近年来, CO₂驱油技术成为低渗透储层提高采收率的热点技术，低渗储层CO₂驱油相对渗透率的研究也被纳入重点研究内容。鄂尔多斯盆地陕北地区延长组的致密砂岩油藏作为中国低渗透油藏的典型代表，学者对其开展了大量CO₂驱油相关技术的研究和矿场实践探索^[1-3]。然而遗憾的是, 目前尚无合适的相渗模型能真实反映该区CO₂驱油的渗流特征。已实施CO₂驱油试验的油区，虽进行过相渗测试, 但其数据处理仍沿用常规岩心的处理方法。以期用传统的“JBN”法获得的“拟相对渗透率曲线”代替真实低渗透储层CO₂驱相渗曲线, 不具有针对性和适用性, 不能真实反映该区储层CO₂驱相渗特征, 对科研及生产产生了误导。在此背景下，本研究立足于低渗透砂岩储层的非线性渗流规律, 全面考虑CO₂驱油过程中溶解降黏、气相滑脱等影响因素, 建立低渗透砂岩油藏CO₂非混相驱相渗计算模型, 以期更加真实地揭示陕北地区致密砂岩储层CO₂驱相渗特征，以提高油藏工程及数值模拟等研究的可靠性。

1 非线性渗流规律与数学描述

CO₂注入低渗透储层后, 由于低渗透储层的复杂孔喉特征、流体的非牛顿流变学特征以及流体/骨架间的相互作用等内在因素, 导致流体在低渗透砂岩储层中产生非线性渗流^[4-8]。

不同学者先后通过不同的数学模型对低速非线性渗流进行了描述, 且呈现了形式由简单到复杂、参数由少到多、精度由低到高的发展历程(见表 1)。拟启动压力模型认为，低于拟启动压力梯度时流体不流动, 同时将高于拟启动压力梯度的低速非线性流简化为通过拟启动压力梯度表征的线性流, 形式过于简单, 得出的分析结果与生产实际相差甚远。直线逼近、二次函数逼近、幂函数逼近等分段模型, 通过不同的函数形式拟合非线性段, 得到了较启动压力模型更接近实际情况的分析结果, 但这种分段模型需要确定线性与非线性渗流的临界点, 增加了应用的难度^[9]。邓英尔等^[10]、杨清立等^[11]、姜瑞忠等^[12]提出的连续函数模型, 通过一个连续方程分析低渗储层介质非线性渗流规律, 解决了分段模型需要确定临界点的问题。但是，这些模型中包含的参数并不能直接来源于实验, 需要通过庞大的计算过程得出, 且无明确的物理意义。黄延章等^[13]提出的三参数连续模型，将最小启动压力梯度(λ_a)、拟启动压力梯度(λ_b)、线性段斜率($\frac{K}{\mu }$)等参数纳入非线性渗流方程，3个参数由实验直接获得, 且物理含义明确, 能很好地反映非线性渗流特征。

2 低渗透砂岩油藏CO₂非混相驱相对渗透率计算模型

基于黄延章等^[13]提出的三参数非线性渗流方程, 考虑CO₂非混相驱过程中的气相滑脱、CO₂在原油中的溶解及溶解后流体黏度的变化等, 综合Buckley-Leverett^[14-15]驱油机理, 建立低渗透砂岩油藏CO₂非混相驱非稳态测算相渗曲线的计算模型。

2.1 不混溶气驱相渗模型的建立 2.1.1 运动方程

根据三参数非线性渗流的规律, 为简化计算, 油、气相的非线性运动方程改写为:

$ \begin{array}{l} - \frac{{\partial {p_{\rm{o}}}}}{{\partial x}} = \frac{{{Q_{\rm{o}}}{\mu _{\rm{o}}}}}{{AK{K_{r{\rm{o}}}}\left[ {1 - \frac{{{\lambda _{b{\rm{o}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{o}}}} - {\lambda _{a{\rm{o}}}}} \right)}}} \right]}} = \\ \frac{{u{f_{\rm{o}}}{\mu _{\rm{o}}}}}{{K{K_{r{\rm{o}}}}\left[ {1 - \frac{{{\lambda _{b{\rm{o}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{o}}}} - {\lambda _{a{\rm{o}}}}} \right)}}} \right]}}, \end{array} $

(1)

$ \begin{array}{l} - \frac{{\partial {p_{\rm{g}}}}}{{\partial x}} = \frac{{{Q_{\rm{g}}}{\mu _{\rm{g}}}}}{{AK{K_{r{\rm{g}}}}\left[ {1 - \frac{{{\lambda _{b{\rm{g}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{g}}}} - {\lambda _{a{\rm{g}}}}} \right)}}} \right]}} = \\ \frac{{u{f_{\rm{g}}}{\mu _{\rm{g}}}}}{{K{K_{r{\rm{g}}}}\left[ {1 - \frac{{{\lambda _{b{\rm{g}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{g}}}} - {\lambda _{a{\rm{g}}}}} \right)}}} \right]}}。\end{array} $

(2)

式中:P为压力, MPa; x为沿驱替方向的流动距离, m; μ为流体黏度, mPa·s; Q为流量, cm³/s; u为流速, cm/s; A为岩心截面积, cm²; f为流体分流量, 无因次; |▽P|为岩心两端的压力梯度, MPa/m; K为绝对渗透率, ×10^-3μm²; K_r为相对渗透率, 无因次; λ_a，λ_b分别为启动压力梯度和拟启动压力梯度, MPa/m; 下标o，g分别代表油相、气相。

2.1.2 等饱和度面移动方程

根据Buckley-Leverett驱油机理, CO₂驱油岩样沿程等饱和度面移动方程为

$ x = \frac{{{{f'}_{\rm{g}}}\left( {{S_{\rm{g}}}} \right)}}{{\varphi A}}\int_0^t {Q\left( t \right){\rm{d}}t} , $

(3)

出口端L处为

$ L = \frac{{{{f'}_{\rm{g}}}\left( {{S_{{\rm{g}}e}}} \right)}}{{\varphi A}}\int_0^t {Q\left( t \right){\rm{d}}t} , $

(4)

$ \frac{x}{L} = \frac{{{{f'}_{\rm{g}}}\left( {{S_{\rm{g}}}} \right)}}{{{{f'}_{\rm{g}}}\left( {{S_{{\rm{g}}e}}} \right)}}。$

(5)

引入无因次注气量V(t)，并代入L表达式:

$ \bar V\left( t \right) = \frac{{\int_0^t {Q\left( t \right){\rm{d}}t} }}{{\varphi AL}} = \frac{1}{{{{f'}_{\rm{g}}}\left( {{S_{{\rm{g}}e}}} \right)}}。$

(6)

式中: S_g为含气饱和度; S_ge为岩心出口端含气饱和度; f_g^′ (S_g)为气相分流量对含气饱和度的导数; f_g^′ (S_ge)为气相在岩心出口端的分流量对含气饱和度的导数; φ为岩样孔隙度。

2.1.3 油相相对渗透率

根据运动方程得岩心两端压差Δp为：

$ \begin{array}{l} \Delta p = {p_1} - {p_2} = - \int_0^L {\frac{{\partial p}}{{\partial x}}{\rm{d}}x} = \\ \frac{{u{\mu _o}}}{K}\int_0^L {\frac{{{f_o}}}{{{K_{r{\rm{o}}}}\left[ {1 - \frac{{{\lambda _b}}}{{\left( {\left| {\nabla P} \right| + {\lambda _b} - {\lambda _a}} \right)}}} \right]}}{\rm{d}}x} 。\end{array} $

(7)

根据等饱和度面移动方程有:

$ {\rm{d}}x = \frac{L}{{{{f'}_{\rm{g}}}\left( {{S_{{\rm{g}}e}}} \right)}}{\rm{d}}{{f'}_{\rm{g}}}\left( {{S_{\rm{g}}}} \right). $

(8)

将dx，V(t)带入Δp中得:

$ \begin{array}{l} \Delta p = \frac{{u{\mu _o}}}{K}\frac{L}{{{{f'}_{\rm{g}}}\left( {{S_{{\rm{g}}e}}} \right)}},\\ \int_0^{{{f'}_{\rm{g}}}\left( {{S_{\rm{g}}}} \right)} {\frac{{{f_{\rm{o}}}}}{{{K_{r{\rm{o}}}}\left[ {1 - \frac{{{\lambda _b}}}{{\left( {\left| {\nabla P} \right| + {\lambda _b} - {\lambda _a}} \right)}}} \right]}}{\rm{d}}{{f'}_{\rm{g}}}\left( {{S_{\rm{g}}}} \right)} , \end{array} $

(9)

整理得:

$ \begin{array}{l} \int_0^{{{f'}_{\rm{g}}}\left( {{S_{\rm{g}}}} \right)} {\frac{{{f_{\rm{o}}}}}{{{K_{r{\rm{o}}}}\left[ {1 - \frac{{{\lambda _b}}}{{\left( {\left| {\nabla P} \right| + {\lambda _b} - {\lambda _a}} \right)}}} \right]}}{\rm{d}}{{f'}_{\rm{g}}}\left( {{S_{\rm{g}}}} \right)} = \\ \;\;\;\;\;\frac{{\Delta pK}}{{u{\mu _{\rm{o}}}L\bar V\left( t \right)}} = \frac{1}{{\bar V\left( t \right)I}} \end{array} $

(10)

其中I为注入能力比, $I = \frac{{u{\mu _{\rm{o}}}L}}{{K\Delta p}}$。

公式(10)两边对f_g^′ (S_g)求导并整理得

$ \begin{array}{l} {K_{r{\rm{o}}}}\left( {{S_{{\rm{g}}e}}} \right) = \\ \frac{{{f_{\rm{o}}}\left( {{S_{{\rm{g}}e}}} \right)}}{{\left[ {1 - \frac{{{\lambda _b}}}{{\left( {\left| {\nabla P} \right| + {\lambda _b} - {\lambda _a}} \right)}}} \right]}}\frac{{{\rm{d}}\left[ {\frac{1}{{\bar V\left( t \right)}}} \right]}}{{{\rm{d}}\left[ {\frac{1}{{\bar V\left( t \right)I}}} \right]}}。\end{array} $

(11)

2.1.4 气相相对渗透率

由非线性渗流方程得岩心出口端气、油流量比为

$ \begin{array}{l} \frac{{{Q_{\rm{g}}}}}{{{Q_{\rm{o}}}}}\frac{{{K_{r{\rm{g}}}}/{\mu _{\rm{g}}}}}{{{K_{r{\rm{o}}}}/{\mu _{\rm{o}}}}}\frac{{1 - \frac{{{\lambda _{b{\rm{g}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{g}}}} - {\lambda _{a{\rm{g}}}}} \right)}}}}{{1 - \frac{{{\lambda _{b{\rm{o}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{o}}}} - {\lambda _{a{\rm{o}}}}} \right)}}}} = \\ \frac{{Q{f_{\rm{g}}}}}{{Q{f_{\rm{o}}}}} = \frac{{1 - {f_{\rm{o}}}}}{{{f_{\rm{o}}}}}。\end{array} $

(12)

因此, 气相相对渗透率为

$ \begin{array}{l} {K_{r{\rm{g}}}} = \frac{{1 - {f_{\rm{o}}}}}{{{f_{\rm{o}}}}}\frac{{{K_{r{\rm{o}}}}\left( {{S_{{\rm{g}}e}}} \right){\mu _{\rm{g}}}}}{{{\mu _{\rm{o}}}}} \times \\ \;\;\;\;\;\;\;\;\;\frac{{1 - \frac{{{\lambda _{b{\rm{o}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{o}}}} - {\lambda _{a{\rm{o}}}}} \right)}}}}{{1 - \frac{{{\lambda _{b{\rm{g}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{g}}}} - {\lambda _{a{\rm{g}}}}} \right)}}}}。\end{array} $

(13)

2.1.5 末端含油饱和度

根据Buckley-Leverett驱油机理, 岩样出口端含气、含油饱和度方程为

$ \begin{array}{l} {S_{{\rm{g}}e}} = \overline {{S_{\rm{g}}}} - \bar V\left( t \right) \cdot {f_{\rm{o}}}\left( {{S_{{\rm{g}}e}}} \right) = \\ \;\;\;\;\;\;\;\;{{\bar V}_{\rm{o}}}\left( t \right) - {f_{\rm{o}}}\left( {{S_{{\rm{g}}e}}} \right)\bar V\left( t \right),\\ {S_{\rm{o}}} = 1 - {S_{wc}} - {S_{{\rm{g}}e}}。\end{array} $

(14)

方程(11)，(13)，(14)构成了满足非线性渗流规律的不混溶气驱相对渗透率的计算模型。

2.2 考虑溶解及滑脱效应的相渗模型

上述相对渗透率计算模型是建立在注入气体与油相不混溶的基础之上的。实际中, 要真实准确地描述CO₂驱各相相对渗透率, 必须考虑CO₂在原油中的溶解、降黏及气相滑脱效应等因素的影响。

2.2.1 CO₂溶解作用

CO₂在原油中有极强的溶解性。根据陕北某致密砂岩油藏油样溶胀实验的结果, CO₂在地层条件下的溶解度达到60%(见图 1)。油气两相在地层中渗流时流体存在的状态有:油、束缚水、自由气、溶解气。而气相的状态取决于地层条件及其在油相中的溶解度大小。不同生产阶段的产出流体在地层压力P_f至计量器压力P_a条件下的状态如图 2所示。

1) 初期不产气阶段:CO₂尚未突破, 产出物仅为油。

2) 中期产少量气阶段:CO₂为分散相, 计量器中产气少而产油多, 在地层条件下CO₂以溶解气形式存在, 且地层压力下的产油量需经CO₂的溶解度α_o和原油的体积膨胀系数β_o进行矫正。

3) 后期产少量油阶段:CO₂突破后呈连续相分布, 地层条件下CO₂主要以自由气形式存在, 同时少量的气体溶于原油。此时产气量需通过地层平均压力进行矫正, 而地层压力下的产油量需经体积膨胀系数β_o进行矫正^[16]。

3个生产阶段的产出物在地层压力下的产量矫正为:

$ \begin{array}{l} {Q_{{\rm{g}}f}} = \left\{ \begin{array}{l} 0,{Q_{{\rm{g}}s}} = 0\\ 0,0 < {Q_{{\rm{g}}s}} < {\alpha _{\rm{o}}} \cdot {Q_{{\rm{o}}s}}\\ \frac{{2{p_a}}}{{\Delta p + 2{p_a}}}\left( {{Q_{{\rm{g}}s}} - {\alpha _{\rm{o}}} \cdot {Q_{{\rm{o}}s}}} \right),{Q_{{\rm{g}}s}} \ge {\alpha _{\rm{o}}} \cdot {Q_{{\rm{o}}s}} \end{array} \right.,\\ {Q_{{\rm{o}}f}} = \\ \left\{ \begin{array}{l} {Q_{{\rm{o}}s}},{Q_{{\rm{g}}s}} = 0\\ \frac{{{Q_{{\rm{g}}s}}}}{{{\alpha _{\rm{o}}}}} \cdot {\beta _{\rm{o}}} + \left( {{Q_{{\rm{o}}s}} - \frac{{{Q_{{\rm{g}}s}}}}{{{\alpha _{\rm{o}}}}}} \right),0 < {Q_{{\rm{g}}s}} < {\alpha _{\rm{o}}} \cdot {Q_{{\rm{o}}s}}\\ {Q_{{\rm{o}}s}} \cdot {\beta _{\rm{o}}},{Q_{{\rm{g}}s}} \ge {\alpha _{\rm{o}}} \cdot {Q_{{\rm{o}}s}} \end{array} \right.。\end{array} $

(15)

2.2.2 溶解降黏作用

CO₂溶于原油后可以有效降低原油黏度, 这是CO₂提高采收率的重要原理之一。图 3中反映CO₂对陕北地区某致密砂岩油藏原油的降黏幅度达到60%。因此, CO₂驱油相对渗透率测算过程中, 流体黏度的变化是必须要考虑的因素。

CO₂对原油的降黏作用是通过CO₂气体与原油组分间的相互作用实现的。组分模型计算原油/CO₂混合体系黏度需要的参数多, 且较为繁琐, 而经典的Arrhenius方程^[17]可以满足一般混合物黏度计算精度的需要。因此，计算过程中采用该方程计算流体黏度。其计算公式为

$ {\mu _{{\rm{mix}}}} = \mu _1^{{x_1}} \cdot \mu _2^{{x_2}}, $

(16)

式中:μ_mix为混合物黏度, mPa·s; x₁，x₂为不同相所占的比例; μ₁^x₁，μ₂^x₂为各相流体黏度, mPa·s。

2.2.3 气相滑脱作用

罗瑞兰等^[18]指出，当储层渗透率低于1×10^-3μm²时, 滑脱效应对气藏产量有明显影响, 且影响程度由渗透率和压力共同决定：渗透率越低, 气藏压力越低, 滑脱效应越显著。因此, 在低渗透储层气相渗流过程中, 滑脱效应是不可忽视的因素。

Klinkenberg^[19]最先提出气体滑脱效应, 并给出考虑滑脱因子和平均压力的气测渗透率表达式:

$ {K_{\rm{g}}} = K\left( {1 + \frac{b}{p}} \right)。$

(17)

气体滑脱因子b与分子平均自由程及毛管半径有关。为了更好地进行工程应用, Sampath等^[20]回归实验数据, 提出滑脱因子与孔隙度和渗透率存在以下关系式:

$ b = 0.0955{\left( {\frac{k}{\varphi }} \right)^{ - 0.53}}。$

(18)

应用该计算方法, 对上述气相相对渗透率计算模型进行修正得到:

$ \begin{array}{l} {K_{r{\rm{g}}}} = \frac{{1 - {f_{\rm{o}}}}}{{{f_{\rm{o}}}}}\frac{{{K_{r{\rm{o}}}}\left( {{S_{{\rm{g}}e}}} \right){\mu _{\rm{g}}}}}{{{\mu _{\rm{o}}}}} \times \\ \frac{{1 - \frac{{{\lambda _{b{\rm{o}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{o}}}} - {\lambda _{a{\rm{o}}}}} \right)}}}}{{1 - \frac{{{\lambda _{b{\rm{g}}}}}}{{\left( {\left| {\nabla P} \right| + {\lambda _{b{\rm{g}}}} - {\lambda _{a{\rm{g}}}}} \right)}}}}\frac{1}{{1 + \frac{b}{{\bar p}}}}。\end{array} $

(19)

式(11)，(14)，(15)，(16)，(19)即构成了低渗透砂岩油藏CO₂非混相驱过程中相对渗透率计算模型。

3 实例计算分析 3.1 非稳态CO₂非混相驱实验描述

实验样品主要包括岩心、地层水、原油及CO₂。岩样源自陕北地区致密砂岩油藏某油区的天然岩心。实验用地层水、原油均为按实际地层样品组成配置的模拟流体。CO₂纯度99.9%。样品详细参数见表 2~4。

实验装置主要包括动力泵、中间容器、岩心夹持器、压力监控器、流体计量装置等。实验流程图见图 4。

实验步骤据行业标准《SY/T5345-2007油水相对渗透率测定方法—非稳态法》执行, 具体包括：①依据地层流体组成对实验流体进行配置; ②岩样洗油、烘干, 几何尺寸及气测渗透率测量; ③抽真空, 饱和地层水并测算孔隙度; ④已饱和地层水的岩样置入驱替装置, 空气浴采用梯度升温的方式升至地层温度60℃, 并恒温24小时; ⑤定流速0.8 cm³/min油驱水, 饱和油, 建立束缚水饱和度, 并保证围压始终比驱替压力高出2 MPa, 待出口流量稳定, 测定束缚水饱和度下的油相有效渗透率; ⑥定流速0.8 cm³/min CO₂驱油, 并记录各个时刻的驱替压力、产油量及产气量; ⑦气驱至不产油, 测定残余油状态下的气相有效渗透率后结束实验。

3.2 相渗曲线计算及分析

利用建立的模型对上述实验数据进行处理, 得到陕北地区致密砂岩油藏某油区基于非线性渗流的, 考虑CO₂溶解、降黏、气相滑脱的CO₂非混相驱油相渗曲线(见图 5)。结果整体上表现出以下特点:①低含气饱和度下, 油相相对渗透率随含气饱和度的增加快速下降, 而气相相对渗透率则缓慢上升; 高含气饱和度下, 气相相对渗透率快速上升。②气相相渗端点值较低, 不足0.1, 表明气体渗流能力相对较弱。③束缚水饱和度、残余油饱和度等特征值均较高, 分别为35%，30%, 说明CO₂非混相驱替能力有限。④两相区范围窄, 约35%, 说明两相共渗能力和空间有限, 具有低黏、低渗油藏相对渗透率曲线的特征。

与传统的“JBN”方法计算的相渗结果相比, 本模型计算结果表现出明显的差异(见图 5):油相相对渗透率有较大幅度的提高; 气相渗透率有较大幅度的降低; 束缚气饱和度有所降低。本模型所反映出的这种差异在油藏工程分析中起到了举足轻重的作用, 因此渗流模式、CO₂在原油中的溶解降黏及气相滑脱效应是相渗测算中不可忽略的因素。

为分析渗流模式、CO₂在原油中的溶解降黏及气相滑脱效应等因素对致密砂岩储层CO₂非混相驱相渗的影响程度, 本研究进行单因素影响的相渗曲线计算, 结果如图 6~8所示。其中，图 7和图 8均基于非线性渗流模型进行计算。

相对于线性渗流模型, 基于三参数非线性渗流的致密砂岩油藏CO₂驱油相渗曲线表现出较大的不同：油相相对渗透率上升, 而气相相对渗透率下降, 等渗点明显右移。如果运用线性渗流模型进行分析和预测, 会得出相对于非线性渗流较差的驱油效果。

由图 7可以看出, 考虑CO₂气体在原油中的溶解降黏作用所计算的相渗曲线与忽略CO₂的溶解作用相比有较大不同。由于忽略CO₂气体的溶解, 气体均视为以自由气的形式产出, 致气相渗透率值偏高、油相渗透率值偏低。同时, 溶解了CO₂气体的部分原油黏度降低, 流动能力增加, 也可引起油相相对渗透率增加。

从图 8可以看出, 考虑CO₂气体的滑脱效应时, 油相渗透率基本不变, 气相相对渗透率有小幅度下降。CO₂气体在致密多孔介质中的流动包括滑脱流和渗流, 未经滑脱效应矫正时, 滑脱产生的气体流量同时被记入渗流流量, 导致气相渗透率偏大。

4 结论

1) 建立了低渗透砂岩油藏非稳态CO₂非混相驱相渗计算模型。该模型以低渗透储层的三参数非线性渗流规律为基础, 同时考虑了CO₂驱油过程中溶解、降黏、气相滑脱效应等因素的影响。与传统“JBN”方法计算的相渗结果相比, 本模型计算结果表现出了明显的差异。

2) 进行了非稳态CO₂驱油相渗实验, 得到了陕北地区某致密砂岩油藏CO₂非混相驱典型相渗曲线。计算并对比了渗流模式、CO₂的溶解降黏作用以及CO₂气体滑脱效应等单因素对相渗曲线的影响。结果表明, 非线性渗流及CO₂的溶解降黏均引起油相相对渗透率上升, 气相相对渗透率下降。考虑气相滑脱时, 气相相对渗透率有小幅度下降, 油相相对渗透率基本无变化。