Approximation of a Second-order Elliptic Equation with Discontinuous and Highly Oscillating Coefficients by Finite Volume Methods

In this paper we consider the numerical approximation of a class of second order elliptic boundary value problems with discontinuous and highly periodically oscillating coefficients. We apply both classical and modified finite volume methods for the approximate solution of this problem. Error estimates depending on ε the parameter involved in the periodic homogenization are established. Numerical simulations for one-dimensional problem confirm the theorical results and also show that the modified scheme has a smaller constant of convergence than the classical scheme based on harmonic averaging for this class of equations.


Introduction
There are many practical computational problems with highly oscillatory solutions e.g. computation of flow in heterogeneous porous media for petroleum and groundwater reservoir simulation (see, e.g., (Hornung, 1997) and the bibliographies therein). If a porous medium with a periodic structure is considered, with the size of the period is small enough compared to the size of the reservoir, and denoting their ratio by ε (0 < ε << 1) an asymptotic analysis, as ε −→ 0, is required. In this paper we will consider problems that are described by a linear elliptic equation in divergence form with highly periodically oscillating coefficients. Especially we will consider the following model problem: Ω ⊂ R n (n = 1, 2, 3) is a bounded polygonal convex domain with a periodic structure and smooth boundary Γ, K ε (x) = K(x /ε), K is a symmetric and uniformly positive definite matrix in Ω which has jumps discontinuities across a given interface. The case of piecewise constant coefficient K ε is very important for the applications. In porous medium flow, the problem (P ε ) results from Darcy's law and continuity for a single phase, incompressible flow through a horizontal heterogeneous porous medium with periodic structure.
Whenever effective equations are applicable they are very useful for computational purposes. There are however many situations for which ε is not sufficiently small so that the effective equations are not practical. In this cases the original equation has to be approximated directly. The numerical approximation partial differential equations with highly oscillating coefficients has been a problem of interest for many years and many methods have been developed (see, e.g., ( Amaziane & Ondami, 1999), (Chen & Hou, 2002), (Versieux & Sarkis, 2008) and (Ondami, 2001(Ondami, , 2015 and the bibliographies therein).
The case where K ε has continuous coeffients is the most studied. The case of discontinuous coefficients has been addressed in some work as in (Bourgat, 1978) and (Bourgat & Dervieux, 1978) where Authors use a double-scale asymptotic expansion. In (Amaziane & Ondami, 1999) and (Ondami, 2001), the numerical approximation of the problem, in the case of discontinuous coefficients was done by finite elements methods. However no error estimate has been established and is still, to our knowledge, an opened issue.
In this paper the approximation will be done by finite volume methods (see, e.g., (Eymard, Gallouët & Herbin, 2000) and (Chernogorova, Ewing, Iliev & Lazarov, 2001)) and the study is limited to one-dimensional problem (1-D problem). This 1-D problem illustrates very clearly the dependence of numerical results to ε, the parameter of homogenization.
Error estimates are established. The obtained results can be generalized in the two-dimensional and three-dimensional problems. The paper is organized as follows. In section 2, a description of methods used is presented as well as the obtained error estimates. Section 3 is devoted to numerical simulations. Lastly, some concluding remarks are presented in section 4.

Methods
Our study will focus on the one-dimensional problem and we assume (without loss of generality) that Ω =]0, 1[. In this case the problem (P ε ) is written simply (1) where k ε (x) = k (x /ε) = k(y), with y = x /ε, k is a discontinuous and periodic function of period 1, on ]0, 1[ . In all this paper we make the following assumptions The assumptions (A1) and (A2) ensure the existence and uniqueness of the solution of the problem (1). From homogenization theory (see, e.g., (Sanchez-Palancia, 1980), (Bensoussan, Lions & Papanicolaous, 1987) and (Jikov, Kozlov & Oleinik, 1994)) follows u ε ⇀ u in H 1 0 (Ω) (consists of functions in Sobolev space H 1 (Ω) that vanish on 0 and 1) weakly, where u (homogenized solution) satisfies the following homogenized problem, and the constant k * is the mean harmonic value of k(y) on (0, 1).
This one-dimensional problem helps to clarify the eventual dependency of numerical results to the parameter ε. The error estimates obtained can be generalized in the two-dimensional and three-dimensional cases which uses, for instance, simplices or parallelepipedes mesh. Two finite volume schemes will be compared: The classical scheme (see, e.g., (Eymard, Gallouët & Herbin, 2000)) and a modified scheme (see, (Ewing, Iliev & Lazarov, 2001)).
In order to compute a numerical approximation to the solution u ε , let us define a mesh, denoted by T , of the interval (0, 1) consisting of N cells (or control volumes), denoted by V i , i = 1, ..., N, and N points of (0, 1), denoted by x i , i = 1, ..., N, satisfying the following assumptions: Definition 1. An admissible mesh of (0, 1), denoted by T , is given by a family , and a family ( x i ) i=0,...,N , such that

Classical Finite Volume Scheme
..,N be an admissible mesh, in the sense of Definition 1 , such that the discontinuities of k ε coincide with the interfaces of the mesh. Classical finite volume scheme consiste to integrate the first equation of the problem (1) over V i (see,e.g., (Eymard, Gallouët & Herbin, 2000)). So we have: ..,N be the discrete unknows and let In order that the scheme be conservative, the discretization of the flux −k ε (x) d u ε (x) d x at x i+ 1 2 should have the same value on V i and V i+1 . To this purpose, one introduces the auxiliary unknown u ε d x may be performe on each side of x i+ 1 2 by using the finite difference principe: with u ε 1/2 = 0 and u ε N+1/2 = 0, for the boundary conditions. Requiring the two above approximation of to be equal (conservativity of the flux) yields the value of u ε which, in turn, allows to give expression of the approximation H ε with τ ε and therefore the mean harmonic value of k ε is involved.
Remark 2. The fact that k ε is discontinuous, periodic (period ε) and with discontinuities that coincide with the interfaces of the mesh T leads to say that: Note. Throughout the paper, we will denote by c generic constants, even if they take different values at different places.
We now state the main result of this section.
..,N , and the discontinuities of k ε coincide with the interfaces of the mesh, ..,N is the solution of (4)- (8). Then there exists a constant c independent of ε and h such that where τ ε i+ 1 2 is defined in (7), and Proof.
The proof of the Theorem 1 is obtained by using the same gait as in (Eymard, Gallouët, & Herbin, 2000). Let and H * ,ε Let us first show there exists c independent of ε and h such that In order to show this, let us introduce . By developing the first equation of (1) in V i (i = 1, ..., N), one obtains According to (Amaziane & Ondami, 1999) and (Ondami, 2001), one has and By using (A1) and (14)- (17), one deduces that there exists c independent of ε and h such that This yields (13) for i = 0 and i = N.
The following equality: yields that where S ε Using (20), this implies that H * ,ε Using (18) and (19), this last inequality yields that there exists c, independent of ε and h such that Therefore (13) is proved.
Remark that Applying once again the Cauchy-Schwarz inequatity, one obtains e ε i AB, which yields Estimation (10).

Modified Finite Volume Scheme
The modified finite volume approach (see, Ewing, Iliev & Lazarov, 2001) is to rewrite the problem (1) into its mixed form: q ε (x) is the flux dependent variable. Conditions for continuity of the fonction and the flux through interface points ξ are added: Here [u ε ] deontes the difference of the right and left limits of u ε at the point of discontinuity. We introduce a standard uniform cell-centered grid x 0 = 0, Note, that the endpoints x 0 = 0 and x N+1 = 1 are part of the grid, but they are at h 2 distance from their neighboring grid points.The internal grid points can be considered as centered around the volumes V i = The values of a funtion f defined at the grid points x i are denoted by f i . Non-uniform grids can be treated in a similar way. The finite volume method exploits the idea of writing the balance equation over the finite volume V i , i.e. integrating the first equation of Problem (1) over each volume V i .
Next, we rewrite the flux equation in the form and integrate this expression over the interval (x i , x i+1 ) : One assumes that the flux is two times continuously differentiable on the interface, so it can be expanded around the point x i+ 1 2 in th Taylor series After replacing the first derivative of the flux at x i+ 1 2 by a two-point backward difference one gets the following approximation of (30).
Finally, by the same gait as in (Ewing, Iliev & Lazarov, 2001) we get the following finite difference approximation of the differential problem (27): where where k εh and u εh i denotes the approximation values of the exact solution. k εh i+ 1 2 is the well known harmonic averaging of the coefficient K ε (x) over the cell (x i , x i+1 ) , which has played a fundamental role in deriving accurate schemes for discontinuous coefficients (see e.g, (Samarskii, 1977) and (Samarskii & Andréev, 1978 )).
We now state the main result of this section.
Theorem 2. Assume that the coefficient k ε (x) is a piecewise C 1 −function and has a finite number of jump discontinuities, the grid is such that the discontinuities are at the points x i+ 1 2 , and the source term f (x) is a C 1 −function on (0, 1). Then the following estimate is valid: where k εh i+ 1 2 is given by (35), and c is a constant independent of ε and h.
Proof. The proof of Theorem 2 is the same as that of Theorem 1.

Numerical Results
In this section, one presents numerical results, comparing the approximations described in this paper and an example of exact solution. More especially, we shall present numerical results obtained with following data of Problem (1): where n p is a positive integer; 0 < p < n p − 1 and the source function is f = 1. Therefore the exact solution u ε is 2k 2 + x 2k 2 + (k 1 −k 2 )εx 4k 2 (k 1 +k 2 ) − (k 1 −k 2 )(p+1)ε 2 4k 2 (k 1 +k 2 ) − (k 1 −k 2 )(p+1)ε 4k 1 k 2 and the homogenized solution is All the simulations presented have been done with uniform grids.

Concluding remarks
The purpose of this paper was to resolve, by finite volume methods, a class of second-order elliptic problems, with discontinuous and highly oscillating coefficients, in order to evaluate the effect of ε the parameter involved in the periodic homogenization on the approximate solution. Study was limited to the academic Dirichlet problem and enabled to clarify the dependence of numerical experiments to ε. Error estimates were obtained. Numerical simulations confirm theorical results and also show that the modified finite volume scheme is much more accurate than the classical scheme in solving these problems. The extension of these results to two-dimensional and three-dimensional problems is currently underway.