说实话看到“水热力三场耦合”这个标题很多老玩家第一反应是“又是那种算到半夜突然不收敛的模型”。但这类问题在岩土、地热、冻土工程里实在太常碰到了躲不开。比如地埋管周围土体受热后的水分迁移、核废料处置库中的缓冲层膨胀、季节性冻融路基的变形预测核心都是水、热、力三个物理场在同一个多孔介质里互相咬合。二维轴对称模型是这类问题最划算的解法。几何上是旋转体比如一根竖直的加热棒插在土里或者一个圆形热源埋在地层中沿着中心轴切开一个截面把截面当成二维区域来算得到的位移、温度、孔压结果在360度方向上完全一样却能把三维网格规模砍掉一个数量级。COMSOL里选“二维轴对称”空间维度本质就是在柱坐标系里求解物理上等于默认了环向没有变化把三维问题降成二维但结果又保留圆柱真实的体积效应。这个设定对圆柱形试件、钻孔热源、轴对称井筒这些工况非常合适。这篇文章不打算给你念手册直接说清楚怎么搭模型、哪些参数是绕不开的坑、求解时怎么避免跟收敛性搏斗到凌晨。案例取材是个很常见的场景一根垂直热源插入饱和-非饱和土柱热源持续放热土体同时发生水分蒸发迁移、孔压变化和热膨胀变形。我会把建模步骤、耦合方程、边界条件、网格策略、求解器配置一条条拆开适合刚做完单场分析、准备迈入多场耦合的人也适合被THM模型折磨过的老手回来查漏补缺。1. 案例背景与建模思路拆解1.1 三场耦合到底在耦合什么很多初学者把“水热力耦合”想成三个方程放在一起算这是误区。耦合的本质是“一个场的状态变量作为另一个场的输入参数”。水分场里含水率的变化会改变土体导热系数和热容这是水对热的作用温度升高导致水的黏度和密度变化进而改变渗透率和渗流速度这是热对水的作用孔压和温度的变化又会引起有效应力改变、土体膨胀或收缩这是力和前两场的相互作用反过来孔隙比的变化会改变渗透率和热物性参数形成闭环。以土柱热源问题为例典型物理过程是这样的热源开启后周围土体温度梯度形成靠近热源的孔隙水受热密度减小、黏度下降水分开始向低温区迁移同时部分水蒸发含水率重新分布含水率一变土体导热系数跟着变温度场又受影响孔隙水压力在热源附近升高有效应力下降土体产生膨胀变形膨胀变形再改变孔隙率。整个过程互相嵌套时间尺度从几小时到几年必须用瞬态计算才能看到真实演化。这个案例里我建议用三个标准物理场接口来完成地下水流模块的“理查兹方程”接口算非饱和渗流传热模块的“固体传热”接口算温度场结构力学模块的“固体力学”接口算变形。三者通过COMSOL的多物理场耦合节点手动连接。版本在6.0以上的朋友也可以直接找预定义的THM耦合但底层数学本质跟我们手动搭是一样的了解原理之后调起多物理场节点会更有底。1.2 为什么用二维轴对称而不是完整三维建模第一步就该想清楚空间维度。如果热源是竖直圆柱周围土体性质在环向没有差异边界条件也是轴对称的那完整三维模型有接近一半的自由度是纯浪费。二维轴对称模型在COMSOL里坐标是(r, z)两个分量求解域是一个截面但物理上它代表一个完整的圆柱体。所有方程中的梯度散度算子都会自带圆柱坐标的几何项比如径向的拉普拉斯算子里有个(1/r)的项这不是人为加的是柱坐标系下本来就有的COMSOL会自动处理。对比一下工作量三维模型哪怕用较粗网格30万到50万自由度很常见二维轴对称同等精度只需要2万到5万自由度。瞬态多物理场耦合本来就是非线性迭代加时间推进自由度少一个数量级单步求解时间和内存占用都大幅下降而且网格可以放心加密到热源附近捕捉到更陡的温度梯度和水压力梯度。所以我给所有做这类问题的人第一个建议就是先确认几何和荷载能不能轴对称能就绝对别碰三维。网格在轴对称模型里还有一层特殊意义r方向靠近对称轴时单元面积趋向于零网格质量天然变差。后面网格部分我会讲怎么处理这里先记住这个隐患就行。2. COMSOL建模与物理场配置2.1 从零搭建模块选择与全局定义打开COMSOL新建模型时选择空间维度“二维轴对称”然后依次添加物理场。拿6.4版本界面为例“添加物理场”搜索框里输入Richards会出现“理查兹方程dl”接口传热选择“固体传热ht”力学选择“固体力学solid”。这里有个细节要注意理查兹方程的因变量是压力水头p单位是米不是帕很多人第一次用在这里卡住。固体传热因变量是温度T固体力学因变量是位移u和v这里v代表z向位移别搞混。添加完物理场后第一件事是去“全局定义”里把材料参数写成参数和变量。不要直接在材料节点里硬填数字因为后面这些参数会随着含水率、温度变化必须用变量表达式。比如渗透率不能是个常数而是要写成饱和渗透率乘以相对渗透率系数导热系数也要写成含水率的函数。这样做的另一个好处是后处理时要画参数云图、做参数化扫描都很方便。我常用的做法是把所有物理参数先放进“参数”节点把随场变化的中间量写进“变量”节点计算时COMSOL每个求解步都会自动重新求值。单位问题再说一句几何尺寸用米时间用秒压力水头用米温度用开尔文。压力水头乘以水的重度ρg才能换算成帕斯卡而固体力学中孔隙水压力的单位必须用帕斯卡所以耦合到有效应力公式时要乘上ρg这个换算漏了的话应力场结果会差好几个数量级典型的低级错误但坑过无数人。2.2 几何简化、边界条件与初始条件的工程处理几何我用一个矩形来表示土柱截面竖直方向深度20米径向半径10米。热源不是实体模型的一部分而是简化成一个边界条件或者一个很小的几何区域。这里有个建模习惯问题如果热源本身刚度温度很高把整个热源几何画出来会让网格在热源内部加密大大增加计算量。工程中如果只关心土体响应热源可以按面热源或体热源载荷施加避免不必要的网格开销。我把热源处理为轴对称轴上r0到0.1米的一段小区域在域上施加热源功率密度模拟一根半径0.1米的加热棒。边界条件分三场分别设置。轴对称轴上r0处三个场都有对称条件力学上r方向位移为零热和水的法向通量为零。土柱外边界r10米处力学上可以设为辊支撑即法向位移被约束但切向可以自由滑动模拟周围土体对柱体的侧向约束温度边界设为初始地层温度如果是无限大区域也可以用热绝缘后再加大半径来近似但直接给定温度更稳定渗流边界设为常水头模拟远处水分补给。土柱底部z0处全部固定绝热不透水。顶部z20米处是自由面力学上自由温度设为地表温度渗流上如果模拟地表蒸发可以给一个通量边界如果模拟封闭系统就设为零通量。初始条件也不能乱拍。温度场初值给一个均匀的地层温度T0水场初值根据地下水位位置给一个静水压力水头分布位移场初值设为零。一个常被忽略的坑是初值不一致会导致求解器在第一步疯狂迭代甚至直接发散。比如初始孔压如果和初始有效应力对不上力学场在t0就会产生一个虚假的瞬态响应。稳妥的办法是先跑一个稳态求解把稳态结果当瞬态的初始值后面求解设置部分我会细讲。2.3 用多物理场耦合节点把三场绑在一起COMSOL里多物理场节点的作用是在各物理场方程里自动添加耦合项。这里推荐用两个预定义耦合再加一个手动变量耦合。第一个是“固体传热-理查兹方程”之间的非等温流动耦合。这个耦合会让水的黏度、密度随温度变化同时把水的流速带入传热方程的对流项。第二个是“固体力学-理查兹方程”的孔隙弹性耦合核心是有效应力原理总应力等于有效应力加上孔隙水压力贡献BioT系数默认取1时对应土颗粒不可压缩的假设实际黏土中取0.6到1之间砂土接近1。第三个是热膨胀耦合在固体力学中把温度变化产生的热应变加上。这三个耦合节点设置好三场的相互作用链就闭环了。如果用的是6.4版本多物理场节点里可以直接搜到这些预定义耦合。如果是老版本就要自己在各物理场的域方程里手动修改在传热方程里加对流项在力学方程里加孔压项和热应变项。手动改域方程对新手不太友好但理解原理后会发现也没什么神秘的本质上就是把耦合项当成源项或材料参数写进去。3. 耦合参数详解与核心变量实现3.1 土水特征曲线与相对渗透率水分场的灵魂非饱和渗流和饱和渗流最大的区别就是渗透率不再是一个常数。土中含水率降低时大孔隙先排水水流只能走小孔隙和薄膜水渗透能力指数级下降。这个关系通常分成两部分描述土水特征曲线SWCC把基质吸力和含水率对应起来相对渗透率曲线kr(Se)把有效饱和度对应到渗透率折减系数。工程中最常用的是van Genuchten模型表达式为Se (θ-θr)/(θs-θr) [1(α|h|)^n]^(-m)其中m 1 - 1/n相对渗透率kr Se^0.5 * [1 - (1 - Se^(1/m))^m]^2参数α和n不是随遍取的拟合值如果手头没有试验数据可以参考典型值砂土α约1.5到51/mn约2到4黏土α约0.05到0.51/mn约1.1到2。这个模型比较灵活一段代码就能覆盖大多数土壤类型。在COMSOL里我把θr、θs、α、n、m都放进参数表在“变量”节点里定义Se和kr然后渗透率表达式写成Ks*kr(Se)理查兹方程接口的材料设置里直接引用这个变量。这里有个大部分人都会踩的坑理查兹方程可以对含水率做限制防止数值振荡导致Se超过1或小于0。但直接在变量里写Seclamp(...)会让表达式不可导影响牛顿迭代的收敛性。更好的做法是用连续可导的平滑函数近似或者把极限值设成θr和θs附近的微小偏移避免求解器在饱和线附近反复震荡。3.2 热物性参数的温度效应热场不是单纯的热传导土的热物性参数会随着含水率变化这一点经常被忽略。干砂的导热系数大概0.3 W/(m·K)饱和砂可以到2.5左右差了快一个数量级。所以温度场和水分场之间不仅是对流换热材料参数本身就在强烈耦合。工程上常用加权法估算λ λdry (Se)^0.5 * (λsat - λdry)体积热容也类似C C_dry Se * (C_sat - C_dry)。这样处理简单有效实测粗糙度可以接受。温度影响水的物性参数也很关键。水的黏度从20度到90度大概下降三分之二渗透率会跟着变大水密度在4度前后有个反常区但这在工程温度范围内不太敏感主要考虑温度对黏度的影响。COMSOL的“非等温流动”耦合节点会自动把水的密度和黏度设为温度函数如果纯手动建模需要自己写黏度公式比如用Vogel方程拟合温度区间。对流项是热场里另一个容易出问题的地方。热水注入导致的水分迁移通道里的热量输运往往比单纯热传导还重要。传热接口默认只算热传导如果没有在物理场里启用对流贡献温度场会明显偏慢。因此在固体传热接口的多孔介质设置里要确认把达西速度或理查兹方程的流速作为输运速度引入对流项。这一步不设置的话你的计算结果可能比实际慢了一倍还不自知。3.3 有效应力原理与热应变力场如何接住前两场力学场在本案例中的作用是计算土柱的膨胀、压缩和应力变化。核心方程是有效应力原理σ_total σ_eff αBiot * p * δ_ij其中p是孔隙水压力单位是帕斯卡。前面说过理查兹方程的因变量是压力水头h单位米所以要写成p ρw * g * h。αBiot在饱和土中通常取1但很多土力学家喜欢用0.8到1之间的值因为它实际上和土颗粒压缩性有关。如果你暂时没有实测值取0.9是个合理的中间值。热应变项为ε_th αT * (T - T_ref)其中αT是热膨胀系数T_ref是参考温度。土的线膨胀系数一般在10^-6到10^-5量级比金属低一两个数量级。但三场耦合问题的有趣之处在于热膨胀本身产生的变形虽然不大但它会改变孔隙比孔隙比一变渗透率和热物性参数又变了最终影响水分和温度分布。这就是所谓的全耦合不是简单的单向驱动。孔隙率与渗透率之间的关系工程上常用Kozeny-Carman公式修正k k0 * (e/e0)^3 * (1e0)/(1e)。这个公式非常好用把力学场算出来的孔隙比变化直接反馈到渗流场。在COMSOL里实现也不复杂定义一个变量e e0 (εvol - 0)再用Kozeny-Carman公式更新渗透率参数。我们在后处理时最喜欢画这一项因为一眼就能看出热-水-力耦合的“反馈”到底有没有起作用。4. 网格、求解器设置与加速技巧4.1 网格策略给温度梯度和水压梯度足够的分辨率网格划分直接决定多物理场耦合计算的成败。网格太粗温度前锋和水压前锋数值抹平结果毫无物理意义网格太细每个时间步的非线性迭代都要算半天瞬态推进上百步根本没法跑。平衡点在哪儿我一般先用一个粗网格试算大致看看温度梯度和水压力梯度集中在哪个区域然后针对性地加密。在这个案例里热源周围0到2米范围内温度梯度最陡水分迁移最活跃网格尺寸安排0.05米就够用中间区域可以放到0.2米外边界附近温度几乎没变化0.5米也完全可以。径向和轴向都按这个逻辑做分布。COMSOL里可以用“尺寸”节点配合“分布”子节点设置单元数量的比例也可以用“边界层”节点在热源边界附近添加边界层网格用来捕捉热边界层内的陡梯度。还有一个很实用的技巧设置完网格后先跑一个短时间的瞬态比如100秒停下来画一下温度场看看热源附近温度梯度是不是平滑如果出现锯齿状说明网格够不上梯度立刻加密。网格质量是轴对称模型里必须检查的项。r0附近因为单元体积趋近零质量最容易出现红色报警特别是三角形单元拉得很长时。解决方案是把对称轴附近的网格做成结构化四边形网格或者至少保证对称轴周边单元的长宽比不要超过5。网格划分完成后点“统计”按钮把最小单元质量目标定在0.1以上低于这个值建议重新调几何或网格不然算到一半遇到网格雅可比奇异直接报错。4.2 瞬态求解器与时间步控制水热力耦合本质上是一个强非线性瞬态问题。求解器建议用“瞬态”研究时间步格式选BDF向后差分公式阶数2。BDF对刚性方程比较友好不像显式格式那样受限于CFL条件。时间步不要用固定步长而是选“自由”模式让求解器根据局部截断误差自动调整步长。但自动时间步有个问题在物理场剧烈变化的起始阶段求解器可能会试图用非常大的步长跳跃导致非线性迭代失败。所以要设置初始步长我一般设成总计算时长的千分之一比如总共模拟1e7秒初始步长就设1e4秒之后让求解器自己加速。非线性迭代设置里有两个关键的容差。绝对容差和相对容差默认值都比较松多物理场计算建议收紧相对容差1e-3绝对容差1e-5。更重要的是对每个物理场的容差单独控制因为压力水头的量级和温度不同位移更是小几个数量级。在“因变量值”的容差设置里可以分别设置水头、温度、位移的绝对容差比如水头0.01米、温度0.1开尔文、位移1e-6米。这一步很多人不做结果就是求解器觉得“已经收敛了”但实际上位移场和孔压场还在明显抖动。阻尼因子是另一个值得关注的旋钮。BDF格式中求解器默认启用阻尼因子最小值0.1。如果遇到剧烈非线性可以把最小阻尼因子降到0.01甚至0.001代价是收敛速度变慢但稳定性会大幅提升。我见过太多人一碰到不收敛就把网格加密其实先调阻尼因子往往更有效而且不用重新剖分网格迭代几次就试出来了。4.3 “先解耦后耦合”的分步求解法这是整个案例里最值得先做的一步。众所周知最稳的三场耦合策略是“分阶段耦合”而不是一上来就全耦合求解。我的做法分三步第一步关闭传热和力学接口只保留理查兹方程先算水分场在常温下的稳态渗流分布。这样能得到一个满足边界条件和初始条件的水压场。第二步打开传热接口保持力学接口关闭以第一步算出的水压场作为初始值算一个热-水双场耦合的瞬态过程。这个过程里可以用一个简化的恒定孔压假设让热场先演化起来。第三步打开力学接口保留前两步的所有结果作为初始值才算完整的三场耦合瞬态。这样每一步只有一个或两个新物理场参与迭代非线性耦合的强度逐级提升收敛难度成倍下降而且每一阶段都能对照单场或双场结果做验证出错了也知道错在哪一段。如果计算资源紧还有一个更近一步的加速技巧先用粗网格跑一个完整耦合过程得到大致的响应包络和关键时间点然后用细网格只计算影响最大的前三步时间范围。对于参数敏感性分析这一步几乎可以省一半时间。5. 常见问题与排查实录5.1 不收敛、负压力和振荡三个高频故障的处置顺序做多物理场耦合遇到不收敛基本是常态。关键是怎么系统地排查而不是抓瞎。第一步看求解器日志里的“残差”。残差一直不下降解决方案中很大可能是初始条件不合适。典型案例是力学场初始孔压为0但渗流场算出边界上孔压非零这一步应该在初始化时就报错。解决方式是执行“研究”前的“初始化”按钮或者用稳态预计算。第二步看变量是否超出物性函数的定义域。经常出现的情况是Se算出来小于0导致相对渗透率kr表达式里出现负数开方或负数幂直接NaN。这时需要用平滑截断函数包裹Se而不是简单用max/min因为max/min不可导。COMSOL里可以定义Se_smooth 0.5*((Se^2eps)^0.5 Se)这种方式既平滑又基本等于原始值。第三步检查时间步长是否太激进。如果是BDF格式求解器在收敛失败时会自动减半时间步但减步长不是无限度的减到最小步长还失败就会停止。这时候不要暴力增加最大时间步数而是先调大“非线性求解器”里的最大迭代次数从默认25加到50同时降低阻尼因子下限。这套组合拳能解决八成收敛问题。最后一个常被忽略的原因模型本身病态比如边界条件冲突或几何有尖角。检查一下几何角落处是否同时施加了多点约束和集中载荷往往会在这些局部点造成应力奇异。这类问题无法通过调求解器解决只能修改边界条件或几何。5.2 网格畸变与大变形movement mesh该何时登场三场耦合里力学场如果产生大位移比如膨胀性土吸水后体积膨胀百分之十几固体力学接口的拉格朗日网格就会跟着变形网格扭曲到一定程度后雅可比矩阵恶化计算发散。此时有两个方案一是使用“移动网格”接口动网格让网格节点跟随材料变形移动同时用自动重剖分当网格质量低于阈值时重新划分。二是用小变形假设忽略几何非线性只在材料本构里引入“体应变”对渗透率和孔隙率的影响这是我在大多数工程案例中的首选因为热致膨胀大多数情况下在5%以内小变形假设完全够用而且计算稳定得多。判别是否需要动网格的标准很简单看计算结果中最大位移量级相对特征尺寸的比例。如果最大位移小于特征尺寸的1%完全没有必要用动网格纯属给自己找麻烦。如果达到5%以上动网格或网格重剖分就得考虑了。我这里案例热膨胀变形很小所以用固定网格就足够。真正必须用移动网格的场景是模拟注浆、压裂这类局部应变很大的工艺过程那种情况里r0附近和裂缝尖端的网格畸变控制才真正考验功力。5.3 结果后处理与验证别让你的漂亮云图骗了你计算收敛了、云图画出来了不等于结果一定正确。我自己有个硬性习惯后处理第一件事不是看云图而是查守恒量。比如在域上对水的质量平衡做积分入口通量减去内部累积量差值应该在数值噪声范围内。COMSOL后处理里可以添加“体积积分”和“边界积分”分别算累积量和通量两者对比一下如果差了超过5%说明时间步或网格精度不足结果只能算定性不能定量。第二件事是并行验证单场解。至少要做一次“去掉耦合”的对照模型只算热传导检查热源周围温度随时间的半解析解比如一维柱坐标系下的径向传热公式。把耦合模型在第一个短时间内的温度分布与纯导热解对比如果偏差在可接受范围内说明热场部分没有打错。这个方法花不了半小时但能帮你排除90%的建模错误。最后建议后处理里不要只画默认的“表面图”多画几个“截面图”和“点探针图”。在热源边界、距热源0.5米、1米等位置放几个探针盯着这些点的温度、孔压、位移时间曲线。这些曲线最能反映耦合动态的细节比如温度爬升过程中孔压有没有先升后降的“不排水热膨胀”特征有的话说明力学-流体耦合真的在起作用。如果有实验数据或文献曲线用探针数据直接叠合上去对比比云图对比更直观。写在最后再分享一个小经验如果你准备把这个案例扩展到自己工程里我的建议是先做“参数敏感性扫描”不要一上来就做精细实测值建模。用COMSOL的扫描研究功能对渗透率、热容量、热膨胀系数这几个核心参数各跑几条曲线看看输出对哪个参数最敏感。这一步能让你知道手头的数据哪项要优先补实验哪项用工程经验估算就够了。我在实际项目里用这招省下过不少时间别人的论文里给了一堆实测参数但不是每个参数都对结果有显著影响把所有精度都堆在对结果无关紧要的参数上是奢侈也是技术上的不成熟。水热力三场耦合模型难的不是数学和软件操作而是理解每个参数在链条里扮演的角色把这个链条想通了模型自然就稳了。