Computational considerations

CHAPTER 4. DISCRETIZATION OF THE WEAK FORMULATION 48













A⁽¹⁾₁₁ A⁽¹⁾₁₂ (B⁽¹⁾₁ )^T A^(I)₁₁ A^(I)₁₂ 0 0 A⁽¹⁾₂₁ A⁽¹⁾₂₂ (B⁽¹⁾₂ )^T A^(I)₂₁ A^(I)₂₂ 0 0 B₁⁽¹⁾ B₂⁽¹⁾ 0 B₁^(I) B₂^(I) 0 0 (A^(I)₁₁)^T (A^(I)₂₁)^T (B₁^(I))^T C₁₁^(I) C₁₂^(I) 0 0 (A^(I)₁₂)^T (A^(I)₂₂)^T (B₂^(I))^T C₂₁^(I) C₂₂^(I)+D^∗ (A^∗)^T B^∗

0 0 0 0 A^∗ A⁽²⁾ B⁽²⁾

0 0 0 0 (B^∗)^T (B⁽²⁾)^T 0













| {z }

K^h:=











 u⁽¹⁾₁ u⁽¹⁾₂ p⁽¹⁾ u^(I)₁ u^(I)₂ c p⁽²⁾













=













 f⁽¹⁾

11

f⁽¹⁾₁₂ f⁽¹⁾

2

f^(I)

11

f^(I)

12

f⁽²⁾

1

f⁽²⁾

2















(4.79) This matrix form can be displayed more clearly in the following way:









A_Stokes B_Stokes^T A^(I)_Stokes 0 0

B_Stokes 0 B_Stokes^(I) 0 0

(A^(I)_Stokes)^T (B_Stokes^(I) )^T C^(I) [0A^∗]^T [0 (B^∗)^T]^T

0 0 [0A^∗] ADarcy BDarcy

0 0 [0 (B^∗)^T] B_Darcy^T 0

















u_Stokes pStokes

u_{Interf ace} c_Darcy pDarcy







=









f1,Stokes

f2,Stokes

f^(I) f1,Darcy

f2,Darcy









(4.80) Here we want to point out that in the left upper part of our discretization matrix the 2×2 block discretization matrix for the usual Stokes problem occurs and in the right lower part the2×2block discretization matrix for the Darcy equations emerges. The other appearing matrices are referable to the coupling of the Stokes and the Darcy equations.

y

x 1

0

-1

1

bb

b bb

bbb b b

b b bcbcbcbcbc

bc

bc

bc

bcbcbcbcbc bbc

degrees of freedom for the velocity inΩ1

dirichlet nodes for the velocity onΓ₁ degrees of freedom for the velocity inΩ₂

known normal components of the velocity onΓ₂

Figure 4.3: Degrees of freedom for approximating the velocity u.

For the integral approximation the following Gauss quadrature rule of degree 2 on the reference triangle is used:

∫ 1 0

∫ 1−y 0

f(x, y)dx dy= 1

6(f(0,1

2) +f(1

2,0) +f(1 2,1

2)), (4.81) see [8]. This quadrature formula is exact for arbitrary polynomials of degree 2.

For the approximation of the appearing integrals over the edges on Γ_I in the bilinear form a₁ we rst transform them to an integral over the reference edge with length1and then we use the Newton Cotes formula (see, e.g. [11]) of degree n = 4 (Milne rule)

∫ ₁

0

f(t)dt =h

∑4 i=0

α_if(x_i) (4.82)

with weights ¹⁴₄₅,⁶⁴₄₅,²⁴₄₅,⁶⁴₄₅ and ¹⁴₄₅ and sampling points x_i =ih, i= 0, . . .4, with h= _n¹. We use the Milne rule because it is exact for polynomials of degree 4. The linear shape functions on the reference triangle are given by

N₁ = 1−x−y, (4.83)

N₂ = x, (4.84)

N₃ = y. (4.85)

CHAPTER 4. DISCRETIZATION OF THE WEAK FORMULATION 50 y

x 1

0

-1

1

bb bb

bbb b b b

*

* * *

*

*

* *

*

degrees of freedom for the pressure inΩ₁ degrees of freedom for the pressure inΩ₂

Figure 4.4: Degrees of freedom for approximating the pressure p.

The quadratic shape functions on the reference triangle are given by

N₁ = 1−3x−3y+ 4xy+ 2x²+ 2y², (4.86)

N₂ = −x+ 2x², (4.87)

N₃ = −y+ 2y², (4.88)

N₄ = 4x−4xy−4x², (4.89)

N5 = 4xy, (4.90)

N₆ = 4y−4xy−4y². (4.91)

Before we write down the explicit form of our basis functionsϕE, E∈ E2, we want to clarify the notation for elements with a common edge E ∈ E2 and determine the orientation of n_E by the next denition.

Denition 4.4. [2, Def 4.2] Let T_± =conv{E∪ {P_±}} for the vertex P_± opposite to E of T_± such that the edge E = conv{A, B} orients from A to B. Then n_E points outward fromT₊ to T₋ (see Figure 4.5) with a positive sign. If E is an exterior edge, i.e., E ∈ Eh(Γ₂)∪ Eh(Γ_I), then n₂ =n_E is the exterior normal and E ⊂∂T₊ denes T₊ (and T₋ is undened).

Note that the normal direction depends on the orientation of the edge E = conv{A, B}. Hence, we choose the orientation of E ∈ E2 in such a way that

for all interior edges it is arbitrary but xed and for the edges on the boundary it is chosen such thatn_E coincides with the exterior normal vector n₂.

Now we can move on to the denition of the edge-basis functions for the lowest-order Raviart-Thomas nite element space.

Denition 4.5. [2, Def 4.3] Given an edge E ∈ E2 there are either two elements T₊ and T₋ in T2^h with the joint edge E =∂T₊∩∂T₋ or there is exactly one element T₊ in T₂^h with E ⊂ ∂T₊. Then if T_± = conv{E∪ {P_±}} for the vertex P_± opposite to E of T_± set

ϕ_E(x, y) :=

{±₂^|_|^E_T±^|_|((x, y)−P_±), for (x, y)∈T_±,

0, elsewhere. (4.92)

Lemma 4.6. [2, Lemma 4.1] There hold 1. ϕ_E ·n_E =

{

0, along (∪E2)\E, 1, along E;

2. ϕ_E ∈H(div,Ω₂);

3. {ϕ_E :E ∈ E2} is a basis of RT₀(T);

4. divϕ_E =

{±_|^|_T^E_±^|_|, on T_±, 0, elsewhere.

Proof. For a proof see [2].

B

P

₊

A

P

₋

E n

_E

T

₊

T

₋

Figure 4.5: Orientation ofn_E.

The assembling of the matrices in K^h associated to the Darcy region and the Matlab implementation of right-hand sides and Neumann conditions is done following section 4 and 6 of [2]. Note that before we solve the algebraic problem we project the

CHAPTER 4. DISCRETIZATION OF THE WEAK FORMULATION 52 right hand side into the orthogonal space of the kernel. Since we approximated the appearing integrals on the right hand side with quadrature formulas, the computation was not exact. We use projection to guarantee orthogonality to the kernel, which implies solvability of the algebraic problem. Let f be the right hand side of the algebraic formulation of the coupled Stokes-Darcy problem and let x₀ be an element of the kernel of the stiness matrix K^h, i.e. x₀ ∈ ker{K^h}, then the projection f˜:=f+αx₀ shall fulll

f +αx₀ ⊥x₀ (4.93)

which is equivalent to

f ·x₀ +αx₀ ·x₀ = 0. (4.94)

Therefore we get

α=−f·x₀

x₀·x₀. (4.95)

Note that for the vector x₀ ∈ ker{K^h} we choose a vector consisting of one and zero entries. For the velocity nodes the corresponding entries are zero and for pressure nodes the corresponding entries are equal 1.

After the assembling of the stiness matrix K^h and the right hand side we x one degree of freedom for the pressure. This means that we delete the corresponding row and column ofK^h and the corresponding entry in the right hand side. The so obtained system of equations is solved with the backslash operator in Matlab. Note that Matlab uses a LU factorization for that. Remember that as a last step we have to subtract the mean value ofp^h from the computed solutionp^h to get the seeked solution in M^h.

Convergence rate:

First we indicate that u = u_g+ω, where u_g ∈ V^g and ω ∈ V is the solution of the homogenized problem.

Before we discuss the convergence of the discrete solution to the Stokes-Darcy problem using Taylor-Hood and Raviart-Thomas elements we collect some denitions of interpolation operators and estimates from [9] and the references therein.

First let I_V^h :W →V^h,0 be the interpolation operator where W is the subspace of suciently smooth functions

W := {v = (v₁, v₂)∈X :v_i ∈W_i :=X_i∩(H^si(Ω_i))², i= 1,2 and v1·n2|ΓI =v2·n2|ΓI inL²(ΓI)},

wherebys_i is chosen suciently large. We deneI_V^h := (I_X^h

1v₁, I_X^h

2v₂−δ₂^h) :W →V^h,0, whereI_X^h_i :Wi →X_i^h,0, i= 1,2are nite element interpolation operators andδ₂^h is the correction such that the continuity of the normal velocities is guaranteed in a discrete sense across the interface Γ_I. We dene δ₂^h := I_X^h₂v₁ −I_X^h₂(I_X^h₁v₁). For more detailed information of the construction of the correction δ₂^h see [9].

In our case the interpolation operator I_V^h can be interpreted in the following way.

We constructI_X^h₁ω₁ inX₁^h,0as a linear combination of the basis functions inX₁^h,0where

we use the function values of the exact solutionω₁ in the corresponding interior nodes as coecients. The interpolated functionI_X^h₂ω₂−δ^h₂ is a linear combination of the edge basis functions in X₂^h,0. The corresponding coecients are the normal components of the exact solutionω₂in the midpoints of the edges except for the edges on the interface.

There we use equation (4.63) to compute the normal components of ω2.

Due to Proposition 4.2 in the paper of Layton, Schieweck an Yotov [9] we know that for all ω∈W ⊂V

∥ω−I_V^hω∥=O(hT). (4.96) Remember that I_M^h

1 :M₁∩H²(Ω₁)→M₁^h is the interpolation operator such that

∥q₁−I_M^h

1q₁∥0,T ≤Ch²_T|q₁|2,T (4.97) and I_M^h

2 :M₂ →M₂^h is the L²-orthogonal projection such that for all q₂ ∈M₂ (I_M^h

2q₂−q₂, q^h₂)_0,Ω2 = 0 for all q₂^h ∈M₂^h and

∥q2−I_M^h₂q2∥0,T ≤ChT|q2|1,T. (4.98) Next, let I_M^h := (I_M^h

1p₁, I_M^h

2p₂) :M →M₁^h×M₂^h be the interpolation operator for the pressure functions inM.

From (4.97) and (4.98) we get that ∥p1 − I_M^h₁p1∥M1 = O(h²_T) and

∥p₂−I_M^h

2p₂∥M2 =O(h_T), respectively. This implies

∥p−I_M^h p∥M = (∥p₁−I_M^h

1p₁∥²M1 +∥p₂−I_M^h

2p₂∥²M2)^1/2 =O(h_T). (4.99) Since X₁^h and M₁^h contain polynomials of degree r₁ = 2 and r₁ −1 = 1 and X₂^h andM₂^h contain polynomials of degreer₂ = 0 and l₂ = 0, respectively we know due to Theorem 4.3 that

∥ω−ω^h∥X +∥p−p^h∥M =O(h^r_T¹ +h^r_T²⁺¹+h^l_T²⁺¹) = O(h_T), (4.100) which means that the convergence rate is equal one. If we halve the mesh width h then also the discretization error halves itself.

Note that all statements and estimates which hold for the homogeneous problem can also be shown for the inhomogeneous case.

Chapter 5

Numerical results

In this chapter we want to verify the error estimate in Theorem 4.3, which is based on Layton, Schieweck and Yotov [9], with a numerical experiment.

We consider the following concrete problem: Find u₁, u₂, p₁ and p₂ such that the equations

−divT(u₁, p₁) = f₁ in Ω₁ (balance of momentum), (5.1) divu1 = 0 in Ω1 (conservation of mass), (5.2) u₁ =g₁ on Γ₁ (boundary condition), (5.3) u₂ =−K∇p₂ in Ω₂ (Darcy's law), (5.4) divu₂ =f₂ in Ω₂ (balance of mass), (5.5) u₂·n₂ =g₂ on Γ₂ (boundary condition), (5.6) u₁·n₁+u₂ ·n₂ = 0 on Γ_I (conservation of mass), (5.7) p₁−2µD(u₁)n₁·n₁ =p₂ on Γ_I (balance of normal forces), (5.8) u₁·τ₁ =

√k˜₁

µα (−⃗t(u₁, p₁))·τ₁, on Γ_I (Beavers-Joseph-Saman cond.), (5.9) are fullled for the uid region Ω₁ = (0,1) × (0,1), the porous medium region Ω2 = (0,1)×(−1,0), the interface ΓI = (0,1)× {0} and for the functions g1, g2, f1

and f₂, which are given by g₁ :=

( −cos(πx) sin(πy) sin(πx) cos(πy)

)

, g₂ :=





−sin(πy) on {0} ×(0,−1) sin(πx) on (0,1)× {−1}

−sin(πy) on {1} ×(0,−1) ,

f₁ :=

( −(4π²+ 1) cos(πx) sin(πy)

−sin(πx) cos(πy)

)

, f₂ :=−2πsin(πx) sin(πy), Furthermore we setα =µ= 1 and K =I. Note thatτ₁ =

( 1 0

) . The exact solution to this problem is

u₁ =

( −cos(πx) sin(πy) sin(πx) cos(πy)

)

, p₁ =−sin(πx) sin(πy)(1

π + 2π) + 4 π, 54

Mesh width h ₁₀¹ ₂₀¹ ₄₀¹ ₈₀¹ Error: ∥I_V^hu−u^h∥X +∥I_M^hp−p^h∥M 0.0854 0.021 0.0052 0.0013 Table 5.1: The error∥I_V^hu−u^h∥X +∥I_M^h p−p^h∥M computed for dierent meshes.

u₂ =

( cos(πx) sin(πy) sin(πx) cos(πy)

)

, p₂ =−1

πsin(πx) sin(πy) + 4 π. The approximated solutions are depicted in Figure 5.1 and Figure 5.2.

Remember that we know due to Theorem 4.3 that

∥u−u^h∥X +∥p−p^h∥M =O(h_T). (5.10) Now we want to illustrate that our numerical results are in good agreement with this error estimation. Note that we computed the error between the interpolated solutionsu, pand the approximated solutions u^h, p^h since it is easier to compute than

∥u−u^h∥X +∥p−p^h∥M. Let h be the short side of the triangles in our uniform grid.

In Table 5.1 we can observe that if we halve the mesh width h then the error between the interpolated solutions u, p and the approximated solutions u^h, p^h is divided by a factor around 4. This means that ∥I_V^hu−u^h∥+∥I_M^h p−p^h∥= O(h²). However, from the triangular inequality and from (4.96) and (4.99) we get that the disretization error satises

∥u−u^h∥X+∥p−p^h∥M ≤ ∥|u−{zI_V^hu∥X}

=O(h)

+∥|p−{zI_M^h p∥M}

=O(h)

+∥|I_V^hu−u^h∥X +{z∥I_M^h p−p^h∥M}

=O(h²)

.

This means that the discretization error for our model problem

∥u−u^h∥X +∥p−p^h∥M is equal toO(h)which coincides with the error estimate (5.10).

CHAPTER 5. NUMERICAL RESULTS 56

Figure 5.1: Computed solution for the velocity u plotted in the midpoints of the elements.

Figure 5.2: Computed solution for the pressure p in the Stokes and the Darcy region.

Chapter 6 Conclusion

In this thesis we considered a Stokes ow and a porous media ow coupled across an interface through appropriate interface conditions. We gured out that the in-terface condition by Beavers and Joseph does not t in our problem setting but the Beavers-Joseph-Saman condition leads us to a well posed problem. We derived a weak formulation for our model problem and also other possible formulations were discussed. We showed boundedness, ellipticity and the inf-sup condition for the re-spective bilinear forms and we concluded from Brezzi's theorem the well-posedness of the Stokes-Darcy problem.

Then we treated the discrete problem for which we also veried existence and uniqueness of a solution. Note that while the verication of the boundedness and ellipticity was rather simple, the proof of the discrete inf-sup condition was quite elaborate. In the next step we exhibited a possible nite element discretization, where we used Taylor-Hood elements in the uid region and Raviart-Thomas elements in the porous medium region. We took a closer look at the characterization of the discrete velocity and how the degrees of freedom are linked on the interface. Then we showed the resulting matrix form of the problem and gave short computational explanations and hints.

In the end we presented numerical results and conrmed that the behaviour of the error between the exact and approximated solution is in good agreement with the error estimation obtained from theory.

58

Bibliography

[1] R. A. Adams. Sobolev Spaces. Academic Press, 1975.

[2] C. Bahriawati and C. Carstensen. Three matlab implementations of the lowest-order Raviart-Thomas MFEM with a posteriori error control. Computational methods in applied mathematics, Vol. 5, No. 4, pp. 333-361, 2005.

[3] D. Braess. Finite Elemente: Theorie, schnelle Löser und Anwendungen in der Elastizitätstheorie. Springer-Verlag Berlin Heidelberg, 2007.

[4] F. Brezzi and M. Fortin. Mixed and hybrid nite element methods. Springer-Verlag, New York, 1991.

[5] M. Discacciati and A. Quarteroni. Navier-Stokes/Darcy coupling: modeling, anal-ysis, and numerical approximation. Rev. Mat. Complut. 22, No. 2, 315-426, 2009.

[6] J. Galvis and M. Sarkis. Non-matching mortar discretization analysis for the coupling Stokes-Darcy equations. Electronic Transactions on Numerical Analysis, Vol. 26, pp. 350-384, 2007.

[7] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations.

Theory and algorithms. Springer-Verlag, Berlin, 1986.

[8] P. C. Hammer and A. H. Stroud. Numerical Integration over Simplexes. Math.

Tables and other Aids to Computation, Vol. 10, No 55, pp. 137-139, 1956.

[9] W. J. Layton, F. Schiewek and I. Yotov. Coupling Fluid Flow with Porous Media Flow. SIAM Journal on Numerical Analysis, Vol. 40, No. 6, pp. 2195-2218, 2003.

[10] T. P. Mathew. Domain Decomposition and Iterative Renement Methods for Mixed Finite Element Discretizations of Elliptic Problems. Ph.D. thesis, Courant Institute of Mathematical Sciences, New York University, New York, 1989.

[11] J. Stoer. Numerische Mathematik 1. Springer-Verlag, Berlin, 2005 (9. Auage).

59

Im Dokument A finite element method for solving the Stokes-Darcy problem / eingereicht von: Manuela Redl (Seite 53-65)