% \documentclass[12pt]{article}
% \pagestyle{headings}
% %\RequirePackage[top=2.54cm,right=2.5cm,bottom=2.54cm,left=2.5cm]{geometry}
% 
% 
% \usepackage{graphicx}
% \usepackage{fixltx2e}
% \usepackage{caption}
% \usepackage{subcaption}
% \usepackage{float}
% 
% \usepackage[numbers,round,colon,authoryear]{natbib}
% 
% 
% \usepackage{url}
% \usepackage[english]{babel}
% \usepackage{amsmath}
% \usepackage{amssymb}
% \usepackage{url}
% \usepackage{multirow}
% %\usepackage{fullpage}
% \usepackage{lscape}
% \usepackage{color}
% %\usepackage[font=small,hang]{caption}
% \graphicspath{{./Figures/}}
% 
% % %  
% % % Tweak the graphics placement.
% % \renewcommand{\textfraction}{.1}
% % \renewcommand{\bottomfraction}{.6}
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %12pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %12pt% this is for 12pt
%% \setlength{\textwidth}{16.4cm}
% \setlengt0.3\textwidth}{16.4cm}
% \setlength{\textheight}{22.8cm}
% \setlength{\hoffset}{-1.3cm}
% \setlength{\voffset}{-1.6cm}
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %11pt%11pt
% %11pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %11pt % this is for 11pt
% %11pt\setlength{\textwidth}{16.5cm}
% %11pt\setlength{\textheight}{22.7cm}
% %11pt\setlength{\hoffset}{-2.0cm}
% %11pt\setlength{\voffset}{-2.3cm}
% %11pt%11pt%11pt
% 
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %  razmak izmedju redova 
% %
% %baselinestretch
% %
% \renewcommand{\baselinestretch}{1.4}
% \small\normalsize  % trick from Kopka and Daly p47.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 


\documentclass[12pt]{article}
%\usepackage{subfigure} % was: subfignarms.sty
%\usepackage{psfig}
\usepackage{graphicx}
%\usepackage{timesmt}   % use times typeface
%\usepackage{chicnarm} % <---------------------------------- MOD
%\usepackage{chicago}   % <---------------------------------- MOD

\usepackage[numbers,round,colon,authoryear]{natbib}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{fancyheadings}


\usepackage{url}
\usepackage[english]{babel}
% \usepackage{amsmath}
% \usepackage{amssymb}
\usepackage{url}
% \usepackage{multirow}
%\usepackage{fullpage}
% \usepackage{lscape}




\graphicspath{{./Figures/}}
%\graphicspath{Figures}
%\graphicspath{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/}



\usepackage{fixltx2e}
%\usepackage{caption}
% \usepackage{subcaption}
% \usepackage{float}
% 
% \usepackage[numbers,round,colon,authoryear]{natbib}
% 
% 
% \usepackage{url}
% \usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
% \usepackage{url}
% \usepackage{multirow}
% %\usepackage{fullpage}
% \usepackage{lscape}
% \usepackage{color}
% %\usepackage[font=small,hang]{caption}

%  
% Tweak the graphics placement.
\renewcommand{\textfraction}{.1}
\renewcommand{\bottomfraction}{.85}
\renewcommand{\floatpagefraction}{0.75}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%12pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%12pt% this is for 12pt
\setlength{\textwidth}{16.4cm}
\setlength{\textheight}{22.8cm}
\setlength{\hoffset}{-1.3cm}
\setlength{\voffset}{-1.6cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%11pt%11pt
%11pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%11pt % this is for 11pt
%11pt\setlength{\textwidth}{16.5cm}
%11pt\setlength{\textheight}{22.7cm}
%11pt\setlength{\hoffset}{-2.0cm}
%11pt\setlength{\voffset}{-2.3cm}
%11pt%11pt%11pt


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \lhead{}
% \rhead{}
% \chead{}
% \lfoot{}
% \rfoot{}
% %\cfoot{\pagenumber}
% \cfoot{}
% \addtolength{\headheight}{14pt}
% 
% \setlength{\headrulewidth}{0mm}
% \setlength{\headrulewidth}{0mm}
% %\setlength{\footrulewidth}{0.1mm}
% %\setlength{\footrulewidth}{0.1mm}
% 
% \pagestyle{fancy}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FEI  LOGO definition
\newcommand{\FEI}
{{\sf \uppercase{F}\kern-0.20em\uppercase{E}\kern-0.20em\uppercase{I}}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\today
{\number\day.\space \ifcase\month\or
January\or
February\or
March\or
April\or
May\or
June\or
July\or
August\or
September\or
October\or
November\or
December\fi,\space
\number\year}





% ----------------------------------------------------------------------
%
% TIME OF DAY
%
%
\newcount\hh
\newcount\mm
\mm=\time
\hh=\time
\divide\hh by 60
\divide\mm by 60
\multiply\mm by 60
\mm=-\mm
\advance\mm by \time
\def\hhmm{\number\hh:\ifnum\mm<10{}0\fi\number\mm}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  razmak izmedju redova 
%
%baselinestretch
%
\renewcommand{\baselinestretch}{1.2}
\small\normalsize  % trick from Kopka and Daly p47.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\title{Simulating stiffness degradation and damping in soils via
a simple visco-elastic-plastic model}

%\title{Coupling Multiaxial Elasto-Plasticity  and Viscous Dissipation 
%for Cyclically Loaded Granular Materials}

\author{Federico Pisan{\`o}\textsuperscript{1} and 
        Boris Jeremi{\'c}\textsuperscript{1,2}\\ 
         \\ 
        \textsuperscript{1} {\small University of California, Davis, CA} \\
        \textsuperscript{2}   {\small  Lawrence  Berkeley  National  Laboratory,
        Berkeley, CA}}

\date{}

\maketitle



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
Stiffness degradation and damping represent some of the most well-known aspects of cyclic soil behavior. While standard equivalent linear approaches reproduce these features by (separately) prescribing stiffness reduction  and damping curves, in this paper a multiaxial visco-elastic-plastic model is developed for the simultaneous simulation of both cyclic curves over a wide cyclic shear strain range.

The proposed  constitutive relationship is based on two parallel resisting/dissipative mechanisms,
purely  frictional  (elastic-plastic)  and viscous. The frictional mechanism is formulated as a bounding surface plasticity model with vanishing elastic domain, including  pressure-sensitive failure locus  and  non-associative plastic flow -- which are essential for effective stress analyses. At the same time, the use of the parallel viscous  mechanism  is shown to be  especially beneficial to improve the simulation of the overall dissipative performance.

In order to enable model calibration  on  few standard experimental data, the constitutive equations are purposely kept as simple as possible with a low number of material parameters. Although the model performance is here explored with reference to pure shear cyclic tests, the multiaxial formulation is appropriate for general loading conditions. 
\end{abstract}


\paragraph{Keywords}  stiffness degradation, damping, plasticity,  bounding  surface, viscosity, cyclic
loading




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

Modeling  soil  behavior  under  cyclic/dynamic loading is
crucial   in   most  Geotechnical  Earthquake  Engineering  (GEE)  applications,
including e.g. site response analyses and soil structure interaction (SSI) problems.
In    the    last    decades,   a     number   of   experimental   studies
\citep{Ishihara_1996,  diPrisco_Wood_2012}  pointed  out  the complexity of such
behavior  -- especially in the presence of pore fluid(s) -- characterized
by   non-linearity,   irreversibility,   anisotropy,   barotropy,   picnotropy,
rate-sensitivity,  etc.  In  principle,  a comprehensive soil model should be
capable  of  reproducing all the aspects of the mechanical response for any loading condition,
as well as predicting the occurrence of liquefaction and cyclic mobility, 
distinguishing  the  conditions for shakedown or ratcheting under repeated loads
and  so forth. However, such a perfect model is expected to require too many data
for calibration, which are hardly available in most practical situations.

Conversely,  many  GEE  problems  are  traditionally  solved in  the frequency domain by using
 1D  (equivalent)  linear visco-elastic models, to be calibrated on standard  stiffness degradation ($G/G_{max}$) and
damping ($\zeta$) curves.  Owing to the availability of computer programs for 1D site response analysis (SHAKE
\citep{Schnabel_1972},  EERA  \citep{Bardet_2000}, DEEPSOIL \citep{Hashash2001})
and  SSI  problems  (SASSI  \citep{Lysmer1988}), the visco-elastic approach has
become  more  and  more popular among practitioners, regardless of the following drawbacks:  
\begin{itemize} 

\item [--] most  energy  dissipation 
in  soils comes from frictional inter-granular mechanisms, rather than viscous flow;

\item [--] $G/G_{max}$  and    $\zeta$    curves do not allow to evaluate
irreversible deformations, nor the influence of pore fluid(s);

\item   [--] adopting  1D  shear  constitutive relationships  has  poor  mechanical 
soundness,  since  soil  behavior  exhibits  a pronounced deviatoric-volumetric
 coupling under general multiaxial loading conditions;

\item  [--]  the  meaning  of  cyclic  shear  strain  amplitude  for  the  choice of
$G/G_{max}$  and  $\zeta$  values  is  not evident in the presence of
irregular seismic loads.

\end{itemize}

From the  above  observations    the  need stems for more
physically      consistent    soil    models.    In the last decades,  several  approaches to cyclic modeling have been explored and gradually
refined in the framework of elasto-plasticity,   including  e.g.  ``multi-surface  plasticity'',  ``bounding  surface
plasticity'',  ``generalized  plasticity'' and ``hypoplasticity''.   A number of valuable contributions are worth citing, such as -- to mention
only  a few -- \cite{Mroz78, Prevost1985, wang1990, Borja_1994, Manzari97, Gajo_Wood_1999,
Papadimitriou2002,  Elgamal_etal_2003, dafalias2004, Taiebat2007a, andrianopoulos2010}; recently, it has been also
shown  how  a  good simulation of dynamic properties can be achieved by means of
even  elastic-perfectly  plastic  models,  as  long  as  formulated  in  a probabilistic
elastic-plastic   framework   \citep{Sett2009b}.  Comprehensive
overviews  on  cyclic soil modeling are given by  \cite{Prevost96}, \cite{Zienkiewicz99} and
\cite{diPrisco_Wood_2012}. 

As is well-known, the  major  issues  about  the
practical  use  of elastic-plastic models concern the complexity of the mathematical
formulations and the possible high number of material parameters. For a model to
appeal to practicing engineers,  a  tradeoff is needed between the overall accuracy and the
number  of  parameters  to be calibrated, particularly provided the frequent lack of detailed
\textit{in situ} or laboratory data. This observation led the authors to set up a constitutive model with main following characteristics:
\begin{itemize}
\item [--] multiaxial formulation to cope with any general stress loading condition; 
\item [--] capability of reproducing the main features of soil behavior, such as stiffness reduction and damping over a wide strain range, as well as frictional failure, irreversible deformation and dilatancy;
\item [--] low number of material parameters, to be calibrated on standard experimental data and especially $G/G_{max}$  and  $\zeta$ curves.
\end{itemize}
A constitutive relationship fulfilling the above requirements is hereafter presented as the combination of two resisting/dissipative mechanism, purely frictional (elastic-plastic) and viscous. The  former mechanism has been formulated starting from the work by \cite{Borja_1994}, who proposed a  kinematic-hardening bounding surface von Mises
model for the seismic total-stress analysis of the clayey deposits  at
Lotung  site  in  Taiwan  \citep{Borja99,  Borja_etal_2000}.  In particular, based  on  the former idea by
\citet{Dafalias_Popov_1977} and \citet{Dafalias_1979}, the \citeauthor{Borja_1994}'s  multiaxial
model  is  characterized  by  the  assumption  of vanishing elastic domain, implying  soil  plastification  at  any  load level and a redefinition of the
standard  loading/unloading  criterion. Recently, the vanishing elastic region approach has been also employed by \cite{andrianopoulos2010}, in order to improve the previous model by \cite{Papadimitriou2002} with respect to numerical implementation and integration.

In  this  paper, a similar  bounding surface approach with vanishing elastic region 
is  exploited  to  derive the frictional component of the overall model in the form of a Drucker-Prager effective-stress relationship, incorporating
pressure-sensitive failure  and    non-associative plastic flow.
In
addition, the  model is endowed with 
a  further  viscous  mechanism,  which  can  be  wisely exploited to improve the
simulation  of  the  experimental  damping. Although numerical convenience often
motivates   the   use   of   viscous   dissipation  into  elastic-plastic
computations,  it has a real physical origin, coming from
the time-dependent   processes  taking place at  both  inter-granular  contacts  and grain/pore fluid interfaces.


While most model features are in fact inherited from other previous works, the proposed formulation should be considered as an attempt 
at reconciling traditional (linear equivalent) and advanced (elastic-plastic) cyclic modeling within an effective-stress plasticity framework. In particular, the model is shown to possess reasonable accuracy in reproducing standard modulus reduction and damping curves over a wide strain range, and it is particularly user-friendly because of the low  number of material parameters. The preliminary literature survey put evidence that these advantages are not easily found in most previous effective-stress models. 

In order to clearly illustrate the main modeling ingredients, the basic version of the model is hereafter presented and tested under symmetric cyclic shear loading.  Its convenient mathematical structure will enable in the near future to easily remove the simplifying assumptions that, depending on the specific application, may lead to excessive inaccuracies.

 






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Frictional and viscous dissipative mechanisms}

The  time-domain  finite  element  (FE) solution of dynamic problems is usually
carried  out  by  solving  an  incremental discrete system of the following form
\citep{local-80, FEM-IV-1}:
%
\begin{equation}
\label{eq:FEsystem}
\mathbf{M}\Delta\ddot{\mathbf{U}}+\mathbf{C}\Delta\dot{\mathbf{U}}+\mathbf{K}^{t}\Delta \mathbf{U}=\Delta \mathbf{F}^{ext}
\end{equation}
%
where  $\Delta$  and  dots  stand  respectively  for  step  increment  and  time
derivative,  $\mathbf{U}$  is the generalized DOF vector (nodal displacement for
example),  $\mathbf{F}^{ext}$  the nodal external force vector and $\mathbf{M}$,
$\mathbf{C}$,  $\mathbf{K}^t$  are  the  mass,  damping  and (tangent) stiffness
matrices, respectively.

In  system \eqref{eq:FEsystem} two dissipative sources are readily recognizable,
namely    the    viscous    (velocity-proportional)    and    the    frictional
(displacement-proportional)  terms  \citep{local-87}. While the latter is given
by  the  elastic-plastic tangent stiffness $\mathbf{K}^t$, the
viscous    term    related  to  the  damping  matrix  $\mathbf{C}$ can represent
 interaction of solid skeleton and pore
fluid, and constitutive   rate-sensitiveness    of    the   soil   skeleton.
The  above  combination  of
frictional  and  viscous dissipation can be interpreted in terms of two distinct
effective stress components acting on the soil skeleton:
%
\begin{equation}
\label{eq:stress_split}
\sigma_{ij}=\sigma_{ij}^{f} + \sigma_{ij}^{v} 
\end{equation}
%
where  the  effective stress tensor $\sigma_{ij}$ has been split into frictional
(elastic-plastic)  and viscous stresses\footnote{Henceforth, effective stresses
are  exclusively accounted for}. From a rheological point of view, the resulting
scheme can be defined as visco-elastic-plastic. In what follows, the frictional component
is  first  specified  via the formulation of the bounding surface
model  with  vanishing elastic domain; then, the role played by the linear viscous term is discussed.

Index  tensor  notation  is  used,  along with the standard Einstein
convention  for  repeated indices; the norm of any second-order tensor $x_{ij}$
is  defined  as  $\Vert x_{ij}\Vert=\sqrt{x_{ij}x_{ij}}$, whereas the deviatoric
component              can              be              extracted             as
$x_{ij}^{dev}=x_{ij}-x_{kk}\delta_{ij}/3   $  ($\delta_{ij}$  is  the
Kronecker delta). In accordance with usual Solid Mechanics conventions, positive
tensile  stresses/strains  are  considered,  whereas  --  as  is  done  in Fluid
Mechanics -- only the isotropic mean pressure is positive if compressive.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Bounding surface frictional model with vanishing elastic domain}

The frictional component of the  model being proposed  is formulated by generalizing the previous constitutive relationship  by
\cite{Borja_1994}. 
%As was expected, the introduction of pressure-dependence into the
%constitutive  equations  implies somewhat more involved %derivations, so that the
%analytical details skipped in this section are reported in Appendix \ref{app:A};
In what follows, the  superscript  $f$  referring  to  the  frictional  component  of  the global
effective  stress  (Equation \eqref{eq:stress_split}) is avoided for the sake of
brevity.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Elastic relationship}

Provided  the  usual additive combination of (incremental)  elastic and plastic strain components
$d \epsilon_{ij}= d \epsilon_{ij}^e   +   d \epsilon_{ij}^p$,  the  incremental  linear
elastic Hooke's law is expressed as follows:
%
\begin{equation}
\label{eq:elastic}
d\sigma_{ij}=D_{ijhk}^e\left ( d\epsilon_{hk}-d\epsilon_{hk}^p \right )\Longrightarrow
\left\{\begin{matrix}
ds_{ij}=2G_{max}\left ( de_{hk}-de_{hk}^p \right )\\
\\ 
dp=-K\left ( d\epsilon_{vol}-d\epsilon_{vol}^p \right )\\
\end{matrix}\right.
\end{equation}
%
where  $d$  stands  for a differentially small increments and  $D_{ijhk}^e$ is the fourth-order
elastic  stiffness  tensor.  Equation \eqref{eq:elastic} also points out the elastic deviatoric/volumetric decoupling, in which $p=-\sigma_{kk}/3$, $\epsilon_{vol}=\epsilon_{kk}$, $s_{ij}=\sigma_{ij}^{dev}$ and $e_{ij}=\epsilon_{ij}^{dev}$ stand for mean
stress,    volumetric  strain, stress      deviator and
 strain deviator, respectively.
The shear modulus $G_{max}=E/2\left(1+\nu\right)$
and  the  bulk  modulus  $K=E/3\left(1-2\nu\right)$  are  derived from the Young
modulus  $E$ and the Poisson's ratio $\nu$. Henceforth, $G_{max}$ will be always
used  for  the  elastic  small-strain  shear modulus, whereas the secant cyclic
shear stiffness will be referred to as $G$.

Although soil elasticity is known to be non-linearly pressure-dependent, a classical linear formulation has been here maintained to simulate constant-pressure cyclic shear tests, exclusively.


%
%\begin{equation}
%\label{eq:elastic_dev}
%ds_{ij}=2G_{max}\left ( de_{hk}-de_{hk}^p \right )\\
%\end{equation}
%

%
%\begin{equation}
%\label{eq:elastic_vol}
%dp=-K\left ( d\epsilon_{vol}-d\epsilon_{vol}^p \right )\\
%\end{equation}
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Drucker-Prager yield and bounding loci}

A    conical   Drucker-Prager  yield  locus is first introduced, similar to
what used by  \citet{Prevost1985} and \citet{Manzari97}:

%
\begin{equation}
\label{eq:yieldlocus}
f_y=\dfrac{3}{2} \left(s_{ij}-p\alpha _{ij}\right)\left(s_{ij}-p\alpha _{ij}\right)-k^2 p^2=0
\end{equation}
%
where   $\alpha_{ij}$   is   the   so called   deviatoric   back-stress  ratio
($\alpha_{kk}=0$)  governing  the  kinematic  hardening  of the yield
surface;  $k$  is  a  parameter determining the opening angle of the cone. The  variation  of  the  back-stress  ratio
$\alpha_{ij}$  in \eqref{eq:yieldlocus} leads to a rigid rotation of the yield
locus and, therefore, a rotational kinematic hardening.

In the spirit of standard bounding surface plasticity, the  yield locus must always reside within an outer surface (the so-called bounding surface), here assumed in the form of a further Drucker-Prager cone (non kinematically hardening, fixed
in size):
%
\begin{equation}
\label{eq:boundingsurface}
f_B=\dfrac{3}{2} s_{ij}s_{ij}-M^2p^2=0
\end{equation}\\
%
where $M$ provides the bounding cone opening and, as a consequence, the material
shear strength.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Plastic flow and translation rules}

The plastic flow of soils is in general  non-associative  \citep{local-38} and   gives rise to  plastic  contractancy  or  dilatancy
depending  on  whether  loose or dense materials are considered. Here, the plastic
flow rule is borrowed from \cite{Manzari97}:

%
\begin{equation}
\label{eq:flow_rule}
d\epsilon_{ij}^p=d\lambda\left(n^{dev}_{ij} - \dfrac{1}{3} D\delta_{ij} \right)
\end{equation}
%
where  $d\lambda$ is the plastic multiplier, $n_{ij}^{dev}$ is a deviatoric unit
tensor  ($\Vert n_{ij}^{dev}\Vert=1$) and $D$ is a dilatancy coefficient defined
as \citep{Manzari97}:

%
\begin{equation}
\label{eq:dil_coeff}
D = \xi \left(\alpha^d_{ij}-\alpha_{ij}\right)n^{dev}_{ij}=\xi \left(\sqrt{\dfrac{2}{3}}k_dn^{dev}_{ij}-\alpha_{ij}\right)n^{dev}_{ij}
\end{equation}
%
in  which  $\xi$  and  $k_d$ are two positive constitutive parameters. While the
former  controls  the amount of volumetric plastic strain, the latter determines
the  position  of  the so called ``dilatancy surface'' and rules the transition
from contractive ($D>0$) to dilative ($D<0$) behavior under undrained triaxial conditions\footnote{Under different loading conditions this transition is not ``exactly'' governed only by the location of the current stress state with respect to the dilatancy surface}.

For the sake of simplicity (and regardless of some experimental evidences), the kinematic hardening evolution of the yield locus and the direction of the deviatoric plastic strain increment are related through the standard
Prager's    rule \citep{Borja_1994}:

%
\begin{equation}
\label{eq:prager_rule}
d\alpha_{ij}=\Vert d\alpha_{ij}\Vert n_{ij}^{dev}
\end{equation}
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Vanishing elastic region}

As  previously mentioned, the main feature of the frictional model concerns
the    vanishing  elastic  domain, corresponding with the limit $k\rightarrow 0$
in  Equation~\eqref{eq:yieldlocus}.    Accordingly,   the  Drucker-Prager  cone 
reduces to its axis, so that:
%
\begin{equation}
\label{eq:vanish_locus}
\lim_{k\rightarrow 0} f_y= 0        
\Rightarrow
\dfrac{s_{ij}}{p} =r_{ij}=\alpha_{ij}
\Rightarrow 
d\alpha_{ij}=dr_{ij}=\dfrac{ds_{ij}}{p}-\dfrac{s_{ij}}{p^2}dp
\end{equation}
%
where $r_{ij}$ is the deviatoric stress ratio tensor \citep{Manzari97}. After  substituting  the  Prager's  rule  \eqref{eq:prager_rule}  into \eqref{eq:vanish_locus} it results:
%
\begin{equation}
\label{eq:from_vanishlocus}
n^{dev}_{ij}=\dfrac{dr_{ij}}{\Vert dr_{ij}\Vert}=\dfrac{ds_{ij}-\alpha_{ij}dp}{\Vert ds_{ij}-\alpha_{ij}dp\Vert}
\end{equation}
%
From Equation \eqref{eq:from_vanishlocus} it can be inferred
that:
\begin{enumerate}
\item  since the direction of the  plastic
strain  increment    overall depends  on the stress increment $d{\sigma}_{ij}$,  the resulting
constitutive  formulation is by definition ``hypoplastic'' \citep{Dafalias86};
\item   hydrostatic   stress  increments  ($d{s}_{ij}=0$)  from  initial
hydrostatic      states   ($\alpha_{ij}=0$)  produce $n_{ij}^{dev}=0\rightarrow \epsilon_{ij}^{dev}=0$;
\item the deviatoric plastic strain increment is along the direction of the deviatoric stress ratio increment. This directly comes from the use of the above Prager's rule in combination with the vanishing elastic region. Although this finding is not in full agreement with general experimental evidences, it will not prevent satisfactory cyclic $G/G_{max}$ and $\zeta$ curves to be obtained.
\end{enumerate}

 Starting from Equation \eqref{eq:vanish_locus},
the  norm  $\Vert d{\alpha}_{ij}\Vert=\Vert dr_{ij}\Vert$  can be  specified for the case of radial
loading  paths  in  the deviatoric plane, which are characterized by $d{p}=0$
and the coaxiality of $s_{ij}$ and $ds_{ij}$. Under these loading conditions, simple
manipulations  lead to find:
%
\begin{equation}
\label{eq:alpha_incr_rad}
\Vert d\alpha_{ij}\Vert=\sqrt{\dfrac{2}{3}}\dfrac{dq}{p}
\end{equation}\\
%
where  $q=\sqrt{3/2}\Vert  s_{ij}  \Vert$ stands for the usual deviatoric stress invariant.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Hardening relationship and plastic multiplier}

An incremental hardening relationship is directly established \citep{Borja_1994}:
%
\begin{equation}
\label{eq:hard_relation}
d{q}= \sqrt{\dfrac{2}{3}}H\Vert d{e}_{ij}^p \Vert
\end{equation}
%
where $H$ is the hardening modulus. Then, the substitution of
both   the   flow  rule  \eqref{eq:flow_rule}  and  the  hardening  relationship
\eqref{eq:hard_relation} into \eqref{eq:alpha_incr_rad} leads to:

%
\begin{equation}
\label{eq:norm_dalpha}
\Vert d{\alpha}_{ij} \Vert=\dfrac{2}{3}\frac{Hd{\lambda}}{p}
\end{equation}
%
%
By  equaling  the two definitions of $ds_{ij}$ arising from  Equations \eqref{eq:elastic}-deviatoric  and
\eqref{eq:vanish_locus}, and then using Equations \eqref{eq:elastic}-volumetric  and \eqref{eq:norm_dalpha} the following relationship is obtained:

%
\begin{equation}
2G_{max}\left(d{e}_{ij}-d{\lambda}n^{dev}_{ij}\right)=
\Vert                      d{\alpha}_{ij}\Vert                      n^{dev}_{ij}
p+\alpha_{ij}d{p}=\dfrac{2}{3}\dfrac{Hd{\lambda}}{p}n^{dev}_{ij}p-\alpha_{ij}K 
\left(d{\epsilon}_{vol}+d{\lambda}D\right)
\end{equation}
%
whence:
%
\begin{equation}
\label{eq:plast_mult}
d{\lambda}=\dfrac{ 2G_{max} d{e}_{ij}n_{ij}^{dev}
+Kd{\epsilon}_{vol}\alpha_{ij}n_{ij}^{dev}  }{2G + \dfrac{2}{3}H -KD\alpha_{ij}n_{ij}^{dev} }
\end{equation}
%
Equation   \eqref{eq:plast_mult}   represents   the   consistent  frictional
generalization  of  Equation  (18) in \cite{Borja_1994}, as well as the limit of
Equation (12) in \cite{Manzari97}\footnote{Different signs result because of the
opposite  sign conventions adopted by these authors} for a vanishing size of the yield locus.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Projection rule, hardening modulus and unloading criterion}

The  theory of bounding  surface  plasticity    relies on the basic concept that the
plastic  modulus  explicitly  depends on the distance between the current stress
state  and an \textit{ad hoc} stress projection onto the bounding surface. On this issue,  benefits and pitfalls of different projection rules have been thoroughly discussed by \cite{andrianopoulos2005}. 

Here, the stress projection in the $\pi$-plane (deviatoric stress ratio plane) is assumed to be along the direction of $\alpha_{ij}-\alpha_{ij}^0 $:
%
\begin{equation}
\label{eq:projection}
{\alpha}_{ij}^b=\alpha_{ij}+\beta \left ( \alpha_{ij}-\alpha_{ij}^0 \right )
\end{equation}\\
%
where $\beta$ is a scalar distance coefficient, while $\alpha_{ij}^0 $ is the back-stress ratio at the last loading reversal (Figure \ref{fig:projection}). As observed by \cite{andrianopoulos2010}, projection rules of the type \eqref{eq:projection} are ``stable'', in the sense that small
perturbations of the loading direction do not  affect severely   the location of the stress projection on the bounding surface.


The coefficient  $\beta$  is obtained by enforcing  the projected stress ${\sigma}_{ij}^b$
to lie on the bounding surface (Equation~\eqref{eq:boundingsurface}):
%
\begin{figure}  [!htb]
 \centering
\includegraphics [width=.6\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/projection.pdf}
{projection2.pdf}
 \caption{Representation of the projection rule \eqref{eq:projection}}
 \label{fig:projection}
\end{figure}
%
\begin{equation}
\label{eq:proj_condition}
\dfrac{3}{2}{s}_{ij}^b{s}_{ij}^b=M^2\left({p}^b\right)^2\Longleftrightarrow
\dfrac{3}{2}{\alpha}_{ij}^b{\alpha}_{ij}^b=M^2
\end{equation}
%
%
that is by substituting \eqref{eq:projection} into \eqref{eq:proj_condition} and then
 deriving  the positive root of the following algebraic equation\footnote{The adoption of a simple Drucker-Prager-type bounding surface allows $\beta$ to be analytically obtained. This would not be the case in the presence of more complex $\pi$-sections}:
%
\begin{equation}
\label{eq:beta}
\left \| \alpha_{ij}-\alpha_{ij}^0 \right \|^2\beta^2+2\alpha_{ij}\left ( \alpha_{ij}-\alpha_{ij}^0 \right )\beta+\left \| \alpha_{ij} \right \|^2-\dfrac{2}{3}M^2=0
\end{equation}
%
The  analytical  relationship  between  $H$ and
$\beta$ is chosen to fulfill the following requirements:
\begin{enumerate}
\item $H\left(\beta=0\right)=0$, that is full mobilization of the material strength  when the stress image point lies on the bounding surface;
\item  $H\left(\beta\rightarrow\infty\right)\rightarrow\infty$, that is (instantaneous) recover of the elastic stiffness at the onset of load reversal (i.e. $\alpha_{ij}=\alpha_{ij}^0$).
\end{enumerate} 
In  this  case,  the  following $H-\beta$ relationship has been introduced because of its simplicity:
%
\begin{equation}
\label{eq:hard_modulus}
H=ph\beta^m
\end{equation}
%
in which $h$ and $m$ are two additional constitutive parameters.

In the absence of a finite-sized elastic region, the occurrence of unloading cannot be checked by comparing the loading direction and the local normal to the yield locus. This requires the difinition of an alternative unloading criterion, coinciding here with providing a rule for updating $\alpha_{ij}^0$.  \cite{Borja_1994} proposed a criterion
based  on  the observation that the hardening modulus $H$ increases at the onset of  unloading, so that -- as long as $H\left(\beta\right)$ is a monotonically
increasing  function --  instantaneous  unloading  is  assumed   whenever
$d{H}>0$, i.e. $d{\beta}>0$. However, the authors experienced that such an unloading criterion lacks robustness in numerical computations under complex/irregular loading paths, because the small values assumed by $\beta$ (especially close to the bounding surface) can be easily corrupted even by numerical inaccuracies; as a result,  unrealistic unloading is often likely to arise. While different proposals are available in literature to overcome this problem \citep{andrianopoulos2005, andrianopoulos2010}, the following unloading criterion has been here preferred:
%
\begin{equation}
\label{eq:unloading}
\left(\alpha_{ij}-\alpha_{ij}^0\right)dr_{ij}<0\Longleftrightarrow
\left(\alpha_{ij}-\alpha_{ij}^0\right)n_{ij}^{dev}<0
\end{equation}
%
in which the coaxiality of $d\alpha_{ij}=dr_{ij}$ and $n_{ij}^{dev}$ (Prager's rule  \eqref{eq:prager_rule}) has been exploited. 


It should be also borne in mind that this kind of updating criterion is likely to give rise to ``overshooting'' phenomena under general loading paths \citep{Dafalias86}. Overshooting takes place when, after loading along a given direction, a very small unloading implies an updated $\alpha_{ij}^0$ before reloading along the original direction. The updated small distance between the current and the reversal stresses determines an unrealistically high reloading stiffness, so that ``the corresponding stress-strain curve will overshoot the continuation of the previous curve which would have occurred if no unloading/reverse loading/reloading had taken place'' \citep{Dafalias86}. While this shortcoming can be usually observed under irregular seismic loading, it will not be detected in the results being presented, exclusively concerning symmetric sinusoidal cyclic loading. A recent discussion on the remediation of overshooting is provided by \cite{Kan2014}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Possible refinements}

The  frictional  model  has been developed trying to keep the number of material
parameters  as  low as possible, even with a non-associated flow rule. However,
it  is  worth  mentioning  which  kind  of  improvements  might be introduced if
required by the problem under examination.

It  should  be first noted that, as a Drucker-Prager type bounding surface has
been  adopted,  the  material shear strength is unaffected by the Lode angle, so
that  for  instance  the  same  failure  obliquity  is  predicted  for  triaxial
compression  and  extension. This drawback could be easily remedied by modifying
the  deviatoric  cross-section of the bounding surface itself, e.g. by adopting
the      well  known     Mohr-Coulomb  deviatoric  locus  or  other  smooth  loci
 \citep{Matsuoka_Nakai_1974, local-103, Lade77}. A change
in the
deviatoric cross-section would negligibly influence the overall formulation, as
just  the  evaluation  of  the  projection distance $\beta$ (Equation  \eqref{eq:proj_condition}).

Secondly,  the  present  version  of the model cannot predict a possible softening
behavior  of  the  soil,  usually  taking  place in the case of dense materials.
Softening  could  be  accounted  for  by incorporating a further
isotropic  hardening  mechanism  at  the  bounding surface level, allowing for a
gradual shrinkage of the outer surface during plastifications.

Another  relevant  point  concerns  the  fact that different parameters must be calibrated for different void ratios of the same granular material, as if
distinct  materials  were  indeed  considered.  As  a matter of fact, continuous transitions  from loose to dense states (and vice versa) spontaneously take place
during straining: this aspect has been successfully addressed and reproduced via
the   concept   of   ``state   parameter''  \citep{Been_Jefferies_1985,  Wood94,
Manzari97}, which could be also introduced into a critical-state version of the
proposed model. Also, the state parameter concept represents a natural way to reproduce softening as an effect of the material evolution from peak strength to critical state conditions.

While the above aspects concern both monotonic and cyclic loading conditions, further issues could be  addressed  with more specific reference to  cyclic/dynamic loading, as for instance fabric effects \citep{Papadimitriou2002, dafalias2004} and anisotropy, the evolution toward shakedown or ratcheting under a large number of loading cycles \citep{diPrisco_Wood_2012}, the occurrence of cyclic mobility \citep{Elgamal_etal_2003},  rate-sensitiveness and related frequency effects. 

Apparently, refining the model formulation in the light of the above observations would result in more  accurate predictions/simulations, implying though higher difficulties in
terms   of   calibration,   implementation  and,  as  a  consequence,  practical
employment. Conceiving a model reasonably accurate but still ``user-friendly'' is the main goal of the present study, in order  to provide engineers with a tool fitting standard cyclic modeling concepts (modulus reduction and damping curves) in a 3D elastic-plastic fashion. Accordingly, while looking for a compromise between accuracy and simplicity,  some experimental evidences have been purposely disregarded on the modeling side.

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{The role of linear viscous damping}

An  additional  viscous  mechanism  (Equation  \eqref{eq:stress_split})  is
usually  available exploited  in  most finite  element  (FE)  codes,  even  though it is not
directly  included in the constitutive model. Indeed, many numerical programs solve
discrete  systems  with  a  viscous damping term (Equation \eqref{eq:FEsystem}),
usually assembled as a linear combination of the mass and the (elastic) stiffness
matrices (Rayleigh formulation \citep{local-87, Chopra01}):
%
\begin{equation}
\label{eq:dampingmatrix}
\mathbf{C}=a_0\mathbf{M}+a_1\mathbf{K}^e
\end{equation}
%
where $a_0$ and $a_1$ are two constant parameters, related to the $n^{th}$ modal damping ratio $\zeta_n$ of the discrete structural system.

It could be easily shown that a constitutive viscosity in the form:
%
\begin{equation}
\label{eq:viscous_stress}
\sigma_{ij}^v = D_{ijhk}^v\dot{\epsilon}_{hk}
\end{equation}
%
gives   rise   to   a  stiffness-proportional  damping  matrix,  which  can  be
equivalently  reproduced  through  the  following  calibration  of  the Rayleigh
damping parameters \citep{Borja_etal_2000, Hashash2002}:
%
\begin{equation}
\label{eq:Rayleigh_calibration}
a_0=0 \qquad a_1=\dfrac{2\zeta_0}{\omega}
\end{equation} 
%
The  calibration  \eqref{eq:Rayleigh_calibration}  establishes  the  same  ratio
between    tangential/bulk   elastic   and   the   viscous   moduli,   that   is
$G^e_{max}/G^v_{max}=K^e/K^v$.  More  importantly,  a damping ratio $\zeta_0$ is
ensured  for  a  given  circular  frequency  $\omega$,  as  long as the parallel
resisting  mechanism  ($\sigma_{ij}^f$)  is  purely  elastic;  as a consequence,
provided the $a_1$ value at the beginning of the analysis, modal frequencies and
the corresponding damping ratios are linearly related.

It  is also worth remarking some further points about the implications of linear
viscous  damping  in conjunction with non-linear soil models. If a soil element
undergoes  an  imposed  shear  strain  history,  the overall shear stress/strain
cycles    $\tau-\gamma$    differ   from   the   purely   frictional   component
$\tau^f-\gamma$, this difference being due to the viscous shear stress $\tau^v$.
As  will be shown in next section, the viscous component implies smoother cycles
and  avoid  the sharp transitions at stress reversal usually exhibited by purely
elastic-plastic   responses   \citep{Borja_etal_2000}.   However,  the  overall
$G/G_{max}$  ratio  between  the  average  cyclic  stiffness  and the
elastic shear modulus is unaffected by viscosity.

As   far   as   the   damping   ratio  is  concerned,  its  standard  definition
\citep{Kramer96}  can  be  easily  adapted  to  point out the frictional/viscous
splitting of the energy $\Delta W$ dissipated in a loading cycle:

%
\begin{equation}
\label{eq:damp_def}
\zeta=\dfrac{\Delta W}{2\pi G\gamma^2_{max}}=\dfrac{\Delta W^{f}+\Delta W^{v}}{2\pi G\gamma^2_{max}}=\zeta^f+\zeta^v
\end{equation}
%
where  $\gamma_{max}$  is  the imposed cyclic shear strain amplitude and $G$ the
corresponding  (secant)  cyclic  shear  stiffness.  As $\gamma_{max}$ approaches
zero,  the  plastic  dissipation tends to zero as well, so that $\zeta=\zeta^v$;
therefore,   the   Rayleigh   parameter   $a_1$  can  be  calibrated  to  obtain
$\zeta\left(\gamma_{max}\rightarrow0   \right)=\zeta_0$  for  a  given  circular
frequency  $\omega$  (see  Equation  \eqref{eq:Rayleigh_calibration}). This is a
desirable  feature  of  the model, as natural soils are well known to dissipate
energy at even very small strain amplitudes.

At  progressively  larger  strains,  both  the frictional and viscous components
contribute   to   the   global   damping,  although  the  relative  quantitative
significance  is  hard to assess \textit{a priori}. In addition, 
the   viscous   component   of   the
$\zeta-\gamma_{max}$  curve  is  not  constant,  since  $\zeta^v$ depends on the
strain-dependent  secant  modulus $G\left(\gamma_{max}\right)$ and, implicitly,
on   the   strain  rate.  This  is  different    to  what  has  been  argued  by
\citet{Borja_etal_2000}.

As  an  example,  consider  the  response  of  an elastic-perfectly  plastic  model with
additional   viscosity   under   a  sinusoidal  shear  excitation  $\gamma\left(
t\right)=\gamma_{max}\sin\left(\omega     t\right)$.    The  simplicity  of  the
elastic-perfectly  plastic  response allows derivation of instructive analytical
formulas for the
$G/G_{max}$  and  the damping ratios, even in the presence of viscous
dissipation. While $\gamma_{max}<\gamma_y$ (yielding shear strain), the material
behavior is linear elastic, so that $G/G_{max}=1$
and  $\zeta$  equals  the  purely viscous contribution at $\gamma_{max}=0$, i.e.
$\zeta=\zeta_0$;  $\gamma_y$  depends  on  the  elastic  stiffness and the shear
strength  of  the material, $\gamma_y=\tau_{lim}/G_{max}$, where $\tau_{lim}$ is
the  limit  (frictional)  shear  stress  for  a  given  confining  pressure. For
$\gamma_{max}>\gamma_y$    plastifications    take    place    with    a    flat
elastic-perfectly plastic $\tau^f-\gamma$ branch, and the following expressions
can be easily derived:

%
\begin{equation}
\dfrac{G}{G_{max}}=\frac{\tau_{lim}}{G_{max}\gamma_{max}}
\end{equation}
%

%
\begin{equation}
\label{eq:damp_PP}
\zeta=\dfrac{\Delta   W^{f}+\Delta  W^{v}}{2\pi  G\gamma^2_{max}}  =\underbrace{
\dfrac{2}{\pi}\left   (   1-\dfrac{\tau_{lim}}{G_{max}\gamma_{max}}   \right   )
}_{\zeta^f}
+\underbrace{\zeta_0\dfrac{G_{max}\gamma_{max}}{\tau_{lim}}}_{\zeta^v}
\end{equation}
%
In Figure \ref{fig:damping_PP} the  $G/G_{max}$ and $\zeta$ ratios are
plotted  for  increasing  $\zeta_0$  values.  As  $\gamma_{max}$  increases, the
frictional  damping  tends  to  $2/\pi\approx0.63$,  while the viscous one keeps
increasing  because of the reduction in the secant stiffness and the increase in
the  shear  strain rate (depending on the strain amplitude). Hence, the value of
$\zeta_0$  is  to  be  carefully chosen, in order to avoid excessive dissipation
when medium/large strains are induced by the loading process.

\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/analytical.pdf}
{analytical.pdf}
 \caption{$G/G_{max}$  and  damping  curves  for  a elastic-perfectly plastic
 model  with  linear viscous damping at varying $\zeta_0$ ($\tau_{lim}$=100 kPa,
 $G_{max}$=100 MPa)}
 \label{fig:damping_PP}
\end{figure}

The   fact   that  the  viscous  mechanism  can  modify  the  purely  frictional
$\zeta-\gamma_{max}$ curve without altering the cyclic stiffness degradation can
be   fruitfully   exploited   to  improve 
experimental-numerical  agreement  in  terms  of  energy
dissipation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\section{Model performance and calibration}

The frictional mechanism of the above model is characterized by a very low number of material parameters, namely the following seven:

\begin{itemize}

\item  two  elastic  parameters,  the  Young  modulus  $E$ (or the shear modulus
$G_{max}$) and the Poisson's ratio $\nu$;

\item  the  shear  strength  parameter  $M$  for  the definition of the bounding
surface (Equation \eqref{eq:boundingsurface});

\item  the flow rule parameters, $\xi$ and $k_d$, governing the increment of the
volumetric  plastic strain under shearing and the size of the dilatancy surface,
respectively (Equation \eqref{eq:flow_rule});

\item  the  hardening parameters $h$ and $m$ for the dependence of the hardening
modulus  on the distance coefficient $\beta$ (Equation \eqref{eq:hard_modulus}),
affecting the pre-failure deformational behavior and, in overall, the resulting
dynamic properties ($G/G_{max}$ and damping curves).

\end{itemize}

Provided  a  reasonable  value  for  the  Poisson's  ratio (usually in the range
$0.25 - 0.4$),  the  small-strain  elastic stiffness can be evaluated from dynamic
laboratory  (RC  tests) or \textit{in situ} (seismic geophysical surveys) tests.
As  far  as the shear strength is concerned, the parameter $M$ can be related to
the friction angle $\phi$ as follows:
%
\begin{equation}
\label{eq:M}
M=\dfrac{6\sin\phi}{3\pm\sin\phi}
\end{equation}
%
to reproduce triaxial compression (sign $-$ in \eqref{eq:M}) or extension
(sign  $+$  in  \eqref{eq:M}) failure conditions. While different  bounding deviatoric  sections would easily capture both
compressive and extensive limits \citep{Manzari97}, the calibration of a circular deviatoric locus can be tackled -- as usual -- by setting a trade-off  $M$ in between the bounding values in Equation \eqref{eq:M} (this is appropriate, for instance, for plane strain problems). Also, since no strain-softening is reproduced, the peak strength or the residual $\phi$ is to be considered on the basis of the specific purpose of the analysis, i.e. depending on whether optimistic or safe assessments are needed. In any case, the present version of the model is not well suited for problems where simulating failure is particularly relevant (e.g. slope stability analyses).

The  calibration of the flow rule parameters, $\xi$ and $k_d$, requires at least a  triaxial  test  to be performed, in order to obtain some information about the volumetric behavior. In particular, under undrained triaxial conditions, $k_d$ coincides with the stress ratio ($\eta=q/p'$) characterizing the so-called ``phase transformation line'' \citep{Ishihara_etal_1975} and determining the compactive/dilative transition. In the case of loose (compactive) materials, since fixed bounding and dilatancy surfaces are considered in this version of the model, $k_d=M$ is set  to ensure a compactive behavior vanishing as the limit external locus is achieved.  Figure \ref{fig:res1} shows the predicted triaxial response
for  three  different  values of $k_d$ (and fixed $\xi$), that is by varying the
opening  angle of the dilatancy surface (the employed parameters are reported in
the figure caption, where $p_0$ stands for the initial mean pressure).

\begin{figure}  [!htb]
 \centering
\includegraphics [width=.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res1a.pdf}
{res1a.pdf}
\includegraphics [width=.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res1b.pdf}
{res1b.pdf}
 \hfill
 \caption{Predicted  triaxial  responses  for  different  dilatancy  surfaces  (
 $p_0$=100   kPa,   $G_{max}$   =   4   MPa,   $\nu$=0.25,   $M$=1.2,   $\xi$=1,
 $h$=$G$/(1.5$p_0$), $m$=1)}
 \label{fig:res1}
\end{figure}

While   the  limit  stress  deviator  $q$  is  exclusively  given  by  $M$,  the
pre-failure  behavior  is influenced by the plastic deformability and therefore
by  $k_d$.  The model possesses sufficient flexibility to reproduce contractive,
dilative  or contractive/dilative behavior; also, such a feature is necessary to
reproduce  undrained conditions (liquefying and non-liquefying responses),
this being a further motivation for non-associativeness when dealing with sandy
materials.


Figure  \ref{fig:res2} exemplifies the  response predicted under pure shear (PS)
cyclic    loading,    applied    as    a   sinusoidal   shear   strain   history
($\gamma_{max}=0.2\%,  20 \%$ , period $T$=2$\pi$ s) at constant normal stresses
(and  thus  constant mean pressure $p_0$ as well). This corresponds with a radial
loading  path  on the deviatoric plane); for the sake of clarity, the volumetric
plastic  response  has  been  inhibited  ($\xi$=0),  in  order  to  evaluate the
deviatoric  mechanism  exclusively.  Both  purely  frictional  (solid  line) and
frictional/viscous (dashed line) responses are plotted.

\begin{figure}  [!htb]
 \centering
\includegraphics [width=.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res2a.pdf}
{res2a.pdf}
\includegraphics [width=.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res2b.pdf}
{res2b.pdf}
 \hfill
 \caption{Predicted pure shear response at two different shear strain amplitudes ($p_0$=100 kPa, $T$=2$\pi$ s, $\zeta_0$ = 0.006,
$G_{max}$= 4 MPa, $\nu$=0.25, $M$=1.2, $k_d$=$\xi$=0, $h$=$G$/(1.5$p_0$), $m$=1)}
 \label{fig:res2}
\end{figure}

\begin{figure}  [!htb]
 \centering
\includegraphics [width=.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res2c.pdf}
{res2c.pdf}
 \caption{Detail of stress reversals for the pure shear response in Figure \ref{fig:res2} ($\gamma_{max}=20\%$)}
 \label{fig:res2c}
\end{figure}

Owing  to  the  kinematic  hardening  of  the vanished yield locus, the model can  reproduce both the Bauschinger and the Masing effects, the latter
implying  the stabilization of the cyclic response to take place after more than
one  loading  cycle.  As  expected, the additional viscous damping increases the
area  of  the  cyclic loop and therefore the overall dissipated energy; however,
the  effect  of the viscous dissipation becomes significant only at medium-high
shear  strains,  corresponding  --  for a given loading frequency -- with higher
strain  rates.  Further,  viscosity  causes  the aforementioned ``smoothing'' of
stress  reversals,  as  it  can  be  noticed in Figure \ref{fig:res2c} by
comparing the purely frictional and the frictional/viscous responses.

Given  the  elastic  stiffness  and  the  strength of the soil, the shape of the
resulting    loading  cycles is totally governed by the hardening properties, by
$h$ and $m$ in Equation~\eqref{eq:hard_modulus}: this directly affects the
simulation  of experimental $G/G_{max}$ and damping curves, which can
be therefore exploited for the calibration of both $h$ and $m$. As is proven in  Appendix  \ref{app:A},  the following
equality  holds  under  PS  loading  conditions:
%
\begin{equation}
\label{eq:hm_calibration}
1    =    \dfrac{G}{G_{max}}\left[    1   +   \dfrac{6G_{max}}{hp_0\gamma_{max}}
\int_{0}^{\gamma_{max}}                         \left(\dfrac{\gamma}{\tau_{lim}/
G-2\gamma+\gamma_{max}}\right)^md\gamma\right]
\end{equation}
%
where  $\tau_{lim}=Mp_0/\sqrt{3}$.  Relationship  \eqref{eq:hm_calibration}  has
been  obtained  by  integrating the constitutive equations over the first loading
cycle,   and   represents   the   frictional  counterpart  of  Equation  (6)  in
\cite{Borja_etal_2000}  --  as  is  testified  by  the explicit influence of the
confining  pressure  $p_0$. The proper use of Equation \eqref{eq:hm_calibration}
requires first the choice of two meaningful points on the $G/G_{max}$
experimental  curve,  i.e.  two  $\left(\gamma_{max},G/G_{max}\right)$  couples;
then,  the  unknowns  $h$  and  $m$  are obtained by solving the integral system
arising  from  the  specification of Equation \eqref{eq:hm_calibration} for both
selected $\left(\gamma_{max},G/G_{max}\right)$ couples.

Figure \ref{fig:res3} illustrates the result of the above calibration procedure,
applied  on  the $G/G_{max}$ and $\zeta$ curves for sands implemented
into    the   code   EERA   \citep{Bardet_2000}   and   formerly   obtained   by
\cite{Seed_Idriss_1970}.

\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res3.pdf}
{res3.pdf}
 \hfill
 \caption{Comparison  between  experimental and simulated $G/G_{max}$
 and  damping  curves ($p_0$=100 kPa, T=2$\pi$ s, $\zeta$ = 0.003, $G_{max}$ = 4
 MPa, $\nu$=0.25, $M$=1.2, $k_d$=$\xi$=0, $h$=$G$/(112$p_0$), $m$=1.38)}
 \label{fig:res3}
\end{figure}

Since     Equation  \eqref{eq:hm_calibration}    exclusively  accounts  for  the
$G/G_{max}$  curve,  the  very  satisfactory  agreement  in  terms of
stiffness  degradation  (viscosity  has no effect on it) should not surprise. On
the other hand, once $h$ and $m$ are set, the predicted damping curve may or may
not match the experimental outcome irrespective of the calibration procedure. In
this  respect,  Figure  \ref{fig:res3}  also presents the comparison between the
damping  curve  by  Seed  and  Idriss  and  the model prediction. The frictional
$\zeta$  curve  lies in the same experimental range, even though the  accuracy at $\gamma_{max}=0.03-1\%$ is not as good as
for  the  $G/G_{max}$  ratio.   In  this case, the contribution of the
viscous  mechanism  is  practically non-existent,  as  it only increases the total
$\zeta$ ratio for $\gamma_{max}>0.1\%$.

Depending  on  the specific application, a ``trial and error'' calibration might
be    preferable,    sacrificing    some   of   the   accuracy   in   terms   of
$G/G_{max}$  ratio  to  improve  the  damping performance. A possible
outcome of a manual calibration is plotted in Figure \ref{fig:res4}: apparently,
while  the  simulation  of  the stiffness curve is still acceptable, the damping
curve  appears  to  be much better than the previous one. The use of the viscous
mechanism  seems to be highly beneficial, since it remedies the lack of accuracy
in the frictional curve at medium/large cyclic strains.

\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res4.pdf}
{res4.pdf}
 \hfill
 \caption{Comparison  between  experimental and simulated $G/G_{max}$
 and  damping  curves ($p_0$=100 kPa, T=2$\pi$ s, $\zeta$ = 0.003, $G_{max}$ = 4
 MPa, $\nu$=0.25, $M$=1.2, $k_d$=$\xi$=0, $h$=$G_{max}$/(15$p_0$), $m$=1)}
 \label{fig:res4}
\end{figure}

It  is   also worth noting that the experimental/numerical agreement is satisfactory  up
to even   $\gamma_{max}=10\%$, where substantial  plasticity  occurs and, in any case, the  extrapolation  of cyclic curves from experimental data is -- to say the  least -- questionable.  

Besides,   if  the  experimental  data  under
examination are unsatisfactorily reproduced for any $h$ and $m$ combination, the
user   still   has   the  chance  of  substituting  the  interpolation  function
\eqref{eq:hard_modulus} with no further changes in the model formulation. In particular, the present model is as flexible as the one by \cite{Borja_1994} in reproducing, for a given initial mean pressure, usual 1D non-linear laws for soils, such as the exponential, the hyperbolic, the Davidenkov and the Ramberg-Osgood models (for this latter a hardening bounding surface would be needed as well). Matching the aforementioned 1D laws ensures sufficient capability of reproducing experimental curves of usual  shape.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\section{Parametric analysis}

In  this section the influence of some relevant input
parameters on the model predictions is parametrically investigated.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Influence of the confining pressure}

Figure  \ref{fig:res5}  illustrates  the  sensitivity, under PS loading, of both
$G/G_{max}$  and  damping  frictional curves to the initial confining
pressure.   As   can   be   noticed,   increasing   $p_0$   does   enlarge   the
``pseudo-elastic''  range,  that  is  the  strain  interval  within  which  the
deviation  by  the  elastic  behavior  is negligible even with a vanishing yield
locus.  It  is  also  noted that the variations in the confining pressure do not
imply appreciable changes in the shape of the curves.

\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res7.pdf}
{res7.pdf}
 \hfill
 \caption{Simulated  $G/G_{max}$  and  damping  curves  at  varying
 confining  pressure  (T=2$\pi$  s,  $G_{max}$  =  4  MPa,  $\nu$=0.25, $M$=1.2,
 $k_d$=$\xi$=0, $h$=$G$/(15$p_0$), $m$=1)}
 \label{fig:res5}
\end{figure}

Although the present version of the model is apparently pressure-dependent, it cannot quantitatively reproduce the pressure-sensitiveness of $G/G_{max}$ and $\xi$ curves arising from real experiments and incorporated into some analytical formulas \citep{Ishibashi_Zhang_1993, Darendeli_2001}. This stems from the fact that the influence of the mean pressure only concerns the plastic component of the model, while constant elastic (and viscous) moduli have been adopted for the initial simulation of only PS loading tests. From this point of view, two easy improvements are possible and mutually compatible: 

\begin{enumerate}
\item use of hyper- or hypo-elastic laws with variable moduli \citep{Papadimitriou2002,andrianopoulos2010};
\item adoption of $p$-dependent hardening parameters (i.e. $h$ and $m$).
\end{enumerate}
In particular, the latter point does not introduce any further difficulty in terms of analytical/numerical treatment, since it only affects the interpolation rule \eqref{eq:hard_modulus}. Appropriate $h\left(p\right)$ and $m\left(p\right)$ relationships could be easily obtained by first calibrating $h$ and $m$  on  experimental or analytical cyclic curves for different $p$ values, and then analytically interpolating the parameters values over a meaningful pressure range.

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Influence of the hardening parameters}

Figures  \ref{fig:res6}  and  \ref{fig:res7} show the influence of the hardening
parameters  $h$  and  $m$  on  the  predicted cyclic curves. In particular, a decrease in either $h$  or  $m$ implies a faster development of plastic strains, so that the
closely-elastic range tends to disappear and $G/G_{max}<1$ and $\zeta>0$ at even  $\gamma_{max}=10^{-4}\%$; conversely, an extended pseudo-elastic behavior
can  be  obtained  over  a  large  strain  range  by  increasing  the  hardening
parameters.  Apparently,  the  model ensures high flexibility in terms of cyclic
curve  shapes,  so  that  the  response of standard elastic-plastic models (i.e. with
non-vanishing  elastic  region)  can  be  smoothly  approximated  (compare  for
instance   the   $m=3$  curves  in  Figure  \ref{fig:res7}  and  the  analytical
elastic-perfectly plastic frictional curves in Figure \ref{fig:damping_PP}).

\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res5.pdf}
{res5.pdf}
 \hfill
 \caption{Simulated  $G/G_{max}$ and  damping  curves  at  varying h
 ($p_0$=100   kPa,  $T$=2$\pi$  s,  $G_{max}$  =  4  MPa,  $\nu$=0.25,  $M$=1.2,
 $k_d$=$\xi$=0, $m$=1)}
 \label{fig:res6}
\end{figure}


\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res6.pdf}
{res6.pdf}
 \hfill
 \caption{Simulated  $G/G_{max}$  and  damping  curves  at  varying m
 ($p_0$=100   kPa,  $T$=2$\pi$  s,  $G_{max}$  =  4  MPa,  $\nu$=0.25,  $M$=1.2,
 $k_d$=$\xi$=0, $h$=$G_{max}$/(15$p_0$))}
 \label{fig:res7}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Influence of the viscous mechanism}

The   influence   of   the   viscous   parameter   $\zeta_0$  on  the  resulting
frictional/viscous  damping  curve  is illustrated in Figure \ref{fig:res9} (the
$G/G_{max}$  is  not  affected by the parallel viscous mechanism). As
was   expected,   an   increase   in   $\zeta_0$   induce   larger   values   of
$\zeta\left(\gamma_{max}\rightarrow0\right)$,  as  well  as a faster increase of
the  $\zeta$ curve at medium/high cyclic strains. Figure \ref{fig:res9} confirms
the usefulness of the viscous mechanism, which can be exploited as an additional degree of freedom for
reproducing the cyclic dissipative soil behavior.


\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.5\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res9.pdf}
{res9.pdf}
 \hfill
 \caption{Damping   curves   simulated  at  varying  $\zeta_0$  ($p_0$=100  kPa,
 $T$=2$\pi$   s,   $G_{max}$   =  4  MPa,  $\nu$=0.25,  $M$=1.2,  $k_d$=$\xi$=0,
 $h$=$G_{max}$/(15$p_0$), $m$=1)}
 \label{fig:res9}
\end{figure}

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Interaction between volumetric behavior and kinematic constraints}
%
All the above simulations have been performed by inhibiting the elastic-plastic soil
dilatancy      ($\xi=0$),    which   in most cases cannot be done to represent real soil
behavior. As previously shown for triaxial loading conditions (Figure
\ref{fig:res1}),  in  the absence of kinematic boundary constraints, a variation
in  the volumetric behavior slightly affects only the hardening evolution of the
stress-strain response toward the limit shear strength; a similar consideration
applies to PS loading conditions, since even in this case the normal confinement
is statically determined.

However,  computational  (FE)  models often contain kinematic constraints arising from
certain  symmetries  (consider e.g.  plane  strain or one-dimensional schemes)
\citep{Prevost_1989,     Borja99,  diPrisco2011}. In addition, for SSI problems,
where  soil  interacts  with a (stiff) structural foundations and wall, the soil
volume change plays an important role. 
The presence of kinematic constraints implies that the value of some stress components directly derives from compatibility requirements (e.g. prevented lateral expansion).
That  means that the local mean confinement is affected by the tendency
of  the  material  to  dilate  or  contract.  In particular, dilative frictional
materials  will  increase  the  limit  shear  stress (with respect to unconfined
conditions),     while  compactive frictional materials will  decrease the limit
shear stress.  Further,  not  only  the  limit  shear  stress,  but  also the whole
pre-failure  response  depends  on  the  plastic  flow  rule whenever kinematic
constraints are imposed \citep{diPrisco_Pisano_2011, diPrisco2011}.

The  above  considerations  suggest that both experimental and numerical results
are  certainly affected by the kinematics of the system, even though this effect
is  not  easy  to  be  \textit{a  priori}  quantified in terms of
$G/G_{max}$  and  $\zeta$  curves.  The  kinematic  conditions  of an
infinite  soil  layer  during  1D  shear  wave  propagation  are  experimentally
approximated   through   the   well known   ``simple   shear  (SS)  apparatus''
\citep{Wood2004},  in  which  the  soil  specimen  is cyclically sheared with no
lateral  expansion  allowed.  In  order  to assess how the kinematic confinement
influences  the  cyclic  response,  stiffness degradation and damping curves are
hereafter  simulated  under  SS conditions by varying the volumetric response of
the  soil;  in particular, three different calibrations of the plastic flow rule
\eqref{eq:flow_rule}  are  considered,  namely (i) isochoric ($k_d=\xi=0$), (ii)
compactive ($k_d=M$, $\xi=1$) and (iii) dilative ($k_d=0.4$, $\xi=1$)

The  results  reported  in  Figure  \ref{fig:res8}  provide  an insight into the
possible  effect  of  the  volumetric  response  in combination with constrained
loading  conditions.  In  the isochoric case, the PS and the SS curves perfectly
match (compare e.g. with the $p_0=$100 kPa curves in Figure \ref{fig:res5}), as,
with  no  plastic  expansion  (or  contraction),  the lateral constraints do not
affect  the  mean pressure during the shear loading; conversely, non-negligible
SS-PS differences arise when dilative or contractive materials are considered.
As  is  evident  in Figure \ref{fig:res8}, the discrepancy between isochoric and
non-isochoric curves becomes evident at medium/high cyclic strains, i.e. when  significant plastifications take place. Indeed, while the mechanical response is
barely  inelastic,  the  deviatoric and the volumetric responses are practically
decoupled, so that no variation of the normal confinement takes place.

Apparently, the quantitative relevance of this effect  strictly relates  to  the actual dilational properties of the material: soils undergoing significant volume changes under unconfined shear will exhibit a high sensitiveness to boundary constraints. The cyclic interaction between volumetric behavior and kinematic constraints seems to be poorly investigated in literature and is worth remarking for both theoretical  and  practical motivations. Indeed, the cyclic behavior measured through  certain experimental  devices (triaxial,  biaxial, simple shear,
torsional  shear,  etc) can differ from the mechanical response characterizing other kinematic conditions in boundary value problems, so that the employment of volume-insensitive models (such as the linear equivalent) may lead to inaccurate predictions. In particular, Equation \eqref{eq:hard_modulus} clearly shows that if any relevant $p$ variation  arises from the interaction between dilatancy and boundary constraints, a variation in the hardening modulus $H$ and, therefore, the resulting stiffness will also take place. As a consequence, an influence on the global (non-linear) dynamic response is expected in terms of both amplitude amplification and frequency content \citep{roten2013}. While this aspect is totally disregarded by most modeling approaches in GEE, further work is currently ongoing to quantitatively investigate the role of soil dilatancy in affecting the outcomes of seismic site response (elastic-plastic) analyses.
%
% 
\begin{figure}  [!htb]
 \centering
\includegraphics [width=0.8\textwidth] %{/home/jeremic/tex/works/Papers/2012/Bounding_Surface_Frictional_Model/Figures/res8.pdf}
{res8.pdf}
 \hfill
 \caption{$G/G_{max}$ and damping curves simulated under SS conditions and different volumetric responses (
 $p_0$=100 kPa, $T$=2$\pi$ s, $G_{max}$ = 4 MPa, $\nu$=0.25, $M$=1.2, $k_d$=[1.2, 0.4], $\xi$=[0,1], $h$=$G_{max}$/(15$p_0$), $m$=1)}
 \label{fig:res8}
\end{figure}
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\section{Concluding remarks}

An incremental  3D  visco-elastic-plastic  constitutive  model  was  developed  to  simulate stiffness degradation and damping in soils under cyclic/dynamic loading. The model is based on
an  effective-stress  formulation  with  two  parallel  dissipative mechanisms,
purely frictional (elastic-plastic) and viscous.

As  far as the frictional mechanism is concerned, a bounding surface formulation
with   vanishing  elastic region   was   adopted,   extending  to  
pressure-sensitive  non-associative  soils  the  previous  cohesive  model  by
\citet{Borja_1994} for total-stress analyses, but maintaining higher simplicity than later works, such as that e.g. by \cite{andrianopoulos2010}. The main features of the frictional
model are: (i) the vanishing yield locus implies an elastic-plastic response at
any  load levels, as is observed in real experiments; (ii) a minimum
number  of  physically  meaningful parameters, which can be easily calibrated on
  few  experimental  data;  (iii)  excellent  performance  and
flexibility  in  reproducing  in  the  elastic-plastic  framework  the  standard
stiffness  degradation  and  damping curves. With reference to these latter, the
parallel  viscous  mechanism  --  easy  to be introduced in FE computations -- was
shown  to  provide  an additional degree of freedom to improve the simulation of
the  cyclic  energy  dissipation,  as  long as the viscous parameter is properly
calibrated.  As  a  matter  of  fact,  the  viscous mechanism, used here,  does 
physically  exist  in  the  form  of  viscous interaction between the soil solid
skeleton  and  the  pore  fluid(s),  and  needs to be taken into account (as for
example done here).

Future  work  will concern the investigation of the model performance under undrained conditions and/or in the presence of non-symmetric loading. Further improvements will be possibly introduced, still  with purpose of keeping the model  simple and with few material parameters -- to be all calibrated on standard experimental data.



\section*{Acknowledgment}
Funding from and collaboration with the US NRC and funding from US DOE for this research is greatly appreciated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%
\bibliographystyle{chicago} 
%%\bibliographystyle{plain}
\bibliography{refmech,refFederico_small}
%\bibliography{} % These is my BiBTeX reference list.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\appendix
%
\section{Derivation of Equation \eqref{eq:hm_calibration}}
\label{app:A}

Under    PS    loading    conditions    (pure shear at constant   mean   pressure),   Equation
\eqref{eq:projection} produces  the following simple expression for $\beta$:

%
\begin{equation}
\label{eq_betaPS}
\beta = \dfrac{\tau_{lim} - \tau}{\tau - \tau_0}
\end{equation}
%
where  $\tau_{lim}=Mp_0/\sqrt{3}$ and $p_0$ is the initial (and constant) mean pressure. Equation \eqref{eq:hard_relation} can be easily specified for PS loading:
\begin{equation}
d{q}=\sqrt{3}d\tau= \sqrt{\dfrac{2}{3}}H\Vert d{e}_{ij}^p \Vert
= \sqrt{\dfrac{2}{3}}H\dfrac{d\gamma^p}{\sqrt{2}}
\Longrightarrow d\tau = \dfrac{H}{3}d\gamma^p
\end{equation}
so that the following form of the PS elastic-plastic response results:
%
\begin{equation}
d\gamma=d\gamma^e+d\gamma^p=\dfrac{d\tau}{G_{max}}+\dfrac{3d\tau}{H}
\end{equation}
%
and, after substituting \eqref{eq_betaPS}  into \eqref{eq:hard_modulus}:
\begin{equation}
d\gamma=\dfrac{d\tau}{G_{max}}+\dfrac{3d\tau}{p_0h\beta^m}=
\dfrac{d\tau}{G_{max}}+\left(\dfrac{\tau -
\tau_0}{\tau_{lim} - \tau}\right)^m\dfrac{3d\tau}{p_0h}
\end{equation}
%
Integration  over  a  strain  interval  between  two  stress  reversals ($\gamma
\in\left[-\gamma_{max};\gamma_{max}\right]$) yields:
%
\begin{equation}
\label{eq:integration_app}
2\gamma_{max}   =   \dfrac{2\tau}{G_{max}}+\dfrac{3}{hp_0}   \int_{-\tau}^{\tau}
\left(\dfrac{\tau' + \tau}{\tau_{lim} - \tau'}\right)^md\tau'
\end{equation}
%
where $\tau_0=-\tau$ has been set. Straightforward variable changes lead to:
%
\begin{equation}
1  = \dfrac{G}{G_{max}} + \dfrac{3}{2hp_0\gamma_{max}} \int_{0}^{2G\gamma_{max}}
\left(\dfrac{\tau''}{\tau_{lim}-\tau''+ G\gamma_{max}}\right)^md\tau''
\end{equation}
%


%
\begin{equation}
\label{eq:calibration_app}
1    =    \dfrac{G}{G_{max}}\left[    1   +   \dfrac{6G_{max}}{hp_0\gamma_{max}}
\int_{0}^{\gamma_{max}}
\left(\dfrac{\gamma}{\tau_{lim}/G-2\gamma+\gamma_{max}}\right)^md\gamma\right]
\end{equation}\\
%

It  is  worth  highlighting  that two approximations are implicitly contained in
Equation  \eqref{eq:calibration_app}: (i) the integration over the first loading
cycle  does not exactly reproduce the stabilized cyclic response (because of the
aforementioned  Masing effect); (ii) a symmetric loading cycle in terms of shear
strain does not in general ensure the symmetry of the corresponding shear stress
range  (as  it is assumed in Equation \eqref{eq:integration_app}). However, such
approximations do not prevent reasonable values for the hardening parameters $h$
and $m$ to be obtained.








\end{document}
