# Convergence Analysis of H(div)-Conforming Finite Element Methods for a Nonlinear Poroelasticity Problem

Authors: Yuping Zeng; Zhifeng Weng; Fen Liang
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In this paper, we introduce and analyze H(div)-conforming finite element methods for a nonlinear model in poroelasticity. More precisely, the flow variables are discretized by H(div)-conforming mixed finite elements, while the elastic displacement is approximated by the H(div)-conforming finite element with the interior penalty discontinuous Galerkin formulation. Optimal a priori error estimates are derived for both semidiscrete and fully discrete schemes.


---

## Body

## 1. Introduction

The Biot poroelasticity model [1] describes the phenomena of coupled mechanics and flow in porous media, and it is of increasing importance today in a divergence of science and engineering applications. Applications of poroelasticity have been made to areas that include carbon sequestration in environment engineering, predictive ability in earthquake engineering, surface subsidence in field phenomena, logging technologies in reservoir engineering, and pathological condition in biomechanics.In the context of linear poroelasticity systems, a large number of numerical methods have been developed in recent years. Among them, finite element methods (FEMs) are commonly used approaches for such coupled system. In [2–4], Taylor–Hood elements were applied for the approximation of the displacement field and fluid pressure. Recently, some techniques based on mixed formulations have been developed. A method coupling standard conforming FEM for the displacement with mixed FEMs for flow variables has been provided in [5, 6]. However, it is well known that a conforming FEM approximation of displacement may give rise to locking or nonphysical oscillations [7]. Some remedies of locking include nonconforming FEMs [8–10], mixed FEMs [11, 12], discontinuous Galerkin (DG) methods [13–15], hybrid high-order methods [16], and weak Galerkin methods [17–20]. We also refer the interested reader to [21–24] for further study of locking-free FEMs on linear poroelasticity. Moreover, some preconditioners that are robust with respect to the physical and discretization parameters have been developed in this direction (refer to [25–28] and the references therein).Compared with linear model in poroelasticity, the numerical methods on nonlinear problem are still very rare. The goal of this paper is to develop H(div)-conforming finite element methods for a nonlinear poroelasticity model. The nonlinearity arises from the dependence of the hydraulic conductivityκ on the dilation ∇⋅u, that is, κ=κ∇⋅u. Such model has been studied in [29], where the existence and uniqueness of the weak solution were proved. Moreover, the authors in [29] introduce a FEM approximation and give the optimal order error estimate. We also find that similar problem has been treated in [30], where a mixed method combining conforming and mixed FEMs was proposed and analyzed. Some other works of numerical methods for nonlinear poroelasticity model can be found in [31, 32]. In the current work, we will propose and analyze H(div)-conforming methods to solve the nonlinear problem studied in [29]. More precisely, the flow variables are discretized by H(div)-conforming mixed FEMs, while the displacement is approximated by the H(div)-conforming FEM with interior penalty DG method. Combining H(div)-conforming finite elements with DG methods was initially proposed in [33, 34] (see also [35–37]); they mainly intended to solve Stokes and Navier–Stokes equations in fluid mechanics. Later, this method was extended to more complex Darcy–Stokes interface problems [38–40], Brinkman problems [41], a magnetic induction model [42], and linear Biot model in poroelasticity [43]. We also refer the interested reader to [44] on robust numerical simulations of H(div)-conforming FEMs for Stokes, coupled Darcy–Stokes, and Brinkman problems. Our work is to extend H(div)-conforming FEMs of linear Biot model to the nonlinear case. With the help of the framework presented in [30], we will give a detailed a priori error analysis of the proposed numerical scheme. It is worth mentioning that we shall address some difficulties arising from the inherent nonlinearity (see the terms (43), (44), (68), and (69) below). The main advantage of our method is that we present a unified treatment for both flow variables and the displacement. In addition, we relax H1 conformity by using H(div) space to approximate the displacement, and thus the normal components are continuous and the tangential components are treated by interior penalty DG approach; this implies that our method is locally conservative.The rest of this paper is organized as follows. In Section2, we describe the nonlinear poroelasticity model and present its mixed weak formulation. We propose a spatial semidiscrete scheme which is based on H(div)-conforming finite elements in Section 3, where a priori error estimate is also provided. In Section 4, a fully discrete scheme with the backward Euler time stepping is proposed and analyzed. Some conclusions are made in Section 5.


## 2. Preliminaries and Model Problem

### 2.1. The Model Problem

LetΩ⊂ℝ2 be a bounded convex polygonal domain, with Lipschitz boundary ∂Ω. We consider the following nonlinear poroelastic model in Ω over a time interval 0,T (see [29]):(1)c0∂tp+α∇⋅∂tu−∇⋅κ∇⋅u∇p=g,inΩ×0,T,−∇⋅σ=f,inΩ×0,T,where(2)σ=λtrϵuI+2μϵu−αpI,ϵu=12∇u+∇uT,in which ux,t is the displacement of the solid phase, px,t is the fluid pressure, σx,t is the total stress tensor, c0>0 is the storage coefficient, 0<α≤1 is the Biot–Willis constant, and λ and μ are Lamé constants; κ⋅ is the hydraulic conductivity that satisfies assumptions (6) and (7) below.Denote byΓd (resp. Γt) the Dirichlet (resp. traction) boundary for the elastic variables. Similarly, denote by Γp (resp. Γf) the pressure Dirichlet (resp. fluid normal flux) boundary. We assume that ∂Ω=Γd∪Γt with Γd>0, and ∂Ω=Γp∪Γf. The boundary conditions and initial conditions for the above system read as(3)u=0,onΓd×0,T,σn=0,onΓt×0,T,p=0,onΓp×0,T,κ∇⋅u∇p⋅n=0,onΓf×0,T,p⋅,0=p0,inΩ,u⋅,0=u0,inΩ,where n denotes the outward unit normal vector.Settingq=−κ∇⋅u∇p, we can reformulate equation (1) by(4)c0∂tp+α∇⋅∂tu+∇⋅q=g,inΩ×0,T,κ−1∇⋅uq+∇p=0,inΩ×0,T,−∇⋅2μεu−λ∇∇⋅u+α∇p=f,inΩ×0,T.


### 2.2. Weak Formulation

We start by introducing some notation. LetHsDs≥0 be the standard Sobolev space, equipped with norm ⋅s,D and seminorm ⋅s,D. If s=0, H0D is understood as L2D. Denote w,vD=∫Dwvdx and w,v∂D=∫∂DwvdF. When D=Ω, we will omit the index Ω. For the space HsD=HsD2, its norm is still denoted by ⋅s,D. A subspace of H1Ω with vanishing trace on Γd is given by H0,Γd1Ω=v∈H1Ω:vΓd=0. Additionally, we define Hdiv;Ω=v∈L2Ω2:∇⋅v∈L2Ω with its graph norm vdiv=v02+∇⋅v021/2. Two subspaces are then introduced by H0,Γfdiv;Ω=v∈Hdiv;Ω:v⋅nΓf=0 and H0,Γddiv;Ω=v∈Hdiv;Ω:v⋅nΓd=0. For convenience, we set P=L2Ω, Q=H0,Γfdiv;Ω and V=H0,Γd1Ω2.Integrating by parts, we deduce that the mixed weak formulation for (4) reads as follows: find p,q,u∈P×Q×V for every t∈0,T such that(5)c0∂tp,w+α∇⋅∂tu,w+∇⋅q,w=g,w,κ−1∇⋅uq,z−p,∇⋅z=0,au,v−αp,∇⋅v=f,v,for any w,z,v∈P×Q×V, where au,v=2μϵu,ϵv+λ∇⋅u,∇⋅v. Additionally, we assume that there exist positive ka and kb such that(6)ka≤κ−1ϑ≤kb,∀ϑ∈ℝ,and κ−1 is Lipschitz continuous:(7)κ−1ϑ−κ−1η≤Lϑ−η,forϑ,η∈ℝ.To deal with functions of time and space, we introduce the standard Bochner spaceLp0,T;HsΩ, which consists of all functions u:0,T⟶HsΩ with norm(8)uLp0,T;HsΩ=∫0Tutspdt1/p,for 1≤p<∞. When p=∞, the norm is defined as(9)uL∞0,T;HsΩ=sup0≤t≤Tuts.


## 2.1. The Model Problem

LetΩ⊂ℝ2 be a bounded convex polygonal domain, with Lipschitz boundary ∂Ω. We consider the following nonlinear poroelastic model in Ω over a time interval 0,T (see [29]):(1)c0∂tp+α∇⋅∂tu−∇⋅κ∇⋅u∇p=g,inΩ×0,T,−∇⋅σ=f,inΩ×0,T,where(2)σ=λtrϵuI+2μϵu−αpI,ϵu=12∇u+∇uT,in which ux,t is the displacement of the solid phase, px,t is the fluid pressure, σx,t is the total stress tensor, c0>0 is the storage coefficient, 0<α≤1 is the Biot–Willis constant, and λ and μ are Lamé constants; κ⋅ is the hydraulic conductivity that satisfies assumptions (6) and (7) below.Denote byΓd (resp. Γt) the Dirichlet (resp. traction) boundary for the elastic variables. Similarly, denote by Γp (resp. Γf) the pressure Dirichlet (resp. fluid normal flux) boundary. We assume that ∂Ω=Γd∪Γt with Γd>0, and ∂Ω=Γp∪Γf. The boundary conditions and initial conditions for the above system read as(3)u=0,onΓd×0,T,σn=0,onΓt×0,T,p=0,onΓp×0,T,κ∇⋅u∇p⋅n=0,onΓf×0,T,p⋅,0=p0,inΩ,u⋅,0=u0,inΩ,where n denotes the outward unit normal vector.Settingq=−κ∇⋅u∇p, we can reformulate equation (1) by(4)c0∂tp+α∇⋅∂tu+∇⋅q=g,inΩ×0,T,κ−1∇⋅uq+∇p=0,inΩ×0,T,−∇⋅2μεu−λ∇∇⋅u+α∇p=f,inΩ×0,T.


## 2.2. Weak Formulation

We start by introducing some notation. LetHsDs≥0 be the standard Sobolev space, equipped with norm ⋅s,D and seminorm ⋅s,D. If s=0, H0D is understood as L2D. Denote w,vD=∫Dwvdx and w,v∂D=∫∂DwvdF. When D=Ω, we will omit the index Ω. For the space HsD=HsD2, its norm is still denoted by ⋅s,D. A subspace of H1Ω with vanishing trace on Γd is given by H0,Γd1Ω=v∈H1Ω:vΓd=0. Additionally, we define Hdiv;Ω=v∈L2Ω2:∇⋅v∈L2Ω with its graph norm vdiv=v02+∇⋅v021/2. Two subspaces are then introduced by H0,Γfdiv;Ω=v∈Hdiv;Ω:v⋅nΓf=0 and H0,Γddiv;Ω=v∈Hdiv;Ω:v⋅nΓd=0. For convenience, we set P=L2Ω, Q=H0,Γfdiv;Ω and V=H0,Γd1Ω2.Integrating by parts, we deduce that the mixed weak formulation for (4) reads as follows: find p,q,u∈P×Q×V for every t∈0,T such that(5)c0∂tp,w+α∇⋅∂tu,w+∇⋅q,w=g,w,κ−1∇⋅uq,z−p,∇⋅z=0,au,v−αp,∇⋅v=f,v,for any w,z,v∈P×Q×V, where au,v=2μϵu,ϵv+λ∇⋅u,∇⋅v. Additionally, we assume that there exist positive ka and kb such that(6)ka≤κ−1ϑ≤kb,∀ϑ∈ℝ,and κ−1 is Lipschitz continuous:(7)κ−1ϑ−κ−1η≤Lϑ−η,forϑ,η∈ℝ.To deal with functions of time and space, we introduce the standard Bochner spaceLp0,T;HsΩ, which consists of all functions u:0,T⟶HsΩ with norm(8)uLp0,T;HsΩ=∫0Tutspdt1/p,for 1≤p<∞. When p=∞, the norm is defined as(9)uL∞0,T;HsΩ=sup0≤t≤Tuts.


## 3. Error Estimates for the Semidiscrete Scheme

In this section, we aim to present the semidiscrete numerical method which focuses on discretizing the spatial variables. LetTh be a shape-regular decomposition of Ω into triangles K. hK denotes the diameter of K and h=maxK∈ThhK. Denote by ℰhI the set of interior edges of elements in Th, ℰhd the set of boundary edges on Γd, and ℰht the set of boundary edges on Γt. Therefore, the set of all edges ℰh=ℰhI∪ℰhd∪ℰht. The length of the edge e∈ℰh is denoted by he. For every edge e∈ℰh, we fix a unit normal ne such that for edges on the boundary ∂Ω, ne is the outward unit normal.Fore∈ℰhI shared by elements K+ and K−, let v be a scalar or vector piecewise smooth function and set v±=ve∩∂K±, and we define the average and jump by(10)v=12v++v−andv=v+−v−.On a boundary edge, we set(11)v=vandv=v.Consider the following finite element spaces:(12)Qh=q∈H0,Γfdiv;Ω:qK∈BDMkK,(13)Vh=v∈H0,Γddiv;Ω:vK∈BDMkK,(14)Ph=w∈L2Ω:wK∈Pk−1K,where BDMkk≥1 is the H(div)-conforming space introduced by Brezzi, Douglas, and Marini and PkK denotes the space of polynomials of degree less than or equal to k on K.LetΠh:Qresp.V⟶Qhresp.Vh be the BDM interpolation and Ph:P⟶Ph be the L2 projection. It is well known that they satisfy the following properties [45, 46]:(15)z−Phz,w=0,∀w∈Ph,(16)z−Phz0,K≤Chlzl,K,0≤l≤k,(17)∇⋅v−Πhv,w=0,∀w∈Ph,(18)v−Πhvs,K≤Chl−svl,K,s=0,1,1≤l≤k+1,(19)∇⋅v−Πhv0,K≤Chl∇⋅vl,K,0≤l≤k.Here and throughout the paper, we utilizeC to represent a positive generic constant that is independent of h and Δt but may take different values at different occurrences.


### 3.1. Semidiscrete H(div)-Conforming Methods with DG Formulation

The expression ofah⋅,⋅ below is similar to the statement in Section 3.1 in the work [43]. For completeness, we derive its scheme as follows.Multiplying the third equation in (4) by any v∈Vh on each element K, integrating by parts, and then summing over all elements in Th, we arrive at(20)−∑T∈Th∇⋅σ,vT=∑T∈Thσ,∇vT−∑T∈ThσnT,v∂T=f,v,where nT is outward normal of T.We first have(21)∑T∈ThσnT,v∂T=∑e∈ℰhIσne,ve+∑e∈ℰhσne,ve=∑e∈ℰhσne,ve=∑e∈ℰhI∪ℰhdσne,ve=2μ∑e∈ℰhI∪ℰhdεune,ve+λ∑e∈ℰhI∪ℰhd∇⋅une,ve+α∑e∈ℰhI∪ℰhdpne,ve=2μ∑e∈ℰhI∪ℰhdεune,ve.Here in the third line, we have usedσne=0 in ℰht, and in the last line, we have used the fact that v⋅ne=0 on e∈ℰhI∪ℰhd and the regularity of the exact solutions.In addition, it is easy to check that(22)∑T∈Thσ,∇vT=2μ∑T∈Thϵu,ϵvT+λ∑T∈Th∇⋅u,∇⋅vT−α∑T∈Th∇⋅v,pT=2μ∑T∈Thϵu,ϵvT+λ∇⋅u,∇⋅v−α∇⋅v,p.Substituting (21) and (22) into (20) yields(23)2μ∑K∈Thϵu,ϵvT−2μ∑e∈εhI∪εhdϵune,ve+λ∇⋅u,∇⋅v−α∇⋅v,p=f,v,∀v∈Vh.Letτe be the unit tangential direction to the edge e so that ne and τe form a right-hand side coordinate system; it follows that for any w, we have(24)w=w⋅nene+w⋅τeτe.Therefore,(25)ϵune⋅v=ϵune⋅nene+ϵune⋅τeτe⋅v⋅nene+v⋅τeτe=ϵune⋅nev⋅ne+ϵune⋅τev⋅τe.Combining (18) and (19) and the regularity of the exact solutions implies that(26)2μ∑e∈ℰhI∪ℰhdϵune,ve=2μ∑e∈ℰhI∪ℰhd∫eϵune⋅vdF=2μ∑e∈ℰhI∪ℰhd∫eϵune⋅nev⋅nedF+2μ∑e∈ℰhI∪ℰhd∫eϵune⋅τev⋅τedF=2μ∑e∈ℰhI∪ℰhd∫eϵune⋅τev⋅τedF=2μ∑e∈ℰhI∪ℰhdϵune⋅τe,v⋅τee.Here in the third line, we have used the fact thatv⋅ne=0 on e∈ℰhI∪ℰhd. Therefore, (23) is reduced to(27)2μ∑K∈Thϵu,ϵvT−2μ∑e∈ℰhI∪ℰhdϵune⋅τe,v⋅τee+λ∇⋅u,∇⋅v−α∇⋅v,p=f,v,∀v∈Vh.Adding some stabilized terms as in [47], we propose the following DG method for the third equation in (4):(28)ahu,v−αp,∇⋅v=f,v,with(29)ahu,v=2μ∑K∈Thϵu,ϵvT+λ∇⋅u,∇⋅v−2μ∑e∈ℰhI∪ℰhdϵune⋅τe,v⋅τee−2μ∑e∈ℰhI∪ℰhdϵvne⋅τe,u⋅τee+∑e∈ℰhI∪ℰhd2μβheu⋅τe,v⋅τee.Consequently, after integrating by parts, we find that the exact solutions of (4) satisfy(30)c0∂tp,w+α∇⋅∂tu,w+∇⋅q,w=g,w,κ−1∇⋅uq,z−p,∇⋅z=0,ahu,v−αp,∇⋅v=f,v,for any w,z,v∈Ph×Qh×Vh.Thus, the corresponding semidiscrete H(div)-conforming FEMs for (4) can be designed as follows: find ph,qh,uh∈Ph×Qh×Vh such that(31)c0∂tph,w+α∇⋅∂tuh,w+∇⋅qh,w=g,w,κ−1∇⋅uhqh,z−ph,∇⋅z=0,ahuh,v−αph,∇⋅v=f,v,for any w,z,v∈Ph×Qh×Vh, with the initial conditions given by(32)ph0=Php0,uh0=Πhu0.


### 3.2. Error Estimates

We begin by defining the mesh-dependent norm⋅h on Vh=V+Vh:(33)vh=∑K∈Th2μϵv0,K2+∑e∈ℰhI∪ℰhd2μhe−1v⋅τe0,e2+∑e∈ℰhI∪ℰhd2μheϵvne⋅τe0,e2+λ∇⋅v021/2.Then, from [43], we can obtain the following result which shows the continuity and coercivity of the bilinear ah⋅,⋅.Lemma 1.
There exists a constantCa>0 such that(34)ahx,v≤Caxhvh,∀x,v∈Vh.Moreover, for sufficiently large penalty parameterβ, it holds that(35)ahv,v≥Cbvh2,∀v∈Vh.We are now in a position to state the following error estimate, which is the main result of this section.Theorem 1.
Letp,q,u∈P×Q×V and ph,qh,uh∈Ph×Qh×Vh be the solutions of (4) and (31), respectively. Then, the following holds:(36)u−uhL∞0,T;Eh+p−phL∞0,T;L2Ω+q−qhL20,T;L2Ω≤Chk,where xL∞0,T;Eh=sup0≤s≤Txsh and C depends on the regularity of u,ut,p,pt, and q.Proof.
Subtracting (30) from (31) yields(37)c0∂tp−ph,w+α∇⋅∂tu−uh,w+∇⋅q−qh,w=0,κ−1∇⋅uq−κ−1∇⋅uhqh,z−p−ph,∇⋅z=0,ahu−uh,v−αp−ph,∇⋅v=0,for any w,z,v∈Ph×Qh×Vh.
We then split the errorp−ph into p−ph=ξp+θp with ξp=p−Php and θp=Php−ph. Similarly, q−qh=ξq+θq with ξq=q−Πhq and θq=Πhq−qh; u−uh=ξu+θu with ξu=u−Πhu and θu=Πhu−uh. Since the estimates for ξp, ξq, and ξu can be derived by the interpolation error bounds in (16), (18), and (19), it remains to estimate θp, θq, and θu. To this end, using (15) and (17), we can rewrite (37) as(38)c0∂tθp,w+α∇⋅∂tθu,w+∇⋅θq,w=0,κ−1∇⋅uq−κ−1∇⋅uhqh,z−θp,∇⋅z=0,ahθu,v−αθp,∇⋅v=−ahξu,v,for any w,z,v∈Ph×Qh×Vh.
Settingw=θp, z=θq, and v=∂tθu in the above equations, we infer that(39)c02ddtθp,θp+α∇⋅∂tθu,θp+∇⋅θq,θp=0,κ−1∇⋅uq−κ−1∇⋅uhqh,θq−θp,∇⋅θq=0,12ddtahθu,θu−αθp,∇⋅∂tθu=−ahξu,∂tθu.
Summing equation (39), integrating in time from 0 to t≤T, and noting that θp0=0, θu0=0, and(40)κ−1∇⋅uq−κ−1∇⋅uhqh=κ−1∇⋅uq−Πhq+κ−1∇⋅u−κ−1∇⋅uhΠhq+κ−1∇⋅uhΠhq−qh,we obtain(41)12ahθut,θut+c02θpt02+∫0tκ−1∇⋅uhθqs,θqsds=A1+A2+A3,with(42)A1=−∫0tκ−1∇⋅uξqs,θqsds,A2=−∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs,θqsds,A3=−∫0tahξus,∂tθusds.
For the first termA1, using (6) and Cauchy–Schwarz and Young’s inequalities, we find that(43)A1=−∫0tκ−1∇⋅uξqs,θqsds≤∫0tκ−1∇⋅uξqs0θqs0ds≤14ϵ1∫0t∥κ−1∇⋅uξqs∥02ds+ɛ1∫0tθqs02ds≤kb24ε1∫0tξqs02ds+ɛ1∫0tθqs02ds,where ɛ1 is an arbitrarily small number.
For the second termA2, with the aid of (6) and (7) and Cauchy–Schwarz and Young’s inequalities, we deduce that(44)A2=−∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs,θqsds≤∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs0θqs0ds≤14ɛ2∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs02ds+ɛ2∫0tθqs02ds≤14ɛ2∫0tκ−1∇⋅u−κ−1∇⋅uh02Πhqs∞2ds+ɛ2∫0tθqs02ds≤L4ɛ2∫0t∇⋅u−∇⋅uh02Πhqs∞2ds+ɛ2∫0tθqs02ds≤LC14ɛ2∫0tus−uhsh2Πhqs∞2ds+ɛ2∫0tθqs02ds≤LC12ɛ2∫0tus−uhsh2Πhqs∞2ds+LC12ɛ2∫0tΠhus−uhsh2Πhqs∞2ds+ɛ2∫0tθqs02ds≤LC1Cq22ɛ2∫0tus−Πhush2ds+LC1Cq22ɛ2∫0tΠhus−uhsh2ds+ɛ2∫0tθqs02ds=LC1Cq22ɛ2∫0tξush2ds+LC1Cq22ɛ2∫0tθus02dsɛ2∫0tθqs02ds,where ϵ2 is an arbitrarily small number and Cq=max0≤s≤tΠhqs∞.
To estimate the third termA3, integrating by parts and noting that θu0=0, we first obtain(45)A3=∫0tah∂tξus,θusds−ahξut,θut,which combined with (34) and Young’s inequality implies that(46)A3≤Ca2∫0t∂tξush2+θush2ds+Ca24ɛ3ξuth2+ɛ3θuth2,where ɛ3 is an arbitrarily small number.
Combining the bounds in (41)–(46) and using (6) and (35), we have(47)Cb2−ɛ3θuth2+c02θpt02+ka−ɛ1−ɛ2∫0tθqs02ds≤LC1Cq22ɛ2∫0tθush2ds+LC1Cq22ɛ2∫0tξush2ds+kb24ɛ1∫0tξqs02ds+Ca2∫0t∂tξush2ds+Ca24ɛ3ξuth2.
Now, we choose appropriateɛii=1,2,3 such that Cmin=minCb/2−ɛ3,c0/2,ka−ɛ1−ɛ2>0. Consequently, the above inequality still holds if we replace its left-hand side term by Cminθuth2+θpt02+∫0tθqs02ds. Then, dividing both sides of the above inequality by Cmin and using Gronwall’s inequality, we have(48)θuth2+θpt02+∫0tθqs02ds≤C∫0tξqs02+∂tξush2+ξush2ds+ξuth2.
Noting that the above estimate holds for all0≤t≤T and using appropriate approximation properties stated in (16), (18), and (19), we obtain(49)sup0≤s≤Tθush2+sup0≤s≤Tθps02+∫0T∫0Tθqs02ds≤Ch2k∫0Tqsk2+∂tusk+12+usk+12ds+h2ksup0≤s≤Tusk+12.
This can be stated by the following equivalent formulation:(50)θusL∞0,T;Eh2+θpsL∞0,T;L2Ω2+θqsL20,T;L2Ω2≤Ch2k∫0Tqsk2+∂tusk+12+usk+12ds+sup0≤s≤Tusk+12.
This together with the standard interpolation error estimates forξp, ξq, and ξu and the triangle inequality yields the desired estimate (36).


## 3.1. Semidiscrete H(div)-Conforming Methods with DG Formulation

The expression ofah⋅,⋅ below is similar to the statement in Section 3.1 in the work [43]. For completeness, we derive its scheme as follows.Multiplying the third equation in (4) by any v∈Vh on each element K, integrating by parts, and then summing over all elements in Th, we arrive at(20)−∑T∈Th∇⋅σ,vT=∑T∈Thσ,∇vT−∑T∈ThσnT,v∂T=f,v,where nT is outward normal of T.We first have(21)∑T∈ThσnT,v∂T=∑e∈ℰhIσne,ve+∑e∈ℰhσne,ve=∑e∈ℰhσne,ve=∑e∈ℰhI∪ℰhdσne,ve=2μ∑e∈ℰhI∪ℰhdεune,ve+λ∑e∈ℰhI∪ℰhd∇⋅une,ve+α∑e∈ℰhI∪ℰhdpne,ve=2μ∑e∈ℰhI∪ℰhdεune,ve.Here in the third line, we have usedσne=0 in ℰht, and in the last line, we have used the fact that v⋅ne=0 on e∈ℰhI∪ℰhd and the regularity of the exact solutions.In addition, it is easy to check that(22)∑T∈Thσ,∇vT=2μ∑T∈Thϵu,ϵvT+λ∑T∈Th∇⋅u,∇⋅vT−α∑T∈Th∇⋅v,pT=2μ∑T∈Thϵu,ϵvT+λ∇⋅u,∇⋅v−α∇⋅v,p.Substituting (21) and (22) into (20) yields(23)2μ∑K∈Thϵu,ϵvT−2μ∑e∈εhI∪εhdϵune,ve+λ∇⋅u,∇⋅v−α∇⋅v,p=f,v,∀v∈Vh.Letτe be the unit tangential direction to the edge e so that ne and τe form a right-hand side coordinate system; it follows that for any w, we have(24)w=w⋅nene+w⋅τeτe.Therefore,(25)ϵune⋅v=ϵune⋅nene+ϵune⋅τeτe⋅v⋅nene+v⋅τeτe=ϵune⋅nev⋅ne+ϵune⋅τev⋅τe.Combining (18) and (19) and the regularity of the exact solutions implies that(26)2μ∑e∈ℰhI∪ℰhdϵune,ve=2μ∑e∈ℰhI∪ℰhd∫eϵune⋅vdF=2μ∑e∈ℰhI∪ℰhd∫eϵune⋅nev⋅nedF+2μ∑e∈ℰhI∪ℰhd∫eϵune⋅τev⋅τedF=2μ∑e∈ℰhI∪ℰhd∫eϵune⋅τev⋅τedF=2μ∑e∈ℰhI∪ℰhdϵune⋅τe,v⋅τee.Here in the third line, we have used the fact thatv⋅ne=0 on e∈ℰhI∪ℰhd. Therefore, (23) is reduced to(27)2μ∑K∈Thϵu,ϵvT−2μ∑e∈ℰhI∪ℰhdϵune⋅τe,v⋅τee+λ∇⋅u,∇⋅v−α∇⋅v,p=f,v,∀v∈Vh.Adding some stabilized terms as in [47], we propose the following DG method for the third equation in (4):(28)ahu,v−αp,∇⋅v=f,v,with(29)ahu,v=2μ∑K∈Thϵu,ϵvT+λ∇⋅u,∇⋅v−2μ∑e∈ℰhI∪ℰhdϵune⋅τe,v⋅τee−2μ∑e∈ℰhI∪ℰhdϵvne⋅τe,u⋅τee+∑e∈ℰhI∪ℰhd2μβheu⋅τe,v⋅τee.Consequently, after integrating by parts, we find that the exact solutions of (4) satisfy(30)c0∂tp,w+α∇⋅∂tu,w+∇⋅q,w=g,w,κ−1∇⋅uq,z−p,∇⋅z=0,ahu,v−αp,∇⋅v=f,v,for any w,z,v∈Ph×Qh×Vh.Thus, the corresponding semidiscrete H(div)-conforming FEMs for (4) can be designed as follows: find ph,qh,uh∈Ph×Qh×Vh such that(31)c0∂tph,w+α∇⋅∂tuh,w+∇⋅qh,w=g,w,κ−1∇⋅uhqh,z−ph,∇⋅z=0,ahuh,v−αph,∇⋅v=f,v,for any w,z,v∈Ph×Qh×Vh, with the initial conditions given by(32)ph0=Php0,uh0=Πhu0.


## 3.2. Error Estimates

We begin by defining the mesh-dependent norm⋅h on Vh=V+Vh:(33)vh=∑K∈Th2μϵv0,K2+∑e∈ℰhI∪ℰhd2μhe−1v⋅τe0,e2+∑e∈ℰhI∪ℰhd2μheϵvne⋅τe0,e2+λ∇⋅v021/2.Then, from [43], we can obtain the following result which shows the continuity and coercivity of the bilinear ah⋅,⋅.Lemma 1.
There exists a constantCa>0 such that(34)ahx,v≤Caxhvh,∀x,v∈Vh.Moreover, for sufficiently large penalty parameterβ, it holds that(35)ahv,v≥Cbvh2,∀v∈Vh.We are now in a position to state the following error estimate, which is the main result of this section.Theorem 1.
Letp,q,u∈P×Q×V and ph,qh,uh∈Ph×Qh×Vh be the solutions of (4) and (31), respectively. Then, the following holds:(36)u−uhL∞0,T;Eh+p−phL∞0,T;L2Ω+q−qhL20,T;L2Ω≤Chk,where xL∞0,T;Eh=sup0≤s≤Txsh and C depends on the regularity of u,ut,p,pt, and q.Proof.
Subtracting (30) from (31) yields(37)c0∂tp−ph,w+α∇⋅∂tu−uh,w+∇⋅q−qh,w=0,κ−1∇⋅uq−κ−1∇⋅uhqh,z−p−ph,∇⋅z=0,ahu−uh,v−αp−ph,∇⋅v=0,for any w,z,v∈Ph×Qh×Vh.
We then split the errorp−ph into p−ph=ξp+θp with ξp=p−Php and θp=Php−ph. Similarly, q−qh=ξq+θq with ξq=q−Πhq and θq=Πhq−qh; u−uh=ξu+θu with ξu=u−Πhu and θu=Πhu−uh. Since the estimates for ξp, ξq, and ξu can be derived by the interpolation error bounds in (16), (18), and (19), it remains to estimate θp, θq, and θu. To this end, using (15) and (17), we can rewrite (37) as(38)c0∂tθp,w+α∇⋅∂tθu,w+∇⋅θq,w=0,κ−1∇⋅uq−κ−1∇⋅uhqh,z−θp,∇⋅z=0,ahθu,v−αθp,∇⋅v=−ahξu,v,for any w,z,v∈Ph×Qh×Vh.
Settingw=θp, z=θq, and v=∂tθu in the above equations, we infer that(39)c02ddtθp,θp+α∇⋅∂tθu,θp+∇⋅θq,θp=0,κ−1∇⋅uq−κ−1∇⋅uhqh,θq−θp,∇⋅θq=0,12ddtahθu,θu−αθp,∇⋅∂tθu=−ahξu,∂tθu.
Summing equation (39), integrating in time from 0 to t≤T, and noting that θp0=0, θu0=0, and(40)κ−1∇⋅uq−κ−1∇⋅uhqh=κ−1∇⋅uq−Πhq+κ−1∇⋅u−κ−1∇⋅uhΠhq+κ−1∇⋅uhΠhq−qh,we obtain(41)12ahθut,θut+c02θpt02+∫0tκ−1∇⋅uhθqs,θqsds=A1+A2+A3,with(42)A1=−∫0tκ−1∇⋅uξqs,θqsds,A2=−∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs,θqsds,A3=−∫0tahξus,∂tθusds.
For the first termA1, using (6) and Cauchy–Schwarz and Young’s inequalities, we find that(43)A1=−∫0tκ−1∇⋅uξqs,θqsds≤∫0tκ−1∇⋅uξqs0θqs0ds≤14ϵ1∫0t∥κ−1∇⋅uξqs∥02ds+ɛ1∫0tθqs02ds≤kb24ε1∫0tξqs02ds+ɛ1∫0tθqs02ds,where ɛ1 is an arbitrarily small number.
For the second termA2, with the aid of (6) and (7) and Cauchy–Schwarz and Young’s inequalities, we deduce that(44)A2=−∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs,θqsds≤∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs0θqs0ds≤14ɛ2∫0tκ−1∇⋅u−κ−1∇⋅uhΠhqs02ds+ɛ2∫0tθqs02ds≤14ɛ2∫0tκ−1∇⋅u−κ−1∇⋅uh02Πhqs∞2ds+ɛ2∫0tθqs02ds≤L4ɛ2∫0t∇⋅u−∇⋅uh02Πhqs∞2ds+ɛ2∫0tθqs02ds≤LC14ɛ2∫0tus−uhsh2Πhqs∞2ds+ɛ2∫0tθqs02ds≤LC12ɛ2∫0tus−uhsh2Πhqs∞2ds+LC12ɛ2∫0tΠhus−uhsh2Πhqs∞2ds+ɛ2∫0tθqs02ds≤LC1Cq22ɛ2∫0tus−Πhush2ds+LC1Cq22ɛ2∫0tΠhus−uhsh2ds+ɛ2∫0tθqs02ds=LC1Cq22ɛ2∫0tξush2ds+LC1Cq22ɛ2∫0tθus02dsɛ2∫0tθqs02ds,where ϵ2 is an arbitrarily small number and Cq=max0≤s≤tΠhqs∞.
To estimate the third termA3, integrating by parts and noting that θu0=0, we first obtain(45)A3=∫0tah∂tξus,θusds−ahξut,θut,which combined with (34) and Young’s inequality implies that(46)A3≤Ca2∫0t∂tξush2+θush2ds+Ca24ɛ3ξuth2+ɛ3θuth2,where ɛ3 is an arbitrarily small number.
Combining the bounds in (41)–(46) and using (6) and (35), we have(47)Cb2−ɛ3θuth2+c02θpt02+ka−ɛ1−ɛ2∫0tθqs02ds≤LC1Cq22ɛ2∫0tθush2ds+LC1Cq22ɛ2∫0tξush2ds+kb24ɛ1∫0tξqs02ds+Ca2∫0t∂tξush2ds+Ca24ɛ3ξuth2.
Now, we choose appropriateɛii=1,2,3 such that Cmin=minCb/2−ɛ3,c0/2,ka−ɛ1−ɛ2>0. Consequently, the above inequality still holds if we replace its left-hand side term by Cminθuth2+θpt02+∫0tθqs02ds. Then, dividing both sides of the above inequality by Cmin and using Gronwall’s inequality, we have(48)θuth2+θpt02+∫0tθqs02ds≤C∫0tξqs02+∂tξush2+ξush2ds+ξuth2.
Noting that the above estimate holds for all0≤t≤T and using appropriate approximation properties stated in (16), (18), and (19), we obtain(49)sup0≤s≤Tθush2+sup0≤s≤Tθps02+∫0T∫0Tθqs02ds≤Ch2k∫0Tqsk2+∂tusk+12+usk+12ds+h2ksup0≤s≤Tusk+12.
This can be stated by the following equivalent formulation:(50)θusL∞0,T;Eh2+θpsL∞0,T;L2Ω2+θqsL20,T;L2Ω2≤Ch2k∫0Tqsk2+∂tusk+12+usk+12ds+sup0≤s≤Tusk+12.
This together with the standard interpolation error estimates forξp, ξq, and ξu and the triangle inequality yields the desired estimate (36).


## 4. Error Estimates for the Fully Discrete Scheme

### 4.1. Fully Discrete Finite Element Method

LetN be a positive integer and set Δt=T/N and tn=nΔt1≤n≤N. For the function ft,x, set fn=fx,tn,∀x∈Ω. The fully discrete mixed element method with back Euler time stepping reads as follows: at each time t=tn, find phn,qhn,uhn∈Ph×Qh×Vh such that(51)c0∂t¯phn,w+α∇⋅∂t¯uhn,w+∇⋅qhn,w=gn,w,κ−1∇⋅uhnqhn,z−phn,∇⋅z=0,ahuhn,v−αphn,∇⋅v=fn,v,for any w,z,v∈Ph×Qh×Vh, where ∂t¯phn=phn−phn−1/Δt and ∂t¯uhn=uhn−uhn−1/Δt with initial condition(52)ph0=Php0,uh0=Πhu0.


### 4.2. Error Estimates

The standard Taylor expansion leads to(53)fn−fn−1Δt=∂tfn+1Δt∫tn−1tntn−1−s∂ttfsds,where ∂tfn=∂tftn.Now, we prove the following error estimate which is the main result of this section.Theorem 2.
Letp,q,u∈P×Q×V and phn,qhn,uhn∈Ph×Qh×Vh be the solutions of (4) and (51), respectively. Then, the following finite element error estimate holds:(54)max1≤n≤Nun−uhnh+max1≤n≤Npn−phn0+Δt∑n=1Nqn−qhn0≤Chk+Δt,where C depends on the regularity of u,ut,utt,p,pt,ptt, and q.Proof.
Since (30) holds for the exact solution u at any time t=tn, this together with Taylor expansion (53) yields(55)c0∂t¯pn,w+α∇⋅∂t¯un,w+∇⋅qn,w=gn,w+c0Δt∫tn−1tntn−1−s∂ttpsds,w+αΔt∫tn−1tntn−1−s∇⋅∂ttusds,w,κ−1∇⋅unqn,z−pn,∇⋅z=0,ahun,v−αpn,∇⋅v=fn,v,for every w,z,v∈Ph×Qh×Vh.
Subtracting (51) from (55), we obtain(56)c0pn−phn−pn−1−phn−1Δt,w+α∇⋅un−uhn−un−1−uhn−1Δt,w+∇⋅qn−qhn,w=c0Δt∫tn−1tntn−1−s∂ttpsds,w+αΔt∫tn−1tntn−1−s∇⋅∂ttusds,w,κ−1∇⋅unqn−κ−1∇⋅uhnqhn,z−pn−phn,∇⋅z=0,ahun−uhn,v−αpn−phn,∇⋅v=0,for every w,z,v∈Ph×Qh×Vh.
We then split the errorpn−phn into pn−phn=ξpn+θpn, with ξpn=pn−Phpn and θpn=Phpn−phn. Similarly, qn−qhn=ξqn+θqn with ξqn=qn−Πhqn and θqn=Πhqn−qhn. un−uhn=ξun+θun with ξun=un−Πhun and θun=Πhun−uhn. Since the estimates for ξpn, ξqn, and ξun can be derived by the interpolation error bounds, it leaves us to estimate θpn, θqn, and θun. To this end, using (15) and (17), we can rewrite (56) by(57)c0θpn−θpn−1Δt,w+α∇⋅θun−θun−1Δt,w+∇⋅θqn,w=c0Δt∫tn−1tntn−1−s∂ttpsds,w+αΔt∫tn−1tntn−1−s∇⋅∂ttusds,w,κ−1∇⋅unqn−κ−1∇⋅uhnqhn,z−θpn,∇⋅z=0,ahθun,v−αθpn,∇⋅v=−ahξun,v,for every w,z,v∈Ph×Qh×Vh.
Settingw=θpn, z=θqn, and v=θun−θun−1/Δt in the above equations, adding these equations, and using(58)κ−1∇⋅unqn−κ−1∇⋅uhnqhn=κ−1∇⋅unqn−Πhqn+κ−1∇⋅un−κ−1∇⋅uhnΠhqn+κ−1∇⋅uhnΠhqn−qhn,we infer that(59)ahθun,θun+c0θpn02+Δtκ−1∇⋅unθqn,θqn=ahθun,θun−1+c0θpn−1,θpn+c0∫tn−1tntn−1−spttsds,θpn+α∫tn−1tntn−1−s∇⋅uttssds,θpn−Δtκ−1∇⋅uhnξqn,θqn−1−Δtκ−1∇⋅un−κ−1∇⋅uhnΠhqn,θqn−ahξun,θun−θun−1.
To give bounds for the above error equation, we introduce the following useful inequalities:(60)ahθun,θun−1≤12ahθun−1,θun−1+ahθun,θun,(61)c0θpn−1,θpn≤c02θpn−102+θpn02.
Summing (59) from 1 to m≤N, applying (60) and (61), and noting that θu0=0 and θp0=0, we obtain(62)12ahθum,θum+c0θpm02+Δt∑n=1mκ−1∇⋅unθqn,θqn≤B1+B2+B3+B4+B5,where(63)B1=c0∑n=1m∫tn−1tntn−1−s∂ttpsds,θpn,B2=α∑n=1m∫tn−1tntn−1−s∇⋅∂ttusds,θpn,B3=−∑n=1mΔtκ−1∇⋅uhnξqn,θqn,B4=−∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn,θqn,B5=−∑n=1mahξun,θun−θun−1.
The first termB1 can be bounded by(64)B1=c0∑n=1m∫tn−1tntn−1−s∂ttpsds,θpn≤c0∑n=1m∫tn−1tntn−1−s∂ttpsds0θpn0.
Since(65)∫tn−1tntn−1−spttsds0≤Δt3/2∫tn−1tn∂ttps02ds1/2,then we further obtain(66)B1≤CΔt∑n=1mθpn02+Δt2∫0tm∂ttps02ds.
Similarly, the second termB2 can be estimated by(67)B2≤CΔt∑n=1mθpn02+Δt2∫0tm∇⋅∂ttus02ds.
Combining (6) and Cauchy–Schwarz and Young’s inequalities, we have(68)B3=−∑n=1mΔtκ−1∇⋅uhnξqn,θqn≤∑n=1mΔtθqn0κ−1∇⋅uhnξqn0≤ε4Δt∑n=1mθqn02+14ε4Δt∑n=1mκ−1∇⋅uhnξqn02≤ε4Δt∑n=1mθqn02+kb24ε4Δt∑n=1mξqn02,where ε4 is an arbitrarily small number.
Similarly, we can boundB4 by(69)B4=−∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn,θqn≤∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn0θqn0≤14ε5∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn02+ε5∑n=1mΔtθqn02≤14ε5Δt∑n=1mκ−1∇⋅un−κ−1∇⋅uhn02Πhqn∞2+ε5Δt∑n=1mθqn02≤L4ε5Δt∑n=1m∇⋅un−∇⋅uhn02Πhqn∞2+ε5Δt∑n=1mθqn02≤LC14ε5Δt∑n=1mun−uhnh2Πhqn∞2+ε5Δt∑n=1mθqn02≤LC12ε5Δt∑n=1mun−Πhunh2Πhqn∞2+LC12ε5Δt∑n=1mΠhun−uhnh2Πhqn∞2+ε5Δt∑n=1mθqn02≤LC1Cm22ε5Δt∑n=1mun−Πhunh2+LC1Cm22ε5Δt∑n=1mΠhun−uhnh2+ε5Δt∑n=1mθqn02,where ε5 is an arbitrarily small number and Cm=max1≤n≤mΠhqn∞.
To bound the last termB5, we need the following equalities:(70)∑n=1mfn−fn−1gn−1=fmgm−f0g0−∑n=1mfngn−gn−1,and(71)ξun−ξun−1Δt=∂tξun+1Δt∫tn−1tntn−1−s∂ttξusds.
Here and in the following,∂tξun=∂tξutn.
Then, using (70), (71), and (34) and Young’s inequality and noting that θu0=0, we infer that(72)B5=−∑n=1mahξun,θun−θun−1=−ahξum,θum+∑n=1mahξun−ξun−1,θun−1≤ε6θumh2+Cξumh2+Δt2∫0tm∂ttξush2ds+Δt∑n=1m∂tξunh2+Δt∑n=1m−1θunh2,where ε6 is an arbitrarily small number.
Combining the bounds above and using (6) and (35), we have(73)Cb2−ε6θumh2+c02θpm02+ka−ε4−ε5∑n=1mθqn02≤CΔt∑n=1mθpn02+θunh2+Δt2∫0tm∂ttps02ds+Δt2∫0tm∇⋅∂ttus02ds+Δt∑n=1mξqn02+Δt∑n=1mξunh2+ξumh2+Δt2∫0T∂ttξush2ds+Δt∑n=1m∂tξunh2.
Choose small enoughεii=4,5,6 to make Cmin=minCb/2−ε6,c0/2,ka−ε4−ε5>0. Then, the above inequality still holds if we replace its left-hand side term by Cminθumh2+θpm02+Δt∑n=1mθqn02. Using the discrete Gronwall’s inequality and some approximation properties and noting that (73) holds for any 1≤m≤N, we obtain(74)max1≤n≤Nθunh2+max1≤n≤Nθpn02+Δt∑n=1Nθqn02≤CΔt2∫0T∂ttps02ds+Δt2∫0T∇⋅∂ttus02ds+h2kmax1≤n≤Nqnk2+h2kmax1≤n≤Nunk+12+h2kΔt2∫0T∂ttusk+12ds+h2kmax1≤n≤N∂tunk+12.
Combining the above estimate with the interpolation error estimates forξpn, ξqn, and ξun and using the triangle inequality, we obtain the desired assertion in (54).


## 4.1. Fully Discrete Finite Element Method

LetN be a positive integer and set Δt=T/N and tn=nΔt1≤n≤N. For the function ft,x, set fn=fx,tn,∀x∈Ω. The fully discrete mixed element method with back Euler time stepping reads as follows: at each time t=tn, find phn,qhn,uhn∈Ph×Qh×Vh such that(51)c0∂t¯phn,w+α∇⋅∂t¯uhn,w+∇⋅qhn,w=gn,w,κ−1∇⋅uhnqhn,z−phn,∇⋅z=0,ahuhn,v−αphn,∇⋅v=fn,v,for any w,z,v∈Ph×Qh×Vh, where ∂t¯phn=phn−phn−1/Δt and ∂t¯uhn=uhn−uhn−1/Δt with initial condition(52)ph0=Php0,uh0=Πhu0.


## 4.2. Error Estimates

The standard Taylor expansion leads to(53)fn−fn−1Δt=∂tfn+1Δt∫tn−1tntn−1−s∂ttfsds,where ∂tfn=∂tftn.Now, we prove the following error estimate which is the main result of this section.Theorem 2.
Letp,q,u∈P×Q×V and phn,qhn,uhn∈Ph×Qh×Vh be the solutions of (4) and (51), respectively. Then, the following finite element error estimate holds:(54)max1≤n≤Nun−uhnh+max1≤n≤Npn−phn0+Δt∑n=1Nqn−qhn0≤Chk+Δt,where C depends on the regularity of u,ut,utt,p,pt,ptt, and q.Proof.
Since (30) holds for the exact solution u at any time t=tn, this together with Taylor expansion (53) yields(55)c0∂t¯pn,w+α∇⋅∂t¯un,w+∇⋅qn,w=gn,w+c0Δt∫tn−1tntn−1−s∂ttpsds,w+αΔt∫tn−1tntn−1−s∇⋅∂ttusds,w,κ−1∇⋅unqn,z−pn,∇⋅z=0,ahun,v−αpn,∇⋅v=fn,v,for every w,z,v∈Ph×Qh×Vh.
Subtracting (51) from (55), we obtain(56)c0pn−phn−pn−1−phn−1Δt,w+α∇⋅un−uhn−un−1−uhn−1Δt,w+∇⋅qn−qhn,w=c0Δt∫tn−1tntn−1−s∂ttpsds,w+αΔt∫tn−1tntn−1−s∇⋅∂ttusds,w,κ−1∇⋅unqn−κ−1∇⋅uhnqhn,z−pn−phn,∇⋅z=0,ahun−uhn,v−αpn−phn,∇⋅v=0,for every w,z,v∈Ph×Qh×Vh.
We then split the errorpn−phn into pn−phn=ξpn+θpn, with ξpn=pn−Phpn and θpn=Phpn−phn. Similarly, qn−qhn=ξqn+θqn with ξqn=qn−Πhqn and θqn=Πhqn−qhn. un−uhn=ξun+θun with ξun=un−Πhun and θun=Πhun−uhn. Since the estimates for ξpn, ξqn, and ξun can be derived by the interpolation error bounds, it leaves us to estimate θpn, θqn, and θun. To this end, using (15) and (17), we can rewrite (56) by(57)c0θpn−θpn−1Δt,w+α∇⋅θun−θun−1Δt,w+∇⋅θqn,w=c0Δt∫tn−1tntn−1−s∂ttpsds,w+αΔt∫tn−1tntn−1−s∇⋅∂ttusds,w,κ−1∇⋅unqn−κ−1∇⋅uhnqhn,z−θpn,∇⋅z=0,ahθun,v−αθpn,∇⋅v=−ahξun,v,for every w,z,v∈Ph×Qh×Vh.
Settingw=θpn, z=θqn, and v=θun−θun−1/Δt in the above equations, adding these equations, and using(58)κ−1∇⋅unqn−κ−1∇⋅uhnqhn=κ−1∇⋅unqn−Πhqn+κ−1∇⋅un−κ−1∇⋅uhnΠhqn+κ−1∇⋅uhnΠhqn−qhn,we infer that(59)ahθun,θun+c0θpn02+Δtκ−1∇⋅unθqn,θqn=ahθun,θun−1+c0θpn−1,θpn+c0∫tn−1tntn−1−spttsds,θpn+α∫tn−1tntn−1−s∇⋅uttssds,θpn−Δtκ−1∇⋅uhnξqn,θqn−1−Δtκ−1∇⋅un−κ−1∇⋅uhnΠhqn,θqn−ahξun,θun−θun−1.
To give bounds for the above error equation, we introduce the following useful inequalities:(60)ahθun,θun−1≤12ahθun−1,θun−1+ahθun,θun,(61)c0θpn−1,θpn≤c02θpn−102+θpn02.
Summing (59) from 1 to m≤N, applying (60) and (61), and noting that θu0=0 and θp0=0, we obtain(62)12ahθum,θum+c0θpm02+Δt∑n=1mκ−1∇⋅unθqn,θqn≤B1+B2+B3+B4+B5,where(63)B1=c0∑n=1m∫tn−1tntn−1−s∂ttpsds,θpn,B2=α∑n=1m∫tn−1tntn−1−s∇⋅∂ttusds,θpn,B3=−∑n=1mΔtκ−1∇⋅uhnξqn,θqn,B4=−∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn,θqn,B5=−∑n=1mahξun,θun−θun−1.
The first termB1 can be bounded by(64)B1=c0∑n=1m∫tn−1tntn−1−s∂ttpsds,θpn≤c0∑n=1m∫tn−1tntn−1−s∂ttpsds0θpn0.
Since(65)∫tn−1tntn−1−spttsds0≤Δt3/2∫tn−1tn∂ttps02ds1/2,then we further obtain(66)B1≤CΔt∑n=1mθpn02+Δt2∫0tm∂ttps02ds.
Similarly, the second termB2 can be estimated by(67)B2≤CΔt∑n=1mθpn02+Δt2∫0tm∇⋅∂ttus02ds.
Combining (6) and Cauchy–Schwarz and Young’s inequalities, we have(68)B3=−∑n=1mΔtκ−1∇⋅uhnξqn,θqn≤∑n=1mΔtθqn0κ−1∇⋅uhnξqn0≤ε4Δt∑n=1mθqn02+14ε4Δt∑n=1mκ−1∇⋅uhnξqn02≤ε4Δt∑n=1mθqn02+kb24ε4Δt∑n=1mξqn02,where ε4 is an arbitrarily small number.
Similarly, we can boundB4 by(69)B4=−∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn,θqn≤∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn0θqn0≤14ε5∑n=1mΔtκ−1∇⋅un−κ−1∇⋅uhnΠhqn02+ε5∑n=1mΔtθqn02≤14ε5Δt∑n=1mκ−1∇⋅un−κ−1∇⋅uhn02Πhqn∞2+ε5Δt∑n=1mθqn02≤L4ε5Δt∑n=1m∇⋅un−∇⋅uhn02Πhqn∞2+ε5Δt∑n=1mθqn02≤LC14ε5Δt∑n=1mun−uhnh2Πhqn∞2+ε5Δt∑n=1mθqn02≤LC12ε5Δt∑n=1mun−Πhunh2Πhqn∞2+LC12ε5Δt∑n=1mΠhun−uhnh2Πhqn∞2+ε5Δt∑n=1mθqn02≤LC1Cm22ε5Δt∑n=1mun−Πhunh2+LC1Cm22ε5Δt∑n=1mΠhun−uhnh2+ε5Δt∑n=1mθqn02,where ε5 is an arbitrarily small number and Cm=max1≤n≤mΠhqn∞.
To bound the last termB5, we need the following equalities:(70)∑n=1mfn−fn−1gn−1=fmgm−f0g0−∑n=1mfngn−gn−1,and(71)ξun−ξun−1Δt=∂tξun+1Δt∫tn−1tntn−1−s∂ttξusds.
Here and in the following,∂tξun=∂tξutn.
Then, using (70), (71), and (34) and Young’s inequality and noting that θu0=0, we infer that(72)B5=−∑n=1mahξun,θun−θun−1=−ahξum,θum+∑n=1mahξun−ξun−1,θun−1≤ε6θumh2+Cξumh2+Δt2∫0tm∂ttξush2ds+Δt∑n=1m∂tξunh2+Δt∑n=1m−1θunh2,where ε6 is an arbitrarily small number.
Combining the bounds above and using (6) and (35), we have(73)Cb2−ε6θumh2+c02θpm02+ka−ε4−ε5∑n=1mθqn02≤CΔt∑n=1mθpn02+θunh2+Δt2∫0tm∂ttps02ds+Δt2∫0tm∇⋅∂ttus02ds+Δt∑n=1mξqn02+Δt∑n=1mξunh2+ξumh2+Δt2∫0T∂ttξush2ds+Δt∑n=1m∂tξunh2.
Choose small enoughεii=4,5,6 to make Cmin=minCb/2−ε6,c0/2,ka−ε4−ε5>0. Then, the above inequality still holds if we replace its left-hand side term by Cminθumh2+θpm02+Δt∑n=1mθqn02. Using the discrete Gronwall’s inequality and some approximation properties and noting that (73) holds for any 1≤m≤N, we obtain(74)max1≤n≤Nθunh2+max1≤n≤Nθpn02+Δt∑n=1Nθqn02≤CΔt2∫0T∂ttps02ds+Δt2∫0T∇⋅∂ttus02ds+h2kmax1≤n≤Nqnk2+h2kmax1≤n≤Nunk+12+h2kΔt2∫0T∂ttusk+12ds+h2kmax1≤n≤N∂tunk+12.
Combining the above estimate with the interpolation error estimates forξpn, ξqn, and ξun and using the triangle inequality, we obtain the desired assertion in (54).


## 5. Conclusion and Future Work

We propose and analyze H(div)-conforming FEMs for a nonlinear poroelasticity model. Optimal a priori error estimates are derived for both semidiscrete and fully discrete schemes. After finite element discretizations, we shall solve a nonlinear discrete system. Thus, there is a need to design efficient iterative methods to carry out numerical simulations. One natural iterative solver is the two-grid method [48] that will be carefully investigated in the future work.



---

*Source: 9464389-2020-09-19.xml*