多孔弹性介质内部不连续性的偏移界面方法

《Computers and Geotechnics》：A shifted interface approach for internal discontinuities in poroelastic media

　　

摘要
含有裂缝、断层或内部不连续性的多孔介质广泛存在于地下工程力学、生物力学和材料科学中。由于压力场和位移场通过 Biot 方程紧密耦合，需要稳定的混合格式，因此耦合流固响应的数值模拟本身就极具挑战性。当存在裂缝时，这些困难会加剧，因为标准的网格拟合方法需要耗费大量人力和成本的非结构化网格划分，而非拟合方法通常需要切割单元积分、富集函数或额外的稳定化措施。在本工作中，我们采用了一种替代方法，即将偏移界面法 (SIM) 适配到具有嵌入界面的耦合瞬态多孔弹性中。该方法通过局部展开将界面条件转移，用代理近似替代真实裂缝。统一推导得出了水力传输和机械牵引耦合的偏移形式。我们广泛比较了两种施加策略：弱（积分）施加和强（逐点）施加。四个几何复杂度递增的测试案例（偏移网格对齐、边界相交倾斜、嵌入倾斜以及多裂缝配置）验证了该公式。在裂缝尖端以外，界面残差以 ??(h) 的速度收敛；在尖端附近，局部后处理伪影降低了全局收敛率，但当排除小的尖端区域时，可恢复一阶收敛。一个包含四个同时嵌入的不同几何形状和界面特性的裂缝的多裂缝演示证实了该框架的实际适用性。这些结果支持偏移界面法作为在非体拟合网格上对具有复杂几何嵌入界面的多孔弹性裂缝建模的实用框架。

引言
贯穿裂缝、断层和其他内部不连续性的多孔材料广泛存在于工程和自然科学中。在地下工程力学中，断层和断裂岩体在储层和盖层的流固耦合响应中起着关键作用 (Caine et al., 1996, Bense et al., 2013, Berre et al., 2019)；在生物力学中，软骨、骨骼和软组织是流体饱和的多孔介质，其中表面裂缝和界面分层严重影响承载和修复 (Cowin, 1999, Ateshian, 2009)；在材料科学中，混凝土的耐久性 (Ulm et al., 2004) 以及水凝胶 (Bouklas et al., 2015) 和金属泡沫 (Iaccarino et al., 2026) 的断裂同样受流体传输与断裂之间相互作用的控制。在所有这些设置中，内部不连续性对流体流动和机械应力传递都施加了控制性影响，并可能同时通过滑移、张开或弹性顺从来容纳位移不连续性。
对于具有广泛社会意义的应用，捕捉这种耦合的流固行为尤为迫切，包括二氧化碳封存 (Rutqvist, 2012, Pan et al., 2016)、地热能开采和诱发地震风险评估 (Ellsworth, 2013, Foulger et al., 2018, Kivi et al., 2023) 以及核废料处置库的长期安全性 (Tsang et al., 2005)。因此，能够表征本体多孔弹性响应与嵌入界面物理相互作用的高效数值方法，而无需网格拟合于可能复杂的界面几何形状，具有相当大的实用价值。
从数值模拟的角度来看，建模不连续性具有挑战性。有两种常见策略。第一种是将不连续性解析为具有极低（或高度各向异性）材料特性的薄 2D/3D 区域。虽然概念上简单且与标准多孔弹性求解器兼容，但随着屏障厚度减小，该方法在数值上性能下降：极端纵横比和小单元导致单元畸变、病态以及在瞬态模拟中限制时间步长。此外，对于复杂几何形状，网格生成变得脆弱，甚至不可能 (Starnoni et al., 2021)。第二种策略是降低屏障的维度，并将其建模为内部界面（即 (???1)D 表面），在拟合于界面的网格上赋予适当的本构定律。
这种方法还使得在不直接划分薄界面层或在复杂几何形状中承受大单元畸变的情况下对断裂进行建模成为可能。事实上，当内部界面几何形状复杂时，非拟合（或浸入式或嵌入式）方法变得具有吸引力，因为界面几何形状独立于计算网格：网格无需拟合于界面，而是将界面条件施加于可能切割或靠近真实界面的背景网格上。XFEM 引入富集函数以在不网格对齐的情况下捕捉嵌入界面处的场跳跃 (Mo?s et al., 1999, Flemisch et al., 2016)。其灵活性是以代价换来的：需要对被界面切割的单元进行特殊积分、混合校正、线性代数病态问题（特别是对于相交或终止特征）以及额外的实现复杂性。CutFEM（及相关非拟合 Nitsche 方法）则在背景网格上保持标准有限元空间，并通过稳定化（例如，幽灵惩罚）弱施加界面条件，以正则化小的切割单元 (Hansbo and Hansbo, 2002, Burman, 2010, Burman et al., 2015)。最近的工作将 CutFEM 扩展到了具有复杂外部边界的准静态多孔弹性 (Berre et al., 2025)、自由流体流与多孔弹性介质的耦合 (Ager et al., 2019)，以及具有材料界面的非均匀多孔弹性域（场保持连续但性质跳跃的弱不连续性）(Liu et al., 2025)。最近，提出了调和方法有限元 (TFEM)，用于包含几乎退化或零测度单元且在完全坍缩极限下的网格；它对应于砌缝 (Quiriny et al., 2025)。然而，其耦合是几何诱导的，而不是通过显式的断裂本构定律引入的，因此它不能直接提供此处考虑的通用水力学和机械界面关系。尽管如此，这些方法尚未关注具有位移和/或压力跳跃的内部断裂，而这是本贡献的重点。相场模型提供了另一条途径，将离散断裂正则化为连续损伤场并与多孔弹性耦合 (Mikeli? et al., 2015, Wilson and Landis, 2016, Chukwudozie et al., 2019)。虽然相场方法自然地处理复杂的裂缝拓扑和扩展，但它们引入了一个全局辅助场，必须在足够精细的网格上解析以捕捉正则化长度，导致显著的 计算成本。
我们采用了一种替代方法，即偏移界面法 (SIM)，它源于偏移边界法 (SBM) (Main and Scovazzi, 2018, Scovazzi et al., 2018)。该方法的关键思想是用拟合于背景网格（例如，单元面的并集）的附近代理界面来替代真实内部界面，并通过局部展开将内部界面条件在位置和值上转移到真实界面条件 (Li et al., 2020, Atallah et al., 2022)。最近的发展包括用于内聚裂缝模拟的偏移断裂法 (SFM) (Li et al., 2020, Li et al., 2023)、用于机械接触的偏移公式 (Li et al., 2025)，以及扩展到耦合瞬态热弹性 (Li et al., 2024)。这些发展共同表明，SIM 为具有复杂几何嵌入界面的复杂多物理场问题提供了一个有前景的框架，同时保持数值基础设施接近标准网格拟合有限元，并避免了 CutFEM 和 XFEM 典型的切割单元积分、富集和稳定化机制。此外，这种方法是非侵入式的，意味着它仅需要在变分形式中增加一些附加项，而不引入新的自由度，也不改变标准 FEM 数据结构或需要特设的稳定化项。该方法在逼近能力和收敛性方面已达到理论成熟度，为我们在此建立和扩展奠定了基础。
在拟合于不连续性的网格上，刚度型机械界面定律与由压力跳跃驱动的横向水力导水率的组合已通过零厚度流固耦合界面单元 (Segura and Carol, 2008, Cerfontaine et al., 2015, Cui et al., 2019) 得到充分建立，通常节理定律以有效界面应力提出，通道型界面流动模型携带储存和孔径相关的导水率。据我们所知，本工作是首次将此类组合界面闭合带入非体拟合网格；此处采用的闭合是总牵引弹簧定律配对具有恒定导水率的锐利零储存膜，这是所考虑半透屏障的适当理想化。现有的 Biot 系统非拟合处理方法处理外部域边界或完美粘合的材料界面；嵌入流固耦合断裂公式将内聚或摩擦力学与通道流或漏失水力配对 (Damirchi et al., 2022)；而本构界面定律的偏移处理迄今仍纯粹是机械性的 (Li et al., 2023, Li et al., 2025)。
我们将 SIM 扩展到具有嵌入界面的耦合瞬态多孔弹性。该公式考虑了水力传输（通过将法向通量与压力跳跃相关的 Robin 型条件）和机械响应（通过将牵引与法向和切向位移跳跃相关的牵引关系）。基于将通用界面通量分解为真实法向和切向失配分量的统一推导，得出了两个场的偏移形式。推导出的公式适用于原始场方法（例如，非连续和连续 Galerkin 法），但本文中我们在沿代理界面分裂连接性的网格上应用连续 Galerkin 法，结合了界面处远处拟合单元的效率以及表示跨界面跳跃的灵活性。
本文结构如下。在第 2 节中，我们介绍了控制偏微分方程、相关的外部 and 内部边界条件，以及计算域上的变分公式。在第 3 节中，提供了给定 PDE 的偏移界面形式的推导，包括界面本构闭合和两种施加策略（弱和强）。第 4 节描述了数值实现：有限元空间、两种施加策略的半离散问题，以及用于评估真实裂缝几何上界面残差的后处理流程。第 5 节展示了四个复杂度递增测试案例的数值结果：偏移网格对齐裂缝、边界相交倾斜裂缝、嵌入倾斜裂缝，以及具有四个同时嵌入的不同几何形状和界面特性的多裂缝配置。

章节片段
问题描述
设 ?????? (??=1,2,3) 为多孔体，??????? 为嵌入的 Lipschitz 界面，其单位法向为 ????（方向从“?"到“＋"）。界面可在 ?? 内终止或延伸至 ????，使得 ??????? 要么连通，要么由不连通子域 ??± 组成。在任何情况下，对于 ??????? 上的场 ??，我们在 ???? 上定义单侧迹 ??± 以及跳跃/平均算子 [[??]]???+????，{{??}}?12?(??++???)。
偏移界面法
该方法的关键特征是真实裂缝 ???? 不与计算网格对齐。偏移界面法通过将界面条件从 ???? 转移到代理 ????? 来解决这个问题。以下，我们展示偏移界面法的推导过程（参见 (Li et al., 2023)）。
数值实现
我们现在规定本工作中使用的空间有限维空间，并澄清半离散仅指空间离散：在有限元基展开后，问题变为关于系数向量的时间连续系统（此处为 DAE）。本体场由连续 Galerkin 有限元近似。然而，代理界面 ????? 必须允许来自每一侧的独立迹，以表示压力和位移跳跃。这是通过...实现的
数值结果
所有数值实验均在平面应变假设下于 2D 中进行。变分形式方程 (55)、(58) 在 scikit-fem 内实现 (Gustafsson and McBain, 2020)。对于时间积分，我们使用了 solve_nivp (Riley, 2025)，采用自适应 SDIRK2 时间积分器。完整的工作流（组装、网格分裂、界面记录管理和后处理）在我们的 Python 代码 PoroSIM 中实现（见数据可用性声明）。梯度修正块（离散...

结论
本工作将偏移界面法适配到了具有嵌入界面的耦合瞬态多孔弹性中。从描述被内部不连续性贯穿的多孔介质中流动和变形的标准 Biot 方程出发，我们推导了一个统一框架，将水力传输（Robin 型通量 - 压力耦合）和机械响应（弹簧型牵引 - 位移耦合）从真实裂缝转移到拟合于背景网格的代理界面。

CRediT 作者贡献声明
David Michael Riley: 撰写 – 初稿，可视化，验证，软件，正式分析。Guglielmo Scovazzi: 撰写 – 审阅与编辑，监督，方法论。Ioannis Stefanou: 撰写 – 审阅与编辑，监督，方法论，资金获取，概念化。

利益冲突声明
作者声明他们没有已知的竞争财务利益或个人关系，这些利益或关系可能影响了本文中报告的工作。

致谢
由欧盟资助。

然而，文中表达的观点和意见仅代表作者，不一定反映欧盟或欧洲研究委员会执行机构的观点。欧盟及资助机构均不承担任何责任。本研究得到 ERC 资助项目 INJECT 的支持，编号 101087771，doi: 10.3030/101087771。
David Michael Riley|Guglielmo Scovazzi|Ioannis Stefanou
IMSIA (UMR 9219), 法国国家科学研究中心，法国电力公司，巴黎国立高等技术与科学学院，巴黎综合理工学院，帕莱索，91120，法国