引 言

由于低渗透多孔介质孔隙几何形状的复杂性、流体性质的特殊性和同时受到流动条件的制约，使得流体在低渗透多孔介质中渗流特征不再符合传统的达西渗流规律^[1-4]，表现为明显的非线性渗流特征，同时具有一定的启动压力梯度。很多学者基于实验和理论研究分析流体性质及流体与多孔介质内表面相互作用机制对低渗透油藏非线性渗流规律的影响。实验研究结果表明低渗透多孔介质固液界面间存在强烈的分子作用力^[5]。黄延章^[6]认为，吸附在多孔介质内表面的边界流体是导致低渗透油藏流体非线性特征和启动压力梯度存在的主要原因；低渗透油藏尤其是致密低渗透油藏的平均喉道半径与原油边界层在同一个数量级上，必须克服固液界面间的分子作用力才能使边界流体参与流动^[7]；李兆敏^[8]指出，低渗透多孔介质的流体在流动过程中不能保持恒定的特性，表现为Bingham型非牛顿流体特征，即施加的驱动压力梯度大于某个阀值压力梯度后流体才开始流动。在低渗透油藏非线性渗流模型表征方面通过拟合物理实验结果和理论推导建立相应的渗流模型，但其相应的参数在实验室难以获得，物理意义不明确，进而限制了其在油藏工程及数值模拟中的应用。

本文以流体边界层理论为基础，提出一种低渗透油藏流体非线性渗流新模型，在模型推导过程中获得动态启动压力梯度的连续变化函数，并以此模型为基础建立一种分析低渗透油藏动边界变化规律的数值分析方法，与线性渗流模型和拟线性渗流模型下地层压力的变化特征进行了对比；最后，研究了非线性系数对低渗透油藏动边界的扩展特征。

1 低渗透油藏流体非线性渗流模型

在低渗透多孔介质中，孔隙半径与流体边界层厚度处于同一个数量级上，流体与多孔介质内表面间的相互作用强烈，影响了流体通过低渗透多孔介质的渗流特征。毛细管实验结果表明^[9-10]，对于同一根毛细管，驱替压力梯度与流体边界层厚度成反比关系，压力越大，壁面黏滞层越薄，呈现指数衰减的规律。基于上述研究结果，假设流体边界层厚度与驱替压力梯度的变化关系为

$ \delta = r_0 \cdot \delta _0 \cdot {\rm{e}}^{ - c\left| {\nabla p} \right|} $

(1)

式中：

$\delta$-边界流体厚度，m；

$r_0 $-毛细管半径，m；

$\delta _0$-驱替压力梯度为零时的边界层厚度与毛细管半径之比，即静态边界层厚度，且$0 < \delta _0 < 1$；

c-流体边界层参数，m·Pa^-1；

$\nabla p$-驱替压力梯度，Pa·m^-1。

1.1 流体在毛细管中的速度及流量公式

低渗透油藏，孔隙及喉道半径细小，边界层效应使流体通过喉道时表现为非牛顿Bingham流体的流动特征。对于Bingham流体，当剪切应力$\tau$大于屈服应力$\tau_0$才开始流动的流体。$\tau_0$值通常是由实验测定的，它是Bingham流体的固有特性。因此，在多孔介质中只有在驱替压力梯度大干某个启动压力梯度的情形下才能实现渗流。Bingham流体应力本构方程^[11-12]，为

$ \tau = \tau _0 + \mu _{\rm{B}} \dot \gamma $

(2)

式中：

$\tau$-流体所受的切应力，Pa；

$\tau _0$-流体所受的屈服应力，Pa；

$\mu _{\rm{B}} $-Bingham流体的塑性系数，Pa·s；

$\dot \gamma$-流体的应变，s^-1。

毛细管模型被广泛用来推导多孔介质的渗流模型^[13-15]，甚至是非牛顿Bingham流体的渗流模型^[16-17]。假设一根半径为$r_0$的毛细管，方程(2)可以描述为

$ \dfrac{{{\rm{d}}\vec v}}{{{\rm{d}}r}} = - \dfrac{1}{{\mu _{\rm{B}} }}\left( {\dfrac{r}{2}\left| {\nabla p} \right| + \tau _0 } \right) $

(3)

式中：

$\vec v$-流体渗流速度，m/s；

r-流体通过毛细管的有效渗流半径，m，$r=r_0-\delta$。

假设边界流体的厚度为$\delta$，相当于有效流动半径减小$\delta$，对方程(3)从r到$r_0 - \delta $进行积分

$ \vec v = - \dfrac{1}{{\mu _{\rm{B}} }}\int_r^{r_0 - \delta } {\left( {\dfrac{r}{2}\left| {\nabla p} \right| + \tau _0 } \right){\rm{d}}r} $

(4)

则得到Bingham流体沿半径方向的速度分布

$ \vec v\!=\!- \dfrac{1}{{\mu _{\rm{B}} }}\left\{ {\left[{\dfrac{{(r_0-\delta )^2\!-\!r^2 }}{4}}\right]\left| {\nabla p} \right|\!+\!\tau _0 \left( {r_0\!-\!\delta\!-\!r} \right)} \right\} $

(5)

再对式(5)进行过流断面上的积分，得到单个毛管中的流量公式

$ q = \int_0^{r_0 - \delta } {\vec v \cdot 2\pi r{\rm{d}}r} = \dfrac{{\pi (r_0 - \delta )^4 }}{{8\mu _{\rm{B}} }}\left| {\nabla p} \right| - \\\hspace{4em}\dfrac{1}{{\mu _{\rm{B}} }}\dfrac{1}{3}\pi \tau _0 (r_0 - \delta )^3 $

(6)

式(6)中流体在驱动压力梯度$\nabla p$的作用下通过半径为$r_0$ 的毛细管，获得的流量由两个部分构成，第一部分为驱动压力作用得到的体积流速，第二部分表示流体在屈服应力$\tau _0$作用下引起的流量损失。在毛细管最小流量为零的前提下，得到Bingham流体屈服应力$\tau _0$对应启动压力梯度

$ G = \dfrac{{8\tau _0 }}{{3(r_0 - \delta )}} = \dfrac{{8\tau _0 }}{{3r}} $

(7)

式(7)表明：(1)流体通过低渗透多孔介质的启动压力梯度本质上是固-液界面间的屈服应力造成的，数值上与剪切应力大小成正比，与毛细管有效渗流半径成反比；(2)启动压力具有动态特征，随驱替压力梯度的增大，启动压力梯度越小。

特别的，当毛细管中静态边界层厚度为零时，得到毛细管中最小启动压力梯度

$ G_{\min } = \dfrac{{8\tau _0 }}{{3r_0 }} $

(8)

流体通过单个毛细管内的最小启动压力梯度可以应用静态法测量获得^[18]。

将式(1)、式(8)代入到式(7)中，得到任意驱替压力梯度下的动态启动压力梯度

$ G = \dfrac{{G_{\min } }}{{1- \delta _0 \cdot {\rm{e}}^{- c\left| {\nabla p} \right|} }} $

(9)

将式(7)代入到式(6)中，得到考虑流体边界层效应的毛管中的流量公式

$ q = \dfrac{{\pi (1- \delta _0 {\rm{e}}^{ - c\left| {\nabla p - G} \right|} )^4 r_0 ^4 }}{{8\mu _{\rm{B}} }}\left( {1 - \dfrac{G}{{\left| {\nabla p} \right|}}} \right)\nabla p $

(10)

式中：

G-流体在驱替压力梯度$\nabla p$下的启动压力梯度，Pa/m；

q-体积流速，m³/s。

1.2 通过低渗透多孔介质的非线性渗流模型

流体通过低渗透多孔介质的渗流模型，可借助毛管束模型推导。将油藏岩芯假设为一组半径为$r_i$的毛细管，总个数为N，则通过截面A的总流量

$ Q\!=\!\sum\limits_{i\!=\!1}^N {q\left({r_i }\right)}\!=\!\sum\limits_{i\!=\!1}^N {\dfrac{{\pi (1\!-\!\delta_{0i}{\rm{e}}^{\!-\!c_i \left|{\nabla p\!-\!G_i }\right|})^4 r_i ^4 }}{{8\mu_{\rm{B}}}}}\left({1\!-\!\dfrac{{G_i }}{{\left|{\nabla p}\right|}}}\right)\nabla p $

(11)

式中：

Q-流量，m³/s；

$r_i$-第i个毛细管半径，m；

$c_i$-第i个毛细管中的流体边界层参数，m/Pa；

$G_i$-流体通过第i个毛细管时需要克服的启动压力梯度，Pa/m；

$\delta _{0i}$-流体在第i个毛细管中的静态边界层厚度，无因次。

定义渗流速度

$ v = \dfrac{Q}{A} = \dfrac{{\sum\limits_{i = 1}^N {q\left( {r_i } \right)} }}{{\left( {\sum\limits_{i = 1}^N {\pi r_i ^2 } } \right)/\phi }} $

(12)

得到多孔介质非牛顿Bingham流体渗流模型为

$ \vec v=\left\{\begin{array}{l} -\dfrac{{K_0 }}{{\mu _{\rm{B}}}}\left[{1\!-\!\delta _{0\phi} {\rm{e}}^{-c_{\phi}\left( {\left| {\nabla p} \right|-G_{\phi}} \right)}}\right]^4 \left({1-\dfrac{G}{{\left|{\nabla p} \right|}}} \right)\nabla p, \left|{\nabla p} \right|>G_{\phi}\\[12pt] 0, \left|{\nabla p}\right|\leqslant G_{\phi} \\ \end{array} \right. $

(13)

式中：

A-多孔介质的横截面积，m²；

$\phi$-多孔介质孔隙度，无因次；

$c_{\phi}$-流体通过多孔介质的平均边界层参数，m/Pa；

$G_{\phi}$-流体通过多孔介质时需要克服的启动压力梯度，Pa/m；

$\delta _{0\phi}$-流体通过多孔介质的平均静态边界层厚度，无因次。

$K_0 = \dfrac{\phi }{8}\dfrac{{\sum\limits_{i = 1}^N {r_i ^4 } }}{{\sum\limits_{i = 1}^N {\pi r_i ^2 } }}$

$\phi = \dfrac{{\sum\limits_{i = 1}^N {\pi r_i ^2 } }}{A}$

$A = \left( {\sum\limits_{i = 1}^N {\pi r_i ^2 } } \right)/\phi$

$K_0$-多孔介质渗透率，D。

上式中，当$\delta_{0\phi}=0$时，式(13)即为传统的达西渗流模型；当$\delta_{0\phi}=0$，$G_{\phi}$取最小启动压力梯度时，式(13)即为考虑拟启动压力梯度作用的拟线性渗流模型。因此，传统的达西渗流模型和拟线性渗流模型是本文建立的非线性渗流模型的一种特殊形式。

1.3 渗流模型及参数物理意义

由上述推导过程可看出，低渗透多孔介质中Bingham流体渗流规律表现为非线性渗流特征，同时又存在动态启动压力梯度；边界层的动态变化是Bingham流体呈现非线性渗流特征的主要原因；低渗透油藏启动压力梯度是流体屈服应力形成的。

公式(13)中$G_{\phi}$和$\mu_{\rm{B}}$称为Bingham流体的流变学参数，物理意义明确^{[11-12, 16-17]}。$c_{\phi}$和$\delta_{0\phi}$为流体边界层参数。其中，$\delta_{0\phi}$为相对静态边界层厚度，属性由多孔介质的润湿性决定，表征在驱替压力梯度为零时边界层流体厚度与储层孔喉半径之比。对于单个毛细管，可通过文献[7, 9-10]中提供的实验方法获得。$c_{\phi}$为非线性渗流系数，表征多孔介质孔喉结构和润湿性影响下的非线性特征，物理意义为在毛细管中流体厚度随驱替压力梯度变化的程度。当流体性质给定的情况下，对于任意多孔介质，其数值的大小仅表明其非线性程度。对于工程应用上流体边界层参数$c_{\phi}$和$\delta_{0\phi}$可通过回归岩芯压力梯度与流量曲线获得。

1.4 模型适应性评价

用黏度2 mPa·s的模拟油在某岩芯中进行了压力梯度与流量的关系实验研究，岩芯及流体参数：空气渗透率20 mD，岩芯截面积0.001 m²。应用式(13)对实验结果进行了拟合，所得的结果如图 1所示。拟合得非线性系数$c_{\phi}$=9.537 m/MPa，相对静态边界层厚度$\delta_{0\phi}$=0.417，应用静态法测得最小启动压力梯度$G_{\rm{min}}$=0.029 8 MPa/m。实验与拟合结果对比可知，本文建立的非线性渗流模型可以有效描述低渗透油藏的流动特征。

图 2表明：随着驱替压力梯度的增大，低渗透多孔介质的启动压力梯度是一个连续变化的过程，当驱替压力梯度增大到一定时，启动压力梯度为常数，此时$G_{\rm{min}}$=0.029 8 MPa/m。

2 数学模型的建立

对低渗透油藏作如下假设：

(1) 油藏均质且等温；

(2) 流动属于一维水平流动，忽略重力及毛管力的影响；

(3) 流体为非牛顿Bingham流体，其流动方程满足式(13)；

(4) 岩石和流体微可压缩，满足以下状态方程

$ \phi (p) = \phi _0 \left[{1 + c_{\rm{r}} (p-p_0 )} \right] $

(14)

式中：

$\phi _0$-初始孔隙压力$p_0$对应的孔隙度，无因次；

$c_{\rm{r}}$-地层孔隙压缩系数，Pa^-1；

$p，p_0$-孔隙压力，初始孔隙压力，Pa。

地层中流体的密度是压力的函数，表示为

$ \rho (p) = \rho _0 {\rm{e}}^{ - \alpha _\rho \left( {p_0 - p} \right)} $

(15)

式中：

$\rho$-孔隙压力p对应的流体密度，kg/m³；

$\rho _0$-流体初始密度，kg/m³；

$\alpha _\rho$-密度修正系数，Pa^-1。

对于式(13)，并根据文献^[12]，得到Bingham流体的视黏度流变模型为

$ \left\{ \begin{array}{l} \mu = \mu _{\rm{B}} /\left[{1-\dfrac{G_\phi}{{\left| {\nabla p} \right|}}} \right];\left| {\nabla p} \right| > G_\phi \\ \mu = \infty ;\left| {\nabla p} \right| \leqslant G_\phi \\ \end{array} \right. $

(16)

根据上述控制方程(13)、式(14)~式(16)，一维低渗透多孔介质不稳定渗流的数学模型可表示为

$ \dfrac{\partial }{{\partial x}}\left( {\rho v} \right) + q = \dfrac{\partial }{{\partial t}}\left( {\rho \phi } \right) $

(17)

初始条件

$ p(x)|_{t = 0} = p_0 $

(18)

内边界条件

$ q = \dfrac{{2\pi r_{\rm{w}} Kh}}{\mu }\dfrac{{{\rm{d}}p}}{{{\rm{d}}r}}|_{r = r_{\rm{w}} } = p_0 $

(19)

式中：

t-时间，s；

$r_{\rm{w}}$-生产井半径，m；

K-地层渗透率，D；

h-地层厚度，m；

$\mu$-原油黏度，Pa·s；

q-油井地面产量，kg/s。

式(14)~式(19)组成了一维低渗透油藏非稳态渗流数学模型。

3 数值求解

上述建立的一维多孔介质流体流动模型中涉及到的流体密度、体积系数、原油黏度、启动压力梯度均为压力的函数，是高度非线性的，需要进行数值求解。在数值求解时需要对方程中的空间传导项、时间累积项和源/汇项进行差分离散。为保证离散后的代数方程组的稳定性，空间传导项采用二阶中心差分近似，时间累积项用一阶向后差分近似。流动方程(17)在节点i的离散格式可写为

$ \sum\limits_{j = 1}^M {F_{ij}^{n + 1} } + q_i^{n + 1} = \dfrac{{V_i }}{{\Delta t}}\left[{\left( {\rho \phi } \right)^{n + 1}-\left( {\rho \phi } \right)^n } \right] $

(20)

式中：

$q_i$-节点单元i的源/汇项的质量流量，kg/s；

$V_i$-节点单元i的体积，m³；

$\Delta t$-时间间隔，s；

$F_{ij}$-临近节点单元j向本地节点单元i的质量流量，kg/s；表示为

$ F_{ij} = T_{ij} \cdot \Delta \varPhi _{ij} $

(21)

式中：$T_{ij}$-相邻节点单元i和j之间的传导率，kg/(Pa·s)；与达西流动条件下传导率的不同在于考虑边界层动态变化特征，即

$ T_{ij} = \dfrac{{K_i \rho _i }}{{\mu _{{\rm{B, }}i} }}\left( {\dfrac{A}{d}} \right)_{ij} \left( {1 - \delta _{0\phi} {\rm{e}}^{ - c——{\phi}\left| {\dfrac{{ \Delta \varPhi }}{d}} \right|_{ij} } } \right)^4 $

(22)

式中：

$\Delta \varPhi _{ij}$-相邻节点单元i和j之间的势差，Pa；定义为

$ \Delta \varPhi _{ij} = (\varPhi _j - \varPhi _i )_{ij} = \Delta p_{ij} - G_{ij} \cdot d_{ij} $

(23)

式中：

$\Delta p_{ij}$-相邻节点单元之间的压力差，Pa；

$G_{ij}$-相邻节点单元之间的启动压力梯度，Pa/m；

$K_i$-节点单元i的渗透率，×10¹² D；

$\rho_i$-节点单元i的密度，kg/m³；

$d_{ij}$-相邻节点单元之间的中心距离，m。

为保证数值求解的稳定性和精度，应用全隐式方法线性化离散差分方程(20)，得到任意控制体网格单元i的残差方程为

$ r_i^{u, n + 1} {\rm{\!=\!}}\sum\limits_{j \ne i} {F_{ij}^{u, n + 1} }\!+\!q_i^{u, n + 1} \!-\!\dfrac{{V_i }}{{\Delta t}}\left[{\left( {\rho \phi } \right)_i^{u, n + 1}\!-\!\left( {\rho \phi } \right)_i^n } \right] $

(i)

应用Newton-Raphson迭代法求解方程(20)。对于一维单相流动，有M个质量守恒方程和1个井方程，有M个压力变量和1个井底流压变量，加上初边值条件，方程(20)封闭可解。对于上述方程，Newton-Raphson迭代格式为

$ \sum {\dfrac{{\partial r_i^{u, n + 1} }}{{\partial p_i }}\left( {\delta p_i^{u + 1, n + 1} } \right)}\!+\! \dfrac{{\partial r_i^{u, n + 1} }}{{\partial p_{{\rm{wf}}} }}\left( {\delta p_{{\rm{wf}}}^{u + 1, n + 1} } \right) \!=\!- r_i^{u, n + 1} $

(25)

每次迭代求解后，其压力变量更新形式为

$ p_i^{u + 1, n + 1} = p_i^{u, n + 1} + \delta p_i^{u + 1, n + 1} $

(26)

式中：

$\delta p_i$-网格单元i的相邻两次迭代的压力变化量，Pa；

i，j-节点单元编号；

M-节点单元总数；

n-当前时间步；

n+1-下一时间步；

u-迭代步；

$A_{ij}$-相邻节点i和j的单元共用面积，m²。

4 数值结果分析

应用一维低渗透油藏渗流模型研究非线性渗流参数对地层压力分布特征的影响。所使用的岩石及流体参数为L=100 m，h=10 m，$G_{\rm{min}}$=0.029 8 MPa/m，$K_0$=20 mD，$\mu$=2.0 mPa·s，q=5 m³/d，B=1.1 m³/m³，$p_{\rm{wf}}$=7 MPa，$p_0$=32 MPa，$\phi_0$=0.3，$c_\phi$=9.54 m/MPa，$\delta_{0\phi}$=0.417，$\lambda$=0.08 MPa/m，$\Delta d$=1 m。

在边界条件和井底流压相同的条件下采用考虑动态启动压力梯度的非线性渗流、考虑最小启动压力梯度的非线性和线性流动规律分别计算低渗透地层压力变化曲线，如图 3所示。在近井区域，采用线性渗流模型时地层压力下降缓慢；考虑动态启动压力梯度的非线性渗流模型时地层压力下降幅度最大；采用考虑最小启动压力梯度的非线性渗流模型时地层压力介于二者之间。在远井地带(流动边界附近)的一个相同泄油半径处，考虑非线性渗流的地层压力大于拟线性渗流模型和线性渗流情况。出现这种现象的主要原因是启动压力梯度的动态变化特征形成的，在近井区域随着驱替压力的增大，启动压力梯度逐渐减小至最小，而远井区域的驱替压力梯度逐渐减小，降低启动压力梯度的能力减弱，减小了地层压力的下降范围。

不同非线性系数下动边界变化曲线如图 4所示。从图 4中可以看出，非线性系数越小，动边界扩展到模型边界需要的时间就越长，任意时刻动边界的距离越小。

5 结 论

(1) 基于流体边界层理论，提出了一种可用于描述在多孔介质中流动特征非线性渗流模型，获得了动态启动压力梯度的连续变化函数。

(2) 应用考虑动态启动压力梯度的非线性渗流模型，建立了一维单相流体流动的数学模型，并应用全隐式方法对离散方程进行了线性化处理；应用数值方法对比研究了线性渗流模型、拟线性和本文非线性渗流模型下地层压力的分布特征和流体非线性渗流系数对动边界扩展的影响

(3) 数值结果分析表明：与线性渗流和拟线性流动规律相比，反映启动压力梯度的动态变化特征的非线性渗流规律下在近井区域地层压力下降幅度更大，地层压力下降范围小。非线性系数是反映多孔介质渗流阻力大小的综合参数，其数值越小，流体流动的非线性程度越强，动边界的扩展速度越慢。