E Khan

The role of superfluidity in nuclear incompressibilities

Keywords: numbers: 21.10.Re, 24.30.Cz, 21.60.Jz

Nuclei are propitious tools to investigate the role of the superfluidity in the compressibility of a Fermionic system. The centroid of the Giant Monopole Resonance (GMR) in Tin isotopes is predicted using a constrained Hartree-Fock Bogoliubov approach, ensuring a full self-consistent treatment. Superfluidity is found to favour the compressibitily of nuclei. Pairing correlations explain why doubly magic nuclei such as 208 Pb are stiffer compared to open-shell nuclei. Fully self-consistent predictions of the GMR on an isotopic chain should be the way to microscopically extract both the incompressibility and the density dependence of a given energy functional. The macroscopic extraction of Ksym, the asymmetry incompressibility, is questioned. Investigations of the GMR in unstable nuclei are called for. Pairing gap dependence of the nuclear matter incompressibility should also be investigated.

The effect of the superfluidity on the compressibility of a Fermionic system remains an open question. Superfluidity initially referred to a system with a dramatic drop of its viscosity [1]: it could be suspected that a super-fluid would be easier to compress than a normal fluid. This question has been investigated in Fermionic atoms traps, by studying the frequency of the compression mode with respect to the scattering length. Experimentally some increase of the frequency may be observed in the weak pairing regime [2], but this signal remains to be confirmed. Theoretically both microscopic and hydrodynamical investigations show no variation of the compression mode between the normal and the superfluid phases [3], but the analysis is complicated by the temperature change between the two phases. In nuclear physics, the study of the role of superfluidity in the compressibility can also be performed: the isoscalar Giant Monopole Resonance (GMR) is a compression mode, allowing to probe for related superfluid effects. Ideal tools are especially isotopic chains, where pairing effects are evolving from normal (doubly-magic) nuclei to superfluid (open-shell) ones [START_REF] Bohr | Nuclear Structure[END_REF]. Moreover, the incompressibility of nuclear matter is a basic parameter in calculations describing neutron stars or supernovae, where superfluid effects are known to occur [START_REF] Pines | [END_REF].

Constraining the nuclear incompressibility modulus K ∞ with experimental data on the Giant Monopole Resonance is a longstanding problem. The first relevant approaches have been performed following the work of Blaizot and Pearson [6,7]: only microscopic predictions of the GMR compared with the data could validate the K ∞ value of the functional which was used. Using a fully self-consistent approach such as the Random Phase Approximation (RPA) [START_REF] Ring | The nuclear many-body problem[END_REF] or the constrained Hartree-Fock (CHF) [START_REF] Bohigas | [END_REF] method, the GMR data on 208 Pb currently provides K ∞ ≃ 230 MeV in non-relativistic approaches [10] whereas K ∞ ≃ 260 MeV is obtained for relativistic one [11]. The roots of this puzzle between the two approaches is still an open question, but it has been shown that the neutron-proton asymmetry dependence of the incompressibility, denoted as K sym , plays also a role, as well as the density dependence of the functional [12]: a value of K ∞ ≃ 250 MeV could be extracted from 208 Pb data, using a non-relativistic CHF method with a modified density dependence of the functional. Therefore, it should be noted that K ∞ cannot be extracted from the measurement of the GMR in a single nucleus: several parts of the functionnal are tested simultaneously, namely K ∞ , K sym and its density dependence.

A previous study on the role of superfluidity on nuclear incompressibility has been performed, finding a negligible effect [13], but the theoretical approach was not self-consistent. Indeed self-consistency is crucial since pairing effects are expected to be small in the GMR: this high energy mode is mainly built from particle-hole configurations located far from the Fermi level, where pairing do not play a major role. However giant resonances are known to be very collective [START_REF] Harakeh | Giant Resonances[END_REF] and pairing can still have a sizable effect on the GMR properties: around 10 % on the centroid, which is the level of accuracy of present analysis on the extraction of K ∞ [10,11]. This requires the advent of accurate microscopic models in the pairing channel, such as fully self-consistent Quasiparticle Random Phase Approximation (QRPA) [START_REF] Li | [END_REF]16], achieved only recently. Experimentally, the measurement of the GMR on an isotopic chain facilitates the study of superfluidity on the GMR properties [17], and the possibility to measure the GMR in unstable nuclei emphasizes this feature [18].

It is therefore necessary to go towards the measurement of the GMR on several nuclei, such as an isotopic chain. The overused method of precise GMR measurements in a single nucleus, such as 208 Pb, may not be the relevant approach. Other nuclei have been used such as 90 Zr and 144 Sm. Indeed when considering the available GMR data from which the K ∞ value has been extracted, 208 Pb is stiffer than the Sn, Zr and Sm nuclei: K ∞ is about 20 MeV larger, both in non-relativistic and in relativistic approaches [10,11]. The question may not be "why Tin are so soft

in2p3-00386181, version 1 -20 May 2009

?" [19] but rather "why 208 Pb is so stiff ?".

Recently the GMR was measured on the stable Tin istopic chain (from 112 Sn to 124 Sn) [17]. Once again it has been noticed that it is not possible to describe the GMR both in Sn and in Pb with the same functional, Tin beeing softer than Pb [START_REF] Li | [END_REF]19]. In the non relativistic case, fully self-consistent QRPA calculations on Sn isotopes lead to K ∞ ≃ 215 MeV. The relativistic DDME2 paremeterisation using QRPA describes well the GMR in the Sn isotopes [START_REF] Lalazissis | Proc. of Nuclear Structure and Dynamics[END_REF], but predict a low value of the GMR in 208 Pb compared to the experiment, as can be seen on Fig. 8 of Ref. [11]. On the contrary, a recent relativistic functional describes well the 208 Pb GMR, but systematically overestimates the Sn GMR values of about 1 MeV [START_REF] Niksic | [END_REF]. Finally, attempts have been performed in order to describe Sn GMR data with relativistic fuctionals having a lower incompressibility and hence different density dependence and K sym value [22]. Once again the 208 Pb and Tin GMR cannot be described at the same time: this puzzling situation is due to the higher value of K ∞ extracted from 208 Pb, compared to Tin, Sm and Zr nuclei.

In Ref. [START_REF] Li | [END_REF] it has been found that including pairing effects in the description of the GMR allows to explain part of the Sn softness : pairing may decrease the predicted centroid of the GMR of few hundreds of keV, located at ∼ 16 MeV. This is sufficient to change by about 10 MeV the extracted value of the incompressibility of nuclei K A , defined as [6]:

E GMR = 2 K A m r 2 (1)
where m is the nucleon mass and r 2 is the ground-state mean square radius.

In this work we follow up this idea and show that the consequences of superfluidity on nuclear incompressibility may solve the above mentioned puzzle. It should be noted that pairing is vanishing in the doubly magic 208 Pb nucleus, unlikely the other nuclei. It is necessary to use a fully microsocopic method including an accurate pairing approach. In order to predict the GMR in a microscopic way we use the constrained HF method, extended to the full Bogoliubov pairing treatment (CHFB). The CHF(B) method has the advantage to very precisely predict the centroid of the GMR using the m -1 sumrule [START_REF] Bohigas | [END_REF]12]. The whole residual interaction (including spin-orbit and Coulomb terms) is taken into account and this method is by construction the best to predict the GMR centroid [12]. Introducing the monopole operator as a constraint, the m -1 value is obtained from the derivative of the mean value of this operator. The m 1 sumrule is extracted from the usual double commutator, using the Thouless theorem [23]. Finally the GMR centroid is given by E GMR = m 1 /m -1 . All details on the CHF method can be found in [START_REF] Bohigas | [END_REF]10].

The extension of the CH method to the CHFB case uses the HFB approach in coordinate space [24]. Skyrme functionals and a zero-range surface pairing interaction are used in this work. The magnitude of the pairing interaction is adjusted so to describe the trend of the neutron pairing gap in Tin isotopes. This interaction is known to describe a large variety of pairing effects in nuclei [25].

Fig. 1 displays the GMR energy obtained from the Sn measurements (times A 1/3 to correct for the slow lowering of the GMR with the nuclear mass [START_REF] Harakeh | Giant Resonances[END_REF]). Microscopic CHFB predictions using two functionals are also shown: SLy4 [26] (K ∞ =230 MeV, which describes well the Pb GMR data), and SkM* [27] (K ∞ =215 MeV). Without pairing, the SLy4 interaction overestimates the Sn GMR data. Pairing effects (CHFB calculations) decrease the centroid of the GMR, getting closer to the data. This confirms the results of [START_REF] Li | [END_REF], where a self-consistent HFB+QRPA approach was used to describe the Tin data. It should be noted that using a less general BCS approach for the pairing channel increases the GMR energy [13,[START_REF] Li | [END_REF]28]. This is due to the problematic treatment of high energy single-particle states in the BCS model. It is therefore important to use the full HFB approach to correctly describe pairing effects on the GMR, especially for the pairing residual interaction. Fig. 1 shows, as expected, that SkM* predictions with pairing effects are in better agreement with the Sn data.

The striking feature of Fig. 1 is the peak of the GMR centroid, located at the doubly magic 132 Sn nucleus, using the CHFB predictions. This indicates that pairing effects should be considered to describe the behaviour of nuclear incompressibility, and that vanishing of pairing make the nuclei stiffer to compress, confirming our previous statement on the stiffness of 208 Pb. The importance of pairing effect can be understood in a simple way: since the nuclear incompressibility is defined as the second derivative of the energy functional at saturation density [6], there is no obvious reason why the pairing terms of the functional would play no role in the nuclear incompressibility. Nuclear incompressibility is indeed very close from a residual interaction (as a second derivative of the energy functionals with respect to the density), and it is known that pairing effects are relevant in residual interactions [25]. This is straightforward in nuclear matter where K ∞ is expressed from the F 0 Landau parameter [29]. However, on Fig. 1, the GMR centroid is shifted to lower energies for more neutron-rich nuclei than 132 Sn, also because of the appearance of a soft L=0 mode, predicted in QRPA calculations beyond 132 Sn [30].

To further investigate the role of pairing on nuclear incompressibility, Fig. 2 displays K A (defined by Eq. ( 1)) with respect to the average pairing gap calculated using the HFB approach, from 112 Sn to 132 Sn. A clear correlation FIG. 1: Top: Excitations energies of the GMR in stable 112-124 Sn isotopes calculated with constrained HF and constrained HFB methods and the SLy4 interaction, compared to the data. Bottom: Same for 112-136 Sn isotopes, using the CHFB method and the SLy4 and SkM* interactions is observed: the more superfluid the nuclei, the lower the incompressibility. Hence it may be easier to compress superfluid nuclei. This may be the first evidence of the role of superfluidity on the compressibility of a Fermionic system. A possible interpretation is that Cooper pairs can modify bulk properties, as known from nuclear physics phenomenology [START_REF] Bohr | Nuclear Structure[END_REF].

The decrease of incompressibility in superfluid nuclei raises the question of a similar effect in infinite nuclear matter: for now, incompressibility is given independently from the pairing part of the functional. However, considering present results, equations of state used for neutron star and supernovae predictions should take into account pairing to provide their incompressibility value. The comparison with GMR data shows, as mentioned above, that the functional as a whole (including pairing effects) is probed. The question of the behaviour of K ∞ with respect to the pairing gap is raised: it seems clear from nuclear data that nuclei incompressibility K A decreases with increasing pairing gap. This should be investigated in nuclear matter.

It should be noted that recent attempts have been performed to extract the K sym value, and its corresponding quantity in nuclei (K τ ): the incompressibility in the Sn isotopic chain has been studied using the macroscopic (liquid drop) formula of nuclei incompressibility K A derived by Blaizot [6,17,31]. However the effect of pairing demonstrated above shows that the current macroscopic approach may not be well designed: on Fig. 2, pairing effects induce ∼ 10 MeV change on the nuclear incompressibility K A . Hence the macroscopic expression of nuclei incompressibility K A should be extended to these terms, since the appropriate definition of nuclei incompressibility is the second derivative of the microscopic energy functional [6]. Presently, the microscopic approach is more relevant to extract K sym : the extraction of K ∞ , K sym and the density dependence of the functional are related, as stated in Ref. [12]. The GMR data on isotopic chains should be used on this purpose.

It is not possible to describe the GMR centroid of both 208 Pb and other nuclei with the same functional, as stated above. The puzzle of the stiffness of 208 Pb may come from its doubly magic behaviour. In Fig. 1 there is a sharp peak at doubly magic 132 Sn, and it would be very interesting to measure the GMR in this unstable nucleus. It should be noted that such experiments are now feasible [18]. A possible explanation of the 208 Pb stiffness is that the experimental data of E GMR is especially increased in the case of doubly magic nuclei, as observed in 208 Pb compared to the GMR data available in other nuclei (such as the Tin isotopic chain). This difficulty to describe with a single functional both doubly magic and other nuclei has already been observed on the masses, namely the so-called "mutually enhancement magicity" (MEM), described in [32,33]: functionals designed to describe masses of open-shell nuclei cannot predict the masses of doubly magic nuclei such as 132 Sn and 208 Pb, which are systematically more bound that predicted. In order to consider MEM, it may be necessary to take into account quadrupole correlation effects due to the flatness of the potentials for open-shell nuclei [34]. K A being related to the second derivative of the energy with respect to the density, it would be useful to find a way to predict the GMR beyond QRPA by taking into account such quadrupole correlations. This may solve the current puzzle of the stiffness of 208 Pb. It should be noted that it would also be relevant to measure the GMR on the Pb isotopic chain in order to provide a similar analysis than the one on the Sn nuclei.

In conclusion, it is shown that superfluidity favours the compressibility of nuclei, using a fully microscopic CHFB approach on the Tin isotopic chain. This may be the first evidence of a sizable effect of superfludity on the compressibility of a Fermionic system. Pairing effects should be described using a full microscopic HFB treatment. Doubly magic nuclei exhibit a specific increase of the GMR energy, due to the collapse of pairing. 208 Pb is therefore the "anomalous" data compared to the others. It is not possible to disentangle pairing interaction from the equation of state when providing the nuclear matter value of K ∞ . Indeed the pairing gap dependence on the nuclear matter incompressibility should be investigated, since it is shown that incompressibility decreases with increasing pairing gap in nuclei. Additional theoretical investigations are called for in order to predict the GMR including the mutually enhancement magicity effect. The macroscopic extraction of K sym may be ill-defined and should be extended to include pairing effects. Experimentally, measurements of the GMR in unstable nuclei should be performed in doubly magic 132 Sn, as well as extending the measurement on the Sn and Pb isotopic chains.

The author thanks G. Colò, M. Grasso, J. Margueron, Nguyen Van Giai, P. Ring, H. Sagawa and M. Urban for fruitful discussions about the results of this work.

FIG. 2 :

 2 FIG.2: Nuclear incompressibilities KA in 112-132 Sn isotopes calculated with the CHFB method and the SkM* interaction, as a function of the pairing gap ∆ predicted by the HFB calculation.


A Njifenjou

Finite volume methods for one-dimensional flow problems with discontinuous coefficients: A new approach based up on second order convergence

Keywords: Diffusion-reaction problems, New Finite Volume Scheme, Diffusion coefficient discontinuities, Approximate solutions, Second order convergence, Discrete energy norm, Stability, Error estimates

The finite volume methods are well known as powerful tools to address system of conservation equations. Over the past two decades a lot of effort has been made to put in place a mathematical framework for the theoretical analysis of Finite Volume Methods for second order differential operators. A perfect illustration of this progress is the development of Discrete Duality Finite Volume Methods (DDFV, for short). Following that trend we expose in this work a new approach of Finite Volume Schemes for second order elliptic problems in one dimension space, involving discontinuities in the diffusion coefficients. That new approach is based up on a system of two grids: a primary grid for having a control on localization of discontinuity points of the diffusion coefficient and a control-volume grid (which is not the dual of the primary grid) defined in such a way to get a second order approximation of the fluxes. The algebraic structure of the discrete problem we have got shows that this new approach is not a 1-D version of the DDFV even that both of them are cell-centered and vertex-centered. The proposed scheme leads to a symmetric, positive definite algebraic system and its solution meets the maximum principle. We have shown the second order convergence of the proposed scheme for pure diffusion problems on any primary grid (respecting quasi-uniformity of grid blocs). The second order convergence still holds for diffusion-reaction problems if the primary mesh elements are uniform, even if the diffusion coefficient gets discontinuity points. The ongoing work on the extension of the proposed method to two-dimensional problems is promising in terms of avoiding the computation of the equivalent diffusion coefficient to allocate to the diamond mesh elements as required from application of DDFV.

The Model problem and the data main assumptions

Let Ω =]a; b[ be an interval from the set of real numbers R, where a, b ∈ R with a < b are given. Also given as data are real-valued functions λ, f and µ.

We are interested in investigating a finite volume approximation of the solution u to the following diffusion-reaction problem:

Find u in an adequate function space such that:

-λ(x)u (x) + µ(x)u(x) = f (x) in Ω (1.1)

u(a) = u(b) = 0. (1.2)
We will proceed in such a way that the discrete problem involves discrete unknowns located at cell points and at vertex points of a primary coarse mesh. These points being cell centers of a finer mesh associated with the primary mesh. For obtaining the discrete equations, the numerical technique used is not that of the Discrete Duality Finite Volumes as we will see later.

Main assumptions on the data:

• On one hand the function λ is supposed to meet the following assumptions:

     • λ(x) = S
s=1 λ s 1 Os (x) a.e. in Ω

• ∃λ -, λ + ∈ R * + such that λ -≤ λ(x) ≤ λ + a.e. in Ω, (1.3) where R * + stands for the subset of R made up of nonnegative real numbers and where {O s } is a partition of Ω made up of nonempty open subintervals of Ω, while {1 Os } is a family of functions defined almost everywhere (a.e. for short) in Ω as:

1 Os (x) =    1 if x ∈ Int(O s ) = O s 0 if x ∈ Ext(O s ) (1.4)
where Int( ) and Ext( ) are respectively the interior and the exterior of a subinterval of R.

• On the other hand the functions f and µ are supposed to satisfy the following assumptions:

     • f ∈ L 2 (Ω)
• µ ∈ L ∞ (Ω), with µ(x) ≥ 0 a.e. in Ω. (1.5) Based on the previous assumptions, the Lax-Milgram theorem ensures existence and uniqueness of a weak solution to the problem (1.1)-(1.2), that is,    There exists a unique function u ∈ H 1 0 (Ω) such that:

Ω λu v dx + Ω µuvdx = Ω f vdx ∀ v ∈ H 1 0 (Ω) (1.6)
where the derivatives u and v are understood in the sense of distributions and where H 1 0 (Ω) is a well-known Sobolev space defined as

H 1 0 (Ω) = {v ∈ L 2 (Ω); v ∈ L 2 (Ω), v(a) = v(b) = 0}. (1.7)
The space H 1 0 (Ω) is endowed with its standard norm defined as

v H 1 0 = v L 2 (Ω) ∀ v ∈ H 1 0 (Ω). (1.8)
Since Ω is a bounded interval the previous norm is equivalent to the following one (direct consequence of Poincaré-Friedrichs inequality that we recall at the end of the current subsection):

v = [ v 2 L 2 (Ω) + v 2 L 2 (Ω) ] 1 2 ∀ v ∈ H 1 0 (Ω). (1.9)
This last norm is in fact the standard norm of the Sobolev space H 1 (Ω) defined as follows:

H 1 (Ω) = {v ∈ L 2 (Ω); v ∈ L 2 (Ω)}. (1.10)
It is then obvious that H 1 0 (Ω) is a subspace of the space H 1 (Ω). From now The standard norm of H 1 (Ω) is denoted by . H 1 (Ω) or simply . H 1 if there is no risk of confusion.

Note that the previous weak solution to the system (1.1)-(1.2) satisfies the following stability inequality:

u H 1 0 (Ω) ≤ C f L 2 (Ω) (1.11)
where C is a nonnegative real number. Let us end this subsection with the following important result the proof of which is a delicious exercise (see for instance [START_REF] Brézis | Functional Analysis, Sobolev spaces and Partial Differential Equations[END_REF]).

∃C > 0 such that ∀ v ∈ H 1 0 (Ω) v L 2 (D) ≤ C gradv L 2 (Ω) .
(1.12)

2 System of meshes for the new Finite Volume Scheme and related discrete function spaces

The new finite volume method we are developing in this work is based on the concept of primary "coarse" mesh associated with a control-volume mesh (to be extensively defined later).

2.1 The system of meshes required for the new finite volume scheme

The new finite volume method we are going to expose requires a system made up of two classes of meshes described in what follows.

• Primary coarse mesh

The primary mesh is the first mesh we define over Ω and any mesh refinement initiative is operated exclusively on that mesh, not anywhere else. Its main role is to give a precise delimitation of the different homogeneous subdomains of Ω in the context of piecewise constant diffusion coefficient. This assumption on the diffusion coefficient is very realistic for many complex engineering problems as subsurface multi-phase flows (see for instance [START_REF] Marle | Cours de production: Les écoulements polyphasiques en milieu poreux, Tome IV[END_REF][2] [START_REF] Galindez | Numerical simulation of 1-D oil and water displacements in petroleum reservoirs using the correction procedure via reconstruction (CPR) method[END_REF]).

Let N ∈ N be given, where N is the set of positive integer, and let x i+ 1 2 N i=0 be an increasing sequence made up of points from Ω = [a; b], i.e.

a = x 1 2 < x 3 2 < • • • < x i+ 1 2 < • • • < x N -1 2 < x N + 1 2 = b. (2.1)
Let us set what follows:

• Ω i def = x i-1 2 ; x i+ 1 2 ∀ i = 1, • • • , N (2.2)
• x i def = x i-1 2 + x i+ 1 2 2 ∀ i = 1, • • • , N (2.3)
• h i def = mes(Ω i ) ∀ i = 0, • • • , N + 1 (2.4)
where mes(.) stands for the Lebesgue measure in space of one dimension, with the following conventions:

h 0 = h N +1 = 0 (2.5)
and

x 0 = x 1 2 = a and x N +1 = x N + 1 2 = b. (2.6) Definition 2.

(Primary mesh)

The family {Ω i } N i=1 ; {x i } N i=1 defines a primary mesh over Ω. The mesh elements Ω i , i = 1, • • • , N , are called primary mesh elements in what follows.

We introduce a parameter h defined by h def = max h i 1≤i≤N with the assumption that there exists a mesh independent nonnegative constant ρ such that

h h i ≤ ρ ∀1 ≤ i ≤ N.
(2.7)

Figure 1: Primary mesh elements in black color, control-volumes in blue and red colors.

As usual this parameter is called the mesh size in what follows and is intended to tend to zero.

Important Assumption: The discontinuity points of the diffusion coefficient λ are located at the boundaries of mesh elements Ω i ] N i=1 .

• Computation finer mesh For computational purposes the primary mesh is associated with a finer mesh made up of two families of control-volumes as indicated in Figure 1 and in the following definition.

Definition 2.2 (Control-volumes)

The first family of control-volumes is defined as follows:

K i = x i - h i 4 ; x i + h i 4 i = 1, • • • , N (2.8)
with cell centers x i , for all i = 1, ..., N.

The second family of control-volumes is defined as follows:

K i+ 1 2 = x i+ 1 2 - h i 4 ; x i+ 1 2 + h i+1 4 i = 0, • • • , N (2.9)
with cell "centers" x i+ 1

2

, for all i = 0, ..., N and accounting with the conventions (2.5).

Note that the points x i+ 1 2 , for i ∈ {0, ..., N }, are not necessarily the midpoints of intervals K i+ 1 2 for i ∈ {0, ..., N } (see also Figure 1).

Definition 2.3 (Computation finer meshes)

The families

{K i } N i=1 , {x i } N i=1 and K i+ 1 2 N i=0 , x i+ 1 2 N i=0
define a computation finer mesh over Ω, simply denoted by M .

The points {x i } N i=1 and

x i+ 1 2 N i=0
are the place where are located the discrete

unknowns {u(x i ) ≡ u i } N i=1 and {u(x i+ 1 2 ) ≡ u i+ 1 2 } N -1 i=1 .
The new method of finite volumes we expose below indicate a way to get second order convergence for diffusion operators involving discontinuous coefficients. The finite volume solution to (1.1)-(1.2) will be denoted by

({u i } N i=1 , {u i+ 1 2 } N -1 i=1 )

Main theoretical tools

We present here some important theoretical tools for numerical analysis of the system (1.1)-(1.2). Among them are discrete function spaces, discrete gradient, adequate inner products and associated discrete norms, a discrete version of Poincaré-Friedrichs inequality, projection and interpolation operators. Such tools were introduced around 2005 by Omnes's team [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF] for numerical analysis of the Discrete Duality Finite Volume method discovered earlier (around 2000) by Hermeline [START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF] and Njifenjou -Moukouop [START_REF] Njifenjou | Traitement des anisotropies de perméabilité en simulation d'écoulement en milieu poreux par les volumes finis. Proceedings of an International Conference on "Systèmes Informatiques pour la Gestion de l'Environnement[END_REF].

Fundamental discrete function spaces and their properties

Let us start with introducing a basic function space denoted by S M,0 (called "discrete function space" because of that it depends on the computation mesh M) and defined by:

S M,0 def = {v h / v h (x) = N i=1 v i 1 K i (x) + N j=0 v j+ 1 2 1 K j+ 1 2 (x) a.e. in Ω, with v i , v j+ 1 2 ∈ R ∀1 ≤ i ≤ N ∀0 ≤ j ≤ N } (2.10)
where we have set, for any subinterval T of Ω and for almost every x ∈ Ω:

1 T (x) =    1 if x ∈ Int(T ) 0 if x ∈ Ext(T ) (2.11)
where Int( ) and Ext( ) respectively stand for interior and exterior of with respect to the standard topology of R. Let us mention the following obvious result.

Proposition 2.4 The family of functions ({1

K i } N i=1 , {1 K j+ 1 2 } N j=0
) is a (natural) basis of the discrete function space S M,0 . So the vector space S M,0 is of dimension 2N + 1. Let us introduce now an important subspace of S M,0 denoted by S M,0 0 and defined as :

S M,0 0 def = {v h ∈ S M,0 / v 1 2 = v N + 1 2 = 0}.
(2.12)

This subspace, with dimension 2N -1, involves the homogeneous Dirichlet boundary conditions (1.2). As we will see later it plays a key role in the finite volume reconstruction of the exact solution to (1.1)-(1.2), with second order convergence rate.

S M,0 α,β def = {v h ∈ S M,0 / v h (x 1 2 ) = α and v h (x N + 1 2 ) = β}.
(2.13)

We will come back to nonhomogeneous Dirichlet boundary conditions for further comments.

On the other hand let us set:

S M,0 0 def = {v h / v h (x) = N i=1 v i 1 K i (x) a.e. inΩ, with v i ∈ R ∀1 ≤ i ≤ N } (2.14)
and • Moreover we have the following relations:

S M,0 0 def = {v h / v h (x) = N -1 i=1 v i+ 1 2 1 K i+ 1 2 (x) a.e. in Ω, with v i+ 1 2 ∈ R ∀1 ≤ i ≤ N -1}.
S M,0 0 ≡ N i=1 P 0 (K i ) and S M,0 0 ≡ N -1 i=1 P 0 (K i+ 1 2
).

Discrete gradient operator

The definition of discrete gradient operator requires the introduction of the mesh D made up of elements from one of the forms:

D i+ 1 4 = [x i , x i+ 1 2 ] or D i-1 4 = [x i-1 2 , x i ], for i ∈ {1, 2, ..., N }. So we have D = {D 1-1 4 , D 1+ 1 4 , D 2-1 4 , D 2+ 1 4 , ..., D i-1 4 , D i+ 1 4 , ..., D N -1 4 , D N + 1 4 }.
Another ingredient is the discrete function space S D,0 defined as

S D,0 = {ζ h / ∃{(ζ i-1 4 , ζ i+ 1 4 )} N i=1 ⊂ R 2 such that ζ h (x) = (2.16) = N i=1 ζ i-1 4
1 D i-1 4 (x) + ζ i+ 1 4 1 D i+ 1 4 (x)
a.e. in Ω} (2.17)

We equip the space S D,0 with the following scalar product:

Definition 2.7 (Scalar Product)
We define over S D,0 a scalar product (., .) L 2 (Ω),D as follows :

(ζ h , ξ h ) L 2 (Ω),D = N i=1 h i 2 ζ i+ 1 4 ξ i+ 1 4 + ζ i-1 4 ξ i-1 4 ∀ ζ h , ξ h ∈ S D,0 .
The previous scalar product is associated with the following norm :

ζ h L 2 (Ω),D = N i=1 h i 2 ζ 2 i+ 1 4 + ζ 2 i-1 4 1 2 ∀ ζ h ∈ S D,0 . (2.18) Definition 2.

(Discrete Gradient operator)

A linear operator ∇ D from S M,0 to S D,0 is called a discrete gradient operator if for all v h ∈ S M,0 the following identity holds :

∇ D v h def = N i=1 2 h i (v i -v i-1 2 )1 D i-1 4 + (v i+ 1 2 -v i )1 D i+ 1 4 . (2.19)
Setting:

[∇ D v h ] i+ 1 4 = v i+ 1 2 -v i h i /2 and [∇ D v h ] i-1 4 = v i -v i-1 2 h i /2 (2.20)
the following obvious result holds:

Proposition 2.9 The mapping

v h -→ ∇ D v h L 2 (Ω),D def = N i=1 h i 2 ([∇ D v h ] i+ 1 4 ) 2 + ([∇ D v h ] i-1 4 ) 2 1 2 ∀ v h ∈ S M,0
is a semi-norm over the discrete function space S M,0 . remark 2.10 It is to notice that the mesh D could be defined as the family {D z j } 2N j=1 , where we have set for j ∈ {1, 2, 3, ..., 2N }:

z j = 3 4 + 1 2 (j -1) and D z j = [x z j -1 4 , x z j + 1 4
].

It then follows from what precedes that the discrete gradient ∇ D of a discrete function v h from the space S M,0 could be also defined as:

∇ D v h = 2N j=1 1 x z j + 1 4 -x z j -1 4 v z j + 1 4 -v z j -1 4 1 Dz j .

2.2.3

Inner products defined on S M,0 and S M,0 0 Recall that S M,0 0 is a subspace of S M,0 . Let us define over the space S M,0 the following inner products (where mes[.] is the Lebesgue measure in one-dimensional space ):

• First inner product on S M,0 and particular case of S M,0 0 :

(v h , w h ) L 2 (Ω),M def = N i=1 mes[K i ]v i w i + N i=0 mes[K i+ 1 2 ]v i+ 1 2 w i+ 1 2 ∀ v h , w h ∈ S M,0
(2.21) with associated norm . L 2 (Ω),M naturally defined by:

v h L 2 (Ω),M def = N i=1 mes[K i ]v 2 i + N i=0 mes[K i+ 1 2 ]v 2 i+ 1 2 1 2 ∀ v h ∈ S M,0 . (2.22)
The restriction of the norm . L 2 (Ω),M to the space S M,0 0 reads as follows:

v h * L 2 (Ω),M def = N i=1 mes[K i ]v 2 i + N -1 i=1 mes[K i+ 1 2 ]v 2 i+ 1 2 1 2 ∀ v h ∈ S M,0 0 . (2.23) since v 1 2 = v N + 1 2
= 0 for all v h ∈ S M,0 0 (see relation (2.12)).

• Second inner product on S M,0 :

(v h , w h ) H 1 (Ω),M,D def = (v h , w h ) L 2 (Ω),M + (∇ D v h , ∇ D w h ) L 2 (Ω),D ∀ v h , w h ∈ S M,0 .
(2.24) . The discrete norm associated with this inner product is denoted by . H 1 ,M,D and naturally defined by

v h 2 H 1 ,M,D = (v h , v h ) L 2 (Ω),M + (∇ D v h , ∇ D v h ) L 2 (Ω),D ∀ v h ∈ S M,0 . (2.25)
• Discrete version of Poincaré-Friedrichs inequality :

We recall that the "continuous" version of Poincaré-Friedrichs inequality is given by (1.12). Let us start with introducing some projection operators at least useful for the proof of the discrete version of Poincaré-Friedrichs inequality.

P M : S M,0 0 -→ S M,0 0 such that P M (v h ) = N i=1 v i 1 K i (2.26) P M : S M,0 0 -→ S M,0 0 such that P M (v h ) = N -1 i=1 v i+ 1 2 1 K i+ 1 2 .
(2.27) Lemma 2.12 The projection operators P M and P M satisfy the following continuity properties:

∀ v h ∈ S M,0 0      P M (v h ) L 2 (Ω),M ≤ mes[Ω] ∇ D v h L 2 (Ω),D P M (v h ) L 2 (Ω),M ≤ mes[Ω] ∇ D v h L 2 (Ω),D (2.28)
Proof.

• Let us prove the first inequality given by the previous Lemma i.e.

∀ v h ∈ S M,0 0 P M (v h ) L 2 (Ω),M ≤ mes[Ω] ∇ D v h L 2 (Ω),D .
Let us choose arbitrarily i ∈ {1, 2, ..., N } and v h ∈ S M,0 0 . So for a.e. x ∈ K i we have what follows:

v h (x) ≡ v i = [v 1 -v 1-1 2 ] + [-v 1 + v 1+ 1 2 ] + [v 2 -v 2-1 2 ] + [-v 2 + v 2+ 1 2 ] + ... + [v i -v i-1 2 ] (2.29)
We deduce that for a.e. x ∈ K i the following holds:

| v h (x) |=| v i |≤| v 1 -v 1-1 2 | + | -v 1 + v 1+ 1 2 | + | v 2 -v 2-1 2 | + | -v 2 + v 2+ 1 2 | +...+ + | v i -v i-1 2 | ≤ N j=1 | v j -v j-1 2 | + | v j -v j+ 1 2 | ≤ N j=1 h j /2 | v j -v j-1 2 | + | v j -v j+ 1 2 | 1 h j /2 (2.30)
It follows from Cauchy-Schawrz inequality that for a.e. x ∈ K i we have:

| v h (x) | 2 =| v i | 2 ≤ mes[Ω] N j=1 2 h j v j -v j-1 2 2 + v j -v j+ 1 2 2 = = mes[Ω] ∇ D v h 2 L 2 (Ω),D (2.31)
Integrating the both sides of the previous inequality in K i and summing on i ∈ {1, 2, ..., N } leads straightly to the investigated inequality.

• Let us prove the second inequality given by the previous Lemma i.e.

∀ v h ∈ S M,0 0 P M (v h ) L 2 (Ω),M ≤ mes[Ω] ∇ D v h L 2 (Ω),D .
As for the first inequality, let us choose arbitrarily i ∈ {1, 2, ..., N } and v h ∈ S M,0 0 . So for a.e. x ∈ K i+ 1 2

we have what follows:

v h (x) = v i+ 1 2 = (v i+ 1 2 -v i )+(v i -v i-1 2 )+(v i-1 2 -v i-1 )+(v i-1 -v i-3 2 )+(v i-3 2 -v i-2 )+...+ +(v 1 -v 1 2 ).
It follows that for a.e. x ∈ K i+ 1 2 we have

| v h (x) |=| v i+ 1 2 |≤ N j=1 h j /2 | v j -v j+ 1 2 | + | v j -v j-1 2 | 1 h j /2
Thanks to Cauchy-Schwarz we get for a.e.

x ∈ K i+ 1 2 | v h (x) | 2 =| v i+ 1 2 | 2 ≤ mes[Ω] ∇ D v h 2 L 2 (Ω),D
Integrating the both sides of the previous inequality in K i+ 1 2 and summing on i ∈ {1, 2, ..., N } lead to the second inequality of the previous Lemma.

We have the following discrete version of Poincaré-Friedrichs inequality.

Proposition 2.13 (Discrete Poincaré-Friedrichs inequality )

There exists a mesh independent nonnegative real number γ such that: The following mapping defined over S M,0 by

v h L 2 (Ω),M ≤ γ ∇ D v h L 2 (Ω),D ∀ v h ∈ S M,0 0 . ( 2
v h -→ ∇ D v h L 2 (Ω),D (2.33)
is a norm for S M,0 0 . Moreover this norm is equivalent on S M,0 0 to the standard norm of S M,0 introduced above and associated with the scalar product (., .) H 1 (Ω),M,D .

Terminology: The norm defined on S M,0 0 by (2.33) is called in the sequel "Discrete H 1 0 -norm" and considered as the standard norm of this discrete functional space. Note that the following bilinear form:

(v h , w h ) H 1 0 (Ω),M,D = (∇ D v h , ∇ D w h ) L 2 (Ω),D ∀ v h , w h ∈ S M,0 0 (2.34)
defines a scalar product on S M,0 0 , called in the sequel "standard scalar product" of S M,0 0 . This scalar product is associated with the "Discrete H 1 0 -norm". Let us give now a result that plays a key role in the proof of the stability of the new finite volume scheme that we will be exposing in the next section.

Proposition 2.15 (A fundamental result) Let Π M be a linear mapping defined over L 2 (Ω) with values in the space S M,0 ≡ S M,0 0 ⊕ S M,0 as follows :

Π M (v) = N i=1 v i 1 K i + N i=0 v i+ 1 2 1 K i+ 1 2 ∀ v ∈ L 2 (Ω) (2.35)
where S M,0 def = S M,0

0 ⊕ V ect(1 K 1 2 , 1 K N + 1 2
) and where we have set:

v i = 1 mes(K i ) K i v(x)dx ∀ 1 ≤ i ≤ N (2.36) and v i+ 1 2 = 1 mes(K i+ 1 2 ) K i+ 1 2 v(x)dx ∀ 0 ≤ i ≤ N. (2.37)
Then Π M is continuous with respect to the norms . L 2 (Ω) and . L 2 (Ω),M respectively defined in L 2 (Ω) and S M,0 , i.e. there exists , a mesh-independent nonnegative real number, such that

Π M (v) L 2 (Ω),M ≤ v L 2 (Ω) ∀ v ∈ L 2 (Ω).
Proof. Let v be a function from the space L 2 (Ω). So according to the definition of the operator Π M (see relations (2.35)-(2.37)) we have :

Π M (v) 2 L 2 (Ω),M = N i=1 h i v 2 i + N i=0 h i+ 1 2 v 2 i+ 1 2 = (2.38) = N i=1 1 h i [ K i v(x)dx] 2 + N i=0 1 h i+ 1 2 [ K i+ 1 2 v(x)dx] 2 ≤ (2.39) ≤ N i=1 K i |v(x)| 2 dx + N i=0 K i+ 1 2 |v(x)| 2 dx = (2.40) = v 2 L 2 (Ω) (2.41)
According to what precedes we have established that

Π M (v) L 2 (Ω),M ≤ v L 2 (Ω) .
(2.42)

So ends the proof.

3 A new finite volume Scheme for the system (1.1)-(1.2)

We expose here the main steps leading to a new finite volume approximation of the balance equation (1.1), with prescribed homogeneous Dirichlet boundary conditions (1.2).

Notations and definitions

Let ϕ(.) be a (real-valued or vector-valued) function defined over Ω and x ∈ Ω. We denote by ϕ(x + ) and ϕ(x -) respectively the limits of ϕ(s) when s 2) and u is its gradient (in one dimension space). Let us define the flow velocity q and the corresponding flux F at the right-side and the left-side of every point x ∈ Ω as follows:

     q(x τ , u ) def = -[λ u ](x τ ) F (x τ , u ) def = q(x τ , u ) • ν(x τ ) (3.1)
where we have set

ν(x τ ) def =      -1 if τ = + +1 if τ = - (3.2)

Discrete mass conservation per control-volume

We proceed in several steps as it follows.

• Step one: Mass conservation principle for control-volumes

Integrating the balance equation (1.1) in the control-volumes K i and K i+ 1 2 , leads to

         F (x - i+ 1 4 , u ) + F (x + i-1 4 , u ) + K i µ(x)u(x)dx = h i 2 f i i = 1, • • • , N F (x - i+ 3 4 , u ) + F (x + i+ 1 4 , u ) + K i+ 1 2 µ(x)u(x)dx = h i +h i+1 4 f i+ 1 2 i = 1, • • • , N -1.
(3.3) where we have set

         f i def = 1 mes(K i ) K i f (x)dx f i+ 1 2 def = 1 mes(K i+ 1 2 ) K i+ 1 2 f (x)dx (3.4)
and where

x i±

1 4 def = x i ± h i 4 , i = 1, 2, ..., N. (3.5)
Remark that the term

x i+ 3 4 def = x i+ 1 2 + h i+1 4 (3.6)
appearing in the second equation of the system (3.3) is of the form (3.5). Indeed we have

x i+ 3 4 = x (i+1)-1 4 i = 0, 1, ..., N -1. (3.7)
It follows from definitions (3.1) that

∀ 1 ≤ i ≤ N        F (x - i+ 1 4 , u ) = -λ i u (x i + h i 4 ) F (x + i-1 4 , u ) = λ i u (x i -h i 4 ) (3.8)
and

∀ 1 ≤ i ≤ N -1        F (x - i+ 3 4 , u ) = -λ i+1 u (x i+ 1 2 + h i+1 4 ) F (x + i+ 1 4 , u ) = λ i u (x i+ 1 2 -h i 4 ).
(3.9)

• Step two: Definition of discrete flux function

Let us investigate now second order approximations of the following flux terms: [ such that

F (x - i+ 1 4 , u ), F (x + i-1 4 , u ), F (x -
u(x i+ 1 2 ) = u(x i + h i 4 ) + h i 4 u (x i + h i 4 ) + h 2 i 32 u (x i + h i 4 ) + ( h i 4 ) 3 u (z d i+ 1 4 ) (3.10) u(x i ) = u(x i + h i 4 ) - h i 4 u (x i + h i 4 ) + h 2 i 32 u (x i + h i 4 ) -( h i 4 ) 3 u (z g i+ 1 4
). (3.11) Subtracting (3.11) from (3.10) side by side leads to

u (x i + h i 4 ) = u(x i+ 1 2 ) -u(x i ) h i /2 + h 2 i 32 u (z i+ 1 4 ) ∀i = 1, • • • , N. (3.12)
where

x i < z g i+ 1 4 < z i+ 1 4 < z d i+ 1 4 < x i+ 1 2
. A similar development yields what follows:

u (x i - h i 4 ) = u(x i ) -u(x i- 1 2
)

h i /2 + h 2 i 32 u (r i-1 4 ), ∀i = 1, • • • , N. (3.13) with x i-1 2 < r i+ 1 4 < x i .
Therefore, thanks to (3.8), we can deduce that

∀i = 1, • • • , N      F (x - i+ 1 4 , u ) = F (x - i+ 1 4 , ∇ D u h ) -R - i+ 1 4 (h, u ) F (x + i-1 4 , u ) = F (x + i-1 4 , ∇ D u h ) -R + i-1 4 (h, u ) (3.14)
where, according to Definition 2.8, we have

∇ D u h def = N i=1 1 h i /2 (u i -u i-1 2 )1 D i-1 4 + (u i+ 1 2 -u i )1 D i+ 1 4
, and where we have set

∀i = 1, • • • , N      F (x - i+ 1 4 , ∇ D u h ) = -λ i [∇ D u h ] i+ 1 4 F (x + i-1 4 , ∇ D u h ) = λ i [∇ D u h ] i-1 4 , (3.15
)

with [∇ D u h ] i+ 1 4 and [∇ D u h ] i-1 4
defined as it follows:

Definition 3.2 ∀i = 1, • • • , N            [∇ D u h ] i+ 1 4 def = u(x i+ 1 2 ) -u(x i ) h i /2 [∇ D u h ] i-1 4 def = u(x i ) -u(x i-1 2 ) h i /2 . (3.16)
Using the same arguments as for the computation of the fluxes F (x -

i+ 1 4
, u ) and

F (x + i-1 4 , u ) (see (3.14) above) leads to ∀i = 1, • • • , N -1            u (x i+ 1 2 + h i+1 4 ) = u(x i+1 )-u(x i+ 1 2 ) h i+1 /2 + h 2 i+1 32 u ( z i+ 3 4 ) u (x i+ 1 2 -h i 4 ) = u(x i+ 1 2 )-u(x i ) h i /2 + h 2 i 32 u ( r i+ 1 4
) (3.17) Therefore, thanks to (3.9), we deduce that

∀i = 1, • • • , N -1      F (x - i+ 3 4 , u ) = F (x - i+ 3 4 , ∇ D u h ) -R - i+ 3 4 (h, u ) F (x + i+ 1 4 , u ) = F (x + i+ 1 4 , ∇ D u h ) -R + i+ 1 4 (h, u ) (3.18)
where F (x -

i+ 3 4 , ∇ D u h ) and F (x + i+ 1 4
, ∇ D u h ) are defined as it follows:

Definition 3.3 ∀i = 1, • • • , N -1        F (x - i+ 3 4 , ∇ D u h ) def = -λ i+1 [∇ D u h ] i+ 3 4 F (x + i+ 1 4 , ∇ D u h ) def = λ i [∇ D u h ] i+ 1 4 (3.19) with [∇ D u h ] i+ 3 4
defined as it follows:

Definition 3.4 ∀i = 1, • • • , N -1 [∇ D u h ] i+ 3 4 def = u(x i+1 ) -u(x i+ 1 2 ) h i+1 /2 (3.20)
Recall that {O s } S s=1 is a family of open (nonempty) intervals defining a partition of Ω and associated with the diffusion coefficient λ in the sense that:

λ(x) = S s=1 λ s 1 Os (x)
a.e. in Ω.

It follows from the previous development that:

Proposition 3.5 If the exact solution u to the problem (1.1)-(1.2) is such that the restriction of u to O s , for all 1 ≤ s ≤ S, lies in C 3 (O s ) then |R ± i+ 1 4 (h, u )| ≤ C h 2 ∀ i ∈ { 1 2 , 1, 3 2 , 2, ..., N - 1 2 , N } (3.21)
where C is a nonnegative mesh-independent number. Moreover the following holds:

R - i+ 1 4 (h, u ) + R + i+ 1 4 (h, u ) = 0 ∀ i ∈ { 1 2 , 1, 3 2 , 2, ..., N - 1 2 , N }. (3.22)
•

Step three: Approximation of reaction term integrals

The term µ(x)u(x) in the left-hand side of the balance equation (1.1) is named the reaction-term. We should look for (at least) a second order approximation of the reaction-term contribution in the integral formulation of the balance equation per control-volume (see system of equations (3.3)). So we should perform these integral Finite Volume for diffusion problems: New approach approximations for the control-volumes K i and K i+ 1 2 , for i ∈ {1, 2, ..., N }. Using the Rectangle-centered quadrature leads to

K i µ(x)u(x)dx = h i 2 [µ(x i )u(x i )] -E i (h, u ) ∀i ∈ {1, 2, ..., N } (3.23) with | E i (h, u ) | ≤ Ch 3 ∀i ∈ {1,
K i+ 1 2 µ(x)u(x)dx = h i + h i+1 4 [µ(x i+ 1 2 )u(x i + 1 2 )]-E i+ 1 2 (h, u ) ∀i ∈ {1, 2, ..., N -1}. (3.25) with | E i+ 1 2 (h, u ) | ≤ Ch 2 ∀i ∈ {1, 2, ..., N -1} (3.26)
as soon as the exact solution u to (1.1)-(1.2) lies in C 1 (Ω s ), for all s ∈ {1, 2, ..., S}.

•

Step four (and the Last one): Equations of the discrete mass conservation principle per control-volume Let us start with introducing the following simplified notations:

u i-1 2 ≡ u(x i-1 2 ), u i ≡ u(x i ), u i+ 1 2 ≡ u(x i+ 1 2 ), u i+1 ≡ u(x i+1 ), µ i ≡ µ(x i ), µ i+ 1 2 ≡ µ(x i+ 1 2
), ... We deduce from the system of equations (3.3) and from the steps two and three above that in the control-volumes K i and K i+ 1 2 we have the following discrete mass conservation principle :

∀ i = 1, • • • , N        λ i [u i -u i+ 1 2 ] h i /2 + λ i [u i -u i-1 2 ] h i /2 + h i 2 µ i u i = = h i 2 f i + R - i+ 1 4 (h, u ) + R + i-1 4 (h, u ) + E i (h, u ) (3.27)
and

∀ i = 1, • • • , N -1        λ i+1 [u i+ 1 2 -u i+1 ] h i+1 /2 + λ i [u i+ 1 2 -u i ] h i /2 + [h i +h i+1 ] 4 µ i+ 1 2 u i+ 1 2 = = [h i +h i+1 ] 4 f i+ 1 2 + R - i+ 3 4 (h, u ) + R + i+ 1 4 (h, u ) + E i+ 1 2 (h, u ) (3.28) where the following truncation errors R ± i± 1 4 (h, u ), R - i+ 3 4 (h, u ), E i (h, u ) and E i+ 1 2 (h, u )
are introduced in what precedes with their respective estimates as well.

Definition of a New Finite Volume Scheme

From the system of equations (3.27)-(3.28) someone can easily see that when the truncation errors are neglected the discrete unknowns {u i } N i=1 and {u i+ 1 2 } N -1 i=1 satisfy the following system of "equations":

λ i [u i -u i+ 1 2 ] h i /2 + λ i [u i -u i-1 2 ] h i /2 + h i 2 µ i u i ≈ h i 2 f i ∀ i = 1, • • • , N (3.29)
and

λ i+1 [u i+ 1 2 -u i+1 ] h i+1 /2 + λ i [u i+ 1 2 -u i ] h i /2 + [h i + h i+1 ] 4 µ i+ 1 2 u i+ 1 2 ≈ ≈ [h i + h i+1 ] 4 f i+ 1 2 ∀ i = 1, • • • , N -1. (3.30)
Inspired by the previous system, let us introduce the following discrete problem:

Find in the discrete space S M,0 0 a function

u h = N i=1 u i 1 K i + N -1 i=1 u i+ 1 2 1 K i+ 1 2 (3.31)
such that its components

({u i } N i=1 , {u j+ 1 2 } N -1 j=1
), in the natural basis of S M,0 0 , are characterized as a solution of the following system of equations:

λ i [u i -u i+ 1 2 ] h i /2 + λ i [u i -u i-1 2 ] h i /2 + h i 2 µ i u i = h i 2 f i ∀ i = 1, • • • , N (3.32)
and

λ i+1 [u i+ 1 2 -u i+1 ] h i+1 /2 + λ i [u i+ 1 2 -u i ] h i /2 + [h i + h i+1 ] 4 µ i+ 1 2 u i+ 1 2 = = [h i + h i+1 ] 4 f i+ 1 2 ∀ i = 1, • • • , N -1. ( 3

Theoretical analysis of the New Finite Volume Scheme

We intend to expose in the current section some important mathematical properties of the New Finite Volume Scheme. Let us start with proving existence, uniqueness and stability of u h in S M,0 0 .

Existence, Uniqueness and Stability for the solution to the system (3.31)-(3.33)

Let v h be a discrete function from the space S M,0 0 , chosen arbitrarily, with components

({v i } N i=1 , {v j+ 1 2 } N -1 j=1 }) in the basis ({1 K i } N i=1 , {1 K j+ 1 2 } N -1 j=1 }).
Multiplying the two sides of equations (3.32) and (3.33) with respectively v i and v j+ 1 2 , summing on i ∈ {1, 2, ..., N } and on j ∈ {1, 2, ..., N -1}, and doing some re-organization of terms lead to the following variational problem:

Find u h ∈ S M,0 0 such that: N i=1 λ i h i /2 u i -u i-1 2 v i -v i-1 2 + u i -u i+ 1 2 v i -v i+ 1 2 + + N i=1 h i 2 µ i u i v i + h i + h i+1 4 µ i+ 1 2 u i+ 1 2 v i+ 1 2 = = N i=1 h i 2 f i v i + h i + h i+1 4 f i+ 1 2 v i+ 1 2 ∀v h ∈ S M,0 0 . (4.1)
In terms of discrete gradient the previous variational problem can be re-written as Find u h ∈ S M,0 0 such that:

N i=1 λ i h i 2 [∇ D u h ] i-1 4 [∇ D v h ] i-1 4 + [∇ D u h ] i+ 1 4 [∇ D v h ] i+ 1 4 + + N i=1 h i 2 µ i u i v i + h i + h i+1 4 µ i+ 1 2 u i+ 1 2 v i+ 1 2 = = N i=1 h i 2 f i v i + h i + h i+1 4 f i+ 1 2 v i+ 1 2 ∀v h ∈ S M,0 0 . (4.2)
First of all, remark that the following equivalence holds: Proof. Follow closely the arguments developed for a similar result in [START_REF] Njifenjou | Overview on conventional finite volumes for elliptic problems involving discontinuous diffusion coefficients . Part I: Focus on the one dimension space models[END_REF], pages 23-24.

u h H 1 0 (Ω),M,D ≤ C f L 2 (Ω) . (4.3)
where C is a mesh independent nonnegative real number.

Proof. Recall that the space S M,0 0 is a closed subspace of the Hilbert space S M,0 which is equipped with its standard inner product (., .) H 1 (Ω),M,D . So S M,0 0 is also a Hilbert space with respect to this inner product. That being said, let us set:

B(w h , v h ) def = N i=1 λ i h i 2 [∇ D w h ] i-1 4 [∇ D v h ] i-1 4 + [∇ D w h ] i+ 1 4 [∇ D v h ] i+ 1 4 + + N i=1 h i 2 µ i w i v i + h i + h i+1 4 µ i+ 1 2 w i+ 1 2 v i+ 1 2 ∀w h , v h ∈ S M,0 0 . (4.4)
and

L(v h ) = N i=1 h i 2 f i v i + h i + h i+1 4 f i+ 1 2 v i+ 1 2 ∀v h ∈ S M,0 0 . (4.5)
• Let us prove Existence and Uniqueness of the solution to the variational equation (4.2) by an application of the Lax-Milgram theorem. It is obviously seen that B(., .) and L(.) are respectively a bilinear form and a linear form over S M,0 0 . Let us check if the conditions of Lax-Milgram are satisfied by these two forms (see for instance [START_REF] Brézis | Functional Analysis, Sobolev spaces and Partial Differential Equations[END_REF] to learn more about theoretical aspects of Lax-Milgram theorem, notably the proof of this theorem; see for instance [START_REF] Njifenjou | Overview on conventional finite volumes for elliptic problems involving discontinuous diffusion coefficients . Part I: Focus on the one dimension space models[END_REF] [START_REF] Njifenjou | Introduction to Finite Element Methods[END_REF] for application of Lax-Milgram theorem to Numerical Analysis of discrete models for diffusion problems).

• Continuity of B(., .): From the fact that (A + B) 2 ≤ 2[A 2 + B 2 ] for all real numbers A and B, we easily get that:

| B(w h , v h ) | 2 ≤ ≤ 2 N i=1 λ i h i 2 | [∇ D w h ] i-1 4 | + | [∇ D w h ] i+ 1 4 | | [∇ D v h ] i-1 4 | + | [∇ D v h ] i+ 1 4 | 2 + + 2 Ω µ(x) | w h (x) || v h (x) | dx 2 ∀w h , v h ∈ S M,0 0 . (4.6)
Thanks to assumptions in (1.3) and (1.5) we can get what follows from the previous inequality.

| B(w h , v h ) | 2 ≤ 2λ 2 + N i=1 h i 2 | [∇ D w h ] i-1 4 | + | [∇ D w h ] i+ 1 4 | | [∇ D v h ] i-1 4 | + | [∇ D v h ] i+ 1 4 | 2 + + 2 µ 2 L ∞ (Ω) Ω | w h || v h | dx 2 ∀w h , v h ∈ S M,0 0 . (4.7)
From a double application of Cauchy-Schwartz inequality in the right-hand side of the previous inequality and thanks to Proposition 2.13 (discrete version of Poincaré-Friedrichs) we obtain the continuity of the bilinear form B(., .).

• Coercivity of B(., .): Let w h be an arbitrarily chosen function from the space S M,0 0 . So we have

B(w h , w h ) = = N i=1 λ i h i 2 ([∇ D w h ] i-1 4 ) 2 + ([∇ D w h ] i+ 1 4 ) 2 + + N i=1 h i µ i 2 w 2 i + h i + h i+1 4 µ i+ 1 2 w 2 i+ 1 2 ) ≥ ≥ λ - N i=1 h i 2 ([∇ D w h ] i-1 4 ) 2 + ([∇ D w h ] i+ 1 4 ) 2 = λ -∇ D w h 2 L(Ω),D . (4.8)
where λ -is a mesh independent nonnegative number coming from the assumption (1.3).

• Continuity of the linear form L(.): Let w h be an arbitrarily chosen function from the space S M,0 0 . So we have to show that there exists a mesh independent nonnegative real number β such that

| L(w h ) | ≤ β ∇ D w h L 2 (Ω),D ∀w h ∈ S M,0 0 . (4.9)
There exist (at least) two ways to prove the inequality (4.9). * First way: Remark that

L(w h ) = N i=1 K i f (x)w h (x)dx + N i=1 K i+ 1 2 f (x)w h (x)dx = (4.10) = Ω f (x)w h (x)dx (4.11)
Since the norm . L 2 (Ω),M is the restriction to S M,0 of the standard norm of the well-known Hilbert space L 2 (Ω), the Cauchy-Schawrz inequality applies and we get

| L(w h ) | ≤ f L 2 (Ω) w h L 2 (Ω),M ∀w h ∈ S M,0 0 . (4.
12)

It follows from the discrete version of Poincaré-Friedrichs inequality that there exists a mesh independent nonnegative real number C such that

| L(w h ) | ≤ C ∇ D w h L 2 (Ω),D ∀w h ∈ S M,0 0 . (4.13)
The continuity of L(.) is proven. * Second way: In virtue of the identity (A + B) 2 ≤ 2A 2 + 2B 2 that holds for all A, B ∈ R, we have what follows for all w ∈ S M,0 0 :

| L(w h ) | 2 ≤ 2 N i=1 w i K i f (x)dx 2 + 2   N i=1 w i+ 1 2 K i+ 1 2 f (x)dx   2 ≤ (from discrete Cauchy-Schwarz inequality we get) ≤ 2 N i=1 h i 2 w 2 i N i=1 1 h i /2 K i f (x)dx 2 + + 2 N i=1 h i + h i+1 4 w 2 i+ 1 2   N i=1 4 h i + h i+1   K i+ 1 2 f (x)dx   2   ≤
(by Cauchy-Schwarz applied to the L 2 -scalar product we have)

≤ 2 N i=1 h i 2 w 2 i N i=1 K i [f (x)] 2 dx + + 2 N i=1 h i + h i+1 4 w 2 i+ 1 2   N i=1 K i+ 1 2 [f (x)] 2 dx   Therefore | L(w h ) | 2 ≤ 2 f 2 L 2 (Ω) w h 2 L 2 (Ω),M . (4.14)
Thanks to the discrete version of Poincaré-Friedrichs we could conclude that L(.) is a continuous linear form.

• The Stability of the discrete solution u h straightly follows from the coercivity of the bilinear form B(., .) and the continuity of the linear form L(.).

Matrix properties of the Scheme (3.32)-(3.33)

The matrix form of the new Finite Volume Scheme (3.32)-(3.33) may be expressed as follows:

A h B h (B h ) t C h U h cc U h vc = F h cc F h vc (4.15)
We can now define the error function e h as follows

e h (x) = N i=1 e i 1 K i (x) + N -1 i=1 e i+ 1 2 1 K i+ 1 2 (x)
a.e. in Ω.

It is then clear that e h lies in S M,0 0 . It is easily seen that this discrete function is solution of a system of the same type as (3.27)-(3.28). Indeed, subtracting side by side equation (3.32) from equation (3.27) and equation (3.33) from (3.28) lead to the so-called Error function system that reads as Concluding remark: Note that P. Omnes in [18] has obtained a second-order convergence for a function reconstructed from Finite Volume approximation of the Laplace operator. The previous theorem asserts that one can get a second order convergence for the finite volume solution of any 1-D diffusion problem involving discontinuous coefficients (under reasonable assumptions for the rest of data).

             2λ i h i /2 e i -e i+ 1 2 + λ i h i /2 e i -e i-1 2
+ h i 2 µ i e i = R - i+ 1 4 + R + i-1 4 + E i ∀i = 1, • • • , N λ i+1 h i+1 /2 e i+ 1 2 -e i+1 + λ i h i /2 e i+ 1 2 -e i + h i +h i+1 4 µ i+ 1 2 e i+ 1 2 = = R - i+ 3 4 + R + i+ 1 4 + E i+ 1 2 ∀i = 1, • • • , N -1.
B(e h , e h ) = N i=1 R - i+ 1 4 (e i -e i+ 1 2 ) + N i=1 R + i-1 4 (e i -e i-1 2 ) + N i=1 E i e i + N -1 i=1 E i+ 1 2 e i+ 1 2 . ( 4

Proposition 1 . 1 (

 11 Poincaré-Friedrichs) Let d be a given space dimension and D a nonempty open subset of R d . If D is bounded in at least one direction then :

remark 2 . 5

 25 Note that when one replaces the boundary conditions (1.2) with the following one: u(a) = α and u(b) = β the discrete function framework S M,0 0 should be replaced with what follows:

( 2 . 2 = v N + 1 2 = 0 .

 2220 [START_REF] Njifenjou | A discrete duality finite volume method for flow problems with prescribed periodic boundary conditions[END_REF] accounting with the boundary conditions v 1 Note that the function spaces S M,0 0 and S M,0 0 are obviously subspaces of S M,0 0 . Moreover the following properties hold: Proposition 2.6 (Identification of spaces)• The following identities hold:

Definition 2 . 11 (

 211 Projection operators) Let P M and P M be two projection operators defined on the space S M,0 0 with values in S M,0 0 and S M,0 0 respectively (see relations (2.12) (2.14) and (2.15) for the definition of these spaces):

>-→

  x and when s < -→ x. In the same order of ideas if ζ(., .) is a (real-valued or vector-valued) function defined over Ω × R, we denote by ζ(x + , z) and ζ(x -, z) respectively the limits of ζ(s, z) when s > -→ x and when s < -→ x for all z ∈ R. Definition 3.1 (Flow Velocity and Flux) Recall that u is the exact solution to (1.1)-(1.

i+ 3 4 , 4 , 1 4 and z g i+ 1 4 4 , x i+ 1 2 [

 441142 u ) and F (x + i+ 1 u ). Under the assumption that the exact solution u to the continuous problem (1)-(2) lies in C 3 (Ω i ), for all i ∈ {1, 2, ..., N }, the Taylor-Lagrange theorem ensures that there exist z d i+ respectively in the intervals ]x i+ 1 and ]x i , x i+ 1 4

Proposition 4 . 1 (

 41 Equivalence between the two formulations) Any function from the discrete function space S M,0 0 is a solution of the system of equations (3.32)-(3.33) if and only if it is a solution of the variational equation (4.2).

Proposition 4 . 2 (

 42 Existence, Uniqueness and Stability Results) (i) There exists a unique function u h in the space S M,0 0 solving the variational equation (4.2). (ii) Moreover the solution u h to (4.2) satisfies the following inequality (named Stability inequality):

( 4 .

 4 23)Let us now investigate some estimates of the Error function with respect to the norm . H 1 0 (Ω),M also denoted by . L 2 (Ω),D and the standard norm of L 2 (Ω). For that purpose let us multiply the two sides of the first equation of the system (4.23) by e i and sum on i ∈ {1, 2, ..., N }. Let us repeat the same operations with the second equation of the same system, but with e i+ 1 2 instead of e i , and i ∈ {1, 2, ..., N -1} instead of i ∈ {1, 2, ..., N }. Adding side by side the two previous sums and re-ordering the terms, accounting with(3.22), lead to what follows:

4 (e i -e i+ 1 2 ) 2 ≤ 4 (e i -e i-1 2 ) 2 ≤ C h 4 N i=1 1 h i / 2 (e i -e i-1 2 ) 2 . 4 N i=1 E i e i 2 ≤ C h 4 N i=1 h i e 2 i = C h 4 P M (e h ) 2 L 2 ( 2 ≤ C h 2 N - 1 i=1 h i + h i+1 4 e 2 i+ 1 2 = 2 L 2 2 L 2 ( 2 L 2 (E i e i 2 ≤ C h 4 ∇ D (e h ) 2 L 2 ( 2 ≤ 2 L 2 4 (e i -e i-1 2 ) 2 ≤

 422422412224i=124242222122222222422222422 .24) From the discrete version of Cauchy-Schwarz inequality we get [B(e h , e h )] 2 ≤ 4 discrete version of Cauchy-Schwarz to the square of the first two sums from the right-hand side of the previous inequality yields, accounting with (3.21) from Proposition 3.5 (note that in what follows C represents diverse mesh independent nonnegative numbers): (4.27) Repeating the same exercise with the square of the two last sums from the right-hand side of (4.25) and accounting with (2.7, (3.24) and (3.26), leads to C h 2 P M (e h ) (Ω),M (4.29)Applying the Lemma 2.12 to the termsC h 4 P M (e h ) Ω),M and C h 2 P M (e h ) C h 2 ∇ D (e h ) (Ω),D(4.31)On the other hand, it follows from inequalities (4.26)-(4.27) that 4 C h 4 ∇ D (e h ) 2 L 2 (Ω),D . (4.32) From the inequalities (4.25), (4.30), (4.31), (4.32) and thanks to the coercivity of the bilinear form B(., .) we can state what follows:

Theorem 4 . 4 (

 44 Error Estimates) Recall that the finite family of non empty subintervals {O s } s∈S (associated with the diffusion coefficient λ(.)) defines a partition of Ω. Consider the assumptions (1.3), (1.5), (2.7) and the following ones:(i) The discontinuity points of the diffusion coefficient λ(.) are part of the set {x i+ 1 2 } N i=0 associated with the primary (relatively coarse) mesh, (ii) The exact solution u is such that the restriction u | Os of u to O s honors the following condition:u | Os ∈ C 3 (O s ) ∀ s ∈ S.Then the finite volume approximation u h of the exact solution u to (1.1)-(1.2) is such that:•In the context of non uniform primary mesh elements Ω i ] N i=1 combined with non negligible reaction effects, the Error function e h = u h -u h satisfies the following estimates (first order convergence):∇ D e h L 2 (Ω),D ≤ C h and e h L 2 (Ω),M ≤ C h;• In the context of uniform primary mesh elements Ω i ] N i=1 combined with non negligible reaction effects, the Error function e h meets what follows (second order convergence):∇ D e h L 2 (Ω),D ≤ C h 2and e h L 2 (Ω),M ≤ C h 2 ;• In the context of pure diffusion problems (i.e. reaction effects negligible), over non-uniform primary meshes, the Error function e h satisfies the following estimates (second order convergence):∇ D e h L 2 (Ω),D ≤ C h 2and e h L 2 (Ω),M ≤ C h 2 .

  2, ..., N } (3.24) if the exact solution u to (1.1)-(1.2) lies in C 2 (Ω s ), for all s ∈ {1, 2, ..., S}.

	Since x i+ 1 2 centered quadrature does not apply. Nevertheless a simple Rectangle quadrature still the Rectangle-is not necessarily the center-point for the interval K i+ 1 2
	applies and leads to what follows:

  .33)

	Definition 3.6 (New Finite Volume Scheme)
	The system of relations (3.32)-(3.33) is what we call a New Finite Volume Scheme
	for one dimension space Diffusion Problems involving discontinuous diffusion coeffi-
	cients.

Note that the previous discrete scheme is not a 1-D version of Discrete Duality Finite Volumes (DDFV, for short). Concerning literature on DDFV see for instance

[START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF]

[10]

[START_REF] Njifenjou | Analysis on general meshes of a discrete duality finite volume method for subsurface flow problems[END_REF]

[12][13][14]

[START_REF] Njifenjou | A discrete duality finite volume method for flow problems with prescribed periodic boundary conditions[END_REF]

for a conventional formulation and

[START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF]

[16]

[START_REF] Delcourte | A Discrete Duality Finite Volume Approach to Hodge Decomposition and div-curl Problems on Almost Arbitrary[END_REF]

for formulations based on discrete gradient and discrete divergence. See also references therein.

where we have set :

and U h vc = u i+ 1 2 1≤i≤N -1 (4.16) and where: F h cc is a sub-vector with N components defined only by the right hand side of (3.32) and F h vc a sub-vector with (N -1) components defined only by the right hand side of (3.33) as the boundary conditions are of homogeneous Dirichlet ones. Concerning the sub-matrices A h , C h and B h , note that the first two ones are respectively N × N and (N -1)×(N -1) diagonal matrices while B h is a N ×(N -1) matrix with a maximum of two coefficients different from 0 per line, So the symmetry and the sparse structure of the matrix associated with the new finite volume scheme are established. Proof. Since the symmetry of M h is obvious let us concentrate on positive definiteness and monotonicity.

• Positive definiteness: It follows straightly from the coercivity of the bilinear form B(., .) introduced in the proof of Proposition 4.2.

• Monotonicity: Recall that M h is monotone if the components of any solution u h to the system (3.32)-(3.33) are positive as soon as the components of Π M f are positive (see Proposition 2.15 for the definition of the operator Π M ). Recall that M h is monotone if and only if M h is nonsingular and the coefficients of its inverse are positive. At least there are two ways to prove that M h is monotone.

• First way is a Classical technique: We start with setting

Let us suppose that 1 2 < r ≤ N. (4.17)

From the previous assumption there are two possibilities:

First possibility: r ∈ {1, 2, ..., i-1, i, i+1, ..., N -1, N }. In this case the following inequalities hold:

Second possibility: r ∈ { 3 2 , ..., i -1 2 , i + 1 2 , ..., N -1 2 }, with i < r < i + 1, where i is an integer from the set {1, 2, ..., N -1}. In this case the following inequalities hold:

This is also absurd. So r ∈ { 3 2 , ..., i -1 2 , i + 1 2 , ..., N -1 2 } is not a possibility. We conclude that the assumption (4.17) is wrong. In consequence r = 1 2 . This implies that all the components of u h are positive since u 1 2 = 0.

• Second way is a geometrical technique: Exposed for the first time, according to our knowledge, in [START_REF] Njifenjou | Geometric arguments for proving the discrete discrete maximum principle met by conventional finite volume schemes in the context of isotropic diffusion problems[END_REF] for two-dimensional diffusion problems, we are going right now to apply it to the new Finite Volume Scheme (3.32)-(3.33).

Let us suppose that the right-hand side of the system (3.32)-(3.33) is positive. We should deduce that all the components

of its solution u h are positive.

Let us denote by u σ the smallest of the quantities u j , with j ∈ { 1 2 , 1, 3 2 , 2, 5 2 , ..., N -

If σ = 1 or σ = N the proof is ended. Otherwise, consider the straight semi-line ∆ σ with origin the node x σ and passing through x 1 2

. Writing the discrete balance for the node x σ-1 2 and accounting with (4.21) we see that

= 1 the proof is ended, otherwise repeat the procedure until the equality [σ -1] = 1 holds. This should happen after a finite number of iterations. The case where σ is a non-integer rational number is analyzed in a similar manner as the case where σ is an integer number.

Error estimates in some discrete energy norm and L

2

. In consequence we have the following discrete boundary conditions: [18] P. Omnes, A second-order convergence of a function reconstructed from finite volume approximation of the Laplace equation on Delaunay-Voronoi meshes, ESAIM:M2AN, Vol.45, N o 4(2011),pp 627-650.