FR 2858444 A1 20050204 FR 0309360 A 20030729 La présente invention concerne une méthode pour simuler les transferts compositionnels et/ou polyphasiques entre la matrice poreuse et les fractures d'un milieu poreux multicouches fracturé. La méthode selon l'invention trouve des applications par exemple pour simuler la production d'un gisement d'hydrocarbures fracturé dans le cas où des transferts polyphasiques se produisent entre des blocs matriciels initialement saturés en huile et les fractures limitrophes balayées par du gaz frais injecté, dans le cadre d'opérations de récupération assistée. Etat de la technique L'état de la technique auquel on va faire ci-après référence, se trouve défini et expliqué dans les publications suivantes : - Chen, W. H., M.L. Wasserman and R. E. Fitzmorris 1987. A Thermal Simulator for Naturally Fractured Reservoirs. Paper SPE 16008 presented at the 9th SPE Symposium on Réservoir Simulation held in San Antonio, Tx, Feb. 1-4,1987 ; - Gilman, J. R. 1986. An Efficient Finite-Difference Method for Simulating Phase Segregation in the Matrix Blocks in Double-Porosity Reservoirs. SPE Reservoir Engineering, July 1986. Pages 403-413 ; - Kazemi, H., Merrill, L. S., Porterfield, K. L. and Zeman, P. R. 1976. Numerical Simulation of Water-Oil Flow in Naturally Fractured Reservoirs. SPE Journal, Dec. 1976,317; - Pruess, K. and T. N. Narasimhan 1985. A Practical Method for Modelling Fluid and Heat Flow in Fractured Porous Media. Society of Petroleum Engineers Journal, Feb. 1985. Pages 14-26 ; - Quintard, M. and Whitaker, S. 1996. Transport in Chemically and Mechanically Heterogeneous Porous Media. Advances in Water Resources, 19(1), 29-60; - Sabathier, J.C., B.J. Bourbiaux, M. C. Cacas and S. Sarda 1998. A New Approach of Fractured Reservoirs. Paper SPE 39825 presented at the SPE International Petroleum Conference and Exhibition of Mexico held in Villahermosa, Mexico, 3-5 March 1998 ; Saïdi, A.M. 1983. Simulation of Naturally Fractured Reservoirs. Paper SPE 12270 presented at the 7th SPE Symposium on Réservoir Simulation held in San Francisco, CA, Nov. 15-18,1983; ou - Warren, J. E. and P. J. Root 1963. The Behavior of Naturally Fractured Reservoirs. Society of Petroleum Engineers Journal, Sept. 1963. Pages 245-255. En outre, différentes méthodes permettant de simuler des flux dans des milieux fracturés ont fait l'objet par le demandeur, des demandes de brevet suivants : FR 2 809 894, FR 02/03 436 et FR 03/01 090. Des estimations fiables concernant un réservoir fracturé en termes de productivité et de récupération nécessitent une approche multi-disciplinaire afin de minimiser les incertitudes liées à la modélisation du réseau de fractures et à la simulation de l'écoulement polyphasique. En particulier les trois étapes suivantes jouent un rôle déterminant : 1. Construction d'un modèle représentatif de réseaux de fractures naturelles à partir des informations de terrain disponibles, 2. Conversion de ce modèle géologique en un modèle double porosité équivalent, 3. Obtention d'une bonne estimation des mécanismes de production physiques en jeu et reproduction de ces derniers au moyen du simulateur double porosité polyphasique. La méthode qui va être décrite, porte sur la troisième étape, et plus précisément sur un type de formulation adaptée pour simuler les transferts matrice-fracture complexes, c'est-à-dire polyphasiques et/ou compositionnels et parfois thermiques en jeu dans la plupart des procédés de récupération par injection de gaz. En effet, les prédictions d'un simulateur double porosité sont extrêmement sensibles à la formulation des transferts matrice-fracture dans la mesure où la majeure partie de l'huile en place se situe dans le milieu matriciel. Le modèle conceptuel double porosité représente le réservoir fracturé sous la forme d'un ensemble de blocs matriciels parallélépipédiques limités par un ensemble de fractures orthogonales uniformes (Fig.l). Deux grilles superposées, représentant les deux milieux fracture et matrice, sont utilisées pour le calcul des écoulements. Les écoulements de fracture sont calculés entre mailles de la grille fracture, les transferts matrice-fracture sont calculés pour chaque couple de mailles matrice et fracture superposées, et les écoulements au sein de la matrice sont également pris en compte dans l'option double porosité de ces simulateurs. Warren et Root ont proposé l'expression suivante pour le flux matrice-fracture (par unité de volume de matrice) en se référant à un transfert monophasique en régime quasi-stationnaire contrôlé par la diffusivité de pression : où pm et pf sont les pressions moyennes de matrice et de fracture. Cette équation comporte une constante de proportionnalité # appelée facteur de forme (dimension : 1/L2), qui reflète la géométrie du bloc matriciel non discrétisé de dimensions (Lx,Ly,Lz). Kazemi a proposé l'expression suivante pour # : Le transfert de masse (flux massique par unité de volume de matrice) dû à la diffusion moléculaire est exprimé d'une manière similaire : où D est le coefficient de diffusion, ()) la porosité, # le facteur de tortuosité, p la masse volumique du fluide et C la concentration des composants au sein du fluide. L'extension de ces expressions à des écoulements polyphasiques donne lieu à des difficultés liées à la représentation d'un phénomène transitoire qui peut difficilement être reproduit sans discrétiser le bloc matriciel. C'est la raison pour laquelle une discrétisation des blocs matriciels a été proposée (cf. Saïdi, 1983) et introduite dans certains simulateurs (cf. Pruess et al., (1985) ou Gilman, (1986) ou Chen et al., (1987). Cette approche est satisfaisante mais elle comporte des besoins en matière de calcul qui ne peuvent pas être satisfaits pour des modèles de réservoirs complexes ou de grande taille. C'est pourquoi, dans le cadre de la présente invention, nous conservons la situation pratique dans laquelle les blocs matriciels ne sont pas discrétisés et les transferts matrice-fracture sont représentés par un seul terme source exprimé en fonction des propriétés matrice/fracture et de variables (pressions, saturations, compositions). Les difficultés s'accroissent encore dans le cas de transferts polyphasiques de fluides en déséquilibre thermodynamique car le transfert de masse des composants s'effectue à l'interface entre les phases et s'ajoute aux transferts convectifs et diffusifs qui ont lieu au sein de chaque phase. Ainsi, il n'existe à ce jour, à notre connaissance, aucun simulateur double porosité capable de simuler de manière fiable les transferts compositionnels polyphasiques entre un bloc matriciel rempli d'huile et une fracture remplie de gaz en déséquilibre, sauf si un processus convectif (la gravité par exemple) déclenche le transfert du gaz dans la matrice. En effet, dans de telles situations, les composants gazeux doivent d'abord être dissous dans la phase huile de la matrice jusqu'à ce que cette dernière soit entièrement saturée de gaz et qu'une phase gazeuse apparaisse dans le bloc matriciel. Ce n'est qu'à ce stade que le processus de diffusion moléculaire monophasique fracture-matrice peut être initié dans la simulation. C'est l'objectif de la présente invention de présenter une méthode permettant de simuler de façon réaliste ces phénomènes de transfert compositionnel. LA METHODE SELON L'INVENTION La méthode selon l'invention permet de simuler les transferts compositionnels et/ou polyphasiques entre la matrice poreuse d'une part et les fractures d'autre part d'un milieu poreux multicouches fracturé, saturées respectivement par une phase liquide et une phase gazeuse distinctes. Elle comporte la construction d'un modèle géologique représentatif des réseaux de fractures naturelles du milieu et la conversion de ce modèle géologique en un modèle double porosité équivalent, la modélisation de la diffusion de la phase gazeuse entre les fractures et la matrice en considérant qu'un front gaz-liquide présentant un équilibre quasi-statique aux interfaces s'établit à partir d'un instant initial et progresse à une vitesse contrôlée par les transferts de diffusion au sein de chaque phase, lesquels transferts font intervenir la différence entre la concentration au sein de la phase et la concentration d'équilibre entre les deux phases, ainsi qu'une distance d'échange variable exprimée en fonction de la saturation de matrice. PRESENTATION 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, en se référant aux dessins annexés où : - la figure 1 montre une représentation classique d'un milieu fracturé dans le modèle à double milieu selon Warren et Root, Lx, Ly, Lz représentant les dimensions du bloc matriciel équivalent et kx, ky et kz, les perméabilités de fracture équivalentes ; - la figure 2 montre la distribution des compositions d'un constituant dans deux phases avec équilibre thermodynamique à l'interface (un gradient de composition existe entre l'interface et le milieu de la phase) ; - la figure 3 montre un exemple d'évolution de la quantité d'huile totale au cours du temps d'un mélange binaire; - la figure 4 montre un exemple d'évolution de la masse totale de Cl au cours du temps d'un mélange binaire; - la figure 5 montre l'évolution de la masse totale de C5 au cours du temps d'un mélange binaire ; - la figure 6 : un exemple d'évolution de la quantité d'huile totale au cours du temps d'un mélange ternaire; - la figure 7 montre un exemple d'évolution de la masse totale de Cl au cours du temps d'un mélange ternaire; - la figure 8 montre un exemple d'évolution de la masse totale de C5 au cours du temps d'un mélange ternaire; - la figure 9 montre un exemple d'évolution de la masse totale de C16 au cours du temps d'un mélange ternaire; - les figures 10 et 11 montrent respectivement un exemple de profil de saturation en gaz pour les mélanges binaire et ternaire obtenu au cours du temps (Fig.10) et et un exemple de profil en fonction de la profondeur (Fig.l 1) ;. et la figure 12 montre des profils de saturation à différentes profondeurs dans un échantillon de roche. DESCRIPTION DETAILLEE On va décrire ci-après la formulation de la méthode, ainsi que son extension à des conditions d'échange 2D et 3D, et sa validation par rapport à des simulations de référence sur modèles finement maillés. Formulation de transferts d'huile et de gaz hors d'équilibre entre une fracture et un bloc matriciel non maillé Principes sous-jacents Les transferts entre une fracture balayée en continu par du gaz d'injection frais et un bloc matriciel initialement saturé d'huile reposent sur la modélisation suivante des phénomènes physiques en jeu : # Un transfert de composants s'effectue à travers l'interface entre l'huile et le gaz de manière à établir un équilibre local entre les deux phases : cet équilibre local est obtenu instantanément par rapport au temps requis pour l'établissement d'un équilibre global entre les fluides de fracture et de matrice, grâce à des phénomènes de diffusion s'opérant à l'échelle du bloc, # Une diffusion moléculaire s'effectue au sein de chaque phase du fait de la différence de composition entre l'interface entre phases et la phase prise dans son ensemble, # En raison de changements volumétriques au sein de chacune des phases, l'interface vapeur-liquide progresse (gonflement) ou recule (vaporisation) au sein du milieu matriciel. Le modèle simplifié de ces phénomènes consiste à considérer qu'un front vapeurliquide présentant un équilibre quasi-statique à l'interface s'établit à partir de l'instant initial et progresse à une vitesse contrôlée par les transferts de diffusion au sein de chaque phase. Un tel modèle permet de simuler les flux d'échange matrice-fracture lorsque la fracture est entièrement saturée de gaz et la matrice est entièrement saturée d'huile, à condition de mettre en #uvre la formulation décrite en détail ci-après. Formulation d'échange 1D Pour un bloc matriciel saturé de liquide et de dimensions latérales Lx, Ly, Lz, en contact avec une fracture saturée de gaz sur 2 faces opposées (échange 1D dans la direction X), l'expression (équation 3) de la formulation classique du flux de diffusion pour chaque phase (vapeur ou liquide) peut être réécrite comme suit : où V = Lx Ly Lz (volume du bloc matriciel), A = 2LyLz , est la section d'échange pour les 2 faces opposées, et d#, est la distance sur laquelle s'effectue l'échange de diffusion (dimension d du demi-bloc, c'est-à-dire Lx/2 dans le cas présent, multipliée par la tortuosité # prise égale à la tortuosité de la matrice #m). Si la matrice et la fracture sont saturées avec des phases différentes (huile et gaz par exemple), la formulation précédente écrite pour la phase vapeur ne prédit aucun flux car la concentration Cm est nulle. De même, aucun flux de diffusion n'est prédit dans la phase liquide car Cf est égal à zéro. De plus, la distance d'échange est fixée (Lx/2) et ne prend pas en compte la dynamique de fluide décrite dans la section précédente. Les deux améliorations majeures apportées à cette formulation sont les suivantes : A. Pour calculer un flux de diffusion de gaz entre fracture et matrice, une concentration fictive Cm est adoptée tant qu'il y a absence d'une phase vapeur dans le milieu matriciel. Cette concentration fictive est celle à l'équilibre avec la composition de la phase liquide de la matrice. Elle est obtenue par calcul d'équilibre entre la phase liquide matrice et la phase vapeur fracture. Dès qu'une phase vapeur est présente dans le milieu matriciel, la composition de la phase vapeur effectivement présente dans la matrice est ré-utilisée dans l'équation de flux de diffusion. B. La distance d'échange (d#) n'est plus constante mais variable de manière à prendre en compte la progression de la phase vapeur dans le bloc matriciel : d# est remplacé par la somme d'une distance fixée représentant le trajet de diffusion au sein de la fracture (c'est-à-dire la moitié de l'épaisseur de la fracture, une très courte distance mais qui n'est pas nulle à l'instant initial) plus le trajet de diffusion variable au sein du bloc matriciel. En admettant que la phase vapeur se développe sous la forme d'un front qui progresse vers le c#ur du bloc, le trajet de diffusion variable peut être exprimé en fonction de la saturation en phase vapeur Sg au sein du bloc matriciel. On obtient finalement, pour un transfert 1D à partir des 2 faces opposées d'un bloc de longueur Lx : d#= ef+ [gamma](Lx/2).#m.Sg ......................... (5) ef est la moitié de l'épaisseur de la fracture (nous admettons une tortuosité égale à 1 dans la fracture, mais une valeur supérieure à 1 pourrait également être introduite). y(Lx/2) # Sg est le trajet effectif de diffusion au sein du milieu matriciel, qui représente la longueur de diffusion moyenne au sein de la partie saturée de vapeur du milieu matriciel. En supposantle déplacement d'un front, la longueur saturée de vapeur est (Lx/2).Sg à partir de l'une quelconque des 2 faces ouvertes du bloc, et (Lx/2).#.Sg si nous prenons en compte la tortuosité des pores de la matrice. Si l'on considère les résultats obtenus avec les techniques de prise de moyenne volumique ( cf. Quintard & Whitaker), la longueur moyenne de diffusion du gaz au sein de cette portion saturée de vapeur du bloc est égale à [gamma](Lx/2).#.Sg, y étant un facteur constant égal à 3 dans le cas présent d'échange 1D. Conformément à ces modifications, la porosité (}) de l'équation (4) est égale à la porosité intrinsèque de la fracture (généralement égale à 1) tant qu'il y a absence d'une phase vapeur dans le bloc matriciel, puis à la porosité de la matrice. Extension de la formulation aux échanges 3D Nous considérons maintenant que toutes les faces du bloc matriciel sont limitées par des fractures. Là encore, nous admettons un bloc matriciel saturé de liquide et des fractures remplies de gaz. En ce qui concerne la concentration fictive de la phase absente, on applique la même procédure que celle décrite ci-avant (≈1.2, A). L'extension décrite ci-dessous concerne le calcul de la distance et de la section d'échange. Nous réécrivons l'expression classique généralisée du flux de diffusion : Si Tx,Ty et/ou #z diffèrent les uns des autres, on peut définir les longueurs suivantes pour obtenir une forme isotrope de l'équation (6) : alors : Comme précédemment, nous considérons un front de vaporisation qui progresse dans le bloc. C'est pourquoi, contrairement à la formulation (6), la longueur et la section d'échange ne sont plus constantes. La longueur d'échange est exprimée sous la forme d=ef+[gamma]#i.X, avec : - X : une fonction de la saturation du bloc matriciel, - i : direction suivant laquelle on a la distance minimale d'échange dans le problème rendu isotrope (i=i[Min(L'x,L'y,L'z)=L'j]), - d : distance d'échange variant d'une valeur initiale égale à l'épaisseur de la fracture, ef, à une valeur finale égale à e+yr.Min(L' x,L' y,L' z)= e+yr. L'i. De plus, pour cette géométrie 3D, en accord avec l'hypothèse de l'existence d'un front, la section d'échange au sein du milieu poreux à un instant t quelconque n'est plus constante mais égale à une section moyenne d'écoulement comprise entre la section des 6 faces externes du bloc (2(L'yL'z + L'XL'Z + L'xL'y)) et la section d'écoulement au front, c'est-à-dire : La section moyenne d'écoulement (A) peut être calculée à partir de la somme suivante entre 0 et X (moyenne harmonique des sections locales dans la zone vaporisée) : La position du front, X, est reliée à la saturation normalisée en gaz du bloc matriciel, Sg , comme suit : où 0 En réalité, des simulations de référence réalisées sur des modèles de bloc finement maillés ont révélé que la prise en compte d'une section d'échange variable et non plus constante n'était pas nécessaire et même moins représentative de la physique réelle du transfert pour les deux raisons suivantes : (a) tout d'abord, en raison de la relation de dépendance polynomiale entre X and Sg (Eq. 10), la section d'échange moyenne demeure très proche de la section des faces latérales tant que la saturation en gaz du bloc matriciel n'excède pas 40% ; (b) l'hypothèse d'un front de vaporisation raide progressant à l'intérieur du bloc est physiquement représentative au cours de la période initiale de transfert mais devient moins adaptée au fur et à mesure que le temps s'écoule : en effet, l'étalement du front de vaporisation avec le temps compense l'effet adverse d'une réduction de la section d'échange. Ainsi, à 3 dimensions, l'usage d'une concentration fictive de phase absente et l'adoption d'une distance d'échange variable demeurent les principales innovations permettant de prévoir les transferts gaz-huile par diffusion au moyen d'un simulateur à double porosité. Exemples de validation Les figures 3 à 10 illustrent les comparaisons entre la formulation de référence simple milieu, la modélisation double milieu standard et les modélisations double milieu améliorées avec section constante et variable sur un cas de mélange binaire (Cldans la fissure et C5 dans la matrice) et ternaire (Cl dans la fissure et C1-C5-C16 dans la matrice). On constate que grâce à la formulation améliorée, le double milieu est en mesure de reproduire le comportement du simple milieu contrairement à la formulation standard qui est totalement inadaptée. 1) Méthode pour simuler les transferts compositionnels et/ou polyphasiques entre la matrice poreuse d'une part et les fractures d'autre part d'un milieu poreux multicouches, saturées respectivement par une phase liquide et une phase gazeuse distinctes, comportant la construction d'un modèle géologique représentatif des réseaux de fractures naturelles du milieu et la conversion de ce modèle géologique en un modèle double porosité équivalent, caractérisé en ce que l'on modélise la diffusion de la phase gazeuse entre les fractures et la matrice en considérant qu'un front gaz-liquide présentant un équilibre quasi-statique aux interfaces s'établit à partir d'un instant initial et progresse à une vitesse contrôlée par les transferts de diffusion au sein de chaque phase, ces derniers faisant intervenir la différence entre la concentration au sein de la phase et la concentration d'équilibre entre les deux phases, ainsi qu'une distance d'échange variable exprimée en fonction de la saturation de matrice.