基于Level Set有限元方法的微观水驱油数值模拟

高亚军，姜汉桥，王硕亮，刘传斌，常元昊,王依诚

(1. 中国石油大学(北京)石油工程教育部重点实验室，北京102249；2. 中国地质大学(北京)能源学院)



基于Level Set有限元方法的微观水驱油数值模拟

高亚军1，姜汉桥1，王硕亮2，刘传斌1，常元昊1,王依诚1

(1. 中国石油大学(北京)石油工程教育部重点实验室，北京102249；2. 中国地质大学(北京)能源学院)

在微观水驱油实验的基础上，建立了二维微观孔隙模型；引入Level Set数学方法，结合N-S方程,建立了微观两相渗流数学模型，借助有限元方法，进行水驱油两相数值模拟，研究微观水驱油动态特征。研究了孔喉非均质性、润湿非均质性、油水两相黏度以及驱替压差对驱油效果的影响；通过对模拟结果中孔隙内流体的波及效率、两相界面的运移速度以及优势渗流通道的分析发现，Level Set方法能很好地处理各因素下油水两相界面的拓扑变化，为微观水驱油机理的研究提供了一种新的技术方法。

水驱实验；数值模拟；微观驱替特征

许多学者在微观水驱油机理和两相流动特征的研究方面做了大量的实验，实验方法主要有岩心薄片平面驱替和玻璃刻蚀物理模拟等二维方法[1-3]，以及驱替岩心CT扫描技术的三维方法[4-5]，但以上方法实验模拟仿真成本高，操作较繁琐和困难,并且实验监测技术和图像获取精度有限，导致实验数据误差较大。数值模拟方法不需要实际岩心参与，利用计算机对微观岩心孔隙模型进行刻画，模拟实际水驱油整个过程的演化规律，研究油水两相动态变化特征，具有便捷和可视化等优点。目前应用较多的是孔隙网络模型[6-7]，但此方法一方面无法更进一步清晰呈现孔道内流体两相界面运移状态，另一方面无法直观获取各流体质点在任何时刻的压力和流速变化。

针对以上方法所存在的问题，本文提出基于Level Set(水平集)方法的微观水驱油两相数值模拟，对孔隙内流动特征进行可视化研究。水平集方法是由Osher和Sethian[8]提出的，最初主要应用于智能控制、图像处理等方面，近年来经国外学者对求解方法的不断改进[9-10]，逐渐被用于气液或液液两相流动的数值模拟研究领域中。有学者用此方法做过简单模型的油水两相模拟[11-12]，但局限于理想孔隙结构或油气两相驱替，均未与实验进行对比验证，且未涉及到微观孔隙内的油水两相模拟。为此笔者根据岩心CT扫描图像，进一步处理后，建立二维微观孔隙模型，结合N-S方程，引入Level Set方法来对油水两相界面的动态运移进行追踪，并用有限元方法求解。

1　油水两相驱数学模型

1.1Level Set方法

假设地层中岩石孔隙内水油两相流体均为不可压缩流体，流动状态为层流。在层流两相流中，水平集界面自动生成界面运移的方程，在Level Set方程中，当φ=0时的水平集轮廓线代表两相流体界面。对于此油水两相，在油中φ=0，水中φ=1。Level Set方程可以被看作是表示水油两相流中水的体积分数，因此，两相流体的界面运移方程可以用下式表示：

(1)

式中：φ为油水两相界面轮廓线；γ为方程求解中的重新初始化参数；t为两相作用时间；u为流速；ε为界面厚度，设置值一般小于计算模型中网格剖分的最小单元。

借助于Level Set函数 ，流体的物性可用下述方程表示，即：

(2)

式中：ρo、ρw分别为油和水的密度；μo、μw分别为油、水的黏度。经过这样定义后，流体的密度和黏度就可以在两相接触面处平滑过渡，不会形成突变，保证计算的稳定收敛。

1.2控制方程

(3)

式中：ρ为密度，kg/cm3；μ为动态黏度，Pa·s；u为速度，m/s；I为单位矩阵；p为入口压力，Pa；Fst为油水界面表面张力，N/m 。

表面张力项可以通过下式计算[13]：

Fst=·T

(4)

T=σ(I-(nnT))δ

(5)

式中：n为界面法向量；σ为表面张力系数；δ表示狄克拉函数，在两相界面以外区域此函数值均为零。函数被近似为：

(6)

界面法向量可表示为：

(7)

当用有限元方法求解方程时，先使方程变形，然后在计算区域内分部积分。层流两相流中，为了表示壁面与流体的相互作用力，在计算域上求积分可以用加上一个边界积分的形式表示如下：

(8)

式中：θ为壁面接触角，当使用无滑移边界条件即非润湿条件时，则无此边界项，因为在边界上时test(u)=0，不能设置接触角。然而，如果允许流体细微的边界滑移(可以忽略不计)，就可以设置接触角。上式方程中给出了润湿壁面所增加的边界条件项，因此可以设置两相的壁面接触角。

1.3边界条件

入口、出口边界上的压力由模型设置，两侧边界为壁面封闭边界，初始界面设置为t=0时水相和油相的接触面，壁面条件为可润湿壁面，垂直壁面方向流速为零。即：

(1)入口边界：p=p0,n·μ(u+(u)T)=0

(2)出口边界：p=p1,n·μ(u+(u)T)=0

式中：β是滑移长度，其中，β值的大小为边界壁面部位剖分网格的尺寸大小；Frf为壁面作用力。

2　水驱油仿真模拟

油水两相驱仿真模拟研究表明，构建仿真模型是研究驱替机理的重要环节[13-14]。本文用CT扫描岩心切片处理后建立的玻璃刻蚀模型进行仿真模拟。模型的构建需要经过多次反复试验、验证、修缮，才能最终建立可行的模型。本文水驱油特征动态模拟的基本步骤如下：

(1)对实际储层岩心进行CT扫描，经图像处理后，建立微观孔喉模型。

(2)根据流体运动方程，结合Level Set数学方法，建立数学模型。

(3)对实际物理模型的水驱油实验中压力、黏滞力、温度、润湿性等因素所反映出的两相动态特征进行分析，建立数学模型的边界控制条件。

(4)改变仿真模拟中两相流体参数以及边界控制条件，根据整个仿真图像中油水运移状态，研究不同条件下水驱机理及剩余油的控制因素。

根据上述步骤，首先需要对实际所建立的微观孔喉模型进行水驱实验，为此笔者选取根据实际岩心CT扫描孔喉特征所建立的微观非均质性不同的玻璃刻蚀模型A和B进行驱替实验，对所建立的微观两相驱仿真模型进行对比验证。

2.1水驱油实验

实验中采用光化学刻蚀的玻璃模型A和B，模型尺寸分别为17.45 mm×14.31 mm和15.28 mm×12.13 mm。为了能够较好地体现出油水两相在驱替过程中相界面的动态变化特征，实验所用原油黏度较大，饱和原油黏度分别为50 mPa·s和20 mPa·s，孔隙度分别为为30.96%和28.91%。实验用盐水密度为1.05 g/cm3，黏度为1 mPa·s。实验步骤如下：

(1)安装固定玻璃刻蚀模型，将模型抽空饱和油，然后再进行水驱油；

(2)按实验方案中压力进行水驱油实验，用显微放大系统和录像机录取实验过程的图像，并用计算机图象处理系统对录像和图片进行处理；

(3)对整个驱替实验的录像进行分析，统计不同时刻的驱油效率，分析微观水驱相界面特征变化并解释剩余油的形成。

2.2数值模拟与实验的对比分析

针对玻璃刻蚀模型，将CT扫描的孔隙结构经过二值化转换处理可得到二维孔隙轮廓模型，对孔隙区域进行自由三角网格剖分，建立相关有限元模型进行求解。实验和仿真模拟均从右端注入，左端采出，各参数如表1所示。

表1　数值模拟模型参数设置

数值模拟条件较为理想，而实验条件下影响因素较多，相同含水饱和度下，实际驱替时间与数值模拟的驱替时间不相等。因此，将各阶段记录时间与数值模拟的时间在(0,1)作归一化处理进行对比分析。记录玻璃刻蚀实验条件下各时刻的的含水饱和度，与数值模拟得到的含水饱和度在相同归一化时间下进行对比，发现吻合度很高(图1)。刚注入时刻，由于数值模拟不存在实验条件下注入速度的稳定过程所带来的误差，因此初始时刻的拟合存在一定的相对误差(27.2%)，随着驱替的进行，相对误差越来越小(5.5%)。

图1　实验与数值模拟的含水饱和度对比

将驱替实验A和B所录取的图像与数值模拟结果进行等含水饱和度对比发现(图2)，驱替前缘两相界面吻合度很高，其中，红色为油相，蓝色为水相。A左端见水时的含水饱和度为31.7%，与图1中模型A曲线拐点处的含水饱和度非常接近。相对于模型B，模型A的非均质性较严重，左端见水后，水相突破后流速通道阻力减小，注入水很难波及到其它孔隙，驱油效率降低，最终采收率为41.2%。模型B左端初始见水的含水饱和度为66.8%(图2)，同样与图1中曲线拐点基本一致，见水后驱油效率基本不再增加，整个驱替过程中，油水两相界面推进较均匀，波及范围大，最终采收率达到72.8%。由于模型A的原油黏度远大于B，因此驱替过程中，水相在进入细小孔喉时的阻力非常大，模型A的指进现象较为严重。

在与实验的驱替过程对比之下，发现本文所建

图2　实验与数值模拟的结果对比

立的数学有限元仿真模拟方法能够较好地刻画水驱油整体过程的相界面动态分布，因此可以对其它特征进行进一步研究。

3　微观水驱油特征

3.1微观孔喉非均质性

由于孔隙介质存在微观非均质性，不同半径孔喉的毛管力差异较大，驱替相所受阻力也不同，此时绕流和指进为主要的微观驱油机理。对整个二维平面模型来看，水相由入口端进入，沿着一条或数条阻力较小、孔径相对较大的孔道向前突进并首先到达出口，即微观指进现象。从实验和仿真模拟中均可以看到当，当模型出口端一旦见水，注入水大部分便会只沿着突破通道流动，模型A从左端初期见水时采收率31.7%到最终稳定采收率41.2%，驱油效率增加幅度很小(图1)。模型B的注入水推进较为均匀，指进现象较弱，水驱效果较好(图2b)。对比模型A和B，可以发现A的孔喉非均质性比B严重，模型A中部存在大孔道，模型B孔喉分布相对较为均质，且B的孔喉比A小，虽然驱替需要更大的压力，但是压力在平行驱替方向孔喉内分布较均匀，油水前缘均匀推进，波及效率较高，驱油效果较好。

3.2润湿非均质性

孔隙内油水两相的流动状态与壁面润湿性紧密相关。为了体现润湿性对水驱油时流动孔道的选择特点，在本文所建立的仿真模型中，以壁面接触角表示润湿角。设置部分孔隙水湿接触角为60°(图3，b1)，其余条件不变，在相同驱替时间下，将水湿条件下的模拟结果与孔隙壁面为中性润湿条件下的模拟结果进行对比(图3)，以同一个压力系统相同驱替条件的a2和b2区域为参考基准，发现在水湿条件下，孔隙中油水两相界面运移更快(图3，b1)。部分壁面水湿加快了水相的指进，导致大孔隙绕流现象更严重；大孔隙油水界面的迅速推移，导致小孔隙内憋压不足，油水界面停止向前推移(图3，b2)。经过以上分析，可以发现本文所建立的微观仿真数学模型能够较好地处理壁面润湿性所引起的油水两相运移动态变化特征。

3.3油水黏度比

原油的黏度在流动过程中反映为原油内部的摩擦阻力，决定两相与壁面的摩擦力大小，从而影响两相驱替的流动速度和驱油特征。取同一驱替时刻不同黏度下模拟结果(图4)可见，在相同的驱替时间下，当油水黏度比较小时，驱替阻力更小，两相界面运移更快。油水黏度比越大，水相进入小孔隙的阻力越大，水沿着大孔道的指进现象越来越明显。由图4可看出，level set方法能够较好地处理油水两相黏度变化所引起的相界面的拓扑变化。

图3　模型A不同润湿角下模拟结果

3.4驱替压力

高压驱替时，水沿大孔道中心迅速突破，指进明显，截断细小孔道，使油相从连续相变为非连续相，形成残余油，但高压驱替可将部分可动剩余油驱替出来(图5和图6)。选取其中典型局部区域放大进行详细分析如下：

(1)在油水前缘到达同等距离时，对比区域a和区域b中驱替现象(图5和图6)，发现高压力驱替下的区域b中，油的体积分数达到50%左右，油水过渡带很长，指进现象明显，驱替效率明显低于压力较小时的驱替结果。

图4　模型A不同油相黏度下模拟结果

图5　P=50 000 Pa时模型模拟结果

图6　P=100 000 Pa时模型B模拟结果

(2)对比图中c、d区域，发现驱替压差过大(d区域)，即水相推进速度过快时，水主要沿着大孔道中心流动，绕过连通性不好、孔喉非均质严重的区域，将孔隙角落的原油封堵死，形成角状残余油。

(3)对于一些小孔道(e区域)，毛管力较大，当驱替压力小于原油的流动阻力时，不能驱替小孔喉内的原油，形成可动残余油，残余油主要以塞状或柱状形式存在。

3.5流速可视化分析

在数值模拟流速示意图中，可以清楚地看出存在主要流速通道(图7)。对整个驱替过程中的压力和流速分析表明，水通过优势通道到达出口端时，整个流域的流动阻力降低，其它小孔隙压力相对降低，波及效率降低，有利于剩余油的形成。流体在孔隙通道中心的流速最大，越靠近孔隙壁面，流速越小。那些微小孔隙、封闭空隙中，油水几乎不参与流动，流速较小，孔隙中压力降落较快，坡度较陡；而对于那些大孔隙高流速的空隙中，压力在流动方向上降落较慢。

图7　模型B流速示意图

沿y轴方向对注入端的7个孔隙编号依次为y1，y2，y3，y4，y5，y6，y7，统计同一驱替时间步长下入口端的驱替流速，做出入口端整个纵剖面在固定压力下的注入速度，可以清楚看到在7个孔隙入口处流速大小(Vy2>Vy5>Vy4>Vy7>Vy1>Vy6>Vy3)，入口速度与图7中的主要流速通道的分布相对应。但由于模型的绝对理想化，在与实际驱替规律不变的情况下，模拟驱替流速大于实际的驱替流速。在孔隙y5入口中可以发现，孔喉中心流速最大，越靠近孔隙壁面越小(图8)。其中在孔道中心处可以看到流速并非最大，呈下凹现象，主要是因为y5孔喉入口内部中心岩石颗粒阻挡形成分流所致。在驱替压力梯度不变的情况下，随着驱替时间的进行，孔隙中流速逐渐增大。主要是因为在微尺度驱替下，惯性力可以忽略不计，黏滞阻力占主导地位，当低黏度液体驱替高黏度液体时，黏滞阻力逐渐减小[15]，因此驱替速率逐渐增大。

图8　单孔隙入口流速随时间变化规律

4　结论

(1)引入Level Set数学方法，进行微观二维平面水驱油仿真模拟，数值模拟的结果与实验得到的结果吻合度很高，为微观水驱油的特征研究提供了新的可视化技术手段。

(2)Level Set方法在处理油水两相黏度、润湿性以及驱替压力等引起的相界面拓扑变化时，具有较高的仿真度，可用此方法来进行更为复杂的微观流体流动机理研究。

(3)Level Set数值模拟方法具备常规物理仿真和数值仿真所不具备的优点，能够清楚研究孔隙内任何部位流体的压力、流速等变化特征。油水两相流在孔道中心流速最大，越靠近孔隙壁面流速越小；孔径越小，压力降落越快；在驱替压力梯度和油水黏度比恒定时(油相黏度大于水相)，孔道内油水两相驱替速度逐渐加快。

[1]高辉，孙卫，路勇，等.低渗透砂岩储层油水微观渗流通道与驱替特征实验研究——以鄂尔多斯盆地延长组为例[J].油气地质与采收率，2011，18(1)：58-62.

[2]李洪玺，刘全稳，何家雄，等.物理模拟研究微观剩余油分布[J].新疆石油地质，2006，27(3)：351-353.

[3]刘太勋，徐怀民.扇三角洲储层微观剩余油分布模拟试验[J].中国石油大学学报：自然科学版，2011，35(4)：20-26.

[4]吕伟峰，冷振鹏，张祖波，等.应用CT扫描技术研究低渗透岩心水驱油机理[J].油气地质与采收率，2013，20(2)：87-90.

[5]曹永娜.CT扫描技术在微观驱替实验及剩余油分析中的应用[J].CT理论与应用研究，2015，24(1)：47-56.

[6]高慧梅，姜汉桥，陈民锋，等.储集层微观参数对油水相对渗透率影响的微观模拟研究[J].石油勘探与开发，2006，33(6)：734-737.

[7]陈民锋，姜汉桥.基于孔隙网络模型的微观水驱油驱替特征变化规律研究[J].石油天然气学报，2006，28(5)：91-95.

[8]Osher S, Sethian J A.Fronts propagating with curvature-dependent speed:Algorithms based on Hamilton-Jacobi formulations[J].Journal of Computational Physics,1988,79(1):12-49.

[9]Sussman M.A sharp interface method for incompressible two-phase flows[J].Journal of Computational Physics, 2007,221(2):469-505.

[10]Olsson E, Kreiss G.A conservative level set method for two phase flow[J].Journal of Computational Physics, 2005,210(1):225-246.

[11]王琳琳，田辉，李国君.基于Level Set方法对油水和气水两相界面的数值模拟[J].应用力学学报，2010，27(2)：298-302.

[12]Qianlin Zhu,Qianlong Zhou,Xiaochun Li.Numerical simulation of displacement characteristics of CO2injected in pore-scale porous media，2016,81:87-92.

[13]孙焕泉，孙国，程会明，等.胜坨油田特高含水期剩余油分布仿真模型[J].石油勘探与开发，2002，29(3)：66-68.

[14]徐守余，朱连章，王德军.微观剩余油动态演化仿真模型研究[J].石油学报，2005，26(2)：69-72.

[15]A.R Ershov,Z M Zorin,V D Sobolev, et al.Displacement of silicone oil from the hydrophobic surface by aqueous trisiloxane solutions[J].Colloid Journal，2001,63:290-295.

编辑：李金华

1673-8217(2016)05-0091-06

2016-05-01

高亚军，1991年生，2015年毕业于中国地质大学(北京)石油工程专业，在读硕士研究生，从事油气田开发方向研究。基金项目：国家重点技术研究发展计划(973计划)项目“致密油高效开发油藏工程理论与方法研究”(2015CB250905)。

TE311

A