引言

考虑启动压力梯度的多孔介质非达西渗流模型研究^{[1,2,3,4,5,6]}在低渗透油气藏开发、稠油油藏中非牛顿宾汉流体流动规律、聚合物提高原油采收率、黏性土壤固结以及水力压裂措施中支撑剂悬浮等方面具有广泛的实际应用背景.考虑启动压力梯度的多孔介质非达西渗流模型属于强非线性动边界问题^[7]，除非数学模型进行等价转化，否则动边界模型的精确解析解将较难求出；而渗流力学领域应用广泛的线性变换例如：拉普拉斯变换、傅里叶变换等不再适用.目前，该类动边界问题研究方法主要包括积分近似解析方法^[7,8,9,10]、基于格林函数的数值近似方法^[11,12]、无网格方法^[13]、(基于空间坐标变换的)有限差分方法^{[14,15,16,17,18,19]}和精确解析方法^[20,21](此处称作精确解析方法是为区别近似解析方法).

相似变量变换方法^[22]是精确解析求解非线性偏微分方程的有效方法，本文将通过广泛应用于传热学中求解经典斯蒂芬(Stefan)动边界问题精确解析解的相似变量变换法，对内边界变压力下、考虑启动压力梯度的一维低渗透多孔介质非达西渗流动边界模型进行解析求解(本文模型与文献^[7]所研究模型的内边界条件不同)；同时，采用基于空间坐标变换的有限差分数值方法^[19]对该动边界模型进行数值求解，并利用所求得的精确解析解对数值解进行了严格验证.

本文的研究内容对低渗透、特低渗透油气藏开发过程所涉及的低渗透多孔介质非达西渗流力学理论进行了进一步完善，为中国非常规油气藏(如：致密油、页岩油气、稠油油藏等)油田开发的试井解释与油藏数值模拟技术提供了一定的理论基础.

1 低渗透多孔介质非达西渗流动边界数学模型 1.1 数学模型建立

本文研究对象为考虑启动压力梯度的一维(半无限长)低渗透多孔介质的非达西渗流问题，定义如下假设条件^[20]：

(1) 流体流动为单相流动;

(2) 多孔介质为均质、各向同性且等温;

(3) 忽略重力的影响，设地层的压力梯度比较小;

(4) 流体和多孔介质是微可压缩的.

流体密度的状态方程为^[8]

$\rho = {\rho _i}\exp \left[{ - {C_f}\left( {{p_i} - p} \right)} \right]{\text{ }}$

(1)

式中，ρ为流体密度；ρ _i为流体初始密度； p为地层压力； p_i为原始地层压力； C_f为流体压缩系数.

岩石孔隙度的状态方程为^[8]

$\phi = {\phi _i}\exp \left[{ - {C_\phi }\left( {{p_i} - p} \right)} \right]{\text{ }}$

(2)

式中，Φ为孔隙度； Φ _i为初始孔隙度；C _Φ为孔隙压缩系数.

考虑启动压力梯度的低渗透多孔介质非达西渗流运动方程为^[1\hbox{-}2]

$\upsilon = \left\{ \begin{gathered} - \frac{k}{\mu } \cdot \frac{{\partial p}}{{\partial x}} \cdot \left( {1 - \frac{\lambda }{{\left| {\frac{{\partial p}}{{\partial x}}} \right|}}} \right),amp;\left| {\frac{{\partial p}}{{\partial x}}} \right| > \lambda \hfill \\ 0,0\left| {\frac{{\partial p}}{{\partial x}}} \right| \leqslant \lambda \hfill \\ \end{gathered} \right.$

(3)

式中， k为渗透率；μ为流体黏度；λ 为启动压力梯度； x为距离； v为渗流速度.

一维多孔介质渗流的连续性方程为^[20]

$ - \frac{{\partial \left( {\rho \upsilon } \right)}}{{\partial x}} = \frac{{\partial \left( {\rho \phi } \right)}}{{\partial t}}{\mkern 1mu} ,\;0 \leqslant x \leqslant s(t){\text{ }}$

(4)

式中， t为时间, s为动边界.

将状态方程式(1)~式(2)和运动方程式(3)代入连续性方程式(4)，推得数学模型的控制方程为^[20]

$\frac{{{\partial }^{2}}p}{\partial {{x}^{2}}}=\frac{\mu \phi _iC_t}{k}\frac{\partial p}{\partial t},0\le x\le s(t)\text{ }$

(5)

式中， C _t为综合压缩系数.

初始条件为

$\text{s (0) = 0 }$

(6)

${{\left. p \right|}_{t}}_{=0}={{p}_{i}}$

(7)

变压力下的内边界条件为

${{\left. p \right|}_{x=0}}=f(t)\text{ }$

(8)

式中， f代表随时间变化的井底压力函数.

动边界条件为^[7,15,19,20]

${{\left. p \right|}_{x=s(t)}}={{p}_{i}}$

(9)

${{\left. \frac{\partial p}{\partial x} \right|}_{x=s(t)}}=\lambda \text{ }$

(10)

动边界条件所表达的物理意义为^{[7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]}：渗流只发生在近井区域地层压力梯度大于启动压力梯度的范围内；而动边界以外的区域，地层压力梯度小于启动压力梯度，渗流不能发生，保持为原始地层压力.此为考虑启动压力梯度影响的非达西渗流模型与经典达西渗流模型之间的主要区别.

为转化为一般模型，定义如下无因次变量^[20]

$\begin{align} & {{x}_{D}}=\frac{x}{{{x}_{w}}},{{t}_{D}}=\frac{k}{\mu {{\phi }_{i}}{{C}_{t}}x_{w}^{2}}t,\delta =\frac{s}{{{x}_{w}}},\\ & {{P}_{D}}=\frac{k}{{{\upsilon }_{w}}{{x}_{w}}\mu }\left( {{p}_{i}}-p \right),{{\lambda }_{D}}=\frac{k\lambda }{{{\upsilon }_{w}}\mu } \\ \end{align}$

式中， x _w为距离常量； x _D为无因次距离； t _D为无因次时间；$\delta$为无因次动边界； P _D为无因次地层压力； $\lambda$ _D为无因次启动压力梯度.

式(5)~式(10)可转化为

$\frac{{{\partial }^{2}}{{P}_{D}}}{\partial {{x}_{D}}^{2}}=\frac{\partial {{P}_{D}}}{\partial {{t}_{D}}},0\le {{x}_{D}}\le \delta \left( {{t}_{D}} \right)$

(11)

$\delta \left( 0 \right) = 0$

(12)

$\left. {P_D } \right|_{t_D = 0} = 0 $

(13)

$\left. {P_D } \right|_{x_D = 0} = f(t_D )$

(14)

$\left. {P_D } \right|_{x_D = \delta \left( {t_D } \right)} = 0$

(15)

$\left. {\dfrac{\partial P_D }{\partial x_D }} \right|_{x_D = \delta \left( {t_D } \right)} = - \lambda _D$

(16)

式(11)~式(16)共同构成了内边界变压力下、考虑启动压力梯度的一维低渗透多孔介质非达西渗流的无因次动边界数学模型.

引入内边界变压力函数如下

$f({{t}_{D}})={{U}_{0}}\sqrt{{{t}_{D}}}\text{ }$

(17)

式中，U ₀为无因次的生产压力数据拟合正参数.

事实上，相对于生产过程中内边界定生产压力较难控制，内边界变生产压力式(17)更符合油井生产的实际情况；特别是生产初期，随着生产时间的增加，井底压力不断降低，无因次井底压力不断增加，试井解释表明一维流动^[23](例如有限导流裂缝井的线性流阶段)无因次压力关于无因次时间的双对数曲线斜率为1/2；而且，式(17)还可以让数学模型保持较好的相似特征.该内边界压力条件独立于与压力梯度相关的非达西运动方程式(3).

1.2 动边界移动规律

由式(15)，可得

$P_D \left( {\delta \left( {t_D } \right)\,,\ \ t_D } \right) = 0$

(18)

对式(18)两边关于时间t_D进行求导，可得

$\left. {\dfrac{\partial P_D }{\partial t_D }} \right|_{x_D = \delta \left( {t_D } \right)} + \left. {\dfrac{\partial P_D }{\partial x_D }} \right|_{x_D = \delta \left( {t_D } \right)} \cdot \dfrac{\partial \delta }{\partial t_D } = 0$

(19)

上式动边界处无因次压力关于无因次距离的导数为左导数.

将动边界条件即式(16)代入式(19)，可得

$\left. {\dfrac{\partial P_D }{\partial t_D }} \right|_{x_D = \delta \left( {t_D } \right)} = \lambda _D \cdot \dfrac{\partial \delta }{\partial t_D }$

(20)

式(11)两边取x _D= $\delta$ (t _D)，可得

${{\left. \frac{\partial {{P}_{D}}}{\partial {{t}_{D}}} \right|}_{{{x}_{D}}=\delta \left( {{t}_{D}} \right)}}={{\left. \frac{{{\partial }^{2}}{{P}_{D}}}{\partial x_{D}^{2}} \right|}_{{{x}_{D}}=\delta \left( {{t}_{D}} \right)}}$

(21)

上式动边界处无因次压力关于无因次距离的二阶导数为二阶左导数.

由式(20)和式(21)联立，可得动边界移动速度公式为^[17,19,20]

$\frac{\partial \delta }{\partial {{t}_{D}}}=\frac{1}{{{\lambda }_{D}}}{{\left. \cdot \frac{{{\partial }^{2}}{{P}_{D}}}{\partial x_{D}^{2}} \right|}_{{{x}_{D}}=\delta \left( {{t}_{D}} \right)}}$

(22)

从式(22)可以看出，考虑启动压力梯度的一维低渗透多孔介质非达西渗流的动边界移动速度与无因次启动压力梯度成反比，而与动边界上无因次压力关于空间距离的二次导数成正比；这与传热学中的经典斯蒂芬动边界问题^[22]有明显区别(动边界速度与一阶导数相关).

2 精确解析求解:相似变量变换方法

通过对比传热学中用于研究单相融化斯蒂芬动边界问题精确解析解的相似变量变换方法，相应地引入如下相似变量^[20,21,22]

$\eta = \dfrac{P_D }{2\sqrt{t_D }}$

(23)

$\xi = \dfrac{x_D }{2\sqrt{t_D }}$

(24)

$\theta = \dfrac{\delta }{2\sqrt{t_D }}$

(25)

于是，由式(23)~式(25)，P_D关于无因次时间和无因次空间的导数可转化为^[20]

$\dfrac{\partial P_D }{\partial t_D } = t_D^{- \tfrac{1}{2}}\eta - \dfrac{x_D }{2}t_D^{-1}\dfrac{\partial \eta }{\partial \xi }$

(26)

$\dfrac{\partial P_D }{\partial x_D }=2 t_D^{\tfrac{1}{2}} \dfrac{\partial \eta }{\partial \xi }\dfrac{1}{2\sqrt{t_D}} = \dfrac{\partial \eta }{\partial \xi }$

(27)

$\dfrac{\partial ^2P_D }{\partial x_D ^2} = \dfrac{\partial ^2\eta }{\partial \xi ^2}\dfrac{1}{2\sqrt{t_D}}$

(28)

利用式(23)~式(28)，式(11)可转化为

$\dfrac{\partial ^2\eta }{\partial \xi ^2}\dfrac{1}{2\sqrt{t_D}} = t_D ^{ - \tfrac{1}{2}}\eta - \dfrac{x_D }{2}t_D^{- 1}\dfrac{\partial \eta }{\partial \xi }\,,\ \ 0 \leqslant \xi \leqslant \theta$

(29)

消去式(29)方程两边相同的变量，可进一步简化为

$\dfrac{1}{2}\dfrac{\partial ^2\eta }{\partial \xi ^2} + \xi \cdot \dfrac{\partial \eta }{\partial \xi } - \eta = 0\,,\ \ 0 \leqslant \xi \leqslant \theta$

(30)

式(14)~式(16)可转化为

${{\left. \eta \right|}_{\xi }}=0=\frac{{{U}_{0}}}{2}$

(31)

$\left. \eta \right|_\xi = \theta = 0$

(32)

$\left. {\dfrac{\partial \eta }{\partial \xi }} \right|_\xi = \theta = - \lambda _D $

(33)

通过引入相似变量即式(23)~式(25)，考虑启动压力梯度的一维低渗透多孔介质非达西渗流的无因次动边界数学模型，可转为由式(30)~式(33)共同构成的一个封闭的常微分方程组，且具有固定的定边界条件，较容易进行解析求解.

常微分方程组的变量包括函数$ \eta (\xi )$和$\theta$ ；首先，由方程式(30)~式(32)构成的线性常微分方程组，可求得该方程组的精确解析解为

$\eta (\xi )=\frac{{{U}_{0}}}{2}\left( {{e}^{-{{\xi }^{2}}}}+\text{ }{{\pi }^{\tfrac{1}{2}}}\xi erf\left( \xi \right)-\frac{\xi }{{{e}^{{{\theta }^{2}}}}\theta }-\text{ }{{\pi }^{\tfrac{1}{2}}}\xi erf\left( \theta \right) \right)$

(34)

将式(34)两边关于$ \xi $ 变量求导，可得

${\eta }'\left( \xi \right)=\frac{{{U}_{0}}}{2}\left( \text{ }{{\pi }^{\tfrac{1}{2}}}erf\left( \xi \right)-\frac{1}{{{e}^{{{\theta }^{2}}}}\theta }-\text{ }{{\pi }^{\tfrac{1}{2}}}erf(\theta ) \right)$

(35)

由式(35)和式(33)联立，可得

${\eta }'\left( \theta \right)=\frac{{{U}_{0}}}{2}\left( \text{ }{{\pi }^{\tfrac{1}{2}}}erf\left( \theta \right)-\frac{1}{{{e}^{{{\theta }^{2}}}}\theta }-\text{ }{{\pi }^{\tfrac{1}{2}}}erf\left( \theta \right) \right)=-{{\lambda }_{D}}$

(36)

化简式(36)可得

${e}^{\theta ^2}\theta = \dfrac{U_0 }{2}\dfrac{1}{\lambda _D }$

(37)

由式(37)可以看出，由于${e}^{\theta ^2}\theta $的严格单调性，对于给定的无因次启动压力梯度$\lambda_D$和U₀>0的值，\theta 的正值可由超越代数 方程式(37)惟一确定.

于是，由定义的相似变量即式(23) \sim 式(25)，可得内边界变压力下、考虑启动压力梯度的一维低渗透多孔介质非达 西渗流动边界模型的精确解析解为

$\begin{align} & {{P}_{D}}=\frac{{{U}_{0}}}{2}\left[2t_{D}^{\tfrac{1}{2}}{{e}^{-\tfrac{x_{D}^{2}}{4{{t}_{D}}}}}+{{\pi }^{\tfrac{1}{2}}}{{x}_{D}}erf\left( \frac{{{x}_{D}}}{2\sqrt{{{t}_{D}}}} \right)- \right.\left. \frac{{{x}_{D}}}{{{e}^{{{\theta }^{2}}}}\theta }-{{\pi }^{\tfrac{1}{2}}}{{x}_{D}}erf\left( \theta \right) \right],\\ & {{x}_{D}}\in \left[0,\delta \right] \\ \end{align}$

(38)

其中，erf为误差函数； $\theta$ >0的值由式(37)确定. 无因次动边界 $\delta$ 的表达式为

$\delta = 2 \theta \sqrt {t_D }$

(39)

由 $\theta$ 正值的惟一性证明了一维低渗透多孔介质渗流动边界模型存在惟一的精确解析解.

此外，为与达西渗流模式下的解析解结果进行对比，进而分析考虑启动压力梯度的多孔介质非达西渗流模型考虑动边界条件的必要性，内边界变压力下、一维多孔介质达西渗流的无因次动边界数学模型的精确解析解也可通过相似变量法求出.

达西渗流模式下，无因次启动压力梯度 $\lambda$_D =0，动边界条件不再有效；其外边界条件应为

$\left. {P_D } \right|_{x_D = \infty } = 0$

(40)

于是，式(11)、式(13)、式(14)和式(40)共同构成了内边界变压力下、一维多孔介质达西渗流的无因次数学模型. 采用相同的相似变量式(23)~式(24)，该模型进行转化；式(11)和式(14)分别转为式(30)和式(31)；式(13)和式(40)可转化为

$\left. \eta \right|_\xi = \infty = 0$

(41)

于是由式(30)、式(31)和式(41)共同构成的封闭线性常微分方程组解析解可表达为

$\eta \left( \xi \right)=\frac{{{U}_{0}}}{2}\left( {{e}^{-{{\xi }^{2}}}}+{{\pi }^{\tfrac{1}{2}}}\xi erf\left( \xi \right)-{{\pi }^{\tfrac{1}{2}}}\xi \right)$

(42)

进而由相似变量式(23)~式(24)，可以求得内边界变压力下、一维多孔介质达西渗流无因次数学模型的精确解析解为

$\begin{align} & {{P}_{D}}\left( {{x}_{D}},{{t}_{D}} \right)=\frac{{{U}_{0}}}{2}\left[2t_{D}^{\tfrac{1}{2}}{{e}^{-\tfrac{x_{D}^{2}}{4{{t}_{D}}}}}+ \right.\left. {{\pi }^{\tfrac{1}{2}}}{{x}_{D}}erf\left( \frac{{{x}_{D}}}{2\sqrt{{{t}_{D}}}} \right)-{{\pi }^{\tfrac{1}{2}}}{{x}_{D}} \right],\\ & {{x}_{D}}\in [0,\infty ) \\ \end{align}$

(43)

非达西渗流动边界模型与达西渗流模型的精确解析解是统一的. 由式(37)可知，当 $\lambda$_D值趋于0 (即达西渗流)时，$\theta$ 值将趋于+ $\infty$ ，而由式(39)，$\delta$ 也趋于+$\infty$ .

式(38)的两边令 $\theta$ 趋于+$\infty$ ，可得

$\begin{gathered} amp;\mathop {\lim }\limits_{\theta \to + \infty } {\mkern 1mu} {P_D} = \mathop {\lim }\limits_{\theta \to + \infty } {\mkern 1mu} \left\{ \begin{gathered} amp;\frac{{{U_0}}}{2}[2t_D^{\tfrac{1}{2}}{e^{ - \tfrac{{x_D^2}}{{4{t_D}}}}} + {\pi ^{\tfrac{1}{2}}}{x_D}erf\left( {\frac{{{x_D}}}{{2\sqrt {{t_D}} }}} \right) \hfill \\ amp;\left. { - \frac{{{x_D}}}{{{e^{{\theta ^2}}}\theta }} - {\pi ^{\tfrac{1}{2}}}{x_D}erf\left( \theta \right)]} \right\} \hfill \\ \end{gathered} \right. \hfill \\ amp; = \frac{{{U_0}}}{2}[2t_D^{\tfrac{1}{2}}{e^{ - \tfrac{{x_D^2}}{{4{t_D}}}}} + {\pi ^{\tfrac{1}{2}}}{x_D}erf\left( {\frac{{{x_D}}}{{2\sqrt {{t_D}} }}} \right) - \hfill \\ amp;{x_D}\mathop {\lim }\limits_{\theta \to + \infty } {\mkern 1mu} \frac{1}{{{e^{{\theta ^2}}}\theta }} - {\pi ^{\tfrac{1}{2}}}{x_D}\mathop {\lim }\limits_\theta {\mkern 1mu} \to + \infty erf\left( \theta \right)] \hfill \\ \end{gathered} $

(44)

由于

$\underset{\theta \to +\infty }{\mathop{\lim }}\,\frac{1}{{{e}^{{{\theta }^{2}}}}\theta }=0$

(45)

$\underset{\theta \to +\infty }{\mathop{\lim }}\,erf\left( \theta \right)=1$

(46)

于是

$\underset{\theta \to +\infty }{\mathop{\lim }}\,{{P}_{D}}=\frac{{{U}_{0}}}{2}\left[2t_{D}^{\tfrac{1}{2}}{{e}^{-\tfrac{x_{D}^{2}}{4{{t}_{D}}}}}+{{\pi }^{\tfrac{1}{2}}}{{x}_{D}}erf\left( \frac{{{x}_{D}}}{2\sqrt{{{t}_{D}}}} \right)-{{\pi }^{\tfrac{1}{2}}}{{x}_{D}} \right]$

(47)

式(47)与式(43)保持一致，即当启动压力梯度值趋于零时，非达西渗流动边界模型的精确解析解将退化为达西渗流情况下的精确解析解.

3 数值求解: 基于空间坐标变换的有限差分方法

由于数学模型中动边界的存在，一维流动区域并不是固定的，而是随着时间逐渐向外扩展；在数值求解时，为克服瞬时流动区域网格划分的难题，引入如下空间坐标变换，可将动边界的数学模型首先转化为定边界的数学模型^[17,19]

${{y}_{D}}\left( {{x}_{D}},{{t}_{D}} \right)=\frac{{{x}_{D}}}{\delta \left( {{t}_{D}} \right)},{{x}_{D}}\in \left[0,\delta \left( {{t}_{D}} \right) \right]$

(48)

通过x_D到y_D的坐标变换即式(48)，可得

${{y}_{D}}\left( 0,{{t}_{D}} \right)=\frac{0}{\delta \left( {{t}_{D}} \right)}=0$

(49)

${{y}_{D}}\left( \delta \left( {{t}_{D}} \right),{{t}_{D}} \right)=\frac{\delta \left( {{t}_{D}} \right)}{\delta \left( {{t}_{D}} \right)}=1$

(50)

由式(48)~式(50)可知，动边界模型瞬时流动区域[0,$\delta$ (t_D)]可以转化到固定流动区域[0,1]上.由式(48)，无因次数学模型即式(11) \sim 式(16)中的微分变量可依次转化为^[17,19]

$\frac{\partial {{P}_{D}}}{\partial {{x}_{D}}}=\frac{\partial {{P}_{D}}}{\partial {{y}_{D}}}\cdot \frac{1}{\delta }$

(51)

$\frac{{{\partial }^{2}}{{P}_{D}}}{\partial x_{D}^{2}}=\frac{{{\partial }^{2}}{{P}_{D}}}{\partial y_{D}^{2}}\cdot \frac{1}{{{\delta }^{2}}}$

(52)

$\begin{align} & \frac{\partial {{P}_{D}}}{\partial {{t}_{D}}}=\frac{\partial {{P}_{D}}}{\partial {{t}_{D}}}+\frac{\partial {{P}_{D}}}{\partial {{y}_{D}}}\cdot \left( -\frac{{{x}_{D}}}{{{\delta }^{2}}} \right)\cdot \frac{\partial \delta }{\partial {{t}_{D}}}= \\ & 1cm\frac{\partial {{P}_{D}}}{\partial {{t}_{D}}}-\frac{\partial {{P}_{D}}}{\partial {{y}_{D}}}\cdot \frac{\partial \delta }{\partial {{t}_{D}}}\cdot \frac{{{y}_{D}}}{\delta } \\ \end{align}$

(53)

将式(51)~式(53)代入式(11)可得

$\frac{{{\partial }^{2}}{{P}_{D}}}{\partial y_{D}^{2}}\cdot \frac{1}{{{\delta }^{2}}}=\frac{\partial {{P}_{D}}}{\partial {{t}_{D}}}-\frac{\partial {{P}_{D}}}{\partial {{y}_{D}}}\cdot \frac{\partial \delta }{\partial {{t}_{D}}}\cdot \frac{{{y}_{D}}}{\delta },0\le {{y}_{D}}\le 1$

(54)

将式(51)~式(53)代入式(12) \sim 式(16)和式(22)可得

$\left. {P_D \left( {y_D ,t_D } \right)} \right|_{t_D = 0} = 0$

(55)

$\left. {P_D } \right|_{y_D = 0} = U_0 \sqrt {t_D }$

(56)

${{\left. \frac{\partial {{P}_{D}}}{\partial {{y}_{D}}}\cdot \frac{1}{\delta } \right|}_{{{y}_{D}}=1}}=-{{\lambda }_{D}}$

(57)

$\left. {P_D } \right|_{y_D = 1} = 0$

(58)

$\dfrac{\partial \delta }{\partial t_D } = \dfrac{1}{\lambda _D } \cdot \left. {\dfrac{\partial ^2P_D }{\partial y_D ^2}} \right|_{y_D = 1} \cdot \dfrac{1}{\delta ^2}$

(59)

将式(59)代入式(54)可得

$\frac{{{\partial }^{2}}{{P}_{D}}}{\partial y_{D}^{2}}\cdot \delta =\frac{\partial {{P}_{D}}}{\partial {{t}_{D}}}\cdot {{\delta }^{3}}-\frac{\partial {{P}_{D}}}{\partial {{y}_{D}}}\cdot \frac{{{y}_{D}}}{{{\lambda }_{D}}}\cdot {{\left. \frac{{{\partial }^{2}}{{P}_{D}}}{\partial y_{D}^{2}} \right|}_{{{y}_{D}}=1}}$

(60)

由式(57)可得

$\delta =-{{\left. \frac{\partial {{P}_{D}}}{\partial {{y}_{D}}} \right|}_{{{y}_{D}}=1}}\cdot \frac{1}{{{\lambda }_{D}}}$

(61)

将式(61)代入变换后的控制方程式即式(60)，可得

$\frac{{{\partial }^{2}}{{P}_{D}}}{\partial y_{D}^{2}}\cdot {{\left. \frac{\partial {{P}_{D}}}{\partial {{y}_{D}}} \right|}_{{{y}_{D}}=1}}-\frac{\partial {{P}_{D}}}{\partial {{t}_{D}}}\cdot {{\left( {{\left. \frac{\partial {{P}_{D}}}{\partial {{y}_{D}}} \right|}_{{{y}_{D}}=1}} \right)}^{3}}\cdot \frac{1}{\lambda _{D}^{2}}-\frac{\partial {{P}_{D}}}{\partial {{y}_{D}}}\cdot {{y}_{D}}\cdot {{\left. \frac{{{\partial }^{2}}{{P}_{D}}}{\partial y_{D}^{2}} \right|}_{{{y}_{D}}=1}}=0$

(62)

式(62)、式(55)、式(56)和式(58)共同组成了内边界变生产压力情况下、考虑启动压力梯度的一维低渗透多孔介质非达西渗流无因次动边界数学模型等价变换后的定边界偏微分方程组.由式(62)可以看出，偏微分方程具有很强的非线性，间接反映出了原动边界问题的强非线性.鉴于变换后模型的强非线性，采用全隐式的有限差分方法进行数值求\linebreak解^[19]；其中，一阶导数采用一阶向前差分格式，二阶导数采用二阶中心差分格式.

控制方程式(62)的全隐式有限差分方程可表示为

$\begin{align} & -\frac{P_{Di+1}^{j+1}-2P_{Di}^{j+1}+P_{Di-1}^{j+1}}{{{\left( \Delta {{y}_{D}} \right)}^{2}}}\cdot \frac{P_{DN-1}^{j+1}}{\Delta {{y}_{D}}}+ \\ & \frac{P_{Di}^{j+1}-P_{Di}^{j}}{\Delta {{t}_{D}}}\cdot {{\left( \frac{P_{DN-1}^{j+1}}{\Delta {{y}_{D}}} \right)}^{3}}\cdot \frac{1}{\lambda _{D}^{2}}- \\ & \frac{P_{Di+1}^{j+1}-P_{Di}^{j+1}}{\Delta {{y}_{D}}}\cdot i\cdot \Delta {{y}_{D}}\cdot \frac{P_{DN-2}^{j+1}-2P_{DN-1}^{j+1}}{{{\left( \Delta {{y}_{D}} \right)}^{2}}}=0\left( i=1,2,\cdots ,N-2 \right) \\ \end{align}$

(63)

式(58)的差分方程为

$P_{{D}N}^{j + 1} = 0$

(64)

由式(64)，第( N -1)个空间网格所对应控制方程的差分方程为

$\begin{align} & -\frac{-2P_{DN-1}^{j+1}+P_{DN-2}^{j+1}}{{{\left( \Delta {{y}_{D}} \right)}^{2}}}\cdot \frac{P_{DN-1}^{j+1}}{\Delta {{y}_{D}}}+\frac{P_{DN-1}^{j+1}-P_{DN-1}^{j}}{\Delta {{t}_{D}}}\cdot \\ & {{\left( \frac{P_{DN-1}^{j+1}}{\Delta {{y}_{D}}} \right)}^{3}}\cdot \frac{1}{\lambda _{D}^{2}}+ \\ & P_{DN-1}^{j+1}\cdot \left( N-1 \right)\cdot \frac{P_{DN-2}^{j+1}-2P_{DN-1}^{j+1}}{{{\left( \Delta {{y}_{D}} \right)}^{2}}}=0 \\ \end{align}$

(65)

式(56)的差分方程为

$P_{{D}0}^{j + 1} = U_0 \sqrt {\left( {j + 1} \right)\Delta t_D }$

(66)

由式(55)得到无因次地层压力的初始值，差分格式为

$P_{{D}i}^0 = 0\,,\ \ \ i = 0,1,2,\cdots,N-1$

(67)

由式(63)、式(65)和式(66)共同构成了第j+1时间层的N个非线性差分方程，含有N个未知变 量$P_{{D}i}^{j + 1} \left( {i = 0,1,2,\cdots,N - 1} \right)$；采用牛顿\!-\!-\!拉夫逊迭代方法^[22]进行数值求解.

4 数值解与精确解析解对比

通过上述数值求解方法对内边界变生产压力下、考虑启动压力梯度的一维低渗透多孔介质非达西渗流的无因次动边界数学模型进行数值求解；其中，空间步长与时间步长设定为\Delta y_D=1/160，$\Delta$ t_D=10.为了对比分析达西渗流的模型解，即 $\lambda$_D=0时，利用内边界变生产压力、一维半无限长多孔介质达西渗流无因次数学模型的精确解析解即式(43)进行计算.

图1和图2分别为U₀=1.0时、不同无因次启动压力梯度 $\lambda$_D=0.852,$\lambda$_D=0.553,$\lambda$_D=0.242\linebreak下^[2]的一维多孔介质非达西渗流的无因次动边界模型数值解与精确解析解关于瞬时无因次动边界距离和t_D=5 000时刻地层压力分布的对比曲线.

图3和图4分别为 $\lambda$_D=0.553时,U₀=0.8,U₀=1.0 ,U₀=1.2,U₀=1.4下的一维多孔介 质非达西渗流无因次动边界模型数值解与精确解析解关于瞬时无因次动边界距离和t_D=5 000,时刻地层压力分布的对比曲线.

由图1~图4可以看出，求解内边界变生产压力下、考虑启动压力梯度的一维低渗透多孔介质非达西渗流动边界模型的数值解与精确解析解吻合程度较好；由此严格验证了动边界模型数值求解方法的正确性和有效性.

从图1和图2还可看出,启动压力梯度对一维多孔介质渗流模型解产生较大影响;无因次启动压力梯度越大,达西渗流模型与考虑启动压力梯度的动边界数学模型解的差异越大.此外,从图2还可看出,非零无因次启动压力梯度值对应的地层压力分布在动边界以内大于零,而动边界以外地层压力保持为零值,曲线表现出紧支性特点,其与达西渗流模型对应的地层压力分布曲线明显不同,达西渗流模型压力降可瞬时传至无穷远;考虑动边界条件与不考虑动边界条件对模型的计算结果影响较大(对于一维渗流模型,考虑启动压力梯度与不考虑启动压力梯度模型的控制方程是相同的^[3],区别在于动边界条件);因此,研究低渗透多孔介质中非稳态渗流问题时,应该考虑动边界的影响,尤其致密多孔介质(例如页岩油、致密油等)中启动压力梯度值较大时.

5 结论

(1)考虑启动压力梯度的一维低渗透多孔介质非达西渗流的动边界移动速度与无因次启动压力梯度成反比，而与动边界上无因次压力关于空间距离的二次导数成正比；其与传热学中的经典斯蒂芬动边界问题有明显区别(动边界速度与一阶导数相关).

(2)该一维低渗透多孔介质渗流动边界模型存在惟一的精确解析解；精确解析解可以用来严格验证模型数值求解方法的正确性和有效性.另外，当启动压力梯度值趋于零时，非达西渗流动边界模型的精确解析解将退化为达西渗流情况下的精确解析解.

(3)考虑启动压力梯度的一维低渗透多孔介质非达西渗流动边界模型的地层压力分布曲线表现出紧支性特点，其与达西渗流模型的有显著不同；研究低渗透多孔介质中非稳态渗流问题时，应该考虑动边界的影响，尤其针对致密多孔介质(例如页岩油、致密油等)启动压力梯度值较大时的情况.

(4)针对于封闭边界的情况，在动边界到达封闭边界以后，可用定边界的渗流模型进行描述，模型初始条件则为动边界到达封闭边界时刻的空间压力分布；该压力分布可通过本文所求得的动边界模型精确解析解求出.