两相流THM热流固耦合建模：COMSOL收敛与求解器配置指南

简介基于COMSOL构建的两相流THM热流固耦合模型面向石油工程、地热开发及二氧化碳地质封存等复杂多物理场耦合问题适合从事渗流力学与岩土工程仿真的科研人员、工程师与研究生。模型以饱和水和油两相流动为核心支持温度场—渗流场—应力场三场耦合求解也能按需求扩展至水与超临界二氧化碳体系用户可根据自身研究对象调整材料参数、边界条件与几何模型计算收敛性好、速度快、鲁棒性强。整套资料共14个文件压缩包约3.91MB其中8个docx文档系统讲解模型构建过程、控制方程、参数设置与结果分析思路4个jpg图片直观展示模型几何造型、网格划分及典型场分布结果另有2个html文件便于交互式快速浏览关键要点适合作为THM耦合建模的参考模板与二次开发基础。该资源已有133人浏览学习可帮助读者快速上手两相流热流固耦合仿真的建模与分析。1. 两相流 THM 热流固耦合给地热与油藏模拟省下三个月调试时间做地热回灌或者油藏注采模拟的人对 THM热-流-固耦合应该不陌生温度场改变流体黏度和骨架应力孔隙压力反过来影响有效应力变形又改变孔隙度和渗透率。三场互相咬合最头疼的不是物理模型本身而是让这一组高度非线性的方程组在一个模型里稳定收敛。我以前搭过一套水-油两相 THM 模型光调初始值和时间步就折腾了快两周。这套 COMSOL 两相流 THM 热流固耦合模型解决的就是这个阶段的问题把三场耦合的方程、模块和求解器预配置做到可直接运行的程度拿到的模型能先算通再按自己的工况换材料参数、换边界条件。适合做地热、油藏、CO2 封存类项目的研究生和工程师。2. THM 三场耦合建模控制方程、相态选择与模块组合策略2.1 热-流-固三场耦合的骨架孔隙弹性、热膨胀与达西流THM 耦合的核心是三个场共享同一套几何和材料属性每一步求解都做数据交换。固体力学接口负责位移场孔隙压力 p 通过 Biot 系数对有效应力做折减流体模块计算达西速度速度由绝对渗透率、相对渗透率、黏度和压力梯度共同决定传热接口求解温度 T而 T 既通过热膨胀系数影响应变又改变水和油的密度、黏度反过来影响流速和压力分布。三个场就这样闭环咬合在一起。实际计算中压力、温度、位移建议同时求解也就是全耦合法而不是顺序耦合。顺序耦合虽然每一场单独看都是线性的但两相流带来的强非线性会让场与场之间信息滞后容易出现压力波和饱和度波错位。全耦合的代价是单步计算量更大但稳定性明显更好这也是这套模型标称“收敛性好、计算速度快”的原因之一——求解器预配置到位省去了大量手动切换的麻烦。骨架参数里最要紧的是有效应力法则的取法。常规孔隙弹性模型中有效应力等于总应力减去 Biot 系数乘以孔隙压力。Biot 系数默认 1 可以理解为颗粒不可压缩这对致密储层是合理的但如果是松散砂岩或裂缝介质取 0.7~0.9 更贴近真实。热膨胀系数建议只给骨架本体流体热膨胀通过密度函数体现不要在流体材料节点里再叠一套热膨胀否则温度应力会重复计算。2.2 两相流驱动饱和水和油的相输运模型怎么选这套模型声明支持饱和水和油两相流动。COMSOL 里两相流有两个完全不同的方向一种是自由流动的水平集/相场法用于追踪界面形态另一种是多孔介质里的两相达西流。对储层和地热这类渗流问题正确选择是多孔介质模块里的两相达西流动接口主变量是压力 p 和海相饱和度 S_w油饱和度 S_o 1 - S_w不需要追踪相界面。两相达西流的关键输入是相对渗透率曲线和毛细压力曲线。常见做法是用 Corey 型曲线水相相对渗透率取归一化饱和度的平方油相取一减归一化饱和度的平方。模型文件里一般内置这套曲线但不同岩性的指数差异很大砂岩水相指数 2~3裂缝里可能接近 1。改模型时如果只替换绝对渗透率而不换相对渗透率曲线驱替前沿形态会差很多这点后面避坑章节还会展开。油水密度和黏度不要用固定常数至少要做成温度或压力的插值函数。原油黏度随温度下降可以差几倍到十几倍固定常数会让温度场计算失真等于砍掉了 THM 耦合里最重要的一条路径。模型文件里如果内置的是常数建议自己在材料节点里把“用户定义”改成“插值”读入一张温度-黏度表格计算会慢一点但结果可信度完全不一样。2.3 模块组合策略固体力学、流体流动与传热接口的变量维度模块组合不是把三个物理接口简单堆进模型树就完事要检查每个接口的因变量和耦合节点。标准组合应该是固体力学sol求解位移三个分量两相达西流dsw求解压力和饱和度传热ht求解温度。多物理场节点里会自动生成“热膨胀”“孔隙弹性”“流体热属性”等耦合项。这里有个常见坑固体力学本身不求解温度只是调用传热接口算出的 T。如果固体力学节点里直接填一个温度场表达式而不是引用 T 变量热膨胀子节点就没有被真正激活。判断方法很简单展开多物理场节点看有没有出现“热膨胀”和“热-流”耦合项没有就是没连上。传热接口也建议用“多孔介质传热”它会把骨架和流体的等效热容、等效导热系数一起算进去而不是只算固体的导热。提示正式运算前先确认求解器是瞬态。THM 地热和注采问题本质是瞬态的注入前沿推进、热突破、压力扩散都需要时间。稳态求解器只适合做初始地应力平衡或者用于求初始压力场。3. 模型落地几何、边界条件与求解器参数设置清单3.1 几何与网格影响收敛性的第一道关卡模型文件的几何通常按储层概念模型构建导入后别直接拿默认网格算。先看几何单位COMSOL 里米和毫米最容易出错储层尺度常用几百米到几公里而井筒半径可能只有 0.1 米量级跨尺度问题必须局部细化。网格策略上我一般用自由四面体为主近井区域加局部加密加密半径至少 5~10 倍井径。如果几何有明确分层比如盖层、储层、底板结构用扫掠网格生成六面体收敛速度会明显提升尤其是对饱和度前沿推进这类双曲型问题。注入面和热源面加 2~3 层边界层避免热通量在第一个单元里就被数值抹平。网格数不是越多越好。两相流 THM 模型网格数建议控制在 10 万到 50 万之间超过 100 万对多数机器会显著拖慢。如果几何允许做轴对称直接用二维轴对称模型网格控制在 2~5 万计算速度快一个量级做参数敏感性分析完全够说明问题。先跑通再换三维这是我从这套模型上得到的最直接经验。3.2 边界条件与初始条件注采井、远场与热源的赋值逻辑边界条件赋值要先分清哪些是压力边界、哪些是通量边界。混合搭配是最常见的建模错误下面这张表可以直接照着设位置推荐条件参考值注入井流量边界或压力边界二选一流量按 m³/s 给定或压力高于初始地层压力 2~5 MPa采出井压力边界定压开采井底压力低于初始地层压力远场边界封闭或定压模型范围够大时远场定初始压力顶底边界热通量或绝热浅层可给土地温梯度深层按绝热初始条件静水压力初始温度油水饱和度压力按深度梯度温度按地温梯度注入流量和压力不能同时设成第一类边界COMSOL 会直接报过约束。井筒附近如果用点源加载到某个节点网格粗时会形成奇异压力场建议用一个小圆柱或球体区域代替点源半径 0.5~1 m压力场会平滑很多。初始条件里最容易被忽略的是饱和度初始值。很多模型给全区域 S_w0.2没考虑束缚水饱和度。如果初始饱和度低于束缚水饱和度相对渗透率曲线在水相起点处就是 0体系里实际没有水相流动两相流方程直接退化后面很难纠回来。初始饱和度应该设在束缚水饱和度之上比如 0.25~0.35。3.3 求解器配置直接法、阻尼 Newton 与时间步约束COMSOL 默认求解器对单相流动问题很好用对 THM 却不一定。这套模型预配置好后通常已经选了全耦合。需要留意的参数有线性系统求解器推荐 PARDISO。MUMPS 内存更大但更稳内存吃紧时换 PARDISO。非线性方法自动阻尼 Newton最大迭代次数设 10~15 次。超过 20 次基本说明模型设置有问题不要盲目等它跑完。时间步建议设最大时间步约束注采初期压力波传播快时间步 0.1 天或更小稳定后放宽到 1~5 天。相对容差默认 0.01 够用做定量敏感性分析建议改 0.001计算时间大约增加 30%~50%。求解器参数建议值说明线性系统PARDISO内存占用中等速度快最大非线性迭代10~15超过 20 基本是设置有问题初始阻尼1e-3 ~ 1e-2收敛困难时先降初始阻尼最大时间步注入段 0.1 d之后 5 d防止一步跨过饱和度前沿相对容差0.01高精度需求时改 0.001日志里出现 Inf 或 NaN 时九成是材料参数有除零重点查渗透率、黏度和孔隙度是否在某个区域被赋成 0。盖层区域孔隙度低但别设 0设 0.001~0.01 既符合物理又避开数值爆炸。4. 避坑清单两相流 THM 模拟最常翻车的六个坑两相流 THM 模型翻车概率最高的环节集中在两个阶段求解器阶段和物理量阶段。下面这六个坑是我复盘多次失败后整理出来的每条按现象到原因到解决记录。4.1 求解器阶段初始值、时间步与异常的三种排查现场现象点计算就弹“找不到一致的初始值”求解器日志停在 0 时刻。原因初始条件与边界条件冲突。最典型的是初始压力与注采井压力差几个量级COMSOL 在组装初始刚度矩阵时得不到一致初场。解决做两步走。第一步把两相流接口暂时设为单相只保留压力变量用稳态求解器算出一个平衡压力场第二步把稳态结果作为瞬态的初始值重新启用两相流。如果还过不去把初始阻尼从默认 1 改成 1e-3并打开求解器设置里的自动稳定选项。现象时间步越缩越小日志里布满“所需步长小于最小步长”。原因非线性残差不收敛自适应时间步被一步步逼到下限。两相流最容易触发这个状况的是饱和度前沿推进和毛细压力突变。解决限制最大时间步并把 BDF 阶数上限从 2 降到 1牺牲一点精度换稳定性。同时检查相对渗透率曲线有没有急剧转折我一般把 Corey 指数大于 3 的曲线做平滑过渡处理用分段线性连接计算稳定性立竿见影。现象计算能跑完但结果里某个点的温度或位移是无穷值。原因材料属性被除零。比如流体储水系数引用了 1/孔隙度而盖层区域孔隙度正好为 0。解决给所有可能归零的材料参数加下界孔隙度用 0.001 作为下限。这种问题用探针检查具体节点就能定位不必怀疑求解器设置。4.2 物理量阶段饱和度震荡、热膨胀混淆与压力突变现象饱和度出现负值或者高于 1空间分布呈棋盘状震荡。原因时间和空间离散不匹配本质是 CFL 类条件被突破。两相流饱和度方程是双曲型问题对流项占优时中心差分格式必然震荡。解决加密近井网格同时把最大时间步限制在 0.01~0.1 天量级观察变化。如果仍然震荡把流动方程的对流项离散格式改为迎风格式COMSOL 里在物理场接口的离散化设置中可以直接切换。现象温度场和应力场同时偏大且注采井之间热突破时间明显早于预期。原因把热对流和热膨胀混为一谈。热对流是流体携带热量迁移需要传热接口与流动接口建立速度耦合热膨胀是骨架因温度变化产生应变。如果只在固体力学节点里加了热膨胀系数却没有让传热接口从达西流里取速度场等于只算了固体的热胀冷缩流体带来的热量搬运完全没参与计算。解决检查多物理场节点确认“传热-流动”耦合项存在。传热接口必须用多孔介质传热而不是固体传热。多孔介质传热会把骨架和流体的等效热容、等效导热系数一起加权这才是 THM 模型该用的传热方程。现象孔隙压力在井筒附近单点突变加密网格后压力剖面相去甚远。原因点源或线源直接加载到了单一网格节点上造成压力奇异性。这种问题在网格敏感性分析里最明显网格一加密井筒附近压力峰值大幅变化说明数值解还没收敛到物理解。解决改用小体积源区域或直接在井筒边界上给压力条件。如果模型几何里有井筒就在井筒壁面边界上施加压力而不是在中心线节点上打点。点源只适合离研究区足够远的辅助计算。5. 模型改造与结果验证油水体系切到水-CO2哪些参数必须动5.1 结果验证先做两件事质量守恒与网格敏感性拿到模型后先不要大改按默认参数跑一个基准算例验证两件事。第一流体质量守恒进出口累积流量与模型内累积量之间的差值绝对值要控制在 1% 以内。两相流 THM 里质量不守恒基本可以定位到时间步过大或数值弥散过强。第二网格敏感性把近井网格加密一倍观察压力剖面和饱和度剖面的变化幅度小于 5% 说明网格达到了工程可接受水平。这两步做完模型才算真正“被验证过”。5.2 油水切到水-CO2 的最小改动状态方程、相渗与注入条件这套模型也有水-超临界 CO2 两相流版本但只有油水版本时自己改造也不难三个地方必须动。第一CO2 的密度和黏度不能填常数超临界 CO2 密度随压力变化非常剧烈材料节点要用真实气体密度或状态方程插值表直接用常数会导致重力驱替和压力扩散全部失真。第二相对渗透率曲线要换油-水体系和 CO2-水体系的界面张力、润湿性都不同Corey 指数和残余饱和度需要重新标定沿用油水参数算出的气驱前沿形态不可信。第三注入条件从注油改成注 CO2 时注意 CO2 在储层温压下可能处于超临界态密度和黏度随局部压力变化会影响重力舌进二维剖面模型里看到明显的重力舌是正常物理现象不要当成不收敛。这套模型的改动空间很大材料、边界、求解器都能动但每改一处都要重新验证。我现在的习惯是每次拿到别人给的耦合模型先跑基准算例确认质量守恒和网格无关性再替换材料参数最后才动求解器。特别是换流体体系时强制自己把相渗曲线和状态方程全部核对一遍跳过任何一步模拟结果就可能变成一只漂亮但不可信的黑匣子。希望帮到你。本文还有配套的精品资源点击获取