Regular Article
Study of compositional multiphase flow formulation using complementarity conditions
IFP Energies nouvelles, Applied Mathematics Department, 14, avenue de BoisPréau, 92852
RueilMalmaison Cedex, France
^{*} Corresponding author: eric.flauraud@ifpen.fr
Received:
6
March
2018
Accepted:
1
March
2019
In this article, two formulations of multiphase compositional Darcy flows taking into account phase transitions are compared. The first formulation is the socalled natural variable formulation commonly used in reservoir simulation, the second has been introduced by Lauser et al. and uses the phase pressures, saturations and component fugacities as main unknowns. We will discuss how the Coats and the Lauser approaches can be used to solve a compositional multiphase flow problem with cubic equations of state of Peng and Robinson. Then, we will study the results of several synthetic cases that are representative of petroleum reservoir engineering problems and we will compare their numerical behavior.
© I. Ben Gharbia & E. Flauraud, published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
In reservoir simulations, compositional multiphase flow coupled with detailed physical laws based on equations of state [1–4] remains an important and challenging problem. This issue consists in solving a large system of nonlinear equations coupling the conservation of the components with the thermodynamical equilibrium constraints. Many formulations have been proposed; we refer to, for example [5–9]. The difficulty lies in handling the appearance and disappearance of phases assumed to be at thermodynamical equilibrium.
The traditional dynamic approach [10], called “variable switching”, considers only the unknowns of the present phases and the equations for them. It is heavy to implement and costly in CPU time, because the “switching” occurs constantly, even from a Newton iteration to another. An alternative approach, called “unified formulation”, enables one to keep a fixed set of unknowns and equations regardless of the context using nonlinear complementarity conditions [11, 12]. These constraints provide elegant models for complex problems and lead to efficient methods to solve them numerically. Interesting results have already been obtained in several fields like solid or fluid mechanics and economics[12, 13]. Thus, this paper will concentrate on the physical validation of this new approach in reservoir simulation industry comparing two formulations for compositional multiphase flow with cubic equations of state.
The first formulation we study is from the variable switching approach, called Natural Variables Formulation (NVF) and introduced by Coats [5, 10]. It is based on the natural unknowns (pressures, saturations, molar fractions) and on phase apparition detection through a flash calculation [14]. This model is stable in regard to phase transitions but can quickly become complex to manage the set of present phases and the associated unknowns/equations at each point of the timespace domain.
For the second formulation, we examine a recent unified formulation introduced by Lauser et al. [15], which uses the phase pressures, saturations and component fugacities as main unknowns with complementarity conditions for handling phase transitions. These complementarity constraints can be equivalently reformulated as Karush–Kuhn–Tucker (KKT) conditions of Gibbs free energy minimization problem and this is known as the Gibbs tangent plane. The resulting system leads to a fixed set of unknowns and a fixed set of equations whatever the present phases, and it allows to avoid the flash calculation. On the other hand, it requires to compute molar fractions and their derivatives as a function of pressures, temperature and fugacities. From the practical viewpoint, as the new formalism involves several nonsmooth “complementarity” equations, it is necessary after discretization to resort to “semismooth” Newton methods, called Newtonmin algorithm [16–18]. This Complementarity Condition Formulation will be denoted by CCF in the following.
The paper is organized as follows. In Section 2, we describe the system of equations for the multiphase compositional model, relevant physical properties along with their dependencies and relations of thermodynamical equilibrium between phases. In Sections 3 and 4, the two formulations are detailed and their advantages and drawbacks are further discussed. In Section 5, we present some numerical aspects of the resolution of both systems and in particular for the nonlinear solver. The specific treatment of aquifer cells is presented in Section 6. Finally, in the last Section 7, the two formulations are compared into several synthetic cases that are representative of petroleum reservoir engineering problems.
2 Multiphase compositional model
For the sake of simplicity, we assume the medium is isotherm with fixed temperature T, hence in the following the dependence of the physical laws on temperature will not be shown. We consider a compositional model for a threephase flow: water (w), oil (o) and gas (g). Let A = {w, o, g} be the set of phases. We suppose that the oil and gas phases are represented by a mixture of n _{ c } components. Moreover, we assume that the water phase is pure, that means that it is composed only of H_{2}O, and that it is immiscible with the other phases. Therefore only the hydrocarbons oil and gas phase are mixable and compositional.
2.1 Conservation laws
The governing partial differential equations are obtained by enforcing the molar conservation law for each component. In particular, the velocity of each phase is given by Darcy’s law. For the water component, we have
(1)and, for each hydrocarbon component i = 1,…,n _{ c }, we have
where:

ϕ: porosity of the medium [–],

S _{ α }: saturation of phase α ∈ A [–],

ρ _{ α }: molar density of phase α ∈ A ,

: molar fraction of component i in phase α ∈ A/{w} [–],

: Darcy’s velocity of phase α ∈ A ,

δ: the Dirac delta function ,

q _{ w }, q _{ i } are the molar flow rates of each component produced or injected at the well .
The molar flow rates have the following form
where Q _{ α } represents the flow rate of phase α and depends on the nature of the associated well. The velocities of the phases are computed through Darcy’s law

K: permeability of the medium [m^{2}],

S _{ α }: saturation of phase α ∈ A [–],

: relative permeability of phase α ∈ A [–],

μ _{ α }: viscosity of phase α ∈ A [Pa s],

P _{ α }: pressure of phase α ∈ A [Pa],

: volumetric mass density of phase α ∈ A ,

: gravity .
The units of the physical quantities defined above are given in the international unit system.
The thermodynamical aspect is constantly present in the simulations of flows in a porous medium. It is necessary, for example, to be able to describe the physical evolution of the mixtures present in the reservoir. This is important not only for being able to simulate flows as well as possible but also for determining which recovery process is best suited. In this section, we will introduce some notions to better understand the thermodynamical aspect involved.
2.2 Equilibrium equations
In addition to the molar conservation equations, we have to consider the relations of thermodynamical equilibrium between different phases. For every component, it can be expressed as the equality of the fugacity of oil and gas phase, that is:
Let be the vector of components molar fractions in phase α, then the fugacities can be computed through
where P is a reference pressure and corresponds to the pressure of a given phase α and is the fugacity coefficient of the component i in the phase α ∈ {o, g}. The fugacity coefficient is computed using an equation of state [19]. We can rewrite equation (3) as
where K _{ i } = K _{ i } (P, C ^{ o }, C ^{ g }) is the equilibrium constant, it is defined as the ratio between the fugacity coefficients of oil and gas phase
and it represents the distribution of component i between these two phases.
2.2.1 Negative flash calculation
The solution of the thermodynamical equilibrium (3) is computed using a negative flash calculation [14]. We define the total molar fraction of the components in the hydrocarbon mixture such as
For a given reference pressure P and a fixed total molar fraction Z, the negative flash is to find the molar fraction of the gas phase and , the molar fractions of the components in the oil and gas phase at equilibrium, which are the solutions of the nonlinear system
The solutions V ≥ 1 and V ≤ 0 correspond to a stable single phase system. In the following, we denote by Flash (P, Z) the negative flash calculation.
2.2.2 From phase stability test to complementarity condition
The necessary and sufficient condition for a phase of a n _{ c } component mixture with composition to be stable at some given pressure and temperature [19, 20] is given by following inequality
(4)where the value of the vector containing the molar fractions of an incipient phase that would separate from the phase with composition Z in case of instability. The function D(C) is called the tangent plane distance function. Hence, the problem of the thermodynamical stability analysis consists in resolving the following constrained minimization problem
The chemical potential μ _{ i } [J mol^{−1}] of a component i is given by

f _{ i }: fugacity of the component i [Pa],

P: reference pressure [Pa],

T: temperature [K],

R: universal gas constant [J mol^{−1} K^{−1}],

: chemical potential of pure component i in the state of perfect gas at the same pressure P and temperature T [J mol^{−1}].
Using the fugacity coefficients
and the relation (6), the function (4) can be rewritten as follow
(7)We define the Lagrangian for the constrained problem (5) as
(8)where λ is the Lagrange multiplier. It is easy to write the optimality conditions from (8)
Replacing these optimality conditions in the equation (7), one can show that λ is nonnegative. Introducing a variable transformation , we can proof that the stationary points indicate stability when
Consequently, this inequality will be used in the following as a term of the complementarity conditions (11) to represent the equilibrium of the oil and gas phases.
2.3 Closure equations
To solve any system of equations, we have to make sure that it is mathematically welldefined. This means that the equations are independent, i.e., that no equation can be expressed in term of the others, that a unique solution exists and that the number of unknowns is equal to the number of equations. In order to close the system, we need to add other equations. Firstly we can consider the relations between capillary pressure
Moreover, since there must be conservation of volume, we have
At last, if the oil phase (resp. gas phase) is present, the following relation holds
In the next section, two formulations are detailed and their advantages/drawbacks are discussed.
3 Natural variable formulation
The reservoir simulation industry commonly uses the formulation introduced by K.H. Coats [10] based on natural variables that are: pressures, saturations and components molar fractions and on the phase apparition detection by a flash calculation. The main advantage of this formulation is the use of the natural set of unknowns for the hydrodynamical and thermodynamical laws. On the other hand, its biggest disadvantage is that the set of unknowns and equations depends on the present phases. For this reason, we have to compute the set of present phases at each point of the spacetime domain and then switch unknowns and equation according to which phases are present. We use a negative flash calculation [14] and the saturation sign to determine phases appearances and disappearances respectively. Let A _{ p } be the set of present phases, thus, depending on the point in the timespace domain we are considering, we can have
We summarize below the principle of the treatment of the phase appearance and disappearance:

The test of phase appearance is done only for the contexts where saturated oil or gas phase is present (i.e., A _{ p } = {w,o} or A _{ p } = {w, g}). A negative flash calculation for the current reference pressure P and overall composition z returns the gas molar fraction V, then

if A _{ p } = {w, o} and V > 0 the oil phase becomes saturated and the gas phase appears. A _{ p } becomes A _{ p } = {w, o, g},

if A _{ p } = {w, g} and V < 1 the gas phase becomes saturated and oil phase appears. A _{ p } becomes A _{ p } = {w, o, g}.


The test of phase disappearance is done only for the context where both oil and gas phases are present (i.e. A _{ p } = {w, o, g}). We use the saturation’s sign:

if A _{ p } = {w, o, g} and S _{ o } ≤ 0, the oil phase disappears and then A _{ p } becomes A _{ p } = {w, g},

if A _{ p } = {w, o, g} and S _{ g } ≤ 0, the gas phase disappears and then A _{ p } becomes A _{ p } = {w, o}.

As in this model, we suppose that the water phase is immiscible with the other phases the context A _{ p } = {w} needs a particular treatment (see Sect. 6).
Hence, in the case of no aquifer cells and if we choose oil pressure as the reference pressure, the problem we have to solve is to find for every α ∈ A _{p}, P, S _{ α }, C ^{ α } satisfying the following system of equations
The following table shows the number of equations we have to solve depending on the context.
Context A _{ p }  {w, o, g}  {w, o}  {w, g} 
Conservation equations  n _{ c } + 1  n _{ c } + 1  n _{ c } + 1 
Closure equations  3  2  2 
Thermodynamical equilibrium  n _{ c }  0  0 
Total number of equations  2n _{ c } + 4  n _{ c } + 3  n _{ c } + 3 
Similarly, the unknowns of the system are given by the following table.
Context A _{ p }  {w, o, g}  {w, o}  {w, g} 
Pressure  1  1  1 
Saturations  3  2  2 
Molar fractions  2n _{ c }  n _{ c }  n _{ c } 
Total number of unknowns  2n _{ c } + 4  n _{ c } + 3  n _{ c } + 3 
We can notice that in all the cases the number of equations corresponds to the number of unknowns. In the next section, we present a new formulation introduced in [15], which is based on nonlinear complementarity conditions [12, 16, 21].
4 Complementarity condition formulation
Recently, new formulations have been introduced for modeling systems which have to switch continuously between different states. A new way of incorporating phase transitions into the simulation of multiphase, multicomponent processes in porous media have been developed. The idea is that a fluid phase appears or vanishes if a physical quantity exceeds a given threshold. Based on this observation we formulate the conditions for the local presence of fluid phases as a set of socalled complementarity conditions. This approach uses the fact that these constraints can be reformulated equivalently as a nondifferentiable but semismooth equation, called nonlinear complementarity equation. The complementarity problems were first manifested in optimality conditions for optimization problems [12]. The methods for the resolution of complementarity problems have spread in the last years following the work of the numerics [9, 16–18, 21, 22]. We focus on the formulation introduced by Lauser et al. [15] which uses the pressures, the phase saturations and the components fugacities as main unknowns. The biggest advantage of this formulation with respect to the previous one is the use of a fixed set of equations and unknowns independent on the present phases and, therefore, the flash calculation is avoided. The molar fractions C ^{ α } of the components are computed starting from the fugacities and the reference pressure. They are defined as the solution of the nonlinear system
If the phase α is present, the quantity coincides with the molar fraction C ^{ α }. On the contrary, if α is absent, represents the molar fraction that is at thermodynamical equilibrium with the ones of the present phase. In this case doesn’t have a physical meaning but it is still included in the system of equations. For this reason the quantities are called extended molar fractions. Thanks to the thermodynamical equilibrium, we have
and then the molar fractions of a component i can be computed from the same fugacity f _{ i } for both the phases. With this formulation, phase transitions are managed using complementarity conditions. The idea is that a phase α ∈ {o, g} is not present if its saturation is zero and the sum of its extended molar fraction is less than one. On the other hand, the phase is present if its saturation is greater than zero and the sum of its extended molar fractions is equal to one. Therefore for one phase α ∈ {o, g} we have the following relations:

If phase α is present S _{ α } > 0 and .

If phase α is absent S _{ α } = 0 and .
These constraints are equivalent to the following complementarity conditions
As said before, they can be reformulated using a complementarity function. In this work, we use the minimum function as complementarity function. Therefore the conditions (11) can be rewritten in the following way
Let A = {w, o, g} be the set of phases. Then, the problem we have to solve with this formulation is as follows. Find satisfying the following system of equations
In the following tables, we can find the number of unknowns and of equations of this formulation.
Conservation equations  n _{ c } + 1 
Closure equation  1 
Complementarity conditions  2 
Total number of equations  n _{ c } + 4 
Pressure  1 
Saturations  3 
Fugacities  n _{ c } 
Total number of unknowns  n _{ c } + 4 
We can notice that, in all the cases, the number of equations corresponds to the number of unknowns. Hence, the resulting system of mass conservation equations and equilibrium conditions is fully free of inequalities (pure set of equations).
5 Discretization and numerical resolution
Let be an admissible finite volume mesh of the reservoir given by a family of control volumes noted K. We also introduce an increasing sequence of discret times such that t ^{0} = 0 and t ^{ N } = T and we denote by Δt ^{ n } the time step such as Δt ^{ n } = t ^{ n }−t ^{ n−1} for n = 1,…, N. The systems (9) and (4) are discretized using fully implicit Euler integration in time and a cellcentered finite volume scheme with a twopoint flux discretization. The mobility terms are upwinded with respect to the sign of the phase Darcy flux [1]. The vectors of the discrete unknowns in each finite volume K and on the whole mesh are denoted respectively by and .
5.1 Natural variable formulation
For NVF, the set of discrete unknowns is depending on the context in each finite volume K and is given at each time t ^{ n } by:

if

if

if
We obtain at each time step n = 1, 2, …, N a nonlinear system denoted by
where is the set of discrete conservation equations and is the set of discrete closure equations and thermodynamical equilibrium equations depending in the context in each cell . In the following, we omit the subscript n for the sake of clarity. The nonlinear system (13) is solved at each time step by a Newton algorithm. At each Newton iteration (l), we have to solve the following linear system.
where is the Jacobian of the vector R _{ K } (resp. F _{ K }) and .
In order to reduce the size of the system (14) in a system of only n _{ c } + 1 equations per cell, we split the set of unknowns X _{ K } into n _{ c } + 1 primary unknowns and a set of remaining secondary unknowns . This splitting is done cell by cell depending on the set of present phases:

if A _{ p,K } = {w, o, g}

if A _{ p,K } = {w, o}

if A _{ p,K } = {w, g}
The linear system (14) can be rewritten in each cell as follows
where , and V(K) is the set of cells which are involved in the calculation of the flux through the edges of the cell K.
The secondary unknowns are eliminated using the closure equations
where is assumed to be non singular.
Finally, the linear system we have to solve at each Newton’s iteration (l) consists of the following equations
The solution of this system is used to update the next iterates as follow
Then, we update the context in each cell using the algorithm presented at the beginning of this section:

if

if then and

if then and


if

Flash

if then , and

else and


if

Flash

if then and

else and

The stopping criteria of the Newton algorithm is given by
where ⋅ is a given norm, ε a specified tolerance and is an initial guess for the Newton algorithm. In general, the initial guess corresponds to the values of the unknowns at the previous time step . If the Newton algorithm does not converge, the time step is restarted with a value divided by two .
5.2 Complementarity condition formulation
One of the advantages of this formulation is that the set of discret unknowns is fixed and it is not depending of the phase state in each finite volume K. It is given by
After discretization, we obtain at each time step n = 1, 2,…, N a nonlinear system denoted by
(15)where is the set of discret conservation equations with closure equations and is the set of discret complementarity equations in each cell . The only drawback of the introduction of the minimum function is that the problem (15) is no longer C ^{1} but we can solve it using the semismooth Newton’s method, called the Newtonmin method [9, 11]. At each Newtonmin iteration (l), we have to solve the following linear system
(16)where and denotes the generalized Jacobian of ϕ _{ K } at a point . This method has locally a superlinear convergence. Globally, when the starting point is not close enough to a solution, the algorithm may diverge, even when the functions are linear, as shown and discussed in [16–18]. The size of the linear system may also be reduced by using a similar splitting strategy as in the Coats formulation for primary and secondary variables but this has not been implemented yet in our prototype. For the stopping criteria and the time step control, we use the same strategy as the one presented for NVF.
The advantage of this formulation is to lead to a fixed set of unknowns and equations, which is easier to deal with the phase appearance and disappearance than the Coats formulation. Nevertheless, we have to resolve a local nonlinear system (10) for oil and gas phases and in each finite volume K to obtain the component molar fractions and propagate its derivatives in the density and viscosity laws by using the chain rule lemma. When the fugacity coefficients do not depend on concentrations, typically with K _{ i }(T, P) equilibrium models, it only occurs a linear scaling of phase concentrations which is easy to follow. In our case, where the fugacity coefficients are computed by an Equation of State (EoS), it results in a nonlinear scaling and artificial coupling between the concentrations which is more challenging. This is the price to pay for the elimination of fugacity equality constraints from the system.
As regards the initial conditions, we have to set the values of pressure, saturation, and composition of each phase at the initial state. In addition, for the complementarity formulation, we have to compute the initial fugacity starting from the reference pressure and compositions of the phases. If in a cell , the gas phase is present and the oil phase is absent, we have
On the other hand, if the oil phase is present and the gas phase is absent, we have
At last, if both oil and gas phases are present we can compute the fugacities using indifferently (17) or (18) because the two phases are supposed to be at thermodynamical equilibrium and then
It is important to remark that also with this formulation, the case of aquifer cells needs a particular treatment (see Sect. 6).
6 Aquifer cells
The case of cells where only the water phase and only the water component are present requires major consideration because the water phase is immiscible with other phases. In fact, in this case, for both the formulations, using the equations introduced in the previous sections the system diverges. Therefore we have to change the system of equations to solve. We define a new set of independent unknowns which are pressure and number of moles per unit volume of each component denoted by n _{ i } for all i = 1,…, n _{ c } and introduced in [23]. In these cells we still consider the conservation equation for the water component, however, the conservation equations for the hydrocarbon components become
where q _{ i } δ represents the source term. Therefore in the case of aquifer cells, for both the formulations, the system we have to solve is the following. If A _{ p } = {w}, find such that
In the following tables, we can find the number of unknowns and of equations of this particular case.
Equations  n _{ c } + 1 
Unknowns  n _{ c } + 1 
For the aquifer cells, if at least one n _{ i } becomes greater than zero, we must test the appearance of the oil phase and/or the gas phase as follows. We first compute the total molar fraction of each component
Knowing the overall composition and the pressure, we can perform a negative flash calculation to obtain the molar fraction of the gas phase V.

if V ≥ 1 the gas phase appears and then A _{ p } becomes A _{ p } = {w, g},

if V ≤ 0 the oil phase appears and then A _{ p } becomes A _{ p } = {w, o}.

if 0 < V < 1 the oil and gas phases appear and then A _{ p } becomes A _{ p } = {w, o, g}.
7 Numerical results
In this section, we compare the behavior of the natural variable formulation and the complementarity condition formulation by means of selected test settings. In order to model the different physical phenomena present in the study of flows in porous media, more choices are possible. The models chosen in our cases are the following.
7.1 Triphasic fluid flow model
The relative permeabilities are modeled using the Brooks and Corey’s correlation:

S _{ w }: water saturation,

S _{wi}: irreducible water saturation,

S _{orw}: residual saturation of oil in the oil–water flow,

S _{org}: residual saturation of oil in the oil–gas flow,

S _{ w }: gas saturation,

S _{gc}: saturation of critical gas,

: maximum values of the relative permeabilities,

nw, now, ng, nog parameters related to the size of the pores.
In our simulation we choose for the parameters the values presented in Table 1, Table 2 and Table 3.
Residual saturation.
End point of relative permeability.
Exponent of relative permeability.
The relative permeability of the oil is computed from the previous functions and from the Stone II model [24]
The other physical properties of oil and gas such as the fugacities and the densities are computed using the Peng–Robinson equation of state and for the computation of the viscosities the LohrenzBrayClark model [25] is used. The properties of water (density and viscosity) are computed using data from [26].
7.2 Peng–Robinson equation of state
For each phase α ∈ {o, g}, we write the Peng–Robinson equation of state [27] as follows

P: the reference pressure,

R: universal gas constant,

T: temperature,

v _{ α }: molar volume of phase α,

a _{ α }: Van der Waals potential of phase α,

b _{ α }: covolume parameter of phase α.
In a mixture, the factors a _{ α } and b _{ α } are defined by
where is the molar fraction of component i in phase α,
and δ _{ ij } is a binary interaction parameter between components i and j, and a _{ i } and b _{ i } are empirical factors for pure component i. The interaction parameters account for molecular interactions between two unlike molecules. By definition, δ _{ ij } is zero when i and j represent the same component, small when i and j represent components that do not differ much (e.g., when components i and j are both alkanes), and large when i and j represent components that are substantially different. Ideally, δ _{ ij } depends on pressure, temperature and on the identities of components i and j [28]. The factors a _{ i } and b _{ i } can be computed from
where T _{ci} and P _{ci} are the critical temperature and the critical pressure of component i, the EOS parameters Ω_{ ia } and Ω_{ ib } are given by
and ω _{ i } is the acentric factor of component i. Let
Introducing the compressibility factor
equation (21) can be expressed as a cubic equation in Z _{ α }:
This equation has three roots. When only one root is real, it is selected. If there are three real roots, the selection of the right one depends on the dominance of the liquid phase or the vapor phase. Now, for i = 1, …, n _{ c } and α = o, g, the fugacity coefficient of component i in a mixture can be obtained from
At last, the fugacity of component i is
7.3 Test of CO_{2} injection in a threecomponent system
The first case is a miscible gas (CO_{2}) injection in a quarter of fivespot saturated with oil. The domain has a size of 100 m in both x and ydirection and it is discretized using a 20 × 20 regular grid blocks. The reservoir model is homogeneous: the permeability is equal to 500 mD and the porosity is 0.3. The gas, composed only of CO_{2}, is injected with a constant rate that is equal to 80 m^{3}/day and the pressure at the producer is fixed to 55 bar. The temperature is assumed to be constant at 80 °C and the initial pressure is equal to 55 bar.
The total simulation time is 300 days, the initial time step is 0.05 day and the minimum and maximum time step are respectively 10^{−5} day and 20 days. The initial water saturation is given by S _{ w } = S _{wi} = 0.25 (S _{wi} is the irreducible water saturation) and the oil saturation is equal to 1 − S _{wi}. The oil and gas phases are a mixture of three components C _{1}, C _{6} and CO_{2} and the initial oil composition is given by C _{1} (20%), C _{6} (80%) and CO_{2} (0%).
Figures 1a and 1b show the distribution of the gas saturation in the domain for NVF formulation (Fig. 1a) and CCF formulation (Fig. 1b) at 300 days. Figures 1c and 1d show the distribution of the CO_{2} molar fraction in the gas phase for NVF formulation (Fig. 1c) and CCF formulation at 300 days (Fig. 1d). In Figures 1e and 1f, we report respectively the cumulative oil and gas production. The results are similar for the two formulations.
Fig. 1 Gas saturation and CO_{2} molar fraction in gas phase after 300 days for Natural Variables Formulation (NVF) (a) and (c) and for Complementarity Condition Formulation (CCF) (b) and (d). Oil and gas cumulative production (e) and (f). 
In Figures 2a–2d, we can see the oil and gas saturation and composition in the middle cell of the domain. The results are the same except for the gas composition when it is missing in the cell. Indeed, for the NVF, the composition of the absent phase is given by the negative flash solution and the sum of the molar fraction in this phase is equal to one, whereas for the complementarity condition formulation, it is obtained by solving the local nonlinear system (10) and the sum is less than one. It is also important to note that the gas phase appears exactly at the same time for both formulations (Fig. 2b).
Fig. 2 Oil and gas saturation and composition versus time in the middle cell. 
Table 4 summarizes the numerical results in terms of number of time steps, number of Newton iterations and number of restarted time steps. For this simulation, the NVF and the complementarity condition formulation performed with the same behavior.
Numerical results for both formulations NVF and CCF.
7.4 Test of CO_{2} injection in a sevencomponent system
The second case study still simulates a CO_{2} injection in a threedimensional quarter five spot saturated with oil. The reservoir size is 100 × 100 × 20 m and we use a 20 × 20 × 4 grid blocks to discretize the reservoir model. The fluid is a seven components model (C_{1}N_{2}, C_{23}, CO_{2}, C_{46}, C_{712}, C_{1319} and ), with the following initial composition: C_{1}N_{2} (38.8209%), C_{23} (14.5821%), CO_{2} (2.2685%), C_{46} (11.9334%), C_{712} (19.4598%), C_{1319} (8.7079%) and (4.2274%). The initial pressure and temperature are respectively 200 bar and 132.77 °C (above the bubble point). The CO_{2} is injected with a fixed rate of 300 m^{3}/day and the production pressure is 150 bar (below the bubble point). The other model’s properties are the same as in the first test.
The results given in Figures 3a–3f show that the agreement between the two formulations is quite good. Nevertheless, we observe slight differences for the oil and gas cumulative production.
Fig. 3 Gas saturation and CO_{2} molar fraction in gas phase after 300 days for Natural Variables Formulation (NVF) (a) and (c) and for Complementarity Condition Formulation (CCF) (b) and (d). Oil and gas cumulative production (e) and (f). 
In Figures 4a–4d, we plot the oil and gas saturation and composition (only for the first four components) in the cell at the middle and at the top of the reservoir. Once again, we observe that the results are identical for both formulations except for the gas composition at the very beginning of the simulation when this phase is absent in the cell.
Fig. 4 Oil and gas saturation and composition versus time in the top middle cell. 
Table 5 indicates that the NVF ran with a lower number of Newton’s iterations than the complementarity condition formulation.
Numerical results for both formulations NVF and CCF.
7.5 Test of condensate gas problem
The third case study is a condensate gas reservoir with a unique producer well in the center. The geometry of the reservoir is the same than in the first test and the mesh used is composed of 21 × 21 grid blocks. The initial gas saturation is 1 − S _{wi} with S _{wi} equals to 0.25. The fluid model is the seven components mixture used in the second test but with a different initial global composition: C_{1}N_{2} (66.8503%), C_{23} (14.6099%), CO_{2} (3.075%), C_{46} (6.8369%), C_{712} (6.0289%), C_{1319} (1.814%) and (0.785%).
We study two different initial points: the first one corresponds to an initial pressure and temperature equal to 290 bar and 196.85 °C (470 K), and the second point is near the critical point with an initial pressure and temperature equal to 310 bar and 132.78 °C (405.93 K) (Fig. 5). The production pressure is 240 bar (below the dew point).
Fig. 5 Phase diagram of the fluid. 
The total simulation time is 1000 days, the initial time step is 10^{−4} day and the minimum and maximum time step are respectively 10^{−6} day and 20 days.
The results are the same for both formulations and we report in Tables 6 and 7 only the numerical behavior of these formulations (the total number of time steps and the number of Newton’s iterations). In this test case, the complementarity formulation has better behavior than the NVF and especially when starting near the critical point for which the complementarity formulation requires fewer time steps and less Newton’s iterations than the NVF.
Numerical results for both formulations NVF and CCF and for the first initial point (290 bar–470 K).
Numerical results for both formulations NVF and CCF and for the second initial point (310 bar–405.93 K).
7.6 Test of gas injection in aquifer cells
For the last test, we model an injection of gas in an aquifer. This is a wellknown case in the reservoir industry and it is interesting because we suppose that the water phase is pure and immiscible with the other phases. We consider a twodimensional domain initially saturated only by water (S _{ w } = 1). The reservoir model and the fluid model are the same as in the first test.
Gas is injected, with a constant rate 30 m^{3}/day, in the aquifer through an injection well located at the lower left corner of the domain. The injected gas is composed by C_{1} (20%), C_{6} (5%) and CO_{2} (75%). A production well is located in the upper right corner of the domain and it works with an imposed pressure equal to 50 bar. The temperature is set at 80 °C and the initial pressure is 50 bar.
Once again, the two formulations lead to similar physical results. Regarding the numerical behavior, as we can see in Table 8 they perform the simulation with the same number of time steps and without restart. However, we can notice that the CCF formulation effectuates slightly more Newton’s iterations than the NVF formulation.
Numerical results for both formulations NVF and CCF.
8 Summary and conclusions
In this work, we develop two isothermal compositional formulations in the same framework and compare them for four different synthetic problems. The first formulation is the well known NVF which is commonly used in reservoir simulator, and the second one is the more recent complementarity condition formulation. The main advantages of this second formulation are to lead to a fixed set of unknowns and equations whatever the present phases and it allows to avoid the flash calculation. Nevertheless, the use of pressures, saturations, and fugacities as main unknowns induces more nonlinearities in the system, and we have to solve two local nonlinear systems in each cell and at each Newton’s iteration to obtain the component molar fractions in the oil and gas phases in function of the main unknowns.
Our first numerical tests indicate that the two formulations lead to similar physical results with quite similar behavior except for the condensate gas test cases and near the critical point for which the complementarity formulation outperform the NVF in terms of the number of time step and the number of Newton iterations.
Then in order to compare fairly the CPU times, we have to reduce the size of the linearized system at each Newton iteration by adding a splitting strategy (eliminating the closure equations and the complementarity constraints as it is done in the NVF formulation). But the choice of the primary unknowns and remaining secondary unknowns depends on the context in each cell. This work is under investigation.
A study of the semiimplicit time scheme for the complementarity condition formulation is a work in progress. This scheme consists in calculating the flux terms implicitly with respect to the pressure and the saturations, and explicitly with respect to the fugacities. More complex and realistic test cases will also be treated to compare the two formulations.
References
 Aziz K., Settari A. (1979) Petroleum reservoir simulation, Vol. 476, Applied Science Publishers, London, UK. [Google Scholar]
 Chavent G., Jaffre J. (1986) Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media, Elsevier Science. [Google Scholar]
 Danesh A. (1998) PVT and phase behavior of petroleum reservoir fluids, Vol. 47, Elsevier. [Google Scholar]
 Firoozabadi A. (1999) Thermodynamics of hydrocarbon reservoirs, McGrawHill, New York, NY. [Google Scholar]
 Coats K.H. (1989) Implicit compositional simulation of singleporosity and dualporosity reservoirs, Paper SPE 18427, 6–8. [Google Scholar]
 Abadpour A., Panfilov M. (2009) Method of negative saturations for modeling twophase compositional flow with oversaturated zones, Transp. Porous Media. 79, 2, 197–214. [Google Scholar]
 Angelini O., Chavant C., Chenier E., Eymard R., Granet S. (2011) Finite volume approximation of a diffusiondissolution model and application to nuclear waste storage, Math. Comput. Simul. 81, 10, 2001–2017. [Google Scholar]
 Masson R., Trenty L., Zhang Y. (2014) Formulations of two phase liquid gas compositional Darcy flows with phase transitions, Int. J. Finite Vol. 11, 1–34. [Google Scholar]
 Gharbia I.B., Jaffre J. (2014) Gas phase appearance and disappearance as a problem with complementarity constraints, Math. Comput. Simul. 99, 28–36. [Google Scholar]
 Coats K.H. (1980) An equation of state compositional model, Soc. Petrol. Eng. J. 20, 5, 363–376. [CrossRef] [Google Scholar]
 Qi L., Sun J. (1993) A nonsmooth version of Newton’s method, Math. Progr. 58, 1–3, 353–367. [CrossRef] [Google Scholar]
 Facchinei F., Pang J.S. (2003) Finitedimensional variational inequalities and complementarity problems (two volumes), Springer Series in Operations Research, Springer. [Google Scholar]
 Hager C., Wohlmuth B. (2010) Semismooth Newton methods for variational problems with inequality constraints, GAMMMitteilungen 33, 1, 8–24. [CrossRef] [Google Scholar]
 Whitson C.H., Michelsen M.L. (1989) The negative flash, Fluid Phase Equilib. 53, 51–71. [Google Scholar]
 Lauser A., Hager C., Helmig R., Wohlmuth B. (2011) A new approach for phase transitions in miscible multiphase flow in porous media, Adv. Water Resour. 34, 8, 957–966. [Google Scholar]
 Ben Gharbia I., Gilbert J.Ch. (2012) Nonconvergence of the plain Newtonmin algorithm for linear complementarity problems with a Pmatrix, Math. Progr. 134, 2, 349–364. [CrossRef] [Google Scholar]
 Ben Gharbia I., Gilbert J.Ch. (2013) An algorithmic characterization of Pmatricity, SIAM J. Matrix Anal. Appl. 34, 3, 904–916. [Google Scholar]
 Ben Gharbia I., Gilbert J.Ch. (2018) An algorithmic characterization of Pmatricity II: adjustments, refinements, and validation, SIAM J. Matrix Anal. Appl. (in revision). [Google Scholar]
 Michelsen L.M. (1986) Simplified flash calculations for Cubic Equations of State, Ind. Eng. Chem. Process Des. Dev. 25, 1, 184–188. [CrossRef] [Google Scholar]
 Sun A.C., Seider W.D. (1995) Homotopycontinuation method for stability analysis in the global minimization of the Gibbs free energy, Fluid Phase Equilib. 103, 4, 213–249. [Google Scholar]
 Hintermüller M., Ito K., Kunisch K. (2003) The primaldual active set strategy as a semismooth Newton method, SIAM J. Optim. 13, 865–888. [Google Scholar]
 Qi L. (1993) Convergence analysis of some algorithms for solving nonsmooth equations, Math. Oper. Res. 18, 18, 227–244. [CrossRef] [Google Scholar]
 Eymard R., Guichard C., Herbin R., Masson R. (2012) Vertexcentred discretization of multiphase compositional Darcy flows on general meshes, Comput. Geosci. 16, 987–1005. [Google Scholar]
 Stone H.L. (1973) Estimation of threephase relative permeability and residual oil data, J. Can. Petrol. Technol. 12, 4. [CrossRef] [Google Scholar]
 Lohrenz J., Bray B.G., Clark C.R. (1964) Calculating viscosities of reservoir fluids from their compositions, J. Pet. Technol. 16, 10, 1–171. [CrossRef] [Google Scholar]
 Schmidt E. (1969) Properties of water and steam in SIunits, SpringerVerlag, Berlin, Germany. [Google Scholar]
 Peng D.Y., Robinson D.B. (1976) A new twoconstant equation of state, Ind. Eng. Chem. Fundam. 15, 1, 59–64. [CrossRef] [Google Scholar]
 Chen Z., Huan G., Ma Y. (2006) Computational methods for multiphase flows in porous media, Computational Science and Engineering, SIAM, Philadelphia, PA. [CrossRef] [Google Scholar]
All Tables
Numerical results for both formulations NVF and CCF and for the first initial point (290 bar–470 K).
Numerical results for both formulations NVF and CCF and for the second initial point (310 bar–405.93 K).
All Figures
Fig. 1 Gas saturation and CO_{2} molar fraction in gas phase after 300 days for Natural Variables Formulation (NVF) (a) and (c) and for Complementarity Condition Formulation (CCF) (b) and (d). Oil and gas cumulative production (e) and (f). 

In the text 
Fig. 2 Oil and gas saturation and composition versus time in the middle cell. 

In the text 
Fig. 3 Gas saturation and CO_{2} molar fraction in gas phase after 300 days for Natural Variables Formulation (NVF) (a) and (c) and for Complementarity Condition Formulation (CCF) (b) and (d). Oil and gas cumulative production (e) and (f). 

In the text 
Fig. 4 Oil and gas saturation and composition versus time in the top middle cell. 

In the text 
Fig. 5 Phase diagram of the fluid. 

In the text 