- Email: [email protected]
A finite volume scheme for an elliptic equation with heterogeneous coefficients. Application to a homogenization problem
Applied Numerical Mathematics 38 (2001) 427–444 www.elsevier.com/locate/apnum
A finite volume scheme for an elliptic equation with heterogeneous coefficients. Application to a homogenization problem R. Luce, S. Perez Laboratoire de mathématiques appliquées, Université de Pau et des Pays de l’Adour, 64000 Pau, France
Abstract The aim of this paper is to develop a finite volume scheme to numerically solve an elliptic equation with heterogeneous coefficients, by using a regular mesh, independent of coefficients discontinuities. The construction of this scheme stems from a primal–dual mixed variational formulation with the solution decomposed into a regular part whose approximation belongs to a standard finite element space and an irregular part that is approximated by bubble functions. The irregular part is eliminated in terms of the regular part by imposing the continuity of the numerical flux across the lines of discontinuities of the coefficient. Several examples are presented and numerical experiments prove the efficiency of the developed method. 2001 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction In many applications, the data describing the properties of the medium (permeability, porosity) involve high discontinuities or numerous heterogeneities. Therefore, elliptic equations with discontinuous coefficients arise naturally in many mathematical modeling processes. In order to solve numerically such equations, we commonly use finite element or finite volume methods. Then, if we work with meshes adapted to the spatial distribution of the coefficients, we obtain a system of large size, that is prohitively expensive to solve. In order to avoid this, we would want to use any reasonable structured mesh. Thereby, some elements of the mesh will be crossed by the coefficient’s discontinuity surfaces and this implies to replace on these elements the true coefficient by coefficients that generate nearly equivalent description of the physical properties. This problem is often called interface problem. Upscaling or homogenization problems raise similar difficulties. Indeed, geophysical methods are continually being developed and improved. Therefore, they generate data describing on a fine scale the spatial distribution and the variability of the medium properties. Unfortunately, the level of these details are often out of scope for today’s computers. Thus, the data have to be upscaled or coarsened before the model could be used in a simulator. The upscaling is nontrivial because heterogeneities at all scales can have a significant effect. 0168-9274/01/$ – see front matter 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 0 4 0 - X
428
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
Such problems have been of much interest for many years. The first idea is to replace the true coefficient by a common average as harmonic, arithmetic or geometric mean. Warren and Price [11] demonstrated numerically that the geometric mean provides a good average for three dimensional systems comprising random permeabilities with no spatial structure. Bourgeat [2] applied homogenization theory to determine permeability but his results are only available for periodic systems. Another method is based on solving the flow equation on each coarse grid block (cf. [3]). A no flow boundary condition is applied to all faces of the block except for two opposite faces where constant Dirichlet boundary conditions are applied. The coarse grid block permeability is then determined such that the total flow across the block is preserved. This approach is more accurate but it is computationally more expensive. Nielsen and Tveito proposed an upscaling method based on a Weighted Output Least Squares approach. However, they obtain equivalent coefficients which do not characterize the medium, since they depend upon the boundary conditions and the source term. Ewing et al. have recently developed a modified finite volume approximation in [4]. The essence of the modification falls in the simultaneous discretization of any two normal components of the flux at the opposite faces of the finite volume. The continuous normal component of the flux through an interface is then approximated by finite difference. The purpose of this paper is to derive a discrete formulation which will allows us to obtain good approximations of the physical variables from non adapted meshes. The sought coefficients are effective in the way that they represent a property of the medium and do not vary with the flow conditions to which the medium is subjected. We describe the method in the context of an elliptic equation with Dirichlet boundary condition:
−div(a gradu) = f u=0
in Ω, on ∂Ω.
(1)
For the sake of simplicity, we take here Ω a rectangle of R2 but the results still hold on any bounded polygonal convex domain of R2 . The function f (x) is the given source term. The diffusion coefficient a is a piecewise constant function and hence, it has jump discontinuities across a given line Σ that is often called interface. In this paper, we will be only interested in two kinds of interfaces: • straight lines parallel to the grid lines, • transversal straight lines parallel among themselves. But, many other configurations can be considered without difficulty. The analogous 1D problem is wellknown by engineers. By assuming the continuity of the flux across the interfaces, we can easily prove that the harmonic mean is the sought coefficient. In order to derive a discrete formulation for the 2D problem, we will introduce a similar hypothesis and, so, we will need a formulation that impose regularity on the u variable and on the flux. The paper is organized as follows. In Section 2, we propose a primal–dual mixed formulation, which imposes some regularity on both the flux and the pressure. In order to take the coefficient’s discontinuities into account, we decompose, in Section 3, the variable u into a regular part and an irregular one, and we introduce a constraint in the discrete formulation. A finite volume scheme for which we state an existence and uniqueness result, is then built in Section 4. Finally, in Section 5, we give some numerical results and we conclude on the efficiency of our method. The notation that we use is standard. Unless otherwise specified, all functions spaces are considered over the domain Ω. We use · 0,Ω and | · |1,Ω to denote respectively the L2 -norm and the H 1 -seminorm. For other norms, we use an index, as, for instance, · H (div,Ω) denotes the usual norm in H (div, Ω).
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
429
2. Primal–dual mixed formulation We present below a primal–dual mixed formulation of (1) which is the following one: (p, u) ∈ H (div, Ω)× H01 (Ω), 2 2 pq dx − a graduq dx = 0, ∀q ∈ L (Ω) , ∀v ∈ L2 (Ω),
Ω
Ω
vdivp dx = −
Ω
(2)
f v dx. Ω
The diffusion coefficients belongs to the space of admissible parameters
Aad = a ∈ L∞ (Ω): 0 < α a β < ∞ , where α and β are positive constants. One can remark that this is not a classical mixed form; see, for instance, [8]. Besides, this formulation brings in the primal variable u and the dual variablep and we cannot straight eliminate the unknown p in order to obtain a u problem. However, using a numerical integration scheme on the two terms of the first equation in (2) allows to localize the problem and to obtain a standard five point finite volume scheme where only the variable u appears. Moreover, we choose this formulation instead of a more standard primal or dual mixed formulation because it allows us to impose the regularity on both the primal variable and the dual variable. The following result holds, cf. Thomas and Trujillo [9]. Theorem 1. The problem (2) has a unique solution (p, u) ∈ H (div, Ω) × H01 (Ω). When the coefficient a is smooth, the solution u of (1) belongs to the space H 2 (Ω). Here, since we only have a ∈ Aad , we globally lose this regularity but we locally preserve it on the parts where a is smooth. So, in order to take this remark into account, we decompose the primal variable u into a regular , i.e., part u and an irregular part u u = u + u .
3. The discrete formulation While the coefficient a has jump across the interface Σ , the average flux a∂u/∂η remains continuous across Σ . On the other hand, if the solution u is approximated in a standard piecewise polynomial finite space, the discrete flux a∂uh /∂η is discontinuous across the lines of discontinuities of the parameter. That’s why we impose on the discrete variables a constraint of the type
a Σ
+ + ∂uh
∂η
−a
− − ∂uh
∂η
ξh dΣ = 0 for every ξh ∈ Yh ,
(3)
where Σ represents the lines of discontinuities of a, Yh is a subspace of finite dimension of H 1/2(Σ) − − and a + ∂u+ h /∂η and a ∂uh /∂η represent the discrete flux on both sides of Σ . However, this relation is global; in order to localize it, we will introduce later some numerical integration operators.
430
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
3.1. Meshes and discretization spaces We define below three meshes on Ω, which will allow us to introduce the finite volume method for (2). Let (Uh )h be a regular family of triangulations formed by rectangles whose edges are parallel to the coordinates axis. Such a triangulation will be associated to the approximation of the primal variable u. The number of internal vertices of Uh is noted nv. By connecting the centers of gravity of the elements of Uh by straight line segments, we build a triangulation Ph of Ω, which is called the dual one of Uh . It will be associated to the dual variable p (see Fig. 1). The family (Ph )h is also a regular family of triangulations of Ω. The number of internal edges of Ph is denoted by ne. By connecting next the midpoints of the horizontal edges of the elements of Uh by straight line segments, we obtain the mesh Q1h . Doing likewise with the midpoints of the vertical edges gives us the mesh Q2h . We next denote by Qh the mesh whose elements belong either to Q1h or to Q2h . It will be associated to the test function q. Moreover, we introduce (see also Fig. 1), for every K of Qh , the piecewise constant vectorial function eK defined on K, equal to 01 if K ∈ Q1h and equal to 10 if K ∈ Q2h . Let the points xi , i = 1, . . . , n, be the intersections of Σ, the lines of discontinuity of a, with the edges of the elements of Uh . Next, for every such point xj , j = 1, . . . , nd, which is internal to Ω, we define the intervals Σj in the following way: Σj contains xj and its extremities are the midpoints of [xk , xj ] and [xj , xl ] with xk , xl (1 k, l n) the neighboring points situated on both sides of xj (see Fig. 2).
Fig. 1.
Fig. 2.
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
431
We next denote Σh = {Σ1 , Σ2 , . . . , Σnd }. Finally, we define the bubble functions wj in the following way: • the wj support is constituted by the two elements of Uh whose closure contains the point xj , • the function wj is a piecewise linear polynomial function, • the function wj is equal to 1 at the point xj and vanishes at any xk , with k = j , as well as at any vertex of Uh . We can now introduce the discretization spaces. First, we choose the space associated with the vectorial unknown ph
W1h = ph ∈ H (div, Ω): ∀K ∈ Ph , p h|K ∈ P1,0 (K) × P0,1 (K) , where P1,0 is defined as the space of polynomial functions linear in the first variable, constant in the second variable. P0,1 is defined in a similar way. Then, we take for the vectorial test functions q h the following space:
2
W2h = q h ∈ L2 (Ω) : ∀K ∈ Qh , q h|K ∈ E(K) , where E(K) is the space of vectorial functions proportional to eK . As regards the scalar unknown uh , we firstly introduce two spaces associated with the regular and the irregular components:
1h = u h ∈ H01 (Ω): ∀K ∈ Uh , u h |K ∈ Q1 (K) Z
and
1h = u h ∈ H01 (Ω): u h ∈ W(Ω) , Z
where Q1 is the space of functions which are linear in each variable and W(Ω) is spanned by the bubble functions wj . 1h and Z 1h : We denote Z1h the direct sum of Z 1h ⊕ Z 1h . Z1h = Z
Then, we choose for the test functions ξh the space
Yh = ξh ∈ H 1/2 (Σ): ∀Σi ∈ Σh , ξh |Σi ∈ P0 (Σi ) . P0 (Σi ) is the space of constant functions on Σi . Finally, we associate to the scalar unknown uh the space
M1h = uh ∈ Z1h :
R a (uh )ηξh dΣ = 0 for every ξh ∈ Yh ,
Σ
where R a is an operator which allows to integrate numerically the constraint (3). It will be define in the following subsection. We still have to choose a discretization space for the scalar test functions vh , which we take as
M2h = vh ∈ L2 (Ω): ∀K ∈ Ph , vh |K ∈ P0 (K) , where P0 (K) is the space of constant functions on K.
432
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
With all these notations, the discrete problem is now written as below: (ph , uh ) ∈ W1h× M1h , ∀q h ∈ W2h ,
p h q h dx −
Ω
∀vh ∈ M2h ,
a graduh q h dx = 0,
Ω
vh divph dx = −
Ω
(4)
f vh dx. Ω
3.2. Characterization of M1h The constraint (3) included in the space M1h is global. In order to localize it, we introduce an integration operator Ra . More precisely, we define an operator R a for the regular part and an operator R a for the irregular part. In order to simplify the definitions, we suppose that every element K of Qh is crossed by at the most one straight line of discontinuity of a. We denote by K + and K − the parts of K placed on both sides of this straight line and by a + and a − the respective values of a. Then, we define the two integration operators in the following way. 1h , we denote R a (u h ) the element of W2h defined by Definition 2. For every u h ∈ Z
∀K ∈ Qh ,
a u h |K = R
a + gradu h eK |γK eK + − a gradu h eK |γK eK −
on K + , on K − ,
where eK + and eK − are the restrictions of eK to K + and K − , respectively. γK is the edge of an element of Uh , contained in K. 1h , we denote R a (u h ) the element of W2h defined by: ∀K ∈ Qh , Definition 3. For every u h ∈ Z
R a u h |K =
a + gradu h |K + eK + |γK eK + − a gradu h |K − eK − |γK eK −
on K + , on K − ,
where eK + and eK − are the restrictions of eK to K + and K − , respectively. γK is the edge of an element of Uh , contained in K. We do not present a general characterization of M1h but we are only interested here in two examples. We firstly study a homogenization problem and secondly we consider the case of a transversal discontinuity which crosses the elements of the triangulation (interface problem). 3.2.1. Homogenization problem We suppose that every element K of Uh is formed by four identical subdomains Ki , on which the coefficient a takes different values (cf. Fig. 3). We denote by ui , i = 1, . . . , 6, the values of the regular h at the vertices of two elements K1 and K2 of Uh , and u 1 the value of the irregular part at the part u intersection point of ∂K1 ∩ ∂K2 with the line of discontinuity Σ of a. The real numbers h1K and h2K represent the dimensions of the rectangle K. Then 1 − h2K u3 − u4 u − (a3 + a5 ) R a u h η + R a u h η w1 dΣ = + h2K (a3 + a5 ) 2 h1K h1K Σ
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
433
Fig. 3.
and
+
+
R a u h η + R a u h η w1 dΣ =
Σ
1 h2K u3 − u4 u (a4 + a6 ) − h2K (a4 + a6 ) . 2 h1K h1K
So, by imposing the constraint (3), we obtain u 1 =
1 a3 + a5 − a4 − a6 (u3 − u4 ). 2 a3 + a4 + a5 + a6
We prove in a similar way that any values u i can be expressed as a linear combination, of ui . Moreover, if we evaluate the flux through Σ1 , we obtain
∂uh (a3 + a5 )(a4 + a6 ) u3 − u4 a = . ∂η Σ1 a3 + a4 + a5 + a6 h1K
This expression leads us to introduce the following equivalent coefficient a∗ =
(a3 + a5 )(a4 + a6 ) a3 + a4 + a5 + a6
(5)
which verifies all the constraints imposed. One can notice that by taking a3 = a5 and a4 = a6 , the expression (5) represents exactly the harmonic mean; this particular case actually privileges one direction and we thus obtain a result already known in 1D. Moreover, when taking a3 = a6 and a4 = a5 , one find the arithmetic mean for a ∗ . i , By doing likewise on every edge of any element K of Uh , there appear four coefficients aK ∗ i = 1, . . . , 4, which have the same form as a . Moreover, for every uh ∈ M1h and every K ∈ Uh , we can write a graduh under the form 1 2 aK (u1 − u2 ) + aK (u3 − u4 ) 1 u2 − u1 a + y K h1K h1K h2K (6) (a graduh )|K = 3 4 aK (u1 − u4 ) + aK (u3 − u2 ) 3 u4 − u1 + x aK h2K h1K h2K where h1K and h2K are the dimensions of the rectangle K.
434
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
Fig. 4.
If we now suppose that every element of Uh is formed by sixteen identical subdomains Ki , as represented in Fig. 4, we have to express three values of the irregular part (u i , i = 1, . . . , 3) by using two values of the regular part. By imposing the constraint (3), we obtain the system 1 − β u 2 = 14 (α − β)(u3 − u4 ), (α + β)u −β u 1 + (β + δ)u 2 − γ u 3 = 14 (β − γ )(u3 − u4 ), 1
−γ u2 + (γ + δ)u3 = 4 (γ − δ)(u3 − u4 )
with
α = 14 (a9 + a13 + a17 + a21 ), β = 1 (a10 + a14 + a18 + a22 ), 4
γ = 14 (a11 + a15 + a19 + a23 ), 1
δ = 4 (a12 + a16 + a20 + a24 ).
This system has a unique solution 1 αβγ + αβδ + αγ δ − 3βγ δ (u3 − u4 ), 4 αβγ + αβδ + αγ δ + βγ δ 1 αβγ + αβδ − αγ δ − βγ δ 2 = (u3 − u4 ), u 2 αβγ + αβδ + αγ δ + βγ δ 1 αβδ + αγ δ + βγ δ − 3αβγ 3 = (u3 − u4 ). u 4 αβγ + αβδ + αγ δ + βγ δ 1 = u
If we evaluate the flux through Σ2 , we obtain
a
∂uh 4αβγ δ (u3 − u4 ). = ∂η Σ2 αβγ + αβδ + αγ δ + βγ δ
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
435
Fig. 5.
This expression leads us to introduce the following equivalent coefficient: a∗ =
4αβγ δ . αβγ + αβδ + αγ δ + βγ δ
We can remark that a ∗ is the harmonic mean of the arithmetic means α, β, γ and δ. Doing likewise on every edge of each element K ∈ Uh gives us four equivalent coefficients. Thus, we can write a graduh under the form (6). 3.2.2. Interface problem The expression of the irregular part depends on the position of Σ, and in particular, on its inclination relatively to the coordinates axis. We do not make a complete study, we are only interested in the case presented in Fig. 5. In order to simplify the description of the method, we suppose here that the line of discontinuity Σ crosses the edges of the elements of Uh at the third of their length; one has to use exactly the same technique when the position of Σ is different. In this case, we only use an integration operator for the irregular part and we evaluate exactly the regular component. We obtain
a−
Σ1
a Σ1
3 1 ∂ u − − h + R a u h η dΣ = a1 2(u8 − u6 ) + (u5 − u3 ) − a1 u 1 , ∂η 6 2
+ + ∂ uh
∂η
3 1 + + R a u h η dΣ = a2 2(u8 − u6 ) + (u5 − u3 ) + a2 u 1 . 6 4
Next, imposing the constraint (3) gives us u 1 =
2 a1 − a2 2(u8 − u6 ) + (u5 − u3 ) . 9 2a1 + a2
We can always express the u i as a linear combination of the uj . Moreover, by calculating the flux through the line of discontinuity Σ of a, we obtain on each element of Uh four coefficients which can be used to express a graduh under the form (6).
436
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
3.3. Concluding remark Let us end this section by drawing a conclusion on the general case. On the one hand, we remark that we can always express the irregular component of the solution from the regular one. So we necessarily have the compatibility of the discrete spaces dimensions: dim W1h = dim W2h = ne, dim M1h = dim M2h = nv. On the other hand, one can see that the elements uh ∈ M1h can be expressed in every element K of Uh , under the form 1 2 aK (u1 − u2 ) + aK (u3 − u4 ) 1 u2 − u1 a + y K h1K h1K h2K , (a graduh )|K = 3 4 a u − u (u − u ) + a (u − u ) 1 1 4 2 K 3 3 4 + K x aK h2K h1K h2K 1 2 3 , aK , aK where u1 , u2 , u3 , u4 are the values of the regular part u h at the four vertices of K and where aK 4 and aK are the different values of a on K obtained by the means of the constraint (3). The real numbers h1K and h2K represent the dimensions of the rectangle K. In order to numerically solve this discrete problem, we will build a finite volume scheme from these remarks.
4. Finite volume scheme 4.1. Numerical integrations In order to eliminate the variable p h in (4), we define below two integration operators for both the primal variable uh and the dual variable ph . For every element K of Ph , ηK denotes the unit exterior normal of K and (p h ηK )i , i = 1, . . . , 4 represent the flux across the four edges of K. We can remark that (p h eK )i = ±(p h ηK )i .
(7)
Definition 4. For every ph ∈ W1h , we denote R(p h ) the element of W2h defined by
R(p h ) =
(p h eK )|γK eK ,
K∈Qh
where γK is the edge of an element of Ph , enclosed in K. Lemma 5. There exist two constants c1 and c2 independent of h such that for every p h ∈ W1h , one has
c1 p h 0,Ω R(p h )0,Ω c2 p h 0,Ω . Proof. For every K ∈ Ph , we have R(p )2
h1K h2K (ph ηK )21 + (ph ηK )22 + (ph ηK )23 + (ph ηK )24 2 and also, using the Simpson formula gives h
0,K
=
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
p h 20,K =
437
2 h1K h2K (ph ηK )21 + (ph ηK )1 + (ph ηK )2 + (p h ηK )22 6 2 h1K h2K + (p h ηK )23 + (p h ηK )3 + (p h ηK )4 + (ph ηK )24 . 6
So, by summing on every element K of Ph , we obtain the result with c1 = 1 and c2 =
√
✷
3.
Definition 6. For every uh ∈ M1h , we denote R ∗a (uh ) the element of W2h defined by
R ∗a (uh ) =
(a graduh )|K eK |γK eK
K∈Qh
where γK is the edge of an element of Uh , enclosed in K and (a graduh )|K is determined from the constraint (3). Lemma 7. There exist two constants c3 and c4 independent of h such that for every uh ∈ M1h ,
c3 a graduh 0,Ω R ∗a (uh )0,Ω c4 a graduh 0,Ω . Proof. For every element K of Uh , we have a graduh 20,K =
2 1 2 2 2
h2K 1 2 aK (u2 − u1 ) + aK (u2 − u1 ) + aK (u3 − u4 ) + aK (u3 − u4 ) 6h1K 2 3 2 4 2
h1K 3 4 + aK (u4 − u1 ) + aK (u4 − u1 ) + aK (u3 − u2 ) + aK (u3 − u2 ) , 6h2K
and ∗ R (uh ) a
0,K
=
2 2 2
h2K 1 aK (u2 − u1 ) + aK (u3 − u4 ) 2h1K 2 4 2
h1K 3 + aK (u4 − u1 ) + aK (u3 − u2 ) . 2h2K
So, using the inequality (a + b)2 2(a 2 + b2 ) gives the result with c3 = 1 and c4 =
√
3.
✷
In order to establish inf–sup conditions for the discrete problem (4), we give the following discrete Green formula for the operators R and R ∗a . Lemma 8. For every uh ∈ M1h and p h ∈ W1h , we have
R ∗a (uh )R(ph ) dx = −
Ω
a divph Π0h (uh ) dx, Ω
is defined on each element K ∈ Ph by the value of uh at the center of K, weighted by the where equivalent coefficient a obtained by the means of (3). a Π0h
Proof. The definitions of the operators R and R ∗a provide the equality
Ω
R ∗a (uh )R(ph ) dx =
K∈Qh
mes(K) (a graduh )|K eK σK (ph eK )|γK ,
(8)
438
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
where σK is the edge of an element of Uh , enclosed in K, and γK is the edge of an element of Ph , enclosed in K. Moreover, an integration by parts on every element K of Ph leads us to −
a divp h Π0h (uh ) dx = −
uh (mK ) h2K a1K (ph ηK )1 + a2K (p h ηK )2
K∈Qh
Ω
+ h1K a3K (ph ηK )3 + a4K (p h ηK )4 ,
(9)
where uh (mK ) is the value taken by uh at the center of K ∈ Ph and ηK is the exterior normal to K. Then, thanks to the relations (7), we reorganize (9) to obtain (8). For more details, we refer to [7]. ✷ By introducing these two operators R and R ∗a in the discrete formulation (4), we obtain the following formulation: (ph , uh ) ∈ W1h× M1h , ∀q h ∈ W2h , ∀vh ∈ M2h ,
R ∗a (uh )q h dx = 0,
R(ph )q h dx −
Ω
Ω
vh divph dx = −
Ω
(10)
f vh dx. Ω
Concerning the numerical solving of this system, let us remark that one can calculate p h from the second equation. By replacing it into the first one, one thus obtain a 5 points scheme in the unknown uh . 4.2. Inf–sup conditions In order to establish the well-posedness of the problem (10), we define the spaces V1h = {ph ∈ W1h : divp h = 0},
V2h = q h ∈ W2h :
R ∗a (uh )q h dx = 0 ∀uh ∈ M1h ,
Ω
and we state some results. Lemma 9. There exists a constant β1 > 0 independent of h such that
vh divp h dx β1 > 0. vh ∈M2h p ∈W1h vh 0,Ω p h H (div,Ω) h inf
Ω
sup
Proof. Let us consider ϕv ∈ H01 (Ω) the unique solution of the auxiliary problem −2ϕv = vh with vh ∈ M2h . We then denote p h = −gradϕv ∈ H (div, Ω) and thus we obtain the equality
vh divp h dx = −
Ω
vh div gradϕv dx = vh 20,Ω .
Ω
Moreover, we have
p h H (div,Ω) = gradϕv 20,Ω + vh 20,Ω
1/2
= |ϕv |21,Ω + vh 20,Ω
1/2
Cvh 0,Ω .
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
439
Therefore, we obtain
vh divp h dx = vh 20,Ω
Ω
1 vh 0,Ω p h H (div,Ω) C
and the stated result follows with β1 = 1/C. One can find a similar result in [8].
✷
Lemma 10. There exists a constant β2 > 0 independent of h such that
inf
sup
uh ∈M1h q ∈W2h h
R ∗a (uh )q h dx β2 > 0. |uh |1,Ω q h 0,Ω Ω
Proof. Let us consider an element uh of M1h and let us denote q h = R ∗a (uh ) ∈ W2h . From Lemma 7 and the definition of the class Aad , we derive the following inequality:
R ∗a (uh )q h dx = R ∗a (uh )0,Ω q h 0,Ω C|uh |1,Ω q h 0,Ω
Ω
where the positive constant C is equal to αc3 .
✷
Lemma 11. There exists a constant β3 > 0 independent of h such that
R(ph )q h dx β3 > 0. ph ∈V1h q ∈V2h p h H (div,Ω) q h 0,Ω h inf
sup
Ω
Proof. Let us consider an element p h ∈ V1h . We denote q h = R(ph ) ∈ W2h . Using the discrete Green formula (Lemma 8) gives, for every uh ∈ M1h , the equality
R ∗a (uh )R(ph ) dx = −
Ω
a divp h Π0h (uh ) dx. Ω
But p h is an element of the space V1h . So divp h = 0 and the following result holds:
R ∗a (uh )R(ph ) dx = 0.
Ω
Therefore q h really is an element of V2h and Lemma 5 involves
R(ph )q h dx = R(ph )0,Ω q h 0,Ω c1 p h 0,Ω q h 0,Ω .
Ω
In this way, we prove the result with β3 = c1 .
✷
Now, it comes the main statement of the section. Theorem 12. The problem has a unique solution (p h , uh ) ∈ W1h × M1h . Proof. We have to verify three hypotheses: (1) The dimensions of the discrete spaces are compatible: dim W1h + dim M1h = dim W2h + dim M2h .
440
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
(2) The Babu˘ska–Brezzi conditions are verified (cf. Lemma 9–11). (3) The continuity of the bilinear forms is satisfied: √ R(ph )q h dx 3p h H (div,Ω) q h 0,Ω , Ω
vh divph dx p h H (div,Ω) vh 0,Ω ,
Ω
R ∗a (uh )q h dx
√
3uh H 1 (Ω) q h 0,Ω . 0
Ω
Then, thanks to [1] or [5], we obtain the result.
✷
5. Numerical experiments In the examples that follow, we shall present some computational results obtained by using our method. The experiments were performed on a UNIX-station in double precision arithmetic. The source term is chosen as
f (x, y) = 100 sin(2π x) sin(2πy) and the integral
K
for every (x, y) ∈ (0, 1) × (0, 1),
f dx is evaluated thanks to the midpoint formula.
Example 1. In this first example, the coefficient a is equal to 1 all over the domain Ω except on a barrier R described by (x, y) ∈ ]0.5375, 0.5875[×]0.125, 0.875[, where a takes the value 100. On the elements which contain the vertices of R, the true coefficient is replaced by the constant c = 2(1 + 100)/(3 + 100). On the elements crossed by the edges of R, the effective coefficient is equal to the harmonic mean. The solution obtained with an adapted mesh and the solution evaluated with a coarser non-adapted mesh are represented in Fig. 6.
Fig. 6.
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
441
Fig. 7 represents on the coarse mesh the absolute difference between the two approximated solutions described in Fig. 6. We can remark that the main error is made on a par with the vertices of the barrier. The relative L∞ -error is of the order of 10−2 .
Example 2. Here, the parameter a is laid out in checkerboard (1 or 100). Fig. 8 represents the solution obtained with an adapted mesh and the solution evaluated with a coarser mesh whose elements contain four elements of the refined mesh. All over the domain Ω, the true coefficient is replaced by the arithmetic mean which is equal to the constant c = 50.5.
Fig. 7.
Fig. 8.
442
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
Fig. 9.
The absolute difference between the two approximated solution is represented in Fig. 9. Here, the relative L∞ -error is of the order of 10−3 . Example 3. In this example, we are interested in the case of a transversal discontinuity described in Section 3.2.2. We consider the true coefficient a defined by
a(x, y) =
c1 c2
1 if y > x − 120 , 1 if y x − 120 ,
where c1 and c2 are strictly positive constants. Fig. 10 represents the solution obtained with an adapted mesh including 31000 elements and a structured non adapted mesh of size 80 × 80. The constants are chosen as c1 = 1000 and c2 = 1. According to the Table 1, we can remark that the L2 -error between the true solution and the approximated one is divided by four when we use an 80 × 80 mesh instead of a 40 × 40 mesh. So, we can expect a convergence of the order of h2 . For all the examples, the results are the same when we change the source term f and they are always satisfactory when the coefficients are modified. As a conclusion, we have constructed in this paper a five points scheme to numerically solve an elliptic equation with discontinuous coefficients. But we can also define a nine points scheme if we choose to evaluate exactly the uh term and to approximate the p h term. Moreover, we can make additional remarks. First, we were only interested in interfaces parallel to the grid lines and in a particular transversal interface on Cartesian meshes. However, many others configurations can be considered but they will require calculations whose complexity depends on the position of Σ . In the same way, using nonstructured meshes is thinkable but it makes difficult the elimination of the irregular part of the solution. On the other hand, the theoretical results (inf–sup conditions, existence and uniqueness theorem, etc.) remain unchanged. Next, only scalar coefficients have been considered. But the presented method can easily fit
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
443
Fig. 10. Table 1 L2 -error for example 3 c1 = 100 and c2 = 1
c1 = 1000 and c2 = 1
Mesh 40 × 40
3.77 · 10−2
4.26 · 10−2
Mesh 80 × 80
1.22 · 10−2
1.33 · 10−2
diagonal matrix coefficient and our intention is to study full matrix coefficient. Lastly, this type of finite volume method can be used to evaluate the solution also at some internal points of any element K of Uh ; this is done by studying a Dirichlet problem set on K. Therefore, our results provide a starting point for developing multigrid methods.
References [1] C. Bernardi, C. Canuto, Y. Maday, Generalized inf–sup conditions for Chebyshev spectral approximation of the Stokes problem, SIAM J. Numer. Anal. 25 (1988) 1237–1271. [2] A. Bourgeat, Homogeneized behavior of two-phase flows in naturally fractured reservoirs with uniform fractures distribution, Comput. Methods Appl. Mech. Eng. 47 (1984) 205–216. [3] L.J. Durlofsky, Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media, Water Resources Res. 27 (1991) 699–708. [4] R. Ewing, O. Iliev, R. Lazarov, A modified finite volume approximation of second-order elliptic equations with discontinuous coefficients, Technical Report, ISC-99-01-MATH, Texas A&M University, 1999. [5] R.A. Nicolaïdes, Existence, uniqueness and approximation for generalized saddle-point problems, SIAM J. Numer. Anal. 19 (1982) 349–357. [6] B.F. Nielsen, A. Tveito, An upscaling method for one-phase flow in heterogeneous reservoirs. A weighted output least squares (WOLS) approach, Comput. Geosci. 2 (1998) 93–123. [7] S. Perez, Identification and homogénéisation de paramètres dans des équations aux dérivées partielles Thèse, Université de Pau and des Pays de l’Adour, 1999.
444
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
[8] J.E. Roberts, J.M. Thomas, Mixed and Hybrid Methods, Handbook of Numerical Analysis, Vol. II, NorthHolland, Amsterdam. [9] J.M. Thomas, D. Trujillo, Mixed finite volume methods, Internat. J. Numer. Methods Eng. 46 (1999) 1351– 1366. [10] D. Trujillo, Couplage espace temps de schémas numériques en simulation de réservoir, Thèse, Université de Pau and des Pays de l’Adour, 1994. [11] J.E. Warren, H.S. Price, Flow in heterogeneous porous media, Soc. Petroleum Eng. J. (1961) 153–169.
A finite volume scheme for an elliptic equation with heterogeneous coefficients. Application to a homogenization problem R. Luce, S. Perez Laboratoire de mathématiques appliquées, Université de Pau et des Pays de l’Adour, 64000 Pau, France
Abstract The aim of this paper is to develop a finite volume scheme to numerically solve an elliptic equation with heterogeneous coefficients, by using a regular mesh, independent of coefficients discontinuities. The construction of this scheme stems from a primal–dual mixed variational formulation with the solution decomposed into a regular part whose approximation belongs to a standard finite element space and an irregular part that is approximated by bubble functions. The irregular part is eliminated in terms of the regular part by imposing the continuity of the numerical flux across the lines of discontinuities of the coefficient. Several examples are presented and numerical experiments prove the efficiency of the developed method. 2001 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction In many applications, the data describing the properties of the medium (permeability, porosity) involve high discontinuities or numerous heterogeneities. Therefore, elliptic equations with discontinuous coefficients arise naturally in many mathematical modeling processes. In order to solve numerically such equations, we commonly use finite element or finite volume methods. Then, if we work with meshes adapted to the spatial distribution of the coefficients, we obtain a system of large size, that is prohitively expensive to solve. In order to avoid this, we would want to use any reasonable structured mesh. Thereby, some elements of the mesh will be crossed by the coefficient’s discontinuity surfaces and this implies to replace on these elements the true coefficient by coefficients that generate nearly equivalent description of the physical properties. This problem is often called interface problem. Upscaling or homogenization problems raise similar difficulties. Indeed, geophysical methods are continually being developed and improved. Therefore, they generate data describing on a fine scale the spatial distribution and the variability of the medium properties. Unfortunately, the level of these details are often out of scope for today’s computers. Thus, the data have to be upscaled or coarsened before the model could be used in a simulator. The upscaling is nontrivial because heterogeneities at all scales can have a significant effect. 0168-9274/01/$ – see front matter 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 0 4 0 - X
428
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
Such problems have been of much interest for many years. The first idea is to replace the true coefficient by a common average as harmonic, arithmetic or geometric mean. Warren and Price [11] demonstrated numerically that the geometric mean provides a good average for three dimensional systems comprising random permeabilities with no spatial structure. Bourgeat [2] applied homogenization theory to determine permeability but his results are only available for periodic systems. Another method is based on solving the flow equation on each coarse grid block (cf. [3]). A no flow boundary condition is applied to all faces of the block except for two opposite faces where constant Dirichlet boundary conditions are applied. The coarse grid block permeability is then determined such that the total flow across the block is preserved. This approach is more accurate but it is computationally more expensive. Nielsen and Tveito proposed an upscaling method based on a Weighted Output Least Squares approach. However, they obtain equivalent coefficients which do not characterize the medium, since they depend upon the boundary conditions and the source term. Ewing et al. have recently developed a modified finite volume approximation in [4]. The essence of the modification falls in the simultaneous discretization of any two normal components of the flux at the opposite faces of the finite volume. The continuous normal component of the flux through an interface is then approximated by finite difference. The purpose of this paper is to derive a discrete formulation which will allows us to obtain good approximations of the physical variables from non adapted meshes. The sought coefficients are effective in the way that they represent a property of the medium and do not vary with the flow conditions to which the medium is subjected. We describe the method in the context of an elliptic equation with Dirichlet boundary condition:
−div(a gradu) = f u=0
in Ω, on ∂Ω.
(1)
For the sake of simplicity, we take here Ω a rectangle of R2 but the results still hold on any bounded polygonal convex domain of R2 . The function f (x) is the given source term. The diffusion coefficient a is a piecewise constant function and hence, it has jump discontinuities across a given line Σ that is often called interface. In this paper, we will be only interested in two kinds of interfaces: • straight lines parallel to the grid lines, • transversal straight lines parallel among themselves. But, many other configurations can be considered without difficulty. The analogous 1D problem is wellknown by engineers. By assuming the continuity of the flux across the interfaces, we can easily prove that the harmonic mean is the sought coefficient. In order to derive a discrete formulation for the 2D problem, we will introduce a similar hypothesis and, so, we will need a formulation that impose regularity on the u variable and on the flux. The paper is organized as follows. In Section 2, we propose a primal–dual mixed formulation, which imposes some regularity on both the flux and the pressure. In order to take the coefficient’s discontinuities into account, we decompose, in Section 3, the variable u into a regular part and an irregular one, and we introduce a constraint in the discrete formulation. A finite volume scheme for which we state an existence and uniqueness result, is then built in Section 4. Finally, in Section 5, we give some numerical results and we conclude on the efficiency of our method. The notation that we use is standard. Unless otherwise specified, all functions spaces are considered over the domain Ω. We use · 0,Ω and | · |1,Ω to denote respectively the L2 -norm and the H 1 -seminorm. For other norms, we use an index, as, for instance, · H (div,Ω) denotes the usual norm in H (div, Ω).
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
429
2. Primal–dual mixed formulation We present below a primal–dual mixed formulation of (1) which is the following one: (p, u) ∈ H (div, Ω)× H01 (Ω), 2 2 pq dx − a graduq dx = 0, ∀q ∈ L (Ω) , ∀v ∈ L2 (Ω),
Ω
Ω
vdivp dx = −
Ω
(2)
f v dx. Ω
The diffusion coefficients belongs to the space of admissible parameters
Aad = a ∈ L∞ (Ω): 0 < α a β < ∞ , where α and β are positive constants. One can remark that this is not a classical mixed form; see, for instance, [8]. Besides, this formulation brings in the primal variable u and the dual variablep and we cannot straight eliminate the unknown p in order to obtain a u problem. However, using a numerical integration scheme on the two terms of the first equation in (2) allows to localize the problem and to obtain a standard five point finite volume scheme where only the variable u appears. Moreover, we choose this formulation instead of a more standard primal or dual mixed formulation because it allows us to impose the regularity on both the primal variable and the dual variable. The following result holds, cf. Thomas and Trujillo [9]. Theorem 1. The problem (2) has a unique solution (p, u) ∈ H (div, Ω) × H01 (Ω). When the coefficient a is smooth, the solution u of (1) belongs to the space H 2 (Ω). Here, since we only have a ∈ Aad , we globally lose this regularity but we locally preserve it on the parts where a is smooth. So, in order to take this remark into account, we decompose the primal variable u into a regular , i.e., part u and an irregular part u u = u + u .
3. The discrete formulation While the coefficient a has jump across the interface Σ , the average flux a∂u/∂η remains continuous across Σ . On the other hand, if the solution u is approximated in a standard piecewise polynomial finite space, the discrete flux a∂uh /∂η is discontinuous across the lines of discontinuities of the parameter. That’s why we impose on the discrete variables a constraint of the type
a Σ
+ + ∂uh
∂η
−a
− − ∂uh
∂η
ξh dΣ = 0 for every ξh ∈ Yh ,
(3)
where Σ represents the lines of discontinuities of a, Yh is a subspace of finite dimension of H 1/2(Σ) − − and a + ∂u+ h /∂η and a ∂uh /∂η represent the discrete flux on both sides of Σ . However, this relation is global; in order to localize it, we will introduce later some numerical integration operators.
430
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
3.1. Meshes and discretization spaces We define below three meshes on Ω, which will allow us to introduce the finite volume method for (2). Let (Uh )h be a regular family of triangulations formed by rectangles whose edges are parallel to the coordinates axis. Such a triangulation will be associated to the approximation of the primal variable u. The number of internal vertices of Uh is noted nv. By connecting the centers of gravity of the elements of Uh by straight line segments, we build a triangulation Ph of Ω, which is called the dual one of Uh . It will be associated to the dual variable p (see Fig. 1). The family (Ph )h is also a regular family of triangulations of Ω. The number of internal edges of Ph is denoted by ne. By connecting next the midpoints of the horizontal edges of the elements of Uh by straight line segments, we obtain the mesh Q1h . Doing likewise with the midpoints of the vertical edges gives us the mesh Q2h . We next denote by Qh the mesh whose elements belong either to Q1h or to Q2h . It will be associated to the test function q. Moreover, we introduce (see also Fig. 1), for every K of Qh , the piecewise constant vectorial function eK defined on K, equal to 01 if K ∈ Q1h and equal to 10 if K ∈ Q2h . Let the points xi , i = 1, . . . , n, be the intersections of Σ, the lines of discontinuity of a, with the edges of the elements of Uh . Next, for every such point xj , j = 1, . . . , nd, which is internal to Ω, we define the intervals Σj in the following way: Σj contains xj and its extremities are the midpoints of [xk , xj ] and [xj , xl ] with xk , xl (1 k, l n) the neighboring points situated on both sides of xj (see Fig. 2).
Fig. 1.
Fig. 2.
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
431
We next denote Σh = {Σ1 , Σ2 , . . . , Σnd }. Finally, we define the bubble functions wj in the following way: • the wj support is constituted by the two elements of Uh whose closure contains the point xj , • the function wj is a piecewise linear polynomial function, • the function wj is equal to 1 at the point xj and vanishes at any xk , with k = j , as well as at any vertex of Uh . We can now introduce the discretization spaces. First, we choose the space associated with the vectorial unknown ph
W1h = ph ∈ H (div, Ω): ∀K ∈ Ph , p h|K ∈ P1,0 (K) × P0,1 (K) , where P1,0 is defined as the space of polynomial functions linear in the first variable, constant in the second variable. P0,1 is defined in a similar way. Then, we take for the vectorial test functions q h the following space:
2
W2h = q h ∈ L2 (Ω) : ∀K ∈ Qh , q h|K ∈ E(K) , where E(K) is the space of vectorial functions proportional to eK . As regards the scalar unknown uh , we firstly introduce two spaces associated with the regular and the irregular components:
1h = u h ∈ H01 (Ω): ∀K ∈ Uh , u h |K ∈ Q1 (K) Z
and
1h = u h ∈ H01 (Ω): u h ∈ W(Ω) , Z
where Q1 is the space of functions which are linear in each variable and W(Ω) is spanned by the bubble functions wj . 1h and Z 1h : We denote Z1h the direct sum of Z 1h ⊕ Z 1h . Z1h = Z
Then, we choose for the test functions ξh the space
Yh = ξh ∈ H 1/2 (Σ): ∀Σi ∈ Σh , ξh |Σi ∈ P0 (Σi ) . P0 (Σi ) is the space of constant functions on Σi . Finally, we associate to the scalar unknown uh the space
M1h = uh ∈ Z1h :
R a (uh )ηξh dΣ = 0 for every ξh ∈ Yh ,
Σ
where R a is an operator which allows to integrate numerically the constraint (3). It will be define in the following subsection. We still have to choose a discretization space for the scalar test functions vh , which we take as
M2h = vh ∈ L2 (Ω): ∀K ∈ Ph , vh |K ∈ P0 (K) , where P0 (K) is the space of constant functions on K.
432
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
With all these notations, the discrete problem is now written as below: (ph , uh ) ∈ W1h× M1h , ∀q h ∈ W2h ,
p h q h dx −
Ω
∀vh ∈ M2h ,
a graduh q h dx = 0,
Ω
vh divph dx = −
Ω
(4)
f vh dx. Ω
3.2. Characterization of M1h The constraint (3) included in the space M1h is global. In order to localize it, we introduce an integration operator Ra . More precisely, we define an operator R a for the regular part and an operator R a for the irregular part. In order to simplify the definitions, we suppose that every element K of Qh is crossed by at the most one straight line of discontinuity of a. We denote by K + and K − the parts of K placed on both sides of this straight line and by a + and a − the respective values of a. Then, we define the two integration operators in the following way. 1h , we denote R a (u h ) the element of W2h defined by Definition 2. For every u h ∈ Z
∀K ∈ Qh ,
a u h |K = R
a + gradu h eK |γK eK + − a gradu h eK |γK eK −
on K + , on K − ,
where eK + and eK − are the restrictions of eK to K + and K − , respectively. γK is the edge of an element of Uh , contained in K. 1h , we denote R a (u h ) the element of W2h defined by: ∀K ∈ Qh , Definition 3. For every u h ∈ Z
R a u h |K =
a + gradu h |K + eK + |γK eK + − a gradu h |K − eK − |γK eK −
on K + , on K − ,
where eK + and eK − are the restrictions of eK to K + and K − , respectively. γK is the edge of an element of Uh , contained in K. We do not present a general characterization of M1h but we are only interested here in two examples. We firstly study a homogenization problem and secondly we consider the case of a transversal discontinuity which crosses the elements of the triangulation (interface problem). 3.2.1. Homogenization problem We suppose that every element K of Uh is formed by four identical subdomains Ki , on which the coefficient a takes different values (cf. Fig. 3). We denote by ui , i = 1, . . . , 6, the values of the regular h at the vertices of two elements K1 and K2 of Uh , and u 1 the value of the irregular part at the part u intersection point of ∂K1 ∩ ∂K2 with the line of discontinuity Σ of a. The real numbers h1K and h2K represent the dimensions of the rectangle K. Then 1 − h2K u3 − u4 u − (a3 + a5 ) R a u h η + R a u h η w1 dΣ = + h2K (a3 + a5 ) 2 h1K h1K Σ
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
433
Fig. 3.
and
+
+
R a u h η + R a u h η w1 dΣ =
Σ
1 h2K u3 − u4 u (a4 + a6 ) − h2K (a4 + a6 ) . 2 h1K h1K
So, by imposing the constraint (3), we obtain u 1 =
1 a3 + a5 − a4 − a6 (u3 − u4 ). 2 a3 + a4 + a5 + a6
We prove in a similar way that any values u i can be expressed as a linear combination, of ui . Moreover, if we evaluate the flux through Σ1 , we obtain
∂uh (a3 + a5 )(a4 + a6 ) u3 − u4 a = . ∂η Σ1 a3 + a4 + a5 + a6 h1K
This expression leads us to introduce the following equivalent coefficient a∗ =
(a3 + a5 )(a4 + a6 ) a3 + a4 + a5 + a6
(5)
which verifies all the constraints imposed. One can notice that by taking a3 = a5 and a4 = a6 , the expression (5) represents exactly the harmonic mean; this particular case actually privileges one direction and we thus obtain a result already known in 1D. Moreover, when taking a3 = a6 and a4 = a5 , one find the arithmetic mean for a ∗ . i , By doing likewise on every edge of any element K of Uh , there appear four coefficients aK ∗ i = 1, . . . , 4, which have the same form as a . Moreover, for every uh ∈ M1h and every K ∈ Uh , we can write a graduh under the form 1 2 aK (u1 − u2 ) + aK (u3 − u4 ) 1 u2 − u1 a + y K h1K h1K h2K (6) (a graduh )|K = 3 4 aK (u1 − u4 ) + aK (u3 − u2 ) 3 u4 − u1 + x aK h2K h1K h2K where h1K and h2K are the dimensions of the rectangle K.
434
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
Fig. 4.
If we now suppose that every element of Uh is formed by sixteen identical subdomains Ki , as represented in Fig. 4, we have to express three values of the irregular part (u i , i = 1, . . . , 3) by using two values of the regular part. By imposing the constraint (3), we obtain the system 1 − β u 2 = 14 (α − β)(u3 − u4 ), (α + β)u −β u 1 + (β + δ)u 2 − γ u 3 = 14 (β − γ )(u3 − u4 ), 1
−γ u2 + (γ + δ)u3 = 4 (γ − δ)(u3 − u4 )
with
α = 14 (a9 + a13 + a17 + a21 ), β = 1 (a10 + a14 + a18 + a22 ), 4
γ = 14 (a11 + a15 + a19 + a23 ), 1
δ = 4 (a12 + a16 + a20 + a24 ).
This system has a unique solution 1 αβγ + αβδ + αγ δ − 3βγ δ (u3 − u4 ), 4 αβγ + αβδ + αγ δ + βγ δ 1 αβγ + αβδ − αγ δ − βγ δ 2 = (u3 − u4 ), u 2 αβγ + αβδ + αγ δ + βγ δ 1 αβδ + αγ δ + βγ δ − 3αβγ 3 = (u3 − u4 ). u 4 αβγ + αβδ + αγ δ + βγ δ 1 = u
If we evaluate the flux through Σ2 , we obtain
a
∂uh 4αβγ δ (u3 − u4 ). = ∂η Σ2 αβγ + αβδ + αγ δ + βγ δ
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
435
Fig. 5.
This expression leads us to introduce the following equivalent coefficient: a∗ =
4αβγ δ . αβγ + αβδ + αγ δ + βγ δ
We can remark that a ∗ is the harmonic mean of the arithmetic means α, β, γ and δ. Doing likewise on every edge of each element K ∈ Uh gives us four equivalent coefficients. Thus, we can write a graduh under the form (6). 3.2.2. Interface problem The expression of the irregular part depends on the position of Σ, and in particular, on its inclination relatively to the coordinates axis. We do not make a complete study, we are only interested in the case presented in Fig. 5. In order to simplify the description of the method, we suppose here that the line of discontinuity Σ crosses the edges of the elements of Uh at the third of their length; one has to use exactly the same technique when the position of Σ is different. In this case, we only use an integration operator for the irregular part and we evaluate exactly the regular component. We obtain
a−
Σ1
a Σ1
3 1 ∂ u − − h + R a u h η dΣ = a1 2(u8 − u6 ) + (u5 − u3 ) − a1 u 1 , ∂η 6 2
+ + ∂ uh
∂η
3 1 + + R a u h η dΣ = a2 2(u8 − u6 ) + (u5 − u3 ) + a2 u 1 . 6 4
Next, imposing the constraint (3) gives us u 1 =
2 a1 − a2 2(u8 − u6 ) + (u5 − u3 ) . 9 2a1 + a2
We can always express the u i as a linear combination of the uj . Moreover, by calculating the flux through the line of discontinuity Σ of a, we obtain on each element of Uh four coefficients which can be used to express a graduh under the form (6).
436
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
3.3. Concluding remark Let us end this section by drawing a conclusion on the general case. On the one hand, we remark that we can always express the irregular component of the solution from the regular one. So we necessarily have the compatibility of the discrete spaces dimensions: dim W1h = dim W2h = ne, dim M1h = dim M2h = nv. On the other hand, one can see that the elements uh ∈ M1h can be expressed in every element K of Uh , under the form 1 2 aK (u1 − u2 ) + aK (u3 − u4 ) 1 u2 − u1 a + y K h1K h1K h2K , (a graduh )|K = 3 4 a u − u (u − u ) + a (u − u ) 1 1 4 2 K 3 3 4 + K x aK h2K h1K h2K 1 2 3 , aK , aK where u1 , u2 , u3 , u4 are the values of the regular part u h at the four vertices of K and where aK 4 and aK are the different values of a on K obtained by the means of the constraint (3). The real numbers h1K and h2K represent the dimensions of the rectangle K. In order to numerically solve this discrete problem, we will build a finite volume scheme from these remarks.
4. Finite volume scheme 4.1. Numerical integrations In order to eliminate the variable p h in (4), we define below two integration operators for both the primal variable uh and the dual variable ph . For every element K of Ph , ηK denotes the unit exterior normal of K and (p h ηK )i , i = 1, . . . , 4 represent the flux across the four edges of K. We can remark that (p h eK )i = ±(p h ηK )i .
(7)
Definition 4. For every ph ∈ W1h , we denote R(p h ) the element of W2h defined by
R(p h ) =
(p h eK )|γK eK ,
K∈Qh
where γK is the edge of an element of Ph , enclosed in K. Lemma 5. There exist two constants c1 and c2 independent of h such that for every p h ∈ W1h , one has
c1 p h 0,Ω R(p h )0,Ω c2 p h 0,Ω . Proof. For every K ∈ Ph , we have R(p )2
h1K h2K (ph ηK )21 + (ph ηK )22 + (ph ηK )23 + (ph ηK )24 2 and also, using the Simpson formula gives h
0,K
=
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
p h 20,K =
437
2 h1K h2K (ph ηK )21 + (ph ηK )1 + (ph ηK )2 + (p h ηK )22 6 2 h1K h2K + (p h ηK )23 + (p h ηK )3 + (p h ηK )4 + (ph ηK )24 . 6
So, by summing on every element K of Ph , we obtain the result with c1 = 1 and c2 =
√
✷
3.
Definition 6. For every uh ∈ M1h , we denote R ∗a (uh ) the element of W2h defined by
R ∗a (uh ) =
(a graduh )|K eK |γK eK
K∈Qh
where γK is the edge of an element of Uh , enclosed in K and (a graduh )|K is determined from the constraint (3). Lemma 7. There exist two constants c3 and c4 independent of h such that for every uh ∈ M1h ,
c3 a graduh 0,Ω R ∗a (uh )0,Ω c4 a graduh 0,Ω . Proof. For every element K of Uh , we have a graduh 20,K =
2 1 2 2 2
h2K 1 2 aK (u2 − u1 ) + aK (u2 − u1 ) + aK (u3 − u4 ) + aK (u3 − u4 ) 6h1K 2 3 2 4 2
h1K 3 4 + aK (u4 − u1 ) + aK (u4 − u1 ) + aK (u3 − u2 ) + aK (u3 − u2 ) , 6h2K
and ∗ R (uh ) a
0,K
=
2 2 2
h2K 1 aK (u2 − u1 ) + aK (u3 − u4 ) 2h1K 2 4 2
h1K 3 + aK (u4 − u1 ) + aK (u3 − u2 ) . 2h2K
So, using the inequality (a + b)2 2(a 2 + b2 ) gives the result with c3 = 1 and c4 =
√
3.
✷
In order to establish inf–sup conditions for the discrete problem (4), we give the following discrete Green formula for the operators R and R ∗a . Lemma 8. For every uh ∈ M1h and p h ∈ W1h , we have
R ∗a (uh )R(ph ) dx = −
Ω
a divph Π0h (uh ) dx, Ω
is defined on each element K ∈ Ph by the value of uh at the center of K, weighted by the where equivalent coefficient a obtained by the means of (3). a Π0h
Proof. The definitions of the operators R and R ∗a provide the equality
Ω
R ∗a (uh )R(ph ) dx =
K∈Qh
mes(K) (a graduh )|K eK σK (ph eK )|γK ,
(8)
438
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
where σK is the edge of an element of Uh , enclosed in K, and γK is the edge of an element of Ph , enclosed in K. Moreover, an integration by parts on every element K of Ph leads us to −
a divp h Π0h (uh ) dx = −
uh (mK ) h2K a1K (ph ηK )1 + a2K (p h ηK )2
K∈Qh
Ω
+ h1K a3K (ph ηK )3 + a4K (p h ηK )4 ,
(9)
where uh (mK ) is the value taken by uh at the center of K ∈ Ph and ηK is the exterior normal to K. Then, thanks to the relations (7), we reorganize (9) to obtain (8). For more details, we refer to [7]. ✷ By introducing these two operators R and R ∗a in the discrete formulation (4), we obtain the following formulation: (ph , uh ) ∈ W1h× M1h , ∀q h ∈ W2h , ∀vh ∈ M2h ,
R ∗a (uh )q h dx = 0,
R(ph )q h dx −
Ω
Ω
vh divph dx = −
Ω
(10)
f vh dx. Ω
Concerning the numerical solving of this system, let us remark that one can calculate p h from the second equation. By replacing it into the first one, one thus obtain a 5 points scheme in the unknown uh . 4.2. Inf–sup conditions In order to establish the well-posedness of the problem (10), we define the spaces V1h = {ph ∈ W1h : divp h = 0},
V2h = q h ∈ W2h :
R ∗a (uh )q h dx = 0 ∀uh ∈ M1h ,
Ω
and we state some results. Lemma 9. There exists a constant β1 > 0 independent of h such that
vh divp h dx β1 > 0. vh ∈M2h p ∈W1h vh 0,Ω p h H (div,Ω) h inf
Ω
sup
Proof. Let us consider ϕv ∈ H01 (Ω) the unique solution of the auxiliary problem −2ϕv = vh with vh ∈ M2h . We then denote p h = −gradϕv ∈ H (div, Ω) and thus we obtain the equality
vh divp h dx = −
Ω
vh div gradϕv dx = vh 20,Ω .
Ω
Moreover, we have
p h H (div,Ω) = gradϕv 20,Ω + vh 20,Ω
1/2
= |ϕv |21,Ω + vh 20,Ω
1/2
Cvh 0,Ω .
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
439
Therefore, we obtain
vh divp h dx = vh 20,Ω
Ω
1 vh 0,Ω p h H (div,Ω) C
and the stated result follows with β1 = 1/C. One can find a similar result in [8].
✷
Lemma 10. There exists a constant β2 > 0 independent of h such that
inf
sup
uh ∈M1h q ∈W2h h
R ∗a (uh )q h dx β2 > 0. |uh |1,Ω q h 0,Ω Ω
Proof. Let us consider an element uh of M1h and let us denote q h = R ∗a (uh ) ∈ W2h . From Lemma 7 and the definition of the class Aad , we derive the following inequality:
R ∗a (uh )q h dx = R ∗a (uh )0,Ω q h 0,Ω C|uh |1,Ω q h 0,Ω
Ω
where the positive constant C is equal to αc3 .
✷
Lemma 11. There exists a constant β3 > 0 independent of h such that
R(ph )q h dx β3 > 0. ph ∈V1h q ∈V2h p h H (div,Ω) q h 0,Ω h inf
sup
Ω
Proof. Let us consider an element p h ∈ V1h . We denote q h = R(ph ) ∈ W2h . Using the discrete Green formula (Lemma 8) gives, for every uh ∈ M1h , the equality
R ∗a (uh )R(ph ) dx = −
Ω
a divp h Π0h (uh ) dx. Ω
But p h is an element of the space V1h . So divp h = 0 and the following result holds:
R ∗a (uh )R(ph ) dx = 0.
Ω
Therefore q h really is an element of V2h and Lemma 5 involves
R(ph )q h dx = R(ph )0,Ω q h 0,Ω c1 p h 0,Ω q h 0,Ω .
Ω
In this way, we prove the result with β3 = c1 .
✷
Now, it comes the main statement of the section. Theorem 12. The problem has a unique solution (p h , uh ) ∈ W1h × M1h . Proof. We have to verify three hypotheses: (1) The dimensions of the discrete spaces are compatible: dim W1h + dim M1h = dim W2h + dim M2h .
440
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
(2) The Babu˘ska–Brezzi conditions are verified (cf. Lemma 9–11). (3) The continuity of the bilinear forms is satisfied: √ R(ph )q h dx 3p h H (div,Ω) q h 0,Ω , Ω
vh divph dx p h H (div,Ω) vh 0,Ω ,
Ω
R ∗a (uh )q h dx
√
3uh H 1 (Ω) q h 0,Ω . 0
Ω
Then, thanks to [1] or [5], we obtain the result.
✷
5. Numerical experiments In the examples that follow, we shall present some computational results obtained by using our method. The experiments were performed on a UNIX-station in double precision arithmetic. The source term is chosen as
f (x, y) = 100 sin(2π x) sin(2πy) and the integral
K
for every (x, y) ∈ (0, 1) × (0, 1),
f dx is evaluated thanks to the midpoint formula.
Example 1. In this first example, the coefficient a is equal to 1 all over the domain Ω except on a barrier R described by (x, y) ∈ ]0.5375, 0.5875[×]0.125, 0.875[, where a takes the value 100. On the elements which contain the vertices of R, the true coefficient is replaced by the constant c = 2(1 + 100)/(3 + 100). On the elements crossed by the edges of R, the effective coefficient is equal to the harmonic mean. The solution obtained with an adapted mesh and the solution evaluated with a coarser non-adapted mesh are represented in Fig. 6.
Fig. 6.
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
441
Fig. 7 represents on the coarse mesh the absolute difference between the two approximated solutions described in Fig. 6. We can remark that the main error is made on a par with the vertices of the barrier. The relative L∞ -error is of the order of 10−2 .
Example 2. Here, the parameter a is laid out in checkerboard (1 or 100). Fig. 8 represents the solution obtained with an adapted mesh and the solution evaluated with a coarser mesh whose elements contain four elements of the refined mesh. All over the domain Ω, the true coefficient is replaced by the arithmetic mean which is equal to the constant c = 50.5.
Fig. 7.
Fig. 8.
442
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
Fig. 9.
The absolute difference between the two approximated solution is represented in Fig. 9. Here, the relative L∞ -error is of the order of 10−3 . Example 3. In this example, we are interested in the case of a transversal discontinuity described in Section 3.2.2. We consider the true coefficient a defined by
a(x, y) =
c1 c2
1 if y > x − 120 , 1 if y x − 120 ,
where c1 and c2 are strictly positive constants. Fig. 10 represents the solution obtained with an adapted mesh including 31000 elements and a structured non adapted mesh of size 80 × 80. The constants are chosen as c1 = 1000 and c2 = 1. According to the Table 1, we can remark that the L2 -error between the true solution and the approximated one is divided by four when we use an 80 × 80 mesh instead of a 40 × 40 mesh. So, we can expect a convergence of the order of h2 . For all the examples, the results are the same when we change the source term f and they are always satisfactory when the coefficients are modified. As a conclusion, we have constructed in this paper a five points scheme to numerically solve an elliptic equation with discontinuous coefficients. But we can also define a nine points scheme if we choose to evaluate exactly the uh term and to approximate the p h term. Moreover, we can make additional remarks. First, we were only interested in interfaces parallel to the grid lines and in a particular transversal interface on Cartesian meshes. However, many others configurations can be considered but they will require calculations whose complexity depends on the position of Σ . In the same way, using nonstructured meshes is thinkable but it makes difficult the elimination of the irregular part of the solution. On the other hand, the theoretical results (inf–sup conditions, existence and uniqueness theorem, etc.) remain unchanged. Next, only scalar coefficients have been considered. But the presented method can easily fit
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
443
Fig. 10. Table 1 L2 -error for example 3 c1 = 100 and c2 = 1
c1 = 1000 and c2 = 1
Mesh 40 × 40
3.77 · 10−2
4.26 · 10−2
Mesh 80 × 80
1.22 · 10−2
1.33 · 10−2
diagonal matrix coefficient and our intention is to study full matrix coefficient. Lastly, this type of finite volume method can be used to evaluate the solution also at some internal points of any element K of Uh ; this is done by studying a Dirichlet problem set on K. Therefore, our results provide a starting point for developing multigrid methods.
References [1] C. Bernardi, C. Canuto, Y. Maday, Generalized inf–sup conditions for Chebyshev spectral approximation of the Stokes problem, SIAM J. Numer. Anal. 25 (1988) 1237–1271. [2] A. Bourgeat, Homogeneized behavior of two-phase flows in naturally fractured reservoirs with uniform fractures distribution, Comput. Methods Appl. Mech. Eng. 47 (1984) 205–216. [3] L.J. Durlofsky, Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media, Water Resources Res. 27 (1991) 699–708. [4] R. Ewing, O. Iliev, R. Lazarov, A modified finite volume approximation of second-order elliptic equations with discontinuous coefficients, Technical Report, ISC-99-01-MATH, Texas A&M University, 1999. [5] R.A. Nicolaïdes, Existence, uniqueness and approximation for generalized saddle-point problems, SIAM J. Numer. Anal. 19 (1982) 349–357. [6] B.F. Nielsen, A. Tveito, An upscaling method for one-phase flow in heterogeneous reservoirs. A weighted output least squares (WOLS) approach, Comput. Geosci. 2 (1998) 93–123. [7] S. Perez, Identification and homogénéisation de paramètres dans des équations aux dérivées partielles Thèse, Université de Pau and des Pays de l’Adour, 1999.
444
R. Luce, S. Perez / Applied Numerical Mathematics 38 (2001) 427–444
[8] J.E. Roberts, J.M. Thomas, Mixed and Hybrid Methods, Handbook of Numerical Analysis, Vol. II, NorthHolland, Amsterdam. [9] J.M. Thomas, D. Trujillo, Mixed finite volume methods, Internat. J. Numer. Methods Eng. 46 (1999) 1351– 1366. [10] D. Trujillo, Couplage espace temps de schémas numériques en simulation de réservoir, Thèse, Université de Pau and des Pays de l’Adour, 1994. [11] J.E. Warren, H.S. Price, Flow in heterogeneous porous media, Soc. Petroleum Eng. J. (1961) 153–169.