FR2850773A1

FR2850773A1 - Fluid e.g. oil, flow in multilayer porous medium modeling method, involves discretizing each fractured layer by grid pattern, and calculating local flows between each fracture grid cell and associated matrix volume in pseudo-steady state

Info

Publication number: FR2850773A1
Application number: FR0301090A
Authority: FR
Inventors: Remy Basquet; Laurent Jeannin; Bernard Bourbiaux
Original assignee: IFP Energies Nouvelles IFPEN
Current assignee: IFP Energies Nouvelles IFPEN
Priority date: 2003-01-31
Filing date: 2003-01-31
Publication date: 2004-08-06
Anticipated expiration: 2023-01-31
Also published as: FR2850773B1

Abstract

Modeling compressible fluid flows in a multilayer porous medium involves discretizing each fractured layer by a grid pattern comprising fracture rigid cells that are centered on nodes each node being associated with a matrix block, and calculating the local flows between each fracture grid cell and associated matrix volume in pseudosteady state. Modeling compressible fluid flows in a multilayer porous medium crossed by a network of fractures of given geometry and by a well involves: (a) discretizing each fractured layer by a grid pattern comprising fracture rigid cells that are centered on nodes at the fracture intersections or at the fracture ends, each node being associated with a matrix block including all the points that are closer than to neighboring nodes; (b) calculating the local flows between each fracture grid cell and associated matrix volume in a pseudosteady state, the transmissivity value being obtained by considering a linear variation of the pressure as a function of the distance between any point of the block and the fracture grid cell; (c) determining the direct flows between the fracture grid cells; (d) determining the direct flow between the matrix volumes through the common edges of the grid cells; (e) determining the direct flows between the fracture grid cells and the matrix grid cells on the one hand, and the volume of the well on the other hand; and (f) simulating the response of the well for imposed flow rate variations by using a diffusion equation that connects the interactions between the pressure and flow rate variations that can be observed in the well.

Description

La présente invention concerne une méthode pour modéliser des flux deThe present invention relates to a method for modeling fluxes of

fluidesfluid

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. highly compressible such as gases in a porous multilayer medium crossed by a network of fractures of given geometry, unevenly distributed in the medium, and some of which communicate with each other.

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. The method according to the invention finds applications in particular in the field of the production of hydrocarbons. It allows reservoir engineers to test by simulation the production of a deposit by one or more wells crossing a porous permeable underground zone comprising two contrasting media, a matrix medium generally containing the largest part of the fluid in place and having a low permeability, and, crossing it in part or in whole, networks of fractures much more conductive in most cases.

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). STATE OF THE ART During a gas well test, the flow conditions imposed on the well cause the gas contained in the tank to flow towards the well. It is a single-phase, mono-component flow (the only mobile phase is gas) and highly compressible (the rock in the reservoir and the fluid are compressible).

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. For any elementary volume of the tank, the gas pressure contained in this volume then obeys the following equation: D CTa = (VK) + PQ (1) with: 2: pore volume CT: total compressibility (fluid + rock) K : permeability of rock g: viscosity of fluid Z: compressibility factor of gas Q: incoming flow, and iv: pseudopressure function defined by yv = 2 | P -dp o p is the pressure of the fluid Po The viscosity and the compressibility factor of the gas are variable as a function of the pressure and are given for a certain number of values in a table called PVT (Pressure Volume Temperature). From this table we can also deduce the function 10 pseudo-pressure defined above and the compressibility of the gas deduced from the equation of state.

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). The PVT table is therefore a 5-column table (pressure, pseudo-pressure, viscosity, compressibility factor and compressibility of the gas) from which, at a given pressure, the corresponding pseudo-pressure, viscosity and compressibility can be obtained by linear interpolation (conversely, from a pseudo-pressure you can obtain the other 4 data by interpolation).

Le débit entrant Q est nul partout sauf aux endroits o le puits est en communication avec le réservoir. The incoming flow Q is zero everywhere except where the well is in communication with the reservoir.

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. To simulate a well test, in any environment, we must solve this equation in space and in time. One thus carries out a discretization of the reservoir (mesh) and the resolution of the problem consists in finding the pressures of the meshes in the course of time, itself discretized in a certain number of time steps.

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). There are known tools for modeling single-phase flows which are applied not to the real geological medium in all its complexity but to a homogenized representation, according to the so-called double medium reservoir model described for example by 25 Warren and Root, in "The Behavior of Naturally Fractured Reservoirs ", SPE Journal September 1963. Any elementary volume of the cracked reservoir is thus modeled as a set of identical parallelepipedic blocks bounded by an orthogonal system of continuous uniform fractures oriented in the direction of one of the three main directions of flow (model called "sugar box"). The flow of fluids on the reservoir scale takes place mainly through the crack medium and exchanges of fluids occur locally between the cracks and the matrix blocks. This representation, which does not account for the complexity of the network of fractures in a reservoir, turns out to be effective nonetheless, but on the scale of a reservoir mesh (typically 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é. However, it is preferable that reservoir engineers have a flow simulator based on a "real" geological model of the environment and not on an equivalent homogeneous model, this tool making it possible to: * validate the geological image of the reservoir constructed by the geologist from all the information he has been able to collect on the fracturing of the reservoir (this validation is carried out by comparison with actual well test data); * reliably test the sensitivity of the hydraulic behavior of the medium to uncertainties on the geological image of the cracked medium.

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. By the patents FR 2 787 219 of the applicant, a method is known for modeling the flows of weakly compressible fluids (oil) in a cracked multilayer porous medium taking into account the real geometry of the network of fractures and local exchanges within the porous matrix and between the porous matrix and the fractures at each node of the network. One discretizes the cracked medium by a mesh, one associates the meshs of fracture with nodes placed either at the intersections of the fractures or at the ends of the fractures. Each node is associated with a volume or block of matrix, the fluxes between each fracture mesh and the volume of matrix are determined in a semi-permanent flow regime. In the case where the fracture network is weakly connected, that is to say for which the fracture network does not cross the whole of the matrix volume considered, the direct flows between the matrix volumes are also determined through the edges common meshes. By taking into account the transmissivities between matrix blocks associated with nodes of discretization of the medium, one can simulate in a more realistic way the response of a well to variations of imposed flow.

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. The object of the method according to the invention is to model flows of highly compressible fluids in a porous multilayer medium traversed by a network of fractures of given geometry and by a well. It includes the following stages: a) one discretizes each fractured layer by a mesh comprising meshes of 5 fracture which one centers on nodes either at the intersections of fractures or at the ends of the fractures, each node being associated with a block of matrix grouping all the points closer to it than from neighboring nodes; b) the local fluxes between each fracture mesh and the associated matrix volume are calculated in a pseudo-permanent flux regime, the value of the transmissivity being obtained by considering a linear variation of the pressure as a function of the distance between any point block and fracture mesh; c) determining the direct fluxes between the fracture meshes; d) the direct flows between the volumes of matrix are determined through the common edges of the meshes; e) direct fluxes between the fracture meshes and the matrix meshes, and on the other hand the volume of the well, are determined from transmissivity values adapted to turbulent flows; and f) simulating the response of the well for variations in flow imposed by means of a diffusion equation modified to take account of the compressibility of the fluid, which links the interactions between the variations in pressure and flow observable in the well.

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. According to one mode of implementation, the diffusion equation is modified by the introduction of a term depending on the pressure, the viscosity and a compressibility factor of the compressible fluid.

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é. According to one mode of implementation, direct flows are determined by introducing a difference term proportional to the imposed fluid flow rate.

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. According to an embodiment, the matrix block associated with each fracture mesh is determined by discretizing each fractured layer into a set of pixels and by calculating the distance from each pixel to the nearest fracture mesh, to determine the 'location of the edges between the meshes and one uses the distances calculated to deduce therefrom the value of transmissivities between each mesh of fracture and the volume of matrix associated on the one hand, and between the common edges of the meshes.

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. For example, the position of the edges between the meshes is identified by moving two elongated windows in the image formed by the pixels, oriented in two perpendicular directions.

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. Thus, taking into account the high compressibility of the gas and adapting the transmissivities between the matrix blocks and the well and between the cracks and the well, it is possible to simulate realistically the response of a gas well to imposed flow variations.

Présentation sommaire des figuresSummary presentation of the 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. Other characteristics and advantages of the method according to the invention will appear on reading the description below, with reference to the accompanying drawings: o - Figure 1 shows an example of a 2D fracture mesh; - Figure 2 shows an example of a matrix block obtained in 2D; FIG. 3 shows an example of two matrix blocks MB communicating via a fracture element FE; - Figure 4 shows an example of two matrix blocks MB each crossed by a fracture element and communicating by means of exchanges through their common edges A; - Figure 5 shows examples of windows oriented along two perpendicular axes, allowing the identification of the edges between the matrix blocks; - Figure 6 shows a well W intersected by two cracks; - Figure 7 shows an example of variation as a function of time t of flow D imposed in a well test; - Figure 8 shows two types of linear (LF) and radial (RF) flow from the matrix to the well; - Figure 9 shows a first part of the flow diagram for implementing the modified modeling in the case considered where the fluids are highly compressible and - Figure 10 shows a complementary part of the flow diagram for implementing the modified modeling in the case considered where fluids are highly compressible.

DESCRIPTION DETAILLEEDETAILED DESCRIPTION

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é. In the case of a cracked medium, if one wants to take into account all the complexity of the network of fractures, it is necessary to carry out an explicit mesh of this network. Furthermore, in such an environment, the K permeabilities are generally much higher in the fractures than in the matrix and the flows are consequently rapid in the fractures and slow in the matrix. It is therefore a question of representing the cracked medium using the concept of double medium (matrix and crack) described in the aforementioned patent FR 2 787 219.

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. The stages of the approach adopted in the case where the fluid produced is oil are the following: * Explicit mesh of the network of fractures * association with each fracture mesh of a single matrix block; * treatment of the flows between each fracture mesh and its associated block in pseudo-permanent regime o the flow from one to the other is proportional to the pressure difference between the two; * taking into account flows between matrix blocks; and * dynamic simulation of flows within the fracture network and the matrix environment.

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. In the case where the fluid produced is highly compressible, typically a gas, 25 we specify or modify: * the system of equations of modeling of flows to find a differential equation of diffusion of the same form as for weakly compressible flows; * the terms describing the flows of fluids between the wells mesh and the adjacent fracture mesh as well as between the wells mesh and the adjacent matrix mesh; and * dynamic simulation to update during the simulation the viscosity and total compressibility parameters at each iteration in time by an interpolation technique, from a table of so-called PVT data.

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. The unknowns are the pressures of the fracture meshes and the pressures of the associated matrix blocks. Since there are as many matrix blocks as there are fracture meshes, the number of unknowns is equal to twice this number.

* 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). * The diffusion equation The diffusion equation that we solve is deduced from the mass conservation equation, from Darcy's law and from an equation of state which depends on the compressibility of the mobile fluid. We consider two different equations of state depending on whether the fluid is weakly compressible (oil, case a) or highly compressible (gas, case 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) at. For weakly compressible fluids (oil) The equation of state translating the equivalent compressibility of the mobile fluid (oil) is as follows: Ch = ap (2) PP) T Considering the system composed of the three equations or law that we have just to quote (equations of conservation of mass and state, Darcy's law), for any elementary volume of the reservoir, the pressure of the oil contained in this volume then obeys the following differential equation: (* CT P = div (- .Vp) + Q (3) ot t with:): porosity CT: total compressibility (fluid + rock) K: permeability of rock gt: viscosity of the fluid Q: incoming flow p: density of the fluid p: unknown pressure 10 b. For highly compressible fluids (gas) In the case of gas, the density which varies considerably depending on the pressure is expressed as follows: pM (4)

ZRTZRT

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é. with M molar mass of the gas, R constant of the ideal gases and Z the compressibility factor of the gas which translates the deviation of the real gas compared to the ideal gases. The compressibility of the gas can then be written: cg P Za (5) PZ [0 PiT The differential diffusion equation then becomes: 20 div kgradp + ICt Z - = q'z (6) / LZ tZ pto Z with C , the overall compressibility of a unitary element of the saturated pore volume.

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). To find the form of a more classical diffusion equation; this equation is expressed using a pseudo-pressure function defined by: y =: 2 {-Ldp (7) Po t The diffusion equation expressed in pseudopression is finally written: k V = 2 pq ( 8) The differential equation to be solved then presents a form similar to that used in the case of weakly compressible fluids, the unknown being pseudopressure Y- (connected to pressure by the PVT table).

* 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. * Discretization of the fracture network In 2D, the calculation nodes are positioned at the intersections between fractures and at the ends of the fractures. The needs of 3D modeling lead to adding additional nodes in the upper and lower layers for fractures crossing several layers.

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. In the mesh of the present method, the nodes of calculation thus defined constitute the centers of the meshs of fracture. In addition, the simulation of highly compressible flows requires knowing the volume of the meshes (p in equation 8). One thus defines limits of meshes positioned in the middle of the segments connecting the nodes of computation. The diagram in FIG. i show an example of a 2D fracture mesh. In terms of number of computation nodes, the mesh of the fracture network is thus made 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é. * Crack-to-crack transmissivities The calculation of the connections between fracture meshes (transmissivities), used for the calculation of the fluxes between these fracture meshes, is identical to that presented in the method 25 described in the aforementioned patent FR 2 757 947.

* 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é. * Discretization of the matrix medium According to the principle of the method, the discretization of the matrix medium consists in assigning a volume of rock to each mesh of fracture defined as indicated in the preceding paragraph. During the dynamic simulation, the matrix-fracture exchanges 5 will be calculated in permanent pseudo between each fracture mesh and its single associated matrix block.

Pour le calcul de ces volumes de matrice, on traite le problème couche par couche. For the calculation of these volumes of matrix, one deals with the problem layer by layer.

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. For a given layer, the blocks are defined as follows: * the block heights are equal and fixed at the height of the layer; * in the layer plane, the surface of the block associated with a fracture mesh is the set of points closer to this fracture mesh than another.

Physiquement, cette définition suppose que les limites à flux nul dans la matrice sont les lignes d'équidistance entre les mailles de fracture. Physically, this definition assumes that the zero flux limits in the matrix are the equidistance lines between the fracture meshes.

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é. To determine the geometry of the blocks in a given layer, it is therefore a question of solving a two-dimensional problem which consists in finding, for each fracture mesh, the points closest to this mesh than another. This problem is solved by using the geometric method exposed in the aforementioned patent FR 2047 957.

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. The method in question makes it possible, by discretizing the fractured layer into a set of pixels and by applying an image processing algorithm, to determine the distance of each pixel to the nearest fracture. It is recalled that in the initialization phase of this algorithm, the value 0 (zero distance) is assigned to the pixels belonging to a fracture and a large value to the others. If in this phase, each fracture pixel is given the number of fracture cells to which it belongs, the same algorithm finally makes it possible to determine, for each pixel: * the distance from this pixel to the nearest fracture cell; * the number of the nearest fracture mesh.

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. The set of pixels having the same number of fracture mesh constitutes the surface of the matrix block associated with this mesh. By multiplying this surface by the height of the layer, we obtain the volume of the block associated with the mesh. An example of a matrix block thus obtained in 2D is shown diagrammatically in 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. By construction, the sum of the surfaces of the blocks is equal to the total surface of the layer, which guarantees the conservation of the volume of the layer.

* 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. * Matrix-rissure transmissivities For each matrix block, the image processing algorithm also gives the distance of each pixel of the block to the associated fracture mesh. This information is used to calculate the transmissivity between the fracture mesh and the matrix block.

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. The hypothesis of a pseudo-permanent regime between a fracture mesh and a matrix block amounts to considering that the flow from one to the other is proportional to the difference in the pressures of the one and the other. We have the following relation: Fmf = - (Pm-pf) (9) with: Fn: matrix-fracture flux Tmf: matrix-fracture transmissivity g: viscosity Pm: pressure of the matrix block pf: pressure of the associated fracture mesh In assuming the pseudo-permanent regime, the value of the transmissivity Tmf is constant during the simulation. In the present method, this value is calculated for each of the fracture mesh - matrix block pairs.

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. For this computation, one makes the assumption that in the matrix block, the pressure varies linearly according to the distance from the point considered to the mesh of fracture associated with the block. We define the transmissivity as follows 21 H- K Tmf D (10) with: l: length of the fractures in the fracture mesh H: height of the layer (and of the matrix block) K: permeability of the matrix D: distance of the fractures at which the pressure is the average pressure of the block We know 1, H and K (the factor 2 comes from the fact that the fracture has two faces). The calculation of D is made from information given by the image processing algorithm on the distances from the pixels to the nearest fractures. Indeed, according to the hypothesis of linear variation of the pressure as a function of the distance to the fracture, we have: 1 - (1 D = Jd (s) ds (11) o S is, in 2D, the surface of the block matrix and d (s) is the function giving the distance to the fractures from the points of the matrix block to the associated fracture mesh.

En terme de pixels, on a donc simplement: D= cIdn (12) In terms of pixels, we therefore simply have: D = cIdn (12)

N NN N

o N est le nombre de pixels du bloc matriciel et dj la distance du pixel n à la maille de 20 fracture. o N is the number of pixels in the matrix block and dj the distance from pixel n to the fracture mesh.

* 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. * Matrix-matrix transmissivity The calculation of exchanges between matrix blocks is based on a numerical diagram of 2-point flux approximation. This approach assumes that the flow between two matrix blocks is perpendicular to the edge separating these two blocks. This numerical diagram is therefore generally used for orthogonal meshes such as Vorono meshes.

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. In our problem, the geometry of the matrix blocks is arbitrary. In addition, the presence of fractures leads to consider two possible geometric cases.

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. In the first case (Fig. 3), the edge between the two matrix blocks is crossed by a fracture. In the second case, it is not.

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. In the second case (Fig. 4), we assume that in each of the two fractures of the problem, the pressure can be considered constant before the pressure variation in the associated matrix block (because the conductivity of the fracture is very large compared to the permeability of the matrix). A two-point transmissivity between the two fracture segments is calculated through the edge separating the two matrix blocks.

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. This calculation is done by an image processing algorithm which uses the digital image defined during the determination of the geometry of the matrix blocks and which therefore informs, for each pixel of the image, on the distance from this pixel to the closest fracture and on the matrix block to which this pixel belongs.

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. The image processing algorithm consists in locating the elementary interfaces between blocks, that is to say the edges separating two pixels associated with two different blocks.

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. The algorithm calculates, for each of these interfaces, an elementary transmissivity between the two blocks and updates the total transmissivity between these blocks according to the following formula: T = k MH p2- (di -dfn) (13) 2 - 2di, where Tmm is the matrix-matrix transmissivity between two matrix blocks, M the number of elementary interfaces between the two blocks, p the size of the side of the pixels and di, and df1 the distances of the two ends of the interface n to nearest fractures. As we can see, the total transmissivity is a sum of elementary transmissivities calculated for each of the interfaces between the two blocks.

La Fig.5 montre un ensemble de pixels appartenant aux mailles notées 1, 2 et 3. Fig.5 shows a set of pixels belonging to the meshes noted 1, 2 and 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. To locate the interfaces between the pixels of these meshes, the image processing algorithm moves two windows made up of 2 times 1 pixel along the lines and columns of the image. There is a vertical window for locating the horizontal interfaces and a horizontal window for locating the vertical interfaces.

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) In the case where a fracture crosses the edge between two matrix blocks, the transmissivity between these two blocks is calculated simply with the following expression: T. = .km H * L (14)

LABLAB

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. where I is the length of the edge and LAB the length of the fracture connecting the centers A and B. The fracture length is known by construction and the length of the edge is calculated by the image processing algorithm from the elementary interfaces between the matrix blocks of centers A and B. * Representation of wells A well is represented by its geometry on the one hand and the flow constraints which are imposed on it on the other hand.

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. A well is made up of a set of connected segments which intersect network fractures (Fig. 6). The geometric representation of a well is therefore a broken line in 3D. The communication points between the well and the network are the points of intersection of this line with the 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. The flow constraints are imposed at the point where the well enters the fractured medium. They are entered by the user in the form of a plateau function giving the flow as a function of time. The well flow change times indicate the different periods to be simulated. For a conventional well test, a flow is imposed during a first period at the end of which the well is closed (zero flow) and the pressure rise in the well is observed. The flow curve to be given is shown diagrammatically in 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. * Well-reservoir connection a) Crack-well transmissivity The calculation of the exchanges between the well and a fracture is based on the pseudopermanent regime hypothesis, that is to say to consider that the flux from one to the other is proportional unlike the pressures from one to the other. We have the following relation: Tp Fpf = f * (p- pf) (15) with Fpf: flow well-fracture Tpf: transmissivity well-fracture u: viscosity pp: pressure of the well pf: pressure of the associated fracture mesh In the pseudopermanent regime hypothesis, the value of the transmissivity Tp, is constant during the simulation. In the present method, this value is calculated for each of the fracture mesh - well segment pairs.

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. Darcy's law used for the establishment of the differential equation which governs the flow of fluids in the porous medium, assumes that the flow regime is laminar. In the case of gases and in the vicinity of the well, the speed of the fluid increases sharply so that one leaves the laminar regime (turbulent regime). In this specific case, Darcy's law is not applicable. To circumvent this problem, it has been shown that by introducing a term called deviation from Darcy's law (dependent on the imposed flow rate) to the formulation established for laminar regimes, the flow of the fluid is correctly represented. This term of deviation from Darcy's law takes the following form: S '= S0 + D.Q (16) S' has the dimension of an additional pressure drop. So represents the 5 modifications of the hydraulic properties around the well (damage during drilling, injection of acids into the well) and the coefficient D, which is a characteristic characteristic of the turbulent flow regime.

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. In the present approach, we introduce these coefficients So and D in the expression of the crack-well transmissivity used for the representation of the flux between the crack and the connected well.

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. In the case where the intersection between a crack and the well is substantially circular, the transmissivity is defined as follows: 27rCf (17) 0ln () + SO + D (q. + Ql with: r0 = 0. 14 (1f) 2 + (h) 2, for isotropic media Cf = kf ef, kf being the intrinsic permeability of the fracture and ef its opening 1f: the fracture length h: the thickness of the layer qn: the flow rate of the well segment at the flow step 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. in the case where the crack-well intersection is arbitrary (elliptical), the calculation of the transmissivity is deduced from the generalized Peaceman function well known to those skilled in the art.

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. b) Matrix-well transmissivity The calculation of the exchanges between the well and an associated matrix mesh is also based on the pseudopermanent regime hypothesis. To express these exchanges, we combine the formulations relating to a radial flow around the well and a linear flow towards the well (Fig. 8). The transmissivity is then defined as follows: Tmw = 02lrkml (18) ln (2h) + 2X (lm -h) + S + Dq (tj) dw, ho: li: the length of the well segment lm: the distance from the well to which the pressure is the average pressure of the matrix block 10 * Dynamic simulation During dynamic simulation, the simulated time period is broken down into time steps whose duration is included for a well test between lx10-3 s (just after l opening or closing of the well) and a few hours.

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. The passage of the instant n (for which we know all the depression values) and the instant n + l, is carried out by solving the following equation in all the meshes and blocks of the reservoir. This equation corresponds to discretized equation 8. Contrary to the case where the fluid is weakly compressible, the equation that we present is a nonlinear differential equation insofar as the compressibility and the viscosity of the fluid 20 vary according to the pressure. The discretized differential equation is written: nv + l -n (TAp n D .. (CT). Ij +. (Jn + l - n + y = 2 Zn Qi (19) t n + 1 t "n [ n J * + with i: number of the mesh or of the block j: numbers of the neighbors of i tn: pseudopression at the instant n pn: pressure at the instant n tn time at the instant n (CT), r ' : the total compressibility of the fluid for a pressure at time n + 1 IL ': the viscosity of the fluid for a pressure at time n + 1 at the interface between the meshes i and j zn: the compressibility factor of the fluid for a pressure at time n / n: the viscosity of the fluid for a pressure at time n Qi: flow entering i at the bottom of the well Tij transmissivity of connection between i and j To solve this nonlinear equation, we use a classic 15 Newton method which consists in linearizing equation (19) by expressing: n + l, k + l = n + l, k ji Yi + 8yi (20) The technique then consists, at each iteration k, in solve the linear system described below such so that the value of the unknown 6Ei becomes negligible. After _.l n + 1, k + 1 convergence with k + 1 iterations, the variable Y Y + 1 + k + i is determined and by linear interpolation on the table PVT, all the mesh pressures, compressibility and viscosity fluid will be refreshed.

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. Constitution of the system to be solved Calculation of the contribution of the fracture meshes and matrix accumulation term: A (CT) VfnOIk + (dCT, tn + l _ tn I n + Ik (VJi -'lYi link between fracture meshes: second member of the system: p Qi tz VI 'Yj _____ ((N + ik yn +, k) - /'% + hIk ó, r Resolution of the linear system Resolution by the iterative conjugate gradient algorithm (CGS) well known to people of the art with updating of the pressures of the fracture meshes and of the matrix blocks Management of the time step In a well test simulation, the time steps are generally very short 15 after the opening or the closing of the well , because it is at these moments that the pressure variations are strong. Later, the time steps can be longer, when the pressure variations are weaker.

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- (<1) A partir de ces valeurs, la gestion du pas de temps est réalisée comme suit: A l'instant tn, on calcule le majorant des variations de pressions dans les mailles de 10 fracture et les blocs matriciels (c'est à dire entre les instants tnet tn). Suivant la valeur de cette variation de pression dp, le pas de temps est soit: * augmenté d'un rapport rdt+ si dp < dpmin, * diminué d'un rapport rdt- si dp > dpmax, * inchangé dans les autres cas. Managing the time step consists in relating the duration of the time step n + 1 to the maximum pressure variation observed during the time step n. To do this, user 20 must provide the following data: Initial time step: dtini Minimum time step: dtmin Maximum time step: dtmax Lower pressure variation terminal: dpmin Upper pressure variation terminal: dpmax Factor d time step increase: rdt + (> 1) Time step reduction factor: rdt- (<1) From these values, the time step is managed as follows: At time tn, we computes the majorizing of the variations of pressures in the meshs of fracture and the matrix blocks (ie between the instants tnet tn). Depending on the value of this pressure variation dp, the time step is either: * increased by a ratio rdt + if dp <dpmin, * decreased by a ratio rdt- if dp> dpmax, * unchanged in the other cases.

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). After each change of flow condition at the well, the time step is reset to the value dtini. In addition, an optimization is carried out at the end of the simulation period: if tf is the time to be simulated to reach the end of the period and dt the expected time step, the time step chosi is 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é. Input data of the simulation The geometry of the fracture network as well as the attributes of the fractures (conductivity, opening) are given in a file in the format defined in the aforementioned patent FR 2 757 947.

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é. The petrophysical data to be provided are the following (the unit is indicated in parenthesis): - * compressibility of the network of fractures (Cf in bar "l) * compressibility of the matrix (cm in barf ') * permeability of matrix (K in mD) * matrix porosity (Jm dimensionless) The data in question are considered to be homogeneous on the cracked medium considered.

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. The data concerning the fluids contained in the cracked medium relate to the gas (mobile) and the water (stationary) which is always present in residual quantity in the matrix: * residual saturation in water in the matrix (Sw without dimension) * compressibility of water (c, on lease) * a PVT table including the compressibility factor of the gas (Z) and the viscosity of the gas (it in cpo) as a function of the pressure.

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. The data concerning the well are: * geometry in the form of a series of connected segments * imposed flows in the form of a curve giving the imposed flow as a function of time.

* le coefficient So d'endommagement de la formation autour du puits, et le coefficient D d'écart à la loi de Darcy. * the coefficient So of damage to the formation around the well, and the coefficient D of deviation from Darcy's law.

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-. The values to provide are the values already defined in the chapter on the management of the time step: dtini, dtmin, dtmax, dpmin, dpmax, rdt + and rdt-.

L'organigramme de la figure 9 récapitule les différentes étapes de la modélisation réalisée. The flowchart in Figure 9 summarizes the different stages of the modeling performed.

Claims

1) Method for modeling flows of compressible fluids in a multilayer porous medium crossed by a network of fractures of given geometry and by a well, comprising the following stages a) one discretizes each fractured layer by a mesh comprising meshs of fracture that l 'we center on nodes either at the intersections of fractures or at the ends of the fractures, each node being associated with a matrix block grouping all the points closer to it than from neighboring nodes; b) the local fluxes between each fracture mesh and the associated volume of matrix 10 are calculated in a pseudo-permanent flux regime, the value of the transmissivity being obtained by considering a linear variation of the pressure as a function of the distance between any point of the block and the fracture mesh; c) determining the direct fluxes between the fracture meshes; d) direct flows between the matrix volumes are determined through the common edges of the meshes; e) direct fluxes between the fracture meshes and the matrix meshes, and on the other hand the volume of the well, are determined from transmissivity values adapted to turbulent flows; and f) simulating the response of the well for variations in flow imposed by means of a modified diffusion equation to take account of the compressibility of the fluid, which links the interactions between the variations in pressure and flow observable in the well.

2) Method according to claim 1, in which the diffusion equation is modified by the introduction of a term depending on the pressure, the viscosity and the compressibility factor of the compressible fluid.

3) Method according to claim i or 2, in which the diffusion equation is modified by the introduction of a pseudopressure taking into account the variation in viscosity and compressibility of the fluids as a function of the pressure.

23 2850773 - -23

4) Method according to one of the preceding claims, characterized in that the direct flows are determined by introducing a term of deviation from Darcy's law proportional to the imposed flow of fluids, reflecting the pressure drops located in the vicinity of the wells linked to the high speed of said fluids.

5) Method according to one of the preceding claims, characterized in that the matrix block associated with each fracture mesh is determined by discretizing each fractured layer into a set of pixels and by calculating the distance from each pixel to the mesh of nearest fracture, to determine the location of the edges between the meshes and the distances calculated are used to deduce therefrom the value of the transmissivities 10 between each fracture mesh and the volume of associated matrix on the one hand, and between the common edges of the meshes.

6) Method according to claim 5, characterized in that one locates the position of the edges between the meshes by moving two elongated windows in the image formed by the pixels, oriented in two perpendicular directions.