干微观驱替模拟的人都知道真正影响结果走向的往往不是软件操作本身而是对物理过程的拆解和对数值稳定性的控制。我最近在COMSOL里做了一件事从岩石扫描图像中重建出真实的孔隙几何再模拟两种油相在这个孔隙结构里的扩散与驱替过程盯着前缘怎么运动、剩余油会在哪里滞留、不同润湿条件下驱替效率差多少。这套流程听起来不算宏大实际完整跑下来涉及的东西非常多从几何清理到相场参数标定任何一个环节掉链子结果都没法看。这篇稿子把整套过程完整写下来包括模型设计思路、关键参数设定、求解器选择、后处理分析以及我在调试过程中踩过的几个坑。适用对象是正在做多孔介质微观渗流、两相流数值模拟或者刚入坑COMSOL流体仿真想找参考的同行。下面直接进入正题。1. 内容整体设计与思路拆解1.1 “两油扩散与驱替”在物理上到底指什么把这个前提说清楚后面的参数才有讨论的意义。所谓两油其实就指两种物性不同的流体相最常见的是黏度、密度各不相同的两相油品也可以把体系换成油和水方法完全通用。我的项目里设定油相1为低黏度侵入流体油相2为高黏度被驱流体用来模拟注气或注低黏介质驱高黏原油的场景。从物理本质看驱替就是把一种流体“推”着取代另一种流体在孔隙中的位置。我们可以做个生活类比用吸管朝杯子里吹气气流把杯中果汁推向杯壁往外溢吹气的能量通过流体的压差传递到相界面上这就是一次微型驱替。发生在岩石孔隙里的驱替要复杂得多因为通道不是一根直管而是粗细交替、坑洼不平、带有死角的复杂网络。流体走到细喉道时界面会变形并产生毛细阻力过去之后又在孔隙里急剧膨胀前缘形状会变得非常不规则这也是微观驱替模拟比宏观渗流计算有趣的地方。标题里的“扩散”同样不是可有可无的字眼。COMSOL的相场法在求解界面运动时会天然包含扩散机制更直白的场景是在驱替基本结束后关闭入口流动让相场继续演化一段时间观察滞留油是否因分子扩散产生缓慢的组分迁移。这一步能把“驱替”和“扩散”两件事结合起来体现多场耦合的价值。1.2 为什么选COMSOL而不是OpenFOAM或其他工具我在做方案选型时纠结过一阵。开源CFD工具里的多相流求解器也能做这类模拟功能并不弱但COMSOL在这个场景下有几点优势非常明显参数化建模方便。材料属性、边界条件、控制方程耦合都在同一个界面里管理改一个变量全局自动更新相场接口成熟。自带的“两相流相场”物理接口集成度高不需要自己去拼代码里的方程组后处理灵活。积分算子、派生值、动画导出全都在一个工作区完成不用在求解器和后处理软件之间反复倒数据。代价是求解速度会比专门的多相流求解器慢一些。但微观孔隙模型规模通常不大百万级网格在单机上是可以接受的这种便捷性换一点计算时间非常划算。1.3 为什么非要“实际孔隙”简化模型不行吗很多初学者习惯用平行毛管束或单根变径管代替岩石孔隙算出漂亮的曲线后直接下结论。但这类简化模型有一个致命问题它会把真实结构最核心的特征也就是孔喉比、迂曲度和盲端孔隙全部抹平。比如一根直管里只要压差存在被驱流体几乎可以被推干净可在真实孔隙里流体沿大孔隙走捷径把小孔隙里的油绕过甚至反向夹住形成非连续的残余油。这种残余油靠简单几何是算不出来的。我做过一个对比相同边界条件下直管模型的驱替效率超过90%而同一参数放进真实CT重建的孔隙模型里驱替效率只有65%左右。这个差距足以说明为什么“实际孔隙”四个字是整个课题的基石。这也是项目名称里特别强调“实际孔隙”的根本原因——不还原结构就谈不上研究真实驱替机理。2. 几何模型从切片到可计算域2.1 孔隙结构从哪里来CT扫描与数值重建两条路要还原实际孔隙第一步是拿到孔隙空间的几何表达。主流途径有两条第一是CT扫描重建。把岩心柱放入显微CT设备得到一组二维切片序列经过滤波、阈值分割、二值化等图像处理后把岩石骨架和孔隙空间分离开。孔隙空间最终导出为STL或面网格文件COMSOL直接导入即可。CT扫描的优点是结构真实缺点是需要设备而且扫描分辨率直接决定可以分辨的最小孔喉尺寸。第二是数值重建。如果没有扫描数据可以用颗粒堆积法、四参数随机生长法等算法生成和真实岩心统计特征相近的数字岩心。这种方法生成的结构不是真实岩心但能反映孔隙尺寸分布、配位数等关键统计特征对于机理研究完全够用而且不用处理噪声数据入坑速度更快。我这次优先用了CT扫描数据主要原因是需要在真实结构里验证驱替过程。但为了前期调试物理场参数我先用数值重建生成了一块小尺寸的孔隙模型参数跑通后才把CT几何替换进来。这个策略节省了大量时间强烈建议参考。2.2 几何清理与导入别让STL的锯齿毁掉你的网格从图像里提取出的面网格通常有很多问题表面锯齿、破损面、相互穿透的面片、退化单元。这些问题如果不处理导入COMSOL后会造成几何修复失败或者后续网格质量极差。在COMSOL里做导入我的操作路径是在“几何”节点下选择“导入”载入STL文件然后开启“修复容差”参数。修复容差决定了多少范围内的缝隙会被自动缝合一般先取全局尺寸的1%左右试算看修复报告里的剩余破损面数量。如果破损面过多千万不要靠无限调大容差硬来那样会把细小孔喉直接填死结构失真。正确做法是回到图像处理阶段把STL重新做一次表面光滑和三角形简化。导入成功后做几件事删除退化面把面积过小、边长比例极端的表面片去掉布尔并集运算把多个独立导入的实体合并成一个整体域避免后续物理场装配的边界混乱检查连通域数量用“几何属性”查看孔隙域是否整体连通。实际岩心中本来就可能存在孤立孔隙但如果孤立孔隙太多到难以统计可能是阈值分割时把骨架和孔隙分错了。2.3 网格划分孔隙尺度模拟成败的第一道关口网格对相场模拟的影响比常规单相流要大得多。因为相场方法需要用一个有限厚度的过渡层去描述相界面这个厚度必须由网格解析。界面厚度参数ε一般取2到3个网格尺寸比如我设定ε5μm那么最大单元尺寸要控制在2.5μm以下。我最终的网格方案是在孔隙主体使用最大单元2.5μm的自由四面体同时在孔喉收缩段和预计界面扫过区域做局部加密壁面附近加12层边界层网格。整体网格量在百万级从孔隙尺度来看相当于用约40个网格分辨一个典型孔径。这个密度不算奢侈但已经很稳定。网格划完后一定要做一次网格无关性验证这也是我在项目里坚持做的事。方法是分别用最大单元尺寸1μm、2.5μm、5μm跑同样的算例对比驱替效率曲线。如果最粗和最细的结果差异在5%以内取中间尺寸就能满足计算精度要是差异明显就不能盲目信任当前网格。这一步看似费时实际上能帮你避开无数“假收敛、真误差”的尴尬。3. 控制方程与关键参数设定3.1 流体运动方程纳维-斯托克斯方程打底孔隙尺度两相流的流动本质由不可压缩纳维-斯托克斯方程控制再叠加相场方程去追踪界面。在COMSOL的“层流两相流相场”接口中这两组方程被耦合在一起求解。不可压缩NS方程的标准形式这里不写全式了关键要理解各物理项在孔隙尺度的角色。惯性力项在这种微米尺度流动中通常可以弱化因为雷诺数非常低。以我用的参数为例特征流速1mm/s特征孔径100μm油相黏度0.1 Pa·s密度850 kg/m³雷诺数约为0.001。在这样的量级下流动几乎完全由黏性力和毛细力主导前缘运动非常缓慢且平滑不会出现宏观流动里的湍流涡旋。这也决定了求解器的一个重要特征你不需要把时间步压到极小去追赶物理上的高速变化但会因为相场方程的强非线性存在数值稳定方面的挑战。3.2 Cahn-Hilliard相场方程界面厚度与迁移率的经典博弈相场方程比NS方程更容易让人翻车因为里面两个关键参数的取舍直接影响界面形态和计算稳定性。Cahn-Hilliard方程的基本形式可以理解为相场变量φ的对流输运加上一个由化学势梯度驱动的扩散项。在这个方程里表面张力σ和混合能密度λ之间存在换算关系σ (2√2/3) · λ / ε也就是说给定目标界面张力和界面厚度后混合能密度就确定了。反过来说如果随便填一个λ算出来的界面张力会和实验值差出好几倍。界面厚度ε不能太大否则界面过度模糊会把孔隙里本应清晰的前缘演变成一大片过渡区驱替特征完全失真也不能太小否则网格解析代价爆炸。实践中ε取1.5到3倍网格尺寸。迁移率γ同样关键。它控制着相场界面的数值扩散速度。γ太小质量守恒好但界面演化迟缓有时在复杂孔喉处会卡出假残留γ太大界面被快速抹开前缘锐度消失。我的经验是从1e-8 m³·s/kg起步试算如果界面前缘出现细碎“毛刺”小幅降低γ如果计算时间步长持续萎缩且界面不前进说明乘积项不合适再逐步回升。3.3 参数表一个可以直接抄作业的案例下面这组参数是我在项目中实际使用过一轮的基准值单位全部统一为国际单位制。你可以直接用这套数据做复现再根据自己的体系调整。参数名 表达式 描述 D 100[um] 特征孔隙直径 U_in 0.001[m/s] 入口速度 mu1 0.1[Pa*s] 油相1动力黏度高黏 mu2 0.001[Pa*s] 油相2动力黏度低黏 rho1 850[kg/m^3] 油相1密度 rho2 1000[kg/m^3] 油相2密度 sigma 0.03[N/m] 界面张力 epsilon 5[um] 相场界面厚度 gamma 1e-8[m^3*s/kg] 相场迁移率 theta_w 60[deg] 壁面接触角这组数据对应的毛细数为Ca μU/σ代入后约为0.0033刚好落在毛细驱替向黏性驱替过渡的区间。如果你希望不同驱替模式对比调整入口速度或黏度即可让Ca数量级发生变化。3.4 边界条件与初始条件最容易被忽视的物理陷阱几何域确定后边界条件设置直接决定流动模式是否正确。我的设置方式如下入口边界给定速度U_in同时指定相场变量φ1代表油相1进入出口边界给定压力为0并开启“抑制回流”选项防止数值回流把出口重新污染壁面边界设置无滑移速度条件同时设定接触角θ_w这是模拟润湿性对驱替形态影响的关键参数初始条件整个孔隙域初始化为φ-1代表初始完全充满油相2入口附近做一段缓冲让油相1可以从边界持续进入。这里有个经验必须单独强调不要让入口速度从0直接跳到目标值。比如U_in1mm/s的入口条件直接加载到一个充满高黏流体的静止域上压力场会瞬间冲击界面导致第一波计算步崩溃。我通常加一个时间相关的过渡函数让入口速度在0.01s内线性爬升到目标值。这个过程能显著减少初始震荡而且对最终驱替结果没有影响因为0.01s相比总模拟时间0.2s只占一个零头。4. 求解设置与后处理分析4.1 瞬态求解器全耦合还是分离式COMSOL里多相流问题一般选择瞬态研究。求解器有两种路线全耦合和分离式。我试过全耦合对小模型收敛很稳但在网格较大或时间步较长的阶段雅可比矩阵规模大会让单步时间显著拉长。后来我改用分离式求解器按“速度场-压力场”和“相场变量”分组迭代更新。分离式每次只解一个子块内存压力小而且可以通过调整子步迭代次数来稳定界面演变。分离式需要多迭代几轮才能达到同样的收敛残差但每个迭代步又快又稳整体算下来通常比全耦合省时。时间步长设置为自适应初始步长1e-4s最大步长0.001s后续交给软件自动调整。Courant数是一个重要的稳定性参考量它表示一个时间步内流体前缘移动的距离与网格尺寸的比值通常希望它保持在1附近。用我的参数算U_in1mm/s最大时间步0.001s前缘单步移动1μm与网格尺寸2.5μm同量级因此自适应时间步不会触发稳定性崩溃这点在做初始设置时可以自己估算。4.2 收敛判断与质量守恒检查判断计算是否真的收敛不能只看求解器报告“已收敛”。我最常做的检查有两项第一项是监察残差曲线。COMSOL求解器面板里可以勾选残差历史横坐标迭代步数、纵坐标残差范数。合格曲线应当持续下降或震荡下降至收敛容差以下而不是反弹后停在某个值附近不动。如果残差在1e-2附近长时间徘徊多半是界面相关条件设置出了问题。第二项是质量守恒。在域内积分相场变量或者更简单地对比入口与出口的体积流量差。如果累计流量差超过入口总流量的2%就要回头检查收敛和网格。我遇到过看起来完全收敛、但出口体积流量和入口差5%的算例最后发现是出口边界上的“抑制回流”选项没有打开数值回流导致假“增油”。这类问题大多可以靠边界条件修复而非调小容差硬撑。4.3 看动画比看曲线更能发现问题这是最想分享的经验之一。定量曲线当然重要但在模型调试阶段界面运动的动画才是真正帮你发现问题的手段。把几组时间截面的相场云图按时间序列播放前缘是否稳定、是否产生指进、油相1是否在孔喉处被截断一眼就能看出来。我在项目里导出了每隔0.005s一帧的相场分布图然后合成动画。第一次播放就发现被驱替流体在一个大孔隙的盲端被反复夹断形成孤立油滴。这个现象在饱和度曲线里只能看到平台期但动画里看机理一目了然。强烈建议每一组参数组合都保存动画很多后续论文或汇报用的机理分析素材就是这么来的。4.4 定量提取指标饱和度、前缘位置与残余油占比后处理阶段不能只看云图需要提取可比较的量化指标。平均饱和度在“派生值”中用积分算子对相场变量φ做域积分再除以区域体积。因为φ从-1到1表示不同相可以换算成油相1的体积分数前缘位置沿主流方向取相场φ0的等值面位置记录每时刻的推进距离残余油占比模拟结束后统计φ0区域中仍构成连通体的体积占比。配合“连通域分析”可以区分哪些残余油是盲端孤立油哪些是可动但未驱走的油。这些指标做出来后不同接触角、不同毛细数之间的对比就有依据了。用表格把静默数据变成趋势描述也方便后续项目总结。5. 算例分析润湿性、毛细数与扩散效应的面面观5.1 不同润湿条件下的驱替差异接触角是润湿性最直接的量化指标。我分别算了θ_w30°、60°、90°、120°四组。接触角30°代表强亲水界面张力对壁面有很强的“贴壁”倾向流体沿壁面铺展更快前缘推进相对均匀总体驱替效率最高。接触角120°后孔隙几何里很多曲面位置出现界面钉扎油相1的前缘推进非常不均衡产生明显的指进形态驱替效率下降显著。这组对比比较直观在不同方向工作时为润湿性调控方向的选择提供了依据。5.2 不同毛细数下的驱替模式转变毛细数的改变可以通过调整入口速度实现本质是在改变黏性力和毛细力的比值。我做了从0.001到0.1的量级扫描。低毛细数区界面张力占主导孔隙孔喉对前缘的束缚作用强前缘容易在细喉处被“锁住”驱替表现为非连续的跳跃式推进高毛细数区黏性力压制界面张力效应前缘推进更连续、更均匀但容易在非优势通道形成纵向指进。在孔隙尺度模拟中重现这种驱替模式的转变是从定性描述走向定量分析的关键一步。后续可以进一步绘制驱替效率随毛细数变化的散点图找到本模型中的过渡区间。5.3 驱替结束后的扩散驰豫不可忽视的“二次效应”驱替主过程结束后我并没有直接终止计算而是把入口速度降为0继续演化0.5s的物理时间观察相场继续扩散驰豫的状态。这时流体不再整体推进但相界面的梯度驱动扩散仍然存在高黏流体在盲端中的组分分布会发生缓慢重分布。这个算例的发现是一部分在驱替过程中被“夹断”的孤立油滴相邻的浅色区域主流通道路内浓度梯度持续存在经过扩散演化油滴界面发生了轻微移动部分死角中的油相浓度被稀释。虽然这个过程在短时间内对采出程度贡献有限但对于长时间尺度的采收率分析却有参考价值。这也是将“扩散”纳入项目标题的意义所在——区块的采收率不仅取决于驱替压差和界面性质还取决于组分在孔隙内的扩散能力。6. 常见问题与排查技巧实录6.1 界面前缘不光滑出现细碎“毛刺”这种情况在初期调参时几乎都会遇到。表面看是网格不够细仔细排查后才发现是迁移率γ取值偏大导致相场界面临近振荡。解决办法将γ降低半个数量级从1e-8降到3e-9然后再看是否出现毛刺。如果降低后前缘变得平滑但推进速度变慢说明γ仍处在合理范围内只是还没有收敛到我需要的形态。关键是要找到形态稳定与扩散速度之间的平衡点。6.2 时间步长越算越小最后算不下去这是相场模型最头痛的故障。出现这种问题先看提示中是否出现界面尖角或压力剧烈变化。常见触发原因是入口速度阶跃启动造成前缘位置的压力尖峰。对策是前面提到的入口速度爬升函数把0.01s的线性过渡加到所有初始参数上。另外要留意界面厚度ε与网格尺寸的比值如果ε小于网格尺寸约1.5倍相场方程对网格的灵敏度会急剧上升时间步也随之萎缩。此时稍微增大ε或加密局部网格问题往往会缓解。6.3 出口回流导致饱和度曲线异常前面提过出口边界如果不开启“抑制回流”求解器会在出口位置产生轻微回流以满足压力条件。这种回流会让出口处的相场波动导致累计流量统计不准确。排查方法很简单绘制出口边界上的速度分布如果边界某区域速度方向为负那就是回流。开启抑制回流选项后大多数情况都能解决。如果开了抑制回流仍异常需要检查出口是否紧邻一个特别狭窄的孔喉段可以在出口与孔喉之间人为加一段延长区。6.4 剩余油驱不走是数值问题还是物理现象模拟运行后把油相2“锁”在某个盲端第一反应经常是“网格不够、参数不对”。我的判断方法是做三个交叉验证一是把最大网格尺寸缩小一倍重新算二是把迁移率改变一个数量级重跑三是把入口速度提高10%再试。如果三种改变下残余油仍然稳定存在于同一位置那基本可以认定这是真实的物理现象是孔隙形态和润湿条件共同决定的残余油模式。如果位置随网格或参数漂移则是数值伪影需要回到网格与参数设置继续调。这个方法帮我区分了大量假问题也让后续分析更有底气。6.5 避坑清单速查入口速度一定要做爬坡过渡不要阶跃加载ε与网格尺寸比值保持1.5倍以上界面才不会被网格离散破坏迁移率过大是界面毛刺的头号元凶出现毛刺先降γ每次参数变更前先保存一份当前算例多参数对比时命名按“接触角_毛细数”的格式避免后期数据混乱求导入几何时修复合并不等于最终几何是正确的一定要用“几何属性”检查域体积是否合理饱和度的积分计算要选对相场变量符号φ1与φ-1表示两相符号反向结果完全不同动画输出分辨率不需要太高云图加流线叠加的模式最实用兼顾信息量和文件大小。我个人这段周期的体验是COMSOL跑这类孔隙尺度两相流并不难难的是永远知道当前结果是在反映真实物理还是放大了数值噪声。要掌握这个判断力我的经验是先从二维切片模型入手把网格、参数、边界和求解器全部摸清再上三维CT结构。盲算三维会浪费大量时间在排查上。最后建议每个算例都保留完整的动画和时间序列数据界面形态的动态变化往往比最终云图更能体现驱替机理做技术汇报和写分析报告时这些素材也会派上大用场。
干微观驱替模拟的人都知道真正影响结果走向的往往不是软件操作本身而是对物理过程的拆解和对数值稳定性的控制。我最近在COMSOL里做了一件事从岩石扫描图像中重建出真实的孔隙几何再模拟两种油相在这个孔隙结构里的扩散与驱替过程盯着前缘怎么运动、剩余油会在哪里滞留、不同润湿条件下驱替效率差多少。这套流程听起来不算宏大实际完整跑下来涉及的东西非常多从几何清理到相场参数标定任何一个环节掉链子结果都没法看。这篇稿子把整套过程完整写下来包括模型设计思路、关键参数设定、求解器选择、后处理分析以及我在调试过程中踩过的几个坑。适用对象是正在做多孔介质微观渗流、两相流数值模拟或者刚入坑COMSOL流体仿真想找参考的同行。下面直接进入正题。1. 内容整体设计与思路拆解1.1 “两油扩散与驱替”在物理上到底指什么把这个前提说清楚后面的参数才有讨论的意义。所谓两油其实就指两种物性不同的流体相最常见的是黏度、密度各不相同的两相油品也可以把体系换成油和水方法完全通用。我的项目里设定油相1为低黏度侵入流体油相2为高黏度被驱流体用来模拟注气或注低黏介质驱高黏原油的场景。从物理本质看驱替就是把一种流体“推”着取代另一种流体在孔隙中的位置。我们可以做个生活类比用吸管朝杯子里吹气气流把杯中果汁推向杯壁往外溢吹气的能量通过流体的压差传递到相界面上这就是一次微型驱替。发生在岩石孔隙里的驱替要复杂得多因为通道不是一根直管而是粗细交替、坑洼不平、带有死角的复杂网络。流体走到细喉道时界面会变形并产生毛细阻力过去之后又在孔隙里急剧膨胀前缘形状会变得非常不规则这也是微观驱替模拟比宏观渗流计算有趣的地方。标题里的“扩散”同样不是可有可无的字眼。COMSOL的相场法在求解界面运动时会天然包含扩散机制更直白的场景是在驱替基本结束后关闭入口流动让相场继续演化一段时间观察滞留油是否因分子扩散产生缓慢的组分迁移。这一步能把“驱替”和“扩散”两件事结合起来体现多场耦合的价值。1.2 为什么选COMSOL而不是OpenFOAM或其他工具我在做方案选型时纠结过一阵。开源CFD工具里的多相流求解器也能做这类模拟功能并不弱但COMSOL在这个场景下有几点优势非常明显参数化建模方便。材料属性、边界条件、控制方程耦合都在同一个界面里管理改一个变量全局自动更新相场接口成熟。自带的“两相流相场”物理接口集成度高不需要自己去拼代码里的方程组后处理灵活。积分算子、派生值、动画导出全都在一个工作区完成不用在求解器和后处理软件之间反复倒数据。代价是求解速度会比专门的多相流求解器慢一些。但微观孔隙模型规模通常不大百万级网格在单机上是可以接受的这种便捷性换一点计算时间非常划算。1.3 为什么非要“实际孔隙”简化模型不行吗很多初学者习惯用平行毛管束或单根变径管代替岩石孔隙算出漂亮的曲线后直接下结论。但这类简化模型有一个致命问题它会把真实结构最核心的特征也就是孔喉比、迂曲度和盲端孔隙全部抹平。比如一根直管里只要压差存在被驱流体几乎可以被推干净可在真实孔隙里流体沿大孔隙走捷径把小孔隙里的油绕过甚至反向夹住形成非连续的残余油。这种残余油靠简单几何是算不出来的。我做过一个对比相同边界条件下直管模型的驱替效率超过90%而同一参数放进真实CT重建的孔隙模型里驱替效率只有65%左右。这个差距足以说明为什么“实际孔隙”四个字是整个课题的基石。这也是项目名称里特别强调“实际孔隙”的根本原因——不还原结构就谈不上研究真实驱替机理。2. 几何模型从切片到可计算域2.1 孔隙结构从哪里来CT扫描与数值重建两条路要还原实际孔隙第一步是拿到孔隙空间的几何表达。主流途径有两条第一是CT扫描重建。把岩心柱放入显微CT设备得到一组二维切片序列经过滤波、阈值分割、二值化等图像处理后把岩石骨架和孔隙空间分离开。孔隙空间最终导出为STL或面网格文件COMSOL直接导入即可。CT扫描的优点是结构真实缺点是需要设备而且扫描分辨率直接决定可以分辨的最小孔喉尺寸。第二是数值重建。如果没有扫描数据可以用颗粒堆积法、四参数随机生长法等算法生成和真实岩心统计特征相近的数字岩心。这种方法生成的结构不是真实岩心但能反映孔隙尺寸分布、配位数等关键统计特征对于机理研究完全够用而且不用处理噪声数据入坑速度更快。我这次优先用了CT扫描数据主要原因是需要在真实结构里验证驱替过程。但为了前期调试物理场参数我先用数值重建生成了一块小尺寸的孔隙模型参数跑通后才把CT几何替换进来。这个策略节省了大量时间强烈建议参考。2.2 几何清理与导入别让STL的锯齿毁掉你的网格从图像里提取出的面网格通常有很多问题表面锯齿、破损面、相互穿透的面片、退化单元。这些问题如果不处理导入COMSOL后会造成几何修复失败或者后续网格质量极差。在COMSOL里做导入我的操作路径是在“几何”节点下选择“导入”载入STL文件然后开启“修复容差”参数。修复容差决定了多少范围内的缝隙会被自动缝合一般先取全局尺寸的1%左右试算看修复报告里的剩余破损面数量。如果破损面过多千万不要靠无限调大容差硬来那样会把细小孔喉直接填死结构失真。正确做法是回到图像处理阶段把STL重新做一次表面光滑和三角形简化。导入成功后做几件事删除退化面把面积过小、边长比例极端的表面片去掉布尔并集运算把多个独立导入的实体合并成一个整体域避免后续物理场装配的边界混乱检查连通域数量用“几何属性”查看孔隙域是否整体连通。实际岩心中本来就可能存在孤立孔隙但如果孤立孔隙太多到难以统计可能是阈值分割时把骨架和孔隙分错了。2.3 网格划分孔隙尺度模拟成败的第一道关口网格对相场模拟的影响比常规单相流要大得多。因为相场方法需要用一个有限厚度的过渡层去描述相界面这个厚度必须由网格解析。界面厚度参数ε一般取2到3个网格尺寸比如我设定ε5μm那么最大单元尺寸要控制在2.5μm以下。我最终的网格方案是在孔隙主体使用最大单元2.5μm的自由四面体同时在孔喉收缩段和预计界面扫过区域做局部加密壁面附近加12层边界层网格。整体网格量在百万级从孔隙尺度来看相当于用约40个网格分辨一个典型孔径。这个密度不算奢侈但已经很稳定。网格划完后一定要做一次网格无关性验证这也是我在项目里坚持做的事。方法是分别用最大单元尺寸1μm、2.5μm、5μm跑同样的算例对比驱替效率曲线。如果最粗和最细的结果差异在5%以内取中间尺寸就能满足计算精度要是差异明显就不能盲目信任当前网格。这一步看似费时实际上能帮你避开无数“假收敛、真误差”的尴尬。3. 控制方程与关键参数设定3.1 流体运动方程纳维-斯托克斯方程打底孔隙尺度两相流的流动本质由不可压缩纳维-斯托克斯方程控制再叠加相场方程去追踪界面。在COMSOL的“层流两相流相场”接口中这两组方程被耦合在一起求解。不可压缩NS方程的标准形式这里不写全式了关键要理解各物理项在孔隙尺度的角色。惯性力项在这种微米尺度流动中通常可以弱化因为雷诺数非常低。以我用的参数为例特征流速1mm/s特征孔径100μm油相黏度0.1 Pa·s密度850 kg/m³雷诺数约为0.001。在这样的量级下流动几乎完全由黏性力和毛细力主导前缘运动非常缓慢且平滑不会出现宏观流动里的湍流涡旋。这也决定了求解器的一个重要特征你不需要把时间步压到极小去追赶物理上的高速变化但会因为相场方程的强非线性存在数值稳定方面的挑战。3.2 Cahn-Hilliard相场方程界面厚度与迁移率的经典博弈相场方程比NS方程更容易让人翻车因为里面两个关键参数的取舍直接影响界面形态和计算稳定性。Cahn-Hilliard方程的基本形式可以理解为相场变量φ的对流输运加上一个由化学势梯度驱动的扩散项。在这个方程里表面张力σ和混合能密度λ之间存在换算关系σ (2√2/3) · λ / ε也就是说给定目标界面张力和界面厚度后混合能密度就确定了。反过来说如果随便填一个λ算出来的界面张力会和实验值差出好几倍。界面厚度ε不能太大否则界面过度模糊会把孔隙里本应清晰的前缘演变成一大片过渡区驱替特征完全失真也不能太小否则网格解析代价爆炸。实践中ε取1.5到3倍网格尺寸。迁移率γ同样关键。它控制着相场界面的数值扩散速度。γ太小质量守恒好但界面演化迟缓有时在复杂孔喉处会卡出假残留γ太大界面被快速抹开前缘锐度消失。我的经验是从1e-8 m³·s/kg起步试算如果界面前缘出现细碎“毛刺”小幅降低γ如果计算时间步长持续萎缩且界面不前进说明乘积项不合适再逐步回升。3.3 参数表一个可以直接抄作业的案例下面这组参数是我在项目中实际使用过一轮的基准值单位全部统一为国际单位制。你可以直接用这套数据做复现再根据自己的体系调整。参数名 表达式 描述 D 100[um] 特征孔隙直径 U_in 0.001[m/s] 入口速度 mu1 0.1[Pa*s] 油相1动力黏度高黏 mu2 0.001[Pa*s] 油相2动力黏度低黏 rho1 850[kg/m^3] 油相1密度 rho2 1000[kg/m^3] 油相2密度 sigma 0.03[N/m] 界面张力 epsilon 5[um] 相场界面厚度 gamma 1e-8[m^3*s/kg] 相场迁移率 theta_w 60[deg] 壁面接触角这组数据对应的毛细数为Ca μU/σ代入后约为0.0033刚好落在毛细驱替向黏性驱替过渡的区间。如果你希望不同驱替模式对比调整入口速度或黏度即可让Ca数量级发生变化。3.4 边界条件与初始条件最容易被忽视的物理陷阱几何域确定后边界条件设置直接决定流动模式是否正确。我的设置方式如下入口边界给定速度U_in同时指定相场变量φ1代表油相1进入出口边界给定压力为0并开启“抑制回流”选项防止数值回流把出口重新污染壁面边界设置无滑移速度条件同时设定接触角θ_w这是模拟润湿性对驱替形态影响的关键参数初始条件整个孔隙域初始化为φ-1代表初始完全充满油相2入口附近做一段缓冲让油相1可以从边界持续进入。这里有个经验必须单独强调不要让入口速度从0直接跳到目标值。比如U_in1mm/s的入口条件直接加载到一个充满高黏流体的静止域上压力场会瞬间冲击界面导致第一波计算步崩溃。我通常加一个时间相关的过渡函数让入口速度在0.01s内线性爬升到目标值。这个过程能显著减少初始震荡而且对最终驱替结果没有影响因为0.01s相比总模拟时间0.2s只占一个零头。4. 求解设置与后处理分析4.1 瞬态求解器全耦合还是分离式COMSOL里多相流问题一般选择瞬态研究。求解器有两种路线全耦合和分离式。我试过全耦合对小模型收敛很稳但在网格较大或时间步较长的阶段雅可比矩阵规模大会让单步时间显著拉长。后来我改用分离式求解器按“速度场-压力场”和“相场变量”分组迭代更新。分离式每次只解一个子块内存压力小而且可以通过调整子步迭代次数来稳定界面演变。分离式需要多迭代几轮才能达到同样的收敛残差但每个迭代步又快又稳整体算下来通常比全耦合省时。时间步长设置为自适应初始步长1e-4s最大步长0.001s后续交给软件自动调整。Courant数是一个重要的稳定性参考量它表示一个时间步内流体前缘移动的距离与网格尺寸的比值通常希望它保持在1附近。用我的参数算U_in1mm/s最大时间步0.001s前缘单步移动1μm与网格尺寸2.5μm同量级因此自适应时间步不会触发稳定性崩溃这点在做初始设置时可以自己估算。4.2 收敛判断与质量守恒检查判断计算是否真的收敛不能只看求解器报告“已收敛”。我最常做的检查有两项第一项是监察残差曲线。COMSOL求解器面板里可以勾选残差历史横坐标迭代步数、纵坐标残差范数。合格曲线应当持续下降或震荡下降至收敛容差以下而不是反弹后停在某个值附近不动。如果残差在1e-2附近长时间徘徊多半是界面相关条件设置出了问题。第二项是质量守恒。在域内积分相场变量或者更简单地对比入口与出口的体积流量差。如果累计流量差超过入口总流量的2%就要回头检查收敛和网格。我遇到过看起来完全收敛、但出口体积流量和入口差5%的算例最后发现是出口边界上的“抑制回流”选项没有打开数值回流导致假“增油”。这类问题大多可以靠边界条件修复而非调小容差硬撑。4.3 看动画比看曲线更能发现问题这是最想分享的经验之一。定量曲线当然重要但在模型调试阶段界面运动的动画才是真正帮你发现问题的手段。把几组时间截面的相场云图按时间序列播放前缘是否稳定、是否产生指进、油相1是否在孔喉处被截断一眼就能看出来。我在项目里导出了每隔0.005s一帧的相场分布图然后合成动画。第一次播放就发现被驱替流体在一个大孔隙的盲端被反复夹断形成孤立油滴。这个现象在饱和度曲线里只能看到平台期但动画里看机理一目了然。强烈建议每一组参数组合都保存动画很多后续论文或汇报用的机理分析素材就是这么来的。4.4 定量提取指标饱和度、前缘位置与残余油占比后处理阶段不能只看云图需要提取可比较的量化指标。平均饱和度在“派生值”中用积分算子对相场变量φ做域积分再除以区域体积。因为φ从-1到1表示不同相可以换算成油相1的体积分数前缘位置沿主流方向取相场φ0的等值面位置记录每时刻的推进距离残余油占比模拟结束后统计φ0区域中仍构成连通体的体积占比。配合“连通域分析”可以区分哪些残余油是盲端孤立油哪些是可动但未驱走的油。这些指标做出来后不同接触角、不同毛细数之间的对比就有依据了。用表格把静默数据变成趋势描述也方便后续项目总结。5. 算例分析润湿性、毛细数与扩散效应的面面观5.1 不同润湿条件下的驱替差异接触角是润湿性最直接的量化指标。我分别算了θ_w30°、60°、90°、120°四组。接触角30°代表强亲水界面张力对壁面有很强的“贴壁”倾向流体沿壁面铺展更快前缘推进相对均匀总体驱替效率最高。接触角120°后孔隙几何里很多曲面位置出现界面钉扎油相1的前缘推进非常不均衡产生明显的指进形态驱替效率下降显著。这组对比比较直观在不同方向工作时为润湿性调控方向的选择提供了依据。5.2 不同毛细数下的驱替模式转变毛细数的改变可以通过调整入口速度实现本质是在改变黏性力和毛细力的比值。我做了从0.001到0.1的量级扫描。低毛细数区界面张力占主导孔隙孔喉对前缘的束缚作用强前缘容易在细喉处被“锁住”驱替表现为非连续的跳跃式推进高毛细数区黏性力压制界面张力效应前缘推进更连续、更均匀但容易在非优势通道形成纵向指进。在孔隙尺度模拟中重现这种驱替模式的转变是从定性描述走向定量分析的关键一步。后续可以进一步绘制驱替效率随毛细数变化的散点图找到本模型中的过渡区间。5.3 驱替结束后的扩散驰豫不可忽视的“二次效应”驱替主过程结束后我并没有直接终止计算而是把入口速度降为0继续演化0.5s的物理时间观察相场继续扩散驰豫的状态。这时流体不再整体推进但相界面的梯度驱动扩散仍然存在高黏流体在盲端中的组分分布会发生缓慢重分布。这个算例的发现是一部分在驱替过程中被“夹断”的孤立油滴相邻的浅色区域主流通道路内浓度梯度持续存在经过扩散演化油滴界面发生了轻微移动部分死角中的油相浓度被稀释。虽然这个过程在短时间内对采出程度贡献有限但对于长时间尺度的采收率分析却有参考价值。这也是将“扩散”纳入项目标题的意义所在——区块的采收率不仅取决于驱替压差和界面性质还取决于组分在孔隙内的扩散能力。6. 常见问题与排查技巧实录6.1 界面前缘不光滑出现细碎“毛刺”这种情况在初期调参时几乎都会遇到。表面看是网格不够细仔细排查后才发现是迁移率γ取值偏大导致相场界面临近振荡。解决办法将γ降低半个数量级从1e-8降到3e-9然后再看是否出现毛刺。如果降低后前缘变得平滑但推进速度变慢说明γ仍处在合理范围内只是还没有收敛到我需要的形态。关键是要找到形态稳定与扩散速度之间的平衡点。6.2 时间步长越算越小最后算不下去这是相场模型最头痛的故障。出现这种问题先看提示中是否出现界面尖角或压力剧烈变化。常见触发原因是入口速度阶跃启动造成前缘位置的压力尖峰。对策是前面提到的入口速度爬升函数把0.01s的线性过渡加到所有初始参数上。另外要留意界面厚度ε与网格尺寸的比值如果ε小于网格尺寸约1.5倍相场方程对网格的灵敏度会急剧上升时间步也随之萎缩。此时稍微增大ε或加密局部网格问题往往会缓解。6.3 出口回流导致饱和度曲线异常前面提过出口边界如果不开启“抑制回流”求解器会在出口位置产生轻微回流以满足压力条件。这种回流会让出口处的相场波动导致累计流量统计不准确。排查方法很简单绘制出口边界上的速度分布如果边界某区域速度方向为负那就是回流。开启抑制回流选项后大多数情况都能解决。如果开了抑制回流仍异常需要检查出口是否紧邻一个特别狭窄的孔喉段可以在出口与孔喉之间人为加一段延长区。6.4 剩余油驱不走是数值问题还是物理现象模拟运行后把油相2“锁”在某个盲端第一反应经常是“网格不够、参数不对”。我的判断方法是做三个交叉验证一是把最大网格尺寸缩小一倍重新算二是把迁移率改变一个数量级重跑三是把入口速度提高10%再试。如果三种改变下残余油仍然稳定存在于同一位置那基本可以认定这是真实的物理现象是孔隙形态和润湿条件共同决定的残余油模式。如果位置随网格或参数漂移则是数值伪影需要回到网格与参数设置继续调。这个方法帮我区分了大量假问题也让后续分析更有底气。6.5 避坑清单速查入口速度一定要做爬坡过渡不要阶跃加载ε与网格尺寸比值保持1.5倍以上界面才不会被网格离散破坏迁移率过大是界面毛刺的头号元凶出现毛刺先降γ每次参数变更前先保存一份当前算例多参数对比时命名按“接触角_毛细数”的格式避免后期数据混乱求导入几何时修复合并不等于最终几何是正确的一定要用“几何属性”检查域体积是否合理饱和度的积分计算要选对相场变量符号φ1与φ-1表示两相符号反向结果完全不同动画输出分辨率不需要太高云图加流线叠加的模式最实用兼顾信息量和文件大小。我个人这段周期的体验是COMSOL跑这类孔隙尺度两相流并不难难的是永远知道当前结果是在反映真实物理还是放大了数值噪声。要掌握这个判断力我的经验是先从二维切片模型入手把网格、参数、边界和求解器全部摸清再上三维CT结构。盲算三维会浪费大量时间在排查上。最后建议每个算例都保留完整的动画和时间序列数据界面形态的动态变化往往比最终云图更能体现驱替机理做技术汇报和写分析报告时这些素材也会派上大用场。
