徐維錚，吳衛(wèi)國 ，況正

1高性能艦船技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，湖北武漢430063

2武漢理工大學(xué)交通學(xué)院，湖北武漢430063

0 引言

針對(duì)約束空間內(nèi)爆炸問題，爆炸后燃燒效應(yīng)是需要重點(diǎn)關(guān)注的內(nèi)爆炸物理現(xiàn)象。大多數(shù)炸藥（例如TNT、溫壓炸藥、SDF（Shock-Dispered-Fuel）混合型炸藥［1］等）的爆轟產(chǎn)物具有負(fù)氧性，在膨脹的過程中會(huì)與周圍空氣中氧氣發(fā)生劇烈的燃燒反應(yīng)，釋放出大量能量，形成爆炸后燃燒效應(yīng)［2］。當(dāng)爆炸發(fā)生在約束空間內(nèi)時(shí)，因壁面限制，相比空爆會(huì)進(jìn)一步加劇后燃燒效應(yīng)，從而對(duì)沖擊波的壁面反射壓力、準(zhǔn)靜態(tài)壓力等產(chǎn)生一定的增強(qiáng)效應(yīng)。

迄今，針對(duì)約束空間內(nèi)后燃燒效應(yīng)對(duì)爆炸載荷的影響規(guī)律，國內(nèi)外學(xué)者開展了大量的實(shí)驗(yàn)和數(shù)值模擬研究。國內(nèi)的相關(guān)研究主要集中在西安近代化學(xué)所和西北核技術(shù)研究所。例如，姬建榮等［3］在小型爆炸容器內(nèi)，實(shí)驗(yàn)測量了在密閉和半密閉爆炸工況下TNT-AL炸藥的后燃燒性能，結(jié)果表明當(dāng)鋁粉質(zhì)量分?jǐn)?shù)為20%時(shí)，TNT-AL炸藥在小型爆炸容器中較易發(fā)生后燃燒反應(yīng)。研究人員用高速攝影儀觀測到了后燃燒反應(yīng)及其發(fā)展過程，發(fā)現(xiàn)后燃燒反應(yīng)在爆炸腔內(nèi)形成了一定的準(zhǔn)靜態(tài)壓力，持續(xù)時(shí)間達(dá)數(shù)十毫秒。鐘巍和田宙［4-5］基于理想氣體狀態(tài)方程和能量守恒定律，將準(zhǔn)靜態(tài)壓力劃分為爆炸壓力（等熵假設(shè)/等壓假設(shè)）和化學(xué)反應(yīng)壓力2個(gè)部分進(jìn)行計(jì)算，推導(dǎo)了考慮爆炸產(chǎn)物發(fā)生化學(xué)反應(yīng)后在密閉空間內(nèi)的準(zhǔn)靜態(tài)壓力計(jì)算公式；金朋剛等［6］在密閉容器內(nèi)通過實(shí)驗(yàn)測量了含有3種粒度（0，13，130 μm）的鋁粉HMX基炸藥的準(zhǔn)靜態(tài)壓力，結(jié)果表明，鋁粉可以提高HMX基炸藥在密閉空間爆炸的準(zhǔn)靜態(tài)壓力，小顆粒鋁粉相比大顆粒鋁粉更有利于提高準(zhǔn)靜態(tài)壓力。在國外，美國勞倫斯·利弗莫爾國家重點(diǎn)實(shí)驗(yàn)室的 Kuhl等［1，7-11］采用自主開發(fā)的三維自適應(yīng)網(wǎng)格加密程序（AMR3D），數(shù)值模擬了1.5 g SDF炸藥在不同容積柱形、矩形密閉容器內(nèi)的爆炸過程。該程序采用兩相模型構(gòu)建SDF炸藥模型，即將PETN或TNT炸藥爆轟產(chǎn)物視為連續(xù)介質(zhì)氣相進(jìn)行模擬，將片狀鋁粉視為連續(xù)介質(zhì)固相進(jìn)行模擬，氣相和固相間通過簡化經(jīng)驗(yàn)?zāi)Ｐ蛠砻枋鲑|(zhì)量、動(dòng)量及能量的轉(zhuǎn)換。Balakrishnan等［12-14］提出了氣相和固相耦合的數(shù)值計(jì)算模型，研究了開敞環(huán)境下TNT-AL炸藥爆炸產(chǎn)物的界面失穩(wěn)及其與鋁粒子相互作用的過程，并考慮了鋁粒子與多組分爆轟產(chǎn)物的后燃燒反應(yīng)，其中氣相被視為連續(xù)介質(zhì)來處理，固相被視為離散介質(zhì)并采用拉格朗日方法來捕捉。為近似描述內(nèi)爆炸的后燃燒效應(yīng)，徐維錚等［15］提出了一種應(yīng)用于內(nèi)爆炸數(shù)值計(jì)算的簡化反應(yīng)率模型，探討了模型中能量釋放常數(shù)及后燃燒能量釋放大小對(duì)爆炸超壓載荷的影響規(guī)律。然而，該反應(yīng)率模型中能量釋放常數(shù)需根據(jù)實(shí)際的內(nèi)爆炸過程來確定，而文獻(xiàn)［15］中尚未給出具體的確定方法。為此，本文將提出一種近似預(yù)估反應(yīng)率時(shí)間歷程的理論公式，在此基礎(chǔ)上進(jìn)一步根據(jù)內(nèi)爆炸化學(xué)反應(yīng)的時(shí)間確定能量釋放常數(shù)。

1 控制方程及數(shù)值求解

爆炸后燃燒涉及爆轟產(chǎn)物與周圍空氣復(fù)雜的多組分化學(xué)反應(yīng)過程，若詳細(xì)考慮各組分的反應(yīng)流動(dòng)，不僅程序編寫復(fù)雜，且網(wǎng)格尺寸需設(shè)置得很密，計(jì)算時(shí)間也很長，這些缺點(diǎn)均不便于在實(shí)際工程中快速計(jì)算。為近似考慮爆炸的后燃燒效應(yīng)，文獻(xiàn)［15］中提出了一種簡化反應(yīng)率模型，即

式中：a為能量釋放常數(shù)；α為內(nèi)爆炸后燃燒反應(yīng)率（初始時(shí)α=1，反應(yīng)完成后α=0）；p為氣體壓力；l，n分別為反應(yīng)率指數(shù)和壓力指數(shù)，l=1 2，n=1/6［16］。

將上述反應(yīng)率模型耦合到二維可壓縮歐拉方程，即

其中，

式中：U為守恒變量；F，G分別為 x，y方向上的數(shù)值通量；S為能量源項(xiàng)；ρ為氣體密度；u，v分別為 x，y方向上的速度分量；E為單位體積氣體的總能量；e為氣體內(nèi)能；Qaf為爆炸后燃燒過程中單位質(zhì)量釋放能量；γ為氣體絕熱指數(shù)，本文中γ=1.4。

本文依次采用五階WENO格式和三階TVD龍格—庫塔（TVD-RK）法對(duì)上述方程中的對(duì)流項(xiàng)和時(shí)間項(xiàng)分別進(jìn)行數(shù)值重構(gòu)及數(shù)值離散［17］。

2 內(nèi)爆炸后燃燒反應(yīng)率數(shù)值計(jì)算

本節(jié)將采用自主開發(fā)的程序?qū)s束空間內(nèi)爆炸后燃燒反應(yīng)率的時(shí)間歷程進(jìn)行數(shù)值計(jì)算，主要探討能量加入對(duì)反應(yīng)率的時(shí)間歷程的影響規(guī)律。

2.1 艙室尺寸及測點(diǎn)布置

計(jì)算艙室為長方形，具體尺寸如圖1（a）所示。在炸藥中心處設(shè)置了4個(gè)典型測點(diǎn)來監(jiān)測爆炸載荷及反應(yīng)率的時(shí)間歷程。均勻網(wǎng)格步長取為10 mm，如圖1（b）所示。

圖1 艙室尺寸及網(wǎng)格劃分Fig.1 Sizes of cabin and mesh division

2.2 初始條件和邊界條件

基于瞬時(shí)爆轟假定［18］，將方形炸藥等效為均勻的高壓氣團(tuán)，具體參數(shù)如下：邊長44.4 mm，密度 1 630 kg/m3，壓力 3.057 9×109Pa。炸藥設(shè)置在艙室中間，周圍區(qū)域?yàn)榭諝庥?，空氣密度?.0 kg/m3，壓力為1.0×105Pa。初始時(shí)刻反應(yīng)率分布如下：在炸藥所在區(qū)域，α=1；在艙室其他區(qū)域，α=0。壁面邊界條件設(shè)置為剛性邊界。由于沖擊波與結(jié)構(gòu)變形耦合效應(yīng)十分復(fù)雜，這里暫時(shí)不考慮艙室結(jié)構(gòu)的變形。

2.3 計(jì)算工況

為了探討能量釋放常數(shù)a和后燃燒過程中單位質(zhì)量釋放能量Qaf對(duì)反應(yīng)率時(shí)間歷程的影響規(guī)律，本文選取了2種典型工況進(jìn)行計(jì)算。這里沒有針對(duì)具體的實(shí)際爆炸工況，而是任意選取系列參數(shù)。對(duì)于工況1，Qaf=4.69×106J/kg，并保持其他參數(shù)不變，取4種不同的能量釋放常數(shù)（a=10，20，30，40）。對(duì)于工況 2，a=10，并保持其他參數(shù)不變，取4種爆炸后燃燒過程單位質(zhì)量釋放能量（Qaf=1.0×106，2.0×106，3.0×106，4.69×106J/kg）。

2.4 能量加入對(duì)反應(yīng)率時(shí)間歷程的影響

圖2（a）和圖2（b）分別給出了工況1和工況2下測點(diǎn)4處反應(yīng)率的時(shí)間歷程曲線。

圖2 不同工況下測點(diǎn)4處的反應(yīng)率時(shí)間歷程曲線Fig.2 Time history curves of reaction rate at gauging point-4 for different cases

由圖2（a）可看出，在后燃燒能量Qaf一定時(shí)，能量釋放常數(shù)a對(duì)反應(yīng)率時(shí)間歷程曲線的影響顯著，反應(yīng)率衰減速度隨著a的增大而增大。

由圖2（b）可看出，當(dāng)能量釋放常數(shù)a一定時(shí)，Qaf對(duì)反應(yīng)率時(shí)間歷程的影響較小，在爆炸初期，各反應(yīng)率時(shí)間歷程曲線基本重合，僅在爆炸后期存在一定的差異，且隨著Qaf的增大，反應(yīng)率衰減速度加快。

3 反應(yīng)率時(shí)間歷程理論預(yù)估

由2.4節(jié)的算例可知，后燃燒能量的加入對(duì)反應(yīng)率時(shí)間歷程具有一定的影響規(guī)律，那么能否采用簡化理論對(duì)上述影響規(guī)律進(jìn)行闡釋和預(yù)估？為了對(duì)炸藥中心處測點(diǎn)4的反應(yīng)率時(shí)間歷程曲線進(jìn)行理論預(yù)估，對(duì)式（1）進(jìn)行變量分離，得

由于氣體壓力 p在內(nèi)爆炸過程中隨著時(shí)間的變化而變得非常復(fù)雜，為了近似預(yù)估反應(yīng)率的時(shí)間歷程，并考慮到密閉空間內(nèi)爆炸載荷中準(zhǔn)靜態(tài)壓力占主要成分，得到如圖3所示結(jié)果。

圖3 密閉艙室內(nèi)爆炸超壓載荷時(shí)間歷程曲線（a=10，Qaf=4.69×106）Fig.3 Overpressure time history curves of explosion in closed cabin（a=10，Qaf=4.69×106）

這里選取準(zhǔn)靜態(tài)壓力峰值 pQS代替 p對(duì)式（6）進(jìn)行積分，可得

若令式（7）中的α=0，可得反應(yīng)結(jié)束時(shí)刻tf：

其中，pQS的計(jì)算公式為［19］

式中：m為炸藥質(zhì)量；V為艙室體積；Qtol=QTNT+Qaf，為爆炸過程中釋放的總能量，其中 QTNT為炸藥爆熱，QTNT=4.69×106J/kg；p0為初始大氣壓力；ρE為炸藥密度，ρE=1 630 kg/m3。

由式（7）～式（9）得到密閉空間內(nèi)爆炸炸藥中心處測點(diǎn)4的反應(yīng)率時(shí)間歷程理論預(yù)估公式。主要計(jì)算流程如下：對(duì)某爆炸工況，首先根據(jù)式（9）計(jì)算 pQS，然后根據(jù)式（8）計(jì)算反應(yīng)結(jié)束時(shí)刻tf，最后，根據(jù)式（7）得到反應(yīng)率的理論時(shí)間歷程。

為了驗(yàn)證理論預(yù)估公式的可靠性，對(duì)工況1和工況2炸藥中心處測點(diǎn)4的內(nèi)爆炸后燃燒反應(yīng)率進(jìn)行數(shù)值計(jì)算，得到了如圖4和圖5所示時(shí)間歷程曲線的數(shù)值與理論計(jì)算對(duì)比結(jié)果。由圖4和圖5可知，雖然在爆炸初期和爆炸后期理論值與數(shù)值計(jì)算結(jié)果存在一定的差異（pQS代替 p造成的），但總體趨勢吻合較好。

圖4 工況1測點(diǎn)4的反應(yīng)率時(shí)間歷程曲線數(shù)值與理論值對(duì)比Fig.4 Comparisons of reaction rate time history curves at gauging point-4 between numerical and theoretical results in explosion case 1

圖5 工況2測點(diǎn)4的反應(yīng)率時(shí)間歷程曲線數(shù)值與理論值對(duì)比Fig.5 Comparisons of reaction rate time history curves at gauging point-4 between numerical and theoretical results in explosion case 2

4 能量釋放常數(shù)的確定

在2.3節(jié)計(jì)算工況中，能量釋放常數(shù)a是任意選取的，然而若要模擬真實(shí)爆炸過程中的后燃燒效應(yīng)，還需確定簡化反應(yīng)率模型中的a值。本節(jié)給出一種近似確定方法。根據(jù)式（8），可進(jìn)一步得到a的表達(dá)式為

由式（10）可知，若可以近似預(yù)估密閉艙室內(nèi)爆炸的化學(xué)反應(yīng)結(jié)束時(shí)刻tf，則可確定a值。文獻(xiàn)［20］中給出了基于化學(xué)反應(yīng)動(dòng)力學(xué)確定tf的方法。這里以TNT炸藥艙室內(nèi)爆炸為例，簡要說明確定tf的方法。

TNT炸藥爆炸的化學(xué)反應(yīng)式為

根據(jù)式（11）及艙室內(nèi)爆炸時(shí)TNT的質(zhì)量m和摩爾質(zhì)量MTNT，可得到如下爆轟產(chǎn)物中各組分的濃度：

式中，c為組分濃度，mol/m3。

排除炸藥占據(jù)的體積，可得周圍空氣所占體積VKQ為

根據(jù)艙室內(nèi)VKQ及對(duì) N2，O2這2種組分體積分?jǐn)?shù)（79%，21%）的基本假定，可得到空氣中各組分濃度為

式中，Vmol為氣體的摩爾體積，Vmol=0.024 m3/L。

文獻(xiàn)［20］中給出的典型算例：藥量體積比m/V=0.6 kg/m3、炸藥密度 ρE=1 640 kg/m3。根據(jù)式（12）和式（14）的計(jì)算，可得艙室內(nèi)各組分的初始時(shí)刻濃度如下：

根據(jù)艙室內(nèi)藥量體積比對(duì)化學(xué)反應(yīng)的劃分［4］，可知在此爆炸算例下（m/V=0.6 kg/m3），艙室內(nèi)將依次發(fā)生如下3種化學(xué)反應(yīng)：

針對(duì)式（15）的化學(xué)反應(yīng)1，其反應(yīng)速率常數(shù)k=0.072 6（mol·m-3）-2/ms。這里，反應(yīng)速率常數(shù)根據(jù)Arrhenius定理k=Ae-Ea/RT（其中Ea為活化能，R為理想氣體常數(shù)，R=8.314（J·mol-1）/K，T為熱力學(xué)溫度，A為指數(shù)前因子）計(jì)算得到。

進(jìn)一步，可列出各組分化學(xué)動(dòng)力學(xué)控制方程為

其中，各組分濃度初始條件的計(jì)算結(jié)果為：

通過采用龍格—庫塔法進(jìn)行數(shù)值求解，得到如圖6所示結(jié)果。由圖6可知，該反應(yīng)結(jié)束時(shí)刻tf近似為15 ms。

圖6 反應(yīng)1（式（15））各組分濃度隨時(shí)間變化曲線Fig.6 Concentration of each component time history curves for reaction one（Eq.（15））

同理，針對(duì)式（16）的化學(xué)反應(yīng)2，將其劃分為兩步反應(yīng)進(jìn)行計(jì)算。

式中，反應(yīng)速率常數(shù) k1=0.021 4（mol·m-3）-2/ms，k2=0.022 4（mol·m-3）-2/ms。

進(jìn)一步，可列出各組分化學(xué)動(dòng)力學(xué)控制方程為

式中，nCO，nO2，nO分別為組分CO，O2，O所對(duì)應(yīng)的化學(xué)反應(yīng)級(jí)數(shù)。其中，各組分的化學(xué)反應(yīng)級(jí)數(shù)為

根據(jù)式（18），反應(yīng)后各組分濃度初始條件的計(jì)算結(jié)果為：

通過采用龍格—庫塔法進(jìn)行數(shù)值求解，得到如圖7所示結(jié)果。由圖7可知，該反應(yīng)結(jié)束時(shí)刻tf近似為15 ms。

圖7 反應(yīng)2（式（16））各組分濃度隨時(shí)間變化曲線Fig.7 Concentration of each component time history curves for reaction two（Eq.（16））

同理，針對(duì)式（17）的化學(xué)反應(yīng)3，其反應(yīng)速率常數(shù) k=0.222 4（mol·m-3）-2/ms，進(jìn)一步可列出各組分化學(xué)動(dòng)力學(xué)控制方程為

根據(jù)式（20），反應(yīng)后各組分濃度初始條件的計(jì)算結(jié)果為：

通過采用龍格—庫塔法進(jìn)行數(shù)值求解，得到如圖8所示結(jié)果。由圖8可知，該反應(yīng)結(jié)束時(shí)刻tf近似為5.0 ms。

圖8 反應(yīng)3（式（17））各組分濃度隨時(shí)間變化曲線Fig.8 Concentration of each component time history curves for reaction three（Eq.（17））

由于以上3種反應(yīng)在真實(shí)的內(nèi)爆炸過程中不一定嚴(yán)格按照上述反應(yīng)順序進(jìn)行（可能各個(gè)反應(yīng)同時(shí)在進(jìn)行中），作為對(duì)內(nèi)爆炸時(shí)間尺度的一種估算，將以上3種反應(yīng)時(shí)刻的平均值作為總的化學(xué)反應(yīng)結(jié)束時(shí)刻。由上述結(jié)果可知，在m/V=0.6 kg/m3藥量體積比工況下，tf≈11.7 ms。文獻(xiàn)［20］針對(duì)此工況下的TNT內(nèi)爆炸開展了實(shí)驗(yàn)測量研究，實(shí)驗(yàn)測量的準(zhǔn)靜態(tài)超壓曲線如圖9所示。

圖9 在m/V=0.6 kg/m3下反應(yīng)3的準(zhǔn)靜態(tài)超壓時(shí)間曲線［20］Fig.9 Quasi-static overpressure time history curves of reaction three at m/V=0.6 kg/m3［20］

由圖9可知，在爆炸初期沖擊波壁面反射形成多峰值載荷特征；在爆炸中期，爆炸產(chǎn)物與艙室內(nèi)部氣體充分混合后開始發(fā)生劇烈的化學(xué)反應(yīng)并釋放能量，導(dǎo)致準(zhǔn)靜態(tài)壓力回升；在爆炸后期，化學(xué)反應(yīng)逐漸完成，準(zhǔn)靜態(tài)壓力也趨于穩(wěn)定。此外，化學(xué)反應(yīng)起主要作用的時(shí)間段約為10 ms，這也初步說明了按照上述化學(xué)反應(yīng)動(dòng)力學(xué)分析方法近似確定的化學(xué)反應(yīng)結(jié)束時(shí)刻tf具有一定的合理性和可靠性。

這里需要注意的是，由于確定化學(xué)反應(yīng)時(shí)間具有一定的難度，它涉及詳細(xì)的化學(xué)反應(yīng)分析及其與爆炸流場的耦合；同時(shí)，實(shí)驗(yàn)手段并不能直接測量得到化學(xué)反應(yīng)的時(shí)間歷程，只能通過可測物理量（例如，壓力和溫度）等數(shù)據(jù)進(jìn)行反推或估算；用上述基于簡化的化學(xué)反應(yīng)分析來近似預(yù)估化學(xué)反應(yīng)時(shí)間，仍需要更多的內(nèi)爆炸實(shí)驗(yàn)數(shù)據(jù)進(jìn)行系列驗(yàn)證及分析。

5 結(jié) 論

通過本文的研究，得到如下主要結(jié)論：

1）當(dāng)內(nèi)爆炸后燃燒能量大小一定時(shí)，能量釋放常數(shù)a對(duì)反應(yīng)率時(shí)間歷程曲線的影響顯著，反應(yīng)率衰減速度隨著a的增大而增大。在a一定時(shí)，后燃燒能量大小對(duì)反應(yīng)率時(shí)間歷程的影響較小，在爆炸初期，各反應(yīng)率時(shí)間歷程曲線基本重合，僅在爆炸后期存在一定的差異。

2）所提的用于近似預(yù)估反應(yīng)率時(shí)間歷程的理論公式具有一定的可靠性，作為一種近似方法，在數(shù)值計(jì)算前可近似預(yù)估反應(yīng)率時(shí)間歷程。

3）根據(jù)化學(xué)反應(yīng)動(dòng)力學(xué)分析方法，近似確定化學(xué)反應(yīng)結(jié)束時(shí)刻tf具有一定的合理性和可靠性，按照反應(yīng)率理論預(yù)估公式，可進(jìn)一步確定反應(yīng)率模型中的能量釋放常數(shù)a，這可為模擬真實(shí)的內(nèi)爆炸后燃燒效應(yīng)提供一定的參考和指導(dǎo)。