崔立堃，葉 偉

(陜西理工大學 機械工程學院， 陜西 漢中 723001)

固體火箭發(fā)動機羽流流場的結構分布特性對固體火箭的總體設計、火箭發(fā)射裝置設計、結構材料的選擇等方面均會產生很大的影響[1-2]，特別是其紅外輻射對發(fā)動機故障監(jiān)控和軍事偵察有著非常重要的作用。為此，國內外許多學者對此進行了研究：H.B.Ebrahimi等[3]運用GPACT模型計算分析了Titan II 液體火箭發(fā)動機紅外波段輻射場；K.Berer等[4]比較了運用逐線模型和分子波帶模型計算彈道導彈羽流紅外輻射的結果；Kim等使用改進的離散坐標法[5]和有限體積法[6]計算了固體發(fā)動機羽流輻射和噴管探照燈(searchlight)發(fā)射的底部加熱效應；張小英等[7]研究了某型固體火箭發(fā)動機在地面試車和不同高度工況下發(fā)動機紅外波段2～6 μm的輻射情況；金云飛等[8]基于羽流輻射的監(jiān)測和光譜特性，對發(fā)動機的工作狀態(tài)、羽流特性等進行了在軌觀測試驗，試驗結果良好；陳杰等[9]采用差分求解N-S方程與直接模擬蒙特卡洛耦合的方法對火星環(huán)繞器姿軌控發(fā)動機羽流力、熱以及污染效應進行了仿真研究。這些研究主要集中在羽流紅外輻射的計算分析方法和特性探究方面，而對羽流場的流場結構及其特性的研究較少。而羽流場的流場結構及特性對包括羽流紅外輻射等其他羽流效應有著重要的意義。因此，本文基于三維N-S方程，采用Spalart-Allmaras湍流模型，通過大型CFD軟件FLUENT對某固體火箭發(fā)動機羽流場進行了數(shù)值分析，并對不同飛行參數(shù)條件下、純氣相作用下的發(fā)動機羽流場進行了綜合分析比較，擬為深入研究固體火箭發(fā)動機羽流效應提供一定的指導。

1 數(shù)值計算模型

可壓縮流動的強守恒型N-S方程在笛卡爾坐標系中可以表示為：

(1)

式(1)中各系數(shù)的具體含義及表達式見文獻[10]。

在動量方程中，涉及到湍流粘性系數(shù)μ1，它必須通過湍流模型求得，本文采用工程上廣泛使用的Spalart-Allmaras模型，該模型是一個相對簡單的模型，只求解一個有關渦粘性的輸運方程，計算量相對較小。此模型是專門應用于涉及到束縛壁面流動的航空領域，對于受到反壓梯度作用的邊界層，能給出很好的模擬結果。湍流粘性系數(shù)μ1按下式計算：

(2)

式(2)中，各系數(shù)的具體含義及表達式見文獻[11]。

1.1 幾何模型及網(wǎng)格劃分

本文以某艦載導彈固體火箭發(fā)動機為例，選定發(fā)動機幾何尺寸。整個彈體結構簡圖及其網(wǎng)格劃分如圖1所示。發(fā)動機噴管主要幾何尺寸參數(shù)見表1。

表1 噴管尺寸參數(shù)

高質量的網(wǎng)格劃分能夠提高計算精度、加快計算速度，加大計算穩(wěn)定性。本文中計算域形狀不規(guī)則，因此為了提高網(wǎng)格質量將其分割為10個部分，

分別對其進行結構化網(wǎng)格劃分。由于粘性計算和激波捕捉的需要對噴管壁面附近和尾流場軸線附近進行了網(wǎng)格加密。總的網(wǎng)格數(shù)為359 577個。網(wǎng)格劃分情況如圖2所示。

模擬介質為某推進劑燃燒產物的成分，推進劑按照鋁18.2%、過氯酸銨67.2%、丁羥11%、葵二酸二辛酯3.1%、甲苯二異氰酸酯0.5%配比，噴管擴張比是9.5，燃燒組分的計算采用布林克萊法，詳細的計算公式和方法及程序參見文獻[11]。表2中列出了計算后得到的主要組分及其摩爾數(shù)。

表2 燃氣中各組分含量

流場計算中取所占比重最大的前七種組份，分別為：CO、HCL、H2、N2、H2O、CO2、H，這七種物質占氣態(tài)產物總量的98.461%。燃燒產物中Al2O3雖然質量占比比較大，但它是凝相，在仿真發(fā)動機羽流場波系結構時可忽略。

1.2 邊界條件

本文數(shù)值計算中涉及到的邊界條件有壓力入口邊界、壁面邊界、對稱邊界和壓力遠場邊界。其中，噴管入口采用壓力入口邊界條件，總壓5 MPa，總溫3 400 K，流動方向與入口邊界垂直；噴管壁面采用絕熱無滑移固體壁面邊界條件；羽流場徑向截面采用對稱邊界條件。外邊界采用壓力遠場邊界條件，壓力與環(huán)境壓力一致，采用標準大氣相應高度下的參數(shù)(見表3)。

很多人心里都不舒服，說來說去的很長一段時間，之后又都閉了口。試問一下，別人除過自己的工作，去掏廁所把自己弄得臭烘烘連老婆孩子都不愿近身了？去牽著驢頂著毒日頭口干舌燥滿山溝地給它找吃的去了？還是在假期本該休息時一個人寡兮兮地在學校修理那些被淘氣孩子損壞的課桌？還是給同事和村里人的家具義務上油漆，給學校充大工修理房子，當農工苦嗆嗆地去墾荒田了？

表3 大氣壓強和高度的關系

2 計算結果分析

2.1 相同飛行高度和攻角，不同飛行馬赫數(shù)羽流場數(shù)值分析

選取飛行高度為0 km，攻角為0°，分別對飛行馬赫數(shù)為0.6、0.9、1.2、1.5這4種情況下的發(fā)動機羽流場進行數(shù)值計算。其羽流場馬赫數(shù)等值線圖如圖3所示。

由圖3可以得出，在地面狀態(tài)時不同飛行馬赫數(shù)下羽流流場結構分布的變化規(guī)律：當飛行馬赫數(shù)為0.6時，燃氣在噴管后方羽流區(qū)域形成了五個明顯的膨脹-壓縮波系，每個波系的長度約為0.8 m，波系的強度逐次降低；當飛行馬赫數(shù)增加到0.9時，噴管后方羽流區(qū)域仍然有五個明顯的膨脹-壓縮波系，波系的膨脹半徑略有減小，波系的長度、強度等方面變化不大；當飛行馬赫數(shù)增加到1.2時，噴管后方羽流區(qū)域形成了四個明顯的膨脹-壓縮波系，波系的膨脹半徑明顯減小，波系的長度略有增加，約為0.9 m，羽流核心區(qū)的長度有所增加；當飛行馬赫數(shù)增加到1.5時，流場結構發(fā)生了明顯的變化，噴管后方羽流區(qū)域只有一個明顯的膨脹-壓縮波系，波系膨脹半徑較亞聲速時明顯降低，第二個和第三個波系未能充分發(fā)展，羽流核心區(qū)的長度繼續(xù)增加。

2.2 相同飛行馬赫數(shù)和攻角，不同飛行高度羽流場數(shù)值分析

選取飛行馬赫數(shù)Ma為0.9，攻角α為0°，分別對飛行高度H為0 km、5 km、15 km這3種情況下的發(fā)動機羽流場進行數(shù)值計算，結果如圖4所示。

由圖4可得出，在亞音速狀態(tài)時，不同飛行高度下羽流流場結構分布的變化規(guī)律如下：

當飛行高度為0 km時，燃氣在噴管后方羽流區(qū)域形成了五個明顯的膨脹-壓縮波系，每個波系的長度約為0.8 m，波系的強度逐次降低；當飛行高度為5 km時，噴管后方羽流區(qū)域形成了四個明顯的膨脹-壓縮波系，波系的膨脹半徑有所增大，波系的長度有所增大，約為1.4 m，波系的強度逐次降低，羽流場邊界和羽流核心區(qū)的長度都有所增加；當飛行高度為15 km時，噴管后方羽流區(qū)域有兩個明顯的膨脹-壓縮波系，第三個波系未能充分發(fā)展，波系的膨脹半徑和波的長度明顯增大，羽流場邊界和羽流核心區(qū)的長度較地面、中低空狀態(tài)時明顯增加。

圖5、圖6和圖7分別為Ma=0.9、攻角α=0°，飛行高度H=0為0 km、5 km、15 km時羽流場軸線上馬赫數(shù)、溫度和壓強的分布對比。

由圖5～7可以得出，由于膨脹-壓縮波系的存在，在不同飛行高度下軸線上的馬赫數(shù)、溫度和壓強這些流場結構參數(shù)都隨著膨脹-壓縮波系個數(shù)、長度和強度的不同而發(fā)生著較為明顯的變化。

2.3 相同飛行高度和飛行馬赫數(shù)，不同攻角羽流場數(shù)值分析

圖8為Ma=0.9、H=0 km，攻角α分別為0°、4°、8°、12°時的羽流場馬赫數(shù)等值線圖。

由圖8可以得出，亞音速狀態(tài)時，不同攻角下羽流流場結構分布的變化規(guī)律如下：

當攻角為0°時，燃氣在噴管后方羽流區(qū)域形成了五個明顯的膨脹-壓縮波系，每個波系的長度約為0.8 m，波系的強度逐次降低；當攻角為4°時，噴管后方羽流區(qū)域仍然有五個明顯的膨脹-壓縮波系，波系的膨脹半徑、長度和強度的變化不大，可以明顯的看出膨脹-壓縮波系形成上揚趨勢，不再具有軸對稱性，羽流核心區(qū)和羽流場邊界偏離軸線；當攻角為8°時，噴管后方羽流區(qū)域仍然有五個明顯的膨脹-壓縮波系，波系的膨脹半徑、長度仍有所減小，膨脹-壓縮波系偏離軸線上揚角度加大，羽流核心區(qū)長度和羽流場邊界減??；當攻角為12°時，噴管后方羽流區(qū)域仍然有五個膨脹-壓縮波系，波系的膨脹半徑、長度繼續(xù)減小，膨脹-壓縮波系較攻角為0°時與軸線形成明顯的上揚角度，羽流核心區(qū)長度繼續(xù)減小，在第四個膨脹-壓縮波系后，羽流核心區(qū)和羽流場邊界完全偏離軸線。

圖9為Ma=0.9、H=0 km，攻角α分別為0°、4°、8°、12°，與噴管出口的距離為0 m、2 m、4 m、6 m時的羽流場截面溫度云圖。

由圖9可以得出，亞音速狀態(tài)時，不同攻角情況下，與噴管出口距離不同的羽流場截面溫度分布的變化規(guī)律如下：

噴管出口截面，攻角為0°時溫度為軸對稱分布，隨著攻角的增大，截面上溫度的分布產生變化，攻角為12°時，可以明顯看到溫度分布不再具有軸對稱性；距離噴管出口2 m處的截面，攻角為0°時溫度仍為軸對稱分布，隨著攻角的增大，可以明顯看到溫度分布偏離中心，向Z軸方向移動，產生上揚；距離噴管出口4 m處的截面，攻角為0°時溫度仍為軸對稱分布，隨著攻角的增大，可以明顯看到溫度分布向Z軸方向移動上揚的距離增大；距離噴管出口6 m處的截面，攻角為0°時溫度仍為軸對稱分布，由于羽流場膨脹-壓縮波系不在形成，此時隨著攻角的增大，溫度分布上揚程度較羽流場過渡段更加明顯，上揚程度加強，攻角為12°時，可以明顯看到溫度分布完全離開截面中心。通過對亞音速狀態(tài)，不同攻角情況下，與噴管出口距離不同截面上溫度分布的觀察，我們進一步驗證了攻角能夠使膨脹-壓縮波系、羽流核心區(qū)偏離軸線，產生上揚，隨著攻角的增大，上揚的角度增大，羽流核心區(qū)后端會完全偏離軸線。

3 結論

1) 隨著飛行馬赫數(shù)的增加，羽流場膨脹壓縮波系的個數(shù)逐漸減少，膨脹壓縮波系的膨脹半徑減小，羽流核心區(qū)的長度逐漸增大。

2) 隨著飛行高度的增加，大氣壓力逐漸減小，噴管出口壓力比增加，導致噴管出口膨脹后的普朗特-梅耶角增大，羽流受到外界的作用變小，能夠膨脹得更完全，羽流場膨脹壓縮波系的個數(shù)逐漸減少，膨脹壓縮波系的膨脹半徑增大，膨脹壓縮波系的長度增加，羽流場邊界和羽流核心區(qū)的長度都逐漸增大。

3) 隨著攻角的增加，羽流場膨脹壓縮波系偏離軸線上揚，不再具有軸對稱性，膨脹壓縮波系的上揚角度逐漸增大，膨脹壓縮波系的個數(shù)逐漸減少，膨脹壓縮波系的膨脹半徑和長度減小，羽流場邊界和尾流核心區(qū)的長度逐漸減小。