王成，方鈞宇，楊同會(huì)

（北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100081）

大多數(shù)凝聚相炸藥在持續(xù)熱源的作用下會(huì)發(fā)生熱分解，可能導(dǎo)致燃燒、爆燃甚至爆轟的發(fā)生，造成嚴(yán)重的損失.一般用烤燃試驗(yàn)[1-2]來(lái)測(cè)試并研究炸藥在這種情況下的宏觀響應(yīng)規(guī)律，以減少意外事件的發(fā)生.但烤燃試驗(yàn)的缺點(diǎn)在于成本較高、周期較長(zhǎng)，部分物理量測(cè)量難度較大，因此結(jié)合數(shù)值模擬可以對(duì)烤燃現(xiàn)象進(jìn)行更全面的分析.

現(xiàn)階段對(duì)烤燃過(guò)程的計(jì)算大多依賴(lài)商業(yè)軟件.MAIENSCHEIN 等[3]使用ALE3D 對(duì)B 炸藥和PBXN-109 炸藥的烤燃進(jìn)行了計(jì)算.常雙君等[4]采用了Fluent 進(jìn)行了LLM-105 基戰(zhàn)斗部裝藥烤燃試驗(yàn)的計(jì)算.孫培培等[5]開(kāi)展了烤燃試驗(yàn)，結(jié)果認(rèn)為殼體厚度對(duì)烤燃的點(diǎn)火臨界溫度無(wú)顯著影響；陳科全等[6]基于ABAQUS 軟件針對(duì)不同殼體厚度的烤燃彈進(jìn)行了數(shù)值模擬研究，得出了和孫培培等類(lèi)似的結(jié)論.常見(jiàn)的凝聚相炸藥多為低熔點(diǎn)炸藥，比如TNT 炸藥熔點(diǎn)為352～355 K，遠(yuǎn)低于其熱分解溫度，因此在烤燃過(guò)程中會(huì)發(fā)生從固相到液相的相變[7].陳朗等[8]、寇永鋒等[9]參考了MCGUIRE 和TARVER[10]在1981 提出的用吸熱描述炸藥相變的簡(jiǎn)化模型，并結(jié)合觀察烤燃試驗(yàn)中的升溫曲線，發(fā)現(xiàn)相變過(guò)程不可忽略，并使用流體力學(xué)計(jì)算軟件Fluent 進(jìn)行了計(jì)算.商業(yè)軟件為烤燃過(guò)程的計(jì)算提供了諸多便利，但大多商業(yè)軟件無(wú)法在算法層面優(yōu)化模擬過(guò)程，且這些軟件在計(jì)算炸藥點(diǎn)火后起爆過(guò)程時(shí)稍有欠缺，比如LS-DYNA 在計(jì)算時(shí)需要將點(diǎn)火區(qū)域換算成等當(dāng)量點(diǎn)火藥，與實(shí)際物理過(guò)程有所偏差.

為了能完整且精確地模擬出炸藥慢烤熱點(diǎn)起爆過(guò)程，形成完善且符合實(shí)際物理過(guò)程的慢烤熱點(diǎn)起爆計(jì)算軟件，本文針對(duì)柱形TNT 烤燃彈算例，基于其軸對(duì)稱(chēng)性質(zhì)，發(fā)展了柱坐標(biāo)下關(guān)于熱傳導(dǎo)方程的交替方向隱格式，對(duì)慢烤升溫過(guò)程進(jìn)行模擬.在炸藥熱點(diǎn)火后，采用化學(xué)反應(yīng)流體動(dòng)力學(xué)Euler 方程組描述起爆過(guò)程，使用5 階WENO 格式，編寫(xiě)了炸藥起爆與爆轟仿真程序.并且針對(duì)算例工況進(jìn)行了慢烤試驗(yàn)，對(duì)計(jì)算得出的升溫曲線、點(diǎn)火溫度、點(diǎn)火時(shí)間進(jìn)行了驗(yàn)證.

1 慢烤熱點(diǎn)火起爆控制方程與數(shù)值格式

在慢烤過(guò)程中，凝聚相炸藥主要經(jīng)歷兩大階段：升溫階段以及熱點(diǎn)火后起爆、爆轟階段.劃分兩個(gè)階段的關(guān)鍵物理量為炸藥反應(yīng)速率.在慢烤升溫階段，炸藥的熱分解十分微弱[11]，反應(yīng)速率幾乎為0，炸藥內(nèi)部為低壓零速的簡(jiǎn)單溫度場(chǎng)，采用熱傳導(dǎo)方程即可描述這個(gè)過(guò)程；當(dāng)溫度到達(dá)一定閾值時(shí)，炸藥發(fā)生熱點(diǎn)火，反應(yīng)速率就會(huì)急劇增長(zhǎng)[12], 強(qiáng)烈的熱分解反應(yīng)導(dǎo)致炸藥內(nèi)部出現(xiàn)高溫高壓，粒子高速移動(dòng)，熱傳導(dǎo)方程無(wú)法描述熱點(diǎn)火后炸藥內(nèi)部演化過(guò)程，因此采用化學(xué)反應(yīng)流體動(dòng)力學(xué)Euler 方程組描述熱點(diǎn)火后起爆、爆轟過(guò)程.

1.1 慢烤升溫過(guò)程

1.1.1 慢烤升溫過(guò)程計(jì)算模型

針對(duì)柱形裝藥，凝聚相炸藥慢烤升溫過(guò)程遵循的熱傳導(dǎo)方程為

dλ/dt實(shí)際上為Arrhenius 公式，表示化學(xué)反應(yīng)速率，m為反應(yīng)級(jí)數(shù)，本文采用一級(jí)反應(yīng)模型進(jìn)行計(jì)算：

1.1.2 慢烤升溫過(guò)程計(jì)算方法

為實(shí)現(xiàn)慢烤升溫過(guò)程與熱點(diǎn)火后響應(yīng)過(guò)程計(jì)算的耦合，保證數(shù)據(jù)交換的準(zhǔn)確性，在兩過(guò)程的計(jì)算中空間網(wǎng)格密度應(yīng)當(dāng)保持一致.慢烤升溫過(guò)程持續(xù)時(shí)間較長(zhǎng)，通常為幾十分鐘至數(shù)小時(shí)，若采用顯式格式進(jìn)行計(jì)算，計(jì)算過(guò)程會(huì)受限于CFL 條件，在高密度空間網(wǎng)格下，會(huì)產(chǎn)生極多的無(wú)用計(jì)算，浪費(fèi)計(jì)算時(shí)間與儲(chǔ)存空間.部分無(wú)條件穩(wěn)定的隱式格式，在求解過(guò)程中需要求解大型方程組，計(jì)算效率低.

因此綜合考慮計(jì)算精度、算法穩(wěn)定性和求解效率，發(fā)展了如下針對(duì)慢烤升溫過(guò)程的柱坐標(biāo)下的交替方向隱格式：

以步長(zhǎng)h劃分空間，取時(shí)間步長(zhǎng) τ，定義導(dǎo)溫系數(shù)D=λT/ρCv，將i、j、n分別作為劃分z、r、t的標(biāo)記，其中 0 ≤i≤imax，0 ≤j≤jmax，0 ≤n≤nmax，定義記號(hào)

該求解格式關(guān)于時(shí)間空間都為2 階精度，無(wú)條件穩(wěn)定，可任意設(shè)置時(shí)間步長(zhǎng)，且在求解過(guò)程中產(chǎn)生的方程組全為三對(duì)角方程組，能有效地減少計(jì)算量，提高計(jì)算效率.

TNT 炸藥等低熔點(diǎn)炸藥在慢烤過(guò)程中會(huì)發(fā)生相變，在計(jì)算這類(lèi)炸藥的慢烤升溫過(guò)程時(shí)，導(dǎo)溫系數(shù)D會(huì)隨著溫度而改變，為此，定義液相率 β[8]，并根據(jù)焓值與溫度、相變潛熱的關(guān)系進(jìn)行相變處理.

式中：TS為 炸藥固-液相變起始溫度；TL為炸藥固-液相 變 終 止 溫 度；λTL為炸藥液態(tài)時(shí)的熱導(dǎo)率；λTS為 炸藥固態(tài)時(shí)的熱導(dǎo)率；Qβ為炸藥相變潛熱.

1.2 熱點(diǎn)火后起爆和爆轟過(guò)程

1.2.1 熱點(diǎn)火后起爆和爆轟過(guò)程計(jì)算模型

針對(duì)柱形裝藥，凝聚相炸藥熱點(diǎn)火后起爆、爆轟過(guò)程遵循柱坐標(biāo)系下Euler 方程組，描述為

式中：vr、vz分 別為徑向和軸向的速度；p為壓力；E為單位質(zhì)量的總能；e為比內(nèi)能.

對(duì)于凝聚相炸藥與爆轟產(chǎn)物，通常采用JWL 狀態(tài)方程來(lái)描述壓強(qiáng)、密度、溫度等物理量之間的關(guān)系：

式中：A、B、R1、R2、ω為 常數(shù)；V=ν/ν0表示當(dāng)前比容與初始比容的比值.

在炸藥熱點(diǎn)火后的一段時(shí)間內(nèi)，由于炸藥內(nèi)部壓力與粒子速度較小，反應(yīng)仍由溫度主導(dǎo)，因此仍可采用Arrhenius 公式描述反應(yīng)速率.當(dāng)炸藥內(nèi)部壓力成長(zhǎng)到一定值時(shí)，壓力將替代溫度成為主導(dǎo)化學(xué)反應(yīng)的因素，此時(shí)采用LEE 和TARVER 提出的點(diǎn)火增長(zhǎng)模型描述反應(yīng)速率：

式中I、a、b、c、d、e、g、x、y、z、G1、G2、λ1、λ2、λ3均為常數(shù).3 個(gè)開(kāi)關(guān)函數(shù)分別控制三項(xiàng)式點(diǎn)火增長(zhǎng)模型中的熱點(diǎn)成核項(xiàng)、慢速反應(yīng)項(xiàng)、快速反應(yīng)項(xiàng).TNT 炸藥等低熔點(diǎn)炸藥在慢烤升溫過(guò)程中會(huì)發(fā)生固-液相變，導(dǎo)致炸藥內(nèi)部無(wú)孔穴等結(jié)構(gòu)缺陷，熱點(diǎn)成核項(xiàng)不起作用，因此在計(jì)算這類(lèi)問(wèn)題時(shí)直接關(guān)閉熱點(diǎn)成核項(xiàng)，取 Φ1=0.

因此采用如下修正后的反應(yīng)模型計(jì)算凝聚相炸藥烤燃熱起爆和爆轟過(guò)程：

在該模型中，當(dāng)溫度主導(dǎo)化學(xué)反應(yīng)時(shí)，采用Arrhenius 公式計(jì)算反應(yīng)速率；當(dāng)壓力主導(dǎo)化學(xué)反應(yīng)時(shí)，通過(guò)采用點(diǎn)火增長(zhǎng)方程計(jì)算反應(yīng)速率.

1.2.2 熱點(diǎn)火后起爆和爆轟過(guò)程計(jì)算方法

在對(duì)熱點(diǎn)火后起爆與爆轟過(guò)程的計(jì)算中，空間離散采用5 階WENO 格式[13-14]，通過(guò)引入自適應(yīng)權(quán)重，在強(qiáng)間斷附近自動(dòng)選擇最光滑的模版進(jìn)行通量重構(gòu)，有效抑制了高精度計(jì)算中的非物理振蕩，時(shí)間離散采用總變差遞減的3 階Runge-Kutta 法，進(jìn)一步增強(qiáng)數(shù)值格式的穩(wěn)定性.

熱傳導(dǎo)方程主要對(duì)溫度進(jìn)行了計(jì)算，而Euler方程組涉及到壓力、速度、能量等物理量，因此在計(jì)算熱點(diǎn)火后起爆與爆轟過(guò)程前，應(yīng)當(dāng)對(duì)升溫過(guò)程的計(jì)算結(jié)果數(shù)據(jù)進(jìn)行預(yù)處理.為此，設(shè)計(jì)如下數(shù)據(jù)接口：

式中：下標(biāo)Euler 與heat 分別表示該變量在Euler 方程組與熱傳導(dǎo)方程中的取值；ε為小量，一般取0.01；αg為爆轟產(chǎn)物與炸藥比熱容的比值.

表1[8-9,15]給出了TNT 炸藥計(jì)算過(guò)程中涉及到的的所有相關(guān)參數(shù).

表1 TNT 炸藥相關(guān)參數(shù)Tab.1 Related parameters of TNT

2 TNT 炸藥烤燃試驗(yàn)

為驗(yàn)證算法可行性，設(shè)計(jì)并開(kāi)展TNT 炸藥烤燃試驗(yàn).試驗(yàn)采用K 型WRNK-191 溫度傳感器進(jìn)行測(cè)溫，HIOKI LR8431-30 數(shù)據(jù)記錄儀用于溫度顯示與采集.在殼體表面均勻纏繞電熱線圈，外層包裹石棉保溫層，圖1 為T(mén)NT 烤燃彈加工前后圖.

圖1 TNT 烤燃彈彈體Fig.1 TNT cook-off projectile

圖2 為T(mén)NT 烤燃彈彈體柱截面示意圖，藥柱直徑為60 mm，高為240 mm，內(nèi)部炸藥為T(mén)NT 炸藥，裝藥密度為1 624 kg/m3，外部殼體厚度為3 mm，材料為45 號(hào)鋼.

圖2 烤燃彈算例截面尺寸圖Fig.2 Example section size diagram

在測(cè)溫點(diǎn)A、B、C處放置溫度傳感器，其中A處溫度傳感器測(cè)得的溫度數(shù)據(jù)用于分析溫度變化規(guī)律，B、C處溫度傳感器測(cè)得的溫度數(shù)據(jù)與溫控系統(tǒng)進(jìn)行實(shí)時(shí)交互，以控制升溫速率.

烤燃試驗(yàn)起始環(huán)境溫度為284 K，升溫速率為1 K/min.位于測(cè)溫點(diǎn)A處的溫度傳感器記錄的溫度數(shù)據(jù)于236 min 終止，記錄的最高溫度為509 K.此時(shí)從監(jiān)控畫(huà)面觀測(cè)，TNT 烤燃彈已開(kāi)始冒煙，溫控系統(tǒng)顯示殼體溫度為520 K，隨后TNT 烤燃彈發(fā)生劇烈反應(yīng)，如圖3 所示.

3 計(jì)算結(jié)果與分析

3.1 慢烤升溫過(guò)程計(jì)算結(jié)果與對(duì)比分析

針對(duì)TNT 炸藥慢烤升溫過(guò)程進(jìn)行數(shù)值模擬計(jì)算，并與試驗(yàn)測(cè)得數(shù)據(jù)進(jìn)行對(duì)比，圖4 為慢烤試驗(yàn)數(shù)據(jù)與計(jì)算結(jié)果對(duì)比.從圖4(a)所示的計(jì)算結(jié)果曲線可以看出，測(cè)溫點(diǎn)A處炸藥升溫曲線在112 ～122 min出現(xiàn)相變平臺(tái)，升溫速率減小至0.46 K/min，溫度由352 K 升至355 K；隨后升溫速率躍升至2.7 K/min；直至146 min 時(shí)，升溫速率逐漸減小至1.0 K/min.測(cè)溫點(diǎn)A處試驗(yàn)結(jié)果曲線中，108 min 時(shí)，柱心傳感器處的升溫速率減小至0.5 K/min 左右；試驗(yàn)進(jìn)行至125 min時(shí)，柱心傳感器處的升溫速率升高至6.0 K/min 左右；試驗(yàn)進(jìn)行至133 min 時(shí)，柱心傳感器處的升溫速率恢復(fù)至1.0 K/min 左右直至紀(jì)錄終止.可以看出測(cè)溫點(diǎn)A處的計(jì)算結(jié)果在升溫趨勢(shì)、相變平臺(tái)上與試驗(yàn)高度一致，在相變段的最大相對(duì)誤差為5%左右，在非相變段的最大相對(duì)誤差為3.5%左右.

圖4(c)所示的計(jì)算結(jié)果中，殼體的升溫速率保持在1 K/min 直至233 min，此時(shí)殼體溫度為525 K；隨后，殼體升溫速率不斷增大，在236 min 左右，殼體溫度直線上升，說(shuō)明炸藥在236 min 左右發(fā)生熱點(diǎn)火并起爆.在試驗(yàn)中，殼體升溫速率在235 min 時(shí)由520 K 開(kāi)始快速升高，與計(jì)算結(jié)果吻合度較高.

圖4 TNT 烤燃試驗(yàn)數(shù)據(jù)與計(jì)算結(jié)果對(duì)比圖Fig.4 Comparison of experimental data and calculated results

圖5 為236 min 時(shí)TNT 烤燃彈溫度分布圖，此時(shí)由于TNT 炸藥熱分解產(chǎn)生的熱量累計(jì)，彈體內(nèi)部2個(gè)圓環(huán)區(qū)域的溫度高于殼體，TNT 炸藥發(fā)生熱點(diǎn)火，點(diǎn)火時(shí)彈體內(nèi)最高溫度大約為550 K.點(diǎn)火區(qū)域是z=10 mm 與z=230 mm 位 置 處 半 徑 為20 mm 的2 個(gè)點(diǎn)火圓環(huán)，圓環(huán)內(nèi)徑大約為3 mm.

圖5 236 min 時(shí)TNT 烤燃彈內(nèi)部溫度云圖Fig.5 Cloud diagram of temperature in TNT projectile at 236 min

3.2 熱點(diǎn)火后起爆和爆轟過(guò)程計(jì)算結(jié)果與分析

以圖5 對(duì)應(yīng)的慢烤模擬結(jié)果為初始溫度場(chǎng)，通過(guò)1.2.2 中的數(shù)據(jù)接口設(shè)置初始條件，同時(shí)設(shè)置剛性壁面條件，開(kāi)展熱點(diǎn)火后起爆、爆轟過(guò)程數(shù)值模擬.圖6、圖7 分別為T(mén)NT 炸藥熱點(diǎn)火后典型時(shí)刻內(nèi)部溫度、壓力云圖，高溫區(qū)域逐漸擴(kuò)大，且溫度由600 K持續(xù)上升至1 800 K.高壓區(qū)域伴隨高溫區(qū)域產(chǎn)生，在2.83 μs 時(shí)，高溫區(qū)域產(chǎn)生了2.6 GPa 的高壓，隨后產(chǎn)生壓力波向柱心傳播，在3.80 μs 時(shí)，壓力波達(dá)到6.5 GPa，并接近柱心.圖8 為在4.28 μs 時(shí)彈體內(nèi)部壓力云圖，此時(shí)沖擊波在柱心匯集，產(chǎn)生26 GPa 的高壓.

圖6 熱點(diǎn)火后典型時(shí)刻烤燃彈體內(nèi)溫度云圖Fig.6 Cloud diagram of temperature in projectile at typical time

圖7 熱點(diǎn)火后典型時(shí)刻烤燃彈體內(nèi)壓力云圖Fig.7 Cloud diagram of pressure in projectile at typical time

圖8 4.28 μs 時(shí)烤燃彈體內(nèi)壓力云圖Fig.8 Cloud diagram of pressure in projectile at 4.28 μs

圖9 為T(mén)NT 烤燃彈熱點(diǎn)火后壓力、反應(yīng)率隨時(shí)間變化曲線.曲線對(duì)應(yīng)的點(diǎn)處于柱心對(duì)稱(chēng)軸上z=10、30、50、70、90、110 mm 位置處.從計(jì)算結(jié)果可以看出，柱心對(duì)稱(chēng)軸上z=10 mm 位置處的TNT 炸藥在3.85 μs 時(shí)達(dá)到壓力峰值，峰值大約為15 GPa，此時(shí)該位置反應(yīng)率接近100%，能量完全釋放.結(jié)合點(diǎn)火圓環(huán)位置，推測(cè)出由點(diǎn)火圓環(huán)產(chǎn)生的壓力波在柱心完成匯聚的時(shí)間點(diǎn)在熱點(diǎn)火后3.85 μs 左右.柱心對(duì)稱(chēng)軸上z=30 mm位置處的TNT 炸藥在4.94 μs 時(shí)達(dá)到壓力峰值，峰值大約為33.7 GPa；柱心對(duì)稱(chēng)軸上z=50、70、90、110 mm 位置處的TNT 炸藥幾乎同時(shí)在5.74 μs 時(shí)達(dá)到壓力峰值，且峰值幾乎一致，大約為33.7 GPa.這是因?yàn)楦邏翰嚸婧笙∈璨ㄅc反應(yīng)流達(dá)到平衡，峰值保持穩(wěn)定，說(shuō)明TNT 炸藥已形成穩(wěn)定爆轟，起爆時(shí)間大約為4.94 μs.

圖9 熱點(diǎn)火后典型位置壓力、反應(yīng)率變化曲線Fig.9 Pressure and reaction rate curves at typical locations

4 結(jié) 論

結(jié)合試驗(yàn)與模擬的結(jié)果，得出以下主要結(jié)論：

①本文推導(dǎo)的交替方向隱格式在計(jì)算TNT 炸藥慢烤時(shí)，具有較高精度且實(shí)現(xiàn)了相變，柱心傳感器處的試驗(yàn)與模擬結(jié)果相對(duì)比，在非相變段的最大相對(duì)誤差為3.5%左右，在相變段的最大相對(duì)誤差為5.0%左右.計(jì)算得出點(diǎn)火區(qū)域?yàn)閦=10 mm 與z=230 mm位置處半徑為20 mm 的2 個(gè)點(diǎn)火圓環(huán)，圓環(huán)內(nèi)徑大約為3 mm.

②TNT 烤燃彈熱點(diǎn)火起爆前殼體臨界溫度為520 K，柱心傳感器處溫度為509 K.

③熱點(diǎn)火后，點(diǎn)火區(qū)域產(chǎn)生壓力波向柱心傳播，于3.85 μs 左右在柱心匯聚，產(chǎn)生15 GPa 的高壓；TNT烤燃彈起爆時(shí)間大約為4.94 μs，形成穩(wěn)定爆轟后，中心軸上峰值壓力大約為33.7 GPa.