陸浩然 孫海亮 馬強 李海濤 于麗晶 馬明輝 孫宇新

（1 北京宇航系統(tǒng)工程研究所，北京，100076；2 北方華安工業(yè)集團有限責任公司，黑龍江齊齊哈爾，161046；3 南京理工大學，江蘇省南京市，210094）

0 引言

爆炸水池是開展炸藥毀傷評估的必要平臺，也是完成戰(zhàn)斗部破碎實驗、進行水下環(huán)境爆炸毀傷能力參數(shù)測定和毀傷技術(shù)研究的重要場所[1,2]。數(shù)值仿真技術(shù)應用范圍涉及各個研究領域，仿真準確性和精確性也越來越高，科學研究正逐步采用數(shù)值仿真代替部分試驗。水下環(huán)境爆炸問題屬于較為復雜、較難進行試驗的科學研究，而且經(jīng)驗公式對水下環(huán)境爆炸相關問題分析能力有限，因此數(shù)值仿真方法顯得尤為重要。

國內(nèi)外專家對水下環(huán)境爆炸開展了大量的數(shù)值仿真：Huang H等[3]對水下爆炸問題進行了一維數(shù)值仿真，研究發(fā)現(xiàn)水的本構(gòu)粘度系數(shù)對仿真結(jié)果有影響；鐘東望等[4]將典型水下爆炸沖擊波傳播的試驗數(shù)據(jù)、經(jīng)驗公式計算數(shù)據(jù)、仿真模擬數(shù)據(jù)三者進行對比，證明了LS-DYNA在該研究條件下的準確性；張斐等[5]研究表明總藥量不變條件下，多次水下爆炸載荷作用下，焊接板與鋼板的塑性變形形貌呈類球冠形，二者的厚度減薄率從中心位置到邊界處呈先減小再增大的趨勢，但焊接板的撓度小于鋼板的撓度；黃興中等[6]仿真研究建議仿真計算時的水域尺寸應該至少大于炸藥尺寸的50倍；Amir Javad Moradloo等[7]采用流固耦合方法模擬水下爆炸條件下混凝土的非線性行為，分析得到混凝土壩受水下爆炸影響大小因素；Gaohui W等[8]研究了混凝土在爆炸沖擊下的響應和應變率效應；焦曉龍等[9]使用ALE方法模擬了多艙室內(nèi)部爆炸；Kira A等[10]將試驗測得的壓力與仿真所得的壓力進行比對，研究表明爆炸和測點的距離與測點沖擊波大小呈負指數(shù)相關性。邵旭東等[11]從理論分析、經(jīng)驗外推以及數(shù)值仿真、試驗研究等角度，系統(tǒng)總結(jié)了國內(nèi)外在航天產(chǎn)品沖擊響應分析與沖擊損傷/失效兩個方面的研究進展，并綜述分析了國內(nèi)外眾多學者提出的各種沖擊損傷/失效評估方法在航天產(chǎn)品上的實用性。陳新發(fā)等[12]總結(jié)和分析了爆炸沖擊波從產(chǎn)生到傳播的各個階段中，壓力、沖量、持續(xù)時間等影響參數(shù)的理論和經(jīng)驗公式，獲得不同比例的入射壓、反射壓等沖擊波參數(shù)，計算了Kingery和Bulmash經(jīng)驗模型參數(shù)，重構(gòu)了爆炸波模型。皮本樓等[13]提出了非接觸式電磁激勵響應板方式的爆炸分離沖擊環(huán)境模擬試驗方法，揭示了沖擊環(huán)境預示與控制的機理，為實現(xiàn)電磁激勵載荷作用下的結(jié)構(gòu)響應的可預測提供了研究基礎。目前，大部分爆炸試驗水池是基于靜力學理論基礎設計其結(jié)構(gòu)強度，但是水下爆炸作為一種高強度的爆炸沖擊波，顯然無法真實給出水池受力狀態(tài)。本文擬采用數(shù)值模擬方法，通過模擬3kg TNT藥包，在較大水池（9000mm×5000mm×5000mm）內(nèi)模擬爆炸沖擊過程，分析水池中沖擊/應力波傳播機理，揭示爆炸沖擊對水池壁面作用力規(guī)律，仿真計算與分析驗證了本文提出的模型方法的有效性。

1 計算模型

本文研究對象是一典型鋼壁水池水下爆炸，水池混凝土墻和鋼板壁使用Lagrange網(wǎng)格劃分方法，水池內(nèi)水體和空氣采用Eular網(wǎng)格離散。

1.1 幾何模型

本文的研究對象如圖1所示，其中藍色水域尺寸為8000mm×4400mm×3900mm，空水界面空氣厚度為100mm，長寬尺寸與水域相同。鋼池壁緊貼水域，厚度為40mm，混凝土外圍尺寸為9000mm×5000mm×5000mm。模型中未顯示空水界面的空氣。

圖1 典型爆炸水池整體計算模型Fig.1 Simulation model of typical explosion pool

由于本文研究沖擊波在爆炸水池中的傳播，計算時通過設置對稱，計算模型如圖 2所示。模型網(wǎng)格尺寸設置為20mm，建立四分之一模型4500mm× 2500mm×5000mm的空（Void）歐拉域，將水和空氣兩種材料對應的模型填充進歐拉域中；鋼板和混凝土使用Lagrange網(wǎng)格劃分，網(wǎng)格尺寸均為20mm。將爆炸沖擊波映射到爆炸水池中，圖2中所示炸藥尺寸已經(jīng)膨脹。

水池上方對歐拉域添加flow-out流出邊界。鋼板和混凝土通過Join連接在一起，歐拉域和拉格朗日之間使用自動流固耦合處理計算。在混凝土外圍三個面（圖2示右、后、下面）設置無反射邊界，視為無限混凝土域。

圖2 爆炸水池計算模型主視圖及網(wǎng)格劃分Fig.2 Front view and grid of typical explosion pool simulation model

1.2 材料模型

本文中涉及材料包括水、TNT、混凝土、鋼、空氣，均直接使用軟件材料庫中的材料本構(gòu)方程和相應材料參數(shù)。

1.2.1 混凝土本構(gòu)模型

混凝土本構(gòu)模型包括失效面、彈性極限面、應變硬化、殘余失效面和損傷模型五個部分。RHT的本構(gòu)模型采用的Polynomial狀態(tài)方程為

其中，e為比內(nèi)能，A1為體積模量，A2、A3、B0、土，本構(gòu)相關參數(shù)如表 1所示（其中除了參數(shù)B1、T1、T2為材料參數(shù)，本文選用35MPa混凝B0、B1，其他參數(shù)的單位均為kPa）。其中失效面的函數(shù)為

表1 混凝土材料參數(shù)Table1 Material parameter of concrete

式中，P為所受壓力大小，P*是通過fc歸一后的壓力值，是P*關于 的函數(shù)，fc為單軸壓縮強度，ft為拉伸強度，A為失效面的相關常數(shù)，θ是加載角度，N是失效面指數(shù)，ε˙是材料的應變率，D和α分別為材料受壓縮和拉伸時的應變率常數(shù)。

其中，felastic是彈性強度和失效面強度的比值，可以根據(jù)壓縮強度cf和拉伸強度tf兩參數(shù)來計算得出，F(xiàn)CAP(P)為常數(shù)，由具體受壓力值大小來計算確定。

應變硬化屈服面參數(shù)由Y*上述失效面failY和彈性極限面Yelastic計算得出

其中G代表各狀態(tài)下的剪切模量。殘余失效面通過殘余失效面常數(shù)B和指數(shù)M確定

損傷模型相關計算如下

其中，Di為損傷經(jīng)驗常數(shù)，為最小失效應變，同時可以進一步得到損傷后的失效面和剪切模量Gfractured為

1.2.2 鋼板材料本構(gòu)模型

鋼板使用4340鋼，鋼板的屈服應力yσ可以表示為

式中，εp為等效塑性應變；ε*是塑性應變率；A是屈服應力；B是硬化常數(shù)；n是硬化指數(shù)；C是鋼材應變率常數(shù)；m是熱軟化指數(shù)；Th通過室溫tT和熔化溫度mT確定

Linear狀態(tài)方程如下

其中，K為體積模量。綜上，鋼材料模型的部分參數(shù)如表2所示（除了指數(shù)n為常數(shù)，其他參數(shù)單位均為kPa）。

表2 鋼板材料參數(shù)Table2 Material parameter of steel plate

1.2.3 空氣的狀態(tài)方程

空氣的狀態(tài)方程為以下理想狀態(tài)方程

其中，P表示壓力，γ表示絕熱指數(shù)，0e是比內(nèi)能。具體采用的參數(shù)如表3所示。

表3 空氣的材料參數(shù)Table3 Material parameter of air

2 仿真計算與分析

在適用范圍內(nèi)經(jīng)驗公式所得結(jié)果高度可信，本文將以由水下自由場爆炸沖擊波經(jīng)驗公式為基準，校驗本數(shù)值模擬可信性。對于球形的TNT裝藥，其沖擊波壓力峰值為

式中Pm表示沖擊波壓力峰值，單位MPa；W是炸藥質(zhì)量，單位kg ；R是距爆心距離，單位是m；0R則是裝藥的初始半徑，單位也是m。

由于本文研究對象尺寸較大，采用Remap方法將一維楔形爆炸模型映射到三維球體模型中進行計算。當炸藥位于水下的深度與裝藥的半徑之比大于10～20，且炸藥在水下的深度與裝藥的半徑比值超過5～10時，該藥量炸藥在水下爆炸時產(chǎn)生的沖擊波壓力峰值可以視為不受外界影響，包括空水界面和其他介質(zhì)反射等，而本文所討論的裝藥位置與裝藥半徑比值大于上述值，因此炸藥爆炸短時間內(nèi)可以視為自由場水下爆炸。炸藥采用球形TNT炸藥，炸藥當量為3kg。

2.1 3kg TNT藥包水下起爆后第一階段仿真

以藥包起爆開始后0.5ms內(nèi)為第一階段。3kg TNT水下爆炸模型如圖 3所示，因為炸藥起爆位置距離空水界面2000mm，所以楔形模型中水域長2000mm，同樣添加無反射邊界視為無限水域，炸藥半徑為76.02mm。在距爆心700mm到1100mm處，每間隔100mm設置一個，由于爆距范圍為700～1100mm，裝藥初始半徑76.02mm，根據(jù)經(jīng)驗公式（16）和（17）進行計算，仿真所得結(jié)果與經(jīng)驗公式計算結(jié)果如表4所示，數(shù)值仿真所得值與經(jīng)驗公式所得偏差不大于2.4%，表明計算可信。

表4 3kg TNT水下起爆壓力峰值Table4 Underwater explosion peak pressure of 3kg TNT

圖3 3kg TNT水下爆炸一維模型Fig.3 One dimension model of 3kg TNT underwater explosion

由于本文研究的是有限區(qū)域爆炸，無限水域中沖擊波計算距離不能太長：楔形模型計算到0.5ms時，輸出文件并寫入爆炸水池計算模型，0.5ms時刻的壓力波分布如圖 4所示，此時沖擊波還遠未到達2000mm處 。

圖4 3kgTNT水下爆炸0.5ms時刻壓力云圖Fig.4 Pressure nephogram of 3kg TNT underwater explosion at 0.5ms

2.2 3kg TNT藥包水下起爆后第二階段仿真

由于2.1節(jié)使用的映射是0.5ms之后的爆炸沖擊波數(shù)據(jù)，因此下文所提到的時間均默認加上映射前的0.5ms，所有相關曲線圖的時間坐標也以0.5ms為起始坐標。

從水中壓力云圖可以看出炸藥從起爆到0.9ms時刻，基本屬于無限水域爆炸情形，爆炸沖擊波呈球形向外傳播；沖擊波在1.1ms左右傳播到空水交界面和水池底部鋼板，因為空氣介質(zhì)的密度和聲速都遠小于水介質(zhì)，所以空氣介質(zhì)的波阻抗遠遠小于水介質(zhì)的波阻抗，那么此時主要是產(chǎn)生與入射波方向相反的反射波，而透射波數(shù)值較小，這從圖中可以觀察到，沖擊波傳播到空水界面時發(fā)生反射產(chǎn)生方向相反的反射波，而由于壓力顯示區(qū)間問題，導致較小的透射波從云圖中不易觀察到；而水池底部的交界面處，水介質(zhì)的波阻抗小于鋼介質(zhì)的波阻抗，可以看到鋼壁與水交界處的確出現(xiàn)了紅色高壓區(qū)，而繼續(xù)往下傳播的是透射波，此情況下，透射波大小大于入射波。圖 5（d）的三處紅色高壓區(qū)也是因為上述波阻抗不同的原因形成的，但是從池底處依然可以發(fā)現(xiàn)出現(xiàn)了與入射波方向相反的反射波傳播軌跡，這是因為水下爆炸是很多方向沖擊波的疊加。圖 5（e）-（g）主要是沖擊波傳播到較近鋼壁和池底的反射與向較遠鋼壁的類球形傳播過程，并且圖示左方與下方的沖擊波相遇形成高壓力區(qū)，而后又沿著各自方向繼續(xù)傳播。從圖 5（h）開始，主要是反射波之間的碰撞疊加，隨后繼續(xù)傳播又反射，循環(huán)進行，直至壓力最終消耗殆盡。

圖5 3kgTNT水池中爆炸壓力云圖Fig.5 Stress nephogram of 3kg TNT underwater explosion at test pool

圖6是將水域隱藏后的池壁等效應力云圖，目的是觀察鋼壁在炸藥爆炸短時間內(nèi)的應力變化。炸藥起爆后1.1ms，離爆心最近的池底最先開始產(chǎn)生應變；隨著爆炸沖擊波呈球形向四周傳播，鋼壁上對應的應力峰值也呈圓形向外傳播，如圖6（b）、（c）所示；圖6（d）和（e）呈現(xiàn)了沖擊波在兩鋼壁夾角處的疊加同樣使拐角交界處的應力產(chǎn)生疊加，產(chǎn)生了應力集中，如圖6（e）中紅色區(qū)域所示；圖6（f）的應力云圖變化，說明沖擊波的反射同樣體現(xiàn)在了鋼壁的應力變化趨勢上；從最后兩張圖中可以看出，后部分時間鋼壁應力體現(xiàn)在較遠壁處，且應力最大的是三鋼板交界處，因為在這些地方會發(fā)生應力集中。另外，從圖6（g）和（e）中可以發(fā)現(xiàn)鋼壁距離池底越近處應力越大，這主要是沖擊波在空水界面發(fā)生卸載導致的。

圖6 3kgTNT水池中爆炸池壁等效應力云圖Fig.6 Equivalent stress nephogram at pool wall of 3kg TNT underwater explosion

2.3 各觀測點應力分析

本文研究對象所設置觀測點如圖 7所示：其中1～5號高斯點為鋼壁靠水表面觀測點，1號點與裝藥中心在同一水平線，1、2、3點處于同一水平線，間隔1000mm。1、4、5號點在同一豎直線上，同樣也間隔1000mm。6～8號點為水域中的定觀測點，三個點與爆炸中心在同一水平線上，與爆心距離2000、3000和4000mm；9-12號觀測點同樣在鋼壁上隨動，水平位置和間距與1～3號點一致。

圖7 計算模型高斯點設置Fig.7 Gauss points of simulation model

2.3.1 水中壓力傳播分析

沖擊波在水池中傳播的大致規(guī)律，6～8號觀測點對應的時間—壓力擬合曲線如圖8（a）所示。圖（a）和（b）的不同點在于，（a）是在爆炸水池中傳播、（b）在無限水域中傳播，圖 8（a）的6、7、8點對應圖8（b）的1、3、5號測點。因為爆炸沖擊波傳播到池壁時會產(chǎn)生各個方向的反射波，且會產(chǎn)生多次反射，因此圖8（a）情況的后面時間段的沖擊波壓力變化曲線錯綜復雜。但是相同的是，6～8點在3.5ms左右第一段波基本也趨于0，與圖 8（b）相同，且6號與1號、7號與3號點第一次到達波峰的時間與大小都幾乎相同，但是8號與5號觀測點雖然到達波峰時間一樣，但是8號觀測點的壓力大小明顯大于5號點，甚至超過了爆距更近的7號點處沖擊波壓力，從圖8（g）中可以發(fā)現(xiàn)，是因為2.43ms左右在8號測點產(chǎn)生了反射波，沖擊波疊加，使得該峰值大小“異常”。

圖8 （a）6-8號高斯點壓力變化擬合曲線（b）無限水域情況壓力曲線Fig.8 (a) Pressure fitting curve at Gauss point 6～8 (b) Pressure of infinite fluid

2.3.2 池壁各點處等效應力曲線分析

圖9是離爆心較遠處的鋼壁上觀測點隨時間變化的等效應力（MIS.STRESS）曲線，由于3、5點處應力集中，應力峰值較其他觀測點較大，為了便于觀察單獨擬合，下文對1-5號測點的應力變化曲線也分開擬合。

圖9 無防護措施下1-5號觀測點處等效應力曲線Fig.9 Equivalent stress nephogram at point 1～5 without protective measure

1、2和4 號觀測點第一個波峰到達時間相差不大，數(shù)值上1處最小是因為他的峰值到達最早，此瞬時的沖擊波對鋼壁的應力影響還未傳遞到2和4，所以沖擊波隨著時間推移傳播到2和4點 時，會產(chǎn)生一個疊加后的波峰，自然就比前一時刻1處的應力峰值高。而對于每個測點處的最大應力來說，4＞1＞2，可以推斷出爆炸水池爆炸時，同一水平線上，爆距越近，鋼壁的最大應力越大， 而豎直方向上，因為空水界面對沖擊波的卸載，以及池底的反射作用，反而距池底較近但離爆心較遠的4號點處的最大應力更大。5號觀測點第一個波峰略高于3號觀測點，兩點處的最高波峰可以視為有兩處，圖中均已標注出來，且總體壓力峰值遠大于3號位置，是因為3號位置距離水面較5號近，許多沖擊波傳播到空水界面發(fā)生卸載，同時沖擊波傳到池底會波及到5號測點。

根據(jù)3.2節(jié)對云圖的初步分析，距離爆心較近鋼壁上的應力到6ms之前就已經(jīng)達到峰值了，9-12號觀測點等效應力擬合曲線如圖 10所示，結(jié)合圖6中不同時刻鋼壁應力云圖，各位置第一次波峰到達時間是匹配的。同樣地，爆距越大，第一個波峰應力大小越大（圖示中11和12號觀測點的第二個波峰視為第一個波峰，可能是因為應力疊加的時間相差較大出現(xiàn)兩個波峰），這是因為沖擊波初步傳播過程，隨著時間的增加，越晚到達的位置受到應力疊加影響越大，從而使得應力越高；第二個標注出的峰值是較近池壁與池底沖擊波疊加后傳播到各觀測點產(chǎn)生的，如圖 5的最后三張圖所示，第二個壓力峰值即是各點處的最大應力，對于各點處最大應力的大小而言，觀測點位置離爆心越近，鋼壁的應力峰值也就越大。

圖10 無防護措施下9-12號觀測點等效應力曲線Fig.10 Equivalent stress nephogram at point 9～12 without protective measure

3 結(jié)論

本文開展了爆炸沖擊對水池壁面作用力數(shù)值仿真，揭示了爆炸沖擊波在水池中會發(fā)生多次反射機制，得出了水池壁面受力規(guī)律：各池壁交界處，由于應力集中，應力明顯大于鋼壁中心區(qū)域；不考慮應力集中區(qū)，同一水平線上，距爆心越近，應力越大，而同一豎直線上，因為水池上的空水界面對沖擊波的卸載作用，加上沖擊波傳播到池底發(fā)生反射后繼續(xù)向側(cè)面池壁傳播，因此距爆心越近但是離池底越遠處的應力反而越??；較遠池壁整體X方向最大位移是0.65mm，最大的正向加速度是 1.01e3m/s2，最大反向加速度是0.64e3m/s2。