Brevet FR2759473A1 — publié le 19980814

L'invention a pour objet une méthode pour simplifier la réalisation d'un modèle permettant de simuler un processus physique dans un milieu matériel et notamment des écoulements de fluides à l'intérieur d'un gisement souterrain. La méthode s'applique par exemple à la modélisation d'écoulements d'hydrocarbures dans des gisements. La caractérisation d'un réservoir pétrolier est un problème délicat auquel sont confrontés les géophysiciens, les géologues et les ingénieurs de gisements. Lorsqu'on cherche à préciser ou à valider la description d'un modèle géologique, I'interprétation de tests de puits, de tests d'interférences et de résultats de production joue un très grand rôle. Certaines de ces techniques d'interprétation sont basées sur des procédures d'inversion où, à partir d'un modèle géologique maillé, on essaye de caler des simulations numériques sur des mesures de terrain en faisant varier les paramètres de ce modèle. Le rôle des simulations numériques est ainsi devenu de plus en plus important et la description géologique devenant de plus en plus fine, la taille des modèles simulés s'accroît et les temps de calculs qui croissent corrélativement, restent très élevés malgré l'amélioration des performances des ordinateurs. De la même façon, lors d'un processus d'inversion utilisant des simulations sur des modèles maillés, il est capital d'obtenir des solutions très précises n'introduisant pas de biais numériques lorsqu'on les compare à des mesures de terrains. Pour modéliser des écoulements dans un gisement pétrolier, il est classique d'utiliser les équations généralisées de Darcy, la loi de conservation de la masse, et de leur associer des conditions aux puits, des conditions aux limites et des conditions initiales. On obtient ainsi un système d'équations non linéaires aux dérivées partielles dont les inconnues sont les pressions P et les saturations S de chaque phase et, éventuellement, les concentrations de chaque constituant dans les phases. où Q est la porosité; Cjp la concentration du constituant i , Sp la saturation de la phase p, , pp la masse volumique et Sqi le débit de constituant injecté ou produit aux puits. La vitesse du fluide et la pression sont reliées par la loi de Darcy où K est la perméabilité absolue, krp la perméabilité relative à la phase p ,up la viscosité relative à la phase p. A ces équations, il convient d'ajouter les bilans de saturations et de constituants, et les relations de capillarité. Dans le cas d'un fluide monophasique et peu compressible, le système d'équations se réduit à K cvx a - Vdiv(- grad(P) +pgz) = q1 > (3) It avec des conditions de flux nul ou de pression imposée aux bords et des conditions initiales. C représente la compressibilité totale (pores et fluide), supposée constante dans le modèle, et VO le volume poreux. Après discrétisation du milieu souterrain étudié par un maillage, en utilisant une méthode de résolution par volumes finis, le système d'équations à résoudre pour chaque maille du modèle maillé est avec: i : indice de maille v(i) : I'ensemble des mailles voisines de la maille i = = I : si le puits p est dans la maille i T"(" : la transmissivité entre les mailles i et v(i) Si on écrit le système sous forme matricielle en notant A la matrice des transmissivités, D, la matrice diagonale des CV et Q, le second membre, on obtient où P est un vecteur de dimension n le nombre d'inconnues qui sont les pressions aux puits et dans chaque maille du modèle maillé, et Po représente l'état initial du gisement. La matrice-A est une matrice carrée nxn symétrique définie positive. La matrice diagonale D étant positive et indépendante du temps, on peut écrire le système matriciel sous la forme suivante La résolution pour toutes les mailles du système d'équations s'avère très lourde et très coûteuse en temps de calcul si le nombre de mailles du modèle est très élevé, comme c'est généralement le cas. La méthode selon l'invention permet de simplifier la réalisation d'un modèle permettant de simuler des écoulements de fluides dans un milieu matériel hétérogène, (tel que des écoulements d'hydrocarbures dans des gisements par exemple), que l'on discrétise par un maillage en un grand nombre n de cellules élémentaires ou mailles, ce modèle étant formé par la résolution d'un système d'équations de modélisation des écoulements de fluides, à partir de valeurs initiales connues d'un certain nombre de paramètres physiques tels que des pressions, des saturations de différentes phases, des débits, des porosités etc., valeurs initiales obtenues par différents moyens connus: carottages, logs, interprétations, etc. On réalise ce modèle pour un coût relativement faible et dans un temps de calcul raisonnable, si on les rapporte à la qualité et la précision des simulations obtenues, essentiellement par une réduction de la taille du problème de résolution, en utilisant un espace de Krylov de dimension d très inférieure au nombre n de mailles. La méthode selon l'invention s'applique en particulier pour constituer un modèle pour simuler des écoulements de fluides monophasiques. D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après de modes de réalisation décrits à titre d'exemples non limitatifs, en se référant aux dessins annexés où - la Fig. I montre les évolutions comparées de la dérivée par rapport au logarithme du temps de la pression au niveau d'un puits de production foré au travers d'un gisement d'hydrocarbures, obtenues par la méthode selon l'invention et une méthode de simulation classique dans le cas d'un gisement homogène monocouche d'épaisseur variable; - la Fig. 2 montre les évolutions comparées de la dérivée de la pression au puits par rapport au logarithme du temps, obtenues par la méthode selon l'invention et une méthode de simulation classique, dans le cas d'un gisement hétérogène monocouche organisé en chenal; et - la Fig. 3 montre les évolutions comparées de la dérivée de la pression par rapport au logarithme du temps, obtenues par la méthode selon l'invention et une méthode de simulation classique, dans le cas d'un gisement hétérogène multicouches. I) Description de la méthode L'objectif de la méthode est d'obtenir beaucoup plus rapidement un modèle permettant la simulation d'écoulements (monophasiques et polyphasiques) dans un milieu souterrain, ceci en réduisant la taille d'un problème complexe de modélisation de dimension n considérable (n est de l'ordre de 106) à un système d'équations différentielles de dimension d petite (d de l'ordre de 10 par exemple). Elle permet aussi de déterminer les gradients de la pression par rapport à des paramètres du milieu. Comme on le verra dans la description ci-après, le modèle est obtenu en extrayant de l'espace de dimension n, un sous-espace pertinent de dimension d. Ce sous-espace est un espace de Krylov dont la construction est basée sur le calcul de dérivées d'ordre supérieur en temps, ou sur une autre technique (sous-espaces propres, vecteurs de bases calculés avec un schéma d'intégration du problème de taille n, etc.). La méthode s'applique à la dérivation par rapport au temps et à la dérivation par rapport à des paramètres du milieu. La dérivation par rapport au temps permet de résoudre le problème direct et la dérivation par rapport aux paramètres permet de traiter le problème inverse. Le problème à résoudre est la résolution d'un système de la forme D. -A.U(t)+Q D. #t U(O) Uo et des conditions aux limites où (Ui )i== 1,n est un vecteur obtenu par discrétisation de l'espace étudié en un maillage de n éléments, D est une matrice diagonale, A une matrice symétrique et Q un vecteur de dimension n. La simplification apportée par la méthode repose sur l'utilisation de l'espace dit de Krylov, décrit dans la littérature, notamment par - Sonneveld, P., CSG: A Fast Lanczos-type Solver for Nonsymmetric Liner Systems, SIAM J. Sci. Comp, Vol 10, No 1, 36-52, 1989. - Saad, Y., et al.; GMRES: A generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, Siam J. Sci Stat. Comp., vol 7, No 3, 856-869, 1986; ou - Golub, G.H. et al; Matrix Computations, The Johns Hopkins University Press, Baltimore, Second Edition, 1989. On rappelle qu'un espace de Krylov est un espace vectoriel de dimension d constitué par la suite de vecteurs Eô La matrice A étant creuse, ces produits matrice-vecteur ont une complexité linéaire en fonction de n . La simplification apportée par la méthode passe par la projection du système d'équations dans un espace de Krylov de dimension d+l, avec d beaucoup plus petit que n, engendré par la suite de vecteurs (u(O), U (1)..., u(d)) où U(i) = AU('-' et U(()) = U(0) On cherche ensuite une solution approchée Ua(t) dans cet espace de Krylov sous la forme Le problème consiste maintenant à calculer les (vj(t)) Od ce qui est un problème de dimension réduite. II est bien connu que la méthode du gradient conjugué donne à la dième itération la solution qui minimise le résidu du système à résoudre dans un espace de Krylov de dimension d. D'où la convergence en n itérations au plus. On n'utilise que des vecteurs de Krylov indépendants en considérant le plus grand indice d l'espace Ed, et donc que l'on peut travailler dans cet espace. Dans la pratique, on va utiliser un espace de Krylov particulier en raison de la présence du terme source Q que l'on cherche à éliminer. A cet effet, on écrit le système d'équations (6) sous la forme suivante et on en cherche la solution approchée dans l'espace de Krylov E, = #yl. Yd sous la forme A cet effet, on utilise par exemple la méthode de Galerkin décrite notamment dans: - Lascaux, P et Théodor, R ; Analyse numérique matricielle appliquée à l'art de l'ingénieur, Masson, Paris, 1987, Tome 2. - Strang, G ; Introduction to Applied Mathematics, Wellesley - Cambridge Press, Wellesley, Massachusetts, 1986. La méthode de Galerkin conduit à: az,, (r) Az,(t)y, > ,lild et en tenant compte de la relation (8), on obtient En posant M = ( )i.,=I.J et B = ( )i ,= le système s'écrit: La condition initiale découle de la relation (8) en écrivant que zu(0)= z(O) et en utilisant le fait que les vecteurs de Krylov sont indépendants. Le système d'équations (9) de dimension d peut se résoudre facilement et on obtient les expressions des vi en fonction du temps. On en déduit ensuitc une solution approchée du système d'équations (6) en calculant II) Résolution du problème réduit Lorsque ci est petit, on peut résoudre le système (9) par une méthode directe. On considère une matrice S telle que M = 5Y5,' (S est une racine carré de M) et on pose w(t) = Sy2v(t) . Le système d'équations (9) devient alors et on peut le résoudre par diagonalisation de la matrice symétrique S-y-2BS-y2. Cette méthode d'intégration "exacte" ne peut pas être utilisée dans le cas où le système d'équations (5) a un second membre qui varie en fonction du temps et ni dans le cas non linéaire. Dans le cas général, on peut utiliser un schéma numérique pour l'intégration du problème réduit. III) Mise en oeuvre numérique En rappelant que M = ( )ij,=,1d, , nous avons, en utilisant la symétrie de la matrice des transmissivités A , = = > . Cela veut dire que les termes de M sont constants sur la codiagonale i+j = constante. On montre également que B est obtenue par translation de la matrice M b, bi,j = bij=m1+1. Ces remarques sont de nature à accélérer les calculs. Mais dans la pratique, les itérés successifs de Ai ont tendance à augmenter fortement en norme. Une démarche possible, consiste à appliquer le procédé d'orthogonalisation de Schmidt aux vecteurs' de la base de Krylov. Cette technique présente deux avantages elle évite la génération d'éléments de la base de norme très grande, elle permet de détecter l'existence de vecteurs de Krylov dépendants et ainsi de déterminer la taille optimale d du petit système. IV) Résultats La méthode a été appliquée à quelques cas de référence en comparant les solutions obtenues aux solutions données par une simulation classique consistant à discrétiser le temps en une suite d'instant tk et à résoudre le système de grande taille à chacun de ses instants. 1) Exemple 1: Cas homogène On considère le cas d'un gisement d'hydrocarbures comportant une couche de perméabilité et porosité homogènes, sur laquelle on a formé un maillage comprenant 58x67 mailles horizontalement. Le maillage est plus fin à proximité du puits qui est placé au centre du gisement modélisé. La perméabilité est de 700 mD et la porosité est de 30%. L'épaisseur du gisement modélisé est variable et comprise entre 40 m et lm et la taille des mailles varie de 50m loin du puits à Sm autour du puits. On a simulé un pompage durant 86400 s (1 jour) avec un débit de 30 m3/jour. On constate sur la Fig. 1 que les valeurs obtenues par la méthode selon l'invention sont semblables aux valeurs obtenues par une méthode classique du type discrétisation du temps et résolution des systèmes linéaires par une méthode de gradients conjugués. Mais les résultats obtenus par la méthode proposée coïncident avec le gradient conjugué lorsque le pas de discrétisation en temps est petit en début de simulation (0,1) et que sa progression géométrique est faible (1,01). Lorsque les pas de temps augmentent avec une raison géométrique plus importante, les résultats obtenus en dérivée divergent. Les pas de temps initiaux et maximaux ainsi que la raison de la progression géométrique des différents tests sont présentés dans le tableau ci-après. A qualité de simulation égale, les temps de calcul nécessaires à la construction du modèle par la méthode selon l'invention sont beaucoup plus faibles que ceux nécessaires à une simulation classique par une méthode du type gradients conjugués. Toutes les simulations présentées ont été réalisées sur une station de travail. Cas Pas de temps Progression du Pas de temps Temps de calcul Min pas de temps Max. (secondes CPU) 1 I 1,3 10000 20 2 I 1,1 1000 68 3 0,1 1,01 1000 268 Méthode Base de Krylov de 25 éléments 85 proposee Tableau 1: Comparaison des temps de calculs du cas homogène d'épaisseur variable 2) Exemple 2: Cas d'un chenal On considère le cas d'un gisement d'hydrocarbures formé d'une couche comportant 2 perméabilités différentes organisées en chenal et de porosité homogène. Le maillage comprend 91x91 mailles horizontalement toutes identiques de 10m sur 10m. Le puits est placé au centre du gisement modélisé et à l'intérieur du chenal de perméabilité est de 100 mD et de porosité 20%. Le chenal est parallèle à l'axe des X et de largeur 310m centré sur le puits. L'épaisseur du gisement modélisé est constante est égale à 10 m. La partie hors chenal a une perméabilité de ImD et une porosité de 20cil. On a simulé un pompage durant 300.000 s (environ 3,5 jours) avec un débit de 50 m3/jour. On constate sur la Fig.2 que, comme dans le cas précédent, les valeurs obtenues par la méthode selon l'invention sont semblables aux valeurs obtenues par une méthode classique du type discrétisation du temps et résolution des systèmes linéaires par une méthode de gradients conjugués. Mais les résultats obtenus par la méthode proposée coïncident avec le gradient conjugué lorsque le pas de discrétisation en temps est petit en début de simulation (0,1) et que sa progression géométrique est faible (1,01). Lorsque les pas de temps augmentent avec une raison géométrique plus importante, les résultats obtenus en dérivée divergent. Les pas de temps initiaux et maximaux ainsi que la raison de la progression géométrique des différents tests sont présentés dans le tableau ci-après. Il est à noter que dans toutes les configurations testées la méthode proposée est plus rapide que la simulation classique et ceci avec des résultats plus précis. Cas Pas de temps Progression du Pas de temps Temps de calcul Min pas de temps Max. (secondes CPU) 1,3 10000 31 2 l 1,1 1000 71 3 0,1 1,01 1000 389 Méthode Base Base de Krylov de 20 éléments 25 proposée 3) Exemple 3: Cas hétérogène multicouche On considère le cas d'un gisement d'hydrocarbures comportant 3 couches formées subdivisée chacune par un maillage de 55 mailles dans les 2 directions horizontales. Le maillage est plus fin près du puits placé au centre du gisement modélisé. La perméabilité horizontale des trois couches est de bas en haut respectivement: l OmD, 200mD et lmD. Les perméabilités verticales valent respectivement lmD, 5mD, 0.imD ct la porosité est homogène (20nô). L'épaisseur de chaque couche est de 5m et la taille horizontale des mailles varie de 100m loin du puits à Sm autour du puits. On a simulé un pompage durant 200 000 s avec un débit de 50m3/jouI-. Comme dans les exemples précédents, on vérifie que les valeurs obtenues par la méthode selon l'invention sont semblables aux valeurs obtenues par une méthode classique. En se référant à la Fig.3, on note que les résultats obtenus par la méthode proposée, coïncident avec le gradient conjugué lorsque le pas de discrétisation en temps est petit au départ (0.1) et augmente lentement avec une raison de 1.01. Les solutions obtenues différent de plus en plus lorsque le pas de discrétisation en temps augmente avec une raison géométrique de plus en plus grande. Les caractéristiques des tests présentés et les temps de calculs respectifs sont regroupés dans le tableau ciapyres Cas Pas de temps Progression du Pas de temps Temps de calcul Min pas de temps Max. (secondes CPU) 1 I 1,3 l0000 74 2 1 1,1 1000 347 3 0,1 1,01 1000 890 Méthode proposée Base de Krylov de 25 éléments I 85 Les trois exemples présentés ci-dessus parmi tous ceux qui ont permis de valider la méthode de modélisation proposée avec projection dans un espace de Krylov, illustrent bien ses avantages. Le problème d'évolution de grande taille est réduit à la résolution d'un système d'équations différentielles de faible dimension, et les résultats numériques obtenus sont plus précis qu'avec une simulation classique pour un coût moindre. Dans l'optique d'une utilisation des simulations numériques pour interpréter les essais de puits, l'utilisation des dérivées d'ordre supérieur en temps pour construire la base de Krylov présente un avantage très important car elle permet la simulation directe de la dérivée de la pression par rapport au temps. Cette dérivée est disponible pour la pression au puits mais également pour chaque maille modélisant le gisement et peut être calculée de façon continue à chaque instant contrairement au cas discret. Il faut noter enfin que le coût de stockage des informations nécessaires en cours de calcul est identique à celui d'une simulation classique avec discrétisation du temps et résolution d'un système linéaire à chaque pas de temps. REVENDICATIONS 1) Méthode pour simplifier la réalisation d'un modèle permettant de simuler des écoulements de fluides dans un milieu matériel hétérogène, que l'on discrétise par un maillage en un grand nombre n de cellules élémentaires ou mailles, le modèle étant formé par la résolution dans chacune des n mailles, d'un système d'équations de modélisation des écoulements de fluides, à partir de valeurs initiales connues d'un certain nombre de paramètres physiques, caractérisée en ce que l'on simplifie la construction du modèle en utilisant un espace de Krylov de dimension d très inférieure au nombre n de mailles pour réduire le volume des opérations de résolution. 2) Méthode selon la revendication 2, caractérisée en ce que l'on modélise la circulation de fluides monophasiques. 3) Méthode selon l'une des revendications 1 ou 2, caractérisée en ce que pour générer l'espace de Krylov, on utilise des dérivées d'ordre supérieur en temps. 4) Application de la méthode selon l'une des revendications précédentes à la réalisation d'un modèle de simulation d'écoulements d'hydrocarbures à l'intérieur d'un gisement souterrain.