黏彈性固體中地下爆炸輻射地震波能量的演化*

盧 強，丁 洋，劉赟哲，唐仕英，郭志昀，王占江

（西北核技術研究所，陜西 西安 710024）

研究地下封閉爆炸時通常把爆炸源簡化為球形爆源，這樣地下爆炸問題可簡化為球面波問題。地下爆炸能量與源區介質發生強烈的耦合，最終以波的形式向周圍介質輻射地震波能量。為表征地下爆炸或地震激發的震源強度，學術界提出了多種模型[1-3]。

根據Mueller等[4-5]提出的二階震源模型，地下爆炸震源可用作用在彈性半徑上的穩態壓力疊加一個指數衰減的動態壓力組成。在地下爆炸地震效應研究中，一般用折合位移勢 ψ (τ) 或折合速度勢 γ (τ)的振幅譜表征震源強度[6]，折合位移勢穩態值 ψ∞和爆后介質彈性半徑處的穩態壓力（也稱準靜態壓力）呈線性關系。同時，折合位移勢穩態值 ψ∞還 可以和地震學中反映源區不可恢復非彈性形變的地震矩M0聯系起來。結合Murphy[7]給出M0∝ψ∞的 線性經驗公式，則爆后彈性半徑處的穩態壓力和地震矩M0也呈線性關系。在研究波傳播路徑上的震動效應時，地下爆炸空腔中的準靜態壓力和動態壓力對震動效應都有貢獻，僅用和準靜態壓力相關的地震矩或折合位移勢穩態值表征地下爆炸震源強度是不完備的，不能全面反映爆炸空腔內動態壓力在地介質內引起的震動效應。

Choy等[2]和Boatwright等[3]利用大量的地震事件，基于能量的概念建立了地震波能量和能量震級之間的關系，以此表征震源強度。袁乃榮等[8]和李贊等[9]利用國家地震臺網和全球地震臺網的寬頻地震波數據對地震波的能量震級進行研究，指出地震波對傳播路徑區域造成的破壞直接與地震釋放的地震波能量相關，可見地下爆炸向介質中耦合的地震波能量是表征地下爆炸震源強度的另外一個重要參數。

為研究地下爆炸與介質的能量耦合效應，Sanchidriár等[10]利用多次野外巖石爆破，結合地震儀記錄圖、高速攝像機拍攝的巖石爆破面的運動速度以及爆破后巖石碎塊尺寸的分布等實驗結果，從地震波能量、動能、斷裂能等方面研究了爆炸能量的組成。田振農等[11]以塊體離散元計算方法為基礎，研究了巖體中爆腔內的壓力脈動及爆炸能量分布特征，指出爆源近區巖體結構性質對爆腔內壓力脈動規律影響不大，爆腔內壓力峰值與炸藥量無關，給出了巖石斷裂破壞消耗能量、變形能、用于拋擲的能量的比例。郭家豪等[12]基于反映高程放大的質點振動公式建立了侵徹鉆地彈爆炸地震震源能量的計算模型，并結合經驗公式推算爆炸總能量，從而評估鉆地彈威力。吳亮等[13]對耦合裝藥條件下柱狀炸藥爆炸的能量分布進行了研究，分析了爆炸能量和巖石介質的耦合機制，計算給出了爆炸能量與介質耦合過程中不同力學過程的能量占比。肖衛國等[14]基于地下爆炸試驗實測的地表速度，以折合速度勢低頻穩態值（即折合位移穩態值）表示的地震耦合強度，研究了黃土、砂礫石和花崗巖場地地下爆炸的地震耦合效應，指出對于相同當量、相當比例埋深下的地下爆炸，地震耦合強度和源區介質特性有著強烈的依賴關系。

本文中，利用廣義Maxwell黏彈性體中球面波Laplace域的解，反演計算黃土中地下爆炸的速度場、位移場、應力場和應變場，以任意兩個不同半徑處球面之間區域作為觀測區域，建立地下爆炸自由場輻射地震波能量的計算方法，分析地震波能量的傳播演化規律，相關結果對于理解地下爆炸地震波能量在黏彈性介質中的衰減規律以及分析地下爆炸中實測的自由場數據具有一定參考作用。

1 地下爆炸黏彈性區輻射能量的組成

球面波條件下，某個半徑r的球面上，其徑向應力為 σr(r,t)，則徑向面力Fr(r,t) 可以表示為 4 πr2σr(r,t)。球面波下材料僅有徑向位移ur，則面力在徑向位移上的做功W可以寫為[10]：

式中：vr為徑向粒子速度，r為爆心距。

力Fr(r,t)在 半徑r的球面上做的功W(r,t) 即 是流入此球面的總能量。根據能量守恒，能量流經r0≤r≤r1的區域時，流入能量W(r0,t) 轉化為此區域內所有質點的動能Ek(r0,r1,t)（即輻射動能）、介質應變能（也稱內能或勢能）Ep(r0,r1,t) 、材料耗散能Ed(r0,r1,t) 以及從半徑r1的 球面上流出的能量W(r1,t)，如圖1所示。則有：

圖1 黏彈性固體中爆炸地震波能量的組成Fig.1 The composition of the energy of explosion seismic wave in visco-elastic solids

式中：Ek(r0,r1,t)和Ep(r0,r1,t) 分別為t時刻整個區域r0≤r≤r1內所有質點動能和應變能密度函數的空間積分，則有：

式中： ρ 為材料密度， εr為徑向應變， εθ為切向應變， σθ為切向應力。

式（2）是基于能量守恒獲得的，適用于任意固體中球面波能量傳播演化的分析。若要求解球面波能量流經r0≤r≤r1的 區域時材料的耗散能Ed(r0,r1,t)，則有：

考慮t→∞ 的特殊情況，此時波完全離開r0≤r≤r1區域，則此區域內質點動能Ek(r0,r1,∞)=0 。按照黏彈性球面波的解，式（4）可以寫為[15-16]：

式中： σ0s為穩態空腔壓力，Ge為黏彈性材料的靜態剪切模量。從式（6）可以看出，穩態空腔壓力 σ0s形成的勢能主要分布在幾倍空腔半徑的范圍內，比如當r1=10r0時，此觀測區域儲存了99.9%的勢能。

由式（3）～（6）可給出球面波能量流經r0≤r≤r1的 區域時材料的耗散能Ed(r0,r1,∞)，有：

2 黏彈性固體中地下爆炸輻射地震波場的計算方法

根據無限介質中黏彈性球面波理論，在Laplace域，以波傳播系數表征的黏彈性球面波粒子速度(r,s)、粒子位移(r,s)、徑向應力(r,s)、切向應力(r,s)、徑向應變(r,s)、切向應變(r,s)可寫為[15-16]：

式中： μ 為材料的泊松比，s為Laplace變量， β(s)和(r0,s) 分別為波傳播系數和彈性半徑為r0的球腔上壓力邊界條件的Laplace表示。

由式（8）～（13）可知，若 ρ 、r0、 μ 、 β(s)和(r0,s) 已知，則方程組完備。下面確定波傳播系數 β(s)和彈性空腔半徑上的壓力邊界條件(r0,s)。廣義Maxwell黏彈性模型可以用多個Maxwell體并聯描述，如圖2所示，其Laplace域表示可寫為[15]：

圖2 廣義Maxwell體模型Fig.2 The generalized Maxwell element model

式中： θi=ηi/Ei、ηi和Ei（i=0,1,2,···,N）分別表示第i個Maxwell體松弛時間、黏度系數和彈性模量，(s) 為Laplace域的彈性模量，彈性元件的彈性模量E0、松弛時間θ0=∞ 。

波傳播系數 β(s)和E(s)之間滿足[17]：

Muller等[4-5]通過對地下爆炸自由場測量結果的分析，認為球腔上壓力邊界條件的(r0,s)可寫為：

式 中： σd為 動 態 壓 力 峰 值，σ0s為 穩 態 空 腔 壓 力，σd+σ0s為空腔壓力峰值，t0為壓力衰減的時間常數，H(t) 為 亥維賽階躍函數（Heaviside function），L表示Laplace算符。

至此，把式（16）～（17）代入式（8）～（13），即可給出黏彈性球面波各參數的Laplace域的形式；通過對式（8）～（13）進行Laplace逆變換，即可求解出時域的粒子速度vr(r,t) 、粒子位移、ur(r,t) 徑向應力σr(r,t) 、切向應力σθ(r,t) 、徑向應變εr(r,t)、切向應變εθ(r,t) ；把這些時域的解代入式（1）～（5），即可給出某個觀測區域r0≤r≤r1地下爆炸的輻射能量。

3 地下爆炸輻射地震波能量的計算實例

以黃土中空腔爆炸為例說明地下爆炸輻射地震波能量的特征。參照文獻[16]，取空腔半徑為0.025 m，空腔邊界壓力為 σr0=?(e?t/t0+1) MPa，其中t0=θ1。地震波的觀測區域分別為r0=0.025 m，r1=0.096 5，0.146 5，0.200 0 m，黃土黏彈性參數見表1。

表1 黃土黏彈性參數[18]Table1 The visco-elastic parameters of loess[18]

圖3給出黃土中地下爆炸輻射地震波能量的傳播演化。對于某個特定的r0和r1，在r0≤r≤r1的觀測區域內，r0邊界上流入的能量最終趨于穩定值W(r0,∞) ；在r1邊界上，當波到達此處，能量開始流出，最終趨于穩定值W(r1,∞) ；自r0邊界上開始流入能量后，勢能Ep(r0,r1,t) 和動能Ek(r0,r1,t)即逐漸增加，當波開始流出r1邊界時，Ep(r0,r1,t)和Ek(r0,r1,t) 突然減少，最終Ep(r0,r1,t) 趨于穩定值（如式（6）），Ek(r0,r1,t)趨于0，Ek(r0,r1,t) 為0時說明波已完整地從r0≤r≤r1觀測區域離開；自r0邊界上開始流入能量后，Ed(r0,r1,t)即逐漸增加，最終趨于穩定值Ed(r0,r1,∞)。Ep(r0,r1,∞)和Ed(r0,r1,∞) 分別表示波完整通過r0≤r≤r1區域后留在該區域介質的勢能和材料吸收的耗散能。

圖3 黃土中地下爆炸輻射地震波能量隨觀測區域的變化Fig.3 The variation of the energy of the radiated seismic wave from underground explosion in loess with the observation region

當觀測區域變大后，比如r0不變，隨著r1由0.096 5 m增大至0.2 m，r1邊界上流出的能量W(r1,∞)越來越低；勢能Ep(r0,r1,t) 、動能Ek(r0,r1,t) 和耗散能Ed(r0,r1,t)在觀測區域重合的部分其變化曲線也是重合的，觀測區域越大則獲得的能量傳播信息越多；Ed(r0,r1,t) 隨著觀測區域的變大，Ed(r0,r1,∞)隨之增大。

在地下爆炸實驗研究中，一般是通過測量不同爆心距處粒子速度、應力等參量表征應力波的傳播特性[7,19-22]，而地下爆炸輻射的地震波能量不能直接進行測量，但可以通過測得的粒子速度、應力等物理量計算。比如，通過同時測量某球面處徑向應力和徑向粒子速度波形，計算得到該球面處流入能量W(r,t)；通過測量r0≤r≤r1區域內不同位置的有限個粒子速度波形，利用一定的物理和數學處理方法[23]，獲得該區域的粒子速度場，從而可由式（3）計算得到Ek(r0,r1,t) ；通過測量r0≤r≤r1區域內不同位置有限個徑向應力、切向應力和徑向粒子速度波形，結合一定的物理和數學處理方法構造徑向應力 σr(r,t)、切向應力σθ(r,t) 、徑 向 應 變εr(r,t) 和 切 向 應 變εθ(r,t) 的 場，從 而 利 用 式（4）計 算 得 到Ep(r0,r1,t) ；在 獲 得W(r,t)、Ek(r0,r1,t)、Ep(r0,r1,t) 之后，即可求得耗散能Ed(r0,r1,t)。

4 地下爆炸輻射地震波能量的傳播特征

4.1 不同球面處流入能量的傳播特征

為研究黏彈性固體中地下爆炸在不同球面處流入能量的傳播特征，把觀測區域取為r=0.025～10 m。圖4給出黃土中不同半徑球面處最終流入能量W(r,∞) 的變化。可看出，W(r,∞) 隨半徑r增加而減小，并呈現出不同的變化規律。結合數據的變化特點，人為把0.025～10 m的空間范圍劃分為4個區域，并分段對W(r,∞)進行擬合。可以發現，在圖中R1、R3、R4區域可近似用冪函數進行擬合，而在R2區域用指數函數擬合更符合W(r,∞)的變化規律。若用R1～R4這4個區域數據的擬合公式分別反推空腔壁上的能量參數，則預估的能量分別為真實輸入能量的0.93、0.60、189、2.71倍，可見若用某局部觀測區域的實測數據預估源區或遠區波動參數時仍需慎重。

圖4 黃土不同半徑球面處最終流入能量W (r,∞)的變化Fig.4 The changes in the inflow energyW (r,∞)at different radii of the sphere in loess

在無黏性的理想彈性固體中，其耗散能為零，則可通過式（7）得到理想彈性材料中地下爆炸不同球面處流入能量：

圖4同時給出了不考慮介質黏性時W(r,∞)的傳播圖像，可看出，理想彈性固體中不同球面處流入的能量隨著波傳播距離的增加而逐漸降低，并在傳播至幾倍空腔半徑處即基本保持不變。可見，黏彈性材料中地下爆炸在不同球面處流入能量的傳播特征要比理性彈性材料復雜的多。

4.2 地震波輻射動能的傳播特征

由公式（3）可知，為求得r0≤r≤r1區域內的輻 射 動 能Ek(r0,r1,t) ，需 先 求 得 此 區 域 內vr(r,t)的空間分布。按照第2節的方法，圖5給出了黃土中不同時刻粒子速度vr(r,t)的空間分布。可以看出，當t=100 μs時，在空腔壁上激發的應力波還未完全進入觀測區域，隨著時間的增加，波陣面一直向遠離爆心方向傳播，同時波長逐漸增加、幅度逐漸降低；當t=300 μs時，在觀測區域內可以觀察到完整的應力波。

圖5 黃土中不同時刻粒子速度的空間分布Fig.5 The spatial distribution of the particle velocity at different times in loess

利用式（3）計算給出不同時刻整個空間中總動能的空間分布，如圖6所示，圖中平臺對應的即是地震波輻射動能穩態值Ek。當在觀測區域內能觀測到完整應力波時（取t=300～9 000 μs），把波陣面對應位置rWF設為橫坐標、地震波輻射動能穩態值Ek設為縱坐標，則可以給出地震波輻射動能穩態值Ek隨波前位置rWF的變化曲線，如圖7所示。與W(r,∞) 類似，在能觀測到完整應力波空間分布的區域內，地震波輻射動能穩態值Ek隨著波傳播距離rWF的增加而減少。結合數據變化特征，在離爆心較近的區域，地震波輻射動能穩態值隨波前位置的變化曲線可近似用指數函數擬合，而在較遠區域，可用分段冪函數進行擬合。

圖 6 黃土中不同時刻輻射動能的空間分布Fig.6 The spatial distribution of the radiated kinetic energy at different times in loess

圖7 黃土中地震波輻射動能隨波前位置的變化Fig.7 The variation of the radiated kinetic energy of seismic wave with the position of wave front in loess

5 結 論

由上述分析得到以下幾點結論：

（1）利用建立的黏彈性固體中地下爆炸輻射地震波場的計算方法，給出了地下爆炸有限黏彈性觀測區域中流入（出）能量、動能、應變能、耗散能的結果，數值計算結果和理論分析結果相符，說明本文計算方法的正確性。

（2）在黏彈性介質中，某球面處流入的能量W(r,∞)隨半徑增加而逐漸降低，同時呈現出復雜的衰減規律，大部分區域W(r,∞) 的 衰減規律可用冪函數擬合，而部分區域W(r,∞)的衰減更符合指數函數衰減特征；在理想彈性介質中，W(r,∞)在幾倍彈性半徑外即可穩定到某一定值。

（3）在某一固定的有限觀測區域內，當觀測時間足夠長時，勢能和耗散能均趨于某一定值，輻射動能趨于零。

（4）當有限的觀測區域能完整地容納地震波時，利用粒子速度的空間分布可以計算給出此時地震波輻射動能的空間分布，地震波輻射動能的穩態值隨波傳播距離的增加而減少，且衰減規律比較復雜，總體上可以用指數函數和冪函數進行分段擬合。