\documentclass[dvips]{nagauth}
\usepackage{graphicx}
\usepackage{lscape}
%\usepackage[authoryear]{natbib}
\usepackage[numbers,round,colon,authoryear]{natbib}
%\usepackage[numbers,round,colon]{natbib}
\usepackage{amsmath}
\usepackage{ulem}
\usepackage{url}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\newcommand{\ud}{{\rm d}}
\def\baselinestretch{2}
\begin{document}
\NAG{1}{6}{01}{28}{01}
\runningheads{Jeremi{\'c}, Cheng, Taiebat and Dafalias}
{Numerical Simulation of Fully Saturated Porous Materials}
\title{Numerical Simulation of Fully Saturated Porous Materials}
\author{
Boris Jeremi{\'c}
(\raisebox{1.0mm}{\includegraphics[width=2.7cm]{BorisJeremic.eps}})\affilnum{1}
\comma\corrauth \\
%%%%%%%%
Zhao Cheng
(\raisebox{1.2mm}{\includegraphics[width=1.1cm]{Zhao.eps}})\affilnum{2} \\
%%%%%%%%
Mahdi Taiebat
(\raisebox{1.5mm}{\includegraphics[width=2.0cm]{mahdi.eps}})\affilnum{3} \\
%%%%%%
and
%%%%%%%%
Yannis Dafalias
(\raisebox{1.0mm}{\includegraphics[width=2.7cm]{yannis.eps}})\affilnum{4}
}
\address{
\affilnum{1}\ {Associate Professor,
Department of Civil and Environmental Engineering, University of
California, Davis, CA 95616. \texttt{jeremic@ucdavis.edu}} \\
%
\affilnum{2}\ {Staff Engineer, Earth Mechanics Inc. Oakland, CA 94621}\\
%
\affilnum{3}\ {Graduate Student,
Department of Civil and Environmental Engineering, University of
California, Davis, CA 95616.}\\
%
\affilnum{4}\ {Professor,
Department of Civil and Environmental Engineering, University of
California, Davis, CA 95616 and Department of Mechanics,
National Technical University of Athens, Zographou 15780, Greece.}
}
\corraddr{ Boris Jeremi{\'c}, Department of Civil and Environmental Engineering,
University of California,
One Shields Ave., Davis, CA 95616, \texttt{jeremic@ucdavis.edu}}
\cgsn{NSF}{EEC9701568}
\cgsn{NSF}{CMS0337811}
\noreceived{}
\norevised{}
\noaccepted{}
\newpage
\paragraph{Abstract}
Fully coupled, porous solidfluid formulation, implementation and related
modeling and simulation issues are presented in this work. To this end, coupled
dynamic field equations with $upU$ formulation are used to simulate pore fluid
and soil skeleton (elasticplastic porous solid) responses. Present formulation
allows, among other features, for water accelerations to be taken into account.
This proves useful in modeling dynamic interaction of media of different
stiffness (as in soilfoundationstructure interaction). Fluid compressibility
is also explicitly taken into account, thus allowing excursions into modeling of
limited cases of nonsaturated porous media. In addition to those feature,
present formulation and implementation models in a realistic way the physical
damping, which dissipates energy. In particular, the velocity proportional
damping is appropriately modeled and simulated by taking into account
interaction of pore fluid and solid skeleton. Similarly, the displacement
proportional damping is physically modeled trough elasticplastic processes in
soil skeleton. An advanced material model for sand is used in presented work and
is discussed at some length. Also explored in this paper are verification and
validation issues related to fully coupled modeling and simulations of porous
media.
Illustrative examples describing dynamical behavior of porous media (saturated
soils) are presented. The verified and validated methods and material models are
used to {\it predict} behavior of level and sloping grounds subjected to seismic
shaking.
\keywords{Elastoplastic,
Finite elements,
Saturated porous solidfluid simulations,
Soil dynamics}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{Background}
Modeling and simulations of soil behavior induced by
earthquakes, that can be
related to liquefaction and cyclic mobility,
continue to challenge engineering research and practice. Such behavior
commonly occurs
in loose to medium dense sands which are fully saturated and will results the almost
complete loss of strength (liquefaction), and partial loss of strength
(cyclic mobility).
%
To model these complex phenomena, a consistent and efficient coupled
formulation must be utilized, including an accurate singlephase
constitutive model for soil.
%
Three general continuum formulations \citep{Zienkiewicz84} are possible for modeling of
the fully coupled problem (soil
skeleton  pore fluid) in geomechanics, namely the (a) $up$, (b) $uU$, and
(c) $upU$ formulations. Here, the unknowns are the soil skeleton displacements
$u$; the pore fluid (water) pressure $p$; and the pore fluid (water)
displacements $U$. The $up$ formulation captures the
movements of the soil skeleton and the change of the pore pressure, and is the most
simplistic one of the three mentioned above. This formulation neglects the
accelerations of the pore fluid (except for combined (same) acceleration of pore fluid
and solid),
and in one version neglects the compressibility
of the fluid (assuming complete incompressibility of the pore fluid). In the
case of incompressible pore fluid and very low permeability, the formulation requires
special treatment of the approximation (shape) function for pore fluid to
prevent the volumetric locking \citep{Ziekienwicz2000a}. This formulation also
must rely on Rayleigh damping to model velocity proportional energy dissipation
(damping).
The majority of the currently available implementations are based on this
formulation. For example \citet{Elgamal2002} and \citet{Elgamal2003}
developed an implementation of the $up$ formulation with the multisurface
plasticity model of \citet{Prevost85}, while \citet{Chan88} and
\citet{Zienkiewicz99}
used generalized theory of plasticity by \citet{Pastor90} in their simulations.
In addition to that \citet{Taiebat2007} used the constitutive model by
\citet{Manzari97} in their $up$
formulation and showed good model performance for capturing the response in
a boundary value problem from
VELACS project \cite{Arulanandan93}.
%
We also note earlier works by \cite{Zienkiewicz1980}, \cite{Simon1984,
Simon1986}, \cite{Sandhu1990} and \cite{Gajo1994}.
%
%
The $uU$ formulations tracks the movements of both the soil skeleton and the
pore fluid. This formulation is complete in the sense of basic variable, but
might still experience numerical problems (volumetric locking) if the difference
in volumetric compressibility of the pore fluid and the solid skeleton is large.
The $upU$ formulation resolves the issues of volumetric locking by including
the
displacements of both the solid skeleton and the pore fluid, and the pore
fluid pressure as well.
This formulation uses additional (dependent) unknown field of pore fluid
pressures to stabilize the solution of the coupled system. The pore fluid
pressures are connected to (dependent on) displacements of pore fluid, as, with known
volumetric compressibility of the pore fluid, pressure can be calculated.
All three formulations were originally derived by \citet{Zienkiewicz84} and are
critically discussed in \citet{Zienkiewicz99}.
% %
Despite it's power, the $upU$ formulation has rarely been implemented into
finite element code.
%
In this paper, complete $upU$ formulation and implementation is presented.
A set of verification examples are used to determine that formulation and the
implementation accurately represent conceptual description of porous media
behavior\footnote{Verification is a mathematics issue. It provides evidence that
the model is solved correctly.}.
%
In
addition to that, a recent critical state twosurface plasticity model
accounting for the fabric dilation effects \citep{Dafalias04}
is described in some details as well.
The material model behavior is validated using a number of constitutive tests from
literature. The validation process is important in determining the degree to which a
model is accurate representation of the real world from the perspective of
the intended uses of the model\footnote{Validation is a physics issue.
It provides
evidence that the correct model is solved.}.
Issues of verification and validation are fundamental to the process of
prediction of mechanical behavior in computational science and engineering
\citep{Oberkampf2002,Roach1998}.
A number of predictive examples are presented in a later part of the paper. The
examples rely on verified formulation and implementation of
behavior of fully coupled porous media and on validated
constitutive material modeling.
%
Predictive examples consist of level and sloping grounds excited by seismic
shaking at the base. Examples explore behavior (displacements of solid
skeleton and pore fluid, excess pore pressure, stressstrain) for two
different types (dense and loose) of Toyoura sand.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Formulation, Finite Element Discretization and Implementation}
This section briefly reviews equations describing behavior of solid skeleton and
pore fluid and their interaction.
%
It is important to note the sign convention used in developments described in
this paper.
The compression stress is assumed negative here (as is common in mechanics of
materials), and the effect of this sign convention has been considered on the
model equations. For example this is the reason for having couple of different signs in
the equations comparing to the original work by \citet{Dafalias04}.
%
The mean effective stress $p$ is defined as $p=\sigma^{'}_{ii}/3 $ and
therefore is positive in compression.
The deviatoric stress tensor $s_{ij}$ then can be obtained using
$s_{ij}=\sigma^{'}_{ij}+p\delta_{ij}$.
%
As for the (small) strain tensor $\epsilon_{ij} = (u_{i,j}+u_{j,i})/2$,
the volumetric component is defined as
$\epsilon_{v}=\epsilon_{ii}$, while the deviatoric
strain tensor $e_{ij}$ can be obtained from
$e_{ij}=\epsilon_{ij}\epsilon_{v}\delta_{ij}/3$.
\subsection{Governing Equations of Porous Media}
\label{Governing Equations of Porous Media}
The material behavior of soil skeleton is fully dependent on the pore fluid
pressures. The behavior of pore fluid is assumed to be elastic and thus all the
material nonlinearity is concentrated in the soil skeleton.
The soil behavior (mix of soil skeleton and pore fluid) can thus be described
using singlephase constitutive analysis approach for skeleton combined with the full
coupling with pore fluid.
The concept of effective stress of the saturated mixture, that is,
the relationship between effective stress,
total stress and pore pressure \citep{Zienkiewicz99} can be written as
$\sigma^{'}_{ij}= \sigma_{ij}+\alpha \delta_{ij} p$,
where $\sigma^{'}_{ij}$ is the effective stress tensor,
$\sigma_{ij}$ is total stress tensor, $\delta_{ij}$ is Kronecker delta and
$\alpha$ is the Biot constant that depends on the geometry of material voids.
For the most part, in soil mechanics problems,
$\alpha \approx 1$ can be assumed.
The relation between total and effective stress becomes
$\sigma^{'}_{ij} = \sigma_{ij}+\delta_{ij} p$,
which corresponds to the classical effective stress definition by \cite{Terzaghi43}.
The overall equilibrium or momentum balance equation
for the soilfluid 'mixture' can be written as
%%%
\begin{equation}
\sigma_{ij,j}  \rho \ddot{u}_i  \rho_f \ddot{w}_i + \rho b_i = 0
\label{upU_GEequilibrium}
\end{equation}
%%%
where $\ddot{u}_i$ is the acceleration of the solid part, $b_i$ is the body
force per unit mass, $\ddot{w}_i$ is the fluid acceleration relative to the
solid part.
For fully saturated porous media (no air trapped inside),
density is equal to $\rho = n \rho_f + (1  n) \rho_s$,
where $n$ is the porosity,
$\rho_s$ and $\rho_f$ are the soil particle and water density respectively.
For the pore fluid, the equation of momentum balance can be written as
%%%
\begin{eqnarray}
p_{,i}  {R_i}  \rho_f \ddot{w}_i  \frac{\rho_f \ddot{w}_i}{n} +\rho_f b_i = 0
\label{GEequilibriumFluid}
\end{eqnarray}
%%%
where $R$ is the viscous drag force. According to the Darcy's seepage law,
the viscous drag forces $R$ between soil matrix and pore fluid (water) can be written as
${R_i} = k_{ij}^{1} \dot{w}_j$,
where
$k_{ij}$ is the tensor of anisotropic Darcy permeability coefficients. For
simple case of isotropic permeability, scalar value of permeability $k$ is used.
The permeability ${k}$ used here with dimension of [$L^{3}TM^{1}$] is different from the permeability
used in the usual soil
mechanics ($K$) which has the same dimension of velocity, i.e. [$LT^{1}$].
Their values are related by $k = K / g \rho_f$,
where $g$ is the gravitational acceleration and the permeability $K$ is
measured in an experiment.
The final equation is the mass conservation of the fluid flow expressed by
%%
\begin{eqnarray}
\dot{w}_{i,i} + \alpha \dot{\varepsilon}_{ii} + \frac{\dot{p}}{Q} = 0
\label{GEconservationFlow}
\end{eqnarray}
%%
where bulk stiffness of the mixture $Q$ is expressed as
${1}/{Q} = {n}/{K_f} + {(\alphan)}/{K_s}$
and $K_s$ and $K_f$ are the bulk moduli of the solid and fluid phases respectively.
In the above governing equations,
convective and terms of lower order are omitted \citep{Zienkiewicz99}.
A change of variable is performed by introducing an alternative variable $U_i$, defined as
$ U_i = u_i + U_i^R = u_i + {w_i}/{n}$, that represents absolute displacement
of the pore fluid.
The basic set of unknowns is then comprised of the soil skeleton displacements $u_i$,
the water pore pressure $p$,
and the water displacements $U_i$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Finite Element Formulation}
\label{Finite Element Formulation}
Standard finite element discretization \citep{FEMIV1,FEMIV2} uses shape
functions to describe
each of the unknown fields (upU) in terms of nodal values
(solid displacements $\bar{u}_{Ki}$, pore fluid pressures $\bar{p}_K$, and pore
fluid displacements $\bar{U}_{Ki}$):
%
\begin{equation}
u_i = N_K^u \bar{u}_{Ki}, \;\;\;
p = N_K^p \bar{p}_K, \;\;\;
U_i = N_K^U \bar{U}_{Ki}
\label{39}
\end{equation}
where $N_K^u$, $N_K^p$ and $N_K^U$ are (same) shape functions for solid displacement,
pore pressure and fluid displacement respectively. Each node of the (upU)
element has thus seven degrees of freedoms in 3 dimensions (three for solid displacements, one
for pore fluid pressures and three or pore fluid displacements). It should be
noted that it is possible to use same shape
functions for both displacement and pore pressure unknown field as upU
formulation with compressible fluid allows that without volumetric locking
\citep{Zienkiewicz84}.
By using finite element discretization and after some tensor algebra and manipulations,
the weak form of governing equations can be obtained from the strong form
described by Equations~(\ref{upU_GEequilibrium}), (\ref{GEequilibriumFluid}) and
(\ref{GEconservationFlow}).
% %%
The weak form is presented in tensor notation (see \cite{FEMIV1}
chapter 6.) as
%%
%
%
%
%
\begin{eqnarray}
& &
\left[ \begin{array}{ccc}
(M_s)_{KijL} & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & (M_f)_{KijL}
\end{array} \right]
\left[ \begin{array}{c}
\ddot{\bar{u}}_{Lj} \\
\ddot{\bar{p}}_N \\
\ddot{\bar{U}}_{Lj}
\end{array} \right]
+
\left[ \begin{array}{ccc}
(C_1)_{KijL} & 0 & (C_2)_{KijL} \\
0 & 0 & 0 \\
(C_2)_{LjiK} & 0 & (C_3)_{KijL} \\
\end{array} \right]
\left[ \begin{array}{c}
\dot{\bar{u}}_{Lj} \\
\dot{\bar{p}}_N \\
\dot{\bar{U}}_{Lj}
\end{array} \right] \nonumber \\
&+&
\left[ \begin{array}{ccc}
0 & (G_1)_{KiM} & 0 \\
(G_1)_{LjM} & P_{MN} & (G_2)_{LjM} \\
0 & (G_2)_{KiL} & 0
\end{array} \right]
\left[ \begin{array}{c}
\bar{u}_{Lj} \\
\bar{p}_M \\
\bar{U}_{Lj}
\end{array} \right] \nonumber \\
&+&
\left[ \begin{array}{c}
\int_{\Omega} N_{K,j}^u \sigma^{'}_{ij} \ud \Omega \\
0 \\
0
\end{array} \right]

\left[ \begin{array}{c}
\bar{f}^u_{Ki} \\
0 \\
\bar{f}^U_{Ki}
\end{array} \right]
= 0
\label{68}
\end{eqnarray}
%%
The left hand side components of the above matrix equation are given as: %%
% %%%
%%%
\begin{eqnarray}
(M_s)_{KijL}= \int_{\Omega} N_K^u (1n) \rho_s \delta_{ij} N_L^u \ud \Omega
\;\;\;\; &\mbox{;}& \;\;\;\;
(M_f)_{KijL}= \int_{\Omega} N_K^U n \rho_f \delta_{ij} N_L^U \ud \Omega \nonumber\\
(C_1)_{KijL}= \int_{\Omega} N_K^u n^2 k_{ij}^{1} N_L^u \ud \Omega
\;\;\;\; &\mbox{;}& \;\;\;\;
(C_2)_{KijL}= \int_{\Omega} N_K^u n^2 k_{ij}^{1} N_L^U \ud \Omega \nonumber\\
(C_3)_{KijL}= \int_{\Omega} N_K^U n^2 k_{ij}^{1} N_L^U \ud \Omega
\;\;\;\; &\mbox{;}& \;\;\;\;
(G_1)_{KiM} = \int_{\Omega} N_{K,i}^u (\alphan) N_M^p \ud \Omega \nonumber\\
(G_2)_{KiM} = \int_{\Omega} n N_{K,i}^U N_M^p \ud \Omega
\;\;\;\; &\mbox{;}& \;\;\;\;
P_{NM} = \int_{\Omega} N_N^p \frac{1}{Q} N_M^p \ud \Omega
\label{691}
\end{eqnarray}
%%%
%
%
while the right hand side components are given as:
% %%%
\begin{eqnarray}
(\bar{f}_s)_{Ki} &=& \int_{\Gamma_t} N_K^u \sigma^{'}_{ij} n_j \ud \Gamma

\int_{\Gamma_p} N_K^u (\alphan) p n_i \ud \Gamma
+
\int_{\Omega} N_K^u (1n) \rho_s b_i \ud \Omega \nonumber\\
(\bar{f}_f)_{Ki} &=&  \int_{\Gamma_p} N_K^U n p n_i \ud \Gamma
+
\int_{\Omega} N_K^U n \rho_f b_i \ud \Omega
\label{692}
\end{eqnarray}
It is very important to note that the velocity proportional damping was
introduced directly through the damping tensor (or matrix) with components
$(C_1)_{KijL}$, $(C_2)_{KijL}$ and $(C_3)_{KijL}$, which are functions of
porosity and permeability of the skeleton. This damping provides for physically
based energy dissipation due to interaction of pore fluid and the solid (soil)
skeleton. It is also emphasized again that presented formulation and
implementation do not use Rayleigh damping.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Constitutive Integration}
\label{Constitutive_Integration}
In work presented, an explicit version of constitutive driver was used for
DafaliasManzari model. While the implicit constitutive integrations
\citep{Jeremic96b} are preferred for local interactions, some
material models have a highly nonlinear regions in the yield surface,
plastic flow directions and/or evolution laws.
The DafaliasManzari material model \citep{Dafalias04}
used in this study, exhibits such
behavior, particularly in the very low mean effective
stress region, which is of much interest here.
Implicit constitutive integration
in such regions are highly unstable as the stress and internal variable state
is far from the so called trust region of the underlying Newton method (see
more in \cite{Schnabel83}). One possible solution for every
failed implicit step is to try the line search algorithm (\cite{Jeremic2001}).
However, while line search (eventually) automatically achieves solution, it
also tends to create very short steps. This stems from the fact that line
search keeps cutting in half Newton steps that fail to converge. In addition
to that, implicit integrations feature much higher computational load per step
compared to explicit integrations, leading to slower computations. Eventual
benefit of implicit integration in faster global Newton solution, was not
pursued in this work but will be investigated in future publications.
%
All of the above reasons led us to use explicit integration algorithm in this
work.
The Explicit (forward Euler) method is based on starting point of plastic flow
in the stress and internal variable space
for finding all the relevant derivatives. The increment of
stress and internal variables is given by
%
\begin{equation}
\Delta \sigma_{mn}^{'} =
{}^{c} E_{mnpq} \; \Delta \epsilon_{pq}^{ep} 
{}^{c} E_{mnpq} \; \frac{ {}^{c} n_{rs} \; {}^{c} E_{rstu} \; \Delta \epsilon_{tu}^{ep}}
{ {}^{c} n_{ab}\; {}^{c} E_{abcd} \; {}^{c} m_{cd}

{}^{n} \xi_{A} {}^{n} h_{A} } \;
{}^{c} m_{pq}
\label{predictor10}
\end{equation}
%
\begin{eqnarray}
\Delta q_{A} =
\left(
\frac{ {}^{c} n_{mn} \; {}^{c} E_{mnpq} \; \Delta \epsilon_{pq}^{ep}}
{\; {}^{c} n_{mn}\; {}^{c} E_{mnpq} \; {}^{c} m_{pq}
 {}^{n} \xi_A {}^{n} h_A }
\right)
h_A
\label{FE_predictor1}
\end{eqnarray}
%%
It is important to note that the explicit algorithm performs only one step of the
computation and does not check on the convergence of the provided solutions.
This usually results in the slow drift of the stressinternal variable point
from the yield surface for monotonic loading. The use of Explicit
integration can also result in spurious elasticplastic
deformations during elastic
unloading in cyclic loadingunloading \citep{Jeremic2000f}. In addition
to that, the Explicit integrations can lead
to completely erroneous step if the load reversal step is large enough that
the stress path misses the elastic region and lands in elasticplastic region
on the opposite side of yield surface.
In this case, derivatives of yield surface, plastic flow directions and
derivatives of evolution law from the starting point will point one way, while
the actual starting point should be on the other side of the yield surface.
Both spurious elasticplastic step and inconsistency in starting point need to
be carefully watched for during computations.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Time Integration}
\label{Time_Integration}
In order to develop integration of dynamic finite element equation in the time
domain, Equation~(\ref{68}) is rewritten in a residual matrix form \citep{local87}
%%%
\begin{eqnarray}
\bf{ R = M \ddot{x} + C \dot{x} + K' x + F(x)  f = 0 }
\label{71}
\end{eqnarray}
%%%
where $\bf{x = \{\bar{u}, \bar{p}, \bar{U}\}^T}$ represent a vector of
generalized unknown variables.
Equation~(\ref{71}) represents the general nonlinear form \citep{local87} for which
the usual tangent stiffness $\bf K$ does not equal to $\bf K'$,
but instead,
\begin{eqnarray}
\bf{ K = \frac{\partial R}{\partial x} = K' + \frac{\partial F(x)}{\partial x} }
\label{72}
\end{eqnarray}
In the specific case of $upU$ formulation of interest here one can write
matrix form
%%%
\begin{eqnarray}
\bf{ M =
\left[ \begin{array}{ccc}
\bf{M_s} & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & \bf{M_f}
\end{array} \right] }, \;\;\;
\bf{ C =
\left[ \begin{array}{ccc}
\bf{C_1} & 0 & \bf{C_2} \\
0 & 0 & 0 \\
\bf{C_2^T} & 0 & \bf{C_3} \\
\end{array} \right] }
\label{upU_MC}
\end{eqnarray}
%%
\begin{eqnarray}
\bf{ K' =
\left[ \begin{array}{ccc}
0 & \bf{G_1} & 0 \\
\bf{G_1^T} & \bf{P} & \bf{G_2^T}\\
0 & \bf{G_2} & 0
\end{array} \right]},\;\;\;
\bf{ K =
\left[ \begin{array}{ccc}
\bf{K^{ep}} & \bf{G_1} & 0 \\
\bf{G_1^T} & \bf{P} & \bf{G_2^T}\\
0 & \bf{G_2} & 0
\end{array} \right] }
\label{upU_K}
\end{eqnarray}
where
\begin{eqnarray}
{\bf K^{ep}} = (K^{ep})_{KijL} = \int_{\Omega} N_{K,m}^u E^{ep}_{imjn} N_{L,n}^u \ud \Omega
\label{upU_Kep}
\end{eqnarray}
The above set of residual (nonlinear) dynamic equations is solved using,
usually, one of the two procedure, namely the
Newmark \citep{Newmark1959} and the
HilberHughesTaylor (HHT) $\alpha$method
\citep{Hilber77,Hughes78a,Hughes78b}.
%
\paragraph{Newmark Integrator.}
The Newmark time integration method \citep{Newmark1959} has two parameters,
$\beta$ and $\gamma$, and is described by the following two equations:
%%%%
\begin{eqnarray}
{}^{n+1} x &=& {}^n x + \Delta t \; {}^n \dot{x} + \Delta^2 t \; [(\frac{1}{2}  \beta) \; {}^n \ddot{x} + \beta \; {}^{n+1} \ddot{x}]
\label{upU_Newmark_EQ1} \\
{}^{n+1} \dot{x} &=& {}^n \dot{x} + \Delta t \; [(1  \gamma) \; {}^n \ddot{x} + \gamma \; {}^{n+1} \ddot{x}]
\label{upU_Newmark_EQ2}
\end{eqnarray}
%%%%
which give the relations between the time step $n$ to the next time step $n+1$.
The method is an implicit, except when both $\beta$ and $\gamma$ are zero.
If the parameters $\beta$ and $\gamma$ satisfy the following conditions
%
\begin{equation}
\gamma \ge \frac{1}{2}, \;\;\; \beta = \frac{1}{4}(\gamma+\frac{1}{2})^2
\label{upU_Newmark_ab}
\end{equation}
%
the time integration method is unconditionally stable. Any $\gamma$ value
greater than $0.5$ will introduce numerical damping. Wellknown members of the
Newmark time integration method family include: trapezoidal rule or average
acceleration method for $\beta = 1/4$ and $\gamma = 1/2$, linear acceleration
method for $\beta=1/6$ and $\gamma = 1/2$, and (explicit) central difference
method for $\beta = 0$ and $\gamma = 1/2$. If and only if $\gamma = 1/2$, the
accuracy is secondorder \citep{local86}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{HilberHughesTaylor Integrator.}
The HHT time integration method uses an alternative residual
form
%(alternative to Newmark method, see \cite{local87})
by introducing an addition parameter $\alpha$ to improve the performance:
%%%%
\begin{eqnarray}
{}^{n+1} R = M \; {}^{n+1} \ddot{x}
+ (1 + \alpha) F({}^{n+1} \dot{x}, {}^{n+1} x)
 \alpha F({}^{n} \dot{x}, {}^{n} x)
 \; {}^{n+1} f
\label{upU_HHT_R}
\end{eqnarray}
%%%%
while retaining the rest of Newmark equations
(\ref{upU_Newmark_EQ1}) and (\ref{upU_Newmark_EQ2}) and its parameters
$\beta$ and $\gamma$.
%
If the parameters $\alpha$, $\beta$ and $\gamma$ satisfy
%%%%
\begin{equation}
1/3 \le \alpha \le 0, \;\;\; \gamma = \frac{1}{2}(12\alpha), \;\;\; \beta = \frac{1}{4}(1\alpha)^2
\label{upU_HHT_ab}
\end{equation}
%%%
the HHT method is unconditionally stable and secondorder accurate \citep{local86}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Program Implementation}
\label{PI}
Implementation of the described algorithms and procedures was performed using a
number of numerical libraries. Parts of OpenSees framework
\citep{McKenna97} were used to connect the finite element domain. In
particular, Finite Element Model Classes from OpenSees (namely, class
abstractions for Node, Element, Constraint, Load and Domain) where used to describe
the finite element model and to store the results of the analysis
performed on the model. In addition to these, Analysis Classes were used to
drive the global level finite element analysis, i.e., to form and
solve the global system of equations.
%
As for the Geomechanics modules, a number of elements, algorithms and material models from UCD
Computational Geomechanics toolset
are used. In particular, set of NewTemplate3Dep numerical libraries were
used for constitutive level integrations, nDarray numerical libraries
\citep{Jeremic97d} were used to handle vector, matrix
and tensor manipulations, while FEMtools element libraries
were used coupled finite elements
(upU) and for element level computations. Finally, solvers from the uMfPACK
set of libraries \citep{Davis97b} were used to solve the
nonsymmetric global (finite element level) system of equations.
All of the above libraries are available either through their original
developers or through first Author's web site
\url{http://geomechanics.ucdavis.edu}.
\vspace*{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Material Model}
\label{Material_Model}
The constitutive model plays a very important role in numerical simulation of the
soil response. Within the critical state soil mechanics framework,
\cite{Manzari97} proposed a twosurface critical state plasticity model for
sands. Employing the effects of the state parameter on the behavior of sand,
this model presents wellestablished mechanisms for prediction of
softening/hardening as well as dilatancy/contractancy response in different
loose and dense states of sand. \cite{Dafalias04} later presented an
improved version of the model. This version introduced the fabricdilatancy
tensor which has a significant effect on the contraction unloading response.
%
This version also considered the effect of a modified Lode angle on the flow rule,
which produces more realistic responses in nontriaxial conditions.
%
For
the sake of completeness of the paper, the essential elements of this
plasticity model are summarized and presented here with index notation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Critical State Line}
The constitutive model proposed by \cite{Dafalias04} is essentially based on
the critical state soil mechanics framework and the defining the correct
position of the CSL is important in this model. Instead of using the common
linear relation of the critical void ratio $e_{c}$ v.s. the logarithm of the
critical mean effective stress $p_{c}$, they have adopted a power relation,
suggested by \cite{Li98}, which gives a greater freedom in representing the
Critical State Line:
\begin{equation}
e_{c}=e_{c,r}\lambda_{c}\left( \frac{p_{c}}{p_{at}}\right) ^{\xi}
\label{DM04CS}%
\end{equation}
where $e_{c,r}$ is the reference critical void ratio, and $\lambda_{c}$ and
$\xi$ are two other material constants (for most sands $\xi=0.7$) while $p_{at}$
refers to the atmospheric pressure used for normalization.
The `state parameter' $\psi=ee_{c}$ can now be used as a measure of how far
the material state $e$, $p$ is from the critical state.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Elasticity}
The isotropic hypoelasticity assumption is adopted in this model and is
defined by
%
\begin{equation}
\dot{e}_{ij}^{e}=\frac{\dot{s}_{ij}}{2G},\;\;\dot{\epsilon}_{v}^{e}%
=\frac{\dot{p}}{K}%
\end{equation}
The elastic shear modulus $G$ is adopted from \cite{Richart70}. Assuming a
constant Poisson's ratio $\nu$, the elastic bulk modulus $K$ can be found from
the corresponding elastic shear modulus
%
\begin{equation}
G=G_{0}\dfrac{(2.97e)^{2}}{(1+e)}\left( \dfrac{p}{p_{at}}\right)
^{1/2}p_{at}
\;\;\;\; \mbox{;} \;\;\;\;
K=\dfrac{2(1+\nu)}{3(12\nu)}G
\end{equation}
%
where $G_{0}$ is a material constant, $e$ is the void ratio, and $p_{at}$ is
the atmospheric pressure used for normalization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Yield Function}
The yield function is defined by
%
\begin{equation}
f=[(s_{ij}p\alpha_{ij})(s_{ij}p\alpha_{ij})]^{1/2}\sqrt{\frac{2}{3}}mp=0
\label{DM04YF}%
\end{equation}
%
where $\alpha_{ij}$ is the deviatoric back stressratio tensor, and $m$ is a
material constant. The above equation describes geometrically a `cone' in the
multiaxial stress space. The trace of the cone on the stress ratio $\pi$plane
is a circular deviatoric yield surface with center $\alpha_{ij}$ and radius
$\sqrt{2/3}m$ as shown in Figure (\ref{DM_Pi_Plane}).
%
The cone type yield surface implies that only changes of the stress ratio can
cause plastic deformations. In a more recent version of this family of
models,
\cite{Taiebat2007a} have introduced a closed yield surface in the form
of a modified eightcurve equation with proper hardening mechanism for capturing
the plastic strain under constant stress to circumvent the mentioned limitation.
In order to incorporate
the effect of the third invariant of the stress tensor in the model, an
effective Lode angle $\theta$ is defined using
`unit' gradient tensor to the yield surface on the deviatoric plane,
$n_{ij}$, as
%
\begin{equation}
\cos{3\theta}=\sqrt{6}n_{ij}n_{jk}n_{ki}
\;\;\;\; \mbox{;} \;\;\;\;
n_{ij}=\frac{s_{ij}p\alpha_{ij}}{[(s_{ij}p\alpha_{ij})(s_{ij}p\alpha
_{ij})]^{1/2}}
\label{DM04Lode}%
\end{equation}
%
and $0\leq\theta\leq\pi/3$. This equation results in $\theta=0$ for triaxial
compression and $\theta_{n}=\pi/3$ for triaxial extension.
The critical stress ratio $M$ at any effective Lode angle $\theta$ can be
interpolated between $M_{c}$ and $M_{e}$, i.e. the triaxial compression and
extension critical stress ratios, respectively.
%
\begin{equation}
M=M_{c}g(\theta,c)
\;\;\;\; \mbox{;} \;\;\;\;
g(\theta,c)=\frac{2c}{(1+c)(1c)\cos{3\theta}%
}
\;\;\;\; \mbox{;} \;\;\;\;
c=\frac{M_{e}}{M_{c}}\label{DM04Argris}%
\end{equation}
The model employs two everchanging surfaces with the state parameter $\psi$,
namely bounding and dilatancy surfaces, in order to account for the
hardening/softening and dilatancy/contractancy response of sand. Upon shearing
these two surfaces change (move) toward a fixed critical state surface to make the
model fully compatible with critical state soil mechanics requirement. The
line from the origin of the stress ratio $\pi$plane parallel to $n_{ij}$ will
intersect the three concentric and homologous bounding, critical and dilatancy
surfaces at the socalled `image' backstress ratio tensor $\alpha_{ij}^{b}$,
$\alpha_{ij}^{c}$, and $\alpha_{ij}^{d}$ respectively as illustrated in
Figure~\ref{DM_Pi_Plane}. The corresponding values of the aforementioned tensors
are
%
\begin{subequations}
\begin{align}
\alpha_{ij}^{b} & =\sqrt{\frac{2}{3}}[M\exp{(n^{b}\psi)}m]n_{ij}=\left(
\sqrt{\frac{2}{3}}\alpha_{\theta}^{b}\right) n_{ij}\label{DM04Alpha_b}\\
\alpha_{ij}^{c} & =\sqrt{\frac{2}{3}}[Mm]n_{ij}=\left( \sqrt{\frac{2}{3}%
}\alpha_{\theta}^{c}\right) n_{ij}\label{DM04Alpha_c}\\
\alpha_{ij}^{d} & =\sqrt{\frac{2}{3}}[M\exp{(n^{d}\psi)}m]n_{ij}=\left(
\sqrt{\frac{2}{3}}\alpha_{\theta}^{d}\right) n_{ij}\label{DM04Alpha_d}%
\end{align}
%
where $\psi=ee_{c}$ is the state parameter and $n^{b}$ and $n^{d}$ are material constants.%
\begin{figure}
[ptb]
\begin{center}
\includegraphics[height=5cm]{piplane.eps}%
\caption{ Schematic illustration of the yield, critical, dilatancy, and
bounding surfaces in the $\pi$plane of deviatoric stress ratio space (after
\cite{Dafalias04}). }%
\label{DM_Pi_Plane}%
\end{center}
\end{figure}
%EndExpansion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Plastic Flow}
In general the model employs a nonassociative flow rule.
The plastic strain
is given by
%
\end{subequations}
\begin{equation}
\Delta{\epsilon}_{ij}^{p}=\Delta{\lambda}R_{ij}=\Delta{\lambda}(R_{ij}^{\prime
}\frac{1}{3}D\delta_{ij})\label{DM04PF1}%
\end{equation}
%
where $\Delta{\lambda}$ is the nonnegative loading index. The deviatoric part
$R_{ij}^{\prime}$ of the flow direction $R_{ij}$ is defined as the deviatoric
part of the normal to the critical surface at the image point $\alpha_{ij}%
^{c}$.
%
\begin{subequations}
\begin{align}
R_{ij}^{\prime} & =Bn_{ij}+C(n_{ik}n_{kj}\frac{1}{3}\delta_{ij}%
)\label{DM04PF2}\\
B & =1+\frac{3}{2}\frac{1c}{c}g\cos{3\theta},\;\;C=3\sqrt{\frac{3}{2}}%
\frac{1c}{c}g\label{DM04PF3}%
\end{align}
The dilatancy coefficient $D$ in Equation~(\ref{DM04PF1}) is defined by
%
\end{subequations}
\begin{equation}
D=A_{d}(\alpha_{ij}^{d}\alpha_{ij})n_{ij}=A_{d}\left( \sqrt{\frac{2}{3}%
}\alpha_{\theta}^{d}\alpha_{ij}n_{ij}\right) \label{DM04PF4}%
\end{equation}
%
where parameter $A_{d}$ is a function of the state variables. In order to account
for the effect of fabric change during dilatancy, the socalled
fabricdilatancy internal variable $z_{ij}$ has been introduced with an
evolution law which will be presented later. The parameter $A_{d}$ in
Equation~(\ref{DM04PF4}) now can be affected by the aforementioned tensor as
%
\begin{equation}
A_{d}=A_{0}(1+\left\langle z_{ij}n_{ij}\right\rangle )\label{DM04PF5}%
\end{equation}
%
where $A_{0}$ is a material constant. The MacCauley brackets $\left\langle
{}\right\rangle $ operate according to $\left\langle x\right\rangle =x$, if
$x>0$ and $\left\langle x\right\rangle =0$, if $x\leq0$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Evolution Laws}
This model has two tensorial internal variables, namely, the back stressratio
tensor $\alpha_{ij}$ and the fabricdilatancy tensor $z_{ij}$. The evolution
law for the back stressratio tensor $\alpha_{ij}$ is function of the distance
between bounding and current back stress ratio tensor in the form of
%
\begin{equation}
\dot{\alpha}_{ij}=\dot{\lambda}[\frac{2}{3}h(\alpha_{ij}^{b}\alpha
_{ij})]\label{DM04Evo1}%
\end{equation}
%
with $h$ a positive function of state parameters. It is defined as a function
of $e$, $p$, and $\eta$ as follows in order to increase the efficiency of the
model in handling the nonlinear response and reverse loading.
%
\begin{subequations}
\begin{align}
h & =\frac{b_{0}}{(\alpha_{ij}\tilde{{\alpha}}_{ij})n_{ij}}\label{DM04Evo1a}%
\\
b_{0} & =G_{0}h_{0}(1c_{h}e)\left( \frac{p}{p_{at}}\right) ^{1/2}%
\label{DM04Evo1b}%
\end{align}
%
where $\tilde{{\alpha}}_{ij}$ is the initial value of $\alpha_{ij}$ at
initiation of a new loading process and is updated to the new value when the
denominator of Equation~(\ref{DM04Evo1a}) becomes negative. The $h_{0}$ and
$c_{h}$ are material constants. Finally the evolution law for the
fabricdilatancy tensor $z_{ij}$ is introduced as
%
\end{subequations}
\begin{equation}
\dot{z}_{ij}=c_{z}\left\langle \dot{\lambda}D\right\rangle (z_{max}%
n_{ij}+z_{ij})\label{DM04Evo2}%
\end{equation}
%
with $c_{z}$ and $z_{max}$ as the material constants that control the maximum
value of $z_{ij}$ and its pace of evolution. Equations~(\ref{DM04PF5}) and
(\ref{DM04Evo2}) are introduced to enhance the subsequent contraction in
unloading after the dilation in loading to generate the wellknown butterfly
shape in the undrained cyclic stress path.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Verification and Validation Examples}
\label{Examples}
Prediction of mechanical behavior comprises use of computational model to foretell the
state of a physical system under consideration under conditions for which the
computational model has not been validated
\citep{Oberkampf2002}. Confidence in predictions relies heavily on proper
Verification and Validation (V\&V) process.
\paragraph{\bf Verification} is the process of determining that a
model implementation accurately represents the developer's conceptual
description and specification. It is a Mathematics issue. Verification provides
evidence that the model is solved correctly. Verification is also meant to
identify and remove errors in computer coding and verify numerical algorithms
and is desirable in quantifying numerical errors in computed solution.
\paragraph{\bf Validation} is the process of determining the degree to which a
model is accurate representation of the real world from the perspective of the
intended uses of the model. It is a Physics issue. {Validation provides evidence
that the correct model is solved.} Validation serves two goals, namely, (a)
tactical goal in identification and minimization of uncertainties and errors
in the computational model, and (b) strategic goal in increasing confidence
in the quantitative predictive capability of the computational model.
Verification and Validation (V\&V) procedures are the primary means of assessing
accuracy in modeling and computational simulations. The V\&V procedures are
essential tools in building confidence and credibility in modeling and
computational simulations.
Figure (\ref{VerifValidFund01}) depicts relationships between verification,
validation and the computational solution.
%
\begin{figure}[!ht]
\begin{center}
{\includegraphics[width=12cm]{VerifValidFund01.eps}}
\label{VerifValidFund01}
\caption{Schematic description of Verification and Validation procedures and
the computational solution. (adopted from \cite{Oberkampf2002}).}
\end{center}
%\vspace*{1cm}
\end{figure}
%
%
In order to verify the upU formulation, a
number of closed form or very accurate solutions were
used.
%\citep{Coussy95}.
Presented bellow are select examples used for verification. In particular,
vertical (1D) consolidation, line injection of fluid in a reservoir
and shock wave propagation in porous medium examples are used to verify the formulation and implementation.
It should be noted that comparison of numerical and closed
form
solutions (verification) is satisfactory within the limitations of numerical
accuracy whereby errors are introduced by finite element discretization, finite
time step size and other artifacts of finite precision arithmetic.
%
In addition to formulation verification, numerical implementation was verified
using a comprehensive set of software tools available at the Computational
Geomechanics Lab at UCD.
%
Validation was performed on set of physical tests on Toyoura sand and is also
presented below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Verification: Vertical Consolidations}
\label{Vertical_Consolidations_upU}
In this section, a vertical consolidation solution is used to verify the upU
finite element formulation.
For the case of linear isotropic elastic soil, \cite{Terzaghi43} developed an
analytic solution (based on infinite trigonometric series). This solution is
used to verify low frequency solidfluid coupling problem.
The consolidation model is assumed to have a depth of $H_0$, Young's modulus $E$,
Poisson's ratio $\nu$, permeability $k$, and the vertical load at the surface is
$p_0$. Finite element model is represented with 10 brick elements with
appropriate boundary conditions to mimic 1D consolidation.
It is important to note that the analytic solutions for the vertical
consolidation is based on assumption that both the soil particles and the
pore fluid (water) are completely incompressible.
Developed $upU$ finite element model, can simulate realistic compressibility
of soil particles and pore fluid. However, in this case, for the purpose of
verification the bulk modulus of soil particles and pore fluid was input as a
large number to mimic their respective incompressibility.
The parameters used for this model are shown in
Table~(\ref{ElasticConsolidation}).
%
%
\begin{table}[!htbp]
\caption{ \label{ElasticConsolidation}
Parameters for the 1D consolidation soil model.}
\begin{center}
\begin{tabular}{ccc}
\hline
parameter & symbol & value \\
\hline
gravity acceleration & $g$ & 9.8~$m/s^2$ \\
soil matrix Young's Modulus & $E$ & $2.0\times 10^7$~$kN/m^2$ \\
soil matrix Poisson's ratio & $v$ & 0.2 \\
soil density & $\rho_s$ & $2.0\times 10^3$~$kg/m^3$ \\
water density & $\rho_f$ & $1.0\times 10^3$~$kg/m^3$ \\
permeability & $k$ & $1.0\times 10^{7}$~$m/s$ \\
solid particle bulk modulus & $K_s$ & $1.0\times 10^{20}$~$kN/m^2$ \\
fluid bulk modulus & $K_f$ & $2.2\times 10^9$~$kN/m^2$ \\
porosity & $n$ & 0.4 \\
% Biot coefficient & $\alpha$ & $1.0$ \\
\hline
\end{tabular}
\end{center}
\end{table}
As an example, a comparison between the numerical and closedform solution (equation (8.171)
from \cite{Coussy95}) for the excess pore pressure isochrones is shown in
Figure~(\ref{OneDComparison}).
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.45\textheight]{ClosedForm_Numerical.eps}
\caption{
\label{OneDComparison}
{ Numerical and closedform excess pore pressure dissipation isochrones.}
}
\end{center}
\end{figure}
%
The agreement is excellent with the note that for
initial stages of loading (when time is close to 0) some differences become
evident. This is to be expected and is an artifact of the time integration of
FEM model.
This situation is
easily resolved using shorter time steps in FEM modeling. However, as this increase in
precision of time integration is needed only at the very beginning of the
analysis, we opt to use large time steps in order to have efficient long term
solutions (which does not require such short time steps).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Verification: Line Injection of a fluid in a Reservoir}
\label{Line Injection of a fluid in a Reservoir}
The analytical solution for the problem of line injection of fluid into
a reservoir (porous medium) is also available from \cite{Coussy95}.
The problem is comprised of a reservoir of infinite extent composed of an
isotropic, homogeneous
and saturated poroelastic material. The cylindrical well of negligible
dimensions, is used to inject a fluid in all directions
orthogonal to the well axis ($z$). As a result of
the axisymmetry, all quantities depends on distance $r$ from the well only.
The injection is assumed instantaneous at time $t=\Gamma$.
The flow rate of fluid mass injection is constant and equal to $q$.
The finite element model used was made up of three dimensional brick elements.
The axisymmetry was mimicked by creating $1/4$ of the completely circular
(axisymmetric) mesh. The use of 3D brick elements was done in order to verify
their formulation and implementation,
while one quarter of the mesh was used in order to facilitate application of
boundary conditions, and since the problems is axisymmetric, to save on
simulation requirements. In theory, any finite slice of the axisymmetric space
cold have been chosen for the model, the alignment of radial boundaries with
coordinate axes proves very useful in defining boundary conditions.
%
Figure~(\ref{fig32}) shows the finite element model.
%
\begin{figure}[!hbpt]
\begin{center}
\includegraphics[width=10cm]{Inject1.eps}
\end{center}
\caption{\label{fig32} {The mesh used for the study of line injection
problem.}}
\end{figure}
Table~(\ref{table8}) presents material parameters and constants used in the
analysis.
%
\begin{table}
\begin{center}
\caption{\label{table8}{Material Properties used to study the line injection problem}}
\begin{tabular}{ccc} \hline
parameter & symbol & value \\ \hline
Poisson ratio & $\nu$ & 0.2 \\
Young's modulus & E & $1.2 \times 10^6~kN/m^2$ \\
Solid particle bulk modulus & $K_s$ & $3.6 \times 10^7~kN/m^2$ \\
Fluid bulk modulus & $K_f$ & $1.0 \times 10^{17}~kN/m^2$ \\
%Undrained bulk modulus & $K_u$ & $6.0 \times 10^7~kN/m^2$ \\
%Skeleton bulk modulus & $K$ & $6.7 \times 10^5~$kN/m^2 \\
Solid density & $\rho_s$ & $2700~kg/m^3$ \\
Fluid density & $\rho_f$ & $1000~kg/m^3$ \\
%Fluid diffusivity coefficient & $c_f$ & $0.4973~m^2/s$ \\
Porosity & n & $0.4$ \\
Permeability & k & $3.6 \times 10^{6}~m/s$ \\ \hline
\end{tabular}
\end{center}
\end{table}
The results are recorded at three points at the radii of $10$~cm, $50$~cm and
$85$~cm. Comparison of analytic and numerical results are shown in
Fig.~(\ref{fig33}).
%
\begin{figure}[!hbpt]
\begin{center}
\includegraphics[width=0.49\textwidth]{InjectionPressure.new.ps.eps}
\includegraphics[width=0.475\textwidth]{InjectionDisplacement.new.ps.eps}
\end{center}
\caption{\label{fig33} The comparison between analytical and
computational solution for pore pressures (left) and radial displacements (right).}
\end{figure}
%
The comparison starts at 1 second after the injection.
Comparison shows very good match between analytic and numerical solution. The
(expected) mismatch is observed in the region close to the point of injection of
fluid. In our case this point represents a singularity and for the finite
element mesh, a very small opening was left in the mesh to mimic this (almost)
singular point.
%
As the distance from the singular injection point increases, the match between
analytic and numerical solution becomes excellent. The match is also improving
in time, that is, at 2 seconds, the match for both the pore pressures and the
displacements is almost perfect, and it continues to improve as the time goes
by.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Verification: Shock Wave Propagation in Saturated Porous Medium}
\label{ShockWave}
In order to verify the dynamic behavior of the system,
an analytic solution developed by Gajo
\cite{Gajo1995} and Gajo and Mongiovi \cite{Gajo1995b} for 1D shock wave
propagation in elastic porous medium was used.
A model was
developed consisting of 1000 eight node brick elements, with boundary conditions
that mimic 1D behavior. In particular, no displacement of solid ($u_x = 0 , u_y
= 0$) and fluid ($U_x = 0 , U_y = 0$) in $x$ and $y$ directions is allowed along
the height of the model. Bottom nodes have full fixity for solid ($u_i = 0$) and
fluid ($U_i=0$) displacements while all the nodes above base are free to move in
$z$ direction for both solid and fluid. Pore fluid pressures are free to develop
along the model. Loads to the model
consist of a unit step function (Heaviside) applied as (compressive) displacements
to both solid and fluid phases of the model, with an amplitude of $0.001$~cm.
The upU model dynamic system of equations was integrated using Newmark algorithm
(see section~\ref{Time_Integration}). Table~\ref{table_shock_wave} gives
relevant parameters for this verification.
%
%
\begin{table}
\begin{center}
\caption{\label{table_shock_wave}{Simulation parameters used for the shock wave
propagation verification problem.}}
\begin{tabular}{ccc} \hline
Parameter & Symbol & Value \\ \hline
Poisson ratio & $\nu$ & 0.3 \\
Young's modulus & E & $1.2 \times 10^6~kN/m^2$ \\
Solid particle bulk modulus & $K_s$ & $3.6 \times 10^7~kN/m^2$ \\
Fluid bulk modulus & $K_f$ & $2.17 \times 10^{6}~kN/m^2$ \\
Solid density & $\rho_s$ & $2700~kg/m^3$ \\
Fluid density & $\rho_f$ & $1000~kg/m^3$ \\
Porosity & n & $0.4$ \\
Newmark parameter & $\gamma$ & 0.6 \\ \hline
%
\end{tabular}
\end{center}
\end{table}
%
Two set of permeability of material were used in our verification. The first
model had permeability set $k=10^{6}$~cm/s which creates very high coupling
between porous solid and pore fluid. The second model had permeability set to
$k=10^{2}$~cm/s which, on the other hand creates a low coupling between porous
solid and pore fluid.
Comparison of simulations and the analytical solution are presented in
Figure~\ref{DynamicVerification}.
%
\begin{figure}[!hbpt]
\begin{center}
\includegraphics[width=0.99\textwidth]{NewMark_gamma0.6_K102and106cms.eps}
\end{center}
\caption{\label{DynamicVerification} Compressional wave in both solid and fluid,
comparison with closed form solution.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Validation: Material Model}
Validation is performed by comparing experimental (physical) results and
numerical (constitutive) simulations for the Toyoura sand.
It should be noted that we have not done validation against 2D or 3D tests (say
centrifuge tests) as characterization of sand used in centrifuge experiments is usually
less than complete for use with advanced constitutive models. Moreover, as our
approach seeks to make predictions of prototype behavior, scaling down
models (and using them for comparison with numerical predictions) brings forward
issues of physics of scaling which we would rather stay out of.
%
The material parameters
used are from \cite{Dafalias04} and
are listed in Table ~(\ref{DMTest_parameters}). Several simulated results are
compared with the experimental data published by \cite{Verdugo1996}.
\begin{table}[tbh]
\caption{ Material parameters of DafaliasManzari model. }
\begin{center}
% {\normalsize
\begin{tabular}
[c]{lcclcc}\hline
\multicolumn{2}{l}{material parameter} & value &
\multicolumn{2}{l}{material parameter} & value\\\hline
Elasticity & $G_{0}$ & 125 kPa & Plastic modulus & $h_{0}$ & 7.05\\
& $v$ & 0.05 & & $c_{h}$ & 0.968\\
Critical sate & $M$ & 1.25 & & $n_{b}$ & 1.1\\
& $c$ & 0.712 & Dilatancy & $A_{0}$ & 0.704\\
& $\lambda_{c}$ & 0.019 & & $n_{d}$ & 3.5\\
& $\xi$ & 0.7 & Fabricdilatancy & $z_{max}$ & 4.0\\
& $e_{r}$ & 0.934 & & $c_{z}$ & 600.0\\
Yield surface & $m$ & 0.01 & & & \\\hline
\end{tabular}
% }
\end{center}
\label{DMTest_parameters}%
\end{table}
Figure~(\ref{drainedTest}) presents both loading and unloading triaxial drained
test simulation results for a relatively low triaxial isotropic pressure of
$100$~kPa but with different void ratios of $e_0 = 0.831, 0.917, 0.996$ at the end of
isotropic compression.
%
\begin{figure}[!htbp]
\begin{center}
\begin{minipage}{0.45\textwidth}
\includegraphics[width=\textwidth]{p100sq.eps} \\
\vspace{0.02\textheight}
\includegraphics[width=\textwidth]{p100eq.eps} \\
\vspace{0.01\textheight}
\includegraphics[width=\textwidth]{p500sq.eps} \\
\vspace{0.01\textheight}
\includegraphics[width=\textwidth]{p500eq.eps}
\vspace{0.05\textheight}
\end{minipage}
\hspace{0.05\textwidth}
\begin{minipage}{0.45\textwidth}
\includegraphics[width=0.95\textwidth]{P100.eps} \\
\vspace{0.08\textheight}
\includegraphics[width=0.95\textwidth]{P500.eps}
\end{minipage}
\vspace*{1.5cm}
\caption{
\label{drainedTest}
Left: Experimental data;
Right: Simulated results.
}
\end{center}
\end{figure}
%
During the loading stage, one can observes the hardening and then softening
together with the contraction and then dilation for the denser sand, while only
hardening together with contraction for the looser sand.
The significance of
the state parameter to the soil modeling is clear from the very different
responses with different void ratios at the same triaxial isotropic pressure.
Figure~(\ref{drainedTest}) also shows both loading and unloading triaxial drained
test simulation results for a relatively high triaxial isotropic pressure of
$500$~kPa but with different void ratio of $e_0$=0.810, 0.886, 0.960 at the end of
isotropic compression. Similar phenomenon are observed as with tests
(physical and numerical) for relatively low triaxial isotropic pressure.
However, due to the higher confinement
pressure, the stressstrain responses are higher at the same strain, which
proves the significant pressure dependent feature for the sands.
Figure~(\ref{undrainedTest}) presents both loading and unloading triaxial
undrained test simulation results for a dense sand with the void ratio of $e_0$=0.735
at the end of isotropic compression but with different isotropic compression
pressures of $p_0$ = 100, 1000, 2000, 3000 kPa.
%
\begin{figure}[!htbp]
\begin{center}
\begin{minipage}{0.45\textwidth}
\includegraphics[width=\textwidth]{e735sq.eps} \\
\vspace{0.02\textheight}
\includegraphics[width=\textwidth]{e735pq.eps} \\
\vspace{0.01\textheight}
\includegraphics[width=\textwidth]{e833sq.eps} \\
\vspace{0.01\textheight}
\includegraphics[width=\textwidth]{e833pq.eps}
\vspace{0.05\textheight}
\end{minipage}
\hspace{0.05\textwidth}
\begin{minipage}{0.45\textwidth}
\includegraphics[width=0.95\textwidth]{e735.eps} \\
\vspace{0.08\textheight}
\includegraphics[width=0.95\textwidth]{e833.eps}
\end{minipage}
\vspace*{1.5cm}
\caption{
\label{undrainedTest}
Left: Experimental data;
Right: Simulated results.
}
\end{center}
\end{figure}
%
During the loading stage, one
observes that each of responses are close to the critical state line for the very
various range of isotropic compression pressures.
%
For the higher isotropic compression pressure,
the contraction response with softening is clearly observed, while for the
smaller isotropic compression pressure,
it is a dilation response without softening.
%%
Close matching of physical testing data with constitutive predictions represents
a satisfactory validation of our material model. This validation with previous
verification gives us confidence that predictions (presented in next section)
represent well the real, prototype behavior.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Liquefaction of Level and Sloping Grounds}
\label{upU_SD_Examples}
Liquefaction of level and sloping grounds represents a very common behavior
during earthquakes. Of interest is to estimate settlement for level ground, and
horizontal movements for sloping grounds. In next few sections, presented are
results for a 1D, vertical (level ground) and sloping ground cases for dense and
loose sand behavior during seismic shaking.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Model Description}
\label{One_Soil_Column_Subjected_to_Earthquake_Ground_Motion}
Vertical soil column consists of a multipleelements subjected to an earthquake shaking.
The soil is assumed to be Toyoura sand and the calibrated parameters are from
\cite{Dafalias04}, and were given in the
Table~(\ref{DMTest_parameters}). The other parameters, related to the boundary
value problem are given in table~(\ref{OneC_Material_Parameter}).
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}
\begin{center}
\caption{ \label{OneC_Material_Parameter}
Additional parameters used in boundary value problem simulations (other than
material parameters from the Table~(\ref{DMTest_parameters})).}
\begin{tabular}{ccc} \hline
Parameter & Symbol & Value \\ \hline
Solid density & $\rho_s$ & $2700~kg/m^3$ \\
Fluid density & $\rho_f$ & $1000~kg/m^3$ \\
Solid particle bulk modulus & $K_s$ & $3.6 \times 10^{7}~kN/m^2$ \\
Fluid bulk modulus & $K_f$ & $2.2 \times 10^{6}~kN/m^2$ \\
%Biot coefficient & $\alpha$ & $1.0$ \\
permeability & $k$ & $5.0\times 10^{4}$~$m/s$ \\
HHT parameter & $\alpha$ & 0.2 \\ \hline
%
\end{tabular}
\end{center}
\end{table}
%
For tracking convenience, the mesh elements are labeled from E01 (bottom) to E10
(surface) and nodes at each layers are labeled from A (bottom) to K (surface).
A static application of gravity analysis is performed before seismic
excitation.
%
The resulting fluid hydrostatic pressures and soil stress states
along the soil column serve as initial conditions for the subsequent dynamic
analysis.
It should be noted that the self weight loading is performed on an initially
zero stress (unloaded) soil column and that the material model and
numerical integration algorithms are powerful
enough to follow through this early loading with proper evolution.
%
%
The boundary conditions are such that the soil and water displacement degree of
freedom (DOF)
at the bottom surface are fixed, while the pore pressure DOFs are free;
the soil and water displacement DOFs
at the upper surface are free upwards to simulate the upward drainage. The pore
pressure DOFs are fixed at surface thus setting pore pressure to zero. On the sides,
soil skeleton and water are prevented from moving in horizontal directions while
vertical movement of both is free.
It is emphasized that those displacements (of soil skeleton and pore fluid)
are different.
In order to simulate the 1D behavior,
all DOFs at the same depth level are connected in a masterslave fashion.
Modeling of sloping ground is done by creating a constant horizontal load, sine of
inclination angle, multiplied by the self weight of soil column, to mimic
sloping ground. In addition to that, for a sloping ground, there
should be a constant flow (slow) downhill, however this is neglected in our
modeling.
%
The permeability is assumed to be isotropic $k=5.0{\times}10^{4}$~m/s.
The input acceleration time history (Figure~(\ref{OneCoulmnGM2}))
is taken from the recorded horizontal
acceleration of Model No.1 of VELACS project \cite{Arulanandan93} by Rensselaer
Polytechnic Institute (http://geoinfo.usc.edu/gees/velacs/).
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=16cm]{SinglePileGM.ps.eps}
\caption{
\label{OneCoulmnGM2}
{Input earthquake ground motion for the soil column.}
}
\end{center}
\end{figure}
%
The magnitude of the motion is close to $0.2$~g, while main shaking lasts
for about 12 seconds (from $1$~s to $13$~s).
For the sloping ground model a slope of \%~3 was considered.
It should be emphasized that the soil parameters are related to Toyoura sand,
not Nevada sand which is used in VELACS project. The purpose of presented simulation
is to show the predictive performance using verified and validated formulation,
algorithms, implementation and models.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Behavior of Saturated Level Ground}
\label{LevelLiq}
Figure~(\ref{Level_Loose}) describes the response of the sample with
loose sand
$e_{0}=0.85$. This figure shows the typical mechanism of cyclic decrease in effective
vertical stress due to pore pressure build up as expected for the looser than
critical granular material.
The lower layers show only the reduction of
effective vertical stress from the beginning. Once the effective vertical (and
therefore confining) stress approaches the smaller values, signs of the
socalled butterfly shape can be observes in the stress path. Similar
observation can made in the upper layers which have the smaller confining pressure
comparing to the lower layers from the very beginning due to lower surcharge.
The upper
layers have lower confining pressure (lower surcharge) at the beginning of the
shaking, hence less contractive response is expected in these layer; however,
soon after the initiation of the shaking these top layers start showing the
liquefaction state and that type of response continues even after the end of the shaking.
%
The top section of the model has remained liquefied well past the end of
shaking. This is explained by the large supply of pore fluid from lower layers,
for which the dissipation starts earlier. For example, for the lowest layer, the
observable drop in excess pore pressure start as soon as the shaking ends,
while, the upper layers then receive this dissipated pore fluid from lower
layers and do liquefy (or continue being liquefied) well past end of shaking
(which happens at approximately $13$~seconds).
%
It is very important to note the significance of this incoming pore water flux on the
pore water pressure of the top layers. Despite the less contractive response of
soil skeleton at the top elements, the transient pore water flux, that enters
these elements from the bottom, forces those to a liquefaction state. In other
words, the top elements have not liquefied only due to their loose state but also
because of the water flow coming from the bottom layers.
%
The maximum horizontal
strains can be observed in the middle layers due to liquefaction and prevents
upper layers from experiencing larger strains. The displacements of water and
soil are presented in the last column. It shows that in all layers the upward
displacement of water is larger than the downward displacement of soil. This
behavior reflects soil densification during shaking.
Figure~(\ref{Level_Dense}) describes the response of the sample with dense
sand ($e_{0}=0.75$). This figure also shows the typical mechanism of cyclic
decrease in effective vertical stress. However, in case of this dense sample the
decreasing rate of the effective confining pressure is much smaller than what
was observed in the loose
sample. Signs of the partial butterfly shape in the effective stress path can be
observed from early stages of shaking. The butterfly is more evident in the
upper layers with the lower confining pressure, i.e. more dilative response. In
later stages of the shaking, i.e. when the confining pressure reduces to smaller
value the butterfly shape of the stress path gets more pronounced due to
having more dilative response in the lower confining pressure based on
CSSM concept. In comparing this dense case to the case of shaking the
loose sand column, the
current case does not show any major sign of liquefaction (when stress
ratio $r_{u}=1$). This is due to the less contractive (more dilative) response of the
sand in this case, which is coming from the the denser state of the sample.
Because of having partial segments of dilative response, the whole column of the
sand has not loosed its strength to the extent that happened for the case of
loose sand and therefore smaller values of horizontal strains has been observed
in the results. The absolute values of soil and water vertical displacements are
also smaller than the case of loose sand which can be again referred to the less
overall contractive response in this case.
Overall, it can be noted that the response in the case of loose sand
($e_{0}=0.85$) is mainly
below the dilatancy surface (phase transformation surface) while the denser
sand sample with $e_{0}=0.75$ shows partially dilative response referring to the
denser than critical state.
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=11cm]{e85L.eps}
\includegraphics[width=1.6cm,angle=90]{OneCShaking.eps}
\caption{ \label{Level_Loose} {Seismic results for (loose sand) soil
column in level ground ($e_0 = 0.85$).} }
\end{center}
\end{figure}
%
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=11cm]{e75L.eps}
\includegraphics[width=1.6cm,angle=90]{OneCShaking.eps}
\caption{ \label{Level_Dense} {Seismic results for (dense sand) soil
column in level ground ($e_0 = 0.75$).} }
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Behavior of Saturated Sloping Ground}
\label{SlopeLiq}
Figures~(\ref{Sloping_Loose}) and (\ref{Sloping_Dense}) present the result of the
numerical simulations for shaking the inclined soil columns (toward right) with
loose and dense sand samples, respectively. The inclination of the soil column
results in presence
of the offset shear stress to the right side. This essentially poses asymmetric
horizontal shear stresses (toward the direction of inclination) during cycles of
shaking. On one hand, this offset shear stress makes the sample more dilative in
the parts of shaking toward the right side (think about the state distance from
the phase transformation line or dilatancy line in the $pq$ space). As a result
asymmetric butterfly loops will be induced causing the soil to regain its
stiffness and strength ($p$) in the dilative parts of the corresponding cycles,
therefore only instantaneous spikes of $r_{u}=1$ can be observed in case of
the sloped columns of soil. There is also a permanent liquefaction in terms of
having stationary portions of $r_{u}=1$ in this case. On the other hand, the
offset shear stress results generation of more horizontal strains in the
portions of loading which are directed toward the right side than those which
are directed back toward the left side. As a result the horizontal shear strains
will accumulative toward the right side and create larger permanent horizontal
displacement comparing to the case of level ground soil column. Since the
overall dilative response of the dense sample, i.e.
Figure~(\ref{Sloping_Dense}),
is larger than that of the loose sample, i.e. Figure~(\ref{Sloping_Loose}), the
dense sample shows stiffer response and therefore less accumulative horizontal
shear strains than the loose sample. The difference in predicted horizontal
displacements is almost three times, that is, for the dense sample the final,
maximum horizontal displacement is approx. $0.5$~m, while for the loose sand
sample, it almost $1.5$~m,
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=11cm]{e85s.eps}
\includegraphics[width=1.6cm,angle=87]{OneCShaking.eps}
\caption{ \label{Sloping_Loose} {Seismic results for (loose sand)
soil column in sloping ground ($e_0 = 0.85$).} }
\end{center}
\end{figure}
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=11cm]{e75s.eps}
\includegraphics[width=1.6cm,angle=87]{OneCShaking.eps}
\caption{ \label{Sloping_Dense} {Seismic results for (dense sand)
soil column in sloping ground ($e_0 = 0.75$).} }
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Summary}
A numerical modeling approach was presented accounting for fully coupled
(elasticplastic porous solid  fluid) behavior of soils using $upU$ formulation. A critical
state elastoplastic model accounting for the fabric dilatancy
effects was used for all stages of loading, starting from zero stress state,
through self weight and finally to seismic shaking.
The formulation and implementation takes into account velocity proportional
damping (usually called viscous damping) by proper modeling of coupling of
pore fluid and solid skeleton, while
the displacement proportional damping is appropriately modeled using
elastoplasticity and a powerful material model. No additional damping was
used in FEM model.
The formulation and resulting implementation
were verified and validate using a number of unit tests and were then applied to
the
problem of sand liquefaction and cyclic mobility phenomena.
Presented verified formulation and implementation, in
conjunction with a validated material model give us confidence in our ability to
realistically predict behavior of saturated soils during dynamic loading events.
%
Moreover, the developed system, that is closely following the physics
(mechanics) of the problem, allows for investigation of energy dissipation
mechanisms during seismic events.
%
As such it allows for high fidelity simulations that will allow designers to
make a transition from empirically based to rational mechanics based methods.
\begin{thebibliography}{42}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi: #1}\else
\providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi
\bibitem[Argyris and Mlejnek(1991)]{local87}
John Argyris and HansPeter Mlejnek.
\newblock \emph{Dynamics of Structures}.
\newblock North Holland in USA Elsevier, 1991.
\bibitem[Arulanandan and Scott(1993)]{Arulanandan93}
Kandiah Arulanandan and Ronald~F. Scott, editors.
\newblock \emph{Verification of Numerical Procedures for the Analysis of Soil
Liquefaction Problems}.
\newblock A. A. Balkema, 1993.
\bibitem[Chan(1988)]{Chan88}
Andrew {Hin}{Cheong} Chan.
\newblock \emph{A Unified Finite Element Solution to Static and Dynamic
Problems in Geomecahnics}.
\newblock PhD thesis, Department of Civil Engineering, University College of
Swansea, February 1988.
\bibitem[Coussy(1995)]{Coussy95}
Olivier Coussy.
\newblock \emph{Mechanics of Porous Continua}.
\newblock John Wiley and Sons, 1995.
\newblock ISBN 471 95267 2.
\bibitem[Dafalias and Manzari(2004)]{Dafalias04}
Yannis~F. Dafalias and Majid~T. Manzari.
\newblock Simple plasticity sand model accounting for fabric change effects.
\newblock \emph{Journal of Engineering Mechanics}, 130\penalty0 (6):\penalty0
622634, 2004.
\bibitem[Davis and Duff(1997)]{Davis97b}
Timothy~A. Davis and Iain~S. Duff.
\newblock A combined unifrontal/multifrontal method for unsymmetric sparse
matrices.
\newblock Technical Report TR97016 (UofF) and TR97046 (RAL), University of
Florida and Rutherford Appleton Laboratory, 1997.
\bibitem[Dennis and Schnabel(1983)]{Schnabel83}
J.~E. Dennis, Jr. and Robert~B. Schnabel.
\newblock \emph{Numerical Methods for Unconstrained Optimization and Nonlinear
Equations}.
\newblock Prentice Hall , Engelwood Cliffs, New Jersey 07632., 1983.
\bibitem[Elgamal et~al.(2002)Elgamal, Yang, and Parra]{Elgamal2002}
Ahmed Elgamal, Zhaohui Yang, and Ender Parra.
\newblock Computational modeling of cyclic mobility and postliquefaction site
response.
\newblock \emph{Soil Dynamics and Earthquake Engineering}, 22:\penalty0
259271, 2002.
\bibitem[Elgamal et~al.(2003)Elgamal, Yang, Parra, and Ragheb]{Elgamal2003}
Ahmed Elgamal, Zhaohui Yang, Ender Parra, and Ahmed Ragheb.
\newblock Modeling of cyclic mobility in saturated cohesionless soils.
\newblock \emph{International Journal of Plasticity}, 19\penalty0 (6):\penalty0
883905, 2003.
\bibitem[Gajo(1995)]{Gajo1995}
A.~Gajo.
\newblock Influence of viscous coupling in propagation of elastic waves in
saturated soil.
\newblock \emph{{ASCE} Journal of Geotechnical Engineering}, 121\penalty0
(9):\penalty0 636644, September 1995.
\bibitem[Gajo and Mongiovi(1995)]{Gajo1995b}
A.~Gajo and L.~Mongiovi.
\newblock An analytical solution for the transient response of saturated linear
elastic porous media.
\newblock \emph{International Journal for Numerical and Analytical Methods in
Geomechanics}, 19\penalty0 (6):\penalty0 399413, 1995.
\bibitem[Gajo et~al.(1994)Gajo, Saetta, and Vitaliani]{Gajo1994}
A.~Gajo, A.~Saetta, and R.~Vitaliani.
\newblock Evaluation of three and twofield finite element methods for the
dynamic response of saturated soil.
\newblock \emph{International Journal for Numerical Methods in Engineering},
37:\penalty0 12311247, 1994.
\bibitem[Hilber et~al.(1977)Hilber, Hughes, and Taylor]{Hilber77}
H.~M. Hilber, T.~J.~R. Hughes, and R.~L. Taylor.
\newblock Improved numerical dissipation for time integration algorithms in
structural dynamics.
\newblock \emph{Earthquake Engineering and Structural Dynamics}, 5:\penalty0
283292, 1977.
\bibitem[Hughes and Liu(1978{\natexlab{a}})]{Hughes78a}
T.~J.~R. Hughes and W.~K. Liu.
\newblock Implicitexplicit finite elements in transient analysis: Stability
theory.
\newblock \emph{{ASME} Journal of Applied Mechanics}, 45:\penalty0 371374,
June 1978{\natexlab{a}}.
\bibitem[Hughes and Liu(1978{\natexlab{b}})]{Hughes78b}
T.~J.~R. Hughes and W.~K. Liu.
\newblock Implicitexplicit finite elements in transient analysis:
Implementation and numerical examples.
\newblock \emph{{ASME} Journal of Applied Mechanics}, 45:\penalty0 375378,
June 1978{\natexlab{b}}.
\bibitem[Hughes(1987)]{local86}
{Thomas J. R.} Hughes.
\newblock \emph{The Finite Element Method ; Linear Static and Dynamic Finite
Element Analysis}.
\newblock Prentice Hall Inc., 1987.
\bibitem[Jeremi{\' c}(2001)]{Jeremic2001}
Boris Jeremi{\' c}.
\newblock Line search techniques in elasticplastic finite element
computations in geomechanics.
\newblock \emph{Communications in Numerical Methods in Engineering},
17\penalty0 (2):\penalty0 115125, January 2001.
\bibitem[Jeremi{\' c} and Sture(1997)]{Jeremic96b}
Boris Jeremi{\' c} and Stein Sture.
\newblock Implicit integrations in elastoplastic geotechnics.
\newblock \emph{International Journal of Mechanics of CohesiveFrictional
Materials}, 2:\penalty0 165183, 1997.
\bibitem[Jeremi{\' c} and Sture(1998)]{Jeremic97d}
Boris Jeremi{\' c} and Stein Sture.
\newblock Tensor data objects in finite element programming.
\newblock \emph{International Journal for Numerical Methods in Engineering},
41:\penalty0 113126, 1998.
\bibitem[Jeremi{\' c} and Yang(2002)]{Jeremic2000f}
Boris Jeremi{\' c} and Zhaohui Yang.
\newblock Template elasticplastic computations in geomechanics.
\newblock \emph{International Journal for Numerical and Analytical Methods in
Geomechanics}, 26\penalty0 (14):\penalty0 14071427, 2002.
\bibitem[Li and Wang(1998)]{Li98}
X.~S. Li and Y.~Wang.
\newblock Linear representation of steadystate line for sand.
\newblock \emph{Journal of Geotechnical and Geoenvironmental Engineering},
124\penalty0 (12):\penalty0 12151217, 1998.
\bibitem[Manzari and Dafalias(1997)]{Manzari97}
M.~T. Manzari and Y.~F. Dafalias.
\newblock A critical state twosurface plasticity model for sands.
\newblock \emph{G{\' e}otechnique}, 47\penalty0 (2):\penalty0 255272, 1997.
\bibitem[Mc{K}enna(1997)]{McKenna97}
Francis~Thomas Mc{K}enna.
\newblock \emph{Object Oriented Finite Element Programming: Framework for
Analysis, Algorithms and Parallel Computing}.
\newblock PhD thesis, University of California, Berkeley, 1997.
\bibitem[Newmark(1959)]{Newmark1959}
N.~M. Newmark.
\newblock A method of computation for structural dynamics.
\newblock \emph{{ASCE} Journal pf the Engineering Mechanics Division},
85\penalty0 (6794), 1959.
\bibitem[Oberkampf et~al.(2002)Oberkampf, Trucano, and Hirsch]{Oberkampf2002}
William~L. Oberkampf, Timothy~G. Trucano, and Charles Hirsch.
\newblock Verification, validation and predictive capability in computational
engineering and physics.
\newblock In \emph{Proceedings of the Foundations for Verification and
Validation on the 21st Century Workshop}, pages 174, Laurel, Maryland,
October 2223 2002. Johns Hopkins University / Applied Physics Laboratory.
\bibitem[Pastor et~al.(1990)Pastor, Zienkiewicz, and Chan]{Pastor90}
M.~Pastor, O.~C. Zienkiewicz, and A.~H.~C. Chan.
\newblock Generalized plasticity and the modeling of soil behaviour.
\newblock \emph{International Journal for Numerical and Analytical Methods in
Geomechanics}, 14:\penalty0 151190, 1990.
\bibitem[Prevost(1985)]{Prevost85}
Jean~H. Prevost.
\newblock Wave propagation in fluidsaturated porous media: an efficient finite
element procedure.
\newblock \emph{Soil Dynamics and Earthquake Engineering}, 4\penalty0
(4):\penalty0 183 202, 1985.
\bibitem[Richart et~al.(1970)Richart, Hall, and Woods]{Richart70}
F.~E. Richart, J.~R. Hall, and R.~D. Woods.
\newblock \emph{Vibration of Soils and Foundations, International Series in
Theoretical and Applied Mechanics}.
\newblock PrenticeHall, Englewood Cliffs, N.J., 1970.
\bibitem[Roache(1998)]{Roach1998}
Patrick~J. Roache.
\newblock \emph{Verif\/ication and Validation in Computational Science and
Engineering}.
\newblock Hermosa Publishers, Albuquerque, New Mexico, 1998.
\newblock ISBN 0913478083.
\bibitem[Simon et~al.(1984{\natexlab{a}})Simon, Wu, Zienkiewicz, and
Paul]{Simon1986}
B.~R. Simon, J.~{S.S.} Wu, O.~C. Zienkiewicz, and D.~K. Paul.
\newblock Evaluation of $uw$ and $u\pi$ finite element methods for the
dynamic response of saturated porous media using onedimensional models.
\newblock \emph{International Journal for Numerical and analytical methods in
Geomechanics}, 10:\penalty0 461482, 1984{\natexlab{a}}.
\bibitem[Simon et~al.(1984{\natexlab{b}})Simon, Zienkiewicz, and
Paul]{Simon1984}
B.~R. Simon, O.~C. Zienkiewicz, and D.~K. Paul.
\newblock An analytical solution for the transient response of saturated porous
elastic solids.
\newblock \emph{International Journal for Numerical and analytical methods in
Geomechanics}, 8:\penalty0 381398, 1984{\natexlab{b}}.
\bibitem[S.Sandhu et~al.(1990)S.Sandhu, Shaw, and Hong]{Sandhu1990}
Ranbir S.Sandhu, H.L. Shaw, and S.J. Hong.
\newblock A threefield finite element procedure for analysis of elastic wave
propagation through fluidsaturated soils.
\newblock \emph{Soil Dynamics and Earthquake Engineering}, 9:\penalty0 5865,
1990.
\bibitem[Taiebat and Dafalias(2008)]{Taiebat2007a}
Mahdi Taiebat and Yannis~F. Dafalias.
\newblock {SANISAND}: Simple anisotropic sand plasticity model.
\newblock \emph{International Journal for Numerical and Analytical Methods in
Geomechanics}, 2008.
\newblock (in print, available in earlyview).
\bibitem[Taiebat et~al.(2007)Taiebat, Shahir, and Pak]{Taiebat2007}
Mahdi Taiebat, Hadi Shahir, and Ali Pak.
\newblock Study of pore pressure variation during liquefaction using two
constitutive models for sand.
\newblock \emph{Soil Dynamics and Earthquake Engineering}, 27:\penalty0 6072,
2007.
\bibitem[Terzaghi(1943)]{Terzaghi43}
Karl Terzaghi.
\newblock \emph{Theoretical Soil Mechanics}.
\newblock Wiley, New York, 1943.
\bibitem[Verdugo and Ishihara(1996)]{Verdugo1996}
Ramon Verdugo and Kenji Ishihara.
\newblock The steady state of sandy soils.
\newblock \emph{Soils and Foundations}, 36\penalty0 (2):\penalty0 8191, 1996.
\bibitem[Zienkiewicz and Shiomi(1984)]{Zienkiewicz84}
O.~C. Zienkiewicz and T.~Shiomi.
\newblock Dynamic behaviour of saturated porous media; the generalized {B}iot
formulation and its numerical solution.
\newblock \emph{International Journal for Numerical Methods in Engineering},
8:\penalty0 7196, 1984.
\bibitem[Zienkiewicz and Taylor(2000)]{Ziekienwicz2000a}
O.~C. Zienkiewicz and R.~L. Taylor.
\newblock \emph{The Finite Element Method, Volume 1, The Basis, 5th Edition}.
\newblock Butterworth Heinemann, London, 2000.
\bibitem[Zienkiewicz et~al.(1999)Zienkiewicz, Chan, Pastor, Schrefler, and
Shiomi]{Zienkiewicz99}
O.~C. Zienkiewicz, A.~H.~C. Chan, M.~Pastor, B.~A. Schrefler, and T.~Shiomi.
\newblock \emph{Computational Geomechanics with Special Reference to Earthquake
Engineering}.
\newblock John Wiley and Sons., 1999.
\newblock ISBN 0471982857.
\bibitem[Zienkiewicz et~al.(1980)Zienkiewicz, Chang, and
Bettess]{Zienkiewicz1980}
O.C. Zienkiewicz, C.~T. Chang, and P.~Bettess.
\newblock Drained, undrained, consolidating and dynamic behaviour assumptions
in soils.
\newblock \emph{G{\'e}otechnique}, 30\penalty0 (4):\penalty0 385395, 1980.
\bibitem[Zienkiewicz and Taylor(1991{\natexlab{a}})]{FEMIV1}
Olgierd~Cecil Zienkiewicz and Robert~L. Taylor.
\newblock \emph{The Finite Element Method}, volume~1.
\newblock McGraw  Hill Book Company, fourth edition, 1991{\natexlab{a}}.
\bibitem[Zienkiewicz and Taylor(1991{\natexlab{b}})]{FEMIV2}
Olgierd~Cecil Zienkiewicz and Robert~L. Taylor.
\newblock \emph{The Finite Element Method}, volume~2.
\newblock McGraw  Hill Book Company, {F}ourth edition, 1991{\natexlab{b}}.
\end{thebibliography}
%
\end{document}
\bye