鋼箱內(nèi)部爆炸破壞的SPH 數(shù)值模擬*

強(qiáng)洪夫，孫新亞，王 廣，黃拳章

（火箭軍工程大學(xué)，陜西 西安 710025）

橋梁是公路、鐵路等交通網(wǎng)的咽喉。為了減小恒載，對(duì)斜拉橋、懸索橋等大跨度橋廣泛采用鋼箱梁作其結(jié)構(gòu)支撐的主梁。隨著恐怖襲擊的不斷演化，戰(zhàn)略性橋梁成為恐怖主義爆炸破壞的重點(diǎn)目標(biāo)，因此，有必要對(duì)鋼箱結(jié)構(gòu)內(nèi)部爆炸過(guò)程中鋼箱的變形規(guī)律進(jìn)行研究[1-4]。Tang 等[5]運(yùn)用LS-DYNA 軟件研究了汽車炸彈爆炸沖擊作用下橋塔、橋墩、橋面結(jié)構(gòu)的局部破壞模式；姚術(shù)健等[2,6]通過(guò)實(shí)驗(yàn)，研究了鋼箱內(nèi)部爆炸過(guò)程中內(nèi)壁失效變形的模式, 并利用ANSYS 軟件分析了箱體的損傷特性。鋼箱內(nèi)部爆炸相比外部爆炸而言，爆炸沖擊波在約束空間內(nèi)多次反射、疊加、匯聚，該物理過(guò)程相當(dāng)復(fù)雜，目前的研究大多局限于爆炸過(guò)程中鋼箱變形的宏觀現(xiàn)象。光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics,SPH)作為一種無(wú)網(wǎng)格的粒子方法，因其在計(jì)算過(guò)程中可避免網(wǎng)格重分和算法耦合，被廣泛用于研究爆炸沖擊問(wèn)題。自Monaghan 等[7]修正的人工黏度解決了SPH 計(jì)算結(jié)果的非物理震蕩問(wèn)題后，強(qiáng)洪夫等[8]、Qiang 等[9]和Liu 等[10]利用該方法對(duì)爆炸沖擊波沖擊鋼板的過(guò)程成功地進(jìn)行了數(shù)值模擬。本文中，擬利用自編的SPH 程序?qū)︿撓鋬?nèi)部爆炸過(guò)程進(jìn)行數(shù)值模擬，并與實(shí)驗(yàn)對(duì)比，驗(yàn)證SPH 算法在模擬鋼箱內(nèi)部爆炸問(wèn)題上的有效性；分析爆炸過(guò)程中鋼箱表面壓力、應(yīng)力波的傳播規(guī)律以及特定點(diǎn)處速度變化趨勢(shì)；得到鋼箱表面的von Mises 應(yīng)力分布；進(jìn)一步對(duì)鋼箱內(nèi)不同位置爆炸的過(guò)程進(jìn)行數(shù)值模擬，得到不同工況下鋼箱的失效變形模式，以期通過(guò)研究為鋼箱類結(jié)構(gòu)的爆炸毀傷提供一種評(píng)估方法。

1 模型與算法

1.1 控制方程和狀態(tài)方程

1.1.1 控制方程

鋼箱的內(nèi)部爆炸是絕熱過(guò)程，外力相對(duì)于內(nèi)部驅(qū)動(dòng)力可忽略不計(jì)，即可不考慮重力的影響，在流體控制方程中忽略黏性項(xiàng)。因此，其控制方程組表述為：

式中：ρ、v、p、e、x 分別為密度、速度、壓強(qiáng)、能量、位移。

1.1.2 狀態(tài)方程

數(shù)值模擬過(guò)程中，TNT 炸藥的爆轟氣體壓力通過(guò)標(biāo)準(zhǔn)的JWL 狀態(tài)方程[11]獲得：

對(duì)于鋼箱材料，選用Q235 鋼。爆炸過(guò)程中鋼箱的各向同性壓力p 利用Mie-Grüneisen 狀態(tài)方程[12]求得:

1.2 本構(gòu)方程

鋼箱材料模型采用B?rvik 等[13]、Nilson[14]提出的修正Johnson-Cook 本構(gòu)模型，該模型通過(guò)流動(dòng)應(yīng)力和失效應(yīng)變描述金屬材料的大變形、高應(yīng)變率等動(dòng)態(tài)力學(xué)特性。其流動(dòng)應(yīng)力的表達(dá)式為：

1.3 控制方程的SPH 離散

SPH 算法是模擬流體流動(dòng)的一種拉格朗日型粒子方法，通過(guò)使用一系列任意分布的粒子來(lái)求解具有各種邊界條件的積分方程或偏微分方程。SPH 方法通常通過(guò)核函數(shù)插值實(shí)現(xiàn)場(chǎng)變量或場(chǎng)變量梯度的插值，通過(guò)粒子近似實(shí)現(xiàn)對(duì)核函數(shù)估計(jì)積分表達(dá)式的粒子離散[9]。

本文中利用SPH 方法模擬鋼箱內(nèi)部爆炸過(guò)程屬于該方法在具有大變形和高度非均勻的動(dòng)力學(xué)極端情形下的應(yīng)用。為了解決爆炸域內(nèi)求解結(jié)果的非物理振蕩，且將沖擊面內(nèi)的動(dòng)能以熱能的形式耗散，引入人工黏性[7]，同時(shí)為了更好地處理爆炸過(guò)程中鋼箱大變形、大扭曲等引起的密度和光滑長(zhǎng)度變化劇烈的問(wèn)題，采用Qiang 等[9]提出的完全變光滑長(zhǎng)度SPH 方法對(duì)鋼箱內(nèi)部爆炸過(guò)程進(jìn)行數(shù)值模擬，結(jié)合后的方程組如下：

式中： fi為 修正系數(shù)； Wij=W(xi-xj,h)為核函數(shù)，一般選用三次樣條函數(shù)；Πij為人工黏度；N 為求解區(qū)域內(nèi)粒子總數(shù)；α、β 表示空間坐標(biāo)軸方向； h 為光滑長(zhǎng)度，在連續(xù)性方程中將其對(duì)時(shí)間求導(dǎo)，得：

式中：d 為空間維數(shù)。

2 鋼箱內(nèi)部爆炸過(guò)程的數(shù)值模擬

2.1 數(shù)值模型

為了證明SPH 方法在數(shù)值模擬鋼箱結(jié)構(gòu)內(nèi)部爆炸過(guò)程時(shí)的有效性，本文中設(shè)定了與姚術(shù)健等[2]的實(shí)驗(yàn)一致的工況條件，如圖1(a)所示，鋼箱模型為立方體形狀，邊長(zhǎng)為600 mm，且采用了120 mm 的延伸結(jié)構(gòu)，Q235 鋼板厚度為4 mm；TNT 炸藥為立方體形狀，邊長(zhǎng)為40 mm。具體粒子配置如圖1（b）所示，粒子間距為4 mm，炸藥粒子數(shù)為1 331，鋼箱粒子數(shù)為265 600。

圖 1 鋼箱和TNT 炸藥試件模型圖Fig.1 Model diagrams for the steel box and TNT dynamite specimen

為保證箱體各壁面受力均衡，TNT 炸藥位于鋼箱內(nèi)部正中心，且各邊與箱體對(duì)應(yīng)邊平行，如圖2 所示。計(jì)算中，TNT 炸藥密度為1 500 kg/m3, 采用點(diǎn)起爆方式，起爆點(diǎn)設(shè)置在炸藥上表面中心點(diǎn)（-73.9，5.6，50.7），爆速為6 930 m/s，初始內(nèi)能為4 290 kJ/kg，狀態(tài)方程選用JWL 狀態(tài)方程；鋼箱密度為7 800 kg/m3,使用Johnson-Cook(J-C)本構(gòu)模型，具體材料參數(shù)A=369 MPa, B=684 MPa, n=0.73, C=0.083, m=1.7,Tm=775 K,=0.1 s-1, CV=875 J/(kg·K), D1=0.13, D2=0.13, D3=-1.5, D4=0.011, D5=0。在組成鋼箱的各板之間施加焊接錨固作用，由文獻(xiàn)[15]可知焊縫金屬?gòu)?qiáng)度不應(yīng)低于母材。為簡(jiǎn)化計(jì)算，忽略焊縫附近區(qū)域母板材料性能的變化，將焊縫與母材視為同一結(jié)構(gòu)，即認(rèn)為鋼箱是一體的。

圖 2 鋼箱內(nèi)部爆炸算例模型Fig.2 The case model for the explosion inside the steel box

2.2 鋼箱表面變形情況

數(shù)值模擬中，核函數(shù)選用三次樣條函數(shù)，光滑長(zhǎng)度取1.5 倍粒子間距，時(shí)間積分采用蛙跳格式，時(shí)間步長(zhǎng)取0.1 μs。鋼箱內(nèi)部爆炸的數(shù)值模擬中，鋼箱表面的變形情況如圖3 所示。起始時(shí)刻，鋼箱表面無(wú)變形，1.97 ms 時(shí)中心處產(chǎn)生明顯的向外鼓起，并向四周傳播變形；3.34 ms 時(shí)變形首次到達(dá)表面邊沿處，且表面角隅處發(fā)生變形；12.03 ms時(shí)，表面變形達(dá)到最大，并保持不變。

圖4 給出了鋼箱內(nèi)部爆炸實(shí)驗(yàn)和數(shù)值模擬的對(duì)比，圖中A 處鋼板發(fā)生卷邊破壞，B 處鋼板向外鼓起，C 處鋼板發(fā)生扭曲變形，通過(guò)對(duì)比發(fā)現(xiàn)，數(shù)值模擬和實(shí)驗(yàn)所得鋼箱壁面失效變形模式吻合較好。圖5 給出了鋼箱內(nèi)部爆炸過(guò)程中箱體壁面撓度的變化情況。從圖5 可以看出，85 μs 時(shí)鋼箱壁面開始變形，撓度值隨時(shí)間逐漸增大，12.03 ms 后基本保持不變。數(shù)值模擬得到的鋼箱壁面最終撓度值為31.80 mm，而實(shí)驗(yàn)值為32.93 mm，實(shí)驗(yàn)和數(shù)值模擬的相對(duì)誤差為3.5%。產(chǎn)生這一誤差的主要原因可能是：(1)炸藥理論設(shè)計(jì)的質(zhì)量為96 g，由于密度不均勻等原因?qū)嶒?yàn)中炸藥的真實(shí)質(zhì)量卻為98.443 g；(2)本文數(shù)值模擬過(guò)程中Q235 鋼材料密度是均勻的，沒(méi)有考慮制作中的工藝誤差。綜上所述，可得SPH方法能夠數(shù)值模擬鋼箱內(nèi)部爆炸過(guò)程中鋼箱壁面的變形情況。

圖 3 鋼箱內(nèi)部爆炸過(guò)程中不同時(shí)刻箱體內(nèi)壁的變形情況Fig.3 Its inner wall deformation at different times for the steel box subjected to inner explosion loading

圖 4 鋼箱壁面破壞的數(shù)值模擬和實(shí)驗(yàn)對(duì)比Fig.4 Comparison of steel box wall damage between numerical simulation and experiment

圖 5 鋼箱壁面撓度變化趨勢(shì)Fig.5 Deflection in the center of the steel box wall varying with time

2.3 鋼箱表面壓力和von Mises 應(yīng)力分布

圖 6 鋼箱內(nèi)部爆炸過(guò)程中，壁面von Mises 應(yīng)力分布Fig.6 Von Mises stress inside the steel wall during the explosion in the steel box

圖 7 鋼箱內(nèi)部爆炸過(guò)程中，壁面的壓力分布Fig.7 Pressure distribution of the wall during the explosion in the steel box

裝藥在完全封閉的鋼箱內(nèi)部爆炸時(shí)，爆炸沖擊波在鋼箱四周壁面上發(fā)生多次反射，沖擊波強(qiáng)度逐漸降低。圖6 給出了鋼箱內(nèi)部爆炸時(shí)，鋼箱表面的von Mises 應(yīng)力分布；圖7 給出了鋼箱內(nèi)部爆炸過(guò)程中，鋼箱表面壓力分布情況。從圖6 和圖7 可得：鋼箱內(nèi)部炸藥爆炸后，很快產(chǎn)生爆炸沖擊波并迅速（85.0 μs）到達(dá)鋼箱表面，爆心投影點(diǎn)附近出現(xiàn)壓強(qiáng)和應(yīng)力最大值；隨后沖擊波沿著鋼箱的內(nèi)表面向板與板的連接處（角隅位置）運(yùn)動(dòng)，同時(shí)爆心投影點(diǎn)處超壓峰值逐漸降低，436.1 μs 時(shí)各板向角隅運(yùn)動(dòng)的沖擊波在角隅處匯聚，形成較強(qiáng)的沖擊波壓力，并出現(xiàn)應(yīng)力集中和應(yīng)力增強(qiáng)現(xiàn)象；隨后各板反射的沖擊波與角隅處匯聚加強(qiáng)了的沖擊波向箱室中心匯聚，使箱室中心形成較高超壓，箱室中心超壓又會(huì)向四周傳播，此過(guò)程與起爆后傳播過(guò)程類似，鋼箱表面承受多次類似的沖擊載荷作用后，鋼箱表面壓強(qiáng)和von Mises 應(yīng)力逐漸消失，且箱內(nèi)維持較長(zhǎng)時(shí)間的余壓。

2.4 鋼箱表面中心點(diǎn)處速度和壓強(qiáng)

圖8 給出了鋼箱內(nèi)部爆炸過(guò)程中，鋼箱表面中心點(diǎn)處速度的變化趨勢(shì)；圖9 給出了鋼箱內(nèi)部爆炸過(guò)程中，鋼箱表面中心點(diǎn)處壓強(qiáng)的變化趨勢(shì)。從圖8 可以看出：101.2 μs 時(shí)（如點(diǎn)A 所示），表面中心點(diǎn)的速度跳躍式增加至22.3 m/s，其后速度值又迅速減小，最終減小至0，且速度值在B 區(qū)域內(nèi)出現(xiàn)上下微小浮動(dòng)。從圖9 可以看出：85.0 μs 時(shí)，鋼箱表面中心點(diǎn)處的壓強(qiáng)迅速增大，達(dá)到峰值487 MPa； 131.6 μs 時(shí)表面壓強(qiáng)達(dá)到第二次峰值394 MPa；隨后壓強(qiáng)值上下擾動(dòng)并逐漸減弱，出現(xiàn)擾動(dòng)的主要原因是爆炸沖擊波在鋼箱表面上的來(lái)回反射；2 000.0 μs 時(shí)壓強(qiáng)值逐漸減弱直至消失。由文獻(xiàn)[16]可得：壓力和速度時(shí)程曲線形狀與實(shí)驗(yàn)及經(jīng)典沖擊波理論相吻合，即初次峰值較大，然后發(fā)生多次反射，速度和壓強(qiáng)逐漸減弱。

圖 8 鋼箱壁面中心點(diǎn)處的速度變化趨勢(shì)Fig.8 Velocity varying with time at the wall center of the steel box

圖 9 鋼箱壁面中心點(diǎn)處的壓強(qiáng)變化趨勢(shì)Fig.9 Pressure varying with time at the wall center of the steel box

3 鋼箱內(nèi)部不同位置爆炸失效模式

為了進(jìn)一步探討同一爆炸當(dāng)量下，鋼箱內(nèi)部不同位置爆炸時(shí)，鋼箱結(jié)構(gòu)的破壞程度。本文設(shè)定了以下3 種計(jì)算工況（如圖10 所示）：工況1，炸藥位于鋼箱底部正中心；工況2，炸藥在底部中心正上方60 mm 處；工況3，炸藥位于底板與側(cè)板之間的角隅處，距離各板60 mm。

圖11 給出了3 種工況下鋼箱表面的失效變形模式。圖11(a)為工況1 鋼箱的變形情況，鋼箱底板局部特征破壞明顯，產(chǎn)生圓形的沖切破壞，斷裂口處相對(duì)比較平滑，且底板產(chǎn)生裂紋（A）；同時(shí)底板與側(cè)板的角隅處產(chǎn)生了長(zhǎng)條棗核形狀的裂口（B）。圖11(b)為工況2 鋼箱的變形情況，鋼箱底板劇烈破壞，產(chǎn)生不規(guī)則的撕裂口，斷裂口呈鋸齒狀，且底板沒(méi)有產(chǎn)生比較規(guī)則的裂紋（A）；同時(shí)底板與側(cè)板的角隅處也產(chǎn)生了長(zhǎng)條狀棗核形狀的裂口，裂口周圍向外凸起（B）。圖11(c)為工況3 鋼箱的變形情況，距離炸藥較近的角隅處發(fā)生嚴(yán)重破壞，斷裂口周圍鋼板發(fā)生扭曲變形（A）；同時(shí)發(fā)現(xiàn)，在遠(yuǎn)離炸藥的3 個(gè)鋼板的角隅處發(fā)生了撕裂破壞（B）。結(jié)合各工況下鋼箱失效變形的情況，可得炸藥在鋼箱角隅處爆炸時(shí)破壞最為嚴(yán)重。

圖 10 鋼箱中炸藥位置示意圖(a=60 mm)Fig.10 The explosive location diagram in the steel box(a=60 mm)

圖 11 3 種工況下，鋼箱表面局部破壞示意圖Fig.11 Surface partial destruction diagrams of the steel box under three different conditions

4 結(jié) 論

采用光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics, SPH)方法對(duì)鋼箱內(nèi)部爆炸的過(guò)程進(jìn)行數(shù)值模擬，研究了爆炸過(guò)程中鋼箱表面的變形，分析了相同當(dāng)量炸藥在鋼箱內(nèi)部不同位置爆炸時(shí)鋼箱表面的失效模式。

（1）通過(guò)比較實(shí)驗(yàn)與SPH 方法數(shù)值模擬獲得的結(jié)果，數(shù)值模擬與實(shí)驗(yàn)所得的失效變形模式吻合較好，實(shí)驗(yàn)和數(shù)值模擬得到的鋼箱表面撓度值相對(duì)誤差為3.5%，驗(yàn)證了本文SPH 方法在模擬鋼箱內(nèi)部爆炸問(wèn)題上的有效性；

（2）采用SPH 方法數(shù)值模擬鋼箱內(nèi)部爆炸的過(guò)程，得出與實(shí)驗(yàn)一致的變形結(jié)果。同時(shí)也發(fā)現(xiàn)，在鋼箱內(nèi)部爆炸過(guò)程中，鋼箱角隅處出現(xiàn)超壓和應(yīng)力集中現(xiàn)象，這與實(shí)驗(yàn)吻合較好。

（3）進(jìn)一步研究了相同當(dāng)量炸藥在鋼箱內(nèi)部不同位置爆炸時(shí)鋼箱失效破壞情況，發(fā)現(xiàn)炸藥在鋼箱角隅處爆炸時(shí)鋼箱破壞最嚴(yán)重。

以上研究結(jié)果可為工事結(jié)構(gòu)的防護(hù)設(shè)計(jì)提供參考。