JuhoK¨onn¨oRolfStenberg ANALYSISOFH(DIV)-CONFORMINGFINITEELEMENTSFORTHEBRINKMANPROBLEM HelsinkiUniversityofTechnologyInstituteofMathematicsResearchReports

Helsinki University of Technology Institute of Mathematics Research Reports

Espoo 2010 A582

ANALYSIS OF H(DIV)-CONFORMING FINITE ELEMENTS FOR THE BRINKMAN PROBLEM

Juho K ¨onn ¨o Rolf Stenberg

AB

HELSINKI UNIVERSITY OF TECHNOLOGY

Helsinki University of Technology Institute of Mathematics Research Reports

Espoo 2010 A582

ANALYSIS OF H(DIV)-CONFORMING FINITE ELEMENTS FOR THE BRINKMAN PROBLEM

Juho K ¨onn ¨o Rolf Stenberg

Helsinki University of Technology

Juho K¨onn¨o, Rolf Stenberg: Analysis of H(div)-conforming finite elements for the Brinkman problem

; Helsinki University of Technology Institute of Mathematics Research Reports A582 (2010).

Abstract: The Brinkman equations describe the flow of a viscous fluid in a porous matrix. Mathematically the Brinkman model is a parameter- dependent combination of the Darcy and Stokes models. We introduce a dual mixed framework for the problem, and use H(div)-conforming finite elements together with Nitsche’s method to obtain a stable formulation. We show the formulation to be stable in a mesh-dependent norm for all values of the pa- rameter. We also introduce a postprocessing scheme for the pressure along with a residual-based a posteriori estimator, which is shown to be efficient and reliable for all parameter values.

AMS subject classifications:

65N30

Keywords:

Brinkman equations, Nitsche’s method, mixed finite element meth- ods, a posteriori

Correspondence

Juho K¨onn¨o, Rolf Stenberg Aalto University

Department of Mathematics and Systems Analysis P.O. Box 11100

00076 Aalto Finland

firstname.lastname@tkk.fi

Received 2010-01-02

ISBN 978-952-248-282-2 (print) ISSN 0784-3143 (print) ISBN 978-952-248-283-9 (PDF) ISSN 1797-5867 (PDF) Helsinki University of Technology

Faculty of Information and Natural Sciences Department of Mathematics and Systems Analysis P.O. Box 1100, FI-02015 TKK, Finland

email: math@tkk.fi

Analysis of H (div)-conforming nite elements for the Brinkman problem

Juho Könnö

Rolf Stenberg

December 18, 2009

Abstract

The Brinkman equations describe the ow of a viscous uid in a porous matrix. Mathematically the Brinkman model is a parameter-dependent combination of the Darcy and Stokes models. We introduce a dual mixed framework for the problem, and use H(div)-conforming nite elements together with Nitsche's method to obtain a stable formulation. We show the formulation to be stable in a mesh-dependent norm for all values of the parameter. We also introduce a postprocessing scheme for the pressure along with a residual-based a posteriori estimator, which is shown to be ecient and reliable for all parameter values.

1 Introduction

In soil mechanics, the Brinkman equation describes the ow of a viscous uid in a very porous medium. For a derivation of and details on the equations we refer to [17, 1, 2, 3, 21]. As opposed to the Darcy model widely used in soil mechan- ics, the Brinkman model adds an eective viscosity to the equations. Typical applications of the equations lie in oil exploration, groundwater modelling, and some special applications, such as heat pipes [15]. The Brinkman model is also often used as a coupling layer between a free surface ow and a porous Darcy ow [10]. Mathematically, the Brinkman equations are a parameter-dependent combination of the Darcy and Stokes equations.

We study the application of H(div)-conforming nite elements designed for the Darcy problem to the more complicated Brinkman problem. H(div)- elements have been considered for the closely related Stokes problem in [9, 26, 14]. Our model constitutes an approximation with non-conforming basis func- tions, since in the discretizations of theH(div)-space only the normal component of the velocity is continuous on interelement boundaries. To enforce the tangen- tial continuity, we use the so-called Nitsche's method introduced in [20]. This in turn requires the use of a mesh-dependent bilinear form. The method has a strong resemblance to totally discontinuous Galerkin methods for the Stokes

∗Department of Mathematics and Systems Analysis, Helsinki University of Technology, P.O. Box 1100, FIN-02015 TKK, Espoo, Finland (jkonno@math.tkk.fi). Supported by the Finnish Cultural Foundation.

†Department of Mathematics and Systems Analysis, Helsinki University of Technology, P.O. Box 1100, FIN-02015 TKK, Espoo, Finland (rolf.stenberg@tkk.fi).

equation, cf. [12]. The motivation for using this non-conforming approximation is the fact thatH(div)-conforming elements are widely used in industry for solv- ing the Darcy equations, and we want to derive a way of easily incorporating viscosity to the existing solvers, along with a rigorous mathematical analysis.

In Part II of this paper we will report on numerical benchmark studies with our methods.

Methods based on Stokes elements are studied in e.g. [13, 5].

The structure of the article is as follows. First, we formulate the Brinkman problem mathematically and introduce the corresponding variational formula- tion. Then we introduce the nite element spaces used for the discretization along with the mesh-dependent energy norms in which the error is measured.

We derive a priori convergence results for this discretization, which implies the possibility of improving the accuracy of the pressure approximation via post- processing. This is addressed in a separate section, where we provide optimal a priori convergence results for the postprocessed solution. Finally, we end the paper with the a posteriori error analysis. A residual-based a posteriori error estimator is introduced. It is then shown that the estimator is both reliable and ecient for all values of the viscosity parameter.

We use standard notation throughout the paper. We denote by C, C₁, C₂ etc. generic constants, that are not always identical in value, but are always independent of the parametertand the mesh size h.

2 The Brinkman model

LetΩ⊂Rⁿ, withn= 2,3, be a domain with a polygonal or polyhedral bound- ary. We denote by u the velocity eld of the uid and by p the pore pres- sure. The equations are scaled as presented in [11], with the single parameter t representing the eective viscosity of the uid, which is assumed constant for simplicity. With this notation, the Brinkman equations are

−t²∆u+u− ∇p=f, in Ω (2.1)

divu=g, in Ω (2.2)

For t > 0, the equations are formally a Stokes problem. The solution (u, p) is sought in V ×Q = [H₀¹(Ω)]ⁿ ×L²₀(Ω). For the case t = 0 we get the Darcy problem, and accordingly the solution space can be chosen as V ×Q =H(div,Ω)×L²₀(Ω). For simplicity of the mathematical analysis, we consider for the case t > 0 homogenous Dirichlet boundary conditions for the velocity eld. Thus the boundary conditions are

u=0. (2.3)

For the limiting Darcy caset= 0we assume Neumann conditions

u·n= 0. (2.4)

In the following, we denote by(·,·)_D the standard L² inner product over a set D ⊂Rⁿ. If D = Ω, the subscript is dropped for convenience. Similarly,

h·,· i_B is theL² inner product over a set B ⊂Rⁿ⁻¹. We dene the following bilinear forms

a(u,v) =t²(∇u,∇v) + (u,v), (2.5)

b(v, p) = (divv, p), (2.6)

and B(u, p;v, q) =a(u,v) +b(v, p) +b(u, q). (2.7) The weak formulation of the Brinkman problem then reads: Find (u, p)∈ V ×Qsuch that

B(u, p;v, q) = (f,v) + (g, q), ∀(v, q)∈V ×Q. (2.8)

3 Solution by mixed nite elements

LetKhbe a shape-regular partition ofΩinto simplices. As usual, the diameter of an element K is denoted by hK, and the global mesh size h is dened as h= maxK∈KhhK. We denote byEh the set of all faces ofKh. We writehE for the diameter of a faceE.

We introduce the jump and average of a piecewise smooth scalar functionf as follows. Let E = ∂K∩∂K⁰ be an interior face shared by two elements K andK⁰. Then the jump off overE is dened by

[[f]] =f|K−f|K⁰. (3.1) and the average as

{f}= 1

2(f|K+f|K⁰). (3.2) For vector valued functions, we dene the jumps and averages analogously.

3.1 The mixed method and the norms

Mixed nite element discretization of the problem is based on nite element spaces V_h×Q_h ⊂ H(div,Ω)×L²₀(Ω) of piecewise polynomial functions with respect to K_h. We will focus here on the Raviart-Thomas (RT) and Brezzi- Douglas-Marini (BDM) families of elements [8]. In three dimensions the coun- terparts are the Nédélec elements [19] and the BDDF elements [7].That is, for an approximation of orderk≥1, the ux spaceV_h is taken as one of the following two spaces

V_h^RT ={v∈H(div,Ω)|v|K ∈[Pk−1(K)]ⁿ⊕xP˜k−1(K)∀K∈ Kh}, (3.3) V_h^BDM ={v∈H(div,Ω)|v|_K ∈[P_k(K)]ⁿ ∀K∈ K_h}, (3.4) where P˜k−1(K) denotes the homogeneous polynomials of degree k−1. The pressure is approximated in the space

Q_h={q∈L²₀(Ω)|q|_K ∈P_k−1(K)∀K∈ K_h}. (3.5) Notice that V_h^RT ⊂ V_h^BDM and Q^BDM_h = Q^RT_h . The combination of spaces satises the following equilibrium property:

div V_h⊂Q_h. (3.6)

To assure the stability of the non-conforming approximation, we use Nitsche's method [13, 20] with a suitably chosen stabilization parameterα. We dene the following mesh-dependent bilinear form

B_h(u, p;v, q) =a_h(u,v) +b(v, p) +b(u, q), (3.7) in which

ah(u,v) = (u,v) +t² X

K∈Kh

(∇u,∇v)K (3.8)

+t² X

E∈Eh

{ α

hEh[[u]],[[v]]i_E− h{∂u

∂n},[[v]]i_E− h{∂v

∂n},[[u]]i_E}.

Then the discrete problem is to nduh∈Vh andph∈Qhsuch that

Bh(uh, ph;v, q) = (f,v) + (g, q), ∀(v, q)∈Vh×Qh. (3.9) We introduce the following mesh-dependent norms for the problem. For the velocity we use

kuk²_t,h=kuk²+t² X

K∈Kh

k∇uk²_0,K+ X

E∈Eh

1

h_Ek[[u·τ]]k²_0,E

!

, (3.10)

and for the pressure

|||p|||²_t,h= X

K∈Kh

h²_K

h²_K+t²k∇pk²_0,K+ X

E∈Eh

h_E

h²_E+t²k[[p]]k²_0,E. (3.11) Note that both of the norms are also parameter dependent.

3.2 A priori analysis

First we prove the consistency of the modied method. For the exact solution (u, p) it holds [[u]]|_E = 0. Inserting u into the modied part of the bilinear form, we have using element-by-element partial integration

ah(u,v) = (u,v) +t² X

K∈Kh

(∇u,∇v)K− X

E∈Eh

h{∂u

∂n},[[v]]iE

!

= (u,v) +t² X

K∈Kh

{(∇u,∇v)_K− h∂u

∂n,vi_∂K}

= (u,v) +t² X

K∈Kh

(−∆u,v)K

= (−t²∆u+u,v).

This gives us the following result.

Theorem 3.1. The exact solution (u, p)∈V ×Q satises

B_h(u, p;v, q) = (f,v) + (g, q), ∀(v, q)∈V_h×Q_h. (3.12)

Next we prove the stability ofa_h(·,·)in the mesh-dependent norm (3.10).

The stability only holds in the discrete space V_h. First we recall, as shown in [24], that the normal derivative can be estimated as

h_Ek∂v

∂nk²_0,E≤C_Ik∇vk²_0,K, ∀v∈V_h. (3.13) Lemma 3.2. The bilinear form a_h(·,·) is coercive in the discrete space V_h; there exists a positive constant C such that

a_h(v,v)≥Ckvk²_t,h, ∀v∈V_h. (3.14) Proof. First we note that by Young's inequality

−2 X

E∈Eh

h{∂v

∂n},[[v]]iE=− X

K∈Kh

h∂v

∂n,[[v]]i∂K

≥ − X

K∈Kh

k∂v

∂nk0,∂Kk[[v]]k0,∂K

≥ −( X

K∈Kh

h_Kk∂v

∂nk²_0,∂K)^1/2( X

K∈Kh

h⁻¹_K k[[v]]k²_0,∂K)^1/2

≥ −1 2

X

K∈Kh

hKk∂v

∂nk²_0,∂K− 2

X

K∈Kh

h⁻¹_K k[[v]]k²_0,∂K

≥ −CI

2 X

K∈Kh

k∇vk²_0,K− 2

X

K∈Kh

h⁻¹_K k[[v]]k²_0,∂K,

for an arbitraryε >0. This immediately gives ah(v,v) =kvk²₀+t² X

K∈Kh

k∇vk²_0,E+t² X

E∈Eh

( α

hKk[[v]]k²_0,E−2h∂v

∂n,[[v]]iE)

≥min{1−CI

2, α−

2}kvk²_t,h. (3.15)

Here C_I is the constant from the discrete trace inequality (3.13). Since and α are free parameters, choosing > CI/2 and α > /2, we have the desired result.

Next, we prove the discrete Brezzi-Babuska stability condition. Recall that we only have to prove the condition in the Raviart-Thomas case since V_h^RT ⊂ V_h^BDM.

Lemma 3.3. There exists a positive constant C such that

v∈Vsuph

b(v, q)

kvkt,h ≥C|||q|||t,h, ∀q∈Qh. (3.16) Proof. We recall that the local degrees of freedom for the RT family are

hv·n, zi_E, ∀z∈P_k−1(E), (3.17) (v,z)K, ∀z∈[Pk−2(K)]ⁿ. (3.18)

Thus, for a given q∈Q_h we denev by hv·n, zi_E= h_K

h²_K+t²h[[q]], zi_E, ∀z∈P_k−1(E), (3.19) (v,z)_K =− h²_K

h²_K+t²(∇q,z)_K ∀z∈[P_k−2(K)]ⁿ. Choosingz= [[q]]∈P_k−1(E)andz=∇q∈[P_k−2(K)]ⁿ gives

hv·n,[[q]]i_E= hK

h²_K+t²k[[q]]k²_0,E, (3.20) (v,∇q)K=− h²_K

h²_K+t²k∇qk²_0,K. (3.21) An explicit inspection of the degrees of freedom yields the relation

h²_K+t²

h²_K kvk²_0,K≤ hK

h²_K+t²k[[q]]k²_0,E+ h²_K

h²_K+t²k∇qk²_0,K. (3.22) Thus, using scaling arguments, we have

kvk²_t,h≤C X

K∈Kh

h²_K+t²

h²_K kvk²_0,K≤C|||q|||²_t,h. (3.23) Next we use element-by-element partial integration on b(v, q), and apply the denitions (3.19) to get

b(v, q) = X

K∈Kh

−(v,∇q)_K+ X

E∈Eh

hv·n, qi_E (3.24)

= X

K∈Kh

h²_K

h²_K+t²(∇q,∇q)_K+ X

E∈Eh

h_K

h²_K+t²h[[q]],[[q]]i_E (3.25)

=|||q|||²_t,h. (3.26)

Combining (3.23) and (3.24) gives (3.16).

By the above stability results forah(·,·)and b(·,·) the following stability result holds, see e.g. [8].

Lemma 3.4. For some positive constant C it holds

(v,q)∈Vsuph×Qh

B_h(r, s;v, q)

kvkt,h+|||q|||t,h ≥C(krk_t,h+|||s|||_t,h), ∀(r, s)∈V_h×Q_h. (3.27) For interpolation inH(div), a special interpolation operator is required. We use the interpolation operator R_h : H(div,Ω) → V_h introduced by Schöberl in [22] satisfying

(div (v−R_hv), q) = 0, ∀q∈Q_h. (3.28) The interpolant satises the following properties. We denote by P_h:L²(Ω)→ Q_h theL²-projection. The equilibrium property (3.6) implies

(div v, q−P_hq) = 0, ∀v ∈V_h. (3.29)

Furthermore, we have the commuting diagram property:

divRh=Phdiv. (3.30)

Traditionally the interpolation operator has been dened based on the moments of the normal component of the velocity on the boundaries. By using the alter- native interpolation operator of [22] we need not the extra regularity assumption of the edge-based interpolation operators. In the following, the standard Sobolev norm of orderkis denotedk· k_k. The main result of the chapter is the following quasioptimal a priori result:

Theorem 3.5. There is a positive constant C such that

ku−uhkt,h+|||Php−ph|||t,h≤Cku−Rhukt,h. (3.31) Proof. By Lemma 3.4 there exists functions (v, q) ∈ Vh ×Qh such that kvkt,h+|||q|||t,h≤C, and

ku_h−R_huk_t,h+|||p_h−P_hp|||_t,h≤ B_h(u_h−R_hu, p_h−P_hp;v, q) (3.32)

=a_h(u_h−R_hu,v) + (div v, p_h−P_hp) + (div (u_h−R_hu), q)

=ah(u−Rhu,v) + (divv, p−Php) + (div (u−Rhu), q),

where the last line follows from the consistency of the method given by The- orem 3.1. By using the interpolation properties (3.28) and (3.29), we arrive at

kuh−Rhukt,h+|||ph−Php|||t,h≤ah(u−Rhu,v)≤Cku−Rhukt,h. (3.33) Using the triangle inequality yields the result of the theorem.

Analogously to the dual mixed formulation of the Poisson problem discussed e.g. in [4, 23, 18], we have a superconvergence result for |||p_h−P_hp|||_t,h. This implies that the pressure solution can be improved by local postprocessing.

Assuming full regularity, we conclude the chapter with the following a priori result:

ku−u_hk_t,h+|||P_hp−p_h|||_t,h≤

(C(h^k+th^k−1)kukk, for RT,

C(h^k+1+th^k)kuk_k+1, for BDM. (3.34)

4 Postprocessing method

In this section we present a postprocessing method for the pressure in the spirit of [18]. We seek the postprocessed pressure in an augmented spaceQ^∗_h ⊃Q_h, dened as

Q^∗_h=

({q∈L²₀(Ω) |q|K ∈Pk(K)∀K∈ Kh}, for RT,

{q∈L²₀(Ω) |q|_K ∈P_k+1(K)∀K∈ K_h}, for BDM. (4.1) The postprosessing method is: ndp^∗_h∈Q^∗_h such that

Php^∗_h=ph (4.2)

(∇p^∗_h,∇q)K= (−t²∆uh+uh−f,∇q)K, ∀q∈(I−Ph)Q^∗_h|K. (4.3)

The method can be compactly treated as an integral part of the problem by embedding it into the bilinear form. We introduce the modied bilinear form B^∗_h(u, p^∗;v, q^∗) =Bh(u, p^∗;v, q^∗)+ X

K∈Kh

h²_K

h²_K+t²(−∇p^∗+u−t²∆u,∇(I−Ph)q^∗)K. (4.4) The postprocessed problem is then: nd(uh, p^∗_h)∈Vh×Q^∗_hsuch that for every pair(v, q^∗)∈V_h×Q^∗_hit holds

B_h^∗(u_h, p^∗_h;v, q^∗) =L_h(f, P_hg;v, q^∗), (4.5) in which

Lh(f, g;v, q^∗) = (f,v) + (g, q^∗) + X

K∈Kh

h²_K

h²_K+t²(f,∇(I−Ph)q^∗)K. (4.6) We have the following theorem relating the solution of the postprocessed prob- lem to the original problem.

Theorem 4.1. Let (u_h, p^∗_h) ∈ V_h×Q^∗_h be the solution of the problem (4.5) and set p_h = P_hp^∗_h. Then (u_h, p_h) ∈ V_h×Q_h is the solution of the original problem (3.9). Conversely, if (u_h, p_h)∈V_h×Q_h is the solution of the original problem (3.9) andp^∗_h is dened as above, then(u_h, p^∗_h)∈V_h×Q^∗_his the solution to (4.5).

Proof. Testing with(v,0)∈V_h×Q^∗_h and using the equilibrium property (3.6) yields

B_h^∗(uh, p^∗_h;v,0) =ah(uh,v) + (divv, p^∗_h)

=a_h(u_h,v) + (divv, P_hp^∗_h)

=ah(uh,v) + (divv, ph) = (f,v).

On the other hand, testing with(0, Phq^∗)∈Vh×Qh⊂Vh×Q^∗_hgives B_h^∗(u_h, p^∗_h;0, P_hq^∗) = (divu_h, P_hq^∗) = (g, P_hq^∗).

Combining the above two equations yields the original problem (3.9) and rst part of the assertion is proved. Next take (u_h, p_h)to be the solution of (3.9), and p^∗_h the postprocessed pressure dened above. Using the denition of the

postprocessed pressure and the equilibrium property, we have B^∗_h(u_h, p^∗_h;v, q^∗) =B_h^∗(u_h, p^∗_h;v, P_hq^∗) +B_h^∗(u_h, p^∗_h;0,(I−P_h)q^∗)

=ah(uh,v) + (divv, p^∗_h) + (divuh, Phq^∗)

+ X

K∈Kh

h²_K

h²_K+t²(−t²∆u_h+u_h− ∇p^∗_h,∇(I−P_h)P_hq^∗)_K + (divuh,(I−Ph)q^∗) + X

K∈Kh

h²_K

h²_K+t²(f,∇(I−Ph)²q^∗)K

+ h²_K

h²_K+t²(−t²∆uh+uh− ∇p^∗_h−f,∇(I−Ph)²q^∗)K

=ah(uh,v) + (divv, Php^∗_h) + (div uh, Phq^∗) + (divuh, Ph(I−Ph)q^∗) + X

K∈Kh

h²_K

h²_K+t²(f,∇(I−Ph)²q^∗)K

=L_h(f, P_hg;v, q^∗)

for arbritrary (v, q^∗) ∈ Vh ×Q^∗_h. Thus the second part of the assertion is valid.

Next we show that the postprocessed method is stable in the discrete spaces.

For this we need the following lemma, essentially proved in [18].

Lemma 4.2. There exists positive constantsC1, C2such that for everyq^∗∈Q^∗_h it holds

|||q^∗|||_t,h≤ |||P_hq^∗|||_t,h+|||(I−P_h)q^∗|||_t,h≤C₂|||q^∗|||_t,h, (4.7) C₁|||q^∗|||_t,h≤ |||P_hq^∗|||_t,h+ X

K∈Kh

k∇(I−P_h)q^∗k²_0,K

!_1/2

≤C₂|||q^∗|||_t,h. (4.8) Since (I−P_h)q^∗ is L²-orthogonal to piecewise constant functions, we further- more have the following estimate, with C3>0,

|||(I−Ph)q^∗|||t,h≤C3

X

K∈Kh

k∇(I−Ph)q^∗k²_0,K

!_1/2

. (4.9)

We are now ready to prove the main stability result.

Theorem 4.3. There exists C >0 such that for every (u, p^∗) ∈ V_h×Q^∗_h it holds

(v,q^∗)∈Vsuph×Q^∗_h

B_h^∗(u, p^∗;v, q^∗)

kvk_t,h+|||q^∗|||_t,h ≥C(kukt,h+|||p^∗|||t,h). (4.10) Proof. Let(u, p^∗)∈Vh×Q^∗_h be arbitrary. Choosingq^∗=q∈Qhwe have

B_h^∗(u, p^∗;v, q) =ah(u,v) + (divv, p^∗) + (divu, q)

=a_h(u,v) + (divv, P_hp^∗) + (div u, q)

=Bh(u, Php^∗;v, q).

Thus Lemma 3.4 guarantees that there exists(v, q)∈V_h×Q_h such that kvk_t,h+|||q|||_t,h≤C(kuk_t,h+|||P_hp^∗|||_t,h)and

B_h(u, p^∗;v, q)≥C(kuk_t,h+|||P_hp^∗|||_t,h). (4.11) Next, we choose(v, q^∗) = (0,(I−P_h)p^∗)∈V_h×Q^∗_h.

B_h^∗(u, p^∗;0,(I−Ph)p^∗)

= (divu,(I−Ph)p^∗) + X

K∈Kh

h²_K

h²_K+t²(−∇p^∗+u−t²∆u,∇(I−Ph)p^∗)K

≥ −C(kukt,h|||(I−Ph)p^∗|||t,h) + X

K∈Kh

h²_K

h²_K+t²k∇(I−Ph)p^∗k²_0,K

+ X

K∈Kh

h²_K

h²_K+t²(−∇P_hp^∗+u−t²∆u,∇(I−P_h)p^∗)_K

≥ −C(kuk_t,h+|||P_hp^∗|||_t,h)|||(I−P_h)p^∗|||_t,h+ X

K∈Kh

h²_K

h²_K+t²k∇(I−P_h)p^∗k²_0,K

− X

K∈Kh

h²_K

h²_K+t²kt²∆uk²_0,K

!_1/2

|||(I−P_h)p^∗|||_t,h.

We can estimate the term containing the Laplacian using the inverse inequality as follows:

X

K∈Kh

h²_K

h²_K+t²kt²∆uk²_0,K≤ X

K∈Kh

t²h²_K h²_K+t²

1

h²_Kkt∇uk²_0,K≤t² X

K∈Kh

k∇uk²_0,K.

Using Young's inequality and the norm equivalence (4.9), we have for any >0 B_h^∗(u, p^∗;0,(I−Ph)p^∗)≥ −C1

2(kukt,h+|||Php^∗|||t,h)²−

2|||(I−Ph)p^∗|||²_t,h

+ X

K∈Kh

h²_K

h²_K+t²k∇(I−P_h)p^∗k²_0,K

≥ −C₁

2(kukt,h+|||Php^∗|||t,h)² + 1−C˜2

2

! X

K∈Kh

h²_K

h²_K+t²k∇(I−Ph)p^∗k²_0,K. Choosing= 1/C1 yields, withC3=C1C2/2,

B_h^∗(u, p^∗;0,(I−P_h)p^∗)≥ −C₃(kuk_t,h+|||P_hp^∗|||_t,h)² (4.12) +1

2 X

K∈Kh

h²_K

h²_K+t²k∇(I−P_h)p^∗k²_0,K.

Combining estimates (4.11) and (4.12), the norm equivalence (4.8), and choosing

δsuciently small, we have

B_h^∗(u, p^∗;v, q+δ(I−Ph)p^∗)≥(1−δC3)(kukt,h+|||Php^∗|||t,h)² δC3

X

K∈Kh

k∇(I−Ph)p^∗k²_0,K

≥C(kukt,h+|||p^∗|||t,h)². Furthermore,

kvkt,h+|||q+δ(I−Ph)p^∗|||t,h≤ kvkt,h+|||q|||t,h+δ|||(I−Ph)p^∗|||t,h

≤C(kuk_t,h+|||P_hp^∗|||_t,h) +δ|||(I−P_h)p^∗|||_t,h

≤C(kukt,h+|||p^∗|||t,h), yielding the desired result.

We have the following a priori result, which shows that given sucient reg- ularity, the postprocessed displacement converges with an optimal rate.

Theorem 4.4. For the postprocessed solution(u_h, p^∗_h)it holds ku−uhkt,h+|||p−p^∗_h|||t,h≤C inf

q^∗∈Q^∗_h

nku−Rhukt,h+|||p−q^∗|||t,h (4.13)

+ ( X

K∈Kh

h²_K

h²_K+t²k − ∇q^∗+Rhu−t²∆Rhu−fk²_0,K)^1/2o .

Proof. Letq^∗∈Q^∗_h. From Theorem 4.3 it follows that we have a pair(v, r^∗)∈ V_h×Q^∗_h such thatkvk_t,h+|||r^∗|||_t,h≤Cand

kuh−Rhukt,h+|||p^∗_h−q^∗|||t,h≤CB_h^∗(uh−Rhu, p^∗_h−q^∗;v, r^∗).

Combining the denition of the postprocessed problem and the consistency re- sult 3.1 gives

kuh−Rhukt,h+|||p^∗_h−q^∗|||t,h≤CB_h^∗(u−Rhu, p−q^∗;v, r^∗)−(g−Phg, r^∗)

=a_h(u−R_hu,v) + (divv, p−q^∗) + (div (u−R_hu), r^∗)−(g−P_hg, r^∗) X

K∈Kh

h²_K

h²_K+t²(−∇(p−q^∗) + (u−Rhu)−t²∆(u−Rhu),∇(I−Ph)r^∗)K. The last two terms on the second line cancel by the commuting diagram prop- erty (3.30). Insertingf into the last equation we have

ku_h−R_huk_t,h+|||p^∗_h−q^∗|||_t,h≤C{ku−R_huk_t,hkvk_t,h+|||p−q^∗|||_t,hkvk_t,h + ( X

K∈Kh

h²_K

h²_K+t²k∇q^∗−Rhu+t²∆Rhu+fk²_0,K)^1/2|||r^∗|||t,h. Thus the assertion is proved using the triangle equality and the above result.

Assuming full regularity, we have the following optimal a priori result for the postprocessed problem.

Theorem 4.5. For the solution (u_h, p^∗_h) of the postprocessed problem (4.5), assuming sucient regularity of the solution, it holds

ku−u_hk_t,h+|||p−p^∗_h|||_t,h≤

(C(h^k+th^k−1)kuk_k, for RT,

C(h^k+1+th^k)kukk+1, for BDM. (4.14)

5 A posteriori estimates

In this section we derive a residual-based a posteriori estimator for the post- processed solution. It should be noted that the postprocessing is vital for a properly functioning estimator. Our derivation of the a posteriori estimator is based on the following saturation assumption. Let K_h/2 be a uniformly rened subtriangulation ofKh, then the saturation assumption is [6, 18]

ku−u_h/2k_t,h/2+|||p−p^∗_h/2|||_t,h/2≤β(ku−uhkt,h+|||p−p^∗_h|||t,h), (5.1) where the constantβ <1. The triangle inequality gives

ku−u_hk_t,h+|||p−p^∗_h|||_t,h≤ 1

1−β(ku_h/2−u_hk_t,h/2+|||p^∗_h/2−p^∗_h|||_t,h/2). (5.2) We divide the estimator into two distinct parts, one dened over the elements and one over the edges of the mesh. The elementwise and edgewise estimators are dened as

η²_K= h²_K

h²_K+t²k −t²∆uh+uh− ∇p^∗_h−fk²_0,K+ (t²+h²_K)kg−Phgk²_0,K, (5.3) η_E² = t²

h_Ek[[u_h]]k²_0,E+ h_E

h²_E+t²k[[p^∗_h]]k²_0,E+ h_E

h²_E+t²k[[t²∂u_h

∂n]]k²_0,E. (5.4) The global estimator is

η= X

K∈Kh

η_K² + X

E∈Eh

η²_E

!_1/2

. (5.5)

Note that settingt= 0gives the standard estimator for the Darcy problem, see e.g. [16, 18]. In the following, we address the reliability and eciency of the estimator and show the terms of the estimator to be properly matched to one another.

5.1 Reliability

First we focus on the reliability and prove the following theorem. Note that the upper bound holds uniformly for all values of the parametertwith the constant C independent oft.

Theorem 5.1. Suppose that the saturation assumption holds. Then there exists a constant C >0 such that

ku−u_hk_t,h+|||p−p^∗_h|||_t,h≤Cη. (5.6) Proof. Due to (5.2) we only have to prove the result

ku_h/2−uhk_t,h/2+|||p^∗_h/2−p^∗_h|||_t,h/2≤Cη. (5.7) By the stability result of Theorem 4.3 we can nd (v, q^∗)∈V_h/2×Q^∗_h/2 such that kvk_t,h/2+|||q^∗|||_t,h/2≤Cfor which it holds

ku_h/2−u_hk_t,h/2+|||p^∗_h/2−p^∗_h|||_t,h/2≤ B_h/2^∗ (u_h/2−u_h, p^∗_h/2−p^∗_h;v, q^∗). (5.8)

Next, we add and subtractR_hv∈V_h andP_hq^∗∈Q^∗_h from the test functionsv and q^∗. For the projections P_h/2 and P_h it holds P_h/2P_h =P_h, which implies (I−P_h/2)(I−P_h) = (I−P_h/2). This property will be important in the analysis to follow. We introduce the following notation:

B_h/2^∗ (u_h/2−u_h, p^∗_h/2−p^∗_h;v, q^∗) =I+II, (5.9) in which

I=B^∗_h/2(u_h/2−uh, p^∗_h/2−p^∗_h;v−Rhv, q^∗−Phq^∗), (5.10) II=B^∗_h/2(u_h/2−u_h, p^∗_h/2−p^∗_h;R_hv, P_hq^∗). (5.11) Keeping in mind that(u_h/2, p^∗_h/2)is the solution on the rened mesh, we have I=L_h/2(f, P_h/2g;v−R_hv, q^∗−P_hq^∗)− B_h/2(u_h, p^∗_h;v−R_hv, q^∗−P_hq^∗)

=Ia+Ib+Ic, in which

Ia= (f,v−Rhv)−ah/2(uh,v−Rhv)−(div(v−Rhv), p^∗_h), (5.12) Ib= (P_h/2g, q^∗−Phq^∗)−(divuh, q^∗−Phq^∗), (5.13) Ic= X

K∈Kh

h²_K

h²_K+ 4t²(∇p^∗_h−uh+t²∆uh+f,∇(I−P_h/2)q^∗). (5.14) We have the following interpolation estimate forR_hv∈V_h:

h²_K+t²

h²_K kv−Rhvk²_0,K ≤C(kvk²_0,K+t²k∇vk²_0,K)≤ kvk²_t,h. (5.15) The following shorthand notation for the residual is used

R_K = (−t²∆u_h+u_h− ∇p^∗_h−f)|_K. (5.16) To estimate the termI_a we rst integrate by parts in the rst and second term.

This gives Ia= X

K∈Kh

{(t²∆uh−uh+∇p^∗_h+f,v−Rhv)K−t² 2h[[∂u_h

∂n]],v−Rhvi∂K}

+ X

E∈Eh

{t²h∂(v−Rhv)

∂n ,[[uh]]iE−αt²

h_Kh[[uh]],[[v−Rhv]]iE

− h(v−R_hv)·n,[[p^∗_h]]i_E}.

Using the inequality (3.13), scaling arguments, and the inequality (5.15) the

termI_a can be estimated as Ia≤ X

K∈Kh

{ h²_K

h²_K+t² 1/2

kRKk0,K

h²_K+t² h²_K

1/2

kv−Rhvk0,K

+ hK

h²_K+t² _1/2

k[[t²∂uh

∂n ]]k_0,∂K

h²_K+t² h_K

_1/2

h^−1/2_K kv−R_hvk_0,K

+tk∇(v−R_hv)k_0,K t

h^1/2k[[u_h]]k_0,∂K}

+ X

E∈Eh

{ t

h^1/2_K k[[uh]]k0,E t

h^1/2_K k[[v−Rhv]]k0,E

+ h_K

h²_K+t² _1/2

k[[p^∗_h]]k_0,E

h²_K+t² h_K

_1/2

h^−1/2_K kv−R_hvk_0,K}

≤C X

K∈Kh

h²_K

h²_K+t²kRKk²_0,K+ X

E∈Eh

{ t²

h_Kk[[uh]]k²_0,E

+ hK

h²_K+t²k[[t²∂uh

∂n]]k²_0,E+ hK

h²_K+t²k[[p^∗_h]]k²_0,E} _1/2

kvkt,h.

Turning to the termI_b, we have by the equilibrium property (3.6) the result divu_h=P_hg.

Adding and subtracting the loading ggives

I_b= (P_h/2g−div u_h, q^∗−P_hq^∗) = (P_h/2g−g+g−P_hg,(I−P_h)q^∗)

≤ X

K∈Kh

{kPh/2g−gk0,K+kPhg−gk0,K}k(I−Ph)q^∗k0,K

≤C X

K∈Kh

(t²+h²_K)^1/2kPhg−gk0,K

h²_K h²_K+t²

_1/2

k∇q^∗k0,K

≤C X

K∈Kh

(t²+h²_K)kP_hg−gk²_0,K

!_1/2

|||q^∗|||_t,h.

Finally, for the termIc we have by straightforward estimation and the inequal- ity (4.7)

Ic≤ X

K∈Kh

h²_K h²_K+t²

1/2

kRKk0,K

h²_K h²_K+t²

1/2

k∇(I−P_h/2)q^∗k0,K

≤ X

K∈Kh

h²_K

h²_K+t²kRKk²_0,K

!_1/2

|||(I−P_h/2)q^∗|||_t,h/2

≤C X

K∈Kh

h²_K

h²_K+t²kR_Kk²_0,K

!_1/2

|||q^∗|||_t,h/2.

Combining the above results we have

I=I_a+I_b+I_c≤Cη. (5.17)

Employing the fact that(I−P_h/2)P_h= 0, we have for the second term II =B_h/2^∗ (u_h/2, p^∗_h/2;Rhv, Phq^∗)− B_h/2^∗ (uh, p^∗_h;Rhv, Phq^∗)

=B_h^∗(uh, p^∗_h;Rhv, Phq^∗)− B_h/2^∗ (uh, p^∗_h;Rhv, Phq^∗)

− L_h(P_hg,f;R_hv, P_hq^∗) +L_h/2(P_h/2g,f;R_hv, P_hq^∗)

=t² X

E∈Eh

α

hKh[[u_h]],[[R_hv]]i_E−t² X

E∈Eh

2α

hKh[[u_h]],[[R_hv]]i_E

+ X

K∈Kh

h²_K

h²_K+t²(R_K,∇(I−P_h)P_hq^∗)_K

− X

K∈Kh

h²_K

h²_K+ 4t²(RK,∇(I−P_h/2)Phq^∗)K

+ (P_hg, P_hq^∗)−(P_h/2g, P_hq^∗)

=−t² X

E∈Eh

α

h_Kh[[uh]],[[Rhv]]iE+ (g, P_h²q^∗−P_h/2Phq^∗)

≤ X

E∈Eh

αt

h^1/2_K k[[u_h]]k_0,E t

h^1/2_K k[[R_hv]]k_0,E

≤C X

E∈Eh

t²

h_Kk[[u_h]]k²_0,E

!_1/2 kvk_t,h.

Combining the estimates for partsI andII gives

B_h/2^∗ (u_h/2−uh, p^∗_h/2−p^∗_h;v, q^∗)≤Cη, (5.18) and thus the theorem holds.

5.2 Eciency

Showing the estimator to be ecient proves to be more tedious than for the case of the pure Darcy ow treated in [18]. Indeed, we have to resort to the standard bubble function techniques to obtain the desired result. In the following we denote byω_E the union of elements sharing an edge or faceE. In addition, two cut-o functionsΨ_K andΨ_E are introduced. The functionΨ_K has its support in Kand0≤ΨK≤1, whilst ΨE is supported inωE and0≤ΨE≤1. Finally, we need an extension operatorχ:L²(E)→L²(ωE)such that on the edgeE it coincides with the identity operator. We have the following lemma [25].

Lemma 5.2. For an element K with an edge E it holds, for any polynomials pandσ,

kΨKpk0,K≤ kpk0,K≤CkΨ^1/2_K pk0,K, k∇(Ψ_Kp)k_0,K≤Ch⁻¹_K kΨ_Kpk_0,K,

kσk_0,E ≤CkΨ^1/2_E pk_0,E,

Ch^1/2_E kσk_0,E ≤ kΨ_Eχσk_0,K≤Ch^1/2_E kσk_0,E, k∇(ΨEχσ)k0,K≤Ch⁻¹_K kΨEχσk0,K.