涡轮叶片气热耦合壁温及压力分布计算与分析

Conjugate heat transfer simulation and analysis of wall temperature and pressure distribution on a turbine blade

摘要

针对气冷涡轮叶片建立了耦合分析数值模拟计算平台，采用气热耦合的方法，对径向气冷MARK Ⅱ 型叶片进行三维气热耦合数值模拟。分析了该内冷涡轮叶片的多场耦合特性，并将计算结果与试验值进行对比。对比结果表明：湍流模型的选择影响叶片表面的温度分布，但对叶片表面压力分布影响较小；在选择的6种湍流模型中，SST k - \omega湍流模型对该流场状态模拟的效果较好，叶片表面温度和压力与试验结果最接近，满足工程上的计算要求。考虑涡轮进口总温径向不均匀时，会在叶片尾缘形成局部高温区，增加了叶身的温度梯度，可适当改进冷却方式，以提高叶片强度。

关键词

涡轮叶片 / 气热耦合 / 流场结构 / 数值仿真 / 湍流模型

Key words

turbine blade / conjugate heat transfer / flow field structure / numerical simulation / turbulence model

涡轮叶片在高温、高压的流场中工作，涡轮进口的高温燃气在通道中膨胀，同时受到冷却气体的冷却作用，流场中的流动复杂，存在湍流及激波作用，涡轮部件的温度场和流场耦合在一起，单独研究某一方面不能准确地获得涡轮表面的温度分布。预测发动机热端部件的寿命及强度对叶片表面温度分布的精确计算十分必要。

Kusterer等^［1］用共轭计算程序对气冷涡轮叶片在实际热流条件下进行数值方法的研究，考虑流场传热和金属热传导的共轭模拟计算，其结果和测量结果一致，得到的壁面温度与试验结果吻合。Montomoli等^［2］基于共轭模拟技术，研究了某气膜冷却涡轮叶片的温度梯度分布，得到的流场分布与试验数据一致。郭兆元等^［3-4］采用三维数值仿真试验台，应用有限差分法搭建了气热耦合数值仿真计算平台，对气冷涡轮叶片进行仿真计算，与试验值进行比较并提出精确耦合计算的改善方法。Silieti等^［5］采用3种不同的湍流模型，研究绝热和共轭传热条件下的三维燃气轮机壁面气膜冷却效率，讨论了流场和温度场，对不同情况的预测能力进行了比较，并与试验值对比。朱卫兵等^［6］针对过高温度的高压涡轮导向叶片采用的冲击冷却进行气热三维耦合计算分析，探索并改进不同的冲击冷却的方式，以达到更好的冷却效果。尹钊等^［7］对带有径向圆孔内冷却通道的超跨音直列叶栅采用流热耦合数值模拟分析，与试验值对比并验证了流热耦合计算的准确性，讨论了耦合分析的必要性；同时采用上述方法对某涡轮导叶环形叶栅开展研究，对内、外部流场、换热及温度场相互作用展开讨论，说明流热耦合计算时，合理的湍流模型是模拟温度场分布的前提。Dyson等^［8］采用流热耦合的方法研究了气膜冷却涡轮叶片的冷却效率，选择基于RANS的SST k - \omega稳态湍流模型，并与试验值进行了比较，讨论了误差产生的原因。卢少鹏等^［9］、罗磊等^［10-12］针对冷却结构及叶型通过气热耦合的分析方式对某涡轮的动叶、导叶进行了优化。于飞龙等^［13］用气热耦合方法，以某燃气轮机涡轮动叶为研究对象开展数值研究，获得温度分布和流场分布，结果表明，该方法准确地预测了涡轮叶片的烧蚀情况。Hwang等^［14］对高压涡轮叶片的流动及传热特性进行研究，将稳态、非稳态的气热耦合数值分析结果与C3X叶型的试验结果进行了对比验证。Scholl等^［15］对某个带有冷却通道的涡轮叶片进行湍流和传热数值计算，采用的是气热耦合方法和大涡模拟求解器，数值模拟与试验结果一致。结果表明，热边界条件对固体区域传热有重要影响。李心语等^［16］采用数值模拟仿真分析多场耦合的特性，采用气热耦合热耦合的方法将具有不同冷却结构的涡轮导叶作为算例，进一步开展冷却结构设计。

本文应用Fluent气热耦合数值分析对MARK Ⅱ型叶片在文献［17］中的5411工况进行气热耦合计算，与试验结果进行对比验证，并考虑涡轮进口总温径向不均匀对叶片表面温度、压力分布的影响。

1 理论方法

1.1 涡轮叶片气热耦合控制方程

流体连续性方程为

\frac{\partial \rho }{\partial t} + \nabla \left(\rho \stackrel{⃑}{\mathit{U}}\right) = \mathrm{0} （1）

流体在每个方向的速度分量满足动量方程为

\begin{array}{l}\frac{\partial \left(\rho \stackrel{⃑}{\mathit{U}}\right)}{\partial t}+\nabla \left(\rho \stackrel{⃑}{\mathit{U}}\times \stackrel{⃑}{\mathit{U}}\right)=\\ \nabla \left(-\rho \delta +\mu \left(\nabla \stackrel{⃑}{\mathit{U}}+{\left(\nabla \stackrel{⃑}{\mathit{U}}\right)}^{\mathrm{T}}\right)\right)+S_M=\mathrm{0}\end{array} （2）

式中： S_M为动量源项。

不考虑流体的黏性耗散，满足能量守恒方程为

\begin{array}{l}\frac{\partial \left(\rho h^{\mathrm{*}}\right)}{\partial t}-\frac{\partial \rho }{\partial t}+\nabla \left(\rho \stackrel{⃑}{\mathit{U}}{h}^{\mathrm{*}}\right)=\\ \nabla \left(\lambda \nabla T_f\right)+\nabla \left(\mu \nabla \stackrel{⃑}{\mathit{U}}+\nabla {\stackrel{⃑}{\mathit{U}}}^{\mathrm{T}}-\frac{\mathrm{2}}{\mathrm{3}}\nabla \stackrel{⃑}{\mathit{U}}\delta \stackrel{⃑}{\mathit{U}}\right)+S_E\end{array} （3）

h^{\mathrm{*}} = h + \frac{\mathrm{1}}{\mathrm{2}} u^{\mathrm{2}} , h = h \left(p,T_f\right) （4）

式中： S_E为能量源项； h为比焓。

理想气体状态方程为

\rho = \frac{P}{R{T}_f} （5）

式中： R为通用气体常数。

耦合传热能量方程为

\begin{array}{l}\frac{\partial \left(\rho c_p{T}_f\right)}{\partial t}=\frac{\partial }{\partial x}\left(\lambda \frac{\partial T_f}{\partial x}\right)+\frac{\partial }{\partial y}\left(\lambda \frac{\partial T_f}{\partial y}\right)+\\ \frac{\partial }{\partial z}\left(\lambda \frac{\partial T_f}{\partial z}\right)+S_E\end{array} （6）

式中： \rho为固体密度； \lambda为热导率； S_E为能量源项，涡轮叶片内部传热问题无内热源， S_E为0。

1.2 湍流模型

湍流是一种复杂的三维、非稳态的流动过程。燃气涡轮内流体流动非常复杂，但以湍流为主。湍流模型在工程计算中主要有标准 k - \omega、SST k - \omega、BSL k - \omega、标准 k - \epsilon、Realizable k - \epsilon、RNG k - \epsilon 6种模型。

SST k - \omega模型是剪切应力输运 k - \omega模型的简称，与标准 k - \omega模型相比，SST k - \omega模型中融合了横向耗散导数项。SST k - \omega应用范围更广的原因是该湍流黏度的公式中考虑了湍流剪切应力的输运过程，更适合模拟激波、逆压梯度下的流动分离等现象。描述方程^［18］为

\begin{array}{l}\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho u_{i}k\right)}{\partial x_i}={\tilde{P}}_k-{\beta }^{\mathrm{*}}\rho k\omega +\\ \frac{\partial }{\partial x_i}\left[\left(\mu +{\sigma }_k{\mu }_t\right)\frac{\partial k}{\partial x_i}\right]\end{array} （7）

\begin{array}{l}\frac{\partial \left(\rho \omega \right)}{\partial t}+\frac{\partial \left(\rho u_i\omega \right)}{\partial x_i}=\alpha \frac{\mathrm{1}}{v_t}{\tilde{P}}_k-{\beta }^{\mathrm{*}}\rho {\omega }^{\mathrm{2}}+\\ \frac{\partial }{\partial x_i}\left[\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\frac{\partial \omega }{\partial x_i}\right]+\mathrm{2}\left(\mathrm{1}-F_{\mathrm{1}}\right)\rho {\sigma }_{\omega \mathrm{2}}\frac{\mathrm{1}}{\omega }\frac{\partial k}{\partial x_i}\frac{\partial \omega }{\partial x_i}\end{array} （8）

2 计算模型

以MARK Ⅱ叶片模型为算例，该叶片是气冷涡轮叶片，共有10个径向圆管冷却孔，径向由底部向上流过冷却空气，单个涡轮叶片及流道模型如图1所示，叶片及流道的几何参数如表1所示。

试验的涡轮叶片材料为310不锈钢，其密度为7 900 kg/m³，质量热容为585.1 J/（kg·K），导热系数是其与温度的函数，计算公式为

\lambda ( T_f ) = \mathrm{6.811} + \mathrm{0.020} \mathrm{176} T_f （9）

涡轮叶片选择文献［17］中5411工况，气热耦合计算的边界条件如下：工质选择理想气体，主流入口总压为337 097 Pa，总温为788 K，湍流度为6.5％，主流出口静压为175 713 Pa，主流流道两侧设置为平移周期性边界条件，其他壁面为无滑移、绝热边界。

10个冷却通道的边界条件如表2所示。对叶片进行气热耦合模拟时，湍流模型采用工程上常用的模型分别对流场进行模拟计算，将温度、压力分布结果与文献［17］中的试验值进行对比。

网格划分如图2所示。流体域和固体域采用Fluent Meshing划分的非结构化多面体网格，可在保证网格质量的基础上减少计算工作量。流体域、固体域共229万个网格，边界层划分为20层，近壁面第一层网格 y +小于1，保证边界层内流动的预测准确性。

3 壁温及压力分布计算结果

3.1 涡轮叶片气热耦合计算结果

图3为MARK Ⅱ型叶片内、外流场中径截面位置的马赫数分布。由图3可见，外部的流动为跨声速流动，激波发生在叶片的吸力面。由于涡轮叶片的外流道逐渐收缩，主流流经前缘驻点沿着吸力面流动时不断膨胀加速，在距前缘驻点约0.43倍相对弦长的位置产生了强激波，马赫数达到1.6。在此逆压梯度下，激波后的边界层产生分离现象，此后的流动为亚音速流动，但在下游又很快附着。在收缩的流道影响下，流体继续在顺压梯度下膨胀，加速流动。在尾缘附近，距前缘0.93倍相对弦长位置出现一道弱激波，马赫数约为1.1，之后向外扩展。而在压力面，马赫数缓慢增长。

图4为中径截面温度分布云图。低温区出现在2号和3号冷却孔附近，原因是这一区域冷却孔的位置比较密集，冷气流量较大，冷却效果较其他区域突出。在叶片的尾缘区域出现了高温区域，主要由于此区域的冷却孔数量较少，冷气量减少及尾缘处复杂的湍流流动，涡量大，发生能量耗散，造成边界层内涡动，产生局部高温区域。在叶片的吸力面，外流场温度在激波前出现一个局部最小值，在激波后温度又快速升高。

3.2 壁温及压力计算结果与试验值的对比

叶片外表面中径截面位置不同湍流模型的温度、压力分布与试验值对比^［17］如图5所示。横坐标表示叶片相对弦长，纵坐标分别表示无量纲温度和无量纲压力；横坐标负半轴表示压力面，正半轴表示吸力面。与试验值相比，温度分布在前缘驻点至吸力面0.43倍相对弦长的区域出现区别，原因是在吸力面激波位置后，复杂的湍流流动条件下，当前选择的湍流模型对边界层内流动情况的预测能力不同，出现较大波动。在较大的逆压梯度下对流动的过度预测，产生较大的计算误差，每种湍流模型在此位置的温度均高于试验值。

温度最高的点为叶片前缘滞止点和尾缘。在叶片前缘，高温流体在此处发生较强的气流冲击，在前缘处滞止，形成了局部的高温区域；在靠近尾缘流动方向上的后半段，有明显的波动，因为其受到了冷却通道低温气体的影响。随着流动沿流道发展，叶片表面边界层逐渐增厚，受到前端冷却通道的降温作用，温度在吸力面距前缘0.43倍相对弦长位置及压力面0.18倍相对弦长的位置达到最小值。发生在吸力面的激波导致其后的主流温度迅速上升，从而导致叶片温度上升，造成了图5a中吸力面上的较高温度。

压力面的气流从前缘驻点沿表面缓慢加速，压力逐渐降低，趋势趋于平缓；气流沿压力面速度快速升高，造成压力大幅度降低。每种湍流模型在中径位置的压力分布基本一致，表明湍流模型的选取对流场气动压力影响较小。

SST k - \omega湍流模型在中径截面的压力分布与试验值^［17］基本相符，温度分布与其他湍流模型相比，在前缘点与试验值最接近，且在吸力面的最大差值为30.90K，相对误差为6.09%，满足工程上的计算要求，后续计算选择与试验值最符合的SST k - \omega模型。

4 进口温度分布不均匀的影响

由于航空发动机和燃气轮机的燃烧室中气流的不均匀性和掺混作用，燃烧室出口即涡轮进口的温度分布不均匀，对涡轮叶片的温度分布及流动产生影响。为了模拟涡轮进口的不均匀温度分布，进口温度的径向分布依据文献［19］设置，如图6所示。进口总温最高值为872 K，位于0.6倍叶高的位置，最低值为621 K，位于叶根。

图7为涡轮进口总温分布均匀与不均匀时，叶片表面温度、压力分布及叶片截面温度。进口总温的径向不均匀分布导致叶身中部的温度显著升高，在叶片尾缘中部形成一个局部高温区域，温度增加约14K；叶根区域的温度大面积降低，降低约6%。

进口总温均匀分布时，叶片表面整体温度随着叶高的增长而升高，叶片内部出现局部低温区，尾缘温度最高。叶片内部的低温区沿叶高逐渐减小，表示冷却气体在给叶片降温的同时，其温度升高对叶片顶端的降温作用逐渐减弱。进口总温径向不均匀时，由于叶身的局部高温区域，中径截面的温度明显升高，特别是前缘和尾缘部位达到局部温度的最高点，叶片

内部温度梯度增大。而叶片表面压力分布没有发生明显变化，说明进口总温的径向不均匀分布对叶片的气动性能没有影响。

5 结论

针对气冷涡轮叶片，选用6种湍流模型对MARK II叶片某工况进行了气热耦合数值模拟，得到以下结论：

（1） MARK II叶片表面流场为跨音速流动，首先在吸力面的前部产生一个强激波，其后又在吸力面接近尾缘处形成一道较弱的激波，由此产生的流动分离使流场的流动状态更加复杂，流场预测的难度变大。

（2） 每种湍流模型对叶片外表面的压力预测与试验值吻合较好；叶片表面的温度分布受湍流模型选择的影响较大，一定程度上影响后续的强度校核及疲劳预估，选择合适的湍流模型进行流场分析十分必要；SST k - \omega湍流模型与叶片表面温度试验值吻合较好，证明涡轮叶片采用气热耦合方法计算的合理性。

（3） 进口总温径向不均匀分布时，在叶身形成局部高温区域，加剧了气冷涡轮叶片内部的温度梯度，必要时可改进冷却方式，改善叶片受力状态，提高使用寿命。

参考文献

[1] Kusterer K， Bohn D， Sugimoto T，et al.Conjugate calculations for a film-cooled blade under different operating conditions［C］//Proceedings of ASME Turbo Expo 2004：Power for Land，Sea，and Air.Vienna：ASME，2004：675-684.

[2] Montomoli F， Adami P， Della Gatta S，et al.Conjugate heat transfer modelling in film cooled blades［C］//Proceedings of ASME Turbo Expo 2004：Power for Land，Sea，and Air.Vienna：ASME，2004：147-156.

[3] 郭兆元，王强，冯国泰.涡轮气热弹耦合计算模型与算例［J］.航空学报，2009，30（2）：213-219.

[4] 郭兆元，王强，黄家骅，涡轮叶片气热耦合数值模拟方法［J］.推进技术，2009，30（2）：159-164.

[5] Silieti M， Kassab A J， Divo E.Film cooling effectiveness：comparison of adiabatic and conjugate heat transfer CFD models［J］.International Journal of Thermal Sciences，2009，48（12）：2237-2248.

[6] 朱卫兵，孙润鹏，徐凌志，高压涡轮导向器叶片冲击冷却数值研究［J］.推进技术，2012，33（5）：710-718.

[7] 尹钊，方祥军，刘思永，超跨音内冷涡轮流热耦合数值模拟研究［J］.中国电机工程学报，2013，33（14）：114-120，1.

[8] Dyson T E， Bogard D G， Bradshaw S D.Evaluation of CFD simulations of film cooling performance on a turbine vane including conjugate heat transfer effects［J］.International Journal of Heat and Fluid Flow，2014，50：279-286.

[9] 卢少鹏，迟重然，罗磊，气热耦合条件下涡轮静叶三维优化［J］.推进技术，2014，35（3）：356-364.

[10] 罗磊，迟重然，卢少鹏，燃气涡轮静叶考虑叶型及冷却结构的气热耦合优化［J］.工程热物理学报，2014，35（6）：1079-1082.

[11] 罗磊，卢少鹏，迟重然，气热耦合条件下涡轮动叶叶型与冷却结构优化［J］.推进技术，2014，35（5）：603-609.

[12] 罗磊，陈朔，温风波，燃气涡轮无气膜冷却动叶参数化设计系统应用［J］.推进技术，2015，36（6）：864-875.

[13] 于飞龙，上官博，李园园，重型燃气轮机透平第一级动叶复合冷却数值研究［J］.中国电机工程学报，2016，36（1）：179-186.

[14] Hwang S， Son C， Seo D，et al.Comparative study on steady and unsteady conjugate heat transfer analysis of a high pressure turbine blade［J］.App-lied Thermal Engineering，2016，99：765-775.

[15] Scholl S， Verstraete T， Duchaine F，et al.Conjugate heat transfer of a rib-roughened internal turbine blade cooling channel using large eddy simulation［J］.International Journal of Heat and Fluid Flow，2016，61（7）：650-664.

[16] 李心语，刘火星.气冷涡轮导叶流热耦合计算及机理［J］.北京航空航天大学学报，2021，47（11）：2378-2386.

[17] Hylton L D， Mihelc M S， Turner E R，et al.Analytical and experimental evaluation of the heat transfer distribution over the surfaces of turbine vanes O：NASA- CR-168015［R］.Washington D.C.：NASA，1983.

[18] Menter F R.Review of the shear-stress transport turbulence model experience from an industrial perspective［J］.International Journal of Computational Fluid Dynamics，2009，23（4）：305-316.

[19] Wang T Y， Xuan Y M， Han X S.Analysis of transient radiation-convection coupled effects on the HP turbine blade heat load with hot streak inlet［J］.Applied Thermal Engineering，2018，138：705-721.

基金资助

两机专项子项(JH2021001G)