徐自強，王麗娟，朱 潔，李 闖，宋長杰，時小凡

(1.中煤科工集團 沈陽設(shè)計研究院，沈陽 110168；2.沈陽建筑大學，沈陽 110168)

爆炸的成因直接決定爆炸的形式和沖擊作用，因此，對爆炸的內(nèi)部機理和發(fā)展過程進行研究顯得尤為重要。在研究遠場爆炸時，可以假設(shè)荷載為自由場中的無限域，進一步利用強沖擊荷載作用下的基本方程進行求解，國內(nèi)外研究人員相繼開展大量研究試驗，得到較多成果[1]；但研究近場爆炸時，這種假設(shè)難以成立，由于其邊界情況復雜，導致理論求解極為困難[2,3]，此外，由于近場爆炸破壞力強，約束條件苛刻，受到技術(shù)和經(jīng)濟等因素的制約，導致相關(guān)試驗研究較少[4,5]。為降低經(jīng)費，并且達到研究目的，采用數(shù)值模擬是現(xiàn)在較好的研究手段之一[6,7]。

近場爆炸數(shù)值模擬大多采用任意拉格朗日-歐拉算法(ALE)和流固耦合算法。ALE既能保持合理的網(wǎng)格形狀，又能準確地描述物體運動，但同時也使計算更加復雜，對方程求解造成困難[8,9]；流固耦合法不需要耦合面上的流體單元和結(jié)構(gòu)單元一一對應(yīng)，能夠大幅減少工作難度，但會造成精度偏低[10]。近年來，逐漸出現(xiàn)新型數(shù)值模擬方法，如爆炸粒子算法(Particle blast method，PBM)。PBM基于分子運動理論，能夠更好地表征高溫下的氣體行為，在一定條件下可以更準確的模擬爆炸產(chǎn)物、空氣和結(jié)構(gòu)之間的相互作用[11,12]。

為探究PBM算法在近場爆炸問題中的適用性，仿真模擬以T Brvik的AL-6XN不銹鋼板近場爆炸試驗為依據(jù)，對比了ALE、S-ALE和PBM方法的計算精確度和運算效率；采用PBM算法對不同爆心距的鋼板近場爆炸工況進行仿真計算，分析了粒子總數(shù)和粒子個數(shù)比對PBM算法計算精度的影響；通過對150 mm爆心距、球狀C4炸藥近場爆炸場景進行仿真，分析了不同算法運算時效占比規(guī)律。

1 PBM方法理論

PBM是在微粒子法(Corpuscular particle method，CPM)的基礎(chǔ)上發(fā)展而來的。CPM是一種基于動力學分子理論(Kinetic molecular theory，KMT)的粗粒度多尺寸粒子方法，多用于氣體動力學模擬。其中，KMT理論是研究氣體粒子及其微觀相互作用進而得出宏觀氣體規(guī)律的理論，它源于Daniel Bernoulli的活塞理論，后因麥克斯韋方程的建立，得以廣泛運用[13]。

KMT理論基于以下假設(shè)：①氣體由大量分子構(gòu)成，這些分子是剛性球狀的物體，始終處于隨機運動狀態(tài)；②分子尺寸遠小于分子間的平均距離；③分子遵循牛頓運動定律；④分子之間及分子與結(jié)構(gòu)間僅存在彈性相互作用；⑤分子之間不存在排斥力和吸引力。

KMT理論可以準確描述理想氣體的各項特性，但計算量巨大，難以在實際模擬中得到運用。為彌補這一缺陷，降低運算量，Olovsson博士提出了將大量分子聚集成一個粒子的方法，即CPM法。但CPM法的下述特征使其無法直接用于爆炸問題模擬：①理想氣體定律無法得到絕熱膨脹條件下爆炸進程中的急劇壓降；②系統(tǒng)處于熱平衡的假設(shè)無法滿足高溫高壓氣體動力學的要求。

PBM改進了CPM，使其可以模擬高溫高壓下的真實氣體規(guī)律。PBM算法引入了協(xié)體積效應(yīng)參數(shù)，能夠更好地反映氣體在極端壓力下的行為。PBM將許多分子合并為一個微粒，使系統(tǒng)的計算對象數(shù)量降低了十幾個數(shù)量級，見圖1。結(jié)合增加的計算臨界時間步長，這種方法使得模擬宏觀系統(tǒng)和長時間問題成為可能。

圖1 爆炸粒子法(PBM)(出自文獻[13])Fig. 1 Particle Blast Method(PBM)

2 近距離爆炸試驗仿真模擬

2.1 試驗場景及分析模型

圖2 試驗裝置(出自文獻[14])Fig. 2 Test device

圖3 試驗簡圖(作者自繪)Fig. 3 Test diagram

表1 AL-6XN鋼板近場爆炸試驗工況Table 1 Near field explosion test conditions of AL-6XN steel plate

仿真計算采用LS-Dyna有限元軟件，分別采用傳統(tǒng)ALE法、改進的S-ALE法及PBM法開展仿真模擬。根據(jù)試驗場景建立1∶1仿真模型，如圖4所示。AL-6XN鋼板采用全積分殼體單元，網(wǎng)格單元尺寸為7 mm，考慮應(yīng)變率效應(yīng)，材料選用關(guān)鍵字MAT_MODIFIED_JOHNSON_COOK來定義，主要參數(shù)見表2。夾具為剛體，采用常應(yīng)力實體單元，網(wǎng)格單元尺寸為10 mm，材料選用關(guān)鍵字MAT_RIGID來定義。

圖4 不同的模擬方式Fig. 4 Different simulation methods

表2 AL-6XN鋼板材料參數(shù)Table 2 Material parameters of AL-6XN steel plate

2.2 PBM仿真計算

根據(jù)試驗工況，由模型坐標計算得出爆心位置，采用DEFINE_PBLAST_GEOMETRY關(guān)鍵字定義炸藥位置及形狀特征；采用DEFINE_PBLAST_AIRGEO關(guān)鍵字定義空氣域特征；采用DEFINE_PARTICLE_BLAST關(guān)鍵字設(shè)置C4炸藥材料參數(shù)并定義粒子與結(jié)構(gòu)的耦合關(guān)系，C4炸藥材料參數(shù)見表3?；赑BM算法的有限元模型見圖3(a)。

表3 C4炸藥材料參數(shù)Table 3 Material parameters of C4

2.3 ALE與S-ALE仿真計算

作為對照組，建立了基于傳統(tǒng)ALE和改進S-ALE的完整有限元模型，見圖3(b)及(c)。相較于傳統(tǒng)ALE，LS-DYNA中改進S-ALE使用關(guān)鍵字創(chuàng)建流體單元網(wǎng)格，即建模時僅需創(chuàng)建結(jié)構(gòu)有限元網(wǎng)格，簡化了建模流程，便于模型的修改處理。ALE和S-ALE的模型尺寸、網(wǎng)格劃分方式和單元總數(shù)均相同。流體域采用單點積分的ALE多物質(zhì)單元，輪廓尺寸均為900 mm×900 mm×540 mm。其中，空氣視為初始壓強100 kPa的理想氣體，網(wǎng)格單元尺寸為9 mm，材料選用關(guān)鍵字MAT_NULL定義，密度為1.293 kg/m2，狀態(tài)方程選用關(guān)鍵字EOS_LINEAR_POLYNOMIAL定義，參數(shù)見表4。C4炸藥材料選用關(guān)鍵字MAT_HIGH_EXPLOSIVE_BURN定義，狀態(tài)方程選用EOS_JWL定義，材料參數(shù)同表3，狀態(tài)方程見表5。

表4 空氣狀態(tài)方程參數(shù)Table 4 Parameters of air equation of state

表5 C4炸藥狀態(tài)方程參數(shù)Table 5 Parameters of C4 equation of state

3 計算結(jié)果

3.1 PBM算法計算結(jié)果

PBM算法精度與粒子總數(shù)N和粒子數(shù)量比γ相關(guān)，其中N=NHE+NAIR，γ=NHE/NAIR，NHE為炸藥粒子個數(shù)，NAIR為空氣粒子個數(shù)。依據(jù)模擬經(jīng)驗，每種工況選用6種粒子總數(shù)進行模擬，模擬結(jié)果見圖5。可知當γ不變時，基本上PBM算法精度隨著N值增加而提高。但對于不同模型，都存在精度的增幅率隨N值增加而衰減的現(xiàn)象。

圖5 不同粒子總數(shù)下的跨中撓度Fig. 5 Deflection at different particle counts

為研究粒子數(shù)量比γ對PBM算法精度的影響，考慮計算精度和時間消耗的情況下，以下述三種工況為基礎(chǔ)：①L=150 mm、N=47.5w；②L=200 mm、N=47.5w；③L=250 mm、N=47.5w。每種工況選用5組γ進行仿真模擬。由試算得知三種工況的粒子質(zhì)量比均為MHE∶MAIR=1∶8.5，其中MHE為初始爆炸粒子總質(zhì)量、MAIR為初始空氣粒子總質(zhì)量。參考上述粒子質(zhì)量比設(shè)置γ值，模擬結(jié)果見圖6，可以看出γ對模擬精度的影響十分明顯，當γ值接近粒子質(zhì)量比時，模擬精度較高，γ值偏離粒子質(zhì)量比時，模擬精度較低。

圖6 不同粒子個數(shù)比下的跨中撓度Fig. 6 Deflection at different Particle number ratio

3.2 不同模擬算法結(jié)果對比

3.2.1 模擬精度方面

對比在PBM、ALE和S-ALE算法下，鋼板跨中最終撓度Z、鋼板總能量E及CPU運算耗時t，不同模擬工況結(jié)果見表6。

表6 不同模擬工況結(jié)果Table 6 Results of different simulation conditions

由表6數(shù)據(jù)可以得出，ALE、S-ALE和PBM算法模擬結(jié)果的離散系數(shù)分別為0.16、0.20和0.12。由此可知，試驗中S-ALE算法模擬結(jié)果離散性最大，模擬效果較差；ALE算法次之，但出現(xiàn)了模擬結(jié)果小于試驗值的情況；PBM算法離散性最小，最大誤差不超過20%，模擬效果最好。

此外，爆心距分別為150 mm、200 mm和250 mm時，由PBM算法模擬得到的鋼板總能量E(包括動能及變形產(chǎn)生的內(nèi)能)最大，約為S-ALE算法的1.15倍，約為ALE算法的1.5倍。三種算法得到的E值曲線趨勢大致相同，但是，在PBM和ALE算法E值曲線均趨于平穩(wěn)后，S-ALE算法仍表現(xiàn)出緩慢下降趨勢，如圖7所示。

圖7 鋼板總能量Fig. 7 Total energy of steel plate

比例距離不同的三種工況下，ALE算法的平均CPU運算耗時為21 193 s；S_ALE算法的平均CPU運算耗時為19 220 s；PBM算法的平均CPU運算耗時為1919 s。ALE與S_ALE算法CPU耗時相差不多，PBM算法的CPU耗時約為ALE和S_ALE算法的十分之一。由此可知，PBM算法的運算耗時遠小于ALE算法和S-ALE算法，具有更高的運算效率。

3.2.2 模擬時效方面

模擬采用戴爾臺式工作站，CPU單核基頻3.60 GHz，多核心同時工作可超頻至4.3 GHz。LS-DYNA求解器為R12_MPP求解器，授權(quán)算力為8核芯。L=150 mm工況下，不同算法CPU耗時的帕累托圖見圖8。

圖8 L=150 mm時不同算法CPU耗時占比Fig. 8 CPU time consumption of different algorithms when L=150 mm

對比ALE和S-ALE的CPU總耗時及各處理環(huán)節(jié)耗時占比，發(fā)現(xiàn)S-ALE與ALE的CPU總耗時差距并不明顯，但算法相關(guān)耗時占比上，S-ALE顯著優(yōu)于ALE，可見改進后的S-ALE算法在算法處理機制上進行了優(yōu)化，但由于授權(quán)核心較少并未有效體現(xiàn)S-ALE在多核心運算上的優(yōu)勢。

對PBM算法而言，粒子總數(shù)N對CPU總耗時t影響巨大，兩者基本上呈線性關(guān)系。PBM算法相關(guān)耗時隨著N的增加不斷增加，占比不斷提高，但單元處理和其他耗時基本未發(fā)生變化。在可接受的仿真精度范圍內(nèi)，PBM算法效率比ALE和S-ALE算法效率提高10倍以上。

4 結(jié)論

(1)PBM方法適用于近場爆炸問題的仿真計算。不同工況下，PBM方法仿真得到的鋼板跨中最終撓度與T.Brvik試驗結(jié)果相近，在最大誤差和離散系數(shù)上均優(yōu)于ALE和S-ALE方法的計算結(jié)果。

(2)PBM方法精度受粒子量值影響較大。不同工況下，粒子個數(shù)比一定時，粒子總數(shù)越多，PBM算法模擬精度越高；粒子總數(shù)一定時，粒子個數(shù)比越接近粒子質(zhì)量比，PBM算法精度越高。

(3)相同條件下，PBM算法運算效率更高。PBM算法CPU運算用時與粒子總數(shù)基本呈線性遞增關(guān)系。結(jié)構(gòu)單元劃分方式與數(shù)量相同時，PBM算法的計算用時僅為ALE和S-ALE算法的十分之一，便能得到相同精度的模擬結(jié)果。