% \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}