GUARANTEED LOCKING-FREE FINITE ELEMENT METHODS FOR BIOT’S CONSOLIDATION MODEL IN POROELASTICITY
JEONGHUN J. LEE
Abstract. We propose a new finite element method for the three-field formu- lation of Biot’s consolidation model in poroelasticity and prove a priori error estimates. Uniform-in-time error estimates of all the unknowns are obtained for both semidiscrete solutions and fully discrete solutions with the backward Euler time discretization. The novelty of our method is that the error analysis does not require the assumption that the constrained specific storage coefficient is uniformly positive. Therefore the method is guaranteed to be locking-free without the additional assumption on material parameters.
1. Introduction
The Biot’s consolidation model simultaneously describes the deformation of a saturated elastic porous medium and the viscous fluid flow inside [2]. Since this model naturally arises from classical engineering applications such as geomechanics and petrolium engineering, the study of its numerical solutions has been of great interest.
There are extensive literature on numerical schemes for Biot’s consolidation model with finite element methods. In early studies of the problem with contin- uous Galerkin finite elements, nonphysical pressure oscillations of numerical solu- tions, called poroelasticity locking, were observed for certain ranges of material parameters and small time step sizes [20, 23, 27].
In order to avoid the poroelasticity locking, various numerical methods for the problem with different formulations were considered. In a series of papers [13, 14, 15], Murad and his collaborators studied a two-field formulation of Biot’s model in incompressible porous media, with displacement and pressure as un- knowns, using mixed finite elements for Stokes equations. A discontinuous Galerkin method for the two-field formulation was studied recently [7]. A Galerkin least square method [9] was proposed for a four-field formulation of the problem, which has displacement, stress, fluid flux, and pressure as unknowns. A three-field formu- lation was studied with various couplings of continuous and discontinuous Galerkin methods, and mixed finite element methods in [16, 17, 18]. A coupling of noncon- forming and mixed finite element methods for the formulation was recently studied in [24]. For more information on previous studies we refer to [7, 10, 24] and the references therein.
In the present paper, we are interested in the three-field formulation of the prob- lem. Compared to the two-field formulation approach, the one more unknown seems to be a disadvantage but there is a reason that the three-field formulation approach may be preferred. An exact solution of Biot’s model typically has a large jump of pressure, so when continuous finite elements are used for the numerical pressure,
1
there is an overshoot for this large jump. Unfortunately, both the Taylor–Hood and MINI elements, the most popular mixed finite element families for Stokes problems, use continuous finite elements for pressure, so the overshooting is inevitable. How- ever, in the methods for three-field formulation, discontinuous finite elements are eligible for the pressure, so there is no overshooting.
To ensure that a method is locking-free, a uniform-in-time error bound of the pressure is needed. In addition, the heuristic analysis and numerical experiments in [19] insinuate that the locking typically occurs when the constrained specific storage coefficientc0≥0 of the model is very close to 0. Thus the uniform-in-time error bound should be robust for very small or even vanishing c0. However, the assumption that c0 is uniformly positive is essential to obtain a uniform-in-time error bound of pressure in the error analysis of aforementioned methods for the three-field formulation [16, 17, 24]. Although numerical tests of those methods suggest that they are good candidates to avoid the poroelasticity locking, their error analysis is not satisfactory from the mathematical point of view.
The goal of this paper is to develop a new locking-free finite element method for the three-field formulation of Biot’s model. To confirm that the method is locking- free we show a priori error estimates of uniform-in-time errors of all the unknowns without the assumption thatc0is uniformly positive. To the best of our knowledge, this result has not been achieved previously.
The paper is organized as follows. In section 2 we introduce notations and preliminaries of the problem. In section 3 we present our numerical method for spatial discretization of the problem and show a priori error analysis of semidiscrete solutions. In section 4 we show error analysis of fully discrete solutions with the backward Euler time discretization. Finally, we conclude with remarks on extending the method to higher order methods and to rectangular meshes.
2. Preliminaries
2.1. Notations. Let Ω be a bounded Lipschitz domain inRn withn= 2 or 3. For a nonnegative integerm,Hm(Ω),Hm(Ω;Rn) denote the standardRandRn-valued Sobolev spaces based on the L2 norm (H0 =L2). For a set G ⊂ Ω, k · km,G is theHmnorm on G. IfG= Ω, then we simply use k · km. The set of functions in L2(Ω;Rn) whose divergence are inL2(Ω), and the corresponding norm are denoted byH(div,Ω) andk · kdiv.
For a reflexive Banach space X and 0 < T0 < ∞, C0([0, T0];X) denotes the set of functionsf : [0, T0]→ X which are continuous in t∈[0, T0]. For an integer m≥1 we define
Cm([0, T0];X) ={f|∂if /∂ti ∈C0([0, T0];X),0≤i≤m},
where ∂if /∂ti is the i-th time derivative in the sense of the Fr´echet derivative in X (see e.g., [25]). For a functionf : [a, b]→ X, we define the space-time norm
kfkLr([a,b];X)=
Rb
akfkrXds1/r
, 1≤r <∞, esssupt∈[a,b]kfkX, r=∞.
If the time interval is fixed as [0, T0], then we useLrX to denoteLr([0, T0];X) for simplicity. We define the space-time Sobolev spacesWk,r([0, T0];X) for nonnegative integerkand 1≤r≤ ∞as the closure ofCk([0, T0];X) with the normkfkWk,rX =
Pk
i=0k∂if /∂tikLrX. The Sobolev embedding [6] gives Wk+1,1X ,→Wk,∞X. (2.1)
We adopt a convention that kf, gkX = kfkX +kgkX for the norm of a Banach spaceX. For two Banach spacesX andY, and forf ∈ X ∩ Y,kfkX ∩Y will stand forkfkX+kfkY. For simplicity of notations, ˙f, ¨f, ...
f will be used to denote time derivatives∂f /∂t,∂2f /∂t2,∂3f /∂t3.
A shape-regular triangulation of Ω will be denoted by Th for which h is the maximum diameter of triangles (or tetrahedra) andEh is the corresponding set of edges (faces), respectively. The interior edges/facesEh◦ is the set{E∈ Eh|E⊂Ω}.
ForE∈ Eh and functionsf,g:Eh→Rn we define hf,giE=
Z
E
f·gds, hf,gi= X
E∈Eh
hf,giE.
ForE∈ Eh and an element-wiseH1functionv,JvKis defined by JvK|E=
(the jump ofv onE, ifE∈ Eh◦,
v, ifE⊂∂Ω.
For an integerk≥0, and G⊂Rn,Pk(G) is the space of polynomials defined onG of degree≤k. We use Pk(Th) to denote the space of piecewise polynomials on Th
of degree≤k. For a vector spaceX, we usePk(G;X) andPk(Th;X) to denote the space ofX-valued polynomials with same conditions.
Throughout this paper we use X . Y to denote the inequality X ≤ cY with a generic constant c >0 which is independent of the mesh size, and X ∼Y will stand for X .Y and Y .X. If needed, we will use c to denote generic positive constants in inequalities and it can be a different constant in every line.
2.2. The Biot’s consolidation model. In this subsection we review the Biot’s consolidation model in poroelasticity. In an elastic porous medium saturated with a fluid, fluid flow and deformation of the porous medium are intimately related and their simultaneous behaviors are described by Biot’s model.
Throughout this paper we restrict our interest to quasistatic consolidation prob- lems. In other words, we assume that the consolidation process is slow and the acceleration term is ignored. In our description of the model,uis the displacement of the porous medium, pis the fluid pressure, f is the body force, and g is the source/sink density function of the fluid. The governing equations of the model are
−divC(u) +α∇p=f, (2.2)
c0p˙+αdiv ˙u−div(κ e
∇p) =g, (2.3)
where C is the elastic stiffness tensor, c0 ≥ 0 is the constrained specific storage coefficient, κ
e
is the hydraulic conductivity tensor, and α > 0 is the Biot–Willis constant which is close to 1. In (2.2), the div is the row-wise divergence of the Rn×n-valued functionC(u).
For isotropic elastic porous media, the elasticity tensorC has the form Cτ
e
= 2µτ e
+λtr(τ e
)I e
, τ
e
∈L2(Ω;Rn×nsym), where the constants µ, λ >0 are Lam´e coefficients, I
e
is the identity matrix, and Rn×nsym is the space of symmetricn×nmatrices. We assume thatµ,λare bounded
from above and below. The coefficientc0≥0 is determined by the permeability of the porous medium, and the bulk moduli of the solid and the fluid. The hydraulic conductivity tensor κ
e
is defined by the permeability tensor of the solid divided by the fluid viscosity and it is positive definite. We assume that κ
e
is uniformly bounded from above and below. For derivation of these equations from physical modeling, we refer to standard porous media references, for instance, [1].
In order to be a well-posed problem, the equations (2.2–2.3) need appropriate boundary and initial conditions. We assume that there are two partitions of∂Ω,
∂Ω = Γp∪Γf, ∂Ω = Γd∪Γt,
with |Γp|,|Γd| >0, i.e., the n−1-dimensional measure of Γp and Γd are positive.
Boundary conditions are given by
p(t) = 0 on Γp, −κ e
∇p(t)·n= 0 on Γf, u(t) = 0 on Γd, σ
e
(t)n= 0 on Γt, (2.4)
for all t ∈ [0, T0], in which nis the outward unit normal vector field on ∂Ω and σ
e
(t) := C(u(t))−αp(t)I e
. Here we only consider this homogeneous boundary condition for simplicity but our method, which will be introduced later, can be extended readily to problems with inhomogeneous boundary conditions. We also assume that given initial datap(0),u(0) and initial body forcef(0) satisfy (2.2).
Regularity of the solutions of the problem (2.2–2.3) with suitable boundary conditions were thoroughly studied in [21]. When we claim a priori error estimates we assume that exact solutions are sufficiently regular to obtain the claimed error bounds. For simplicity, we assume that α= 1 andκ
e
=I e
in the rest of the paper.
However, we assume thatc0≥0 is only bounded from above.
2.3. Variational formulation. Assuming α= 1,κ e
=I e
, and introducing a new unknownz:=κ
e
∇pin (2.2–2.3), we have
−divC(u) +∇p=f, (2.5)
z− ∇p= 0, (2.6)
c0p˙+ div ˙u−divz=g.
(2.7) Let
ΣΓd={u∈H1(Ω;Rn)|u|Γd= 0}, VΓf ={z∈H(div,Ω)|z·n|Γf = 0},
W =L2(Ω).
and define a bilinear form
a(u,v) = (C(u), (v)) = 2µ((u), (v)) +λ(divu,divv), u,v∈H1(Ω;Rn).
Then a variational formulation of (2.5–2.7) with boundary conditions (2.4) is to seek (u, p)∈C1([0, T0]; ΣΓd×W) andz∈C0([0, T0];VΓf) such that
a(u,v)−(p,divv) = (f,v), v∈ΣΓd, (2.8)
(z,w) + (p,divw) = 0, w∈VΓf, (2.9)
(c0p, q) + (div ˙˙ u, q)−(divz, q) = (g, q), q∈W.
(2.10)
2.4. Finite element spaces. For discretization of (2.8–2.10) we need three finite element spaces Σh,Vh,Wh for unknownsu,z,p, respectively.
We defineVh and Wh as the lowest order Raviart–Thomas–Ned´el´ec space and piecewise constant finite element space
Vh={w∈VΓf|w|T ∈(P0(T)n+xP0(T)), ∀T ∈ Th}, Wh={q∈L2(Ω)|q|T ∈ P0(T), ∀T ∈ Th}.
Here xis the vector field (x1 x2)T and (x1 x2 x3)T in two and three dimensions, respectively. Let ΠRTh be the canonical Raviart–Thomas interpolation operator into Vh, and Qh be the orthogonal L2 projection into Wh. It is well-known that, for w∈H1(Ω;Rn) andq∈H1(Ω),
divVh=Wh, div ΠRTh w=Qhdivw, (2.11)
kw−ΠRTh wk0.hkwk1, kq−Qhqk0.hkqk1, (2.12)
hold. Furthermore, there exists aβ >0, independent of the mesh size, such that
06=q∈Winf h
sup
06=w∈Vh
(q,divw)
kqk0kwkdiv ≥β >0.
(2.13)
For Σh we use vector-valued nonconforming H1 elements. In two dimensions we use the Mardal–Tai–Winther element [11] which has local shape functions on T ∈ Th as
ΣT ={v ∈(P3(T))2|divv∈ P0(T), v·n|E ∈ P1(E), E⊂∂T}, and DOFs
v 7→
Z
E
v·nEr ds, ∀r∈ P1(E), v7→
Z
E
v·tEds,
in whichnE andtE are the unit normal and tangential vectors on an edgeEofT. Let Σhbe the Mardal–Tai–Winther element with an additional condition that
all DOFs associated to the edges in Γdvanish.
(2.14)
By definition one can see that Σh⊂H(div,Ω) and div Σh⊂Wh. In three dimen- sions we use the element developed in [22] which is a three dimensional analogue of the Mardal–Tai–Winther element.
The rest of this section will be devoted to present properties of Σh which will be important for well-posedness and the a priori error analysis of our methods. Let us first define a discrete semi-norm for element-wiseH1 functions by
kvk21,h= X
T∈Th
k∇vk20,T.
By the Poincar´e inequality this is in fact a norm on ΣΓd. By (2.14), the interelement continuity of Σh, and a discrete Poincar´e inequality [4],k · k1,his a norm on Σh as well. Let
ΣΓd+ Σh={v|v=v1+v2 forv1∈ΣΓd,v2∈Σh}.
Then we are able to prove the following discrete Korn’s inequality.
Lemma 2.1. Forv∈ΣΓd+ Σh and the element-wise symmetric gradienth, kvk1,h∼ kh(v)k0.
Proof. By the definition ofk · k1,h, kh(v)k0≤ kvk1,his obvious.
Proof of the other direction is essentially same to the proof of Theorem 3.1 in [12], so we only sketch it. The inequality (1.12) in [5] gives
kvk21,h.kh(v)k20+ X
E∈Eh,E⊂Ω
h−1E kΠEJvKk20,E
+ sup
m∈RM(Ω)
kmk0,Γd=1,R
Γdmds=0
Z
Γd
v·mds 2
,
where ΠE is theL2projection toP1(E;R2) andRM(Ω) is the space of rigid body motions on Ω. Actually, the last term in the above vanishes becausev ∈ΣΓd+ Σh. In the second term ΠEJvK can be reduced to the jumps projected to the space of traces of rigid body motions onE under the help of kh(v)k20 as in [12]. However, there is no jumps in traces of rigid body motions due to the interelement continuity
of Σh, so the desired inequality follows.
It is proven in [11] that there exists an interpolation Πh : H2(Ω;Rn) → Σh, satisfying
div Πhv=Qhdivv, (2.15)
kv−Πhvk0.hkvk1, kΠhvk1,h.kvk1, kΠhv−vk1,h.hkvk2. (2.16)
Finally, there exists aβ0>0 such that
06=q∈Winf h sup
06=v∈Σh
(q,divv)
kqk0kvk1,h ≥β0>0.
(2.17)
3. Error analysis of semidiscrete solutions
In this section we discuss well-posedness of the semidiscrete problems of (2.5–
2.7) and the a priori error analysis of the semidiscrete solutions. In the rest of this paper Σh×Vh×Whwill always be the finite element spaces introduced in the previous section.
3.1. Semidiscrete problem and its well-posedness. The semidiscrete problem of (2.8–2.10) is to seek (uh,zh, ph) : [0, T0]→Σh×Vh×Wh such that
ah(uh,v)−(ph,divv) = (f,v), v∈Σh, (3.1)
(zh,w) + (ph,divw) = 0, w∈Vh, (3.2)
(c0p˙h, q) + (div ˙uh, q)−(divzh, q) = (g, q), q∈Wh, (3.3)
where
ah(u,v) = 2µ(h(u), h(v)) +λ(divu,divv), u,v∈Σh.
For the existence and uniqueness of solutions of (3.1–3.3), we first point out that this is not a system of ordinary differential equations (ODE) but a differential algebraic equation (DAE). Thus the theory of ODEs for existence of solutions is not available.
Unfortunately, there is a difficulty adopting a standard DAE theory to prove existence of solutions for the problem. To see it we review the basic theory of linear DAE problems. Let C be the complex field and m be a positive integer.
For E0, E1 ∈ Cm×m, F ∈ C0([0,∞);Cm), X0 ∈ Cm, a linear DAE is to seek X ∈C1([0,∞),Cm) such that
E0X˙ +E1X=F, X(0) =X0.
In general E0 can be singular. In standard DAE theory, this DAE is calledwell- posed if the matrix pencil E0+λE1 is regular for some λ ∈ C. However, this definition is easy to mislead because it doesnotguarantee the existence of solutions for arbitrary initial data. It is known that a well-posed DAE has a unique solution if given initial dataX0 is compatible, i.e., the initial data satisfies some intrinsic algebraic equations of the DAE [3]. To derive those intrinsic algebraic equations from a DAE, we need to know the Jordan canonical form of the augmented matrix [E0E1], and it is certainly impractical ifmis large. Another difficulty arises from the variance of c0 because finding the Jordan canonical form of the augmented matrix changes drastically whenc0 changes. Therefore we will prove the existence and uniqueness of solutions of (3.1–3.3) directly.
Theorem 3.1. For initial data ph(0)∈Wh and givenf ∈C1([0, T0];L2(Ω;Rn)), g∈C0([0, T0];L2(Ω)) there exists a unique solution(uh,zh, ph)∈C1([0, T0]; Σh× Vh×Wh)of(3.1–3.3).
Proof. We show that (3.1–3.3) is equivalent to an ODE system and check well- posedness of the ODE system.
Let {φi}, {ψi}, {χi} be bases of Σh, Vh, andWh, respectively. We use Auu, Azz,App,Bup,Bzp to denote the matrices whose (i, j)-entries are
ah(φj, φi), (ψj, ψi), (c0χj, χi), (χi,divφj), (χi,divψj), respectively. We writeuh=P
iαiφi,zh=P
iβiψi,ph=P
iγiχi,Phf =P
iζiψi, Qhg=P
iξiχi with (time-dependent) coefficient vectorsα,β,γ,ζ,ξ, wherePhis theL2 projection intoVh. Then (3.1–3.3) can be written as a matrix equation
0 0 0
0 0 0
BTup 0 App
˙ α β˙
˙ γ
+
Auu 0 −Bup
0 Azz Bzp
0 −BTzp 0
α β γ
=
ζ 0 ξ
. (3.4)
It is obvious thatAzzis symmetric positive definite andApp is symmetric positive semidefinite becausec0≥0. By Lemma 2.1,Auuis also symmetric positive definite.
The first and second rows of (3.4) give
α=A−1uuBupγ+A−1uuζ, β=−A−1zzBzpγ.
(3.5)
The third row of (3.4) givesBTupα˙ +Appγ˙ −BTzpβ =ξ. Substituting α and β in this equation using (3.5), one obtains
(BTupA−1uuBup+App) ˙γ+BTzpA−1zzBzpγ=−BTupA−1uuζ˙ +ξ,
which is an ODE system of γ. The inf-sup condition (2.16) implies that the ma- trix Bup is injective. Since BTupA−1uuBup is positive definite and App is positive semidefinite, their sum is positive definite. Thus the above system has a unique solution γ for given initial data γ(0) by a standard ODE theory. Furthermore, BTupA−1uuBup+App is still invertible asApp→0, so this well-posedness is not influ- enced by small or even vanishingc0. Note thatα,βare uniquely determined from γ by (3.5). Now the assertion follows from equivalence of (3.4) and (3.1–3.3).
In this theorem only initial dataph(0) is given butuh(0),zh(0) are determined by (3.1) and (3.2), which we will call compatibility conditions.
Definition 3.2. A triple(v0,w0, q0)∈Σh×Vh×Wh is called a compatible initial data of (3.1–3.3) if
ah(v0,v)−(q0,divv) = (f(0),v), v ∈Σh, (3.6)
(w0,w) + (q0,divw) = 0, w∈Vh, (3.7)
hold.
If the backward Euler scheme is used for time discretization, then the compat- ibility of initial data may not be crucial because the compatibility conditions are parts of the equations, so the numerical solution after one time step satisfies the compatibility conditions. However, incompatibility of initial data in DAE problems may generate a spurious numerical solution with the Crank–Nicolson scheme even if it is a very stable time discretization scheme in general.
3.2. Error analysis. We state and prove main results of semidiscrete error anal- ysis. We denote the discreteH1space with the normk · k1,hbyHh1.
Theorem 3.3. Suppose that (u,z, p) is a solution of (2.5–2.7) which is suffi- ciently regular and (uh,zh, ph) is a solution of (3.1–3.3) with compatible initial data (uh(0),zh(0), ph(0)) such that
ku(0)−uh(0)k1,h+kp(0)−ph(0)k0.h(ku(0)k2+kp(0)k1), kz(0)−zh(0)k0.hkz(0)k1.
(3.8) Then
ku−uhkL∞Hh1+kp−phkL∞L2
(3.9)
.hmax{kukW1,1H2+kpkW1,1H1,kukL∞H2+kpkL2H2+kpkL∞H1}.
and
kz−zhkL∞L2 .hmax{kpkW1,1H2,kpkL∞H2+kukW1,2H2+kpkW1,2H1}.
Remark 3.4. Compatible initial data satisfying (3.8) can be found by setting ph(0) =Qhp(0) and findinguh(0) andzh(0)using(3.6–3.7).
For the proof of Theorem 3.3 we need some preliminary results.
Lemma 3.5. Forτ e
∈H1(Ω;Rn×n) andv∈Σh define Eh(τ
e
,v) := X
E∈Eh
hτ e
nE,JvKiE.
Then
|Eh(τ e
,v)|.hkτ e
k1kvk1,h.
Proof. Let E ∈ Eh be an interior edge/face in Ω and T+, T− be two distinct triangles/tetrahedra sharing E as the common boundary. The shape regularity ofTh yieldshT− ∼hT+ ∼hE where hT−,hT+,hE are the diameters ofT−,T+,E, respectively.
SinceJvK|E is perpendicular toP0(E;Rn), hτ
e
nE,JvKiE=hτ e
nE−c,JvKiE≤ inf
c∈Rnkτ e
nE−ck0,EkJvKk0,E.
By a standard scaling argument and shape regularity of meshes one can see that kJvKk0,E .hE12k∇vk0,T+∪T−,
c∈infRn
kτ e
nE−ck0,E .h
1 2
Ek∇τk0,T+∪T−.
For an edge/face E on boundary and the triangle T ∈ Th containing E in its boundary, a similar argument gives
kJvKk0,E .h
1 2
Ek∇vk0,T, inf
c∈Rnkτ e
nE−ck0,E .hE12k∇τk0,T.
By the Cauchy–Schwarz inequality and the above results, we have Eh(τ
e
,v).hkτ e
k1kvk1,h.
as desired.
The following simple lemma will be useful in our error analysis.
Lemma 3.6. Suppose that A, B, C, D >0 satisfy A2+B2≤CA+D.
Then either A+B≤4C orA+B ≤2√
D holds.
Proof. Since eitherCA≤D orCA≥D is true, one of the followings holds.
A2+B2≤2D, A2+B2≤2CA.
(3.10)
If the first inequality in (3.10) holds, then
(A+B)2≤2(A2+B2)≤4D, which impliesA+B≤2√
D.
Suppose that the second inequality in (3.10) holds. IfB≥A, then dividing the inequality byAgives
A+B ≤(A2+B2)/A≤2C.
If B ≤A, then dividing the second inequality in (3.10) by A givesA ≤2C, and
therefore we haveA+B≤2A≤4C as desired.
Now we are ready to prove Theorem 3.3.
of Theorem 3.3. We denote the errors by
eu=u−uh, ez=z−zh, ep=p−ph.
From the left-hand side of (2.5), through the integration by parts, we obtain
−(divC(u) +∇p,v) = (C(u), h(v))−(p,divv) +Eh(σ e
,v), v ∈Σh. The difference of this and (3.1) gives
ah(eu,v)−(ep,divv) =−Eh(σ e
,v), v∈Σh. (3.11)
In addition, the differences of (2.9–2.10) and (3.2–3.3) yield
(ez,w) + (ep,divw) = 0, w∈Vh, (3.12)
(c0e˙p, q) + (div ˙eu, q)−(divez, q) = 0, q∈Wh. (3.13)
For the error analysis we split the errors into
eu=eIu+eAu := (u−Πhu) + (Πhu−uh), (3.14)
ez=eIz+eAz := (z−ΠRTh z) + (ΠRTh z−zh), (3.15)
ep=eIp+eAp := (p−Qhp) + (Qhp−ph).
(3.16)
By (2.12) and (2.16)
keIu(t)k1,h.hku(t)k2, keIz(t)k0.hkz(t)k1, keIp(t)k0.hkp(t)k1. (3.17)
Similar inequalities also hold for the time derivatives ofeIu,eIz,eIp. These inequali- ties, (3.8), and the triangle inequality give
keAu(0)k1,h+keAp(0)k0.h(ku(0)k2+kp(0)k1) keAz(0)k0.hkz(0)k1.
(3.18)
Furthermore, by the definitions ofeIu,eIz, eIpand the properties in (2.11), div ˙eIu⊥Wh, diveIz⊥Wh, eIp⊥div Σh.
(3.19)
Rewriting (3.11–3.13) using (3.14–3.16) and the orthogonalities in (3.19), we have ah(eAu,v)−(eAp,divv) =−ah(eIu,v)−Eh(σ
e ,v), (3.20)
(eAz,w) + (eAp,divw) =−(eIz,w), (3.21)
(c0e˙Ap, q) + (div ˙eAu, q)−(diveAz, q) =−(c0e˙Ip, q), (3.22)
for any (v,w, q)∈Σh×Vh×Wh. The rest of proof consists of three parts: estimates ofku−uhkL∞Hh1,kz−zhkL∞L2, andkp−phkL∞L2.
Estimate of ku−uhkL∞H1h : LetX(t), Y(t) be X2(t) =keAu(t)k2a,h+keAp(t)k2c0, Y2(t) =
Z t 0
keAzk0ds,
where kvk2a,h := ah(v,v) and kqk2c
0 := (c0q, q). Since k · k1,h ∼ k · ka,h, we will estimate max0≤t≤T0X(t) instead of max0≤t≤T0keAu(t)k.
If we take v = ˙eAu, w = eAz, q = eAp in (3.20–3.22), and add those equations together, then
1 2
d
dt(X(t))2+keAz(t)k20=−ah(eIu,e˙Au)−Eh(σ e
,e˙Au)−(eIz, eAz)−(c0e˙Ip, eAp).
Integrating this from 0 tot in time, (X(t))2+ 2(Y(t))2
(3.23)
= (X(0))2+ 2 Z t
0
−ah(eIu,e˙Au)−Eh(σ e
,e˙Au)−(eIz, eAz)−(c0e˙Ip, eAp) ds,
=: (X(0))2+ Φ1(t) + Φ2(t) + Φ3(t) + Φ4(t).
Defining ¯tby
X(¯t) = sup
0≤t≤T0
X(t), it suffices to estimateX(¯t). Moreover, since
either Φ1(¯t) + Φ2(¯t) + Φ4(¯t)≤Φ3(¯t) or Φ1(¯t) + Φ2(¯t) + Φ4(¯t)≥Φ3(¯t),
is always true, we have, from (3.23), that
either (X(¯t))2+ 2(Y(¯t))2≤(X(0))2+ 2Φ3(¯t), (3.24)
or (X(¯t))2+ 2(Y(¯t))2≤(X(0))2+ 2(Φ1(¯t) + Φ2(¯t) + Φ4(¯t)).
(3.25)
Thus it is enough to prove an estimate ofX(¯t) for these two cases.
Case I : Suppose that (3.24) is true. Note that 2Φ2(¯t) =−4
Z t¯ 0
(eIz, eAz)ds≤2 Z t¯
0
keIzk20ds+ 2(Y(¯t))2,
by Young’s inequality and the definition ofY(¯t). Combining it with (3.24) gives X2(¯t)≤X2(0) + 2
Z ¯t 0
keIzk20ds≤X(0)2+ch2kzk2L2H1. EstimatingX(0) using (3.18) and taking square roots,
X(¯t).h(ku(0)k2+kp(0)k1+kzkL2H1) (3.26)
.h(ku(0)k2+kp(0)k1+kpkL2H2).
Case II: Suppose that (3.25) is true. To estimate Φ1(¯t), we use the integration by parts in time, which gives
Φ1(¯t) =−2 Z ¯t
0
ah(eIu,e˙Au)ds, (3.27)
= 2ah(eIu(¯t), eAu(¯t))−2ah(eIu(0), eAu(0)) + 2 Z t¯
0
ah( ˙eIu, eAu)ds
≤2 keIu(¯t)ka,h+keIu(0)ka,h+ Z t¯
0
ke˙Iu(0)ka,hds
! X(¯t)
≤chkukW1,1H2X(¯t).
Similarly,
Φ2(¯t) =−2 Z ¯t
0
Eh(σ e
,e˙Au)ds (3.28)
= 2
−Eh(σ e
(¯t), eAu(¯t)) +Eh(σ e
(0), eAu(0)) + 2
Z ¯t 0
Eh( ˙σ e
, eAu)ds
≤ch kσ e
(¯t)k1+kσ e
(0)k1+ Z t¯
0
kσ˙ e
k1ds
! X(¯t)
≤chkσ e
kW1,1H1X(¯t).
For Φ4(¯t),
Φ4(¯t) =−2 Z ¯t
0
( ˙eIp, eAp)c0ds≤2 Z ¯t
0
ke˙Ipkc0dsX(¯t)≤chkpkW1,1H1X(¯t).
(3.29)
Combining (3.25), (3.27), (3.28), and (3.29), we have an inequality of the form in Lemma 3.6 withA=X(¯t),B =√
2Y(¯t),D=X(0)2, and C=ch(kσ
e
, pkW1,1H1+kukW1,1H2).
Note that kσ e
kW1,1H1 . kukW1,1H2 +kpkW1,1H1 because σ e
= C(u)−pI e
. By Lemma 3.6 and (3.18)
X(¯t).ch(kσ e
, pkW1,1H1+kukW1,1H2).h(kpkW1,1H1+kukW1,1H2) or X(¯t).h(ku(0)k2+kp(0)k1).
(3.30)
Now we complete the proof of the estimate of ku−uhkL∞Hh1. By the triangle inequality and the inequalitykeAu(t)k1,h.X(¯t),
ku(t)−uh(t)k1,h≤ keIu(t)k1,h+keAu(t)k1,h.keIu(t)k1,h+X(¯t).
By (3.17), the triangle inequality, (3.26), (3.30), and (2.1), ku(t)−uh(t)k1,h
.hmax{kukW1,1H2+kpkW1,1H1,kukL∞H2+kpkL2H2+kpkL∞H1}, for any 0≤t≤T0, which is the estimate for ku−uhkL∞Hh1 in (3.9).
Estimate ofkp−phkL∞L2: By the inf-sup condition (2.13), for any 06=q∈Wh, there exists av∈Σh such that
(divv, q0) = (q, q0), ∀q0 ∈Wh, kvk1,h.kqk0. If we use thisv in (3.20) withq=eAp(t), then
keAp(t)k20=ah(u(t)−uh(t),v) +Eh(σ e
(t),v) .(ku(t)−uh(t)k1,h+hkσ
e
(t)k1)kvk1,h,
where the last inequality is due to the estimate ofku−uhkL∞H1h. By the triangle inequality, the above estimate, (3.17), and the estimate ofku−uhkL∞H1h,
kp−phkL∞L2
.hmax{kukW1,1H2+kpkW1,1H1,kukL∞H2+kpkL2H2+kpkL∞H1}, so (3.9) is proven.
Estimate of kz−zhkL∞L2 : The time derivatives of the equations (3.20–3.21) give
ah( ˙eAu,v)−( ˙eAp,divv) =−ah( ˙eIu,v)−Eh( ˙σ e
,v), v∈Σh, ( ˙eAz,w) + ( ˙eAp,divw) =−( ˙eIz,w), w∈Vh,
Taking v = ˙eAu, w = eAz in the above, taking q = ˙eAp in (3.22), and adding these three equations together yield
1 2
d
dtkeAzk20+ke˙Auk2a,h+ke˙Apk2c0 =−ah( ˙eIu,e˙Au)−Eh( ˙σ e
,e˙Au)−( ˙eIz, eAz) + (c0e˙Ip,e˙Ap).
LetkeAz(¯t)k0= max0≤t≤T0keAz(t)k0. By integrating the above from 0 to ¯t, keAz(¯t)k20+
Z ¯t 0
{ke˙Auk2a,h+ke˙Apk2c
0}ds
=keAz(0)k20+ Z ¯t
0
{−ah( ˙eIu,e˙Au)−Eh( ˙σ e
,e˙Au)−( ˙eIz, eAz) + (c0e˙Ip,e˙Ap)}ds
=:keAz(0)k20+ Ψ1(¯t) + Ψ2(¯t) + Ψ3(¯t) + Ψ4(¯t).
If Ψ1(¯t) + Ψ2(¯t) + Ψ4(¯t)≤Ψ3(¯t), then
keAz(¯t)k20≤ keAz(0)k20+ 2Ψ3(¯t).
Using the Cauchy–Schwarz inequality and dividing both sides bykeAz(¯t)k0 yield keAz(¯t)k0≤ keAz(0)k0+ 2
Z t¯ 0
ke˙Izk0ds.
If we use (3.18) and (3.17) to estimatekeAz(0)k0 and the integral term, then keAz(¯t)k0.h(kz(0)k1+kzkW1,1H1).hkzkW1,1H1,
(3.31)
where the last inequality is due to (2.1).
On the other hand, if Ψ1(¯t) + Ψ2(¯t) + Ψ4(¯t)≥Ψ3(¯t), then keAz(¯t)k20+
Z t¯ 0
{ke˙Auk2a,h+ke˙Apk2c0}ds≤ keAz(0)k20+ 2(Ψ1(¯t) + Ψ2(¯t) + Ψ4(¯t)).
The Cauchy–Schwarz and Young’s inequalities yield
keAz(¯t)k20≤ keAz(0)k20+ch2(kuk2W1,2([0,¯t];H2)+kσ e
, pk2W1,2([0,¯t];H1)) (3.32)
.h2(kz(0)k21+kuk2W1,2H2+kpk2W1,2H1).
By the triangle inequality, (3.17) and estimates (3.31), (3.32), kz−zhkL∞L2 .keIzkL∞L2+keAzkL∞L2
.hmax{kzkW1,1H1,kzkL∞H1+kukW1,2H2+kpkW1,2H1}.
4. Error analysis of fully discrete solutions
In this section we consider the error analysis of the fully discrete solutions with the backward Euler time discretization. As in our error analysis for the semidiscrete solutions we do not use Gr¨onwall’s inequality, so our error bounds do not contain exponentially growing factors.
4.1. Fully-discrete problem and its well-posedness. Let ∆t >0 be the time step size such that T0 =N∆t for an integerN, andtj =j∆tfor j= 0,1,· · · , N.
For a continuous functionf on [0, T0], we definefj =f(tj). For a sequence{fj}j≥0, define
∂tfj+1= fj+1−fj
∆t . (4.1)
Suppose that (U0,Z0, P0) is a numerical initial data satisfying (3.8) but not nec- essarily compatible. In the backward Euler scheme, (Uj+1,Zj+1, Pj+1), the nu- merical solution at the (j+ 1)-th time step is defined inductively by
ah(Uj+1,v)−(Pj+1,divv) = (fj+1,v), (4.2)
(Zj+1,w) + (Pj+1,divw) = 0, (4.3)
(c0∂tPj+1, q) + (div(∂tUj+1), q)−(divZj+1, q) = (gj+1, q), (4.4)
for (v,w, q)∈Σh×Vh×Whand j≥0.
In order to prove that the fully discrete solution is well-defined we show that (Uj+1,Zj+1, Pj+1) is uniquely determined by the linear system (4.2–4.4) when Uj, Pj,fj+1, gj+1 are given. Rewriting (4.2–4.4),
ah(Uj+1,v)−(Pj+1,divv) = (fj+1,v), (Zj+1,w) + (Pj+1,divw) = 0,
(c0Pj+1, q) + (divUj+1, q)−∆t(divZj+1, q) = (c0Pj, q) + (divUj, q) + ∆t(gj+1, q),
for (v,w, q)∈Σh×Vh×Wh. RegardingUj+1,Zj+1, Pj+1as unknowns, the above is a system of linear equations with the same number of equations and unknowns.
Suppose that Uj = Zj = Pj = fj+1 = gj+1 = 0 and we want to show that Uj+1=Zj+1=Pj+1= 0. If we take v=Uj+1,w=Zj+1,q=Pj+1 and add all equations together, then we have
0≤ah(Uj+1,Uj+1) + (Zj+1,Zj+1) + (c0Pj+1, Pj+1) = 0,
which impliesUj+1=Zj+1= 0. Note that (4.2) is now given as (Pj+1,divv) = 0 for any v ∈Σh. ThenPj+1 = 0 because div Σh = Wh. Hence the fully discrete solution is well-defined.
4.2. Error analysis. Let us split the error (uj−Uj,zj−Zj, pj−Pj) as uj−Uj= (uj−Πhuj) + (Πhuj−uj) =:eI,ju +θuj,
(4.5)
zj−Zj= (zj−ΠRTh zj) + (ΠRTh zj−Zj) =:eI,jz +θzj, (4.6)
pj−Pj= (pj−Qhpj) + (Qhpj−Pj) =:eI,jp +θjp. (4.7)
Applying a similar argument used to obtain (3.20–3.22) yields
ah(θuj+1,v)−(θj+1p ,divv) =−ah(eI,j+1u ,v)−Eh(σ e
j+1,v), (4.8)
(θj+1z ,w) + (θj+1p ,divw) =−(eI,j+1z ,w), (4.9)
(c0∂tθj+1p , q) + (div∂tθj+1u , q)−(divθzj+1, q) = (ω1j+1+ωj+12 , q), (4.10)
where
ω1j+1=c0( ¯∂tpj+1−p˙j+1−∂¯teI,j+1p ), ω2j+1= div( ¯∂tuj+1−u˙j+1−∂¯teI,j+1u ).
Theorem 4.1. Suppose that(u,z, p)is an exact solution of(2.5–2.7) with sufficient regularity and a fully discrete solution{(Uj,Zj, Pj)}1≤j≤N is defined by (4.2–4.4) with initial data satisfying (3.8). LetM1 be the maximum ofkukW2,1H1∩W1,1H2+ kpkW2,1L2∩W1,1H1 and kukL∞H2 +kpkW1,2H2∩L∞H1 and M2 be the maximum of kukW1,2H2∩W2,2H1+kpkW1,2H1∩W2,2L2 and kpkW1,1H2. Then
1≤i≤Nmax kui−Uik1,h+ max
1≤i≤Nkpi−Pik0≤c(∆t+h)M1,
1≤i≤Nmax kzi−Zik0≤c(∆t+h)M2, with constantsc independent ofT0.
