0 引言

低渗透致密储层喉道半径差异较大、孔隙结构复杂、非均质性强，注水开发效果较差^[1-3]。由于该类储层孔隙尺度为纳米—微米级别，连通性较差，开发过程中存在大量剩余油^[4-5]，CO₂驱在低渗透油藏可大幅驱替剩余油，具有较好的驱油效果^[6]。CO₂驱分为混相驱和非混相驱。当注入压力大于最小混相压力时为混相驱，反之为非混相驱。当非混相驱压力增大时，采收率曲线会出现向上升高的转折点，将此时的注入压力与最小混相压力之间的区域称为近混相驱压力范围^[7]，本文将近混相驱归为非混相驱范畴。非混相驱的驱替过程中CO₂与原油之间因存在界面张力而不能形成混相区域，且CO₂指进窜流现象明显^[8]，所以研究CO₂在孔隙中的渗流特征对进一步认识CO₂非混相驱的驱油机理和气窜机制具有重要意义。

目前对于非混相驱渗流特征的研究主要通过室内实验^[9-11]和数值模拟^[12-15]进行。Xiao等^[16]通过CO₂驱替实验和核磁共振测量技术，研究了气水交替驱、非混相驱、混相驱和循环注CO₂等在特低渗透油藏中的驱油效果，结果表明，注CO₂可改善驱油效果。Zhu等^[17]采用数值模拟方法在非均质圆形多孔介质中研究了CO₂在油湿多孔介质中的驱替过程。Ma等^[18]建立了非均质圆形多孔介质模型，分析了CO₂和原油在非均质多孔介质中的渗流现象和CO₂驱采收率。以往研究多基于理想非均质圆形多孔介质模型，无法表征流体在实际孔喉中的渗流特征。

针对上述问题，利用水平集法数值模拟来表征CO₂与原油的界面变化，建立CO₂非混相驱数值模型，模拟扩散吸附作用下CO₂非混相驱和近混相驱的微观渗流规律，对CO₂在孔隙中的微观渗流特征和扩散吸附特征展开研究，并选取注入速度、扩散系数、吸附反应速率常数等对近混相驱微观渗流特征的影响进行分析。

1 CO₂非混相驱数值模拟 1.1 数值模型建立

模拟储层多孔介质内CO₂非混相驱的微观渗流特征，需做出以下假设：①多孔介质内流体为不可压缩流体，流体之间及流体与孔隙壁之间不发生化学反应，多孔介质在驱替过程中不发生形变；②CO₂在原油中的扩散系数在驱替过程中不产生变化，扩散现象符合菲克定律；③CO₂和原油被视为牛顿流体；④非混相驱过程中CO₂与原油之间的界面张力为常数；⑤非混相驱过程中多孔介质润湿性及两相接触角不发生变化。

基于假设，不可压缩流体的质量平衡用如下连续性方程表示：

$ \nabla \cdot \vec{u}=0 $

(1)

式中：▽为哈密顿算子；$\vec{u}$为流体矢量速度，m/s。

对于不可压缩流体，采用Navier-Stokes方程对流体的动量平衡^[17]进行控制：

$ \begin{gathered} \rho \frac{\partial \vec{u}}{\partial t}+\rho(\vec{u} \cdot \nabla) \vec{u}= \\ \nabla \cdot\left\{-p I+\mu\left[\nabla \vec{u}+(\nabla \vec{u})^T\right]\right\}+F_{\mathrm{st}} \end{gathered} $

(2)

式中：ρ为流体（CO₂或油）密度，kg/m³；t为时间，s；p为压力，Pa；I为单位矩阵；μ为流体（CO₂或油）黏度，Pa·s；F_st为用于CO₂或油界面的表面张力，N/m³。

驱替过程中，采用水平集法数值模拟确定CO₂和原油的界面。两相流体界面由水平集函数φ的等值线值0.5来表示。原油中φ = 0，CO₂中φ = 1。在靠近界面的过渡层中，φ从0平滑过渡至1。分隔2个相态的流体界面的传递由下式^[19]表示：

$ \frac{\partial \varphi}{\partial t}+\vec{u} \cdot \nabla \varphi=\gamma \nabla \cdot\left[\varepsilon \nabla \varphi-\varphi (1-\varphi) \frac{\nabla \varphi}{|\nabla \varphi|}\right] $

(3)

式中：φ为水平集函数；γ为重新初始化参数，m/s；ε为界面厚度控制参数，m。

两相界面上的密度突变和黏度突变可进行平滑处理：

$ \rho=\rho_1+\left (\rho_2-\rho_1\right) \varphi $

(4)

$ \mu=\mu_1+\left (\mu_2-\mu_1\right) \varphi $

(5)

式中：ρ₁为原油密度，kg/m³；ρ₂为CO₂密度，kg/m³；μ₁为原油黏度，Pa·s；μ₂为CO₂黏度，Pa·s。

结合水平集函数，表面张力由以下公式进行计算：

$ F_{\mathrm{st}}=\sigma \delta \kappa \vec{n} $

(6)

$ \delta=6|\varphi (1-\varphi)||\nabla \varphi| $

(7)

$ \kappa=-\nabla \cdot \vec{n} $

(8)

$ \vec{n}=\frac{\nabla \varphi}{|\nabla \varphi|} $

(9)

式中：σ为表面张力系数，N/m；δ为狄拉克δ函数，仅在流体界面为非零值；κ为曲率；$\vec{n}$为界面单位法向。

CO₂在多孔介质中向原油的扩散可通过下式表示：

$ \frac{\partial c}{\partial t}+\nabla \cdot (-D \nabla c)+\vec{u} \cdot \nabla c=0 $

(10)

式中：c为CO₂浓度，mol/m³；D为扩散系数，m²/s。

CO₂在多孔介质壁面的吸附方程为

$ \frac{\partial c_{\mathrm{s}}}{\partial t}+\nabla \cdot\left (-D_{\mathrm{s}} \nabla c_{\mathrm{s}}\right)=k_{\mathrm{ads}} c\left (\Gamma_{\mathrm{s}}-c_{\mathrm{s}}\right)-k_{\mathrm{des}} c_{\mathrm{s}} $

(11)

式中：c_s为流体表面浓度，mol/m²；D_s为表面扩散速率，m²/s；k_ads为吸附速率常数，m³/（mol·s）；k_des为解吸速率常数，1/s；Γ_s为吸附区总表面浓度，mol/m²。

CO₂表面吸附浓度的求解还需与多孔介质内壁的通量方程进行耦合：

$ -\vec{n} \cdot (-D \nabla c)=k_{\mathrm{des}} c_{\mathrm{s}}-k_{\mathrm{ads}} c\left (\Gamma_{\mathrm{s}}-c_{\mathrm{s}}\right) $

(12)

1.2 数值模型验证

基于非均质多孔介质模型，将水平集法CO₂驱数值模型与Ma等^[18]等采用相场法得到的CO₂驱数值模型的采收率曲线进行对比。结果显示2种模拟方法的采收率曲线随时间变化的趋势较为接近，注入3 s时，相场法CO₂驱数值模拟采收率为51.29%，水平集法CO₂驱数值模拟采收率为53.60%。水平集法CO₂驱数值模拟采用自由边界等值面曲率计算表面张力，且以Navier滑移处理基质边界处的流体速度；相场法CO₂驱数值模拟利用化学势计算表面张力，基质边界无滑移（图 1）。由此可知，水平集法数值模拟更适用于CO₂非混相驱的渗流过程模拟。

1.3 二维孔隙结构建立

数值模拟所采用的实际岩心孔隙度为12%，渗透率为1.2 mD，根据图 2a所示的实际岩心铸体薄片的电镜扫描结果，利用AutoCAD软件处理后可获得实际岩心的微观二维孔隙结构（图 2b），制成可供CO₂驱微观渗流模拟的模型。

1.4 材料参数

水平集法数值模拟采用的原油与CO₂在343 K温度下的最小混相压力为12.4 MPa^[18]。原油与CO₂在不同注入压力下的参数如表 1所列^[20]。

1.5 计算参数

水平集法数值模型耦合了层流物理场、水平集物理场和稀物质传递物理场。在层流物理场中注入速度为0.04 m/s，出口压力为0；在水平集物理场中，水平集变量指定相为原油，入口处注入相为CO₂，气相接触角为60°；在稀物质传递物理场中初始CO₂浓度为0，入口浓度在非混相驱和近混相驱时分别为3 950.3 mol/m³和7 269.1 mol/m³，CO₂在壁面吸附时的表面扩散速率为1×10^-11 m²/s。在求解数值模型时，选用瞬态探究，步长为0.000 1 s，计算时长为0.04 s。

1.6 网格划分参数

将图 2（b）的二维孔隙结构基于物理场控制网格方法对孔隙域进行网格划分。孔隙结构尺寸为350 μm×350 μm，网格形状为三角形，三角形数92 029个，网格顶点数52 410个，边单元数13 051个，顶点单元数2 934个，域单元数92 029个，最小单元质量为0.345 4，平均单元质量为0.822 2，单元面积比为7.768×10^-4，网格面积为43 940 μm²。二维孔隙结构网格划分如图 3所示。

2 模拟结果分析 2.1 不同驱替条件下的微观渗流特征

非混相驱时，CO₂注入后与原油之间存在清晰界面，CO₂前缘在模型中部的运移速度最快，形成指进现象，CO₂到达采出端后，模型上、下部原油被波及驱替至采出端，驱替末期，靠近采出端的模型边部及小孔隙中存在滞留剩余油（图 4a）。近混相驱的毛细数为9.68×10^-4，是非混相驱的7.29倍，在黏性力增大和界面张力减小的作用下，两相接触面积和油气过渡区域均增大；当CO₂注入0.012 s时，模型中下部的波及面积大于非混相驱，注入压力的增大使CO₂的密度和黏度增大，改善了气油流度比，有利于驱替前缘进入小孔隙；驱替末期，与非混相驱相比，部分孔隙中滞留的剩余油被驱出，剩余油面积减小（图 4b）。

驱替初期，非混相驱和近混相驱的采收率增长趋势较为一致（图 5），随驱替时间增加，非混相驱和近混相驱采收率曲线的增长趋势均不断减缓，但近混相驱采收率数值高于非混相驱。在驱替至模拟结束时（0.04 s），非混相驱采收率为87.7%，近混相驱采收率为91.1%。

2.2 不同驱替条件下CO₂浓度

当驱替温度为343 K时，CO₂在驱替压力为8.0 MPa时的入口浓度为3 950.3 mol/m³，扩散系数为8.85×10^-9 m²/s；在驱替压力为11.5 MPa时的入口浓度为7 269.1 mol/m³，扩散系数为1.19×10^-8 m²/s。驱替过程中孔隙内CO₂的浓度分布如图 6所示。非混相驱条件下，CO₂优先向模型中部孔隙扩散，注入前期浓度分布呈锥形，CO₂前缘向原油中扩散，使外围边部的浓度低于孔隙中部，随驱替时间增加，CO₂逐渐向未波及的小孔隙扩散，在驱替末期，仅靠近采出端的模型边部有少量孔隙未被波及（图 6a）。近混相驱条件下，孔隙内波及区域的CO₂浓度通常大于非混相驱，CO₂更易向小孔隙扩散，有利于扩大与原油的接触面积（图 6b）。

由图 7可知，驱替压力为8.0 MPa和11.5 MPa的见气时间分别为0.002 7 s和0.002 5 s。驱替中后期曲线逐渐平稳后，近混相驱的气体体积分数高于非混相驱；驱替结束时，近混相驱的出口气体体积分数为97.01%，大于非混相驱的71.60%，说明近混相驱的驱油速度较非混相驱快。

2.3 不同驱替条件下CO₂吸附浓度

由图 8（a）可知，孔隙内壁CO₂吸附浓度随注入时间增加而增大，吸附浓度的增加方向为先中间后两边，与孔隙内CO₂浓度变化趋势一致。非混相驱替结束时，孔隙内壁CO₂最大表面吸附浓度为3.16×10^-4 mol/m²。由图 8（b）可知，同时刻下近混相驱孔隙内壁CO₂最大吸附浓度大于非混相驱，近混相驱替结束时，孔隙内壁CO₂最大表面吸附浓度为5.81×10^-4 mol/m²。

取孔隙内壁上某一个点S分析CO₂吸附浓度随时间的变化（图 9）。非混相驱和近混相驱条件下，S点见气后吸附浓度曲线均呈上升趋势，且近混相驱的吸附浓度曲线上升幅度大于非混相驱；在驱替结束时（0.04 s），非混相驱的表面吸附浓度为2.66×10^-4 mol/m²，近混相驱的表面吸附浓度为5.21×10^-4 mol/m²。

3 近混相驱微观渗流特征影响因素 3.1 注入速度

选取注入速度为0.01 m/s和0.07 m/s进行驱替模拟。注入速度的增大可使同时刻下CO₂的波及面积增大，从而加快了孔隙内CO₂的扩散。不同注入速度下的采收率曲线均呈现为先快速上升，之后上升幅度不断减小，最后趋于平缓的趋势（图 10）。注入速度增大可使驱替力增大，从而使更多CO₂进入小孔隙，使剩余油面积不断减少，提高了采收率。驱替至模拟结束时（0.04 s），注入速度为0.01 m/s的采收率为68.66%，注入速度为0.07 m/s的采收率为97.02%。

由图 11可见，注入速度增大可使见气时间提前，含气率曲线达到平稳的时间也提前，并使驱替结束时出口含气率较高。注入速度0.01 m/s和0.07 m/s的出口见气时间分别为0.010 1 s和0.001 5 s。驱替结束时，注入速度0.01 m/s和0.07 m/s的出口气体体积分数分别为84.4% 和98.4%。CO₂注入速度的增大可使前缘移动速度加快，使出口见气时间提前，含气率趋于稳定后形成气窜通道，提高了出口含气率。

3.2 扩散系数

随注入时间增加，扩散系数增大，CO₂易于向复杂小孔隙中的原油扩散。由图 12可知，不同扩散系数的采收率曲线随注入时间增加均呈先快速上升后趋于平缓趋势。扩散系数为1.19×10^-9 m²/s和1.19×10^-7 m²/s时的采收率分别为91.65% 和92.27%。扩散系数的增大有利于增大波及面积，因扩散系数的增大对CO₂运移速度影响较小，所以对波及范围的影响仅限于在原油中的扩散溶解作用。

不同扩散系数的含气率曲线变化趋势接近（图 13），注入后期曲线接近平稳时，扩散系数较大时出口含气率较低。不同扩散系数的采出端见气时间相同，均为0.002 5 s，驱替结束时扩散系数1.19×10^-9 m²/s和1.19×10^-7 m²/s的出口气体体积分数分别为96.67% 和94.93%。扩散系数的增大会增加CO₂与原油的接触面积，使出口的产油量略高，所以出口含气率值会稍有降低。

3.3 吸附反应速率常数

吸附反应速率常数的增大有利于中后期模型边部区域及壁面剩余油的驱替，从而增大波及面积。由图 14可见，CO₂注入0.02 s之前2条曲线的变化趋势相似，均为先快速上升，之后上升幅度逐渐减小。CO₂注入0.03 s后，吸附反应速率常数数值较大时的采收率曲线比数值较小时上升幅度大。驱替至模拟结束时（0.04 s），吸附反应速率常数2×10^-7 m³/（mol·s）的采收率为90.95%，2×10^-5 m³/（mol·s）的采收率为92.08%。CO₂吸附反应速率常数的增加有利于驱替后期剩余油从壁面剥离。

由图 15可见，不同吸附反应速率常数对应的含气率曲线变化趋势接近，在出口见气后均急速上升，之后上升幅度出现波动，最后趋于平稳。不同吸附反应速率常数的采出端见气时间相同，均为0.002 5 s。驱替结束时，吸附反应速率常数2×10^-7 m³/（mol·s）和2×10^-5 m³/（mol·s）的出口气体体积分数分别为97.47% 和97.14%。吸附反应速率常数的变化对CO₂运移速度和出口含气率的影响均较小。驱替结束时，吸附反应速率常数的增大会使孔隙内CO₂浓度减小，减少了气体分子数量，从而使出口含气率略有降低。

4 结论

（1）非混相驱条件下，CO₂优先向模型中部孔隙扩散，注入前期CO₂浓度分布呈锥形；近混相驱条件下，CO₂更易向小孔隙扩散，有利于扩大与原油的接触面积。驱替结束时，近混相驱的最终采收率和出口气体体积分数分别较非混相驱高3.4% 和25.41%，驱油速度较快。

（2）近混相驱条件下，注入速度的增大使同时刻下CO₂的波及面积增大，使孔隙内CO₂扩散加快，从而使见气时间提前，驱替结束时出口气体体积分数较高。驱替结束时，注入速度的增大使采收率增加了28.36%，出口气体体积分数增加了14%。

（3）近混相驱条件下，随注入时间增加，扩散系数增大时，CO₂易于向复杂小孔隙中的原油扩散；吸附反应速率常数的增大会降低CO₂在孔隙内的扩散浓度。驱替结束时，扩散系数的增大使采收率增加了0.62%，出口气体体积分数降低了1.74%；吸附反应速率常数的增大使采收率增加了1.13%，出口气体体积分数降低了0.33%。