分形复合非牛顿幂律流体渗流模型解的研究

Research on the Solution of Fractal Composite Non-Newtonian Power-Law Fluid Seepage Model

摘要: 在实际油气田开发中，分形复合油藏由于储层孔隙结构复杂、非均质性较强，地层流体多表现为非牛顿幂律流体特性，针对此问题，本文建立了分形复合非牛顿幂律流体渗流模型，该模型结合了分形油藏的地质特征和非牛顿流体的流变特征，考虑了井筒储集与表皮因子的内边界情形，以便更精确地贴合实际油藏边界的弹性变形特性。在求解过程中，本文推广了相似构造法，结合无量纲化、Laplace变换及数值反演，获得了统一相似连分数结构的实空间解，从而简化了求解过程。

Abstract: In practical oil and gas field development, fractured composite reservoirs exhibit complex pore structures and strong heterogeneity, with formation fluids often displaying non-Newtonian power-law fluid characteristics. To address this problem, this paper establishes a seepage model for fractal composite non-Newtonian power-law fluids. This model combines the geological characteristics of fractal reservoirs and the rheological properties of non-Newtonian fluids, taking into account the inner boundary conditions of wellbore storage and skin factor, so as to more accurately reflect the elastic deformation characteristics of actual reservoir boundaries. In the solution process, this study extends the similarity construction method, combining dimensionless analysis, Laplace transform, and numerical inversion to obtain a unified similarity continued fraction structure for the real space solution, thus simplifying the solution process.

1. 引言

石油作为不可再生重要矿产资源，其油藏合理开发与高效利用对经济社会发展意义重大。随着易开采油藏枯竭，复杂地质结构与流体性质给传统渗流理论带来挑战，分形理论可有效描述地层非均质特征，复合油藏能表征多介质油藏系统，但目前相关研究多集中于牛顿流体，对非牛顿幂律流体的研究较少，而实际油藏中常存在牛顿与非牛顿渗流区混杂现象，因此建立分形复合非牛顿幂律渗流模型具有重要实际意义。相较于封闭、定压、无穷大等理想外边界，弹性外边界条件更具一般性，可将三种理想边界作为特殊情形，基于此，本文引入弹性外边界条件，建立考虑井筒储集和表皮效应的分形复合油藏非牛顿幂律流体模型，以探究油藏渗流规律、拓展渗流力学理论，为油藏开发方案优化、采油效率提升提供理论支撑与计算工具。

1961年Mandelbrot [1]建立分形理论基础框架，该理论随后渗透到多学科领域，成为重要跨学科工具。1990年Katz和Lee [2]建立分形流动模型，提出其可有效描述复杂油藏流体流动的非均质性。2003年Flamenco-Lopez和Camacho-Velázquez [3]提出基于压力瞬态数据反演裂缝网络分形参数的方法，为裂缝油藏流动建模提供理论支持。2004年Jordan [4]通过理论与数值模拟，对比复合与双层油藏压力瞬态响应，为油藏类型判别提供参考。2001~2004年间，向开理[5] [6]等人结合分形几何与渗流力学，构建三种理想外边界下分形复合油藏不稳定渗流模型并求得解析解，推动了相关试井分析技术发展。2006年徐昌学、李顺初[7]等人引入相似结构理论，实现分形复合油藏试井解析解的规范化表达，简化求解分析过程。2007年孔祥言[8]等人推导分形渗流基本方程，建立压力动态求解方法并绘制试井样板曲线，奠定相关研究理论基础。2013年Santos [9]通过实验与数值模拟，发现油藏分形维数对气驱、水驱效果的影响，为开发方案调整提供依据。2014年李顺初[10]等人提出复合油藏球向渗流的相似结构解法，提供了新的求解思路。2017~2019年Abbasi [11] [12]等人求解非线性扩散方程解析解，建立考虑多因素的裂缝性储层瞬态形状因子模型，提升计算贴合度。2020年Su [13]将分形理论应用于低渗透油藏油水两相流动分析，实现油水相对渗透率定量刻画。2023年罗静[14]引入弹性概念建立渗流模型，构建考虑多因素的分形复合油藏模型。2024年Dong [15]等人提出基于分形几何的复合油藏建模方法，更贴合实际油藏特征。目前分形模型应用仍有瓶颈，分形参数提取及油藏边界复杂性刻画仍是后续研究重点。

非牛顿流体指黏度不恒定，且依赖于剪切速率、温度、应力等因素的流体，与牛顿流体相比，其流动性质复杂多变，通常表现出剪切变稀、剪切增稠、时变性等特性。1969年Van Poollen [16]等人率先开展非牛顿幂律流体在多孔介质内瞬态流动行为的系统性研究，通过理论推导得出相关偏微分方程，为后续研究搭建基础理论框架；1979年Ikoku [17]基于稳态流量假设，给出非牛顿幂律流体的渗流控制方程，并研究其在油藏中的不稳定渗流特性，推动其在油藏工程中的应用。1980年Ikoku [18]扩展先前研究，重点考虑井眼储存效应对非牛顿流体渗流行为的影响，揭示其对压力瞬态的作用机制；1987年Vongvuthipornchai [19]针对外边界无穷大的均质储层，构建非牛顿幂律流体试井解释模型，并提出通过压力导数方法改进试井数据分析的方案。1996年刘泽俊等[20]研究有界地层中的非牛顿幂律流体试井模型，给出不同边界条件下的解析解；1997年宋考平等[21]提出考虑多个非牛顿渗流区的复合油藏模型，对聚合物驱等非牛顿流体驱油的试井分析具有重要意义。2000年徐建平[22]提出适用于球向渗流的剪切速率计算公式，构建考虑多因素的非牛顿幂律流体不稳定球形流动模型；2009年王雅茹等[23]采用Duhamel褶积法，构建非牛顿幂律流体现代试井典型理论图版，为试井解释提供直观参考。2016年龚伟等[24]总结非牛顿流体试井分析的研究现状，展望未来研究方向并强调过渡区附加阻力的影响；2019年刘文涛等[25]建立三区复合油藏试井解释模型，绘制聚驱过程渗流特征双对数理论曲线并识别出四个流动阶段。2022年Oluogun等[26]推导二维解析解，构建通用型复合数学模型，进一步完善非牛顿流体在复合油藏中的渗流理论体系。总体而言，非牛顿流体在试井分析中的应用逐步深入，未来将继续聚焦相关影响因素，提升试井解释的准确性与实用性。

2. 预备知识

变型Bessel函数的定义与性质

变型bessel函数作为一类特殊函数的统称，其本质是满足如下二阶常微分方程的解[27]，

$x^2\frac{{\text{d}}^{2}y}{\text{d}{x}^2}+x\frac{\text{d}y}{\text{d}x}+\left(x-v^2\right)y=0,$ (2.1)

变型bessel函数可分为第一类变型bessel函数$I_v\left(x\right)$ 和第二类变型bessel函数$K_v\left(x\right)$ ，具体定义如下：

$I_v\left(x\right)={\displaystyle \sum _{m=0}^{\infty }\frac{1}{m!\Gamma \left(m+v+1\right)}}{\left(\frac{x}{2}\right)}^{2m+v},$ (2.2)

其中，$\Gamma (.)$ 表示Gamma函数。

$K_v\left(x\right)=\underset{\mu \to v}{\lim }\frac{J_{\mu }\left(x\right)\cos \left(\mu \pi \right)-J_{-\mu }\left(x\right)}{\sin \left(\mu \pi \right)},$ (2.3)

在渗流模型的分析与求解过程中，需用到以下变型Bessel函数的组合

${\psi }_{m,n}\left(x,y,z\right)=K_m\left(xz\right)I_n\left(yz\right)+{\left(-1\right)}^{m+n-1}{I}_m\left(xz\right)K_n\left(yz\right),$

其中，$m$ ，$n$ 为常实数，$x$ ，$\chi$ ，$\delta$ 为实变量，$I_m(\cdot )$ 和$K_m(\cdot )$ 分别为第一类、第二类$m$ 阶变型bessel函数，$I_n(\cdot )$ 和$K_n(\cdot )$ 分别为第一类、第二类$n$ 阶变型bessel函数，为了方便研究渗流模型的解，给出如下变型Bessel函数的性质：

1) 微分性质

$\{\begin{array}{l}\frac{\text{d}}{\text{d}x}{I}_v\left(x\right)=\frac{v}{x}{I}_v\left(x\right)+I_{v+1}\left(x\right),\\ \frac{\text{d}}{\text{d}x}{K}_v\left(x\right)=\frac{v}{x}{K}_v\left(x\right)-K_{v+1}\left(x\right).\end{array}$

2) ${\psi }_{m,n}\left(x,\chi ,\delta \right)$ 的微分性质

$\frac{\partial }{\partial x}{\psi }_{m,n}\left(x^{\theta },y^{\theta },z\right)=\frac{\theta m}{x}{\psi }_{m,n}\left(x^{\theta },y^{\theta },z\right)-\theta z{x}^{\theta -1}{\psi }_{m+1,n}\left(x^{\theta },y^{\theta },z\right)$

$\frac{\partial }{\partial y}{\psi }_{m,n}\left(x^{\theta },y^{\theta },z\right)=\frac{\theta n}{y}{\psi }_{m,n}\left(x^{\theta },y^{\theta },z\right)+\theta z{y}^{\theta -1}{\psi }_{m,n+1}\left(x^{\theta },y^{\theta },z\right)$

$\frac{\partial }{\partial x}{\psi }_{m,n+1}\left(x^{\theta },y^{\theta },z\right)=\frac{\theta m}{x}{\psi }_{m,n+1}\left(x^{\theta },y^{\theta },z\right)-\theta z{x}^{\theta -1}{\psi }_{m+1,n+1}\left(x^{\theta },y^{\theta },z\right)$

$\frac{\partial }{\partial y}{\psi }_{m+1,n}\left(x^{\theta },y^{\theta },z\right)=\frac{\theta n}{y}{\psi }_{m+1,n}\left(x^{\theta },y^{\theta },z\right)+\theta z{y}^{\theta -1}{\psi }_{m+1,n+1}\left(x^{\theta },y^{\theta },z\right)$

3) 渐近公式

当$x\to \infty$ 时，$I_v\left(x\right)$ 和$K_v\left(x\right)$ 有如下渐近公式

$\{\begin{array}{l}{I}_v\left(x\right)=\frac{{\text{e}}^x}{\sqrt{2\pi x}}\left[1+O\left(x^{-1}\right)\right],\\ K_v\left(x\right)=\sqrt{\frac{\pi }{2x}}{\text{e}}^{-x}\left[1+O\left(x^{-1}\right)\right].\end{array}$

通过该渐进公式可知，当$x\to \infty$ 时$I_v\left(x\right)$ 按${\text{e}}^x$ 型指数趋于无穷大，而$K_v\left(x\right)$ 按${\text{e}}^{-x}$ 型指数趋近于零。

3. 一类复合变型Bessel微分方程的边值问题的解

$\{\begin{array}{l}\frac{{\partial }^2{y}_1\left(x,z\right)}{\partial x^2}+\frac{A_1}{x}\frac{\partial y_1\left(x,z\right)}{\partial x}-B_1{x}^{C_1}{y}_1\left(x,z\right)=0a<x<b,\\ \frac{{\partial }^2{y}_2\left(x,z\right)}{\partial x^2}+\frac{A_2}{x}\frac{\partial y_2\left(x,z\right)}{\partial x}-B_2{x}^{C_2}{y}_2\left(x,z\right)=0b<x<c,\\ {\left[E{y}_1\left(x,z\right)+\left(1+EF\right)\frac{\partial y_1\left(x,z\right)}{\partial x}\right]|}_{x=a}=D,\\ y_1\left(x,z\right){|}_{x=b}-y_2\left(x,z\right){|}_{x=b}=0,\\ {\frac{\partial y_1\left(x,z\right)}{\partial x}|}_{x=b}-\lambda x^{\theta }{\frac{\partial y_2\left(x,z\right)}{\partial x}|}_{x=b}=0,\\ {\left[\left(N_1+N_2{x}^2\right)y_2\left(x,z\right)+N_3\frac{{\partial }^2{y}_2\left(x,z\right)}{\partial z^2}+x\frac{\partial y_2\left(x,z\right)}{\partial x}\right]|}_{x=c}=0.\end{array}$ (3.1)

其中，$y_i\left(x,z\right)$ 是关于$x,z$ 的二元函数，$A_i,B_i,C_i,D,E,F,a,b,c,\lambda ,\theta ,N_1,N_2,N_3\in R$ ，且$D\ne 0$ ，$a<b<c$ ，$C_i\ne 0,\left(i=1,2\right)$ 。

若上述边值问题方程存在唯一解，那么方程在内区$\left(a<x<b\right)$ 的解可以表示为

$y_1\left(x,z\right)=D\cdot \frac{1}{E+\frac{1}{F+{\Phi }_1\left(a\right)}}\cdot \frac{1}{F+{\Phi }_1\left(a\right)}\cdot {\Phi }_1\left(x\right)$ (3.2)

在外区$\left(b<x<c\right)$ 的解可以表示为

$y_2\left(x,z\right)=D\cdot \frac{1}{E+\frac{1}{F+{\Phi }_1\left(a\right)}}\cdot \frac{1}{F+{\Phi }_1\left(a\right)}\cdot \frac{{\phi }_{0,1}^1\left(b,b\right)}{{\Phi }_2\left(b\right){\phi }_{1,1}^1\left(a,b\right)-\lambda b^{\theta }{\phi }_{1,0}^1\left(a,b\right)}\cdot {\Phi }_2\left(x\right)$ (3.3)

式(3.2)、(3.3)中的${\Phi }_1\left(x\right)$ 称为内区的相似核函数，${\Phi }_2\left(x\right)$ 称为外区的相似核函数，可分别表示为

${\Phi }_1\left(x\right)=\frac{{\Phi }_2\left(b\right){\phi }_{0,1}^1\left(x,b\right)-\lambda b^{\theta }{\phi }_{0,0}^1\left(x,b\right)}{{\Phi }_2\left(b\right){\phi }_{1,1}^1\left(a,b\right)-\lambda b^{\theta }{\phi }_{1,0}^1\left(a,b\right)}$ (3.4)

${\Phi }_2\left(x\right)=\frac{\left(N_1+N_2{c}^2\right){\phi }_{0,0}^2\left(x,c\right)+N_3{\phi }_{0,1}^{*}\left(x,c\right)+c{\phi }_{0,1}^2\left(x,c\right)}{\left(N_1+N_2{c}^2\right){\phi }_{1,0}^2\left(b,c\right)+N_3{\phi }_{0,1}^{\#}\left(b,c\right)+c{\phi }_{1,1}^2\left(b,c\right)}$ (3.5)

其中${\phi }_{m,n}^i\left(x,\xi \right)\left(i=1,2;m,n=0,1\right)$ 被称为引解函数，定义如下

${\phi }_{0,0}^i\left(x,\xi \right)=y_{i1}\left(x,z\right)y_{i2}\left(\xi ,z\right)-y_{i2}\left(x,z\right)y_{i1}\left(\xi ,z\right)$ (3.6)

${\phi }_{1,0}^i\left(x,\xi \right)=\frac{\partial y_{i1}\left(x,z\right)}{\partial x}{y}_{i2}\left(\xi ,z\right)-\frac{\partial y_{i2}\left(x,z\right)}{\partial x}{y}_{i1}\left(\xi ,z\right)$ (3.7)

${\phi }_{0,1}^i\left(x,\xi \right)=y_{i1}\left(x,z\right)\frac{\partial y_{i2}\left(\xi ,z\right)}{\partial \xi }-y_{i2}\left(x,z\right)\frac{\partial y_{i1}\left(\xi ,z\right)}{\partial \xi }$ (3.8)

${\phi }_{1,1}^i\left(x,\xi \right)=\frac{\partial y_{i1}\left(x,\xi \right)}{\partial x}\frac{\partial y_{i2}\left(\xi ,z\right)}{\partial \xi }-\frac{\partial y_{i2}\left(x,z\right)}{\partial x}\frac{\partial y_{i1}\left(\xi ,z\right)}{\partial \xi }$ (3.9)

${\phi }_{0,1}^{*}\left(x,\xi \right)=y_{21}\left(x,z\right)\frac{{\partial }^2{y}_{22}\left(\xi ,z\right)}{\partial z^2}-y_{22}\left(x,z\right)\frac{{\partial }^2{y}_{21}\left(\xi ,z\right)}{\partial z^2}$ (3.10)

${\phi }_{1,1}^{\#}\left(x,\xi \right)=\frac{y_{21}\left(x,z\right)}{\partial x}\frac{{\partial }^2{y}_{22}\left(\xi ,z\right)}{\partial z^2}-\frac{\partial y_{22}\left(x,z\right)}{\partial x}\frac{{\partial }^2{y}_{21}\left(\xi ,z\right)}{\partial z^2}$ (3.11)

其中$y_1{}_1\left(x,z\right)$ ，$y_{12}\left(x,z\right)$ ，$y_2{}_1\left(x,z\right)$ ，$y_{22}\left(x,z\right)$ 为(3.1)中两个内外区方程的解。

4. 分形复合非牛顿幂律流体渗流模型的建立与求解

4.1. 模型的建立

在分形油藏中，由于经典达西定律中的渗透率与介质特征长度紧密相关，该定律无法直接应用于分形结构的渗流研究。为此，研究者将分形几何理论与达西定律结合，建立了适用于分形渗流网络的理论框架。此外，考虑到实际油藏中的非牛顿流体特性，模型引入了幂律流体流变特性。通过流动行为指数和稠度系数，能够描述流体的剪切变稀或剪切增稠现象，从而更准确地反映油藏中复合流体的渗流规律。即将分形维数为$d_f$ 的分形渗流网络嵌入到$d$ 维的欧几里得空间中，则可得到内、外区储集层的孔隙度和渗透率的分布函数，即：

${\Phi }_i\left(r\right)={\varphi }_i{\left(\frac{r}{r_0}\right)}^{d_{f_i}-d},$

$K_i\left(r\right)=k_i{\left(\frac{r}{r_0}\right)}^{d_{f_i}-d-{\theta }_i}\left(i=1,2\right)$

其中，${\Phi }_i\left(r\right)$ 表示孔隙度，$K_i\left(r\right)$ 表示渗透率；$r_0$ 表示参考半径，本文取井筒半径$r_w$ 为参考半径。内外区均为非牛顿幂律流体，其中$i=1$ 表示内区，$i=2$ 表示外区；${\varphi }_i$ 和$k_i$ 表示$r_0$ 处内外区孔隙度和渗透率；$d_{f_i}$ 是分形维数，${\theta }_i$ 是分形指数。

1) 连续性方程：

由质量守恒定理，可得流体的连续性方程为

$\frac{1}{r}\frac{\partial \left(r{\rho }_i{v}_i\right)}{\partial r}=-\frac{\partial \left({\varphi }_i{\rho }_i\right)}{\partial t}\left(i=1,2\right)$ (4.1)

2) 运动方程：

内外区为非牛顿幂律流体渗流区，且非牛顿流体为幂律型，则运动方程为

$v_i={\left(-\frac{k_i}{{\mu }_e}\frac{\partial p_i}{\partial r}\right)}^{\frac{1}{n}}\left(i=1,2\right)$ (4.2)

3) 状态方程：

${\rho }_i={\rho }_{0i}{\text{e}}^{C_{{\rho }_i}(p_i-p_0)}\left(i=1,2\right)$ (4.3)

${\varphi }_i={\varphi }_{i0}{\text{e}}^{C_{{\varphi }_i}(p_i-p_0)}\left(i=1,2\right)$ (4.4)

将等式(4.2)，(4.3)，(4.4)带入等式(4.1)得到内区渗流微分方程

$\frac{{\partial }^2{p}_1}{\partial r^2}+\frac{d_{f_1}-d-{\theta }_1+n}{r}\frac{\partial p_1}{\partial r}=n{G}_1{\varphi }_1{C}_{t_1}{\left(\frac{r}{r_w}\right)}^{{}^{{\theta }_1+1-n}}\frac{\partial p_1}{\partial t}{r}_w<r<\alpha r_w,t>0,$ (4.5)

同理可得外区渗流微分方程

$\frac{{\partial }^2{p}_2}{\partial r^2}+\frac{d_{f_2}-d-{\theta }_2+n}{r}\frac{\partial p_2}{\partial r}=n{G}_2{\varphi }_2{C}_{t_2}{\left(\frac{r}{r_w}\right)}^{{\theta }_2+1-n}\frac{\partial p_2}{\partial t}{r}_w<r<\alpha r_w,t>0,$ (4.6)

其中，$C_{t_i}=C_{{\rho }_i}+C_{{\varphi }_i}$ ，$G_i=\frac{{\mu }_e}{k_i}{\left(\frac{q}{2\pi h}\right)}^{n-1}\left(i=1,2\right)$ 。

4) 初始条件：

$p_1\left(r,0\right)=p_2\left(r,0\right)=p_0$ (4.7)

5) 考虑井筒储集系数C和表皮因子S的内边界条件：

$\{\begin{array}{l}{p}_w\left(t\right)={\left[p_1\left(r,t\right)-Sr\frac{\partial p_1}{\partial r}\right]|}_{r=r_w}\\ {\left(r\frac{\partial p_1}{\partial r}\right)|}_{r=r_w}={\frac{{\mu }_1}{2\pi k_{1}h}\left(Bq+C\frac{\text{d}{p}_w}{\text{d}t}\right)|}_{{}_{r=r_w}}\end{array}$ (4.8)

交界面条件：

${\{\begin{array}{l}{p}_1\left(\alpha r_w,t\right)=p_2\left(\alpha r_w,t\right)\\ {\frac{\partial p_1}{\partial r}|}_{r=\alpha r_w}=\frac{{\mu }_1{k}_2}{{\mu }_2{k}_1}\left[{\left(\frac{r}{r_w}\right)}^{d_{f_2}-d_{f_1}-\left({\theta }_2-{\theta }_1\right)}\frac{\partial p_2}{\partial r}\right]|\end{array}}_{r=\alpha r_w}$ (4.9)

6) 弹性外边界条件

${\left[\left({\chi }_1+\frac{{\chi }_2}{r_w^2}{r}^2+\frac{{\chi }_3}{{\left(n{\varphi }_1{C}_{t_1}{G}_1{r}_w^2\right)}^2}{t}^2\right)\left(p_0-p\right)+r\frac{\partial \left(p_0-p\right)}{\partial r}\right]|}_{r=R}=0$ (4.10)

其中，${\chi }_1,{\chi }_2,{\chi }_3$ 为任意常数。因此，得到考虑井筒储集和表皮系数的分形复合非牛顿幂律流体渗流模型

$\{\begin{array}{l}\frac{{\partial }^2{p}_1}{\partial r^2}+\frac{d_{f_1}-d-{\theta }_1+n}{r}\frac{\partial p_1}{\partial r}=n{G}_1{\varphi }_1{C}_{t_1}{\left(\frac{r}{r_w}\right)}^{{}^{{\theta }_1+1-n}}\frac{\partial p_1}{\partial t},r_w<r<\alpha r_w,t>0,\\ \frac{{\partial }^2{p}_2}{\partial r^2}+\frac{d_{f_2}-d-{\theta }_2+n}{r}\frac{\partial p_2}{\partial r}=n{G}_2{\varphi }_2{C}_{t_2}{\left(\frac{r}{r_w}\right)}^{{}^{{}^{{\theta }_2+1-n}}}\frac{\partial p_2}{\partial t},\alpha r_w<r<R,t>0,\\ p_1\left(r,0\right)=p_2\left(r,0\right)=p_0,p_w\left(t\right)={\left[p_1\left(r,t\right)-Sr\frac{\partial p_1}{\partial r}\right]|}_{r=r_w},\\ {\left(r\frac{\partial p_1}{\partial r}\right)|}_{r=r_w}={\frac{{\mu }_1}{2\pi k_{1}h}\left(Bq+C\frac{\text{d}{p}_w}{\text{d}t}\right)|}_{{}_{r=r_w}},\\ p_1\left(\alpha r_w,t\right)=p_2\left(\alpha r_w,t\right)\\ {\frac{\partial p_1}{\partial r}|}_{r=\alpha r_w}={\frac{{\mu }_1{k}_2}{{\mu }_2{k}_1}\left[{\left(\frac{r}{r_w}\right)}^{d_{f_2}-d_{f_1}-({\theta }_2-{\theta }_1)}\frac{\partial p_2}{\partial r}\right]|}_{r=\alpha r_w},\\ {\left[\left({\chi }_1+\frac{{\chi }_2}{r_w^2}{r}^2+\frac{{\chi }_3}{{\left(n{\varphi }_1{C}_{t_1}{G}_1{r}_w^2\right)}^2}{t}^2\right)p_2+r\frac{\partial p_2}{\partial r}\right]|}_{r=R}=0.\end{array}$ (4.11)

4.2. 模型的求解

1) 无量纲化

引入下列无量纲变换：

$p_{iD}=\frac{2\pi k_{1}h}{Bq{\mu }_1}\left(p_0-p_i\right)\left(i=1,2\right)$ ，$p_{wD}=\frac{2\pi k_{1}h}{Bq{\mu }_1}\left(p_0-p_w\right)$ ，$t_D=\frac{t}{n{\varphi }_1{C}_{t_1}{G}_1{r}_w^2}$ ，$r_D=\frac{r}{r_w}$ ，$R_D=\frac{R}{r_w}$ ，$C_D=\frac{C}{2\pi {\varphi }_1{C}_{t_1}h{r}_w^2}$ ，$\lambda =\frac{k_2{\mu }_1}{k_1{\mu }_2}$ ，$\omega =\frac{G_2{\varphi }_2{C}_{t_2}}{G_1{\varphi }_1{C}_{t_1}}$ 。

将式(4.11)无量纲化后，模型可化简为

$\{\begin{array}{l}\frac{{\partial }^2{p}_{1D}}{\partial r_D^2}+\frac{d_{f_1}-d-{\theta }_1+n}{r_D}\frac{\partial p_{1D}}{\partial r_D}=r_D^{{\theta }_1+1-n}\frac{\partial p_{1D}}{\partial t_D},1<r_D<\alpha ,t_D>0,\\ \frac{{\partial }^2{p}_{2D}}{\partial r_D^2}+\frac{d_{f_2}-d-{\theta }_2+n}{r_D}\frac{\partial p_{2D}}{\partial r_D}=\omega r_D^{{\theta }_2+1-n}\frac{\partial p_{2D}}{\partial t_D},\alpha <r_D<R_D,t_D>0,\\ p_{1D}\left(r_D,0\right)=p_{2D}\left(r_D,0\right)=0,p_{wD}\left(t_D\right)={\left[p_{1D}\left(r_D,t_D\right)-S{r}_D\frac{\partial p_{1D}}{\partial r_D}\right]|}_{r_D=1}\\ {\left(r_D\frac{\partial p_{1D}}{\partial r_D}\right)|}_{r_D=1}=-1+C_D\frac{\text{d}{p}_{wD}}{\text{d}{t}_D}\\ p_{1D}\left(\alpha ,t_D\right)=p_{2D}\left(\alpha ,t_D\right)\\ {\frac{\partial p_{1D}}{\partial r_D}|}_{r=\alpha }={\lambda \left[{\left(r_D\right)}^{d_{f_2}-d_{f_1}-({\theta }_2-{\theta }_1)}\frac{\partial p_{2D}}{\partial r_D}\right]|}_{r_D=\alpha }\\ {\left[\left({\chi }_1+{\chi }_2{r}_D^2+{\chi }_{{}_3}{t}_D^2\right)p_{2D}+r_D\frac{\partial p_{2D}}{\partial r_D}\right]|}_{r_D=R_D}=0\end{array}$ (4.12)

2) Laplace变换

对式(4.12)进行如下关于无量纲时间$t_D$ 的laplace变换

${\overline{p}}_{1D}\left(r_D,z\right)={\displaystyle {\int }_0^{+\infty }{\text{e}}^{-z{t}_D}}{p}_{1D}\left(r_D,t_D\right){\text{d}}_{t_D},$

${\overline{p}}_{2D}\left(r_D,z\right)={\displaystyle {\int }_0^{+\infty }{\text{e}}^{-z{t}_D}}{p}_{2D}\left(r_D,t_D\right){\text{d}}_{t_D},$

${\overline{p}}_{wD}\left(z\right)={\displaystyle {\int }_0^{+\infty }{\text{e}}^{-z{t}_D}}{p}_{wD}\left(t_D\right){\text{d}}_{t_D},$

即可得到laplace空间中的渗流模型：

$\{\begin{array}{l}\frac{{\text{d}}^2{\overline{p}}_{1D}}{\text{d}{r}_D^2}+\frac{d_{f_1}-d-{\theta }_1+n}{r_D}\frac{\text{d}{\overline{p}}_{1D}}{\text{d}{r}_D}=z{r}_D^{{\theta }_1+1-n}{\overline{p}}_{1D},1<r_D<\alpha ,t_D>0,\\ \frac{{\text{d}}^2{\overline{p}}_{2D}}{\text{d}{r}_D^2}+\frac{d_{f_2}-d-{\theta }_2+n}{r_D}\frac{\text{d}{\overline{p}}_{2D}}{\text{d}{r}_D}=\omega z{r}_D^{{\theta }_2+1-n}{\overline{p}}_{2D},\alpha <r_D<R_D,t_D>0,\\ {\left[-C_{D}z{\overline{p}}_{1D}+\left(1+C_{D}Sz\frac{\text{d}{\overline{p}}_{1D}}{\text{d}{r}_D}\right)\right]|}_{r_D=1}=-\frac{1}{z}\\ {\overline{p}}_{1D}\left(\alpha ,z\right)={\overline{p}}_{2D}\left(\alpha ,z\right)\\ {\frac{\text{d}{\overline{p}}_{1D}}{\text{d}{r}_D}|}_{r_D=\alpha }={\lambda \left[{\left(r_D\right)}^{d_{f_2}-d_{f_1}-({\theta }_2-{\theta }_1)}\frac{\text{d}{\overline{p}}_{2D}}{\text{d}{r}_D}\right]|}_{r_D=\alpha }\\ {\left[\left({\chi }_1+{\chi }_2{r}_D^2\right){\overline{p}}_{2D}+{\chi }_3\frac{{\text{d}}^2{\overline{p}}_{2D}}{\text{d}{z}^2}+r_D\frac{\text{d}{\overline{p}}_{2D}}{\text{d}{r}_D}\right]|}_{r=R_D}=0\end{array}$(4.13)

3) 相似构造法

模型(4.13)为一类复合变型bessel方程的边值问题，可用本文第二章阐述的相似构造法求解，在其边界条件中，$D=-\frac{1}{z}$ ，$E=-C_{D}z$ ，$F=-S$ ；$a=1,b=\alpha ,c=R_D$ 。

第一步：求内、外区定解方程的基础解系

$\frac{{\text{d}}^2{\overline{p}}_{1D}}{\text{d}{r}_D^2}+\frac{d_{f_1}-d-{\theta }_1+n}{r_D}\frac{\text{d}{\overline{p}}_{1D}}{\text{d}{r}_D}=z{r}_D^{{\theta }_1-n+1}{\overline{p}}_{1D}$的两个线性无关解为

$r_D^{\frac{d+{\theta }_1-n-d_{f_1}+1}{2}}{K}_{v_1}\left(\frac{2\sqrt{z}}{{\theta }_1-n+3}{r}_D^{\frac{{\theta }_1-n+3}{2}}\right)$ ，

$r_D^{\frac{d+{\theta }_1-n-d_{f_1}+1}{2}}{I}_{v_1}\left(\frac{2\sqrt{z}}{{\theta }_1-n+3}{r}_D^{\frac{{\theta }_1-n+3}{2}}\right)$

$\frac{{\text{d}}^2{\overline{p}}_{2D}}{\text{d}{r}_D^2}+\frac{d_{f_2}-d-{\theta }_2+n}{r_D}\frac{\text{d}{\overline{p}}_{2D}}{\text{d}{r}_D}=\omega z{r}_D^{{\theta }_2-n+1}{\overline{p}}_{2D}$的两个线性无关解为

$r_D^{\frac{d+{\theta }_2-n-d_{f_2}+1}{2}}{K}_{v_2}\left(\frac{2\sqrt{\omega z}}{{\theta }_2-n+3}{r}_D^{\frac{{\theta }_2-n+3}{2}}\right)$ ，

$r_D^{\frac{d+{\theta }_2-n-d_{f_2}+1}{2}}{I}_{v_2}\left(\frac{2\sqrt{\omega z}}{{\theta }_2-n+3}{r}_D^{\frac{{\theta }_2-n+3}{2}}\right)$

其中，$v_i=1+\frac{d-d_{f_i}-2}{{\theta }_i-n+3},i=1,2.$

第二步：根据基础解系构造引解函数：

$\begin{array}{l}{G}_{0,0}^i\left(r_D,R_D\right)=r_D^{\frac{d+{\theta }_i-n-d_{f_i}+1}{2}}{K}_{v_i}\left(\frac{2\sqrt{z_i}}{{\theta }_i-n+3}{r}_D^{\frac{{\theta }_i-n+3}{2}}\right)R_D^{\frac{d+{\theta }_i-n-d_{f_i}+1}{2}}{I}_{v_i}\left(\frac{2\sqrt{z_i}}{{\theta }_i-n+3}{R}_D^{\frac{{\theta }_i-n+3}{2}}\right)\\ -R_D^{\frac{d+{\theta }_i-n-d_{f_i}+1}{2}}{K}_{v_i}\left(\frac{2\sqrt{z_i}}{{\theta }_i-n+3}{R}_D^{\frac{{\theta }_i-n+3}{2}}\right)r_D^{\frac{d+{\theta }_i-n-d_{f_i}+1}{2}}{I}_{v_i}\left(\frac{2\sqrt{z_i}}{{\theta }_i-n+3}{r}_D^{\frac{{\theta }_i-n+3}{2}}\right)\\ ={\left(r_D{R}_D\right)}^{\frac{d+{\theta }_i-n-d_{f_i}+1}{2}}{\psi }_{v_i,v_i}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)\end{array}$

$\begin{array}{c}{G}_{0,1}^i\left(r_D,R_D\right)=\frac{\partial }{\partial R_D}{G}_{0,0}^i\left(r_D,R_D\right)\\ =r_D^{\frac{d+{\theta }_i-n-d_{f_i}+1}{2}}{R}_D^{\frac{d+{\theta }_i-n-d_{f_i}-1}{2}}\\ [\left(d+{\theta }_i-n-d_{f_i}+1\right){\psi }_{v_i,v_i}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)+\sqrt{z_i}{R}_D^{\frac{{\theta }_i-n+3}{2}}{\psi }_{v_i,v_i+1}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)]\end{array}$

$\begin{array}{c}{G}_{1,0}^i\left(r_D,R_D\right)=\frac{\partial }{\partial r_D}{G}_{0,0}^i\left(r_D,R_D\right)\\ =r_D^{\frac{d+{\theta }_i-n-d_{f_i}-1}{2}}{R}_D^{\frac{d+{\theta }_i-n-d_{f_i}+1}{2}}\\ [\left(d+{\theta }_i-n-d_{f_i}+1\right){\psi }_{v_i,v_i}\left(r_D^{{}^{\frac{{\theta }_i-n+3}{2}}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)+\sqrt{z_i}{r}_D^{\frac{{\theta }_i-n+3}{2}}{\psi }_{v_i+1,v_i}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)]\end{array}$

$\begin{array}{c}{G}_{1,1}^i\left(r_D,R_D\right)=\frac{\partial }{\partial r_D}{G}_{0,1}^i\left(r_D,R_D\right)=\frac{\partial }{\partial R_D}{G}_{1,0}^i\left(r_D,R_D\right)\\ ={\left(r_D{R}_D\right)}^{\frac{d+{\theta }_i-n-d_{f_i}-1}{2}}\left[{\left(d+{\theta }_i-n-d_{f_i}+1\right)}^2{\psi }_{v_i,v_i}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)\right)\\ +\sqrt{z_i}\left(d+{\theta }_i-n-d_{f_i}+1\right)R_D^{\frac{{\theta }_i-n+3}{2}}{\psi }_{v_i,v_i+1}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)\\ -\sqrt{z_i}\left(d+{\theta }_i-n-d_{f_i}+1\right)r_D^{\frac{{\theta }_i-n+3}{2}}{\psi }_{v_i+1,v_i}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)\\ -z_i{\left(R_D{r}_D\right)}^{\frac{{\theta }_i-n+3}{2}}{\psi }_{v_i+1,v_i+1}\left(r_D^{\frac{{\theta }_i-n+3}{2}},R_D^{\frac{{\theta }_i-n+3}{2}},\frac{2\sqrt{z_i}}{{\theta }_i-n+3}\right)]\end{array}$

$\begin{array}{c}{G}_{0,1}^{*}\left(r_D,R_D\right)=\frac{\partial }{\partial R_D}{G}_{0,1}^2\left(r_D,R_D\right)\\ ={\left(r_D{R}_D\right)}^{{}^{\frac{d+{\theta }_2-n-d_{f_2}+1}{2}}}[\left(d_{f_2}-d-{\theta }_2+n\right)\left(d_{f_2}-d-{\theta }_2+n-1\right)R_D^{-2}{\psi }_{v_2,v_2}\left(r_D^{\frac{{\theta }_2-n+3}{2}},R_D^{\frac{{\theta }_2-n+3}{2}},\frac{2\sqrt{z_2}}{{\theta }_2-n+3}\right)\\ +\left(3{\theta }_2-3n-2d_{f_2}+2d+4\right)\sqrt{z_2}{R}_D^{\frac{{\theta }_2-n-1}{2}}{\psi }_{v_2,v_2+1}\left(r_D^{{}^{\frac{{\theta }_2-n+3}{2}}},R_D^{\frac{{\theta }_2-n+3}{2}},\frac{2\sqrt{z_2}}{{\theta }_2-n+3}\right)\\ +z_2{R}_D^{{\theta }_2-n-1}{\psi }_{v_2,v_2+2}\left(r_D^{\frac{{\theta }_2-n+3}{2}},R_D^{\frac{{\theta }_2-n+3}{2}},\frac{2\sqrt{z_2}}{{\theta }_2-n+3}\right)]\end{array}$

$\begin{array}{l}{G}_{1,1}^{\#}\left(r_D,R_D\right)=\frac{\partial }{\partial R_D}{G}_{1,1}^2\left(r_D,R_D\right)\\ =\frac{{\partial }^3\left[r_D^{\frac{d+{\theta }_2-n-d_{f_2}+1}{2}}{K}_{v_2}\left(\frac{2\sqrt{z_2}}{{\theta }_2-n+3}{r}_D^{{}^{\frac{{\theta }_2-n+3}{2}}}\right)R_D^{\frac{d+{\theta }_2-n-d_{f_2}+1}{2}}{I}_{v_2}\left(\frac{2\sqrt{z_2}}{{\theta }_2-n+3}{R}_D^{{}^{\frac{{\theta }_2-n+3}{2}}}\right)\right]}{\partial r_D\partial R_D{}^2}\\ -\frac{{\partial }^3\left[R_D^{\frac{d+{\theta }_2-n-d_{f_2}+1}{2}}{K}_{v_2}\left(\frac{2\sqrt{z_2}}{{\theta }_2-n+3}{R}_D^{\frac{{\theta }_2-n+3}{2}}\right)r_D^{\frac{d+{\theta }_2-n-d_{f_2}+1}{2}}{I}_{v_2}\left(\frac{2\sqrt{z_2}}{{\theta }_2-n+3}{r}_D^{\frac{{\theta }_2-n+3}{2}}\right)\right]}{\partial r_D\partial R_D^2}\end{array}$

其中

$z_1=z$ ，$z_2=\omega z$ ，${\psi }_{m,n}\left(x,y,z\right)=K_m\left(xz\right)I_n\left(yz\right)+{\left(-1\right)}^{m-n+1}{I}_m\left(xz\right)K_n\left(yz\right)$ ，

且${\psi }_{m,n}\left(x,y,z\right)$ 是变型的$n$ 阶的Bessel函数$K_n\left(yz\right)$ ，$I_n\left(yz\right)$ 的线性组合，$m$ ，$n$ 是常数。

第三步：由边界条件构造相似核函数

边界条件可以转换为以下形式

${\left[-C_{D}z{\overline{p}}_{1D}+\left(1+C_{D}Sz\frac{\text{d}{\overline{p}}_{1D}}{\text{d}{r}_D}\right)\right]|}_{r_D=1}=-\frac{1}{z}$

结合引解函数、交界面条件、外边界条件可构造出外区相似核函数${\Phi }_2\left(r_D\right)$ 和内区相似核函数${\Phi }_1\left(r_D\right)$ ，即

${\Phi }_1\left(r_D\right)=\frac{{\Phi }_2\left(\alpha ,z\right)G_{0,1}^1\left(r_D,\alpha \right)-\lambda {\alpha }^{{}^{d_{f_2}-d_{f_1}-({\theta }_2-{\theta }_1)}}{G}_{0,0}^1\left(r_D,\alpha \right)}{{\Phi }_2\left(\alpha ,z\right)G_{1,1}^1\left(1,\alpha \right)-\lambda {\alpha }^{{}^{d_{f_2}-d_{f_1}-({\theta }_2-{\theta }_1)}}{G}_{1,0}^1\left(1,\alpha \right)}$

${\Phi }_2\left(r_D\right)=\frac{\left({\chi }_1+{\chi }_2{r}_D^2\right)G_{0,0}^2\left(r_D,R_D\right)+{\chi }_3{G}_{0,1}^{*}\left(r_D,R_D\right)+R_D{G}_{0,1}^2\left(r_D,R_D\right)}{\left({\chi }_1+{\chi }_2{r}_D^2\right)G_{1,0}^2\left(\alpha ,R_D\right)+{\chi }_3{G}_{1,1}^{\#}\left(\alpha ,R_D\right)+R_D{G}_{1,1}^2\left(\alpha ,R_D\right)}$

第四步：得到内、外区的解

结合相似核函数、内边界条件可得到内区的解${\overline{p}}_{1D}\left(r_D,z\right)$ 和外区的解${\overline{p}}_{2D}\left(r_D,R_D\right)$ ，即

${\overline{p}}_{1D}\left(r_D,z\right)=-\frac{1}{z}\cdot \frac{1}{-C_{D}z+\frac{1}{-S+{\Phi }_1\left(1,z\right)}}\cdot \frac{1}{-S+{\Phi }_1\left(1,z\right)}\cdot {\Phi }_1\left(r_D,z\right)$ (4.14)

$\begin{array}{l}{\overline{p}}_{2D}\left(r_D,R_D\right)=-\frac{1}{z}\cdot \frac{1}{-C_{D}z+\frac{1}{-S+{\Phi }_1\left(1\right)}}\cdot \frac{1}{-S+{\Phi }_1\left(1\right)}\cdot \\ =\frac{G_{0,0}^1\left(\alpha ,\alpha \right)}{{\Phi }_2\left(\alpha \right)G_{1,1}^1\left(1,\alpha \right)-\chi {\alpha }^{d_{f_2}-d_{f_1}-({\theta }_2-{\theta }_1)}{G}_{1,0}^1\left(1,\alpha \right)}\cdot {\Phi }_2\left(r_D,z\right)\\ \end{array}$(4.15)

根据模型里的内边界条件${\overline{p}}_{wD}\left(z\right)={\left[{\overline{p}}_{1D}\left(r_D,t_D\right)-S{r}_D\frac{\text{d}{p}_{1D}}{\text{d}{r}_D}\right]|}_{r_D=1}$，可以得到无量纲井底压力在laplace空间中的解，即

${\overline{p}}_{wD}\left(r_D,z\right)=\frac{1}{z}\cdot \frac{1}{C_{D}z-\frac{1}{{\Phi }_1\left(1,z\right)-S}}$ (4.16)

特别地，当$C_D=0$ ，即为不考虑井筒储集的情况；当$S=0$ 为不考虑表皮效应的情况。

4) Stehfest数值反演

利用Gaver-Stehfest数值反演，可得到无量纲井底压力在实空间中的解。

$P_{wD}\left(t_D\right)=\frac{\ln 2}{t_D}{\displaystyle \sum _{j=1}^N{V}_j}{\overline{P}}_{wD}\left(z\right)=\frac{\ln 2}{t_D}{V}_j{\overline{P}}_{wD}\left(\frac{j\ln 2}{t}\right)$ (4.17)

其中

$V_j={\left(-1\right)}^{N/2+j}{\displaystyle \sum _{k=\left[\frac{j+1}{2}\right]}^{\min \left(j,\frac{N}{2}\right)}\frac{k^{\frac{N}{2}}\left(2k+1\right)!}{k!\left(k+1\right)！\left(\frac{N}{2}-k+1\right)!\left(j-k+1\right)!\left(2k-j+1\right)!}}$ (4.18)

5. 结论

本文建立的弹性外边界分形复合非牛顿幂律流体渗流模型，通过纳入井筒储集与表皮因子，有效弥补了传统模型忽略地层弹性、简化边界假设的局限，显著提升了对实际储层渗流状态的适配性。基于模型结果，得出关键新发现：弹性外边界会使油藏压力导数曲线晚期呈现“上翘后趋平”特征，与传统理想边界模型的单调变化规律存在本质差异，且分形特性与非牛顿流体流变特征会加剧该差异，为油藏边界识别及试井解释提供了全新理论支撑。

该模型可精准揭示此类油藏的真实渗流规律，为其试井解释、开发方案优化提供可靠理论依据，彰显出较高的学术价值与工程应用潜力。

基金项目

四川省科技厅科技计划项目(2015JY0245)。