基于波传播系数构建填实爆炸下花岗岩中运动及变形场*

卢 强，王占江，朱玉荣，丁 洋，郭志昀

（1. 西北核技术研究院，陕西 西安 710024；2. 西北核技术研究院强动载与效应重点实验室，陕西 西安 710024）

波传播系数法的研究始于Kolsky[1]及Hunter[2]的工作，他们采用傅里叶变换法研究了线黏弹性杆中的一维波传播问题，由杆衰减系数和波数（或波速）表示线黏弹性材料的复模量，并以频率衰减因子和波数作为波传播分析的参数。自此之后，众多学者对波传播系数分析方法在杆中的应用进行了研究。Zhao 等[3-4]给出了考虑三维效应的无限长杆中纵波的通解，并进行了一些材料力学行为动态测试实验验证，结果表明利用考虑三维效应后的波传播系数去修正实验结果能够提高分析精度。Bacon 等[5-8]将前人的研究成果加以整理和完善，利用入射波和反射波的分离技术，提出了波传播系数的一点应变测试方法。Casem 等[9]及Mousavi 等[10-11]将波传播系数法在低密度泡沫材料和聚丙烯材料动态力学性能研究中进行了初步的尝试。Benatar 等[12]简化了Pochhammer-Chree 频率方程[13]，采用 ∅12 mm 和 ∅6.4 mm 的PMMA（polymethyl methacrylate）杆进行了一维波传播实验，修正了黏弹性杆中应力波传播的几何效应，并把两种杆径下实测的随频率变化的相速曲线、衰减曲线进行了对比，结果表明黏弹性杆理论可以在较宽的频率范围内确定材料的黏弹性特性。Ahonsi 等[14]采用钢球撞击作为杆中产生应力波的源，理论分析中采用了一个弹性元件和一个Maxwell 元件并联的模型（标准线性固体模型）对传播系数进行了分析。Butt 等[15-16]采用一个弹性元件和两个Maxwell 元件并联的模型（5 参数模型）分析了PMMA 杆中的传播系数，并反演了此模型对应的材料参数。Fan 等[17]采用一个非线性弹性元件和一个Maxwell 元件并联的模型分析了混凝土材料的黏弹性特性，给出了混凝土的衰减系数和波数，通过参数识别给出了混凝土材料非线性本构参数。Othman[18]采用尼龙材料作为输入输出杆对泡沫铝材料进行了SHPB 测试，施加在泡沫铝两个端面上的载荷由尼龙杆的传播系数校正（标准线性固体模型），指出若按弹性假设处理软材料的SHPB 实验数据，在小应变、高应变率时会引入不可忽略的误差。

上述研究均为杆中黏弹性波传播相关的内容，目的是为结合黏弹性霍普金森压杆实验技术来研究低阻抗材料的动态力学性能。本文中从黏弹性球面波的频率方程出发，利用花岗岩中有限个实测的球面波径向粒子速度频谱信息，给出球面波传播系数的求解方法，分析花岗岩球面波传播系数的变化，提出一种构建地下爆炸介质运动及变形场的方法。

1 黏弹性球面波的频率方程

根据黏弹性球面波的分析，地下爆炸自由场径向粒子速度和震源函数在频域内满足如下关系[19-20]：

任意两个位置 r1和 r2处 粒子速度的理论频谱比 HT(r1,r2,ω)可写为：

因此理论频谱比 HT(r1,r2,ω)可以写为：

由公式（4）可以看出，理论频谱比 HT(r1,r2,ω)是 爆心距参数 r1、 r2和 波传播系数 β(ω) 的 函数。 β (ω)是控制波在传播过程中形状改变的重要参数，客观上反映介质的黏性对波传播演化的影响。若爆心距 r1、r2处 的径向粒子速度已知，则可通过公式（4）求解出波传播系数β (ω)。

2 球面波传播系数的求解方法

利用球面波实验技术和圆环型粒子速度测试技术，可以得到有限个不同爆心距位置处的径向粒子速度[21]。利用这些测得的粒子速度中的任意两个，可以给出相应的实验频谱比 HE(r1,r2,ω)，即：

式中： ∆t 为径向粒子速度的采样时间间隔， M 和 N 分别为 r1和 r2处粒子速度的有效采样点数。

若假设岩土是黏弹性介质，则理论频谱比与实验频谱比一致，通过公式（4）定义函数g (r1,r2,β(ω))：

从公式（6）可以看出，仅有 β(ω)是需要求解的。利用Newton 迭代法，有：

式中： βn(ω)的 下标表示第 n次迭代。

数值迭代求解公式（7）的关键是确定波传播系数的初值 β0(ω)， 即要分别确定衰减因子初值 α0(ω)和波数 k0(ω)。把公式（4）展开，可得到：

式中： |HT(r1,r2,ω)|和 φT(r1,r2,ω)分 别为 HT(r1,r2,ω)的模和辐角。

把公式（9）进行简化，忽略公式右边的两项，得到波数的近似值作为其初值，即：

再把 k0(ω)代 入公式（8），即可求出衰减因子初值 α0(ω)，即：

由此给出波传播系数的初值 β0(ω)=α0(ω)+k0(ω)i，按照公式（7）经过有限次迭代即可给出收敛的波传播系数 β(ω)。

3 填实爆炸下花岗岩中传播系数的实例与讨论

本文中以王占江等[21]在0.125 g TNT 填实爆炸下实测的花岗岩中径向粒子速度为基础（如图1 所示），按照前述方法对花岗岩的波传播系数进行了分析。由于爆心距10 mm 处粒子速度计在冲击下损坏、爆心距大于60 mm 的粒子速度计受样品边界反射波的影响，这些传感器获得的粒子速度信号相对不够完整，因此在进行波传播系数分析时不予以考虑。

利用爆心距15～50 mm 处的粒子速度频谱，依次选取两个相邻测点来计算相应的实验频谱比HE(r1,r2,ω) 。 图2 给出了 r1= 15 mm、 r2= 20 mm 和 r1= 40 mm、 r2= 50 mm 时的辐角曲线 φE(r1,r2,ω)，可以看出，当 r1= 15 mm、 r2= 20 mm 时， ω≈ 1.5×107rad/s（对应频率 f≈2.38 MHz）时波数曲线开始出现突然增大或减小，这是违背物理规律的。同样，当 r1= 40 mm、 r2= 50 mm 时， φE(r1,r2,ω)出现类似的现象。本文中把这些开始出现违背物理规律的频率点视为波传播系数有效频段的上限 ωmax。因此，由实测数据计算得到的衰减因子 α(ω)和 波数 k (ω)只在有限频段内是可信的。

图3～4 分别给出了相邻测点之间区域内花岗岩传播系数中的衰减因子 α(ω)和 波数k (ω)，图5 给出了和波数曲线对应的相速度曲线。从图3～4 可以看出，测点距爆心越远，传播系数有效频段的上限ωmax越低。这是因为波传播过程中，由于介质耗散和几何发散的影响，粒子速度的高频成分衰减快，导致远区的高频信息较弱，高频信号成分的信噪比较低，从而造成波传播系数有效频段上限 ωmax的降低。

另外，由于样品尺寸小，导致信号低频成分未能充分发展即受到样品边界反射波的影响，因此衰减因子 α(ω)和 波数 k(ω)的低频结果是不可信的。按照王占江等[22]的结果，本文中使用的花岗岩用超声测得的波速为2 700 m/s。把图5 中相速度低于2 700 m/s 的部分进行标示，可以近似对衰减因子 α(ω)和波数k(ω)有 效频段的下限 ωmin进行估计，即本文中给出的在几十kHz 以下的波传播系数是不可信的。

从图5 给出的相速度曲线还可以看出，在有效频段内，相速度曲线有一个平台值，距爆心越远，这个平台值越小。按照前述的黏弹性假设，理论上获得的波传播系数应具有一致性，但从本文中处理的结果看，在15～50 mm 区域所处的应力条件下花岗岩没有体现出理想的黏弹性行为。利用卢强等[23]给出的利用球面波径向粒子速度波形反推有机玻璃力学参数的方法，图6 给出了0.125 g TNT 填实爆炸下花岗岩中等效应力峰值 τmax随爆心距 r 的变化。可以看出，在爆心距35 mm 处的等效应力峰值 τmax约为158 MPa，略大于花岗岩的单轴压缩强度154 MPa[22]，可近似认为0.125 g TNT 填实爆炸下，花岗岩弹性区的半径约为35 mm。因此，本文中所处理的区域中，爆心距15～35 mm 范围属于塑性区，35～50 mm 区域属于黏弹性区。从黏弹性区计算得到的波传播系数看，即使是低幅值的弱波，花岗岩表现出来的也不是理想的黏弹性力学行为。

图 1 0.125 g TNT 填实爆炸下花岗岩中实测的径向粒子速度[21]Fig. 1 Measured radial particle velocities in granite under the tamped explosion of 0.125 g TNT

图 2 花岗岩中实验频谱比 HE(r1,r2,ω)的辐角 φE(r1,r2,ω)随 ω的变化Fig. 2 Argument φ E(r1,r2,ω) of the experimental spectrum ratio HE(r1,r2,ω) in granite vs the circular frequency ω

图 3 利用花岗岩中相邻测点数据计算的衰减因子α(ω)Fig. 3 Attenuation factor α( ω) calculated from the data of adjacent measuring points in granite

图 4 利用花岗岩中相邻测点数据计算的波数k(ω)Fig. 4 Wave number k (ω) calculated from the data of adjacent measuring points in granite

4 利用有限个实测粒子速度构建花岗岩中运动及变形场

4.1 局部黏弹性等效方法

前面按照黏弹性假设计算了相邻测点之间花岗岩的波传播系数，衰减因子α (ω)和 波数 k (ω)基本反映出了爆炸应力波从近区的高压状态演化到相对远区的低压状态时花岗岩对波吸收和弥散的频率相关性。下面对上述不同区域的波传播系数作进一步的应用分析。

如图7 所示，爆炸应力波由 r1处传播至 r2，由两个位置处粒子速度的频谱比，可以求得一个局部黏弹性等效的波传播系数 βvisco(r1,r2,ω)。假设局部黏弹性等效成立，由公式（2）、（4）可以得到 r1和 r2之间任意位置 r处的频谱比，即：

图 5 利用花岗岩中相邻测点数据计算的相速度c(ω)Fig. 5 Phase velocity c (ω) calculated from the data of adjacent measuring points in granite

图 6 0.125 g TNT 填实爆炸下花岗岩中等效应力峰值τmax随爆心距 r的变化Fig. 6 Peak value of the equivalent stress τ max vs. r under the tamped explosion of 0.125 g TNT in granite

图 7 局部黏弹性等效下粒子速度场的构建方法Fig. 7 Method for constructing particle velocity field under local viscoelastic equivalence

若以局部理想弹性等效处理，即忽略 β(r1,r2,ωk) 中 的频率衰减因子 α (r1,r2,ωk) ，并把 β (r1,r2,ωk)中的波数k(r1,r2,ωk)以 爆炸应力波由 r1处 传播至 r2处 的平均波速 c (r1,r2)表示：

局部弹性等效条件下， r1和 r2之间任意位置 r 处的粒子速度 vr(r,t)可写为：

图8 通过花岗岩 r1= 15 mm 和 r2=25 mm 处粒子速度计算的局部黏弹性等效波传播系数βvisco(r1,r2,ω)以及局部弹性等效波传播系数βelastic(r1,r2,ω)，分别计算了r=20 mm 和25 mm 处的粒子速度。可以看出，采用局部黏弹性等效方法计算的粒子速度在r1和r2区域两端精度很高，中间位置(r=20 mm)处计算的粒子速度在峰值以及形状方面均和实验结果保持较高的相似性。以局部弹性等效方法计算的r=20，25 mm 处的粒子速度同实验结果的差异较大，无论是粒子速度峰值还是波形形状均不能很好地反映出当地粒子速度波形本来的特点。这里强调，0.125 g TNT填实爆炸下花岗岩中15、25 mm 处还是塑性区的范围，但从基于局部黏弹性等效方法计算 r1和 r2之间区域的粒子速度看，其精度远高于局部弹性等效方法，这也说明虽然局部黏弹性等效的波传播系数是基于黏弹性理论给出的结果，但在塑性区应用时仍有较好的表现。

图 8 局部黏弹性等效和局部弹性等效方法计算的粒子速度波形的比较Fig. 8 Comparison of particle velocity waveforms calculated by local viscoelastic with that by elastic equivalence method

4.2 花岗岩中变形场的分析

采用局部黏弹性等效方法，利用相邻测点获得的径向粒子速度可给出相邻测点区域内任意位置的粒子速度，即给出粒子速度的时间-空间场 vr(r,t)， 如图9 所示。对粒子速度场 vr(r,t)积分可得粒子位移场ur(r,t)， 如图10 所示。由位移场 ur(r,t)可得径向和切向应变（率）场：

图 9 采用局部黏弹性等效方法构建的粒子速度场vr(r,t)Fig. 9 Particle velocity field vr( r,t) constructed by local viscoelastic equivalence method

图 10 采用局部黏弹性等效方法构建的粒子速度场ur(r,t)Fig. 10 Particle displacement field u r(r,t) constructed by local viscoelastic equivalence method

与0.125 g TNT 填实爆炸下花岗岩中粒子速度实测位置相对应，图11～12 分别给出了花岗岩中不同位置的径向应变 εr(r,t)和 切向应变 εθ(r,t)（以压为负）。可以看出，半径15～50 mm 范围内，花岗岩中径向应变峰值由−1.7×10−2下降为−2.1×10−3，切向应变峰值由4.7×10−3下降为0.4×10−3。另外，从图11～12 还可看出，径向和切向应变达到峰值后降低一段时间，而后又发生一定的上升，这和卢强等[24]给出的理论模拟结果体现的变化规律一致。

图 11 花岗岩中的径向应变Fig. 11 Radial strain in granite at different radii

图 12 花岗岩中的切向应变Fig. 12 Tangential strain in granite at different radii

图13～14 分别给出了花岗岩中不同位置的径向应变率 ε˙r(r,t) 和 切向应变率 ε˙θ(r,t)。由图13 可以看出，径向应变率 ε˙r(r,t)在µs 级时间内由压缩加载转变为拉伸卸载。随着波传播距离的增加，径向压缩加载应变率峰值由−5.1×104s−1下降为−2.5×103s−1，径向拉伸卸载应变率峰值由3.5×104s−1下降为5.0×102s−1。由图14 可以看出，随着波传播距离的增加，切向拉伸加载的应变率峰值由5.0×103s−1下降为1.4×102s−1，切向压缩卸载的应变峰率值由−2.0×102s−1下降为−4.0×101s−1。从上述这些数据可以看出，在半径15～50 mm 区域内应变（率）峰值约有一个数量级的变化，涵盖了高应变（率）到中低应变（率）加、卸载的全过程。

图 13 花岗岩中的径向应变率Fig. 13 Radial strain rates in granite at different radii

图 14 花岗岩中的切向应变率Fig. 14 Tangential strain rates in granite at different radii

图15 给出了花岗岩中不同位置的应变状态。可以看出，爆炸近区应变状态 εr- εθ主要为压拉模式，随着波传播距离的增加，开始逐渐出现拉拉、拉压、压压模式。由静态分析结果可知，当填实爆炸激发的爆腔压力稳定时，介质的应变状态 εr- εθ为压拉模式[24-25]。这里需指出，图15 中给出的花岗岩不同位置的应变状态 εr- εθ最终会稳定在压拉模式。由于样品边界反射波的影响，图15 中远离爆心的几个位置的应变状态并不完整，其最终应变状态 εr-εθ没有处于压拉模式。

图 15 花岗岩中不同位置的应变状态Fig. 15 Strain states in granite at different radii

5 结 论

由上述分析得到以下几点结论：

（1）利用花岗岩中实测的粒子速度频谱信息，计算得到的衰减因子 α(ω)和 波数 k (ω)反映出了爆炸应力波从近区的高压状态演化到相对远区的低压状态时花岗岩对波吸收和弥散的频率相关性；

（2）花岗岩中波传播系数随爆心距的增加而变化，即使是在确定的弹性区内传播的低幅值弱波，花岗岩表现出来的也不是理想的黏弹性力学行为。换言之，花岗岩的波传播系数对其所处的应力应变状态敏感；

（3）以花岗岩中相邻测点之间区域内局部黏弹性等效假设为基础，分区域构建了填实爆炸下花岗岩介质运动和变形的时空分布，其处理精度高于局部弹性等效方法；

（4）0.125 g TNT 填实爆炸下，在半径15～50 mm 区域内：花岗岩的应变状态 εr-εθ主要为压拉模式，随着波传播距离的增加，开始逐渐出现拉拉、拉压、压压模式；花岗岩的径向应变率很快由压缩加载转变为拉伸卸载，而切向应力率则由拉伸加载转变为压缩卸载；应变（率）峰值约有一个数量级的变化，涵盖了高应变（率）到中低应变（率）加、卸载的全过程；

（5）球面波传播过程中其频率成分不断发生变化，部分频段的粒子速度信息由于粒子速度计无法响应（或响应精度降低）、测试记录设备精度不足、样品尺寸小导致信号低频成分未能充分发展、空间电磁干扰等一系列原因将影响数据的分析精度；

（6）根据本文中提出的构建地下爆炸介质运动及变形场的新方法，可进一步丰富对地下爆炸复杂应力应变状态下介质变形特征的认识。