© I. Ben Gharbia & E. Flauraud, published by IFP Energies nouvelles, 2019

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 time-space 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 Newton-min 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 three-phase 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₂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

(2)

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

where

• K: permeability of the medium [m²],

• 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:

(3)

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

(5)

The chemical potential μ _i [J mol⁻¹] of a component i is given by

(6)where

• f _i : fugacity of the component i [Pa],

• P: reference pressure [Pa],

• T: temperature [K],

• R: universal gas constant [J mol⁻¹ K⁻¹],

• : chemical potential of pure component i in the state of perfect gas at the same pressure P and temperature T [J mol⁻¹].

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 well-defined. 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 space-time 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 time-space 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

(9)

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 so-called complementarity conditions. This approach uses the fact that these constraints can be reformulated equivalently as a non-differentiable 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

(10)

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

(11)

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

(12)

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

where α ∈ A.

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 and t ^N  = T and we denote by Δt ⁿ the time step such as Δt ⁿ  = t ⁿ −t ⁿ⁻¹ for n = 1,…, N. The systems (9) and (4) are discretized using fully implicit Euler integration in time and a cell-centered finite volume scheme with a two-point flux discretization. The mobility terms are up-winded 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 ⁿ by:

• if

• if

• if

We obtain at each time step n = 1, 2, …, N a nonlinear system denoted by

(13)

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.

(14)

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

and

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

for all .

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 ¹ but we can solve it using the semismooth Newton’s method, called the Newton-min method [9, 11]. At each Newton-min 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 non-linear 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

(17)

On the other hand, if the oil phase is present and the gas phase is absent, we have

(18)

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

(19)

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:

where:

• 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.

The relative permeability of the oil is computed from the previous functions and from the Stone II model [24]

(20)

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 Lohrenz-Bray-Clark 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

(21)where

• 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₂ injection in a three-component system

The first case is a miscible gas (CO₂) injection in a quarter of five-spot saturated with oil. The domain has a size of 100 m in both x and y-direction 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₂, is injected with a constant rate that is equal to 80 m³/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⁻⁵ 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 ₁, C ₆ and CO₂ and the initial oil composition is given by C ₁ (20%), C ₆ (80%) and CO₂ (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₂ 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.

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).

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.

7.4 Test of CO₂ injection in a seven-component system

The second case study still simulates a CO₂ injection in a three-dimensional 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₁N₂, C₂₃, CO₂, C₄₆, C₇₁₂, C₁₃₁₉ and ), with the following initial composition: C₁N₂ (38.8209%), C₂₃ (14.5821%), CO₂ (2.2685%), C₄₆ (11.9334%), C₇₁₂ (19.4598%), C₁₃₁₉ (8.7079%) and (4.2274%). The initial pressure and temperature are respectively 200 bar and 132.77 °C (above the bubble point). The CO₂ is injected with a fixed rate of 300 m³/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.

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.

Table 5 indicates that the NVF ran with a lower number of Newton’s iterations than the complementarity condition formulation.

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₁N₂ (66.8503%), C₂₃ (14.6099%), CO₂ (3.075%), C₄₆ (6.8369%), C₇₁₂ (6.0289%), C₁₃₁₉ (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).

The total simulation time is 1000 days, the initial time step is 10⁻⁴ day and the minimum and maximum time step are respectively 10⁻⁶ 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.

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 well-known 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 two-dimensional 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³/day, in the aquifer through an injection well located at the lower left corner of the domain. The injected gas is composed by C₁ (20%), C₆ (5%) and CO₂ (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.

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 non-linearities 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 semi-implicit 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

All Tables

All Figures