两相流THM热流固耦合建模实践：COMSOL实现与收敛性避坑

简介COMSOL两相流THM热流固耦合模型面向从事石油工程、地热开发、二氧化碳封存等领域的研究者与工程师可模拟饱和水和油两相流动也可扩展至水超临界二氧化碳两相流情况。模型实现热-流-固三场耦合能够在考虑温度、流体压力和固体变形相互作用的前提下分析多孔介质中的渗流与传热过程。压缩包共14个文件大小3.91MB其中8个Word文档提供模型理论基础、控制方程、参数设置及结果解读的详细说明4张JPEG图片展示几何建模、网格划分和温度/压力分布等典型结果2个HTML文件便于浏览器直接查看模型思路概览。已有130人学习使用说明其具备一定参考价值。模型收敛性好、计算速度快、鲁棒性强且结构清晰用户可根据研究需求修改流体参数、边界条件和耦合项适配不同场景便于快速开展两相流与热流固耦合仿真。1. 两相流THM热流固耦合模型为什么值得你花时间搭建两相流THM热流固耦合模型听起来是一套层层嵌套的黑匣子又要算饱和水和油的流动又要算温度场还要算骨架变形。我这几年的项目验证里有个反直觉的经验三场耦合不一定比两场更难收敛真正让模型翻车的往往是变量量纲和离散顺序没对齐。这个方向的典型价值藏在工程预判里——注入井附近温度、压力、应力同步演化先看趋势再定方案比一次性堆满细节可靠得多。适合正在做地热回灌、稠油热采、CO2封存方案预判的工程师也适合做热-流-固耦合研究的学生作为起步基线。读完后你能判断它适不适合你的问题并能照着一套可复现的路径搭出模型。2. 把THM三场耦合拆开看控制方程、耦合路径与选型理由2.1 三场各自管什么耦合点在哪几个方程里热流固耦合里三个场负责的事边界其实很清楚。流动场用达西定律描述多孔介质中的渗流饱和水和油各自有压力、饱和度和速度温度场通过能量守恒方程把热传导和热对流放在一起固体场则是孔隙弹性介质在外力、热应变和孔隙压力作用下的变形。关键不是单独解这三个方程而是它们互相咬合的几个点温度改变流体粘度和密度粘度密度改变流动速度流动速度再反过来改变热量输运孔隙压力和热应变改变有效应力有效应力又改变骨架变形而骨架变形会挤压缩小孔隙导致孔隙度和渗透率实时变化。以两相达西为例一个常见形式是每个相写质量守恒∇·(ρᵢ λᵢ K (∇pᵢ - ρᵢ g)) qᵢ ∂(φ Sᵢ ρᵢ)/∂t其中 λᵢ 是相对渗透率相关的流度。能量方程里最容易被漏掉的是对流项里的流体速度应该用达西速度而不是实际孔隙速度漏了会把热前锋算晚。固体力学方程用有效应力形式σ_eff σ_total - α_B pα_B 是 Biot 系数一般取 0.7-1.0。耦合路径可以用一句话记温度影响流动是通过物性压力影响力学是通过有效应力变形影响流动是通过孔隙度和渗透率。这三条链路就是多物理场节点里要挂的东西。很多从单一力学转来做多场的人最容易犯的错是只做了单向耦合比如先算流动拿到压力再把这个压力当不变的载荷去算力学。这在稳态弱耦合问题里偶尔能凑合但只要注入温度变化明显单向耦合的前提就不成立。2.2 为什么选 COMSOL 而不是手写有限元代码这里不是替商业软件说话而是从工程账上来算两相流THM模型要在代码层面同时协调非线性材料、移动的饱和度前锋、随时变化的渗透率场和解耦困难的压力-变形关系如果是自己写有限元代码调试周期长得离谱。常见做法是用 COMSOL 这类带多物理场装配的有限元环境把每个场的方程以接口形式摆出来多物理场节点自动生成耦合项不用手工组装刚度矩阵和非线性雅可比。COMSOL 模型优势中最实际的一点是求解器栈完整。非线性牛顿迭代带阻尼和线搜索分离式求解器可以把压力和饱和度放一组、温度放一组、力学放一组按物理强耦合程度分组求解这对条件数很高的多孔介质问题非常友好。如果换成自编程序光处理压力-变形交叉项的前置条件和预处理就要花掉一到两周。当然 COMSOL 也有它的边界它对网格数量的扩展性一般超大规模三维模型跑起来内存占用明显高而且求解器里那些看似自动的默认选项比如阻尼因子初始值、非线性残差收敛阈值在强非线性两相流问题里经常默认值过于乐观。我的经验是不要全信默认设置打开求解器细节把阻尼上限从 1.0 降到 0.8 或者 0.5收敛性立刻稳下来。这一点留到第五章展开。2.3 模型边界和材料参数怎么定量纲最容易在这里翻车材料参数和边界条件直接决定两相流THM模型收敛难度和物理正确性。我一般先定一张参数表把所有量统一到 SI 制再去考虑单位换算。下面这张表是从热采场景里抽出来的初值模板实际值按你的储层换但量级可以参考。参数示例值作用与说明初始孔隙度 φ₀0.25控制储集能力也参与有效应力修正初始渗透率 K⁰1e-14 m²流动能力基线约等于 1 md 量级水相动力粘度 μ_w0.001 Pa·s随温度按指数或 Arrhenius 关系修正油相动力粘度 μ_o0.02 Pa·s热采里随温度下降通常能降一个量级岩石密度 ρ_s2600 kg/m³用于力学重力载荷岩石热导率 k_s2.5 W/(m·K)热传导主导项岩石比热容 C_s900 J/(kg·K)能量方程储热项水比热容 C_w4180 J/(kg·K)对流热载能力油比热容 C_o2200 J/(kg·K)对流热载能力Biot 系数 α_B0.8有效应力计算关键系数入口流速0.01 m/s用达西速度给不是孔隙流速出口压力1 atm参考压力边界反射风险低单位换算在这里要特别小心。渗透率常用毫达西 md 表示1 md 约等于 1e-15 m²如果直接导入数值不换算就会出现量级离谱的问题。压力也一样COMSOL 默认压强单位是 Pa你习惯用 bar 或 MPa 时可以在参数里写成10[MPa]求解器会按 COMSOL 语法自动换算但如果你在表达式里混写10e6和1e7最后排查起来非常折腾。我通常用带方括号单位的形式这是 COMSOL 参数表达式里非常实用的做法。边界条件上入口不要直接给压力脉冲最好给流量或流速边界让压力和饱和度自己建立起来这样初期的非线性冲击小得多。出口给定压边界模拟外边界连通如果模型四周全封闭压力会随时间持续升高物理上没错但数值上容易出现压力漂移。另外初始压力一定要和流体密度设置自洽否则一开算就产生人为流动。3. 从零搭一个饱和水油两相流THM模型接口、参数与求解器3.1 几何、网格与初始值粗网格反而更好收敛一个典型起步模型可以简化成二维长方形域比如长度 50 m、高度 20 m。左侧代表注入井附近右侧代表生产端或者开放边界中间是储层。先别着急按真实井网建复杂几何二维单注入单采出的概念模型足以验证耦合机制和参数敏感性。网格策略上我常用的做法是整域最大单元尺寸控制在 1-2 m在入口和出口附近做一层 0.3-0.5 m 的局部加密。重点提醒饱和度锋面区域不需要刻意加密尤其在隐式时间步下过度加密会让局部库朗数变大非线性迭代反而更容易震荡。这个坑我踩得很深具体现象写在了第五章。初始值设置要尽量靠近稳态。如果初始孔隙压力给 10 MPa而模型内部还有个压力突变则第一步迭代就要消化这个大梯度过大的残差很容易发散。我的习惯是用一个非常小的初始时间步比如 1e-2 秒先让压力场扩散再做正式时间推进但更常用的做法是直接用 COMSOL 的辅助扫掠功能从静水压力解起步相当于先用小步长瞬态算一个初值。这里给一个几何和网格的基本框架用 COMSOL 的 Java API 示意模型定义过程实际用图形界面等价操作// 设置基础参数 model.param().set(L, 50[m]); model.param().set(H, 20[m]); model.param().set(phi0, 0.25); model.param().set(k0, 1e-14[m^2]); model.param().set(mu_w, 0.001[Pa*s]); model.param().set(mu_o, 0.02[Pa*s]); // 建立二维几何 model.component().create(comp1, true); model.geom().create(geom1, 2); model.geom(geom1).feature().create(r1, Rectangle); model.geom(geom1).feature(r1).set(size, new double[]{50, 20}); // 划分网格整域最大尺寸 2m入口边界加密到 0.3m model.component(comp1).mesh().create(mesh1); model.component(comp1).mesh(mesh1).autoMeshSize(2);这段代码最关键的是autoMeshSize(2)它会按比例给全局网格一个参考尺寸。实际操作中我会把这个数值调大而不是调小经验是两相流THM模型网格太细会让饱和度前锋处产生压力小尖峰进而拖慢迭代先把粗网格跑通再局部加密而不是一开始就把网格拉满。3.2 物理场接口与耦合变量把渗透率随有效应力变化写进去用图形界面建模型时物理场接口按顺序添加固体力学、传热、达西定律如果你的 License 里有地下水流模块可以直接选两相流达西没有也别慌用饱和度和压力两个因变量的系数型 PDE 也能搭。我这里为了说清楚耦合变量写法用 Java API 展示接口初始化顺序// 创建三个物理场接口 model.physics().create(solid, SolidMechanics, geom1); model.physics().create(ht, HeatTransfer, geom1); model.physics().create(darc, DarcyFlow, geom1); // 定义多物理场耦合变量渗透率随有效应力变化 model.variable().create(var1); model.variable(var1).set(k_eff, k0*exp(-3*(sigma_eff/1[MPa]))); model.variable(var1).set(mu_eff, (S_w*mu_w (1-S_w)*mu_o));这段代码里sigma_eff是固体力学模块算出的有效应力除以 1 MPa 是为了归一化前面的系数 3 控制了应力敏感度。这个指数型渗透率衰减模型在孔隙弹性里很常见工程上通常取 3-5 之间配合少看绝对渗透率、多看归一化趋势来标定。mu_eff用了饱和度加权平均是一种平滑近似实际项目里用对数平均更接近相对渗透率曲线但初版用加权平均更容易收敛。饱和度变量的处理是两相流模型的关键。如果用的是两相流达西接口饱和度方程已经封装好不用额外写如果是用系数型 PDE 自己搭那要对应一个因变量S_wmodel.physics().create(c, CoefficientPDE, geom1); // 系数型 PDE 设置d 是阻尼系数gamma 是对流通量 model.physics(c).feature(cfe1).set(d, phi0*rho_w); model.physics(c).feature(cfe1).set(gamma, rho_w*k_eff*lambda_w*(grad(p)-rho_w*g)); model.physics(c).feature(cfe1).set(f, 0);这里的lambda_w是水相的流度包含相对渗透率除以粘度的计算。系数型 PDE 写法很自由代价是把控要求高一旦系数写错很难从报错信息里找到方向。我的建议优先用现成的两相流达西接口自建 PDE 只用于补充附加组分方程否则调试成本成倍上升。3.3 求解器配置分离式与全耦合怎么选求解器配置直接决定这个模型是十几分钟跑完还是跑一夜还发散。常见做法是对两相流THM模型优先选择分离式求解器把求解分成两步或三步——第 1 组算压力和饱和度两相流强耦合第 2 组算温度第 3 组算固体力学。每组内部用阻尼牛顿迭代组间用缓和的顺序耦合。为什么这么分压力和饱和度之间通过相对渗透率强耦合必须放一组里同时迭代温度对流动的影响通过物性间接传递可以放到下一轮力学场与流动场之间有 Bohm 效应但成熟的迭代策略下往主循环里塞全耦合牛顿反而容易让雅可比矩阵过于稠密内存和时间双双爆炸。只有当你的模型是弱非线性、几乎没有饱和度前锋的时候全耦合才划算。瞬态求解时间步的设置也应该分组。初版不要直接设固定大时间步我用的是自适应时间步加一个合理的初始步长比如model.study().create(std1); model.study(std1).feature().create(time, Transient); model.study(std1).feature(time).set(tlist, range(0, 3600, 60*24*30)); model.sol().create(sol1); model.sol(sol1).study(std1); model.sol(sol1).feature().create(st1, Time); model.sol(sol1).feature(st1).set(tsteps, strict); model.sol(sol1).feature().create(seg1, Segregated);tlist里的写法是从 0 秒开始每 3600 秒输出一个点总时长约 30 天。tsteps设成strict意味着严格按输出步长推进如果你的模型在前峰阶段容易发散可以改为不严格模式让求解器在内部自动加密时间步。分离式求解器在内存上比全耦合省得多对个人电脑上的三维模型尤其明显。3.4 后处理怎么证明模型物理上没算错求解完成后第一件事不是看云图好不好看而是做四个物理合理性检查。第一个是饱和度越界检查用全域探针统计 S_w 的最小和最大值只要超出 [0,1] 范围说明离散或时间步设置出问题。第二个是压力连续性在任意剖面上看压力是否平滑过渡如果有孤立跳变点先怀疑网格质量。第三个是能量守恒在全局计算界面分别积分入口和出口的能量通量再对比域内总能变化差值应该收敛到时间步容差量级。这个检查能区分“温度场在缓慢变化是正常的物理现象”和“数值泄漏”避免你把收敛误判成发散。第四个是探针回看在注入点附近放一个压力探针如果压力响应曲线出现锯齿说明局部时间步精度不够需要调小时间步容差或者改阻尼因子。一个实用小技巧以 30 天为总时长时先把总时长缩短到 3 天把初始时间步放到 1e-3 天验证 3 天内的趋势是否平滑再放宽到完整时长。这样做能提前暴露出大概率会发散的现象而不是浪费一整夜跑完后发现第 7 天就崩了。4. 从水-油切到水超临界CO2物性、相态与参数修改清单4.1 超临界CO2的物性为什么不能当常数处理标题里提到“水超临界二氧化碳两相流的版本另议”这其实是比水-油更贴近碳封存和驱替场景的改造方向。超临界CO2与油的物性相比有一个巨大差异它的密度不是常数而是在临界点附近随压力和温度剧烈变化。超临界CO2在 8-20 MPa、60-80 °C 之间的密度可能在 200-700 kg/m³ 之间波动粘度也可能上下浮动数倍。如果你直接把油相的物理性质替换成一段常数压力解释和饱和度推进都会失真。比较靠谱的做法是做一张随压力和温度变化的插值表把密度、粘度、比热容、热导率整理成四列数据然后在 COMSOL 里用插值函数加载。数据来源用物性数据手册或公开物性数据库都行关键是覆盖你计算压力范围内的全部区间不能只查一个点。4.2 相态模型互溶与不互溶的取舍水与超临界CO2并不像水-油那样几乎不互溶CO2在水中有一定溶解度反过来水在CO2相中也微量存在。严格处理需要组分输运方程、状态方程和局部平衡假设模型的非线性程度显著提高。标题说“另议”说明默认走简化路线假设两相不互溶、局部平衡饱和度方程和相对渗透率曲线沿用两相流框架。这个简化在初版方案里完全够看因为你要先抓住压力和温度响应的主趋势而不是细扣溶解度带来的质量汇。但要注意简化后不能用于CO2矿化封存相关的组分分析只能用于物理封存里的压力与温度影响评估。如果后面必须考虑互溶我会把溶解量处理成一个与溶解度相关的源项而不是直接改成双组分全耦合方程。在 COMSOL 里物性插值函数的建立方式不复杂图形界面在“全局定义-函数-插值”里选好数据表自变量列表写p和T因变量写rho和mu。Java 片段则是这样的示意model.func().create(table_cp, Interpolation); model.func(table_cp).set(table, new double[][]{{10,323,550,0.00002}, {10,353,400,0.00002}, {15,323,700,0.00003}, {15,353,520,0.00003}, {20,323,780,0.00004}, {20,353,600,0.00004}}); model.func(table_cp).set(argnames, new String[]{p,T}); model.func(table_cp).set(funcs, new String[]{rho,mu});这段插值数据只作为占位示意真正建模时要放一个覆盖完整压力和温度网格的表至少 30-50 个数据点插值才够平滑。插值表选线性插值即可不要用三次样条因为CO2物性在临界点附近本来就有高压梯度样条容易产生伪震荡。4.3 从水-油切换到水-CO2版本的参数修改清单参数水-油模型典型值水-CO2模型建议值修改原因初始压力5 MPa10 MPa需高于CO2临界压力 7.38 MPa初始温度293.15 K323.15 K需高于CO2临界温度 304.13 KCO2密度油相密度 800 kg/m³插值表 400-800 kg/m³随压力温度剧烈变化CO2粘度油相粘度 0.02 Pa·s约 0.00002-0.00005 Pa·s超临界态粘度接近气体量级差 2 位以上CO2比热容油相 2200 J/(kg·K)1000-1200 J/(kg·K)对流热容改变CO2热导率油相 0.13 W/(m·K)0.02-0.05 W/(m·K)热扩散比例变化注入流率0.01 m/s0.005 m/sCO2粘度低流速快过大的注入流率容易数值发散界面张力系数与油体系标定重新标定相对渗透率曲线和毛管压力随界面张力变化切换成水-CO2两相流版本时最容易翻车的是 CO2粘度太小导致流速变大在相同的界面通量下前锋推进速度远快于水-油体系。如果时间步还沿用原来的设置前锋会在几个时间步内越过几个单元尺寸出现时间和空间上的振铃。解决的办法不是把网格加密而是降低达西速度或时间步让每个时间步的前锋推进距离小于一个单元尺寸。这里注意一个单位陷阱水的粘度是 0.001 Pa·s超临界CO2粘度可能是 0.00003 Pa·s数值上差了约 30 倍。这个量级差距会直接反映到达西速度上如果不修改求解器阻尼参数第二步迭代残差就可能放大几十倍然后发散。5. 收敛性与鲁棒性避坑排查五个高频踩坑点5.1 细网格让饱和度前锋周围压力震荡现象模型在饱和度前锋附近出现逐点振荡的压力残差下降曲线像锯齿时间步越缩越小。原因锋面宽度远小于局部网格尺寸时隐式时间积分和中心差分会产生伪振荡过度加密让局部库朗数变大而库朗数条件被打破。解决优先保证前锋区域的单元尺寸与每个时间步内饱和度推进距离接近而不是一味加密。如果必须加密请同时把时间步长对应缩短并检查是否用了迎风格式。一个可复现的做法是把全局最大单元尺寸设为 2 m入口局部加密 0.5 m时间步初始值设为 0.01 天跑通后再逐步加大。5.2 瞬态能算稳态求解器却不收敛现象同一套模型改成稳态研究后报“求解失败”残差不下降甚至上升。原因稳态问题在给定强非线性条件下可能不存在唯一解或者初始猜测太差。两相流THM里的温度场和饱和度推进往往天然是瞬态过程强行稳态会抹掉物理上的延迟响应。解决先用瞬态求解到接近稳定的时间点把结尾时刻的场保存为初始值然后再切到稳态求解器。也可以对入口流量做一个从 0 到目标值的参数扫描让求解器先在一个温和的载荷下建立解再逐步提高流量。5.3 温度场长期慢漂分不清是物理还是数值现象压力饱和度和力学场都看着正常唯独温度场在很长时间内持续缓慢上升或下降残差已经小于容差但解仍在漂移。原因可能来自绝热边界导致热量在域内累积也可能是能量方程里的对流项顺序写错或时间积分容差太松导致累积误差。解决做一个整体能量平衡检查——分别积分所有边界上的热流入和热流出对比域内热储量的变化率。如果相对误差在容差量级比如相对值小于 1e-3这个慢漂是物理的不用管如果误差持续单向增大就把非线性和时间步的容差同时收紧并核查温度边界条件是否存在重复。5.4 参数跨越数量级导致雅可比矩阵病态现象报错里出现“奇异矩阵”或迭代长期卡在同一个残差值减少时间步也没有明显帮助。原因渗透率 1e-14 m²、粘度 1e-3 Pa·s、岩石刚度 1e10 Pa 这几个数值放在一个矩阵里条件数极高。如果单位不统一比如压力用 bar 写成1e5的倍数更容易让矩阵病态。解决全部用 SI 单位表达关键参数所有自定义表达式里的参数都带上单位方括号。然后在求解器配置里打开自动缩放功能并给每组的因变量设置合理的初始值。经验上这一步解决了一半以上的两相流THM收敛困难属于成本最低、收益最高的操作。5.5 饱和度越界S_w 出现大于 1 或小于 0现象计算报告里饱和度场出现绿蓝色斑块探针显示 S_w 超过 [0,1] 范围随后时间步一砍再砍直到求解失败。原因两相流饱和度方程是双曲型方程默认的中心差分离散格式不适合处理强对流前锋数值振荡让饱和度越过物理边界。解决换成迎风格式或给饱和度方程加一个极小的数值扩散项来驯化高阶振荡。另外在变量定义里加一个饱和度的投影限制表达式比如min(1, max(0, S_w))虽然这是一层保险但能让求解器不至于因为越过界限而崩溃。真正的根子还在离散格式和时间步匹配上。6. 让模型跑得更快又稳的进阶技巧6.1 用对称性把模型缩到四分之一两相流THM模型最大的成本在非线性迭代而不是方程数量本身。如果你的几何是规则的矩形或轴对称井网先用对称性把域缩到四分之一或一半网格量直接砍一半以上。边界条件里对称边界用零通量效果和完整域一致。这个技巧能让你在同样的机器上把三维的模型从跑不动变成几分钟一轮。6.2 参数化扫描做初值预热不要直接在最严苛的参数组合下起算。常见做法是先把入口流量降到目标值的十分之一跑一个短时间得到相对平滑的场再把流量作为参数扫描目标值逐步提高每一步以上一步的解作为初值。这比任何阻尼设置都有效因为非线性求解器最怕的不是强非线性而是从一个毫无根据的初值开始硬来。6.3 分离式求解器加阻尼系数协同调分离式求解器不是选了就完事。把“完全耦合”改成“分离式”后还要在每个求解步里打开高级设置看看阻尼类型。我一般会把阻尼上限从 1.0 降到 0.6阻尼因子设为 1.2让每次迭代不要走太远代价是迭代次数增加一点但稳定性大幅提升。配合时间步的严格模式关闭求解器可以在前锋处自己加密时间步整体时间反而比一路跌跌撞撞更快。我印象最深的一次翻车是我为了追求云图精细把油水前锋附近的网格加密到 0.2 m结果跑了一夜都没收敛。后来把网格粗化到 1 m反而半小时跑完饱和度锋面平滑得离谱。那一轮之后我就养成了习惯任何两相流THM模型都先跑一个 1-2 m 量级的粗网格作为预演确认趋势正常再做局部加密。这个习惯帮我躲过了后面好几次时间上的血亏希望帮到你。本文还有配套的精品资源点击获取