A discontinuous finite volume method for a coupled fracture model

Computers and Mathematics with Applications 78 (2019) 3429–3449

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A discontinuous finite volume method for a coupled fracture model Shuangshuang Chen Beijing Institute for Scientific and Engineering Computing (BISEC), Beijing University of Technology, Beijing, 100124, China

article

info

Article history: Received 21 July 2018 Received in revised form 4 May 2019 Accepted 14 May 2019 Available online 31 May 2019 Keywords: Fracture model Discontinuous finite volume method Error estimates Numerical experiments

a b s t r a c t In this paper, we present a discontinuous finite volume method for a hybrid-dimensional coupled fracture model on triangular meshes. In matrix and fracture, the pressure is approximated by the discontinuous and continuous piecewise linear elements, and in order to retain local conservation property, the velocity is sought in the discontinuous lowest-order Raviart–Thomas element space and piecewise linear element space. The proposed numerical scheme is only associated with the pressure, so we avoid solving an indefinite saddle-point problem. In addition, the discontinuous method eliminates the continuity of approximate functions across the interelement boundary, so it has the advantage of high localizability and easy handling of complicated geometries. Optimal order error estimates for the pressure and velocity are proved, and numerical experiments are carried out to confirm our theoretical analysis. © 2019 Published by Elsevier Ltd.

1. Introduction In this paper, we are concerned with a single phase fluid governed by Darcy’s law in a domain containing fractures. Fracture models have become a very active research area based on their wide variety of applications in petroleum reservoir, nuclear waste storage monitoring, groundwater waste problem, and so on. Due to the complex geometries involved and the high contrast of length scales between matrix and fractures, the presence of fractures in porous media greatly complicates the modeling of flow and increases difficulty of numerical simulation. Thus, providing accurate and efficient numerical solutions for such fracture models is not only the requirement of numerical mathematics, but also the demand of real world applications. When it comes to very thin fracture, an idea of treating fractures as (n − 1)dimensional interfaces in the n-dimensional domain was proposed in [1], but in the paper fractures were supposed to be of high permeability so pressure was continuous across the fracture, and flux was supposed to be discontinuous since fluid could flow into and out of as well as along fractures. In [2], a new model generalizing the earlier models to handle fractures with large or small permeability was presented, in which the flow equations in the fracture were averaged along the normal direction to obtain a new set of reduced models suited for the hybrid-dimensional description of fluid flow, and the Robin boundary condition was proposed on the interface to describe the interaction between the fracture and surrounding domains. We also refer to [3,4] for similar fracture models. The reduction of dimensionality in fractures overcomes most of the difficulties associated to the problem, and fracture aperture becomes a coefficient in the equations and not a geometrical constrain for the grid generation, hence it has been a hot topic in studying the fluid flow process in fractured media. The lowest order Raviart–Thomas mixed finite element method was applied to fracture models in [2], and the existence and uniqueness of the numerical solutions were analyzed. Then the mixed Raviart–Thomas element method was extended to the nonmatching meshes between fractures and surrounding domains in [5]. The extended finite E-mail address: [email protected]. https://doi.org/10.1016/j.camwa.2019.05.014 0898-1221/© 2019 Published by Elsevier Ltd.

3430

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

element method was used in [6] for a fracture model coupling two phase flows in surrounding media and single phase flow in the fracture. A block-centered finite difference method was developed in [7] for a coupled fracture model, in which both the pressure and velocity were approximated with second-order accuracy. A node-centered finite volume method based on conforming piecewise linear elements was presented in [8]. In this paper, we will derive the discontinuous finite volume method (DFVM) for the hybrid-dimensional fracture model, in which the pressure is approximated by discontinuous piecewise linear functions in matrix and the continuity across the interelement boundary is eliminated, hence it is more suitable to handle the complicate porous domain with fractures. The finite volume method (FVM) [9,10] as a discretization technique is widely used to solve various partial differential equations (PDEs) due to its local conservation property of original problems, such as the elliptic problems [11,12], the parabolic problems [13,14], the advection–diffusion equation [15], the Stokes problem [16], and so on. Unlike the standard conforming FVM [17,18] and nonconforming FVM [19–21] with Crouzeix–Raviart element [22], the DFVM [23–25] is inspired by the discontinuous Galerkin (DG) methods [26,27]. The discontinuous piecewise polynomials and piecewise constant elements are used for trial functions and test functions in the DFVM. The continuity of the numerical solution across the interelement boundary is eliminated, and a penalty term is added to enforce connections between elements. The DFVM combines the benefits of DG methods and FVMs, such as the flexibility, highly order accuracy of DG methods, and conservative property, simplicity of FVMs. In [23], Ye first proposed the DFVM for second order elliptic problems, and an optimal order error estimate for the discrete H 1 norm of the pressure was analyzed. Then the idea was extended to solve Stokes equations on both triangular and rectangular meshes in [24]. In [25], Chou and Ye established a general framework for analyzing the class of finite volume methods, in which continuous or totally discontinuous polynomials were used as trial functions, and the optimal order convergence in the L2 norm was proved under their framework. A discontinuous finite volume method with different penalty term from [23,24] was developed in [28], in which the DFVM was treated as a perturbation of the interior penalty method and the closer relationship between the DFVM and the interior penalty method was revealed. In this paper, we will apply the DFVM on triangular meshes to the hybrid-dimensional coupled fracture model. The pressure in surrounding domains and fracture is approximated by discontinuous and continuous piecewise linear polynomials. For the velocity, in order to retain local conservation property at the element level, we use the conservative method proposed in [29] to calculate flux by the discontinuous lowest-order Raviart–Thomas elements [30] in matrix and by piecewise linear elements in the fracture. The piecewise constant elements are used as test functions, so in order to keep the same dimension for the spaces of the trial and test element, two different partitions of the computational domain in the DFVM are required, one is called the primary partition corresponding to the trial space, and the other is called the dual partition corresponding to the test space. Here we use the standard triangular elements in the primal partition. The dual partition is constructed as follows: each triangle is divided into three triangles by connecting the barycenter and three vertices, and each small triangle is regarded as a control volume corresponding to a degree of freedom on an edge. Compared to the classical conforming and nonconforming FVMs, elements in the dual partition of the DFVM have the best local property because only one element in the primary partition is needed. Therefore, the DFVM has the advantage of high localizability and flexibility to handle complicate geometries. In addition, since low-order polynomials are used as trial functions, the DFVM is easy to code, but the desired high order of accuracy is still achieved. Error estimates show that optimal order error estimates in the discrete H 1 semi-norm and L2 norm of the pressure p as well as the (L2 )2 norm and discrete H(div ) norm of the velocity u are obtained in general triangular meshes. We also carry out numerical experiments to illustrate the accuracy of our theoretical analysis. The rest of the paper is organized as follows. In Section 2, we describe the hybrid-dimensional coupled fracture model, and define the triangulation and the corresponding dual partition. The discontinuous finite volume scheme is established in Section 3. In Section 4, the stability and convergence for the proposed method are analyzed. In Section 5, numerical experiments are provided to show errors and convergence rates of the proposed scheme for the coupled fracture model. Finally, in Section 6, we give a conclusion for our work. Through this paper, we will use the letter C to denote a positive constant with independence of mesh size and indicate difference values at its difference appearances. 2. Description of the fracture model In this article, we are concerned with the hybrid-dimensional coupled fracture model developed in [2]. The fracture is regarded as one-dimensional interface in the two-dimensional domain. As shown in Fig. 1, for simplicity, we denote by Ω = (a1 , a2 ) × (b1 , b2 ), a rectangle domain in R2 with the boundary Γ . Let γ = {x = xf } × (b1 , b2 ) ⊂ Ω be the fracture, and Ω is separated by the fracture γ into two bounded sub-domains:

Ω1 = (a1 , xf ) × (b1 , b2 ), and Ω2 = (xf , a2 ) × (b1 , b2 ), Ω ⊂ R2 , Ω \γ = Ω1 ∪ Ω2 , Γ = ∂ Ω . Denote by Γi the part of the boundary of Ωi in common with the boundary of Ω , i.e., Γi = ∂ Ωi ∩ Γ , i = 1, 2. Let u1 , u2 and uf be the velocity and p1 , p2 and pf be the pressure defined in Ω1 , Ω2 and γ , respectively. Suppose that fluid flow in every subdomain is governed by Darcy’s law relating the gradient of pressure to the flow velocity together

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3431

Fig. 1. The domain of the fracture model.

with a conservation equation. The pressure is no longer assumed to have the continuity across the fracture. The flux exchange between the fracture and surrounding domains is added to the equation in the fracture as the coupling term. The coupled fracture model is detailedly described as follows. in Ωi , in Ωi ,

ui div ui

= =

gi

uf

=

−kf ,y d

=

gf + (u1 · n1 |γ +u2 · n2 |γ )

in γ ,

= = =

αf pf − (1 − ξ )ui+1 · ni+1 p¯ i p¯ f

in γ , on Γi , on ∂γ ,

∂ uf ∂y −ξ ui · ni + αf pi pi pf

−Ki ∇ pi ∂ pf ∂y

i = 1, 2, i = 1, 2,

in γ , (2.1)

i = 1, 2,

where the unit vectors ni , i = 1, 2 are the normal directions to γ outward from Ωi , and n1 = −n2 . We suppose that ¯ i ), 1 ≤ s, t ≤ 2, and the diagonal matrix Kf = diag(kf ,x , kf ,y ) the symmetric matrices Ki = (kist )2s,t =1 , i = 1, 2, kist ∈ C 1 (Ω respectively represent permeability tensors of the surrounding domains Ωi and the fracture γ , satisfying the following property. There exist constants kmax > kmin > 0 such that kmin ζ ′ ζ ≤ ζ ′ Ki (x)ζ ≤ kmax ζ ′ ζ , ∀ζ ∈ R2 , ∀x ∈ Ωi , i = 1, 2, kmin ≤ (kf ,y d), (kf ,x /d) ≤ kmax .

(2.2)

αf = 2kf ,x /d is the exchange coefficient, and the parameter d denotes the thickness of the fracture. For common application, the parameter ξ should be such that ξ ∈ (1/2, 1]. The index i varies in Z/2Z satisfying 2 + 1 = 1. For simplicity, we will consider the homogeneous Dirichlet boundary conditions in this paper, i.e., p¯ i = 0 and p¯ f = 0. Similar to [2], we rewrite the Robin boundary conditions in γ , the fifth equation in (2.1), in the following form. αf u2 · n2 |γ +u1 · n1 |γ = (p2 |γ +p1 |γ −2pf ), (2.3) 2ξ − 1 u2 · n2 |γ −u1 · n1 |γ = αf (p2 |γ −p1 |γ ). Summing and subtracting the above two equations, we can obtain ui · ni |γ =

αf (ξ pi |γ +(1 − ξ )pi+1 |γ −pf ), 2ξ − 1

i = 1, 2.

(2.4)

In this paper, we will use the standard definitions of the Sobolev spaces W m,s , for any non-negative integer m and number s ≥ 1, equipped with the standard semi-norm |·|W m,s and norm ∥ · ∥W m,s . When s = 2, this is the Hilbert space H m . We also introduce the associate inner product (·, ·)s , and we drop s when s = 2. 3. Discontinuous finite volume method In this section, we will construct and analysis the discontinuous finite volume scheme for the coupled fracture model (2.1). In the surrounding matrix Ωi , i = 1, 2, the pressure is solved by discontinuous piecewise linear elements. The continuity of the approximate functions is eliminated but instead a penalty term is added to enforce the connection between elements. In the fracture γ , the pressure is approximated by continuous piecewise linear functions. A numerical scheme only related to the pressure is formed, so we avoid solving the indefinite saddle-point problem resulting from the mixed method. When the approximate pressure is obtained, the velocity will be calculated in the discontinuous lowest-order Raviart–Thomas element space and piecewise linear element space in matrix and fracture element by element.

3432

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

Let Rhi be a triangulation of the sub-domain Ωi , i = 1, 2 with the longest diameter of elements hi , and Rh = ∪i=1,2 Rhi with h = maxi=1,2 hi . We also assume that Rh1 and Rh2 match on the fracture γ to form a partition for γ , Rhf : b1 = y1/2 < y3/2 < · · · < yn+1/2 = b2 .

Define yℓ =

yℓ−1/2 + yℓ+1/2 2

,

ℓ = 1, . . . , n.

Let σℓ , ℓ = 1, . . . , n be the line segment located in γ connecting yℓ−1/2 and yℓ+1/2 , and let σℓ+1/2 , ℓ = 0, . . . , n be the line segment connecting yℓ and yℓ+1 with y0 = y1/2 and yn+1 = yn+1/2 . We assume that Rh is quasi-uniform and the following hypotheses are true. There exist two positive constants α2 > α1 > 0 such that

α1 h2 ≤ |K | ≤ α2 h2 , and α1 h ≤ |σ | ≤ α2 h, where |K | and |σ | denote the area of any triangle K and the length of any edge σ . In Ωi , i = 1, 2, we define the dual partition Ti h for the test function space. For every triangle K ∈ Rhi , connect the barycenter with three vertices, and we divide K into three small triangles, and every small triangle is regarded as a control volume. The dual partition Ti h consists of all these small triangles, and T h = T1h ∪ T2h . Besides, denote by Tf h the dual partition of Rhf in the fracture γ , Tf h : b1 < y1 < · · · < yn < b2 .

Denote by EK the set of three edges of the triangle K . Let Ei◦ denote the sets of all interior edges of triangle K ∈ Rhi and let EiΓ denote all edges located on the boundary Γi , and Ei = Ei◦ ∪ EiΓ . We define Eγ to be the union of the edges σℓ , ℓ = 1, . . . , n lying in the fracture γ . For any interior edge σ ∈ Ei◦ , shared by two triangles K1 and K2 , we use the same definitions of the average {·} and jump [·] as those in [27]. For any scalar qi and any vector vi , 1

{qi } =

2 1

{vi } =

2

(qi |∂ K1 +qi |∂ K2 ),

[qi ] = qi |∂ K1 nK1 + qi |∂ K2 nK2 ,

(vi |∂ K1 +vi |∂ K2 ),

[vi ] = vi |∂ K1 ·nK1 + vi |∂ K2 ·nK2 ,

where nK1 and nK2 are unit normal vectors on σ pointing exterior to K1 and K2 . For any edge σ ∈ EiΓ or σ ∈ Eγ , we define

{qi } = qi ,

[vi ] = vi · nK ,

where nK denotes the unit outward normal vector of triangle K containing σ . The following identity can be obtained by a direct computation, shown in [23],

∑ ∫ K ∈Rhi

∂K

vi · nqi ds =

∑∫ σ ∈Eio

σ

∑ ∫

[qi ] · {vi }ds +

σ ∈Ei ∪Eγ

σ

{qi }[vi ]ds.

(3.1)

Based on triangulations, we define the finite element spaces Wih , i = 1, 2 and Wfh for the pressure in Ωi and γ . Wih = {wi ∈ L2 (Ωi ) : wi |K ∈ P1 (K ), ∀K ∈ Rhi },

i = 1, 2,

Wfh = {wf ∈ H01 (γ ) : wf |σℓ ∈ P1 (σℓ ), ℓ = 1, . . . , n}, where Ps (D) is the space containing all polynomials with degree less than or equal to s in domain D. The associated discrete H 1 semi-norm in Ωi is denoted by

|wi |21,h,Ωi =

∑

|wi |2H 1 (K ) .

K ∈Rih

We define the global space for p, W h = w = (w1 , w2 , wf ) ∈ W1h × W2h × Wfh ,

{

}

equipped with the following semi-norm

|w|21,h =

∑ i=1,2

|wi |21,h,Ωi + |wf |2H 1 (γ ) .

(3.2)

For the velocity, we introduce the discontinuous lowest-order Raviart–Thomas element space Vih in matrix Ωi , i = 1, 2, Vih = {vi ∈ L2 (Ωi ) : vi |K ∈ RT0 (K ), ∀K ∈ Rhi },

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3433

where RT0 (K ) = {φ = (φ 1 , φ 2 ) : φ 1 = a + bx, φ 2 = c + by, in K }, with the discontinuous H(div ) norm,

|vi |2H(div),h,Ωi =

∑

∥div vi ∥2L2 (K ) .

K ∈Rhi

We introduce the discontinuous piecewise linear element space Vfh in the fracture γ , Vfh = {vf ∈ L2 (γ ) : vf |σℓ ∈ P1 (σℓ ), ℓ = 1, . . . , n}, with the discrete H 1 (γ ) norm

|vf |21,h,γ =

n   ∑  ∂vf 2 .   ∂ y L2 (σℓ )

ℓ=1

The global approximate space for the velocity Vh = {v = (v1 , v2 , vf ) ∈ V1h × V2h × Vfh }, with the following discrete H(div ) norm,

|v|2H(div),h =

∑ i=1,2

|vi |2H(div),h,Ωi + |vf |21,h,γ .

We define the finite element spaces Q h = Q1h × Q2h × Qfh for test functions. Qih = {qi ∈ L2 (Ωi ) : qi |T ∈ P0 (T ), ∀T ∈ Ti h }, Qfh

i = 1, 2,

= {qf ∈ L (γ ) : qf |σℓ+1/2 ∈ P0 (σℓ+1/2 ), ℓ = 0, . . . , n}. 2

Let Wi (h) = Wih + H 2 (Ωi ) ∩ H01 (Ωi ), i = 1, 2, and similar to [23,24], there exist mappings Πi : Wi (h) ↦ → Qih

Πi wi |T =

∫

1

|σ |

σ

wi |T ds,

T ∈ Ti h ,

where σ is the edge sharing by triangles T ∈ Ti h and K ∈ Rhi consisting of T . Let Wf (h) = Wfh + H 2 (γ ) ∩ H01 (γ ) and we define a mapping Πf : Wf (h) ↦ → Qfh such that

Πf wf |σℓ+1/2 = wf (yℓ+1/2 ),

ℓ = 0, . . . , n.

Denote by W (h) = W1 (h) × W2 (h) × Wf (h), and let Π = (Π1 , Π2 , Πf ) be a mapping from W (h) to Q h , which has the following set of properties. For any K ∈ Rhi , let σ be an edge of K , and σℓ , σℓ+1/2 are defined as before in fracture γ , we have

Lemma 3.1.

∫ (1)

(wi − Πi wi )dx = 0,

∀wi ∈ Wih ,

(wi − Πi wi )ds = 0,

∀wi ∈ Wih ,

K

∫ (2) σ

(3) (4)

∥wih − Πi wi ∥L2 (K ) ≤ Ch|wi |H 1 (K ) , ∀wi ∈ Wih , ∫ [wi − Πi wi ]ds = 0, ∀wi ∈ H 2 (Rhi ),

(3.3)

σ

(5)

if [wi ] = 0, then [Πi wi ] = 0,

∫ (6) σℓ

(7)

(wf − Πf wf )ds = 0,

∀wi ∈ H 2 (Rhi ),

∀wf ∈ Wfh ,

∥wf − Πf wf ∥L2 (σℓ ) ≤ Ch|wf |H 1 (σℓ ) ,

∀wf ∈ Wfh ,

where H 2 (Rhi ) := {wi ∈ L2 (Ωi ) : wi |K ∈ H 2 (K ), ∀K ∈ Rhi }. Proof. Properties (1)–(5) have been proved in [25,28]. For the properties of Πf , (6)–(7), noting that wf ∈ Wfh is a piecewise linear function and yℓ is the midpoint of edge σℓ ,

∫ σℓ

wf ds =

|σℓ | ( 2

) |σℓ | ( ) wf (yℓ+1/2 ) + wf (yℓ−1/2 ) = Πf wf |σℓ+1/2 +Πf wf |σℓ−1/2 = 2

∫ σℓ

Πf wf ds.

3434

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

Fig. 2. The triangulation and the dual partition.

∫ σℓ

yℓ

∫

(wf − Πf wf )2 ds =

( )2 wf − wf (yℓ−1/2 ) ds +

yℓ+1/2

∫

)2 wf − wf (yℓ+1/2 ) ds ≤ Ch2 |wf |2H 1 (σ ) . □ ℓ

yℓ

yℓ−1/2

(

We then derive the formula of the DFVM for the coupled fracture model. Multiply the second equation in (2.1) by a test function qi ∈ Qih , integrate in control volumes, and the area integrals will disappear using integrating by part.

∫

∫ ∂T

gi qi dx,

ui · nqi ds = T

where n is the unit outward normal vector on ∂ T . By the first equation of (2.1), we get

−

∑∫ ∂T

T ∈Ti h

∑∫

Ki ∇ pi · nqi ds =

T ∈Ti h

gi qi dx.

(3.4)

T

As shown in Fig. 2, Tj ∈ Ti h (j = 1, 2, 3) are three triangles in every triangle K ∈ Rhi . We can get the same equation as that in [23].

∑∫ T ∈Ti

=

h

∂T

K ∈Rhi j=1

3 ∫ ∑ ∑ K ∈Rhi j=1

3 ∫ ∑ ∑

Ki ∇ pi · nqi ds =

Ki ∇ pi · nqi ds + Aj+1 CAj

∂ Tj

Ki ∇ pi · nqi ds (3.5)

∑ ∫ K ∈Rhi

∂K

Ki ∇ pi · nqi ds,

where A4 = A1 . Using the identity (3.1) and the Robin boundary condition (2.4) on γ , and by the fact [Ki ∇ pi ] = 0 in any interior edge, we have

∑ ∫ ∂K

K ∈Rhi

Ki ∇ pi · nqi ds =

∑∫ σ ∈Ei

σ

[qi ]{Ki ∇ pi }ds −

αf

∑ ∫ σℓ ∈Eγ

σ

2ξ − 1

(ξ pi |γ +(1 − ξ )pi+1 |γ −pf )qi ds.

(3.6)

Hence, we obtain the following equation in Ωi ,

−

3 ∫ ∑ ∑ K ∈Rhi j=1

+

∑ ∫ σℓ ∈Eγ

σℓ

Ki ∇ pi · nqi ds − Aj+1 CAj

∑∫ σ ∈Ei

σ

[qi ]{Ki ∇ pi }ds

αf (ξ pi |γ +(1 − ξ )pi+1 |γ −pf )qi ds = (gi , qi )Ωi . 2ξ − 1

(3.7)

In the fracture γ , we multiply the fourth equation in (2.1) by the test function qf ∈ Qfh and denote by qf ,ℓ+1/2 the value of qf on the edge σℓ+1/2 ,

−

n−1 (( ∑

kf ,y d

ℓ=1

∂ pf ∂y

)⏐ ⏐ ⏐

yℓ+1

( ) ) ∫ ∫ ∂ pf ⏐⏐ αf − k f ,y d qf ,ℓ+1/2 = gf qf ds + (p1 |γ +p2 |γ −2pf )qf ds, ⏐ yℓ ∂y γ γ 2ξ − 1

where φ|yℓ = φ (yℓ ) for any function φ .

(3.8)

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3435

For simplicity, we give the following definition Dℓ qf =

qf (yℓ+1/2 ) − qf (yℓ−1/2 )

|σℓ |

.

By direct calculation, (3.8) can be rewritten as n ( ∑

ℓ=1

∂ pf kf ,y d ∂y

)⏐ ∫ ⏐ ⏐ Dℓ qf |σℓ | − yℓ

γ

αf (p1 |γ +p2 |γ −2pf )qf ds = 2ξ − 1

∫ γ

gf qf ds.

(3.9)

Keeping (3.7)–(3.9) in mind, the DFVM for the coupled fracture model (2.1) reads as follows.

∑

A(ph , w ) =

(gi , Πi wi )Ωi + (gf , Πf wf )γ ,

∀w = (w1 , w2 , wf ) ∈ W h ,

(3.10)

i=1,2

with A(ph , w ) = −

3 ∫ ∑ ∑ ∑ i=1,2 K ∈Rh j=1 i

∑ ∑∫

−

i=1,2 σ ∈Ei

σ

Ki ∇ phi · nΠi wi ds − Aj+1 CAj

∑ ∑∫ i=1,2 σ ∈Ei

[Πi phi ]{Ki ∇wi }ds + α

∑∑

σ

[Πi wi ]{Ki ∇ phi }ds

[Πi phi ]σ [Πi wi ]σ

i=1,2 σ ∈Ei

(3.11)

αf + (ξ phi + (1 − ξ )phi+1 − phf )Πi wi ds 2ξ − 1 γ i=1,2 ) ( ∫ n ∑ ∂ phf ⏐⏐ αf D Π w |σ | − (ph1 + ph2 − 2phf )Πf wf ds, + kf ,y d ⏐ ℓ f f ℓ yℓ ∂y 2 ξ − 1 γ

∑∫

ℓ=1

where α is a real number determined in the following analysis. From the numerical scheme (3.10)–(3.11), the approximate pressure ph = (ph1 , ph2 , phf ) is obtained, and then we use the method proposed in [29] to calculate the velocity uhi and uhf in the discontinuous lowest-order Raviart–Thomas element space and piecewise linear element space. In every triangle K ∈ Rhi , similar to equation (2.5) in [29], uhi is expressed as follows. uhi |K = −Ai,K ∇ phi |K +gi,K PK (x), ∀x ∈ K , i = 1, 2,

with PK (x) =

where (xB , yB ) is the coordinates of the barycenter of K , and Ai,K = tensor Ki and the right side function gi over K . On every edge σℓ ∈

1

1

(

2

∫

K Rhf , uhf

|K |

x − xB y − yB

)

,

Ki dx and gi,K =

(3.12) 1

|K |

∫ K

gi dx are averages of the

is computed by

) ( ∂ phf ⏐⏐ αf ,ℓ h h h (3.13) (p + p2,ℓ − 2pf ,ℓ ) (y − yℓ ), | = −Af ,ℓ ⏐ + gf ,ℓ − ∂ y σℓ 2ξ − 1 1,ℓ ∫ where Af ,ℓ = |σ1 | σ (kf ,y d)ds is the average of the product kf ,y d over the edge σℓ , similar to αf ,ℓ , phi,ℓ , i = 1, 2, phf,ℓ . ℓ ℓ For the exact solution p of the coupled fracture model, [Πi pi ] = 0 by Lemma 3.1, and from (3.7) and (3.9), it is easy to see that the bilinear form A(·, ·) is consistent with the fracture model (2.1), i.e., ∑ A(p, w ) = (gi , Πi wi )Ωi + (gf , Πf wf )γ , ∀w = (w1 , w2 , wf ) ∈ W h . (3.14) uhf σℓ

i=1,2

In order to analyze the proposed numerical scheme, we define the following |||·||| norm over W (h),

|||p|||2 =

∑ i=1,2

|pi |21,h,Ωi + |pf |2H 1 (γ ) +

∑∑ i=1,2 σ ∈Ei

[Πi pi ]2σ +

∑ i=1,2

∥pi − pf ∥2L2 (γ ) .

(3.15)

4. Error estimate In this section, we will prove the stability and convergence of the proposed DFVM for the coupled fracture model. To this end, We first list some results used in the following analysis. The first lemma had been proved in [23].

3436

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

For any p = (p1 , p2 , pf ), w = (w1 , w2 , wf ) ∈ W (h), we have

Lemma 4.1.

3 ∫ ∑ ∑

Ki ∇ pi · nΠi wi ds = (Ki ∇h pi , ∇h wi )Ωi Aj+1 CAj

K ∈Rhi j=1

∑ ∫

+

K ∈Rhi

∂K

(4.1)

(Πi wi − wi )Ki ∇ pi · nds +

∑

(∇ · Ki ∇ pi , wi − Πi wi )K .

K ∈Rhi

Furthermore, if p, w ∈ W h , we have 3 ∫ ∑ ∑ K ∈Rhi j=1

Aj+1 CAj

Ki ∇ pi · nΠi wi ds ≥ (Ki ∇h pi , ∇h wi )Ωi − C1 h|pi |1,h,Ωi |wi |1,h,Ωi ,

where ∇h pi |K = ∇ pi , ∀K ∈ Rhi . If α is large enough and h is small enough, there exists a constant C independent of h such that

Lemma 4.2.

A(p, p) ≥ C |||p|||,

∀p ∈ W h .

Proof. We first introduce the trace inequality, ∀w ∈ H 2 (K ),

( ) ∥w∥2L2 (e) ≤ C h−1 ∥w∥2L2 (K ) + h|w|2H 1 (K ) .

(4.2)

Then for any p ∈ W h , we get

∑ ∑∫ i=1,2 σ ∈Ei

σ

[Πi pi ]{Ki ∇h pi }ds

⎞ 12 ⎛ ⎞ 21 ∫ ( ) ∑ ∑ ⎟ ⎜ [Πi pi ]2 ds⎠ h−1 |pi |2H 1 (K )2 + h|pi |2H 2 (K ) ⎠ ⎝ ∥Ki ∥1,∞ ⎝ ⎛

∑

≤ C

i=1,2

( ∑

≤ C

σ ∈Ei

K ∈Rhi

σ

(4.3)

⎞ 21 ) 12 ⎛ ∑∑ ⎝ |pi |21,h,Ωi [Πi pi ]2σ ⎠ .

i=1,2

i=1,2 σ ∈Ei

For the sixth term in the bilinear form A(p, p), since pf ∈ Wfh is a piecewise linear function, it is easy to check that Dℓ Πf pf =

Πf pf (yℓ+1/2 ) − Πf pf (yℓ−1/2 ) pf (yℓ+1/2 ) − pf (yℓ−1/2 ) ∂ pf ⏐⏐ = = ⏐ . |σℓ | |σℓ | ∂ y yℓ

Hence, using the assumption of the coefficient kf ,y d ≥ kmin , we have n ( ∑

ℓ=1

∂ pf k f ,y d ∂y

)⏐ )2 n ( ∑ ∂ pf ⏐⏐ ⏐ |σℓ | = kmin |pf |2H 1 (γ ) . ⏐ Dℓ Πf pf |σℓ | ≥ kmin ⏐ yℓ ∂ y yℓ

(4.4)

ℓ=1

Combining the fifth term with the seventh term in A(p, p), we obtain

∫ αf αf (ξ pi + (1 − ξ )pi+1 − pf )Πi pi ds − (p1 + p2 − 2pf )Πf pf ds 2 ξ − 1 2 ξ −1 γ i=1,2 ∫γ ∑ ) αf ( = ξ (pi − pf ) + (1 − ξ )(pi+1 − pf ) (Πi pi − Πf pf )ds 2ξ − 1 i=1,2 ∫γ ∑ ) αf ( = ξ (pi − pf ) + (1 − ξ )(pi+1 − pf ) (pi − pf )ds 2 ξ − 1 i=1,2 γ ∫ ∑ )( ) αf ( + ξ (pi − pf ) + (1 − ξ )(pi+1 − pf ) (Πi pi − pi ) − (Πf pf − pf ) ds γ 2ξ − 1 ∑∫

i=1,2

:= S1 + S2 .

(4.5)

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449 2kf ,x

By the assumption of αf =

≥ 2kmin , it follows that

d

) αf ξ (pi − pf ) + (1 − ξ )(pi+1 − pf ) (pi − pf )ds 2ξ − 1 i=1,2 ∫γ ∫ ∑ αf αf ξ (pi − pf )2 ds + 2(1 − ξ ) (p1 − pf )(p2 − pf )ds = 2ξ − 1 γ 2ξ − 1 i=1,2 γ ∑ 2kmin ≥ (2ξ − 1) ∥pi − pf ∥2L2 (γ ) 2ξ − 1 i=1,2 ∑ = 2kmin ∥pi − pf ∥2L2 (γ ) . ∑∫

S1 =

3437

(

(4.6)

i=1,2

Using the Cauchy–Schwarz inequality, we have

αf

∑∫

S2 =

γ

i=1,2

( ≤C

2ξ − 1

)( ) ξ (pi − pf ) + (1 − ξ )(pi+1 − pf ) (Πi pi − pi ) − (Πf pf − pf ) ds )1/2 (

∑∫

2

γ

i=1,2

(

(pi − pf ) ds

∑∫ i=1,2

(Πi pi − pi ) ds + 2

γ

)1/2

∫

(Πf pf − pf ) ds 2

γ

.

By the trace inequality (4.2), it is easily obtained that

⎞

⎛ ∫

(Πi pi − pi )2 ds ≤ C ⎝h−1

⎜

γ

∑

∑

∥Πi pi − pi ∥2L2 (K ) + h

K ∈Rhi

∑ ⎟ |pi |2H 1 (K ) . |Πi pi − pi |2H 1 (K ) ⎠ ≤ Ch

(4.7)

K ∈Rhi

K ∈Rhi

According to (4.7) and the inequality (7) in Lemma 3.1, we get

)1/2 (

( 1/2

∑

S2 ≤ Ch

pi 21,h,Ωi

| |

+| |

i=1,2

)1/2 ∑

pf 2H 1 (γ )

2 L2 (γ )

∥pi − pf ∥

i=1,2

.

(4.8)

Taking (4.3)–(4.4), (4.6) and (4.8) back into the expression of A(p, p) and by Lemma 4.1, if α is large enough, we have

∑

A(p, p) ≥ kmin

|pi |21,h,Ωi + kmin |pf |2H 1 (γ ) + α

i=1,2

∑∑

[Πi pi ]2σ

i=1,2 σ∈ Ei

⎞ 12

⎛ −Ch

∑

|pi |21,h,Ωi − C

i=1,2

∑

|pi |1,Ωi ,h ⎝

∑ i=1,2

≥ C⎝

i=1,2

|pi |21,h,Ωi + |pf |21,h,γ

i=1,2

(4.9)

∥pi − pf ∥2L2 (γ )

⎞ ∑∑ ∑ + [Πi pi ]2σ + ∥pi − pf ∥2L2 (γ ) ⎠ i=1,2 σ∈ Ei

i=1,2

( − C1 h

∥pi − pf ∥2L2 (γ )

)1/2 ∑

|pi |21,h,Ωi + |pf |2H 1 (γ )

∑ i=1,2

)1/2 (

⎛ ∑

[Πi pi ]2σ ⎠ + 2kmin

σ ∈Ei

i=1,2

( −Ch1/2

∑

) ∑

pi 21,h,Ωi

+h

| |

1/2

i=1,2

∑

pi 21,h,Ωi

| |

i=1,2

+h

1/2

pf 2H 1 (γ )

| |

.

Hence, if h is small enough, the desired result holds. □ Lemma 4.3.

For any p, w ∈ W (h), we have

⎞ 12 ⎞ ⎜ ⎜∑ ∑ 2 2 ⎟ ⎟ A(p, w ) ≤ C ⎜ h |pi |H 2 (K ) + h2 |pf |2H 2 (γ ) ⎠ ⎟ ⎠ ⎝|||p|||+ ⎝ ⎛

⎛

i=1,2 K ∈Rh i

⎞ 12 ⎞ ∑ ∑ ⎜ ⎟ ⎟ ⎜|||w|||+ ⎜ h2 |wi |2H 2 (K ) + h2 |wf |2H 2 (γ ) ⎠ ⎟ ⎝ ⎝ ⎠ ⎛

⎛

i=1,2 K ∈Rh i

(4.10)

3438

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

Furthermore, if p, w ∈ W h , A(p, w ) ≤ C |||p||||||w|||.

(4.11)

Proof. According to Lemma 4.1, and using the trace inequality (4.2), it yields

⏐ ⏐ 3 ∫ ⏐∑ ∑ ∑ ⏐ ⏐ ⏐ Ki ∇ pi · nΠi wi ds⏐ ⏐ ⏐ ⏐ Aj+1 CAj h i=1,2 K ∈R

j=1

i

⏐ ⏐ ⏐∑ ⏐ ⏐∑ ∑ ∫ ⏐ ⏐ ⏐ ⏐ ⏐ ≤ ⏐ (Ki ∇h pi , ∇h wi )Ωi ⏐ + ⏐ (Πi wi − wi )Ki ∇ pi · nds⏐ ⏐ ⏐ ∂ K h i=1,2

i=1,2 K ∈R

i

⏐ ⏐ ⏐∑ ∑ ⏐ ⏐ ⏐ +⏐ (∇ · Ki ∇ pi , wi − Πi wi )K ⏐ ⏐ ⏐ h i=1,2 K ∈R

≤ C

∑

i

∑ ∑

|pi |1,h,Ωi |wi |1,h,Ωi + C

i=1,2

+C

(4.12)

)

(

h |pi |H 1 (K ) + |pi |H 2 (K ) |wi |H 1 (K )

i=1,2 K ∈Rh i

∑ ∑ (

h−1 ∥wi − Πi wi ∥2L2 (K ) + h|wi − Πi wi |2H 1 (K )

) 12 ) 21 ( · h−1 |pi |2H 1 (K ) + h|pi |2H 2 (K )

i=1,2 K ∈Rh i

⎛

⎞ 21 ⎞ ∑ ⎜∑ ∑ 2 2 ⎟ ⎟ |wi |1,h,Ωi . +⎝ h |pi |H 2 (K ) ⎠ ⎟ ⎠ ⎛

⎜∑ |pi |1,h,Ωi ≤ C⎜ ⎝ i=1,2

i=1,2 K ∈Rh

i=1,2

i

Using the Cauchy–Schwarz inequality and the trace inequality (4.2), we have

−

∑ ∑∫ i=1,2 σ ∈Ei

σ

{Ki ∇ pi }[Πi wi ]ds −

∑ ∑∫ i=1,2 σ ∈Ei

σ

{Ki ∇wi }[Πi pi ]ds + α

∑∑

[Πi pi ]σ [Πi wi ]σ

i=1,2 σ ∈Ei

⎞1/2 ⎛ ⎞1/2 ) ∑∑ ⎜∑ ∑ ⎟ ≤ C⎝ h h−1 |pi |2H 1 (K ) + h|pi |2H 2 (K ) ⎠ ⎝ [Πi wi ]2σ ⎠ ⎛

(

i=1,2 σ ∈Ei

i=1,2 K ∈Rh i

⎞1/2 ⎛ ⎞1/2 ( ) ∑ ∑ ⎜ ⎟ ⎝∑ ∑ −1 2 2 2⎠ +C ⎝ h h |wi |H 1 (K ) + h|wi |H 2 (K ) ⎠ [Πi pi ]σ ⎛

i=1,2 σ ∈Ei

i=1,2 K ∈Rh i

⎞1/2 ⎛

⎛ +α⎝

∑∑

[Πi pi ]2σ ⎠

⎞1/2 ∑∑

⎝

i=1,2 σ ∈Ei

(4.13)

[Πi wi ]2σ ⎠

i=1,2 σ ∈Ei

⎛ ⎞ 12 ⎞ ⎞ 12 ⎞ ⎛ ∑⎜∑ ∑⎜∑ ⎜ ⎟ ⎟ ⎟ ⎟⎜ ⎜ h2 |wi |2H 2 (K ) ⎠ ⎟ ≤ C⎜ h2 |pi |2H 2 (K ) ⎠ ⎟ ⎝ ⎝ ⎠ ⎝|||w|||+ ⎠. ⎝|||p|||+ ⎛

⎛

i=1,2

i=1,2

K ∈Rih

K ∈Rih

For the term in the fracture γ , we have n ( ∑

)⏐ ) n ∫ ( ∑ ∂ pf ⏐⏐ wf ,ℓ+1/2 − wf ,ℓ−1/2 ⏐ D Π w |σ | = k d ds ⏐ ℓ f f ℓ ⏐ f ,y yℓ yℓ ∂y σℓ ℓ=1 σℓ ℓ=1 ( )1/2 ( )1/2 n n ∑ ∑ 2 2 2 2 2 2 ≤ C |pf |H 1 (γ ) + h |pf |H 2 (σ ) |wf |H 1 (γ ) + h |wf |H 2 (σ ) . kf ,y d

∂ pf ∂y

ℓ=1

ℓ

ℓ=1

ℓ

(4.14)

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3439

For the coupled items on interface, similar to (4.8), it is obtained that

∫ αf αf (ξ pi + (1 − ξ )pi+1 − pf )Πi wi ds − (p1 + p2 − 2pf )Πf wf ds 2ξ − 1 2ξ − 1 γ γ i=1,2 ∫ ∑ ) αf ( = ξ (pi − pf ) + (1 − ξ )(pi+1 − pf ) (wi − wf )ds 2ξ − 1 i=1,2 γ ∫ ∑ )( ) αf ( + ξ (pi − pf ) + (1 − ξ )(pi+1 − pf ) (Πi wi − wi ) − (Πf wf − wf ) ds 2 ξ − 1 γ (i=1,2 )1/2 ( )1/2 ∑ ∑ ≤ C ∥pi − pf ∥2L2 (γ ) ∥wi − wf ∥2L2 (γ ) ∑∫

i=1,2

( ∑

+C

∥pi − pf ∥2L2 (γ )

i=1,2

=1,2 )1/2 i( ∑

i=1,2

(4.15)

)1/2 |wi |21,h,Ωi + |wf |2H 1 (γ )

.

Thanks to (4.12)–(4.15) and the definition of the bilinear form A(·, ·), we get (4.10). If p, w ∈ W h , we have |pi |H 2 (K ) = 0 and |wi |H 2 (K ) = 0, ∀K ∈ Rhi , so the continuity over W h is obtained. □ Lemma 4.4.

There exists a constant C with independence of h such that

∀w ∈ W h .

∥w∥L2 (Ω ) ≤ C |||w|||,

Proof. For any w = (w1 , w2 , wf ) ∈ W h , by Lemma 3.1 in [23],

⎛

⎞

∥wi ∥2L2 (Ω ) ≤ ⎝|wi |21,h,Ωi + i

∑

[Πi wi ]2σ ⎠ ,

i = 1, 2.

(4.16)

σ ∈Ei

By using the Poincaré inequality,

∥wf ∥2L2 (γ ) ≤ C |wf |2H 1 (γ ) .

(4.17)

From (4.16)–(4.17), the result holds. Let w = (w , w , w following property. I 2

I f)

|wi − wiI |H s (K ) |wf − wfI |H s (σ )

≤ ≤

I 1

I

ℓ

□

∈ W be the usual continuous piecewise interpolation of any function w, which satisfies the h

Ch2−s |wi |H 2 (K ) , 2−s

Ch

∀K ∈ Rih , i = 1, 2, s = 0, 1, |wf |H 2 (σℓ ) , ℓ = 1, . . . , n, s = 0, 1.

(4.18)

By the definition of the norm |||·|||,

|||p − pI |||

∑

=

i=1,2

+

2

2

|pi − pIi |1,h,Ωi + |pf − pIf |H 1 (γ ) +

∑ i=1,2

∑∑

[Πi (pi − pIi )]2σ

i=1,2 σ ∈Ei

(4.19)

∥(pi − pIi ) − (pf − pIf )∥2L2 (γ ) .

Using the trace inequality and (4.18), for any edge σ ∈ Ei ,

Πi (pi − pIi )|2σ ≤ C

1

∫

|σ |

(

σ

2

(pi − pIi )2 ds ≤ C h−2 ∥pi − pIi ∥2L2 (K ) + |pi − pIi |H 1 (K )

)

≤ Ch2 |pi |2H 2 (K ) .

(4.20)

Similarly, for any σℓ , ℓ = 1, . . . , n,

( ) 2 ∥pi − pIi ∥2L2 (σ ) ≤ C h−1 ∥pi − pIi ∥2L2 (K ) + h|pi − pIi |H 1 (K ) ≤ Ch3 |pi |2H 2 (K ) ,

(4.21)

∥pf − pIf ∥L2 (γ ) ≤ Ch2 |pf |H 2 (γ ) .

(4.22)

ℓ

and

Hence, we get the following property about the interpolation pI ,

( I

2

|||p − p ||| ≤ Ch

2

) ∑ i=1,2

pi 2H 2 (Ω ) i

| |

pf 2H 2 (γ )

+| |

.

(4.23)

3440

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

Theorem 4.1. Suppose that p = (p1 , p2 , pf ) is the exact solution for the fracture model (2.1) and ph = (ph1 , ph2 , phf ) is the approximation. Then there is a positive constant C independent of mesh size h such that

( h 2

|||p − p ||| ≤ Ch

) ∑

2

pi 2H 2 (Ω ) i

| |

i=1,2

.

pf 2H 2 (γ )

+| |

Proof. Comparing (3.10) with (3.14), we get the following error equation, A(p − ph , w ) = 0,

∀w = (w1 , w2 , wf ) ∈ W h .

(4.24)

By using Lemma 4.2, we get

|||pI − ph |||2 ≤ CA(pI − ph , pI − ph ) = C (A(pI − p, pI − ph ) + A(p − ph , pI − ph )).

(4.25)

From Lemma 4.3, it is obtained that A(pI − p, pI − ph )

⎞ 21 ⎞ ⎜ I 2 2 ⎟ ⎟ I ⎜∑ ∑ 2 I h h |pi − pi |H 2 (K ) + h2 |pIf − pf |H 2 (γ ) ⎠ ⎟ C⎜ ⎠ |||p − p ||| ⎝|||p − p|||+ ⎝ ⎛

≤

⎛

i=1,2 K ∈Ri

(4.26)

h

)1/2

( ∑

Ch|||pI − ph |||

≤

i=1,2

.

|pi |2H 2 (Ω ) + |pf |2H 2 (γ ) i

Using (4.24), we have A(p − ph , pI − ph ) = 0.

(4.27)

Taking (4.26)–(4.27) back into (4.25), we have

)1/2

( I

h

|||p − p |||≤ Ch

∑

pi 2H 2 (Ω ) i

| |

i=1,2

pf 2H 2 (γ )

+| |

.

(4.28)

From (4.23) and (4.28), by using the triangle inequality, the result holds.

□

Next, we turn our attention to the L2 error estimate of the pressure p for the proposed DFVM. Consider the following dual problem of coupled fracture model: Find φi ∈ H 2 (Ωi ) ∩ H01 (Ωi ), i = 1, 2 and φf ∈ H 2 (γ ) ∩ H01 (γ ) satisfying

−div(Ki ∇φi ) = pi − phi , in Ωi , i = 1, 2 ( ) ) ( ∂ ∂φf h − kf ,y d = pf − pf − K1 ∇φ1 · n1 |γ +K2 ∇φ2 · n2 |γ , in γ ∂y ∂y ξ Ki ∇φi · ni |γ +αf φi = αf φf + (1 − ξ )Ki+1 ∇φi+1 · ni+1 |γ ,

(4.29)

in γ

φi = 0,

on Γi , i = 1, 2

φf = 0 ,

on ∂γ .

Similar to (2.4), we have the Robin boundary condition in γ

− Ki ∇φi · ni |γ =

αf

2ξ − 1

(ξ φi |γ +(1 − ξ )φi+1 |γ −φf ),

i = 1, 2.

(4.30)

We assume that the solution φ = (φ1 , φ2 , φf ) has the following regularity property,

( ∑

∥φi ∥H 2 (Ωi ) + ∥φf ∥H 2 (γ ) ≤ C

i=1,2

) ∑

∥pi −

phi L2 (Ωi )

∥

+ ∥ pf −

phf L2 (γ )

∥

.

(4.31)

i=1,2

Theorem 4.2. Let (p1 , p2 , pf ) and (ph1 , ph2 , phf ) be respectively solutions of (2.1) and (3.10), and we assume that pi |γ ∈ H 1 (γ ), i = 1, 2. Then there exists a constant C with independence of h such that

∑

∥pi − phi ∥L2 (Ωi ) + ∥pf − phf ∥L2 (γ ) i=1,2( ) ∑ ∑ ∑ 2 ≤ Ch ∥pi ∥H 2 (Ωi ) + ∥pf ∥H 2 (γ ) + ∥ p i ∥ H 1 (γ ) + ∥gi ∥H 1 (Ωi ) + ∥gf ∥H 1 (γ ) . i=1,2

i=1,2

i=1,2

(4.32)

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3441

Proof. In surrounding domains Ωi , i = 1, 2, from the first equation in (4.29) and (4.30), we get

) ( ∥pi − phi ∥2L2 (Ω ) = pi − phi , −div(Ki ∇φi ) Ω i i ∑ ∫ ) ( h = ∇h (pi − pi ), Ki ∇φi Ω − i

K ∈Rhi

∂K

Ki ∇φi · n(pi − phi )ds

∑∫ ) ( {Ki ∇φi }[pi − phi ]ds = ∇h (pi − phi ), Ki ∇φi Ω − i σ σ ∈Ei ∫ ) αf ( + ξ φi |γ +(1 − ξ )φi+1 |γ −φf (pi − phi )ds, 2 ξ − 1 γ

(4.33)

where using identity (3.1) and by the fact [Ki ∇φi ] = 0, the last equation holds. In the fracture γ , using the second equation in (4.29) and (4.30), it yields

( ) ) ( ) ∂ ∂φf − kf ,y d , pf − phf + K1 ∇φ1 · n1 |γ +K2 ∇φ2 · n2 |γ , pf − phf γ ∂y γ ( ∂y ) ∫ h ) ∂φf ∂ (pf − pf ) αf ( − = kf ,y d , φ1 |γ +φ2 |γ −2φf (pf − phf )ds. ∂y ∂y γ 2ξ − 1

∥pf − phf ∥2L2 (γ ) =

(

(4.34)

γ

Combining the last term in (4.33) with the last term in (4.34), we have

∫ ) αf ( αf h ξ φi |γ +(1 − ξ )φi+1 |γ −φf (pi − pi )ds − (φ1 |γ +φ2 |γ −2φf )(pf − phf )ds 2 ξ − 1 2 ξ −1 γ γ i=1,2 ∫ ∑ ) αf ( ξ (pi − phi ) + (1 − ξ )(pi+1 − phi+1 ) − (pf − phf ) φi ds = 2ξ − 1 i=1∫ ,2 γ ) αf ( − (p1 − ph1 ) + (p2 − ph2 ) − 2(pf − phf ) φf ds. γ 2ξ − 1 ∑∫

(4.35)

From (4.33)–(4.35), it is obtained that

∑ i=1,2

∥pi − phi ∥2L2 (Ω ) + ∥pf − phf ∥2L2 (γ ) i

( ) h ∑ ∑∫ ∑( ) ∂φf ∂ (pf − pf ) h h {Ki ∇φi }[pi − pi ]ds + kf ,y d = ∇h (pi − pi ), Ki ∇h φi Ω − , i ∂y ∂y σ σ ∈ E i = 1 , 2 i=1,2 γ i ∑∫ ) αf ( ξ (pi − phi ) + (1 − ξ )(pi+1 − phi+1 ) − (pf − phf ) φi ds + 2ξ − 1 i∫ =1,2 γ ) αf ( − (p1 − ph1 ) + (p2 − ph2 ) − 2(pf − phf ) φf ds. 2 ξ − 1 γ

(4.36)

Let φ I = (φ1I , φ2I , φfI ) ∈ W h be the usual continuous piecewise linear interpolation of φ , and we have A(p − ph , φ I ) = 0 by (4.24), i.e.

−

3 ∫ ∑ ∑ ∑ i=1,2 K ∈Rh j=1

Ki ∇ (pi −

phi )

· nΠi φ

Aj+1 CAj

i

−

∑ ∑∫ i=1,2 σ ∈Ei

σ

[Πi (pi − phi )]{Ki ∇φiI }ds + α

I i ds

−

∑ ∑∫ i=1,2 σ ∈Ei

∑∑ i=1,2 σ ∈Ei

σ

[Πi φiI ]{Ki ∇ (pi − phi )}ds

[Πi (pi − phi )]σ [Πi φiI ]σ

) αf ( ξ (pi − phi ) + (1 − ξ )(pi+1 − phi+1 ) − (pf − phf ) Πi φiI ds + 2 ξ − 1 i=1,2( γ ) ∫ n ∑ ∂ (pf − phf ) ⏐⏐ ) αf ( + kf ,y d (p1 − ph1 ) + (p2 − ph2 ) − 2(pf − phf ) Πf φfI ds = 0. ⏐ Dℓ Πf φfI |σℓ | − yℓ ∂y γ 2ξ − 1 ∑∫

ℓ=1

(4.37)

3442

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

By Lemma 4.1 and the fact [Πi φiI ]σ = 0, ∀σ ∈ Ei , we obtain

∑

phi )

(Ki ∇h (pi −

, ∇h φ

I i )Ωi

+

∑ ∑ ∫ i=1,2 K ∈Rh

i=1,2

( ∂K

Ki ∇ (pi − phi ) · n(Πi φiI − φiI )ds

)

i

∑ ∑∫ ∑ ∑ ( ( ) ) [Πi (pi − phi )]{Ki ∇φiI }ds + ∇ · Ki ∇ (pi − phi ) , φiI − Πi φiI K − σ

i=1,2 σ ∈Ei

i=1,2 K ∈Rh i

) αf ( + ξ (pi − phi ) + (1 − ξ )(pi+1 − phi+1 ) − (pf − phf ) Πi φiI ds 2 ξ − 1 i=1,2( γ ) ∫ n ∑ ∂ (pf − phf ) ⏐⏐ ) αf ( (p1 − ph1 ) + (p2 − ph2 ) − 2(pf − phf ) Πf φfI ds = 0. + kf ,y d ⏐ Dℓ Πf φfI |σℓ | − yℓ ∂y γ 2ξ − 1 ∑∫

(4.38)

ℓ=1

Subtracting (4.38) from (4.36), it yields

∑ i=1,2

∥pi − phi ∥2L2 (Ω ) + ∥pf − phf ∥2L2 (γ ) i

∑( ∑ ∑ ( ( ) ) ) = ∇h (pi − phi ), Ki ∇h (φi − φiI ) Ω − ∇ · Ki ∇ (pi − phi ) , φiI − Πi φiI K i

i=1,2

⎛ +⎝

σ

i=1,2 σ ∈Ei

−

i=1,2 K ∈Rh i

∑ ∑∫

∑ ∑ ∫ (∫ kf ,y d

I i

]{Ki ∇φ }ds −

i=1,2 σ ∈Ei

(

∂K

i=1,2 K ∈Rh i

[Πi (pi −

phi )

⎞

∑ ∑∫ σ

{Ki ∇φi }[pi −

phi

]ds⎠

Ki ∇ (pi − phi ) · n(Πi φiI − φiI )ds

)

∂ (pf − phf ) ∂φf

n ∑

(4.39)

∂ (pf − phf )

(

)

) ⏐ ⏐ I ⏐ Dℓ Πf φf |σℓ |

ds − kf ,y d yℓ ∂y ∂y ∂y ℓ=1 ( ) αf + ξ (pi − phi ) + (1 − ξ )(pi+1 − phi+1 ) − (pf − phf ) (φi − Πi φiI )ds 2ξ − 1 i∫ =1,2 γ ) αf ( − (p1 − ph1 ) + (p2 − ph2 ) − 2(pf − phf ) (φf − Πf φfI )ds. 2 ξ − 1 γ

+

γ

∑∫

We rewrite the fourth item in the above equation as follows,

−

∑ ∑ ∫ ∂K

i=1,2 K ∈Rh

= −

i ∑ ∑ ∫

Ki ∇ (pi − phi ) · n(Πi φiI − φiI )ds

(

)

(4.40)

i

∑∫ γ

i=1,2

)

∂ K \γ

i=1,2 K ∈Rh

+

Ki ∇ (pi − phi ) · n(Πi φiI − φiI )ds

(

) αf ( ξ (pi − phi ) + (1 − ξ )(pi+1 − phi+1 ) − (pf − phf ) (Πi φiI − φiI )ds. 2ξ − 1

For the sixth item in (4.39), we have

∫ γ

kf ,y d

∫ =

kf ,y d

∂ (pf − phf ) ∂φf ∂y ∂ (pf −

∂y phf )

ds −

∂ (φf − φ

(

n ∫ ∑

ℓ=1 I f)

kf ,y d

σℓ

∫ ds +

k f ,y d

∂y ∂) y γ n h ∑ ∂ (pf − pf ) ⏐⏐ − kf ,y d ⏐ Dℓ Πf φfI |σℓ |, yℓ ∂y γ

∂ (pf − phf )

)

∂y ∂ (pf − ∂y

⏐ ⏐ ⏐ Dℓ Πf φfI ds yℓ

phf )

∂φfI ∂y

ds

(4.41)

(

ℓ=1

and

∫ ∫ I ∂ pf ∂φf αf I kf ,y d ds = gf φf ds + (p1 + p2 − 2pf )φfI ds, ∂ y ∂ y 2 ξ −1 γ γ γ ) ∫ ∫ n ( ∑ ∂ pf ⏐⏐ αf − kf ,y d (p1 + p2 − 2pf )Πf φfI ds. ⏐ Dℓ Πf φfI |σℓ | = − gf Πf φfI ds − yℓ ∂y γ γ 2ξ − 1 ∫

ℓ=1

(4.42)

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3443

Combining (4.41) with (4.42), it follows that

∫ kf ,y d

γ

∫ kf ,y d

= γ

−

∂ (pf − phf ) ∂φf ∂y ∂y

σℓ

k f ,y d

( kf ,y d

∂ (pf − phf )

∂y ∂ phf ∂φfI ∂y ∂y

ds +

∫ ds + γ

n ∑

)

∂y

ℓ=1

∂ (pf − phf ) ∂ (φf − φfI )

n ∫ ∑

ℓ=1

∂y

ds −

n ∑

⏐ ⏐ ⏐ Dℓ Πf φfI |σℓ | yℓ

gf (φfI − Πf φfI )ds

( kf ,y d

∂ phf

ℓ=1

)

∂y

⏐ ⏐ ⏐ Dℓ Πf φfI |σℓ |

(4.43)

yℓ

) αf ( (p1 − ph1 ) + (p2 − ph2 ) − 2(pf − phf ) (φfI − Πf φfI )ds 2ξ − 1

∫ + γ

∫ + γ

) αf ( h p1 + ph2 − 2phf (φfI − Πf φfI )ds. 2ξ − 1

Taking (4.40) and (4.43) back into (4.39), the L2 norm of p − ph reads as follows.

∑

∥pi − phi ∥2L2 (Ω ) + ∥pf − phf ∥2L2 (γ ) i

i=1,2

∫ ∑( ∂ (pf − phf ) ∂ (φf − φfI ) ) = ∇h (pi − phi ), Ki ∇h (φi − φiI ) Ω + kf ,y d ds i ∂y ∂y γ i=1,2 ∑ ∑ ( ( ) ) − ∇ · Ki ∇ (pi − phi ) , φiI − Πi φiI K (

)

i=1,2 K ∈Rh i

⎛

∑ ∑∫

+⎝

σ

i=1,2 σ ∈Ei

−

[Πi (pi − phi )]{Ki ∇φiI }ds −

i=1,2 σ ∈Ei

∑ ∑ ∫ i=1,2 K ∈Rh

∑ ∑∫

{Ki ∇φi }[pi − phi ]ds⎠

Ki ∇ (pi − phi ) · n(Πi φiI − φiI )ds

(

)

∂ phf ∂φfI

(

∂ K \γ

σ

⎞

(4.44)

i

−

n ∫ ∑

σℓ

ℓ=1

( +

k f ,y d

αf

∑∫ i=1,2

∂y ∂y

γ

2ξ − 1

(

ds +

n ∑

kf ,y d

ℓ=1

∂ phf ∂y

)

⏐ ⏐ ⏐ Dℓ Πf φfI |σℓ | yℓ

) ξ (pi − phi ) + (1 − ξ )(pi+1 − phi+1 ) − (pf − phf ) (φi − φiI )ds

) ) αf ( h h h I − (p1 − p1 ) + (p2 − p2 ) − 2(pf − pf ) (φf − φf )ds γ 2ξ − 1 (∫ ) ∫ ) αf ( h + gf (φfI − Πf φfI )ds + p1 + ph2 − 2phf (φfI − Πf φfI )ds γ γ 2ξ − 1 := I1 + I2 + I3 + I4 + I5 + I6 + I7 . ∫

We will estimate I1 − I7 defined in the above equation. Using Cauchy–Schwarz inequality, it is easy to see that

|I1 | ≤

∑ i=1,2

≤ Ch2

|pi − phi |1,h,Ωi |φi − φiI |1,h,Ωi + |pf − phf |H 1 (γ ) |φf − φfI |H 1 (γ )

∑ i=1,2

|pi |H 2 (Ωi ) |φi |H 2 (Ωi ) + Ch2 |pf |H 2 (γ ) |φf |H 2 (γ )

( ≤ Ch2

)( ∑

i=1,2

|pi |H 2 (Ωi ) + |pf |H 2 (γ )

(4.45)

) ∑ i=1,2

|φi |H 2 (Ωi ) + |φf |H 2 (γ ) .

3444

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

According to conclusions in Theorem 2.2 in [25], we have

( |I2 | ≤ Ch

)( ∑

2

∥pi ∥H 2 (Ωi ) +

i=1,2

|I3 | ≤ Ch

|I4 | ≤ Ch2

∥gi ∥H 1 (Ωi )

(i=1,2 ∑

) ∑

∥pi ∥H 2 (Ωi ) ∥pi ∥H 2 (Ωi )

) (i=1,2 ∑

i=1,2

∥φi ∥H 2 (Ωi ) ,

i=1,2

)( ∑

) ∑

i=1,2

( 2

∑

∥φi ∥H 2 (Ωi ) ,

(4.46)

) ∥φi ∥H 2 (Ωi ) .

i=1,2

We estimate I5 , since the approximation phf , φfI ∈ Wfh are all piecewise linear functions in γ , and using the second-order Taylor expression of kf ,y , we have

) ( n n ∫ ⏐ ⏐ ∑ ∑ ∂ phf ⏐⏐ ∂ phf ∂φfI ⏐ ⏐ ds + |I5 | = ⏐ − kf ,y d kf ,y d ⏐ Dℓ Πf φfI |σℓ |⏐ y ∂ y ∂ y ∂ y ℓ ℓ=1 ℓ=1 σℓ n ∫ h I ⏐∑ ⏐ ∂ p ∂φ ⏐ ⏐ f f (kf ,y − kf ,y |yℓ )d ≤ ⏐ ds⏐ ∂y ∂y σ ℓ ℓ=1 ( ) ( ) n ∫ n ⏐ ⏐∑ ∑ ∂ phf ⏐⏐ ∂φfI ∂ phf ⏐⏐ ⏐ ⏐ k f ,y d +⏐ ds − kf ,y d ⏐ ⏐ Dℓ Πf φfI |σℓ |⏐ y y ∂ y ∂ y ∂ y ℓ ℓ σℓ ℓ=1

(4.47)

ℓ=1

≤ Ch2 ∥pf ∥H 2 (γ ) ∥φf ∥H 2 (γ ) , where the last inequality is established due to n ∫ n ⏐∫ ⏐∑ ⏐⏐ ∂ ph ⏐ ⏐⏐ ∂φ I ⏐ ⏐ ∑ ∂ phf ∂φfI ⏐⏐ ⏐ ⏐⏐ f ⏐ ⏐⏐ f ⏐ ⏐ ⏐ (kf ,y − kf ,y |yℓ )d ds⏐ ≤ C ⏐ ⏐ (kf ,y − kf ,y |yℓ )ds⏐⏐ ⏐ ⏐⏐ ⏐ ⏐ ∂ y ∂ y ∂ y yℓ ∂ y yℓ σℓ ℓ=1 σℓ ℓ=1 n ⏐ h ⏐ ⏐⏐ I ⏐ ⏐ ∑ ⏐ ∂ pf ⏐ ⏐⏐ ∂φf ⏐ ⏐ ≤ Ch2 ∥kf ,y ∥W 2,∞ (γ ) ⏐ ⏐ ⏐⏐ ⏐ ⏐|σℓ | ≤ Ch2 ∥pf ∥H 2 (γ ) ∥φf ∥H 2 (γ ) , ∂ y yℓ ∂ y yℓ

ℓ=1

and

∂φfI ∂y

=

φfI (yℓ+1/2 ) − φfI (yℓ−1/2 ) |σℓ |

and Dℓ Πf φfI =

φfI (yℓ+1/2 ) − φfI (yℓ−1/2 ) |σℓ |

,

in σℓ .

For I6 , by Theorem 4.1 and (4.23), we obtain

⏐ ⏐∑∫ ) ( )) αf ( ( ⏐ |I6 | = ⏐ ξ (pi − phi ) − (pf − phf ) + (1 − ξ ) (pi+1 − phi+1 ) − (pf − phf ) ⏐ 2ξ − 1 i=1,2 γ ⏐ ( ) ⏐⏐ I I (φi − φi ) − (φf − φf ) ds⏐ ⏐ ( )1/2 ( )1/2 ∑ ∑ ≤ C ∥(pi − phi ) − (pf − phf )∥2L2 (γ ) ∥(φi − φiI ) − (φf − φfI )∥2L2 (γ ) i=1,2

(4.48)

i=1,2

≤ C |||p(− ph ||||||φ − φ I ||| )( ) ∑ ∑ 2 ≤ Ch ∥pi ∥H 2 (Ωi ) + ∥pf ∥H 2 (γ ) ∥φi ∥H 2 (Ωi ) + ∥φf ∥H 2 (γ ) . i=1,2

i=1,2

For the last term I7 , by identity (6) in Lemma 3.1, we get n ∫ ⏐ ⏐∑ ⏐∫ ⏐ ⏐ ⏐ ⏐ ⏐ (gf − gf (yℓ ))(φfI − Πf φfI )ds⏐ ⏐ gf (φfI − Πf φfI )ds⏐ = ⏐

γ

≤

n ∑

ℓ=1

ℓ=1

σℓ

∥gf − gf (yℓ )∥L2 (σℓ ) ∥φfI − Πf φfI ∥L2 (σℓ ) ≤ Ch2 ∥gf ∥H 1 (γ ) ∥φf ∥H 2 (γ ) .

(4.49)

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3445

Similarly, for the second term in I7 , it is obtained that

⏐ ) αf ( h ⏐ p1 + ph2 − 2phf (φfI − Πf φfI )ds⏐ γ 2ξ − 1 n ∫ ⏐ ⏐ 1 ⏐∑∑ ⏐ αf (phi − phf )(φfI − Πf φfI )ds⏐ ⏐ 2ξ − 1 i=1,2 ℓ=1 σℓ n ∫ ⏐ ⏐ (( h ) ( )) 1 ⏐∑∑ ⏐ αf pi − (αf pi )|yℓ − αf phf − (αf pf )|yℓ (φfI − Πf φfI )ds⏐ ⏐ 2ξ − 1 σ ℓ i=1,2 ℓ=1 n ∫ ⏐ ) ⏐⏐ )) ( I ) ( (( h 1 ⏐∑∑ φf − Πf φfI ds⏐ αf pi − αf pi − αf phf − αf pf ⏐ 2ξ − 1 i=1,2 ℓ=1 σℓ n ∫ ⏐ ⏐ ( ) 1 ⏐∑∑ ⏐ + (αf pi − (αf pi )|yℓ ) − (αf pf − (αf pf )|yℓ ) (φfI − Πf φfI )ds⏐ ⏐ 2ξ − 1 σ ℓ i=1,2 ℓ=1 ( ) ∑ h I I C |||p − p |||∥φf − Πf φf ∥L2 (γ ) + Ch∥αf ∥W 1,∞ (γ ) |pi |H 1 (γ ) + |pf |H 1 (γ ) ∥φfI − Πf φfI ∥L2 (γ ) i = 1 , 2 ( ) ∑ ∑ 2 Ch |pi |H 2 (Ωi ) + ∥pf ∥H 2 (γ ) + |pi |H 1 (γ ) ∥φf ∥H 2 (γ ) . ⏐∫ ⏐ ⏐

= = =

≤ ≤

i=1,2

(4.50)

i=1,2

Taking (4.49) and (4.50) into I7 , it follows that

( |I7 | ≤ Ch

2

) ∑

|pi |H 2 (Ωi ) + ∥pf ∥H 2 (γ ) +

i=1,2

∑

|pi |H 1 (γ ) + ∥gf ∥H 1 (γ ) ∥φf ∥H 2 (γ ) .

(4.51)

i=1,2

Taking the aforementioned seven terms back into (4.44), and using the regularity property (4.31), the error estimate for the pressure p on the L2 norm is obtained. □ Finally, we prove the convergence of the velocity u in (L2 )2 norm and the discrete H(div ) norm. Theorem 4.3. Let u = (u1 , u2 , uf ) be the exact solution for the fracture model (2.1) and uh = (uh1 , uh2 , uhf ) be the discontinuous finite volume approximation. There exists a positive constant C independent of mesh size h such that

( h

∥u − u ∥(L2 (Ω ))2 ≤ Ch ( h

|u − u |H(div),h ≤ Ch

) ∑

∥gi ∥L2 (Ωi ) + ∥gf ∥L2 (γ ) +

i=1,2

∑ i=1,2

∑

∥pi ∥H 2 (Ωi ) + ∥pf ∥H 2 (γ ) +

i=1,2

|gi |H 1 (Ωi ) + |gf |H 1 (γ ) +

∑

∥pi ∥H 2 (Ωi ) + ∥pf ∥H 2 (γ ) +

i=1,2

∑ i=1,2

∑

∥pi ∥L2 (γ ) , )

(4.52)

∥pi ∥H 1 (γ ) .

i=1,2

Proof. We first consider the (L2 (Ωi ))2 norm of ui − uhi . In matrix Ωi , from the first equation of (2.1) and (3.12), we have ui |K −uhi |K = −Ki ∇ pi |K +Ai,K ∇ phi |K −gi,K PK (x),

∀K ∈ Rhi .

(4.53)

By using Theorem 4.4 in [29], we obtain

( ) ∥ui − uhi ∥(L2 (Ωi ))2 ≤ C h∥gi ∥L2 (Ωi ) + |pi − phi |1,h,Ωi + h∥pi ∥H 1 (Ωi ) .

(4.54)

In the fracture γ , from the third equation of (2.1) and (3.13), we get uf |σℓ −uhf |σℓ = −kf ,y d

) ( ∂ phf ⏐⏐ ∂ pf ⏐⏐ αf ,ℓ (ph1,ℓ + ph2,ℓ − 2phf,ℓ ) (y − yℓ ). ⏐ + Af ,ℓ ⏐ − gf ,ℓ − ∂ y σℓ ∂ y σℓ 2ξ − 1

(4.55)

Hence,

 ⏐ ⏐ ∂ phf  ∂ pf αf ,ℓ  ⏐ ⏐  ∥uf − uhf ∥L2 (σℓ ) ≤ kf ,y d − Af ,ℓ (ph1,ℓ + ph2,ℓ − 2phf,ℓ )⏐∥y − yℓ ∥L2 (σℓ ) .  2 + ⏐gf ,ℓ − ∂y ∂ y L (σℓ ) 2ξ − 1

(4.56)

By the triangle inequality and assumption of kf ,y d, it yields that

 ∂ phf  ∂ pf   − Af ,ℓ kf ,y d  ∂y ∂ y L2 (σℓ )

≤ ≤

  ( ∂p ∂ phf ) ∂ pf ∂ pf      f − Af ,ℓ − kf ,y d  2 + Af ,ℓ  ∂y ∂ y L (σℓ ) ∂y ∂ y L2 (σℓ ) ( ) C |σℓ |∥pf ∥H 1 (σℓ ) + ∥pf − phf ∥H 1 (σℓ ) .

(4.57)

3446

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

We then consider the last term in (4.56), it is easy to check that

|gf ,ℓ | =

)1/2 (∫ )1/2 (∫ ⏐∫ ⏐ 1 1 ⏐ 2 1ds = g ds ∥gf ∥L2 (σℓ ) . g ds ≤ ⏐ ⏐ f f 1 |σℓ | σℓ |σℓ | σℓ σℓ |σℓ | 2 1 ⏐

(4.58)

Similarly, we have

|phi,ℓ − phf,ℓ | ≤

1

|σℓ |

3

1 2

∥phi − phf ∥L2 (σℓ ) , i = 1, 2. and ∥y − yℓ ∥L2 (σℓ ) ≤ |σℓ | 2 .

(4.59)

Using the triangle inequality, we have the following estimate for the last term in (4.56),

≤ ≤

⏐ ⏐ αf ,ℓ ⏐ ⏐ (ph1,ℓ + ph2,ℓ − 2phf,ℓ )⏐∥y − yℓ ∥L2 (σℓ ) ⏐gf ,ℓ − 2 ξ − 1 ( ) ∑ h h ∥pi − pf ∥L2 (σℓ ) C |σℓ | ∥gf ∥L2 (σℓ ) + i=1,2 ( ) ∑ ∑ ∑ h h C |σℓ | ∥gf ∥L2 (σℓ ) + ∥(pi − pi ) − (pf − pf )∥L2 (σℓ ) + ∥pi ∥L2 (σℓ ) + ∥pf ∥L2 (σℓ ) . i=1,2

i=1,2

(4.60)

i=1,2

Taking (4.57) and (4.60) back into (4.56), we obtain

( ∥ uf −

uhf L2 (σℓ )

∥

C

≤

|σℓ |∥pf ∥H 1 (σℓ ) + |σℓ |∥gf ∥L2 (σℓ ) + |σℓ |

∑ i=1,2

∥pi ∥L2 (σℓ ) + ∥pf − phf ∥H 1 (σℓ ) (4.61)

) ∑

+

∥(pi −

phi )

− (pf −

i=1,2

.

phf ) L2 (σℓ )

∥

Summing over ℓ = 1, . . . , n, we have

( ∥uf −

uhf L2 (γ )

∥

≤ C

h∥pf ∥H 1 (γ ) + h∥gf ∥L2 (γ ) + h

∑

∥pi ∥L2 (γ ) + ∥pf − phf ∥H 1 (γ )

i=1,2

(4.62)

) +

∑

∥(pi −

phi )

− (pf −

phf ) L2 (γ )

∥

.

i=1,2

Combining (4.62) with (4.54) and using Theorem 4.1, it is obtained that

( h

∥u − u ∥(L2 (Ω ))2 ≤ Ch

) ∑

∑

∑

∥gi ∥L2 (Ωi ) + ∥gf ∥L2 (γ ) + ∥pi ∥H 1 (Ωi ) + ∥pf ∥H 1 (γ ) + ∥pi ∥L2 (γ ) i=1,2 i=1,2 i=1,2 ( ) ∑ ∑ h h h h +C |pi − pi |1,h,Ωi + ∥pf − pf ∥H 1 (γ ) + ∥(pi − pi ) − (pf − pf )∥L2 (γ ) i = 1 , 2 i = 1 , 2 ( ) ∑ ∑ ∑ ≤ Ch ∥gi ∥L2 (Ωi ) + ∥gf ∥L2 (γ ) + ∥pi ∥H 2 (Ωi ) + ∥pf ∥H 2 (γ ) + ∥pi ∥L2 (γ ) . i=1,2

i=1,2

(4.63)

i=1,2

Thus, the desired result of the (L2 )2 norm holds. We then consider the discrete H(div ) norm for the velocity. In matrix Ωi , from the second equation of (2.1) and (3.12), it is obtained that

∥div ui − div uhi ∥L2 (K ) = ∥gi − gi,T ∥L2 (K ) ≤ Ch|gi |H 1 (K ) ,

K ∈ Rhi .

(4.64)

Summing over K ∈ Rhi , we have

|ui − uhi |H(div),Ωi ,h ≤ Ch|gi |H 1 (Ωi ) .

(4.65)

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3447

In the fracture γ , from the fourth equation of (2.1) and (3.13), we obtain

 ∂u ∂ uhf   f  −   ∂y ∂ y L2 (σℓ )

  α αf ,ℓ   f (p1 + p2 − 2pf ) − (ph1,ℓ + ph2,ℓ − 2phf,ℓ ) ≤ ∥gf − gf ,ℓ ∥L2 (σℓ ) +  L2 (σℓ ) 2ξ − 1 2ξ − 1 α − α   f  f ,ℓ ≤ ∥gf − gf ,ℓ ∥L2 (σℓ ) +  (p1 + p2 − 2pf ) L 2 (σ ℓ ) 2ξ − 1  ⏐ ⏐ α  ⏐ f ,ℓ ⏐ +⏐ ⏐(p1 − ph1,ℓ ) + (p2 − ph2,ℓ ) − 2(pf − phf,ℓ ) 2 2ξ ( −1 ) (L (σℓ ) ) ∑ ∑ ≤ C |σℓ | |gf |H 1 (σℓ ) + ∥pi ∥L2 (σℓ ) + ∥pf ∥L2 (σℓ ) + C ∥pi − pi,ℓ ∥L2 (σℓ ) + ∥pf − pf ,ℓ ∥L2 (σℓ ) i=1,2

+C

∑

∥(pi −

phi )

(4.66)

i=1,2

− (pf −

i=1,2

phf ) L2 (σℓ )

∥

.

Summing over ℓ = 1, . . . , n, we get

( |uf −

uhf 1,h,γ

|

)

≤ Ch |gf |H 1 (γ ) +

∑

∥pi ∥H 1 (γ ) + ∥pf ∥H 1 (γ ) +

i=1,2

∑

∥pi ∥H 2 (Ωi ) .

(4.67)

i=1,2

Combining (4.65) with (4.67), it yields that

(

) ∑

h

|u − u |H(div),h ≤ Ch

i=1,2

|gi |H 1 (Ωi ) + |gf |H 1 (γ ) +

∑

∥pi ∥H 2 (Ωi ) +

i=1,2

∑

∥pi ∥H 1 (γ ) + ∥pf ∥H 1 (γ ) . □

(4.68)

i=1,2

5. Numerical experiment In this section, numerical experiments are carried out to confirm the accuracy of the proposed discontinuous finite volume scheme. For simplicity, the domain for all examples is a rectangle, Ω = [−1, 1] × [0, 1], and the fracture γ = {x = 0}×[0, 1]. We form the triangulation of Ω by bisecting rectangles of the same shapes, and the finer triangulation can be obtained by connecting all midpoints of edges in the coarser grid. The integral of the right side functions gi over triangle K is computed by the midpoint rule using the three edges, and the integral of gf over edge σℓ is also computed by midpoint rule. We take the penalty parameter α = 10 in all of the experiments. In addition to testing the convergence rates for the L2 norm and |||·||| norm of the pressure p and the (L2 )2 norm of the velocity u, we also define the same discrete L2 norm of the flux as that in [29] as follows.

∥(u − uh ) · n∥20,h =

∑ ∑ ∑ (∫ i=1,2 K ∈Rh σ ∈EK

σ

(ui − uhi ) · nds

)2

.

i

Example 1. In the first example, the permeability in the fracture is much smaller than that in the surrounding domain so that the fracture can be seen as geologic barrier and the fluid barely flows through the fracture. The tensors of permeability in the surrounding domains and the fracture are all isotropic and constant, K1 = K2 = 10I, Kf = 10−3 I, where I is the two dimensional identity matrix. The width of the fracture d = 10−3 and ξ = 4/5. The right-hand side functions and Dirichlet boundary conditions are determined by the following exact solutions. We list the prior error estimates and convergence rates for p and u in Table 1. p1 = (x + 1) sin(2π y) p2 = −(4x + 1/2) sin(2π y) pf = sin(2π y).

{

Example 2. In the second example, a model with anisotropic fracture permeability is tested to verify the accuracy of the discontinuous finite volume method. The width of the fracture d = 10−3 , and ξ = 2/3. Dirichlet conditions are imposed on boundaries of fractures and surrounding domains. −3

K1 = K2 = 10

I,

( and Kf =

10−2 0

0 10

)

.

In this case, the tangential permeability in the fracture kf ,y = 10 and kf ,y d = 10−2 . The fracture is more permeable than surrounding matrix so fluid tends to fluid rapidly along the fracture. The exact solutions are shown as follows, and we

3448

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449 Table 1 Errors and convergence rates of DFVM for Example 1. h 1/8 1/16 1/32 1/64 1/128 1/256 1/512

∥p − ph ∥L2

Rate

||p − ph ||

Rate

∥u − uh ∥(L2 )2

Rate

∥(u − uh ) · n∥0,h

Rate

3.735e−02 9.704e−03 2.451e−03 6.143e−04 1.536e−04 3.841e−05 9.602e−06

– 1.94 1.98 1.99 2.00 2.00 2.00

2.209e−01 1.091e−01 5.431e−02 2.709e−02 1.353e−02 6.764e−03 3.381e-o3

– 1.02 1.01 1.00 1.00 1.00 1.00

1.798e−01 9.185e−02 4.619e−02 2.312e−02 1.156e−02 5.784e−03 2.892e−03

– 0.97 0.99 0.99 1.00 1.00 1.00

4.189e−02 1.105e−02 2.806e−03 7.045e−04 1.763e−04 4.410e−05 1.102e−05

– 1.92 1.97 1.99 2.00 2.00 2.00

Table 2 Errors and convergence rates of DFVM for Example 2. h 1/8 1/16 1/32 1/64 1/128 1/256 1/512

∥p − ph ∥L2

Rate

||p − ph ||

Rate

∥u − uh ∥(L2 )2

Rate

∥(u − uh ) · n∥0,h

Rate

1.294e−01 2.429e−02 4.264e−03 8.303e−04 1.860e−04 4.433e−05 1.084e−05

– 2.41 2.51 2.36 2.15 2.06 2.03

3.548e−01 1.653e−01 8.088e−02 4.010e−02 1.996e−02 9.960e−03 4.974e-o3

– 1.10 1.03 1.01 1.01 1.00 1.00

4.493e−01 2.085e−02 1.009e−02 4.971e−02 2.466e−02 1.228e−02 6.129e−03

– 1.10 1.04 1.02 1.01 1.00 1.00

1.735e−01 3.700e−02 6.157e−03 1.188e−03 2.619e−04 6.170e−05 1.499e−05

– 2.22 2.58 2.37 2.18 2.08 2.04

Table 3 Errors and convergence rates of DFVM for Example 3. h 1/8 1/16 1/32 1/64 1/128 1/256 1/512

∥p − ph ∥L2

Rate

||p − ph ||

Rate

∥u − uh ∥(L2 )2

Rate

∥(u − uh ) · n∥0,h

Rate

4.082e−02 1.038e−02 2.609e−03 6.531e−04 1.633e−04 4.083e−05 1.020e−05

– 1.97 1.99 1.99 2.00 2.00 2.00

1.961e−01 9.757e−02 4.870e−02 2.433e−02 1.216e−02 6.080e−03 3.040e-o3

– 1.08 1.01 1.00 1.00 1.00 1.00

1.822e−01 9.216e−02 4.623e−02 2.313e−02 1.156e−02 5.784e−03 2.892e−03

– 0.98 0.99 0.99 1.00 1.00 1.00

3.744e−02 9.728e−03 2.459e−03 6.166e−04 1.542e−04 3.857e−05 9.645e−06

– 1.94 1.98 1.99 1.99 2.00 2.00

show the numerical results in Table 2. p1 = (x + 1) sin(2π y) p2 = − (2x − (1 − 1/20000)) sin(2π y) pf = sin(2π y).

{

Example 3. In the last example, we are also concerned with a model with anisotropic fracture permeability, and there exists great difference between the tangential and normal permeability. d = 10−3 , ξ = 2/3, and K1 = K2 = I, kf ,y = 103 , kf ,y d = 1.

( Kf =

10−4 0

0 103

)

.

The Dirichlet boundary conditions in the fracture and surrounding domains are determined by the following exact solutions, and we have shown numerical results in Table 3. p1 = (x + 1) sin(2π y) p2 = −(2x + 4) sin(2π y) pf = sin(2π y).

{

Remark 5.1. According to the above three numerical examples, one can easily see that the proposed DFVM performs well for coupled fracture models in which fracture has a higher or lower permeability than the surrounding domains, or even has an anisotropic permeability with a big difference. Numerical results in Tables 1–3 show that errors for pressure p on the L2 norm and |||·||| norm hold O(h2 ) and O(h) accuracy, and the first-order convergence rate for velocity u on the (L2 )2 norm is obtained, which confirm our theoretical analysis in the last section. In addition, experiments also suggest a superconvergence property for the discrete L2 norm of the flux. Hence, we prove the accuracy and efficiency of the presented DFVM from two aspects of theory and practice.

S. Chen / Computers and Mathematics with Applications 78 (2019) 3429–3449

3449

6. Conclusion In this article, a DFVM for the hybrid-dimensional coupled fracture model with homogeneous Dirichlet boundary condition is proposed and analyzed. The pressure in matrix and fracture is approximated respectively by discontinuous and continuous piecewise linear functions. The DFVM eliminates the continuity of the numerical solution across the interelement boundary, hence it can be easily applied to the complex domain with fractures. The velocity is calculated in the discontinuous lowest-order Raviart–Thomas element space and piecewise linear element space in matrix and fracture. Hence, the presented method retains the local conservation property at the element level, similar to the mixed finite element method, but we avoid solving the indefinite saddle-point problem. Optimal order error estimates for the pressure p on the L2 norm and |||·||| norm, and for the velocity u on the (L2 )2 norm and the discrete H(div ) norm are proved. Numerical experiments with higher, lower and anisotropic fracture permeability are carried out and results show the same order accuracy for p and u as theoretical analysis, illustrating the accuracy and efficiency of the DFVM for the fracture model. Acknowledgment The author gratefully acknowledges financial support from the China Postdoctoral Science Foundation funded project (2019M650400). References [1] C. Alboin, J. Jaffre, J.E. Roberts, C. Serres, Domain decomposition for flow in fractured porous media, Domain Decompos. Methods Sci. Eng. (1999) 365–373. [2] V. Martin, J. Jaffre, J.E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, SIAM J. Sci. Comput. 26 (5) (2005) 1667–1691. [3] T.T.P. Hoang, C. Japhet, M. Kern, J.E. Roberts, Space-time domain decomposition for reduced fracture models in mixed formulation, SIAM J. Numer. Anal. 54 (1) (2016) 288–316. [4] M. Lesinigo, C. DAngelo, A. Quarteroni, A multiscale Darcy-Brinkman model for fluid flow in fractured porous media, Numer. Math. 117 (4) (2011) 717–752. [5] N. Frih, V. Martin, J.E. Roberts, A. Saâda, Modeling fractures as interfaces with nonmatching grids, Comput. Geosci. 16 (4) (2012) 1043–1060. [6] A. Fumagalli, A. Scotti, Numerical modelling of multiphase subsurface flow in the presence of fractures, Commun. Appl. Ind. Math. 3 (1) (2011). [7] W. Liu, Z. Sun, A block-centered finite difference method for reduced fracture model in karst aquifer system, Comput. Math. Appl. 74 (6) (2017) 1455–1470. [8] S. Chen, H. Rui, A node-centered finite volume method for a fracture model on triangulations, Appl. Math. Comput. 327 (2018) 55–69. [9] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, Comput. Math. Appl. 7 (6) (2000) 713–1018. [10] H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows (Ph.D. thesis), University of London, 1996. [11] S.H. Chou, D.Y. Kwak, K. Kim, Mixed finite volume methods on nonstaggered quadrilateral grids for elliptic problems, Math. Comput. 72 (242) (2003) 525–539. [12] H. Rui, T. Lu, An expanded mixed covolume method for elliptic problems, Numer. Methods Partial Differential Equations 21 (1) (2005) 8–23. [13] H. Rui, Symmetric modified finite volume element methods for self-adjoint elliptic and parabolic problems, J. Comput. Appl. Math. 146 (2) (2002) 373–386. [14] H. Rui, Symmetric mixed covolume methods for parabolic problems, Numer. Methods Partial Differential Equations 18 (5) (2002) 561–583. [15] H. Rui, A conservative characteristic finite volume element method for solution of the advection-diffusion equation, Comput. Methods Appl. Mech. Engrg. 197 (45–48) (2008) 3862–3869. [16] H. Rui, Analysis on a finite volume element method for Stokes problems, Acta Math. Appl. Sin. 21 (3) (2005) 359–372. [17] Z. Cai, On the finite volume element method, Numer. Math. 58 (1) (1991) 713–735. [18] Z. Cai, J. Mandel, S. McCormick, The finite volume element method for diffusion equations on general triangulations, SIAM J. Numer. Anal. 28 (2) (1991) 392–402. [19] P. Chatzipantelidis, A finite volume method based on the Crouzeix-Raviart element for elliptic pde’s in two dimensions, Numer. Math. 82 (3) (1999) 409–432. [20] C. Bi, L. Li, The mortar finite volume method with the crouzeix-raviart element for elliptic problems, Comput. Methods Appl. Mech. Engrg. 192 (1) (2003) 15–31. [21] H. Rui, C. Bi, Convergence analysis of an upwind finite volume element method with Crouzeix-Raviart element for non-selfadjoint and indefinite problems, Front. Math. China 3 (4) (2008) 563–579. [22] M. Crouzeix, P.A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations, Rev. Francaise Dautomatique Inform. Rech. Oper. 7 (3) (1973) 33–75. [23] X. Ye, A new discontinuous finite volume method for elliptic problems, SIAM J. Numer. Anal. 42 (3) (2005) 1062–1072. [24] X. Ye, A discontinuous finite volume method for the Stokes problems, SIAM J. Numer. Anal. 44 (1) (2006) 183–198. [25] S.H. Chou, X. Ye, Unified analysis of finite volume methods for second order elliptic problems, SIAM J. Numer. Anal. 45 (4) (2007) 1639–1653. [26] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (4) (1982) 742–760. [27] D.N. Arnold, F. Brezzi, B. Cockburn, D.L. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2002) 1749–1779. [28] C. Bi, M. Liu, A discontinuous finite volume element method for second-order elliptic problems, Numer. Methods Partial Differential Equations 28 (2) (2012) 425–440. [29] S.H. Chou, S. Tang, Conservative P1 conforming and nonconforming Galerkin FEMs: effective flux evaluation via a nonmixed method approach, SIAM J. Numer. Anal. 38 (2) (2000) 660–680. [30] P.A. Raviart, J.M. Thomas, A finite element method for 2nd order elliptic problem, Math. Asepects Finite Elem. Methods 606 (1977) 292–315.