A new divergence-free interpolation operator with ... - Semantic Scholar

Comment

Report 1 Downloads 50 Views

A new divergence-free interpolation operator with applications to the Darcy-Stokes-Brinkman equations Xuejun Xu and Shangyou Zhang LSEC, Institute of Computational Mathematics, Chinese Academy of Sciences, PO Box 2719, Beijing 100080, China Department of Mathematical Sciences, University of Delaware, DE 19716, USA. [email protected]

July 13, 2009 Abstract A local interpolation operator, preserving the divergence, is constructed explicitly, which is the first ever one for the divergence-free elements. A divergence-free finite element method is applied to the Darcy-Stokes-Brinkman flow in a mixed region of both free and porous medium. The method is of optimal order in contrast to the traditional Hdiv and H1 mixed finite elements, which work for either the Darcy flow or the Stokes flow but not the mixture. Comparing to the existing non-conforming elements, the divergence-free element method provides a continuous solution for the velocity which is also an orthogonal projection within a Hilbert subspace of the true velocity. Numerical tests supporting the theory are presented. Keywords. mixed finite element, Stokes equations, Darcy’s law, Brinkman equations, divergence-free element. AMS subject classifications (2000). 65M60, 65N30, 76M10, 76D07.

1

Introduction

We consider a model Darcy-Stokes-Brinkman mixed flow problem (cf. [10, 19]): Find the velocity u and the pressure p such that  2   −ǫ ∆u + u + ∇p = f in Ω, (1.1) div u = g in Ω,   u = 0 on ∂Ω. R where Ω is a bounded polygonal domain in R2 , and Ω g = 0. When ǫ = 0, (1.1) becomes the Darcy’s law for porous medium flow. When ǫ = 1, (1.1) is a form of Stokes equations, modeling free fluid flow. When the singular perturbation parameter ǫ varies between 0 and 1 in the domain, the partial differential equations (1.1) change the type. The Darcy-StokesBrinkman equations can be applied to fuel cell dynamics, for example, which models the flow of gas and water in a multi-domain by one set of equations on a single domain, cf. [19]. As pointed out in [10], the traditional Hdiv mixed elements for 2nd order elliptic equations (when ǫ = 0) and H1 mixed elements for the incompressible fluids and nearly incompressible 1

materials (when ǫ = 1) do not work for ǫ both close to 0 and 1. In most mixed finite elements for the Stokes equations, to satisfy the inf-sup condition, the finite element space for the velocity is chosen too big in the sense that there is a (lot) vh the finite element velocity space Vh such that k div vh kL2 6= 0 but (div vh , qh ) = 0 (1.2) for all qh in the finite element pressure space Ph . In this case, a non-divergence free part of finite element solution uh is not controlled by the second equation of (1.1) due to (1.2), while not controlled well by the first equation of (1.1) either if ǫ = 0. It is controlled only by the lower order term in the first equation of (1.1). This results to k div uh − gkL2 6→ 0 and even k div uh kL2 → ∞, usually. Two such examples are provided in [10]. On the other side, in most mixed element methods for the Darcy’s law, the finite element velocity space is not even an H1 space, neither an H1 -nonconforming space, but only an Hdiv space. Such mixed elements would obviously not work for the Stokes equations. For a pair of finite elements spaces to work uniformly for the problem (1.1) with ǫ ∈ [0, 1], a compatible condition is, implied by [10] too, div Vh = Ph .

(1.3)

Of course, the div in (1.3) can be interpreted in a weak sense, i.e., div is a discrete divergence operator. In other words, (1.2) would not happen. So far, very few such element pairs are constructed [10, 19], satisfying (1.3), for the Darcy-Stokes-Brinkman system (1.1), and they are all H 1 non-conforming elements. For the Stokes or Navier-Stokes equations, rewritten in variational forms, the primitive unknowns, the velocity and the pressure, belong to Sobolev spaces H1 and L2 , respectively. dc element which approximates the veNaturally, a finite element method would be the Pk -Pk−1 1 locity in an H -subspace of continuous piecewise Pk polynomials and approximates the pressure in an L2 -subspace of discontinuous Pk−1 polynomials. This is a truly conforming element as the incompressibility condition is satisfied pointwise and the discrete solution for the velocity is a projection within the space of divergence-free functions. We call such H1 finite elements divergence-free elements. A fundamental study on the method was done by Scott and Vogelius ([14, 15]) that the method is stable and consequently of the optimal order on 2D triangular grids for any k ≥ 4, provided the grids have no nearly-singular vertex. Some other divergencefree elements are also discovered, cf. [1, 3, 4, 12, 20, 21, 22]. To be divergence-free elements, the necessary and sufficient condition is (1.3). In principle, all divergence-free elements would approximate the solutions of the Darcy-Stokes-Brinkman system in optimal order, uniformly for all ǫ ∈ [0, 1]. In this work, however, we only analyze one type of divergence-free element, the 2D P2 Hsieh-Clough-Tocher divergence-free element. This element was analyzed by Qin and shown to be stable for the Stokes equations, in an unpublished work [12]. Nevertheless, our analyze is not based on Qin’s work. We constructed a first-ever divergence-preserving interpolation operator, extending the Scott-Zhang operator [16]. Here the divergence-preserving is meant “pointwise” in the following sense: P0 div u(x) = div Ih u(x),

(1.4)

where P0 is the L2 -orthogonal projection operator to the the image space of divergence of discrete velocity space. In particular, div Ih u = 0 if div u = 0. Previously, such operators were constructed by Girault and Scott [9], and some others [10, 19], for some finite elements but not for divergence-free elements. A divergence-preserving operator would imply the infsup condition for the finite element pair (Vh , Ph ). Therefore, we re-produce Qin’s result [12] 2

as a byproduct. With this divergence-preserving operator, the uniform convergence for (1.1) by the P2 Hsieh-Clough-Tocher divergence-free element would be deduced. Using the same technique, one can construct divergence-preserving operators for higher order Hsieh-CloughTocher divergence-free elements on triangular grids, for example, the P3 HCT element in [12], as well as to Pk HCT divergence-free elements on tetrahedral grids for any k ≥ 3, cf. [20]. But it seems this kind of construction would be possible only for such a few macro-elements. Finally, we provide numerical results supporting the theory.

2

The divergence-free element

Multiplying the Darcy-Stokes-Brinkman equations (1.1) by test functions and integrating by parts, the following variational form for the problem is derived: Find u ∈ H10 (Ω) and p ∈ L20 (Ω) such that aǫ (u, v) − (p, div v) = (f , v) ∀v ∈ H10 (Ω), (2.1) (q, div u) = (g, q) ∀q ∈ L20 (Ω), where H10 = H01 (Ω)d (d = 2, 3) is the vector Sobolev space with 0 trace on ∂Ω, L20 (Ω) is the R space of square integrable functions with mean Ω (·) dx = 0, and Z aǫ (u, v) = ǫ2 ∇u∇vdx + (u, v) ∀u, v ∈ H 1 (Ω)d , d = 2, 3, Z Ω uvdx ∀u, v ∈ L2 (Ω)d , d = 1, 2, 3. (u, v) = Ω

The problem (2.1) is well posed (cf. [10, 13]) when ǫ ∈ (0, 1]. By [10], it is well posed for ǫ = 0 too. Let Ω be a 2D polygonal domain, triangulated into quasiuniform grids T˜h (cf. [7]) of triangles, with a grid size h. We consider one type divergence-free element. The theory established here would be generalized to all divergence-free elements in 2D and 3D, though. In 2D, cf. Figure 1, we connect the bary-center of each triangle M in T˜h to its 3 vertices. The refined grid is called a Hsieh-Clough-Tocher grid (cf. [7, 12]): Th = {KM,i | M = ∪3i=1 KM,i

∀M ∈ T˜h }.

Then the P2 -P1dc mixed finite element space for the 2D problem (2.1) is defined as follows, to be referred by P2 Hsieh-Clough-Tocher divergence-free element, (2.2) Vh = vh ∈ H10 (Ω) | vh |K ∈ P2 (K)2 ∀K ∈ Th , 2 (2.3) Ph = qh ∈ L0 (Ω) | qh |K ∈ P1 (K) ∀K ∈ Th . The following uniform (in h) inf-sup condition is shown, via the macro-element technique of Stenberg [18], by Qin in [12]: inf

sup

qh ∈Ph vh ∈Vh

(qh , div vh ) ≥ C. kvh kH1 kqh kL2

(2.4)

However, we will show a stronger version inf-sup condition in this manuscript, for our singular perturbation problem (2.1), by constructing a divergence-preserving, polynomial preserving, 3

@ @

@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ ⇒ @ @ @ @ @ @

H H A HH A@HH @ H A@ H A@ A @ AA A @ AA AH @ AH @ H A H A H@ H@ A H H A @ @ H H A@H A H @ A@HH A@HH A A @ A A @ AA AH @ AH @ H A H A H@ H@ A H H A @ @

T˜h

Th

Figure 1: A uniform triangulation T˜h splits into a Hsieh-Clough-Tocher grid Th . and boundary-condition preserving interpolation operator. For Vh -Ph defined in (2.2)–(2.3), by the inf-sup condition (2.4), it follows that, i.e., (1.3) holds, Ph = div Vh = {div wh | wh ∈ Vh }.

(2.5)

For any qh = div wh ∈ Ph , we have, by the divergence theorem, Z Z Z (wh · ns )ds = 0, div wh dx = qh dx = Ω

∂Ω

Ω

as wh ⊂ H10 . Here ns denotes the unit normal to ∂Ω. Therefore, the mixed pair of finite element spaces is conforming: (2.6) Vh ⊂ H10 , Ph ⊂ L20 . Furthermore, the divergence free finite element space is also a subspace of the continuous divergence free space, Zh ⊂ Z, (2.7) where Zh = {vh ∈ Vh | div vh = 0} = {vh ∈ Vh |

Z

div vh qh dx = 0 ∀qh ∈ Ph },

Ω

Z = {v ∈ H10 (Ω) | div v = 0}.

The finite element approximation problem for (2.1) is then: Find uh ∈ Vh and p ∈ Ph such that aǫ (uh , vh ) − (ph , div vh ) = (f , vh ) ∀vh ∈ Vh , (2.8) (qh , div uh ) = (g, qh ) ∀qh ∈ Ph .

3

Divergence-preserving interpolation

The main work is in this section, constructing an interpolation operator from H10 space to the discrete velocity space Vh , which preserves the divergence. Such operators were constructed before by Girault and Scott [9] for Taylor-Hood elements and (P2 +B3 )-P1 triangular elements. Straightforwardly, one would hope to interpolate a velocity u = hu1 , u2 i to preserve its gradient. But it is not possible. In general, it is not even possible to interpolate u to preserve ∂x u1 and ∂y u2 only. Here, we preserve ∂x u1 + ∂y u2 pointwise in the sense of (1.4). As we are unable to prove the existence of such an Ih abstractly, we give a constructive proof where the explicit 4

dual basis functions are calculated. For example, instead of showing the invertibility of an 8 × 8 matrix A, we compute its inverse and listed the entries in Figure 4. The mixed finite element spaces Vh and Ph in (2.2) and (2.3) are decomposed into global and local spaces: ˜h ⊕V ¯ h, Vh = V

Ph = P˜h ⊕ P¯h ,

where n o ˜ h = vh ∈ H01 (Ω)2 | vh |M ∈ P2 (M ) ∀M ∈ T˜h , V n o ¯ h = vh ∈ Vh | vh |∂M = 0 ∀M ∈ T˜h , V n o P˜h = qh ∈ L20 (Ω) | qh |M ∈ P0 (M ) ∀M ∈ T˜h , Z qh dx = 0 ∀M ∈ T˜h . P¯h = qh ∈ Ph |

(3.1) (3.2) (3.3) (3.4)

M

˜ h is the C 0 -P2 space on the original grid T˜h , V ¯ h is the space local macroIn other words, V ˜ P2 bubbles on each big triangle M ∈ Th (which is refined into 3 triangles of Th ), P˜h is the traditional P0dc space on T˜h , and P¯h is the P1dc space of functions with mean value 0 on each macro-triangle M of the base grid T˜h . Here P1dc stands for the space of discontinuous, piecewise linear polynomials. We define the divergence-preserving interpolation operator in two steps. First, we define ˜ h . This would be standard as described in [6], an operator ˜Ih from H01 (Ω)2 to the subspace V for example, if the function to be interpolated is slightly regular than H1 so that its nodal values are well defined. To be general, we define ˜Ih by the Scott-Zhang operator [16] which uses the edge averages as the nodal values of interpolant so that the homogeneous boundary condition and polynomial functions are preserved. For each vertex xi of triangulation T˜h , we choose an edge Exi from T˜h which is any one edge if xi is an internal vertex, or one boundary edge, cf. Figure 2, (3.5) Exi ⊂ ∂Ω if xi ∈ ∂Ω. ˜ h on Let {φj,xi (x), j = 1, 2, 3} be the three standard Lagrange P2 nodal basis functions of V the edge Exi . In particular, φ1,xi is numbered to be the P2 nodal basis at the vertex xi . Then on the edge Exi we have a dual basis {φj,xi ∈ P2 (Exi ), j = 1, 2, 3} which are defined uniquely by ( Z 0 if j 6= k, φj,xi (xs )φk,xi (xs )ds = 1 if j = k. Exi In fact, it is easy to get that φ1,xi (xs ) = h−1 30(s/h)2 − 36(s/h) + 9 , 2,xi

φ

3,xi

φ

−1

(xs ) = h

−1

(xs ) = h

2

−15(s/h) + 15(s/h) − 3/2 , 2

30(1 − s/h) − 36(1 − s/h) + 9 ,

Exi = {xs = xi +

(s/h)(x′i

− xi ) | 0 ≤ s ≤ h},

(3.6) (3.7) (3.8) (3.9)

where x′i is the other end point of edge Exi , s is an arc-length parameter and h is the length of the edge, cf. Figure 2. Therefore, we uniquely define the nodal value (a vector) of ˜Ih v for 5

xi = x0 , x′i = xh

2 xi A E : m i A w A x1 mi A i A QQ Q Q Q Q Q Q Q

′ xi E : A xi A x i A w A Q Q Q Q Q Q Q Ex = {xs } Q i xi w Q x′i

A A

A

˜ h node, cf. (3.9) and (3.11) Figure 2: Choose an edge for each V v ∈ H01 (Ω)2 at node xi : (˜Ih v)(xi ) =

Z

v(xs )φ1,xi (xs )ds.

(3.10)

Exi

Next we define the mid-edge values of ˜Ih v(mi ) so that the divergence of v is preserved on each macro-triangle M ∈ T˜h . For each mid-edge point mi of grid T˜h , let Emi be the corresponding edge, cf. Figure 2. Using the same notations above, we let {φj,mi (x), j = 1, 2, 3} ˜ h on the edge Em , where φ2,m be the three standard Lagrange P2 nodal basis functions of V i i is the P2 nodal basis at the mid-point mi . We define, for each v ∈ H01 (Ω)2 , R 1 )φ 1 (x ) − (˜ 2 )φ 2 (x ) φ ˜ v(x ) − ( I v)(x I v)(x s s s 2,mi (xs )ds h h 1,x 1,xi i i E mi R i , (3.11) (˜Ih v)(mi ) = Em φ2,mi (xs )ds i

where x1i and x2i are the two end points of the edge Emi , cf. Figure 2. Here (3.11) is well defined, as the integral in the denominator is nonzero.

Lemma 3.1 The following linear interpolation operator preserves the boundary condition, the ˜ h functions and the divergence element-wise, and it is H1 -stable: V X X ˜Ih v = (˜Ih v)(mi )φ2,mi , v ∈ H01 (Ω)2 , (3.12) (˜Ih v)(xi )φ1,xi + mi ∈T˜h

xi ∈T˜h

see (3.10) and (3.11). That is,

Z

M

˜ h, ˜Ih : H1 (Ω) → V Z0 div vdx ∀M ∈ T˜h , div ˜Ih vdx =

˜Ih vh = vh

˜ h, ∀vh ∈ V

k˜Ih vkH1 ≤ CkvkH1 .

M

Proof. Since we choose a boundary edge for each boundary vertex and for each mid-edge node on the boundary in (3.5), the averaged values are (˜Ih v)(xi ) = 0, (˜Ih v)(mj ) = 0,

6

when xi , mj ∈ ∂Ω.

˜ h . Next, the dual basis preserves the nodal values of P2 Hence ˜Ih v|∂Ω = 0 and ˜Ih v ∈ V polynomial, cf. (3.6)–(3.8), i.e.,   Z X X ˜Ih vh (xi ) = φ1,xi  vh (xj )φxj + vh (mj )φmj  ds = vh (xi ), Exi

˜Ih vh (mi ) =

Z

E mi

Thus we have

xj ∈T˜h



φ2,mi 

X

mj ∈T˜h

vh (xj )φxj +

xj ∈T˜h

X

mj ∈T˜h

˜Ih vh = vh



vh (mj )φmj  ds = vh (mi ).

˜ h. ∀vh ∈ V

By mapping the dual basis to [0, 1] and taking the averaging process on [0, 1], it is straightforward to show the stability, cf. [16], that X X k˜Ih vk2H 1 (Ω)2 ≤ (˜Ih v)(xi )2 kφ1,xi k2H 1 (M )2 M ∈T˜h

+

xi ∈T˜h

(˜Ih v)(mi )2 kφ2,mi k2H 1 (M )2

X

kφ1,xi k2L2 (Ex ) kvk2L2 (Ex ) + C

X

kvk2L2 (Ex ) + C

X

kvkL2 (M ) |v|H1 (M ) ≤ CkvkL2 (Ω) |v|H1 (Ω) ≤ Ckvk2H 1 (Ω)2

mi ∈T˜h

≤C

i

xi ∈T˜h

≤C

i

i

xi ∈T˜h

≤C

X

X

X

kφ2,mi k2L2 (Em ) kvk2L2 (Em

mi ∈T˜h

kvk2L2 (Em

mi ∈T˜h

i

i)

i)

(3.13)

M ∈T˜h

for all v ∈ H01 (Ω)2 . Finally, for element-wise divergence preserving, it follows by the divergence theorem and (3.11): Z = =

div(˜Ih v)dx = M

3 Z X

i=1 Emi 3 Z X i=1

E mi

Z

(˜Ih v) · nds = ∂M

3 Z X i=1

(˜Ih v) · ni ds E mi

[(˜Ih v)(x1i ) · ni ]φx1i + [(˜Ih v)(mi ) · ni ]φmi + [(˜Ih v)(x2i ) · ni ]φx2i ds Z v(xs ) · ni ds =

div vdx. M

The second step in constructing the interpolation operator is to define the internal values for a piecewise P2 function on a macro Hsieh-Clough-Tocher triangle, which consists of three triangles in Th . This is done first on the reference Hsieh-Clough-Tocher triangle, i.e., the unit right triangle at the origin, shown in Figure 3. We need to define eight internal nodal values of (Ih v). There are two components of (Ih v) at 4 internal P2 nodes, see Figure 3. On the other 7

x3

A@ A@ A @ x7As @ A @ @ A @ A As x4 @ HH @ HHs @ x5 s x6 HH@ HH x1 @ x2

Figure 3: The reference Hsieh-Clough-Tocher triangle and its nodes. side, there are precisely eight conditions to preserve the divergence of v against Ph functions, since the dimension of the pressure space P¯h on one macro triangle is 8 (3 P1 functions would give a dimension of 9, but the condition of zero mean value would eliminate one dimension.) This is how Qin showed the inf-sup condition for the P2 Hsieh-Clough-Tocher P2 -P1 divergence free element in [12]. But we would build directly the interpolation operator by constructing explicitly a dual basis. ˆ be numbered Let the internal nodes of C 0 -P2 polynomials on the reference macro-triangle K as in Figure 3, and {φˆi , i = 4, 5, 6, 7} the four internal P2 nodal basis functions. One such nodal basis function is plotted in Figure 5. Let A be the 8×R 8 matrix representing the vector semi-H1 inner product: its 2×2 blocks of 4×4 submatrics are Kˆ ∂l φˆi ∂m φˆj dˆ x, i, j = 4, 5, 6, 7, l, m = 1, 2. By the analysis in [12], i.e., the inf-sup condition, it follows that A is invertible. But here we compute A−1 explicitly, in order to construct a semi-H1 dual basis for the four vector nodal ˆ The entries of A−1 are displayed in Figure 4. basis functions on K. ˆi x)] by, cf. Figure 3 and Figure 4, ˆ i (x) = [ˆbi (ˆ We defined the dual basis b 1 x), b2 (ˆ ! 7 ˆi,x X φ i j ˆ (ˆ φˆj (ˆ x), i = 4, 5, 6, 7, (3.14) b x) = φˆi,y j=4

j

where the coefficients are displayed in Figure 4. Two such dual basis functions are depicted in ˆ i is a dual basis, i.e., Figure 5. One can verify that b Z

ˆ K

∂ φˆj ˆ i (ˆ (ˆ x) div b x)dˆ x = δji . ∂x ˆ

(3.15)

We note that we only need half of the matrix A−1 . In (3.14), the entries of first four columns of A−1 are used. We may use the second four columns of A−1 to define another set of dual basis, then the partial directive in x ˆ direction in (3.15) would be replaced by ∂/∂ yˆ if the other dual basis is used. In fact, as the reference element is symmetric in x ˆ and yˆ, by switching i ˆ the two components of b , we would get the other set of dual basis functions. As we used the inverse matrix A−1 , we get immediately that Z Z ∂yˆφˆj (∂xˆˆbi2 + ∂yˆˆbi1 )dˆ x = δij . (3.16) ∂xˆ φˆj (∂xˆˆbi1 + ∂yˆˆbi2 )dˆ x = δij , ˆ K

ˆ K

8

7,x 6,x A@ A@ @ @ φˆ i φˆ i A A 7,y @ @ A A φˆ φˆ6,y A @ i A @ i At @ At @ @ @ A A 5 1 @ @ A A 12 6 −1 A −7 A 2@ 5@ 24 A A 5,x 4,x 24 9 @ 9 @ ˆ −5 −1 At −5 @ 1 @ φˆ i φ t A i H 18 H 18 A@ ˆ5,y A@ ˆ4,y @ @ @ @ 12 6 H H A A φi φi −1 1 H @ H @ H H @ @ A A 24 t 8 t 1 Ht 11 Ht HH@ HH@ @ 6 12 A @ A −1 −11 @ H H@ @ @ At At HH HH 8 24 @ @ @ @ A A −1 2 @ @ 12 A 9 A 1 −5 A 8@ A 2@ 8 18 A A 9 @ 9 @ 5 At 1 @ 2 At −4 @ 18 HH HH9 @ @ 12 9 1 1 HH @ HH @ 24 t 18 t 1 Ht 5 Ht H @ H @ 6 9 H HH@ −1 −5 H@ HH HH 24 18 @ @ Figure 4: The nodal values of semi-H1 dual basis, cf. (3.14) and Figure 5.

Nodal basis φ4

4,y

4,x

Dual basis φ

Dual basis φ

0.2 1 0

1 0.5

−0.2 0.5

0

0

−0.5 1

−0.4 −0.6 1

0.8

−0.5 1

0.8 1

0.8 0.8

0.6

0.6

0.4

0.6

0.4 0.2

1

0.4

0.8 0.2 0

0

0.4

1 0.8 0.6

0.2

0.4

0.2 0

0.6

0.6 0.2

0.4 0

0

0.2 0

Figure 5: A nodal basis and two dual basis functions, cf. (3.14). ˆ → M be the reference mapping, known as Piola For a general triangle M ∈ T˜h , let F : K transformation ([6]), i.e, x = F (x) =

3 X

xj φ¯j (ˆ x),

ˆ x ∈ M, ˆ ∈ K, x

j=1

where xj are the three vertices of macro-triangle M and φ¯j are the P1 nodal basis functions ˆ Let the Jacobian matrix and Jacobian of the affine transformation F be, respectively, on K. ∂x t11 t12 T = = , J = det(T ). (3.17) t21 t22 ˆ ∂x ˆ to integrals In order to find the dual basis on M , we map the averaging integrals (3.16) on K

9

on M : δij =

Z

ˆ K

∂xˆ φˆi (∂xˆˆbj1 + ∂yˆˆbj2 )dˆ x=

Z

M

(t11 ∂x φi + t21 ∂y φi )(t11 ∂x bj1

+ t21 ∂y bj1 + t12 ∂x bj2 + t22 ∂y bj2 )J −1 dx Z Z j j ˆ ˆ ˆ (t12 ∂x φi + t22 ∂y φi )(t11 ∂x bj1 ∂yˆφi (∂xˆ b1 + ∂yˆb2 )dˆ x= 0= ˆ K

M

+ t22 ∂y bj2 )J −1 dx, Z j j ˆ ˆ ˆ (t12 ∂x φi + t22 ∂y φi )(t11 ∂x bj2 ∂yˆφi (∂xˆ b2 + ∂yˆb1 )dˆ x= +

δij =

Z

ˆ K

t21 ∂y bj1

+

t12 ∂x bj2

M

+ t22 ∂y bj1 )J −1 dx, Z Z j j ˆ ˆ ˆ (t11 ∂x φi + t21 ∂y φi )(t11 ∂x bj2 ∂xˆ φi (∂xˆ b2 + ∂yˆb1 )dˆ x= 0= +

t21 ∂y bj2

+

t21 ∂y bj2

+

t12 ∂x bj1

+

t12 ∂x bj1

ˆ K

M

+ t22 ∂y bj1 )J −1 dx.

Hence, we can linearly combine them to get Z ∂x φi [∂x (t11 bj1 + t12 bj2 ) + ∂y (t21 bj1 + t22 bj2 )]dx, t22 δij = ZM ∂y φi [∂x (t11 bj1 + t12 bj2 ) + ∂y (t21 bj1 + t22 bj2 )]dx, −t12 δij = ZM ∂x φi [∂x (t11 bj2 + t12 bj1 ) + ∂y (t21 bj2 + t22 bj1 )]dx, −t21 δij = ZM ∂y φi [∂x (t11 bj2 + t12 bj1 ) + ∂y (t21 bj2 + t22 bj1 )]dx. t11 δij = M

Therefore, we conclude the calculation with the definition of of the dual basis on M : ψj,x (x) ψj,y (x) bj,1 (x) = , bj,2 (x) = , ψj,y (x) ψj,x (x)

(3.18)

where !  7 j,x j,y ˆ ˆ X  t φ + t φ 1 11 k 12 k   φk (x),  j,x  t22 ˆ ˆj,y t φ + t 21 k 22 φk k=4 ! ψj,x (x) = 7 ˆj,x + t11 φˆj,y X  t φ −1  12 k k   φk (x), t ˆj,x + t21 φˆj,y t φ 21 22 k k k=4 !  7 j,y j,x  1 X t11 φˆk + t12 φˆk   φk (x),   t11 ˆj,y + t22 φˆj,x t φ 21 k k k=4 ! ψj,y (x) = 7 ˆj,y + t11 φˆj,x X  t φ −1  12 k k   φk (x), t ˆj,y + t21 φˆj,x t φ 12 22 k k k=4

if |t11 t22 | ≥ J/2, (3.19) if |t11 t22 | < J/2, if |t11 t22 | ≥ J/2, (3.20) if |t11 t22 | < J/2,

ˆj,y where φˆj,x k and φk are defined in Figure 4, cf. Figure 3. The second step in constructing the

10

interpolation operator is ¯Ih : H 1 (Ω)2 → V ¯ h, 0 7 XX ¯Ih : v = v1 7→ ¯Ih v (xj )φj (x), v2 ˜ j=4 K∈Ω Z ∂x v1 (x) div bj,1 (x) ¯Ih v (xj ) = dx, K ∂y v2 (x) div bj,2 (x)

(3.21)

¯ h is defined in (3.2) and bj,i are defined by (3.18)–(3.20). where V

Theorem 3.1 The following linear, H1 -stable interpolation operator preserves the boundary condition, the Vh finite element functions and the divergence inside the finite element space Vh : Ih v = ˜Ih v + ¯Ih (v − ˜Ih v), v ∈ H01 (Ω)2 , (3.22) where ˜Ih and ¯Ih are defined in (3.12) and (3.21) respectively. Proof. The idea is to use the eight degrees of freedom of Ih v inside each macro-triangle M ∈ T˜h , cf. Figure 3, to match the eight nodal values of P0 div(v − ˜Ih v), where P0 is the L2 -orthogonal projector to Ph . Note that each linear P1dc function in Ph has 3 × 3 coefficients on a macro-triangle M which consists of three triangles. But div ˜Ih v matches the mean value of div v on M so that P0 div(v − ˜Ih v) has 8 independent coefficients. Since the support of ¯Ih v is inside each macro-triangle K of T˜h , Ih preserves the homogeneous ¯ h, boundary condition of v, by Lemma 3.1. For any vh ∈ Vh , by Lemma 3.1, (˜Ih vh − vh ) ∈ V defined in (3.2), i.e., (˜Ih vh − vh ) = 0 on ∂M for all M ∈ T˜h . Further ˜Ih vh (xj ) = 0 at 4 internal nodes xj , j = 4, 5, 6, 7, cf. Figure 3. Therefore (vh − ˜Ih vh )(x) =

7 X X

vh (xj )φj (x) =

7 X X CM,j φj (x), DM,j

M ∈T˜h j=4

M ∈T˜h j=4

for some constants CM,j and DM,j . Now, inside each triangle M ∈ T˜h , by (3.19) and (3.20), ¯Ih (vh − ˜Ih vh )(xi ) =

Z

7 X CM,j ∂x φj (x) div bi,1 (x)

M j=4

Hence

DM,j ∂y φj (x) div bi,2 (x)

dx =

CM,i . DM,i

Ih vh = ˜Ih vh + ¯Ih (vh − ˜Ih vh ) = ˜Ih vh + (vh − ˜Ih vh ) = vh .

(3.23)

(3.24)

Next, for showing the divergence preserving, we need to notice that for any q¯h ∈ P¯h , there is ¯ h , cf. (3.2) and (3.4), such that ¯h ∈ V aw ¯ h = q¯h . div w

(3.25)

ˆ As the matrix To show (3.25) on each element M ∈ T˜h , we map M to the reference element K. ˆ {div w ¯ h |M } is an 8-dimensional subspace of P¯h A for defining the dual basis is invertible, on K, restricted on M . On the other side, the dimension of itself, P¯h restricted on M , is 8. So (3.25)

11

¯ h restricted on each M , we have the following holds. Then, as {bj,1 , bj,2 } form a basis for V linear combination, cf. (3.18) and (3.20), ¯h = w

7 X X

w1M,j bj,1 (x) + w2M,j bj,2 (x).

(3.26)

M ∈T˜h j=4

In fact, the coefficients wiM,j in (3.26) can be easily found by (3.23) and (3.19)–(3.20), as ¯Ih w ¯h = w ¯ h . For all v ∈ H01 (Ω)2 and wh ∈ Vh , denoting div wh = q˜h ⊕ q¯h for some q˜h ∈ P˜h ¯ h ∈ P¯h , we have, by Lemma 3.1 and (3.26), and q¯h = div w Z div(v − Ih v) div wh dx Ω

=

=

=

Z Z Z

div(v − ˜Ih v) div wh dx −

Ω

Ω

div(v − ˜Ih v)¯ qh dx − Ω

div(v − ˜Ih v)¯ qh dx − Ω

=

Z

Ω

7 X X

M ∈T˜h j=4

cM,j,1 ∂x φj (x) + cM,j,2 ∂y φj (x) div wh dx

7 X X ¯ h dx cM,j,1 ∂x φj (x) + cM,j,2 ∂y φj (x) div w

M ∈T˜h j=4

7 7 X X X

M ∈T˜h i=4 j=4

cM,j,1 Z

Z

Z

M

∂x φj (x)w1M,i

div(v − ˜Ih v)¯ qh dx − Ω

div bi,1 (x)dx + cM,j,2

X

7 X

M ∈T˜h j=4

Z

M

∂y φj (x)w2M,i div bi,2 (x)dx

cM,j,1 w1M,j + cM,j,2 w2M,i ,

where

cM,j,1 cM,j,2

= ¯Ih (v − ˜Ih v)(xj ) =

Hence, Z

div(v − Ih v) div wh dx

Z

div(v − ˜Ih v)¯ qh dx −

R ˜Ih v)1 (y) div bj,1 (y)dy ∂ (v − x M R . ˜ M ∂y (v − Ih v)2 (y) div bj,2 (y)dy

Ω

=

Ω

7 Z X X

M ∈T˜h j=4

Z

M

∂x (v − ˜Ih v)1 (y) div bj,1 (y)dyw1M,j

∂y (v − ˜Ih v)2 (y) div bj,2 (y)dyw2M,j M Z Z ˜ ¯ h dx = 0. div(v − ˜Ih v) div w qh dx − = div(v − Ih v)¯ +

Ω

Ω

Finally, by mapping Ih v to the reference element, the same way in proving (3.13), it is standard to show the stability that kIh vkH1 ≤ CkvkH1 for all v ∈ H01 (Ω)2 , cf. [16].

12

4

The convergence analysis

We will show that the divergence-free element solution is the best approximation in the finite element subspaces, under appropriate norms, to the exact solutions of the Darcy-StokesBrinkman system (2.1), uniformly for ǫ ∈ [0, 1]. Following [10], the following parameter dependent norm is introduced to the analysis: kuk2ǫ = aǫ (u, u) + (div u, div u) ∀u ∈ H10 (Ω).

(4.1)

When ǫ = 1, k · kǫ is equivalent to the H1 -norm, and equivalent to the Hdiv -norm when ǫ = 0. We show the inf-sup condition again with respect this ǫ-dependent norm in the next lemma. Lemma 4.1 Both bilinear forms are coercive: aǫ (zh , zh ) ≥ kzh k2ǫ ∀zh ∈ Zh , (div vh , qh ) sup ≥ Ckqh kL2 ∀qh ∈ Ph . kvh kǫ vh ∈Vh

(4.2) (4.3)

Proof. We note that (4.2) holds in fact as an equality, by (2.7) and (4.1). For (4.3), cf. [13, 2], there is a v ∈ H10 (Ω) such that div v = qh and kvkH1 ≤ Ckqh kL2 . Let vh = Ih v. By Theorem 3.1, we get div vh = qh and consequently the inf-sup condition (div vh , qh ) = kqh k2L2 ≥ CkvkH1 ≥ Ckvh kH1 ≥

C kvh k2ǫ . 1+ǫ

We are ready to show the main theory, the quasi-optimality of the Hsieh-Clough-Tocher P2 divergence-free element. Theorem 4.1 Let Vh and Ph be defined in (2.2) and (2.3), respectively. The finite element solutions in (2.8) approximate the solution in the Darcy-Stokes-Brinkman system (2.1) in optimal order: k div(u − uh )kL2 = inf k div u − qh kL2 , qh ∈Ph inf ku − vh ka + inf kp − qh kL2 . ku − uh ka + kp − ph kL2 ≤ C vh ∈Vh

qh ∈Ph

(4.4)

Further, if the solution (u, p) is regular with some 0 < r ≤ 2 in (2.1), k div(u − uh )kL2 ≤ Chr kukH r+1 , r

ku − uh ka + kp − ph kL2 ≤ Ch (ǫ|u|H r+1 + |u|H r + |p|H r ) .

13

(4.5) (4.6)

Proof. By the inf-sup condition, there is a unique solution (uh , p): aǫ (u − uh , vh ) − (p − ph , div vh ) = 0 ∀vh ∈ Vh , (qh , div(u − uh )) = 0 ∀qh ∈ Ph .

(4.7) (4.8)

By the second equation, we get the best approximation for the divergence: k div(u − uh )k2L2 = (div u − qh , div(u − uh )) ≤ k div u − qh kL2 k div(u − uh )kL2 , k div(u − uh )kL2 = inf k div u − qh kL2 .

(4.9) (4.10)

qh ∈Ph

Let P0 p be the L2 orthogonal projection of the pressure p in space Ph , i.e., (qh , P0 p) = (qh , p) for all qh ∈ Ph . By (4.8), we have (P0 p − ph , div(u − uh )) = 0.

(4.11)

Since (qh , div(u − Ih u)) = 0, shown in Theorem 3.1, we get (P0 p − ph , div(Ih u − uh )) = 0. Letting vh = (Ih u − uh ) ∈ Vh in (4.7), we get aǫ (u − uh , Ih u − uh ) = 0.

(4.12)

Using the Schwartz inequality, we get the following estimates for the velocity ku − uh k2a = aǫ (u − uh , u − Ih u) ≤ ku − uh ka ku − Ih uka , ku − uh ka ≤ ku − Ih uka .

(4.13)

For the pressure approximation, we have kp − ph kL2 = kp − P0 pkL2 + k P0 p − ph kL2 .

(4.14)

By the continuous inf-sup condition, i.e., the maximal right inverse of the divergence operator, cf. [13, 2], there is a v ∈ H01 (Ω)2 such that div v = P0 p − ph

and

kvkH1 ≤ Ck P0 p − ph kL2 .

Let vh = Ih v in (4.7). By the preservation of divergence of Ih in Theorem 3.1, we get k P0 p − ph k2L2 = aǫ (u − uh , Ih v) ≤ ku − uh ka (ǫ|Ih v|H1 + kIh vkL2 ) ≤ ku − uh ka kIh vkH1 ≤ Cku − uh ka kvkH1 ≤ Cku − uh ka k P0 p − ph kL2 . Thus, by (4.14), (4.13) and Theorem 3.1, we get kp − ph kL2 ≤ kp − P0 pkL2 + Cku − uh ka ≤ inf kp − qh kL2 + C inf ku − vh ka . qh ∈Ph

vh ∈Vh

14

(4.15)

Normally, the solution of (1.1) depends on ǫ. So as ǫ approaches zero, the convergence estimates given in Theorem 4.1 will be deteriorate, especially in the case that the solution of (1.1) has boundary layers. Let Z |g(x)|2 1 1 2 dx < ∞, j = 1, 2, ..., N } H+ = {g ∈ H ∩ L0 | 2 Ω |x − xj | with the norm kgk21,+

=

kgk21

+

N Z X j=1

Ω

|g(x)|2 dx, |x − xj |2

where x1 , x2 , ..., xN denote the vertices of polygonal domain Ω. We assume the domain Ω and the data in (1.1) are regular enough so that the solutions of (1.1) satisfy ǫ2 kukH2 + ǫkukH1 +ku − u0 kL2 + kp − p0 kH 1 + 1

1

1

ǫ 2 ku0 kH1 + ǫ 2 kp0 kH 1 ≤ Cǫ 2 (kf kHcurl + kgk1,+ ),

(4.16)

where (u0 , p0 ) is the solution of (1.1) when ǫ = 0. (4.16) is proved in [10] for convex polygonal domain Ω, using the regularity results of [2]. Applying the same techniques developed in [10], we have the following uniform convergence estimate. We note that, comparing to the convergence rate in Theorem 4.1 (assuming there the regularity constant r = 1), we lose half an order of convergence in next theorem. Nevertheless, the rate for the divergence remains optimal by (4.5) of Theorem 4.1. Theorem 4.2 Let Vh and Ph be defined in (2.2) and (2.3), respectively. If (4.16) holds, then there is a constant C, independent of ǫ, such that ku − uh kL2 + ǫkcurl(u − uh )kL2 + kp − ph kL2 ≤ Ch1/2 (kf kHcurl + kgk1,+ ). Proof. We repeat the proof in [10], with Theorem 3.1. By (4.4) and (4.16), we get ku − uh kL2 ≤ Cku − Ih ukL2 ≤ Ck(I − Ih )(u − u0 )kL2 + Ck(I − Ih )u0 kL2 1/2

1/2

≤ Ch1/2 (ku − u0 kL2 ku − u0 kH1 + h1/2 ku0 kH1 ) ≤ Ch1/2 (kf kHcurl + kgk1,+ ). By (4.4), (4.5) and (4.16), 1/2

1/2

ǫkcurl(u − uh )kL2 ≤ Cǫ|u − Ih u|H1 ≤ Cǫh1/2 kukH1 kukH2 ≤ Ch1/2 (kf kHcurl + kgk1,+ ). Finally, by (4.15), (4.16) and the two estimates above, it follows that kp − ph kL2 ≤ kp − P0 pkL2 + Cku − uh ka ≤ Chkp − p0 kH 1 + Chkp0 kH 1 + Cku − uh ka ≤ Ch1/2 (kf kHcurl + kgk1,+ ).

15

(4.17)

HH HH A@ A@ H A@HH A@ H A @ A A @ A A @ A A @ A HH@A HH@A HH HH A A @ @ H H A@H A@H A@HH A@HH A @ A A @ A A @ A A @ A HH@A HH@A HH HH A A @ @

HH A@ HH @ H HH A@ A H A@ A@A A@A A@A A@A H@ A H@ A H@ H H HA H@ HA H H H HH A@H A@H @ A H A@ A@A A@A A@A A@A H@ A H@ H H HA H@ HA A H@ H H H HH A@H A@H @ A H A@ A@A A@A A@A A@A H@ A H@ A H@ H H HA H@ HA H H H HH A@H A@H @ A H A@ A@A A@A A@A A@A H@ A H@ A H@ H H HA H@ HA

Figure 6: The level 3, 4 and 5 Hsieh-Clough-Tocher grid Th in Table 1.

5

Numerical tests

We apply the Hsieh-Clough-Tocher P2 divergence-free element to two model problems of DarcyStokes-Brinkman equations, where one is regular while the other has a boundary layer. The numerical results confirm the analysis for both cases. We continue to apply the method to a chancel flow problem where the viscosity has a big jump. We repeat the computation for a model problem in [10]. We solve (1.1) on the unit square Ω = (0, 1)2 with the exact solution: sin2 (πx1 ) sin(2πx2 ) u = curl sin2 (πx1 ) sin2 (πx2 ) = π , (5.1) − sin2 (πx2 ) sin(2πx1 ) 2 (5.2) p = sin(πx1 ) − . π That is, g = 0 in (1.1) and 2 2

2 3

f = (1 + 4ǫ π )u + 2ǫ π

cos(2πx1 ) sin(2πx2 ) π cos(πx1 ) − . − cos(2πx2 ) sin(2πx1 ) 0

We take a sequence of uniform grids on the domain Ω as the base grids. Then each base grid is refined into a Hsieh-Clough-Tocher grid by connecting the barycenter of each triangle to the three vertices. Three levels of computational grids are shown in Figure 6. We define multi-level finite element spaces Vh and Ph ((2.2) and (2.3)) on such Hsieh-Clough-Tocher grids. For the resulting linear systems of equations (2.8), we solve them by the iterated penalty method, cf. [8, 6, 5, 17]. In Table 1, we list the errors of finite element solutions in various norms, and the convergence orders for the solutions. They seem to be of optimal order in convergence, confirming the theory in Theorem 4.1. That is, order 3 for the velocity solutions in maximum norm and in L2 norm while order 2 in H1 norm, and order 2 for the pressure solutions in maximum and L2 norms. These orders of convergence are independent of the perturbation parameter ǫ, because the solutions are regular here, independent of ǫ. Next, we apply the Hsieh-Clough-Tocher P2 divergence-free element to the Darcy-StokesBrinkman equations, where the solutions are not regular, i.e., with boundary layers. We solve (1.1) on the unit square Ω = (0, 1)2 with the exact solution: u = curlf1 ,

p = ǫ(1 − e−1/ǫ ) − e−x1 /ǫ , 16

(5.3)

Table 1: The errors and orders of convergence for P2 HCT divergence-free element. Level

|eu |l∞

hn

|eu |L2

3 4 5 6

1.75780 0.32259 0.04646 0.00603

− 2.4 2.8 2.9

1.11846 0.18254 0.02586 0.00316

3 4 5 6

1.72846 0.31621 0.04583 0.00596

− 2.5 2.8 2.9

1.10515 0.17821 0.02531 0.00312

3 4 5 6

1.62753 0.29508 0.04368 0.00594

− 2.5 2.8 2.9

1.06127 0.16490 0.02364 0.00299

3 4 5 6

1.44164 0.27507 0.03917 0.00590

− 2.4 2.8 2.7

0.88100 0.11814 0.01814 0.00257

3 4 5 6

1.55194 0.29178 0.04053 0.00594

− 2.4 2.8 2.8

0.83451 0.10507 0.01628 0.00243

hn

|eu |H1 hn ǫ = 1 in (2.8) − 13.53769 − 2.6 4.51836 1.6 2.8 1.44821 1.6 3.0 0.42560 1.8 ǫ = e−1 in (2.8) − 13.54040 − 2.6 4.51816 1.6 2.8 1.44822 1.6 3.0 0.42560 1.8 ǫ = 4−2 in (2.8) − 13.57410 − 2.7 4.52095 1.6 2.8 1.44842 1.6 3.0 0.42560 1.8 ǫ = 4−4 in (2.8) − 14.53164 − 2.9 4.62103 1.7 2.7 1.45392 1.7 2.8 0.42580 1.8 ǫ = 0 in (2.8) − 16.83536 − 3.0 4.91961 1.8 2.7 1.49852 1.7 2.7 0.43134 1.8

|ep |l∞

hn

kep kL2

hn

53.55211 31.45767 12.24402 3.71159

− 0.8 1.4 1.7

18.50145 8.16259 3.39324 1.19218

− 1.2 1.3 1.5

13.94648 7.93903 3.06569 0.92926

− 0.8 1.4 1.7

4.78990 2.05655 0.84939 0.29811

− 1.2 1.3 1.5

4.01395 2.05128 0.77066 0.23366

− 1.0 1.4 1.7

1.35218 0.52846 0.21333 0.07458

− 1.3 1.3 1.5

0.72949 0.18672 0.05206 0.01641

− 2.0 1.8 1.7

0.22448 0.04422 0.01426 0.00474

− 2.3 1.6 1.6

0.41060 0.04987 0.00876 0.00188

− 3.0 2.5 2.2

0.12466 0.01088 0.00209 0.00050

− 3.5 2.3 2.0

where f1 = 28 e−x1 x2 /ǫ x21 (1 − x1 )2 x22 (1 − x2 )2 . That is, g = 0 in (1.1) and f = (1 − ǫ∆)u − grad p. We note that the two components of u have each a boundary layer in each direction, shown in Figure 7. Comparing to the convergence orders in Table 1, we can see that we lose some orders of convergence here in Table 2. However, the errors and the orders of convergence shown in Table 2 are consistent to the theory shown in Theorem 4.2. Our third numerical test is on computing a flow driven by a boundary force, not by an internal force, i.e., for the problem (1.1) where f = 0, g = 0, but u|∂Ω 6= 0. The domain is shown in Figure 8, where the in-flow and out-flow boundary condition are, respectively, 0 0 u= , u= . (5.4) x(2 − x) 8x(1 − x) 17

Figure 7: The exact solution u of test problem (5.3) with boundary layers. In Figure 8, we plot the computed velocity field uh and the second component of uh , for ǫ = 1 in (1.1). The corresponding graphs for ǫ = 4−4 are nearly identical. In Figure 9, the computed pressure ph is shown for ǫ = 1 and ǫ = 4−4 . The are slightly different.

Figure 8: The uh filed and the second component of uh , cf. (5.4).

References [1] D. N. Arnold and J. Qin, Quadratic velocity/linear pressure Stokes elements, in Advances in Computer Methods for Partial Differential Equations VII, ed. R. Vichnevetsky and R.S. Steplemen, 1992. [2] D. Arnold, L.R. Scott and M. Vogelius, Regular inversion of the divergence operator with Dirichlet conditions on a polygon, Ann. Sc. Norm. Super Pisa, C1. Sci., IV Ser., 15(1988), 169–192 [3] G. Baker, W. Jureidini and A. Karakashian, Piecewise solenoidal vector fields and the Stokes problem, SIAM J. Num. Anal. 27 (1990), 1466–1485. [4] B. Bernstein, K.A. Feigl and T. Olsen, A first-order exactly incompressible finite element for Axisymmetric Fluid Flow, SIAM J. Num. Anal., 33-5 (1996), 1736–1758. 18

Table 2: The P2 divergence-free HCT element for test problem (5.3). Level

|eu |L2

2 3 4 5 6

1.0892 0.9074 0.1697 0.0223 0.0026

2 3 4 5 6

0.6322 0.5535 0.1263 0.0168 0.0019

2 3 4 5 6

0.1634 0.2300 0.1078 0.0218 0.0033

2 3 4 5 6

1.9040 0.9475 0.4342 0.1798 0.0723

hn |eu |H1 hn ǫ = 1 in (5.3) 12.01 0.26 11.09 0.12 2.42 4.16 1.41 2.93 1.31 1.67 3.12 0.37 1.84 −1 ǫ = 4 in (5.3) 6.696 0.19 7.822 2.13 3.426 1.19 2.90 1.098 1.64 3.08 0.304 1.84 ǫ = 4−2 in (5.3) 2.584 4.106 1.09 3.610 0.18 2.30 1.767 1.03 2.72 0.580 1.60 −4 ǫ = 4 in (5.3) 36.87 1.00 19.08 0.95 1.12 12.17 0.64 1.27 5.82 1.06 1.31 2.40 1.27

kep kL2

hn

24.23 15.86 9.11 3.56 1.17

0.61 0.79 1.35 1.61

4.124 2.957 1.888 0.875 0.553

0.47 0.64 1.10 0.66

0.9806 0.5823 0.4875 0.3750 0.3363

0.75 0.25 0.37 0.15

3.0790 1.2450 0.5166 0.2123 0.1035

1.30 1.26 1.28 1.03

[5] S.C. Brenner and L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 1994. [6] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer, 1991. [7] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [8] M. Fortin and R. Glowinski, Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-value Problems, North Holland, Amsterdam, 1983. [9] V. Girault and L.R. Scott, A quasi-local interpolation operator preserving the discrete divergence, Calcolo 40 (2003), 1–19. [10] K. A. Mardal, X.-C. Tai and R. Winther, A robust finite element method for Darcy-Stokes flow, SIAM J. Numer. Anal. 40 (2002), 1605–1631. [11] D.A. Nield and A. Bejan, Convection in Porous Media, Springer-Verlag, 2006. [12] J. Qin , On the convergence of some low order mixed finite elements for incompressible fluids, Thesis, Pennsylvania State University, 1994. 19

0.98E-02 0.91E-02 0.87E-02 0.83E-02 0.78E-02 0.74E-02 0.69E-02 0.65E-02 0.61E-02 0.56E-02 0.52E-02 0.47E-02 0.43E-02 0.39E-02 0.34E-02 0.30E-02 0.25E-02 0.21E-02 0.16E-02 0.12E-02 0.76E-03 0.32E-03 -0.12E-03 -0.56E-03 -0.10E-02 -0.14E-02 -0.19E-02 -0.23E-02 -0.28E-02 -0.32E-02 -0.37E-02 -0.41E-02 -0.45E-02 -0.50E-02 -0.54E-02 -0.59E-02 -0.63E-02 -0.67E-02 -0.72E-02 -0.76E-02 -0.81E-02 -0.85E-02 -0.94E-02 -0.98E-02 -0.10E-01 -0.11E-01 -0.12E-01

0.56E-04 0.51E-04 0.49E-04 0.47E-04 0.45E-04 0.43E-04 0.41E-04 0.39E-04 0.37E-04 0.35E-04 0.33E-04 0.31E-04 0.28E-04 0.26E-04 0.24E-04 0.22E-04 0.20E-04 0.18E-04 0.16E-04 0.14E-04 0.12E-04 0.10E-04 0.82E-05 0.61E-05 0.41E-05 0.21E-05 0.47E-07 -0.20E-05 -0.40E-05 -0.60E-05 -0.81E-05 -0.10E-04 -0.12E-04 -0.14E-04 -0.16E-04 -0.18E-04 -0.20E-04 -0.22E-04 -0.24E-04 -0.26E-04 -0.28E-04 -0.30E-04 -0.32E-04 -0.34E-04 -0.37E-04 -0.39E-04 -0.46E-04

Figure 9: ph plots for ǫ = 1 and ǫ = 4−4 , cf. (5.4). [13] P. A. Raviart and V. Girault, Finite element methods for Navier-Stokes equations, Springer, 1986. [14] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, RAIRO, Modelisation Math. Anal. Numer. 19 (1985), 111–143. [15] L. R. Scott and M. Vogelius, Conforming finite element methods for incompressible and nearly incompressible continua, in Lectures in Applied Mathematics 22, 1985, 221–244. [16] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions , Math. Comp. 54 (1990), 483–493. [17] L. R. Scott and S. Zhang, Multilevel Iterated Penalty Method for Mixed Elements, the Proceedings for the Ninth International Conference on Domain Decomposition Methods, 133-139, Bergen, 1998. [18] R. Stenberg, Analysis of mixed finite element methods for the Stokes problem: a unified approach, Math. Comp. 42 (1984), 9–23. [19] X. P. Xie, J. C. Xu and G. R. Xue, Uniformly-stable finite element methods for DarcyStokes-Brinkman models, J. Comput. Math. 26 (2008) 437–455. [20] S. Zhang, A new family of stable mixed finite elements for 3D Stokes equations, Math. Comp. 74 (2005), 250, 543–554. [21] S. Zhang, On the P1 Powell-Sabin divergence-free finite element for the Stokes equations, J. Comp. Math., 26 (2008), 456-470. [22] S. Zhang, A family of Qk+1,k ×Qk,k+1 divergence-free finite elements on rectangular grids, SIAM Num. Anal., 47 (2009), 2090-2107.

20