徐自強(qiáng)，王麗娟，朱 潔，李 闖，宋長杰，時(shí)小凡

(1.中煤科工集團(tuán) 沈陽設(shè)計(jì)研究院，沈陽 110168；2.沈陽建筑大學(xué)，沈陽 110168)

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

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

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

1 PBM方法理論

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

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

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

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

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

2 近距離爆炸試驗(yàn)仿真模擬

2.1 試驗(yàn)場(chǎng)景及分析模型

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

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

表1 AL-6XN鋼板近場(chǎng)爆炸試驗(yàn)工況Table 1 Near field explosion test conditions of AL-6XN steel plate

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

根據(jù)試驗(yàn)工況，由模型坐標(biāo)計(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。基于PBM算法的有限元模型見圖3(a)。

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

2.3 ALE與S-ALE仿真計(jì)算

作為對(duì)照組，建立了基于傳統(tǒng)ALE和改進(jìn)S-ALE的完整有限元模型，見圖3(b)及(c)。相較于傳統(tǒng)ALE，LS-DYNA中改進(jìn)S-ALE使用關(guān)鍵字創(chuàng)建流體單元網(wǎng)格，即建模時(shí)僅需創(chuàng)建結(jié)構(gòu)有限元網(wǎng)格，簡化了建模流程，便于模型的修改處理。ALE和S-ALE的模型尺寸、網(wǎng)格劃分方式和單元總數(shù)均相同。流體域采用單點(diǎn)積分的ALE多物質(zhì)單元，輪廓尺寸均為900 mm×900 mm×540 mm。其中，空氣視為初始?jí)簭?qiáng)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 計(jì)算結(jié)果

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

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

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

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

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

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

3.2.1 模擬精度方面

對(duì)比在PBM、ALE和S-ALE算法下，鋼板跨中最終撓度Z、鋼板總能量E及CPU運(yùn)算耗時(shí)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。由此可知，試驗(yàn)中S-ALE算法模擬結(jié)果離散性最大，模擬效果較差；ALE算法次之，但出現(xiàn)了模擬結(jié)果小于試驗(yàn)值的情況；PBM算法離散性最小，最大誤差不超過20%，模擬效果最好。

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

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

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

3.2.2 模擬時(shí)效方面

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

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

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

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

4 結(jié)論

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

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

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