纵波在地下含流体岩石中传播时会引起岩石孔隙内流体的流动, 按照尺度的不同, 可以将流体的流动大致分为宏观尺度、微观尺度和中观尺度的流动。同时, 不同尺度的流体流动都会对纵波的传播产生影响, 使纵波出现衰减和速度频散现象^[1]。

BIOT针对孔隙内流体的宏观流动对弹性波的影响给出了比较完整的论述, 形成了一套研究双相介质弹性波传播问题的Biot理论^[2-4]。Biot理论认为, 当弹性波在各向同性均匀多孔介质内传播时, 孔隙内的流体会产生宏观尺度上的流动。由于流体的惯性作用, 使得流体的运动速度小于固体的运动速度, 流体流动的滞后使得流体与固体之间产生了内摩擦力, 从而导致弹性波的衰减和速度频散。虽然Biot理论成功预测了慢纵波的存在, 但由于理论假设的局限性, Biot理论得到的弹性波衰减和速度频散程度都远小于实验结果^[5]。

MAVKO等^[6]认为含流体多孔岩石中还存在一种微观尺度的喷射流动, 且这种流动会造成弹性波更强的衰减和速度频散。造成流体产生喷射流动的主要原因是当弹性波在介质内传播时, 由于介质内孔隙可压缩性的差异, 使流体从裂缝中进入邻近孔隙; 或由于介质部分饱和, 使流体进入含高压缩性气体的孔隙。

对于部分饱和岩石, 由于存在局部非均匀性, 介质内部还会存在一种中观尺度的流体流动。很多学者都对部分饱和岩石中弹性波的传播规律进行了研究。White理论详细论述了弹性波在具有局部非均匀性的部分饱和模型中的传播规律^[7-8]; DUTTA等^[9-10]通过求解Biot方程, 分析了White模型内的衰减及速度频散现象; 杜伟等^[11]构建了一个裂隙孔隙介质, 给出了垂直于周期性层状孔隙介质层面传播的纵波的频散方程, 再结合线性滑动裂隙模型推导出垂直于裂隙传播的纵波的等效介质模型; 为了研究纵波的耗散过程, SUN等^[12]提出了一种饱含两种不相容流体的3层部分饱和模型; CHEN等^[13]通过分析地震波衰减和反射行为, 认为波致流体流动受多相流体的影响; 未晛^[14]基于随机斑块饱和孔隙介质模型(CRM), 建立了改进的随机斑块饱和孔隙介质模型(MCRM), 研究了中观尺度斑块饱和对地震波衰减和速度频散特征的影响。

Biot理论、喷射流机制和White理论都是从单一尺度流体流动的角度出发, 对双相介质内的弹性波传播规律进行研究。但当弹性波在地下实际岩石中传播时, 孔隙内流体的流动形式十分复杂, 用单一尺度的流动模型进行描述不合理。

DVORKIN等^[15]认为, 当纵波在含流体多孔介质内传播时, 流体的宏观流动和微观喷射流动同时存在。他们同时考虑了这两种尺度的流动, 建立了一种统一的一维Biot-喷射流(Squirt)模型(BISQ模型), 并推导了纵波相速度和品质因子的表达式, 研究了纵波的衰减和速度频散现象。BISQ模型是一种典型的孔隙弹性理论, 很多学者都以BISQ模型为基础, 研究孔隙介质内流体的宏观和微观流动。DIALLO等^[16]对BISQ模型进行了改进, 重新推导了介质内的平均流体压力项, 使其不含喷射长度; 杨顶辉等^[17]首先对各向同性介质的BISQ模型进行了详细的回顾, 并提出了各向异性介质的BISQ模型; 杜艺可等^[18]以Biot理论为基础, 利用平面波分析推导了纵波相速度和衰减因子的解析表达式, 并与BISQ理论进行了深入对比; YANG等^[19]建立了一个含多相流的黏弹性BISQ模型, 研究了低孔低渗条件下的波衰减和速度频散现象; LIU等^[20]对BISQ模型中喷射长度为负值或具有非零虚部的情况进行了研究; 杨盈盈等^[21]以BISQ模型为基础, 改进了Pride震电耦合理论, 研究了震电耦合波的传播特性; 凌云等^[22]将Zener线性体和BISQ模型的作用机理相结合, 建立了Zener-BISQ模型, 研究了纵波在新模型中的衰减和速度频散现象。

与单一尺度的Biot理论相比, BISQ模型虽然与实际情况更为吻合, 但其仍无法解释弹性波在地震频段内的衰减和频散现象。造成弹性波在地震频段内衰减和频散的主要原因是流体的中观流动, 因此很多学者将流体的中观流动与宏观流动相结合, 研究弹性波在含流体孔隙介质内的传播规律。BA等^[23]同时考虑了流体的宏观流动和中观流动, 建立了一种双孔模型, 推导了该模型的波动方程, 研究了双孔模型中的局部流体流动过程。BA等^[24]在双孔模型的基础上, 考虑了球状内含物内部的流体动能, 重新推导了模型的波动方程; BA等^[25]提出了一种新的双孔模型, 研究了含黏土孔隙的致密粉砂岩中的强速度频散现象; SONG等^[26]提出了有效动态剪切模量的解析解, 并通过分析S波衰减和速度频散, 研究了波致流体流动与剪切模量之间的关系; BA等^[27]建立了一种双重双孔模型并推导了模型的控制方程, 研究了构造和饱和度的非均匀性对弹性波衰减和速度频散的影响; ZHAO等^[28]将White层状部分饱和模型与双孔模型相结合, 建立了层状双孔模型, 研究了纵波的衰减和速度频散。为了预测双孔介质中的衰减和速度频散, ZHENG等^[29]推导了一种新的波动方程, 并利用解析法求解该方程。很多学者以双孔模型为基础, 研究了两种尺度流体流动的共同作用对弹性波反射系数的影响。DAI等^[30]计算了流体饱和单孔和双孔介质分界面处的反射和透射系数, 认为反射和透射系数与频率和入射角有关。为了分析波致流体流动对反射系数的影响, ZHAO等^[31]将PRIDE等^[32]建立的双孔模型作为研究背景, 计算了弹性波在模型分界面处的反射和透射系数; 贺鹏飞等^[33]认为P波在半空间自由边界处的反射系数与入射角和频率有关, 并利用双孔介质模型理论对这种设想进行了验证。

BISQ理论同时考虑了流体的宏观流动和微观喷射流动, 双孔模型将宏观流动和中观局部流动相结合, 但目前很少有文献提到将3种尺度流体流动同时考虑的模型。本文在ZHAO等^[28]建立的层状双孔模型的基础上, 考虑微观尺度喷射流的影响, 建立含横向喷射流的层状部分饱和模型, 研究纵波的衰减和速度频散特性。首先建立含横向喷射流的层状部分饱和模型; 然后通过求解模型中的平均流体压力推导波动方程, 计算纵波的相速度和品质因子; 最后, 利用数值模型分析了3种尺度共同作用下纵波的衰减和速度频散特性, 同时研究了岩石参数对衰减和速度频散的影响。

1 基础理论 1.1 含横向喷射流的层状部分饱和模型

以ZHAO等^[28]建立的层状双孔模型为基础, 引入孔隙内流体的横向喷射流, 建立含横向喷射流的层状部分饱和模型, 如图 1所示。

WHITE等^[8]建立了一个含气层和含水层交替排列的周期性层状部分饱和模型, 如图 1a所示, 图中h₁和h₂表示模型每一层的厚度。我们以图 1a中的黑框为纵向切面, 在周期层状部分饱和模型中选取一个两层的圆柱形单元体, 如图 1b所示, 其中上层介质内饱含水, 下层介质内饱含气, X₁和X₂分别表示单元体第一层的厚度和模型的总厚度。当纵波沿图中黑色箭头所示方向在模型介质内传播时, 孔隙内的流体会产生宏观尺度的流动, 如图中蓝色箭头所示。由于模型上、下两层介质内饱含了两种不同的流体, 所以在两层介质间还会产生中观尺度的局部流体流动, 如图中绿色箭头所示。在实际地下岩石中, 流体微观尺度的喷射流动与宏观尺度的流动同时发生, 所以在模型中再引入垂直于纵波传播方向的横向喷射流, 如图 1b中的红色箭头所示, 这样就在同一个模型中同时考虑了宏观、微观和中观尺度的流体流动。图 1b中的R表示横向喷射长度, 同时表示圆柱形单元体的半径。

1.2 波动方程

BIOT利用拉格朗日方程推导了各向同性多孔介质的波动方程^[4], 本文将此方法推广到含横向喷射流的层状部分饱和模型中, 可以得到一维拉格朗日方程为:

$ \begin{array}{l} \frac{{\partial \tau }}{{\partial x}} = \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial T}}{{\partial {{\dot u}_x}}}} \right) + \frac{{\partial D}}{{\partial {{\dot u}_x}}}\\ \frac{{\partial {\tau ^{\left( m \right)}}}}{{\partial x}} = - {\varphi _m}\frac{{\partial p_{\rm{f}}^{\left( m \right)}}}{{\partial x}} = \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial T}}{{\partial \dot U_x^{\left( m \right)}}}} \right) + \frac{{\partial D}}{{\partial \dot U_x^{\left( m \right)}}} \end{array} $

(1)

式中:u_x和U_x^(m)分别表示波传播方向上的固体和流体位移, 变量上面的点表示对时间求导数; T和D分别表示模型的动能和耗散函数(具体表达式见附录A); φ_m和p_f^(m)(m=1, 2)分别表示上、下两层介质的绝对孔隙度和流体压力; τ和τ^(m)分别表示模型介质固体和流体部分的应力, 具体表达式为:

$ \begin{array}{l} \tau = \left( {A + 2N} \right)e + \sum\limits_m {{Q_m}{\xi _m}} \\ {\tau ^{\left( m \right)}} = {Q_m}e + {R_m}{\xi _m} \end{array} $

(2)

其中, A, N, Q_m和R_m表示刚度系数; e和ξ_m分别表示介质固体和流体部分的体应变。

将(2)式和(A1)式代入(1)式, 可以得到含横向喷射流的层状部分饱和模型的一维波动方程:

$ \begin{array}{l} \left( {A + 2N - \frac{{Q_1^2}}{{{R_1}}} - \frac{{Q_2^2}}{{{R_2}}}} \right)\frac{{{\partial ^2}{u_x}}}{{\partial {x^2}}} - \frac{{{\varphi _1}{Q_1}}}{{{R_1}}}\frac{{\partial p_{\rm{f}}^{\left( 1 \right)}}}{{\partial x}} - \\ \;\;\;\;\frac{{{\varphi _2}{Q_2}}}{{{R_2}}}\frac{{\partial p_{\rm{f}}^{\left( 2 \right)}}}{{\partial x}} = {\rho _{00}}{{\ddot u}_x} + {\rho _{01}}\ddot U_x^{\left( 1 \right)} + {\rho _{02}}\ddot U_x^{\left( 2 \right)} + \\ \;\;\;\;{b_1}\left( {{{\dot u}_x} - \dot U_x^{\left( 1 \right)}} \right) + {b_2}\left( {{{\dot u}_x} - \dot U_x^{\left( 2 \right)}} \right) \end{array} $

(3a)

$ - {\varphi _1}\frac{{\partial p_{\rm{f}}^{\left( 1 \right)}}}{{\partial x}} = {\rho _{01}}{{\ddot u}_x} + {\rho _{11}}\ddot U_x^{\left( 1 \right)} - {b_1}\left( {{{\dot u}_x} - \dot U_x^{\left( 1 \right)}} \right) $

(3b)

$ - {\varphi _2}\frac{{\partial p_{\rm{f}}^{\left( 2 \right)}}}{{\partial x}} = {\rho _{02}}{{\ddot u}_x} + {\rho _{22}}\ddot U_x^{\left( 2 \right)} - {b_2}\left( {{{\dot u}_x} - \dot U_x^{\left( 2 \right)}} \right) $

(3c)

其中, ρ₀₀, ρ₀₁, ρ₀₂, ρ₁₁, ρ₂₂表示模型介质的密度系数; b₁和b₂分别表示上层和下层介质的耗散系数; 变量上面两点表示对时间求二阶导数。

由于模型上、下两层间存在流体的局部流动, 故还需用一个控制方程来描述这种流动。本文将ZHAO等^[28]给出的控制方程推广到饱含两种流体的圆柱形单元体模型中, 可以得到:

$ \begin{array}{l} {\varphi _1}{Q_2}e + {\varphi _1}{R_2}\left( {{\xi _2} + {\varphi _1}\zeta } \right) - {\varphi _2}{Q_1}e - {\varphi _2}{R_1}\left( {{\xi _1} - } \right.\\ \;\;\;\;\left. {{\varphi _2}\zeta } \right) = \varphi _2^2X_1^2\left[ {\rho _{\rm{f}}^{\left( 1 \right)}{\varphi _1} + 3\rho _{\rm{f}}^{\left( 2 \right)}{\varphi _{20}}{\rm{ \mathsf{ π} }}{R^2}\left( {{X_2} - } \right.} \right.\\ \;\;\;\;\left. {\left. {\frac{{{\varphi _1}}}{{{\rm{ \mathsf{ π} }}{R^2}{\varphi _{10}}}}} \right)} \right]\ddot \zeta + \varphi _2^2X_1^2\left[ {\frac{1}{3}\frac{{{\eta _1}}}{\kappa }{\varphi _1}{\varphi _{10}} + \frac{{{\eta _2}}}{\kappa }\varphi _{20}^2{\rm{ \mathsf{ π} }}{R^2} \cdot } \right.\\ \;\;\;\;\left. {\left( {{X_2} - \frac{{{\varphi _1}}}{{{\rm{ \mathsf{ π} }}{R^2}{\varphi _{10}}}}} \right)} \right]\dot \zeta \end{array} $

(4)

其中, ζ表示模型介质中流体的局部流动造成的流体应变增量; ρ_f⁽¹⁾和ρ_f⁽²⁾分别表示上层和下层介质的流体密度; φ₁₀和φ₂₀分别表示上层和下层介质的局部孔隙度(φ₁₀=φ₂₀); η₁和η₂分别表示上层和下层介质内流体的黏滞系数; κ表示介质的渗透率。(3)式和(4)式即为含横向喷射流的层状部分饱和模型的波动方程, 但此时的(3)式中仍含有未知量p_f^(m), 下面我们将求解(3)式中的平均流体压力p_f^(m)。

1.3 3种尺度流体流动同时存在情况下的平均流体压力

DVORKIN等^[15]推导了含横向喷射流的平均流体压力方程, 本文将其推导过程推广到3种尺度流体流动同时存在的情况, 即加入中观尺度局部流体流动的影响。上层介质的流体质量守恒方程为:

$ \begin{array}{l} \frac{{\partial \left( {{\varphi _1}\rho _{\rm{f}}^{\left( 1 \right)}} \right)}}{{\partial t}} + \frac{{\partial \left[ {\rho _{\rm{f}}^{\left( 1 \right)}{\varphi _1}\left( {\dot U_x^{\left( 1 \right)} - {{\dot u}_x}} \right)} \right]}}{{\partial x}} + \\ \frac{{\partial \left( {\rho _{\rm{f}}^{\left( 1 \right)}{\varphi _1}{{\dot \nu }^{\left( 1 \right)}}} \right)}}{{\partial r}} + \frac{{\rho _{\rm{f}}^{\left( 1 \right)}{\varphi _1}{{\dot \nu }^{\left( 1 \right)}}}}{r} - \rho _{\rm{f}}^{\left( 1 \right)}{\varphi _1}{\varphi _2}\dot \zeta = 0 \end{array} $

(5)

式中:ν⁽¹⁾表示上层介质内流体的径向位移。根据DVORKIN等^[15]的推导, 孔隙度与流体密度的时间导数可以写为:

$ \frac{{\partial {\varphi _1}}}{{\partial t}} = {\alpha _1}\frac{{{\partial ^2}{u_x}}}{{\partial x\partial t}} + \frac{1}{{{{Q'}_1}}}\frac{{\partial \rho _{\rm{f}}^{\left( 1 \right)}}}{{\partial t}} $

(6)

$ \frac{{\partial \rho _{\rm{f}}^{\left( 1 \right)}}}{{\partial t}} = \frac{1}{{c_{01}^2}}\frac{{\partial \rho _{\rm{f}}^{\left( 1 \right)}}}{{\partial t}} $

(7)

其中,

$ {\alpha _1} = 1 - \frac{{{K_{{\rm{dl}}}}}}{{{K_{\rm{s}}}}} $

(8)

$ \frac{1}{{{{Q'}_1}}} = \frac{1}{{{K_{\rm{s}}}}}\left( {1 - {\varphi _1} - \frac{{{K_{{\rm{dl}}}}}}{{{K_{\rm{s}}}}}} \right) $

(9)

式中:K_s和K_d1分别表示模型介质的颗粒体积模量和上层介质的干燥体积模量; c₀₁表示上层介质内流体的声波速度。将(6)式、(7)式代入(5)式, 得到:

$ \begin{array}{l} \frac{{\partial \rho _{\rm{f}}^{\left( 1 \right)}}}{{\partial t}} = - {F_1}\left( {\frac{{{\partial ^2}U_x^{\left( 1 \right)}}}{{\partial x\partial t}} + \frac{{{\gamma _1}}}{{{\varphi _1}}}\frac{{{\partial ^2}{u_x}}}{{\partial x\partial t}} + \frac{{{\partial ^2}{\nu ^{\left( 1 \right)}}}}{{\partial r\partial t}} + } \right.\\ \;\;\;\;\;\;\;\;\;\left. {\frac{1}{r} - \frac{{\partial {\nu ^{\left( 1 \right)}}}}{{\partial t}} - {\varphi _2}\frac{{\partial \zeta }}{{\partial t}}} \right) \end{array} $

(10)

其中,

$ {F_1} = {\left( {\frac{1}{{c_{01}^2\rho _{\rm{f}}^{\left( 1 \right)}}} + \frac{1}{{{\varphi _1}{{Q'}_1}}}} \right)^{ - 1}} $

(11)

$ {\gamma _1} = {\alpha _1} - {\varphi _1} $

(12)

流体的横向喷射流动同样满足波动方程, 且固体部分在径向上没有移动, 故根据(3b)式可以得到:

$ {\rho _{11}}\frac{{{\partial ^2}{\nu ^{\left( 1 \right)}}}}{{\partial {t^2}}} + {b_1}\frac{{\partial {\nu ^{\left( 1 \right)}}}}{{\partial t}} = - {\varphi _1}\frac{{\partial \rho _{\rm{f}}^{\left( 1 \right)}}}{{\partial r}} $

(13)

将流体径向位移ν⁽¹⁾和流体压力p_f⁽¹⁾写成平面波解的形式:

$ {\nu ^{\left( 1 \right)}} = \nu _0^{\left( 1 \right)}{{\rm{e}}^{{\rm{i}}\left( {\omega t - kx} \right)}} $

(14)

$ p_{\rm{f}}^{\left( 1 \right)} = p_{f{\rm{0}}}^{\left( 1 \right)}{{\rm{e}}^{{\rm{i}}\left( {\omega t - kx} \right)}} $

(15)

并将此平面波解代入(13)式, 可以得到:

$ \frac{{\partial \rho _{f0}^{\left( 1 \right)}}}{{\partial r}} = \nu _0^{\left( 1 \right)}{\omega ^2}\left( {\frac{{{\rho _{11}}}}{{{\varphi _1}}} - {\rm{i}}\frac{{{b_1}}}{{\omega {\varphi _1}}}} \right) $

(16)

同样将固体位移u_x, 流体垂向位移U_x⁽¹⁾和流体应变增量ζ写成平面波解的形式:

$ {u_x} = {u_{x0}}{{\rm{e}}^{{\rm{i}}\left( {\omega t - kx} \right)}} $

(17)

$ U_x^{\left( 1 \right)} = U_{x0}^{\left( 1 \right)}{{\rm{e}}^{{\rm{i}}\left( {\omega t - kx} \right)}} $

(18)

$ \zeta = {\zeta _0}{{\rm{e}}^{{\rm{i}}\left( {\omega t - kx} \right)}} $

(19)

将(14)式到(19)式代入(10)式, 得到:

$ \begin{array}{l} \frac{{{\partial ^2}\rho _{f0}^{\left( 1 \right)}}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial \rho _{f0}^{\left( 1 \right)}}}{{\partial r}} + \frac{{{\omega ^2}\left( {\frac{{{\rho _{11}}}}{{{\varphi _1}}} - i\frac{{{b_1}}}{{\omega {\varphi _1}}}} \right)}}{{{F_1}}}\rho _{{\rm{f0}}}^{\left( 1 \right)}\\ = \left[ {{\rm{i}}k\left( {U_{x0}^{\left( 1 \right)} + \frac{{{\gamma _1}}}{{{\varphi _1}}}{u_{x0}}} \right) + {\varphi _2}{\zeta _0}} \right]{\omega ^2}\left( {\frac{{{\rho _{11}}}}{{{\varphi _1}}} - {\rm{i}}\frac{{{b_1}}}{{\omega {\varphi _1}}}} \right) \end{array} $

(20)

求解(20)式, 可以得到:

$ \begin{array}{l} \rho _{\rm{f}}^{\left( 1 \right)} = - {F_1}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _1}R} \right)}}{{{\lambda _1}R{J_0}\left( {{\lambda _1}R} \right)}}} \right]\left( {\frac{{\partial U_x^{\left( 1 \right)}}}{{\partial x}} + } \right.\\ \;\;\;\;\;\;\;\;\;\left. {\frac{{{\gamma _1}}}{{{\varphi _1}}}\frac{{\partial {u_x}}}{{\partial x}} - {\varphi _2}\zeta } \right) \end{array} $

(21)

式中:J₀和J₁分别表示零阶和一阶贝塞尔函数;

$ \lambda _1^2 = \frac{{{\omega ^2}\left( {\frac{{{\rho _{11}}}}{{{\varphi _1}}} - {\rm{i}}\frac{{{b_1}}}{{\omega {\varphi _1}}}} \right)}}{{{F_1}}} $

(22)

(21) 式即为同时存在3种尺度流体流动的上层介质内流体的平均压力。

类似的, 下层介质内流体的平均压力可以写为:

$ \begin{array}{l} \rho _{\rm{f}}^{\left( 2 \right)} = - {F_2}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _2}R} \right)}}{{{\lambda _2}R{J_0}\left( {{\lambda _2}R} \right)}}} \right]\left( {\frac{{\partial U_x^{\left( 2 \right)}}}{{\partial x}} + } \right.\\ \;\;\;\;\;\;\;\;\;\left. {\frac{{{\gamma _2}}}{{{\varphi _2}}}\frac{{\partial {u_x}}}{{\partial x}} + {\varphi _1}\zeta } \right) \end{array} $

(23)

式中:

$ {F_2} = {\left( {\frac{1}{{c_{02}^2\rho _{\rm{f}}^{\left( 2 \right)}}} + \frac{1}{{{\varphi _2}{{Q'}_2}}}} \right)^{ - 1}} $

(24)

$ \lambda _2^2 = \frac{{{\omega ^2}\left( {\frac{{{\rho _{22}}}}{{{\varphi _2}}} - {\rm{i}}\frac{{{b_2}}}{{\omega {\varphi _2}}}} \right)}}{{{F_2}}} $

(25)

$ {\gamma _2} = {\alpha _2} - {\varphi _2} $

(26)

$ \frac{1}{{{{Q'}_2}}} = \frac{1}{{{K_{\rm{s}}}}}\left( {1 - {\varphi _2} - \frac{{{K_{{\rm{d2}}}}}}{{{K_{\rm{s}}}}}} \right) $

(27)

$ {\alpha _2} = 1 - \frac{{{K_{{\rm{d2}}}}}}{{{K_{\rm{s}}}}} $

(28)

$ U_x^{\left( 2 \right)} = U_{x0}^{\left( 2 \right)}{{\rm{e}}^{{\rm{i}}\left( {\omega t - kx} \right)}} $

(29)

其中, c₀₂表示下层介质内流体的声波速度, K_d2表示下层介质的干燥体积模量。

将(21)式和(23)式代入(3)式, 就可以得到含横向喷射流的层状部分饱和模型的最终波动方程(具体过程见附录B)。

1.4 相速度和品质因子

为了研究3种尺度流体流动同时存在时纵波的衰减和速度频散现象, 利用平面波分析求解(3)式和(4)式, 计算纵波的相速度和品质因子。

将(17)式、(18)式、(19)式和(29)式代入(3)式和(4)式, 可以得到3类纵波(快纵波和第1、第2类慢纵波)的波数解(具体过程见附录C)。3类纵波的相速度和品质因子可以写为:

$ \begin{array}{*{20}{c}} {V = {{\left[ {{\mathop{\rm Re}\nolimits} \left( {\frac{k}{\omega }} \right)} \right]}^{ - 1}}}\\ {Q = \frac{{{\mathop{\rm Re}\nolimits} \left( {\frac{\omega }{k}} \right)}}{{2{\mathop{\rm Im}\nolimits} \left( {\frac{\omega }{k}} \right)}}} \end{array} $

(30)

其中, k表示3类纵波的复波数, ω表示圆频率。

2 数值模拟

根据(30)式可以求得3类纵波的相速度和品质因子, 分析纵波的衰减和速度频散特性。计算过程中需要的岩石参数如表 1所示, 其中μ_s和ρ_s分别表示固体颗粒的剪切模量和密度; K_f1和K_f2分别表示上、下两层介质内流体的体积模量。

图 2给出了利用5种孔隙弹性理论计算得到的快纵波相速度和衰减因子曲线。本文所建立的含横向喷射流的层状部分饱和模型以层状双孔模型为基础, 从图 5中可以看出, 在低频范围内, 这两种模型的频散和衰减曲线差别不大, 且都与采用White理论得到的频散和衰减曲线吻合较好。在高频范围内, 这两种模型的频散和衰减曲线差别很大, 引入了横向喷射流后, 频散程度大幅提高, 产生频散的频率范围增大, 衰减增大, 衰减峰向低频方向移动, 且高频范围内的频散和衰减曲线与采用BISQ理论得到的频散和衰减曲线吻合较好。层状双孔模型只考虑了孔隙流体的宏观和中观流动, 这种模型虽然能较好地描述纵波在低频范围内的频散和衰减现象, 但由于Biot理论得到的结果远小于实际情况, 故层状双孔模型无法较好地描述纵波在高频范围内的频散和衰减现象。含横向喷射流的层状部分饱和模型同时考虑了3种尺度的流体流动, 不仅保留了纵波在低频范围内的频散和衰减, 而且很好地描述了纵波在高频范围内的高频散和强衰减现象。

图 3给出了模型内含气饱和度不同时快纵波相速度和衰减因子曲线。从图 3中可以看出, 当模型内的含气饱和度逐渐增加时, 低频范围内的相速度先减小再逐渐增加, 频散程度逐渐减小, 低频衰减先增大再逐渐减小; 高频范围内的相速度逐渐减小, 频散程度也逐渐减小, 高频衰减变化不大。从图 3a中也可以看出, 当含气饱和度从0变化到10%时, 相速度明显减小, 说明相速度对模型中是否存在气体较为敏感。此外, 在低频范围内, 当含气饱和度在较低范围内变化时相速度曲线变化较为明显; 在高频范围内, 当含气饱和度在较高范围内变化时相速度曲线变化较为明显, 说明相速度曲线对模型的含气饱和度有较好的指示作用。

图 4是模型内孔隙度不同时快纵波相速度和衰减因子曲线图。从图 4中可以看出, 随着模型孔隙度的减小, 频散曲线向速度增加的方向移动, 衰减峰值略有下降。此外, 随着孔隙度的减小, 低频范围内的频散曲线和衰减峰向低频方向移动, 这说明在低孔低渗岩石中, 纵波更容易在地震频段内产生衰减和速度频散现象。

图 5是模型内喷射长度不同时快纵波相速度和衰减因子曲线图。从图 5中可以看出, 随着喷射长度的增加, 高频范围内的频散曲线和衰减峰向低频方向移动, 同时由于喷射长度与中观尺度流体流动过程无关, 所以低频范围内的频散曲线和衰减峰无明显变化。

模型内不仅存在快纵波, 还存在两类慢纵波, 两类慢纵波相速度和衰减因子曲线如图 6所示。从图 6中可以看出, 在低频范围内, 两类慢纵波相速度基本为0, 在高频范围内相速度随着频率的增加而增大, 但第一类慢纵波相速度增加更明显。第一类慢纵波的衰减因子在低频范围内随着频率的增加而增大, 在高频范围内随着频率的增加而减小; 第二类慢纵波在低频范围内衰减因子几乎为0, 在高频范围内随着频率的增加而减小。

3 结论

本文在层状双孔模型的基础上, 引入横向喷射流, 同时考虑了3种尺度流体流动对纵波传播的影响, 建立了含横向喷射流的层状部分饱和模型, 推导了模型的波动方程, 并通过分析平面波, 计算了纵波的相速度和品质因子, 研究了纵波的衰减和速度频散特性。数值模拟发现, 与层状双孔模型相比, 本文建立的模型不仅可以在低频范围内描述纵波的衰减和速度频散现象, 还可以在高频范围内很好地描述纵波的高频散和强衰减现象, 说明含横向喷射流的层状部分饱和模型是一个较为合理的跨尺度理论模型。本文同时分析了部分岩石参数对纵波衰减和速度频散特性的影响, 结果表明, 频散曲线对模型的含气饱和度有较好的指示作用; 随着模型孔隙度的减小, 频散曲线向速度增加的方向移动, 且低频范围内的频散曲线和衰减峰向低频方向移动; 随着喷射长度的增加, 高频范围内的频散曲线和衰减峰向低频方向移动, 低频范围内的频散曲线和衰减峰无明显变化。

附录A

将ZHAO等^[28]推导的动能和耗散函数推广到含横向喷射流的层状部分饱和模型中, 可以写为:

$ 2T = {\rho _{00}}\sum\limits_i {\dot u_x^2} + 2\sum\limits_m {{\rho _{0m}}} \sum\limits_i {{{\dot u}_x}\dot U_x^{\left( m \right)}} + \sum\limits_m {{\rho _{mm}}} \sum\limits_i {{{\left( {\dot U_x^{\left( m \right)}} \right)}^2}} + 2{T_{\rm{L}}} $

(A1a)

$ 2D = \sum\limits_m {{b_m}} \sum\limits_i {{{\left( {\dot U_x^{\left( m \right)} - {{\dot u}_x}} \right)}^2}} + 2{D_{\rm{L}}} $

(A1b)

其中, T_L和D_L分别表示动能和耗散函数的局部流体流动项, 具体表达式为:

$ {T_{\rm{L}}} = \frac{1}{2}\rho _f^{\left( 1 \right)}{\varphi _1}\varphi _2^2X_1^2{{\dot \zeta }^2} + \frac{3}{2}{\rm{ \mathsf{ π} }}{R^2}\rho _f^{\left( 2 \right)}\varphi _2^2{\varphi _{20}}X_1^2\left( {{X_2} - \frac{{{\varphi _1}}}{{{\rm{ \mathsf{ π} }}{R^2}{\varphi _{10}}}}} \right){{\dot \zeta }^2} $

(A2a)

$ {D_{\rm{L}}} = \frac{1}{6}\left( {\frac{{{\eta _1}}}{\kappa }} \right){\varphi _1}{\varphi _{10}}\varphi _2^2X_1^2{{\dot \zeta }^2} + \frac{1}{2}{\rm{ \mathsf{ π} }}{R^2}\left( {\frac{{{\eta ^2}}}{\kappa }} \right)\varphi _{20}^2\varphi _2^2X_1^2\left( {{X_2} - \frac{{{\varphi _1}}}{{{\rm{ \mathsf{ π} }}{R^2}{\varphi _{10}}}}} \right){{\dot \zeta }^2} $

(A2b)

附录B

将(21)式和(23)式代入(3)式, 可以得到:

$ \begin{array}{l} {Z_1}\frac{{{\partial ^2}{u_x}}}{{\partial {x^2}}} + {Z_2}\frac{{{\partial ^2}U_x^{\left( 1 \right)}}}{{\partial {x^2}}} + {Z_3}\frac{{{\partial ^2}U_x^{\left( 2 \right)}}}{{\partial {x^2}}} - {Z_4}\frac{{\partial \zeta }}{{\partial x}}\\ \;\;\;\;\;\; = {\rho _{00}}{{\ddot u}_x} + {\rho _{01}}\ddot U_x^{\left( 1 \right)} + {\rho _{02}}\ddot U_x^{\left( 2 \right)} + {b_1}\left( {{{\dot u}_x} - \dot U_x^{\left( 1 \right)}} \right) + {b_2}\left( {{{\dot u}_x} - \dot U_x^{\left( 2 \right)}} \right) \end{array} $

(B1a)

$ {Z_5}\frac{{{\partial ^2}{u_x}}}{{\partial {x^2}}} + {Z_6}\frac{{{\partial ^2}U_x^{\left( 1 \right)}}}{{\partial {x^2}}} - {Z_7}\frac{{\partial \zeta }}{{\partial x}} = {\rho _{01}}{{\ddot u}_x} + {\rho _{11}}\ddot U_x^{\left( 1 \right)} - {b_1}\left( {{{\dot u}_x} - \dot U_x^{\left( 1 \right)}} \right) $

(B1b)

$ {Z_8}\frac{{{\partial ^2}{u_x}}}{{\partial {x^2}}} + {Z_9}\frac{{{\partial ^2}U_x^{\left( 2 \right)}}}{{\partial {x^2}}} + {Z_{10}}\frac{{\partial \zeta }}{{\partial x}} = {\rho _{02}}{{\ddot u}_x} + {\rho _{22}}\ddot U_x^{\left( 2 \right)} - {b_2}\left( {{{\dot u}_x} - \dot U_x^{\left( 2 \right)}} \right) $

(B1c)

其中,

$ \begin{array}{l} {Z_1} = A + 2N - \frac{{Q_1^2}}{{{R_1}}} - \frac{{Q_2^2}}{{{R_2}}} + \frac{{{\gamma _1}{Q_1}}}{{{R_1}}}{F_1}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _1}R} \right)}}{{{\lambda _1}R{J_0}\left( {{\lambda _1}R} \right)}}} \right]\\ \;\;\;\;\;\;\; + \frac{{{\gamma _2}{Q_2}}}{{{R_2}}}{F_2}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _2}R} \right)}}{{{\lambda _2}R{J_0}\left( {{\lambda _2}R} \right)}}} \right] \end{array} $

(B2a)

$ {Z_2} = \frac{{{\varphi _1}{Q_1}}}{{{R_1}}}{F_1}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _1}R} \right)}}{{{\lambda _1}R{J_0}\left( {{\lambda _1}R} \right)}}} \right] $

(B2b)

$ {Z_3} = \frac{{{\varphi _2}{Q_2}}}{{{R_2}}}{F_2}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _2}R} \right)}}{{{\lambda _2}R{J_0}\left( {{\lambda _2}R} \right)}}} \right] $

(B2c)

$ {Z_4} = {\varphi _1}{\varphi _2}\left\{ {\frac{{{Q_1}}}{{{R_1}}}{F_1}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _1}R} \right)}}{{{\lambda _1}R{J_0}\left( {{\lambda _1}R} \right)}}} \right] - \frac{{{Q_2}}}{{{R_2}}}{F_2}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _2}R} \right)}}{{{\lambda _2}R{J_0}\left( {{\lambda _2}R} \right)}}} \right]} \right\} $

(B2d)

$ {Z_5} = {\gamma _1}{F_1}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _1}R} \right)}}{{{\lambda _1}R{J_0}\left( {{\lambda _1}R} \right)}}} \right] $

(B2e)

$ {Z_6} = {\varphi _1}{F_1}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _1}R} \right)}}{{{\lambda _1}R{J_0}\left( {{\lambda _1}R} \right)}}} \right] $

(B2f)

$ {Z_7} = {\varphi _1}{\varphi _2}{F_1}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _1}R} \right)}}{{{\lambda _1}R{J_0}\left( {{\lambda _1}R} \right)}}} \right] $

(B2g)

$ {Z_8} = {\gamma _2}{F_2}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _2}R} \right)}}{{{\lambda _2}R{J_0}\left( {{\lambda _2}R} \right)}}} \right] $

(B2h)

$ {Z_9} = {\varphi _2}{F_2}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _2}R} \right)}}{{{\lambda _2}R{J_0}\left( {{\lambda _2}R} \right)}}} \right] $

(B2i)

$ {Z_{10}} = {\varphi _1}{\varphi _2}{F_2}\left[ {1 - \frac{{2{J_1}\left( {{\lambda _2}R} \right)}}{{{\lambda _2}R{J_0}\left( {{\lambda _2}R} \right)}}} \right] $

(B2j)

附录C

将(4)式进行整理, 并将整理后的(4)式、(17)式、(18)式、(19)式和(29)式代入(3)式, 可以得到:

$ \left( {\begin{array}{*{20}{c}} {{a_{11}}{k^2} + {b_{11}}}&{{a_{12}}{k^2} + {b_{12}}}&{{a_{13}}{k^2} + {b_{13}}}\\ {{a_{21}}{k^2} + {b_{21}}}&{{a_{22}}{k^2} + {b_{22}}}&{{a_{23}}{k^2} + {b_{23}}}\\ {{a_{31}}{k^2} + {b_{31}}}&{{a_{32}}{k^2} + {b_{32}}}&{{a_{33}}{k^2} + {b_{33}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{u_{10}}}\\ {U_{10}^{\left( 1 \right)}}\\ {U_{10}^{\left( 2 \right)}} \end{array}} \right) = 0 $

(C1)

其中,

$ \begin{array}{l} {a_{11}} = {Z_1} - {Z_4}\frac{{{\varphi _1}{Q_2} - {\varphi _2}{Q_1}}}{S},{a_{12}} = {Z_2} + \frac{{{Z_4}{\varphi _2}{R_1}}}{S},{a_{13}} = {Z_3} - \frac{{{Z_4}{\varphi _1}{R_2}}}{S}\\ {a_{21}} = {Z_5} - {Z_7}\frac{{{\varphi _1}{Q_2} - {\varphi _2}{Q_1}}}{S},{a_{22}} = {Z_6} + \frac{{{Z_7}{\varphi _2}{R_1}}}{S},{a_{23}} = - \frac{{{Z_7}{\varphi _1}{R_2}}}{S}\\ {a_{31}} = {Z_8} + {Z_{10}}\frac{{{\varphi _1}{Q_2} - {\varphi _2}{Q_1}}}{S},{a_{32}} = - \frac{{{Z_{10}}{\varphi _2}{R_1}}}{S},{a_{33}} = {Z_9} + \frac{{{Z_{10}}{\varphi _1}{R_2}}}{S} \end{array} $

(C2a)

$ \begin{array}{l} {b_{11}} = - {\rho _{00}}{\omega ^2} + {\rm{i}}\omega \left( {{b_1} + {b_2}} \right),\;\;\;{b_{12}} = - {\rho _{01}}{\omega ^2} - {\rm{i}}\omega {b_1},\;\;\;{b_{13}} = - {\rho _{02}}{\omega ^2} - {\rm{i}}\omega {b_2}\\ {b_{21}} = {b_{12}},{b_{22}} = - {\rho _{11}}{\omega ^2} + {\rm{i}}\omega {b_1},\;\;\;\;{b_{23}} = 0\\ {b_{31}} = {b_{13}},{b_{32}} = {b_{23}},{b_{33}} = - {\rho _{22}}{\omega ^2} + {\rm{i}}\omega {b_2} \end{array} $

(C2b)

(C1) 式有解的条件为系数矩阵的行列式为0, 即:

$ \left| {\begin{array}{*{20}{c}} {{a_{11}}{k^2} + {b_{11}}}&{{a_{12}}{k^2} + {b_{12}}}&{{a_{13}}{k^2} + {b_{13}}}\\ {{a_{21}}{k^2} + {b_{21}}}&{{a_{22}}{k^2} + {b_{22}}}&{{a_{23}}{k^2} + {b_{23}}}\\ {{a_{31}}{k^2} + {b_{31}}}&{{a_{32}}{k^2} + {b_{32}}}&{{a_{33}}{k^2} + {b_{33}}} \end{array}} \right| = 0 $

(C3)

求解(C3)式, 可以得到3类纵波的波数解。