气液非混相驱替过程中的卡断机理及模拟研究1)

张晟庭 * 李 靖 *, 陈掌星 *, 张 涛 ** 吴克柳 * 冯 东 * 毕剑飞 * 李相方 *

* (中国石油大学 (北京) 油气资源与探测国家重点实验室,北京 102249)

† (卡尔加里大学化学与石油工程系,加拿大阿尔伯塔 T2N1N4)

** (西南石油大学油气藏地质及开发工程国家重点实验室,成都 610500)

引言

复杂多孔介质中的非混相驱替在两相流动中占据重要地位,并在气驱,气水交替及泡沫驱等提高油气采收率领域具有重要研究意义[1-5].在气驱提高油气采收率过程中,连续气相注入储层后,将会形成复杂的气液非混相驱替模式;注入储层中的气相一般为非润湿相,在非润湿相驱替润湿相的过程中,孔喉壁面上将会滞留一层液膜,液膜的存在可能会造成驱替过程中发生卡断甚至水锁现象.因此,在孔隙尺度研究气液非混相驱替过程中的卡断机理有利于人们认识微观孔喉中的气液流动状态.

研究人员通过理论研究,实验及数值模拟等方法,对气液两相驱替过程中的卡断机理展开了大量研究工作.在理论方面,Roof[6]基于孔喉毛管压力不平衡这一机制,建立了预测卡断模型的静态准则,认为当最小喉道半径小于孔隙半径的一半时,润湿相将沿着孔隙壁面发生回流,也即导致非润湿相发生卡断.Gauglitz 等[7]给出了收缩圆柱形孔喉系统的卡断准则,认为毛管数Ca对卡断发生的时间有一定影响,并且对于不同宽度比的孔喉系统影响规律不同.Ransohoff 等[8]给出了方形喉道截面的以及三角形喉道截面的静态卡断准则,并认为在非圆形截面的孔喉中更容易发生卡断现象.然而,上述卡断准则并没有考虑动态参数对卡断影响,Tsai 和Miksis[9]通过数值方法发现毛管数过大或者过小时,气泡不会发生卡断现象.Deng 等[10-11]扩展了圆形及非圆形缩颈孔喉系统的静态卡断准则,发现之前的静态准则对卡断的预测过于保守,同时发现即使喉道系统达到了卡断发生的静态准则,高毛管数也会抑制卡断的发生.在实验方面,随着微流控模型的发展,大量学者采用微流控模型对两相流动的卡断现象进行可视化实验,Tian 等[12]利用2D 微模型研究了气水流动过程中的卡断现象,并验证了Roof 等提出的静态准则.Cha 等[13]利用2D 矩形截面微模型研究了气泡卡断机理,并与提出的卡断理论模型得到了较好的一致性.Tetteh 等[14]利用微流控模型观察到了低矿化度水驱过程中的油相卡断.Wu 等[15]利用实验研究非缩颈孔喉的气泡卡断现象,并且认为喉道长度及喉道宽度对卡断有一定影响.

相比于实验,数值模拟方法可以更好的帮助理解及预测孔隙尺度气液两相流动状态及卡断现象.在孔隙尺度上模拟多相流的方法有孔隙网络模型(PNM)[16]、光滑粒子动力学方法(SPH)[17]、密度泛函流体动力学方法(DFH)[18]、流体体积方法(VOF)[19]以及格子玻尔兹曼方法[20]等.部分学者采用VOF方法研究了缩颈孔喉系统中的卡断现象,Raeini 等[21]利用VOF 方法研究了星形喉道截面的卡断现象,并考虑了动态因素的影响.Starnoni 和Pokrajac[22]采用VOF 方法考虑接触角及黏度比对卡断影响,研究表明,喉道截面圆度越高发生卡断的临界接触角越小,黏度比的增大也会降低发生卡断的临界接触角.Zhang 等[23]采用VOF 方法研究了方形截面喉道的卡断现象,并认为高毛管数及高黏度比可以有效抑制卡断的发生,Cha 等[13]利用VOF 方法模拟了矩形截面的卡断现象,并给出了卡断发生的孔喉宽度比条件.然而,VOF 方法很难准确计算界面的法向和曲率,且界面重构方法复杂,向高维推广困难[24].格子玻尔兹曼方法是介于微观分子动力学模拟方法和宏观传统数值模拟方法之间的一种介观数值模拟方法,与VOF 方法等界面跟踪或界面捕捉方法相比,格子玻尔兹曼方法中的气液界面可以自动产生,演化与迁移,无须采用任何界面跟踪或界面捕捉技术[25],并能直接处理流体与流体间的作用力及流体与固壁间的作用力[26].目前多相格子玻尔兹曼模型可以分为四类,包括颜色梯度模型[27]、伪势模型[28]、自由能模型[29]以及相场模型[30].张磊等[31]基于颜色梯度模型在孔喉系统中模拟非混相驱替过程,并认为多孔介质的非均质性越强,流体越容易发生卡断.赵玉龙等[32]基于颜色梯度模型,模拟地层高温高压条件下致密气驱替地层水的流动过程,并发现在多孔介质中大量连通的微小通道被卡断水占据.Alpak等[33]采用相场模型模拟3D 缩颈孔隙中的两相驱替过程,观察到了卡断现象并与Roof[6]静态准则保持较好的一致性,但是模拟两相流体的黏度相差不大,不能反应气液两相的黏度比.Wei 等[34]采用伪势模型模拟单个气泡通过喉道过程,并观察到了卡断现象.Zhang等[35]采用改进的双组分伪势模型模拟气驱水过程,成功在喉道壁面捕捉到了水膜的存在,并认为微纳尺度气水两相流动过程中壁面液膜的厚度不可忽略.

但是,目前对卡断现象的研究往往建立在缩颈渐变孔喉系统中,弱化了壁面液膜厚度的影响,同时宏观VOF 方法不能充分表征壁面和液膜相互作用.为此,本文在原始的伪势格子玻尔兹曼方法的基础上,改进了流体之间作用力格式并添加流固作用力以捕捉气液两相驱替过程中喉道壁面上的液膜变化规律及其对气液两相流动状态的影响.并进一步研究了驱替压差(毛管数Ca),喉道长度及宽度等因素对液膜厚度及流动状态的影响规律.本文的研究内容为表征复杂多孔介质的卡断机理提供了理论及模拟基础.

1 单组分多相伪势模型

1.1 经典BGK 碰撞模型

在格子玻尔兹曼方法中,采用离散的密度分布函数表征流体的运动规律,标准BGK 碰撞过程可以表示为[36]

式中,fi(x,t)为i方向密度分布函数,Δx和Δt分别表示网格步长和时间步长,τ表示松弛时间,其取值与流体的运动黏度ν相关,表示为

ei为各方向的离散速度,对于D2Q9 模型表示为[37]

1.2 流体之间作用力及流固作用力

原始伪势模型中流体之间的作用力定义为[28]

式中,G为格林函数,ψ(x,t)为有效质量,其定义为ψ(ρ)=ρ0[1-exp(-ρ/ρ0)],ρ0=1 为参考密度.进一步,考虑第一层粒子之间的作用力,上述作用力格式可以简化为

式中,g为流体之间作用强度的控制参数,w(|ei|2)为权重系数,|ei|2=1,w(|ei|2)=1/3;|ei|2=2,w(|ei|2)=1/12.采用原始的伪势模型的作用力格式会在两相界面处存在较大的虚假速度.因此,本文采用Gong 等[38]提出的作用力改进格式以降低虚假速度,具体表示为

进一步考虑最近及次近的节点之间作用力,Chen 等[39]提出了上述作用力的简化格式

式中,β为调节参数,通过改变β大小可以改善模型的热力学一致性并降低虚假速度.若仅考虑第一层节点的作用力,可将式(9)进一步扩展为[40]

式中,γ1和γ2为作用力格式的离散系数,其中γ1=1/3,γ2=1/12.上述离散作用力格式也称为E4 作用力格式.读者可以参考文献[40]的研究,将作用力扩展至E6 或E8 作用力格式获得更低的虚假速度及更精确的界面张力.

流体与壁面的作用力类似于流体之间的作用力格式,其定义如下[41]

式中,s(x+eiΔt,t)为判断流体节点和壁面的标识函数,如果x+eiΔt是固体节点其值为1,x+eiΔt是流体节点则其值为0.ρw为壁面密度,通过调节ρw的大小可以模拟不同的接触角.

1.3 状态方程

原始的伪势模型的状态方程可以通过泰勒展开获得[28],具体表示为

当采用上述的状态方程时,会导致严重的热力学不一致性,同时模拟的最大密度比和黏度比均在10 以内.文献[42]提出将真实的流体状态方程与伪势模型结合,大大提高了伪势模型的热力学一致性.本文采用常用的Redlich-Kwong (RK)状态方程表征气液两相,RK 状态方程的定义为

式中,pEOS为实际流体的压力,R为气体常数,T为系统的温度,ρ为实际流体的密度,a,b为临界参数,其中a=0.427 48R2Tc2.5/pc,b=0.086 62RTc/pc.Huang等[43]给出了伪势模型中的参数取值,a=2/49,b=2/21,R=1,临界温度Tc=0.196 1,临界压力pc=0.178 4,临界密度ρc=2.988 7.在添加RK 状态方程后可以采用下述方法计算伪势

值得注意的是,当采用上述格式计算伪势时,g将会在计算压力p的时候消掉,不再具有实际物理意义,仅作为保证根号内的符号为正的作用,故本文的模拟中取g=-1.

1.4 外力项处理格式

在添加了流体-流体作用力及流固作用力后,本文采用Kupershtokh 等[44]提出的精确差分方法(EDM)格式处理外力项,其定义如下

式中,Ftotal表示流体之间作用力与流固作用力之和,EDM 格式能够适应较宽的温度范围,是目前最常见的外力处理格式之一.

1.5 边界条件

格子玻尔兹曼方法中的边界条件对数值计算精度,计算稳定性以及计算效率具有很大影响[45].常用于模拟两相流动边界条件包括周期边界以及无滑移反弹边界等.其中,周期边界指流体粒子从一侧边界离开流场时,在下一时间步会从流场的另一侧流入流场,周期边界可以保证整个系统的质量和动量守恒.传统的周期边界采用增加虚拟流体节点(x0,xN+1)的方法实现,可以表示为[46]

对于静止的固体边界,常用的处理方法就是对边界上的粒子进行反弹处理,反弹边界可以严格保证边界上速度无滑移,反弹边界可以表示为[46]

式中,ybottom和ytop分别为位于底部及顶部的固体边界.

2 模型验证

2.1 热力学一致性验证

热力学不一致性一直是伪势模型中的主要问题之一,将直接影响气液密度比及界面张力的计算精度.因此,需要首先评价模型的热力学一致性.考虑到气液弯液面引起的毛管压力会对气液密度有一定影响,采用假想的平直界面气液模型来研究气液共存密度,并与麦克斯韦理论计算结果进行对比来评价模型的热力学一致性.

首先,建立31 × 201 的格子空间,在x和y方向均设置为周期边界,只改变温度,以模拟不同温度下的气液共存密度,采用Huang 等[47]密度初始化方法提高数值稳定性,该方法可以表示为

式中,W为相界面宽度,一般为2～ 5 个格子,ρv和ρl分别为麦克斯韦理论的气相和液相密度.采用上述密度初始化方法之后,在50 ≤y≤150 的位置为液相,其他位置为气相,气液界面为直线,可以忽略毛管压力作用[47].对比β=1 和β=1.125 及β=1.25 三种不同参数条件下的气液共存密度和麦克斯韦理论计算的气液共存密度,结果如图1 所示,β=1.125 时,LB 模型与麦克斯韦理论所得的气液共存密度基本一致,据此采用β=1.125.

图1 气液共存密度曲线对比Fig.1 Comparison of gas-liquid coexistence density curves

2.2 界面张力验证

在验证热力学一致性之后,通过模拟静态圆形液滴来计算界面张力.首先,构建3 种不同网格尺寸(101 × 101,151 × 151,201 × 201)的格子空间,模拟静态圆形液滴密度分布并进行网格独立性检验.采用Huang 等[48]提出的密度初始化方法

式中,R0为液滴初始化半径,(x0,y0)为液滴初始圆心位置,模拟条件选取为:T=0.8Tc,ρv=0.342,ρl=6.60,四周采用周期边界,其余参数与2.1 节中保持一致.为了便于比较,将不同网格尺寸中液滴的直径与网格尺寸的比值保持一致,模拟时间步选为50 000步,足以使模拟结果达到稳定.模拟结果如图2 所示,不同网格尺寸下中心线处(x=Nx/2)对应的气相及液相密度分布几乎一致,说明模拟结果具有网格独立性.模拟的密度分布结果仅在两相过渡区域存在微小偏差,当网格密度超过151 × 151 时,密度分布的模拟结果几乎一致,足以满足模拟精度要求,为了节约计算资源量,选取网格尺寸为151 × 151 的格子空间进行界面张力模拟验证.

图2 中心线处密度分布Fig.2 Density distribution of center line

进一步模拟T=0.7Tc,T=0.75Tc以及T=0.8Tc三种不同温度条件下圆形液滴的静态行为,模拟达到平衡后根据等式(14)计算实际的气液压差,利用Young-Laplace 方程ΔP=Pin-Pout=σ/R确定界面张力[49].本文认为气液的分界面处的密度为(ρv+ρl)/2,并以此为根据确定液滴直径大小,通过改变液滴的不同直径,获得一系列的压差和液滴半径的对应关系.如图3 所示,压差与半径的倒数呈现良好的线性关系(R2均大于0.999),通过线性拟合得到不同温度下气液界面张力大小分别为0.463 lu,0.361 lu,0.264 lu (其中lu 代表格子单位).

图3 界面张力验证Fig.3 Verification of interfacial tension

2.3 静态接触角验证

本节在考虑流体-流体的作用力的基础上,添加流体与固体壁面之间的作用力,模拟液滴与壁面间的不同接触角.Chen 等[39]给出了一种液滴在壁面上的接触角计算方法

式中,θ为液滴的润湿角,ξ1为液滴与固体壁面接触线的长度,ξ2为液滴的高度.为了模拟验证不同的润湿角,建立151 × 151 的格子区域,在x方向设置周期边界,在y=0 及y=Ny设置为固体壁面,并采用反弹边界[50],模拟温度T=0.8Tc,密度初始化方法采用Li 等[51]提出的密度初始化方法,其具体格式如下

式中,R0为液滴初始化半径,本文取15 lu,x0,y0为液滴初始化圆心位置,本文设置(x0,y0)=(75,135).模拟结果如图4 所示,通过改变不同的壁面密度ρw可以模拟0°～ 180°的润湿角.

图4 不同静态接触角模拟Fig.4 Simulation of different static contact angles

2.4 角隅液相滞留特征

在方形截面孔隙及喉道中,气相驱替液相之后,液相将会在孔隙及喉道的角隅处滞留,液相的滞留特征(曲率半径)对非圆形孔隙中的卡断机理有一定的影响[52].为此,本文采用提出的伪势模型模拟液相在角隅的曲率半径与液相饱和度Sl的关系,验证模型的准确性.建立151 × 151 的格子空间,模拟温度T=0.8Tc,四周均为固体壁面并采用反弹边界,润湿角设置为0°,其余参数与2.3 节中一致.在文献[51]基础上提出一种新的初始化密度方法以表征液相在角隅中的不同曲率半径,具体格式为

式中,(x0,y0)=(75,75) 为气泡初始圆心坐标位置,R1为初始气泡半径,将R1的取值范围设为75～之间即可表征角隅液的不同曲率半径(0～ 75).如图5 所示,所提模型可以捕捉到壁面的液膜及液相在角隅处的曲率半径.壁面上静态液膜厚度的表征方法可以参考文献[53-57],本文主要聚焦评价液相在角隅处的不同曲率直径d及饱和度的Sl关系.

图5 角隅液相滞留特征Fig.5 Corner liquid retention characteristics

通过将液相在角隅的曲率半径与液相饱和度模拟结果与Li 等[58]提出的方形孔隙角隅液的曲率半径与液相饱和度的关系相对比来验证模型的准确性,若不考虑壁面液膜厚度的影响,理论模型可以简化为

式中,ε为无量纲曲率直径,ε=d/(L1-2h),L1为方形孔隙的长度,h为液膜厚度,d为角隅液的曲率直径.图6 给出了模拟结果与理论结果的对比,可以看出,模拟结果与理论结果基本一致,说明所提模型可以准确评价液相在角隅的赋存规律以及确定角隅液的曲率半径.

图6 角隅液曲率直径与饱和度的关系Fig.6 Relationship between corner liquid curvature diameter and saturation

2.5 网格独立性检验

通过模拟不同网格尺寸条件下单管内气相驱替液相过程(弯液面顶点位置随时间变化),对模型进行网格独立性检验.建立4 种不同的网格尺寸(241 ×41,301 × 51,361 × 61,421 × 71),将y方向的底部及顶部设置为固体边界并采用无滑移反弹格式.在入口和出口分别设置为5～10 格子的气相缓冲区域,其余位置设置为液相,在每个格子节点添加外力驱动流动,在流动方向(x方向)采用周期边界,当液相被驱替至出口的气相缓冲区域中时,液相密度在下一个时间步被自动替换为气相密度,实现缓冲区域内液相“消失”,从而实现驱替过程[35].模拟条件选取为:T=0.8Tc,ρv=0.342,ρl=6.60,将模拟的驱替压差Δp以及模拟时间步长与取点时间间隔的比值保持一致.模拟结果如图7 所示,对于不同的网格尺寸,模拟结果几乎一致.因此,在后续驱替的模拟中,选取的模拟条件网格尺寸为301 × 51.

图7 网格独立性验证Fig.7 Grid independence verification

3 动态卡断机理及影响因素分析

3.1 静态卡断准则评价

Ransohoff 等[8]基于Roof 提出的卡断判断机制,(即喉道毛管力Pc-neck大于前缘毛管力Pc-front),提出了方形孔喉截面的静态卡断准则.其核心假设是在方形截面孔喉系统中,液相(润湿相)在角隅处以角流为主要流动方式,如图8 所示,忽略壁面液膜的影响,缩颈方形截面孔喉系统的卡断判断准则可以表示为

图8 气液两相流动示意图Fig.8 Schematic diagram of gas-liquid two-phase flow

式中,Rc为孔隙曲率半径,Rt为喉道曲率半径,Rzt为喉道横向曲率半径,Cm为无量纲曲率,与孔隙截面形状有关.基于最小表面能模型,Ransohoff 等[8]给出了方形孔隙中无量纲曲率Cm=1.89.当孔隙收缩性较平缓时1/Rzt可以忽略,式(27)进一步简化为Rc＞1.89Rt.

Ransohoff 等[8]提出的静态准则一般适用于驱替压差较小的准静态情况.为了验证该静态准则的适用性,本文采用改进的伪势模型进行模拟验证.为了方便将模拟结果与Ransohoff 静态准则进行对比,首先定义无量纲孔喉长度比L*及孔喉宽度比R*表征孔隙结构的变化,L*和R*可以表示为

式中,Lp为整个孔喉系统的长度,L为喉道长度.基于此,可以得到卡断静态准则为R*≤0.53.

建立301 × 51 的格子空间,在120 ＜x＜ 180 及10 ＜y＜ 40 的位置处设置宽度为30 lu 的喉道(L*=0.2,R*=0.6).x方向的边界设置与2.5 节中一致,在y=0 及y=50 及喉道壁面设置为固体壁面,润湿角设置为0°,并采用反弹边界.模拟的两相流体性质为:气相密度ρv=0.53 lu、液相密度ρl=6.08 lu、气相黏度μv=0.088 lu、液相黏度μl=1.013 lu、两相界面张力为σ=0.178 lu.选取的模拟密度比(ρr=11.5)与15 MPa,350 K 温度条件下的水(980 kg/m3)与甲烷(90 kg/m3)的密度比(ρr=11)接近(NIST 化学数据库).

为了更清楚的表示驱替过程,笔者将整个驱替过程分为三个阶段(以时间为分界线).第一阶段,气相在左端的孔隙中的驱替阶段;第二阶段,气相侵入喉道的驱替阶段;第三阶段,气相突破喉道进入右端孔隙的驱替阶段.在这三个阶段中,选取气相饱和度Sg以及滞留液相饱和度Srl作为量化表征液相变化的参数.气相饱和度Sg是指在驱替过程中气相体积分数与整个区域体积分数之比;滞留液相饱和度Srl是指在驱替过程中,以两相弯液面的顶点为分界线(垂直线),将整个区域分为两部分,分界线以内的液相体积分数与整个区域体积分数的比值.换句话说,Sg+Srl等效为石油工程领域中活塞式驱替(界面为平直界面)过程中气驱波及的全部面积.捕捉了驱替过程中的Sg以及Sg+Srl随时间的变化规律以及三种阶段某一时刻下的两相流体的密度分布(图9标注).第一阶段(t＜ 9000),气相在左端的孔隙中的驱替阶段,驱替过程中形成相对稳定的弯液面,Srl保持较低的数值;第二阶段(9000≤t≤13 000),气相侵入喉道的驱替阶段,此时在孔隙角隅处及喉道壁面上存在滞留的液相,随着模拟时间步的增大Srl逐渐增大;第三阶段(t＞ 13 000),气相突破喉道进入右端孔隙的驱替阶段,此时由于壁面的强润湿性,气相无法与壁面接触,导致右端孔隙壁面上滞留大量液相,Srl迅速增大,并且孔隙右侧滞留的液相(润湿相)随时间逐渐发生回流到喉道中,最终卡断相界面.此外,所提模型孔喉条件为(R*=0.6),并不满足卡断发生的静态准则条件(R*≤ 0.53),但是通过模拟可以观察到卡断现象.这是因为在驱替的第二和第三阶段,气相无法与喉道壁面以及右端孔隙壁面接触,导致气相突破喉道后的前缘曲率半径始终小于孔隙半径,使得Cm取值偏大,从而导致基于角流推导的静态准则预测结果偏保守.

图9 不同时间步下的Sg 和Srl 变化Fig.9 Sg and Srl changes at different time steps

进一步计算了驱替过程中的第三阶段内的Pc-neck及Pc-front数值变化.如图10 所示,随着时间步的增大Pc-neck不断增大,Pc-front不断减小,这是因为Srl增大使得Pc-neck不断增大,气相突破喉道之后前缘曲率半径不断增大导致Pc-front不断减小,最终Pc-neck＞Pc-front,使得滞留的液相发生回流从而卡断相界面.

图10 不同时间步下Pc-neck 及Pc-front 随时间变化Fig.10 Pc-neck and Pc-front changes at different time steps

上述结果与之前静态准则判断发生卡断的条件一致,说明Roof[6]基于毛管压力不平衡的卡断判断机制在突变孔喉系统中同样适用.同时,上述2D 模拟结果与Wu 等[15]在3D 孔喉结构中气泡卡断实验结果基本一致,说明模拟结果对实际的3D 情况仍然有效.然而,大量研究表明,动态因素对卡断现象有较大影响[11],但是目前的理论模型都是准静态几何准则,无法考虑动态因素的影响.为此,利用提出的伪势模型,采用数值模拟方法评价毛管数Ca,L*及R*对气液两相驱替中相界面卡断的影响规律.

3.2 毛管数的影响

毛管数Ca是表征两相流动的重要参数,其定义为

式中,ρn,νn分别代表非润湿相的密度及运动黏度,在本文的模拟中与气相的密度及运动黏度相同,uave为两相驱替的平均速度[9].

通过改变驱替压差大小进而改变毛管数,研究毛管数对于卡断的影响.建立301 × 51 的格子空间,在100 ＜x＜ 200 及10 ＜y＜ 40 的位置处设置长度为100 lu,宽度为30 lu 的喉道(L*=0.33,R*=0.6),驱替压差Δp=0.069,其余条件均与3.1 中相同.每隔500 时间步,捕捉相界面前缘位置,通过相界面位置与时间关系曲线获得两相驱替平均速度(uave=Δs/Δt).需要指出的是,Δs和Δt表示发生卡断现象之前的位移和时间.图11(a)给出了Ca=4.8 × 10-3时两相驱替过程中的Sg以及Sg+Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布,第一阶段的模拟结果与3.1 节中一致;当气相进入喉道之后,孔隙角隅及喉道壁面始终存在滞留的液相,Srl逐渐增大.由于驱替压差较小,气相不能突破喉道,最终气液相界面在喉道内保持平衡不再延伸,Srl保持不变(Srl=0.11).造成上述现象的原因是气液界面的毛管压力以及固液界面的作用力形成较大的流动阻力,这种现象称为毛管“钉扎”效应[59].因此,发生卡断的首要条件即为驱替压差可以克服毛管“钉扎”效应,形成有效驱替.

增大驱替压差至Δp=0.078,图11(b)给出了Ca=5.6 × 10-3时两相驱替过程的Sg以及Sg+Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布.驱替的第一阶段与前文模拟结果一致;在气相进入喉道之后,孔隙角隅及喉道壁面始终存在滞留的液相,Srl逐渐增大,且滞留在喉道壁面上的液相厚度呈非均匀分布.随着模拟时间步的增大,滞留的液相逐渐回流并卡断相界面.值得一提的是,发生卡断时间属于驱替的第二阶段,说明在非渐变孔喉系统中卡断可以在气相突破喉道之前发生.

继续增大驱替压差至Δp=0.105,其余条件保持不变,图11(c)给出了Ca=8.0 × 10-3时两相驱替过程中的Sg以及Sg+Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布.此时毛管数相对较大,相比于之前的模拟结果,在第二阶段,喉道壁面上滞留的液相饱和度Srl较小;随着模拟的时间步的增大,气相不断向前流动,在突破喉道前并未发生卡断,直至气相到达出口端,随着模拟时间步的继续增大,滞留的液相发生回流并形成卡断.并且随着毛管数的增大,发生卡断的位置向喉道出口端移动,这与Wei 等[60]采用3D 孔喉结构研究不同毛管数条件下气泡卡断的实验结果一致.

图11 不同时间步下的Sg 和Srl 变化Fig.11 Sg and Srl changes at different time steps

进一步增大驱替压差至Δp=0.114,其余条件保持不变,图11(d)给出了Ca=8.6 × 10-3时两相驱替过程中的Sg以及Sg+Srl随时间的变化规律以及三个阶段某一时刻下的两相流体的密度分布.此时驱替压差很大,相比于之前的模拟结果Srl较小;随着时间步的增大,气相从喉道中突破,并最终到达出口端.在气相到达出口端以后,随着模拟时间步的增大,并未发生卡断现象.这说明对于固定的孔喉结构,存在一个临界毛管数,当系统的毛管数大于临界毛管数时,将会抑制滞留的液相在喉道内发生回流,从而抑制卡断现象的发生.

图12 给出了50 000 时间步内,上述四种条件下相界面前缘顶点位置随时间的变化.发生卡断之后,把回退到喉道中的弯液面顶点作为实际的相界面位置,这样当相界面位置突然变化时,即可判断为发生卡断.可以看出,气相为连续相时,驱替过程中可以发生连续卡断,并且每次发生卡断的位置几乎不会变,该过程类似于液滴/气泡的形成过程[61-62].同时,对于特定的孔喉系统,只有在满足一定的毛管数范围内才会发生卡断,较大的毛管数会抑制驱替过程中滞留液相的回流从而抑制卡断的发生,较小的毛管数无法克服毛管“钉扎”效应形成有效驱替.模拟结果与文献[9]在渐变孔喉中研究气泡在压力梯度作用下的卡断行为所得出的结论一致.

图12 相界面位置随时间的变化Fig.12 Change of phase interface position with time

3.3 孔喉长度比的影响

孔喉结构对卡断同样具有重要影响[9-10],为此本节利用提出的伪势模型研究喉道长度对卡断的影响.建立301 × 51 的格子空间,在10 ＜y＜ 40 的位置处设置相同宽度(R*=0.6),不同长度的喉道来研究不同L*对卡断的影响.对于不同的喉道长度,驱替压差比毛管数更能体现动态因素对卡断的影响,因此采用驱替压差代替毛管数表征动态参数的影响.

在模拟过程中,改变驱替压差从0.036 至0.123,改变L*从0.08 至0.4 来确定卡断发生的临界条件;这里仅用Δp=0.084,L*=0.3 和Δp=0.084,L*=0.2 两种模拟条件为例,展示不同孔喉长度比的驱替效果.图13 给出了Δp=0.084,L*=0.3 及L*=0.2 条件下,不同时间步的驱替过程.随着喉道长度的减小,流动阻力减小,两相驱替平均速度增大,驱替过程中滞留在孔隙角隅及喉道壁面处的液相饱和度Srl较低.值得一提的是,设置的孔喉宽度比R*=0.6,不满足静态准则预测发生卡断的条件(R*≤0.53),但是模拟结果可以观察到卡断现象,说明即使不满足静态准则预测发生卡断的孔喉宽度比,足够长的喉道长度也会促进卡断的发生.因此,在相同驱替压差的作用下,不同喉道长度的孔喉系统会呈现不同的流动状态.

图13 两相驱替过程对比Fig.13 Comparison of two-phase displacement processes

为了确定不同孔喉长度比的卡断发生条件,开展了一系列模拟,得到不同孔喉长度比下发生卡断的临界驱替压差,如图14 所示.对于固定的喉道宽度的孔喉系统,存在临界孔喉长度比(L*=0.08),使得无论驱替压差如何变化,气液两相驱替过程中不会出现卡断现象,这与Deng 等[10]提出的短波长的渐变孔喉中卡断会被抑制的结论一致.随着喉道长度增大,不同的驱替压差会使得气液驱替过程中出现卡断,连续流动以及无效驱替三种流动状态.当驱替压差小于临界驱替压差下界时,无法克服毛管“钉扎”效应,形成无效驱替;当驱替压差大于临界驱替压差上界时,气相将保持连续流动的状态,抑制卡断的发生.

图14 不同L*条件下气液流动状态与驱替压差的关系Fig.14 The relationship between gas-liquid flow state and displacement pressure difference with different L*

3.4 孔喉宽度比的影响

除了喉道长度,喉道的宽度对于卡断现象也有一定影响,本节通过固定L*改变R*,研究不同喉道宽度对卡断的影响.建立301 × 51 的格子空间,在125 ＜x＜ 175 设置长度为50 lu 的喉道(L*=0.167),其余参数设置与3.3 节相同.在模拟过程中,改变驱替压差从0.033 至0.1,改变R*大小从0.52 至0.68 来确定卡断临界条件;这里仅用Δp=0.084,R*=0.52 和Δp=0.084,R*=0.6 两组条件为例,展示不同孔喉宽度比对气液两相驱替过程中的流动状态的影响.图15 给出了Δp=0.084,R*=0.52 及R*=0.6 两种条件下的驱替过程.随着喉道宽度的减小,流动阻力显著增大,两相驱替平均速度减小,驱替过程中滞留在孔隙角隅及喉道壁面处的液相饱和度Srl较高.值得一提的是,R*=0.52 满足静态准则预测发生卡断的条件(R*≤ 0.53),模拟结果同样观察到了卡断现象,验证了模拟结果的准确性.对于R*=0.6 的情况下,驱替过程中始终不会发生卡断现象,这也说明了喉道宽度增大将会抑制卡断的发生.

图15 两相驱替过程对比Fig.15 Comparison of two-phase displacement processes

为了确定不同孔喉宽度比对卡断现象的影响,开展了一系列模拟,得到了不同孔喉宽度比下发生卡断的临界驱替压差.如图16 所示,与3.3 节类似,对于固定喉道长度的孔喉系统,气液两相在不同驱替压差的作用下形成气相连续流动,卡断及无效驱替3 类流动状态.随着喉道宽度的增加,发生卡断的驱替压差范围逐渐减小,当孔喉宽度比增大至R*=0.68 时,无论施加多大的驱替压差,喉道内都不会发生卡断现象.值得一提的是,对于R*=0.52 的孔喉结构(满足静态准则预测发生卡断的条件),当驱替压差大于0.1 时,驱替过程中将不会发生卡断现象,这说明即使达到静态准则所预测的卡断条件,较大的驱压差也可以抑制卡断的发生.同时,基于等式(27)预测临界孔喉宽度比为(R*=0.53),相比于通过模拟得出的临界孔喉宽度比(R*=0.68)降低28.3%.

图16 不同R*条件下气液流动状态与驱替压差的关系Fig.16 The relationship between gas-liquid flow state and displacement pressure difference with different R*

4 结论

本文在原始单组分伪势格子玻尔兹曼模型基础上改进流体-流体作用力格式,建立了改进的伪势模型,并在此基础上模拟了孔喉系统中的气液两相在不同驱替条件下的流动状态,得出以下结论.

(1) 明确了非渐变孔喉结构卡断机理,即孔喉毛管压力不平衡使得驱替过程中滞留的液相(润湿相)随时间逐渐发生回流是造成相界面卡断的根本原因.但是,模拟过程中在孔隙角隅及喉道壁面观察到大量滞留的液相,这导致基于角流的假设思想建立的静态卡断判断准则将会高估喉道右侧气泡的曲率半径从而低估卡断发生的条件.

(2) 揭示了动态因素驱替压差(毛管数)对气液两相驱替的影响规律.在非渐变孔喉系统中,只有 驱替压差(毛管数)大小在一定范围时,喉道内才会发生卡断;超过该毛管数上界,即使满足静态准则条件,卡断也会被抑制;低于该下界,无法完成驱替;同时毛管数的增大使得发生卡断的位置向喉道出口端移动.

(3) 揭示了孔喉长度比对气液两相驱替的影响规律.对于固定宽度的孔喉系统,即使不满足静态卡断发生的孔喉宽度比(R*≤ 0.53),足够长的喉道也可以促进卡断的发生,随着喉道长度增大,发生卡断的驱替压差范围越大.并且存在一个临界喉道长度,使得无论驱替压差如何变化,喉道内都不会发生卡断,对于本文的模型,临界孔喉长度比L*=0.08.

(4) 揭示了孔喉宽度比对气液两相驱替的影响规律.对于固定长度的孔喉系统,喉道宽度越大,发生卡断的驱替压差范围越小;并且存在一个临界喉道宽度,使得无论驱替压差如何变化,喉道内都不会发生卡断.对于本文的模型,临界孔喉宽度比R*=0.68.若采用静态准则预测该模型发生卡断的条件(R*=0.53),将会低估28.3%.