冻融循环威胁高海拔寒冷地区岩石斜坡的稳定性：SPH模拟技术的进步

《INTERNATIONAL JOURNAL OF ROCK MECHANICS AND MINING SCIENCES》：Freeze-thaw cycles threaten high-altitude cold-region rock slopes stability: Advances in SPH simulation

　　成志夏|沈振明|李航|彭明|张立民
香港科技大学土木与环境工程系，沿海城市气候韧性国家重点实验室，999077，香港

**摘要**
在全球变暖的背景下，寒冷地区日益严重的冻融循环引发了频繁的岩石滑坡，对基础设施安全构成了重大威胁。本文提出了一种以工程为导向的热-机械-损伤耦合平滑粒子流体动力学（SPH）方法，用于模拟冻融引起的岩石滑坡的渐进性破坏过程，该方法更注重实际应用而非理论复杂性。通过将相变与多孔介质热传递相结合，分析了断裂岩体内的冻融演变过程。冻胀压力和热应力被纳入柯西应力方程中，有效耦合了岩石裂隙与完整基质之间的冻胀和融化沉降现象。所提出的方法经过严格验证，并应用于研究青藏高原不同工程条件下的节理岩坡破坏机制：冻融循环次数、温度范围、节理倾角、冻胀系数（βf）和弹性模量（E）。结果表明，改进的SPH框架成功捕捉了完整的破坏过程，包括热融化弱化、冻胀开裂、冰楔传播、滑动倾倒破坏以及大变形。节理倾角、温度幅度、冻融循环模式、βf和E共同控制着破坏模式的转变、不稳定时间、损伤程度和滑动幅度。值得注意的是，矩形冻融循环造成的损伤最为严重，而βf和E的协同效应与坡体位移和损伤强度呈线性相关。这项工作克服了传统数值方法在模拟冻融岩石滑坡机制方面的关键局限性，为关键寒冷地区工程项目中的地质灾害风险评估和灾害缓解提供了可靠且实用的工具。

**引言**
在全球气候变化的背景下，高海拔寒冷地区的岩坡面临着日益严重的冻融循环威胁。受昼夜温度变化和季节性温度波动的影响，岩体内的孔隙水经历反复的相变，导致复杂的物理-力学效应，如冻胀和断裂损伤。典型案例包括中国哈尔滨-大连高铁沿线切割坡体的冻融崩塌、中国青海门源县山区公路沿线岩石坡体的冻融崩塌，以及青藏高原益公藏布河沿线的本都奥山滑坡。通过水和冰的相变，冻融循环驱动岩体内的岩-水-冰多相多场耦合，导致岩石强度显著下降，安全系数降低。这些现象不仅威胁当地社区的生命和财产安全，还对寒冷地区的生态环境造成不可逆的破坏。然而，现有研究缺乏对岩坡冻融损伤内在机制的系统理解，以及冻融恶化→坡体开裂→岩体损伤→岩石滑坡灾害过程的了解。此外，从工程角度来看，热应力与冻胀对岩体行为的耦合效应尚未得到充分阐明。

**领域规模的冰川相关岩石滑坡**
以白格、 Sedongpu和益公事件为代表的冰川相关岩石滑坡，已通过数值模拟进行了广泛研究，以表征其破坏行为。对于白格滑坡，众多研究详细考察了其破坏后的动态、沉积特征和运动参数，深入了解了其崩塌后的演变过程。在Sedongpu滑坡案例中，数值模型主要关注冰川内的热-水-力（THM）耦合过程及随后碎屑流的动态传播，而控制岩体破坏的潜在启动机制通常被简化处理。至于益公滑坡，基于连续介质、非连续介质和混合数值方法均成功再现了其高速长距离运动模式。然而，对这些研究而言，关键的逐步破坏前演变阶段（即岩体在循环冻融作用下的渐进性退化过程）并未得到足够重视。这些数值研究的一个普遍局限是，它们对循环冻融作用下岩-冰系统逐步破坏演变的关注不足。具体而言，从水-冰相变引发的岩体完整性初始退化，到冻胀引起的裂纹起始和传播，直至最终发生灾难性坡体破坏的耦合热-机械-损伤（TMD）演变过程，在当前的数值建模框架中仍未得到充分解决。

**平滑粒子流体动力学（SPH）**
平滑粒子流体动力学（SPH）是一种拉格朗日无网格数值方法。其核心思想是摒弃传统基于网格的方法对计算域的依赖性，而是将连续介质离散成大量携带物理属性（如质量、密度、速度）的粒子，以数值求解连续介质的控制方程。在温度场模拟方面，SPH已应用于多个研究，例如高精度瞬态热传导问题模拟、大应变热弹性和热粘塑性模拟。在岩石断裂模拟领域，SPH通过使用弹塑性断裂模型来模拟岩体的渐进性破坏过程。SPH方法中提出了伪弹簧策略，通过虚拟弹簧追踪岩石断裂的传播路径，实现动态断裂的介观尺度表征；Zhou和Yu基于改进的核函数方法实现了岩石断裂传播模拟；Xia开发了考虑残余强度的改进损伤框架来模拟岩石断裂传播；Ma利用基于统一孪生剪切强度准则的弹塑性损伤模型预测岩石断裂。然而，这些研究均未考虑岩体断裂后的接触行为。对于大规模岩石滑坡，定义明确的断裂后接触行为对于有效捕捉块体分离、摩擦滑动和倾倒等不连续行为至关重要。因此，Xia开发了一种非连续SPH方法，有效模拟了岩体从断裂到大变形接触的整个过程。在此基础上，为了解决断裂传播与岩体温度场演变之间的强耦合问题，提出了总拉格朗日SPH（TLSPH）耦合模型，以有效解决热断裂过程中的拉伸不稳定性。同时，Yu使用改进的SPH算法模拟了岩体中的热-水-力-损伤耦合问题。为了模拟冻融条件下的坡体渐进性破坏，Bui率先引入了基于固相应力的关键状态弹塑性粘土-砂模型（CASM-F），该模型专门用于捕捉受冻融循环影响的冻结土壤的状态转换机制和大变形特性。Yu提出了基于不可压缩SPH理论的数值框架，将断裂标记嵌入热应力耦合方程中，从而准确模拟由反复冻融循环引起的岩体退化。

**当前对冻融岩体整个演变过程的理解**
目前对冻融岩体从连续变形到冻融和冻胀效应下的断裂起始和块体破碎的整个演变过程的理解仍缺乏理论基础。特别是以下挑战尚未解决：（1）使用SPH方法模拟冻融相变的热力学方程缺乏严格的理论推导，水-冰相变过程中的热应力和冻胀特性及其对岩坡稳定性的影响需要量化。（2）尚未建立多场耦合机制和水-冰相变过程的统一框架，现有模型难以准确描述冻胀、损伤、位移和接触之间的相互作用。（3）使用SPH对冻融恶化→坡体开裂→岩体损伤→岩石滑坡灾害过程进行定性和定量评估的方法尚不明确。

**为解决上述研究空白，我们提出了一种统一的相变热-机械SPH（PT-TM-SPH）框架**
该框架用于模拟岩体在冻融恶化条件下的整个破坏过程，从初始裂纹到大变形。通过冻胀和热应力的耦合，可以揭示岩石滑坡渐进性破坏过程中的水-冰-岩相互作用。为了捕捉水-冰相变-损伤的相互耦合效应，将伪弹簧裂纹策略和具有拉伸截止值的Mohr-Coulomb准则结合到弹性-脆性本构模型中，从而实现裂纹起始和传播的模拟。使用块体-块体接触和库仑摩擦定律捕捉破坏后的动态接触行为变化。本文通过温度模型实验和冻融冻胀模拟实验验证了所提算法的准确性和鲁棒性。此外，还利用泰顺县山区公路沿线岩石坡体的冻融崩塌案例进行了模拟验证，并考虑了温度参数和强度参数变化对冻融滑坡的起始机制、破坏模式和运动机制的影响。同时，评估了冻胀系数、弹性模量和摩擦系数对冻融冻胀响应的影响。最后，讨论了所提模拟方法的优点和缺点，并对相关参数进行了敏感性分析。

**计算框架和工作流程**
增强的冻融-应力-损伤SPH模拟框架包括六个核心模块：输入、热传导、冻融-应力耦合、损伤和裂纹、接触以及输出。计算工作流程如图3所示。
（1）输入模块：导入岩体粒子的初始数据，包括位置、速度、材料参数、初始应力和边界条件。
（2）热传导模块：计算基质和裂隙中粒子的温度。

**结论**
本研究开发了PT-TM-SPH方法，用于模拟青藏高原冻融岩石坡体的变形和破坏过程。此外，还研究了各种参数（包括节理倾角、循环温度幅度、循环模式、粒子分辨率、核半径、冰的摩擦系数、弹性模量和冻胀系数）对模拟结果的影响。主要结论如下：

**作者贡献声明**
成志夏：撰写——原始草案、方法论、概念化；
沈振明：撰写——审阅与编辑、资金获取；
李航：验证、数据管理；
彭明：资源获取、调查、形式分析；
张立民：撰写——审阅与编辑、监督、资金获取。

**利益冲突声明**
作者声明没有已知的竞争性财务利益或个人关系可能影响本文所述的工作。

**致谢**
本研究得到了国家自然科学基金博士学者项目（项目编号424B2055）、国家重点研发计划（项目编号2023YFC3008300）、NSFC/RGC联合研究计划（N_HKUST620/20）以及香港特别行政区政府研究资助委员会（16206923）的财政支持。