LOCKING-FREE FINITE ELEMENT METHODS FOR BIOT’S CONSOLIDATION MODEL IN POROELASTICITY

(1)

JEONGHUN LEE

Abstract. In this paper we propose a new locking-free finite element scheme for Biot’s consolidation model in poroelasticity. The main advantage of our method is that uniform-in-time bounds of the pressure and the flux errors are obtained with an analytic proof. We show a priori error estimates of semidiscrete solutions and of fully discrete solutions with the backward Euler scheme.

1. Introduction

The Biot’s consolidation model describes deformation of saturated elastic porous media and viscous fluid flow inside, simultaneously ^MR0066874[2]. Since this model naturally arises from classical engineering applications such as geomechanics and petrolium engineering, the study of numerical solutions of it has been of great interest.

There is an extensive literature on numerical schemes for Biot’s consolidation model with finite element methods. In researches using continuous Galerkin (CG) finite elements, nonphysical pressure oscillations of numerical solutions were ob-

served in certain ranges of material parameters and time step sizes MR603175,NAG:NAG1610080304,NAG:NAG1610080106

[18, 15, 22], which is called poroelasticity locking. In studies of Biot’s model in incompressible porous media using mixed elements for Stokes equations MR1156589,MR1257948,MR1393902

[8, 9, 10], Murad et.al.

proved decay of numerical pressure oscillation in time.

In order to overcome the locking problem, various other numerical methods were suggested, for example, a least square method[5], a stabilized CG method (Wan),^MR2177147 coupling mixed and CG mehtods MR2327964,MR2327966

[11, 12], discontinuous Galerkin (DG) methods (Liu), coupling mixed and DG methods [13]. Recently, coupling nonconforming^MR2461315 and mixed methods NUM:NUM21775

[19], and a stabilized DG method ^MR3047799[3] are proposed. For more information of previous studies we refer NUM:NUM21775,MR3047799, lewis1998finite

[19, 3, 6] and the references therein.

According to the heuristic analysis in phillips-wheeler

[14] the locking phenomena occurs when the constrained specific storage coefficient, denoted by c0 ≥ 0, is very close to 0. It is somewhat surprising that many of the aforementioned methods do not provide complete analytic proof that they are locking-free but only show numerical evidences. More precisely, whenc0is not uniformly positive, a priori error estimates of pressure is obtained only inL²([0, T0];L²) norm, not inL^∞([0, T0];L²) norm.

In this paper we aim to provide a new locking-free finite element scheme for Biot’s consolidation model. To the best of the author’s knowledge, our method gives the first analytic result on a priori error bound of pressure inL^∞([0, T0];L²) norm, without any time weights and an assumption of uniformly positive c0. As in NUM:NUM21775

[19] a coupling of nonconforming and mixed finite elements is considered but we use more sophisticated nonconforming elements, say Mardal–Tai–Winther type

1

(2)

elements. These are nonconforming vector-valued H¹ elements, originally devel-

oped for Brinkman (or Stokes–Darcy flow) problems MR1950614,MR2283095,MR2991835,MR2495068

[7, 17, 4, 21]. It turns out that they have nice features to give locking-free error numerical schemes for linear poroelasticity problems.

This paper is organized as follows. In section 2 we introduce notations and^sec:prelim^sec:prelim review preliminary knowledges to define our numerical scheme. We show a priori

error analysis of semidiscrete and fully discrete solutions in sectionssec:semidiscrete-analysissec:semidiscrete-analysis

3sec:fullydiscrete-analysissec:fullydiscrete-analysis

4, respectively.

Finally, we give remarks for extension of the method to higher order and rectangular elements in section sec:extensionsec:extension

5.

2. Preliminaries sec:prelim

2.1. Notations. Let Ω be a bounded polygonal domain inRⁿwithn= 2 or 3. For a nonnegative integerm,H^m(Ω),H^m(Ω;Rⁿ) denote the standardRandRⁿ-valued Sobolev spaces based onL² norm.

For a Banach space X and 0 < T₀ < ∞, C⁰([0, T₀];X) denotes the set of functions f : [0, T₀]→X which are continuous in t∈[0, T₀]. For an integerm≥1 we define

C^m([0, T0];X) ={f|∂ⁱf /∂tⁱ ∈C⁰([0, T0];X),0≤i≤m},

where∂f /∂t is the time derivative in the sense of the Fr´echet derivative inX (see e.g.,Yosida-book

[20]). For a functionf : [a, b]→X, we define the space-time norm kfk_Lp([a,b];X)=





 Rb

akfk^p_Xds1/p

, 1≤p <∞, esssup_t∈[a,b]kfkX, p=∞.

If the time interval is fixed as [0, T0], then we useL^pX to denoteL^p([0, T0];X) for simplicity. We define the space-time Sobolev spacesW^k,p([0, T0];X) for nonnegative integerkand 1≤p≤ ∞as the closure ofC^k([0, T0];X) with the normkuk_Wk,pX= Pk

i=0k∂ⁱu/∂tⁱkL^pX. We adopt a convention thatkf, gkX =kfkX+kgkX for the norm of a Banach space X. For simplicity of notations, ˙σ is used to denote the time derivative ofσ.

For a triangulation of Ω, Th is used to denote a shape-regular triangulation for which h is the maximum diameter of triangles (or tetrahedra) and Eh is the corresponding set of edges (faces), respectively. For E ∈ Eh and functionsf,g : Eh→Rⁿ we define

hf,giE= Z

E

f·gds, hf,gi= X

E∈E_h

hf,giE.

For an integerk≥0, a set G⊂Rⁿ, Pk(G) is the space of polynomials defined on Gof degree ≤k. We usePk(Th;X) to denote the space of piecewise polynomials adapted on Th of degree less ≤ k. For a vector space X, we use Pk(G;X) and Pk(Th;X) to denote the space ofX-valued polynomials with same conditions.

2.2. The Biot’s consolidation model. In this section we review the Biot’s consolidation model in poroelasticity. In an elastic porous medium, fluid flow and deformation of porous medium are intimately related when permeability of the medium is small, i.e., deformation of the medium and retardation of fluid flow occur simultaneously. Biot’s consolidation model describes this phenomena in saturated porous media.

(3)

Throughout this paper we restrict our interest on quasistatic consolidation problems. In other words, we assume that consolidation process is slow and acceleration term is ignored. In our description of the model,u is the displacement of porous media,pis fluid pressure,f is the body force,gis the source/sink density function of fluid. The governing equations of Biot’s consolidation model are

−divC(u) +α∇p=f, eq:strong-eq1

eq:strong-eq1 (2.1)

c₀p˙+αdiv ˙u−div(κ e

∇p) =g, eq:strong-eq2

eq:strong-eq2 (2.2)

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 generalC is a rank 4 tensor giving a symmetric positive definite linear map onL²(Ω;R^n×nsym) into itself, whereR^n×nsym is the space of symmetricn×nmatrices. In the rest of this paper, we assume that the elastic medium is homogeneous isotropic, so the elasticity tensorC has the form

Cτ e

= 2µτ e

+λtr(τ e

)I, τ

e

∈L²(Ω;R^n×nsym), eq:homog-elasticity

eq:homog-elasticity (2.3)

where the constants µ, λ >0 are Lam´e coefficients, I is the identity matrix. We assume that µ, λ are bounded above and below, in particular, λ does not tend to infinity. In (eq:strong-eq2eq:strong-eq2

2.2) the quantity c0p measures the amount of fluid that can be injected into a fixed material volume andc0≥0 is determined by the permeability of porous media and bulk moduli of solid and fluid. The hydraulic conductivity tensorκ

e

is defined by the permeability tensor of solid divided by fluid viscosity and it is positive definite. We assume that κ

e

is uniformly bounded above and below.

On details of deriving these equations from physical modelling, we refer standard porous media texts, e.g., anandarajah2010computational

[1].

In order to be a well-posed problem, the equations (eq:strong-eq1eq:strong-eq1

2.1–eq:strong-eq2eq:strong-eq2

2.2) need appropriate boundary and initial conditions. We assume that there are partitions of∂Ω which are

∂Ω = Γp∪Γf, ∂Ω = Γd∪Γt, |Γd|>0.

We also assume that boundary conditions are given to be p(t) = 0 on Γp, −κ

e

∇p(t)·n= 0 on Γf, u(t) = 0 on Γd, σ e

(t)n= 0 on Γt, eq:bc

eq:bc (2.4)

for all t ∈ [0, T₀] where n is the outward unit normal vector field on ∂Ω and σ

e

:=C(u)+αpI, the Cauchy stress tensor. Here we only consider the homogeneous boundary condition for simplicity but our method can be extended naturally to problems with nonhomogeneous boundary conditions. We also assume that given initial datap(0),u(0) andf(0) satisfy the compatibility condition (eq:strong-eq1eq:strong-eq1

2.1).

Well-posedness of this problem and regularity of solutions were thoroughly stud- ied in[16]. When we claim a priori error estimates we assume that exact solutions^MR1790411 are sufficiently regular to obtain the error bounds. In the rest of this paper the parameters α,κ

e

are assumed to be 1 and the identity for simplicity. In contrast, c0 ≥0 is not necessarily uniformly positive because we are mainly concerned with robust error bounds for generalc0≥0. In previous studies it is observed that the poroelasticity locking occurs when thec₀, the hydraulic conductivity, and time step are very small.

(4)

2.3. Variational formulation. Assumingα= 1,κ e

= identity, and introducing a new unknownz:=κ

e

∇pin (eq:strong-eq1eq:strong-eq1

2.1–eq:strong-eq2eq:strong-eq2

2.2), we have a new system

−divC(u) +∇p=f, eq:new-strong-eq1

eq:new-strong-eq1 (2.5)

z− ∇p= 0, eq:new-strong-eq2

eq:new-strong-eq2 (2.6)

c₀p˙+ div ˙u−divz=g.

eq:new-strong-eq3

eq:new-strong-eq3 (2.7) Let

Σ ={u∈H¹(Ω;Rⁿ)|u|Γ_d= 0}, V ={z∈H(div,Ω;Rⁿ)|z·n|Γ_t = 0}, W =L²(Ω).

A variational formulation of (eq:new-strong-eq1eq:new-strong-eq1

2.5-eq:new-strong-eq3eq:new-strong-eq3

2.7) is to seek (u,z, p)∈Σ×V ×W such that a(u,v)−(p,divv) = (f,v), v∈Σ,

eq:weak-eq1

eq:weak-eq1 (2.8)

(z,w) + (p,divw) = 0, w∈V, eq:weak-eq2

eq:weak-eq2 (2.9)

(c0p, q) + (div ˙˙ u, q)−(divz, q) = (g, q), q∈W, eq:weak-eq3

eq:weak-eq3 (2.10)

where, regarding (eq:homog-elasticityeq:homog-elasticity

2.3) and the integration by parts,

a(u,v) = (C(u), (v)) = 2µ((u), (v)) +λ(divu,divv), u,v∈Σ.

(2.11)

Assuming that|Γd|>0, the coercivity ofa(·,·) on Σ is obtained by Korn’s inequality.

2.4. finite element spaces. In our numerical scheme based on (eq:weak-eq1eq:weak-eq1

2.8–eq:weak-eq3eq:weak-eq3

2.10) we need three finite element spaces Σh,Vh,Wh for unknownsu, z, p, respectively.

We chooseVh,Wh as the lowest order Raviart–Thomas and discontinuous finite element spaces

V_h:={w∈H(div,Ω)|w|T ∈(P0(T)ⁿ+~xP0(T)), ∀T ∈ Th, w·n|Γf = 0}, W_h:={q∈L²(Ω)|q|T ∈ P0(T)}.

Let Π⁰_h be the canonical interpolation operator intoVh, andQh be the orthogonal L² projection intoWh. It is well-known that

divVh=Wh, div Π⁰_hw=Qhdivw, forw∈H¹(Ω;Rⁿ), eq:commute1

eq:commute1 (2.12)

kw−Π⁰_hwk0.hkwk1, kq−Qhqk0.hkqk1. eq:approx1

eq:approx1 (2.13)

Furthermore, the following inf-sup condition hold.

06=q∈Winf _h sup

w∈Vh

(q,divw) kqk0kwkdiv

≥C >0.

eq:inf-sup1

eq:inf-sup1 (2.14)

For Σ_h we use vector-valued nonconformingH¹ finite elements. Before we de- scribe necessary properties of Σ_h we need some notations. Let us define

|v|²_1,h= X

T∈Th

k∇vk²_0,T, kvk²_1,h:=|v|²_1,h+ X

E∈Eh

1

h_EkJvKk²_0,E,

for a piecewise polynomial functionv where JvKis the jump term on edges/faces.

LetRM(G) be the space of rigid body motions on a setG⊂Rⁿ. LetT₁, T₂ ∈ Th

be any two adjacent triangles/tetrahedra and M := T₁∪T₂. We suppose that

(5)

there exists a finite element space Σh, with an interpolation Πh:H¹(Ω;Rⁿ)→Σh, satisfying

Σh⊂H(div,Ω), div Σh=Wh, div Πhv=Qhdivv, eq:sigma-property1

eq:sigma-property1 (2.15)

there existsC >0 such that inf

06=q∈W_h sup

06=v∈Σh

(q,divv)

kqk₀kvk_1,h ≥C >0, eq:sigma-property2

ifv∈Σ_hsatisfiesv|T1∈RM(T₁),v|T2 ∈RM(T₂), thenv|M ∈RM(M), eq:sigma-property3

eq:sigma-property3 (2.17) and

kv−Πhvk0.hkvk1, kΠhvk1,h.kvk1, kΠhv−vk1,h.hkvk2. eq:sigma-property4

There are some known elements satisfying the above conditions. For example, two and three dimensional triangular elements are constructed in [7] and in^MR1950614 ^MR2283095[17], respectively. Rectangular elements are also constructed in [3]. In^MR3047799^MR2991835[4] some higher order triangular elements are developed.

From now on we assume an additional condition for Σ_h, which is that all DOFs associated to edges/faces in Γdare vanishing.

eq:sigma-property5

With this assumption, it is not difficult to see that k · k1,h is a norm on Σh due to the boundary condition on Γ_d and the jump terms. Furthermore, since JvK is mean-value zero on each interior edgeE∈ Eh from the definition of Σ_h, a standard scaling argument leads tok · k_1,h∼ | · |_1,hon Σ_h.

lemma:korn Lemma 2.1. The following inequality holds.

kvk1,h.kh(v)k0, v∈Σh, where_h is the element-wise symmetric gradient.

Proof. By the definition of ah, kh(v)k²₀ .ah(v,v) holds forv ∈Σh. Due to the equivalencek · k1,h∼ | · |1,hit is enough to prove

|v|1,h.kh(v)k0, v∈Σh.

Suppose that there is v∈Σh such that kh(v)k0= 0. Thenv is a piecewise rigid body motion on each triangle T ∈ Th. The interelement continuity of Σh implies thatJvKis vanishing on all edgesE in the interior of Ω. Thereforev is a rigid body motion on Ω, and we conclude that v = 0 due to the boundary condition on Γd. Since both of| · |1,hand kh(·)k0 are norms on a finite dimensional space Σh they are equivalent. Furthermore, these two norms have same scaling, so equivalence

constants do not depend onh. Proof is completed.

As a consequence, the following equivalence holds.

kvk_1,h∼ |v|_1,h∼ k_h(v)k₀, v∈Σ_h. 3. Error analysis of semidiscrete solutions sec:semidiscrete-analysis

In this section we discuss well-posedness of semidiscrete problems of (eq:new-strong-eq1eq:new-strong-eq1

2.5–eq:new-strong-eq3eq:new-strong-eq3

2.7) and the a priori error analysis. In the rest of this paper Σ_h×V_h×W_h will always denote the finite element spaces introduced in the previous section.

(6)

3.1. Semidiscrete problem. We consider a semidiscrete problem of (eq:weak-eq1eq:weak-eq1

2.10), which is seeking (uh,zh, ph) : [0, T0]→Σh×Vh×Whsuch that

ah(uh,v)−(ph,divv) = (f,v), v∈Σh, eq:weak-eq1-disc

eq:weak-eq1-disc (3.1)

(κ⁻¹zh,w) + (ph,divw) = 0, w∈Vh, eq:weak-eq2-disc

eq:weak-eq2-disc (3.2)

(c0p˙h, q) + (div ˙uh, q)−(divzh, q) = (g, q), q∈Wh, eq:weak-eq3-disc

eq:weak-eq3-disc (3.3) with

a_h(u,v) := 2µ(_h(u), _h(v)) +λ(divu,divv), u,v∈Σ_h. eq:ah

eq:ah (3.4)

Note that_h(u) =(u) ifu∈H¹(Ω;Rⁿ).

3.2. Well-posedness of semidiscrete problem. In this section we discuss existence and uniqueness of semidiscrete solutions. We first point out that the system (eq:weak-eq1-disceq:weak-eq1-disc

3.1–eq:weak-eq3-disceq:weak-eq3-disc

3.3) is not a system of ordinary differential equations (ODE) but a differential algebraic equation (DAE), so a theory of ODE is not directly applicable. In general proof of existence and uniqueness of DAE solutions, is not a consequence of a standard DAE theory if the system is large.

To explain a difficulty of applying standard DAE theory, here we discuss more on well-posedness of linear DAE problems. LetCbe the complex field andmbe a positive integer. ForE0, E1∈C^m×m,F∈C⁰([0,∞);C^m),X0∈C^m, a linear DAE is to seekX ∈C¹([0,∞),C^m) such that

E0X˙ +E1X=F, X(0) =X0.

Note thatE0 is not necessarily nonsingular. In a standard DAE theory, this DAE is called awell-posed DAEif the matrix pencil (E0, E1) is regular, i.e., there exists λ∈Csuch thatE0+λE1is nonsingular. However, this definition is easy to mislead because it does not guarantee existence of solutions for arbitrary initial data. In fact, it is known that a well-posed DAE has a unique solution only when initial data is compatible, i.e., the initial dataX0satisfies some intrinsic algebraic equations of the DAE. Furthermore, to recognize those intrinsic algebraic equations from a DAE, one needs to compute Jordan canonical forms of the augmented matrix [E₀E₁] in general, and it is certainly not practical for large systems. Therefore, instead of adopting a standard theory of DAEs, we clarify compatible initial data of (eq:weak-eq1-disceq:weak-eq1-disc

3.3), and prove existence and uniqueness of solutions directly.

Definition 3.1. A triple(v⁰,w⁰, q⁰)∈Σ_h×V_h×W_h is called a compatible initial data of (eq:weak-eq1-disceq:weak-eq1-disc

3.3) if

a_h(v⁰,v)−(q⁰,divv) = (f(0),v), v ∈Σ_h, eq:ic-comp1

eq:ic-comp1 (3.5)

(w⁰,w) + (q⁰,divw) = 0, w∈V_h, eq:ic-comp2

eq:ic-comp2 (3.6) hold.

The following lemma allows us to find a compatible initial data for numerical schemes.

lemma:ic-comp Lemma 3.2. For any givenq⁰ ∈Wh one can find a unique pair(v⁰,w⁰)∈Σh×Vh

such that (v⁰,w⁰, q⁰)is a compatible initial data of (eq:weak-eq1-disceq:weak-eq1-disc

3.3).

Proof. Suppose that q⁰ ∈ W_h is given. Since the bilinear form a_h is symmetric positive definite by its definition in (^eq:ah^eq:ah3.4) and Lemma2.1. Hence (^lemma:korn^lemma:korn eq:ic-comp1eq:ic-comp1

3.5) has a unique

(7)

solution v⁰ ∈ Σh. It is obvious that (eq:ic-comp2eq:ic-comp2

3.6) has a unique solution w⁰ ∈ Vh. By definition (v⁰,w⁰, q⁰) is a compatible initial data of (eq:weak-eq1-disceq:weak-eq1-disc

3.3).

Theorem 3.3. For given initial dataph(0)∈Wh andf ∈C¹([0, T0];L²(Ω;Rⁿ)), g∈C⁰([0, T0];L²(Ω)) there exists a unique solution(uh,zh, ph)∈C¹([0, T0]; Σh× Vh×Wh)of(eq:weak-eq1-disceq:weak-eq1-disc

3.3).

Proof. We will show that (eq:weak-eq1-disceq:weak-eq1-disc

3.3) is equivalent to an ODE system and check well- posedness of the ODE system. In order for that we rewrite (eq:weak-eq1-disceq:weak-eq1-disc

3.3) as a matrix equation using bases of finite element spaces.

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 write u_h = P

iα_iφ_i, z_h = P

iβ_iψ_i, p_h = P

iγ_iχ_i, P_hf = ζ_iψ_i,

Q_hg=ξ_iχ_i with (time-dependent) coefficient vectorsα,β,γ,ζ,ξ. Then (eq:weak-eq1-disceq:weak-eq1-disc

3.3) can be written as a matrix equation, which is





0 0 0

B^Tup 0 App









˙ α β˙

˙ γ



+





Auu 0 −Bup

0 Azz Bzp

0 −B^Tzp 0







 α β γ



=



 ζ 0 ξ



. eq:matrix-eq

eq:matrix-eq (3.7)

It is obvious thatAzzis symmetric positive definite andApp is symmetric positive semidefinite due toc0≥0. By Lemma2.1,^lemma:korn^lemma:kornAuu is also symmetric positive definite.

The first and second rows of (eq:matrix-eqeq:matrix-eq

3.7) give

α=A⁻¹uuBupγ+A⁻¹uuζ, β=−A⁻¹zzBzpγ.

eq:uz-replace

eq:uz-replace (3.8)

The third row of (eq:matrix-eqeq:matrix-eq

3.7) gives B^Tupα˙ +Appγ˙ −B^Tzpβ =ξ. Substituting α and β in this equation using (eq:uz-replaceeq:uz-replace

3.8) one obtains

(B^TupA⁻¹uuBup+App) ˙γ+B^TzpA⁻¹zzBzpγ=−B^TupA⁻¹uuζ˙ +ξ, which is an ODE system ofγ. The inf-sup condition (eq:sigma-property2eq:sigma-property2

2.16) implies that the matrix Bup is injective. SinceB^TupA⁻¹uuBup is positive definite andApp 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, the invert- ibility of B^TupA⁻¹uuBup+App is stable as App ≈0, so this well-posedness is robust for c0 ≈0. Note that α, β are uniquely determined by (eq:uz-replaceeq:uz-replace

3.8). Now the assertion follows by considering equivalence of (eq:matrix-eqeq:matrix-eq

3.7) and (eq:weak-eq1-disceq:weak-eq1-disc

3.3).

3.3. Error analysis. We state and prove main results of semidiscrete error analysis. We denote the discreteH¹space byH_h¹ with the normk · k1,h.

thm:semidiscrete-error Theorem 3.4. Suppose that u, p∈W^1,1([0, T];H²)and(u_h,z_h, p_h)is a solution of (eq:weak-eq1-disceq:weak-eq1-disc

3.3) with a compatible initial data such thatph(0) =Qhp(0). Then ku−uhk_L∞H¹_h+kp−phkL^∞L²≤ch(kukW^1,1H²+kpkW^1,1H¹),

holds withcindependent ofh. Furthermore, if we assume thatu, p∈W^2,1([0, T];H²), then

kz−z_hk_L∞L² ≤ch(kuk_W2,1H²+kpk_W2,1H¹).

For the proof of it we need some preliminary results.

(8)

lemma:consistency-err-bound Lemma 3.5. Forτ e

∈H¹(Ω;R^n×n) andv∈Σh define Eh(τ

e

,v) := X

E∈Eh

hτ e

nE,JvKiE.

Then the following inequality holds.

|Eh(τ e

,v)|.hkτ e

k1kvk1,h.

Proof. Let E ∈ E_h be an interior edge/face in Ω and T₊, T₋ be two distinct triangles/tetrahedra sharing E as the common boundary. The shape regularity of meshes yieldsh_T₋∼h_T₊∼h_E.

SinceJvKis perpendicular to constants inL²(E;Rⁿ), hτ

e

nE,JvKiE=hτ e

nE−c,JvKiE≤ inf

c∈Rⁿkτ e

nE−ck0,EkJvKk0,E.

By a standard scaling argument and shape regularity of meshes one can see that kJvKk0,E .h_E¹²k∇vk0,T₊∪T−,

c∈infRⁿ

kτ e

nE−ck0,E .h

1 2

Ek∇τk0,T₊∪T−.

For an edgeE on boundary and the triangleT ∈ ThcontainingEin the boundary, a similar argument gives

kJvKk_0,E .h_E¹²k∇vk_0,T,

c∈infRⁿ

kτ e

nE−ck0,E .h

1 2

Ek∇τk0,T.

By the Cauchy–Schwarz inequality and the results obtained above, we have E_h(τ

e

,v).hkτ e

k1|v|1,h,

as desired.

The following simple lemma will be useful in our error analysis.

lemma:basic-ineq Lemma 3.6. Suppose that A, B, C, D >0 satisfy A²+B²≤CA+D.

Then either A+B≤4C orA+B ≤2√

D holds.

Proof. Since eitherCA≤D orCA≥D one of the following inequalities holds.

A²+B²≤2D, A²+B²≤2CA.

eq:two-ineq

eq:two-ineq (3.9)

If the first inequality in (eq:two-ineqeq:two-ineq

3.9) holds, then

(A+B)²≤2(A²+B²)≤4D, which impliesA+B≤2√

D.

Suppose that the second inequality in (eq:two-ineqeq:two-ineq

3.9) holds. IfB ≥A, then dividing the inequality byAgives

A+B ≤(A²+B²)/A≤2C.

If B ≤A, then dividing the second inequality in (eq:two-ineqeq:two-ineq

3.9) by A gives A ≤ 2C, and

therefore we haveA+B≤2A≤4C as desired.

Now we have all preliminaries to prove Theorem thm:semidiscrete-errorthm:semidiscrete-error

3.4.

(9)

Proof of Theorem thm:semidiscrete-errorthm:semidiscrete-error

3.4. The proof consists of three parts: estimates ofku−uhk_L∞H_h¹, kz−z_hk_L2L², andkp−p_hk_L∞L².

Estimate of ku−u_hk_L∞H¹_h : We denote the error terms with eu=u−uh, ez=z−zh, ep=p−ph. From the left-hand side of (eq:new-strong-eq1eq:new-strong-eq1

2.5) we can obtain

−(divC(u) +∇p,v) = (C(u), h(v))−(p,divv) +Eh(σ e

,v), eq:consistency-form

eq:consistency-form (3.10)

forv∈Σ_h, by the integration by parts. The difference of this and (eq:weak-eq1-disceq:weak-eq1-disc

3.1) gives ah(eu,v)−(ep,divv) =−Eh(σ

e

,v), v∈Σh. eq:err-eq1

eq:err-eq1 (3.11)

In addition, the differences of (eq:weak-eq2eq:weak-eq2

2.10) and (eq:weak-eq2-disceq:weak-eq2-disc

3.3) yield

(ez,w) + (ep,divw) = 0, w∈Vh, eq:err-eq2

eq:err-eq2 (3.12)

(c₀e˙_p, q) + (div ˙e_u, q)−(dive_z, q) = 0, q∈W_h. eq:err-eq3

eq:err-eq3 (3.13)

In order for the error analysis we define splitting of error terms as in (eq:err-split1eq:err-split1

3.14–eq:err-split3eq:err-split3

3.16) and In order for the error analysis we split the errors into

eu=e^I_u+e^A_u := (u−Π˜hu) + ( ˜Πhu−uh), eq:err-split1

eq:err-split1 (3.14)

ez=e^I_z+e^A_z := (z−Π⁰_hz) + (Π⁰_hz−zh), eq:err-split2

e_p=e^I_p+e^A_p := (p−Q_hp) + (Q_hp−p_h).

eq:err-split3

The definitions ofeÎ_u,eÎ_z,eÎ_p and properties in (eq:commute1eq:commute1

2.12), yield div ˙eÎ_u⊥W_h, diveÎ_z⊥W_h, eÎ_p⊥div Σ_h. eq:err-cancel2

eq:err-cancel2 (3.17)

Rewriting (3.11–êq:err-eq1êq:err-eq13.13) using (êq:err-eq3êq:err-eq3 eq:err-split1eq:err-split1

3.14–eq:err-split3eq:err-split3

3.16) in consideration of the orthogonalities in (eq:err-cancel1eq:err-cancel1

??), (eq:err-cancel2eq:err-cancel2

3.17), we have

ah(eÂ_u,v)−(eÂ_p,divv) =−ah(eÎ_u,v)−Eh(σ e

,v), v∈Σh, eq:new-err-eq1

eq:new-err-eq1 (3.18)

(eÂ_z,w) + (eÂ_p,divw) =−(eÎ_z,w), w∈Vh, eq:new-err-eq2

(c0e˙Â_p, q) + (div ˙eÂ_u, q)−(diveÂ_z, q) =−(c0e˙Î_p, q) q∈Wh. eq:new-err-eq3

LetX(t), Y(t) be defined by

X²(t) =ah(eÂ_u(t), eÂ_u(t)) +keÂ_p(t)k²_c₀, Y²(t) = Z t

0

ke^A_zkds.

If we takev = ˙eÂ_u,w=eÂ_z,q=eÂ_p in (eq:new-err-eq1eq:new-err-eq1

3.18–eq:new-err-eq3eq:new-err-eq3

3.20), and add those equations, then 1

2 d

dt(X(t))²+keÂ_z(t)k²₀=−a_h(eÎ_u,e˙Â_u)−E_h(σ e

,e˙Â_u)−(eÎ_z, eÂ_z)−(c₀e˙Î_p, eÂ_p).

Integration of this from 0 tot in time, gives (X(t))²+ 2(Y(t))²

eq:interm-err-eq1

eq:interm-err-eq1 (3.21)

= (X(0))²+ 2 Z t

0

−ah(e^I_u,e˙^A_u)−Eh(σ e

,e˙Â_u)−(eÎ_z, eÂ_z)−(c0e˙Î_p, eÂ_p) ds,

=: (X(0))²+ Φ₁(t) + Φ₂(t) + Φ₃(t) + Φ₄(t).

Defining ¯tby

X(¯t) = sup

0≤t≤T0

X(t), eq:tmax-def

eq:tmax-def (3.22)

(10)

it suffices to estimateX(¯t). Moreover, since one of the followings is always true, Φ1(¯t) + Φ2(¯t) + Φ4(¯t)≤Φ3(¯t) or Φ1(¯t) + Φ2(¯t) + Φ4(¯t)≥Φ3(¯t), we always have

either (X(¯t))²+ 2(Y(¯t))²≤(X(0))²+ 2Φ₃(¯t), eq:estm-case1

eq:estm-case1 (3.23)

or (X(¯t))²+ 2(Y(¯t))²≤(X(0))²+ 2(Φ₁(¯t) + Φ₂(¯t) + Φ₄(¯t)).

eq:estm-case2

eq:estm-case2 (3.24)

Thus it is enough to prove estimates ofX(¯t) for these two cases separately.

As a preliminary result, we show an estimate ofX(0). Note first thateÂ_p(0) = 0, so X(0)²=keÂ_u(0)k²_a,h. To estimatekeÂ_u(0)k²_a,h, takev=eÂ_u(0) in (eq:new-err-eq1eq:new-err-eq1

3.18) att= 0, and have

keÂ_u(0)k²_a,h=−ah(eÎ_u(0), eÂ_u(0))−Eh(σ e

(0), e^A_u(0))

=−2µ(h(e^I_u(0)), _h(e^A_u(0)))−E_h(σ e

(0), e^A_u(0)), (∵dive^I_u(0)⊥W_h) By the Cauchy–Schwarz inequality, Lemma^lemma:korn^lemma:korn2.1 and Lemma lemma:consistency-err-boundlemma:consistency-err-bound

3.5, X(0) =ke^A_u(0)ka,h.h(ku(0))k2+kσ

e (0)k1).

eq:X0-estm

eq:X0-estm (3.25)

Now we show estimates ofX(¯t) for the cases (eq:estm-case1eq:estm-case1

3.23) and (eq:estm-case2eq:estm-case2

3.24).

Case I: Suppose that (eq:estm-case1eq:estm-case1

3.23) is true. Note that 2Φ2(¯t) =−4

Z ^¯t 0

(e^I_z, e^A_z)ds≤4 Z ^¯t

0

ke^I_zk²₀ds+ (Y(¯t))², eq:young-ineq1

eq:young-ineq1 (3.26)

by Young’s inequality and the definition ofY(¯t). Combining it with (eq:estm-case1eq:estm-case1

3.23) gives X²(¯t) +Y²(¯t)≤X²(0) + 4

Z ^¯t 0

ke^I_zk²₀ds eq:semi-xy-ineq1

eq:semi-xy-ineq1 (3.27)

≤X(0)²+ 4ke^I_zk²_L2L²

≤X(0)²+ch²kzk²_L2H¹. (Bramble–Hilbert lemma) Case II: Suppose that (eq:estm-case2eq:estm-case2

3.24) is true.

To estimate Φ1(¯t), we use the integration by parts in time, which gives Φ₁(¯t) =−2

Z ^¯t 0

a_h(eÎ_u,e˙Â_u)ds (∵diveÎ_u⊥W_h), eq:semi-Phi1-estm

eq:semi-Phi1-estm (3.28)

= 2ah(eÎ_u(¯t), eÂ_u(¯t))−2ah(eÎ_u(0), eÂ_u(0)) + 2 Z t^¯

0

ah( ˙e^I_u, e^A_u)ds

≤2 ke^I_u(¯t)ka,h+ke^I_u(0)ka,h+ Z t^¯

0

ke˙^I_u(0)ka,hds

! X(¯t)

≤chkuk_W1,1H²X(¯t).

(11)

Similarly,

Φ2(¯t) =−2 Z ^¯t

0

Eh(σ e

,e˙^A_u)ds eq:semi-Phi2-estm

= 2

−Eh(σ e

(¯t), e^A_u(¯t)) +Eh(σ e

(0), e^A_u(0)) + 2

Z ^¯t 0

Eh( ˙σ e

, e^A_u)ds

≤ch kσ e

(¯t)k1+kσ e

(0)k1+ Z t^¯

0

kσ˙ e

k1ds

! X(¯t)

≤chkσ e

k_W1,1H¹X(¯t).

For Φ4(¯t),

Φ₄(¯t) =−2 Z ^¯t

0

( ˙e^I_p, e^A_p)_c₀ds≤2 Z ^¯t

0

ke˙^I_pk_c₀dsX(¯t)≤chkpk_W1,1H¹X(¯t).

eq:semi-Phi4-estm

Combining (eq:estm-case2eq:estm-case2

3.24), (eq:semi-Phi1-estmeq:semi-Phi1-estm

3.28), (eq:semi-Phi2-estmeq:semi-Phi2-estm

3.29), and (eq:semi-Phi4-estmeq:semi-Phi4-estm

3.30), we have an inequality of the form in Lemma lemma:basic-ineqlemma:basic-ineq

3.6 withA=X(¯t),B =√

2Y(¯t),D=X(0), and C=ch(kσ

e

, pk_W1,1H¹+kuk_W1,1H²).

By Lemma (lemma:basic-ineqlemma:basic-ineq

3.6),

X(¯t)≤2X(0).h(ku(0)k2+kσ e

k1), eq:semi-xy-ineq2

eq:semi-xy-ineq2 (3.31)

or X(¯t)≤4C.h(kσ e

, pk_W1,1H¹+kuk_W1,1H²).

Note thatkpkL²H² ∼ kzkL²H¹ by the definition ofz. Then for any 0≤t≤T₀, by the triangle inequality, (eq:semi-xy-ineq1eq:semi-xy-ineq1

3.27), and (eq:semi-xy-ineq2eq:semi-xy-ineq2

3.31), ku(t)−u_h(t)k_H1

h ≤ ke^I_u(t)k_H1

h+ke^A_u(t)k_H1

eq:semi-u-estm h

eq:semi-u-estm (3.32)

.ke^I_u(t)k_H¹

h+X(¯t) (3.33)

.hmax{kukW^1,1H²+kσ e

, pkW^1,1H¹,kpkL²H²}, (3.34)

which completes the estimation ofku−uhk_L∞H_h¹.

Estimate of kp−p_hkL^∞L² : Note that the estimate ofX(t) is enough to give an estimate ofkp−p_hkL^∞L² whenc₀is uniformly greater than a positive constant.

We aim to show that anL^∞L² estimate ofp−p_h can be obtained and the error bound constants do not grow unboundedly asc₀tends to 0 andλtends to infinity.

By the triangle inequality it suffices to estimatee^A_p(t). By the inf-sup condition (eq:inf-sup1eq:inf-sup1

2.14), for any 06=q∈Wh, there existsv∈Σh such that

(divv, q⁰) = (q, q⁰), ∀q⁰∈Wh, kvk1,h≤ckqk0. If we use thisv to (eq:new-err-eq1eq:new-err-eq1

3.18) withq=e^A_p(t), then ke^A_p(t)k²₀=a_h(u(t)−u_h(t),v) +E_h(σ

e (t),v) .(ku(t)−uh(t)k1,h+hkσ

e

(t)k1)kvk1,h

.hmax{kuk_W1,1H²+kσ e

, pk_W1,1H¹,kpk_L2H²}ke^A_p(t)k₀. By the triangle inequality, the above estimate, and the Bramble–Hilbert lemma,

kp−phk_L∞L² .hmax{kuk_W1,1H²+kσ e

, pk_W1,1H¹,kpk_L2H²}.

(12)

Estimate ofkz−zhkL^∞L²: If we take time derivative of the equations (eq:new-err-eq1eq:new-err-eq1

3.18–eq:new-err-eq3eq:new-err-eq3

3.20) and repeat similar argument in the estimates ofku−uhk_L∞H_h¹ andkp−phkL^∞L², then we have

ku˙ −u˙hk_L∞H_h¹+kp˙−p˙hk_L∞L² .h(kuk_W2,1H²+kpk_W2,1H¹).

To estimate kz −z_hkL^∞L², by the triangle inequality, it is enough to estimate ke^A_zkL^∞L². To do it, takew =e^A_z in (eq:new-err-eq2eq:new-err-eq2

3.19), q=e^A_p in (eq:new-err-eq3eq:new-err-eq3

3.20), and add these two equations. Then we get

keÂ_zk²₀=−(eÎ_z, eÂ_z)−(c0e˙p, eÂ_p)−(div ˙eÂ_u, eÂ_p).

By the Cauchy–Schwarz inequality and the arithmetic-geometric mean inequality keÎ_zk0keÂ_zk0≤(keÎ_zk²₀+keÂ_zk²₀)/2, we have

keÂ_zk²₀.keÎ_zk²₀+ke˙_pk₀keÂ_pk₀+ke˙Â_uk_akeÂ_pk₀ .h²(kzk²_L∞H¹+

4. Error analysis of fully discrete solutions sec:fullydiscrete-analysis

In this section we consider error analysis of fully discrete solutions with the backward Euler scheme for time discretization. As in our error analysis for semidiscrete solutions we do not use Gronwall’s inequality, so the error bound does not contain exponentially growing factors. An additional technical difficulty arises from time discretization error terms which vanish in semidiscrete case because we want error bounds which are robust forc₀≥0. (kpk_c₀ does not boundkpk₀)

4.1. Well-posedness of fully-discrete problem. Let ∆t > 0 such that T₀ = N∆tfor an integerN, andt_j=j∆tforj= 0,1,· · ·, N. For a continuous function f defined on [0, T₀], we definef^j =f(t_j) andf^j+1/2=f(t_j+∆t/2). For a sequence {f^j}_j≥0, we define

∂¯_tf^j+1= f^j+1−f^j

∆t , time-quotient

time-quotient (4.1)

Note that forf defined on [0, T₀], ˆf^j+1/26=f^j+1/2in general.

Suppose that (U⁰,Z⁰, P⁰) is a compatible initial data defined in Lemma lemma:ic-complemma:ic-comp

3.2.

In the backward Euler scheme, (U^j+1,Z^j+1, P^j+1), the numerical solution at time t_j+1 is defined inductively by

ah(U^j+1,v)−(P^j+1,divv) = (f^j+1,v), eq:fulldisc1

eq:fulldisc1 (4.2)

(Z^j+1,w) + (P^j+1,divw) = 0, eq:fulldisc2

eq:fulldisc2 (4.3)

(c0∂¯tP^j+1, q) + (div( ¯∂tU^j+1), q)−(divZ^j+1, q) = (g^j+1, q), eq:fulldisc3

eq:fulldisc3 (4.4)

for (v,w, q)∈Σ_h×V_h×W_hand j≥0.

(13)

4.2. Well-definedness. First of all, we show that the fully discrete solution is well-defined. We need to check that (U^j+1,Z^j+1, P^j+1) is uniquely determined by the linear system (eq:fulldisc1eq:fulldisc1

4.2–eq:fulldisc3eq:fulldisc3

4.4) whenU^j, P^j,f^j+1, g^j+1 are given. Rewriting (eq:fulldisc1eq:fulldisc1

4.2–eq:fulldisc3eq:fulldisc3

4.4), a_h(U^j+1,v)−(P^j+1,divv) = (f^j+1,v),

(Z^j+1,w) + (P^j+1,divw) = 0,

(c₀P^j+1, q) + (divU^j+1, q)−∆t(divZ^j+1, q) = (c₀P^j, q) + (divU^j, q) + ∆t(g^j+1, q),

for (v,w, q)∈Σh×Vh×Wh. RegardingU^j+1,Z^j+1, P^j+1as unknowns, the above is a system of linear equations with same number of equations and unknowns. In order to show that it is well-defined, we assumeU^j =Z^j=P^j=f^j+1=g^j+1= 0 and will show thatU^j+1 =Z^j+1 =P^j+1 = 0. If we takev =U^j+1, w =Z^j+1, q=P^j+1 and add all equations, then we have

0≤a_h(U^j+1,U^j+1) + (Z^j+1,Z^j+1) + (c₀P^j+1, P^j+1) = 0, which impliesU^j+1=Z^j+1= 0. Note that (eq:fulldisc1eq:fulldisc1

4.2) is now given as (P^j+1,divv) = 0 for anyv ∈Σ_h. ThenP^j+1= 0 follows since div Σ_h=W_h. Hence the fully discrete solution is well-defined by induction.

Note that the backward Euler scheme enforces numerical solutions satisfy the compatibility conditions at every time step, so it does not need numerical initial data satisfying the compatibility conditions. However, initial data satisfying the compatibility conditions can be essential when one wants to use other numerical schemes, for instance, the Crank–Nicolson scheme. add reference of counterexam- ple.

4.3. Convergence. Let us denote the error (u^j−U^j,z^j−Z^j, p^j−P^j) by E_u^j =u^j−U^j = (u^j−Πhu^j) + (Πhu^j−u^j) =:e^I,j_u +θ^j_u, eq:err-decomp1

eq:err-decomp1 (4.5)

E_z^j =z^j−Z^j = (z^j−Π⁰_hz^j) + (Π⁰_hz^j−Z^j) =:e^I,j_z +θ_z^j, eq:err-decomp2

E_p^j=p^j−P^j= (p^j−Qhp^j) + (Qhp^j−P^j) =:e^I,j_p +θ^j_p. eq:err-decomp3

In the semidiscrete error analysis of the previous section, we already obtained error bounds ofke^I,j_u ka andke^I,j_p k0. By the triangle inequality, we only need to consider the a priori estimates ofkθ_u^jka andkθ_p^jk0.

Lemma 4.1. Forθ^j_u,θ^j_p defined in (eq:err-decomp1eq:err-decomp1

4.5)and (eq:err-decomp3eq:err-decomp3

4.7), kθ_p^jk0.kθ^j_uka+ke^I,j_u ka+hkσ

e

jk1. eq:thetap-estm1

eq:thetap-estm1 (4.8)

Proof. From (eq:consistency-formeq:consistency-form

3.10) att=t_j we have ah(u^j,v)−(p^j,divv) +Eh(σ

e

j,v) = (f^j,v), v∈Σh. The difference between this and the equation (eq:fulldisc1eq:fulldisc1

4.2) at the (j−1)-st time step, gives a_h(E_u^j,v)−(E_p^j,divv) +E_h(σ

e

j,v) = 0.

Rewriting it using (eq:err-decomp1eq:err-decomp1

4.5–eq:err-decomp3eq:err-decomp3

4.7) and the orthogonalitye^I,j_p ⊥divv forv∈Vh, we have (θ^j_p,divv) =ah(θ^j_u+e^I,j_u ,v) +Eh(σ

e

j,v).

eq:thetap-estm2

eq:thetap-estm2 (4.9)