一种用于冻土中颗粒尺度相变与液态水流的热力学FEM-DEM算法

《Computers and Geotechnics》：A thermodynamic FEM-DEM algorithm for granular-scale phase transition and liquid-water flow in frozen soil

　　

摘要

本研究开发了一种耦合有限元-离散元（FEM-DEM）模型。基于颗粒热力学理论，该模型在COMSOL中求解连续介质热-水-力（THM）耦合过程，并在PFC中模拟颗粒运动和冰胶结演化。场间交互通过体积含冰量、未冻水含量、孔隙水压力和等效导热系数实现。水-冰相变以冰颗粒的生成、融化和再分布来表示。孔隙水压力被转化为等效的颗粒尺度渗流驱动力，同时利用颗粒温度来表征非平衡耗散区。实验验证和数值模拟表明，冻结锋是一个活跃的颗粒结构区。未冻水持续向该区域迁移，为冰晶生长提供主要水源。水力边界效应主要集中在底部附近，对整体冻结演化的影响有限。孔隙水压力对冻胀具有阶段依赖的放大效应。等效颗粒半径主要影响热传导和冻胀速率。所提出的模型为冻土中多尺度相变、水迁移和冻胀机制的捕捉提供了一种数值方法。

引言

土壤冻结是由温度降低驱动的多尺度耦合过程，涉及水迁移、水-冰相变、潜热释放和土壤结构演化。在冻结过程中，孔隙水不仅在原地冻结，还在温度梯度和基质吸力的协同作用下持续向冻结锋迁移（Bai等人，2020，Tang等人，2021）。未冻水含量和土壤冻结特性曲线（SFCC）是连接温度场、水迁移和水-冰相变的核心状态变量，从而构成了描述冻土水热行为的本构基础（Zhou等人，2020，Vu等人，2022，Bai等人，2025b）。在寒冷地区路堤、地基处理和人工冻土等工程应用中，水迁移和冰晶形成改变了土壤的孔隙结构、颗粒间接触和宏观变形，从而影响工程稳定性和服役性能（Yu等人，2024）。在宏观尺度上，总含水率表征了冻结过程中水的整体再分配，而需要分别演化未冻水和冰才能揭示局部相平衡、水源补给通道和活跃冰晶生长区域。水迁移、相变和冰胶结结构形成不仅受宏观温度场或总含水率场的控制，还受颗粒排列、接触状态和局部传热条件的控制（Zhang等人，2024）。因此，对这些耦合过程进行多尺度解释对于理解冻土的力学行为演化至关重要。

冻土中耦合热-水力过程的理论基础已经变得相对系统化。Dall'Amico等人（2011）基于能量方程和Richards方程开发了变饱和度冻土的守恒型热-水力模型，在冻结过程中实现了能量守恒和稳定解。Kurylyk和Watanabe（2013）系统地建立了水迁移和相变过程，阐明了修正Richards方程、热传导方程和相平衡关系在数值模拟中的作用。Painter和Karra（2014）提出了亚零条件下的未冻水含量本构模型，为描述温度、吸力和水-冰相态的分布提供了基础。近年来，非局部热-水力输运模型（Nikolaev等人，2024）、基于热力学的SFCC方程（Li和Vanapalli，2023）、数据驱动的SFCC预测方法（Park等人，2025）、热-水-力-化学耦合框架（Li和Yin，2025）以及冰胶结颗粒介质的多尺度模型（Yu等人，2024）不断涌现，推动研究从连续介质描述向多场和多尺度视角发展。然而，大多数现有模型依赖等效参数来表征孔隙尺度和相变过程，难以捕捉颗粒排列、接触传热、孔隙喉道水流通道和冰胶结结构对冻结行为的控制作用。

与理论进展并行，冻土的数值模拟正从连续介质计算向跨尺度耦合演进。基于FEM的连续介质模型可以求解非饱和冻土中的耦合THM问题（Yin等人，2018），也可以与SFCC曲线结合模拟人工冻土过程中的温度和渗流变化（Li等人，2023）。这些方法适用于工程尺度预测，但无法直接捕捉颗粒尺度上的控制机制。为解决这一局限性，DEM模型已被用于描述颗粒接触、局部变形和热-力响应（Sun等人，2024）。CFD-DEM方法进一步将连续介质流场与颗粒运动耦合，揭示了渗流下的颗粒迁移和孔隙结构演化（Zhang等人，2023）。与上述策略不同，COMSOL可以将颗粒热力学的控制方程嵌入连续介质多物理场框架中，求解颗粒温度和颗粒熵等变量，从而在连续介质尺度上反映局部传热和颗粒重排效应（Bai等人，2025a）。FEM-DEM接口已成功应用于多孔介质的跨尺度建模（Pirnia等人，2019），PFC-COMSOL耦合在煤层水力压裂和渗透率演化分析中也显示出良好的适用性（Wang等人，2022）。因此，将COMSOL的多物理场求解能力与PFC的颗粒尺度表征能力相结合，实现了统一框架内连续介质场与颗粒尺度响应的跨尺度耦合。这为冻土多尺度建模提供了一条具有物理解释力的数值途径。

基于上述分析，本研究开发了一种用于冻土THM过程的非对称耦合FEM-DEM多尺度框架。FEM组件在COMSOL中求解温度场、渗流场、相变过程和土骨架位移。DEM组件在PFC中模拟颗粒运动、接触演化、冰胶结形成和颗粒结构重排。在该模型中，水-冰相变以颗粒间冰胶结的演化来等效表示。孔隙水压力梯度被引入作为颗粒尺度渗流力，以追踪未冻水向冻结锋的方向性补给。颗粒温度用于识别冻结锋附近的局部耗散和结构活化。通过模型验证、单向冻结模拟、水力边界对比以及孔隙压力和等效颗粒尺度的敏感性分析，本研究旨在揭示冻结锋内颗粒尺度下水累积、相变再分布、冰胶结演化和冻胀发展的耦合机制。

章节片段

控制方程

非饱和冻土可以视为由土颗粒（S）、液态水（W）、冰（I）和水蒸气（G）组成的四相系统。颗粒温度最初由Haff（1983）引入，用于量化颗粒相对于局部平均运动的波动动能。在颗粒热力学中，颗粒温度进一步被引入作为颗粒熵的共轭变量，作为描述颗粒尺度运动、结构重排和

水热响应演化的验证

为评估模型在再现冻结过程中温度和未冻水含量耦合演化方面的能力，选取Zhong等人（2025）报道的粉质黏土快速冷却试验进行对比。试样为直径10 cm、高10 cm的圆柱体，初始干密度为1.7 g/cm³，初始体积含水率为0.316。试验前，试样在0.5 °C下平衡。随后环境温度迅速降低

水力边界条件的影响

为阐明宏观场在单向冻结下如何传递到颗粒尺度响应，仅在顶部施加单向冷却，产生自顶向下的冻结锋。与多边界冷却相比，这种配置减少了面内对称性，并产生了更清晰的温度梯度、冰晶累积和孔隙水压力梯度的垂直分层。

为评估底部水力边界对冻结过程的影响，采用无流量边界

结论

本研究开发了一种热力学FEM-DEM算法来模拟土壤冻结过程，捕捉了连续介质尺度THM场与颗粒尺度相变、液态水流动和结构演化之间的相互作用。根据模型验证和单向冻结模拟的结果，主要结论如下：

• (1)
  所建立的热力学FEM-DEM算法实现了连续介质THM场与颗粒尺度响应之间的跨尺度耦合。

CRediT作者贡献声明

Bixia Zhang：撰写 – 初稿，调查研究，形式化分析，数据整理。Bing Bai：撰写 – 审阅与编辑，方法论，调查研究，概念构思。

利益冲突声明

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

致谢

本研究由中央高校基本科研业务费（2026YJS163）和中央高校基本科研业务费（2026JBZD015）资助。

Bixia Zhang|Bing Bai

教育部城市地下工程重点实验室，北京交通大学，北京 100044，中国