企业官网建设流程全解析

项目标题里“水热力三场耦合”这几个字，放到工程模拟里其实就是固体力学、流体渗流、温度传导三个物理场互相咬合的问题。我前阵子用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)，软件默认r=0处为旋转对称轴。建模时只需要画一个矩形代表地层剖面，然后用“物理场接口”在域内定义固体力学、达西定律、多孔介质传热，边界条件在r=0处自动按对称处理。

1.3 物理场接口的取舍：用内置多物理场还是自定义PDE

COMSOL里做THM耦合有两条路：一条是直接用物理场接口加多物理场耦合节点，另一条是自己在“数学模块”里用系数型PDE甚至一般形式PDE去写三套方程，再用耦合项链接。

我建议绝大多数人走第一条路。因为COMSOL 6.4里已经有相当成熟的“达西定律接口”、“多孔介质传热接口”、“固体力学接口”，再加上多物理场节点里的“多孔弹性”和“热膨胀”，已经能把主要物理勾连起来。内置物理场的好处是边界条件语义清晰，比如达西接口里可以直接设置“压力”边界，传热接口里可以直接设“热通量”，后处理结果也有现成的压力云图、温度云图，省去自己写方程时边界的表达问题。

只有在遇到内置接口没有覆盖的特殊本构，比如考虑渗透率随应力动态变化、或者水合物相变影响孔隙比时，才需要动用自定义PDE。自定义PDE灵活性高，但麻烦也很大，尤其是初边值条件一旦写错，数值发散连报错都找不到方向。后文我会专门给一个渗透率随有效应力变化的扩展思路。

2. 几何建模与材料参数准备

2.1 几何建模步骤与单位统一

模拟对象我取了一个半径500米、深度1000米的地层圆柱体，井筒半径0.1米，回注段取顶部20米。几何画起来就是个4.9999米×1000米的矩形，r方向从0.1到500，z方向从-1000到0。

COMSOL默认单位我们统一用SI制，几何单位设为米。创建几何时，直接添加一个矩形，宽设置成499.9，高设置成1000，左下角坐标(r=0.1, z=-1000)。然后在这个平面域上构建网格。虽然这个矩形长宽比接近500倍，但因为轴对称模型本来就是细长柱体剖面，所以只要网格划分得当，不会影响解的正确性。

这里有一点值得提醒：二维轴对称模型的几何单位会直接影响表达式里的r值。写对的是一个从0.1米开始的域，如果画几何时把半径起点画成0，然后想用“边界条件-井筒”去模拟井筒，就会出现一个半径为零的奇异点，尤其是在达西定律的径向压力梯度里，r在分母处会出现奇异性。所以正确做法是让几何域从井筒半径开始，把井筒作为内边界处理，而不是把旋转轴本身当成井筒。

2.2 材料参数取值与依据

材料参数是这类模拟最容易被质疑的地方，一定要挨个交代清楚。我用的基本参数如下：

参数里最容易引起争议的是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)。取地表压力为0，g=9.81，在地下1000米处，p0数值就是9.81 MPa。很多初次做的人直接在“达西定律-初始值”里填“1e7 Pa”，结果整个压力场从第一秒开始就在向静水压力平衡调整，产生一个人为压力脉冲，干扰了真实注水信号。

温度初始场按地温梯度设：T0(z) = 30 + 0.03×(−z)，在地表是30℃，在地下1000米是60℃。这个梯度大约每百米3℃，符合大陆地壳常见地温梯度。

边界条件分成几组：

• 井筒边界r=0.1：达西压力固定为p_in=12 MPa，高于初始静水压力，形成注水压差；温度固定为T_in=10℃，模拟冷水的注入。
• 远场边界r=500：达西压力按初始静水压力设置，温度按初始地温设置，岩石位移设置为自由或滚动支撑边界，防止刚体漂移。
• 顶部z=0：设为排水自由边界，压力固定为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 网格划分的实战心得

对于轴对称细长模型，最忌讳的是用默认自由三角形网格，因为长宽比太大，剖出来密密麻麻而且单元质量差。正确做法是先选择“映射”网格，把整个矩形剖成四边形，再在r=0.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_ref=20℃或者T_ref=平均储层温度，再在热膨胀特征里引用这个参数。如果地层温度随深度差异很大，更精细的做法是给T_ref设成位置函数，但这个在常规案例里不适用。

5.4 渗透率固定还是随应力变化

大多数教材案例里渗透率是常数，但真实地层在注水压差下会发生孔隙压力变化，进而改变有效应力，渗透率随有效应力减小而升高，这就是应力敏感渗透率现象。如果只做教学实例，固定渗透率没问题；如果做工程预测，建议把渗透率写成有效应力的函数，比如k = k0·exp(a·Δσe)，其中σe是有效应力，a是应力敏感系数。这个可以用“变量”功能定义，并在达西接口里把渗透率填成表达式k(effective stress)。这样三场耦合就升级成了四场联动，不过收敛难度也会上一个台阶。

6. 后处理、验证与批量参数化的经验

6.1 结果验证与网格无关性检查

模拟做完第一步不是急着画云图，而是做网格无关性验证。我习惯取三种网格：粗网格1万单元、中网格10万单元、细网格50万单元，对比井底压力随时间曲线和r=100米处温度随时间曲线。当两条曲线的差异缩小到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模拟不像电热模拟那样“装上模块就能跑”，它是真正需要逐场验证耦合项的。先把压力场单独算准，再把温度场单独算准，最后才把力学场挂上去，每一步都有清晰的物理对照，这样出了问题也能快速定位。二维轴对称模型绝不是三维模型的降级，它是把复杂耦合先盘活、再理清物理逻辑的最佳试验台。