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

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

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

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

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

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

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

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

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

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

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

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

若假設(shè)巖土是黏彈性介質(zhì)，則理論頻譜比與實驗頻譜比一致，通過公式（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)行簡化，忽略公式右邊的兩項，得到波數(shù)的近似值作為其初值，即：

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

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

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

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

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

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

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

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

圖 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 利用花崗巖中相鄰測點數(shù)據(jù)計算的衰減因子α(ω)Fig. 3 Attenuation factor α( ω) calculated from the data of adjacent measuring points in granite

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

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

4.1 局部黏彈性等效方法

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

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

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

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

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

4.2 花崗巖中變形場的分析

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

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

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

與0.125 g TNT 填實爆炸下花崗巖中粒子速度實測位置相對應(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á)到峰值后降低一段時間，而后又發(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 級時間內(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)變（率）峰值約有一個數(shù)量級的變化，涵蓋了高應(yīng)變（率）到中低應(yīng)變（率）加、卸載的全過程。

圖 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)填實爆炸激發(fā)的爆腔壓力穩(wěn)定時，介質(zhì)的應(yīng)變狀態(tài) εr- εθ為壓拉模式[24-25]。這里需指出，圖15 中給出的花崗巖不同位置的應(yīng)變狀態(tài) εr- εθ最終會穩(wěn)定在壓拉模式。由于樣品邊界反射波的影響，圖15 中遠(yuǎn)離爆心的幾個位置的應(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é) 論

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

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

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

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

（4）0.125 g TNT 填實爆炸下，在半徑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)變（率）峰值約有一個數(shù)量級的變化，涵蓋了高應(yīng)變（率）到中低應(yīng)變（率）加、卸載的全過程；

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

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