一种用于冻土冻胀裂缝热 - 水 - 力建模的混合近场动力学 - 有限元方法

《Computers and Geotechnics》：A hybrid peridynamics-finite element method for thermo-hydro-mechanical modeling of frost cracking in frozen soils

　　

**摘要**

寒冷地区土壤冻裂的模拟涉及多种物理过程的耦合，包括热传导、水分迁移、相变和机械变形，这对传统数值方法构成了巨大挑战。本研究基于经典冻土理论的热-水-力（THM）耦合模型，结合混合邻域动力学-有限元方法（PD-FEM），建立了一个模拟土壤冻结过程中损伤与开裂行为的数值框架。邻域动力学方法用于描述裂纹演化过程，无需预设裂纹路径，而有限元方法则用于求解耦合的温度与水分迁移方程。本研究模拟了单向冻结条件下土柱的冻胀开裂过程。模拟结果与实验数据吻合良好，表明所提出的模型能够有效捕捉冻结过程中土壤的损伤、裂纹萌生与扩展，同时高效地耦合多种物理过程。此外，本研究还探讨了裂纹扩展及裂纹内充填材料对冻结过程的影响。

**引言**

多年冻土和季节性冻土覆盖了北半球约50%的陆地表面（Li and Yin, 2024, Dobinski, 2011）。在这些地区，地表通常在冬季冻结，夏季融化。在冻结和冻融循环作用下，土壤开裂、冻胀和融化沉降现象频繁发生，进而降低土壤的承载力（Chen et al., 2024, Li et al., 2018, Zheng et al., 2024）。其中，冻土开裂不仅破坏土体的结构完整性，还可能改变其力学和水力性质（Cheng et al., 2020）。裂纹的存在为水分迁移和渗漏提供了通道，从而增加水力传导率并降低不透水性（Cherkauer and Lettenmaier, 1999）。此外，开裂可能导致剪切强度降低，使土壤结构和地基更容易发生破坏和承载力丧失（Tang et al., 2021a, Lu et al., 2016, Albrecht and Benson, 2001），导致结构稳定性下降（Fan and Kong, 2012），并加剧崩塌、滑坡和泥石流等侵蚀相关灾害的风险（Goodarzi and Salimi, 2015, Moravej et al., 2018）。因此，深入研究冻结条件下冻土的开裂行为对于阐明冻裂机理、预测结构变形和破坏模式、指导寒冷地区土木工程设计及防护具有重要意义。

当环境温度降至冰点以下时，土壤孔隙水结冰，体积膨胀约9%。在约束条件下，这种体积膨胀将产生显著的冻胀应力。当拉应力超过土壤抗拉强度时，微裂纹形成并逐渐扩展，最终发育成宏观裂纹（Benson and Othman, 1993, Liu et al., 2022, Zhang et al., 1916）。研究人员通过广泛的实验室和野外试验对土壤冻胀机理及其相关模型进行了研究（Koopmans and Miller, 1966, O'Neill, 1983, Konrad, 1989, Konrad and Duquennoi, 1993）。基于这些研究，提出了一系列理论模型，包括毛细管理论（Tsytovich, 1975）、冰透镜分离理论（Miller, 1972）、刚性冰模型（O'Neill and Miller, 2008）和分离潜势模型（Konrad and Morgenstern, 1981），逐步丰富了冻胀理论理解。尽管冻胀机理的理解经历了百余年的发展深化，但尚未形成统一的模型。普遍认为冻胀开裂并非单一过程，而是热-水-力（THM）耦合的复杂行为。

基于上述理论，研究人员采用数值模拟方法来研究冻胀现象，以弥补实验研究的局限性。与实验室试验相比，数值模拟可以避免高昂的人工成本、仪器精度限制以及人为操作干扰的影响，从而为冻胀机理和参数敏感性分析提供更灵活高效的方法（Tang et al., 2021a）。William Harlan是最早建立多孔介质一维热-水（TH）耦合模型的研究者之一，系统研究了土壤冻结和融化过程中温度场和含水率场的演化（Harlan, 1973）。随后，许多研究者发展了各种TH耦合模型（Harris et al., 2009, He et al., 2021, Heinze, 2021, Wang et al., 2011, Li et al., 2019, Tang et al., 2021b, Wang et al., 2023）。基于不同的控制方程组和参数化策略，采用有限元法、有限差分法及相关数值算法系统地模拟和分析了冻土中的水分迁移和热量传递过程。此外，基于Terzaghi（Nash, 1961）和Biot（1956）建立的理论基础，研究人员发展并完善了适用于冻土问题的各种THM耦合模型。Mu和Ladanyi（1987）开发了冻结土壤的耦合THM模型，成功预测了温度、总冻胀量和应力演化。Selvadurai等（2000）开发了一个数值模型来模拟土壤与埋地管道相互作用过程中的冻胀行为。Lai等（2001）基于传热理论、渗流理论和冻土力学基本原理建立了THM耦合分析框架，并利用该框架对寒冷地区冻结条件下现有隧道和支护结构的力学响应进行了数值研究。Yang等（2006）采用半耦合THM方法研究了人工冻土过程中的冻胀机理，但其工作中使用经验关系简化了水分扩散过程。Zhou、Li和Pei建立了包含温度、孔隙率和位移变量的冻胀数学模型，其中采用Clausius-Clapeyron方程作为冰-水相之间的平衡条件（Li et al., 2018, Zhou and Li, 2012）。Li等（2022b）通过解析推导和有限元方法研究了冻结条件下有限均质土层中的THM响应及水力边界条件的影响。研究表明，在表面冰积情况下，孔隙水压力基本保持为负值，由于低温吸力的作用，水压将引起地表消散和沉降，水和热耦合效应下的沉降分析是可靠且必要的。Li等（2026）构建了考虑冰增强效应的THM弹塑性本构模型，有效模拟了冻融条件下土壤的变形行为。三轴试验的数值模拟结果表明，所提出的模型很好地反映了强度增长与冰含量增长之间的典型关系，揭示了强度的衰减变化规律，并进一步证明了冰在土体中的增强机理。然而，现有冻土数值研究主要集中于模拟冻胀变形特征和水-热耦合迁移过程，而关于冻裂萌生和扩展的数值研究相对有限，有待进一步深入探索。

冻土冻裂过程涉及水-冰相变、THM耦合场以及裂纹萌生和扩展等不连续行为。然而，基于偏微分方程的经典理论和数值方法难以处理与裂纹扩展相关的奇异性问题（Zheng et al., 2024, Silling et al., 2007）。Silling（2000）提出的基于积分的非局部邻域动力学（PD）理论及其相关数值方法为分析固体材料和结构的损伤与破坏开辟了新的途径。该理论在处理涉及非局部长程相互作用、多裂纹扩展和多点损伤碎裂的不连续破坏问题方面具有显著优势（Madenci and Oterkus, 2014, Gerstle and Herbert, 2015, Bobaru et al., 2016）。PD理论已广泛应用于模拟不同加载条件下各种材料的传热传质过程及破坏行为（Diehl et al., 2019, Zeleke et al., 2021, Oterkus and Oterkus, 2024, Liu et al., 2025, Dorduncu et al., 2024, Li et al., 2022a, Li et al., 2024），取得了重大研究进展。

近年来，由于其表征裂纹扩展的独特优势，PD理论已被广泛应用于模拟土壤开裂。例如，Menon和Song（2019）、Yan等（2023）、Zhang等（2023）、Qiao等（2025a）、Liu等（2024）、Liu等（2023a）以及Qiao等（2025b）均利用PD方法研究了土壤的干缩开裂，取得了良好效果。关于土壤冻裂纹，Zheng等（2024）提出了一种考虑材料相变的耦合PD模型，分析了单向冻结条件下土壤冻裂纹的演化。然而，该模型未考虑水分迁移机制，也无法反映冻裂过程中THM的耦合效应。鉴于PD在模拟冻土裂纹萌生和扩展方面的显著优势，该方法为冻土裂纹行为的数值分析提供了有效途径，具有广阔的应用前景。

基于Li等（2018）、Li等（2014a）和Li等（2014b）建立的冻土THM耦合模型，本研究提出了一种邻域动力学-有限元方法（PD-FEM）混合建模方法，用于模拟冻土中冻胀引起的开裂行为。与完全基于局部连续介质力学的传统THM方法不同，本研究采用PD方法描述冻土的变形和损伤演化，同时利用有限元方法离散和求解温度及水分迁移方程。这一混合策略将PD在模拟损伤和裂纹扩展方面的优势与有限元在多物理场计算方面的效率和成熟性相结合，在统一的计算框架内一致地处理相变、水分迁移和机械变形过程。该方法能够自然地模拟裂纹的萌生、扩展和相互作用，而无需额外的裂纹追踪算法或网格重划分技术。此外，利用COMSOL Multiphysics与MATLAB的二次开发接口，实现了混合建模方案的数值求解与耦合计算过程，显著提高了THM问题的计算效率。在此基础上，对冻土冻胀开裂过程进行了系统数值模拟。

本文第二节介绍了经典冻土THM耦合模型和PD基本理论。第三节描述了PD-FEM混合建模方法的实现过程。第四节给出了土壤冻裂的计算模型及具体参数设置。第五节通过模拟单向冻结条件下土柱的冻结过程验证了模型和方法的准确性，分析了裂纹萌生与扩展过程以及裂纹充填介质对冻胀开裂的影响。最后，第六节对研究进行了总结与展望。

**热传导模型**

土壤冻结过程中，由于孔隙水流动速度相对缓慢，对流换热对温度分布的影响可以忽略不计，仅考虑热传导。根据傅里叶定律，将相变潜热视为热源。冻土的热传导微分方程为：

ρC(?T/?t) = λ?2T + L·ρ_I(?θ_I/?t)

其中，T表示土壤的瞬态温度，ρ和ρ_I分别为土壤和冰的密度，

**PD-FEM混合建模**

由于PD的非局部相互作用特性，其计算效率通常低于传统有限元方法（FEM），特别是在土壤冻结等长期加载过程中尤为明显。

为平衡PD在表征裂纹萌生与扩展行为方面的优势以及FEM在求解效率上的优势，本研究采用PD-FEM混合建模策略对土壤冻胀开裂计算模型中的THM耦合问题进行数值分析。本研究采用第2节提出的THM耦合模型以及第3节开发的PD-FEM混合建模方法和算法，对单向冻结条件下黏土土柱中的裂纹扩展过程进行了数值模拟（见图5），与Li等人（2018）开展的实验相一致。该模拟旨在验证模型和方法在模拟冻土开裂方面的准确性和有效性。详细实验信息请参阅Li的学位论文（Li等人，2018）。

温度分析

为进一步研究土壤冻结过程中的温度场演化并验证本研究中数值模型的可靠性，将四个特征位置T1–T4（如图6所示）的温度–时间曲线与Li等人（2018）的实验数据进行了对比分析。结果如图8所示。需要注意的是，正如Li等人所述，实验过程中土柱的初始温度未达到……

结论

寒冷地区土壤的冻胀开裂过程极为复杂，其数值模拟涉及水分与热量的迁移、冰的形成与积聚、THM耦合以及裂纹扩展，所有这些都涉及不连续行为。这对现有的计算模型和数值方法构成了重大挑战。本研究基于多孔介质理论中的冻土胀裂模型，开发了一种PD-FEM混合模型，并提出了一种……

作者贡献声明：

魏晓丽：撰写原稿、可视化、方法论。周志伟：监督指导、资源提供、调查研究。顾鑫：监督指导、方法论、概念构思。李星：撰写原稿、方法论、形式分析。张青：项目管理、资金获取、概念构思。

利益冲突声明：

作者声明不存在已知的竞争性财务利益或个人关系，这些关系可能影响本文所报告的工作。

致谢

本研究由国家重点研发计划（2022YFC2903903）、河西学院博士研究基金（项目编号：KYQD2025037）、甘肃省教育厅高校教师创新基金项目（项目编号：2025B-192）以及国家自然科学基金（编号：12472200、42271153、11932006）资助。

魏晓丽 | 周志伟 | 顾鑫 | 李星 | 张青

河西学院土木工程学院，张掖，734000，中国