液体火箭发动机多级跨声速涡轮流面迭代计算

雷宗琪 梁国柱

(北京航空航天大学 宇航学院，北京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