FR 2850773 A1 20040806 FR 0301090 A 20030131 La présente invention concerne une méthode pour modéliser des flux de fluides fortement compressibles tels que les gaz dans un milieu poreux multicouches traversé par un réseau de fractures de géométrie donnée, inégalement réparties dans le milieu, et dont certaines communiquent les unes avec les autres. La méthode selon l'invention trouve des applications notamment dans le domaine 10 de la production d'hydrocarbures. Elle permet aux ingénieurs de réservoir de tester par simulation la mise en production d'un gisement par un ou plusieurs puits traversant une zone souterraine poreuse perméable comportant deux milieux contrastés, un milieu matriciel contenant généralement la plus grande part du fluide en place et présentant une faible perméabilité, et, le traversant en partie ou en totalité, des réseaux de fractures 15 beaucoup plus conductrices dans la plupart des cas. Etat de la technique Lors d'un essai de puits à gaz, les conditions de débit imposées au puits amènent le gaz contenu dans le réservoir à s'écouler vers le puits. Il s'agit d'un écoulement monophasique, mono constituant (la seule phase mobile est le gaz) et fortement 20 compressible (la roche du réservoir et le fluide sont compressibles). Pour tout volume élémentaire du réservoir, la pression du gaz contenue dans ce volume obéit alors à l'équation suivante: D CTa= ( VK) + P Q (1) avec: 2: volume poreux CT: compressibilité totale (fluide + roche) K: perméabilité de la roche g: viscosité du fluide Z: facteur de compressibilité du gaz Q: débit entrant, et iv: fonction pseudopression définie par yv = 2| P -dp o p est la pression du fluide Po La viscosité et le facteur de compressibilité du gaz sont variables en fonction de la pression et sont donnés pour un certain nombre de valeurs dans un table dite PVT (Pression Volume Température). A partir de cette table nous pouvons également déduire la fonction 10 pseudo-pression définie cidessus et la compressibilité du gaz déduite de l'équation d'état. La table PVT est donc un tableau à 5 colonnes (pression, pseudo-pression, viscosité, facteur de compressibilité et compressibilité du gaz) à partir duquel, à une pression donnée on peut obtenir les pseudo-pression, viscosité et compressibilité correspondantes par interpolation linéaire (inversement, à partir d'une pseudo-pression on peut obtenir par 15 interpolation les 4 autres données). Le débit entrant Q est nul partout sauf aux endroits o le puits est en communication avec le réservoir. Pour simuler un essai de puits, dans quelque milieu que ce soit, on doit résoudre cette équation dans l'espace et dans le temps. On réalise donc une discrétisation du 20 réservoir (maillage) et la résolution du problème consiste à trouver les pressions des mailles au cours du temps, lui-même discrétisé en un certain nombre de pas de temps. On connaît des outils de modélisation d'écoulements monophasiques qui sont appliqués non au milieu géologique réel dans toute sa complexité mais à une représentation homogénéisée, suivant le modèle de réservoir dit de double milieu décrit par exemple par 25 Warren et Root, dans " The Behavior of Naturally Fractured Reservoirs ", SPE Journal Septembre 1963. Tout volume élémentaire du réservoir fissuré est ainsi modélisé sous la forme d'un ensemble de blocs parallélépipédiques identiques limités par un système orthogonal de fractures uniformes continues orientées dans la direction de l'un des trois principaux sens d'écoulement (modèle dit de "boite à sucre"). L'écoulement des fluides à l'échelle du réservoir s'effectue principalement à travers le milieu fissure et des échanges de 5 fluides surviennent localement entre les fissures et les blocs matriciels. Cette représentation qui ne rend pas compte de la complexité du réseau de fractures dans un réservoir, se révèle efficace néanmoins mais à l'échelle d'une maille de réservoir (typiquement 100m x 100m). Cependant il est préférable que les ingénieurs réservoir disposent d'un simulateur d'écoulement basé sur un modèle géologique " réel " du milieu et non sur un modèle 10 homogène équivalent, cet outil permettant de: * valider l'image géologique du réservoir construite par le géologue à partir de l'ensemble des informations qu'il a pu recueillir sur la fracturation du réservoir (cette validation s'effectuant par comparaison aux données réelles d'essais de puits); * tester de manière fiable la sensibilité du comportement hydraulique du milieu 15 aux incertitudes sur l'image géologique du milieu fissuré. Par les brevets FR 2 787 219 du demandeur, on connaît une méthode pour modéliser les flux de fluides faiblement compressibles (huile) dans un milieu poreux multicouches fissuré en tenant compte de la géométrie réelle du réseau de fractures et des échanges locaux au sein de la matrice poreuse et entre la matrice poreuse et les fractures en 20 chaque noeud du réseau. On discrétise le milieu fissuré par un maillage, on associe les mailles de fracture à des noeuds placés soit aux intersections des fractures soit aux extrémités des fractures. Chaque noeud est associé à un volume ou bloc de matrice, on détermine les flux entre chaque maille de fracture et le volume de matrice associé dans un régime de flux semi-permanent. Dans le cas ou le réseau de fracture est faiblement 25 connecté, c'est à dire pour lequel le réseau de fractures ne traverse pas l'ensemble du volume matrice considéré, on détermine aussi les flux directs entre les volumes de matrice au travers des arêtes communes des mailles. En tenant compte des transmissivités entre blocs matriciels associés à des noeuds de discrétisation du milieu, on peut simuler de façon plus réaliste la réponse d'un puits à des variations de débit imposées. La méthode selon l'invention a pour objet de modéliser des flux de fluides fortement compressibles dans un milieu poreux multicouches traversé par un réseau dé fractures de géométrie donnée et par un puits. Elle comprend les étapes suivantes: a) on discrétise chaque couche fracturée par un maillage comprenant des mailles de 5 fracture que l'on centre sur des noeuds soit aux intersections de fractures soit aux extrémités des fractures, chaque noeud étant associé à un bloc de matrice regroupant tous les points plus proches de lui que de noeuds voisins; b) on calcule les flux locaux entre chaque maille de fracture et le volume de matrice associé dans un régime de flux pseudo permanent, la valeur de la transmissivité étant 10 obtenue en considérant une variation linéaire de la pression en fonction de la distance entre tout point du bloc et la maille de fracture; c) on détermine les flux directs entre les mailles de fracture; d) on détermine les flux directs entre les volumes de matrice au travers des arêtes communes des mailles; e) on détermine les flux directs entre d'une part les mailles de fracture et les mailles de matrice, et d'autre part le volume du puits, à partir de valeurs de transmissivité adaptées à des écoulements turbulents; et f) on simule la réponse du puits pour des variations de débit imposées au moyen d'une équation de diffusion modifiée pour tenir compte de la compressibilité du fluide, qui 20 relie les interactions entre les variations de la pression et du débit observables dans le puits. Suivant un mode de mise en oeuvre, on modifie l'équation de diffusion par introduction d'un terme dépendant de la pression, de la viscosité et d'un facteur de compressibilité du fluide compressible. Suivant un mode de mise en oeuvre, on détermine les flux directs en introduisant un 25 terme d'écart proportionnel au débit de fluides imposé. Suivant un mode de mise en oeuvre, on détermine le bloc de matrice associé à chaque maille de fracture en discrétisant chaque couche fracturée en un ensemble de pixels et en calculant la distance de chaque pixel à la maille de fracture la plus proche, pour déterminer l'emplacement des arêtes entre les mailles et l'on utilise les distances calculées pour en déduire la valeur des transmissivités entre chaque maille de fracture et le volume de matrice associé d'une part, et entre les arêtes communes des mailles. On repère par exemple la position des arêtes entre les mailles en déplaçant deux 5 fenêtres allongées dans l'image formée par les pixels, orientées suivant deux directions perpendiculaires. Ainsi, en tenant compte de la forte compressibilité du gaz et adaptant les transmissivités entre les blocs matriciels et le puits et entre les fissures et le puits on peut simuler de façon réaliste la réponse d'un puits à gaz à des variations de débit imposées. Présentation sommaire des figures D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après, en se référant aux dessins annexés o: - la figure 1 montre un exemple de maille de fracture en 2D; - la figure 2 montre un exemple de bloc matriciel obtenu en 2D; la figure 3 montre un exemple de deux blocs matriciels MB communiquant par l'intermédiaire d'un élément de fracture FE; - la figure 4 montre un exemple de deux blocs matriciels MB traversés chacun par un élément de fracture et communiquant par le biais d'échanges au travers de leurs arêtes communes A; - la figure 5 montre des exemples de fenêtres orientées suivant deux axes perpendiculaires, permettant le repérage des arêtes entre les blocs matriciels; - la figure 6 montre un puits W intersecté par deux fissures; - la figure 7 montre un exemple de variation en fonction du temps t de débit D imposé dans un essai de puits; - la figure 8 montre deux types d'écoulement linéaire (LF) et radial (RF) de la matrice vers le puits; - la figure 9 montre une première partie de l'organigramme de mise en oeuvre de la modélisation modifié dans le cas considéré o les fluides sont fortement compressibles et - la figure 10 montre une partie complémentaire de l'organigramme de mise en oeuvre de 5 la modélisation modifié dans le cas considéré o les fluides sont fortement compressibles. DESCRIPTION DETAILLEE Dans le cas d'un milieu fissuré, si l'on veut tenir compte de toute la complexité du réseau de fractures, il est nécessaire de réaliser un maillage explicite de ce réseau. Par 10 ailleurs, dans un tel milieu, les perméabilités K sont en général beaucoup plus élevées dans les fractures que dans la matrice et les écoulements sont par conséquents rapides dans les fractures et lents dans la matrice. Il s'agit donc de représenter le milieu fissuré à l'aide du concept de double milieu (matrice et fissure) décrit dans le brevet FR 2 787 219 précité. Les étapes de l'approche adoptée dans le cas o le fluide produit est de l'huile sont les 15 suivantes: * Maillage explicite du réseau de fractures * association à chaque maille de fracture d'un bloc matriciel unique; * traitement des écoulements entre chaque maille de fracture et son bloc associé en régime pseudo permanent o le flux de l'un à l'autre est proportionnel à la différence de 20 pression entre les deux; * prise en compte des écoulements entre blocs matriciels; et * simulation dynamique des écoulements au sein du réseau de fractures et du milieu matriciel. Dans le cas o le fluide produit est fortement compressible, typiquement un gaz, 25 nous précisons ou modifions: * le système d'équations de modélisation des écoulements pour retrouver une équation différentielle de diffusion de même forme que pour les écoulements faiblement compressibles; * les termes décrivant les flux de fluides entre les mailles puits et la maille de fracture adjacente ainsi qu'entre les mailles puits et la maille matrice adjacente; et * la simulation dynamique pour réactualiser au cours de la simulation les paramètres de viscosité et de compressibilité totale à chaque itération en temps par une technique d'interpolation, à partir d'une table de données dites de PVT. Les inconnues sont les pressions des mailles de fracture et les pressions des blocs 10 matriciels associés. Puisqu'il y a autant de blocs matriciels que de mailles de fracture, le nombre d'inconnues est égal à deux fois ce nombre. * L'équation de diffusion L'équation de diffusion que nous résolvons est déduite de l'équation de conservation de la masse, de la loi de Darcy et d'une équation d'état qui dépend de la compressibilité du 15 fluide mobile. Nous considérons deux équations d'état différentes selon que le fluide est faiblement compressible (huile, cas a) ou fortement compressible (gaz, cas b). a. Pour lesfluides faiblement compressibles (huile) L'équation d'état traduisant la compressibilité équivalente du fluide mobile (l'huile) est la suivante: Ch= ap (2) PP)T Considérant le système composé des trois équations ou loi que nous venons de citer (équations de conservation de la masse et d'état, loi de Darcy), pour tout volume élémentaire du réservoir, la pression de l'huile contenue dans ce volume obéit alors à l'équation différentielle suivante: ( *CT P = div(- .Vp) +Q (3) ot t avec: ): porosité CT: compressibilité totale (fluide + roche) K: perméabilité de la roche gt: viscosité du fluide Q: débit entrant p: densité du fluide p: inconnue pression 10 b. Pour les fluides fortement compressibles (gaz) Dans le cas du gaz, la densité qui varie considérablement en fonction de la pression s'exprime comme suit: pM (4) ZRT avec M masse molaire du gaz, R constante des gaz parfaits et Z le facteur de compressibilité du gaz qui traduit l'écart du gaz réel par rapport aux gaz parfaits. La compressibilité du gaz peut alors s'écrire: cg P Za (5) P Z[ 0 PiT L'équation différentielle de diffusion devient alors: 20 div kgradp + ICt Z - = q'z (6) / LZ tZ pto Z avec C, la compressibilité globale d'un élément unitaire du volume poreux saturé. Pour retrouver la forme d'une équation de diffusion plus classique; cette équation est exprimée à l'aide d'une fonction de pseudo-pression définie par: y =:2{ -Ldp (7) Po t L'équation de diffusion exprimée en pseudopression s'écrit finalement: k V=2 p q (8) L'équation différentielle à résoudre présente alors une forme similaire à celle utilisée dans le cas des fluides faiblement compressibles, l'inconnue étant la pseudopression Y- (reliée à la pression par la table PVT). * Discrétisation du réseau de fractures En 2D, les noeuds de calcul sont positionnés aux intersections entre fractures et aux extrémités des fractures. Les besoins de la modélisation 3D conduisent à rajouter des noeuds supplémentaires dans les couches supérieures et inférieures pour les fractures traversant plusieurs couches. Dans le maillage de la présente méthode, les noeuds de calcul ainsi définis constituent les centres des mailles de fracture. De plus, la simulation d'écoulements fortement compressibles nécessite de connaître le volume des mailles (p dans l'équation 8). On définit donc des limites de mailles positionnées au milieu des segments reliant les noeuds de calcul. Le schéma de la Fig. i montre un exemple de maille de fracture en 2D. En 20 terme de nombre de noeuds de calcul, le maillage du réseau de fractures est ainsi rendu optimal. * Transmissivités fissure-fissure Le calcul des liaisons entre mailles de fracture (transmissivités), utilisées pour le calcul des flux entre ces mailles de fracture, est identique à celui présenté dans la méthode 25 décrite dans le brevet FR 2 757 947 précité. * Discrétisation du milieu matriciel Suivant le principe de la méthode, la discrétisation du milieu matriciel consiste à affecter un volume de roche à chaque maille de fracture définie comme indiqué dans le paragraphe précédent. Lors de la simulation dynamique, les échanges matrice-fracture 5 seront calculés en pseudo permanent entre chaque maille de fracture et son unique bloc matriciel associé. Pour le calcul de ces volumes de matrice, on traite le problème couche par couche. Pour une couche donnée, on définit les blocs de la manière suivante: * les hauteurs de blocs sont égales et fixées à la hauteur de la couche; * dans le plan de la couche, la surface du bloc associé à une maille de fracture est l'ensemble des points plus proches de cette maille de fracture que d'une autre. Physiquement, cette définition suppose que les limites à flux nul dans la matrice sont les lignes d'équidistance entre les mailles de fracture. Pour déterminer la géométrie des blocs dans une couche donnée, il s'agit donc de 15 résoudre un problème bidimensionnel qui consiste à trouver, pour chaque maille de fracture les points les plus proches de cette maille que d'une autre. Ce problème est résolu en utilisant la méthode géométrique exposée dans le brevet FR 2047 957 précité. La méthode en question permet, en discrétisant la couche fracturée en un ensemble de pixels et en appliquant un algorithme de traitement d'image, de déterminer la distance de 20 chaque pixel à la fracture la plus proche. On rappel que dans la phase d'initialisation de cet algorithme, on affecte la valeur 0 (distance nulle) aux pixels appartenant à une fracture et une grande valeur aux autres. Si dans cette phase, on donne de plus à chaque pixel fracture le numéro de maille de fracture auquel il appartient, le même algorithme permet finalement de déterminer, pour chaque pixel: * la distance de ce pixel à la maille de fracture la plus proche; * le numéro de la maille de fracture la plus proche. L'ensemble des pixels ayant le même numéro de maille de fracture constitue la surface du bloc matriciel associé à cette maille. En multipliant cette surface par la hauteur de la couche, on obtient le volume du bloc associé à la maille. Un exemple de bloc matriciel ainsi obtenu en 2D est schématisé à la Fig.2. Par construction, la somme des surfaces des blocs est égale à la surface totale de la couche, ce qui garantit la conservation du volume de la couche. * Transmissivités matrice-rissure Pour chaque bloc matriciel, l'algorithme de traitement d'image donne aussi la distance de chaque pixel du bloc à la maille de fracture associée. On se sert de cette 10 information pour calculer la transmissivité entre la maille de fracture et le bloc matriciel. L'hypothèse de régime pseudo permanent entre une maille de fracture et un bloc matriciel revient à considérer que le flux de l'un à l'autre est proportionnel à la différence des pressions de l'un et l'autre. On a la relation suivante: Fmf = -(Pm-pf) (9) avec: Fn: flux matrice-fracture Tmf: transmissivité matrice-fracture g: viscosité Pm: pression du bloc matriciel pf: pression de la maille de fracture associée Dans l'hypothèse du régime pseudo permanent, la valeur de la transmissivité Tmf est constante au cours de la simulation. Dans la présente méthode, cette valeur est calculée pour chacun des couples maille de fracture - bloc matriciel. Pour ce calcul, on fait l'hypothèse que dans le bloc matriciel, la pression varie linéairement en fonction de la distance du point considéré à la maille de fracture associée au bloc. On définit la transmissivité comme suit 21 H- K Tmf D (10) avec: l: longueur des fractures dans la maille de fracture H: hauteur de la couche (et du bloc matriciel) K: perméabilité de la matrice D: distance des fractures à laquelle la pression est la pression moyenne du bloc On connaît 1, H et K (le facteur 2 vient du fait que la fracture à deux faces). Le calcul de D est fait à partir des informations données par l'algorithme de traitement d'image sur les distances des pixels aux fractures les plus proches. En effet, suivant l'hypothèse de variation linéaire de la pression en fonction de la distance à la fracture, on a: 1 -(1 D = Jd(s) ds (11) o S est, en 2D, la surface du bloc matriciel et d(s) est la fonction donnant la distance aux fractures des points du bloc matriciel à la maille de fracture associée. En terme de pixels, on a donc simplement: D= cIdn (12) N N o N est le nombre de pixels du bloc matriciel et dj la distance du pixel n à la maille de 20 fracture. * Transmissivité matrice-matrice Le calcul des échanges entre blocs matriciels repose sur un schéma numérique d'approximation des flux à 2 points. Cette approche suppose que le flux entre deux blocs matriciels est perpendiculaire à l'arête séparant ces deux blocs. Ce schéma numérique est 5 donc généralement utilisé pour les maillages orthogonaux comme les maillages de Vorono. Dans notre problème, la géométrie des blocs matriciels est quelconque. De plus, la présence des fractures amène à considérer deux cas géométriques possibles. Dans le premier cas (Fig.3), l'arête entre les deux blocs matriciels est traversée par 10 une fracture. Dans le deuxième cas, elle ne l'est pas. Dans le deuxième cas (Fig.4), on fait l'hypothèse que dans chacune des deux fractures du problème, la pression peut être considérée comme constante devant la variation de pression dans le bloc matriciel associé (car la conductivité de la fracture est très grande devant la perméabilité de la matrice). On calcule une transmissivité à deux 15 points entre les deux segments de fractures à travers l'arête séparant les deux blocs matriciels. Ce calcul est fait par un algorithme de traitement d'image qui utilise l'image numérique définie lors de la détermination de la géométrie des blocs matriciels et qui renseigne donc, pour chaque pixel de l'image, sur la distance de ce pixel à la fracture la 20 plus proche et sur le bloc matriciel auquel appartient ce pixel. L'algorithme de traitement d'image consiste à repérer les interfaces élémentaires entre blocs, c'est à dire les arêtes séparant deux pixels associés à deux blocs différents. L'algorithme calcule, pour chacune de ces interfaces, une transmissivité élémentaire entre les deux blocs et met à jour la transmissivité totale entre ces blocs suivant la formule 25 suivante: T =k M H p2-(di -dfn) (13) 2 - 2di, o Tmm est la transmissivité matrice-matrice entre deux blocs matriciels, M le nombre d'interfaces élémentaires entre les deux blocs, p la taille du coté des pixels et di, et df1 les distances des deux extrémités de l'interface n aux fractures les plus proches. Comme on le voit, la transmissivité totale est une somme de transmissivités élémentaires calculées pour chacune des interfaces entre les deux blocs. La Fig.5 montre un ensemble de pixels appartenant aux mailles notées 1, 2 et 3. Pour repérer les interfaces entre les pixels de ces mailles, l'algorithme de traitement d'image déplace deux fenêtres constituées de 2 fois 1 pixel le long des lignes et des colonnes de l'image. Il y a une fenêtre verticale pour repérer les interfaces horizontales et 10 une fenêtre horizontale pour repérer les interfaces verticales. Dans le cas o une fracture traverse l'arête entre deux blocs matriciels, la transmissivité entre ces deux blocs est calculée simplement avec l'expression suivante: T. = .km H* L (14) LAB o I est la longueur de l'arête et LAB la longueur de la fracture reliant les centres A et B. La 15 longueur de fracture est connue par construction et la longueur de l'arête est calculée par l'algorithme de traitement d'image à partir des interfaces élémentaires entre les blocs matriciels de centres A et B. * Représentation des puits Un puits est représenté par sa géométrie d'une part et les contraintes de débit qui lui 20 sont imposées au cours du temps d'autre part. Un puits est constitué par un ensemble de segments connectés qui intersectent les fractures du réseau (Fig. 6). La représentation géométrique d'un puits est donc une ligne brisée en 3D. Les points de communication entre le puits et le réseau sont les points d'intersection de cette ligne avec les fractures. Les contraintes de débit sont imposées au point o le puits entre dans le milieu fracturé. Elles sont entrées par l'utilisateur sous la forme d'une fonction en plateau donnant le débit en fonction du temps. Les temps de changement de débit du puits indiquent les différentes périodes à simuler. Pour un essai de puits classique, on impose un débit pendant une première période au terme de laquelle on ferme le puits (débit nul) et l'on observe la remontée de pression dans le puits. La courbe de débit à donner est schématisée sur la Fig. 7. * Liaison puits-réservoir a) Transmissivité fissure-puits Le calcul des échanges entre le puits et une fracture repose sur l'hypothèse de régime pseudopermanent, c'est à dire à considérer que le flux de l'un à l'autre est proportionnel à la différence des pressions de l'un à l'autre. On a la relation suivante: Tp Fpf = f *(p- pf) (15) avec Fpf: flux puits- fracture Tpf: transmissivité puits-fracture u: viscosité pp: pression du puits pf: pression de la maille de fracture associée Dans l'hypothèse de régime pseudopermanent, la valeur de la transmissivité Tp, est constante au cours de la simulation. Dans la présente méthode, cette valeur est calculée pour chacun des couples maille de fracture - segment puits. La loi de Darcy utilisée pour l'établissement de l'équation différentielle qui régie l'écoulement des fluides dans le milieu poreux, suppose que le régime d'écoulement est laminaire. Dans le cas des gaz et au voisinage du puits, la vitesse du fluide augmente fortement de telle sorte que l'on quitte le régime laminaire (régime turbulent). Dans ce cas précis, la loi de Darcy n'est pas applicable. Pour contourner ce problème, il a été montré 25 qu'en introduisant un terme dit d'écart à la loi de Darcy (dépendant du débit imposé) à la formulation établie pour les régimes laminaire, l'écoulement du fluide est correctement représenté. Ce terme d'écart à la loi de Darcy prend la forme suivante: S'= S0 + D.Q (16) S' a la dimension d'une perte de charge additionnelle. So représente les 5 modifications des propriétés hydrauliques autour du puits (endommagement lors du forage, injection d'acides dans le puits) et le coefficient D, qui est une donnée caractéristique du régime turbulent de l'écoulement. Dans la présente approche, nous introduisons ces coefficients So et D dans l'expression de la transmissivité fissure-puits utilisée pour la représentation du flux entre la 10 fissure et le puits connectés. Dans le cas o l'intersection entre une fissure et le puits est sensiblement circulaire, on définit la transmissivité comme suit: 27rCf (17) 0ln( ) + SO + D(q. + ql avec: r0 = 0. 14 (1f)2 + (h)2, pour les milieux isotropes Cf = kf ef, kf étant la perméabilité intrinsèque de la fracture et ef son ouverture 1f: la longueur de fracture h: l'épaisseur de la couche qn: le débit du segment de puits à l'échelon de débit n. dans le cas o l'intersection fissure-puits est quelconque (elliptique), le calcul de la 20 transmissivité est déduit de la fonction de Peaceman généralisée bien connue des gens de métier. b) Transmissivité matrice-puits Le calcul des échanges entre le puits et une maille matrice associée repose aussi sur l'hypothèse de régime pseudopermanent. Pour exprimer ces échanges, nous combinons les formulations relatives à un écoulement radial autour du puits et un écoulement linéaire vers le puits (Fig. 8). La transmissivité se définit alors comme suit: Tmw = 02lrkml (18) ln(2h)+ 2X(lm -h)+S +Dq (tj) dw, h o: li: la longueur du segment puits lm: la distance du puits à laquelle la pression est la pression moyenne du bloc 10 matriciel * Simulation dynamique Lors de la simulation dynamique, la période de temps simulée est décomposée en pas de temps dont la durée est comprise pour un essai de puits entre lx10-3 s (juste après l'ouverture ou la fermeture du puits) et quelques heures. Le passage de l'instant n (auquel on connaît toutes les valeurs depressions) et l'instant n+l, est réalisé en résolvant l'équation qui suit dans toutes les mailles et blocs du réservoir. Cette équation correspond à l'équation 8 discrétisée. Contrairement au cas o le fluide est faiblement compressible, l'équation que nous présentons est une équation différentielle non linéaire dans la mesure o la compressibilité et la viscosité du fluide 20 varient en fonction de la pression. L'équation différentielle discrétisée s'écrit: nv+l -n (TAp n D..(CT) . i j + .(Jn+l - n+y = 2 Zn Qi (19) t n+1 t"n [ n J *+ avec i: numéro de la maille ou du bloc j: numéros des voisins de i tn: pseudopression à l'instant n pn: pression à l'instant n tn temps à l'instant n (CT),r': la compressibilité totale du fluide pour une pression à l'instant n+1 IL' :la viscosité du fluide pour une pression à l'instant n+ 1 à l'interface entre les mailles i et j zn: le facteur de compressibilité du fluide pour une pression à l'instant n /n: la viscosité du fluide pour une pression à l'instant n Qi: débit entrant dans i en fond de puits Tij transmissivité de liaison entre i et j Pour résoudre cette équation non linéaire, nous utilisons une méthode classique de 15 Newton qui consiste à linéariser l'équation (19) en exprimant: n+l,k+l = n+l,k ji Yi + 8yi (20) La technique consiste alors, à chaque itération k, à résoudre le système linéaire décrit ci-après de telle sorte que la valeur de l'inconnue 6Èi devienne négligeable. Après _.l n+1,k+1 convergence avec k+1 itérations, la variable Y Y+1+ k+i est déterminée et par interpolation linéaire sur la table PVT, toutes les pressions de mailles, la compressibilité et la viscosité du fluide seront réactualisées. Constitution du système à résoudre Calcul de la contribution des mailles de fracture et matrice terme d'accumulation: A (CT)VfnOIk + ( dCT, tn+l _ tn I n+Ik (VJi -'lYi liaison entre mailles de fractures: second membre du système: p Qi tz VI' Yj _____ ( (N +ik yn+,k) -/'% + hIk ó,r Résolution du système linéaire Résolution par l'algorithme itératif du gradient conjugué (CGS) bien connu des gens de l'art avec mise à jour des pressions des mailles de fractures et des blocs matriciels Gestion du pas de temps Dans une simulation d'essai de puits, les pas de temps sont en général très courts 15 après l'ouverture ou la fermeture du puits, parce que c'est à ces moments que les variations de pressions sont fortes. Plus tard, les pas de temps peuvent être plus longs, lorsque les variations de pressions sont plus faibles. La gestion du pas de temps consiste à relier la durée du pas de temps n+1 à la variation de pression maximum observée lors du pas de temps n. Pour ce faire, l'utilisateur 20 doit fournir les données suivantes: Pas de temps initial: dtini Pas de temps minimal: dtmin Pas de temps maximum: dtmax Borne inférieure de variation de pression: dpmin Borne supérieure de variation de pression: dpmax Facteur d'augmentation du pas de temps: rdt+ (>1) Facteur de réduction du pas de temps: rdt- ( dpmax, * inchangé dans les autres cas. Après chaque changement de condition de débit au puits, le pas de temps est réinitialisé à la valeur dtini. De plus, une optimisation est réalisée en fin de période de simulation: si tf est le temps à simuler pour atteindre la fin de période et dt le pas de temps prévu, le pas de temps chosi est min(tf/2,dt). Données d'entrée de la simulation La géométrie du réseau de fractures ainsi que les attributs des fractures (conductivité, ouverture) sont donnés dans un fichier au format défini dans le brevet FR 2 757 947 précité. Les données pétrophysiques à fournir sont les suivantes (on indique l'unité entre parenthèse) : - * compressibilité du réseau de fractures (Cf en bar"l) * compressibilité de la matrice (cm en barf') * perméabilité de matrice (K en mD) * porosité de matrice (Jm sans dimension) Les données en question sont considérées comme homogènes sur le milieu fissuré considéré. Les données concernant les fluides contenues dans le milieu fissuré concernent le gaz (mobile) et l'eau (immobile) qui est toujours présente en quantité résiduelle dans la matrice: * saturation résiduelle en eau dans la matrice (Sw sans dimension) * compressibilité de l'eau (c, en bail) * une table PVT comprenant le facteur de compressibilité du gaz (Z) et la viscosité du gaz (it en cpo) en fonction de la pression. Les données concernant le puits sont: * géométrie sous la forme d'une série de segments connectés * débits imposés sous la forme d'une courbe donnant le débit imposé en fonction du temps. * le coefficient So d'endommagement de la formation autour du puits, et le coefficient D d'écart à la loi de Darcy. Les valeurs à fournir sont les valeurs déjà définies dans le chapitre sur la gestion du pas de temps: dtini, dtmin, dtmax, dpmin, dpmax, rdt+ et rdt-. L'organigramme de la figure 9 récapitule les différentes étapes de la modélisation réalisée. REVENDICATIONS 1) Méthode pour modéliser des flux de fluides compressibles dans un milieu poreux multicouches traversé par un réseau de fractures de géométrie donnée et par un puits, comprenant les étapes suivantes a) on discrétise chaque couche fracturée par un maillage comprenant des mailles de fracture que l'on centre sur des noeuds soit aux intersections de fractures soit aux extrémités des fractures, chaque noeud étant associé à un bloc de matrice regroupant tous les points plus proches de lui que de noeuds voisins; b) on calcule les flux locaux entre chaque hmaille de fracture et le volume de matrice 10 associé dans un régime de flux pseudo permanent, la valeur de la transm-iissivité étant obtenue en considérant une variation linéaire de la pression en fonction de la distance entre tout point du bloc et la maille de fracture; c) on détermine les flux directs entre les mailles de fracture; d) on détermine les flux directs entre les volumes de matrice au travers des arêtes 15 communes des mailles; e) on détermine les flux directs entre d'une part les mailles de fracture et les mailles de matrice, et d'autre part le volume du puits, à partir de valeurs de transmissivité adaptées à des écoulements turbulents; et f) on simule la réponse du puits pour des variations de débit imposées au moyen 20 d'une équation de diffusion modifiée pour tenir cokpte de la compressibilité du fluide, qui relie les interactions entre les variations de la pression et du débit observables dans le puits. 2) Méthode selon la revendication 1, dans laquelle on modifie l'équation de diffusion par introduction d'un terme dépendant de la pression, de la viscosité et du facteur de compressibilité du fluide compressible. 3) Méthode selon la revendication i ou 2, dans laquelle on modifie l'équation de diffusion par introduction d'une pseudopression tenant compte de la variation de viscosité et de compressibilité des fluides en fonction de la pression. 23 2850773 - -23 4) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on détermine les flux directs en introduisant un terme d'écart à la loi de Darcy proportionnel au débit de fluides imposé, traduisant les pertes de charge localisées aux abords du puits liés à la vitesse élevée des dits fluides. 5) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on détermine le bloc de matrice associé à chaque maille de fracture en discrétisant chaque couche fracturée en un ensemble de pixels et en calculant la distance de chaque pixel à la maille de fracture la plus proche, pour déterminer l'emplacement des arêtes entre les mailles et l'on utilise les distances calculées pour en déduire la valeur des transmissivités 10 entre chaque maille de fracture et le volume de matrice associé d'une part, et entre les arêtes communes des mailles. 6) Méthode selon la revendication 5, caractérisée en ce que l'on repère la position des arêtes entre les mailles en déplaçant deux fenêtres allongées dans l'image formée par les pixels, orientées suivant deux directions perpendiculaires.