A Discontinuous Finite Volume Method for the Darcy-Stokes Equations

Abstract

This paper proposes a discontinuous finite volume method for the Darcy-Stokes equations. An optimal error estimate for the approximation of velocity is obtained in a mesh-dependent norm. First-order 𝐿²-error estimates are derived for the approximations of both velocity and pressure. Some numerical examples verifying the theoretical predictions are presented.

1. Introduction

The study of discontinuous Galerkin methods has been a very active research field since they were proposed by Reed and Hill [1] in 1973. Discontinuous Galerkin methods use discontinuous functions as finite element approximation and enforce the connections of the approximate solutions between elements by adding some penalty terms. The flexibility of discontinuous functions gives discontinuous Galerkin methods many advantages, such as high parallelizability and localizability. Arnold et al. [2] provided a framework for the analysis of a large class of discontinuous Galerkin methods for second-order elliptic problems.

Based on the advantages of using discontinuous functions for approximation in discontinuous Galerkin methods, it is natural to consider using discontinuous functions as trial functions in the finite volume method, which is called the discontinuous finite volume method. Such a method has the flexibility of the discontinuous Galerkin method and the simplicity and conservative properties of the finite volume method. Ye [3] developed a new discontinuous finite volume method and analyzed it for the second-order elliptic problem. Bi and Geng [4] proposed the semidiscrete and the backward Euler fully discrete discontinuous finite volume element methods for the second-order parabolic problems. Ye [5] considered the discontinuous finite volume method for solving the Stokes problems on both triangular and rectangular meshes and derived an optimal order error estimate for the approximation of velocity in a mesh-dependent norm and first-order 𝐿²-error estimates for the approximations of both velocity and pressure.

The Darcy-Stokes problem is interesting for a variety of reasons. Apart from being a modeling tool in its own right, it also appears, less obviously, in time-stepping methods for Stokes and for high Reynolds number flows (where of course the convective term causes additional difficulties). In [6], the nonconforming Crouzeix-Raviart element is stabilized for the Darcy-Stokes problem with terms motivated by a discontinuous Galerkin approach. In [7], a new stabilized mixed finite element method is presented for the Darcy-Stokes equations.

In this paper, we will extend the discontinuous finite volume methods to solve the Darcy-Stokes equations. In our methods, velocity is approximated by discontinuous piecewise linear functions on triangular meshes and by discontinuous piecewise rotated bilinear functions on rectangular meshes. Piecewise constant functions are used as the test functions for velocity in the discontinuous finite volume methods. We obtained an optimal error estimate for the approximation of velocity in a mesh-dependent norm. First-order 𝐿²-error estimates are derived for the approximations of both velocity and pressure. For the sake of simplicity and easy presentation of the main ideas of our method, we restrict ourselves to the model problem.

We consider the Darcy-Stokes equations𝜎𝐮−𝜇Δ𝐮+∇𝑝=𝐟,𝑥∈Ω,∇⋅𝐮=0,𝑥∈Ω,𝐮=𝟎,𝑥∈𝜕Ω, (1.1a)(1.1b)(1.1c)where Ω is a bounded polygonal domain in 𝑅² with boundary 𝜕Ω. 𝐮=(𝑢₁,𝑢₂) is the velocity, 𝑝 is the pressure, and 𝐟 is a given force term. We assume 𝜎=1,𝜇=1.

2. Discontinuous Finite Volume Formulation

Let ℛ_ℎ be a triangular or rectangular partition of Ω. The triangles or rectangles in ℛ_ℎ are divided into three or four subtriangles by connecting the barycenter of the triangle or the center of the rectangles to their corner nodes, respectively. Then we define the dual partition 𝒯_ℎ of the primal partition ℛ_ℎ to be the union of the triangles shown in Figures 1 and 2 for both triangular and rectangular meshes.

Figure 1 

Element 𝑇∈𝒯ℎ for triangular mesh.

Figure 2 

Element 𝑇∈𝒯ℎ for rectangular mesh.

Let 𝑃_𝑘(𝑇) consist of all the polynomials with degree less than or equal to 𝑘 defined on 𝑇. We define the finite dimensional trial function space for velocity on a triangular partition by 𝑉_ℎ=𝐯∈𝐿²(Ω)²∶𝐯|_𝐾∈𝑃₁(𝐾)²,∀𝐾∈ℛ_ℎ (2.1) and on rectangular partition by 𝑉_ℎ=𝐯∈𝐿²(Ω)²∶𝐯|_𝐾∈𝑄₁(𝐾)²,∀𝐾∈ℛ_ℎ, (2.2) where 𝑄₁ denotes the space of functions of the form 𝑎+𝑏𝑥₁+𝑐𝑥₂+𝑑(𝑥²₁−𝑥²₂) on 𝐾.

Let 𝑄_ℎ be the finite dimensional space for pressure 𝑄_ℎ=𝑞∈𝐿²₀(Ω)∶𝑞_𝐾∈𝑃₀(𝐾),∀𝐾∈ℛ_ℎ, (2.3) where 𝐿²₀(Ω)=𝑞∈𝐿²(Ω)∶_Ω𝑞𝑑𝑥=0. (2.4) Define the finite dimensional test function space 𝑊_ℎ for velocity associated with the dual partition 𝒯_ℎ as 𝑊_ℎ=𝝃∈𝐿²(Ω)²∶𝝃|_𝑇∈𝑃₀(𝑇)²,∀𝑇∈𝒯_ℎ. (2.5)

Multiplying (1.1a) and (1.1b) by 𝜉∈𝑊_ℎ and 𝑞∈𝑄_ℎ, respectively, we have (𝐮,𝝃)−_{𝑇∈𝒯ℎ𝜕𝑇}𝜕𝐮𝜕𝐧⋅𝝃𝑑𝑠+_{𝑇∈𝒯ℎ𝜕𝑇}𝑝𝝃⋅𝐧𝑑𝑠=(𝐟,𝝃_{𝐾∈ℛℎ𝐾}),∇⋅𝐮𝑞𝑑𝑥=0, (2.6) where 𝐧 is the unit outward normal vector on 𝜕𝑇.

Let 𝑇_𝑗∈𝒯_ℎ(𝑗=1,…,𝑡) be the triangles in 𝐾∈ℛ_ℎ, where 𝑡=3 for triangular meshes and 𝑡=4 for rectangular meshes, as shown as Figures 3 and 4. Then we have _{𝑇∈𝒯ℎ𝜕𝑇}𝜕𝐮𝜕𝐧⋅𝝃𝑑𝑠=_{𝐾∈ℛℎ}^𝑡_{𝑗=1𝐴𝑗+1𝐶𝐴𝑗}𝜕𝐮𝜕𝐧⋅𝝃𝑑𝑠+_{𝐾∈ℛℎ𝜕𝐾}𝜕𝐮𝜕𝐧⋅𝝃𝑑𝑠, (2.7) where 𝐴_𝑡+1=𝐴₁.

Figure 3 

Triangular partition and its dual.

Figure 4 

Rectangular partition and its dual.

For vectors 𝐯=(𝑣₁,𝑣₂) and 𝐧=(𝑛₁,𝑛₂), let 𝐯⊗𝐧 denote the matrix whose 𝑖𝑗th component is 𝑣_𝑖⋅𝑛_𝑗 as in [5]. For two matrix valued variables 𝜎 and 𝜏, we define 𝜎∶𝜏=∑²_𝑖,𝑗=1𝜎_𝑖,𝑗𝜏_𝑖,𝑗. Let Γ=∑_{𝐾∈ℛℎ}𝜕𝐾,Γ₀=Γ⧵𝜕Ω. Let 𝑒 be an interior edge shared by two elements 𝐾₁ and 𝐾₂ in ℛ_ℎ. We define the average {⋅} and jump [⋅] on 𝑒 for scalar 𝑞, vector 𝑤, and matrix 𝜏, respectively. If 𝑒∈Γ₀, 𝑞=12𝑞|_𝜕𝐾1+𝑞|_𝜕𝐾2,{𝐰}=12𝐰|_𝜕𝐾1+𝐰|_𝜕𝐾2,{𝜏}=12𝜏|_𝜕𝐾1+𝜏|_𝜕𝐾2,𝑞=𝑞|_𝜕𝐾1𝐧₁+𝑞|_𝜕𝐾2𝐧₂,[𝐰]=𝐰|_𝜕𝐾1⋅𝐧₁+𝐰|_𝜕𝐾2⋅𝐧₂,[𝜏]=𝜏|_𝜕𝐾1⋅𝐧₁+𝜏|_𝜕𝐾2⋅𝐧₂, (2.8) where 𝐧₁ and 𝐧₂ are unit normal vectors on e pointing exterior to 𝐾₁ and 𝐾₂, respectively. We also define a matrix valued jump ⟦⋅⟧ for a vector 𝐰 as ⟦𝐰⟧=𝐰|_𝜕𝐾1⊗𝐧₁+𝐰|_𝜕𝐾2⊗𝐧₂. (2.9) If 𝑒∈𝜕Ω, define 𝑞=𝑞,[𝐰]=𝐰⋅𝐧,{𝜏}=𝜏,⟦𝐰⟧=𝐰⊗𝐧. (2.10) A straightforward computation gives _{𝐾∈ℛℎ𝜕𝐾}𝑞𝐯⋅𝐧𝑑𝑠=_{𝑒∈Γ0𝑒}𝑞⋅{𝐯}𝑑𝑠+_{𝑒∈Γ𝑒}𝑞[𝐯]𝑑𝑠,_{𝐾∈ℛℎ𝜕𝐾}𝐯⋅𝜏𝐧𝑑𝑠=_{𝑒∈Γ0𝑒}[𝜏]⋅{𝐯}𝑑𝑠+_{𝑒∈Γ𝑒}{𝜏}∶⟦𝐯⟧𝑑𝑠. (2.11)(2.12) Let ∫_Γ𝑞𝑑𝑠=∑_𝑒∈Γ∫_𝑒𝑞𝑑𝑠. Using (2.7), (2.12), and the fact that [∇𝐮]=0 for 𝐮∈(𝐻¹₀(Ω)∩𝐻²(Ω))² on Γ₀, (2.7) becomes _{𝑇∈𝒯ℎ𝜕𝑇}𝜕𝐮𝜕𝐧⋅𝝃𝑑𝑠=_{𝐾∈ℛℎ}^𝑡_{𝑗=1𝐴𝑗+1𝐶𝐴𝑗}𝜕𝐮𝜕𝐧⋅𝝃𝑑𝑠+_Γ⟦𝝃⟧∶{∇𝐮}𝑑𝑠. (2.13) Since [𝑝]=0 for 𝑝∈𝐻¹(Ω) on Γ₀, we also have _{𝑇∈𝒯ℎ𝜕𝑇}𝑝𝝃⋅𝐧𝑑𝑠=_{𝐾∈ℛℎ}^𝑡_{𝑗=1𝐴𝑗+1𝐶𝐴𝑗}𝑝𝝃⋅𝐧𝑑𝑠+_Γ𝑝[𝝃]𝑑𝑠. (2.14) Let 𝑉(ℎ)=𝑉_ℎ+(𝐻²(Ω)∩𝐻¹₀(Ω))². Define a mapping 𝛾∶𝑉(ℎ)→𝑊_ℎ, 𝛾𝐯|_𝑇=1ℎ_𝑒𝑒𝐯|_𝑇𝑑𝑠,∀𝑇∈𝒯_ℎ, (2.15) where ℎ_𝑒 is the length of the edge 𝑒.

We define two norms for 𝑉(ℎ) as follows: ‖|𝐯‖|²₁=‖𝐯‖²_1,ℎ+_𝑒∈Γ𝛾𝐯²_𝑒,‖|𝐯‖|²=‖|𝐯‖|²₁+_{𝐾∈ℛℎ}ℎ²_𝐾|𝐯|²_2,𝐾, (2.16) where ‖𝐯‖²_1,ℎ=|𝐯|²_0,ℎ+|𝐯|²_1,ℎ,  |𝐯|²_0,ℎ=∑_{𝐾∈ℛℎ}|𝐯|²_0,𝐾,|𝐯|²_1,ℎ=∑_{𝐾∈ℛℎ}|𝐯|²_1,𝐾, and ℎ_𝐾= diameter of 𝐾.

As in [5], the standard inverse inequality implies that there is a constant 𝐶 such that ‖|𝐯‖|≤𝐶‖|𝐯‖|₁,∀𝐯∈𝑉_ℎ. (2.17)

Lemma 2.1. There exists a positive constant 𝐶 independent of ℎ such that ℎ‖|𝐯‖|≤𝐶‖𝐯‖,‖𝐯‖≤𝐶‖|𝐯‖|,∀𝐯∈𝑉_ℎ. (2.18)

Proof. As in [4], ℎ‖|𝐯|‖_1,ℎ≤𝐶‖𝐯‖,‖𝐯‖≤𝐶‖|𝐯|‖_1,ℎ,∀𝐯∈𝑉_ℎ, (2.19) where ‖|𝐯‖|²_1,ℎ=|𝐯|²_1,ℎ+∑_𝑒∈Γ⟦𝛾𝐯⟧²_𝑒+∑_{𝐾∈ℛℎ}ℎ²_𝐾|𝐯|²_2,𝐾. Since ‖|𝐯‖|_1,ℎ≤‖|𝐯‖|, we have ‖𝐯‖≤𝐶‖|𝐯‖|. Note that 𝐯∈𝑉_ℎ is a piecewise linear function, and ℎ²‖|𝐯‖|²=ℎ²|𝐯|²_0,ℎ+ℎ²|𝐯|²_1,ℎ+ℎ²∑_𝑒∈Γ⟦𝛾𝐯⟧²_𝑒=𝐼₁+𝐼₂+𝐼₃. By Lemma 3.6 in [4], 𝐼₂≤𝐶‖𝐯‖²,𝐼₃≤𝐶‖𝐯‖², we have ℎ‖|𝐯‖|≤𝐶‖𝐯‖.

Lemma 2.2 (see [4]). The operator 𝛾 is self-adjoint with respect to the 𝐿²-inner product, (𝐮,𝛾𝐯)=(𝐯,𝛾𝐮),∀𝐮,𝐯∈𝑉_ℎ. Define ‖|𝐯‖|₀=(𝐯,𝛾𝐯)^1/2. Then ‖|⋅‖|₀ and ‖⋅‖ are equivalent; here the equivalence constants are independent of ℎ. And ‖𝛾𝐯‖=‖𝐯‖,∀𝐯∈𝑉_ℎ.

Let 𝑎₀(𝐯,𝜉)=(𝐯,𝜉)−_{𝐾∈ℛℎ}^𝑡_{𝑗=1𝐴𝑗+1𝐶𝐴𝑗}𝜕𝐯𝜕𝐧⋅𝜉𝑑𝑠−_Γ⟦𝝃⟧∶{}𝑐∇𝐯𝑑𝑠,𝜉,𝑞=_{𝐾∈ℛℎ}^𝑡_{𝑗=1𝐴𝑗+1𝐶𝐴𝑗}𝑞𝜉⋅𝐧𝑑𝑠+_Γ𝑞𝜉𝑏𝑑𝑠,₀𝐯,𝑞=_{𝐾∈ℛℎ𝐾}∇⋅𝐯𝑞𝑑𝑥. (2.20) It is clear that the solutions (𝐮,𝑝) of the Darcy-Stokes equations (1.1a)–(1.1c) satisfy the following: 𝑎₀(𝐮,𝜉)+𝑐𝜉,𝑝=(𝐟,𝜉),∀𝜉∈𝑊_ℎ,𝑏₀𝐮,𝑞=0,∀𝑞∈𝑄_ℎ. (2.21) Define the following bilinear forms: 𝐴₀(𝐯,𝐰)=𝑎₀(ℎ),𝐵𝐯,𝛾𝐰,∀𝐰,𝐯∈𝑉₀𝐯,𝑞=𝑏₀𝐯,𝑞,∀𝐯∈𝑉(ℎ),∀𝑞∈𝐿²₀(Ω),𝐶𝐯,𝑞=𝑐𝛾𝐯,𝑞,∀𝐯∈𝑉(ℎ),∀𝑞∈𝐿²₀(Ω). (2.22) Then systems (2.21) are equivalent to 𝐴₀(𝐮,𝐯)+𝐶𝐯,𝑝=𝐟,𝛾𝐯,∀𝐯∈𝑉_ℎ,𝐵₀𝐮,𝑞=0,∀𝑞∈𝑄_ℎ. (2.23)

We propose two discontinuous finite volume formulations based on modification of the weak formulation (2.23) for Darcy-Stokes problem (1.1a)–(1.1c). Let us introduce the bilinear forms as follows: 𝐴₁(𝐯,𝐰)=𝐴₀(𝐯,𝐰)+𝛼_𝑒∈Γ𝛾𝐯_𝑒∶𝛾𝐰_𝑒,𝐵𝐯,𝑞=𝐵₀𝐯,𝑞−_Γ𝑞𝛾𝐯𝑑𝑠, (2.24) where 𝛼>0 is a parameter to be determined later. For the exact solution (𝐮,𝑝) of (1.1a)–(1.1c), we have 𝐴₀(𝐮,𝐯)=𝐴₁(𝐮,𝐯),∀𝐯∈𝑉_ℎ,𝐵₀𝐮,𝑞=𝐵𝐮,𝑞,∀𝑞∈𝑄_ℎ. (2.25) Therefore, it follows from (2.23) that 𝐴₁(𝐮,𝐯)+𝐶𝐯,𝑝=𝐟,𝛾𝐯,∀𝐯∈𝑉_ℎ,𝐵𝐮,𝑞=0,∀𝑞∈𝑄_ℎ. (2.26) The corresponding discontinuous finite volume scheme seeks (𝐮_ℎ,𝑝_ℎ)∈𝑉_ℎ×𝑄_ℎ, such that 𝐴₁𝐮_ℎ,𝐯+𝐶𝐯,𝑝_ℎ=𝐟,𝛾𝐯,∀𝐯∈𝑉_ℎ,𝐵𝐮_ℎ,𝑞=0,∀𝑞∈𝑄_ℎ. (2.27)

Let 𝑒 be an edge of element 𝐾. It is well known (see [2]) that there exists a constant 𝐶 such that for any function 𝑔∈𝐻²(𝐾), 𝑔²_𝑒≤𝐶ℎ_𝐾⁻¹𝑔²_𝐾+ℎ_𝐾𝑔²_1,𝐾,𝜕𝑔𝜕𝐧²_𝑒≤𝐶ℎ_𝐾⁻¹𝑔²_1,𝐾+ℎ_𝐾𝑔²_2,𝐾, (2.28)(2.29) where 𝐶 depends only on the minimum angle of 𝐾.

Let ∇_ℎ𝐯 and ∇_ℎ⋅𝐯 be the functions whose restriction to each element ∀𝐾∈ℛ_ℎ is equal to ∇𝐯 and ∇⋅𝐯, respectively.

Lemma 2.3. For 𝐯,𝐰∈𝑉(ℎ), there exists a positive constant 𝐶 independent of ℎ such that 𝐴₁(𝐯,𝐰)≤𝐶‖|𝐯‖|‖|𝐰‖|. (2.30)

Proof. Let 𝐴_∗∗(𝐯,𝐰)=(𝐯,𝛾𝐰)+𝐴_∗(𝐯,𝐰), 𝐴_∗(𝐯,𝐰)=−_{𝐾∈ℛℎ}^𝑡_{𝑗=1𝐴𝑗+1𝐶𝐴𝑗}𝜕𝐯𝜕𝐧⋅𝛾𝐰𝑑𝑠. (2.31) By Lemma 3.1 in [5], 𝐴_∗(𝐯,𝐰)=∇_ℎ𝐯,∇_ℎ𝐰+_{𝐾∈ℛℎ𝜕𝐾}𝛾𝐰−𝐰𝜕𝐯𝜕𝐧𝑑𝑠+_{𝐾∈ℛℎ}Δ𝐯,𝐰−𝛾𝐰_𝐾.𝐴_∗∗(𝐯,𝐰)≤𝐯,𝛾𝐰+∇_ℎ𝐯,∇_ℎ𝐰+_{𝐾∈ℛℎ𝜕𝐾}𝛾𝐰−𝐰𝜕𝐯+𝜕𝐧𝑑𝑠_{𝐾∈ℛℎ}Δ𝐯,𝐰−𝛾𝐰_𝐾≤𝐶|𝐯|_0,ℎ|𝐰|_0,ℎ+|𝐯|_1,ℎ|𝐰|_1,ℎ+_{𝐾∈ℛℎ}ℎ_𝐾⁻¹𝐰−𝛾𝐰²_𝐾+ℎ_𝐾𝐰−𝛾𝐰²_1,𝐾^1/2×ℎ_𝐾⁻¹|𝐯|²_1,𝐾+ℎ_𝐾|𝐯|²_2,𝐾^1/2+_{𝐾∈ℛℎ}ℎ_𝐾|𝐯|_2,𝐾|𝐰|_1,𝐾≤𝐶|𝐯|_0,ℎ|𝐰|_0,ℎ+|𝐯|_1,ℎ|𝐰|_1,ℎ+_{𝐾∈ℛℎ}|𝐰|²_1,𝐾^1/2×|𝐯|_1,ℎ+_{𝐾∈ℛℎ}ℎ²_𝐾|𝐯|²_2,𝐾^1/2+_{𝐾∈ℛℎ}ℎ²_𝐾|𝐯|²_2,𝐾^1/2|𝐰|_1,ℎ𝐯𝐰,𝐴≤𝐶‖|‖|‖|‖|₁(𝐯,𝐰)=𝐴_∗∗(𝐯,𝐰)−_Γ𝛾𝐰∶{∇𝐯}𝑑𝑠+𝛼_𝑒∈Γ𝛾𝐯_𝑒∶𝛾𝐰_𝑒,≤𝐶‖|𝐯‖|‖|𝐰‖|+_{𝐾∈ℛℎ}|𝐯|²_1,𝐾+ℎ²_𝐾|𝐯|²_2,𝐾^1/2_𝑒∈Γ𝛾𝐰²_𝑒^1/2+𝛼_𝑒∈Γ𝛾𝐯²_𝑒^1/2_𝑒∈Γ𝛾𝐰²_𝑒^1/2≤𝐶‖|𝐯‖|‖|𝐰‖|. (2.32)

Lemma 2.4 (see [5]). For any (𝐯,𝑞)∈𝑉(ℎ)×𝑄_ℎ, one has 𝐶𝐯,𝑞=−𝐵𝐯,𝑞. (2.33)

Lemma 2.5 (see [5]). For (𝐯,𝑞)∈𝑉(ℎ)×𝐿²₀(Ω), there exists a positive constant 𝑀 independent of ℎ such that 𝐶𝐯,𝑞≤𝑀‖|𝐯‖|𝑞+_{𝐾∈ℛℎ}ℎ²_𝐾𝑞²_1,𝐾^1/2. (2.34) If (𝐯,𝑞)∈𝑉_ℎ×𝑄_ℎ, then 𝐶𝐯,𝑞≤𝑀‖|𝐯‖|𝑞. (2.35)

Lemma 2.6. For any 𝐯∈𝑉_ℎ, there is a constant 𝐶₀ independent of ℎ such that for 𝛼 large enough 𝐴₁(𝐯,𝐯)≥𝐶₀‖|𝐯‖|². (2.36)

Proof. Using the proof of Lemmas 3.1 and 3.5 in [5], for 𝐯∈𝑉_ℎ, 𝛾𝐯∶⟦_Γ∇𝐯⟧𝑑𝑠≤𝐶‖|𝐯‖|_1𝑒∈Γ𝛾𝐯²_𝑒^1/2,𝐴_∗(𝐯,𝐰)=∇_ℎ𝐯,∇_ℎ𝐰,∀𝐯,𝐰∈𝑉_ℎ, (2.37) we have 𝐴₁(𝐯,𝐯)=𝐯,𝛾𝐯+∇_ℎ𝐯,∇_ℎ𝐯+𝛼_𝑒∈Γ𝛾𝐯²_𝑒−_Γ∶{}≥|𝐯|𝛾𝐯∇𝐯𝑑𝑠,_20,ℎ+|𝐯|²_1,ℎ+𝛼_𝑒∈Γ𝛾𝐯²_𝑒−𝐶‖|𝐯‖|_1𝑒∈Γ𝛾𝐯²_𝑒^1/2≥𝐶‖|𝐯‖|²₁≥𝐶₀‖|𝐯‖|², (2.38) when 𝛼 is large enough.

The value of 𝛼 depends on the constant in the inverse inequality. Therefore, the value of 𝛼 for which 𝐴₁(⋅,⋅) is coercive is mesh dependent. We introduce a second discontinuous finite volume scheme which is parameter insensitive. Define a bilinear form as follows: 𝐴₂(𝐯,𝐰)=𝐴₁(𝐯,𝐰)+_Γ𝛾𝐯∶{∇𝐰}𝑑𝑠. (2.39) Similar to the bilinear form 𝐴₁(⋅,⋅), for the exact solution (𝐮,𝑝) of the Darcy-Stokes problem we have 𝐴₂(𝐮,𝐯)=𝐴₀(𝐮,𝐯),∀𝐯∈𝑉_ℎ. (2.40) Consequently, the solution of the Darcy-Stokes problem satisfies the following variational equations: 𝐴₂(𝐮,𝐯)+𝐶𝐯,𝑝=𝐟,𝛾𝐯,∀𝐯∈𝑉_ℎ,𝐵𝐮,𝑞=0,∀𝑞∈𝑄_ℎ. (2.41) Our second discontinuous finite volume scheme for (1.1a)–(1.1c) seeks (𝐮_ℎ,𝑝_ℎ)∈𝑉_ℎ×𝑄_ℎ, such that 𝐴₂𝐮_ℎ,𝐯+𝐶𝐯,𝑝_ℎ=𝐟,𝛾𝐯,∀𝐯∈𝑉_ℎ,𝐵𝐮_ℎ,𝑞=0,∀𝑞∈𝑄_ℎ. (2.42)

For any value of 𝛼>0, we have 𝐴₂(𝐯,𝐯)=𝐯,𝛾𝐯+∇_ℎ𝐯,∇_ℎ𝐯+𝛼_𝑒∈Γ𝛾𝐯²_𝑒≥𝐶‖|𝐯‖|²₁≥𝐶₀‖|𝐯‖|²,∀𝐯∈𝑉_ℎ. (2.43) Similarly, we can prove that 𝐴₂(𝐯,𝐰)≤𝐶‖|𝐰‖|‖|𝐯‖|,∀𝐯,𝐰∈𝑉(ℎ). (2.44) Let 𝐴(𝐯,𝐰)=𝐴₁(𝐯,𝐰) or 𝐴(𝐯,𝐰)=𝐴₂(𝐯,𝐰). In the rest of the paper, we assume that the following is true: 𝐴(𝐯,𝐯)≥𝐶₀‖|𝐯‖|². (2.45) If 𝐴(𝐯,𝐰)=𝐴₂(𝐯,𝐰), (2.45) holds for any 𝛼>0. If 𝐴(𝐯,𝐰)=𝐴₁(𝐯,𝐰), (2.45) holds for only 𝛼 large enough.

3. Error Estimates

We will derive optimal error estimates for velocity in the norm ‖|⋅‖| and for pressure in the 𝐿²-norm. A first-order error estimate for velocity in 𝐿²-norm will be obtained.

Let 𝑒 be an interior edge shared by two elements 𝐾₁ and 𝐾₂ in ℛ_ℎ. If ∫_𝑒𝐯|_𝐾1𝑑𝑠=∫_𝑒𝐯|_𝐾2𝑑𝑠, we say that 𝐯 is continuous on 𝑒. We say that 𝐯 is zero at 𝑒∈𝜕Ω if ∫_𝑒𝐯𝑑𝑠=0. Define a subspace 𝑉_ℎ of 𝑉_ℎ by 𝑉_ℎ=𝐯∈𝐿²(Ω)²∶𝐯|_𝐾∈𝑄₁(𝐾)²∀𝐾∈ℛ_ℎiscontinuousat𝑒∈Γ₀andiszeroat𝑒∈𝜕Ω (3.1) for rectangular meshes and by 𝑉_ℎ=𝐯∈𝐿²(Ω)²∶𝐯|_𝐾∈𝑃₁(𝐾)²∀𝐾∈ℛ_ℎiscontinuousat𝑒∈Γ₀andiszeroat𝑒∈𝜕Ω (3.2) for triangular mesh.

It has been proven in [8, 9] that the following discrete inf-sup condition is satisfied; that is, there exists a positive constant 𝛽₀ such that sup_𝐯∈𝑉_ℎ∇_ℎ⋅𝐯,𝑞|𝐯|_1,ℎ≥𝛽₀𝑞,∀𝑞∈𝑄_ℎ. (3.3)

Lemma 3.1. The bilinear form 𝐵(⋅,⋅) satisfies the discrete inf-sup condition sup_{𝐯∈𝑉ℎ}𝐵𝐯,𝑞𝐯𝑞‖|‖|≥𝛽,∀𝑞∈𝑄_ℎ, (3.4) where 𝛽 is a positive constant independent of the mesh size ℎ.

Proof. For 𝐯∈𝑉_ℎ⊂𝑉_ℎ and 𝑞∈𝑄_ℎ, we have 𝐵(𝐯,𝑞)=(∇_ℎ⋅𝐯,𝑞), and ‖|𝐯‖|₁=‖𝐯‖_1,ℎ. By Poincare-Friedrichs ‖𝐯‖_1,ℎ≤𝐶|𝐯|_1,ℎ, with (3.3), and (2.17) we get for any 𝑞∈𝑄_ℎ𝛽₀𝑞≤sup_𝐯∈𝑉_ℎ∇⋅𝐯,𝑞|𝐯|_1,ℎ≤𝐶sup_𝐯∈𝑉_ℎ𝐵𝐯,𝑞‖|𝐯‖|₁≤𝐶₁sup_{𝐯∈𝑉ℎ}𝐵𝐯,𝑞𝐯.‖|‖| (3.5) With 𝛽=𝛽₀/𝐶₁, we have proven (3.4).

Define an operator 𝜋_𝐾∶𝐻¹(𝐾)→𝑃₁(𝐾) or 𝑄₁(𝐾). For all 𝑣∈𝐻¹(𝐾), _𝑒𝑖𝜋_𝐾𝑣𝑑𝑠=_𝑒𝑖𝑣𝑑𝑠,𝑖=1,…,𝑡, (3.6) where 𝑒_𝑖,𝑖=1,…,𝑡, are the 𝑡 sides of the element 𝐾. 𝑡=3 if 𝐾 is a triangle and 𝑡=4 if 𝐾 is a rectangle. It was proven in [8] that 𝜋_𝐾𝑣−𝑣_𝑠,𝐾≤𝐶ℎ^2−𝑠|𝑣|_2,𝐾,𝑠=0,1,2. (3.7) For all 𝐯=(𝑣₁,𝑣₂)∈𝐻¹₀(Ω)², define Π₁𝐯=(Π₁𝐯₁,Π₁𝐯₂)∈𝑉_ℎ by Π₁𝑣_𝑖|_𝐾=𝜋_𝐾𝑣_𝑖,∀𝐾∈ℛ_ℎ,𝑖=1,2. (3.8) Using the definition of Π₁ and integration by parts, we can show that 𝐵𝐯−Π₁𝐯,𝑞=0,∀𝑞∈𝑄_ℎ. (3.9) The Cauchy-Schwarz inequality implies 𝛾𝐯²_𝑒=1ℎ_𝑒𝑒⟦𝐯⟧𝑑𝑠²≤1ℎ_𝑒²_𝑒⟦𝐯⟧²𝑑𝑠_𝑒𝑑𝑠=_𝑒1ℎ_𝑒⟦𝐯⟧²𝑑𝑠. (3.10) Equations (2.28) and (3.8) imply that _𝑒∈Γ𝛾𝐮−Π₁𝐮²_𝑒≤𝐶𝐮−Π₁𝐮²_1,ℎ+_{𝐾∈ℛℎ}ℎ⁻²𝐮−Π₁𝐮²_𝐾≤𝐶ℎ²‖𝐮‖²₂. (3.11) The definitions of the norm ‖|⋅‖|, (3.7), and (3.11) give |𝐮−Π₁𝐮‖²=𝐮−Π₁𝐮²_0,ℎ+𝐮−Π₁𝐮²_1,ℎ+_𝑒∈Γ𝛾𝐮−Π₁𝐮²_𝑒+_{𝐾∈ℛℎ}ℎ²𝐮−Π₁𝐮²_2,𝐾≤𝐶ℎ²‖𝐮‖²₂. (3.12)

Theorem 3.2. Let (𝐮_ℎ,𝑝_ℎ)∈𝑉_ℎ×𝑄_ℎ be the solution of (2.27), and let (𝐮,𝑝)∈(𝐻²(Ω)∩𝐻¹₀(Ω))²×(𝐿²₀(Ω)∩𝐻¹(Ω)) be the solution of (1.1a)–(1.1c). Then there exists a constant 𝐶 independent of ℎ such that 𝐮−𝐮_ℎ+𝑝−𝑝_ℎ≤𝐶ℎ‖𝐮‖₂+𝑝₁,𝐮−𝐮_ℎ≤𝐶ℎ‖𝐮‖₂+𝑝₁. (3.13)(3.14)

Proof. Let 𝜀=𝐮−Π₁𝐮,𝜀_ℎ=𝐮_ℎ−Π₁𝐮,𝜂=𝑝−Π₂𝑝,𝜂_ℎ=𝑝_ℎ−Π₂𝑝, where Π₂ is 𝐿² projection from 𝐿²₀(Ω)→𝑄_ℎ. Then 𝐮−𝐮_ℎ=𝜀−𝜀_ℎ,𝑝−𝑝_ℎ=𝜂−𝜂_ℎ. Subtracting (2.26) from (2.27) and using Lemma 2.4, we get error equations𝐴𝜀_ℎ,𝐯−𝐵𝐯,𝜂_ℎ=𝐴(𝜺,𝐯)+𝐶𝐯,𝜂,∀𝐯∈𝑉_ℎ,𝐵𝜺_ℎ,𝑞=𝐵𝜺,𝑞=0,∀𝑞∈𝑄_ℎ. (3.15a)(3.15b)By letting 𝐯=𝜀_ℎ in (3.15a) and 𝑞=𝜂_ℎ in (3.15b), the sum of (3.15a) and (3.15b) gives 𝐴𝜺_ℎ,𝜺_ℎ=𝐴𝜺,𝜺_ℎ+𝐶𝜺_ℎ,𝜂. (3.16) Thus, it follows from the coercivity (2.45), the boundedness (2.30), (2.44), and (2.34) that |𝜺_ℎ‖²≤𝐶‖|𝜺‖|𝜺_ℎ+𝜂+_{𝐾∈ℛℎ}ℎ²_𝐾𝜂²_1,𝐾^1/2𝜺_ℎ, (3.17) which implies the following: 𝜺_ℎ≤𝐶𝜺+𝜂+_{𝐾∈ℛℎ}ℎ²_𝐾𝜂²_1,𝐾^1/2. (3.18) The previous estimate can be rewritten as 𝐮_ℎ−Π₁𝐮≤𝐶𝐮−Π₁𝐮+𝑝−Π₂𝑝+_{𝐾∈ℛℎ}ℎ²_𝐾𝑝−Π₂𝑝²_1,𝐾^1/2. (3.19) Now using the triangle inequality, (3.7), the definition of Π₂, and the inequality mentioned previously, we get 𝐮−𝐮_ℎ≤𝐶𝐮−Π₁𝐮+𝐮_ℎ−Π₁𝐮≤𝐶ℎ‖𝐮‖₂+𝑝₁, (3.20) which completes the estimate for the velocity approximation.
Discrete inf-sup condition (3.4), (3.15a), (3.15b), Lemmas 2.5, 2.4, and inverse inequality give 𝑝_ℎ−Π₂𝑝≤1𝛽sup_{𝐯∈𝑉ℎ}𝐵𝐯,Π₂𝑝−𝑝_ℎ‖|𝐯‖|_1,ℎ=1𝛽sup_{𝐯∈𝑉ℎ}𝐶𝐯,𝑝_ℎ−Π₂𝑝‖|𝐯‖|_1,ℎ=1𝛽sup_{𝐯∈𝑉ℎ}𝐶𝐯,𝑝_ℎ−𝑝+𝐶𝐯,𝑝−Π₂𝑝𝐯=1‖|‖|𝛽sup_{𝐯∈𝑉ℎ}𝐴𝐮−𝐮_ℎ,𝐯+𝐶𝐯,𝑝−Π₂𝑝‖|𝐯‖|≤𝐶𝐮−𝐮_ℎ+𝑝−Π₂𝑝+_{𝐾∈ℛℎ}ℎ²_𝐾𝑝−Π₂𝑝²_1,𝐾^1/2≤𝐶ℎ‖𝐮‖₂+𝑝₁. (3.21) Using the previous inequality and the triangle inequality, we have completed the proof of (3.13).
Using Lemma 2.1, (3.12), and (3.13), we have 𝐮_ℎ−Π₁𝐮_ℎ≤𝐶𝐮_ℎ−Π₁𝐮_ℎ≤𝐶𝐮−𝐮_ℎ+𝐮−Π₁𝐮_ℎ≤𝐶ℎ‖𝐮‖₂+𝑝₁. (3.22) Equations (3.22) and (3.7) and the triangle inequality imply (3.14). We have completed the proof.

4. Numerical Experiments

In this section, we present a numerical example for solving the problems (1.1a)–(1.1c) by using the discontinuous finite volume element method presented with (2.27) and (2.42). Let Ω=(0,1)×(0,1), ℛ_ℎ be the Delaunay triangulation generated by EasyMesh [10] over Ω with mesh size ℎ as shown in Figure 5. We consider the case of 𝜎=1,𝜇=1, the exact velocity 𝑢₁(𝑥,𝑦)=−𝑥²(𝑥−1)²𝑦(𝑦−1)(2𝑦−1), 𝑢₂(𝑥,𝑦)=−𝑢₁(𝑦,𝑥) and the pressure 𝑝(𝑥,𝑦)=(𝑥−0.5)(𝑦−0.5). Denote the numerical solution as 𝐮_ℎ and 𝑝_ℎ with step ℎ_𝑑 which is used to generate the mesh data in the EasyMesh input file, and ℎ=max{ℎ_𝑒∶𝑒∈Γ}. For 𝛼=2, the numerical results are presented in Tables 1 and 2. It is observed from the tables that the numerical results support our theory.

Table 1 

Error behavior for scheme (2.27).

Table 2 

Error behavior for scheme (2.42).

Figure 5 

Triangular and its dual partition of (0,1)×(0,1).

Acknowledgments

This paper is supported by the Excellent Young and Middle-Aged Scientists Research Fund of Shandong Province (2008BS09026), National Natural Science Foundation of China (11171193), and National Natural Science Foundation of Shandong Province (ZR2011AM016).