李 猛，張曉宏，孫 美，王長健

（西安近代化學(xué)研究所，西安 710065）

0 引言

固體推進(jìn)劑羽流是一種含有離子和自由電子的多組分含化學(xué)反應(yīng)的高溫混合物，也是一種含有兩相流的弱電離的等離子體場(chǎng)。由于它的存在以及與周圍環(huán)境相互作用，會(huì)造成噪聲、煙霧、熱輻射、環(huán)境污染以及信號(hào)衰減等效應(yīng)，而這些效應(yīng)的研究都以固體推進(jìn)劑羽流特性研究為基礎(chǔ)。

為評(píng)估固體推進(jìn)劑羽流特性，世界上一些發(fā)達(dá)國家都建立了相關(guān)測(cè)試評(píng)估方法，建設(shè)了各種實(shí)驗(yàn)測(cè)試設(shè)施來對(duì)推進(jìn)劑羽流特性進(jìn)行檢測(cè)及表征，同時(shí)也建立了各種理論模型預(yù)示推進(jìn)劑羽流特性。國內(nèi)在對(duì)國外測(cè)試設(shè)施研究的基礎(chǔ)上，已建立起一套具有獨(dú)立知識(shí)產(chǎn)權(quán)的低成本、易控制的固體推進(jìn)劑羽流特性檢測(cè)系統(tǒng)，起草完成了相關(guān)行業(yè)標(biāo)準(zhǔn)，發(fā)表了系列的研究論文，申請(qǐng)了一些國防專利。而在羽流特性理論模型預(yù)估方面，國內(nèi)研究還不完善，尚需進(jìn)一步研究。

1 化學(xué)平衡計(jì)算模型

若系統(tǒng)內(nèi)各組分的摩爾數(shù)沒有自發(fā)的變化趨勢(shì)，稱之為達(dá)到化學(xué)平衡。在質(zhì)量守恒的約束條件下，采用吉布斯自由能最小化為化學(xué)平衡判據(jù)，建立體系的化學(xué)平衡方程為：

其中：λk是拉格朗日乘子。根據(jù)變分原理，對(duì)方程兩邊取變分得：

因?yàn)槿魏尾坏扔诹愕男×喀膎i、δλk是相互獨(dú)立的，因此可得：

根據(jù)所獲得的方程（3）和（4），在給定溫度和壓強(qiáng)的情況下，就可以求得熱力學(xué)狀態(tài)的平衡組成。NASA－CEA正是基于以上控制方程及邊界條件，采用Newton－Raphson迭代法來修正組分、拉格朗日乘子、摩爾數(shù)及溫度的初始值，將有約束問題轉(zhuǎn)化為無約束問題，并利用數(shù)學(xué)中的極值條件重新建立方程組來對(duì)問題進(jìn)行求解。

2 羽流場(chǎng)計(jì)算模型

2.1 氣相多組分輸運(yùn)模型

氣相羽流基本控制方程可以表示為以下統(tǒng)一張量形式：

式中：Γφ是對(duì)應(yīng)于φ的輸運(yùn)系數(shù)；Sφ是相應(yīng)的源項(xiàng)；φ表示流場(chǎng)通用變量，當(dāng)φ分別為l、u、v、w、T、k、ε時(shí)，上述方程分別表示連續(xù)性方程、3個(gè)坐標(biāo)方向的動(dòng)量方程、能量方程、湍流動(dòng)能和湍流耗散率方程。

羽流多組分輸運(yùn)方程為：

式中：Ri為與化學(xué)反應(yīng)相關(guān)的第i種組分的凈生成率，即單位體積的質(zhì)量生成率；Si為離散相及用戶定義源項(xiàng)所額外產(chǎn)生的質(zhì)量生成率。Vi為第i種組分的擴(kuò)散速度。

針對(duì)推進(jìn)劑羽流特性評(píng)估的特點(diǎn)，采用二維軸對(duì)稱簡化模型及非穩(wěn)態(tài)歐拉方程求解，湍流模型采用兩方程的修正的k－ε模型（Realizable k－εturbulent model）。

2.2 有限速率化學(xué)反應(yīng)模型

處理化學(xué)反應(yīng)采用有限速率的化學(xué)反應(yīng)模型，其在本問題中的適用性已由Rodionov等人進(jìn)行了驗(yàn)證[3]，反應(yīng)速度系數(shù)是溫度的強(qiáng)烈非線性函數(shù)，對(duì)具有N個(gè)基元反應(yīng)的某反應(yīng)，其當(dāng)量表達(dá)式可以寫為：

其中：v′I、v″I分別為基元反應(yīng)中反應(yīng)物和生成物的當(dāng)量反應(yīng)系數(shù)；用Arrhenius定律表示正、逆反應(yīng)速率常數(shù)如下式所示：

其中：E為活化能；R為通用氣體常數(shù)；A稱作指數(shù)前因子或頻率因子。從上式可以看出，反應(yīng)溫度T是決定反應(yīng)速率的主要因素，但當(dāng)流動(dòng)存在湍流脈動(dòng)時(shí)，反應(yīng)速率不但受化學(xué)反應(yīng)動(dòng)力學(xué)影響，也受到湍流脈動(dòng)的影響，因此可以根據(jù)化學(xué)反應(yīng)動(dòng)力學(xué)和湍流脈動(dòng)性質(zhì)來確定組分的化學(xué)反應(yīng)速率。

2.3 離散相模型

采用拉格朗日兩相流模型對(duì)混合燃?xì)?、三氧化二鋁等顆粒兩相流在噴管及排氣羽流場(chǎng)內(nèi)的流動(dòng)進(jìn)行仿真。其實(shí)質(zhì)是用顆粒間的動(dòng)量交換來模擬顆粒間的碰撞過程，是一種重要的兩相流模型，對(duì)顆粒相采用單顆粒尺度上的跟蹤描述，而對(duì)氣相采用連續(xù)介質(zhì)假設(shè)模型，即對(duì)模型中氣相流場(chǎng)中的計(jì)算采用歐拉模型，而固相顆粒作為離散相，其動(dòng)量守恒方程為：

其中：Fdr為流體拖曳力；Fp為壓強(qiáng)梯度力；Fam為虛擬質(zhì)量力，F(xiàn)b為體積力。

3 計(jì)算實(shí)例及結(jié)果分析

3.1 推進(jìn)劑羽流計(jì)算前處理

3.1.1 化學(xué)平衡計(jì)算及化學(xué)反應(yīng)動(dòng)力學(xué)選取

基礎(chǔ)配方（質(zhì)量分?jǐn)?shù)）：NC＋NG（45%～60%）；RDX（20%～35%）；AL（0～5%）；燃燒催化劑（3.5%～4.5%）；其它助劑（4.5%～6.5%）。

應(yīng)用NASA－CEA計(jì)算7MPa下某改性雙基推進(jìn)劑化學(xué)平衡性能后，噴管入口主要參數(shù)以及主要組分質(zhì)量百分含量如表1所示。根據(jù)所得組分選取化學(xué)反應(yīng)模型為9組分10反應(yīng)的化學(xué)動(dòng)力學(xué)模型，所用反應(yīng)機(jī)理數(shù)據(jù)如表2所示。

表1 噴管入口參數(shù)及組分質(zhì)量分?jǐn)?shù)

表2 化學(xué)反應(yīng)模型

3.1.2 計(jì)算區(qū)域及網(wǎng)格結(jié)構(gòu)

選取燃燒室末端為燃?xì)馊肟?，將噴管流?chǎng)和羽流場(chǎng)進(jìn)行一體化處理，計(jì)算區(qū)域如圖1所示，其中區(qū)域1為噴管流動(dòng)區(qū)域，區(qū)域2和區(qū)域3為羽流流動(dòng)區(qū)域，對(duì)稱半?yún)^(qū)域尺寸為x方向6m，y方向2m。噴管內(nèi)區(qū)域靠近壁面處進(jìn)行了網(wǎng)格加密，羽流區(qū)域在靠近軸線處進(jìn)行了網(wǎng)格加密，網(wǎng)格數(shù)量在兩萬左右（見圖2）。在羽流計(jì)算中，噴管長度為30mm；羽流部分長6m，寬3m。

圖1 羽流計(jì)算區(qū)域分布簡圖

圖2 羽流計(jì)算網(wǎng)格圖

3.2 推進(jìn)劑羽流計(jì)算后處理

3.2.1 化學(xué)反應(yīng)動(dòng)力學(xué)影響分析

圖3顯示出了考慮化學(xué)反應(yīng)模型前后，羽流場(chǎng)各組分質(zhì)量分?jǐn)?shù)分布圖形，當(dāng)考慮羽流為氣相多組分流場(chǎng)而不考慮羽流的化學(xué)反應(yīng)時(shí)，羽流場(chǎng)中各組分具有相同的流場(chǎng)分布，如圖3（a）所示。考慮化學(xué)反應(yīng)動(dòng)力學(xué)后，如圖3（b）～圖3（d），組分分布圖產(chǎn)生了明顯變化，一些原無或很少的組分（OH、H、O）隨著化學(xué)反應(yīng)的進(jìn)行，這些物質(zhì)濃度逐漸增加并參加化學(xué)反應(yīng)，對(duì)比是否考慮化學(xué)模型組分分布圖可以看出，二次燃燒主要發(fā)生在燃?xì)饬髋c空氣混合的邊界區(qū)域，原因是二次燃燒的化學(xué)反應(yīng)產(chǎn)生熱量，各組分之間發(fā)生了質(zhì)量和能量的擴(kuò)散。二次燃燒主要發(fā)生在尾焰的外邊界區(qū)域和距離噴管較遠(yuǎn)的尾焰區(qū)，而在尾焰核心區(qū)內(nèi)，二次燃燒基本不發(fā)生，隨著尾焰核心區(qū)的終止，外流空氣與燃?xì)饣旌?，化學(xué)反應(yīng)加劇。

3.2.2 離散相影響分析

圖4中的x軸與發(fā)動(dòng)機(jī)軸線重合，y軸代表流場(chǎng)中速度，取噴管入口中心點(diǎn)為坐標(biāo)原點(diǎn)。

考慮離散相模型對(duì)羽流的影響，在兩相流場(chǎng)中，由于相間作用，顆粒一方面受燃?xì)庾饔米冯S燃?xì)膺\(yùn)動(dòng)，一方面受慣性作用保持原有運(yùn)動(dòng)狀態(tài)，因此在燃?xì)馍淞鞯膬上嗔鲌?chǎng)中，顆粒運(yùn)動(dòng)對(duì)燃?xì)獾臏笞饔煤苊黠@。顆粒相對(duì)燃?xì)膺\(yùn)動(dòng)的滯后，使得顆粒與燃?xì)庵g存在速度梯度和溫度梯度，產(chǎn)生相間傳熱，從而影響兩相流情況下的燃?xì)馍淞髁鲌?chǎng)。與純氣相流場(chǎng)相比，兩相流情況下燃?xì)馑俣鹊戎稻€會(huì)出現(xiàn)滯后。

圖3 羽流組分質(zhì)量分?jǐn)?shù)分布圖

圖4 速度分布對(duì)比曲線

4 結(jié)論

文中建立了改性雙基推進(jìn)劑化學(xué)平衡性能計(jì)算模型，使用有限速率化學(xué)反應(yīng)模型、離散相模型對(duì)改性雙基推進(jìn)劑羽流場(chǎng)進(jìn)行了計(jì)算，獲得了羽流結(jié)構(gòu)及其組分分布情況。并對(duì)考慮有限速率化學(xué)反應(yīng)模型前后羽流組分分布情況進(jìn)行了研究，結(jié)果表明不考慮化學(xué)反應(yīng)模型時(shí)，羽流場(chǎng)中各組分具有相同的流場(chǎng)分布，而在考慮化學(xué)反應(yīng)模型后，組分分布發(fā)生了明顯變化，很重要的火箭發(fā)動(dòng)機(jī)二次燃燒現(xiàn)象發(fā)生。進(jìn)一步對(duì)比分析了考慮離散相模型前后羽流軸線上的速度分布情況，結(jié)果顯示速度數(shù)值差別很大，顆粒運(yùn)動(dòng)對(duì)燃?xì)獾臏笮?yīng)很明顯。因此，在改性雙基推進(jìn)劑羽流計(jì)算中，需要考慮有限速率化學(xué)反應(yīng)模型和離散相模型的共同影響。

[1]Sanford Gordon，Bonnie J McBride.Computer program for calculation chemical eguilibrium compositions and applications：I Analysis，NASA RP－1311[R].1994.

[2]Bonnie J MeBride，Sanford Gordon.Computer program for calculation chemical eguilibrium compositions and applications：ⅡUsers Mawual and profram description，NASA－RP－1311[R].1996.

[3]J Troyes，I Dubois，V Borie，et al.Multi－phase reactive numerical simulations of a model solid rocket motor exhaust jet，AIAA2006－4414[J].2006.

[4]Xiao－Yen J Wang，James R Yuko.Thermal analysis on plume heating of the main engine on the crew exploration vehicle service module，NASA／TM－2007－215049[R].2007.

[5]NATO Research and Technology Organisation.Terminology and assessment methods of solid propellant rocket exhaust signature，AGARD－AR－287[R].1993.

[6]徐義華，胡春波，張勝敏，等.固體火箭發(fā)動(dòng)機(jī)羽流紅外輻射特性研究[J].固體火箭技術(shù)，2010，33（2）：176－181.

[7]李猛，王宏.火箭發(fā)動(dòng)機(jī)排氣羽流特征信號(hào)預(yù)示技術(shù)研究綜述[J].飛航導(dǎo)彈，2010（3）：54－57.

[8]李猛，王宏，陳雪莉.復(fù)雜化學(xué)平衡應(yīng)用計(jì)算程序[J].四川兵工學(xué)報(bào)，2010，31（9）：132－134.

[9]王偉臣，魏志軍，等.后燃對(duì)火箭發(fā)動(dòng)機(jī)羽流紅外特性的影響[J].航空動(dòng)力學(xué)報(bào)，2010，25（11）：2612－2618.