徐維錚 ，吳衛(wèi)國(guó)

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

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

0 引言

針對(duì)約束空間內(nèi)的爆炸問(wèn)題，爆炸后燃燒效應(yīng)是需要重點(diǎn)關(guān)注的內(nèi)爆炸物理現(xiàn)象。大多數(shù)炸藥（TNT炸藥、溫壓炸藥、SDF混合型炸藥［1］等）的爆轟產(chǎn)物具有負(fù)氧性，高溫高壓爆轟產(chǎn)物在膨脹過(guò)程中會(huì)與周?chē)諝庵械难鯕膺M(jìn)行劇烈的燃燒反應(yīng)并釋放大量能量，這一物理現(xiàn)象為后燃燒效應(yīng)［2］。后燃燒效應(yīng)對(duì)爆炸沖擊波的傳播過(guò)程、沖擊波壁面反射壓力、準(zhǔn)靜態(tài)壓力都會(huì)產(chǎn)生一定的影響。

為了研究約束空間內(nèi)后燃燒效應(yīng)對(duì)爆炸載荷的影響規(guī)律，國(guó)內(nèi)外學(xué)者開(kāi)展了大量的實(shí)驗(yàn)和數(shù)值模擬研究。金朋剛等［3］采用壓力傳感器測(cè)量了在密閉罐體內(nèi)分別存在氮?dú)?、空氣和氧?種條件下TNT炸藥在后燃燒過(guò)程中產(chǎn)生的準(zhǔn)靜態(tài)壓力，結(jié)果表明，在氧氣環(huán)境下，TNT炸藥的后燃燒效應(yīng)較明顯，能夠產(chǎn)生更大的準(zhǔn)靜態(tài)壓力。李鴻賓等［4］在容積為500 L的密閉爆炸罐中進(jìn)行了TNT炸藥爆炸實(shí)驗(yàn)，發(fā)現(xiàn)隨著環(huán)境中氧氣量的增大，準(zhǔn)靜態(tài)壓力增大，說(shuō)明提高環(huán)境中的氧氣含量能夠提高爆轟產(chǎn)物的反應(yīng)率。李芝絨等［5］通過(guò)實(shí)驗(yàn)測(cè)量了在圓柱形密閉爆炸罐內(nèi)分別存在空氣和氮?dú)膺@2種條件下溫壓炸藥爆炸的沖擊波峰值以及罐體內(nèi)的準(zhǔn)靜態(tài)壓力，結(jié)果表明，在爆轟反應(yīng)階段，氧氣與超細(xì)鋁粉發(fā)生氧化反應(yīng)，空氣環(huán)境中的沖擊波峰值和沖量比氮?dú)猸h(huán)境的略高；在后燃燒階段，氧氣與鋁粉混合產(chǎn)生后燃燒反應(yīng)，釋放出大量的熱量，與氮?dú)猸h(huán)境相比，空氣環(huán)境中的準(zhǔn)靜態(tài)壓力和熱響應(yīng)溫度峰值顯著增大。Kuhl等［6］提出采用多流體模型描述考慮后燃燒效應(yīng)的爆炸場(chǎng)，將爆炸組分劃分為燃料、空氣和產(chǎn)物3種成分，分別對(duì)這3種成分建立質(zhì)量、動(dòng)量、能量守恒方程，在爆炸過(guò)程中，采用狄拉克函數(shù)進(jìn)行燃燒陣面的捕捉，釋放后燃燒能量，并開(kāi)發(fā)了網(wǎng)格自適應(yīng)數(shù)值計(jì)算程序。Togashi等［7-8］將自主開(kāi)發(fā)的三維爆炸波數(shù)值計(jì)算程序與勞倫斯利弗莫爾實(shí)驗(yàn)室開(kāi)發(fā)的Cheetah編碼銜接，利用Cheetah編碼計(jì)算后燃燒能量，再將后燃燒能量添加到自主開(kāi)發(fā)的程序，實(shí)現(xiàn)爆炸過(guò)程中后燃燒效應(yīng)的模擬。

爆炸流場(chǎng)包含高密度比和高壓力比強(qiáng)間斷等復(fù)雜流場(chǎng)結(jié)構(gòu)，對(duì)爆炸過(guò)程進(jìn)行數(shù)值模擬時(shí)，需要高精度、強(qiáng)穩(wěn)定性激波捕捉格式。Liu等［9］提出了加權(quán)本質(zhì)無(wú)振蕩格式（Weighted Essentially Non-Oscillation Scheme，WENO），Jiang等［10-11］發(fā)展了該格式并擴(kuò)展了其應(yīng)用。目前，WENO格式作為一種典型的高精度激波捕捉格式，對(duì)流場(chǎng)內(nèi)的激波間斷具有較高的分辨率，適用于求解包含激波、膨脹波以及接觸間斷等復(fù)雜結(jié)構(gòu)的流場(chǎng)。

由于后燃燒過(guò)程涉及復(fù)雜的多組分燃燒化學(xué)反應(yīng)，如果考慮詳細(xì)的化學(xué)反應(yīng)過(guò)程，不僅程序編寫(xiě)復(fù)雜，且由于化學(xué)反應(yīng)時(shí)間尺度與流場(chǎng)時(shí)間尺度存在差別，計(jì)算量大，難以應(yīng)用于實(shí)際工程問(wèn)題計(jì)算中。為此，本文擬提出一種考慮后燃燒效應(yīng)的簡(jiǎn)化反應(yīng)率模型，考慮到WENO格式精度高及穩(wěn)定性較好的優(yōu)勢(shì)，基于Fortran平臺(tái)，采用五階WENO有限差分格式，自主開(kāi)發(fā)約束空間內(nèi)考慮爆炸后燃燒效應(yīng)的二維數(shù)值計(jì)算程序，并探討反應(yīng)速率及后燃燒能量大小對(duì)約束空間內(nèi)爆炸載荷的影響規(guī)律。

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

為了近似考慮爆炸后燃燒效應(yīng)，在Miller反應(yīng)率模型［12］思想的啟發(fā)下，提出一種考慮炸藥爆炸后燃燒效應(yīng)的數(shù)值計(jì)算方法，初步探討后燃燒效應(yīng)對(duì)約束空間內(nèi)爆炸載荷的影響規(guī)律。

Miller反應(yīng)率模型最初提出的目的是近似描述含鋁炸藥爆轟過(guò)程中鋁粒子在爆轟波陣面后的反應(yīng)過(guò)程［12］。我們嘗試將模型思想推廣應(yīng)用到約束空間內(nèi)爆炸后燃燒效應(yīng)的數(shù)值計(jì)算中。本文采用的反應(yīng)率模型如下：

式中：a為反應(yīng)速率常數(shù)；p為流體壓力；α為后燃燒過(guò)程中反應(yīng)率（初始時(shí)α=1，反應(yīng)完成后α=0）。

不考慮后燃燒效應(yīng)的可壓縮歐拉方程為

其中，

將反應(yīng)率模型式（1）耦合到式（2）中，并以源項(xiàng)的形式進(jìn)行后燃燒能量的添加，可得考慮后燃燒效應(yīng)的可壓縮歐拉方程為

其中，

以上式中：ρ為密度；u，v分別為x，y方向上的速度分量；E為單位體積流體的總能量；e為比內(nèi)能；Qaf為爆炸后燃燒過(guò)程中單位質(zhì)量釋放的能量，J/kg；γ為氣體的絕熱指數(shù)，文中γ統(tǒng)一取為1.4。

式（5）在每個(gè)方向上均可以看成是一個(gè)帶有源項(xiàng)的雙曲守恒律方程：

例如，針對(duì)x方向，式（8）的半離散守恒型格式為

2 約束空間內(nèi)爆炸后燃燒效應(yīng)數(shù)值計(jì)算

采用五階WENO有限差分格式，自主開(kāi)發(fā)了約束空間內(nèi)考慮爆炸后燃燒效應(yīng)的二維數(shù)值計(jì)算程序。自主程序?qū)τ诒úǖ臄?shù)值計(jì)算可靠性在文獻(xiàn)［14-15］中已經(jīng)得到驗(yàn)證。本節(jié)采用自主程序?qū)s束空間內(nèi)炸藥爆炸過(guò)程進(jìn)行數(shù)值模擬研究，主要探討反應(yīng)率演化過(guò)程及后燃燒能量的加入對(duì)爆炸載荷的影響規(guī)律。

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

長(zhǎng)方形艙室的尺寸如圖1（a）所示，圖中數(shù)值單位為mm。在艙室壁面上設(shè)置了3個(gè)典型測(cè)點(diǎn)，用于監(jiān)測(cè)爆炸載荷時(shí)間歷程和反應(yīng)率時(shí)間歷程。均勻網(wǎng)格步長(zhǎng)取為10 mm，如圖1（b）所示。

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

2.2 初始條件和邊界條件

基于瞬時(shí)爆轟假定，將方形炸藥等效為均勻高壓氣團(tuán)，具體參數(shù)為：邊長(zhǎng)22.2 mm，密度1 630 kg/m3，壓力 3.057 9×109Pa。炸藥設(shè)置在艙室中間，即圖2（a）中的紅色區(qū)域，周?chē){(lán)色區(qū)域?yàn)榭諝庥颍諝饷芏葹?.0 kg/m3，壓力為1.0×105Pa。初始時(shí)刻反應(yīng)率分布見(jiàn)圖2（b）中的紅色區(qū)域。初始時(shí)刻反應(yīng)還沒(méi)發(fā)生，炸藥所在區(qū)域反應(yīng)率為1，艙室其他區(qū)域反應(yīng)率為0。壁面邊界條件設(shè)置為剛性邊界，由于沖擊波與結(jié)構(gòu)變形耦合效應(yīng)十分復(fù)雜，這里暫時(shí)不考慮艙室結(jié)構(gòu)的變形。

2.3 計(jì)算工況

圖2 爆炸初始條件Fig.2 Initial condition for the explosion

為了探討反應(yīng)速率和后燃燒能量大小對(duì)爆炸過(guò)程的影響規(guī)律，選取2類(lèi)典型工況進(jìn)行計(jì)算。工況1：后燃燒單位質(zhì)量釋放的能量Qaf=4.69×106J/kg，保持其他參數(shù)不變，反應(yīng)速率常數(shù)分別取為a=0，10，40。工況2：反應(yīng)速率常數(shù)為a=10，保持其他參數(shù)不變，后燃燒單位質(zhì)量釋放的能量分別取為Qaf=0，3.0×106，4.69×106J/kg。

2.4 反應(yīng)率演化過(guò)程分析

為了探討爆炸后反應(yīng)率在艙室內(nèi)部的演化過(guò)程，給出了2種工況在不同時(shí)刻的反應(yīng)率分布圖，如圖3所示，其中左圖為a=0，Qaf=4.69×106J/kg；右圖為a=10，Qaf=4.69×106J/kg。由圖3可知，2種工況下，反應(yīng)率在不同時(shí)刻總體呈現(xiàn)相似的空間分布形態(tài)。在爆炸工況a=0，Qaf=4.69×106J/kg中反應(yīng)速率常數(shù)為零，即在整個(gè)爆炸過(guò)程中沒(méi)有發(fā)生反應(yīng)，沒(méi)有后燃燒能量的加入。根據(jù)式（4）可知，該工況求解的是爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)的演化過(guò)程；由于爆炸工況a=10，Qaf=4.69×106J/kg中反應(yīng)速率常數(shù)不為零，反應(yīng)速度快，從而導(dǎo)致艙室內(nèi)部的反應(yīng)率快速下降，最終趨近于0。

圖3 不同時(shí)刻的反應(yīng)率分布圖（左：a=0，Qaf=4.69×106J/kg；右：a=10，Qaf=4.69×106J/kg）Fig.3 Reaction rate distribution at different moments（left：a=0，Qaf=4.69× 106J/kg；right：a=10，Qaf=4.69 × 106J/kg）

為定量顯示艙室內(nèi)部反應(yīng)率的變化，圖4給出了工況1（a=0，10，40；Qaf=4.69×106J/kg）壁面測(cè)點(diǎn)1和2的反應(yīng)率時(shí)間歷程曲線(xiàn)。由圖4可見(jiàn)，隨著反應(yīng)速率常數(shù)a的增大，反應(yīng)率迅速下降，反應(yīng)完成后，反應(yīng)率等于0。這里需要說(shuō)明的是，當(dāng)反應(yīng)速率常數(shù)a=0時(shí)，反應(yīng)沒(méi)有發(fā)生，根據(jù)式（4）可知，反應(yīng)率的演化過(guò)程就是爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)的演化過(guò)程，由于爆炸后爆轟產(chǎn)物與艙室內(nèi)部的空氣進(jìn)行了復(fù)雜的摻混過(guò)程，因此艙室內(nèi)部爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)數(shù)值將趨近于大于0而小于1。

圖4 工況1壁面測(cè)點(diǎn)1和2的反應(yīng)率時(shí)間歷程曲線(xiàn)Fig.4 Time history curves of reaction rate for gauging points No.1 and No.2 in explosion case one

圖5給出了工況2（Qaf=0，3.0×106，4.69×106J/kg；a=10）壁面測(cè)點(diǎn)1和2的反應(yīng)率時(shí)間歷程曲線(xiàn)。由圖5可以看出，在反應(yīng)速率常數(shù)一定的情況下，后燃燒能量的加入對(duì)反應(yīng)速率時(shí)間歷程的影響顯著，然而隨著后燃燒能量Qaf的增大，其大小對(duì)反應(yīng)速率時(shí)間歷程的影響較小。

圖5 工況2壁面測(cè)點(diǎn)1和2的反應(yīng)率時(shí)間歷程曲線(xiàn)Fig.5 Time history curves of reaction rate for gauging points No.1 and No.2 in explosion case two

2.5 后燃燒能量的加入對(duì)爆炸載荷的影響

為了探討不同反應(yīng)速率常數(shù)a對(duì)艙室內(nèi)部爆炸載荷的影響規(guī)律，給出了工況1（a=0，10，40;Qaf=4.69×106J/kg）壁面典型測(cè)點(diǎn)3的超壓時(shí)間歷程曲線(xiàn)（圖6）。由圖6可以看出，后燃燒能量的加入明顯增強(qiáng)了沖擊波載荷和沖量；隨著反應(yīng)速率常數(shù)a的增大，沖擊波到達(dá)時(shí)間提前，沖擊波峰值增大，沖量增大。由于后燃燒能量的加入，與工況a=0，Qaf=4.69×106J/kg相比，工況a=10，40；Qaf=4.69×106J/kg的準(zhǔn)靜態(tài)超壓峰值增大。由于后燃燒能量相同，工況a=10，Qaf=4.69×106J/kg和工況a=40，Qaf=4.69×106J/kg具有相同的準(zhǔn)靜態(tài)超壓峰值。

圖6 不同反應(yīng)速率常數(shù)下測(cè)點(diǎn)3爆炸載荷的時(shí)間歷程曲線(xiàn)Fig.6 Blast load time histories with different reaction rate constants for gauging point No.3

為了探討不同后燃燒能量大小Qaf對(duì)艙室內(nèi)部爆炸載荷的影響規(guī)律，給出了工況2（Qaf=0，3.0×106，4.69×106；a=10）壁面典型測(cè)點(diǎn) 3的超壓時(shí)間歷程曲線(xiàn)（圖7）。由圖7可以明顯看出，在反應(yīng)速率常數(shù)一定的情況下，隨著后燃燒能量Qaf的增大，載荷強(qiáng)度增大，最終的準(zhǔn)靜態(tài)超壓峰值增大。

為了初步驗(yàn)證本文后燃燒能量加入的可靠性，從準(zhǔn)靜態(tài)超壓峰值的角度進(jìn)行了對(duì)比分析。根據(jù)文獻(xiàn)［16］可知，封閉艙室內(nèi)部準(zhǔn)靜態(tài)超壓峰值的計(jì)算公式為

圖7 不同后燃燒能量下測(cè)點(diǎn)3的爆炸載荷時(shí)間歷程曲線(xiàn)Fig.7 Blast load time histories with different afterburning energy for gauging point No.3

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

將工況1，2中的后燃燒能量Qaf代入式（10），得到準(zhǔn)靜態(tài)超壓峰值的理論計(jì)算值，再將理論計(jì)算結(jié)果與數(shù)值模擬結(jié)果進(jìn)行對(duì)比，結(jié)果如圖6（a）和圖7（a）所示。由圖可以明顯看出，理論計(jì)算結(jié)果與數(shù)值模擬結(jié)果吻合較好，證明了后燃燒能量加入的可靠性。

3 結(jié) 論

通過(guò)研究得到如下主要結(jié)論：

1）在后燃燒能量大小一定的情況下，反應(yīng)速率常數(shù)增大時(shí)，沖擊波到達(dá)時(shí)間提前，沖擊波峰值、沖量均增大，準(zhǔn)靜態(tài)超壓峰值保持不變。

2）在反應(yīng)速率常數(shù)一定的情況下，隨著后燃燒能量的增大，沖擊波峰值、沖量及準(zhǔn)靜態(tài)超壓峰值均增大，后燃燒能量的加入能顯著增強(qiáng)爆炸載荷強(qiáng)度。

3）本文針對(duì)爆炸后燃燒過(guò)程的數(shù)值計(jì)算，盡管沒(méi)有考慮復(fù)雜的多組分反應(yīng)過(guò)程，然而采用一種簡(jiǎn)化的反應(yīng)率模型近似描述后燃燒能量對(duì)爆炸載荷的影響，不失為一種滿(mǎn)足實(shí)際工程計(jì)算的有效方法。

本文的研究方法及結(jié)果可為進(jìn)一步研究?jī)?nèi)爆炸復(fù)雜多組分后燃燒效應(yīng)及抗爆結(jié)構(gòu)設(shè)計(jì)提供參考和借鑒。