\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}{EEC-9701568}
\cgsn{NSF}{CMS--0337811}

\noreceived{}
\norevised{}
\noaccepted{}


\newpage

\paragraph{Abstract}


Fully   coupled,    porous  solid-fluid  formulation, implementation and related
modeling  and simulation issues are presented in this work. To this end, coupled
dynamic field equations with $u-p-U$ formulation are used to simulate pore fluid
and soil skeleton (elastic--plastic 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 soil--foundation--structure interaction). Fluid compressibility
is also explicitly taken into account, thus allowing excursions into modeling of
limited  cases  of  non-saturated  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 elastic--plastic 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{Elasto--plastic, 
          Finite elements,
          Saturated porous solid--fluid 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 single-phase
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) $u-p$, (b) $u-U$, and
(c) $u-p-U$ formulations. Here, the unknowns are  the soil skeleton displacements
$u$;    the  pore  fluid  (water)    pressure  $p$;   and the pore fluid (water)
displacements $U$. The $u-p$ 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  $u-p$ formulation with the multi-surface
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 $u-p$
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  $u-U$  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  $u-p-U$  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  $u-p-U$ formulation has rarely been implemented into
finite element code. 
%
In   this  paper,  complete $u-p-U$ 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 two-surface 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,    stress-strain)  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  single--phase  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 soil-fluid '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} + {(\alpha-n)}/{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{FEM-IV-1,FEM-IV-2} uses shape
functions to describe 
each  of  the  unknown  fields  (u--p--U) 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 (u--p--U)
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 u--p--U
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{FEM-IV-1}
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 (1-n) \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 (\alpha-n) 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{69-1}  
\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 (\alpha-n) p n_i \ud \Gamma  
                     + 
                     \int_{\Omega} N_K^u (1-n) \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{69-2}  
\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
Dafalias--Manzari     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  Dafalias--Manzari  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 stress--internal variable point
from        the    yield   surface  for  monotonic  loading. The use of Explicit
integration can also result in spurious elastic--plastic
deformations during elastic
unloading    in cyclic loading-unloading   \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 elastic--plastic 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 elastic--plastic 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{local-87}
%%%
\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 non-linear form  \citep{local-87} 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  $u-p-U$  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
Hilber--Hughes--Taylor        (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. Well-known 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 second-order \citep{local-86}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Hilber--Hughes--Taylor Integrator.} 

The   HHT  time integration method uses an alternative residual
form 
%(alternative to Newmark method, see  \cite{local-87})
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}(1-2\alpha), \;\;\; \beta = \frac{1}{4}(1-\alpha)^2
\label{upU_HHT_ab}     
\end{equation}
%%%
the HHT method  is unconditionally stable and second-order accurate \citep{local-86}. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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
(u-p-U)  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 two-surface critical state plasticity model for
sands. Employing the effects of the state parameter on the behavior of sand,
this model presents well-established 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 fabric-dilatancy
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 non-triaxial 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=e-e_{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.97-e)^{2}}{(1+e)}\left(  \dfrac{p}{p_{at}}\right)
^{1/2}p_{at}
\;\;\;\; \mbox{;} \;\;\;\;
K=\dfrac{2(1+\nu)}{3(1-2\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 stress-ratio 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 eight-curve 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)-(1-c)\cos{3\theta}%
}
\;\;\;\; \mbox{;} \;\;\;\;
c=\frac{M_{e}}{M_{c}}\label{DM04Argris}%
\end{equation}


The model employs two ever--changing 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 so-called `image' back-stress 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}}[M-m]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=e-e_{c}$ is the state parameter and  $n^{b}$ and $n^{d}$ are material constants.%

\begin{figure}
[ptb]
\begin{center}
\includegraphics[height=5cm]{pi-plane.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 non-associative 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 non-negative 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{1-c}{c}g\cos{3\theta},\;\;C=3\sqrt{\frac{3}{2}}%
\frac{1-c}{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 so-called
fabric-dilatancy 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 stress-ratio
tensor $\alpha_{ij}$ and the fabric-dilatancy tensor $z_{ij}$. The evolution
law for the back stress-ratio 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}(1-c_{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
fabric-dilatancy 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 well-known 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 u--p--U 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 u-p-U
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 solid-fluid 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 1-D 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  $u-p-U$  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}{|c|c|c|}
\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 closed-form 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 closed-form 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}{|c|c|c|} \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 u--p--U 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}{|c|c|c|} \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_K10-2and10-6cms.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 Dafalias-Manzari model. }
\begin{center}
% {\normalsize
\begin{tabular}
[c]{|l|c|c||l|c|c|}\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 & Fabric-dilatancy & $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  stress-strain  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{u-p-U_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 multiple-elements 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}{|c|c|c|} \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 master--slave 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
so-called    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 $p-q$ 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
(elastic--plastic porous solid -- fluid) behavior of soils using $u-p-U$ formulation. A critical
state elasto--plastic 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
elasto--plasticity  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)]{local-87}
John Argyris and Hans-Peter 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
  622--634, 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 TR97-016 (UofF) and TR-97-046 (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 post--liquefaction site
  response.
\newblock \emph{Soil Dynamics and Earthquake Engineering}, 22:\penalty0
  259--271, 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
  883--905, 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 636--644, 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 399--413, 1995.

\bibitem[Gajo et~al.(1994)Gajo, Saetta, and Vitaliani]{Gajo1994}
A.~Gajo, A.~Saetta, and R.~Vitaliani.
\newblock Evaluation of three-- and two--field finite element methods for the
  dynamic response of saturated soil.
\newblock \emph{International Journal for Numerical Methods in Engineering},
  37:\penalty0 1231--1247, 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
  283--292, 1977.

\bibitem[Hughes and Liu(1978{\natexlab{a}})]{Hughes78a}
T.~J.~R. Hughes and W.~K. Liu.
\newblock Implicit--explicit finite elements in transient analysis: Stability
  theory.
\newblock \emph{{ASME} Journal of Applied Mechanics}, 45:\penalty0 371--374,
  June 1978{\natexlab{a}}.

\bibitem[Hughes and Liu(1978{\natexlab{b}})]{Hughes78b}
T.~J.~R. Hughes and W.~K. Liu.
\newblock Implicit--explicit finite elements in transient analysis:
  Implementation and numerical examples.
\newblock \emph{{ASME} Journal of Applied Mechanics}, 45:\penalty0 375--378,
  June 1978{\natexlab{b}}.

\bibitem[Hughes(1987)]{local-86}
{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 elastic--plastic finite element
  computations in geomechanics.
\newblock \emph{Communications in Numerical Methods in Engineering},
  17\penalty0 (2):\penalty0 115--125, January 2001.

\bibitem[Jeremi{\' c} and Sture(1997)]{Jeremic96b}
Boris Jeremi{\' c} and Stein Sture.
\newblock Implicit integrations in elasto--plastic geotechnics.
\newblock \emph{International Journal of Mechanics of Cohesive--Frictional
  Materials}, 2:\penalty0 165--183, 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 113--126, 1998.

\bibitem[Jeremi{\' c} and Yang(2002)]{Jeremic2000f}
Boris Jeremi{\' c} and Zhaohui Yang.
\newblock Template elastic--plastic computations in geomechanics.
\newblock \emph{International Journal for Numerical and Analytical Methods in
  Geomechanics}, 26\penalty0 (14):\penalty0 1407--1427, 2002.

\bibitem[Li and Wang(1998)]{Li98}
X.~S. Li and Y.~Wang.
\newblock Linear representation of steady-state line for sand.
\newblock \emph{Journal of Geotechnical and Geoenvironmental Engineering},
  124\penalty0 (12):\penalty0 1215--1217, 1998.

\bibitem[Manzari and Dafalias(1997)]{Manzari97}
M.~T. Manzari and Y.~F. Dafalias.
\newblock A critical state two--surface plasticity model for sands.
\newblock \emph{G{\' e}otechnique}, 47\penalty0 (2):\penalty0 255--272, 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 (67-94), 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 1--74, Laurel, Maryland,
  October 22-23 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 151--190, 1990.

\bibitem[Prevost(1985)]{Prevost85}
Jean~H. Prevost.
\newblock Wave propagation in fluid-saturated 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 Prentice-Hall, 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 0-913478-08-3.

\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 $u-w$ and $u-\pi$ finite element methods for the
  dynamic response of saturated porous media using one--dimensional models.
\newblock \emph{International Journal for Numerical and analytical methods in
  Geomechanics}, 10:\penalty0 461--482, 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 381--398, 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 three-field finite element procedure for analysis of elastic wave
  propagation through fluid-saturated soils.
\newblock \emph{Soil Dynamics and Earthquake Engineering}, 9:\penalty0 58--65,
  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 60--72,
  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 81--91, 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 71--96, 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 0-471-98285-7.

\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 385--395, 1980.

\bibitem[Zienkiewicz and Taylor(1991{\natexlab{a}})]{FEM-IV-1}
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}})]{FEM-IV-2}
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



