FR 3014581 A1 20150612 FR 1362311 A 20131209 PROCEDE ET DISPOSITIF DE DETERMINATION D'UN ENDOMMAGEMENT D'UN MODELE GEOLOGIQUE La présente invention concerne le domaine de la détermination d'un endommagement d'un élément de formation géologique susceptible de modéliser un materiau de type roche subissant des évolutions structurelles telles que des déformations, des fissurations, des fragmentations avec interactions entre les fragments, dans le cadre de la recherche et de exploitation d'hydrocarbures, liquides ou gazeux issus de la pyrolyse de kérogène. Le matériau modélisé est par exemple un schiste bitumineux qui est un matériau hétérogène constitué de roche mère (encore appelée kérogène) et de matières minérales. Le kérogène peut être amené à subir une pyrolyse. naturelle ou provoquée, au cours de laquelle il se transforme partiellement en hydrocarbures, liquides et gazeux. associés à d'autres composé chimiques tels que CO2, H20, etc.. De tels changements de phase du kérogène, principalement par libération de gaz sous pression, génèrent des contraintes importantes dans le matériau et sont la. source de déformations, de fissurations et de fragmentations du matériau. La modélisation de l'endommagement d'un matériau subissant de telles évolutions structurelles dans une formation géologique donnée est utile dans le cadre de la recherche et de la production d'hydrocarbures ou de gaz dans les schistes bitumineux car celle-ci permet d'estimer l'effet d'une pyrolyse provoquée sur le kérogène contenu dans la formation, et ainsi les ressources disponibles associées à ladite formation, la production future et, donc, la valeur économique du champ. Il existe plusieurs approches pour modéliser le comportement des matériaux sous contraintes que l'on peut regrouper en deux grandes familles : les approches continues et les approches discrètes. Les approches continues prennent en compte le comportement des différentes phases du matériau et la physique du problème sous la forme de lois de comportement formulées en équations aux dérivées partielles. Des méthodes de résolution par cléments finis permettent de résoudre ces équations. Les inconvénients des approches continues sont nombreux lorsqu'on l'on cherche à les appliquer aux matériaux hétérogènes comme les schistes bitumineux. Ces approches nécessitent un maillage extrêmement fin pour pouvoir prendre en compte les hétérogénéités micro-structurelles des schistes bitumineux et donc un temps de calcul conséquent. En outre, les lois de comportement dans les éléments finis ne peuvent pas répondre correctement au problème d'endommagement des matériaux dès lors que les phénomènes tels que l'ouverture et la propagation de fissures jouent un râle prépondérant, ce qui est en particulier le cas au cours de l'évolution structurelle des schistes bitumineux. Les approches discrètes ont été développées dans le but de modéliser des matériaux divisés non cohésifs tels que les matériaux granulaires. Dans ces méthodes, le matériau est décrit comme une collection de corps rigides et indépendants. les particules, qui interagissent entre eux à travers des lois de forces reliant la force de contact à une variable de déformation associée au contact (définie à travers des degrés de liberté rigides des particules). L'évolution du système est alors obtenue par l'intégration des équations du mouvement. Les interactions considérées. telles que la cohésion, le contact ou le frottement, diffèrent selon les méthodes. leurs implantations et l'application visée. Les approches discrètes se regroupent sous l'appellation de méthode « DEM » acronyme anglo-saxon pour « Discret Element Method -. par opposition aux méthodes des éléments finis. Les méthodes « DEM » sont de façon générale peu adaptées à la modélisation des matériaux continus. En effet, la déformation et l'endommagement des particules ne sont pas considérés et la propagation des fissures se fait uniquement par la rupture des liens cohésifs et par contournement des particules. De plus. le comportement global du matériau est le résultat d'un grand nombre d'interactions à l'échelle des particules. Les propriétés mécaniques du matériau sont donc des propriétés résultantes et ne peuvent être directement introduites dans la formulation discrète comme c'est le cas dans les méthodes continues. En outre. en trois dimensions, l'utilisation de particules de forme sphérique offre des performances numériques optimales en raison de !a facilité de gestion des contacts mais introduit une phase de vide artificiel dans le matériau. L'utilisation d'autres formes de particules nécessite la gestion de contacts plus complexes et donc un temps de calcul plus important. Enfin, la configuration géométrique initiale du domaine discret peut avoir un grand impact sur le comportement mécanique du matériau et des phénomènes tels que la propagation de fissures suivant des chemins préférentiels ou encore la diffusion non uniforme des ondes élastiques peuvent etre observé. Une troisième approche a été proposée pour surmonter ces inconvénients, appelée discrétisation sur réseau (en anglais (, lattice model »). Dans cette approche, le matériau n'est pas représenté par un ensemble d'éléments volumiques, comme dans l'approche continue, ou par un ensemble de corps rigides étendus en contact. comme dans l'approche discrète, mais par une distribution de points, les noeuds, reliés entre eux par des lois d'interaction, les liens. Seules les lois locales, comme l'équilibre des forces et des moments. sont considérées, et leur implémentation se fait au niveau de chaque noeud relié à un nombre limite et défini de voisins. Cette approche numérique permet de discrétiser !es phases du matériau et éventuellement leurs interfaces sur un réseau d'éléments unidimensionnels, les liens, pouvant avoir différents comportements rhéologiques. La complexité de la microstructure et le comportement de ses phases sont ainsi pris en compte. Les liens qui relient les noeuds permettent de modéliser le comportement mécanique des différentes phases géologiques du matériau et ont usuellement un comportement rhéologique élastique fragile caractérisé par une raideur liée au module élastique de la phase géologique en question et un seuil de rupture en force lié à la contrainte de rupture de ladite phase géologique. La formation de fissures dans le modèle géologique est alors prise en compte à travers la suppression de liens entre des noeuds lorsqu'une force appliquée à un lien est supérieure à une valeur seuil de force du lien. Cette méthode présente cependant des inconvénients et ne permet pas de modéliser de façon satisfaisante l'endommagement du modèle géologique au-delà de la simple ouverture d'une fissure. Ainsi en particulier, une fois un lien rompu, les deux noeuds reliés à l'origine par ledit lien, et formant alors les deux lèvres de la fissure, n'interagissent plus entre eux. Cette méthode n'est ainsi pas à même de simuler la propagation des fissures les déformations du matériau. les fissurations multiples et les fragmentations avec interactions entre les fragments du matériau. Il existe ainsi un besoin pour une méthode de détermination d'un endommagement d'un modèle géologique capable de prendre en compte des hétérogénéités micro-structurelles, de modéliser de façon fiable la déformation des phases géologiques du matériau. et de simuler des phénomènes tels que l'ouverture et la propagation de fissures qui jouent un rôle prépondérant au cours de l'évolution structurelle des schistes bitumineux. La présente invention vient améliorer la situation. A cet effet. la présente invention propose d'utiliser une approche mixte fondée sur la discrétisation sur réseau et sur la dynamique moléculaire pour modéliser le comportement mécanique et l'endommagement d'un modèle géologique. 10 La présente invention vise alors un procédé de détermination d'un endommagement d'un modèle géologique, ledit modèle géologique comportant une pluralité de noeuds ainsi qu'une pluralité de liens, chaque lien de la pluralité de liens reliant ensemble deux noeuds de la pluralité de noeud. Le procédé comprend !es étapes, 15 - pour chaque lien de la pluralité de liens : * détermination d'une force appliquée audit lien, * si ladite force est supérieure à une valeur seuil de force : - pour chaque noeud des deux noeuds de la pluralité de noeud reliés ensemble par ledit lien, création d'une sphère de contact 20 associée audit noeud, et - suppression dudit lien, - détermination d'une pluralité de contacts entre sphères comportant au moins un contact entre sphères chaque contact entre sphères de la pluralité de contacts entre sphères associant deux sphères de contact parmi la pluralité de sphères de 25 contact. Les « sphères de contact » sont des corps ayant une certaine extension spatiale et aptes à interagir deux à deux par des interactions de contact : les « contacts entre sphères ». Ainsi, contrairement aux noeuds, les voisins de chaque sphère de contact, avec lesquels ladite sphère de contact interagit, ne sont pas déterminés a priori mais peuvent évoluer dans le temps au gré de l'endommagement du modèle géologique. De cette façon, il est par exemple possible de modéliser précisément le déplacement relatif des lèvres d'une fissure. Il est également possible de modéliser l'interaction de fragments du modèle géologique entre eux. Dans un mode de réalisation, pour chaque contact entre sphères de la pluralité de contacts entre sphères. le procédé peut comprendre - détermination d'une force d'interaction entre les deux sphères de contact associées par ledit contact entre sphères. - pour chaque lien de la pluralité de liens tel que, parmi les deux noeuds de la pluralité de noeud reliés ensemble par ledit lien, au moins un noeud est associé à une sphère de contact associée audit contact entre sphères, transmission de ladite force d'interaction audit lien. Dans un mode de réalisation du procédé, à chaque contact entre sphères de la pluralité de contacts entre sphères est associé - un plan de contact entre sphères, - une composante normale de la force d'interaction. normale au plan de contact entre sphères, - une composante tangentielle de la force d'interaction, comprise dans le plan de contact entre sphères, et - la composante normale de la force d'interaction est une composante proportionnelle à une distance d'interpénétration des deux sphères de contact associées par ledit contact entre sphères. et - la composante tangentielle de la force d'interaction est une composante de direction opposée à un déplacement relatif tangentiel des deux sphères de contact associées par ledit contact entre sphères et de module proportionnel audit déplacement relatif tangentiel. Dans un mode de réalisation du procédé, une valeur de raideur, une valeur seuil de force et une valeur de longueur de référence sont associées à chaque lien de la pluralité de liens. lesdites valeurs de raideur et valeurs seuil de force étant déterminées a partir de valeurs de module élastique E et de valeurs de seuil de résistance à la traction o d'au moins une phase géologique modélisée par le modèle géologique. Dans un mode de réalisation du procédé, les liens de la pluralité de liens et les noeuds de la pluralité de noeuds sont disposés en un réseau régulier désordonné, de préférence en un réseau tel qu'une distribution statistique de valeurs de longueur de référence ou une distribution statistique de valeurs de raideur, associées aux liens de la pluralité de liens, est une distribution statistique unimodale. Dans un mode de réalisation. le procédé comprend en outre les étapes, - application d'une pluralité de conditions aux limites à une pluralité de noeuds de limites parmi la pluralité de noeuds. - pour chaque noeud de la pluralité de noeuds, détermination d'une position spatiale dudit noeud. Dans un mode de réalisation du procédé, la détermination d'une position spatiale pour chaque noeud de la pluralité de noeuds comprend soit une minimisation d'une énergie potentielle du modèle géologique, soit une intégration d'une pluralité d'équations du mouvement associées aux noeuds. Dans un mode de réalisation du procédé : - une première pluralité de noeuds parmi la pluralité de noeuds du modèle géologique est associée à une première phase géologique modélisée par le modèle 25 géologique, - une deuxième pluralité de noeuds parmi la pluralité de noeuds du modèle géologique est associée à une deuxième phase géologique modélisée par le modèle géologique, - une pluralité de liens de première phase géologique parmi la pluralité de liens du modèle géologique est constituée de liens de !a pluralité de liens reliant ensemble deux noeuds de la première pluralité de noeuds, à chaque lien parmi la pluralité de liens de première phase géologique est associe une valeur de raideur de première phase géologique et une valeur seuil de force de première phase géologique déterminées à partir de valeurs de module élastique E et de valeurs de seuil de résistance à la traction o de la première phase géologique modélisée par le modèle géologique, - une pluralité de liens de deuxième phase géologique parmi la pluralité de liens du modèle géologique est constituée de liens de la pluralité de liens reliant ensemble deux noeuds de la deuxième pluralité de noeuds. à chaque lien parmi la pluralité de liens de deuxième phase géologique est associé une valeur de raideur de deuxième phase géologique et une valeur seuil de force de deuxième phase géologique déterminées à partir de valeurs de module élastique E et de valeurs de seuil de résistance à la traction o de la deuxième phase géologique modélisée par le modèle géologique, et - une pluralité de liens d'interface parmi la pluralité de liens du modèle géologique est constituée de liens de la pluralité de liens reliant ensemble un noeud de la première pluralité de noeuds et un noeud de la deuxième pluralité de noeuds, chaque lien parmi la pluralité de liens d'interface est associé une valeur de raideur d'interface et une valeur seuil de force d'interface , et une valeur de raideur de première phase géologique est différente d'une valeur de raideur de deuxième phase géologique et/ou une valeur de seuil de force de première phase géologique est différente d'une valeur de seuil de force de deuxième 30 phase géologique. Un dispositif destiné à déterminer un endommagement d'un modèle géologique peut être avantageux, enlui-rnerne. dès lors qu'il permet e déterminer un endommagement d'un modèle géologique. Ainsi. la présente invention vise également un dispositif de détermination d'un 5 endommagement d'un modèle géologique. ledit modèle géologique comportant une pluralité de noeuds ainsi qu'une pluralité de liens, chaque lien de la pluralité de liens reliant ensemble deux noeuds de la pluralité de noeud. Le dispositif comporte un circuit adapté pour permettre, - pour chaque lien de la pluralité de liens 10 la détermination d'une force appliquée audit lien, * si ladite force est supérieure à une valeur seuil de force: - pour chaque noeud des deux noeuds de la pluralité de noeud reliés ensemble par ledit lien, la création d'une sphère de contact associée audit noeud, et 15 - la suppression dudit lien, - la détermination d'une pluralité de contacts entre sphères comportant au moins un contact entre sphères, ledit contact entre sphères associant deux sphères de contact parmi la pluralité de sphères de contact. 20 Un programme informatique, mettant en oeuvre tout ou partie du procédé décrit ci- avant, installé sur un équipement préexistant, est enlui-même avantageux, dès lors qu'il permet de déterminer un endommagement d'un modèle géologique. Ainsi, la présente invention vise également un programme informatique comportant des instructions pour la mise en oeuvre du procédé précédemment décrit, 25 lorsque ce programme est exécuté par un processeur. Ce programme peut utiliser n'importe quel angage de programmation (par exemple, un langage objet ou autre). et être sous ia forme d'un code source interprétable, d'un code partiellement compilé ou d'un code totalement compilé. La figure 2 décrite en détails ci-après, peut former l'organigramme de l'algorithme général d'un tel programme informatique. D'autres caractéristiques et avantages de l'invention apparaîtront encore à la lecture de la description qui va suivre. Celle-ci est purement illustrative et doit être lue en regard des dessins annexés sur lesquels - la figure 1 illustre un modèle géologique pour le fonctionnement d'un procédé possible selon ['invention ; - la figure 2 illustre un lien reliant deux noeud d'un modèle géologique pour le fonctionnement d'un procédé possible selon l'invention ; - la figure 3a illustre un dét d'un modèle géologique comportant un lien à rompre au cours de la mise en oeuvre d'un procédé possible selon l'invention; - la figure 3d illustre le détail d'un modèle géologique de la figure 3a comportant un lien à rompre et deux sphères de contact au cours de la mise en oeuvre d'un procédé possible selon l'invention: - la figure 3c illustre le détail d'un modèle géologique de la figure 3b comportant deux sphères de contact et dans lequel le lien à rompre a été supprimé au cours de la mise en oeuvre d'un procédé possible selon l'invention; - la figure 3d illustre un modèle géologique au cours de la mise en oeuvre d'un procédé possible selon l'invention dans lequel plusieurs liens à rompre ont été supprimés; - la figure 4 illustre un détail d'un modèle géologique comportant un contact entre sphère au cours de la mise en oeuvre d'un procédé possible selon l'invention; - la figure 5 illustre un diagramme de fonctionnemen d n procédé dans un mode de réalisation selon l'invention - la figure 6 illustre un dispos tif de calcul permettant la mise en oeuvre selon un mode de réalisation de l'invention. La figure 1 illustre un modèle géologique 100 pour le fonctionnement d'un procédé possible selon !invention. Un tel modèle géologique 100 est en particulier apte à modéliser, ou représenter, un matériau 1 de type roche, originellement consolidée, homogène ou hétérogène, et vouée à subir d'importantes évolutions structurelles telles que des déformations, des fissurations multiples. des fragmentations avec interactions entre les fragments. Le modèle géologique 100 est en particulier apte à représenter un matériau 1 comportant plusieurs phases géologiques la. lb de caractéristiques mécaniques distinctes. Par « un matériau comportant plusieurs phases géologiques -, on entend un matériau hétérogène constitué de plusieurs portions qui peuvent différer par exemple par leurs compositions chimiques, leurs caractéristiques mécaniques et/ou leur état (solide. fluide. gazeux). Le modèle géologique peut ainsi présenter des hétérogénéités de dimensions plus ou moins importantes. Le modèle géologique 100 comporte une pluralité de noeuds 200 ainsi qu'une pluralité de liens 300. Comme illustré sur les figures 1 et 2, chaque lien i de a pluralité de liens 300 relie ensemble deux noeuds i,j de la pluralité de noeud 200. Ainsi, le matériau 1 n'est pas représenté par un ensemble d'éléments volumiques, comme en éléments finis. mais par une distribution de points, les noeuds 200, reliés entre eux par des lois d'interaction, les liens 300. Les liens 300 permettent de transmettre les forces entre noeuds 200. Dans un mode de réalisation de invention, les liens 300 transmettent uniquement une force normale I; entre les noeuds 200. Les liens peuvent ainsi avoir par exemple un comportement de ressort fragile. On associe alors à chaque lien 1 de la pluralité de liens 300 une valeur de raideur ki, une valeur seuil de force Ï! et une valeur de longueur de référence lit,. La valeur seuil de force défini une force maximale applicable au lien i avant que ledit lien ne rompe comme il sera décrit plus détail ci-après. Les valeurs de longueur de référence lo sont par exemple définies en fonction de la plus petite hétérogénéité du modèle géologique 100, et de sorte à discrétiser de façon convenable ladite plus petite hétérogénéité, par exemple de sorte que ladite plus petite hétérogénéité comporte un nombre de noeud 200 prédéfini, de préférence supérieur ou égal à un noeud 200. La valeur de raideur ki et la valeur de longueur de référence ( permettent d'associer à chaque lien une force normale errn proportionnelle à un allongement Ai' du lien 1 tel que - Dans une variante, la force norma e Mon, peut être une fonction non linéaire de l'allongement Ar du lien, de la valeur de raideur ki et la valeur de longueur de référence Dans un autre mode de réalisation de invention les liens peuvent transmettre, en plus de la force normale fun,,,,, une force de cisaillement et un moment. Dans un tel mode de réalisation, les liens représentent alors des « poutres reliant les noeuds. Les valeurs de raideur k' et les valeurs seuil de force fr associées à chaque lien i de la pluralité de liens peuvent avantageusement être déterminées à partir des valeurs de module élastique L' et des valeurs de seuil de résistance à la traction 0- des phases géologiques la, lb modélisés par le modèle géologique 100. Comme illustré sur la figure 1, les noeuds de la pluralité de noeuds 200 et les liens de la pluralité de liens 300 sont disposés en un réseau 120. Dans un mode de réalisation de l'invention. ledit réseau 120 est un réseau tridimensionnel. Dans une variante, ledit réseau 120 peut être un réseau bidimensionnel ou même un réseau unidimensionnel. En particulier, les liens de la pluralité de liens 300 peuvent être disposés en un réseau 120 régulier désordonné. On entend par « réseau régulier désordonné » par exemple un réseau tel que les liens de la pluralité de liens 300 ont en moyenne la même valeur de longueur de référence /0 avec un écart type non nul. C'est-à-dire un réseau tel qu'une distribution statistique de valeurs de longueur de référence 1(') associées aux liens ide la pluralité de liens 300 est une distribution statistique unimodale, par exemple avec un mode unique 10, en particulier une distribution statistique unimodale avec un écart-type non nul. On entend également par « réseau régulier désordonné » un réseau tel que les liens de !a pluralité de liens 300 ont en moyenne la même valeur de raideur k avec un écart type non nu!. C'est-à-dire un réseau te! qu'une distribution statistique de valeurs de raideur k' associées aux liens ide la pluralité de liens 300 est une distribution statistique unimodale, par exemple avec un mode unique k, en particulier une distribution statistique unimodale avec un écart-type non nul. Pour cela, on peut par exemple disposer !es noeuds de la pluralité de noeuds 200 et les liens de la pluralité de liens 300 en un réseau 120 parfaitement régulier (avec un écart-type nul} puis appliquer un désordre contrôle aux positions desdits noeuds, aux valeurs de longueur de référence /0 des liens ou aux valeurs de raideur k des liens de sorte à constituer un réseau régulier désordonné. L'application d'un désordre contrôlé permet de supprimer, ou au moins de réduire, les directions privilégiées oie fracturation du modèle géologique 100 liées à la discrétisation en réseau 120. Ledit réseau 120 régulier tridimensionnel peut par exemple être un réseau tétraédrique. Comme mentionné ci-avant. le modèle géologique 100 peut en particulier modéliser un matériau 1 comportant plusieurs phases géologiques et comportant par exemple au moins une première phase géologique 1 a et une deuxième phase géologique I b. Comme illustré sur la figure 1, une première pluralité de noeuds 210 parmi la pluralité de noeuds 200 du modèle géologique 100 est ainsi associée à la première phase géologique la modélisée par le modèle géologique. Une deuxième pluralité de noeuds 220 parmi la pluralité de noeuds 200 du modèle géologique est ainsi associée à la deuxième phase géologique 1 b modélisée par le modèle géologique. Les mêmes considérations s'appliquent aux liens 300 de sorte qu'une pluralité de liens de première phase géologique 310 parmi la pluralité de liens du modèle géologique est constituée de liens de la pluralité de liens reliant ensemble deux noeuds de la première pluralité de noeuds 210. A chaque lien parmi la pluralité de liens de première phase géologique 310 est associé une valeur de raideur de première phase géologique ic4)' et une valeur seuil de force de première phase géologique La valeur de raideur de première phase géologique et la valeur seuil de force de première phase géologique 1,-1) peuvent en particulier être définies par une valeur de module élastique de la première phase géologique L'`' modélisée par le modèle géologique et une valeur de seuil de résistance à la traction de la première phase géologique 044 modélisée par le modèle géologique, De façon analogue, une pluralité de liens de deuxième phase géologique 320 parmi la pluralité de liens 300 du modèle géologique est constituée de liens de la pluralité de liens reliant ensemble deux noeuds de la deuxième pluralité de noeuds 220. A chaque lien parmi la pluralité de liens de deuxième phase géologique 320 est associé une valeur de raideur de deuxième phase géologique et une valeur seuil de force de deuxième phase géologique ez qui sont définies par une valeur de module élastique de la deuxième phase géologique Ll'Hodelisée par le modèle géologique et une valeur de seuil de résistance à la traction de la deuxième phase géologique modélisée par le modèle géologique. Enfin, on peut également définir des liens d'interface 330. Une pluralité de liens d'interface 330 parmi la pluralité de liens 300 du modèle géologique est ainsi constituée de liens de la pluralité de liens reliant ensemble un noeud de la première pluralité de noeuds 210 et un noeud de la deuxième pluralité de noeuds 220. A chaque lien parmi la pluralité de liens d'interface 330 est associé une valeur de raideur d'interface Int''Tince et une valeur seuil de force d'interface r f w-c Les caractéristiques mécaniques de la première phase géologique la et de la deuxième phase géologique 1 b peuvent différer. Ainsi. !a première phase géologique peut par exemple modéliser une roche mère tandis que la deuxième phase géologique peut modéliser du kérogène, un fluide ou même un gaz. Une valeur de raideur de première phase géologique et et/ou une valeur de seuil de force de première phase géologique peu(ven)t donc être différente(s) d'une valeur de raideur de deuxième phase géologique 1::" et/ou d'une valeur de seuil de force de deuxième phase géologique f ;72 Les caractéristiques mécaniques des interfaces, c'est-à-dire en particulier les valeurs ktnterf 'ire et "ner-f an' associées liens d'interface 330, peuvent être différentes des valeurs associées aux liens de la première phase géologique 210 et de la deuxième phase géologique 220 ou être identiques aux valeurs de l'une desdites phases géologiques la, 1 b. Les caractéristiques mécaniques des interfaces peuvent par exemple être calibrées par des mesures sur des échantillons test de matériau. La figure 5, en relation avec les figures 3a à 3d et 4. illustre un diagramme de fonctionnement d'un procédé dans un mode de réalisation selon l'invention. Le mode de fonctionnement décrit dans cette figure n'est qu'un des modes possibles de fonctionnement de I nvention. Une étape préliminaire du procédé (étape 1100, DEF_MOD) comprend la définition du modèle géologique 100 en associant à chaque noeud I de la pluralité de noeuds une position spatiale initiale r;'it et en définissant les valeurs de raideur, de seuil de force et de longueur de référence associées aux liens de la pluralité de liens. Ce modèle géologique 100 tel que défini est un modèle géologique initial qui va être modifié au cours de la mise en oeuvre du procédé qui cela va maintenant être détaillé. En se référant à la figure 2, une étape du procédé comporte ensuite l'application d'une pluralité de conditions aux limites au modèle géologique (étape 1200, APP_COND). Dans un mode de réalisation, les conditions aux limites peuvent être appliquées à une pluralité de noeuds de limites 240 parmi la pluralité de noeuds 200. Dans un premier mode de réalisation, es noeuds de limites 240 peuvent être des noeuds formant le contour spatial externe 110 du modèle géologique 100. Les noeuds de limites 240 ont alors usuellement un nombre de voisins inférieur aux autres noeuds du modèle géologique. Dans un deuxième mode de réalisation, les noeuds de limites 240 peuvent être des noeuds disposés à l'intérieur du modèle géologique 100, en particulier n'importe quels noeuds du modèle géologique 100. Dans ce deuxième mode de réalisation, on peut alors appliquer des conditions aux limites à une portion intérieure du matériau. A titre d'exemple. ce mode de réalisation permet ainsi de déterminer l'endommagement d'un matériau soumis à des efforts particuliers suite à un changement de phase d'une portion intérieure du matériau ou à une dilation ou contraction de ladite portion intérieure. par exemple en réponse à une élévation ou une diminution de température de ladite portion intérieure. L'application desdites conditions aux limites peut par exemple comporter un déplacement des positions spatiales desdits noeuds de limites 240. Dans un autre mode de réalisation, l'application desdites conditions aux limites peut par exemple comporterl'application d'une force auxdits noeuds de limites 24 Une étape ultérieure du procédé comporte la détermination d'une position spatiale des noeuds de la pluralité de noeuds (étape 1300, DET_POS). Une telle étape de détermination d'une position spatiale peut être mise en oeuvre de différentes façons. Dans un premier mode de réalisation, on détermine la position spatiale ri de chaque noeud I de la pluralité de noeuds 200 par minimisation d'une énergie potentielle El, du modèle géologique 100. Dans le repère local d'un noeud l'énergie potentielle d'un lien reliant le noeu un noeud j est donnée par : 1 '112 -2kiki oùr1 et r1' sont respectivement la position initiale et la position après un déplacement du noeud j. Une énergie potentielle totale est alors obtenue par sommation On peut par exemple minimiser une telle énergie potentielle Ep par une méthode des gradients conjugués telle que la méthode d'Hestenes et Stiefel. Dans un deuxième mode de réalisation, on détermine la position spatiale ri de chaque noeud t de la pluralité de noeuds 200 par intégration d'une pluralité 5 d'équations du mouvement associées aux noeuds. Dans cette méthode, on discretise temporellement révolution du modèle géologique 100. La dynamique des noeuds 200 peut alors être calculée à chaque instant t par mise en oeuvre d'un algorithme d'intégration qui consiste à déterminer, au temps t + At. pour chacun des noeuds i, la position rt. !a vitesse vi et 10 l'accélération ai dudit noeud t dans un repère global du modèle géologique 100, en fonction de la position. la vitesse et l'accélération des noeuds au temps t, ainsi que des forces appliquées et transmises aux noeuds par les liens. A cette fin, on peut par exemple mettre en oeuvre un algorithme tel que 15 l'intégration de Verlet. Plus précisément, pour chaque noeud i de la pluralité de noeuds la position ri et la vitesse yr, au temps t + At, peuvent être déterminées par intégration numérique des équations du mouvement associée audit noeud t. Pour chaque noeud L de la pluralité de noeuds, l'accélération ai, au temps t + At, 20 peut être déterminée par application de la seconde loi de Newton, c'est-à-dire en sommant les forces exercées sur ledit noeud t par les liens reliant ledit noeud à d'autres noeuds i de la pluralité de noeuds 200 du modèle géologique. Par exemple ntiai = ZJ fil où la somme est effectuée sur l'ensemble des noeuds ./ de a pluralité de noeuds reliés au noeud I par un lien de la pluralité de liens et où f" 25 est la force exercée par le noeud j sur le noeud i, et où ait est la masse du noeud t. La stabilité numérique du calcul est assurée si le pas de temps est petit devant un temps de réponse élastique ?-7717, où mi est ici la masse des noeuds 200 du modèle géologique 100, par exemple de la phase géologique 1 a du matériau 1 la plus rigide. La masse m des noeuds 200 d'une phase géologique du modèle géologique 100 peut être aisément obtenue à partie de la masse volumique de ladite phase géologique. Dans ce deuxième mode de réalisation, on peut en outre inclure un terme d'amortissement visqueux T permettant de dissiper les oscillations autoentretenues. Le terme d'amortissement visqueux T est par exemple obtenu à partir du coefficient de restitution du matériau ou de la phase géologique modélisé(e) par le modèle géologique. Le procédé comprend ensuite, pour chaque lien i de la pluralité de liens 300, la détermination d'une force appliquée au lien fc:7,2,. Dans le mode de réalisation dans lequel les liens sent des ressorts fragiles transmettant uniquement une force normale la force appliquée au lien i, est par exemple donnée par la force normale décrite ci-dessus fi1711 01'7? = kt Al' Sila force Lipp appliquée au lien i est supérieure à la valeur seuil de force f,? du lien e Hen r est un lien à rompre 305 et on met alors en oeuvre les deux opérations suivantes. En se référant aux figures 3a à 3d, tout d'abord, pour chacun des deux noeuds 201, 202 reliés ensemble par le lien à rompre 305. on réalise une opération de création d'une sphère de contact 400 associée audit noeud (étape 1510, CRT SPH). _ Une première sphère de contact 410 et une deuxième sphère de contact 420 sont ainsi créées et associées respectivement au premier noeud 201 et au deuxième noeud 202 reliés ensembles par le lien à rompre 305. Les sphères de contact 400 forment une pluralité de sphères de contact. Les sphères de contact 400 comportent chacune un centre 430 ayant une position spatiale de centre par exemple similaire à la position spatiale r1 du noeud i associé à la sphère de contact 400. Les sphères de contact 400 comportent aussi un rayon de référence prédéfini I?. Le rayon de référence R peut par exemple être égal à la moitié de la longueur de référence du lien à rompre. Le rayon de référence peut également être inférieur à la moitié de !a longueur de référence Io du lien à rompre 350, par exemple inférieur à la moitié de la longueur de référence lo du lien à rompre 350 d'une valeur prédéterminée LII? ou présentant un ratio prédéterminé avec la moitié de la longueur de référence 10 du lien à rompre 350. Dans une variante de l'invention, les sphères de contact 400 peuvent ne pas avoir la forme d'une sphère parfaite mais une forme de polyèdre plus ou moins proche d'une sphère. De façon générale les « sphères de contact » sont ainsi des corps ayant !a dimensionnalité du réseau et une certaine extension spatiale. Dans l'exemple d'un réseau tridimensionnel, les sphères de contact sont ainsi des corps tridimensionnels. Dans d'autres exemples de réalisations, les sphères de contact pourront être des corps bidimensionnel comme des disques voir des corps unidimensionnel comme des segments. Une autre opération comprend la suppression du lien à rompre 350 de la pluralité de lien du modèle géologique 100 (étape 1520. SPR_LN). Comme l'illustre la figure 3d, ces opérations permettent la création de fissure et de fragments dans le modèle géologique à l'endroit des liens soumis à une force appliquée trop élevée. En se référant à présent plus précisément à la figure 4. une fois les liens à rompre supprimés et les sphères de contact créées. le procédé comporte la détermination d'une pluralité de contacts entre sphères 500 (étape 1600. DET_CTCT). La pluralité de contacts entre sphères 500 comporte au moins un contact entre sphères et peut en particulier être réduite à un contact entre sphères. Chaque contact entre sphères 500 associe deux sphères de contact 410, 420 parmi !a pluralité de sphères de contact 400. Un contact entre sphères 500 associe en particulier deux sphères de contact 410, 420 telles que les centres 430 desdites deux sphères de contact 410. 420 sont éloignés d'une distance d'éloignement d inférieure ou égale à la somme des rayons de référence 12 desdites deux sphères de contact 410, 420. Une fois la pluralité de contacts entre sphères 500 déterminée, pour chaque contact entre sphères de la pluralité de contacts entre sphères. on procède à la détermination d'une force d'interaction 1;,',ter entre les deux sphères de contact 410, 420 de la pluralité de sphères de contact 400 associées par ledit contact entre sphères / (étape 1700, DET_FINTER). A cette fin, on peut par exemple associer à chaque contact entre sphères de la pluralité de contacts entre sphères 500: - un plan de contact P, entre sphères, - une composante normale de la force d'interaction normale au plan de contact P,, et - une composante tangentielle de la force d'interaction f;i7;;,., comprise dans le plan de contact P,. La force d'interaction est alors la somme de la composante normale de la force d'interaction et de la composante tangentielle de la force d'interaction rat t Pi' - Plus précisément, la composante normale de la force d'interaction t,,, peut-être par exemple une composante élastique proportionnelle à une distance d'interpénétration M't' des deux sphères de contact 410. 420 associées par ledit contact entre sphères Une valeur de constante normale d'interaction peut être associée à ladite composante normale de la force d'interaction 1,7te, de s e a ce que la composante normale de la force d'interaction so t lt 17 tinr- En outre la composante tangentielle de la force d'interaction ;,V',. peut-être une composante de dissipation. Ainsi, la composante tangentielle de la force d'interaction py?. peut avoir une direction opposée à un déplacement relatif tangentiel des deux sphères de contact 410, 420 associées par ledit contact entre sphères r et un module proportionnel au module dudit déplacement relatif tangentiel f Ainsi la composante tangentielle de ia force d'interaction 7.7C7 peut- par exemple être donnée pa ' 'f. ,f,,7_ est une valeur de constante t 777 r? tangentielle d'interaction. Dans des variantes de l'invention, les composantes tangentielles et normales de la force d'interaction ftiliti,' pourront prendre des formes différentes. Pour chaque contact entre sphère de la pluralité de contacts entre sphères 500, le procédé peut encore comporter une étape de transmission de ladite force d'interaction fi77te7- à chaque lient de la pluralité de liens 300 tel que, parmi les deux noeuds de la pluralité de noeud 200 reliés ensemble par ledit lien 1, au moins un noeud est associé à une sphère de contact 400 associée audit contact entre sphères (étape 1800, TRSM_FINTER). De cette façon, les forces d'interaction générées à l'endroit des contacts entre sphères 500 peuvent être redistribuées dans le modèle géologique 100. L'étape de la détermination d'une position spatiale des noeuds de la pluralité de noeuds détaillée ci-avant (étape 1300. DET POS) peut avantageusement prendre en compte la pluralité de contacts entre sphères 500. Ainsi. pour chaque noeud i associé à une sphère de contact 400 associée à un contact entre sphères 500, la mise en oeuvre de l'algorithme d'intégration peut prendre en compte la force d'interaction entre les deux sphères de contact 410, 420 associées par ledit contact entre sphères 500, par exemple de la façon suivante. Pour chaque noeud i associé à une sphère de contact ni, l'accélération ai, au temps t + At. peut être déterminée par application de la seconde loi de Newton, en sommant les forces exercées sur ledit noeud i par les liens / reliant ledit noeud i à d'autres noeuds de la pluralité de noeuds 200 du modèle géologique et en sommant également les forces f't, exercées par les sphères de contact n associées à la sphère de contact in par un contact entre sphères 500. Ainsi par exemple l'accélération ai peut prendre la forme m 1 f" + Y,' f=e,- où m est la masse du noeud I. où la première somme est effectuée sur l'ensemble des noeuds j de la pluralité de noeuds 200 reliés au noeud i par un lien de la pluralité de liens 300, où fil est la force exercée par le noeud j sur le noeud où la deuxième somme est effectuée sur l'ensemble des sphères de contact in de la pluralité de sphères de contact 400 telles que la sphères de contact n est associé à la sphère de contact in associée au noeud t par un contact entre sphères 500 et oùest la force d'interaction entre la sphère de contact n et la sphère de contact in assoc ee au noeud 1. Les étapes du procédé décrites ci-dessus peuvent être réitérées au moins une fois, hormis l'étape de définition du modèle géologique. Avantageusement, les étapes du procédé décrites ci-dessus, exceptée l'étape de définition du modèle géologique, peuvent être réitérées tant qu'une condition de test n'est pas remplie, par exemple tant qu'un déplacement de la position spatiale d'au moins un noeud de la pluralité de noeud entre une itération n et une itération n + 1 du procédé est supérieur à un seuil de déplacement spatial prédéfini Le procédé fourni alors en sortie un modèle géologique endommagé. Par « modèle géologique endommagé » on entend ici en particulier un modèle géologique dans lequel les positions spatiales des noeuds et les liens sont modifiés par rapport au modèle géologique initial. La figure 4 représente un exemple de dispositif de détermination d'un endommagement d'un modèle géologique dans un mode de réalisation de l'invention. Dans ce mode de réalisation, le dispositif comporte un ordinateur 400, comprenant une mémoire 404 pour stocker des instructions permettant la mise en oeuvre du procédé, les données comme ou caractéristiques détaillées ci-avant des noeuds, des liens, des sphères de contacts, des contacts entre sphères, par exemple les positions, vitesses et accélérations des noeuds, et des données temporaires pour réaliser les différentes étapes du procédé tel que décrit précédemment. L'ordinateur comporte en outre un circuit de détermination d'un endommagement d'un modèle géologique. Ce circuit peut être, par exemple : - un processeur apte à interpréter des instructions sous la forme programme informatique, ou - une carte électronique dont les étapes du procédé de l'invention sont décrites dans le silicium, ou encore - une puce électronique programmable comme une puce FPGA (pour Field-Programmable Gate Array » en anglais). Cet ordinateur comporte une interface d'entrée 402 pour la réception de données comme cité précédemment. et une interface de sortie 403 pour la fourniture de données. Enfin, l'ordinateur comporte, pour permettre une interaction aisée avec un utilisateur. un écran 405. une souris et un davier 406. L'interface de sortie 406 peut être confondue avec une sortie vidéo permettant un affichage sur l'écran 401. Par ailleurs, le schéma fonctionnel présenté sur la figure 2 est un exemple typique 15 d'un programme dont certaines instructions peuvent être réalisées auprès de l'équipement décrit. A ce titre, la figure 2 peut correspondre à l'organigramme de l'algorithme général d'un programme informatique au sens de l'invention. Bien entendu, la présente invention ne se limite pas aux formes de réalisation décrites ci-avant à titre d'exemples ; elle s'étend à d'autres variantes. 20 D'autres réalisations sont possibles. REVENDICATIONS1. Procédé de détermination d'un endommagement d'un modèle géologique (100), ledit modèle géologique comportant une pluralité de noeuds (200) ainsi qu'une pluralité de liens (300). chaque lien de la pluralité de liens reliant ensemble deux noeuds de la pluralité de noeud, le procédé comprenant les étapes, - pour chaque lien de la pluralité de liens : * détermination (1400) d'une force appliquée audit lien, * si ladite force est supérieure (1500) à une valeur seuil de force - pour chaque noeud j des deux noeuds de la pluralité de noeud (200) reliés ensemble par ledit lien /, création (1510) d'une sphère de contact (400) associée audit noeud j, - suppression (1520) dudit lien i - détermination (1600) d'une pluralité de contacts entre sphères (500) comportant au moins un contact entre sphères, chaque contact entre sphères de la pluralité de contacts entre sphères associant deux sphères de contact (410, 420) parmi la pluralité de sphères de contact (400). 2. Procédé selon la revendication 1, dans lequel, pour chaque contact entre sphères i de la pluralité de contacts entre sphères (500), - détermination d'une force d'interaction entre les deux sphères de contact (410, 420) associées par ledit contact entre sphères 1, - pour chaque lien I de la pluralité de liens (300) tel que, parmi les deux noeuds de la pluralité de noeud reliés ensemble par ledit lien 1, au moins un noeud est associé à une sphère de contact (400) associée audit contact entre sphères 1, transmission de ladite force d'interaction fii'7ter audit lien I. 3. Procédé selon la revendication 2, dans lequel, à chaque contact entre sphères t de la pluralité de contacts entre sphères (500) est associé - un plan de contact P, entre sphères,N - une composante normale de la force d'interaction 117tCr normale au plan de contact et - une composante tangentielle de la force d'interaction ;;;;,,,r, comprise dans le plan de contact P. et dans lequel la composante normale de la force d'interaction est une composante proportionnelle à une distance d'interpénétration Dmt,,, des deux sphères de contact (410.420) associées par ledit contact entre sphères i, et - la composante tangentielle de la force d interaction r est une nte composante de direction opposée à un déplacement relatif tangentiel des deux sphères de contact (410, 420) associées par ledit contact entre sphères 1, et de module proportionnel audit déplacement relatif tangentiel 4. Procédé selon l'une des revendications précédentes, dans lequel une valeur de raideur k', une valeur seuil de force J! et une valeur de longueur de référence sont associées à chaque lien ï de la pluralité de liens (200), lesdites valeurs de raideur ki et valeurs seuil de force / étant déterminées à partir de valeurs de module élastique L' et de valeurs de seuil de résistance à la traction 0- d'au moins une phase géologique modélisée par le modèle géologique (100). 5. Procédé selon l'une des revendications précédentes, dans lequel les liens de la pluralité de liens (200) et les noeuds de la pluralité de noeuds (300) sont disposés en un réseau (120) régulier désordonné, de préférence en un réseau (120) tel qu'une distribution statistique de valeurs de longueur de référence 1 ou une distribution statistique de valeurs de raideur k, associées aux liens de la pluralité de liens (300), est une distribution statistique unimodale. 6. Procédé selon une des revendications précédentes, comprenant en outre les étapes, - application 1200) d'une pluralité de conditions aux limites une pluralité de noeuds de limites (240) parmi la pluralité de noeuds (200), - pour chaque noeud i de la pluralité de noeuds (200), détermination d'une position spatiale ri dudit noeud i. 7. Procédé selon l'une des revendications précédentes. dans lequel la détermination d'une position spatiale ri pour chaque noeud de la pluralité de noeuds (200) comprend soit une minimisation d'une énergie potentielle E7, du modèle géologique (100), soit une intégration d'une pluralité d'équations du mouvement associées aux noeuds (200). 8. Procédé selon l'une des revendications précédentes, dans lequel - une première pluralité de noeuds (210) parmi la pluralité de noeuds (200) modèle géologique (100) est associée à une première phase géologique (1 a) modélisée par le modèle géologique, - une deuxième pluralité de noeuds (220) parmi la pluralité de noeuds du modèle géologique est associée à une deuxième phase géologique (lb) modélisée par le modèle géologique, - une pluralité de liens de première phase géologique (310) parmi a pluralité de liens (300) du modèle géologique est constituée de liens de la pluralité de liens reliant ensemble deux noeuds de la première pluralité de noeuds (210), à chaque lien parmi la pluralité de liens de première phase géologique est associé une valeur de raideur de première phase géologique et une valeur seuil de force de première phase géologique f» déterminées à partir d'une valeur de module élastique de ta première phase géologique E`t modélisée par le modèle géologique et d'une valeur de seuil de résistance a la traction de la première phase géologique 0-In modélisée par le modèlegéologique, - une pluralité de liens de deuxième phase géologique (320) parmi la pluralité de liens (300) du modèle géologique est constituée de liens de la pluralité de liens reliant ensemble deux noeuds de la deuxième pluralité de noeuds (220), à chaque lien parmi la pluralité de liens de deuxième phase géologique est associé une valeur de raideur de deuxième phase géologique kl) et une valeur seuil de force de deuxième phase géologique déterminées à partir de valeurs de module élastique et de valeurs de seuil de résistance à la traction o de la deuxième phase géologique modélisée par le modèle géologique et - une pluralité de liens d'interface (330) parmi la pluralité de liens du modèle géologique est constituée de liens de la pluralité de liens (300) reliant ensemble un noeud de la première pluralité de noeuds (210) et un noeud de la deuxième pluralité de noeuds (220), à chaque lien parmi pluralité de liens d'interface est associé une valeur de raideur d'interfaceii71-cïr,yrc et une valeur seuil de force d'interface, -u?.Yr°a f et dans lequel une valeur de raideur de première phase géologique est différente d'une valeur de raideur de deuxième phase géologique et/ou une valeur de seuil de force de première phase géologique C est différente d'une valeur 20 de seuil de force de deuxième phase géologique f72. 9. Dispositif de détermination d'un endommagement d'un modèle géologique (100), ledit modèle géologique comportant une pluralité de noeuds (200) ainsi qu'une pluralité de liens (300). chaque lien de la pluralité de liens reliant 25 ensemble deux noeuds de la pluralité de noeud, le dispositif comportant un circuit (404) adapté pour permettre, - pour chaque lien de la pluralité de liens : a détermination d'une force appliquée audit lien, * si ladite force est supérieure à une valeur seuil de force :- pour chaque noeud des deux noeuds de la pluralité de noeud reliés ensemble par ledit lien. la création d'une sphère de contact associée audit noeud, et - la suppression dudit lien, - la détermination d'une pluralité de contacts entre sphères comportant au moins un contact entre sphères. ledit contact entre sphères associant deux sphères de contact parmi la pluralité de sphères de contact. 10. Produit programme informatique comportant des instructions pour la mise en oeuvre, lorsque ce programme est exécuté par un processeur. de la détermination d'un endommagement d'un modèle géologique, ledit modèle géologique comportant une pluralité de noeuds ainsi qu'une pluralité de liens, chaque lien de la pluralité de liens reliant ensemble deux noeuds de la pluralité de noeud, comportant les étapes : - pour chaque lien de la pluralité de liens : * détermination d'une force appliquée audit lien, ladite force est supérieure à une valeur seuil de force : - pour chaque noeud des deux noeuds de la pluralité de noeud relies ensemble par ledit lien, création d'une sphère de contact associée audit noeud, et - suppression dudit lien. - détermination d'une pluralité de contacts entre sphères comportant au moins un contact entre sphères, ledit contact entre sphères associant deux sphères de contact parmi la pluralité de sphères de contact.