0 引言

一直以来饱和流体多孔介质中的弹性波传播问题是工程界和力学界的重要研究课题，Biot(1956a, b)首次考虑流固耦合作用，建立了低频和高频弹性波传播流固耦合双相介质模型，即著名的Biot多孔介质理论，用以预测均匀各向同性岩石中快纵波、慢纵波和剪切波的传播规律.此后，Biot(1962)引入黏性项，建立了最早的黏弹双相介质弹性波模型.基于Biot理论，Brutsaert(1964)、Santos等(1990, 2004)和蔡袁强等(2006)等先后建立了描述饱和两非混相流体多孔介质中弹性波传播规律模型，对弹性波传播和衰减机制进行宏观解释.近60年来，Biot理论得到不断发展，在未固结、高渗沉积岩等各类饱和流体多孔介质中亦得到应用.工程应用中，由于Biot弹性波理论对弹性波衰减和频散的预测值与实际值差距较大，为进一步揭示造成该现象的原因，人们逐渐转向研究微观机制模型和介观尺度模型.许多学者认为微观局部喷射流是引起弹性波衰减和频散的重要因素，但研究发现喷射流机制受到单个孔隙结构的影响，并不适用于复杂多孔介质分析(Mavko and Jizba, 1991，1994；Gurevich et al., 2009).由于喷射流模型和宏观模型都无法单独解释复杂孔隙介质中弹性波的耗散机制，因此Dvorkin和Nur(1993)结合Biot理论和喷射流机制，提出了BISQ模型.朱建伟等(2001)建立起含油水两相流体多孔介质中的弹性波传播模型，将固体位移、流体的流动压力与孔隙度、渗透率、饱和度和密度相联系.Wei和Muraleetharan(2002)建立的模型中考虑了动态相容性条件，提出了第三类纵波的概念.针对介观尺度的弹性波传播模型研究，White(1975)最先提出球状斑块饱和模型，模拟了介观尺度下含气水多孔介质中弹性波的传播，此后，诸多研究人员(Norris，1993；Johnson，2001；Pride et al., 2004；Ciz and Gurevich, 2005；Ba et al., 2011；Sun et al., 2015)运用各种方法对White模型进行了扩展和修正，研究了不同的流体斑块模型中的渗流状态，得到了关于弹性波的衰减和波速度耗散的主要机理.Ba等(2016, 2017)分析了黏土射流强化岩石基质引起的纵波频散，研究了双重孔隙介质弹性波传播，由于片状饱和、结构非均质等因素表现出的岩石滞弹性.以上模型中饱和流体均处于静止状态，该假设广泛应用于地震勘探和地球物理学等领域.

应用于油田开发过程中的低频人工地震波辅助强化采油技术亦涉及弹性波的传播问题，由于储层流体在注采压差下始终处于渗流状态，关于该条件下波场与渗流场耦合作用的波动耦合渗流力学问题(刘静等, 2014a, b ；郑黎明等, 2015, 2016, 2017)，与目前饱和静态流体多孔介质中弹性波传播理论存在明显区别，直接利用静态弹性波传播理论解释其动力学原理与设计、优化工艺参数，会出现与现场工程技术应用背景不统一的问题.因此，需进一步考虑储层渗流场，对静态模型进行改进和发展.

为了建立反映弹性波场与流固耦合渗流场相互作用的控制方程，本文基于等效流体理论和Biot理论假设条件，考虑油水两相渗流、流固压缩性、初始条件和边界条件，建立了各向同性、均匀油水两相渗流储层多孔介质弹性波传播动力学模型，通过引入宏观渗流和位移势函数，对渗流介质中纵/横波波速和品质因子进行求解，分析了考虑与忽略动态饱和度梯度两种情形下弹性波传播的物理意义，研究成果对深化低频人工地震波辅助强化采油技术作用机理认识和优化工艺参数具有重要意义.

1 模型假设与弹性波方程建立 1.1 假设条件

(1) 储层为均匀、各向同性多孔介质；(2)流体和骨架可压缩；(3)频率低于特征频率时，弹性波引起的诱导流动呈poiseuille型流动；(4)等温渗流过程，忽略热耗散作用(蒲春生等，2017)；(5)油藏流体初始时刻处于渗流状态；(6)油、水两相互不混相，用体积平均物性反映孔隙流体整体物性；(7)忽略重力作用对渗流的影响.

目前饱和静态流体多孔介质中弹性波传播理论主要考虑恒定流体饱和情形，油藏条件和储层物性保持不变；而饱和油水两相渗流储层多孔介质弹性波传播理论模型，由于用于储层开发过程中，油藏条件和储层物性随开发过程发生变化.因此，模型建立考虑了饱和静态流体多孔介质中弹性波传播理论与多相渗流力学的耦合，将渗流力学辅助方程和迭代变化的耦合物性方程，在饱和静态流体多孔介质中弹性波传播理论运动方程、连续性方程中同时考虑迭代变化的耦合物性参数(饱和度、压力、孔隙度、导流能力等的变化均不忽略)引入饱和静态流体多孔介质中弹性波传播理论模型.弹性波场求解过程中，引入表示初始状态流体流动的一维线性势函数，将其与谐函数形式的波诱导位移势函数耦合.

1.2 弹性波传播动力学模型建立

(1) 运动方程

多孔介质中含油水两相流体时，用下标w表示水相，o表示油相，系统运动方程如式(1)所示；受储层开发压力变化和渗流场影响，流体饱和度、孔隙介质导流能力、流体密度与固体密度均为变量.

(1)

式中u、U、W分别为固体位移、水相相对位移和油相相对位移；根据Dupuit-Forchheimer渗流公式，令水相、油相分别相对于固体骨架的渗流速度为分别为水相和油相的真实位移，U=U_ff+U_fs，W=W_ff+W_fs；U_ff为水相初位移；U_fs为水相波诱导位移；W_ff为油相初位移；W_fs为油相波诱导位移；τ_ij为作用于孔隙介质的有效应力；p_w为水相压力；p_o为油相压力；ρ_w、ρ_o、ρ为水相密度、油相密度和多孔介质平均密度，ρ=(1－ϕ)ρ_g+ϕs_wρ_w+ϕs_oρ_o(余俊等，2014)，ρ_g为固体颗粒密度；η_w、η_o分别为水相黏度系数和油相黏度系数，视为定值；k_w、k_o分别为水和油的相渗透率；s_w、s_o分别为水相饱和度和油相饱和度；ϕ为孔隙度.

(2) 物性方程

油气藏开发过程中，储层压力变化较大，地下流体体积和储层孔隙体积受油水压缩性的影响不应被忽略.文中用体积模量来表示孔隙流体的体积应变.孔隙流体体积应变和岩石骨架、颗粒固体体积应变分别为

(2)

式中K_w、K_o、K_s、K_g为水相体积模量、油相体积模量、骨架体积模量和颗粒体积模量；p′为作用在骨架上的平均有效压力；p是孔隙压力，p=s_wp_w+s_op_o(体积平均法)，dp=s_wdp_w+s_odp_o－p_cds_w；p_c为油水相间毛管力；ρ_d为固体骨架密度，ρ_d=(1－ϕ)ρ_g.

对低渗透多孔介质(如岩石)的有效压力进行修正(杜修力等，2018；Bishop，1959)，表示为式(3)，

(3)

式中λ和μ为Lame常数；α为等效孔隙压力系数，α=1－K_s/K_g∈(0, 1]，τ′_ij为总应力，ε_v为固体骨架体应变，ε_v=∇·u；δ为Kronecker符号，δ=0(i≠j)或1(i=j).

(3) 状态方程

受开发储层场影响，运动方程中需考虑渗透率、油水相间毛管力随饱和度的变化.引入描述流体相对渗透率的Mualem模型，当非润湿相为油时，相对渗透率与饱和度关系表示为式(4)所示(Lo et al., 2005)：

(4)

毛管力曲线表示为式(5)所示：

(5)

式中η、n、χ、m为毛管力曲线拟合参数，m=1-1/n.

(4) 连续性方程

本文所研究的弹性波波动引起的流动属于非稳态无源(汇)流动，在多孔介质空间内任意取一个有界体积；根据质量守恒定律，可得到含油水两相渗流流体多孔介质中的连续方程如式(6)所示(孔祥言，2010).

(6)

(5) 饱和非混相渗流流体多孔介质弹性波方程

展开式(6)并舍去二阶小量u·∇，联立其他基础方程，推导得到饱和两非混相渗流流体多孔介质中的波动方程，见式(7)，考虑到波动采油对渗流速度和孔隙物性有较大改善，因此连续性方程推导过程中未考虑孔隙度与压力的变化.

其中，参数.

2 弹性波传播动力学模型求解

在均匀各向同性多孔介质中，纵波与横波独立存在，不存在耦合，遵从独立的传播规律，此时可将速度按照纵波(P)和横波(S)进行分解.本文引入了宏观渗流，将水相和油相的位移表示为初始位移项(U_ff、W_ff)和波诱导位移项(U_fs、W_fs)的形式.对固体和流体位移进行Helmholtz分解，表示为标量势的形式(姜领发等，2007)：

式中下标P和S表示纵波和横波，φ是标量，ψ是矢量；下标ff和fs表示初始状态和波诱导状态.

将上述标量势代入波动方程(7)，并用三角函数表示波场下的流固位移势：

式中ω为角频率，r/s；l_P、l_S分别为P波和S波的波数；r为波的位置矢量；v_of和v_wf分别为油相和水相的初始渗流速度，v_of=v_f(1－f_w)，v_wf=v_ff_w，v_f为初始渗流速度；A_i和B_i为系数，且有：

2.1 忽略grads_w项情形求解

当原始油藏中渗流速度非常小(但不等于0)时，近似于饱和静态流体的情况，对一维均匀各向同性油藏而言，可认为grads_w=0.适用于油藏开发初期(图 1).

将平面波解代入波动方程中，整理可得到纵波和横波的波矢量方程，分别如式(8)和式(9)所示.

对式(8)和式(9)进行整理，得到纵波波矢量方程为

(10)

得到横波波矢量方程为

(11)

求解以上波矢量特征方程，由于此模型中有限制条件Re(l_j)>0和Im(l_j) < 0(Re和Im分别为复数的实部和虚部)，所以可得到纵波波数l_P1、l_P2、l_P3和横波波数l_S，进而根据式(12)求得四种体波的传播速度v_j及品质因子Q_j：

(12)

2.2 考虑grads_w项情形求解

当油藏饱和渗流流体时，必须考虑饱和度在各方向上的变化，即grads_w≠0.以一维均匀活塞驱替为例(图 2)，水驱前缘位置处grads_w∈(0, ∞)，其中s_w=1－s_o；水区：grads_w=0且s_w=1－s_or；油区：grads_w=0且s_w=s_wi；适用于和油藏开发过程匹配的波动采油技术.

(1) 纵波波矢量方程

(13)

式(13)用系数性方程表示如下：

(14)

由以上矩阵可知，当波速l_P已知时，系数矩阵为已知，由此可求出势系数A₁、A₂、A₃，如果势系数A₁、A₂、A₃已知时，可以求出流固位移u、W、U；当l_P未知时，矩阵C_i包含原始渗流的影响，是A₁、A₂、A₃的函数.

(2) 横波波矢量方程

计算发现，考虑grads_w时的横波波矢量方程与不考虑grads_w时的横波波矢量方程相同.

3 算例分析与对比 3.1 算例分析

为便于与前人理论进行对比，流固物性参数的选取与Lo(2005)的Columbia砂岩保持一致，如表 1所示.

本文数学模型中给出了渗流存在条件下，考虑grads_w=0和grads_w≠0时多孔介质中波数的求解矩阵，给出了grads_w=0的情形下的实际算例，适用于油田开发前期油藏内部稳定渗流阶段弹性波传播规律分析.

求解式(9)、(11)可得到4种体波(第一类纵波P₁，第二类纵波P₂，第三类纵波P₃和横波S)的波矢量值.进而根据式(14)求得波速和品质因子.4种体波的产生机理为：P₁波的传播由岩石骨架和流体的耦合运动产生，P₂波传播由岩石骨架和两相孔隙流体的反相位运动产生，P₃波传播由于毛管力波动而产生，S波沿液流垂直方向作剪切运动.

(1) 波速变化规律分析

由图 3可知：(1)s_w一定时，纵波和横波的波速均随振动频率的增加而增大，符合波的传播定理；(2)频率一定，随着含水饱和度的增加，纵波波速先逐渐增大后趋于平稳，而横波波速随饱和度增加而增大的幅度较小.

分析认为，随s_w的增加，孔隙流体密度随之增加，系统压缩系数减小，体积模量增加，而且密度变化引起的体积模量变化要比密度本身变化大得多(水相饱和度为40%时，岩石平均密度为1.815 g·cm^-3，岩石体积模量M为1.34 GPa；水相饱和度为60%时，岩石平均密度为1.837 g·cm^-3，岩石体积模量M变化为1.39 GPa)，因此，纵波的波速随s_w的增大而增大；与之相比，由于多孔介质剪切模量受孔隙流体影响甚微，所以横波的波速随含水饱和度变化较小，当含水饱和度大于80%时，横波波速略微增加，基于拟合相渗曲线下此时水相占优，油相接近残余油状态，此时由饱和双相流体情形弹性波传播趋近于饱和水相弹性波传播规律.

(2) 品质因子的变化规律分析

图 4(a-b)中，P₁波品质因子随含水饱和度s_w增大先缓慢减小后小幅增加，P₂波的品质因子随含水饱和度s_w增大而增加，其中P₂波品质因子变化幅度大于P₁波品质因子.

图 4c中，P₃波品质因子随s_w的增大先降低后趋于稳定.图 4d中，不同频率下S波品质因子均随s_w增加而降低.分析造成图 4c与图 4(a-b)三类纵波品质因子不同变化趋势的原因为：P₁和P₂波是流固相惯性耦合的结果，而P₃波认为是两流体相惯性耦合产生，产生机理的不同导致波衰减规律不同.且实验也已证明，在低振动频率下，只有P₁波能明显观测到，P₂、P₃和横波相对能量较弱，衰减严重，较难测量到(杜艺可，2012；丁卫等，2014).除弹性波产生原理和常规传播规律影响，在双相渗流条件下含水饱和度的变化会影响整体渗透率(k_w+k_o)的数值，弹性波波速和衰减规律会进一步受到影响.

通过以上算例分析，可分别揭示低频人工地震波辅助强化采油技术中波速、品质因子受频率和含水饱和度的影响规律.

(3) 模型验证

本模型考虑孔隙度随时间变化、忽略grads_w的变化时，P₁波波速随饱和度的增加而增加.对模型进行简化，忽略孔隙度随时间的变化和grads_w的变化，模型可退化为类似于饱和静态流体弹性波传播模型的变化趋势(图 5所示)，参考Murphy无因次波速变化实验结果(Fort Union sandstone，孔隙度8.5%，频率4~6 Hz)进行拟合(Murphy Ⅲ，1982)，当含水饱和度s_w较小时，P1波波速呈下凸型减小趋势，当含水饱和度s_w增加(s_w大于30%时)，波速呈逐渐增加的趋势，本文模型计算值(s_w处于30%及以上)与实验结果基本吻合.

3.2 宏观渗流场对弹性波传播影响特征分析

(1) 宏观渗流场对波速影响规律分析

由图 6可知：频率为50 Hz时，宏观渗流的存在使P₁波的波速下降了10~10²m·s^-1左右，使P₂波和P₃波的波速大幅增加，S波波速小于静态模型，随s_w的增加逐渐趋近于静态模型值.宏观渗流的存在影响固、液相间的交错运动，根据不同种类波的产生机制，流固交错运动程度变化必然对诱导波的性质产生影响，当主要宏观渗流方向与波诱导液体运动方向一致时，本模型中P₁波的波速相对静态模型中P₁波的波速减弱，而P₂和P₃波波速却相对增强，S波波速由于垂直于宏观渗流方向而受交错运动变化影响小.由于储层开发中开展的波场检测与分析，主要采用快纵波，文中关于P₁波波速变化对于开发储层物性参数分析具有较强指导意义，关于P₂和P₃波在开发储层中的应用需进一步探索.

(2) 宏观渗流场对品质因子影响规律分析

由图 7可知，本模型中各纵波和横波的品质因子相对于饱和两相静态流体模型结果均有很大差异，P₁波的品质因子与静态模型值在同一数量级，P₂波、P₃波和S波的衰减程度相对于饱和静态流体模型更加显著，数量级增加了10²左右，但当含水饱和度逐渐增大时，考虑宏观渗流饱和双相流体模型中弹性波波速数值与静态情形弹性波波速数值差距减小.计算结果差值反映油藏渗流场对弹性波传播的影响程度，初始渗流场的存在增强了油水两相流体的黏滞作用，不同流相与固相的交错运动反复受到激励与抑制，且由于模型假设孔隙可压缩，波诱导液体流动在垂直渗流方向发生射流机制，平行渗流方向的弹性诱发流动被削弱，因此纵、横波品质因子均相对饱和两相静态流体模型结果发生变化.

4 结论

(1) 本文基于等效流体理论和Biot理论基本原理，引入宏观渗流场和位移势函数，建立了适用于油藏开发过程的油水两相渗流条件下储层多孔介质弹性波传播动力学模型.通过算例求解与分析，发现含油水两相渗流储层多孔介质中同时存在着3种纵波P₁、P₂、P₃和1种横波S，为揭示低频人工地震波辅助强化采油技术动力学机理和优化设计工艺参数奠定了重要的理论基础.

(2) 纵、横波波速和品质因子受到频率和储层含水饱和度s_w的影响呈现不同变化规律，当含水饱和度s_w一定时，纵波和横波的波速都随着振动频率的增加而增大，频率一定时，纵、横波的波速均随含水饱和度s_w的增加而增大.第一纵波P₁、第二纵波P₂的品质因子随频率和含水饱和度s_w增大而增强，第三纵波P₃品质因子随含水饱和度s_w的增大先降低后逐渐趋于平稳，横波S品质因子随含水饱和度s_w的增大而降低.

(3) 初始渗流场的存在，增强了油水两相流体的黏滞作用，使本模型弹性波传播规律与统静态模型存在明显差异，不同相态间的交错运动受到反复激励与抑制，其中渗流场对“流-固”系统的正向弹性运动以抑制作用为主，对其反向弹性运动以促进作用为主，两相流体的惯性作用交错增强.

附录

式(10)、(11)方程中的系数定义如下：