\documentclass[10pt]{article}
\usepackage{amsmath}
\newcommand\bm[1] {{\bf #1}}
\newcommand\Div {{\rm div\ }}
\newcommand\grad {{\rm grad\ }}
\newcommand\Dp {\delta p}
\newcommand\dphi {\delta \varphi}
\newcommand\dpsi {\delta \psi}
\newcommand\du {\delta \bm u}
\newcommand\odt {O(\Delta t^2)}
\newcommand\odd {O(\Delta t^3)}
\newcommand\nal {\par\medskip\noindent}

\title{Short note on the Biot soil consolidation}
\author{Jos van Kan}
\begin{document}
\maketitle

\begin{abstract}
We introduce briefly a computational scheme for the time dependent Biot
soil consolidation model. The scheme is equivalent to solving 6 Laplacian
equations per time step. We argue that the scheme is $O(\Delta t^2)$ in the pore
pressure $p$ and the {\em dilatation} $\Div \bm u$ in which $\bm u$ is the soil
displacement vector.
\end{abstract}

\section{Introduction}
The Biot soil consolidation model is given by
\begin{align}
-\mu\tilde\Delta \bm u -(\lambda + \mu)\grad\Div \bm u +\grad p &=0\\
\frac{\partial{\ }}{\partial t}(\Div \bm u) - \frac{\kappa}{\eta} \Delta p &=
f(\bm x, t)
\end{align}
with boundary conditions $p=0$ on $\Gamma$ and $\bm u=0$ on $\Gamma$. Initial
conditions are given for the dilatation $\Div  \bm u$. To facilitate the
calculations we introduce a special symbol $\varphi$ for the dilatation $\Div
\bm u$ which transforms the system to:
\begin{align}
 \frac{\partial \varphi}{\partial t} -\frac{\kappa}{\eta} \Delta p &= f(\bm x,
t)\\
-\mu\tilde\Delta \bm u - (\lambda + \mu)\grad \varphi +\grad p &=0\\
\Div \bm u &= \varphi
\end{align}
This also introduces an extra boundary condition 
\begin{equation}
\varphi = \frac{\partial u_n }{\partial n} \quad {\rm on\ } \Gamma.
\end{equation} 


\section {Numerical model}
We assume that we have a natural staggered space discretization. ($\varphi$- and
$p$-points must coincide). We also assume that we have the quantities $p$, $\bm
u$ and $\varphi$ on time level $n$  with sufficient accuracy and the
straightforward way to calculate an update $\Dp$, $\du$ and $\dphi$ would be the
solution of the following system of simultaneous linear equations:
\begin {align}
\dphi - \frac12 \Delta t \frac{\kappa}{\eta} \Delta\ \Dp &= \Delta t
[\frac{\kappa}{\eta} \Delta\ p^n + f^{n+1/2}]\\
\label{eq:elastic}-\mu\tilde\Delta\ \du - (\lambda + \mu) \grad \dphi + \grad
\Dp &=0\\
\dphi - \Div \du = 0 
\end {align}
The accuracy of such an update would be $O(\Delta t^3)$ at best (that is the
error of the mid point integration) and for that to happen the error in $p^n$
should be no worse than $O(\Delta t^2)$. Fortunately these two characteristics
are compatible.
\nal
This simultaneous system is a bit unwieldy and we therefore try a
simplification that won't lose us the overall accuracy. First of all we note,
by taking the divergence of the second equation that
\begin{equation}\label{eq:phi-p}
-(\lambda+2\mu) \Delta \dphi + \Delta \Dp =0.
\end{equation}
Therefore we may write the {\em first} equation of the system as:
\begin{equation}\label{eq:time}
\dphi - \frac12 \Delta t \frac{\kappa (\lambda+2\mu)}{\eta} \Delta\ \dphi =
\Delta t[\frac{\kappa}{\eta} \Delta\ p^n + f^{n+1/2}]
\end{equation}
with boundary condition
\begin{equation}
\dphi=\frac{\partial \du_n }{\partial n} \quad {\rm on\ } \Gamma.
\end{equation}
Unfortunately, we do not know $\du_n$ yet, but we can be sure, that it is of
order $O(\Delta t)$. So if we take $\dphi^*=0$ on $\Gamma$ instead, we have an
error of $O(\Delta t)$ in the boundary conditions, but only an error $O(\Delta
t^2)$  in the solution $\dphi^*$ of system \eqref{eq:time} on the interior
points, since those boundary condtions are multiplied by $\Delta t$ in the
equations. 
\nal
We define $\dpsi= \Dp-(\lambda+2\mu)\dphi$. From equation \eqref{eq:phi-p} we
now could calulate $\dpsi^*$ with an accuracy
of $\odt$, provided we had boundary values to $\dphi^*$ with that accuracy.
We know  $\dphi^*$ in the interior with an accuracy of $\odt$, hence we
can extrapolate to the boundary without loss of accuracy. From this solution
$\dpsi*$ we can obtain $\Dp^*$ with $\odt$ accuracy.
\nal
We plug $\Dp^*$ and $\dphi^*$ into
equation \eqref{eq:elastic} to obtain $\du^*$ with $\odt$ accuracy. Note that
this system is uncoupled in the $\bm u$ components, hence is equivalent to two
Laplacian equations. This yields $\du^*$ with $\odt$ accuracy. Using this in the
boundary condition of equation \eqref{eq:time} we obtain an $\odd$ accurate 
solution $\dphi$. From this we extrapolate the boundary contions
for $\dphi$ once more to obtain $\odd$ accurate boundary conditions and an $\odd$
accurate solution $\dpsi$ to \eqref {eq:phi-p}. From $\dpsi$ we get an $\odd$
accurate $\Dp$.
\nal
This suffices to obtain the pressure $p$ and the dilatation $\varphi$ with an
global $\odt$ accuracy for the cost of 6 Laplacian like equations per time step.
To get the displacement $\bm u$ to the same accuracy would require an
additional two Laplacian equations.
\nal
There may be a stability problem. 


\end{document}
