FR 2850187 A1 20040723 FR 0300430 A 20030116 La présente invention concerne une méthode pour déterminer des zones d'un milieu poreux 10 stratifié dans lesquelles le front ou interface entre des fluides en place et des fluides de balayage se déplace de façon stationnaire, c'est-à-dire sans déformation et à vitesse constante. La méthode trouve des applications notamment quand on a à balayer par des fluides injectés (de l'eau par exemple) un gisement d'hydrocarbures stratifié. Cette détermination permet d'obtenir très rapidement les zones pouvant être agrégées et constituer des unités 15 hydrodynamiques lors de la mise à l'échelle d'un modèle de réservoir. Etat de la technique La simulation des écoulements diphasiques en milieu poreux hétérogène peut 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 20 domaine pétrolier. Afin de pouvoir procéder à la simulation d'écoulements avec des moyens raisonnables, il faut utiliser une description réduite du réservoir. Une technique consiste à agréger les mailles du modèle numérique de façon à obtenir un maillage plus grossier, composé de mailles aux valeurs effectives moyennes. Cette technique, dite des pseudosfonctions, a initialement été proposée pour des réservoirs stratifiés de façon à agréger les 25 mailles dans la direction verticale par - Coats, K.H., Nielsen, R.L. and Terhune, M.H. "Simulation of threedimensional, twophase flow in oil and gas reservoirs." SPE 1961, 1967 Cette méthode, reprise et étendue par de nombreux auteurs, repose sur des hypothèses d'équilibre vertical capillaire ou d'équilibre gravitaire o les forces capillaires ou gravitaires 30 sont prédominantes devant les effets dus à la viscosité, et la distribution des saturations dans le réservoir peut être connue sans recourir à une simulation sur le maillage précis (simulation fine). On peut en tirer facilement des valeurs moyennes suivant la direction verticale. Une seule couche aux valeurs moyennes peut finalement suffire à décrire l'évolution des saturations dans le réservoir. Une autre méthode conduisant à obtenir des valeurs moyennes est proposée par: - Hearn, C.L. "Simulation of stratified waterflooding by pseudo-relative permeability curve." Journal of Petroleum Technology, pages 805-813, July 1971. Elle repose sur l'hypothèse d'équilibre vertical d aux effets visqueux. Dans ce cas, les rééquilibrages visqueux dans le réservoirs sont très rapides par rapport aux variations de saturation, de telle sorte que des valeurs moyennes peuvent toujours être calculées sans 10 recourir à la simulation fine, car les fluides se propagent à vitesse constante dans chaque couche. La méthode de Hearn, ultérieurement étendue par, entre autres par: - Simon, A.D. and Koederitz, L.F. "An improved method for the determination of pseudorelative permeability data for stratified systems. " SPE 10975, 1982, se base sur un calcul itératif conduisant à ordonner les couches du réservoir en fonction de la 15 vitesse de propagation des fluides dans celles-ci. On peut montrer que ce calcul itératif peut ne pas converger. On est alors en présence d'un front stationnaire dans au moins deux couches du réservoir. Les auteurs proposent alors de ne prendre qu'une seule couche aux propriétés moyennes pour ces zones particulières. Une étude théorique de l'équilibre vertical et de ses implications peut être trouvée chez - Yortsos, Y.C. "Analytical studies for processes at vertical equilibrium. " SPE 26022, 1992. Lorsque la capillarité ou la gravité ne peuvent être considérées comme prédominantes sur les effets visqueux, mais que l'équilibre vertical visqueux ne peut se mettre en place, les méthodes proposées jusqu'à présent ne peuvent se passer d'une simulation fine sur tout ou partie du réservoir pour calculer les propriétés moyennes pouvant être affectées aux blocs 25 grossiers. On parle alors de pseudo fonctions dynamiques. Celles-ci furent introduites notamment par: - Kyte, J.R. and Berry, D.W. "New pseudo functions to control numerical dispersion." SPE 5105, 1975. Ces méthodes dynamiques posent toutefois de nombreux problèmes théoriques et pratiques. Des revues de ces méthodes et des problèmes associés peuvent être trouvées dans les publications suivantes - Archer, R. "Pseudo function generation", Master of Science thesis, Departement of Petroleum Engineering of Stanford University, 1996. - Ahmadi, A. "Utilisation des propriétés équivalentes dans les modèles de réservoir: cas des écoulements diphasiques incompressibles." Thèse de Doctorat, Université de Bordeaux I, 1992; - Barker J.W. and Thibeau S. "A critical review ofthe use ofpseudo relative permeabilities 10 for upscaling ". SPE 35491, 1996. Les méthodes dynamiques ont en général recours à une agrégation systématique et uniforme des mailles. Afin de prendre en compte les influences locales des hétérogénéités sur l'écoulement, des méthodes d'agrégation non-uniformes ont été proposées notamment par: - Darman, N.H. and Durlofsky, L.J. "Upscaling immiscible gas displacements: Quantitative 15 Use of Fine Grid Flow Data in Grid Coarsening Schemes." SPE 59452, 2000. Ces méthodes ont l'avantage de n'agréger que des zones préférentielles et laissent en état initial les parties du réservoir à caractères hydrodynamiques plus complexes. Elles ne peuvent cependant s'affranchir d'une simulation fine. Des méthodes permettant d'agréger de façon non uniformes les strates d'un réservoir, sans 20 avoir recours à une simulation détaillée, ont été proposées. On peut citer des méthodes entièrement statiques, comme: - Li, D. and Beckner, B. "Optimal uplayering for scaleup of multimillion- cell geologic models", SPE 62927, 2000 ou basées sur les résultats d'une simulation d'écoulement monophasique (qui est beaucoup 25 rapide qu'une simulation d'écoulement diphasique), comme: - Stern, D. and Dawson, A.G. "A technique for generating reservoir simulation grids to preserve geologic heterogeneity", SPE 51942, 1999. De telles méthodes, par définition, ne prennent donc pas en compte les effets dynamiques et visqueux des écoulements diphasiques. L'évolution du front dans le réservoir au cours de l'écoulement est considérablement influencée par le couplage visqueux entre le champ de pression et le champ de saturation. 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 5 l'écoulement des fluides dans les couches les plus perméables. 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. Si, dans au moins une partie du réservoir, il existe un tel front stationnaire, tout se passe comme si l'écoulement n'y voyait qu'une seule couche au sens hydrodynamique. On peut donc raisonnablement penser qu'une seule couche aux propriétés moyennes dans cette zone suffit à la modélisation lors de la simulation numérique de l'écoulement. Ceci est d'importance capitale lorsque l'on passe du modèle géologique de réservoir avec toutes les couches, au 15 modèle de simulation, pour lequel seules les zones les plus influentes sur l'écoulement suffisent à la description. Deux problèmes se posent alors. Connaissant la viscosité du fluide injecté, il faut pouvoir déterminer a priori l'existence et la localisation de zones o l'écoulement est stationnaire, sans recourir pour autant à une simulation complète d'écoulement sur le modèle détaillé. Il faut 20 également être en mesure avec la même contrainte, de choisir la viscosité du fluide injecté de telle sorte qu'un maximum de couches aient un comportement stationnaire. La méthode que nous proposons permet de reconnaître les zones d'un réservoir stratifié se comportant comme une seule couche au sens hydrodynamique sans recourir à une simulation fine et sans invoquer les hypothèses d'équilibre vertical. Elle facilite le choix des zones à 25 agréger dans les réservoirs stratifiés. Elle permet donc de prendre en compte les effets dynamiques et visqueux tout en permettant une détermination des couches grossières très rapide par rapport aux anciennes solutions. La méthode selon l'invention La méthode selon l'invention permet de déterminer dans un milieu stratifié dont les propriétés 30 physiques sont connues ou estimées, au moins une zone o une interface entre un fluide en place dans le milieu et un fluide de balayage, de viscosités et densités différentes connues, injectés dans le milieu, se déplace de façon stationnaire, dans le but de simplifier la construction d'un modèle de simulation des écoulements dans le milieu. Elle comporte essentiellement les étapes suivantes - on considère une forme d'interface a priori en supposant que les déplacements de fluides en tous ses points sont stationnaires; - on détermine le champ de pression de part et d'autre de l'interface a priori; - on modifie itérativement la forme de l'interface jusqu'à obtenir que les pressions de part et d'autre d'au moins une partie de l'interface s'égalisent en tous points de cette partie; et - on attribue des propriétés hydrodynamiques moyennes uniformément à chaque zone du 10 milieu délimitée par chaque partie d'interface, quand la dite égalisation est atteinte. Suivant un mode de mise en oeuvre, quand on ne peut obtenir une égalisation des pressions de part et d'autre de l'interface tout le long de celle-ci, on segmente l'interface en plusieurs parties et on modifie itérativement et séparément la forme de ces différentes parties, jusqu'à obtenir une égalisation des pressions de part et d'autre de chacune d'elles, l'étendue de chaque partie 15 d'interface, quand la dite égalisation est atteinte, délimitant une zone favorable à laquelle on attribue uniformément des propriétés hydrodynamiques moyennes. Suivant un exemple de mise en oeuvre, on détermine la forme d'au moins une zone du milieu délimitée par une interface à déplacement stationnaire, qui correspond à des valeurs différentes de la viscosité des fluides de balayage, et on choisit la viscosité pour laquelle on 20 optimise les déplacements stationnaires dans le dit milieu. 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 mise en oeuvre, en se référant aux dessins annexés o: - la figure 1 montre une boucle de détermination itérative de l'existence et de la forme d'un front stationnaire pour un écoulement diphasique de rapport de viscosités connu dans un milieu stratifié donné - la figure 2 montre une boucle de construction itérative de l'évolution des couches hydrodynamiques en fonction du rapport des viscosités, pour des écoulements diphasiques dans un milieu stratifié donné; - la figure 3 montre un exemple de réservoir stratifié incliné par rapport à la verticale. Cette figure a pour but d'expliciter les notations employées dans la description de la méthode; la figure 4 montre la valeur des perméabilités absolues en fonction de la profondeur, employées dans l'exemple synthétique de réservoir illustrant la méthode de façon nonlimitative; - la figure 5 montre l'évolution de l'interface au cours du processus itératif de recherche de 10 la forme d'un front stationnaire, sur l'exemple décrit par la figure 4; - la figure 6 montre l'évolution de la solution stationnaire en fonction du rapport des viscosités employées, sur l'exemple décrit par la figure 4; et - la figure 7 montre l'évolution du nombre de couches hydrodynamiques en fonction du rapport des viscosités employées, sur l'exemple décrit par la figure 4. Description détaillée de la méthode On considère un balayage de fluides de viscosités connues dans un milieu stratifié. Supposons qu'un front stationnaire de forme connue se développe dans le milieu considéré. Ceci impose une condition de flux au niveau de l'interface. Le flux est également connu sur les bords du réservoir (o il est en général nul), en entrée et en sortie (o il est imposé par les conditions 20 d'injection et de récupération). Comme la pression est une variable continue, le saut de pression de part et d'autre de l'interface est nul en tout point de celle-ci. On cherche à présent l'existence et la forme d'un front stationnaire dans le milieu. On procède de façon itérative en partant d'une forme de front a priori, en supposant que celle-ci correspond à la solution stationnaire. Comme le front est supposé stationnaire, on connaît le 25 flux à l'interface et l'on a donc un problème de Neumann parfaitement défini à résoudre pour chacune des deux zones que sépare l'interface. Le champ de pression peut donc être résolu indépendamment dans les deux régions. On évalue alors le saut de pression de part et d'autre de l'interface, en tout point de celle-ci. Si le saut de pression est nul en tout point, on possède alors la forme correcte du front correspondant à un écoulement stationnaire. Si le saut n'est pas nul partout, la condition de continuité des pressions n'est pas respectée et la forme de front ne correspond pas à une solution stationnaire. On doit alors déformer l'interface et recommencer le processus jusqu'à la convergence. Deux cas apparaissent alors: - Si la convergence est atteinte, on possède alors la forme du front correspondant à une 5 solution stationnaire de l'écoulement. Toutes les couches sont traversées par le fluide à la même vitesse, et l'ensemble du système se comporte comme une seule unité hydrodynamique. On peut donc ne considérer qu'une seule couche aux propriétés moyennes lors de la mise à l'échelle. - Si le processus itératif ne converge pas, c'est qu'un front stationnaire ne peut pas exister 10 sur l'ensemble du milieu. Lors de la déformation du front, il apparaît alors en général deux zones s'éloignant perpétuellement l'une de l'autre. On découpe alors le milieu en deux parties correspondants aux zones précédemment identifiées. Puis le processus de recherche d'une solution stationnaire reprend dans chacune de ces deux parties. En procédant ainsi de manière itérative, on finira par obtenir toutes les zones dans lesquelles 15 des fronts stationnaires peuvent se développer, et qui peuvent être considérées comme des unités hydrodynamiques. La méthode peut aussi être appliquée de façon itérative en faisant varier le rapport des viscosités des fluides considérés. On part d'un fluide injecté beaucoup plus visqueux que celui en place. On doit alors trouver rapidement une solution stationnaire sur toute l'épaisseur du 20 milieu. Lorsque l'on diminue le rapport des viscosités entre le fluide injecté et le fluide en place, la forme du front stationnaire, si celui-ci existe, doit être modifiée. On reprend alors la méthode décrite ci-dessus pour trouver la nouvelle forme du front, en utilisant comme forme de front initiale du processus itératif, la forme finale obtenue pour le précédent rapport des mobilités. Ceci permet de converger beaucoup plus rapidement vers la solution stationnaire. On peut ainsi procéder de façon itérative sur le rapport des viscosités. Au fur et à mesure que l'on décroît ce rapport, le couplage visqueux ne sera plus suffisant pour compenser les contrastes de perméabilité et une solution stationnaire ne pourra plus être trouvée sur toute l'épaisseur du milieu. On découpe alors le milieux en deux portions suivant la méthode décrite ci-dessus en cas de non-convergence, et la recherche se poursuit dans les deux portions de 30 milieu. En procédant ainsi de façon itérative sur le rapport des viscosités, on obtient le nombre de couches hydrodynamiques du milieu en fonction des viscosités considérées. Ce résultat pourrait permettre par exemple de connaître la viscosité du fluide d'injection qui permet d'optimiser la récupération. On considère une section verticale d'un réservoir stratifié selon une direction x inclinée par rapport à l'horizontale (figure 3). Les valeurs des perméabilités absolues sont connues dans 5 toutes les strates. L'épaisseur du réservoir est selon la direction z. A l'extrémité amont (x = 0) on injecte, à débit constant Qo imposé, un fluide incompressible de viscosité pl et de densité pi, qui pousse de façon immiscible le fluide initialement présent dans le réservoir, incompressible et de viscosité M2 et de densité P2. On suppose qu'il n'y a pas de pression capillaire entre les deux fluides. La récupération se fait à l'extrémité aval (x = L). Le flux est 10 nul à travers le toit (z = 0) et le mur du réservoir (z,= A). On suppose également que les courbes de perméabilités relatives sont uniformes dans le réservoir et que chacune de ces deux zones est à saturation constante en fluides considérés. La zone non encore balayée est donc à saturation irréductible Sir en fluide injecté et la zone balayée à saturation maximum Sm Ou saturation de front Sf selon les hypothèses de l'utilisateur. Le front de saturation est 15 arbitrairement placé au centre du réservoir, avec une forme présupposée qui correspond à notre a priori sur la forme du front stationnaire. Cette forme sera par la suite modifiée jusqu'à trouver la "vraie" solution stationnaire, si celle-ci existe. On doit donc résoudre un problème aux frontières libres. Les conditions aux limites et la condition de front stationnaire permettent de définir un 20 problème de Neumann dans les deux zones que sépare le front: - A l'extrémité amont: la condition d'injection impose un flux constant sur l'épaisseur du réservoir en entrée. Si le milieu est supposé suffisamment long et le profil de saturations déjà bien développé, on peut supposer que le gradient de pression en entrée du système est identique dans toutes les couches: al = Y-, +pgX ax =O K1i p1g. ap = Plg, Z az x=O o p est la pression au point considéré et K1 est la moyenne des perméabilités effectives en entrée sur toute l'épaisseur du milieu, et Q0 est le débit d'injection. Les perméabilités effectives prennent en compte l'effet des perméabilités relatives, en fonction des saturations choisies par l'utilisateur dans les zones amont et aval. - A l'extrémité aval: la condition d'incompressibilité impose que le flux sortant du système soit le même que le flux entrant. On obtient donc une condition similaire à la tranche de '5 récupération: ap = 2Qo + P2g'x ax. =L K K2 ap -x = p2g-z azx=L o K2est la moyenne des perméabilités effectives en sortie sur toute l'épaisseur du milieu. - de part de d'autre de l'interface: Comme le front est supposé stationnaire, il se déplace 10 horizontalement sans déformation à la vitesse de filtration U0, définie telle que: U0 = Q0 x A O x est le vecteur unitaire dans la direction horizontale. Cette condition impose le flux et donc le gradient de pression de part et d'autre de l'interface: K (plg-Vpl). n _ K (Z) (P2g -VP2) n K1P(z) =K2(z) 2 ]M1]'2 o n est le vecteur unitaire orthogonal à l'interface au point considéré, Vp, est le gradient de pression à l'interface du côté amont et Vp2 est le gradient de pression à l'interface du côté aval. - Equation en pression: L'équation en pression est donnée classiquement de part et d'autre de l'interface. Pour la zone aval: V IK(z)(pig_-Vp)j= M Pour la zone amont: K2(z) (P2g-VP)$=O U2 Les équations en pression et les conditions aux limites permettent de définir pour chaque zone de part et d'autre du front un problème de Neumann. Le champ de pression dans chaque zone 5 peut donc être calculé numériquement à une constante près pour chaque zone. On élimine ce problème en posant la condition supplémentaire suivante, permettant d'éliminer les constantes en recalant les pressions le long du front: X p1 (z)dz = P2 (z)dz = C o pl et p2 sont les pressions évaluée au front respectivement pour les zones amont et aval, et 10 C est définie arbitrairement par l'utilisateur et peut être nulle. Une fois la pression p calculée dans chaque région et recalée le long de l'interface, le saut de pression ôpa (z) à travers le front en tout point de ce dernier peut être calculé: a5pf (Z) = P2 (Z) -p1 (z) - Ap(g.x + g.z) Si le saut de pression est partout nul le long du front, ou du moins en deçà de la valeur 15 critique définie par l'utilisateur, on dispose de la forme du front correspondant à un déplacement stationnaire. Sinon, il faut déformer l'interface de façon à faire diminuer les sauts de pression jusqu'à obtenir la forme stationnaire. Ceci revient alors à résoudre un problème d'optimisation, pour lequel de nombreuses méthodes sont applicables. Ce processus de minimisation peut ne pas converger. Dans ce cas, il n'existe pas de front 20 stationnaire sur toute l'épaisseur du milieu. La méthode décrite doit donc être reprise sur des parties plus petites du réservoir. On peut ainsi trouver des fronts stationnaires dans plusieurs zones du réservoir, qui se déplacent à des vitesses différentes. Exemple On décrit ici un exemple synthétique illustrant la méthode de façon nonlimitative. On considère un réservoir stratifié avec quatre strates horizontales de perméabilités différentes, indiquées sur la figure 4. Pour simplifier au maximum, on négligera les effets de la gravité et on considérera que les perméabilités relatives des strates sont telles que la zone non balayée est saturée en huile seule et qu'en amont de l'interface, seule l'eau est présente à saturation S=1. La porosité est constante dans le réservoir. L'eau est injectée sur la face amont à flux constant connu. L'huile est à viscosité de 1. On commence par chercher l'existence et la forme d'un front stationnaire dans le milieu lorsque la viscosité de l'eau est de 12. Comme on ne dispose pas encore d'a priori, on part d'une interface plane, que l'on déforme de façon itérative jusqu'à atteindre la convergence vers une solution stationnaire (figure 5). Une fois la solution obtenue pour cette viscosité, on peut rechercher la solution stationnaire 10 pour une viscosité de l'eau inférieure, par exemple une viscosité de 10. La solution stationnaire obtenue précédemment peu servir d'a priori en entrée du processus itératif de déformation du front. On converge alors vers une nouvelle solution. En procédant ainsi itérativement sur différentes valeurs de la viscosité de l'eau, on obtient les formes du front stationnaire associées à ces valeurs (figure 6). Comme le rapport des mobilités entre le fluide 15 ne place et le fluide injecté diminue chaque fois que l'on diminue la viscosité de l'eau, le front s'étale au fur et à mesure, car les effets visqueux ont de plus en plus de mal à compenser les différences de perméabilité d'une strate à l'autre. Lorsque la viscosité de l'eau diminue en deçà de la valeur critique 5, on s'aperçoit que le processus de recherche d'un front stationnaire ne converge plus. Les positions du front dans les deuxième et troisième strates s'éloignent 20 perpétuellement l'un de l'autre. On considère alors qu'un front stationnaire ne peut se former sur toute l'épaisseur du milieu. Ce dernier doit être découpé en deux parties, la séparation se situant entre la deuxième et la troisième couche. On poursuit ensuite le processus de recherche de front stationnaire dans chacune des parties précédemment identifiées. En itérant ainsi sur la valeur de la viscosité du fluide injecté, on peut construire un diagramme 25 représentant la réponse du milieu du point de vue hydrodynamique en fonction du rapport des mobilités M entre le fluide en place et le fluide injecté (figure 7). On voit sur ce diagramme que lorsque le rapport des mobilités est supérieur à 5, le milieu se comporte comme une seule couche hydrodynamique, puisqu'un front stationnaire existe sur toute son épaisseur. Lorsque M passe entre 4,5 et 5 deux couches hydrodynamiques apparaissent: la première est formée 30 des deux premières couches "géologiques", et la seconde des deux autres couches géologiques". Puis les deux premières couches géologiques se séparent au sens hydrodynamique en deçà d'un rapport M de 4,5. On a alors trois couches au sens hydrodynamique dans le milieu. Les deux dernières couches géologiques se séparent en deçà de M=3. On a alors quatre couches hydrodynamiques, qui correspondent aux quatre couches géologiques présentes dans le milieu. REVENDICATIONS 1) Méthode pour déterminer dans un milieu stratifié dont les propriétés physiques sont connues ou estimées, au moins une zone o une interface entre un fluide en place dans le milieu et un fluide de balayage, de viscosités et densités différentes connues, injectés dans le 5 milieu, se déplace de façon stationnaire, dans le but de simplifier la construction d'un modèle de simulation des écoulements dans le milieu, caractérisée en ce que: - on considère une forme d'interface a priori en supposant que les déplacements de fluides en tous ses points sont stationnaires; - on détermine le champ de pression de part et d'autre de l'interface a priori; et - on modifie itérativement la forme de l'interface jusqu'à obtenir que les pressions de part et d'autre d'au moins une partie de l'interface s'égalisent en tous points de cette partie; et - on attribue des propriétés hydrodynamiques moyennes uniformément à chaque zone du milieu délimitée par chaque partie d'interface, quand la dite égalisation est atteinte. 2) Méthode selon la revendication 1, caractérisée en ce que, à défaut d'obtenir une égalisation 15 des pressions de part et d'autre de l'interface tout le long de celle-ci, on segmente l'interface en plusieurs parties et on modifie itérativement et séparément la forme de ces différentes parties, jusqu'à obtenir une égalisation des pressions de part et d'autre de chacune d'elles, l'étendue de chaque partie d'interface, quand la dite égalisation est atteinte, délimitant une zone favorable à laquelle on attribue uniformément des propriétés hydrodynamiques moyennes. 3) Méthode selon la revendication 1 ou 2, caractérisée en ce que l'on détermine la forme d'au moins une zone du milieu délimitée par une interface à déplacement stationnaire, qui correspond à des valeurs différentes de la viscosité des fluides de balayage, et on choisit la viscosité pour laquelle on optimise les déplacements stationnaires dans le dit milieu.