Discontinuous Galerkin Immersed Finite Volume Element Method for Anisotropic Flow Models in Porous Medium

Abstract

By choosing the trial function space to the immersed finite element space and the test function space to be piecewise constant function space, we develop a discontinuous Galerkin immersed finite volume element method to solve numerically a kind of anisotropic diffusion models governed by the elliptic interface problems with discontinuous tensor-conductivity. The existence and uniqueness of the discrete scheme are proved, and an optimal-order energy-norm estimate and 𝐿²-norm estimate for the numerical solution are derived.

1. Introduction

Let us consider the following elliptic interface problems in a convex domain Ω⊂𝑅²: −∇⋅(𝔹∇𝑢)=𝑓,𝑥,𝑦∈Ω,𝑢=0,𝑥,𝑦∈𝜕Ω, (1) where Ω is separated into two subdomains Ω⁺ and Ω⁻ by interface Γ∈𝐶², see Figure 1 for an illustration, and 𝑓∈𝐿²(Ω); 𝑢 satisfies the following homogenous jump conditions on the interface Γ: [𝑢],𝔹𝜕𝑢𝜕𝐧=0. (2)

Figure 1 

Equation (1) describes many real diffusion processes in fluid dynamics and engineering applications, such as the miscible displacement with discontinues conductivity due to complex strata or multiphase flux. It is significant to seek efficiently the numerical solution to the interface problems for better understanding of the mechanism of the flow process and conducting engineering practice.

When 𝔹(𝑥) is a scale function, which corresponds to an isotropic flow case, two classes of numerical methods were developed to approximate (1) in terms of the meshes. One is the fitted finite element or fitted difference method [1–3], which restricts the mesh to be aligned with the smooth interface Γ. Consequently, the fitted methods are costly for more complicated time dependent problems in which the interface moves with time and repeated grid generation is called for. The other one is the immersed interface difference or finite element methods in which the Cartesian grid is naturally used even though it cannot match a curved interface. Although the immersed difference methods [4, 5] were demonstrated to be very effective, convergence analysis of related finite difference methods is extremely difficult and is still open. On the other hand, the immersed finite element method (IFE) was developed, which allows the interface to go through the interior of the element; see the references [6–9] and the references therein. Numerical experiments demonstrated an optimal order of the errors. Once again, it is not easy to analyze this method. Further, to preserve the conservative characteristics of the interface model (1), [10] developed an immersed finite volume element (IFVE) method by combining the finite volume element method [11–16] and the immersed finite element method.

In realistic diffusion processes, the interface problem (1) displays much often its anisotropic type. That is, the conductivity 𝔹(𝑥) becomes a tenser-formed function. The goals of this paper are as follows: (1) to develop a discontinuous Galerkin-immersed interface-finite volume element (DGIFVE) method for the second-order elliptic problems with tensor-formed conductivity 𝔹(𝐱) defined by 𝔹(𝐱)=𝔹^𝑙(𝐱)=𝑚^𝑙𝑠^𝑙𝑠^𝑙𝑛^𝑙, (3) where 𝐱∈Ω^𝑙,𝑚^𝑙>0,𝑛^𝑙>0,𝑚^𝑙𝑛^𝑙>(𝑠^𝑙)²,𝑙=+,−; by doing so, we can use the ability of the penalty term in discontinuous Galerkin method to control the integrals on an element boundary, in order to prove the solvability of the scheme and derive easily an optimal-order error analysis, and we can use the conservation characteristics of the finite volume element method to construct a conservation-preserved numerical method; (2) to prove the existence and uniqueness of the proposed discontinuous Galerkin-immersed interface-finite volume element procedure based on the nonconforming interface finite element space for anisotropic flow model [17]; (3) to establish its optimal-order energy-norm estimate and 𝐿²-norm estimate.

This paper is organized as follows. In the next section, we will introduce the trial function space and its approximation properties on primal triangulation. In Section 3, we will formulate the DGIFVE procedure. In Section 4, we will introduce some important lemmas. In Section 5, we will prove the existence and uniqueness of the solution of the discrete scheme. In Section 6, we will derive the convergence analysis.

Throughout this paper, the symbol 𝐶 will be used as a generic positive constant independent of ℎ and may have different values at different places.

2. The Construction of the Trial Function Space

In this section, we recall the definitions of IFE spaces discussed in [7]. Let 𝒯_ℎ={𝐾} be a regular triangulation of Ω with the diameter size ℎ. We can separate the triangles on a partition into two classes:(1)interface element: the interface Γ passes through the interior of 𝐾;(2)noninterface element: the interface does not intersect with this triangle, or it intersects with this triangle but does not separate its interior into two nontrivial subsets. Let 𝒯^𝑛_ℎ be the collection of all noninterface elements and let 𝒯^𝑚_ℎ be the collection of all interface elements. We will use 𝐴_𝑖=(𝑥_𝑖,𝑦_𝑖),  𝑖=1,2,3 to denote the vertices of 𝑇, and we will use 𝐷𝐸 to denote the line segment connecting the intersection of the interface and the edges of a triangle 𝐾. This line segment 𝐷𝐸 divides 𝑇 into two parts 𝐾⁺ and 𝐾⁻ with 𝐾=𝐾⁺∪𝐾⁻∪𝐷𝐸 (see Figure 2).

Figure 2 

For the analysis, we introduce the spaces 𝐻²(Ω)=𝑢∈𝐻¹₀(Ω)∶𝑢∈𝐻²Ω^𝑠,𝐻,𝑠=+,−₂(𝐾)=𝑢∈𝐻¹(𝐾)∶𝑢∈𝐻²𝐾∩Ω^𝑠,𝑠=+,− (4) equipped with the norms ‖𝑢‖²_𝐻2(Ω)=‖𝑢‖²_𝐻2(Ω+)+‖𝑢‖²_𝐻2(Ω−),‖𝑢‖²_𝐻2(𝐾)=‖𝑢‖²_{𝐻2(𝐾∩Ω+)}+‖𝑢‖²_{𝐻2(𝐾∩Ω−)}, (5) where 𝐻^𝑚=𝑊^𝑚₂(𝑚=1,2) is the standard Sobolev spaces. In order to define the bilinear formulation, we introduce the broken Sobolev space 𝐻¹(𝒯_ℎ)={𝑢∈𝐿²(Ω)∶∀𝐾∈𝒯_ℎ,𝑢|_𝐾∈𝐻¹(𝐾) and ∀𝑒⊂𝜕𝐾∩𝜕Ω≠0/,𝑢|_𝑒=0}.

For a noninterface element 𝐾, we use the standard linear shape functions on 𝐾 whose degrees of freedom are functional values on the vertices of 𝐾, and we use 𝑆_ℎ(𝐾) to denote the linear spaces spanned by the three nodal basis functions on 𝐾 as follows: 𝑆_ℎ(𝐾)=span𝜙_𝑖∈𝑃₁(𝐾)∶𝑖=1,2,3, (6) where 𝑃₁(𝐾) is the linear space on 𝐾. For this space, we have the following estimate of the interpolant: 𝑢−Π_ℎ𝑢_𝐿2(𝐾)+ℎ𝑢−Π_ℎ𝑢_𝐻1(𝐾)⩽𝐶ℎ²‖𝑢‖_𝐻²_(𝐾), (7) where Π_ℎ∶𝐻²(𝐾)→𝑆_ℎ(𝐾) is the interpolation operator. We use 𝑆_ℎ(Ω) to denote the space of the conforming piecewise linear polynomials on the domain 𝒯^𝑛_ℎ.

For an interface element 𝐾 whose geometric configuration is given in Figure 3 in which 𝐴₁=(0,0)^𝑇,𝐴₂=(ℎ₁,0)^𝑇,𝐴₃=(0,ℎ₂)^𝑇, the interface points 𝐷=(𝑏ℎ₁,0)^𝑇 and 𝐸=(0,𝑎ℎ₂)^𝑇, where 0<𝑎⩽1 and 0<𝑏⩽1. Let 𝜙_𝑖,  𝑖=1,2,3, denote the usual Lagrange nodal basis function associated with the vertex 𝐴_𝑖,  𝑖=1,2,3, respectively. Here we assume that the ratio 𝑟=ℎ₁/ℎ₂ is bounded below and above by some constants.

Figure 3 

By 𝜙(𝐴_𝑖)=𝑉_𝑖,  𝑖=1,2,3, we can construct the basis function 𝜙(𝐱) on an interface element 𝐾 as follows: 𝜙(𝐱)=𝜙⁻(𝐱)=𝑉₁𝜙₁+𝐶₁𝜙₂+𝐶₂𝜙₂,𝐱=𝑥,𝑦∈𝐾⁻,𝜙⁺(𝐱)=𝐶₃𝜙₁+𝑉₂𝜙₂+𝑉₃𝜙₂,𝐱=𝑥,𝑦∈𝐾⁺. (8) Satisfy 𝜙⁺(𝐷)=𝜙⁻(𝐷),𝜙⁺(𝐸)=𝜙⁻(𝐸),𝔹⁺𝜕𝜙⁺𝜕𝐧_𝐷𝐸=𝔹⁻𝜕𝜙⁻𝜕𝐧_𝐷𝐸, (9) where 𝐧_𝐷𝐸 is the unit normal vector on the line segment 𝐷𝐸.

By [17, 18], we have the following conclusions.

Lemma 1. When 𝑠^𝑙⩾0,𝑙=+,−, the piecewise linear function 𝜙(𝑥,𝑦) defined by (8) is uniquely decided by three conditions in (9).

Remark 2. By [17], the condition 𝑠^𝑙≥0,𝑙=+,− is necessary in Lemma 1. For some specially selected entries of 𝔹^𝑙 and the intersection points of the interface with the edges on 𝐾, 𝜙(𝐱) satisfying (9) is uniquely undetermined by 𝜙(𝐴_𝑖)=𝑉_𝑖,𝑖=1,2,3.
Based on the above results, the finite element space 𝑆_ℎ(𝐾) on a typical interface element 𝐾∈𝒯^𝑚_ℎ is defined by 𝑆_ℎ(𝐾)=𝜙∶𝜙ispiecewiselinearandsatisfies(9). (10) We call 𝑆_ℎ(𝐾) the immersed interface element space. For any 𝑢∈𝐻²(Ω) and 𝐾∈𝒯^𝑚_ℎ, we define (Π_ℎ𝑢)|_𝐾=Π_ℎ(𝑢|_𝐾)∈𝑆_ℎ(𝐾) by Π_ℎ𝑢𝐴_𝑖=𝑢𝐴_𝑖,𝑖=1,2,3, (11) and we call Π_ℎ𝑢 the interpolant of 𝑢 in 𝑆_ℎ(𝐾). Similar to [7], we have an estimate of the interpolant given in the following theorem.

Theorem 3. For ∀𝐾∈𝒯^𝑚_ℎ, there exists a constant 𝐶>0 such that the interpolation operator Π_ℎ∶𝐻²(𝐾)→𝑆_ℎ(𝐾) satisfies 𝑢−Π_ℎ𝑢_𝐿2(𝐾)+ℎ𝑢−Π_ℎ𝑢_𝐻1(𝐾)⩽𝐶ℎ²‖𝑢‖_𝐻²_(𝐾)∀𝑢∈𝐻²(𝐾). (12)

Finally, we define trial function space 𝑆_ℎ(Ω) as the collection of functions such that 𝜙_𝐾∈𝑆_ℎ(𝐾)𝜙,𝐾isanoninterfaceelement,_𝐾∈𝑆_ℎ(𝐾),𝐾isaninterfaceelement. (13) The space 𝑆_ℎ(Ω)⊂𝐿²(Ω) is a subspace of 𝐻¹(𝒯_ℎ). We also use the space 𝑆_0ℎ(Ω)={𝑣_ℎ∈𝑆_ℎ(Ω),𝑣_ℎ|_𝜕Ω=0}.

3. DGIFVE Procedure

In this section, we will construct a dual grid 𝒯^∗_ℎ based on 𝒯_ℎ. Assume that the triangulation 𝒯_ℎ is quasi-uniform. For a given triangle 𝐾∈𝒯_ℎ, we divide 𝐾∈𝒯_ℎ into three triangles by connecting the barycenter 𝑄 and the three corners of the triangle as shown in Figure 4. Let 𝒯^∗_ℎ consist of all these triangles 𝑇.

(a)

(b)

Figure 4 

(a) 𝐾 is a noninterface element and (b) 𝐾 is an interface element.

For the 𝒯^∗_ℎ, we define the test function space as follows: 𝑆^∗_ℎ(Ω)=𝜙(𝐱)𝜙(𝐱)_𝑇=constant,∀𝑇∈𝒯^∗_ℎ. (14) Analogous to the operator Π_ℎ, we introduce the interpolation operator 𝛾_ℎ∶𝐻_ℎ=𝑆_ℎ(Ω)+𝐻²(Ω)∩𝐻¹₀(Ω)→𝑆^∗_ℎ(Ω) defined by, for ∀𝑣∈𝐻_ℎ, 𝛾_ℎ𝑣|_𝑇=1|𝑒|_𝑒𝑣(𝐱)𝑑𝐱,𝐱∈𝑇. (15)

Let 𝑒 be an interior edge shared by two elements 𝐾₁ and 𝐾₂ in 𝒯_ℎ. Define the unite normal vectors 𝐧₁ and  𝐧₂ on 𝑒 pointing exterior to 𝐾₁ and 𝐾₂, respectively. For scalar 𝑣 function and vector function 𝐪, we define their average {⋅} and jump [⋅] on 𝑒, as follows (see [19]): {𝑣}|_𝑒=12𝑣|_𝜕𝐾1+𝑣|_𝜕𝐾2,[𝑣]|_𝑒=𝑣|_𝜕𝐾1⋅𝐧₁+𝑣|_𝜕𝐾2⋅𝐧₂,{𝐪}|_𝑒=12𝐪|_𝜕𝐾1+𝐪|_𝜕𝐾2,[𝐪]|_𝑒=𝐪|_𝜕𝐾1⋅𝐧₁+𝐪|_𝜕𝐾2⋅𝐧₂. (16) If 𝑒 is an edge on the boundary of Ω, we define {𝑣}|_𝑒=𝑣,[𝐪]|_𝑒=𝐪⋅𝐧. (17) Let 𝜀_ℎ denote the union of the boundaries of the triangle 𝐾 of 𝒯_ℎ and let 𝜀⁰_ℎ∶=𝜀_ℎ⧵𝜕Ω, 𝜀^∗_ℎ be the union of the boundaries cutting by the Γ. A straightforward computation gives _{𝐾∈𝒯ℎ𝜕𝐾}𝑣𝐪⋅𝐧𝑑𝑠=_{𝑒∈𝜀ℎ𝑒}[𝑣]{𝐪}𝑑𝑠+_{𝑒∈𝜀0ℎ𝑒}{𝑣}[𝐪]𝑑𝑠. (18)

We multiply (1) by 𝑣_ℎ∈𝑆^∗_ℎ(Ω); using [𝔹∇𝑢]|_Γ=0 and Green's formula, we have −_{𝑇∈𝒯∗ℎ𝜕𝑇}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠=𝑓,𝑣_ℎ, (19) where 𝐧 is the unit outward normal vector on 𝜕𝑇. Let 𝑇_𝑗∈𝒯^∗_ℎ(𝑗=1,2,3) be three triangles in 𝐾∈𝒯_ℎ. Then, we have _{𝑇∈𝒯∗ℎ𝜕𝑇}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠=_{𝐾∈𝒯ℎ𝑗=1𝜕𝑇𝑗}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠=_{𝐾∈𝒯ℎ𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠+_{𝐾∈𝒯ℎ𝜕𝐾}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠. (20) Using (18) and the fact that [𝔹∇𝑢]=0, (20) becomes _{𝑇∈𝒯∗ℎ𝜕𝑇}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠=_{𝐾∈𝒯ℎ𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠+_{𝑒∈𝜀ℎ𝑒}{𝔹∇𝑢}𝑣_ℎ𝑑𝑠. (21) By (19) and (21), we can get −_{𝐾∈𝒯ℎ𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢⋅𝐧𝑣_ℎ𝑑𝑠−_{𝑒∈𝜀ℎ𝑒}{𝔹∇𝑢}𝑣_ℎ𝑑𝑠=𝑓,𝑣_ℎ. (22)

By the definition of 𝛾_ℎ, the discontinuous Galerkin immersed finite volume element formulation is equivalent to finding 𝑢_ℎ∈𝑆_0ℎ(Ω) such that 𝑎^∗_ℎ𝑢_ℎ,𝛾_ℎ𝜔_ℎ=𝑓,𝛾_ℎ𝜔_ℎ,𝜔_ℎ∈𝑆_0ℎ(Ω), (23) where 𝑎^∗_ℎ𝑢_ℎ,𝛾_ℎ𝜔_ℎ=−_{𝐾∈𝒯ℎ𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ𝑑𝑠−_{𝑒∈𝜀ℎ𝑒}𝔹∇𝑢_ℎ𝛾_ℎ𝜔_ℎ𝑑𝑠+_{𝑒∈𝜀ℎ}𝜎₀ℎ_𝑒𝑒𝑢_ℎ𝜔_ℎ𝑑𝑠 (24) is the bilinear formulation defined on 𝑆_0ℎ(Ω)×𝑆_0ℎ(Ω), and in addition to penalty term ∑_{𝑒∈𝜀ℎ}(𝜎₀/ℎ_𝑒)∫_𝑒[𝑢_ℎ][𝜔_ℎ]𝑑𝑠, the penalty parameter 𝜎₀>0. Since [𝛾_ℎ𝑢]_𝑒=0, it is easy to see that 𝑢 satisfies the solution of (1) as follows: 𝑎^∗_ℎ𝑢,𝛾_ℎ𝜔_ℎ=𝑓𝑢,𝛾_ℎ𝜔_ℎ. (25)

Let 𝐴𝑢_ℎ,𝜔_ℎ=−_{𝐾∈𝒯ℎ}³_{𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ𝑑𝑠. (26) If 𝐾∈𝒯^𝑛_ℎ, we have by ∇⋅(𝔹∇𝑢_ℎ)=0,𝑢_ℎ∈𝑆_0ℎ(Ω)0=_𝐾∇⋅𝔹∇𝑢_ℎ𝛾_ℎ𝜔_ℎ=𝑑𝑥_{3𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ+𝑑𝑠_{3𝑗=1𝐴𝑗𝐴𝑗+1}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ𝑑𝑠. (27) Then, −³_{𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ𝑑𝑠=³_{𝑗=1𝐴𝑗𝐴𝑗+1}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ𝑑𝑠. (28) Similarly, 0=−_𝐾∇⋅𝔹∇𝑢_ℎ𝜔_ℎ𝑑𝑥=−³_{𝑗=1𝐴𝑗𝐴𝑗+1}𝔹∇𝑢_ℎ⋅𝐧𝜔_ℎ𝑑𝑠+_𝐾𝔹∇𝑢_ℎ⋅∇𝜔_ℎ𝑑𝑥. (29) We find that _𝜕𝐾𝔹∇𝑢_ℎ⋅𝐧𝜔_ℎ−𝛾_ℎ𝜔_ℎ𝑑𝑠=0, (30) due to the fact that 𝔹∇𝑢_ℎ is a constant vector on each edge and the definition of 𝛾_ℎ. Thus, we can get _𝐾𝔹∇𝑢_ℎ∇𝜔_ℎ𝑑𝑥=−³_{𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ𝑑𝑠 (31) following from (28) and (29). For 𝐾∈𝒯^𝑚_ℎ (see Figure 4(b)), it follows from the same arguments above and the [𝐵∇𝑢_ℎ]|_𝐷𝐸=0 that _𝐾𝔹∇𝑢_ℎ∇𝜔_ℎ𝑑𝑥=−³_{𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ+𝑑𝑠_𝐴1𝐷𝔹⁺∇𝑢⁺_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ−𝜔⁺_ℎ𝑑𝑠+_𝐷𝐴2𝔹⁻∇𝑢⁻_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ−𝜔⁻_ℎ+𝑑𝑠_𝐴2𝐸𝔹⁻∇𝑢⁻_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ−𝜔⁻_ℎ𝑑𝑠+_𝐸𝐴3𝔹⁺∇𝑢⁺_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ−𝜔⁺_ℎ𝑑𝑠. (32) Summarizing the results above, we have 𝐴𝑢_ℎ,𝜔_ℎ=−_𝐾³_{𝑗=1𝐴𝑗+1𝑄𝐴𝑗}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝜔_ℎ𝑑𝑠=_{𝐾∈𝒯ℎ𝐾}𝔹∇𝑢_ℎ∇𝜔_ℎ𝑑𝑥−_{𝐾∈𝒯𝑚ℎ𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ⋅𝐧𝜔_ℎ−𝛾_ℎ𝜔_ℎ=𝑑𝑠_{𝐾∈𝒯ℎ𝐾}𝔹∇𝑢_ℎ∇𝜔_ℎ−𝑑𝑥_{𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ𝜔_ℎ−𝛾_ℎ𝜔_ℎ−𝑑𝑠_{𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ𝜔_ℎ−𝛾_ℎ𝜔_ℎ𝑑𝑠. (33) Thus, (24) can be written by 𝑎^∗_ℎ𝑢_ℎ,𝛾_ℎ𝜔_ℎ=𝐴𝑢_ℎ,𝜔_ℎ−_{𝑒∈𝜀ℎ𝑒}𝔹∇𝑢_ℎ𝛾_ℎ𝜔_ℎ+𝑑𝑠_{𝑒∈𝜀ℎ}𝜎₀|𝑒|^𝛼0_𝑒𝑢_ℎ𝜔_ℎ𝑑𝑠. (34)

4. Some Lemmas

We define a norm |||⋅||| for 𝐻_ℎ as follows: |||𝑣|||²=|𝑣|²_1,ℎ+_𝑒1ℎ_𝑒𝑒[𝑣]²𝑑𝑠+_{𝑒∈𝜀∗ℎ}ℎ_𝑒𝑒𝛽∇𝑣²𝑑𝑠. (35) In order to prove the existence and uniqueness of the solution to (24) and conduct its convergence analysis, we need the following lemmas.

Lemma 4. The operator 𝛾_ℎ in (15) has the following properties: _𝑒𝜔−𝛾_ℎ𝜔𝑑𝑠=0,∀𝜔∈𝐻_ℎ,∀𝑒∈𝜀_ℎ;[𝜔]|_𝑒=0⟹[𝛾_ℎ𝜔]_𝑒=0,∀𝜔∈𝐻_ℎ;𝛾_ℎ𝜔−𝜔_0,𝐾⩽𝐶ℎ_𝐾|𝜔|_1,𝐾,∀𝐾∈𝒯_ℎ,𝜔∈𝐻_ℎ. (36)(37)(38)

Proof. Obviously, (36) and (37) follow from the definition of 𝛾 in (15). We only prove (38) below.
Let 𝐾 be a noninterface element; we have the conclusion (38) by [12]. Therefore, we focus (38) on interface element 𝐾 (Figure 4(b)). For ∀𝜔∈𝑆_ℎ(Ω), we have the following form: 𝜔𝜔=⁺=𝑎₀+𝑏₀𝑥+𝑐₀𝑦,𝐱=𝑥,𝑦∈𝐾⁺,𝜔⁻=𝑎₁+𝑏₁𝑥+𝑐₁𝑦,𝐱=𝑥,𝑦∈𝐾⁻. (39) The jump conditions on 𝐷𝐸 lead to (see [7]) ∇𝜔⁺=𝐧²_𝑦+𝜌𝐧²_𝑥𝜌−1𝐧_𝑥𝐧_𝑦𝜌−1𝐧_𝑥𝐧_𝑦𝐧²_𝑥+𝜌𝐧²_𝑦∇𝜔⁻, (40) or ∇𝜔⁺=𝑁⁻_𝐷𝐸∇𝜔⁻, (41) where 𝐧_𝐷𝐸=(𝐧_𝑥,𝐧_𝑦)^𝑇 and 𝜌=(𝛽⁻/𝛽⁺). We know that 𝜔−𝛾_ℎ𝜔²_0,𝐾=³_{𝑗=1𝑇𝑗}𝜔−𝛾_ℎ𝜔²𝑑𝐱, (42) where _𝑇𝑖𝜔−𝛾_ℎ𝜔²=𝑑𝐱_𝑇+𝑖𝜔⁺−𝛾_ℎ𝜔²𝑑𝐱+_𝑇−𝑖𝜔⁻−𝛾_ℎ𝜔²𝑑𝐱,𝑖=1,2. (43) Since 𝜔 is continuous on 𝐴₁𝐴₂, there exists a point 𝜉 such that 𝛾_ℎ𝜔|_𝐴1𝐴2=𝜔(𝜉). (44) We suppose that 𝜉 fall on 𝐴₁𝐷; then, we have _𝑇1𝜔−𝛾_ℎ𝜔²𝑑𝐱=_𝑇+1𝜔⁺−𝜔⁺(𝜉)²𝑑𝐱+_𝑇−1𝜔⁻−𝜔⁺(𝜉)²=𝑑𝐱_𝑇+1𝜔⁺−𝜔⁺(𝜉)²+𝑑𝐱_𝑇−1𝜔⁻−𝜔⁻(𝐷)+𝜔⁺_ℎ(𝐷)−𝜔⁺(𝜉)²𝑑𝐱, (45) where we used 𝜔⁻(𝐷)=𝜔⁺(𝐷). Because 𝜔⁺(𝐱) and 𝜔⁻(𝐱) are linear polynomial, we have 𝜔⁺(𝐱)=𝜔⁺(𝜉)+∇𝜔⁺(),𝜔𝐱−𝜉₋(𝐱)=𝜔⁻(𝐷)+∇𝜔⁻(𝐱−𝐷). (46) From these expansions of 𝜔 and (41), we have _𝑇1𝜔−𝛾_ℎ𝜔²𝑑𝐱=_𝑇+1∇𝜔⁺(𝐱−𝜉)²+𝑑𝐱_𝑇−1∇𝜔⁻(𝐱−𝐷)+𝑁⁻_𝐷𝐸∇𝜔⁻(𝐷−𝜉)²𝑑𝐱. (47) Then, _𝑇1𝜔−𝛾_ℎ𝜔²𝑑𝐱⩽𝐶ℎ²|𝜔|²_1,𝑇1. (48) If 𝜉 is on 𝐷^→𝐴₂, similarly, we also have (48). Analogously, we can have the following inequality: _𝑇𝑖𝜔−𝛾_ℎ𝜔²𝑑𝐱⩽𝐶ℎ²|𝜔|²_1,𝑇𝑖,𝑖=2,3. (49)
This completes the proof of (38) by (48) and (49).

Lemma 5. For any 𝜔∈𝐻_ℎ and 𝑒∈𝜀_ℎ, one has [𝛾_ℎ𝜔]_𝐿2(𝑒)⩽‖[𝜔]‖_𝐿²_(𝑒);[𝜔−𝛾_ℎ𝜔]_𝐿2(𝑒)⩽2‖[𝜔]‖_𝐿²_(𝑒). (50)(51)

Proof. If 𝑒 is the common side of 𝑇₁, 𝑇₂∈𝒯^∗_ℎ, and 𝑇₁⊂𝐾₁, 𝑇₂⊂𝐾₂, by the definition of 𝛾_ℎ, we have [𝛾_ℎ𝜔]²_𝐿2(𝑒)=𝛾_ℎ𝜔|_𝜕𝑇1⋅𝐧₁+𝛾_ℎ𝜔|_𝜕𝑇2⋅𝐧₂²_𝐿2(𝑒)=ℎ_𝑒⁻¹_𝑒𝜔|_𝑇1𝑑𝑠⋅𝐧₁+_𝑒𝜔|_𝑇2𝑑𝑠⋅𝐧₂²_𝐿2(𝑒)=ℎ_𝑒⁻²_𝑒[𝜔]𝑑𝑠²_𝐿2(𝑒)=ℎ_𝑒⁻¹_𝑒[𝜔]𝑑𝑠². (52) Using the Hölder inequality, we can get 𝛾_ℎ𝜔²_𝐿2(𝑒)⩽_𝑒[𝜔]²𝑑𝑠=‖[𝜔]‖²_𝐿2(𝑒). (53) If 𝑒∈𝜕Ω, there exists a 𝑇∈𝒯^∗_ℎ, such that 𝑒∈𝜕𝑇 and 𝑇⊂𝐾; we have [𝛾_ℎ𝜔]²_𝐿2(𝑒)=𝛾_ℎ𝜔⋅𝐧²_𝐿2(𝑒)=ℎ_𝑒⁻¹_𝑒𝜔|_𝑇𝑑𝑠⋅𝐧²_𝐿2(𝑒)=ℎ_𝑒⁻¹_𝑒[𝜔]𝑑𝑠²⩽‖[𝜔]‖²_𝐿2(𝑒). (54) Thus, (50) is valid.
For (51), we have [𝜔−𝛾_ℎ𝜔]_𝐿2(𝑒)=[𝜔]−[𝛾_ℎ𝜔]_𝐿2(𝑒)⩽2‖[𝜔]‖_𝐿²_(𝑒) (55) following from (50).

Lemma 6 (see [20]). Let 𝒯 be a regular triangulation; then, there exists a constant 𝐶>0 independent of ℎ_𝐾 such that, for 𝑤∈𝐻¹(𝐾) and 𝐾∈𝒯_ℎ, the following inequality holds: _𝜕𝐾|𝑤|²𝑑𝑠⩽𝐶ℎ_𝐾⁻¹‖𝑤‖²_0,𝐾+ℎ_𝐾|𝑤|²_1,𝐾. (56)

5. Existence, Uniqueness, and Convergence of DIFVE Solution

In this section, we will prove the existence and uniqueness of the solution to (24) and conduct its convergence analysis in the broken |||⋅||| norm.

Lemma 7. There is a constant 𝐶 independent of  ℎ such that for 𝜎₀ large enough and 𝜖 small enough 𝑎^∗_ℎ𝑢_ℎ,𝛾_ℎ𝑢_ℎ⩾𝐶𝑢_ℎ². (57)

Proof. By Cauchy-Schwarz's inequality, we have _{𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ⋅𝐧𝑢_ℎ−𝛾_ℎ𝑢_ℎ⩽𝑑𝑠_{𝑒∈𝜀∗ℎ}𝔹∇𝑢_ℎ⋅𝐧²_𝐿2(𝑒)^1/2_{𝑒∈𝜀∗ℎ}𝑢_ℎ−𝛾_ℎ𝑢_ℎ²_𝐿2(𝑒)^1/2. (58)
Using the trace inequality (56), we have {𝑢_ℎ−𝛾_ℎ𝑢_ℎ}²_𝐿2(𝑒)⩽12𝑢_ℎ−𝛾_ℎ𝑢_ℎ|_𝐾^𝑒₁²_𝐿2(𝑒)+𝑢_ℎ−𝛾_ℎ𝑢_ℎ|_𝐾^𝑒₂²_𝐿2(𝑒)⩽𝐶ℎ_𝐾^𝑒₁𝑢_ℎ²_1,𝐾1+ℎ_𝐾^𝑒₂𝑢_ℎ²_1,𝐾2⩽𝐶ℎ𝑢_ℎ²_1,𝐾1+𝑢_ℎ²_1,𝐾2, (59) where edge 𝑒 is shared by the elements 𝐾^𝑒₁ and 𝐾^𝑒₂. Therefore, _{𝑒∈𝜀∗ℎ}𝑢_ℎ−𝛾_ℎ𝑢_ℎ²_𝐿2(𝑒)^1/2⩽𝐶ℎ^1/2_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝑘^1/2. (60) By Young's inequality, we have, for 𝜀>0, _{𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ⋅𝐧𝑢_ℎ−𝛾_ℎ𝑢_ℎ𝑑𝑠⩽𝐶ℎ^1/2_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝑘^1/2_{𝑒∈𝜀∗ℎ}[𝔹∇𝑢_ℎ⋅𝐧]²_𝐿2(𝑒)^1/2⩽𝐶_{0𝑒∈𝜀∗ℎ}ℎ[𝔹∇𝑢_ℎ⋅𝐧]²_𝐿2(𝑒)+𝜀3_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝑘. (61) Similarly, we obtain _{𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ⋅𝐧𝑢_ℎ−𝛾_ℎ𝑢_ℎ⩽𝑑𝑠_{𝑒∈𝜀∗ℎ}𝔹∇𝑢_ℎ⋅𝐧²_𝐿2(𝑒)^1/2×_{𝑒∈𝜀∗ℎ}𝑢_ℎ−𝛾_ℎ𝑢_ℎ²_𝐿2(𝑒)^1/2. (62) On the one hand, we have {𝔹∇𝑢_ℎ⋅𝐧}²_𝐿2(𝑒)⩽𝐶ℎ⁻¹𝑢_ℎ²_1,𝐾1+𝑢_ℎ²_1,𝐾2, (63) and thus _{𝑒∈𝜀∗ℎ}𝔹∇𝑢_ℎ⋅𝐧²_𝐿2(𝑒)^1/2⩽𝐶ℎ^−1/2_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝑘^1/2. (64) On the other hand, we get _{𝑒∈𝜀∗ℎ}𝑢_ℎ−𝛾_ℎ𝑢_ℎ²_𝐿2(𝑒)^1/2⩽2_{𝑒∈𝜀∗ℎ}𝑢_ℎ²_𝐿2(𝑒)^1/2 (65) following from (51). Therefore, _{𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ⋅𝐧𝑢_ℎ−𝛾_ℎ𝑢_ℎ𝑑𝑠⩽𝐶ℎ^−1/2_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝐾^1/2_{𝑒∈𝜀∗ℎ}𝑢_ℎ²_𝐿2(𝑒)^1/2⩽𝐶ℎ^{(𝛼0−1)/2}_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝐾^1/2×_{𝑒∈𝜀∗ℎ}1|𝑒|^𝛼0𝑢_ℎ²_𝐿2(𝑒)^1/2, (66) where 𝛼₀>1. By 𝜀-inequality, we have _{𝑒∈𝜀∗ℎ𝑒}𝔹∇𝑢_ℎ⋅𝐧𝑢_ℎ−𝛾_ℎ𝑢_ℎ⩽𝜖𝑑𝑠3_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝐾+𝐶_{1𝑒∈𝜀∗ℎ}1|𝑒|^𝛼0𝑢_ℎ²_𝐿2(𝑒). (67) Similarly, we can get _{𝑒∈𝜀ℎ𝑒}𝔹∇𝑢_ℎ⋅𝐧𝛾_ℎ𝑢_ℎ⩽𝜖𝑑𝑠3_{𝐾∈𝒯ℎ}𝑢_ℎ²_1,𝐾+𝐶_{2𝑒∈𝜀ℎ}1|𝑒|^𝛼0𝑢_ℎ²_𝐿2(𝑒). (68) Combining (61), (67), and (68), we obtain 𝑎^∗_ℎ𝑢_ℎ,𝛾_ℎ𝑢_ℎ⩾𝜖𝐶−3_{𝐾∈𝒯𝑛ℎ}𝑢_ℎ²_1,𝐾+(𝐶−𝜖)_{𝐾∈𝒯𝑚ℎ}𝑢_ℎ²_1,𝐾+_{𝑒∈𝜀∗ℎ}𝜎₀−𝐶₁−𝐶₂|𝑒|^𝛼0𝑢_ℎ²_𝐿2(𝑒)+_{𝑒∈𝜀ℎ⧵𝜀∗ℎ}𝜎₀−𝐶₂|𝑒|^𝛼0[𝑢_ℎ]²_𝐿2(𝑒)+_{𝑒∈𝜀∗ℎ}𝐶₀ℎ[𝔹∇𝑢_ℎ⋅𝐧]²_𝐿2(𝑒). (69) Choosing 𝜎₀ large enough and 𝜖 small enough, we have 𝑎^∗_ℎ𝑢_ℎ,𝛾_ℎ𝑢_ℎ⩾𝐶𝑢_ℎ². (70)

Lemma 8. For 𝑢, 𝜔∈𝐻_ℎ, one has 𝑎^∗_ℎ𝑢,𝛾_ℎ𝜔⩽𝐶+|||𝑢|||+_{𝐾∈𝒯ℎ}ℎ²|𝑢|²_𝐻2(𝐾)^1/2|||𝜔|||. (71) If 𝑢_ℎ, 𝜔_ℎ∈𝑆_ℎ(Ω), then 𝑎^∗_ℎ𝑢_ℎ,𝛾_ℎ𝜔_ℎ⩽𝐶𝑢_ℎ𝜔_ℎ. (72)

Proof. By (38) and the trace inequality (56), we have, for any 𝑒^𝑠∈𝜀^∗_ℎ,𝑠=+,−, _𝑒𝑠𝔹∇𝑢⋅𝐧𝜔−𝛾_ℎ𝜔ℎ𝑑𝑠⩽𝐶⁻¹|𝑢|²_1,𝐾^𝑠+ℎ|𝑢|²_2,𝐾^𝑠1/2ℎ⁻¹𝜔−𝛾_ℎ𝜔²_{𝐿2(𝐾𝑠)}+ℎ𝜔−𝛾_ℎ𝜔²_1,𝐾𝑠^1/2⩽𝐶|𝑢|²_1,𝐾^𝑠+ℎ²|𝑢|²_2,𝐾^𝑠1/2|𝜔|_1,𝐾^𝑠⩽𝐶|𝑢|²_1,𝐾+ℎ²|𝑢|²_𝐻2(𝐾)^1/2|𝜔|_1,𝐾. (73) By the Cauchy-Schwarz inequality, we obtain |𝐴(⩽𝑢,𝜔)|_{𝐾∈𝒯ℎ𝐾}+𝔹∇𝑢∇𝜔𝑑𝑥_{𝐾∈𝒯𝑚ℎ𝑒∈𝜀∗ℎ𝑒}𝐵∇𝑢⋅𝐧𝜔−𝛾_ℎ𝜔|𝑢|𝑑𝑠⩽𝐶_1,ℎ+_{𝐾∈𝒯ℎ}|𝑢|²_1,𝐾+ℎ²|𝑢|²_𝐻2(𝐾)^1/2|𝜔|_1,ℎ. (74) The definition of 𝑎_ℎ(𝑢,𝛾_ℎ𝜔) and the inequality above imply that 𝑎_ℎ𝑢,𝛾_ℎ𝜔⩽𝐶|𝑢|_1,ℎ|𝜔|_1,ℎ+_{𝐾∈𝒯ℎ}|𝑢|²_1,𝐾+ℎ²|𝑢|²_𝐻2(𝐾)^1/2|𝜔|_1,ℎ+_{𝐾∈𝒯ℎ}|𝑢|²_1,𝐾+ℎ²|𝑢|²_𝐻2(𝐾)^1/2×_{𝑒∈𝜀ℎ}‖[𝜔]‖²_𝐿2(𝑒)^1/2+_{𝑒∈𝜀ℎ}‖[𝑢]‖²_𝐿2(𝑒)^1/2_{𝑒∈𝜀ℎ}‖[𝜔]‖²_𝐿2(𝑒)^1/2⩽𝐶|||𝑢|||+_{𝐾∈𝒯ℎ}ℎ²|𝑢|²_𝐻2(𝐾)^1/2|||𝜔|||. (75) This completes the proof of (71). For (72), we can get the following from (71) and |𝑢_ℎ|_𝐻²_(𝐾)=0, ∀𝑢_ℎ∈𝑆_ℎ(Ω).

Lemmas 7 and 8 guarantee the existence and uniqueness of the discontinuous immersed finite volume element solution to (23) when choosing 𝜎₀ large enough.

6. Error Estimates in the Energy Norm

We will derive an optimal-order error estimate in the norm |||⋅||| defined in (35) and a first order error estimate in 𝐿²-norm. We start with the following lemmas.

Lemma 9. Let u be the solutions of (1); one has the conclusion 𝑢−Π_ℎ𝑢⩽𝐶ℎ‖𝑢‖_𝐻²_(Ω). (76)

Proof. By the definition of norm |||⋅|||, we have 𝑢−Π_ℎ𝑢²=𝑢−Π_ℎ𝑢²_1,ℎ+_{𝑒∈𝜀ℎ}1|𝑒|^𝛼0[𝑢−Π_ℎ𝑢]²_𝐿2(𝑒)+_{𝑒∈𝜀∗ℎ}ℎ[∇(𝑢−Π_ℎ𝑢)⋅𝐧]²_𝐿2(𝑒). (77) Using trace inequality and (12), we have [𝑢−Π_ℎ𝑢]²_𝐿2(𝑒)⩽𝐶𝑢−Π_ℎ𝑢|_𝐾^𝑒₁²_𝐿2(𝑒)+𝑢−Π_ℎ𝑢|_𝐾^𝑒₂²_𝐿2(𝑒)⩽𝐶ℎ⁻¹𝑢−Π_ℎ𝑢²_𝐿2𝐾1+ℎ𝑢−Π_ℎ𝑢²_1,𝐾1+ℎ⁻¹𝑢−Π_ℎ𝑢²_𝐿2𝐾2+ℎ𝑢−Π_ℎ𝑢²_1,𝐾2⩽𝐶ℎ²‖𝑢‖²_𝐻2𝐾1+‖𝑢‖²_𝐻2𝐾2, (78) where 𝑒 is shared by the elements 𝐾₁ and 𝐾₂. Thus, we obtain _{𝑒∈𝜀ℎ}[𝑢−Π_ℎ𝑢]²_𝐿2(𝑒)⩽𝐶ℎ²‖𝑢‖²_𝐻2(Ω). (79) Analogously, we can get _{𝑒∈𝜀∗ℎ}ℎ[∇(𝑢−Π_ℎ𝑢)⋅𝐧]²_𝐿2(𝑒)⩽𝐶ℎ²‖𝑢‖²_𝐻2(Ω). (80) By (12), (79), and (80), we obtain (76).

Lemma 10 (see [21]). There exists a constant 𝐶 independent of ℎ such that ‖𝜔‖⩽𝐶|||𝜔|||,∀𝜔∈𝑆_ℎ(Ω). (81)

Theorem 11. Let 𝑢_ℎ∈𝑆_0ℎ(Ω) and 𝑢∈𝐻²(Ω)∩𝐻¹₀(Ω) be the solutions of (23) and (1), respectively; then, there exists a constant 𝐶 independent of ℎ such that 𝑢−𝑢_ℎ⩽𝐶ℎ‖𝑢‖_𝐻²_(Ω),𝑢−𝑢_ℎ⩽𝐶ℎ‖𝑢‖_𝐻²_(Ω). (82)(83)

Proof. Subtracting (25) from (23) gives 𝑎^∗_ℎ𝑢−𝑢_ℎ,𝛾_ℎ𝜔_ℎ=0,∀𝜔_ℎ∈𝑆_0ℎ(Ω). (84) Using (71), (76), and (84), we have 𝑢_ℎ−Π_ℎ𝑢²⩽𝐶𝑎^∗_ℎ𝑢_ℎ−Π_ℎ𝑢_ℎ,𝛾_ℎ𝑢_ℎ−Π_ℎ𝑢_ℎ=𝐶𝑎^∗_ℎ𝑢−Π_ℎ𝑢,𝛾_ℎ𝑢_ℎ−Π_ℎ𝑢_ℎ⩽𝐶𝑢−Π_ℎ𝑢+_{𝐾∈𝒯ℎ}ℎ²𝑢−Π_ℎ𝑢²_𝐻²_(𝐾)^1/2×𝑢_ℎ−Π_ℎ𝑢⩽𝐶ℎ²‖𝑢‖²_𝐻2(Ω). (85) Thus, we can get 𝑢−𝑢_ℎ⩽𝐶ℎ‖𝑢‖_𝐻²_(Ω) (86) following from triangle inequality, (76), and (85).
For (83), we can get by (12), (81) and triangle inequality. We have completed the proof.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the NSF of China under Grants 10971254 and 10926100 and the Development projects of Shandong Province Science and Technology 2012GGB02298.