固體火箭發(fā)動機(jī)尾焰中顆粒速度的粒子軌跡法測量①

相恒升,程 博,張成飛,喬文生,徐海濤

(1.清華大學(xué) 能源與動力工程系及燃燒能源中心,北京 100084; 2.內(nèi)蒙航天動力機(jī)械測試所,呼和浩特 010010;3.清華大學(xué) 航天航空學(xué)院及燃燒能源中心,北京 100084)

0 引言

導(dǎo)彈及某些空間飛行器多使用具有特定性能的復(fù)合固體推進(jìn)劑,使用這類推進(jìn)劑時一般會在其中加入金屬或非金屬顆粒以提高推進(jìn)劑的比沖以及抑制震蕩燃燒[1-3]。然而,添加金屬或非金屬顆粒也會帶來一些不利影響:推進(jìn)劑燃燒產(chǎn)生的凝聚相燃燒產(chǎn)物導(dǎo)致發(fā)動機(jī)出現(xiàn)兩相流損失[4-5];金屬顆粒也會加重射流對噴管的侵蝕[6-9],如 1993 年美國 Titan IV 運(yùn)載火箭因為顆粒對噴管的侵蝕在點火后約 100 s發(fā)生了爆炸[10]。發(fā)動機(jī)尾焰中顆粒速度的實測數(shù)據(jù)可以用來驗證固體火箭發(fā)動機(jī)中兩相流的理論分析和數(shù)值計算[11],為新一代發(fā)動機(jī)能量損失模型構(gòu)建提供數(shù)據(jù)支撐。因為顆粒慣性正比于粒徑的平方,大顆粒的慣性更大,所以大顆粒對于噴管的侵蝕作用也更明顯。此外,發(fā)動機(jī)尾焰流場中存在激波引起的氣流速度突變,這種突變對小顆粒和大顆粒的影響也可能完全不同。因此,測量尾焰中大尺寸顆粒的速度對于正確認(rèn)識發(fā)動機(jī)尾焰流動具有重要的實際意義。

然而,測量固體火箭發(fā)動機(jī)尾焰中的顆粒速度存在諸多困難,主要體現(xiàn)在以下幾點:(1)流場溫度高,固體火箭發(fā)動機(jī)尾焰中的溫度超過2000 K[12-13];(2)尾焰流場速度高[14];(3)尾焰中顆粒濃密[15];(4)火箭發(fā)動機(jī)尾焰為超音速氣流并伴有強(qiáng)烈的溫度及壓強(qiáng)變化,尾焰流場中氣體密度變化較大[16],并且尾焰邊緣存在強(qiáng)烈的密度梯度,這導(dǎo)致的折射率變化使得經(jīng)過尾焰的光線發(fā)生復(fù)雜的偏轉(zhuǎn)。這些因素給發(fā)動機(jī)尾焰測量帶來了很大困難,目前能夠測量的固體火箭發(fā)動機(jī)內(nèi)外流場參數(shù)比較少,幾乎只能測量發(fā)動機(jī)推力曲線和發(fā)動機(jī)前端壓力隨時間變化的曲線[16]。發(fā)動機(jī)尾焰中顆粒速度測量受制于這些困難[17-18],常見的粒子速度測量技術(shù),如PIV、LDA、LTV等難以應(yīng)用到真實固體火箭發(fā)動機(jī)尾焰流場中顆粒速度的測量[19]。PIV雖然可以實現(xiàn)全場顆粒速度測量,但在發(fā)動機(jī)尾焰中因顆粒數(shù)密度高PIV將失效;LDA因為入射激光經(jīng)過高溫尾焰后產(chǎn)生不可控制的偏轉(zhuǎn)而無法形成測量區(qū),導(dǎo)致LDA無法測量;LTV的測量區(qū)域每次只能有一個顆粒,同樣因為尾焰中顆粒濃密導(dǎo)致該技術(shù)失效。此外,發(fā)動機(jī)尾焰中高溫輻射對LTV是強(qiáng)烈的信號噪音,這也將會導(dǎo)致LTV無法測量。

鑒于固體火箭發(fā)動機(jī)尾焰中顆粒速度測量比較困難,發(fā)動機(jī)尾焰中顆粒速度測量的相關(guān)文獻(xiàn)比較少。PAGLIAROLI等[14]利用兩個光電探頭采集尾焰中顆粒散射信號,通過對兩個信號做條件互相關(guān)計算得到顆粒通過兩個光電探頭探測區(qū)所用時間,兩光電探頭探測區(qū)間距已知,可求得顆粒速度。BALAKUMAR與ADRIAN[16]利用 PIV 技術(shù)測量了噴喉直徑為 2.7 mm 的模型固體火箭發(fā)動機(jī)尾焰流場,流場中最高速度為 630 m/s。SAILE等[17]利用 PIV 技術(shù)對一系列不同噴喉直徑(噴喉直徑分別為 7.75、6.75、5.74 mm)、使用不同配方推進(jìn)劑及不同燃燒室壓力的縮比固體火箭發(fā)動機(jī)尾焰流場進(jìn)行了測量,PIV 測量區(qū)域約為 30 mm×30 mm,實驗發(fā)現(xiàn),即使對如此小尺寸的發(fā)動機(jī)進(jìn)行 PIV 尾焰流場測量,對某些工況(推進(jìn)劑 HTPB1814),PIV 測量仍然會失敗,作者將 PIV 測量失敗的原因歸結(jié)為流場中顆粒濃度太高及背景輻射太強(qiáng)。

由于尾焰中大顆粒的輻射強(qiáng)度也遠(yuǎn)高于小顆粒。本文提出利用大顆粒輻射強(qiáng)這一特點,使用高速相機(jī)對火箭發(fā)動機(jī)尾焰進(jìn)行較長曝光時間的直接拍攝,獲得相機(jī)曝光時間內(nèi)大尺寸顆粒形成的粒子軌跡圖像,采用粒子軌跡法[20-21]對尾焰中的大尺寸顆粒的速度進(jìn)行測量。粒子軌跡法測速的關(guān)鍵是準(zhǔn)確識別粒子軌跡長度,為此本文提出了一種一維擬合算法,以獲得具有亞像素精度的粒子軌跡長度,并采用合成粒子軌跡圖像來檢驗此算法的精度與計算效率。之后使用本文發(fā)展的圖像處理程序采用粒子軌跡法對某固體火箭發(fā)動機(jī)尾焰中大尺寸顆粒的速度進(jìn)行了全場測量,分析了尾焰中大尺寸顆粒速度特點。

1 測量原理及數(shù)據(jù)處理

1.1 粒子軌跡法速度測量原理

一般情況下,使用相機(jī)對發(fā)光顆粒攝像記錄時,顆粒在相機(jī)成像平面上的像為圓斑。然而當(dāng)顆粒高速運(yùn)動時,即使曝光時間很短,顆粒的像在相機(jī)成像平面上的位移可能超過相機(jī)的最小成像單元(對于數(shù)字相機(jī)而言即一個像素,pixel)。此時粒子在相機(jī)中的像不是圓斑,而是圓斑在成像平面上不斷移動形成的拖影,即粒子軌跡,如圖1所示。

圖1 粒子軌跡法測速原理Fig.1 Principle of particle streak velocimetry(PSV)

設(shè)相機(jī)成像的空間分辨率為Rs(m/pixel)、曝光時間為Δt(s),粒子軌跡法即為利用數(shù)字圖像處理技術(shù)識別出圖像中的獨立的粒子軌跡并計算出其長度lp(pixel),因為曝光時間很短,可以假設(shè)曝光時間內(nèi)粒子速度不變,則粒子軌跡為直線,所以粒子速度vp(m/s)可以表示為

(1)

因此,粒子軌跡法特別適合于測量高速運(yùn)動的粒子速度。對于粒子軌跡法而言,相機(jī)的曝光時間可以從相機(jī)參數(shù)中準(zhǔn)確讀取,分辨率由成像系統(tǒng)決定,可以在測量之前標(biāo)定,因此提高粒子軌跡法速度測量精度的關(guān)鍵是提高粒子軌跡長度lp計算精度。

1.2 亞像素精度粒子軌跡長度算法及精度檢驗

1.2.1 亞像素精度粒子軌跡長度算法

在得到含有粒子軌跡的圖像后,首先需要判別出粒子軌跡。對于數(shù)字圖像而言,這意味著判別出哪些像素屬于某個特定的粒子軌跡,即所謂的圖像分塊(segmentation)。通常實驗所得的圖像為單色圖像,因此判別的依據(jù)是粒子軌跡的亮度要高于背景流場,可用的算法包括亮度閾值以及邊緣梯度等經(jīng)過檢驗的分塊算法。之后,需要根據(jù)分塊得到的粒子軌跡圖像來確定粒子軌跡長度。僅根據(jù)圖像所包含的像素數(shù)目,只能得到像素級精度軌跡長度。為進(jìn)一步實現(xiàn)亞像素級精度,需要借助相關(guān)模型,對粒子圖像進(jìn)行擬合。對大多數(shù)實際情況而言,粒子的圖像主要來自于粒子所發(fā)光(或散射光)經(jīng)成像系統(tǒng)衍射后的光斑[22],此類粒子圖像的光強(qiáng)分布近似為高斯分布[23-24],則粒子軌跡的光強(qiáng)分布gs(x′)為[23-24]

(2)

其中,α為灰度值系數(shù);x′為以軌跡灰度權(quán)重中心為原點的粒子軌跡坐標(biāo);n為沿粒子運(yùn)動方向的單位矢量;σ為粒子散射光強(qiáng)分布參數(shù);I0為與圖像背景亮度有關(guān)的常數(shù);G(x)為高斯函數(shù):

(3)

根據(jù)式(2)對粒子軌跡圖像的光強(qiáng)值進(jìn)行擬合,可以得到亞像素精度粒子軌跡長度。本文中稱此方法(或與之類似的方法,如采用不同的光強(qiáng)分布函數(shù))為二維擬合(2D fit)算法。

二維擬合方法需要一次同時對α、lp、σ、I0及n共6個參數(shù)(n含有2個未知量)進(jìn)行擬合,尤其是方向n的確定與其他參數(shù)耦合,通常需要多次迭代求解,計算速度較低。為此,提出了一種新的,同樣具有亞像素精度,但計算較為簡便的粒子軌跡長度算法。新算法的計算步驟如下:

(1)圖像分塊,得到單個粒子軌跡(詳見1.2節(jié)中介紹的粒子軌跡識別程序)。

(2)得到每個粒子軌跡所在的直線方程。如果記粒子軌跡的直線方程為

y=kx+b

(4)

則此步目的是得到參數(shù)k與b,該參數(shù)組合通過對每個粒子軌跡采用帶亮度值權(quán)重的最小二乘法進(jìn)行直線擬合得到。為此,定義擬合誤差函數(shù)Eerror為

(5)

式中wi為該粒子軌跡圖像中像素亮度值。

根據(jù)最小二乘法可以求得式(4)中的k與b的值分別為

b=Sy-kSx

(6)

其中,

(7)

式(7)中的求和是針對這個粒子軌跡圖像所含的所有像素。

(3)由于尾焰中粒子速度相對于軸線方向偏離很小,圖像上粒子軌跡基本都為水平方向,因此可以沿豎直方向?qū)壽E圖像降維,即計算粒子軌跡每列像素的亮度值之和I1D(x),其中x為列位置,如圖2所示。

圖2 亞像素精度軌跡長度一維擬合算法示意圖Fig.2 Schematic diagram of sub-pixel accuracy particle streak length of 1D fitting algorithm

(4)對一維亮度函數(shù)I1D(x)進(jìn)行擬合,得到粒子軌跡長度。將式(2)給出的粒子軌跡圖像模型沿豎直方向積分可得一維亮度函數(shù)的模型函數(shù):

(8)

其中,c1D、a1D、I0為常數(shù);x0為粒子軌跡起點在x軸上的坐標(biāo);lp_x為粒子軌跡在x軸方向上的投影長度;erf(x)為誤差函數(shù),定義為

(9)

利用式(8)對I1D(x)進(jìn)行擬合可得x0及l(fā)p_x。

(5)根據(jù)式(4),確定粒子軌跡長度lp為

(10)

如果需要,可以更近一步確定粒子軌跡的起點坐標(biāo)(x0,kx0+b)以及終點坐標(biāo)(x0+lp_x,k(x0+lp_x)+b)。將該方法被稱為一維擬合(1D fit)算法。相比于原二維擬合算法,此方法需擬合5個參數(shù)c1D、a1D、x0、lp_x及I0,雖然看起來僅減少了1個擬合參數(shù),但是剩余5個參數(shù)間耦合程度低,擬合收斂快。為更進(jìn)一步節(jié)省計算時間,對I1D(x)擬合時可以將式(8)中的誤差函數(shù)erf(x)用指數(shù)函數(shù)來近似,即用式(11)來代替式(8)。

(11)

將I1D(x)進(jìn)行歸一化后,使用式(8)與式(11)對某個顆粒軌跡的I1D進(jìn)行擬合,結(jié)果如圖3所示。從圖3中可以發(fā)現(xiàn)式(8)與式(11)都可以很好地對I1D進(jìn)行擬合。但采用式(11)進(jìn)行擬合的優(yōu)點是擬合算法更為簡潔,無需計算誤差函數(shù)erf(x),進(jìn)一步減少計算量。

圖3 兩種形式擬合結(jié)果Fig.3 Fitting results of two equations

1.2.2 含有模擬粒子軌跡的合成圖像的生成

為分析上述兩種粒子軌跡長度算法(二維及一維擬合算法)的精度及計算效率,生成了各個參數(shù)已知的模擬粒子軌跡的計算合成圖像(synthetic image),利用合成圖像來檢驗上述算法。為生成模擬粒子軌跡圖像,考慮一個自身發(fā)光或散射光的粒子在相機(jī)成像面對應(yīng)的光強(qiáng)分布函數(shù)P(x,y)為

P(x,y)=f(x-xc,y-yc)

(12)

式中 (xc,yc)為粒子中心在相機(jī)成像平面的坐標(biāo)。

理想情況下,球形顆粒在成像平面的光強(qiáng)分布函數(shù)P(x,y)為艾瑞斑(Airy disk),即

(13)

其中,A為系數(shù);J1(r)為第一類貝塞爾函數(shù);r為成像面上任意一點到粒子中心的距離:

(14)

對于實際成像系統(tǒng),粒子散射光經(jīng)過光學(xué)系統(tǒng)中一系列光學(xué)器件后其在像平面的光強(qiáng)分布近似為二維高斯函數(shù)[25],即

(15)

通常可以假設(shè)σx=σy=σ。

假設(shè)在相機(jī)曝光時間Δt內(nèi)粒子以速度vxi+vyj(i與j是單位方向矢量)做勻速直線運(yùn)動,則將光強(qiáng)分布對相機(jī)曝光時間Δt及像素區(qū)域積分后得到圖像上任意像素(x1,y1)處的灰度值I(x1,y1)為

(16)

式中β為與成像系統(tǒng)有關(guān)的能量-灰度值轉(zhuǎn)換系數(shù)。

在實際測量過程中,采集到的圖像不可避免地會含有由于背景雜光引起的圖像噪音。為此在生成的粒子軌跡模擬圖像上添加了強(qiáng)度為γnoise的高斯白噪音,即圖像上每一像素的實際亮度值為

(17)

式中 floor為取整函數(shù);γnoise為代表噪音強(qiáng)度的參數(shù);N(0,1)為服從均值為0,方差為1的高斯分布的隨機(jī)變量;H代表圖像位深,H=12。

計算生成的合成圖像大小為1280 pixel × 384 pixel,每幅圖像中有20個模擬粒子軌跡,軌跡的初始位置在圖像中隨機(jī)分布。為檢驗粒子軌跡擬合算法在不同工況下的精度,改變了模擬粒子軌跡的參數(shù),包括:軌跡寬度σ、粒子軌跡長度lp、噪音強(qiáng)度γnoise,以及粒子軌跡方向Γ=vy/vx。各參數(shù)取值見表1。

表1 模擬粒子軌跡圖像參數(shù)Table 1 Parameters of the generated particle streaks

對應(yīng)于表1中每種參數(shù)組合,生成了4000幅圖像,即每種模擬工況下有80 000條模擬粒子軌跡。同時,為檢驗光強(qiáng)分布函數(shù)的影響,在每種模擬工況下,分別生成了光強(qiáng)分布為艾瑞斑與高斯分布的粒子軌跡。為保證可比性,假設(shè)在給定工況下無論粒子圖像光強(qiáng)為何種分布,其所含輻射能量相同。由此易得在對應(yīng)工況下高斯分布的參數(shù)σ與艾瑞斑參數(shù)A之間滿足以下關(guān)系:

Aσ≈2

(18)

1.2.3 亞像素精度粒子軌跡長度算法精度

將擬合得到的粒子軌跡長度與合成圖像時輸入的模擬參數(shù)進(jìn)行比較,可以得到擬合算法的粒子軌跡長度計算誤差。結(jié)果表明,任意給定工況下,二維與一維擬合算法得到的粒子軌跡誤差均值都為零。為定量評價誤差的散布并避免極少數(shù)離群值的影響,定義粒子軌跡長度誤差散布Δlp為[24, 26]

(19)

式中φ(x)為粒子軌跡長度誤差的概率密度分布;Δlp為擬合得到的粒子軌跡長度80%的概率在真值附近[-Δlp/2,Δlp/2]的區(qū)間內(nèi)。

不同工況下統(tǒng)計得到的誤差散布Δlp如圖4所示。

(a)γnoise=64,Γ=0.2,lp=20 pixel (b)γnoise=64,σ=0.6 pixel,lp=20 pixel

(c)γnoise=64,σ=0.6 pixel,Γ=0.2 (d)σ=0.6 pixel,Γ=0.2,lp=20 pixel圖4 不同工況下算法精度Fig.4 Algorithm accuracy for different conditions

從圖4(a)可以發(fā)現(xiàn),對于大多數(shù)測試的粒子軌跡寬度,二維擬合算法精度要高于一維擬合精度,這是符合預(yù)期的。然而,對于一維擬合算法,隨著圖像寬度增加擬合精度略有提高,對于二維擬合算法,隨著粒子軌跡寬度增加擬合精度略有降低。

尤為特別的是,當(dāng)粒子軌跡寬度σ<0.4 pixel時,不論是粒子圖像光強(qiáng)分布是高斯分布還是艾瑞斑分布,一維擬合算法精度都比二維擬合的精度高。這是因為當(dāng)σ<0.4 pixel時(算例中σ=0.2 pixel)散射能量在垂直于粒子運(yùn)動方向上集中在一個窄于一個像素的寬度內(nèi),由于圖像是以像素為單位離散化的,因此會出現(xiàn)沿粒子運(yùn)動方向,同一個粒子軌跡在不同位置的光強(qiáng)分布在一個或兩個像素的情況,而相關(guān)像素的亮度值會劇烈變化,這導(dǎo)致二維擬合算法精度降低。對于一維擬合算法,由于將粒子軌跡所涉及像素的亮度在水平方向(接近粒子運(yùn)動方向)加和,相當(dāng)于將粒子圖像所含的能量在垂直于運(yùn)動方向進(jìn)行了積分,這消除了粒子圖像在垂直于運(yùn)動方向分布于幾個像素的影響,所以一維擬合精度幾乎不受影響,故當(dāng)粒子軌跡寬度較小時,一維擬合算法的精度高于二維擬合算法精度。

此外,粒子圖像的光強(qiáng)分布形式對一維擬合算法精度有一定影響:在相同粒子軌跡寬度下,一維擬合算法對粒子光強(qiáng)分布為高斯分布的軌跡長度計算精度略高于光強(qiáng)分布為艾瑞斑時的計算精度。這應(yīng)該是由于擬合所用的算法,即式(8)或其近似式(11),是基于高斯光強(qiáng)分布導(dǎo)出的。二維擬合算法的精度則對粒子圖像的光強(qiáng)分布不敏感,對于高斯分布和艾瑞斑分布得到的結(jié)果幾乎一樣。這是由于雖然二維擬合時也采用了粒子光強(qiáng)分布為高斯分布的假設(shè),但是二維擬合時更多利用到圖像的空間分布信息。從圖4(b)～(d)可以發(fā)現(xiàn),這一規(guī)律對其他參數(shù)變化時也同樣成立,因此后文中不再贅述。值得指出的是,對于實際的光學(xué)成像系統(tǒng),粒子圖像的光強(qiáng)分布更接近于高斯分布,因此采用式(8)或式(11)進(jìn)行擬合是合理的。

圖4(b)表明,粒子運(yùn)動方向無論對一維擬合算法還是二維擬合算法的精度都只有非常微弱的影響。

粒子軌跡長度對算法精度的影響體現(xiàn)在圖4(c)中:一維和二維擬合算法的精度都先隨著粒子軌跡長度的增加而增加,但當(dāng)粒子軌跡長度大于大約20 pixel 后,軌跡長度進(jìn)一步增加對計算精度就不再有影響,無論是采用一維還是二維擬合都是如此。這個結(jié)果對于實際實驗測量具有重要的指導(dǎo)意義,可用于調(diào)整相機(jī)曝光時間以獲得最佳的粒子軌跡長度。因為如果曝光時間過長形成過長的粒子軌跡的話,出現(xiàn)粒子軌跡交叉的可能性也增加,這會給圖像處理帶來額外的困難。實際上,圖4(c)的結(jié)果表明,如果實際測量時粒子過濃,為避免粒子軌跡交叉,可以考慮通過縮短曝光時間將粒子軌跡長度降至10 pixel左右,這只帶來軌跡長度測量精度的略微降低。

圖4(d)給出了圖像噪音水平對軌跡長度測量精度的影響:兩種算法的誤差散布都隨著圖像噪音水平增大而幾乎線性增大。注意圖中橫軸的噪音強(qiáng)度為背景高斯噪音的標(biāo)準(zhǔn)差(standard deviation),對于生成的12-bit合成圖像,噪音強(qiáng)度為128(即7-bit)已是非常強(qiáng)烈的背景噪聲了。即便是這種情況下,兩種擬合算法的精度依然是亞像素級。實際情況下,通過對原始圖像進(jìn)行去噪濾波,可以進(jìn)一步降低背景噪音的影響。

通過以上分析發(fā)現(xiàn),通常情況下,二維擬合算法的精度都略優(yōu)于一維擬合算法,但只要粒子軌跡長度不過短,如不低于20 pixel,一維擬合算法計算得到單個粒子軌跡長度的相對誤差即使在背景噪音很強(qiáng)烈的情況下也不超過±2%。尤其值得指出的是,這是指單個粒子軌跡長度的測量誤差,對于統(tǒng)計量,如粒子的平均速度,因為算法的無偏性,其誤差更小。

除了分析不同工況下兩種算法的精度,還比較了兩種方法計算軌跡長度所用時間,統(tǒng)計而言,對同一條粒子軌跡,二維擬合算法所用時間約是一維擬合算法的8倍。綜合考慮算法精度和算法計算效率,在后續(xù)實驗測量中采用一維擬合算法進(jìn)行粒子軌跡長度計算。

1.3 固體火箭發(fā)動機(jī)尾焰中顆粒速度測量

利用前文介紹的粒子軌跡法,對一真實固體火箭發(fā)動機(jī)尾焰中顆粒的速度進(jìn)行了實驗測量。在本小節(jié)中,對實驗裝置布置和所用計算程序進(jìn)行簡要介紹。具體實驗結(jié)果及討論則在第2節(jié)中詳細(xì)敘述。

1.3.1 實驗裝置

固體火箭發(fā)動機(jī)尾焰中較大的顆粒其亮度通常比背景亮度高,因此,實驗時采用對尾焰直接攝影的方法,無需設(shè)置額外的照明光源,如圖5所示。

圖5 固體火箭發(fā)動機(jī)尾焰中顆粒速度測量原理示意圖Fig.5 Schematic diagram of particle velocity measurement in solid rocket motor

實驗時,采用高速相機(jī)正對地面熱試車中的固體火箭發(fā)動機(jī)尾焰區(qū)拍攝。本文中使用的高速相機(jī)型號為Phantom 710L單色相機(jī),鏡頭焦距105 mm,圖像幅面設(shè)置為1280 pixel×384 pixel。在實驗開始前先標(biāo)定得到的成像系統(tǒng)空間分辨率Rs=0.37 mm/pixel,標(biāo)定時采用高速相機(jī)拍攝表面加工有規(guī)則圓斑陣列的標(biāo)定板,對拍攝得到的圖像處理后得到標(biāo)定圖像中兩圓斑的圖像間距(單位為像素),則圓斑間物理間距除以對應(yīng)圖像間距即為空間分辨率Rs,經(jīng)估算標(biāo)定誤差不超過1%。高速相機(jī)距離發(fā)動機(jī)尾焰大約1.5 m。通過對尾焰中粒子速度的預(yù)估,結(jié)合標(biāo)定得到的成像系統(tǒng)空間分辨率,設(shè)置高速相機(jī)曝光時間Δt為9.41 μs,使得圖像上粒子軌跡長度大約為25 pixel。相機(jī)拍攝幀速為15 000 fps,接近該幅面下相機(jī)所能支持的最高拍攝速度。高速相機(jī)的拍攝開始時刻由外部觸發(fā)信號控制,該外部觸發(fā)信號可以是火箭發(fā)動機(jī)的點火信號或者由實驗者設(shè)置的其他信號。

1.3.2 發(fā)動機(jī)尾焰圖像中粒子軌跡識別

利用高速相機(jī)拍攝得到發(fā)動機(jī)尾焰圖像后,需要識別出尾焰圖像中的粒子軌跡。由于尾焰中存在細(xì)小的顆粒云團(tuán),這使得尾焰圖像背景噪音較大,粒子軌跡識別程序需要從背景噪音中分辨出單個的粒子軌跡。性能良好的粒子軌跡識別程序是進(jìn)行尾焰中顆粒速度測量的重要步驟。本文用MATLAB語言編寫了對數(shù)字圖像的粒子軌跡識別程序,其主要流程如圖6所示。

圖6 圖像中粒子軌跡識別主要過程Fig.6 Main process of particle streak recognition in image

在讀入圖像后首先對圖像去噪,本文中使用均值濾波去噪音。圖像去噪后利用圖像灰度值梯度閾值進(jìn)行粒子軌跡邊緣識別得到二值圖(binary image),找到二值圖中的8鄰域連通區(qū),并去掉小的連通區(qū),因為這些小的區(qū)域?qū)?yīng)的是圖像中殘余噪聲。之后對通過篩選的二值圖像進(jìn)行圖像膨脹(dilation),圖像膨脹用的卷積核為3×3。圖像膨脹的目的是進(jìn)行圖像填充,以補(bǔ)充軌跡中間灰度梯度較小而未能被梯度閾值法識別部分,但圖像膨脹也使得粒子軌跡圖像向四周擴(kuò)展,因此再用3×3的卷積核對圖像進(jìn)行腐蝕(erosion),以恢復(fù)粒子軌跡的大小[27]。

經(jīng)過上述操作后可找到疑似粒子軌跡,之后以疑似粒子軌跡為中心,劃分較小的粒子軌跡檢測區(qū),進(jìn)行第二次粒子軌跡識別。對這些小的粒子軌跡檢測區(qū)利用局部直方圖法進(jìn)行局部圖像增強(qiáng)[28],對增強(qiáng)后的檢測區(qū)再進(jìn)行圖像膨脹、腐蝕、邊緣檢測等一系列處理,其中邊緣檢測使用Sobel算法。局部圖像增強(qiáng)一方面可以避免將一個顆粒軌跡誤識別為多個軌跡,另一方面可以提高顆粒軌跡起始點位置識別精度。此處得到的軌跡起始點位置作為亞像素精度軌跡長度算法擬合的初值,較高精度的擬合初值可以減少擬合計算所用時間。

上述處理后,圖像中的單個噪點不會影響軌跡識別結(jié)果,但亮度變化劇烈的區(qū)域,如較亮的顆粒云團(tuán)邊緣仍可能被錯誤識別為粒子軌跡。為此引入粒子軌跡的判別條件以去除誤識別的軌跡。粒子軌跡判別條件包括:(1)顆粒的軌跡應(yīng)該是直線狀;(2)單幅圖像上粒子軌跡長度(即在給定的曝光時間內(nèi)顆粒飛行的距離)不可能超過某一上限(該上限可以通過對顆粒速度最大值的估計來確定);(3)顆粒速度矢量指向不應(yīng)違反物理規(guī)律,如速度矢量應(yīng)該遠(yuǎn)離中心線,而非指向中心線。利用此判別準(zhǔn)則得到最終粒子軌跡,記錄粒子軌跡長度、端點位置、出現(xiàn)時刻等,用于后續(xù)顆粒速度統(tǒng)計分析。

經(jīng)過以上圖像處理過程,得到了圖像中代表單個粒子軌跡圖像的一組連續(xù)的像素,對此應(yīng)用前文介紹的一維擬合算法得到具有亞像素精度的粒子軌跡長度。

2 結(jié)果與討論

在前文介紹的對某固體火箭發(fā)動機(jī)測試過程的現(xiàn)場測量中,共得到25 910幅圖像。對其利用本文所介紹的粒子識別程序?qū)D像中的粒子軌跡進(jìn)行識別,共識別出粒子軌跡約400 000個。圖7(a)所示為某時刻的該固體火箭發(fā)動機(jī)尾焰圖像及識別出的尾焰中粒子軌跡,圖中黃色框表示本文中粒子識別程序處理的圖像范圍,紫色箭頭從程序識別出的粒子軌跡起點指向粒子軌跡終點,圖7(b)為粒子軌跡局部放大圖,圖中粒子軌跡寬度約為3 pixel,其他粒子軌跡圖像多與之相似?？紤]粒子成像系統(tǒng)的放大比、光圈大小以及衍射效應(yīng),根據(jù)文獻(xiàn)[29]可以估計粒子直徑約為100～200 μm。

識別單個粒子軌跡,并利用前文介紹的一維擬合算法對粒子軌跡長度進(jìn)行計算,根據(jù)式(1)得到粒子瞬時速度。圖8中給出了在拍攝時間段內(nèi),距離噴管82、228、430 mm處粒子軸向平均速度隨時間的變化。該平均速度為在相應(yīng)位置處每300幀圖像(對應(yīng)的時間為 0.02 s)中所有粒子速度的平均值。如前所述,利用粒子軌跡法計算顆粒速度時,對于顆粒軌跡長度平均值,粒子軌跡長度算法是無偏的。實驗結(jié)果表明,在拍攝時間段內(nèi)三個不同位置處的顆粒平均速度隨時間沒有明顯變化。因此接下來對得到的粒子速度在整個拍攝時間段內(nèi)平均,以計算尾焰流場內(nèi)不同位置處顆粒的平均速度。

(a)Recognition results of particle streaks in SRM plume for one image (b)Enlarged view of particle streaks圖7 尾焰流場及粒子軌跡局部放大圖Fig.7 Image of plume field and enlarged view of particle streaks

圖8 尾焰中3個不同位置處顆粒平均速度隨時間變化關(guān)系Fig.8 Mean particles velocityin three different locations varies with time

圖9所示為尾焰不同位置處顆粒軸向平均速度。從圖9可以看出尾焰中的粒子速度在靠近噴管出口處最大,超過1000 m/s,隨著距噴口距離的增加尾焰中的粒子速度降低,但粒子速度隨著距離衰減特性與氣相湍流射流的速度衰減特性不同:大尺寸粒子速度衰減弱于氣相湍流射流的速度衰減,這是因為大粒子的斯托克數(shù)較大造成的。

圖9 平均顆粒軸向速度隨距離噴口位置的變化Fig.9 Mean particles axial velocity varies with distance from the nozzle

被測發(fā)動機(jī)尾焰為超音速氣流,根據(jù)發(fā)動機(jī)設(shè)計參數(shù)仿真計算表明尾焰中存在馬赫盤。為觀察是否有大尺寸顆粒穿過激波后發(fā)生速度突變,圖10中給出了二維流動平面內(nèi)粒子速度幅值分布的等值線圖(圖中左下角處白色區(qū)域是因為在該處沒有拍攝到粒子軌跡,故該處沒有粒子速度數(shù)據(jù))。

圖10 全流場中顆粒平均速度分布Fig.10 Mean particles velocity distribution in full plume field

與前面軸向平均速度類似,從圖10中可以發(fā)現(xiàn):靠近噴管軸線位置處粒子速度最高,隨著遠(yuǎn)離軸線,粒子速度降低。因為穩(wěn)態(tài)射流流場中激波位置是固定的,如果粒子速度受到流場中激波明顯影響的話,從該圖中應(yīng)該可以觀察到在空間中某些位置出現(xiàn)粒子速度劇烈變化的區(qū)域。然而圖中并不存在類似區(qū)域,表明尾焰中激波對于觀察到的大尺寸粒子的速度沒有顯著影響。

3 結(jié)論

本文采用合成圖像對粒子軌跡測速法中的粒子軌跡長度擬合算法進(jìn)行了仿真分析,并開發(fā)了完整的粒子軌跡測速法程序,對某固體火箭發(fā)動機(jī)尾焰中大尺寸顆粒速度進(jìn)行了實驗研究,得到以下結(jié)論:

(1)模擬顆粒軌跡計算結(jié)果表明,本文提出的一維擬合粒子軌跡亞像素精度長度算法計算效率高,單個粒子軌跡長度計算的不確定性都在亞像素精度。為減小單個粒子軌跡長度的相對誤差,最有效的方法是使得實驗圖像中粒子軌跡的長度約為20 pixel,至少不短于10 pixel。此外,擬合算法得到的粒子軌跡長度平均值是無偏的。

(2)該發(fā)動機(jī)尾焰流場中顆粒速度測量誤差約為3%。粒子軌跡法可以實現(xiàn)尾焰流場中大尺寸顆粒速度測量,該技術(shù)為尾焰中大尺寸顆粒的速度測量提供了新的思路。

(3)測量結(jié)果表明,固體火箭發(fā)動機(jī)尾焰中的大尺寸顆粒速度隨著距噴管中心距離增加而減小,也隨距噴管出口的軸向距離增加而衰減,但是衰減速度小于氣相湍流射流中的速度衰減。結(jié)果中未觀察到尾焰中激波對大尺寸顆粒速度產(chǎn)生影響。

致謝：感謝基礎(chǔ)科研項(JCKY2019204B020)的支持。