Parameter-robust stability of classical three-field formulation of Biot's consolidation model

TL;DR

This paper establishes parameter-robust stability for a classical three-field formulation of Biot's consolidation model, enabling the development of effective preconditioners and stable discretizations that ensure mass conservation.

Contribution

It introduces parameter-dependent norms that prove full stability of the model uniformly across all parameters, including the Lamé parameter, facilitating robust numerical methods.

Findings

Proves uniform inf-sup stability with respect to all model parameters.

Constructs a block diagonal preconditioner based on the stability analysis.

Designs stable discretizations with optimal error estimates and mass conservation.

Abstract

This paper is devoted to the stability analysis of a classical three-field formulation of Biot's consolidation model where the unknown variables are the displacements, fluid flux (Darcy velocity), and pore pressure. Specific parameter-dependent norms provide the key in establishing the full parameter-robust inf-sup stability of the continuous problem. Therefore, stability results presented here are uniform not only with respect to the Lam\'e parameter $\lambda$, but also with respect to all the other model parameters. This allows for the construction of a uniform block diagonal preconditioner within the framework of operator preconditioning. Stable discretizations that meet the required conditions for full robustness and guarantee mass conservation, both locally and globally, are discussed and corresponding optimal error estimates proven.

Equations

$$-\mathrm{div}\mathbf{\sigma }+c_{up}\nabla p$$

$$\mathbf{v}$$

$$c_{pu}\mathrm{div}\dot{\mathbf{u}}-\mathrm{div}\mathbf{v}+c_{pp}\dot{p}$$

$$\mathbf{ϵ}(\mathbf{u})$$

$$\mathbf{\sigma }$$

$$\lambda :=\frac{\nu E}{(1+\nu )(1-2\nu )},\mu :=\frac{E}{2(1+\nu )}$$

$$p(\mathbf{x},t)$$

$$\mathbf{v}(\mathbf{x},t)\cdot \mathbf{n}(\mathbf{x})$$

$$\mathbf{u}(\mathbf{x},t)$$

$$\mathbf{\sigma }(\mathbf{x},t)\mathbf{n}(\mathbf{x})$$

$$p(\mathbf{x},0)$$

$$\mathbf{u}(\mathbf{x},0)$$

$$\mathcal{A}\left[\begin{array}[]{c}{\bm{u}}^{k}\\ {\bm{v}}^{k}\\ p^{k}\end{array}\right]=\left[\begin{array}[]{c}{\bm{f}}^{k}\\ {\bm{0}}\\ \tilde{g}^{k}\end{array}\right]$$

$$\mathcal{A}:=\left[\begin{array}[]{ccc}-2\mu\operatorname{div}{\bm{\epsilon}}-\lambda\nabla\operatorname{div}&0&c_{up}\nabla\\ 0&\tau K^{-1}&\tau\nabla\\ -c_{pu}\operatorname{div}&-\tau\operatorname{div}&-c_{pp}I\end{array}\right]$$

$${\mathbf{u}}^k=\mathbf{u}(\mathbf{x},t_k)\in \mathbf{U}$$

$${\mathbf{v}}^k=\mathbf{v}(\mathbf{x},t_k)\in \mathbf{V}$$

$$p^k=p(\mathbf{x},t_k)\in P$$

$$2\mu \to 1,\lambda /2\mu \to \lambda ,c_{up}/2\mu \to c_{up},\mathbf{f}/2\mu \to \mathbf{f},\tau /2\mu \to \tau ,c_{pp}/2\mu \to c_{pp},g/2\mu \to g.$$

$$-\mathrm{div}\mathbf{ϵ}(\mathbf{u})-\lambda \nabla \mathrm{div}\mathbf{u}+c_{up}\nabla p$$

$$\tau K^{-1}\mathbf{v}+\tau \nabla p$$

$$-c_{pu}\mathrm{div}\mathbf{u}-\tau \mathrm{div}\mathbf{v}-c_{pp}p$$

$$-\mathrm{div}\mathbf{ϵ}(\tilde{\mathbf{u}})-\lambda \nabla \mathrm{div}\tilde{\mathbf{u}}+\nabla \tilde{p}$$

$${\alpha }^2{\tau }^{-1}{K}^{-1}\tilde{\mathbf{v}}+\nabla \tilde{p}$$

$$-\mathrm{div}\tilde{\mathbf{u}}-\mathrm{div}\tilde{\mathbf{v}}-c_{pp}{\alpha }^{-2}\tilde{p}$$

$$R_p^{-1}={\alpha }^2{\tau }^{-1}{K}^{-1},{\alpha }_p=c_{pp}{\alpha }^{-2},$$

$$\lambda \ge 1,R_p^{-1}\ge 0,{\alpha }_p\ge 0.$$

$$-\mathrm{div}\mathbf{ϵ}(\mathbf{u})-\lambda \nabla \mathrm{div}\mathbf{u}+\nabla p$$

$$R_p^{-1}\mathbf{v}+\nabla p$$

$$-\mathrm{div}\mathbf{u}-\mathrm{div}\mathbf{v}-{\alpha }_{p}p$$

$$A\left[\begin{array}[]{c}{\bm{u}}\\ {\bm{v}}\\ p\end{array}\right]=\left[\begin{array}[]{c}\bm{f}\\ {\bm{0}}\\ \tilde{g}\end{array}\right]$$

$$A:=\left[\begin{array}[]{ccc}-\operatorname{div}{\bm{\epsilon}}-\lambda\nabla\operatorname{div}&0&\nabla\\ 0&R_{p}^{-1}I&\nabla\\ -\operatorname{div}&-\operatorname{div}&-\alpha_{p}I\end{array}\right].$$

$$(\mathbf{ϵ}(\mathbf{u}),\mathbf{ϵ}(\mathbf{w}))+\lambda (\mathrm{div}\mathbf{u},\mathrm{div}\mathbf{w})-(p,\mathrm{div}\mathbf{w})$$

$$(R_p^{-1}\mathbf{v},\mathbf{z})-(p,\mathrm{div}\mathbf{z})$$

$$-(\mathrm{div}\mathbf{u},q)-(\mathrm{div}\mathbf{v},q)-{\alpha }_p(p,q)$$

$$(\mathbf{u},\mathbf{w}{)}_{\mathbf{U}}$$

$$(\mathbf{v},\mathbf{z}{)}_{\mathbf{V}}$$

$$(p,q{)}_P$$

$$\underset{q\in P}{\inf }\underset{\mathbf{v}\in \mathbf{V}}{\sup }\frac{(\mathrm{div}\mathbf{v},q)}{\Vert \mathbf{v}{\Vert }_{\mathrm{div}}\Vert q\Vert }\ge {\beta }_d.$$

$$\underset{q\in P}{\inf }\underset{\mathbf{u}\in \mathbf{U}}{\sup }\frac{(\mathrm{div}\mathbf{u},q)}{\Vert \mathbf{u}{\Vert }_1\Vert q\Vert }\ge {\beta }_s.$$

$$\underset{q\in P}{\inf }\underset{(\mathbf{u},\mathbf{v})\in \mathbf{U}\times \mathbf{V}}{\sup }\frac{(\mathrm{div}\mathbf{u},q)+(\mathrm{div}\mathbf{v},q)}{(\Vert \mathbf{u}{\Vert }_{\mathbf{U}}+\Vert \mathbf{v}{\Vert }_{\mathbf{V}})\Vert q{\Vert }_P}\ge {\beta }_{dv}>0.$$

Full text

1. Introduction: Biot’s consolidation model

Poroelastic models describe mechanical deformation and fluid flow in porous media. They have a wide range of applications in medicine, biophysics and geosciences such as the computation of intracranial pressure, trabecular bone stiffness under different loading conditions, reservoir simulation, waste repository performance, ${\rm CO}_{2}$ sequestration, consolidation of soil under surface loads, subsidence due to fluid withdrawal and many others, see, e.g. [1, 2, 3, 4].

A classical and widely used model has been introduced by Biot [5, 6] and is based on the following assumptions:

• (i)

the porous medium is saturated by fluid and the temperature is constant,

• (ii)

the fluid in the porous medium is (nearly) incompressible,

• (iii)

the solid skeleton (matrix) is formed by an elastic material and deformations and strains are relatively small,

• (iv)

the fluid flow is driven by Darcy’s law (laminar flow).

For homogeneous isotropic linear elastic porous media, the Biot model in an open domain $\Omega\subset I\!\!R^{d}$, $d=2,3$, comprises the following system of partial differential equations (PDEs):

[TABLE]

Here $\lambda$ and $\mu$ denote the Lamé parameters which are defined by

[TABLE]

in terms of the modulus of elasticity (Young’s modulus) $E$ and the Poisson ratio $\nu\in[0,1/2)$.

The constant $c_{up}=c_{pu}=\alpha$ coupling the pore pressure $p$ and the displacement variable ${\bm{u}}$ is the Biot-Willis constant, $K$ is the hydraulic conductivity, given by the quotient between the permeability of the porous medium $\kappa$ and the viscosity of the fluid $\eta$; $\bm{I}$ denotes the identity tensor and ${\bm{\sigma}}$ and ${\bm{\epsilon}}$ are the effective stress and strain tensors, respectively, which are related to each other via the constitutive equation (1e); The strain tensor ${\bm{\epsilon}}({\bm{u}})$ is given by the symmetric part of the gradient of the displacement field as defined in the compatibility condition (1d). The time derivatives of ${\bm{u}}$ and $p$ in the continuity equation (1c) are denoted by $\dot{{\bm{u}}}$ and $\dot{p}$. Finally, ${\bm{v}}$ denotes the fluid flux, sometimes also called percolation velocity of the fluid, which is assumed to be proportional to the (negative) pressure gradient as expressed by Darcy’s law (1b). The right hand side $\bm{f}$ in the equilibrium equation (1a) represents the density of the applied body forces and the source term $g$ in (1c) a forced fluid extraction or injection.

The system (1) is completed by proper boundary and initial conditions, e.g.,

[TABLE]

where $\Gamma_{p,D}\cap\Gamma_{p,N}=\emptyset$, $\overline{\Gamma}_{p,D}\cup\overline{\Gamma}_{p,N}=\Gamma=\partial{\Omega}$ and $\Gamma_{\bm{u},D}\cap\Gamma_{\bm{u},N}=\emptyset$, $\overline{\Gamma}_{\bm{u},D}\cup\overline{\Gamma}_{\bm{u},N}=\Gamma$; Initial conditions at the time $t=0$ to complement the boundary conditions (2) are given by

[TABLE]

Making use of the constitutive equation (1e) to eliminate the stress variable from the system (1) results in the classical three-field formulation of the Biot model.

A common way to solve the time-dependent problem numerically, is then to discretize it in time and solve a static problem in each time step. Using the backward Euler method for time discretization, in this way one obtains a three-by-three block system of time-step equations

[TABLE]

where

[TABLE]

for the unknown time-step functions

[TABLE]

and the right hand side time-step functions ${\bm{f}}^{k}={\bm{f}}(\bm{x},t_{k})$ and $\tilde{g}^{k}=-\tau g(\bm{x},t_{k})-c_{pu}\operatorname{div}({\bm{u}}^{k-1})-c_{pp}p^{k-1}$ at any given time moment $t=t_{k}=t_{k-1}+\tau$. As in the remainder of this paper we will consider the static problem (4)–(5), for convenience, we will drop the superscript for the time-step functions, that is, denote ${\bm{u}}^{k}$, ${\bm{v}}^{k}$ and $p^{k}$ by ${\bm{u}}$, ${\bm{v}}$ and $p$, respectively.

Following the standard notation, $L^{2}(\Omega)$ denotes the space of square Lebesgue integrable functions equipped with the standard $L^{2}$ norm $\|\cdot\|$ and $H^{1}(\Omega)^{d}$ the space of vector-valued $H^{1}$-functions equipped with the norm $\|\cdot\|_{1}$ defined by $\|\bm{u}\|^{2}_{1}:=\|\bm{u}\|^{2}+\|\nabla\bm{u}\|^{2}$. Moreover, $H(\operatorname{div};\Omega):=\{\bm{v}\in L^{2}(\Omega)^{d}:\operatorname{div}\bm{v}\in L^{2}(\Omega)\}$ where the standard Sobolev norm $\|\cdot\|_{\operatorname{div}}$ is defined by $\|\bm{v}\|^{2}_{\operatorname{div}}:=\|\bm{v}\|^{2}+\|\operatorname{div}\bm{v}\|^{2}$.

Frequently we will consider the case in which $\Gamma_{{\bm{u}},D}=\Gamma_{p,N}=\Gamma$ and ${{\bm{u}}}_{D}={\bm{0}}$, $q_{N}=0$ in which we write ${\bm{U}}=H^{1}_{0}(\Omega)^{d}$ and ${\bm{V}}=H_{0}(\operatorname{div},\Omega)$. In order to determine the solution for the pressure variable $p$ uniquely one can set $P=L^{2}_{0}(\Omega):=\{p\in L^{2}(\Omega):\int_{\Omega}p\,d\bm{x}=0\}$.

In many applications the variations of the model parameters are quite large. In geophysical applications, the permeability typically varies in the range from $10^{-9}$ to $10^{-21}m^{2}$ whereas Young’s modulus is typically in the order of GPa and the Poisson ratio in the range $0.1-0.3$, see [3, 7, 8]. Soft tissue of the central nervous system on the other hand has a permeability of about $10^{-14}$ to $10^{-16}m^{2}$ whereas Young’s modulus is typically in the order of kPa and the Poisson ratio in the range $0.3$ to almost $0.5$, see [1, 2]. For that reason it is important that not only the formulation of the problem but also the numerical methods for its solution are stable over the whole range of values of the parameters in the model.

The stability of the time discretization has been studied in [9] and will not be addressed here. Instead we will focus on the issue of inf-sup stable finite element discretizations of the static problem (4)–(5). It is a well known fact that the LBB condition, see [10, 11], plays the crucial role in the well-posedness analysis of the continuous problem and its discrete counterparts arising from mixed finite element discretizations. It is also the key tool in deriving a priori error estimates. Inf-sup stability for the Darcy problem as well as for the Stokes and linear elasticity problems is well understood and various stable mixed discretizations of either of these systems of PDEs have been proposed over the years, see, e.g. [12] and the references therein.

Biot’s model of poroelasticity combines theses equations and the parameter-robust stability of its three-field formulation becomes more delicate as we will see in the next sections. Alternative formulations that can be proven to be stable independently of the model parameters (in certain norms) include a two-field formulation for the displacements and pore pressure, see [7, 13], and a new three-field formulation that–besides the displacements–introduces two pressure unknowns, one for the fluid pressure and one for the total pressure defined as a weighted sum of fluid and solid pressure [7].

Compared to the new three-field formulation presented in [7], the classic three-field formulation of Biot’s consolidation model keeps Darcy’s law in order to guarantee fluid mass conservation. Recently, a four-field formulation has been proposed in which the stress tensor is kept as a variable in the system, see [14], and the error analysis there is robust with respect to $\lambda$, but not uniform with respect to the other parameters such as $\tau$ and $K$.

Nonconforming finite elements have been shown to be beneficial with regard to reducing pressure oscillations in computations based on the classical three-field formulation, see [15]. The lowest approximation order, consisting of Crouzeix-Raviart finite elements for the displacements, lowest-order Raviart-Thomas elements for the Darcy velocity, and a piecewise constant approximation of the pressure unknown, in combination with a mass-lumping technique for the Raviart-Thomas elements results in a computationally efficient method. However, the norms defined in this work, and in many others, do not allow to establish the full parameter-robust stability for which we aim.

In the present paper, we establish this full parameter-robust stability for the classic three-field formulation of Biot’s consolidation model. Crucial in our analysis is the definition of proper norms for which we prove that the constants in the related inf-sup conditions do not depend on any of the model parameters. Further, we propose a discretization that preserves fluid mass conservation at a discrete level. We also prove the full parameter-robust stability of the discretized problem and of course the related optimal error estimates. The remainder of the paper is organized as follows.

In Section 2, we briefly revisit non-uniform stability results, and make some useful observations which motivate the subsequent analysis. In Section 3 we introduce the parameter-dependent norms based on which we establish the parameter-robust stability of the weak formulation of the continuous problem (4)–(5). Section 4 analyzes mixed finite element discretizations that provide discrete parameter-robust inf-sup stability and full fluid mass conservation. Applying the theory of operator preconditioning, see [16], the results from Sections 3 and 4 imply the uniformity (parameter-robustness) of the (canonical) norm-equivalent block-diagonal preconditioners. In Section 5 we use our findings to derive robust optimal a priori error estimates. Finally, Section 6 gives some concluding remarks.

Throughout this paper, the hidden constants in $\lesssim,~{}\gtrsim$ and $\eqsim$ are independent of the parameters $\mu,\lambda,c_{up},\tau,K,c_{pu},c_{pp}$ and the mesh size $h$. Hence the hidden constants in $\lesssim,~{}\gtrsim$ and $\eqsim$ are independent of $\lambda,R_{p},\alpha_{p}$ and the mesh size $h$.

2. A revisit of non-uniform stability results

We begin our stability analysis with recasting the equations (4)–(5). We first divide all the parameters in the model (4)–(5) by $2\mu$ herewith eliminating the parameter $\mu$. That is, we make the substitutions

[TABLE]

Herewith, equation (4) becomes

[TABLE]

Now let $\tilde{\bm{u}}=c_{pu}\bm{u},\tilde{\bm{v}}=\tau\bm{v},\tilde{p}=c_{pu}^{2}p,\tilde{\bm{f}}=c_{pu}\bm{f}$ and divide equation (7b) by $\tau$ to get

[TABLE]

where we have also used that $c_{up}=c_{pu}=\alpha$. For convenience we denote

[TABLE]

and further assume that the range of the parameters is

[TABLE]

These assumptions are very general and reasonable.

In the subsequent, in order to simplify the notation, we skip the “tilde” symbol, i.e., we make the substitutions $\bm{\tilde{u}}\rightarrow\bm{u},\bm{\tilde{v}}\rightarrow\bm{v},\bm{\tilde{f}}\rightarrow\bm{f},\tilde{p}\rightarrow p$, and finally write the system, without loss of generality, in the form

[TABLE]

or, in short notation,

[TABLE]

where

[TABLE]

The weak formulation of (9a)–(9c) reads: Find $(\bm{u},\bm{v},p)\in\bm{U}\times\bm{V}\times P$, such that for any $(\bm{w},\bm{z},q)\in\bm{U}\times\bm{V}\times P$

[TABLE]

Motivated by the work [17], let us first consider the Hilbert spaces ${\bm{U}}=H^{1}_{0}(\Omega)^{d},{\bm{V}}=H_{0}(\operatorname{div},\Omega),$ and $P=L^{2}_{0}(\Omega)$ with the natural norms defined by

[TABLE]

Before we study the Biot’s equations, we recall the following well known results, see, e.g. [11, 12].

Lemma 1**.**

There exists a positive constant $\beta_{v}>0$, such that

[TABLE]

Lemma 2**.**

There exists a positive constant $\beta_{s}>0$, such that

[TABLE]

Let us now turn to the stability of the formulation (12a)–(12c). We consider two different cases. In the first case, for some nearly incompressible materials, we have that the Lamé parameter $\lambda$ tends to infinity. If we assume that $\lambda\gg 1$, $R_{p}^{-1}\eqsim 1$ and $0\leq\alpha_{p}\lesssim 1$ then, defining the norms according to (13)–(15), the boundedness of both $({\bm{\epsilon}}(\bm{u}),{\bm{\epsilon}}(\bm{w}))+\lambda(\operatorname{div}\bm{u},\operatorname{div}\bm{w})+(R_{p}^{-1}\bm{v},\bm{v})$ and $(\operatorname{div}\bm{u},q)+(\operatorname{div}\bm{v},q)$ are obvious. Further, for $0\leq\alpha_{p}\lesssim 1$, we obtain the boundedness of $\alpha_{p}(p,q)$. Moreover, defining the norms by (13)–(15), using Lemma 1, and choosing $(\bm{u},\bm{v})=(\bm{0},\bm{v})\in\bm{U}\times\bm{V}$, we obtain the inf-sup condition

[TABLE]

Finally, the coercivity of $({\bm{\epsilon}}(\bm{u}),{\bm{\epsilon}}(\bm{u}))+\lambda(\operatorname{div}\bm{u},\operatorname{div}\bm{u})+(R_{p}^{-1}\bm{v},\bm{v})$ on the kernel set

[TABLE]

can also be verified. In fact, since $(\bm{u},\bm{v})\in\bm{Z}$ means $\operatorname{div}\bm{v}=-\operatorname{div}\bm{u}$, it follows that

[TABLE]

where the last inequality comes from the assumption $\lambda\gg 1$, $R_{p}^{-1}\eqsim 1$.

On the other hand, in many practical applications one has $\lambda\eqsim 1$, $c_{up}\eqsim 1$, $c_{pp}=c_{up}/{\lambda}$, and $K~{}\hbox{or}~{}\tau$ tending to zero. Then, since in this case we have $R_{p}^{-1}\gg 1$ and $0\leq\alpha_{p}\lesssim 1$, defining the norms according to (13)–(15), the boundedness of both $({\bm{\epsilon}}(\bm{u}),{\bm{\epsilon}}(\bm{w}))+\lambda(\operatorname{div}\bm{u},\operatorname{div}\bm{w})+(R_{p}^{-1}\bm{v},\bm{v})$ and $(\operatorname{div}\bm{u},q)+(\operatorname{div}\bm{v},q)$ are again obvious. Further, the assumption $0\leq\alpha_{p}\lesssim 1$ implies the boundedness of $\alpha_{p}(p,q)$.

Next, using the definition of the norms (13)–(15), Lemma 2, and choosing $(\bm{u},\bm{v})=(\bm{u},\bm{0})\in\bm{U}\times\bm{V}$, we obtain the inf-sup condition

[TABLE]

We see, however, that in this case the coercivity of $({\bm{\epsilon}}(\bm{u}),{\bm{\epsilon}}(\bm{u}))+\lambda(\operatorname{div}\bm{u},\operatorname{div}\bm{u})+(R_{p}^{-1}\bm{v},\bm{v})$ can not be valid any more on the kernel set $\bm{Z}$, where $\bm{Z}$ is defined by (19). In fact, since $(\bm{u},\bm{v})\in\bm{Z}$ means $\operatorname{div}\bm{v}=-\operatorname{div}\bm{u}$, it follows that for any $M>0$, there exits $(\bm{u},\bm{v})$, where, e.g., $\operatorname{div}\bm{v}\neq 0$, and $R_{p}^{-1}$ is large enough, such that

[TABLE]

where the second inequality comes from $\lambda\lesssim 1$.

Using the norms defined in (13)–(14), the estimate (22) implies that for any $M>0$, there exists $(\bm{u},\bm{v})\in\bm{Z}$ (and $R_{p}^{-1}$ large enough) such that

[TABLE]

Therefore the system (12a)–(12c) is not uniformly stable with respect to the parameter $R_{p}^{-1}$ under the norms (13)–(15).

From this observation we conclude that we have to define proper norms (as we do below in (26a)–(26c)) in order to establish the coercivity of $({\bm{\epsilon}}(\bm{u}),{\bm{\epsilon}}(\bm{u}))+\lambda(\operatorname{div}\bm{u},\operatorname{div}\bm{u})+(R_{p}^{-1}\bm{v},\bm{v})$ on $\bm{Z}$ in both above cases.

3. Parameter-robust stability of the model

In this section, we first define proper parameter-dependent norms for the spaces $\bm{U},\bm{V}$ and $P$ based on which we establish then the parameter-robust stability of the Biot’s model (12a)–(12c) for parameters in the ranges

[TABLE]

Let us denote

[TABLE]

and consider the Hilbert spaces ${\bm{U}}=H^{1}_{0}(\Omega)^{d},{\bm{V}}=H_{0}(\operatorname{div},\Omega),P=L^{2}_{0}(\Omega)$ with parameter-dependent norms $\|\cdot\|_{{\bm{U}}}$, $\|\cdot\|_{{\bm{V}}}$, $\|\cdot\|_{P}$ induced by the inner products

[TABLE]

The above norms are the key to establish the parameter-robust stability of the model.

Directly related to problem (12a)–(12c) we introduce the bilinear form

[TABLE]

In view of the definition of the norms (26a)–(26c), the boundedness of $A((\bm{u},\bm{v},p),(\bm{w},\bm{z},q))$ is obvious. We come to our first main result.

Theorem 1**.**

There exists a constant $\beta>0$ independent of the parameters $\lambda,~{}R_{p}^{-1},~{}\alpha_{p}$, such that

[TABLE]

Proof.

Case I:

[TABLE]

For any $(\bm{u},\bm{v},p)\in\bm{U}\times\bm{V}\times P$, by Lemma 2, there exists

[TABLE]

Choose

[TABLE]

where $\delta$ is a positive constant that will be determined later.

First we verify the boundedness of $(\bm{w},\bm{z},q)$ by $(\bm{u},\bm{v},p)$.

By (30), and noting that $\gamma^{-1}\leq\rho=\lambda$, we have

[TABLE]

Next, since $1\leq\lambda$, we get the boundedness of $\bm{w}$, i.e.,

[TABLE]

It is obvious that $\|\bm{z}\|_{\bm{V}}=\delta\|\bm{v}\|_{\bm{V}}$. We still need to bound $q$. Using (31) and

[TABLE]

we obtain

[TABLE]

Collecting the estimates above we get

[TABLE]

Next we show the coercivity of $A((\bm{u},\bm{v},p),(\bm{w},\bm{z},q))$. Using the definition of $(\bm{w},\bm{z},q)$ and (30), we find

[TABLE]

Applying Cauchy’s inequality and using (30), we therefore obtain

[TABLE]

Now, letting $\epsilon_{1}=\frac{\beta_{s}^{2}}{2},\epsilon_{2}=\epsilon_{3}=\epsilon_{4}=\frac{1}{2}$ and noting that $\rho^{-1}\leq\gamma,\lambda\geq\gamma^{-1}>0$ and $\lambda\geq 1$, we further conclude that

[TABLE]

Next, letting $\delta:=\max\{\beta_{s}^{-2}+\frac{1}{2},2+\frac{1}{2}\}$ and noting that $\alpha_{p}\leq\gamma,\lambda\geq 1$, we arrive at the following coervicity estimate

[TABLE]

Case II:

[TABLE]

For any $(\bm{u},\bm{v},p)\in\bm{U}\times\bm{V}\times P$, by Lemma 1, there exists

[TABLE]

Choose

[TABLE]

where $\delta$ is a constant which we will specify later.

Again we verify the boundedness of $(\bm{w},\bm{z},q)$ by $(\bm{u},\bm{v},p)$ first. It is obvious that $\|\bm{w}\|_{\bm{U}}=\delta\|\bm{u}\|_{\bm{U}}$.

Moreover, by (40) and noting that $\gamma^{-1}\leq\rho=R_{p}^{-1}$, we have

[TABLE]

Hence, we get the boundedness of $\bm{z}$, i.e.,

[TABLE]

The boundedness of $q$ follows as in (35).

Next we verify the coercivity of $A((\bm{u},\bm{v},p),(\bm{w},\bm{z},q))$ in Case II.

Using the definition of $(\bm{w},\bm{z},q)$ and (39), we find

[TABLE]

Applying Cauchy’s inequality and using (39), we get

[TABLE]

Now, letting $\epsilon_{1}=\displaystyle\beta_{d}^{2},\epsilon_{2}=\epsilon_{4}=\displaystyle\frac{1}{2}$ and noting that $\rho^{-1}\leq\gamma,\lambda\geq\rho\geq\gamma^{-1}>0$ it follows that

[TABLE]

Next, we set $\delta=\max\left\{\frac{1}{2}\beta_{d}^{-2}+\displaystyle\frac{1}{2},1+\displaystyle\frac{1}{2}\right\}$, observe that $\alpha_{p}\leq\gamma$, and finally obtain the coercivity estimate

[TABLE]

which completes the proof. ∎

The above theorem implies the following stability estimate.

Corollary 2**.**

Let $(\bm{u},\bm{v},p)\in\bm{U}\times\bm{V}\times P$ be the solution of (12a)–(12c). Then there holds the estimate

[TABLE]

where $C$ is a constant independent of $\lambda,R_{p}^{-1},\alpha_{p}$ and $\|\bm{f}\|_{\bm{U}^{*}}=\sup\limits_{\bm{w}\in\bm{U}}\frac{(\bm{f},\bm{w})}{\|\bm{w}\|_{\bm{U}}},\|g\|_{P^{*}}=\sup\limits_{q\in P}\frac{(g,q)}{\|q\|_{P}}=\gamma^{-1}\|g\|$.

Remark 3**.**

We want to emphasize that the parameter ranges as specified in (24) are indeed relevant since the variations of the model parameters are quite large in many applications. For that reason, Theorem 1 is a very important and basic result that provides the parameter-robust stability of the model (12a)–(12c).

Remark 4**.**

Define

[TABLE]

Due to the theory presented in [16], Theorem 1 implies that the operator $B$ in (44) defines a norm-equivalent (canonical) block-diagonal preconditioner for the operator $A$ in (11) which is robust in all model parameters.

Remark 5**.**

Note that if $\lambda\ll R_{p}^{-1}$ and $0\leq\alpha_{p}\leq\frac{1}{\lambda}$ then $\rho=\min\{\lambda,R^{-1}_{p}\}=\lambda$ and the norms defined in (26a)–(26c) are given by

[TABLE]

Then the coercivity of $({\bm{\epsilon}}(\bm{u}),{\bm{\epsilon}}(\bm{u}))+\lambda(\operatorname{div}\bm{u},\operatorname{div}\bm{u})+(R_{p}^{-1}\bm{v},\bm{v})$ on the kernel set

[TABLE]

can be verified by direct computation. However, in this case the $H(\operatorname{div})$ inf-sup condition

[TABLE]

fails if $R_{p}^{-1}\gg\lambda$, for instance, $\lambda\eqsim 1$ and $R_{p}^{-1}\gg 1$.

Hence, in order to obtain the inf-sup condition for $(\operatorname{div}\bm{u},q)+(\operatorname{div}\bm{v},q)$, namely,

[TABLE]

the Stokes inf-sup condition from Lemma 2 has to be satisfied at the discrete level, as we can see by choosing $(\bm{u},\bm{v})=(\bm{u},0)\in\bm{U}\times\bm{V}$ in (49).

From the above observation, we conclude that we need to choose a proper space for the approximation of the displacement field $\bm{u}$, even if $\lambda$ is small, in order to satisfy the Stokes inf-sup condition stated in Lemma 2 at the discrete level. This shows that the discretization of (12a)–(12c) using on the spaces $P_{1}\times RT_{0}\times P_{0}$ can not be uniformly stable with respect to all model parameters!

Therefore, in our paper, we use $H(\operatorname{div})$ conforming spaces such as $BDM,RT$ and $BDFM$ to replace $P_{1}$. Further, in this way, we can preserve the divergence condition exactly, which means the fluid mass conservation is fully preserved. Details will be presented in the following Section 4.

4. Uniformly stable discretizations of the model

There are various discretizations that meet the requirements for the proof of the full parameter-robust stability that will be presented in this section. They include conforming as well as nonconforming methods. In general, whenever $U_{h}/P_{h}$ is a Stokes-stable pair and $V_{h}/P_{h}$ satisfies the $H(\operatorname{div})$ inf-sup condition similar to (63), the norm that we have proposed in Section 3 allows to prove the full parameter-robust stability using similar arguments as in the proof of Theorem 6. To give two popular examples, the triplet $CR/RT_{0}/P_{0}$ together with stabilization does result in a parameter-robust stable discretization of the Biot model if the norms are defined as in Section 3. The same is true for the conforming discretization based on the spaces $P_{2}/RT_{0}/P_{0}$. However, these finite element methods can not preserve the fluid mass conservation exactly and globally although they can achieve parameter-robustness.

Therefore, in this section, motivated by the works [18, 19, 20], we propose discretizations of the Biot’s model problem (12a)–(12c). These discretizations preserve the divergence condition (namely equation (9c)) exactly and globally, which means an exact conservation of the fluid mass. Furthermore, the discretizations are also locking-free when the Lamé parameter $\lambda$ tends to $\infty$. First we introduce some notations.

4.1. Preliminaries and notation

By $\mathcal{T}_{h}$ we denote a shape-regular triangulation of mesh-size $h$ of the domain $\Omega$ into triangles $\{K\}$. We further denote by $\mathcal{E}_{h}^{I}$ the set of all interior edges (or faces) of $\mathcal{T}_{h}$ and by $\mathcal{E}_{h}^{B}$ the set of all boundary edges (or faces); we set $\mathcal{E}_{h}=\mathcal{E}_{h}^{I}\cup\mathcal{E}_{h}^{B}$.

For $s\geq 1$, we define

[TABLE]

The vector functions are represented column-wise.

As we consider discontinuous Galerkin (DG) discretizations, we define some trace operators next. Let $e=\partial K_{1}\cap\partial K_{2}$ be the common boundary (interface) of two subdomains $K_{1}$ and $K_{2}$ in $\mathcal{T}_{h}$ , and $\bm{n}_{1}$ and $\bm{n}_{2}$ be unit normal vectors to $e$ pointing to the exterior of $K_{1}$ and $K_{2}$, respectively. For any edge (or face) $e\in\mathcal{E}_{h}^{I}$ and a scalar $q\in H^{1}(\mathcal{T}_{h})$, vector $\bm{v}\in H^{1}(\mathcal{T}_{h})^{d}$ and tensor $\bm{\tau}\in H^{1}(\mathcal{T}_{h})^{d\times d}$, we define the averages

[TABLE]

and jumps

[TABLE]

where $\bm{v}\odot\bm{n}=\frac{1}{2}(\bm{v}\bm{n}^{T}+\bm{n}\bm{v}^{T})$ is the symmetric part of the tensor product of $\bm{v}$ and $\bm{n}$.

When $e\in\mathcal{E}_{h}^{B}$ then the above quantities are defined as

[TABLE]

If $\bm{n}_{K}$ is the outward unit normal to $\partial K$, it is easy to check that

[TABLE]

Also, for $\bm{\tau}\in H^{1}(\Omega)^{d\times d}$ and for all $\bm{v}\in H^{1}(\mathcal{T}_{h})^{d}$, we have

[TABLE]

The finite element spaces are denoted by

[TABLE]

[TABLE]

[TABLE]

The discretizations that we consider here, define the local spaces $\bm{U}(K)/\bm{V}(K)/Q(K)$ via the triplets $RT_{l}(K)/RT_{l-1}(K)/P_{l-1}(K)$, or $BDM_{l}(K)/RT_{l-1}(K)/P_{l-1}(K)$, or $BDFM_{l}(K)/RT_{l-1}(K)/P_{l-1}(K)$. Note that for all these choices the important condition $\operatorname{div}\bm{U}(K)=\operatorname{div}\bm{V}(K)=Q(K)$ is satisfied.

We recall the basic approximation properties of these spaces: for all $K\in\mathcal{T}_{h}$ and for all $\bm{u}\in H^{s}(K)^{d}$, there exists $\bm{u}_{I}\in\bm{U}(K)$ such that

[TABLE]

4.2. DG discretization

We note that according to the definition of $\bm{U}_{h}$, the normal component of any $\bm{u}\in\bm{U}_{h}$ is continuous on the internal edges and vanishes on the boundary edges. Therefore, by splitting a vector $\bm{u}\in\bm{U}_{h}$ into its normal and tangential components $\bm{u}_{n}$ and $\bm{u}_{t}$

[TABLE]

we have

[TABLE]

implying that

[TABLE]

A direct computation shows that

[TABLE]

Therefore, the discretization of the variational problem (12a)–(12c) is given as follows. Find $(\bm{u}_{h},\bm{v}_{h},p_{h})\in\bm{U}_{h}\times\bm{V}_{h}\times P_{h}$, such that for any $(\bm{w}_{h},\bm{z}_{h},q_{h})\in\bm{U}_{h}\times\bm{V}_{h}\times P_{h}$

[TABLE]

where

[TABLE]

and $\eta$ is a stabilization parameter which is independent of $h,\lambda,R_{p}^{-1},\alpha_{p}$.

For any $\bm{u}\in H^{1}(\mathcal{T}_{h})^{d}$, we introduce the mesh dependent norms:

[TABLE]

Next, for $\bm{u}\in H^{2}(\mathcal{T}_{h})^{d}$, we define the “DG”-norm

[TABLE]

and, finally, the mesh-dependent norm $\|\bm{\cdot}\|_{\bm{U}_{h}}$ by

[TABLE]

We now summarize several results on well-posedness and approximation properties of the DG formulation, see, e.g. [19, 20].

• •

From the discrete version of Korn’s inequality we have that the norms $\|\cdot\|_{DG}$, $\|\cdot\|_{h}$, and $\|\cdot\|_{1,h}$ are equivalent on $\bm{U}_{h}$, namely,

[TABLE]

• •

The bilinear form $a_{h}(\cdot,\cdot)$, introduced in (4.2) is continuous and we have

[TABLE]

• •

For our choice of the finite element spaces $\bm{V}_{h}$ and $P_{h}$ we have the following inf-sup conditions

[TABLE]

where $\beta_{sd}$ and $\beta_{dd}$ are constant independent of the parameters $\lambda,R_{p}^{-1},\alpha_{p}$ and the mesh size $h$.

• •

We also have that $a_{h}(\cdot,\cdot)$ is coercive, and the proof of this fact parallels the proofs of similar results.

[TABLE]

where $\alpha_{a}$ is a positive constant independent of parameters $\lambda,R_{p}^{-1},\alpha_{p}$ and the mesh size $h$.

Related to the discrete problem (57a)-(57c) we introduce the bilinear form

[TABLE]

In view of the definitions of the norms $\|\cdot\|_{\bm{U}_{h}},\|\cdot\|_{\bm{V}}$ and $\|\cdot\|_{P}$, the boundedness of the $A_{h}((\bm{u}_{h},\bm{v}_{h},p_{h}),(\bm{w}_{h},\bm{z}_{h},q_{h}))$ is obvious, and we come to our second main result.

Theorem 6**.**

There exits a constant $\beta_{0}>0$ independent of the parameters $\lambda,~{}R_{p}^{-1},~{}\alpha_{p}$ and the mesh size $h$, such that

[TABLE]

Proof.

Case I:

[TABLE]

For any $(\bm{u}_{h},\bm{v}_{h},p_{h})\in\bm{U}_{h}\times\bm{V}_{h}\times P_{h}$, by the first inequality in (63), there exists

[TABLE]

Choose

[TABLE]

where the constant $\delta$ will be determined later.

We verify first the boundedness of $(\bm{w}_{h},\bm{z}_{h},q_{h})$ by $(\bm{u}_{h},\bm{v}_{h},p_{h})$.

By (68), the equivalence between norms $\|\cdot\|_{DG}$ and $\|\cdot\|_{1,h}$, namely (61), and noting that $\gamma^{-1}\leq\rho=\lambda$, we have

[TABLE]

Therefore, by taking into account that $1\leq\lambda$ we get for $\bm{w}_{h}$ the estimate

[TABLE]

Obviously we have $\|\bm{z}_{h}\|_{\bm{V}}=\delta\|\bm{v}_{h}\|_{\bm{V}}$ and it remains to bound $q_{h}$. From

[TABLE]

it follows that

[TABLE]

Next we establish the coercivity of $A_{h}((\bm{u}_{h},\bm{v}_{h},p_{h}),(\bm{w}_{h},\bm{z}_{h},q_{h}))$.

Using the definition of $(\bm{w}_{h},\bm{z}_{h},q_{h})$ and (68), we find

[TABLE]

Next we apply Cauchy’s inequality, use the coercivity and the continuity of $a_{h}(\cdot,\cdot)$, the equivalence of the norms $\|\cdot\|_{DG}$ and $\|\cdot\|_{1,h}$, and (68), to get

[TABLE]

Now letting $\epsilon_{1}=\frac{1}{2}C^{-2}_{1}C_{0}^{-2}\beta_{sd}^{2},\epsilon_{2}=\epsilon_{3}=\epsilon_{4}=\frac{1}{2}$, and noting that $\rho^{-1}\leq\gamma,\lambda\geq\gamma^{-1}>0$, we obtain

[TABLE]

Next, setting $\delta:=\max\{C^{2}_{1}C^{2}_{0}\beta_{sd}^{-2}+\frac{1}{2}\alpha^{-1}_{a},2+\frac{1}{2}\}$ and noting that $\alpha_{p}\leq\gamma,\lambda\geq 1$, we derive the coervicity estimate

[TABLE]

Case II:

[TABLE]

For any $(\bm{u}_{h},\bm{v}_{h},p_{h})\in\bm{U}_{h}\times\bm{V}_{h}\times P_{h}$, by the second inequality in (63), there exists

[TABLE]

Choose

[TABLE]

where $\delta$ is a constant which will be specified later.

Again we first verify the boundedness of $(\bm{w}_{h},\bm{z}_{h},q_{h})$ by $(\bm{u}_{h},\bm{v}_{h},p_{h})$. We note that $\|\bm{w}_{h}\|_{\bm{U}_{h}}=\delta\|\bm{u}_{h}\|_{\bm{U}_{h}}$.

From (75) and noting that $\gamma^{-1}\leq\rho=R_{p}^{-1}$ we have

[TABLE]

Hence, we get the boundedness of $\bm{z}_{h}$, that is

[TABLE]

Again we have the boundedness for $q_{h}$ according to (73).

In what follows we show the coercivity of $A_{h}((\bm{u}_{h},\bm{v}_{h},p_{h}),(\bm{w}_{h},\bm{z}_{h},q_{h}))$ in Case II.

Using the definition of $(\bm{w}_{h},\bm{z}_{h},q_{h})$ and (75), we find

[TABLE]

Next, we apply Cauchy’s inequality, use (75), and the coercivity of $a_{h}(\cdot,\cdot)$, to get

[TABLE]

Now, letting $\epsilon_{1}=\beta_{dd}^{2},\epsilon_{2}=\epsilon_{3}=\frac{1}{2}$, and noting that $\rho^{-1}\leq\gamma,\lambda\geq\rho\geq\gamma^{-1}>0$, we obtain

[TABLE]

Finally, we choose $\delta:=\max\{\frac{1}{2}\alpha_{a}^{-1},\frac{1}{2}\beta_{dd}^{-2}+\frac{1}{2},1+\frac{1}{2}\}$, note that $\alpha_{p}\leq\gamma$, and conclude the coervicity of the bilinear form, i.e.,

[TABLE]

This completes the proof. ∎

From the above theorem, we get the following stability estimate.

Corollary 7**.**

Let $(\bm{u}_{h},\bm{v}_{h},p_{h})\in\bm{U}_{h}\times\bm{V}_{h}\times P_{h}$ be the solution of (57a)-(57c), then we have the estimate

[TABLE]

where $\|\bm{f}\|_{\bm{U}_{h}^{*}}=\sup\limits_{\bm{w}_{h}\in\bm{U}_{h}}\frac{(\bm{f},\bm{w}_{h})}{\|\bm{w}_{h}\|_{\bm{U}_{h}}},\|g\|_{P^{*}}=\sup\limits_{q_{h}\in P_{h}}\frac{(g,q_{h})}{\|q_{h}\|_{P}}$ and $C_{2}$ is a constant independent of $\lambda,R_{p}^{-1},\alpha_{p}$ and mesh size $h$.

Remark 8**.**

Denote by $A_{h}$ the operator induced by the bilinear form (65), namely

[TABLE]

and define

[TABLE]

Then due to the theory presented in [16], Theorem 6 implies that the norm-equivalent (canonical) block-diagonal preconditioner $B_{h}$ for $A_{h}$ is parameter-robust, which means that the condition number $\kappa(B_{h}A_{h})$ is uniformly bounded with respect to the parameters $\lambda,R_{p}^{-1},\alpha_{p}$ in the ranges (24) and with respect to the mesh size $h$.

5. Error estimates

In this section, we derive the error estimates that follow from the results presented in Section 4. Let $\Pi_{B}^{\operatorname{div}}:H^{1}(\Omega)^{d}\mapsto\bm{U}_{h}$ be the canonical interpolation operator. We also denote the $L^{2}$-projection on $P_{h}$ by $Q_{h}$. The following Lemma, see [19], summarizes some of the properties of $\Pi_{B}^{\operatorname{div}}$ and $Q_{h}$ needed for our proof.

Lemma 3**.**

For all $\bm{w}\in H^{1}(K)^{d}$ we have

[TABLE]

Theorem 9**.**

Let $(\bm{u},\bm{v},p)$ be the solution of (12a)–(12c) and $(\bm{u}_{h},\bm{v}_{h},p_{h})$ be the solution of (57a)–(57c). Then the error estimates

[TABLE]

and

[TABLE]

hold, where $C_{e,u},C_{e,p}$ are constants independent of $\lambda,R_{p}^{-1},\alpha_{p}$ and the mesh size $h$.

Proof.

Subtracting (57a)–(57c) from (12a)–(12c) and noting the consistency of $a_{h}(\cdot,\cdot)$, we have that for any $(\bm{w}_{h},\bm{z}_{h},q_{h})\in\bm{U}_{h}\times\bm{V}_{h}\times P_{h}$

[TABLE]

Let $\bm{u}_{I}=\Pi_{B}^{div}\bm{u}\in\bm{U}_{h},p_{I}=Q_{h}p\in P_{h}$. Now for arbitrary $\bm{v}_{I}\in\bm{V}_{h}$, from (84)–(86), noting that $\operatorname{div}\Pi_{B}^{div}=Q_{h}\operatorname{div}$ and $\operatorname{div}\bm{U}_{h}=\operatorname{div}\bm{V}_{h}=P_{h}$, we conclude

[TABLE]

Next, since $(\bm{u}_{I}-\bm{u}_{h})\in\bm{U}_{h},(\bm{v}_{I}-\bm{v}_{h})\in V_{h},(p_{I}-p_{h})\in P_{h}$, by the stability result (66) for the discrete problem (57a)–(57c), we obtain

[TABLE]

and

[TABLE]

Hence, using the boundedness of $a_{h}(\cdot,\cdot)$, the second inequality in Lemma 3, and triangle inequality, we arrive at

[TABLE]

and

[TABLE]

∎

Remark 10**.**

From the above theorem, we can see that the discretizations are locking-free.

6. Conclusions

This paper presents the stability analysis of a classical three-field formulation of Biot’s consolidation model where the unknown variables are the displacements, fluid flux (Darcy velocity), and pore pressure. Specific parameter-dependent norms provide the key to establish the parameter-robust stability of the continuous problem. This allows for the construction of a parameter-robust block diagonal preconditioner in the framework of operator preconditioning. Discretizations that fully preserve the fluid mass conservation are designed. Further, both discrete parameter-robust stability and locking-free error estimates are proved.