目前油气藏中普遍应用的试井分析模型, 包括均质模型、双重介质模型和三重介质模型等, 都是基于达西定律建立的连续介质渗流理论模型。Bruce等^[1]建立的单重介质数学模型, 解决了连续多孔介质渗流机理的表征问题。Warren和Root^[2]研究了裂缝-孔隙型储集层中的渗流问题, 建立了Warren-Root模型。Camacho等^[3]对Warren-Root模型进行了扩展, 建立了适用于裂缝-溶洞-基质型碳酸盐岩储集层的三重孔隙介质渗流模型。这些以达西渗流理论和连续性假设为基础的传统渗流力学模型在常规油气藏开发中发挥了重要作用。

然而, 达西渗流理论有一定的适用范围, 超过了这些范围, 渗流规律将变为非线性流动规律^[4]。用于描述非线性流动规律的典型公式有：考虑启动压力梯度的低速非达西公式^[5], 考虑速度平方项的Forchheimer高速非达西公式^[6], 考虑流体流变特性的幂律非牛顿公式^[7], 考虑渗透率随压力变化的指数公式^[8], 考虑均匀流沿程水头损失的Darcy-Weisbach公式^[9]以及考虑空腔自由流的Navier-Stokes方程^[10]等。对于线性流动规律, 除达西定律外, 常用的描述方法还包括泊肃叶定律^[11], 泊肃叶定律描述流体沿细管压力损失。

上述线性和非线性流动规律在整个油藏中根据实际情况可以单独使用, 也可以组合使用。单独使用的情况相对比较简单, 如常规均质砂岩油藏可直接使用达西渗流模型对整个储集层进行表征^[1]。均质稠油油藏可使用幂律模型对整个储集层进行表征^[12]。组合使用的情况相对复杂, 如在先水驱后聚驱的油藏中将形成以注聚井为中心的聚合物、水、原油3个分布区。聚合物区的流动规律用非牛顿流体模型表征, 而水区和原油区则用达西渗流模型表征^[13], 不同区域之间是一种渗流-渗流复合流动。又如在疏松砂岩油藏注水开发过程中形成的大孔道和热采稠油油藏蒸汽驱阶段形成的大孔道, 都同时存在管流和渗流复合流动情况。缝洞型油藏空间形态复杂, 介质尺度跨度大且非均质性严重^{[14, 15]}, 水平井生产过程中, 井筒内管流与地层渗流互相制约、互相影响、互为边界条件, 储集层与井筒之间是一种渗流-管流耦合流动^{[16, 17]}。流体在这类油气藏中的流动规律异常复杂, 既有多孔介质渗流又有大空间的自由流动, 是一个复杂的渗流-自由流耦合流动^[18]。渗流-管流耦合流动及渗流-自由流耦合流动基本都是采用数值方法对模型进行求解。尽管目前已有比较成熟的数值求解算法, 但如果不经过任何的简化处理直接求解这类耦合流动模型, 将存在网格剖分困难、计算量巨大的问题。

针对渗流-管流耦合流动问题：Collins等^[19]和吴淑红等^[20]基于等效渗透率建立了水平井管流的简化模型。在水利地质工程领域, 陈崇希等^[21]基于等效渗透系数建立了渗流-管流耦合模型; 陈崇希^[22]在渗流-管流耦合模型的基础上, 提出岩溶管道-裂隙-孔隙三重空隙介质地下水流模型; 赵延林等^[23]基于折算渗透系数建立承压溶洞突水的非线性渗流-管流耦合模型。这类方法的核心思想是基于Darcy-Weisbach公式^[9]定义渗流、层流管流和紊流管流3种情况的等效渗透率或者等效渗透系数, 实现渗流、层流管流和紊流管流的统一。其优点是可以将储集层中的渗流区、层流区和紊流区的流体流动规律用一个统一的运动方程来描述, 有利于模型方程的构建及求解。但长期以来很少用于地下流体管流（自由流）-渗流耦合流动模拟, 特别是试井分析领域未见文献报道。

针对渗流-自由流耦合流动问题：万义钊等^[24]利用高渗透率、高孔隙度的区域（块）来描述大尺度的裂缝或溶洞。流体在裂缝或溶洞中的流动仍然符合达西渗流, 由此构建了缝洞型油藏的半解析解流动模拟模型。这种方法的优点是模型求解速度快, 适合大规模油藏计算, 缺点是洞穴计算精度大大降低。段宝江等^[25]认为可将溶洞视为等势体, 在此假设条件下建立了缝、洞试井分析模型, 其中流体在裂缝中的流动为拟稳态窜流。Wu Yonghui等^[26]将缝洞型油藏划分为岩块系统（包括基岩、微裂缝和微小溶洞）、裂缝、溶洞。裂缝和溶洞镶嵌在岩块系统中, 并且相互连成网络。流体在岩块系统中的流动采用三重孔隙介质模型进行刻画; 流体在裂缝中的流动规律采用达西渗流模型表征; 流体在溶洞中的流动为拟稳态流, 采用等势体模型进行表征。这类将溶洞视为等势体的方法的优点是模型构建和求解方便, 计算速度快, 可以计算出缝洞的体积等。缺点是无法确定缝洞的几何尺寸, 对长条状或条带型流动系统的适应性差。Popov等^[27]认为溶洞往往伴随着不同程度的充填, 提出应用Stokes-Brinkman方程来^[28]描述缝洞介质中的流动。在模型求解过程中, 通过选择适当的参数使得Stokes-Brinkman方程可以简化为Stokes方程或者达西渗流方程。同时, 由于采用了统一的控制方程, 避免了在自由流动区和渗流区交界面上的耦合计算, 使得问题求解简单化。Jie等^[29]认为Stokes-Brinkman方程既适用于稳态流也适用于非稳态流, 因此基于广义质量守恒方程和Stokes-Brinkman方程构建了管流-渗流耦合的不稳定流动模拟分析方法。黄朝琴等^[30]基于Stokes-Brinkman方程建立了离散缝洞网络宏观流动数学模型。

本文引入“ 广义流度” ^[31], 以实现管流和渗流在形式上的统一、实现常用线性和非线性流动规律在形式上的统一。油气藏在不同区域或不同尺度上能够使用相同形式的运动方程构建统一的控制方程, 进而能够将线性和非线性以及非线性和非线性的复杂耦合流动问题换化为复合流动问题。在广义流度的定义下, 层流管流和达西渗流在形式上没有任何差异。在此基础上, 提出基于广义流度的模型体系, 建立广义管流-渗流耦合储集体的基本控制方程。构建两个试井分析示例模型, 并采用Laplace变换方法对模型进行解析求解, 给出不同模型参数变化对应的典型特征图。结合具体实例, 验证本文方法的正确性和合理性。

图1给出了一种典型的复杂储集体管流-渗流耦合物理模型示意图。整个储集空间Ω 包含4个子空间：井筒空间Ω ₀、多孔介质空间Ω ₁、孔道空间Ω ₂及溶洞空间Ω ₃。流体在4个子空间中的流动分别符合不同的流动规律。

图1

	Figure Option ViewDownloadNew Window
	图1 管流-渗流耦合物理模型示意图

为了解决图1所示的管流-渗流耦合问题, 在储集体Ω 上定义一个形式上完全统一的运动方程：

$v=-\lambda \nabla p$ （1）

其中

$\lambda =\text{diag}\left[ {{\lambda }_{X}}(x, t), {{\lambda }_{Y}}(x, t), {{\lambda }_{Z}}(x, t) \right]$ $x\in \Omega $ $t\in \left( 0, T \right]$

（1）式中, λ 称为广义流度, 其为对角矩阵。λ _X(x, t), λ _Y(x, t), λ _Z(x, t)分别为X, Y, Z方向的广义流度分量。广义流度分量在子空间Ω ₀, Ω ₁, Ω ₂, Ω ₃内部是关于x的连续函数, 但并不要求在子空间边界$\partial {{\Omega }_{0}}$、$\partial {{\Omega }_{2}}$以及$\partial {{\Omega }_{3}}$上也是关于x的连续函数。

基于（1）式可构建统一的管流-渗流耦合不稳定流动控制方程：

$\operatorname{div}\left[ \rho \left( -\lambda \nabla p \right) \right]+\frac{\partial \left( \rho \phi \right)}{\partial t}=0$ （2）

初始条件为：

$p\left( x, 0 \right)={{h}_{0}}\left( x \right)$ $x\in \Omega $ （3）

边界条件为：

${{c}_{1}}p+{{c}_{2}}\lambda \frac{\partial p}{\partial n}={{h}_{1}}\left( x, t \right)$$x\in \partial {{\Omega }_{1, \text{out}}}$$t\in \left( 0, T \right]$ （4）

（2）— （4）式组成的方程组一般需要采用数值方法进行求解：先将整个储集空间Ω 剖分为一系列互不重叠的四面体, 然后利用有限体积法将（2）— （4）式组成的方程组进行离散化形成离散方程组, 最后求解离散方程组可获得求解区域内任意时刻压力分布。

（1）式同时描述了管流和渗流规律, 对一维流动, 根据不同的流体流动情况, 广义流度可以定义为表1的不同形式。为了便于理解和应用, 表1中列出了一部分常见单相流动规律对应的一维广义流度公式。这些一维广义流度可直接推广至二维和三维, 也可以推广到多相流。图1管流-渗流耦合物理模型中的不同部位可以使用表1中的不同公式形成复合模型。

表1

表1

表1 一维广义流度公式 流动规律 广义流度（依据已有研究成果推导得到） 线性 流动 达西渗流[4] $\lambda \text{=}K/\mu $ 层流管流[11] $\lambda ={{d}^{2}}/(32\mu )$ 非线性 流动 低速非达西流[5] $\lambda \text{=(}K/\mu )(1-G/|\partial p/\partial x|)$ 高速非达西流[6] $\lambda \text{=1/}\left( \mu /K+\beta \rho \left| v \right| \right)$ 幂律非牛顿流[7] $\lambda \text{= }\!\!|\!\!\text{ }v{{\text{ }\!\!|\!\!\text{ }}^{1-n}}K/{{\mu }_{\text{e}}}$ 应力敏感非达西流[8] $\lambda \text{=}\left( K/\mu \right){{\text{e}}^{\gamma \left( p-{{p}_{\text{i}}} \right)}}$ 紊流管流[9] $\lambda \text{=}2{{d}^{2}}\text{/}\left[ f(Re, e/d)\mu Re \right]$ 注：管流中的水力直径d用$\sqrt{32K}$替换, 以便于将一维管流推广至二维和三维, 也便于引入多相流的相渗透率

表1 一维广义流度公式

利用广义流度构建的复杂管流-渗流耦合模型在不同的区域具有完全相同的形式, 因而采用业界熟知的数值计算方法能够更方便地对模型进行离散化, 能够降低耦合问题的复杂性, 使问题的求解变得简单统一。

对于不稳定流动问题, 为了能够获得解析解, 下文将基于上述管流-渗流耦合模型, 对油藏中微可压缩单相流体符合线性流动规律的不稳定流动情况建立两种典型示例模型。

在存在层流管流和渗流耦合的一维不稳定流的管道形复合储集体（见图2）中, 井筒以定产量生产或者注入, 且开井前储集体中压力均为原始地层压力。同时, 假设储集体及流体满足：①储集体由两个管流-渗流耦合的管道形储集体复合而成; ②储集体1的左端与井筒相连, 右端与储集体2相连, 两个储集体中的流体流动均符合线性流动规律; ③岩石微可压缩; ④储集体1和2中流体为微可压缩单相流（使用拟压力函数将可压缩单相流及多相流问题转化为微可压缩单相流问题^{[32, 33]}）; ⑤考虑井筒储集效应及表皮效应; ⑥管道考虑为一维流动, 不考虑具体的井型, 储集体1左端的流出（入）流量等于井筒产出（注入）量。

图2

	Figure Option ViewDownloadNew Window
	图2 管道形复合储集体物理模型示意图

根据以上假设, 建立如下无因次模型。

$\frac{{{\partial }^{2}}{{p}_{\text{1D}}}}{\partial x_{\text{D}}^{2}}=\frac{\partial {{p}_{\text{1D}}}}{\partial {{t}_{\text{D}}}}$ （5）

其中 ${{t}_{\text{D}}}=\frac{{{\lambda }_{1}}{{t}_{\text{a}}}\left( t \right)}{{{\phi }_{1}}{{C}_{\text{t1}}}{{L}^{2}}}$ ${{x}_{\text{D}}}=\frac{x}{L}$

${{p}_{\text{1D}}}=\frac{{{\lambda }_{1}}{{A}_{1}}}{\text{2 }\!\!\pi\!\!\text{ }qL}\left[ M\left( {{p}_{\text{i}}} \right)-M\left( {{p}_{1}} \right) \right]$

$\frac{{{\partial }^{2}}{{p}_{\text{2D}}}}{\partial x_{\text{D}}^{2}}=\delta \frac{\partial {{p}_{\text{2D}}}}{\partial {{t}_{\text{D}}}}$ （6）

其中 $\delta =\frac{{{\lambda }_{1}}/\left( {{\phi }_{1}}{{C}_{\text{t1}}} \right)}{{{\lambda }_{2}}/\left( {{\phi }_{2}}{{C}_{\text{t2}}} \right)}$ ${{p}_{\text{2D}}}=\frac{{{\lambda }_{1}}{{A}_{1}}}{\text{2 }\!\!\pi\!\!\text{ }qL}\left[ M\left( {{p}_{\text{i}}} \right)-M\left( {{p}_{2}} \right) \right]$

${{p}_{\text{1D}}}\left( {{x}_{\text{D}}}, {{t}_{\text{D}}}=0 \right)={{p}_{\text{2D}}}\left( {{x}_{\text{D}}}, {{t}_{\text{D}}}=0 \right)=0$ （7）

$\frac{\partial p_{1D}}{\partial x_{D}}|_{x_{D}=x_{1D}} =m\frac{\partial p_{2D}}{\partial x_{D}}|_{x_{D}=x_{1D}} $ （8）

其中 ${{x}_{\text{1D}}}=\frac{{{x}_{1}}}{L}$ $m=\frac{{{A}_{2}}{{\lambda }_{2}}}{{{A}_{1}}{{\lambda }_{1}}}$

${{p}_{\text{1D}}}|_{{x}_{\text{D}}={x}_{\text{1D}}} ={{p}_{\text{2D}}}|_{{x}_{\text{D}}={x}_{\text{1D}}}$ （9）

${{p}_{\text{2D}}}\left| \begin{align} & \\ & {{x}_{\text{D}}}=\infty \\ \end{align} \right.=0$ （10）

$\frac{\partial {{p}_{\text{1D}}}}{\partial {{x}_{\text{D}}}}\left| \begin{align} & \\ & {{x}_{\text{D}}}=0 \\ \end{align} \right.=-1$ （11）

${{p}_{\text{wD}}}={{p}_{\text{1D}}}\left( {{x}_{\text{D}}}=0, {{t}_{\text{D}}} \right)$ （12）

基于t_D对（5）— （12）式作Laplace变换, 进行求解可得Laplace空间无因次井底压力解：

${{\bar{p}}_{\text{wD}}}=\frac{m\sqrt{\delta }\tanh \left( \sqrt{u}{{x}_{\text{1D}}} \right)+1}{{{u}^{3/2}}\left[ \tanh \left( \sqrt{u}{{x}_{1D}} \right)+m\sqrt{\delta } \right]}$ （13）

其中 ${{p}_{\text{wD}}}=\frac{{{\lambda }_{1}}{{A}_{1}}}{\text{2 }\!\!\pi\!\!\text{ }qL}\left[ M\left( {{p}_{\text{i}}} \right)-M\left( {{p}_{\text{w}}} \right) \right]$

对（13）式利用杜哈美原理^[4]可得Laplace空间考虑井筒储集效应和表皮效应的无因次井底压力解：

${{\bar{p}}_{\text{wD}}}=\frac{u{{{\bar{p}}}_{\text{wD}}}+S}{u+{{u}^{2}}{{C}_{\text{D}}}\left( u{{{\bar{p}}}_{\text{wD}}}+S \right)}$ （14）

其中 ${{C}_{\text{D}}}=\frac{C}{{{\phi }_{1}}{{C}_{\text{t1}}}{{A}_{1}}L}$

采用Euler数值反演算法^[34]对（14）式进行反演获得实空间的井底压力, 并绘制部分参数变化的压力降落曲线（见图3、图4）和压力恢复曲线（见图5）。压力恢复曲线的绘制按照如下计算公式实施^[35]：

图3

	Figure Option ViewDownloadNew Window
	图3 组合参数mδ 1/2变化响应特征（压力降落）图

图4

	Figure Option ViewDownloadNew Window
	图4 储集体1长度变化响应特征（压力降落）图

图5

	Figure Option ViewDownloadNew Window
	图5 生产时间对压力恢复曲线的影响

${{p}_{\text{sD}}}={{p}_{\text{wD}}}\left( {{t}_{\text{D}}} \right)+{{p}_{\text{wD}}}\left( {{t}_{\text{pD}}} \right)-{{p}_{\text{wD}}}\left( {{t}_{\text{pD}}}+{{t}_{\text{D}}} \right)$ （15）

其中 ${{p}_{\text{sD}}}=\frac{{{\lambda }_{1}}{{A}_{1}}}{\text{2 }\!\!\pi\!\!\text{ }qL}\left[ M\left( {{p}_{\text{i}}} \right)-M\left( {{p}_{\text{s}}} \right) \right]$

图3为组合参数mδ ^1/2变化响应特征（压力降落）图。大体上可分为4个流动段, Ⅰ 为纯井储段及早期过渡段, 压力和压力导数曲线前期基本呈斜率为1的直线, 之后出现驼峰（正表皮作用）; Ⅱ 、Ⅳ 主要为储集体1、2的线性流段, 后期压力和压力导数曲线基本呈斜率为0.5的平行线; Ⅲ 为储集体1、2线性流之间的过渡段。随着mδ ^1/2的增大, 压力和压力导数曲线位置下移; 当mδ ^1/2=1时, 两个线性流段之间无过渡段响应特征。

图4为储集体1长度变化响应（压力降落）特征图。显然, 储集体1越长, 其线性流响应时间越长。

图5展示了生产时间对压力恢复曲线的影响。早期井筒储集段压力恢复导数与压力降落导数曲线重合, 二者无明显差异。随着生产时间的延长, 压力恢复导数曲线可以无限地接近压力降落导数曲线。生产时间越短, 压力恢复导数曲线开始偏离压力降落导数曲线的时间就越早, 与压力降落曲线偏差也越大。生产时间一定时, 只要测试时间足够长, 压力恢复导数曲线后期是下坠的。

在存在层流管流和渗流耦合的三维不稳定流的圆柱形储集体（见图6）中, 生产井以定产量生产或者注入, 且开井前储集体中压力均为原始地层压力, 同时假设储集体及流体满足：①储集体为单个管流与渗流耦合的圆柱形储集体, 储集体垂向和径向的广义流度具有差异性, 储集体水平等厚, 储集体中流体流动符合线性流动规律; ②岩石微可压缩; ③储集体中流体为微可压缩单相流（使用拟压力函数将可压缩单相流及多相流问题转化为微可压缩单相流问题^{[32, 33]}）; ④井筒位于圆柱形储集体的中轴线上, 且只打开一部分; ⑤考虑表皮效应及井筒储集效应。

图6

	Figure Option ViewDownloadNew Window
	图6 圆柱形储集体物理模型示意图

根据模型假设, 在柱坐标中可建立如下无因次模型：

$\frac{1}{{{r}_{\text{D}}}}\frac{\partial }{\partial {{r}_{\text{D}}}}\left( {{r}_{\text{D}}}\frac{\partial {{p}_{\text{D}}}}{\partial {{r}_{\text{D}}}} \right)+\frac{{{\lambda }_{\text{z}}}}{{{\lambda }_{\text{r}}}}\frac{{{\partial }^{2}}{{p}_{\text{D}}}}{\partial z_{\text{D}}^{2}}=\frac{\partial {{p}_{\text{D}}}}{\partial {{t}_{\text{D}}}}$ （16）

其中 ${{t}_{\text{D}}}=\frac{{{\lambda }_{\text{r}}}{{t}_{\text{a}}}\left( t \right)}{\phi {{C}_{\text{t}}}{{L}^{2}}}$ ${{p}_{\text{D}}}=\frac{{{\lambda }_{\text{r}}}h}{q}\left[ M\left( {{p}_{\text{i}}} \right)-M\left( p \right) \right]$

${{p}_{\text{D}}}({{r}_{\text{D}}}, {{z}_{\text{D}}}, {{t}_{\text{D}}}=0)={{p}_{\text{wD}}}=0$ （17）

其中 ${{r}_{\text{D}}}=\frac{r}{L}$ ${{z}_{\text{D}}}=\frac{z}{L}$ ${{p}_{\text{wD}}}=\frac{{{\lambda }_{\text{r}}}h}{q}\left[ M\left( {{p}_{\text{i}}} \right)-M\left( {{p}_{\text{w}}} \right) \right]$

$\frac{{{\varepsilon }_{\text{D}}}}{{{h}_{\text{D}}}}{{r}_{\text{D}}}\frac{\partial {{p}_{\text{D}}}}{\partial {{r}_{\text{D}}}}|_{r_{D}=1}=\left\{ \begin{align} & -1\quad \quad \left| {{z}_{\text{D}}}-{{z}_{\text{wD}}} \right|\le \frac{{{\varepsilon }_{\text{D}}}}{2} \\ & 0\quad \ \ \quad \left| {{z}_{\text{D}}}-{{z}_{\text{wD}}} \right|\frac{{{\varepsilon }_{\text{D}}}}{2} \\ \end{align} \right.$ （18）

其中 ${{z}_{\text{wD}}}=\frac{{{z}_{\text{w}}}}{L}$ ${{h}_{\text{D}}}=\frac{h}{L}$ ${{\varepsilon }_{\text{D}}}=\frac{\varepsilon }{L}$

${{p}_{\text{wD}}}=\frac{1}{{{\varepsilon }_{\text{D}}}}\int\limits_{{{z}_{\text{wD}}}-{{\varepsilon }_{\text{D}}}/2}^{{{z}_{\text{wD}}}+{{\varepsilon }_{\text{D}}}/2}{{{p}_{\text{D}}}({{r}_{\text{D}}}\text{=}1, {{z}_{\text{D}}}, {{t}_{\text{D}}})\operatorname{d}{{z}_{\text{D}}}}$ （19）

$\frac{\partial {{p}_{\text{D}}}}{\partial {{z}_{\text{D}}}}\left| \begin{align} & \\ & {{z}_{\text{D}}}=0 \\ \end{align} \right.=0$ （20）

$\frac{\partial {{p}_{\text{D}}}}{\partial {{z}_{\text{D}}}}\left| \begin{align} & \\ & {{z}_{\text{D}}}={{h}_{\text{D}}} \\ \end{align} \right.=0$ （21）

$\frac{\partial {{p}_{\text{D}}}}{\partial {{r}_{\text{D}}}}\left| \begin{align} & \\ & {{r}_{\text{D}}}={{r}_{\text{eD}}} \\ \end{align} \right.=0$ （22）

其中 ${{r}_{\text{eD}}}=\frac{{{r}_{\text{e}}}}{L}$

基于t_D对（16）— （22）式作Laplace变换, 并利用分离变量法求得Laplace空间无因次井底压力解：

${{\bar{p}}_{\text{wD}}}=\frac{1}{u}\sum\limits_{i=0}^{\infty }{{{a}_{i}}\frac{{{\text{K}}_{0}}\left( {{b}_{i}} \right)+\frac{{{\text{K}}_{1}}\left( {{r}_{\text{eD}}}{{b}_{i}} \right)}{{{\text{I}}_{1}}\left( {{r}_{\text{eD}}}{{b}_{i}} \right)}{{\text{I}}_{0}}\left( {{b}_{i}} \right)}{{{b}_{i}}{{\text{K}}_{1}}\left( {{b}_{i}} \right)-\frac{{{\text{K}}_{1}}\left( {{r}_{\text{eD}}}{{b}_{i}} \right)}{{{\text{I}}_{1}}\left( {{r}_{\text{eD}}}{{b}_{i}} \right)}{{b}_{i}}{{\text{I}}_{1}}\left( {{b}_{i}} \right)}}$ （23）

其中 ${{a}_{i}}=\left\{ \begin{array}{* {35}{l}} 1 & \quad i=0 \\ 8{{\left( \frac{\sin \frac{i\text{ }\!\!\pi\!\!\text{ }{{\varepsilon }_{\text{D}}}}{2{{h}_{\text{D}}}}\cos \frac{i\text{ }\!\!\pi\!\!\text{ }{{z}_{\text{wD}}}}{{{h}_{\text{D}}}}}{\frac{i\text{ }\!\!\pi\!\!\text{ }{{\varepsilon }_{\text{D}}}}{{{h}_{\text{D}}}}} \right)}^{2}} & \quad i=1, 2, \cdots , \infty \\\end{array} \right.$

${{b}_{i}}=\sqrt{u+{{\lambda }_{\text{z}}}/{{\lambda }_{\text{r}}}{{\left( i\text{ }\!\!\pi\!\!\text{ }/{{h}_{\text{D}}} \right)}^{2}}}$

对（23）式利用杜哈美原理^[4]可得Laplace空间考虑井筒储集效应和表皮效应的无因次井底压力解：

${{\bar{p}}_{\text{wD}}}=\frac{u{{{\bar{p}}}_{\text{wD}}}+S}{u+{{u}^{2}}{{C}_{\text{D}}}\left( u{{{\bar{p}}}_{\text{wD}}}+S \right)}$ （24）

其中 ${{C}_{\text{D}}}=\frac{C}{2\text{ }\!\!\pi\!\!\text{ }\phi {{C}_{\text{t}}}h{{L}^{2}}}$

采用Euler数值反演算法^[34]对（24）式进行反演获得实空间的井底压力, 并绘制部分参数变化的压力降落图（见图7— 图10）和压力恢复图（见图11— 图12）。压力恢复图的绘制按照（15）式实施, 此时：

${{p}_{\text{sD}}}=\frac{{{\lambda }_{\text{r}}}h}{q}\left[ M\left( {{p}_{\text{i}}} \right)-M\left( {{p}_{\text{s}}} \right) \right]$ （25）

图7为圆柱形储集体半径变化响应特征（压力降落）图。在垂向广义流度与径向广义流度相等, 打开段厚度相对储集体厚度较小, 储集体厚度一定条件下, 随着储集体半径的增大压力导数曲线后期拟稳定流出现时间延迟。除此之外, 当r_eD小于h_D时, 压力导数曲线中部可能出现球形流和线性流两个特征（图7红色曲线）; 随着r_eD的增大压力导数曲线中部球型流特征越来越明显, 而线性流特征逐渐消失（图7黑色曲线）; 当r_eD继续增大时, 压力导数曲线可能出现球形流和径向流两个特征（图7蓝色曲线）。

图7

	Figure Option ViewDownloadNew Window
	图7 储集体半径变化特征（压力降落）图

图8为储集体厚度变化响应特征（压力降落）图。在垂向广义流度与径向广义流度相等, 打开段厚度相对储集体厚度较小且二者比值一定, 储集体半径大小一定条件下, 随着储集体厚度的增大压力导数曲线中部下凹变浅。

图8

	Figure Option ViewDownloadNew Window
	图8 储集体厚度变化特征（压力降落）图

图9为垂向与径向广义流度比变化响应特征（压力降落）图。在储集体厚度、储集体半径及打开段厚度一定条件下, 随着广义流度比的增大, 压力导数曲线中部下凹变深, 线性流持续时间变短。

图9

	Figure Option ViewDownloadNew Window
	图9 垂向与径向广义流度比变化特征（压力降落）图

图10为打开段比例变化响应特征（压力降落）图。在储集体厚度、储集体半径及垂向与径向广义流度比一定条件下, 随着打开段厚度的增大, 压力导数曲线过渡段位置变低, 过渡段可能出现环绕打开段的早期径向流（一般较难出现）。一种特殊情况是当ε _D=h_D时, 压力导数曲线中部无球形流和线性流特征响应。

图10

	Figure Option ViewDownloadNew Window
	图10 打开段比例变化特征（压力降落）图

图11和图12为生产时间对压力恢复曲线的影响, 其中图11为垂向流度相对较大情况, 图12为垂向广义流度相对较小情况。早期井筒储集段压力恢复导数与压力降落导数曲线重合, 二者无明显差异。随着生产时间的增大, 压力恢复导数曲线可以更加接近压力降落导数曲线。但当生产时间增大到一定值时, 压力恢复导数曲线的形状基本不再发生变化。生产时间越小, 压力恢复导数曲线开始偏离压力降落导数曲线的时间就越早, 与压力降落导数曲线偏差也越大。生产时间一定时, 只要测试时间足够长, 压力恢复导数曲线后期是下坠的。

图11

	Figure Option ViewDownloadNew Window
	图11 垂向广义流度相对较大时生产时间对压力恢复曲线的影响

图12

	Figure Option ViewDownloadNew Window
	图12 垂向广义流度相对较小时生产时间对压力恢复曲线的影响

通过定义广义流度将不同的流体运动方程进行了统一。给出了单相流管渗耦合的一般试井分析模型, 该模型能够用于复杂管渗耦合流动模拟。

在广义流度意义下, 整个油气藏在不同区域或不同尺度上均能使用相同形式的运动方程构建统一的控制方程。因此, 可以更方便地对模型进行离散化求解, 降低耦合问题的复杂性, 使模型求解变得简单统一。

构建了两种典型管流-渗流耦合的试井分析示例模型。管道形复合储集体模型压力降落图可出现两个线性流特征。圆柱形储集体模型压力降落图压力导数曲线中部可出现球形流和线性流两个特征或可出现球形流和径向流两个特征。通过应用实例分析验证了模型的实用性和可靠性。

目前给出的两个示例模型为解析解, 为了获得更丰富油气藏参数, 下一步将重点研究广义流度为非线性情况的一般性数值求解方法, 实现对油气藏精细特征更准确的解释和描述。

符号注释：

A₁, A₂— — 储集体1和2过流面积, m²; c₁— — 压力函数项系数, Pa^-1; c₂— — 压力梯度函数项系数, s/m; C— — 井筒储集系数, m³/Pa; C_t— — 综合压缩系数, Pa^-1; C_t1, C_t2— — 储集体1和2综合压缩系数, Pa^-1; d— — 水力直径, m; e— — 粗糙度, m; f(· )— — 关于Re和e/d的函数; G— — 启动压力梯度, Pa/m; h— — 圆柱形储集体高度, m; h₀— — 初始压力分布函数, Pa; h₁— — 边界条件函数, 无因次; I₀(· ), I₁(· ) — — 0阶和1阶第1类修正Bessel函数; K— — 渗透率, m²; K₀(· ), K₁(· )— — 0阶和1阶第2类修正Bessel函数; L— — 参考长度（可任选）, m; m— — 过流面积比与广义流度比的乘积; M(· )— — 拟压力函数, Pa; n— — 流动特性指数; n— — 储集体边界的法向量, m; p— — 压力, Pa; p_i— — 初始压力, Pa; p_s— — 关井恢复时的井底压力, Pa; p_w— — 开井生产时的井底压力, Pa; ${{\bar{p}}_{\text{w}}}$— — Laplace空间井底压力, Pa; p₁, p₂— — 储集体1和储集体2沿水平方向的压力分布, Pa; q— — 流量, m³/s; r— — 径向距离, m; r_e, r_w— — 圆柱形储集体半径和井筒半径, m; Re— — 雷诺数; S— — 井筒与储集体连接处的表皮因子, 在管道形复合储集体中以储集体1的过流面积为参考, 在圆柱形储集体中以井筒打开段面积为参考; t— — 时间, s; t_a(· )— — 拟时间函数, s; t_pD— — 无因次生产时间, s; T— — 时间最大值, s; u— — 无因次时间对应的Laplace变量; v— — 速度, m/s; v— — 速度矢量, m/s; x— — 距离, m; x— — 距离矢量, m; x₁— — 储集体1的长度, m; X, Y, Z— — 直角坐标系, m; z— — 垂向距离, m; z_w— — 打开段中心点垂向位置, m; β — — 高速非达西系数, m^-1; γ — — 渗透率模量, Pa^-1; δ — — 导压系数比; ε — — 打开段厚度, m; λ — — 广义流度, m²/(Pa· s); λ — — 广义流度矩阵, m²/(Pa· s); λ _r, λ _z— — 储集体径向和垂向广义流度, m²/(Pa· s); λ _X(x, t), λ _Y(x, t), λ _Z(x, t)— — X, Y, Z方向的广义流度分量, m²/(Pa· s); λ ₁, λ ₂— — 储集体1和2的广义流度, m²/(Pa· s); μ — — 流体黏度, Pa· s; μ _e— — 有效黏度, (Pa· s)· (m/s)^1-n; ρ — — 密度, kg/m³; ϕ — — 孔隙度, %; ϕ ₁, ϕ ₂— — 储集体1和2的孔隙度, %。下标：D— — 无因次。

（编辑 胡苇玮）