企业官网建设流程全解析

简介：这是一份基于COMSOL的两相流THM热流固耦合模型资源，面向从事油气渗流、地热开发或CO2封存模拟的工程师与研究生，可模拟饱和水与油两相流动下的温度场、流场和应力应变三场耦合响应，也可扩展为水与超临界二氧化碳体系，并支持按自身需求修改。压缩包共14个文件，大小3.91MB，其中8个docx提供模型原理、方程推导、操作步骤与不同场景说明，4个jpg直观展示几何建模、网格划分及结果后处理效果，2个html便于在浏览器中快速查阅相关图文内容。资源目前已有130人学习。模型具有收敛性好、计算速度快、鲁棒性强的特点，文档中对两相流热流固耦合的构建思路、参数设置及常见问题均有展开，读者可结合文档与示例快速上手，显著减少自行调试模型的时间。

1. 一个能直接跑通的 COMSOL 两相流 THM 模型：饱和水与油的热流固耦合，稳定不玄学

地下储层模拟里最常翻车的组合是什么？就是两相渗流、温度场和应力场三个物理场叠在一起。油田热采、地热回灌、干热岩换热这类问题，实际跑 COMSOL 时往往不是方程搭不起来，而是算到一半压力振荡、饱和度前锋数值扩散，或者换一组边界条件后直接发散，一改参数就全盘崩掉。这套两相流 THM 热流固耦合 COMSOL 模型，把饱和水和油的两相流动、温度场、应力场放进同一个求解框架里，收敛性好、计算速度快、鲁棒性好，岩石骨架和流体物性都留了改动余地，也支持替换成水与超临界二氧化碳两相流版本。适合做储层注采模拟、敏感性分析和方案对比的从业者，拿到手不需要从零搭物理场，直接跑通再二次开发。

2. 从物理场到数值构型：THM 到底在解哪几组方程

2.1 三场耦合的作用关系：压力、温度和应力谁也离不开谁

这类模型的核心不是“三个场同时算”，而是三个场之间存在真实物理反馈。温度升高会让水的粘度下降、密度减小，流速分布跟着变；流动压力变化改变有效应力，骨架受力状态又反过来挤压孔隙，造成孔隙度和渗透率变化；骨架变形加上热膨胀，位移场和应力场也就不再是单纯力学问题。比如热采注热水时，温度前锋先到，先改写粘度，近井地带的油相流动性变好，压力扩散速度也变了，而压力变化又让有效应力下降，可能诱发骨架微变形，渗透率再变。这个反馈环如果不建全，模型和现场实测数据对不上。

在 COMSOL 里实现时，常见做法是让三个物理场接口同时参与计算。流体流动部分由 Darcy 定律接口提供流速和压力场，传热部分用多孔介质传热接口，固体力学接口计算位移和应力。多物理场耦合节点把三者接起来，温度进流场是“粘度温度相关”，流场进力学是“孔隙压力作为体内载荷”，力学进流场是“孔隙度随应变更新”。这套模型的搭建思路就是这种强耦合，不是把三个场分开算再手工传递结果，那样只能做单步近似。

2.2 油水两相对饱和度场的处理：为什么常走“压力方程加系数型 PDE”

COMSOL 内置的 Darcy 定律接口默认是单相流动，直接拿它算饱和水和油的两相驱替并不够用。常见做法是保留 Darcy 接口解孔隙压力，再额外解一个饱和度输运方程，把油相和水相的速度分别用相对渗透率折算出来。方程形式大致是：

∂(φρ_w S_w)/∂t + ∇·(ρ_w v_w) = 0

油相类似，总速度是两相速度之和，两相压力差由毛管压力关联。这个实现方式的好处是自由度开放，你想换相对渗透率模型、毛管压力曲线都直接改表达式，不用重新搭物理场。模型里默认用的相对渗透率一般是 Corey 曲线，毛管压力可以用 Brooks-Corey 或 van Genuchten 表达式，矿场有实测曲线时替换成插值表也很方便。

需要注意饱和度有一个硬约束：S_w + S_o = 1，且 S_w 必须夹在束缚水饱和度和残余油饱和度之间。初值给不好，求解器第一轮迭代就会碰到物理量越界，随后发散。模型文件里对饱和度做了限制器处理，这也是它稳的关键之一。

2.3 求解策略：全耦合牛顿迭代与时间步选择的边界

逐场交替求解在这个问题上很容易卡住，因为非线性太强。这个模型走的是全耦合牛顿迭代，也就是把压力、饱和度、温度、位移放在同一个未知向量里做雅可比矩阵求解。全耦合求解内存占用更高，但胜在鲁棒。查阅模型求解器设置时会看到阻尼因子默认大概在 0.9 附近，这个值意味着每步迭代做适度阻尼，防止压力场或者饱和度场一步跨过头。

时间步进方面，模型采用 BDF 格式，常见是二阶，自动时间步进开启后能根据非线性残差调整步长。对于注采模拟来说时间步长策略比网格尺寸更容易影响速度：初始步长给得太大，注入早期饱和度前沿会瞬间冲过两三个网格距离，造成数值震荡。实际使用时我会把初始时间步控制在峰值注入量到达一个网格所需时间的四分之一以内，后面再放开自动步长。下一章会展开说具体参数改哪里。

3. 上手实操：从打开模型到改出属于自己的一套算例

3.1 打开模型后先看全局参数，不要急着点运行

拿到这套 COMSOL 模型文件，第一件事不是运行，而是打开全局参数列表。全局参数是这个模型能否二次开发的关键。参数以名称、表达式、单位为一组列出来，储层性质、流体物性、注入条件全部集中在这里，不需要钻进每个物理场的边界条件里一个个找。我习惯按表 3-1 逐项过一遍，确认默认值代表的工况是否贴近自己的问题。

需要注意参数之间不是独立的。孔隙度改了，渗透率大概率要跟着改，否则 P 波传播或者压力扩散特征会失真。骨架的杨氏模量和热膨胀系数改完之后，应力场量级会明显变化，这时要检查网格变形是否过大。

3.2 参数修改三步走：全局参数、物理场边界、求解器

修改参数的路线应该是“全局参数优先，物理场次之，求解器兜底”。绝大多数工况调整只需改全局参数。我一般用 COMSOL 的 Java API 做批量修改，比在图形界面里逐个点快很多，尤其是要做多组算例的时候。参考脚本如下：

// 用 COMSOL Java API 批量修改全局参数并运行 model.param().set("poro", "0.22"); // 孔隙度从0.20改为0.22 model.param().set("perm", "5e-15[m^2]"); // 渗透率改为5mD量级 model.param().set("Sw_init", "0.25"); // 初始含水饱和度0.25 model.param().set("T_inj", "373.15[K]"); // 注入温度100摄氏度 model.study("std1").run();

set 方法第一个参数是全局参数名，第二个参数是含单位的表达式。COMSOL 里单位必须写进表达式字符串，写成 "1e-15" 会被当成无量纲数，压力场量级会直接出错。std1 是对应的瞬态研究，如果文件里研究标识不同，要按实际名称替换。改完参数先不要急着看结果，跑完第一步后回到收敛日志确认残差曲线是否平滑。

3.3 先跑稳态给初值，再开瞬态注入，否则第一轮就翻车

这类两相 THM 模型最容易被忽略的一个环节是初始条件初始化。直接设置好注入流量就运行瞬态，注入井附近会出现压力突变，饱和度场在第一个时间步就会剧烈震荡。正确的流程是先跑一个稳态研究或者“无注入”的预处理步骤，让孔压场和应力场在初始载荷下平衡，再以此为初值打开注入。

模型文件里如果已经内置了这种两步研究策略，直接按顺序执行即可。如果没有，我通常手动加一个稳态研究步骤，关闭注入边界条件，让它只计算地应力与初始孔压平衡。这一步跑完再换成瞬态，后继就很顺。初始饱和度场不能简单设成整个区域一个常数，最好在井附近设置一个小范围的润湿相分布，模拟毛管力作用下近井含水率提高的现象，可以有效防止饱和度初值触碰限制器边界。

3.4 怎么判断“跑通了”：读取收敛日志和关键点动态

COMSOL 求解过程中的“日志窗口”不是摆设。全耦合 Newton 迭代每一轮的阻尼因子、迭代次数、残差大小都在里面。经验是阻尼因子稳定在 0.9 以上、单时间步迭代次数不超过四次，说明当前时间步相对宽松；如果阻尼因子掉到 0.1 以下还不停，基本是边界条件或物性表达式出了问题。除了收敛日志，还可以在模型里预先设置注入点压力探针和采出端温度探针，运行结束后立即查看这段时间序列曲线。压力曲线如果呈平滑上升后趋缓，温度曲线出现明显的热突破拐点，结果基本是可信的。曲线出现锯齿状抖动，就要回到时间步设置或饱和度限制器上找原因。

4. 二次开发与自由修改：从油水两相换到超临界 CO₂

4.1 可以放心改的部分，和需要谨慎对待的红线

这个标称“可以根据需求自由修改”的模型，修改区域分两类。第一类是全局参数、边界流量、注入温度、初始饱和度，直接改值就行，不触及模型结构。第二类是改造物理场接口，比如把油相换成超临界二氧化碳、换相对渗透率模型、把二维模型改三维，这些操作是在原有骨架上升级，需要保证新增表达式里的所有变量都有定义。红线在于饱和度方程自带的限制器设置不能取消，人工扩散系数也不要随意清零，否则饱和度前锋会变成陡峭间断，数值直接失稳。

4.2 改成水+超临界二氧化碳两相流版本：物性表与压缩性差异

标题里提到“水+超临界二氧化碳两相流版本另议”，说明这种替换是常见需求。把油相替换为超临界 CO₂ 并不只是改一个粘度值。CO₂ 在储层温度和压力范围内密度变化剧烈，超临界条件下可能从 100 kg/m³ 跳到 700 kg/m³ 以上，无法用常数密度搞定。实现上我会在全局参数里换掉 mu_o 的表达式，改成温度压力的插值函数。常见做法是用密度表和粘度表导入 COMSOL 插值节点，温度和压力作为自变量，密度、粘度为因变量。毛管压力曲线也要换，CO₂-水系统的界面张力比油水系统低很多，毛管压力数值要按实验数据修正。

模型里原来的油相方程不需要大改，因为控制方程形式都是“饱和度输运+相对渗透率”，换的只是物性表达式和密度计算方式。但要注意 CO₂ 的压缩性问题，压力方程的求解难度会上升。建议把压力求解的阻尼因子调低到 0.7 左右，初始时间步再缩小一半，先让压力场稳定再逐步恢复步长。

4.3 几何替换与网格重建：先做尺寸检查再做局部加密

如果想把默认模型改到自己的工区上，几何可能需要替换。模型一般由参数化尺寸控制，比如井半径、储层长度、模拟区域厚度。改几何参数之后必须重新检查网格质量，特别是注入井附近的边界层网格。COMSOL 自动重新划分网格时，如果几何尺寸变化过大，原来加密区的网格密度可能不再匹配流速场分辨率。我一般先跑一版较粗网格快速检验模型是否稳定，确认没问题再做边界层加密和网格尺寸细化。

4.4 用参数扫描批量产出敏感性曲线

这个模型结构干净，非常适合直接挂参数扫描研究。比如考察注入温度对热突破时间的影响，定义 T_inj 从 353.15 K 到 393.15 K 取五个值，扫描输出累计产油量或热突破时间。使用 COMSOL 的 Parametric Sweep 功能时，不建议把每组扫描拆成多个独立研究，那会重复编译模型，效率低。正确的方式是在研究设置里添加参数扫描节点，扫描列表填参数名和取值列表，一次运行，结果存在同一个数据组里。后处理时用一维绘图组把所有曲线叠在一张图上，比单独跑五次再导数据省太多时间。

5. 避坑记录：从这套模型里踩出的五个常见问题

5.1 第一轮瞬态计算发散，残差直线上升

现象：设置好注入条件后点运行，迭代日志显示残差在第一到第三步内快速增大，求解器报错退出。原因通常是初始条件与边界条件冲突。比如初始含水饱和度全区域 0.15，但井边界立即注入纯水，饱和度突变过大。解决：先运行稳态预处理步骤平衡孔压和应力场，再修改 Sw_init 使得井附近有一个渐变的初始饱和度分布，最后才打开瞬态注入。从那以后我拿到任何新模型，第一步都是检查有没有初始化研究。

5.2 注入井压力曲线出现锯齿状振荡

现象：结果文件里注入点压力不是平滑上升，而是每个时间步出现上下跳动，幅度随时间增大。原因：时间步长过大，饱和度前锋在单个时间步内穿越了多个网格单元，相对渗透率在相邻两个时间步之间跳跃。解决：把初始时间步缩小到原来的五分之一，时间步进方法保持 BDF，回到一阶先跑一段再切换二阶。此问题不影响最终稳态值，但过程曲线失真，不能用来做压力历史拟合。

5.3 把油相改成 CO₂ 后求解器每一步都卡满迭代次数

现象：替换 CO₂ 物性后，全耦合迭代始终无法降到残差阈值以下，每个时间步都在“最大迭代次数达到”边缘警告。原因：CO₂ 的密度和粘度随压力的导数过大，数值上的源项刚度过高。解决：改用插值表提供物性，避免原始表达式在临界点附近出现剧烈梯度；同时把求解器的阻尼因子降到 0.7，并打开“自动步长”下限限制。这一步做完，迭代次数基本回到油相工况的水平。

5.4 网格加密后结果反而比粗网格更不平稳

现象：把整体网格从粗换细之后，饱和度场出现局部扰动，看起来数值噪音更多。原因：网格加密改变了局部单元尺寸，但时间步没有同步缩小，前锋通过细网格时单位时间步跨越的网格数变多，Courant 条件被破坏。解决：正确做法是网格加密的同时按比例缩小初始时间步，再跑一遍。这个属于网格与时间步的联合收敛，不能只加密网格不缩步长。

5.5 同时修改多个全局参数后结果“跳变”

现象：一次性把孔隙度、渗透率、杨氏模量、注入温度全部改掉，输出结果与默认工况比不是渐进偏移，而是数量级错乱。原因：参数之间存在隐式耦合，孔隙度变化影响渗透率经验公式，杨氏模量变化影响热应力项，同时改会让中间物理量短暂越过合理范围。解决：逐项修改并运行观察，确认每组参数变化都在合理区间后再做多参数扫描。这个习惯后来帮我省了大量复查时间，参数一多，黑匣子就越深，一次动太多变量很难定位是谁把结果带偏了。

6. 结果验证：网格无关性检验与解析解对照的方法

拿到任何 COMSOL 模型，我第一件事永远是验证，而不是信它的默认输出。这套 THM 模型也不例外。验证分两步走：网格无关性检验和解析解对照。网格无关性以注入点压力为目标，分别用基础网格、半尺寸网格和四分之一尺寸网格跑同一算例，比较注入点压力曲线差异。工程经验是三套网格在同一时刻的最大压力差小于 2% 即认为网格密度足够，这个判据比“看起来差不多”靠谱得多。

解析解对照可以用单相压力扩散的 Theis 解做一个近似检验。把模型关掉两相和温度场，只保留单相流动和固体力学，给定恒定注入流量，取出压力降落曲线，和 Theis 公式计算的理论曲线放一起比较。误差来源一般在边界反射，模拟区域尺寸不够大时压力波传到边界再反射回来，后期曲线会明显偏离解析解。这时就要扩大模型区域半径，直到目标时间窗口内压力波未到达边界。这一招对任何渗流模型的修正都很关键，它能告诉你模型在时间维度上到底能跑多远不出边界效应。我在验证这一套模型时就是靠这个流程确认了默认尺寸的设置合理，后面换实际工区参数时也更有底气。从那以后我每次跑新算例前，都强制自己先做一遍网格无关性和解析解对照，不跳过验证直接信结果。希望这个方法对你也能有参考价值，能帮你少走我当年走过的弯路。

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