孔隙内填充单一固体的固-固孔隙介质中的声波传播*

2022-05-26刘琳张秀梅王秀明

物理学报 2022年9期

刘琳张秀梅† 王秀明

1)(中国科学院声学研究所,声场声信息国家重点实验室,北京 100190)

2)(中国科学院大学,北京 100149)

3)(中国科学院声学研究所,北京市海洋深部钻探研究中心,北京 100190)

连通孔隙空间被有别于骨架固体的另一种固相介质充填而形成的双组分连通固体孔隙介质,简称固-固孔隙介质.推导出了固-固孔隙介质的声波动力学方程和本构关系,利用平面波分析的方法研究了各种波的频散和衰减特性.在此基础上,提出了基于一阶速度-应力方程的时间分裂的高阶交错网格有限差分算法,对其中的声场演化特点进行了模拟计算,分析了各种波的产生机制和能量分布,并详细讨论了两种固体间的摩擦系数和声源频率对各种波传播特性的影响.数值模拟表明:固-固孔隙介质中存在两种纵波(P1 和P2)和两种横波(S1 和S2),其中P1 和S1 波能量主要在骨架固体中传播;P2 和S2 波是骨架固体和孔隙固体之间相对运动产生的慢波,能量主要在孔隙固体中传播.固体骨架和孔隙固体之间的摩擦主要影响慢波(P2和S2)的衰减,且低频时衰减大于高频.

1 引言

Biot[1−4]于二十世纪五六十年代首次提出流体饱和孔隙介质的弹性波传播理论,奠定了孔隙介质声学的基础.Biot 理论预言了在流体饱和孔隙介质中存在着两类纵波(快纵波和慢纵波)和一类横波.Plona[5]从实验室中观测到Biot 理论中预言的慢纵波,从而推动了对Biot 理论的进一步研究和应用[6−8].自Biot 理论建立以来,多相孔隙介质理论及其波动问题得到了广泛而持续的研究[9−13].当流-固孔隙介质中的流体相被另一种固体替代时,固体填满整个孔隙空间,就形成了固-固孔隙介质,例如,在能源勘探开发领域,地下石油与天然气储层是典型的由固体颗粒骨架和内部孔隙空间组成的多孔介质,储层中悬浮着可随孔隙流体自由运移的颗粒,并会发生聚集,进一步可能会完全填充原有的孔隙空间形成固-固孔隙介质[14,15];在天然气水合物的开采中,孔隙中完全填充天然气水合物的水合物储层也可以看成是固-固孔隙介质.其次,在声学压电材料的研究中,由压电材料及聚合物构成的双固体压电复合材料可以提高材料的机械性能以及电性能,其构成包含两种规则或者不规则分布的孔隙压电介质[16,17].此外,在医学领域,骨质骨中的骨松质结构疏松,孔隙中存在骨髓,可以看成是具有孔隙结构的介质,孔隙中的骨髓应该看作是剪切模量很小的软组织.为此,研究固-固孔隙介质中的声波传播机理和规律具有重要意义.

自Biot 理论建立以来,国内外学者对孔隙介质中波的传播规律的数值模拟研究有了较大进展.Zhu 和McMechan[18]利用以位移表述的有限差分算法模拟了非均匀孔隙弹性介质中声波的传播.Dai 等[19]提出了用于模拟孔隙介质声波传播的速度-应力有限差分算法.Carcione 等[20]利用解析式分析了孔隙度、渗透率、黏滞性等对声波的影响特点.Atalla 等[21,22]利用基于Biot 理论的有限元法模拟了孔隙介质声波的传播.杜启振等[23]采用伪谱法正演模拟了各向异性黏弹性孔隙介质地震波,并对其波场特征进行了分析.王秀明等[6]利用高阶交错网格有限差分算法模拟了波在非均匀孔隙介质中的传播.在此基础上,赵海波等[24]基于Santos理论[25]发展了方程存在刚性部分情况下的时间分裂[20]的高阶交错网格有限差分法,解决了方程存在刚性的问题,清楚地展示了不饱和孔隙介质中存在的三类纵波和一类横波.为了研究孔隙中填充一种固体和一种流体的孔隙介质中弹性波传播规律,Carcione 和Seriani[26]利用伪谱法模拟出冷冻多孔介质(冻土或永冻土)的波场,指出该孔隙介质中存在三种纵波和两种横波,但其数值计算的声场存在一定程度的频散,不能清晰地展示每种波的传播和相互作用情况.Gao 等[27,28]使用变阶数时域有限差分模拟孔隙弹性方程.Zhan 等[29]使用间断有限元模拟孔隙弹性波动方程.刘琳等[30]基于Carcione 理论[26]利用时间分裂的高阶交错网格有限差分法模拟了弹性波在孔隙中填充一种固体和一种流体的孔隙介质中的传播,清楚地显示出该类孔隙介质中存在的5 种波.以上工作表明,针对多相孔隙介质中声波传播问题目前已经开展了较多工作,然而对于固-固孔隙介质中的声波传播规律认识较少.因此,开展声波在固-固孔隙介质中的传播研究,是进一步进行相关的实际应用的物理基础.

本文首先根据动能密度、势能密度和耗散势函数,结合Lagrange 方程得到固-固孔隙介质动力学方程和本构关系.在此基础上,利用平面波分析的方法研究了各种波的频散和衰减特性.研究声波在孔隙介质中的传播过程,必须进行数值模拟,有限差分法是进行声波数值模拟的有效手段.孔隙介质中各部分之间的相对运动会产生一定的摩擦.由于摩擦系数的存在,其运动方程传播矩阵的特征值差别很大,在此情况下,方程被称之为刚性方程.可以利用时间分裂法解决方程的刚性问题[20,24,30,31].因此,利用时间分裂的交错网格有限差分法可以实现对孔隙介质波场的高精度数值模拟.本文根据一阶速度-应力方程,构建了时间分裂的高阶交错网格有限差分算法.以孔隙中完全填充天然气水合物的水合物储层为例,对该类介质中的波场进行了高精度的数值模拟,得到了波场快照图以及质点振动速度时间变化曲线,分析了各类波的能量分布与传播特征.基于此,研究了骨架固体和填充固体之间的摩擦系数和频率对声波传播特征的影响.

2 固-固孔隙介质声传播理论

当流体饱和孔隙介质中的流体相被另一种固体替代时,就形成了由两种固体组成的固-固孔隙介质,且固体骨架和孔隙固体之间存在相对运动.假设孔隙均匀分布,其尺度远小于声波波长.根据Biot 理论[1−4],孔隙介质的一阶速度-应力方程可根据其动能密度、势能密度以及耗散能密度结合Lagrange 定理得到.与Biot 理论描述的双相介质中的动能密度[1−4]类似,孔隙中填充一种固体的固-固孔隙介质的动能密度可以表示为

式中,上角标1 和2 分别表示骨架固体和填充固体,ρ11和ρ22分别为骨架固体和填充固体的有效质量密度,ρ12为两种固体的耦合质量密度,和分别表示骨架固体和填充固体的质点振动位移分量.

骨架固体和填充固体间的相对运动会导致能量的耗散,耗散力表现为两种固体之间的相互摩擦作用力,可以近似成运动速度差的一次函数.耗散力是一种特殊的非保守力,可由显含速度的势函数给出[32].对于各向同性介质,固-固孔隙介质中的耗散势D可以用Biot 理论中对饱和流体孔隙介质中耗散函数的描述来表示[1−4]:

式中,b12为两种固体之间的摩擦系数.根据流体饱和孔隙介质中骨架固体与孔隙流体之间的摩擦系数以及Berryman 和Wang[33]对阻力系数的定义,Guerin 和Goldberg[34]给出了两固一流孔隙介质中固体骨架和孔隙固体间的摩擦系数.在此基础上,定义固-固孔隙介质中两种固体之间的摩擦系数为

固-固孔隙介质的势能密度V只与固体骨架的应变张量和孔隙固体的应变张量有关.忽略二阶以上的项,其Taylor 展开可以表示为

由固体骨架势能密度、孔隙固体势能密度以及二者耦合的势能密度项构成,其中,m=1,2.Leclaire 等[32]对孔隙中填充一种固体和一种流体的孔隙介质势能密度的描述忽略了固体骨架与孔隙固体之间的接触.考虑到固体骨架与孔隙固体之间的接触,可以得到固-固孔隙介质的势能密度为

式中,K1和K2分别是骨架固体和填充固体的体积模量,µ1和µ2分别是骨架固体和填充固体的剪切模量,C12是骨架固体和填充固体的弹性耦合系数,µ12是由于骨架固体和填充固体耦合形成的剪切模量,θm和dm为应变张量的不变量.θm=,,其中δ是克罗内克符号.由于应力分量可以表示为

因此可以得到固-固孔隙介质的本构方程为

固-固孔隙介质的运动方程可以由Hamilton原理得到.对于保守系统的Lagrange 密度可以表示为

可以用基于Hamilton 最小作用原理的Lagrange方程来描述保守系统的运动.带耗散势的以位移为广义坐标的Lagrange 方程[1,9,32]可以写成:

式中上标m=1,2 分别表示骨架固体和填充固体.由于动能密度是速度的函数,势能密度是应变的函数,因此可以得到

将动能密度、势能密度和耗散势函数代入Lagrange方程中可以得到固-固孔隙介质的运动方程为

将(7a)—(7f)式代入到(11a)式和(11b)式中,可以得到

式中,ρ,b,R和µ分别为质量密度、摩擦系数、刚度和剪切模量矩阵:

其中,

这里,c1和c2分别代表骨架固体和孔隙固体骨架的固结系数,具体表示见附录A .将(12)式分别进行梯度和散度运算,得到纵横波的运动方程为

纵波和横波的位移势函数分别为

式中,A(1)和A(2)分别为固体骨架和孔隙固体纵波势函数振幅分量,B(1)和B(2)分别为固体骨架和孔隙固体横波势函数振幅分量,kP和kS分别表示纵波和横波波矢量.将纵波和横波的位移势函数分别代入纵横波运动方程(15)中得到

由于固体骨架和孔隙固体势函数振幅存在非零解,因此可以得到纵波和横波的频散方程为

纵横波的传播系数和衰减因子分别为

用来评价衰减的逆品质因子可以表示为

式中j=P1,P2,S1 和S2.

3 数值算法

本文采用基于时间分裂的高阶交错网格有限差分法实现对固-固孔隙介质波场的模拟.分裂法[20,24,30,31]是将运动方程分为刚性和非刚性两部分,刚性部分可以得到解析解,将刚性部分得到的解析解作为初始值输入到方程的非刚性部分,再使用交错网格有限差分算法求解.

将方程(11a)和(11b)两边对时间求导,在二维情况下为

将运动方程(11a)和(11b)写成如下矩阵形式:

式中,

这样得到了如下速度-应力关系的运动平衡方程,该表述方式有利于之后的数值模拟算法的建立,即:

式中γij和Π=ij分别为矩阵γ和Π中元素.

(24a)—(24f)式,(26a)式和(26b)式构成了两种固体组成的固-固孔隙介质中弹性波在二维平面内传播的一阶速度-应力方程组.

将方程(25)的刚性部分写成如下矩阵形式:

其中S=.在二维平面内,刚性方程(25)的解析解为

可利用Sylvester 公式求取 exp(SΔt),即

式中λ为矩阵S的特征值,λ=b12(γ12−γ11+γ21−γ22).

方程(25)中非刚性部分写成如下矩阵形式:

式中B=.

定义中间过渡向量V′=,作为初始量输入到方程非刚性部分.过渡向量中质点振动速度为刚性部分解析解.二维情况下非刚性部分的交错网格有限差分形式为(以固体颗粒骨架为例)

式中Dx和Dz表示离散化的微分算子,以x方向的差分算子Dx作用于为例,表示为

其中L为差分算子的长度,本文中L=4,即采用空间阶数为8 的高阶交错网格有限差分算法.二维情况下各场量在交错网格有限差分法中的位置分布如图1 所示.

图1 各场量在交错网格有限差分法中的位置分布Fig.1.Distribution of the field components in the staggeredgrid finite-difference algorithm.

4 数值模拟

本节以孔隙中完全填充天然气水合物的水合物储层为例,首先求解固-固孔隙介质频散方程,分析各种波的相速度和衰减与频率的关系.其次,将时间分裂的有限差分算法运用到孔隙中填充固体的固-固孔隙介质中的声波数值模拟中,讨论了各种波的产生机制,分析了多孔介质中波的能量分布与传播特征,研究了两种固体之间的摩擦系数和频率对声波传播特性的影响.

4.1 传播特性分析

现针对一个均匀的孔隙中充满固体的固-固孔隙介质模型(以孔隙中完全填充天然气水合物的水合物储层为例),两种固体之间存在耦合.模型的物性参数见表1[34].

表1 模型物性参数表Table 1. Physical parameters of different phases.

通过求解纵横波频散方程(19)和(20)得到固-固孔隙介质中存在两种纵波和两种横波:第一类纵波(P1)、第二类纵波(P2)、第一类横波(S1)和第二类横波(S2).在此基础上进一步计算出各种波相速度和衰减.图2(a)—(d)分析了体波的相速度和衰减与频率的关系.从图2(a)和图2(b)可以看出,P1 和S1 波的速度随频率增大有微小的变化,这一变化是由固体骨架和孔隙固体之间的摩擦系数引起的,且影响程度与频率有关.P1 和S1 波的衰减较小,随频率的增大先增大后减小,存在衰减峰.从图2(c)和图2(d)可以发现,P2 和S2 波的速度在低频时随频率的增大而增大,其衰减远大于两种快波的衰减.另外,由于两种慢波的衰减主要受固体骨架和孔隙固体之间的摩擦影响,因此两种慢波的衰减随频率的增大而减小.结合频散方程可以发现,固-固孔隙介质体波传播特性与两种固体之间的摩擦系数参考值和频率之比有关,且两种慢波高频时衰减比低频情况下小.

图2 相速度和衰减随频率的变化曲线 (φ=0.3,=2.2×108)(a) P1 波;(b) S1 波;(c) P2 波;(d) S2波Fig.2.Effects of frequency on the phase velocities and attenuation of waves (φ=0.3,=2.2×108):(a) P1;(b) S1;(c) P2;(d) S2.

为了进一步研究各种波的传播特性,以天然气水合物储层为例,利用时间分裂的有限差分法模拟了波在固-固孔隙介质中的传播.此时假设两种固体之间的摩擦系数b12=0 .在二维x-z平面内,取x轴正方向水平向右,z轴正方向垂直向下.选用模型网格大小为1600×1600,模型x方向和z方向空间步长均为2 m,时间步长为0.4×10–3s.模型的物性参数如表1[34]所列.声源位于模型中心位置(1600 m,1600 m)处.数值模拟时选用雷克子波作为声源函数,其时域形式为

式中,f0为声源主频,f0=20 Hz.声源能量按比例分配到骨架和孔隙固体质点振动速度水平方向分量上,具体形式为

图3 显示了0.45 s 时刻的质点振动速度水平和垂直分量的波场快照图,图3(a)—(d)的幅度显示比例为1∶1,分别对应于骨架固体质点振动速度水平分量、填充固体质点振动速度水平分量、骨架固体质点振动速度垂直分量和填充固体质点振动速度垂直分量质点的情况.在图中可以清晰地观察到4 种体波,从外到内依次为P1,S1,P2 和S2 波.为了进一步分析各种波的产生机制,图4 分别计算和分析了骨架固体和孔隙固体的质点振动速度水平分量归一化的波形.模型以及参数与图3 选用模型和参数一致,接收点位于(200 m,2600 m)位置处.其中蓝色实线和红色虚线分别表示骨架固体和填充固体传播的波形.从图4 可以观察到,P1 和S1 波在骨架固体质点振动速度波形中幅度最大,骨架固体和填充固体质点运动同相.P2 和S2 波在填充固体质点振动速度波形中幅度最大,骨架固体和填充固体质点运动反相.类似于Biot理论中由于骨架固体与孔隙流体相对运动产生的慢纵波,P2 波和S2 波分别是由于骨架固体和填充固体之间相对运动产生的慢纵波和慢横波.

图3 频率为20 Hz 且 b12=0 时,0.45 s 时刻质点振动速度水平及垂直分量波场快照 (a) 骨架固体质点振动速度水平分量;(b) 填充固体质点振动速度水平分量;(c) 骨架固体质点振动速度垂直分量;(d) 填充固体质点振动速度垂直分量Fig.3.Snapshots of the (a),(b) horizontal and (c),(d) vertical particle velocity component of (a),(c) solid grain frame and (b),(d) pore solid when the frequency is 20 Hz and b12=0 at 0.45 s.

图4 频率为20 Hz 且 b12=0 时,各相质点振动速度水平分量的波形Fig.4.Waveform comparisons of the horizontal particle velocity component between different phases when the frequency is 20 Hz and b12=0 (normalized).

4.2 参数影响

多相孔隙介质中孔隙中固体或流体与固体颗粒骨架的摩擦会使慢波发生衰减.对于固-固孔隙介质,分析了两种固体间的摩擦系数对各种波传播的影响.图5(a)和图5(b)分别为0.45 s 时刻两种介质之间存在摩擦时(=2.2×108)骨架固体和填充固体的质点振动速度水平分量的波场快照,图5(a)和图5(b)的幅度显示比例为1∶1.图6 显示了两种介质之间存在摩擦时骨架固体和填充固体质点振动速度水平分量的波形情况.从图5 和图6 可以看出,两种介质之间存在摩擦时骨架固体和填充固体中都已经完全观察不到P2 和S2 波.说明两种固体之间的摩擦对P2 和S2 波的衰减作用非常明显.由于P2 和S2 波的能量主要在填充固体中传播,骨架固体和填充固体之间的摩擦使P2 与S2 波衰减,导致存在摩擦时填充固体质点振动速度水平及垂直分量波场快照与骨架固体相同.

图5 频率为20 Hz 且=2.2×108 时,0.45 s 时刻质点振动速度水平分量波场快照 (a) 骨架固体;(b) 填充固体Fig.5.Snapshots of the horizontal particle velocity component when the frequency is 20 Hz and =2.2×108 at 0.45 s:(a) Solid grain frame;(b) pore solid.

图6 频率为20 Hz 且=2.2×108 时,各相质点振动速度水平分量的波形Fig.6.Waveform comparisons of the horizontal particle velocity component between different phases when the frequency is 20 Hz and =2.2×108 (normalized).

图7 频率为20 Hz 时,=0 和=2.2×108 两种情况下各相质点振动速度水平分量的波形 (a) 骨架固体;(b) 填充固体Fig.7.Waveforms of the horizontal velocity component of each phase when =0 and =2.2×108 (the frequency is 20 Hz):(a) Solid grain frame;(b) pore solid.

孔隙介质中波的传播速度和衰减与频率有关.地震和声波测井是和水合物开发有关的两个重要手段,为了说明不同频率下两种固体之间的摩擦系数b12对波的传播特性的影响,本文还研究了声源频率为10 kHz(声波测井典型频率)的情况下b12对各种波传播速度和衰减的影响.现针对一个均匀的孔隙中充满固体的固-固孔隙介质模型,两种固体之间存在耦合.在二维x-z平面内,取x轴正方向水平向右,z轴正方向垂直向下.选用模型网格大小为800 × 800,模型x方向和z方向空间步长均为0.005 m,时间步长为2×10–7s.模型的物性参数见表1.计算时选用的雷克子波声源主频为10 kHz,声源位于模型中心位置处(对应x和z的网格点坐标分别为400 和400),接收点位于模型网格(100,780)位置处.图8 显示了频率为10 kHz时,=0 和=2.2×108两种情况下骨架固体和填充固体质点振动速度水平分量的波形.从图8可以看出,两种固体间摩擦系数的存在使P2 和S2 波均产生一定程度的衰减,但衰减比低频情况小得多(与图7 相比).此分析进一步证明了高频情况下摩擦系数对P2 和S2 波的衰减作用比低频情况小的多这一特性.

图8 频率为10 kHz 时,=0 和 =2.2×108 两种情况下各相质点振动速度水平分量的波形 (a) 骨架固体;(b) 孔隙固体Fig.8.Waveforms of the horizontal velocity component of each phase when =0 and =2.2×108 (the frequency is 10 kHz):(a) Solid grain frame;(b) pore solid.

5 结论

1)为了描述孔隙中填充一种固体的固-固孔隙介质中波的传播特性,本文由动能密度、势能密度以及耗散势函数出发,结合Lagrange 方程建立固-固孔隙介质的基本声学理论;2) 本文利用平面波分析的方法研究了固-固孔隙介质中各种波的传播特性,分析了两种纵波(P1 和P2 波)和两种横波(S1 和S2 波)的频散和衰减特性,得到P2 和S2波在低频时频散,且衰减比P1 和S1 波大得多,两种慢波的衰减随频率的增大而减小;3) 本文利用分裂法解决方程中的刚性问题,结合高阶交错网格有限差分法,提出了适用于固-固孔隙介质的时间分裂的高阶交错网格有限差分算法,模拟了声波在固-固孔隙介质中的传播.P1和S1 波的能量主要集中在骨架固体中,P2 和S2 波类似于Biot 理论中的慢纵波,是由于骨架固体和填充固体间的相对运动产生,能量主要在孔隙固体中传播;4) 本文研究了两种固体间摩擦系数和频率对各种波传播特性的影响.摩擦系数b12对P1 和S1 波的传播速度和幅度几乎无影响,主要影响P2 和S2 波的衰减.这就是在实际情况下很难观测到慢波的原因,但由于慢波也携带了一部分能量,在研究声波在介质中能量的衰减时是不能忽略的.声源的频率对激发模式的影响主要体现在幅度上,高频情况下摩擦系数对P2 和S2 波的幅度衰减较低频情况小的多.

本文的研究成果主要针对天然气水合物的特殊填充情况,显然本文方法可以推广到其他的固-固孔隙介质中的声波传播研究中.

附录A 物理量的表示及意义[26,35]

固体骨架有效质量密度:ρ11=a12(1−φ)ρs+(a21−1)φρps,其中下标s,ps 分别代表固体骨架和孔隙固体;

孔隙固体有效质量密度:ρ22=a21φρps+(a12−1)(1−φ)ρs;

两种固体耦合质量密度:ρ12=−(a12−1)(1−φ)ρs−(a21−1)φρps;

固体骨架通过孔隙固体的弯曲度:a12=1+r12φ[(1−φ)ρs+φρps]/(1−φ)ρs;

孔隙固体通过固体骨架的弯曲度:a21=1+r21(1−φ)[(1−φ)ρs+φρps]/φρps,其中a12和a21为弯曲度量级的惯性阻力参数,r12和r21表示骨架固体和孔隙固体边界的几何特征,对于球形颗粒r12=1/2[1,26,32].

固体骨架有效体积模量:K1=[(1−c1)(1−φ)]2Kav+Ksm,其中Ksm为固体骨架体积模量;

孔隙固体有效体积模量:K3=[(1−c2)φ]2Kav+Kpsm,其中Kpsm为孔隙固体骨架体积模量;

固体骨架有效剪切模量:µ1=[(1−g1)(1−φ)]2µav+µsm,其中µsm为固体骨架剪切模量;

孔隙固体有效体积模量:µ3=[(1−g2)φ]2µav+µpsm,其中µpsm为孔隙固体骨架剪切模量;

两种固体间耦合剪切模量:µ12=(1−g1)(1−g2)µav;