EP1889999B1 - Method of optimising the assisted recovery of a fluid placed in a porous medium by front follow-up - Google Patents

Description

The present invention relates to a method for optimizing the exploitation of a heterogeneous porous medium in the assisted fluid recovery framework in place, by determining the position of the front separating a sweeping fluid and the fluid in place.

In particular, the method according to the invention makes it possible to determine the position of a front separating two immiscible fluids in motion, without systematic updating of the pressure field, so as to obtain a simulation of the multiphase flows in the porous medium sufficiently fast and precise to obtain quantitative information allowing an optimal exploitation of the medium.

The method finds applications in particular for the exploitation of oil or gas deposits, or the exploitation of underground storage of gas for example.

Presentation of the prior art

All the notations used to describe the prior art and the invention are defined at the end of the description.

Knowing how to describe and simulate multiphase flows in underground reservoirs is the basis of the job of reservoir engineers in oil and gas companies (and similarly in water supply companies). The presence of heterogeneities in the subsoil and the increasing complexity of drainage systems make any simple analytical solution impossible, which imposes the development of digital solutions using a mesh model. The simulation of multiphase flows in a heterogeneous porous medium can then require large computer resources, in particular when the numerical model of the medium in question is highly detailed. This is particularly the case in tank engineering, in the oil field. This cost comes mainly from the resolution of the large linear systems resulting from the equation governing the pressure, which must be updated by following the displacement of the fluids, if one wishes to have a solution having a good precision.

$${\mathbf{u}}_{\mathit{w}}\mathit{=}\mathit{-}{\mathit{\lambda }}_{\mathit{w}}\mathrm{\nabla }{\mathit{p}}_{\mathit{w}},{\mathit{\lambda }}_{\mathit{nw}}=\frac{\mathbf{K}{k}_{\mathit{rw}}}{{\mathit{\mu }}_{\mathit{w}}}$$

$$\mathbf{u}={\mathbf{u}}_{\mathit{nw}}+{\mathbf{u}}_w$$

où les indices nw et w note respectivement les fluides non mouillant et mouillant. La différence des pressions entre les deux fluides est notée p_c (S) = p_nw - p_w. On a la relation : $$S_{\mathit{nw}}+S_w=1$$

In the context of a two-phase flow in a heterogeneous porous medium, it is considered the displacement of a fluid in place (for example the oil) under the effect of the injection of another fluid (for example from the water). The generalized Darcy law governing the movement of fluids is written, using standard assumptions and notation, as follows: $${\mathbf{u}}_{\mathit{nw}}\mathit{=}\mathit{-}{\mathit{\lambda }}_{\mathit{nw}}\mathrm{\nabla }{\mathit{p}}_{\mathit{nw}},{\mathit{\lambda }}_{\mathit{nw}}=\frac{\mathbf{K}{k}_{\mathit{rnw}}}{{\mathit{\mu }}_{\mathit{nw}}}$$

$${\mathbf{u}}_{\mathit{w}}\mathit{=}\mathit{-}{\mathit{\lambda }}_{\mathit{w}}\mathrm{\nabla }{\mathit{p}}_{\mathit{w}},{\mathit{\lambda }}_{\mathit{nw}}=\frac{\mathbf{K}{k}_{\mathit{rw}}}{{\mathit{\mu }}_{\mathit{w}}}$$

$$\mathbf{u}={\mathbf{u}}_{\mathit{nw}}+{\mathbf{u}}_w$$

where the indices nw and w respectively note the non-wetting and wetting fluids. The difference in the pressures between the two fluids is noted p _c ( S ) = p _nw - p _w . We have the relationship: $$S_{\mathit{nw}}+S_w=1$$

As a result, one can write all the functions dependent on the saturation S in terms of the saturation S _w .

Neglecting the capillary pressure, we have $$p_{\mathit{nw}}=p_w=p$$

avec $$\nabla \cdot \mathbf{u}=0$$

The total speed then takes the form $$\mathbf{u}=-\mathit{\lambda }\left(S_w\right)\nabla p$$

with $$\nabla \cdot \mathbf{u}=0$$

This system of equations (6) and (7) defines the pressure equation.

In equations (1) and (2), the functions k _r are the relative permeabilities. We have k _r = k _r ( S _w ). The functions λ _nw , λ _w are the respective mobilities of the fluids, and λ ( S _w ) = λ _nw ( S _w ) + λ _w ( S _w ) is the total mobility, which therefore depends explicitly on S _w .

$$\varphi \frac{\partial S_{\mathit{w}}}{\partial t}+\nabla \cdot {\mathbf{u}}_{\mathit{w}}=0$$

The laws of conservation of the mass relating to each of the fluids are written $$\phi \frac{\partial S_{\mathit{nw}}}{\partial t}+\nabla \cdot {\mathbf{u}}_{\mathit{nw}}=0$$

$$\phi \frac{\partial S_{\mathit{w}}}{\partial t}+\nabla \cdot {\mathbf{u}}_{\mathit{w}}=0$$

This system of equations (8) and (9) defines the saturation equation. Here φ is the porosity (assumed to be uniform in the tank), t is the time, and u _{n, nw} is the velocity of the fluid considered.

One of the difficulties in solving these equations comes from the coupling between equations (6) to (9) which couple the saturation to the pressure field. These effects are well known and drive the development of possible viscous instabilities. In particular, in the case 1D, considering the scanning of a medium initially saturated with oil by water, this system of equations can be solved by the method of the characteristics, whose solutions are characterized by the existence of oscillations, corresponding to saturation leaps propagating in the medium.

The discontinuity is characterized by a saturation in water passing from S _wi to S _f , the saturation at the front whose value is given by a Rankine-Hugoniot type condition. In the oil context, the construction of Welge's tangent makes it possible to determine S _f graphically. Behind this front, a "rarefaction wave" ensures the transition between S _f and the water saturation at the injection.

Several techniques are known for determining the simulation of multiphase flows in a heterogeneous porous medium. These methods all rely on the resolution of the following system of equations: $$\{\begin{array}{l}\begin{array}{l}\mathrm{Pressure\; equation}:\\ \mathbf{u}=-\mathit{\lambda }\left(S_w\right)\nabla p\left(10\right)\\ \nabla \cdot \mathbf{u}=0\left(11\right)\end{array}\\ \begin{array}{l}\mathrm{Saturation\; equation}:\\ \phi \frac{\partial S_{\mathit{nw}}}{\partial t}+\nabla \cdot {\mathbf{u}}_{\mathit{nw}}=0\left(12\right)\\ \phi \frac{\partial S_{\mathit{w}}}{\partial t}+\nabla \cdot {\mathbf{u}}_{\mathit{w}}=0\left(13\right)\end{array}\end{array}$$

The techniques differ in the method used to solve this system. In the general case, in a complex heterogeneous medium, this type of system of equations must be solved numerically, a reference of which is: B. Noetinger, PE Spesivtsev and E. Teodorovich / Stochastic Analysis of a Displacement Front in a Randomly Heterogeneous Medium, Fluid Dynamics Vol. 41, No. 5, 2006 pp 830-842 .

The oil companies most often use conventional discretization techniques, such as finished volumes ( Aziz, Kb., Settary, A.: Petroleum reservoir simulation. Applied science publishers, London, 1979 ). Regarding temporal discretization, the choices are many: if we favor digital robustness, we can be totally implicit in pressure as in saturation. If one prefers to avoid oscillations of the solution, and a better precision, one will choose a schema implicit in pressure, explicit in saturation. We generally observe the existence of fronts characterized by the same discontinuity of saturation whose value passes abruptly from the saturation irreducible in water, S _wi , to the saturation of the front, denoted S _f . These discontinuities are due to the hyperbolic character of the saturation transport equation, and are well described analytically in the case of 1D displacement: this is the Buckley-Leverett theory, a detailed description of which can be found in Charles-Michel Marle. Multiphase flows in porous media. Second edition revised and expanded 1972. Technip Publishing, Par est.

Whatever the method chosen, regular updates of the equation in pressure are necessary, so as to estimate the speeds with precision.

In another class of techniques, which are becoming more and more used in the operating companies, we find the so-called current-line techniques, a reference of which is: RP Batycky, MJ Blunt, and MR Thiele, A 3d field-scale streamline-based reservoir simulator, SPERE, 11, 246-254 (1997) .

In these methods the saturations are updated by following the movement of fluids in their movement along the streamlines, parametric curves defined by dx / dt = u (x, t). A change of variable can be reduced to the Buckley-Leverett theory previously mentioned. However, in order to correctly calculate the velocity u ( x , t) the pressure is solved as many times as precision requires, on a Cartesian grid. These methods are faster than conventional methods. Nor do they possess the versatility and robustness to deal with complex compositional problems, where the traditional method remains the only possible solution. Expensive updates of the pressure field are not avoided.

Finally, there are the front-end techniques ( R Juanes a, d KA Lie "A front tracking method for efficient simulation of miscible gas injection." paper SPE 92298, Houston Jan. 31 to 2 feb. 2005 ), whose idea is to follow the front by a Lagrangian approach, which requires the same updates of the pressure. In this type of method, the description of the multiphase flows in the underground reservoirs consists in determining the position of the front (also called interface) separating two immiscible fluids in motion: a fluid in place (oil for example) and a fluid of sweeping (of water for example), also called injected fluid

The evolution of the front in the reservoir during the flow is considerably influenced by the coupling (called viscous coupling by the reservoir engineers) between the pressure field and the saturation field ( Saffman PG and G. Taylor. The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid. Proc. Royal Society of London, A245: 312-329, 1958 , Noetinger B., V. Artus, and L. Ricard. Dynamics of the water-front for two-phase, immiscible flows in heterogeneous porous media. 2-isotropic media. Transport in porous media, 56: 305-328, 2004 ).

In particular, when the injected fluid is less viscous and therefore more mobile at the forehead than the fluid in place, the viscous instabilities will always favor the flow of fluids in the most permeable layers of the medium. The breakthrough time through these is much faster than in the rest of the tank. On the contrary, if the injected fluid is less mobile, the viscous coupling can slow it down in the initially faster layers, thus compensating for the differences in permeability due to stratification. A stationary front appears then ( V. Artus, B. Noetinger, and L. Ricard. Dynamics of the water-front for two-phase, immiscible flows in heterogeneous porous media. 1-layered media. Transportation in Porous Media, 56: 283-303, 2004 ).

It is known that the stability of the front following small disturbances is determined by the ratio of the viscosities of the fluids. If the less viscous fluid is displaced by the more viscous fluid, then the front is stable, and vice versa. If the situation is unstable, then phenomena, called "viscous fingers", develop over time.

The effects of relative permeability are taken into account when considering Bucley-Leverett displacements in a homogeneous medium. Such displacements are described by a hyperbolic saturation equation. A solution to this equation by the method of the provided characteristics of solutions that have oscillations characterized by a condition of Rankine Hugoniot. It has been shown that the stability of the front where the saturation has a "shock" is determined by a frontal mobility ratio, noted M _f : if M _f <1, then the front is stable with regard to small perturbations, and vice versa . This frontal mobility ratio can be approximated to the viscosity ratio using the relative permeability function.

où $$c_0=\frac{u_0}{\varphi }\frac{f_w\left(S_f\right)-f_w\left(S_{\mathit{wr}}\right)}{S_f-S_{\mathit{wr}}}=\frac{u_0}{\varphi }{f}_w^{ʹ}\left(S_f\right)$$

Avec:

M_f :

le rapport de mobilité frontal

u ₀ :

la vitesse de filtration moyenne le long de la direction de l'axe X

φ :

la porosité

f_w :

la fonction de Buckley-Leverett représentant le flux fractionnaire en eau.

S_f :

Saturation en eau au front

S_wr :

Saturation maximale en eau

More recently, this result has been confirmed in linearized stability analysis of the water injection process. In this context, it has been shown that the evolution of the Fourier transform of the fluctuations of the position of the front, fluctuations noted δ h ( q, t ), in a homogeneous medium is given by: $${\partial }_t\mathit{.DELTA.h}\left(qt\right)={\mathrm{vs}}_0\left|q\right|\frac{M_f-1}{M_f+1}\mathit{.DELTA.h}\left(qt\right)$$

or $${\mathrm{vs}}_0=\frac{u_0}{\phi }\frac{f_w\left(S_f\right)-f_w\left(S_{\mathit{wr}}\right)}{S_f-S_{\mathit{wr}}}=\frac{u_0}{\phi }{f}_w^{\text{'}}\left(S_f\right)$$

With:

M _f :

the frontal mobility report

u ₀ :

the average filtration rate along the direction of the X axis

φ:

porosity

f _w :

the Buckley-Leverett function representing the fractional flow of water.

S _f :

Water saturation at the front

S _wr :

Maximum water saturation

To summarize, all the aforementioned numerical techniques have a very high cost of computation time which comes mainly from the large number of resolutions of large linear systems from the equation governing the pressure, which must be updated by following the displacement fluids with good accuracy.

1. 1 Échantillonner rapidement l'espace des paramètres d'étude sur lesquels porte la modélisation, de façon à optimiser un ou plusieurs critères économiques et techniques d'exploitation (décision d'exploitation, choix de la position des puits, d'un procédé de récupération etc.).
2. 2 Savoir modifier de façon cohérente le modèle géologique afin de caler au mieux les données de production observées, y compris des données de sismique répétée. Ceci peut permettre de forer en s'adaptant aux hétérogénéités géologiques, de localiser les fluides, et donc de mieux contrôler le scénario de récupération.
3. 3 Pouvoir estimer des incertitudes en effectuant des simulations de Monte Carlo, de façon à tester le rôle d'hétérogénéités, ou de paramètres mal connus caractérisant le sous-sol. Ceci permet de quantifier le niveau de risque lié à l'exploitation d'un réservoir et donc de redéfinir des paramètres d'exploitation voire les paramètres économiques.

As a result, the existing solutions do not make it possible to obtain quantitative answers to practical questions that the petroleum engineer poses as soon as tank simulation studies are carried out, namely:

1. 1 Quickly sample the space of the study parameters being modeled, in order to optimize one or more economic and technical criteria exploitation (exploitation decision, choice of the position of the wells, a recovery process, etc.).
2. 2 Knowing how to modify the geological model in a coherent way in order to better calibrate the observed production data, including repeated seismic data. This can make it possible to drill by adapting to the geological heterogeneities, to locate the fluids, and thus to better control the recovery scenario.
3. 3 To be able to estimate uncertainties by carrying out Monte Carlo simulations, so as to test the role of heterogeneities, or poorly known parameters characterizing the subsoil. This makes it possible to quantify the level of risk associated with the operation of a reservoir and thus to redefine operating parameters or even the economic parameters.

An exhaustive study incorporating these various aspects (sensitivity study, data calibration and estimation of uncertainties) can then potentially require numerous calls to a flow simulator (dedicated software), hence the interest of having tools fast calculation, accuracy compatible with the accuracy of the input data.

The method according to the invention makes it possible to determine the position of a front separating two immiscible fluids in motion in a heterogeneous porous medium, without systematic updating of the pressure field so as to obtain a simulation of the multiphase flows in the sufficiently fast porous medium. to obtain quantitative information allowing an optimal exploitation of the reservoir, by the determination of technical operating parameters and / or economic parameters.

The method according to the invention

• on détermine une vitesse et une direction d'écoulement d'au moins un desdits fluides, en utilisant un simulateur d'écoulement pour résoudre une équation en pression ;
• on définit une relation décrivant une position d'un front séparant lesdits fluides, en s'affranchissant d'un couplage visqueux à l'aide de la théorie des perturbations, et en prenant en compte des fluctuations de vitesse dans la direction d'avancée du front et des fluctuations de vitesse dans la direction perpendiculaire à celle d'avancée du front ; et pour différents pas de temps,
• on reconstruit la position du front dans ladite grille à l'aide d'une discrétisation de ladite relation et d'une transformée de Fourier rapide ; et
• on optimise l'injection dudit fluide de balayage en fonction de la position dudit front.

The invention relates to a method for optimizing the exploitation of a heterogeneous porous medium containing a fluid, said fluid in place. The operation is performed by injection into the medium of a second fluid, said scanning fluid, to cause a fluid flow in place. These two fluids are immiscible. According to the method, the medium is discretized into a grid consisting of a set of meshes. The method has the following steps:

• determining a velocity and a flow direction of at least one of said fluids, using a flow simulator to solve a pressure equation;
• defining a relationship describing a position of a front separating said fluids, dispensing with a viscous coupling using the perturbation theory, and taking into account speed fluctuations in the direction of advance of the forehead and velocity fluctuations in the direction perpendicular to that of advancing the front; and for different time steps,
• the position of the front in said grid is reconstructed by means of a discretization of said relation and a fast Fourier transform; and
• the injection of said scanning fluid is optimized as a function of the position of said front.

• un premier terme représentant un couplage visqueux décrivant la perturbation due à la variation de saturation ;
• un second terme décrivant les perturbations de la position du front provoquées par les hétérogénéités du milieu.
  Ce premier terme peut être obtenu en considérant le milieu homogène et les déplacements du front de type Bucley-Leverett. Le second terme peut quant à lui, être obtenu en utilisant le fait que le front est une surface matérielle, et en représentant la vitesse comme une somme d'une vitesse moyenne avec des fluctuations de vitesses dues aux hétérogénéités.
  Le premier terme peut être fonction d'au moins les paramètres suivants : un rapport de mobilité frontal, une vitesse de filtration moyenne le long d'une direction d'avancée du front, une porosité du milieu, une fonction de Buckley-Leverett représentant un flux fractionnaire en eau, une saturation en eau au front, une saturation maximale en eau, le vecteur d'onde dans l'espace de Fourier.
  Le second terme peut être fonction d'au moins les paramètres suivants : une vitesse du front non perturbé, une vitesse totale moyenne, des perturbations des composantes de la vitesse totale dans la direction d'avancée du front et perpendiculaire à cette direction.
  Enfin pour optimiser l'injection, on peut déterminer en chaque maille et pour différents pas de temps la saturation d'au moins un desdits fluides à partir de la position dudit front.
  D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après d'exemples non limitatifs de réalisations, en se référant aux figures annexées et décrites ci-après.

Depending on the method, the defined relationship may include the following two terms:

• a first term representing a viscous coupling describing the disturbance due to saturation variation;
• a second term describing the disturbances of the position of the front caused by the heterogeneities of the medium.
  This first term can be obtained by considering the homogeneous medium and displacements of the Bucley-Leverett type front. The second term can be obtained by using the fact that the front is a material surface, and representing the speed as a sum of an average speed with velocity fluctuations due to heterogeneities.
  The first term may be a function of at least the following parameters: a frontal mobility ratio, an average filtration rate along a direction of advance of the front, a medium porosity, a Buckley-Leverett function representing a fractional flow of water, water saturation at the front, maximum saturation in water, the wave vector in the Fourier space.
  The second term may be a function of at least the following parameters: an undisturbed front velocity, an average total velocity, disturbances of the components of the total velocity in the forward direction of the front and perpendicular to this direction.
  Finally, to optimize the injection, it is possible to determine in each mesh and for different time steps the saturation of at least one of said fluids from the position of said front.
  Other characteristics and advantages of the method according to the invention will appear on reading the following description of nonlimiting examples of embodiments, with reference to the appended figures and described below.

Summary presentation of figures

• The Figures 1A and 1B show a comparison between the result provided by the method according to the invention ( Fig. 1A ), and a current line simulation ( Fig. 1B ).

Detailed description of the method

The method according to the invention makes it possible to optimize the exploitation of a heterogeneous porous medium in the framework of enhanced oil recovery (EOR) for example. This technique consists in injecting, in a petroleum tank, a sweeping fluid (CO ₂ , water) to cause a flow of a fluid in place (oil) contained in the reservoir and that one seeks to extract. These two fluids must of course be immiscible.

In this framework of assisted recovery, we consider a situation where there is an interface, called a front, separating these two fluids. At the scale of the tank the width of this front is small compared to the size of the oil tank studied, and for the sake of simplicity it is assumed that this front is narrow relative to the size of the tank. We focus the study on the dynamics of this front in a heterogeneous porous medium, neglecting capillary effects and gravity.

• on discrétise le milieu en un ensemble de maille ;
• on détermine le champ de vitesse pour chacun des fluides ;
• on définit une relation décrivant la position du front en milieu hétérogène ;

et à chaque pas de temps,

• on estime la position du front à l'aide de la relation et des vitesses ;
• on optimise l'injection du fluide de balayage en fonction de la position du front.

The method according to the invention makes it possible to optimize the injection by determining, for different instants, called time steps, the position of the front separating the two fluids in motion in a heterogeneous porous medium (the reservoir for example). The method is based on an original technique of systematic updating of the pressure field which has the advantage of being both accurate and fast, since the method requires only a very limited call to a simulator. performing both long and complex calculations. To do this, the update of the pressure field, due to the change of geometry of the front, is performed semi-analytically using Fast Fourier Transform (FFT) techniques. This saves most of the calculation cost due to these updates of the pressure field, hence the speed of the method. The method mainly comprises the following steps:

• the medium is discretized into a set of stitches;
• the velocity field is determined for each of the fluids;
• a relation describing the position of the front in a heterogeneous medium is defined;

and at each time step,

• the position of the forehead is estimated using relation and velocities;
• the injection of the sweeping fluid is optimized as a function of the position of the front.

Discretization of the environment

The discretization of the medium is a classic step and well known to specialists. It makes it possible to represent the structure of the medium in a set of meshes, characterized by their dimensions and their geographical positions. The calculations are done for each mesh. This step is essential to use flow simulators (the simulator is a dedicated software). This set of mesh is called "grid". This grid, discretizing the medium can be rectangular size L _x × L _y meshes.

Determination of a heterogeneous speed field

During this step, a heterogeneous velocity field is determined in the context of a monophasic flow. To do this, we solve numerically the pressure equation (equations (6) and (7)), which is obtained by combining Darcy's law and the incompressibility equation using a monophasic flow simulator. . The velocity vector of each of the fluids, i.e., their velocity and direction of flow, can be determined at a time t = 0.

For a 2D domain (the 3D case is treated in exactly the same way.) Of size L _x × L _y whose permeability map is K ( x, y ) and for which the main flow is parallel to an X axis , y identifying the transverse coordinate (Y-axis), the resolution of the pressure equation is to solve the following equation: $$\nabla \cdot (\mathit{\lambda }\left(S_w\left(x,\mathrm{there},t=0\right)\cdot \nabla P\left(x\mathrm{there}\right)\right)=0$$

From these equations, we deduce the total velocity field u : u _nw , the velocity vector of the non-wetting fluid, and u _w , the velocity vector of the wetting fluid.

This resolution, which uses a simulator, is performed only once according to the method. The pressure equation is solved at the beginning, that is to say for t = 0, it is no longer solved at subsequent time steps ( t > 0).

We can use a standard five-point finite volume discretization to solve this equation using the simulator: the pressure is estimated at the center of the meshes, the fluxes between meshes involving transmissivities calculated as the harmonic mean of the permeabilities of the two meshes. Darcy's law then makes it possible to determine fluxes and speeds.

• Le Ravalec, M., Noetinger, B., and Hu, L.Y. [2000] The FFT moving average (FFT-MA) generator: An efficient numerical method for generating and conditioning Gaussian simulations. Mathematical Geology 32(6), 701-723 .

The permeability map K ( x, y ) can be generated using a random field generator, as described in the following document for example:

• The Ravalec, M., Noetinger, B., and Hu, LY [2000] The FFT moving average (FFT-MA) generator: An efficient numerical method for generating and conditioning Gaussian simulations. Mathematical Geology 32 (6), 701-723 .

Definition of a relation describing the position of the front in heterogeneous medium

The basic idea is therefore to avoid recalculating at each iteration in time the flow in the whole computational domain, by restricting itself to the description of the movement of the front, described by the function δ h ( y , t ). This is made possible by an approximation making it possible to relate the modification of the pressure field δ p ( x, y, t ), and thus of the velocities, to the distortion of the front δ h ( y , t ). We can eliminate the pressure of the problem and we are then brought to the resolution of a differential system involving δ h ( y, t ).

où $$c_0=\frac{u_0}{\varphi }\frac{f_w\left(S_f\right)-f_w\left(S_{\mathit{wr}}\right)}{S_f-S_{\mathit{wr}}}=\frac{u_0}{\varphi }{f}_w^{ʹ}\left(S_f\right)$$

Avec :

M_f :

le rapport de mobilité frontal

u ₀ :

la vitesse de filtration moyenne le long de la direction de l'axe x

φ :

la porosité

f_w :

la fonction de Buckley-Leverett représentant le flux fractionnaire en eau.

S_f :

Saturation en eau au front

S_wr :

Saturation maximale en eau

q :

vecteur d'onde dans l'espace de Fourier

We have seen that in a homogeneous medium, the evolution of the Fourier transform of the fluctuations of the position of the front, fluctuations noted δ h ( q, t ), is given by: $${\partial }_t\mathit{.DELTA.h}\left(qt\right)={\mathrm{vs}}_0\left|q\right|\frac{M_f-1}{M_f+1}\mathit{.DELTA.h}\left(qt\right)$$

or $${\mathrm{vs}}_0=\frac{u_0}{\phi }\frac{f_w\left(S_f\right)-f_w\left(S_{\mathit{wr}}\right)}{S_f-S_{\mathit{wr}}}=\frac{u_0}{\phi }{f}_w^{\text{'}}\left(S_f\right)$$

With:

M _f :

the frontal mobility report

u ₀ :

the average filtration rate along the direction of the x- axis

φ:

porosity

f _w :

the Buckley-Leverett function representing the fractional flow of water.

S _f :

Water saturation at the front

S _wr :

Maximum water saturation

q :

wave vector in Fourier space

In a heterogeneous medium, the problem is more complex. The heterogeneities of the permeability field initiate perturbations of the front, that is to say that the front is deformed according to the permeability of the strata. Thus, these disturbances can increase or decrease depending on the ratio of viscosities. This is the problem of viscous coupling. The problem is non-linear and a mathematical analysis proves very complicated.

According to the invention, it is proposed to take into account the heterogeneities, on the one hand by avoiding the viscous coupling with the aid of the perturbation theory, and on the other hand taking into account the transverse fluctuations of the speed, that is to say the fluctuations in a Y axis perpendicular to the X axis of advance of the front.

Decoupling by the technique of disturbances

The method according to the invention estimates the coupling between the disturbances of velocities δ u and the fluctuation of the position of the front δ h semi analytically. In this sense, the method uses the perturbation theory.

From a heuristic point of view, the perturbation theory is a general mathematical method that allows to find an approximate solution of a mathematical equation ( E _λ ) dependent on a parameter λ when the solution of the equation ( E ₀ ) , corresponding to the value λ = 0, is exactly known. The mathematical equation ( E _λ ) can be an algebraic equation, a differential equation, an equation with eigenvalues, ... The method consists in looking for the approximate solution of the equation ( E _λ ) in the form of a development in series of the powers of the parameter λ, this approximate solution being supposed to be a better approximation of the exact solution, but unknown, that the absolute value of the parameter λ is "smaller".

Thus, according to the invention, the perturbation theory leads to consider locally the speed perturbation as a sum of two disturbances: $$\delta \mathbf{u}\left(x\mathrm{there}t\right)=\delta {\mathbf{u}}_{\mathrm{at}}\left(x\mathrm{there}t\right)+\delta {\mathbf{u}}_b\left(x\mathrm{there}t\right)$$

The first term, δ u _, describes the disturbance due to the saturation variation, and the second term δ u _b describes the perturbation due to the heterogeneity of the medium.

Avec : $$A=\frac{M_f-1}{M_f+1}$$

δu_bx :

les fluctuations de la vitesse de filtration dans la direction d'avancée du front (axe X)

x_f :

l'abscisse de la position du front sur l'axe X

According to this decomposition, we can describe the position of the front during its advance in a heterogeneous porous medium by the following equation: $${\partial }_t\mathit{.DELTA.h}\left(qt\right)={\mathrm{vs}}_0\mathrm{AT}\left|q\right|\mathit{.DELTA.h}\left(qt\right)+\frac{{\mathrm{vs}}_0}{u_0}{\mathit{\delta u}}_{\mathit{bx}}\left(x_{f}qt\right)$$

With: $$\mathrm{AT}=\frac{M_f-1}{M_f+1}$$

δ u _bx :

the fluctuations of the filtration speed in the forward direction of the front ( X axis)

x _f :

the abscissa of the position of the front on the X axis

Thus, using this formulation of the position of the front, it is possible to decouple the problem. In this way, for patency fields with small variances, the stability of the front undergoing small disturbances is determined by the frontal mobility ratio.

Taking into account transversal fluctuations in speed

According to the invention, the propagation of the front in a heterogeneous medium, that is to say in a heterogeneous speed field, is simulated by modeling the influence of the viscous effects analytically. The heterogeneous velocity field was obtained using monophasic flow simulation (step 2 of the method).

To simplify the explanations, we place ourselves in two dimensions. The front separating the two immiscible fluids is then described by a 2D function, x ( y , t ) = h ( y, t ). We can write this equation in the form: $$F\left(x\mathrm{there}t\right)=0$$

By solving it with respect to x, we have: $$F_1\left(x\mathrm{there}t\right)=h\left(\mathrm{there}t\right)-x=0$$

Then, using the fact that the surface in question is material, ie composed of fluid elements that follow the local movement of the fluid, and using the expression for the velocity of propagation of the front for Buckley-Leverett displacements , we obtain : $${\partial }_t\mathit{h}\left(\mathrm{there}t\right)={\mathrm{vs}}_0\left(\mathbf{u}\cdot \nabla \mathit{\phi }\right){|}_{\mathit{\phi }=0}$$

Avec : r = {x,y}, u ₀ = {u ₀,0} et δu={δu_x ,δu_y }Now considering the heterogeneity of the velocity field, velocity is represented as a sum of two terms. The first represents the average, and the second the fluctuation due to heterogeneities: $$\mathbf{u}\left(\mathrm{r}\right)={\mathbf{u}}_0+\delta \mathbf{u}\left(\mathrm{r}\right)$$

With: r = { x, y }, u ₀ = { u ₀ , 0} and δ u = {δ u _x , δ u _y }

The perturbations of the velocity field cause the perturbations of the front h ( y, t ) = h ₀ ( t ) + δ h ( y, t ). Substituting this decomposition in equation (19) and extracting the average term using the relation ∂ _t h ₀ ( t ) = c ₀ , which is valid for Buckley-Leverett type displacements, we obtain the following equation , describing the disturbances of the position of the forehead caused by the heterogeneities of the medium: $${\partial }_t\mathit{h}\left(\mathrm{there}t\right)=\frac{{\mathrm{vs}}_0}{u_0}\left(\delta u_x\left(x_f\mathrm{there}\right)-\delta u_{\mathrm{there}}\left(x_f\mathrm{there}\right){\partial }_{\mathrm{there}}\mathit{.DELTA.h}\left(\mathrm{there}t\right)\right)$$

Le premier terme (à gauche) est l'expression du couplage visqueux et peut être calculé avantageusement dans l'espace de Fourier. Il n'y a donc plus lieu d'estimer ce couplage par une coûteuse résolution de système linéaire. Dans tout ce qui suit, les perturbations de vitesses δu_x et δu_y peuvent maintenant être supposées découplées de l'équation en saturation. Le second terme (à droite) est obtenu à l'aide d'un simulateur d'écoulement monophasique (étape 2 de l'invention).Thus, by applying the theory of perturbations, the theory of material surfaces, and using the expression of the propagation velocity of the front for the Buckley-Leverett displacement, we can describe the position of the front in a porous medium heterogeneous in the following way (δ h has been substituted by h in this expression for the sake of clarity): $${\partial }_t\mathit{h}\left(\mathrm{there}t\right)={\mathrm{vs}}_0\mathrm{AT}\int \frac{\mathit{dq}}{2\mathit{\pi }}\left|q\right|h\left(qt\right)e^{\mathit{iqy}}+\frac{{\mathrm{vs}}_0}{u_0}\left(\delta u_x\left(x_f\mathrm{there}\right)-\delta u_{\mathrm{there}}\left(x_f\mathrm{there}\right){\partial }_{\mathrm{there}}\mathit{h}\left(\mathrm{there}t\right)\right)$$

The first term (left) is the expression of the viscous coupling and can be calculated advantageously in the Fourier space. It is therefore no longer necessary to estimate this coupling by an expensive linear system resolution. In all that follows, the velocity perturbations δu _x and δu _y can now be assumed to be uncoupled from the saturation equation. The second term (right) is obtained using a monophasic flow simulator (step 2 of the invention).

The next step in the method is to propose a technique for solving the full-frontal equation.

Estimation of the position of the front

A heterogeneous velocity field in the context of a monophasic flow was established in step 2 of the method according to the invention by numerical resolution (with the aid of a monophasic flow simulator) of the equation in pressure.

Then, equation (21) is discretized on the grid discretizing the medium. Then this discretized equation is used to obtain the perturbations of the front for each step of time. An explicit numerical scheme is then used, described below.

The index i is introduced to number the meshes in the direction of the X axis, the index j to number the meshes along the Y axis, and n to number the time steps. We denote Δ x the size of a mesh in the X direction, Δ y the size of a mesh in the Y direction, and Δ t the time step.

1. 1- Tout d'abord on effectue une transformée de Fourier rapide pour les perturbations du front. En effet, l'approximation par un développement en série de perturbation possède une formulation explicite simple dans l'espace de Fourier. Ainsi, selon la méthode, on effectue le calcul dans le domaine où celui-ci est le plus simple, puisque le coût d'une Transformée de Fourier Rapide (FFT) est considéré comme négligeable. Un algorithme de Transformée de Fourier Rapide (FFT) en y permet d'estimer rapidement le premier membre de l'équation (21). $${\mathit{\delta h}}_k^n={\displaystyle \sum _{j=1}^{N_y}}{\mathit{\delta h}}_j^n\exp \left(\frac{2\mathit{\pi i}\left(\mathrm{k}-1\right)\left(n-1\right)}{N_y}\right)$$
   
2. 2- Puis on multiplie les modes de Fourier des fluctuations de la position du front par un module du vecteur d'onde dans l'espace de Fourier, et l'on utilise une transformée de Fourier inverse pour revenir dans l'espace réel. Pour obtenir une approximation de la dérivée de la vitesse selon y, on utilise la différence centrale, bien connue des spécialistes. On obtient donc le schéma explicite suivant : $$\begin{array}{ll}{\mathit{\delta h}}_j^{n+1}& ={\mathit{\delta h}}_j^n+c_{0}A\mathrm{\Delta }t\frac{1}{N_y}\frac{2\mathit{\pi }\left(k-1\right)}{N_y\mathrm{\Delta }y}{\mathit{\delta h}}_k^n\exp \left(\frac{2\mathit{\pi i}\left(k-1\right)\left(n-1\right)}{N_y}\right)+\frac{c_0}{u_0}\mathrm{\Delta }t\left({\mathit{\delta u}}_{i,j}^n-{\mathit{\delta v}}_{i,j}^n\frac{{\mathit{\delta h}}_{j+1}^n-{\mathit{\delta h}}_{j-1}^n}{2\mathrm{\Delta }y}\right)\\ {\mathit{\delta h}}_k^n& ={\displaystyle \sum _{j=1}^{N_y}}{\mathit{\delta h}}_j^n\exp \left(-\frac{2\mathit{\pi i}\left(k-1\right)\left(n-1\right)}{N_y}\right)\end{array}$$
   

où l'index i est défini par la position du front au pas de temps précèdent n. Ici le terme ${\mathit{\delta h}}_{i,k}^n$

représente la transformée de Fourier du front à la fréquence k, δu représente la fluctuation de la vitesse longitudinale selon l'axe X (δu_x ), et δv représente la fluctuation de la vitesse longitudinale selon l'axe Y (δu_y ). La transformation de Fourier discrète et son inverse peuvent être implémentée en utilisant un algorithme de transformée de Fourier rapide (FFT) pour accélérer les calculs.The method of solving the discretized equation (21) is then as follows.

1. 1- Firstly, a fast Fourier transform is carried out for the disturbances of the front. Indeed, the approximation by a serial disturbance development has a simple explicit formulation in the Fourier space. Thus, according to the method, the calculation is carried out in the field where it is the simplest, since the cost of a Fast Fourier Transform (FFT) is considered negligible. A transform Fourier algorithm (FFT) y allows to quickly estimate the first member of equation (21). $${\mathit{.DELTA.h}}_k^{\mathrm{not}}={\displaystyle \underset{j=1}{\overset{{\mathrm{NOT}}_{\mathrm{there}}}{\Sigma }}}{\mathit{.DELTA.h}}_j^{\mathrm{not}}\exp \left(\frac{2\mathit{\pi i}\left(\mathrm{k}-1\right)\left(\mathrm{not}-1\right)}{{\mathrm{NOT}}_{\mathrm{there}}}\right)$$
   
2. 2- Then we multiply the Fourier modes of the fluctuations of the position of the front by a modulus of the wave vector in the Fourier space, and we use an inverse Fourier transform to return to the real space. To get an approximation of the derivative of the velocity according to y , we use the central difference, well known to the specialists. We thus obtain the following explicit scheme: $$\begin{array}{ll}{\mathit{.DELTA.h}}_j^{\mathrm{not}+1}& ={\mathit{.DELTA.h}}_j^{\mathrm{not}}+{\mathrm{vs}}_0\mathrm{AT}\mathrm{\Delta }t\frac{1}{{\mathrm{NOT}}_{\mathrm{there}}}\frac{2\mathit{\pi }\left(k-1\right)}{{\mathrm{NOT}}_{\mathrm{there}}\mathrm{\Delta }\mathrm{there}}{\mathit{.DELTA.h}}_k^{\mathrm{not}}\exp \left(\frac{2\mathit{\pi i}\left(k-1\right)\left(\mathrm{not}-1\right)}{{\mathrm{NOT}}_{\mathrm{there}}}\right)+\frac{{\mathrm{vs}}_0}{u_0}\mathrm{\Delta }t\left({\mathit{\delta u}}_{i,j}^{\mathrm{not}}-{\mathit{.DELTA.V}}_{i,j}^{\mathrm{not}}\frac{{\mathit{.DELTA.h}}_{j+1}^{\mathrm{not}}-{\mathit{.DELTA.h}}_{j-1}^{\mathrm{not}}}{2\mathrm{\Delta }\mathrm{there}}\right)\\ {\mathit{.DELTA.h}}_k^{\mathrm{not}}& ={\displaystyle \underset{j=1}{\overset{{\mathrm{NOT}}_{\mathrm{there}}}{\Sigma }}}{\mathit{.DELTA.h}}_j^{\mathrm{not}}\exp \left(-\frac{2\mathit{\pi i}\left(k-1\right)\left(\mathrm{not}-1\right)}{{\mathrm{NOT}}_{\mathrm{there}}}\right)\end{array}$$
   

where the index i is defined by the position of the front at the preceding time step n . Here the term ${\mathit{.DELTA.h}}_{i,k}^{\mathrm{not}}$

represents the Fourier transform of the front at the frequency k , δ u represents the fluctuation of the longitudinal velocity along the X axis (δ u _x ), and δ v represents the fluctuation of the longitudinal velocity along the Y axis (δ u _y ). The discrete Fourier transform and its inverse can be implemented using a Fast Fourier Transform (FFT) algorithm to accelerate computations.

Having calculated the fluctuations of the front for a given time step N _t , we reconstruct the position of the front in the discretization grid of the medium using the following formula: $${\mathit{h}}_j^{{\mathrm{NOT}}_t}=\mathrm{\Delta }t{\mathrm{NOT}}_t{\mathrm{vs}}_0+\mathit{\delta }{\mathit{h}}_j^{{\mathrm{NOT}}_t}$$

It is recalled that the approximation of the velocity u was carried out once and for all, since the viscous coupling is modeled via the convolution term (first term of the equation (21)).

It is therefore possible to determine explicitly the position of the front for each time step, after having determined once, at time t = 0, the velocities of the fluids by solving the pressure equation with the help of a flow simulator.

It should be noted that the proposed technique requires only one complete numerical resolution of the pressure equation, unlike the market techniques that require as much as there is no time. A comparison between the method according to the invention ( Figure 1A ) and a simulation using a market-current simulator ( Figure 1B ) illustrates that despite the approximation made to save computing time, the precision of the results remain comparable to other methods. The Figure 1A shows the simulated interface front with the method, while the Figure 1B shows the front of the simulated interface using a market-current simulator. These figures represent a section of the subsoil, the x-axis corresponds to a horizontal coordinate x , while the y-axis represents the depth y.

Optimization of oil recovery by injection

From the position of the front, it is known to determine the saturations in each mesh of the medium and for each time step. One can for example use the current tube method (well known to specialists) far from the front, and interpolate saturations near the front whose position is now well determined.

By following the geographical and temporal evolution of the saturations one can determine when and where the fluid injected to facilitate the recovery will reach the well. When water is used to "push" the oil in place of a reservoir, the knowledge of the saturations makes it possible to determine the time of breakthrough of the water, which a key data of the exploitation of a oil field.

1. 1 D'échantillonner rapidement l'espace des paramètres d'étude sur lesquels porte la modélisation, de façon à optimiser un ou plusieurs critères économiques d'exploitation (décision d'exploitation, choix de la position des puits, d'un procédé de récupération etc.).
2. 2 De modifier de façon cohérente le modèle géologique afin de caler au mieux les données de production observées, y compris des données de sismiques répétées. Ceci peut permettre de forer en s'adaptant aux hétérogénéités géologiques, de localiser les fluides, et donc de mieux contrôler le scénario de récupération.
3. 3 D'estimer des incertitudes en effectuant des simulations de Monte Carlo, de façon à tester le rôle d'hétérogénéités, ou de paramètres mal connus caractérisant le sous-sol. Ceci permet de quantifier le niveau de risque lié à l'exploitation d'un réservoir et donc de redéfinir des paramètres d'exploitation voire les paramètres économiques.

Due to a rapid determination of saturations, those skilled in the art are able to:

1. 1 To rapidly sample the space of the study parameters to which the modeling relates, so as to optimize one or more economic criteria of exploitation (exploitation decision, choice of the position of the wells, a recovery process etc.).
2. 2 To consistently modify the geological model to best fit the observed production data, including repeated seismic data. This can make it possible to drill by adapting to the geological heterogeneities, to locate the fluids, and thus to better control the recovery scenario.
3. 3 To estimate uncertainties by carrying out Monte Carlo simulations, in order to test the role of heterogeneities, or poorly known parameters characterizing the subsoil. This makes it possible to quantify the level of risk associated with the operation of a reservoir and thus to redefine operating parameters or even the economic parameters.

u _nw :

Vecteur vitesse du fluide non-mouillant

u _w :

Vecteur vitesse du fluide mouillant

u :

Vitesse totale des fluides.

λ _nw :

Mobilité du fluide non-mouillant

λ _w :

Mobilité du fluide mouillant

p_nw :

Pression du fluide non-mouillant

p_w :

Pression du fluide mouillant

K :

Tenseur perméabilité absolue du milieu

k _rnw :

Perméabilité relative du fluide non-mouillant

k _rw :

Perméabilité relative du fluide mouillant

µ _nw :

Viscosité du fluide non-mouillant

µ _w :

Viscosité du fluide mouillant

p_c (S) :

Pression capillaire

S_nw :

Saturation du fluide non-mouillant

S_w :

Saturation du fluide mouillant

p :

Pression commune des deux fluides, si p_c (S) = 0.

φ :

Porosité (supposée uniforme dans le réservoir)

t :

Temps

S _wi :

Saturation irréductible en eau

S_f :

Saturation en eau au front

M_f :

Rapport de mobilité totale au front.

q :

vecteur d'onde dans l'espace de Fourier

c₀ :

Vitesse du front non perturbé

u ₀ :

Vitesse constante de référence

S_wr :

Saturation maximale en eau

f_w' :

Dérivée du flux fractionnaire en eau

h ₀ :

Position de référence du front

h(y, t) :

Position du front

δh(y, t) :

Fluctuation de la position du front.

n :

Vecteur normal au front

u ₀ :

Vitesse totale moyenne.

δu _nw :

Perturbation du vecteur vitesse du fluide non-mouillant

δu _w :

Perturbation du vecteur vitesse du fluide mouillant

δu :

Perturbation du vecteur vitesse totale des fluides.

δu_x :

Perturbation de la composante x de la vitesse totale.

δu_y :

Perturbation de la composante y de la vitesse totale.

F :

Équation du front.

F_l :

Équation du front sous forme canonique

L_x, L_y :

Taille en x et en y du domaine.

V_n :

Vitesse normale au front

f(S) :

Flux fractionnaire en eau

l₁ :

Aire d'injection

Q :

Débit total injecté.

V :

Volume du milieu poreux

x :

Coordonnée longitudinale.

y :

Coordonnée transverse.

k :

Vecteur d'onde discrétisé

t_D :

Temps adimensionné par le temps de percée moyen.

An exhaustive study incorporating these various aspects (sensitivity study, data calibration and estimation of uncertainties) potentially requires numerous calls to a flow simulator (dedicated software). As a result, the existing solutions do not make it possible to obtain quantitative answers to practical questions that the petroleum engineer asks himself when tank simulation studies are carried out. Hence the interest of the method according to the invention which proposes a rapid calculation tool, of precision compatible with the accuracy of the input data.
The following notations are used during the description of the invention:

u _nw :

Vector velocity of non-wetting fluid

u _w :

Vector speed of the wetting fluid

u :

Total speed of the fluids.

λ _nw :

Mobility of the non-wetting fluid

λ _w :

Mobility of the wetting fluid

p _nw :

Non-wetting fluid pressure

p _w :

Wetting fluid pressure

K :

Tensor absolute permeability of the medium

k _rnw :

Relative permeability of the non-wetting fluid

k _rw :

Relative permeability of the wetting fluid

μ _nw :

Viscosity of the non-wetting fluid

μ _w :

Viscosity of the wetting fluid

p _c ( S ):

Hair pressure

S _nw :

Saturation of the non-wetting fluid

S _w :

Saturation of the wetting fluid

p:

Common pressure of the two fluids, if p _c ( S ) = 0.

φ:

Porosity (supposedly uniform in the tank)

t :

Time

S _wi :

Irreducible saturation in water

S _f :

Water saturation at the front

M _f :

Total mobility ratio at the front.

q:

wave vector in Fourier space

c ₀ :

Front speed undisturbed

u ₀ :

Constant reference speed

S _wr :

Maximum water saturation

f _w ' :

Derived from the fractional flow of water

h ₀ :

Reference position of the front

h (y, t) :

Front position

δ h (y, t) :

Fluctuation of the position of the forehead.

n :

Normal vector at the front

u ₀ :

Average total speed.

δ u _nw :

Disturbance of the velocity vector of the non-wetting fluid

δ u _w :

Disturbance of the speed vector of the wetting fluid

δ u :

Disturbance of the vector total velocity of the fluids.

δ u _x :

Disturbance of the x component of the total speed.

δ u _y :

Disturbance of the y component of the total speed.

F :

Forehead equation.

F _l :

Equation of the front in canonical form

L _x , L _y :

Size in x and y of the domain.

V _n :

Normal speed at the front

f (S) :

Fractional flow of water

l ₁ :

Injection area

Q :

Total flow injected.

V :

Volume of the porous medium

x:

Longitudinal coordinate.

y:

Transverse coordinate.

k:

Discretized wave vector

t _D :

Time scaled by average breakthrough time.

Claims (6)

1. A method of optimizing the recovery of a fluid in place in a heterogeneous porous medium, by injecting into said medium a sweeping fluid intended to cause the fluid in place to flow, said fluids being immiscible, wherein said medium is discretized into a grid consisting of a set of cells, the method being characterized in that it comprises the following stages:

   - determining a velocity and a direction of flow for at least one of said fluids, using a flow simulator to solve a pressure equation,

   - defining a relation describing a position of a front separating said fluids, by disregarding a viscous coupling using the perturbation theory, and by taking into account the velocity fluctuations in the front advance direction and the velocity fluctuations in the direction perpendicular to the front advance direction, and for different time intervals, and wherein said relation comprises a first term representing a viscous coupling describing the perturbation due to the saturation variation, and a second term describing the front position perturbations caused by the medium heterogeneities,

   - reconstructing the position of the front in said grid by means of a discretization of said relation and of a fast Fourier transform ; and

   - optimizing the injection of said sweeping fluid according to the position of said front.
2. A method as claimed in claim 1, wherein said first term is obtained by considering said homogeneous medium and the Buckley-Leverett type displacements of said front.
3. A method as claimed in claim 1, wherein said second term is obtained using the fact that the front is a material surface, and by representing the velocity as a sum of a mean velocity with velocity fluctuations due to the heterogeneities.
4. A method as claimed in claim 2, wherein said first term depends on at least the following parameters: a frontal mobility ratio, a mean filtration rate along a front advance direction, a porosity of the medium, a Buckley-Leverett function representing a fractional water flow, a water saturation at the front, a maximum water saturation, the wave vector in the Fourier space.
5. A method as claimed in claim 3, wherein said second term depends on at least the following parameters: a nonperturbed front velocity, a mean total velocity, perturbations of the components of the total velocity in the front advance direction and in the direction perpendicular thereto.
6. A method as claimed in any one of the previous claims, wherein injection is optimized by determining at each cell and for different time intervals the saturation of at least one of said fluids from the position of said front.

Images