引言

非线性渗流是中国低渗透油气田开发的重要研究内容之一，启动压力梯度更是其中的研究热点。低渗透井在试油试采阶段的瞬态压力响应的特点是：关井所测的压力及其导数在后期呈平行上翘趋势。众多学者对此进行了研究，提出启动压力梯度的概念，并从实验、机理、应用等方面对其进行了研究。多数学者认为：低速非达西渗流是由边界层、细微的喉道、非牛顿流体等因素共同形成的。其表现形式为启动压力梯度。当压力梯度大于启动压力梯度时，并非所有喉道内的流体都参与流动，而是随着驱动压力梯度的增加，不断有更多小喉道的流体参与流动，显现为非线性。这常称为拟启动压力梯度。

众多学者在实验中发现了启动压力梯度，并研究了其流动规律。闫庆来等^[1]通过室内实验发现流体在低渗透储层内渗流时，存在非线性段和启动压力梯度。黄延章、葛家理提出了利用带拟启动压力梯度的分段函数描述低渗透油藏非线性渗流^[2-3]。邓英尔等提出了带有3个参数的单一函数非线性渗流模型^[4]。李爱芬等^[5]用实验研究了水在特低渗岩芯中的启动压力梯度。同时，给出了几组由实验获得的启动压力梯度，如当渗透率为0.5560 mD时，最小启动压力梯度为0.04543 MPa/m；当渗透率为4.1229 mD时，最小启动压力梯度为0.00609 MPa/m。阮敏等^[6]基于实验数据用蒙特卡洛数值解法给出了临界压力梯度与流体黏度、岩石渗透率间的关系图。

同时，启动压力梯度的机理研究也受到了重视。姚约东等^[7]研究了线性与非线性的判别标准，并用边界层理论进行分析。徐绍良等^[8]以离子水在半径为2.5 μm的微圆管中流动进行实验，认为边界层是致密低渗油气藏非线性流的主要原因之一。时宇等^[9]利用毛管模型及边界层理论，建立了含启动压力梯度的非线性渗流模型，并给出了几组由实验获得的启动压力梯度，如渗透率为1.010 0 mD时，启动压力梯度为0.042 85 MPa/m；渗透率为3.970 0 mD时，启动压力梯度为0.015 50 MPa/m。

近年来，陈元千^[10]认为，启动压力梯度不能用于平面径向流方程。王晓冬等^[11]认为启动压力梯度能够简明扼要地表述低速非Darcy渗流，是低渗透储层有效开发的理论基础；并认为启动压力梯度决定于流体的性质、介质的表面作用和孔隙结构，并不受渗流方式的影响。

在压力导数特征方面，人们普遍认为，启动压力梯度使压力及其导数曲线后期呈上翘趋势, 而且启动压力梯度越大, 压力及其导数曲线抬升越高，上升得也越早。通过叠加原理得到的关井段的井底压力导数曲线与生产段的井底压力导数曲线特征相同，即在关井段，压力与压力导数也后翘，且意味着关井后的压力可以恢复到初始压力(若油藏边界很大)。蔡明金等^[12]研究了低渗透双重介质的试井模型，并进行了实例解释。王晓冬等^[13]研究了考虑启动压力梯度的非稳态压力分析，并运用Duhamel褶积建立了相应的现代试井分析典型曲线图版，并基于现场数据进行了试井解释，算例中的启动压力梯度为0.157 00 MPa/m。姚军等^[14]则用数值方法研究了启动压力梯度的低渗透油藏的试井解释模型。另外，尹芝林等^[15]通过油水两相数值模拟研究了启动压力梯度对压力场的影响。王道富等^[16]则研究了启动压力梯度在超前注水方面的应用。

然而，低渗储层的非达西渗流理论自20世纪80年代提出至今，仍受到部分学者的质疑。有人认为低渗非达西渗流是由实验误差和错误产生的。例如，通常认为黏土中的渗流是低渗非达西渗流。但是王秀艳等^[17]通过实验发现，黏性土中孔隙水运动不存在绝对意义上的初始水力梯度和临界水力梯度。李传亮等^[18]认为启动压力梯度属于实验假象，其实并不存在，低渗透储层与中高渗透储层没有本质区别；认为启动压力梯度主要是由于流量测试困难和测试时间太短所致，表皮效应加剧了启动压力存在的假象，注水井存在启动压力是由于井底污染和润湿滞后所致，流体在地层中流动并不需要启动压力梯度^[19]。窦宏恩^[20]对此提出异议，指出低渗(特低渗)储层渗流存在启动压力梯度，而且启动压力梯度对低渗透油藏注水井网优化设计、压裂设计及注水压力的确定都有很高的指导意义和价值。窦宏恩同时指出，不是所有低渗透储层流体渗流都存在“启动压力梯度”，只有当储层压力系数低于1时，低渗透储层才会表现出一种特殊的渗流特征。同样属于低渗透致密储层，如果成藏原始压力较高，低渗透致密储层流体渗流就根本不需要“启动压力梯度。

由于含启动压力梯度的渗流方程没有解析解，现有的含启动压力梯度的渗流方程都是利用数值方式进行求解的。然而，在进行瞬态压力特征分析中，很多文献都是基于叠加原理对关井段的井底压力进行求解。含启动压力梯度的渗流方程已不是齐次方程，已不满足叠加原理的适用条件。为此，本文基于PEBI网格，直接利用数值方式求解含启动压力梯度的单相渗流方程，通过数值模拟研究启动压力梯度的瞬态压力响应特征，并与叠加原理的解进行对比，发现叠加原理不适合求解含启动压力梯度的渗流方程。

1 基本方程

对低渗透油藏作如下假设：油井以定产量生产；流体弱可压缩，且压缩系数为常数；考虑启动压力梯度；考虑井筒储集和表皮效应；在本文的算例中，储层是水平的，为简洁起见，忽略下面公式中的重力项。含启动压力梯度的Darcy方程有多种方式，但最常用的形式如下

$ \left\{ {\begin{array}{*{20}{l}} {\boldsymbol{u} =-\frac{K}{\mu }\left( {\nabla p-\lambda } \right)\ \ \ \nabla p > \lambda }\\ {\boldsymbol{u} = 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \nabla p < \lambda } \end{array}} \right. $

(1)

式中：

u-速度，m/s；

K-绝对渗透率，mD；

p-压力，Pa；

λ-压力梯度，Pa/m；

μ-黏度，Pa·s。

方程(1)中的启动压力梯度不是拟启动压力梯度，本文称之为绝对启动压力梯度。由于启动压力梯度总是阻止流动的，一般地，方程(1)可写为

$ \left\{ {\begin{array}{*{20}{l}} {\boldsymbol{u} =-\frac{K}{\mu }\left( {\nabla p-f\lambda } \right)\left| {\nabla p} \right| > \lambda }\\ {\boldsymbol{u} = 0\left| {\nabla p} \right| < \lambda } \end{array}} \right. $

式中：

$ \left\{ {\begin{array}{*{20}{l}} {f = 1{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \nabla p > \lambda }\\ {f =-1{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}-\nabla p > \lambda } \end{array}} \right. $

为简洁起见，在后面的叙述中，仍用方程(1)的形式。将方程(1)代入质量守恒方程，用地面密度去除等是两边，并引入体积系数B，有

$ \nabla \cdot \left[{\frac{K}{{\mu B}}\left( {\nabla p-\lambda } \right)} \right] = \frac{\partial }{{\partial t}}\left( {\frac{\phi }{B}} \right) + {q_1} $

(2)

式中：

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

q₁-质量流量与地面密度的比值，1/s；

B-体积系数，无因次。

2 含启动压力梯度方程的数值求解 2.1 PEBI网格

矩形网格很不灵活，主要表现在：不能描述复杂油藏的边界形状；对区域较大、井数较多的油藏，油井不会都位于网格中心；对水平井或斜井，网格很难与井的方向保持一致；存在严重的网格取向效应。非结构网格PEBI(Perpendicular bisection)是一种局部正交网格，任意两个相邻网格块的交界面一定垂直平分相应网格节点的连线，如图 1所示。Heinemann Z E等首次将其应用到油藏模拟中^[21]。研究表明，PEBI网格具有如下优点^[22-23]：比结构网格灵活；近井处可以局部加密并且粗细网格过渡较为平滑；PEBI网格取向效应小；易于构造断层等。因而，PEBI网格很快被引入到数值试井领域，并在数值试井领域率先得到了商业化应用。

2.2 数值求解格式

数值计算中，对方程(2)进行积分，考虑井筒污染(S)后有

$ {q_{{\rm{sc}}}} = {J_{\rm{w}}}(p-{p_{{\rm{wf}}}}) $

(3)

式中：

${q_{{\rm{sc}}}} $-体积流量，m³/s；

${J_{\rm{w}}} = \frac{{\theta Kh}}{{\mu B\ln (r/{r_{\rm{w}}}) + S}} $；

$p_{{\rm{wf}}} $-井底压力，Pa；

$\theta $-井的射开角度，(°)；

h-地层厚度，m；

r-井与相邻网格的距离，m；

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

S-表皮系数，无因次。

假定井口产量为Q，储层为单层，井为内边界，则可得耦合井储C的井控制方程

$ Q = \sum\limits_k {{q_{{\rm{sc}}, k}}}-\frac{C}{{\Delta t}}\left( {p_{{\rm{wf}}}^{n + 1}-p_{{\rm{wf}}}^n} \right) $

(4)

式中：

Q-体积流量，m³/s；

C-井储，m³/Pa；

下标k-与井相邻的网格；

$q_{{\rm{sc}}, k} $-与井相邻的网格k与井间的流量，m³/s；

$\sum\limits_k {{q_{{\rm{sc}}, k}}} $-整个井的流量，m³/s。

将方程(4)对网格i进行体分数(图 1)，经过若干推导，其隐式的离散格式为

$ \sum\limits_j {T_{ij}^{n + 1}} {\rm{ }}\left[{p_j^{n + 1}-p_i^{n + 1}-{d_{ij}}({\lambda _j}-{\lambda _i})} \right] =\\ \frac{{{V_i}}}{{\Delta t}}\left[{\left( {\frac{\phi }{B}} \right)_i^{n + 1}-\left( {\frac{\phi }{B}} \right)_i^n} \right] + q_{{\rm{sc}}}^{n + 1} $

(5)

式中：

$p_i $，$p_j $-网格i与网格j的压力，Pa；

$d_{ij} $-网格中心距离，m；

$\lambda_i $，$\lambda_j $-网格i、网格j的启动压力梯度，Pa/m；

$V_i $-网格i的体积，m³；

$T_{ij} $-传导系数，在径向流动区域和井的T_ij要进行修正，相关内容可参见文献[22]。

流动项的传导率采用空间上游加权，采用全隐式线性化方法，写成迭代的方式，对流动项，有

$ \sum\limits_j {{{\left\{ {{T_{ij}}[\Delta p-{d_{ij}}({\lambda _j}-{\lambda _i})]} \right\}}^{n + 1}}} = \\ \sum\limits_j {T_{ij + }^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}\left[{p_j^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}-p_i^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}-{d_{ij}}({\lambda _j}-{\lambda _i})} \right]} + \sum\limits_j {T_{ij + }^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}\left( {\delta p_j^{\begin{array}{*{20}{c}} {{\rm{ }}}\\ {{\rm{ }}} \end{array}} - \delta p_i^{\begin{array}{*{20}{c}} {{\rm{ }}}\\ {{\rm{ }}} \end{array}}} \right)} + \\ \sum\limits_j {\left( {\frac{{\partial {T_{ij}}}}{{\partial p}}} \right)_ + ^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}\left[{p_j^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}-p_i^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}-{d_{ij}}({\lambda _j}-{\lambda _i})} \right]\delta p_ + ^{\begin{array}{*{20}{c}} {{\rm{ }}}\\ {{\rm{ }}} \end{array}}} $

(6)

其中，下标“＋”表示上游加权。

累积项采用守恒方式展开，写成迭代形式，有

$ \frac{{{V_i}}}{{\Delta t}}\left[{\left( {\frac{\phi }{B}} \right)_i^{n + 1}-\left( {\frac{\phi }{B}} \right)_i^n} \right] = \\ C_{{\rm{op}}}^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}\delta {p_i} + C_{{\rm{op}}}^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}\left( {p_i^v -p_i^n} \right) $

(7)

式中：

$ \delta {p_i} = p_i^{\begin{array}{*{20}{c}} {(v + 1){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}-p_i^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}, $

$ C_{{\rm{op}}}^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}} = \frac{{{V_i}}}{{\Delta t}}\left[{\frac{{{\phi _{{\rm{ref}}}}{C_{\rm{r}}}}}{{B_{\rm{o}}^n}} + \frac{{\phi _i^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}C_{{\rm{o}}, {\rm{ref}}}^{}}}{{B_{{\rm{o}}, {\rm{ref}}}^{{\rm{ref}}}}}} \right] $

井既可以使用Peaceman模型处理成源汇，也可把其当成内边界。这里将井作为对内边界处理。对方程(3)进行全隐式线性化，写成迭代形式，有

$ q_{{\rm{sc}}}^{n + 1} = q_{{\rm{sc}}}^v + \frac{{\partial J_{\rm{w}}^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}}}{{\partial {p_i}}}\left( {p_i^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}-p_{{\rm{wf}}}^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}} \right)\delta p_i^{\begin{array}{*{20}{c}} {{\rm{ }}}\\ {{\rm{ }}} \end{array}} +\\ J_{\rm{w}}^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1} \end{array}}\delta p_i^{\begin{array}{*{20}{c}} {{\rm{ }}}\\ {{\rm{ }}} \end{array}}-J_{\rm{w}}^{\begin{array}{*{20}{c}} {(v){\rm{ }}}\\ {n + 1{\rm{ }}} \end{array}}\delta p_{{\rm{wf}}}^{} $

(8)

对方程(4)进行线性化，有

$ Q = \sum\limits_k {q_{{\rm{sc}}, k}^v}-\frac{C}{{\Delta t}}\left( {p_{{\rm{wf}}}^{\begin{array}{*{20}{c}} {(v)}\\ {n + 1} \end{array}}-p_{{\rm{wf}}}^n} \right) +\\ \sum\limits_k {\frac{{\partial J_{{\rm{w}}, k}^{\begin{array}{*{20}{c}} {(v)}\\ {n + 1} \end{array}}}}{{\partial {p_i}}}\left( {p_i^{\begin{array}{*{20}{c}} {(v)}\\ {n + 1} \end{array}}-p_{{\rm{wf}}}^{\begin{array}{*{20}{c}} {(v)}\\ {n + 1} \end{array}}} \right)\delta p_i^{\begin{array}{*{20}{c}} {}\\ {} \end{array}}} +\\ \sum\limits_k {J_{{\rm{w}}, k}^{\begin{array}{*{20}{c}} {(v)}\\ {n + 1} \end{array}}\delta p_i^{\begin{array}{*{20}{c}} {}\\ {} \end{array}}} - \left( {\sum\limits_k {J_{{\rm{w}}, k}^{\begin{array}{*{20}{c}} {(v)}\\ {n + 1} \end{array}}} + \frac{C}{{\Delta t}}} \right)\delta p_{{\rm{wf}}}^{\begin{array}{*{20}{c}} {}\\ {} \end{array}} $

(9)

将方程(6)、方程(7)和方程(8)代入方程(5)，就可得线性化的方程组，再加上井控制方程(9)、边界条件与初始条件，就可对网格的压力及井底流压进行耦合求解。

对单相渗流问题，虽然可采用更为简单的线性化方法、井变量也可采用解耦的方法进行求解，但我们在黑油模型的数值求解中统一使用的是全隐式耦合求解的方法。

需要提及的是，由于启动压力梯度没有方向性，在程序中需要判断网格间的压力梯度的正负号，并作相应调整。相关细节这里不详细讨论。

3 基于叠加原理的含启动压力梯度方程求解 3.1 无因次方程

对于均质情形，含启动压力梯度的平面径向流渗流方程可表示为

$ \frac{{{\partial ^2}p}}{{\partial {r^2}}} + \left( {\frac{1}{r} + \lambda {{\bar c}_{\rm{f}}}} \right)\frac{{\partial p}}{{\partial r}}-\frac{\lambda }{r} = \frac{1}{\chi }\frac{{\partial p}}{{\partial t}} $

(10)

式中：$\chi = \frac{K}{{\phi \mu {{\bar c}_{\rm{t}}}}} $；

$ {\bar c_{\rm{f}}} = \frac{{{c_{\rm{f}}}}}{{1 + {c_{\rm{f}}}\Delta \bar p}}; $

$ {\bar c_{\rm{t}}} = \frac{{{c_{\rm{t}}}}}{{1 + {c_{\rm{f}}}\Delta \bar p}}; $

$c_{\rm{f}} $-流体的压缩系数，MPa^-1；

$c_\phi $-岩石压缩系数，MPa^-1；

$c_{\rm{t}}=c_{\rm{f}}+c_\phi $-综合压缩系数，MPa^-1；

$\bar p $-压力恢复期间的地层平均压力，MPa。

引进无因次量

$ {r_{\rm{D}}} = \frac{r}{{{r_{\rm{w}}}}}, $

$ {t_{\rm{D}}} = \frac{{Kt}}{{\phi \mu c{r_{\rm{w}}}}}, $

$ {\lambda _{\rm{D}}} = \frac{{2{\rm{\pi }}Kh{r_{\rm{w}}}\lambda }}{{Q\mu }}, $

$ {c_{{\rm{fD}}}} = \frac{{Q\mu }}{{2{\rm{\pi }}Kh}}{c_{\rm{f}}}, $

$ {p_{\rm{D}}} = \frac{{2{\rm{\pi }}Kh({p_{\rm{i}}}-p)}}{{BQ\mu }}, $

$ {r_{{\rm{De}}}} = {r_{\rm{D}}}{{\rm{e}}^S}, $

$ {P_{{\rm{De}}}} = {P_{\rm{D}}}, $

$ {R_{{\rm{De}}}} = {R_{\rm{D}}}{{\rm{e}}^S}$，且分别用${t_{\rm{D}}}{{\rm{e}}^{2S}} $和${C_{\rm{D}}}{{\rm{e}}^{2S}} $代替$t_{\rm{D}} $和$C_{\rm{D}} $，式\(10\)变为

$ \frac{{{\partial ^2}{p_{{\rm{De}}}}}}{{\partial {r_{{\rm{De}}}}^2}} + \frac{1}{{{r_{{\rm{De}}}}}}\frac{{\partial {p_{{\rm{De}}}}}}{{\partial {r_{{\rm{De}}}}}} + \frac{{{\lambda _{\rm{D}}}{{\rm{e}}^{-S}}}}{{{r_{{\rm{De}}}}}} = \\ \frac{1}{{{C_{\rm{D}}}{{\rm{e}}^{2S}}}}\frac{{\partial {p_{{\rm{De}}}}}}{{\partial ({t_{\rm{D}}}/{C_{\rm{D}}})}} $

(11)

式中：

$p_{\rm{i}} $-原始地层压力，MPa。

内边界条件为

$ {\left. {\frac{{{\rm{d}}{p_{{\rm{De}}}}}}{{{\rm{d}}({t_{\rm{D}}}/{C_{\rm{D}}})}}-\frac{{\partial {p_{{\rm{De}}}}}}{{\partial {r_{{\rm{De}}}}}}} \right|_{{r_{{\rm{De}}}} = 1}} = 1 + {\lambda _{\rm{D}}}{{\rm{e}}^{-S}} $

(12)

运动外边界条件

$ {\left. {\frac{{\partial {p_{{\rm{De}}}}}}{{\partial {r_{{\rm{De}}}}}}} \right|_{{r_{{\rm{De}}}} = {r_{{\rm{cD}}}}}} =-{\lambda _{\rm{D}}}{{\rm{e}}^{-S}} $

(13)

封闭外边界

$ {\left. {\frac{{\partial {p_{{\rm{De}}}}}}{{\partial {r_{{\rm{De}}}}}}} \right|_{{r_{{\rm{De}}}} = {R_{{\rm{De}}}}}} = 0 $

(14)

初始条件

$ {\left. {{p_{{\rm{De}}}}({r_{{\rm{De}}}}, {t_{\rm{D}}}/{C_{\rm{D}}})} \right|_{{t_{\rm{D}}}/{C_{\rm{D}}} = 0}} = 0 $

(15)

式(12)~(16)就是最终的数学模型。上述方程的数值求解格式较为简单，这里不详细讨论。

3.2 叠加原理

在渗流力学中，叠加原理常用来求解关井段的井底流压随时间的变化情况，其示意图如图 2所示。这意味着，关井段的试井曲线同生产段的试井曲线是类似的。

实际生产史为(图 2a)：以产量q进行生产，生产t_p时间后关井，关井时间为$\Delta t $。在上一节的数值计算中，是以产量q进行生产$t_{\rm p} $时间后的压力分布为关井期间的油藏压力分布的初始值，然后再进行关井段的流动模拟计算。

叠加原理的处理方式为(图 2b)：以产量q进行生产，生产的时间为$t_{\rm p}+\Delta t $；在$t_{\rm p} $时间后以-q的产量进行注入，注入时间为$\Delta t $。因而，使用叠加原理进行求解时，在计算以-q的产量注入时，认为地层的压力是处处相等的，注入时间为$\Delta t $。这种处理方式会使压力漏斗恢复。

叠加原理要求方程是齐次的，而含启动压力梯度的渗流方程(12)不是齐次的。因而，叠加原理的适用性条件不满足。若仍使用叠加原理进行求解，就会引起误差。然而，该误差没能引起众多研究者的注意。其中很重要的原因是，在生产段，启动压力梯度使油藏流动的面积逐渐变大，出现“动边界效应”，从而压力导数曲线上翘。根据叠加原理计算出来的压力导数曲线上翘，而很多低渗透关井压力恢复导数曲线也上翘。计算曲线与实测曲线特征相符，这让人们误认为，虽然方程不是齐次的，但叠加原理引起的误差小，其结果仍可用。

4 计算与讨论

若无特别说明，油藏为1 000 m×1 000 m的矩形封闭地层(图 3)，直井位于油藏中心，产量、生产时间、关井时间、网格大小见具体算例。

其他参数取值如下：

初始压力=20 MPa；

油藏厚度=10 m；

水平渗透率=5.922 mD；

孔隙度=0.2；

岩石压缩系数=0.000 15 MPa^-1；

黏度=0.002 mPa·s；

体积系数=1.05；

井储=0.2 m³/MPa；

表皮因子=5；

井半径=0.1 m。

4.1 数值解与基于叠加原理解的对比

图 4给出了流量为5 m³/d的直井生产120 d，然后关井30 d，启动压力梯度为0.05 MPa/m的数值解与基于叠加原理解的井底流压对比。在生产段，数值解与基于无量纲方程的解吻合的很好。但在关井段，二者的井底压力差别很大。基于叠加原理的井底压力不仅升的高，升的快，而且还没有稳定。数值解的井底压力很快稳定，压力不再升高。这矛盾的结果是不同的计算方法造成的。

叠加原理认为，关井后的井底压力可用一个产量为q的生产井和一个产量为-q注入井的二者井底压力之和来获得。因而，关井后，井保持以q产量继续生产；同时虚拟一个-q注入量，在油藏压力为初始压力的前提下进行注入。虚拟的注入井的井底压力会很快升高，二者之和当然也会升高。

在数值解中，计算方式与叠加原理完全不同。在生产段的最后时刻，压降漏斗已经形成，井附近的压力梯度大，井远的地方的压力梯度小。关井后，以该压力分布为初值进行计算。因而关井后，离井远的区域的压力梯度很快降低，并很快低于启动压力梯度。随着关井时间的增加，更多的区域会停止流动，压力恢复的区域变小，压力恢复困难。而且，由于启动压力梯度所形成的压降漏斗会一直存在，没有外力干扰的情形下不会消失。这就是数值解的井底压力升高慢，很快稳定的原因。

由上分析可知，由于含启动压力梯度的渗流方程不是齐次的，叠加原理的使用条件不满足，若强行使用叠加原理进行求解会导致完全错误的结果。

图 5给出了数值解与基于叠加原理的试井曲线对比图。

基于叠加原理的试井曲线后期上翘(图 5a)，而数值解的试井曲线后期下掉(图 5b)。二者的趋势完全不同。试井曲线后期下掉(图 5b)的特征容易解释。由于启动压力梯度的存在，小于启动压力梯度的区域没有流动，因而整个生产期间都是一个封闭油藏，只是油藏边界在不断扩大。试井曲线也应表现出封闭油藏曲线特征，即试井曲线后期下掉。并且，随着时间的增加，在边界附近的区域的压力梯度会很快下降，并且很快小于启动压力梯度，停止流动，所以，油藏可流动的区域很快变小，试井曲线应下掉更快。

图 6给出了启动压力梯度为0.05 MPa/m，关井4 800 d时的井底压力曲线。图 6表明，即使关井达十多年，井底压力仍在上升，仍没稳定。这显然难以解释。启动压力梯度应使流动更加困难，关井后稳定更快，但基于叠加原理的解却表明，流动难以稳定。显然这是由于叠加原理的不适用性造成的。

4.2 讨论

首先，本文的研究结果不能回答“启动压力梯度是否存在”的问题。本文所采用的启动压力梯度是绝对启动压力梯度，即压力梯度小于某个压力梯度，流动就不存在。事实上，众多研究者认为，实际的低渗油藏中不存在绝对的启动压力梯度，存在一个非线性段。在实际的岩芯中，孔喉直径存在一个分布。小的孔喉对应大的启动压力梯度，大的孔喉对应小的启动压力梯度。因而，当油藏开始生产时，大孔喉的孔道启动压力梯度小，先开始流动；随着压力梯度的增加，具有中等尺寸孔喉的孔道开始流动，直至最后小的喉道里的流体开始流动。针对这种情形，很多学者都提出了拟启动压力梯度的概念。理论上拟启动压力梯度会使关井压力恢复试井曲线就会上翘。这部分内容将另文论述。因而本文的结果仅能证明含启动压力梯度的渗流方程不能使用叠加原理进行求解。

其次，本文的研究还表明，实际的低渗油田中不存在绝对的启动压力梯度。在实际的试油试采阶段，关井末点压力值往往很高，有时仅略低于初始压力。显然，绝对启动压力梯度不能解释这种现象。当存在绝对启动压力梯度时，压降漏斗一旦形成，在没有外力、边水、注入等干扰的情形下，井周围的压力漏洞当永远存在。

5 结论

(1) 叠加原理不能用于求解含启动压力梯度的渗流方程。含启动压力梯度的渗透方程不再是齐次方程，叠加原理的适用条件不成立。由叠加原理导致的误差，会使结果失真。

(2) 绝对启动压力梯度使得生产段的井底压力导数曲线上翘，关井段的井底压力导数曲线下掉；启动压力梯度降低关井恢复的末点压力值。这与实际低渗油藏的曲线特征不同。因而实际的低渗油藏不存在绝对的启动压力梯度。