A Stabilized Hybrid Mixed Finite Element Method for Poroelasticity

TL;DR

This paper introduces a stabilized hybrid mixed finite element method for Biot's poroelasticity model, improving stability and efficiency through face bubble functions, elimination techniques, and robust preconditioners, with confirmed numerical results.

Contribution

It develops a stabilized hybrid scheme with face bubble functions for Biot's model, ensuring stability, optimal error estimates, and efficient preconditioning, addressing previous stability issues.

Findings

Stable hybrid scheme with face bubble functions

Elimination of bubble functions and Darcy velocity

Numerical experiments confirm theoretical stability and efficiency

Abstract

In this work, we consider a hybrid mixed finite element method for Biot's model. The hybrid P1-RT0-P0 discretization of the displacement-pressure-Darcy's velocity system of Biot's model presented in \cite{C. Niu} is not uniformly stable with respect to the physical parameters, resulting in some issues in numerical simulations. To alleviate such problems, following \cite{V. Girault}, we stabilize the hybrid scheme with face bubble functions and show the well-posedness with respect to physical and discretization parameters, which provide optimal error estimates of the stabilized method. We introduce a perturbation of the bilinear form of the displacement which allows for the elimination of the bubble functions. Together with eliminating Darcy's velocity by hybridization, we obtain an eliminated system whose size is the same as the classical P1-RT0-P0 discretization. Based on the well-posedness of the eliminated system, we design block preconditioners that are parameter-robust. Numerical experiments are presented to confirm the theoretical results of the stabilized scheme as well as the block preconditioners.

Tables

Table 1. Table 1: Error and convergence rate of the hybrid mixed finite element method ( 5 ).

		$N=4$	$N=8$	$N=16$	$N=32$	$N=64$
K=1e-4	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0509	0.0270	0.0135	0.0068	0.0034
	$\Vert p-p_h\Vert$	0.1160	0.0535	0.0088	0.0015	0.0003
K=1e-6	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0570	0.0543	0.0314	0.0081	0.0034
	$\Vert p-p_h\Vert$	0.1587	0.3277	0.3199	0.0763	0.0099
K=1e-8	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0571	0.0571	0.0565	0.0478	0.0169
	$\Vert p-p_h\Vert$	0.1591	0.3553	0.7157	1.1509	0.6537
K=1e-10	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0571	0.0571	0.0571	0.0570	0.0550
	$\Vert p-p_h\Vert$	0.1588	0.3550	0.7271	1.4616	2.9182

Table 2. Table 2: Error and convergence rate of the stabilized hybrid mixed finite element method ( 11 ) with perturbation

		$N=4$	$N=8$	$N=16$	$N=32$	$N=64$
K=1e-4	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0369	0.0183	0.0093	0.0047	0.0024
	$\Vert p-p_h\Vert$	0.0511	0.0185	0.0034	0.0006	0.0001
K=1e-6	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0377	0.0189	0.0091	0.0045	0.0022
	$\Vert p-p_h\Vert$	0.0593	0.0346	0.0155	0.0062	0.0019
K=1e-8	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0377	0.0189	0.0092	0.0045	0.0023
	$\Vert p-p_h\Vert$	0.0594	0.0349	0.0162	0.0074	0.0035
K=1e-10	${\Vert 𝒖-{𝒖}_h\Vert }_a$	0.0377	0.0189	0.0092	0.0045	0.0023
	$\Vert p-p_h\Vert$	0.0594	0.0349	0.0162	0.0074	0.0035

Table 3. Table 3: Iteration counts of the block preconditioners for the eliminated system 𝒜 E superscript 𝒜 𝐸 \mathcal{A}^{E} discretized from the first example in Section 5.1 with varying physical parameters K 𝐾 K and ν 𝜈 \nu .

	$\nu =0$ and varying $K$
	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$	${10}^{-10}$	${10}^{-12}$
${ℬ}_D^E$	21	28	38	40	40	38
${ℬ}_U^E$	12	13	14	15	15	15
${ℬ}_L^E$	13	14	14	15	15	15
$\widehat{{ℬ}_D^E}$	29	38	44	46	44	43
$\widehat{{ℬ}_U^E}$	16	18	23	22	23	21
$\widehat{{ℬ}_L^E}$	20	22	21	22	20	20
	$K={10}^{-6}$ and varying $\nu$
	$0$	$0.1$	$0.2$	$0.4$	$0.45$	$0.49$
${ℬ}_D^E$	38	38	38	36	33	29
${ℬ}_U^E$	14	14	14	13	11	8
${ℬ}_L^E$	14	14	14	13	11	8
$\widehat{{ℬ}_D^E}$	44	44	44	45	44	40
$\widehat{{ℬ}_U^E}$	23	21	20	19	15	12
$\widehat{{ℬ}_L^E}$	21	21	20	16	15	12

Table 4. Table 4: Iteration counts of the block preconditioners for the eliminated system 𝒜 E superscript 𝒜 𝐸 \mathcal{A}^{E} discretized from the first example in Section  5.1 with varying discretization parameters h ℎ h and τ 𝜏 \tau .

	${ℬ}_D^E$
	$\frac{1}{4}$	$\frac{1}{8}$	$\frac{1}{16}$	$\frac{1}{32}$	$\frac{1}{64}$
$1$	35	39	40	38	35
$0.1$	32	39	40	40	39
$0.01$	34	39	40	40	39
$0.001$	33	37	38	38	38
$0.0001$	33	38	38	38	38

Table 5. Table 5: Iteration counts for the block preconditioners for the eliminated system 𝒜 E superscript 𝒜 𝐸 \mathcal{A}^{E} discretized from the cantilever bracket problem with varying physical parameters K 𝐾 K and ν 𝜈 \nu .

	$\nu =0.45$ and varying $K$
	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$	${10}^{-10}$	${10}^{-12}$
${ℬ}_D^E$	4	4	5	14	21	25
${ℬ}_U^E$	2	2	2	2	2	2
${ℬ}_L^E$	3	3	4	4	3	3
$\widehat{{ℬ}_D^E}$	5	6	10	29	36	38
$\widehat{{ℬ}_U^E}$	4	4	4	6	9	16
$\widehat{{ℬ}_L^E}$	5	5	6	8	8	11
	$K={10}^{-7}$ and varying $\nu$
	$0$	$0.1$	$0.2$	$0.4$	$0.45$	$0.49$
${ℬ}_D^E$	13	13	13	10	9	6
${ℬ}_U^E$	2	2	2	2	2	2
${ℬ}_L^E$	5	5	5	4	4	3
$\widehat{{ℬ}_D^E}$	23	23	23	20	18	11
$\widehat{{ℬ}_U^E}$	7	7	7	5	5	5
$\widehat{{ℬ}_L^E}$	9	9	9	7	7	7

Table 6. Table 6: Iteration counts for the block preconditioners for the eliminated system 𝒜 E superscript 𝒜 𝐸 \mathcal{A}^{E} discretized from the cantilever bracket problem with varying discretization parameters h ℎ h and τ 𝜏 \tau .

	${ℬ}_D^E$
	$\frac{1}{4}$	$\frac{1}{8}$	$\frac{1}{16}$	$\frac{1}{32}$	$\frac{1}{64}$
$1$	15	15	12	10	8
$0.1$	22	20	18	16	14
$0.01$	25	23	21	20	18
$0.001$	25	25	24	22	21
$0.0001$	26	26	26	25	23

Equations

$$\begin{array}[]{l}-\nabla\cdot\bm{\sigma}^{\prime}+\alpha\nabla p=\bm{f},\hskip 9.95845pt\text{in}\ \Omega,\\ \bm{\sigma}^{\prime}=2\mu\epsilon(\bm{u})+\lambda\ (\nabla\cdot\bm{u})\bm{I},\hskip 9.95845pt\text{in}\ \Omega,\\ \displaystyle{\frac{\partial}{\partial t}(\frac{1}{M}p+\alpha\nabla\cdot\bm{u})+\nabla\cdot\bm{w}=g},\hskip 9.95845pt\text{in}\ \Omega,\\ \bm{w}=-\kappa\nabla p,\hskip 9.95845pt\text{in}\ \Omega,\\ \end{array}$$

$$p=0,\text{ on}{\overline{\Gamma }}_t,-\kappa \nabla p\cdot \text{n}=0,\text{on}{\Gamma }_c,$$

$$\mathbf{u}=0,\text{on}{\overline{\Gamma }}_c,{\mathbf{\sigma }}^{\prime }\text{n}=0,\text{on}{\Gamma }_t,$$

$$p(\mathbf{x},0)=p_0,\mathbf{x}\in \Omega ,\mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0,\mathbf{x}\in \Omega .$$

$$\mathbf{V}=\{\mathbf{v}\in {\mathbf{H}}^1(\Omega ):\mathbf{v}{|}_{{\overline{\Gamma }}_c}=0\},$$

$$Q=L^2(\Omega ),$$

$$\mathbf{W}=\{\mathbf{w}\in \mathbf{H}(\text{div},\Omega ):\mathbf{w}\cdot \text{n}{|}_{{\Gamma }_c}=0\},$$

$$a(\mathbf{u},\mathbf{v})-\alpha (p,\nabla \cdot \mathbf{v})=(\mathbf{f},\mathbf{v}),\forall \mathbf{v}\in \mathbf{V},$$

$$\alpha (\nabla \cdot \frac{\partial \mathbf{u}}{\partial t},q)+\frac{1}{M}(\frac{\partial p}{\partial t},q)+(\nabla \cdot \mathbf{w},q)=(g,q),\forall q\in Q,$$

$$({\kappa }^{-1}\mathbf{w},\mathbf{r})-(p,\nabla \cdot \mathbf{r})=0,\forall \mathbf{r}\in \mathbf{W},$$

$$L^2({\mathcal{T}}_h)=\{p\in L^2(\Omega ):p{|}_T\in L^2(T),\forall T\in {\mathcal{T}}_h\},$$

$$H^1({\mathcal{T}}_h)=\{p\in L^2(\Omega ):p{|}_T\in H^1(T),\forall T\in {\mathcal{T}}_h\},$$

$$L^2(\partial {\mathcal{T}}_h)=\{p\in L^2(e),\forall e\in \partial {\mathcal{T}}_h\},$$

$$(p,q{)}_h=\sum _{T\in {\mathcal{T}}_h}(p,q{)}_T,(p,q{)}_T={\int }_{T}pqdx.$$

$$(\beta ,\rho {)}_{\partial {\mathcal{T}}_h}=\sum _{T\in {\mathcal{T}}_h}(\beta ,\rho {)}_{\partial T},(\beta ,\rho {)}_{\partial T}={\int }_{\partial T}\beta \rho ds,$$

$$\begin{array}[]{l}\bm{V}_{h,1}=\{\bm{v}_{h}\in\bm{V}:\bm{v}_{h}|_{T}\in[P_{1}(T)]^{d},\ \bm{v}_{h}|_{\partial\mathcal{T}_{B}}=\bm{0},\ \forall\ T\in\mathcal{T}_{h}\},\\ Q_{h}=\{q_{h}\in L^{2}(\mathcal{T}_{h}):q_{h}|_{T}\in P_{0}(T),\forall\ T\in\mathcal{T}_{h}\},\\ B_{h}=\{\rho_{h}\in L^{2}(\partial\mathcal{T}_{h}):\rho_{h}|_{e}\in P_{0}(e),\forall\ e\in\partial\mathcal{T}_{h}\},\\ \bm{W}_{h}=\{\bm{w}_{h}\in[H^{1}(\mathcal{T}_{h})]^{d}:\bm{w}_{h}|_{T}\in RT_{0}(T),\bm{w}_{h}\cdot\bm{n}_{e}|_{\partial\mathcal{T}_{B}}=0,\ \forall\ T\in\mathcal{T}_{h}\},\end{array}$$

$${\mathbf{V}}_h={\mathbf{V}}_{h,1}\oplus {\mathbf{V}}_b,{\mathbf{V}}_b=\text{span}\{{\mathbf{\Phi }}_e{\}}_{e\in \partial {\mathcal{T}}_{0,t}},$$

$$\begin{array}[]{l}\displaystyle{\Pi\bm{v}=\Pi_{1}\bm{v}+\sum_{e\in\partial\mathcal{T}_{0,t}}\nu_{e}\bm{\Phi}_{e}},\end{array}$$

$$\begin{array}[]{l}a(\bm{u}_{h}^{n},\bm{v}_{h})-\alpha(p_{h}^{n},\nabla\cdot\bm{v}_{h})=(\bm{f},\bm{v}_{h}),\ \forall\bm{v}_{h}\in\bm{V}_{h},\\ \displaystyle{\alpha(\nabla\cdot\bm{u}_{h}^{n},q_{h})+\frac{1}{M}(p_{h}^{n},q_{h})+\tau(\nabla\cdot\bm{w}_{h}^{n},q_{h})=\tau(\widetilde{g},q_{h})},\ \forall q_{h}\in Q_{h},\\ \tau(\bm{w}_{h}^{n}\cdot\bm{n}_{e},\rho_{h})_{\partial\mathcal{T}_{h}}=0,\ \forall\rho_{h}\in B_{h},\\ -\tau(p_{h}^{n},\nabla\cdot\bm{r}_{h})-\tau(\beta_{h}^{n},\bm{r}_{h}\cdot\bm{n}_{e})_{\partial\mathcal{T}_{h}}+\tau(\kappa^{-1}\bm{w}_{h}^{n},\bm{r}_{h})_{h}=0,\ \forall\bm{r}_{h}\in\bm{W}_{h},\end{array}$$

$$\begin{array}[]{l}a(\bm{u}_{h}^{n},\bm{v}_{h})-\alpha(p_{h}^{n},\nabla\cdot\bm{v}_{h})=(\bm{f},\bm{v}_{h}),\ \forall\bm{v}_{h}\in\bm{V}_{h,1},\\ \displaystyle{\alpha(\nabla\cdot\bm{u}_{h}^{n},q_{h})+\frac{1}{M}(p_{h}^{n},q_{h})+\tau(\nabla\cdot\bm{w}_{h}^{n},q_{h})=(\widetilde{g},q_{h})},\ \forall q_{h}\in Q_{h},\\ \tau(\bm{w}_{h}^{n}\cdot\bm{n}_{e},\rho_{h})_{\partial\mathcal{T}_{h}}=0,\ \forall\rho_{h}\in B_{h},\\ -\tau(p_{h}^{n},\nabla\cdot\bm{r}_{h})-\tau(\beta_{h}^{n},\bm{r}_{h}\cdot\bm{n}_{e})_{\partial\mathcal{T}_{h}}+\tau(\kappa^{-1}\bm{w}_{h}^{n},\bm{r}_{h})_{h}=0,\ \forall\bm{r}_{h}\in\bm{W}_{h}.\end{array}$$

$$\begin{array}[]{l}\mathcal{B}((\bm{u}_{h},p_{h},\beta_{h},\bm{w}_{h}),(\bm{v}_{h},q_{h},\rho_{h},\bm{r}_{h}))\\ :=a(\bm{u}_{h},\bm{v}_{h})-\alpha(p_{h},\nabla\cdot\bm{v}_{h})+\alpha(\nabla\cdot\bm{u}_{h},q_{h})+\displaystyle{\frac{1}{M}(p_{h},q_{h})}+\tau(\nabla\cdot\bm{w}_{h},q_{h})\\ \hskip 11.38109pt+\tau(\bm{w}_{h}\cdot\bm{n}_{e},\rho_{h})_{\partial\mathcal{T}_{h}}-\tau(p_{h},\nabla\cdot\bm{r}_{h})-\tau(\beta_{h},\bm{r}_{h}\cdot\bm{n}_{e})_{\partial\mathcal{T}_{h}}+\tau(\kappa^{-1}\bm{w}_{h},\bm{r}_{h})_{h},\end{array}$$

$$\begin{array}[]{l}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\bm{u}_{h},p_{h},\beta_{h},\bm{w}_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}:=\displaystyle{\biggl{[}\|\bm{u}_{h}\|_{a}^{2}+\delta\|p_{h}\|^{2}+\frac{1}{h}|\beta_{h}|_{\partial\mathcal{T}_{h}}^{2}+\tau\|\bm{w}_{h}\|_{h,\kappa^{-1}}^{2}}\\ \hskip 105.27519pt+\displaystyle{\frac{\tau^{2}}{\delta}\|\nabla\cdot\bm{w}_{h}\|_{h}^{2}+\tau^{2}h\|\bm{w}_{h}\cdot\bm{n}_{e}\|_{h}^{2}\biggr{]}^{1/2}},\\ \end{array}$$

$$\begin{array}[]{l}\displaystyle{\sup_{(\bm{v}_{h},q_{h},\rho_{h},\bm{r}_{h})\in S_{h}}}\frac{\mathcal{B}((\bm{u}_{h},p_{h},\beta_{h},\bm{w}_{h}),(\bm{v}_{h},q_{h},\rho_{h},\bm{r}_{h}))}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\bm{v}_{h},q_{h},\rho_{h},\bm{r}_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}}\geq\gamma_{\ast}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\bm{u}_{h},p_{h},\beta_{h},\bm{w}_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},\end{array}$$

$$\begin{array}[]{l}|\mathcal{B}((\bm{u}_{h},p_{h},\beta_{h},\bm{w}_{h}),(\bm{v}_{h},q_{h},\rho_{h},\bm{r}_{h}))|\leq\gamma^{\ast}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\bm{u}_{h},p_{h},\beta_{h},\bm{w}_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\bm{v}_{h},q_{h},\rho_{h},\bm{r}_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}.\end{array}$$

$$\begin{array}[]{l}\displaystyle{(\frac{2\mu}{d}+\lambda)\|\nabla\cdot\bm{v}\|^{2}\leq\|\bm{v}\|_{a}^{2}\Rightarrow\|\nabla\cdot\bm{v}\|\leq\frac{1}{\zeta}\|\bm{v}\|_{a}}.\\ \end{array}$$

$$\frac{1}{\tau }\Vert \nabla u{\Vert }_{\kappa }^2+\frac{{\theta }^2}{{\tau }^2}\Vert u{\Vert }^2\le C\Vert {\mathbf{z}}_{\mathbf{n}}{\Vert }_{H^{-1/2}(T)}^2,$$

$$\frac{1}{\tau }\Vert {\kappa }^{-1}\mathbf{z}{\Vert }_{\kappa }^2+\frac{{\theta }^2}{{\tau }^2}\Vert \frac{\tau }{{\theta }^2}\nabla \cdot \mathbf{z}{\Vert }^2\le C\Vert \mathbf{z}\cdot {\mathbf{n}}_e{\Vert }_{H^{-1/2}(T)}^2.$$

$$\frac{1}{\tau^{2}}\biggl{(}\tau\|\bm{z}\|^{2}_{h,\kappa^{-1}}+\frac{\tau^{2}}{\theta^{2}}\|\nabla\cdot\bm{z}\|^{2}\biggr{)}\leq C\|\bm{z}\cdot\bm{n}_{e}\|^{2}_{H^{-1/2}(T)}.$$

$$\begin{array}[]{l}\displaystyle{\tau\|\bm{z}_{h}\|^{2}_{h,\kappa^{-1}}+\frac{\tau^{2}}{\theta^{2}}\|\nabla\cdot\bm{z}_{h}\|_{h}^{2}}\\ \displaystyle{=\sum_{T\in\mathcal{T}_{h}}\biggl{(}\tau\|\Pi^{\text{div}}_{T}\bm{z}\|^{2}_{T,\kappa^{-1}}+\frac{\tau^{2}}{\theta^{2}}\|\Pi^{L^{2}}_{T}\nabla\cdot\bm{z}\|_{T}^{2}\biggr{)}}\\ \displaystyle{\leq\sum_{T\in\mathcal{T}_{h}}\biggl{(}\tau\|\bm{z}\|^{2}_{T,\kappa^{-1}}+\frac{\tau^{2}}{\theta^{2}}\|\nabla\cdot\bm{z}\|_{T}^{2}\biggr{)}}\\ \displaystyle{\leq\sum_{T\in\mathcal{T}_{h}}\frac{C}{h}\|\beta_{h}\|_{L^{2}(T)}^{2}}\\ \displaystyle{\leq\frac{C}{h}|\beta_{h}|_{\partial\mathcal{T}_{h}}^{2}},\end{array}$$

$$\mathcal{B}(({\mathbf{u}}_h,p_h,{\beta }_h,{\mathbf{w}}_h),({\mathbf{v}}_h,q_h,{\rho }_h,{\mathbf{r}}_h))$$

$$=\Vert {\mathbf{u}}_h{\Vert }_a^2-\theta a({\mathbf{u}}_h,\mathbf{h})+\theta \alpha (p_h,\nabla \cdot \mathbf{h})+{\theta }_1\alpha \tau (\nabla \cdot {\mathbf{u}}_h,\nabla \cdot {\mathbf{w}}_h)+\frac{1}{M}\Vert p_h{\Vert }^2$$

$$+{\theta }_1\frac{\tau }{M}(p_h,\nabla \cdot {\mathbf{w}}_h)+{\theta }_1{\tau }^2\Vert \nabla \cdot {\mathbf{w}}_h{\Vert }_h^2+{\theta }_2{\tau }^{2}h\Vert {\mathbf{w}}_h\cdot {\mathbf{n}}_e{\Vert }_h^2+{\theta }_3\tau (p_h,\nabla \cdot {\mathbf{z}}_h)$$

$$+{\theta }_3\tau ({\beta }_h,{\mathbf{z}}_h\cdot {\mathbf{n}}_e{)}_{\partial {\mathcal{T}}_h}+\tau \Vert {\mathbf{w}}_h{\Vert }_{h,{\kappa }^{-1}}^2-{\theta }_3\tau ({\kappa }^{-1}{\mathbf{w}}_h,{\mathbf{z}}_h{)}_h$$

$$\ge \Vert {\mathbf{u}}_h{\Vert }_a^2-\frac{1}{2}\Vert {\mathbf{u}}_h{\Vert }_a^2-\frac{{\theta }^2}{2}\Vert \mathbf{h}{\Vert }_a^2+\theta \frac{\alpha {\gamma }_B}{\zeta }\Vert p_h{\Vert }^2-\frac{{\theta }_1{\alpha }^2}{2}\Vert \nabla \cdot {\mathbf{u}}_h{\Vert }^2$$

$$-\frac{{\theta }_1{\tau }^2}{2}\Vert \nabla \cdot {\mathbf{w}}_h{\Vert }_h^2+\frac{1}{M}\Vert p_h{\Vert }^2-\frac{3{\theta }_1}{2M^2}\Vert p_h{\Vert }^2-\frac{{\theta }_1{\tau }^2}{6}\Vert \nabla \cdot {\mathbf{w}}_h{\Vert }_h^2$$

$$+{\theta }_1{\tau }^2\Vert \nabla \cdot {\mathbf{w}}_h{\Vert }_h^2+{\theta }_2{\tau }^{2}h\Vert {\mathbf{w}}_h\cdot {\mathbf{n}}_e{\Vert }_h^2-\frac{{\theta }^2}{2}\Vert p_h{\Vert }^2-\frac{{\theta }_3^2{\tau }^2}{2{\theta }^2}\Vert \nabla \cdot {\mathbf{z}}_h{\Vert }_h^2$$

Full text

††journal: Computational Geosciences

1 Introduction

Poroelasticity theory simultaneously describes the interaction between the deformation and fluid flow in a fluid-saturated porous medium. Nowadays, poroelasticity theory is well developed, dating back to the pioneering work of Terzaghi K. Terzaghi , who studied a one-dimensional consolidation problem. A general three-dimensional mathematical model was established by Biot M.A. Biot41 , who later extended it to anisotropic and nonlinear materials M.A. Biot55 ; M.A. Biot3 . Results on the existence and uniqueness of the solution for these models have been investigated by Showalter in R. Showalter and Ženíšek in A. Zenisek . The well-posedness for nonlinear poroelastic models is considered in R.Z. Dautov . Lewis and Schrefler R.W. Lewis1 presented a complete work on theoretical, practical and numerical aspects of geomechanical problems. The analysis and numerical simulation of Biot’s model has became increasingly popular due to the wide range of applications such biomechanics, petroleum engineering, and environmental engineering.

Numerically solving the poroelasticity model is a challenging task due to its complex coupled nature. There is extensive literature on numerical methods for poroelasticity. For the two-field displacement-pressure formulation, Lewis and Schrefler R.W. Lewis2 analyzed one and two-dimensional problems in consolidation and the problem of subsidence in Venice by using a continuous Galerkin (CG) method. Furthermore, Liu R. Liu implemented a discontinuous Galerkin (DG) method for both displacement and pressure. A weak Galerkin (WG) method for the two-field formulations was considered in X. Hu1 ; Y. Chen . Since Phillips and Wheeler P.J. Phillips1 ; P.J. Phillips2 ; P.J. Phillips3 developed a method that couples a CG method for the displacement with a mixed finite element method for the pressure and velocity, more and more researchers have begun to focus on the three-field formulation, which has displacement, Darcy’s velocity, and pressure as unknowns. Yi S.Y. Yi1 developed a nonconforming finite element method, which has been used to overcome nonphysical oscillations in the pressure variable. Hu et.al. X. Hu2 also presented a nonconforming method based on the Crouzeix-Raviart elements for the displacements, lowest order Raviart-Thomas-Nédélec elements for Darcy’s velocity P.A. Raviart ; J.C. Nedelec , and piecewise constant elements for the pressure. A DG scheme was proposed in HongKraus2018 for the three-field formulation. A WG scheme which couples a WG method for the displacement with a standard mixed finite element method for the pressure and Darcy’s velocity, was proposed in M. Sun to avoid the locking. Other numerical schemes, such as least squares mixed finite element methods, were proposed in Korsawe ; Tchonkova for the four-field formulation (displacement, stress, Darcy’s velocity and pressure). And a new mixed finite element method for both the flow subproblem and the mechanical subproblem was introduced by Yi S.Y. Yi14 , and then S.Y. Yi17 presented the iteratively coupled solution strategies for the four-field formulation. More recently, an element-based finite volume formulation was proposed in Herminio T.Honorio18 to avoid pressure instabilities in poromechanics, and I. Sokolova19 presented a multiscale finite volume method.

For the three-field formulation, one standard approximation scheme, piecewise constant approximation for the pressure and piecewise linear approximation for the displacements and fluid flux may cause instability and non-physical oscillations, a local pressure jump stabilization term was introduced in L. Berger15 to ensure stability. One of the most frequently considered low-order schemes for the three-field formulation is P1-RT0-P0, i.e., piecewise linear elements for the displacement, Raviart-Thomas-Nédélec elements for the Darcy’s velocity, and piecewise constant elements for the pressure. However, the triple P1-RT0-P0 does not satisfy the Biot-Stokes stability condition uniformly with respect to the discretization and physical parameters of the problem V. Girault ; C. Rodrigo ; HongKraus2018 . For example, when the permeability is small with respect to the mesh size, volumetric locking may occur C. Rodrigo . Due to the same reason, the hybridized P1-RT0-P0 scheme presented in C. Niu has the same stability issue. Recently, a stabilization strategy based on the macro-element theory and the local pressure jump approach for the hybridized P1-RT0-P0 formulation has been proposed in M. Frigo20 , and an efficient solver was introduced to improve computational efficiency.

In this paper, we apply the stabilization technique introduced in C. Rodrigo and enriching the piecewise linear finite element space with edge/face bubble functions for displacement to stabilize the hybrid P1-RT0-P0 scheme developed in C. Niu . As shown in F. Brezzi ; D.N. Arnold ; H. Egger , the hybridization technique, which was proposed in T.H.H. Pian , is characterized by the removal of the continuity of the normal component of Darcy’s velocity along each element interface. Therefore, one particular advantage is that the degrees of freedom for Darcy’s velocity can be eliminated by static condensation, which reduces the computational cost. For our proposed stabilized hybrid scheme, a perturbation of the bilinear form allows for the elimination of the bubble functions. Therefore, both the unknowns of Darcy’s velocity and the bubble functions can be eliminated in our case, which results in a system with the same number of degrees of freedom as the standard P1-RT0-P0 discretization. The eliminated system is proved to be well-posed with respect to the discretization as well as the physical parameters, and thus, based on the framework developed in D. Loghin ; K.A. Mardal ; James , we develop robust block preconditioners to solve the resulting linear systems of equations efficiently. Moreover, due to the introduction of hybridization and eliminations of the fluid flux, in the implementation of the block preconditioners, we avoid solving complicated $\operatorname{grad}\operatorname{div}$-type subproblem, which requires special solvers James , and only need to solve standard $\operatorname{div}\operatorname{grad}$-type subproblem, which can be efficiently handled by standard multigrid methods.

The rest of the paper is organized as follows. In Section 2, we briefly recall the Biot’s model, and mixed variational formulation of the three-field formulation, and then devote to a stabilization technique with edge bubble functions and the perturbation of the bilinear form. The well-posedness of the resulting scheme, as well as the corresponding error analysis are also provided. In Section 3 we discuss the elimination of bubbles and Darcy’s velocity, and show the well-posedness of the eliminated system. Both the block diagonal and triangular preconditioners are presented in Section 4. Numerical results are presented in Section 5 to validate the accuracy and efficiency of the stabilization method and demonstrate the robustness and effectiveness of the preconditioners. Finally, our conclusion are presented in Section 7.

2 Stabilized Hybrid Mixed Finite Element Method

In this section, we first recall the Biot’s consolidation model M.A. Biot41 in a domain $\Omega\subset R^{d}$ $(d=2,3)$, and then analyze the stabilized hybrid mixed finite element method.

2.1 The Biot’s Model

In this subsection, we review Biot’s consolidation model and present a three-field formulation of which the primary variables are displacement $\bm{u}$, pressure $p$ and Darcy’s velocity $\bm{w}$.

[TABLE]

where $\kappa$ stands for the permeability tensor, $M$ is the Biot modulus, and $\alpha$ is the Biot-Willis constant. $\lambda$ and $\mu$ are Lamé constants. Here, the effective stress tensor and the strain tensor are denoted by $\bm{\sigma}^{\prime}$ and $\epsilon(\bm{u})=\frac{1}{2}(\nabla\bm{u}+\nabla\bm{u}^{\top})$, and $\bm{I}$ is the identity tensor. Let $g$ be the volumetric source/sink term, and $\bm{f}$ represents the body force. The boundary conditions are

[TABLE]

[TABLE]

where n represents the outward unit normal vector to the boundary, $\overline{\Gamma}_{t}\cup\overline{\Gamma}_{c}=\overline{\partial\Omega}$ with $\Gamma_{c}$ and $\Gamma_{t}$ being open (with respect to $\partial\Omega$) subsets of $\partial\Omega$ with nonzero measure. The initial condition at $t=0$ is given by,

[TABLE]

Next, we give the function spaces which are used in the variational form:

[TABLE]

where $\bm{H}^{1}(\Omega)$ is the space of square integrable vector-valued functions whose first derivatives are also square integrable, and $\bm{H}(\textup{div},\Omega)$ contains the square integrable vector-valued functions with square integrable divergence.

Finally, we give the mixed variational formulation: For all $t\in(0,t_{max}]$, find $(\bm{u},p,\bm{w})\in\bm{V}\times Q\times\bm{W}$, such that,

[TABLE]

where $a(\bm{u},\bm{v})=2\mu\int_{\Omega}\epsilon(\bm{u}):\epsilon(\bm{v})+\lambda\int_{\Omega}\nabla\cdot\bm{u}\nabla\cdot\bm{v}.$ Well-posedness of the continuous problem for this three field formulation was established by Lipnikov K. Lipnikov .

2.2 Stabilized Hybrid Mixed Finite Element Method

Let $\{\mathcal{T}_{h}\}$ be a quasi-uniform regular partition of $\Omega$, the element $T\in{\mathcal{T}_{h}}$ is triangular (d=2) or tetrahedral (d=3). Denote $\partial\mathcal{T}_{I}$ the set of interior edges, i.e., the set of common edges of $T^{+}$ $\cap$ $T^{-}$ for all neighboring $T^{+},T^{-}\in{\mathcal{T}_{h}}$. Correspondingly, let $\partial\mathcal{T}_{B}$ denote the set of boundary edges, i.e., the set of the common edges of $T\cap\partial\Omega$ for all $T\in{\mathcal{T}_{h}}$, and denote $\partial{\mathcal{T}_{h}}=\partial\mathcal{T}_{I}\cup\partial\mathcal{T}_{B}$. For every edge $e\in\partial\mathcal{T}_{h}$, we associate it with a unit outer normal $\bm{n}_{e}$. For the boundary faces $e\in\partial\mathcal{T}_{B}$, we set $\bm{n}_{e}=\pm\bm{n}_{e,T}$, where $\bm{n}_{e,T}$ is the outward (with respect to $T\in\mathcal{T}_{h}$) unit normal vector to the edge $e$. For the interior edges $e\in\partial\mathcal{T}_{I}$, the direction of $\bm{n}_{e}$ is fixed, while the particular direction of $\bm{n}_{e}$ is not important. Hence, we define $\bm{n}_{e}=\bm{n}_{e,T^{+}}=-\bm{n}_{e,T^{-}}$.

Now we introduce the piecewise Sobolev spaces:

[TABLE]

where $L^{2}(T)$ denotes the set of square integrable functions on element $T$, $H^{1}(T)$ is the space of square integrable functions on element $T$ whose first derivatives are also square integrable, and $L^{2}(e)$ the set of square integrable functions on edge $e$. The norm endowed with $L^{2}(\mathcal{T}_{h})$ is defined as $\|p\|_{h}:=\sqrt{(p,p)_{h}}$, and the corresponding inner product is defined as

[TABLE]

For $\beta,\rho\in L^{2}(\partial\mathcal{T}_{h})$, we define the inner product on element interfaces as follows,

[TABLE]

and the corresponding norms are denoted by $|\rho|_{\partial\mathcal{T}_{h}}:=\sqrt{(\rho,\rho)_{\partial\mathcal{T}_{h}}}$. For the vector-valued functions, the norm and inner product are defined in a similar manner with corresponding modifications.

Following the standard idea of hybrid mixed finite element method, we use completely discontinuous piecewise polynomial functions for the Darcy’s velocity and ensure the continuity of the normal fluxes over element interfaces by adding appropriate constraints and introducing a Lagrange multiplier. We introduce the discrete finite element spaces for the displacement, pressure, Lagrange multiplier and Darcy’s velocity,

[TABLE]

where $P_{0}(T)$ denotes the set of piecewise constant functions restricted to $T$ for each $T\in\mathcal{T}_{h}$ and $RT_{0}(T)=[P_{0}(T)]^{d}\oplus\textup{span}(\bm{x}P_{0}(T))$ denotes the standard lowest order Raviart-Thomas-Nédélec space.

For uniformly positive definite permeability tensor $\kappa$, the choice of spaces $\bm{V}_{h,1}$ has been successfully employed for numerical simulations of Biot’s consolidation model (see P.J. Phillips2 ; K. Lipnikov ). However, the heuristic considerations that expose some of the issues with this discretization are observed when $\kappa\rightarrow\bm{0}$. In such case, the discrete problem approaches a P1-P0 discretization of the Stokes’ equation. As it is well known, the finite element pair, $\bm{V}_{h,1}\times Q_{h}$, does not satisfy the inf-sup condition and is unstable for the Stokes’ problem.

Following V. Girault ; C. Rodrigo , we design a stabilized hybrid mixed finite element method by enriching the piecewise linear finite element space $\bm{V}_{h,1}$, with edge bubble functions in 2D or face bubble functions in 3D. Now we introduce the stabilized finite element space $\bm{V}_{h}$ as

[TABLE]

where $\bm{\Phi}_{e}=\varphi_{e}\bm{n}_{e}$, and $\varphi_{e}|_{T^{\pm}}=\varphi_{e,T^{\pm}}=\displaystyle{\prod_{l=1,l\neq j^{\pm}}^{d+1}\lambda_{l,T^{\pm}}}$, for every face $e\in\partial\mathcal{T}_{h}$, $e=T^{+}\cap T^{-}$. Here, $\lambda_{l,T^{\pm}},l=1,...,(d+1)$ are barycentric coordinates on $T^{\pm}$ and $j^{\pm}$ is the vertex opposite to face $e$ in $T^{\pm}$. Note that $\bm{\Phi}_{e}\in\bm{V}$ is a continuous piecewise polynomial function of degree $d$.

The degree of freedom associated with $\bm{V}_{h}$ are the values at the vertices of $\mathcal{T}_{h}$ and the total flux through $e\in\partial\mathcal{T}_{0,t}$. The canonical interpolant $\Pi:C(\overline{\Omega})\mapsto\bm{V}_{h}$ is defined as:

[TABLE]

where $\Pi_{1}:C(\overline{\Omega})\mapsto\bm{V}_{h,1}$ is the standard piecewise linear interpolant, and $\displaystyle{\nu_{e}=\frac{1}{|e|}\int_{e}(I-\Pi_{1})\bm{v}}$.

Finally, by using backward Euler time discretization with constant time step size $\tau$, the stabilized hybrid mixed finite element method is, for all $t_{n}=n\tau,n=1,2,...,$ find $(\bm{u}_{h}^{n},p_{h}^{n},\beta_{h}^{n},\bm{w}_{h}^{n})\in\bm{V}_{h}\times Q_{h}\times B_{h}\times\bm{W}_{h}$, such that,

[TABLE]

where $(\widetilde{g},q_{h})=\tau(g,q_{h})+\frac{1}{M}(p_{h}^{n-1},q_{h})+\alpha(\nabla\cdot\bm{u}_{h}^{n-1},q_{h})$, $(\bm{u}_{h}^{n},p_{h}^{n},\beta_{h}^{n},\bm{w}_{h}^{n})\approx(\bm{u}(\cdot,t_{n}),p(\cdot,t_{n}),\beta(\cdot,t_{n}),\bm{w}(\cdot,t_{n}))$, $t_{n}=n\tau$, $n=1,2,...$, and the initial conditions are $\bm{u}_{h}^{1}=\bm{u}_{0}$, $p_{h}^{1}=p_{0}$.

Remark 2.1**.**

When the space $\bm{V}_{h,1}$ is used for displacement, (4) is just the hybrid mixed finite element method presented in C. Niu : For all $t_{n}=n\tau,n=1,2,...,$ find $(\bm{u}_{h}^{n},p_{h}^{n},\beta_{h}^{n},\bm{w}_{h}^{n})\in\bm{V}_{h,1}\times Q_{h}\times B_{h}\times\bm{W}_{h}$, such that,

[TABLE]

2.3 Well-posedness

To show that the discrete problem (4) is well-posed, we introduce a bilinear form on the space $S_{h}:=\bm{V}_{h}\times Q_{h}\times B_{h}\times\bm{W}_{h}$,

[TABLE]

and a norm,

[TABLE]

where $\displaystyle{\zeta=\sqrt{\lambda+2\mu/d}}$, $\displaystyle{\delta=\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}}$, and $\|\bm{w}\|_{\mathcal{T}_{h},\kappa^{-1}}=(\kappa^{-1}\bm{w},\bm{w})^{1/2}$. Then, we have the following theorem to show the bilinear form $\mathcal{B}(\cdot,\cdot)$ is inf-sup stable and continuous, which shows that the hybrid mixed finite element method is well-posed at each time step.

Theorem 2.2**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), the following inf-sup condition holds

[TABLE]

and, $\mathcal{B}(\cdot,\cdot)$ is continuous with respect to ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$, i.e.,

[TABLE]

Here, the constants $\gamma_{\ast}$, $\gamma^{\ast}>0$ are independent of mesh size $h$, time step size $\tau$, and the physical parameters.

Proof.

First, for $d=2,3$, $2\mu(\epsilon(\bm{u}):\epsilon(\bm{v}))\leq a(\bm{u},\bm{v})\leq(2\mu+d\lambda)(\epsilon(\bm{u}):\epsilon(\bm{v}))$, then $(\nabla\cdot\bm{v},\nabla\cdot\bm{v})\leq d(\epsilon(\bm{u}):\epsilon(\bm{v}))$, and

[TABLE]

According to the definition of $\zeta$, we have $\|\bm{v}\|_{a}\leq\sqrt{d}\zeta\|\bm{v}\|_{1}$. Due to our choice of finite element spaces, for a given $p\in Q_{h}$, there exists an $\bm{h}\in\bm{V}_{h}$ such that $\displaystyle{(p,\nabla\cdot\bm{h})\geq\frac{\gamma_{B}}{\zeta}\|p\|^{2}}$ and $\|\bm{h}\|_{a}=\|p\|$ V. Girault .

Next, following J. Schoberl , let $\bm{z}_{\bm{n}}\in H^{-1/2}(T)$, we solve the weak form of the scalar equation $\displaystyle{-\frac{1}{\tau}\nabla\cdot\kappa\nabla u+\frac{\theta^{2}}{\tau^{2}}u=0}$ with boundary conditions $\displaystyle{\kappa\frac{\partial u}{\partial n}}$ $\displaystyle{=\bm{z}_{\bm{n}}}$, and there exists a unique solution $u\in H^{1}(T)$ such that,

[TABLE]

where, $\|\nabla u\|_{\kappa}=(\kappa\nabla u,\nabla u)^{1/2}$. Set $\bm{z}=\kappa\nabla u$, then $\nabla\cdot\bm{z}=\kappa\triangle u=\displaystyle{\frac{\theta^{2}}{\tau}u}$ and $\bm{z}_{\bm{n}}=\bm{z}\cdot\bm{n}_{e}$, thus,

[TABLE]

Immediately, we have,

[TABLE]

Introduce the interpolation operator $\Pi^{\text{div}}_{T}$ F. Brezzi : $H(\text{div},T)\rightarrow RT_{0}(T)$ as $\bm{z}_{h}|_{T}=\Pi^{\text{div}}_{T}\bm{z}|_{T}$. Let $\Pi^{L^{2}}_{T}$ be the $L^{2}$ projection on $P_{0}(T)$, we have $\nabla\cdot\bm{z}_{h}|_{T}=\nabla\cdot\Pi^{\text{div}}_{T}\bm{z}|_{T}=\Pi^{L^{2}}_{T}\nabla\cdot\bm{z}|_{T}$. Let $\bm{z}_{h}\cdot\bm{n}_{e}=\displaystyle{\frac{1}{\tau h}\beta_{h}\in B_{h}}$, then we have

[TABLE]

where the properties of the operators $\Pi^{\text{div}}_{T}$ and $\Pi^{L^{2}}_{T}$ R. Hiptmair derives the first inequality, the second inequality comes from (6) directly and the third inequality is derived from the property of space $H^{-1/2}$ J. Schoberl .

Let $\bm{v}_{h}=\bm{u}_{h}-\theta\bm{h}$, $q_{h}=p_{h}+\theta_{1}\tau\nabla\cdot\bm{w}_{h}$, $\rho_{h}=\beta_{h}+\theta_{2}\tau h\bm{w}_{h}\cdot\bm{n}_{e}$, and $\bm{r}_{h}=\bm{w}_{h}-\theta_{3}\bm{z}_{h}$, then by the Cauchy-Schwarz inequality and (7), we have,

[TABLE]

where $\displaystyle{\gamma_{1}=\min\{\frac{\gamma_{B}^{2}}{4},\frac{1}{2C},\frac{1}{6}\}}$, and we choose $\displaystyle{\theta=\frac{\alpha\gamma_{B}}{2\zeta}}$, $\displaystyle{\theta_{1}=\frac{1}{2\delta}}$, $\theta_{2}=1$, and $\displaystyle{\theta_{3}=\frac{1}{C}}$. Then by the Cauchy-Schwarz inequality, we have,

[TABLE]

where $\displaystyle{\gamma_{2}^{2}=\max\{2+\frac{\gamma_{B}^{2}}{2},\frac{2}{C}+\frac{\gamma_{B}^{2}}{2C}+\frac{2}{C^{2}}+2,4\}}$. Therefore, the inf-sup condition holds with $\displaystyle{\gamma_{\ast}=\frac{\gamma_{1}}{\gamma_{2}}}$.

For the boundedness of $\mathcal{B}(\cdot,\cdot)$, using the Cauchy-Schwarz inequality for all the terms on right hand side and $\displaystyle{\|\nabla\cdot\bm{v}\|\leq\frac{1}{\zeta}\|\bm{v}\|_{a}}$, it is easy to show that $\mathcal{B}(\cdot,\cdot)$ is continuous with respect to norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$. ∎

2.4 Error Estimates

To derive the error analysis of the fully discrete scheme, following the error analysis in F. Brezzi ; C. Rodrigo ; V. Thomee , we first define the following elliptic projections $\overline{\bm{u}}_{h}\in\bm{V}_{h}$, $\overline{p}_{h}\in Q_{h}$, $\overline{\beta}_{h}\in B_{h}$ and $\overline{\bm{w}}_{h}\in\bm{W}_{h}$ for $t>0$,

[TABLE]

We split the errors as follows,

[TABLE]

And the following lemma shows the error estimates for the elliptic projection defined in (8).

Lemma 2.3**.**

The following error estimates for the elliptic projections defined in (8) hold for $t>0$,

[TABLE]

Proof.

The proof is same as the proof of Lemma 4.4 C. Rodrigo , based on the error analysis of the hybrid mixed formulation of Poisson problems. ∎

Similarly, we define the elliptic projections, $\overline{\partial_{t}\bm{u}_{h}}$, $\overline{\partial_{t}p_{h}}$ and $\overline{\partial_{t}\bm{w}_{h}}$ of $\overline{\partial}_{t}\bm{u}_{h}$, $\overline{\partial}_{t}p_{h}$ and $\overline{\partial}_{t}\bm{w}_{h}$ respectively. This gives similar estimates as above for $\partial_{t}\rho_{\bm{u}}$, $\partial_{t}\rho_{p}$ and $\partial_{t}\rho_{\bm{w}}$, where on the right-hand side of the inequalities we use norms of $\overline{\partial}_{t}\bm{u}_{h}$, $\overline{\partial}_{t}p_{h}$ and $\overline{\partial}_{t}\bm{w}_{h}$ instead of the norms of $\bm{u}_{h}$, $p_{h}$ and $\bm{w}_{h}$ respectively. Then, we estimate the errors using the following norm,

[TABLE]

where $\|\bm{w}_{h}\|^{2}_{h,\kappa^{-1}}:=(\kappa^{-1}\bm{w}_{h},\bm{w}_{h})_{h}.$

Theorem 2.4**.**

Let $\bm{u}$, $p$ and $\bm{w}$ be the solution of (1), and $\bm{u}_{h}$, $p_{h}$ and $\bm{w}_{h}$ be the solution of (4). If they satisfy the following regularity assumptions,

[TABLE]

then,

[TABLE]

Proof.

See Appendix. ∎

2.5 Perturbation of the Bilinear Form $a(\cdot,\cdot)$

In general, using edge/face bubbles leads to a prohibitively large linear system. To resolve this, following C. Rodrigo we introduce a perturbation of $a(\cdot,\cdot)$, which has a diagonal matrix representation. It is then easy to eliminate the unknowns corresponding to the bubble functions in $\bm{V}_{b}$.

First, consider a natural decomposition of $\bm{u}\in\bm{V}_{h}$:

[TABLE]

where $\bm{u}^{l}\in\bm{V}_{h,1}$ is the linear part and $\bm{u}^{b}\in\bm{V}_{b}$ is the bubble part. The local bilinear form for $T\in\mathcal{T}_{h}$, $\bm{u}\in\bm{V}_{h}$ and $\bm{v}\in\bm{V}_{h}$ is,

[TABLE]

On each element $T$, introduce

[TABLE]

Then we define a perturbed bilinear form $a^{D}(\cdot,\cdot)$ of $a(\cdot,\cdot)$ as follows,

[TABLE]

As shown in C. Rodrigo , the bilinear forms $a(\cdot,\cdot)$ and $a^{D}(\cdot,\cdot)$ are spectrally equivalent.

Lemma 2.5**.**

(C. Rodrigo Lemma 4.3) The following inequalities hold,

[TABLE]

where $\eta$ depends on the shape regularity of the mesh.

Now, we consider the variational problem with diagonal bubble terms: Find $(\bm{u}_{h}^{n},p_{h}^{n},\beta_{h}^{n},\bm{w}_{h}^{n})\in\bm{V}_{h}\times Q_{h}\times B_{h}\times\bm{W}_{h}$, such that,

[TABLE]

To show (11) is well-posed, we define a bilinear form with $a^{D}(\cdot,\cdot)$ replacing $a(\cdot,\cdot)$ in $\mathcal{B}(\cdot,\cdot)$,

[TABLE]

Similarly, we have the following theorem,

Theorem 2.6**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), the following inf-sup condition holds

[TABLE]

and $\mathcal{B}^{D}(\cdot,\cdot)$ is continuous with respect to ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$, i.e.,

[TABLE]

Here, the constants $\gamma_{\ast\ast}$, $\gamma^{\ast\ast}>0$ are independent of mesh size $h$, time step size $\tau$, and the physical parameters.

Proof.

Based on the equivalence between $a^{D}(\cdot,\cdot)$ and $a(\cdot,\cdot)$, the proof is similar to that of Theorem 2.2. ∎

3 Elimination of Bubbles and Darcy’s Velocity

The perturbation of the bilinear form $a(\cdot,\cdot)$ allows for local elimination of the bubble functions. On the other hand, the normal component of Darcy’s velocity is discontinuous across the interior edges, so the mass matrix $M_{\bm{w}}$ corresponding to Darcy’s velocity is block diagonal. Therefore, both the bubble unknowns and Darcy’s velocity can be eliminated by static condensation. In this section, we discuss such an elimination and the well-posedness of the resulting eliminated linear system.

3.1 Elimination

The variational form (11) can be represented in block matrix form,

[TABLE]

where $\bm{U}_{b}$, $\bm{U}_{l}$, $P$, $B$ and $\bm{W}$ are the unknown vectors for the bubble component of the displacement, the piecewise linear component of the displacement, the pressure, the Lagrange multiplier, and Darcy’s velocity, respectively. The blocks in $\mathcal{A}^{D}$ correspond to the following bilinear forms:

[TABLE]

[TABLE]

[TABLE]

After eliminating the unknowns corresponding to the bubbles and Darcy’s velocity from $\mathcal{A}^{D}$, we arrive at a smaller size matrix as follows,

[TABLE]

where $\displaystyle{B_{p}^{E}=\frac{1}{M}M_{p}+\alpha^{2}B_{b}D_{bb}^{-1}B_{b}^{\top}+\tau B_{\bm{w}}M_{\bm{w}}^{-1}B_{\bm{w}}^{\top}}$, $A_{\bm{u}}^{E}=A_{ll}-A_{bl}^{\top}D_{bb}^{-1}A_{bl}$, $B_{\bm{u}}^{E}=B_{l}-B_{b}D_{bb}^{-1}A_{bl}$. Note that, the size of $\mathcal{A}^{E}$ is same as the classical P1-RT0-P0 discretization.

3.2 Well-posedness of the Eliminated System

We have shown the well-posedness of system (4), and the well-posedness of (12) follows directly since it is obtained by static condensation. However, for the purpose of developing preconditioners for the linear system $\mathcal{A}^{E}$, we show the well-posedness of (12) explicitly with proper chosen weighted norms.

In order to show that the system (12) is well-posed, we group the multiplier and the pressure together because the multiplier is considered to be the trace of the pressure on the element boundaries. $\mathcal{A}^{E}$ can be rewritten as in the following two-by-two block form,

[TABLE]

where

[TABLE]

In order to present the well-posedness theorem, we first give some useful lemmas.

Lemma 3.1**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), the matrix $B_{p\beta}^{E}$ is SPD.

Proof.

For any $\bm{x}=(p,\beta)^{\top}\in Q_{h}\times B_{h}$,

[TABLE]

where

[TABLE]

The matrix $\mathcal{M}$ is symmetric positive semi-definite. The eigenvector corresponding to zero eigenvalue is $(0,1,1)^{\top}$, but the vector $([\frac{1}{M}\|p\|^{2}_{M_{p}}+\alpha^{2}\|B_{b}^{\top}p\|^{2}_{D_{bb}^{-1}}]^{1/2}$, $\sqrt{\tau}\|B_{\bm{w}}^{\top}p\|_{M_{\bm{w}}^{-1}},\sqrt{\tau}\|B_{\beta}^{\top}\beta\|_{M_{\bm{w}}^{-1}})^{\top}$ can not be a multiple of $(0,1,1)^{\top}$, because when $\frac{1}{M}\|p\|^{2}_{M_{p}}+\alpha^{2}\|B_{b}^{\top}p\|^{2}_{D_{bb}^{-1}}=0$, we have $p=0$, which implies $\|B_{\bm{w}}^{\top}p\|_{M_{\bm{w}}^{-1}}=0$. Therefore, $B_{p\beta}^{E}$ is SPD. ∎

Corollary 3.2**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), the matrix $A_{p\beta}^{E}$ is SPD, where

[TABLE]

and $\displaystyle{A_{p}^{E}=B_{p}^{E}+\frac{\alpha^{2}}{\zeta^{2}}M_{p}}$.

Proof.

The proof is similar to that of Lemma 3.1. ∎

Now, we introduce the weighted norm

[TABLE]

The following lemma is useful for proving of the well-posedness of the system (13) with respect to the norm (14).

Lemma 3.3**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), then

[TABLE]

[TABLE]

[TABLE]

Proof.

See (4.19) and (4.22) in James for the proofs of (15) and (16). We only show (17).

[TABLE]

So, $\|(B_{\bm{u}}^{E})^{\top}p\|_{(A_{\bm{u}}^{E})^{-1}}\leq\displaystyle{\frac{1}{\zeta}}\|p\|_{M_{p}}$. ∎

Theorem 3.4**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), then the eliminated system (13) satisfies the following inf-sup condition,

[TABLE]

and the continuity condition,

[TABLE]

Thus, (13) is well-posed with respect to the weighted norm (14).

Proof.

(15) gives a weak inf-sup condition,

[TABLE]

Then for a given $p\in Q_{h}$, there exists $\bm{h}\in\bm{V}_{h}$ such that ((B_{\bm{u}}^{E})^{\top}p,\bm{h})$$\geq\displaystyle{\frac{\gamma^{2}_{B}}{\eta^{2}\zeta^{2}}\|p\|^{2}_{M_{p}}-(D_{bb}^{-1}B_{b}^{\top}p,B_{b}^{\top}p)} and $\|\bm{h}\|^{2}_{A_{\bm{u}}^{E}}=\displaystyle{\frac{\gamma^{2}_{B}}{\eta^{2}\zeta^{2}}\|p\|^{2}_{M_{p}}-(D_{bb}^{-1}B_{b}^{\top}p,B_{b}^{\top}p)}$. For $\bm{x}=(\bm{u},p,\beta)^{\top}$, let $\bm{y}=(\bm{u}+\vartheta\bm{h},p,\beta)^{\top}$, then,

[TABLE]

where $\vartheta_{1}=\displaystyle{\min\{\frac{3\gamma^{2}_{B}}{8\eta^{2}},\frac{1}{2}\}}$, and we choose $\vartheta=\displaystyle{\frac{\alpha}{2}}$. On the other hand, by Cauchy-Schwarz inequality, we have,

[TABLE]

where $\vartheta_{2}=\displaystyle{\max\{\frac{1}{2}+\frac{\gamma^{2}_{B}}{2\eta^{2}},2\}}$. Therefore, (18) holds with $\beta_{1}=\vartheta_{1}/\sqrt{\vartheta_{2}}$. Following from (16), (17) and the Cauchy-Schwarz inequality, we have

[TABLE]

Thus, (19) holds with $\beta_{2}=2$, which concludes the proof. ∎

4 Block Preconditioner

In this section, we use the well-posedness to develop block preconditioners for the linear system $\mathcal{A}^{E}$ (13). Following the general framework developed in D. Loghin ; K.A. Mardal , we first consider block diagonal preconditioners (norm-equivalent preconditioners), then we discuss block triangular preconditioners following the framework developed in D. Loghin ; G. Starke ; A. Klawonn ; Y. Ma for (FOV) equivalent preconditioners. We theoretically show that their performance is robust with respect to the discretization and physical parameters.

4.1 Block Diagonal Preconditioners

Based on the framework proposed in D. Loghin ; K.A. Mardal , a natural choice of a norm-equivalent preconditioner is the Riesz operator with respect to the inner product that induces the weighted norm (14). The Riesz operator for (14) takes the following block diagonal form,

[TABLE]

Then, we have the following theorem on the condition number.

Theorem 4.1**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), then

[TABLE]

Proof.

This result follows from Theorem 3.4 and the theory of condition number in D. Loghin ; K.A. Mardal . ∎

The action of inverting the diagonal blocks is expensive and sometimes infeasible in practice. Thus, we use spectrally equivalent SPD approximations to replace the diagonal blocks, i.e.,

[TABLE]

where $S_{\bm{u}}^{E}$ and $S_{p\beta}^{E}$ are spectrally equivalent to the action of the inverse of the diagonal blocks $A_{\bm{u}}^{E}$ and $A_{p\beta}^{E}$, respectively, i.e.,

[TABLE]

[TABLE]

where the constants $c_{1,\bm{u}}^{E}$, $c_{2,\bm{u}}^{E}$, $c_{1,p\beta}^{E}$ and $c_{2,p\beta}^{E}$ are independent of discretization and physical parameters.

Similarly to Theorem 4.1, we have the following result about the condition number.

Theorem 4.2**.**

If the finite element spaces $Q_{h},B_{h},\bm{W}_{h}$ are chosen as (2), and $\bm{V}_{h}$ as (3), then

[TABLE]

where $\displaystyle{c_{1}=\min\{\frac{1}{c_{2,\bm{u}}^{E}},\frac{1}{c_{2,p\beta}^{E}}\}}$ and $\displaystyle{c_{2}=\max\{\frac{1}{c_{1,\bm{u}}^{E}},\frac{1}{c_{1,p\beta}^{E}}\}}$.

4.2 Block Triangular Preconditioner

Now, we consider more general preconditioners, in particular, block upper and lower triangular preconditioners for the linear system $\mathcal{A}^{E}$. Following D. Loghin , based on the weighted norm (14) and Riesz operator (20), the block lower triangular preconditioner takes the following block form,

[TABLE]

And the inexact block lower triangular preconditioner is,

[TABLE]

Next theorem shows that (23) and $\mathcal{A}^{E}$ are FOV-equivalent and, therefore, $\mathcal{B}_{L}^{E}$ provides a preconditioner for general minimal residual (GMRES) method as suggested in D. Loghin ; G. Starke ; A. Klawonn ; Y. Ma

Theorem 4.3**.**

Assuming a shape regular mesh and the discretization described above, there exists constants $\Sigma_{L}$ and $\Upsilon_{L}$, independent of discretization and physical parameters, such that, for any $\bm{x}=(\bm{u},p,\beta)^{\top}\neq 0$,

[TABLE]

Proof.

By direct computation, we have,

[TABLE]

The matrix in the middle is SPD. Thus, there exists a $\varrho_{1}>0$ such that

[TABLE]

where $\Sigma_{L}=\displaystyle{\varrho_{1}\min\{\frac{\gamma^{2}_{B}}{2\eta^{2}},\frac{1}{2}\}}$. Using inequalities (16), (17) and the Cauchy-Schwarz inequality, we have $\|\mathcal{B}_{L}^{E}\mathcal{A}^{E}\bm{x}\|_{(\mathcal{B}_{D}^{E})^{-1}}\leq\Upsilon_{L}\|\bm{x}\|_{(\mathcal{B}_{D}^{E})^{-1}}$ with $\Upsilon_{L}=2$, which concludes the proof. ∎

Now we prove the inexact block lower triangular preconditioner (24) satisfies the requirements to be an FOV-equivalent preconditioner for the $\mathcal{A}^{E}$ system as well, when the diagonal blocks are solved sufficiently accurately.

Theorem 4.4**.**

Assuming a shape regular mesh, the discretization described above, and that (21) and (22) hold with $c_{2,\bm{u}}^{E}>\frac{1}{2}$, then there exists constants $\Sigma^{I}_{L}$ and $\Upsilon^{I}_{L}$, independent of discretization and physical parameters, such that, for any $\bm{x}=(\bm{u},p,\beta)^{\top}\neq 0$,

[TABLE]

Proof.

By direct computation, we have,

[TABLE]

The matrix in the middle is SPD. Thus, there exists a $\varrho_{2}>0$ such that

[TABLE]

where $\Sigma^{I}_{L}=\displaystyle{\frac{\varrho_{2}}{2c_{2,p\beta}^{E}c_{2,\bm{u}}^{E}}\min\{2c_{2,p\beta}^{E},\frac{\gamma^{2}_{B}}{\eta^{2}},2c_{2,\bm{u}}^{E}-1\}}$. The upper bound follows from the inequalities (16), (17), and (21), and Cauchy-Schwarz inequality, with $\Upsilon^{I}_{L}=\displaystyle{\max\{\frac{3}{(c_{1,\bm{u}}^{E})},\frac{1}{(c_{1,p\beta}^{E})}(1+\frac{2}{c_{1,\bm{u}}^{E}})\}}$. ∎

Similarly, we can also consider the following block upper triangular preconditioner for $\mathcal{A}^{E}$,

[TABLE]

and its corresponding inexact version,

[TABLE]

The following theorems show that these block upper triangular preconditioners are also parameter robust. The proofs are similar to the proofs for Theorem 4.3 and 4.4 and, therefore, are omitted.

Theorem 4.5**.**

Assuming a shape regular mesh and the discretization described above, then there exists constants $\Sigma_{U}$ and $\Upsilon_{U}$, independent of discretization and physical parameters, such that, for any $\bm{x}^{\prime}=(\mathcal{B}_{U}^{E})^{-1}\bm{x}$, with $\bm{x}=(\bm{u},p,\beta)^{\top}\neq 0$,

[TABLE]

Theorem 4.6**.**

Assuming a shape regular mesh, the discretization described above, and that (21) and (22) hold with $c_{2,\bm{u}}^{E}>\frac{1}{2}$, then there exists constants $\Sigma^{I}_{U}$ and $\Upsilon^{I}_{U}$, independent of discretization or physical parameters, such that, for any $\bm{x}^{\prime}=(\widehat{\mathcal{B}_{U}^{E}})^{-1}\bm{x}$, with $\bm{x}=(\bm{u},p,\beta)^{\top}\neq 0$,

[TABLE]

5 Numerical Experiments

In this section, we give some numerical examples. The examples in Section 6.1 are used to investigate the accuracy of the stabilized hybrid method. In addition, the cantilever bracket problem is proposed to show the effectiveness of our discretization. In Section 6.2, we demonstrate the robustness of the preconditioners presented in Section 5. In all test cases we consider a diagonal permeability tensor $\kappa=K\bm{I}$ with constant $K$.

5.1 The Accuracy and Efficiency of the Stabilization Method

Firstly, we demonstrate numerically that for Biot’s model, the error in the finite element approximation does not decrease when the permeability is small relative to the mesh size. We consider $\Omega=(0,1)\times(0,1)$ with homogeneous Dirichlet boundary conditions for $\bm{u}$, and Neumann boundary conditions for $p$. We cover $\Omega$ with a uniform triangular grid by dividing an $N\times N$ uniform square mesh into right triangles, where the mesh spacing is defined by $h=1/N$. The material parameters are $\lambda=2$, $\mu=1$, $\alpha=1$, and $M=10^{6}$. We consider a diagonal permeability tensor $\kappa=K\bm{I}$ with constant $K$. The data is set so that the exact solution is given by

[TABLE]

Note that the solution is designed to satisfy $\nabla\cdot\bm{u}=0$ at any time $t$. We set $\tau=1$ and $t_{max}=1$.

As shown in Table 1, the results for the energy norm ($\|\bm{v}\|^{2}_{a}:=a(\bm{v},\bm{v})$) of the displacement errors and the $L^{2}$-norm for pressure errors convergence when $K$ is relatively big. When $K$ becomes smaller, the performance becomes worse, and eventually divergences when $K$ is really small.

Then, we use the proposed stabilized hybrid discretization to solve it. Table 2 shows good convergence results even when $K\rightarrow 0$.

We also compare the errors obtained by (11) with diagonal bubble functions with those provided by the bubble enriched system (4), in order to see that the same error reduction is achieved. Fig. 1 displays a comparison of the displacement and pressure errors in the energy and $L^{2}$ norms, respectively, for different grid sizes. We choose $K=10^{-8}$ here, though similar results are obtained for different values of $K$. We observe both schemes produces the same convergence rate although the scheme corresponding to the diagonal version provides slightly worse errors. However, this scheme, when the bubble block is eliminated, uses fewer degrees of freedom, thus computationally more efficient.

We also consider the cantilever bracket problem P. J. Phillips4 illustrated in Fig. 2. The domain is the unit square $[0,1]\times[0,1]$. For the flow problem we impose a no-flow boundary condition along the entire boundary. For the elasticity problem, we assume that the left side boundary of the domain is fixed, so a no-displacement boundary condition is imposed. We impose a downward traction at the top of the domain and a traction-free boundary condition at the right and bottom of the domain. The initial displacement and pressure are assumed to be zero.

For the material properties, the Lamé coefficients are computed in terms of the Young modulus, $E$, and the Poisson ratio $\nu$: $\lambda=\frac{E\nu}{(1-2\nu)(1+\nu)}$, and $\mu=\frac{E}{1+2\nu}$ with $\nu=0.45$, $E=10^{5}$, and $K=10^{-7},\ \alpha=0.93,\ M=10^{10}$. Setting $\tau=0.001$, we show the approximation for the pressure field obtained without stabilization terms in Fig. 3 (left), and observe nonphysical oscillations. We clearly observe that the stability scheme (13) removes the nonphysical pressure oscillations, see Fig. 3 (right).

5.2 Robustness of the Preconditioners

In this section, we demonstrate the robustness of the block preconditioners. We use flexible GMRES to solve the eliminated system $\mathcal{A}^{E}$ (13). For each test we use flexible GMRES to solve the linear system $\mathcal{A}^{E}$ and a stopping tolerance of $10^{-8}$ for the relative residual is used. Different numerical tests are presented to show the robustness of the proposed preconditioners with respect to the physical parameters (mainly the permeability, $\kappa$ and the poisson ratio $\nu$) and the discretization parameters (mesh size $h$ and time step size $\tau$). In the exact block preconditioners $\mathcal{B}_{D}^{E}$, $\mathcal{B}_{L}^{E}$, and $\mathcal{B}_{U}^{E}$, the diagonal blocks are inverted by the UMFPACK library T. A. Davis1 ; T. A. Davis2 ; T. A. Davis3 ; T. A. Davis04 . For the inexact block preconditioners $\widehat{\mathcal{B}_{D}^{E}}$, $\widehat{\mathcal{B}_{L}^{E}}$, and $\widehat{\mathcal{B}_{U}^{E}}$, $S_{\bm{u}}^{E}$ and $S_{p,\beta}^{E}$ are implemented using conjugate gradient (CG) method preconditioned with a V-cycle unsmoothed aggregation AMG (UA-AMG) and solved to a tolerance that the relative residual is less than $10^{-3}$. We use zero right hand side and random initial guess in all the tests. Each test is repeated $5$ times and the average number of iterations are reported.

We consider the first example presented in Section 5.1. Choosing the physical parameters as $M=10^{6}$, $\alpha=1$, $E=1$, Table 3 shows iteration counts for the block preconditioners for the eliminated system when the physical values of $\nu$ and $K$ are varying. The mesh size is fixed as $h=\frac{1}{64}$, and the time step size is $\tau=1$. When $K$ is varying, the number of iterations grows slightly at the beginning but stabilizes when $K$ approached $10^{-12}$ for block diagonal preconditioners $\mathcal{B}_{D}^{E}$ and $\widehat{\mathcal{B}_{D}^{E}}$. Block triangular preconditioners perform better than the block diagonal preconditioners as expected and are stable with respect to $K$. Overall, from the relatively consistent iteration counts, we can conclude that the block preconditioners are indeed robust with respect to the physical parameters. Table 4 shows the number of iterations for the block preconditioners for solving $\mathcal{A}^{E}$ for different mesh sizes and time step sizes with fixed $\lambda=2$, $\mu=1$, and $K=10^{-6}$. Again, we observe robustness with respect to the discretization parameters. The block upper and lower triangular preconditioners contain more coupling information than the block diagonal preconditioners, and, as a result, we can see that they perform better than the block diagonal preconditioners. In addition, from both tables, we can see that the number of iterations of the inexact block preconditioners is slightly higher than that of the exact block preconditioners. This is expected and the difference is by no means significant. Therefore, we would suggest using the inexact version in practice to avoid expensive direct solves used in the exact version.

Next, we illustrate that the block preconditioners are robust for the cantilever bracket problem as well. Table 5 shows iteration counts for the block preconditioners when the values of $\nu$ and $K$ varies for the eliminated system (13). The mesh size is fixed as $h=1/64$ and the time step size is $\tau=1$. When $K$ is varying, the number of iterations of the diagonal block preconditioners grows as $K$ approaches $10^{-12}$ but stabilize towards the end. The exact block triangular preconditioners are stable with respect to $K$ and the number of iterations basically stays the same. But the iteration counts of the inexact block triangular preconditioners increases as $K$ approached to $10^{-12}$. This is due to the inexactness and, when we solve the diagonal block more accurately, the number iterations will stabilize. When $\nu$ is varying, the number iterations decreases as $\nu$ gets closer to $0.5$ for all block preconditioners. However, we want to point out that, when $\nu$ approaches $0.5$, the diagonal block corresponding to the elasticity part becomes more and more difficult to solve and, therefore, the inner preconditioned CG methods might need more iterations since we are using simple UA-AMG method. Overall, the proposed block preconditioners are robust with respect to the physical parameters which confirms the theoretical results.

Table 6 shows the iteration counts for the block preconditioners on the eliminated system for different mesh sizes and time step sizes with $\nu=0.45$ and $K=10^{-7}$. As shown in Table 6, We saw the same behavior as the previous experiments. The block triangular preconditioners work better than the block diagonal preconditioners. The number of iterations of the inexact block preconditioners is higher than the exact block preconditioners. Overall, the number of iterations of all the block preconditioners are stable when we varying $h$ and $\tau$, which demonstrates the robustness of the block preconditioners. In this example, the exact block triangular preconditioners only need $2$-$4$ iterations to achieve tolerance. This is due to the choice of the physical parameters, the block triangular preconditioned linear system has eigenvalues are clustered near $-1$ and $1$ and, therefore, the Krylov iterative method converges quickly in this case.

6 Conclusion

In this work, we considered a stabilized hybrid mixed method for the three-field Biot’s model. The hybrid mixed finite element method based on the triple P1-RT0-P0 presented in C. Niu is not uniformly stable with respect to the physical parameters, for instance, when the permeability is small with respect to the mesh size, it does not converge. To overcome such problem, we presented a stabilization technique with bubble functions following V. Girault ; C. Rodrigo . The well-posedness of the stabilized scheme with respect to the energy norm is given independent of the physical and discretization parameters. We then use a spectrally equivalent perturbed bilinear form as suggested in C. Rodrigo . The perturbation of the bilinear form allows for elimination of the bubble functions; The normal component of Darcy’s velocity is discontinuous across the interior edges, so the mass matrix corresponding to Darcy’s velocity is block diagonal. Therefore, the unknowns of the bubble functions and Darcy’s velocity can be eliminated by static condensation. In fact, the eliminated system is the same size as the P1-RT0-P0 discretization, which is widely used in practice. We also prove that the eliminated system is well-posed and efficient preconditioners are developed to solve the resulting linear system. We show theoretically that the block preconditioners are robust with respect to physical and discretization parameters. Finally, numerical experiments validate the accuracy and efficiency of the stabilization method and also demonstrate the robustness of the block preconditioners.

Appendix. Proof of Theorem 2.4

In this appendix, we provide the proof of Theorem 2.4.

Proof.

First of all, we introduce the notation $\beta\in L^{2}(\partial\Omega)$ for pressure on the boundary of domain $\Omega$. Here, $L^{2}(\partial\Omega)$ denotes the set of square integrable functions on the boundary of $\Omega$. Then, for the model problem (1), integrating by parts, we have,

[TABLE]

Choosing $\bm{v}=\bm{v}_{h}$, $q=q_{h}$, $\rho=\rho_{h}$ and $\bm{r}=\bm{r}_{h}$ in (27) and subtracting these equations from (4) and using the definition of elliptic projections given in (8), we obtain

[TABLE]

where $(R^{n},q_{h})=\displaystyle{\alpha(\nabla\cdot R_{\bm{u}}^{n},q_{h})+\frac{1}{M}(R_{p}^{n},q_{h})}$, with

[TABLE]

Then, choosing $\bm{v}_{h}=\overline{\partial}_{t}e_{\bm{u}}^{n}$, $q_{h}=e_{p}^{n}$, $\rho_{h}=e_{\beta}^{n}$, and $\bm{r}_{h}=e_{\bm{w}}^{n}$, respectively, and adding all the equations of (28), we have,

[TABLE]

Due to our choice of finite element spaces and the first equality in (28), we have

[TABLE]

Therefore, we have

[TABLE]

This implies,

[TABLE]

and

[TABLE]

By recursion, we get,

[TABLE]

and

[TABLE]

Combining (29) and (30), we have the estimate

[TABLE]

Now, following the same procedures of Lemma 8 in C. Rodrigo2 , we have

[TABLE]

and similarly,

[TABLE]

Then, the error estimate (10) follows from (31), (32), (33), (9) and triangle inequality. ∎