EP1889999B1

EP1889999B1 - Méthode pour optimiser la récupération assistée d'un fluide en place dans un milieu poreux par suivi de front

Info

Publication number: EP1889999B1
Application number: EP07290895A
Authority: EP
Inventors: Benoît Noetinger; Pavel Spesivtsev
Original assignee: IFP Energies Nouvelles IFPEN
Current assignee: IFP Energies Nouvelles IFPEN
Priority date: 2006-08-16
Filing date: 2007-07-16
Publication date: 2011-09-14
Anticipated expiration: 2027-07-16
Also published as: FR2904982B1; US20080046223A1; FR2904982A1; US7912686B2; NO20073881L; EP1889999A1

Description

La présente invention concerne une méthode pour optimiser l'exploitation d'un milieu poreux hétérogène dans le cadre de récupération assistée d'un fluide en place, par détermination de la position du front séparant un fluide de balayage et le fluide en place.
En particulier la méthode selon l'invention permet de déterminer la position d'un front séparant deux fluides immiscibles en mouvement, sans remise à jour systématique du champ de pression, de façon à obtenir une simulation des écoulements polyphasiques dans le milieu poreux suffisamment rapide et précise pour obtenir des informations quantitatives permettant une exploitation optimale du milieu.
La méthode trouve des applications notamment pour l'exploitation de gisements de pétrole ou de gaz, ou l'exploitation de stockages souterrains de gaz par exemple.

Présentation de l'art antérieur

Toutes les notations utilisées pour décrire l'art antérieur et l'invention sont définies en fin de description.
Savoir décrire et simuler les écoulements polyphasiques dans les réservoirs souterrains est à la base du métier des ingénieurs de gisement dans les compagnies pétrolières ou gazières (et de façon analogue dans des compagnies d'adduction d'eau). La présence des hétérogénéités du sous-sol, et la complexité croissante des systèmes de drainage rendent impossible toute solution analytique simple, ce qui impose le développement de solutions numériques utilisant un modèle maillé. La simulation des écoulements polyphasiques en milieu poreux hétérogène peut alors requérir d'importantes ressources informatiques, en particulier lorsque le modèle numérique du milieu considéré est fortement détaillé. C'est le cas notamment en ingénierie de réservoirs, dans le domaine pétrolier. Ce coût provient principalement de la résolution des systèmes linéaires de grande taille issus de l'équation gouvernant la pression, qu'il convient de réactualiser en suivant le déplacement des fluides, si l'on désire disposer d'une solution possédant une bonne précision.
Dans le cadre d'un écoulement diphasique dans un milieu poreux hétérogène, on considère le déplacement d'un fluide en place (par exemple l'huile) sous l'effet de l'injection d'un autre fluide (par exemple de l'eau). La loi de Darcy généralisée gouvernant le mouvement des fluides s'écrit, moyennant des hypothèses et notations standard, de la façon suivante : $u_{nw} = - λ_{nw} \nabla p_{nw}, λ_{nw} = \frac{K k_{rnw}}{μ_{nw}}$
$u_{w} = - λ_{w} \nabla p_{w}, λ_{nw} = \frac{K k_{rw}}{μ_{w}}$
$u = u_{nw} + 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_{nw} + S_{w} = 1$
De ce fait, on peut écrire toutes les fonctions dépendant de la saturation S en terme de la saturation S_w .
En négligeant la pression capillaire, on a $p_{nw} = p_{w} = p$
La vitesse totale prend alors la forme $u = - λ (S_{w}) \nabla p$

avec $\nabla\cdot u = 0$
Ce système d'équations (6) et (7) définit l'équation en pression.
Dans les équations (1) et (2), les fonctions k_r sont les perméabilités relatives. On a k_r = k_r (S_w ). Les fonctions λ _nw , λ _w sont les mobilités respectives des fluides, et λ(S_w ) = λ _nw (S_w ) + λ _w (S_w ) est la mobilité totale, qui dépend donc explicitement de S_w .
Les lois de conservation de la masse relatives à chacun des fluides s'écrivent $ϕ \frac{\partial S_{nw}}{\partial t} + \nabla\cdot u_{nw} = 0$
$ϕ \frac{\partial S_{w}}{\partial t} + \nabla\cdot u_{w} = 0$
Ce système d'équations (8) et (9) définit l'équation en saturation. Ici φ est la porosité (supposée uniforme dans le réservoir), t le temps, et u _n,nw est la vitesse du fluide considéré.
Une des difficultés pour la résolution de ces équations vient du couplage entre les équations (6) à (9) qui couplent la saturation au champ de pression. Ces effets sont bien connus et pilotent le développement d'éventuelles instabilités visqueuses. En particulier, dans le cas 1D, en considérant le balayage d'un milieu initialement saturé en huile par de l'eau, ce système d'équations peut être résolu par la méthode des caractéristiques, dont les solutions sont caractérisées par l'existence d'oscillations, correspondant à des saut de saturation se propageant dans le milieu.
La discontinuité est caractérisée par une saturation en eau passant de S_wi à S_f, la saturation au front dont la valeur est donnée par une condition de type Rankine-Hugoniot. Dans le contexte pétrolier, la construction de la tangente de Welge permet de déterminer S_f graphiquement. En arrière de ce front, une "onde de raréfaction" assure la transition entre S_f et la saturation en eau à l'injection.
On connaît plusieurs techniques pour déterminer la simulation des écoulements polyphasiques en milieu poreux hétérogène. Ces méthodes reposent toutes sur la résolution du système d'équations suivant : ${\begin{cases} \begin{array}{l} Équation en pression : \\ u = - λ (S_{w}) \nabla p (10) \\ \nabla\cdot u = 0 (11) \end{array} \\ \begin{array}{l} Équation en saturation : \\ ϕ \frac{\partial S_{nw}}{\partial t} + \nabla\cdot u_{nw} = 0 (12) \\ ϕ \frac{\partial S_{w}}{\partial t} + \nabla\cdot u_{w} = 0 (13) \end{array} \end{cases}$
Les techniques se différencient par la méthode utilisée pour résoudre ce système. Dans le cas général, en milieu hétérogène complexe, ce type de système d'équations doit être résolu numériquement dont une référence est : B. Noetinger, P.E. 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.
Les compagnies pétrolières utilisent le plus souvent les techniques de discrétisation classiques, de type volumes finis (Aziz, Kb., Settary, A. : Petroleum réservoir simulation. Applied science publishers, London, 1979). Concernant la discrétisation temporelle, les choix sont multiples : si on privilégie la robustesse numérique, on peut être totalement implicite en pression comme en saturation. Si on préfère éviter des oscillations de la solution, et une meilleure précision, on choisira un schéma implicite en pression, explicite en saturation. On constate en général l'existence de fronts caractérisés par la même discontinuité de la saturation dont la valeur passe brutalement de la saturation irréductible en eau, S_wi , à la saturation du front, notée S_f . Ces discontinuités sont dues au caractère hyperbolique de l'équation de transport en saturation, et sont bien décrites analytiquement dans le cas du déplacement 1D : c'est la théorie de Buckley-Leverett dont on trouvera une description détaillée dans l'ouvrage de Charles-Michel Marle. Les écoulements polyphasiques en milieu poreux. Seconde édition revue et augmentée 1972. Editions Technip, Paris.
Quelque soit la méthode retenue, des remises à jour régulières de l'équation en pression sont nécessaires, de façon à estimer les vitesses avec précision.
Dans une autre classe de techniques, qui deviennent de plus en plus employées dans les compagnies opératrices, on trouve les techniques dites de lignes de courant, dont une référence est : R. P. Batycky, M.J. Blunt, and M.R. Thiele, A 3d field-scale streamline-based reservoir simulator, SPERE, 11, 246-254 (1997).
Dans ces méthodes les saturations sont mises à jour en suivant le déplacement des fluides dans leur mouvement le long des lignes de courant, courbes paramétrées définies par dx/dt = u(x,t). Un changement de variable permet de se ramener à la théorie de Buckley-Leverett précédemment évoquée. Toutefois, afin de bien calculer la vitesse u(x, t) la pression est résolue autant de fois que la précision l'exige, sur une grille cartésienne. Ces méthodes sont plus rapides que les méthodes classiques. Elles n'en possèdent pas non plus la versatilité et la robustesse pour traiter des problèmes compositionnels complexes, où la méthode traditionnelle reste la seule issue possible. Les mises à jour onéreuses du champ de pression ne sont pas évitées.
Enfin, il existe les techniques de suivi de front (R Juanes a,d KA Lie "A front tracking method for efficient simulation of miscible gas injection." paper SPE 92298, Houston 31 jan. to 2 feb. 2005), dont l'idée est de suivre le front par une approche lagrangienne, qui nécessite les mêmes mises à jour de la pression. Dans ce type de méthode, la description des écoulements polyphasiques dans les réservoirs souterrains consiste à déterminer la position du front (appelé aussi interface) séparant deux fluides immiscibles en mouvement : un fluide en place (de l'huile par exemple) et un fluide de balayage (de l'eau par exemple), appelé également fluide injecté
L'évolution du front dans le réservoir au cours de l'écoulement est considérablement influencée par le couplage (appelé couplage visqueux par les ingénieurs réservoirs) entre le champ de pression et le champ de saturation (Saffman P.G. 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-oil front for two-phase, immiscible flows in heterogeneous porous media. 2-isotropic media. Transport in porous media, 56:305-328, 2004).
En particulier, lorsque le fluide injecté est moins visqueux et par conséquent plus mobile au niveau du front que le fluide en place, les instabilités visqueuses favoriseront toujours l'écoulement des fluides dans les couches les plus perméables du milieu. Le temps de percée au travers de celles-ci est beaucoup plus rapide que dans le reste du réservoir. Au contraire, si le fluide injecté est moins mobile, le couplage visqueux peut le ralentir dans les couches initialement plus rapides, compensant ainsi les différences de perméabilité dues à la stratification. Un front stationnaire apparaît alors (Artus V., B. Noetinger, and L. Ricard. Dynamics of the water-oil front for two-phase, immiscible flows in heterogeneous porous media. 1-stratified media. Transport in Porous Media, 56:283-303, 2004).
Il est connu que la stabilité du front suite à de petites perturbations est déterminée par le rapport des viscosités des fluides. Si le fluide le moins visqueux est déplacé par le fluide le plus visqueux, alors le front est stable, et vice versa. Si la situation est instable, alors des phénomènes, appelés « doigts visqueux », se développent avec le temps.
Les effets de la perméabilité relative sont pris en compte en considérant les déplacements de type Bucley-Leverett dans un milieu homogène. De tels déplacements sont décrits par une équation en saturation hyperbolique. Une solution à cette équation par la méthode des caractéristiques fournis des solutions qui ont des oscillations caractérisées par une condition de Rankine Hugoniot. Il a été montré que la stabilité du front où la saturation a un « choc » est déterminée par un rapport de mobilité frontal, noté M_f : si M_f < 1, alors le front est stable au regard des petites perturbations, et vice versa. Ce rapport de mobilité frontal peut être rapproché du rapport de viscosité en utilisant la fonction de perméabilité relative.
Plus récemment, ce résultat a été confirmé dans le cadre d'analyses de stabilité linéarisées du procédé d'injection d'eau. Dans ce cadre, il a été montré que l'évolution de la transformée de Fourier des fluctuations de la position du front, fluctuations notées δh(q, t), dans un milieu homogène est donnée par : $\partial_{t} δh (q t) = c_{0} |q| \frac{M_{f} - 1}{M_{f} + 1} δh (q t)$

où $c_{0} = \frac{u_{0}}{ϕ} \frac{f_{w} (S_{f}) - f_{w} (S_{wr})}{S_{f} - S_{wr}} = \frac{u_{0}}{ϕ} f_{w}^{ʹ} (S_{f})$
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

Pour résumer, toutes les techniques numériques précitées ont un coût en temps de calcul très important qui provient principalement du grand nombre de résolutions des systèmes linéaires de grande taille issus de l'équation gouvernant la pression, qu'il convient de réactualiser en suivant le déplacement des fluides avec une bonne précision.
De ce fait, les solutions existantes ne permettent pas d'obtenir des réponses quantitatives à des questions pratiques que l'ingénieur pétrolier se pose dès lors que l'on procède à des études de simulation de réservoir, à savoir :

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

Une étude exhaustive incorporant ces divers aspects (étude de sensibilité, calage aux données et estimation d'incertitudes) peut alors nécessiter potentiellement de nombreux appels à un simulateur d'écoulement (logiciel dédié), d'où l'intérêt de disposer d'outils de calcul rapides, de précision compatible avec la précision des données d'entrée.
La méthode selon l'invention permet de déterminer la position d'un front séparant deux fluides immiscibles en mouvement en milieu poreux hétérogène, sans remise à jour systématique du champ de pression de façon à obtenir une simulation des écoulements polyphasiques dans le milieu poreux suffisamment rapide pour obtenir des informations quantitatives permettant une exploitation optimale du réservoir, par la détermination de paramètres techniques d'exploitation et/ou de paramètres économiques.

La méthode selon l'invention

L'invention concerne une méthode pour optimiser l'exploitation d'un milieu poreux hétérogène contenant un fluide, dit fluide en place. L'exploitation est réalisée par injection dans le milieu d'un second fluide, dit fluide de balayage, pour provoquer un écoulement du fluide en place. Ces deux fluides sont immiscibles. Selon la méthode, on discrétise le milieu en une grille constituée d'un ensemble de mailles. La méthode comporte les étapes suivantes :

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.

Selon la méthode, la relation définie peut comporter notamment les deux termes suivants :

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.

Présentation sommaire des figures

Les figures 1A et 1B montrent une comparaison entre le résultat fourni par la méthode selon l'invention (fig. 1A), et une simulation par lignes de courant (fig. 1B).

Description détaillée de la méthode

La méthode selon l'invention permet d'optimiser l'exploitation d'un milieu poreux hétérogène dans le cadre de récupération assistée de pétrole (EOR) par exemple. Cette technique consiste à injecter, dans un réservoir pétrolier, un fluide de balayage (CO₂, eau) pour provoquer un écoulement d'un fluide en place (huile) contenu dans le réservoir et que l'on cherche à extraire. Ces deux fluides doivent bien évidemment être immiscibles.
Dans ce cadre de récupération assistée, on considère une situation où il existe une interface, appelée front, séparant ces deux fluides. A l'échelle du réservoir la largeur de ce front est petite comparée à la taille du réservoir pétrolier étudié, et par souci de simplicité on suppose que ce front est étroit relativement à la taille du réservoir. On focalise l'étude sur la dynamique de ce front dans un milieu poreux hétérogène, en négligeant les effets capillaires et la gravité.
La méthode selon l'invention permet d'optimiser l'injection en déterminant, pour différents instants, appelés pas de temps, la position du front séparant les deux fluides en mouvement en milieu poreux hétérogène (le réservoir par exemple). La méthode s'appuie sur une technique originale de remise à jour systématique du champ de pression qui procure l'avantage d'être à la fois précise et rapide, dans la mesure où la méthode ne nécessite qu'un appel très limité à un simulateur effectuant des calculs à la fois longs et complexes. Pour ce faire, la mise à jour du champ de pression, due au changement de géométrie du front, est effectuée de façon semi analytique à l'aide de techniques de Transformées de Fourier Rapide (FFT). On économise ainsi l'essentiel du coût de calcul dû à ces mises à jour du champ de pression, d'où la rapidité de la méthode. La méthode comporte principalement les étapes suivantes :

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 ;

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.

Discrétisation du milieu

La discrétisation du milieu est une étape classique et bien connue des spécialistes. Elle permet de représenter la structure du milieu en un ensemble de mailles, caractérisées par leurs dimensions et leurs positions géographiques. Les calculs sont effectués pour chaque maille. Cette étape est indispensable pour utiliser des simulateurs d'écoulement (le simulateur est un logiciel dédié). Cet ensemble de maille est appelé « grille ». Cette grille, discrétisant le milieu peut être rectangulaire de taille L_x × L_y mailles.

Détermination d'un champ de vitesse hétérogène

Au cours de cette étape, on détermine un champ de vitesse hétérogène dans le cadre d'un écoulement monophasique. Pour ce faire on résout numériquement l'équation en pression (équations (6) et (7)), qui est obtenue en combinant la loi de Darcy et l'équation d'incompressibilité à l'aide d'un simulateur d'écoulement monophasique. On peut déterminer ainsi le vecteur vitesse de chacun des fluides, c'est-à-dire, leur vitesse et leur direction d'écoulement, à un instant t = 0.
Pour un domaine 2D (le cas 3D se traite exactement de la même façon.) de taille L_x × L_y dont la carte de perméabilité est K(x, y) et pour lequel l'écoulement principal est parallèle à un axe X, y repérant la coordonnée transverse (sur l'axe Y), la résolution de l'équation en pression consiste à résoudre l'équation suivante : $\nabla\cdot (λ (S_{w} (x, y, t = 0) \cdot \nabla P (x y)) = 0$
On déduit donc de ces équations le champ de vitesse total u : u _nw , le vecteur vitesse du fluide non-mouillant, et u _w , le vecteur vitesse du fluide mouillant.
Cette résolution, qui fait appel à un simulateur, est effectuée une seule fois selon la méthode. L'équation en pression est résolue au début, c'est-à-dire pour t = 0, elle ne l'est plus aux pas de temps ultérieurs (t>0).
On peut utiliser une discrétisation de type volumes finis cinq points standard pour résoudre cette équation à l'aide du simulateur : la pression est estimée au centre des mailles, les flux entre mailles faisant intervenir des transmissivités calculées comme moyenne harmonique des perméabilités des deux mailles. La loi de Darcy permet alors de déterminer les flux et les vitesses.
La carte de perméabilité K(x, y) peut être générée à l'aide d'un générateur de champs aléatoire, comme décrit dans le document suivant par exemple :

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.

Définition d'une relation décrivant la position du front en milieu hétérogène

L'idée de base consiste donc à éviter de recalculer à chaque itération en temps l'écoulement dans tout le domaine de calcul, en se restreignant à la description du mouvement du front, décrit par la fonction δh(y, t). Ceci est rendu possible par une approximation permettant de relier la modification du champ de pression δp(x, y, t), et donc des vitesses, à la distorsion du front δh(y, t). On peut ainsi éliminer la pression du problème et on est alors amenés à la résolution d'un système différentiel impliquant δh(y, t).
On a vu que dans un milieu homogène, l'évolution de la transformée de Fourier des fluctuations de la position du front, fluctuations notées δh(q, t), est donnée par : $\partial_{t} δh (q t) = c_{0} |q| \frac{M_{f} - 1}{M_{f} + 1} δh (q t)$

où $c_{0} = \frac{u_{0}}{ϕ} \frac{f_{w} (S_{f}) - f_{w} (S_{wr})}{S_{f} - S_{wr}} = \frac{u_{0}}{ϕ} f_{w}^{ʹ} (S_{f})$

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

En milieu hétérogène, le problème est plus complexe. Les hétérogénéités du champ de perméabilité initient des perturbations du front, c'est-à-dire que le front se déforme en fonction de la perméabilité des strates. Ainsi, ces perturbations peuvent grossir ou diminuer en fonction du rapport des viscosités. C'est le problème du couplage visqueux. Le problème est non linéaire et une analyse mathématique s'avère très compliquée.
Selon l'invention, on propose de tenir compte des hétérogénéités, d'une part en s'affranchissant du couplage visqueux à l'aide de la théorie des perturbations, et d'autre part en prenant en compte les fluctuations transversales de la vitesse, c'est-à-dire les fluctuations dans un axe Y perpendiculaire à l'axe X d'avancée du front.

Découplage par la technique des perturbations

La méthode selon l'invention estime le couplage entre les perturbations de vitesses δu et la fluctuation de la position du front δh de façon semi analytique. En ce sens, la méthode utilise la théorie des perturbations.
D'un point de vue heuristique, la théorie des perturbations est une méthode mathématique générale qui permet de trouver une solution approchée d'une équation mathématique (E _λ) dépendante d'un paramètre λ lorsque la solution de l'équation (E ₀), correspondant à la valeur λ = 0, est connue exactement. L'équation mathématique (E _λ) peut être une équation algébrique, une équation différentielle, une équation aux valeurs propres, ... La méthode consiste à chercher la solution approchée de l'équation (E _λ) sous la forme d'un développement en série des puissances du paramètre λ, cette solution approchée étant supposée être une approximation d'autant meilleure de la solution exacte, mais inconnue, que la valeur absolue du paramètre λ est plus « petite ».
Ainsi, selon l'invention, la théorie des perturbations conduit à considérer localement la perturbation de vitesse comme une somme de deux perturbations : $δ u (x y t) = δ u_{a} (x y t) + δ u_{b} (x y t)$
Le premier terme, δu _a décrit la perturbation due à la variation de saturation, et le second terme δu _b décrit la perturbation due à l'hétérogénéité du milieu.
Selon cette décomposition, on peut décrire la position du front au cours de son avancée dans un milieu poreux hétérogène par l'équation suivante : $\partial_{t} δh (q t) = c_{0} A |q| δh (q t) + \frac{c_{0}}{u_{0}} {δu}_{bx} (x_{f} q t)$
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

Ainsi, à l'aide de cette formulation de la position du front, il est possible de découpler le problème. De cette façon, pour des champs de perméabilité à petites variances, la stabilité du front subissant de petites perturbations est déterminée par le rapport de mobilité frontal.

Prise en compte des fluctuations transversales de la vitesse

Selon l'invention, on simule la propagation du front dans un milieu hétérogène, c'est-à-dire dans un champ de vitesse hétérogène, en modélisant l'influence des effets visqueux de façon analytique. Le champ de vitesse hétérogène a été obtenu en utilisant une simulation d'écoulement monophasique (étape 2 de la méthode).
Pour simplifier les explications, on se place en deux dimensions. Le front séparant les deux fluides immiscibles est alors décrit par une fonction 2D, x(y, t) = h(y, t). On peut écrire cette équation sous la forme : $F (x y t) = 0$
En la résolvant par rapport à x, on a : $F_{1} (x y t) = h (y t) - x = 0$
Puis, en utilisant le fait que la surface en question est matérielle, i.e. composée d'éléments de fluide qui suivent le mouvement local du fluide, et en utilisant l'expression pour la vitesse de propagation du front pour les déplacements de type Buckley-Leverett, on obtient : $\partial_{t} h (y t) = c_{0} (u \cdot \nabla φ) |_{φ = 0}$
En considérant maintenant l'hétérogénéité du champ de vitesse, on représente la vitesse comme une somme de deux termes. Le premier représente la moyenne, et le second les fluctuation dues aux hétérogénéités : $u (r) = u_{0} + δ u (r)$
Avec : r = {x,y}, u ₀ = {u ₀,0} et δu={δu_x ,δu_y }
Les perturbations du champ de vitesse provoquent les perturbations du front h(y,t) = h ₀(t)+δh(y,t). En substituant cette décomposition dans l'équation (19) et en extrayant le terme moyen en utilisant la relation ∂ _th ₀(t) = c ₀, qui est valable pour les déplacement de type Buckley-Leverett, on obtient l'équation suivante, décrivant les perturbations de la position du front provoquées par les hétérogénéités du milieu : $\partial_{t} h (y t) = \frac{c_{0}}{u_{0}} (δ u_{x} (x_{f} y) - δ u_{y} (x_{f} y) \partial_{y} δh (y t))$
Ainsi, par application de la théorie des perturbations, de la théorie des surfaces matérielles, et en utilisant l'expression de la vitesse de propagation du front pour les déplacement de type Buckley-Leverett, on peut décrire la position du front dans un milieu poreux hétérogène de la façon suivante (δh a été substituée par h dans cette expression par soucis de clarté) : $\partial_{t} h (y t) = c_{0} A \int \frac{dq}{2 π} |q| h (q t) e^{iqy} + \frac{c_{0}}{u_{0}} (δ u_{x} (x_{f} y) - δ u_{y} (x_{f} y) \partial_{y} h (y t))$
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).
L'étape suivante de la méthode consiste à proposer une technique de résolution de l'équation du front complet.

Estimation de la position du front

Un champ de vitesse hétérogène dans le cadre d'un écoulement monophasique, a été établit à l'étape 2 de la méthode selon l'invention par résolution numérique (à l'aide d'un simulateur d'écoulement monophasique) de l'équation en pression.
Ensuite, l'équation (21) est discrétisée sur la grille discrétisant le milieu. Puis cette équation discrétisée est utilisée pour obtenir les perturbations du front pour chaque pas de temps. On utilise alors un schéma numérique explicite, décrit ci-après.
On introduit l'indexe i pour numéroter les mailles selon la direction de l'axe X, l'indexe j pour numéroter les mailles selon l'axe Y, et n pour numéroter les pas de temps. On note Δx la taille d'une maille selon la direction X, Δy la taille d'une maille selon la direction Y, et Δt le pas de temps.
La méthode de résolution de l'équation (21) discrétisée se présente alors comme suit.

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). ${δh}_{k}^{n} = \sum_{j = 1}^{N_{y}} {δh}_{j}^{n} \exp (\frac{2 πi (k - 1) (n - 1)}{N_{y}})$
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}{l} {δh}_{j}^{n + 1} & = {δh}_{j}^{n} + c_{0} A Δ t \frac{1}{N_{y}} \frac{2 π (k - 1)}{N_{y} Δ y} {δh}_{k}^{n} \exp (\frac{2 πi (k - 1) (n - 1)}{N_{y}}) + \frac{c_{0}}{u_{0}} Δ t ({δu}_{i, j}^{n} - {δv}_{i, j}^{n} \frac{{δh}_{j + 1}^{n} - {δh}_{j - 1}^{n}}{2 Δ y}) \\ {δh}_{k}^{n} & = \sum_{j = 1}^{N_{y}} {δh}_{j}^{n} \exp (- \frac{2 πi (k - 1) (n - 1)}{N_{y}}) \end{array}$

i

n

{δh}_{i, k}^{n}

k

u

X

u_x

v

Y

u_y

Ayant calculé les fluctuations du front pour un pas de temps donné N_t , on reconstruit la position du front dans la grille de discrétisation du milieu à l'aide de la formule suivante : $h_{j}^{N_{t}} = Δ t N_{t} c_{0} + δ h_{j}^{N_{t}}$
On rappelle que l'approximation de la vitesse u a été réalisée une fois pour toute, puisque le couplage visqueux est modélisé via le terme de convolution (premier terme de l'équation (21)).
On peut donc déterminer de façon explicite la position du front pour chaque pas de temps, après avoir déterminer une seule fois, à l'instant t = 0, les vitesses des fluides par résolution de l'équation en pression à l'aide d'un simulateur d'écoulement.
On notera que la technique proposée n'a besoin que d'une seule résolution numérique complète de l'équation en pression, contrairement aux techniques du marché qui en demande autant qu'il y a de pas de temps. Une comparaison entre la méthode selon l'invention (figure 1A) et une simulation utilisant un simulateur à ligne de courant du marché (figure 1B) illustre que malgré l'approximation réalisée pour gagner en temps de calcul, la précision des résultats restent comparables aux autres méthodes. La figure 1A montre le front de l'interface simulé avec la méthode, tandis que la figure 1B montre le front de l'interface simulé à l'aide d'un simulateur à ligne de courant du marché. Ces figures représentent une coupe du sous-sol, l'axe des abscisses correspond à une coordonnée géographique horizontale x, tandis que l'axe des ordonnées représente la profondeur y.

Optimisation de la récupération d'huile par injection

A partir de la position du front, il est connu de déterminer les saturations en chaque maille du milieu et pour chaque pas de temps. On peut par exemple utiliser la méthode des tubes de courant (bien connue des spécialistes) loin du front, et interpoler les saturations près du front dont la position est maintenant bien déterminée.
En suivant l'évolution géographique et temporelle des saturations on peut déterminer à quel moment et où le fluide injecté pour faciliter la récupération va atteindre le puits. Lorsque l'on utilise de l'eau pour « pousser » l'huile en place d'un réservoir, la connaissance des saturations permet de déterminer le temps de percée de l'eau, qui une donnée clef de l'exploitation d'un champ pétrolier.
Du fait d'une détermination rapide des saturations, l'homme du métier est en mesure :

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

Une étude exhaustive incorporant ces divers aspects (étude de sensibilité, calage aux données et estimation d'incertitudes) nécessite potentiellement de nombreux appels à un simulateur d'écoulement (logiciel dédié). De ce fait, les solutions existantes ne permettent pas d'obtenir des réponses quantitatives à des questions pratiques que l'ingénieur pétrolier se pose dès lors que l'on procède à des études de simulation de réservoir. D'où l'intérêt de la méthode selon l'invention qui propose un d'outil de calcul rapide, de précision compatible avec la précision des données d'entrée.
Les notations suivantes sont utilisées au cours de la description de l'invention :

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.

Claims

Méthode pour optimiser la récupération d'un fluide en place dans un milieu poreux hétérogène, par injection dans ledit milieu d'un fluide de balayage pour provoquer un écoulement du fluide en place, lesdits fluides étant immiscibles, dans laquelle on discrétise ledit milieu en une grille constituée d'un ensemble de mailles, la méthode étant caractérisée en ce qu'elle comporte les étapes suivantes:
- 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, et dans laquelle ladite relation comporte un premier terme représentant un couplage visqueux décrivant la perturbation due à la variation de saturation, et un second terme décrivant les perturbations de la position du front provoquées par les hétérogénéités du milieu.

- 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.
Méthode selon la revendication 1, dans laquelle ledit premier terme est obtenu en considérant ledit milieu homogène et les déplacements dudit front de type Bucley-Leverett.
Méthode selon la revendication 1, dans laquelle ledit second terme est 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.
Méthode selon la revendication 2, dans laquelle ledit premier terme est 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.
Méthode selon la revendication 3, dans laquelle ledit second terme est 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.
Méthode selon l'une des revendications précédentes, dans laquelle on optimise l'injection en déterminant en chaque maille et pour différents pas de temps la saturation d'au moins un desdits fluides à partir de la position dudit front.