A New Finite Element Method for Darcy-Stokes-Brinkman Equations

Abstract

We present a new finite element method for Darcy-Stokes-Brinkman equations using primal and dual meshes for the velocity and the pressure, respectively. Using an orthogonal basis for the discrete space for the pressure, we use an efficiently computable stabilization to obtain a uniform convergence of the finite element approximation for both limiting cases.

1. Introduction

Recently developing an efficient finite element method for Darcy-Stokes-Brinkman problem has become an active area of research [1–5]. This mathematical model is often called a generalized Stokes problem [6] and is also a special case of Oseen's equations when the convective term is absent [1]. This model incorporates both Darcy and Stokes problem. If the viscosity vanishes, the model becomes the Darcy problem, and if another parameter is set to zero, it reduces to Stokes problem. Therefore, this is important when studying the coupling between Stokes and Darcy equations. Some applications of Darcy-Stokes-Brinkman problem are given in [4]. This model also arises from the time discretization of time-dependent Stokes equations.

In this paper we propose a new finite element method for this problem using stabilization based on a local projection. Our technique is the combination of the local projection stabilization for Stokes equations [7] and the grad-div stabilization proposed earlier for the Stokes equations [8]. We note that we mainly focus on the Darcy-Stokes-Brinkman model so that the divergence of the velocity field is not identically zero but is given as a function in 𝐿²-space as in [3, 4]. Using a standard continuous linear, bilinear, or trilinear finite element space for the velocity, the grad-div stabilization [8] is to be modified in this case. Since the modification involves projecting the divergence of the velocity field onto the finite element space for the pressure, we first use a discontinuous pressure space in the finite element scheme. This leads to an efficient realization of the stabilization. The finite element space for the pressure is the piecewise constant function space on the dual mesh. The inf-sup condition is satisfied using an element-wise defined bubble function in each element.

This paper is organized as follows. We present the boundary value problem for the Darcy-Stokes-Brinkman problem in the next section. We present our finite element method, propose stabilization, and analyze the discrete problem in Section 3. We prove an optimal a priori error estimate for both cases in this section. Finally, a conclusion is drawn in the last section.

2. Darcy-Stokes-Brinkman Problem

This section is devoted to the introduction of the boundary value problem of the Darcy-Stokes-Brinkman flow. Let Ω in ℝ^𝑑, 𝑑∈{2,3}, be a bounded domain with polygonal or polyhedral boundary Γ. For a prescribed body force 𝐟∈[𝐿²(Ω)]^𝑑 and another function 𝑔∈𝐿²(Ω), the Darcy-Stokes-Brinkman equations with homogeneous Dirichlet boundary condition in Γ read −𝜈Δ𝐮+∇𝑝+𝜎𝐮=𝐟inΩ,div𝐮=𝑔inΩ, (1) with 𝐮=0 on Γ, where 𝐮 is the velocity, 𝑝 is the pressure, and 𝜈 denotes the viscosity of the fluid. For simplicity we assume that 𝜈≥0 and 𝜎 are constant. However, the theory also can be easily extended to cover the case when 𝜈 and 𝜎 depend on the spatial coordinates 𝐱∈Ω. When the viscosity of the fluid 𝜈 tends to zero, this problem reduces to the Darcy problem, and when the parameter 𝜎 approaches zero, this problem reduces to the Stokes problem. Here we aim at developing a finite element method which is stable in both limit cases.

Here we use standard notations 𝐿²(Ω), 𝐻¹(Ω), and 𝐻¹₀(Ω) for Sobolev spaces; see [9, 10] for details. Let 𝐖∶=[𝐻¹₀(Ω)]^𝑑 be the vector Sobolev space with inner product (⋅,⋅)_1,Ω and norm ‖⋅‖_1,Ω defined in the standard way: (𝐮,𝐯)_1,Ω∶=∑^𝑑_𝑖=1(𝑢_𝑖,𝑣_𝑖)_1,Ω, with the norm being induced by this inner product. We also define another subspace 𝑀 of 𝐿²(Ω) as 𝑀=𝑞∈𝐿²(Ω)∶_Ω𝑞𝑑𝑥=0. (2)

The weak formulation of the Darcy-Stokes-Brinkman problem is to find (𝐮,𝑝)∈𝐖×𝑀 such that 𝜈_Ω∇𝐮∶∇𝐯𝑑𝑥+_Ωdiv𝑣𝑝𝑑𝑥+𝜎_Ω𝐮𝐯𝑑𝑥=ℓ(𝐯),𝐯∈𝐖,_Ωdiv𝐮𝑞𝑑𝑥=𝑔𝑞,𝑞∈𝑀, (3) where ℓ(𝐯)=∫_Ω𝐟⋅𝐯𝑑𝑥. It is well known that the weak formulation of this Darcy-Stokes-Brinkman problem is well-posed [3, 11].

We note that, when the viscosity 𝜈→0, the appropriate space to seek a solution of the variational problem (8) is 𝐻₀(div,Ω) defined as 𝐻₀(div,Ω)={𝝈∈𝐻(div,Ω)∶𝝈⋅𝐧=0on𝜕Ω}, (4) where 𝐻(div,Ω)=𝐮∈𝐿²(Ω)^𝑑∶div𝐮∈𝐿²(Ω) (5) (see [3]). A standard norm on the space 𝐻(div,Ω) is ‖𝐮‖²_div,Ω=‖𝐮‖²_0,Ω+‖div𝐮‖²_0,Ω. (6)

With this change of the character of the solution, when the viscosity 𝜈 vanishes, we define an appropriate norm for the problem (8) as ‖𝐯‖²_𝑎∶=𝜈‖𝐯‖²_1,Ω+𝜎‖𝐯‖²_div,Ω,𝐯∈[𝜈𝐖∩𝜎𝐻(div,Ω)]. (7)

Therefore now we use the space 𝐕∶=𝜈𝐖∩𝜎𝐻(div,Ω) as the space for the velocity 𝐮. Due to the occurrences of the parameters in the problem it is appropriate to use the suitable norm defined above although it is clear that 𝐻(div,Ω)⊂[𝐻¹(Ω)]^𝑑. Thus we consider the problem of finding (𝐮,𝑝)∈𝐕×𝑀 such that 𝜈_Ω∇𝐮∶∇𝐯𝑑𝑥+_Ωdiv𝑣𝑝𝑑𝑥+𝜎_Ω𝐮𝐯𝑑𝑥=ℓ(𝐯),𝐯∈𝐖,_Ωdiv𝐮𝑞𝑑𝑥=𝑔𝑞,𝑞∈𝑀. (8)

3. Finite Element Discretizations

We consider a quasiuniform triangulation 𝒯_ℎ of the polygonal or polyhedral domain Ω, where 𝒯_ℎ consists of simplices (triangles or tetrahedra), where ℎ denotes the mesh-size. Note that 𝒯_ℎ denotes the set of elements.

For nonnegative integer 𝑛 and an element 𝐾∈𝒯_ℎ, we let 𝒫_𝑛(𝐾) denote the space of polynomials in two variables of total degree less than or equal to 𝑛, where 𝐾 is a triangle or tetrahedron. Let 𝐾 be a reference triangle or tetrahedron, where the reference triangle is defined as 𝐾∶=𝑥,𝑦∈ℝ²∶𝑥>0,𝑦>0,𝑥+𝑦<1, (9) and the reference tetrahedron is 𝐾∶=𝑥,𝑦,𝑧∈ℝ³∶𝑥>0,𝑦>0,𝑧>0,𝑥+𝑦+𝑧<1. (10)

A typical element 𝐾∈𝒯_ℎ is generated by an affine map 𝐹_𝐾 from the reference element 𝐾, in which 𝐹_𝐾 is defined using the basis functions corresponding to 𝒫₁.

The finite element space of displacement is taken to be the space of continuous functions whose restrictions to an element 𝐾 are obtained by maps of linear functions from the reference element: 𝑆_ℎ∶=𝑣_ℎ∈𝐻¹₀(Ω),𝑣_ℎ∣𝐾=𝑣_ℎ̂∘𝐹_𝐾⁻¹,𝐯_ℎ∈𝒫₁𝐾,𝐾∈𝒯_ℎ. (11)

Let the set of all vertices in 𝒯_ℎ be denoted by 𝒱_ℎ∶={𝐱_𝑖}^𝑁_𝑖=1 and 𝒩_ℎ∶={1,2,…,𝑁}.

A dual mesh 𝒯^∗_ℎ is introduced based on the primal mesh 𝒯_ℎ so that the elements of 𝒯^∗_ℎ are called control volumes. Each control volume element 𝑉_𝑖∈𝒯^∗_ℎ is associated with a vertex 𝐱_𝑖∈𝒱_ℎ. For simplicity we explain the construction of the dual mesh for triangular meshes [12, 13]. The construction is easily extended to a tetrahedral mesh. For a vertex 𝐱_𝑖∈𝒱_ℎ let 𝒯^𝑖_ℎ and ℰ^𝑖_ℎ be the set of elements and edges touching 𝐱_𝑖, respectively. Then the volume element 𝑉_𝑖 corresponding to the vertex 𝐱_𝑖∈𝒱_ℎ is the polygonal domain joining centroids of all elements in 𝒯^𝑖_ℎ and centroids of all edges in ℰ^𝑖_ℎ. The set of all volume elements in 𝒯^∗_ℎ will form a nonoverlapping decomposition of the polygonal domain Ω; see also [12, 13]. A dual mesh for a triangular mesh is shown in Figure 1.

Figure 1 

Primal and dual meshes for a triangular grid.

In the following we use a generic constant 𝐶, which takes different values at different places but is always independent of the mesh-size ℎ and the parameters 𝜈 and 𝜎.

We call the control volume mesh 𝒯^∗_ℎ regular or quasiuniform if there exists a positive constant 𝐶>0 such that 𝐶ℎ^𝑑≤𝑉_𝑖≤ℎ^𝑑,𝑉_𝑖∈𝒯^∗_ℎ, (12) where ℎ is the maximum diameter of all elements 𝑇∈𝒯_ℎ. It can be shown that, if 𝒯_ℎ is locally regular, that is, there is a constant 𝐶 such that 𝐶ℎ^𝑑_𝑇≤|𝑇|≤ℎ^𝑑_𝑇,𝑇∈𝒯_ℎ (13) with diam(𝑇)=ℎ_𝑇 for all elements 𝑇∈𝒯_ℎ, then this dual mesh 𝒯^∗_ℎ is also locally regular. The dual volume element space 𝑆^∗_ℎ to discretize the pressure is now defined by 𝑆^∗_ℎ∶=𝑝∈𝐿²₀(Ω)∶𝑝_∣𝑉∈𝒫₀(𝑉),𝑉∈𝒯^∗_ℎ. (14)

Now 𝑝_ℎ∈𝑆^∗_ℎ and 𝑢_ℎ∈𝑆_ℎ can be written as 𝑢_ℎ=∑^𝑁_𝑖=1𝑢_𝑖𝜙_𝑖 and 𝑝_ℎ=∑^𝑁_𝑖=1𝑝_𝑖𝜒_𝑖, where 𝜙_𝑖 are the standard nodal basis functions associated with the vertex 𝐱_𝑖, and 𝜒_𝑖 are the characteristic functions of the volume 𝑉_𝑖.

For the reference element 𝐾, we define a bubble function𝑏𝐾=27𝑥𝑦1−𝑥−𝑦,if𝐾isthereferencetriangle,256𝑥𝑦𝑧1−𝑥−𝑦−𝑧,if𝐾isthereferencetetrahedron. (15)

Then defining the space of bubble functions 𝐵_ℎ as 𝐵_ℎ∶=𝑏_ℎ∈𝐻¹₀(Ω),𝑏_ℎ∣_𝑘𝐾=𝑏∘𝐹_𝐾⁻¹𝐾,𝑏isabubblefunctionon𝐾,𝐾∈𝒯_ℎ, (16) we introduce our finite element space for velocity as 𝐕_ℎ=[𝑆_ℎ⊕𝐵_ℎ]^𝑑.

With this definition of the spaces 𝐕_ℎ and 𝑆^∗_ℎ we have the following lemma [13].

Lemma 1. There exists a constant 𝛽>0 independent of the mesh-size ℎ such that inf_{𝑞ℎ∈𝑆∗ℎ⧵{0}}sup_{𝐯ℎ∈𝐕ℎ⧵{0}}∫_Ωdiv𝐯_ℎ𝑞_ℎ𝐯_ℎ1,Ω𝑞_ℎ0,Ω≥𝛽. (17)

We now consider the limiting case when 𝜈→0. The main difficulty in this case using the conforming approximation is to get the coercivity of the bilinear form 𝑎𝐮_ℎ,𝐯_ℎ=𝜎_Ω𝐮_ℎ⋅𝐯_ℎ𝑑𝑥 (18) on the kernel space 𝐊_ℎ∶=𝐯_ℎ∈𝐕_ℎ∶_Ωdiv𝐯_ℎ𝑞_ℎ𝑑𝑥=0,𝑞_ℎ∈𝑆^∗_ℎ. (19)

When 𝜈→0, the coercivity of the bilinear form 𝑎(⋅,⋅) on 𝐊_ℎ demands the existence of a constant 𝐶>0 such that 𝑎𝐯_ℎ,𝐯_ℎ≥𝐶𝐯_ℎ²_𝑎=𝐶𝜎𝐯_ℎ²_div,Ω. (20)

Assuming that 𝜎>0 this condition requires the existence of a constant 𝐶>0 such that 𝐯_ℎ0,Ω≤𝐶div𝐯_ℎ0,Ω. (21)

Since this condition is not satisfied for the space 𝐕_ℎ, we add a stabilising term to the bilinear form 𝑎(⋅,⋅) to obtain another bilinear form 𝑎(⋅,⋅)𝑎𝐮_ℎ,𝐯_ℎ∶=𝑎𝐮_ℎ,𝐯_ℎ+𝜎_Ωdiv𝐮_ℎ−𝑄_ℎdiv𝐮_ℎ×div𝐯_ℎ−𝑄_ℎdiv𝐯_ℎ𝑑𝑥, (22) where 𝑄_ℎ is an orthogonal projection onto the space 𝑆^∗_ℎ. With this definition of the bilinear form, our discrete problem is to find [𝐮_ℎ,𝑝_ℎ]∈𝐕_ℎ×𝑆^∗_ℎ such that 𝑎𝐮_ℎ,𝐯_ℎ+𝑏𝐯_ℎ,𝑝_ℎ=ℓ𝐯_ℎ,𝐯_ℎ∈𝐕_ℎ,𝑏𝐮_ℎ,𝑞_ℎ=𝑔𝑞_ℎ,𝑞_ℎ∈𝑆^∗_ℎ. (23)

Remark 2. Due to the fact that the basis functions for the space 𝑆^∗_ℎ are orthogonal with respect to 𝐿²-inner product, computing the action of the operator 𝑄_ℎ is very easy and efficient. Since the basis functions of 𝑆^∗_ℎ are piecewise constant, the stabilization is also local. In fact, similar stabilization can be applied using mini type finite element method [14] at least for a simplicial grid, where computation of the operator 𝑄_ℎ will be more difficult, expensive, and nonlocal.

Remark 3. If the source term 𝑔 is identically zero, the stabilization is simplified as 𝑄_ℎdiv𝐯_ℎ=0. In fact, such stabilization is called grad-div stabilization and was used before for Stokes equations in [8].

We now show that the bilinear form 𝑎(⋅,⋅) is uniformly coercive on 𝐊_ℎ with respect to the norm ‖⋅‖_𝑎.

Lemma 4. One has 𝑎𝐯_ℎ,𝐯_ℎ=𝜈𝐯_ℎ²_1,Ω+𝜎𝐯_ℎ²_div,Ω,𝐯_ℎ∈𝐊_ℎ. (24)

Proof. Note that 𝑎𝐯_ℎ,𝐯_ℎ=𝜈𝐯_ℎ²_1,Ω+𝜎𝐯_ℎ²_0,Ω+𝜎div𝐯_ℎ−𝑄_ℎdiv𝐯_ℎ²_0,Ω=𝜈𝐯_ℎ²_1,Ω+𝜎𝐯_ℎ²_0,Ω+𝜎div𝐯_ℎ²_0,Ω=𝜈𝐯_ℎ²_1,Ω+𝜎𝐯_ℎ²_div,Ω, (25) as 𝑄_ℎdiv𝐯_ℎ=0 for 𝐯_ℎ∈𝐾_ℎ.

Clearly, we see that the bilinear forms 𝑎(⋅,⋅) and 𝑏(⋅,⋅), and the linear form ℓ(⋅) are also continuous on the spaces where they are defined. Moreover, since for a constant 𝐶>0 we have 𝐯_ℎ𝑎≤𝐶𝐯_ℎ1,Ω, (26) we have the inf-sup condition from Lemma 1 when the ‖⋅‖_1,Ω is replaced by ‖⋅‖_𝑎. Thus we have the following theorem.

Theorem 5. The discrete saddle point problem (23) has a unique solution [𝐮_ℎ,𝑝_ℎ]∈𝐕_ℎ×𝑆^∗_ℎ such that 𝐮_ℎ𝑎+𝑝_ℎ0,Ω≤𝐶‖ℓ‖_0,Ω. (27)

The error estimate is given in the following theorem.

Theorem 6. Let [𝐮,𝑝]∈𝐕×𝑀 be the unique solution of (8) and [𝐮_ℎ,𝑝_ℎ]∈𝐕_ℎ×𝑆^∗_ℎ that of the discrete (23). One has the following error estimate uniformly with respect to 𝜎 and 𝜈: 𝐮−𝐮_ℎ𝑎+𝑝−𝑝_ℎ0,Ω≤𝐶inf_{𝐯ℎ∈𝐕ℎ}𝐮−𝐯_ℎ𝑎+inf_{𝑞ℎ∈𝑆∗ℎ}𝑝−𝑞_ℎ0,Ω. (28)

4. Conclusion

We have presented a very efficient finite element technique for approximating the solution of Darcy-Stokes-Brinkman equations. Using local projection stabilization for the limiting case of vanishing viscosity we have shown that the error estimate is uniform for both limiting cases.