煙火藥燃燒粒子流場的實驗研究與數(shù)學(xué)模型

薛 銳, 許厚謙, 李 燕, 朱晨光

(1. 南京理工大學(xué) 能源與動力工程學(xué)院, 江蘇 南京 210094; 2. 南京理工大學(xué) 化工學(xué)院, 江蘇 南京 210094; 3. 南京工程學(xué)院 能源與動力工程學(xué)院, 江蘇 南京 211167)

1 引 言

煙火藥通常是由可燃劑、氧化劑、粘合劑等粉末材料經(jīng)過機械混合、造粒、壓制等工藝制成的復(fù)合型含能材料,是一種典型的非均勻多孔介質(zhì),這種含能材料的燃燒具有高溫、多相和復(fù)雜的中間成分,因此建立煙火藥火焰結(jié)構(gòu)的研究歷來是個難題[1-3]。為此,不少學(xué)者對于該問題進行了深入的研究與探討,埃利·弗里曼研究認為煙火藥燃燒是從初始溫度開始到轉(zhuǎn)變生成火焰溫度時的氣體、液體和固體燃燒產(chǎn)物的復(fù)雜過程[4]; 希洛夫認為煙火藥燃燒火焰結(jié)構(gòu)分還原層、完全燃燒層、熱輻射層、氧氣層和煙層等[5]; 還有柏木(Kashiwagi)建立的均勻物系固體燃料的燃燒模型以及正田強半經(jīng)驗式延期藥燃燒模型[6],都是經(jīng)驗的物理模型,缺乏具體實施數(shù)據(jù)描述,忽略了煙火藥這種粉劑性材料燃燒時火焰中的燃燒粒子,有很大的局限性。

近些年,隨著科學(xué)儀器的發(fā)展進步,人們對煙火藥燃燒的認識發(fā)生較大改變。一方面,使用了大量計算軟件,以化學(xué)反應(yīng)和熱力學(xué)計算出發(fā)描述煙火藥性能,如Farid等[7]通過不同側(cè)面利用熱力學(xué)理論對鎂粉和聚四氟乙烯體系的燃燒機理和紅外輻射效應(yīng)進行了研究; 另一方面,將流體力學(xué)兩相流動理論引入煙火藥燃燒火焰的計算與分析,并借助粒子圖像速度場儀測試設(shè)備(Particle Image Velocimetry，PIV)與高速攝影儀等設(shè)備,發(fā)現(xiàn)煙火藥燃燒火焰由氣體流場與離散燃燒粒子組成,且離散燃燒粒子是不可忽略的重要因素[8-9]。改變了原有觀點,即認為煙火藥燃燒火焰是一種與碳氫燃料一樣的連續(xù)火焰,分成內(nèi)外焰結(jié)構(gòu),炙熱煙塵的影響很難分清火焰的觀點?；鹧娼Y(jié)構(gòu)診斷設(shè)備的不斷更新使火焰結(jié)構(gòu)研究取得了突破性進展,如電荷藕合器件圖像傳感器(Charge Coupled Device,CCD)、平面激光誘導(dǎo)熒光(Planar Laser Induced Fluorescence,PLIF)、可調(diào)諧半導(dǎo)體激光吸收光譜(Tunable Diode Laser Absorption Spectroscopy,TDLAS)等的使用。但煙火藥燃燒產(chǎn)生的大量煙塵,明顯干擾這些設(shè)備對煙火火焰診斷。感光設(shè)備所記錄的圖像及其強度與曝光時間直接相關(guān)。曝光時間長,對于不動點光源,體現(xiàn)為光強累加,對于動點光源,體現(xiàn)為跡線效應(yīng),即在圖像上表現(xiàn)為一個光線段(在曝光時間內(nèi)動點光源移動距離的體現(xiàn))。煙火藥燃燒火焰含有大量高速運動的燃燒粒子,在普通光學(xué)設(shè)備所成圖像表現(xiàn)為連續(xù)的。先進高速攝影儀,可以設(shè)置超短曝光時間,不僅可以濾掉不發(fā)光的燃燒產(chǎn)物顆粒,而且可消除動光源的跡線效應(yīng)[9]。盡管采用高速攝影儀(High-Speed Camera, HSC)可以拍攝火焰的大量燃燒粒子圖像,但這種圖像并不能體現(xiàn)出粒子的速度與空間分布等火焰流場的重要參數(shù)。

為此,本研究對煙火藥火焰中的正在燃燒粒子進行實驗與統(tǒng)計,通過圖片的二值化處理得出正在燃燒粒子的坐標(biāo)位置,通過連續(xù)圖片的粒子追蹤法得到速度分布,獲得其行為軌跡,進而確立火焰粒子速度分布的數(shù)學(xué)模型。

2 實驗

2.1 藥劑與設(shè)備

所用的煙火藥劑由高氯酸鉀(KClO4)、鎂粉(Mg)和硝化棉組成。其中,KClO4為分析純,粒徑為80～100目,分解溫度為400 ℃,國藥集團生產(chǎn); 可燃劑為霧化鎂粉,粒徑為200～325目,熔點648.5 ℃,沸點1090 ℃,唐山威豪公司生產(chǎn); 硝化棉,四川瀘州245廠。

采用丙酮將硝化棉溶解,(按照質(zhì)量比為8∶92混合); 將高氯酸鉀和鎂粉混合均勻后,與硝化棉溶液濕混,并經(jīng)過造粒、干燥,稱量后用30噸油壓機分別壓制在內(nèi)徑為18 mm和10 mm、高度為20 mm的單向開口的柱狀鐵殼內(nèi)用于實驗。其中高氯酸鉀50%～55%,鎂粉35%～40%,硝化棉約10%。

HG100K高速攝影系統(tǒng),柯達REDLAKE公司,170萬像素CMOS感光器,后觸發(fā)連續(xù)記錄,可提供25～100000幅/s的幀比率,圖像最大分辨率為1504×1128。

2.2 實驗

實驗在暗室中進行,裝置如圖1所示。實驗通過調(diào)整HSC與燃燒火焰距離,保證燃燒火焰的燃燒質(zhì)點充滿高速攝影儀視場。根據(jù)燃燒質(zhì)點的速度、亮度等因素,確定曝光時間為15 μs,高速攝影儀的幀頻設(shè)置為500幀/s。

圖1實驗裝置示意圖

Fig.1Schematic diagram of experimental apparatus

實驗采用普通攝像機(最低照度0.1Lux)與高速攝像儀對直徑為18 mm的煙火藥柱的燃燒火焰進行研究,對比火焰圖像中粒子的解析效果,結(jié)果如圖2所示。由圖2a可以看出,普通攝像機的曝光時間處在十幾毫秒到幾百毫秒之間,能記錄火焰氣體的光輻射與發(fā)光粒子輻射,以及那些通過反射和散射的不發(fā)光粒子,如圖像中煙縷與背景,記錄的是整個火焰的外形圖像,但難以分辨火焰中的發(fā)光粒子。由圖2b～圖2e可以看出,隨著曝光時間的縮短,不發(fā)光物質(zhì)與發(fā)光粒子的跡線效應(yīng)被消除,煙火藥燃燒火焰逐漸表現(xiàn)出高密度粒子特征。

圖2中煙火藥柱直徑為18 mm,燃燒粒子密集度偏高,為方便圖像的處理,選擇直徑為10 mm藥柱進行處理,取一組曝光時間為15 μs的圖像作為示例進行研究分析,結(jié)果如圖3所示。

從圖2b～圖2e與圖3可以看出,曝光時間為15 μs時,火焰大部分為粒子簇,僅在噴口高溫區(qū)域仍有顯示一點火焰形狀,對其二值化處理,適當(dāng)選取閾值后,這部分小火焰也變?yōu)榱Ｗ訄D像。

a. the common camerab. 100 μsc. 80 μsd. 50 μse. 10 μs

圖2普通攝像機與高速攝像機(不同曝光時間下)所得煙火藥火焰圖像

Fig.2Images of pyrotechnic flames obtained from the common camera and the high-speed camera under different exposure time

a. the first sampling frameb. the second sampling frame

圖3曝光時間為15 μs的HSC所得煙火藥燃燒火焰圖

Fig.3Succession images of pyrotechnic flames from high-speed camera at exposure time of 15 μs

3 數(shù)學(xué)模型

對圖3進行分析可以發(fā)現(xiàn),圖3中的亮點部分(連續(xù)噴射的過程圖片)為煙火藥燃燒粒子。通過高速攝影儀慢放觀察可知,煙火藥火焰中心軸線且在噴口位置的粒子較快,周邊速度隨著直徑的增大而逐漸降低。隨著燃燒的進行,這些正在燃燒的粒子將向上運動,呈噴射狀,當(dāng)粒子上升速度為零時,由于重力的作用開始逐漸下降。由于高速攝影儀每張圖片只能夠給出每個燃燒粒子二維平面坐標(biāo)和形狀信息,要想獲得燃燒粒子的速度、軌跡及其流場,還需要根據(jù)給出圖面尺寸和連續(xù)拍攝時間間隔,通過建立粒子追蹤模型和軌跡模型,可以求出粒子運動速度的變化過程,進而建立該燃燒過程的數(shù)學(xué)模型。

3.2 建模思路

建模思路如圖4所示。

圖4建模過程

Fig.4The modeling process

3.3 問題分析

(1)二值化處理: 依據(jù)文獻[10]的運動軌跡提取算法對圖片進行二值化處理,高通濾波之后,選擇0.9為二值化的閾值比較合適,影像的灰度高于0.9的令為1,其余為0。直接使用MATLAB 中DIP工具箱函數(shù)im2bw,使用閾值變換法把灰度圖像轉(zhuǎn)換成二值圖像,為下一步用灰度值加權(quán)質(zhì)心計算法[10]求燃燒粒子的坐標(biāo)奠定基礎(chǔ)。

(2)灰度值加權(quán)質(zhì)心計算法: 本研究的粒子定位采用灰度值加權(quán)質(zhì)心計算法,該方法不僅適用于對稱圖像的中心計算,而且其優(yōu)點還在于能夠充分利用圖像中每一點的灰度信息。加權(quán)質(zhì)心計算公式如下:

(1)

式中,x和y分別為質(zhì)心的橫坐標(biāo)和縱坐標(biāo);m為粒子占據(jù)的像素數(shù),且m≥2;f(xi,yi) 為第i個像素的灰度值。基于式(1)的灰度值加權(quán)質(zhì)心計算法,獲取每張圖片粒子的坐標(biāo),并將其標(biāo)注,如圖5所示。同時,對于二值化圖像,我們還得到了其中每個粒子的面積(pixel2),這為粒子軌跡的求取提供了一組用以辨識的參數(shù)。通過文件批處理命令,每張圖片上粒子的特性矩陣D以*.mat格式輸出,便于下面MATLAB調(diào)用。

圖5粒子標(biāo)定及軌跡提取

Fig.5Particle calibration and trajectory extraction

(3)粒子標(biāo)定與軌跡提取: 粒子的標(biāo)定是關(guān)鍵,找到每兩幅圖中的對應(yīng)粒子,采用一組坐標(biāo)值序列標(biāo)定得到軌跡。這里,采用面積和坐標(biāo)作為特征量,對粒子進行標(biāo)定。

首先對第一第二幅圖做粒子的初始標(biāo)定,從第一幅圖片的粒子跟蹤第二幅圖片中的粒子,并以此作為后續(xù)標(biāo)定的初值。對于第一幅圖中的某一粒子坐標(biāo)為(x0,y0),在其幾何臨域(這里采用扇形加矩形區(qū)域)內(nèi)進行搜索,下一幅圖中所有處在該粒子鄰域內(nèi)的粒子均選為候選粒子; 對候選粒子逐個進行分析,在候選粒子中標(biāo)記變量不為1的且選擇面積最相近的標(biāo)定為該粒子下一時刻的狀態(tài),并且置此兩點標(biāo)記變量flag為1以保證不會再次被標(biāo)定; 由兩組坐標(biāo)分別相減得到粒子當(dāng)前速度(pixel/f); 對第一張圖中所有例子重復(fù)上述步驟。

初始時的大范圍搜索圖示,此鄰域判定條件的數(shù)學(xué)表達式:

(2)

這里把圖像進行倒序處理,從后往前進行標(biāo)定,由于在粒子數(shù)密度較低的散射區(qū)域內(nèi),粒子標(biāo)定比較準(zhǔn)確,以此作為初值對高密度的火焰中心區(qū)域進行標(biāo)定,可以有效降低在高密度區(qū)內(nèi)標(biāo)定的偏差。

粒子標(biāo)定之后得到了大量粒子坐標(biāo)的時間序列,把該坐標(biāo)序在坐標(biāo)系下按時序連接,就可以得到粒子的運動軌跡。本研究通過MATLAB程序把這800張圖中所有標(biāo)記過的點的坐標(biāo)、速度,記錄在矩陣?yán)?共記錄了27313組數(shù)據(jù)。

3.4 建立模型

3.4.1 建模條件

(1)假設(shè)所有燃燒粒子形狀為球形。

(2)火焰各區(qū)域溫度變化較大,但是燃燒粒子溫度遠高于其周圍溫度,不僅燃燒粒子不斷產(chǎn)生氣體,而且燃燒粒子呈現(xiàn)非剛性,故假設(shè)忽略熱泳力。

(3)假設(shè)某取樣時間內(nèi)火焰中心單位時間發(fā)出的燃燒粒子數(shù)目恒定。

(4)假設(shè)火焰中心發(fā)射的所有粒子,在以火焰中心為頂點,火焰中心軸線為軸線的一定圓錐角范圍內(nèi),且各個方向的單位時間內(nèi)粒子發(fā)射數(shù)目相同。

(5)假設(shè)經(jīng)過一段時間,系統(tǒng)為相對穩(wěn)態(tài)。

(6)假設(shè)流場氣相流速與燃燒粒子速度差值為一常量。

3.4.2 計算結(jié)果

對燃燒粒子進行受力分析: 主要有重力、浮力和曳力及其他影響。由牛頓第二定律可得下一方程:

(3)

式中，up為粒子速度,m·s-1;t為時間,s;FD,曳力,N;u,氣相速度,m·s-1;gx,重力加速度,m·s-2;ρp,燃燒粒子密度,kg·m-3;ρ,氣相密度,kg·m-3;Fx,其它影響力,N。

式中，dp,燃燒粒子直徑,m;Re,雷諾數(shù);CD,a1,a2,a3為常數(shù)。

根據(jù)假設(shè)二,燃燒粒子在流場中忽略升力,布朗力,熱泳力對運動的影響,故Fx=0。由于流場氣相流速明顯高于燃燒粒子,所以可作假設(shè),流場氣相流速與燃燒粒子速度差值為一常量,記為C=u-up。

經(jīng)過上述假設(shè),可得方程:

(4)

根據(jù)上述受力分析,結(jié)合牛頓第二定律,得出下列常微分方程組:

(5)

3.4.3 問題求解

根據(jù)前文粒子的標(biāo)定與軌跡的提取,可以通過色彩、等高線、粒子數(shù)量等方式表述顯示的粒子速度與密度分布,如圖6所示。

a. particles velocity profile

b. velocity contour map

c. the number density of particles

d. 3D figure of the number density of particles

圖6數(shù)學(xué)模型的煙火藥火焰流場參數(shù)示意圖

Fig.6Schematic diagram of flow field parameters of pyrotechnic flame

從圖6a和圖6b可以看出,粒子在火焰中心區(qū)域速度最大; 伴隨著粒子的上升,粒子的運動速度明顯減小?；鹧嬷姓谌紵Ｗ拥乃俣确植加袀€梭型平衡區(qū),如圖6b內(nèi)環(huán)區(qū)域內(nèi)(每個等速等高線形成的近似圖形),而非預(yù)想粒子速度隨著高度的增加,速度迅速衰減。

從圖6c和圖6d可以看出,粒子在火焰中心區(qū)域最集中,數(shù)量所占比例最大; 隨著燃燒粒子的上升,區(qū)域粒子數(shù)量明顯逐漸減少; 燃燒的粒子上升主要集中在一個以火焰為中心的錐角范圍內(nèi),在此范圍外,幾乎沒有粒子出現(xiàn)。

4 結(jié) 論

基于HSC實驗數(shù)據(jù),采用二值化圖形處理方法,以及粒子運動軌跡模型,得出了煙火藥燃燒粒子的粒子速度、粒子分布的計算模擬圖,并通過實驗數(shù)據(jù)對數(shù)學(xué)模型進行了驗證,并得出如下結(jié)論:

(1)采用灰度值加權(quán)質(zhì)心計算法,并通過倒序處理數(shù)據(jù),可以很好地解決煙火藥粒子的燃燒軌跡;

(2)噴口位置并不是粒子速度的最高區(qū)域。速度分布的等高線在內(nèi)環(huán)區(qū)有個梭型平衡區(qū),證明正在燃燒粒子離開噴口后仍處在加速運動。

(3)由于HSC可拍攝從點火到熄滅,故模型可進行煙火藥燃燒火焰的全過程分析。

參考文獻:

[1] Takeo Shimizu. Study on the reaction mechanism of black powder and its application ballistics of firework shells[J].SelectedPyrotechnicPublicationofDr.TakeoShimizu, Part 2, 45-58.

[2] Kwon Y S, Gromov A A, Ilyin A P, et al. The mechanism of combustion of superfine aluminum powders[J].CombustionandFlame, 2003, 133(4): 385-391.

[3] Zarko V, Gusachenko L K. Simulation of energetic materials combustion[R]. Institute of Chemical Kinetics and Combustion. Russia, 31-January-2000. ADA378994. July, 2002.

[4] Yetter R A, Risha G A, Son S F. Metal particle combustion and nanotechnology[J].ProceedingsoftheCombustionInstitute, 2009, 32(2): 1819-1838.

[5] Shoshin Y L, Dreizin E L. Particle combustion rates for mechanically alloyed Al-Ti and aluminum powders burning in air[J].CombustionandFlame, 2006, 145(4): 714-722.

[6] Christo F C. Thermochemistry and Kinetics models for Magnesium/Teflon/Viton Pyrotechnic Compositions, ADA377423[R]. Weapons Systems Division Aeronautical and Maritime Research Laboratory, 2000.

[7] De yong LV, Smit J J. A theoretical study of the combustion of Magnesium/Teflon/Viton pyrotechnic compositions, ADA243244[R]. Materials Research Laboratory of USA, 2000.

[8] Duraes L, Costa B F O, Santos R, et al. Fe2O3/aluminum thermite reaction intermediate and final products haracterization[J].MaterialsScienceandEngineering, 2007, A 465:199-210.

[9] Zhu C G, Xu C G, Xue R. Study of the spatial distribution of burning particles in a pyrotechnic flame based on particle velocity[J].JournalofEnergeticMaterials, 2014, 32(4): 252-263.

[10] 張健, 周曉東, 張春華. 空間目標(biāo)運動軌跡提取算法研究[J]. 紅外技術(shù), 2007, 29(8): 459-462.

ZHANG Jian, ZHOU Xiao-dong, ZHANG Chun-hua. The research on the algorithm for space target motion trajectory extraction[J].InfraredTechnology, 2007, 29(8): 459-462.