李　猛，孫　美，趙鳳起，軒春雷，儀建華，楊燕京，孫志華，王長健，許　毅

(西安近代化學研究所燃燒與爆炸技術重點實驗室，陜西 西安 710065)

引　言

固體推進劑排氣羽流是一種超音速排出噴管的氣相和凝聚相粒子混合物，并伴隨有燃燒產物的二次燃燒、湍流、電子激發(fā)和電離等，在噴管出口處下游常形成羽毛狀的發(fā)光火焰流場[1-3]。排氣羽流的氣動力結構不僅導致速度和壓力的非連續(xù)性，還影響氣流與周圍空氣的混合過程[3-4]。如機載導彈發(fā)射時排氣羽流所形成的壓力擾動使得進氣道增壓，產生的溫度與壓力畸變引起壓力機失速或停車；艦載導彈對排氣羽流的防護要求較高，由于艦面范圍有限，其上層建筑和雷達外接設備等都可能受到羽流的危害；排氣羽流對飛行中的火箭導彈的空氣動力特性和飛行穩(wěn)定性等都有較大影響。流速和溫度是排氣羽流重要的特性參數，國外學者對其開展了數值模擬及實驗驗證，并將模擬結果和實驗數據進行了對比分析[5-7]；國內近年在固體推進劑排氣羽流數值模擬方面也開展了一些研究[8-9]，建立了排氣羽流特性測試表征手段[10-12]，但對于排氣羽流場數值模擬的準確性和可靠性未進行驗證，使得模擬計算與工程實際嚴重脫節(jié)。

本研究采用基于最小自由能原理的能量計算星程序[13-15]和流體仿真商業(yè)程序ANSYS-Fluent對雙基(DB)推進劑、復合改性雙基(CMDB)推進劑及硝酸酯增塑聚醚(NEPE)推進劑進行了流場模擬計算，獲得了3種推進劑裝藥排氣羽流特性參數，采用紅外熱像儀和TDLAS流速測試系統(tǒng)對其進行了實驗驗證，將模擬結果與實驗結果進行對比，為數值模擬技術用于工程實踐及數值模型修正提供參考。

1　數值模擬

1.1　熱力學計算

根據推進劑燃燒產物組分沒有自發(fā)的變化趨勢時體系自由能最小的原理建立自由能函數方程，聯立質量和能量守恒方程，獲得給定壓強下燃燒產物組分、絕熱燃燒溫度等熱力學參數[13-15]。

在燃燒室壓強為7MPa、噴管出口壓強為0.1MPa的假設條件下，應用能量計算星程序計算獲得了雙基(DB)推進劑、復合改性雙基(CMDB)推進劑及硝酸酯增塑聚醚(NEPE)推進劑3種推進劑配方的燃燒溫度、燃燒產物組分及含量等熱力學參數，并將獲得的熱力學參數作為羽流場計算的初始條件。

1.2　計算區(qū)域及網格結構

對于推進劑裝藥排氣羽流來說，流場中各處溫度、速度等特性分布都具有軸對稱性。為簡化前處理、提高計算速度及滿足計算精度，采用二維軸對稱計算區(qū)域對羽流場進行計算，圖1為排氣羽流計算區(qū)域示意圖。

邊界條件：噴管入口處設置為壓力入口邊界，其總壓為恒定值；噴管壁面為壁面邊界條件，假設流體和壁面之間沒有熱量傳遞，固體壁面靜止不動，且為無滑移、絕熱壁面；噴管及羽流軸線為軸對稱邊界條件；其余邊界為壓力出口邊界，出口處的氣體參數按照試驗過程中的大氣條件給出。由推進劑羽流特性參數實驗可知，其采用Ф50mm標準固體火箭發(fā)動機，噴管幾何型面均為收斂-擴張的拉瓦爾型，不同的是噴管喉部半徑及出口半徑和擴張段長度。在Gambit中編寫并固化了Journal文件對計算區(qū)域進行建模和對稱結構化網格劃分，在壁面處和噴管出口后緣速度剪切層進行加密，網格數量約為兩萬。因此只用在Journal文件中改變噴喉半徑、出口半徑和擴張段長度，即可實現3種推進劑裝藥羽流特性計算的網格自動劃分。圖2為采用上述方法獲得的雙基推進劑裝藥羽流計算網格。

1.3　排氣羽流場計算

為獲得較高的計算精度，采用二階迎風格式的有限體積法對羽流流場的氣凝兩相多組分控制方程進行離散，并采用密度基求解器來模擬高馬赫數流場，湍流模型選用RNGk—ε模型[16]，然后采用有限反應速率模型和DPM模型分別研究排氣羽流多組分的復燃和凝聚相顆粒的運動規(guī)律。

由于高速運動的氣流對凝聚相Al2O3顆粒的剪切作用，排氣羽流中顆粒的直徑通常小于10～15μm，且是多種直徑組成的顆粒群，考慮到工程實際情況及保證相應計算精度的同時盡可能簡化，在數值計算過程中認為鋁顆粒直徑分布在1～10μm之間，使用Rosin-Rammler分布規(guī)律獲得羽流中凝聚相顆粒的粒度分布情況，即：

(1)

式中：F為分布函數；d為粒徑尺寸；D43為平均粒徑；n為指數系數。

對于某些包含原子及自由基的反應，需要考慮第三種物質以接受過剩的反應能量，從而使反應得以穩(wěn)定進行。為更真實地反映排氣羽流二次燃燒所帶來的溫度場變化，在分析上述3種固體推進劑的燃燒產物組分特點及含量的基礎上，對DB推進劑和CMDB推進劑，主要采用C /H /O 組分系統(tǒng)9種組分、10個基元反應的有限速率化學反應模型來模擬羽流中的二次燃燒現象；由于NEPE推進劑的燃燒產物包含Cl元素，故采用C /H /O /Cl組分系統(tǒng)11種組分，14個方程的有限速率化學反應模型來研究羽流二次燃燒。所用反應機理數據如表1所示。

表1　3種典型推進劑排氣羽流中的化學反應

為節(jié)省計算資源，首先計算單組分的氣相流場，待計算收斂且結果符合理論解后，在現有流場結果的基礎上開啟組分反應模型，設置各組分在各個邊界處的含量分布和組分之間的反應模型，從而進行組分二次燃燒的模擬。最后再開啟DPM模型，設置顆粒的材料、速度、流量等屬性，并加載UDF來控制顆粒的運動、二次燃燒及氣相、顆粒相之間的耦合來完成最終的模擬。

2　實　驗

2.1　推進劑及裝藥

選取DB、CMDB、NEPE 3種固體推進劑，其中DB推進劑和CMDB推進劑采用螺旋壓伸工藝制備，NEPE推進劑采用淤漿澆鑄工藝制備。推進劑配方見表2。

表2　3種典型推進劑的配方

根據推進劑裝藥排氣羽流檢測需要，選定Ф50mm標準固體火箭發(fā)動機，制備3種推進劑單孔管狀裝藥，其中外徑為45mm，內徑為8mm；單端包覆裝藥。同時按GJB770A方法706.1靶線法要求制備藥條試樣，并測試燃速，由3種推進劑在7MPa下的燃速、推進劑密度以及計算獲得的特征速度，根據內彈道理論設計出3種推進劑裝藥地面靜止發(fā)動機實驗所需要的噴管尺寸，見表3。

表3　實驗噴管構型尺寸

注：R為噴喉直徑；L為噴喉長度；r為入口半徑；θ1、θ2分別為收斂角和擴張角；η為擴張比

2.2　實驗裝置

羽流特性參數紅外輻射溫度采用德國IRCAM327型紅外熱像儀測試[17-18]。羽流流速采用可調諧半導體激光吸收光譜技術(TDLAS)測量[12]，即利用半導體激光器的窄線寬和可調諧特性對氣體分子的吸收譜線進行快速掃描和測量，實時反演出氣流速度。

2.3　實驗原理及方法

2.3.1羽流溫度的測量

利用紅外熱像儀對羽流流場進行實時溫度測量，根據羽流發(fā)射率計算羽流溫度場的紅外輻射溫度平均值和最高值的分布[17-18]。以火箭發(fā)動機羽流場作為紅外輻射源，在一定位置上布設紅外熱像儀和比色計對輻射溫度進行測量，利用比色溫度對被測目標發(fā)射率不依賴的特性，計算獲得了發(fā)動機羽流場發(fā)射率與時間的關系，對紅外熱像儀測得的溫度分布進行實時修正，使得測量結果的準確度大幅度提高。接通點火電源前觸發(fā)紅外熱像儀采集數據，在發(fā)動機工作完后停止采集數據。紅外熱像儀放置位置如圖3所示。

2.3.2羽流速度的測量

利用TDLAS進行流場速度的非接觸式測量，其實質是利用多普勒效應，利用頻移來實現流場速度測量[10-12]。將TDLAS流速測量裝置放置于固體火箭發(fā)動機噴口后方，TDLAS羽流流速測試系統(tǒng)測量點(激光交叉點)距發(fā)動機噴口30cm，以噴口軸心為原點，則TDLAS的激光發(fā)射點與接收點的坐標分別為(0.05，0.25)、(0.55，0.25)、(0.05，-0.25)、(0.55，-0.25)。在接通點火電源使發(fā)動機點火時，觸發(fā)TDLAS流速測量系統(tǒng)采集數據，在發(fā)動機工作完畢后停止采集數據。系統(tǒng)示意圖如圖4所示。

3　實驗與數值模擬結果對比分析

3.1　排氣羽流溫度結果的對比

采用紅外熱像儀在實驗室條件下對3種推進劑裝藥的羽流紅外輻射溫度進行測試，結果見圖5。

由圖5可以看出，CMDB推進劑羽流紅外輻射面積最大，NEPE推進劑羽流紅外輻射面積次之，且輻射溫度最高，DB推進劑羽流紅外輻射面積最小，紅外輻射溫度最低。

表4列出了3種推進劑羽流紅外輻射溫度測試結果的最大值和平均值、能量星程序計算的燃燒室溫度及模擬計算的羽流溫度。由表4可知，3種推進劑紅外輻射溫度最大值和平均值的大小順序一致，NEPE推進劑最高，CMDB推進劑次之，DB推進劑最小。通過與能量星程序計算獲得的燃燒室和完全膨脹后羽流的溫度值進行對比，發(fā)現3種推進劑羽流溫度值之間大小順序一致，能量星程序計算結果可以對羽流溫度進行定性判斷。

表4　3種典型推進劑羽流紅外輻射溫度數據

注：t1和t2分別為紅外輻射溫度的測試值和計算值。

圖6為3種推進劑裝藥排氣羽流數值模擬溫度場云圖。

由圖6的溫度云圖可以看出，由于出口壓強高于環(huán)境壓強，燃氣出噴管即開始膨脹解壓，溫度降低，進而截斷斜激波造成溫度間斷升高，激波遇到混合層邊界反射成膨脹波，從而導致新一輪的膨脹開始，如此往復，馬赫盤和溫度間斷不斷出現，直到射流與外界環(huán)境壓力相同。同時，二次燃燒主要發(fā)生在排氣羽流混合層內，并且隨著排氣羽流距離的增加而逐漸與主流混合，由于混合層二次燃燒的大量放熱，混合層的溫度比核心區(qū)溫度要高幾百攝氏度。由于DB推進劑的氧平衡較低，使得DB推進劑二次燃燒最劇烈，CMDB推進劑次之，NEPE推進劑二次燃燒相對最弱。圖7為3種推進劑裝藥發(fā)動機軸線上的溫度分布曲線。

由圖7可看出，NEPE推進劑羽流溫度最高，CMDB推進劑次之，DB推進劑最小，流場計算獲得的羽流溫度與紅外熱像儀測試結果一致。由于排氣羽流的膨脹壓縮波系，使得排氣羽流溫度出現振蕩，但經過馬赫盤后，燃氣的溫度突然升高。由于NEPE推進劑的能量較高，燃燒溫度高，噴管噴出燃氣的溫度也相應較高，同時由于溫度的升高，在混合層處化學反應速率加快，復燃產生的熱量使混合層局部溫度更高，從而使此區(qū)域等溫線分布更密。CMDB推進劑次之，DB推進劑能量最小，推進劑燃燒溫度也最小。

綜上所述，羽流溫度作為排氣羽流特性重要參數之一，本實驗測試采用的紅外熱像儀是非接觸式測量，由于輻射率等因素的影響，測試值比真實值偏小；所采用的羽流場數值模型是基于能量守恒定律的，相對于外界是絕熱狀態(tài)，與實際排氣羽流溫度有差別，計算值比真實值要偏高。

3.2　排氣羽流流速結果對比

采用TDLAS對3種推進劑裝藥的羽流流速進行測試，結果見圖8。

由圖8可知，羽流流速測試曲線與推進劑裝藥工作過程p—t曲線一致，經分析和數據處理，獲得3種推進劑裝藥羽流流速平均值分別為879.0、1103.0、1457.1m/s，總體來看，由于NEPE推進劑能量高、CMDB推進劑次之，DB推進劑能量最小，由一維等熵流動假設可知，羽流流速與推進劑能量性能大小順序一致，本研究測試結果可靠。在羽流場計算數據中，選取TDLAS的激光發(fā)射點、接收點及激光交叉點坐標處的流速數據，經數學處理獲得羽流流速模擬值為1002.38、1279.01、1243.16m/s，模擬值與計算值之間偏差為14.04%、15.9%、-14.7%。

圖9為3種推進劑裝藥排氣羽流速度場云圖。

由圖9可知，氣流在噴喉處被加速到音速，隨后流通截面增大，氣流繼續(xù)加速到超音速，由于噴管出口處的壓強小于外界環(huán)境壓強，處于高度過膨脹狀態(tài)，隨著氣流經過噴管出口，伴隨著壓強升高、激波產生，激波形狀是從噴口邊緣發(fā)出的收縮截錐形，截錐頂部為一盤形正激波即馬赫盤，經過馬赫盤后，速度迅速降低。此后，羽流再次經過一個膨脹波系，速度升高，直到下一個馬赫盤出現，如此反復，羽流核心區(qū)和邊界層最終混合均勻。

圖10為3種推進劑在發(fā)動機軸線上的羽流流速。

由圖10可看出，NEPE推進劑羽流流速最高，CMDB次之，DB推進劑最小。流場計算獲得的羽流流速與TDLAS羽流流速測試結果一致。單就一種推進劑羽流流速軸線分布情況來看，羽流流速數據經過多次振蕩后達到平衡狀態(tài)，隨著遠離噴口速度逐漸衰減到零，其中NEPE推進劑羽流流速振蕩次數最多，其次為CMDB推進劑，DB推進劑振蕩次數最少，這與排氣羽流馬赫盤的形成分析一致。

3種推進劑羽流流速模擬結果與實驗測試結果大小順序一致，分析其原因主要是：(1)TDLAS實驗技術選取的是1392nm附近的H2O吸收譜線，并且在均勻氣相風洞中進行標定，而本研究的推進劑排氣羽流是包含氣、凝兩相多組分、含化學反應的流場，其速度呈現明顯的滯后，根據文獻[11-12]分析，測量值與理論值偏差可能是流場的不均勻導致；(2)TDLAS實驗測試羽流流速是在測量光學路徑上的平均氣流流速，由于考慮到實驗裝置的安全性，TDLAS裝置架設位置不可能完全貼合排氣羽流場，造成排氣羽流場外緣很多靜止的空氣也被考慮到羽流流速中，導致測試結果偏低。

4　結　論

(1)通過日志文件驅動Gambit和ANSYS-Fluent程序讀入Journal文件來自動生成3種推進劑羽流對稱結構化網格和開展氣凝兩相化學反應流場模擬，采用C/H/O 9種組分、10個基元反應模擬DB、CMDB推進劑羽流二次燃燒，采用C/H/O/Cl 11種組分、14個方程來模擬NEPE推進劑羽流二次燃燒，只需輸入燃氣成分及含量即可獲得羽流溫度和速度分布。

(2)采用紅外熱像儀和可調諧半導體激光吸收光譜技術(TDLAS)對排氣羽流特性參數羽流溫度和羽流速度進行了實驗驗證，將模擬結果與測試結果進行對比，對產生偏差的原因進行了分析，分析結果可對模型修正提供一些途徑，為實驗系統(tǒng)的設置及實驗方法的改進提供參考。

參考文獻:

[1]Thorn L B, Wharton W W. Terminology and assessment methods of solid propellant rocket exhaust signature, AGARD-AR-287[R]. Loughton: AGARD, 1993.

[2]Chastenet J C, Kessel P A , Lawrence R, et al. Rocket motor plume technology, AGARD-LS-188[R]. Loughton: AGARD, 1993.

[3]李猛,王宏.火箭發(fā)動機排氣羽流特征信號預示技術研究綜述[J].飛航導彈,2010(3):54-57.

LI Meng, WANG Hong. A review studies on the technology for prediction of rocket motor plume signature[J]. Aerodynamic Missile Journal,2010(3):54-57.

[4]Guery J F, Chang I S, Shimada T, et al. Solid propulsion for space applications: an updated road map[J]. Acta Astronautics, 2010, 66: 201-219.

[5]Troyes J, Dubois I, Borie V, et al. Multi-phase reactive numerical simulations of a model solid rocket motor exhaust jet[C]∥42th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit. New York: AIAA, 2006:1-13.

[6]Wang X Y ，Yuko J R.Thermal analysis on plume heating of the main engine on the crew exploration vehicle service module, NASA/TM-215049[R].Washington D C:NASA, 2007.

[7]Gueyffier D, Fromentin-Denoziere B, Simon J, et al. Numerical simulation of ionized rocket plumes[C]∥44th AIAA Plasmadynamics and Laser Conference. New York: AIAA, 2013:1-13.

[8]李猛, 張曉宏, 孫美, 等. 改性雙基推進劑兩相化學反應羽流特性研究[J]. 彈箭與制導學報, 2012, 32(1): 123-126.

LI Meng, ZHANG Xiao-hong, SUN Mei, et al. The study on plume characteristic of two phase chemical reaction of modified double-base propellant[J]. Journal of Projectiles, Rockets, Missiles Guidance, 2012, 32(1): 123-126.

[9]李崢, 向紅軍, 張小英. 復合推進劑固體火箭發(fā)動機噴流流場數值模擬[J]. 固體火箭技術, 2014, 37(1): 37-42.

LI Zheng, XIANG Hong-jun, ZHANG Xiao-ying. Numerical simulation of composite solid propellant rocket motor exhaust plume[J]. Journal of Solid Rocket Technology, 2014, 37(1): 37-42.

[10] 楊燕京, 趙鳳起, 孫美, 等. 基于TDLAS技術的雙基系推進劑裝藥排氣羽流流速特性研究[J]. 火炸藥學報, 2016，39(3)：70-74.

YANG Yan-jing, ZHAO Feng-qi, SUN Mei, et al. Investigations on the exhaust plume velocity characteristics of double base and CMDB propellants charge based on TDLAS technique[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao), 2016，39(3)：70-74.

[11] 楊斌, 黃斌, 劉佩進, 等. 利用TDLAS技術評估火箭基組合循環(huán)發(fā)動機試驗性能[J]. 宇航學報, 2015，36(7)：840-848.

YANG Bin, HUANG Bin, LIU Pei-jin, et al. Test performance evaluation of rocket based combined cycle engine using TDLAS technique[J]. Journal of Astronautics, 2015，36(7)：840-848.

[12] 賈良權, 劉文清, 闞瑞峰, 等. 波長調制-TDLAS 技術測量風洞中氧氣流速方法研究[J]. 中國激光, 2015 (7): 304-311.

JIA Liang-quan, LIU Wen-qing, KAN Rui-feng, et al. Study on oxygen velocity measurement in wind tunnel by wavelength modulation-TDLAS technology[J]. Chinese Journal of Lasers, 2015 (7): 304-311.

[13] McBride B J, Gordon S. Computer program for calculation of chemical equilibrium compositions and applications: II users manual and program description, NASA RP-1311[R]. Washington D C:NASA, 1996.

[14] 李猛, 趙鳳起, 徐司雨, 等. 三種能量計算程序在推進劑配方設計中的比較[J]. 火炸藥學報, 2013, 36(3): 73-77.

LI Meng, ZHAO Feng-qi, XU Si-yu, et al. Comparison of three kinds of energy calculation programs in formulation design of solid propellants[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao), 2013, 36(3): 73-77.

[15] 田德余,趙鳳起,劉劍洪.含能材料及相關物手冊[M].北京:國防工業(yè)出版社,2011.

[16] Yakhot V , Orszag S A. Renormalization group analysis of turbulence: I. basic theory[J]. Journal of Scientific Computing, 1986,1(1):1-51.

[17] 儀建華. 火箭發(fā)動機羽流紅外輻射溫度動態(tài)測量方法：CN，201610846494.4[P]. 2016-09-23.

[18] 張波, 呂事桂. 基于雙響應波段工作的紅外熱像儀測溫原理與誤差分析[J]. 宇航計測技術, 2010, 30(6): 11-15.

ZHANG Bo, Lü Shi-gui. Temperature measurement of the thermal infrared imager based on two response working wavebands and its error analysis[J]. Journal of Astronautic Metrology and Measurement, 2010, 30(6): 11-15.