«上一篇
文章快速检索 高级检索
下一篇»
  力学学报  2015, Vol. 47 Issue (4): 605-612  DOI: 10.6052/0459-1879-14-385
0

研究论文

引用本文 [复制中英文]

刘文超, 姚军, 陈掌星, 刘曰武, 孙海. 低渗透多孔介质渗流动边界模型的解析与数值解[J]. 力学学报, 2015, 47(4): 605-612. DOI: 10.6052/0459-1879-14-385.
[复制中文]
Liu Wenchao, Yao Juny, Chen Zhangxin, Liu Yuewu, Sun Haiy. RESEARCH ON ANALYTICAL AND NUMERICAL SOLUTIONS OF A MOVING BOUNDARY MODEL OF SEEPAGE FLOWIN LOW-PERMEABLE POROUS MEDIA[J]. Chinese Journal of Ship Research, 2015, 47(4): 605-612. DOI: 10.6052/0459-1879-14-385.
[复制英文]

基金项目

国家自然科学基金项目(51404232)、国家留学基金委员会项目、中国博士后科学基金项目(2014M561074) 和国家科技重大专项(2011ZX05038-003) 资助.

作者简介

姚军, 教授, 主要研究方向:油气渗流力学及油气藏工程. E-mail: RCOGFR_UPC@126.com

文章历史

2014-12-04收稿
2015-03-30录用
2015-04-21网络版发表
Contents              Abstract              Full text              Figures/Tables              PDF
低渗透多孔介质渗流动边界模型的解析与数值解
刘文超1, 2, 3, 姚军2, 陈掌星3, 刘曰武1, 孙海2    
1. 中国科学院力学研究所流固耦合系统力学重点实验室, 北京100190;
2. 中国石油大学(华东), 青岛266580;
3. 卡尔加里大学, 卡尔加里T2N 1N4, 加拿大
2014-12-04 收稿，2015-03-30 录用，2015-04-21 网络版发表
基金项目：国家自然科学基金项目(51404232)、国家留学基金委员会项目、中国博士后科学基金项目(2014M561074) 和国家科技重大专项(2011ZX05038-003) 资助.
作者简介：姚军, 教授, 主要研究方向:油气渗流力学及油气藏工程. E-mail: RCOGFR_UPC@126.com
摘要：考虑启动压力梯度的低渗透多孔介质非达西渗流模型属于强非线性动边界问题, 分别利用相似变量变换方法和基于空间坐标变换的有限差分方法, 对内边界变压力情况下、考虑启动压力梯度的一维低渗透多孔介质非达西渗流动边界模型进行了精确解析与数值求解研究. 研究结果表明:该动边界模型存在唯一的精确解析解, 且所求得的精确解析解可严格验证数值解的正确性;且当启动压力梯度值趋于零时, 非达西渗流动边界模型的精确解析解将退化为达西渗流情况下的精确解析解. 由求解结果作出的非零无因次启动压力梯度下的地层压力分布曲线表现出紧支性特点, 其与达西渗流模型的有显著不同. 因此, 研究低渗透多孔介质中非稳态渗流问题时, 应该考虑动边界的影响. 研究内容完善了低渗透多孔介质的非达西渗流力学理论, 为低渗透油气藏开发的试井解释与油藏数值模拟技术提供了理论基础.
关键词： 低渗透油藏    启动压力梯度    动边界    精确解析解    数值解    
引言

考虑启动压力梯度的多孔介质非达西渗流模型研究[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为地层压力；pi为原始地层压力；Cf为流体压缩系数.

岩石孔隙度的状态方程为[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)
式中，Ct为综合压缩系数.

初始条件为

$\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}$
式中，xw为距离常量；xD为无因次距离；tD为无因次时间；$\delta$为无因次动边界；PD为无因次地层压力； $\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)
式中，U0为无因次的生产压力数据拟合正参数.

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

1.2 动边界移动规律

由式(15)，可得

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

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

$\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)两边取xD = $\delta$ (tD)，可得
${{\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$和U0 >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)

通过xD到yD的坐标变换即式(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$ (tD)]可以转化到固定流动区域[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 yD =1/160，$\Delta$ tD =10.为了对比分析达西渗流的模型解，即 $\lambda$D=0时，利用内边界变生产压力、一维半无限长多孔介质达西渗流无因次数学模型的精确解析解即式(43)进行计算.

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

图 1 不同无因次启动压力梯度下的瞬时无因次动边界距离对比 Fig.1 Comparison of transient distance of moving boundary under different values of dimensionless threshold pressure gradient
图 2 不同无因次启动压力梯度下的地层压力分布对比 Fig.2 Comparison of formation pressure distribution under different values of dimensionless threshold pressure gradient

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

图 3 不同U0值下的瞬时无因次动边界距离对比 Fig.3 Comparison of transient distance of moving boundary under different values of U0
图 4 不同U0值下的地层压力分布对比 Fig.4 Comparison of formation pressure distribution under different values of U0

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

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

5 结论

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

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

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

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

参考文献
[1] 黄延章. 低渗透油层非线性渗流特征. 特种油气藏, 1997, 4(1): 9-14 (Huang Yanzhang. Nonlinear percolation feature in low permeability reservoir. Special Oil & Gas Reservoirs, 1997, 4(1): 9-14 (in Chinese))
[2] Prada A, Civan F. Modification of Darcy's law for the threshold pressure gradient. Journal of Petroleum Science and Engineering, 1999, 22(4): 237-240
[3] 孔祥言. 高等渗流力学. 中国科学技术大学出版社, 2010 (Kong Xiangyan. Advanced Mechanics of Fluid in Porous Media. Hefei: Press of University of Science and Technology of China, 2010 (in Chinese))
[4] Zhu WY, Song HQ, Huang XH, et al. Pressure characteristics and effective deployment in a water-bearing tight gas reservoir with low-velocity non-Darcy flow. Energy & Fuels, 2011, 25(3): 1111-1117
[5] 姚同玉, 黄延章, 李继山. 孔隙介质中稠油流体非线性渗流方程.力学学报, 2012, 44(1): 106-110 (Yao Tongyu, Huang Yanzhang, Li Jishan. Nonliner flow equations for heavy oil in porous media. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(1): 106-110 (in Chinese))
[6] Cai JC. A fractal approach to low velocity non-Darcy flow in a low permeability porous medium. Chinese Physics B, 2014, 23(4): Article ID 044701
[7] Pascal H. Nonsteady flow through porous media in the presence of a threshold gradient. Acta Mechanica, 1981, 39(3-4): 207-224
[8] Wu YS, Pruess K, Witherspoon PA. Flow and displacement of Bingham non-Newtonian fluids in porous media. SPE Reservoir Engineering, 1993, 7(3): 369-376
[9] 祝春生, 程林松, 张淑娟. 低渗透油藏不稳定渗流的近似解法.西南石油大学学报(自然科学版), 2008, 30(4): 69-72 (Zhu Chunsheng, Cheng Linsong, Zhang Shujuan. Approximate solution for non-steady flow in low-permeability reservoir. Journal of Southwest Petroleum University (Science m & Technology Edition), 2008, 30(4): 69-72 (in Chinese))
[10] 程时清, 张盛宗, 黄延章等. 低速非达西渗流动边界问题的积分解.力学与实践, 2002, 24(3): 15-17 (Cheng Shiqing, Zhang Zongsheng, Huang Yanzhang, et al. An integral solution of free-boundary problem of non-Darcy flow behavior. Mechanics in Engineering, 2002, 24(3): 15-17 (in Chinese))
[11] Ghedan S, Lu J. Pressure behavior of vertical wells in low permeability reservoirs with threshold pressure gradient. Special Topics and Reviews in Porous Media-An International Journal, 2011, 2(3): 157-169
[12] Lu J. Pressure behavior of uniform flux hydraulically fractured wells in low permeability reservoirs with threshold pressure gradient. Special Topics and Reviews in Porous Media-An International Journal, 2012, 3(4): 307-320
[13] 郭永存, 曾清红, 王仲勋. 动边界低渗透油藏的无网格方法数值模拟.工程力学, 2006, 23(11): 188-192 (Guo Yongcun, Zeng Qinghong, Wang Zhongxun. Numerical simulation of low permeability flow with moving-boundary using meshless methods. Engineering Mechanics, 2006, 23(11): 188-192 (in Chinese))
[14] Li FH, Liu CQ. Pressure transient analysis for unsteady porous flow with start-up pressure derivative. Well Testing, 1997, 6: 1-11
[15] Song FQ, Liu CQ, Li F H. Transient pressure of percolation through one dimension porous media with threshold pressure gradient. Applied Mathematics and Mechanics, 1999, 20(1): 27-35
[16] 郭永存, 卢德唐, 马凌霄.低渗透油藏渗流的差分法数值模拟.水动力学研究与进展, 2004, 19(3): 288-293 (Guo Yongcun, Lu Detang, Ma Lingxiao. Numerical simulation of fluid flow in low permeability reservoir using finite difference method. Journal of Hydrodynamics, 2004, 19(3): 288-293 (in Chinese))
[17] 姚军, 刘文超. 低渗透油藏非达西渗流动边界模型求解新方法.力学季刊, 2012, 33(4): 597-601 (Yao Jun, Liu Wenchao. New method for solution of the model of non-Darcy seepage flow in low-permeability reservoirs with moving boundary. Chinese Quarterly of Mechanics, 2012, 33(4): 597-601 (in Chinese))
[18] 刘文超, 姚军, 孙致学等. 基于渗透率连续变化的低渗透多孔介质非线性渗流模型.计算力学学报, 2012, 29(6): 885-892 (Liu Wenchao, Yao Jun, Sun Zhixue, et al. Model of nonlinear seepage flow in low-permeability porous media based on the permeability gradual recovery. Chinese Journal of Computational Mechanics, 2012, 29(6): 885-892 (in Chinese))
[19] Yao J, Liu W C, Chen ZX. Numerical solution of a moving boundary problem of one-dimensional flow in semi-infinite long porous media with threshold pressure gradient. Mathematical Problems in Engineering, 2013, Article ID 384246
[20] Liu W, Yao J, Wang Y. Exact analytical solutions of moving boundary problems of one-dimensional flow in semi-infinite long porous media with threshold pressure gradient. International Journal of Heat and Mass Transfer, 2012, 55(21-22): 6017-6022
[21] Liu WC, Yao J, Chen ZX, et al. Analytical solution of a double moving boundary problem for nonlinear flows in one-dimensional semi-infinite long porous media with low permeability. Acta Mechanica Sinica, 2014, 30(1): 50-58
[22] Crank J. Free and Moving Boundary Problems. Oxford: Clarendon Press, 1984
[23] 刘能强. 实用现代试井解释方法. 北京:石油工业出版社, 2008 (Liu Nengqiang. Practical Modern Well Testing Interpretation Method. Beijing: Petroleum Industry Press, 2008 (in Chinese))
RESEARCH ON ANALYTICAL AND NUMERICAL SOLUTIONS OF A MOVING BOUNDARY MODEL OF SEEPAGE FLOWIN LOW-PERMEABLE POROUS MEDIA
Liu Wenchao1,2,3, Yao Juny2, Chen Zhangxin3, Liu Yuewu1, Sun Haiy2    
1. Key Laboratory for Mechanics in Fluid Solid Coupling System, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China;
2. China University of Petroleum, Qingdao 266580, China;
3. University of Calgary, Alberta Calgary T2N 1N4, Canada
Fund: The project was supported by the National Natural Science Foundation of China (51404232),China Scholarship Council (CSC),China Postdoctoral Science Foundation (2014M561074) and the National Science and Technology Major Project (2011ZX05038-003).
Abstract: The modes of non-Darcy seepage flow in low-permeable porous media with threshold pressure gradient belong to the moving boundary problems with strong nonlinearity. In this paper, a moving boundary model of one-dimensional non-Darcy seepage flow in low-permeable porous media with threshold pressure gradient is studied for the case of a variable pressure at the inner boundary. A similarity transformation method and a spatial coordinate transformation based finite difference method are applied to obtain the exact analytical solution and numerical solution of the moving boundary model, respectively. Research results show that the moving boundary model has a unique exact analytical solution, which also strictly verifies the accuracy of the numerical solution; and when the value of threshold pressure gradient tends to zero, the exact analytical solution of the moving boundary model of non-Darcy seepage flow can reduce to the one corresponding to Darcy's seepage flow. The formation pressure distribution curves with non-zero dimensionless threshold pressure gradient from these solutions exhibit the characteristics of compact support, which is obviously different from those corresponding to Darcy's flow model. Therefore, the study on the unsteady seepage flow problem in low-permeable porous media should take into account the effect of moving boundary condition. The presented research improves the theory of non-Darcy seepage flow in low-permeable porous media, and also builds theoretical foundation for the technologies of well testing interpretation and reservoir numerical simulation involved in the development of low-permeable oil and gas reservoirs.
Key words: low-permeable reservoir    threshold pressure gradient    moving boundary    exact analytical solution    numerical solution