简介面向石油开采、地热开发、二氧化碳地质封存等工程场景这份COMSOL两相流THM热-流-固耦合模型资源适用于需要开展多物理场数值模拟的研究人员与工程师。模型支持饱和水-油两相流动同时可延伸至水与超临界二氧化碳两种流体的耦合工况热流固三场框架允许依据实际需求修改物性参数、边界条件与网格设置具备收敛性能好、计算效率高、鲁棒性强的优势能够明显降低复杂耦合问题的建模与调试成本。压缩包共14个文件以8份docx技术文档为核心辅以4张jpg示意图与2个html版本页面文档内容覆盖模型原理、控制方程、耦合策略与结果分析图片直观展示模型几何、网格划分及关键场分布整体包体约3.91MB。目前已吸引133人学习使用。读者可借助完整模型文件与技术文档快速复现算例深入理解两相流体流动与热流固耦合的实现过程并在现有基础上自由扩展至不同流体组合或实际工况是进行THM耦合仿真研究的可靠参考资料。1. 两相流THM热流固耦合这套模型到底在解决什么在饱和水-油两相流动的储层里真正难算的不只是某个单一过程而是流动、传热、变形同时发生并且互相反馈注入热水让温度场变化温度变化改变油水的黏度和岩石的热应力孔隙压力变化改变有效应力岩石变形反过来又改写孔隙度与渗透率渗透率一变饱和度前缘和压力分布全部重算。这就是标题里“两相流THM热流固耦合模型”在做的事——把饱和水和油两相达西流动、固体传热、孔隙弹性力学放进同一个方程组求解目标不是让每个场单独算准而是让场与场之间的反馈路径在数值上收敛得住。它凭什么能做到“收敛性好、计算速度快、鲁棒性好”一句话物理建模时砍掉多余自由度初始化与载荷加载给足保护求解器用全耦合阻尼牛顿而非盲目的分离迭代。适合谁做稠油热采、地热储层改造、二氧化碳驱替或封存前期评价的工程师这套建模思路都值得完整过一遍。2. 三场耦合怎么搭THM控制方程与每条耦合路径2.1 饱和水-油两相流动饱和度方程、达西速度与毛管压力两相流部分的起点是水和油各自的质量守恒方程假设岩石骨架变形不改变流体质量源项则水相和油相分别写为水相∂/∂t(φρwSw) ∇·(ρw vw) qw油相∂/∂t(φρoSo) ∇·(ρo vo) qo饱和度约束为 Sw So 1。这里的主变量取水相压力 pw 和含水饱和度 Sw油相压力由毛管压力关系给出po pw pc(Sw)。这种取变量的方式不是随便定的它直接决定了 COMSOL 里因变量列表的写法也决定了边界条件好不好给。两相的速度各自服从达西定律但每相都要乘上自己的相对渗透率vw -(K k_rw(Sw))/μw(T) (∇pw - ρw g)vo -(K k_ro(Sw))/μo(T) (∇po - ρo g)注意 k_rw 和 k_ro 都是 Sw 的非线性函数最常用的是 Corey 形式k_rw S^2k_ro (1-S)^2但这只是最粗的近似真实岩心实验出来的曲线往往带残相饱和度截断残油饱和度高时 k_ro 在 S 接近0 之前就已经归零这种截断直接影响压力和饱和度前缘。温度在这组方程里不是直接出现的而是藏在 μw(T) 和 μo(T) 里。水和油的黏度随温度都可能下降一个数量级以上温度升高、黏度下降、流速上升这就是热场对流动场最明显的驱动路径。做稠油热采时这条路径往往比压力梯度更重要。2.2 传热场能量守恒与对流-导热竞争传热场采用局部热平衡假设即流体与岩石骨架在同一位置共享同一温度这在油藏尺度的 THM 模型里是公认可用的近似。能量守恒方程写为((ρC)_eff) ∂T/∂t ∇·(-k_eff ∇T ρw Cw vw T ρo Co vo T) Q这里的等效热容和等效导热系数按体积平均(ρC)_eff φ(Sw ρw Cw So ρo Co) (1-φ)ρs Cs关键点在于对流通度项油和水各自以自身的达西速度搬运热量不能简单用总流速乘平均热容代替。如果模型里用了一个总的“流速矢量”做对流两相各自的热容差异会被抹掉热前缘的位置会偏移进而影响热应力分布。从数值角度看热前缘和压力前缘的速度往往相差一到两个数量级这给时间步长选择造成了很大麻烦。压力波传播快、热前缘传播慢如果统一用热前缘的时间尺度定步长压力方程会白白吃掉大量计算量反过来用压力波的时间尺度定步长热前缘在模拟周期里几乎不动。后面求解器章节说的 BDF 格式和分阶段加载很大程度上就是为了同时容纳这两个不同尺度。2.3 固体力学场有效应力、热应变与孔隙弹性力学场的控制方程是准静态应力平衡不考虑惯性项因为储层变形的时间尺度远大于波传播时间∇·σ F 0增量形式的本构关系写成dσ D : (dε - αT dT·I) - αb dp̄ ·I其中 p̄ 是等效孔隙压力。两相同时存在时工程上通行的取法是按饱和度加权p̄ Sw pw So po。这个加权对结果影响很大如果只取水相压力高含油区域的骨架载荷会被低估如果只取油相压力含水区域的载荷会被高估。更严格的做法是把每相压力分别作为各自占据孔隙体积的载荷但线性加权已经能满足绝大多数热采和驱替评估。热应变项 αT·ΔT 直接写在力学接口的“热膨胀”设置里。这里的参考温度是储层初始温度而不是地表温度或者室温。对深部储层来说初始温度可能超过一百摄氏度一旦参考温度写错整个初始应力场就会平白多出一个巨大的热应力分量这在后面避坑章里是最容易翻车的点之一。边界条件方面底边固定位移侧边界给法向位移约束辊支撑顶面可以是自由面或按上覆岩层压力加载。特别注意井壁周围流体注入或采出会在这里产生压力集中而井径又很小网格必须局部加密否则力学响应会严重失真。2.4 耦合路径与闭合关系渗透率怎么随变形变化把三套方程放在一起耦合路径可以列成五条流到热通过达西速度的对流换热热到流通过黏度和密度的温度依赖流到力通过等效孔压载荷力到热通过热膨胀反过来影响孔隙度和热容力到流则通过孔隙度和渗透率的更新。前四条在 COMSOL 里都是通过“因变量之间的直接表达式”实现的真正难处理的是最后一条。最常见的闭合关系是幂律型K/K0 (φ/φ0)^nn 在 Kozeny-Carman 推导里取 3很多储层模型也用 3。这条关系式物理上说得通数值上却埋了一个正反馈隐患孔压上升→有效应力下降→骨架膨胀→孔隙度上升→渗透率上升→流体更容易流到该区域→孔压继续上升。反过来也一样孔压下降→骨架压实→渗透率下降→孔压更难消散→继续下降。后者正是压裂或衰竭开发中常见的发散路径。所谓“鲁棒性好”本质上是把这些正反馈闭环加了安全边界。工程上的做法是给渗透率设置下限或者在指数关系外面再用平滑函数截断而不是写硬 if 判断。硬 if 会让雅可比矩阵出现不连续点牛顿法在这种点附近经常直接失步。这里先立住结论后面建模章节会给出具体的表达式写法。3. 在 COMSOL 里落地两相流 THM接口选择、参数清单与耦合表达式3.1 物理场接口选择的两种方案对比如果许可证里有多孔介质相关的模块最省事的方案是直接用官方的两相达西物理场接口再叠加固体传热和固体力学。这个接口默认的因变量就是水相压力和含水饱和度自带迎风机制和相对渗透率模型选择质量守恒性质比自建 PDE 好一个量级。另一个方案是只用一个通用的达西接口加两个系数型 PDE自己写质量守恒方程和所有非线项。这个方案的优势是完全透明缺点是数值稳定项全部要自己控制而且饱和度方程是双曲型不做额外稳定化的话细网格下数值振荡能把求解器拖死。两种方案没有绝对好坏取决于你要拿模型做什么。如果目标是快速给一个热采区块的热前缘和应力场做敏感性分析用官方两相达西接口省下的时间是实打实的如果目标是研究自己实验室测出来的相对渗透率曲线怎么影响缝网扩展那自建 PDE 才能把每个数值细节握在手里。我的习惯是先用方案 A 跑通一个粗算版本确认物理趋势合理之后再决定要不要用方案 B 重写。3.2 参数清单油水两相与岩石热力学参数怎么给参数清单是模型里最容易抄作业的部分也是最容易抄出问题的地方。下面这组参数是我在热采类 THM 模型里常用的基准值单位完全是 COMSOL 默认的 SI 风格直接照搬可以跑通但务必按你的岩心数据替换参数名基准值单位物理含义phi00.251参考孔隙度k01e-13m²参考渗透率约 100 mDSw00.301初始含水饱和度T060degC储层初始温度muW06e-4Pa·s水相参考黏度muO02e-3Pa·s油相参考黏度c_perm81渗透率耦合指数nw21水相 Corey 指数no21油相 Corey 指数E5e9Pa岩石弹性模量nu0.251泊松比alpha_biot0.81Biot 系数alphaT1e-51/K线膨胀系数rho_s2600kg/m³骨架密度Cp_s900J/(kg·K)骨架比热容注意几个容易出错的参数。E 和 nu 应该取排水条件下的岩石骨架模量而不是岩块模量如果用了实验室无水条件下的杨氏模量又同时给 Biot 系数取了 0.8力学响应会偏硬。alphaT 对砂岩通常取 (0.5~1.5)e-5 这个量级但实际值对温度场的应力计算结果非常敏感有条件的话用变温三轴实验数据而不是查手册默认值。温度和孔隙度都会反过来影响这些参数的取值。严格做法是把有效应力、温度都设为参数的函数但第一版模型建议把 E、nu、alphaT 当常数处理先看整体响应不要第一轮就把所有物性都做成非线性否则排查收敛问题时将分不清是哪个非线性项在捣乱。3.3 相对渗透率与毛管压力的表达式写法下面这组表达式是 COMSOL“定义-变量”面板里直接输入的内容变量名以官方接口实际生成为准我这里用 S 表示含水饱和度T 表示温度phi 表示孔隙度k0 和 phi0 是前述参数。这个写法兼顾了物理合理性和数值平滑性// 水相相对渗透率Corey 形式 k_rw S^2 // 油相相对渗透率残油饱和度影响截断 k_ro (1-S)^2 // 毛管压力第一版先给 0稳定后再加 pc 0 // 渗透率随孔隙度变化用上限截断防正反馈发散 Kx k0 * min(exp(c_perm*(phi/phi0-1)), 10) Ky Kx Kz Kx/3 // 油相黏度随温度指数衰减 muO muO0 * exp(-alpha_mu*(T-T0))直接看可能觉得 k_ro (1-S)^2 太粗糙确实如此但第一版模型的价值不是精度而是先让耦合方程组在数值上稳定。残油饱和度、残余水饱和度、端点相对渗透率这些细节应该在稳定模型之后逐个加进去每一个非线性项的加入都对应一次重新收敛测试。Kx 表达式里 max/min 直接用了 COMSOL 的内置函数用来限制渗透率不超过原始值的 10 倍。这个 10 倍上限是经验值压缩导致的渗透率下降不需要上限因为物理上 K 不会变负但膨胀方向如果没有约束正反馈会让 K 无上限增长最终把压力场算出一个不合理的尖峰。类似的思路也可以用在毛管压力曲线上——如果加 van Genuchten 形式注意让它在饱和度趋近端点时保持有限斜率否则雅可比矩阵会出现很大跳变。4. 收敛性与计算速度从哪来初始值、阻尼牛顿与逐步加载4.1 一致性初始化和“先单场后三场”的加载次序一个宣称“收敛性好”的模型最不显眼但最关键的设置往往在初始值这一步。两相流 THM 的初始条件不是简单填两个常数初始水饱和度分布应该由初始油水界面的毛管压力平衡给出初始孔隙压力分布应该满足静水压力梯度初始温度场应该满足地温梯度而初始应力场还要与上覆岩层和侧向应力匹配。如果这些初值只是各自填一个常数求解器第一步就要在同一套方程组里自己把水力、热力、力学三个场“拉”到一致这一步往往直接让牛顿迭代发散。我常用的做法是建一个只包含渗流和传热的稳态预求解模型把因此得到的压力和温度分布作为瞬态 THM 模型的初始值力学场则单独用“初始应力和应变”节点写入地应力位移初值保持全零。这类“分阶段加载”不是偷懒而是把强非线性问题拆成一系列线性程度逐步提高的子问题每一步都从前一步已经收敛的解出发而不是让求解器从零猜起。4.2 全耦合牛顿法与时间步控制参数求解器设置直接决定那句“收敛性好、计算速度快”能不能兑现。我的推荐配置是全耦合牛顿法加 PARDISO 直接求解器雅可比更新频率设为每次牛顿迭代都更新。下面是常用的一组参数可作为起步模板设置项推荐值说明求解器PARDISO强多场耦合下比迭代法稳妥雅可比更新每个牛顿步保证收敛速率不退化阻尼因子0.75振荡时降到 0.5最大迭代次数25超过直接报错不拖时间初始步长特征时间的 1e-4宁小勿大最大步长总模拟时长的 1/100限制跳跃BDF 阶数2兼顾时间精度与稳定性相对容差1e-2 起步1e-4 终算先快后准阻尼因子是牛顿法在非线性区存活的关键。全耦合雅可比矩阵中跨场耦合项的数值可能比对角项大直接牛顿步会过冲阻尼因子把这步收缩到原来七成五代价是迭代次数略增收益是几乎不出现“第一步就发散”这种事故。如果日志里看到残差先下降后反弹把阻尼因子降到 0.5 再试比换稳定化方法快得多。时间步方面BDF 是隐式格式对刚度问题天然免疫但固定步长会浪费大量算力。COMSOL 会根据局部截断误差自动调整步长因此不要手动把时间步设死只限制最大步长即可。若日志显示时间步反复被强制缩短到初始步长附近说明某个非线性项产生了高频振荡优先检查相对渗透率曲线和毛管压力曲线而不是继续折腾求解器。4.3 网格和稳定化两相饱和度方程的数值扩散两相饱和度方程在数学性质上与对流方程更接近流动前缘的后方会出现饱和度间断。这种间断在网格上表现为两种数值现象细网格无稳定化时产生前缘振荡粗网格带迎风时产生前缘抹平。官方两相达西接口自带迎风机制一般不会出现严重振荡但如果你用自建 PDE 版本就必须手动把流线扩散和交叉扩散两个稳定化开关打开。网格尺度和时间步长之间有一个实际的耦合约束每个时间步内饱和度前缘移动的距离不应超过一个网格单元。如果网格太粗而时间步太大前缘会一步跨越多个单元饱和度出现大量负值或大于 1 的异常值。粗网格下这部分误差会被数值扩散掩盖看起来“稳定”了但压力响应和后续的力学响应都会失真。所以验证网格无关性时除了看压力曲线和温度曲线务必把含水饱和度的剖面也画出来检查前缘是否随网格加密产生明显位移。5. THM 两相流建模的五个常见坑从压力振荡到质量不守恒5.1 相对渗透率函数过陡压力振荡与时间步崩盘现象是模型在早期阶段跑得好好的一旦饱和度前缘进入某一区域求解器日志里牛顿迭代次数突然从三四次跳到二十次以上同时时间步被压到初始值的万分之一最终报“找不到收敛解”。原因是相对渗透率曲线的斜率太大。Corey 指数取 3 或 4 时k_rw 在 S0.5 附近的斜率已经非常陡饱和度只要变化一点点两相流速差异就急剧放大雅可比矩阵里这些大梯度元素让牛顿步在迭代空间里反复横跳阻尼因子收缩也救不回来。解决方法是先降低 Corey 指数到 2 或者更平缓的端点截断形式把模型跑稳定再逐步提高指数的过程其实就是为这个参数做敏感性分析。另一个有效的调整是给毛管压力设一个非零但不大的值pc 的引入能在物理上平滑饱和度间断避免纯双曲型方程带来的尖锐前缘。5.2 参考温度给错未加载就先“热胀”了一大截现象是瞬态模型第一步的结果里位移就已经达到上万个微应变孔压场同时出现一个与注入无关的异常分布整个应力云图看起来像被整体加热过。原因几乎永远是固体力学接口里的热膨胀参考温度与固体传热的初始温度不一致。储层初始温度 60 度力学参考温度却用了默认的 20 度方程在 t0 就自动给整个模型加了 40 度的热膨胀荷载初始位移自然异常。解决方法是把力学的热膨胀参考温度显式设为 T0并且在固体传热接口中设置初始值 T0。排查时可以直接画 t0 时刻的应力分布如果初场出现明显的量级异常九成是参考温度没配对。5.3 井点源和渗透率正反馈质量亏与发散同时出现现象是用点源表示注水井时井筒附近的压力先出现一个难以收敛的尖峰然后整个压力场像被“顶住”了一样停滞同时监测面上的累计流量与注入流量对不上账质量守恒误差超过 5%。原因是点源的奇异性和渗透率正反馈叠加在同一个位置。点源处的压力在数学上趋于无穷网格加密只是让奇异的极限更逼近与此同时注入引起的孔隙度增大又抬高了局部渗透率更多的流体涌向井点形成正反馈闭环。解决方法是把井处理成一小段圆柱面边界施加法向通量或压力边界条件而不是几何上的点。渗透率表达式加下限保护例如把 Kx 的 min 上限从 10 改成 3或者整个换成对数关系的 K k0 * exp(c_perm*ln(phi/phi0))正反馈的放大效应会温和很多。质量守恒的核对要等井边界改建完成后再做点源版本的 5% 误差是网格奇异性导致的不是你方程写错了。5.4 力学约束不够刚体位移导致矩阵奇异现象是求解器报矩阵奇异而且报错位置在力学接口的某一行但压力场、温度场看起来都正常有时干脆给出一堆时间为 nan 的结果云图空白。原因是力学场存在刚体自由度。储层模型如果顶部全是自由面底部又没有固定约束整个实体可以整体滑动或旋转刚度矩阵奇异牛顿法自然失效。这个坑在二维模型里尤其常见二维只约束水平位移和垂直位移还不够还要约束绕面外轴的转动自由度。解决方法是底部给固定约束侧边界给法向位移为零的辊支撑必要时在模型角落或对称轴上再加弱约束把最后几个零特征值消除。加完约束后顺手算一下特征值确认最小特征值不在机器精度量级这一步比肉眼检查边界条件可靠得多。5.5 排查习惯每加一个物理场都应该有一个“对照组”现象是模型昨天还能收敛今天加了一个新的物理场或者改了一条材料曲线之后突然全面发散而且根本分不清是哪个场出了问题只能一版一版回退。原因是缺少可对照的中间版本。全耦合 THM 模型里任何一个环节的非线性增强都可能引发发散没有对照版本排错就是在黑匣子里猜。我现在的做法是严格按照“等温单相流动 → 两相流动 → 开启传热 → 最后开启力学”四个阶段保存版本每一步都确认残差曲线和关键物理量没有量级跳变再进入下一步。哪个版本开始发散发散就出在这个版本新增的环节里。这个习惯本身比任何一条具体参数改进都值钱能省掉大半排查时间。6. 从油水版本扩展到水超临界 CO2 两相流你要改的四个关键点6.1 物性从常数到温压插值把油相换成超临界 CO2最大的结构性变化是物性从“一组常数”变成“随温度和压力剧烈变化的插值函数”。CO2 在临界点附近密度可以随压力变化好几个量级必须用插值表或公开状态方程导出的数据在 COMSOL 里以插值函数接入。注意插值数据的温压网格要有足够密度尤其跨临界点附近否则一阶导数不连续牛顿法会在临界区反复失步。6.2 相平衡改法先忽略溶解再补亨利项水和 CO2 并非完全不互溶但第一版模型强烈建议先按完全不互溶处理把两相流动结构原样保留只换物性。等整个 THM 耦合在 CO2 物性下重新收敛、验证通过之后再加入亨利溶解项把它作为能量方程和质量方程的附加源项。一上来就加溶解会让方程组额外增加两个未知数常见做法是先“冻结”它拿到稳定基线再说。6.3 验证入口等温注气与饱和度前缘对照验证 CO2 版本我会先做一个等温、单相、恒压注气的小模型核对密度场和压力场与状态方程计算值一致再做等温两相注气把含水饱和度剖面与一维解析驱替解对照确认前缘速度和形状合理最后才把这些子模型合并进完整的 THM 模型。这样把“物性错误”和“耦合错误”分开排查最大程度避免两者混在一起时无解可查。把油换成 CO2 的坑主要藏在物性导数的不光滑里所以“另议”从来不是几十行的改动而是一个要求重新调参的迭代过程。我自己接手这类改造时至少预留两轮完整调优周期第一轮只求不发散第二轮才追求结果精度。希望这篇整理能帮你在遇到两相流 THM 模型时少走弯路直接把精力花在参数校验和物理分析上。本文还有配套的精品资源点击获取