波传播系数法的研究始于Kolsky^[1]及Hunter^[2]的工作，他们采用傅里叶变换法研究了线黏弹性杆中的一维波传播问题，由杆衰减系数和波数（或波速）表示线黏弹性材料的复模量，并以频率衰减因子和波数作为波传播分析的参数。自此之后，众多学者对波传播系数分析方法在杆中的应用进行了研究。Zhao等^[3-4]给出了考虑三维效应的无限长杆中纵波的通解，并进行了一些材料力学行为动态测试实验验证，结果表明利用考虑三维效应后的波传播系数去修正实验结果能够提高分析精度。Bacon等^[5-8]将前人的研究成果加以整理和完善，利用入射波和反射波的分离技术，提出了波传播系数的一点应变测试方法。Casem等^[9]及Mousavi等^[10-11]将波传播系数法在低密度泡沫材料和聚丙烯材料动态力学性能研究中进行了初步的尝试。Benatar等^[12]简化了Pochhammer-Chree频率方程^[13]，采用$\varnothing$12 mm和$\varnothing$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]：

式中：$\omega = 2{{{\text{π}} }}f$为圆频率，$f$为频率，$r$为爆心距，${\overline v _{\rm{r}}}(r,\omega )$和$\overline \gamma (r,\omega )$分别为径向粒子速度${v_{\rm{r}}}(r,t)$和震源函数$\gamma (r,t)$的频谱，$\beta (\omega ) = \alpha (\omega ) + k(\omega ){\rm{i}}$为球面波传播系数，$\alpha (\omega )$为频率衰减因子，$k(\omega ) = \omega /c(\omega )$为波数，$c(\omega )$为相速度，${\rm{i = }}\sqrt { - 1}$为虚数单位。

任意两个位置${r_1}$和${r_2}$处粒子速度的理论频谱比${H_{\rm{T}}}({r_1},{r_2},\omega )$可写为：

由于$\overline \gamma ({r_1},\omega )$和$\overline \gamma ({r_2},\omega )$之间满足：

因此理论频谱比${H_{\rm{T}}}({r_1},{r_2},\omega )$可以写为：

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

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

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

式中：$\Delta t$为径向粒子速度的采样时间间隔，$M$和$N$分别为${r_1}$和${r_2}$处粒子速度的有效采样点数。

若假设岩土是黏弹性介质，则理论频谱比与实验频谱比一致，通过公式（4）定义函数$g({r_1},{r_2},\beta (\omega ))$：

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

式中：${\beta _{{n}}}(\omega )$的下标表示第$n$次迭代。

数值迭代求解公式（7）的关键是确定波传播系数的初值${\beta _0}(\omega )$，即要分别确定衰减因子初值${\alpha _0}(\omega )$和波数${k_0}(\omega )$。把公式（4）展开，可得到：

式中：$\left| {{H_{\rm{T}}}({r_1},{r_2},\omega )} \right|$和${\varphi _{\rm{T}}}({r_1},{r_2},\omega )$分别为${H_{\rm{T}}}({r_1},{r_2},\omega )$的模和辐角。

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

再把${k_0}(\omega )$代入公式（8），即可求出衰减因子初值${\alpha _0}(\omega )$，即：

由此给出波传播系数的初值${\beta _0}(\omega ) = {\alpha _0}(\omega ) + {k_0}(\omega )i$，按照公式（7）经过有限次迭代即可给出收敛的波传播系数$\beta (\omega )$。

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

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

图  1  0.125 g TNT填实爆炸下花岗岩中实测的径向粒子速度^[21]

Figure  1.  Measured radial particle velocities in granite under the tamped explosion of 0.125 g TNT

下载: 全尺寸图片 幻灯片

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

图  2  花岗岩中实验频谱比${{H_{\rm{E}}}({r_1},{r_2},\omega )}$的辐角${{\varphi _{\rm{E}}}({r_1},{r_2},\omega )}$随${\omega }$的变化

Figure  2.  Argument ${{\varphi _{\rm{E}}}({r_1},{r_2},\omega )}$ of the experimental spectrum ratio ${{H_{\rm{E}}}({r_1},{r_2},\omega )}$ in granite vs the circular frequency ${\omega }$

下载: 全尺寸图片 幻灯片

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

图  3  利用花岗岩中相邻测点数据计算的衰减因子$\alpha (\omega )$

Figure  3.  Attenuation factor $\alpha (\omega )$ calculated from the data of adjacent measuring points in granite

下载: 全尺寸图片 幻灯片

图  4  利用花岗岩中相邻测点数据计算的波数${k(\omega )}$

Figure  4.  Wave number ${k(\omega )}$ calculated from the data of adjacent measuring points in granite

下载: 全尺寸图片 幻灯片

图  5  利用花岗岩中相邻测点数据计算的相速度$c(\omega )$

Figure  5.  Phase velocity $c(\omega )$ calculated from the data of adjacent measuring points in granite

下载: 全尺寸图片 幻灯片

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

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

图  6  0.125 g TNT填实爆炸下花岗岩中等效应力峰值${\tau _{\max }}$随爆心距$r$的变化

Figure  6.  Peak value of the equivalent stress ${\tau _{\max }}$ vs. r under the tamped explosion of 0.125 g TNT in granite

下载: 全尺寸图片 幻灯片

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

4.1   局部黏弹性等效方法

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

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

图  7  局部黏弹性等效下粒子速度场的构建方法

Figure  7.  Method for constructing particle velocity field under local viscoelastic equivalence

下载: 全尺寸图片 幻灯片

对上式进行整理，并对${\overline v _{\rm{r}}}(r,\omega )$进行傅里叶逆变换，可得$r$处的粒子速度${v_{\rm{r}}}(r,t)$，即：

式中：$L$为${\overline v _{\rm{r}}}({r_1},{\omega _{{k}}})$的有效长度。

若以局部理想弹性等效处理，即忽略$\beta ({r_1},{r_2},{\omega _{{k}}})$中的频率衰减因子$\alpha ({r_1},{r_2},{\omega _{{k}}})$，并把$\beta ({r_1},{r_2},{\omega _{{k}}})$中的波数$k({r_1},{r_2},{\omega _{{k}}})$以爆炸应力波由${r_1}$处传播至${r_2}$处的平均波速$c({r_1},{r_2})$表示：

局部弹性等效条件下，${r_1}$和${r_2}$之间任意位置$r$处的粒子速度${v_{\rm r}}(r,t)$可写为：

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

图  8  局部黏弹性等效和局部弹性等效方法计算的粒子速度波形的比较

Figure  8.  Comparison of particle velocity waveforms calculated by local viscoelastic with that by elastic equivalence method

下载: 全尺寸图片 幻灯片

4.2   花岗岩中变形场的分析

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

图  9  采用局部黏弹性等效方法构建的粒子速度场${v_{\rm{r}}}(r,t)$

Figure  9.  Particle velocity field ${v_{\rm{r}}}(r,t)$ constructed by local viscoelastic equivalence method

下载: 全尺寸图片 幻灯片

图  10  采用局部黏弹性等效方法构建的粒子速度场${u_{\rm{r}}}(r,t)$

Figure  10.  Particle displacement field ${u_{\rm{r}}}(r,t)$ constructed by local viscoelastic equivalence method

下载: 全尺寸图片 幻灯片

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

图  11  花岗岩中的径向应变

Figure  11.  Radial strain in granite at different radii

下载: 全尺寸图片 幻灯片

图  12  花岗岩中的切向应变

Figure  12.  Tangential strain in granite at different radii

下载: 全尺寸图片 幻灯片

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

图  13  花岗岩中的径向应变率

Figure  13.  Radial strain rates in granite at different radii

下载: 全尺寸图片 幻灯片

图  14  花岗岩中的切向应变率

Figure  14.  Tangential strain rates in granite at different radii

下载: 全尺寸图片 幻灯片

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

图  15  花岗岩中不同位置的应变状态

Figure  15.  Strain states in granite at different radii

下载: 全尺寸图片 幻灯片

5.   结　论

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

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

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

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

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

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

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