基于孔隙网络模型的非牛顿流体驱替模式研究

2023-12-04杨鑫李星甫唐雁冰李闽BernabYves李晨曦赵金洲杜翔宇

地球物理学报 2023年12期

杨鑫, 李星甫, 唐雁冰*, 李闽, Bernabé Yves,李晨曦, 赵金洲, 杜翔宇

1 西南石油大学油气藏地质及开发工程国家重点实验室, 成都 6105002 美国麻省理工学院地球、大气和行星科学系, 剑桥 021383 日本早稻田大学地球科学、资源与环境工程系, 东京 1698555

0 引言

研究多孔介质中的流体流动规律对提高地下油气采收率(Bikkina et al., 2016)、清理地下污染水层(Ji et al., 2008)和二氧化碳地质封存(Bakhshian et al., 2020;吕鹏飞等,2023)等许多工程技术都具有重要意义.随注入速度增加,一种流体驱替另一种黏度更大的流体时表现为毛管指进向黏性指进转换,并且两者之间存在一个大范围、低驱替效率的过渡区(Rodríguez de Castro et al., 2015).指进是毛管力、黏性力、重力和岩石表面润湿性等因素的复杂相互作用导致的优势通道或指状流动(Detwiler et al., 2009),可以被描述为不同流体之间相互作用的能量和无序系统的约束之间的竞争(Dias et al., 2012).指进常常是一种不利的现象,例如在石油工业中,当人们试图通过将水泵入油藏(即注水驱替)来置换原油时,指进会使得注入流体绕过和捕集原油,导致较早的突破(Chen et al., 2017).突破时的流体驱替效率常常会影响最终驱替效果(Xu et al., 2014),尤其是在非均质性很强的储层中,突破后很难再通过注水的方式大幅度提高原油采收率(Mosavat and Torabi, 2016).

流体驱替模式受孔喉微观结构,流体与流体或流体与岩石的耦合机制影响,预测注入流体的驱替模式是十分困难的(Bultreys et al., 2016).Lenormand等(1988)用简单的二维多孔介质模型进行实验,绘制了毛管数Ca和两相流体黏度比M的经典相图,划分了紧凑型驱替(CD)、毛管指进(CF)、黏性指进(VF)和过渡区(CZ)等不同驱替模式所处的区域.不同驱替模式的驱替效率存在差异,研究人员可以根据相图判断不同流速条件下的驱替模式,从而调整合理的驱替速度来达到提高原油采收率的目的.研究人员利用实验和数值模拟方法,结合分形理论等手段建立了不同条件下和经典相图类似的相图,并分析了相图中不同驱替模式的影响因素(Cottin et al., 2010; Yang et al., 2019).例如,在Lenormand等(1988)的基础上,Chen等(2017)使用二维平板裂缝物理模型进行水驱油实验,将经典相图的适用范围拓展至适用于复杂多孔介质.Holtzman和Segre(2015)、Zhao等(2016)和Lan等(2020)分析了不同注入毛管数情况下润湿角θ对驱替模式的影响,并提出了用于区分不同驱替模式范围的理论公式,从而建立了[Ca,θ]空间中发生毛管指进、黏性指进和过渡区的相图.同时,储层岩石孔喉非均质性也会对驱替模式产生影响(Toussaint et al., 2005),使得毛管指进和黏性指进在相图中所处的区域和范围发生变化.

上述研究中采用的被驱替相均为牛顿流体,而地下原油常常表现出非牛顿流体的特性(即由剪切应力引起的黏度变化)(Wang et al., 2019;刘文山等,2021).Shah和 Yortsos(1993)的研究指出非牛顿流体(如Ellis流体和幂律流体)驱替也存在一个相图,其中毛管指进,黏性指进,过渡区所处的区域均随剪切流变性改变而发生不同程度的移动,但没有明确给出不同驱替模式之间的分界线.Tian和Yao(1999)的研究表明,非牛顿流体的剪切流变性会影响驱替过程,但他们并没有研究毛管指进、黏性指进和过渡区的特征,也没有比较牛顿流体和非牛顿流体性质差异对不同驱替模式的影响.此外,储层流体是可压缩的,现有的准静态和动态网络模拟技术是假设流量与压力传递过程只发生在孔隙网络模型的喉道上,节点间的流动满足基尔霍夫定律(Tang et al., 2020;王猛等,2022),多用于研究较小尺度模型中的油水两相流(Zhao et al., 2021; Yang et al., 2022),未考虑油水的压缩性,故而忽略了当模型较大时,压力波传播和时间等因素对驱替过程的影响(Li et al., 2017).因此,此研究在以前的孔隙网络模型中引入适用于可压缩流体的非稳态渗流理论,即通过将稳态渗流方法(Tang et al., 2021)中的拉普拉斯方程转换为与压力波传播和时间相关的非稳态渗流的压力扩散方程,形成了非稳态油水两相孔隙网络模拟方法.

本文提出了一种考虑了流体可压缩性的非稳态孔隙网络模型,并在模型中增加了非牛顿流体的剪切流变性,用以模拟跨越多个数量级的毛管数以及不同黏度比条件下的水驱油过程.分析了毛管指进、黏性指进和过渡区各自的特征,并采用特征前缘流量对各区域进行量化.最后比较了被驱替相为牛顿流体和非牛顿流体时不同驱替模式的差异.这项工作为深入了解地下复杂流体黏性力和毛管力的作用机制提供了新思路,在石油工程应用中具有实际意义.

1 数值模拟方法

1.1 构建孔隙网络模型

动态网络模拟方法的基本原理是假设能量耗散或熵产生过程(如黏性流与扩散混合等流体流动行为)只发生在孔隙网络模型的喉道上,孔隙节点流量遵守质量守恒定律,因此孔隙网络模型可以看作是真实多孔介质的简化(David, 1993).而孔隙网络模拟器结合数学、物理控制方程,可以合理地模拟再现流体在孔喉内的流动(Blunt et al., 2002; Xiong et al., 2020).

本文根据Bernabé等(2011)和Li等(2015)研究中构建孔隙网络模型的方法建立了圆盘状模型,模型是由二维正方形网格组成的网络(Seeburger and Nur, 1984),如图1a所示.模型直径方向节点数为301,每个节点i与相邻四个节点j(j=1,2,3,4)由4个不同半径的管道连接.管道长度l=300 μm,半径r按对数均匀分布(如图1b)随机赋值,最小值和最大值分别为1.35 μm和78.65 μm,半径分布标准差σ=20 μm,平均值〈r〉=19 μm,平均水力半径rH=40 μm,赋值后的管道半径均恒定不变.计算得到模型半径R=4.5 cm,总孔隙体积Vp=0.1×10-6m3.通过上述步骤设置的模型,网络节点为孔隙,管道为喉道,其他部分均视为被固体颗粒充填.该模型常被用于研究非混相驱替过程中的黏性力与毛管力对流动的影响以及渗流过程中的临界路径和优势通道等现象(Tang et al., 2019, 2020).

图1 (a) 多孔介质中水驱油示意图; (b) 孔隙网络模型中喉道半径分布频率和累计曲线图

1.2 非稳态两相流动

模型初始被润湿相流体(原油,黏度为μo)充分饱和,模拟过程中,将非润湿相流体(水,黏度为μw)以恒定的流量qin从中心注入,保持模型边界出口端的压力恒定为大气压(0.1 MPa),使注入流体由模型中心沿网络向四周边界流动(如图1a所示).当两相流体存在同一喉道中时,存在如下假设条件:(1)所有的流体都被认为包含在喉道和孔隙中,流体压降只存在于孔隙之间的喉道中;(2)孔隙网络喉道中的两相流体间只存在一个界面(无扩散);(3)喉道中只有活塞式驱替(Zhao et al., 2016).两相流体界面引起的毛管压力pcij=-2γcosθ/rij(Teige et al., 2006),其中γ为界面张力(20×mN/m),θ为润湿角(150°).利用扩展Hagen-Poiseuille方程(Aker et al., 1998)计算喉道内体积流量qij:

qij=-gijeff(Δpij-pcij),

(1)

多孔介质内部流体是可压缩的,每个控制体内部的中心孔隙i与相邻孔隙j(假设j的编号为1,2,3,4)之间的流动满足(推导见附录A):

(2)

其中Vp为孔隙i所代表的控制体总体积(控制体由孔隙i与相邻四个长度一半的喉道组成),Δpi/Δt表示控制体的平均压力随时间的变化,Cρeff表示混合流体的有效压缩系数,其值与具体流体的性质有关,受温度和压力等因素的影响,本文研究对象为油水,设置为1×10-10Pa-1(He, 1994).用非稳态孔隙网络模型模拟了水驱油(牛顿流体),流体突破时的驱替形态和分形维数与Chen和Wilkinson(1985)的实验与模拟结果吻合性很好(见附录B).

(3)

参照Sadowski(1963)的对Ellis模型的实验研究,τ1/2和α取值分别为3.5 Pa和2.4.

数值模拟计算过程中,式(1)中毛管压力项可写为:

(4)

(5)

从而将式(2)改写为:

(6)

对式(6)进行离散处理,同时增加汇源项Qi:

(7)

其中,t和t+Δt分别对应流动状态变化前后的时刻,式(7)可写为:

(8)

(9a)

(9b)

在给定边界条件时,未知数为各个孔隙处的流体压力.通过上述步骤迭代求解更新压力场,从而利用孔隙与孔隙间的压力计算每个喉道中的流量,并通过驱替流体完全占据的喉道体积计算当前时刻的饱和度.在每个时间步内,由于非牛顿流体黏度随压差变化,且两相流体界面实时移动,每个时间步长内都需要更新传导率,对计算机算力有很高的要求.本文使用GPU加速算法(Naumov et al., 2015)进行模拟以加快模拟速度.

采用注入流体的毛管数Ca和黏度比M这两个无量纲数来表征流体驱替模式,并用以比较黏性力和毛管力的大小(Reynolds et al., 2017).Ca=uμw/γ(u=qin/Ad,Ad=4π〈r〉2),Ca的范围设置为lgCa=-8.5～lgCa=-2.5;M=μw/μo,μw=1 mPa·s,μo被分别设置为1000、500、300、200、100、50 mPa·s.

2 结果与分析

2.1 流体驱替模式

模拟结果可观察到不同黏度比M和大范围变化的毛管数Ca在注入流体突破边界时的流体驱替模式,如图2a所示,用不同颜色表示注入流体占据的不同孔喉.

图2 (a)不同黏度比M和毛管数Ca情况下注水突破时的驱替模式示意图; (b) 毛管指进示意图(lgCa=-8.5); (c) 过渡区示意图(lgCa=-3.5); (d) 黏性指进示意图(lgCa=-2.5)

毛管力与黏性阻力的相互作用机制会导致水驱油时产生不同形态的手指(Martys et al., 1991).低流速条件下,毛管力占据主导作用,两种流体的黏性阻力很小,如图2b为lgCa=-8.5,M=1/50时的侵入剖面图,毛管指进会产生较粗形态的手指,手指沿各个路径方向生长(横向生长(标记1),向入口处生长(标记2),形成圈闭(标记3)捕集原油,圈闭数量和大小与网格非均质性和模型尺度有关(Lenormand et al., 1988)).而在高流速条件下(lgCa>-3.5),原油黏性阻力占据主导作用,高注入压力使得驱替流体在多个流动路径中产生大量的“细手指状”的指进,并向出口处不断延伸,驱替模式表现为黏性指进(如图2d).黏性手指形态细窄,呈树状,且难以见到树状手指向入口处生长的现象.随着Ca由lgCa=-7.5逐渐增加到-3.5,黏性力的影响逐渐增加.如lgCa=-4.5时,黏性力与毛管力的作用相当,在这种状态下,多孔介质中同时发生毛管指进和黏性指进,流动路径显著减少(如图2c),大量原油被滞留,波及效率(Holtzman and Segre, 2015)降低,这是流动状态由毛管力主导向黏性力主导转变的过渡区.

2.2 过渡区的特征与划分依据

图3 计算分形维数示意图

前人的研究表明(Lan et al., 2020),毛管指进、黏性指进和过渡区处的Sw和D的大小均不一样,过渡区中的Sw和D最小(Singh and Mohanty, 2003),驱油效率最低.通过模拟结果绘制了Sw和D随毛管数Ca和黏度比M变化的关系图(分别对应图4a和图4b).结果表明,不同M时,随Ca增加,D和Sw两者都先减小后增加,在中间Ca处的值最小,这也表明非稳态孔隙网络模型能很好地模拟从毛管指进到黏性指进的过渡区.

图4 (a)、(b) 不同黏度比lgM和不同毛管数lgCa时驱替流体突破时的饱和度Sw和分形维数D的变化情况; (c) 饱和度Sw和分形维数D关系图

使用不同材料与不同结构的物理模型或注入不同类型的流体介质,也会使得过渡区处的分形维数和饱和度产生差异.例如,Wang等(2013)在多孔介质中利用超临界二氧化碳驱替盐水,过渡区中,注入流体饱和度小于0.30.Zhao等(2016)认为较高的Ca时(如图4a过渡区左边,lgCa>-4),D>1.62为黏性指进;较低的Ca时(如图4a过渡区右边,lgCa<-6),D>1.82为毛管指进.Chen等(2017)的研究中,过渡区处的D<1.62,Sw<0.25.由图4a和图4b可以得到如图4c所示的不同Ca和M情况下D和Sw关系图,表明D随Sw增加而增加.结合图3和图4c,本文将分形维数D<1.61时的驱替模式划分为过渡区(图4c中阴影部分),划分依据与Holtzman和Juanes(2010)的研究一致.过渡区内,饱和度Sw<0.10,该值更小是因为本文建立的模型非均质性更强.阴影部分外的其他值均落在Måløy等(1985)和King(1987)划分的二维多孔介质中经典毛管指进和黏性指进的各自值范围内.将划分的不同驱替模式对应于图3a中,红色框区域以下部分为毛管指进区,红色框区域内为过渡区,红色框区域以上部分为黏性指进区.

图5为Lenormand等(1988)提出的经典相图,灰色区域表示毛管指进、黏性指进和紧凑型驱替,空白处为这三者之间的过渡区.图中虚线部分为Zhang等(2011)提出的不同驱替模式之间的分界线,其过渡区范围更小.结合分形维数和饱和度的大小,将模拟结果得到的毛管指进(红色圆圈)、过渡区(黑色菱形)和黏性指进(蓝色方形)对应到相图的位置并作比较,本文绝大部分预测给出的驱替模式与他们的相图显示出了良好的一致性.Lenormand等(1988)证明了过渡区的点满足公式lgCa=lgM+c,其中c是与多孔介质和流体性质有关的系数.过渡区处Sw最小,因此定义临界毛管数Ca*为每个M情况下最低Sw时对应的毛管数.每个M情况下可以计算得到一个c值,利用6组(Ca*,M)(图5中的菱形点中标有红色线)得到的不同的c值后取平均值,得到分界线(在图5中用红色线表示,其下部分表示毛管指进向过渡区的转变,其上部分表示过渡区向黏性指进转变)表达式:

图5 lgCa-lgM平面相图中不同驱替模式所处的位置

log10Ca*=log10M-2.3.

(10)

分界线斜率为1是由于使用了驱替流体定义的毛管数.在注入速度恒定的情况下,更大黏度的被驱替相会产生更大的压降,也可以使用基于压力梯度(Δp/R)的毛管数来代表不同的驱替模式(Lenormand et al., 1988),如Chatzis和Morrow(1984)提出的Ca=(kw/γ)⋅(Δp/R),其中kw为驱替相渗透率.但为了与经典相图作比较,本文采用了保持统一定义的Ca来绘制图表,这在Chen等(2018)的实验研究中已经得到了应用.

2.3 不同驱替模式的饱和度和压力特征

图6 驱替流体饱和度Sw与归一化前缘位置演化关系曲线

数值模拟结果也提供了流体突破时的压力分布信息.如图7所示,显示了使用黏度比为1/200时,毛管指进(图7a)、过渡区(图7b)和黏性指进(图7c)的压力分布.当毛管数非常小时(图7a,lgCa=-8.5),压力分布均匀且较低(不同的色条刻度表示不同大小的压力),压力梯度很小.参照压力分布色带,随毛管数增加,压力梯度增加(图7b和图7c).为观察流体压力在整个驱替过程中的演变,本文也绘制了注入压力Pin随归一化驱替相饱和度Swin/Swmax(Swin为实时变化的饱和度,Swmax为最大饱和度)增加而变化的关系图.在低毛管数情况下(图7d),初始入口压力的波动最大,这是因为以恒定速度注入时,必须不断调整节点间的压差以满足克服孔径差异引起的毛管压力的变化.随毛管数增加,毛管力的影响减弱,节点间的压差调整频率降低,入口压力仅会出现偶尔的波动(图7e),此时黏性阻力的主导优势逐步显现,达到可以和毛管力比较的级别,驱替流体需要同时克服毛管力和较大的黏性阻力.在高毛管数情况下,两相流体快速流动会引起黏性阻力显著增加,需要较高的初始入口压力以满足流动条件.此外,随饱和度增加,注入压力逐渐降低,快突破时,入口压力快速下降(图7f).这是因为随两相界面推进,流体前缘到出口的距离逐渐缩短,被驱替相饱和度降低,黏性阻力急剧减小.因此,还可以通过判断压差的波动规律和入口压力是否大幅下降来定性区分不同的驱替模式.

图7 (a)、(b)、(c) 毛管指进,过渡区和黏性指进的压力分布; (d)、(e)、(f) 毛管指进,过渡区和黏性指进的中心注入压力Pin和归一化饱和度Swin/Swmax的演化关系; (g) 平均注入压力Pm和不同毛管数lgCa的关系

图7g为整个模拟过程中不同Ca和M情况下注入压力的平均值Pm.相同Ca时,被驱替相黏度越大,黏性阻力越大,注入压力越高.如lgCa=-2.5时,M=1/1000时的入口压力明显高于M=1/50时的入口压力(6.5倍).高Pm导致中间Ca至高Ca情况下M=1/1000时的部分驱替效率比M=1/50时高,这与前人(Chen et al., 2017, 2018)在二维物理模型实验中观察到的现象一致.

2.4 驱替流体在多孔介质中的占比

图8为M=1/500,不同Ca时,侵入半径小于平均半径的总数(小喉道)和大于平均半径的总数(大喉道)与占据总侵入喉道数比例.当较小的孔喉充满高黏度原油时,注入压力难以克服小孔喉中的毛管力与黏性阻力,而更倾向于驱替阻力更小的大孔喉中的原油,故被驱替的路径中大孔喉占比更多.由图3可以看出,过渡区流体驱替路径中的孔喉数量比毛管指进和黏性指进少,波及区域更小.而图8表明,虽然过渡区的驱替路径少,但是路径中的大孔喉占比更多,这是因为随Ca增加,所有孔喉中的黏性力阻力增加,除小孔喉外,注入压力也无法克服中等孔喉中的毛管力与黏性力,流体更倾向于进入注入压力与阻力大小更接近的更大孔喉.同时,本文也绘制了突破时的饱和度Sw与所占据喉道的平均半径〈rH〉的关系图(如图9所示),也表明〈rH〉在过渡区处(图中灰色区域)具有较大值.

图8 不同Ca时小喉道和大喉道占据总侵入喉道数的比例分布柱状图

图9 平均侵入半径〈rH〉与饱和度Sw的关系图

2.5 驱替前缘的特征流量

(11)

其中t*是标定时间,表示为t*=tqin/(AdR),t是注入时间,突破前的饱和度Sw与t的关系为Sw=tqin/Vp.

(12)

则由式(11)可得到:

(13)

(14a)

(14b)

(15)

式(15)表明,特征前缘流量与划分的时间间隔与圆环数量无关,而与多孔介质结构参数(即Ad,R,rH和Vp)相关,其值随饱和度减小而增加,随侵入平均半径的增加而增加.

图10 (a) 平均的特征前缘流量和毛管数lgCa关系曲线;(b)和lg(Ca/Ca*)关系曲线

3 讨论

图11 (a)不同毛管数和黏度比时水驱牛顿流体剖面图; (b) 不同流体注入压力Pin与归一化饱和度Swin/Swmax演化关系; (c) 牛顿流体与非牛顿流体过渡区比较图

在低渗储层中,注水驱替黏度普遍较高的具有非牛顿流体性质的原油需要极大的注入压力,对实验设备要求高,难度大.本文通过对比前人的实验结果,验证了非稳态孔隙网络模型能用于水驱牛顿流体研究,并通过牛顿流体与非牛顿流体的注入压力差异也表明了该模型能体现剪切流变性对驱替过程的影响,能解释二者不同驱替模式所处范围差异的原因.

4 结论

(1)非稳态孔隙网络模型能够模拟再现由毛管力和黏性力的作用机制而产生的不同驱替模式,且不同驱替模式所处的区域与经典相图给出的各区域吻合.

(2)本文提出的特征前缘流量可以通过饱和度和平均侵入半径计算得到,能用于区分不同黏度比与毛管数下的毛管指进、黏性指进和二者之间的过渡区.

(3)非牛顿流体具有比牛顿流体更宽的过渡区.

致谢感谢审稿人的宝贵意见和建议.

附录A 两相非稳态模型

流体具有压缩性,其密度ρ(kg·m-3)随压力p(Pa)变化满足:

ρ=ρ0eCρ(p-p0),

(A1)

其中,p0和ρ0分别为初始压力和密度,其值会在下一个时间步更新;Cρ为流体的压缩系数,单位为Pa-1.将公式(A1)按麦克劳林级数展开并取前两项可得ρ≈ρ0[1+Cρ(p-p0)].

图A1 控制体示意图

流体在控制体中的流动过程遵循质量守恒定律,即:

(A2)

(A3)

在控制体中:

(A4)

则控制体内任意一个方向:

(A5)

所以:

(A6)

其中〈pij〉=(pi+pj)/2,表示连接相邻节点的管道中间压力,pj为节点j的压力.在同一管道中:

(A7)

由式(A2)、(A4)、(A6)、(A7)得:

(A8)

控制体内两相流体流动同样满足公式(A2).任意圆管内,注入流体I和被驱替流体II的混合流体的密度计算式为:

(A9)

任意节点i为中心的控制体内,混合流体的密度计算式为:

(A10)

其中,ρI为注入流体I的密度,ρII为被驱替流体II的密度,VijI和VijII分别为节点i和j之间孔喉通道内注入流体和被驱替流体的体积,ViI和ViII为以节点i为中心的控制体内注入流体和被驱替流体的体积.由于油水密度差异较小,为简化计算过程,认为以节点i为中心的控制体内混合流体的密度与节点i和j之间孔喉通道内混合流体密度相近,则控制体内混合流体的有效密度ρeff满足:

ρij≈ρi=ρeff=ρISI+ρIISII,

(A11)

其中,SI和SII分别为以节点i为中心的控制体内注入流体和被驱替流体的饱和度.任意控制体内混合流体的压缩系数满足:

(A12a)

ρeff=ρeff0eCeffρ(p-p0),

(A12b)

Ceffρ=CISI+CIISII,

(A12c)

其中,Ceffρ为混合流体的有效压缩系数,ρeff0为初始条件下混合流体的密度,CI和CII分别为注入流体和被驱替流体的压缩系数.与非稳态单相渗流模型推导中采用的方法类似,可以得到两相渗流数学模型:

(A13)

其中,gijeff为两相流体的等效传导率,pcij为节点i和j之间孔喉中两相流体的毛管力,μeff为节点i和j之间混合流体的有效黏度.

公式(A13)可适用于二维三角形、正方形和六边形网络模型以及三维的SC、FCC、BCC等规则网络模型以及其他复杂不规则网络模型.