A four-field mixed finite element method for Biot's consolidation problems

1.   Introduction

In [4,5], Biot's equations coupling mechanical and flow problem were introduced to describe solid consolidation phenomena by coupling the solid and fluid variables. As an example, squeezing water out of an elastic porous medium can be understood as a consolidation process. In [20,15], fluid movement within soft tissue was considered with the corresponding consolidation models, which were crucial for treatment in solid tumors, such as, delivery of drugs and nutrients, elastography.

The Biot's consolidation equations have been solved numerically by various finite element methods in literature. In [17], an error analysis for semi-discrete and fully-discrete finite element method of Biot's models was established, and the backward Euler discretization was used in fully-discrete scheme. The formulation involved displacement and pressure as its unknown variables. A least-squares four-field finite element method was introduced in [10], and the unknown variables were displacement, stress tensor, pressure and flux. The stress tensor and flux were approximated using Raviart-Thomas elements. With Raviart-Thomas approximating space for flux and continuous Galerkin method for displacement, a three field formulation was derived to approximate the Biot's consolidation models in [23,24]. In [28], a four-field mixed finite element method was proposed and analyzed where the unknowns functions include the total stress tensor, displacement, fluid flux, and pore pressure. By introducing the stabilization term, a stabilized lowest order finite element method for three-field poroelasticity was developed in [3] by using piecewise constant element for pressure and piecewise linear element for displacement and flux. A three field finite element method using Crouzeix-Raviart element for displacement and Raviart-Thomas element for flux was presented in [9], with optimal order of convergence for backward Euler fully discrete scheme. A four-field mixed finite element formulation with displacement, total pressure, pressure and flux as unknowns was discussed in [11] where a finite volume method was designed for the displacement.

For simplicity, we consider the Biot's consolidation model that seeks the displacement ${\bf{u}}$ and the pore pressure $p$ such that

where $\Omega$ is an open bounded polygonal or polyhedral domain in $\mathbb{R}^d,\ d = 2,3$ with boundary $\Gamma_D\cup\Gamma_N = \partial \Omega$. Let ${\bf{f}}$ and $g$ be known functions, and $\beta$ be the data on boundary $\Gamma_N$ with positive measure in $\mathbb{R}^{d-1}$. Here, $\varepsilon({\bf{u}}) = \frac{1}{2}(\nabla{\bf{u}}+(\nabla{\bf{u}})^T)$ is the strain tensor, $\mu$ and $\lambda$ are the elastic parameters, $\kappa$ is the permeability and $\chi$ is the average microfiltration coefficient, and ${\bf{n}}$ denotes the unit outward normal vector on $\partial \Omega$. Assume the parameters $\mu, \ \lambda, \ \kappa$ are strictly positive and $\chi \ge 0$ is non-negative. For $x\in\Omega$, the initial conditions are given by

By introducing the total stress $z = -\lambda\nabla\cdot{\bf{u}}+p$ and the fluid flux (Darcy velocity) ${\bf{q}} = -\kappa\nabla p$, the initial condition for the total stress would be given by

The total stress $z$ has been introduced for the Biot's model with three fields in [22,13,14]. In particular, we have used the displacement, total stress and pressure as the unknown variables for the Biot's model in [25] and presented a finite element method and a convergence theory for the fully-discrete scheme with backward Euler discretization in time.

The goal of this paper is to devise and analyze a four-field finite element method for the Biot's model (1)-(5) which consists of displacement, total stress, pressure and the Darcy flux for the system. To describe a weak formulation with four field variables, we denote two Sobolev spaces as follows

For simplicity, we shall adopt the notation and the definition of Sobolev spaces in [1,26].

The weak form of (1)-(5) reads as follows: find ${\bf{u}}\in[H^1_{0,D}(\Omega)]^d$, $z\in L^2(\Omega)$, $p \in L^2(\Omega)$ and ${\bf{q}}\in H(div; \Omega)$ such that

Our four-field mixed finite element method is designed by using the Raviart-Thomas element $(RT_k)$ [27,19] for the flux variable ${\bf{q}}$ and the pressure $p$, with an employment of conforming finite elements for the displacement ${\bf{u}}$ and the total stress variable $z$ satisfying the inf-sup condition of Babuška [2] and Brezzi [7]. The time direction is discretized by using the Crank-Nicolson (CN) scheme. The satisfaction of the inf-sup condition for the displacement ${\bf{u}}$ and the total stress variable $z$ ensures a locking-free nature of the numerical scheme.

The paper is organized as follows. In Section 2, we describe the semi-discrete and fully-discrete schemes for a four-field mixed finite element method. In Section 3, we derive some error estimates for the semi-discrete numerical scheme, while the main result of error estimates for the fully-discrete scheme shall be established in Section 4. In section 5, we present some numerical results for several benchmark problems with various boundary conditions. Our numerical experiments include an implementation of the Hood-Taylor element for the approximation of the displacement and total stress.

2.   Four-field mixed finite element method

Let $\mathcal{T}_h$ be a partition of the domain $\Omega$ into triangular or tetrahedral elements satisfying the quasi-uniform condition. For each element $T\in\mathcal{T}_h$, denote by $h_T$ its diameter and $h = \max_{T\in\mathcal{T}_h}h_T$ the mesh size of the triangulation $\mathcal{T}_h$. Denote by $\mathcal{E}_{I}$ the set of interior edges/faces and $\partial T$ the boundary of the element $T$. Let ${P}_{k}(T)$ be the space of polynomials of degree less than or equal to $k$ in $T$. We consider the following finite element spaces for $k\geq 0$

Define the Stokes projection operator $Q_h^{{\bf{u}}}\times Q_h^z: \ [H^1_{0,D}(\Omega)]^d\times L_{2}(\Omega) \rightarrow {\bf{U}}_h\times Z_h$ by

Next, denote by $Q_h^{{\bf{q}}}: H_{0,\Gamma_N}(div; \Omega) \rightarrow {\bf{X}}_{h}$ the Fortin projection operator and $Q_h^{p}: L_{2}(\Omega) \rightarrow P_h$ the $L^2$ projection operator satisfying the following properties:

Consider the following two lemmas that we will need to prove error estimates in the next section.

Lemma 2.1. [18] Let $Q_h^{{\bf{u}}}\times Q_h^z$ be the Stokes projection operator defined by (8). For $i = 0,1$ and convex $\Omega$, there exists a constant $C$ such that

Lemma 2.2. (see Proposition 3.6 in [8]) For the Fortin projection operator $Q_h^{{\bf{q}}}$ and the $L^2$ projection operator $Q_h^p$ satisfying (9), there exists a constant $C$ such that

Given a suitable approximation of the initial conditions (6) ${\bf{u}}_h(0) = Q_h^{{\bf{u}}}\varphi$, $z_h(0) = Q_h^z(-\lambda\nabla\cdot\varphi+\phi)$, $p_h(0) = Q_h^p\phi$ and ${\bf{q}}_h(0) = -Q_h^{{\bf{q}}}(\kappa\nabla\phi)$, the semi-discrete scheme for the Biot's model problem (1)-(5) seeks ${\bf{u}}_h\in {\bf{U}}_{h}$, $z_h\in Z_{h}$, $p_h\in P_{h}$ and ${\bf{q}}_h\in {\bf{X}}_h$ such that

To describe a fully-discrete numerical method with Crank-Nicolson discretization in time, denote by $\tau$ the time step, and $t^n = \tau n$ the time level, where $n$ is a non-negative integer. Then, given a suitable approximation of the initial conditions (6) ${\bf{u}} ^0 = Q_h^{{\bf{u}}}\varphi$, $z^0 = Q_h^z(-\lambda\nabla\cdot\varphi+\phi)$, $p^0 = Q_h^p\phi$ and ${\bf{q}}^0 = -Q_h^q(\kappa\nabla\phi)$, the fully-discrete scheme is to find ${\bf{u}}^n\in {\bf{U}}_{h}$, $z^n\in Z_{h}$, $p^n\in P_{h}$ and ${\bf{q}}^n\in {\bf{X}}_h$ such that

where $\bar{\partial}_t z^n: = \frac{z^n-z^{n-1}}{\tau}$ and $f^n: = f(x, t^n),\ x\in\Omega$.

Lemma 2.3. The numerical solution $({\bf{u}}_h, z_h, p_h, {\bf{q}}_h)$ for semi-discrete scheme (10) is existence and uniqueness.

Proof. It suffices to show that the trivial functions are the only solution to the homogeneous problem. Observe that the semi-discrete scheme for the homogeneous problem seeks ${\bf{u}}_h\in {\bf{U}}_{h}$, $z_h\in Z_{h}$, $p_h\in P_{h}$ and ${\bf{q}}_h\in {\bf{X}}_h$ such that

By taking the time derivative for the second equation and then choosing ${\bf{v}} = {\bf{u}}_{ht}, w_z = z_h, w_p = p_h, w_{{\bf{q}}} = {\bf{q}}_h$ in (12), we have from summing up all the equations

Integrating over $(0, t)$, then from ${\bf{u}}_h(0) = 0$ and $z_h(0)-p_h(0) = 0$, we have

Since $\mu$, $\lambda$, and $\kappa$ are strictly positive and $\chi\ge 0$, it follows that $\varepsilon({\bf{u}}_h(t)) = 0$, $z_h(t)-p_h(t) = 0$, and ${\bf{q}}_h = 0$. With the use of Korn's inequality [21,6], we have ${\bf{u}}_h(t) = 0$. Finally, from the inf-sup condition for the Raviart-Thomas element and the fourth equation in (12) we have $p_h = 0$.

3.   Error estimates

We note that the finite element pair $Z_h\times{\bf{U}}_h$ satisfies the inf-sup condition [2,8], i.e. there exists a constant $C>0$ such that for any $w_z \in Z_h$ one has

Next, the Raviart-Thomas space $P_h\times{\bf{X}}_h$ satisfies the following inf-sup condition [26]

3.1.   Error estimates for semi-discrete scheme

Denote the error of the displacement for the semi-discrete scheme by $e_h^{{\bf{u}}} = {\bf{u}}-{\bf{u}}_h = \eta_h^{{\bf{u}}}+\xi_h^{{\bf{u}}}$, where $\eta_h^{{\bf{u}}} = {\bf{u}}-Q_h^{{\bf{u}}}{\bf{u}}$ and $\xi_h^{{\bf{u}}} = Q_h^{{\bf{u}}}{\bf{u}}-{\bf{u}}_h$. Similarly, we denote the errors of the total stress, the pressure, the flux by $e_h^z, \ e_h^p,\ e_h^{{\bf{q}}}$, respectively.

Theorem 3.1. Let $({\bf{u}}, z, p, {\bf{q}})$ be the solution of (7) and $({\bf{u}}_h, z_h, p_h, {\bf{q}}_h)$ the numerical solution of (10). There exists a constant $C$ such that for each $t\in(0, \bar{T}]$

and

Proof. Substituting (10) into (7) leads to

Using the properties of projections and differentiating the second term with respect to time $t$, we deduce that

where we have used the property of $\nabla\cdot w_{\bf{q}}\in P_h$. By choosing ${\bf{v}} = \xi^{{\bf{u}}}_{ht}, w_z = \xi^z_h, w_p = \xi^p_h, w_{{\bf{q}}} = \xi^{{\bf{q}}}_h$ in (15) and summing up, we get

Then, integrating over $(0, t)$ yields

Using the Gronwall Lemma [26], we obtain

Moreover, from the first equation of (15), we find $2\mu(\varepsilon(\xi^{{\bf{u}}}_h), \varepsilon({\bf{v}})) = (\xi^z_h,\nabla\cdot{\bf{v}})$. Hence, combining with inf-sup condition (13), it follows

On the other hand, from differentiating the fourth equation of (15) with respect to $t$ and taking ${\bf{v}} = \xi^{{\bf{u}}}_{ht}, w_z = \xi^z_{ht}, w_p = \xi^p_{ht}, w_{{\bf{q}}} = \xi^{{\bf{q}}}_h$, we arrive at

Integrating the equation over $(0, t)$, we have from $\xi_h^p(0) = 0$, $\xi_h^{\bf{q}}(0) = 0$

Using the Gronwall Lemma [26], we get

The desired error estimates now stem from (16)-(18).

Using the inf-sup condition (13) similar to the process in (17), we have

Thus, in view of (16)-(19), we have derived the following result.

Corollary 1. Let $({\bf{u}}, p, z, {\bf{q}})$ be the solution of (7) and $({\bf{u}}_h, p_h, z_h, {\bf{q}}_h)$ the solution of (10). There exists a constant $C$ such that

and

Finally, combining Theorem 3.1 and Lemmas 2.1, 2.2, we arrive at the following error estimates.

Theorem 3.2. Let $({\bf{u}}, z, p, {\bf{q}})$ be the solution of (7) and $({\bf{u}}_h, z_h, p_h, {\bf{q}}_h)$ the numerical solution of (10). Assume ${\bf{u}}\in L^{\infty}(0, \bar{T}; [H^{k+3}(\Omega)]^d)$, $z\in H^{1}(0, \bar{T}; H^{k+1}(\Omega))$, $p\in H^{1}(0, \bar{T}; H^{k+1}(\Omega))$ and ${\bf{q}}\in L^{2}(0, \bar{T}; [H^{k+1}(\Omega)]^d)$. Then, there exists a constant $C$ such that for each $t\in(0, \bar{T}]$

and

3.2.   Error estimates for fully-discrete scheme

To derive an error estimate for the fully-discrete scheme (11), we denote the error $e_{{\bf{u}}}^n = {\bf{u}}(t^n)-{\bf{u}}^n = \eta_{{\bf{u}}}^n+\xi_{{\bf{u}}}^n$, where $\eta_{{\bf{u}}}^n = {\bf{u}}(t^n)-Q_h^{{\bf{u}}}{\bf{u}}(t^n)$ and $\xi_{{\bf{u}}}^n = Q_h^{{\bf{u}}}{\bf{u}}(t^n)-{\bf{u}}^n$. The error for other variables are denoted in a similar way.

Theorem 3.3. Let $({\bf{u}}, p, z, {\bf{q}})$ be the solution of (7) and $({\bf{u}}^n, p^n, z^n, {\bf{q}}^n)$ the solution of (11). There exists a constant $C$ such that

and

Here we denote $RF = \int_{0}^{t^N}\Vert\eta_{ht}^z-\eta_{ht}^p\Vert^2ds+\bar{T}\max\limits_{0\leq n\leq N}\Vert\eta_{{\bf{q}}}^n\Vert^2$.

Proof. By subtracting (7) from (11), for each ${\bf{v}}_h\in{\bf{U}}_h$, $w_z\in Z_h$, $w_p\in P_h$ and $w_{{\bf{q}}}\in{\bf{X}}_h$, we have

Using $\frac{z_t(t^n)-z_t(t^{n-1})}{2}-\bar{\partial}_t z(t^n) = -\frac{1}{2\tau}\int_{t^{n-1}}^{t^n}(t^n-s)(t^{n-1}-s)z_{ttt}ds$, the third one of the above equations can be rewritten as

Hence, combining the above equations with the properties of projections (8) and (9), we get

Therefore, by taking ${\bf{v}}_h = \bar{\partial}_t\xi_{{\bf{u}}}^n$, $w_z = \frac{\xi_z^n+\xi_z^{n-1}}{2}$, $w_p = \frac{\xi_p^n+\xi_p^{n-1}}{2}$ and $w_{\bf{q}} = \frac{\xi_{\bf{q}}^n+\xi_{\bf{q}}^{n-1}}{2}$ in (20), we deduce

Adding the above equations leads to

Here, we have used the identity

Choosing $\epsilon_1 = \frac{1}{2}$ and $\epsilon_2\lambda^{-1} = \frac{\chi}{2}$ leads to

Adding from $1$ to $N$, we find from $\xi_{{\bf{u}}}^0 = 0$ and $\xi_z^0-\xi_p^0 = 0$

where we have used

Now from the discrete Gronwall Lemma [26], we conclude

Furthermore, we rewrite the above inequality as

Similar to the proceed in the error estimate for the semi-discrete scheme, from the inf-sup condition (13) we have

Therefore, with the help of (22) and (23), the error estimates for $e_{{\bf{u}}}^N$ and $e_z^N$ are completed for their derivation. Note that, from the triangle inequality and (23), we can deduce that

which, by combining (22) and (24), yields the error estimate for $e_p^N$.

On the other side, taking ${\bf{v}}_h = \bar{\partial}_t\xi_{{\bf{u}}}^n$, $w_z = \bar{\partial}_t\xi_z^n$, $w_p = \bar{\partial}_t\xi_p^n$ and $w_{\bf{q}} = \frac{\xi_{\bf{q}}^n+\xi_{\bf{q}}^{n-1}}{2}$ in (20), one finds

Adding the above equations and using (21) lead to

Combining the first equation in (24) with the inf-sup condition (13), we get

Hence, by taking $\lambda^{-1}\epsilon_1(1+C)\leq \min\{2\mu,\frac{\lambda^{-1}u}{2}\}$, i.e $\epsilon_1\leq \frac{1}{1+C}\min\{2\mu\lambda^{-1},\frac{1}{2}\}$, the above inequality can be reformulated as

Adding from $1$ to $N$, we obtain

With the help of the error estimate (22), we obtain the desired error estimate for the flux $e^N_{{\bf{q}}}$.

Combining Theorem 3.3 with Lemmas 2.1, 2.2 leads to the following results.

Theorem 3.4. Let $({\bf{u}}, p, z, {\bf{q}})$ be the solution of (7) and $({\bf{u}}^n, p^n, z^n, {\bf{q}}^n)$ the numerical solution of (11). Assume that the displacement ${\bf{u}}\in L^{\infty}(0, \bar{T}; [H^{k+3}(\Omega)]^d)$, the total stress $z\in L^{\infty}(0, \bar{T}; H^{k+1}(\Omega))$, $z_{ttt}\in L^{2}(0, \bar{T}; L^{2}(\Omega))$, the pressure $p\in L^{\infty}(0, \bar{T}; H^{k+1}(\Omega))$, $p_{ttt}\in L^{2}(0, \bar{T}; L^{2}(\Omega))$ and the flux ${\bf{q}}\in H^{1}(0, \bar{T}; [H^{k+1}(\Omega)]^d)$. Then, there exists a constant $C$ such that

and

where $R1 = \int_{0}^{t^N}(\Vert {z_t}\Vert_{k+1}^2+\Vert{p_t}\Vert_{k+1}^2)ds +\max\limits_{0\leq n\leq N}\Vert{\bf{q}}(t^n)\Vert^2_{k+1}$.

Remark 1. It should be noted that the finite element pairs $(Z_h, {\bf{U}}_h)$ can be replaced by any finite element spaces that satisfy the inf-sup condition (13). In particular, the Hood-Taylor element $([P_{k+2}]^d,P_{k+1})$ could be employed in the place of $(Z_h, {\bf{U}}_h)$ for which all the error estimates can be derived without any difficulty. Analogously, the Raviart-Thomas element $(P_k, RT_k)$ can be substituted by other stable elements in the mixed finite element literature, and the theory remains true. Details are left to interested readers as an exercise.

4.   Numerical examples

In this section, we shall present some numerical results for the four-field mixed finite element method implemented on uniform triangulations for the Biot's consolidation model in the two dimensional space. The Darcy velocity and the fluid pressure are approximated by the Raviart-Thomas element of lowest order. Two examples of stable elements are employed for the elasticity equation: the first one makes use of $C^0$-conforming piecewise quadratic functions for the displacement variable and piecewise constant for the total stress. The second such example is based on the Hood-Taylor element. For simplicity, the first example shall be referred as $([P_2]^2,P_0,P_0, RT_0)$, and the second one as $([P_2]^2,P_1,P_0, RT_0)$. Denote by $|||e_{\bf{u}}^n|||^2 = 2\mu\Vert \varepsilon(e_{\bf{u}}^n)\Vert^2$ the energy norm for the displacement.

Example 1. The domain is given by $\Omega = (0,1)^2$ and the terminal time is $\bar{T} = 1$. The boundary condition is of full Dirichlet (i.e., $\Gamma_D = \partial \Omega$). The exact solution is given by

and

The body force function ${\bf{f}}$, the forced fluid $g$, the initial value and Dirichlet boundary value are determined by the exact solution. The parameters are set as $\kappa = 10^{-2}$, $\mu = 10^{-2}$ and $\lambda = 10^{-2}$, $\chi = 1$ in this example.

Table 1 reveals that the convergence for the finite element solution is of order $O(h)$ for the displacement in energy norm, and $O(h^2)$ in the $L^2$ norm. We also see from Table 1 that the convergence for the total stress, pressure and flux in $L^2$ norm appear to be $O(h^2)$. Note that, when the mesh size $h$ is sufficiently small, Table 2 shows a second order convergence $O(\tau^2)$ in time. Overall, these numerical results are in good agreement with our error estimates.

Example 2. The domain is given by $\Omega = (0,1)^2$ and the final time is $\bar{T} = 1$. The exact solution for the displacement and the pore pressure is given as [12]

and

The body force function ${\bf{f}}$, the forced fluid $g$ and initial conditions are determined by the exact solution. The parameters are set as $\kappa = 10$, $\mu = 10$ and $\lambda = 1$, $\chi = 0$ in this example. Full Dirichlet boundary condition is imposed on the boundary, which is non-homogeneous.

Table 3 illustrates the numerical performance of the four-field mixed finite element method with Example 2. It can be seen that the convergence for the corresponding numerical displacement is of order $O(h)$ and $O(h^2)$ in energy norm and $L^2$ norm, respectively. The almost $O(h^2)$ convergence for the total stress, pressure and flux in $L^2$ norm is also showed in Table 3. These numerical results confirm the theoretical convergence for problems with non-homogeneous Dirichlet boundary conditions.

We further tested the numerical performance of the four-field mixed finite element method for Example 2 with mixed Dirichlet and Neumann boundary conditions. In this test, the Neumann boundary condition is set at $\Gamma_N = \{(x,y) \mid x = 1,y\in (0,1)\}$ and the rest is of Dirichlet type (i.e., $\Gamma_D = \partial\Omega\setminus\Gamma_N$). The numerical results are presented in Table 4, which shows that our method is accurate and efficient for the cases with mixed Dirichlet and Neumann boundary conditions. It can be seen that a quadratic convergence for the total stress, the pressure, and the flux in $L^2$ norm was achieved in Table 4.

The Hood-Taylor element of $(P_1, [P_2]^2)$ was also implemented for the elastic part of the Biot's consolidation model in our four-field mixed finite element method. In other words, in the numerical scheme (11), the finite element pair $(Z_h, {\bf{U}}_h)$ is chosen as the lowest order Hood-Taylor element $(P_1, [P_2]^2)$ and the mixed finite element pair $(P_h, {\bf{X}}_h)$ is kept as the lowest order Raviat-Thomas element $(P_0, RT_0)$. The exact solutions and the parameters are the same as in Example 2. Table 5 shows the numerical results when full Dirichlet boundary condition is imposed in the model; i.e., $\Gamma_D = \partial \Omega$. It can be seen that the convergence for the displacement in energy norm is of order $O(h^2)$, which is an expected improvement over the $([P_2]^2,P_0,P_0, RT_0)$ element.

The numerical results with the use of Hood-Taylor element for the mixed Dirichlet and Neumann boundary condition are shown in Table 6. In this numerical test, the Neumann boundary condition is imposed on $\Gamma_N = \{(x, y) \mid x = 1,y\in (0,1)\}$ and the of the boundary assumed Dirichlet data. It can be seen in Table 6 that the convergence of the displacement approximation in the energy norm is of the order $O(h^2)$ while the order of convergence for the total stress, the pressure, and the Darcy flux in $L^2$ norm is higher than two.

5.   Concluding remarks

In the paper, we presented and analyzed a four-field mixed finite element method using the Raviart-Thomas element for Biot's consolidation model. The new numerical scheme retains the important mass conservation property for the flow equation. The novelty of the four-field mixed finite element method lies in the use of the total stress for the elastic equation, the mixed form of the flow equation, and the Crank-Nicolson scheme in the time discretization for the consolidation model. Although each of the total stress and Darcy flux have been considered on its own in other applications, the combined four-field method with the total stress and flux and Crank-Nicolson scheme is novel for the Biot's consolidation equations.

Acknowledgments

The research of Wenya Qi was partially supported by China Scholarship Council, Grant Number: 201906180039 and National Natural Science Foundation of China (Grant No.11471150). The research of Junping Wang was supported by the NSF IR/D program, while working at National Science Foundation. However, any opinion, finding, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.