盧 強(qiáng)，王占江，朱玉榮，丁 洋，郭志昀

（1. 西北核技術(shù)研究院，陜西 西安 710024；2. 西北核技術(shù)研究院強(qiáng)動(dòng)載與效應(yīng)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710024）

波傳播系數(shù)法的研究始于Kolsky[1]及Hunter[2]的工作，他們采用傅里葉變換法研究了線黏彈性桿中的一維波傳播問(wèn)題，由桿衰減系數(shù)和波數(shù)（或波速）表示線黏彈性材料的復(fù)模量，并以頻率衰減因子和波數(shù)作為波傳播分析的參數(shù)。自此之后，眾多學(xué)者對(duì)波傳播系數(shù)分析方法在桿中的應(yīng)用進(jìn)行了研究。Zhao 等[3-4]給出了考慮三維效應(yīng)的無(wú)限長(zhǎng)桿中縱波的通解，并進(jìn)行了一些材料力學(xué)行為動(dòng)態(tài)測(cè)試實(shí)驗(yàn)驗(yàn)證，結(jié)果表明利用考慮三維效應(yīng)后的波傳播系數(shù)去修正實(shí)驗(yàn)結(jié)果能夠提高分析精度。Bacon 等[5-8]將前人的研究成果加以整理和完善，利用入射波和反射波的分離技術(shù)，提出了波傳播系數(shù)的一點(diǎn)應(yīng)變測(cè)試方法。Casem 等[9]及Mousavi 等[10-11]將波傳播系數(shù)法在低密度泡沫材料和聚丙烯材料動(dòng)態(tài)力學(xué)性能研究中進(jìn)行了初步的嘗試。Benatar 等[12]簡(jiǎn)化了Pochhammer-Chree 頻率方程[13]，采用 ?12 mm 和 ?6.4 mm 的PMMA（polymethyl methacrylate）桿進(jìn)行了一維波傳播實(shí)驗(yàn)，修正了黏彈性桿中應(yīng)力波傳播的幾何效應(yīng)，并把兩種桿徑下實(shí)測(cè)的隨頻率變化的相速曲線、衰減曲線進(jìn)行了對(duì)比，結(jié)果表明黏彈性桿理論可以在較寬的頻率范圍內(nèi)確定材料的黏彈性特性。Ahonsi 等[14]采用鋼球撞擊作為桿中產(chǎn)生應(yīng)力波的源，理論分析中采用了一個(gè)彈性元件和一個(gè)Maxwell 元件并聯(lián)的模型（標(biāo)準(zhǔn)線性固體模型）對(duì)傳播系數(shù)進(jìn)行了分析。Butt 等[15-16]采用一個(gè)彈性元件和兩個(gè)Maxwell 元件并聯(lián)的模型（5 參數(shù)模型）分析了PMMA 桿中的傳播系數(shù)，并反演了此模型對(duì)應(yīng)的材料參數(shù)。Fan 等[17]采用一個(gè)非線性彈性元件和一個(gè)Maxwell 元件并聯(lián)的模型分析了混凝土材料的黏彈性特性，給出了混凝土的衰減系數(shù)和波數(shù)，通過(guò)參數(shù)識(shí)別給出了混凝土材料非線性本構(gòu)參數(shù)。Othman[18]采用尼龍材料作為輸入輸出桿對(duì)泡沫鋁材料進(jìn)行了SHPB 測(cè)試，施加在泡沫鋁兩個(gè)端面上的載荷由尼龍桿的傳播系數(shù)校正（標(biāo)準(zhǔn)線性固體模型），指出若按彈性假設(shè)處理軟材料的SHPB 實(shí)驗(yàn)數(shù)據(jù)，在小應(yīng)變、高應(yīng)變率時(shí)會(huì)引入不可忽略的誤差。

上述研究均為桿中黏彈性波傳播相關(guān)的內(nèi)容，目的是為結(jié)合黏彈性霍普金森壓桿實(shí)驗(yàn)技術(shù)來(lái)研究低阻抗材料的動(dòng)態(tài)力學(xué)性能。本文中從黏彈性球面波的頻率方程出發(fā)，利用花崗巖中有限個(gè)實(shí)測(cè)的球面波徑向粒子速度頻譜信息，給出球面波傳播系數(shù)的求解方法，分析花崗巖球面波傳播系數(shù)的變化，提出一種構(gòu)建地下爆炸介質(zhì)運(yùn)動(dòng)及變形場(chǎng)的方法。

1 黏彈性球面波的頻率方程

根據(jù)黏彈性球面波的分析，地下爆炸自由場(chǎng)徑向粒子速度和震源函數(shù)在頻域內(nèi)滿足如下關(guān)系[19-20]：

任意兩個(gè)位置 r1和 r2處 粒子速度的理論頻譜比 HT(r1,r2,ω)可寫為：

因此理論頻譜比 HT(r1,r2,ω)可以寫為：

由公式（4）可以看出，理論頻譜比 HT(r1,r2,ω)是 爆心距參數(shù) r1、 r2和 波傳播系數(shù) β(ω) 的 函數(shù)。 β (ω)是控制波在傳播過(guò)程中形狀改變的重要參數(shù)，客觀上反映介質(zhì)的黏性對(duì)波傳播演化的影響。若爆心距 r1、r2處 的徑向粒子速度已知，則可通過(guò)公式（4）求解出波傳播系數(shù)β (ω)。

2 球面波傳播系數(shù)的求解方法

利用球面波實(shí)驗(yàn)技術(shù)和圓環(huán)型粒子速度測(cè)試技術(shù)，可以得到有限個(gè)不同爆心距位置處的徑向粒子速度[21]。利用這些測(cè)得的粒子速度中的任意兩個(gè)，可以給出相應(yīng)的實(shí)驗(yàn)頻譜比 HE(r1,r2,ω)，即：

式中： ?t 為徑向粒子速度的采樣時(shí)間間隔， M 和 N 分別為 r1和 r2處粒子速度的有效采樣點(diǎn)數(shù)。

若假設(shè)巖土是黏彈性介質(zhì)，則理論頻譜比與實(shí)驗(yàn)頻譜比一致，通過(guò)公式（4）定義函數(shù)g (r1,r2,β(ω))：

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

式中： βn(ω)的 下標(biāo)表示第 n次迭代。

數(shù)值迭代求解公式（7）的關(guān)鍵是確定波傳播系數(shù)的初值 β0(ω)， 即要分別確定衰減因子初值 α0(ω)和波數(shù) k0(ω)。把公式（4）展開，可得到：

式中： |HT(r1,r2,ω)|和 φT(r1,r2,ω)分 別為 HT(r1,r2,ω)的模和輻角。

把公式（9）進(jìn)行簡(jiǎn)化，忽略公式右邊的兩項(xiàng)，得到波數(shù)的近似值作為其初值，即：

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

由此給出波傳播系數(shù)的初值 β0(ω)=α0(ω)+k0(ω)i，按照公式（7）經(jīng)過(guò)有限次迭代即可給出收斂的波傳播系數(shù) β(ω)。

3 填實(shí)爆炸下花崗巖中傳播系數(shù)的實(shí)例與討論

本文中以王占江等[21]在0.125 g TNT 填實(shí)爆炸下實(shí)測(cè)的花崗巖中徑向粒子速度為基礎(chǔ)（如圖1 所示），按照前述方法對(duì)花崗巖的波傳播系數(shù)進(jìn)行了分析。由于爆心距10 mm 處粒子速度計(jì)在沖擊下?lián)p壞、爆心距大于60 mm 的粒子速度計(jì)受樣品邊界反射波的影響，這些傳感器獲得的粒子速度信號(hào)相對(duì)不夠完整，因此在進(jìn)行波傳播系數(shù)分析時(shí)不予以考慮。

利用爆心距15～50 mm 處的粒子速度頻譜，依次選取兩個(gè)相鄰測(cè)點(diǎn)來(lái)計(jì)算相應(yīng)的實(shí)驗(yàn)頻譜比HE(r1,r2,ω) 。 圖2 給出了 r1= 15 mm、 r2= 20 mm 和 r1= 40 mm、 r2= 50 mm 時(shí)的輻角曲線 φE(r1,r2,ω)，可以看出，當(dāng) r1= 15 mm、 r2= 20 mm 時(shí)， ω≈ 1.5×107rad/s（對(duì)應(yīng)頻率 f≈2.38 MHz）時(shí)波數(shù)曲線開始出現(xiàn)突然增大或減小，這是違背物理規(guī)律的。同樣，當(dāng) r1= 40 mm、 r2= 50 mm 時(shí)， φE(r1,r2,ω)出現(xiàn)類似的現(xiàn)象。本文中把這些開始出現(xiàn)違背物理規(guī)律的頻率點(diǎn)視為波傳播系數(shù)有效頻段的上限 ωmax。因此，由實(shí)測(cè)數(shù)據(jù)計(jì)算得到的衰減因子 α(ω)和 波數(shù) k (ω)只在有限頻段內(nèi)是可信的。

圖3～4 分別給出了相鄰測(cè)點(diǎn)之間區(qū)域內(nèi)花崗巖傳播系數(shù)中的衰減因子 α(ω)和 波數(shù)k (ω)，圖5 給出了和波數(shù)曲線對(duì)應(yīng)的相速度曲線。從圖3～4 可以看出，測(cè)點(diǎn)距爆心越遠(yuǎn)，傳播系數(shù)有效頻段的上限ωmax越低。這是因?yàn)椴▊鞑ミ^(guò)程中，由于介質(zhì)耗散和幾何發(fā)散的影響，粒子速度的高頻成分衰減快，導(dǎo)致遠(yuǎn)區(qū)的高頻信息較弱，高頻信號(hào)成分的信噪比較低，從而造成波傳播系數(shù)有效頻段上限 ωmax的降低。

另外，由于樣品尺寸小，導(dǎo)致信號(hào)低頻成分未能充分發(fā)展即受到樣品邊界反射波的影響，因此衰減因子 α(ω)和 波數(shù) k(ω)的低頻結(jié)果是不可信的。按照王占江等[22]的結(jié)果，本文中使用的花崗巖用超聲測(cè)得的波速為2 700 m/s。把圖5 中相速度低于2 700 m/s 的部分進(jìn)行標(biāo)示，可以近似對(duì)衰減因子 α(ω)和波數(shù)k(ω)有 效頻段的下限 ωmin進(jìn)行估計(jì)，即本文中給出的在幾十kHz 以下的波傳播系數(shù)是不可信的。

從圖5 給出的相速度曲線還可以看出，在有效頻段內(nèi)，相速度曲線有一個(gè)平臺(tái)值，距爆心越遠(yuǎn)，這個(gè)平臺(tái)值越小。按照前述的黏彈性假設(shè)，理論上獲得的波傳播系數(shù)應(yīng)具有一致性，但從本文中處理的結(jié)果看，在15～50 mm 區(qū)域所處的應(yīng)力條件下花崗巖沒有體現(xiàn)出理想的黏彈性行為。利用盧強(qiáng)等[23]給出的利用球面波徑向粒子速度波形反推有機(jī)玻璃力學(xué)參數(shù)的方法，圖6 給出了0.125 g TNT 填實(shí)爆炸下花崗巖中等效應(yīng)力峰值 τmax隨爆心距 r 的變化?？梢钥闯觯诒木?5 mm 處的等效應(yīng)力峰值 τmax約為158 MPa，略大于花崗巖的單軸壓縮強(qiáng)度154 MPa[22]，可近似認(rèn)為0.125 g TNT 填實(shí)爆炸下，花崗巖彈性區(qū)的半徑約為35 mm。因此，本文中所處理的區(qū)域中，爆心距15～35 mm 范圍屬于塑性區(qū)，35～50 mm 區(qū)域?qū)儆陴椥詤^(qū)。從黏彈性區(qū)計(jì)算得到的波傳播系數(shù)看，即使是低幅值的弱波，花崗巖表現(xiàn)出來(lái)的也不是理想的黏彈性力學(xué)行為。

圖 1 0.125 g TNT 填實(shí)爆炸下花崗巖中實(shí)測(cè)的徑向粒子速度[21]Fig. 1 Measured radial particle velocities in granite under the tamped explosion of 0.125 g TNT

圖 2 花崗巖中實(shí)驗(yàn)頻譜比 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 利用花崗巖中相鄰測(cè)點(diǎn)數(shù)據(jù)計(jì)算的衰減因子α(ω)Fig. 3 Attenuation factor α( ω) calculated from the data of adjacent measuring points in granite

圖 4 利用花崗巖中相鄰測(cè)點(diǎn)數(shù)據(jù)計(jì)算的波數(shù)k(ω)Fig. 4 Wave number k (ω) calculated from the data of adjacent measuring points in granite

4 利用有限個(gè)實(shí)測(cè)粒子速度構(gòu)建花崗巖中運(yùn)動(dòng)及變形場(chǎng)

4.1 局部黏彈性等效方法

前面按照黏彈性假設(shè)計(jì)算了相鄰測(cè)點(diǎn)之間花崗巖的波傳播系數(shù)，衰減因子α (ω)和 波數(shù) k (ω)基本反映出了爆炸應(yīng)力波從近區(qū)的高壓狀態(tài)演化到相對(duì)遠(yuǎn)區(qū)的低壓狀態(tài)時(shí)花崗巖對(duì)波吸收和彌散的頻率相關(guān)性。下面對(duì)上述不同區(qū)域的波傳播系數(shù)作進(jìn)一步的應(yīng)用分析。

如圖7 所示，爆炸應(yīng)力波由 r1處傳播至 r2，由兩個(gè)位置處粒子速度的頻譜比，可以求得一個(gè)局部黏彈性等效的波傳播系數(shù) βvisco(r1,r2,ω)。假設(shè)局部黏彈性等效成立，由公式（2）、（4）可以得到 r1和 r2之間任意位置 r處的頻譜比，即：

圖 5 利用花崗巖中相鄰測(cè)點(diǎn)數(shù)據(jù)計(jì)算的相速度c(ω)Fig. 5 Phase velocity c (ω) calculated from the data of adjacent measuring points in granite

圖 6 0.125 g TNT 填實(shí)爆炸下花崗巖中等效應(yīng)力峰值τ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 局部黏彈性等效下粒子速度場(chǎng)的構(gòu)建方法Fig. 7 Method for constructing particle velocity field under local viscoelastic equivalence

若以局部理想彈性等效處理，即忽略 β(r1,r2,ωk) 中 的頻率衰減因子 α (r1,r2,ωk) ，并把 β (r1,r2,ωk)中的波數(shù)k(r1,r2,ωk)以 爆炸應(yīng)力波由 r1處 傳播至 r2處 的平均波速 c (r1,r2)表示：

局部彈性等效條件下， r1和 r2之間任意位置 r 處的粒子速度 vr(r,t)可寫為：

圖8 通過(guò)花崗巖 r1= 15 mm 和 r2=25 mm 處粒子速度計(jì)算的局部黏彈性等效波傳播系數(shù)βvisco(r1,r2,ω)以及局部彈性等效波傳播系數(shù)βelastic(r1,r2,ω)，分別計(jì)算了r=20 mm 和25 mm 處的粒子速度?？梢钥闯?，采用局部黏彈性等效方法計(jì)算的粒子速度在r1和r2區(qū)域兩端精度很高，中間位置(r=20 mm)處計(jì)算的粒子速度在峰值以及形狀方面均和實(shí)驗(yàn)結(jié)果保持較高的相似性。以局部彈性等效方法計(jì)算的r=20，25 mm 處的粒子速度同實(shí)驗(yàn)結(jié)果的差異較大，無(wú)論是粒子速度峰值還是波形形狀均不能很好地反映出當(dāng)?shù)亓Ｗ铀俣炔ㄐ伪緛?lái)的特點(diǎn)。這里強(qiáng)調(diào)，0.125 g TNT填實(shí)爆炸下花崗巖中15、25 mm 處還是塑性區(qū)的范圍，但從基于局部黏彈性等效方法計(jì)算 r1和 r2之間區(qū)域的粒子速度看，其精度遠(yuǎn)高于局部彈性等效方法，這也說(shuō)明雖然局部黏彈性等效的波傳播系數(shù)是基于黏彈性理論給出的結(jié)果，但在塑性區(qū)應(yīng)用時(shí)仍有較好的表現(xiàn)。

圖 8 局部黏彈性等效和局部彈性等效方法計(jì)算的粒子速度波形的比較Fig. 8 Comparison of particle velocity waveforms calculated by local viscoelastic with that by elastic equivalence method

4.2 花崗巖中變形場(chǎng)的分析

采用局部黏彈性等效方法，利用相鄰測(cè)點(diǎn)獲得的徑向粒子速度可給出相鄰測(cè)點(diǎn)區(qū)域內(nèi)任意位置的粒子速度，即給出粒子速度的時(shí)間-空間場(chǎng) vr(r,t)， 如圖9 所示。對(duì)粒子速度場(chǎng) vr(r,t)積分可得粒子位移場(chǎng)ur(r,t)， 如圖10 所示。由位移場(chǎng) ur(r,t)可得徑向和切向應(yīng)變（率）場(chǎng)：

圖 9 采用局部黏彈性等效方法構(gòu)建的粒子速度場(chǎng)vr(r,t)Fig. 9 Particle velocity field vr( r,t) constructed by local viscoelastic equivalence method

圖 10 采用局部黏彈性等效方法構(gòu)建的粒子速度場(chǎng)ur(r,t)Fig. 10 Particle displacement field u r(r,t) constructed by local viscoelastic equivalence method

與0.125 g TNT 填實(shí)爆炸下花崗巖中粒子速度實(shí)測(cè)位置相對(duì)應(yīng)，圖11～12 分別給出了花崗巖中不同位置的徑向應(yīng)變 εr(r,t)和 切向應(yīng)變 εθ(r,t)（以壓為負(fù)）?？梢钥闯?，半徑15～50 mm 范圍內(nèi)，花崗巖中徑向應(yīng)變峰值由?1.7×10?2下降為?2.1×10?3，切向應(yīng)變峰值由4.7×10?3下降為0.4×10?3。另外，從圖11～12 還可看出，徑向和切向應(yīng)變達(dá)到峰值后降低一段時(shí)間，而后又發(fā)生一定的上升，這和盧強(qiáng)等[24]給出的理論模擬結(jié)果體現(xiàn)的變化規(guī)律一致。

圖 11 花崗巖中的徑向應(yīng)變Fig. 11 Radial strain in granite at different radii

圖 12 花崗巖中的切向應(yīng)變Fig. 12 Tangential strain in granite at different radii

圖13～14 分別給出了花崗巖中不同位置的徑向應(yīng)變率 ε˙r(r,t) 和 切向應(yīng)變率 ε˙θ(r,t)。由圖13 可以看出，徑向應(yīng)變率 ε˙r(r,t)在μs 級(jí)時(shí)間內(nèi)由壓縮加載轉(zhuǎn)變?yōu)槔煨遁d。隨著波傳播距離的增加，徑向壓縮加載應(yīng)變率峰值由?5.1×104s?1下降為?2.5×103s?1，徑向拉伸卸載應(yīng)變率峰值由3.5×104s?1下降為5.0×102s?1。由圖14 可以看出，隨著波傳播距離的增加，切向拉伸加載的應(yīng)變率峰值由5.0×103s?1下降為1.4×102s?1，切向壓縮卸載的應(yīng)變峰率值由?2.0×102s?1下降為?4.0×101s?1。從上述這些數(shù)據(jù)可以看出，在半徑15～50 mm 區(qū)域內(nèi)應(yīng)變（率）峰值約有一個(gè)數(shù)量級(jí)的變化，涵蓋了高應(yīng)變（率）到中低應(yīng)變（率）加、卸載的全過(guò)程。

圖 13 花崗巖中的徑向應(yīng)變率Fig. 13 Radial strain rates in granite at different radii

圖 14 花崗巖中的切向應(yīng)變率Fig. 14 Tangential strain rates in granite at different radii

圖15 給出了花崗巖中不同位置的應(yīng)變狀態(tài)。可以看出，爆炸近區(qū)應(yīng)變狀態(tài) εr- εθ主要為壓拉模式，隨著波傳播距離的增加，開始逐漸出現(xiàn)拉拉、拉壓、壓壓模式。由靜態(tài)分析結(jié)果可知，當(dāng)填實(shí)爆炸激發(fā)的爆腔壓力穩(wěn)定時(shí)，介質(zhì)的應(yīng)變狀態(tài) εr- εθ為壓拉模式[24-25]。這里需指出，圖15 中給出的花崗巖不同位置的應(yīng)變狀態(tài) εr- εθ最終會(huì)穩(wěn)定在壓拉模式。由于樣品邊界反射波的影響，圖15 中遠(yuǎn)離爆心的幾個(gè)位置的應(yīng)變狀態(tài)并不完整，其最終應(yīng)變狀態(tài) εr-εθ沒有處于壓拉模式。

圖 15 花崗巖中不同位置的應(yīng)變狀態(tài)Fig. 15 Strain states in granite at different radii

5 結(jié) 論

由上述分析得到以下幾點(diǎn)結(jié)論：

（1）利用花崗巖中實(shí)測(cè)的粒子速度頻譜信息，計(jì)算得到的衰減因子 α(ω)和 波數(shù) k (ω)反映出了爆炸應(yīng)力波從近區(qū)的高壓狀態(tài)演化到相對(duì)遠(yuǎn)區(qū)的低壓狀態(tài)時(shí)花崗巖對(duì)波吸收和彌散的頻率相關(guān)性；

（2）花崗巖中波傳播系數(shù)隨爆心距的增加而變化，即使是在確定的彈性區(qū)內(nèi)傳播的低幅值弱波，花崗巖表現(xiàn)出來(lái)的也不是理想的黏彈性力學(xué)行為。換言之，花崗巖的波傳播系數(shù)對(duì)其所處的應(yīng)力應(yīng)變狀態(tài)敏感；

（3）以花崗巖中相鄰測(cè)點(diǎn)之間區(qū)域內(nèi)局部黏彈性等效假設(shè)為基礎(chǔ)，分區(qū)域構(gòu)建了填實(shí)爆炸下花崗巖介質(zhì)運(yùn)動(dòng)和變形的時(shí)空分布，其處理精度高于局部彈性等效方法；

（4）0.125 g TNT 填實(shí)爆炸下，在半徑15～50 mm 區(qū)域內(nèi)：花崗巖的應(yīng)變狀態(tài) εr-εθ主要為壓拉模式，隨著波傳播距離的增加，開始逐漸出現(xiàn)拉拉、拉壓、壓壓模式；花崗巖的徑向應(yīng)變率很快由壓縮加載轉(zhuǎn)變?yōu)槔煨遁d，而切向應(yīng)力率則由拉伸加載轉(zhuǎn)變?yōu)閴嚎s卸載；應(yīng)變（率）峰值約有一個(gè)數(shù)量級(jí)的變化，涵蓋了高應(yīng)變（率）到中低應(yīng)變（率）加、卸載的全過(guò)程；

（5）球面波傳播過(guò)程中其頻率成分不斷發(fā)生變化，部分頻段的粒子速度信息由于粒子速度計(jì)無(wú)法響應(yīng)（或響應(yīng)精度降低）、測(cè)試記錄設(shè)備精度不足、樣品尺寸小導(dǎo)致信號(hào)低頻成分未能充分發(fā)展、空間電磁干擾等一系列原因?qū)⒂绊憯?shù)據(jù)的分析精度；

（6）根據(jù)本文中提出的構(gòu)建地下爆炸介質(zhì)運(yùn)動(dòng)及變形場(chǎng)的新方法，可進(jìn)一步豐富對(duì)地下爆炸復(fù)雜應(yīng)力應(yīng)變狀態(tài)下介質(zhì)變形特征的認(rèn)識(shí)。