寒区块石路基水热力耦合与冻胀变形的COMSOL模拟要诀

一提到寒区路基，干这行的人第一反应通常都是那六个字——冻胀、融沉、翻浆。我前两年接了一个某寒区道路工程的数值评估任务，核心就是把一段块石路基的水热力耦合行为在COMSOL里算清楚，最后落到冻胀变形预测上。这类项目看着是纯仿真，实际跑起来你会发现，真正难的从来不是建模，而是把冻土物理机制翻译成COMSOL能算的东西。这篇文章就围绕我当时踩过的坑、试过的方法、最终能出可靠结果的操作流程展开，适合正在做冻土路基温度场、水分场或冻胀变形模拟的同行参考，也适合刚接触COMSOL多物理场耦合、想找个真实工程案例练手的读者。

1. 块石路基为什么会冻胀？先搞清机理再谈仿真

很多教程上来就教建模，我觉得这是个误区。冻土区的仿真，物理机理没理顺，后面所有参数和耦合关系都是在瞎填。所以先把块石路基的工作机制和冻胀的形成过程讲透，后面你在COMSOL里的每一个设置才会知道是为什么。

1.1 块石层的"热半导体"效应是如何起作用的

块石路基的基本构造，就是路基面以下铺一层大粒径块石，厚度通常在1米左右。它抗冻的核心机理不是什么保温，恰恰相反，是靠冬天主动"降温"。

块石层里孔隙率很大，通常做到0.35到0.45，块石之间形成了大量连通的大孔隙。冬天大气温度远低于土体温度，块石层内的空气因为顶部受冷变重而下沉，底部相对温暖的空气变轻上升，形成一种自然对流循环。这个循环非常高效——它把地表冷量源源不断地"泵"进路基深处，让基底土体温度持续降低，从而抬高冻结面、保护下卧冻土。

夏天反过来，地表温度高于路基内部温度，空气下冷上热形成稳定层结，对流被抑制，热量主要靠传导缓慢向下传递，进入土体的热量就远小于冬天被"泵"走的热量。这一正一反，就产生了类似半导体整流那样的单向导热特性，工程上把它叫块石层的"热半导体效应"。

搞清这个机理对仿真非常关键：如果你在模型里只给了块石的导热系数，没有考虑冬季对流传热带来的等效导热增强，那一整年的温度场算出来都会偏保守，块石层的作用根本表现不出来。所以模型里必须显式考虑块石层的孔隙对流，或者用等效导热系数随温差变化的方式来覆盖这个效应。

1.2 单算温度场的问题出在哪

常见的错误做法是只做温度场模拟，算完给出冻结深度就完事。这在粗看宏观季节变化时凑合能用，但一落到冻胀变形就必须把水加进来，原因有两方面。

首先是相变潜热。水在0°C附近结冰要释放大量热量——1立方米水完全冻结大约释放3.34×10^5千焦。假设某层土体体积含水率20%，那每立方米土体完全冻结释放的潜热约6.7万千焦，这个数字比土体本身降温释放的显热高出一个量级。潜热会显著拖慢冻结锋面的推进速度，还会让0°C附近出现一个明显的"温度平台期"。不做水热耦合的话，冻结深度推进速率是算不准的，有时甚至会偏差几十厘米。

其次是水分迁移。土体冻结过程中，冻结锋面附近的未冻水会向冷端迁移、在锋面附近富集并冻结成冰，形成分凝冰。这个过程的直接后果是：某个位置的总含水率会从初始20%左右升高到30%甚至更高，冻胀量随之剧增。这正好解释了为什么野外路基冻胀常集中在特定深度而不是均匀分布，也解释了为什么有的路段表面看着没变化、底下却已经拱起一大块。只算温度场，永远解释不了这种局部剧烈变形。

1.3 水热力三场的相互作用链

冻土路基里，水、热、力三场是强耦合关系，我习惯把它画成一条因果链：

温度变化引起水相变和水分迁移，这是热到水的传递；水分迁移和相变改变土体的热容与导热系数、释放或消耗潜热，这是水到热的反馈；温度变化和水分重分布改变土体的应力状态、诱发冻胀应变，这是热和水到力的作用；而应力状态变化又反过来改变土体的孔隙比和渗透系数，这是力到水的反馈，虽然这一项在多数工程模拟里量级较小、常被忽略，但严格做研究时不能完全不管。

在COMSOL里实现这条链，就是接下来要说的核心内容。先记住结论：三场耦合是全链路的，任何一环砍掉，最终冻胀变形结果都不完整。

2. COMSOL里怎么把三场物理过程搭起来

模型搭建是纯操作，但里面每个选择背后都有原因。我按当时项目的推进顺序来写，你可以直接照着排布，能省很多试错时间。

2.1 几何简化：二维断面先跑通，三维再扩展

我的建议是先用二维平面应变模型把方法链跑通，不要一上来就建三维。冻土路基模拟本质是瞬态非线性多场耦合，三维模型自由度动辄上百万，单次瞬态计算可能跑十几个小时还不收敛，调试阶段效率极低。二维模型自由度少一个量级，修改参数后几分钟就能看到结果，非常适合做机理分析和参数敏感性研究。

几何上用一个典型路基横断面：地基土往下取10到12米（保证底部恒温边界的影响不会干扰路基范围），路基填土按设计断面形状，块石层夹在路面结构与路基填土之间。层与层之间用共形网格，先在几何建模里画矩形和梯形，再通过"布尔分割"切出不同材料域。块石层内的大孔隙结构不需要逐块画出，几何上按均质多孔介质层处理，孔隙效应通过等效参数体现。

2.2 物理场接口与多物理场耦合入口的选择

COMSOL里对应三场的接口我建议这样选：

• 温度场用"固体传热"接口，或"多孔介质传热"接口。两者的差别在于后者自带了多孔介质有效导热系数的体积平均算法，对有孔隙率的块石层和土体更合适。
• 水分场用"达西定律"接口描述孔隙水压力驱动的渗流。冻土模拟中饱和-非饱和状态切换很常见，达西接口加"孔隙水压力"因变量，配合非饱和渗透系数方程就能覆盖。
• 力学场用"固体力学"接口，因变量是位移。块石层和土体都按弹性或弹塑性材料处理，冻胀应变作为初始应变或外部应变叠加。

耦合关系不需要你手动写复杂方程，大部分通过多物理场节点完成。我实际项目里用到的耦合设置大概如下：

这里有个细节：COMSOL的"多物理场"节点里并没有一个叫"冻土水热力耦合"的一键选项，它给你的是"固体传热+达西定律"的双向耦合和"热膨胀"这类通用耦合。剩下的冻胀应变、渗透系数折减、潜热源项，都需要你以"变量"或"PDE弱形式"的方式补进去。这也是做这类模拟真正的技术门槛。

2.3 缺什么就补什么：未冻水含量关系的自定义实现

冻土水热耦合里有一个任何现成接口都不提供的核心关系——未冻水含量曲线。它的含义是：土体在负温下仍然保留一部分未冻水，温度越低，未冻水含量越少，但永远不会降到零。这个关系决定了潜热释放的分布、渗透系数的折减幅度，是整个水热耦合的灵魂。

常见的经验关系式是：

w_u = a × |T|^(-b)（T < 0°C时）

其中w_u为未冻水质量含水率，T为摄氏温度，a和b是与土质有关的经验系数。粉质黏土大致a=0.12、b=0.36，砂土大约a=0.05、b=0.24。这个公式在COMSOL里可以直接作为变量定义，温度场求出后实时计算未冻水含量。

未冻水含量确定之后，两个关键派生量都能算出来：

第一个是相变潜热源项。单位体积土体在温度变化dT时发生相变的水量为ρ_d×d(w_u)，释放的潜热就是L×ρ_d×d(w_u)。在传热方程里加一个等效热容增量项即可实现。为避免数值震荡，我习惯把未冻水含量对温度的导数做平滑处理，让潜热在-0.5°C到0°C之间均匀释放，而不是集中在单点。

第二个是渗透系数折减。温度降到负温后，孔隙水部分结冰堵塞渗流通道，渗透系数急剧下降。工程上常用k = k_0×(w_u/w_0)^n来近似，n通常取3左右。把这个关系写进达西定律接口的渗透系数变量，水分场就会自动表现出"冻结后基本不流动、融化后恢复"的行为。

系数型PDE接口在这个项目里主要用来实现一些高级本构，比如考虑水分场对温度场的非线性反馈时，把它作为一个附加方程域使用。基础项目用不到，但知道这条路总没坏处。

3. 材料参数才是模拟的命门——怎么取值才靠谱

我见过不少同行在COMSOL里参数填得跟猜谜似的，算完结果漂亮但完全对不上实测。冻土模拟里，参数的可靠性直接决定结果可信度。这块我把参数分成三类，分别说说它们为什么重要、怎么取值。

3.1 块石层参数：孔隙率和当量粒径决定了对流强度

块石层的两项关键参数是孔隙率和当量粒径。

孔隙率决定空气占据的体积比例，也直接影响有效导热系数和孔隙空气的流动空间。工程上块石层孔隙率一般在0.35到0.45之间，我这个项目取0.40。当量粒径决定孔隙通道的水力直径，是计算孔隙内自然对流强弱的重要参数。块石粒径越大、孔隙通道越通畅，冬季对流越强，降温效果越明显。粒径取0.2米左右比较典型。

如果不用显式对流传热，而是采用等效导热系数法，就要按温差修正有效导热系数。冬季大温差时等效导热系数可达到静态导热系数的两倍以上，这是块石层"热半导体效应"在数学上的直接体现。

块石本身的导热系数约2.5 W/(m·K)，体积热容约2.0×10^6 J/(m³·K)。这些数据在一般岩石手册里都能查到，可靠度较高。

3.2 未冻水含量曲线：冻土水热耦合的灵魂参数

未冻水含量曲线对粉质黏土和砂土的差异非常明显。粉质黏土比表面积大，结合水含量高，负温下仍有较多未冻水；砂土则相反，温度一降到零下，绝大部分孔隙水迅速冻结。

这条曲线不仅用于潜热计算，还直接控制负温区的渗透折减和冻胀应变激活量。我做项目时没有试验数据，就参照同类土的文献取值，并做了-30%到+30%的参数敏感性分析，确认最终变形结果对a和b的取值不是特别敏感，才敢把结果作为工程参考。如果条件允许，强烈建议做一组冻土未冻水含量试验，哪怕只有两个温度点（比如-1°C和-5°C），都能把参数标定精度提高一大截。

3.3 冻土力学参数的温度相关性

冻土的弹性模量随温度降低显著增大，这是因为冰胶结作用增强了土颗粒间的连接。常温压实填土弹性模量可能在20到50 MPa左右，到-5°C时可能升到100 MPa以上。如果整个模型用固定模量，变形的应力分布就不准。

冻胀率是另一个关键参数，它表示土体在冻结过程中的体积膨胀程度。原位冻结条件下，体应变近似为可冻水体积乘以9%（水冰相变体积膨胀率）。这里要注意，不是所有含水率都参与冻胀，负温下保留下来的未冻水就不产生冻胀。

内聚力和内摩擦角在负温下也会升高，如果采用摩尔-库仑本构，建议把强度参数写成温度的分段线性函数，实测数据没有的话就参考工程冻土试验手册的典型值。

3.4 一组可以直接参考的示例参数

我基于当时项目经验整理了一组典型参数，可以用于模型初步标定。不同工程土质差异很大，直接用之前务必对照你的勘察报告校正。

注意：上表中的渗透系数是常温饱和状态值。负温区的实际渗透系数需要乘上未冻水比例的三次方左右，这一折减必须写进达西定律的渗透系数表达式，否则水分场的冻结特征完全失真。

这里放一个我当时的实际经历：第一次算出来的冻胀集中在路基中心正下方，跟监理给的实测点位完全对不上。后来排查发现是地基土渗透系数填得太保守，冻结锋面附近水分来不及向冷端迁移，凝冰量不足，冻胀当然集中在原位冻结区域。把渗透系数按实测结果调大一档、并加入未冻水折减后，冻胀峰值位置才跟实测吻合。参数之间的联动关系，比单点参数的绝对值更重要。

4. 冻胀变形怎么算：从"水冻成冰"到"位移场"

这是整个模拟里最核心、也最容易绕晕的一步。做了水热耦合，得到了温度场和水分场的时空分布，接下来怎么把"冻胀"翻译成力学计算里的变形？

4.1 三种冻胀实现方式，怎么选

我梳理过冻胀在数值模型中的三种主流实现方式：

第一种是等效线膨胀系数法。把冻胀应变表示成温度的函数，温度每下降1°C对应一个等效线膨胀系数。优点是简单直接，不需要复杂的本构推导，适合工程初步评估。

第二种是应变驱动法。根据水分场的计算结果，确定每个单元内参与冻结的水量，乘以9%的体积膨胀率得到体应变，再按各向同性分配到三个方向。这个方法物理意义明确，能反映水分迁移引起的局部冻胀差异，我最终采用的就是这个方法。

第三种是全耦合法，在土体本构中引入冰晶生长与应力耦合的弹塑性模型。学术研究常用，但参数标定难度极大，工程阶段不建议作为首选。

4.2 应变驱动法的具体操作

实际在COMSOL里，我的实现路径是这样的：

第一步，在"全局定义"里用一个污染片函数定义单位体积土体在温度T下可冻水的比例。它等于总含水率减去当前温度对应的未冻水含量，再除以总含水率。这个比例直接决定了该单元"有多少水会变成冰"。

第二步，定义冻胀体应变变量ε_vol = 0.09 × θ_freezable，其中θ_freezable是单元内参与冻结的水的体积含量。假设初始体积含水率20%，温度降到-3°C时未冻水含量约5%，那么参与冻结的水约15%，体应变就是0.09×0.15=0.0135，对应线应变约0.0045，也就是0.45%。

第三步，在"固体力学"接口中，把冻胀体应变以初始应变的方式赋给每个单元。初始应变的三个对角线分量取ε_vol的三分之一，切应变为零。这个操作相当于告诉求解器：这些单元在冻结过程中产生了各向同性的膨胀。

要注意的是，冻胀应变只有在负温区才会激活。所以我给应变变量乘了一个温度开关函数，温度低于0°C时激活，高于0°C时应变为零。开关函数不能从0跳到1，要在-0.2°C到0°C之间做平滑过渡，否则非线性求解器很容易因为突变而发散。

4.3 冻结深度与冻胀应变的空间分布

算完你会发现，冻胀并非均匀分布。它有几个明显的空间特征：

• 靠近路基表面和边坡的区域冻结深度大，冻胀量大，因为冷空气直接接触、降温最快。
• 路基中心下方因为有块石层的"热半导体效应"，基底温度被压低，冻结锋面往深推进，但这个区域的初始含水率可能低于边坡区，冻胀反而不一定最大。
• 边坡脚和路基两侧排水条件不同，含水率差异大，导致不均匀冻胀。这种差异最终反映在路面上就是纵向裂缝和横向错台。

所以结果的判读，重点不是看某一个点的位移数值，而要看冻胀量沿路基横向的分布梯度。梯度大的位置，就是路面开裂风险最高的位置。这个思路直接从模拟结果中能找到，也是你说的"把模拟变成工程判断"的关键一环。

5. 边界条件、求解器与收敛性控制——跑不动是常态

设置完物理场和材料，第一次点"计算"，大概率会收到一个"求解器未收敛"的报错。这一步几乎所有人都会经历，区别只在于谁能更快地找出问题并调整好。我在这个项目上花了差不多一周时间解决收敛问题，总结出几套行之有效的对策。

5.1 三类边界条件怎么施加才符合实际

温度边界：路基顶面是最重要的热交换面。它同时受对流换热、太阳辐射、地面长波辐射和蒸发散热的影响。简化处理时，用等效对流换热系数加外部气温函数来描述。气温函数建议直接用实测日平均气温序列，如果只有月平均，就用正弦波拟合并加密采样。路基底部取恒温边界，数值取当地年均地温，比如-1°C到-3°C。左右两侧边界在足够远的距离上按绝热或恒温处理，关键是边界离路基足够远，不会干扰核心区。

水分边界：路基顶面按降水入渗和蒸发边界处理，但冻土路基的蒸发量在冬季几乎为零，可以用分段函数控制。底部若是地下水位较深，按零通量边界处理；若较浅，则给定定水头。侧向边界按对称或零通量处理，保证水分迁移主要是竖向的。

力学边界：底部固定约束，左右两侧约束水平位移，路基表面自由变形。要注意的是，路面结构层和路基填土之间如果设置绑定约束，会限制两者的相对滑移，实际工程中这种接触面往往有薄弱层，可以做一次"接触力学"的敏感性检验，看冻胀位移差多大。

5.2 瞬态时间步与初始场的处理

这类模拟通常要算全年或多年的温度、水分、变形演化，属于瞬态分析。时间步的设置直接影响计算量和收敛性。我当时的做法是：

• 以一年为周期，先做"预处理"计算，让温度场和水分场达到一个与当地气候平衡的初始状态。具体做法是连续计算两到三个气象年的循环，取最后一个周期的结果作为正式计算的初值。
• 时间步采用自适应，并限制最大步长不超过1天。冻结期和融化期（大概每年10月到次年5月）是潜热释放最剧烈的时段，最大步长收紧到0.5天。夏季稳定期可以放宽到2到3天。
• 绝对容差的设置不宜过严，我取默认值偏松一档，因为多场强耦合问题对容差极其敏感，太严会导致求解器反复迭代甚至不收敛。

5.3 收敛失败时的排查与解决顺序

COMSOL报"未收敛"的原因很多，按概率从高到低排查：

第一，检查是否有变量突变。未冻水含量对温度的导数、渗透系数折减、冻胀应变开关，只要有阶跃跳变，非线性迭代就很容易折戟。解决办法就是做平滑处理，把突变区间扩展到0.3到0.5°C范围。

第二，检查初始值是否合理。瞬态计算之前建议先跑一个稳态温度场，把结果作为瞬态初始场，能大幅减少第一个时间步的发散概率。

第三，检查网格质量。块石层和路面交界处建议细化网格，单元质量要高于0.6。其实冻土问题用映射网格分区划分，比自由三角形网格稳定很多。

第四，检查时间步。如果连续好几次都卡在同一时间点，把该时间点的最大步长缩小一个数量级，或者改用BDF（向后差分）求解器的更高阶格式。

第五，实在不行就分步解耦调试：先只算温度场，再开水分场，最后加力学场，逐步排查耦合哪里出了问题。我调试过最头疼的一次问题，最后发现是某条耦合变量里写错了变量名，导致传热源项的物理单位差了一个量级。

6. 结果怎么读：从云图到工程结论

算完后工作其实只完成了一半，另一半是结果的解读和验证。云图不会直接告诉你"这儿要开裂"，你需要自己从场分布里提炼工程结论。

6.1 温度场：看块石层是否真的在"降温"

先看温度场冬夏两季的对比。一个有热半导体效应的块石层，冬天块石层内部温度会明显低于同深度普通路基土体，夏季则差距不大。如果在温度云图上看不到这种"单向导热"的痕迹，说明你的等效导热系数或对流传热设置没有生效。

我习惯在后处理里画一条中轴线的温度随深度变化曲线，比较1月、4月、7月、10月四个时刻。块石层的作用会以"冬季曲线明显下凹"的形式展示出来。冻结深度则提取0°C等温线的最大深度，与设计值或实测值对比，偏差一般应控制在10%到15%以内。

6.2 水分场与变形场的对应关系

水分场要重点关注冻结锋面附近的含水率增量。如果水分迁移设置正确，你会看到冻结锋面靠冷端一侧含水率显著升高，形成一条"富冰带"。这个富冰带的位置与冻胀最大的区域高度对应，两者是可以互相印证的。

变形场则要提取竖向位移云图，标注最大冻胀量及其位置。我当时的计算结果最大冻胀量约80毫米，出现在路基边缘下方约1米深度处，与实测点位一致。同时把横向不均匀冻胀量（横向位移差）单独画出来，它直接对应路面开裂风险。

6.3 工程判据与实测验证

工程上对路基变形有明确的控制指标，比如冻胀量超过20毫米/年就会明显影响路面平整度，超过50毫米就必须采取处置措施。你可以把模拟得到的冻胀量与此类指标对比，判断是否需要加强块石层厚度或添加保温层。

实测验证这块，有条件的话在路基内埋设温度探头和冻胀计，一年后对比模拟与实测的月变化曲线。我们当时把温度场计算结果与实测做对比，全年最大偏差约1.2°C，主要集中在融化期——这是因为融化过程涉及非饱和水的复杂相变和水分重分布，模型作了简化处理。这个偏差范围在工程上是可接受的，但你要在报告里明确说明简化的影响。

我个人的体会是，这套"水热力耦合+冻胀变形"的模拟流程，最大的价值不在于预测绝对数值，而在于方案对比。同一套模型，改块石层厚度、粒径、含水率条件，能很快看出哪个方案冻胀控制效果更好。这种"相对比较"的结论，在方案设计阶段比绝对数值更有说服力。

最后分享一个小经验：做冻土模拟，你一定会在某些参数的取值上反复斟酌。我的建议是别追求每一个参数都"绝对精确"，而是把重点放在两件事上——第一，所有关键参数必须有依据，哪怕是文献值或经验公式，不能拍脑袋；第二，对核心参数做敏感性分析，用图表说明结果对参数的依赖程度。这样审模型的人会问的问题，你都能给出有依据的回答。这个思路在工程报告里极其加分，也是这个项目最后能被认可的关键。