A control volume method to solve an elliptic equation on a two-dimensional irregular mesh

Computer Methods in Applied Mechanics and Engineering 100 (1992) 275-290 North-Holland CMA 268

A control volume method to solve an elliptic equation on a two-dimensional irregular mesh I. Faille lnstitut Frangais du Pdtrole, B.P. 311,, 92506 Rueil Malmaison cedex France

Received 9 November 1991

The approximation of an elliptic partial differential operator on a two-dimensional irregular grid by a control volume method is studied. A conservative control volume scheme with a consistent approximation of the fluxes is shown to satisfy a property of weak consistency and to be convergent under appropriate stability assumptions. Numerical experiments are presented.

Introduction

The present paper deals with the approximation by a cell centered control volume method, of an elliptic partial differential operator such as div(~ grad u) on a 2D irregular meshing where gr is a discontinuous tensor. This problem arises in oil recovery [1] or basin modelling [2] for instance, when modelling 2D diphasic fluid flow in a heterogeneous porous medium. Indeed, in order to take into account the heterogeneities of the medium, the grid used to represent the medium follows the geological discontinuities and is therefore composed of irregular quadrangles. As the set of equations governing the flow (in particular conservation laws) couples unknowns of non-linear hyperbolic behaviour and parabolic behaviour, a control volume method is chosen to discretize the equations. Applying a control volume method to solve a conservation equation requires the following main steps: - the exact solution is approximated by a piecewise constant function, defined by its value in each cell; the value in one cell approximates the mean value of the solution in the cell. - the discrete problem is obtained by integrating the conservation equation over each grid block. Using Green's formula, the integral over one cell reduces to the sum of the fluxes over each edge of the cell. Difficulties in applying this method to approximate div(Yf grad u) are of two kinds: (1) The first is related to the irregularity of the meshing: convergence properties of control volume methods are usually stated for uniform grids, or irregular rectangular grids. (2) The second comes from the anisotropy and discontinuity of the tensor ~Tf.This discontinuity implies the discontinuity of the gradient of the solution, which has to be taken into account. Correspondence to: Dr. I. Faille, Institut Franqais du P6trole, B.P. 311, 92506 Rueil Malmaison cedex, France.

0045-7825/92/$05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved

1. Faille, A control volume method

276

Before going more deeply into the study, let us make a remark about another kind of control volume method. The method mentioned above is called a cell centered control volume method as opposed to a point centered control volume method. A point centered control volume method differs from a cell centered one in the way the grid used to approximate the solution is built: an initial grid allows us to define the points where the solution is to be approximated and a secondary one allows us to define the control volumes over which the equation is integrated. Recently, several papers have appeared in the literature dealing with such methods on irregular grids [3-5]. These methods work well with irregular meshing due to complex domain geometry. However, when the irregularity of the meshing comes from the description of heterogeneities, a point centered control volume method would lead to the integration of the equation over control volumes that overlap heterogeneous domains, which does not seem suitable. The paper is organized as follows. In the first part, we focus on the irregularity of the grid by looking at a more simple problem, i.e. the approximation of the Laplace operator. The usual finite difference criteria of consistency is not necessary for a control volume scheme to be convergent, thus a non-consistent scheme is not necessarily a bad scheme. A study about standard cell centered control volume schemes was carried out; it enhances two criteria that seem to be more appropriate to qualify a control volume scheme, These two criteria are a conservative scheme with a consistent approximation of the fluxes. Indeed a scheme satisfying these two criteria can be shown to satisfy a property of weak consistency which under stability assumptions implies the convergence of the scheme. A method of building such a scheme on a 2D irregular grid is given. In the second part, this scheme is extended to the approximation of div(~ g r a d u): the discontinuity of K is taken into account by writing the continuity of the fluxes on each cell boundary. Finally, in the third part numerical experiments are presented.

1, Approximation of the

Laplace

operator

The problem under consideration is -Au=f

inn,

u -- u~

on F ,

(1)

where D is an open bounded subset of R: and F its boundary. We first introduce the following notation: Mh = { Vk, k = 1 , . . . , K} is a meshing covering t~, h-is the grid size, i.e. a scalal ~such that there exist two constants c and d, independent of h, satisfying the inequality ch <~I <~dh for any length I of the grid (distance between two neighbouring cell centers or length of an edge). - % is the centre, i.e. the isobarycentre and S~ the surface of the cell Vk. - n~ is the normal to an edge 8 orientated in an arbitrary way. - n~ is the outward normal to the boundary of cell Vk on the edge 6, i.e. n~ -n~. - D(~f~) is the space of infinitely smooth functions with compact support i n / ] . - li is the exact solution of the problem under consideration.

!. Faille, A control volume method

277

u h is the approximated piecewise constant solution. - Uh is the vector whose components are the values of u h in each cell. - Uh is the vector whose components are the mean values cf ~ in each cell. - U h is the vector whose components are the values of t/at each cell center. Let us now approximate problem (1). Following the control volume rule, the values of the approximated function u h are a solution of -

± Sk ~ edge

f

L

o f Vk

jo

(-grad u •

.i)

=

f doJ = fk ,

(2)

where f~ ( - g r a d u. nk)do- has to be approximated with the values of u h. We denote this approximation by L ~ ( u h ) . If we assume that this approximation is linear, the discrete problem can be written as a linear system whose unknowns are the values of u h"

A,,U,, = Z,,.

(3)

The control volume method is defined as soon as the approximation L [ ( u h ) of (f8 grad un k do-) is specified. Several possibilities can be thought of for this approximation. Usually, in finite difference schemes, the approximation is chosen in order to get a consistent scheme: the finite difference truncation error Th, i.e. the vector whose component k is obtained by substituting the values of the appro~mate solution Uh by the values of the exact solution U h in the discrete equation, ( T h ) k = ( A h U h - B h ) k , is such that

IIThll

o

when h--, 0.

After a study about cell centered control volume methods on irregular 1D and 2D rectangular meshing, this criteria of consistency appears to be unnecessary for a control volume scheme to be convergent. The main results of this study are given in the first section. Consequently to this study, two criteria which seem to be more appropriate to qualify a control volume scheme have been enhanced. A scheme fulfilling these two criteria satisfies a property of weak consistency, which under stability assumptions leads to the convergence of the scheme. These results are detailed in the second paragraph. Then, in the last section, a method of building ~ control volume scheme satisfying the two criteria is suggested. 1.1. Usual schemes on 1 D a n d 2 D irregular m e s h i n g

On a 1D irregular grid, a usual control volume scheme is the three points scheme, which uses for the approximation of the flux on the extremity of a cell, the difference between the values of the approximated function in the two cells adjacent to that extremity, divided by the distance between the centers of these two cells. This scheme has been studied by several authors [6, 7]. They have proven that although it is not consistent, this scheme is convergent. Indeed, the part of the truncation error which does not converge towards zero with the grid size can be seen as the product of a vector Wh by the matrix A h, which reduces to zero with the grid size, i.e.

I. Faille, A control volume method

278

r , = T°+

7",',,

where T o has a good behaviour, i.e.

IIr 2 II

0

when h

---->

0

and T~ = A , I~h and II II ~ 0 when h --->0. Once such a decomposition ls obtained, the convergence of the scheme derives from the L ~ stability of the scheme (the inverse of matrix A h is bounded in L ~ norm independently of the grid size) and the following equality: U h - U, = A h ' T ° + W h .

These results were extended to the case of a 2D irregular rectangular grid, but on a general 2D irregular grid, no result on control volume schemes is available. However, on 2D grids composed of irregular quadrangles, the five points scheme, which uses the values of u h in the two adjacent cells for the approximation of the flux on an edge, leads to a conservative scheme but does not have a co~sistent approximation of the flux. It can be shown by studying this scheme on a particular irregular grid, that the computed solution does not converge towards the exact solution [9]. 1.2. New criteria for a control volume scheme

In control volume methods, it is more natural to speak in terms of piecewise constant functions rather than in terms of vectors. Therefore, we first introduce the following entities: E, is the set of functions, constant in each cell of M,. - For any function @ continuous on •, @, is the function of Eh associated with M, and O. Its value in one cell V, is Oh such that

-

-

-

~h is the linear operator from E, into Eh, defined by the matrix A,. The truncation error 3"h is the function of E h whose value in the cell Vk is (3"h)k =

(A,,Uh- B,) k.

1.2.1. The two criteria

Taking into account the results in the first section, two criteria appear as good criteria for a control volume scheme: (1) A conservative scheme. A control volume solution satisfies conservation over any cell because of the way it is obtained. In order to get a solution which satisfies conservation over any group of cells and so over the whole domain, a unique expression of fn grad u. n6 do" has to be used, i.e. the same expression (except for the sign) is used when writing conservation over the two cells adjacent to edge & If 8 is an edge of the cell Vk, we denote the approximation of j'~ (-grad u. ha)de" by L n ( u , ) and we have

!. Faille, A control volume method

279

This is a natural criteria for a control volume scheme, and from now on, we assume that a control volume scheme always satisfies this criteria. (2) A consistent approximation of the flux. The approximation chosen for the flux on each edge is such that

L,(~h) = -/,[grad @- n(co,) + as] V~ in C2(~),

(4)

where l, denotes the length of edge 3, ca8 its centre and a 8 a scalar such that [ a , l < ah, a is a constant independent of h. The importance of these two criteria is strengthened by the following propositions. The first is related to the consistency of a control volume scheme (conservative) with a consistent approximation of the fluxes. Although it is not consistent, such a scheme satisfies a property of weak consistency. Moreover, under stability assumptions, weak consistency implies convergence of the scheme, i.e. Proposition 2. Thus, Proposition 2 shows that a control volume scheme (conservative) with a consistent approximation of the fluxes is convergent under suitable stability hypothesis. For simplicity, the assumption u b = 0 is made in these two propositions.

PROPOSITION I. A cell centered control volume scheme (conservative) with a consistent approximation of the fluxes is such that its truncation error ~, satisfies the following: 3"h converges in L ®for the weak • topology towards zero fa 3"hO dca--*O when h--cO VO in L~(I~).

(5)

The scheme is said to satisfy the property of weak-consistency. PROOF. The value of the truncation error 3"h in the cell Vk is (6)

(ffh)k ffi (AhOi, -- Bh)k

1 (

m ~k

~

8edgeofVk

L~(uh))--fk.

(7)

The proof of Proposition 1 needs two steps: Step 1. We first notice that the truncation error is the sum of two terms:

where 1( = T,

X

l,[-gradu'nk(<°~)])-fk,

(9)

. od,o o, 2

(n~ . nk)lras .

(10)

The first term (3"h) is related to the control volume formulation and the use of Green's formula. It behaves well, i.e. II3"°11="-> 0 when h - > 0 . Indeed we have

I, Faille, A control voh~me method

280

/8[-grad t~. nk(~aa)] = f6 (-grad a. nk) d¢o+ O(h3),

(~)

and then 3 h satisfies

o =~ (~'h)~

(-Aa- f) doo) + O(h).

02)

Therefore (~h)k is a O(h) term and we have

119',II1=~0

when h--,0.

(13)

The second term ( ~ ) is related to the approximation of the flux on each edge. As a~ - O(h), [[~r], [[~, is bounded independently of h, i.e. there exists a constant C, independent of h, such that II~~ [[~ < c .

(14)

Step 2. We now prove that the function 3.~, given by (10) satisfies: ~'~--,0 in L ~ for the

weak • topology, i.e. f# 3"~,O dto---,0 when h--,0

¥4, in L J(gt).

(15)

As (W],) is bounded in L ~ (14), there exists a subsequence of (3.~,), still denoted by (3"~,), which converges towards a limit 3.J in L '~ for the weak • topology. Let us prove that this limit satisfies j'# 3" )O dto -- 0, for any • in D(I']). Multiplying the value of 3":, in each cell by the value of O~, (the function of Ej, associated with O and M~,) and summing over all the cells of Mp,, we obtain

k

= '~k

( S edgeZ of V'~(.~. ,,~)t,,~,)~,~.

(17)

As the scheme is conservative, a s depends only on edge 8 and not on the cell Vk. If V~ and Vk. k0 are the two adjacent cells to the edge 8, we have (n s • n~) --- - ( n s . n s ) and we can integrate by parts in (17). This leads to

f, ~,s'~ do, = ~ (n,. ,,)t~,(~P~- ~,.).

(18)

6

Finally, as Is, ~Pk - Ok, are O(h) terms and as the approximation of the flux is consistent, a# is

also an O(h) term, so (18) is reduced to

S

281

!. Faille, A control vohm~e method

The number of edges is 0 ( 1 / h " ) and therefore

f~ ~h •~o7"1 i, d~o---*O when h---~O.

(20)

Furthermore @~,--->@ in L '(g~) and ~07-1 ~,--->O"! in L ~ for the weak * topology, so

fl qbl,'~ o71h dto---->fl ~ ~$- I dto when h--->O.

(21)

Hence we deduce from (20) and (21) that

f @~"do=()V @ i n D ( n ) ,

(22)

0"~ = 0 .

(23)

and then

Then Proposition 1 is deduced from (13) and (23).

El

PROPOSITION 2. Consider ul, E El,, given by a control volume scheme (conservative) with a consistent approximation of the fluxes for problem (I). Under the following stability assumptions: (l ) there exists a constant C such that Ilu,, II= < c, (2) ul°--+u a.e. and u E H ~ ( I ] ) , (3) Jbr any function @ ~ D(I~), there exists constant C, such that IIM~,q~l,ll~ < C,, where @h is the function of E h associated with @ and M n, Ml* is ate adjoint operator of Mh for the L 2 scalar product in E h, then the limit u of the sequence (ul,) is the solution of problem (l). PROOF. Let us first establish the following iemma about M~*. LEMMA. Under the assumptions of Proposition 2, the subsequence (M~@h) is such that (M~, @h) ~ -- A ~ in L ~ ]'or the weak • topology when h --, O. PROOF OF LEMMA. Under the third assumption of Proposition 2, the sequence (M~,~h) is bounded in L ~ so there exists a subsequence, still denoted (M~,@h) which converges towards a limit L(@) in L ~ for the weak • topology. Let us prove that this limit satisfies

ft

g'L(@) dto = - ~ g' A@ do,

We consider a function I/' in have

D(O) and ~, the function of

fg, ~,,ff~ @h dto = f~ .fib ~,tOh dto.

(24)

for any 1/' • D ( O ) . Eh

associated with 1/' and

M h.

We

(25)

!. Faille, A control volume method

282

As the scheme is weakly consistent, we can write

fo

do-. fo

when h-->0.

(26)

Moreover, qrh --->qt in L l(/]) and (M~ ~,)---> L(W) in L ~ for the weak • topology so, from (25) and (26), we have

f,, ~L(,O)do =

-

fa A ~ ¢ d~o

(27)

= - (. qt A@ do~.

(28)

ja

We deduce from (28) that L ( @ ) = - A @ and then (a~,a,,)~-A¢, in L ~ for the weak • topology when h---, 0. I-1 Let us now go back to the proof of Proposition 2. The solution u h, given by the control volume scheme, satisfies A,U, = B h. So for any function • in D(.O), we have

f Mi,u,,@h

do = f~ fl, q~h do

(29)

f u,,a,*, q~h do = fs f~' @~'d o .

(30)

From the first and second hypothesis, we obtain u,---, u in L 102 ) and from the lemma, we obtain (A *J,Oh)--*-AO in L for the w e a k , topology. Then reducing h to zero in equality (30) gives

f u(-AC,)d,o=f, fcPdoo.

(3t)

Finally, as u E H~(/~), u satisfies

f. (-grad u grad O ) do - f. fO d o , and hence u is a solution of problem (1).

for any • in D(/~),

(32)

12

1.3. A control volume scheme ~,~e are going to build a conservative control volume scheme with a consistent approximation of the fluxes. Hence, we have to specify how a consistent approximation of the flux on each edge can be obtained. Let us consider an edge ~ (Fig. I). To obtain an approximation of the gradient of u in the

I. Faille, A control volume method

283

Fig. 1. Points used in the approximation of the flux on the edge 6.

direction normal to edge 6, the basic idea is to approximate u at two points H, B located on both sides of 6 and on the straight line A~, going through to~ in the direction orthogonal to 8. An expression for L~ is then given by

UH- Un

L~ ffi - 1~ d(H, B) '

(33)

where U n, U~ are approximations of u at points H and B and d(H, B) is the distance between H and B. In order to obtain a simple expression for U n, H is chosen on the straight line joining the centers of two adjacent cells and then a linear combination between the values of u~, in these two cells gives an expression for U H. For instance, in the case of Fig. 1, H is the intersection between line A~ and line (M, M') and U~s = a UM + (1 - ~)UM,,

(34)

where a ffi d(H, M')/d(M, M') and UM, UM. are the values of u~, in the two cells M and M'. B and Utj are obtained in a similar way. Hence, the expression of the flux on an edge is a linear combination between four of the six values of u h in its neighbouring cells, the choice of the cells depends on the direction of the edge. When (2) is written for a particular cell, we obtain a relation between at most nine values of u h (the value in the considered cell and the values in the eight neighbouring cells) so the scheme is a nine points one. In the sequel, we call it VF9. The approximation of the flux on an edge located on the boundary of the domain, requires a particular treatment. This is presented in Appendix A.

1.3.1. Properties of the scheme The scheme satisfies the hypothesis of Proposition 1 so it satisfies the property of weak consistency. Concerning its stability, it can be established that under the assumption that the grid is not too much distorted, the matrix A h is irreducibly diagonally dominant with positive real entries on the diagonal and negative real entries outside the diagonal [9]. Therefore, the

I, Faille, A comul volume method

284

discrete solution satisfies a discrete maximum principle and the inverse of A,, is bounded independently of h. There exists a constant C > 0, independent of h, such that

More stability results are needed to be under the hypothesis of Proposition 2, however numerical experiments presented in the third part show that the scheme behaves well. 2. Approximation of div(X grad u) The problem under consideration is now

-div(X grad u) = f U =U

h

in 0

,

135)

onr.

Following the method presented in the first section, the construction of a scheme is reduced to the approximation &(u,,) of the fiux on each edge:

IR(-XgradEc~n,)de. As the medium is heter eneous, X can hav different values the two cells adjacent to 6. This discontinuity impli a discontinuity of rad u on the ed , the continuous quantity is ad u e nA, In order to take this discontinui o account, we proceed in the followin we approximate the flux on both sides of 6 (this can be done as soon as U *, th value of u at the middle of e 6 is introduced), then we use the continuity to eliminate 1:’ between the two expressis the Wuxto obtain a unique expression, independent of UY recisely, the followi ps are required: aider tP an appro ion of u at the middle CL)~ of e this value to approximate X rsrd u 9 n, in the two cells nt to 6. Let us consider, for instance, a vertical edge 6. If XR is the permeability tensor in the right cell and SC; its transpose, we have the following equality: = grad u (X’,n,) . l

Then an expression for the approximation of L, on the right side of 6, Lt9 is (36)

where 14is a point leciated on the straight line d,, going through o4 and in the direction (X’,n,) and K, = 11 Kit& /12.As in the case of the Laplace operator, R is chosen on the straight line joining the centres of two adjacent cells (Fig. 2) and UR is a linear

!. Faille, A control volume method

285

AR L~

Fig. 2. Points used in the approximation of the flux on the edge 6.

approximation between the values of u at these two cell centres. In a similar way, an expression for L L (left side) is obtained. (3) Write the continuity of Y/"grad u.n~ on the edge, to obtain an expression for U*, i.e.

LI'=L ,

[ K.

KL

R)

+

d(~%.L)

]

KR

U*=

d(~.

KL

R) U. +

d(~-, L)

UL"

(4) Substitute this expression of U* in one of the two expressions of the flux to obtain the final expression, KL KR [ KL KR ] -~ L s = -l~ d(toa, L) d ( ~ R') d ( ~ L) + d ( ~ . R) (Ua

- -

UL).

(37)

When the cells are rectangular and when the tensor ~/" is diagonal in the reference of the coordinate axes, AL and AR are simply the line orthogonal to the edge 8. The scheme then reduces to the expression of the usual five points scheme with harmonic average of the permeabilities [ 1]. 3. Numerical results

In this section, we present the results of some numerical experiments. The first three illustrate different aspects of the control volume method for the Laplace equation. The interest of Test 1 lies in the grid used as it is very distorted; in the second test the control volume scheme is compared to a bilinear finite element method and the third test shows some of the convergence properties of VF9 by solving the same problem on a coarse and on a refined grid. In the last test, the scheme i:: applied to the resolution of div(~/" grad u). TEST I: gershaw's test. This test was proposed by Kershaw and concerns the resolution on a very distorted grid of the following problem:

i. Faille, A control volume method

286 ? ....

~ ....

,4 . . . .

~ ....

? ....

0l l

l,o

~

, . . I , , . , I

4 . . . . 6I ,

8

, . , 1 . , , , 1

10

l

ll

"i

1

] 10 Fig. 4. Test 1 - isovalue curves.

Fig. 3. Test I - V F 9 solution. Each cell is coloured in accordance with the value of the approximated function in that cell.

- A u = 0 in [0, 10] x [0, 10],

u(x, o) = o,

(38)

u(x, 10)-- 10, Ouu (0, y)= Ou

ox

~ (1o, y) = o.

The exact solution is t~(x, y) -~ y. Figure 3 gives the solution u~, computed by VF9 whereas Fig. 4 gives the isovalue curves of that solution. It can be noticed that these curves are not altered by the distortion of the grid as they would be with a five points scheme (see Section 1.1). TEST 2: Comparison between VF9 and a finite element method. The problem

-Au - 0

in [0, 100l × [0, 10],

u(x, o) = u(O, y)= u(1oo, y)= 1oo, u(x, lo)= ~2500 (1 - c o s ( 2 - ~ ) )

(39)

+ 100

is solved with VF9 (Fig. 5) and with a bilinear finite element method. Since the exact solution is known [8], the relative errors between the approximated solution and the exact solution have been computed for VF9 (Fig. 6) and for the finite element method (Fig. 7). We can observe that these errors are of the same order, although they are distributed differently (let us recall that the finite element method gives the exact solution on the boundary of the domain as the boundary conditions are of Dirichlet type).

!. Faille, A control volume method

287

nnumnun!llnllnnnnnnu I iungiii iiii!nnnnu n '-~::iii :':~~lll iiiiiil

liniiii,~=: III ili

' ililll

| n i i i iiili

i.iii ! ! l l

unnu ii;i

iilmm

+

980 870• 760nn 6 5 0 N 540I 430mN 3 2 0 210-

--I-

1090 1090 980 870 760 650 540 430 320 210

Fig. 5. Test 2 - VF9 solution.

I. l i A . I . . I

. I.,

..

I....

I....I...*1

....

I*.*

.1.*..

I*''LI

II U

+ 0.17 0.15 - 0.17 0.13 - 0.15 0.11 - 0.13 0.09 - 0.11

6

...... 0.07 0.05 0.03 0.01 ~ -

B

8 |0

- 0.09 -- 0.07 - 0.05 - 0.03 0.01

Fig. 6. Test 2 - relative errors for VF9. The colour o f each cell corresponds to the

Z

value of the relative error at the

cell centre.

.l_[j,

, I = t=LI.,

,= I,,,,

I,

, , ,i

,,

,,l,,

, , I ....

I,,,,L,,,,I

0

+ 0.15 0.13 0.11 0.090.07 0.05-

2

4 $

0.17 0.17 0.15 0.13 0.11

0.09

0.07 0 . 0 3 - 0.05 0 . o , - 0.03 0.01

8 tO

Fig. 7. Test 2 - relative errors for the finite element method. The colour of one cell corresponds to the mean of the relative errors at the four nodes of the cell.

T E S T 3. The problem under consideration is

- A u = - 2 [ x ( x - 10) + y ( y - 10)1

in [0, 10] x [0, 10], (40)

u(x, O) = u(x, 10) = u(0, y) = u(10, y) = 0.

288

!. Faille, A control volume method ,o ....

o .... 2 .... 4 .... ~ .... ~ .... to

i ~

~

+ 3.2 2.8 2.4 2.01.6-1.20.80.4 -

,2 ....

4 ....

e, ....

9....t° +

3.6 3.6 3.2 2.8 2.4 2.0 1.6 1.2 0.8 0.4

3.2

5.6 3.6

-

i I I

2 . 8 - 3.2 2 . 4 - 2.8 2 . 0 - 2.4 1.6 -2.o

!~

1.2 o.8

I

0.4

-

g

1.6

-

1.2

0.8 0.4

-

I--I-

10

Fig. 8. Test 3 - relative errors on the coarse grid.

Fig. 9. Test 3 - relative error on the refined grid.

It is solved on a coarse and on a refined grid, obtained when dividing the edges of the coarse grid by two. The relative errors between the VF9 solution and the exact solution have been computed on the coarse grid (Fig. 8) and on the refined grid (Fig. 9). These errors are much smaller for the refined grid, in fact they are almost divided by four.

TEST 4. In this test VF9 is used tO solve div(A" grad u) on a 2D domain. This problem comes from basin modelling: a monophasic fluid flow is modelled in a 2D porous medium composed of lithologic layers (Fig. 11). The grid used to represent the medium follows these layers and is made of irregular trapeziums. The tensor Y/' can be seen as a permeability tensor, it is diagonal in a orthonormal reference (s, a) related to the geometry of the layers (Fig. 10). The problem solved is

=div(~ grad u) = 0 in [0, I0] x [0, I0],

(41)

u(x. O) = O,

u(x, 10)= 100, no flux on lateral boundaries,

+

90

80 70 -

J

60

-

90 80 70

50

-

60 50

i

s...

30 - 40 20 -

30

20 10

10

Fig, 10, Reference (s, a) in a layer,

Fig. 11, Test 4 - VF9 solution.

!. Faille, A control volume method

289

Fig. 12. Points used in the approximation of the flux on the edge &

where the values of Yf, K s and Lo in the reference (s, a) are K s = 1 , K a = 10 -3 in the first four rows and the last three rows. - K s = 1, K,, = 10 -5 in the three other rows. Figure 11 gives the VF9 solution. It reflects well the main features of the permeability tensor: - anisotropy: the isovalue curves follow the stratigraphic layers as K s is greater than K~. - discontinuity: the gradient of the solution is more important in the three rows where K a is lower.

Appendix A. Approximation of the flux on the boundary of the domain Let us consider an edge 8 located on the boundary of the domain. We have to approximate the gradient of u in the direction given by z~, the straight line going through to~ and orthogonal to 8. The boundary condition u = u h on F gives us a value u~ of u at the centre of edge 3: u~ -- u~,(to~). To approximate grad u. n, it is sufficient to approximate u at a point B located on a~. Indeed, an expression for L~ is then given by --

//,s - / ¢ / t

B)' As in the case of an interior e d g e , point B is chosen on the straight line joining the centres of two adjacent cells in order to obtain a simple expression for un; u n is a linear combination between the values of u~, in the two adjacent cells.

References [1] K. Aziz and A. Settari, Petroleum Reservoir Simulation (Elsevier, London, 1979). [2] P. Ungerer, J. Burrus, B. Doligez, P.Y. Chenet and F. Bessis, Basin evaluation by integrated two dimensional modelling of heat transfer, fluid flow. hydrocarbon generation and migration, AAPG Bull. 14 (1990) 309-335. [3] Z. Cai, J. Mandel and S. McCormick, The finite volume element method for diffusion equations on general triangulations, SIAM J. Numer. Anal., to appear. [4] P.A, Forsyth, A control finite element method for local mesh refinement, SPE 18415, February 1989. [5] B. Rozon, A generalized finite volume discretization mt,thod for reservoir simulation, SPE 1814, February 1989. [6] P.A. Forsyth Jr. and P.H. Sammon, Quadratic convergence for cell centered grids, Appl. Numer. Math. (1988) 377-394.

290

!. Faille, A control volume method

[7] T.A. Manteuffei and A.B. White, The numerical solution of second order boundary value problem on non uniform meshes, Math. Comput. 47 (1986) 511-536. [8] H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, 2nd Edition (Clarendon Press, Oxford). [9] I. Faille, Les volumes finis en vue de la modelisation bidimensionnelle de la genese et la migration des hydrocarbures, Rapport IFP no. 37763, 1990. [10] D.S. Kershaw, Differencing of the diffusion equation in Lagrangian hydrodynamic codes, J. Comput. Phys. 39 (1987) 375-395.