趙象潤，嚴 楠，李朝振，閆利偉，趙變玲，王德君

（1.北京理工大學 爆炸科學與技術(shù)國家重點實驗室，北京100081；2.遼寧北方華豐特種化工有限公司研發(fā)中心，遼寧 撫順113003）

1 引言

隨著數(shù)值仿真技術(shù)的發(fā)展，采用非線性顯式動力學軟件對煙火藥作用過程的可視化建模及分析提供了必要手段［1］。目前常用的非線性顯式動力學仿真軟件（如LS-DYNA，Autodyn，ABAQUS 等）中內(nèi)置的與藥劑有關(guān)的狀態(tài)方程諸如JWL 狀態(tài)方程、γ 律狀態(tài)方程（Gamma Law Equation of State）、點火增長模型（Ignition and Growth Explosive Material Model）等都是描述炸藥高速爆轟的作功過程［2-4］。這些模型并不適合描述煙火藥爆燃的作功過程，因此現(xiàn)階段對于煙火藥爆燃過程的模擬大多采用基于經(jīng)典內(nèi)彈道方程組的理論體系，利用MATLAB 進行數(shù)值計算，獲得燃燒室內(nèi)的壓力時間歷程（p-t 曲線）［5-7］。在進行以煙火藥為驅(qū)動能源的作動裝置（如分離螺母）作用過程的動力學仿真時，將此p-t 曲線作為動力學仿真的輸入條件，代替藥劑的作用過程［8-10］。這種方法涉及兩個不同軟件的聯(lián)合仿真，過程較為繁瑣，而且這種方法不能獲得壓力在燃燒室內(nèi)的空間分布。

采用非線性顯式動力學軟件MSC.Dytran 新開發(fā)的爆燃狀態(tài)方程（Deflagration equation of state）可直接用來模擬煙火藥作功過程。本研究中采用燃速較慢熱穩(wěn)定性良好的炭黑/硝酸鉀（C/KNO3）和燃速較快的鋁/高氯酸鉀（Al/KClO4）兩種煙火藥，利用爆燃狀態(tài)方程模擬定容條件下的p-t 曲線［11-13］。并對比了仿真結(jié)果與密閉爆發(fā)器試驗結(jié)果，驗證仿真方法的準確性和適用性［14］。以期為煙火藥的內(nèi)彈道仿真以及煙火藥驅(qū)動的分離裝置作功過程的動力學仿真提供一種新的方法。

2 實驗部分

2.1 試劑與儀器

炭黑：工業(yè)級，GB3778-2011，天津億博瑞化工有限公司；硝酸鉀：分析純，GB/T 647-2011，北京化工廠；鋁粉：GB/T 2085.2-2007，沈陽東輕鋁材經(jīng)銷有限公司；高氯酸鉀：HG 3-1395-1981，分析純，天津市風船化學試劑科技有限公司。

裝藥量：200 mg；裝藥方式：壓裝。密閉爆發(fā)器：（10±0.5）mL。

壓力傳感器：KISTLER601A 電荷型壓電式壓力傳感器；量程：0～25 MPa；線性誤差：±0.27%；靈敏度：150 pC/MPa；工作溫度：（-196～200）℃；固有頻率：150 kHz；抗沖擊10000 g；過載：50 MPa。

信號調(diào)理器：丹麥戴維創(chuàng)DEWE-PACK-16。動態(tài)分析儀：丹麥戴維創(chuàng)DEWE-5001-TR，8 個通道，每個通道最高采樣速率為10 MS/s/CH，單通道最高采樣速率為25 MS/s/CH，12bit A/D。

壓力藥筒p-t 曲線測試裝置如圖1 所示。

圖1 p-t 曲線測試裝置Fig.1 Test device of p-t curves

2.2 仿真模型的建立

MSC.Dytran 內(nèi)置的爆燃狀態(tài)方程是用于模擬煙火藥或推進劑爆燃過程的數(shù)學模型，可模擬煙火藥燃燒時產(chǎn)生燃氣的壓力和反應(yīng)速率，還能獲得燃燒過程的能量特征（如溫度、能量值等）［15］。該爆燃狀態(tài)方程基于經(jīng)典內(nèi)彈道學，不考慮熱散失，并假定密閉爆發(fā)器不變形、不漏氣。在該模型中，煙火藥作用產(chǎn)生的燃氣的壓力由諾貝爾-阿貝爾狀態(tài)方程定義［16-17］：

式中，p 為壓力，Pa；R 為與氣體組分有關(guān)的氣體常數(shù)，J·（kg·K）-1；T 為爆溫，K；v 為比容，m3·kg-1；b 為余容，表示與單位質(zhì)量氣體分子體積有關(guān)的修正量，m3·kg-1；γ 為比熱容比；E 為單位質(zhì)量藥劑內(nèi)能，J·kg-1。

為描述火藥燃燒過程，假設(shè)：①裝藥的物理化學性質(zhì)均一，幾何形狀相同；②所有藥粒表面同時著火，且燃燒過程按藥粒表面平行層或同心層逐層燃燒的；③所有藥粒的燃燒環(huán)境和條件相同，各方向的燃速相同；④不論火藥燃燒期間還是燃燒結(jié)束后，生成物的成分始終保持不變［18］。設(shè)mp為裝藥初始質(zhì)量，kg；mb為一段時間內(nèi)火藥的燃去量，kg；Vp為裝藥初始體積，m3；V 為一段時間內(nèi)裝藥已燃去體積m3；n 為裝藥單體數(shù)目；ρp為裝藥密度，kg·m-3。則可以得到式（3）［19］：

式中，ψ 為裝藥已燃去質(zhì)量與裝藥初始質(zhì)量之比，為也稱為火藥燃去百分比。將ψ 對時間微分可得式（4）：

為了導出dV/dt，設(shè)某時刻單體藥粒的表面積為Sb，m2，在dt 時間內(nèi)燃去的厚度為de。則得到式（5）：

式中，e 為燃去厚度，m。將e 對時間微分，則得到火藥的線性燃速表達式（6）：

式中，u 為火藥燃燒的線性燃速，m·s-1；聯(lián)立式（5）～（6）則得到式（7）：

由于燃燒涉及多種復雜化學反應(yīng)，實際上難以根據(jù)火藥的化學和物理性質(zhì)預測火藥的線性燃速［20］。在內(nèi)彈道學上，火藥的線性燃速常用圣羅伯特方程（Saint-Robert equation）表示［21］：

式中，w 為燃速系數(shù)，m·s-1·Pa-β，與火藥性能、燃燒壓強范圍、裝藥初溫等有關(guān)；β 為燃速壓強指數(shù)，為量綱為1 的量，對于特定的火藥，在一定壓強范圍內(nèi)是常數(shù)［22］。

火藥在燃燒時快速產(chǎn)生氣體的能力由活度（vivacity）表示為式（9）［23］：

式中，ξ 為藥劑活度，s-1；SA是藥劑初始表面積（Sp）與Vp的比率。在燃燒過程中，隨著碎裂程度的增加，火藥燃速更快，活度更大。SA可以基于顆粒幾何形狀計算初始表面積與體積比，每種幾何形狀都具有特定的SA值［24］。

定義相對燃燒表面σ=Sb/Sp，其物理意義是燃燒期間當前表面積與初始表面積之比，它提供火藥顆粒幾何形狀以影響燃燒部分的燃燒速率［25］，得到式（10）：

因為在分析過程中當前的顆粒表面積是未知的，為了簡化輸入，使用兩個參數(shù)X 和s 輸入到式（11）來近似替代典型的顆粒幾何燃燒表面［26］：

式中，X、s 均為與形狀相關(guān)的參數(shù)，軟件幫助文檔提供了各種形狀的X，s 值供選用。

密閉爆發(fā)器燃燒室簡化為圓柱形，共有15180 個節(jié)點和13664 個網(wǎng)格單元。按照密閉爆發(fā)器的形狀設(shè)置歐拉域，再使用EOSDEF 模型定義炭黑/硝酸鉀的物理和化學屬性。在燃燒期間，燃氣流入歐拉單元。仿真模型所需藥劑輸入?yún)?shù)通過公開文獻或試驗獲取，部分能量參數(shù)應(yīng)用Real 熱力學程序進行求解得到。軟件中每個參數(shù)都有對應(yīng)的域名稱（field heading），仿真輸入的參數(shù)如表1 所示。

3 結(jié)果與討論

3.1 炭黑/硝酸鉀裝藥壓力藥筒升壓過程

炭黑/硝酸鉀的壓力云圖、p-t 曲線如圖2～圖3 所示。由于仿真計算不考慮橋絲升溫點火過程，因此燃燒室從0 s 開始就產(chǎn)生壓力。

圖2 表明密閉爆發(fā)器內(nèi)的壓力傳遞過程和空間分布狀態(tài)，整體上壓力是近似逐層分布的。圖3 可以看出仿真與試驗的p-t 曲線之間的差異。二者的量化對比主要考察峰值壓力（pm）和對應(yīng)的作用時間t 的接近程度。在壓力快速上升階段（從0 到0.5pm），仿真結(jié)果和實測結(jié)果的升壓速率（曲線的斜率）幾乎重合。壓力上升到0.5pm時，仿真與試驗結(jié)果的作用時間僅相差0.03 ms，相對誤差為3.4%；壓力到達pm的時刻，時間差擴大到0.3 ms，相對誤差為4.8%。實測pm比仿真結(jié)果高約0.5 MPa，壓力偏高的原因是由于實測的壓力包括了始發(fā)裝藥（鎂點火藥）產(chǎn)生的壓力。通過單獨的壓力測試得到鎂點火藥自身產(chǎn)生的壓力約為1.1 MPa，可見仿真的結(jié)果實際上偏高約0.6 MPa。這主要是由于仿真分析沒有考慮熱散失，且仿真分析中輸入的能量參數(shù)是由Real 軟件按照化學反應(yīng)完全進行的理想條件下計算得出的。而實際作用過程中，由于制藥工藝、裝藥條件等因素的影響，總會有能量損失。扣除發(fā)火藥壓力貢獻的因素，仿真得到的峰值壓力相對誤差為10%?？梢奙SC.Dytran 對于火藥燃燒過程的模擬準確度較好。

表1 炭黑/硝酸鉀仿真模型輸入?yún)?shù)Table 1 Input parameters of the simulation model of the carbon black/potassium nitrate

圖2 炭黑/硝酸鉀裝藥壓力藥筒在密閉爆發(fā)器內(nèi)的壓力-時間云圖Fig.2 Pressure-time fringe/deformation diagram of the pressure cartridge with carbon black/potassium nitrate charge in closed bomb

仿真和測試p-t 曲線的明顯不同之處在于實測曲線在壓力達到峰值后開始下降，而仿真曲線并沒有下降趨勢，但壓力波動減小。這是由于仿真運算是把整個反應(yīng)過程按照絕熱過程進行的，沒有考慮密閉爆發(fā)器器壁的熱散失，而測試結(jié)果的壓力降一般認為是熱散失導致的［28-30］。仿真曲線在壓力達到峰值后震蕩減小是因為軟件在運算過程中，達到峰值后反應(yīng)已經(jīng)完全結(jié)束，沒有燃氣產(chǎn)生，因而密閉爆發(fā)器內(nèi)的氣體流動趨于穩(wěn)定。

圖3 炭黑/硝酸鉀裝藥壓力藥筒p-t 曲線Fig.3 p-t curves of the pressure cartridge with carbon black/potassium nitrate charge

3.2 鋁/高氯酸鉀裝藥壓力藥筒升壓過程

為驗證仿真方法的準確性和對不同藥劑的適用性，在相同條件下對鋁/高氯酸鉀進行了試驗和仿真（試驗條件與炭黑/硝酸鉀相同）并進行分析，結(jié)果如圖4 所示。

圖4 鋁/高氯酸鉀裝藥壓力藥筒p-t 曲線Fig.4 p-t curves of the pressure cartridge with aluminum/potassium perchlorate charge

圖4 表明，鋁/高氯酸鉀到達壓力峰值的時間約為炭黑/硝酸鉀的1/5，壓力峰值也要比炭黑/硝酸鉀高4.3 MPa。壓力上升到1/2pm時，仿真結(jié)果的作用時間與實測結(jié)果相差0.02 ms，相對誤差為22%；壓力達到pm時，時間差擴大到0.1 ms，相對誤差為7.7%。鋁/高氯酸鉀相對誤差偏大的原因是其反應(yīng)速率比炭黑/硝酸鉀快很多，因此計算時的基數(shù)小導致?？鄢l(fā)火藥壓力貢獻的因素，仿真所得峰值壓力比實測峰值壓力高0.4 MPa，相對誤差為4%。由鋁/高氯酸鉀的仿真也可以看出，MSC.Dytran 對于升壓階段的仿真結(jié)果與實測結(jié)果一致性很好。

上述分析表明，利用MSC.Dytran 的爆燃狀態(tài)方程來模擬定容條件下火藥的作用過程（僅指升壓過程）有較好的準確性。事實上，利用該方法對于變?nèi)葸^程也是適用的，因為軟件本身也可以解算不同體積時的藥劑作用過程并進行迭代。從航天系統(tǒng)的實際應(yīng)用情況來看，使用MSC.Dytran 軟件進行分離螺母作用過程的動力學仿真本身是有較好的準確性的，整個過程也無需進行多軟件的聯(lián)合仿真［31］。

4 結(jié)論

通過對炭黑/硝酸鉀和鋁/高氯酸鉀鋁/高氯酸鉀作用過程的密閉爆發(fā)器試驗和內(nèi)彈道仿真分析，結(jié)果表明：

（1）炭黑/硝酸鉀裝藥壓力藥筒的壓力上升到1/2pm時，仿真結(jié)果與實測結(jié)果的時間差為0.03 ms，壓力達到pm時的時間差為0.3 ms，峰值壓力相對誤差為10%。

（2）鋁/高氯酸鉀鋁/高氯酸鉀裝藥壓力藥筒的壓力上升到1/2pm時，仿真結(jié)果的作用時間與實測結(jié)果時間差為0.02 ms，達到pm時的時間差為0.1 ms，峰值壓力的相對誤差為4%。

使用MSC.Dytran 軟件的爆燃狀態(tài)方程來模擬煙火藥在定容條件下的p-t 曲線有較好的準確性，且可以得到壓力在任意時刻的空間分布。本研究為火藥爆燃過程的內(nèi)彈道仿真分析提供了一種新的方法，該方法可以簡化對于以火藥驅(qū)動的分離裝置的動力學仿真過程，也可以用于同類裝置的p-t 曲線仿真分析。