核爆炸對地沖擊作用下土體運動特性等效模擬

榮吉利， 宋逸博,2， 王璽， 郭振， 項大林， 吳志培,2

(1.北京理工大學 宇航學院, 北京 100081; 2.中國運載火箭技術研究院, 北京 100076;3.北京航天發射技術研究所, 北京 100076; 4.北京宇航系統工程研究所, 北京 100076)

0 引言

核武器由于其巨大的毀傷能力及戰略價值，得到了快速的發展。以往對于地面及地下結構在核爆襲擊下的強度校核，主要采用沖擊波超壓峰值和正壓作用時間或沖量值作為創傷指示參數[1]，該方法過于單一且不能完全反映出結構所遭受的威脅。隨著武器技術的發展以及對核爆效應的了解日益深入，由爆炸產生的對地沖擊效應引起了廣泛的關注，并被認為是對地下結構最大的威脅。1968年4月26日，在美國內華達地下進行的1 200 kt當量的Boxcar核爆炸產生了強大的對地沖擊效應及巖土的不可逆變形，引起最遠震中距6 100 m處地表出現斷層破裂，斷層錯動的最大垂直相對位移為1 m[2]. 美軍在1964至1976年間對土體介質及巖石場地進行了一系列大當量爆炸試驗以對核爆炸對地沖擊效應進行論證分析[3-7]，并將研究成果擬合成經驗公式，整理在TM 5-858設施抗核武器效應設計手冊中的TM 5-858-2武器效應分冊中[8]。Pathak等[9]利用一維彈性波傳播模型，對核武器空爆引起的地運動位移進行了研究，為工程提出了一種預測地面位移的方法。我國學者也對爆炸地沖擊效應進行了大量的研究[10-12]，形成了較為完備的計算方法。但由于對爆炸引起的地運動規律缺乏有效且統一的評估方法，不同計算方式對地沖擊效應的計算結果存在較大差異，因此難以得到一致結論。此外大當量的對地沖擊試驗工程巨大，難以大范圍開展獲取充足的試驗數據，因此對地沖擊效應的認識仍存在一定的局限性。

隨著計算機模擬技術的發展，對核爆炸對地沖擊效應進行準確的數值仿真成為了可能，林大超等[13]提出地面隨機平穩噪聲的時變曲線是一條分型曲線的假設，給出了爆炸所致地面運動的仿真方法和過程；唐廷等[14]利用有限元分析軟件LS-DYNA對地沖擊后地下預置裂紋拱結構進行了抗爆炸分析，以研究受損的硐室抵抗二次爆炸的能力；Ma等[15]、李永勝等[16]分析了地下硐室小當量化爆下巖土響應，其結果與試驗結果擬合較好，但相關工作研究的爆炸當量較小，無法應用于核爆等大當量爆炸仿真；王小盾等[17]、苗青等[18]針對天津8·12爆炸事故，考慮了爆炸空氣沖擊波與爆炸對地沖擊波聯合作用，利用有限元分析軟件Abaqus對網架結構動力響應進行分析，證明了爆炸對地沖擊效應對結構的影響不可忽視，但由于缺少相關實測數據，仿真結果缺乏對比驗證，模型的準確性和可靠性難以保證，且沒有考慮空氣沖擊波激發的感生沖擊波的作用。

綜上所述，對于大當量爆炸下土體運動準確仿真方法的建立與模型的選取亟需探討。本文考慮了感生沖擊波與直接沖擊波的疊加作用，通過使用流體- 固體耦合技術，得到了由空氣沖擊波激發的感生沖擊波，研究了二者的貢獻程度，證明了爆炸引起的對地沖擊效應是感生沖擊波和直接沖擊波聯合作用的產物。提出了以Shock Hugoniot絕熱關系描述的Mie-Gruneisen狀態方程以及Drucker-Prager(D-P)彈塑性準則來描述黏土的行為，使用JWL狀態方程描述爆轟產物的分析方法，通過有限元分析軟件AUTODYN以美軍代號LN302的對地沖擊試驗[3]為例進行數值仿真模擬。將分析結果與美軍試驗數據和美軍手冊TM5-858-2的公式計算結果進行對比，吻合良好，驗證了本模型的準確性和有效性，為核爆炸對地沖擊效應的仿真模擬提供了有效的方法。

1 核爆炸對地沖擊試驗

1968年，美軍實施了LN302對地沖擊試驗中代號為PRAIRIE FLAT的子試驗[3]，以研究黏土地質下核爆炸對地沖擊效應。試驗在廣闊的平原區域進行，該地區由河床沉積而成，沒有巖層，主要地質成分為黏土。由于核爆炸與化學爆炸主要區別在于爆炸反應發生的初期階段經歷了不同的物理過程以及有無后續輻射產物，而對于爆炸對地沖擊波的產生幾乎沒有差別，因此美軍使用高能炸藥梯恩梯(TNT)作為核爆炸的等效爆炸源[12]。在該試驗中，等效TNT當量為500 t，成球狀堆疊置于地表。根據預測的超壓峰值，在距離爆心25.603～350.520 m的距離范圍內布置了地運動測量儀器，深度0.457～9.144 m. 儀器封裝在鋁制圓柱罐內，固定置于地下，其具體布置方案如表1所示。

表1 地運動測量儀器布置陣列[3]

核武器在地表爆炸會在空氣中產生向外高速擴張的球面空氣沖擊波，空氣沖擊波掠過地表，其激波前緣會黏連并牽引地表的土體與間隙的空氣產生應力波擾動，稱為感生沖擊波，而由于爆炸釋放的能量直接向地面耦合，在地下產生從彈坑處向外擴散的應力波稱為直接沖擊波[11]。美軍通過理論分析與試驗驗證，將對地沖擊作用區域從爆心向外劃分為超震區和超前區。距離爆心較近的區域為超震區，由于爆心附近空氣沖擊波速度遠大于爆炸引起的直接沖擊波速度，因此位于超震區內的物體先接收到空氣沖擊波及其激發的感生沖擊波，之后接收到直接沖擊波。隨著距離的增加，空氣沖擊波速度迅速衰減，而直接沖擊波在地下傳播速度不變，直接沖擊波逐漸追趕并超越感生沖擊波，以二者同時到達的距離為界，之后進入超前區。位于超前區的物體先感受到直接沖擊波，后感受到沖擊波和感生沖擊波。

試驗測得的正向速度峰值數據如表2所示，表中垂直速度和水平速度的正方向相對于爆心分別為向上和向外。

表2 對地沖擊正向速度峰值試驗數據[3]

2 美軍手冊TM 5-858-2經驗公式

美軍將一系列大當量爆炸對地沖擊試驗的數據結果擬合出用于計算核爆炸對地沖擊效應的經驗公式，作為手冊TM 5-858-2[8]中核爆炸對地沖擊效應的計算方法。在計算對地沖擊速度峰值響應之前，需要根據爆炸工況判斷目標位于超震區還是超前區，其判定方法如下：

(1)

式中：pso為沖擊波壓力最大值與大氣壓力po的差值(Pa)；L為與爆心的直線距離(m)；W為有效當量(kt)；u為空爆超壓沖擊波速度(m/s)；tf為從爆炸發生到直接沖擊波到達目標位置所需的時間(s)；s為目標位置距離爆心的直線距離(m)；t為時間(s)；Range為判據，當Range≥1時使用超震區公式計算，當Range<1時使用超前區公式計算。

根據PRAIRIE FLAT試驗報告中的設置工況和黏土的實測參數，確定公式中基本參數為當量W=0.5 kt，載荷下土壤壓縮波速cL=130 m/s，s為布置有速度計的測點爆心距，分別為25.6 m、42.7 m、67.1 m、100.6 m與121.9 m，經計算均處于超震區，故使用超震區的速度峰值計算公式進行計算：

(2)

(3)

式中：uV為土壤垂直速度峰值(m/s)；uH為土壤水平速度峰值(m/s)；Z為目標深度(m)。對地下3.048 m深度的速度響應進行計算，得到的結果如表3所示，表中垂直速度和水平速度的正方向相對于爆心分別為向上和向外。

表3 3.048 m深度速度峰值手冊結果[8]

3 仿真分析模型

3.1 材料模型

3.1.1 空氣

空氣使用理想氣體狀態方程描述，

p=(γ-1)ρe，

(4)

式中：p為氣體壓力；γ為絕熱指數，γ=1+R/cv，R為一般氣體常數R0除以特定氣體的有效分子量，cv為氣體體積恒定時的比熱容；ρ為密度；e為比內能。其材料參數如表4所示。

表4 空氣材料參數[19]

3.1.2 TNT

TNT使用JWL狀態方程來描述，

(5)

式中：pe為爆轟產物壓力；v=1/ρe為比容，ρe為裝藥密度；C1、r1、C2、r2為常數；ω為絕熱常數；ee為爆轟產物的比內能。材料參數如表5所示。

表5 TNT材料參數[19]

3.1.3 土體

由于土體在應力作用下表現出較強的流動性，必須考慮其剪切應力與剪切破壞，因此本文以沖擊狀態方程結合基于摩爾- 庫倫準則的D-P彈塑性強度模型對土體進行建模，并通過控制點定義屈服應力隨壓力分段線性變化來模擬描述土體抗壓不抗拉的特性，當土體受拉時，屈服應力迅速降為0 Pa. 失效使用Hydro模型，當材料壓力小于定義的最大拉伸壓力時判定失效，參數如表6所示。

表6 土體材料參數[20]

3.2 模型建立與網格劃分

根據文獻[3]中的試驗布置情況(見圖1)建立如圖2所示有限尺寸的二維軸對稱模型，在試驗速度計布置處設置測點。模型需要考慮空氣與大地的耦合效應，建立150 m×400 m的歐拉網格模擬空氣，建立150 m×400 m的拉格朗日網格模擬土體，為與經驗公式結果進行對比，對土體進行簡化，不考慮河床線對土體水飽和度的影響。由于測點分布在全距離下，因此使用均勻網格劃分，網格尺寸為0.5 m×0.5 m. 二者設置50 m重疊以實現空氣與土體的耦合作用，在地面以上的歐拉網格填充500 t球型TNT，對球心單點起爆以模擬核彈觸地爆炸。對稱軸處為對稱邊界，將歐拉區域邊界條件設置為流出邊界，將拉格朗日區域除地表面外的邊界均設置為透射邊界以消除邊界反射，防止反射波與后續沖擊波疊加。定義拉格朗日域自接觸以模擬彈坑附近土體回落堆積，定義歐拉/拉格朗日耦合，實現應力從空氣到大地的傳播。

圖1 試驗測點布置Fig.1 Arrangement of test measuring points

圖2 模型及測點布置Fig.2 Model and arrangement of measuring points

4 結果與分析

4.1 物質運動

基于上述建模方法與材料參數建立了500 t TNT在地面爆炸的地沖擊效應仿真模型。為了觀察完整的地沖擊傳遞過程，仿真時間設置為2 s，輸出時間間隔為0.002 s，截取爆炸過程部分時間的物質運動狀態云圖如圖3所示。圖3中深藍色部分為空氣，淺藍色為黏土地質，土中呈陣列分布的點狀物為測點，中心綠色部位為TNT。提取0.02 s時局部地沖擊應力云圖及速度云圖如圖4、圖5所示。

圖3 500 t TNT觸地爆炸物質運動云圖Fig.3 Material motion contour of 500 t TNT contact explosion

圖4 0.02 s時對地沖擊應力云圖Fig.4 Stress contour of ground shock at 0.02 s

圖5 0.02 s時對地沖擊速度云圖Fig.5 Velocity contour of ground shock at 0.02 s

通過圖3可以看到，爆炸初期爆轟產物迅速膨脹形成的沖擊波會帶走爆心附近土體并形成半球形彈坑，后續的沖擊波將彈坑周圍的土體壓實、拋擲和沖蝕，使彈坑繼續擴大并形成凸起的唇緣。周圍的土體在爆炸的擠壓作用下產生了向外、向上的運動。距離爆心25.6 m的測點由于位于彈坑邊緣附近的破裂區和塑性區，其位置在土體的流動下被大幅抬升并高于地面，隨著與爆心距離的增加，土體位移幅度逐漸衰減。

4.2 對地沖擊波

根據0.02 s時對地沖擊應力云圖可知，在靠近爆心的位置，由爆炸直接形成并從彈坑處向外擴散的球面應力波前方出現了一段與空氣沖擊波前緣相連的應力波。球面應力波為直接沖擊波，在土體中以縱波波速傳播，前方的應力波是感生沖擊波。地沖擊速度云圖顯示了土體中粒子速度的變化，測點先接收到感生沖擊波帶來的地沖擊速度，直接沖擊波隨后到達。該現象與文獻[8]中結論相互印證，也說明了對于地面建筑、設備和人員來說，在沖擊波來臨的同時還有其激發的感生沖擊波在起作用，對于結構的強度校核不應只局限于超壓，還應考慮與對地沖擊作用的組合工況。

對直接沖擊波與空氣沖擊波分離情況進行分析，結果如圖6所示。

圖6 超震區與超前區分界圖Fig.6 Boundary of superseismic region and outrunning region

通過圖6可以看到，在距離爆心約158 m處直接沖擊波與空氣沖擊波前緣分離，與報告中超震區與超前區的預測分界點(約爆心距171 m處)十分接近，誤差小于10%. 隨著空氣沖擊波速度的衰減，先前生成的感生地震波與激波前緣分離，在地下以縱波波速繼續向前傳播并逐漸衰減，導致170.7 m處的測點仍能先接收到較強的感生沖擊信號，因此實際的分界點更接近預測分界點。

分別選取位于超震區67 m和超前區253 m處的典型位置測點，對比分析其壓力變化以研究直接地沖擊波和感生地沖擊波的貢獻，結果如圖7、圖8所示。由圖7、圖8可以看出：位于超震區小于3.048 m的深度上，感生沖擊波的壓力峰值超過直接沖擊波壓力峰值，成為地沖擊效應的主導部分，隨著深度的增加，直接沖擊波的壓力峰值逐漸增大，感生沖擊波的貢獻比例逐漸降低；在超前區，直接沖擊波起主導作用，壓力幅值在距離上的衰減程度遠小于感生沖擊波，且隨著深度的增加壓力峰值逐漸增大。產生這種現象的原因是由于巖土中接近地表的直接沖擊波會被地面運動釋放，而地下的壓力耗散緩慢，因此隨著深度的增加強度衰減緩慢，而感生沖擊波的壓力是由空氣沖擊波提供，導致其峰值在深度上無明顯變化。

圖7 超震區對地沖擊壓力變化曲線Fig.7 Ground shock pressure curves in superseismic region

圖8 超前區對地沖擊壓力曲線Fig.8 Ground shock pressure curves in outrunning region

4.3 速度時程

由于試驗TNT當量為500 t，根據物質云圖顯示靠近爆心的淺層測點可能存在測量失真，因此選取25.6 m遠、9.144 m深測點的速度波形與美軍手冊指導波形進行對比。使用該測點的歸一化爆距在手冊TM 5-858-2中選取波形，根據手冊指導通過匹配峰值與上升時間的方式進行擬合，二者對比結果如圖9、圖10所示。通過對比可以發現二者吻合良好，對于垂直速度波形，其先出現方向向下，幅值與向上的速度峰值相近的脈沖，隨后向上運動，之后速度緩慢下降。水平速度波形向外運動，之后速度迅速下降。波形對比的吻合印證了本模型得到的速度波形正確性，側面證明了本仿真方法的準確性。

圖9 垂直速度波形對比Fig.9 Comparison of vertical velocity waveforms

圖10 水平速度波形對比Fig.10 Comparison of horizontal velocity waveforms

由于試驗TNT當量為500 t，根據物質云圖顯示靠近爆心的淺層測點可能存在測量失真，因此選取25.6 m遠、9.144 m深測點的速度波形與美軍手冊指導波形進行對比。經計算得到該測點的歸一化爆距并在手冊中選取波形進行擬合，二者對比結果如圖9、圖10所示。通過對比可以發現二者吻合良好，對于垂直速度波形，其先出現方向向下，幅值與向上的速度峰值相近的脈沖，隨后向上運動，之后速度緩慢下降。水平速度波形向外運動，之后速度迅速下降。波形對比的吻合印證了本模型得到的速度波形正確性，側面證明了本仿真方法的準確性。

將仿真得到的各測點正向速度峰值與美軍試驗數據使用B樣條方法進行對比，結果如圖11～圖16所示。

圖11 地面下0.457 2 m深垂直速度Fig.11 Vertical velocity at 0.457 2 m depth

圖12 地面下0.457 2 m深的水平速度Fig.12 Horizontal velocity at 0.457 2 m depth

圖13 地面下1.524 m深的垂直速度Fig.13 Vertical velocity at 1.524 m depth

圖14 地面下1.524 m深的水平速度Fig.14 Horizontal velocity at 1.524 m depth

圖15 地面下3.048 m深的垂直速度Fig.15 Vertical velocity at 3.048 m depth

圖16 地面下3.048 m深的水平速度Fig.16 Horizontal velocity at 3.048 m depth

通過曲線可以看到土體的正向速度峰值隨爆心距的增加呈現出指數衰減的趨勢，同一距離不同深度的速度峰值也隨深度而減小。通過對比可以發現仿真結果與試驗數據擬合很好，由于試驗的淺層測點表面覆土由人工回填，與周圍土體可能存在一定的非連續性，且靠近爆心的測點位于彈坑附近，容易在爆炸作用下被發掘出地面，結合其他測點的整體趨勢判斷圖11中42.7 m處試驗數據可能存在失真。

使用美軍防核設施設計手冊武器效應分冊的對地沖擊效應計算公式，根據試驗條件及當地黏土介質特性測量結果，對地下3.048 m的測點進行計算，并與上述結果進行對比，結果如圖17、圖18所示。通過對比可以看出三者整體擬合情況良好，仿真結果與試驗數據更接近，垂直速度擬合優于水平速度，且均呈現出隨爆心距的增加指數衰減的特征。其中防核設施設計手冊水平分量與仿真結果的近場計算結果相比于試驗數據均存在較大誤差，計算結果高于試驗數據，而仿真結果低于試驗數據。原因可能是由于美軍的計算公式基于大量地沖擊試驗數據進行擬合得到，由于爆炸效應的強非線性、測量數據的離散性以及擬合方式的選取等原因，使得擬合公式在強非線性段存在一定的誤差，其特征是爆炸當量越大，爆心距離越小，深度越小，誤差越大。此外有限元分析基于連續介質力學理論建立，真實試驗場地的巖土介質實際上是非連續的，介質之間存在孔隙、節理等使得介質流動性更大、連接更松散，對于接近爆心的表層巖土亦存在拋擲及氣化蒸發的可能，導致接近爆心處的巖土實際運動幅度比有限元仿真所得到的結果更大，造成了上述誤差。

圖17 地面下3.048 m深的垂直速度Fig.17 Vertical velocity at 3.048 m depth

圖18 地面下3.048 m深的水平速度Fig.18 Horizontal velocity at 3.048 m depth

5 結論

核爆炸引起的對地沖擊效應成分復雜，難以進行準確地仿真及預測。本文提出了一種綜合考慮感生沖擊波與直接沖擊波疊加效應的核爆炸對地沖擊準確描述方法，通過流體- 固體耦合技術還原了感生沖擊波，使用Shock Hugoniot絕熱關系描述的Mie-Gruneisen狀態方程以及D-P彈塑性準則來描述黏土在對地沖擊下的流動行為，并以美軍核爆炸對地沖擊試驗為例對核爆炸對地沖擊作用下土體運動特性進行了研究。得到如下結論：

1)本文所提出的仿真分析方法準確地反映了核爆炸對地沖擊作用下黏土介質的運動特性，與美軍地沖擊試驗數據以及美軍手冊計算公式結果吻合良好，可用于核爆炸對地沖擊效應計算。

2)仿真結果表明核爆炸引起的空氣沖擊波會在地面激發感生沖擊波，完整的核爆炸對地沖擊效應是感生沖擊波和直接沖擊波兩部分疊加作用的結果。在超震區淺層地表范圍內感生沖擊波起主導作用，在深層地下及超前區范圍直接沖擊波起主導作用。

3)在核爆炸打擊下地面物體會同時受到空氣沖擊波及感生沖擊波的共同作用，對于地面及淺層地表的目標，除了要經受空氣沖擊波的沖刷，還會經歷地面瞬時大范圍運動，因此在核爆炸襲擊下對于建筑強度的校核不應只局限于沖擊波超壓峰值，還應考慮其抗震性能。

4)由于大當量爆炸的強非線性、測量數據的離散性以及擬合方式的選取等原因導致在接近爆心處仿真結果、試驗數據與經驗公式存在較大誤差，因此難以準確描述爆心附近的土體運動情況。