COMSOL相场方法模拟多孔介质两相驱替：毛管数与饱和度计算

两相驱替这个词在不同行业里叫法不一样油气行业叫“水驱油”地下水修复那边叫“非水相液体迁移”燃料电池方向则习惯叫“水气管理”。但说到底都是同一件事——一种流体把另一种不混溶流体从多孔介质里挤出去。我这次要分享的项目就是用COMSOL的相场方法把多孔介质里的两相驱替过程完整模拟出来同时把毛管数、饱和度这两个最常用的评价指标的计算方法也一并整理好并且做了一些参数拓展研究。这套东西很适合正在做渗流模拟、微流控实验对照、油藏数值模拟和CCUS相关研究的人参考。COMSOL里跑两相流物理上最麻烦的就是那个油水界面。界面在孔隙里会不断变形、断裂、合并传统的网格追踪方法在这种拓扑变化面前非常吃力。相场方法好就好在它不直接追踪界面而是用一个连续的序参数场把界面隐式地“糊”在变量场里界面变化自然就处理掉了。这也正是我为什么最终选择相场而不是水平集或者VOF的原因。下面我把整个项目从思路、建模、算例到拓展和避坑完整拆开讲。1. 项目思路与核心物理背景1.1 多孔介质两相驱替的应用场景多孔介质两相驱替的应用范围其实比很多人想象的要广。最常见的自然是油藏水驱开发水注入地层把原油推向生产井这个过程是否高效直接决定了最终采收率能到多少。再比如地下水污染修复把表面活性剂溶液注入含水层去冲刷残留在孔隙里的有机污染物本质上也是驱替过程。还有燃料电池气体扩散层里的水管理——反应产生的水必须被及时排出否则会堵住气体通道这也是典型的两相驱替问题。在这些场景里驱替效果好坏的核心评价指标有两个一个是驱替效率也就是被驱替相最终能采出多少通常用饱和度变化来衡量另一个是驱替过程中的形态特征流体会不会出现指进、绕流、卡断这些微观现象直接影响宏观采出程度。这两个指标看起来简单实际模拟起来牵扯的东西却非常多。1.2 相场方法如何解决界面追踪问题界面追踪在数值模拟里一直是个难点。早期的界面追踪类方法比如VOF和水平集都需要专门维护一个界面位置。VOF要靠几何重构来保持界面锐利水平集要定期重新初始化符号距离函数。在孔隙尺度这种复杂拓扑环境里界面频繁合并和断裂这两种方法的维护成本很高而且在分支点容易出现数值上的质量问题。相场方法的思路完全不一样。它引入一个相场变量比如在COMSOL的层流两相流相场接口中变量 φ 的物理意义可以理解成“归一化的流体浓度场”。φ 1 表示一种流体φ -1 表示另一种流体两者之间用一个很薄的过渡层连续连接起来这个过渡层就是“扩散界面”。界面处的表面张力效应被转化成一个分布力施加在过渡层上这样就不需要显式追踪界面位置了。扩散界面的做法带来的直接好处是拓扑变化随便发生。界面想分叉就分叉想合并就合并数值格式不需要做任何特殊处理。而且COMSOL把相场方法和Navier-Stokes方程耦合得很成熟拿来就能用。缺点也有后面我会专门讲比如界面厚度这个参数要取合适值网格分辨率要跟上计算量比纯VOF可能更大。1.3 毛管数与饱和度从微观形态到宏观评价做完驱替模拟如果你只是盯着动画看界面怎么弯弯曲曲地走那是远远不够的。你得有定量指标来比较不同工况的优劣。毛管数和饱和度就是这个项目里我认为最核心的两个定量指标。毛管数定义为粘滞力与毛管力的比值宏观定义一般是 Ca μU/σ其中 μ 是流体黏度U 是特征流速σ 是两相界面张力。它控制的是驱替前沿的稳定性Ca 很低的时候毛管力占主导驱替过程呈毛细指进特征界面推进比较慢容易在细孔隙处发生卡断Ca 升高到一定程度粘滞力占主导驱替前沿变得更稳定波及效率可以提高但Ca过高又可能在局部形成粘性指进。饱和度则直接给出定量结果湿相饱和度 Sw 表示湿相占据的孔隙体积比例残余油饱和度 S_or 是评价驱替效率最常用的指标。把这两个参数量出来驱替过程才能从“看动画”变成“比数字”。所以这个项目的整体逻辑就是先用相场方法把孔隙尺度的两相驱替过程复现出来再通过毛管数和饱和度这两个指标对结果做定量评估最后利用参数扫描探索不同工况下的驱替规律。下面我按模型搭建、求解设置、算例实操、拓展研究的顺序展开。2. 多孔介质几何与相场模型搭建2.1 几何建模随机球堆积和网络模型怎么选多孔介质的几何构建方式直接决定模拟结果是否可靠这个环节一定要花时间认真做。我在这个项目里用过的几何生成方案主要有三种各有各的适用场景。第一种是随机球堆积。在COMSOL里先用“随机分布”功能生成一系列圆去除掉相互重叠的部分然后布在矩形模拟域内。这种方式的优点是孔隙结构具有一定随机性接近真实岩心的拓扑特征适合做研究性算例。孔隙率可以通过控制圆的半径和数量来粗调比如目标孔隙率0.4左右时通常可以随机放置半径在10到30微米之间的圆形颗粒重叠处理后再通过布尔运算得到孔隙域。第二种是规则微通道网络。用周期性排列的柱体或通道交叉形成网络结构高度规则重复性好适合做物理机理研究。这类几何更容易获得干净的网格计算效率高而且便于验证结果比如饱和度曲线能跟解析解或者文献数据对比。第三种是把高精度CT扫描图像导入COMSOL。这是最接近真实多孔介质的方式但建模工作量大需要做图像分割、生成曲面网格计算成本也会成倍增加。如果项目刚刚起步我建议先从第一和第二种方式做起先把相场参数和求解策略摸透再考虑真实岩心图像。无论哪种方式都需要在建模时计算孔隙体积因为后处理里算饱和度要用。COMSOL提供了体积积分算子建立好几何后可以先积分一次得到孔隙体积作为基准值。2.2 COMSOL物理场与控制方程两相驱替的物理场搭建我用的是COMSOL内置的“层流两相流”物理场接口中的相场分支。这个接口本质上是由两个方程耦合组成的一个是Navier-Stokes方程描述流体的动量守恒另一个是Cahn-Hilliard方程描述相场变量 φ 的输运。Cahn-Hilliard方程可以拆解成两个子方程来理解。第一个是相场变量的对流扩散方程∂φ/∂t u·∇φ ∇·(M∇η)其中 η 是相场辅助变量M 是迁移率。第二个子方程给出 η 与 φ 的关系η λ(-∇²φ φ(φ²-1)/ε²)其中 λ 是混合能密度ε 是界面厚度控制参数。这个方程组把界面的扩散和移动统一描述在一个框架下。表面张力通过体积力项施加在界面上COMSOL里会将其转化为包含相场梯度的体积力表达式。实际使用中你不需要手工推导这些方程COMSOL已经把参数化工作做好了。核心参数包括界面厚度 ε控制扩散界面的宽度需要远小于孔隙特征尺寸通常取孔隙直径的1/20到1/50迁移率 M控制界面处相场的扩散速度影响数值稳定性COMSOL里有默认值建议但要根据特征时间调整表面张力 σ直接从物性参数设定接触角 θ在三相交界处设定决定壁面的润湿性。需要特别注意界面厚度和迁移率是一对“数值敏感性参数”它们不像物性参数那样可以直接从手册上查必须根据网格尺寸和特征时间合理估算。我一般会把界面厚度设置为网格特征尺寸的2到3倍这样既能保证界面处的连续性也能避免需要极端细网格。2.3 边界条件与润湿性设置边界条件牵涉到入口、出口和壁面三类。入口边界我用的是定速度入口给定一个驱替流速速度方向与流动方向一致。出口边界设为定压通常是0或者某个参考压力。壁面设为无滑移壁面并且在壁面边界上额外设置润湿接触角。接触角在相场方法里的处理是润湿边界条件。COMSOL的“润湿壁”边界可以指定接触角 θ原理是让相场变量的法向梯度满足一个与接触角相关的边界条件。这个设置对指进形态影响非常大我在早期测试时跳过这个设置结果所有工况都变成了中性润湿驱替形态跟真实情况差得很远。所以这个边界条件一定不能漏。初始条件方面通常在模拟域内先填充被驱替相比如油让相场变量初始值在大部分区域为 -1入口附近设置为 1表示注入水速度场初始为0。压力场初始为0或者静水压分布即可。3. 数值参数与求解策略3.1 特征时间与时间步长选择相场两相流模拟中选择合适的时间步长是项目成败的决定性因素之一。COMSOL虽然支持自适应步长但自适应算法在很多复杂孔隙几何下并不够稳我更倾向于自己设定一个合理的时间步长范围。时间步长的估算基准是特征时间。假设模拟域长度为 L 500 μm注入速度 U 1 mm/s那么对流特征时间 t_c L/U 0.5 s。如果你对整个过程需要分辨100个时间点那么初始时间步长取 t_c/1000 0.0005 s 是比较合理的起点。另外一个必须考虑的是相场方程本身的扩散时间尺度。相场界面厚度 ε 取1 μm时相场扩散特征时间约为 ε²/MM 取1×10^-10 m²/s量级的话这个时间是微秒级导致显式时间推进很受限制。实际计算中COMSOL的隐式求解器允许较大的时间步长但需要监测界面位置处是否出现数值振荡一旦振荡就要减小步长。我自己习惯的做法是先把时间步长设得很小跑一小段比如前0.01秒观察界面是否平滑移动然后逐步放大步长到对流时间尺度的百分之一左右。这样可以兼顾计算速度和数值稳定。3.2 网格加密和界面厚度关系网格是相场模拟里最容易踩坑的地方。相场方法在界面附近有一个厚度为 2ε 的过渡层这个过渡层必须至少覆盖两到三个网格单元否则界面处的相场梯度无法被准确计算会导致表面张力计算失真甚至产生虚假的“液滴分裂”现象。网格策略大致可以这样定。先确定界面厚度 ε比如取 ε 5 μm那么界面附近的网格尺寸需要控制在 ε/2 ≈ 2.5 μm 左右。孔隙直径如果是 30 μm那么一个孔隙直径方向大约需要12个网格对于二维模拟来说这就是数百倍甚至上千倍于最小网格的网格数量。整个计算域通常要布置5万到20万个三角形单元网格数多少取决于模型大小和精度要求。COMSOL里可以使用“自适应网格细化”功能让细化集中在相场界面运动轨迹可能经过的区域。但需要注意相场界面不是静止的它会在模拟过程中移动所以自适应网格必须能跟随界面移动否则你会看到界面冲进粗网格区后形态明显异常。如果网格实在加密不起来一个折中方案是把界面厚度 ε 适当调大但必须在合理范围过大会导致界面弯曲区域被抹平。3.3 收敛性调整与常用求解器设置COMSOL的默认求解器设置对大多数简单流动问题足够但加上相场耦合以后默认配置经常会出现“收敛失败”或者“时间步长持续减小”的情况。我的经验是把求解器设置从默认配置改为“拟牛顿”或“阻尼牛顿”迭代并开启“松弛因子”来稳定非线性迭代。具体操作上我通常用分离式求解器先解Navier-Stokes方程更新速度场和压力场再解相场输运方程更新 φ两个子问题交替迭代。这样做的好处是内存占用低、收敛相对稳定。缺点是需要设置好松弛系数速度场的松弛系数一般设在0.7左右相场的松弛系数可以低一些比如0.5尤其在接触角接近90度时更需要放松相场迭代。还有一个很实用的小技巧先给一个很低的注入速度跑通整个模型得到稳定的驱替流场作为初始化再恢复到目标速度继续计算。这种“分阶段加载”的做法能有效避免初始时刻流场突变导致的求解失败处理复杂几何时尤其管用。另外COMSOL的相场接口里有一个“初始所有变量解耦”的选项开启后能让初始时刻的迭代更稳定但也会稍微增加前期迭代次数对严重非线性问题很有效。4. 算例实操与毛管数、饱和度计算4.1 算例设计模型参数速查直接上一个完整的算例参数表方便后面复现。这里我用了随机球堆积几何模拟域尺寸取1000 μm × 500 μm孔隙率控制在0.38左右。两相流体设定为水和油具体参数如下参数数值说明模拟域尺寸1000 μm × 500 μm二维颗粒半径范围15~30 μm随机分布孔隙率0.38通过颗粒覆盖面积估算水相黏度1 mPa·s注入相油相黏度10 mPa·s被驱替相水相密度1000 kg/m³油相密度800 kg/m³界面张力0.03 N/m油水界面接触角45°水湿注入速度1 mm/s入口恒定界面厚度 ε3 μm约为孔隙直径的1/10迁移率 M2×10^-10 m²/s经验值出口压力0 Pa参考压力利用这个算例计算全局毛管数Ca μ_非湿相 × U / σ。这里被驱替相是油黏度10 mPa·s速度1 mm/s界面张力0.03 N/m得到 Ca ≈ 10 × 0.001 × 0.001 / 0.03 ≈ 3.3×10^-4。这个数量级处在粘性指进和毛细指进的过渡区间模拟结果会非常容易看到指进结构很适合做相场方法的演示算例。4.2 驱替过程后处理饱和度场与界面前沿跟踪算例跑完后最直观的结果就是相场变量 φ 的分布图。把 φ 用颜色标出来正负分别对应水和油中间过渡层就是界面位置。在COMSOL结果模块里可以制作动画查看驱替前沿随时间推进的过程。你会看到水相从入口进入孔隙网络遇到颗粒阻碍时绕过颗粒流动在局部窄孔处出现卡断形成被困住的油滴。跟踪界面前沿我常用的是等值线图也就是 φ 0 的等值线。直接在结果模块里绘制“φ0”的等值线可以看清界面形状随时间的演化。如果配合速度场云图还能分析出高流速通道和死区的位置关系这对接下来的毛管数分析很有帮助。4.3 毛管数计算从达西速度到局部毛管数毛管数的计算看似简单实际操作中有几个容易出错的点。全局毛管数用的是达西尺度特征速度一般取入口速度这个积分平均意义下的毛管数只能反映驱替过程的整体操作条件并不能描述局部的指进行为。因为孔隙网络内部流速分布不均匀窄孔处的局部速度可能是入口速度的数十倍局部毛管数比全局毛管数大得多。局部毛管数的定义需要考虑该处的特征剪切速率。更严格的说法是使用局部速度梯度和黏度的乘积替代 μU即 Ca_local μ × |∇u| × d_p / σ其中 d_p 是局部孔隙直径。我在COMSOL里计算时会先在孔隙中心线上布置若干探针点输出这些点处的速度、速度梯度和压力数据导入到数据处理软件里计算出局部毛管数分布。值得注意的是对于粘性指进现象使用全局毛管数往往低估了局部不稳定性。我在模拟结果中看到过全局毛管数只有3×10^-4但局部毛管数已经接近10^-2这直接导致界面出现了明显的分叉指进结构。所以在结果分析时一定要分辨清楚讨论的是全局毛管数还是局部毛管数否则很容易得出互相矛盾的结论。4.4 饱和度计算体积积分与无因次曲线饱和度是指湿相体积占孔隙体积的比例。在COMSOL的后处理中可以直接定义变量 Sw (φ 1)/2它表示水相体积分数取值范围从0到1。然后使用体积积分算子对整个模拟域积分 Sw除以孔隙体积就得到当前的饱和度。这里有个容易忽略的细节积分域必须限制在孔隙区域而不是整个模拟域。因为模拟域里还有固体颗粒区域如果对整个矩形域积分结果会把孔隙率混进去。正确做法是在“定义”里创建一个“积分算子”积分域选择“孔隙域”组件几何实体这样后续每次评估饱和度都直接用这个算子。随着驱替过程推进饱和度会逐渐上升。把饱和度对无因次时间画出来就是经典的驱替曲线。无因次时间可以定义为累积注入体积除以孔隙体积。这个曲线能告诉你驱替的早期快速推进段、中期的稳定段以及后期的低速段分别在哪里。在模拟结果中通常能看到驱替前期饱和度上升很快因为水相沿着主干通道快速推进后期出现指进和绕流后上升速度明显放缓。4.5 结果解读粘性指进与毛细指进判别最后要会判断模拟结果属于哪种驱替模式。我的判断依据主要看两个因素全局毛管数和接触角。在强水湿条件下如果全局毛管数较小比如低于10^-5量级驱替过程以毛细力主导流体倾向于优先进入大孔隙界面推进缓慢均匀形成类似“分层”的推进模式。如果全局毛管数升高超过临界值粘性力开始主导界面前沿会出现明显的指状结构。如果黏度比小于1即驱替相黏度低于被驱替相那么即使毛管数不高也容易出现粘性指进。拿上面的算例来看全局毛管数约3×10^-4已经是过渡区间偏粘性主导的位置模拟结果确实出现了明显的分支状指进。这种指进形态通常是入口速度偏快或者油水黏度比过大导致的。如果减小入口速度比如降到0.1 mm/sCa降到3×10^-5驱替过程就会过渡为更均匀的活塞式推进。这个对比做起来非常直观也适合拿来做敏感性分析。5. 拓展研究思路5.1 参数扫描Ca数、黏度比与润湿性的系统影响这个项目做完基础算例之后自然要往前走一步做拓展研究。最核心的拓展方向就是对影响驱替的参数做系统扫描找出各参数对毛管饱和度曲线和指进模式的定量影响。参数扫描可以选好几个维度。第一个是注入速度。把入口速度从0.05 mm/s到2 mm/s扫描一遍获得对应的毛管数从10^-5到10^-3量级范围观察驱替形态从毛细指进到粘性指进的转变。这里需要跑多个瞬态模拟每个算例的网格和时间步长基本保持一致保证对比的公平性。第二个维度是黏度比。固定毛管数不变把黏度比从0.1扫到10观察相同毛管数下不同黏度比对界面前沿稳定性的影响。这个参数其实比毛管数更直观地决定指进模式。第三个维度是润湿性。通过改变接触角从30°到120°可以考察水湿、中性润湿、油湿条件下的驱替差异。油湿条件下水相难以进入小孔隙残余油饱和度会显著升高这在地下CO2封存和油藏水驱方案设计里都有重要意义。参数扫描数据整理出来后可以画成驱替效率随毛管数变化的经典曲线。通常这个曲线是S形低Ca时驱替效率很低高Ca时驱替效率升高到某个平台。这个平台对应的饱和度就是不同工况下的残余油饱和度对油藏工程有直接参考价值。5.2 三维建模与随机实时生成多孔介质二维模拟是很好的起步但驱替过程毕竟是三维的。二维模型会严重低估孔隙网络中的连通性对卡断现象的模拟也可能跟实际不相符。所以更进一步的研究方向就是三维建模。三维模型跟前文提到的几何生成方式是一样的只是把圆换成球模拟域变成三维立方体。但计算代价会大幅上升。一个二维算例往往只要10分钟到半小时就能跑完同样的物理参数在三维中网格数可能增加两个数量级计算时间会从小时级别增长到几天级别。我实际做过的三维算例在COMSOL中使用了共轭梯度求解器和并行计算效率提升还是明显的。如果追求更真实的孔隙结构可以考虑用随机填充和图像重建方法生成三维网络模型。COMSOL支持从图像导入并识别孔喉结构这也是一条可行的路线。三维模拟的意义在于能更好地评估驱替效率和残余油饱和度的代表值这些结果更接近真实岩心实验的数据。5.3 多物理场耦合反应流与热-流耦合相场两相驱替模拟的另一个拓展方向是与其他物理场的耦合这恰好是COMSOL最具优势的地方。比如考虑驱替过程中的传热。界面张力和流体黏度都是温度的函数将温度场方程与流动方程耦合可以模拟热采或蒸汽驱过程。实现起来就是在流体传热接口中添加温度依赖的表达式把σ(T)和μ(T)写成温度插值函数。再比如考虑溶质输运与化学反应的耦合。在地下CO2封存或污染修复中驱替液与孔隙壁会发生化学反应改变孔隙结构并影响润湿性。可以把化学反应产物的浓度场与相场边界条件耦合让接触角随表面活性剂浓度变化这样就能模拟表面活性剂驱油的详细过程。这种耦合计算在纯理论推导中很难做但COMSOL的模块化接口让你能比较容易地把这些物理过程叠上去改起来可控性好。6. 常见问题与排查实录6.1 收敛困难与初始条件设置我在做这个项目的过程中遇到最多的问题就是求解器不收敛具体表现是时间步长一次一次缩小最后停在极小值报错。排查下来这类问题一般有三个原因。第一个原因是初始时刻流场突变。入口速度从零瞬间加到目标值压力场也要瞬间建立数值上很难适应。解决方法是分阶段加载或使用一个较短的斜坡函数来加载入口速度比如让入口速度在0.05秒内线性升高到目标值这样前期迭代就稳定得多。第二个原因是相场界面厚度与网格不匹配。如果ε取得太小相场过渡层里只有一层网格数值格式很难解析正确的曲率。我建议用试算确定固定ε逐步加密网格对比界面长度和饱和度的变化直到结果不再显著依赖于网格。第三个原因是迁移率过大数值上引起界面的非物理扰动。最简单的办法是降低迁移率一个量级观察结果是否变化如果基本不变就说明已经收敛。6.2 界面模糊与质量守恒检查很多新手做相场模拟跑完之后发现界面区特别宽整个孔隙都快被“糊”成一片这时候首先检查三件事ε是否过大、网格是否过粗、迁移率是否过高。ε过大的直接表现是界面宽度接近孔隙半径这时候表面张力效应会被“摊薄”界面运动形态失真。网格过粗则是分辨率跟不上界面附近无法准确计算相场梯度。注意这两者会互相影响界面厚度越薄需要的网格数越多计算量越大所以需要折中。质量守恒检查是每个算例跑完必做的一项。方法很简单统计总的水相体积φ积分随时间的演化理论上应该等于入口累计注入流量加上初始水相体积。如果体积明显不守恒先检查润湿壁边界条件是否启用了接触角边界对质量平衡有影响再检查迁移率是否偏大导致数值扩散最后检查求解器是否使用了过大的时间步长导致相场输运被数值抹平。这里我踩过最深的坑就是忘记计算时间步长上限把步长拉得太快结果饱和度看似上升很快实际上全是数值扩散的假象。6.3 毛管数和饱和度计算的坑毛管数计算常见的坑主要有两个。一个是粘度取错了相。全局毛管数在油藏工程里通常使用被驱替相的黏度也有文献用驱替相的黏度两者数值可能差十倍以上。建议在报告里明确写明采用哪种定义否则对比文献时极易混淆。另一个是速度取值问题。如果模拟结果中入口速度不是恒定的比如用压力驱动毛管数计算就不能直接取入口速度而应该取整个孔隙网络内速度的达西平均。我一般用体积平均值也就是用总流量除以孔隙域截面积这样算出来的毛管数才和达西尺度实验结果可比。饱和度计算容易踩的坑是积分域和权重。前面已经说过要排除固体区域另外不能用算术平均替代体积积分。如果直接用相场变量在全部网格上的平均值再除孔隙率在高非均质孔隙结构下会产生不小的偏差。务必用积分算子按体积加权计算。最后再分享一个我做了很久之后才总结出来的小技巧在COMSOL里可以建一个“全局评估”节点把饱和度、注入体积、毛管数都定义成全局变量每次算完一个算例后自动输出一张画好的驱替曲线这样参数扫描跑几十个工况时能轻松对比。如果你也想拿相场方法做两相驱替我建议先把本文的算例完整跑通一遍把饱和度曲线和毛管数都算出来再往三维和多物理场方向拓展。那段从“能跑出来”到“能解释清楚”的过程比实现本身更有价值。