Modéliser les cavités rocheuses revêtues pour le stockage souterrain d’hydrogène

Dans cet article de blog, Qinghua Lei, auteur invité, présente un cadre de modélisation permettant d’analyser la mise sous pression cyclique à l’hydrogène d’une LRC intégrée dans une masse rocheuse fracturée.

Le stockage souterrain de l’hydrogène devient essentiel pour la transition énergétique mondiale. Les cavités rocheuses revêtues (LRC, pour Lined Rock Caverns) offrent une solution flexible et adaptable sur le plan géographique, mais leur sécurité dépend d’interactions complexes entre l’hydrogène gazeux, les revêtements structurels et les masses rocheuses fracturées. La compréhension de ces effets couplés nécessite une modélisation numérique avancée. Dans cet article de blog, nous montrons comment le logiciel COMSOL Multiphysics^® peut être utilisé pour modéliser le comportement des LRC lors de l’injection d’hydrogène sous pression.

Pourquoi utiliser COMSOL Multiphysics^® pour modéliser le stockage d’hydrogène dans les LRC?

La modélisation de cavités rocheuses revêtues destinées au stockage de l’hydrogène pose d’importants défis de calcul. Le système implique plusieurs matériaux en interaction (hydrogène gazeux, acier, béton armé et roche fracturée), chacun régi par un comportement physique distinct. La masse rocheuse environnante contient de nombreuses fractures préexistantes, ce qui introduit des discontinuités géométriques marquées et des non-linéarités constitutives qui influencent fortement la réponse du système. La déformation de la masse rocheuse fracturée interagit étroitement avec la structure de la caverne, régissant la redistribution des contraintes et la stabilité globale du système. Sous l’effet d’une exposition à l’hydrogène sous haute pression, le revêtement en acier peut être sujet à une fragilisation, tandis que le revêtement en béton peut présenter des fissures qui redistribuent les contraintes et influencent la déformation du revêtement en acier. Ces processus étroitement couplés s’étendent sur plusieurs échelles spatiales et temporelles, faisant de l’analyse réaliste un problème multiphysique complexe.

COMSOL Multiphysics^® est particulièrement adapté à la résolution de ce type de problèmes multiphysiques grâce à ses capacités exceptionnelles pour:

• Résoudre simultanément des équations multiphysiques entièrement couplées, permettant des interactions directes entre les phénomènes mécanique, hydraulique et de transport
• Définir des paramètres du modèle en fonction d’autres variables de champ, permettant ainsi des couplages indirects tels que des propriétés matériaux fonction des contraintes ou une évolution de la rigidité induite par l’endommagement
• Représenter explicitement des fractures discrètes au sein des masses rocheuses et modéliser des processus hydromécaniques non linéaires au sein de réseaux de fractures complexes
• Gérer plusieurs matériaux en interaction au sein d’un cadre de travail unifié, grâce à une vaste bibliothèque de lois de comportement prédéfinies, ainsi qu’à des formulations de matériaux personnalisables et flexibles
• Implémenter des équations et des modèles constitutifs personnalisées, spécifiques à chaque application, permettant d’intégrer directement dans l’analyse des formulations avancées telles que des modèles de fragilisation par l’hydrogène

Nous abordons ci-dessous la procédure de création de modèles numériques dans COMSOL Multiphysics^® pour l’analyse des LRC et présentons un exemple de simulation.

Aperçu des étapes de modélisation

Lors de l’utilisation de COMSOL Multiphysics^®, la simulation comporte trois étapes clés:

1. Génération de la géométrie et du maillage
2. Paramètres de couplage et définition des propriétés matériaux, des conditions aux limites, etc.
3. Calcul de la solution

Examinons ces étapes plus en détail.

Géométrie et maillage

Tout d’abord, la géométrie du modèle est construite selon une stratégie multi-échelle. À grande échelle (figure 1a), un domaine bidimensionnel est défini pour représenter une coupe transversale horizontale de la cavité rocheuse revêtue, intégrée dans une masse rocheuse fracturée. La géométrie de la cavité, comprenant le revêtement en acier, la couche de béton armé et le béton projeté environnant, peut être créée directement dans COMSOL Multiphysics^® ou importée à partir d’un logiciel de CAO. Les réseaux de fractures discrets sont représentés géométriquement sous forme de segments de droite au sein du domaine rocheux et peuvent être générés dans COMSOL^® ou à l’aide d’outils externes tels que MATLAB^® ou un logiciel de CAO, puis exportés sous forme de fichiers DXF™ pour être importés directement dans COMSOL^®. À petite échelle, un modèle dédié au revêtement en acier est construit afin de modéliser les processus de diffusion de l’hydrogène et de fragilisation. Les deux modèles sont couplés à l’interface acier-béton en imposant la compatibilité des déplacements, ce qui permet une interaction cohérente entre la déformation structurelle et l’endommagement du matériau.

Une fois la géométrie définie ou importée, le domaine à grande échelle est discrétisé à l’aide d’un maillage non structuré d’éléments finis triangulaires généré par tessellation de Delaunay (figure 2). Un raffinement du maillage est appliqué à proximité de la frontière de la cavité et autour des intersections de fractures afin de modéliser avec précision les concentrations de contraintes et l’évolution de l’endommagement. Pour représenter les fractures naturelles, des éléments d’interface entre éléments finis voisins sont mis en oeuvre, ce qui permet une analyse explicite du glissement et de l’ouverture non linéaires des fractures au sein de la masse rocheuse. Les composants structurels minces, tels que le revêtement en acier et la couche de béton projeté, sont modélisés à l’aide d’éléments d’interface dans le modèle à grande échelle afin de préserver l’efficacité du calcul. Dans le modèle à petite échelle, en revanche, l’épaisseur du revêtement en acier est représentée explicitement (figure 1b) et discrétisée à l’aide d’éléments solides afin de modéliser la diffusion de l’hydrogène et la fragilisation à travers toute l’épaisseur. Cette stratégie de discrétisation multi-échelle permet de rendre compte efficacement de la réponse structurelle globale tout en modélisant localement les mécanismes de dégradation du matériau au sein du revêtement en acier.

Figure 1. Schéma conceptuel et conditions aux limites du modèle multi-échelle, comprenant (a) un modèle à grande échelle représentant une LRC située dans un massif rocheux fracturé et (b) un modèle à petite échelle décrivant la réponse du revêtement en acier soumis à une pressurisation interne cyclique, à des conditions aux limites de déplacement imposé, ainsi qu’à la diffusion de l’hydrogène. Les dimensions et les échelles indiquées suivent la configuration adoptée dans notre étude précédente (réf. 1) et sont présentées ici afin d’illustrer la stratégie de mise en place du modèle. La configuration réelle peut varier selon les applications et les conditions du site.

Figure 2. Discrétisation du maillage pour le modèle de LRC (réf. 2). La couche de béton armé est maillée à l’aide d’éléments quadrangles structurés, tandis que la matrice rocheuse est maillée à l’aide d’éléments triangulaires libres. Le revêtement en acier et le béton projeté sont discrétisés à l’aide d’éléments d’arêtes, tandis que les fractures le sont à l’aide d’éléments d’interface.

Interfaces et couplages

Nous utilisons l’interface Mécanique du solide de COMSOL Multiphysics^® pour simuler la déformation de la masse rocheuse fracturée, du revêtement en béton et des composants en acier sous l’effet de la pression interne exercée par l’hydrogène. La matrice rocheuse est modélisée à l’aide d’une formulation d’endommagement scalaire afin de rendre compte de l’amorçage et de la propagation des fissures, tandis que les fractures discrètes sont modélisées sous forme d’éléments d’interface dotés de relations constitutives non linéaires normales et de cisaillement pour représenter le glissement et l’ouverture des fractures.

Le revêtement en béton armé est modélisé à l’aide du modèle d’endommagement de Mazars, disponible dans COMSOL^®, qui prend en compte la fissuration en traction et la perte progressive de rigidité. L’influence de l’armature est prise en compte par le biais d’un module d’élasticité amélioré et de paramètres de résistance résiduelle, conformément à la classe de béton retenue.

Dans le modèle à grande échelle, le revêtement en acier est représenté à l’aide d’éléments d’interface en raison de sa faible épaisseur. Son comportement mécanique obéit à une loi constitutive élastoplastique avec écrouissage exponentiel. Afin de prendre en compte les effets de la fragilisation par l’hydrogène, un modèle distinct à petite échelle représente explicitement l’épaisseur du revêtement en acier à l’aide d’éléments solides. Dans ce modèle, une interface personnalisée a été développée à l’aide du Constructeur de physiques de COMSOL^® afin de simuler la diffusion de l’hydrogène à travers le revêtement tout en tenant compte du couplage avec la mécanique du solide. Le comportement mécanique de l’acier est couplé à la concentration en hydrogène en définissant la limite d’élasticité comme une variable dépendant de la concentration, ce qui permet d’analyser la perte de résistance induite par l’hydrogène.

Lorsque l’on tient compte des effets de la pression interstitielle dans la masse rocheuse fracturée, l’interface Loi de Darcy du module Subsurface Flow est ajoutée pour simuler l’écoulement fluide au sein de la matrice rocheuse et des fractures. Le couplage Poroélasticité peut être activé pour obtenir un couplage direct entre la déformation mécanique et l’évolution de la pression interstitielle. Dans ce cas, les propriétés hydrauliques telles que la perméabilité ou l’ouverture des fractures peuvent être définies comme des variables dépendant des contraintes, ce qui permet un couplage hydromécanique indirect compatible avec notre cadre de modélisation précédent.

Les propriétés matériaux et les équations constitutives sont définies séparément pour la matrice rocheuse, les fractures, le béton et l’acier. Un couplage multiphysique direct garantit une interaction cohérente entre la déformation et le transport d’hydrogène (ainsi que la pression interstitielle lorsqu’elle est activée), tandis que des couplages indirects supplémentaires sont introduits en définissant les paramètres du modèle comme des fonctions de variables de champ telles que la contrainte, l’endommagement ou la concentration en hydrogène. Les conditions aux limites mécaniques comprennent les contraintes in situ appliquées sur les frontières extérieures de la roche et la pression interne d’hydrogène appliquée sur la paroi de la cavité.

Calcul de la solution

L’analyse se déroule en deux étapes: lors de la première étape, le modèle à grande échelle est mis à l’équilibre statique, en appliquant progressivement les contraintes in situ à l’aide d’une procédure de ramping. Lors de la deuxième étape, une mise sous pression d’hydrogène est appliquée aux frontières de la cavité, de façon monotone ou cyclique, afin de simuler l’exploitation du site de stockage. Le champ de déplacement résultant du modèle à grande échelle est imposé au modèle en acier à petite échelle, où la diffusion de l’hydrogène et l’endommagement mécanique dépendant de la concentration sont résolus en temporel. Des méthodes de résolution non linéaires sont utilisées pour traiter la réactivation des fractures, l’évolution de l’endommagement et la plasticité.

Exemple de simulation d’une LRC

Nous appliquons le modèle multi-échelle pour simuler la mise sous pression cyclique à l’hydrogène d’une LRC encastrée dans une masse rocheuse fracturée (réf. 1). Les résultats de la modélisation à grande échelle montrent un déplacement radial cyclique du revêtement en béton et un développement progressif de l’endommagement tant dans le béton que dans la roche environnante (figure 3). Il se localise principalement dans les zones soumises à des contraintes de traction et à proximité des intersections de fractures (réf. 1-2), ce qui met en évidence l’influence déterminante de la répartition des fractures dans la roche sur la réponse structurelle de la LRC (figure 4).

Le modèle à petite échelle reproduit la diffusion de l’hydrogène et la fragilisation au sein du revêtement en acier (figure 5) (réf. 1). La concentration en hydrogène augmente depuis la surface interne et évolue au fil des cycles de charge, entraînant une dégradation locale de la résistance qui est en corrélation avec les zones de concentration de contraintes. Ces résultats mettent en évidence l’interaction couplée entre la mise sous pression cyclique, la réactivation des fissures, la redistribution des contraintes et l’endommagement induit par l’hydrogène à différentes échelles.

Figure 3. Résultats de simulation (réf. 1) illustrant la répartition et l’évolution (a) du déplacement radial dans le revêtement en béton et (b) de l’endommagement dans le béton et la roche environnante au cours de plusieurs cycles de mise sous pression à l’hydrogène (durée du cycle T0 = 24 h).

Figure 4. Résultats de simulation (réf. 1) illustrant la répartition de l’endommagement et la contrainte principale locale maximale à proximité de la LRC lors d’une mise sous pression cyclique à l’hydrogène gazeux (période de cycle T0 = 24 h).

Figure 5. Résultats de simulation (réf. 1) illustrant la répartition spatiale (a) de la concentration en hydrogène, (b) de la variation de la contrainte principale maximale et (c) de la perte de résistance du revêtement en acier à différents stades de mise en charge.

De plus, le modèle a été étendu pour prendre en compte le couplage hydromécanique dans la masse rocheuse fracturée (réf. 3) et le fluage de la roche en fonction du temps (réf. 4). Ces améliorations nous permettent d’évaluer dans quelle mesure la diffusion de la pression des fluides et la déformation viscoélastique de la masse rocheuse environnante influencent les performances à long terme de la LRC.

Références

1. C. Zhao et al., “Modelling lined rock caverns subject to hydrogen embrittlement and cyclic pressurisation in fractured rock masses,” International Journal of Hydrogen Energy, 2025; 152: 150027.
2. C. Zhao, Q. Lei, Z. Zhang, “Impact of fracture networks on the structural deformation of lined rock caverns under high internal gas pressure,” Underground Space, 2025; 21: 252-269.
3. C. Zhao, Z. Zhang, Q. Lei, “Coupled hydro-mechanical simulation of the interaction between adjacent lined rock caverns subject to internal gas pressurisation,” Geomechanics for Energy and the Environment, vol. 43, 2025: 100701.
4. C. Zhao et al., “Influence of rock creep on the performance of lined caverns under cyclic pressurization and hydrogen embrittlement,” International Journal of Rock Mechanics and Mining Sciences, vol. 199, 2026; 106401.

À propos de l’auteur

Qinghua Lei est maître de conférences à l’université d’Uppsala, en Suède. Il a obtenu sa licence (2009) et son master (2012) en génie civil à l’université de Tongji, en Chine, ainsi que son doctorat (2016) en mécanique des roches à l’Imperial College London, au Royaume-Uni. Il a occupé un poste de chercheur postdoctoral en mécanique des fluides à l’Imperial College London (2016-2018), puis de chercheur senior et de maître de conférences en géologie de l’ingénierie à l’ETH Zurich, en Suisse (2018-2023). Il est lauréat de la bourse ERC Consolidator Grant 2025, du prix Chin-Fu Tsang Coupled Processes Award 2024, de la médaille Rocha 2019, du prix NGW Cook PhD Dissertation Award 2016 et du prix Rock Mechanics Research Award 2015. Il est membre de la Young Academy of Europe. Ses domaines de recherche comprennent la mécanique des roches, l’hydrogéologie, la géophysique, les risques naturels et le génie géotechnique.
 

MATLAB is a registered trademark of The MathWorks, Inc. Autodesk, the Autodesk logo, AutoCAD, and DXF are registered trademarks or trademarks of Autodesk, Inc., and/or its subsidiaries and/or affiliates in the USA and/or other countries.