致密油藏存在低孔低渗、储层非均质性、各向异性等特征^[1-2], 其储层改造区域(SRV)的渗流规律研究也遇到众多挑战。研究致密油藏SRV的渗流规律, 首先需要考虑由基质孔隙(微米、纳米级)-天然裂缝(微米、毫米级)-水力裂缝(毫米级)组成的空间多尺度体系^[3-7], 该体系渗透率、孔隙度差异大, 因此, 流体在不同介质中的流速差异大, 而且不同介质之间存在窜流。其次, 空间多尺度体系导致了流动的时间滞后效应和窜流的不平衡效应。其中时间滞后是指因孔隙、裂缝介质迂曲, 流体在介质中的运移时间较长, 因此某时刻的流动压力会对后续的流动产生影响, 但产生的影响会随着流动的进行而衰减。通过时间滞后效应^[8]表征SRV渗流可更好地修正等效连续介质模型, 同时可以有效考虑孔隙、裂缝几何形态对井底流动的影响。而不平衡效应指的是因不同介质流动截面的几何形态、空间尺度不同, 介质交界面前后的流速差异引起的窜流现象, 这样的窜流不仅受介质间压差的影响, 还受交界面附近的储容^{[5, 9]}流体影响。此外, 空间多尺度体系导致了流动的时间多尺度现象, 该体系的流动压力主要产生于基质孔隙、天然裂缝、人工裂缝等不同空间尺度的介质之中。随着生产的进行, 每种介质产生的流动压力降低的快慢相差较大, 这就体现了时间多尺度性。目前, 致密油藏SRV的渗流规律研究尚未考虑流体流动的时间多尺度现象^[10-11]。笔者研究SRV渗流时间多尺度问题, 建立考虑时间滞后的致密油藏SRV渗流模型, 分析不同渗透率、孔隙度分布下的流动时间尺度, 描述介质的非连续性、储层的非均质性、窜流的不平衡效应对井底流动压力的影响规律。

1 物理模型

图 1(a)为SRV流动的物理模型, 由SRV、水力裂缝2个区域组成。除水力裂缝、井筒交点外, 其余边界均封闭, 流动为单相流。

将SRV中天然裂缝和基质等效为非均质连续介质。其中, SRV渗透率k_m(m²)在水力裂缝、井筒交点处为k₀(m²), 沿x方向线性^[12-13]递减至k_x^* (m²), 沿y方向线性递减至k_y^*(m²):

$ {k_{\rm{m}}}\left( {x,y} \right) = {k_0}\left( {1 - \left( {1 - \frac{{k_x^ * }}{{{k_0}}}} \right)\frac{x}{{{l_{\rm{f}}}}}} \right)\left( {1 - \left( {1 - \frac{{k_y^ * }}{{{k_0}}}} \right)\frac{y}{{{l_{\rm{f}}}}}} \right), $

(1)

孔隙度φ_m在水力裂缝、井筒交点处为φ₀, 沿x方向线性递减至φ_x^*, 沿y方向线性递减至φ_y^*:

$ {\varphi _{\rm{m}}}\left( {x,y} \right) = {\varphi _0}\left( {1 - \left( {1 - \frac{{\varphi _x^ * }}{{{\varphi _0}}}} \right)\frac{x}{{{l_{\rm{f}}}}}} \right)\left( {1 - \left( {1 - \frac{{\varphi _y^ * }}{{{\varphi _0}}}} \right)\frac{y}{{{l_{\rm{f}}}}}} \right). $

(2)

式中, x为井筒方向, m; y为射孔方向, m; z为SRV厚度方向, m; l_f为水力裂缝半长, m。

水力裂缝与SRV中的流速差异极大, 两者交界面附近的储容流体也流向水力裂缝, 因此水力裂缝与SRV之间的窜流存在不平衡效应(图 1(b))。另外, 储层中的孔隙迂曲, 天然裂缝、人工裂缝也迂曲, 且不断分叉, 流动截面随空间变化, 流动过程还存在着时间滞后效应(图 1(c))。

2 数学模型 2.1 单条孔隙微观流动模型

SRV中单条孔隙流动微观上符合Navier-Stokes方程, 流体沿孔隙方向流动, 流动方程为

$ \frac{{\partial \mathit{\boldsymbol{v}}}}{{\partial t}} - \frac{\mu }{\rho }\nabla _n^2\mathit{\boldsymbol{v}} = - \frac{1}{\rho }\nabla p, $

(3)

边界条件为

$ \mathit{\boldsymbol{v}}\left| {_{\partial \sigma }} \right. = 0, $

(4)

$ \int\limits_\sigma {\mathit{\boldsymbol{v}}{\rm{d}}\sigma } = Q\left( t \right), $

(5)

初始条件为

$ \mathit{\boldsymbol{v}}\left| {_{t = 0}} \right. = 0. $

(6)

式中, t为时间, s; v为流速, m/s; ρ为流体密度, kg/m³; p为流体产生的压力, Pa, 即流动压力; μ为流体黏度, Pa·s。沿流动方向取自然坐标, ∇_n为自然坐标下的梯度算子, ∇为真实空间直角坐标下的梯度算子, x_n为自然坐标, 垂直于x_n方向的孔隙截面为σ。如图 2(a)所示。

孔隙、天然裂缝迂曲, 截面几何形态随空间变化(图 2(b))。流动压力不仅包括克服流动阻力的压力, 还包括因流动截面、流速方向改变, 孔隙壁对流体产生的附加压力, 因此要想确定流动压力梯度与流速分布之间的关系, 应以通过流动截面σ的流量Q(m³/s)作为已知条件, 而如果以流动压力梯度作为已知条件, 得到的关系不够准确。

由文献[14], t时刻的流量受时间段[0, t]内所有u(0≤u≤t)时刻的流动压力梯度的影响:

$ \frac{{\partial \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}{{\partial t}} - \frac{\mu }{\rho }\nabla _n^2\mathit{\boldsymbol{ \boldsymbol{\varTheta} }} = 0, $

(7)

式中, Θ为u时刻的流动压力梯度对后续时刻流动的影响。边界条件为

$ \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left| {_{\partial \sigma }} \right. = 0, $

(8)

初始条件为

$ \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left| {_{t = 0}} \right. = \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}. $

(9)

式中, Ψ表示u时刻的流动压力梯度对u时刻流速分布的影响:

$ - \frac{\mu }{\rho }\nabla _n^2\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = - \frac{1}{\rho }\nabla p\left| {_{t = u}} \right., $

(10)

边界条件为

$ \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}\left| {_{\partial \sigma }} \right. = 0. $

(11)

根据文献[14]中的方法可得出流量Q与流动压力梯度之间的关系:

$ \int\limits_0^t {H'\left( {t - u} \right)\left| {\nabla p} \right|{\rm{d}}u} = Q\left( t \right). $

(12)

其中

$ H = \int\limits_\sigma {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}/\left| {\nabla p} \right|{\rm{d}}\sigma } . $

(13)

为得到式(7)的解Θ, 将方程(7)、(10)中的自然坐标变换为真实空间直角坐标:

$ \frac{{\partial \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}{{\partial t}} - \frac{\mu }{{\rho J}}{\nabla ^2}\mathit{\boldsymbol{ \boldsymbol{\varTheta} }} = 0, $

(14)

$ \frac{\mu }{{\rho J}}{\nabla ^2}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = \frac{1}{\rho }\nabla p\left| {_{t = u}} \right.. $

(15)

式中, J为自然坐标到真实空间直角坐标的坐标变换系数。孔隙越迂曲, J越大。Θ表示为

$ \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} = \sum\limits_{j = 1}^m {{f_j}\left( {x,y,z} \right){{\rm{e}}^{ - \frac{t}{{{\tau _j}}}}}} . $

(16)

式中, τ_j为孔隙、裂缝迂曲引起的滞后时间, s, 表示u时刻流动压力梯度对流动的影响衰减所需的时间。τ_j随着流动阻力的增加而减小, 但随着孔隙、裂缝迂曲程度的增加而增加。t=0时, 将式(16)代入(13)中得

$ H\left( 0 \right) = \int\limits_\sigma {\sum\limits_{j = 1}^m {{\mathit{\boldsymbol{f}}_j}\left( {x,y,z} \right)} /\left| {\nabla p} \right|{\rm{d}}\sigma } . $

(17)

H(0)表示流动截面的几何形态。根据Poiseuille定律, 有

$ H\left( 0 \right) = \frac{{{\rm{ \mathsf{ π} }}{r^4}}}{{8\mu }}. $

(18)

式中, r为孔隙半径, m。

滞后时间越小, 流动压力梯度的影响消失得越迅速, 流动也就越接近于达西流动, 较小的滞后时间无法很好地描述孔隙、裂缝的微观形态对油藏流动的影响。本文中忽略了较小的滞后时间, 仅考虑最大的滞后时间, 滞后时间个数m=1, 式(12)表示为

$ Q\left( t \right) = - H\left( 0 \right)\int\limits_0^t {\frac{{{{\rm{e}}^{ - \frac{{t - u}}{\tau }}}}}{\tau }\left| {\nabla p} \right|{\rm{d}}u} . $

(19)

2.2 对微观流动的尺度升级

在SRV中, 取长dx, 宽dy, 高h的单元体, 以x方向流量为例进行分析:

$ h{\mathit{\boldsymbol{v}}_x}{\rm{d}}y = hQ\left( t \right){n_\varphi }{\rm{d}}y = - hH\left( 0 \right){n_\varphi }\left( {\int\limits_0^t {\frac{{{{\rm{e}}^{ - \frac{{t - u}}{\tau }}}}}{\tau }\left| {\nabla p} \right|{\rm{d}}u} } \right){\rm{d}}y. $

(20)

式中, n_φ为单位面积上的孔隙数, m^-2。可以看出, 通过流动的时间滞后效应, 可对微观上迂曲的孔隙、裂缝中流动进行尺度升级。SRV渗透率k_m可表示为

$ {k_{\rm{m}}} = \mu H\left( 0 \right){n_\varphi }. $

(21)

单元体流动连续性方程为

$ \nabla \rho \mathit{\boldsymbol{v}} + \frac{{\partial {\varphi _{\rm{m}}}\rho }}{{\partial t}} = 0. $

(22)

采用文献[15]中的方法, 将式(20)代入式(22)中得到考虑时间滞后的SRV渗流数学模型:

$ \frac{1}{\tau }\int\limits_0^t {{{\rm{e}}^{ - \frac{{t - u}}{\tau }}}\nabla \frac{{{k_{\rm{m}}}}}{\mu }\left( {x,y} \right)\nabla p{\rm{d}}u} = {C_t}{\varphi _{\rm{m}}}\left( {x,y} \right)\frac{{\partial p}}{{\partial t}}. $

(23)

式中, B为流体体积分数; C_t为SRV综合压缩系数, Pa^-1; p为流动压力, Pa; p₀为油藏初始压力, Pa; p_w为井底流动压力, Pa; Q为产量, m³/s; t为时间, s; τ为滞后时间, s。

将式(23)无因次化, 定义如下无因次参数：

$ {y_{\rm{D}}} = \frac{y}{{{l_{\rm{f}}}}}, $

(24)

$ {x_{\rm{D}}} = \frac{x}{{{l_{\rm{f}}}}}, $

(25)

$ {w_{\rm{D}}} = \frac{w}{{{l_{\rm{f}}}}}, $

(26)

$ {p_{\rm{D}}} = \frac{{2{\rm{ \mathsf{ π} }}nh{k_0}\left( {{p_0} - p} \right)}}{{QB\mu }}, $

(27)

$ {t_{\rm{D}}} = \frac{{{k_0}t}}{{{\varphi _0}\mu {C_t}l_{\rm{f}}^2}}, $

(28)

$ {k_{{\rm{mD}}}} = \frac{{{k_m}}}{{{k_0}}}, $

(29)

$ {\varphi _{{\rm{mD}}}} = \frac{{{\varphi _m}}}{{{\varphi _0}}}, $

(30)

$ {\tau _{\rm{D}}} = \frac{{{k_0}\tau }}{{{\varphi _0}\mu {C_t}l_{\rm{f}}^2}}. $

(31)

式中, p_D表示无因次压力; t_D表示无因次时间; τ_D表示无因次滞后时间; SRV无因次渗透率k_mD、无因次孔隙度φ_mD沿x_D, y_D轴递减; (0, 0)为水力裂缝与井筒交点, 此处SRV渗透率为k₀, 孔隙度为φ₀。

无因次SRV渗流数学模型为

$ \frac{1}{{{\tau _{\rm{D}}}}}\int\limits_0^{{t_{\rm{D}}}} {{{\rm{e}}^{ - \frac{{{t_{\rm{D}}} - u}}{{{\tau _{\rm{D}}}}}}}\nabla {k_{{\rm{mD}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}}} \right)\nabla {p_{\rm{D}}}{\rm{d}}u} = {\varphi _{{\rm{mD}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}}} \right)\frac{{\partial {p_{\rm{D}}}}}{{\partial {t_{\rm{D}}}}}. $

(32)

2.3 水力裂缝流动模型

在水力裂缝中, 取长dy, 宽w, 高h的单元体, 单元体流动连续性方程为

$ \nabla \rho \mathit{\boldsymbol{v}} + \frac{{\partial {\varphi _{\rm{f}}}\rho }}{{\partial t}} + \rho {q^ * } = 0. $

(33)

式中, h为地层厚度, m; w为水力裂缝宽度, m; q^*为SRV到水力裂缝的窜流量。水力裂缝与SRV之间的窜流存在不平衡效应, 窜流量如下:

$ {q^ * } = - \frac{1}{\tau }\int\limits_0^t {{{\rm{e}}^{ - \frac{{t - u}}{\tau }}}\left( {\frac{{2{k_{\rm{m}}}}}{{\mu w}}\frac{{\partial p}}{{\partial x}}\left| {_{y = 0}} \right.} \right){\rm{d}}u} + \frac{{2{C_{{\rm{gd}}}}}}{{hw{l_{\rm{f}}}}}\frac{{\partial p}}{{\partial t}}. $

(34)

式中, C_gd为窜流不平衡效应因子, m³/Pa, 表示流动截面变化处的流动压力减小时, 截面附近的储容流体的流出量。可以看出, 窜流量由两部分组成:SRV中流体流出量(式(34)的第一项)和SRV、水力裂缝交界面附近储容流体流出量(式(34)的第二项)。

人工裂缝平面迂曲, 且不断分叉, 如图 2(b)所示，因此人工裂缝中流动也存在时间滞后效应。类似于式(12), 通过单元体的流量可表示为

$ {\mathit{\boldsymbol{v}}_y}hw = Q\left( t \right) = - H\left( 0 \right)\left( {\int\limits_0^t {\frac{{{{\rm{e}}^{ - \frac{{t - u}}{\tau }}}}}{\tau }\left| {\nabla p} \right|{\rm{d}}u} } \right). $

(35)

根据立方定律, 其中

$ H\left( 0 \right) = \frac{{h{w^3}}}{{24\mu }}, $

(36)

水力裂缝渗透率k_f可表示为

$ {k_{\rm{f}}} = \frac{{2\mu H\left( 0 \right)}}{{hw}} = \frac{{{w^2}}}{{12}}. $

(37)

采用文献[15]中的方法, 将式(34)、(35)代入式(33)中得到水力裂缝渗流数学模型:

$ \begin{array}{l} \frac{1}{\tau }\int\limits_0^t {{{\rm{e}}^{ - \frac{{t - u}}{\tau }}}\left( {\nabla \frac{{{k_{\rm{f}}}}}{\mu }\nabla p + \frac{{2{k_{\rm{m}}}}}{{\mu w}}\frac{{\partial p}}{{\partial x}}\left| {_{y = 0}} \right.} \right){\rm{d}}u} = \\ {C_{{\rm{lf}}}}{\varphi _{\rm{f}}}\frac{{\partial p}}{{\partial t}} + \frac{{2{C_{{\rm{gd}}}}}}{{hw{l_{\rm{f}}}}}\frac{{\partial p}}{{\partial t}}, \end{array} $

(38)

井筒存在储容效应:

$ - \frac{{hw{k_{\rm{f}}}}}{{2\mu }}\frac{{\partial p}}{{\partial y}}\left| {_{x = 0}} \right. = \frac{{BQ}}{{2n}} + 2{C_{\rm{w}}}\frac{{\partial p}}{{\partial t}}. $

(39)

式中, φ_f为水力裂缝孔隙度; n为水力裂缝级数; l_f为水力裂缝半长, m; C_lf为水力裂缝综合压缩系数, Pa^-1; C_w为井筒储容系数, m³/Pa。

将式(38)、(39)无因次化, 定义如下无因次参数:

$ {k_{{\rm{fD}}}} = \frac{{{w^2}}}{{12{k_0}}}, $

(40)

$ {\varphi _{{\rm{fD}}}} = \frac{{{\varphi _{\rm{f}}}{C_{{\rm{lf}}}}}}{{{\varphi _0}{C_t}}}, $

(41)

$ {C_{{\rm{wD}}}} = \frac{{2{C_{\rm{w}}}}}{{{\rm{ \mathsf{ π} }}{\varphi _0}{C_t}hl_{\rm{f}}^2}}, $

(42)

$ {C_{{\rm{Dgd}}}} = \frac{{2{C_{{\rm{gd}}}}}}{{hw{\varphi _0}{C_t}l_{\rm{f}}^3}}. $

(43)

式中，C_wD为无因次井筒储容系数; C_Dgd为无因次窜流不平衡效应因子; φ_fD为无因次水力裂缝孔隙度; k_fD为无因次水力裂缝渗透率; w_D为无因次水力裂缝宽度。

得到无因次水力裂缝渗流数学模型:

$ \begin{array}{l} \frac{1}{{{\tau _{\rm{D}}}}}\int\limits_0^{{t_{\rm{D}}}} {{{\rm{e}}^{ - \frac{{{t_{\rm{D}}} - u}}{{{\tau _{\rm{D}}}}}}}\left( {\nabla {k_{{\rm{fD}}}}\nabla {p_{\rm{D}}} + \frac{{2{k_{\rm{D}}}}}{{{w_{\rm{D}}}}}\frac{{\partial {p_{\rm{D}}}}}{{\partial {x_{\rm{D}}}}}\left| {_{{y_{\rm{D}}} = 0}} \right.} \right){\rm{d}}u} = \\ {\varphi _{{\rm{fD}}}}\frac{{\partial {p_{\rm{D}}}}}{{\partial {t_{\rm{D}}}}} + {C_{{\rm{Dgd}}}}\frac{{\partial {p_{\rm{D}}}}}{{\partial {t_{\rm{D}}}}}, \end{array} $

(44)

无因次井筒流动模型:

$ \frac{{{w_{\rm{D}}}{k_{{\rm{fD}}}}}}{2}\frac{{\partial {p_{\rm{D}}}}}{{\partial {y_{\rm{D}}}}}\left| {_{{x_{\rm{D}}} = 0}} \right. = {\rm{ \mathsf{ π} }}\left( {1 - {C_{{\rm{wD}}}}\frac{{\partial {p_{\rm{D}}}}}{{\partial {t_{\rm{D}}}}}} \right). $

(45)

3 模型求解

对式(32)、(44)、(45)进行Laplace变换得

$ \frac{{\nabla {k_{{\rm{mD}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}}} \right)\nabla \overline {{p_{\rm{D}}}} }}{{1 + s{\tau _{\rm{D}}}}} = s{\varphi _{{\rm{mD}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}}} \right)\overline {{p_{\rm{D}}}} , $

(46)

$ \frac{{\nabla {k_{{\rm{fD}}}}\nabla \overline {{p_{\rm{D}}}} + \frac{{2{k_{{\rm{mD}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}}} \right)\overline {{p_{\rm{D}}}} }}{{{w_{\rm{D}}}}}\left| {_{{y_{\rm{D}}} = 0}} \right.}}{{1 + s{\tau _{\rm{D}}}}} = s\left( {{\varphi _{{\rm{fD}}}} + {C_{{\rm{Dgd}}}}} \right)\frac{{\partial \overline {{p_{\rm{D}}}} }}{{\partial {t_{\rm{D}}}}}, $

(47)

$ \frac{{{w_{\rm{D}}}{k_{\rm{D}}}}}{2}\frac{{\partial {p_{\rm{D}}}}}{{\partial {y_{\rm{D}}}}}\left| {_{{x_{\rm{D}}} = 0}} \right. = {\rm{ \mathsf{ π} }}\left( {1 - {C_{{\rm{wD}}}}s\overline {{p_{\rm{D}}}} } \right). $

(48)

式中, $\overline {{p_{\rm{D}}}} $为p_D在Laplace空间中的像, 为s的独立函数; s为Laplace变换参数。

将无因次井底流动压力分解为若干部分, 不同压力部分由不同介质产生, 其变化的快慢由相应的时间尺度反映, 这就是对无因次井底流动压力的时间尺度分析。为确保上述分解的顺利进行, 对式(46)~(48)进行有限元离散得

$ \left( {\frac{K}{{1 + s{\tau _{\rm{D}}}}} + Cs} \right)\overline {{\mathit{\boldsymbol{p}}_{\rm{D}}}} = \frac{{\rm{ \mathsf{ π} }}}{{2s\left( {1 + s{\tau _{\rm{D}}}} \right)}}\mathit{\boldsymbol{B}}, $

(49)

其中

$ \begin{array}{l} K = \frac{{{w_{\rm{D}}}}}{2}\int\limits_\mathit{\Gamma } {\left( {\nabla \mathit{\boldsymbol{N}}} \right){k_{{\rm{fD}}}}\left( {\nabla {\mathit{\boldsymbol{N}}^{\rm{T}}}} \right){\rm{d}}\mathit{\Gamma }} + \\ \int\limits_A {\left( {\nabla \mathit{\boldsymbol{N}}} \right){k_{{\rm{mD}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}}} \right)\left( {\nabla {\mathit{\boldsymbol{N}}^{\rm{T}}}} \right){\rm{d}}A} , \end{array} $

(50)

$ \begin{array}{l} C = \frac{{{w_{\rm{D}}}}}{2}\int\limits_\mathit{\Gamma } {\mathit{\boldsymbol{N}}\left( {{\varphi _{{\rm{fD}}}} + {C_{{\rm{Dgd}}}}} \right){\mathit{\boldsymbol{N}}^{\rm{T}}}{\rm{d}}\mathit{\Gamma }} + \\ \int\limits_A {\mathit{\boldsymbol{N}}{\varphi _{{\rm{mD}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}}} \right){\mathit{\boldsymbol{N}}^{\rm{T}}}{\rm{d}}A} + \frac{{{\rm{ \mathsf{ π} }}{C_{{\rm{wD}}}}}}{2}. \end{array} $

(51)

式中, A表示SRV; Γ为水力裂缝区域; $\overline{\mathit{\boldsymbol{p}}_{\rm{D}}} $为Laplace空间中各点无因次流动压力组成的向量; N为形函数组成的向量; 列向量B除第一个元素为1外, 其余元素都为0。

求得矩阵K和C的广义特征值和广义特征向量, 使得

$ PK{P^{\rm{T}}} = \Lambda , $

(52)

$ PC{P^{\rm{T}}} = I. $

(53)

式中, Λ为对角矩阵; P为由矩阵K、C广义特征向量组成的矩阵的逆。

Laplace空间中各点无因次流动压力为

$ \overline {{\mathit{\boldsymbol{p}}_{\rm{D}}}} = \frac{{\rm{ \mathsf{ π} }}}{2}{P^{\rm{T}}}{\left( {{\tau _{\rm{D}}}{s^3}I + {s^2}I + \Lambda s} \right)^{ - 1}}P\mathit{\boldsymbol{B}}, $

(54)

取式(54)中与(0, 0)点对应的行, 井筒无因次流动压力可表示为

$ \overline {{p_{{\rm{wD}}}}} = \sum\limits_{i = 1}^n {{\alpha _i}\overline {{p_i}} } , $

(55)

$ \overline {{p_i}} = \frac{{{\lambda _i}}}{{{\tau _{\rm{D}}}{s^3} + {s^2} + {\lambda _i}s}}. $

(56)

式中, n为离散得到的单元的数目; p_i描述了与时间尺度λ_i对应的无因次流动压力部分随时间变化的规律; $\overline {{p_i}} $为p_i在Laplace空间中的像。每一部分无因次流动压力先从0增加到最终值, 而后不再随时间变化。时间尺度λ_i为矩阵K和C的广义特征值, 即矩阵Λ的对角元素, 其倒数为时间常数^[16-17], 时间常数为某种介质产生的流动压力升到其最终值的68.3%所需的时间; α_i为时间尺度λ_i对井底流动压力的贡献度, 表示与时间尺度λ_i相对应的压力部分变化的幅度, 数值上等于该无因次压力部分增加后的最终值, 在与时间尺度λ_i相对应的储层介质中, α_i随着流动阻力的增加而增加, 并随着流体采出量的增加而增加。以λ_i为横轴, 以α_i为纵轴, 进而可以作出时间尺度分布图。

通过对式(56)解析反演, 求得p_i, λ_i≠0时, 有

$ {p_i} = \frac{{\frac{{1 - {{\rm{e}}^{ - {x_1}{t_{\rm{D}}}}}}}{{{x_1}}} - \frac{{1 - {{\rm{e}}^{ - {x_2}{t_{\rm{D}}}}}}}{{{x_2}}}}}{\Delta }, $

(57)

其中

$ \Delta = \sqrt {1 - 4\tau {\lambda _i}} , $

(58)

$ {x_1} = \frac{{1 + \delta }}{{2\tau }}, $

(59)

$ {x_2} = \frac{{1 + \delta }}{{2\tau }}. $

(60)

λ_i=0时，有

$ {p_i} = {t_{\rm{D}}} - \tau + \tau {{\rm{e}}^{ - \frac{{{t_{\rm{D}}}}}{\tau }}}. $

(61)

最后得到真实空间中井底无因次流动压力:

$ {p_{{\rm{wD}}}} = \sum\limits_{i = 1}^n {{\alpha _i}{p_i}} . $

(62)

4 计算实例 4.1 模型验证

对于图 1(a)所示的物理模型, 假定其x方向长度为100 m, y方向长度为100 m, 厚度为10 m。SRV渗透率k_m沿x_D方向从8.333×10^-14 m²线性衰减到8.333×10^-15 m², 如图 3(a)所示。其中x_D、y_D方向分别与x、y方向一致, 孔隙度为0.1, 综合压缩系数为4×10^-9 Pa^-1。井筒储容系数为3.14×10^-8 m³/Pa, 窜流不平衡效应因子为1×10^-11 m³/Pa。水力裂缝共有10级, 宽度为0. 001 m, 孔隙度为0.5, 综合压缩系数为4×10^-9Pa^-1。流体黏度为0.001 Pa·s, 体积系数为1.2, 井筒以产量5.787×10^-3m³/s生产。将模型划分为70×70网格。通过对Laplace空间中的解析解进行数值反演验证模型, 井底无因次压力随无因次时间变化的曲线如图 3(b)所示。可以看出, 2种方法计算得到的井底无因次流动压力曲线重合, 说明本文研究的时间尺度分析方法的正确性, 并且解决了时间上的离散带来的不收敛问题。此外, 在空间上的数值求解还避免了因为孔隙度、渗透率存在强烈的非均质性、各向异性, 模型在Laplace空间中的解析解难以求得, 甚至不能求得。得到的流动压力解是时间的连续函数, 可分解成若干部分, 每部分由不同介质产生。压裂方案的选择和生产参数的调整取决于每部分的变化快慢, 否则容易错过调整生产参数的时机, 而本方法可以揭示每部分的变化快慢, 对致密油藏的高效开发具有十分重要的指导作用。

4.2 时间滞后对井底无因次压力的影响

图 4(a)为不同滞后时间下井底无因次压力随无因次时间的变化。可以看出, 流动的滞后时间主要影响早期的无因次井底压力, 其中τ=960 s时井底早期流动压力最低, τ=96 s时井底早期流动压力次之, 不考虑时间滞后的井底早期流动压力最高。因为随着裂缝、孔隙迂曲度增大, τ也相应增大, 流动压力传播所需的时间增大, 生产早期流动压力传播范围相应地减小。然而, 当流动压力传遍整个SRV时, 无因次井底压力较不考虑时间滞后的井底流动压力增加得更快, 最后, 逐渐接近无时间滞后的井底流动压力, τ越小, 接近越快。因此水力压裂形成的迂曲、分叉裂缝可减慢流动压力在水力裂缝中的传播, 使SRV中的流体更早地采出。

图 4(b)为时间尺度分布图。可以看出:SRV渗透率连续变化导致了时间尺度的连续分布, 最大时间尺度的贡献度最小, 对应近水力裂缝地带; 最小时间尺度的贡献度最大, 对应远水力裂缝地带, 说明远离人工裂缝的区域如同封闭边界, 流动阻力较大。另外, 时间尺度分布呈两个条带, 它们之间的距离说明人工裂缝和SRV之间的不连续性。上面的条带对应着SRV渗流, 时间尺度跨度较大, 条带左侧较右侧稀疏, 说明人工裂缝附近的区域对井底流动影响较大; 而下面的条带对应着人工裂缝, 时间尺度跨度较小。流动的时间滞后不影响时间尺度分布, 只影响各时间尺度对应的压力部分随时间变化的规律。

4.3 非均质性对井底无因次压力的影响

考虑SRV孔隙度、渗透率沿x_D、y_D两个方向均发生变化。假设SRV孔隙度从0.1沿x_D方向线性递减至0.001, 沿y_D方向线性递减至0.001。渗透率存在各向异性, x_D方向的无因次渗透率沿x_D方向从8.333×10^-14 m²线性递减至1.667×10^-14 m², 沿y_D方向线性递减至1.667×10^-17 m², 而y_D方向的无因次渗透率从8.333×10^-15 m²沿x_D方向线性递减至8.333×10^-16 m², 沿y_D方向线性递减至8.333×10^-19 m²。井底无因次压力随无因次时间变化见图 5(a)。可以看出, 考虑y_D方向非均质性和各向异性, 早期井底流动压力略微减小, 说明近井地带x_D方向渗流阻力较小; 晚期井底流动压力略微增加, 说明远井地带渗流阻力较大, 孔隙流体压力较小, 因此远井地带向井底的供液较少。定产量生产时, 流动压力难以传至远井低孔低渗带, 因此渗透率、孔隙度分布对井底流动的影响极其微弱。

渗流不存在时间滞后时, y_D方向非均质性和各向异性不影响早期井底流动, 仅在流动压力传遍水力裂缝后, 井底压力略增加。因水力裂缝不再迂曲、分叉, 流动压力早期主要沿水力裂缝传播, 而不沿x_D方向在SRV内传播, 因此早期井底压力不再因为SRV渗透率各向异性而减小。然而晚期井底流动压力较考虑时间滞后的井底压力略小, 说明介质迂曲度低时, 早期压力梯度不再影响后续流动，接近达西定律。时间滞后增加流动压力在水力裂缝中传播的时间, 进而增加y_D方向非均质性和各向异性对井底流动的影响。

图 5(b)为时间尺度分布图。考虑y_D方向非均质性和各向异性, 最大、最小时间尺度都增加, 说明储层孔隙度越低, 流动压力随时间变化越剧烈, 供液量随时间递减也就越快。水力裂缝对应的条带向右移动, 时间尺度范围变窄, 进一步说明了流动压力只能在近井高孔高渗地带传播, 难以传播到远井地带。另外, 时间尺度分布图中上下两条带之间的距离减小, 说明y_D方向渗透率和远井渗透率减小, 近井地带采出的流体量相应地增加, 因此该地带对井底流动影响较大。

图 6为不同的时间尺度下, 时间尺度贡献度随空间的分布。某介质产生的压力部分的最终值反映了从不同区域采出的流体量。不考虑y_D方向非均质性和各向异性时, 从图 6(a)看出, 较小的时间尺度在远人工裂缝地带贡献度大, 说明该地带流动阻力大; 而在远井地带贡献度小, 该地带供液量减少, 说明相应的流动压力部分难以在低渗透介质中传播。然而从图 6(b)看出, 较大的时间尺度在远井、远人工裂缝地带贡献较小, 说明相应的流动压力部分尚未传到该地带, 该地带供液相应减少; 在远井地带, 该时间尺度贡献度在近人工裂缝地带较大, 说明储层渗透率越高, 流动压力随时间递减越快; 而在近井地带, 流动压力未经人工裂缝就传到了近井、远人工裂缝地带, 近井地带供液较多, 该时间尺度贡献度也随着流动阻力增加而增加, 同图 6(b)下部分趋势一致。若考虑y_D方向非均质性和各向异性, 从图 6(c)看出, 较小时间尺度的贡献度沿x_D方向等势分布, 分布条带平行于x_D方向, 说明y_D方向渗透率较小, 再加上远井地带流动阻力大, 因此流动压力主要沿x_D方向传播。同时, 近井、近人工裂缝地带采出流体量略增加, 进一步说明了流动压力难以在低渗透介质中传播。从图 6(d)看出, 较大时间尺度的贡献度沿y_D方向先增加后减小, 最后增加, 中间的减小阶段说明相应的压力部分不再随时间增加时, 流动压力尚未传遍人工裂缝, 远井地带供液相应减少; 时间尺度贡献度沿x_D方向递减, 而不再等势分布, 说明流动压力在较短时间内无法传到远人工裂缝地带。

对比井底无因次压力随无因次时间变化和时间尺度贡献度随时间尺度、空间的分布, 可以看出时间尺度分析可以准确表明储层非均质性、各向异性对井底流动的影响, 因为时间尺度贡献度随时间尺度、空间的分布在不同渗透率、孔隙度储层中差异明显。另外, 储层强烈的非均质性是导致多个时间尺度在整个SRV内共存的原因，而且时间尺度和空间区域互相对应, 因为在流动压力传播到与某时间尺度对应的区域, 该时间尺度贡献度在该区域较大。

4.4 窜流的不平衡效应对井底无因次压力的影响

图 7(a)为不同窜流不平衡效应因子下, 井底流动压力随无因次时间变化的曲线。可以看出, C_gd主要影响晚期的无因次井底流动压力, 因为早期仅水力裂缝中的流体大量采出, 然而SRV中的流体未进入人工裂缝。晚期SRV中的流体开始进入水力裂缝, 流动截面于SRV、水力裂缝接触面发生了变化, 导致接触面附近的储容流体流向井筒, 因此井底流动压力轻微降低。窜流的不平衡效应对井底流动的影响晚于且小于流动的时间滞后效应产生的影响, 即流动压力通过迂曲的水力裂缝才能传到SRV, 接触面附近的储容流体才开始流向井筒。

渗流不存在时间滞后时, 窜流的不平衡效应使井底流动压力轻微降低, 而且降低的时间早于考虑时间滞后的。因水力裂缝不再迂曲, 流动压力很快传到SRV。然而流动存在时间滞后时, 迂曲的水力裂缝与SRV接触面积较大, 且流动压力在水力裂缝中传播缓慢, 因此在生产早期, 从接触面附近流出的流体较流动不存在时间滞后时少。

图 7(b)为时间尺度分布图。介质间耦合存在不平衡效应时, 最大、最小时间尺度都减小, 时间尺度分布图上方还比不存在不平衡效应的多一条连续曲线, 曲线右边比左边陡, 说明接触面附近的流体比SRV流体先采出, 因此SRV流动压力变化晚于人工裂缝流动压力。上述曲线高于不存在不平衡效应的时间尺度分布带, 说明压力传到远人工裂缝地带时, 体积改造区域流出量接近人工裂缝流入量, 因此SRV流出的流体较多。但时间尺度分布带下方比不存在不平衡效应的密集, 说明流动压力刚传到SRV时, 流入人工裂缝的流体不仅包括SRV流出的流体, 还包括从交界面附近流出的储容流体, 加上流动压力在近井地带传播较快, 近井和近人工裂缝地带流出的流体相应减少。另外, 介质间耦合存在不平衡效应时, 最大时间尺度贡献度高于不存在不平衡效应的时间尺度贡献度。对比井底无因次压力随无因次时间变化和时间尺度贡献度随时间尺度的分布, 可以得出时间尺度分析方法可以系统地表明不同介质间的渗流速度的差异以及耦合的不平衡效应。

5 结论

(1) 时间尺度与空间区域相对应, 并且储层强烈的非均质性是导致多个时间尺度共存的原因; 近井压裂高渗区域对井底流动的影响较大, 而远井低渗区域对井底流动的影响较小。

(2) 流动的时间滞后效应影响各介质产生的压力部分随时间变化的规律, 使早期井底流压降低, 迂曲、分叉的水力裂缝可减慢压力的传播, 进而促使SRV中的流体更早地采出; 而窜流的不平衡效应对晚期渗流影响较大, 同时可减弱远井低渗地带对井底流动的影响, 并且能够更好地描述SRV、水力裂缝之间的渗流速度差异。

(3) 水力裂缝的不连续性使得时间尺度分布图呈现两个条带; 渗透率、孔隙度非均质性加剧不同区域采出流体量的差异, 但时间尺度分布图中上、下两个条带之间的距离相应地减小; 耦合存在不平衡效应的时间尺度较不考虑不平衡效应的小, SRV流出的流体较少。对井底流动压力的时间尺度分析能够清晰描述不同渗透率、孔隙度储层对井底流动压力的影响规律, 进而对生产制度的制定具有重要的指导作用。