爆破地震的數(shù)值模擬及爆破振動規(guī)律分析

張智超，劉漢龍，陳育民，王維國

（1.河海大學(xué) 巖土力學(xué)與堤壩工程教育部重點實驗室，江蘇 南京210098；2.河海大學(xué) 土木與交通學(xué)院，江蘇 南京210098）

0 引言

土體地震動力響應(yīng)一直是巖土領(lǐng)域的研究重點.以往所開展的地震動力試驗，主要是通過動三軸等試驗的結(jié)果來分析土體的地震液化特性和動力響應(yīng)［1］，以及振動臺、剪切模型箱和離心機等配合使用的手段來還原土工建筑或土體結(jié)構(gòu)遭受天然地震時的受力狀態(tài)［2－5］.然而，常規(guī)1 g室內(nèi)試驗難以滿足土體原型應(yīng)力水平，而離心機試驗雖然能夠還原土體的實際應(yīng)力狀態(tài)，但也存在一定缺陷，如離心應(yīng)力場不均勻，在沿著離心力作用的半徑方向上，離心力與半徑大小成正比，這就給原位應(yīng)力場的模擬帶來一定誤差，且離心機造價昂貴，技術(shù)復(fù)雜，很大程度上制約了其推廣使用.

因此，王蘭民［6］采用現(xiàn)場微差爆破試驗的方法，創(chuàng)造出與天然地震相似的持續(xù)振動狀態(tài)，以此來分析場地的地震震陷特征，取得了不錯的效果；此外，通過爆破手段來模擬天然地震的試驗研究也不斷地為學(xué)者們所嘗試［7－10］，以試圖避免前述室內(nèi)實驗和模型試驗的缺陷，從原位的角度來研究土體力學(xué)行為.

目前，對此類利用微差爆破模擬天然地震的研究，主要集中在現(xiàn)場試驗方面.然而，現(xiàn)場試驗的一些弊端，如對場地要求較高，數(shù)據(jù)監(jiān)測范圍有限，較為耗費人力、物力，難以進行大量重復(fù)等，都在一定程度上制約了爆破模擬天然地震的研究.

作者利用數(shù)值手段，對文獻［6］中的微差爆破試驗進行了模擬，將計算結(jié)果與實測值進行比較，結(jié)合理論分析，為數(shù)值結(jié)果的正確性提供保障；在此基礎(chǔ)上，對地層中的爆破振動規(guī)律進行了分析，豐富了爆破地震動的研究，對將來利用現(xiàn)場爆破試驗來模擬天然地震的研究工作具有一定指導(dǎo)意義.

1 現(xiàn)場試驗

試驗［6］位于甘肅某黃土場地，場地深度0～1 m為植耕土，1～15 m為風成黃土，15～28 m為次生黃土，28 m以下為卵礫石層.

炮點布局如圖1所示，作為爆源的30個藥包埋設(shè)在以場地中心為圓心、半徑為15 m的圓周上，等間距分布，埋深23 m.每個藥包裝有40 kg二號巖石炸藥，軸對稱上的兩個藥包作為一對同時起爆.如圖1所示，炮點邊上的標序表示引爆順序，而炮點邊上無標序則表示此處藥包引爆失敗，模擬中將忽略之.

圖1 藥包布局及引爆順序［6］Fig.1 Distribution of explosive charges and their detonation or der［6］

場地中心的地表布置了三分相標準強震儀，用于捕捉試驗過程中的爆破振動加速度時程.

2 數(shù)值模型

2.1 幾何模型

根據(jù)試驗現(xiàn)場情況，數(shù)值模型采取直徑50 m、高度30 m的圓柱體計算區(qū)域，如圖2所示.分別將深度0～1 m、1～15 m、15～28 m和28～30 m的土體劃分為不同材料，用不同顏色表示.x向?qū)?yīng)現(xiàn)場的WE（東西）向，y向?qū)?yīng)于現(xiàn)場的NS（南北）向，z向?qū)?yīng)于現(xiàn)場的UD（豎直）向.30個藥包在深度23 m、半徑15 m的圓環(huán)上等間距分布，每個藥包質(zhì)量為40 kg.

圖2 有限元模型Fig.2 FEMmodel

2.1.1 有限元網(wǎng)格劃分

模型中的所有單元均為六面體八節(jié)點實體單元，并在藥包附近進行了網(wǎng)格細化，共計527 085個單元和545 184個節(jié)點；同時，毗鄰單元間的網(wǎng)格尺寸比都小于1∶4，并對單元大小進行限制，通過適當?shù)木W(wǎng)格漸變，將單元尺寸由炸藥單元的0.2 m逐漸過度到最外圍土體單元的1 m，即遠離炸藥的單元網(wǎng)格逐漸變得稀疏.這在一定程度上保證了計算精度，防止單元退化.

2.1.2 邊界條件

模型上表面設(shè)置為自由邊界，以示與空氣接觸，四周與底部設(shè)置為無反射邊界，以真實地反映波在此的透射情況.此外，藥包與底部邊界距離為7 m，與圓周邊界距離10 m，在一定程度上也保證了足夠的計算范圍，消除了應(yīng)力波的反射.

2.1.3 求解時間

孫軍杰［10］對實測加速度結(jié)果進行了處理，合理地截掉了6對未爆爆源所導(dǎo)致的振動中斷時段，得到了疊續(xù)合成的振動時程.因此，數(shù)值計算中也將對這6段未爆爆源的微差時段進行忽略，只考慮成功引爆的微差時段，總計算時間為6 s.根據(jù)試驗中的微差時間，計算爆破時刻如表1所示.

表1 數(shù)值計算的爆破順序Tab.1 Blasting sequence in the numerical simulation

式中：P為壓力；φ為內(nèi)摩擦角；J2為偏應(yīng)力張量第二不變量；K（θ）為張量平面角的函數(shù)；c為黏聚力；A HYP為決定修正后的Mohr－Coulo mb屈服面和標準Mohr－Coulomb屈服面相似度的參數(shù).

當AHYP＝0時，上式表示標準 Mohr－Coulomb屈服面，當AHYP采用較大值時，修正后的屈服面明顯偏離標準屈服面，對于數(shù)值模擬來講，A HYP的取值應(yīng)該小于c cotφ，一般按照下式進行選?。?/p>

2.2 物理模型

2.2.1 土體本構(gòu)

LS－DYNA中的 MAT＿FH WA＿SOIL［11］本構(gòu)可以模擬諸多土體特性，如爆炸荷載下的應(yīng)變率效應(yīng)、應(yīng)變硬（軟）化以及孔壓等，是較為合適的計算土體爆炸響應(yīng)的模型.該模型在爆炸分析方面的優(yōu)越性已在文獻［7］中得到充分驗證.因此，作者選取其來描述現(xiàn)場爆破試驗中土體的力學(xué)行為.

該模型將標準Mohr－Coulo mb屈服面修正為光滑的曲面，并且垂直于壓力軸，提高了數(shù)值計算的有效性和穩(wěn)定性.屈服函數(shù)為：

圖3為標準Mohr－Coulomb屈服面與修正的屈服面對比.除了在低應(yīng)力區(qū)外，二者幾乎一致.

深度0～28 m的黃土層采用MAT＿FH WA＿SOIL本構(gòu).根據(jù)現(xiàn)場探測的剪切波速和壓縮波速［6］，按 照 公 式算出各土層的剪切模量和體積模量；由于黃土含水率較低，將孔壓系數(shù)Ksk設(shè)置為0，即不考慮孔壓發(fā)展，進行總應(yīng)力分析；其余參數(shù)根據(jù)現(xiàn)場實測數(shù)據(jù)并結(jié)合相關(guān)經(jīng)驗得出［11，13］，參數(shù)值見表2.

圖3 標準與修正后的Mohr－Coulomb屈服面對比Fig.3 The standard and modified Mohr－Coulomb yield surface

表2 MAT＿FHWA＿SOIL模型參數(shù)Tab.2 Parameters of MAT＿FHWA＿SOIL

對于第28～30 m的卵礫石層，則采用線彈性本構(gòu)，密度2 000 kg／m3，楊氏模量5e8 Pa，泊松比0.25.

2.2.2 炸藥材料模型

采用高能炸藥燃燒模型和J WL狀態(tài)方程來模擬現(xiàn)場試驗中所使用的二號巖石炸藥［13］.JWL狀態(tài)方程能夠精確描述爆炸過程中爆轟產(chǎn)物的壓力、體積、能量特性，其關(guān)系式為

式中：A，B，R1，R2，ω 為材料參數(shù)；E01為爆轟產(chǎn)物單位體積的內(nèi)能；V為單位體積裝藥產(chǎn)生的爆轟產(chǎn)物的體積；P為爆炸爆轟壓力.炸藥材料參數(shù)取值分別為［14］：密度1 000 kg／m3，爆轟速度3 600 m／s，CJ壓力3.24 GPa，A＝214 GPa，B＝0.182 GPa，R1＝4.15，R2＝0.95，ω＝0.38，E01＝4.5 GPa，V＝1.

2.3 結(jié)果分析

在爆破計算前，利用LS－DYNA本身的動力松弛算法對土體初始自重靜應(yīng)力進行計算，結(jié)果如圖4所示.可以看出，自重結(jié)果合乎規(guī)律，為其后的爆破分析提供了準確的初始應(yīng)力條件.

圖4 初始自重應(yīng)力云圖側(cè)視圖Fig.4 Side view of initial gravity stress nephogram

2.3.2 模擬與實測地表加速度時程的比較

為了驗證模擬結(jié)果的準確性，分別將模擬得到的x、y和z向加速度時程與所對應(yīng)E W、SN和UD向加速度時程實測值［10］進行比較，圖5是對比圖.由圖可以看出，模擬的加速度時程出現(xiàn)了9個較為明顯的峰值群，與現(xiàn)場9對藥包依次爆破后激發(fā)的地震動相對應(yīng).由于地震本身就是一種雜亂無章的隨機振動現(xiàn)象，加之現(xiàn)場監(jiān)測數(shù)據(jù)的離散性，以及數(shù)值計算的近似性，無法做到計算曲線與實測曲線的完全一致，但由圖仍可看出，模擬的x、y向加速度時程與現(xiàn)場實測的EW、SN向加速度時程在幅值、曲線形狀和峰值出現(xiàn)的時刻都吻合得較好，從一定程度上說明了數(shù)值手段的準確性；同時還可以看出，微差爆破在場地中心所產(chǎn)生的連續(xù)的爆破地震波時程曲線在宏觀形態(tài)上與典型的天然地震波曲線具有比較好的相似性，因此，可以通過爆破手段來對原位土體場地的地震響應(yīng)進行分析，其反映出的總體振動趨勢大體上可以被用來模擬天然地震.

值得一提的是，對于模擬的z向加速度時程與實測豎向加速度時程，雖然在曲線形狀和峰值出現(xiàn)時刻上較為吻合，但峰值大小差別相對較大，模擬的峰值普遍接近2 g，且在第一組峰值群有著不合常理的躍遷；而實測峰值普遍不到1 g，模擬值幾乎是實測峰值的2倍.分析原因可能是由于計算參數(shù)不夠精確以及有限元計算方法的近似性等造成的，加之實際環(huán)境的干擾，如土質(zhì)不均勻、監(jiān)測設(shè)備精度有限等，與數(shù)值計算中完全理想的情況有所不同，就更容易導(dǎo)致實測值與模擬值的差異，但我們?nèi)钥梢酝ㄟ^數(shù)值手段預(yù)測爆破震動的大致規(guī)律.

2.3.3 爆破振動隨地層深度的變化規(guī)律

基于上述結(jié)論，選取了從深度23 m——即炸藥埋深至地表共10個不同深度的節(jié)點A～J作為參考，分析爆破振動變化規(guī)律.參考節(jié)點位置見圖6.

圖6 參考節(jié)點示意圖Fig.6 Location of reference nodes

圖7 為不同深度x、y和z向加速度時程，子圖中的上圖為時程全紀錄，下圖將加速度特征局部放大，以便更為清晰地展示其變化規(guī)律.由圖7可以看出，隨著埋深的減小、爆心距的增加，各節(jié)點x、y向（水平向）加速度幅值逐漸減小，而z向（豎向）加速度幅值則先增大、后減小，到了近地表又趨于穩(wěn)定.之所以產(chǎn)生這種現(xiàn)象，原因分析如下：爆炸引起的球形應(yīng)力波在土體空間中擴張，爆炸振動最為劇烈的方向近似在爆心與參考質(zhì)點的連線方向上，即球形波陣面的外法線方向，因此，對于位于一對同時引爆的藥包連線中點上的質(zhì)點來說，z向為球形波陣面的切線方向，振動強度趨于最小，所以，當參考節(jié)點埋深減小，即節(jié)點上移、離開了藥包的連心線中點后，z向的振動強度分量增大，這導(dǎo)致了z向加速度幅值的首次增大；進而，隨著節(jié)點深度的繼續(xù)減小，爆心距也逐漸增大，由于爆破振動強度與爆心距呈反關(guān)系，這就導(dǎo)致了z向（豎向）加速度幅值產(chǎn)生了減小的趨勢；但藥包附近存在自由面時，土體中的應(yīng)力狀態(tài)將變得較為復(fù)雜［15］，自由面的存在，會增強爆破效果，較之自由面切線方向（x、y向）來說，自由面反射的增強作用對自由面法線方向（z向）的影響更大，因此，這就導(dǎo)致臨近自由面時，反射波的放大作用和爆心距增加而導(dǎo)致的減弱作用大致相互抵消，所以雖然爆心距增加了，但質(zhì)點豎向振動強度并不隨著爆心距的增加而減小，而是基本維持穩(wěn)定；另一方面，自由面切向的質(zhì)點運動特征受到自由面反射的影響相對較小，因而切向的振動幅值并未因為臨近自由面而增強，而僅是隨著爆心距的增加呈減小的趨勢.

質(zhì)點速度與加速度的變化規(guī)律呈現(xiàn)一致性，即從炸藥所在的埋置深度向上統(tǒng)計，隨著深度的減小、爆心距的增大，質(zhì)點x、y向速度幅值呈現(xiàn)逐漸減小的趨勢，而質(zhì)點z向速度幅值則先增大、后減小，最后在臨近自由面處又趨于穩(wěn)定.分析原因與前述類似，在此不再贅述.

3 結(jié)論

對微差爆破模擬天然地震的現(xiàn)場試驗進行了數(shù)值分析，可以得到以下幾點結(jié)論.

（1）數(shù)值手段可以很好地還原微差爆破作用下的土體質(zhì)點振動特征.數(shù)值手段簡便易行，可靈活地針對不同的研究目標和內(nèi)容進行調(diào)整，因此，開展爆破地震動的數(shù)值分析工作，對將來利用現(xiàn)場爆破來模擬天然地震的試驗設(shè)計、施工及爆破振動預(yù)測都具有一定指導(dǎo)意義.

（2）從數(shù)值計算的角度驗證了采用爆破手段來模擬天然地震在一定程度上是可行的，可以利用爆破地震波創(chuàng)造出與天然地震相似的震動環(huán)境，進而較為符合實際地進行一系列原位土體地震響應(yīng)研究，避免其他室內(nèi)試驗和模型試驗的缺陷.

（3）在本文的算例中，自炸藥埋深向上，隨著深度的減小、爆心距的增加，爆破水平向加速度幅值和速度幅值逐漸減小，而豎直向加速度幅值和速度幅值則先增大、后減小，到了近地表又因自由面反射效應(yīng)的影響而基本維持穩(wěn)定.

［1］ 何昌榮.動模量和阻尼的動三軸試驗研究［J］.巖土工程學(xué)報，1997，19（2）：39－48.

［2］ 吳俊賢，倪至寬，高漢棪.土石壩的動態(tài)反應(yīng)：離心機模型試驗與數(shù)值模擬［J］.巖石力學(xué)與工程學(xué)報，2007，26（1）：1－14.

［3］ 汪明武，Susu mu IAI，Test uo TOBITA.液化場地堤壩地震響應(yīng)動態(tài)土工離心試驗及模擬［J］.水利學(xué)報，2008，39（12）：1346－1352.

［4］ 黃春霞，張鴻儒.大型疊層剪切變形模型箱的研制［J］.巖石力學(xué)與工程學(xué)報，2006，25（10）：2128－2134.

［5］ 李敏霞，楊澤群，陳建秋.地震模擬振動臺技術(shù)的開發(fā)與應(yīng)用［J］.世界地震工程，1996，2（1）：49－54.

［6］ 王蘭民，孫軍杰，徐舜華，等.爆破模擬地震動條件下黃土場地震陷研究［J］.巖石力學(xué)與工程學(xué)報，2008，27（5）：913－921.

［7］ Lee，WAYNE Y.NUMERICAL MODELING OF BLAST－INDUCED LIQUEFACTION ［D］.Brigham Young University，Ph D Thesis.August 2006.

［8］ PAT HIRAGE K S.Critical assessment of the CANLEX blast experi ment to facilitate a development of an in－situ liquefaction methodology using explosives［D］.Depart ment of Civil Engineering，Degree of Master of Science，2000.

［9］ ASHFORD S A，ROLLINS K M.Blast－Induced Liquefaction for Full－Scale Foundation Testing［J］.Jour nal of Geotechnical and Geoenviron mental Engineering.2004，130（8）：798－806.

［10］孫軍杰.黃土場地震陷與樁基負摩阻力現(xiàn)場試驗研究 ［D］.蘭州：蘭州大學(xué)土木工程與力學(xué)學(xué)院博士學(xué)位論文，2010.

［11］Lewis B A.Manual for LS－DYNA Soil Material Model 147［R］.Federal Highway Administration，2004.

［12］徐芝綸.彈性力學(xué) ［M］.北京：人民教育出版社，1983.

［13］Anon.Ls－dyna key word users’menual（Version 971／Release 4）［M］.s.l.Liver more Soft ware Technology Cor poration，2009.

［14］劉紅巖，王新生.地表垂直爆破震動速度的數(shù)值計算［J］.爆炸與沖擊，2008，28（3）：255－270.

［15］亨利奇.爆炸動力學(xué)及其應(yīng)用 ［M］.熊建國，譯.北京：科學(xué)出版社，1987.