FR 2809494 A1 20011130 FR 0006874 A 20000526 La présente invention concerne une méthode pour modéliser des écoulements de fluides dans un milieu fissuré traversé par des fractures relativement grandes, au sens que l'on va préciser ci-après. La méthode peut être mise en oeuvre par exemple dans le domaine de la production pétrolière par des ingénieurs de gisement en vue d'obtenir des prédictions d'écoulement fiables dans des réservoirs pétroliers comportant des discontinuités structurales ou sédimentaires dont les propriétés d'écoulement sont contrastées par rapport au milieu encaissant. Par cette modélisation, il s'agit de simuler des écoulements dans un milieu poreux perméable (réservoir) traversé par un réseau de fissures et/ou de couches fines beaucoup plus conductrices que la matrice poreuse. Les réservoirs fissurés constituent un type extrême de réservoirs hétérogènes 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 un milieu fissuré représentant moins de 1% de l'huile en place et hautement conducteur. Le milieu fissuré lui-même peut être complexe, avec différents ensembles de fissures caractérisés par leur densité, longueur, orientation, inclinaison et ouverture respectives. Dans un réservoir en production, les conditions de pression et de débit imposées aux puits producteurs et injecteurs induisent des écoulements des fluides en place dans le réservoir (huile, gaz et eau). La simulation de ces écoulements consiste à déterminer l'évolution des pressions et des saturations dans le réservoir au cours du temps. De par la nature très hétérogène des réservoirs fracturés, les fluides se déplacent relativement vite au travers du réseau de fractures et beaucoup plus lentement dans la matrice. La simulation d'écoulements dans le contexte d'un réservoir fracturé nécessite donc une très bonne représentation des hétérogénéités majeures que constituent les fractures. La précision de cette représentation dépend du type de simulation mise en oeuvre et de la taille des fractures à modéliser. Ainsi, pour une simulation très précise d'un essai de puits, on sera amené à représenter la géométrie exacte du réseau de fractures. A I l'inverse, pour la prise en compte d'un réseau dense de petites fractures à l'échelle d'un champ, on cherchera une représentation simplifiée équivalente. La présente invention permet de simuler précisément les écoulements à l'échelle d'un champ en présence de nombreuses grandes fractures. ETAT DE LA TECHNIQUE Actuellement dans l'industrie pétrolière notamment, on utilise les modèles à double porosité (double milieu) pour simuler les écoulements dans des milieux fracturés, mais ils 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 Warren et Root, dans à une modélisation explicite des deux milieux. Par le brevet FR-A-2.757.947 (US 6 023 656) du demandeur, on connaît une technique de détermination de la perméabilité de fissure équivalente d'un réseau de fissures dans un milieu multi-couche souterrain à partir d'une représentation connue de ce réseau. Elle permet 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. Par le brevet FR-A-2.757.957 du demandeur, on connaît une technique pour obtenir une modélisation simplifiée d'un milieu géologique hétérogène poreux (tel qu'un réservoir traversé par un réseau irrégulier de fissures 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 98/15.727 du demandeur, on connaît également une méthode pour modéliser les flux de fluides dans un milieu poreux multi-couches fissuré en tenant compte de la géométrie réelle du réseau de fissures et des échanges locaux entre la matrice poreuse et les fissures en chaque noeud du réseau. On discrétise le milieu fissuré par un maillage, on centre les mailles de fissure sur des noeuds aux différentes intersections des fissures, chaque noeud étant associé à un volume de matrice, et l'on détermine les flux entre chaque maille de fissure et le volume de matrice associé dans un régime de flux semi-permanent. Il existe des cas o les techniques précédentes sont difficiles à mettre en oeuvre, o l'on se trouve en présence d'un milieu traversé par de grandes fractures ou failles sub-sismiques dont le comportement hydraulique n'est pas homogénéisable à l'échelle de la maille. Une modélisation explicite de ces objets est alors nécessaire a priori mais leur nombre trop élevé rendrait une telle approche prohibitive à l'échelle d'un champ (nombre de mailles trop élevé et contraintes numériques). Le même problème se pose pour des réservoirs renfermant des bancs minces et très perméables dont le comportement s'apparente à celui de grandes fractures horizontales. La méthode selon l'invention La méthode de modélisation selon l'invention permet de simuler des écoulements de fluides dans un milieu géologique poreux fissuré de structure connue que l'on discrétise par un maillage et que l'on modélise en considérant que chaque volume élémentaire du milieu géologique fissuré est constitué d'un milieu fissure et d'un milieu matrice équivalents à l'échelle de chaque maille entre lesquels on détermine les échanges de fluides, ce milieu géologique étant traversé par un réseau de d'objets conducteurs de fluides de géométrie définie mais non homogénéisables à 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.). La méthode comporte une détermination des échanges intervenant entre le milieu matrice et le milieu fissure, 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 cet objet conducteur. 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. Dans toutes les mailles traversées par des objets conducteurs géométriquement définis (blocs matriciels de formes et de tailles irrégulières), on détermine un milieu transposé comportant un ensemble de blocs régulièrement disposés et séparés par un maillage régulier de fissures donnant sensiblement la même fonction de récupération de fluide lors d'un processus d'imbibition capillaire que le milieu réel. On calcule la dimension verticale des blocs du milieu transposé à partir des positions des couches sédimentaires très perméables dans la maille et on obtient les dimensions horizontales des blocs de ce milieu transposé, à partir d'une image à deux dimensions (2D) du milieu géologique sous la forme d'une série de pixels, par: - la détermination, pour chaque pixel, de la distance minimale qui le sépare de la fissure la plus proche; - la formation d'une distribution du nombre de pixels par rapport à la distance minimale au milieu fissuré et la détermination, à partir de cette distribution, de la fonction de récupération (R) dudit ensemble de blocs; et - la détermination de dimensions des blocs réguliers équivalents du milieu transposé à partir de la récupération (R) et de la récupération du bloc équivalent. Présentation des figures D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après d'un exemple non limitatif de réalisation, en se référant aux dessins annexés o: - la Fig.1 montre schématiquement deux mailles adjacentes d'une même couche d'un maillage de réservoir o la présence de niveaux sédimentaires fins et très perméables induit une transmissivité horizontale dans le milieu, fissure ", - la Fig.2 montre un maillage de réservoir traversé par un réseau de fissures; - la Fig.3 montre des mailles O, A, B, C d'un réservoir traversée par une fissure que l'on modélise par des transmissivités fissure- fissure; - la Fig.4 illustre le mode de calcul de la transmissivité entre deux mailles A et B traversées par une fissure; - la Fig.5 montre par comparaison le trajet en marches d'escalier au travers de mailles traversées par une fracture oblique, que l'on prend en compte pour les besoins de la simulation, et dont la transmissivité, selon le mode de calcul choisi, est cependant équivalente à la transmissivité réelle directement le long de la fracture; et - la Fig.6 illustre le mode de calcul de la taille d'un bloc équivalent en fonction du nombre de niveaux sédimentaires très perméable traversant une maille; - la Fig.7 présente un exemple de pixels voisins servant au calcul de la valeur attribuée à un pixel; et - la Fig. 8 présente une variation possible d'une zone envahie normée en fonction de la distance aux fissures. DESCRIPTION DETAILLEE On va considérer ci-après l'exemple d'un réservoir poreux traversé par un réseau de fractures F (Fig.2) que l'on suppose verticales pour simplifier et par des niveaux sédimentaires fins (sub- horizontaux) L (Fig.1) dont les propriétés pétrophysiques (perméabilité notamment) sont contrastées par rapport à l'encaissant matriciel. Ce réservoir est modélisé sous forme de deux grilles ,. Les échanges entre les mailles matrice sont calculées de manière standard bien connue des gens de l'art. I-Transmissivités entre mailles > 1-1 Transmissivités associées aux niveaux sédimentaires fins et perméables Les niveaux sédimentaires fins et perméables sont inclus dans le milieu > du modèle double milieu. Dans une maille traversée par de tels niveaux sédimentaires, les propriétés pétrophysiques de ces niveaux (porosité, perméabilité, saturation en eau) sont affectées au milieu fissure de la maille et les propriétés du reste de la roche contenue dans la maille sont affectées au milieu matrice de la même maille. La présence de niveaux sédimentaires fins et très perméable dans deux mailles adjacentes induit une transmissivité horizontale dans le milieu "fissure ", entre les deux mailles (transmissivité , induite puisque les niveaux sont horizontaux. Dans cet exemple, la transmissivité horizontale fissure " entre les mailles i et i+1 est calculée comme suit: Ks.Esii1+l AY Tsij, - AX - Ks est la perméabilité des niveaux sédimentaires très perméables, *Y est la taille des mailles en Y, *X est la taille en X, et Esi,+, est l'épaisseur des contacts entre niveaux sédimentaires fins et très perméables des deux mailles adjacentes i et i+1. Cette épaisseur est nulle si les niveaux des deux mailles ne sont i5 pas connectés. A l'inverse, s'ils sont totalement connectés, elle peut être égale à la plus petite des épaisseurs cumulées de niveaux dans les deux mailles. 1-2 Transmissivités associées aux fractures Le réseau de fractures verticales est aussi pris en compte dans le milieu , de l'approche double milieu. Dans chaque couche du modèle réservoir (Fig.2), ce réseau peut être représenté par une série de segments fissure qui traversent le maillage réservoir, comme sur le schéma suivant: Les données relatives à ces segments fissure sont: - la perméabilité fissure Kf, I'épaisseur de la fissure Ef, et - la longueur de fissure Lf. La porosité de fissure.f, calculée par la formule suivante Cf Lf Ef AX AY La communication entre les mailles du réservoir à travers le réseau de fractures est modélisée par des transmissivités fissure- fissure. Sur l'exemple de la Fig.3, des transmissivités fissure-fissure sont calculées entre les mailles traversées par la fracture, c'est à dire pour les couples O et A, A et Bet Bet C: TFOA, TFAB et TFBC. Cet exemple montre que le chemin réel d'un flux à travers une is fracture peut être éloigné du trajet imposé par la modélisation qui passe par les centres des mailles O, A, B et C. La solution qui consisterait à remplacer le segment de fracture par une ligne brisée passant par les centres de mailles conduirait à une mauvaise simulation des flux à travers ces mailles. Aussi, on détermine la transmissivité horizontale " fissure- fissure >, entre deux mailles traversées par la même fracture (cf. Fig.4), comme suit: AZ TfAB = Kf Ef - o: Lab Kf est la perméabilité intrinsèque de la fracture, Ef est l'épaisseur de la fracture, ÀZ est l'épaisseur de la couche, Lab est la longueur du segment ab, a est le milieu du segment de fracture traversant la maille A, et b est le milieu du segment de fracture traversant la maille B On note que plus la longueur Lab est petite, plus la transmissivité entre les mailles A et B est grande. Ainsi, on corrige numériquement l'effet de grille qui impose un chemin en marches d'escalier au flux. Ainsi le flux entre deux mailles distantes reliées par une fracture peut être correctement simulé en dépit du trajet en escalier (Fig.5) imposé par la grille. En effet, sur l'exemple de la Fig. 5, on peut vérifier que: y1 _ 1 o les Tf sont les transmissivités, fissure-fissure " entre TT f fi TfMN mailles réservoir cartésiennes le long du trajet en marches d'escalier reliant M et N, et TMN est la transmissivité réelle de la fracture entre M et N. Concernant la perméabilité verticale, fissurefissure " induite par la présence d'une fracture traversant plusieurs couches du réservoir, on a: Lf TfV =Kf Ef -- oU: AZ Kf est la perméabilité intrinsèque de la fracture, Ef est l'épaisseur de la fracture, *Z est l'épaisseur de la couche, Lf est la longueur du segment de fracture dans la maille considérée. Enfin, quand plusieurs fractures traversent une maille, les transmissivités calculées pour ces fractures prises individuellement sont ajoutées. 1-3 Transmissivité résultante Les transmissivités relatives aux niveaux sédimentaires très perméables (Ts) et les transmissivités relatives aux fracture (T) sont calculées séparément suivant les méthodes exposées précédemment. Les transmissivités sont calculées simplement par addition des T, et Tf. De même, la porosité finale du milieu " fissure " dans chaque maille est la somme des porosités dues aux niveaux sédimentaires très perméables d'une part et aux fractures d'autre part. Il- Dimensions de bloc matriciel équivalent 11-1 Dimensions horizontales Les dimensions horizontales du bloc équivalent (a,b) sont contrôlées par les fractures verticales présentes dans le réservoir. En effet, les échanges, matrice-fissure, dues aux fractures verticales sont uniquement horizontaux. Dans chaque maille traversée par au moins une fracture, ces dimensions horizontales sont calculées par la méthode décrite dans le brevet FR 2. 757.957 précité. Pour les mailles dans lesquelles il n'y a pas de fracture, les dimensions horizontales de bloc équivalent sont infinies. Pratiquement, on leur affecte une valeur très grande (par exemple 10 km). Selon cette méthode on détermine les dimensions de bloc équivalent par identification des comportements du milieu fracturé réel et du milieu équivalent pour un mécanisme d'imbibition diphasique eau-huile. Cela consiste en la mise en adéquation de la fonction de 1o récupération d'huile R(t) (du milieu fracturé réel) obtenue par une méthode de traitement d'image (décrite plus bas) et de la fonction de récupération Req(t) du milieu équivalent dont l'expression analytique est connue et dépend des dimensions du bloc équivalent. 11-1-a) Formulation géométrique Les fissures étant définies par les coordonnées de leurs points limites sur une coupe bidimensionnelle XY d'une couche, le processus d'imbibition par lequel de l'eau est présente dans les fissures et de l'huile est présente dans les blocs matriciels doit être déterminé. On suppose que l'invasion de la matrice par de l'eau est du type piston. On20 admet que la fonction x=f(t) qui relie l'avancée du front d'eau au temps est la même pour tous les blocs matriciels, quelle que soit leur forme, et pour tous les blocs élémentaires. Par conséquent, la mise en adéquation des fonctions R(t) et Req(t) équivaut à la mise en adéquation des fonctions R(x) et Req(x). Ces fonctions définissent physiquement des zones normées envahies par l'eau en fonction de l'avancée du front d'imbibition dans le milieu fissuré. En 2D, I'expression analytique de Req(x) est la suivante: I{Req(x) = 1(a-2x)(b-2x) 2a+ a xb x x 0, min 2 Req(x) = 1, x> min4 2'J o a et b sont les dimensions du bloc rectangulaire ou carré équivalent (a et b > 0): La fonction R(x) n'a pas d'expression analytique. Elle est calculée à partir d'une discrétisation de la coupe XY de la couche étudiée suivant l'algorithme défini ci-après. 11-1-b) Algorithme de calcul de la fonction R(x) La coupe XY de la couche étudiée est considérée comme une image dont chaque pixel représente un élément de surface. Ces pixels sont espacés régulièrement d'un pas dX dans la direction X et dY dans la direction Y (Fig.7). L'algorithme mis en oeuvre vise à déterminer, pour chaque pixel de cette image, la distance minimale qui le sépare de la fissure la plus proche. L'image est traduite par un tableau de nombres réels à deux dimensions: Pict[0:nx+1,0:ny+1] o nx et ny sont les nombres de pixels de l'image dans les directions X et Y. En pratique, le nombre total de pixels (nx.ny) est par exemple de l'ordre du million. Les valeurs des éléments du tableau Pict sont les distances recherchées. Initialisation: Tous les pixels par lesquels passe une fissure sont à une distance nulle de la fissure la plus proche. Pour ces pixels, le tableau Pict est donc initialisé à la valeur O. Ceci est fait par un algorithme connu en soi (l'algorithme de Bresline par exemple) auquel on donne les coordonnées des pixels correspondants aux deux extrémités d'une fissure considérée comme un segment de droite et qui initialise (à 0 dans le cas présent) les pixels les plus proches. Les autres éléments de Pict sont initialisés à une valeur supérieure à la plus grande distance existant entre deux pixels de l'image. Cette valeur est par exemple nx.dX +ny.dY. Calcul: Pour un pixel donné, le calcul de la distance recherchée 1o à la fissure la plus proche se fait à partir des valeurs de distance déjà calculées pour les pixels voisins. On lui affecte une valeur qui, si elle s'avère inférieure à la valeur qui lui a été assignée initialement, est le minimum des valeurs des pixels voisins auxquelles on ajoute la distance de ces pixels à celui considéré. Ce calcul est réalisé en deux phases successives. Lors de la passe descendante, on parcourt l'image ligne par ligne, de haut en bas et de gauche à droite (de Pict [1,1] à Pict [nx,ny]). Les pixels dont on tient compte sont différents selon que l'on est dans une passe descendante ou une passe montante. Comme le montre la Fig. 7, les pixels en noir et en grisé sont ceux que l'on prend en compte respectivement durant les passes descendantes et les passes ascendantes, pour le pixel Px. L'écart oblique dxy étant défini comme: dxy = dx2 +dy2, l'algorithme s'écrit: pour j=1 aà ny I pour i=1 à nx I I Pict[i,j] = min Pict[i-1l,j] + dx,: passe descendante I IPict[i-1,j-1] + dxy, 1I Pict[i,j-1] + dy, I IPict[i+1,j-1] + dxy, I IPict[i,j] I fin de boucle sur i fin de boucle sur j pour j=ny à 1, I pouri=lxà 1, I I Pict[i,j] = min Pict[i+1,j] + dx, : passe descendante I IPict[i+1,j+l] + dxy, I I Pict[i,j+1] + dy, I Pict[i-1,j+1] + dxy, II Pict[i,j] I fin de boucle sur i fin de boucle sur j Histogramme: A partir du tableau Pict ainsi calculé, on peut tracer un histogramme en classant les valeurs non nulles (celles affectées aux pixels hors des fissures) par ordre croissant. Le cumulé de cet histogramme donne, pour toute distance délimitant deux intervalles de l'histogramme, le nombre de pixels non nuls dont la valeur est inférieure à cette distance. Dans l'application décrite à un milieu poreux fissuré o cette distance correspond à l'avancée du front d'eau, le cumulé de l'histogramme indique donc l'aire envahie par l'eau. La courbe R(x) est obtenue en divisant ce cumulé par le nombre total de pixels non nuls (pour la normer). Le nombre d'intervalles utilisés en abscisse pour l'histogramme correspond au nombre de points de discrétisation de la courbe R(x). On le choisit égal à 500 par exemple. 11-1-c) Recherche des dimensions du bloc équivalent A ce stade, on connaît la fonction R(x) et l'on cherche les paramètres (a,b) (dimensions du bloc équivalent qui minimisent la fonctionnelle): N to J(a,b) =,(R(xi)- Req- a,b,xi))2 i=l o N est le nombre de points de discrétisation de R(x) et (xi) sont les abscisses de ces points de discrétisation. Discrétisation suivant l'ordonnée de R(x): Pour donner autant de poids à tous les volumes d'huile récupérés lors de l'imbibition, la courbe R(x) est re-discrétisée suivant un pas constant sur l'axe des ordonnées (Fig. 8). La suite (xi) utilisée par la fonctionnelle est déduite de cette discrétisation. Minimisation de la fonctionnelle: Comme a et b jouent des rôles symétriques dans l'expression Req(a,b,x), on utilise en fait la fonctionnelle N J(u, v) = (R(xi) - Rê q(,v,ixi)) i =] Req(u,v,x)= uxx+vxb avec RE q(u,v,x) axb Minimiser cette fonctionnelle revient à trouver le couple (ii,) pour lequel J'(,v) = 0. Ceci est fait par un algorithme de Newton. Ensuite, le couple (a,b) recherché est déduit de (i,v). Trois cas peuvent se présenter: 1) v > 0 signifie qu'une des valeurs du couple (a,b) est négative, ce qui n'a aucun sens physique. On pose alors v=0 dans l'expression de Réq(u,v,x), ce qui implique que les fissures sont parallèles. L'opération est recommencée et le couple (, b) est calculé comme suit 2 2 aU lbinf ini 2) Le cas U2+4v résultants de ces niveaux sont en effet uniquement verticaux. Dans chaque maille traversée par au moins un niveau très perméable, la dimension verticale du bloc équivalent est calculée par la to formule suivante AZ C=- Ns o *Z est l'épaisseur de la maille et Ns est le nombre de niveaux sédimentaires très perméables distincts dans la maille. Par exemple, sur le schéma suivant, Ns vaut 2: Pour les mailles dans lesquelles il n'y a pas de niveaux fins très perméables, la hauteur de bloc équivalent est infinie. Pratiquement, on leur affecte aussi une très grande valeur (par exemple 10 km). REVENDICATIONS 1) Méthode de modélisation permettant de simuler des écoulements de fluides dans un milieu géologique poreux fissuré de structure connue que l'on discrétise par un maillage et que l'on modélise en considérant que chaque volume élémentaire du milieu géologique fissuré est constitué d'un milieu fissure et d'un milieu matrice équivalents à l'échelle de chaque maille entre lesquels on détermine les échanges de fluides, ce milieu géologique étant traversé 1o par un réseau d'objets conducteurs de fluides de géométrie définie mais non homogénéisables à l'échelle de chaque maille du modèle, comportant une détermination des échanges intervenant entre le milieu matrice et le milieu fissure, caractérisée en ce que l'on modélise les 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 cet objet conducteur. 2) Méthode selon la revendication 1, caractérisée en ce que, dans le cas o les objets conducteurs sont des couches sédimentaires très perméables, on attribue à la transmissivité entre mailles adjacentes traversées par chaque couche très perméable une valeur dépendant des dimensions des mailles, de l'aire de contact commune entre couches à la jonction des mailles adjacentes et de la perméabilité des couches très perméables. 3) Méthode selon la revendication 1, caractérisée en ce que, dans le cas o les objets conducteurs sont des fractures, on attribue à la transmissivité entre mailles adjacentes traversées par chaque fracture, une transmissivité dépendant des dimensions des mailles, de l'épaisseur de la fracture et de sa perméabilité intrinsèque. 4) Méthode selon l'une des revendications précédentes, caractérisée en ce que dans les mailles traversées par des objets conducteurs géométriquement définis, on détermine un milieu transposé comportant un ensemble de blocs régulièrement disposés et séparés par un maillage régulier de fissures, ledit milieu transposé donnant sensiblement la même récupération (Req) de fluide lors d'un processus d'imbibition capillaire que le milieu réel, la dimension verticale des blocs du milieu transposé étant déterminée à partir des positions des couches sédimentaires très perméables dans la maille et les dimensions horizontales des blocs de ce milieu transposé étant obtenues, à partir d'une image à deux dimensions (2D) du milieu géologique sous la forme d'une série de pixels, par: - la détermination, pour chaque pixel, de la distance minimale le séparant de la fissure la plus proche; - la formation d'une distribution du nombre de pixels par rapport à la distance minimale au milieu fissuré et la détermination, à partir de cette distribution, de la fonction de récupération (R) dudit ensemble de blocs; et - la détermination de dimensions (a,b) des blocs réguliers équivalents du milieu transposé à partir de la récupération (R) et de la récupération (Req) du bloc équivalent.