★ FEATURED ARTICLE 一提到寒区路基干这行的人第一反应通常都是那六个字——冻胀、融沉、翻浆。我前两年接了一个某寒区道路工程的数值评估任务核心就是把一段块石路基的水热力耦合行为在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实现方式关键设置温度→水分达西定律中渗透系数和孔隙度定义为温度函数或直接给未冻水含量-温度关系负温区渗透系数按未冻水比例折减水分→温度传热方程中追加相变潜热源项用温度导数形式表示潜热在-0.5°C到0°C区间平滑释放温度/水分→力学固体力学中定义温度应变和冻胀应变用插值函数驱动应变激活范围与负温区间对应这里有个细节COMSOL的多物理场节点里并没有一个叫冻土水热力耦合的一键选项它给你的是固体传热达西定律的双向耦合和热膨胀这类通用耦合。剩下的冻胀应变、渗透系数折减、潜热源项都需要你以变量或PDE弱形式的方式补进去。这也是做这类模拟真正的技术门槛。2.3 缺什么就补什么未冻水含量关系的自定义实现冻土水热耦合里有一个任何现成接口都不提供的核心关系——未冻水含量曲线。它的含义是土体在负温下仍然保留一部分未冻水温度越低未冻水含量越少但永远不会降到零。这个关系决定了潜热释放的分布、渗透系数的折减幅度是整个水热耦合的灵魂。常见的经验关系式是w_u a × |T|^(-b)T 0°C时其中w_u为未冻水质量含水率T为摄氏温度a和b是与土质有关的经验系数。粉质黏土大致a0.12、b0.36砂土大约a0.05、b0.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 一组可以直接参考的示例参数我基于当时项目经验整理了一组典型参数可以用于模型初步标定。不同工程土质差异很大直接用之前务必对照你的勘察报告校正。参数块石层路基填土地基土孔隙率0.400.300.35初始体积含水率0.050.200.25导热系数 W/(m·K)2.51.51.4体积热容 J/(m³·K)2.0e62.4e62.6e6渗透系数 m/s1.0e-41.0e-61.0e-7弹性模量 MPa804030泊松比0.250.300.35注意上表中的渗透系数是常温饱和状态值。负温区的实际渗透系数需要乘上未冻水比例的三次方左右这一折减必须写进达西定律的渗透系数表达式否则水分场的冻结特征完全失真。这里放一个我当时的实际经历第一次算出来的冻胀集中在路基中心正下方跟监理给的实测点位完全对不上。后来排查发现是地基土渗透系数填得太保守冻结锋面附近水分来不及向冷端迁移凝冰量不足冻胀当然集中在原位冻结区域。把渗透系数按实测结果调大一档、并加入未冻水折减后冻胀峰值位置才跟实测吻合。参数之间的联动关系比单点参数的绝对值更重要。4. 冻胀变形怎么算从水冻成冰到位移场这是整个模拟里最核心、也最容易绕晕的一步。做了水热耦合得到了温度场和水分场的时空分布接下来怎么把冻胀翻译成力学计算里的变形4.1 三种冻胀实现方式怎么选我梳理过冻胀在数值模型中的三种主流实现方式第一种是等效线膨胀系数法。把冻胀应变表示成温度的函数温度每下降1°C对应一个等效线膨胀系数。优点是简单直接不需要复杂的本构推导适合工程初步评估。第二种是应变驱动法。根据水分场的计算结果确定每个单元内参与冻结的水量乘以9%的体积膨胀率得到体应变再按各向同性分配到三个方向。这个方法物理意义明确能反映水分迁移引起的局部冻胀差异我最终采用的就是这个方法。第三种是全耦合法在土体本构中引入冰晶生长与应力耦合的弹塑性模型。学术研究常用但参数标定难度极大工程阶段不建议作为首选。4.2 应变驱动法的具体操作实际在COMSOL里我的实现路径是这样的第一步在全局定义里用一个污染片函数定义单位体积土体在温度T下可冻水的比例。它等于总含水率减去当前温度对应的未冻水含量再除以总含水率。这个比例直接决定了该单元有多少水会变成冰。第二步定义冻胀体应变变量ε_vol 0.09 × θ_freezable其中θ_freezable是单元内参与冻结的水的体积含量。假设初始体积含水率20%温度降到-3°C时未冻水含量约5%那么参与冻结的水约15%体应变就是0.09×0.150.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主要集中在融化期——这是因为融化过程涉及非饱和水的复杂相变和水分重分布模型作了简化处理。这个偏差范围在工程上是可接受的但你要在报告里明确说明简化的影响。我个人的体会是这套水热力耦合冻胀变形的模拟流程最大的价值不在于预测绝对数值而在于方案对比。同一套模型改块石层厚度、粒径、含水率条件能很快看出哪个方案冻胀控制效果更好。这种相对比较的结论在方案设计阶段比绝对数值更有说服力。最后分享一个小经验做冻土模拟你一定会在某些参数的取值上反复斟酌。我的建议是别追求每一个参数都绝对精确而是把重点放在两件事上——第一所有关键参数必须有依据哪怕是文献值或经验公式不能拍脑袋第二对核心参数做敏感性分析用图表说明结果对参数的依赖程度。这样审模型的人会问的问题你都能给出有依据的回答。这个思路在工程报告里极其加分也是这个项目最后能被认可的关键。
阅读完成 · 觉得有帮助？
