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

盧 強(qiáng)，丁 洋，劉赟哲，唐仕英，郭志昀，王占江

（西北核技術(shù)研究所，陜西 西安 710024）

研究地下封閉爆炸時(shí)通常把爆炸源簡(jiǎn)化為球形爆源，這樣地下爆炸問(wèn)題可簡(jiǎn)化為球面波問(wèn)題。地下爆炸能量與源區(qū)介質(zhì)發(fā)生強(qiáng)烈的耦合，最終以波的形式向周圍介質(zhì)輻射地震波能量。為表征地下爆炸或地震激發(fā)的震源強(qiáng)度，學(xué)術(shù)界提出了多種模型[1-3]。

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

Choy等[2]和Boatwright等[3]利用大量的地震事件，基于能量的概念建立了地震波能量和能量震級(jí)之間的關(guān)系，以此表征震源強(qiáng)度。袁乃榮等[8]和李贊等[9]利用國(guó)家地震臺(tái)網(wǎng)和全球地震臺(tái)網(wǎng)的寬頻地震波數(shù)據(jù)對(duì)地震波的能量震級(jí)進(jìn)行研究，指出地震波對(duì)傳播路徑區(qū)域造成的破壞直接與地震釋放的地震波能量相關(guān)，可見(jiàn)地下爆炸向介質(zhì)中耦合的地震波能量是表征地下爆炸震源強(qiáng)度的另外一個(gè)重要參數(shù)。

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

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

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

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

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

力Fr(r,t)在 半徑r的球面上做的功W(r,t) 即 是流入此球面的總能量。根據(jù)能量守恒，能量流經(jīng)r0≤r≤r1的區(qū)域時(shí)，流入能量W(r0,t) 轉(zhuǎn)化為此區(qū)域內(nèi)所有質(zhì)點(diǎn)的動(dòng)能Ek(r0,r1,t)（即輻射動(dòng)能）、介質(zhì)應(yīng)變能（也稱內(nèi)能或勢(shì)能）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時(shí)刻整個(gè)區(qū)域r0≤r≤r1內(nèi)所有質(zhì)點(diǎn)動(dòng)能和應(yīng)變能密度函數(shù)的空間積分，則有：

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

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

考慮t→∞ 的特殊情況，此時(shí)波完全離開(kāi)r0≤r≤r1區(qū)域，則此區(qū)域內(nèi)質(zhì)點(diǎn)動(dòng)能Ek(r0,r1,∞)=0 。按照黏彈性球面波的解，式（4）可以寫(xiě)為[15-16]：

式中： σ0s為穩(wěn)態(tài)空腔壓力，Ge為黏彈性材料的靜態(tài)剪切模量。從式（6）可以看出，穩(wěn)態(tài)空腔壓力 σ0s形成的勢(shì)能主要分布在幾倍空腔半徑的范圍內(nèi)，比如當(dāng)r1=10r0時(shí)，此觀測(cè)區(qū)域儲(chǔ)存了99.9%的勢(shì)能。

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

2 黏彈性固體中地下爆炸輻射地震波場(chǎng)的計(jì)算方法

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

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

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

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

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

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

Muller等[4-5]通過(guò)對(duì)地下爆炸自由場(chǎng)測(cè)量結(jié)果的分析，認(rèn)為球腔上壓力邊界條件的(r0,s)可寫(xiě)為：

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

至此，把式（16）～（17）代入式（8）～（13），即可給出黏彈性球面波各參數(shù)的Laplace域的形式；通過(guò)對(duì)式（8）～（13）進(jìn)行Laplace逆變換，即可求解出時(shí)域的粒子速度vr(r,t) 、粒子位移、ur(r,t) 徑向應(yīng)力σr(r,t) 、切向應(yīng)力σθ(r,t) 、徑向應(yīng)變?chǔ)舝(r,t)、切向應(yīng)變?chǔ)纽?r,t) ；把這些時(shí)域的解代入式（1）～（5），即可給出某個(gè)觀測(cè)區(qū)域r0≤r≤r1地下爆炸的輻射能量。

3 地下爆炸輻射地震波能量的計(jì)算實(shí)例

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

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

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

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

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

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

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

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

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

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

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

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

4.2 地震波輻射動(dòng)能的傳播特征

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

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

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

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

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

5 結(jié) 論

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

（1）利用建立的黏彈性固體中地下爆炸輻射地震波場(chǎng)的計(jì)算方法，給出了地下爆炸有限黏彈性觀測(cè)區(qū)域中流入（出）能量、動(dòng)能、應(yīng)變能、耗散能的結(jié)果，數(shù)值計(jì)算結(jié)果和理論分析結(jié)果相符，說(shuō)明本文計(jì)算方法的正確性。

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

（3）在某一固定的有限觀測(cè)區(qū)域內(nèi)，當(dāng)觀測(cè)時(shí)間足夠長(zhǎng)時(shí)，勢(shì)能和耗散能均趨于某一定值，輻射動(dòng)能趨于零。

（4）當(dāng)有限的觀測(cè)區(qū)域能完整地容納地震波時(shí)，利用粒子速度的空間分布可以計(jì)算給出此時(shí)地震波輻射動(dòng)能的空間分布，地震波輻射動(dòng)能的穩(wěn)態(tài)值隨波傳播距離的增加而減少，且衰減規(guī)律比較復(fù)雜，總體上可以用指數(shù)函數(shù)和冪函數(shù)進(jìn)行分段擬合。