地下储氢衬砌岩洞的建模仿真

在本篇博文中，客座作者 Qinghua Lei 探讨了一种用于分析嵌入裂隙岩体中的衬砌岩洞（LRC）在氢气循环加压条件下性能表现的仿真框架。

地下储氢正成为全球能源转型的关键一环。衬砌岩洞（Lined rock caverns ,LRC）提供了一种灵活且不受地域限制的解决方案，但其安全性取决于氢气、结构衬砌与裂隙岩体之间复杂的相互作用。要深入理解这些多物理场耦合效应，离不开先进的数值仿真技术。本文，我们将展示如何利用 COMSOL Multiphysics^® 软件对氢气加压过程中 LRC 的力学性能进行仿真分析。

为何使用 COMSOL Multiphysics^® 进行 LRC 储氢建模？

对储氢衬砌岩洞进行仿真面临着巨大的计算挑战。该系统涉及多种相互作用的介质（氢气、钢材、钢筋混凝土以及裂隙岩体），每种介质都遵循不同的物理规律。周围岩体中存在大量天然裂隙，会产生显著的几何不连续性和本构非线性，极大地影响系统的整体响应。裂隙岩体的变形与岩洞结构紧密交互，决定应力的重新分布和系统的整体稳定性。在高压氢气环境下，钢衬砌可能会发生氢脆现象，而混凝土衬砌则可能产生裂缝，改变应力的分布，进而影响钢衬砌的变形。这些紧密耦合的过程跨越了多个空间和时间尺度，使得对其进行精准分析成为一个极具挑战性的多物理场问题。

COMSOL Multiphysics^® 软件非常适合求解此类多物理场问题，因其具备以下卓越功能：

• 同步求解完全耦合的多物理场方程，实现力学、水力及传输过程之间的直接交互；
• 将模型参数定义为其他场变量的函数，实现间接耦合（例如应力相关的材料属性或损伤驱动的刚度演化）；
• 显式表征岩体中的离散裂隙，并解析复杂裂隙网络中的非线性流体-力学过程；
• 在统一框架内处理多种相互作用的材料，内置丰富的本构模型库并支持灵活的用户自定义材料公式；
• 可自定义控制方程和适配各类应用场景的本构模型，支持将氢脆模型等高级公式直接引入分析中。

接下来，我们将介绍在 COMSOL Multiphysics^® 中搭建用于 LRC 分析的数值模型的流程，并给出一则仿真示例。

仿真流程概述

使用 COMSOL Multiphysics^® 进行仿真时，包含三大核心步骤：

1. 构建几何模型并划分网格
2. 定义参数、材料属性、边界条件等
3. 求解计算

下面我们详细介绍各步骤。

几何模型和网格

首先，采用多尺度建模策略构建模型几何结构。在宏观尺度（图 1a）上，定义了一个二维区域来表示嵌入裂隙岩体中的衬砌岩洞的水平截面。岩洞的几何形状（包括钢衬砌、钢筋混凝土层和周围的喷射混凝土）可以直接在 COMSOL Multiphysics^® 中搭建，也可以从 CAD 软件导入。离散裂隙网络在岩体计算域内以线段形式进行几何表征，既可以在 COMSOL^® 中生成，也可以借助 MATLAB^® 或 CAD 软件等外部工具生成后导出为 DXF™ 格式文件，直接导入 COMSOL^®。在微观尺度上，则构建了专门的钢衬砌模型，以解析氢气的扩散和氢脆过程。这两个模型在钢-混凝土界面处通过位移协调条件进行耦合，从而确保结构变形与材料劣化之间能够保持一致的相互作用。

完成几何模型的定义或导入后，宏观区域采用非结构化三角形有限元网格进行离散化，该网格基于 Delaunay 三角剖分生成（图 2）。在岩洞边界附近和裂隙交叉处进行了网格加密，以准确解析应力集中和损伤演化。为了表征天然裂隙，在相邻有限元之间引入了节理单元（joint elements），从而能够对岩体中的非线性裂隙滑移和张开进行显式分析。为了保持计算效率，在宏观模型中，钢衬砌和喷射混凝土层等薄壁结构构件采用界面单元（interface elements）进行建模。然而，在微观模型中，将显式地表征钢衬砌的厚度（图 1b），并使用实体单元进行离散化，以解析氢气在厚度方向上的扩散和氢脆过程。这种多尺度离散化策略使得在高效捕捉全局结构响应的同时，能够在局部解析钢衬砌内部的材料劣化机制。

图 1. 多尺度模型的结构设计与边界条件示意图，包含：(a) 表征位于裂隙岩体中的衬砌岩洞（LRC）的宏观模型；(b) 用于分析钢衬砌在循环内压、边界位移约束作用下响应行为及氢气扩散过程的微观模型。图中标注的尺寸与比例沿用了本团队前期研究采用的模型参数（参考文献 1），此处展示仅为说明建模策略。实际的模型设置可根据具体应用场景与场地条件调整。

图 2. 衬砌岩洞（LRC）模型的网格离散化（参考文献 2）。钢筋混凝土层采用结构化四边形单元进行离散化，而岩石基质采用非结构化三角形单元进行离散化。钢衬砌和喷射混凝土采用线单元进行离散化，裂隙则采用节理单元进行离散化。

接口和耦合

我们使用 COMSOL Multiphysics^® 中的 固体力学 接口，模拟裂隙岩体、混凝土衬砌和钢构件在内部氢气加压作用下的变形。岩石基质采用标量损伤公式进行建模，以捕捉裂缝的萌生与扩展；离散裂隙则作为节理单元引入，并赋予非线性法向和剪切本构关系，以表征裂隙的滑移与张开。

钢筋混凝土衬砌采用了 COMSOL^® 中内置的 Mazars 损伤模型 进行模拟，该模型能够有效捕捉拉伸开裂和刚度的渐进式退化。钢筋的影响通过提高弹性模量和残余强度参数来体现，这些参数与所采用的混凝土等级保持一致。

在宏观模型中，由于钢衬砌几何形状较薄，因此采用界面单元进行表征。其力学行为遵循带有指数硬化特征的弹塑性本构规律。为了捕捉氢脆效应，我们构建了一个独立的微观模型，使用实体单元显式地表示钢衬砌的厚度。在该模型中，我们利用 COMSOL^® 中的 “物理场开发器” 开发了自定义接口，以模拟氢气在衬砌中的扩散，并考虑了其与固体力学的耦合。通过将屈服应力定义为随着氢浓度变化的变量，钢材的力学行为与氢浓度实现了耦合，从而能够分析由氢引起的材料强度退化。

若考虑裂隙岩体中的孔隙压力效应，可使用地下水流模块中的 达西定律 接口，模拟流体在岩石基质和裂隙中的流动。此时可启用 多孔弹性 耦合节点，以实现力学变形与孔隙压力演化之间的直接耦合。在这种情况下，渗透率或裂隙开度等水力属性可被定义为随应力变化的变量，从而实现与我们前期建模框架一致的间接流体-力学耦合。

针对岩石基质、裂隙、混凝土和钢材，我们分别定义了其材料属性和本构方程。直接多物理场耦合可确保变形与氢气传输（以及激活时的孔隙压力）之间的一致性交互；同时，通过将模型参数定义为随应力、损伤或氢浓度等演化场变量变化的函数，可引入额外的间接耦合。在力学边界条件方面，外部岩石边界施加了 原位 应力，岩洞内壁则施加了内部氢气压力。

求解计算

本次分析分为两个阶段进行：第一阶段，采用分级加载方案，使宏观模型在给定的原位应力下达到平衡状态。第二阶段，在岩洞边界施加氢气压力（单调递增或循环加载），以模拟储氢库的实际运行工况。随后，将宏观模型计算得出的位移场施加到微观钢衬砌模型上，并在该模型中以随时间变化的方式（瞬态）求解氢气扩散及与浓度相关的力学退化问题。此外，分析中采用了非线性求解方案，以准确解析裂隙的重新活化、损伤演化以及塑性变形问题。

LRC 仿真示例

我们采用多尺度模型对嵌入裂隙岩体中的衬砌岩洞（LRC）在氢气循环加压下的行为展开模拟（参考文献1）。宏观模型的模拟结果显示，混凝土衬砌呈现出循环径向位移，且混凝土与周围岩体中均出现了渐进式的损伤演化（图 3）。损伤主要集中在受拉区域以及裂隙交叉处附近（参考文献 1-2），这表明岩体中裂隙分布对 LRC 的结构响应具有显著控制作用（图 4）。

微观模型捕捉到了钢衬砌内部的氢气扩散与氢脆现象（图 5）（参考文献 1）。氢浓度从内壁表面开始逐渐升高，并随着加载循环不断演化，导致局部强度发生退化，且该退化与应力集中区域高度相关。这些结果充分证明了在不同尺度下，循环加压、裂隙重新活化、应力重新分布以及氢致损伤之间存在着复杂的耦合作用。

图 3. 仿真结果（参考文献 1）展示了在多次氢气加压循环（循环周期 T0 = 24 h）下：(a) 混凝土衬砌的径向位移分布；以及 (b) 混凝土及周围岩体的损伤演化的分布及发展规律。

图 4. 仿真结果（参考文献 1）展示了在循环氢气加压（循环周期 T0 = 24 h）条件下，LRC 周边的损伤分布以及局部最大主应力分布。

图 5. 仿真结果（参考文献 1）展示了不同加载阶段下钢衬的 (a) 氢浓度空间分布、(b) 最大主应力变化，以及 (c) 强度退化。

此外，该框架已进行了扩展，将裂隙岩体中的流体-力学耦合（参考文献 3）以及岩石随时间变化的蠕变（参考文献 4）纳入其中。通过这些扩展，我们能够评估围岩中的流体压力扩散和黏弹性变形对 LRC 长期性能的影响。

参考文献

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.

关于作者

Qinghua Lei 现任瑞典乌普萨拉大学（Uppsala University）副教授。他先后于 2009 年和 2012 年在中国同济大学获得土木工程专业学士和硕士学位，并于 2016 年在英国帝国理工学院获得岩石力学专业博士学位。2016 年至 2018 年，他在帝国理工学院从事流体力学方向博士后研究；随后于 2018 年至 2023 年在瑞士苏黎世联邦理工学院（ETH Zurich）担任工程地质学高级研究员和讲师。他曾荣获 2025 年欧洲研究委员会基金（ERC Consolidator Grant）、2024 年 Chin-Fu Tsang 耦合过程奖（Chin-Fu Tsang Coupled Processes Award）、2019 年罗哈奖（Rocha Medal）、2016 年 NGW 库克博士论文奖（NGW Cook PhD Dissertation Award）以及 2015 年岩石力学研究奖（Rock Mechanics Research Award）。此外，他还是欧洲青年科学院（Young Academy of Europe）院士，其研究方向包括岩石力学、水文地质学、地球物理学、地质灾害以及岩土工程。
 
MATLAB 是 MathWorks, Inc. 的注册商标。Autodesk、Autodesk 徽标、AutoCAD 和 DXF 是 Autodesk, Inc. 和/或其子公司和/或关联公司在美国和/或其他国家的注册商标或商标。