Brevet FR2870358A1 — publié le 20051118

La présente invention concerne une méthode rapide de simulation des écoulements dans un milieu poreux hétérogène tel qu'un réservoir pétrolier ou un aquifère, efficace pour tous types de contrastes de perméabilité. Cette méthode s'applique à des données pétrophysiques telles que les perméabilités, acquises entre autre par des mesures diagraphiques, sismiques et expérimentales (mesures sur carottes par exemple). La méthode peut être mise en oeuvre par exemple dans le domaine de la production pétrolière par des ingénieurs de gisement en vue d'obtenir rapidement des prédictions d'écoulement dans un milieu peu hétérogène à très hétérogène, voire discontinu comme par exemple un réservoir fissuré. L'une des propriétés les plus importantes du milieu poreux est sa perméabilité K, définie par la loi de Darcy (1856). Du fait de l'hétérogénéité du milieu poreux, la perméabilité K est souvent assimilée à une variable aléatoire. De fait, le milieu poreux est représenté par un champ de perméabilité qui varie dans l'espace. L'écoulement d'un unique fluide en régime permanent dans un milieu poreux est régi par la loi de Darcy: U= Kgrad(P)= KVP 77 77 où U est la vitesse de filtration, K la perméabilité, P la pression du fluide et 17 la viscosité du fluide. Par ailleurs, le fluide étant incompressible, l'équation de conservation de la masse s'écrit: div(U) = 0 En combinant les deux équations précédentes, on obtient l'équation de diffusivité : div(K VP) = 0 17 Quelle que soit la technique utilisée pour simuler des écoulements dans un milieu poreux hétérogène, il faut résoudre cette équation de diffusivité. En d'autres termes, il faut estimer, pour un champ de perméabilité donné et discrétisé sur un maillage de réservoir, les champs de pression et vitesses correspondants. Etat de la technique Les documents suivants que l'on va citer dans le cours de la description ci-après, illustrent l'état de la technique: Darcy, H., Les fontaines publiques de la ville de Dijon, Victor Dalmont, Paris, France, 1856. - Eyre, D.J., et Milton, G.W., A fast numerical scheme for computing the response of composites using grid refinement. The european physical journal applied physics, 6:41-47, 1999. - Michel J.C., Moulinec H., et Suquet P., A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast, Comput. Model. Eng. Sc., 1(2), 79-88, 2000. - Mizell S.A., Gutjahr A.L., et Gelhar L.W., Stochastic analysis of spatial variability in two-dimensional steady groundwater flow assuming stationary and non-stationary heads, Water Resources Research, 18(4), 1053-1067, 1982. - Moulinec, H., et Suquet, P., A numerical method for computing the overail response of nonlinear composites with complex microstructure. Computer methods in applied mechanics and engineering, 157:69-94, 1998. - Van Lent, T. J., Numerical spectral methods applied to flow in highly heterogeneous aquifers, PhD thesis, Stanford University, 1992. De nombreux algorithmes ont été développés dans le passé pour simuler les écoulements dans des milieux poreux hétérogènes. Ils impliquent essentiellement des différences finies, des volumes finis, des éléments finis mixtes hybrides, des éléments finis conformes ou encore des éléments finis mixtes spectraux. Ces techniques sont très générales et permettent de modéliser des écoulements très complexes. Elles se caractérisent en ce qu'elles s'appliquent aux transmissivités plutôt qu'aux perméabilités. A ce jour, il n'existe pas de méthode simple et totalement fiable pour estimer ces transmissivités à partir des perméabilités. Parallèlement à ces différentes approches, des méthodes permettant de travailler directement à partir des perméabilités ont été développées. Ces méthodes privilégiant l'emploi des transformées de Fourier pour simuler les écoulements ont été introduites par Van Lent (1992). Van Lent développe l'équation de diffusivité suivant la méthode des petites perturbations et propose de la résoudre dans le domaine fréquentiel par une approche itérative. Cependant, cette méthode n'est appropriée que pour les milieux poreux continus et peu hétérogènes. Elle n'est pas adaptée à la simulation des écoulements dans un milieu fortement hétérogène du fait des approximations associées à la résolution de l'équation de diffusivité développée suivant la méthode des petites perturbations. Cette méthode ne peut pas, par exemple, s'appliquer à un milieu contenant deux faciès comme les grès et l'argile dont les propriétés pétrophysiques sont très différentes. Ces techniques, qui travaillent directement à partir des perméabilités, c'est à dire sans passer par les transmissivités, sont plus rapides que les méthodes complexes précitées, mais elles sont limitées à la description d'un écoulement monophasique en régime permanent dans un milieu poreux peu hétérogène. De fait, leur utilisation est restée très marginale. La méthode selon l'invention permet de simuler rapidement l'écoulement d'un fluide unique et incompressible dans un milieu poreux quelle que soit son hétérogénéité. Elle s'appuie sur le développement de l'équation de diffusivité dans le domaine spectral selon la formule de Neumann. La méthode permet alors de résoudre cette équation par un modèle itératif, grâce à des algorithmes faisant intervenir des transformées de Fourier rapides. Les schémas de résolution proposés permettent un calcul rapide des champs de pression et vitesse et sont une alternative aux simulateurs développés à ce jour dans l'industrie. La méthode selon l'invention L'invention concerne une méthode pour simuler, de façon simple et rapide, l'écoulement d'un fluide unique et incompressible dans un milieu poreux hétérogène, tel qu'une zone souterraine, quel que soit son degré d'hétérogénéité. La méthode comporte les étapes suivantes: - on discrétise ledit milieu hétérogène par un maillage; - on détermine un champ de perméabilité à partir de mesures obtenues par au moins l'une des ces techniques: mesures diagraphiques, mesures sismiques et mesures expérimentales (mesures sur carottes par exemple); - on traduit dans le domaine fréquentiel l'équation de diffusivité, reliant la perméabilité (K) au gradient de pression (VP) et à la vitesse de filtration (U), puis on développe cette équation selon la formule de Neumann; - on détermine en chaque maille dudit maillage, le gradient de pression (VP) et la vitesse de filtration (U), par résolution itérative de l'équation de diffusivité ainsi traduite, comportant l'utilisation de transformées de Fourier rapides pour accélérer les calculs. Selon l'invention, la détermination en chaque maille dudit gradient de pression et de ladite vitesse de filtration peut être réalisée directement à partir de la valeur dudit champ de perméabilités. Cette détermination en chaque maille du gradient de pression et de la vitesse de filtration peut comporter les étapes suivantes: - on définit, selon le degré d'hétérogénéité dudit champ de perméabilité, une variable d'itération, en fonction d'au moins un des paramètres pétrophysiques suivants: la vitesse de filtration (U) et le gradient de la pression (VP) ; - on détermine la vitesse de filtration (U) et le gradient de la pression (VP) par un calcul itératif basé sur ladite variable d'itération. Selon l'invention, dans le cas d'un champ de perméabilité peu hétérogène, la variable itérative peut être la vitesse de filtration (U). Dans le cas d'un champ de perméabilité très hétérogène, la variable itérative peut être une vitesse (z- ), fonction d'au moins un des paramètres suivants: la perméabilité moyenne (K ), la perméabilité (K) et le gradient de pression (VP). Enfin, dans le cas de contrastes de perméabilité infinis, par exemple dans le cas d'une matrice peu perméable traversée par des fractures très perméables, la variable itérative peut être le multiplicateur de Lagrange X. Présentation succincte des figures Les caractéristiques et avantages de la méthode selon l'invention, apparaîtront plus clairement à la lecture de la description ci-après d'un exemple non limitatif de réalisation, en se référant aux dessins annexés où : - la figure 1 montre un exemple de champ de perméabilité continue; - les figures 2A, 2B et 2C montrent les pressions calculées à partir d'un solveur par différences finies (DF), du premier algorithme proposé (FF1'1) et du deuxième algorithme proposé (FF1'2) ; - les figures 3A, 3B et 3C montrent les vitesses de Darcy suivant l'axe X calculées à partir d'un solveur par différences finies (DF), du premier algorithme proposé (FF1'1) et du deuxième algorithme proposé (FF12) ; - les figures 4A, 4B et 4C montrent les vitesses de Darcy suivant l'axe Y calculées à partir d'un solveur par différences finies (DF), du premier algorithme 15 proposé (FFT1) et du deuxième algorithme proposé (FFT2) ; - la figure 5 montre un exemple de champ de perméabilité en faciès (le faciès noir a une perméabilité de 20 mD et le faciès blanc une perméabilité de 500 mD) ; - les figures 6A, 6B et 6C montrent le champ de pression correspondant à la figure 5 calculé par volumes finis (VF), le champ de pression correspondant calculé par transformées de Fourier rapide (FFT1), et le champ de pression correspondant calculé par transformées de Fourier à base de Lagrangien augmenté (1''F1'3). Description détaillée de la méthode La méthode de simulation selon l'invention permet de prédire de façon rapide l'écoulement des fluides d'un milieu hétérogène. Les différentes étapes de la méthode sont présentées dans le cas où le milieu est un gisement pétrolier, mais la méthode peut s'appliquer à tout type de milieu poreux hétérogène. L'objectif fixé est de fournir à l'ingénieur de gisement un outil lui permettant de rapidement améliorer la description des modèles de réservoir et ainsi la production d'hydrocarbures. Pour ce faire, on se propose de développer un simulateur d'écoulement rapide, limité à une physique simple: on se place dans le cadre d'un fluide unique et incompressible, et dans le cas d'écoulement en régime permanent. Ce simulateur peut s'appliquer à un milieu perméable hétérogène quel que soit son degré d'hétérogénéité. 1- Définition du problème La modélisation des écoulements dans un réservoir pétrolier est fondée essentiellement sur l'application, au réservoir préalablement maillé, de la loi de diffusivité, combinaison de la loi de Darcy et de l'équation de conservation de la masse pour un fluide incompressible. Afin de simplifier les développements qui suivent, on suppose que la viscosité du fluide vaut 1 et que la perméabilité correspond à un scalaire en chaque maille. Ces simplifications n'enlèvent rien à la généralisation des algorithmes qui sont proposés dans la suite. Les équations s'écrivent alors de la façon suivante: U = KVP div(U) = 0 Le gradient de pressions VP, la vitesse de filtration U et la perméabilité K, sont considérées comme des variables aléatoires. La méthode s'affranchit de toute hypothèse réductrice concernant le degré d'hétérogénéité du milieu poreux, en résolvant cette équation de diffusivité dans le domaine spectral développée selon la formule de Neumann. La résolution de cette équation permet alors de définir, directement à partir des perméabilités définies en chaque maille du maillage de réservoir, les valeurs de gradient de pression et de vitesse de filtration. Cette étape fait appel à un algorithme itératif comportant, pour accélérer les calculs, l'utilisation de transformées de Fourier rapides telles que décrites dans le document suivant: Cooley, J.W., et Tukey, J.W., 1965, "An algoritm for the machine computation of complex Fourier series", Math. Comp., 19, 297-301. En effet, la complexité de l'algorithme des transformées de Fourier rapides est de nlog(n) au lieu de n2 pour un algorithme classique de transformées de Fourier, ce qui implique un gain de temps de calcul considérable. La méthode peut alors se décomposer en 4 étapes: - on discrétise le gisement par un maillage; - on détermine un champ de perméabilité à partir de mesures effectuées sur des puits (diagraphies), de mesures sismiques et/ou de mesures expérimentales en laboratoires (mesures sur carottes par exemple). - on traduit l'équation de diffusivité, reliant la perméabilité au gradient de pression et à la vitesse de filtration, dans le domaine fréquentiel, puis on développe cette équation selon la formule de Neumann. - on détermine en chaque maille du maillage de réservoir, le gradient de pression ainsi que la vitesse de filtration par résolution itérative de l'équation de diffusivité ainsi traduite. Trois algorithmes sont proposés, chacun d'eux comportant l'utilisation de transformées de Fourier rapides pour accélérer considérablement les calculs. Les deux premières étapes peuvent être réalisées par des méthodes bien connues de l'homme du métier. On peut par exemple, pour discrétiser le réservoir, utiliser un maillage de type cartésien uniforme. Ce type de maillage est tout à fait classique pour un ingénieur de réservoir. Le champ de perméabilité peut, par exemple, être issu d'estimations géostatistiques effectuées à partir des mesures réalisées aux puits par diagraphies, des mesures de laboratoires permettant de déterminer la perméabilité de carottes issues du réservoir étudié, et de données sismiques pour, par exemple, définir la variabilité spatiale de la perméabilité. De nombreuses techniques, reposant sur les estimations par krigeages et bien connues de l'homme du métier, peuvent être utilisées. La façon d'obtenir un maillage du réservoir renseigné par un champ de perméabilité est donc aisément réalisable par un ingénieur de réservoir. En revanche, l'étape de développement de la formule de diffusivité et l'étape de simulation des écoulements des fluides sont présentées ciaprès. 8 2- Développement de l'équation de diffusivité L'idée est tout d'abord de développer une formulation itérative de l'équation de diffusivité, permettant de tenir compte de tous types d'hétérogénéité, avant de proposer des algorithmes de résolution. Problème auxiliaire On se concentre, dans un premier temps, sur un problème auxiliaire qui est celui de l'écoulement d'un fluide dans un milieu poreux homogène de perméabilité K soumis à un champ de vitesse r. Les équations de base deviennent alors: U = K VP+r div(U) = 0 Ce qui donne dans le domaine fréquentiel: ê = iK P + iê = 0. Le symbole "^" indique que l'on considère la transformée de Fourier, est le vecteur d'onde (discrétisation de l'espace des fréquences) et i est tel que: i2 = -1. En substituant ê dans l'équation de conservation de la masse, on montre que: OP= Cette équation se ramène à : P = Î'0z avec r est appelé "opérateur de Green" ou "solution fondamentale associée à K ". Problème général On utilise à présent le problème auxiliaire pour résoudre le problème de l'écoulement d'un fluide dans un milieu poreux hétérogène de perméabilité K et soumis 25 à un gradient de pression OP. Le gradient de pression est la somme de deux termes: une moyenne, appelée gradient de pression macroscopique (Vpmacro) et une perturbation (Vp). U = KVP = K(Vp +VP') div(U) = 0 Par souci de simplification, on suppose que le gradient VPmar est fixé. 5 Néanmoins, d'autres conditions pourraient être envisagées (Moulinec and Suquet, 1998). On introduit un champ de référence, de perméabilité homogène K : U = KOP + K VP K VP = (K K )VP K VP =r K VP Le champ Test alors défini par z = (K K)VP. Avec cette formulation, la loi de Darcy et l'équation de conservation de la masse deviennent: TU = K0(Vp+VPmacro)+r L div(U) = 0 La solution déterminée pour le problème auxiliaire nous amène à proposer la solution générale suivante: VP(ç) = r î V Ç 0 VP(0) = Vpmacro ce qui donne dans l'espace réel: VP=F0 *r+VPmacro. Formulation itérative En tenant compte de l'expression de a, on obtient pour VP: VP = (K K)F0 *VP+VPmacro 20 On montre finalement que: (I (K K)I'0)*VP =Vpmacro Si, à présent, on définit A tel que: A=(K K) F0, le développement suivant la formule de Neumann nous conduit à : (I A) ' =I+A+A2+A3+... Le gradient de pression devient alors: VP = E ((K * Vpmacro 1=0 Pour calculer cette expression, on introduit: k vp) =1r(Ko _K)r0) *vpmacro k=0 qui fournit finalement la formule itérative suivante: VP' =Vpmacro+(K0 K) I,0 *Vp3-'. En combinant cette expression avec la loi de Darcy, on montre que: VP3 =ppmacro +r *U' 1 +Koto * 7p1' 3- Simulation des écoulements La formulation itérative précédente constitue le c ur des trois algorithmes proposé ci-après. Les algorithmes itératifs proposés permettent de résoudre directement cette équation, c'est à dire de résoudre l'équation de diffusivité sans aucune approximation concernant le type d'hétérogénéité. Premier algorithme Le premier algorithme, qui gère directement les pressions et vitesses, est général. 1 - Initialisation: La valeur initiale du gradient de pression est donnée par le gradient de pression macroscopique: VP = Vpmacro La vitesse initiale vaut: U = KVPmacro. En utilisant un algorithme de transformées de Fourier rapides (FFT), on calcule VP . 2 - Itération j+l: DP' et U' sont connus. 1) en utilisant un algorithme de F14'1, on calcule ê ; 2) on pose vP'+' () + K r ().vP' () quand 0 et OP3+1(0) _ v pmacro 3) on utilise un algorithme de FFI' inverse pour calculer VP +1; 4) on calcule U'+ 5) on calcule la perméabilité effective K eff' = (vp i+1 1 / (U'+i) 6) on pose j = j+l et on retourne à l'étape 1 jusqu'à ce que K;+ eff eff E. Parce que l'on utilise des transformées de Fourier rapides, les calculs sont très rapides. Néanmoins, cela nous oblige à considérer des conditions aux limites périodiques et des grilles, modélisant le milieu poreux hétérogène, régulières. Le critère d'arrêt proposé s'appuie sur le calcul de la perméabilité effective, mais d'autres critères peuvent être envisagés. Le choix du K est crucial pour la convergence. Empiriquement, on montre que la valeur K = mm(K) 2max(K) est judicieuse. Le taux de convergence est donné (Eyre et Milton, 1999) par: max(K) OP' = y; I ppmacro avec: y = min(K) max(K) min(K) Comme yest inférieur à 1, l'algorithme développé converge pour toute valeur de départ. Cependant, plus le contraste de perméabilité est important et plus la convergence est lente. C'est pourquoi un deuxième algorithme est proposé : 2870358 12 Deuxième algorithme Le second algorithme, qui introduit une variable intermédiaire mise à jour à chaque itération, permet d'une part, de supporter des contrastes de perméabilité plus importants que le premier algorithme, et d'autre part, de traiter plus rapidement ces contrastes importants de perméabilité. Pour rendre le processus plus rapide, même dans le cas de champs de perméabilité affichant des contrastes importants de perméabilité, on s'inspire des travaux de Eyre et Milton (1999) développés dans le cadre des permittivités diélectriques. Suivant ces auteurs, la variable r(x) peut s'écrire sous une forme différente: z = 2(1 yz) 1 KOVPmacro avec: I l'identité y = 2K I'0 I K + K z K K A nouveau, en se référant à un développement de Neumann, on met en évidence la formulation itérative suivante pour z: ri' = 2K VPmacro +z;yz A partir de cette expression, on construit l'algorithme suivant: 1 - Initialisation: La valeur initiale du gradient de pression est donnée par le gradient de pression 20 macroscopique: VP = ppmacro. On a alors r0 = 2KoO pmacro 2 - Itération j+l: 1) on calcule Q'+' zr' 2) en utilisant un algorithme de FFT, on calcule Q'+' ; 3) on calcule (4) = z, V 0 et î'+' (0) = î (0) ; 4) en utilisant un algorithme de FFT inverse, on calcule ri' ; 5) on pose j = j+l et on retourne à l'étape 1 jusqu'à ce que Après convergence, on estime VP et U à partir de z et^ de la formule a=(K K) VP. A nouveau, le taux de convergence dépend de K . D'après Eyre et Milton (1999), la convergence est optimale quand K = Vmin(K)max(K) et le taux de convergence est donné par: \7P' vP'-'I = y; vpmacrol avec: y = Comme y est inférieur à 1, cet algorithme converge quel que soit le point de 10 départ. Il converge aussi plus rapidement que le premier algorithme comme le démontreront les applications présentées ultérieurement. Troisième algorithme Les méthodes précédentes permettent une résolution rapide de l'équation de la diffusivité. Nous avons pu montrer que le taux de convergence de ces algorithmes est lié au contraste de perméabilité. Ces deux algorithmes permettent de traiter des contrastes de perméabilité très importants, mais pas infinis. Ce type de milieu peut par exemple correspondre à une matrice peu perméable traversée par des fractures très perméables. En effet les réservoirs fissurés constituent un type extrême de réservoirs hétérogènes comportant deux milieux contrastés, un milieu matriciel contenant la plus grande part de l'huile en place et présentant une faible perméabilité, et un milieu fissuré représentant moins de 1% de l'huile en place et hautement conducteur. Dans ce dernier cas de contrastes infinis, on envisage un nouveau type d'algorithme à base de lagrangien augmenté, parallèlement aux travaux de Michel et al. (2000). Il s'agit à présent de voir le problème de simulation d'écoulement comme un problème de minimisation d'énergie. Ce problème est classiquement abordé en trois temps selon l'algorithme d'Uzawa. Tout d'abord, on résout l'écoulement de façon globale. Puis, on résout une équation de dissipation de l'énergie locale. Enfin, on met à jour le multiplicateur de Lagrange. L'algorithme de simulation d'écoulement se déroule comme suit. 1 - Initialisation: La valeur initiale du gradient de pression est donnée par le gradient de pression macroscopique: VP = VP'0 La valeur initiale du multiplicateur de Lagrange est supposée égale à l'opposée de la vitesse: = K Vp . On introduit la variable Ex, dont la valeur initiale est l'opposée du gradient de pression macroscopique: E _ OP . 2 Itération j+l: OPj, 2j, Ej sont connus. 1) On calcule zj' = .j + K Ex 2) en utilisant un algorithme de FFT, on calcule îJ+1; 3) on calcule OPj+1 ( ) = j+ll) si 0 "7p (o) ll sinon 4) en utilisant un algorithme de FFT inverse, on calcule VPj+1 5) on en déduit Ej+' = VP j+1 1 (2j + KOPi+1) Y K +K 6) on en déduit Xj+1 = 2j _ K (OP E j+l) ; 7) on pose j = j+l et on retourne à l'étape 1 jusqu'à ce que j+l 2j A partir de quelques tests, on montre que K = Vmin(K)max(K) est une 20 valeur optimale. 4- Applications Pour illustrer le potentiel des algorithmes développés, on effectue des tests comparatifs en deux dimensions. On considère dans un premier temps, un champ de perméabilité construit sur une grille bidimensionnelle de 50x50 mailles comme illustré par la Figure 1. La distribution de ce champ de perméabilité est de type log- normale. La moyenne et la variance du logarithme des perméabilités valent respectivement 1 et 3. La distribution des valeurs du logarithme de la perméabilité respecte un variogramme sphérique isotrope d'une portée de 20 mailles. Pour ce champ de perméabilité, on résout le problème d'écoulement décrit plus haut en utilisant d'une part un solveur classique par différences finies et d'autre part les algorithmes à base de transformées de Fourier rapides proposés ici. Les résultats en pression sont reportés sur les figures 2A, 2B et 2C, les résultats en vitesses selon l'axe horizontal X sont reportés sur les figures 3A, 3B et 3C, et enfin, les résultats en vitesses selon l'axe vertical Y sont reportés sur les figures 4A, 4B et 4C. Les figures 2A, 3A et 4A représentent les résultats obtenus avec un solveur de type différences finies (DF). Les figures 2B, 3B et 4B représentent les résultats obtenus avec le premier algorithme (FFT1) et les figures 2C, 3C et 4C représentent les résultats obtenus avec le second algorithme (FFT2). Clairement, les résultats sont très semblables. Il existe cependant de légères différences qui peuvent s'expliquer, au moins en partie, par la nature même des solveurs. Par exemple, le solveur par différences finies fournit les vitesses aux interfaces entre les mailles alors que les algorithmes à base de FFT proposé ici donnent les vitesses au centre des mailles. Il y a donc un décalage en vitesse. En outre, les solveurs par différences finies ne travaillent pas directement avec les perméabilités: ils font intervenir des transmissivités. La transmissivité entre deux mailles dépend des perméabilités de ces deux mailles. A ce jour, il n'existe pas à de technique universelle pour calculer ces transmissivités. En général, on les estime à partir de la moyenne harmonique. Cette approche est justifiée à une dimension, mais très improbable pour des dimensions supérieures. Au contraire, les algorithmes proposés utilisent directement les perméabilités: ils ne nécessitent pas l'utilisation de transmissivités pour lesquelles on ne dispose que d'une valeur approchée. Une analyse plus fine est reportée dans le Tableau 1 pour comparer les deux premiers algorithmes par transformées de Fourier rapides. Les tests sont effectués pour le même champ de perméabilité, celui illustré sur la figure 1, mais des valeurs de e décroissantes. On rappelle que e contrôle le processus d'arrêt des algorithmes. On observe que le deuxième algorithme est nettement plus efficace que le premier. Pour un même e, le deuxième algorithme requiert moins d'itérations. De plus, la perméabilité effective se stabilise plus rapidement avec le deuxième algorithme qu'avec le premier. Pour le cas étudié, la perméabilité effective vaut 15. 71. Le premier algorithme nécessite 290 itérations pour atteindre cette valeur, alors que le deuxième n'en demande que 21. Tableau 1. Comparaison entre les deux premiers algorithmes proposés. E est une constante qui détermine l'arrêt du processus itératif. Keff est la perméabilité effective à la dernière itération. e 1e1 algorithme 2eme algorithme Nombre d'itérations Kef Nombre d'itérations Keff 10.2 27 15.83 21 15.71 10-3 58 15.74 31 15.71 10-4 102 15.72 41 15.71 10-5 158 15.72 49 15.71 10-6 290 15.71 59 15.71 Un deuxième exemple est donné par un réservoir en faciès illustré sur la figure 5. Le milieu est constitué de 2 faciès. Le faciès dominant, présent à 70%, a une perméabilité de 500 mD: il s'agit du faciès réservoir (en blanc sur la figure 5). Le faciès non-réservoir (en noir sur la figure 5), présent à 30%, a une perméabilité de 20 mD. Dans ce cas, seul l'algorithme par FFT à base de Lagrangien augmenté (FFT3 sur la figure 6C) converge. Il conduit à une solution (figure 6C) assez semblable à celle déterminée par volumes finis (VF sur la figure 6A) , même si, comme noté plus haut, des différences persistent. Ces différences peuvent ici encore être liées à la façon dont les perméabilités sont traitées. Dans le cas des volumes finis, on utilise des transmissivités, calculées comme la moyenne harmonique des perméabilités de mailles adjacentes. Ce faisant, on a tendance à lisser les discontinuités en perméabilité. Les algorithmes par FFT (11' 1'1 sur la figure 6B) considèrent directement les perméabilités et sont donc vraisemblablement plus sensibles aux discontinuités de perméabilités. On observe en outre que dans le cas du premier algorithme de résolution parFFT proposé ici, même s'il n'y a pas convergence, la première itération suffit pour approcher assez bien le champ de pression. REVENDICATIONS 1) Méthode pour simuler l'écoulement d'un fluide unique et incompressible dans un milieu poreux hétérogène, tel qu'une zone souterraine, comportant les étapes suivantes: - on discrétise ledit milieu hétérogène par un maillage; - on détermine un champ de perméabilité à partir de mesures obtenues par au moins l'une des ces techniques: mesures diagraphiques, mesures sismiques et mesures expérimentales; - on traduit dans le domaine fréquentiel l'équation de diffusivité, reliant la perméabilité (K) au gradient de pression (VP) et à la vitesse de filtration (U), puis on développe cette équation selon la formule de Neumann; - on détermine en chaque maille dudit maillage, le gradient de pression (VP) et la vitesse de filtration (U), par résolution itérative de l'équation de diffusivité ainsi traduite, comportant l'utilisation de transformées de Fourier rapides pour accélérer les calculs. 2) Méthode selon la revendication 1, dans laquelle la détermination en chaque maille dudit gradient de pression et de ladite vitesse de filtration est réalisée directement à partir de la valeur dudit champ de perméabilités. 3) Méthode selon l'une des revendications 1 et 2, dans laquelle la détermination en chaque maille dudit gradient de pression et de ladite vitesse de filtration comporte les étapes suivantes: - on définit, selon le degré d'hétérogénéité dudit champ de perméabilité, une variable d'itération, en fonction d'au moins un des paramètres pétrophysiques suivants: la vitesse de filtration (U) et le gradient de la pression (VP) ; - on détermine la vitesse de filtration (U) et le gradient de la pression (\7P) par un calcul itératif basé sur ladite variable d'itération. 4) Méthode selon la revendication 3, dans laquelle la variable itérative est la vitesse de filtration (U). 5) Méthode selon la revendication 3, dans laquelle la variable itérative est une vitesse (a), fonction d'au moins un des paramètres suivants: la perméabilité moyenne (K ), la perméabilité (K) et le gradient de pression (OP). 6) Méthode selon la revendication 3, dans laquelle la variable itérative est le multiplicateur de Lagrange X.