这是用户在 2025-3-12 12:07 为 https://app.immersivetranslate.com/pdf-pro/e3f9265e-86cc-4d7d-9a0c-4c39d7283e7f/ 保存的双语快照页面,由 沉浸式翻译 提供双语支持。了解如何保存?


离岸重油储层中多井干扰下倾斜井的压力性能分析


中国海洋石油总公司天津分公司,渤海油田研究院,天津,300459

请将信件发送至王公昌;wanggch12@cnooc.com.cn


收到日期:2022 年 4 月 11 日;修订日期:2022 年 5 月 22 日;接受日期:2022 年 6 月 6 日;发布日期:2022 年 7 月 13 日


学术编辑:刘书扬


版权 © 2022 Kuiqian Ma 等人。本文为开放获取文章,依据知识共享署名许可协议分发,允许在任何媒介中无限制使用、分发和复制,前提是正确引用原作。

海上重油储层主要通过倾斜井开发,经过几次井网调整后,井距不断缩小。井测试数据常常受到邻近井的干扰,导致井测试解释结果不理想,因此有必要对倾斜井的多井干扰井测试问题进行研究。本文综合考虑了重油的阈值压力梯度和对渗透率敏感的应力,建立了多井干扰下双介质储层倾斜井的井测试模型。通过使用格林函数和叠加原理获得了拉普拉斯空间中的解析解。结果表明,典型的压力动态曲线可以分为 12 个流动区域,中央倾斜井的压力导数曲线在邻近井的干扰下向上翘起并形成多个“平台”;邻近井的干扰将加剧压力导数曲线的上升;倾斜井的临界井倾角为 40 40 40^(@)40^{\circ} 。 当井倾角大于 40 40 40^(@)40^{\circ} 时,将出现类似于水平井的垂直径向流动。新模型在渤海湾的 BZ 油田应用中得到了很好的匹配和解释,为类似储层的多井干扰试井提供了理论指导。

  1. 引言


目前常见的井测试解释方法假设在油藏中仅存在一个中心测试井。随着油田开发到后期,油藏的连通性变得更加复杂,邻近井之间的干扰变得更加明显。测试井的压力动态曲线常常受到邻近井的干扰,导致后期径向流动段“上翘”,这通常在常见的单井测试模型中被视为边界影响,从而导致井测试数据的错误处理[1-3]。Onur 等人提出了一种包括多产井系统的压力恢复模型[4]。Marhaendrajana 等人开发了一种通过将“干扰效应”视为区域压力下降来解释多井系统中压力的方法[5]。Adewole 基于干扰下测试井的压力数据评估了井之间的连通性[6]。Deng 等人建立了一种在多井干扰下进行压力恢复的解释方法[7]。Cheng 等人使用多井干扰方法来判断水


水平井的流入方向 [8]。杨等人开发了一种多段水平井(MSHW)的新型干扰测试模型,以更好地理解注入井在地平观测井开放生产时的干扰 [9]。库马尔等人建立了一个数学干扰测试模型,用于裂缝孔隙碳酸盐岩储层,然后用于测试塔里木油田的一个观测井和两个干扰井 [10]。韩等人提出了一种基于示踪剂和压力干扰数据分析的综合方法,以获取多井平台中裂缝水平井之间的干扰程度 [11]。石青等人绘制了通过大裂缝连接的井的干扰测试的压力类型曲线及压力导数,并与干扰测试数据进行了验证 [12]。

由于海上平台的特殊特性,倾斜井能够适应复杂的地质问题,并增加油渗透面积以提高井的产量,因此现在在海上油田中得到广泛应用。Cinco 等人首次分析并绘制了一个的压力动态。


通过建立一个上下闭合的测试井模型来研究倾斜井 [13]。Zhang 等人研究了地层非均匀性对倾斜井压力的影响 [14]。Sousa 研究了均质储层中倾斜井的压力动态 [15]。Li 等人建立了属于复合气藏中具有常规内区和低渗透外区的倾斜井,综合考虑了低渗透复合气藏的应力敏感性和非达西流动特征 [16]。

总之,目前尚无关于相邻井对倾斜井压力动态干扰的研究。然而,由于倾斜井目前是海上油田最主要的开发方法,因此研究适用于倾斜井的相邻井干扰试井模型是非常重要的。本文考虑了重油的渗透率应力敏感性和阈值压力梯度,建立了海上重油储层中多井干扰下的倾斜井试井模型,并分析了多种敏感参数的影响。最后,该模型在渤海湾的 SZ 油田得到了良好的应用。

  2. 模型开发


2.1. 重油的非线性渗流。重油具有高粘度、大渗流阻力以及液固界面和液液界面之间的较大相互作用力[17-19]。因此,孔隙介质中的渗流特性与常规轻油不同,通常表现出非线性渗流特性(阈值压力梯度)[20-23]。只有当位移压力梯度超过阈值压力梯度时,重油才开始流动,其渗流特性如图 1 所示。

渤海 SZ 油田的核心位移实验表明,重油在多孔介质中的阈值压力梯度和流动性符合非线性关系,如图 2 所示。当流动性较小时,压力梯度随着流动性的增加而迅速下降。随着流动性的持续增加,阈值压力梯度的下降速度减缓,并与指数函数相匹配,如公式(1)所示。这种现象的原因在于,随着重油粘度的降低,重油中胶质、沥青质和高分子烃的含量减少,导致重油的结构特征减弱,流动过程中的分子间力减小,流动阻力降低。随着渗透率的增加,阻力减少的速度更快,


图 1:重油渗漏特性动态曲线。


导致随着流动性的增加,重油阈值压力梯度逐渐降低。
λ B = a ( K μ o ) b λ B = a K μ o b lambda_(B)=a((K)/(mu_(o)))^(-b)\lambda_{\mathrm{B}}=a\left(\frac{K}{\mu_{\mathrm{o}}}\right)^{-b}

2.2. 物理模型。如图 3 所示,在无限外边界储层中有一个中央测试倾斜井,周围的邻近井与其具有良好的连通性:


(1) 采用 Warren-Root 模型来描述双孔隙性地层


(2)孔隙间流动是通过伪稳态模型计算的


(3) 考虑了渗透率的应力敏感性


(4) 多孔介质中的重油具有阈值压力梯度的特性


2.3. 数学模型。考虑到储层的渗透率应力敏感性和重油流体的阈值压力梯度,流体运动方程改进如下:
v fr = K fh μ e γ ( p i p f ) ( p f r λ B ) , v fz = K fv μ e γ ( p i p f ) ( p f r λ B ) v fr = K fh μ e γ p i p f p f r λ B , v fz = K fv μ e γ p i p f p f r λ B {:[v_(fr)=(K_(fh))/(mu)e^(-gamma(p_(i)-p_(f)))((delp_(f))/(del r)-lambda_(B))","],[v_(fz)=(K_(fv))/(mu)e^(-gamma(p_(i)-p_(f)))((delp_(f))/(del r)-lambda_(B))]:}\begin{gathered} v_{\mathrm{fr}}=\frac{K_{\mathrm{fh}}}{\mu} e^{-\gamma\left(p_{\mathrm{i}}-p_{\mathrm{f}}\right)}\left(\frac{\partial p_{\mathrm{f}}}{\partial r}-\lambda_{\mathrm{B}}\right), \\ v_{\mathrm{fz}}=\frac{K_{\mathrm{fv}}}{\mu} e^{-\gamma\left(p_{\mathrm{i}}-p_{\mathrm{f}}\right)}\left(\frac{\partial p_{\mathrm{f}}}{\partial r}-\lambda_{\mathrm{B}}\right) \end{gathered}

根据方程(1)和方程(2),得到描述应力敏感储层中重油阈值压力梯度的中央测试倾斜井的渗流微分方程:
{ K ff μ 1 r r [ r e γ ( p i p f ) ( p f r λ B ) ] + K fv μ z [ e γ ( p i p f ) ( p f z λ B ) ] = ( ϕ C t V ) f p f t α K m μ ( p m p f ) ( ϕ C t V ) m p m t + α K m μ ( p m p f ) = 0 K ff μ 1 r r r e γ p i p f p f r λ B + K fv μ z e γ p i p f p f z λ B = ϕ C t V f p f t α K m μ p m p f ϕ C t V m p m t + α K m μ p m p f = 0 {[(K_(ff))/(mu)(1)/(r)(del)/(del r)[re^(-gamma(p_(i)-p_(f)))((delp_(f))/(del r)-lambda_(B))]+(K_(fv))/(mu)(del)/(del z)[e^(-gamma(p_(i)-p_(f)))((delp_(f))/(del z)-lambda_(B))]=(phiC_(t)V)_(f)(delp_(f))/(del t)-(alphaK_(m))/(mu)(p_(m)-p_(f))],[(phiC_(t)V)_(m)(delp_(m))/(del t)+(alphaK_(m))/(mu)(p_(m)-p_(f))=0]:}\left\{\begin{array}{l} \frac{K_{\mathrm{ff}}}{\mu} \frac{1}{r} \frac{\partial}{\partial r}\left[r e^{-\gamma\left(p_{\mathrm{i}}-p_{\mathrm{f}}\right)}\left(\frac{\partial p_{\mathrm{f}}}{\partial r}-\lambda_{\mathrm{B}}\right)\right]+\frac{K_{\mathrm{fv}}}{\mu} \frac{\partial}{\partial z}\left[e^{-\gamma\left(p_{\mathrm{i}}-p_{\mathrm{f}}\right)}\left(\frac{\partial p_{\mathrm{f}}}{\partial z}-\lambda_{\mathrm{B}}\right)\right]=\left(\phi C_{\mathrm{t}} V\right)_{\mathrm{f}} \frac{\partial p_{\mathrm{f}}}{\partial t}-\frac{\alpha K_{\mathrm{m}}}{\mu}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right) \\ \left(\phi C_{\mathrm{t}} V\right)_{\mathrm{m}} \frac{\partial p_{\mathrm{m}}}{\partial t}+\frac{\alpha K_{\mathrm{m}}}{\mu}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right)=0 \end{array}\right.

图 2:压力梯度与流动性的关系曲线。

以下是相邻干扰井的微分方程:
{ K fh μ 1 r r [ r e γ ( p i p f ) ( p f r λ B ) ] = ( ϕ C t V ) f p f t α K m μ ( p m p f ) ( ϕ C t V ) m p m t + α K m μ ( p m p f ) = 0 K fh μ 1 r r r e γ p i p f p f r λ B = ϕ C t V f p f t α K m μ p m p f ϕ C t V m p m t + α K m μ p m p f = 0 {[(K_(fh))/(mu)(1)/(r)(del)/(del r)[re^(-gamma(p_(i)-p_(f)))((delp_(f))/(del r)-lambda_(B))]=(phiC_(t)V)_(f)(delp_(f))/(del t)-(alphaK_(m))/(mu)(p_(m)-p_(f))],[(phiC_(t)V)_(m)(delp_(m))/(del t)+(alphaK_(m))/(mu)(p_(m)-p_(f))=0]:}\left\{\begin{array}{l} \frac{K_{\mathrm{fh}}}{\mu} \frac{1}{r} \frac{\partial}{\partial r}\left[r e^{-\gamma\left(p_{\mathrm{i}}-p_{\mathrm{f}}\right)}\left(\frac{\partial p_{\mathrm{f}}}{\partial r}-\lambda_{\mathrm{B}}\right)\right]=\left(\phi C_{\mathrm{t}} V\right)_{\mathrm{f}} \frac{\partial p_{\mathrm{f}}}{\partial t}-\frac{\alpha K_{\mathrm{m}}}{\mu}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right) \\ \left(\phi C_{\mathrm{t}} V\right)_{\mathrm{m}} \frac{\partial p_{\mathrm{m}}}{\partial t}+\frac{\alpha K_{\mathrm{m}}}{\mu}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right)=0 \end{array}\right.

由于 ( p f / r ) × λ B p f / r × λ B (delp_(f)//del r)xxlambda_(B)\left(\partial p_{\mathrm{f}} / \partial r\right) \times \lambda_{\mathrm{B}} ( p f / z ) × λ B p f / z × λ B (delp_(f)//delz)xxlambda_(B)\left(\partial p_{\mathrm{f}} / \partial \mathrm{z}\right) \times \lambda_{\mathrm{B}} 的值


中央测试倾斜井:小,可以进行四舍五入,上述方程可以简化如下。
{ K fh [ 1 r r ( r p f r ) + γ ( p f r ) 2 λ B r ] + K fv [ 2 p f z 2 + γ ( p f z ) 2 ] = e γ ( p t p f ) [ ( ϕ C t V ) f μ p f t α K m ( p m p f ) ] ( ϕ C t V ) m p m t + α K m μ ( p m p f ) = 0 K fh 1 r r r p f r + γ p f r 2 λ B r + K fv 2 p f z 2 + γ p f z 2 = e γ p t p f ϕ C t V f μ p f t α K m p m p f ϕ C t V m p m t + α K m μ p m p f = 0 {[K_(fh)[(1)/(r)(del)/(del r)(r(delp_(f))/(del r))+gamma((delp_(f))/(del r))^(2)-(lambda_(B))/(r)]+K_(fv)[(del^(2)p_(f))/(delz^(2))+gamma((delp_(f))/(del z))^(2)]=e^(gamma(p_(t)-p_(f)))[(phiC_(t)V)_(f)mu(delp_(f))/(del t)-alphaK_(m)(p_(m)-p_(f))]],[(phiC_(t)V)_(m)(delp_(m))/(del t)+(alphaK_(m))/(mu)(p_(m)-p_(f))=0]:}\left\{\begin{array}{l} K_{\mathrm{fh}}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial p_{\mathrm{f}}}{\partial r}\right)+\gamma\left(\frac{\partial p_{\mathrm{f}}}{\partial r}\right)^{2}-\frac{\lambda_{\mathrm{B}}}{r}\right]+K_{\mathrm{fv}}\left[\frac{\partial^{2} p_{\mathrm{f}}}{\partial z^{2}}+\gamma\left(\frac{\partial p_{\mathrm{f}}}{\partial z}\right)^{2}\right]=e^{\gamma\left(p_{\mathrm{t}}-p_{\mathrm{f}}\right)}\left[\left(\phi C_{\mathrm{t}} V\right)_{\mathrm{f}} \mu \frac{\partial p_{\mathrm{f}}}{\partial t}-\alpha K_{\mathrm{m}}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right)\right] \\ \left(\phi C_{\mathrm{t}} V\right)_{\mathrm{m}} \frac{\partial p_{\mathrm{m}}}{\partial t}+\frac{\alpha K_{\mathrm{m}}}{\mu}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right)=0 \end{array}\right.

邻近干扰井:
{ K fh [ 1 r r ( r p f r ) + γ ( p f r ) 2 λ B r ] = e γ ( p i p f ) [ ( ϕ C t V ) f μ p f t α K m ( p m p f ) ] ( ϕ C t V ) m p m t + α K m μ ( p m p f ) = 0 K fh 1 r r r p f r + γ p f r 2 λ B r = e γ p i p f ϕ C t V f μ p f t α K m p m p f ϕ C t V m p m t + α K m μ p m p f = 0 {[K_(fh)[(1)/(r)(del)/(del r)(r(delp_(f))/(del r))+gamma((delp_(f))/(del r))^(2)-(lambda_(B))/(r)]=e^(gamma(p_(i)-p_(f)))[(phiC_(t)V)_(f)mu(delp_(f))/(del t)-alphaK_(m)(p_(m)-p_(f))]],[(phiC_(t)V)_(m)(delp_(m))/(del t)+(alphaK_(m))/(mu)(p_(m)-p_(f))=0]:}\left\{\begin{array}{l} K_{\mathrm{fh}}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial p_{\mathrm{f}}}{\partial r}\right)+\gamma\left(\frac{\partial p_{\mathrm{f}}}{\partial r}\right)^{2}-\frac{\lambda_{\mathrm{B}}}{r}\right]=e^{\gamma\left(p_{\mathrm{i}}-p_{\mathrm{f}}\right)}\left[\left(\phi C_{\mathrm{t}} V\right)_{\mathrm{f}} \mu \frac{\partial p_{\mathrm{f}}}{\partial t}-\alpha K_{\mathrm{m}}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right)\right] \\ \left(\phi C_{\mathrm{t}} V\right)_{\mathrm{m}} \frac{\partial p_{\mathrm{m}}}{\partial t}+\frac{\alpha K_{\mathrm{m}}}{\mu}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right)=0 \end{array}\right.

  • 邻近干扰井


    _- 中心测试倾斜井


    图 3:物理模型图。

无量纲量定义如下:

p fD = ( 2 π K fl h / q μ ) ( p i p f ) p mD = ( 2 π K m h / q μ ) ( p i p m p fD = 2 π K fl h / q μ p i p f p mD = 2 π K m h / q μ p i p m p_(fD)=(2piK_(fl)h//q mu)(p_(i)-p_(f))p_(mD)=(2piK_(m)h//q mu)(p_(i)-p_(m):}p_{\mathrm{fD}}=\left(2 \pi K_{\mathrm{fl}} h / q \mu\right)\left(p_{\mathrm{i}}-p_{\mathrm{f}}\right) p_{\mathrm{mD}}=\left(2 \pi K_{\mathrm{m}} h / q \mu\right)\left(p_{\mathrm{i}}-p_{\mathrm{m}}\right. ) t D = ( K fh t / μ r w 2 ( V ϕ C t ) f + m ) ω = ( ( V ϕ C t ) f / ( V ϕ C t ) f + m ) λ m t D = K fh t / μ r w 2 V ϕ C t f + m ω = V ϕ C t f / V ϕ C t f + m λ m t_(D)=(K_(fh)t//mur_(w)^(2)(V phiC_(t))_(f+m))omega=((V phiC_(t))_(f)//(V phiC_(t))_(f+m))lambda_(m)t_{\mathrm{D}}=\left(K_{\mathrm{fh}} t / \mu r_{\mathrm{w}}^{2}\left(V \phi C_{t}\right)_{\mathrm{f}+\mathrm{m}}\right) \omega=\left(\left(V \phi C_{\mathrm{t}}\right)_{\mathrm{f}} /\left(V \phi C_{\mathrm{t}}\right)_{\mathrm{f}+\mathrm{m}}\right) \lambda_{\mathrm{m}} = α r w 2 ( K m / K fh ) γ D = ( q μ B / K fh h ) γ λ BD = ( K fh h r w / q μ B ) λ B = α r w 2 K m / K fh γ D = q μ B / K fh h γ λ BD = K fh h r w / q μ B λ B =alphar_(w)^(2)(K_(m)//K_(fh))gamma_(D)=(q mu B//K_(fh)h)gammalambda_(BD)=(K_(fh)hr_(w)//q mu B)lambda_(B)=\alpha r_{\mathrm{w}}^{2}\left(K_{\mathrm{m}} / K_{\mathrm{fh}}\right) \gamma_{D}=\left(q \mu B / K_{\mathrm{fh}} h\right) \gamma \lambda_{\mathrm{BD}}=\left(K_{\mathrm{fh}} h r_{\mathrm{w}} / q \mu B\right) \lambda_{\mathrm{B}}
r D = ( r / r w ) z D = ( z / h ) C D = ( C / ( V ϕ C t ) f + m h r w 2 ) L D = ( L / h ) r D = r / r w z D = ( z / h ) C D = C / V ϕ C t f + m h r w 2 L D = ( L / h ) r_(D)=(r//r_(w))z_(D)=(z//h)C_(D)=(C//(V phiC_(t))_(f+m)hrw^(2))L_(D)=(L//h)r_{\mathrm{D}}=\left(r / r_{\mathrm{w}}\right) z_{\mathrm{D}}=(z / h) C_{\mathrm{D}}=\left(C /\left(V \phi C_{t}\right)_{\mathrm{f}+\mathrm{m}} h r \mathrm{w}^{2}\right) L_{\mathrm{D}}=(L / h) K fh / K fv h D = ( h / r w ) K fh / K fv K fh / K fv h D = h / r w K fh / K fv sqrt(K_(fh)//K_(fv))h_(D)=(h//r_(w))sqrt(K_(fh)//K_(fv))\sqrt{K_{\mathrm{fh}} / K_{\mathrm{fv}}} h_{\mathrm{D}}=\left(h / r_{\mathrm{w}}\right) \sqrt{K_{\mathrm{fh}} / K_{\mathrm{fv}}}.

将上述无量纲量代入公式 (5) 和 (6),得到无量纲渗流微分方程:
{ 1 r D r D ( r D p fD r D ) γ D ( p fD r D ) + λ BD r D + [ 2 p fD z D + γ D ( p fD z D ) 2 ] = e γ D p fD [ ω p fD t D λ m ( p mD p fD ) ] , ω p mD t D + λ m ( p mD p fD ) = 0 , { 1 r jD r jD ( r jD p fD r jD ) γ D ( p fD r jD ) + λ BD r jD = e γ D p fD [ ω p fD t jD λ m ( p mD p fD ) ] , ω p mD t jD + λ m ( p mD p fD ) = 0 . : 1 r D r D r D p fD r D γ D p fD r D + λ BD r D + 2 p fD z D + γ D p fD z D 2 = e γ D p fD ω p fD t D λ m p mD p fD , ω p mD t D + λ m p mD p fD = 0 , 1 r jD r jD r jD p fD r jD γ D p fD r jD + λ BD r jD = e γ D p fD ω p fD t jD λ m p mD p fD , ω p mD t jD + λ m p mD p fD = 0 .  :  {:{[(1)/(r_(D))(del)/(delr_(D))(r_(D)(delp_(fD))/(delr_(D)))-gamma_(D)((delp_(fD))/(delr_(D)))+(lambda_(BD))/(r_(D))+[(del^(2)p_(fD))/(delz_(D))+gamma_(D)((delp_(fD))/(delz_(D)))^(2)]=e^(gamma_(D)p_(fD))[omega(delp_(fD))/(delt_(D))-lambda_(m)(p_(mD)-p_(fD))]","],[omega(delp_(mD))/(delt_(D))+lambda_(m)(p_(mD)-p_(fD))=0","],[{[(1)/(r_(jD))(del)/(delr_(jD))(r_(jD)(delp_(fD))/(delr_(jD)))-gamma_(D)((delp_(fD))/(delr_(jD)))+(lambda_(BD))/(r_(jD))=e^(gamma_(D)p_(fD))[omega(delp_(fD))/(delt_(jD))-lambda_(m)(p_(mD)-p_(fD))]","],[omega(delp_(mD))/(delt_(jD))+lambda_(m)(p_(mD)-p_(fD))=0.]:}]:}" : ":}\begin{aligned} & \left\{\begin{array}{l} \frac{1}{r_{\mathrm{D}}} \frac{\partial}{\partial r_{\mathrm{D}}}\left(r_{\mathrm{D}} \frac{\partial p_{\mathrm{fD}}}{\partial r_{\mathrm{D}}}\right)-\gamma_{D}\left(\frac{\partial p_{\mathrm{fD}}}{\partial r_{\mathrm{D}}}\right)+\frac{\lambda_{\mathrm{BD}}}{r_{\mathrm{D}}}+\left[\frac{\partial^{2} p_{\mathrm{fD}}}{\partial z_{\mathrm{D}}}+\gamma_{\mathrm{D}}\left(\frac{\partial p_{\mathrm{fD}}}{\partial z_{\mathrm{D}}}\right)^{2}\right]=e^{\gamma_{\mathrm{D}} p_{\mathrm{fD}}}\left[\omega \frac{\partial p_{\mathrm{fD}}}{\partial t_{\mathrm{D}}}-\lambda_{\mathrm{m}}\left(p_{\mathrm{mD}}-p_{\mathrm{fD}}\right)\right], \\ \omega \frac{\partial p_{\mathrm{mD}}}{\partial t_{\mathrm{D}}}+\lambda_{\mathrm{m}}\left(p_{\mathrm{mD}}-p_{\mathrm{fD}}\right)=0, \\ \left\{\begin{array}{l} \frac{1}{r_{\mathrm{jD}}} \frac{\partial}{\partial r_{\mathrm{jD}}}\left(r_{\mathrm{jD}} \frac{\partial p_{\mathrm{fD}}}{\partial r_{\mathrm{jD}}}\right)-\gamma_{D}\left(\frac{\partial p_{\mathrm{fD}}}{\partial r_{\mathrm{jD}}}\right)+\frac{\lambda_{\mathrm{BD}}}{r_{\mathrm{jD}}}=e^{\gamma_{\mathrm{D}} p_{\mathrm{fD}}}\left[\omega \frac{\partial p_{\mathrm{fD}}}{\partial t_{\mathrm{jD}}}-\lambda_{\mathrm{m}}\left(p_{\mathrm{mD}}-p_{\mathrm{fD}}\right)\right], \\ \omega \frac{\partial p_{\mathrm{mD}}}{\partial t_{\mathrm{jD}}}+\lambda_{\mathrm{m}}\left(p_{\mathrm{mD}}-p_{\mathrm{fD}}\right)=0 . \end{array}\right. \end{array}\right. \text { : } \end{aligned}

中央测试倾斜井如下:

  初始条件
p fD ( r D , 0 ) = p mD ( r D , 0 ) = 0 p fD r D , 0 = p mD r D , 0 = 0 p_(fD)(r_(D),0)=p_(mD)(r_(D),0)=0p_{\mathrm{fD}}\left(r_{\mathrm{D}}, 0\right)=p_{\mathrm{mD}}\left(r_{\mathrm{D}}, 0\right)=0
  外边界条件
p fD | r D = p mD | r D = 0 p fD r D = p mD r D = 0 p_(fD)|_(r_(D)rarr oo)=p_(mD)|_(r_(D)rarr oo)=0\left.p_{\mathrm{fD}}\right|_{r_{\mathrm{D}} \rightarrow \infty}=\left.p_{\mathrm{mD}}\right|_{r_{\mathrm{D}} \rightarrow \infty}=0
  内边界条件
lim ε 0 [ lim r D 0 z wD ε / 2 z wD + ( ε / 2 ) r D e γ D p fD ( p fD r D + λ BD ) d z D ] = h D | z D z wD | ε 2 lim ε 0 lim r D 0 z wD ε / 2 z wD + ( ε / 2 ) r D e γ D p fD p fD r D + λ BD d z D = h D z D z wD ε 2 lim_(epsi rarr0)[lim_(r_(D)longrightarrow0)int_(z_(wD)-epsi//2)^(z_(wD)+(epsi//2))r_(D)e^(-gamma_(D)p_(fD))((delp_(fD))/(delr_(D))+lambda_(BD))dz_(D)]=-h_(D)|z_(D)-z_(wD)| <= (epsi)/(2)\lim _{\varepsilon \rightarrow 0}\left[\lim _{r_{\mathrm{D}} \longrightarrow 0} \int_{z_{\mathrm{wD}}-\varepsilon / 2}^{z_{\mathrm{wD}}+(\varepsilon / 2)} r_{\mathrm{D}} e^{-\gamma_{\mathrm{D}} p_{\mathrm{fD}}}\left(\frac{\partial p_{\mathrm{fD}}}{\partial r_{\mathrm{D}}}+\lambda_{\mathrm{BD}}\right) d z_{\mathrm{D}}\right]=-h_{\mathrm{D}}\left|z_{\mathrm{D}}-z_{\mathrm{wD}}\right| \leq \frac{\varepsilon}{2}

相邻干扰井 ( j = 1 , 2 , 3 ) ( j = 1 , 2 , 3 ) (j=1,2,3cdots)(j=1,2,3 \cdots) 如下。

  初始条件
p fD ( r jD , 0 ) = p mD ( r jD , 0 ) = 0 . p fD r jD , 0 = p mD r jD , 0 = 0 . p_(fD)(r_(jD),0)=p_(mD)(r_(jD),0)=0.p_{\mathrm{fD}}\left(r_{\mathrm{jD}}, 0\right)=p_{\mathrm{mD}}\left(r_{\mathrm{jD}}, 0\right)=0 .
  外边界条件
p fD | | r jD = p mD | | r jD = 0 p fD | | r jD = p mD | | r jD = 0 p_(fD)||_(r_(jD)rarr oo)=p_(mD)||_(r_(jD)rarr oo)=0p_{\mathrm{fD}}| |_{r_{\mathrm{jD}} \rightarrow \infty}=p_{\mathrm{mD}}| |_{r_{\mathrm{jD}} \rightarrow \infty}=0
  内边界条件
r jD p fD r jD | r jD = 1 λ BD = q jD . r jD p fD r jD r jD = 1 λ BD = q jD . r_(jD)(delp_(fD))/(delr_(jD))|_(r_(jD)=1)-lambda_(BD)=-q_(jD).\left.r_{\mathrm{jD}} \frac{\partial p_{\mathrm{fD}}}{\partial r_{\mathrm{jD}}}\right|_{r_{\mathrm{jD}}=1}-\lambda_{\mathrm{BD}}=-q_{\mathrm{jD}} .

顶部和底部条件为
p fD z D | z D = 0 = p mD z D | z D = 0 = 0 , p fD z D | z D = h D = p mD z D | z D = h D = 0 . p fD z D z D = 0 = p mD z D z D = 0 = 0 , p fD z D z D = h D = p mD z D z D = h D = 0 . {:[(delp_(fD))/(delz_(D))|_(z_(D)=0)=(delp_(mD))/(delz_(D))|_(z_(D)=0)=0","],[(delp_(fD))/(delz_(D))|_(z_(D)=h_(D))=(delp_(mD))/(delz_(D))|_(z_(D)=h_(D))=0.]:}\begin{gathered} \left.\frac{\partial p_{\mathrm{fD}}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=0}=\left.\frac{\partial p_{\mathrm{mD}}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=0}=0, \\ \left.\frac{\partial p_{\mathrm{fD}}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=h_{\mathrm{D}}}=\left.\frac{\partial p_{\mathrm{mD}}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=h_{\mathrm{D}}}=0 . \end{gathered}

  3. 模型解答


为了消除上述无量纲方程中的二次项,采用了 Pedrosa 变量替换和常规扰动方法:
p fD = 1 γ D ln ( 1 γ D ξ D ) p fD = 1 γ D ln 1 γ D ξ D p_(fD)=-(1)/(gamma_(D))ln(1-gamma_(D)xi_(D))p_{\mathrm{fD}}=-\frac{1}{\gamma_{\mathrm{D}}} \ln \left(1-\gamma_{\mathrm{D}} \xi_{\mathrm{D}}\right)

因为 γ D 1 γ D 1 gamma_(D)≪1\gamma_{D} \ll 1 ,所以采用零阶扰动解,然后对其进行拉普拉斯变换。最后,方程(5)和(6)变为如下:

中央测试倾斜井:
1 r D r D ( r D ξ ¯ D 0 r D ) + λ BD s r D + 2 ξ ¯ D 0 z D 2 = ω s ξ ¯ D 0 + ( 1 ω ) s p ¯ mD , p ¯ mD = λ m λ m + ( 1 ω ) s ξ ¯ D 0 , ξ ¯ D 0 | r D = p ¯ mD | r D = 0 , lim ε 0 [ lim r D 0 z wD ε / 2 z wD ( ε / 2 ) r D e γ D p DD ( ξ ¯ D 0 r D + λ BD ) d z D ] = h D , | z D z wD | ε 2 . 1 r D r D r D ξ ¯ D 0 r D + λ BD s r D + 2 ξ ¯ D 0 z D 2 = ω s ξ ¯ D 0 + ( 1 ω ) s p ¯ mD , p ¯ mD = λ m λ m + ( 1 ω ) s ξ ¯ D 0 , ξ ¯ D 0 r D = p ¯ mD r D = 0 , lim ε 0 lim r D 0 z wD ε / 2 z wD ( ε / 2 ) r D e γ D p DD ξ ¯ D 0 r D + λ BD d z D = h D , z D z wD ε 2 . {:[(1)/(r_(D))(del)/(delr_(D))(r_(D)(del bar(xi)_(D0))/(delr_(D)))+(lambda_(BD))/(sr_(D))+(del^(2) bar(xi)_(D0))/(delz_(D)^(2))=omega s bar(xi)_(D0)+(1-omega)s bar(p)_(mD)","],[ bar(p)_(mD)=(lambda_(m))/(lambda_(m)+(1-omega)s) bar(xi)_(D0)","],[ bar(xi)_(D0)|_(r_(D)rarr oo)= bar(p)_(mD)|_(r_(D)rarr oo)=0","],[lim_(epsi rarr0)[lim_(r_(D)longrightarrow0)int_(z_(wD)-epsi//2)^(z_(wD)(epsi//2))r_(D)e^(-gamma_(D)p_(DD))((del bar(xi)_(D0))/(delr_(D))+lambda_(BD))dz_(D)]=-h_(D)","|z_(D)-z_(wD)| <= (epsi)/(2).]:}\begin{gathered} \frac{1}{r_{\mathrm{D}}} \frac{\partial}{\partial r_{\mathrm{D}}}\left(r_{\mathrm{D}} \frac{\partial \bar{\xi}_{\mathrm{D} 0}}{\partial r_{\mathrm{D}}}\right)+\frac{\lambda_{\mathrm{BD}}}{s r_{\mathrm{D}}}+\frac{\partial^{2} \bar{\xi}_{\mathrm{D} 0}}{\partial z_{\mathrm{D}}^{2}}=\omega s \bar{\xi}_{\mathrm{D} 0}+(1-\omega) s \bar{p}_{\mathrm{mD}}, \\ \bar{p}_{\mathrm{mD}}=\frac{\lambda_{\mathrm{m}}}{\lambda_{\mathrm{m}}+(1-\omega) s} \bar{\xi}_{\mathrm{D} 0}, \\ \left.\bar{\xi}_{\mathrm{D} 0}\right|_{r_{\mathrm{D}} \rightarrow \infty}=\left.\bar{p}_{\mathrm{mD}}\right|_{r_{\mathrm{D}} \rightarrow \infty}=0, \\ \lim _{\varepsilon \rightarrow 0}\left[\lim _{r_{\mathrm{D}} \longrightarrow 0} \int_{z_{\mathrm{wD}}-\varepsilon / 2}^{z_{\mathrm{wD}}(\varepsilon / 2)} r_{\mathrm{D}} e^{-\gamma_{\mathrm{D}} p_{\mathrm{DD}}}\left(\frac{\partial \bar{\xi}_{\mathrm{D} 0}}{\partial r_{\mathrm{D}}}+\lambda_{\mathrm{BD}}\right) d z_{\mathrm{D}}\right]=-h_{\mathrm{D}},\left|z_{\mathrm{D}}-z_{\mathrm{wD}}\right| \leq \frac{\varepsilon}{2} . \end{gathered}

相邻干扰井 ( j = 1 , 2 , 3 ) ( j = 1 , 2 , 3 ) (j=1,2,3cdots)(j=1,2,3 \cdots) :
1 r jD r jD ( r jD ξ ¯ D 0 j r jD ) + λ BD s r D = ω s ξ ¯ D 0 j + ( 1 ω ) s p ¯ mD 1 r jD r jD r jD ξ ¯ D 0 j r jD + λ BD s r D = ω s ξ ¯ D 0 j + ( 1 ω ) s p ¯ mD (1)/(r_(jD))(del)/(delr_(jD))(r_(jD)(del bar(xi)_(D0j))/(delr_(jD)))+(lambda_(BD))/(sr_(D))=omega s bar(xi)_(D0j)+(1-omega)s bar(p)_(mD)\frac{1}{r_{\mathrm{jD}}} \frac{\partial}{\partial r_{\mathrm{jD}}}\left(r_{\mathrm{jD}} \frac{\partial \bar{\xi}_{\mathrm{D} 0 \mathrm{j}}}{\partial r_{\mathrm{jD}}}\right)+\frac{\lambda_{\mathrm{BD}}}{s r_{\mathrm{D}}}=\omega s \bar{\xi}_{\mathrm{D} 0 \mathrm{j}}+(1-\omega) s \bar{p}_{\mathrm{mD}}
p ¯ mD = λ m λ m + ( 1 ω ) s ξ ¯ D 0 j ξ ¯ D 0 j | r jD = p ¯ mD | r iD = 0 , ξ ¯ D 0 j r jD + λ BD s = q jD s ξ ¯ D 0 z D | z D = 0 = p mD z D | z D = 0 = 0 , ξ ¯ D 0 z D | z D = h D = p mD z D | z D = h D = 0 p ¯ mD = λ m λ m + ( 1 ω ) s ξ ¯ D 0 j ξ ¯ D 0 j r jD = p ¯ mD r iD = 0 , ξ ¯ D 0 j r jD + λ BD s = q jD s ξ ¯ D 0 z D z D = 0 = p mD z D z D = 0 = 0 , ξ ¯ D 0 z D z D = h D = p mD z D z D = h D = 0 {:[ bar(p)_(mD)=(lambda_(m))/(lambda_(m)+(1-omega)s) bar(xi)_(D0j)],[ bar(xi)_(D0j)|_(r_(jD)rarr oo)= bar(p)_(mD)|_(r_(iD)rarr oo)=0","],[(del bar(xi)_(D0j))/(delr_(jD))+(lambda_(BD))/(s)=-(q_(jD))/(s)],[(del bar(xi)_(D0))/(delz_(D))|_(z_(D)=0)=(delp_(mD))/(delz_(D))|_(z_(D)=0)=0","],[(del bar(xi)_(D0))/(delz_(D))|_(z_(D)=h_(D))=(delp_(mD))/(delz_(D))|_(z_(D)=h_(D))=0]:}\begin{gathered} \bar{p}_{\mathrm{mD}}=\frac{\lambda_{\mathrm{m}}}{\lambda_{\mathrm{m}}+(1-\omega) s} \bar{\xi}_{\mathrm{D} 0 j} \\ \left.\bar{\xi}_{\mathrm{D} 0 \mathrm{j}}\right|_{r_{\mathrm{jD}} \rightarrow \infty}=\left.\bar{p}_{\mathrm{mD}}\right|_{r_{\mathrm{iD}} \rightarrow \infty}=0, \\ \frac{\partial \bar{\xi}_{\mathrm{D} 0 \mathrm{j}}}{\partial r_{\mathrm{jD}}}+\frac{\lambda_{\mathrm{BD}}}{s}=-\frac{q_{\mathrm{jD}}}{s} \\ \left.\frac{\partial \bar{\xi}_{\mathrm{D} 0}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=0}=\left.\frac{\partial p_{\mathrm{mD}}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=0}=0, \\ \left.\frac{\partial \bar{\xi}_{\mathrm{D} 0}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=h_{\mathrm{D}}}=\left.\frac{\partial p_{\mathrm{mD}}}{\partial z_{\mathrm{D}}}\right|_{z_{\mathrm{D}}=h_{\mathrm{D}}}=0 \end{gathered}

根据叠加原理,在拉普拉斯空间中对模型进行有限余弦积分变换,关于 z D z D z_(D)z_{\mathrm{D}} ,然后使用格林函数和零阶再生解,可以找到模型的点源解为
ξ ¯ D ξ ¯ D 0 + ξ ¯ D 0 j = 1 s K 0 [ r D s f ( s ) ] + 2 s n = 1 K 0 ( r D s f ( s ) + n 2 π 2 h D 2 ) × cos ( n π z D h D ) cos ( n π z wD h D ) + j = 1 m [ q Dj K 0 ( r jD s f ( s ) ) s s f ( s ) K 1 ( s f ( s ) ) + 0 + G ( r jD , τ ) d τ ] + 0 + G ( r D , τ ) d τ . ξ ¯ D ξ ¯ D 0 + ξ ¯ D 0 j = 1 s K 0 r D s f ( s ) + 2 s n = 1 K 0 r D s f ( s ) + n 2 π 2 h D 2 × cos n π z D h D cos n π z wD h D + j = 1 m q Dj K 0 r jD s f ( s ) s s f ( s ) K 1 ( s f ( s ) ) + 0 + G r jD , τ d τ + 0 + G r D , τ d τ . {:[ bar(xi)_(D)~~ bar(xi)_(D0)+ bar(xi)_(D0j)=(1)/(s)K_(0)[r_(D)sqrt(sf(s))]],[+(2)/(s)sum_(n=1)^(oo)K_(0)(r_(D)sqrt(sf(s)+(n^(2)pi^(2))/(h_(D)^(2))))xx cos((n piz_(D))/(h_(D)))cos((n piz_(wD))/(h_(D)))],[+sum_(j=1)^(m)[(q_(Dj)K_(0)(r_(jD)sqrt(sf(s))))/(ssqrt(sf(s))K_(1)(sqrt(sf(s))))+int_(0)^(+oo)G(r_(jD),tau)d tau]],[+int_(0)^(+oo)G(r_(D),tau)d tau.]:}\begin{aligned} \bar{\xi}_{\mathrm{D}} \approx & \bar{\xi}_{\mathrm{D} 0}+\bar{\xi}_{\mathrm{D} 0 \mathrm{j}}=\frac{1}{s} K_{0}\left[r_{\mathrm{D}} \sqrt{s f(s)}\right] \\ & +\frac{2}{s} \sum_{n=1}^{\infty} K_{0}\left(r_{\mathrm{D}} \sqrt{s f(s)+\frac{n^{2} \pi^{2}}{h_{\mathrm{D}}^{2}}}\right) \times \cos \left(\frac{n \pi z_{\mathrm{D}}}{h_{\mathrm{D}}}\right) \cos \left(\frac{n \pi z_{\mathrm{wD}}}{h_{\mathrm{D}}}\right) \\ & +\sum_{\mathrm{j}=1}^{m}\left[\frac{q_{\mathrm{Dj}} \mathrm{~K}_{0}\left(r_{\mathrm{jD}} \sqrt{s f(s)}\right)}{s \sqrt{s f(s)} \mathrm{K}_{1}(\sqrt{s f(s)})}+\int_{0}^{+\infty} G\left(r_{\mathrm{jD}}, \tau\right) d \tau\right] \\ & +\int_{0}^{+\infty} G\left(r_{\mathrm{D}}, \tau\right) d \tau . \end{aligned}

假设中央倾斜井的井筒内流动均匀受到干扰,根据点源理论,沿井筒方向对点源解进行积分,可以得到中央倾斜井在多井干扰下的底部压力解。
ξ ¯ wDN = 1 s L D L D / 2 L D / 2 K o [ r ¯ D s f ( s ) ] d η + 2 s L D L D / 2 L D / 2 n = 1 K 0 ( r ¯ D s f ( s ) + n 2 π 2 h D 2 ) × cos ( n π z D h D ) cos ( n π z ¯ wD h D ) d η + 1 L D 0 + L D / 2 L D / 2 G ( r ¯ D , τ ) d η d τ + j = 1 m [ q Dj K 0 ( r jD s f ( s ) ) s s f ( s ) K 1 ( s f ( s ) ) + 0 + G ( r jD , τ ) d τ ] ξ ¯ wDN = 1 s L D L D / 2 L D / 2 K o r ¯ D s f ( s ) d η + 2 s L D L D / 2 L D / 2 n = 1 K 0 r ¯ D s f ( s ) + n 2 π 2 h D 2 × cos n π z D h D cos n π z ¯ wD h D d η + 1 L D 0 + L D / 2 L D / 2 G r ¯ D , τ d η d τ + j = 1 m q Dj K 0 r jD s f ( s ) s s f ( s ) K 1 ( s f ( s ) ) + 0 + G r jD , τ d τ {:[ bar(xi)_(wDN)=(1)/(sL_(D))int_(-L_(D)//2)^(L_(D)//2)K_(o)[ bar(r)_(D)sf(s)]d eta],[+(2)/(sL_(D))int_(-L_(D)//2)^(L_(D)//2)sum_(n=1)^(oo)K_(0)( bar(r)_(D)sqrt(sf(s)+(n^(2)pi^(2))/(h_(D)^(2))))],[ xx cos((n piz_(D))/(h_(D)))cos((n pi bar(z)_(wD))/(h_(D)))d eta+(1)/(L_(D))int_(0)^(+oo)int_(-L_(D)//2)^(L_(D)//2)G( bar(r)_(D),tau)d eta d tau],[+sum_(j=1)^(m)[(q_(Dj)K_(0)(r_(jD)sqrt(sf(s))))/(ssqrt(sf(s))K_(1)(sqrt(sf(s))))+int_(0)^(+oo)G(r_(jD),tau)d tau]]:}\begin{aligned} \bar{\xi}_{\mathrm{wDN}}= & \frac{1}{s L_{\mathrm{D}}} \int_{-L_{\mathrm{D}} / 2}^{L_{\mathrm{D}} / 2} K_{\mathrm{o}}\left[\bar{r}_{\mathrm{D}} s f(s)\right] d \eta \\ & +\frac{2}{s L_{\mathrm{D}}} \int_{-L_{\mathrm{D}} / 2}^{L_{\mathrm{D}} / 2} \sum_{n=1}^{\infty} K_{0}\left(\bar{r}_{\mathrm{D}} \sqrt{s f(s)+\frac{n^{2} \pi^{2}}{h_{\mathrm{D}}^{2}}}\right) \\ & \times \cos \left(\frac{n \pi z_{\mathrm{D}}}{h_{\mathrm{D}}}\right) \cos \left(\frac{n \pi \bar{z}_{\mathrm{wD}}}{h_{\mathrm{D}}}\right) d \eta+\frac{1}{L_{\mathrm{D}}} \int_{0}^{+\infty} \int_{-L_{\mathrm{D}} / 2}^{L_{\mathrm{D}} / 2} G\left(\bar{r}_{\mathrm{D}}, \tau\right) d \eta d \tau \\ & +\sum_{\mathrm{j}=1}^{m}\left[\frac{q_{\mathrm{Dj}} \mathrm{~K}_{0}\left(r_{\mathrm{jD}} \sqrt{s f(s)}\right)}{s \sqrt{s f(s)} \mathrm{K}_{1}(\sqrt{s f(s)})}+\int_{0}^{+\infty} G\left(r_{\mathrm{jD}}, \tau\right) d \tau\right] \end{aligned}
r ¯ D = ( x D x wD η sin θ w ) 2 + ( y D y wD ) 2 z ¯ wD = z wD + η cos θ w r ¯ D = x D x wD η sin θ w 2 + y D y wD 2 z ¯ wD = z wD + η cos θ w {:[ bar(r)_(D)=sqrt((x_(D)-x_(wD)-eta sin theta_(w))^(2)+(y_(D)-y_(wD))^(2))],[ bar(z)_(wD)=z_(wD)+eta cos theta_(w)]:}\begin{gathered} \bar{r}_{\mathrm{D}}=\sqrt{\left(x_{\mathrm{D}}-x_{\mathrm{wD}}-\eta \sin \theta_{\mathrm{w}}\right)^{2}+\left(y_{\mathrm{D}}-y_{\mathrm{wD}}\right)^{2}} \\ \bar{z}_{\mathrm{wD}}=z_{\mathrm{wD}}+\eta \cos \theta_{\mathrm{w}} \end{gathered}
θ w = arctan ( K fv K fh tan θ ) , L D = L r w K fv K fh sin 2 θ + K fv K fh cos 2 θ . θ w = arctan K fv K fh tan θ , L D = L r w K fv K fh sin 2 θ + K fv K fh cos 2 θ . {:[theta_(w)=arctan(sqrt((K_(fv))/(K_(fh)))tan theta)","],[L_(D)=(L)/(r_(w))sqrt((K_(fv))/(K_(fh))sin^(2)theta+(K_(fv))/(K_(fh))cos^(2)theta).]:}\begin{gathered} \theta_{\mathrm{w}}=\arctan \left(\sqrt{\frac{K_{\mathrm{fv}}}{K_{\mathrm{fh}}}} \tan \theta\right), \\ L_{\mathrm{D}}=\frac{L}{r_{\mathrm{w}}} \sqrt{\frac{K_{\mathrm{fv}}}{K_{\mathrm{fh}}} \sin ^{2} \theta+\frac{K_{\mathrm{fv}}}{K_{\mathrm{fh}}} \cos ^{2} \theta} . \end{gathered}

杜哈梅尔原理和压力降叠加原理用于解决中央倾斜井在多井干扰下的井底压力解。
ξ ¯ WD = s ξ ¯ WDN + S s + C D s 2 ( s ξ ¯ WDN + S ) ξ ¯ WD = s ξ ¯ WDN + S s + C D s 2 s ξ ¯ WDN + S bar(xi)_(WD)=(s bar(xi)_(WDN)+S)/(s+C_(D)s^(2)(s bar(xi)_(WDN)+S))\bar{\xi}_{\mathrm{WD}}=\frac{s \bar{\xi}_{\mathrm{WDN}}+S}{s+C_{\mathrm{D}} s^{2}\left(s \bar{\xi}_{\mathrm{WDN}}+S\right)}

通过 Stehfest 数值反演,可以获得实际空间中的井底压力解[24, 25]:
p wD = 1 γ D ln ( 1 γ D ξ ¯ WD ) p wD = 1 γ D ln 1 γ D ξ ¯ WD p_(wD)=-(1)/(gamma_(D))ln(1-gamma_(D) bar(xi)_(WD))p_{\mathrm{wD}}=-\frac{1}{\gamma_{\mathrm{D}}} \ln \left(1-\gamma_{\mathrm{D}} \bar{\xi}_{\mathrm{WD}}\right)


4. 结果与讨论


4.1. 模型验证。当无量纲阈值压力梯度 ( λ BD ) λ BD (lambda_(BD))\left(\lambda_{\mathrm{BD}}\right) 、无量纲应力敏感系数 ( γ D ) γ D (gamma_(D))\left(\gamma_{D}\right) 和井倾角 ( θ ) ( θ ) (theta)(\theta) 取为零,并且没有邻井干扰时,该模型与双介质储层的常规垂直井压力降测试模型相同。为了验证本文中的模型,通过数值方法获得的类型曲线与图 4 中所示的通过解析解获得的常规双介质储层井测试曲线进行比较。

如图 4 所示,当无量纲阈值压力梯度 ( λ BD ) λ BD (lambda_(BD))\left(\lambda_{\mathrm{BD}}\right) 、无量纲应力敏感系数( γ D γ D gamma_(D)\gamma_{\mathrm{D}} )和井倾角 ( θ ) ( θ ) (theta)(\theta) 取为零,并且没有邻井干扰时,通过数值解和解析解获得的压力降测试的类型曲线是相同的。因此,本文中的数值解是可靠的。


4.2. 类型曲线。通过数值反演获得实际空间中的压力解。以中央倾斜井周围存在三个相邻井为例,其压力曲线和压力导数曲线如图 5-12 所示;流动区域被识别。

区域 1。早期井筒储存期:压力曲线和压力导数曲线是斜率等于“1”的线段,反映了井筒储层效应的影响。

区域 2。皮肤效应期:这一时期的主要影响因素是皮肤因子效应,压力导数曲线呈现“隆起”形状。

区域 3。井倾角影响期:压力导数曲线呈现“凹形”,反映了井倾角的影响。当井倾角较大时,显示出特征


在水平井中垂直径向流动的特征;当井倾角较小时,表现出垂直井的径向流动特征

区域 4。裂缝径向流动期:在此期间,主要影响因素是流体从裂缝流向井底,压力导数曲线是一条水平直线。

区域 5。泄漏期:此期间的主要影响因素是基质到裂缝的泄漏能力,压力导数曲线呈现“凹形”。

区域 6。第一次径向流动期:这一时期的主要影响因素是中央倾斜井的晚期径向流动,压力导数曲线是一条水平直线。

区域 7. 过渡流动期:径向流动之间的过渡期

区域 8。第二个径向流动周期:它反映了受到最近邻井干扰的中央倾斜井的径向流动特征,压力导数曲线是一条水平直线。

区域 9. 过渡流动期:径向流动之间的过渡期

区域 10。第三个径向流动周期:它反映了中央倾斜井受到两个最近邻井干扰的径向流动特征,压力导数曲线是一条水平直线。

区域 11. 过渡流动期:径向流动之间的过渡期

区域 12。第四个径向流动周期:它反映了受到三个最近邻井干扰的中央倾斜井的径向流动特征,压力导数曲线是一条水平直线。


4.3. 敏感性分析


4.3.1. 渗透率模量。如图 6 所示,随着渗透率模量的增加,压力动态曲线在第 7 区后逐渐上升,且渗透率模量越大,上升的程度也越大。上述现象的原因是,渗透率模量越大,随着压力的增加,渗透率的下降幅度也越大;因此,流体流动的阻力变大,流动所需的压力也更大。当渗透率模量较大时,储层渗透率在渗流后期的高压下急剧下降,压力动态曲线反映出类似于封闭外边界的特征。


4.3.2. 阈值压力梯度。如图 7 所示,重油的阈值压力梯度越高,压力动态曲线的向上翘曲程度越大,并且在区域 7 之后向上翘曲明显。这是因为重油具有阈值压力梯度,重油在孔隙中的流动能力减弱,重油流动所需的位移压力差也更大。


图 4:数值解与压力降井测试的解析解比较。


图 5:典型压力动态曲线。


4.3.3. 邻近井的生产。如图 8 所示,以中央倾斜井周围存在三口相连的邻近井为例,中央倾斜井的压力动态曲线有四个晚期径向流动阶段, N th ( N > 2 ) N th  ( N > 2 ) N^("th ")(N > 2)N^{\text {th }}(N>2) 径向流动值与第一个径向流动值的比值等于 1 + j = 1 k q j D 1 + j = 1 k q j D 1+sum_(j=1)^(k)q_(jD)1+\sum_{j=1}^{k} q_{j D} k k kk 是对 N N NN th th  ^("th "){ }^{\text {th }} 径向流动阶段有影响的邻近井数量)。以 q 1 D = q 2 D = q 3 D = 1 q 1 D = q 2 D = q 3 D = 1 q_(1D)=q_(2D)=q_(3D)=1q_{1 \mathrm{D}}=q_{2 \mathrm{D}}=q_{3 \mathrm{D}}=1 为例,第三个径向流动阶段是测试井与两口邻近井共同影响的结果,因此第三个径向流动值为 0.5 ( 1 + 2 ) = 1.5 0.5 ( 1 + 2 ) = 1.5 0.5**(1+2)=1.50.5 *(1+2)=1.5


4.3.4. 井倾角。如图 9 所示,井倾角越大,井筒压力降越大,压力动态曲线越低。由于倾斜井完全穿透地层,井倾角越大,井筒长度越长。在假设井筒内流量相等的情况下,井筒长度越长 L L LL ,井筒内的压力降越大。而井倾角主要影响井附近的流动,当井倾角较大 ( θ > 40 ) θ > 40 (theta > 40^(@))\left(\theta>40^{\circ}\right) 时,曲线将显示出类似于水平井的早期垂直径向流动。


图 6:渗透率模量对压力动态曲线的影响。


图 7:阈值压力梯度对压力动态曲线的影响。


当井倾角相对较小时 ( θ < 40 ) θ < 40 (theta < 40^(@))\left(\theta<40^{\circ}\right) ,垂直径向流动的特征消失,这类似于垂直井的径向流动。


4.3.5. 邻近井距离。如图 10 所示,随着邻近井与中央倾斜井之间距离的增加,第二次径向流动的开始时间是


图 8:相邻井的产量对压力动态曲线的影响。


图 9:井倾角对压力动态曲线的影响。


图 10:距离对压力动态曲线的影响。


图 11:井位示意图。


延迟,换班持续时间变得更短。当距离较短时,第二个径向流动段逐渐消失并融入第三个径向流动段。

  5. 现场应用


5.1. 现场测试数据匹配方法。为了减少非唯一性,在将模型应用于现场测试数据的井测试解释时,应遵循以下步骤:


(1) 使用实际的压力恢复数据,绘制压力曲线和压力导数曲线。判断在压力响应曲线的后期是否存在多个径向流动阶段,如果它


符合相邻井的干扰井测试曲线,采用相邻井干扰模型进行参数解释;否则,采用普通单井模型


(2) 根据所选的井测试解释模型(不考虑应力敏感性和阈值压力梯度的影响),通过改变井筒储存系数、皮肤系数、孔隙间流动系数和弹性储存比来匹配初始压力响应曲线


(3)输入相邻井之间的距离和相邻井的产量,以匹配后期的多级径向流动段


(4) 常规油藏中相邻井的解释结果参数作为初始参数,输入到考虑应力敏感性和阈值压力梯度的双介质油藏相邻井干扰井测试模型中。通过改变应力敏感性系数和启动压力梯度,进一步匹配和解释压力响应曲线。最终,获得相关参数。


5.2. 示例说明。SZ 油藏位于中国渤海湾。它是一个重油油藏,平均油粘度为 320 mPa s 320 mPa s 320mPa*s320 \mathrm{mPa} \cdot \mathrm{s} 。X2 井是 SZ 油藏中的一口倾斜井,倾斜角度为 73 73 73^(@)73^{\circ} ,油藏的有效厚度为 31.5 米,孔隙度为 13.5 % 13.5 % 13.5%13.5 \% ,体积系数为 1.06。通过油藏砂


图 12:压力匹配曲线。

表 1:所研究储层的基本参数。
  井名   距离 (m)
每日液体产量 ( m 3 / d ) m 3 / d (m^(3)//d)\left(\mathrm{m}^{3} / \mathrm{d}\right)
X1 173 200
X3 162 150
X4 218 200
Well name Distance (m) Daily liquid production (m^(3)//d) X1 173 200 X3 162 150 X4 218 200| Well name | Distance (m) | Daily liquid production $\left(\mathrm{m}^{3} / \mathrm{d}\right)$ | | :--- | :---: | :---: | | X1 | 173 | 200 | | X3 | 162 | 150 | | X4 | 218 | 200 |

表 2:解释结果。
  参数   

井筒储存系数
0.17 m 3 / MPa 0.17 m 3 / MPa 0.17m^(3)//MPa0.17 \mathrm{~m}^{3} / \mathrm{MPa}
  皮肤系数 0.09

孔隙间流动系数
5.7 × 10 7 5.7 × 10 7 5.7 xx10^(-7)5.7 \times 10^{-7}
  渗透率 1420 mD
  储存比 5.4 × 10 2 5.4 × 10 2 5.4 xx10^(-2)5.4 \times 10^{-2}

阈值压力梯度
2.4 × 10 4 MPa / m 2.4 × 10 4 MPa / m 2.4 xx10^(-4)MPa//m2.4 \times 10^{-4} \mathrm{MPa} / \mathrm{m}
  渗透率模量 1.5 × 10 2 MPa 1 1.5 × 10 2 MPa 1 1.5 xx10^(-2)MPa^(-1)1.5 \times 10^{-2} \mathrm{MPa}^{-1}

X1 的干扰量
72 m 3 / d 72 m 3 / d 72m^(3)//d72 \mathrm{~m}^{3} / \mathrm{d}

X3 的干扰量
10 m 3 / d 10 m 3 / d 10m^(3)//d10 \mathrm{~m}^{3} / \mathrm{d}

X4 的干扰量
0 m 3 / d 0 m 3 / d 0m^(3)//d0 \mathrm{~m}^{3} / \mathrm{d}
Parameters Value Wellbore storage coefficient 0.17m^(3)//MPa Skin factor 0.09 Interporosity flow coefficient 5.7 xx10^(-7) Permeability 1420 mD Storativity ratio 5.4 xx10^(-2) Threshold pressure gradient 2.4 xx10^(-4)MPa//m Permeability modulus 1.5 xx10^(-2)MPa^(-1) Interference quantity of X1 72m^(3)//d Interference quantity of X3 10m^(3)//d Interference quantity of X4 0m^(3)//d| Parameters | Value | | :--- | :---: | | Wellbore storage coefficient | $0.17 \mathrm{~m}^{3} / \mathrm{MPa}$ | | Skin factor | 0.09 | | Interporosity flow coefficient | $5.7 \times 10^{-7}$ | | Permeability | 1420 mD | | Storativity ratio | $5.4 \times 10^{-2}$ | | Threshold pressure gradient | $2.4 \times 10^{-4} \mathrm{MPa} / \mathrm{m}$ | | Permeability modulus | $1.5 \times 10^{-2} \mathrm{MPa}^{-1}$ | | Interference quantity of X1 | $72 \mathrm{~m}^{3} / \mathrm{d}$ | | Interference quantity of X3 | $10 \mathrm{~m}^{3} / \mathrm{d}$ | | Interference quantity of X4 | $0 \mathrm{~m}^{3} / \mathrm{d}$ |

连通性分析,三个相邻的井与其(X1/X3/X4)具有良好的连通性,如图 11 所示。当进行 X2 压力建立测试时,相邻的井(X1/X3/X4)也同时被封闭,X2 周围三个干扰井的位置和参数如表 1 所示。本文提出的模型被选择用于匹配压力动态曲线,如图 12 所示。该模型的匹配精度良好,特别是在压力动态曲线的后期阶段。表 2 显示了解释结果。

  6. 结论


(1) 针对海上重油储层在多井干扰下的倾斜井测试问题,本文综合考虑了储层的渗透率应力敏感性和重油的阈值压力梯度,建立了多井干扰下双介质重油储层的倾斜井测试模型。并利用格林函数和叠加原理,得到了拉普拉斯空间下的解析解。最后,通过数值反演绘制了测试倾斜井的压力动态曲线,并以三个相邻井为例。


(2) 本文分析了渗透率模量、阈值压力梯度、多井生产、井倾角和井距对压力动态曲线的影响规律。在相邻井的影响下,测试倾斜井后期的压力导数曲线向上翘起,产生多个“平台”。通过敏感性分析,倾斜井存在一个临界井倾角 40 40 40^(@)40^{\circ} 。当井倾角大于 40 40 40^(@)40^{\circ} 时,将出现类似于水平井的垂直径向流动。


(3) 新模型应用于渤海湾 SZ 油田的井测试解释,获得了高匹配精度。该模型不仅为类似类型储层的多井干扰井测试提供了理论指导,还为井间连通性的定量表征提供了基础。

  术语表

p f p f p_(f)p_{\mathrm{f}} :
裂缝系统压力 ( MPa )
v fr v fr  v_("fr ")v_{\text {fr }} :
裂缝中流体的径向速度 ( m / s m / s m//s\mathrm{m} / \mathrm{s} )
K fh K fh  K_("fh ")K_{\text {fh }} :
裂缝水平渗透率 ( mD )
v fz v fz v_(fz)v_{\mathrm{fz}} :
裂缝中流体的垂直速度 ( m / s m / s m//s\mathrm{m} / \mathrm{s} )
μ μ mu\mu :
油的表观粘度 ( mPa s mPa s mPa*s\mathrm{mPa} \cdot \mathrm{s} )
K fv K fv  K_("fv ")K_{\text {fv }} :
裂缝垂直渗透率
λ B λ B lambda_(B)\lambda_{B} :
重油的阈值压力梯度 (MPa / m)
p i p i p_(i)p_{\mathrm{i}} :
初始地层压力 (MPa)
r r rr :   距离 (m)
ϕ f ϕ f phi_(f)\phi_{\mathrm{f}} :
裂缝孔隙度(无量纲)
C tf C tf  C_("tf ")C_{\text {tf }} :
裂缝压缩性 ( MPa 1 MPa 1 MPa^(-1)\mathrm{MPa}^{-1} )
C tm C tm  C_("tm ")C_{\text {tm }} :
矩阵可压缩性 ( MPa 1 MPa 1 MPa^(-1)\mathrm{MPa}^{-1} )
V f V f V_(f)V_{\mathrm{f}} :
裂缝体积比(无量纲)
K m K m K_(m)K_{\mathrm{m}} :   矩阵渗透率 (mD)
p m p m p_(m)p_{\mathrm{m}} :
矩阵系统压力 (MPa)
α α alpha\alpha :
矩阵块形状因子(无量纲)
t t tt :   时间 (s)
z z zz :   距离 (m)
q j D q j D q_(jD)q_{j \mathrm{D}} :
相邻井的无量纲产量(无量纲)
r jD r jD  r_("jD ")r_{\text {jD }} :
相邻井距离(无量纲)
ξ D ξ D xi_(D)\xi_{\mathrm{D}} :
扰动变形函数
s s ss :   拉普拉斯因子
ξ D 0 ξ D 0 xi_(D0)\xi_{\mathrm{D} 0} :
中央倾斜井的零阶扰动解
ξ D 0 : ξ 0  :  xi_("D "0" : ")\xi_{\text {D } 0 \text { : }} :
相邻井的零阶扰动解
m : m : m:m:
邻近井造成干扰的数量
η η eta\eta :   井筒积分因子
θ θ theta\theta :   井倾角
L L LL :
储层中倾斜井的长度
x D , y D , z D x D , y D , z D x_(D),y_(D),z_(D)x_{\mathrm{D}}, y_{\mathrm{D}}, z_{\mathrm{D}} :   坐标
S:   皮肤系数。
p_(f) : Fracture system pressure ( MPa ) v_("fr ") : Radial velocity of fluid in fracture ( m//s ) K_("fh ") : Fracture horizontal permeability ( mD ) v_(fz) : Vertical velocity of fluid in fracture ( m//s ) mu : Oil apparent viscosity ( mPa*s ) K_("fv ") : Fracture vertical permeability lambda_(B) : Threshold pressure gradient of heavy oil ( MPa / m) p_(i) : Initial formation pressure (MPa) r : Distance (m) phi_(f) : Fracture porosity (dimensionless) C_("tf ") : Fracture compressibility ( MPa^(-1) ) C_("tm ") : Matrix compressibility ( MPa^(-1) ) V_(f) : Fracture volume ratio (dimensionless) K_(m) : Matrix permeability (mD) p_(m) : Matrix system pressure (MPa) alpha : Matrix block shape factor (dimensionless) t : Time (s) z : Distance (m) q_(jD) : Dimensionless production of adjacent wells (dimensionless) r_("jD ") : Adjacent well distance (dimensionless) xi_(D) : Perturbation deformation function s : Laplace factor xi_(D0) : Zero order perturbation solution of central inclined well xi_("D "0" : ") : Zero order perturbation solution of adjacent well m: Number of adjacent wells causing interference eta : Wellbore integral factor theta : Well inclination angle L : Length of inclined well in reservoir x_(D),y_(D),z_(D) : Coordinate S: Skin factor.| $p_{\mathrm{f}}$ : | Fracture system pressure ( MPa ) | | :---: | :---: | | $v_{\text {fr }}$ : | Radial velocity of fluid in fracture ( $\mathrm{m} / \mathrm{s}$ ) | | $K_{\text {fh }}$ : | Fracture horizontal permeability ( mD ) | | $v_{\mathrm{fz}}$ : | Vertical velocity of fluid in fracture ( $\mathrm{m} / \mathrm{s}$ ) | | $\mu$ : | Oil apparent viscosity ( $\mathrm{mPa} \cdot \mathrm{s}$ ) | | $K_{\text {fv }}$ : | Fracture vertical permeability | | $\lambda_{B}$ : | Threshold pressure gradient of heavy oil ( MPa / m) | | $p_{\mathrm{i}}$ : | Initial formation pressure (MPa) | | $r$ : | Distance (m) | | $\phi_{\mathrm{f}}$ : | Fracture porosity (dimensionless) | | $C_{\text {tf }}$ : | Fracture compressibility ( $\mathrm{MPa}^{-1}$ ) | | $C_{\text {tm }}$ : | Matrix compressibility ( $\mathrm{MPa}^{-1}$ ) | | $V_{\mathrm{f}}$ : | Fracture volume ratio (dimensionless) | | $K_{\mathrm{m}}$ : | Matrix permeability (mD) | | $p_{\mathrm{m}}$ : | Matrix system pressure (MPa) | | $\alpha$ : | Matrix block shape factor (dimensionless) | | $t$ : | Time (s) | | $z$ : | Distance (m) | | $q_{j \mathrm{D}}$ : | Dimensionless production of adjacent wells (dimensionless) | | $r_{\text {jD }}$ : | Adjacent well distance (dimensionless) | | $\xi_{\mathrm{D}}$ : | Perturbation deformation function | | $s$ : | Laplace factor | | $\xi_{\mathrm{D} 0}$ : | Zero order perturbation solution of central inclined well | | $\xi_{\text {D } 0 \text { : }}$ : | Zero order perturbation solution of adjacent well | | $m:$ | Number of adjacent wells causing interference | | $\eta$ : | Wellbore integral factor | | $\theta$ : | Well inclination angle | | $L$ : | Length of inclined well in reservoir | | $x_{\mathrm{D}}, y_{\mathrm{D}}, z_{\mathrm{D}}$ : | Coordinate | | S: | Skin factor. |

  数据可用性


本研究所支持的发现的数据包含在文章中。

  利益冲突


作者声明他们在本论文的发表上没有利益冲突。

  致谢


本工作得到了中国海洋石油有限公司天津分公司的支持(YXKY-2018-TJ-04)。

  参考文献


[1] W. Weilin 和 A. Hongliang,“裂缝水平井的干扰压力分析模型,”Fault-Block Oil & Gas Field,第 26 卷,第 4 期,页 501-505,2019 年。


[2] W. Hongfeng, L. Xiaoping, W. Xiaopei, L. Xiuyu, Z. Songbai, 和 N. Yanbo, “多井干扰测试技术在克深气田勘探与开发中的应用”


领域,”石油地质与采收效率,第 25 卷,第 1 期,页码 100-105,2018 年。


[3] S. Hedong, “相邻井干扰下多井系统的压力建立分析,” 天然气工业, 第 36 卷, 第 5 期, 第 62-68 页, 2016 年。


[4] M. Onur, K. V. Serra, 和 A. C. Reynolds, “多井系统中井的压力建立数据分析,” SPE Formation Evaluation, 第 6 卷, 第 1 期, 页 101-110, 1991.


[5] T. Marhaendrajana, N. J. Kaczorowski, 和 T. A. Blasingame,“阿伦油田井测试性能的分析与解释,”发表于 1999 年休斯顿德克萨斯州的 SPE 年技术会议和展览。


[6] E. S. Adewole,“在双边水驱动的油藏中水平和垂直井组合的干扰测试分析”,发表于 2012 年尼日利亚年度国际会议和展览,尼日利亚拉各斯。


[7] Q. Deng, R.-s. Nie, Y.-l. Jia, X.-l. Wang, Y.-y. Chen, 和 Y. Xiong, “一种多井储层中井的压力建立分析新方法,” 发表在 2015 年埃及开罗 SPE 北非技术会议和展览上。


[8] Y. He, S. Cheng, J. Qin 等,“多段水平井的分析干扰测试分析,”《石油科学与工程杂志》,第 171 卷,页码 919-927,2018 年。


[9] J. Yang, L. I. Qi, F. Chen, Y. Liu, J. Zhang, 和 D. Gao, “裂缝-孔隙碳酸盐岩储层干扰测试研究,” Atlantis Press, vol. 2017, pp. 510-514, 2017.


[10] A. Kumar, P. Seth, K. Shrivastava, R. Manchanda, 和 M. M. Sharma, “综合分析示踪剂和压力干扰测试以识别井干扰,” SPE Journal, 第 25 卷, 第 4 期, 页 1623-1635, 2020.


[11] G. Han, Y. Liu, W. Liu, 和 D. Gao, “关于通过大裂缝连接的井的干扰测试的研究,” 应用科学, 第 9 卷, 第 1 期, 页 206, 2019.


[12] C. Shiqing, L. Meng, H. Youwei 等,“一种多井干扰压力瞬态分析方法用于确定多裂缝水平井的水源方向,”《中国石油大学学报(自然科学版)》,第 42 卷,第 5 期,页码 81-88,2018 年。


[13] H. Cinco, F. G. Miller, 和 H. J. Ramey Jr.,“由定向钻井井产生的非稳态压力分布,”《石油技术杂志》,第 27 卷,第 11 期,页 1392-1400,1975 年。


[14] L. Zhang, M. B. Dusseault, 和 J. A. Franklin,具有渗透率各向异性的介质中的倾斜井生产,1993。


[15] B. R. Sousa, “倾斜井测试分析,”坎皮纳大学(葡萄牙语),2012。


[16] H. Li, Q. Zhang, K. Wei, Y. Zeng, 和 Y. Zhu, “考虑非达西流的低渗透复合气藏倾斜井的井测试分析,” Energies, vol. 15, no. 5, p. 1654, 2022.


[17] Z. Xl, L. Xp, 和 X. Zhang, “斜井测试模型在正交断层层状双孔隙储层中的研究,” 井测试, 第 21 卷, 第 6 期, 2012 年。


[18] J. Ruizhong, S. Zeyang, C. Yongzheng, Z. Fulei, Z. Chunguang, 和 Y. Jianwei, “双介质低渗透储层中倾斜井的动态特性,” 岩性储层, 第 30 卷, 第 6 期, 第 131-137 页, 2018 年。


[19] Y. L. Jia, G. F. Sun, R. S. Nie, J. M. Li, 和 H. K. Li, “四重介质储层的流动模型和井测试曲线,” 岩性储层, 第 28 卷, 第 1 期, 第 123-127 页, 2016 年。


[20] X. Youjie, L. Qiguo, W. Rui, 和 L. Yicheng, “复合储层中具有复杂裂缝分布的裂缝水平井的压力瞬态,” 岩性储层, 第 31 卷, 第 5 期, 页 161-168, 2019.


[21] C. Zhongliang, W. Nutao, C. Huanghui, J. Tao, 和 D. Sen, “关于多层混合和倾斜井的井测试解释,” Well Testing, 第 26 卷, 第 3 期, 页 29-32+76, 2017.


[22] G. Li, X. Xin, G. Yu 等,“沥青质对多孔介质中重油阈值压力梯度的影响,”《先进材料科学》,第 13 卷,第 2 期,页 273-279,2021。


[23] D. Zhang, J. Peng, Y. Gu, 和 Y. Leng, “重油储层阈值压力梯度的实验研究,” 新疆石油地质, 第 33 卷, 第 2 期, 页 201-204, 2012.


[24] H. Stehfest, “算法 368:拉普拉斯变换的数值反演 [D5],” ACM 通讯, 第 13 卷, 第 1 期, 第 47-49 页, 1970 年。


[25] H. Stehfest,“关于算法 368 的备注:拉普拉斯变换的数值反演,”ACM 通讯,第 13 卷,第 10 期,第 624 页,1970 年。