液體火箭發(fā)動機噴管仿真模型

孫得川，楊建文

(1.大連理工大學 航空航天學院，遼寧 大連 116024；2.西安航天動力研究所，陜西 西安 710100)

0 引言

化學火箭發(fā)動機使用拉瓦爾噴管作為能量轉換的部件，用以產生推力，對其流場進行仿真分析是現代噴管優(yōu)化設計和理解各種復雜流動過程必須經歷的環(huán)節(jié)和步驟。隨著計算流體力學(CFD)的發(fā)展，對噴管流動的仿真似乎已經不成問題，流行的CFD軟件幾乎都可以“輕松地”得到計算結果；但是得到結果容易，得到更接近真實值的結果卻并非易事。尤其是現代噴管設計要結合具體的發(fā)動機條件，對性能仿真的準確性提出了更高的要求，以滿足優(yōu)化設計的目標。

噴管性能計算的發(fā)展歷史主要以美國為代表。NASA早在20世紀60年代就開始研究發(fā)動機性能計算的程序，最初基于化學平衡假設，采用最小吉布斯自由能法來計算燃燒室和噴管的化學平衡狀態(tài)，也可以計算透平機械、激波管等流動中的化學平衡問題。目前該程序的化學平衡部分代碼已經基本公開，即近年來廣為流傳的CEA軟件。在該程序的基礎上，NASA主持開發(fā)了液體火箭發(fā)動機的二維化學動力學計算軟件TDK(two-dimensional kinetics)并多次更新，可處理諸如三組元發(fā)動機、超燃沖壓發(fā)動機中的化學反應流動問題。TDK軟件預測的比沖與試車相比偏差在3%以內，已經成為美國JANNAF的標準軟件，并且共享給其盟友使用，在液體發(fā)動機研制中發(fā)揮了重要的作用。

在眾多文獻中，Manski等對液體火箭發(fā)動機性能的研究比較具有代表性，他們利用TDK詳細計算研究了SSME和Vulcain發(fā)動機的性能，重點討論了室壓對噴管效率的影響。其研究指出：對于氫氧一級發(fā)動機，當室壓較低時噴管效率約為97%，而當室壓提高到10 MPa時噴管效率可達98%以上；由噴管二維效應引起的擴張損失是總損失的主要部分，隨著室壓提高，擴張損失從初始的0.6%上升到1.5%；SSME的室壓超過10 MPa后，化學動力學損失幾乎為零，但是在低壓情況下，化學動力學損失可達1%。可見，室壓提高時噴管效率提高的主要貢獻在于化學動力學損失減少。

與西方國家相比，國內在噴管性能的研究方面還比較欠缺，這主要是因為該研究內容偏于“老舊”，似乎缺乏創(chuàng)新。因此對于噴管流動和各種損失的理解還較為淺薄，這不利于建立在仿真基礎上的發(fā)動機優(yōu)化設計。

仿真是通過數學方法模擬真實世界，其“保真”有兩個環(huán)節(jié)，首先是數學模型是否真正反映了物理過程，其次是求解數學模型的過程是否準確。筆者認為目前計算技術的發(fā)展已經可以保證第二個環(huán)節(jié)正確，噴管性能仿真準確度的提高主要取決于所采用的數學模型是否恰當。因此，本文以典型的液體火箭發(fā)動機噴管為例討論氣體模型、化學反應模型等對計算結果的影響，并進一步厘清影響噴管性能的物理因素，以期對發(fā)動機設計有所啟發(fā)和幫助。

1 噴管性能的仿真模型

1.1 理想噴管

噴管性能以理想噴管的性能作為參照。理想噴管是出口壓力等于環(huán)境壓力的一維噴管，能夠產生設計條件下的最大推力，其推力系數表示為

(1)

式(1)完全是根據氣體等熵膨脹假設得出，可見理想推力系數僅與燃氣的比熱比和壓比/相關。圖1以壓比為變量，給出了不同比熱比所對應的曲線。顯然，壓比增大會得到較大的推力系數，即(在不發(fā)生流動分離的情況下)面積比大的噴管推力系數大。而在同樣的壓比條件下，如果燃氣比熱比較小，也可以獲得較大的推力系數。

圖1 理想噴管的推力系數與壓比的關系Fig.1 Relationship between thrust coefficient and pressure ratio forideal nozzle

因為燃氣的比熱比與推進劑相關，并隨化學反應流動而變化(燃氣在噴管中流動時，由于溫度下降和化學組分變化，比熱比也隨之變化，一般隨面積比或壓比增大而增大)，所以工程應用中一般以化學平衡流動所得到的比熱比作為理想值。圖2是采用一維化學平衡計算得到的某空間發(fā)動機噴管中燃氣的比熱比隨面積比的變化。

圖2 假設流動處于不同狀態(tài)下，噴管中燃氣比熱比和主要組分與面積比的關系Fig.2 Relationship of specific heat ratio and main components to area ratio under different flow conditions

該發(fā)動機采用四氧化二氮/一甲基肼作為推進劑，混合比為1.65，室壓為0.85 MPa?？梢钥吹饺艏僭O燃氣處于化學平衡狀態(tài)，則比熱比介于1.24～1.26之間；而若假定燃氣凍結，則由于燃氣組分不再變化，其比熱比隨燃氣膨脹、降溫而增大得較多；若考慮多組分處于化學非平衡，則得到的比熱比略低于凍結流假設計算值。對照主要組分的變化，可見在化學平衡計算中，H和CO逐漸增多而CO和HO逐漸減少，但是在實際的反應流中，組分只在膨脹初期(面積比小于10)有明顯變化，而后基本凍結。所以，噴管擴張段中的燃氣更接近凍結流，其比熱比大于平衡流的值。而由圖1可知較大的比熱比所對應的推力系數較小，故而不容易提高噴管的性能。

1.2 氣體模型

由以上分析可知，比熱比對準確預示噴管性能起決定性作用，而比熱比的計算首先由所采用的氣體模型決定。表1對比了噴管流場仿真中所使用的氣體模型的特點。

表1 不同氣體模型的特點Tab.1 Characteristics of different gas models

各模型的狀態(tài)方程為

=

(2)

式中為氣體常數，且=-?？梢妼τ诹繜嵬耆珰怏w假設，其定壓比熱和比熱比均為常數，而比熱或比熱比的確定則完全靠經驗；對于噴管流動計算，常采用喉部位置的化學平衡計算值。

對于熱完全混合氣體假設，其中各組分的摩爾定壓比熱的計算常采用由JANNAF表擬合的多項式，即

(3)

式中：為通用氣體常數；為組分序號。混合氣的摩爾定壓比熱為

(4)

式中：為組分數量；為組分的摩爾分數。因絕大多數氣體單質的定壓比熱都是隨溫度升高而增大的，所以混合氣體的比熱也具有同樣的特點。

若假設氣體組分凍結，則燃氣成分的含量不變，故定壓比熱只與溫度有關，一般隨溫度降低而減小。若假設氣體組分處于化學平衡狀態(tài)，則組分種類和含量都隨溫度瞬態(tài)變化、而與反應過程無關，通常用吉布斯自由能最小的方法求解組分及其含量。在真實的噴管流動中，燃氣組分處于化學非平衡狀態(tài)，必須采用化學動力學的方法求解。

1.3 二維非平衡流動的控制方程

二維非平衡流動的控制方程為軸對稱N-S方程，各組分的質量守恒方程為

(5)

式中：為組分的凈生成率；坐標和都以喉部半徑進行無量綱化處理。對所有組分求和，得到混合燃氣的連續(xù)方程為

(6)

動量方程為

(7)

能量方程為

(8)

其中

(9)

1.4 化學動力學模型

設混合氣體中含有種組分，描述化學反應過程的機理含有個反應，反應機理表達為

(10)

式中:為反應序號;表示組分名稱；和′分別為各組分在反應前、后的計量系數。

組分的凈生成率為

(11)

(12)

(13)

式中:為指前因子;為溫度系數;a為活化能由反應機理給出。

1.5 湍流模型

對于沒有流動分離的噴管流動，湍流模型對計算結果的影響很小。本文采用經典的-兩方程模型，近壁處采用壁面函數法處理，在此不贅述。

1.6 控制方程求解

因噴管結構相對簡單，可采用結構網格，故在計算程序中對空間差分采用了適用于結構網格的三階WENO格式，對時間推進采用三階精度的Runge-Kutta格式。

對于化學反應源項，有

(14)

當只考慮反應引起的質量變化時，有

(15)

其中

(16)

采用如下半隱式梯形公式求解

(17)

邊界條件按照如下方法給出：入口邊界指定流量，其溫度和化學成分組成由化學平衡計算確定；壁面邊界為無滑移條件，壓強、溫度、組分濃度梯度為零；軸線邊界只有軸向速度，壓強、溫度、組分濃度梯度為零；出口邊界為無反射條件。

2 典型噴管的性能仿真分析

2.1 噴管參數及計算條件

本文針對3個發(fā)動機噴管進行了仿真研究，其參數見表2，其中針對噴管Ⅱ還計算了室壓為22.0 MPa的情況，以評估室壓對化學動力學的影響。

表2 噴管參數Tab.2 Nozzle parameters

圖3給出了噴管Ⅰ的計算網格，其軸向節(jié)點數為151，徑向節(jié)點數為51，在喉部和壁面附近加密。其中壁面第一層網格的≈30～50，與湍流模型適應，不影響摩擦損失的計算。其他兩噴管的網格類似。

圖3 噴管Ⅰ的計算網格Fig.3 Calculation grid of nozzle Ⅰ

噴管入口溫度、組分及其含量由化學平衡計算給出，見表3。采用的反應機理見表4和表5。

表3 噴管入口溫度、組分及質量分數Tab.3 Temperature， component and mass fraction at nozzle inlet

表4 CHON系統(tǒng)的反應速率Tab.4 Reaction rate of CHON system 單位:cc，K，mole，sec

表5 CHON反應的三體系數Tab.5 Coefficients of the third body for CHON reaction

2.2 氣體模型的影響

首先針對噴管Ⅰ用不同氣體模型進行性能計算。采用量熱完全氣體假設時，燃氣的比熱值設為2 156.6 J/(kg·K)(該值為燃燒室化學平衡計算結果，對應比熱比為1.232 7)，對應的溫度為3 042 K。圖4對比了按不同氣體模型所計算的噴管溫度場。

圖4 噴管Ⅰ內的溫度分布Fig.4 Temperature distribution in nozzle Ⅰ

可以看到采用量熱完全氣體假設所計算的噴管出口最低溫度約為500 K，而采用凍結的熱完全氣體假設所計算的出口溫度則明顯偏低，非平衡氣體的出口溫度則介于兩者之間。另外在緊靠噴管喉部的初始膨脹段，采用量熱完全氣體假設和熱完全氣體假設計算得到的溫度等值線接近，呈更加外凸的形式，而化學非平衡計算的等值線曲率則較小。相應地，因采用量熱完全氣體計算的噴管出口溫度較高，所以對應的馬赫數較低(見圖5)。

圖5 噴管Ⅰ內的馬赫數分布Fig.5 Mach number in nozzle Ⅰ

3種氣體模型所計算的軸向速度分布見圖6，其中凍結的熱完全氣體所對應的出口速度明顯低于其他兩種情況。另外，從溫度和馬赫數分布可以觀察到當采用量熱完全氣體和凍結的完全氣體假設時，噴管內從擴張段起始的壓縮波較明顯，等值線的折轉較突然；而實際的非平衡氣體流動中，等值線折轉比較圓滑，說明組分的變化起耗散作用，而且這種影響起源于擴張段的起始部分。

圖6 噴管Ⅰ內的軸向速度分布Fig.6 Axial velocity distribution in nozzle Ⅰ

圖7顯示了反應中自由基團OH的分布，因OH是反應中的活性組分，其濃度可以反映出反應的活躍程度，所以該圖表示喉部附近的流動為非平衡流動，在面積比較大的噴管下游才能視為凍結流動。

圖7 噴管Ⅰ內OH組分的質量分數Fig.7 Mass fraction of OH in nozzle Ⅰ

表6列出了采用不同氣體模型計算的噴管Ⅰ的真空比沖與理論值和試車測量值的對比。比較可知，二維化學動力學計算得到的真空比沖與發(fā)動機實測值最接近，除燃燒效率之外的綜合效率為94%；若考慮實際燃燒效率大于98%，則采用化學動力學計算的發(fā)動機比沖將落在試車值范圍內，這充分說明了化學動力學計算的準確性；而采用量熱完全氣體和熱完全氣體假設所計算的效率都偏高。

表6 計算得到的噴管Ⅰ的性能Tab.6 Calculation performance of nozzle Ⅰ

為了更進一步說明氣體模型對計算結果帶來的影響，表7給出了采用不同模型計算的噴管Ⅱ的比沖。

比較表6和表7可知，采用量熱完全氣體假設所計算的比沖高于用熱完全氣體假設的計算值。這是因為量熱完全氣體假設比熱比不變，而熱完全氣體的比熱比隨溫度降低而增大，這符合式(1)(圖1)對于比熱比對噴管性能影響的說明。但是在噴管Ⅰ的計算中，采用量熱完全氣體或熱完全氣體假設所計算的比沖比二維化學動力學計算的比沖高，而在噴管Ⅱ的計算中恰恰相反。造成這種差異的原因可能有二，其一是兩種發(fā)動機所采用的推進劑組合不同；其二正是化學動力學損失的影響，噴管Ⅰ的室壓較低，化學動力學損失較大，而噴管Ⅱ的室壓較高，化學動力學損失較小。

表7 計算得到的噴管Ⅱ(pc=17.7 MPa)的性能Tab.7 Calculation performance of nozzle Ⅱ(pc=17.7 MPa)

目前噴管性能的CFD計算中，很多都采用量熱完全氣體或熱完全氣體假設，這顯然會造成較大的偏差。

2.3 室壓對化學動力學的影響

理論上講，提高燃燒室壓強有助于提高發(fā)動機性能。這主要有兩方面的原因：其一是提高室壓會提高燃燒效率，其二是高室壓下噴管的化學動力學損失小。

為了說明第二點，對比了噴管Ⅱ在17.7 MPa和22.0 MPa室壓下的化學非平衡流場。計算時給定兩個算例的入口組分含量完全相同，相當于燃燒效率相同。計算得到當室壓提高到22.0 MPa時，比沖為329.8 s，比表7給出的室壓為17.7 MPa時的比沖328.8 s高了1.0 s。

為了解室壓增大提高噴管性能的原因，圖8對比了用兩種室壓條件計算的噴管Ⅱ的流場參數。由該圖可知，從馬赫數是很難分辨出差別的。但是靜溫的等值線分布說明了高室壓下性能略高的原因，高室壓下具有相同值的等值線略偏向下游，這說明在噴管相同位置，高室壓所對應的靜溫略高，即混合燃氣的能量更高。提高室壓可以略提高燃氣溫度的原因在于碳氫燃料燃燒中的重要基元反應，即CO+O=CO，該反應是主要的放熱步驟。因為該反應是聚合反應，反應物的分子數多于產物的分子數，所以增大室壓使該反應向偏向于生成CO的方向移動，即放熱量增大。圖8(c)和圖8(d)中CO和CO質量分數的分布證明了這一點，即提高室壓使CO消耗得更快一些而CO生成得更多一些。

圖8 不同室壓條件下噴管Ⅱ內參數對比Fig.8 Comparison of internal parameters of nozzle Ⅱ under different chamber pressures

2.4 噴管型面對性能的影響

噴管Ⅱ和噴管Ⅲ所對應的發(fā)動機均采用了液氧/煤油推進劑，混合比也都是2.62，但是室壓不同、特征尺寸也不同，噴管內型面曲線均采用Rao方法設計。為了比較兩個噴管性能的差異，將兩個噴管都用喉部尺寸進行無量綱化處理，使喉部之前的無量綱型面重合。圖9顯示噴管Ⅱ、Ⅲ的無量綱型面曲線非常接近，在面積比較小時基本重合。

圖9給出了兩個噴管內的靜溫等值線，與圖8(b)基本一致，這說明噴管Ⅲ在提高室壓后燃氣具有較高的能量，對應更高的性能。

圖9 噴管內溫度場(噴管尺寸無量綱處理)Fig.9 Temperature field in nozzle (dimensionless treatment for nozzle size)

圖10給出了兩個噴管在面積比分別為4、9、16、25、32各截面的比沖。為了討論室壓的影響，圖中還給出了噴管Ⅱ在室壓為22.0 MPa條件下的比沖。對于噴管Ⅱ，當室壓提高到22.0 MPa后，各截面的比沖都提高了1.0 s左右，但是仍然比噴管Ⅲ對應(面積比相同)截面的比沖低。

從圖10看到，兩個噴管的型面在面積比小于4之前基本完全相同，但是計算得到比沖相差0.9 s；在面積比較大的位置，噴管Ⅲ的比沖約高1.0 s左右。產生這種差異的原因在于，針對噴管Ⅱ計算室壓22.0 MPa條件時只是增大了入口流量，其他入口參數(組分和總溫)與17.7 MPa工況相同，而噴管Ⅲ的入口則用的是22.0 MPa的化學平衡計算結果，這相當于噴管Ⅲ的燃燒效率略高?？紤]到這個因素，在相同入口條件下，兩個噴管在相同面積比下的比沖基本相同。這也說明，噴管Ⅲ的型面實際是偏瘦的，即在相同性能要求下采用噴管Ⅱ的型面可以略短；換個角度說，這是因為噴管設計采用了Rao方法而沒有考慮黏性效應，進行附面層修正后應可以略提高性能。

圖10 噴管各截面的比沖(虛線上方的值對應17.7 MPa，下方的值對應22.0 MPa)Fig.10 Specific impulse of each section of nozzle (the values above the dashed line corresponds to 17.7 MPa and the values below corresponds to 22.0 MPa)

3 結論

本文對液體火箭發(fā)動機噴管仿真的計算模型進行了詳細的分析，并通過3個噴管、4種工況的仿真詳細討論了氣體模型對仿真結果的影響，以及室壓、噴管擴張段型面對性能的影響。根據分析，本文得到如下主要結論：

1)采用量熱完全氣體假設所計算的比沖高于用熱完全氣體假設的計算值；這兩種氣體模型給出的計算結果偏差較大，可能高于、也可能低于二維化學動力學計算結果，沒有規(guī)律性；二維化學動力學計算結果比較準確，接近試車測試值。

2)提高室壓不僅能提高燃燒效率，也能促進聚合反應、減小流動過程中的化學動力學損失，使噴管性能提高。

3)Rao方法設計的噴管型面偏“瘦”，進行附面層修正可略提高性能，或在同樣性能要求下略減小噴管長度。