1 引言

为了大幅度提高涡轮负荷，涡轮气动设计广泛采用跨声速设计；同时燃气涡轮发动机进口温度越来越高，受材料许用温度的限制，有效的涡轮叶片冷却设计越来越重要^{[1, 2]}。而涡轮叶栅通道内为复杂的三维流动，会出现转捩、二次流，激波及附面层相互干扰等现象^{[3, 4]}。刘聪等^[5]对跨声速涡轮导叶的数值研究表明，激波诱导边界层分离和转捩，换热系数会有大幅度变化。因此，准确预测转捩位置和叶片表面对流换热系数分布对设计高效的冷却方案具有重要意义。

叶栅通道内的流动转捩会严重影响叶片表面对流换热系数，工程常用的转捩预测方法主要有：经验关联方法^[6]、以e^N为代表的线性稳定性方法、基于间歇因子的预测方法^[7]，以及Menter和Langtry等发展的$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型^[8]。前三种工程方法难以与CFD方法相匹配。$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型完全基于当地变量信息，适于在CFD中应用。该模型考虑了包括自然转捩、旁路转捩、分离泡转捩和再层流化（Relaminarization）在内的各类算例^{[9, 10]}，能较为准确地预测转捩位置。国内外学者如张晓东等^[11]，Tomac等^[12]、郑赟等^[13]对转捩模型中的经验关联式进行了研究，得到的关联函数同Langtry类似，均依赖于当地自由流湍流度。王刚等^[14]采用转捩模型对较低来流湍流度（0.03%~0.3%）的平板、翼型流动进行了计算，预测的转捩点位置较为准确，但并未涉及高湍流度情况。陈立立等^[15]运用Michel转捩判据，分别采用Langtry和Tomac提出的3个经验关联式，对低雷诺数翼型计算；该方法降低了对湍流度等参数的敏感性，能得到较接近试验值的结果，但对中高雷诺数以及分离预测并未验证。

董平等^{[16, 17]}、李涛等^[18]采用$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型，数值模拟了跨声速涡轮叶栅Mark Ⅱ内来流湍流度6.5%的流动换热情况，结果表明该转捩模型能较准确地获得转捩点位置，但计算的叶栅压力面换热系数大大低于实验值，而吸力面转捩点附近换热系数又远高于实验值。曾军和卿雄杰^[19]研究发现进口湍流度对外换热系数影响较大。孟德虹^[20]和张玉伦^[21]在TRIP软件平台上研究了$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型在二维低速问题中的应用，文中指出：进口湍流度沿程会迅速衰减，无法确保物面附近有正确的自由流湍流度分布；他们采用了Spalart等^[22]建议的环境源项法来控制湍流度衰减，在远场区域保持环境湍流度，在物面附近让其自由衰减，但文中未给出远场区域的具体确定方法。为了准确地模拟入口至物面前缘的湍流度分布，Langtry等^[10]提出用反复试验入口粘性比的方法来匹配不同位置处的实验湍流度水平。显然这种方法需要耗费大量的计算资源和时间。

基于上述问题，本文主要参考零压力梯度时自由流湍流度衰减率实验结果，给出了一种快速估算来流粘性比的方法，避免多次试算的繁复。通过在入口处设置合适的粘性比来控制湍流度衰减，确保叶栅前缘具有准确的湍流度和期望的衰减率。并对MARK Ⅱ与VKI涡轮叶栅通道在亚、跨声速下的流动换热情况进行了算例验证，结果表明，按本文方法设置入口粘性比能很好地改善$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型对涡轮叶栅外换热的预测精度，尤其是较高入口湍流度时叶栅压力面的换热精度。

2 $ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型

$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型通过引入当地转捩动量厚度雷诺数（$ \tilde R{e_{\theta {{\rm t}}}}$）的输运方程将自由流的影响扩散到边界层内，并引入间歇因子（γ）输运方程来控制边界层内的湍动能，实现转捩的模拟。

间歇因子输运方程^[10]为

$\frac{{\partial \left( {\rho \gamma } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {U_{{\rm j}}}\gamma } \right)}}{{\partial {x_{{\rm j}}}}} = {P_{{\rm \mathsf{ γ} }}} - {E_{{\rm \mathsf{ γ} }}} + \frac{\partial }{{\partial {x_{{\rm j}}}}}\left[ {\left( {\mu + \frac{{{\mu _{{\rm t}}}}}{{{\sigma _{{\rm f}}}}}} \right)\frac{{\partial \gamma }}{{\partial {x_{{\rm j}}}}}} \right]$

(1)

其中，转捩源项

${P_{{\rm \mathsf{ γ} }}} = {F_{{{\rm length}}}}{c_{{{\rm a}}1}}\rho S{\left[ {\gamma {F_{{{\rm onset}}}}} \right]^{0.5}}\left( {1 - {c_{{{\rm e}}1}}\gamma } \right)$

(2)

${E_{{\rm \mathsf{ γ} }}} = {c_{{{\rm a}}2}}\rho \mathit{\Omega }\gamma {F_{{{\rm turb}}}}\left( {{c_{{{\rm e}}2}}\gamma - 1} \right)$

(3)

${F_{{{\rm onset}}}} = {{\rm max}}\left( {{F_{{{\rm onset}}2}} - {F_{{{\rm onset}}3}}, 0} \right)$

(4)

${F_{{{\rm turb}}}} = {{{\rm e}}^{ - {{\left( {{R_{{\rm T}}}/4} \right)}^4}}}, {R_{{\rm T}}} = \rho k/\mu \omega $

(5)

${F_{{{\rm onset}}1}} = R{e_{{\rm v}}}/\left( {2.193R{e_{\theta {{\rm c}}}}} \right), {{\rm }}R{e_{{\rm v}}} = \rho {y^2}S/\mu $

(6)

${F_{{{\rm onset}}2}} = {{\rm min}}\left( {{{\rm max}}\left( {{F_{{{\rm onset}}1}}, F_{{{\rm onset}}1}^4} \right), 2.0} \right)$

(7)

${F_{{{\rm onset}}3}} = {{\rm max}}\left( {1 - {{\left( {{R_{{\rm T}}}/2.5} \right)}^3}, 0} \right)$

(8)

式中F_length为控制转捩区长度的关键参数；Re_θc为间歇函数在边界层内第一次开始增大的临界雷诺数。F_length和Re_θc都是当地转捩动量厚度雷诺数$ \tilde R{e_{\theta {{\rm t}}}}$的经验函数，具体表达式和相关参数详见文献[10]。

当地转捩动量厚度雷诺数$ \tilde R{e_{\theta {{\rm t}}}}$的输运方程^[10]

$ \begin{array}{*{20}{l}} {\frac{{\partial \left( {\rho \tilde R{e_{\theta {{\rm t}}}}} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {U_{{\rm j}}}\tilde R{e_{\theta {{\rm t}}}}} \right)}}{{\partial {x_{{\rm j}}}}} = }\\ {{P_{\theta {{\rm t}}}} + \frac{\partial }{{\partial {x_j}}}\left[ {{\sigma _{\theta {{\rm t}}}}\left( {\mu + \frac{{{\mu _{{\rm t}}}}}{{{\sigma _{{\rm f}}}}}} \right)\frac{{\partial \tilde R{e_{\theta {{\rm t}}}}}}{{\partial {x_{{\rm j}}}}}} \right]} \end{array} $

(9)

式中，源项

$ {P_{\theta {{\rm t}}}} = {c_{\theta {{\rm t}}}}\frac{\rho }{{{t_{{\rm s}}}}}\left( {R{e_{\theta {{\rm t}}}} - \tilde R{e_{\theta {{\rm t}}}}} \right)\left( {1.0 - {F_{\theta {{\rm t}}}}} \right) $

(10)

这里t_s是为了协调量纲而存在的时间参数。混合函数F_θt用来控制$ \tilde R{e_{\theta {{\rm t}}}}$源项的开关，在0~1之间变化：边界层以外F_θt趋于0，方程受源项控制，只有当$ \tilde R{e_{\theta {{\rm t}}}}$的值接近当地Re_θt的值才能维持平衡；边界层内F_θt趋于1，源项关闭，通过当地转捩动量厚度雷诺数$ \tilde R{e_{\theta {{\rm t}}}}$输运方程将自由流的影响扩散到边界层内，充分考虑自由流对边界层的影响。

转捩动量厚度雷诺数Re_θt直接影响着转捩过程，它是自由流湍流强度Tu，流向压力梯度参数λ_θ的关联函数；Re_θt采用Langtry^[10]给出的形式

$ \begin{array}{*{20}{l}} {R{e_{\theta {{\rm t}}}} = \left[ {1173.51 - 589.428Tu + \frac{{0.2196}}{{T{u^2}}}} \right]F\left( {{\lambda _\theta }} \right), }\\ {Tu \le 1.3} \end{array} $

(11)

$ R{e_{\theta {{\rm t}}}} = 331.50{\left[ {Tu - 0.5658} \right]^{ - 0.671}}F\left( {{\lambda _\theta }} \right), Tu > 1.3 $

(12)

式中F(λ_θ)代表压力梯度的影响。

转捩模型耦合于SST湍流模型中，通过间歇因子γ作用于湍动能方程的生成项和耗散项来模拟转捩过程，即

$ \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {U_{{\rm j}}}k} \right)}}{{\partial {x_{{\rm j}}}}} = \overline {{P_{{\rm \mathsf{ γ} }}}} - \overline {{E_{{\rm \mathsf{ γ} }}}} + \frac{\partial }{{\partial {x_{{\rm j}}}}}\left[ {\left( {\mu + \frac{{{\mu _{{\rm t}}}}}{{{\sigma _{{\rm k}}}}}} \right)\frac{{\partial k}}{{\partial {x_{{\rm j}}}}}} \right] $

(13)

$ \overline {{P_{{\rm \mathsf{ γ} }}}} = {\gamma _{{{\rm eff}}}}{P_{{\rm k}}}\;\;\overline {{E_{{\rm \mathsf{ γ} }}}} = {{\rm min}}\left( {{{\rm max}}\left( {{\gamma _{{{\rm eff}}}}, 0.1} \right), 1.0} \right){D_{{\rm k}}} $

式中P_k和D_k是原来k-ω SST湍流模型的生成项和耗散项；对ω输运方程不做任何修改。

3 自由流湍流度控制

由上节可知，转捩模型中转捩动量厚度雷诺数Re_θt强烈依赖于当地自由流湍流度Tu。研究发现^[22]，在数值模拟中，入口湍流度沿程会迅速衰减，且衰减速度取决于进口粘性比；进口粘性比越大，湍流衰减越慢；但是过大的入口粘性比则会影响壁面摩阻。因此，在数值模拟中需指定进口自由流湍流度并给定合理的入口粘性比，以保证叶栅前缘附近具有合理的湍流度水平和期望的衰减率。然而在实验中一般只测量来流的湍流强度，很少测量湍流长度尺度^[23]，所以难以确定入口的湍流粘性比。

采用SST模型描述湍流，在匀速运动参考系上观察，湍流度强度随时间的衰减率为

$\frac{{{{\rm d}}k}}{{{{\rm d}}t}} = - {\beta ^{{\rm *}}}\omega k, \frac{{{{\rm d}}\omega }}{{{{\rm d}}t}} = - \beta {\omega ^2}$

初始时入口处t=0, k=k_inlet, ω=ω_inlet

求解上式可以得到自由流中的湍动能衰减公式

$k = {k_{{{\rm inlet}}}}{\left( {1 + {\omega _{{{\rm inlet}}}}\beta t} \right)^{ - \frac{{{\beta ^{{\rm *}}}}}{\beta }}}$

(14)

式中β=0.09, β^*=0.0828。

在静止参考系中，湍流度变化是定常的，湍流强度从入口处的值沿流向不断衰减。将时间尺度用空间尺度替换，即t=X/V，X为流向距离，V为流动速度；涡粘性μ_t=ρk/ω，引入粘性比R_T=μ_t/μ=ρk/ωμ；湍流度$Tu = \sqrt {2k/3} /V $，可得到从入口处自由流湍流度随流向距离和粘性比的衰减关系式^[23]

$Tu = {\left( {Tu_{{{\rm inlet}}}^2{{\left[ {1 + \frac{{3\rho Vx\beta Tu_{{{\rm inlet}}}^2}}{{2\mu \left( {{\mu _{{\rm t}}}/\mu } \right)}}} \right]}^{ - \frac{{{\beta ^{{\rm *}}}}}{\beta }}}} \right)^{0.5}}$

(15)

则入口粘性比为

${R_{{\rm T}}}_{{{\rm inlet}}} = 1.5\beta \frac{{\rho {V_{{{\rm inlet}}}}x}}{\mu }\frac{{Tu_{{{\rm inlet}}}^2}}{{{{\left( {Tu/T{u_{{{\rm inlet}}}}} \right)}^{ - 2\beta /{\beta ^{{\rm *}}}}} - 1}}$

(16)

为了高效地获得合理的入口粘性比，提高计算效率，本文参考零压力梯度时自由流湍流强度的衰减规律^[23]，给出了入口粘性比的计算方法，具体方法如下。

图 1为不同入口湍流度下相对湍流度Tu/Tu_inlet沿流向距离分布的实验点^[22]及其拟合曲线。由图可知，随着流向距离增大，湍流衰减趋于平缓; 在入口湍流度较低（Tu_inlet=0.9%）时，与中高入口湍流度之间的衰减程度相差较大，且入口湍流度越大，相对湍流度沿流向衰减越快；在入口湍流度较高时（Tu_inlet > 3.3%），不同入口湍流度下的相对湍流度沿流向衰减程度基本相同。在涡轮叶栅数值模拟中，可由确定的入口与涡轮叶栅前缘之间距离X₁和入口湍流度Tu_inlet，依据图 1来估算出对应X₁处的前缘湍流度Tu_LE及其与入口湍流度Tu_inlet之比，代入到公式（16）中可得到所需要的入口粘性比R_{T inlet}。

4 计算结果及分析

本文采用商业软件ANSYS CFX，针对两个典型的涡轮叶栅MARK Ⅱ和VKI模型，选用$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型进行了涡轮叶栅外换热数值模拟。采用High Resolution离散格式。

计算中，在入口处分别给定$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型中缺省粘性比和按本文方法得到的粘性比，通过数值结果的对比来验证本文方法的有效性。通过与不同网格数的计算结果对比，已验证本文的数值解为网格无关解。

4.1 Mark Ⅱ涡轮叶栅

表 1和图 2分别为Mark Ⅱ叶型的主要几何参数和计算网格，详细叶型数据和冷却通道数据可参见文献[24]。叶栅流道内采用结构化网格，网格数量约为332万，在叶栅壁面附近加密网格，壁面网格最小尺度为5.0 × 10^-7m，保证y⁺ < 1，网格法向增长率为1.1；固体叶栅部分采用非结构化网格，网格数量为26万。计算工况如表 2所示，各冷气通道平均温度和换热系数参见文献[24]。

在数值模拟中采用$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型，入口粘性比分别采用缺省值（R_{T 1}=65，即Tu_inlet×1000）及按本文方法计算的入口粘性比（R_{T inlet}=415）。作为对比，同时还采用层流模型及SST完全湍流模型进行了计算。

图 3给出了跨声速5411工况下涡轮叶栅通道马赫数分布云图。此时在流道中出现了局部超声速区域，在吸力面X/L=0.43位置产生一道较强的激波，同时尾缘处形成一对斜激波，激波与激波相互干涉以及激波与边界层的相互干扰，使流动变得更加复杂。

图 4给出了Mark Ⅱ涡轮叶栅5411工况下涡轮叶栅表面无量纲压力分布，其中p表示叶栅表面压力，p_t表示进口总压，X/L表示轴向位置与轴向弦长之比。可以看出，除层流模型外，数值模拟的压力分布与试验值吻合得较好。由于在吸力面0.43倍轴向弦长处产生了激波，激波后压力骤增，产生较强的逆压梯度，引起层流边界层分离，从图 5的吸力面局部流线图中可以看到分离泡的存在；层流分离触发流动转捩，使得吸力面上对流换热系数突增。

图 6为Mark Ⅱ叶片5411工况下叶片表面对流换热系数分布。叶片表面对流换热系数由下式计算

$h = {q_{{\rm w}}}/\left( {T_1^* - {T_{{\rm w}}}} \right)$

(17)

式中q_w为壁面热流密度；T₁^*为叶栅进口总温；T_w为叶栅壁面温度。

由图 6可知，对于跨声速工况5411，在压力面，边界层对自由流湍流强度敏感，当转捩模型采用缺省进口粘性比时（R_{T 1}=65），此时入口湍流度会迅速衰减，叶片前缘湍流度大大低于入口湍流度，预测的换热系数远低于实验值，与层流结果相近；当设定合理进口粘性比（R_{T inlet}=415）后，使入口到前缘位置具有合适的湍流度衰减率，此时叶栅压力面为完全湍流流动，转捩模型的换热结果与实验值吻合得较好，相对于缺省粘性比情况，对于压力面大部分区域相对误差减小约30%；SST湍流模型的结果稍高于实验值。在吸力面，转捩模型能较好地获得转捩点位置；在转捩点之前，转捩模型计算结果与实验值相吻合，而完全湍流模型SST结果偏大；在转捩点之后，层流模型的换热结果偏低，转捩模型结果与SST模型结果基本相同，均出现了一较大的峰值，过高估计了对流换热系数。这可能是由于此处产生了较强的激波，实际流动压缩性效应显著，而转捩模型本身未考虑流体压缩效应，高估了边界层内的湍动能，从而引起较大误差。对于这种情况进行流动可压缩性修正，可显著减小吸力面传热系数数值计算误差^[25]。总体来看，在$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型中设置合理入口粘性比，显著改善了叶栅表面对流换热系数的计算精度，特别是压力面换热精度。

4.2 VKI高压涡轮叶栅

VKI叶型的主要几何参数和计算网格如表 3和图 7所示，详细叶型数据可参见文献[26]。采用结构化网格，流道内网格数量约为156万。在叶栅壁面加密网格，壁面第一层厚度为7.0 × 10^-7m，y⁺ < 1.3，网格法向增长率为1.1。计算工况如表 4所示，MUR 241和MUR 247分别为出口跨声速、亚声速。

图 8和图 9分别为工况MUR 241下，涡轮叶栅表面压力分布与外换热系数分布图。由图 8可知，除层流模型之外，其它模型得到的压力分布基本相同；在吸力面中部出现弱激波，此时逆压梯度较小，未发生边界层分离。

由图 9可知，工况MUR241下各模型换热结果与实验的对比情况与Mark Ⅱ的跨声速工况5411结果基本类似。在压力面，采用缺省入口粘性比（R_{T 1}=60）时转捩模型预测的换热系数远低于实验值，在压力面大部分区域，相对误差高达50%，与层流模型的结果相当；而当合理设置入口粘性比（R_{T inlet}=636.23）后，整个压力面的对流换热系数与实验值符合得较好，相对误差最大约7%；在此区域，SST完全湍流模型也得到较好的换热结果。而在吸力面，转捩模型能较准确地获得转捩点位置；从滞止点到转捩点之前，层流模型和转捩模型的结果均与试验值符合得较好，给定合理入口粘性比时转捩模型预测结果与实验最为接近。在转捩点之后，层流模型预测的对流换热系数过低；而转捩模型与SST模型一样，均高估了对流换热系数。

在工况MUR247下，VKI叶片表面压力分布如图 10所示。同MUR 241工况类似，几种湍流模型模拟的压力分布基本相同。转捩模型模拟结果显示在叶栅吸力面中部喉道附近达到声速，但并未产生激波。图 11是转捩位置附近局部流线，在靠近吸力面尾缘位置产生了流动分离，引起转捩。

图 12为工况MUR 247下的对流换热系数分布。由图可知，对于SST湍流模型，除前缘区外均高估了对流换热水平，层流模型则在吸力面转捩点之后低估了对流换热水平。对于转捩模型，在压力面，由于入口自由流湍流度较低，流动保持层流，采用缺省入口粘性比（R_{T 1} =10）和设置合理入口粘性比（R_{T inlet}=18）时的换热结果均与试验值相符。而在吸力面，采用缺省入口粘性比使得入口湍流度衰减过大，转捩之后换热结果远低于实验值；当设置合理的入口粘性比时，换热结果与试验值基本相符。

本文还采用给定入口粘性比的转捩模型，对比了VKI叶栅在出口亚声速、超声速时不同入口湍流度下的叶栅表面换热系数分布，如图 13和图 14所示。其中，试验值参见文献[26]，图中FSTI表示入口来流湍流度。由图可知，叶栅压力面的换热模拟结果与实验值均符合得较好，在来流湍流度较低时，压力面保持层流流动；来流湍流度较高时，压力面为湍流流动。

在吸力面，对于出口超声速工况，吸力面中弦附近产生弱激波，激波与边界层相互作用，使流动变得更加复杂。在转捩点之前流动处于层流，换热模拟结果与实验值较为符合，转捩点之后的换热结果与实验值偏差较大。对于高湍流度情况，激波之后触发转捩，转捩模型高估了对流换热系数。对于低湍流度工况，在吸力面中弦附近的微弱激波不足以触发转捩，流动保持层流，直至靠近尾缘处产生流动转捩，换热增强^[10]，此时模拟的换热结果与试验值有较大差异。对于中等湍流度，计算值接近于高湍流度的结果，激波后发生转捩，而实验值接近于低湍流度结果，激波后未发生转捩，这是由于压缩性对转捩有一定抑制作用。表明此时流动可压缩性效应显著，需进一步考虑流动可压缩性修正^[11]。

如图 14，对于出口亚声速工况，流动在喉道截面达到声速，在此之后，由于通道扩张，流动减速，无激波产生。在吸力面转捩点之前，换热系数的模拟结果与试验值基本符合；在中低入口湍流度时，转捩点位置较试验值稍有延后，转捩之后的换热系数与试验值相差不大；而在高入口湍流度时，转捩点相对出口超声速工况有所提前（图 13），在S=26.45mm处（喉道之前）发生转捩，其转捩区的范围相对实验结果明显偏窄，换热系数相对实验值偏高。

5 结论

通过研究，得到如下结论：

（1）当采用本文方法对$ \gamma - \tilde R{e_{\theta {{\rm t}}}}$转捩模型给定合理的入口粘性比时，能够显著改善涡轮叶栅表面对流换热系数的计算精度。

（2）在本文的中等雷诺数范围，对于涡轮叶栅压力面，当较低入口湍流度时，流动保持层流；较高入口湍流度时，流动为完全湍流，采用合理入口粘性比时，压力面换热系数模拟值与试验值符合得较好，与缺省粘性比情况相比，传热系数计算相对误差降低30%以上。

（3）在本文的中等雷诺数范围，对涡轮叶栅吸力面，当出口亚声速、无激波产生时，中低湍流度时在近叶栅尾缘处发生转捩，换热强化；在高湍流度时计算模型能够较好地获得转捩位置，但转捩宽度偏窄，转捩后换热系数偏高。当出口超声速时叶栅通道内有弱激波，较高入口湍流度时，激波引发转捩，换热增强；设置合理粘性比的转捩模型能较准确地预测叶栅前缘至吸力面转捩点的换热情况；转捩点之后转捩模型高估了对流换热系数。表明流动可压缩性显著，需进一步考虑流动可压缩性修正。