FR2952210A1 - Methode pour la construction d'un espace de depot correspondant a un domaine geologique - Google Patents

Abstract

La présente invention concerne une méthode pour la construction d'un espace de dépôt correspondant à un domaine géologique, comprenant les étapes de (i) partitionnement du domaine géologique actuel avec au moins un maillage conforme correspondant sensiblement aux limites dudit domaine géologique, (ii) calcul des coordonnées de dépôt définissant un espace de dépôt, dans lequel lesdites coordonnées de dépôt comprennent des calculs de champs de déplacements comprenant des simulations de déformations mécaniques dans le domaine géologique, utilisant un modèle de déformation solide. La présente invention concerne également un programme informatique mettant en œuvre la méthode.

Description

-1- « Méthode pour la construction d'un espace de dépôt correspondant à un domaine géologique »

Domaine de l'invention

Le domaine général de l'invention concerne, sans y être limité, la modélisation de terrains stratifiés de subsurface, notamment dans le domaine des géosciences pétrolières. L'invention décrit une méthode de construction d'un espace de dépôt virtuel et fiable correspondant à l'environnement chronostratigraphique au moment du dépôt des couches de terrains géologiques. Cet espace de dépôt est alors utilisé comme espace de calcul dans lequel de nombreuses applications peuvent être exécutées avantageusement. De telles applications comprennent la modélisation et l'interpolation précises des propriétés des terrains géologiques comme la porosité ou la perméabilité, le contrôle de la qualité de l'interprétation des données de réflexion sismique, ou la construction de grilles satisfaisant la géométrie complexe de structures géologiques.

Contexte de l'invention

De nombreuses applications peuvent potentiellement bénéficier de l'utilisation d'un espace de calcul modélisant l'environnement chronostratigraphique au moment du dépôt des terrains (appelé « espace de dépôt »). L'exemple précis de la modélisation des propriétés physiques de terrains géologiques est avancé ci-dessous comme passage en revue de l'état antérieur de la technique pour donner une vision claire du contexte de l'invention. Les propriétés physiques d'un domaine géologique tridimensionnel telles que la porosité ou la perméabilité sont généralement modélisées par des méthodes géostatistiques. Ces méthodes interpolent les propriétés du sous-sol dans un maillage tridimensionnel haute résolution en fonction de données généralement ponctuelles, telles que les mesures obtenues dans des puits. Ce processus d'interpolation utilise de façon intensive les distances, telles que les distances euclidiennes entre les centres des éléments du maillage, ou les distances de corrélation fournies par des variogrammes.

Les distributions calculées des propriétés physiques doivent refléter le paléo-environnement au moment du dépôt des terrains ; ainsi l'interpolation des données obtenues dans les puits n'est juste que si les distances euclidiennes calculées sont proches des distances équivalentes au moment du dépôt (appelées « distances géodésiques »). Toutefois, les roches ont généralement subi des altérations depuis leur dépôt, soit par érosion, soit par des failles ou des plis provoqués par des contraintes tectoniques. Par conséquent, la géométrie actuelle d'un domaine géologique est généralement très différente de sa géométrie au moment du dépôt. Ainsi, les méthodes géostatistiques produisent des résultats inexacts si elles sont appliquées directement dans l'espace actuel, décrit par le système de coordonnées cartésiennes (x, y, z). Ce problème peut être compensé en appliquant les méthodes géostatistiques de préférence dans un « espace de calcul » ou dans un « espace de dépôt » qui vise à représenter l'environnement chronostratigraphique au moment du dépôt. Cet espace de calcul ou de dépôt est généralement défini par un système de coordonnées curvilinéaires (u, v, w) également appelé « système de coordonnées de calcul » ou « système de coordonnées de dépôt ». Par la suite, une « isochrone » fera toujours référence à une surface joignant des points dans l'espace actuel où les sédiments se sont déposés au même moment. Les horizons sismiques ou le toit des couches géologiques qui ne représentent pas de lacunes dans un profil stratigraphique sont généralement isochrones. Les limites de séquences correspondant à des discordances stratigraphiques ne sont pas isochrones. On parle de « lacunes » lorsque les points de l'espace de dépôt n'ont pas de correspondance dans l'espace actuel, et de « chevauchements » lorsque les points de l'espace de dépôt correspondent à plus d'un emplacement dans l'espace actuel. Les espaces de dépôt, tels que définis dans des documents antérieurs de l'état de la technique, doivent respecter les deux propriétés suivantes : - Propriété P1 : Dans l'espace de dépôt, la géométrie du domaine géologique représente un espace chronostratigraphique dans lequel les isochrones identifiées au sein des séquences stratigraphiques sont essentiellement planes et parallèles, et dans lequel les déformations contemporaines et postérieures aux dépôts (c.-à-d. failles et plis) ont été effacées pour la plupart ; - Propriété P2 : Tout point (u, v, w) de l'espace de dépôt situé au sein d'une séquence stratigraphique n'a qu'un seul et unique emplacement (x, y, z) correspondant dans l'espace actuel.

Autrement dit, la portion de l'espace de dépôt qui correspond à une séquence stratigraphique ne contient ni lacune ni chevauchement, permettant la mise en correspondance des points d'un espace à un autre sans aucune ambiguïté. Jusqu'à présent, deux types de solutions principales existent dans l'état de la technique pour définir un tel espace de dépôt. On appelle le premier type de solutions « solutions d'indexation (i, j, k) »car elles consistent à construire un maillage structuré conforme tridimensionnel dans l'espace actuel en utilisant un système de coordonnées cartésiennes (x, y, z) qui représente la géométrie actuelle du domaine géologique. Ce maillage étant structuré, un indice (i, j, k) peut donc être attribué à chaque noeud et élément du maillage. Cette indexation est telle que les voisins d'un noeud ou d'un élément peuvent être obtenus par simple transformation de l'indice, ces transformations étant identiques pour tous les noeuds et éléments du maillage. Des exemples de telles transformations sont (i-1, j, k), (i, j-1, k), (i, j, k-1), (i+1, j, k), (i, j+1, k) et (i, j, k+1). Une fois l'indexation (i, j, k) terminée, l'espace de dépôt est très simplement défini par un système de coordonnées curvilinéaire (u, v, w) tel que u(x, y, z) = i(x, y, z), v(x, y, z) = j(x, y, z) et w(x, y, z) = k(x, y, z). Le maillage étant conforme, il existe une série de facettes dans le maillage correspondant à chaque horizon. Par conséquent, si tous les noeuds de ces facettes ont le même indice k, la propriété « P1 » est respectée. Autrement dit, le maillage doit être ce qu'on appelle souvent une « grille stratigraphique ». Cette solution a certains inconvénients et limites : - D'une part, la plupart des algorithmes servant à construire de telles grilles stratigraphiques créent des éléments très déformés (c.-à-d. étirés, écrasés, concaves ou de volume nul) dans l'espace (x, y, z) actuel, tout particulièrement à proximité des failles. Ces éléments étant tous des cuboïdes droits ou des parallélépipèdes réguliers dans l'espace de dépôt (u, v, w), cela signifie qu'une distorsion est introduite lors de la mise en correspondance des points d'un espace à l'autre. Autrement dit, les distances mesurées dans l'espace de dépôt devraient être mises à l'échelle de la distorsion afin de correspondre aux distances géodésiques. Cette distorsion peut être très élevée localement et elle diminue l'exactitude de la modélisation des propriétés. Par exemple, elle peut altérer radicalement le calcul du volume poreux des roches, et donc les quantités de pétrole et de gaz contenues dans ces roches. -4- - D'autre part, il existe certaines configurations habituelles pour lesquelles il n'est pas possible de garantir que deux éléments n'auront pas le même indice (i, j, k), ainsi des cellules « nulles » ou « mortes » sont ajoutées artificiellement aux grilles stratigraphiques pour résoudre ces configurations problématiques. Par conséquent, il y aura des lacunes ou des chevauchements dans l'espace de dépôt. Ces deux situations sont inacceptables car elles enfreignent la propriété « P2 » obligatoire. Un deuxième type de solutions, appelées solutions « paramétriques », consistent à calculer les coordonnées de dépôt par une interpolation de certaines quantités suivant une série de règles géométriques. Des exemples d'une telle solution géométrique sont fournis dans les documents WO 03/050766 par Deny et al., et WO 2005/119304 par Dulac et al., lesquels sont tous deux incorporés dans le présent document par référence. Trois fonctions de transfert u(x, y, z), v(x, y, z) et w(x, y, z) sont interpolées aux noeuds d'un maillage conforme 3D représentant la géométrie actuelle du domaine géologique. La fonction de transfert w(x, y, z) est calculée en premier. Afin de satisfaire aux exigences décrites ci-dessus, l'interpolation de la fonction de transfert w(x, y, z) est effectuée de sorte que des iso-valeurs (isosurfaces) particulières de la fonction se rapprochent de la géométrie des horizons, ce qui garantit que la propriété "Pl" est respectée. Puis les deux autres fonctions de transfert u(x, y, z) et v(x, y, z) sont calculées. L'interpolation de ces fonctions de transfert est effectuée de telle sorte que leurs gradients soient orthogonaux à la fois l'un à l'autre et au gradient de la fonction de transfert w(x, y, z) calculée précédemment, et de sorte que les longueurs de leurs gradients soient égales. Par ailleurs, comme condition aux limites supplémentaires, les valeurs de ces deux fonctions de transfert sont calculées selon un horizon de référence, précisant que leurs gradients sont orthogonaux l'un à l'autre et de même longueur. Cet horizon de référence est choisi de sorte à ce qu'il coupe un nombre maximum de failles. Le document WO 2005/119304 donne l'exécution de l'interpolation sur un maillage tétraédrique (maillage dont les éléments sont tous des tétraèdres). Cette solution a certains inconvénients et limites : - En premier lieu, elle tombe dans la catégorie des approches géométriques et cinématiques de modélisation des terrains au cours du temps. Comme il est mentionné dans l'article intitulé « Space-time Mathematical Framework for Sedimentary Geology », de Mallet, Mathematical Geology, Vol. 36, No. 1, 2004 incorporé dans le présent document par référence, le calcul d'un système de coordonnées curvilinéaires (u, v, w) ayant recours aux contraintes sur les gradients telles que -5- décrites ci-dessus revient à supposer que les terrains ont été déformés par un « glissement dans le plan des couches » de style cinématique, c'est-à-dire dans lequel la déformation est facilitée par un glissement le long d'un nombre infini d'interfaces. Par conséquent, l'espace de dépôt résultant ne convient pas à la modélisation des propriétés physiques du domaine géologique si les terrains sont supposés s'être déformés selon un style cinématique différent ou selon toute autre combinaison de styles cinématiques, tels que des cisaillements verticaux, obliques, des écoulements parallèles aux failles, etc. Par ailleurs, il n'est pas clair dans ce cas quel style cinématique ou quelle combinaison de styles cinématiques utiliser pour modéliser la déformation des terrains. En effet, les méthodes cinématiques ne sont pas fondées sur les principes fondamentaux de la conservation de la masse et de la quantité de mouvement, lesquels régissent la déformation des roches. Par ailleurs, seule la déformation, qui dépend fortement du style cinématique ou de la combinaison des styles cinématiques utilisés, peut être calculée en utilisant des méthodes cinématiques, et quoique Mallet déclare dans « Space-time Mathematical Framework for Sedimentary Geology » qu'il est possible d'expliquer le développement des fractures et leur interaction, la mécanique des fractures nous apprend plutôt que c'est l'état des contraintes qui est une donnée requise. En second lieu, afin d'assurer que les points de l'espace de dépôt n'aient qu'une seule correspondance dans l'espace actuel afin d'éviter les lacunes et les chevauchements et donc d'enfreindre à la propriété « P2 », des conditions aux limites précises doivent être appliquées aux failles lorsqu'on a recours à une telle solution paramétrique : ces solutions consistent à imposer des directions et des magnitudes de glissement arbitraires à tous les noeuds des facettes du maillage correspondant aux failles. Ceci peut être considéré comme une limitation sévère car elle impose des directions et des magnitudes de glissement qui ne sont pas connues à priori et qui peuvent ne pas correspondre aux rejets des failles qui résultent réellement de la déformation mécanique des roches. Comme il n'est pas certain quels vecteurs de glissement devraient être utilisés comme conditions aux limites, de tels choix arbitraires peuvent conduire à une distorsion notable lors de la mise en correspondance des points d'un espace à un autre, et par conséquent, la géométrie du domaine géologique dans l'espace de dépôt peut ne pas être proche de sa géométrie au moment du dépôt de celui-ci. Autrement dit, les distances mesurées dans l'espace de dépôt devraient être mis à l'échelle de la distorsion pour correspondre aux distances géodésiques. Cette distorsion peut être très élevée localement et elle diminue l'exactitude de la modélisation des propriétés. Par exemple, elle peut -6- modifier de façon significative le calcul du volume poreux des roches, et donc les quantités de pétrole et de gaz contenues dans ces roches. L'invention a pour objectif de proposer une méthode pour calculer physiquement un système de coordonnées curvilinéaires (u, v, w) définissant un espace de dépôt.

Résumé de l'invention

Afin d'atteindre cet objectif, l'invention propose une méthode pour la construction d'un espace de dépôt correspondant à un domaine géologique, comprenant les étapes suivantes : - partition dudit domaine géologique par au moins un maillage conforme correspondant sensiblement aux limites dudit domaine géologique, - calcul des coordonnées de dépôt définissant un espace de dépôt, dans lequel le calcul desdites coordonnées de dépôt comprend le calcul des champs de déplacement, faisant appel à des simulations de déformations mécaniques dans le domaine géologique, au moyen d'un modèle de déformation solide. L'espace de dépôt tel que calculé par la méthode selon l'invention peut être défini avantageusement (sans vouloir être limitatif) de sorte à respecter les deux propriétés suivantes (appelées propriétés « P1 » et « P2 ») : - Propriété P1 : Dans l'espace de dépôt, la géométrie du domaine géologiquereprésente un espace chronostratigraphique dans lequel les isochrones identifiées au sein des séquences stratigraphiques sont essentiellement planes et parallèles, et dans lequel les déformations contemporaines et postérieures aux dépôts (c.-à-d. failles et plis) ont été effacées pour la plupart ; - Propriété P2 : Tout point de l'espace de dépôt situé au sein d'une séquence stratigraphique n'a qu'un seul emplacement correspondant dans l'espace actuel. Autrement dit, la portion de l'espace de dépôt correspondant à une séquence stratigraphique ne contient ni lacune ni chevauchement, permettant la mise en correspondance des points d'un espace à un autre sans aucune ambiguïté. Ainsi, - ledit maillage conforme étant défini par les coordonnées (x, y, z) des noeuds des éléments, - les coordonnées de dépôt (u, v, w) peuvent être définies comme étant u = u(x, y, z), v = v(x, y, z), et w= w(x, y, z), et, -7- - les champs de déplacement (dx, dy, dz) peuvent être définis de telle sorte que : u(x, y, z) = x + dx(x, y, z), v(x, y, z) = y + dy(x, y, z), w(x, y, z) = z + dz(x, y, z). Le comportement mécanique du solide peut être quelconque, élastique, plastique, visqueux, ou une combinaison quelconque de ces comportements. Le calcul des coordonnées de dépôt (u, v, w) peut comprendre : - une étape de mise en correspondance des limites, comprenant une opération associant chacune des limites du domaine géologique à une série de facettes ou d'éléments du maillage correspondant sensiblement à la géométrie de ladite limite, la série comprenant au moins une facette d'un maillage conforme 3D, ou un élément d'un maillage conforme 2D (défini dans l'espace 3D), - une étape de préparation de la simulation mécanique, comprenant la définition des conditions aux limites visant à contraindre la géométrie calculée du domaine géologique dans l'espace de dépôt (de sorte que les deux propriétés « P1 » et « P2 » soient toutes les deux respectées, sans intention de limitation), - une étape de simulation mécanique, comprenant le calcul dans le domaine géologique (actuel) (x, y, z) d'un champ de déplacements (dx, dy, dz) correspondant auxdites conditions aux limites. Ainsi l'invention permet de calculer un système de coordonnées de dépôt (u, v, w), respectant les propriétés « P1 » et « P2 », dans lequel la géométrie du domaine géologique tridimensionnel est physiquement fiable, permettant ainsi la définition d'un espace de calcul à faible distorsion où les distances euclidiennes mesurées sont proches des distances géodésiques. Lors de leur application dans cet espace de dépôt, les interpolations géostatistiques ou les algorithmes de simulation peuvent modéliser les propriétés physiques du domaine géologique tridimensionnel avec une précision accrue. La génération de maillages respectant les géométries complexes des structures géologiques est également simplifiée lorsqu'effectuée dans cet espace de dépôt. La qualité de l'interprétation de la qualité d'une sismique réflexion 3D peut être plus précisément vérifiée lorsque celle-ci est effectuée dans cet espace de dépôt. Dans tout ce qui suit, le terme « horizon » est utilisé au sens très large, faisant référence à une isochrone précise identifiée au sein d'une séquence du domaine géologique, sans représenter de -8- lacune dans la série stratigraphique. Dans le système de coordonnées de dépôt (u, v, w) calculé par la méthode selon l'invention, les horizons peuvent être plans et parallèles (propriété « P1 »). Par conséquent les couches ne peuvent pas être faillées dans l'espace de dépôt. Par comparaison à un espace de dépôt obtenu par une méthode paramétrique selon un état technique antérieur tel que "Space-time Mathematical Framework for Sedimentary Geology" de Mallet, il existe une différence fondamentale concernant un espace de dépôt calculé par l'invention décrite dans le présent document : - aucune supposition n'est requise en ce qui concerne le style cinématique ou la combinaison de styles cinématiques utilisés (glissement dans le plan des couches, cisaillement vertical, oblique, écoulement parallèle aux failles, etc.), - aucune supposition n'est requise pour les directions de glissement et les magnitudes le long des failles. Avec la présente invention, les failles (ou toute autre limite) sont libres de glisser et de se déformer dans la mesure où le contact mécanique est respecté (propriété « P2 »). Cela signifie que l'équilibre des forces a été atteint, sans chevauchement ou pénétration à travers les failles. Les directions correspondantes de glissement et les magnitudes sont donc des résultats du calcul plutôt que des paramètres d'entrée. Les vecteurs de glissement produits résultent d'une simulation qui respecte les principes fondamentaux à la fois de la mécanique de fracturation et de la physique qui régissent la déformation des roches. Par conséquent, ils produisent une géométrie physique fiable du domaine géologique dans l'espace de dépôt. La même remarque s'applique aux horizons : Avec la méthode paramétrique telle que celle dont il est question dans le document WO 2005/119304, il n'est pas possible de modéliser le glissement/détachement de la couche (ou toute autre forme de glissement le long d'un horizon), qui aurait pu affecter certaines couches du domaine géologique, alors qu'il est très facile de le faire avec la méthode selon l'invention, en définissant simplement une condition aux limites qui permet le contact mécanique avec cisaillement aux interfaces des couches. Par ailleurs, la présente invention ne simplifie par le modèle géologique traité. Toutes les failles sont prises en compte, quelle que soit le type de connexion entre elles. Les failles peuvent avoir des bords libres et n'ont pas besoin d'être prolongées artificiellement pour être complètement connectées à d'autres limites du modèle. La présente invention n'est pas non plus limitée à un contexte tectonique particulier. Par exemple, les structures en extension ou en compression sont -9- gérées de façon identique. Ainsi la présente invention peut gérer toute les complexités des modèles géologiques. Afin de respecter la propriété « PI », l'étape de préparation de la simulation mécanique peut comprendre par ailleurs la définition d'une condition aux limites imposant aux noeuds correspondant à un horizon donné du domaine géologique à avoir sensiblement la même valeur finale w(x, y, z), - la valeur finale w(x, y, z) correspondant à un horizon donné peut être un résultat de l'étape de simulation mécanique, - la valeur finale w(x, y, z) correspondant à un horizon donné peut être précisée explicitement, - la condition aux limites peut également imposer que l'angle entre les paires de facettes adjacentes ou éléments correspondant à l'horizon soit sensiblement égal à un arc de 180 degrés. Afin de respecter la propriété « P2 », l'étape de préparation de la simulation mécanique peut également comprendre la définition de : - une condition aux limites exigeant que les noeuds correspondant à un côté donné d'une faille donnée du domaine géologique, soient libres de mouvement à condition qu'ils restent en contact avec l'autre côté de ladite faille, lequel contact sera respecté mécaniquement, - une condition aux limites exigeant que les noeuds situés à l'intersection entre un horizon donné du domaine géologique et un côté donné d'une faille donnée du domaine géologique, soient libres de mouvement à condition qu'ils restent en contact avec les bords correspondant à l'intersection entre ledit horizon et l'autre côté de ladite faille, le contact étant respecté mécaniquement, - une condition aux limites exigeant que les noeuds correspondant à un côté donné d'un horizon donné du domaine géologique, soient libres de mouvement à condition qu'ils restent en contact avec l'autre côté dudit horizon, lequel contact sera respecté mécaniquement. Ainsi, - les déplacements (dx, dy, dz) peuvent être calculés pour les noeuds du (des) maillage(s) représentant la géométrie (actuelle) du domaine géologique, - l'étape de simulation mécanique n'est pas limitée à l'hypothèse de faibles déplacements et peut gérer les grandes déformations, évitant dans ce cas des effets secondaires indésirables tels que des rotations, - les séquences identifiées dans le domaine géologique peuvent être traitées séparément, 2,952210 -10- - l'étape de simulation mécanique peut également comprendre l'utilisation d'une formulation « maître/esclave »pour contraindre les déplacements (dx, dy, dz) des noeuds à rester mécaniquement en contact avec les séries d'éléments, de facettes ou de bords tels que définis dans les conditions aux limites, - la formulation maître/esclave peut être basée sur une méthode de pénalisation, - la formulation maître/esclave peut être basée sur l'utilisation de multiplicateurs de Lagrange. Dans une première version, appelée la version « FEM », l'invention fournit une méthode dans laquelle le champ des déplacements (dx, dy, dz) est calculé selon la méthode des éléments finis (Finite Element Method - FEM), la géométrie (actuelle) du domaine géologique étant représentée par au moins un maillage conforme 3D, les éléments dudit maillage étant des polyèdres. Voici les aspects spécifiques de cette version : - les éléments des maillages conformes 3D peuvent inclure l'un des éléments suivants : tétraèdres, cuboïdes, - les maillages conformes 3D peuvent être hybrides et comprendre des groupes d'éléments structurés et des groupes d'éléments non structurés - les groupes d'éléments structurés peuvent comprendre des cuboïdes et les groupes d'éléments non structurés peuvent comprendre des tétraèdres et des pyramides à base quadrilatérale, - les éléments des maillages conformes 3D peuvent être associés à des paramètres représentant, les propriétés mécaniques locales du matériau contenu dans la portion du domaine géologique comprise dans ledit élément, - le modèle de déformation solide peut être un modèle de déformation élastique, - les propriétés mécaniques locales des matériaux peuvent comprendre au moins deux des paramètres suivants : coefficient de Poisson, module de Young, coefficients de Lamé - l'étape de mise en correspondance des limites peut également comprendre l'étape de récupération des éléments du maillage contenus dans une couche donnée du domaine géologique qui contient un élément donné, en explorant de façon récursive les éléments adjacents à un élément, et en sélectionnant ceux qui satisfont à la condition que les facettes partagées par des paires d'éléments voisins ne correspondent pas à une limite de cette couche, -11- - les éléments du maillage correspondant à une couche donnée du domaine géologique peuvent être associés à des propriétés mécaniques locales qui sont sensiblement identiques de sorte à définir une couche de matériau sensiblement homogène, - les propriétés mécaniques des matériaux associés aux éléments ou ensembles d'éléments du maillage conforme 3D peuvent varier en fonction du temps et/ou des itérations de calcul, - les déplacements (dx, dy, dz) peuvent être calculés en construisant une matrice de rigidité globale, comme pour une méthode classique basée éléments finis connue de la technique, - les déplacements (dx, dy, dz) peuvent être calculés de façon itérative en construisant et en mettant à jour des matrices de rigidité locales définies aux noeuds du maillage 3D représentant la géométrie (actuelle) du domaine géologique, permettant aux forces d'être transmises de noeud à noeud (ou d'une série de noeuds à d'autres séries de noeuds), jusqu'à ce que l'équilibre soit atteint, selon n'importe quelle méthode de « relaxation dynamique », « multi-agent », de Gauss-Seidel ou de Jacobi connue de la technique, La méthode selon la version « FEM » peut comprendre par ailleurs le calcul à unemplacement (x, y, z) dans le domaine géologique de la valeur d'une propriété de terrain calculée (par exemple simulée selon une méthode géostatistique effectuée dans l'espace de dépôt) à un emplacement (u, v, w) dans l'espace de dépôt, le calcul de ladite propriété de terrain comprenant les étapes suivantes : - récupération de l'élément de maillage contenant l'emplacement (x, y, z), - calcul des coordonnées de dépôt (u, v, w) dudit emplacement (x, y, z) au moyen d'une interpolation analytique basée sur les valeurs des coordonnées de dépôt (u, v, w) correspondant aux coordonnées des noeuds (x, y, z) de l'élément de maillage, - récupération de la valeur de la propriété de terrain correspondant aux coordonnées (u, v, w) dans l'espace de dépôt. Selon certains aspects bien précis : - l'élément contenant peut être un tétraèdre et le mécanisme d'interpolation analytique peut être basé sur l'utilisation de coordonnées barycentriques, - l'élément contenant peut être un cuboïde et le mécanisme d'interpolation analytique peut être une interpolation trilinéaire. Dans une seconde version, appelée la version « BEM », l'invention fournit une méthode dans laquelle le champ des déplacements (dx, dy, dz) est calculé en utilisant la méthode des éléments -12- frontière (Boundary Element Method - BEM), la géométrie (actuelle) du domaine géologique étant représentée par une série de maillages conformes bidimensionnels incorporés en trois dimensions, chacun de ces maillages représentant une limite dudit domaine géologique, les éléments de ces maillages bidimensionnels étant des polygones. Voici les aspects spécifiques de cette version : - les éléments des maillages 2D peuvent comprendre l'un des éléments suivants : triangles, quadrilatères, - certains maillages peuvent être hybrides et comprendre des groupes d'éléments structurés et des groupes d'éléments non structurés - les groupes d'éléments structurés peuvent comprendre des quadrilatères et les groupes d'éléments non structurés peuvent comprendre des polygones, - les groupes d'éléments non structurés peuvent comprendre des triangles, - chaque couche du domaine géologique peut être associée à des paramètres représentant les propriétés mécaniques du matériau contenu dans la portion dudit domaine géologique correspondant à l'intérieur de ladite couche, définissant ainsi un matériau homogène pour cette couche, - le modèle de déformation solide peut être un modèle de déformation élastique, - les propriétés mécaniques des matériaux peuvent comprendre au moins deux des paramètres suivants : coefficient de Poisson, module de Young, coefficients de Lamé - les propriétés mécaniques associées à une couche peuvent varier en fonction du temps et/ou des itérations de calcul, - l'étape de préparation de la simulation mécanique peut également comprendre le calcul de déplacements (dx, dy, dz) à des emplacements (x, y, z) arbitraires appelés « points d'observation » en utilisant les déplacements (dx, dy, dz) calculés aux emplacements des noeuds des maillages bidimensionnels représentant la géométrie du domaine géologique tel que décrit dans toute méthode des éléments frontière connue, La méthode d'analyse par éléments frontière (« BEM ») peut comprendre par ailleurs le calcul de la valeur à unemplacement (x, y, z) d'une propriété de terrain calculée (par exemple simulée selon une méthode géostatistique effectuée dans l'espace de dépôt) à un emplacement (u, v, w) dans l'espace de dépôt, le calcul de ladite propriété de terrain comprenant les étapes suivantes : - calcul des coordonnées de dépôt (u, v, w) dudit emplacement (x, y, z) utilisé comme point d'observation, -13- - récupération de la valeur de la propriété de terrain correspondant aux coordonnées (u, v, w) dans l'espace de dépôt. L'invention a pour objectif supplémentaire de fournir un programme mettant en oeuvre la méthode.

Description des dessins

Les méthodes selon les différentes versions de la présente invention peuvent être mieux comprises en se référant aux dessins, sachant que ces dessins sont fournis à des fins d'illustration uniquement et ne sont pas limitatifs. D'autres aspects, objectifs et avantages de l'invention seront mis en évidence à partir des descriptions fournies ci-dessous. - La figure 1 montre une vue en coupe transversale verticale d'un domaine géologique en 3D composé de quatre séquences, avec horizons, failles et discordances ; - la figure 2 montre une vue en coupe transversale verticale d'un maillage conforme non structuré incorporé au domaine géologique de la figure 1 ; - la figure 3 montre une vue en coupe transversale verticale des espaces de dépôt associés au domaine géologique des figures 1 et 2 (un seul espace de dépôt par séquence), tel que calculé selon la méthode des éléments finis ; - la figure 4 illustre les différences entres les distances euclidiennes et les distances géodésiques lorsque mesurées (a) dans le domaine géologique actuel et (b) dans l'espace de dépôt ; - La figure 5 illustre une séquence d'opérations selon l'invention : (a) dans le domaine géologique actuel dont les limites ont été identifiées, (b), discrétisation des limites et/ou des couches avec maillages 2D ou 3D conformes, et (c) calcul de l'espace de dépôt. (d) montre le détail de la couche figurée en (c) dans l'espace de dépôt dans les environs d'une faille, où un certain glissement horizontal peut être observé suite à la simulation mécanique ; - La figure 6 illustre une méthode de post-traitement, montrant en (a) le maillage haute résolution en 3D dans l'espace de dépôt peuplé de propriétés des matériaux simulées dans cet espace et en (b) une vue interne d'un autre maillage représentant une couche du domaine géologique actuel peuplée des propriétés simulées dans le maillage figuré en (a). -14- - La figure 7 illustre une méthode de post-traitement pour la construction d'une grille en 3D hybride (essentiellement structurée) représentant la géométrie actuelle du domaine géologique, peuplé de propriétés simulées dans l'espace de dépôt produit par la méthode selon l'invention : sont figurés en (a) le domaine géologique en 3D avec les couches et les limites identifiées, en (b) le même domaine géologique maillé avec la grille ainsi créée, peuplé avec les propriétés simulées dans l'espace de dépôt, et en (c) une vue éclatée des séquences du domaine géologique figurées en (b).

Description détaillée de l'invention

En se référant aux figures 1 et 2, sans perte de généralité et à moins d'avoir été précisé autrement, les définitions suivantes relatives à la modélisation de structures géologiques s'appliquent à l'ensemble du présent document : - « Maillages » 9 se réfèrent aux partitions d'un domaine en des sous-ensembles géométriques plus petits appelés « éléments » 11. Ainsi, les maillages 9 sont des modèles discrets du domaine qu'ils représentent. Les éléments de maillages tridimensionnels (3D) sont des volumes tels que des tétraèdres ou des cuboïdes. Les éléments 12 de maillages bidimensionnels (2D) sont plans, tels que des triangles, ou des quadrilatères quoiqu'ils aient aussi une géométrie 3D ; - Chacun des éléments simples 11 d'un maillage 9 relie un nombre variable de points arbitraires (au moins quatre en 3D, trois en 2D) du domaine 3D appelés « noeuds » 10 - un maillage 9 est dit « structuré » Si de telles connexions entre les noeuds 10 suivent le même format élémentaire pour tous les noeuds du maillage. Sinon, un maillage est dit « non structuré », et dans ce cas les noeuds ont un nombre variable de connexions aux autres noeuds ; - les « faces » d'un élément font référence aux intersections non-vides de cet élément avec tout plan dans l'espace 3D, à condition que l'élément soit entièrement compris dans l'un des deux demi-espaces définis par le plan. Les faces de dimension 0 s'appellent les noeuds de l'élément, l'ensemble des faces de dimension 1 s'appellent les « arêtes » de l'élément, et l'ensemble des faces de dimension 2 s'appellent les « facettes » 12 de l'élément (un élément de maillage 2D n'a pas de facettes, seulement des noeuds et des arêtes) ; - à l'échelle d'un bassin ou d'un réservoir, les domaines géologiques 3D 1 sont limités par des limites externes et internes appelées « limites » 4. Une limite est un objet 2D équivalent à une surface, incorporée dans l'espace 3D ; - un maillage 3D 9 représentant un domaine géologique 3D est dit « conforme » s'il existe un sous-ensemble de facettes 12 du maillage 3D correspondant (dans le sens le plus large) à une limite 4 donnée, un tel sous-ensemble de facettes formant une approximation de la géométrie de la limite - un maillage 2D 9 représentant un domaine géologique 3D est dit « conforme » s'il existe un sous-ensemble de facettes du maillage 2D correspondant (dans le sens le plus large) à une limite 4 donnée, un tel sous-ensemble de facettes formant une approximation de la géométrie de la limite ; - Les « failles » 7 se réfèrent aux limites d'un domaine géologique 3D correspondant à des discontinuités dans la roche ; - Les « horizons » 5 se réfèrent à des limites d'un domaine géologique 3D séparant une roche plus jeune d'une roche plus ancienne sans représenter de lacune dans le profil stratigraphique; - par opposition aux horizons 5, les limites 4 d'un domaine géologique 3D qui séparent une roche plus jeune d'une roche plus ancienne et qui représentent une lacune dans le profil stratigraphique, sont appelées des « discordances » 6. - Les « isochrones » se réfèrent à des surfaces reliant des points dans l'espace actuel où les sédiments se sont déposés au même moment. Les horizons 5 sont généralement des isochrones tandis que les discordances 6 ne le sont pas. Dans tout ce qui suit, le terme « horizon » se réfère à des isochrones particulières identifiées au sein d'un domaine géologique. - les régions de l'espace fermées et entièrement limitées par des limites 4 d'un domaine géologique 3D correspondant ainsi à des sous-volumes dudit domaine sont appelées « couches » 3 ; - Les « séquences » 8 se réfèrent à des ensembles de couches limitées par des discordances 6. Cette définition implique que les dépôts sédimentaires d'une séquence 8 se sont produits continuellement dans le temps, sans érosion ou lacune de sédimentation. Dans la description suivante, plusieurs versions de l'invention sont décrites et des exemples précis sont avancés afin de permettre une bonne compréhension de l'invention. Toutefois, il sera évident pour un expert en la matière que d'autres versions sont possibles et que l'invention n'est pas limitée aux exemples décrits dans le présent document. Par ailleurs, les caractéristiques qui sont bien -16- connues d'un expert en la matière peuvent être omises ou simplifiées pour éviter de brouiller la description des différentes versions de l'invention. En se référant aux figures 1, 2 et 3, la méthode selon l'invention est composées des étapes suivantes : - une étape de « mise en correspondance des limites » pour associer chaque limite 4 pertinente du domaine géologique 1 à un « identifiant géologique » unique. Cet identifiant géologique est également associé aux facettes 12 du maillage 3D conforme 9 ou aux éléments 12 d'un maillage 2D qui approchent la géométrie de la limite 4, permettant ainsi la récupération sans équivoque soit de toutes les facettes 12 d'un maillage conforme 3D 9 soit de tous les éléments 12 d'un maillage 2D qui correspondent à une limite 4 donnée du domaine géologique 1. Cette étape ne s'applique qu'aux limites 4 qui sont impliquées dans les étapes suivantes de la méthode ; - une « étape de préparation de simulation mécanique » pour la mise en place de conditions aux limites nécessaires pour assurer que la géométrie résultante du domaine géologique 3D 1 dans l'espace de dépôt 2 est telle que les deux propriétés « P1 » et « P2 » sont respectées (sans intention de limitation) ; - une « étape de simulation mécanique » pour calculer, partout dans l'espace actuel muni de coordonnées (x, y, z), du champ des déplacements (dx, dy, dz) correspondant aux conditions aux limites prescrites, de telle sorte que son calcul respecte les principes physiques qui régissent la déformation des roches. Les déplacements calculés permettent de définir un système de coordonnées de dépôt (u, v, w) tel que : u(x, y, z) = x + dx(x, y, z), v(x, y, z) = y + dy(x, y, z), et w(x, y, z) = z + dz(x, y, z). Le domaine géologique actuel 1 des figures 1 et 2 sont des vues transversales à y constant, ainsi seules les coordonnées x et z sont figurées. La représentation dans l'espace de dépôt de la figure 3 est une vue transversale à v constant, ainsi seules les coordonnées u et w sont figurées. Selon la figure 4, la méthode de l'invention permet de définir un espace de dépôt 2 muni d'un système de coordonnées de dépôt (u, v, w) respectant les deux propriétés « P1 » et « P2 », où la géométrie du domaine géologique 3D est fiable physiquement, permettant ainsi la définition d'un espace de calcul à faible distorsion où les distances euclidiennes 20 mesurées sont proches des 7 -- distances géodésiques 21. Les interpolations géostatistiques ou les algorithmes de simulation, lorsqu'appliqués à cet espace de dépôt 2, peuvent modéliser les propriétés physiques du domaine géologique 3D avec une meilleure précision. La génération de maillages respectant les géométries complexes des structures géologiques est également simplifiée lorsqu'effectuée dans cet espace de dépôt. Comme autre exemple, la qualité de l'interprétation d'une sismique réflexion 3D peut être plus précisément vérifiée lorsque celle-ci est effectuée dans cet espace de dépôt. Il existe de nombreuses méthodes de mises en oeuvre possibles potentielles pour l'étape de correspondance des limites selon les algorithmes et les structures de données utilisées. Généralement, plus la structure des données est complexe, plus l'algorithme sera simple, quoique la simplicité de l'algorithme sera modifiée par la gestion de mises à jour de la structure de données. Par ailleurs, une structure de données simple est facile à mettre à jour mais elle oblige l'algorithme à effectuer une série d'opérations bien précises pour recréer certaines données nécessaires à chaque fois. Par conséquent, le choix final représente souvent un compromis. Les identifiants géologiques associés aux limites du domaine géologique 1, B = {B1, B2, B3, ..., Bn}, où B est l'ensemble restreint des limites 4 utilisées dans les étapes suivantes de la méthode, peuvent être représentés par une série de n nombre entiers arbitraires {bi, b2, b3, ..., bn}. Quoique facultatif et pour des raisons de clarté, B peut être représenté par une liste linéaire de n enregistrements, c'est à dire un ensemble de n enregistrements dont les propriétés structurelles ne concernent essentiellement que la position relative des enregistrements. Les positions des enregistrements peuvent être mises en correspondance une à une avec l'ensemble des n premiers entiers, de sorte que l'identifiant géologique bi, associé au ie enregistrement de la liste B, puisse être égal à i. Selon la figure 5, dans une première version, le domaine géologique peut être représenté par un maillage 3D conforme 9 dont les éléments 11 sont tous des tétraèdres. La génération du maillage ne fait pas partie de l'invention et peut être effectuée selon toute méthode classique quelconque connue de la technique. Chaque tétraèdre est composé de quatre facettes 12, et chacune de ces quatre facettes peut être identifiée simplement par un indice entre 1 et 4. Les facettes dans le maillage tétraédrique peuvent ainsi être représentées sans équivoque par une paire (T, J) dans laquelle T est un tétraèdre et 3 un indice compris entre 1 et 4. Considérons également que la structure de données représentant -18- un tétraèdre t simple est codée de sorte que des cases mémoires sont réservées pour stocker et récupérer un tableau de quatre valeurs entières A(t) = {al, a2, a3, a4}, chaque aj étant associé à la facette (t, j), la je facette de t. les valeurs de aj peuvent être définies comme suit : - si la facette (t, j) ne correspond pas à une limite de l'ensemble B, alors aj = -1, - si la facette (t, j) correspond au ie enregistrement dans l'ensemble B, alors aj est égal à l'identifiant géologique correspondant bi (aj = i par exemple). La méthode suppose que le maillage tétraédrique est conforme, de sorte qu'il y ait au moins une facette qui corresponde à chacune des limites 4 contenues dans l'ensemble B. Les valeurs de aj peuvent être définies pendant une étape d'initialisation prenant place au début de l'étape de mise en correspondance des limites. Ainsi, la récupération de toutes les facettes du maillage tétraédrique qui correspondent à une limite donnée Bi du domaine géologique peut comprendre les étapes suivantes : - d'abord, création d'une liste linéaire vide de facettes Fi, - ensuite, récupération de l'identifiant géologique bi correspondant à Bi (égal à i par exemple), - enfin, exploration de tous les tétraèdres du maillage. Si, pour le tétraèdre donné t, il existe un j de telle sorte que dans A(t), aj = bi, alors la facette (t, j) est insérée dans la liste Fi. La liste des facettes Fi peut être divisée en deux listes Fi+ et Fi-. La liste Fi+ contient toutes les facettes de Fi dont le vecteur normal correspond localement (au sens le plus large) au vecteur normal de la limite Bi, le vecteur normal d'une facette (T, I) étant défini par exemple comme le vecteur normal au plan défini par les trois noeuds de la facette orientés vers l'extérieur de T. La liste Fi- contient toutes les autres facettes de Fi. Ces deux listes peuvent servir à récupérer toutes les facettes qui correspondent à un côté donné (« + » ou « - ») d'une limite donnée du domaine géologique. Dans cette première version, ou dès que la géométrie actuelle du domaine géologique 1 est représentée par un maillage conforme 3D 9 dont les éléments 11 sont des polyèdres, les déplacements inconnus (dx, dy, dz) peuvent être déterminés en utilisant la méthode des éléments finis (FEM). La méthode des éléments finis est une méthode numérique classique, bien connue dans la technique, elle ne sera donc que partiellement décrite ici afin de ne pas brouiller les différentes versions de l'invention. Une description détaillée de la méthode se trouve par exemple dans le livre de -19- Zienkiewicz et Taylor, "The Finite Element Method - Volume 1: The Basis", éditions Butterworth-Heinemann (2000). Dans cette version, appelée méthode FEM, chaque élément 11 du maillage conforme 3D 9 est associé à au moins deux paramètres représentant les propriétés mécaniques du matériau contenu dans la portion du domaine géologique correspondant à l'intérieur de cet élément. Le nombre de paramètres dépend du modèle de déformation solide mis en oeuvre. Un modèle de déformation élastique exige deux paramètres, tandis que des modèles de déformation plastique ou des déformations plus complexes demandent davantage de paramètres. Les valeurs des paramètres peuvent varier d'un élément 11 à un autre élément 11 à travers le maillage 9 et peuvent également changer en fonction du temps pendant l'étape de simulation mécanique. Si le matériau est supposé se déformer élastiquement pendant l'étape de simulation mécanique, ces paramètres sont généralement le coefficient de Poisson et le module de Young (ou les coefficients de Lamé correspondants) pour le matériau contenu dans l'élément du maillage. Afin de respecter la propriété « P1 » une condition aux limites est établie imposant que tous les noeuds 10 des facettes 12 correspondant à un horizon donné 5 du domaine géologique 1 aient la même valeur finale w(x, y, z) de sorte à apparaître dans un même plan dans l'espace de dépôt tel qu'illustré en figure 5(c). Cette valeur peut être un résultat de l'étape de simulation mécanique, ou elle peut être précisée explicitement. Afin de respecter la propriété « P2 », une autre condition aux limites est établie pour imposer que tous les noeuds 10 des facettes 12 correspondant à un côté donné d'une faille donnée 7 ou un horizon 5 du domaine géologique 1 soient libres de mouvement à condition qu'ils restent en contact avec les facettes correspondant à l'autre côté de la faille 7 ou de l'horizon 5. En option, une autre condition aux limites peut être précisée pour imposer que tous les noeuds 10 situés aux intersections entre un horizon 5 donné du domaine géologique 1 et un côté donné d'une faille donnée 7 du domaine géologique 1 soient libres de mouvement à condition qu'ils restent en contact avec les bords correspondant à l'intersection entre l'horizon 5 et l'autre côté de la faille 7. La figure 5(d) illustre un exemple de la façon dont une faille apparaît dans l'espace de dépôt, les deux côtés de celle-ci restant mécaniquement en contact sans interpénétration, mais libres de glisser. Le calcul des déplacements inconnus (dx, dy, dz) à l'emplacement de chaque noeud 10 du maillage 3D 9 représentant la géométrie actuelle du domaine géologique 1 est effectué par la méthode -20- des éléments finis basée sur les conditions aux limites définies dans l'étape précédente de préparation de la simulation mécanique. Chaque séquence 8 identifiée dans le domaine géologique 1 peut être traitée séparément, menant à une représentation de son propre espace de dépôt 2, comme l'illustre la figure 3. Le calcul par la méthode des éléments finis peut être effectué en construisant une matrice de rigidité globale comme dans toute méthode des éléments finis classique implicite quelconque connue dans la matière. - Une autre solution consiste à effectuer le calcul FEM de façon itérative en construisant et en mettant à jour des matrices de rigidité locales définies aux noeuds 10 des mailles 3D 9 représentant la géométrie actuelle du domaine géologique 1, permettant aux forces d'être transmises de noeud à noeud (ou d'une série de noeuds à d'autres séries de noeuds), jusqu'à ce que l'équilibre soit atteint, selon n'importe quelle méthode de « relaxation dynamique », « multi-agent », de Gauss-Seidel ou de Jacobi connue. De telles méthodes sont décrites par exemple dans le livre de Saad « Iterative Methods for Sparse Linear Systems", éditions Siam (2003). Le calcul FEM n'est pas limité à une hypothèse de petits déplacements et peut également gérer de grandes déformations, évitant dans ce cas les effets secondaires indésirables liés aux rotations. Afin de contraindre les noeuds 10 à rester mécaniquement en contact avec des ensembles de facettes 12 ou des bords comme cela peut être imposé par des conditions aux limites définies dans l'étape de préparation de la simulation mécanique, une formulation « maître/esclave » peut être utilisée. Cette formulation peut être basée par exemple sur une méthode de pénalisation classique ou sur l'utilisation des multiplicateurs de Lagrange. Un exemple de l'application de la méthode des éléments finis est donné ci-dessous pour le calcul de la déformation mécanique de solides en fonction d'un modèle de déformation élastique. Seul un aperçu est donné car la méthode des éléments finis est bien connue. Selon l'hypothèse de faible déplacement, le tenseur de déformation de Green-Lagrange E est défini par : -21- où dI,~ dénote la 3th dérivée partielle de la Ith composée du déplacement inconnu d = (dx, dy, dz). Dans le cas de solides élastiques isotropes, le tenseur des contraintes a est donné par la loi de Hooke : r -~ 2Ai 0 0 0 Ey>. À À À + 2p. o 0 0 Enz oxx 0 0 0 0 p 0° 2exa 0 0 0 0 L2sb Dans lequel les coefficients de Lamé A et p sont définis à partir du module de Young E et du coefficient de Poisson v: 1'i À _ (1+v)(1 - E 2(1+'v) L'équation de conservation de la quantité du mouvement dans un cas quasi-statique donne: div a(d) = 0, avec déplacements et conditions aux limites habituelles. En particulier, lorsque les noeuds 10 doivent rester mécaniquement en contact avec des ensembles de facettes 12 ou des bords, une contrainte de non-pénétration est établie, imposant que la distance de projection de tels noeuds sur de tels ensembles de facettes ou de bords au plan normal aux noeuds soit nulle (formulation « maître/esclave »). Selon la méthode des éléments finis, le champ de déplacements inconnus d est recherché dans le champ de déplacement cinématique admissible D. On définit également l'espace V de « test » le champ de déplacement admissible virtuel v. Le principe des travaux virtuels donne la loi de conservation de la quantité de mouvement sous la forme d'une intégrale : And d c D, W' L* E V, c(d.). v d = où S? est le solide considéré défini par ses coordonnées spatiales. Cette dernière équation peut être considérée comme un problème de minimisation d'une fonctionnelle. Ainsi, afin de tenir compte des contraintes de non pénétration, le problème de minimisation est réécrit, remplaçant la -22- fonctionnelle par une fonction de Lagrange L(d, A), où A est un multiplicateur de Lagrange inconnu. Appelons h le maillage utilisé pour la discrétisation du solide considéré. La méthode des éléments finis donne une valeur approximative aux inconnues dh et Ah dans des espaces dimensionnels finis, et un système vectoriel non linéaire R(dh, Ah) peut être obtenu en appliquant : = 0 `(t,t 1d 1 où { gi, ..., gn } et { Al, ..., An } sont la bases de fonctions de forme données approximant des déplacements et des multiplicateurs de Lagrange respectivement. Le système non-linéaire résultant peut être résolu itérativement afin d'annuler le résiduel de premier ordre à l'étape (n+l) à condition que R(dh, Ah) ne soit pas nul à l'étape n. Définissant : t'Adz-1 d"' û d ic, n.+1 n+s ùÀn A Il s'ensuit que : 'c R a(( ) >.d, où Gai(dA, 1.j a(d,À) est la matrice tangente du résiduel R. Cette matrice peut être inversée directement en deux dimensions avec une factorisation LU, ou résolue itérativement en trois dimensions en utilisant un algorithme d'Uzawa ou de facteur résiduel minimum (consulter par exemple le livre de Saad, "Iterative Methods for Sparse Linear Systems", Siam editions (2003)).

Enfin, les déplacements nodaux résolus (dx, dy, dz) sont utilisés pour définir le système de coordonnées de dépôt (u, v, w) à l'emplacement de chaque noeud 10 utilisé dans le maillage conforme en 3D 9. Selon les figures 2 et 3, un point P aux coordonnées cartésiennes (x, y, z) dans le domaine géologique 1 actuel reçoit des coordonnées (u, v, w) dans l'espace de dépôt 2. -23- Le calcul de la valeur à tout emplacement actuel (x, y, z) d'une propriété 30 simulée par n'importe quelle méthode géostatistique effectuée dans le système de coordonnées de dépôt (u, v, w) peut alors être effectué par les opérations suivantes : - récupération de l'élément du maillage contenant l'emplacement actuel donné (x, y, z), - en supposant que les coordonnées de dépôt (u, v, w) dans chacun des éléments de maillage 11 sont entièrement définies par leurs valeurs aux noeuds dudit élément, le calcul des coordonnées de dépôt (u, v, w) d'un emplacement actuel donné (x, y, z) utilisant un système d'interpolation analytique basé sur les valeurs des coordonnées de dépôt aux noeuds 10 de l'élément contenant, - récupération de l'élément du maillage de haute résolution 3D 35 de l'espace de dépôt sur lequel la méthode géostatistique a été appliquée, qui contient l'emplacement calculé (u, v ,w), - récupération de la valeur de la propriété simulée liée à cet élément. Dans une seconde version, le domaine géologique 1 peut être représenté par une série M de maillages 2D indépendants (par exemple des triangles 12), chacun de ces maillages représentant une limite 4 de l'ensemble B. La génération des maillages 2D ne fait pas partie du champ de cette invention et peut être effectuée en suivant toute méthode classique quelconque connue dans la matière. M peut être représentée par exemple par une liste linéaire M = {M1 ,M2 ,M3 ..., Mn). La structure de données représentant un tel maillage triangulaire Mj est codée de sorte que des cases mémoires sont réservées pour stocker et récupérer une valeur entière A(Mj). La valeur de A(Mj) est définie de sorte que si le maillage Mj représente la limite Bi, alors A(Mj) est égal à l'identifiant géologique correspondant bi (égal à i par exemple). La valeur de A(Mj) peut être définie pendant l'étape d'initialisation prenant place au début de l'étape de mise en correspondance des limites. Ainsi, la récupération de tous les triangles du maillage qui correspondent à une limite donnée Bi du domaine géologique peut comprendre les étapes suivantes : - d'abord, création d'une liste linéaire vide d'éléments Ti, - ensuite, récupération de l'identifiant géologique bi correspondant à Bi (égal à i par exemple), - enfin, exploration de tous les maillages de l'ensemble M, et si pour le maillage actuel d'indice j, A(Mj) = bi, alors insertion de tous les éléments de Mj dans la liste Ti. Dans cette seconde version, les déplacements inconnus (dx, dy, dz) sont déterminés grâce à la méthode des éléments frontière (BEM). La méthode des éléments frontière est une méthode -24- numérique classique, bien connue dans la matière, elle ne sera donc que partiellement décrite ici afin de ne pas brouiller les différentes versions de l'invention. Une description détaillée de la méthode se trouve par exemple dans le livre d'Aliabadi, "The Boundary Element Method: Applications in Solids and Structures - Vol. 2", aux éditions Wiley (2002). Dans le contexte de la géomécanique, une représentation par éléments frontière s'applique lorsque les équations aux dérivées partielles gérant la déformation des roches sont homogènes par morceaux, et pour lesquelles une fonction de Green peut être trouvée ou approximée. Dans de tels cas, le déplacement inconnu du système peut être réduit au déplacement des interfaces séparant les couches 3 ayant des propriétés mécaniques homogènes. La méthode des éléments frontière utilise des conditions aux limites données pour faire se conformer l'équation intégrale aux valeurs limites. La méthode BEM suppose que la géométrie actuelle du domaine géologique est représentée par une série de maillages 2D, chacun de ces maillages représentant une limite 4 dudit domaine géologique (voir figure 5), et que les éléments 12 de ces maillages 2D sont des polygones. Les fonctions utilisées pour approximer la solution au sein du domaine géologique étant choisies de sorte à satisfaire exactement aux équations, aucun élément de volume (3D) n'étant requis. Chaque polygone 12 individuel formant les maillages 2D peut être associé à des conditions aux limites de traction et de déplacement, et lorsque de telles conditions aux limites sont précisées, la méthode résout les déplacements inconnus correspondant. Les polygones peuvent être des triangles, mais les méthodes des éléments frontière habituelles peuvent être utilisées avec d'autres types de polygones. Dans cette version, appelée la version BEM, chaque couche du domaine géologique est associée à au moins deux paramètres représentant les propriétés mécaniques du matériau contenu dans la portion du domaine géologique correspondant à l'intérieur de cette couche, définissant ainsi un matériau homogène pour cette couche. Le nombre de paramètres dépend du modèle de déformation solide mis en oeuvre. Un modèle de déformation élastique exige deux paramètres, tandis que des modèles de déformation plastique ou des déformations plus complexes demandent davantage de paramètres. Les valeurs des paramètres peuvent varier en fonction du temps pendant l'étape de simulation mécanique. Si le matériau est supposé se déformer élastiquement pendant l'étape de simulation mécanique, ces paramètres sont généralement le coefficient de Poisson et le module de Young (ou les coefficients de Lamé correspondants) pour le matériau contenu dans la couche. -25- Afin de respecter la propriété « P1 », une condition aux limites est établie pour obliger tous les noeuds 10 des éléments 12 correspondant à un horizon donné 5 du domaine géologique 1 à avoir la même valeur finale w(x, y, z) de sorte à apparaître dans un même plan dans l'espace de dépôt tel qu'illustré en figure 5(c). Cette valeur peut être un résultat de l'étape de simulation mécanique, ou elle peut être précisée explicitement. Afin de respecter la propriété « P2 », une autre condition aux limites est établie imposant que tous les noeuds 10 des éléments 12 correspondant à un côté donné d'une faille donnée 7 ou un horizon 5 du domaine géologique 1 soient libres de mouvement à condition qu'ils restent en contact avec les éléments correspondant à l'autre côté de la faille 7 ou de l'horizon 5. En option, une autre condition aux limites peut être précisée imposant que tous les noeuds 10 situés aux intersections entre un horizon 5 donné du domaine géologique 1 et un côté donné d'une faille donnée 7 du domaine géologique 1 soient libres de mouvement à condition qu'ils restent en contact avec les bords correspondant à l'intersection entre l'horizon 5 et l'autre côté de la faille 7. La figure 5(d) montre un exemple de la façon dont une faille 7 apparaît dans l'espace de dépôt, les deux côtés de celle-ci restant mécaniquement en contact sans interpénétration, mais libres de glisser. Le calcul des déplacements inconnus (dx, dy, dz) à l'emplacement de chaque noeud 10 du maillage 2D représentant la géométrie actuelle du domaine géologique 1 est effectué par la méthode des éléments frontière basée sur les conditions aux limites définies dans l'étape de mise en place de la simulation mécanique précédente. Chaque séquence 8 identifiée dans le domaine géologique 1 peut être traitée séparément. Le calcul BEM n'est pas limité à une hypothèse de petits déplacements et peut également gérer de grandes déformations, évitant dans ce cas les effets secondaires indésirables liés à des rotations. Optionnellement, d'autres déplacement inconnus (dx, dy, dz) peuvent être calculés à des emplacements actuels (x, y, z) arbitraires appelés « points d'observation », utilisant les déplacements calculés aux emplacements de chaque noeud 10 des maillages 2D représentant la géométrie actuelle du domaine géologique 1, tel que décrit dans toute méthode quelconque des éléments frontière connue dans la matière. Un exemple de l'application de la méthode des éléments frontière est donné ci-dessous pour le calcul de la déformation mécanique de solides, en fonction d'un modèle de déformation élastique. Seul un aperçu est donné car la méthode des éléments frontière est bien connue. -26- Cette méthode suppose que les équations de comportement et d'équilibre définissant le problème mécanique utilisent deux opérateurs ou matrices T(p, q) et U(p, q), où T(p, q) est l'opérateur de traction donnant les tractions t au point de « terrain » q causé par une charge unitaire au point « source » p, et U(p, q) est une fonction de Green donnant le déplacement cl au point de « terrain » q causé par une charge unitaire au point de source p. Par exemple, ces deux opérateurs peuvent être donnés par les solutions fondamentales de l'élasticité linéaire. Le théorème de Betti donne le déplacement inconnu d au point p : `.U(p q)ds ù (p, q) Ls où S est la limite du solide considéré. En utilisant une discrétisation de S en éléments polygonaux, on obtient le système d'équations suivant : Mr .d = Mu. t avec t l'ensemble de tractions imposées, et MT et Mu deux matrices denses et non symétriques. Ainsi, la quantité de mémoire nécessaire pour stocker le système d'équations est de l'ordre de n et la complexité du résolveur itératif est de l'ordre de k-n2, avec n le nombre d'inconnues et k le nombre d'itérations requises pour résoudre le problème. Toutefois, des techniques précises, telles que l'expansion multipole (décrite par exemple dans l'oeuvre de Geengard et Rokhlin, "A Fast Algorithm for Particle Simulations", 3. Comput. Phys., Vol. 73, 1987) ou les matrices hiérarchiques (telles que décrites par exemple dans l'oeuvre de Hackbusch, "A Sparse Matrix Arithmetic Based on H-Matrices - Part I -Introduction to H-Matrices", Computing, Vol. 86, No 4, 1999) peuvent être utilisées pour réduire la complexité de l'ordre de n-log(n). Enfin les déplacements nodaux résolus (dx, dy, dz) sont utilisés pour définir le système de coordonnées de dépôt (u, v, w) à l'emplacement de chaque noeud 10 des maillages 2D représentant la géométrie actuelle du domaine géologique 1 et/ou à n'importe quel point d'observation. En se référant aux figures 6 et 7, le calcul de la valeur à tout emplacement actuel (x, y, z) d'une propriété 30 simulée par n'importe quelle méthode géostatistique effectuée dans le système de coordonnées de dépôt calculé (u, v, w) peut alors être effectué avec les opérations suivantes : - calcul des coordonnées de dépôt (u, v, w) d'un emplacement actuel donné (x, y, z), qui a été calculé précédemment en tant que point d'observation, -27- - récupération du maillage 3D haute résolution 35 de l'espace de dépôt sur lequel la méthode géostatistique a été appliquée, et qui contient l'emplacement calculé (u, v, w), - récupération de la valeur de la propriété simulée liée à cet élément. De façon plus générale, l'invention comprend également une méthode de post-traitement pour la « peinture » des éléments d'un maillage 3D 9 représentant la géométrie actuelle d'un domaine géologique 1 avec les propriétés 30 simulées dans l'espace de dépôt qui a été produit par la méthode selon cette invention. Elle comprend les étapes suivantes: - pour chaque élément e d'un maillage 3D 9, partitionnant le volume du domaine géologique actuel, calcul de l'emplacement actuel b de son barycentre, - calcul des coordonnées de dépôt b' correspondant à l'emplacement actuel b décrit précédemment, utilisant soit une méthode FEM ou BEM, - récupération de l'élément c d'un maillage 3D (généralement) à haute-résolution 35, partition du volume de l'espace de dépôt 2 auquel la méthode géostatistique a été appliquée et qui contient l'emplacement calculé b', - récupération de la valeur de la propriété simulée 30 attachée audit élément c, - enfin, association de cette valeur à l'élément e. Des illustrations de ce procédé sont montrées aux figures 6 et 7. Dans la figure 6(b) un maillage 3D 9 (qui peut être celui déjà utilisé pour le calcul de l'espace de dépôt 2) est peuplé de propriétés 30 calculées dans l'espace de dépôt 2 muni d'un maillage haute densité 35 comme illustré à la figure 6(a). Dans la figure 7(b), les propriétés 30 calculées dans l'espace de dépôt 2 sont intégrées dans le domaine géologique suivant un nouveau maillage 3D, essentiellement structuré, qui convient mieux pour l'application. Dans la figure 7(c), les séquences sont affichées séparément pour une meilleure compréhension. Toute version de méthode (ou portion de celle-ci) décrite dans le présent document peut être mise en oeuvre par des instructions programmées. Les instructions programmées peuvent être enregistrées sur des mémoires ordinateurs variées. Les instructions programmées peuvent être lues et exécutées par tout ordinateur ou ensemble d'ordinateurs pour mettre en oeuvre la version de la méthode (ou une portion de celle-ci). Bien que cette invention ait été décrite en plusieurs versions, de nombreuses alternatives, modifications et variations sembleront ou sont évidentes pour ceux qui sont experts dans la matière. -28- Ainsi, il est envisagé d'englober toutes ces alternatives, modifications, équivalences et variations qui sont dans l'esprit et la portée de cette invention, l'objectif essentiel étant d'utiliser la théorie de la mécanique et de la déformation des roches pour bâtir un espace de dépôt, où (i) les isochrones (horizons) identifiées dans un domaine géologique 3D sont essentiellement planes et parallèles, (ii) de telle sorte que tous les points de l'espace de dépôt situé à l'intérieur de la séquence stratigraphique n'ont qu'un et un seul emplacement correspondant dans l'espace actuel, et (iii) où la géométrie du domaine géologique 3D est physiquement fiable, l'espace de dépôt étant ainsi utilisé en tant qu'espace de calcul pour des applications telles que la modélisation et l'interpolation exacte des propriétés physiques du sous-sol, la génération de maillages pour des géométries complexes de structures géologiques, ou la vérification de la qualité d'une interprétation d'une sismique réflexion en 3D.

Claims (25)

1. REVENDICATIONS1. Une méthode pour la construction d'un espace de dépôt (2) correspondant à un domaine géologique (1) comprenant les étapes suivantes : - partition dudit domaine géologique (1) avec au moins un maillage conforme (9) correspondant sensiblement aux limites (4) dudit domaine géologique (1), - calcul des coordonnées de dépôt définissant un espace de dépôt (2), dans lequel le calcul de ces coordonnées de dépôt comprend des calculs de champs de déplacements, comprenant des simulations de déformations mécaniques dans le domaine géologique (1), au moyen d'un modèle de déformation solide.
2. 2. Une méthode selon la déclaration 1, (i) ledit maillage conforme (9) étant muni d'un système de coordonnées (x, y, z), (ii) des coordonnées de dépôt (u, v, w) étant définies comme u = u(x, y, z), v = v(x, y, z), et w= w(x, y, z), et, (iii) les champs de déplacements (dx, dy, dz) étant définis de telle sorte que : u(x, y, z) = x + dx(x, y, z), v(x, y, z) = y + dy(x, y, z), w(x, y, z) = z + dz(x, y, z). où le calcul des coordonnées de dépôt (u, v, w) comprend :-29- - une étape de mise en correspondance des limites, comprenant l'association de chacune des limites (4) du domaine géologique (1) à une série de facettes du maillage ou à des éléments correspondant sensiblement à la géométrie de la limite (4), une telle série comprenant au moins une facette d'un maillage conforme 3D, ou un élément d'un maillage conforme 2D (12), - une étape de préparation de la simulation mécanique, comprenant la définition des conditions aux limites visant à contraindre la géométrie calculée du domaine géologique dans l'espace de dépôt (2), - une étape de simulation mécanique, comprenant le calcul dans le domaine géologique (x, y, z) d'un champ de déplacements (dx, dy, dz) correspondant auxdites conditions aux limites.
3. 3. Une méthode selon la déclaration 2, selon laquelle l'étape de mise en place de la simulation mécanique comprend par ailleurs la définition d'une condition aux limites obligeant les noeuds (10) correspondant à un horizon donné (5) du domaine géologique (1) à avoir sensiblement la même valeur finale w(x, y, z).
4. 4. Une méthode selon la déclaration 3, selon laquelle la valeur w(x, y, z) finale correspondant à un horizon donné (5) est l'une des valeurs suivantes : - le résultat de la simulation mécanique, - une valeur spécifiée explicitement.
5. 5. Une méthode selon les déclarations 3 ou 4, selon laquelle la condition aux limites oblige également à ce que l'angle entre les paires de facettes ou éléments adjacents correspondant à l'horizon (5) soit sensiblement égal à un arc de 180 degrés.
6. 6. Une méthode selon l'une quelconque des déclarations 3 à 5, selon laquelle l'étape de préparation de la simulation mécanique comprend par ailleurs la définition d'une condition aux limites imposant que les noeuds (10) correspondant à un côté donné d'une faille donnée (7) du domaine géologique (1) soient libres de mouvement à condition qu'ils restent en contact avec l'autre côté de ladite faille (7).-30-
7. 7. Une méthode selon la déclaration 6, selon laquelle l'étape de préparation de la simulation mécanique comprend par ailleurs la définition d'une condition aux limites imposant que les noeuds (10) situés à l'intersection entre un horizon donné (5) du domaine géologique (1) et à un côté donné d'une faille donnée (7) du domaine géologique (1) soient libres de mouvement à condition qu'ils restent en contact avec les bords correspondant à l'intersection entre ledit horizon (5) et l'autre côté de ladite faille (7).
8. 8. Une méthode selon l'une quelconque des déclarations 3 à 7, selon laquelle l'étape de préparation de la simulation mécanique comprend par ailleurs la définition d'une condition aux limites imposant que les noeuds (10) correspondant à un côté donné d'un horizon donné (5) du domaine géologique (1) soient libres de mouvement à condition qu'ils restent en contact avec l'autre côté dudit horizon (5).
9. 9. Une méthode selon l'une quelconque des déclarations 3 à 8 selon laquelle les déplacements (dx, dy, dz) sont calculés aux noeuds (10) du (des) maillage(s) (9) représentant la géométrie du domaine géologique (1).
10. 10. Une méthode selon l'une quelconque des déclarations 3 à 9, selon laquelle les séquences (8) identifiées dans le domaine géologique (1) sont traitées séparément.
11. 11. Une méthode selon l'une quelconque des déclarations 3 à 10, selon laquelle l'étape de simulation mécanique comprend également l'utilisation d'une formulation « maître/esclave » pour contraindre les déplacements (dx, dy, dz) des noeuds (10) à rester mécaniquement en contact avec les séries d'éléments, de facettes ou de bords tel que défini dans les conditions aux limites.
12. 12. Une méthode selon la déclaration 11, selon laquelle la formulation maître/esclave est basée sur une méthode de pénalisation.
13. 13. Une méthode selon la déclaration 11, selon laquelle la formulation maître/esclave est basée sur l'utilisation des multiplicateurs de Lagrange.-31-
14. 14. Une méthode selon l'une quelconque des déclarations 2 à 13, selon laquelle le champ des déplacements (dx, dy, dz) est calculé avec la méthode des éléments finis (FEM), la géométrie du domaine géologique étant représentée par au moins un maillage 3D conforme (9).
15. 15. Une méthode FEM selon la déclaration 14, selon laquelle les maillages 3D conformes (9) sont hybrides, et comprennent des groupes d'éléments structurés (11) et des groupes d'éléments non structurés (11).
16. 16. Une méthode FEM selon la déclaration 14 ou 15, selon laquelle les éléments (11) des maillages 3D conformes (9) sont associés à des paramètres représentant les propriétés mécaniques locales du matériau contenu dans la portion du domaine géologique (1) renfermé dans ledit élément (11).
17. 17. Une méthode FEM selon la déclaration 16, selon laquelle le modèle de déformation d'un solide est un modèle de déformation élastique.
18. 18. Une méthode FEM selon l'une quelconque des déclarations 14 à 17, selon laquelle les déplacements (dx, dy, dz) sont calculés selon l'une quelconque des méthodes suivantes : - en construisant une matrice de rigidité globale, - itérativement en construisant et en mettant à jour des matrices de rigidité locales définies aux noeuds (1) du maillage 3D (9) représentant la géométrie du domaine géologique (1).
19. 19. Une méthode FEM selon l'une quelconque des déclarations 14 à 18, comprenant également le calcul à l'emplacement (x, y, z) dans le domaine géologique (1) de la valeur d'une propriété de terrain (30) calculée à l'emplacement (u, v, w) dans l'espace de dépôt (2), le calcul de ladite propriété de terrain comprenant les étapes suivantes : - récupération de l'élément de maillage (11) contenant l'emplacement (x, y, z),-32- - calcul des coordonnées de dépôt (u, v, w) dudit emplacement (x, y, z) au moyen d'un système d'interpolation analytique basé sur les valeurs des coordonnées de dépôt (u, v, w) correspondant aux coordonnées de noeuds (x, y, z) dudit élément de maillage contenant (11), - récupération de la valeur de ladite propriété de terrain (30) correspondant aux coordonnées (u, v, w) dans l'espace de dépôt (2).
20. 20. Une méthode selon l'une quelconque des déclarations 2 à 12, selon laquelle le champ des déplacements (dx, dy, dz) est calculé avec la méthode des éléments frontière (BEM), la géométrie du domaine géologique (1) étant représentée par une série de maillages 2D conformes, chacun de ces maillages représentant une limite (4) dudit domaine géologique (1).
21. 21. Une méthode BEM selon la déclaration 20, selon laquelle les maillages 2D conformes (9) sont hybrides, et comprennent des groupes d'éléments structurés (12) et des groupes d'éléments non structurés (12).
22. 22. Une méthode BEM selon les déclarations 20 ou 21, selon laquelle le modèle de déformation d'un solide est un modèle de déformation élastique.
23. 23. Une méthode BEM selon l'une quelconque des déclarations 20 à 22, selon laquelle l'étape de simulation mécanique comprend par ailleurs le calcul des déplacements (dx, dy, dz) à des emplacements arbitraires (x, y, z) appelés « points d'observation », utilisant les déplacements calculés (dx, dy, dz) à l'emplacement des noeuds (10) des maillages 2D (12) représentant la géométrie du domaine géologique (1).
24. 24. Une méthode BEM selon l'une quelconque des déclarations 20 à 23, comprenant également le calcul à l'emplacement (x, y, z) dans le domaine géologique (1) de la valeur d'une propriété de terrain (30) calculée à l'emplacement (u, v, w) dans l'espace de dépôt (2), le calcul de ladite propriété de terrain comprenant les étapes suivantes : - calcul des coordonnées de dépôt (u, v, w) dudit emplacement (x, y, z) utilisé comme point d'observation,-33- - récupération de la valeur de ladite propriété de terrain (30) correspondant aux coordonnées (u, v, w) dans l'espace de dépôt (2).
25. 25. Un programme informatique mettant en oeuvre la méthode selon l'une quelconque des déclarations précédentes.