張瑞華， 靳建偉， 張鄒鄒， 趙宏立， 王瓊林

(西安近代化學(xué)研究所， 西安 710065)

目前，國內(nèi)外對發(fā)射裝藥引起的膛炸機(jī)理已逐步形成共識，即“發(fā)射裝藥低溫—擠壓—破碎—增面—增燃—增壓—膛炸”，發(fā)射裝藥膛內(nèi)擠壓破碎是藥粒的低溫脆性和彈底發(fā)射裝藥著火前受到擠壓和摩擦作用的共同結(jié)果[1-2]，而發(fā)射藥在發(fā)生破碎后引起增面、增燃的燃?xì)馍梢?guī)律計算方法是模擬火炮膛炸亟待解決的難題，這是實現(xiàn)膛炸模擬的關(guān)鍵環(huán)節(jié)。

正常發(fā)射藥的燃?xì)馍梢?guī)律是由幾何燃燒定律推導(dǎo)出的形狀函數(shù)來描述的[3-4]，但是在膛內(nèi)隨機(jī)發(fā)生破碎后的發(fā)射藥顆粒大小不同、形狀各異，傳統(tǒng)的形狀函數(shù)計算方法已不再適用。大量學(xué)者通過試驗或仿真對發(fā)生破碎后的發(fā)射藥燃?xì)馍梢?guī)律進(jìn)行了研究。Horst等[2]考慮發(fā)射裝藥的破碎對膛內(nèi)壓力異常的影響時，將彈底6%的發(fā)射裝藥燃面人為增加為未破裂前的2～5倍，僅在燃速系數(shù)上乘以系數(shù)，發(fā)射裝藥的形狀不改變，應(yīng)用此方法模擬了M110E2榴彈炮M188E1裝藥結(jié)構(gòu)下炮尾膛炸的現(xiàn)象。Gazanas等[5-6]用落錘和高速液壓伺服裝置研究了M30、JA2發(fā)射裝藥的壓縮和撞擊力學(xué)性能，用小型密閉爆發(fā)器測量破碎藥粒的燃燒規(guī)律，得到了破碎藥粒的燃面是未破碎藥粒的6倍的結(jié)論。翁春生等[7-8]在數(shù)值仿真過程中通過增加燃?xì)馍伤俾蔬_(dá)到增加燃面的效果。張小兵等[9]認(rèn)為火藥破碎后按比例破裂成大塊、小塊和粉末狀3種類型，通過試驗獲得不同的顆粒間應(yīng)力下3種碎藥的質(zhì)量分?jǐn)?shù)。芮筱亭等[1]、Rui等[10]、Li等[11]、陳琪等[12]通過理論與試驗獲得了起始動態(tài)活度比定量表征發(fā)射裝藥破碎程度的方法，并建立了破碎發(fā)射藥等效形狀函數(shù)，用于計算破碎發(fā)射藥的燃?xì)馍梢?guī)律，該方法需要通過試驗獲取參數(shù)。Jiang等[13]建立了基于離散單元法的發(fā)射裝藥擠壓破碎程序，通過統(tǒng)計彈簧斷裂個數(shù)來計算破碎發(fā)射裝藥初始燃面。以上研究是通過數(shù)據(jù)統(tǒng)計或者人為設(shè)定的方法預(yù)測發(fā)射藥發(fā)生破碎后的燃?xì)馍梢?guī)律。

目前，離散單元法已經(jīng)較為成熟的應(yīng)用于發(fā)射裝藥由于擠壓、摩擦作用導(dǎo)致的大規(guī)模破碎仿真研究[13-14]中，且獲得實驗驗證，但是缺乏基于離散單元法發(fā)生破碎后發(fā)射藥的燃燒規(guī)律計算方法研究。

現(xiàn)以廣泛應(yīng)用的花邊19孔發(fā)射藥為研究對象，建立考慮花邊和內(nèi)孔的發(fā)射藥離散單元力學(xué)模型，構(gòu)建正三角形燃燒模型和燃燒函數(shù)，并運用燃燒函數(shù)和傳統(tǒng)形狀函數(shù)對正常和破碎程度相同的發(fā)射藥進(jìn)行密閉爆發(fā)器燃燒規(guī)律數(shù)值仿真，驗證燃燒函數(shù)模擬發(fā)射藥發(fā)生破碎后增面、增燃的準(zhǔn)確性。將該計算方法應(yīng)用于基于離散單元法仿真獲得的彈底破碎發(fā)射裝藥中，以期為模擬由發(fā)射裝藥破碎引起的火炮內(nèi)彈道增壓、膛炸現(xiàn)象提供技術(shù)支撐。

1 發(fā)射藥離散單元力學(xué)模型

如圖1所示，花邊19孔發(fā)射藥在縱向呈柱狀，簡化建立二維離散單元力學(xué)模型。依據(jù)真實發(fā)射藥形狀和尺寸，考慮花邊和19個內(nèi)孔，運用EDEM軟件和顆粒替換方法[14-15]，建立如圖2所示的二維花邊19孔發(fā)射藥離散單元力學(xué)模型。離散小球顆粒半徑r為0.067 mm，離散單元數(shù)為8 460，任意相鄰的小球單元之間由黏結(jié)鍵連接，發(fā)射藥的損傷、破壞通過黏結(jié)鍵的斷裂來體現(xiàn)，且小球單元是發(fā)生破碎的最小單位。

圖1 花邊19孔發(fā)射藥Fig.1 The lace 19-hole propellant

圖2 花邊19孔發(fā)射藥離散單元力學(xué)模型Fig.2 Discrete element mechanical model of the lace 19-hole propellant

2 發(fā)射藥燃?xì)馍梢?guī)律仿真理論

2.1 密閉爆發(fā)器燃燒經(jīng)典理論

密閉爆發(fā)器定容情況下的火藥氣體狀態(tài)方程如式(1)所示[1,3]，火藥已燃百分比ψ計算公式如式(2)所示，火藥燃速公式如式(3)所示。通過仿真獲得p-t曲線，即為發(fā)射藥燃?xì)馍梢?guī)律。針對發(fā)生破碎的發(fā)射藥，其難點在于計算已燃百分比ψ。

(1)

(2)

(3)

2.2 建立基于正三角形模型的燃燒函數(shù)

2.2.1 燃燒函數(shù)建模思路

建立發(fā)射藥燃燒函數(shù)的基本思想如圖3所示。為了清晰地表示建模過程，選取圖2中紅框內(nèi)的局部力學(xué)模型進(jìn)行詳細(xì)說明，局部發(fā)射藥燃燒模型建模過程如圖4所示。具體步驟如下。

圖3 燃燒函數(shù)建模思路圖Fig.3 Thought diagram of combustion function modeling

圖4 局部發(fā)射藥燃燒模型建模過程Fig.4 Modeling process of combustion model of the local propellant

步驟1依據(jù)建立的發(fā)射藥離散單元力學(xué)模型，預(yù)設(shè)小球顆粒離散單元代表一個正六棱柱塊體，由于在縱向呈柱狀，簡化建立平面正六邊形幾何模型。

步驟2平面正六邊形幾何模型由6個正三角形組成，建立正三角形燃燒模型，推導(dǎo)燃燒函數(shù)。

步驟3所有正六棱柱塊體都沿著與空氣接觸的燃燒面的平形層逐層燃燒，統(tǒng)計小球顆粒離散單元之間黏結(jié)鍵斷裂情況，依據(jù)建立的燃燒函數(shù)計算出整個藥床的燃燒剩余體積，再計算發(fā)射藥已燃百分比ψ，從而計算獲得密閉爆發(fā)器定容情況下的發(fā)射藥氣體壓力。

2.2.2 建立燃燒函數(shù)

建立燃燒函數(shù)的基本假設(shè)如下。

(1) 所有的藥粒具有均一的理化性質(zhì)。

(2) 塊體與空氣接觸的表面都同時著火。

(3) 所有藥粒具有相同的燃燒環(huán)境，因此燃燒面各個方向上燃燒速度相同。

根據(jù)小球顆粒離散單元之間黏結(jié)鍵斷裂情況，正三角形的燃燒函數(shù)分為7種，如表1所示。

表1 7種正三角形燃燒模型和燃燒函數(shù)Table 1 Seven equilateral triangle combustion models and combustion functions

根據(jù)建立的7種平面燃燒函數(shù)，某一種燃燒函數(shù)表達(dá)式為Si(e) (i=1,2,…,7)；根據(jù)小球顆粒之間黏結(jié)鍵連接情況，正三角形燃燒模型燃燒層數(shù)為k，某一層燃燒模型對應(yīng)的燃燒函數(shù)統(tǒng)計數(shù)量為Ni,j(i=1,2,…,7；j=1,2,…,k)；正六棱柱塊體縱向黏結(jié)鍵連接情況可分為3類：縱向無黏結(jié)鍵連接、縱向連接1個黏結(jié)鍵和縱向連接2個黏結(jié)鍵，對應(yīng)統(tǒng)計數(shù)量分別為n0,i,j、n1,i,j、n2,i,j，滿足Ni,j=n0,i,j+n1,i,j+n2,i,j；正六棱柱塊體的初始高度均為l，則可計算整個藥床的燃燒剩余體積Vall；計算公式如式(4)所示；從而計算發(fā)射藥發(fā)生破碎的已燃百分比ψ，如式(5)所示，將其代入式(1)中用于發(fā)射藥氣體壓力仿真研究。

n2,i,jl]， 0

(4)

(5)

2.3 形狀函數(shù)

為了驗證新建立燃燒函數(shù)的合理性和準(zhǔn)確性，通過傳統(tǒng)的形狀函數(shù)[3]進(jìn)行驗證，形狀函數(shù)如式(6)所示，將其代入式(1)、式(3)中用于發(fā)射藥氣體壓力仿真，并與燃燒函數(shù)仿真結(jié)果進(jìn)行對比分析。

(6)

式(6)中：χ、λ、μ、χs、λs為發(fā)射藥形狀特征量；Z為發(fā)射藥相對燃燒厚度；Zk為燃燒結(jié)束時的相對燃燒厚度。

2.4 起始動態(tài)活度比理論

為了進(jìn)一步驗證新建立燃燒函數(shù)的實用性，引入發(fā)射藥動態(tài)活度比R和起始動態(tài)活度比R0概念[1]。定義密閉爆發(fā)器中燃燒的發(fā)射藥的動態(tài)活度L為

(7)

式(7)中：p(t)為氣體壓力，Pa；pm為發(fā)射藥在密閉爆發(fā)器內(nèi)產(chǎn)生的最大壓力，Pa；dp(t)/dt為壓力變化率，Pa/s。

動態(tài)活度比R為相等壓力下破碎發(fā)射藥與原未破碎發(fā)射藥動態(tài)活度之比。已理論證明破碎發(fā)射藥的動態(tài)活度比等于相同壓力下破碎發(fā)射藥與原未破碎發(fā)射藥的面積比S/S0[1]，發(fā)射藥起始動態(tài)活度比R0為發(fā)射藥被點燃時刻的發(fā)射藥動態(tài)活度比。計算公式為

(8)

式(8)中：L′為破碎發(fā)射藥動態(tài)活度，(Pa·s)-1；L0為原未破碎發(fā)射藥動態(tài)活度，(Pa·s)-1。研究中選取L′為形狀函數(shù)計算獲得的發(fā)射藥動態(tài)活度，L0為燃燒函數(shù)計算獲得的發(fā)射藥動態(tài)活度。

由式(7)可以看出，由密閉爆發(fā)器仿真獲得壓力時間曲線后，即可求得其動態(tài)活度和起始動態(tài)活度比，用來驗證新建立的燃燒函數(shù)計算發(fā)生破碎發(fā)射藥增加燃燒面積的準(zhǔn)確性。

3 發(fā)生破碎的發(fā)射藥燃?xì)馍梢?guī)律仿真及驗證

3.1 發(fā)生破碎的發(fā)射藥離散單元力學(xué)模型

由于建立的燃燒函數(shù)是二維的，為了檢驗平面燃燒函數(shù)統(tǒng)計判斷是否正確、驗證燃燒規(guī)律計算是否有效，通過預(yù)設(shè)發(fā)射藥破碎形式進(jìn)行研究，而在實際研究過程中為擠壓仿真獲得的隨機(jī)破碎發(fā)射藥。預(yù)設(shè)如圖5所示發(fā)生破碎的發(fā)射藥模型，呈中心十字形裂成4塊，建立了如圖5(a)所示的單層發(fā)射藥離散單元力學(xué)模型；建立的局部[圖5(a)紅框內(nèi)力學(xué)模型]正六邊形幾何模型和正三角形燃燒模型分別如圖5(b)、圖5(c)所示；正常和發(fā)生破碎的發(fā)射藥在縱向高度l均為17 mm，縱向黏結(jié)鍵連接情況只存在一種情況：縱向無黏結(jié)鍵連接，即n0,i,j=Ni,j、n1,i,j=0、n2,i,j=0。

圖5 發(fā)生破碎的發(fā)射藥燃燒模型建立過程Fig.5 Establishment process of combustion model of propellant after fracture

3.2 燃燒函數(shù)統(tǒng)計結(jié)果

建立的正常和發(fā)生破碎的發(fā)射藥離散單元力學(xué)模型中離散單元數(shù)量均為8 460，則正三角形燃燒模型數(shù)量為8 460×6=50 760。表2、表3所示為對正常和發(fā)生破碎的發(fā)射藥的燃燒函數(shù)統(tǒng)計結(jié)果，發(fā)射藥的燃燒層數(shù)k為16層，由第1層到第16層逐層燃燒，Ni(i=1,2,…,7)表示某層燃燒的正三角形燃燒模型數(shù)量對應(yīng)的燃燒函數(shù)統(tǒng)計結(jié)果，第1層中未燃燒和已燃燒的共計數(shù)量為50 760，為初始正三角形燃燒模型數(shù)量；從第2層起共計數(shù)量均為上一層未燃燒數(shù)量，表明燃燒函數(shù)統(tǒng)計結(jié)果正確；對比表2和表3，在相同燃燒層數(shù)下，表3中的未燃燒數(shù)量統(tǒng)計結(jié)果均小于表2，表明燃燒層數(shù)相同時，發(fā)生破碎的發(fā)射藥已燃百分比大于未破碎發(fā)射藥。

表2 正常發(fā)射藥燃燒函數(shù)統(tǒng)計結(jié)果Table 2 Statistical results of the combustion function of normal propellant

表3 發(fā)生破碎的發(fā)射藥燃燒函數(shù)統(tǒng)計結(jié)果Table 3 Statistical results of the combustion function of propellant after fracture

圖6所示為圖5(c)中局部正三角形燃燒模型的燃燒層數(shù)判斷模型(前5層)，每層顏色代表某層正三角形燃燒模型，畫圖時5種顏色(黑色、藍(lán)色、紅色、綠色、黃色)依次循環(huán)，空白部分表示發(fā)射藥未燃燒部分，可以看出發(fā)射藥沿著與空氣接觸的面逐層向內(nèi)燃燒，通過圖6(e)中內(nèi)孔細(xì)節(jié)可以看出燃燒模型判斷正確。正常和發(fā)生破碎的發(fā)射藥燃燒層數(shù)判斷模型如圖7所示，可以看出燃燒層數(shù)為16層，分別對應(yīng)表2、表3中的燃燒函數(shù)統(tǒng)計結(jié)果，充分說明燃燒模型統(tǒng)計結(jié)果判斷正確。將表2、表3的統(tǒng)計結(jié)果代入式(4)即可計算發(fā)射藥的燃燒剩余體積Vall，然后通過式(5)、式(1)、式(3)計算密閉爆發(fā)器燃?xì)鈮毫?。通過燃燒函數(shù)和形狀函數(shù)分別進(jìn)行正常和發(fā)生破碎的發(fā)射藥密閉爆發(fā)器仿真研究。

圖6 發(fā)生破碎的發(fā)射藥燃燒層數(shù)示意圖(前5層)Fig.6 Schematic diagram of the number of combustion layers of propellant after fracture (the first five layers)

圖7 發(fā)射藥燃燒層數(shù)示意圖Fig.7 Schematic diagram of the number of combustion layers of propellant

3.3 燃燒函數(shù)驗證

運用燃燒函數(shù)和形狀函數(shù)進(jìn)行密閉爆發(fā)器燃?xì)鈮毫Ψ抡娉跏紖?shù)均相同，仿真裝填密度均為0.2 g/mL。對于正常發(fā)射藥，密閉爆發(fā)器壓力時間對比曲線如圖8(a)所示，曲線吻合較好；依據(jù)式(7)處理得到圖8(b)所示動態(tài)活度曲線，進(jìn)一步依據(jù)式(8)處理得到燃燒函數(shù)與形狀函數(shù)在仿真過程的動態(tài)活度比曲線，如圖8(c)所示，可以看出燃燒表面積之比在1.0左右，表明燃燒函數(shù)可以描述正常發(fā)射藥的燃燒面積和燃?xì)馍梢?guī)律。

圖8 正常發(fā)射藥燃燒規(guī)律仿真結(jié)果Fig.8 Simulation results of combustion laws of normal propellant

當(dāng)運用燃燒函數(shù)，圖5(a)中發(fā)生破碎的發(fā)射藥燃燒初始表面積是正常發(fā)射藥的1.567倍，即起始動態(tài)活度比為1.567；當(dāng)運用形狀函數(shù)，令發(fā)射藥的破碎方式為取0.219倍原發(fā)射藥高度，起始動態(tài)活度比為1.567，達(dá)到兩種仿真方法用的破碎發(fā)射藥初始表面積相同的目的。分別運用燃燒函數(shù)和形狀函數(shù)計算發(fā)生破碎發(fā)射藥的燃?xì)馍梢?guī)律，仿真結(jié)果如圖9所示，發(fā)現(xiàn)在相同起始動態(tài)活度比下壓力時間曲線吻合較好，且燃燒函數(shù)與形狀函數(shù)在仿真過程的燃燒表面積之比值也在1.0左右，同樣表明燃燒函數(shù)可以描述基于離散單元法仿真發(fā)射藥發(fā)生破碎引起的燃?xì)馍梢?guī)律變化。以上結(jié)果驗證了建立的燃燒函數(shù)的有效性和合理性，且具有較高可實施性。

圖9 發(fā)生破碎的發(fā)射藥燃燒規(guī)律仿真結(jié)果Fig.9 Simulation results of combustion laws of propellant after fracture

4 結(jié) 論

(1) 首次通過預(yù)設(shè)離散單元模型代表正六棱柱幾何模型，進(jìn)一步簡化為正三角形燃燒模型，推導(dǎo)和建立了7種燃燒函數(shù)，具有建模方法簡單易懂、燃燒函數(shù)類型少、適用性強等優(yōu)點。

(2) 運用傳統(tǒng)的形狀函數(shù)和新建立的燃燒函數(shù)對正常和破碎程度相同的發(fā)射藥進(jìn)行密閉爆發(fā)器仿真，結(jié)果表明新建立的燃燒函數(shù)不僅能模擬正常發(fā)射藥的燃?xì)馍梢?guī)律，同時能夠精確計算基于離散單元法仿真獲得的發(fā)射藥破碎引起的燃燒規(guī)律變化。

(3) 研究結(jié)果為科研人員提供了一種基于離散單元法的發(fā)射藥發(fā)生隨機(jī)破碎的燃?xì)馍梢?guī)律計算新方法和新思路；為進(jìn)一步進(jìn)行伴隨發(fā)射裝藥擠壓破碎的內(nèi)彈道兩相流仿真從而預(yù)測膛炸現(xiàn)象的研究提供了技術(shù)支撐。