\documentclass[a4paper,11pt]{article}
%\documentclass[dvips]{article}
%\documentclass[12pt,letterpaper]{article}
%\usepackage{times}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{setspace}
\usepackage[normalsize]{caption}
\usepackage[numbers,round,colon,authoryear]{natbib}
%\usepackage{natbib}
\newcommand{\ud}{{\rm d}}
\usepackage{url}
\setlength{\textwidth}{17.5cm}
\setlength{\textheight}{26.8cm}
\setlength{\hoffset}{2.5cm}
\setlength{\voffset}{2.9cm}
%\usepackage{pause}
%\definecolor{webgreen}{rgb}{0, 0.15, 0} % less intense green
%\definecolor{webblue}{rgb}{0, 0, 0.15} % less intense blue
%\definecolor{webred}{rgb}{0.15, 0, 0} % less intense red
%\usepackage[pdfauthor={Boris Jeremic},
% colorlinks=true,
% linkcolor=webblue,
% citecolor=webred,
% urlcolor=webgreen,
% linktocpage,
% dvips]{hyperref}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ovo je predefinisanje datuma za cirilicu ( vidi Kopka93 strana 175 )
\def\today
{\number\day. \space \ifcase\month\or
January\or
February\or
March\or
April\or
May\or
June\or
July\or
August\or
September\or
October\or
November\or
December\fi,\space
\number\year}
%time of day as found in layout.tex
% 
%
% TIME OF DAY
%
\newcount\hh
\newcount\mm
\mm=\time
\hh=\time
\divide\hh by 60
\divide\mm by 60
\multiply\mm by 60
\mm=\mm
\advance\mm by \time
\def\hhmm{\number\hh:\ifnum\mm<10{}0\fi\number\mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\def\baselinestretch{2}
\def\baselinestretch{1.1}
%Mahdi
%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%\usepackage{lineno}
%\usepackage[left,pagewise,displaymath, mathlines]{lineno}
\usepackage[left,displaymath, mathlines]{lineno}
\begin{document}
\begin{center}
%{\Large \bf Simulation of Fully Coupled Behavior of Porous Media: Verification}
{\Large \bf
Solution Verification Procedures for
Modeling and Simulation of Fully Coupled Porous Media: Static and Dynamic Behavior}
\end{center}
\vspace*{1cm}
\begin{center}
{\large \bf Panagiota Tasiopoulou} \\
{PhD Student,
National Technical University of Athens,
Athens, Greece}\\ ~ \\
%Mahdi Taiebat \scalebox{1}{\includegraphics[viewport=0 4.5 0 0]{mahdi.eps}}
{\large \bf Mahdi Taiebat} \\
Associate Professor,
Department of Civil Engineering, The University of British Columbia,
Vancouver, BC, Canada\\ ~ \\
%%%%%%%%
{\large \bf Nima Tafazzoli} \\
Geotechnical Engineer, Tetra Tech EBA, Vancouver, BC, Canada \\ ~ \\
{\large \bf Boris Jeremi{\'c}}\footnote{Corresponding Author: Boris Jeremi{\'c}, \\
Department of Civil
and Environmental Engineering University of California, Davis, CA,
U.S.A. \texttt{jeremic@ucdavis.edu}, and \\
Earth Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA,
U.S.A. \texttt{bjeremic@lbl.gov}}
\\
Professor,
Department of Civil and Environmental Engineering,
University of California, Davis, CA, U.S.A. \\ and \\
Faculty Scientist,
Earth Science Division,
Lawrence Berkeley National Laboratory, Berkeley, CA, U.S.A.
\end{center}
%%%%%%%%
%and
%%%%%%%%%
%Yiorgos Perikleous\affilnum{4}
%\vspace*{2cm}
%\begin{center}
%{\sc Version: \today{} , \hhmm{} }
%\end{center}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage
%\tableofcontents
%\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\end{abstract}
%\linenumbers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Abstract}
%{\sc \Large I have updated it this}
Numerical prediction of dynamic behavior of fully coupled saturated porous media
is of great importance in many engineering problems. Specifically, static and
dynamic response of soils  porous media with pores filled with fluid, such as
air, water, etc.  can only be modeled properly using fully coupled approaches.
Modeling and simulation of static and dynamic behavior of soils require
significant Verification and Validation (V\&V) procedures in order to build
credibility and increase confidence in numerical results. By definition,
Verification is essentially a mathematics issue and it provides evidence that
the model is solved correctly, while Validation, being a physics issue, provides
evidence that the right model is solved.
%
This paper focuses on Verification procedure for fully coupled modeling and
simulation of porous media. Therefore, a complete Solution Verification suite
has been developed consisting of analytical solutions for both static and
dynamic problems of porous media, in time domain. Verification for fully coupled
modeling and simulation of porous media has been performed through comparison of
the numerical solutions with the analytical ones. Modeling and simulation is
based on the so called, $u$$p$$U$ formulation.
%\citep{Zienkiewicz84, Jeremic2004d,
%Jeremic2007e}.
Of particular interest are numerical dispersion effects which determine the
level of numerical accuracy.
%
These effects are
investigated in detail, in an effort to suggest a compromise between
numerical error and computational cost.
%
%In this paper, a solution verification suite (a subset of the whole
%verification suite, divided into code and solution verification processes),
%consisting of both static and dynamic problems, is developed for fully coupled
%modeling and simulation of porous media. Modeling and simulation is based on the
%so called, $u$$p$$U$ formulation \citep{Zienkiewicz84, Jeremic2004d,
%Jeremic2007e}.
%
%Static and dynamic examples are used to verify the numerical solutions
%through comparison with analytical solutions.
%
%Of particular interest are numerical dispersion effects.
%
%These effects are
%investigated in detail in an effort to suggest a compromise between
%numerical error and the computational cost.
\vspace*{0.3cm}
\noindent{\bf keywords:} verification, finite elements, fully coupled
analysis, porous media, numerical dispersion
%elasticplastic
%Note (mahdi): we are not dealing with elastic plastic materials here in this
%paper  our focus is rather only on elastic material.
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
\section{Introduction}
One of the main questions raised when results of a numerical analysis
are
evaluated is: ``How much can (should) we trust model
implementations
and how much can (should) we trust numerical simulations?"
Trust in numerical implementation is gained through a verification process,
while, trust in results obtained from numerical simulation, is gained through a
validation process \citep{Oberkampf2002, Roy2011}.
%
\begin{itemize}
\item[] {\bf Verification} is a process of determining that a model
implementation
accurately represents the developer's conceptual description and specification.
It is essentially a mathematics issue and it provides evidence that the model is
solved correctly.
%
\item[] {\bf Validation}, on the other hand, is a 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, and it provides evidence that
the correct model is solved.
%
\end{itemize}
\noindent
Verification and Validation (V\&V) are the primary means of assessing accuracy
in modeling and computational simulations in order to build confidence and
credibility in numerical predictions.
\citet{Oberkampf2002}
mainly focus on trying to accurately model the ``real world", aiming to provide
verified and validated
computer simulation of the reality.
%
Slightly different approach is taken by \citet{Oden2010a, Oden2010b}, where
the
actual purpose of numerical simulation is to facilitate the engineering
decision
process, rather than to explicitly represent the ``real world". Despite the
alternative approaches on the objectives of numerical modeling, it is important
to note that both schools of thought
place significant emphasis on the need
for
validation of numerical predictions under uncertainty. In an attempt to describe
how V\&V provide the connection between reality and numerical simulation
(computer implementation) in the process of reaching an engineering decision, a
flow chart has been constructed and shown in Figure \ref{RoleVV}.
%
\begin{figure}[!htb]
\begin{center}
{\includegraphics[width=7.0cm]{FIGS/RoleVV_NEW_knowledge.pdf}}
%\vspace*{0.5cm}
\caption{\label{RoleVV} Schematic representation of the role of Verification and
Validation on numerical modeling (inspired from \protect{\citet{Oberkampf2002,Roy2011}} and
\protect{\citet{Oden2010a}}).}
\end{center}
\end{figure}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%\paragraph{Verification Procedures.}
%\label{VandV_P}
It is evident that Verification and Validation procedures have separate but
complementary objectives by definition: i) the model is solved correctly
(Verification) and ii) the correct model is solved (Validation). On this basis,
they should be dealt as two distinct stages requiring different activities in
order to achieve them. In particular, two main activities related to
verification procedure can be described as \citep{Oberkampf2002}:
\begin{itemize}
\item {\bf Code Verification} which identifies and removes errors in computer
coding. It can further be separated into:
%
(a) {\bf numerical algorithm verification}, used to ensure that all
the numerical
algorithms are correctly implemented in the code, and
%
(b) {\bf software quality assurance}, used to ensure that the program
system is implemented reliably so that it produces repeatable
(stable) results with different compilers and on different computer
architectures.
\item {\bf Solution Verification} which quantifies numerical errors in
computed solutions. Numerical accuracy can only be estimated by direct comparison
with analytical solutions for a given application. Thus, evidence that the model can be solved correctly is provided.
\end{itemize}
On the other hand, Validation procedure involves numerical simulation of
realistic conditions and thus, direct comparison with experiments and/or well
documented case histories is needed. It is common in practice that only
validation procedures are performed, considered misleadingly to show enough
evidence that the numerical code is appropriate to solve correctly realistic
numerical models, while the verification stage is omitted. Therefore, it should
be emphasized that verification is the initial necessary step to ensure that
numerical solution of the model is accurate while validation is the final stage
to ascertain the numerical model can represent realistic conditions and thus,
numerical predictions can be trusted.
The current study focuses on the Verification stage and particularly, on the
solution verification procedure, assuming that code verification, which is a
strictly programming issue, has been successfully performed. The twofold goal of
this paper is to: a) provide a complete solution verification suite for fully
coupled, static and dynamic behavior of porous solidfluid systems, consisting
of analytical solutions for porous media in time domain, available in
literature, which can be used to verify any finite element/finite difference
code and b) to use this particular solution verification suite in order to
verify the UCD Computational Geomechanics Group's implementation of 3D
$u$$p$$U$ finite elements \citep{Jeremic2007e}, which is
available through an open source license. The available analytical solutions are
provided for idealized problems, based on assumptions such as linear elastic
behavior and simplified geometries (mostly 1D). Validation of coupled behavior
using elastoplastic material models, is the subject of another paper
\citet{Tasiopoulou2014b}.
It is important to note that this paper provides, to our knowledge, the
first verification of fully coupled (porous solid  pore fluid) modeling and
simulation. In the past, there were few published papers that used some of
the analytical and numerical solutions, used here, for comparison purposes.
For example a paper by \citet{Gajo1994}, published a while ago, focused on
comparing different methods, namely upU, uU and uw and used two closed
form solutions for that purpose. However, they did not provide details of
numerical algorithms used, since they used the same algorithm and parameters
(whatever they were) comparing different methods (as noted above, upU, uU
and uw). Focus of this paper is entirely on verification of upU method.
For this purpose all the available closed form solutions (to our knowledge)
are used and all the algorithms and parameters that we used are fully
documented. It can
be claimed that the upU formulation developed by \citet{Zienkiewicz84},
together with our implementation (presented verification, and validation provided by
\citep{Tasiopoulou2014b}), is now complete.
%It should be mentioned that the current verification process has been
%performed using the deterministic approach. While this might be appropriate
%for solution verification and code verification, validation must address the
%issue of uncertainty, as noted by \citet{Oden2010a, Oden2010b, Oden2011}.
%Such undertaking is planned for near future and will rely on Probabilistic
%ElasticPlasticity and Stochastic ElasticPlastic Finite Element Method
%\citep{Jeremic2005a, Sett2005a, Sett2009b, Sett2010a}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
% \clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%\section{TheoreticalNumerical Formulation
%%of Porous Media
%and Implementation
%Background}
%\label{Formulation}
%{\sc \Large THIS NEEDS TO BE UPDATED, it is from our old paper... The following
%subsection has been changed and the constitutive integration part has been
%omitted..}
%\subsection{Governing Equations of Porous Media}
\section{Governing Equations of Porous Media}
\label{Governing Equations of Porous Media}
The response of fully saturated soils under transient loads plays an important
role in geotechnical engineering. A literature review on the subject of
twophase materials identifies two major theories which have been developed
during the last century and are currently used, namely the Biot theory and
the theory of Porous Media . Biot's theory
%\citep{Biot1941,Biot1956a,Biot1956b,Biot1956c}
is initially based on
the work of \citet{Terzaghi23}, who proposed a model for
onedimensional
consolidation. Theory of Porous Media is based on the axioms of continuum theory of
mixtures extended by the concept of fractions by
\citet{Bowen1980,Bowen1982,Ehlers1993}.
\citet{Biot1941} presented a theoretical description of linear elastic porous
materials fully saturated by a viscous fluid under quasistatic conditions.
Extensions of his theory to anisotropic cases \citep{Biot1955} and to
poroviscoelasticity \citep{Biot1956a,Biot1956b,Biot1956c} followed. In Biot's
theory a fully saturated material is considered, which consists of two
compressible parts: the solid skeleton and the pore fluid. Thus, the solid
displacements, $u_i$ and the fluid displacements, $U_i$, are introduced. The
equations of motion are formed by stating the balance of momentum, or else the
dynamic equilibrium, for the solid skeleton and the pore fluid separately.
\citet{Zienkiewicz84} presented a slightly different formulation,
which is based on the concept of Biot's theory; however, it is more convenient
for numerical solution. Alternatively to the fluid displacements, a relative
displacement of fluid to the solid, $w_i$ (the socalled seepage displacement)
is used. Moreover, no additional mass density, $\rho_\alpha$, is introduced for
the dynamic interaction between the fluid and the skeleton, in contrast to
Biot's theory.
Formulation presented by \citet{Zienkiewicz84} is the most general one and is
used in our developments.
%
They did not, however, provide verification of
developed formulation that will prove that the formulation is implemented
accurately and that the examples are solved correctly.
%
Detailed solution verification of $upU$ formulation is precisely the goal of
this paper.
Fully coupled ($upU$) formulation is given below in
some detail, and is needed as reference for some of the verification solutions
given later in the paper.
%
We start from the relationship between effective stress,
total stress and pore pressure, written as $\sigma^{'}_{ij}= \sigma_{ij}+\alpha
\delta_{ij} p$, where $\sigma^{'}_{ij}$ is the effective stress
tensor\footnote{stresses are positive in tension}, $\sigma_{ij}$ is total stress
tensor, $\delta_{ij}$ is Kronecker delta, $p$ is the pore fluid pressure\footnote{pore pressure is positive in compression} 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, so that the relationship between total and effective stress becomes
$\sigma^{'}_{ij} = \sigma_{ij}+\delta_{ij} p$, which corresponds to the
classical effective stress definition by \citet{Terzaghi43}.
The material behavior of soil skeleton (porous solid) 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 on the soil skeleton. The
soil behavior ('mixture' 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 overall equilibrium (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
porous solid.
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 fluid density, respectively.
For the pore fluid, the equation of momentum balance can be written as
%%%
\begin{eqnarray}
p_{,i}  {R_i}  \rho_f \ddot{u}_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 can be written as
${R_i} = k_{ij}^{1} \dot{w}_j$,
where
$k_{ij}$ is the tensor of anisotropic intrinsic permeability coefficients.
For
simple case of isotropic permeability, scalar value of permeability $k$ is used,
so as ${R_i} = k^{1} \dot{w}_j$.
The permeability ${k}$ used here with dimension of [$L^{3}TM^{1}$] is
different from Darcy permeability (hydraulic conductivity), ($k_{D}$), which has
the dimension of velocity, i.e. [$LT^{1}$].
Their values are related by $k = k_{D} / g \rho_f$,
where $g$ is the gravitational acceleration and $\rho_f$ is the density of
the pore fluid.
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 (particles, not the skeleton) and fluid phases, respectively.
In the above governing equations, convective terms of lower order are
omitted. A change of variable is performed by introducing an
alternative variable $U_i$, defined as $ U_i = u_i + {w_i}/{n}$, that represents
absolute displacement of the pore fluid. Thus, the basic set of unknowns
comprises the soil skeleton displacements $u_i$, the fluid pore pressure $p$,
and the fluid displacements $U_i$. The set of the modified governing equations
is summarized below:
%
\begin{eqnarray}
\sigma_{ij,j}^{'}(\alphan) p_{,i}+(1n) \rho_s b_i(1n) \rho_s \ddot{u}_i +
n R_i=0
\label{34}
\end{eqnarray}
\begin{eqnarray}
n p_{,i}+n \rho_f b_in \rho_f \ddot{U}_i
 n R_i=0
\label{35}
\end{eqnarray}
\begin{eqnarray}
n \dot{U}_{i,i}=(\alphan) \dot{\varepsilon}_{ii}+\frac{1}{Q} \dot{p}
\label{36}
\end{eqnarray}
From the modified equation set (\ref{34}), (\ref{35}), and (\ref{36}), we can see that only
$\ddot{u}_i$ occurs in the first equation, and only $\ddot{U}_i$ in the second,
thus leading to a convenient diagonal form in discretization.
This theoretical approach describes the interaction of solid skeleton and pore
fluid with three major equations, governing both dynamic and quasistatic
phenomena: i) equilibrium of total solidfluid ``mixture" 
Eqs.~(\ref{upU_GEequilibrium}) and (\ref{34}), ii) equilibrium of pore fluid 
Eqs.~(\ref{GEequilibriumFluid}) and (\ref{35}), and iii) mass conservation of
the fluid flow  Eqs.~(\ref{GEconservationFlow}) and (\ref{36}). The boundary
conditions imposed on these variables will complete the problem. These equations
together with an appropriate constitutive relationship describe the behavior of
solid skeleton accounting for its full interaction with the pore fluid under
both quasistatic and dynamic conditions. Herein, the constitutive relationship
between increments of effective stress and strain of the solid skeleton is given in a general
incremental form:
%
\begin{equation}
d\sigma^{'}_{ij} = D_{ijkl} d\varepsilon_{kl}
\label{constitutive}
\end{equation}
%%%
where $\varepsilon_{ij} = (u_{i,j}+u_{j,i})/2$ is the small strain tensor of the
solid skeleton, and $D_{ijkl}$ is the tangent stiffness that can be elastic or
elastoplastic.
%In order to solve the equation of motion, the degrees of freedom must be determined. In retrospect, nowadays three general continuum formulations \citep{Zienkiewicz84} are common for modeling of the fully coupled problem %(soil skeleton  pore fluid) in geomechanics, namely the:
%
%\begin{enumerate}
%\item A combination of the pore pressure $p$ and the solid displacement $u_i$ with four unknowns in 3D, called $u$$p$ (twofield formulation).
%\item The solid displacement $u_i$ and the fluid displacement $U_i$ with six unknowns in 3D, called $uU$ (twofield formulation).
%\item A combination of the pore pressure $p$, the solid displacement $u_i$ and the fluid displacement $U_i$ with seven unknowns in 3D, called $u$$p$$U$ (threefield formulation).
%\end{enumerate}
%
%
%The $u$$p$ formulation captures the movements of the soil skeleton and the change of the pore pressure. Despite the fact that it is the most simplistic one of the three mentioned above, it can be valuable when the movements of %the water are not important. This formulation neglects the relative 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). This formulation must rely on Rayleigh damping to model velocity proportional energy dissipation (damping). 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 $u$$p$$U$ formulation considers compressible pore fluid and solid grains.
Additionally, due to computation of fluid inertia forces, this formulation is
applicable to any range of frequency for which the hypothesis of constant
permeability is valid \citep{Gajo1994}. Moreover, 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 (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.
An important advantage of $u$$p$$U$ formulation over commonly used
$u$$p$ formulations is that velocity proportional viscous damping is
introduced directly through the damping tensor which is a function of porosity
and permeability of the soil skeleton. This viscous damping provides for
physically based, velocity proportional energy dissipation, stemming from the
interaction of pore fluid and the porous solid (soil) skeleton. In addition to
these advantages, the inclusion of both solid skeleton and pore fluid
displacements in the field of unknowns allows for independent treatment of
accelerations of both constituents (skeleton and fluid) which improves accuracy
of simulations.
%
Some details of the finite element formulation and implementation
is given in the appendix for reference.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\section{Solution Verification Suite}
\label{Verification_Suite}
Solution verification of the $u$$p$$U$ formulation and its
implementation is presented in the following. Analytical solutions for porous
media, found in literature, are used to verify numerical solutions obtained by
application of the developed code. The analytical solutions in time domain,
available in literature, have been developed for linear elastic porous media and
simplified geometries (mostly 1D) and they are divided into two main categories,
based on: a) consolidation theory, where no inertia effects are taken into
account, called quasistatic poroelasticity, and b) dynamic version of the
twophase theory, called poroelastodynamics. Characteristics solutions from both
categories are chosen in the following to verify the $u$$p$$U$ finite
element formulation, as developed by \cite{Jeremic2007e}.
It should be mentioned that the majority of the analytical solutions are
provided for idealized examples, such as 1D linear elastic wave propagation,
which are not always representative of realistic cases in terms of 3D geometry
and nonlinearity. However, while they would seem to be inadequate for
Validation purposes, which require realistic experiments and/or well documented
case histories (see \citet{Tasiopoulou2014b} for more details), they successfully
serve the mathematical objective of Verification, ensuring the accuracy of the
numerical solution provided by the code.
\subsection{Quasistatic Poroelasticity}
\label{Quasistatic Poroelasticity}
This section focuses on analytical solutions based on coupled theory of linear
porous media under quasistatic conditions \citep{Coussy95,Coussy2004}.
Eqs.~(\ref{upU_GEequilibrium}) and (\ref{GEequilibriumFluid}), which are based
on equilibrium of total ``mixture" and fluid respectively, are reformed in order
to apply for quasistatic phenomena (inertia terms are omitted). Thus,
Eq.~(\ref{upU_GEequilibrium}) is written as:
\begin{equation}
\sigma_{ij,j} = 0
\label{equil_mix_quasi}
\end{equation}
and Eq.~(\ref{GEequilibriumFluid}), assuming isotropic permeability, is written as:
\begin{eqnarray}
p_{,i}  \frac{\dot{w}_i}{k} +\rho_f b_i = 0
\label{equil_Fluid_quasi}
\end{eqnarray}
Differentiating Eq.~(\ref{equil_Fluid_quasi}) with respect to spatial coordinates, we obtain:
\begin{eqnarray}
\dot{w}_{i,i} = k p_{,ii}
\label{equil_Fluid_quasi_darcy}
\end{eqnarray}
Combining Eqs.~(\ref{equil_Fluid_quasi_darcy}) and
(\ref{GEconservationFlow}), the mass conservation equation becomes:
\begin{eqnarray}
\alpha \dot{\varepsilon}_{ii} + \frac{\dot{p}}{Q} = k p_{,ii}
\label{mass_cons_quasi}
\end{eqnarray}
A simple linear constitutive relationship between effective stress and strain is
considered,
\begin{eqnarray}
\sigma_{ij} = (K  \frac{2}{3}\mu)\varepsilon_{v}\delta_{ij}+2\mu\varepsilon_{ij}+\alpha p\delta_{ij}
\label{shear_strain}
\end{eqnarray}
where $K$ is the bulk modulus, $\mu$ is the shear modulus, and $\varepsilon_{v}$ is
the volumetric strain of porous solid skeleton. Combination of
Eqs.~(\ref{shear_strain}) and (\ref{equil_mix_quasi}), and subsequent
differentiation with respect to spatial coordinates, leads to $(K + (4/3)\mu)\varepsilon_{v,ii} 
\alpha p_{,ii} = 0$, which can be solved for $\varepsilon_{v}$ as:
\begin{eqnarray}
\varepsilon_{v} = \frac{\alpha p}{K + (4/3)\mu}
\label{strain_p}
\end{eqnarray}
Lastly, differentiating $\varepsilon_{v}$ from Eq.~(\ref{strain_p}) with respect to time and substituting it in Eq.~(\ref{mass_cons_quasi}) in combination with the relationship $K_{u}=K +\alpha^{2}Q$ (where $K_u$ is the undrained bulk modulus of soil), the diffusion equation is derived as:
\begin{eqnarray}
\dot{p} = c_{f} p_{,ii}
\label{diffusion}
\end{eqnarray}
where $c_f$ is the fluid diffusivity coefficient, defined as:
\begin{eqnarray}
c_{f} = kQ\frac{K+(4/3)\mu}{K_{u}+(4/3)\mu}
\label{c_f}
\end{eqnarray}
\citet{Coussy95,Coussy2004} has presented analytical solutions for the diffusion
Eq.~(\ref{diffusion}) for different boundary condition problems which have been
used as verification examples and are studied in the following.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Vertical Consolidation of a Soil Layer}
\label{Vertical_Consolidation_of_a_Soil_Layer}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Analytical Solution by \citet{Coussy2004}} 
%
%
%he consolidation process can be defined as follows: When an elastic
%soil layer is subjected to an external change in mean normal stress, immediately
%the water will alone sustain this increment of mean normal stress and cause the
%buildup the excessive pore water pressure. In the progress of the flow of the
%water to the surface, the load is gradually transferred to the soil skeleton and
%the excessive pore water pressure will dissipate. At the same time, the
%settlement of the soil layer occurs. As settlement is usually a major concern in
%geotechnical engineering, this is a key problem in soil mechanics.
%
A soil layer of thickness $h$ in the $z$ direction, composed of isotropic,
homogeneous and saturated poroelastic material is assumed, as shown in
Figure~\ref{num_model}(a). The top surface of the soil layer, $z=0$, is drained,
while its base, $z=h$, is rigid and impervious.
%
\begin{figure}[!htb]
\begin{LARGE}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.40\textwidth]{FIGS/cons_mod1.pdf}
\vspace*{0.3cm}
\caption{\label{num_model}{(a) The geometry of the soil layer of thickness $h$ under uniform constant vertical pressure $\varpi$ = 400~kPa applied at the surface, and (b) the finite element mesh with the appropriate boundary conditions and loading at top.}}
\end{center}
\end{sffamily}
\end{LARGE}
\end{figure}
%
Since this is a onedimensional problem, the only nonzero displacement is the
vertical one, $u_z$. In particular, the fluid pore pressure, $p$, as well as the displacement, $u_z$, depend only
on $z$ and $t$. Thus, the only nonzero strain component is $\varepsilon_{zz}= {\partial u_z}/{\partial z}$.
The hydraulic boundary conditions require that $p= 0$ for the soil surface ($z= 0$) and ${\partial p}/{\partial z}=0$ for the base of the layer ($z=h$).
The undeformability of the substratum requires that the vertical soil displacement, $u_z$, is zero ($u_z=0$ for $z=h$).
%\begin{center}
%\begin{eqnarray}
%\boldsymbol u=u(z,t)e_z \qquad p=p(z,t)
%\label{99}
%\end{eqnarray}
%\end{center}
%\begin{center}
%\begin{eqnarray}
%\varepsilon_{zz}=\frac{\partial u}{\partial z}
%\label{1099}
%\end{eqnarray}
%\end{center}
%\begin{center}
%\begin{eqnarray}
%z=0 \qquad p= 0 ; \qquad z=h \qquad \frac{\partial p}{\partial z}=0
%\label{95}
%\end{eqnarray}
%\end{center}
%\begin{center}
%\begin{eqnarray}
%z=h \qquad u=0
%\label{97}
%\end{eqnarray}
%\end{center}
At time $t=0$, the instantaneous vertical constant load $\varpi$ is suddenly applied
on the top surface $z=0$. The equilibrium requires that, for $t>0$ and $z=0$, $\sigma_{zz}=\varpi$.
%\begin{center}
%\begin{eqnarray}
%t>0 \rightarrow \qquad z=0 \qquad \sigma_{zz}=\varpi
%\label{96}
%\end{eqnarray}
%\end{center}
Thus, the new fields of fluid pressure and
displacement, induced by the external loading, remain to be determined.
%
Application of the constitutive Eq.~(\ref{shear_strain}) gives:
\begin{eqnarray}
\frac{\partial u_z}{\partial z}=\varepsilon_{zz}=\frac{\alpha p\varpi}{K+(4/3)\mu}
\label{1101}
\end{eqnarray}
The instantaneous response of the porous medium to any external loading is undrained and can be expressed in the form:
\begin{eqnarray}
\varepsilon(t=0_+)=\frac{1}{\alpha M}p(t=0_+)
\label{1103}
\end{eqnarray}
%
which, when combined with Eq.~(\ref{shear_strain}), leads to:
\begin{eqnarray}
p(z,t=0_+)
=
\frac{\alpha M\varpi}{K_u+(4/3)\mu}
\qquad
\frac{\partial \ud_z}{\partial z}(z,t=0_+)
=
\varepsilon_{zz}
=
\frac{\varpi}{K_u+(4/3)\mu}
\label{1104}
\end{eqnarray}
Finally, in case of onedimensional consolidation, the diffusion Eq.~(\ref{diffusion}), needed to be solved, reads:
\begin{eqnarray}
t>0 \qquad \frac{\partial p}{\partial t}=c_f\frac{\partial^2p}{\partial z^2}
\label{102}
\end{eqnarray}
Satisfying the above mentioned boundary conditions, the fluid pressure can be expressed in the form of infinite series:
\begin{eqnarray}
p(z,t)=\frac{\alpha M\varpi}{K_u+(4/3)\mu}\sum_{m=0}^\infty\frac{4}{\pi(2m+1)}\sin\left(\frac{(2m+1)\pi}{2}\frac{z}{h}\right)\exp\left(\frac{(2m+1)^2\pi^2}{4}\frac{t}{\tau}\right)
\label{103}
\end{eqnarray}
Each term of the series decreases exponentially with respect to the ratio
${t}/{\tau}$, with $\tau$ being the characteristic consolidation time, defined as:
\begin{eqnarray}
\tau=\frac{h^2}{c_f}
\label{104}
\end{eqnarray}
By substituting Eq.~(\ref{103}) into Eq.~(\ref{1101}), the series
converges and, after integration term by term, we obtain:
\begin{eqnarray}
u_z(z,t) =
\nonumber \\
\frac{\varpi(hz)}{K+(4/3)\mu} +
\left(\frac{\varpi h}{K_u+(4/3)\mu}\frac{\varpi h}{K+(4/3)\mu}\right)
\nonumber \\
\sum_{m=0}^\infty\frac{8}{\pi^2(2m+1)^2}
\cos\left(\frac{(2m+1)\pi}{2}\frac{z}{h}\right)
\exp\left(\frac{(2m+1)^2\pi^2}{4}\frac{t}{\tau}\right)
%\nonumber \\
\label{106}
\end{eqnarray}
%
Using Eq.~(\ref{106}) and substitute $z=0$, the total settlement of the
soil layer, $s$, can be expressed as:
%
\begin{eqnarray}
s(t)=s_\infty+(s_{0^+}s_\infty)\sum_{m=0}^\infty\frac{8}{\pi^2(2m+1)^2}exp\{\frac{(2m+1)^2\pi^2}{4}\frac{t}{\tau}\}
\end{eqnarray}
%
where
%
\begin{eqnarray}
s_{0^+}=\frac{h\varpi}{K_u+(4/3)\mu}
\;\;\; \mbox{;} \;\;\;
s_\infty=\frac{h\varpi}{K+(4/3)\mu}
\label{107}
\end{eqnarray}
%
{\bf Numerical Analysis} 
%
A soil column of $u$$p$$U$, 8node brick, finite elements is used to
model the horizontal layer. Each element has dimensions
1~m$\times$1~m$\times$1~m, resulting in soilcolumn mesh, 10 m high, as
illustrated in Figure~\ref{num_model}(b). The elastic material properties, shown
in Table~\ref{Soil Properties for quasistatic problems}, are chosen as
indicative values for the natural soil deposit. A constant uniform vertical
pressure of 400~kPa is instantaneously applied on the top of the finite element
mesh.
%
\begin{table}[!htb]
\begin{center}
\caption{ \label{Soil Properties for quasistatic problems}
Soil properties used in numerical analysis of quasistatic problems.}
%\begin{spacing}{1.4}
%\setlength{\extrarowheight}{3pt}
\begin{tabular}{lccc}
Parameter & Symbol [Units] & 1D consolidation & Line injection \\
\hline
gravity acceleration & $g$ [m/s$^2$] & $9.81$ & $9.81$ \\
soil Young's Modulus & $E$ [kN/m$^2$] & $10^4$ & $1.2\times 10^6$ \\
soil Poisson's ratio & $v$ & $0.25$ & $0.2$ \\
solid density & $\rho_s$ [ton/m$^3$] & $2.65$ & $2.7$ \\
fluid density & $\rho_f$ [ton/m$^3$] & $1.0$ & $2.7$ \\
solid bulk modulus & $K_s$ [kN/m$^2$] & $37\times 10^6$ & $3.6\times 10^7$ \\
fluid bulk modulus & $K_f$ [kN/m$^2$] & $2.2\times 10^6$ & $1.0\times 10^{17}$ \\
porosity & $n$ & 0.46 & 0.4 \\
Biot coefficient & $\alpha$ & $1.0$ & $1.0$\\
Darcy permeability & $k_D$ [m/s] & $10^3$ & $3.6\times 10^{6}$ \\
%\hline
soil drained bulk modulus & $K$ [kN/m$^2$] & $6.67\times 10^3$ & $6.7\times 10^5$ \\
soil undrained bulk modulus & $K_u$ [kN/m$^2$] & $4\times 10^6$ & $6.0\times 10^7$ \\
% soil shear modulus & $\mu$ & $4\times 10^3$~kN/m$^2$ \\
fluid diffusivity coefficient & $c_f$ [m$^2$/s] & $1.2$ & $0.5$ \\
% characteristic consolidation time & $\tau$ & $83.3$~s \\
\end{tabular}
%\end{spacing}
\end{center}
\end{table}
%
The following boundary conditions are applied to the model, shown in
Figure~\ref{num_model}(b): the solid and fluid displacements are fixed at the
base and the pore pressure is kept zero at the top surface. In order to simulate
the 1D consolidation problem, the lateral movement of the solid and fluid phase
is suppressed so that the vertical displacement remains the only nonzero
displacement. To remedy any artificial oscillation due to spatial and time
discretization, numerical damping is introduced into the analysis by using
$\gamma=0.6$ and $\beta=0.3025$ for the Newmark time integrator
\citep{Newmark1959,local87}.
The instantaneous response of the soil skeleton is undrained which leads to
sudden increase of pore pressure. The boundary condition of perfect drainage at
the surface ($p=0$) causes an upward flow of the pore water, so that pore
pressure gradually reduces with time, as depicted by
Figure~\ref{dissipation}(a). The time factor, $T_v$ is defined as $T_v=t/\tau$
so that
for $T_v=1$, the process of consolidation is considered to have been completed,
since the rate of dissipation has been practically diminished.
Figure~\ref{dissipation}(b) depicts the distribution of pore pressure with depth
for
different moments in time, $T_v$, confirming the upward movement of pore water,
as the rate of pore pressure dissipation is clearly greater for depths closer to
the surface.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.9\textwidth]{FIGS/cons_pp.pdf}
\vspace*{0.5cm}
\caption{\label{dissipation}{Comparison between numerical analysis and
analytical solution by \citet{Coussy2004}: a) normalized pore pressure at
various depths with respect to real time, $t$, and b) distribution of normalized
pore pressure with normalized depth, $z/h$, for different time factors, $T_v$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
The gradual upward movement of pore water leads to increase of the
vertical soil displacement with time, as shown in Figure~\ref{settle}(a). It can
be concluded that the numerical analysis can effectively demonstrate the process
of the dissipation of the pore pressure in agreement with the analytical
solution by \citet{Coussy2004}.
%
\begin{figure}[tb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.4\textwidth]{FIGS/cons_disp.pdf}
\vspace*{0.5cm}
\caption{\label{settle}{(a) Comparison between numerical analysis and analytical
solution by \citet{Coussy2004} in terms of time histories: a) solidskeleton
displacements, and b) fluid displacements, for various depths.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Moreover, the numerical analysis based on $u$$p$$U$ formulation permits the
calculation of upward fluid movement, $U_z$, plotted with time in
Figure~\ref{settle}(b). The upward displacement of the fluid, $U_z$, at the end
of
consolidation has been computed as 0.39~m, more than the final downward
displacement of the soil skeleton (0.33~m). Thus, an accurate estimation of
the change in porosity due to settlement of the soil and consequent upward
drainage of amount of pore fluid can be performed. The amount of fluid that
escaped is calculated as $\Delta V_f=nU_z=0.46\times 0.39\text{~m}\times 1\text{~m}\times
1\text{~m}=0.1794\text{~m}^3$. The change of solidskeleton volume is equal to $\Delta
V_s=(1n)s=0.54\times 0.33\text{~m}\times 1\text{~m}\times 1\text{~m}=0.1782\text{~m}^3$. If the total initial
volume of soil is $V=10\text{~m}^3$ and the initial volume of pore fluid is $nV=4.6\text{~m}^3$,
the final porosity at the end of consolidation is estimated
$n'={(4.60.1794)}/{(100.1782)}=0.45$, less than the initial one, $n=0.46$, indicating slight compaction of the soil layer.
%It should be mentioned
%that theory does not account for change of porosity throughout the consolidation
%process since it does not involve the fluid displacements.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Line Injection of Fluid in a Reservoir}
\label{Line Injection of a fluid in a Reservoir}
%{\sc \Large I revised everything, for some reason everything was quite out of place, even eqns. Keep in mind when using the thesis you have taken this from...}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Analytical Solution by \citet{Coussy2004}} 
%
A reservoir of infinite extent composed of an isotropic, homogeneous
and saturated poroelastic material is considered. Through a cylindrical well of negligible
dimensions, fluid mass injection is performed in all directions
orthogonal to the well axis of symmetry, forming the $z$ axis of coordinates.
The flow rate of fluid mass injection is constant and equal to $q$. This is a problem of cylindrical symmetry, requiring the use of cylindrical
coordinates ($r$,$\theta$,$z$). The vectors of relative flow of fluid
mass, $\boldsymbol w$, and solidskeleton displacement, $\boldsymbol \xi$ are expressed as $\boldsymbol w=w(r,t)\boldsymbol {e_r}$ and $\boldsymbol \xi=\xi(r,t)\boldsymbol {e_r}$,
%
%
%\begin{center}
%\begin{eqnarray}
%\boldsymbol w=w(r,t)e_r \qquad \boldsymbol{\xi}=\xi(r)e_r
%\label{108}
%\end{eqnarray}
%\end{center}
%
where $\boldsymbol{e_r}$ is the unit vector along the radius. Solidskeleton strain components are derived in the form:
\begin{eqnarray}
\varepsilon_{rr}=\frac{\partial\xi}{\partial r}
\;\;\; \mbox{;} \;\;\;
\varepsilon_{\theta\theta}=\frac{\xi}{r}
\;\;\; \mbox{and} \;\;\;
\varepsilon_{ij}=0
\label{73}
\end{eqnarray}
Thus, Eqs.~(\ref{strain_p}) and (\ref{diffusion}) are reformed in cylindrical coordinates:
\begin{eqnarray}
\varepsilon_v=\frac{1}{r} \frac{\partial (r\xi)}{\partial r}= \frac{\alpha p}{K + (4/3)\mu}
\label{strain_p1}
\end{eqnarray}
\begin{eqnarray}
\frac{\partial p}{\partial t}=c_f \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial p}{\partial r})
\label{diffusion1}
\end{eqnarray}
Requiring that the fluid flow reduces to zero infinitely far from the
well, $rw \rightarrow 0$ as $r \rightarrow \infty$ , the fluid mass balance relationship is applied:
\begin{eqnarray}
%rw \rightarrow 0 \qquad r \rightarrow \infty \qquad t<\infty \nonumber\\
\int_{0}^{\infty}\dot{w}_{,r}(r,t)2\pi rdr=q \qquad
\label{110}
\end{eqnarray}
%
The solution of the set of Eqs.~(\ref{strain_p1}), (\ref{diffusion1}), (\ref{110}), and (\ref{GEconservationFlow}) for the instantaneous injection of a finite volume of fluid, $\Omega$, is given in the following form:
%
\begin{eqnarray}
p=Q\frac{K + (4/3)\mu}{K+u + (4/3)\mu}\frac{\Omega}{4\pi c_{f} t}\exp(\frac{r^2}{4c_ft}) \nonumber
\end{eqnarray}
\begin{eqnarray}
\xi=\frac{\alpha Q\Omega}{2\pi r(K_u +(4/3)\mu)}\left[1\exp(\frac{r^2}{4c_ft})\right]
\label{111}
\end{eqnarray}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Numerical Analysis} 
%
As a result of axisymmetry of the geometry, the model has been
developed for $\pi/2$ segment, 5~cm thick, with radius equal to 1~m. A
cylindrical well is drilled at the center of the model, with 1~cm radius, having
negligible size when compared to the model radius, but still providing
wellshaped elements.
%The final mesh is shown in Figure~\ref{fig32}.
%%
%\begin{figure}[!htb]
%\begin{center}
%\includegraphics[width=6cm]{FIGS/Inject1.pdf}
%\vspace*{0.5cm}
%\caption{\label{fig32} {Mesh generated for line injection problem.}}
%\end{center}
%\end{figure}
%%
Due to
planestrain conditions, all the movements for solid and fluid phases are
suppressed along the vertical axis of symmetry, $z$. The $u$ and $U$ nodes along
the outside perimeter are fixed. It should be noted that $\Omega$ is the volume
of the fluid injected per unit of vertical well length, m$^3$/m. In order to
generate the volume of 1~cm$^3$/m, the corresponding fluid displacement of the
nodes along the well has been calculated and applied as a step function at the
time $t=0$~sec. For simplicity, the initial fluid pressure $p_0$ is set to be
0~kPa. The analytical solution is studied below using the set of parameters
shown in Table~\ref{Soil Properties for quasistatic problems}.
In the analysis, the pore pressure and the radial displacements are recorded at
three points on the radius (10~cm, 50~cm, and 85~cm). The analytical and
numerical results compare well with each other as indicated by
Figure~\ref{fig33}.
%
\begin{figure}[~htb]
\begin{center}
\includegraphics[width=6cm]{FIGS/Injection.pdf}
%\includegraphics[width=8cm]{FIGS/InjectionPressure_new_ps1.pdf}
%\hfill
%\includegraphics[width=7.8cm]{FIGS/InjectionDisplacement_new_ps1.pdf}
\vspace*{0.3cm}
\caption{\label{fig33} {Comparison between numerical analysis and analytical solution by \citet{Coussy2004}: a) radial displacement, and b) pore pressure versus time at different distances, R, from the center of the model.}}
\end{center}
\end{figure}
%
The time step, $\Delta t$, was set equal to 1~sec; thus, the
first data point is recorded at the time of 1~sec.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Poroelastodynamics}
\label{Poroelastodynamics}
One of the significant findings by \citet{Biot1956a,Biot1956b,Biot1956c} was the
identification of two kinds of compressional waves and one shear wave, which has
been the starting point for further research on the mechanics of dispersive
waves through analytical or numerical approaches. Several analytical and
semianalytical solutions for dynamic poroelasticity are available in the
literature, varying with geometry, type of loading, time domain or transformed
domain solution, assumptions related to the compressibility of the fluid and
the formulation used ($u$$p$, $u$$p$$U$, $u$$U$).
%
A review of this wide range of solutions in dynamic poroelasticity is nicely
presented by \citet{Schanz2009}. The current paper focuses on solutions
regarding onedimensional wave propagation in saturated linear elastic porous
media due to step and sinusoidal excitations.
In order to study the accuracy of $u$$p$$U$ finite element procedure in
representing a dispersive wave propagation in saturated soil, numerical
dispersion effects should be taken into account. Dependence of the wavevelocity
on the frequency content of the excitation is the basic concept behind the
theory of dispersive waves. Thus, when wave propagation is modeled in a
numerical scheme, numerical features, such as the element size, integration
scheme, numerical damping, element type etc., should be chosen appropriately, so
that the frequency range of the excitation will be represented in the model.
Otherwise, the wavevelocities developed in the numerical model will diverge
from the theoretical solution due to an artificial dispersion attributed to
numerics and not physics (numerical dispersion).
Various theoretical works concern the effect of numerical dispersion
\citep{Deraemaeker1999,local86,Ihlenburg1995,Semblat2000}, highlighting
spatial discretization as one of the most influential numerical features. In
particular, the effect of numerical dispersion decreases with the increase in wave length,
$\lambda$, and increases with the increase in mesh size $\Delta h$, which can lead to the
necessity to use fine meshes, even uneconomical ones for small wave lengths
\citep{Deraemaeker1999}. A common rule of thumb for typical finite element
procedures is to resolve the wavelength by 10 elements \citep{Ihlenburg1995, local86};
thus the element size, $\Delta h$, should be less that 1/10 of the smallest
wavelength ($\lambda_{\min}$) involved in the application:
\begin{eqnarray}
\Delta h \le \frac{\lambda_{\min}}{10}
\mbox{~~~~~~~~~where~~~~~~~~}
\lambda_{\min}= \frac{V_{\min}}{f_{\max}}
\label{110a}
\end{eqnarray}
%
where $V_{\min}$ is the lowest wave velocity of interest in the
simulation and generally, it is a function of the elastic modulus of the medium,
while $f_{\max}$ is the maximum frequency that is modeled.
Spatial and temporal discretization are linked in a way requiring that mesh
spacing is accompanied by an appropriate time interval. As a wave front
progresses in space it reaches one point after the other. If the time step in
the finite element analysis is too large, the wave front can appear to reach two
consecutive nodes at the same moment. This would violate a fundamental
property of wave propagation and can lead to instability. In order to ensure
that the propagation of wave can be modeled properly, the required time
step, $\Delta t$, must be limited to:
%
\begin{eqnarray}
\Delta t \le \frac{\Delta h}{V_{\max}}
\label{111a}
\end{eqnarray}
%
where $V_{\max}$ is the highest wave velocity involved, usually being a function
of the undrained elastic modulus in case of fully saturated porous media.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Step displacement excitation}
%\label{One dimensional shock wave propagation with a step displacement boundary condition at the surface}
\label{1DDisp_BC}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Analytical Solution by \citet{Gajo1995b}} 
%
An onedimensional analytical solution of the Biot's equations (using $u$$U$
field displacements) is provided by \citet{Gajo1995b} for a general transient
problem in poroelasticity. The analytical solution was obtained using Fourier
series expansion and it is not based on any assumptions with respect to the
inertial, viscous or mechanical coupling. Furthermore, it is applicable to any
type of boundaryinitial value problem.
Since each term of the Fourier series represents a frequency component of the
excitation signal, the analytical solution can describe the behavior of every
single frequency component related to the problem. Thus, it can illustrate the
mechanics of dispersive wave propagation in which higher frequencies propagate
with two waves and lower frequencies with only one wave, as a function of
permeability and travel length. Considering the above, the analytical solution
can provide a useful comparative tool towards the verification of the existing
numerical solutions based on the finite element method \citep{Gajo1994}.
%
\citet{Gajo1995b} present the response of the porous medium only relative to the
first arrival of the waves of the first and second kind. Specifically,
analytical results are shown for a finite soil column subjected to a step
displacement boundary condition (Heaviside function) at the surface applied both
to solid and fluid phases, as shown in Figure~\ref{Fig1001}(a).
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.40\textwidth]{FIGS/Gajo_mod.pdf}
\caption{\label{Fig1001}{a) A representative soil column subjected to a step
vertical displacement equal to $1.0\times10^{3}\text{~cm}$ at the surface, and
b) the finite element mesh and the applied boundary conditions used for
numerical modeling.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
This problem can
demonstrate better the mechanics of dispersive wave propagation, since the step
excitation contains waves of all frequencies. In the following section (3.2.2),
analytical solution is compared with numerical analysis for a soil column
subjected to a single sine pulse with characteristic predominant frequency of
$10^6\text{~Hz}$, producing smoother wavefronts.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Numerical Analysis} 
%
Fine spatial discretization is required due to step boundary condition. This
kind of excitation contains waves of all frequencies leading to meshes with
a large number of finite elements, according to Eq.~(\ref{110a}), which
suggests that the appropriate
element size is inversely proportional to the maximum frequency introduced to
the model. In fact, there is no element size that could satisfy the equation;
thus, numerical dispersion is inevitable. It can only be diminished by
increasing the spatial refinement. Evidently, a compromise needs to be made
between the quantity of numerical error and the computational cost.
Numerical models with three different values of permeability  resulting in
three different viscous coupling scenarios ($k=10^{8}\text{~cm$^3$s/g}$,
$k=10^{6}\text{~cm$^3$s/g}$, $k = 10^{5}\text{~cm$^3$s/g}$)  were solved in
order to verify the $u$$p$$U$ formulation in a wide range of drag forces. At
the top surface of the model, a vertical step displacement of
$1.0\times10^{3}\text{~cm}$ is applied both to the solid and the fluid phases.
No lateral flow or displacement is allowed. The base of the model is rigid and
impervious, as schematically illustrated in Figure~\ref{Fig1001}(a).
According to \citet{Gajo1995}, the two extreme values of wave velocity for
longitudinal propagation in a porous medium with properties shown in the third
column of Table~\ref{Soil Properties_gajo}, are: $V_{\min}=610\text{~m/s}$ and
$V_{\max}=1819\text{~m/s}$, corresponding to the case of significantly low
viscous coupling, where two longitudinal waves exist (one slow and one fast).
%
\begin{table}[!htb]
\begin{center}
\caption{ \label{Soil Properties_gajo}
Soil Properties for 1D wave propagation due
to various boundary conditions and loading.}
%\begin{spacing}{1.4}
\begin{tabular}{lcccc}
Parameter & Symbol [Units] & Step/Sine pulse disp. & Step/Sine loading & Step vel. \\
\hline
gravity acceleration & $g$ [m/s$^2$] & $9.81$ & $9.81$ & $9.81$ \\
soil Young's Modulus & $E$ [kN/m$^2$] & $1.2\times 10^6$ & $20\times 10^3$ & $23.2\times 10^6$ \\
soil Poisson's ratio & $v$ & $0.3$ & $0.2$ & $0.17$ \\
solid density & $\rho_s$ [ton/m$^3$] & $2.7$ & $2.0$ & $2.66$ \\
water density & $\rho_f$ [ton/m$^3$] & $1.0$ & $1.0$ & $1.0$ \\
solid bulk modulus & $K_s$ [kN/m$^2$] & $36\times 10^6$ & $36\times 10^6$ & $36\times 10^6$ \\
fluid bulk modulus & $K_f$ [kN/m$^2$] & $2.2\times 10^6$ & $2.2\times 10^6$ & $2.2\times 10^6$ \\
porosity & $n$ & $0.4$ & $0.33$ & $0.18$ \\
Biot coefficient & $\alpha$ & $1.0$ & $1.0$ & $0.677$ \\
\end{tabular}
%\end{spacing}
\end{center}
\end{table}
%
For infinitely large viscous coupling, the velocity of the single longitudinal
wave is $1773\text{~m/s}$. The chosen finite element mesh consists of 400
$u$$p$$U$ brick finite elements of dimensions $0.01\text{~cm} \times
0.01\text{~cm} \times 0.01\text{~cm}$, creating a soil column 4 cm high, as
illustrated in Figure~\ref{Fig1001}(b). The maximum frequency allowed to be
present in a model with the abovementioned grid spacing is $18.19
\times{10^6}\text{~Hz}$ for the first longitudinal wave (fast) and $6.1
\times{10^6}\text{~Hz}$ for the second one (slow). There is a cutoff of greater
frequencies imposed by the specific computational mesh.
The maximum time step is equal to
$(0.0001\text{~m})/(1819\text{~m/s})=5.5\times{10^{8}}\text{~sec}$, as
suggested by Eq.~(\ref{111a}). The selected temporal integration involves 800
steps of $\Delta t = 2.0\times10^{8}\text{~sec}$, which allows a maximum wave
velocity of $5.0\times10^5\text{~m/s}$.
% In addition, the Newmark time
%integration method was used \citep{Newmark1959}.
% Sets of parameters, assuming
%unconditionally numerical stability, were chosen for both integrators.
Figure~\ref{Fig13} shows the comparative results for all three different values
of viscous coupling using Newmark parameters: $\gamma = 0.6$ and $\beta =
0.3025$, which introduce numerical damping in the model.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Gajo_comp.pdf}
\vspace*{0.5cm}
\caption{\label{Fig13}{ Comparison between numerical analysis and analytical
solution by \citet{Gajo1995b} for three different values of viscous coupling,
$k$: a) solid displacement versus time from 416$\times$10$^{6}$~sec, b) solid
displacement versus time from 56.5$\times$10$^{6}$~sec, c) fluid displacement versus
time from 416$\times$10$^{6}$~sec, and d) fluid displacement versus time from
56.5$\times$10$^{6}$~sec, obtained $1$~cm below the surface.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
In general, the numerical results are in good agreement with analytical
solution for 1D wave propagation in fully saturated, elastic porous media. For
example, numerical analysis demonstrates that during the propagation of the
first wave, the solid and fluid displacement are in phase with each other,
whereas during the propagation of second wave, the displacements are in opposite
phase. In case of high viscous coupling, only one longitudinal wave exists,
since the relative motion between solid and fluid phase is highly constrained.
However, as expected, the rise time of wave fronts obtained from the numerical
analysis is larger than the rise time estimated by the analytical solution, due
to numerical dispersion, observed also by \citet{Gajo1994} and
\citet{Simon1986}. This discrepancy is smaller for high viscous coupling which
physically induces a more dispersive response, in contrast to low
viscous coupling which produces abrupt wavefronts. This dispersive response in case of large viscous coupling is
attributed to the physically based, high viscous damping due to coupling of pore
fluid and porous solid, which can be captured by the $u$$p$$U$ formulation
(see matrices $C_1$, $C_2$, and $C_2$ of Eq.~(\ref{691}), in the appendix).
Furthermore, the rise time of the second wave is even larger; a fact that can be
partly attributed to the diffusive behavior of the second wave in lowfrequency
range. In general, the distortion and the smearing of the wave fronts is a
drawback of all types of numerical solutions, which is linked to the limitation
of the highest frequency allowed by the computational grid and the numerical
damping introduced to the model.
In order to further examine the sensitivity of the numerical solution to
characteristic numerical features, such as a) spatial refinement and b) numerical
damping, an indicative parametric investigation was conducted. For this purpose, a finer spatial discretization of $\Delta h=
0.005\text{~cm}$ was selected, accompanied by a smaller time step equal to
$10^{8}\text{~sec}$. Two different sets of parameters for Newmark integrator
were used: i) $\gamma = 0.5$ and $\beta = 0.25$ corresponding to no numerical
damping and ii) $\gamma = 0.6$ and $\beta = 0.3025$, which introduce numerical
damping. Figure~\ref{NM05} shows
the numerical response versus the analytical solution for two combinations
of spatial and temporal refinement when no numerical damping is introduced in
the model ($\gamma = 0.5$ and $\beta = 0.25$).
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_step_NM05.pdf}
\vspace*{0.3cm}
\caption{\label{NM05}{ Analytical solution by \citet{Gajo1995b} versus numerical
analysis with no numerical damping ($\gamma = 0.5$ and $\beta = 0.25$) for two
different combinations of spatial and temporal discretization:i) $\Delta
h=0.01\text{~cm}$ and $\Delta t = 2 \times 10^{8}\text{~sec}$, ii)
$\Delta
h=0.005\text{~cm}$ and $\Delta t=10^{8}\text{~sec}$, in the case of viscous
coupling $k = 10^{5}\text{~cm$^3$s/g}$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Highfrequency artificial
oscillations develop within the model at the first wave arrival with gradually
decreasing amplitude over time.
%
%
%
Evidently, the oscillations are narrower with equal or smaller amplitude in
case of the finer mesh, while the rise time of abrupt wavefronts diminishes,
better approaching the analytical solution.
Figure~\ref{NM06} shows the comparison of numerical results with analytical
solution for the two selected meshes when numerical damping is used in the
numerical scheme ($\gamma = 0.6$ and $\beta = 0.3025$).
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_step_NM06.pdf}
\vspace*{0.3cm}
\caption{\label{NM06}{ Analytical solution by \citet{Gajo1995b} versus numerical
analysis with no numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$) for two
different combinations of spatial and temporal discretization:i) $\Delta
h=0.01\text{~cm}$ and $\Delta t=2x10^{8}\text{~sec}$, ii) $\Delta
h=0.005\text{~cm}$ and $\Delta t=10^{8}\text{~sec}$, in the case of viscous
coupling $k = 10^{5}\text{~cm$^3$s/g}$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Once more, closer
convergence with analytical solution is achieved in case of finer mesh. It is
interesting to notice that the slight oscillations occurring close to the
wavefronts in case of coarser mesh, have been completely damped out by
numerical damping in case of finer spatial discretization. This fact leads to
the conclusion that numerical damping is, in this case, more effective
in damping out higherfrequency oscillations which develop in case of finer
mesh.
Overall, Figure~\ref{zoom} suggests that the analytical solution is better
approached in terms of the rise time of wavefronts, when no numerical damping is
introduced in the model and finer mesh is used.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_step_zoom.pdf}
\vspace*{0.3cm}
\caption{\label{zoom}{ Analytical solution by \citet{Gajo1995b} versus numerical
analysis for two cases of numerical damping: a) $\gamma = 0.5$ and $\beta =
0.25$ (no numerical damping) and b) $\gamma = 0.6$ and $\beta = 0.3025$ and for
two different combinations of spatial and temporal discretization: i) $\Delta
h=0.01\text{~cm}$ and $\Delta t=2 \times 10^{8}\text{~sec}$, ii) $\Delta
h=0.005\text{~cm}$ and $\Delta t=10^{8}\text{~sec}$, in case of viscous
coupling $k = 10^{5}\text{~cm$^3$s/g}$. Magnified view of the first wavefront.
}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
However, numerical damping is needed to damp out the artificial
oscillations developed due to the discretization of continuum into finite
element model in combination with the range of frequencies imposed by the
excitation.
%
The compromise in choosing the parameters of the numerical features is
directly related to the desired level of accuracy in capturing the details of
response, where both low and high frequency waves are present.
% It should be noted that high
%frequency waves are damped out faster.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Single sine pulse displacement}
\label{1D wave propagation due to single sine pulse displacement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to investigate the efficiency of $u$$p$$U$ finite element under more
realistic excitations with a characteristic predominant frequency, a single sine
pulse displacement was imposed to both solid and fluid phases at the top of a
finite soil column. The characteristic duration of the single sine pulse is
$10^{6}\text{~sec}$, corresponding to a predominant frequency of
$10^6\text{~Hz}$. The soil properties assigned to the model are the same as in
the previous example and are shown in the third column of Table~\ref{Soil Properties_gajo}. The
minimum wavelength, $\lambda_{\min}$, is estimated as $V_{\min}/f_{\max}=
(610\text{~m/s})/(10^6\text{Hz}) = 6.1 \times 10^{4}\text{~m}$. According to
Eq.~(\ref{110a}), the element size is limited to $1/10$ of the minimum
wavelength, so that it is less than $6.1 \times{10^{5}}\text{~m}$.
The selected mesh is comprised of $800$ $u$$p$$U$ brick finite elements
of dimensions $0.005\text{~cm}\times{0.005\text{~cm}}\times{0.005\text{~cm}}$,
creating a soil column $4\text{~cm}$ high. The time step is limited to
$\Delta t \le \Delta
h/V_{\max} = 2.7 \times{10^{8}}\text{~sec}$ according to Eq.~(\ref{111a}). The
chosen temporal discretization consists of time intervals equal to
$\Delta t = 2 \times 10^{8}\text{~sec}$.
%
Numerical models with three different values of
permeability, resulting in three different viscous coupling cases ($k =
10^{8}\text{~cm$^3$s/g}$, $k = 5 \times 10^{8}\text{~cm$^3$s/g}$, and $k =
10^{5}\text{~cm$^3$s/g}$), were analyzed.
Figure~\ref{gajo_sine} shows the comparative results for all three different
values of viscous coupling using Newmark parameters: $\gamma = 0.5$ and $\beta =
0.25$, introducing no numerical damping.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_sine.pdf}
\vspace*{0.3cm}
\caption{\label{gajo_sine}{ Comparison between numerical analysis and analytical
solution by \citet{Gajo1995b} in terms of solid displacement, obtained
$1\text{~m}$ below the surface, for three different values of viscous coupling,
$k$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
%
%
In general, the numerical results converge well with analytical solution,
especially for higher values of viscous coupling where the response is
physically more dispersive. One can observe that the there is a slight
discrepancy between analytical and numerical response in terms of the rise time
at the arrival of the pulse, for cases of lower viscous coupling. This fact emphasizes the
sensitivity of the response to the spatial refinement even under smoother
excitations with a predominant characteristic frequency, particularly for cases
of low viscous coupling. The numerical error can be eliminated with finer mesh,
as shown in the previous case. However, herein, the effect of numerical
dispersion is significantly smaller compared to the example of step excitation and it could be suppressed in case of completely smooth pulse (e.g. Ricker).
Moreover, no need for numerical damping occurred since no high frequency
artificial oscillations developed as in case of step pulse.
%old
%In general, the numerical results
%compare well with analytical solution, especially for higher values of viscous
%coupling where the response is physically more dispersive. One can observe that
%the there is discrepancy between analytical and numerical response in
%the rise time at the arrival of the pulse due to numerical dispersion. This fact
%emphasizes the sensitivity of the response to the spatial refinement even under
%smoother excitations, particularly for cases with low viscous coupling. The
%numerical error can be eliminated with finer mesh, as shown in the previous
%case. There was no need for numerical damping since there were no artificial
%high frequency oscillations as in case of step loading.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Step loading}
\label{One dimensional shock wave propagation with step loading at the surface}
{\bf Analytical Solution by \citet{deBoer1993}} 
%
An analytical solution for onedimensional wave propagation in
fluidsaturated elastic porous media is provided by \citet{deBoer1993}, based on
theory of mixtures \citep{Bowen1980,Bowen1982,Ehlers1993}. The fluidsaturated
porous material is assumed to be a two phase system composed of incompressible
solid and fluid phases. An exact analytical solution is obtained via Laplace
transform technique exhibiting only one independent compressive wave in both the
solid and fluid phases, as a result of the incompressibility constraint.
The problem configuration consists of an infinitely long column, extracted
from the halfspace of a fluidsaturated porous elastic skeleton material. The
motions of both solid and fluid materials are constrained to occur only in the
vertical direction. Loading is applied to the surface boundary, as a function of
time, $\sigma(z=0,t) = f(t)$. Both solid and fluid displacements, as well as
solid skeleton extra stress and pore pressure are obtained with respect to time and
depth \citep{deBoer1993}. The analytical solution has been
developed for i) sinusoidal, ii) step (presented here) and iii) impulsive
loading.
{\bf Numerical Analysis} 
%
Numerical example for the step loading case was used to verify the
implementation of $u$$p$$U$ formulation. Step loading excitation contains
all frequencies, thus requiring a fine spatial discretization. However, this
example does not examine the details of wavefront; it rather focuses on the long
term response describing the consolidation process under dynamic loading.
Considering the above, the numerical grid used for the simulation of the 1D
shock wave propagation consists of 1000 $u$$p$$U$ brick finite elements of
dimensions $1\text{~cm}\times1\text{~cm}\times1\text{~cm}$ creating a soil
column $10$~m high. Since the numerical simulation of a semi infinite soil
column is not possible, a soil column of finite depth of $10$~m was considered
adequate for the current problem configuration. No lateral flow or displacement
are allowed in the model; thus, only the vertical displacement is free. The pore
pressure is constrained at the top surface to be equal to the atmospheric
pressure; thus, modeling drained condition. The base of the model is rigid and impervious.
Table~\ref{Soil Properties_gajo} shows the soil properties of the numerical
model, which are the same with those used for the analytical results presented
in the paper by \citet{deBoer1993}. The only difference is associated with the
elastic modulus, which was selected to be $20$~MN/m$^2$ for the numerical
solution instead of $30$~MN/m$^2$, as noted by \citet{deBoer1993}. The value of Darcy permeability assigned to the soil skeleton is equal to 0.01 m/s. The numerical
analysis indicates that the analytical results presented in the paper correspond
to a soil column with elastic modulus equal to $20$~MN/m$^2$ instead of
$30$~MN/m$^2$. Moreover, it should be mentioned that the solid and fluid (water)
compressibility were given realistic values (see the fourth column of Table~\ref{Soil Properties_gajo}), which practically means that the two constituents can be
assumed fairly incompressible compared to the soilskeleton compressibility.
At the top surface of the soil column, a step loading of $\sigma(z=0,t) =
3\text{~kN/m$^2$}$ is applied to the solid part through four nodal loads at the
top. The nodal loads are developed as:
%
\begin{eqnarray}
F_N(z=0,t)
=
\frac{\sigma(z=0,t)\times{A}}{4}
=
\frac{3\text{~kN/m$^2$}\times0.01\text{~m}\times0.01\text{~m}}{4}
=
7.5\times10^{5}\text{~kN}
\label{1}
\end{eqnarray}
The maximum wave velocity of the soil can be given approximately by equation:
%
\begin{eqnarray}
V_{\max}
=
\sqrt{\frac{M+\alpha\times{Q}}{(1n)\times{\rho_s}+n\times{\rho_f}}}
=
2000\text{~m/s}
\label{vmax}
\end{eqnarray}
%
where $M$ is the constrained modulus of the solid skeleton, $\alpha$ is Biot
coefficient, while $Q$ is the bulk stiffness of a mixture, defined after Eq.
(\ref{GEconservationFlow}) on page \pageref{GEconservationFlow}.
% %
Thus, according to Eq.~(\ref{111a}) the time step should be
$5\times{10^{6}}\text{~sec}$ or less. The selected temporal integration involves 80000
steps of $\Delta t = 5.0\times10^{6}\text{~sec}$. The Newmark parameters used in the analysis are: $\gamma = 0.7$ and $\beta =
0.4$. This sets of parameters introduce considerable numerical damping to the
model which can increase numerical dispersion, but as mentioned before, the
current study does not focuses on the details of the high frequency components
which are quickly damped out with time.
%
Figures~\ref{Fig_disp}, \ref{Fig_stress}, and \ref{Fig_pp} show the
comparison between analytical and numerical results.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Boer_disp.pdf}
\vspace*{0.5cm}
\caption{\label{Fig_disp}{Comparison between numerical analysis and analytical
solution by \citet{deBoer1993} in terms of a) solid displacements versus time at
different depths, b) fluid displacements versus time at different depths, c)
solid displacements versus depth at different times, and d) fluid displacements
versus depth at different times.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{.5cm}
\end{figure}
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Boer_stress.pdf}
\vspace*{0.5cm}
\caption{\label{Fig_stress}{Comparison between numerical analysis and
analytical solution by \citet{deBoer1993} in terms of a) solidskeleton stress
versus time at different depths, and b) solidskeleton stress versus depth at
different times.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{.5cm}
\end{figure}
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Boer_pp.pdf}
\vspace*{0.5cm}
\caption{\label{Fig_pp}{Comparison between numerical analysis and analytical
solution by \citet{deBoer1993}: a) pore fluid pressure versus time at different
depths, b) pore fluid pressure stress versus depth at time $t=0.01$~sec, c) pore
fluid pressure stress versus depth at time $t=0.05$~sec, and d) pore fluid
pressure stress versus depth at time $t=0.1$~sec.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{.5cm}
\end{figure}
%
In general, the numerical
results are in good agreement, with respect to time and depth, with those
obtained by the analytical solution. During the consolidation process, the solid
skeleton settles and the fluid escapes out of the solid skeleton while the load
is transferred from the pore pressure to the solid skeleton stress. The only
discrepancy between analytical and numerical results is found in the pore
pressures which develop an oscillatory trend in case of numerical analysis. The
existence of oscillatory waves has also been observed and discussed by
\citet{Zienkiewicz84}. These artificial oscillations, attributed to the
combination of numerical features and the nature of dynamic step loading, are
damped out over time. Numerical damping would be more efficient in eliminating
the artificial oscillations if finer mesh was used and thus, higherfrequency
oscillations were developed. This conclusion, indicating that the efficiency of
numerical damping increases as the artificial oscillations become narrower, was
shown in section \ref{1DDisp_BC} and will be also shown in section \ref{1DVelBC}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Sinusoidal loading}
\label{1D wave propagation due to sinusoidal loading}
Numerical example for the sinusoidal loading case was also used to verify
implementation of $u$$p$$U$ formulation. The frequency of the excitation is
equal to $75$ rad/s and the load is given by the function
$3\times{(1\text{cos}(75\times{t}))}$. Table~\ref{Soil Properties_gajo} shows
the soil properties of the numerical model, which are the same as in the
previous section. The value of Darcy permeability assigned to the soil skeleton is again equal to 0.01 m/s.
The minimum wave velocity that can be present in the model can be approximately
estimated by equation:
%
\begin{eqnarray}
V_{\min}=\sqrt{\frac{M}{(1n)\times{\rho_s}}}=128\text{~m/s}
\label{vmin}
\end{eqnarray}
%
corresponding to a dry solid skeleton. Thus, the minimum wavelength can be
estimated by Eq.~(\ref{110a}) as $10.7\text{~m}$ leading to a maximum element
size of $1.07\text{~m}$. The numerical grid used for the simulation of the 1D
shock wave propagation consists of 100 $u$$p$$U$ brick finite elements of
dimensions $10\text{~cm} \times 10\text{~cm} \times 10\text{~cm}$ creating a
soil column $10\text{~m}$ high. The maximum timestep required is estimated by
Eq.~(\ref{111a}) equal to $4\times{10^{6}}\text{~sec}$ approximately for the
selected element size. The Newmark parameters were selected to be: $\gamma = 0.6$ and $\beta = 0.3025$, in order to introduce numerical damping. At the top
surface of the soil column, the sinusoidal loading is applied to the solid part
through four nodal loads at the top.
Figures~\ref{Boer_sine} depicts the timehistories of solid and fluid
displacements, as well as pore pressure and effective stress as obtained by both
numerical analysis and analytical solution.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/deBoer_sine.pdf}
\vspace*{0.5cm}
\caption{\label{Boer_sine}{Comparison between numerical analysis and analytical
solution by \citet{deBoer1993} in terms of a) solid displacements, b) fluid
displacements, c) pore pressure and d) solid skeleton stress versus time for
different depths.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{.5cm}
\end{figure}
%
Numerical analysis captures
the response of fully saturated soil under sinusoidal
loading in a satisfactory manner. In particularly, the solid displacement gradually
accumulates downwards
while the fluid is squeezed out of the skeleton. However, there are moments of
partial recovery, when the load on the surface decreases, which renders the pore
pressure values negative close to surface, as the fluid is absorbed inwards
(suction). The effective stress seems to be more sensitive to loading close to the
surface following the shape of the excitation, whereas the effect tends to
diminish for greater depths due to solidfluid interaction.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Step velocity excitation}
%\label{One dimensional shock wave propagation with a step velocity boundary condition}
\label{1DVelBC}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf SemiAnalytical Solution by \citet{Hiremath88}} 
%
\citet{Hiremath88} present a semianalytical solution of Biot's equations of
motion for onedimensional wave propagation in a fluidsaturated linear elastic
isotropic medium (soil) using Laplace transformation followed by numerical
inversion. This study is considered to be an extension of the exact transient
solution presented by \citet{Garg74} for two limiting cases of infinitely small
and infinitely large viscous coupling. In both cases, a soil column of finite
dimension (50~cm) subjected to velocity boundary conditions at the surface was
analyzed, allowing for reflection of waves at the boundaries.
One of the most important observations, made by both \citet{Garg74}
and \citet{Hiremath88}, is that in case of strong viscous coupling (high drag),
the material behaves as a single continuum with internal dissipation and the two
wave fronts tend to become a single one.
In particular, \citet{Hiremath88} examined two cases allowing for small
and large viscous coupling. Moreover, two different types of excitations were
applied at the boundary surface in terms of solid and fluid velocity. In the
first case, a unit step boundary condition was applied at the top surface for
both solid and fluid phases (to be studied below). In the second case, the fluid
velocity applied at the boundary is slightly different from the
applied solid
velocity, and it is increasing gradually to unity over a short time scale. The
results obtained
from the numerical inversion allowed for six reflections of the fast
compressional wave of first kind and two reflections of the secondary slow
longitudinal wave.
{\bf Numerical Analysis} 
%
Models for two extreme values of viscous coupling: a) large
($k=0.148\times10^{8}\text{~cm$^3$s/g}$) and b) low
($k=0.148\times10^{2}\text{~cm$^3$s/g}$) are developed and simulated. At the
top surface of the soil column, a step velocity of $1.0\times10^{2}\text{~m/s}$
is applied both to the solid and the fluid phase. Only the vertical
translational degrees of freedom are free so that no lateral flow or
displacement are allowed. The base of the model is rigid and impervious. A
step excitation introduces all frequencies in the model,
requiring a very small element size (theoretically, based on Eq.~(\ref{110a}),
the element size should really be infinitely small, in order to accommodate all
introduced frequencies). The
selected finite element mesh consists of 100 $u$$p$$U$ brick finite elements
of dimensions $0.005\text{~m}\times0.005\text{~m}\times0.005\text{~m}$ creating
a soil column 50~cm high, as suggested by \citet{Hiremath88}. Table~\ref{Soil Properties_gajo} shows the soil properties used for this model. The solid
skeleton is so rigid that the wave of the first kind propagates with a velocity
close to 3500m/s according to equation Eq.~(\ref{vmax}). The required time step
is estimated by Eq.~(\ref{111a}) and is equal to $\Delta t =
1.4\times10^{6}\text{~sec}$. The selected temporal integration involves 1972
steps of $\Delta t = 5.0\times10^{7}\text{~sec}$. The Newmark set of parameters
was selected as $\gamma = 0.6$ and $\beta = 0.3025$, introducing numerical
damping in the model.
Figures~\ref{solid_10_high} and \ref{solid_30_high} show the comparison
between numerical and analytical results for both extreme cases of viscous
coupling.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Hiremath_high.pdf}
\vspace*{0.5cm}
\caption{
\label{solid_10_high}
Comparison between numerical analysis and analytical solution
\citet{Hiremath88} for the case of strong viscous coupling
($k=0.148\times10^{8}\text{~cm$^3$s/g}$) in terms of time histories of: a) solid velocity at
$10$~cm below the surface, b) solid velocity at $30$~cm below the surface, c) fluid
velocity at $10$~cm below the surface and, d) fluid velocity at $30$cm below the
surface.}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Hiremath_low.pdf}
\vspace*{0.5cm}
\caption{\label{solid_30_high}{Comparison between numerical analysis and
analytical solution \citet{Hiremath88} for the case of weak viscous coupling
($k=0.148\times10^{2}\text{~cm$^3$s/g}$): a) solid
velocity at $10$~cm below the surface, b) solid velocity at $30$~cm below the
surface, c) fluid velocity at $10$~cm below the surface and d) fluid velocity at
$30$cm below the surface, versus time.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Overall, the finite element solution reproduces correctly the trends
of wave propagation in both limiting cases of viscous coupling. In particular,
numerical analysis demonstrates well that for the case of strong viscous
coupling, solid and fluid velocities are in phase with each other.
%It should be
%noted that the numerical analysis develops oscillations at the instantaneous
%changes in velocity due to reflection of wave fronts.
In order to improve the accuracy of the numerical solution, a finer spatial and
temporal discretization was used, with $\Delta h = 0.001\text{~m}$ and $\Delta t
= 2\times{10^{7}}\text{~sec}$, while Newmark parameters remained the same.
%
Figure~\ref{hiremath_NM06} depicts the analytical solution versus the numerical
results obtained from both meshes.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/hiremath_NM06.pdf}
\vspace*{0.5cm}
\caption{\label{hiremath_NM06}{Analytical solution \citet{Hiremath88} versus
numerical analysis, with two different combinations of spatial and temporal
discretization and numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$), for
the case of weak viscous coupling ($k=0.148\times10^{2}\text{~cm$^3$s/g}$): a)
solid velocity at 10 cm below the surface and b) fluid velocity at 10 cm below
the surface, with time.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Evidently, the finer mesh achieves a greater
accuracy, suppresses oscillations, and gives a more similar response to the analytical
solution, as expected.
%
It is interesting to notice that when numerical damping is used, all artificial
oscillations are suppressed in case of finer mesh, while they persist in case of
coarser mesh, as illustrated in Figure~\ref{hiremath_fine}.
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/hiremath_fine.pdf}
\caption{\label{hiremath_fine}{Analytical solution \citet{Hiremath88} versus
numerical analysis with no numerical damping ($\gamma = 0.6$ and $\beta =
0.3025$), in terms of fluid velocity time history at 10 cm below the surface,
for the case of weak viscous coupling ($k=0.148\times10^{2}\text{~cm$^3$s/g}$):
a) $\Delta h = 0.001\text{~m}$ and $\Delta t = 2 \times 10^{7}\text{~sec}$ and
b) $\Delta h = 0.005\text{~m}$ and $\Delta t = 5 \times 10^{7}\text{~sec}$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
A solution verification suite for modeling and
simulation of fully coupled behavior of saturated porous media, using upU formulation,
have been presented in this paper.
%
Analytical solutions for both static and
dynamic examples have been included in this suite.
%
%
Detailed comparison between numerical and analytical solutions shows close
matching, while common numerical discretization effects are captured and
controlled. In particular, artificiallyintroduced (by discretization)
higherfrequency effects have been observed, as in most numerical simulations.
These numerical dispersion effects have been investigated focusing on the
influence of numerical features such as spatial and temporal refinement, as well
as numerical damping on the obtained response. It is illustrated that the finer
meshes better approach the analytical solutions, as expected; however, such
finer meshes produce even higher frequency (artificial) oscillations. It is also
shown that such oscillations produced by finer meshes can be more efficiently
damped out numerically than those produced by coarser discretization.
%
Presented solution verification suite can be used by any
modeling and simulation effort that deals with the response of fully saturated porous
media. In this particular case, presented solution verification results provide
evidence that the fully coupled models for
saturated porous media, from the UCD Computational Geomechanics Group numerical
libraries, are solved correctly for both static and
dynamic cases.
It is important to observe that any numerical method develop to model and
simulation behavior of solids/structures and fluids needs to be verified (and
validated if quality experimental data is available). Only after such
verification (and validation) process has been successfully completed, a method
can be claimed to be complete. The main aim of this paper, with presented
verification of the upU formulation, is to contribute to such approach to
numerical modeling and simulations.
%Presented solution verification results provide evidence that
%the fully coupled model for saturated porous media was solved correctly by UCD Computational Geomechanics
%Group's
%implementation of 3D $u$$p$$U$ finite elements \citep{Jeremic2008}. The developed solution verification procedure can be used for any finite element or finite difference code oriented to simulate the response of fully saturated porous media.
%The next step would be
%developing a validation suite, providing evidence that the correct model (or set of
%models) is (are) solved. For validation purposes, high quality experimental results are
%required; thus, finding the most appropriate ones, is our next goal.
%\clearpage
\begin{appendix}
%\section{APPENDIX}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Finite Element Formulation}
\label{Finite Element Formulation}
Standard finite element discretization \citep{FEMIV1,FEMIV2}, using shape
functions, is used 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 (in this case same) shape functions
for solid displacement,
pore pressure, and fluid displacement, respectively. Each node of the
($u$$p$$U$) element thus features seven degrees of freedoms in three 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 the
$u$$p$$U$ formulation with compressible fluid allows that without volumetric
locking.
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 Eqs.~(\ref{34}), (\ref{35}), and
(\ref{36}):
% %%
%%
%
%
%
%
\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 V \\
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}
%%
where the left hand side components are: %%
% %%%
%%%
\begin{eqnarray}
(M_s)_{KijL}= \int_{V} N_K^u (1n) \rho_s \delta_{ij} N_L^u \ud V
\;\;\;\; &\mbox{;}& \;\;\;\;
(M_f)_{KijL}= \int_{V} N_K^U n \rho_f \delta_{ij} N_L^U \ud V \nonumber\\
(C_1)_{KijL}= \int_{V} N_K^u n^2 k_{ij}^{1} N_L^u \ud v
\;\;\;\; &\mbox{;}& \;\;\;\;
(C_2)_{KijL}= \int_{V} N_K^u n^2 k_{ij}^{1} N_L^U \ud V \nonumber\\
(C_3)_{KijL}= \int_{V} N_K^U n^2 k_{ij}^{1} N_L^U \ud V
\;\;\;\; &\mbox{;}& \;\;\;\;
(G_1)_{KiM} = \int_{V} N_{K,i}^u (\alphan) N_M^p \ud V \nonumber\\
(G_2)_{KiM} = \int_{V} n N_{K,i}^U N_M^p \ud V
\;\;\;\; &\mbox{;}& \;\;\;\;
P_{NM} = \int_{V} N_N^p \frac{1}{Q} N_M^p \ud V
\label{691}
\end{eqnarray}
%%%
%
%
while the right hand side are:
% %%%
\begin{eqnarray}
(\bar{f}_s)_{Ki} &=& \int_{A_t} N_K^u \sigma^{'}_{ij} n_j \ud A

\int_{A_p} N_K^u (\alphan) p n_i \ud A
+
\int_{V} N_K^u (1n) \rho_s b_i \ud V \nonumber\\
(\bar{f}_f)_{Ki} &=&  \int_{A_p} N_K^U n p n_i \ud A
+
\int_{V} N_K^U n \rho_f b_i \ud V
\label{692}
\end{eqnarray}
%
where $V$ is the volume, and $A_t$ and $A_p$ are the domain boundaries with
traction and the pore fluid pressure defined, respectively.
It is very important to note that the velocity proportional viscous
damping was introduced directly through the damping tensor 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, velocity proportional energy dissipation due to interaction of pore fluid
and the solid (soil) skeleton. It is also emphasized
that presented
formulation and implementation do not (need to) use Rayleigh damping.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Time Integration}
\label{Time_Integration}
In order to develop integration of dynamic finite element equation in the time
domain, Eq.~(\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.
Eq.~(\ref{71}) represents the general nonlinear form \citep{local87} for which
the usual tangent stiffness $\bf K$ obtained from:
\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_{V} N_{K,m}^u E^{ep}_{imjn} N_{L,n}^u \ud V
\label{upU_Kep}
\end{eqnarray}
The above set of residual (nonlinear) dynamic equations is solved using the
Newmark procedure \citep{Newmark1959}. This time integration method has two parameters,
$\beta$ and $\gamma$, and is described by the following two equations:
%%%%
%{\sc \Large I checked this, same with Lecture Notes}
\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 current time step $n$ to the next time step $n+1$.
The method is generally an implicit method, except when both $\beta$ and $\gamma$ are zero.
%
%{\sc \Large I checked this, same with Lecture Notes}
%
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$, it is second order accurate \citep{local86}.
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \section{Program Implementation}
% \label{PI}
%
% %{\sc \Large NEEDS TO BE REWRITTEN}
%
% Implementation of the described algorithms and procedures was performed using a
% number of numerical libraries, algorithms, material models, and finite
% elements available through one of open source licenses.
% %
% %
% The finite element domain is managed using Modified OpenSees Services (MOSS)
% library
% \citep{Jeremic2004d, Jeremic2007d, Jeremic2008a, McKenna97}.
% %
% On a lower functional level, a set of Template3Dep numerical libraries
% \citep{Jeremic2000f}
% are used for constitutive level modeling, while nDarray numerical libraries
% \citep{Jeremic97d} are used to handle vector, matrix
% and tensor manipulations, and FEMtools element libraries from UCD
% CompGeoMech toolset \citep{Jeremic2004d} are used to supply other
% necessary libraries and components.
% %
% %
% %
% The solution of system of equations is provided by UMFPACK solver
% \citep{Davis97b, Davis97, Davis1999, Davis2004b, Davis2004}.
% %
\end{appendix}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\include{Appendix}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\bibliography{refmech,Mahdi,zhao}
%\bibliography{./refmech,./refcomp,./Mahdi,./zhao,./ref_Giota}
\bibliography{refmech,refcomp,Mahdi,zhao,ref_Giota}
\bibliographystyle{abbrvnat}
%\bibliographystyle{agsm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%TABLES
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
% \clearpage
% \begin{table}[!htb]
% \begin{center}
% \caption{ \label{Soil Properties for quasistatic problems}
% Soil properties used in numerical analysis of quasistatic problems.}
% %\begin{spacing}{1.4}
% %\setlength{\extrarowheight}{3pt}
% \begin{tabular}{lccc}
% Parameter & Symbol [Units] & 1D consolidation & Line injection \\
% \hline
% gravity acceleration & $g$ [m/s$^2$] & $9.81$ & $9.81$ \\
% soil Young's Modulus & $E$ [kN/m$^2$] & $10^4$ & $1.2\times 10^6$ \\
% soil Poisson's ratio & $v$ & $0.25$ & $0.2$ \\
% solid density & $\rho_s$ [ton/m$^3$] & $2.65$ & $2.7$ \\
% fluid density & $\rho_f$ [ton/m$^3$] & $1.0$ & $2.7$ \\
% solid bulk modulus & $K_s$ [kN/m$^2$] & $37\times 10^6$ & $3.6\times 10^7$ \\
% fluid bulk modulus & $K_f$ [kN/m$^2$] & $2.2\times 10^6$ & $1.0\times 10^{17}$ \\
% porosity & $n$ & 0.46 & 0.4 \\
% Biot coefficient & $\alpha$ & $1.0$ & $1.0$\\
% Darcy permeability & $k_D$ [m/s] & $10^3$ & $3.6\times 10^{6}$ \\
% %\hline
% soil drained bulk modulus & $K$ [kN/m$^2$] & $6.67\times 10^3$ & $6.7\times 10^5$ \\
% soil undrained bulk modulus & $K_u$ [kN/m$^2$] & $4\times 10^6$ & $6.0\times 10^7$ \\
% % soil shear modulus & $\mu$ & $4\times 10^3$~kN/m$^2$ \\
% fluid diffusivity coefficient & $c_f$ [m$^2$/s] & $1.2$ & $0.5$ \\
% % characteristic consolidation time & $\tau$ & $83.3$~s \\
% \end{tabular}
% %\end{spacing}
% \end{center}
% \end{table}
%
%\begin{table}[!htb]
%\begin{center}
%\caption{\label{table8}{Soil properties used to study the line injection problem.}}
%\begin{spacing}{1.4}
%\setlength{\extrarowheight}{3pt}
%\begin{tabular}{lcc}
% Parameter & Symbol & Value \\
%\hline
%Young's modulus & E & $1.2\times 10^6$ kN/m$^2$ \\
%Poisson ratio & $\nu$ & $0.2$ \\
%solid density & $\rho_s$ & $2.7$ ton/m$^3$ \\
%fluid density & $\rho_f$ & $1.0$ ton/m$^3$ \\
%solid 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$ \\
%porosity & $n$ & 0.4 \\
%Darcy permeability & $k_{D}$ & $3.6\times 10^{6}$ m/s \\
%bulk modulus & $K$ & $6.7\times 10^5$ kN/m$^2$ \\
%undrained bulk modulus & $K_u$ & $6.0\times 10^7$ kN/m$^2$ \\
%fluid diffusivity coefficient & $c_f$ & $0.4973$ m$^2$/s \\
%\end{tabular}
%\end{spacing}
%\end{center}
%\end{table}
% \begin{table}[!htb]
% \begin{center}
% \caption{ \label{Soil Properties_gajo}
% Soil Properties for 1D wave propagation due
% to various boundary conditions and loading.}
% %\begin{spacing}{1.4}
% \begin{tabular}{lcccc}
% Parameter & Symbol [Units] & Step/Sine pulse displacement & Step/Sinusoidal loading & Step velocity \\
% \hline
% gravity acceleration & $g$ [m/s$^2$] & $9.81$ & $9.81$ & $9.81$ \\
% soil Young's Modulus & $E$ [kN/m$^2$] & $1.2\times 10^6$ & $20\times 10^3$ & $23.2\times 10^6$ \\
% soil Poisson's ratio & $v$ & $0.3$ & $0.2$ & $0.17$ \\
% solid density & $\rho_s$ [ton/m$^3$] & $2.7$ & $2.0$ & $2.66$ \\
% water density & $\rho_f$ [ton/m$^3$] & $1.0$ & $1.0$ & $1.0$ \\
% solid bulk modulus & $K_s$ [kN/m$^2$] & $36\times 10^6$ & $36\times 10^6$ & $36\times 10^6$ \\
% fluid bulk modulus & $K_f$ [kN/m$^2$] & $2.2\times 10^6$ & $2.2\times 10^6$ & $2.2\times 10^6$ \\
% porosity & $n$ & $0.4$ & $0.33$ & $0.18$ \\
% Biot coefficient & $\alpha$ & $1.0$ & $1.0$ & $0.677$ \\
% \end{tabular}
% %\end{spacing}
% \end{center}
% \end{table}
%
%\begin{table}[!htb]
%\begin{center}
%\caption{ \label{Soil Properties_boer}
%Soil Properties for 1D wave propagation due
%to step loading.}
%\begin{spacing}{1.4}
%\begin{tabular}{lcc}
% Parameter & Symbol & Value \\
%\hline
% gravity acceleration & $g$ & 9.81~m/s$^2$ \\
% soil matrix Young's Modulus & $E$ & $20\times 10^3$~kN/m$^2$ \\
% soil matrix Poisson's ratio & $v$ & 0.2 \\
% solid particle density & $\rho_s$ & $2.0\times 10^3$~kg/m$^3$ \\
% water density & $\rho_f$ & $1.0\times 10^3$~kg/m$^3$ \\
% solid particle bulk modulus & $K_s$ & $36.0\times 10^6$~kN/m$^2$ \\
% fluid bulk modulus & $K_f$ & $2.177\times 10^6$~kN/m$^2$ \\
% porosity & $n$ & 0.33 \\
% Darcy permeability & $k_D$ & 0.01~m/s \\
%\end{tabular}
%\end{spacing}
%\end{center}
%\end{table}
%\begin{table}[!htb]
%\begin{center}
%\caption{ \label{Soil Properties_Hiremath}
%Soil Properties for 1D wave propagation due to step velocity boundary condition.}
%\begin{spacing}{1.4}
%\begin{tabular}{lcc}
% Parameter & Symbol & Value \\
%\hline
% gravity acceleration & $g$ & 9.81~m/s$^2$ \\
% soil matrix Young's Modulus & $E$ & $23.21\times 10^6$~kN/m$^2$ \\
% soil matrix Poisson's ratio & $v$ & 0.171 \\
% solid particle density & $\rho_s$ & $2.66\times 10^3$~kg/m$^3$ \\
% fluid density & $\rho_f$ & $1.0\times 10^3$~kg/m$^3$ \\
% solid particle bulk modulus & $K_s$ & $36.0\times 10^6$~kN/m$^2$ \\
% fluid bulk modulus & $K_f$ & $2.2\times 10^6$~kN/m$^2$ \\
% porosity & $n$ & 0.18 \\
% Biot coefficient & $\alpha$ & $0.6772$ \\
%\end{tabular}
%\end{spacing}
%\end{center}
%\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%FIGURES
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
% \clearpage
%
% %
% \begin{figure}[!htb]
% \begin{center}
% {\includegraphics[width=9.0cm]{FIGS/RoleVV_NEW_knowledge.pdf}}
% \caption{\label{RoleVV} Schematic representation of the role of Verification and
% Validation on numerical modeling (inspired from \protect{\citet{Oberkampf2002,Roy2011}} and
% \protect{\citet{Oden2010a}}).}
% \end{center}
% \end{figure}
% \begin{figure}[!htb]
% \begin{LARGE}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=0.50\textwidth]{FIGS/cons_mod1.pdf}
% \caption{\label{num_model}{(a) The geometry of the soil layer of thickness $h$ under uniform constant vertical pressure $\varpi$ = 400~kPa applied at the surface, and (b) the finite element mesh with the appropriate boundary conditions and loading at top.}}
% \end{center}
% \end{sffamily}
% \end{LARGE}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/cons_pp.pdf}
% \caption{\label{dissipation}{Comparison between numerical analysis and
% analytical solution by \citet{Coussy2004}: a) normalized pore pressure at
% various depths with respect to real time, $t$, and b) distribution of normalized
% pore pressure with normalized depth, $z/h$, for different time factors, $T_v$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=0.5\textwidth]{FIGS/cons_disp.pdf}
% \caption{\label{settle}{(a) Comparison between numerical analysis and analytical
% solution by \citet{Coussy2004} in terms of time histories: a) solidskeleton
% displacements, and b) fluid displacements, for various depths.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
%
% \begin{figure}[!htb]
% \begin{center}
% \includegraphics[width=8cm]{FIGS/Injection.pdf}
% %\includegraphics[width=8cm]{FIGS/InjectionPressure_new_ps1.pdf}
% %\hfill
% %\includegraphics[width=7.8cm]{FIGS/InjectionDisplacement_new_ps1.pdf}
% \caption{\label{fig33} {Comparison between numerical analysis and analytical solution by \citet{Coussy2004}: a) radial displacement, and b) pore pressure versus time at different distances, R, from the center of the model.}}
% \end{center}
% \end{figure}
% \begin{figure}[!htb]
% \begin{center}
% \includegraphics[width=8cm]{FIGS/Inject1.pdf}
% \caption{\label{fig32} {Mesh generated for line injection problem}}
% \end{center}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=0.50\textwidth]{FIGS/Gajo_mod.pdf}
% \caption{\label{Fig1001}{a) A representative soil column subjected to a step vertical displacement equal to $1.0\times10^{3}\text{~cm}$ at the surface, and b) the finite element mesh and the applied boundary conditions used for numerical modeling.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Gajo_comp.pdf}
% \caption{\label{Fig13}{ Comparison between numerical analysis and analytical
% solution by \citet{Gajo1995b} for three different values of viscous coupling,
% $k$: a) solid displacement versus time from 416$\times$10$^{6}$~sec, b) solid
% displacement versus time from 56.5$\times$10$^{6}$~sec, c) fluid displacement versus
% time from 416$\times$10$^{6}$~sec, and d) fluid displacement versus time from
% 56.5$\times$10$^{6}$~sec, obtained $1$~cm below the surface.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_step_NM05.pdf}
% \caption{\label{NM05}{ Analytical solution by \citet{Gajo1995b} versus numerical
% analysis with no numerical damping ($\gamma = 0.5$ and $\beta = 0.25$) for two
% different combinations of spatial and temporal discretization:i) $\Delta
% h=0.01\text{~cm}$ and $\Delta t = 2 \times 10^{8}\text{~sec}$, ii)
% $\Delta
% h=0.005\text{~cm}$ and $\Delta t=10^{8}\text{~sec}$, in the case of viscous
% coupling $k = 10^{5}\text{~cm$^3$s/g}$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_step_NM06.pdf}
% \caption{\label{NM06}{ Analytical solution by \citet{Gajo1995b} versus numerical
% analysis with no numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$) for two
% different combinations of spatial and temporal discretization:i) $\Delta
% h=0.01\text{~cm}$ and $\Delta t=2x10^{8}\text{~sec}$, ii) $\Delta
% h=0.005\text{~cm}$ and $\Delta t=10^{8}\text{~sec}$, in the case of viscous
% coupling $k = 10^{5}\text{~cm$^3$s/g}$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_step_zoom.pdf}
% \caption{\label{zoom}{ Analytical solution by \citet{Gajo1995b} versus numerical
% analysis for two cases of numerical damping: a) $\gamma = 0.5$ and $\beta =
% 0.25$ (no numerical damping) and b) $\gamma = 0.6$ and $\beta = 0.3025$ and for
% two different combinations of spatial and temporal discretization: i) $\Delta
% h=0.01\text{~cm}$ and $\Delta t=2 \times 10^{8}\text{~sec}$, ii) $\Delta
% h=0.005\text{~cm}$ and $\Delta t=10^{8}\text{~sec}$, in case of viscous
% coupling $k = 10^{5}\text{~cm$^3$s/g}$. Magnified view of the first wavefront.
% }}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_sine.pdf}
% \caption{\label{gajo_sine}{ Comparison between numerical analysis and analytical
% solution by \citet{Gajo1995b} in terms of solid displacement, obtained
% $1\text{~m}$ below the surface, for three different values of viscous coupling,
% $k$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Boer_disp.pdf}
% \caption{\label{Fig_disp}{Comparison between numerical analysis and analytical solution by \citet{deBoer1993} in terms of a) solid displacements versus time at different depths, b) fluid displacements versus time at different depths, c) solid displacements versus depth at different times, and d) fluid displacements versus depth at different times.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{.5cm}
% \end{figure}
%
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Boer_stress.pdf}
% \caption{\label{Fig_stress}{Comparison between numerical analysis and analytical solution by \citet{deBoer1993} in terms of a) solidskeleton stress versus time at different depths, and b) solidskeleton stress versus depth at different times.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{.5cm}
% \end{figure}
%
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Boer_pp.pdf}
% \caption{\label{Fig_pp}{Comparison between numerical analysis and analytical solution by \citet{deBoer1993}: a) pore fluid pressure versus time at different depths, b) pore fluid pressure stress versus depth at time $t=0.01$~sec, c) pore fluid pressure stress versus depth at time $t=0.05$~sec, and d) pore fluid pressure stress versus depth at time $t=0.1$~sec.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{.5cm}
% \end{figure}
%
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/deBoer_sine.pdf}
% \caption{\label{Boer_sine}{Comparison between numerical analysis and analytical
% solution by \citet{deBoer1993} in terms of a) solid displacements, b) fluid
% displacements, c) pore pressure and d) solid skeleton stress versus time for
% different depths.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{.5cm}
% \end{figure}
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Hiremath_high.pdf}
% \caption{
% \label{solid_10_high}
% Comparison between numerical analysis and analytical solution
% \citet{Hiremath88} for the case of strong viscous coupling
% ($k=0.148\times10^{8}\text{~cm$^3$s/g}$) in terms of time histories of: a) solid velocity at
% $10$~cm below the surface, b) solid velocity at $30$~cm below the surface, c) fluid
% velocity at $10$~cm below the surface and, d) fluid velocity at $30$cm below the
% surface.}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
%
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Hiremath_low.pdf}
% \caption{\label{solid_30_high}{Comparison between numerical analysis and
% analytical solution \citet{Hiremath88} for the case of weak viscous coupling
% ($k=0.148\times10^{2}\text{~cm$^3$s/g}$): a) solid
% velocity at $10$~cm below the surface, b) solid velocity at $30$~cm below the
% surface, c) fluid velocity at $10$~cm below the surface and d) fluid velocity at
% $30$cm below the surface, versus time.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
%
% \clearpage
%
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/hiremath_NM06.pdf}
% \caption{\label{hiremath_NM06}{Analytical solution \citet{Hiremath88} versus
% numerical analysis, with two different combinations of spatial and temporal
% discretization and numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$), for
% the case of weak viscous coupling ($k=0.148\times10^{2}\text{~cm$^3$s/g}$): a)
% solid velocity at 10 cm below the surface and b) fluid velocity at 10 cm below
% the surface, with time.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
%
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/hiremath_fine.pdf}
% \caption{\label{hiremath_fine}{Analytical solution \citet{Hiremath88} versus
% numerical analysis with no numerical damping ($\gamma = 0.6$ and $\beta =
% 0.3025$), in terms of fluid velocity time history at 10 cm below the surface,
% for the case of weak viscous coupling ($k=0.148\times10^{2}\text{~cm$^3$s/g}$):
% a) $\Delta h = 0.001\text{~m}$ and $\Delta t = 2 \times 10^{7}\text{~sec}$ and
% b) $\Delta h = 0.005\text{~m}$ and $\Delta t = 5 \times 10^{7}\text{~sec}$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
%
%
\end{document}
\bye