Numerical simulation for two-phase flow in a porous medium

Abstract

In this paper, we introduce a numerical study of the hydrocarbon system used for petroleum reservoir simulations. This system is a simplified model which describes a two-phase flow (oil and gas) with mass transfer in a porous medium, which leads to fluid compressibility. This kind of flow is modeled by a system of parabolic degenerated non-linear convection-diffusion equations. Under certain hypotheses, such as the validity of Darcy’s law, incompressibility of the porous medium, compressibility of the fluids, mass transfer between the oil and the gas, and negligible gravity, the global pressure formulation of Chavent (Mathematical Models and Finite Elements for Reservoir Simulation: Single Phase, Multiphase and Multicomponent Flows Through Porous Media, 1986) is formulated. This formulation allows the establishment of theoretical results on the existence and uniqueness of the solution (Gasmi and Nouri in Appl. Math. Sci. 7(42):2055-2064, 2013). Furthermore, different numerical schemes have been considered by many authors, among others we can refer the reader to (Chen in Finite element methods for the black oil model in petroleum reservoirs, 1994; Chen in Reservoir Simulation: Mathematical Techniques in Oil Recovery, 2007) and (Gagneux et al. in Rev. Mat. Univ. Complut. Madr. 2(1):119-148, 1989). Here we make use of a scheme based on the finite volume method and present numerical results for this simplified system.

1 Introduction

The fluids flow within porous media has an important role in various domains, such as geothermal studies, geotechnics (the mechanics of soils), chemical engineering, ground water storage, hydrocarbon exploitation (see references [1] and [2]), etc. In some cases, there are two or more fluids with different characteristics. We are concerned with the simulation of the displacement of a fluid by another one, within a porous medium, while the displacing fluid is immiscible with the fluid being displaced. Different numerical techniques, mainly based on finite elements, have been used by many authors to solve such problems, for example see [3] and [4].

In this paper we introduce a finite volume method for solving the hydrocarbon system model often used for petroleum reservoir simulations. To prevent consistency defects in the scheme, we propose to modify the mesh where the discontinuity occurs, because of the porous media. We propose to design our new scheme, called the ‘diamond meshes scheme’ (DMS), whose convergence is proved, and which can be used to solve the nonlinear discrete equations. Finally, numerical simulations confirm that the DMS is significantly useful for such difficult problems taking into account their physical properties.

In this case, there are only two chemical species, or components, gas and oil, where this last component may exist in both phases (gas and oil), that is to say, the heavy residual component in the oil phase and the light volatile component in the gas phase. In order to reduce confusion, we need to distinguish carefully between the ‘oil component’ and the ‘oil phase’. This model, called a hydrocarbon system, is a simplified compositional model describing two-phase flow in a porous medium with mass interchange between them. Therefore it can predict compressibility and mass transfer effects, in the sense that it is assumed that there is mass transfer between the oil and the gas phase. In this model the ‘oil component’ (also called stock-tank oil) is the residual liquid at atmospheric pressure left after a differential vaporization, while the ‘gas component’ is the remaining fluid.

2 Mathematical model

One of the fluids wets the porous medium more than the other; we refer to this as the wetting phase fluid and we refer to the other as the non-wetting phase fluid. In an oil-gas system, oil is the wetting phase. Let us consider a bounded connex open Ω of $R^{d}$ ($d=2\mbox{ or }3$), describing the porous medium (the reservoir), with a Lipchitz boundary Γ, and let t be the time variable t in $[0,T[$, $T <\infty$. Let $C_{Gg}$ be the mass fraction of the gas component in the gas phase, $C_{Og}$ the mass fraction of the oil component (the light component) in the gas phase, and $C_{Oo}$ the mass fraction of the oil component (the heavy component) in the oil phase which is equal to 1. While this distribution of the hydrocarbon components between the oil and gas phases plays an important role in a steam drive process, we cannot say that the mass of each phase is conserved because of the possibility of transfer of the oil component between the two phases. Instead, we observe that the total mass of each component must be conserved.

Then the mass flux of the oil and the gas components are

$$\begin{aligned}& C_{Og}\rho_{g}U_{g}+\rho_{o}U_{o}, \\& C_{Gg}\rho_{g}U_{g}. \end{aligned}$$

The oil and gas mass components are

$$\begin{aligned}& \phi ( C_{Og}\rho_{g}S_{g}+ \rho_{o}S_{o} ) , \\& \phi C_{Gg}\rho_{g}S_{g}. \end{aligned}$$

Then, for each fluid, we can write the conservation equations as

$$ \begin{aligned} &\frac{\partial}{\partial t} \bigl[ \phi ( C_{Og} \rho_{g}S_{g}+ \rho _{o}S_{o} ) \bigr] +\nabla\cdot ( C_{Og} \rho_{g}U_{g}+\rho _{o}U_{o} ) =0, \\ &\frac{\partial}{\partial t} \bigl[ \phi ( C_{Gg}\rho _{g}S_{g}) \bigr] +\nabla\cdot ( C_{Gg}\rho_{g}U_{g}) =0, \end{aligned} $$

(1)

where ($i=o,g$) $S_{i}$, $U_{i}$, $\rho_{i}$ represent the saturation, the velocity, and the density of the phase i, respectively. The parameter ϕ is the porosity of the medium. We have

$$\begin{aligned} &C_{Og}+C_{Gg} =1, \\ &\rho_{g} =f_{1} ( P_{g},C_{Gg} ) , \\ &\rho_{o} =f_{2} ( P_{o} ), \end{aligned}$$

where $P_{o}$ and $P_{g}$ are the oil and gaz pressure respectively.

We consider compressible fluids, with constant dynamic viscosities and where the gravity effect is neglected. Under these hypotheses, Darcy’s law combined with the mass conservation equations for each one of the component leads to the following system of partial differential equations of parabolic convection-diffusion type:

$$\begin{aligned}& \phi(x)\frac{\partial}{\partial t} \bigl(\rho_{g}S_{g}+ \rho_{o}\omega _{o}^{l}S_{o} \bigr)+ \nabla\cdot \bigl(\rho_{g}U_{g}+\rho_{o} \omega_{o}^{l}U_{o} \bigr) =0, \end{aligned}$$

(2)

$$\begin{aligned}& \phi(x)\frac{\partial}{\partial t} \bigl(\rho_{o}\omega _{o}^{h}S_{o} \bigr)+\nabla \cdot \bigl(\rho_{o}\omega_{o}^{h}U_{o} \bigr) =0, \end{aligned}$$

(3)

$$\begin{aligned}& U_{o} =-K(x)\frac{k_{ro}}{\mu_{o}}\nabla P_{o}, \end{aligned}$$

(4)

$$\begin{aligned}& U_{g} =-K(x)\frac{k_{rg}}{\mu_{g}}\nabla P_{g}, \end{aligned}$$

(5)

where ($i=o,g$) $P_{i}$, $\mu_{i}$, $k_{ri}$ represent the pressure, the viscosity, and the relative permeability of the phase i, respectively. The parameter K is the absolute permeability of the medium and $\omega_{o}^{c}$, $c=h,l$ is the mass fraction of the component c, denoted by h for the heavy component and by l for the light component in the oil phase. We have

$$\begin{aligned}& \mu_{g} =f_{3} ( P_{g},C_{Gg} ) , \\& \mu_{o} =f_{4} ( P_{o} ) , \\& k_{rg} =f_{5} ( S_{g},S_{o} ) , \\& k_{ro} =f_{6} ( S_{g},S_{o} ) . \end{aligned}$$

Furthermore, we shall use the subscript S to indicate standard conditions, i.e. appropriate conditions for the temperature and the pressure of medium.

Let $\rho_{OS}$, $\rho_{GS}$ be the density (measured at standard conditions) of the oil and the gas components, respectively. The gas formation volume factor, denoted by $B_{G}$, is the ratio of the volume of free gas (all of which is gas component), measured at the reservoir conditions, to the volume of the same gas measured at standard conditions. Thus

$$B_{G} ( P,T ) =\frac{V_{G} ( P,T ) }{V_{GS}}. $$

Let $W_{G}$ be the weight of free gas, since $V_{G}=\frac{W_{G}}{\rho_{g}}$ and $V_{GS}=\frac{W_{G}}{\rho_{gS}}$, then

$$B_{G}=\frac{\rho_{gS}}{\rho_{g}}, $$

so that the volatility of the oil in the gas is expressed by the ratio

$$R_{V}=\frac{V_{OS}}{V_{GS}}. $$

The mass fractions of the two components in the gas phase are

$$\begin{aligned}[b] &C_{Og} =\frac{R_{V}\rho_{oS}}{B_{G}\rho_{g}}, \\ &C_{Gg} =\frac{\rho_{gS}}{B_{G}\rho_{g}}. \end{aligned} $$

By adding the last two equations and noting that $C_{Gg}+C_{Og}=1$, we obtain

$$\rho_{g}=\frac{ ( \rho_{GS}+R_{V}\rho_{OS} ) }{B_{G}}. $$

The substitution of these mass fractions and densities into (1) gives, for the gas and the oil components,

$$ \begin{aligned} &\frac{\partial}{\partial t} \biggl[ \phi \biggl( \frac{R_{V}}{B_{G}}\rho _{oS}S_{g}+\rho_{o}S_{o} \biggr) \biggr] +\nabla\cdot \biggl( \frac {R_{V}}{B_{G}}\rho_{oS}U_{g}+ \rho_{o}U_{o} \biggr) =0, \\ &\frac{\partial}{\partial t} \biggl[ \phi \biggl( \frac{\rho _{gS}S_{g}}{B_{G}} \biggr) \biggr] +\nabla\cdot \biggl( \frac{\rho_{gS}U_{g}}{B_{G}} \biggr) =0. \end{aligned} $$

(6)

We suppose that it is a saturated regime and is expressed by

$$ S_{o}+S_{g}=1. $$

(7)

The capillary pressure is given by

$$ P_{g}-P_{o}=P_{c}(S_{o})=p_{c}(S_{o})p_{cM}, $$

(8)

where

$$ p_{cM}=\sup \bigl\vert P_{c}(S_{o}) \bigr\vert \quad\mbox{and}\quad0\leq p_{c}(S_{o})\leq1. $$

(9)

We define the mobility of each phase by the formula

$$ \lambda_{i}=\frac{k_{ri}}{\mu_{i}},\quad i=o,g, $$

(10)

and the total mobility λ by

$$ \lambda=\lambda_{o}+\lambda_{g}. $$

(11)

To simplify, we set

$$\begin{aligned}& \rho_{o}^{h} =\rho_{o}\omega_{o}^{h}, \end{aligned}$$

(12)

$$\begin{aligned}& \rho=\rho_{g}+\rho_{o}, \end{aligned}$$

(13)

$$\begin{aligned}& b =\rho_{g}\lambda_{g}+\rho_{o} \lambda_{o}, \end{aligned}$$

(14)

$$\begin{aligned}& d =\rho_{g}-\rho_{o}. \end{aligned}$$

(15)

Let us now introduce the new unknowns, namely the reduced saturation and the global pressure in the following way: if we denote by $S_{i,m}$ and $S_{i,M}$, the residual and the maximum saturations of the fluid $i=o,g$, respectively; the reduced saturation is given by

$$\begin{aligned}& S =\frac{S_{o}-S_{o,m}}{1-S_{g,m}-S_{o,m}}, \end{aligned}$$

(16)

$$\begin{aligned}& 0 \leq S\leq1. \end{aligned}$$

(17)

The ‘global pressure’ was first introduced by Chavent and Jaffre [1] in the following form:

$$ P=\frac{1}{2} ( P_{g}+P_{o} ) +\gamma ( S ) , $$

(18)

with

$$ \gamma ( S ) =\frac{1}{2}\int_{S_{o,m}}^{S} \frac {\lambda _{g}-\lambda_{o}}{\lambda}p_{c}^{\prime} ( \xi ) p_{cM}\,d\xi. $$

(19)

Hence

$$ \gamma(S)=\frac{1}{2}\int_{0}^{S}\alpha ( \xi )\,d\xi, $$

(20)

where

$$ \alpha ( S ) =\frac{\lambda_{g} ( S ) -\lambda _{o} ( S ) }{\lambda ( S ) }p_{c}^{\prime} ( S ) P_{cM} $$

(21)

is the capillary diffusion. Therefore, our model is given by the following simplified system:

$$\begin{aligned}& \begin{aligned}[b] &\Phi(x)\frac{\partial}{\partial t} \bigl( \rho_{o}^{h}S \bigr) +\phi (x)S_{o,m}\frac{\partial}{\partial t} \bigl( \rho_{o}^{h} \bigr) \\ &\quad{}-\nabla \cdot \bigl( K(x)\rho_{o}^{h}\lambda_{o}(S) \nabla P \bigr) +\nabla\cdot \bigl( K(x)\rho _{o}^{h}\alpha (S) \nabla(S) \bigr) =0, \end{aligned} \end{aligned}$$

(22)

$$\begin{aligned}& \begin{aligned}[b] &\Phi(x)\frac{\partial}{\partial t} ( \rho S ) +\phi(x)\frac{ \partial}{\partial t} ( \rho S_{o,m}+\rho_{g} ) \\ &\quad{}-\nabla \cdot \bigl(K(x)b(S,P)\nabla P \bigr)+ \nabla\cdot \bigl( K(x)d ( P ) \alpha(S)\nabla S \bigr) =0, \end{aligned} \end{aligned}$$

(23)

$$\begin{aligned}& \lambda ( S ) \nabla P\cdot\eta=0, \qquad\alpha(S)\nabla S=0, \quad\mbox{on }\Gamma\times ( 0,T ) , \end{aligned}$$

(24)

$$\begin{aligned}& S ( x,0 ) =S^{0} ( x ),\qquad P ( x,0 ) =P^{0} ( x ) ,\quad\mbox{in }\Omega, \end{aligned}$$

(25)

where $\phi(x)=\varphi(x)(1-S_{o,m}-S_{o,g})$.

3 Weak formulations

First we introduce the following spaces:

$$\begin{aligned}& H ( \operatorname{div},\Omega ) = \bigl\{ v\in \bigl( L^{2} \bigl( \Omega , (0,T ) \bigr) \bigr) ^{d},\operatorname{div} ( v ) \in L^{2} \bigl( \Omega, ( 0,T )\bigr) ,d=2,3 \bigr\} , \end{aligned}$$

(26)

$$\begin{aligned}& \begin{aligned} &V ( \Omega ) = \bigl\{ v\in H ( \operatorname{div},\Omega ) ,v\cdot\eta =0\mbox{ on } \Gamma \bigr\} ,\\ &W ( \Omega ) = \bigl\{ v\in V ( \Omega ) ,v ( x,T ) =0\mbox{ in }\Omega \bigr\} . \end{aligned} \end{aligned}$$

(27)

The weak formulation of problem (22)-(25) is written as

$$\begin{aligned}& \biggl( \Phi(x)\rho_{o}^{h}S,\frac{\partial v}{\partial t} \biggr) _{\Omega}- \bigl( K(x)\rho_{o}^{h} \lambda_{o}(S)\nabla P,\nabla v \bigr) _{\Omega}+ \bigl( K(x) \rho_{o}^{h}\alpha(S)\nabla S,\nabla v \bigr) _{\Omega}= ( f_{1},v ), \end{aligned}$$

(28)

$$\begin{aligned}& \biggl( \Phi(x)\rho S,\frac{\partial v}{\partial t} \biggr) _{\Omega }- \bigl( K(x)b(S,P)\nabla P,\nabla v \bigr) _{\Omega}+ \bigl( K(x)d ( P ) \alpha(S) \nabla S,\nabla v \bigr) _{\Omega}= ( f_{2},v ), \end{aligned}$$

(29)

where $(\cdot,\cdot)_{\Omega}$ is the inner product defined on $W ( \Omega )$. A theoretical study of the existence and uniqueness of the weak solution was done and the details of the results can be found in [5].

4 Finite volume approximation

In this section we study the approximation of solutions to our model in the discrete finite volume framework. This family of schemes allows very general meshes and deals with the main properties of the physical features of the treated problems. The time interval $[ 0,T [$ is divided into finite sub-intervals $[ t_{n},t_{n+1} [$ of length $\Delta t_{n}$, $n=0,\ldots,M$ with $t_{0}=0$ and $t_{M}=T$. The space domain (the reservoir Ω) is discretized by a non-structured stitching $T_{h}$.

4.1 Notations

We introduce the following notation:

• Let $\vert K\vert$ denote the cell K surface, $N(K)$ the set of triangles having a common side with the cell K.

• Let $e_{K,L}$ be the common side of the triangles K and L and $\overrightarrow{\eta}_{K,L}$ be the normal oriented from K towards L.

• $\overrightarrow{\eta}_{e_{i}}$ is the external normal corresponding to the part of $e_{i}$ at the boundary Γ.

• Let $S_{h}$ be the set of sides of the stitching $T_{h}$ and $S_{h}^{\ast}$ be the set of the interior sides.

• For a given side e, let us denote by S and N the extremities e, by W and E the two triangles where $e=\bar{W}\cap\bar{E}$; by $\chi _{e}$ the diamond cell associated with e given by connecting the centers of gravities of the cells W and E with the extremities N and S of e.

• $( ( \varepsilon_{_{i}} ) _{i=1,4} )$ are the four segments forming the border of $\chi_{e}$.

• $\overrightarrow{\eta}_{\varepsilon}=\frac{1}{\vert \varepsilon _{_{i}}\vert } ( \mu_{x_{i}},\mu_{y_{i}} )$ is the normal on $\varepsilon_{_{i}}$ outgoing of $\chi_{e}$.

• For a given node, $V ( N )$ is the set of triangles with this node in common.

For the numerical resolution of problem (22)-(25), the first two equations will be discretized separately. For more details on finite volume methods, see for example [6] and [7].

4.2 Discretization of the first equation

Let C be a cell of the stitching $T_{h}$, at time $t_{n}$; we integrate (28) on C to get

$$\begin{aligned} &\frac{1}{\Delta t}\int_{C} \bigl[\Phi(x) \rho_{o}^{h,n} \bigl( S^{n+1}-S^{n} \bigr) \bigr]\,dx \\ &\qquad{}+\frac{1}{\Delta t}\int_{C} \bigl[\Phi(x) \bigl( \rho_{o}^{h,n} \bigr)^{\prime}\ S^{n}\ \bigl( P^{n+1}-P^{n} \bigr) \bigr]\,dx \\ &\qquad{}+\frac{1}{\Delta t}\int_{C} \bigl[\phi(x) S_{o,m} \bigl(\rho _{o}^{h,n} \bigr)^{\prime } \bigl( P^{n+1}-P^{n} \bigr) \bigr]\,dx \\ &\qquad{}-\sum_{D\in N ( C ) }\int_{e_{CD}} \bigl[K(x) \rho_{o}^{h,n}\ \lambda _{o}^{n}(S) \nabla P^{n} \eta_{e} \bigr]\,dx \\ &\qquad{}+\sum_{D\in N ( C ) }\int_{e_{CD}} \bigl[K(x) \rho_{o}^{h,n}\ \alpha ^{n}(S) \nabla \bigl(S^{n} \bigr) \eta_{e} \bigr]\,dx \\ &\quad=0. \end{aligned}$$

(30)

Therefore we obtain

$$\begin{aligned} &\frac{\vert C\vert }{\Delta t}\Phi_{C} \rho _{o,C}^{h,n} \bigl( S_{C}^{n+1}-S_{C}^{n} \bigr) \\ &\qquad{}+\frac{\vert C\vert }{\Delta t}\Phi_{C} \bigl(\rho _{o,C}^{h,n} \bigr)^{\prime} S_{C}^{n} \bigl( P_{C}^{n+1}-P_{C}^{n} \bigr) \\ &\qquad{}+\frac{\vert C\vert }{\Delta t}\phi_{C} S_{o,m} \bigl(\rho _{o,C}^{h,n} \bigr)^{\prime} \bigl( P_{C}^{n+1}-P_{C}^{n} \bigr) \\ &\qquad{}-\sum_{D\in N ( C ) }\int_{e_{CD}} \bigl[K \rho_{o}^{h,n} \lambda _{o}^{n}(S) \nabla P^{n} \eta_{e} \bigr]\,dx \\ &\qquad{}+\sum_{D\in N ( C ) }\int_{e_{CD}} \bigl[K \rho_{o}^{h,n} \alpha ^{n}(S)\nabla \bigl(S^{n} \bigr) \eta_{e} \bigr]\,dx \\ &\quad=0. \end{aligned}$$

(31)

This implies that

$$\begin{aligned} &\frac{\vert C\vert }{\Delta t}\Phi_{C} \rho_{o,C}^{h,n} S_{C}^{n+1}+\frac{\vert C\vert }{\Delta t} \bigl(\rho _{o,C}^{h,n} \bigr)^{\prime} \bigl( \Phi_{C} S_{C}^{n}+ \phi_{C}\ S_{o,m} \bigr) P_{C}^{n+1} \\ &\quad=- \frac{\vert C\vert }{\Delta t}\Phi_{C} \rho _{o,C}^{h,n}S_{C}^{n}-\frac{\vert C\vert }{\Delta t} \bigl(\rho_{o,C}^{h,n} \bigr)^{\prime } \bigl( \Phi_{C}S_{C}^{n}- \phi_{C}S_{o,m} \bigr) P_{C}^{n} \\ &\qquad{}+\sum _{D\in N ( C ) }F_{C,e}^{n}- \sum _{D\in N ( C ) }F_{D,e}^{n}, \end{aligned}$$

(32)

where $(\rho_{o,C}^{h,n})^{\prime}$ is the space derivative of $(\rho _{o,C}^{h,n})$, and $F_{C,e}^{n}$ and $F_{D,e}^{n}$ are the numerical flows of convection and diffusion, which are approximated by

$$\begin{aligned}& F_{C,e}^{n} =K_{e}\rho_{o,e}^{h,n} \lambda_{o,e}^{n}\nabla _{e}P^{n} \vert e\vert , \end{aligned}$$

(33)

$$\begin{aligned}& F_{D,e}^{n} =K_{e}\rho_{o,e}^{h,n} \alpha_{e}^{n}\nabla _{e}S^{n}\vert e\vert , \end{aligned}$$

(34)

where e is a side of the stitching, $\overrightarrow{\eta}_{e}$ is the normal of e outgoing from C; $K_{e}$, $\rho_{o,e}^{h,n}$, and $\lambda _{o,e}^{n}$ denote the interpolations on e of the functions K, $\rho _{o}^{h}$, and $\lambda_{o}$, respectively, while $\nabla_{e}P^{n}$ is the approximation of the gradient of the pressure on the side e. The construction of the approached gradient on e is done by the so-called Green-Gauss method. It consists of approaching the gradient by its means on a co-volume in the form of a diamond constructed around the side e. We write

$$\begin{aligned}& \rho_{o,e}^{h,n}\lambda_{o,e}^{n}= \rho_{o}^{h} \bigl( P_{e}^{n} \bigr) \lambda_{o} \bigl( S_{e}^{n} \bigr) , \end{aligned}$$

(35)

$$\begin{aligned}& \rho_{o,e}^{h,n}\alpha_{e}^{n}= \rho_{o}^{h} \bigl( P_{e}^{n} \bigr) \alpha \bigl(S_{e}^{n} \bigr), \end{aligned}$$

(36)

and

$$ \nabla_{e}P=\frac{1}{\vert \chi_{e}\vert }\sum_{\varepsilon\in\partial\chi_{e}} \frac{1}{2} ( P_{N_{1} ( \varepsilon ) }+P_{N_{2} ( \varepsilon ) } ) \vert \varepsilon \vert \overrightarrow{\eta}_{e}, $$

(37)

where $\chi_{e}$ is the diamond cell associated with e, and $N_{1} ( \varepsilon )$ and $N_{2} ( \varepsilon )$ are the extremities of ε, one of the four segments forming $\partial \chi_{e}$ (the boundary of $\chi_{e}$) and $\overrightarrow{\eta }_{e}$ is the unit normal vector. The values of P at the centers W and E are $P_{W}$ and $P_{E}$, while the values at nodes N and S are interpolated on the boundary and denoted by $P_{N}$ and $P_{S}$. For convenience, we omit the indication η. Hence at each node N we have

$$ P_{N}=\sum_{K\in V ( N ) }y_{K} ( N ) P_{K}, $$

(38)

where $V ( N )$ is the set of triangles sharing the node N, $P_{K}$ is the value of P at the center of the cell K and $y_{K} ( N )$ are the interpolation weights. The gradient of the saturation is given by

$$ \nabla_{e}S=\frac{1}{\vert \chi_{e}\vert }\sum_{\varepsilon\in\partial\chi_{e}} \frac{1}{2} ( S_{N_{1} ( \varepsilon ) }+S_{N_{2} ( \varepsilon ) } ) \vert \varepsilon \vert \overrightarrow{\eta}_{e}, $$

(39)

with

$$ S_{N}=\sum_{K\in V ( N ) }y_{K} ( N ) S_{K}. $$

(40)

4.3 Discretization of the second equation

In the same way, integration of (29) implies

$$\begin{aligned} &\frac{\vert C\vert }{\Delta t}\Phi_{C}\rho _{C}^{n}S_{C}^{n+1}+\frac{\vert C\vert }{\Delta t} \bigl( \Phi_{C} \bigl(\rho _{C}^{h,n} \bigr)^{\prime}S_{C}^{n}+\phi_{C} \bigl( \rho_{C}^{h,n} \bigr)^{\prime }S_{o,m}+ \phi_{C} \bigl(\rho_{C}^{h,n} \bigr)^{\prime} \bigr) P_{C}^{n+1} \\ &\quad=-\frac{\vert C\vert }{\Delta t}\Phi_{C}\rho_{C}^{n}\ S_{C}^{n}-\frac{\vert C\vert }{\Delta t}\left ( \Phi_{C}(\rho_{o,C}^{h,n})^{\prime}S_{C}^{n}-\phi_{C}(\rho _{o,C}^{h,n})^{\prime}S_{o,m}+ \phi_{C}(\rho_{o,C}^{h,n})^{\prime} \right ) P_{C}^{n} \\ &\qquad{}+\sum _{D\in N ( C ) }G_{C,e}^{n}-\sum_{D\in N ( C ) }G_{D,e}^{n}, \end{aligned}$$

(41)

where $G_{C,e}^{n}$ and $G_{D,e}^{n}$ are the numerical flows of the convection and diffusion approximated by

$$\begin{aligned}& G_{C,e}^{n} =K_{e}b_{e}^{n} \nabla_{e}P^{n}\vert e\vert , \end{aligned}$$

(42)

$$\begin{aligned}& G_{D,e}^{n} =K_{e}d_{e}^{n} \alpha_{e}^{n}\nabla_{e}S^{n}, \end{aligned}$$

(43)

with

$$\begin{aligned}& b_{e}^{n} =b_{e}^{n} \bigl( S_{e}^{n},P_{e}^{n} \bigr) , \end{aligned}$$

(44)

$$\begin{aligned}& d_{e}^{n}\alpha_{e}^{n} =d_{e}^{n} \bigl( P_{e}^{n} \bigr) \alpha \bigl(S_{e}^{n} \bigr), \end{aligned}$$

(45)

where e is a side of the stitching that limits the cell C, $\nabla _{e}P^{n}$ is given by (42) and (44), and $\nabla_{e}S^{n}$ is given by (43) and (45).

5 Numerical experiments

The problem given in (22)-(25) was said not to be typical of a hydrocarbon system simulation, but it can be designed to test the resolution capabilities of the numerical method for problems involving sharp fronts.

The numerical tests are done on a homogeneous $isotropic$ reservoir. The physical setting in 2-D was as follows. A horizontal rectangular reservoir $\Omega=\,]0,300[\,\times\,]0,60[$, discretized with a structured mesh made up of $3{,}000$ cells, with an intrinsic permeability of 0.001 was initially filled with a mixture of gas and oil, their residual saturations are 0.02 and 0.03, respectively. The initial fluid (oil) distribution was taken to be uniform on the whole reservoir surface, and it had an associated initial pressure $P_{0}=2{,}000~\mbox{bar}$. The porosity is $\Phi=0.02$. The mobilities and the capillary pressure are given by $P_{c}=(s^{2}-1)/2$, $\lambda _{o}=0.5$, $\lambda_{g}= 0.2$, and $\mu_{w}=1$, and $\mu_{o}=3$, the time step is $\Delta t=10~\mbox{s}$. These tests were carried out using the free and open-source simulator for flow and transport processes in porous media, based on the Distributed and Unified Numerics Environment DUNE (see www.dumux.org).

5.1 Results and discussion

The numerical results shown in Figures 1-3 give the pressure for the wetting and non-wetting fluid, and the saturation for the non-wetting one. Because of the complexity of the problem, we have introduced new unknowns, namely the global pressure and the reduced saturation, in order to reduce the number of unknowns from six to two (for more details see [5]). An explicit non-structured finite volume scheme has been used to solve this simplified problem with the new unknowns, and because of the non-structured meshes, we proposed a method based on ‘diamond cells’ to approximate the gradient. These results show that the scheme is very stable.

Figure 1

Non-wetting pressure at $\pmb{t=0, 50, 100, \mbox{and }175}$ .

Figure 2

Wetting pressure at $\pmb{t=0, 50, 100,\mbox{ and }175}$ .

Figure 3

Non-wetting pressure at $\pmb{t=0, 50, 100,\mbox{ and }175}$ .

6 Conclusion

In this paper we proposed a new scheme based on the finite volume method for solving the displacement of a fluid by another one, within a porous medium, while the displacing fluid is immiscible with the fluid being displaced. This gives the hydrocarbon system model often used for petroleum reservoir simulations. To prevent consistency defects in the scheme, we suggested to modify the mesh where the discontinuity occurs, because of the porous media. This convergent scheme, called the ‘diamond meshes scheme’ (DMS) was designed, to solve the associated nonlinear discrete equations. Finally, the numerical results, which are linked to our theoretical ones in [5], confirm that the DMS is significantly useful for such difficult problems (see Figures 1-3).