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

简介一份基于COMSOL的两相流THM热流固耦合数值仿真模型面向油气渗流、地热开发与地下多相流模拟领域的工程师和科研人员也适合用作高校及科研机构数值仿真教学课程中的进阶案例。模型以饱和水和油的两相流动为对象同时耦合热传导、流体渗流与固体力学过程支持按实际需求调整物理场、材料参数、边界条件或几何结构还可进一步扩展至水与超临界二氧化碳的两相流体系。压缩包共14个文件包括8个Word文档、2个HTML文件与4张图片整体大小3.91MBWord文档系统介绍两相流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 材料参数表六个必填项与两个易错项材料参数是多场耦合的数据地基。下面是初学者可直接抄的参数表更新为应用相关的数据参数符号数值说明孔隙度φ0.2均匀场暂不做非均质渗透率K1e-13 m²约100 mD中渗砂岩水相粘度μ_w1e-3 Pa·s20°C时按温度插值油相粘度μ_o2e-2 Pa·s20°C时温度升高则下降水密度ρ_w1000 kg/m³可加微压缩性油密度ρ_o800 kg/m³浮力驱动取决于差值水热导率k_w0.6 W/(m·K)与温度略相关油热导率k_o0.15 W/(m·K)骨架热导率k_s2.5 W/(m·K)砂岩典型值水比热c_p,w4200 J/(kg·K)油比热c_p,o2100 J/(kg·K)骨架比热c_p,s900 J/(kg·K)骨架密度ρ_s2650 kg/m³弹性模量E10 GPa岩石模量泊松比ν0.25比奥系数α_B0.8孔隙弹性耦合热膨胀系数α_T1e-5 1/K骨架线性热膨胀束缚水饱和度S_wr0.2相对渗透率曲线端点残余油饱和度S_or0.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中跑起瞬态研究的流程实际操作流程按以下步骤执行建立矩形几何设置注水井点与生产井点。定义全局参数包括注入压力、地层温度、孔隙度和渗透率。添加三个物理场接口按照2.1到2.3节设置方程来源。建立多物理场耦合节点确认孔隙弹性与热膨胀被激活。设定初始值压力写p0饱和度写S0温度写T0。网格划分并加密井附近区域。研究设置为瞬态时间范围从0到3.15e7秒一年步长设置上限3600秒。求解器用默认的“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后处理里对出入口做质量积分、热量积分确认总累计注入量减产出量等于储量变化量。同时在参数扫描结果中画出“温度突破时间随注入压力”的曲线观察趋势是否单调。有一次我对比扫描结果发现温度突破时间没有随注水量增加而显著提前排查后发现是入口边界误设成了无热流条件热量交换根本没发生。这种问题靠云图很难发现靠趋势曲线一眼就看穿。我的血泪教训是不管多急着出图先把这三个验证动作练成习惯。模型自己没有对错只有你还没发现的不一致。希望帮到你。本文还有配套的精品资源点击获取