花崗巖中實測球面波粒子速度的時域和頻域分析

盧 強，王占江，朱玉榮，丁 洋，郭志昀，張永偉

(西北核技術研究所，西安710024； 強動載與效應重點實驗室，西安710024)

球面波實驗是研究波傳播規(guī)律和3維應力-應變條件下介質(zhì)動力學特性的重要手段之一[1-4]，被國內(nèi)外眾多學者所關注。美國進行了大量的球面波實驗研究，涉及的介質(zhì)有凝灰?guī)r、花崗巖、頁巖、沖積土、混凝土等[5-9]。Larson對球面波加載下巖石特性與爆炸能量之間的耦合關系進行了大量的研究[5-7]。Nagy采用0.375 g的太安球形裝藥對低孔隙率脆性巖石中的球面波傳播特征進行了研究，分析了圍壓、孔隙率、含水量等因素對波傳播演化的影響[8]，指出了實測球面波徑向粒子速度波形與采用傳統(tǒng)的理想彈性假設方法預測的粒子速度波形，在峰值和脈寬上均存在較大差異，但研究中未考慮波傳播衰減過程中對頻率的依賴性。James等報道了φ1 829 mm×1 829 mm混凝土樣品中454 g TNT炸藥球爆炸加載實驗[9]，目的是通過測量球形爆炸下快速膨脹空腔周圍的應力和粒子運動參數(shù)，驗證與評估混凝土中侵徹和爆炸效應的計算模型。1985年，中國科學院武漢巖土力學研究所開展了球面波化爆實驗[10]，建立了以∏形電磁粒子速度計和壓電晶體應力計為測量手段、以微秒電雷管起爆直徑30 mm鑄裝炸藥球作波源的室內(nèi)巖石球面波化爆實驗裝置。王占江等建立了0.1 g TNT的球面波實驗裝置，并獲得多種介質(zhì)中球面波徑向粒子速度數(shù)據(jù)[11]。盧強等利用實測球面波數(shù)據(jù)，結(jié)合黏彈性理論對球面波在黃土及有機玻璃介質(zhì)中的傳播特性進行了分析[12-18]。

本文以0.125 g TNT、直徑為5 mm的炸藥球填實爆炸加載下花崗巖中實測球面波數(shù)據(jù)為基礎，給出了球面波在花崗巖中傳播演化的時域特征，結(jié)合黏彈性球面波傳播理論對球面波傳播的頻率特征進行了近似分析，并對分析結(jié)果進行了評估。

1 花崗巖中球面波實驗簡介

花崗巖是一種典型的介質(zhì)材料，屬于火成巖?；◢弾r主要成分為鉀長石(40%～45%)、酸性斜長石(30%～35%)、石英(18%～22%)、黑云母(2%～5%)、磁鐵礦(0.1%～1%)等，礦物顆粒直徑為0.2～5 mm。用于球面波實驗的花崗巖樣品是從米級尺度的新鮮巖石上切割成型后再由專門的巖石加工設備精加工而成[1]。本文所用的花崗巖取自華山[19]，花崗巖樣品的密度為2.60 g·cm-3，尺寸為φ276 mm×300 mm，爆炸源為0.125 g TNT精密微型炸藥球。考慮到樣品邊界反射波對粒子速度計的影響，離爆心最遠的粒子速度計設置在半徑為120 mm處。本文主要分析填實爆炸條件下，花崗巖半徑為10～120 mm時的球面波數(shù)據(jù)。球面波實驗示意圖，如圖1所示。

圖1球面波實驗示意圖Fig.1Schematic of spherical wave experiment

2 花崗巖中球面波傳播的時域特征

圖2--圖5給出了測點半徑為10～25 mm，30～50 mm，60～90 mm，100～120 mm時實測的球面波徑向粒子速度曲線。

圖2半徑為10～25 mm時，花崗巖中實測的粒子速度曲線Fig.2Measured particle velocities in granite (r=10～25 mm)

圖3半徑為30～50 mm時，花崗巖中實測的粒子速度曲線Fig.3Measured particle velocities in granite(r=30～50 mm)

圖4半徑為60～90 mm時，花崗巖中實測的粒子速度曲線Fig.4Measured particle velocities in granite(r=60～90 mm)

圖5半徑為100～120 mm時，花崗巖中實測的粒子速度曲線Fig.5Measured particle velocities in granite(r=100～120 mm)

圖6和圖7分別給出了花崗巖中球面波粒子速度峰值vmax和比例粒子位移峰值umax隨比距離R(R=r/Q1/3)的變化規(guī)律。同時，對實驗數(shù)據(jù)進行了單段線性擬合和雙段線性擬合，擬合結(jié)果一并繪于圖中。單段線性擬合公式為

(1)

(2)

雙段線性擬合公式為

(3)

(4)

從圖6和圖7可以看出，粒子速度峰值和比例粒子位移峰值均隨波傳播距離的增加而減??；單段線性擬合結(jié)果和實驗結(jié)果偏差較大，而雙段線性擬合結(jié)果與實驗結(jié)果符合較好。從雙段線性擬合結(jié)果可以看出：粒子速度峰值和比例粒子位移峰值在不同區(qū)域內(nèi)的衰減規(guī)律是不同的，離爆心較近的區(qū)域衰減較慢，離爆心較遠的區(qū)域粒子衰減較快。可見，花崗巖中球形波的衰減不僅包含了傳播距離的幾何衰減，還包含其非線性動力學特性引入的黏性衰減。

圖6花崗巖中粒子速度峰值隨比距離的變化Fig.6Peak values of particle velocity vs. scaled distance in granite

圖7花崗巖中比例粒子位移峰值隨比距離的變化Fig.7Peak values of scaled particle displacement vs. scaled distance in granite

圖8為0.125 g TNT炸藥球填實爆炸下，花崗巖中實測的球面波徑向粒子速度波形的時間特征隨比距離的變化關系。可以得到：在測點半徑10 mm處，波形上升沿ΔTr約為0.5 μs；在測點半徑120 mm處，ΔTr約為5.3 μs；二者相比，ΔTr展寬了約10倍。另外，ΔTr、波形半高寬ΔT1/2及波形正向脈寬ΔT+均有隨著比距離增加而增加的趨勢。球面波在花崗巖中傳播的展寬效應也反映了其非線性動力學響應的復雜性，基于理想彈性假設的波傳播理論，無法對其波形的傳播演化規(guī)律進行合理解釋。

圖8花崗巖中粒子速度波形的時間特征隨比距離的變化Fig.8Time features of particle velocity waveforms vs. scaled distance in granite

3 花崗巖中球面波粒子速度的頻率響應函數(shù)

3.1 粒子速度的頻率響應函數(shù)典型特征

按照粒子速度的頻率響應函數(shù)定義，可以給出以測點半徑r1、r2處球面波徑向粒子速度為輸入量的頻率響應函數(shù)的計算公式為[13, 18, 20]

(5)

在對粒子速度的頻率響應函數(shù)進行分析時，如果單獨分析幅頻特性會損失相頻特性的信息，只分析相頻特性同樣也會損失幅頻特性的信息。鑒于粒子速度響應函數(shù)Hvr(r2,r1，ωi)的實部和虛部均包含了幅頻特性及相頻特性信息，因此，將其分別作為分析對象。實部和虛部可以分別寫為

Re[Hvr(r2,r1，ωi)]=

|[Hvr(r2,r1，ωi)| cos[φvr(r2,r1，ω)]

(6)

Im[Hvr(r2,r1，ωi)]=

|[Hvr(r2,r1，ωi)| sin[φvr(r2,r1，ω)]

(7)

式中，|[Hvr(r2,r1，ωi)|描述了幅頻特性；φvr(r2,r1，ω)描述了相頻特性。φvr(r2,r1，ω)的計算公式為

φvr(r2,r1，ω) =unwrap{angle[Hvr(r2,r1，ωi)]}

(8)

式中，angle為輻角函數(shù)；unwrap為相位展開函數(shù)，其功能是檢查數(shù)據(jù)相位跳變并糾正其跳變。

按照式(5)的定義，圖9和圖10分別給出了r1為35 mm，r2為50 mm時，花崗巖中粒子速度的頻率響應函數(shù)的實部和虛部隨圓頻率的變化曲線，并采用幅值按指數(shù)衰減的正弦振蕩函數(shù)曲線進行了擬合，擬合結(jié)果一并在圖中給出。擬合函數(shù)為

Re[Hvr(r2,r1，ωi)]=Ae-ω/ω0sin[π(ω-ωc)/ω1]

(9-1)

Im[Hvr(r2,r1，ωi)]=Ae-ω/ω0sin[π(ω-ωc)/ω1]

(9-2)

式中，A代表擬合函數(shù)的幅值；ω0為延遲因子，rad·s-1；ω1和振蕩周期相關；ωc和波形的初始相位相關。

從圖9和圖10可以看出，實部和虛部的振蕩幅值均隨著圓頻率的增加而逐漸降低；擬合的結(jié)果在ω>1.0×106rad·s-1的高頻衰減部分與實驗曲線符合程度較好，在ω<1.0×106rad·s-1的低頻部分與實驗曲線符合較差。

圖9花崗巖中粒子速度的頻率響應函數(shù)的實部隨圓頻率ω的變化曲線(r1=35 mm，r2=50 mm)Fig.9The real part of the particle velocity response function in granite(r1=35 mm，r2=50 mm)

圖10花崗巖中粒子速度的頻率響應函數(shù)的虛部隨圓頻率ω的變化曲線(r1=35 mm，r2=50 mm)Fig.10The imaginary part of the particle velocity response function in granite(r1=35 mm，r2=50 mm)

3.2 粒子速度的頻率響應函數(shù)有效頻段的估計方法

理論上，任何位置的球面波粒子速度信號均應包含由源激發(fā)的所有頻率信息，受到介質(zhì)黏性以及幾何發(fā)散的影響，粒子速度的部分頻段受到抑制的程度高，而有些頻段受到抑制的程度低，但用任何兩個位置的粒子速度信息均應能給出粒子速度所有頻率對應的響應值。對于球面波實驗，球面波傳播過程中其頻率成分不斷發(fā)生變化，部分頻段的粒子速度信息由于粒子速度計無法響應或響應精度降低、測試記錄設備精度不足、樣品尺寸小導致信號低頻成分未能充分展現(xiàn)等一系列原因，將影響后續(xù)的數(shù)據(jù)分析精度，這就要求對利用實測粒子速度分析粒子速度的頻率響應函數(shù)的有效頻段進行評估。

按式(8)給出的當r1為15 mm，r2為20 mm和r1為40 mm，r2為50 mm時的相頻特性曲線，如圖11所示?？梢钥闯?，當r1為15 mm，r2為20 mm時，相頻特性曲線在ω>1.15×107rad·s-1(對應頻率f>1.83 MHz)時開始出現(xiàn)異常變化，這說明用該位置實測粒子速度計算的頻率響應函數(shù)的有效頻段的上限ωmax為1.15×107rad·s-1(fmax為1.83 MHz)。同理，用r1=40 mm，r2=50 mm位置實測粒子速度計算的ωmax為5.64×106rad·s-1(fmax為0.90 MHz)。對比圖11中近區(qū)和遠區(qū)的相頻特性曲線，可以發(fā)現(xiàn)，遠區(qū)粒子速度的頻率響應函數(shù)的有效頻段上限降低，這是因為隨著波的傳播，遠區(qū)的粒子速度由于介質(zhì)耗散、幾何發(fā)散等原因?qū)е赂哳l粒子速度信息喪失，高頻成分信號弱、空間電磁噪聲、記錄設備精度等均會對數(shù)據(jù)獲取精度造成影響。

圖11花崗巖中粒子速度的典型相頻特性曲線Fig.11Typical phase-frequency characteristic of particle velocity in granite

粒子速度的頻率響應函數(shù)的有效頻段下限ωmin或fmin受多種因素的影響，包括樣品尺寸小導致的低頻信號發(fā)展不充分、空間電磁噪聲對低幅度粒子速度信號的掩蓋等，對其評估較為困難。本文利用粒子速度信號的采樣頻率及實際采樣點數(shù)，近似給出有效頻段下限值的確定方法。

本文花崗巖球面波實驗中的信號采樣頻率S為50 M·s-1，扣除邊界反射波影響后測試的有效持續(xù)時間為50～60 μs，L約為2 500。按照傅里葉理論分析，頻率的分辨率近似為

(10)

由式(10)可推斷，粒子速度的頻率響應函數(shù)有效頻段的下限fmin應高于信號的頻率分辨率Δf。由此可得，在本文的測試區(qū)域內(nèi)，粒子速度的頻率響應函數(shù)的有效頻段應滿足：

10 kHz=Δf≤fmin

(11)

這里需指出的是，fmax和采用哪兩個位置的粒子速度信號直接相關，fmin是較為粗略的估計，總體上看式(11)只是對有效頻段的近似估計，便于工程上參考使用。

3.3 粒子速度的頻率響應函數(shù)的高頻近似分析

考慮黏彈性條件下，式(5)的幅頻及相頻特性的理論公式可寫為[13, 20]：

|Hvr(r2,r1,ωi)|=

(12)

φvr(r2,r1,ω)=-k(ω)(r2-r1)+

(13)

式(12)和(13)中，α(ω)和k(ω)分別為頻率衰減因子和波數(shù)，它們共同構(gòu)成波傳播系數(shù)β(ω)。β(ω)可寫為[13]

(14)

式中，C(ω)=ω/k(ω)為相速度。波傳播系數(shù)β(ω)控制波傳播過程中的衰減及形狀變化。

式(12)和式(13)是非線性方程，為滿足工程需要，可做近似分析。當k(ω)滿足：

且k(ω)?α(ω)

(15)

則可將式(12)和式(13)近似為

(16)

φvr(r2,r1,ω)=-k(ω)(r2-r1)

(17)

以r1=15 mm，r2=20 mm時得到的粒子速度的頻率響應函數(shù)為例，按上述分析可得其有效頻段的估計值為

10 kHz=Δf≤fmin

(18)

按時域給出的波傳播速度估計值c≈5 000 m·s-1計算，波數(shù)的最大值kmax(ω)為

(19)

由式(19)可以看出，在有效頻段接近fmax的頻率區(qū)間，k(ω)滿足式(15)。另外，由文獻[4]給出的方法，利用式(16)及兩個不同位置的粒子速度峰值，估計出高頻段的α(ω)約為21.7，也滿足式(15)。此高頻條件下，粒子速度的頻率響應函數(shù)的實部和虛部可以近似寫為

cos[-k(ω)(r2-r1)]

(20)

sin[-k(ω)(r2-r1)]

(21)

由式(20)和式(21)可知，粒子速度的頻率響應函數(shù)實部和虛部的振蕩特性主要由k(ω)的變化特征和傳播距離r2-r1決定，其衰減特征主要由α(ω)的變化特征和r2-r1決定。在k(ω)和α(ω)確定的條件下，r2-r1控制著粒子速度的頻率響應函數(shù)實部和虛部的衰減和振蕩特性。把式(9)同式(20)及式(21)進行對比，可以給出：

(22)

(23)

從式(22)和式(23)可以看出：頻率衰減因子α(ω)和波數(shù)k(ω)均與圓頻率ω近似成線性關系，其比例系數(shù)分別為κ和λ，則

β(ω)=α(ω)+k(ω)i=κω+λωi

(24)

在測點半徑為15～70 mm內(nèi)布置粒子速度計測量粒子速度信息，并計算κ，λ，fmax，結(jié)果如表1所列。

表1基于實測粒子速度計算的κ，λ，fmax

Tab.1κ，λ，fmaxcalculatedwithmeasuredparticlevelocities

r1/mmr2/mmκ/(10-5 s·m-1)λ/(10-4 s·m-1)fmax/MHz15201.002.191.8320251.051.981.7725300.682.211.6230352.032.201.6135403.362.471.3340501.792.670.9050603.522.750.6960704.102.910.49

從表1可以看出，由不同位置的實測粒子速度給出的fmax隨波傳播距離的增加而逐漸降低，這反映了波傳播過程中高頻成分的耗散。κ和λ在不同頻段內(nèi)的擬合值不同，這表明α(ω)和k(ω)在不同頻段內(nèi)變化的快慢不同，反映了球形應力波在花崗巖中傳播的復雜性。

4 討論

由式(3)可知，在測點半徑為10～120 mm的測試區(qū)域內(nèi)，近似以r為28 mm(R為56 m·kt-1/3)為分界點，靠近爆心的區(qū)域內(nèi)粒子速度峰值的冪衰減指數(shù)為1.46，遠離爆心的區(qū)域內(nèi)粒子速度峰值的冪衰減指數(shù)為2.60。由表1可以看出，由r1和r2均小于30 mm區(qū)域內(nèi)的實測粒子速度給出的頻率衰減因子的比例系數(shù)κ整體上明顯小于r1和r2均大于30 mm區(qū)域內(nèi)的κ值。時域和頻域的分析結(jié)果都表明，在半徑r為30 mm(R為60 m·kt-1/3)附近時，球面波粒子速度峰值的衰減出現(xiàn)了轉(zhuǎn)折。另外，上述結(jié)果也表明：本文提出的球面波頻域分析方法給出的結(jié)果和時域的結(jié)果在物理上是相符的。

為對本文提出的球面波傳播的頻域分析方法進行評估，利用表1給出的波傳播系數(shù)及r1處粒子速度作為輸入量，計算r2處的粒子速度，并將其和實測粒子速度進行比較。按照粒子速度的頻率響應函數(shù)的定義，r2處粒子速度的頻域計算公式可以寫為[20]

(25)

vr(r2，t)=

(26)

式(26)中，N為粒子速度頻域?qū)狞c數(shù)。

將式(26)中的波傳播系數(shù)β(ω)分別取為

β(ω)=κω+λωi

(27)

β(ω)=λωi

(28)

式(27)表示利用r1處粒子速度并結(jié)合式(26)反演粒子速度時，同時考慮了波傳播過程中衰減因子和波數(shù)的影響，而式(28)則沒有考慮衰減因子的影響。

圖12預測的粒子速度波形和實測波形的對比(r1=15 mm，r2=20 mm)Fig.12Comparison between the predicted particle velocity waveform and the measured waveform (r1=15 mm，r2=20 mm)

圖12和圖13分別給出了r1=15 mm，r2=20 mm和r1=20 mm，r2=25 mm時的實測粒子速度曲線同預測的粒子速度曲線的對比?？梢钥闯?，采用式(27)反演的粒子速度峰值比采用式(28)反演的粒子速度峰值更接近實測結(jié)果，這反映了介質(zhì)黏性對波衰減的影響。還可以看出，無論是基于式(27)還是基于式(28)的假設，反演的粒子速度波形和實測粒子速度波形之間均有較大的差別。

圖13預測的粒子速度波形和實測波形的對比(r1=20 mm，r2=25 mm)Fig.13Comparison between the predicted particle velocity waveform and the measured waveform (r1=20 mm，r2=25 mm)

圖14給出了花崗巖中粒子速度峰值預測結(jié)果和實測結(jié)果的對比。可以看出，本文給出的球面波傳播系數(shù)的線性近似計算方法考慮了波傳播過程中衰減因子的影響，與傳統(tǒng)的局部理想彈性假設方法[8]相比，提高了對粒子速度峰值的預測精度，但對粒子速度波形演化過程中形狀的變化，其預測能力不高。解決這個問題的關鍵是利用式(5)直接給出粒子速度的頻率響應函數(shù)的實驗結(jié)果，進而求解式(12)和式(13)構(gòu)成的非線性方程組，確切獲得波傳播系數(shù)β(ω)隨頻率的變化，這部分內(nèi)容將在后續(xù)研究中予以關注。

圖14花崗巖中粒子速度峰值預測結(jié)果和實測結(jié)果的對比Fig.14Comparison of peak values of particle velocity obtained by predicted and measured in granite

5 結(jié)論

本文分析了0.125 g TNT、直徑為5 mm的炸藥球填實爆炸加載下花崗巖中球面波傳播的時域特征，提出了一種基于黏彈性假設的球面波粒子速度頻域分析方法，得到如下結(jié)論：

1)對花崗巖中粒子速度峰值和比例粒子位移峰值進行對數(shù)坐標下的線性擬合分析，結(jié)果表明粒子速度峰值和比例粒子位移峰值在靠近爆心的區(qū)域衰減慢，遠離爆心的區(qū)域衰減快；

2)從花崗巖中實測球面波粒子速度波形的時域和頻域分析結(jié)果可知，在測點半徑約為30 mm(比距離約為60 m·kt-1/3)處，粒子速度峰值的衰減規(guī)律出現(xiàn)由慢變快的轉(zhuǎn)折；

3)在高頻近似分析下，花崗巖中頻率衰減因子α(ω)和波數(shù)k(ω)均與圓頻率ω近似成線性關系，但由不同相鄰位置實測的粒子速度信息得到的α(ω)和k(ω)各不相同；

4)球面波傳播系數(shù)的線性近似處理方法對波傳播演化過程中波形形狀的預測精度不高，但相比于傳統(tǒng)的局部理想彈性假設方法，本文的近似處理考慮了波傳播過程中衰減因子的影響，可以提高粒子速度峰值的預測精度。