液體火箭發動機多級跨聲速渦輪流面迭代計算

雷宗琪 梁國柱

(北京航空航天大學 宇航學院，北京100191)

渦輪是液體火箭發動機渦輪泵的核心部件，對渦輪內部三維流動過程的深刻認識對于提高液體火箭發動機渦輪泵的性能具有重要的意義.目前國內外對于葉輪機械內部三維流動的計算主要分為直接求解三維控制方程組的CFD(Computational Fluid Dynamics)方法和相對流面法兩大類基本方法.其中CFD方法可以直接求解三維Euler方程[1－2]和Navier-Stokes方程[3－4]，具有較高的計算精度，能夠比較真實地反映流動情況，但對計算資源和時間的消耗也過于巨大.而相對流面法的思路是基于吳仲華教授提出的相對流面理論[5]，這一理論把三維流場分解為S1流面和S2流面兩組二維流面，通過流面迭代得到準三維或全三維的解.流面內二維流場的求解常采用流線曲率法.與CFD方法相比，相對流面法和流線曲率法對計算資源和時間的消耗小，可以同時應用于流場計算和葉型設計，而且配合適當的損失模型也能夠得到高精度結果[6]，從而在工程計算中具有很高的實用價值.國內外有大量采用多個S1流面和一個S2流面迭代求解葉輪機械內部準三維流場的研究，但多個S1流面和多個S2流面迭代的全三維方法卻由于復雜的計算流程而沒有得到推廣[7－8]，僅有的一些研究報道也僅局限于亞聲速葉輪機械流場的求解[9].目前關于渦輪的公開研究主要集中于民用和航空領域，針對液體火箭發動機渦輪泵中渦輪部分的研究較為少見.相對于民用渦輪機和航空發動機內的燃氣渦輪而言，液體火箭發動機渦輪泵使用的渦輪具有尺寸小、負荷高、轉速快等特點.馬冬英等[10]曾采用基于正交線的流線曲率法對某型液體火箭發動機跨聲速多級渦輪流場進行了計算并得到了不同流面內的二維無激波結果，但未能將研究拓展至三維流場，也沒有考慮激波對流場的影響.

本文以某型液體火箭發動機二級跨聲速渦輪為例，使用Fortran語言編制計算程序，采用基于任定準正交線的流線曲率法求解二維流面內的流場，并通過流面迭代的方法進一步得到了三維流場.在計算過程中采用混合平面法處理動、靜葉間的非定常干涉，通過Miller正激波模型計算激波損失，用總壓恢復系數計算流動過程中的黏性損失，并對結果以及計算方法的精確性進行了分析.

1 理論模型

相對流面的基本思想是把葉輪機械內部的三維流場分解在兩類相對流面上，分別稱為S1流面和S2流面[5].在軸流葉輪機械中S1流面和S2流面的示意圖和相對關系如圖1所示.

圖1 軸流式葉輪中的流面

1.1 基本假設

本文采用的理論模型基于以下假設:

1)渦輪內的燃氣是理想氣體;

2)渦輪內的流動是相對定常的;

3)渦輪內的流動是絕熱的;

4)忽略渦輪內燃氣的體積力.

需要注意的是，本文的研究方法為全三維方法，認為S1流面是可以翹曲的.

1.2 流線曲率法

流線曲率法在求解時主要涉及以下方程:

1)速度梯度方程

經過一系列變換，連續方程、動量方程和能量方程被整合為可以直接進行數值計算的速度梯度方程:

式中，q表示連接橫向計算站的曲線的弧長;m表示子午流線，即流線在子午面上投影的弧長;β表示相對速度與子午面的夾角;φ表示子午流線與z軸的夾角;Wm表示相對速度在子午流線方向的分量;表示相對滯止轉焓;s表示熵.相關參數的關系如圖2所示.

圖2 流線相關參數示意圖

2)流量方程

S1流面和S2流面內的流量方程具有相似的形式，均為沿計算站進行的流量積分.

S1流面流量方程:

S2流面流量方程:

其中，q表示計算站方向的微元長度;κ，λ表示流線與S1流面和S2流面內計算站的夾角;τ，σ表示S1流面和S2流面的厚度.

2 求解過程

2.1 邊界條件

渦輪葉片流道的邊界由若干曲面組成，如圖3所示.其中下表面A1為輪轂面，上表面A2為輪蓋面，左面A3為入口截面，右面A4為出口截面.流道側表面的中間部分A5，A6為葉片表面，前后兩段A7，A8為周期性邊界.

圖3 單級渦輪的三維流道

1)固體壁面A1，A2，A5，A6

A1，A2，A5，A6處的邊界均為固體壁面，應當滿足不可滲透條件，即相對速度W沿壁面法向n的分量為0:

2)進出口截面A3，A4

由于動靜葉間存在非定常干涉，上一級葉片的總溫和總壓分布不能直接傳遞給下一級葉片，本文采用了混合平面法進行處理.

根據混合平面法的思路，上一級葉片的出口截面A4和下一級葉片的入口截面A3相互重疊，該界面即為混合平面，如圖4所示.對于S1流面的流動計算，由圖4a可見各葉柵排的流道寬度不同，各條流線不能直接連接，所以計算過程中對周向的總溫和總壓進行平均處理.對于S2流面的流動計算，由圖4b可見在混合平面前后沿徑向各條流線可以連接，所以混合平面后一區域S2流面各條流線的入口總溫和總壓分別取前一區域對應各條流線的出口總溫和總壓.

3)周期性邊界A7，A8

圖4 混合平面示意圖

A7和A8是周期性邊界，應當滿足周期性條件.由于A7和A8的初始形狀和位置都是未知的，所以采用迭代修正的方法求解.

2.2 損失模型

渦輪內部的真實流動過程非常復雜，存在各種損失，例如葉型損失、二次流損失、激波損失等，一般采用經驗公式和模型對葉輪機械內部的損失進行計算.

就本文而言，由于缺乏相關數據，在計算黏性損失時無法直接套用已公開的損失模型.目前采用的方法是將這些損失等效為總壓損失，以在每一級計算站上施加總壓恢復系數的方式進行計算.

在入口條件為超聲速的渦輪級中，激波的影響是不可忽略的.本文在計算激波損失時采用了Miller激波模型，該模型在各類跨聲速葉輪機械的流場計算中被廣泛使用，其主要特點是將葉輪流道內的激波系簡化為入口處的一道正激波，并以波前平均馬赫數M'為依據計算通過激波的總壓損失[11]，如圖5所示.

圖5 Miller激波模型

圖中為葉輪流道入口處的來流馬赫數，為激波前吸力面處的馬赫數，波前平均馬赫數的計算方法為

對于激波與附面層的相互作用，本文的處理方法是在激波后的流場計算中加入流量系數.

2.3 跨聲速問題的解決方法

流線曲率法采用流量作為計算站內的流速分布是否正確的判據，在求解跨聲速流場時會遇到雙值問題，即對于同樣的流量m1，同時存在亞聲速w1和超聲速w'1兩種速度與之對應，見圖6.

圖6 雙值問題

針對雙值問題的解決方法是在計算過程中逐個計算站檢查計算流量以確定臨界截面的位置，之后將流場分成亞聲速區和超聲速區分別求解.

2.4 流面迭代過程

流面間的迭代是逐級進行的，求出一級葉片流道的全三維流場之后才開始進行下一級的計算.流面迭代是否收斂的判斷方法是相同空間位置上不同流面的求出的流動參數是否一致.本文中采用的判據是絕對馬赫數的相對誤差ε，其定義為

式中，MS1和MS2分別表示同一計算節點在S1流面和S2流面內的馬赫數.

3 計算示例和結果分析

3.1 搭建模型

以某型液氫液氧火箭發動機速度分級的兩級沖擊式氧渦輪為算例，如圖7所示.其中第2排和第4排是動葉，其余是靜葉.

圖7 渦輪葉型示意圖

在計算過程中認為氣體常數和絕熱指數為常量，計算工況如表1所示.

表1 渦輪工作參數

根據渦輪的流道形狀建立了計算模型，每級渦輪葉片采用了7個S1流面和7個S2流面，如圖8所示.

圖8 算例渦輪流面

3.2 主要參數計算結果

表2列出了各級葉片流道在流面迭代中達到的最小誤差和相應的迭代次數.計算結果和實驗數據見表3，其中實驗數據由工程單位提供，出口壓力和溫度的實驗值是在表1工況下的測量結果，渦輪效率的實驗值則是多種工況下測量結果的范圍.

表2 流面迭代的收斂程度

表3 計算結果和實驗數據

3.3 馬赫數分布

圖9展示了各級流道中的馬赫數分布.從中可以看出氣流在第1級靜葉中從亞聲速變為超聲速，在第1級動葉入口經過一道激波降為亞聲速，之后在流道內重新加速為至馬赫數較低的超聲速.在第2級靜葉入口存在一道弱激波，氣流經過弱激波后降為亞聲速，并在之后的流動過程中始終保持亞聲速的狀態.

圖9 各級流道中的馬赫數分布

3.4 誤差分析

從表3中可以看出計算得到的出口平均壓強比實驗數據高9.6%，出口總溫比實驗數據低2.2%，效率處于實驗數據的范圍之內.產生誤差的原因主要有以下幾點:

首先，流面迭代的最終計算結果是兩組流面上結果的平均值，所以流面迭代的收斂情況會對計算精度產生較大的影響.

其次，在計算過程中采用的控制方程并未包含黏性項，黏性損失是通過總壓恢復系數加入計算過程的.而且在計算激波后的流場時為了考慮附面層的影響引入了流量系數.這在一定程度上影響到了流動參數和渦輪效率的計算精度.

最后，二次流等損失在計算時也沒有被考慮在內，這也影響到了渦輪效率的計算.

4 結論

本文采用流面迭代的方法實現了液體火箭發動機渦輪泵多級渦輪內跨聲速流場的全三維計算，得到如下結論:

1)計算結果較好地描述了多級跨聲速渦輪內部的流動情況，主要性能參數的計算值與實驗值基本吻合，可為進一步改進渦輪設計、提高渦輪性能提供理論依據;

2)采用混合平面法處理動靜葉間的非定常干涉，能夠解決多級葉輪機械葉柵進出口邊界的參數傳遞問題;通過引入簡化的激波模型和由黏性引起的總壓恢復系數，使相對流面法可以用來近似計算渦輪內部的實際流動.

采用無奇點的速度梯度方程可以正確求解跨聲速流場;通過流量以及出口參數的對比可以較為準確地判斷多級跨聲速渦輪內不同區域的流動狀態.

References)

[1]Denton J D.An improved time-marching method for turbomachinery flow calculation[J].Journal of Engineering for Power，1983，105:514－521

[2]Denton J D.Calculation of three dimensional flow through any type of turbomachine[R].AGARD，Lecture Series 140，1985

[3]Hah C A.A navier-stokes analysis of 3d turbulent flows inside turbine blade rows[J].Journal of Engeering for Power，1984(106):421－429

[4]Dawes W N.Development of a three-dimensional navier-stokes solver for application to all types of turbomachinery[R].ASME.Paper No.88-GT-70，1988

[5]Wu Chunghua.A general theory of three-dimensional flow in subsonic and supersonic turbomachines of axial-，radial-，and mixedflow types[R].NACA，TN-2604，1952

[6]Boyer K M，O'Brien W F.Application of an improved streamline curvature approach to a modern，two-stage transonic fan:comparison with data and cfd[R].ASME Turbo Expo 2002.GT2002-30383，2002

[7]華鑫，喬渭陽，盧蕊.基于流線曲率法的航空軸流渦輪損失模型研究[J].機械設計與制造，2005(12):12－14

Hua Xin，Qiao Weiyang，Lu Rui，et al.A study on different loss models in axial-flow turbines based on the methed of streamlinecurvature[J].Machinery Design＆Manufacture，2005(12):12－14(in Chinese)

[8]Denton J D，Dawes W N.Computational fluid dynamics for turbomachinery design[J].Journal of Mechanical Engineering Science，1998，213(107):107－124

[9]汪慶桓，吳文權，朱根興，等.軸流式壓氣機兩類流面迭代三元完全解[J].工程熱物理學報，1983，4(4):336－346

Wang Qinghuan，Wu Wenquan，Zhu Genxing，et al.Solution of three-dimensional compressible flow in an axial-flow compressor stator by the use of s1 and s2 relative surfaces[J].Journal of Engineering Thermophysics，1983，4(4):336－346(in Chinese)

[10]馬冬英，梁國柱.基于S1，S2流面理論的液體火箭發動機渦輪內部流場計算[J].推進技術，2008，29(4):401－406

Ma Dongying，Liang Guozhu.Calculation of the internal now for turbine of liquid rocket engine based on the theory of S1，S2 stream surface[J].JournalofPropulsion Technology，2008，29(4):401－406(in Chinese)

[11]Bloch G S.Flow losses in supersonic compressor cascades[D].Blacksburg，Virgina:Virginia Polytechnic Institute and State U-niversity，1996