COMSOL两相流THM热流固耦合建模：从方程到收敛调试全解析

简介：一份基于COMSOL的两相流THM热流固耦合数值仿真模型，面向油气渗流、地热开发与地下多相流模拟领域的工程师和科研人员，也适合用作高校及科研机构数值仿真教学课程中的进阶案例。模型以饱和水和油的两相流动为对象，同时耦合热传导、流体渗流与固体力学过程，支持按实际需求调整物理场、材料参数、边界条件或几何结构，还可进一步扩展至水与超临界二氧化碳的两相流体系。压缩包共14个文件，包括8个Word文档、2个HTML文件与4张图片，整体大小3.91MB；Word文档系统介绍两相流THM耦合模型的构建流程、控制方程、求解器配置与后处理解读，HTML文件便于在浏览器中快速浏览模型说明，图片则展示了几何建模、网格划分与结果云图等关键环节。资源已有130人浏览学习。文档中对模型收敛性、计算速度与鲁棒性也进行了分析，并提供了参数化可调方案，方便读者移植到自己课题或深入二次开发。

1. 两相流THM热流固耦合模型：为什么说它是储层仿真的“联席解题”

两相流THM热流固耦合模型，指的是在COMSOL中把饱和水与油的渗流、多孔介质传热、岩石或土体的力学响应放进同一个求解框架里同步计算。它能回答“注水驱油时油前缘推到哪、地层温度怎么变、孔压抬升会不会让盖层变形”这类单一物理场模型答不出的问题。地热回灌、聚合物驱、CO2埋存等场景本质上都吃这套耦合：温度改变流体粘度，流动带着热量跑，孔隙压力又改写有效应力，最后应力变形再回头影响渗透率。这篇笔记直接面向要复现模型的人，把物理场拆解、参数表、边界条件、求解器配置和踩坑点逐条写清楚。新手能按步骤搭出收敛的初版，熟手也能在参数与失效模式上找到自己翻过车的影子。

2. 先把三场拆开：油水两相、传热与固体力学各自的方程和耦合通道

2.1 油水两相流动：饱和度、毛细压力与相对渗透率的配合

COMSOL里做储层两相流，常见做法是用多孔介质两相流接口，变量是水相压力和含水饱和度，油相压力通过毛细压力关联。核心方程是两相各自的达西速度，水的速度写为：

$u_w = -\frac{k_{rw}}{\mu_w(T)},K,(\nabla p_w + \rho_w,g,\nabla z)$

油相同理，只是用油的相对渗透率、粘度和密度。饱和度由质量守恒方程推进：

$\phi \frac{\partial S_w}{\partial t} + \nabla \cdot u_w = q_w$

这里最容易忽略的是相对渗透率和毛细压力曲线的端点处理。油水两相模型的默认做法是给定一条从束缚水饱和度到残余油饱和度的归一化曲线，比如Brooks-Corey形式。水的相对渗透率一般写成：

$k_{rw}=S_{we}^{,(2+\lambda)/\lambda}$

其中$S_{we}$是归一化有效水饱和度，$\lambda$是孔尺寸分布指数。油相相对渗透率在形式上与油饱和度相关，但在靠近残余油饱和度时会趋于零。COMSOL允许直接以插值表形式导入实验室测得的曲线，也可以输入解析表达式。对刚起步的模型，我建议先用光滑的解析曲线，不要直接把现场报告的离散点塞进去，否则会在端点处出现折角，非线性求解器容易卡住。

两相流接口里还需要定义毛细压力。它决定油的吸入与排出行为，直接影响饱和度空间分布。注意毛细压力是饱和度的函数，而且存在迟滞，多数简化模型不考虑迟滞；如果你做的是注水驱油这种驱替过程，使用排水主曲线即可，先别加迟滞。另一个常被忽略的点：初始饱和度不要让任何一相刚好位于残余饱和度端点，否则该相相对渗透率是零，雅可比矩阵可能奇异。初始水饱和度建议取在束缚水饱和度以上5%左右，例如束缚水饱和度为0.2，初始就设0.25。

2.2 传热子问题：有效热性质与对流热传递

多孔介质传热的稳态与瞬态方程，在COMSOL中默认按等效体积热容处理。等效热容近似为：

$(\rho c_p){eff} = \phi,\left(S_w \rho_w c{p,w} + S_o \rho_o c_{p,o}\right) + (1-\phi),\rho_s c_{p,s}$

等效热导率类似地是水和油按饱和度的加权，再加上骨架贡献。这里有一个新手容易踩的误区：直接把固相热导率和流体热导率分开设置，然后让界面自动“平均”，但多孔介质传热接口默认用的是体积平均模型。如果你用的是“多孔介质传热”接口，它已经把各相按饱和度混合好了，不要再额外加一个“层流”或“流体传热”接口去算流体域，否则热容会被重复统计。

对流项是热-流耦合的主通道。达西速度由流动接口给出，进入传热方程的对流项$\rho_f c_{p,f} u \cdot \nabla T$，这个耦合在COMSOL的多物理场节点里是自动建立的。如果你是手动搭耦合，这可能是唯一一个两个物理场均需保留速度变量的地方：流动接口输出速度场，传热接口读取它并计算热量输运。

从参数敏感性来看，水比热4200 J/(kg·K)，油通常只有2000多一点，骨架约900。在典型含水饱和度0.3的情况下，等效热容主要由骨架承担。温度场传播速度远慢于压力波速，但会明显滞后于饱和度前缘。因此时间步长要兼顾两个特征时间：压力传播快，需要小步或隐式格式；热前缘慢，允许大步。COMSOL的BDF时间步进会自动调节，但如果初始步长给太大，饱和度前缘可能被“跨过去”，导致温度场和饱和度场的空间分布错位。

2.3 固体力学子问题：有效应力、热应变与孔压的三种耦合

固体力学在THM中承担的是“看到应力变形”的角色。经典做法是采用孔隙弹性材料模型。有效应力写成：

$\sigma'{ij} = \sigma{ij} - \alpha_B, p_w, \delta_{ij}$

这里的$\alpha_B$是比奥系数，介于孔隙度和1之间。对砂岩可取0.6到0.9，对致密岩石取0.3左右。COMSOL材料模块里有“孔隙弹性”子节点，设置比奥系数后，系统会把孔隙压力作为体载荷自动纳入力学方程，不需要你再手写一遍。

热应变项则来自温度变化：$\epsilon_T = \alpha_T,(T - T_{ref})$。这个量在热-力耦合中占主导，尤其当地层温度从80°C降到20°C的冷水回灌场景，冷缩量可能远超注水引起的孔压形变。如果你只算力学不写热应变，位移结果可能小一个量级，那样的THM模型没有实用价值。

三场耦合的通道可以概括为：温度场→流体粘度和密度→流动行为；流动→对流换热→温度场重分布；孔压与温度→有效应力→变形；变形→孔隙度与渗透率→流动参数变化。第四条通道最复杂。初学者建议先忽略变形对渗透率的反馈，把孔隙度和渗透率设为常数。等模型跑通，再按Kozeny-Carman公式或经验关系把渗透率做成体积应变的函数，逐步加复杂度。

2.4 建模前必须决定的三个边界条件

第一个是流动边界。注入端用压力边界还是流量边界，差别很大。如果目标工况是定流量注入，尽量用通量边界，让压力自由响应；如果用高压边界而实际流量还未知，压力可能被抬得过高，力学模块容易崩溃。第二个是热边界，通常注入温度固定，远边界设为地层温度或热绝缘。第三个是力学边界。底部建议固定法向位移，顶部自由或加有效应力载荷；侧向不要全部固定，否则热应变引起的侧向变形被锁死，应力结果过于保守。

这三组边界必须在建模前定下来，因为它们决定模型的刚度矩阵性质。流动压力边界会直接影响孔压载荷，力学约束又决定应力场能不能释放，两者耦合不好，常常出现后续章节要讲的“压力尖峰”问题。

3. COMSOL建模实操：从几何到物理场再到求解流程

3.1 几何简化与网格划分策略

对于THM两相流，几何不必一开始就做复杂层状构造。先用一个矩形剖面或简单的二维平面模型跑通流程，是最省时的做法。以二维矩形剖面为例，长取400米，高50米，左侧模拟一口注水井，右侧为生产边界或定压边界。这样的几何在COMSOL里叫“几何1”，直接用矩形加上点或边来表示井位。井半径在二维模型里通常用一个点源来近似，或者用一个小半径圆弧施加压力。

网格划分是两相流THM的隐形天花板。饱和度前缘是一个移动的锐边界，网格太粗则前缘数值扩散严重，油水过渡带可能被抹掉；网格太细则每个时间步的雅可比组装和线性求解代价高。常见做法是在注水井附近使用加密区，用“边界层”或“扫掠”网格，让井周围网格尺寸在0.5米量级，远端网格放大到5米。对二维模型，网格数控制在3万到5万之间就能兼顾速度和精度。

COMSOL中可以用“网格序列”自定义多级加密。要留意两相流的Peclet数，它表征对流与弥散的相对强度。网格尺寸过大，饱和度方程表现为强对流主导，BDF格式容易产生数值振荡。这与计算流体力学里的Courant数类似；工程经验是让网格尺寸小于压力扩散特征长度的1/5。如果前缘宽度本身只有几米，那网格就必须在那一带给到0.5米以下。

提示：不要在第一次运行时就把网格加密到极致。先用粗网格跑通，确认物理趋势合理后再逐步加密到网格无关解。

3.2 材料参数表：六个必填项与两个易错项

材料参数是多场耦合的数据地基。下面是初学者可直接抄的参数表，更新为应用相关的数据：

两个易错项值得专门提醒。第一个是油相粘度。很多人习惯在材料库里直接选“油”，然后默认粘度是一个定值，结果注入冷水后测温升高但流速几乎不变，完全丢失THM的核心耦合效应。正确处理是给粘度定义温度插值函数，常见做法是使用Arrhenius型或Vogel方程，在COMSOL的插值函数里导入温度和粘度的数据点。第二个是比奥系数。默认设置里如果不启用孔隙弹性，固体力学里根本没有孔压项，输出的位移和应力就跟注水压力没有半点关系，后面排查起来非常费劲。

3.3 物理场接口与多物理场耦合节点

在COMSOL中新建模型时，三个关键物理场接口分别是“多孔介质两相流”“多孔介质传热”和“固体力学”。不同版本的名称可能略有差异，但本质是同一个组合。先把三个接口都加到组件中，再检查右侧的“多物理场耦合”节点。

两相流接口中，你需要指定水相和油相的属性，包括密度、粘度和相对渗透率模型。相对渗透率模型选“解析”或者“插值表”。插值表需要提前准备好饱和度与kr和pc的对应数据。传热接口里，流体和固体的热性质按材料定义，孔隙率会自动从流动接口读取。力学接口中，设置弹性模量、泊松比和热膨胀系数，并在固体材料中钩选“孔隙弹性”子节点，填入比奥系数。

多物理场耦合节点是THM的灵魂。在COMSOL中，多孔介质传热与多孔介质两相流之间存在一个“多孔介质传热”与“流动”的耦合；固体力学与两相流之间通过孔隙弹性材料节点关联；固体力学与传热之间通过热膨胀节点关联。建议按默认的自动耦合建立节点，不要手动去添加边界载荷形式的孔压，那不仅容易遗漏，还可能在力学方程里重复加入孔压项。

如果版本中的自动耦合没有生成孔隙弹性节点，退回半手工方案：在力学接口中增加一个“体载荷”，体积力的表达式设为孔隙压力的梯度，例如$F_x = -\alpha_B \frac{\partial p_w}{\partial x}$，并加上热应变。这是常见做法，但要注意体载荷作用的是单元质心处的压力梯度，在压力梯度剧烈的井附近容易出现数值振荡。

3.4 初始条件与分析设置

初始条件对多场耦合模型的成败影响极大。压力和饱和度的初始条件必须与边界条件协调。常见的错误是初始压力设为静水压力分布，而入口边界是定压注入，两者在第一个时间步就产生一个压力跳变，井点附近的网格承受极大通量，非线性迭代直接发散。

我一般先把初始压力场设成与出口边界相同的均匀压力，或者用一次稳态单相流动计算得到初始压力分布，再以此为初值做两相流THM瞬态计算。饱和度的初始值则设定为束缚水饱和度以上几个百分点，这样水相和油相相对渗透率都是正数。传热初始温度设定为地层温度，入口温度低于或高于地层温度。

时间步长的设置也值得讲究。如果研究时间为一年，合适的最小步长可以设为1秒，最大步长设为3600秒。BDF时间步进会在非线性迭代失败时自动退回，但最大步长设得太大，会跳过早期的饱和度前缘形成阶段。建议先跑前半天时间步，观察饱和度剖面是否随时间单调推进。

3.5 在COMSOL中跑起瞬态研究的流程

实际操作流程按以下步骤执行：

1. 建立矩形几何，设置注水井点与生产井点。
2. 定义全局参数，包括注入压力、地层温度、孔隙度和渗透率。
3. 添加三个物理场接口，按照2.1到2.3节设置方程来源。
4. 建立多物理场耦合节点，确认孔隙弹性与热膨胀被激活。
5. 设定初始值，压力写p0，饱和度写S0，温度写T0。
6. 网格划分并加密井附近区域。
7. 研究设置为瞬态，时间范围从0到3.15e7秒（一年），步长设置上限3600秒。
8. 求解器用默认的“PARDISO”或“MUMPS”，如果内存不够再换“迭代求解器”，并设置预处理为“左预处理”。

第一次跑的时候，先关掉力学场快速验证流动与温度场，再把力学场打开，这样出问题能快速定位。利用参数的复用，可以编写MATLAB脚本批量修改参数，具体做法在第5章展示。

4. COMSOL两相流THM模型排查：四个让人翻车的收敛坑

两相流THM的收敛问题几乎不能用“调大阻尼”四个字一言蔽之。这里记录四个高频故障，原因是它们在项目里反复出现，现象、机理和解决路径都很有代表性。

4.1 非线性残差停在1e-3不动，时间步反复退回

现象：瞬态求解时，非线性迭代前期下降很快，到残差1e-3左右就开始水平蠕动，散步次数达到上限，时间步长被退回。日志显示“找不到一致的初始条件”或“最大迭代次数达到，未收敛”。

原因：最常见是相对渗透率曲线在端点处不平滑。实测曲线数据在饱和度接近束缚点或残余油饱和度时出现跳变，解析式里的指数函数在端点数值急剧变化。另一种原因是某个相饱和度在局部网格中越过端点，导致该相相对渗透率变成负值或NaN，雅可比矩阵出现坏条件。

解决：给相对渗透率和毛细压力曲线添加端点磨平，常见的做法是让曲线在端点的导数连续。检查模型中是否有区域油相饱和度低于残余油饱和度，这种情况不在边界处加流线扩散就很容易振荡。求解器设置方面，把“非线性容差”从默认的1e-3放宽到5e-3，或将“最大迭代次数”从25提高到50，效果有时立竿见影。但要注意，放宽容差不等于放弃计算，要观察累计物质平衡误差是否仍小。

4.2 饱和度在局部变成负值或大于1，出现“饱和度息肉”

现象：后处理中含水饱和度云图出现局部小斑块，值超过1或低于0，时间步不断缩小。日志中提示“检测到负饱和度”，计算速度骤降甚至最后中断。

原因：饱和度方程具有强对流性，使用标准伽辽金有限元离散时在没有流线扩散的情况下，饱和度前缘出现数值振荡。网格太粗或时间步太大，都会加剧。

解决：在求解器稳定化选项中给两相流接口启用“流线扩散”。这是一个经典做法；它可以抑制前缘数值振荡但也会轻微抹平前缘。更本质的解法是加密饱和度前缘网格并把时间步长上限调低到能与网格匹配合适的范围。工程经验：前缘所在区域的网格Peclet数控制在2以下，时间步与网格尺寸要满足“饱和度前缘单步推进不超过一个网格”的经验准则。

4.3 注入井点压力尖峰，第一或第二步即发散

现象：计算刚开始一两步，近井区域压力从1e7 Pa飙升到1e11 Pa，然后求解中断。位移场同时出现异常的凹凸变形。

原因：初始压力与边界条件不协调。入口边界直接设为2e7 Pa的高压，而整个初始场是1e7 Pa，第一个时间步在注水点形成巨大的压力梯度，流动方程的非线性迭代无法解析这种突变。另一个原因是网格过粗，井点附近没有边界层，压力梯度无法被准确解析。

解决：将初始孔隙压力场设置为与入口边界压力近似一致的场，或先做一次稳态单相流动获得压力初始场，再开始两相流THM的瞬态计算。另一种稳妥做法是采用“缓启动”注入：时间函数里前1000秒压力从出口压力线性上升到目标注入压力，避免阶跃。网格层面在井点处切分一个小圆或小矩形区域，施加至少三层边界层网格。

4.4 位移场与注水压力完全无关，力学耦合“形同虚设”

现象：把注入压力从1e7 Pa提高到5e7 Pa，结果位移场几乎不变。检查发现应力分布和孔压没有任何相关趋势，仿佛力学模块没有读入压力数据。

原因：在固体力学接口中没有启用孔隙弹性材料节点，或者多物理场耦合节点里没有把两相流压力接口和力学接口关联起来。很多版本默认不会自动感知，需要你在多物理场节点里手动添加“孔隙弹性”耦合，并指定压力变量。如果材料节点手动覆盖了应力分量，也会发生同样情况。

解决：在固体力学的材料定义中启用“孔隙弹性”子节点，设置比奥系数。同时确认多物理场节点中出现“孔隙弹性耦合”，并在其中选择两相流压力变量。完成后重新求解并检查井点附近孔压与应力梯度的量级。

5. 自由修改与扩展：从油水两相到水+超临界CO2及批量工况

5.1 改造成水+超临界CO2两相流的三个关键改动

标题里提到的水+超临界二氧化碳两相流版本，很多人拿到油水模型后直接替换材料库里的油，这是翻车率最高的操作。超临界CO2的密度和粘度随温度压力变化极大，在储层条件下密度可能从150 kg/m³变到700 kg/m³，不能再用材料的常数描述，甚至不能用简单的线性插值。

第一个改动是密度和粘度。常见做法是使用Span-Wagner状态方程或Peng-Robinson状态方程预先计算出温度压力网格上的物性表，再以插值表导入COMSOL。把油相材料替换为“CO2相”时，不要保留原来油相材料里的常密度与常粘度，否则重力驱替行为完全是错的。第二个改动是毛细压力曲线。CO2与水的界面张力与油水不同，等值吸/排曲线也不同；如果你没有实验数据，可以用基于界面张力的无因次毛细压力模型暂替。第三个改动是初始饱和度的设置。由于CO2密度明显低于水，浮力不可忽略，模型需要竖直剖面才能体现CO2的向上聚集，不能用水平面模型替代。

超临界CO2的流动前缘通常更不稳定，因为粘度比差异大，黏性指进更加明显。要保证网格分辨率足以捕捉前缘的细观结构，时间步长也要相应缩小。如果模型区域较大，建议先做二维剖面验证，再考虑三维化的扩展。

5.2 按自己的需求自由修改：三种典型改法

自由修改绝不是打开源文件改两个参数就交差。要对物理场做增减，用参数化控制更高效。三种典型改法具有普遍意义。

第一种是只保留流动与传热，去掉力学场。直接在“研究”设置中取消“固体力学”对应的求解场，或将饱和度和应力解耦。去掉力学场后，计算量能显著下降，因为刚度矩阵规模大幅缩小。第二种是只保留流动与力学，温度场用固定值不参与求解。做法是把传热接口的方程关闭，同时将粘度设为常数，这样THM退化成HM模型。第三种是增加化学场，考虑矿物溶解沉淀引起的孔隙度和渗透率变化。这种扩展通常需要额外引入溶质输运方程，并将化学反应速率表达式耦合到孔隙度变化项中。

在COMSOL中实现开关，推荐在“参数”表里添加一个控制变量，例如use_mech = 1，然后通过表达式使孔隙弹性耦合项乘以这个开关。这样不需要删除物理场接口，只需要切换参数即可在模型家族里共享同一个几何和网格设置。

如果自己会写COMSOL的“Application Builder”或“LiveLink for MATLAB”，可以用更系统化方式管理。下面是基于LiveLink for MATLAB实现的参数化扫描，适合批量修改注入压力并输出出口温度或累计产油量：

% 假设COMSOL模型已打开并定义了参数 p_inj 和 S_w0 model = mphopen('thm_2phase.mph'); p_inj_list = [1.1e7, 1.2e7, 1.3e7]; % 不同注入压力 T_out = zeros(length(p_inj_list), 1); for i = 1:length(p_inj_list) model.param.set('p_inj', num2str(p_inj_list(i))); model.study('std1').run(); T_out(i) = model.result.numerical('intop_T_out').getReal(); fprintf('p_inj=%.1f MPa, T_out=%.2f K\n', ... p_inj_list(i)/1e6, T_out(i)); end

该脚本的逻辑很直白：对每个注入压力，重设模型参数、运行研究、用积分算子读取出口温度。几个关键参数说明：p_inj必须是模型已定义的全局参数，否则set会报错；intop_T_out需要提前在模型结果中定义好积分算子，并且位于出口边界上；std1是研究的名称，不同模型里可能是std2或自定义名称。运行这个脚本之前要确认COMSOL的LiveLink for MATLAB接口已经正确连接并指向当前打开的模型文件。

5.3 汇总一并列出参数扫描对收敛数据的验证方法

做多工况扫描时，除了关注数值结果，还要让每次运算的结果通过检验，否则一个算例的噪声会被当成物理信号。我的常见做法是把每次扫描的质量平衡偏差、迭代次数和守恒积分记录下来。通过对比同一网格下不同注入压力下的累积注水量与油相产出量，可以快速辨别模型的鲁棒性，而不必逐个查看云图。

6. 验证模型对的三个习惯：从“能跑”到“敢用”

模型能跑通不等于结果能用。我养成了三个固定验证手段，在这里完整提供。

第一个习惯是网格无关性检验。对同一工况，用三套不同密度的网格计算，比如5万、10万和20万单元，对比生产端累计产油量和温度突破时间。如果温度和产油量变化在3%以内，认为网格足够；如果偏差大，以最细网格为准，同时检查是否存在局部加密盲区。这个检验必须在第一次完成正式方案前做掉，不然之后填任何网格都会成为“黑匣子”里的一块绊脚石。

第二个习惯是极限退化测试。把温度场设为常数、力学场关闭，把两相流模型退化成普通的油水两相渗流，看它能否匹配单相或解析解的流动特征。比如一维注水驱油可以用Buckley-Leverett解析结果对比。如果连退化后的结果都对不上，那么THM全耦合的可信度无从谈起。这个步骤看起来很笨，却能在几分钟内帮你看穿模型是否在物理上站得住。

第三个习惯是守恒与趋势检查。在COMSOL后处理里对出入口做质量积分、热量积分，确认总累计注入量减产出量等于储量变化量。同时在参数扫描结果中画出“温度突破时间随注入压力”的曲线，观察趋势是否单调。有一次我对比扫描结果，发现温度突破时间没有随注水量增加而显著提前，排查后发现是入口边界误设成了无热流条件，热量交换根本没发生。这种问题靠云图很难发现，靠趋势曲线一眼就看穿。

我的血泪教训是：不管多急着出图，先把这三个验证动作练成习惯。模型自己没有对错，只有你还没发现的不一致。希望帮到你。

本文还有配套的精品资源，点击获取