同心筒發(fā)射燃氣流二次燃燒數(shù)值研究及導流板結(jié)構改進

胡曉磊，樂貴高，馬大為，朱忠領

(南京理工大學 機械工程學院，江蘇 南京210094)

0 引言

同心筒熱發(fā)射具有獨立排導、結(jié)構簡單、可靠性高、雷達非同步和全方位發(fā)射等特點，在導彈發(fā)射中引起普遍關注［1］。由于固體火箭發(fā)動機噴出的燃氣射流含未完全燃燒的氣體，主要有H2O、CO2、H2、CO、N2、HCl、O2、Cl2、H、OH 等多種組分，在流場分析時需要考慮燃氣射流組分二次燃燒的影響［2］。如果忽略二次燃燒，計算所得的流場有較大失真。

針對同心筒熱發(fā)射問題，研究學者普遍采用計算液體力學數(shù)值模擬的方法和實驗方法對發(fā)射筒內(nèi)燃氣射流流動的規(guī)律進行分析。文獻［3］采用多重網(wǎng)格法研究了內(nèi)外筒間隙、導流錐高度對燃氣流流場的影響。熊永亮等［4］采用動網(wǎng)格技術研究了內(nèi)外筒尺寸對旁泄流的影響。藺翠郎等［5］對同心筒發(fā)射過程中筒內(nèi)的純氣相流場熱效應進行了分析。姜毅等［6-8］對同心筒進行了大量的實驗研究和數(shù)值模擬工作，研究了尾部尺寸、筒口安裝導流裝置對燃氣排導的影響，并且提出了一種新型“引射同心筒”的概念。何朝勛等［9］用數(shù)值模擬的方法研究了導流柵結(jié)構和型面對同心筒筒內(nèi)流場的影響。于勇等［10］提出同心筒氣體動力學過程的簡化模型，給出了估算最小狹縫寬度和平均附加彈射力的方法。以上研究均沒有考慮燃氣射流二次燃燒對流場的影響。

針對某同心筒熱發(fā)射裝置，采用域動分層網(wǎng)格更新方法數(shù)值研究同心筒熱發(fā)射導彈的運動，建立考慮二次燃燒的11 組分12 步化學反應動力學模型，得到典型時刻流場溫度的分布。分析了導流板結(jié)構對導彈出筒過程的影響。研究的數(shù)據(jù)為同心筒熱發(fā)射裝置的熱防護和結(jié)構優(yōu)化提供理論參考依據(jù)。

1 數(shù)學模型

1.1 控制方程

針對同心筒熱發(fā)射過程中二次燃燒特點，采用有限速率化學反應模型建立考慮二次燃燒的二維軸對稱組分擴散控制方程

式中:

(1)式和(2)式中:Ω(t)表示體積可變控制體，Γ(t)表示其外表面，i =1，…，N;N 表示混合氣體的組分總數(shù);ρi為第i 組分的密度;ρ 為混合物的密度;p、T分別為混合物的壓力和溫度;u、v 分別為流體在x、y方向的速度;uw和vw分別為運動網(wǎng)格在x、y 方向的速度;hi為組分i 單位質(zhì)量的焓值;τ 為粘性應力張量;k 為熱擴散系數(shù);Di是組分i 的擴散系數(shù);ω·i為組分i 的質(zhì)量生成率;E 為混合物單位的總能。

利用有限體積法，采用顯示格式求解控制方程組(1)式，其中對流項采用對激波具有較高捕捉精度的3 階MUSCL 格式離散［10］，時間項采用4 階Runge-kutta 法處理，湍流模型選用重整化群RNG k-ε二方程模型［11］，壁面附近采用標準壁面函數(shù)。

1.2 數(shù)值方法

選擇Arrhenius 定律描述同心筒內(nèi)二次燃燒的詳細化學反應機理，忽略湍流脈動對化學反應過程的影響。

考慮以如下形式的第r 個反應:

式中:N 為系統(tǒng)中化學物質(zhì)數(shù)目;v'i，r為反應r 中反應物i 的化學計量系統(tǒng);v″i，r為反應r 中生成物i 化學計量系統(tǒng);Mi為第i 物質(zhì)的符號;kf，r為反應r 的正向速率常數(shù);kb，r為反應r 逆向速率常數(shù)。

反應r 中物質(zhì)i 產(chǎn)生/分解摩爾速度以以下公式給出:

式中:Nr為反應r 化學物質(zhì)數(shù)目;Cj，r為反應r 每種反應物或生成物j 摩爾濃度;η'j.r為反應r 每種反應物或生成物j 正向速度指數(shù);η″j.r為反應r 每種反應物或生成物j 的逆向速度指數(shù);忽略第三體的影響。

反應r 正向速率常數(shù)kf，r通過Arrhenius 公式計算

式中:Ar為指前因子(恒定單位);βr為溫度指數(shù)(無量綱);Er為反應活化能(J/kmol);R 為氣體常數(shù)(J/(kmol·K))。

采用11 組分12 步化學反應動力學模型［13］，化學式及相關參數(shù)見表1. 其中M 為第三體(不參與化學反應)，本例中近似認為是HCl.

表1 化學反應模型Tab.1 Chemical reaction model

1.3 域動分層網(wǎng)格更新方法

根據(jù)導彈在同心筒內(nèi)運動過程，計算模型采用結(jié)構化網(wǎng)格對流場區(qū)域進行空間離散，計算域的變化采用彈簧近似法和動態(tài)分層法相結(jié)合的域動分層網(wǎng)格更新技術進行處理，具體細節(jié)可參考文獻［6］.

2 計算模型和邊界條件

2.1 計算區(qū)域

同心筒熱發(fā)射裝置由外筒、內(nèi)筒、后蓋構成，其中后蓋采用錐形導流型面。圖1(a)為某同心筒熱發(fā)射裝置示意圖，圖1(b)為含有導流板的同心筒發(fā)射裝置示意圖，其中L 為導流板的長度。

針對同心筒發(fā)射裝置的結(jié)構特點，考慮到計算條件的限制，在不影響真實流場的前提下，采用軸對稱計算模型進行流場分析。圖2為噴管附近的網(wǎng)格劃分情況。為了保證數(shù)值研究的精度，分析時流場全部采用結(jié)構化網(wǎng)格。

圖1 同心筒熱發(fā)射裝置示意圖Fig.1 Sketch of concentric canister launcher

圖2 噴管附近網(wǎng)格劃分Fig.2 Grid around the nozzle exit

2.2 邊界條件

噴管入口采用壓力入口邊界條件，總壓為9.119 ×106Pa，總溫為3 387 K. 燃氣入口主要組分和質(zhì)量分數(shù)見表2. 同心筒熱發(fā)射周圍流場為靜止大氣條件;靜壓為1.013 ×105Pa，靜溫為300 K，大氣中N2的質(zhì)量分數(shù)為0.77，O2的質(zhì)量分數(shù)為0.23.

表2 燃氣入口組分和質(zhì)量分數(shù)Tab.2 Species and mass fraction of inlet

3 計算結(jié)果及討論

3.1 計算方法的驗證

為了檢驗所建數(shù)值方法的正確性，采用11 組分12 步化學反應動力學模型對文獻［2］中的固體火箭發(fā)動機尾焰流場進行研究。計算結(jié)果如圖3所示。

圖3為數(shù)值計算得到的溫度場云圖，從圖中可見，在噴管出口處溫度較高。并且數(shù)值計算得到的溫度場分布與圖4中尾焰實驗熱像圖［2］分布趨勢基本一致。表3是軸線上各觀測點溫度計算值與實驗值對比。通過計算值和實驗值［2］的對比發(fā)現(xiàn)最大溫度誤差為5.31%. 表明所建立的數(shù)值方法是可行的，能夠運用到燃氣射流二次燃燒流場分析中。

圖3 數(shù)值計算溫度場云圖Fig.3 Temperature contour of numerical simulation

圖4 尾焰熱像圖［2］Fig.4 Thermal image of exhaust plume

表3 軸線上觀測點溫度計算值與實驗值對比Tab.3 Comparison of calculated and experimental results of temperature along axis

3.2 網(wǎng)格無關性檢測

由于導彈在筒內(nèi)運動時間短暫，燃氣流二次燃燒流場的計算精度對網(wǎng)格密度有較強的依賴性。因此，網(wǎng)格無關性的檢驗至關重要。為了選取合適的計算網(wǎng)格，對3 種網(wǎng)格數(shù)對流場進行網(wǎng)格無關性測試，開始分別取98 380(網(wǎng)格A)，142 270(網(wǎng)格B)，205 690(網(wǎng)格C)網(wǎng)格個數(shù)，其中燃氣溫度為3 190 K，對5 ms 時刻后蓋的溫度峰值及其所在的位置進行對比。計算模型網(wǎng)格無關性測試如表4所示。

從表4可以發(fā)現(xiàn)，網(wǎng)格B 與網(wǎng)格C 在5 ms 時刻峰值溫度差異為3.33%，Y 坐標值相差2.67%，雖然更高分辨率的網(wǎng)格能夠更好地捕捉流動細節(jié)，但考慮計算條件的限制，網(wǎng)格B 已滿足的研究，故采用網(wǎng)格B 分析同心筒燃氣流二次燃燒流場。

表4 5 ms 時3 種網(wǎng)格溫度峰值對比Tab.4 Comparison of maximum temperatures with 3 different meshes at 5 ms

3.3 同心筒化學反應動力學流場分析

圖5 不同時刻流場的速度云圖Fig.5 Velocity contours at different time

圖5(a)～圖5(d)為0.10 s、0.20 s、0.32 s 和0.35 s 時刻流場的速度云圖。從整個導彈運動過程來看，在導彈出筒之前，從筒內(nèi)噴出的燃氣射流沿著導彈兩側(cè)壁面向前流動，并吞噬初始流場、追趕包圍導彈;在導彈側(cè)面和頭部，形成一個向前突出的氣團。導彈出筒之后，燃氣射流開始離開發(fā)射筒，外界空氣與燃氣射流開始發(fā)生化學反應，由于受到筒內(nèi)排除的燃氣流、導彈尾部的燃氣射流和彈底結(jié)構的影響，高度欠膨脹射流向筒口兩邊排去，形成反濺激波。隨著導彈尾部離發(fā)射筒口越來越遠，彈底部空間逐漸增大，反濺激波逐漸減弱。

圖6(a)～圖6(d)為0.10 s、0.20 s、0.32 s 和0.35 s 時刻溫度場云圖。從圖6中溫度場云圖可見，在導彈發(fā)射過程中，筒內(nèi)一直處于高溫燃氣射流作用中。除燃氣射流核心區(qū)域外，筒內(nèi)大部分區(qū)域溫度在3 000 K左右。出筒之前，高溫燃氣向?qū)椷\動方向流動，包圍筒外的彈體。由此可見高溫是同心筒發(fā)射的特點。

圖6 不同時刻溫度場云圖Fig.6 Temperature contours at different time

圖7是不同時刻彈體表面溫度曲線，其中0.10 s時刻和0.20 s 時刻為導彈出筒之前的彈體表面溫度，0.32 s 和0.35 s 為導彈出筒之后溫度。從曲線的變化趨勢可見，導彈出筒之前，彈體表面溫度明顯高于導彈出筒之后。通過比較0.10 s 和0.20 s 時刻彈體表面溫度可見，0.20 s 時刻在0 m ＜x ＜1.65 m區(qū)間內(nèi)彈體的表面溫度高于0.10 s 時刻。這是由于隨著導彈的運動，燃氣在筒內(nèi)積聚，造成彈體表面溫度升高。在0.20 s 時刻1.65 m ＜x ＜1.78 m 區(qū)間內(nèi)，由于導彈頭部開始沖破燃氣高溫區(qū)域，使得此區(qū)域的彈體表面溫度低于0.10 s 時刻。而且隨著導彈離筒內(nèi)聚集的高溫燃氣越來越遠，彈體表面的溫度逐漸下降。

圖7 不同時刻彈體表面溫度曲線Fig.7 Missile body temperatures at different time

圖8為0.20 s 時刻和0.32 s 時刻同心筒發(fā)射流場密度云圖。從云圖可見，導彈在發(fā)射過程中燃氣射流與外界空氣之間有明顯的密度間斷，這是由于燃氣射流區(qū)域的氣體受到激波加熱后溫度升高，密度相應降低。因此在燃氣射流推動空氣向周圍運動的同時，二者之間形成一個形狀不規(guī)則的接觸間斷。

圖8 不同時刻密度云圖Fig.8 Density contours at different time

圖9(a)和圖9(b)是0.32 s 時刻O 和CO2的質(zhì)量分布。針對本文計算情況，O 是原來沒有的組分，因此可以通過分析O 的質(zhì)量分數(shù)分布反映氧化反應的劇烈程度。從圖9(a)可以看出，在0.32 s 時刻筒內(nèi)O 的質(zhì)量分數(shù)幾乎為0，說明筒內(nèi)被高溫燃氣充滿，空氣較少，即缺少發(fā)生氧化反應的O2. 導彈出筒之后，筒內(nèi)的燃氣射流迅速與空氣中的O2發(fā)生化學反應。在導彈兩側(cè)含有較多的O，質(zhì)量分數(shù)達0.004. 可見在導彈出筒之后，氧化反應主要發(fā)生在導彈頭部和兩彈身兩側(cè)壁面附近的混合層。從外邊界層到混合核心區(qū)域質(zhì)量分數(shù)由0.000 5 變化到0.005 2. 尤為明顯的是燃氣射流的核心區(qū)域，O 的質(zhì)量分數(shù)為0，即在此區(qū)域可認為不發(fā)生化學反應。這與文獻［13］中的結(jié)論一致。

圖9 0.32 s 時刻反應產(chǎn)物質(zhì)量分布Fig.9 Mass fraction contours of reaction products at 0.32 s

由于CO2只在燃氣射流中存在，可用其質(zhì)量分數(shù)分布表示該時刻燃氣射流噴出后所覆蓋的流場區(qū)域，如圖9(b)所示。從圖可以看出燃氣射流從筒內(nèi)噴出后所覆蓋的流場區(qū)域與圖6(c)流場中的高溫區(qū)域完全重合，且燃氣射流與空氣之間的界面非常清晰，這主要是由于高溫高速的燃氣射流在推動前方空氣的過程中與外界空氣的熱量和質(zhì)量交換速度相對較慢導致的。

3.4 化學反應動力學對溫度場的影響

圖10(a)給出了多組分凍結(jié)流和多組分反應流溫度場對比。圖10(b)給出了在0.20 s 時刻單一組分燃氣、多組分凍結(jié)流和多組分反應流軸線溫度對比。從圖10(a)凍結(jié)流和反應流的溫度場云圖對比可見，反應流的溫度場波及范圍要大于凍結(jié)流流場。由圖10(b)中可見，在發(fā)射筒內(nèi)考慮/不考慮組分和化學反應對軸線溫度無明顯影響。這說明在同心筒內(nèi)有限的空間內(nèi)，幾乎不發(fā)生化學反應。而在筒外，單一組分燃氣流、多組分凍結(jié)流與含有化學反應反應流軸線溫度相差比較大。這是由于單一組分燃氣流中沒有考慮組分的擴散，而多組分凍結(jié)流流沒有考慮從筒內(nèi)排除的燃氣流與周圍空氣混合發(fā)生強烈的化學反應。從圖10(b)可見筒外多組分反應流流場的溫度比多組分凍結(jié)流流場的溫度高500 K 左右，比單組分燃氣流場溫度最高高1 000 K. 這是由于筒外的富燃燃氣與空氣中的O2在筒口附近混合，混合流場中不僅有膨脹波而且存在壓縮波，壓縮波后燃氣溫度、壓強上升，更易引起化學反應并釋放熱量，導致反應區(qū)域的溫度比無反應和單組分燃氣的溫度高。這與文獻［2］中實驗結(jié)果一致。對比有無化學反應時流場溫度表明，同心筒發(fā)射裝置流場考慮二次燃燒是必要的。如果忽略二次燃燒，將導致同心筒溫度場計算誤差增大。

圖10 0.20 s 時刻不同工況溫度對比Fig.10 Comparison of temperatures at 0.20 s

3.5 筒口導流板對流場的影響

為了研究筒口導流板對導彈出筒時彈體表面溫度的影響，選取3 種形式的導流板，如圖1(b)所示，導流板長度分別是0.5s，1.5s 和3.0s，其中s 為內(nèi)外筒間隙。將這3 種工況分別與無導流板的結(jié)構進行對比。

圖11(a)是含有不同長度的導流板和不含導流板在0.1 s 時刻彈體表面溫度的對比。從圖中可見，當導彈出筒時，導流板對內(nèi)外筒之間的燃氣具有很好的導向作用。從溫度曲線上可以看出，加上導流板能夠使彈體表面溫度下降1 000 K 以上。導流板長度從0.5s 增加到3.0s 時，彈體表面的溫度先降低后升高，在1.5s 時，彈體表面的溫度最低。圖11(b)～圖11(e)為不同導流板長度時筒口流場的矢量圖。從圖中可看出，無論有無導流板，內(nèi)外筒之間都存在倒吸現(xiàn)象。這是由于導彈運動過程中，內(nèi)筒與彈壁之間大部分區(qū)域為負壓區(qū)，從內(nèi)外筒之間排出的燃氣在發(fā)射筒出口處產(chǎn)生抽吸效應進入彈筒間隙。

圖11 不同導流板長度時筒口流場矢量圖Fig.11 Velocity vectors for different guider lengths

當導流板存在時，從內(nèi)外筒之間排除的部分燃氣在筒口與周圍的空氣以一定比例混合形成漩渦后進入導彈與內(nèi)筒之間。導流板長度為0.5s 和3.0s時，漩渦在筒口附近形成。當導流板長度為1.5s時，由于燃氣排導順暢，漩渦離筒口較遠，進入彈筒間隙的高溫燃氣流量要比0.5s 和3.0s 時少，所以導流板長度為1.5s 時彈底壁面溫度最低。當導流板長度為3.0s 時，燃氣出口截面增大，從內(nèi)外筒間隙排出的高溫燃氣徑向速度偏小，沿軸線方向高溫燃氣繞過導流板進入彈筒間隙的流量比0.5s 時要多，所以導流板高度為3.0s 時的彈體壁面溫度高于0.5s 時。筒口無導流板時從內(nèi)外筒間隙排出的高溫燃氣直接與空氣混合進入彈筒間隙，進入筒內(nèi)的高溫燃氣流量最多，因而溫度最高。

根據(jù)以上分析可知對于同心筒熱發(fā)射裝置來說，導流板的長度最優(yōu)值在1.5s 附近。

4 結(jié)論

1)建立了耦合有限速率化學反應和導彈運動的同心筒熱發(fā)射計算模型，采用11 組分12 步基元化學反應動力學模型和域動分層網(wǎng)格更新方法，實現(xiàn)了同心筒發(fā)射過程H2/CO 混合燃氣二次燃燒流場的數(shù)值模擬。

2)導彈出筒之后，由于筒內(nèi)燃氣流、導彈尾部燃氣射流和彈底結(jié)構的綜合影響，高度欠膨脹射流在導彈尾部形成反濺激波。隨著導彈尾部離發(fā)射筒口越來越遠，反濺激波逐漸減弱。在此過程中彈體表面溫度先升高后降低。

3)根據(jù)不同時刻流場中典型組分的分布，說明在同心筒內(nèi)幾乎不發(fā)生化學反應。二次燃燒主要發(fā)生在筒口燃氣射流與空氣混合區(qū)域。

4)同心筒發(fā)射裝置多組分反應流軸線的溫度比凍結(jié)流高500 K 左右，比單一組分燃氣流高1 000 K 左右。可見考慮二次燃燒對同心筒熱發(fā)射流場作用是必要的。

5)同心筒筒口增加導流板能夠降低導彈出筒時彈體表面的溫度。同心筒導流板長度的最優(yōu)值在1.5s 附近。

References)

［1］Yagla J J. Concentric canister launcher［J］. Naval Engineer Journal，1997，1(9):313 -327.

［2］張光喜，周為民，張鋼錘，等.固體火箭發(fā)動機尾焰流場特性研究［J］.固體火箭技術，2008，31(1):19 -23.ZHANG Guang-xi，ZHOU Wei-min，ZHANG Gang-chui，et al.Study on characteristics of flow field in exhaust plu-me for SRM［J］. Journal of Solid Rocket Technology，2008，31(1):19 -23.(in Chinese)

［3］苗佩云，袁曾鳳. 同心筒式發(fā)射時筒內(nèi)流場機理及內(nèi)外筒間隙的影響［J］. 戰(zhàn)術導彈技術，2006，29(1):8 -13.MIAO Pei-yun，YUAN Zeng-feng. The effect of flow mechanism and annular sizes on concentric canister launcher［J］. Tactical Missile Technology，2006，29(1):8 -13. (in Chinese)

［4］熊永亮，郜冶，李燕良. 同心筒發(fā)射中旁泄流影響的數(shù)值研究［J］.彈箭與制導學報，2007，27(4):194 -197.XIONG Yong-liang，GAO Ye，LI Yan-liang. Numerical investigate the influence of the side releasing flow on the concentric canister launch［J］. Journal of Projectiles，Rocket，Missiles and Guidance，2007，27(4):194 -197. (in Chinese)

［5］藺翠郎，畢世華. 同心發(fā)射筒內(nèi)燃氣流溫度場的數(shù)值模擬［J］.彈箭與制導學報，2007，27(5):160 -162.LIN Cui-lang，BI Shi-hua. Numerical simulation of the combustion-gas temperature field in concentric canister launcher［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2007，27(5):160 -162. (in Chinese)

［6］姜毅，郝繼光，劉群. 同心筒垂直發(fā)射裝置排導燃氣的改進［J］. 北京理工大學學報，2007，27(2):95 -99.JIANG Yi，HAO Ji-guang，LIU Qun. Improvement measures of exhausting the jet in the concentric vertical launching equipment［J］. Transaction of Beijing Institute of Technology，2007，27(2):95 -99. (in Chinese)

［7］姜毅，郝繼光，傅德彬. 導彈發(fā)射過程三維非定常數(shù)值模擬［J］.兵工學報，2008，29(8):911 -915.JIANG Yi，HAO Ji-guang，F(xiàn)U De-bin. 3D unsteady numerical simulation of missile launching［J］. Acta Armamentarii，2008，29(8):911 -915. (in Chinese)

［8］姜毅，郝繼光. 新型“引射同心筒”垂直發(fā)射裝置理論及實驗研究［J］. 宇航學報，2008，29(1):236 -241.JIANG Yi，HAO Ji-guang. Study on theory and test of a new type“concentric canisterwith jet flow”vertical launcher［J］. Journal of Astronautics，2008，29(1):236 -241. (in Chinese)

［9］何朝勛，戴宗妙. 同心筒式發(fā)射裝置燃氣排導研究［J］. 艦船科學技術，2007，29(1):71 -75.HE Chao-xun，DAI Zong-miao. Research into the jet exhaust of concentric canister launcher［J］. Ship Science and Technology，2007，29(1):71 -75. (in Chinese)

［10］于勇，徐新文，傅德彬. 同心筒發(fā)射裝置燃氣排導的氣體動力學原理分析［J］. 航空動力學報，2012，27(3):501 -507.YU Yong，XU Xin-wen，F(xiàn)U De-bin. Aerodynamic theory analysis of gas exhausting from concentric canister launcher［J］. Journal of Aerospace Power，2012，27(3):501 -507. (in Chinese)

［11］Shu C W，Osher S. Efficient implementation of essentially nonoscillatory shock capturing schemesⅡ［J］. Journal of Computational Physics，1989，83:32 -78.

［12］Shih T H，Liou W W，Shabbir A，et al. A new k-ε eddy-viscosity model for high reynolds number turbulent flows-model development and validation［J］. Computers ＆ Fluids，1995，24(3):227-238.

［13］Jensen D E，Jones G A. Reaction rate coefficients for Flame calculations［J］. Combustion and Flame，1978，32:1 -34.

［14］Troyes J，Boischot A. Multi-phase reaction numerical simulations of a model solid rocket motor exhaust jet［R］. California:American Institute of Aeronautics and Astronautics，2006:1 -13.