刘 静, 郑黎明, 李悦静, 夏军勇, 蒲春生

(1.中国石油大学(华东) 石油工程学院， 山东 青岛 266580； 2.非常规油气开发教育部重点实验室(中国石油大学(华东))， 山东 青岛 266580; 3.燕山大学 车辆与能源学院， 河北 秦皇岛 066004)

0 引言

低频波动采油技术，是一种利用物理场激励油层达到提高原油采收率目的的增产技术[1,2].揭示低频弹性波作用下储层渗流动力学变化机制是该技术应用和推广的基础.然而，关于低频波动采油机理的系列研究主要集中于实验研究和半定量简化动力学模型分析[3-6].实验模拟主要定性揭示了储层物性受振动参数的影响敏感性，无法定量解释波动对储层渗流过程的影响规律.理论模型研究方面，部分研究人员通过建立半定量简化波动渗流模型(将波动力直接叠加到渗流运动方程)或者简化为单一毛管，研究波动对粘附层、液滴脱落等的影响[7-10].关于初始时刻渗流场与波场的耦合作用问题，研究人员也开展了部分理论分析工作，但主要针对一维、二维单相流体的作用过程，对于油水两相渗流场与波场耦合模型研究目前较少[11,12].

基于此，借鉴经典的孔隙介质弹性波传播理论模型，并引入目前低频波动采油技术已有实验机理认识，联立波动条件下的流体运动方程、固体运动方程、连续性方程和状态方程，利用COMSOL数值模拟软件，建立物理场作用下油水两相渗流的物理模型和数学模型，开展波动耦合渗流规律的研究，揭示低频波动对采油过程的影响及振动参数影响的敏感性.

1 数学模型

低频波动采油对实际地层渗流的影响较为复杂，为研究方便，本文模拟中假设储层岩石满足线弹性变化关系，假设低频采油激振力较小，不引起岩石本身弹塑性质的变化，仍采用线性本构关系.

1.1 运动方程和连续方程

油水渗流过程中，不仅液相在多孔介质中发生运动，固相岩石受到载荷作用也会发生形变，从而产生运动.初始时刻流体流动不考虑惯性作用和流固耦合作用，结合Darcy定律，基于均匀、各向同性、弹性多孔介质的经典Biot弹性波传播模型，考虑流固耦合的固相、液相运动方程为[13-15]：

(1)

固相、液相的连续性方程为：

(2)

式(1)、(2)中：下标w表示水相，o表示油相；τ为作用在单元体上的总应力，MPa；u为固体位移，m；U为水相对于固体的相对位移，U=Ui+Uv，Ui为水相初始宏观位移，Uv为水相波诱导位移，m；W为油相对于固体的相对位移，W=Wi+Wv，Wi为油相初始宏观位移，Wv为油相波诱导位移，m；pw为水相压力，MPa，po为油相压力，MPa；ρw为水相密度，kg/m3；ρo为油相密度，kg/m3；ρ为多孔介质平均密度，kg/m3，ρd为固体颗粒密度，kg/m3；ηw为水相粘度，mPa·s；ηo为油相粘度，mPa·s；kw为水相渗透率，μm2；ko为油相渗透率，μm2；sw为水相饱和度；so为油相饱和度；φ为孔隙度.

1.2 状态方程

弹性波作用下流体密度随压力变化而变化，在渗流力学中常用流体压缩系数和压力变化量的函数表示，此文中用体积模量来表示液体压缩方程.

(3)

式(3)中：Kw为水相体积模量，MPa；Ko为油相的体积模量，MPa.

考虑固体压缩性，用孔隙压力引起固体密度改变和有效应力引起的固体密度变化之和来表示固体压缩性[13].

(4)

式(4)中：φ为孔隙度，εv为固体骨架体应变，εv=·u；α为Biot系数；1/Qc为Biot引入的假设下一截面的代表流固间耦合关系的系数，为正值，其中K为岩石体积模量，MPa；Kdry为干燥固体骨架体积模量，MPa.

其中K的取值可利用斑块饱和模型计算得到，岩石骨架体积模量和流体模的关系可表示为：

(5)

式(5)中：Ki(i=w，o)为第i相流体的体积模量，MPa；Kmi为单独饱和某一相流体时的岩石体积模量，MPa.

岩石骨架的变形以及流体压力的变化会引起孔隙度发生改变，根据Biot理论传播模型，孔隙度、形变和孔隙压力的关系式可表示为[14,15]：

(6)

式(6)中：P为油水两相渗流系统平均孔隙压力，用体积平均法求得p=swpw+sopo，MPa.

波动下粘度变化方程[16]：

ηo=ηo(ω)≈ηiF(ω)

(7)

孔隙变化方程：

(8)

式(8)中：k为绝对渗透率，μm2；k0为初始绝对渗透率，μm2；φ0为初始孔隙度.

相渗模型采用Mualem理论模型，毛管力曲线采用以下方程：

(9)

式(9)中：m、n、l为van Genuchten毛管力曲线隐式方程的拟合系数.

1.3 定解条件

1.3.1 初始压力

在该一维储层中，弹性波作用之前，油水处于渗流平衡状态.在单相流动区，可根据液柱重力梯度算出初始压力分布状况.各相的压力遵循单相流动时的计算方法，同一点不同流相间的压力差为毛管力.初始时刻压力可表示为：

pf(x,0)=pe

(10)

式(10)中：pe为油藏压力，MPa.

1.3.2 初始饱和度

对于多相渗流问题，需要指定油藏的初始饱和度分布，在单相流动区，饱和度为定值；在油水两相区，各相饱和度通过联立运动方程、连续性方程、状态方程等的波动耦合渗流模型求得.初始时刻饱和度可表示为：

sw(x,0)=swc，so(x,0)=1-swc

(11)

式(11)中：swc为束缚水饱和度.

1.3.3 初始位移和速度

假设研究的初始时刻，固相位移为0，流体宏观渗流位移也为0，但在弹性波作用之后的渗流过程中同时存在宏观渗流和波诱导渗流，初始时刻速度不等于0.

流固初始渗流位移为：

u(x,0)=0,U(x,0)=W(x,0)=0

(12)

初始渗流速度为：

(13a)

(13b)

1.4 边界条件

弹性传播方向与研究储层渗流方向一致时，油藏边界包括顶部、底部、注入井边界和采出井边界四部分，其中，低频波动边界可与注入或采出边界条件重叠.

油藏边界处水相压力等于注入端压力，油相压力等于采出端压力：

pw(0,t)=pin,po(L,t)=pout

(14)

震源边界采用最简单的位移边界表示，本文中用谐函数形式表示如下：

u(0,t)=u0(x,t)=u0exp(iωt)

(15)

式(15)中：u0为初始振幅.

无限远或强制位移边界(数值模拟过程中引入吸收边界条件)处：

u(∞,t)=0

(16)

2 模型求解

2.1 物理模型建立

利用COMSOL进行有限元模拟，基于控制方程分别构建固体运动方程、流体运动方程、压力-位移方程(水相和油相压力方程)和饱和度-位移方程的系数型偏微分方程形式，一维渗流模拟时，忽略z方向的控制方程分量，并建立弹性波一维驱替物理模型，弹性波作用下一维驱替模型中，保证弹性波传播的主方向与原渗流场方向一致，且沿x方向，如图1所示.

在研究区域纵向外围各拓展L×d/4，右侧外围拓展d×d/4，作为PML吸收层，以尽可能地降低弹性波在人工边界处虚假反射，提高数值模拟的结果准确性，且PML吸收层上、下边界为不渗透边界.渗流场设定为定压开采，左侧为注入端，右侧为采出端，初始时刻Pf=2.0 MPa/m.在计算区域左侧加载u=u0exp(iωt)的波场.

图1 计算区域与边界条件

2.2 波动耦合渗流与Darcy渗流对比分析

将波动耦合渗流模型模拟结果与Darcy渗流模型模拟结果对比，分析波动对流体渗流过程的影响，基础参数如表1所示.图2为不加载弹性波的Darcy渗流和加载弹性波耦合渗流两种情况下含水饱和度和含油饱和度的分布图.

表1 基础参数表

模拟发现当同样开采500 s时，相较与Darcy渗流水驱前缘(饱和度在x=4 m处发生明显变化)，波动耦合渗流中水驱前缘(含油饱和度在x=8 m处发生明显变化)运动更快，离注入端更远.说明在一维定压开采过程中，弹性波作用于储层，给储层一个反复的振动作用力，改变了储层物性，促进了油水渗流，低频波动强化采油技术有较好的增渗增油效果，从理论上说明波动可以提高水驱波及范围和水驱效率的机理.

(a)Darcy渗流含油饱和度变化

(b)波动耦合渗流含油饱和度变化图2 作用500 s时波动耦合渗流和 Darcy渗流油水分布变化图

2.3 波动参数对渗流的影响

前期研究表明振动频率为固有频率时激励效果最佳，本模型模拟基于参数表，令ω=15 Hz，η0=100 mPa·s，k0=50×10-3μm2，φ0=20%，模拟不同振幅u0=100～350μm对含油分布、油水渗流速度和孔隙压力的影响，为振动采油参数优化提供依据.

2.3.1 振幅对含油饱和度影响规律分析

模拟结果如图3所示，相同模拟时间内(500 s)，当振幅u0在100～350μm之间变化时，随着振幅的增大，水驱前缘向左偏移，即振幅越小时同一位置处油相优先被水相驱替.当u0<300μm时，水驱前缘在空间上随振幅的变化不明显，一方面可能是因为驱替时间较短，另一方面可能是振幅足够强时饱和度变化敏感性变弱，即在一定储层物性特征条件下可选择合适的振幅强度作为工作振幅.

(a)u0=100 μm

(b)u0=150 μm

(c)u0=200 μm

(d)u0=250 μm

(e)u0=300 μm

(f)u0=350 μm图3 500 s时不同振幅下含油 饱和度的空间分布

对不同振幅条件下轴线A-A′含水饱和度变化进行绘图，如图4所示.由于计算时间较小，含水饱和度变化差距较小；弹性波场激励下饱和度沿轴线A-A′发生波动性变化，呈先急剧减小、后波动性增加与减小的趋势；根据急剧减小段(x<10 m)的变化速率可知，振幅愈大，弹性波场激励作用愈强，对含水饱和度的影响愈大，这一定程度上亦反映了实验揭示的波动降低束缚水饱和度的现象.

图4 500 s时不同振幅下A-A′线含水 饱和度变化

2.3.2 振幅对流速影响规律分析

模拟结果如图5所示，在图5(a)中水相流速随振幅增大波动变化，但变化范围不大(0～0.5×10-5m/s之间)，当振幅u0=150μm时，水相流速最小.在图5(b)中油相速度随振幅的增加先增大后减小，当振幅u0=250μm时油相速度与水相速度最接近.油水相速度越接近，流速越小，水驱突进程度越低，波及面积越大，驱油效果越明显.油藏内部相对注入端，振动对油藏内部油水相流速的促进明显降低，说明随传播距离的增加，波动引起的渗流速度变化逐渐降低，反映了波衰减情况下对渗流的影响减弱.综合水相和油相的速度变化规律，最佳振幅为150～250μm之间.

(a)水相速度随振幅的变化规律

(b)油相速度随振幅的变化规律图5 500 s时振幅对A-A′线处液体速度 的影响规律

随着距离的增加，从注入水端面到油水交界处(变化的曲面)到纯油水区，油水饱和度逐渐发生变化，毛管力、油水相对渗透率亦发生变化，水相流速出现先增加后减小的情况；论文分析过程中考虑了波场对渗流场的激励作用，尤其是低渗孔隙介质中激励作用愈强，流固耦合交错下流速变化幅度愈大，因此波动耦合渗流场流速变化比达西流速变化更为剧烈.另外由于弹性波传播过程是一个周期性的变化过程，随距离和时间均会出现相位变化(振幅出现正负之分)，波诱导流速既可实现正向激励，又可表现为反向抑制，当然在考虑毛管力作用下实际孔隙中油水流度以正向激励为主，因此在沿距离上流速出现了波动性的变化趋势，表现为流速速度振动剧烈；且可能在一定位置或时间时反向波诱导流速大于原渗流场流速情形，使得图5(b)中流速出现负值，当孔隙介质渗透率较高、波诱导流速较小，原渗流场流速较高时，油水两相流速即可全部表现为正值.

2.3.3 振幅对压力影响规律分析

图6为弹性波作用下，500 s时A-A′上的压力随振幅的变化规律.在距振动端较近的位置(x<12 m)，压力随振幅的增大先减小后增大，在u0=250μm时达到最大，随后减小，反映了在距波源较近的位置处u0=250μm为最佳振幅.在油藏中部(12 m

综合以上各物性特征的变化规律，建议低频波动强化采油最佳振幅为150～300μm，矿场应用时根据设备和具体井筒、储层情况选择.

图6 500 s时A-A′线处孔隙 压力随振幅的变化规律

3 结论

(1)基于Biot理论和等效流体理论，考虑油水两相初始Darcy渗流、低频波动下孔隙-应变关系、状态方程、连续性方程和辅助方程，引入实际开发油藏边界条件和初始条件，建立低频波作用下一维油水两相流固耦合渗流模型.

(2)通过对模型求解和敏感性参数分析，认为与常规油水两相Darcy渗流规律相比，波动耦合渗流场中储层受到一个反复的振动作用力，改变了储层物性，促进了油水渗流，说明低频波动强化采油技术有较好的增渗增油效果.在地层固有频率15 Hz下，基于振幅对含油饱和度、水相流速和油相流速、压力等的影响规律，优选出低频波动强化采油最佳振幅为150～300μm.