FR 2837592 A1 20030926 FR 0203436 A 20020320 zo de visualisation (2), font partie d'un ensemble spécifique. S La présente invention concerne une méthode pour modéliser des flux de fluides dan s un milieu poreux multi couches traversé par un réseau de fractures de géométrie donnée, inégalement réparties dans le rnilieu, et dont certaines communiquent les unes avec les autres. La méthode selon l'invention trouve des applications notamment dans le domaine de la production pétrolière. 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 un zone souterraine poreuse perméable comportant deux milieux contrastés, un milieu matriciel contenant la plus grande part de l'huile en place et présentant une faible perméabilité, et, le traversant en partie ou en totalité, des réscaux de fractures beaucoup plus conductrices. 1S Etat de la technique Lors d'un essai de puits, les conditions de débit imposées au puits amènent l'huile (le pétrole) contenue dans le réservoir à s'écouler vers le puits. I1 s'agit d'un écoulement monophasique (la seule phase mobile est l'huile) et compressible (la roche du réservoir et le fluide huile sont compressibles). Pour tout volume élémentaire du réservoir, la pres si on de l 'hui le contenue dans ce volume oboit alors à l'équation suivante (on ne tient pas compte de la gravité): CT À at = div( VP) + Q (1) avec: 4: volume poreux CT: compressibilité totale (fluide + roche) K: perméabilité de la roche : viscosité du fluide Q: débit entrant P: inconnue pression 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 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 connat des outils de modélisation d'écoulements monophasiques qui sont appliqués non au milieu géologique rcel dans toute sa complexité mais à une représentation homogéncisoe, suivant le modèle de réservoir dit de double milieu décrit par exemple par 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 fluides surviennent localement entre les fissures et les blocs matriciels. Cette représentation qui ne rend pas compte de la complexité du réscau de fractures dans un réservoir, se révèle efficace néanmoins mais à l'échelon d'une maille de réservoir ayant typiquement pour dimensions 100m x 100m. I1 est bien préférable cependant pour les ingénieurs de réservoir, de disposer d'un simulateur d'écoulement basé sur un modèle géologique " récl " du milieu et non sur un modèle homogène équivalent, de façon à: À 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 donnces rcelles d'essais de puits); À tester de manière fiable la sensibilité du comportement hydraulique du milieu aux incertitudes sur l'image géologique du milieu fissuré. Une méthode de modélisation connue consiste à mailler finement le réseau de fractures et la matrice en ne fai s ant aucune approximati on concern ant l es éch anges de fluides entre les deux milieux. Elle est cependant difficile à mettre en _uvre car la géométrie souvent complexe des espaces entre les fractures rend difficile leur maillage et qu'en tout état de cause, le nombre de mailles à traiter finalement est souvent énorme. La complexité augmente encore s'il s'agit d'un maillage en 3D. Des techniques de modélisation de milieux poreux fissurés sont décrites dans les brevets FR-A-2.757.947 et FR-A-2.757.957 du demandeur. La première est relative à la détermination de la perméabilité de fracture équivalente d'un réseau de fractures dans un milieu multicouches souterrain à partir d'une représentation connue de ce réseau, permettant de relier de manière systématique des modèles de caractérisation de réservoir fissuré à des simulateurs double porosité en vue de réaliser une modélisation plus réaliste d'une structure géologique souterraine fissurée. La deuxième est relative à la modélisation simplifiée d'un milieu géologique hétérogène poreux (tel qu'un réservoir traversé par un réscau irrégulier de fractures par exemple) sous la forme d'un milieu transposé ou équivalent de manière à ce que le milieu transposé soit équivalent au milieu d'origine, relativement à un type déterminé de fonction de transfert physique (connu pour le milieu transposé). Par la demande de brevet FR 00/06 874 du demandeur, on conna^t aussi une méthode de simulation des écoulements de fluides dans un milieu géologique poreux fissuré traversé par un réscau de d' objets conducteurs de fluides de géométrie définie mais non homogéncisables à l'échelle de chaque maille du modèle (des grandes fractures, des failles sub-sismiques par exemple, des couches sédimentaires très perméables, etc.), o l' on détermine les échanges intervenant entre le milieu matrice et le milieu fracture, et une modélisation des transmissivités des différentes mailles traversées par chaque objet conducteur, de façon que la transmissivité résultante corresponde à la transmissivité directe le long de chacun de ces objets. Dans le cas o les objets conducteurs sont des couches sédimentaires très perméables, on attribue à la transmissivité entre mailles traversées par chaque couche très perméable une valeur dépendant des dimensions des mailles et de l'aire de contact commune entre couches à la jonction des mailles adjacentes. Dans le cas o les objets conducteurs sont des fractures, on attribue à la transmissivité entre mailles traversées par chaque fracture, une transmissivité dépendant des dimensions des mailles et de l'aire commune de la fracture à la jonction des mailles adjacentes. Par le brevet FR 2 787 219 du demandeur, on conna^t également une méthode pour modéliser les flux de fluides dans un milieu poreux multicouches fissuré en tenant compte de la géométrie réelle du réscau de fractures et des échanges locaux entre la matrice poreuse et les fractures en chaque n_ud du réscau. On discrétise le milieu fissuré par un maillage, on centre les mailles de fracture sur des n_uds aux différentes intersections des fractures, chaque n_ud étant associé à un volume ou bloc de matrice, et l'on détermine les flux entre chaque maille de fracture et le volume de matrice associé dans un régime de flux semi- permanent. La méthode permet de calculer le volume de chaque bloc ainsi que du coefficient de proportionnalité (appelé transmissivité) qui relie de débit matrice-fracture à la différence de pression matrice-fracture. Cette méthode s'applique bien aux cas o un réseau dense de fractures traverse l'ensemble du volume matrice considéré. Dans ces cas, les écoulements à grande échelle ont lieu dans le réseau de fractures et la matrice n'intervient que dans les échanges entre matrice et fractures locales. La méthode selon l'invention La méthode selon l'invention a pour objet de généraliser la technique décrite dans le brevet FR 2 787 219 précité, dans le cas o les écoulements à grande échelle traversent non seulement des zones avec des réseaux relativement denses de fractures mais également des zones (certaines couches par exemple) que sont peu fracturées. La méthode selon l'invention permet de modéliser des flux de fluides dans un milieu poreux multi couches traversé par un réseau de fractures de géométrie donnce, dont certaines communiquent les unes avec les autres. Elle comprend les étapes suivantes: a) on discrétise chaque couche fracturce par un maillage comprenant des mailles de fracture que l' on centre sur des n_uds soit aux intersections de fractures soit aux extrémités des fractures, chaque n_ud étant associé à un bloc de matrice regroupant tous les points plus proches de lui que de n_uds voisins; s 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 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 fracture; d) on détermine les flux directs entre les volumes de matrice au travers des arêtes communes des mailles; et e) on simule les interactions entre les variations de la pression et du débit observables dans au moins un puits au travers du milieu. On peut déterminer le bloc de matrice associé à chaque maille de fracture en discrétisant chaque couche fracturce en un ensemble de pixels et en calculant la distance de chaque pixcl à la maille de fracture la plus proche, pour déterminer l' emplacement des arêtes entre les mailles. 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 fenêtres allongées dans l'image formée par les pixels, orientées suivant deux directions perpendiculaires. En tenant compte aussi des transmissivités entre blocs matriciels associés à des n_uds 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. Présentation sommaire des figures D'autres caractéristiques et avantages de la méthode selon l'invention, apparâîtront à la lecture de la description ci-après, en se référant aux dessins annexés o: - la Fig. 1 montre un exemple de maille de fracture en 2D; - la Fig.2 montre un exemple de bloc matriciel obtenu en 2D; - la Fig.3 montre un exemple de deux blocs matriciels communiquant par l'intermédiaire d'un élément de fracture; - la Fig.4 montre un exemple de deux blocs matriciels traversés chacun par un élément de fracture et communiquant par le biais d'échanges au travers de leurs arétes communes; - la Fig.S montre des exemples de fenétres orientées suivant deux axes perpendiculaires, S permettant le repérage des arêtes entre les blocs matriciels; - la Fig.6 montre un exemple de variation de déhit imposé dans un essai de puits; et la Fig.7 montre un exemple de couche avec une répartition très inégale de fractures o la méthode s'applique avantageusement. Description détaillée - Dans le cas d'un milieu fissré, si l'on veut tenir compte de toute la complexité du réscau de fractures, il est nocessaire de réaliser un maillage explicite de ce réseau. Par 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. Enfin, dans certains cas, le réscau de fractures n'est pas présent partout dans le volume poreux et l'on doit tenir compte des écoulements à grande échelle à la fois dans les fractures et dans le milieu matrice. Partant de ces constatations, l'approche présentée consiste à: À mailler explicitement le réseau de fractures; À associer à chaque maille fracture un bloc matriciel unique; À traiter les écoulements entre chaque maille fracture et son bloc associé en pseudo permanent, c'est-à-dire que le flux de l'un à l'autre est proportionnel à la différence de pression entre les deux; et À tenir compte des écoulements entre blocs matriciels. Les inconnues du problème sont les pressions des mailles de fracture et les pressions des blocs matriciels associés. Puisqu'il y a autant de blocs matriciels que de mailles de fracture, le nombre d'inconnues est égal à 2 fois ce nombre. Discrétisation du réseau de fractures En 2D, les n_uds 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 n_uds 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 n_uds de calcul ainsi définis constituent le centre des mailles fracture. De plus, la simulation d'écoulements compressibles nécessite de connâitre le volume des mailles ( dans l'équation 2-1). On définit donc des limites de mailles positionnées au milieu des segments reliant les n_uds de calcul. Le schéma de la Fig.1 montre un exemple de maille fracture en 2D. En terme de nombre de n_uds de calcul, le maillage du réscau de fractures est ainsi rendu optimal. Le calcul des liaisons entre mailles fracture (transmissivités), utilisées pour le calcul des flux entre ces mailles fracture, est identique à celui présenté dans la méthode 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 fracture définie comme indiqué dans le paragraphe précédent. Lors de la simulation dynamique, les échanges matrice-fracture seront calculés en pseudo permanent entre chaque maille 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 donnce, 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 fracture est l'ensemble des points plus proches de cette maille 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 fracture. Pour déterminer la géométrie des blocs dans une couche donnce, il s'agit donc de résoudre un problème bidimensionnel qui consiste à trouver, pour chaque maille 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 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 mai l le fracture auquel il appartient, le même al gori thme permet finalement de déterminer, pour chaque pixel: la distance de ce pixcl à la maille fracture la plus proche; Ie numéro de la maille fracture la plus proche. L'ensemble des pixcls ayant le méme numéro de maille 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. Pour chaque bloc matriciel, l'algorithme de traitement d'image donne aussi la distance de chaque pixel du bloc à la maille fracture associce. On se sert de cette information pour calculer la transmissivité entre la maille fracture et le bloc matriciel. L'hypothèse de régime pseudo permanent entre une maille fracture et un bloc matriciel revient à considérer que le flux de 1'un à l'autre est proportionnel à la différence des pressions de l'un et l'autre. On a la relation suivante: TFmf = mf (Pm - Pf) (2) avec: Fmf: flux matricefracture Tmf: transrnissivité matrice-fracture p: viscosité Pm pression du bloc matriciel Pf: pression de la maille fracture associce Dans l'hypothèse du régime pseudo permanent, la valeur de la transmissivité T est constante au cours de la simulation. Dans la présente méthode, cette valeur est calculée pour chacun des couples maille fracture - bloc matriciel. Pour ce calcul, on fait l'hypothèse que dans le bloc matriciel, la pression varie S linéairement en fonction de la distance du point considéré à la maille fracture associée au bloc. On définit la transmissivité comme suit: Tmf = D (3) avec 1: longueur des fractures dans la maille 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 connâî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 donnces par l'algoritUme de traitement d'image 1S sur les distances des pixcls 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: D= S Jd(s) ds (4) 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 fracture associée. En terme de pixcls, on a donc simplement: D = N À, dn o N est le nombre de pixels du bloc matriciel et dn la distance du pixel n à la maille 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 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 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 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 pixcl de l'image, sur la distance de ce pixel à la fracture la plus proche et sur le bloc matriciel auquel appartient ce pixel. L'algorithwe de traitement d'image consiste à repérer les interfaces élémentaires entre blocs, c'est à dire les arêtes séparant deux pixcls 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 suivante: M p2 _ (din-dfn)2 T. = k À H, 2. di n (6) o Tn 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 din et dfn 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. I1 y a une fenêtre verticale pour repérer les interfaces horizontales et une fenêtre horizontale pour repérer les interfaces verticales. Dans le cas de figure 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 (7) AB o I est la longueur de l'arête et LAB la longueur de la fracture reliant les centres A et B. La 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 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éscau. 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 entrces 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 es sai de puits c l as sique, on imp o se un débi t 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ématisoe sur la Fig. 6. Simulation dynamique Lors de la simulation dynamique, la période de temps simulée est décomposée en pas de temps dont la durce est comprise pour un essai de puits entre 1 s (juste après l'ouverture ou la fermeture du puits) et quelques heures. Le passage de l'instant n (auquel on connâît toutes les valeurs de pressions) et I'instant n+l, est réalisé en résolvant l'équation qui suit dans toutes les mailles et blocs du réservoir. Cette équation est l'équation 2-1 discrétisce. pn+1 _ pn T i CT. n+ i +y, J.(pin+1_pn+l) Q (8) avec i: numéro de la maille ou du bloc j: numéros des voisins de i pn: pression à l'instant n tn: temps à l'instant n Qi: débit entrant dans i Tij: transmissivité de liaison entre i et j A part les pressions, les termes de cette équation sont connus. Les volumes poreux des mailles de fracture et des blocs matriciels (qJ) sont connus par le maillage, les transmissivités fracturefracture, matrice-matrice et matrice-fracture (Tij) sont calculées comme décrit plus haut et les débits (Qi) sont nuls partout sauf aux connexions puits réservoir o ils sont imposés. L'algorithme général de résolution d'un pas de temps consiste à constituer un système linéaire correspondant à l'équation 8 écrite pour toute les mailles fracture et mailles matrice i, à résoudre ce système linéaire et à mettre à jour les inconnues Les étapes de cet algorithme sont les suivantes. Constitution du système linéaire (calcul des termes de l'équation 8) Calcul de la contribution des mailles fracture et matrice terme d'accomulation ( i) liaison entre mailles de fractures conditions limites puits (Qi) Résolution du système linéaire (pour trouver les valeurs de pression à l'instant n+1) Résolution par l'algorithme itératif du gradient conjugué (CGS) Mise à jour Mise à jour des pressions des mailles fractures et des blocs matriciels Gestion du pas de temps Dans d'une simulation d'essai de puits, les pas de temps sont en général très courts 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 durce du pas de temps n+1 à la variation de pression maximum observée lors du pas de temps n. Pour ce faire, l'utilisateur 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éscau de fractures ainsi que les attributs des fractures (conductivité, ouverture) sont donnés dans un fichier au format défini dans lé 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) À compressibilité de la matrice (cm en bar) À perméabilité de matrice (K en mD) À porosité de matrice ((Pm sans dimension). Les données en question sont considérées comme homogènes sur le rnilieu fissuré considéré. Les données concernant les fluides contenues dans le milieu fissuré concernent l'huile (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 (cw en bar) À compressibility de l'huile (cO en bard) À viscosité de l'huile (p en cpo) Les données concernant le puits sont: À géométrie sous la forme d'une série de segments connoctés S À débits imposés sous la forme d'une courbe donnant le débit imposé en fonction du temps. 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-. REVENDICATIONS 1) Méthode pour modéliser des flux de fluides dans un rnilieu poreux multicouches traversé par un réseau de fractures de géométrie donnée, dont certaines communiquent les unes avec les autres, caractérisée en ce que: f) on discrétise chaque couche fracturée par un maillage comprenant des mailles de fracture que l' on centre sur des n_uds soit aux intersections de fractures soit aux extrémités des fractures, chaque n_ud étant associé à un bloc de matrice regroupant tous les points plus proches de lui que de n_uds vo1sins; g) 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 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; h) on détermine les flux directs entre les mailles de fracture; i) on détermine les flux directs entre les volumes de matrice au travers des arêtes communes des mailles; et j) on simule les interactions entre les variations de la pression et du débit observables dans au moins un puits au travers du milieu. 2) Méthode selon la revendication 1, caractérisce en ce que l'on détermine le bloc de matrice associé à chaque maille de fracture en discrétisant chaque couche fracturce 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. 3) Méthode selon la revendication 2, caractérisce en ce que 1'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