COMSOL+XFEM破壁之作：对不连续多孔介质中流体流动的热-流-力学建模
一、摘要
本文针对多场耦合的裂隙多孔介质问题，提出并实现了基于扩展有限元法（XFEM）的数值模拟框架，并将其集成于商业软件COMSOL Multiphysics中。研究背景源于地质工程领域对裂隙、断层等不连续结构在热-流-力多场耦合下的精确建模需求。方法上，作者针对COMSOL平台，设计了专属的富集策略和物理界面修改，实现了多场物理过程（力学变形、孔隙流动、热传导）的统一耦合，并利用高阶Gauss积分与特殊变量处理提升计算效率与精度。通过多组2D/3D数值算例（如双材料梁的温度致裂、含边裂纹板的热力分析、双夹持梁的热力接触、坝体中板桩的渗流分析及多断层岩体模拟），验证了算法的准确性与稳健性，结果与理论或实验数据高度吻合。结论指出，该XFEM框架具备处理复杂多场多断层问题的灵活性与高精度，显著降低了自编代码门槛，促进了地质工程中多物理场裂隙问题的实用化高效模拟，对相关领域的数值研究具有重要推动作用。
二、研究背景及意义
该论文旨在填补多场耦合下可变形多孔介质中不连续界面建模的研究空白，基于Biot理论，提出并实现了与COMSOL结构兼容的XFEM热-水-力多场耦合数值框架，利用专属加密策略与MATLAB联动，实现界面追踪、物理场耦合及精确积分，显著提升了商业软件针对复杂裂隙、多物理场问题的通用性与效率，并通过多组二维与三维算例验证了方法的鲁棒性和准确性。
三、创新点
内容描述：作者首次将扩展有限元方法（XFEM）全面集成到通用多物理场仿真平台COMSOL Multiphysics中，实现了对变形多孔介质中热-流-力（THM）多场耦合问题的高效建模。该框架通过标准和富集物理接口的耦合，结合MATLAB外部函数（LiveLink），支持裂缝等强弱间断的独立建模，突破了以往XFEM仅限于自研代码或有限功能商用软件的局限。 重要性：极大提升了多场裂缝和断层问题的建模灵活性，显著拓展了XFEM在工程实际和大规模复杂多物理问题中的应用范围。
内容描述：提出了一种专用于COMSOL结构的富集区处理策略，利用点约束和变量插值方法，有效消除无关自由度，保证富集区与标准区的物理变量一致性（如应力、流速、热通量等），并通过MATLAB函数实现富集区域的自动识别和更新，支持裂缝扩展与复杂间断的动态追踪。
重要性：显著提升了计算效率，避免了因富集DOF错误定义导致的数值不稳定，为多场耦合裂缝模拟的准确性和稳定性奠定了基础。
图2 裂隙积分转换原理
内容描述：作者开发了基于COMSOL内置圆路径积分和域积分算子的裂尖应力强度因子（SIF）计算方法，并将高阶Gauss积分与富集物理接口结合，实现了多场（热、流、力）耦合环境下的高精度数值集成。同时提出了裂尖区辅助位移场与实际位移场的交互积分方法，提升了裂缝扩展判据和SIF计算的准确性。 重要性：推动了多物理场裂缝模拟的数值精度与理论一致性，为复杂地质、岩石等多场环境下的断层动态分析提供了可靠工具。
四、实验过程
实验一：双材料梁的温度诱导裂纹扩展实验
实验指标：观察温度变化引起的双材料（硼硅玻璃/钢）梁试样中的裂纹扩展路径，并与实验结果进行比对。
实验步骤：
原材料/试剂：AISI304不锈钢条与硼硅玻璃条，分别为3.17 mm×300 mm×6 mm和31.7 mm×300 mm×6 mm，通过胶合连接，形成复合梁。材料的力学与热学性质在文献中给出。
仪器/设备：COMSOL Multiphysics 软件，实验中使用14,432个有限元，平面应力条件，裂纹区平均单元大小0.5 mm。
具体操作步骤：
在复合梁的玻璃部分距左端30 mm处预制初始缺口。 对试样整体施加室温（20°C）降温，实现均匀温度降低。 通过差异热膨胀系数产生纯弯曲，触发裂纹在玻璃部分起始并向双材料界面扩展。 使用最大环向应力准则（t =80 MPa）进行准静态裂纹扩展分析，裂纹增长步长为2.5 mm，初始裂纹长度11 mm。 利用COMSOL的“Auxiliary Sweep”功能进行裂纹增长仿真。
实验结论：裂纹从缺口处起始并向钢/玻璃界面扩展，最终裂纹路径与实验观察高度一致，裂纹趋向沿界面平行方向传播。仿真结果准确反映了温度诱导裂纹扩展的实际机制，表明XFEM和COMSOL的联合应用能够精确模拟此类复杂过程[9][10]。
图3双材料梁中的裂纹扩展；（a）基于𝜓的数值裂纹轨迹，（b）格鲁齐克和里迪（2018 年）的实验裂纹形态
实验二：带边裂纹矩形板的热-力耦合分析
实验指标：计算并验证带有绝热边裂纹的矩形板在热载荷下的规范化应力强度因子（SIF）与裂尖温度分布。
实验步骤：
原材料/试剂：二维矩形板（2.0 m ×0.5 m），半宽绝热边裂纹，材料热-力学参数参照文献。
仪器/设备：COMSOL Multiphysics，均匀结构网格160 ×40单元。
具体操作步骤：
板左右两端分别施加+50°C和-50°C的温度边界条件。 进行总时长100 s的热-力耦合仿真。 计算规范化SIF随时间的变化，并与解析解进行对比。 提取裂尖处温度分布，并与理论公式进行比对。
实验结论：在20 s后，SIF达到稳态值0.491，与解析解0.495相对误差小于1%；裂尖温度分布与理论公式高度吻合。表明提出的XFEM策略在模拟非绝热裂纹环境下的应力和温度场具有极高精度[10]。
图4a - SIF时变曲线；图4b -裂尖温度分布
实验三：双端固定梁热-力接触分析
实验指标：分析含中部贯穿裂纹的双端固定梁在热-力复合载荷下的接触行为及温度场连续性。
实验步骤：
原材料/试剂：3 m×10 m梁，中心有全深度垂直裂纹。材料参数详见文献。
仪器/设备：COMSOL Multiphysics，结构网格155 ×47个四边形单元。
具体操作步骤：
左端固定，右端施加6 cm水平位移，左/右端分别施加0°C和10°C温度。 梁顶边中心施加240 MN均匀分布的竖向载荷。 应用无摩擦的法向接触约束（惩罚参数10^9 MN/m^3）。 在接触区采用热接触导热系数h_cont =10^5 W/m•C，保证温度场连续。
实验结论：裂纹导致梁分为两段，接触区温度场连续，热-力耦合行为被精确捕捉。仿真结果与理论预期一致，验证了XFEM框架在复杂接触与热传导场景下的适用性[10]。
图5温度/位移场演化
图6裂隙面温度剖面
实验四：二维孔隙介质中倾斜/多裂隙的热-流-力耦合分析
实验指标：探究不可渗透断层（单个或多个）对孔隙介质中压力、温度场和流动路径的影响。
实验步骤：
原材料/试剂：1 m ×1 m二维孔隙块体，含长0.3 m、不同倾角的不可渗透断层；多裂隙实验含11条长度0.2 m同尺寸断层。材料参数详见文献。
仪器/设备：COMSOL Multiphysics，结构网格91 ×91单元。
具体操作步骤：
底边施加50°C温度和10^-4 m/s法向流体通量；顶边为排水条件，压力和温度均为零。 对不同断层倾角（0°、30°、45°、60°）进行仿真，观察温度跳跃和压力分布。 多裂隙实验仿真时间8000 s，分析终态压力与温度场分布，提取热流与流体流线分布。
实验结论：不可渗透裂隙对压力和温度场均有显著屏障作用，断层倾角越小（越水平）温度跳跃越大。在多裂隙场景下，热流和流体流动被完全阻断于裂隙处，热流集中于裂尖，表现出强烈不连续性。XFEM强不连续性富集策略表现优异[12]。
图7压力/温度场
实验五：三维孔隙介质中圆盘状断层的热-流-力耦合分析
实验指标：评价三维孔隙介质中不可渗透圆盘状断层对位移、压力、温度场及流动路径的影响。
实验步骤：
原材料/试剂：边长50 m的立方体域，中心含直径20 m的圆盘状不可渗透断层。材料参数同前述实验。
仪器/设备：COMSOL Multiphysics，域内用2,312个砖元素（断层附近，平均单元1.5 m）及45,415个四面体元素（其他区域）进行离散。
具体操作步骤：
底面施加50°C温度与10^-4 m/s流体流量，顶面压力和温度均为零，其余面为无流/无热边界。 除顶面外各面法向约束，模拟原位边界条件。 总仿真时长3 ×10^5 s，提取断层附近的垂直位移、压力和温度场分布。 分析流体与热流线在断层附近的偏转与分布变化。
实验结论：不可渗透圆盘断层在三维场景下同样引起显著的场分布不连续，流体与热流线在断层附近发生偏转，断层区表现出明显阻断作用。XFEM框架可灵活应对二维/三维复杂多场耦合断层问题，具有极高的通用性和准确性[13]。
图8三维位移/压力/温度场及流线
实验六：XFEM在COMSOL中的实现与数值方法验证
实验指标：开发并验证XFEM在COMSOL中处理多场（热-流-力）问题的能力，包括裂纹预处理、富集区识别、数值积分与求解策略。
实验步骤：
原材料/试剂：各类含裂纹或断层的孔隙介质试样，具体结构见前述各实验。
仪器/设备：COMSOL Multiphysics与MATLAB（通过联合仿真接口），用于预处理和富集区识别。
具体操作步骤：
通过FEM初步网格划分并导出节点信息，使用MATLAB脚本进行Level-set分析，识别裂纹/断层区。 用Heaviside函数对富集区进行强不连续性处理，更新物理接口相关变量 采用“Solid Mechanics”、“Darcy’s Law”、“Heat Transfer in Porous Media”等模块分别处理结构变形、流体流动和热传导。 通过20 ×20 Gauss积分对富集元素数值积分，提高计算精度。 采用Generalized alpha时域积分与全隐式耦合求解器处理非线性耦合方程。
实验结论：XFEM与COMSOL结合可实现无网格对齐的裂纹/断层多场耦合模拟，数值精度高，且富集区处理灵活。该方法能有效模拟多种复杂实际地质场景，无需自编代码即可应用于工程问题，具有极高的工业与科研推广价值[1][5][6][9]。
图9热流与流线分布
五、所有实验结论汇总
XFEM在COMSOL中的实现能够高精度模拟温度、流体与力学场耦合下的裂纹或断层问题，无需网格与断层对齐，极大提高了建模效率与适用性[1][9][13]。 各类裂纹/断层对压力场、温度场及流体流动具有显著屏障效应，断层几何与分布显著影响场分布和流动路径[12][13]。 XFEM结合Heaviside强不连续性富集策略能精确捕捉裂尖行为、场跳跃及热流/流线集中现象[5][12]。 多场耦合下的接触行为和温度连续性可通过合理的物理参数设定与富集策略进行准确模拟[10]。 通过COMSOL与MATLAB集成，可灵活实现裂纹/断层预处理、富集区识别与物理接口耦合，方法简便，通用性强[6][1]。
综上，本文所述各项实验系统验证了XFEM在多场耦合复杂断层/裂纹模拟中的卓越表现，具有极高的工程应用与科学研究价值。
六、参考文献：
[1] Jafari, A., Vahab, M., Broumand, P., & Khalili, N. (2023). An eXtended finite element method implementation in COMSOL multiphysics: Thermo-hydro-mechanical modeling of fluid flow in discontinuous porous media. Computers and Geotechnics, 159, 105458.
[2] Khoei, A. R. (2014). Extended Finite Element Method: Theory and Applications. John Wiley & Sons.
[3] Grutzik, S. J., & Reedy, E. D. (2018). Crack path selection in thermally loaded borosilicate/steel bibeam specimen. Experimental Mechanics, 58(1), 1–10.
[4] Duflot, M. (2008). The extended finite element method in thermoelastic fracture mechanics. International Journal for Numerical Methods in Engineering, 74(5), 827–847.
[5] Khoei, A., & Bahmani, B. (2018). Application of an enriched FEM technique in thermo-mechanical contact problems. Computational Mechanics, 62(5), 1127–1154.
[6] Vahab, M., Hirmand, M., & Jafari, A. (2021). Numerical an alysis of multiple hydro-fracture growth in layered media based on a non-differentiable energy minimization approach. Engineering Fracture Mechanics, 241, 107361.
[7] Salimzadeh, S., Paluszny, A., Nick, H. M., & Zimmerman, R. W. (2018). A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems. Geothermics, 71, 212–224.
[8] Borja, R. I. (2008). Assumed enhanced strain and the extended finite element methods: A unification of concepts. Computer Methods in Applied Mechanics and Engineering, 197(33–40), 2789–2803.
[9] Khalili, N., & Selvadurai, A. (2003). A fully coupled constitutive model for thermo-hydro-mechanical a nalysis in elastic media with double porosity. Geophysical Research Letters, 30(24).
[10] Biot, M. A. (1956). General solutions of the equations of elasticity and consolidation for a porous material. Journal of Applied Mechanics, 23(1), 91–96.
[11] Dolbow, J., Moës, N., Belytschko, T. (2000). Discontinuous enrichment in finite elements with a partition of unity method. Finite Elem. A nal. Des. 36 (3–4), 235–260.
[12] Duflot, M. (2008). The extended finite element method in thermoelastic fracture mechanics. Internat. J. Numer. Methods Engrg. 74 (5), 827–847.
[13] Flemisch, B., Berre, I., Boon, W., Fumagalli, A., Schwenck, N., Scotti, A., Stefansson, I., Tatomir, A. (2018). Benchmarks for single-phase flow in fractured porous media. Adv. Water Resour. 111, 239–258.
七、内容来源声明
l 本文核心研究成果源自：Jafari, A., Vahab, M., Broumand, P., & Khalili, N. (2023). Computers and Geotechnics, 159, 105458.（DOI: https://doi.org/10.1016/j.compgeo.2023.105458）
l 本文所述模型、算法及结论仅代表原研究团队观点，公 众号平台不承担其工程应用的直接责任
l 读者参考本文内容进行的科研或实践项目，需自行承担技术风险
