Locking free staggered DG method for the Biot system of poroelasticity on general polygonal meshes

The Biot system of poroelasticity [3, 41] models fluid flow within deformable porous media, and has been extensively studied in the literature due to its broad range of applications such as geosciences, carbon sequestration and biomedical applications. The system is a fully coupled system, consisting of an equilibrium equation for the solid and a mass balance equation for the fluid. A great amount of effort has been devoted to developing efficient numerical schemes for the Biot system. The two-field displacement-pressure formulation is studied by numerous numerical methods such as finite difference method [20], finite volume method [35], finite element methods [31, 32, 33, 40], coupling of HHO method with SWIP method [4], weak Galerkin method [22] and HDG method [19]. The three-field displacement-pressure-Darcy velocity formulation with various couplings of continuous and discontinuous Galerkin methods, and mixed finite element methods is studied in [37, 38, 45, 44, 36, 29, 47, 27]. Very recently, a stabilized hybrid mixed finite element method uniformly stable with respect to the physical parameters is proposed in [34]. Alternatively, fully-mixed formulations of the Biot system has also drawn great attention [24, 46, 26, 1]. Typical difficulty encountered in the design of numerical schemes for the Biot system is to avoid locking and to get estimates independent of the storativity coefficient.

Therefore, the purpose of this paper is to develop and analyze a staggered discontinuous Galerkin (DG) method on general polygonal meshes for the Biot system of poroelasticity based on a five-field formulation. As a new generation of numerical schemes, staggered DG method offers many salient features such as the easy treatment of general polygonal elements possibly including hanging nodes, the local and global mass conservation and superconvergence. All these features make it a good candidate for practical applications, thus it has been extensively studied for a wide range of partial differential equations arising from physical applications [12, 13, 15, 14, 23, 9, 8, 10, 11, 16, 18, 51, 54, 53, 52, 48, 50, 49]. Therefore, it is of great importance to develop a staggered DG method for the Biot system of poroelasticity. To this end, we will first extend staggered DG method developed in [53] with piecewise constant approximations to arbitrary polynomial orders, where the symmetry of stress is weakly imposed via the carefully designed approximation space for rotation. Then it is combined with staggered DG method proposed in [51] to discretize the Biot system of poroelasticity.

The resulting formulation is fully mixed with stress tensor, the fluid flux, the rotation, the displacement and the pressure as the unknowns. At a first sight this appears to increase the computational burden by introducing additional variables compared to methods that only involve displacement and pressure. However, the introduction of these variables enables us to achieve all the variables of physical interest directly without resorting to postprocessing, which is actually more accurate. We emphasize that the rotation variable introduced here is a scalar field, which can reduce the size of the global system. Another advantage of our method is that no further stabilization or penalty term is needed in contrast to other DG methods. In addition, local conservation is preserved on the dual mesh. It is worth mentioning that the mass matrix is block diagonal, which allows local elimination of the stress, velocity and rotation. As such, we can exploit local elimination to reduce the global saddle point system to displacement-pressure system. Since this procedure is similar to that proposed in [1], where block diagonal matrices can be generated thanks to the vertex quadrature rule utilized, we will not repeat here. Instead, we will introduce fixed stress splitting scheme to reduce the size of the global system. Fixed stress splitting scheme is a popular method for iteratively solving the Biot system, but has not been addressed for the formulation we exploited. Thus, we will propose a fixed stress splitting scheme for our formulation and prove the linear convergence of the given scheme. The main ingredient is to decompose the global system into two subproblem, then we can solve the flow problem and the mechanics problem separately. In particular, as mentioned before, the mass matrix is block diagonal, thus we can also apply local elimination for each problem, consequently, we can obtain a system which is solely based on displacement and pressure.

We analyze the unique solvability, stability and convergence estimates for the semidiscrete continuous-in-time and the fully discrete methods. Note that a different approach from [28] is pursued to prove the inf-sup condition for the elasticity equation. The optimal convergence rates for all the variables in their natural norms are analyzed. In addition, the error estimates are independent of the anisotropy of the permeability of the porous media. More importantly, the stability and convergence estimates are independent of the storativity coefficient $c_{0}$ and are valid for $c_{0}=0$. Moreover, like the linear elasticity [53], the error bounds are robust for nearly incompressible materials. The numerical experiments including cantilever bracket benchmark problem presented verify the optimal convergence rates and locking free property of the proposed method.

The rest of the paper is organized as follows. In the next section, we describe the model problem and the associated weak formulation. Then in section 3, we derive the discrete formulation for the Biot system of poroelasticity and prove the unique solvability. Convergence analysis for semi-discrete and fully discrete methods is presented in section 4. In section 5, a fixed stress splitting scheme is introduced and analyzed. Several numerical experiments are carried out in section 6 to confirm the proposed theories. Finally, a conclusion is given.

In this section we introduce the poroelasticity system and its fully mixed weak formulation, which lays foundation for the derivation of our staggered DG scheme. Let $\Omega\subset\mathbb{R}^{2}$ be a simply connected bounded domain, occupied by a poroelastic media saturated with fluid. Then the stress-strain constitutive relationship for the poroelastic body is given by

$\displaystyle A\sigma_{e}=\epsilon(u),$

(2.1)

where $\epsilon(u)=\frac{1}{2}(\nabla u+\nabla u^{T})$ and $\sigma_{e}=2\mu\epsilon(u)+\lambda\nabla\cdot uI$. Here $A$ is a symmetric, bounded and uniformly positive definite linear operator representing the compliance tensor, $\sigma_{e}$ is the elastic strain, $u$ is the solid displacement.

In the case of a homogeneous and isotropic body, we have

$\displaystyle A\sigma=\frac{1}{2\mu}(\sigma-\frac{\lambda}{2\mu+2\lambda}tr(\sigma)I),$

where $I$ is the $2\times 2$ identity matrix and $\mu>0$, $\lambda\geq 0$ are the Lamé coefficients. In this case the elastic stress is $\sigma_{e}=2\mu\epsilon(u)+\lambda\nabla\cdot uI$. The poroelastic stress, which includes the effect of the fluid pressure $p$, is given as

$\displaystyle\sigma=\sigma_{e}-\alpha pI,$

(2.2)

where $0<\alpha\leq 1$ is the Biot-Willis constant.

Given a vector field $f$ representing the body forces and a source term $q$, the quasi-static Biot system that governs the fluid flow within the poroelastic media is given as follows:

$\displaystyle\begin{aligned} -\nabla\cdot\sigma&=f\quad&&{\rm in~{}}\Omega\times[0,T],\\ K^{-1}z+\nabla p&=0\quad&&{\rm in~{}}\Omega\times[0,T],\\ \frac{\partial}{\partial t}(c_{0}p+\alpha\nabla\cdot u)+\nabla\cdot z&=q\quad&&{\rm in~{}}\Omega\times[0,T],\end{aligned}$

(2.3)

where $z$ is the Darcy velocity, $c_{0}\geq 0$ is a mass storativity coefficient, $T>0$ is a finite time and $K$ is a symmetric and positive definite tensor representing the permeability of the porous media divided by the fluid viscosity, i.e., there exist two strictly positive real numbers $K_{1}$ and $K_{2}$ satisfying for almost every $x\in\Omega$ and all $\xi\in\mathbb{R}^{2}$ such that $|\xi|=1$

$$0<K_{1}\leq K(x)\xi\cdot\xi\leq K_{2}.$$

Further, we introduce the global anisotropy ratio

$\displaystyle\varrho_{B}=\frac{K_{2}}{K_{1}}.$

(2.4)

For the sake of simplicity, we assume that the system satisfies homogeneous Dirichlet boundary conditions for $u$ and $p$, i.e., $u=0,p=0$ on $\partial\Omega$ and the initial condition $p(x,0)=p^{0}$. Now we are ready to derive the weak formulation, we have from (2.1) and (2.2)

$\displaystyle\nabla\cdot u=tr(\epsilon(u))=tr(A\sigma_{e})=trA(\sigma+\alpha pI),$

which can be substituted in the third equation of (2.3) to get

$\displaystyle\partial_{t}(c_{0}p+\alpha trA(\sigma+\alpha pI))+\nabla\cdot z=q.$

In the weakly symmetric stress formulation, we allow for $\sigma$ to be non-symmetric and introduce the Lagrange multiplier $\gamma=\mbox{rot}(u)/2$, $\mbox{rot}(u)=-\frac{\partial u_{1}}{\partial y}+\frac{\partial u_{2}}{\partial x}$. The constitutive equation (2.1) can be rewritten as

$\displaystyle A(\sigma+\alpha pI)=\nabla u-\gamma\chi,$

where

$\displaystyle\chi=\left(\begin{array}[]{cc}0&-1\\ 1&0\\ \end{array}\right).$

By the preceding arguments, we can propose the following mixed weak formulation for the Biot system of poroelasticity: Find $(\sigma,u,\gamma,z,p)\in L^{2}(\Omega)^{2\times 2}\times H^{1}_{0}(\Omega)^{2}\times L^{2}(\Omega)\times L^{2}(\Omega)^{2}\times H^{1}_{0}(\Omega)$ such that

	$\displaystyle(A(\sigma+\alpha pI),\psi)-(\nabla u,\psi)+(\gamma,as(\psi))$	$\displaystyle=0\hskip 25.6073pt\forall\psi\in L^{2}(\Omega)^{2\times 2},$		(2.5)
	$\displaystyle(\sigma,\nabla v)$	$\displaystyle=(f,v)\quad\forall v\in H^{1}_{0}(\Omega)^{2},$		(2.6)
	$\displaystyle(as(\sigma),\eta)$	$\displaystyle=0\hskip 28.45274pt\forall\eta\in L^{2}(\Omega),$		(2.7)
	$\displaystyle(K^{-1}z,\xi)+(\nabla p,\xi)$	$\displaystyle=0\hskip 28.45274pt\forall\xi\in L^{2}(\Omega)^{2},$		(2.8)
	$\displaystyle c_{0}(\partial_{t}p,w)+\alpha(\partial_{t}A(\sigma+\alpha pI),wI)-(z,\nabla w)$	$\displaystyle=(q,w)\quad\forall w\in H^{1}_{0}(\Omega),$		(2.9)

where we use $(trA\psi,w)=(A\psi,wI)$ and $as(\psi)=-\psi_{12}+\psi_{21}$ for $\psi\in L^{2}(\Omega)^{2\times 2}$.

For the initial condition we can set $z^{0}=-K\nabla p^{0}$ which naturally satisfies (2.8). We can also determine $(\sigma^{0},u^{0},\gamma^{0})$ by solving the elasticity problem (2.5)-(2.7) with $p^{0}$ given as initial data. We refer to the initial data obtained by this procedure as compatible initial data.

Before closing this section, we also introduce some notations that will be used later. Let $D\subset\mathbb{R}^{2}$. By $(\cdot,\cdot)_{D}$, we denote the scalar product in $L^{2}(D):(p,q)_{D}:=\int_{D}p\;q\;dx$. When $D$ coincides with $\Omega$, the subscript $\Omega$ will be dropped. We use the same notation for the scalar product in $L^{2}(D)^{2}$ and in $L^{2}(D)^{2\times 2}$. More precisely, $(\xi,w)_{D}:=\sum_{i=1}^{2}(\xi^{i},w^{i})$ for $\xi,w\in L^{2}(D)^{2}$ and $(\psi,\zeta)_{D}:=\sum_{i=1}^{2}\sum_{j=1}^{2}(\psi^{i,j},\zeta^{i,j})_{D}$ for $\psi,\zeta\in L^{2}(D)^{2\times 2}$. The associated norm is denoted by $\|\cdot\|_{0,D}$. Given an integer $m\geq 0$ and $n\geq 1$, $W^{m,n}(D)$ and $W_{0}^{m,n}(D)$ denote the usual Sobolev space provided the norm and semi-norm $\|v\|_{W^{m,n}(D)}=\{\sum_{|\ell|\leq m}\|D^{\ell}v\|^{n}_{L^{n}(D)}\}^{1/n}$, $|v|_{W^{m,n}(D)}=\{\sum_{|\ell|=m}\|D^{\ell}v\|^{n}_{L^{n}(D)}\}^{1/n}$. If $n=2$ we usually write $H^{m}(D)=W^{m,2}(D)$ and $H_{0}^{m}(D)=W_{0}^{m,2}(D)$, $\|v\|_{H^{m}(D)}=\|v\|_{W^{m,2}(D)}$ and $|v|_{H^{m}(D)}=|v|_{W^{m,2}(D)}$. In addition, for an integer $l\geq 0$ we define spaces of vector valued functions such as $H^{l}(0,T;H^{m}(\Omega))$ and $L^{\infty}(0,T;H^{m}(\Omega))$ with the norms

$\displaystyle\|\psi\|_{H^{l}(0,T;H^{m}(\Omega))}=\sum_{i=0}^{l}\Big{(}\int_{0}^{T}\|d_{t}^{i}\psi(t)\|_{H^{m}(\Omega)}^{2}\;dt\Big{)}^{1/2},\quad\|\psi\|_{L^{\infty}(0,T;H^{m}(\Omega))}=\max_{0\leq t\leq T}\|\psi(t)\|_{H^{m}(\Omega)}.$

In the sequel, we use $C$ to denote a positive constant independent of the meshsize and of the anisotropy ratio $\varrho_{B}$ in (2.4), which may have different values at different occurrences. In addition, for any $\psi\in L^{2}(\Omega)^{2\times 2}$, $(A\psi,\psi)^{1/2}$ is a norm equivalent to $\|\psi\|_{0}$, which will be denoted by $\|A^{1/2}\psi\|_{0}$ throughout the paper.

In this section, we present the staggered DG discretization for the Biot system of poroelasticity. Specifically, a fully mixed formulation will be used, and the symmetry of stress is imposed weakly via the introduction of a lagrange multiplier.

To begin with, we introduce the construction of staggered meshes, which lays foundation for the construction of finite element spaces suitable for our setting. Let $\mathcal{T}_{u}$ be a shape-regular family of matching meshes covering $\Omega$ exactly, where hanging nodes are allowed. Note that the shape-regularity for general polygonal meshes can be described as follows: (1) Every element $M\in\mathcal{T}_{u}$ is star-shaped with respect to a ball of radius $\geq\rho_{B}h_{E}$, where $\rho_{B}$ is a positive constant and $h_{E}$ is the diameter of $M$; (2) For every element $M\in\mathcal{T}_{u}$ and every edge $e\subset\partial M$, it satisfies $h_{e}\geq\rho_{e}h_{M}$, where $\rho_{e}$ is a positive constant and $h_{e}$ represents the length of edge $e$. For an element $M\in\mathcal{T}_{u}$, we select an interior point $\nu$ and connect it to all the vertices of $M$. The resulting sub-triangulation is denoted as $\mathcal{T}_{h}$. We remark that $\nu$ is an arbitrary point interior to $M$ and for simplicity we choose it as center point. The union of all the edges in the primal partition $\mathcal{T}_{u}$ is denoted as $\mathcal{F}_{u}$. We use $\mathcal{F}_{u}^{0}$ to represent the subset of $\mathcal{F}_{u}$, that is the set of edges in $\mathcal{F}_{u}$ that does not lie on $\partial\Omega$. In addition, the edges generated during the subdivision process due to the connection of the interior point $\nu$ to all the vertices are called dual edges, which is denoted as $\mathcal{F}_{p}$. In addition, we define $\mathcal{F}=\mathcal{F}_{u}\cap\mathcal{F}_{p}$ and $\mathcal{F}=\mathcal{F}_{u}^{0}\cap\mathcal{F}_{p}$. For each triangle $\tau\in\mathcal{T}_{h}$, we let $h_{\tau}$ be the diameter of $\tau$ and $h=\max\{h_{\tau},\tau\in\mathcal{T}_{h}\}$. Next, we introduce the dual mesh. For each interior edge $e\in\mathcal{F}_{u}^{0}$, we use $D(e)$ to represent the dual mesh, which is the union of the two triangles in $\mathcal{T}_{h}$ sharing the common edge $e$. For each edge $e\in\mathcal{F}_{u}\backslash\mathcal{F}_{u}^{0}$, we use $D(e)$ to denote the triangle in $\mathcal{T}_{h}$ having the edge $e$, see Figure 1 for an illustration. For each edge $e$, we define a unit normal vector $n_{e}$ as follows: If $e$ is an interior edge, then $n_{e}$ is the unit normal vector of $e$ pointing towards the outside of $\Omega$. If $e$ is an interior edge, we then fix $n_{e}$ as one of the two possible unit normal vectors on $e$. When there is no ambiguity, we use $n$ instead of $n_{e}$ to simplify the notation. Let $k\geq 0$ be the order of approximation, for every $\tau\in\mathcal{T}_{h}$ and $e\in\mathcal{F}$, we define $P^{k}(\tau)$ and $P^{k}(e)$ as the spaces of polynomials of degree less than or equal to $k$ on $\tau$ and $e$, respectively.

For a scalar, vector, or tensor function $v$ that is double-valued on an interior edge $e\in\mathcal{F}^{0}$, its jump and average on $e$ are defined as

$\displaystyle[v]=v\mid_{\tau_{1}}-v\mid_{\tau_{2}},\quad\{v\}=\frac{1}{2}(v\mid_{\tau_{1}}+v\mid_{\tau_{2}}),$

where $\tau_{1}$ and $\tau_{2}$ are the two triangles belonging to $\mathcal{T}_{h}$ which share the common edge $e$. We set $[v]=v\mid_{\tau_{1}}$ and $\{v\}=v\mid_{\tau_{1}}$ on boundary edges.

Now we are ready to introduce the finite element spaces that will be used for the construction of our method. First, the locally $H^{1}(\Omega)-$conforming SDG space $S_{h}$ is defined by

$$S_{h}:=\{w\in L^{2}(\Omega):w\mid_{\tau}\in P^{k}(\tau),[w]\mid_{e}=0\;\forall e\in\mathcal{F}_{u};w\mid_{\partial\Omega}=0\},$$

which is equipped by the following norm

$\displaystyle\|w\|_{Z}^{2}=\sum_{\tau\in\mathcal{T}_{h}}\|\nabla w\|_{0,\tau}^{2}+\sum_{e\in\mathcal{F}_{p}}h_{e}^{-1}\|[w]\|_{0,e}^{2}.$

For a function $v=(v_{1},v_{2})\in[S_{h}]^{2}$ we define $\|v\|_{h}^{2}=\|v_{1}\|_{Z}^{2}+\|v_{2}\|_{Z}^{2}$.

Next, the locally $H(\mbox{div};\Omega)$-conforming SDG space $\Sigma_{h}$ is defined by

$$\Sigma_{h}:=\{\xi\in L^{2}(\Omega)^{2}:\xi\mid_{\tau}\in P^{k}(\tau)^{2}\;\forall\tau\in\mathcal{T}_{h},[\xi\cdot n]\mid_{e}=0\;\forall e\in\mathcal{F}_{p}\}.$$

Finally, the locally $H^{1}(\Omega)-$conforming finite element space is defined as

$$M_{h}:=\{\eta\in L^{2}(\Omega):\eta\mid_{\tau}\in P^{k}(\tau)\;\forall\tau\in\mathcal{T}_{h},[\eta]\mid_{e}=0\;\forall e\in\mathcal{F}_{p}\}.$$

We can derive the staggered DG discretization for the Biot system of poroelasticity by following [51, 53], which reads as follows: Find $(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h})\in[\Sigma_{h}]^{2}\times[S_{h}]^{2}\times M_{h}\times\Sigma_{h}\times S_{h}$ such that

	$\displaystyle(A(\sigma_{h}+\alpha p_{h}I),\psi)-B_{h}^{*}(u_{h},\psi)+(\gamma_{h},as(\psi))$	$\displaystyle=0\hskip 25.6073pt\forall\psi\in[\Sigma_{h}]^{2},$		(3.1)
	$\displaystyle B_{h}(\sigma_{h},v)$	$\displaystyle=(f,v)\quad\forall v\in[S_{h}]^{2},$		(3.2)
	$\displaystyle(as(\sigma_{h}),\eta)$	$\displaystyle=0\hskip 28.45274pt\forall\eta\in M_{h},$		(3.3)
	$\displaystyle(K^{-1}z_{h},\xi)+b_{h}^{*}(p_{h},\xi)$	$\displaystyle=0\hskip 29.87547pt\forall\xi\in\Sigma_{h},$		(3.4)
	$\displaystyle c_{0}(\partial_{t}p_{h},w)+\alpha(\partial_{t}A(\sigma_{h}+\alpha p_{h}I),wI)-b_{h}(z_{h},w)$	$\displaystyle=(q,w)\quad\forall w\in S_{h},$		(3.5)

where

	$\displaystyle B_{h}(\psi,v)$	$\displaystyle=-\sum_{e\in\mathcal{F}_{p}}(\psi n,[v])_{e}+(\psi,\nabla v)\quad\forall(\psi,v)\in[\Sigma_{h}]^{2}\times[S_{h}]^{2},$
	$\displaystyle B_{h}^{*}(v,\psi)$	$\displaystyle=\sum_{e\in\mathcal{F}_{u}^{0}}(v,[\psi n])_{e}-(v,\nabla\cdot\psi)\quad\forall(\psi,v)\in[\Sigma_{h}]^{2}\times[S_{h}]^{2}$

and

	$\displaystyle b_{h}(\xi,w)$	$\displaystyle=-\sum_{e\in\mathcal{F}_{p}}(\xi\cdot n,[w])_{e}+(\xi,\nabla w)\quad\forall(\xi,w)\in\Sigma_{h}\times S_{h},$
	$\displaystyle b_{h}^{*}(w,\xi)$	$\displaystyle=\sum_{e\in\mathcal{F}_{u}^{0}}(w,[\xi\cdot n])_{e}-(w,\nabla\cdot\xi)\quad\forall(\xi,w)\in\Sigma_{h}\times S_{h}.$

The bilinear forms defined above satisfy the following adjoint properties

	$\displaystyle b_{h}(\xi,w)$	$\displaystyle=b_{h}^{*}(w,\xi)\quad\forall(w,\xi)\in S_{h}\times\Sigma_{h},$
	$\displaystyle B_{h}(\psi,v)$	$\displaystyle=B_{h}^{*}(v,\psi)\quad\forall(\psi,v)\in[\Sigma_{h}]^{2}\times[S_{h}]^{2}.$

In addition, the following inf-sup condition holds; cf. [13].

$\displaystyle\|w\|_{Z}\leq C\sup_{\xi\in\Sigma_{h}}\frac{b_{h}^{*}(w,\xi)}{\|\xi\|_{0}}\quad\forall w\in S_{h}$

(3.6)

and

$\displaystyle\|v\|_{h}\leq C\sup_{\psi\in[\Sigma_{h}]^{2}}\frac{-B_{h}^{*}(v,\psi)}{\|\psi\|_{0}}\quad\forall v\in[S_{h}]^{2}.$

(3.7)

Now we prove the inf-sup condition, which is the key to proving the stability and convergence estimates of the proposed method. The methodology exploited here is different from that of [28]. Indeed the norm used here bypasses the construction of divergence free functions in proving the convergence estimates.

In order to prove the well-posedness of (3.1)-(3.5), we will make use of the following theorem; cf. [43].