水驱油这几个字在我们做油气田开发模拟的人眼里其实就是把地层里那些靠天然能量已经采不出来的原油用注入水一点一点顶出来。整个过程听起来简单真正落到数值模拟层面就麻烦不少水进到孔隙里不是均匀推进的而是会沿着高渗条带形成指进油水两相之间还有一个模糊的过渡带再加上毛细压力和相对渗透率的非线性关系这些物理机制叠在一起用单一连续介质方程根本描述不清楚。所以我把目光放在了Comsol Multiphysics上用多孔介质两相流框架去还原这个水驱油过程。这篇内容适合两类人看一类是刚开始接触多孔介质多相流、被一堆饱和度方程和相对渗透率模型劝退的研究生另一类是在做储层/岩土/环境渗流相关项目、想快速搭一个可复现模型来验证自己方案的工程师。这里我不谈教科书上那种从零推导理论的过程而是直接告诉你Comsol里做水驱油模拟需要关注哪几个物理量、怎么设置相场和Darcy接口、最容易出问题的网格和求解器配置在哪以及我踩过的那些坑是怎么填平的。整个模型我从搭建到出结果一共花了两天半中间删掉重建了好几次最后跑通的方案比最初设想的要稳很多。下面从头拆解。1. 水驱油模拟的核心物理图景与Comsol方案选型1.1 多孔介质多相流到底在算什么先捋一下我们面对的对象。多孔介质粗略理解就是一块海绵或一块砂岩里面有大量的微小孔隙和喉道流体只能在连通的空间里流动。原油刚开始是均匀分布在孔隙里的注入井以恒定压力或恒定流速往里面打水水沿着孔隙通道推进把油从生产井那端逼出来。这里面有两个尺度的问题值得先说清楚。在我做的这个模型里采用的是连续介质尺度也就是把岩石骨架中的每个数学单元体看作包含孔隙和骨架的“等效体”用孔隙度、渗透率、饱和度来宏观描述而不是真的把某一条孔隙通道用三维几何画出来。这样做的好处是计算量可以接受而且工程上最关心的是油、水两个相的推进前缘和压力剖面而不是单个孔隙内部的流动细节。多相流的核心方程是基于Darcy定律扩展出来的两相流动方程油相和水相各写一个动量方程和一个质量守恒方程。每个相都有一个独立的压力场而两相之间的压差就是毛细压力它是饱和度的函数。再加上相对渗透率它也随饱和度变化。这几个变量之间互相耦合、互相制约所以求解过程本质上是一个强非线性问题。用大白话描述就是水往前跑的时候孔隙里含水饱和度升高水的相对渗透率变大油必须让出通道反过来某个位置油的饱和度越高油越容易流动水就越难往那边挤。这个“你争我让”的过程就是水驱油模拟的底层逻辑。1.2 为什么选Comsol而不是CFD工具有人可能会问水驱油用Fluent、OpenFOAM这类CFD工具不是也很常规吗确实可以但我选择Comsol的原因有三个。第一多相流Darcy模型在Comsol中被封装成了专门的物理接口不需要像在CFD工具里那样自己耦合多孔介质动量源项和相输运方程。Comsol里有现成的“多孔介质多相流”多物理场耦合底层帮你把Darcy定律、饱和度输运方程、毛细压力模型组合到一起。第二Comsol对各种材料属性和自定义方程支持非常方便。比如相对渗透率曲线你可以用内置的Brooks-Corey或van Genuchten模型直接填参数也可以自己写一个分段函数插入进来。这对于模拟不同岩心实验数据的情况非常友好。第三Comsol后处理高度集成尤其是在观察油水前缘推进、动态调整监测点时比把计算结果导出到其他软件里再画图省事很多。我后来还顺手用它做了Desktop高密度3D打印样件的热-力耦合扩展阅读发现同一套多物理场框架确实能覆盖很多不同的物理过程手感熟悉之后跨度很大也不慌。1.3 物理接口与算法选择两相流Darcy、相场、水平集Comsol里处理多孔介质多相流主要有三条路两相流Darcy接口、相场接口、水平集接口。不少人第一次进入“模型向导”时看到一堆接口名字会懵我用过之后给你一个直接的选择建议。如果模拟的目标是油水两相在孔隙介质中的宏观推进优先考虑“多孔介质多相流”下的两相流Darcy定律接口。这个接口是专门为饱和度方程设计的计算的是Darcy速度场和压力场同时输运水相饱和度。它不需要显式捕捉油水界面适合地层尺度和岩心尺度。如果是想把孔隙尺度上的油水界面边界层效应也模拟出来或者你的模型几何本身是孔隙级微通道那就用“两相流相场”或“两相流水平集”。这两个都是基于界面捕获方法计算量明显更大但能够得到更精细的界面形态适用于数字岩心、微流控芯片这类场景。我做的是长约1 m、高约0.2 m的二维模型属于典型的宏观岩心尺度所以果断选了两相流Darcy接口加Brooks-Corey相对渗透率模型。后面的所有设置都围绕这个方案来展开。2. 模型搭建全流程从几何到边界条件2.1 几何建模与材料参数怎么填打开Comsol后建议直接选择二维模型。先在最左侧的组件节点下创建几何用矩形画一个长1 m、高0.2 m的岩心区域。不需要画注入井孔洞因为在这个尺度上井筒细节没有意义直接定义边界条件就能代表注采过程这也是连续介质模型的常见处理方式。材料块方面我在多孔介质子节点下给整个矩形域赋予同一组多孔介质属性主要包括孔隙度、渗透率和有效扩散张量如果是等温不涉及组分扩散扩散项可以不管。我用的参数组合是孔隙度0.35渗透率100 mD。这里要解释一下单位100 mD换算成国际单位约为1e-13 平方米这个值代表一个中等渗透率的砂岩岩心比较有代表性。如果渗透率太低注水压力就需要给得很大早期算起来容易发散如果太高油相又很容易被一下子全部推出缺乏“指进”效果不方便观察现象。流体的物性参数也很关键。水相密度设为998 kg/m^3动力粘度1e-3 Pa·s油相密度设定稍轻一点800 kg/m^3粘度设为50e-3 Pa·s。这组参数很好的模拟了一种典型的中质原油粘度是水的50倍这样的粘度差最容易在水驱过程中出现粘性指进现象结果直观且具有代表性。我建议你在第一次做的时候不要用太极端的参数先保证容易收敛模型跑通后再把粘度拉大、换渗透率做敏感性分析。2.2 相对渗透率与毛细压力模型的计算两相流Darcy接口里水相和油相各自的流动方程是分开写的但两相的相对渗透率和毛细压力都会统一挂在同一个液相属性节点下。Comsol里比较常用的是Brooks-Corey模型它的数学形式是水相相对渗透率 k_rw (S_e)^( (2 3λ)/λ )油相相对渗透率 k_ro (1 - S_e)^2 * (1 - S_e²)其中S_e是有效水饱和度S_e (S_w - S_wr)/(1 - S_wr - S_or)S_wr是束缚水饱和度S_or是残余油饱和度。我在模型里给的束缚水饱和度S_wr是0.2残余油饱和度S_or是0.15孔径分布指数λ取2。这样算下来当含水饱和度为0.2时水相相对渗透率为0即水完全不流动饱和度为0.85时油相相对渗透率为0即油被锁住无法移动。这些参数看着简单但实际经验是饱和度端点值设置不当会导致某个相在边界面上的计算不收敛尤其是初始时刻饱和度若接近残余油饱和度那油相相对渗透率趋近于零方程数值性质会变得很差。所以我建议把初始含水饱和度设在0.21左右而不是正好等于束缚水饱和度这样既接近实际情况又不至于让油相流动性一开始就为零。毛细压力我用的是van Genuchten形式p_c -p_0 * ((S_e)^(-1/m) - 1)^(1/n)其中p_0取5000 Pam和n由经验公式关联。毛细压力在宏观尺度上的具体值虽然比注水压差小不少但它在饱和度前缘的玫瑰花状分布中起关键作用不能直接设为0否则前缘推进形态会产生偏差。2.3 初始条件与边界条件的设置思路初始条件设置很简单把整个计算域初始水饱和度设为0.21油饱和度是0.79。初始压力场我并没选择零压力开始而是做了一个静水平衡式的估算让注入压力从入口到出口线性分布这样求解器在第一步迭代时就不会出现压力剧烈调整。边界条件我分成左、右、上、下四个边界来处理。左边界是注入端我给定一个速度型条件也就是质量流量边界让水以恒定速率注入换算成Darcy速度大约为1e-4 m/s这个值很小确保渗流处于层流Darcy适用范围内。右边界是生产端设为定压边界压力为0即表压为0的开放出口。上下边界做对称或者无流动边界。有一个细节容易忽略就是在地层/岩心出流边界上。如果直接使用“通量”边界而不限制流体回流一旦局部压力反超出口压力油或水可能会从外面回流进计算域出现非物理震荡。我的做法是加了一个“仅流出”条件这是Comsol提供的一个约束选项相当于单向阀只允许计算域里的流体流出去不允许外部流体倒灌回来效果很稳定。3. 网格划分与求解器设置的实操要点3.1 网格策略粗网格为什么会让前缘糊掉网格这件事可以说是多相流模拟里最让人头疼的一环。我前后试过三种网格密度100 × 20200 × 40和400 × 80。对比下来100 × 20这个极度粗糙的方案算出来饱和度云图虽然也能看出大致的前缘推进但过渡带特别宽含水饱和度从0.2到0.8跨越了接近四分之一的模型长度这种情况物理上不合理属于典型数值弥散。我最终用的是200 × 40的映射网格在注水入口和出口附近做了局部加密。做法是先给上下左右边界定义好分布点击“映射”生成矩形网格然后用一个分布节点把左边界附近40%的网格密度提高一倍。这个做法其实就是在保证精度的同时不让计算量爆炸因为前缘先从左端发育那边的饱和度梯度最大。对于Darcy模型因为用不到边界层的壁面定律不需要像CFD那样去画很薄的边界层网格。实际上在宏观多孔介质Darcy求解中边界处的速度满足的是无穿透条件不存在速度梯度剧烈变化所以强行加密边界层反而只会增加无意义的计算开销。网格质量控制我是这么把握的先跑一个200×40的网格记录入口流量、出口采出量这些关键响应值再加密到400×80如果两次计算结果的采收率曲线偏差小于2%说明网格已经收敛。如果偏差大再回头调整。3.2 求解器与时间步的控制两相流Darcy问题是一个典型的瞬态强非线性问题求解器用全耦合线性解算器我用PARDISO。PARDISO这种直接求解器可以避免Krylov迭代方法在处理不对称系数矩阵时出现的收敛不稳代价是内存占用较高但对于二维中型规模模型2万个自由度以内完全没问题。时间步进方式我建议选BDF最大阶次设到2。BDF在刚性问题中表现比广义alpha更稳定而水驱油前缘移动过程中饱和度方程经常会从平滑变剧烈用BDF能减少振荡。时间步控制是最影响效率的地方。我最初直接用固定时间步长0.01 s跑了一个小时后才推进到一半时间效率实在太低。后来我改成自由时间步并设置一个初始步长为1e-4 s、最大步长5 s配合相对误差10^-4做自适应控制。这样前缘位置变化快时会自动加密时间步前缘移动到中段后步长会逐渐放大整体计算时间能缩短三分之二以上。在后期油水前缘快到出口时建议手动把最大时间步降低到1 s不然采样点曲线会出现台阶状跳动。这个台阶状问题如果发生不要怪物理模型大概率是时间步太大导致饱和度前缘“跳跃”过了监测点位置。3.3 如何判断计算结果已经收敛判断收敛不能只顾着看求解器有没有报错还要观察几个关键量的曲线是否平滑。我主要看两个出口含水率随时间变化曲线以及整个模型区域内的平均油饱和度变化曲线。如果出口含水率曲线出现毛刺、抖动首先排查是不是时间步太大如果含水率曲线根本没有上升趋势大概率是初始条件设置时油和水饱和度给错了导致油相一开始就不能流动。还有一次我碰到了高振荡警告检查后发现是某些网格单元里的饱和度降到负值此时需要在物理接口设置中开启饱和度下限限制再把初始水饱和度从0.2调到0.21问题就消失了。4. 结果后处理与数据提取4.1 观察油水前缘推进形态计算完成后最有成就感的一步就是点开“二维绘图组”选“表面”把表达式切换成水的饱和度sw然后播放时间动画。可以看到水的饱和度从左边界开始一点一点向右侧推进前缘形态不是一条直线而是呈现出细长的指状凸起这就是粘性指进现象。细想为什么会形成指进模型里油相粘度是水相粘度的50倍高粘度油对水产生了很强的流动阻力水更倾向于从已经形成水流的低阻力通道往前突进这种不稳定推进在图上看就是一条条往前伸的“手指”。这里有一个重要的后处理技巧饱和度云图的默认色标范围是全局最小到全局最大但在初始时刻全局最小可能显示成0.19左右这不方便观察。我习惯把色标范围锁定在0.2到1之间这样初期微小变化也能看得清楚。做法是在表面绘图节点的“范围”选项卡中手动设置最大值和最小值才能让动画的颜色映射更稳定。如果想把界面位置更清晰地提取出来可以画一条等值线在等值线表达式里填sw设为0.5这条线就是油水界面的大致位置。观察它在不同时间点的位置差异可以直接看出推进速度是否均匀。4.2 含水率与采出程度曲线怎么画光看云图还不够工程上总要用曲线做定量分析。我通常是定义几个全局计算表达式出口水通量除以出口总通量就是产水率fw平均初始油饱和度减当前平均油饱和度再除以初始油饱和度就是采出程度。绘制这两个量的时间曲线时Comsol的“派生值”菜单下的“全局计算”可以直接处理。先选中全部域在表达式里写mean(sw)然后把它和时间关联起来生成一个采出程度时间曲线。这个曲线倒U型关系很明显刚开始一段时间水没到出口产油量以较高水平持续输出当前缘突破之后产水率快速上升采油速率迅速下降采出程度曲线也随之变得平缓。曲线里最能说明问题的指标是“见水时间”也就是含水率开始上升的时间节点。它和模型里的注入速度、两相粘度比都有直接关系。我后来做了三组不同注入速度的对比模拟明显看到注入速度越快见水时间越早但总采收率反而稍微下降原因是高速驱替更容易把水从高渗通道短路掉低渗区里的油没有足够时间被置换出来。这个结论在油田开发方案中很有参考价值。4.3 残余油是怎么分布的最后再分析一下残余油的分布。查看油相饱和度so等于1-sw的云图可以看到在经过长时间水驱后模型的左下角和右下角区域仍然保留着较高油饱和度。左下角靠近入口的是由于水的推进主要从中央通道通过角落处流速极低油大部分被绕流截留右下角则是靠近出口末端水驱的波及系数有限油无法被动用。为了量化残余油我在整个域内用积分算子算了一遍油饱和度对面积的积分得到残余油体积分数我把这个数据做成柱状图分配到不同分区就能清楚看到主流通道和被绕过区域的残余油占比差异。这些信息对指导后续的“调剖堵水”措施很有帮助知道油主要被截留在哪个区域才能决定是注聚合物调整粘度比还是细化注采井网。5. 常见问题与排查技巧实录5.1 典型报错与解决方法我在做这个模拟时前后遇到了几个典型问题挑有代表性的列出来给你排查时做个参考。求解器报“找不到初始解”是最常见的。我第一次遇到时第一反应是模型方程写错了实际上检查后发现是材料属性的单位没对上。分子级渗透率有mD和m^2两套单位我一开始顺手填了100 m^2这直接让Darcy方程快溢出。把渗透率改成1e-13 m^2后模型立刻正常。第二种情况是收敛迭代次数过多每步都要迭代几十次才能满足容差。这时候先看相对渗透率曲线端值。如果初始饱和度正好设在束缚水饱和度0.2那么水的相对渗透率为0方程雅可比矩阵出现退化特征导致每步迭代都很艰难。解决方法是把初始水饱和度改为0.21让方程水相部分保有正梯度收敛速度明显加快。还有一种情况是在饱和度前缘位置的网格单元上出现振荡表现为云图上有一条条带状噪声。这通常是因为局部时间步太大或者网格太粗。我会优先回到网格设置把前缘可能经过的区域细分然后把BDF阶次降为1时间步最大限制收紧基本可以消除。5.2 排查思路的优先级如果你遇到水驱油模拟不收敛我建议按这个顺序排查先看材料属性和单位再看初始条件然后看边界条件最后才调网格和求解器。很多新手一上来就改网格其实往往问题出在最基本的单位换算上。单位问题用Comsol的“单位检查”工具可以直接排查初始条件则看饱和度是否落在合理区间尤其是油相饱和度必须留出流动空间边界条件重点检查出口是否可能倒流用之前提到的“仅流出”条件可以一劳永逸解决。网格方面我建议用“网格统计”查看最小单元质量如果最小质量低于0.3哪怕是拉格朗日型Darcy方程也会开始出问题。在二维矩形几何下映射网格的单元质量一般能保持在0.9以上如果哪些区域单元质量过低直接局部重画即可。5.3 几个额外的实用细节最后分享三个容易被忽略但很务实的小细节。第一多物理场接口中如果开启了“自动开启一致稳定化”在Darcy类的多孔介质流动中一般不要关闭。这个选项会添加流线扩散和侧风扩散对抑制前缘震荡非常有帮助。有些CFD背景的用户习惯把稳定化关掉但在多孔介质Darcy框架下这会导致解出现严重振荡。第二采样点时不要直接选网格节点。我用点监测功能在整个区域布置了若干监测点后发现位置稍微偏离网格节点时曲线会抖动得很厉害。后来我把监测点全部放在几何顶点或边线上的均匀位置数据就平滑了其实底层原因是插值失真。第三模型跑完后记得多保存几个状态。我通常会在时间推进过半时存一个中间结果文件即处理不同敏感性参数时可以直接以这个状态为初始值继续跑不用从头算起。这个习惯在后来的多工况对比中为我节省了大量时间。做实操类仿真的乐趣就在于从“方程写出来”到“图像动起来”的这段距离里每一个参数的选择背后都有它的理由。水驱油这个经典问题恰好把这些关键环节密集地串在了一起。有了这套思路和排查经验相信你也能顺利完成自己的模拟实验。