基于非结构网格的涡轮气热耦合数值方法研究

汤超，葛宁，SHENG Chun-hua

基于非结构网格的涡轮气热耦合数值方法研究

汤超1，葛宁1，SHENG Chun-hua2

(1.南京航空航天大学能源与动力学院，江苏南京210016；2.College of Engineering，University of Toledo，OH 43606，USA)

利用课题组自主开发的三维非结构隐式N-S计算软件CU_Turbo，采用气热耦合计算方法，对MarkⅡ内冷径向涡轮导向叶片、带气膜冷却涡轮导叶MT1的流场和温度场进行了数值模拟。计算过程中，隐式时间推进中Jacobi⁃ans矩阵采用对Roe通量的一种近似方法求解。结果表明，计算值与试验值吻合良好，验证了气热耦合计算方法的实用性和有效性，为涡轮工程设计提供了一种新的计算分析方法；涡轮叶片通道内附面层的不同流动状态及气膜冷却，对当地换热都有很大的影响。

涡轮叶片；气热耦合；气膜冷却；隐式时间推进；非结构网格

1 引言

涡轮是发动机中热负荷和动力负荷最大的部件，采用有效的冷却措施是发动机安全可靠工作的有力保证，也是降低材料成本的有效措施。近年来，涡轮前燃气温度的逐年提高与涡轮的有效冷却，特别是叶片内部空气冷却技术的迅速提高分不开。有研究表明，涡轮叶片设计中，如果预测的叶片温度超过实际工作温度28℃，叶片寿命将减半[1]。因此准确估算叶片的传热系数和温度，有助于预防热腐蚀和设计高效冷却系统。

涡轮叶片的数值模拟中，常规方法是，用给定的固体表面温度作为边界求出叶片表面的换热系数分布，再以此作为固体域热传导计算的边界条件。该方法中，壁面温度的给定更多的是依靠经验，从而导致计算的不确定性增加，加大了预测热负荷的难度，尤其是对于复杂的几何体。而实际中，涡轮部件的温度场与流场相耦合，单方面研究无法为发动机设计提供准确依据，故提出了气热耦合方法[1]。该方法不需要用壁面的温度分布和换热系数分布作为边界条件，而是在整场计算中获得，很好地解决了常规方法中的不确定性，且可适用于任意复杂的几何体和工作状态。同时，由于复杂的冷却系统导致涡轮叶片造型很不规律，相对于结构网格来说，非结构网格求解器对其有更好的适应性。

本文使用课题组自主开发的基于非结构网格的三维计算软件CU_Turbo，通过对叶型MarkⅡ和MT1进行数值模拟，验证了气热耦合计算方法的实用性和有效性，为涡轮工程设计提供了一种新的计算分析方法。

2 数值计算方法

2.1控制方程和数值方法

笛卡儿坐标系下，守恒变量形式的三维可压缩流雷诺平均N-S方程的形式为：

式中：Q=[ρ ρu ρv ρw ρet]为守恒变量，F为通量矢量。

CU_Turbo软件流体域的计算采用格点格式的有限体积法，无粘通量的离散采用Roe格式，Ven⁃katakrishnan限制器[2]抑制数值振荡，时间推进采用LU-SGS[3]隐式方法并加入当地时间步长和隐式残值光顺来加速收敛，高阶精度由高斯-格林和最小二乘梯度重构方法实现，边界条件依据特征波原理[4]，加入Spalart改进的S-A湍流模型把握近壁区域流动。固体域采用热传导控制方程，计算方法与流体域保持一致。

2.2非结构网格的气热耦合方法

在流固交界面上，必须满足交界面温度的连续性和交界面热通量的连续性两个物理条件[5]：

对于如图1所示的非结构网格，由式(2)可得：式中：kf、ks分别为流体域和固体域的热传导系数，Tf、Ts分别为流体域和固体域边界相邻边对应点的静温，Δnf、Δns分别为流体域和固体域相邻边对应点到边界的距离，m、n分别为交界面上点在流体域和固体域的相邻边数，i、j分别指交接面上点相应在流体域和固体域的相邻边，wi、wj分别为相应边的权重且表达式为

由式(3)可求得壁面温度Tw，再以此作为流体域和固体域的定温边界条件分别进行计算。

2.3隐式推进数值方法

时间推进采用LU-SGS隐式算法，把式(1)有限体积法离散并采用向后差分可得：

式中：Vi为控制体i的体积，Rn+1i为时间第n+1层上残值且代表控制体i的第ij个控制面上的无粘通量，Fv代表该控制面上的粘性通量，n̑代表控制面外法矢，Ωij代表控制面面积。

式(5)为一非线性系统，直接求解非常麻烦，通过泰勒展开将其线化可得：

对于上式中每个时间步的系数矩阵，可分解为对角阵D、上三角阵U和下三角阵L：

其中D为纯对角线矩阵，L、U分别对应编号小和编号大的邻居单元：

将式(8)带入式(7)可进一步改写为：

向前推进

向后推进

求得的修正项Δq更新本时间步：

2.4通量Jacobians矩阵的计算方法

对流通量Jacobians矩阵的求法是隐式求解的关键，区别于Steger-Warming矢通量分裂的解析方法[6]求解Jacobians矩阵，本文采用对Roe通量的一种近似方法求解。该方法简单可靠且易于编程实现。本课题组史万里等对以上两种方法做过比较研究，评估了Jacobians矩阵计算的可靠性和实用性[7]。控制面ij的Roe通值表达式为：

式中：带上标⌒的项采用修正的Roe平均变量计算，认为该项为当地值并假设为常数。求微分后可得任一控制面对于界面左右两侧的对流通量Jacobians矩阵的表达式：

粘性通量Jacobians矩阵的求解采用薄层假设，求解过程参见文献[8]。

3 算例验证及分析

CU_Turbo软件已经过湍流平板气动验证、跨声流凸包气动验证、层流平板气动热力耦合验证[9]，且计算值与试验值吻合良好。本文对径向涡轮导叶MarkII和带有6排圆柱形气膜孔的涡轮导叶MT1进行详细的计算和分析，来验证该软件的实用性和有效性。

3.1MarkII

MarkII叶片是有10个冷却孔的涡轮直叶片，叶高76.20 mm，弦长136.22 mm。图2给出了叶片及计算域几何尺寸，其中左、右表面分别为进口和出口，上、下两边界面为周期性边界面，中间为带冷却孔叶片，顶端和底端面为对称边界。Hylton已在接近工程运行条件下做过大量试验[10]，本文模拟的是Hyl⁃ton试验中的5411工况，进出口参数为：总压332 000 Pa，总温811 K，背压171 282 Pa，进口马赫数0.19，出口马赫数1.04。

图2 MarkⅡ叶片几何尺寸及工况Fig.2 Geometric configuration and boundary conditions of MarkII vane

采用cooper的方法划分计算域网格，在贴近叶片表面的壁面处使用边界层网格，使得第一层y+≤1以保证精确模拟边界层内流动。采用商用软件Gambit生成网格，通过CU_Turbo接口软件GMT[9]实现数据结构转换。

图3为马赫数等值线分布云图。从图中看，由于高进、出口压比和吸力面收缩通道的影响，在吸力面前缘附近流动有一个强加速过程，在x cx=0.45 (x cx为横坐标与轴向弦长之比)位置出现一道强激波，这也将影响到该区域的热负荷；激波后边界层出现分离，之后流动减速到亚声速范围，在更远的下游流体再次加速，在尾缘附近出现一道较弱激波。而在压力面一侧，马赫数增加较缓慢。

图3马赫数等值线分布Fig.3 Mach number contours

图4 所示为温度等值线分布云图。可见，温度最低点发生在第2和第3冷却孔之间，由于这一区域冷却孔位置比较集中，冷却效果要比其它位置的好；吸力面一侧，气流温度在激波前降低到局部最小值，在激波后又迅速抬升到当地最大值，温度最高的位置在叶片尾缘。原因为尾缘部分冷却孔数量没有前缘部分多，冷气流量小，加之叶片尾缘的直角形状，气体流动复杂，各种涡系相互掺混，所以温度在这一区域达到最高。

图4温度等值线分布Fig.4 Temperature contours

图5 给出了叶片表面中径处的性能曲线。从图5(a)中可知，在吸力面，从滞止点开始压力迅速降低，在x cx=0.45处压力减小到最小值，激波与压力分布的位置相对应。

图5(b)中，滞止点温度达到局部最大值。沿压力面一侧，由于第2和第3冷却孔间冷却效果较好，以及叶片轮廓由凸转凹、凸面流动换热量不断减小使得温度降低[11]，压力面温度从滞止点开始逐渐降到局部最小值。沿吸力面一侧，由于与压力面相同的原因，温度先降低到局部最小值，激波后边界层发生分离诱导转捩，导致换热量急剧上升，温度也随之迅速增加，接近尾缘时由于叶片变薄和冷却不足致使温度进一步上升。温度计算值最大误差小于3%，壁面平均误差小于1%。

图5(c)中，换热系数h定义为h=qw(T0-Tw)，qw为叶片表面热流密度，T0为气体来流温度。可见，温度影响换热，换热系数与表面温度分布呈现出很强的相关性，计算值与试验值趋势基本吻合。驻点后压力面和吸力面的热传导系数急剧下降，吸力面激波处诱导分离转捩，导致热传导系数开始迅速上升。由于计算的温度边界层比实际的厚，使得压力面换热系数的计算值比试验值偏低，靠近尾缘区域的误差，还可能是由于当地附面层流动受相邻叶片压力面强激波影响所致。

图5 MarkII叶片表面中径处性能曲线Fig.5 Performance profiles of MarkII vane

3.2MT1

MT1叶片的冷却系统由2个圆柱形冷却腔和6排圆柱形气膜孔组成，冷却气自冷却腔上端进入、经气膜孔流出，与主流掺混在叶片表面形成气膜，冷却腔下端封闭，两腔在吸力面均只有1排气膜孔，在压力面有叉状排列的2排气膜孔[12]。

图6所示为叶片计算域及网格，左侧为进口边界，右侧为出口边界，计算域上、下端设为对称边界，前、后为周期边界，所有壁面均设为定温边界。各边界具体参数为：主流来流总压P*=4.60×105Pa，主流来流总温T*=444 K，前冷却腔进口总压P*1=6.28×105Pa，后冷却腔进口总压P*2=4.88×105Pa，冷却腔进口总温T=286.5 K，Tw=288 K。

图6 MT1叶片计算域及网格Fig.6 Computation domain and grids of MT1 vane

在贴近叶片表面的壁面处使用边界层网格，使得第一层y+在1的量级以保证精确模拟边界层内流动。气膜孔进、出口处也采用了与壁面边界层一样的加密方式，以保证气膜孔与主流通道交界面上网格过渡平滑。流体域网格170万。

图7示出了叶片表面等熵马赫数和叶片表面中径处努赛尔数的分布情况。从图7(a)中可看出，计算值与试验值吻合较好，吸力面后端(x=30 mm)有一道弱激波，气膜孔对吸力面等熵马赫数分布的影响较大，波动明显强于压力面。

图7(b)中，努赛尔数的定义式为Nu=hd λ，d为特征长度，λ为流体域热传导系数。对比有无气膜孔时的努赛尔数分布，可清楚看出气膜孔对叶片表面换热的影响，在每组气膜孔后换热强度均有较大幅度的降低，有效减小了叶片表面的热负荷。对比带气膜的计算值与试验值，压力面由于计算得出的温度附面层偏厚导致努赛尔数的计算值稍偏低，吸力面峰值处则略高于试验值，总体趋势吻合良好。同时，可观察到在吸力面约x cx=0.30处换热量急剧增加，这是由于此处发生边界层转捩所致。

图7 MT1叶片表面等熵马赫数及努赛尔数分布Fig.7 Isentropic Mach number and Nu of MT1 vane

图8 示出了冷却气经气膜孔与主流气掺混的过程，冷却气进入主流通道后在叶片表面形成气膜，有效地把高温气体与叶片隔离开，达到了保护叶片的目的。

图8 冷却气与主流气掺混迹线Fig.8 Mixing stream trace of cooling air and gas

4 结论

CU_Turbo是一个基于非结构网格的三维计算软件，由于其对复杂内冷结构涡轮叶片的适应能力及气热耦合方法的实用性，使得软件非常适用于涡轮叶片的流动分析和换热机理研究，从而优化涡轮叶片的冷却流路设计，有较大的工程应用价值。本文选择对带内冷的涡轮导叶MarKⅡ和带6排气膜孔的涡轮导叶MT1，进行流场和温度场耦合计算，得出如下结论：

(1)CU_Turbo软件能准确预测带内冷的涡轮叶片MarkⅡ和MT1的流场及温度场。MarkⅡ在超声速工况下计算得到的表面温度与试验值符合良好，壁面温度计算值最大误差小于3%；MT1的换热分布同样与试验值吻合较好。

(2)该软件采用隐式时间推进，隐式通量Jaco⁃bians矩阵采用对Roe通量的一种近似方法求解，进一步在涡轮叶片计算中验证了将Jacobians矩阵与Roe格式结合这种方法的实用性和可靠性。

(3)附面层内不同流动状态(层流、湍流和转捩)对当地的传热过程有很大影响，气膜冷却可显著减少气膜孔下游叶片与高温燃气的热交换，有效保护叶片。

[1]冯国泰，黄家骅，李海滨，等.涡轮发动机三维多场耦合数值仿真试验台的数学模型研究[C]//.中国工程热物理学会热机气动热力学学术会议.2000.

[2]Venkatakrishman V.Convergence to Steady State Solu⁃tions of the Euler Equations on Unstructured Grids with Limiters[J].Journal of Computational Physics，1995，118 (1)：120—130.

[3]Luo Hong，Baum J D，Löhner R.A Fast,Matrix-Free Im⁃plicit Method for Computing Low Mach Number Flows on Unstructured Grids[R].AIAA 99-3315，1999.

[4]Sheng C，Wang X.Characteristic Variable Boundary Con⁃ditions for Arbitrary Mach Number Algorithm in Rotating Frame[R].AIAA 2003-3976，2003.

[5]Sheng Chun-hua，Xue Qing-luan.Aerothermal Analysis of Turbine Blades Using an Unstructured Flow Solv⁃er-U2NCLE[R].AIAA 2005-4683，2005.

[6]Wang X.A Preconditioning Algorithm for Turbomachinery Viscous Flow Simulation[D].MS：Mississippi State Univer⁃sity，2005.

[7]史万里，葛宁，Sheng Chunhua.粘性计算中预处理的隐式算法[J].空气动力学学报，2009，27(5)：565—571.

[8]王靖宇.涡轮内非结构网格气动热力耦合数值计算方法研究[D].南京：南京航空航天大学，2012.

[9]赵滋阳.非结构网格下气动热力耦合数值方法研究[D].南京：南京航空航天大学，2010.

[10]Hylton L D，Mihelc M S，Turner E R，et al.Analytical and Experimental Evaluation of the Heat Transfer Distribution overtheSurfacesofTurbineVanes[R].NASA CR-168015，1983.

[11]Bohn D，Heuer T.Conjugate Flow and Heat Transfer Cal⁃culation of a High Pressure Turbine Nozzle Guide Vane [R].AIAA 2001-3304，2001.

[12]Adami P，Martelli F，Chana K S，et al.Numerical Predic⁃tionsofFilmCooledNGVBlades[R].ASME GT2003-38861，2003.

Study of Turbine Conjugate Heat Transfer Simulation Based on Unstructured Grids

TANG Chao1，GE Ning1，SHENG Chun-hua2
(1.College of Energy and Power Engineering，NUAA，Nanjing 210016，China；2.College of Engineering，University of Toledo，OH 43606，USA)

This paper describes an attempt to predict the aerodynamic pressure loading and heat transfer of two cases with an unstructured Navier-Stokes flow solver,which include a NASA turbine vanes-MarkII and NGV MT1.A three-dimensional Reynolds-averaged Navier-Stokes unstructured flow solver,referred to as CU_Turbo,has been modified to solve the governing equations in the flow field and a thermal diffusion equa⁃tion inside the solid region in a coupled manner.During implicit time iterative,an approximate method of Roe flux is used for finding the Jacobians matrix.The results showed that the numerical simulation results were agreement well with the experimental values.The results verified the feasibility and effectiveness of the code.So the code provides a new method of calculation and analysis for the engineering design of tur⁃bine;flow condition in boundary layer and film cooling have great impacts on the local heat transfer in tur⁃bine vane channels.

turbine vane；conjugate heat transfer；film cooling；implicit time iterative；unstructured grids

V231.1

A

1672-2620(2013)02-0033-05

2012-08-06；

2012-12-22

航空推进技术验证计划项目

汤超(1988-)，男，安徽巢湖人，硕士研究生，主要从事叶轮机械气动热力学研究。