黏弹性节理岩体中应力波的传播特性分析

汪书敏，王志亮，贾帅龙，王浩然，王昊辰

(合肥工业大学 土木与水利工程学院，合肥 230009)

节理裂隙作为一种普遍存在于岩体中的缺陷结构，对应力波(爆炸波、地震波等)在岩体中的传播规律产生较大影响。应力波与节理之间的作用关系复杂，常常导致振幅下降、波形变化和能量损失[1-2]。岩体中的天然裂隙常常由饱和砂土、黏土以及风化岩石等具有黏弹性变形特性的介质填充，在自然条件下逐渐演化成黏弹性节理，故研究具有黏弹特性的节理对应力波传播的影响机制在抗震工程、损伤测定和岩石力学等领域都具有重大意义。

国内外学者针对应力波在节理岩体中的传播过程进行了大量研究。王志亮等[3]提出了一种三参数模型用于描述节理岩体的动态响应，基于位移不连续理论(DDM)推导了一维应力波的传播方程并与实验进行对比，发现当应力波持续时间较短时，理论结果与实验结果吻合较好。Li等[4-7]根据波前动量守恒原理提出了一种新的时域递归方法，相继导出了应力波通过线弹性节理、非线弹性节理以及考虑节理切向滑移特性时的传播方程。而对于天然夹层裂隙中常常充填有饱和砂土、黏土等形成的节理，Fehler等[8-9]认为应考虑其黏弹特性，对此，也有不少学者进行了相关研究。Huang等[10]利用矩阵传播法研究了应力波通过一组黏弹性节理的传播规律，发现在Kelvin节理中存在一个特定黏度系数值使得能量耗散率最大。Yi等[11]采用等效介质法将节理充填介质等效为静态的各向同性介质，发现该方法不能准确反映填充介质较疏时的振幅各向异性和波速随频率变化的规律。Zhu等[12]考虑填充介质质量产生的惯性应力，提出了位移与应力不连续理论(DSDM)，发现应力波通过填充节理时的能量耗散率由节理黏度系数和填充介质质量决定。Wu等[13]使用分离式霍普金森岩杆对充填节理进行的动态试验结果表明，DSDM比DDM更适合描述厚度在10 mm以上的填充节理对波传播的影响。此外，Li等[14-16]提出了一种薄层界面模型用于研究应力波入射充填节理时的传播规律，结果表明，在剪切波入射该模型引起剪应力大于薄层界面的抗剪强度而产生滑移时该方法仍有效。

综上，目前基于时域递归方法(TDRM)且考虑应力不连续条件的平行黏弹性节理中波传播特性研究开展较少。因此，本研究先采用TDRM法，引入饱依丁-汤姆逊本构模型描述节理的黏弹特性，并考虑填充节理质量，将DDM扩展至DSDM，导出了应力波在一组平行黏弹性节理中传播的完整解析方程。通过与频域内的封闭解对比验证该方法的正确性。最后,分析了节理模型参数、无量纲节理厚度和平行节理间距等参数对透反射系数和能量耗散率的影响，力求得出有参考价值的结论。

1 理论分析

1.1 问题描述

在离爆炸、地震等扰动源较远的地层中，应力波相对于节理面可看作平面波。当应力波穿过一组等间距平行节理时，在节理处会发生透射和反射现象，分别形成透射P波与S波，反射P波与S波，如图1所示，S为节理间距。图2为应力波通过第X条节理的左右波形图。其中，d为节理厚度，符号“-”和“+”分别代表节理的左侧与右侧，LP-、LS-、RP-和RS-分别为节理左侧的左行P波、左行S波、右行P波和右行S波，LP+、LS+、RP+和RS+分别为节理右侧的左行P波、左行S波、右行P波和右行S波。假设节理两侧的岩石均为线弹性介质，且密度与泊松比相同，P波入射角为α，S波入射角为β。应力波通过含黏弹特性的充填介质时，会产生一定的边际弥散效应而损失部分到达节理右侧的能量。本研究不考虑应力波入射充填节理时产生的散射效应，认为该透反射过程遵循Snell定律。故右侧出射波角度与左侧入射角相同，即透反射P波和S波与法线夹角均分别为α和β[12,15,17]。

图1 应力波通过平行节理Fig.1 Stress wave traversing through parallel joint

图2 第X条节理两侧的左行波和右行波Fig.2 Left-and right-running waves on both sides of the Xth joint

1.2 理论公式推导

(1)

(2)

根据Snell定律

(3)

式中：cp和cs分别为P波和S波在岩石中的波速。图3(b)～3(h)分析过程同图3(a)。

图3 应力波波前与节理两侧应力Fig.3 Stress on wavefront and rock joint surfaces

结合上述结果及波前动量守恒原理，根据Li等[4]推导过程可得节理两侧的法向应力、切向应力、法向速度与切向速度：

(4)

(5)

(6)

(7)

式中：上标e分别代表“-”或“+”。下标r和l分别代表右行波和左行波。

本研究引入饱依丁-汤姆逊本构模型对岩体节理的黏弹特性进行描述，如图4所示。由一个Maxwell体与弹簧元件并联而成，本构方程为

图4 饱依丁-汤姆逊本构模型[18]Fig.4 Poyting-Thomson constitutive model[18]

(8)

节理左右界面的连接如图2所示，考虑节理填充介质质量的应力不连续条件[12,19]：

(9)

(10)

将式(10)作为连接节理左右两侧法向和切向的位移不连续条件可以得

(11)

将式(4)～(7)代入式(9)和(11)，并整理为矩阵形式：

(12)

(13)

式中的相关系数矩阵如下：

假定节理间的岩石为完全弹性，则应力波由第X-1条节理右侧出射到达第X条节理左侧或由第X+1条节理左侧反射至第X条节理右侧的过程中，只相差一个时间差，故时移函数方程表达式为

(14)

(15)

式中：下标X表示第X条节理,S表示相邻节理间的间距。在计算过程中，时间步长为Δt，则两侧速度时移步数差为np=S/(cosα·cp·Δt)或ns=S/(cosβ·cs·Δt)。

式(12)～(15)即为应力波在一组平行节理中传播的完整方程。通过编制MATLAB程序，输入相应已知条件和初始条件并采用满足精度要求的时间步长，对上述方程进行时域递归计算即可得解。

2 理论验证

为方便研究节理对应力波在岩体中传播特征的影响，定义反射和透射系数分别为

(16)

(17)

式中：下标a和b分别表示P波或S波，下标N表示平行节理条数。

引入兼具Maxwell体与Kelvin体黏弹特性的饱依丁-汤姆逊黏弹本构模型作为位移不连续的变形条件。当不考虑节理的黏性特征和应力不连续条件时，可将黏性系数η及惯性质量m均设为0，便退化为线弹性节理，容易发现退化后的方程与Li等[4]一致。当弹簧元件k1趋于正无穷时，原模型演变为Kelvin模型，适合描述饱和土填充节理的动态响应规律[21-22]。而将弹簧元件k2设为0时，原模型便退化为适合描述横波穿过黏土填充的薄层节理情形的Maxwell模型[23]。为验证上述所推导方程的正确性，将计算结果与Zhu等[12]所得的结果进行对比。假定在完整岩石中P波和S波的传播速度分别为6 131和3 830 m/s，节理填充介质中cplate=1 000 m/s，这些参数均可通过岩石和节理填充介质的杨氏模量、密度和泊松比求得[12,19]。根据Snell定律可得入射P波和入射S波的入射临界角分别为90°与38.65°，分析中仅讨论入射角小于临界角的情况。定义两个法向特定刚度系数分别为Kn1=kn1/(zpω)与Kn2=kn2/(zpω)，特定黏度系数Hn=ηn/zp，两个切向特定刚度系数分别为Kτ1=kτ1/(zsω)与Kτ2=kτ2/(zsω)，特定黏度系数为Hτ=ητ/zs。ω为入射波的圆频率，即ω=2πf，zp和zs分别为P波和S波的波阻抗。不同入射角下的透反射系数计算结果如图5、6所示，可以看出，Kelvin节理和Maxwell节理随入射角的变化趋势大致相同，临界角之前透射系数变化较为平稳，临界角附近急剧变化，这与Zhu等[12]所得频域内封闭解具有良好的一致性。

3 参数讨论与分析

节理模型参数往往通过霍普金森压杆(SHPB)对充填介质进行动态冲击实验的方式获取[12]。但由于SHPB实验中S波难以获得，切向参数的确定存在困难。现假设节理法向与切向同性，即取Kτ1=Kn1，Kτ2=Kn2，Hτ=Hn，除待研究参数外，其余参数同上。设入射波为半个周期的正弦P波或S波，如式(18)所示，Is=Ip。

(18)

式中f为正弦波的频率，取1 kHz。

定义法向入射波的能量耗散率为[12]

eloss=1-(R2+T2)

(19)

3.1 节理模型参数的影响分析

为探究模型参数对透反射系数与能量耗散率的影响，以P波通过单条节理为例，特定刚度系数Kn2分别取为0.2、0.4、0.8和1.6。透反射系数及能量耗散率随1/Kn1和1/Hn的变化曲线分别如图7、8所示。可以看出，当其他参数为定值时，Kn2对节理模型的影响较为明显，透射系数随着Kn2的增加而增大，而反射系数与能量耗散率均减小。

由图7可知，当Kn2较小时，随着1/Kn1的增加，反射系数逐渐增加，而透射系数初期有一小段上升，随后缓慢下降；当Kn2较大时，反射和透射系数的初期上升趋势趋弱，且后续随1/Kn1的增大基本持平。能量耗散率随着1/Kn1的增大而逐渐下降，并趋于一定值e0。由于Kn2的增加，模型中Maxwell体的影响逐渐弱化，节理的变形特征由黏弹性逐渐向线弹性演变，因而能量耗散率相应下降。故e0主要由Kn2决定，且随Kn2的增大而减小。

由图8可知，当Kn2较小时，反射系数初期小幅下降，且下降的趋势随Kn2的增加而逐渐趋缓，但总体随1/Hn增加而上升。同时，透射系数随1/Hn的增加而逐渐下降，且下降幅度随Kn2增加而逐渐减小。能量耗散率eloss随1/Hn的增加逐渐上升，而当1/Hn达2～3之后则基本持平。

3.2 无量纲节理厚度的影响分析

应力波在节理岩体传播的过程中，节理厚度同样会对应力波的传播特性产生影响。考虑到透反射系数的波频依赖性，定义无量纲节理厚度h为ωd/(2πcplate)。图9为IS和IP法向入射单条节理时透反射系数及能量耗散率随h的变化曲线。可以看出，当S波法向入射时，透反射系数及能量耗散率的大小随lgh的变化趋势与P波法向入射时基本相同。当h<10-2时，无量纲节理厚度对透反射系数的影响不明显。随着h的增大，透射系数出现明显的下降，并逐渐接近于0。而反射系数随着无量纲节理厚度的增加呈现先降低后增加的现象。此外，随着厚度增加，S波入射时透反射系数的变化趋势要先于P波，而能量耗散率随无量纲节理厚度的增加表现出先上升后下降的趋势。

图9 透反射系数与能量耗散率随无量纲节理厚度的变化Fig.9 Variation of transmission and reflection and energy dissipation rate with non-dimensional joint thickness

3.3 应力波入射角的影响分析

取入射S波表达式同P波，图10、11分别为不同入射角下P波和S波在节理间距为1 m的一组平行等间距节理中传播的透射系数TPP、TPS、TSS和TSP。可以看出，当P波和S波以任意角度入射时，随着节理条数的增加，TPP和TSS均出现了不同程度的下降。在入射P波与S波分别接近临界角αc=90°与βc=arcsin(cs/cp)=38.6°时，TPP和TSS分别出现急剧下降或上升的现象。随着节理条数N的增加，透反射系数在临界角处的“尖点”显著下降，而TPS和TSP随着节理条数的增加而逐渐增大。

图10 一组平行节理对不同入射角度P波透射系数的影响Fig.10 Influence of a set of parallel joints on the transmission coefficients of incident P-wave at different angles

图11 一组平行节理对不同入射角度S波透射系数的影响Fig.11 Influence of a set of parallel joints on the transmission coefficients of incident S-wave at different angles

3.4 无量纲节理间距的影响分析

定义无量纲节理间距ξ为节理间距S与入射应力波波长的比值。本节采用半个周期的正弦入射P波，如式(18)所示。图12为透射系数与无量纲节理间距ξ的关系，易见透射系数随ξ的变化可以分为上升区间(ξ≤ξcri)、下降区间(ξcri<ξ<ξthr)以及不变区间(ξthr≤ξ)。在节理条数较多时，TPP随无量纲节理间距的增加呈线性下降，而黏弹节理则是分为两段——陡峭段和平缓段下降。图13为临界值ξcri和ξthr随节理数N的变化曲线，可以看出,随着节理数的增加，临界值ξcri缓慢下降，且基本保持线性关系；对于临界值ξthr，当N<5时，增长较缓；当N>6时，增长较快。

图12 透射系数与无量纲节理间距ξ的关系Fig.12 Transmission coefficients versus non-dimensional joint spacing (ξ)

图13 ξcri和ξthr随节理数N的变化Fig.13 Variation of ξcri and ξthr with joint number N

4 结 论

1)时域递归方法可用于推导任意角度的入射应力波通过黏弹性节理的传播方程。引入的饱依丁-汤姆逊模型退化为Kelvin模型和Maxwell模型的计算结果不仅验证了所用模型及推导过程的正确性，同时印证了其兼具两种模型的特点，为描述节理复杂的黏弹特性提供了新的模型选择。

2)参数化研究表明，刚度系数、黏度系数、无量纲节理厚度等都影响透反射系数及能量耗散率。其中，特定刚度系数Kn2与透射系数呈正相关，与反射系数及能量耗散率呈负相关。当无量纲节理厚度小于10-2时，其变化对透反射系数及能量耗散率的影响不明显。

3)透反射系数随黏弹性节理数及节理间距的变化而变化。对于入射P波(S波)，透射P波(S波)随节理数增加而逐渐下降，透射S波(P波)随节理数的增加逐渐上升。当节理条数较少时，节理条数的变化对透射系数影响较大；节理条数较多时，其变化对透射系数的影响则较小。此外，ξcri随节理数的增长基本保持线性关系，但ξthr随节理数的增加呈现先缓后陡的增长趋势。