COMSOL两相流THM热流固耦合模型:从物理原理到收敛性优化

📅 发布时间:2026/10/19 17:42:24
COMSOL两相流THM热流固耦合模型:从物理原理到收敛性优化
简介面向COMSOL仿真用户的THM热流固耦合两相流模型资源适用于石油工程、地热开发等领域中饱和水与油流动的数值模拟也可扩展至水与超临界二氧化碳两相流场景。模型完整涵盖热-流-固三场耦合用户可根据实际需求自由修改几何、材料参数及边界条件灵活性高。资源包共含14个文件以8个Word文档为主体详细讲解模型构建原理与结果分析辅以4张示意图直观展示网格与云图另有2个HTML页面方便浏览器快速查看关键结论压缩包整体仅3.91MB轻量易下载。已有133人学习使用收敛性表现优异计算速度快且鲁棒性强能够稳定获取可靠结果。无论是初涉多场耦合的入门研究者还是需要高效仿真工具的工程人员这套模型都能显著节省建模调试时间是实践THM两相流耦合分析的高效参考。1. 为什么说THM两相流模型不只是一个“能跑的算例”COMSOL两相流THM热流固耦合模型听起来像是某个课题的定制算例但它其实是地热开发、油藏工程、二氧化碳封存里最常用到的一类硬骨头。很多工程师在THM耦合上翻车往往不是因为流体模块不会建而是温度、压力、饱和度、应力四个场互相咬合稍有设置不当就发散或结果明显失真。这个模型把饱和水和油的两相流动、热传导、固体变形放在一个耦合框架里声称收敛性好、计算速度快、鲁棒性好实际使用中确实能省去大量调试时间。适合正在做油藏热采、地热储层改造、核废料处置等方向的人也适合想从单相流过渡到多场耦合的仿真用户。2. 先拆物理油水两相流与热流固三场耦合的数学结构在动手设置COMSOL之前我建议你先有控制方程在手边。不是说要手动推导整套公式而是要知道你用的每个接口到底在解什么方程哪些耦合项是内置的哪些需要你手动加。这一章就把油水两相流、传热、固体力学三块拆开再合到一起。2.1 流场饱和水和油的两相Darcy流动怎么进COMSOL饱和水和油的两相流动在低渗多孔介质里通常用两相Darcy流动来描述。每一相的达西速度为v_α - (k_rα K / μ_α) (∇p_α - ρ_α g ∇z)其中α代表水或油k_rα是相对渗透率K是绝对渗透率μ_α是动力黏度。饱和度方程则是∂(φ S_α)/∂t ∇·v_α q_α并且有S_w S_o 1。毛管压力把两相压力联系起来p_o - p_w p_c(S_w)。水相压力通常作为主要求解变量油相压力由毛管压力函数得到。在COMSOL里模型向导中选择“多孔介质多相流”接口在两相情况下会直接以压力和饱和度为因变量。如果你用的版本里找不到这个接口可以用一个“Darcy定律”接口加上一个对流-扩散PDE把饱和度作为因变量再自己定义相对渗透率和毛管压力。我一般会优先用现成的“多孔介质多相流”接口因为它已经内置了饱和度约束、相流量公式和存储项能少写很多容易出错的耦合项。设置时要注意下面几点流体材料水和油分别建材料密度和黏度都可以写成温度的函数。油类物质在温度升高时黏度下降明显这个变化对THM耦合影响很大不要用常数凑合。相对渗透率最常用的是Corey模型。比如水相相对渗透率写成S_w的4次方油相写成(1-S_w)^2具体指数根据岩心实验标定。没有实验数据时用这两个指数作为默认起点是常见的做法。毛管压力可以用van Genuchten模型也可以用简单的指数函数。重点是毛管压力函数要单调导函数不能为0否则饱和度方程会出现退化也是后面积分不收敛的原因之一。初始饱和度不要一拍脑袋给一个常数。最好先做一步“无流动初始化”让饱和度和毛管压力在域内平衡再开始瞬态求解。很多新手第一次跑就报饱和度越界就是忽略了这一步。这里我给出一个典型的参数起点表适合尚未拿到实验数据时的快速验证具体数值要按你的岩心或数值实验调整参数建议值/范围说明水密度1000 kg/m3可按温度修正油密度850 kg/m3可按温度修正水黏度1e-3 Pa·s随温度指数下降油黏度5e-30.1 Pa·s温度敏感建议用插值函数绝对渗透率100 mD约1e-13 m²各向异性时可分列孔隙度0.20.3会随有效应力变化水相残余饱和度0.2经验值油相残余饱和度0.15经验值Corey指数水相4、油相2非油藏实验时可用注意两相Darcy流动在COMSOL中的变量默认可能使用“多孔基体体积平均流速”或“Darcy流速”这会影响你加载的速率边界条件和后处理。我习惯用Darcy流速即单位面积体积流量这样和井的注入量、采出量换算比较直观。后处理时如果发现速度量级不对先检查这个单位约定。2.2 温度场与应力场THM三场耦合的耦合项在哪温度场在多孔介质里的控制方程是能量守恒(ρC_p)_eff ∂T/∂t (ρC_p)_fluid (u·∇)T - ∇·(k_eff ∇T) Q其中(ρC_p)_eff是基体和流体的体积平均热容k_eff是有效导热系数u是达西流速或者体积平均流速。注意这个方程里的对流项不是通常的流速而是多孔介质内的流体热容加权的速度。在COMSOL中使用“多孔介质传热”接口时它会自动把流体的速度场耦合进来你只需要指定有效导热系数和串联/并联模型。应力场部分就是线性弹性或弹塑性固体的力平衡方程∇·σ_eff F 0σ_eff σ - α_Biot p_w I这里采用有效应力原理α_Biot是Biot系数。温度的影响通过热应变进入本构方程线弹性本构里Δε_th α_T (T - T_ref)。热膨胀系数是体积热膨胀系数各向同性情况下等于线膨胀系数的3倍在COMSOL固体力学接口里按线膨胀系数输入。压力耦合则在多物理场“孔隙弹性”节点中自动加到应力张量上。耦合关系实际上有三个方向流→固孔隙压力和饱和度变化引起有效应力变化产生变形。热→流温度影响流体黏度和密度从而影响渗透和饱和度分布。热→固热应变产生位移位移又会通过孔隙弹性改变体积应变进而影响孔隙度和渗透率。固→流这步最容易被忽略。变形改变孔隙度孔隙度变化又改变存储系数对两相Darcy方程产生一个源项。如果在COMSOL里直接用一个“变量”把孔隙度改为有效应力的函数而不在流体接口的“存储”项中增加∂φ/∂p项会导致流体质量不守恒这是后面会讲到的常见坑。在COMSOL里实现三场耦合通常是在同一域上叠加三个物理场接口多孔介质多相流、多孔介质传热、固体力学。然后通过多物理场节点加入“孔隙弹性”和“非等温流动”或“热膨胀”耦合。如果你用的是较新版本的COMSOL多孔介质多相流接口本身支持与固体力学耦合的“变形多孔介质”设置确认在“多孔介质”节点下启用“孔隙率依赖于压力/应力”选项。一个稳妥的设置顺序是先单独跑通流动再跑通传热再挂上固体力学。三场全耦合容易把问题混在一起出问题时不知道是哪个场在发散。先两两耦合找到每个耦合节点的效果再全耦合是省时间的好习惯。另外一个常见疑问是“我能不能用用户自定义的PDE来写两相流避免‘多孔介质多相流’接口的黑匣子”我的看法是如果你需要实现特殊相渗模型或非标准毛管压力用PDE未必比现成接口慢。但PDE接口要自己处理保守形式和数值震荡调试成本高。对于标题里提到的“收敛性好、计算速度快、鲁棒性好”现成接口带来的帮助远大于损失。只有当你的流动方程与标准Darcy模型有本质区别时才值得自写PDE。这个判断同样适用于COMSOL模型中的其他场。3. 搭建模型从几何到求解器的关键设置物理拆完之后我们开始搭模型。这一章按我常用的顺序走几何与边界条件、物理场与耦合节点、求解器配置。每一步都直接影响你能不能真正“跑起来”而不是只在概念上成立。3.1 几何与边界条件简化不是偷工减料THM两相流模型本身已经够复杂几何上没有必要一上来就还原真实储层形态。常见做法是先用二维矩形或轴对称圆筒验证可行性再逐步增加断层、裂缝、夹层。一个经典的简化模型是水平方向100米、垂向50米的矩形域模拟一个注采井对或者只在左侧设注入井右侧设产层边界。这样的几何能让你先看清饱和度前锋、温度锋面和位移分布再谈工程细节。在COMSOL里创建几何时建议直接在矩形域中挖掉井筒区域而不是用点源/汇。因为井筒附近压力梯度大用点源会在后处理和网格加密时带来奇异问题。挖掉一个半径0.1米的小圆作为井孔然后用边界条件或汇项定义流量效果稳定得多。根据地层对称性可以使用轴对称模型把三维问题降到二维计算量少一个数量级这也是“计算速度快”的一个重要来源。边界条件的设置顺序是流动边界注入井壁给定流量或压力产层边界给定初始压力。如果模拟定压生产外边界压力固定如果模拟封闭油藏外边界为零流量。热边界注入温度低于或高于地层温度时井壁边界给固定温度外边界可以用隔热零热通量或固定地层温度。这个选择会影响温度锋面传播。力学边界底部固定左右两侧滚支顶部自由。如果模型包含盖层盖层顶面可以加一个自由面并受上覆压力。需要注意COMSOL中不同物理场对同一个几何实体可以使用不同的边界条件例如井孔边界既是流体的流入边界也是固体的自由表面边界。建立选择时建议为每个物理场单独命名边界选择比如“注入井壁_流体”防止后边调整几何时选错边界。说到几何参数我习惯把所有几何尺寸定义为参数而不是直接在绘图界面输入数字。这样后面做参数化扫描比如改变井距、水油比、注水量会非常方便。COMSOL中的参数扫描可以自动创建一组衍生解避免你手动复制模型。同一批参数还能在导入网格、设置材料时复用减少改一处漏一处的尴尬。3.2 物理场接口选择与耦合节点先建流场再挂温度与应力在模型向导里我一般按以下顺序添加物理场多孔介质多相流两相油和水多孔介质传热固体力学顺序不是必须的但流场先建因为温度和固体力学都需要用到速度场的表达。添加完成后在“组件”下会看到三个物理场节点和一个“多物理场”节点。在多物理场节点中加入以下耦合孔隙弹性耦合固体力学与多孔介质多相流。这里COMSOL会要求选择一个压力变量默认是流体压力。在多相流情况下通常采用水相压力作为等效孔隙压力。如果你的问题要求两相压力的各自贡献比如考虑毛管压力对应力的影响需要自定义一个等效孔隙压力变量比如p_eff S_w * p_w S_o * p_o。这个细节是决定THM模型是否精确的常见分歧点。非等温流动或热-流耦合耦合多孔介质传热与多孔介质多相流。它会自动把温度场中的流体热容项与流动速度结合你只需要确保流体材料的密度和热容随温度变化是合理的。热膨胀耦合固体力学与多孔介质传热把温度场作为热载荷。在设置物理场时最容易出错的是“界面条件”。比如注入井壁上的热边界实际可能同时存在流体注入和热通量。COMSOL允许在同一个边界上叠加多个物理场的边界条件但要避免“定义了流量又定义了固定压力”这会让求解器无所适从。我的习惯是边界上只对流动定义一个条件流量或压力热边界和力学边界各自独立设置不互相影响。材料设置上流体材料建议用内置材料库中的常用流体再修改为温度相关的属性。如果你已经有实测的黏温曲线用插值函数是最佳选择。COMSOL中可以用“插值”函数导入一个文本表格然后在材料属性里引用该函数。对于基体岩石各向同性时只填密度、杨氏模量、泊松比、导热系数、热容和热膨胀系数。如果做正交各向异性渗透率和弹性张量都要定义完整复杂度会立刻翻倍不建议新手一上来就挑战。注意两个物理场共享同一个几何域时要确保网格域选择一致尤其不要在一个物理场中选择“所有域”另一个只选择部分域这会在耦合时产生无效变量。3.3 求解器配置为什么默认求解器会卡手动调什么模型建完后求解器默认会自动选择全耦合或分步。对于THM两相流默认全耦合常常导致内存占用大、迭代发散。我一般会手动切到“分离式求解器”自定义求解顺序先解流体压力和饱和度再解温度最后解位移。这一步在哪里改取决于版本大体上是求解器配置→分离步骤→添加三个分离步分别选择对应的因变量。分离式求解器还有几个参数值得调最大迭代次数默认可能只有几次。两相流加THM耦合经常需要较多分离迭代我通常设到25~50次并在“高级”里开启“一致初始化”或“全局牛顿”选项。阻尼因子如果压力残差震荡把阻尼从1改成0.7~0.9能明显提高稳定性。但阻尼过小会让收敛变慢建议从0.9开始往下试探。因变量标度压力和位移的量纲差异很大。COMSOL默认会自动标度但如果残差提示“因变量标度不可靠”手动把压力标度设为1e6以Pa为单位位移设为1e-3以m为单位流动时间按你的时间尺度标定。时间步长瞬态求解时COMSOL的自适应时间步在THM耦合中容易“冒进”导致某个时刻残差突然增加。建议设置最大时间步为总模拟时间的1/50或更小再开启“恢复早期解”功能。这样可以避免前锋扫过时步长过大。网格方面三相耦合模型对网格质量很敏感。粗网格虽然快但饱和度前锋和应力集中区会被抹平细网格又会导致每个时间步的牛顿迭代代价剧增。我通常先跑一个粗网格短时间比如1小时的模拟检查场分布是否光滑、饱和度是否出界、位移是否合量级再针对注入井附近做一层边界层网格加密。记住网格无关性验证不是最后一刻才做的事它要伴随参数调优的全过程。4. 参数怎么设物性、初始值、收敛判据的实用清单THM模型的“玄学”很多时候不在模型结构而在参数。参数给得不合理求解器再强大也救不回来。这一章把油水两相流、热流固耦合里最关键的那些参数列出来并说明它们对收敛和速度的真实影响。4.1 油水两相流的相渗曲线与毛管压力参数相渗曲线是两相流模型的灵魂。和单相渗流不同两相流中每一相的有效渗透率都随饱和度变化如果相对渗透率为零该相在网格内的传输就关闭容易造成数值不连续。COMSOL里可以用解析函数定义S_w_star (S_w - S_wr)/(1 - S_wr - S_or) k_rw (S_w_star)^nw k_ro (1 - S_w_star)^no其中S_wr和S_or分别是束缚水和残余油饱和度。这个模型的优点是简单参数少数值稳定缺点是没有体现相渗曲线在低饱和度端的缓变特征。如果你手头有实验测得的曲线建议用COMSOL的“插值”函数导入并启用“端点外推”到零。毛管压力函数常用van Genuchten形式p_c (1/α_vG) (S_e^{-1/m} - 1)^(1-n_vG)其中S_e是有效饱和度。需要注意毛管压力函数的导数在饱和度趋近残余饱和度时会变得非常大这会让线性化矩阵病态。一个规避方法是在饱和度接近端点时对p_c函数做平滑截断或者给p_c加一个很小的最小斜率。工程上常说“毛管压力是两相流不收敛的主要来源之一”我深有同感。另一个容易忽略的是初始饱和度与毛管压力的协调。如果你给初始饱和度S_w0.5但按毛管压力函数这个饱和度对应的两相压差为某个正值而你又给初始压力为静水压那么求解器要在一个极小时间步内把压力场和饱和度场重新平衡很容易震荡。解决方案是先跑一个“无流动初始化”求解步关闭注入流量只让压力和饱和度在毛管压力平衡下稳定下来然后用这个稳定解作为瞬态初值。4.2 热-流-固耦合的关键物性孔隙度、渗透率、热膨胀系数三场耦合里孔隙度是承上启下的枢纽。它不仅控制流体存储量还影响导热系数和弹性模量。很多建模者把孔隙度设为常数这简化了问题但当压力变化大时孔隙度的变化会显著影响有效应力反过来变形也会改变储层渗透性。在COMSOL中你可以定义孔隙度随有效应力的线性变化φ φ_0 - c_φ * (σ_eff - σ_eff0)其中c_φ是孔隙度对应力的敏感系数。但要注意这样定义后质量守恒方程中的存储项必须包含∂φ/∂p贡献。更稳妥的方式是利用COMSOL的“孔隙弹性”节点自动生成的孔隙度表达式而不是手动写变量。渗透率与孔隙度的关系最常用的是Kozeny-Carman式K K_0 * (φ/φ_0)^3 * ((1-φ_0)/(1-φ))^2这式子虽是经验关系但它能保证渗透率在孔隙度下降时快速下降用于应力敏感储层建模很顺手。如果你观察到仿真结果中压力分布对渗透率极敏感先检查这个关系是否在合理范围而不是怀疑求解器。热膨胀系数对THM耦合的稳定性影响很大。如果只做常温渗流热膨胀系数对结果的影响可以忽略但如果模拟热采或热流注入温度变化几十度到几百度线膨胀系数每增加1e-5 1/K热应变就可能超过弹性应变单元扭曲就不再是罕见事。所以热膨胀系数不要取“典型岩石值”就完事要根据你的岩性和温度变化范围来评估。另外导热系数在干、湿状态下差异很大建议使用饱和状态下的有效导热系数用体积平均或串联并联平均而不是用干岩值。以下是常用物性参数表数值用于初步估算参数岩石骨架水油密度 kg/m326001000850热导率 W/(m·K)2.50.60.15比热容 J/(kg·K)90042002000线膨胀系数 1/K1e-52e-45e-4注意流体的线膨胀系数在传热里通常不起主导作用但如果你用“热膨胀”耦合固体力学只用固体的热膨胀系数即可。4.3 计算速度与鲁棒性时间步长、迭代容差、网格粗化标题里强调了模型收敛性好、计算速度快、鲁棒性好这其实和求解器参数、网格策略密切相关。时间步长是最直接的调节开关。瞬态THM耦合中饱和度前锋和温度锋面以不同速度传播如果时间步长大于局部传播时间就容易出现振荡或负的饱和度。建议采用最大时间步长限制同时允许自适应步长在内缩。通常把最大时间步设为总模拟时长的1/1000作为起点如果发现太慢再逐步放宽到1/500。迭代容差方面相对容差设为1e-3通常足够工程应用设为1e-4会明显增加迭代次数且不一定让结果更准。不要一开始就追求1e-6那是数值洁癖。我通常使用分离式求解器时压力和饱和度步用1e-3温度和位移步用1e-4但COMSOL的容差设置是全局统一或分物理场具体可以调整因变量的容差加权新手建议保持默认。需要知道的是若模型反复在“牛顿迭代不收敛”上报错先把相对容差放宽到5e-3试跑看是不是容差设得太苛刻。网格粗化是“计算速度快”的重要来源但盲目粗化会让结果失真。常见做法是注入井附近用边界层网格或局部加密远离井的区域可以按比例放粗。一个矩形100m*50m的域井周围0.1m网格用10~20个单元中间区域用0.5~2m的最大单元总单元数控制在几千到一两万之间就足以跑通二维模型。三维模型则要小心几十万单元在THM耦合下内存会爆炸建议先用轴对称或薄层模型验证。5. 避坑THM两相流模型最常见的5个翻车点这一章是血泪经验合集。每个问题我都按现象、原因、解决来写方便你对照排查。翻车不可怕怕的是不知道往哪个方向查。5.1 油水两相流初始饱和度的“黑匣子”问题现象模型一开始就报错或者饱和度在第一个时间步出现负值、大于1。原因初始饱和度与毛管压力、初始压力分布不协调。COMSOL在求解两相流动时饱和度和压力必须同时满足毛管压力关系。如果给一个均匀饱和度值但压力场又是静水压那么在边界附近会强制产生一个压力突变导致初始状态不是平衡态。解决先做“无流动初始化”求解步。具体做法是在瞬态求解前先冻结注入/采出条件用一个小时间步让流动场自行平衡。也可以在“多孔介质多相流”接口中开启“初始化”功能自动求解静水平衡。另一个补救措施是把初始饱和度定义为随压力变化的毛管压力反函数保证初始压力与饱和度在同一毛管压力曲线上。5.2 流固耦合不收敛孔隙度更新导致的水量不平衡现象分离式求解器在“流场”步迭代次数持续增加压力残差不降甚至提示“质量守恒残留超过容差”。原因固体力学更新了孔隙度但流场方程的质量存储项没有同步更新。例如把孔隙度定义成位移的函数却仍使用常数孔隙度的存储系数那么压力变化引起的流体质量变化和力学变形不一致。解决不要在流体接口里手动把孔隙度写成位移函数。正确做法是使用“孔隙弹性”多物理场它会自动引入Biot系数和体积应变对存储项的贡献。如果是新版COMSOL在多孔介质多相流“多孔介质”子节点中开启“孔隙率表达式”并选择“应力相关”选项。改完后要重新检查流体物质平衡方法是在后处理中计算域内平均饱和度变化与边界累计流量的差值。5.3 热应力把单元拉成负体积几何非线性忘开了现象传热求解正常但固体力学在温度变化较大时出现“未定义结果”或单元雅可比为负。原因线性几何假设下当热应变和孔隙压力引起的变形超过小变形范围仍使用未更新几何就会导致应力应变关系失真甚至单元畸变。解决在固体力学接口的设置中把“几何非线性”选为“自动检测”或强制开启。开启大变形后多孔介质流动区域是否需要移动网格取决于变形量级如果变形超过几个百分点最好把“变形”节点加到多孔介质多相流中让流体区域跟随固体变形。否则可以分步加载温度例如每步只加20度让求解器更容易收敛。5.4 超临界CO2版本的特殊坑状态方程与密度突变现象把油相换成超临界CO2后瞬态求解在临界压力附近压力振荡收敛情况和鲁棒性都劣于油相。原因超临界CO2的密度、黏度和焓在临界点附近随温度和压力变化剧烈。若用常数或者简单的线性插值会导致控制方程的非线性过于强烈线性化迭代一越界就发散。解决换用温压相关的真实气体物性。常见做法是在外部查询表中按温度、压力网格取密度和黏度导入COMSOL作为插值函数或者用解析的Peng-Robinson状态方程配合COMSOL“理想气体”之外的“真实气体”模块。如果临界点区域不重要还可以在设置中将工作区限定在远离临界点的压力区间减少非线性。注意CO2在水中的溶解度也会改变流体性质但两相流模型如果不考虑组分溶解先不要加这一项否则模型复杂度会失控。5.5 网格无关性验证不是形式主义现象把网格加密一倍后关键点温度或位移变化较小但饱和度锋面位置差出几个网格宽度甚至压力云图差异明显。原因两相流方程具有迎风特性和锋面在粗网格下数值扩散会自动抹平饱和度前锋。应力场在线弹性问题中往往受网格影响较小所以只看应力结果会让网格无关性验证得到假阳性。解决做网格无关性验证时同时监测井底压力、生产端流量或含水率曲线、剖面饱和度曲线而不是只看某个点。加密网格时采用等比细化并保证井周边界层网格数同步增加。如果饱和度剖面在加密后仍有明显差异检查毛管压力和重力数可能问题处在物理参数而不是网格。网格无关性从来不是一步就能完成先跑粗网格、中等网格、细网格三个级别再决定最终网格密度。6. 验证与进阶用解析解校准模型再改成水CO2模型能跑不代表结果正确。在放大工程应用之前我建议先用解析解或半解析解给模型“验身”。验证通过后再改造成水超临界CO2版本才算把标题里的能力真正握在手里。6.1 用一维Terzaghi固结解做热流固耦合验证对于流固耦合最经典的一维Terzaghi固结问题可以验证孔隙压力和位移的瞬态行为。构造一个宽度很窄的柱状域顶部自由并施加压力载荷底部固定且排水边界时间演化中孔隙压力消散过程与固结解析解对比。在COMSOL中关闭两相流用一个Darcy接口和固体力学耦合按有效应力原理求解观测某点孔隙压力随时间变化。两者偏差在5%以内说明“流固耦合”的核心路径正确。再做传热验证可以加一个恒定热源或边界温度对比一维热扩散解析解。只有这两步都通过才允许自己继续调油水两相流参数。6.2 从油水两相流扩展到水超临界CO2两相流的改动清单扩展时不要复制整个模型再改流体材料。我一般会按以下清单操作将油的材料节点替换为CO2定义密度和黏度作为温度和压力的插值函数或真实气体表达式。检查相对渗透率和毛管压力曲线。水和CO2的相渗曲线与油水不同尤其是在超临界条件下界面张力和润湿性变化明显需要更换参数。修正热物性CO2的比热容、导热系数随温度和压力变化大不能用常数。热容的插值尤其重要否则能量方程可能不守恒。关闭不适用的大变形假设。如果模拟的是封存场景注入压力高变形大几何非线性要保留。重新初始化饱和度场避免用油水模型的初值。CO2密度低重力作用和重力稳定性与油完全不同。我在做这类模型时有一个习惯先不急着上全耦合把每个物理场单独跑通再两两耦合最后三场全耦合。这个习惯让我少踩了很多坑。THM两相流模型不是一次就能跑得漂亮的“一键算例”但只要你按本文的顺序拆解用解析解验证每一环再按避坑清单排查收敛和鲁棒性是可以建立起来的。希望帮到你。本文还有配套的精品资源点击获取