Recuperation assistée d'hydrocarbures en configuration cavernes/fractures



Belerrajoul Mohamed: 3rd year student at ENSEEIHT INP Toulouse, fluid mechanics dept

Kong Gaopan: 2nd year International master at ENSEEIHT INP Toulouse, fluid mechanics dept

Rocher Léa: 3rd year student at ENSEEIHT INP Toulouse, fluid mechanics dept


Supervised by :


Bernard MONTARON: Schlumberger China Petroleum Institute






The BEI (Bureau d'Etudes Industrielles) is a six weeks project which allows last year engineering students to work as a team on a industrial challenge. The project which is presented here has been developped in collaboration with Schlumberger China Petroleum Institute.

The student team has joined the Simkarts Project, initiated in 2012 by Pr. Han Pingchou, Peking University PKU, and Bernard Montaron, Schlumberger China. This project has been created in order to develop a reservoir simulator which would simulate realistic multi-phase flows in the particular geologic context of the Tarim Basin, situated in West China. However, this simulator is not based on fluid dynamics laws and therefore needs to be validated.

In the BEI, three-phase flow (oil/water/gas) simulations are carried out, using  Jadim and Fluent, two CFD softwares, in order to simulate nitrogen injection in the specific fracture/cave configuration of the Tarim Basin.  The simulation method used is the Volum Of Fluid (VOF) method ,based on Navier-Stokes equations.

These simulations will give a good understanding of flow behavior in a case of gas injection and will help to optimized this oil recovery method. Moreover, results will permit to validate and optimized the reservoir simulator developped by Schlumberger.

The team




​    3rd year ENSEEIHT Student, Fluid Mechanics dpt                            



       Léa ROCHER

       3rd year ENSEEIHT student, Fluid Mechanics dpt



    Gaopan KONG

    International master student, Fluid Mechanics dpt



The supervising team

Mr Climent: ENSEEIHT professor and director of IMFT (Institut de Mécanique des Fluides de Toulouse)

Mme Pedrono: member of the IMFT (Institut de Mécanique des Fluides de Toulouse)


The industrial partner          

Schlumberger is the world’s leading supplier of technology, integrated project management and information solutions to customers working in the oil and gas industry worldwide. The company provides the industry’s widest range of products and services from exploration through production.

Present in China since 1980, Schlumberger is actually involved in project which take place in the Tarim Basin. Bernard Montaron has initiated the Simkarts Project in 2012 and will be our interlocutor inside Schlumberger China.


Industrial context

The Simkarts project has been initiated in 2012 by Pr. Han Pingchou, PKU, and Bernard Montaron, Schlumberger China, in order to understand  multi-phase flow dynamics in deep oil reservoirs located in the Tarim basin, West China.

Karst reservoirs of Tarim Basin in Xinjiang Province of People’s Republic of China have been determined as an area with large quantities of oil and gas and are exploited by Sinopec and Petrochina Tarim Oil Company​ since the 1990’s.

This carbonate formation of Ordovician age is buried at more than 6000 meters and represents a complex network of connected caves full of oil and water large of 20 x 20 km.                                                           


                           The Tarim Basin                                                        A fracture/cave reservoir  

In order to optimize oil production, it is necessary to understand and predict oil flow inside the reservoir. Indeed, the fluid dynamic inside the reservoir is complex.  In the Tarim basin, an aquifer is situated below the reservoir. Therefore, as oil is pumped up, water goes up by capillarity effect. Moreover, in many reservoir, water or gas are injected in order to maintain pressure inside the reservoir.

If geologists can find oil deep in the ground by using variety of methods such as sismology, surface features, soil types, sniffers or satellites images, numerical simulators are needed in order to understand fluid dynamic inside reservoirs.

Very few oil reservoirs in the world produce oil from caves. In traditional reservoirs oil is contained in microporous rocks with typical pore size smaller than 30 microns. 

Existing reservoir simulators currently used by the oil and gas industry are designed to simulate fluid flow in porous media governed by Darcy's law.

These simulators are not appropriate for simulating fluid flow in cave networks where Darcy's law does not apply. Moreover, Navier-Stokes model is impracticable because of the complexity of the cave's geometry, the scale of the reservoir and massive calculation and data burden that it would involve.

Thus, the ultimate objective is to design a 3D reservoir simulator that can make realistic multi-phase flow simulations (gas, water, and oil) in large and complex cave networks made of thousands of caves distributed in 3D space and connected according to a given pattern. 


Bernard Montaron, Schlumberger China, proposed to investigate the feasibility of using cellular automata technology for such a simulator, and therefore initiated the SimKARST Project in 2012 with Pr. Han Pingchou, Peking University PKU.

The first part of the project started in October 2012 with two PKU students who developed a prototype cellular automata software to simulate 2D flow of water and oil in a cell with a simple geometry.

A cellular automaton is a discrete model where particules are constrained to move on points of the lattice. Each particle is in a certain state, which is be in this case oil, water, gas or rock. Movement of particules between neighboring lattices points are governed by prescribed rules based on gravity for the vertical state and on probablility for the horizontal state. For instance, a rule stipulate that oil blocks will always move upwards when water blocks are on top of them.

If this method is a fast simlulation technique and give reasonable kinetic results, it is not based on any fluid dynamics and does not respect the physical properties of fluids. Moreover, velocities are limited by lattice spacing and time. Therefore, this method needs to be compared with experiments and simulations based on physics laws, in order to validate its accurency.


                                 Representation of the reservoir by the automata simulator

It order to validate and to improve the automata simulator, two ENSEEIHT students, Nicolas Sobecki and Thibault Moreau, realized a summer intership in China  in order to conduct two-phase flow experiments on a specific geometry. The acquisition of this experimental data allowed to compare the experimental flow behaviour observed with  flow simulated with the cellular automata software.

At the same time, Tinting Han, another ENSEEIHT student, realized CFD studies with fluent in order to simulate  two-phase flow on the same geometry, in order to compare the results with experiments and data from the automata cellular. As CFD simulations are based on the Navier and Stokes equations, they respect physics law and give therefore accurate results.

          Two-phase flow experiment                                                     CFD Simulation with Fluent
As the simulator is still under development, a ENSEEIHT students team had, during their BEI project, carried out CFD simulations using complex geometries. These simulations, realized with two different CFD softwares and based on Navier and Stokes equations, helped Schlumberger to validate and to improve its simulator.

Gas Injection

Schlumberger came to us with another phenomenon to simualtion: the injection of gas in order to optimize oil recovery.

Indeed, oil extraction can be devided into different recovery stage. During the primary recovery stage, which corresponds to the first years of production, reservoir drive comes from a number of natural mechanisms: underground pressure in the oil reservoir is sufficient to force oil to the surface. However, over the lifetime of the well, pressure fall and is insufficient to force oil to the surface.

Therefore, after natural reservoir drive diminishes, secondary recovery methods are applied in order to increase reservoir pressure and thereby stimulate production.

Water or gas injection is a secondary recovery method used by oil and gas companies to increase oil recovery from an existing reservoir. Water or gas is injected to support pressure of the reservoir (also known as voidage replacement), and to sweep or displace oil from the reservoir, and push it towards a well.

Normally only 30% of the oil in a reservoir can be extracted, but this method increases that percentage (known as the recovery factor) and maintains the production rate of a reservoir over a longer period.

In the Tarim Basin, due to the complex geometry of the karst reservoir, a large amount of oil is wedged into the cap of the cavities:  38,000 m3 of oil produced (32%) for 80,000 m3 remaining (68%).


                                          Oil wedged into the cap of the cavities

The solution proposed by Schlumberger is to inject nitrogen, lighter than oil, in the reservoir to replace oil wedged  in the cavities. 

In this BEI, CFD simulations are therefore carried out in order to simulate this phenomenon. These CFD simulations will than help Sclumberger to understand flow behaviour in a case of gas injection and see how the technique can be optimized. Moreover, a comparison of CFD simulation and automata simulator results will permit to validate the reservoir simulator developped by Schlumberger for a three-phase flow.


The objectives of the BEI is to:
1. Understand  three-phase flow dynamics when nitrogen is injected in a fracture/cave configured reservoir
  •  A bibliography of different methods to simulate multiphase flow is realized with a focus on the VOF method.
  • The main caracteristics of Jadim and Fluent, two CFD codes which permit three-phase flow simulation, are analyzed.
  •  CFD simulations of nitrogen injection in a simple cave geometry (three-phase flow) is carried out with Fluent and Jadim.
  • A parametic study is also carried out in order to observe the influence of the gas injection well position and of the outlet position.

2. Compare a commercial software (Fluent) and a research software (Jadim)

3. Compare CFD simulations based on physics Navier and Stokes equation with cellular automata simulations in order to validate the cellular automata simulator


Multiphase-flow simulations

The numerical simulation of flows with interfaces is a vast and complex topic, with applications in domains as varied as environment, geophysics, engineering, and fundamental physics.

Indeed, phenomena considered often happen on scales of space and time where experimental visualizations are difficult or impossible. In such cases, numerical simulation may be a useful tool to validate the intuition of the physicist, the engineer, or the mathematician.

In our case, gas injections in the fracture/cave configuration has not been developped yet. Therefore, the simulations carried out in the BEI will help our industrial partner  to anticipate and to optimize futur gas injections in the reservoir.

Different numerical techniques exit to simulate three-phase flows. If most of them only allow low Reynolds or Weber numbers, the continuous improvement of computational power extends the range of affordable problems.

In the first part of the project, an analyze of different simulation methods based on Navier-Stokes equations is conducted. Understanding problems and limitations of each methods will allow us to choose the most appropriate method to our case and to be aware of its limitations.

Moreover, a comparative study of Jadim and Fluent, two CFD software, is conducted to know how they handle three-phase flow simulations.


Methods for direct numerical simulation of multiphase flow

If numerical simulations of three-phase flow are still a challenge in the CFD field and the purpose of lots of research groups, differents methods already exist.

Three main classes of methods for direct numerical simulations of multiphase flow exists:

The Multiphase particle-in-cell method (mixed Eulerian-Lagrangian method)

lt is an hybrid method where individual particles in a Lagrangian frame are tracked in continuous phase space. This method solves Navier-Stokes equations (Eulerian) for the continuous phase and Lagrangian equation of motion for individual particles.
Interactions between particles are calculated on the Eulerian grid and convection algorithms are highly accurate (Lagrangian advection is non-diffusive)
However, this method has limitations:
  • All particles are assumed to be spherical. If corrections for non-spherical particles can be included in particle drag model, true interactions may therefore not be well represented.
  • The size of particles must be small compared to the Eulerian grid  for accurate interpolation.​
  • The major problem of this method is that only a low volume fraction of the dispersed second phase is acceptable. This last consideration makes the model inappropriate for the modeling of liquid-liquid mixtures, fluidized beds, or any application where the volume fraction of the second phase is not negligible.

Moving grids methods (Lagragian methods)​

The interface is a boundary between two subdomains of the grid. ​The method is mainly used to follow the motion of a rising bubble, small amplitude waves or weakly deformed bubble. However, when the  points move simply in a Lagrangian way, the grid may deform considerably and in case of large deformations of the interface or topology changes, the grid has to be remeshed. This complex method is only sucessful and very accurate for small interface deformation, which is not the case here.


Fixed-grid methods (Eulerian methods)

These methods correspond to a predifined grid that does not move with the interface.
  • The marker method
Tracers or marker particles moved in a Lagrangian way are used to locate the interface on a fixed grid.

The advantage of the method is that it allows the capture of details of interface motion on scales much smaller than the grid spacing  (Sub grid resolution) and a high degree of accuracy that may be achieved by representing the interface through high-order interpolation polynomials.

However, if surface markers are more accurate than volume markers, they cannot handle more than two phases and topology changes are difficult. Control the marker distribution is also problematic and it needs an important CPU time.
  • The Level Set method
It is a numerical technique for tracking interfaces that can perform numerical computations on a fixed  grid without having to parameterize the interface. A new dimension is introduced to the case and defines the interface as a level set of function φ( x, y) which represents the minimum distance from each point to the interface.
This method tracks the evolution of φ and determines the zero level set thanks to the equation:
                                    $\frac{\partial φ}{\partial t}+U. \Delta φ=0$
This method is a great tool for modeling time-varying objects and works in any dimension. Moreover, no special treatment is needed for topological changes ans the interface geometry reconstruction is simple.
However, this method has not yet produced wide range of results, especially in 3D compare to other methods and does not guarantee conservation of the volume.
  • The VOF method (Volume of Fluid)
It is a surface-tracking technique applied to a fixed Eulerian mesh where Navier Stokes equations which describe the motion of the flow have to be solved separately.
The method is based on the solution of a transport equation for variable ‘C’ (often also referred as indicator or colour function) for the liquid phase.
$C_{ij}$ represents the portion of the area of the cell (i, j) filled with liquid phase  and the phase function χ :
where  0 < C < 1 in cells cut by the interface S and C = 0 or 1 away from it.
The VOF method doesn't explicitly track the interface, it reconstructs the interface based on calculation of the volume fraction of fluid . The Color Function also cannot be solved easily. Several methods for the reconstruction of the interface exist, the most popular being PLIC (Piecewise Linear Interface Calculate). In a 3D space, the interface can be described by


where n is the normal vector to the interface and $\alpha$ is a constant line. 

$\alpha$  can be solved by root-finding method or analytical formulas $\alpha=\alpha(C)$ and  $\vec{n}$ has several approaches such as  Parker and Yong's method, Least-squarts method etc.

  • Advantages of the VOF method

The VOF's use, reliability and effectiveness are widespread: the method has been known for several decades, has gone through a continuous process of improvement and is used by many commercially available software programs.

Moreover, the volume conservation is good and no special provision is necessary to perform reconnection or breakup of the interface  as the change of topology is implicit in the algorithm. The VOF method is easy to extend to 3D of space and simple to implement.

Applications of the VOF model include stratified flows or the steady or transient tracking of any liquid-gas interface, which correspond to our case.

  • Limitations of the VOF method

The simplification of nonlinear terms and the fact that high order terms are omitted after discretization can lead to less accurate solutions in some cases. Moreover, C<0 or C>1 is possible, which is not relevant with the physic. Another limit of the method is the transitional region between the phases which has to be at least equal to the grid distance.                                                                                                                    

If the interface geometry reconstruction is challenging, normal interface movements are not straightforwards. Last but not the least, the VOF method involves massive calculations and data burden, which leads to important CPU times.





CFD Softwares: Fluent versus Jadim

Two different software are used to simulate the gas injection: FLUENT and Jadim. 

Fluent is a commercial CFD software which provide a wide range of multiphase flow models to model the behavior and interactions of a large number of phases. The user can not have access to the code; therefore, there are no possibilities to check equation implementation or to had function in the code.This software has been largely used in the industry for decades.

JADIM LES is a research code developped by Jacques Magnaudet and Dominique Legendre's team in the Interface group at IMFT(Institut de Mécanique des Fluides de Toulouse). The code permits to describe in an accurate way such physical mecanisms present
in multi-phase fows. It is especially  able to simulate two and three-phases flow with the VOF method without reconstruction of the interface. Jadim resolves Navier-Stokes equations in 3D for incompressible fluids and instationnarities.
A second order space and time finite volum method with a structured mesh is used (third order Runge-Kutta scheme for non-linear term resolution coupled with a Crank-Nicolson scheme for the semi-implicit part). This code can only be used with structured mesh.

The table below compare  the main caracteristics of Jadim and Fluent:


It is important to be aware about the difference that exits between the two codes in order to understand simulations issues such as lack of convergence or simulations results such as differences in flow features.


More informations about    can be found on this website:

More informations about can be found on this website:



The bibliography carried out leads us to conclusions developped in the table below:


Simulations will be carried out with the VOF method on both Jadim and Fluent codes and results will be compare. However, as the level set method is not available for three-phase flow on Jadim, this method will only be tested with Fluent.


The simulations

In this part of the project, three-phase flow simulations are carried out using a simple geometry which represents oil wedged into the cap of a cave. Simulated with Fluent and Jadim with the VOF method, the results are analyzed and compared in order to conclude on the accuracy of the recovery method.

It is important to note that these results are only relevant for this particular geometry and can not be applied to the geometry of the Tarim Basin. They will allow us to compare Jadim and Fluent and will then be used by Schlumberger to validate the automata software for the gas injection recovery method.


Presentation of the simulation

The geometry

As three-phase flow simulations are complex calculations for Jadim and Fluent, a simple geometry is used. The geometry represents one cave with oil wedged on the top and has same dimensions as the geometry used for experiments carried out in 2013.


The boundary conditions are:

- An inlet velocity which represents  the gas injection, situated on the bottom-right of the domain

- A pressure outlet which represents a connection to another cave, situated on the middle left side of the domain

- Walls otherwize which represent the rocks

The simulation starts with an initiale quantity of oil wedged which represents 37% of the domain.

A solid obstacle is defined in the middle of the rectangular geometry in order to avoid a direct oil flow into the outlet and to force gas flow to go up and therefore replace the oil.


In order to know if the flow is compressible or incompressible, it is important to look at the value of the Mach number. This dimentionaless number is defined by:

$  Ma= \frac {U }{ a} $   where a is the sound velocity (a=340 m/s) and  U is the velocity of the flow.

A flow can only be considered as incompressible if its Mach number is smaller than 1.

As velocities inlet used for the simulation do not exceed 1 m/s, velocities of the three-phase flow will be small enough to consider that a Mach number of the flow smaller than 1.

In the simulations, the flow will be considered as incompressible.

Physical Properties

CFD simulations are based on Navier and Stoles equations. Therefore, it is essential to know the physical properties of the differente phases as they have a direct impact on the flow behavior.

Investigate real properties of phases is difficult. Indeed, due to the large dimension of the Tarim basin, oil properties can change from one well to another. Moreover, fluids and gas properties change in function of the temperature, which is itself variable, especially is function of the depth. Regarding gas viscosity, and as gas injection has not been tested yet in the reservoir, it is not known yet if the gas properties will be modified or not before being used.

After  discussion with our supervising team, and on the advices of Mr Montaron, choices have been made to use:

- The same properties of oil and water used for the experiments carried out in summer 2013

- Gas properties of an unmodified nitrogen gas

- Determinate surface tension value between oil/air/water phases instead of oil/nitrogen water phases.

All the physical properties used for the simulations are recap in the table below:


Gravity is constant and equals to -9.81 m/s².


Multiphase flow simulation is a complex topic, still in development and which need long calculation times, lots of memory and high-performance material. Therefore, for this BEI, we choose small velocities in order to have a laminar flown which would not require to add turbulence models to the simulation.

Indeed, it is difficult to calculate the Reynolds Number of the flow because the velocity inside the domain is variable. In order to limit the turbulence phanomenon in the domain, the Reynolds Number of the gas is calculated, which must be small enough not to have a turbulent jet at the inlet.

                                              $$ Re_{gas}=  \frac{U_{inlet}\times D\times \rho} { \mu }$$

with $  \rho_{gas}=1.251 kg/m³$, $D=0.025m$ and ${ \mu_{gas}=12.10^{-6} kg/ms }$

A flow can be considered turbulent for a Re>2000. For each simulation, the Renolds Number of the gas is calculated in order to be sure that velocities used are not too high.

Moreover, for each simulation, it is important to verify that velocities in the domain are not to high.


Jadim Simulations

The Jadim code: simulation parameters

Jadim is a research code developped by the IMFT to simulate two and three-phases flow with the VOF method without reconstruction of the interface.

The code is in Fortran 90 and the post-treatment can be acheived with Tecplot or Paraview, open source softwares and data can be analyzed with Matlab. In the project, the post-treatment is acheived with Paraview and Matlab.

Jadim needs at least 8 input files to run. The *.para le is necessary to fix the numerical parameters, *.geom, *.bord and *.datr descript the geometry, *.limi gives the boundary conditions. The *.init file permits to defined the initial conditions, the *.phys concerns the physical properties and the last one jcl file gives the unit number of Jadim input and output les.

In order to create the *.geom, *.bord, *.datr  and the *.limi, a mesher is included witch only provide cartesian structured mesh. The geometry and the mesh are created together. The square geometry is first entirely meshed. Then a obstacle is added in order to represent the empty square in the middle of the geometry. In the obstacle, velocities and pressure fields are always equal to zero.

The inlet and the wall conditions are defined as Dirichet conditions where the velocity value is zero for a wall and different from zero for the inlet.


Simulations can be restarted from an old binary and can be run in parallele, using more than one processor to compute and therefore decreasing the CPU time.

Mesh impact study

 The mesh can have a major influence on the results of a simulation. Therefore, a mesh study is carried out in order to investigate the impact of the mesh in the three-phase flow case. As Jadim can only sustain structured mesh, a coarse structured mesh and a refined structured mesh are tested. Two simulations are compared, one using a coarse mesh and one using a refined mesh. The main parameters of the simulations are recap in the table below:

Type of mesh Number of cells  Velocity Inlet (m/s) Gas Reynold Number Time Step max Physical Time  
Refined Mesh     270x352            0.1                260         10-4     1.93 sec  
Coarse Mesh       70x92            0.1                260         10-4     1.93 sec  

As Jadim does not give direct access to the quality of the mesh, an advanced analyzis of the results is conduct in order to appreciate the impact of the mesh on the results.

In the two pictures below, the oil contour is plotted at the dimentional time 0.025. A first observation of these results shows that a refined mesh reduce the phenomenon of numerical diffusion. Indeed, the interface between oil and the other phases obtained with a refined mesh is clearer and more accurate than the one obtained with a coarse mesh. However, oil seems to mixed with other phases in both cases.

The critical analyze of the results must take into account the numerical diffusion phenomenon. Indeed, if a refined mesh reduce its effect, it does not eliminate it.


     Oil volume fraction - Coarse mesh                                Oil volume fraction- Refined mesh

In order to compare the two meshes, the phase volume fraction in the domain and the outlet oil flow rate are plotted and compare.

The analyze of these data is done in the part "analyze of the results". Here, only the impact of the mesh is studied.

The graph below shows the evolution of mean oil, gas and water volume fraction in the domain in function of an dimentionless time. In this case, the physical time of the simulation is 1.93 sec, which correponds to a dimentionless time of 0.036.

The comparison of the simulation shows that after one seconde, the oil, gas and water ratio in the domain is not the same for the both meshes. Moreover, as the gas injection is continuous, the evolution of gas should be linear in function of time. If a refined mesh gives a linear evolution of the gas volume fraction, a coarse mesh does not.

The importante numerical diffusion observed for the coarse mesh makes the phase mixed together and therefore, the percentage of each phase in the domain is less accurate.


The evolution of outlet oil flow rate in function of time is also an interesting graph to analyze. The value is divided by the input gas flow rate in order to work with dimentionless values.

The graph below shows that the maximum flow rate value obtained for a refined mesh is twice the value obtained for a coarse mesh. In the pictures representing the oil volume fraction at time 0.025, it can be observed that at this time, only oil goes out of the domain. Therefore, at time 0.025, the outlet oil flow rate should be equal to the input gas flow rate. The graph shows that for a refined mesh, the flow rate ratio equals to 1 at time 0.025, while a coarse mesh gives a ratio of 0.5.

Therefore, a coase mesh makes a error of 50% on the value of the outlet oil flow rate.


In conclusion, both local flow features and macroscopic data are affected by the mesh quality.  A refined mesh is necessary in order to have accurate results and to limit the numerical diffusion phenomenon.

Convergence issues

For the first simulation carried out, real properties of oil, gas and water were used with a velocity inlet of 1m/s. Quickly, convergence issue appeared: on one hand, as the maximum time step allowed was to high, the simulation stopped due to a too weak time step and on another hand, velocity and pressure values in few cells diverged.

Even after descreasig the maximum time step, convergence issue remained. Two solutions were tried to resolve the problem:

- to multiply oil, gas and water viscosities by 100

- to decrease the velocity to 0.1 m/s

This two cases did not show convergence problem and therefore have been carried out. Results of these simulations are analyzed in the part "analyze of the results".

CPU time

The complexity of resolving three-phase flow induce high CPU time. Moreover, as the mesh is refined and the time step is small, simulations lasts even longer.

Therefore, the choice of the simulations is important as the project last only six weeks.


Fluent Simulations

Simulation parameters

The VOF model proposed by Fluent can model two or more immiscible fluids by solving a single set of momentum equations and tracking the volume fraction of each of the fluids throughout the domain. The VOF formulation relies on the fact that two or more fluids (or phases) are not interpenetrating.

To reconstruct the interface, Fluent offers different schemes. In the geometric reconstruction and donor-acceptor schemes, FLUENT applies a special interpolation treatment to the cells that lie near the interface between two phases.

As small velocities inlet are used in order to avoid turbulence, a laminar model is selected.

Concerning the boundary conditions:

  • The outlet boundary condition is a pressure outlet initialized at 0 Pa.
  • The gas injection is represented by a velocity inlet boundary condition.

The initialization of the domain is an important step in the preparetion of the simulation. The domain is first initialized with water evrywhere. Then, a region is created and which represents the initial position of the oil (Adapt/ Region... Mark). Then, the region is associated to the oil (Solve/ Initialized/ Patch... Oil Volume Fraction=1).

The convergence criterions are fixed at $1\times 10^{-6}$ and the maximimum number of iteration per timestep is 30.

Only the pressure-based solver can be used.The VOF model proposed by Fluent is not available with either of   the density-based solvers. Therefore, a pressure-based solver is used.

Nota Bene:

The VOF model does not allow void regions where no fluid of any type is present. However, this restriction is not an issue in the case as the cave is initially filled up with oil and water.

Moreover, only one of the phases can be defined as a compressible ideal gas. There is no limitation on using compressible liquids using user-defined functions.

First order versus second order

With Fluent, the second-order time-stepping formulation can only be used with the VOF implicit scheme. As Jadim is only available for a second-order formulation, the second-order implicit formulation is used for the comparison between Fluent and Jadim.                                              

However, as explicit schemes are easier to program, with only simple calculations performed at each timestep, they required smaller timestep, and thus reduce the computational time of a simualtion.

A simulation is carried out with first a second order implicit scheme and then a first order explicit scheme and the results appeared to be similar. Thus, for the parametric study which require many heavy simulations, the choice has been made to use a first-order explicit formulation to reduce the computational time.

The implicit scheme is used for time discretization, Fluent's standard finite-difference interpolation schemes, QUICK, Second Order Upwind and First Order Upwind, and the Modified HRIC schemes, are used to obtain the face fluxes for all cells, including those near the interface.

In the explicit approach, Fluent's standard finite-difference interpolation schemes are applied to the volume fraction values that were computed at the previous time step.

Mesh impact study

Geometries and meshes are realized with IcemCFD . 

 As shown previously on a Jadim simulation, the influence of the mesh on the results can be important. Therefore, a mesh study is carried out in order to investigate the impact of the mesh in the three-phase flow case. Two simulations are compared, one using a coarse mesh and one using a refined mesh. The main parameters of the simulations are recap in the table below:

Type of mesh Number of cells  velocity Inlet  Gas Reynold time step max  Physical Time
Coarse Mesh        70x88    0.1 m/s      260     $10^{-4}$      2.41
Refined Mesh      140x176      0.1 m/s      260      $10^{-5}$       2.43

In the two pictures below, the oil, gas and water contour are plotted at the same time. As observed for the Jadim simulations, a refined mesh reduces the phenomenon of numerical diffusion. Indeed, the interface between oil and the other phases obtained with a refined mesh is clearer and more accurate than the one obtained with a coarse mesh.

  Oil volume fraction - Coarse mesh            

        Refined mesh                                                                         Coarse mesh

The graph below shows that the oil and water volume fraction at the outlet are not the same for a coarse and a refined mesh. The coarse mesh  leads to an error of 7% on the results.


In conclusion, both local flow features and macroscopic data are affected by the mesh quality.  A refined mesh is necessary in order to have accurate results and to limit the numerical diffusion phenomenon.

Level Set method versus VOF method

Two simulations are carried out, one with the Level Set method and one with the VOF method. The graph below shows the outlet oil flow rate and oil and water remaining in the domain found. We can see that the results obtained are accurate and similar for the two methods: a maximum oil flow rate ratio of 11 is reach, which shows a good flow rate conservation and evolution of the oil and water volume fraction is relevant.

Therefore, the two methods can be used to simulate the problem. However, as the Level set method for a three-phase flow is not available on Jadim, all the Fluent simulations will be carried out with the VOF method.



Analyze of the results

In this part, results of the simulations carried out with Fluent and Jadim are analyzed and compared.

As explained previously, convergence issues make three-phase flow simulations on Jadim difficult to achieve. On another hand, if Fluent seems to avoid these problems, the fact that the access to the code is not possible leads us to be critical about the results.

Therefore, a comparison of similar simulations carried out with these two softwares is important to validate the results.

Two simulations with different velocities and physical properties are realized:

     - In the first simulation, a velocity of 1 m/s is used , oil, gas and water viscosities are multiplied by 100 and the gas density is multiplied by 10.

     - In the second simulation, real properties of oil, gas and water are used with a velocity of 0.1 m/s.


  • Matlab Post-treatment

To visualize the geometry with matlab, the command  les2asc filename geom is used.
To visualize an animation of a serie of timesteps (for instance every 2 timestep from 4 to 10), the command les2asc filename range 4 2 10 is used.

Then, a matlab script file allows to load output jadim variables in matlab, using the function fscanf. Afterwards, values of variables such as volume fraction of oil and gas, velocity, position or pressure are saved into matrices which contain a value for each cell. The loop time gives access to every timestep.

The matlab script is used to plot macroscopic dimentionless data of interest:

- The repartition of oil/water/gas volume fraction in the domain in function of time, plotted in function of a time:

$$oil_{mean}=\frac {\sum \limits_{i=1}^m \sum \limits_{j=1}^l \tau_{oil,ij}}{m \times l}$$                    $$gas_{mean}=\frac {\sum \limits_{i=1}^m \sum \limits_{j=1}^l \tau_{gas,ij}}{m \times l}$$



$m \times l$: numbers of cells of the domain

$\tau_{oil,i}$: Oil volume fraction of cell $i$

$\tau_{gas,i}$: Gas volume fraction of cell $i  $                                

- The outlet oil flow rate in function of time, plotted in function of time:

$$Q_{outlet}=\frac{1}{n}\sum \limits_{i=1}^n \tau_{oil,i} v_i L_{outlet}$$

$n$: number of outlet cells
$v_i$: outlet velocity of cell $i$
$m \times l$: numbers of cells of the domain
- The ratio of remaining and the ratio of oil recovered:
$oil_{initial}=\sum \limits_{i=1}^m \sum \limits_{i=1}^l \tau_{oil,ij} S_{cell}$ 
and $S_{cell}$: area of a cell
  • Paraview post-treatment

Thanks to the command les2par filename range 1st-timestep interval final-timestep, jadim data can be loaded on paraview and it is possible to vizualize the different variable in function of time. In our BEI, paraview was mainly used to create video and to observe the general flow behaviour.

  • Dimentionless Values

As results of this project will be used by Schlumberger, it is important to work with dimentionless variables:

- the outlet oil flow rate is divided by the inlet gas flow rate

- the time is divided by Tref, the time needed to fill up the domain with gas:  

                                        $ t= \frac{ physical-time}{T_{ref}}$

                                        $ T_{ref}=\frac{S}{V_{Inlet}\times D}$


$S$:  surface of the domain without the obstacle

$V_{inlet}$ : velocity of the gas injection

$D$:  lenght of the inlet

Simulation 1 : High viscosities

In this first simulation, oil, gas and water viscosities are multiplied by 100. Indeed, increasing viscosity values should allows Jadim to avoid convergence issues and to reduced numerical diffusion.

The table below recaps the main caracteristics of the simulation:

  Physical properties Velocity inlet Time step Scheme Order Time simulated Quantity of N2 injected
Fluent Viscosities x100 0.1 m/s 10-4 Second Order 2.5 sec 6.25.10-3 m²
Jadim Viscosities x100 0.1 m/s max. 10-4 Second Order 2.5 sec 6.25.10-3 m²

This simulation is carried out for a velocity of 0.1 m/s, 1m/s and 10 m/s. In term of convergence, CPU time and accuracy of the results, the velocity of 1 m/s appeared to be the suitable velocity for this case, for bth Fluent and Jadim simulations. Therefore, the simulations presented in this page have been carried out with a velocity of 1 m/s. 


The inlet is characterized by a gas Reynolds Number of 260 at the inlet. Therefore, the gas jet at the inlet is not turbulent.

If a flow Reynold Number in the domain  is complicated   to calculate, velocities in function of time can be  observed. In  both Jadim and Fluent  simulations,  flow  velocity  does  not exceed 7 m/s. The picture  on the   right shows the velocity field at 1.5 sec with Jadim.

Thus, it is relevant to conclude that the flow  is laminar and  that turbulence models are not needed.



The two videos below show the same case simulated with Fluent and Jadim. As explained previoulsy, gas is injected on the lower right side of the domain, and the outlet is situated on the middle left side. The physical time simulated is 2.5 seconds, which correpond to an injection of %47 of gas i the domain.

Simulation carried out with JADIM


Simulation carried out with Fluent

In the videos, if befferent local flow features can be observed, the global flow in the domain has the same behaviour with Fluent and Jadim. Moreover, the phenomenon of numerical diffusion is not too important and the position of the norrow gas interface is accurate. From the observation of these videos, Jadim seems to handle the three-phase flow simulation and the results seem accurate and relevant.

In order to confirm these observations, macroscopic data of interest are analyzed with the matlab post-treatment. The figure below represents the evolution of oil, water and gas volume fraction in the domain in function of a dimentionless time. Results of Fluent and Jadim simulations are compared.

The first observation is that for both simulations, oil and water volume fraction decrease as gas volume fraction increase, which corresponds to the injection of gas and the disparition of oil and water into the outlet.

Secondly, it shows that no gas injected have left the domain after an injection of 74% of nitrogen. Indeed, the evolution of the gas volume fraction is a  linear curve with represents a constant evolution in the domain.

Finally, oil, pushed by water, starts to leave the domain after an injection of 14% of gas and at the end of the simulation, after an injection of 47% of gas, 52.7% of oil wedged remains in the domain.

The results obtained with Jadim and Fluent are similar: the same amount of oil has left the domain at the end and the evolution of phases repartition is similar during this period of time. These results confirm that  Jadim predicts well flow features in a case of   high viscosities.



Then, the evolution of the outlet oil flow rate is plotted in function of a dimentionless time.  

When oil starts to leave the domain after 14% of nitrogen injected, the dimentionless oil flow rate reach a maximum of 1. The previous videos show that at the time, only oil goes out of the domain. Theses observations lead to the conclution that flow conservation between the inlet and the outlet is respected.

Moreover, like the evolution of the volume fraction, the comparison of the outlet oil flow rate shows that Jadim and Fluent produce similar results in this case.


Conclusion: If local flow features observed on the videos  have differences in function of the CFD software, the study of macroscopic values of interest shows that in this case, Jadim and Fluent give similar results. Jadim seems to gives more accurate results when it does not have to deal with the small viscosity of nitrogen.

Simulation 2: Real properties of oil, gas and water

A second simulation is carried out with real properties of gas. In order to avoid convergence issues, the timestep is reduced to $10^{-5}$ and a small velocity inlet of 0.1 m/s is used. Because of this small velocities, the computational time of the simulation is really high. Therefore, only 5% of nitrogen is injected in the domain before convergence issues appear in the Jadim simulation.

The table below recaps the main caracteristics of the simulation:

  Physical properties Velocity inlet Time step Scheme Order Time simulated Nitrogen Injected (% of the domain)
Fluent Real properties 0.1 m/s $10^{-5}$ Second Order 2.8 sec 5%
Jadim Real properties 0.1 m/s max. $10^{-5}$ Second Order 2.8 sec 5%

The gas injection is characterized by a Gas Reynolds Number of 26. Moreover, as for the simulation 1, the velocity field in the domain observed on paraview shows that the velocity remains small in the domain in function of time.

Thus, it is relevant to conclude that the flow  can be consider as laminar and that turbulence models are not needed.

The two videos below show the same case simulated with Fluent and Jadim. As explained previoulsy, gas is injected on the lower right side of the domain, and the outlet is situated on the middle left side. The physical time simulated is 2.8 seconds, which correponds to a dimentionless time of 0.4 

                                                                                Simulation carried out with Fluent   

                                                                              Simulation carried out with Jadim

On videos, differencies between local flow features can be observed with the Fluent and Jadim simulations.

In the Jadim simulation, the gas interface is not accurate and it seems that Jadim has difficulties to handle flow features. While bubble of gas should be expected, the gas is diffused and seems to mix with oil and water.

In the Fluent simulation, the interface is more accurate but the gas seems to "desappear", probably because of numerical diffusion.

If  numerical diffusion can explained the diffusion of the gas observed with Jadim, another explanation can be found in the scale of the mesh. Indeed, the size of cells is close to 1 mm. Therefore, if the gas forms bubbles smaller than 1 mm, Jadim is not able to represent the topology changes properly. If Fluent seems to provide a accurate interface, it is important to note that it is a really strong commercial software design to always give a result.

The graph below represents the repartition of  phases in function of time. The comparison between Jadim and Fluent shows that after an injection of 5% of nitrogen, the repartition of gas and water in the domains already starts to be different.  If the simulation has been run longer, the repartirion of phases obtained with Fluent and Jadim would have probably been even more different.


If we plot the outlet oil flow rate in function of time, results obtained with Jadim and Fluent are also different. When only oil goes out of the domain, which corresponds to a maximum oil flow rate, the flow rate ratio of the Fluent simulation as a value of 0.9, which represents an error of 10%.


Due to the small velocity, only 5% of nitrogen is injected in the domain and only 4% of oil is recovered.

To conclude, the results obtained by Fluent and Jadim does not seem as accurate as in the Simulation 1. As the viscosity of the gas is small, bubbles of gas are difficult to represent and the interface between gas and other phases is hard to predict.

However, the computational time of Fluent is smaller than the one with Jadim and Fluent is a strong software which had a good capability to handle topology changes. Therefore, the parametric study will only be simulated with Fluent.


Parametric study

A parametric study is carried out to investigate the influence the inlet and the outlet on the oil recovery ratio. For the same geometry, the parametric study investigates:

- the influence of the position of the inlet/outlet 

- the influence of the size of the inlet/outlet 

- the influence of the initial quantity of oil in the domain

 The objective is to observe if the inlet and the outlet have a major impact on the flow behavior and which configuration would provide the greater oil recovery ratio.

Simulations are carried out with Fluent at a velocity inlet of 0.6m/s, using real physical properties and a first order explicit scheme. 

The oil viscosity value is increased in order to reduce convergence issues and to help Fluent to handle topology changes. Indeed, as a crude oil viscosity value can vary from 1 to 1000mPa.s,  the oil viscosity  is increased to 420mPa.



Optimization of gas injection

The first part of the parametric study investigates how the position and size of the injection well can influence the results.

Two positions are tested, on the bottom and lateral  wall of the domain, with an intlet lenght constant equals to 0.025m. Then the size of the inlet is reduced by two as it is shown in the pictures below:

Simulation parameters are listed in the table below:

  Velocity Inlet  Percentage of gaz injection  Percentage of oil recovered
Inlet 1  0,6 m/s 75% 68,5%
Inlet 2 0,6 m/s 75% 77%
inlet 3 1,2 m/s 75% 66%

The graph below represents the percentage of oil in the domain in function of time. We can observe that after injecting 70% of gas in the domain, the quantity of oil remains constant: the oil recovery ratio is maximum. There is a maximum quantity of gas injected after which oil viscosity ratio can not be increased.

After an injection of 75% of nitrogen, the same amount of oil is recovered in the case of Inlet 1 and 3 . Therefore, the size of the injection well does not  influence the oil recovery ratio.

On the contrary, the percentage of oil recovered by Inlet 2 is increase by 12.6% compared to Inlet 1. The position of Inlet 2 is more efficient than the position of Inlet 1 to recover oil. Thereforethe position of the injection well have an impact on the oil recovery ratio.

We recommand to investigate carefuly the position of the injection well before developping the nitrogen injection recovery method in the Tarim Basin, as it impacts the oil viscosity ratio.




Influence of the outlet

The influence of the oulet is also investigated. Indeed, connexions between caves can be completly differente. Thus, it is important to see how the flow behaves for a different outlet.

Two positions are tested, with an outlet lenght of 0.05m. Then,  the size if the outlet is increased by two as it is shown in the pictures below. 

In these configurations, the inlet is the same used in the reference geometry.

Simulation parameters are listed in the table below:

  Velocity Inlet  Percentage of gaz injection  Percentage of oil recovered
Outlet 1  0.6 m/s 75% 68,5%
Outlet 2 0.6 m/s 75% 69%
Outlet 3 0.6 m/s 75% 70%

The graph below represents the percentage of oil in the domain in function of time. We can observe that after injecting 75% of gas in the domain, the quantity of oil remains constant: the oil recovery ratio is maximum. There is a maximum quantity of gas injected after which oil viscosity ratio can not be increased.

After an injection of 75% of nitrogen, the same amount of oil is recovered in the case of Outlet 1, 2 and 3 . Therefore, the size a,d the position of the injection well does not  influence the oil recovery ratio.



Influence of the initial amount of oil

Finally, the parametric study investigates how the initial quantity of oil in the domain can influence the results.

Two simulations are carried out: one with the size of the inlet is reduced by two as it is shown in the pictures below:

         37% of  oil initialy in the domain                                18% of  oil initialy in the domain

Simulation parameters are listed in the table below:

  Velocity Inlet Percentage of gaz injection Percentage of oil recovered
37% oil initial 0.6 m/s 75% 78%
18% oil initial 0.6 m/s 75% 66.5%

After an injection of 75% of nitrogen, the percentage of oil recovered in function of the initial quantity of oil is higher when the initial quantity of oil is important. 

These results were expected as every recovery method have difficulties to recover oil when the initial quantity is really small.

 Therefore, the initial amount of oil in the domain impact the efficiency of the recovery method. The oil recovery ratio will not be the same for every cave and will depends on the initial amound of oil present in the cave.



Complex Geometry

To conclude this project, the simulation of gas injection have been carried out, using a complex geometry investigated by our partner group.

Indeed, they  simulated the injection of water in a cave full of oil in order to see the amount of oil which can be recovered and where the remaining oil is wedged. The geometry and mesh caracteristics can be found in the website of our partner group:

Here, a injection of gas in this geometry is simulated. The initial amount of oil corresponds to the quantity of oil wedged into the cap of the cavity at the final time of the two-phase flow simulation.

The table below recaps the main characteristics of the simulation:

Inlet velocity Physical properties Mesh Time step Physical time simulated
0.6 m/s real properties refined $10^{-4}$ 5 s

On the video,we can see that while the gas is injected into the domain, water and oil are driven out. Indeed, the graphs below shows that their volume fractions in the domain decrease with time.

Because of its lighter density, gas first occupies the first cave, pushing the oil out of the cavity. The oil is therefore drived out to the outlet, increasing the oil recovery ratio.

After the gas has fulfilled the cave, some of the gas flows over the obstacle between these two caves and drives the oil in the second cave out. Finally, it seems that after a certain time,  a maximum oil recovery ratio has been reached.

The first figure shows the outlet oil and water flow rate in function of time and the second figure shows the oil, water and gas volume fraction in the domain. 

From these figures, we can observe that when 25% of gas is injected in the domain, gas starts to leave the domain. Moreover, the repartition of phase does not change: a stable state is reached where no oil can be recovered anymore.

In the end, the simulation show that 34.3% oil was recovered in this procedure. Therefore, after water injection recovery, gas injection would be a good choice to improve oil recovery efficiency.