用于分段均质多孔介质水力-力学耦合分析的界面增强型广义有限元方法

《Computers and Geotechnics》：Interface-enriched generalized finite element method for the coupled hydro-mechanical analysis of piecewise homogeneous porous media

　　普里亚·皮尔莫拉迪（Pouriya Pirmoradi）| 亚历杭德罗·M·阿拉贡（Alejandro M. Aragón）| 帕亚姆·普尔斯索尔胡伊（Payam Poorsolhjouy）| 阿克·S.J. 苏伊克（Akke S.J. Suiker）
埃因霍温理工大学建筑环境系，De Groene Loper 6，5612 AP 埃因霍温，荷兰

摘要
本文提出了一种增强的界面通用有限元方法（IGFEM），用于耦合水力-力学分析，这种分析针对的是由不同但完美结合的物质相组成的可变形、饱和的多孔介质。使用IGFEM对动量平衡方程和储存方程进行空间离散化，然后通过广义Newmark方法对这些方程进行时间离散化。这导致了一个完全耦合的非线性方程组，该方程组采用整体更新方案进行迭代求解。IGFEM公式在准确捕捉固体相位移场和流体相压力场在材料界面处的弱不连续性方面非常有效，因为它直接在这些界面上放置了增强节点。多个数值示例表明，所提出的IGFEM公式不仅与使用共形网格的标准FEM具有相同的精度，而且性能优于扩展/通用有限元方法（XFEM/GFEM）。此外，它能够在不需要网格对齐的情况下准确捕捉复杂的非平面界面，凸显了该方法在处理具有几何形状复杂边界的多孔介质的水力-力学分析中的灵活性和鲁棒性。总体而言，IGFEM为涉及材料界面的瞬态耦合问题提供了一种高精度且高效的方法。

1. 引言
多孔介质的水力-力学行为受流体和固体相之间的相互作用控制。准确预测这种行为需要建模流体流动与固体基质变形之间的强耦合关系。描述固体骨架与孔隙流体之间准静态相互作用的基础方程最初由Biot（1941年）提出，后来他扩展了该框架以考虑动态波动传播（Biot，1956年；Biot，1962年）。经典混合理论是多孔介质中不同相建模的连续介质力学的基石，最初由Truesdell和Toupin（1960年）提出，该理论基于体积分数的概念。这一理论通过几项关键发展得到了进一步完善：Green和Naghdi（1969年）开发了一个符合热力学第二定律的热力学一致框架；Bowen（1980年、1982年）专注于为可压缩和不可压缩多孔介质制定严格的平衡定律和相互作用力；Whitaker（1977年）引入了体积平均技术以实现微观到宏观的放大；Hassanizadeh和Gray（1980年、1990年）将该框架扩展到多相流动和界面现象的非平衡热力学领域。对于涉及几何和非线性材料特性的非线性问题，Zienkiewicz等人（1980年、1982年、1990a年；Prévost，1980年、1982年）在有限元方法（FEM）框架内推进了迭代求解程序。基于这些开创性工作，已经开发了许多用于饱和多孔介质（Zienkiewicz等人，1990a年；Schrefler等人，1998年；Li等人，2004年；Zhang等人，2009年；Hajiabadi和Khoei，2019年；Scheperboer等人，2022年）和部分饱和多孔介质（Zienkiewicz等人，1990b年；Schrefler和Xiaoyong，1993年；Schrefler和Scotta，2001年；Sheng等人，2003年；Sheng等人，2008年；Hoteit和Firoozabadi，2008年；Wang等人，2009年；Khoei和Mohammadnejad，2011年；Khoei和Saeedmonir，2021年）的水力-力学分析的数值模型。De Boer（2012年）的综述提供了关于多孔介质的各种理论框架的全面概述。

多孔介质中的不连续性，如材料界面、空隙和裂缝，显著影响固体基质的变形行为和孔隙流体的流动模式。因此，在水力-力学分析中准确考虑这些特征至关重要。有限元方法通过使用与不连续性几何形状相匹配的网格来建模不连续多孔介质，使元素边缘与裂缝或界面位置对齐。这种网格适配方法允许结合粘聚性界面元素，从而使用标准有限元对多孔介质的完整区域进行建模，同时通过粘聚性元素显式捕获裂缝（Segura和Carol，2008年；Khoei等人，2011年；Sarris和Papanastasiou，2012年；Carrier和Granet，2012年）。然而，当不连续性演变时（例如在裂缝扩展过程中），这种方法需要网格自适应性（Karimi-Fard和Firoozabadi，2001年；Simoni和Secchi，2003年；Monteagudo和Firoozabadi，2004年；Secchi等人，2007年），这可能导致相当大的计算效率降低。

为了克服传统FEM的网格依赖性限制，开发了增强型有限元方法（e-FEMs），它们通过将有限元离散化与不连续性解耦来实现。在这些方法中，扩展有限元方法（XFEM）和通用有限元方法（GFEM）通过在标准有限元函数空间中添加额外的增强函数来建模强不连续性和弱不连续性，即主场变量及其空间梯度的不连续性（Belytschko和Black，1999年；Mo?s等人，1999年；Daux等人，2000年；Fries和Belytschko，2010年）。这些增强函数能够充分解决与不连续性相关的非光滑行为，从而恢复网格与不连续性不对齐时通常会丢失的精度。Khoei（2014年）和Aragón与Duarte（2023年）提供了关于XFEM/GFEM的理论基础、方法论和应用的全面综述。此外，de Borst（2017b年；de Borst，2017a年）还展示了使用XFEM/GFEM和界面元素的裂缝或多孔介质中流体流动的比较研究。

在过去二十年中，XFEM/GFEM在饱和和部分饱和多孔介质中孔隙流体流动的数值分析中得到了广泛应用。关键应用包括具有材料界面的多孔介质的水力-力学模拟（Khoei和Haghighat，2011年；Mohammadnejad和Khoei，2013a年）、裂缝或多孔介质（Réthoré等人，2007年；Réthoré等人，2008年；Callari等人，2010年；Mohammadnejad和Khoei，2013c年；Khoei等人，2023b年；Khoei等人，2023a年）以及水力压裂传播（Gordeliy和Peirce，2013a年；Gordeliy和Peirce，2013b年；Remij等人，2015年；Mohammadnejad和Khoei，2013b年；Mohammadnejad和Andrade，2016年；Khoei等人，2018年；Iranmanesh和Pak，2023年）。其他研究将XFEM/GFEM应用于经历大变形的裂缝多孔介质中的流体流动（Irzal等人，2013年），模拟裂缝和压裂域中的热-水力-力学（THM）耦合（Khoei等人，2012年；Varnosfaderani等人，2017年；Iranmanesh和Pak，2018年；Khoei和Mortazavi，2020年；Khoei等人，2021年；Mortazavi和Khoei，2024年），分析天然裂缝油藏中的热水注入（Mortazavi等人，2022年），以及研究增强型地热系统中的热提取过程（Li等人，2021年；Mortazavi等人，2023年）。

一类独特的e-FEM通过直接在界面处放置增强节点来处理不连续性，这与标准的XFEM/GFEM方法不同，后者是在现有元素节点上应用增强功能。这些所谓的不连续性增强公式可以被视为XFEM/GFEM的衍生物（Aragón和Duarte，2023年），它们将网格与不连续性解耦，同时恢复了XFEM/GFEM中通常会丢失的标准FEM的某些特性。更具体地说，由于增强函数仅在剪切元素中非零，不连续性增强公式具有以下关键优势：
（i）保持标准形状函数的Kronecker-δ属性：由于增强函数在标准元素节点处消失，数值解与节点值完全匹配，确保与标准元素节点相关的自由度（DOFs）保持其物理意义；
（ii）简化基本边界条件的施加：非零Dirichlet边界条件可以像在标准FEM中一样直接应用（Soghrati等人，2012年），包括在浸没的Dirichlet边界上，其中可以弱施加（Cuba-Ramos等人，2015年）或强施加（van den Boom等人，2019年）；
（iii）恢复平滑的反应牵引力：可以从Dirichlet边界准确检索牵引力（van den Boom等人，2019年），这解决了XFEM/GFEM中的一个长期存在的挑战（即使使用稳定化技术也无法实现）；
（iv）混合元素中不会降低精度：由于增强函数在混合元素中（即那些相邻于剪切元素的元素）中确切为零，因此无需在XFEM/GFEM中使用特殊技术来避免某些增强函数选择导致的精度下降（van den Boom等人，2021年；van Bergen等人，2024年）；
（v）内在的数值稳定性：全局刚度矩阵的条件数（影响数值求解器稳定性和效率的关键因素）按Oh?2比例缩放，其中h表示网格元素大小。这种缩放与使用拟合（或几何形状匹配）网格的标准FEM一致，可以通过缩放增强函数使其在标准节点附近的贡献减弱，或通过应用简单的类Jacobi预处理器来实现（Aragón等人，2020年）；
（vi）易于实现：不连续性增强方法计算效率高且易于实现，因为它们允许重用通常在基于位移的FEM代码中使用的标准数据结构。

基于上述优势，在不连续性增强公式类别中，提出了增强界面通用有限元方法（IGFEM）来处理具有弱不连续性的问题（Soghrati等人，2012年；Cuba-Ramos等人，2015年；van den Boom等人，2019年；Aragón等人，2020年）。IGFEM在多项研究中证明了其在无奇异性的平滑问题上的最优收敛性（Soghrati等人，2012年；Cuba-Ramos等人，2015年；van den Boom等人，2019年；De Lazzari等人，2021年），并且在最近的研究中研究了其在涉及具有奇异性的材料连接处处理更复杂情况的能力（Liu等人，2025年）。该方法进一步发展为分层界面增强有限元方法（HIFEM）（Soghrati，2014年；Aragón等人，2020年），用于处理多个不连续性切割单个元素的情况。它还演变为不连续性增强有限元方法（DE-FEM）（Aragón和Simone，2017年；Zhang等人，2019年；Liu等人，2025年；Zhang等人，2025年），通过单三项增强近似同时捕捉弱不连续性和强不连续性。不连续性增强有限元方法已应用于连续介质学的广泛问题，包括多相材料系统的分析（Soghrati等人，2012年；Aragón等人，2013年；Soghrati，2014年；Aragón等人，2020年）、浸没边界（虚拟域）问题（Cuba-Ramos等人，2015年；van den Boom等人，2019年）、断裂力学问题（Aragón和Simone，2017年；van den Boom等人，2019年；Zhang等人，2019年；De Lazzari等人，2021年；Liu等人，2025年；Zhang等人，2025年）、非协调网格耦合和接触问题（Liu等人，2022年）以及拓扑优化问题（van den Boom等人，2021年；van Bergen等人，2024年；van den Boom等人，2023年；Zhang等人，2022年）。然而，这种方法在分析涉及强耦合偏微分方程的复合物理现象方面的适用性仍有待探索。

在本文中，我们介绍了一种新颖的IGFEM框架，用于耦合水力-力学分析分段均匀多孔介质。饱和多孔介质中的材料界面会在固体位移和孔隙压力场的梯度中引起不连续性，这通常被称为弱不连续性。与需要网格与材料界面对齐的标准FEM不同，IGFEM允许界面任意切割元素，同时准确捕捉固体位移和流体压力场中的这些梯度不连续性。这种在剪切元素中的分层增强是该方法的关键创新，使得无需重新网格化即可准确建模弯曲的、几何形状不匹配或其他不规则的材料界面。控制方程（包括线性动量平衡和储存方程）使用IGFEM函数空间进行离散化，时间积分通过Newmark方法进行，最终得到一个使用Newton–Raphson方案整体求解的完全耦合的非线性系统。该方法在四个实际问题上进行了演示：
（i）分层土壤柱的一维固结；
（ii）水平分层土壤结构的二维固结；
（iii）向分层地下储层注水；
（iv）具有正弦材料界面的水平分层土壤的瞬态响应。所有IGFEM结果都与使用拟合网格的标准FEM进行了比较，以突出该方法的精度和鲁棒性，并且第一个示例还与已发表的XFEM/GFEM结果进行了验证。这些演示证实了IGFEM为涉及材料界面间弱不连续性的水力-力学问题提供了一种新颖、高效且灵活的框架。模型构建和丰富的有限元离散化

本研究考虑的饱和多孔介质具有一个固体基质，其中包含相互连接的孔隙，这些孔隙中充满了可迁移的流体，例如气体、水或油。在模型构建中，假设多孔介质处于等温条件下，并且固体基质与流体相之间不发生质量传递或化学反应。流体-力学问题的控制方程和边界条件在第2.1节中介绍，随后在第2.2节中对问题进行了弱形式的表述。第2.3节介绍了使用IGFEM进行空间离散化的方法。第2.4.2节详细介绍了用于时间积分的代数方程耦合系统的整体更新方案。

2.1 控制方程和边界条件

饱和多孔介质的流体-力学行为受到其固体相和流体相之间相互作用的影响。为了理解这些相互作用，将材料点上的总应力σ分解如下是必要的（Zienkiewicz等人，1999年）：
$$
\sigma = \sigma' - \alpha p I,
$$
其中 $ p = p_x,t $ 是流体相中的压力，$ \sigma' = \sigma'x,t $ 是多孔骨架中的有效应力，$ I $ 是二阶单位张量，$ \alpha $ 是Biot系数，$ x $ 和 $ t $ 分别是空间和时间坐标。注意，在方程（1）中 $ p $ 前有一个负号，因为流体压力通常被定义为正值，而根据固体力学惯例，总应力 $ \sigma $（以及有效应力 $ \sigma' $）在拉伸时为正，在压缩时为负。此外，Biot系数由 $ \alpha = \frac{1 - K_t}{K_s} $ 定义，其中 $ K_t $ 是多孔骨架的体积模量，$ K_s $ 是构成多孔骨架的固体相的体积模量。对于大多数土壤来说，单个土壤颗粒的体积模量远大于整个颗粒骨架的体积模量，因此Biot系数接近于1；相反，对于岩石和混凝土，$ \alpha $ 的值通常在1/3到2/3之间（Zienkiewicz等人，1999年）。

多孔介质的流体-力学响应在小变形框架内描述，假设多孔骨架的本构响应是弹性的，即
$$
\sigma' = D \epsilon,
$$
其中 $ D $ 是四阶弹性刚度张量，$ \epsilon = \epsilon_x,t $ 是无穷小应变张量，定义为位移梯度的对称部分，$\epsilon = \frac{1}{2} \nabla u + (\nabla u)^T$，其中 $ u = u_x,t $ 是固体基质的位移，$ \nabla $ 是梯度运算符。

对于饱和多孔介质中的材料点，固体-流体混合物的整体动量平衡由以下方程表示（Zienkiewicz等人，1999年）：
$$
\nabla \cdot \sigma - \rho u'" - \rho f_w' + \nabla w + \rho b = 0.
$$
这里，$ u'" = u'"x,t $ 是固体基质的加速度，$ w_x,t $ 是达西速度，定义为流体相对于固体基质的相对速度，乘以多孔介质的孔隙率。此外，$ w' = w'x,t $ 是流体相的相对加速度，它是相对于附着在固体基质上的参考框架的时间导数 $ w' = \frac{dw}{dt} $。进一步，$ b $ 是单位质量的体力，$ \rho_f $ 和 $ \rho $ 分别是流体相和多孔介质的密度。饱和多孔介质的密度表示为 $ \rho = n\rho_f + (1 - n)\rho_s $，其中 $ \rho_s $ 是固体相的密度，$ n $ 是孔隙率，表示为孔隙体积与多孔介质总体积之比。

另外，相对于附着在固体基质上的参考系统，流体相的动量平衡由以下方程表示（Zienkiewicz等人，1999年）：
$$
-\nabla p - R - \rho_f u'" - \rho_f n w' + \nabla w + \rho_b = 0,
$$
其中 $ R $ 是单位体积的粘性阻力。粘性阻力由达西渗流定律确定，$ R = \frac{w}{k}$，其中 $ (各向同性的) 渗透率 $ k = \kappa/\mu$，$ \kappa $ 是多孔介质的固有渗透率，$ \mu $ 是流体相的动态粘度。

进一步，从固体相和流体相的质量守恒可以推导出饱和多孔介质的所谓储能方程（Zienkiewicz等人，1999年）：
$$
\nabla \cdot w + \alpha \nabla \cdot u' + p' Q = 0,
$$
其中 $ 1/Q $ 通常被称为水力容量，$ K_f $ 是流体相的体积模量。方程（5）考虑了流体流动（通过 $ w $）和不排水压缩（通过 $ p' $）都会导致多孔介质的体积变化。

这三个控制方程（方程（3）、（4）、（5），加上应力分解方程（1）和本构关系方程（2），描述了在一般静态和动态条件下的饱和多孔介质的流体-力学行为。基于方程（3）到（5）的模型称为u-w-p公式，对应于需要求解的未知主变量，即固体相的位移 $ u $、流体相的达西速度 $ w $ 和流体相的压力 $ p $。

方程（3）、（4）中出现的对流流体加速度项 $ w' + \nabla w $ 通常相对较小，因此可以忽略（Zienkiewicz等人，1999年）。相应地，方程（3）简化为
$$
\nabla \cdot \sigma - \rho u'" + \rho b = 0,
$$
而方程（4）的简化形式可以重写为粘性阻力 $ R $ 的表达式，随后调用达西渗流定律并将结果代入方程（5），得到
$$
\nabla \cdot k - \nabla p - \rho_f u'" + \rho_f b + \alpha \nabla \cdot u' + p' Q = 0.
$$
注意，上述两个控制方程的简化集合不依赖于达西速度 $ w $，因此称为u-p公式（Zienkiewicz等人，1999年）。应当注意的是，在涉及高频、短持续时间动态现象的特定问题中，例如由非常快速的加载引起的现象，u-p公式可能会失去准确性。在这种情况下，位移场和孔隙压力场中可能会发展出强烈的流体-力学耦合和振荡响应，当前公式中采用的简化假设可能不再适用。对于这些情况，更通用的u-w-p公式保留了对流流体加速度项，通常能提供更准确的耦合动力学描述（Zienkiewicz等人，1999年）。

为了完成初始-边界值问题的强形式，u-p公式需要补充适当的初始和边界条件。图1a示意性地显示了一个包含饱和多孔介质的域 $\Omega$。该域由外域边界 $\Gamma$ 包围，其法向外单位法向量是 $ n_G$。因此，$\Omega$ 的封闭表示为 $\Omega = \Omega \cup \Gamma$。外域边界由子边界组成，在这些子边界上分别规定了本质（狄利克雷）和自然（诺伊曼）边界条件，分别表示为 $\Gamma_u$ 和 $\Gamma_t$（对于力学问题），以及 $\Gamma_p$ 和 $\Gamma_q$（对于水文问题）。根据 $\Gamma$ 的定义，边界的组成可以表示为 $\Gamma = \Gamma_u \cup \Gamma_t = \Gamma_p \cup \Gamma_q$，条件是子边界是不相交的，即 $\Gamma_u \cap \Gamma_t = \Gamma_p \cap \Gamma_q = \emptyset$。该域由两个材料相 $\Omega_i$（$ i = 1, 2 $）组成，因此 $\Omega = \cup_i \Omega_i$。因此，两个子域之间的材料界面由 $\Gamma_{12} = \Omega_1 \cap \Omega_2$ 定义，其单位法向量由 $ n_{12}$ 给出。此外，时间 $ t = 0 $ 时的初始条件表示为
$$
u_x,0 = u_{0x,0} \in \Omega,
u'x,0 = u'_{0x,0} \in \Omega,
p_x,0 = p_{0x,0} \in \Omega,
$$
其中 $ u_0$ 和 $ u'_{0}$ 分别表示固体相的位移场和速度场的初始值，$ p_0$ 表示流体压力场的初始值。多孔介质的时间依赖行为将分析到预定义的最大时间 $ T $，即 $ t \in [0, T] $。在该时间域内的力学边界条件由以下方程定义
$$
u_x,t = u?_x, t \in \Gamma_u \times [0, T],
t_x,t = n_G x \cdot \sigma_x,t = t?_x, t \in \Gamma_t \times [0, T],
$$
其中 $ t $ 表示牵引力，$ u? $ 和 $ t? $ 分别是在边界处规定的位移和牵引力。此外，材料界面处的力学连续性条件为
$$
???u? = 0 \in \Gamma_{12} \times [0, T],
????t? = ?n_{12} \cdot \sigma? = 0 \in \Gamma_{12} \times [0, T],
$$
其中 ????? 表示在计算材料界面两侧 $ x+ $ 和 $ x- $ 的值后得到的相应场变量的跳变。为了便于解释，在方程（13）和（14）中省略了变量对空间和时间的依赖性，这与边界条件方程（11）和（12）中呈现的依赖性相当。与力学边界条件类似，水文边界条件表示为
$$
p_x,t = p?_x, t \in \Gamma_p \times [0, T],
q_x,t = n_G x \cdot w_x,t = q?_x, t \in \Gamma_q \times [0, T],
$$
其中 $ q $ 是垂直于外边界的流体通量。此外，材料界面处的水文连续性条件为
$$
???p? = 0 \in \Gamma_{12} \times [0, T],
????q? = ?n_{12} \cdot w? = 0 \in \Gamma_{12} \times [0, T],
$$
在方程（16）、（18）中，u-p公式中的达西速度 $ w $ 由以下公式得出
$$
w = - k \nabla p + \rho_f u'" - \rho_f b.
$$

使用上述公式，我们将寻找域 $\Omega$ 内连续的位移场 $ u(x,t)$ 和压力场 $ p(x,t)$ 的解，同时它们的空间梯度在材料界面 $\Gamma_{12}$ 处显示跳跃，即 $ ???n_{12} \cdot \nabla u? \neq 0 $ 和 $ ???n_{12} \cdot \nabla p? \neq 0 $。这些弱不连续性的准确建模确保了使用IGFEM和非共形网格在材料界面上获得的数值结果（见图1b）能够与使用标准FEM和共形网格在材料界面获得的结果进行适当比较。重要的是要提到，方程（13）假设材料界面是连贯的，即在材料界面处两个固体相完美结合且没有损伤。界面处的损伤（裂缝、空洞）或塑性滑移会导致局部位移场的强烈不连续性，这将改变周围多孔介质的流体-力学响应；然而，考虑这一方面超出了本工作的范围。

**下载：** 下载高分辨率图像（257KB）
**下载：** 下载全尺寸图像

图1. (a) 多孔介质的2D示意图，包含域 $\Omega$ 和外边界 $\Gamma$，其法向外单位法向量为 $ n_G$。多孔介质包含一个材料界面 $\Gamma_{12}$，其单位法向量为 $ n_{12}$，它分隔了两个不同的材料相 $\Omega_i$（$ i = 1, 2 $）。该域受到在 $\Gamma_u$ 上规定的位移 $ u?$、在 $\Gamma_t$ 上规定的牵引力 $ t?$、在 $\Gamma_p$ 上规定的流体压力 $ p?$、在 $\Gamma_q$ 上规定的流体通量 $ q?$ 以及单位质量的体力 $ b $。 (b) 材料界面 $\Gamma_{12}$ 处的非共形网格，非共形连续体的标准节点和富集节点分别用黑色和红色实心圆表示。富集节点放置在材料界面与非共形连续体元素的边缘的交点处，从而引入了所谓的切割元素，随后这些元素被进一步细分为四个三角形积分元素。

2.2. 问题的弱形式

通过将两个控制方程（6）、（7）分别乘以允许的测试函数 $ \delta_u $ 和 $ \delta_p $，并在域 $\Omega$ 上对这些表达式进行积分，可以得到这两个方程的弱形式。试验位移场和压力场及其相应的测试函数假定属于适当的函数空间。具体来说，试验函数和测试函数在域 $\Omega$ 上定义，以满足整个域 $\Omega$ 的平方可积性要求，同时在每个材料子域内要求一阶Sobolev规则性，即我们使用“断裂”的Sobolev空间。数学上，我们定义标量值集合 $ V_{\Omega} = \{v \in L^2(\Omega), v_{\Omega_i} \in H^1(\Omega_i) \}$。然后从向量值空间 $ V_u \Omega = \{uu \in V_{\Omega}, u_Gu = 0\}$ 中取出位移场，从 $ V_p \Omega = \{pp \in V_{\Omega}, p_Gp = 0\}$ 中取出压力场。值得注意的是，尽管我们在处理位移场和压力场的均匀狄利克雷边界条件时使用这些空间，但在处理非均匀条件时也可以使用这些空间，只需简单定义所谓的线性多样性（Aragón和Duarte，2023年第2.1.2.2节）。应用散度定理并使用方程（12）、（16）中的定义后，得到
$$
??∫_{\Omega} \sigma : \delta \epsilon d\Omega + \int_{\Omega} \rho u'" \cdot \delta u d\Omega - \int_{\Omega} \rho b \cdot \delta u d\Omega - \int_{\Gamma_t} t? \cdot \delta u d\Gamma_t + \int_{\Gamma_{12}} ?t \cdot \delta u \cdot d\Gamma_{12} = 0 \forall \delta_u \in V_u,
$$
和
$$
??∫_{\Omega} k \nabla p \cdot \nabla \delta p d\Omega + \int_{\Omega} k \rho_f u'" \cdot \nabla \delta p d\Omega + \int_{\Omega} \alpha \nabla \cdot u' \delta p d\Omega + \int_{\Omega} p' Q \delta p d\Omega - \int_{\Omega} k \rho_f b \cdot \nabla \delta p d\Omega + \int_{\Gamma_{12}} ?q \delta p \cdot d\Gamma_{12} = 0 \forall \delta_p \in V_p,
$$
其中 $ \delta \epsilon = \frac{1}{2} \nabla \delta u + (\nabla \delta u)^T $，并且利用了柯西应力张量的对称性 $ \sigma = \sigma^T $。注意，方程（20）、（21）中与材料界面跳跃相关的项会消失，因为方程（13）、（14）、（17）、（18）提供了界面连续性条件。因此，方程（20）、（21）简化为
$$
??∫_{\Omega} \sigma : \delta \epsilon d\Omega + \int_{\Omega} \rho u'" \cdot \delta u d\Omega = \int_{\Omega} \rho b \cdot \delta u d\Omega + \int_{\Gamma_t} t? \cdot \delta u d\Gamma_t \forall \delta_u \in V_u,
$$
和
$$
??∫_{\Omega} k \nabla p \cdot \nabla \delta p d\Omega + \int_{\Omega} k \rho_f u'" \cdot \nabla \delta p d\Omega + \int_{\Omega} \alpha \nabla \cdot u' \delta p d\Omega + \int_{\Omega} p' Q \delta p d\Omega = \int_{\Omega} k \rho_f b \cdot \nabla \delta p d\Omega - \int_{\Gamma_{12}} ?q \delta p \cdot d\Gamma_{12} = 0 \forall \delta_p \in V_p,
$$
在上述表达式中，左边的项依赖于主场变量 $ u $ 和 $ p $ 以及测试函数 $ \delta_u $ 和 $ \delta_p $，而右边的项仅依赖于测试函数。因此，方程（22）、（23）的耦合集合可以正式表示为：
$$
找到 $ u, p \in V_u \times V_p $ 使得
$$
Bu(u, p, \delta_u) = L_u(\delta_u) \forall \delta_u \in V_u
$$
和
$$
Bp(u, p, \delta_p) = L_p(\delta_p) \forall \delta_p \in V_p,
$$
其中双线性形式 $ Bu $ 和 $ Bp $ 以及线性形式 $ Lu $ 和 $ Lp $ 可以分别从方程（22）、（23）的左边和右边推断出来。

2.3.为了求解方程(22)、(23)的有限维对应问题，将域Ω用不遵循材料界面的四边形有限元网格离散化，并采用界面广义有限元方法（IGFEM）来恢复因此丢失的精度。这在图1b中示意性地展示。在材料界面Γ12与非协调连续体元素边缘的交点处引入了富集节点，形成了所谓的切割元素。这些切割元素随后被细分为三角形积分元素，用于对切割元素上的局部数组进行积分。根据Bubnov–Galerkin程序，与固体矩阵中的位移场相关的试函数和检验函数从相同的IGFEM富集空间VuhΩh?VuΩ中选择。离散位移场表示为（25）uhx,t=∑i∈ιhNixu?it︸标准FEM+∑i∈ιwψixαit︸富集，u?i,αi∈R2，检验函数的定义类似。同样，流体压力场的试函数和检验函数也从有限维空间VphΩh?VpΩ中选取，离散压力场表示为（26）ph(x,t)=∑i∈ιhNixp?it︸标准FEM+∑i∈ιwψixγit︸富集，p?i,γi∈R，检验函数的定义方式也相同。在方程（25）、（26）中，第一项对应于标准FEM近似，其中ιh表示原始网格的标准节点索引集（图1b中的黑色实心圆表示），第二项对应于富集，ιw表示切割元素中的富集节点索引集（图1b中的红色实心圆表示）。在第一项中，Ni是与第i个网格节点相关的标准拉格朗日形状函数，u?i和p?i分别代表固体矩阵位移和流体压力的标准自由度。在第二项中，ψi是富集函数，αi和γi分别代表与第i个富集节点相关的固体位移和流体压力的富集自由度。通过在相邻切割元素之间的交界处共享富集节点来确保C0-连续性，因为这些节点通常位于相邻元素的公共边缘上。富集项通过弱不连续的富集函数ψi增强了标准有限元近似，这些函数包含了沿材料界面Γ12的固体矩阵位移梯度和流体压力梯度的不连续性。考虑材料界面将一个四边形切割元素分成两个四边形子域的情况（见图2），我们称这些子域为子元素。在这种情况下，富集函数ψi是通过这两个四边形子域中的拉格朗日形状函数构建的（Soghrati等人，2012年）。富集函数在相应的富集节点处取最大值1，在切割元素的标准节点处减小到0。因此，在混合元素中（即与切割元素相邻的元素）富集函数完全为零；因此，由于混合元素中的寄生项，不会出现精度降低或收敛速度下降的情况。子元素进一步被划分为积分元素，这些元素专用于对局部元素数组进行高斯求积，如下文更详细解释的那样。值得注意的是，材料界面也可能将一个四边形父元素分成两个三角形子域（界面穿过两个相对的节点），或者创建一个三角形和一个五边形子域（界面切割四边形的两个相邻边缘）。IGFEM也能够处理这种情况，如Aragón等人（2020年）所解释的；在第3.4节介绍的流体力学问题背景下提供了关于四边形父元素分成三角形和五边形子域的详细处理。此外，IGFEM已通过分层界面富集有限元方法（HIFEM）（Soghrati，2014年；Aragón等人，2020年）扩展到处理多个界面与单个元素相交的情况。在HIFEM中，富集函数是分层构建的，由初始界面生成的子元素作为后续由其他界面形成的子元素的父元素，从而在有序的树结构中组织层次结构。

为了求解方程（22）、（23）的有限维对应问题，域Ω用不遵循材料界面的四边形有限元网格离散化，并采用界面广义有限元方法（IGFEM）来恢复因此丢失的精度。这在图1b中示意性地展示。在材料界面Γ12与非协调连续体元素边缘的交点处引入了富集节点，形成了所谓的切割元素。这些切割元素随后被细分为三角形积分元素，用于对切割元素上的局部数组进行积分。根据Bubnov–Galerkin程序，与固体矩阵中的位移场相关的试函数和检验函数从相同的IGFEM富集空间VuhΩh?VuΩ中选择。离散位移场表示为（25）uhx,t=∑i∈ιhNixu?it︸标准FEM+∑i∈ιwψixαit︸富集，u?i,αi∈R2，检验函数的定义类似。同样，流体压力场的试函数和检验函数也从有限维空间VphΩh?VpΩ中选取，离散压力场表示为（26）ph(x,t)=∑i∈ιhNixp?it︸标准FEM+∑i∈ιwψixγit︸富集，p?i,γi∈R，检验函数的定义方式也相同。在方程（25）、（26）中，第一项对应于标准FEM近似，其中ιh表示原始网格的标准节点索引集（图1b中的黑色实心圆表示），第二项对应于富集，其中ιw表示切割元素中的富集节点索引集（图1b中的红色实心圆表示）。在第一项中，Ni是与第i个网格节点相关的标准拉格朗日形状函数，u?i和p?i分别代表固体矩阵位移和流体压力的标准自由度。在第二项中，ψi是富集函数，αi和γi分别代表与第i个富集节点相关的固体位移和流体压力的富集自由度。通过在不同切割元素之间共享富集节点来确保C0-连续性，因为这些节点位于通常为相邻元素共有的边缘上。富集项通过弱不连续的富集函数ψi增强了标准有限元近似，这些函数包含了沿材料界面Γ12的固体矩阵位移梯度和流体压力梯度的不连续性。考虑材料界面将一个四边形切割元素分成两个四边形子域的情况（见图2），我们称这些子域为子元素。在这种情况下，富集函数ψi是通过这两个四边形子域中的拉格朗日形状函数构建的（Soghrati等人，2012年）。富集函数在相应的富集节点处取得最大值1，在切割元素的标准节点处减小到0。因此，在混合元素中（即与切割元素相邻的元素）富集函数完全为零；结果，由于混合元素中的寄生项，不会出现精度降低或收敛速度下降的情况。子元素进一步被划分为积分元素，这些元素专门用于对局部元素数组进行高斯求积，如下文更详细解释的那样。值得注意的是，材料界面也可能将一个四边形父元素分成两个三角形子域（界面穿过两个相对的节点），或者创建一个三角形和一个五边形子域（界面切割四边形的两个相邻边缘）。IGFEM也能够处理这种情况，如Aragón等人（2020年）所解释的；在第3.4节介绍的流体力学问题背景下提供了关于四边形父元素分成三角形和五边形子域的详细处理。此外，IGFEM已通过分层界面富集有限元方法（HIFEM）（Soghrati，2014年；Aragón等人，2020年）扩展到处理多个界面与单个元素相交的情况。在HIFEM中，富集函数是分层构建的，由初始界面生成的子元素作为后续由其他界面形成的子元素的父元素，从而在有序的树结构中组织层次结构。

为了求解方程（22）、（23）的有限维对应问题，将域Ω用不遵循材料界面的四边形有限元网格离散化，并采用界面广义有限元方法（IGFEM）来恢复因此丢失的精度。这在图1b中示意性地展示。在材料界面Γ12与非协调连续体元素边缘的交点处引入了富集节点，形成了所谓的切割元素。这些切割元素随后被细分为三角形积分元素，用于对切割元素上的局部数组进行积分。根据Bubnov–Galerkin程序，与固体矩阵中的位移场相关的试函数和检验函数从相同的IGFEM富集空间VuhΩh?VuΩ中选择。离散位移场表示为（25）uhx,t=∑i∈ιhNixu?it︸标准FEM+∑i∈ιwψixαit︸富集，u?i,αi∈R2，检验函数的定义类似。同样，流体压力场的试函数和检验函数也从有限维空间VphΩh?VpΩ中选取，离散压力场表示为（26）ph(x,t)=∑i∈ιhNixp?it︸标准FEM+∑i∈ιwψixγit︸富集，p?i,γi∈R，检验函数的定义方式也相同。在方程（25）、（26）中，第一项对应于标准FEM近似，其中ιh表示原始网格的标准节点索引集（图1b中的黑色实心圆表示），第二项对应于富集，其中ιw表示切割元素中的富集节点索引集（图1b中的红色实心圆表示）。在第一项中，Ni是与第i个网格节点相关的标准拉格朗日形状函数，u?i和p?i分别代表固体矩阵位移和流体压力的标准自由度。在第二项中，ψi是富集函数，αi和γi分别代表与第i个富集节点相关的固体位移和流体压力的富集自由度。通过在不同切割元素之间共享富集节点来确保C0-连续性，因为这些节点位于通常为相邻元素共有的边缘上。富集项通过弱不连续的富集函数ψi增强了标准有限元近似，这些函数包含了沿材料界面Γ12的固体矩阵位移梯度和流体压力梯度的不连续性。考虑材料界面将一个四边形切割元素分成两个四边形子域的情况（见图2），我们称这些子域为子元素。在这种情况下，富集函数ψi是通过这两个四边形子域中的拉格朗日形状函数构建的（Soghrati等人，2012年）。富集函数在相应的富集节点处取得最大值1，在切割元素的标准节点处减小到0。因此，在混合元素（即与切割元素相邻的元素）中，富集函数完全为零；结果，由于混合元素中的寄生项，不会出现精度降低或收敛速度下降的情况。子元素进一步被划分为积分元素，这些元素专门用于对局部元素数组进行高斯求积，如下文更详细解释的那样。值得注意的是，材料界面也可能将一个四边形父元素分成两个三角形子域（界面穿过两个相对的节点），或者创建一个三角形和一个五边形子域（界面切割四边形的两个相邻边缘）。IGFEM也能够处理这种情况，如Aragón等人（2020年）所解释的；在第3.4节介绍的流体力学问题背景下提供了关于四边形父元素分成三角形和五边形子域的详细处理。此外，IGFEM已通过分层界面富集有限元方法（HIFEM）（Soghrati，2014年；Aragón等人，2020年）扩展到处理多个界面与单个元素相交的情况。在HIFEM中，富集函数是分层构建的，由初始界面生成的子元素作为后续由其他界面形成的子元素的父元素，从而在有序的树结构中组织层次结构。

为了求解方程（22）、（23）的有限维对应问题，域Ω用不遵循材料界面的四边形有限元网格离散化，并采用界面广义有限元方法（IGFEM）来恢复因此丢失的精度。这在图1b中示意性地展示。在材料界面Γ12与非协调连续体元素边缘的交点处引入了富集节点，形成了所谓的切割元素。这些切割元素随后被细分为三角形积分元素，用于对切割元素上的局部数组进行积分。根据Bubnov–Galerkin程序，与固体矩阵中的位移场相关的试函数和检验函数从相同的IGFEM富集空间VuhΩh?VuΩ中选择。离散位移场表示为（25）uhx,t=∑i∈ιhNixu?it︸标准FEM+∑i∈ιwψixαit︸富集，u?i,αi∈R2，检验函数的定义类似。同样，流体压力场的试函数和检验函数也从有限维空间VphΩh?VpΩ中选取，离散压力场表示为（26）ph(x,t)=∑i∈ιhNixp?it︸标准FEM+∑i∈ιwψixγit︸富集，p?i,γi∈R，检验函数的定义方式也相同。在方程（25）、（26）中，第一项对应于标准FEM近似，其中ιh表示原始网格的标准节点索引集（图1b中的黑色实心圆表示），第二项对应于富集，其中ιw表示切割元素中的富集节点索引集（图1b中的红色实心圆表示）。在第一项中，Ni是与第i个网格节点相关的标准拉格朗日形状函数，u?i和p?i分别代表固体矩阵位移和流体压力的标准自由度。在第二项中，ψi是富集函数，αi和γi分别代表与第i个富集节点相关的固体位移和流体压力的富集自由度。通过在不同切割元素之间共享富集节点来确保C0-连续性，因为这些节点位于通常为相邻元素共有的边缘上。富集项通过弱不连续的富集函数ψi增强了标准有限元近似，这些函数包含了沿材料界面Γ12的固体矩阵位移梯度和流体压力梯度的不连续性。考虑材料界面将一个四边形切割元素分成两个四边形子域的情况（见图2），我们称这些子域为子元素。在这种情况下，富集函数ψi是通过这两个四边形子域中的拉格朗日形状函数构建的（Soghrati等人，2012年）。富集函数在相应的富集节点处取得最大值1，在切割元素的标准节点处减小到0。因此，在混合元素中（即与切割元素相邻的元素）富集函数完全为零；结果，由于混合元素中的寄生项，不会出现精度降低或收敛速度下降的情况。子元素进一步被划分为积分元素，这些元素专门用于对局部元素数组进行高斯求积，如下文更详细解释的那样。值得注意的是，材料界面也可能将一个四边形父元素分成两个三角形子域（界面穿过两个相对的节点），或者创建一个三角形和一个五边形子域（界面切割四边形的两个相邻边缘）。IGFEM也能够处理这种情况，如Aragón等人（2020年）所解释的；在第3.4节介绍的流体力学问题背景下提供了关于四边形父元素分成三角形和五边形子域的详细处理。此外，IGFEM已通过分层界面富集有限元方法（HIFEM）（Soghrati，2014年；Aragón等人，2020年）扩展到处理多个界面与单个元素相交的情况。在HIFEM中，富集函数是分层构建的，由初始界面生成的子元素作为后续由其他界面形成的子元素的父元素，从而在有序的树结构中组织层次结构。

界面增强型广义有限元离散化

为了求解方程（22）、（23）的有限维对应问题，域Ω使用不遵循材料界面的四边形有限元网格进行离散化，并采用界面增强型有限元方法（IGFEM）来恢复因此丢失的精度。这如图1b所示。在材料界面Γ12与非协调连续体元素的边缘交点处引入了增强节点，形成了所谓的切割元素。这些切割元素随后被细分为三角形积分元素，用于对切割元素上的局部数组进行积分。遵循Bubnov–Galerkin程序，与固体矩阵中的位移场相关的试函数和检验函数从相同的IGFEM增强空间VuhΩh?VuΩ中选取。离散位移场表示为（25）uhx,t=∑i∈ιhNixu?it︸标准FEM+∑i∈ιwψixαit︸增强，u?i,αi∈R2，检验函数的定义类似。同样，流体压力场的试函数和检验函数也从有限维空间VphΩh?VpΩ中选取，离散压力场表示为（26）ph(x,t)=∑i∈ιhNixp?it︸标准FEM+∑i∈ιwψixγit︸增强，p?i,γi∈R。在方程（25）、（26）中，第一项对应于标准FEM近似，其中ιh表示原始网格的标准节点索引集（图1b中的黑色实心圆），第二项对应于增强，其中ιw表示切割元素的增强节点索引集（图1b中的红色实心圆）。在第一项中，Ni是与第i个网格节点相关的标准拉格朗日形状函数，u?i和p?i分别代表固体矩阵位移和流体压力的标准自由度。在第二项中，ψi是增强函数，αi和γi分别代表与第i个增强节点相关的固体位移和流体压力的增强自由度。通过在不同切割元素之间共享增强节点来确保C0-连续性，因为这些节点位于通常为相邻元素共有的边缘上。增强项通过弱不连续的增强函数ψi增强了标准有限元近似，这些函数包含了沿材料界面Γ12的固体矩阵位移梯度和流体压力梯度的不连续性。考虑材料界面将一个四边形切割元素分成两个四边形子域的情况（见图2），我们将其称为子这可以通过以下两种方法实现：（i）从雅可比矩阵中排除惯性矩阵C（在计算J时可以忽略它，因为其贡献相对较小（Mohammadnejad和Khoei，2013c；Khoei和Mortazavi，2020）；（ii）将方程（36）乘以?β2/(β1Δt)的值。因此，方程（38）变为（40）JsymmδUnδPni=?β2β1ΔtΨn+1UΨn+1Pi，对称化的雅可比矩阵表示为（41）Jsymm=?β2β1Δt1β1Δt2M+Kβ2β1ΔtGsymm.H+1β?1ΔtS。当残差Ψn+1U和Ψn+1P的组合L2范数降到接近零的预定容忍度以下时，迭代更新过程（方程（38）被认为是收敛的。

3. 数值结果
在本节中，我们展示了IGFEM在处理具有弱不连续性的饱和多孔介质的力学分析中的准确性和鲁棒性。第3.1节中的第一个例子涉及一个由两种不同沙材料组成的饱和土柱，在垂直荷载作用下进行一维固结。第3.2节中的第二个例子探讨了由两层不同沙层组成的土壤结构在垂直表面荷载作用下的二维固结行为。第3.3节模拟了水注入由三种不同岩材料组成的地下储层的过程。最后，第3.4节考虑了一个具有正弦形材料界面的水平分层土壤的瞬态响应，以证明IGFEM能够在不需要网格对齐的情况下准确捕捉复杂的非平面界面。所有示例都使用了标准FEM（在材料界面处使用共形有限元素）和IGFEM（在材料界面处使用非共形有限元素）进行分析。FEM结果作为基准解，用于验证使用IGFEM方法获得的结果。数值结果是通过在MATLAB数值计算平台上实现第2节中的模型公式计算得出的。有限元网格使用四节点四边形（Q4）元素构建，对位移和压力场都采用双线性插值，并且每个领域都使用2 × 2高斯积分。在IGFEM分析中，对于富集元素的面积，采用三角形积分元素进行数值积分。
众所周知，标准的u–p公式在违反Ladyzhenskaya–Babu?ka–Brezzi（LBB）条件时会出现数值稳定性问题，特别是在两个极限情况下：（i）不排水极限（w→0）和（ii）刚性骨架极限（u→0），参见Lotfian和Sivaselvan（2018）。在这些情况下，系统行为趋近于不可压缩状态，因此选择用于位移和压力的有限元素空间变得至关重要。LBB条件提供了一个稳定性标准，确保这些场之间的正确耦合；如果不满足该条件，可能会导致checkerboard压力模式、锁定以及收敛不良——这些现象在标准的等阶位移和压力插值中很常见。在本研究中，材料参数和边界条件自然远离这些极限情况，从而放宽了不可压缩性约束，并减轻了LBB条件的严格性。在这种条件下，使用双线性插值对于位移和压力场都能获得稳定且准确的解，没有出现虚假的压力振荡或锁定的现象。

下载：下载高分辨率图像（222KB）
下载：下载完整尺寸图像
图4. 一个由两种不同沙材料“1”和“2”组成的饱和土柱，在中部的材料界面处分开，在其顶面施加的时间依赖性垂直荷载t?y的作用下进行固结。（a）荷载剖面、几何形状和边界条件。（b）在材料界面使用共形有限元素的FEM网格（用红色虚线表示）。（c）在材料界面使用非共形有限元素的IGFEM网格。

3.1. 分层土柱的一维固结
考虑Khoei和Haghighat（2011）提出的一维基准问题，该问题涉及一个由两种不同沙材料“1”和“2”组成的饱和土柱，在中部的材料界面处分开，见图4。这个柱子在顶面施加的时间依赖性垂直荷载作用下进行固结。这个问题使得能够将IGFEM的结果与FEM的结果以及Khoei和Haghighat（2011）得到的XFEM/GFEM结果进行准确比较。问题的几何形状、加载和边界条件显示在图4a中。两种沙材料的杨氏模量分别为E1=60×10^6 N/m^2和E2=30×10^6 N/m^2，固有渗透率分别为κ1=1.02×10^-12 m^2和κ2=1.02×10^-11 m^2。两种土壤的其余属性相同，泊松比为ν=0.2，固体和流体的密度分别为ρs=2000 kg/m^3和ρf=1000 kg/m^3，固体和流体相的体积模量分别为Ks=1.0×10^20 N/m^2和Kf=2.1×10^9 N/m^2，Biot系数为α=1，孔隙率为n=0.3，流体相（即地下水）的动粘度为μ=10^-3 N·s/m^2。为了简化，忽略了重力加载效应。所有参数值与Khoei和Haghighat（2011）报告中相同问题的参数值一致，需要注意的是，上述体积模量Ks的值人为地较高，即假设沙颗粒几乎是“理想刚性”的。尽管如此，这里未展示的额外模拟表明，这一假设并不显著影响数值结果；砂岩的颗粒体积模量通常相对较高，大约在10到60 GPa之间（Qin等人，2022），这验证了分析中采用的“理想刚性”沙颗粒的假设。多孔介质的初始位移、速度和超孔隙压力被规定为零。在土柱顶面施加的垂直牵引力从t?y=0到1 kN/m^2线性增加，持续时间为0.1秒，然后保持不变。底面受到刚性支撑，u?x=u?y=0，两个侧面的法向水平位移也为零，u?x=0。在加载过程中，顶面的孔隙压力被规定为零，p?=0，而在其他三个（不透水的）外部边界，流体通量为零，q?=0。此外，假设平面应变条件。计算域分别使用10个Q4元素进行标准FEM分析，以及9个Q4元素进行IGFEM分析，如图4b和4c所示。时间积分方案中使用的Newmark参数为β1=β?1=0.6和β2=0.605，时间增量从t=0对数增加到t=200秒，共1000步。空间离散化和计算设置与Khoei和Haghighat（2011）中应用的设置相似。

图5显示了在给定荷载作用下土柱的时间依赖性响应，其中IGFEM方法（红色虚线）、标准FEM方法（蓝色实线）以及Khoei和Haghighat（2011）报告的XFEM/GFEM方法（黑色虚线点线）。图5a和5b分别显示了在y= 9 m和27 m处评估的超孔隙压力p，图5c和5d分别显示了在y= 24 m和30 m处的垂直位移uy。图5中的插图阐明了在t= 0.1 s到t= 1 s之间评估的IGFEM和标准FEM分析的初始瞬态响应。由于Khoei和Haghighat（2011）没有详细报告XFEM/GFEM模拟的初始瞬态响应，因此该结果没有出现在插图中。可以看出，IGFEM计算出的结果与标准FEM计算出的结果在初始瞬态阶段和整个时间域内都非常吻合。相比之下，Khoei和Haghighat（2011）报告的XFEM/GFEM结果显示出压力和位移响应存在一些不准确性。这种不准确性的一个潜在原因可能是在线性混合元素中对主要变量的空间近似中出现了寄生非线性项，如第2.3节所讨论的。值得注意的是，Khoei和Haghighat（2011）中展示的XFEM/GFEM结果仅作为参考。由于他们的公式没有在我们的计算框架中实现，因此无法评估或直接与IGFEM结果比较诸如收敛率、矩阵条件性和计算成本等指标。

图5a显示，在大部分考虑的时间内，结构下部的超孔隙压力p（y= 9 m）几乎保持不变，而结构上部的孔隙压力p（y= 27 m）逐渐减小到零。结构上部孔隙水更快地消散显然是由于孔隙水从顶面离开结构，因为在顶面规定的外部压力为零。图5c显示，在y= 24 m处的垂直位移uy最初由于孔隙压力的瞬态效应略有振荡，大约10秒后以几乎线性的方式单调增加。最终在200秒后获得的位移值uy=?0.14mm是在上表面y= 30 m观察到的最终位移值uy=?0.28mm的一半，图5d所示，这说明了位移随深度的快速减小。

下载：下载高分辨率图像（641KB）
下载：下载完整尺寸图像
图5. 使用标准FEM模型（蓝色实线）、IGFEM模型（红色虚线）以及Khoei和Haghighat（2011）提出的XFEM/GFEM模型（黑色虚线线）计算的图4中所描述固结问题的时间依赖性响应。（a）y= 9 m处的超孔隙压力p。（b）y= 27 m处的超孔隙压力p。（c）y= 24 m处的垂直位移uy。（d）y= 30 m处的垂直位移uy。插图详细显示了使用标准FEM模型（蓝色实线）和IGFEM模型（红色虚线）在t= 0.1 s到t= 1 s之间计算的初始瞬态响应。
图6a和6b分别显示了在t=200s时沿柱高度的垂直位移和孔隙压力分布。使用IGFEM获得的结果与使用标准FEM计算的结果非常吻合。在位于土柱中间高度y=15m的材料界面处，观察到两个剖面的斜率有明显变化。虽然垂直位移在界面处保持连续，但其梯度突然改变（见图6a），表明上层土壤的刚度较低，因此沉降较大。同样，孔隙压力剖面在界面处连续，但在界面处有一个明显的分界点（见图6b），这是由于两种土壤层之间固有渗透率的差异造成的。

下载：下载高分辨率图像（280KB）
下载：下载完整尺寸图像
图6. 在t=200s时，使用标准FEM模型（蓝色实线）和IGFEM模型（红色虚线）计算的沿柱高度的垂直位移和孔隙压力分布。（a）垂直位移uy。（b）超孔隙压力p。

3.2. 水平分层土壤结构的二维固结
第二个问题涉及在表面荷载和平面应变条件下，一个饱和的、水平分层的土壤结构的二维固结行为，见图7a。这种分析的实际好处在于表明施加表面荷载对于管理固结至关重要；它有助于以安全且成本效益高的方式减少未来的沉降并控制施工时间表。作为岩土设计中的主动措施，它直接影响建在软土或可压缩土壤上的基础设施的性能和长期耐久性。在数值问题中，土壤结构的垂直截面积为30×30m^2，由两个水平沙层组成，这两个沙层在距结构底部y= 24 m处的材料界面处分开。固结过程在上方表面施加的时间依赖性垂直牵引力t?y的作用下发生，该牵引力在距右侧边界6 m的水平距离上施加。上方结构的其余部分被建模为无牵引力，且上表面的压力等于零（即大气压）。与前一个问题类似，底部表面受到刚性支撑，u?x=u?y=0，两个侧面的法向位移也设置为零，u?x=0。此外，两个侧面和底部表面的通量也设置为零，q?=0。在时间积分方案中使用的Newmark参数以及材料参数与第3.1节中介绍的第一个问题相同，除了两种土壤材料的内禀渗透率；这里的渗透率值稍大一些，分别为κ1=5×10?11 m2和κ2=10?1? m2，以便更清楚地突出固结过程中时间依赖响应的特定特征。多孔介质内部的初始位移、速度和超孔隙压力被设定为零。顶部表面的牵引力从t?y=0线性增加到t?y=350 kN/m2，之后保持恒定。计算域使用结构化的网格进行离散化，标准FEM分析使用729个Q4单元，IGFEM分析使用702个Q4单元，如图7b和7c所示。空间离散化的选择是基于对三种不同分辨率的结构化FEM网格的网格收敛性研究：一个粗网格包含225个Q4单元，一个中间网格包含729个Q4单元，一个细网格包含2916个Q4单元，见图8。时间依赖的流体-力学响应在结构域内的点A处进行评估，该点位于x=27m,y=21m（参见图7a），并比较了三种网格配置下的响应。图9展示了点A处超孔隙压力p和垂直位移uy随时间的演变。图9中的插图阐明了t=1 s到t=5 s之间的响应（以对数时间尺度显示）。如图所示，与中间和细网格相比，粗网格显示出小但明显的偏差，特别是在位移响应方面。相比之下，中间和细网格的结果非常一致，表明网格收敛。基于这一观察，最终分析选用了中间网格，因为它在准确性和计算效率之间提供了良好的平衡。除了空间离散化外，还使用中间网格研究了时间分辨率的影响。在模拟的时间区间t∈[0,200s]内考虑了三种对数增加的时间步进方案：第一种方案有500个时间步长，第二种方案有1000个时间步长，第三种方案有1500个时间步长。对于所有情况，Newmark积分参数都设置为β1=β?1=0.6和β2=0.605，以确保模拟过程中的数值阻尼和稳定性一致。点A处超孔隙压力p和垂直位移uy的相应时间历史在图10中显示。图9的插图更详细地展示了t=1 s到t=5 s之间的响应（使用对数时间尺度）。使用1000个和1500个时间步长得到的结果非常吻合，而使用500个时间步长的结果略有偏差，表明在1000个时间步长时达到了时间收敛。因此，在最终分析中采用了更高效的1000个时间步长的方案。图11a和11b分别展示了在x=27m,y=21m的点A和x=27m,y=27m的点B处，IGFEM和FEM方法计算的超孔隙压力p和垂直位移uy随时间的变化。显然，这两种方法计算的压力和位移响应是相同的。随着表面牵引力增加到最大值，点A和B处的孔隙压力从零上升到峰值。一旦表面牵引力达到并保持在这个最大值，由于土壤结构的固结，两点的孔隙压力逐渐减小到零。点B处的孔隙压力变化比点A处更大，这归因于第二层土壤的渗透率高于第一层土壤。大约在t=0.1s时，点A和B处的垂直位移开始从零增加，在t=1s左右出现轻微的振荡现象，这是由于质量惯性效应造成的。这些效应是由施加的表面牵引力达到最大值时孔隙压力的突然峰值引发的。随着时间的推移，由于固结过程的特点，点A和B处的垂直位移进一步增加。

下载：下载高分辨率图像（469KB）
下载：下载全尺寸图像
图7. 一个饱和的分层土壤结构在其顶部表面受到局部、时间依赖的垂直载荷t?y的作用而发生固结。(a) 载载剖面、几何形状和边界条件。(b) 材料界面处使用共形有限元的结构化FEM网格（由红色虚线表示)。(c) 材料界面处使用非共形有限元的结构化IGFEM网格。
下载：下载高分辨率图像（461KB）
下载：下载全尺寸图像
图8. 用于固结问题网格收敛性研究的有限元网格。(a) 含有225个Q4单元的粗网格。(b) 含有729个Q4单元的中间网格。(c) 含有2916个Q4单元的细网格。
下载：下载高分辨率图像（425KB）
下载：下载全尺寸图像
图9. 图8中显示的三种网格分辨率对应的固结问题的时间依赖FEM响应，即粗网格（黑色实线）、中间网格（蓝色实线）和细网格（粉色实线）。(a) 图7a中指示的点A处的孔隙压力p。(b) 图7a中指示的点A处的垂直位移uy。插图详细展示了t=1 s到t=5 s之间的响应（以对数时间尺度显示）。
下载：下载高分辨率图像（416KB）
下载：下载全尺寸图像
图10. 对于500时间步长（黑色实线）、1000时间步长（蓝色实线）和1500时间步长（粉色实长）的时间依赖FEM响应。(a) 图7a中指示的点A处的孔隙压力p。(b) 图7a中指示的点A处的垂直位移uy。插图详细展示了t=1 s到t=5 s之间的响应（以对数时间尺度显示）。
图12和图13分别展示了t=20 s和t=200 s时超孔隙水压力p和垂直位移uy的总体分布。在两个评估时间点，IGFEM和FEM的模拟结果在整个结构中完全一致。固结过程中的孔隙水压力从非均匀的高压力剖面（t=20 s）发展成整个结构宽度上均匀的低压力剖面（t=200 s）。如图12a和13a所示，随着时间的进一步推移，整个结构的超孔隙压力趋于逐渐接近零（导致固结完成），此时位移场稳定并变得时间不变，其最大值位于施加的表面载荷下方。
下载：下载高分辨率图像（359KB）
下载：下载全尺寸图像
图11. 图7中所示的固结问题的时间依赖响应，使用标准FEM模型（蓝色实线）和IGFEM模型（红色虚线）计算。(a) 图7a中指示的点A和B处的超孔隙压力p。(b) 图7a中指示的点A和B处的垂直位移uy。
下载：下载高分辨率图像（449KB）
下载：下载全尺寸图像
图12. 使用标准FEM模型（左）和IGFEM模型（右）计算的图7中所示的固结问题的场变量等高线图。(a) t=20 s时整个结构上的超孔隙压力p。(b) t=20 s时整个结构上的垂直位移uy。
下载：下载高分辨率图像（391KB）
下载：下载全尺寸图像
图13. 使用标准FEM模型（左）和IGFEM模型（右）计算的图7中所示的固结问题的场变量等高线图。(a) t=200 s时整个结构上的超孔隙压力p。(b) t=200 s时整个结构上的垂直位移uy。
3.3. 向分层地下储层注水
第三个问题考虑向地下储层注水，见图14。这是多孔介质研究中的关键领域（Monteagudo和Firoozabadi，2004；Hoteit和Firoozabadi，2008；Hajiabadi和Khoei，2019；Mortazavi等人，2022；Khoei等人，2023a），此类分析有助于表征储层属性，如渗透率、孔隙度和异质性，从而支持关于井位和生产策略的更明智的决策。此外，检查注水导致的孔隙压力变化对于评估与诱发地震、井筒不稳定性和潜在岩石破坏相关的风险至关重要。在数值示例中，储层的水平截面积为600×600m2，由三个不同的垂直岩层组成，每个厚度为200米，由两个材料界面分隔，如图14a和14b所示。岩石材料的杨氏模量分别为E1=4×101?N/m2、E2=8×101?N/m2和E3=12×101?N/m2，内禀渗透率分别为κ1=1.97×10?1?m2（20mD）、κ2=9.87×10?1?m2（10mD）和κ3=4.93×10?1?m2（5mD）。其余材料属性与第3.1节中处理的第一个问题相同。该问题使用二维模型进行模拟，该模型在深度方向上受到平面应变条件的约束。域使用结构化网格进行离散化，标准FEM分析使用900个Q4单元，IGFEM分析使用840个Q4单元，如图14c和14d所示。进一步细化网格对数值结果的准确性没有显著影响。
在外部边界，储层被刚性支撑，u?x=u?y=0，且不透水，q?=0，见图14b。水以恒定速率1.25 × 10?? m3/s（相当于0.0001 PV/d）在位于二维几何形状左下角的注入点注入。在位于二维几何形状右上角的生产点，水压被规定为零，p?=0。储层内部的初始位移、速度和超孔隙压力被设定为零。时间步进方案采用固定的时间增量Δt=3h，Newmark参数设置为β1=β?1=0.6和β2=0.605，这些参数已被证实可以产生收敛的结果。图15a展示了注入点处流体压力p的时间演变（即注入压力），而图15b展示了注入开始后1500天内储层中心的总水平位移u=ux2+uy2。此外，图15c显示了生产流量与注入流量的比值|Qp/Qi|，图15d展示了同一时期累积产水量与累积注入体积的比值|Vp/Vi|。显然，使用IGFEM获得的结果与标准FEM方法的结果完全一致。在注入点瞬时流动开始导致孔隙压力初始瞬态跳跃后，注入点的压力稳步增加，因为孔隙水缓慢地通过分层储层向生产点迁移。这导致储层中心的水平位移增加，同时生产流量与注入流量的比值和累积生产量与注入体积的比值也增加。随着时间的推移，多孔介质内的流场逐渐趋于稳态。在这个阶段，注入点的压力在100 MPa以上逐渐稳定，见图15a；中心点的位移达到最终值约为13 mm，见图15b；生产流量与注入流量的比值趋近于1，见图15c。累积产水量与注入体积的比值在0.8以上趋于稳定，见图15d，其中与1的偏差归因于多孔介质中保留的水体积。
下载：下载高分辨率图像（863KB）
下载：下载全尺寸图像
图14. 向由三个不同垂直岩层组成的地下储层注水。(a) 储层的几何形状以及注入井和生产井的位置。(b) 二维模型几何形状，包括注水点和生产点以及边界条件。(c) 两个材料界面处使用共形有限元的结构化FEM网格（由红色虚线表示)。(d) 两个材料界面处使用非共形有限元的结构化IGFEM网格。图16a和16b分别展示了t=1500天时储层上的超孔隙压力和位移剖面，这大约对应于最终的稳态流动条件。这些图形清楚地表明IGFEM和标准FEM模型产生了相同的结果。请注意，从注入点到生产点，孔隙压力逐渐减小，最大水平位移大致发生在储层的中心。尽管生产点的流体压力被规定为零，但在等高线图中显示为50MPa，这对应于选定的压力范围的下限，选择该范围是为了优化孔隙压力场的视觉表示。

图15. 图14所示的向地下储层注水问题的时变响应，分别用标准FEM模型（蓝色实线）和IGFEM模型（红色虚线）计算。(a) 注入点的流体压力p。(b) 储层中心的总水平位移u=ux2+uy2。(c) 生产流量与注入流量的比值|Qp/Qi|。(d) 累计产水体积与累计注水体积的比值|Vp/Vi|。

图16. 图14所示的向地下储层注水问题的场变量等高线图，分别用标准FEM模型（左）和IGFEM模型（右）计算。(a) t=1500时刻结构上的孔隙压强p。(b) t=1500时刻结构上的总水平位移u=ux2+uy2。

3.4. 具有正弦形材料界面的水平分层土壤的瞬态响应

第四个问题旨在展示所提出的IGFEM公式在具有几何复杂材料界面的多孔介质水力-力学分析中的准确性和鲁棒性。具体来说，这个例子突显了IGFEM在以下情况下的有效性：标准有限元方法需要复杂的网格对齐来准确表示界面几何形状，而IGFEM使用不遵循界面的结构化网格，能够准确且高效地解决水力-力学问题。

图17展示了一个垂直截面为25×25m2的土壤结构，该结构由两个不同的沙层组成，二者之间通过一个正弦形材料界面分隔。土壤结构在平面应变条件下建模，并且整个上表面受到脉冲垂直牵引力的作用。牵引力从t?y=0到t?y=1500kN/m2线性增加，在1到2秒的间隔内线性减小至t?y=0。上边界处的孔隙压力被设定为零，相当于大气压力。下边界受到刚性支撑，u?x=u?y=0，而两个侧边界处的法向位移通过施加u?x=0进行水平约束。此外，流体流量也被规定为零，q?=0，沿两侧边界和底部边界都如此。时间积分方案中使用的Newmark参数以及材料参数与第3.1节中介绍的第一个问题中的参数相同，除了两个沙层的杨氏模量和固有渗透率不同。杨氏模量分别为E1=3×10? N/m2和E2=3×10? N/m2（即E1/E2=100），固有渗透率分别为κ1=10?12,m2和κ2=10?1?,m2（即κ1/κ2=0.01）。多孔介质内的初始位移、速度和孔隙压力均被设定为零。时间步长从t=0对数增加到t=50s，共分为1000个增量。图18a展示了包含1087个Q4元素的界面拟合FEM网格，在正弦形材料界面附近对元素进行了局部细化，以准确表示界面几何形状。进一步细化对FEM结果的准确性影响可以忽略不计。图18b显示了由25×25=625个四边形元素组成的结构化IGFEM网格，该网格不遵循材料界面。界面几何形状导致一些四边形切割元素被划分为两个四边形子域（图18d），而其他元素则被划分为一个三角形和一个五边形子域（图18c）。在第一种情况下，可以通过使用具有双线性形状函数的四边形子元素对每个四边形子域进行建模，如第2.3节中详细描述并在图3中示意。然而，当界面将一个四边形父元素划分为一个三角形和一个五边形子域时，需要额外的处理来确保IGFEM近似保持最佳准确性。

图17. 由两个沙层组成的土壤结构的载荷剖面、几何形状和边界条件，两个沙层之间通过一个正弦形材料界面分隔。为了处理这个问题，四边形父元素被细分为四个三角形积分元素，如图18c所示。需要注意的是，从四边形父元素的拉格朗日插值开始，使用3节点三角形积分元素来评估富集函数并不能达到最佳准确性。这是因为使用线性三角形元素构建的富集函数不包含四边形父元素插值中的双线性项。为了解决这个问题，采用了Aragón等人（2020年）提出的方法，即通过使用6节点二次三角形积分元素来提高富集函数的多项式阶数。如图18c所示，二次富集公式在四个三角形积分元素的边缘中点引入了额外的富集节点，用红色实心圆表示。因此，基于3节点积分元素的线性富集导致的准确性损失通过使用基于6节点积分元素的二次富集得到补偿。

图18. 为图17所示的具有正弦形材料界面的土壤结构生成的FEM和IGFEM离散化。(a) 在正弦形材料界面周围进行局部细化的共形FEM网格（红色虚线）。(b) 在材料界面处使用非共形有限元素的结构化IGFEM网格。(c) 一个四边形父元素被划分为一个三角形和一个五边形子域。父元素被划分为四个具有二次富集的6节点三角形积分元素；在边缘中心添加了额外的富集节点。(d) 一个四边形父元素被划分为两个四边形子元素。在(c)和(d)中，黑色和红色实心圆分别代表标准节点和富集节点。

图19a和图19b分别展示了点A（x=12.5m,y=17m）和B（x=12.5m,y=23m）处的孔隙压力p和垂直位移uy随时间的变化（以对数时间尺度显示），其位置见图17。这两种方法都使用了FEM和IGFEM进行计算，结果非常一致。在加载阶段（0< />
图19. 使用标准fem模型（蓝色实线）和igfem模型（红色虚线）计算的图17所示的具有正弦形材料界面的土壤结构的时变响应。(a) 图17中所示的a点和b点的孔隙压力p。(b) 图17中所示的a点和b点的垂直位移uy。

图20. 图17所示的具有正弦形材料界面的土壤结构的场变量等高线图，分别用标准fem模型（左）和igfem模型（右）计算。(a) t= 2秒时的孔隙压力p。(b) t= 2秒时的垂直位移uy。

4. 总结和结论性意见

在本文中，我们介绍了用于求解由分段均匀相组成的可变形、饱和多孔介质的完全耦合水力-力学方程的界面增强广义有限元方法（igfem）。不同、完美结合的均匀相之间的连续界面在固相位移和液相压力场中引入了弱不连续性。在标准fem中，准确表示这些梯度不连续性需要元件边缘与界面对齐，因为不对齐的元件无法自然重现主要变量导数的跳跃。igfem通过在切割元素（被材料界面相交的元素）的近似空间中添加函数来克服这一限制，这些函数能够再现弱不连续性的运动学特性，从而无需网格遵循界面即可准确表示梯度跳跃。耦合水力-力学分析的控制方程包括线性动量平衡方程和饱和多孔介质的储存方程，这些方程使用igfem进行空间离散化。随后的时间离散化采用广义newmark方法进行。所得到的耦合方程组使用基于newton-raphson方法的单块更新方案进行迭代求解。

igfem的一个关键优势是它允许材料界面任意地与元件相交，而不需要网格一致。这种灵活性是通过将富集函数结合到切割元素的近似空间中实现的。对于数值积分，这些切割元素进一步使用拉格朗日形状函数细分为积分元素，然后用于高效计算局部元素贡献。值得注意的是，积分网格也可以被视为标准的fem网格，其中切割的父元素在局部被线性三角形替换。这种方法本质上是noble等人（2010年）提出的共形分解fem（cd-fem）的基础。主要区别在于igfem是层次化构建离散空间的：标准fem的基础保持不变，而在其上添加了富集函数。尽管这牺牲了切割元素中的单位分割属性，但两种公式得出的解是相同的。然而，在这两种方法中，细条（几何退化的）元素可能导致不良条件。cd-fem通过预处理来解决这个问题，而igfem除了允许预处理外，还提供另一种补救措施：当界面接近网格节点时，可以缩放富集函数，从而在不修改底层网格的情况下保持稳定性（aragón等人，2020年）。

与其他富集有限元方法（如xfem/gfem）相比，igfem采用了不同的富集策略，其中富集函数直接与插入不连续性的节点相关联，而不是与标准网格节点相关联。因此，在混合元素中，igfem中的富集函数完全消失，即在靠近切割元素的元素中消失，从而避免了可能降低准确性和减少收敛率的寄生项。 e2=100），大部分变形发生在较软的上层。因此，位于较硬下层的B点比A点经历更大的垂直位移。图20a和图20b分别展示了卸载阶段结束时（t = 2秒）的孔隙水压力p和垂直位移uy的等高线图。可以看出，fem和igfem模拟在整个域内的压力和位移场都非常吻合。压力等高线表明，由于快速卸载和土壤骨架内孔隙流体的延迟重新分布，产生的负孔隙压力在域的上部区域最为明显。此外，位移等高线显示位移幅度随深度减小，大部分变形发生在较软的上层。材料界面的正弦形状进一步导致沿着界面的压力和位移的空间变化，反映了两个沙层之间的局部刚度和渗透率差异。 图19. 使用标准fem模型（蓝色实线）和igfem模型（红色虚线）计算的图17所示的具有正弦形材料界面的土壤结构的时变响应。(a) 图17中所示的a点和b点的孔隙压力p。(b) 图17中所示的a点和b点的垂直位移uy。 图20. 图17所示的具有正弦形材料界面的土壤结构的场变量等高线图，分别用标准fem模型（左）和igfem模型（右）计算。(a) t=2秒时的孔隙压力p。(b) t=2秒时的垂直位移uy。 4. 总结和结论性意见 在本文中，我们介绍了用于求解由分段均匀相组成的可变形、饱和多孔介质的完全耦合水力-力学方程的界面增强广义有限元方法（igfem）。不同、完美结合的均匀相之间的连续界面在固相位移和液相压力场中引入了弱不连续性。在标准fem中，准确表示这些梯度不连续性需要元件边缘与界面对齐，因为不对齐的元件无法自然重现主要变量导数的跳跃。igfem通过在切割元素（被材料界面相交的元素）的近似空间中添加函数来克服这一限制，这些函数能够再现弱不连续性的运动学特性，从而无需网格遵循界面即可准确表示梯度跳跃。耦合水力-力学分析的控制方程包括线性动量平衡方程和饱和多孔介质的储存方程，这些方程使用igfem进行空间离散化。随后的时间离散化采用广义newmark方法进行。所得到的耦合方程组使用基于newton-raphson方法的单块更新方案进行迭代求解。 igfem的一个关键优势是它允许材料界面任意地与元件相交，而不需要网格一致。这种灵活性是通过将富集函数结合到切割元素的近似空间中实现的。对于数值积分，这些切割元素进一步使用拉格朗日形状函数细分为积分元素，然后用于高效计算局部元素贡献。值得注意的是，积分网格也可以被视为标准的fem网格，其中切割的父元素在局部被线性三角形替换。这种方法本质上是noble等人（2010年）提出的共形分解fem（cd-fem）的基础。主要区别在于igfem是层次化构建离散空间的：标准fem的基础保持不变，而在其上添加了富集函数。尽管这牺牲了切割元素中的单位分割属性，但两种公式得出的解是相同的。然而，在这两种方法中，细条（几何退化的）元素可能导致不良条件。cd-fem通过预处理来解决这个问题，而igfem除了允许预处理外，还提供另一种补救措施：当界面接近网格节点时，可以缩放富集函数，从而在不修改底层网格的情况下保持稳定性（aragón等人，2020年）。 与其他富集有限元方法（如xfem>
图19. 使用标准fem模型（蓝色实线）和igfem模型（红色虚线）计算的图17所示的具有正弦形材料界面的土壤结构的时变响应。(a) 图17中所示的a点和b点的孔隙压力p。(b) 图17中所示的a点和b点的垂直位移uy。

图20. 图17所示的具有正弦形材料界面的土壤结构的场变量等高线图，分别用标准fem模型（左）和igfem模型（右）计算。(a) t= 2秒时的孔隙压力p。(b) t= 2秒时的垂直位移uy。

4. 总结和结论性意见

在本文中，我们介绍了用于求解由分段均匀相组成的可变形、饱和多孔介质的完全耦合水力-力学方程的界面增强广义有限元方法（igfem）。不同、完美结合的均匀相之间的连续界面在固相位移和液相压力场中引入了弱不连续性。在标准fem中，准确表示这些梯度不连续性需要元件边缘与界面对齐，因为不对齐的元件无法自然重现主要变量导数的跳跃。igfem通过在切割元素（被材料界面相交的元素）的近似空间中添加函数来克服这一限制，这些函数能够再现弱不连续性的运动学特性，从而无需网格遵循界面即可准确表示梯度跳跃。耦合水力-力学分析的控制方程包括线性动量平衡方程和饱和多孔介质的储存方程，这些方程使用igfem进行空间离散化。随后的时间离散化采用广义newmark方法进行。所得到的耦合方程组使用基于newton-raphson方法的单块更新方案进行迭代求解。

igfem的一个关键优势是它允许材料界面任意地与元件相交，而不需要网格一致。这种灵活性是通过将富集函数结合到切割元素的近似空间中实现的。对于数值积分，这些切割元素进一步使用拉格朗日形状函数细分为积分元素，然后用于高效计算局部元素贡献。值得注意的是，积分网格也可以被视为标准的fem网格，其中切割的父元素在局部被线性三角形替换。这种方法本质上是noble等人（2010年）提出的共形分解fem（cd-fem）的基础。主要区别在于igfem是层次化构建离散空间的：标准fem的基础保持不变，而在其上添加了富集函数。尽管这牺牲了切割元素中的单位分割属性，但两种公式得出的解是相同的。然而，在这两种方法中，细条（几何退化的）元素可能导致不良条件。cd-fem通过预处理来解决这个问题，而igfem除了允许预处理外，还提供另一种补救措施：当界面接近网格节点时，可以缩放富集函数，从而在不修改底层网格的情况下保持稳定性（aragón等人，2020年）。

与其他富集有限元方法（如xfem/gfem）相比，igfem采用了不同的富集策略，其中富集函数直接与插入不连续性的节点相关联，而不是与标准网格节点相关联。因此，在混合元素中，igfem中的富集函数完全消失，即在靠近切割元素的元素中消失，从而避免了可能降低准确性和减少收敛率的寄生项。>在XFEM/GFEM中处理这类问题需要使用更高阶的近似方法或修改过的富集函数来消除寄生项，正如Aragón等人（2010年）所解释的那样。由于界面与网格节点距离过近而导致的病态问题在XFEM/GFEM和IGFEM中都会出现。在XFEM/GFEM中，通过Stable GFEM（SGFEM）和Strongly Stable GFEM（SSGFEM）公式（Kergrene等人，2016年；Zhang等人，2020年）来解决这一问题，这些公式修改了富集函数以减少其对系统条件数的不利影响。相比之下，IGFEM提供了一种更为简单的方法：可以通过使用类似Jacobi的对角预处理器或在界面接近网格节点时缩放富集函数来提高稳定性（Aragón等人，2020年）。在这两种情况下，条件数都与网格尺寸h成Oh^-2的平方关系变化，与标准FEM在拟合网格上的行为一致。在精度方面，IGFEM在几何形状符合的网格上能够达到与标准FEM相同的收敛速度（Soghrati等人，2012年；van den Boom等人，2019年），这对于没有奇异性的问题来说是理想的。即使存在奇异性，只要采用适当的特定问题函数进行近似，XFEM/GFEM也能恢复最优收敛性（Aragón和Duarte，2023年）。与标准FEM相比，IGFEM和XFEM/GFEM所需的CPU开销很小，因为由界面富集的自由度只占总自由度的一小部分。本文提出的四个数值示例展示了当前IGFEM公式的有效性及其解决实际岩土问题的能力。具体来说，考虑的示例包括：（i）分层土柱的一维固结；（ii）水平分层土结构的二维固结；（iii）向分层地下储层注水；（iv）具有正弦材料界面的水平分层土的瞬态响应。在所有四个示例中，IGFEM的结果与使用标准FEM和符合网格的方法得到的结果都非常吻合，这验证了所提出方法的准确性和稳健性。通过系统地细化网格和时间步长，分别证明了在空间和时间上都能得到收敛解。在第一个示例中，数值结果还与Khoi和Haghighat（2011年）使用XFEM/GFEM得到的结果进行了比较。虽然IGFEM提供了准确的解，但XFEM/GFEM方法在位移和压力场中表现出明显的误差。XFEM/GFEM中观察到的不准确性可能来源于混合元素中出现的寄生非线性项。如上所述，IGFEM避免了这个问题，因为富集函数在这些元素中恰好消失。此外，第四个示例涉及正弦材料界面，展示了IGFEM在不需要网格对齐的情况下准确捕捉复杂非平面界面的能力，进一步突显了该方法在处理具有几何复杂边界的多孔介质的实用水力学分析中的灵活性和稳健性。本文提出的模型公式和结果为未来关于具有材料界面的饱和多孔介质的IGFEM建模研究提供了有意义的方向。例如，IGFEM方法可以扩展用于进行具有嵌入材料界面的多孔介质的热-水-力学（THM）分析。在这样的问题中，不仅固体位移和流体压力场，温度场也可能在材料界面处表现出较弱的不连续性。另一个未来的研究方向可以是结合Discontinuity-Enriched Finite Element Method（DE-FEM）（Aragón和Simone，2017年）来扩展当前的IGFEM公式，从而能够分析断裂的异质多孔介质或其断裂传播过程（参见Zhang等人，2025年）。