孫 行， 聶萬(wàn)勝*， 蔡紅華， 陳 朋， 黃衛(wèi)東， 石天一

(1.航天工程大學(xué)宇航科學(xué)與技術(shù)系，北京 101416；2.航天員科研訓(xùn)練中心，北京 100094)

推進(jìn)劑在火箭發(fā)動(dòng)機(jī)中燃燒，經(jīng)拉瓦爾噴管?chē)姵?，形成高溫尾焰，尾焰組分發(fā)生電離，形成稠密不均勻的弱電離等離子體[1-2]，對(duì)火箭進(jìn)行測(cè)控時(shí)發(fā)現(xiàn)，電磁波穿過(guò)火箭尾焰時(shí)會(huì)受到嚴(yán)重干擾[3-8]。研究火箭尾焰等離子體的各項(xiàng)特征參數(shù)(等離子體濃度、等離子體頻率)對(duì)于研究火箭尾焰對(duì)電磁波的干擾作用意義重大。

近年來(lái)，中外對(duì)火箭尾焰等離子體特征參數(shù)的研究不斷深入。Fromentin-Denoziere等[7]描述了尾焰中堿金屬元素的電離過(guò)程；Smoot等[8]描述了尾焰中自由電子的分布規(guī)律，預(yù)測(cè)了雷達(dá)信號(hào)在尾焰中的衰減；Troyes等[9]采用數(shù)值計(jì)算的方法對(duì)固體火箭推進(jìn)劑是否含鋁兩種情況的尾焰流場(chǎng)進(jìn)行了研究；Chen[10]對(duì)等離子體特征參數(shù)的計(jì)算方法進(jìn)行了描述。清華大學(xué)的楊任剛等[11]在采用RNGk-ε湍流模型描述氣相的同時(shí)采用離散相模型描述固相，研究了尾焰氣固兩相射流流場(chǎng)；石雁祥等[5-6]在考慮火箭尾焰為塵埃等離子體的基礎(chǔ)上，對(duì)火箭尾焰等離子體特征參數(shù)進(jìn)行了大量的理論研究；郭祥天等[12]建立了固體火箭發(fā)動(dòng)機(jī)尾焰的氣固兩相流計(jì)算模型，研究了固相顆粒Al2O3對(duì)發(fā)動(dòng)機(jī)尾焰流場(chǎng)的影響。研究發(fā)現(xiàn)，相較于液體火箭，固體火箭尾焰對(duì)電磁波的干擾作用更加明顯。這些研究多以火箭尾焰流場(chǎng)參數(shù)的經(jīng)驗(yàn)值為基礎(chǔ)，計(jì)算等離子體特征參數(shù)，本文以某型號(hào)固體火箭為研究對(duì)象，對(duì)其不同工況下的尾焰流場(chǎng)進(jìn)行仿真，獲取準(zhǔn)確的尾焰流場(chǎng)參數(shù)；在此基礎(chǔ)上，計(jì)算尾焰等離子體特征參數(shù)分布，分析其分布規(guī)律，并討論不同工況對(duì)尾焰等離子體特征參數(shù)的影響。

1 物理模型和計(jì)算方法

1.1 控制方程

基于Navier-Stokes方程計(jì)算火箭尾焰流場(chǎng)，添加組分輸運(yùn)方程，得到如下四個(gè)控制方程[13]。

1.1.1 質(zhì)量守恒方程

(1)

式(1)中：ρ為流體的密度；v為流體的速度矢量。

1.1.2 動(dòng)量守恒方程

(2)

1.1.3 能量守恒方程

(3)

式(3)中：cp為定壓比熱容；T為流體的溫度；k為流體的傳熱系數(shù)；ST為流體的內(nèi)熱源及由于黏性作用而使流體機(jī)械能轉(zhuǎn)化為熱能的部分。

1.1.4 組分質(zhì)量守恒方程

(4)

式(4)中：cs、ρcs及Ds分別為組分s的體積濃度、質(zhì)量濃度及擴(kuò)散系數(shù)；Ss為系統(tǒng)內(nèi)單位時(shí)間單位體積通過(guò)化學(xué)反應(yīng)產(chǎn)生的組分s的質(zhì)量。

1.2 化學(xué)反應(yīng)

采用文獻(xiàn)[14]中給出的12組分17步化學(xué)反應(yīng)模型計(jì)算燃燒過(guò)程，以CHEMKIN嵌入計(jì)算流體力學(xué)(CFD)的流場(chǎng)計(jì)算中。參與化學(xué)反應(yīng)的12組分為CO、CO2、H、H2、H2O、OH、O、O2、HCl、Cl、Cl2、N2,反應(yīng)速率由式(5)計(jì)算：

(5)

式(5)中：κ為化學(xué)反應(yīng)速率,mol/(cm·s·k)；A為指前因子，對(duì)應(yīng)于二級(jí)化學(xué)反應(yīng)，cm3/(mol·s)；T為溫度，K；n為溫度指數(shù)1；E為活化能，cal/mol；具體反應(yīng)參數(shù)如表1所示。

表1 多步反應(yīng)Table 1 Multi-step chemical reaction

采用1～4步分支反應(yīng)和5～8步重組/離解反應(yīng)計(jì)算H2/O2系統(tǒng)，9～11步反應(yīng)計(jì)算CO/CO2系統(tǒng)，12～17步反應(yīng)描述抑制劑HCl的作用[14]。

1.3 湍流模型與燃燒模型

尾焰流場(chǎng)流動(dòng)屬于湍流形態(tài)，采用realizable，k-ε雙方程湍流模型計(jì)算可壓縮、雷諾平均N-S方程[15-16]。由于采用多步化學(xué)反應(yīng)機(jī)理，燃燒模型采用渦耗散概念模型。

1.4 等離子體濃度

根據(jù)理想氣體狀態(tài)方程，可得氣體粒子數(shù)密度為

(6)

式(6)中:ne為電離粒子的數(shù)密度,cm-3，nm為未電離粒子的數(shù)密度，cm-3;P為壓強(qiáng)，T為溫度，均采用國(guó)際單位制。

文獻(xiàn)[10]中根據(jù)薩哈-朗繆爾方程導(dǎo)出處于熱平衡狀態(tài)的氣體電離量為

(7)

式(7)中：Ui為氣體的電離能(最外層電子離開(kāi)原子所需的能量，單位為J)，K=1.38×10-23J/K為玻爾茲曼常數(shù)。

(8)

1.5 等離子體頻率

文獻(xiàn)[10]給出了等離子體頻率的計(jì)算公式：

(9)

式(9)中：ωpe為電子的等離子體角頻率；ne為電子的等離子體密度；e為電子所帶電量；ε0為真空介電常數(shù)；me為電子的質(zhì)量。在實(shí)際計(jì)算應(yīng)用中，常采用更簡(jiǎn)潔的表達(dá)方式[10]：

(10)

式(10)中：fp為等離子體角頻率對(duì)應(yīng)的頻率，Hz。

1.6 熱力學(xué)計(jì)算

研究的固體火箭發(fā)動(dòng)機(jī)使用復(fù)合推進(jìn)劑，發(fā)動(dòng)機(jī)燃燒室壓強(qiáng)為9.15 MPa，燃燒室平均質(zhì)量流量333 kg/s，固相燃燒過(guò)程采用熱力計(jì)算簡(jiǎn)化，計(jì)算結(jié)果如表2所示[13]。

表2 熱力學(xué)計(jì)算結(jié)果Table 2 Thermodynamic calculation results

1.7 網(wǎng)格構(gòu)建及邊界條件

采用CFD對(duì)固體火箭尾焰進(jìn)行仿真，以獲取火箭尾焰流場(chǎng)參數(shù)，計(jì)算區(qū)域從噴管喉部開(kāi)始，軸向長(zhǎng)度40 m，徑向長(zhǎng)度15 m，為減小仿真計(jì)算量，采用1/2對(duì)稱模型。喉部采取質(zhì)量流動(dòng)入口邊界條件，噴管壁面采取固壁邊界條件，對(duì)稱面采取對(duì)稱邊界條件，遠(yuǎn)場(chǎng)采取壓力遠(yuǎn)場(chǎng)邊界條件，網(wǎng)格構(gòu)建完成后進(jìn)行無(wú)關(guān)性驗(yàn)證，共采用15×104、44×104和62×104三套網(wǎng)格計(jì)算相同工況下的尾焰流場(chǎng)，15×104網(wǎng)格對(duì)應(yīng)的計(jì)算結(jié)果有一定誤差，44×104和62×104兩套網(wǎng)格對(duì)應(yīng)的計(jì)算結(jié)果相近，因此采用44×104網(wǎng)格方案，如圖1所示。

圖1 計(jì)算網(wǎng)格Fig.1 Computing grid

計(jì)算火箭發(fā)動(dòng)機(jī)在四種情況下的工作情況，四種工況如表3所示。

表3 發(fā)動(dòng)機(jī)工況Table 3 Engine conditions

考慮Al2O3作為離散相參與換熱，粒徑取50 μm，并假設(shè)Al2O3顆粒為球體，溫度相同，且顆粒之間無(wú)相互作用。

2 尾焰流場(chǎng)計(jì)算結(jié)果

2.1 Al2O3含量對(duì)尾焰流場(chǎng)的影響

圖2、圖3是發(fā)動(dòng)機(jī)處于地面，Al2O3含量10%、20%兩種情況下，尾焰的溫度云圖與壓強(qiáng)云圖。圖4為Al2O3含量對(duì)尾焰溫度與壓強(qiáng)沿軸線變化情況的影響圖。

圖2 工況1下的尾焰溫度及壓強(qiáng)云圖Fig.2 Contours of temperature and pressure under condition 1

圖3 工況2下的尾焰溫度及壓強(qiáng)云圖Fig.3 Contours of temperature and pressure under condition 2

圖4 不同推進(jìn)劑下尾焰溫度與壓強(qiáng)沿軸線變化情況Fig.4 Variation of temperature and pressure along the axis under different propellants

固體推進(jìn)劑中加入鋁粉可提高發(fā)動(dòng)機(jī)比沖[17]，鋁粉燃燒生成Al2O3，隨著Al2O3的增加，尾焰溫度明顯升高，推進(jìn)劑燃燒更充分，燃燒產(chǎn)物充分膨脹，噴出噴管后在外界大氣壓的作用下，形成錐形波與倒錐形波的反復(fù)出現(xiàn)，但Al2O3含量從10%增加到20%對(duì)尾焰壓強(qiáng)的影響較小。

2.2 飛行高度對(duì)尾焰流場(chǎng)的影響

圖5、圖6是Al2O3含量10%，發(fā)動(dòng)機(jī)處于不同飛行高度下，尾焰的溫度云圖與壓強(qiáng)云圖。圖7為Al2O3含量10%時(shí)，不同飛行高度對(duì)尾焰溫度與壓強(qiáng)沿軸線變化情況的影響圖。

圖5 工況3下的尾焰溫度及壓強(qiáng)云圖Fig.5 Contours of temperature and pressure under condition 3

圖6 工況4下的尾焰溫度及壓強(qiáng)云圖Fig.6 Contours of temperature and pressure under condition 4

圖7 不同飛行高度下尾焰溫度與壓強(qiáng)沿軸線變化情況Fig.7 Variation of temperature and pressure along the axis under different heights

飛行高度從10 km增加到20 km，外界大氣壓降低，對(duì)過(guò)膨脹燃?xì)獾膲嚎s作用降低，尾焰溫度峰值明顯后移，但其峰值幾乎不變，這也說(shuō)明了確實(shí)是Al2O3的含量對(duì)尾焰溫度起主要影響作用；飛行高度對(duì)尾焰壓強(qiáng)分布影響不大。

3 尾焰等離子體特征參數(shù)計(jì)算結(jié)果

尾焰等離子體濃度是表征尾焰等離子體屬性的重要物理量；當(dāng)電磁波頻率小于尾焰等離子體頻率時(shí)，電磁波將發(fā)生全反射而無(wú)法穿過(guò)尾焰[4，18-20]。計(jì)算以上等離子體特征參數(shù)，對(duì)于計(jì)算尾焰與電磁波的相互作用意義重大。

3.1 等離子體濃度

使用CFD后處理模塊的自定義函數(shù)功能，獲得尾焰等離子體濃度分布云圖如圖8所示。使用編程軟件計(jì)算不同工況下等離子體濃度隨軸線變化趨勢(shì)如圖9所示。二者吻合較好。

圖8 四種工況下尾焰等離子體濃度云圖Fig.8 Contours of plasma concentration under four conditions

圖9 不同工況下尾焰等離子體濃度沿軸線變化情況Fig.9 Variation of plasma concentration along the axis under different working conditions

隨著Al2O3含量的增加，尾焰溫度升高，燃?xì)猥@得更多的電離能，電離程度變高，等離子體濃度顯著增高，高濃度區(qū)域范圍增大，峰值位置后移；隨著發(fā)動(dòng)機(jī)工作高度增加，尾焰高溫區(qū)域后移，等離子體濃度峰值位置后移。

3.2 等離子體頻率

電磁波僅在頻率高于等離子體頻率時(shí)可以穿過(guò)等離子體[18-20]，等離子體頻率對(duì)研究尾焰對(duì)電磁波的衰減有重要作用。使用CFD后處理模塊的自定義函數(shù)功能，獲得尾焰等離子體頻率分布云圖如圖10所示。使用編程軟件計(jì)算不同工況下等離子體頻率隨軸線變化趨勢(shì)如圖11所示。二者吻合較好。

圖10 四種工況下尾焰等離子體頻率云圖Fig.10 Contours of plasma frequency under four conditions

圖11 不同工況下尾焰等離子體頻率沿軸線變化情況Fig.11 Variation of plasma frequency along the axis under different working conditions

根據(jù)簡(jiǎn)化后的等離子體頻率公式[式(10)]可知，在此種簡(jiǎn)化條件下，等離子體頻率為等離子體濃度的函數(shù)，二者正相關(guān)，因此趨勢(shì)變化相似。

4 結(jié)論

研究了不同推進(jìn)劑組成以及不同飛行高度下火箭尾焰的流場(chǎng)分布，計(jì)算了火箭尾焰等離子體濃度、等離子體頻率隨尾焰軸線方向的變化情況，給出了火箭尾焰等離子體特征參數(shù)在尾焰對(duì)稱面上的分布情況，得出以下結(jié)論。

(1)在噴管擴(kuò)張段，熱能向動(dòng)能轉(zhuǎn)化，燃?xì)馑俣壬叩耐瑫r(shí)，溫度下降，等離子體濃度降低；燃?xì)鈬姵鰢姽芎?，在激波位置，溫度迅速升高，等離子體濃度也隨之增大。

(2)尾焰等離子體特征參數(shù)高度依賴尾焰溫度，推進(jìn)劑中鋁粉含量增加，尾焰溫度顯著升高，等離子體特征參數(shù)(等離子體濃度、等離子體頻率)峰值升高；火箭飛行高度增加，尾焰高溫區(qū)域明顯后移，等離子體特征參數(shù)(等離子體濃度、等離子體頻率)峰值位置后移。

以往相關(guān)領(lǐng)域的研究，多以火箭尾焰流場(chǎng)經(jīng)驗(yàn)值為基礎(chǔ)計(jì)算尾焰等離子體參數(shù)，在此基礎(chǔ)上針對(duì)尾焰等離子體對(duì)電磁波的折射以及衰減的研究的準(zhǔn)確性難以得到保證。得到了尾焰等離子體特征參數(shù)在尾焰對(duì)稱面上的分布情況，對(duì)后續(xù)計(jì)算尾焰對(duì)電磁波的折射及衰減提供了重要參考。