项目标题里“水热力三场耦合”这几个字放到工程模拟里其实就是固体力学、流体渗流、温度传导三个物理场互相咬合的问题。我前阵子用COMSOL 6.4做深地工程场景的THM耦合模拟时就选了二维轴对称模型来跑而不是一上来就砸资源上全三维。这个取舍一方面是工程几何本身有对称性另一方面是求解器调试和网格控制都会轻松得多。这篇文章整理的是我当时从物理场接口选择、几何搭建、材料参数标定到耦合节点设置、求解器收敛性调试的完整过程附带踩坑记录希望能给同样在做地热开发、核废料处置或地下储能数值模拟的朋友省点时间。1. 三场耦合的物理本质与轴对称建模思路1.1 水热力三场耦合的物理本质与数学模型先搞清楚一个问题三场耦合到底在耦合什么。力场和水场之间的耦合是岩土工程里最经典的“有效应力原理”土体或岩体骨架承受的总应力一部分由固体颗粒承担一部分由孔隙水承担。写成通俗的式子就是σ总应力 σ有效 α_B × p孔隙水压力也就是说孔隙压力增加有效应力就减小如果有效应力下降到岩石承载极限以下就可能产生压裂、沉降或者剪切破坏。反过来岩石骨架发生压缩或膨胀也会改变孔隙体积从而迫使孔隙水压力改变这就是Skempton效应。温度场也不独立。温度升高引起固体骨架热膨胀直接改变应变场这是热-力耦合温度升高会让水的粘度下降、密度变化从而改变渗流速度这是热-水耦合而水的流动会携带热量对流渗流速度本身又直接影响温度分布这是水-热耦合。三套方程最后串成一条链力学方程管位移达西方程管压力能量方程管温度。从控制方程角度THM耦合可以粗暴地概括成三个偏微分方程第一固体骨架的力平衡。∇·σ ρ_b b 0其中σ必须用有效应力表示并且计入热应变项α_T(T−T_ref)。这里T_ref千万不能随便给给了错的参考温度整个应力场就是歪的后面我会单独讲。第二孔隙水的质量守恒。S ∂p/∂t α_B ∂ε_v/∂t − ∇·(k/μ ∇p ) 0。注意中间这个α_B ∂ε_v/∂t就是力学变形对水压力的反馈项很多人漏了它结果算出来孔隙压力变化离谱。第三能量守恒。(ρC_p)_eq ∂T/∂t ρ_f C_f u·∇T − ∇·(k_eq ∇T) 0。里面ρ_f C_f u·∇T是对流项u用的是达西速度这就是把渗流场和温度场绑在一起的钥匙。1.2 为什么用二维轴对称模型而不是直接上三维我这次算的问题是一个竖直地热井周边地层在注冷水过程中的响应井周围岩体可以近似看作以井轴为中心旋转对称的结构。注入点附近的压力锥、温度锥都是沿径向向外扩展的最大的梯度方向就是径向r和竖向z周向θ方向几乎没有梯度变化。这种条件下直接用二维轴对称坐标系等于把三维问题降了一维同时又完整保留了r方向上的径向扩散和z方向上的分层效应。有人会问三维模型也不难建为什么非省这点算力我的体会是THM三场耦合是非线性的求解起来往往要反复调。如果直接上全三维网格数量至少是轴对称模型的十几倍每次试错都要等很久而且一旦不收敛你很难判断到底是物理设置问题还是网格问题。用轴对称模型先把所有条件摸清楚再决定要不要扩展到三维是效率高很多的路径。COMSOL里实现二维轴对称的操作不复杂新建组件时选择“二维轴对称”空间维度几何平面就是(r, z)软件默认r0处为旋转对称轴。建模时只需要画一个矩形代表地层剖面然后用“物理场接口”在域内定义固体力学、达西定律、多孔介质传热边界条件在r0处自动按对称处理。1.3 物理场接口的取舍用内置多物理场还是自定义PDECOMSOL里做THM耦合有两条路一条是直接用物理场接口加多物理场耦合节点另一条是自己在“数学模块”里用系数型PDE甚至一般形式PDE去写三套方程再用耦合项链接。我建议绝大多数人走第一条路。因为COMSOL 6.4里已经有相当成熟的“达西定律接口”、“多孔介质传热接口”、“固体力学接口”再加上多物理场节点里的“多孔弹性”和“热膨胀”已经能把主要物理勾连起来。内置物理场的好处是边界条件语义清晰比如达西接口里可以直接设置“压力”边界传热接口里可以直接设“热通量”后处理结果也有现成的压力云图、温度云图省去自己写方程时边界的表达问题。只有在遇到内置接口没有覆盖的特殊本构比如考虑渗透率随应力动态变化、或者水合物相变影响孔隙比时才需要动用自定义PDE。自定义PDE灵活性高但麻烦也很大尤其是初边值条件一旦写错数值发散连报错都找不到方向。后文我会专门给一个渗透率随有效应力变化的扩展思路。2. 几何建模与材料参数准备2.1 几何建模步骤与单位统一模拟对象我取了一个半径500米、深度1000米的地层圆柱体井筒半径0.1米回注段取顶部20米。几何画起来就是个4.9999米×1000米的矩形r方向从0.1到500z方向从-1000到0。COMSOL默认单位我们统一用SI制几何单位设为米。创建几何时直接添加一个矩形宽设置成499.9高设置成1000左下角坐标(r0.1, z-1000)。然后在这个平面域上构建网格。虽然这个矩形长宽比接近500倍但因为轴对称模型本来就是细长柱体剖面所以只要网格划分得当不会影响解的正确性。这里有一点值得提醒二维轴对称模型的几何单位会直接影响表达式里的r值。写对的是一个从0.1米开始的域如果画几何时把半径起点画成0然后想用“边界条件-井筒”去模拟井筒就会出现一个半径为零的奇异点尤其是在达西定律的径向压力梯度里r在分母处会出现奇异性。所以正确做法是让几何域从井筒半径开始把井筒作为内边界处理而不是把旋转轴本身当成井筒。2.2 材料参数取值与依据材料参数是这类模拟最容易被质疑的地方一定要挨个交代清楚。我用的基本参数如下参数数值单位备注岩石弹性模量10GPa中等砂岩/泥岩量级泊松比0.251经验值岩石密度2200kg/m³干密度孔隙度0.21有效性有待实测校准渗透率1e-17m²约10 mDBiot-Willis系数0.81非压实砂岩经验值岩石热导率2.5W/(m·K)含孔隙介质等效值岩石体积热容2.2e6J/(m³·K)固体骨架孔隙水混合热膨胀系数1e-51/K岩石骨架线膨胀系数流体密度1000kg/m³常密度简化流体粘度1e-3Pa·s20℃水未建温度依赖流体比热4.2e3J/(kg·K)水的比热参数里最容易引起争议的是Biot-Willis系数α_B。有些默认材料库直接给1那意味着颗粒不可压缩骨架体积变化完全由孔隙体积贡献。实际砂岩不会这么理想0.8是个合理值。另一个要注意的是体积热容COMSOL里传热模块默认可以用多孔介质混合公式ρC_p_eq θ_p·ρ_f·C_p_f (1−θ_p)·ρ_s·C_p_s。如果你直接填了岩石干骨架的热容而忽略水温度场计算就会出现明显偏差尤其在地热这类水占比重的问题里。2.3 初始条件与边界条件的设置初始孔隙水压力要按静水压力分布而不是一个常数。假设孔隙与地表连通p0(z) p_top ρ_f·g·(−z)。取地表压力为0g9.81在地下1000米处p0数值就是9.81 MPa。很多初次做的人直接在“达西定律-初始值”里填“1e7 Pa”结果整个压力场从第一秒开始就在向静水压力平衡调整产生一个人为压力脉冲干扰了真实注水信号。温度初始场按地温梯度设T0(z) 30 0.03×(−z)在地表是30℃在地下1000米是60℃。这个梯度大约每百米3℃符合大陆地壳常见地温梯度。边界条件分成几组井筒边界r0.1达西压力固定为p_in12 MPa高于初始静水压力形成注水压差温度固定为T_in10℃模拟冷水的注入。远场边界r500达西压力按初始静水压力设置温度按初始地温设置岩石位移设置为自由或滚动支撑边界防止刚体漂移。顶部z0设为排水自由边界压力固定为0温度按空气对流边界处理。底部z-1000竖向位移约束为零底部设为绝热和零压力梯度。关键原则初始条件必须尽量落在稳态解附近。也就是说在未开井之前把模型先算一遍稳态确认位移、压力、温度场没有异常“启动响应”再在井筒边界上加上扰动。3. 三大物理场耦合实现细节3.1 固体力学接口中的有效应力与热应力在固体力学接口里默认情况下解的是总应力但本构关系定义时我们关心的是有效应力增量。需要在“线弹性材料”节点下加“热膨胀”子特征材料热膨胀系数填α_T参考温度填T_ref。然后通过“多孔弹性”多物理场耦合节点把达西接口算出来的孔隙压力p传给力学。实际操作时可以点开多物理场节点选择“添加多物理场”找到“多孔弹性耦合”它会自动把固体力学和达西定律配对并且生成一个α_B×p的贡献项。这个节点里有“Biot-Willis系数”输入框就填0.8。再加上“热膨胀”节点后力学场的驱动因素就齐了外力载荷、孔隙压力、热应变三者会对位移同时起作用。需要注意如果把孔隙压力直接加载到固体表面而不是作为体积力处理在井筒边界和表面边界上容易出现应力集中。在多孔弹性耦合里孔隙压力的作用是通过本构关系以α_Bp成出现的体应变实现的一般不需要额外在边界上加压力。如果你用边界载荷代替等于重复计入了孔隙压力误差会放大很多。3.2 达西定律接口中的渗流与孔隙压力达西定律接口的设置相对简单核心参数就是渗透率k和流体粘度μ。方程内部用Darcy速度u −(k/μ)·(∇p ρ_f·g·∇z)。这里流体密度默认是常数如果温度跨度大可以把密度改成温度的函数但那会带来浮力对流也让非线性程度上一个台阶。首次跑通案例时建议先保持常密度等基准结果稳定后再考虑浮力项。达西接口初始值里填p0域上用水文地质参数里的“储水系数S”。S的值很重要岩石的储水系数典型范围是1e-6到1e-4 1/m。填得太大压力扩散速度就慢注水压力传播路径也会失真填太小瞬态过程又容易出现刚性方程很难收敛。我习惯先取1e-6 1/m然后用半解析井底压力历史校验后微调。3.3 多孔介质传热接口中的对流与热传递温度场接口我用的是“多孔介质传热”它会自动根据孔隙度混合固体和流体的热物性。这个接口里有个“达西速度”子特征可以把达西定律计算出来的速度场作为对流速度导入。设置路径是传热接口–热通量特征下把速度场选择为“来自达西定律接口”。这里最关键的物理效应是热对流项。井筒注入冷水后冷锋会以速度u推进。如果忽略对流项则冷锋只会靠热传导缓慢扩散计算结果会跟实际完全脱节。你可以简单估算一下热扩散系数约1e-6 m²/s传导扩散到100米要几十年而达西速度若有1e-6 m/s量级对流推进100米只需要几年时间尺度差一个量级还多。所以“多孔介质传热”里的达西速度必须正确关联。3.4 多物理场节点绑定方式与顺序在COMSOL 6.4里我习惯按下面的顺序建立多物理场链第一个节点多孔弹性耦合固体力学↔达西定律保证孔隙压力和骨架变形的双向耦合。第二个节点热膨胀固体力学↔多孔介质传热让温度场驱动力学应变。第三个节点达西定律 ↔ 多孔介质传热通常传热接口里已经引用了达西速度但还需要在“多孔介质传热”节点里确认勾选了“流动传热”或者“多物理场耦合-多孔介质传热/达西定律”选项。耦合顺序背后有个逻辑先算渗流场得到速度再算温度场得到热对流和热源再把温度和压力带入力学计算位移反过来位移影响储水系数和孔隙度。求解器按这个顺序迭代比三场同时全耦合要稳得多。4. 网格、求解器与收敛性调试4.1 网格划分的实战心得对于轴对称细长模型最忌讳的是用默认自由三角形网格因为长宽比太大剖出来密密麻麻而且单元质量差。正确做法是先选择“映射”网格把整个矩形剖成四边形再在r0.1附近和顶部注水段加密。网格密度怎么判断以“井筒附近的压力梯度”和你关心的“冷锋推进半径”为标准。注水段20米长度内沿z方向网格加密到0.5米左右沿r方向从井筒向外按1.05倍比例增长最外层单元可以放到30到50米。边界建议加两层到三层的边界层单元用来捕捉井筒边上的温度和压力边界层。我最初直接把全区域用均匀1米网格跑了结果网格数飙到近百万求解时间非常长而且很多远场单元纯粹是浪费。换成渐变映射网格后单元数降到十万量级结果反而更稳因为井筒附近的解被充分加密了。4.2 瞬态求解器配置THM耦合这种问题一般是瞬态总时长我设成30年。时间步长列表可以设置成“range(0,0.1,1, 1,5,50)”式的分段格式前期注水初期每0.1年输出一个点后期每5年输出一个点。求解器用“瞬态”在求解器配置里把“高度非线性”选项打开相对容差放宽到1e-3就行太严格的1e-6会让非线性迭代卡死在细枝末节上。在“参数”选项卡里非线性方法建议用“自动牛顿”阻尼因子初始设置0.5如果反复振荡把阻尼因子降到0.2或者0.1。如果模型里渗透率、储水系数这些量存在量级差异再用“分离式逐步求解”顺序是达西压力→温度→固体力学每个子步内做若干次迭代。这个顺序相当重要先求压力再求温度再求位移每一步都有明确的物理依据收敛速度明显好于全耦合牛顿法。4.3 收敛性失败的常见征兆与对策THM模拟最常遇到的情况是前几步还算正常算到第十步突然残差飙升然后就报“找不到解”。我总结下来先要检查这几处如果残差从第一步开始就不降八成是初始条件没有落稳或者边界条件里有数量级冲突。比如压力用了MPa而材料参数里渗透率用了m²导致数值矩阵病态。COMSOL的物理场单位会自动处理但如果通过表达式写了自定义边界就容易把单位搞乱。如果在前几十步收敛之后发散多数是达西速度过大导致对流占主导Courant数超限。这时要用边界层加密并减小时步或者把传热接口改成“全耦合”让对流项和传导项同时迭代。如果网格一变形就出现负雅可比那是几何变形量级超出小变形假设。固态力学模块默认是“小变形”如果你设置了“几何非线性”必须确认大变形理论在实际中是否成立否则就不要勾选。5. 常见问题实录与避坑清单5.1 初始孔隙压力没设成静水第一步就给你回弹这个坑我印象特别深。第一次跑模型时我把初始水压力直接设成了整个域10 MPa的常数结果达西定律一启动整个压力场瞬间往静水压力分布调整在顶部和底部形成了巨大的人为压力波连带着力学场出现初始回弹位移把真正注水的信号全掩盖了。后来把“达西定律-初始值”改成随z坐标线性变化的函数并把初始位移场设为0再配合重力载荷求解一次稳态位移场作为初始条件启动响应基本消失。5.2 边界上的压力振荡与温度怪值在井筒边界上冷水和高温岩石之间温差大如果网格不够密温度场会出现波浪状阶梯。这个在传热问题里非常典型。对策就是在井筒附近加边界层。做了这个之后温度和压力曲线的振荡几乎消失。另外传热接口在“对流”项的Petrov-Galerkin稳定参数上也值得看一眼默认“一致稳定化”建议打开。不打开的话达西速度较高时温度会在边界前后各出现一个鬼峰。5.3 热应力“参考温度”设置的致命错误热膨胀子特征需要一个参考温度T_ref这个值必须是你所设初始温度场在对应位置的值或者是整个模型统一初始温度。我当时为了省事填了293.15K也就是20℃但初始地温在60℃左右于是算出来每个单元都带一个巨大的压缩热应变位移和应力场全部失真。虽然温度场本身没受这个影响但力学结果完全没法解释。正确做法在全局参数里定义一个T_ref20℃或者T_ref平均储层温度再在热膨胀特征里引用这个参数。如果地层温度随深度差异很大更精细的做法是给T_ref设成位置函数但这个在常规案例里不适用。5.4 渗透率固定还是随应力变化大多数教材案例里渗透率是常数但真实地层在注水压差下会发生孔隙压力变化进而改变有效应力渗透率随有效应力减小而升高这就是应力敏感渗透率现象。如果只做教学实例固定渗透率没问题如果做工程预测建议把渗透率写成有效应力的函数比如k k0·exp(a·Δσe)其中σe是有效应力a是应力敏感系数。这个可以用“变量”功能定义并在达西接口里把渗透率填成表达式k(effective stress)。这样三场耦合就升级成了四场联动不过收敛难度也会上一个台阶。6. 后处理、验证与批量参数化的经验6.1 结果验证与网格无关性检查模拟做完第一步不是急着画云图而是做网格无关性验证。我习惯取三种网格粗网格1万单元、中网格10万单元、细网格50万单元对比井底压力随时间曲线和r100米处温度随时间曲线。当两条曲线的差异缩小到2%以内时就把中网格定为标准网格。此外还会找一个半解析解或者简化模型做对标。比如纯地层膨胀问题可以用固体直径与热弹性位移的解析公式渗流部分可以对比Theis井函数解。这些对标看起来繁琐但能快速暴露耦合系数的符号错误和量纲问题。我自己就通过对比Theis解发现渗流模块里储水系数单位写反了一次。6.2 用MATLAB/Python控制COMSOL批量跑参数这个案例里需要扫的参数包括渗透率、注水温度、注水压力、Biot系数。逐一手动改模型参数效率太低我后来直接用LiveLink模块或者COMSOL Java/Python API来批量驱动。在COMSOL中可以把参数扫描做成“参数化扫描”研究步骤但更复杂的多参数组合或者需要外部优化时我会写脚本调用COMSOL的批量计算接口。用Python控制COMSOL的典型思路是在COMSOL里打开一个.mph模型文件通过mphstart或者LiveLink for MATLAB/Python模块连接然后循环修改参数并运行研究最后把结果导出成CSV或图片。前期把模型建立好参数全部设置为全局参数脚本只需要负责循环和收集结果整体自动化程度很高。我也在Linux集群上跑过这类批量任务COMSOL对Linux支持得很成熟可以用命令行comsolbatch -inputfile xxx.mph -batch“solve” -study std1提交作业。批量跑参数时每个模型文件独立计算不会互相干扰。6.3 最后再分享一个小技巧二维轴对称模型里后处理看径向剖面图时有个非常实用的操作在“二维绘图组”里添加“旋转二维数据集”可以围绕对称轴旋转出三维云图效果。虽然本质还是二维解但展示给合作方或写成报告时会直观很多颜色云图一眼就能看出冷锋在径向上的推进形态。另一个技巧是利用“全局定义”里的“非局部耦合”做体积平均。比如我想看注水区半径方向的平均温度随时间变化可以定义range积分或者点平均值省去了在不同时刻截取剖面再导数的重复劳动。特别是做参数扫描时把所有需要的物理量提前定义成全局量导出时就自动带上了。整体走完这套流程后我自己最大的体会是THM模拟不像电热模拟那样“装上模块就能跑”它是真正需要逐场验证耦合项的。先把压力场单独算准再把温度场单独算准最后才把力学场挂上去每一步都有清晰的物理对照这样出了问题也能快速定位。二维轴对称模型绝不是三维模型的降级它是把复杂耦合先盘活、再理清物理逻辑的最佳试验台。