Finite element method with the total stress variable for Biot’s consolidation model

Biot’s fundamental equations for the soil consolidation process that describes the gradual adaptation of the soil to a load variation have been established in [1]. The mechanism of consolidation phenomenon described using a linear isotropic model is identical with the process of squeezing water out of an elastic porous medium in many cases. This solid-fluid coupling was extended to general anisotropy case in [2]. Such poroelasticity models have a lot applications in many areas including geomechanics [3], medicine [4], biomechanics [5], reservoir engineering [6]. The Biot’s consolidation models have also been used to combine transvascular and interstitial fluid movement with the mechanics of soft tissue in [7], which can be applied to improve drug delivery in solid tumors. Existence, uniqueness, and regularity theory were developed in [8] for poroelasticity and quasi-static problem in thermoelasticity. In [9], experiments with finite element method for the model in [7] have demonstrated the effects of fluid flow on the spatio-temporal patterns of soft-tissue elastic strain under a variety of physiological condition, which emphasized simulations relevant to a quasistatic elasticity imaging for the characterization of fluid flow in solid tumors.

Biot’s consolidation model has been considered by many researchers using finite element methods. In [10], a variational principle and the finite element method for a model with applications to a nonhomogeneous, anisotropic soil were developed. The fully discretization with backward Euler time discrete finite element method has been carried out and the existence and uniqueness were proved in [11]. Moreover, the simplest Taylor-Hood finite elements were employed. The stability and error estimates for the semi-discrete and fully-discrete finite element approximations were derived in [12] and [13]. Decay functions and post processing technique also were employed to improve the pore pressure accuracy. With the displacement and pore pressure fields as unknowns, the short and long time behavior of spatially discrete finite element solutions have been studied in [14]. Asymptotic error estimates have been derived for both stable and unstable combinations of the finite element spaces. For the coupled problem, a least squares mixed finite element method was presented and analyzed in [15] with the unknowns fluid displacement, stress tensor, flux, and pressure. In [16], coupling of mixed and continuous Galerkin finite element methods for pressure and displacements have been formulated deriving the convergence error estimates in time continuous setting. Methods for coupling mixed and discontinuous Galerkin have been presented in [17]. The error estimates for a fully-discrete stabilized discontinuous Galerkin method were obtained with the unknowns pressure and displacement in [18]. In [19], a discretization method in irregular domains with general grids for discontinuous full tensor permeabilities was developed. A new mixed finite element method for Biot’s consolidation problem in four variables was proposed in [20] and later, a three field mixed finite element which was free of pressure oscillations and Poisson locking has been proposed [21]. The priori error estimates that were robust for material parameters were provided in [22] with a four-field mixed method formulation. A three field finite element formulation with nonconforming finite element space for the displacements was considered in [23]. Based on the parameter dependent norms, the parameter-robust stability was established in [24], and parameter robust inf-sup stability and strong mass conservation were derived for three field mixed discontinuous Galerkin discretizations. A stabilized finite element method with equal order elements for the unknowns pressure and displacement was proposed in [25] to reduce the effects of non-physical oscillations. Combining the mixed method with symmetric interior penalty discontinuous Galerkin method obtained a H(div) conforming finite element method in [26]. The method achieved strong mass conservation.

Moreover, in [27], the total stress (or the soil pressure) expressed as a combination of the divergence of the velocity and pressure has been introduced coupling the solid and fluid robustly for Biot’s consolidation problem with the unknown displacement, pressure, and volumetric stress. Using a Fredholm argument for a static model, error estimates were derived independently of the Lam$\acute{e}$ constants for both continuous and discrete formulations. A three field formulation of Biot’s model with the total stress variable has also been proposed in [28], which developed a parameter-robust block diagonal preconditioner for the associated discrete systems. Then, in [29], a priori error estimates for semi-discrete scheme has been presented by introducing the total pressure variable for quasi-static multiple-network poroelasticity equations. Our goal in this paper is to emphasize the time dependence of field variables in error estimates, so we choose to consider the fully-discrete scheme. Thus, using the total stress as a new variable for the three field formulation, we give the error estimates for semi-discrete and fully-discrete with backward Euler time discretization schemes.

Let $\Omega$ be an open bounded polygonal or polyhedral domain in $\mathbb{R}^{d},~{}d=2,3$. Biot’s consolidation model is described as follows: Find the displacement vector $\mathbf{u}$ and the fluid pressure $p$ in

where $\Gamma_{D}$ and $\Gamma_{N}$ are disjoint closed subsets with $\Gamma_{D}\cup\Gamma_{N}=\partial\Omega$ and the Dirichlet boundary $|\Gamma_{D}|\neq 0$. Here, $\varepsilon(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})$ is the strain tensor expressed in terms of symmetrized gradient of displacements, $\kappa\in(0,\infty)$ is the permeability of the porous solid and parameters $\mu\in(0,\infty)$, $\lambda\in(0,\infty)$ are the elastic Lam$\acute{e}$ constants. The right hand side term $\mathbf{f}$ in (1) represents the density of the applied body forces, and the source term $g$ in (2) represents a forced fluid extraction or injection. The outward unit normal vector on $\partial\Omega$ is denoted by $\mathbf{n}$.

Next, we introduce the coupling between the solid and fluid using $q=-\lambda\nabla\cdot\mathbf{u}+p$, where $q$ is the total stress, $\mathbf{u}$ is the displacement and $p$ is the fluid pressure [27, 28]. This can be rewritten as

For $x\in\Omega$, the initial conditions are given by

Hence, the initial condition for the total stress is $q(x,0)=-\lambda\nabla\cdot\mathbf{\varphi}+\phi$. Denote the spaces $[H^{1}_{0,D}(\Omega)]^{d}=\{\mathbf{v}\in[H^{1}(\Omega)]^{d},~{}\mathbf{v}=\mathbf{0}~{}\mbox{on}~{}\Gamma_{D}\}$ and $H^{1}_{0,D}(\Omega)=\{v\in H^{1}(\Omega),~{}v=0~{}\mbox{on}~{}\Gamma_{D}\}$.

Multiplying (1), (6) and (2) by test functions, integrating by parts and applying boundary conditions yield the following weak formulation: Find $\mathbf{u}\in[H^{1}_{0,D}(\Omega)]^{d}$, $q\in L^{2}(\Omega)$, $p\in H^{1}_{0,D}(\Omega)$ such that

In Section 2, we present the semi-discrete and fully-discrete finite element formulation for system (8) with the unknown displacement, total stress and pressure and prove uniqueness of the solution for each of these schemes. Section 3 presents the derivation and analysis of the error estimates for both the semi-discrete and fully-discrete schemes. In Section 4 we present computational experiments on benchmark problems that validate the theoretical convergence rates with respect to mesh size $h$ and time step $\tau$.

In the paper, we denote the arbitrary constants by $\epsilon_{i}\geq 0$, where $i$ is positive integer and $C$ is a constant which is independent of time step $\tau$ and mesh size $h$. Let ${P}_{k}$ be the space of polynomials of degree less than or equal to $k$ in all variables. Moreover, let $\|\cdot\|$ be the norm in $L^{2}(\Omega)$ space and $\|\cdot\|_{k}$ be the norm in $H^{k}(\Omega)$ space. Denote the space-time space by $L^{k}(0,\bar{T};V)$ for a Banach space $V$ (see details in [34]).

Let $\mathcal{T}_{h}$ be a regular and quasi-uniform triangulation of domain $\Omega$ into triangular or tetrahedron elements [30, 31]. For each element $T\in\mathcal{T}_{h}$, $h_{T}$ is its diameter and $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ is the mesh size of triangulation $\mathcal{T}_{h}$. We consider the following finite element spaces on $\mathcal{T}_{h}$,

where $k\geq 2,~{}l\geq 0.$ In order to describe the initial conditions of discretization schemes, we define the following two projection operators. Let us define the Stokes projection $Q_{h}^{\mathbf{u}}:[H^{1}_{0,D}(\Omega)]^{d}\rightarrow\mathbf{U}_{h}$ and $Q_{h}^{q}:L^{2}(\Omega)\rightarrow Z_{h}$ by

Also the elliptic projection $Q_{h}^{p}:H^{1}_{0,D}(\Omega)\rightarrow P_{h}$ is defined with the following properties,

Hence, given a suitable approximation of initial conditions $\mathbf{u}_{h}(0)=Q_{h}^{\mathbf{u}}\varphi$, $q_{h}(0)=Q_{h}^{q}(-\lambda\nabla\cdot\varphi+\phi)$ and $p_{h}(0)=Q_{h}^{p}\phi$, the semi-discrete scheme corresponding to a three field formulation (8) is for all $t\in[0,\bar{T}]$, to seek $\mathbf{u}_{h}(t)\in\mathbf{U}_{h}$, $q_{h}(t)\in Z_{h}$, $p_{h}(t)\in P_{h}$ such that

To obtain a fully-discrete formulation, we denote time step by $\tau$, and $t^{n}=n\tau$, where $n$ is non-negative integer. Thus, given a suitable approximation of initial conditions $\mathbf{u}^{0}=Q_{h}^{\mathbf{u}}\varphi$, $q^{0}=Q_{h}^{q}(-\lambda\nabla\cdot\varphi+\phi)$ and $p^{0}=Q_{h}^{p}\phi$, the fully-discrete scheme with backward Euler time discretization corresponding to the three field formulation (8) is to find $\mathbf{u}^{n}\in\mathbf{U}_{h}$, $q^{n}\in Z_{h}$, $p^{n}\in P_{h}$ such that

where $\bar{\partial}_{t}{\mathbf{u}^{n}}:=\displaystyle\frac{\mathbf{u}^{n}-\mathbf{u}^{n-1}}{\tau}$ and $f^{n}:=f(x,t^{n}),~{}x\in\Omega$.

Then, we present an inequality, which will be useful in the uniqueness and convergence analysis.

Next we use Lemma 1 to prove the uniqueness for semi-discrete (11) and fully-discrete (12) schemes.

In order to derive the error estimates, we first show the inf-sup condition for space pair $([H^{1}_{0,D}(\Omega)]^{d},L^{2}(\Omega))$ with $|\Gamma_{N}|>0$. Denote $b(\mathbf{v},w_{q}):=(w_{q},\nabla\cdot\mathbf{v})$. For each $w_{q}\in L^{2}(\Omega)$, there exists $\mathbf{v}\in[H^{1}_{0,D}(\Omega)]^{d}$ satisfying $w_{q}=\nabla\cdot\mathbf{v}$ and $|\mathbf{v}|_{1}\leq C\|w_{q}\|$. Thus, we can deduce that for $\forall w_{q}\in L^{2}(\Omega)$[35, 27]

Since not all the discrete finite element spaces meet the inf-sup condition, we assume that the space pair $(\mathbf{U}_{h},Z_{h})$ satisfy the inf-sup condition, i.e. there exists a positive constant such that for $\forall w_{q_{h}}\in Z_{h}$