某型液體火箭發(fā)動機噴注參數(shù)與尾焰紅外圖像特征研究

申帥帥，張 楠，王大銳

0 引 言

發(fā)動機尾焰紅外圖像的研究主要在以下兩個方面：a）利用尾焰紅外輻射能量大、紅外圖像突出的特點，來識別目標飛行器；b）對尾焰紅外圖像本身特征的研究。胡炳梁等[1]成功研制了中國第1個超聲速機載紅外測量吊艙系統(tǒng)，為測量目標機尾焰紅外特性提供了有效的測量手段；史麗芳等[2]完成了飛機尾焰紅外圖像識別軟件系統(tǒng)，該系統(tǒng)利用低空衛(wèi)星攜帶的紅外熱像儀來探測目標，并根據(jù)各個階段尾焰呈現(xiàn)的不同特性（形體特征、溫度）來檢測與識別來襲飛機，及時預警并主動攻擊，具有較高的辨識率和較好的發(fā)展前景；張碩等[3]開發(fā)了紅外圖像識別系統(tǒng)，將小波分析、中值濾波以及數(shù)學形態(tài)學的方法、密度分割法和頻率閾法應用于固體火箭發(fā)動機排氣羽焰圖像的識別處理中，得到了較好的火焰輪廓及其基本參數(shù)，如火焰量值、長寬比、均值對比度、亮度最大像素點、復雜度等參數(shù)等，具有較好的識別效果。

尾焰紅外圖像中包含豐富的燃燒信息，較好地反映發(fā)動機工作狀態(tài)。但是，目前通過建立尾焰紅外圖像特征與發(fā)動機工作狀態(tài)之間的關系對發(fā)動機工作狀態(tài)進行監(jiān)測診斷的研究還比較少。由于室壓、壁溫等常規(guī)監(jiān)測參數(shù)不對發(fā)動機燃燒過程中所有狀態(tài)變化都敏感，而尾焰流場的亮度、特定區(qū)域的形狀等特征參數(shù)對于不同的燃燒工況反應明顯[4]，且尾焰紅外圖像的測量方法具有非接觸、無干擾、精度高等特點，因此對發(fā)動機尾焰光學測量監(jiān)測的應用研究越來越廣泛。

本文采用多學科仿真技術，耦合蒸發(fā)、燃燒、傳熱過程，建立了三維非穩(wěn)態(tài)推力室蒸發(fā)-燃燒-尾焰場一體化仿真模型；采用尾焰圖像識別方法，開展液體火箭發(fā)動機噴注參數(shù)與尾焰紅外圖像特征之間關系的研究。

1 物理模型和計算方法

1.1 物理模型及邊界條件

本文仿真對象為某型液體火箭發(fā)動機，燃料組合為甲基肼（MMH）和綠色四氧化二氮（NTO）。該仿真模型采用三維非穩(wěn)態(tài)仿真模型，耦合了蒸發(fā)、燃燒、傳熱等多個過程。

三維仿真計算模型網(wǎng)格劃分如圖 1所示，整個計算域采用多層O形網(wǎng)格劃分，壁面附近劃分邊界層網(wǎng)格，網(wǎng)格精度較高。計算域由流體區(qū)域Fluid塊和固體區(qū)域Solid塊組成，內(nèi)壁面采用Coupled壁面條件，外壁面采用固定熱流壁面條件，尾焰場區(qū)域邊界采用Pressure-Out邊界條件。

圖1 三維仿真計算模型示意Fig.1 Schematic Diagram of Three-dimensional Simulation Model

1.2 計算模型

該推力室-尾焰場仿真模型的計算采用Euler-Lagrange方法，控制方程包括 Euler坐標系下的N-S方程和 Lagrange方程坐標系下的離散顆粒軌道方程。氣液兩相之間的耦合通過源項實現(xiàn)。將液滴簡化為質點，用 Lagrange法進行追蹤，液滴狀態(tài)由其中心位置XP、直徑DP、速度UP和密度ρp4個物理量確定；對于氣相流體采用Euler坐標系下的N-S方程進行描述。

為提高計算精度，在計算允許的條件下，應盡量增加噴注器數(shù)量，本次仿真中布置了 352個噴注器，共計1958束液滴流束。噴注器布置面垂直于推力室軸線方向，其示意如圖2所示。

圖2 噴注器布置示意Fig.2 Schematic Diagram of Injector Layout

各個噴注器的位置、液滴噴射速度、粒徑、流量均由相位多普勒粒子動態(tài)分析儀（Phase Doppler Particle Analyzer，PDPA）霧化試驗測量獲得。由于仿真過程中采用的工質為MMH和NTO，而霧化試驗采用的工質是水，所以布置噴注器時，液滴的噴注參數(shù)需要進行換算。

蒸發(fā)模型中采用的顆粒類型為液滴，遵循加熱、冷卻、蒸發(fā)、沸騰規(guī)律，當液滴剩余質量大于揮發(fā)極限質量，且溫度分別達到汽化和沸騰溫度時，液滴便會以相應的速率發(fā)生蒸發(fā)和沸騰，液滴表面質量輸運方程為

式中 mp(t)為 t時刻液滴質量；mp(t+Δt)為 t+Δt時刻液滴質量；Ni為蒸發(fā)速率系數(shù)，與飽和蒸汽壓、蒸汽擴散系數(shù)、運動粘度系數(shù)有關；Ap(t)為t時刻液滴表面積。

化學反應采用MMH和NTO一步總包反應，化學反應模型為Finite-Rate/Eddy-Dissipation燃燒模型[5]，假定化學反應速率由化學動力學模型和渦耗散模型共同控制，取兩模型各自計算反應速率的較小值作為化學反應速率。其中，化學動力學控制的反應速率采用阿累尼烏斯（Arrhenius）公式計算，渦耗散模型假定化學反應速率由渦團破碎的速率計算。

MMH和NTO液體推進劑的物性參數(shù)（密度、比熱容、潛熱、飽和蒸汽壓等）采用《液體推進劑》[6]提供的數(shù)據(jù)；壁面材料的鈮鎢合金物性參數(shù)（熱導率、比熱容等）采用《航天用鈮鎢合金棒材規(guī)范》[7]提供的數(shù)據(jù)；各氣相物質的物性參數(shù)，如熱導率、粘度等為常數(shù)，定壓比熱容為關于溫度的多項式。

1.3 尾焰紅外圖像識別方法

本文采用液體火箭發(fā)動機尾焰紅外圖像特征識別方法[4]作為尾焰場特征識別的工具，該方法具體流程如圖3所示。

基本步驟如下：

a）對圖像進行灰度化處理，并設定合適閾值獲取需要的灰度圖像，過濾掉部分干擾點后得到邊緣點圖像；

b）將邊緣點連接成為線段，過濾掉干擾線段，并為線段排序編號；

c）將各邊緣線段和構造的區(qū)域邊界連接起來，構成指定區(qū)域（尾焰出口區(qū)、尾焰射流區(qū)、尾焰核心區(qū)）；

d）計算各區(qū)域周長、面積、灰度熵，第1馬赫盤位置等特征參數(shù)值。

圖3 尾焰紅外圖像特征識別方法流程Fig.3 Schematic Diagram of the Process of Infrared Flame Image Identification

通過紅外圖像特征識別方法識別出2個區(qū)域，如圖4所示，區(qū)域1為尾焰出口區(qū)，區(qū)域2為尾焰射流區(qū)。

圖4 仿真尾焰紅外圖像識別分區(qū)示意Fig.4 Schematic Diagram of Simulation of Tail Flame Infrared Image Recognition

將2個區(qū)域的周長、面積、灰度熵、第1馬赫盤位置作為尾焰場紅外圖像的特征指標。

假定噴管出口直徑為D ，圖像識別出其像素個數(shù)為n，則單個像素尺寸為D/n。

a）區(qū)域周長Li：識別區(qū)域的邊緣像素點個數(shù)與單個像素尺寸的乘積：

式中 n1i為橫或者豎方向線段數(shù)量；n2i為對角線方向的線段數(shù)；下標i為識別區(qū)域代號，i=1，2。

b）區(qū)域面積Si：識別區(qū)域的像素點個數(shù)與單個像素尺寸的平方的乘積：

式中 Ni為區(qū)域內(nèi)的像素點數(shù)，下標i為識別區(qū)域代號，i=1，2。

c）灰度熵：灰度熵是圖像分割中常用的特征參數(shù)，它表示圖形中像素灰度的不均勻程度或者復雜程度?；叶褥卦叫。砻鲌D像均勻性越好[8]。

對于一幅256級灰度圖像而言，每個像素的灰度值可能是0～255級灰度中的某一級，因此像素灰度值的符號集為{0，1，2，…，255}。假設圖像中灰度值為j的像素個數(shù)為kj，則圖像中灰度值為j的像素出現(xiàn)的概率為

整幅圖像的灰度熵值為

式中 I為灰度。

d）第1馬赫盤位置：如圖5所示，噴管出口射流經(jīng)強烈壓縮形成的第 1道正激波所在的位置即為第 1馬赫盤位置。用第1馬赫盤所在位置與噴口之間的像素點數(shù)與單個像素尺寸的乘積表示，計算表達式如下：

式中 y1為第1馬赫盤所在位置的像素坐標；y0為噴管出口所在位置的像素坐標；y1-y0為第 1馬赫盤距離噴管出口的像素個數(shù)。

圖5 噴管出口附近波系示意Fig.5 Schematic Diagram of the Waves Construction out of the Nozzle

2 仿真和熱試車結果對比分析

2.1 仿真和熱試車室溫、室壓、壁溫、尾焰場對照

多對互擊式發(fā)動機仿真和熱試車外壁面、尾焰場紅外圖像特性如圖6所示。

圖6 仿真和熱試車外壁面、尾焰場紅外特性云圖Fig.6 Simulation and Hot Firing Test of Outer Wall Surface and Tail Flame Infrared Image

續(xù)圖6

由圖6可知：圖6a和圖6b的仿真結果較為一致，仿真尾焰場出現(xiàn)明顯的馬赫盤結構，且第1馬赫盤位置清晰；圖6c和圖6d尾焰輪廓相差較大，分析可知：由于圖6d中的尾焰圖像是發(fā)生復燃反應的發(fā)動機的試車數(shù)據(jù)，所以仿真結果尾焰邊緣未出現(xiàn)外擴現(xiàn)象是正常的，并發(fā)現(xiàn)仿真和熱試車尾焰場溫度波動幅度不一樣，這是因為不同溫度下流場的輻射系數(shù)不同，而圖6d是根據(jù)紅外圖像對全尾焰場取同一輻射系數(shù)求解所得溫度場，導致結果有偏差。

多對互擊式發(fā)動機推力室內(nèi)溫度、壓力、馬赫數(shù)特性如圖7所示。

圖7 推力室中軸面溫度、壓力、馬赫數(shù)仿真云圖Fig.7 Simulation of Axial Plane in the Thrust Chamber Temperature, Pressure and Mach Number Cloud

由于離散化噴嘴布置在噴注面板處，液滴汽化吸熱，使該處溫度較低；推力室直線段壓力變化不大，再經(jīng)過喉部和噴管擴張段，燃氣壓力逐漸減小，馬赫數(shù)逐漸增大，且在喉部附近達到聲速，與試車結果和氣體動力學理論符合較好。仿真與熱試車結果對照如表1所示。

表1 仿真、熱試車結果對照Tab.1 Comparison Table of Simulation and Hot Firing Test

如表1可知，仿真和熱試車的室壓和外壁最高溫度偏差均小于5%，該仿真模型結果較為準確。

仿真和熱試車外壁面溫度沿推力室軸向無量綱位置變化規(guī)律如圖8所示。

圖8 仿真和熱試車外壁面溫度沿軸向變化規(guī)律Fig.8 Figure of Simulation and Hot Firing Test of Outer Wall Surface Temperature Changes Along the Axis

由圖 8可知，在距離噴注面板一定距離內(nèi)（區(qū)域1），由于邊區(qū)液膜、氣膜冷卻作用，壁面溫度較低；距離噴注面板一定距離（區(qū)域2）時，氣膜消失，高溫燃氣直接沖刷內(nèi)壁，壁面溫度迅速升高；隨著距離的增加（區(qū)域3），流動截面積變小，熱流密度增加，壁面溫度繼續(xù)上升（一直持續(xù)到喉部前端）；隨著軸向位置的進一步增加（區(qū)域4），由于流動截面積變小，燃氣內(nèi)能轉換為動能加劇，氣流溫度降低，壁面溫度也隨之降低；過了喉部（區(qū)域5），流動截面積變大，熱流密度變小，而且燃氣內(nèi)能進一步轉化為動能，燃氣溫度迅速降低，壁面溫度隨之急劇降低。

對比分析仿真結果可知：仿真和熱試車外壁面溫度變化曲線相吻合，驗證了推力室蒸發(fā)-燃燒模型的計算精度可用該型發(fā)動機的仿真分析。

2.2 仿真和熱試車尾焰識別參數(shù)對照

為了驗證該仿真模型可以應用上述尾焰紅外圖像特征識別方法進行后續(xù)研究，將仿真和熱試車尾焰紅外圖像識別結果進行對比分析，如表2所示。

表2 額定工況下仿真和熱試車尾焰特征參數(shù)對照Tab.2 Comparison Table of Simulation and Hot Firing Test the Tail Flame Characteristic Parameters

由表2可知，除了灰度熵，尾焰圖像的其他特征，第1馬赫盤位置、區(qū)域周長和面積的偏差均小于10%，由于熱試車采用的是紅外熱像儀拍攝的尾焰紅外圖像，而仿真結果采用的是尾焰場溫度云圖近似代替尾焰紅外圖像，其顏色梯度上有差別，導致灰度熵有所偏差。綜上所述，該尾焰紅外圖像特征識別方法對仿真尾焰紅外圖像的識別結果準確可靠。

3 噴注參數(shù)與尾焰特征關系

采用上述仿真模型和尾焰紅外圖像特征識別方法研究了不同的噴注參數(shù)（液滴噴射速度、液滴平均粒徑、邊區(qū)冷卻流量百分比）條件下尾焰紅外圖像特征的變化規(guī)律。

3.1 液滴平均噴射速度與尾焰特征關系

通過對不同液滴噴射速度下尾焰場紅外圖像進行特征識別，得到不同噴注速度下第1馬赫盤位置和區(qū)域1、2灰度熵隨液滴平均噴注速度變化圖（見圖9）、區(qū)域1、2周長隨液滴平均噴注速度變化圖（見圖10）、區(qū)域1、2面積隨液滴平均噴注速度變化圖（見圖11）。

由圖9～11可知，隨著液滴噴射速度由20 m/s變化到40 m/s，尾焰場紅外圖像中第1馬赫盤位置、尾焰出口區(qū)和尾焰射流區(qū)的灰度熵和面積均出現(xiàn)明顯變化，且出現(xiàn)極值，而尾焰出口區(qū)和尾焰射流區(qū)的周長沒有明顯變化，說明這兩個特征參數(shù)對于液滴噴射速度變化不敏感。

圖9 第1馬赫盤位置和區(qū)域灰度熵隨液滴平均噴注速度變化Fig.9 Figure of the Area of the First Mach Locus and Gray Value Entropy Variation with Droplet Ejection Average Speed

圖10 區(qū)域1、區(qū)域2周長隨液滴平均噴注速度變化Fig.10 Figure of the Zone 1, 2 Perimeter Variation with the Average Jetting Speed

3.2 液滴平均粒徑與尾焰特征關系的研究

通過對不同的液滴平均粒徑下尾焰場紅外圖像進行特征識別，得到不同液滴平均粒徑下第1馬赫盤位置和區(qū)域1、區(qū)域2灰度熵隨液滴平均粒徑變化圖（見圖12）、區(qū)域1、區(qū)域2周長隨液滴平均粒徑變化圖（見圖13）、區(qū)域1、區(qū)域2面積隨液滴平均粒徑變化圖（見圖 14）。

由圖12～14可知，隨著液滴平均粒徑由15 μm變化到35 μm，尾焰場紅外圖像中第1馬赫盤位置、尾焰出口區(qū)和尾焰射流區(qū)的灰度熵和面積均出現(xiàn)明顯變化。在平均粒徑為25 μm附近，尾焰出口區(qū)和尾焰射流區(qū)灰度熵均出現(xiàn)極小值，尾焰出口區(qū)、尾焰射流區(qū)面積和第1馬赫盤位置均出現(xiàn)極大值，而尾焰出口區(qū)和尾焰射流區(qū)的周長沒有明顯變化，說明這兩個特征參數(shù)對于液滴平均粒徑變化不敏感。

圖12 第1馬赫盤位置和區(qū)域灰度熵隨液滴平均粒徑變化Fig.12 Figure of the Area of the First Mach Locus and Gray Value Entropy Variation with Droplet Average Particle Size

圖13 區(qū)域1、區(qū)域2周長隨液滴平均粒徑變化Fig.13 Figure of the Zone 1, 2 Perimeter Variation with roplet Average Particle Size

圖14 區(qū)域1、區(qū)域2面積隨液滴平均粒徑變化Fig.14 Figure of the Zone 1, 2 Area Variation with Droplet Average Particle Size

3.3 邊區(qū)冷卻流量百分比與尾焰特征關系

通過對不同的邊區(qū)冷卻流量百分比下尾焰場紅外圖像進行特征識別，得到不同邊區(qū)冷卻流量百分比下第1馬赫盤位置和區(qū)域1、區(qū)域2灰度熵隨邊區(qū)冷卻流量百分比變化圖（見圖15）、區(qū)域1、區(qū)域2周長隨邊區(qū)冷卻流量百分比變化圖（見圖16）、區(qū)域1、區(qū)域2面積隨邊區(qū)冷卻流量百分比變化圖（見圖17）。

圖15 第1馬赫盤位置和區(qū)域灰度熵隨邊區(qū)冷卻流量百分比變化Fig.15 Figure of the Area of the First Mach Locus and Gray Value Entropy Variation with Side to Cool the Flow Percentage

圖16 區(qū)域1、區(qū)域2周長隨邊區(qū)冷卻流量百分比變化Fig.16 Figure of the Zone 1, 2 Perimeter Variation with Side to Cool the Flow Percentage

圖17 區(qū)域1、區(qū)域2面積隨邊區(qū)冷卻流量百分比變化Fig.17 Figure of the Zone 1, 2 Area Variation with Side to Cool the Flow Percentage

由圖15～17可知，隨著邊區(qū)冷卻流量百分比由22%變化到30%，尾焰場紅外圖像中第1馬赫盤位置、尾焰出口區(qū)和尾焰射流區(qū)的灰度熵和面積均出現(xiàn)明顯變化，且第1馬赫盤位置、尾焰出口區(qū)和尾焰射流區(qū)的灰度熵出現(xiàn)雙峰變化規(guī)律，尾焰出口區(qū)和尾焰射流區(qū)的面積變化明顯，但沒有發(fā)現(xiàn)特定的變化規(guī)律；尾焰出口區(qū)和尾焰射流區(qū)的周長沒有明顯變化，說明這兩個特征參數(shù)對于邊區(qū)冷卻流量百分比變化不敏感。

觀察還發(fā)現(xiàn)，隨著這3個噴注參數(shù)的變化，區(qū)域灰度熵和第1馬赫盤距離噴管出口的位置均呈現(xiàn)出很強的負相關性。分析知，隨著燃燒效率增大，噴管出口馬赫數(shù)增大，同時，燃燒更充分，尾焰紅外圖像灰度熵減小。說明噴管出口馬赫數(shù)與第1馬赫盤距離噴管出口的距離成正相關，與文獻[9]中結論一致。

4 結 論

采用蒸發(fā)、燃燒、傳熱多學科算法和尾焰場紅外圖像識別方法，利用霧化試驗液滴分布規(guī)律，建立的三維非穩(wěn)態(tài)推力室蒸發(fā)-燃燒-尾焰場一體化仿真模型，對某型發(fā)動機推力室噴嘴霧化特性和邊區(qū)冷卻流量與尾焰場紅外圖像特征的關系仿真分析結論如下：

a）本文建立的模型與熱試車結果吻合較好，具有較高的計算準確性；

b）尾焰射流場的周長對噴注參數(shù)不敏感，而灰度熵、馬赫盤位置和圖像面積的敏感性較強，能夠用來監(jiān)測噴注特性；

c）采用馬赫盤位置、尾焰出口區(qū)和尾焰射流區(qū)的灰度熵和面積峰值可以評估推力室噴嘴霧化參數(shù)和邊區(qū)冷卻流量的選取。

[1] 胡炳梁, 劉學斌，等. 超音速機載吊艙紅外測量系統(tǒng)研究[J]. 紅外與激光工程, 2003, 32(5)： 441-444.

Hu Bingliang, Liu Xuebin, et al. Research on an ultrasonic pod infrared measurement system[J]. Infrared and Laser Engineering, 2003, 32(5)：441-444.

[2] 史麗芳. 基于圖像處理的飛機尾焰檢測與識別技術的研究[D]. 西安：西安電子科技大學, 2015.

Shi Lifang. A research of detection and recognition technology of the plane tail flame based on image processing[D]. Xi’an： Xidian University,2015.

[3] 張碩, 王寧飛, 張平. SRM羽焰圖像的檢測與參數(shù)預估[J]. 導彈與航天運載技術, 2008 (4)： 47-50, 55.

Zhang Shuo, Wang Ningfei, Zhang Ping. Image detection and parameter estimation on flame images of solid rocket motor[J]. Missiles and Space Vehicles, 2008(4)： 47-50, 55.

[4] 王大銳, 張楠, 葛明和. 液體火箭發(fā)動機噴管燃氣外流場紅外圖像研究[J]. 導彈與航天運載技術, 2016(2)： 26-30.

Wang Darui, Zhang Nan, Ge Minghe. Infrared image study on gas external flow field of liquid rocket engine nozzle[J]. Missiles and Space Vehicles,2016(2)： 26-30.

[5] 劉濤. 燃燒理論[M]. 長沙： 國防科學技術大學出版社, 2010.

Liu Tao. Combustion Theory[M]. Changsha： University of Defense Science and Technology University Press, 2010.

[6] 李亞裕. 液體推進劑[M]. 北京： 中國宇航出版社, 2011.

Li Yayu. Liquid propellant[M]. Beijing： China Aerospace Press, 2011.

[7] 周小軍. GJB8507-2015. 航天用鈮鎢合金棒材規(guī)范[S]. 北京： 總裝備部軍標出版發(fā)行部, 2015.

Zhou Xiaojun. GJB8507-2015. Specification for niobium-tungsten alloy bars for aerospace[S]. Beijing： General Armament Department Military Standards Publication and Distribution Department, 2015.

[8] 田相軍, 羅琳. 綜合圖像灰度熵和灰度值的人臉識別方法[J]. 現(xiàn)代電子技術, 2005, 29(24)： 46-48.

Tian Xiangjun, Luo Lin. Face recognition of grayhound of entropy and grayhound of value[J]. Modern Electronic Technique, 2005, 29(24)： 46-48.

[9] 胡俊, 邱明勇, 郭紹剛. 激光切割氣體出口馬赫數(shù)對馬赫盤的影響[J].中國激光, 2008, 35(11)： 1789-1794.

Hu Jun, Qiu Mingyong, Guo Shaogang. Effect of gas outlet mach number on mach shock disk in laser cutting[J]. Chinese Journal of Lasers, 2008,35(11)： 1789-1794.