Titre : PROCEDE DE GENERATION DE MAILLAGE HEXAEDRIQUE Domaine technique

[0001] La présente divulgation concerne un procédé de génération d’un maillage hexaédrique structuré à partir d’un modèle géologique. Elle trouve des applications privilégiées pour l’étude de la migration des fluides, de la géomécanique, du dépliage des unités géologiques ainsi que de la diffusion de la chaleur d’un domaine géologique, e.g. un réservoir géologique d’hydrocarbures, de ressources géothermiques ou de stockage de C02

Technique antérieure

[0002] Pour réaliser des simulations numériques, comme des simulations d’écoulements des fluides, dans un réservoir représenté par un modèle géologique, il est connu de réaliser un maillage de ce modèle géologique qui peut ensuite être peuplé de paramètres pétrochimiques, pétrophysiques, et mécaniques permettant de réaliser les simulations numériques.

[0003] Pour faciliter la mise en oeuvre de ces simulations, il est préférable de représenter le modèle géologique par un maillage hexaédrique structuré, c’est-à-dire un maillage dans lequel chaque cellule présente la forme d’un hexaèdre, et dont chaque cellule (sauf les cellules situées aux limites du modèle) a un nombre constant de cellules voisines, ou autrement dit dont chaque sommet d’une cellule appartient à un nombre constant de cellules, typiquement 8. De cette façon, chaque cellule d’un maillage structuré peut être simplement indexée par un tableau en trois dimensions (U,V,W).

[0004] En référence à la figure 1a, on a représenté un exemple simple d’un modèle géologique à représenter, comprenant des horizons H, et des failles F. Pour ce type de modèle, il est connu de réaliser un maillage hexaédrique structuré dont les failles et les frontières du domaine d’étude sont approximées comme s’étendant uniquement selon l’une de deux directions U et V orthogonales, comme représenté sur les figures 1 b et 1c. Sur la figure 1 d, une troisième direction W orthogonale aux deux premières est définie par l’empilement des horizons, de sorte que chaque horizon s’étend selon une surface iso-W. Le maillage obtenu est représenté sur la figure 1e.

[0005] Ce maillage n’est cependant pas applicable à des cas plus complexes, et en particulier lorsqu’il existe des intersections entre des failles, ou entre les failles et les frontières du modèle, ou encore en cas de failles mourantes, c’est-à-dire une faille qui se termine avant d’atteindre une frontière du domaine. On a représenté un tel exemple sur la figure 2 : la faille représentée en pointillés peut être approximée par une faille s’étendant selon la direction U, mais l’autre faille ne peut alors ni être approximée par une faille parallèle (puisque par définition elle est sécante avec la première faille) ni par une faille s’étendant selon la direction V puisqu’elle est également sécante avec la frontière du domaine s’étendant selon la direction V.

[0006] Pour réaliser une représentation précise de ce type de modèle géologique, plusieurs approches ont été proposées. On connaît par exemple des approches, comme par exemple dans le document US 8,150,663, visant à approximer une faille par une représentation en « marches d’escalier » combinant des faces s’étendant selon U et selon V. Ceci permet d’obtenir un maillage structuré qui est compatible avec la plupart des logiciels de simulation d’écoulement actuellement sur le marché. Cependant ce type de représentation n’est pas approprié pour la mise en oeuvre de calculs de simulation mécanique, qui requièrent une grande fidélité dans la représentation des failles afin d’appliquer des contraintes de limite. Donc lorsqu’il est souhaité de réaliser à la fois des calculs d’écoulement et des calculs mécaniques (c’est-à-dire, des calculs couplés), il est nécessaire de réaliser deux maillages, l’un comme indiqué précédemment pour les calculs d’écoulement et un autre, plus fidèle sur la représentation des failles, qui peut par exemple être tétraédrique (et donc non structuré) pour les calculs mécaniques. Cela consomme deux fois plus de mémoire de l’ordinateur. Le couplage nécessite un transfert aller et retour des propriétés pétrophysiques et mécaniques entre les deux maillages, qui est sujet à erreur et chronophage.

[0007] D’autres approches utilisent un maillage mixte qui est majoritairement hexaédrique dans les zones du modèle géométrique exemptes de failles, et qui est localement non- hexaédrique au niveau des failles, en remplissant les zones du maillage adjacentes aux failles de polyèdres comme dans US10, 114,134 ou d’autres éléments non-hexaédriques comme des tétraèdres ou des prismes, dans US 2017/0184760.

[0008] L’inconvénient de ce type d’approche est qu’elle n’est pas compatible avec certains simulateurs d’écoulement disponibles sur le marché qui ne peuvent fonctionner ou ne sont optimaux qu’avec des maillages hexaédriques structurés. Le même problème se pose aussi pour des maillages intégralement non-structurés, comme par exemple des maillages dits de Voronoi, décrits dans Palagi C. L. et Aziz, K : « Use of Voronoi Grid in Réservoir Simulation », SPE Advanced Technology Sériés 2 (02): 69-77. SPE-22889-PA

Résumé

[0009] La présente divulgation vient améliorer la situation.

[0010] En particulier, un but de l’invention est de proposer une méthode permettant de réaliser un maillage hexaédrique structuré à partir de modèles géologiques pour lesquels les solutions actuelles ne fournissent que des maillages localement non structurés ou intégralement non structurés.

[0011] Un but de l’invention est de réaliser un maillage hexaédrique structuré à partir de modèles géologiques comprenant des failles mourantes, ou des intersections sensiblement verticales entre deux failles ou entre une faille et une frontière du domaine.

[0012] Un autre but de l’invention est de permettre de réaliser un maillage majoritairement hexaédrique structuré à partir d’un modèle géométrique pouvant aussi comprendre des intersections sensiblement horizontales entre deux failles ou entre une faille et une frontière du domaine. [0013] A cet égard, il est proposé un procédé de génération d’un maillage hexaédrique à partir d’un modèle géologique délimité par un ensemble de frontières et comprenant un ensemble d’horizons et un ensemble de failles, le procédé étant mis en oeuvre par un calculateur et comprenant : la génération d’un maillage tétraédrique en trois dimensions à partir du modèle géologique, le maillage tétraédrique respectant un ensemble de contraintes surfaciques correspondant aux failles, aux horizons et aux frontières du modèle, l’obtention d’un maillage hexaédrique en trois dimensions à partir du maillage tétraédrique, comprenant : o l’application d’une première déformation au maillage tétraédrique pour rendre les failles et les frontières verticales, o le calcul d’une paramétrisation globale en deux dimensions (U,V) d’un horizon de référence du maillage tétraédrique déformé, o la détermination d’une grille en deux dimensions formée de quadrilatères représentant l’horizon de référence à partir de la paramétrisation globale selon deux axes (U,V) dudit horizon de référence, o l’obtention d’un maillage hexaédrique en trois dimensions par propagation de la grille représentant l’horizon de référence aux autres horizons du modèle géologique, et o l’application, au maillage hexaédrique obtenu, d’une déformation inverse à la première déformation pour restaurer la géométrie initiale des failles, et si le maillage hexaédrique obtenu comprend au moins un point singulier, la mise en oeuvre d’un traitement complémentaire du maillage pour rendre le maillage structuré au niveau de chaque point singulier. [0014] Dans des modes de réalisation, la propagation de la grille en trois dimensions comprend : le calcul d’une paramétrisation en W du maillage tétraédrique en trois dimensions de sorte que chaque horizon corresponde à une iso-surface selon la coordonnée w, et la détermination d’une grille en deux dimensions à base de quadrilatères représentant chaque iso-surface selon la coordonnée W à partir de la grille représentant l’horizon de référence, et la construction d’un ensemble d’hexaèdres entre deux iso-surfaces selon la coordonnée W consécutives, de sorte que deux faces opposées d’un hexaèdre soient formées respectivement par un quadrilatère de chaque iso-surface, ledit ensemble d’hexaèdres formant le maillage hexaédrique.

[0015] Dans des modes de réalisation, la détermination d’une grille en deux dimensions à base de quadrilatères représentant une iso-surface selon la coordonnée W est réalisée par une projection verticale, sur ladite iso-surface, de la grille représentant l’horizon de référence.

[0016] Dans des modes de réalisation, le calcul de la paramétrisation globale en (U,V) de l’horizon de référence est mis en oeuvre de sorte que les bords de l’horizon de référence coïncide avec des lignes d’iso-valeurs entières de U et V, et la détermination de la grille en deux dimensions formée de quadrilatères représentant l’horizon de référence est mise en oeuvre en formant un ensemble de quadrilatères dont les sommets sont définis par des points d’intersection de lignes d’iso-valeurs entières de U et V.

[0017] Dans des modes de réalisation, le calcul d’une paramétrisation globale selon deux axes (U,V) d’un horizon de référence est mis en oeuvre de manière à ce que chaque cellule de l’horizon de référence soit paramétrée par rapport à deux axes U et V, les directions des axes U et V du repère étant variables sur l’horizon de référence.

[0018] Dans des modes de réalisation, la paramétrisation globale (U,V) de l’horizon de référence associe à chaque bord commun entre deux triangles voisin une fonction de transition établissant une relation entre les coordonnées (Uj, Vj) du bord du premier triangle et les coordonnées (Uk,Vk) du bord de l’autre triangle telle que : (Uj,Vj) = R ^a (Uk, Vk) + (b,c), où R ^a = est la rotation du repère de multiple de 90 degrés, (b,c) est la translation des coordonnées , et a, b, c sont tous des valeurs entières.

[0019] Dans des modes de réalisation, la génération d’un maillage tétraédrique à partir du modèle géologique est mise en oeuvre par tessellation de Delaunay contrainte. [0020] Dans des modes de réalisation, le traitement complémentaire du maillage pour rendre le maillage structuré comprend :

Le partitionnement du maillage hexaédrique obtenu en une pluralité de sous- maillages hexaédriques structurés associés chacun à une indexation locale (U,V,W) respective,

L’intégration de chaque sous-maillage hexaédrique structuré dans un maillage hexaédrique structuré associé à une indexation (U,V,W) globale et, L’identification de l’ensemble des paires de mailles voisines dans le maillage hexaédrique obtenu initialement et qui ne le sont plus dans le maillage hexaédrique structuré d’indexation globale.

[0021] Dans des modes de réalisation, le partitionnement du maillage hexaédrique est mis en oeuvre par : le partitionnement de la grille en deux-dimensions à base de quadrilatères représentant l’horizon de référence en sous-grilles, par des lignes de découpe partant d’un point singulier et s’étendant le long d’un bord d’un quadrilatère de la grille dont le point singulier est un sommet, et la propagation selon l’axe W de chaque ligne de découpe pour obtenir un plan de découpe du maillage hexaédrique.

[0022] Dans des modes de réalisation, le procédé peut également comprendre la mise en oeuvre, sur le maillage obtenu, d’un traitement complémentaire parmi le groupe suivant :

Modélisation par des méthodes géostatistiques des propriétés physiques,

Simulation d’écoulement de fluide,

Simulation de la géomécanique,

Simulation de la chaleur,

Simulation chemique,

Simulation multiphysique couplée,

Simulation cinématique.

[0023] Dans des modes de réalisation, le procédé est mis en oeuvre pour la représentation d’un modèle géologique comprenant en outre au moins un puits, dans lequel le puits est représenté par une ligne d’intersection entre deux surfaces traitées comme des failles virtuelles.

[0024] [0025] Dans des modes de réalisation, l’intégration de chaque sous-maillage hexaédrique structuré dans un maillage hexaédrique structuré associé à une indexation (U,V,W) globale comprend :

Le partitionnement d’une grille en deux dimensions représentant un horizon en une pluralité de sous-grilles, l’intégration de l’ensemble des sous-grilles dans une grille en deux dimensions associée à une indexation (U,V) globale, ladite intégration comprenant : o pour chaque sous-grille, la détermination d’un cadre en deux dimensions U et V englobant l’ensemble des indices occupés de la sous-grille dans son indexation locale, et la détermination des dimensions L axj et Vmaxj de ce cadre, o L’initialisation d’une grille en deux dimensions associée à une indexation globale (U,V) dont les dimensions en U et V sont respectivement la somme des dimensions en U et la somme des dimensions en V des cadres déterminés pour toutes les sous-grilles, o Le remplissage de la grille d’indexation globale par chaque sous-grille par ordre décroissant de taille de chaque sous-grille, tel que chaque sous-grille insérée dans la grille d’indexation globale n’occupe que des indices non occupés de la grille d’indexation globale, et qu’elle soit la plus proche possible des sous-grilles précédemment insérées, o Le redimensionnement de la grille d’indexation globale pour correspondre à un cadre minimum englobant toutes les sous-grilles, et

La propagation selon l’axe W de la grille d’indexation globale obtenue pour obtenir un maillage hexaédrique structuré intégrant l’ensemble des sous- maillages inclus dans le maillage hexaédrique initial.

[0026] Selon un autre objet, il est décrit un procédé de génération d’un maillage majoritairement hexaédrique, à partir d’un modèle géologique comprenant au moins deux failles sécantes présentant une intersection s’étendant selon une ligne sensiblement horizontale, et le procédé comprend : la mise en oeuvre du procédé selon la description qui précède à partir d’un modèle géologique modifié dans lequel une faille a été supprimée d’un ensemble de failles sécantes de manière à supprimer l’intersection sensiblement horizontale, pour obtenir un maillage hexaédrique, la suppression d’un ensemble de cellules du maillage hexaédrique situées à l’intersection avec l’emplacement de la faille supprimée, l’intégration de la faille supprimée dans le maillage hexaédrique et le remplissage des cellules supprimées par des cellules polyédriques. [0027] Il est également décrit un produit programme d’ordinateur, comprenant des instructions de code pour la mise en oeuvre de l’un des procédés décrit précédemment, lorsqu’il est mis en oeuvre par un ordinateur.

[0028] Il est également décrit un support non transitoire lisible par ordinateur sur lequel sont stockées des instructions de code pour la mise en oeuvre de l’un des procédés selon la description qui précède, lorsqu’elles sont exécutées par un ordinateur.

[0029] Selon un autre objet, il est également décrit un dispositif de génération d’un maillage hexaédrique à partir d’un modèle géologique, caractérisé en ce qu’il comprend un calculateur configuré pour la mise en oeuvre de l’un des procédés selon la description qui précède. [0030] Le procédé proposé permet d’obtenir un maillage hexaédrique structuré avec une indexation régulière à partir d’un modèle géologique y compris dans le cas où celui-ci comprend des failles mourantes ou des failles présentant une intersection sensiblement verticale.

[0031] En effet, le fait de verticaliser les failles, puis de réaliser une paramétrisation globale d’un horizon en deux dimensions du domaine qui est propagée aux autres horizons, permet de réaliser de façon efficace et robuste un paramétrage global de l’ensemble du domaine en trois dimensions, en produisant un minimum de points singuliers, et en respectant avec fidélité la géométrie des failles et des horizons. Le fait de disposer pour chaque horizon du même paramétrage en deux dimensions permet aussi de retraiter localement chaque point singulier pour rendre le maillage structuré au niveau dudit point et ainsi obtenir un maillage intégralement structuré.

[0032] Dans le cas de failles dites « en X » ou « en Y », c’est-à-dire de deux failles présentant une intersection sensiblement horizontale, ce procédé peut aussi être mis en oeuvre en ignorant l’une des failles de l’intersection, puis en complétant le maillage par des mailles polyédriques, par exemple tétraédriques ou hexaédriques au voisinage de la faille préalablement ignorée afin d’obtenir, pour ce type de modèle géométrique plus complexe, un maillage majoritairement hexaédrique structuré.

Brève description des dessins [0033] D’autres caractéristiques, détails et avantages apparaîtront à la lecture de la description détaillée ci-après, et à l’analyse des dessins annexés, sur lesquels :

Fig. 1a

[0034] [Fig. 1a], déjà décrite, représente un exemple de modèle géologique. Fig. 1b

[0035] [Fig. 1b], déjà décrite, représente des surfaces iso-u d’un paramétrage du modèle géologique.

Fig. 1c

[0036] [Fig. 1c], déjà décrite, représente des surfaces iso-v d’un paramétrage du modèle géologique.

Fig. 1d

[0037] [Fig. 1 d], déjà décrite, représente des surface iso-w d’un paramétrage du modèle géologique.

Fig. 1e [0038] [Fig. 1e], déjà décrite, représente un maillage hexaédrique structuré obtenu pour un modèle géologique simple dans lequel une faille peut être approximée comme étant parallèle à une frontière du modèle.

Fig. 2

[0039] [Fig. 2] déjà décrite également, illustre l'impossibilité d’appliquer le même type de maillage dans le cas de failles présentant une intersection.

Fig. 3a

[0040] [Fig. 3a] représente schématiquement une étape d’obtention d’un maillage tétraédrique d’un modèle géologique.

Fig. 3b [0041] [Fig. 3b] représente schématiquement une étape de verticalisation d’un maillage tétraédrique.

Fig. 3c

[0042] [Fig. 3c] représente schématiquement la représentation d’un horizon du modèle géologique par une grille de quadrilatères. Fig. 3d [0043] [Fig. 3d] représente l’obtention d'une grille de quadrilatères représentant un horizon à partir d’un paramétrisation globale en U et V.

Fig. 3e

[0044] [Fig. 3e] représente un maillage hexaédrique obtenu par propagation de la grille représentée en figure 3c.

Fig. 3f

[0045] [Fig. 3f] est une vue éclatée du maillage représenté sur la figure 3e.

Fig. 4a

[0046] [Fig. 4a] représente un exemple de paramétrisation d’une surface avec un repère en deux dimensions fixe sur toute la surface,

Fig. 4b

[0047] [Fig. 4b] représente un exemple de paramétrisation globale d’une surface avec un repère en deux dimensions, variable sur la surface.

Fig. 5a

[0048] [Fig. 5a] représente un exemple de partitionnement d’une surface d’un maillage hexaédrique localement non structuré.

Fig. 5b

[0049] [Fig. 6b] représente une intégration des sous-grilles obtenues sur la figure 4a dans une grille d’un espace d’indexation globale.

Fig. 6a

[0050] [Fig. 6a] représente un exemple de faille en X ou Y.

Fig. 6b

[0051] [Fig. 6b] représente un exemple d’application d’un procédé d’obtention d’un maillage majoritairement hexaédrique selon un mode de réalisation, à un modèle géologique comprenant des failles en X ou en Y, présentant une intersection sensiblement horizontale.

Fig. 7

[0052] [Fig. 7] représente schématiquement les principales étapes du procédé selon un mode de réalisation.

Fig. 8

[0053] [Fig. 8] représente un maillage d’un modèle géologique intégrant un puits. Description des modes de réalisation

[0054] En référence aux figures 3a à 3e et à la figure 7, on va maintenant décrire un procédé de génération d’un maillage hexaédrique à partir d’un modèle géologique. Ce procédé est mis en oeuvre par un calculateur, par exemple un processeur central ou un processeur graphique, par exécution d’une série d’instructions de code stockées dans une mémoire.

[0055] Le modèle géologique comprend un ensemble d’horizons H ou couches géologiques, les horizons étant, lors de leur formation, sensiblement horizontaux et superposés les uns sur les autres. En raison des déformations géologiques intervenues entre la formation des horizons et le présent, les horizons dans le modèle géologique ne sont pas nécessairement horizontaux. Le modèle géologique peut également comprendre des failles F s’étendant transversalement aux horizons et pouvant également induire une discontinuité de hauteur d’un horizon de part et d’autre de la faille.

[0056] Dans des modes de réalisation, le modèle géologique peut également comprendre un ou plusieurs puits P. Chaque puits peut être représenté par un ou plusieurs segments sensiblement rectilignes, chaque segment étant défini comme l’intersection de deux surfaces F _p qui sont traitées dans la suite du procédé comme des failles virtuelles. Un puit de géométrie simple est représenté en figure 8, ce puits comprenant un segment rectiligne vertical défini par l’intersection de deux surfaces sensiblement verticales qui sont traitées comme des failles. Ces surfaces virtuelles ne sont utilisées que pour la mise en oeuvre des contraintes et ne sont pas représentées dans le maillage final. Seule l’intersection des surface virtuelles approximant la trajectoire du puits est représentée dans le maillage final. Des puits plus complexes peuvent également être représentés, comprenant par exemple un segment sensiblement vertical et un ou plusieurs segments sensiblement horizontaux. [0057] Dans la suite, on considérera séparément les cas où : le modèle géologique ne comprend pas de failles, ou comprend des failles non intersectées ou dont l’intersection s’étend sensiblement verticalement, y compris des surfaces virtuelles F _p pour la représentation de puits ou de segments de puits sensiblement verticaux, et - le cas où le modèle géologique comprend des failles dont l’intersection s’étend sensiblement horizontalement, qui peuvent également être désignées par les termes de failles en X ou failles en Y, et/ou des puits ou segments de puits sensiblement horizontaux, formés par l’intersection de surfaces virtuelles en X ou en Y. [0058] Par ailleurs, le modèle géométrique est délimité par un ensemble de frontières (frontières artificielles qui délimitent la zone de modélisation).

[0059] En référence à la figure 7, dans le cas où le modèle géométrique comprend, comme sur l’exemple de la figure 5a, des failles sécantes dont l’intersection s’étend selon une ligne sensiblement horizontale, c’est-à-dire des failles en X ou en Y, ou des puits ou segments de puits sensiblement horizontaux, le procédé de génération d’un maillage hexaédrique à partir du modèle géologique comprend une étape préliminaire 90 de suppression temporaire d’au moins une faille d’un ensemble de failles sécantes de sorte qu’il n’y ait plus d’intersection horizontale entre deux failles. Cette étape permet de se ramener au premier cas décrit ci-avant, qui ne comprend pas de failles en X ou en Y ou de puits horizontaux, pour la mise en oeuvre de la suite des étapes.

[0060] Les étapes suivantes sont donc décrites en référence à un modèle géologique ne comprenant pas de failles en X ou en Y ou de puits horizontaux, soit parce que le modèle géologique correspond initialement à ce cas, soit parce qu’elles s’appliquent à un modèle géologique modifié suite à la mise en oeuvre de l’étape 90.

[0061] En référence à la figure 3a, le procédé de génération d’un maillage hexaédrique à partir du modèle géologique initial ou modifié selon l’étape 90 comprend une première étape 100 de génération d’un maillage tétraédrique T en trois dimensions à partir du modèle géologique, le maillage tétraédrique respectant les failles, les horizons et les frontières du domaine couvert par le modèle géologique.

[0062] Cette étape peut être réalisée classiquement en mettant en oeuvre une tessellation de Delaunay contrainte dans laquelle les failles, les horizons et les frontières du domaine sont traitées comme des contraintes surfaciques. Une tessellation correspond à une triangulation en dimension 3. A l’issue de cette étape, les failles, les horizons et les frontières du domaine sont représentées sous la forme de surfaces en deux dimensions triangulées.

[0063] Le procédé comprend ensuite des étapes permettant d’obtenir, à partir du maillage tétraédrique, un premier maillage hexaédrique du modèle. Dans un maillage hexaédrique, chaque point du maillage est relié à un nombre constant de huit cellules voisines (sauf aux bords). Un point singulier est un point qui est relié à un nombre de cellules voisines différent, qui peut être inférieur ou supérieur. Les étapes décrites ci-après permettent d’obtenir un maillage hexaédrique qui peut, en fonction de la complexité du modèle géométrique représenté, soit être structuré, c’est-à-dire ne comprendre aucun point singulier, soit comprendre, par la nature des étapes mises en oeuvre, un nombre réduit de points singuliers. [0064] En référence à la figure 3b, une étape 200 comprend la verticalisation des surfaces du maillage tétraédrique formées par les failles et par les frontières du domaine. Pour ce faire, pour chaque surface triangulée représentant une faille ou une frontière du domaine, une nouvelle position (Xv, Yv, Zv) est calculée pour chaque sommet d’un triangle de la surface, de sorte que la surface soit verticalisée. La nouvelle position est calculée en mettant en oeuvre un procédé d’optimisation numérique sous contrainte, minimisant une énergie définie comme une distance entre la géométrie initiale de la surface et la géométrie verticalisée de la surface, sous la contrainte que la géométrie verticalisée de la surface soit effectivement verticale, c’est-à-dire que la composante en Z de la normale à tous les triangles de la surface verticalisée soit nulle. La distance entre la géométrie initiale et la géométrie verticalisée peut par exemple être calculée entre la position de chaque point de la géométrie initiale et la position du même point dans la géométrie verticalisée.

[0065] Dans le cas où le sommet d’un triangle est partagé par plus d’une surface, ce qui est notamment le cas d’un sommet situé au niveau d’une faille, qui dans ce cas appartient à deux surfaces disposées de part et d’autre de la faille, ou d’un sommet situé au niveau d’une intersection entre une faille et une frontière du domaine, alors la position verticalisée de ce sommet est déterminée comme la valeur moyenne calculée pour la position de chaque surface à laquelle appartient ce sommet.

[0066] Les géométries verticalisées des surfaces représentant les failles et les frontières du domaine sont ensuite utilisées pour calculer une déformation W(T) de l’ensemble des positions des sommets du maillage tétraédrique T par interpolation. La déformation est calculée en mettant en oeuvre un procédé d’optimisation numérique sous contrainte, minimisant une distance entre la géométrie initiale du maillage T et sa géométrie déformée (ou verticalisée) W(T), en utilisant comme contraintes les positions (Xv, Yv, Zv) des sommets de chaque surface verticalisée. Ici encore, la distance entre la géométrie initiale et la géométrie verticalisée peut être calculée entre la position initiale de chaque sommet du maillage et la position du même sommet dans la géométrie verticalisée.

[0067] Le procédé comprend ensuite, le calcul d’une paramétrisation globale 300 en deux dimensions (U,V) d’un horizon de référence du maillage tétraédrique verticalisé W(T). Sur la figure 3b, on a représenté par W _ref l’horizon de référence. Sur cette figure, il s’agit de l’horizon correspondant à la surface supérieure du modèle, mais l’horizon de référence peut être quelconque.

[0068] L’horizon de référence W _ref est donc une surface en deux dimensions triangulée, comprenant un ensemble de triangles qui correspondent à des faces des tétraèdres du maillage tétraédrique verticalisé W(T). La frontière de cette surface 5W _ref comprend l’ensemble des sommets et des bords des triangles de cette surface à l’intersection avec des failles, d’autres horizons ou des frontières du domaine.

[0069] Le calcul de la paramétrisation globale en deux dimensions (U,V) est mis en oeuvre de manière à pouvoir paramétrer l’ensemble de l’horizon de référence avec un même repère (U,V). Cependant, contrairement à la représentation de la figure 4a, qui montre un exemple de paramétrisation d’une surface formée de quadrilatère, dans laquelle le repère (U,V) est fixe sur tout l’horizon, la paramétrisation globale mise en oeuvre à l’étape 300 définit un repère (U,V) dont la direction des axes U et V est variable sur l’horizon de référence. [0070] Comme représenté sur la figure 4b, qui représente un exemple de paramétrisation globale obtenue avec un repère (U,V) variable sur la même surface que la figure 4a, la paramétrisation globale est déterminée de sorte que, d’une maille à une maille adjacente, les repères sont à une rotation (par exemple multiple de 90 degrés) près. C’est à dire que les coordonnées U et V du repère peuvent être interchangée (U devient V et inversement) et / ou les directions des axes U et V du repère peuvent être inversées (par exemple : U devient -U). Ceci confère des degrés de liberté supplémentaires dans la paramétrisation de l’horizon de référence, qui permet de respecter des angles qui ne sont pas à angle droit, et en particulier de respecter la géométrie de failles ou d’intersection de failles avec les frontières du domaine qui ne sont pas à angle droit. [0071] La paramétrisation globale F en U et V obtenue est représentée sur l’horizon de référence W _ref par :

Une paire de coordonnées (U,V) à valeurs réelles associée à chaque sommet de chaque triangle de W _ref . Au sein de chaque triangle des coordonnées de n’importe quel point du triangle peuvent être obtenues par une interpolation linéaire à partir des coordonnées des trois sommets du triangle.

Une fonction de transition définie pour chaque frontière entre deux triangles voisins et établissant la relation entre les coordonnées paramétriques (Uj,Vj) du bord d’un triangle et celles (Uk,Vk) du même bord du triangle voisin : (Uj, Vj) = R ^a (Uk, Vk) + (b,c), où R ^a est la rotation du repère de multiple de 90 degrés, (b,c) est la translation des coordonnées et a, b, c sont tous des valeurs entières. Cette fonction de transition n’est définie que s’il ne s’agit pas de la fonction identité (i.e. a=0, b=0, c=0). A titre d’exemple non limitatif une telle fonction peut être définie, entre un triangle d’indice 1 et un triangle d’indice 2 par : {U2=V1 , V2 = -U 1 }, i.e. a=1 , b=0, c=0 [0072] Le calcul de la paramétrisation globale en deux dimensions (U,V) de l’horizon de référence, avec un repère (U,V) variable sur l’horizon de référence, peut être mis en oeuvre par application :

- de la méthode décrite dans la publication de Ray et al. : « Periodic Global Parameterization », ACM Transactions on Graphics, Association for Computing Machinery,

2006, ou

- de la méthode décrite dans la publication de D. Bommes et al. « Mixed-lnteger Quadrangulation », ACM SIGGRAPH 2009, en utilisant la frontière 5W _ref comme une contrainte telle que tous les bords de 5W _ref doivent coïncider avec des lignes d’iso valeurs entières des coordonnées U et V.

[0073] Une fois la paramétrisation globale en U et V obtenue pour l’horizon de référence représenté par un maillage de triangles, cette paramétrisation est utilisée pour mailler 400 le même horizon de référence avec des quadrilatères Qi, c’est-à-dire pour déterminer une grille en deux dimensions, formée de quadrilatères et représentant l’horizon de référence. Un exemple de mise en oeuvre de cette étape est représenté sur la figure 3c.

[0074] En référence à la figure 3d, une telle grille peut être formée par un ensemble de lignes d’iso valeurs de U et V, par exemple, des lignes correspondant à des iso valeurs entières de U et V. Etant donné que lors de la mise en oeuvre de la paramétrisation globale à l’étape précédente, la frontière 5W _ref de l’horizon de référence a été traitée comme une contrainte devant coïncider avec des lignes d’iso valeurs entières des coordonnées U et V, les bords des quadrilatères de la grille obtenue par les lignes d’iso valeurs entières de U et V correspondent nécessairement à la frontière 5W _ref de l’horizon.

[0075] Ainsi, lors de l’étape 400, on crée un sommet de la grille à chaque intersection entre deux lignes d’iso valeurs de U et V. Des bords sont créés pour relier deux sommets s’ils sont la même coordonnée sur un axe et une différence de coordonnée égale à un, en valeur absolue, sur l’autre axe. Les quadrilatères sont ensuite recensés à partir des sommets et des bords ainsi définis de manière à former une grille paramétrée en U et V. La taille des quadrilatères peut être déterminée par un utilisateur en fonction de la densité des lignes d’iso valeurs. Par exemple, des quadrilatères peuvent être définis par deux lignes consécutives en U ou en V, ou en variante un utilisateur peut fixer un pas plus élevé, qui sera, si les lignes d’iso valeurs sont à valeurs entières, un entier également.

[0076] Une fois que la grille est obtenue, elle peut être utilisée pour construire un maillage hexaédrique en trois dimensions sur l’ensemble du modèle géométrique. Préalablement, il faut d’abord calculer une paramétrisation 500 selon un troisième axe W perpendiculaire à U et V. Cette paramétrisation est de préférence calculée à partir du maillage tétraédrique obtenu à l’issue de l’étape 100, en associant une coordonnée en W à tous les sommets de ce maillage, et en utilisant les horizons comme contraintes géométriques. De cette façon, les surfaces à iso valeur de W correspondent naturellement à la géométrie des horizons dans l’espace initial. En variante, cette paramétrisation en W peut également être calculée sur la géométrie déformée W(T).

[0077] Pour ce faire, les horizons sont d’abord triés par ordre stratigraphique, c’est-à-dire dans l’ordre dans lequel les horizons géologiques ont été formés, et à chaque horizon est assignée une valeur réelle fonction de sa profondeur stratigraphique, tous les sommets des mailles appartenant à un même horizon étant assignés à une même valeur. La valeur de chaque horizon h peut par exemple correspondre à l’épaisseur cumulée des couches géologiques empilées jusqu’à cet horizon.

[0078] La paramétrisation est ensuite obtenue par interpolation à partir des valeurs des sommets des mailles appartenant aux différents horizons en résolvant un système d’équations linéaires, de préférence de façon numérique. Des surfaces intermédiaires d’iso valeurs de W peuvent être définies entre les horizons en fonction des préférences de l’utilisateur. Dans la figure 3c, l’exemple fourni ne comprend que cinq surfaces d’iso-valeurs de W qui correspondent aux cinq horizons du modèle.

[0079] Dans des modes de réalisations, l’étape 500 peut être mise en oeuvre immédiatement après l’étape 100, et les étapes visant à élaborer une grille de quadrilatères représentant un horizon de référence peuvent être mises en oeuvre après cette étape.

[0080] Le procédé comprend ensuite une étape 600 d’obtention d’un maillage hexaédrique à partir du maillage tétraédrique, en combinant le maillage en quadrilatères Qi, _re f de la grille représentant l’horizon de référence W _ref , et le paramétrage en W du maillage tétraédrique. Le maillage hexaédrique est obtenu en propageant la grille en deux-dimensions calculée sur un horizon de référence, sur toutes les autres surfaces d’iso valeurs de w, et notamment sur tous les autres horizons, dans le maillage tétraédrique verticalisé W(T). Etant donné que ce maillage a été verticalisé, toutes les failles et les frontières du domaine sont verticales, et il est donc possible de déterminer, pour chaque surface d’iso-valeur de W autre que l’horizon de référence, une grille en deux dimensions formée par des quadrilatères Qi représentant ladite surface, en calculant une projection de la grille calculée pour l’horizon de référence W _ref sur chacune des surfaces d’iso-valeurs de W. On obtient ainsi, pour chacune de ces surfaces, une grille formée de quadrilatères, et chaque grille présentant la même topologie que la grille représentant l’horizon de référence. [0081] En référence à la figure 3e, un maillage hexaédrique peut ensuite être obtenu en construisant un ensemble d’hexaèdres H, entre deux surfaces consécutives d’iso-valeur de W, de sorte que deux faces opposées d’un hexaèdre soient formées respectivement par un quadrilatère Q,, Q, _+i de chaque iso-surface. Ce maillage hexaédrique est obtenu à partir du maillage tétraédrique verticalisé, qui représente un modèle géométrique déformé dans lequel toutes les failles et les frontières du domaine sont verticalisées. Par conséquent, un maillage hexaédrique du modèle géométrique initial est obtenu en appliquant 700 au maillage hexaédrique obtenu la déformation inverse W ¹ de la déformation W ayant permis la verticalisation du maillage tétraédrique à l’étape 200. Un exemple d’un maillage hexaédrique obtenu à ce stade est représenté sur la figure 3e.

[0082] A l’issue de cette étape, on obtient un maillage hexaédrique MH qui, grâce aux étapes ayant permis son obtention, contient peu, voire pas du tout, de points singuliers, en fonction de la complexité initiale du modèle géométrique, et notamment du nombre et de la disposition des failles. En effet, le fait d’utiliser une paramétrisation globale en U et V pour mailler une surface de référence qui est ensuite projetée sur les autres surfaces dans une version verticalisée du domaine, permet d’assurer que tous les horizons et toutes les surfaces d’iso valeur de W ont le même nombre de points singuliers. De plus, la méthode mise en oeuvre pour obtenir la paramétrisation globale en U et V de la surface de référence permet elle-même de minimiser le nombre de points singuliers dans le maillage de cette surface. Comme visible sur la figure 3f, qui représente une vue éclatée du maillage de la figure 3e, la méthode permet également de mailler les frontières du domaine ainsi que les bords des failles en mailles hexaédriques.

[0083] Néanmoins, dans le cas où le maillage hexaédrique obtenu MH comporte tout de même des points singuliers, c’est-à-dire des points où le maillage est localement non structuré, le procédé comprend en outre un traitement complémentaire 800 du maillage, afin de rendre le maillage hexaédrique structuré en chacun de ces points (et cellules).

[0084] En référence à la figure 5a, ce traitement complémentaire 800 du maillage MH comprend le partitionnement 810 du maillage hexaédrique MH en plusieurs sous-maillages dont chacun peut être indexé localement par un triplet (Ui,Vi,Wi) respectif. Dans la suite, le terme « indexé » ou « indexation » désigne une paramétrisation où chaque sommet et cellule du maillage présente des coordonnées entières. Une telle indexation peut être déduite du maillage des horizons de l’étape 400, qui est obtenu à partir de lignes d’iso valeurs entières de U et V.

[0085] D’après la construction du maillage hexaédrique MH décrite précédemment, puisque toutes les surfaces d’iso valeur de W ont le même maillage en deux-dimensions (puisqu’il est obtenu par projection d’une grille d’un horizon de référence sur toutes les autres surfaces d’iso-valeur de W), ce partitionnement en trois dimensions peut être obtenu en partitionnant une grille en deux dimensions représentant un horizon ou toute autre surface d’iso valeur de W. [0086] Pour mettre en oeuvre le partitionnement en deux dimensions, on identifie d’abord tous les points singuliers de la surface considérés, c’est-à-dire, en deux dimensions, tous les points qui sont reliés à un nombre de cellules différent de 4 (sauf aux bords). Pour chaque point singulier, on définit un ensemble de lignes de découpe partant du point singulier et s’étendant selon un bord d’un quadrilatère dont le point singulier est un sommet, chaque ligne de découpe s’arrêtant à l’intersection avec une autre ligne de découpe, ou avec une faille ou encore un bord du domaine. Une ligne de découpes L1 est représentée sur la figure 4b, et d’autres lignes de découpe sont représentées en pointillés sur les figures 5a et 5b. Ces lignes de découpe partitionnent la surface considérée en un ensemble de grilles structurées formées de quadrilatères, ou chaque grille peut être aisément indexée par un couple de coordonnées locales [U, V]

[0087] Ces lignes de découpe, si elles sont propagées selon l’axe W, forment des plans de découpe partitionnant le maillage hexaédrique en sous-maillages structurés SMi qui peuvent tous être indexés par des coordonnées locales [Ui,Vi,Wi] respectives.

[0088] Une fois le partitionnement réalisé, en référence à la figure 4b, tous les sous- maillages SM, obtenus sont intégrés 820 dans un maillage hexaédrique MU d’indexation globale, comprenant un nombre de mailles supérieur au nombre de mailles cumulé total de tous les sous-maillages. L’intégration des sous-maillages dans ce maillage MU est faite de manière à préserver autant que possible les paires de mailles voisines dans le maillage hexaédrique initial MH. Pour les paires de mailles voisines dans le maillage hexaédrique initial et qui ne le sont plus dans le maillage global obtenu à l’issue de cette étape, ces paires sont néanmoins identifiées au cours d’une étape 830 pour pouvoir être indiquées explicitement à un simulateur d’écoulement, et qu’elles soient prises en compte comme des mailles voisines.

[0089] Dans un mode de réalisation particulier, l’intégration des sous-maillages hexaédriques dans un maillage d’indexation globale est mise en oeuvre de manière simplifiée en deux dimensions en intégrant chaque sous-grille obtenue à l’étape précédente dans une grille en deux dimensions d’indexation globale. Cette intégration est ensuite propagée en trois dimensions aux autres surfaces d’iso-valeur de W pour obtenir le maillage hexaédrique structuré global. [0090] L’intégration des grilles en deux-dimensions en une grille d’indexation globale peut être mise en oeuvre comme suit.

[0091] Pour chacune des sous-grilles, on détermine un cadre en deux dimensions englobant l’ensemble des indices occupés par une cellule de la sous-grille dans son indexation locale, et les dimensions de ce cadre déterminées pour obtenir les dimensions du cadre selon U et V : [I axj, Vmaxj]. I axj et Vmaxj correspondent donc à la taille du cadre.

[0092] Une grille en deux dimensions d’indexation globale est ensuite initialisée au cours d’une étape. Il s’agit d’une grille présentant une indexation en U et V qui sera commune à l’ensemble des sous-grilles à l’issue de l’étape d’intégration. Cette grille présente lors de son initialisation des dimensions en U et V correspondant respectivement à la somme des dimensions en U et à la somme des dimensions en V des cadres déterminés pour l’ensemble des sous-grilles. En notant NUi _nit et NVi _nit les dimensions initiales de cet espace global, NUinit est défini par å, Umaxj, et NVimtest défini par ¾ Vmaxj, .

[0093] La grille d’indexation globale est ensuite remplie par les sous-grilles, dans un ordre décroissant de taille de sous-grilles, où la taille de chaque sous-grille est déterminée par la taille Umax, Vmax du cadre englobant chaque sous-grille.

[0094] Pour chaque sous-grille, la sous-grille est insérée dans la grille d’indexation globale en déterminant une translation associant à chaque indice occupé de la sous-grille un indice correspondant dans la grille d’indexation globale, de manière à respecter les deux conditions suivantes :

Chaque sous-grille insérée dans la grille de paramétrage global ne doit occuper que des indices non occupés de cette grille,

Chaque sous-grille doit être insérée de manière à être la plus proche possible des autres sous-grilles, de manière à minimiser le nombre d’indices inoccupés de la grille d’indexation globale.

[0095] Une fois que toutes les sous-grilles ont été insérées, la grille d’indexation globale est redimensionnée de façon à correspondre à un cadre en deux dimensions minimal englobant toutes les sous-grilles. La taille de ce cadre est dans la plupart des cas très inférieure à la taille initiale de l’espace indexé de N Uinit x NVimt. [0096] Un exemple de résultat obtenu est représenté à la figure 5b. Une fois ce remaillage obtenu en deux dimensions, il peut être aisément propagé selon la coordonnée W pour obtenir un maillage hexaédrique complètement structuré.

[0097] Comme indiqué ci-avant, et comme visible sur la figure 5b, certaines mailles géométriquement voisines sont indexées avec des indices non consécutifs dans le nouveau maillage. Par exemple, deux mailles positionnées respectivement à (U,V) et (U+1,V) dans le maillage non structuré initial peuvent se trouver avec des indices (U,V) et (U+10, V-30) dans le nouveau maillage. Toutes les paires de mailles initialement voisines (au moins l’un des indices U et V est consécutif entre les deux mailles dans le maillage initial) le long d’une faille ou d’un plan de découpe (ou une ligne de découpe, en deux dimensions) et qui ne sont plus indexées avec des indices consécutifs à l’issue de l’étape 5b sont identifiées lors d’une étape 830. Sur la figure 5b, on a ainsi repéré par des flèches et des traits de figurés identiques les surfaces de paires de mailles initialement voisines.

[0098] A l’issue de ce traitement, on a donc obtenu un maillage hexaédrique structuré représentant le modèle géométrique initial. Ce maillage est donc compatible avec tous les simulateurs d’écoulement actuellement exploités.

[0099] Dans le cas décrit en référence à la figure 6a où le modèle géométrique initial comprenait des failles en X ou en Y ou des puits ou segments de puits sensiblement horizontaux formés par l’intersection de surfaces virtuelles de ce type, dont certaines failles ont été temporairement supprimées à l’étape préliminaire 90, le procédé comprend en outre les étapes supplémentaires 900 permettant d’introduire à nouveau les failles supprimées dans le maillage.

[0100] Tout d’abord, les cellules du maillage hexaédrique structuré obtenu à l’issue de la mise en oeuvre du procédé décrit ci-avant qui sont situées à l’intersection avec l’emplacement initial de la faille supprimée ou de chaque faille supprimée, sont supprimées 910 de manière à former une ou plusieurs cavités dans le maillage.

[0101] Ensuite, chaque faille supprimée est réintroduite dans le maillage hexaédrique, et les vides restant du maillage, à l’intersection avec chaque faille réintroduite, sont remplis 920 par des cellules polyédriques, qui peuvent par exemple être des cellules tétraédriques ou hexaédriques ou pyramides ou prismes ou polyédrique arbitraire ou un mélange des différents éléments. Ce remplissage peut par exemple être mis en oeuvre par tessellation de Delaunay contrainte, en utilisant les failles réintroduites et les bords de chaque cavité comme contraintes géométriques. On a représenté sur la figure 6b un exemple de maillage obtenu au voisinage d’une faille qui a été supprimée puis réintroduite conformément à ce mode de réalisation.

[0102] De ce fait, même pour des modèles géométriques complexes intégrant des failles en Y ou en X, ou pour représenter des puits complexes, on obtient un maillage qui est majoritairement hexaédrique structuré, avec uniquement quelques zones locales de maillage non structuré autour de certaines failles ou puits complexes. [0103] La mise en œuvre de la méthode décrite ci-avant permet donc d’obtenir une représentation numérique d’un réservoir ou d’un bassin géologique. La représentation numérique obtenue permet la mise en œuvre de la modélisation des propriétés pétrophysiques telles que la perméabilité et la porosité via des méthodes de géostatistiques. Des simulations d’écoulements de fluides peuvent être ensuite réalisées sur cette représentation numérique remplie de propriétés pétrophysiques dans des logiciels de simulation d’écoulements de fluides, par exemple pour, dans les cas de l’exploitation d’hydrocarbures, déterminer le positionnement des puits d’injection et de production ainsi que déterminer une pression optimale d’injection d’un fluide au niveau des puits afin d’obtenir un meilleur rendement de production. Dans les cas du stockage géologique du C02 ou H2 ainsi que l’exploitation de l’énergie géothermique, une simulation de géomécanique peut être menée dans le but de déterminer la déformation du réservoir en conséquence du changement de pression, de perméabilité, de porosité et d’autres propriétés engendrées durant la production/opération, et donc de prédire un risque de subsidence, une réactivation des failles du réservoir ainsi qu’une sismicité induite.