湿燃气透平叶片热流固耦合换热特性的数值研究

2021-10-22符阳春李健武廖贵鄂张志伟

动力工程学报 2021年10期

符阳春，李健武，廖贵鄂，张志伟，李亮,2

(1.西安交通大学能源与动力工程学院，西安 710049；2.陕西省叶轮机械及动力装备工程实验室，西安 710049)

近年来，随着能源短缺和环境污染问题日益严重，燃气轮机开始向“高效率和环保性能兼优”的方向发展。在此背景下，研究人员提出了适用于燃气轮机的湿空气透平循环[1]，该循环通过加湿工质提高燃气轮机的输出功率。根据研究，湿化燃气轮机燃烧室出口的高温燃气中水蒸气含量最高可达40%[2]。由于水蒸气和空气的比热、导热系数等物性参数存在显著差异，湿化燃气轮机透平燃气侧的热通量可增加20%以上。因此，准确预测湿燃气传热系数的分布，对湿化燃气轮机叶片冷却结构的设计和应用有着重要的作用。

Carullo等[3]实验研究了叶片表面换热的影响因素，结果表明叶片外表面的换热与进口湍流强度、出口马赫数和雷诺数以及壁面表面粗糙度等有关。Maffulli等[4]数值模拟了不同壁温对叶片表面换热的影响，结果表明传热系数随着壁面与燃气温比的减小而增大。Schmidt等[5]数值研究了稳态和非稳态条件下，高压涡轮叶片表面耦合换热的换热特性，结果表明非稳态条件下，叶片压力面传热系数明显增加。Shiau等[6]通过TSP技术实验测量了跨声速条件下扇形气膜孔透平叶片吸力面的传热系数。Al-Zurfi等[7]通过大涡模拟方法研究了旋转对透平叶片表面传热系数的影响，结果表明旋转促使吸力面边界层提前发生转捩，且转捩区扩大。郭隽等[8]数值研究了不同湍流模型对透平叶片表面传热系数的预测精度，结果表明带转捩的SST湍流模型对透平叶片表面传热系数的计算明显优于其他模型。曾军等[9]采用气热耦合的方法，数值研究了不同冷气流量比、温比以及压比对低压透平动叶叶片换热的影响，提出了一种有效的透平叶片温度场分析方法。陈大为等[10]采用热色液晶技术测量了尾迹对透平动叶表面传热系数的影响，结果表明尾迹会导致叶片表面传热系数增大，且对吸力面的影响大于压力面。杜长河等[11]通过数值模拟研究了湿蒸汽对透平叶片耦合换热特性的影响，结果表明，相比于空气和过热蒸汽，湿蒸汽的冷却效率更高。

目前，对透平叶片换热特性的研究大多以干燃气为工质，而对湿化燃气轮机中湿燃气换热特性的研究较少。随着湿化燃气轮机循环的逐渐应用，有必要对湿燃气温度和换热特性展开研究。此外，叶片冷却的研究大多以空气或湿蒸汽为冷却介质，而对湿空气冷却特性的研究较少。湿化燃气轮机按照压缩空气加湿运行的特点，本身具备采用湿空气进行透平冷却的条件。湿空气冷却与干空气冷却相比，减少了冷气的使用量，与蒸汽冷却相比，既降低了冷却系统的复杂度，也避免了可能出现的叶片腐蚀、流道堵塞等问题。目前，对透平叶片的研究主要集中在燃气侧传热系数的分布，缺少冷气侧燃气侧综合考虑下的叶片受热负荷的研究。笔者以NASA C3X叶片的实验数据为基础[12]，采用商业软件CFX，研究了湿燃气含湿量对透平叶片表面温度和传热系数的影响，在研究范围内，给出了透平叶片燃气侧传热系数的无量纲关系式，对比分析了干空气冷却与湿空气冷却效果的差异，探究了湿空气含湿量对冷却效果的影响，为透平叶片优化和冷却结构设计提供参考。

1 计算模型和数值方法

1.1 计算模型及边界条件

图1为C3X叶片结构，其10个径向圆柱形冷却孔均匀地分布在叶片内部。表1为C3X叶片的几何参数。笔者根据C3X叶片的实际尺寸建立了数值计算三维模型，并采用ICEM生成了结构化网格，如图2所示。网格由叶片主流通道流体域、叶片固体域和冷却通道流体域三部分构成，并在流体域壁面附近进行加密，以满足k-ω、SST和SST 转捩湍流模型对于y+<1的要求。

图1 C3X叶片结构Fig.1 Structure of C3X blade

表1 C3X叶片几何参数Tab.1 Geometric parameters of C3X blade

(a)主流通道和冷却通道

湍流模型验证和网格无关性验证的计算边界条件按照实验参数[12]给定。其中，主流通道进气燃气温度为818 K，总压为413 286 Pa；出口静压为254 172 Pa；叶片两侧为周期性边界条件，上下端壁为绝热壁面。冷气通道沿轴向进气，各冷却通道进口的总温和质量流量与文献[12]保持一致。

1.2 数值方法

数值模拟采用计算流体动力学(CFD)软件CFX进行。研究时分别以干燃气和湿燃气、干空气和湿空气作为工质，CFX中设置干燃气是由二氧化碳、氮气和氧气组成的混合气体，干空气是由氩气、氮气和氧气组成的混合气体。湿燃气和湿空气则由干燃气和干空气分别与水蒸气混合而得。二氧化碳、氮气、氧气和氩气均设置为实际气体，气体状态方程选择对理想气体状态方程进行修正的The Aunigiter Redilich Kwong Model[13]，以描述上述气体在数值模拟工况条件下的状态。导热系数λ和动力黏度μ则由Sutherland公式[14]确定，其表达式如下：

(1)

(2)

式中：Tref为参考温度，K；μref为参考温度下的动力黏度，kg/(m·s)；λref为参考温度下的导热系数，W/(m·K)；S为Sutherland常数，K；n取1.5;T为计算温度。

根据文献[14]可知，计算温度离参考温度越近，Sutherland公式对动力黏度和导热系数的预测精度越高。笔者将参考温度设置为1 750 K，计算动力黏度的Sutherland常数设置为340 K，计算导热系数的Sutherland常数设置为416 K。经计算，在数值模拟工况的压力和温度范围内，干空气的物性计算值与美国国家标准与技术研究所(NIST)公布的物性数据，在密度、动力黏度、导热系数和比定压热容上的平均误差分别为2.7%、0.029%、0.021%和2.3%。

计算采用热流固耦合方法对叶片的传热过程进行分析，在流体域求解三维N-S方程，在固体域求解Laplace导热方程，并在流固交界面处保证温度和热通量的连续性。计算过程中，监测进口雷诺数和气体进出口质量流量不平衡度等参数。当方程组均方根残差小于10-4、气体进出口质量流量不平衡度小于0.1%以及其他监测变量均保持稳定时，模拟计算获得稳定的数值解。

含湿量d为1 kg干空气所带的水蒸气质量，单位为g/kg。

(3)

其中，mg为水蒸气质量，g；ma为干空气质量，kg。用dmg表示湿燃气的含湿量，用dma表示湿空气的含湿量。

为了定量比较干燃气和湿燃气的换热强度，采用对流传热系数h作为评判参数，表达式如式(4)所示，可以看出h越大，换热强度越高。

h=q/(Tg,t-Tw)

(4)

式中：q为叶片表面热流量,W/m2；Tg,t为主流燃气进口温度,K；Tw为叶片壁面温度,K。

以叶片弦长为特征尺度，定义表征对流换热强度的努塞尔数Nul为：

Nul=q·l/[(Tg,t-Tw)·λg]

(5)

式中：λg为主流工质导热系数,W/(m·K)；l为叶片弦长,m。

为了定量比较干空气和湿空气的冷却效果，采用冷却效率η作为评判参数，其表达式为：

(6)

式中：Tc,t为冷气进口总温,K。

以冷却孔直径为特征尺度，定义表征对流换热强度的努塞尔数Nud为：

Nud=qc,wdc/(Tc,t-Tc,w)λc

(7)

式中：qc,w为冷却通道壁面热流密度,W/m2；Tc,w为冷却通道壁面温度,K；λc为冷却工质导热系数,W/(m·K)；dc为冷却孔直径,m。

1.3 湍流模型验证

对比分析k-ε、k-ω、SST和SST 转捩4种湍流模型对热流固耦合计算结果的影响。图3～图5分别给出了不同湍流模型下叶片中截面的压力p/p0分布、温度T/T0分布和无量纲传热系数h/h0的分布曲线，其中下标0表示气体实验进口值。图6给出了不同湍流模型下冷却通道中径处的努塞尔数分布。图中都选用叶片的轴向弦长L作为基准，以其在X坐标轴上的无量纲参数X/L来表示叶片表面不同位置，X/L为负值表示压力面，为正值则表示吸力面，以叶片前缘滞止点为原点X/L=0，尾缘为顶点。

图3 不同湍流模型下的叶片中截面压力分布Fig.3 Pressure distribution at mid-span section under different turbulence models

图4 不同湍流模型下的叶片中截面温度分布Fig.4 Temperature distribution at mid-span section under different turbulence models

图5 不同湍流模型下的叶片中截面对流传热系数Fig.5 Heat transfer coefficients at mid-span section under different turbulence models

图6 不同湍流模型下冷却通道中径处的努塞尔数分布Fig.6 Nusselt number distribution at mid-span of cooling channels under different turbulence models

由图3可知，4种湍流模型计算的叶片中截面的压力分布基本上均与实验值吻合。由图4可知，与实验值相比，k-ε湍流模型的温度计算结果明显偏大，k-ω、SST湍流模型的计算结果稍有偏大，SST 转捩湍流模型的计算结果与实验值重合度最好。由图5可知，与实验值相比，k-ε湍流模型的传热系数计算结果明显偏大；k-ωSST湍流模型的计算结果在叶片压力面和吸力面尾缘处与实验值吻合较好，但在叶片吸力面前缘处仍有较大偏差；SST 转捩湍流模型计算的传热系数与实验值重合度最好。由图6可知，k-ε湍流模型的努塞尔数计算结果比实验值小，k-ω、SST湍流模型的计算结果比实验值大，SST 转捩湍流模型计算的冷却通道中径处的努塞尔数与实验值重合度最好。

综上所述，SST 转捩湍流模型在透平叶片热流固耦合计算中最为可靠，因此后续均选择SST 转捩湍流模型进行计算。

1.4 网格无关性验证

为了保证计算的准确性，选取网格数分别为120万、170万、230万、320万和410万进行计算，以验证网格无关性。图7给出了不同网格数下的叶片表面平均对流传热系数havg，其中下标s表示吸力面，p表示压力面。从图7可以看出，随着网格数的增多，叶片表面平均对流传热系数有所增加，当网格数从230万增加到410万时，叶片吸力面的表面平均对流传热系数仅增加0.44%。综合考虑计算结果的准确性与计算成本，最终选择网格数为230万。

图7 叶片表面平均对流传热系数的网格无关性验证Fig.7 Independent verification of heat transfer coefficients on blade surface

2 计算结果与讨论

2.1 湿燃气含湿量对燃气侧流动的影响

参照某型燃气轮机透平进气参数，取主流燃气进口温度分别为1 473 K和1 773 K，总压为1.7 MPa，出口静压为1.0 MPa，冷却空气的总温分别为602 K和725 K，冷气的质量流量与文献[12]保持一致。图8给出了不同主流进口温度下，湿燃气和干燃气叶片中截面的压力分布。从图8可以看出，气流沿着压力面不断膨胀，压力不断降低，从叶片前缘到X/L=-0.5处，压力降低比较缓慢，此后压力急剧下降，一直到尾缘达到最小值。在吸力面，压力从前缘点开始就很快下降，随后由于气流速度的变化，压力出现了波动。同样的边界条件下，湿燃气和干燃气叶片中截面处的压力分布差别不大。图9给出了主流进口温度为1 773 K时，叶片吸力面的极限流线分布。由图9可知，在叶片根部和顶部附近存在明显的通道涡，并将主流向叶片中间高度处挤压。这是由于端壁处存在流速低的边界层区域，导致该区域内流体不能抵抗压力面和吸力面之间的压力梯度，压力面的高压气流经端壁边界层流向吸力面，从而卷起了吸力面边界层。当主流从干燃气变为湿燃气后，近壁面通道涡的范围没有变大，因而流阻损失没有增大。

(a)Tg,t=1 473 K

(a)dmg=0 g/kg

图10给出了不同主流进口温度下，湿燃气和干燃气叶片中截面的马赫数Ma分布云图。从图10可以看出，湿燃气和干燃气的流动情况类似，在叶片通道内的流动主要是亚音速流动，在吸力面前缘附近存在一处激波(Ma=1),在吸力面尾缘处存在一块高速区域。C3X型叶片头部比较大且曲率半径变化较小，叶栅的能量损失变化较小。当湿燃气含湿量从0 g/kg增加到150 g/kg时，叶片尾迹处高速区域显著增加。

图10 不同主流进口温度下叶片中截面的马赫数分布云图Fig.10 Mach number distribution at mid-span section under different inlet temperatures

2.2 湿燃气含湿量对燃气侧换热的影响

图11给出了不同主流进口温度下，湿燃气和干燃气叶片中截面的传热系数分布。由图11可知，在压力面，由于压力降低缓慢，湍流边界层迅速增厚，因此传热系数在前缘附近急速下降，随后由于主流速度逐渐增大，边界层减薄，传热系数又逐渐增大；在吸力面，传热系数在前缘附近达到最大，此后由于层流边界层增厚，传热系数减小。在转捩区，传热系数有一定的增大，经过湍流转捩区到达转捩点以后，层流转化为湍流，传热系数又开始骤增，此后随着湍流边界层的增厚，传热系数逐渐减小。随着湿燃气含湿量的增加，由于水蒸气密度小于干燃气，导致湿燃气流速增大，流体与壁面间的换热增强，叶片中截面传热系数逐渐增大。

(a)Tg,t=1 473 K

图12和图13给出了主流进口温度分别为1 473 K和1 773 K，dmg=0 g/kg(即干燃气)和dmg=150 g/kg时，叶片壁面的温度分布云图。由图12和图13可知，干燃气和湿燃气下的叶片表面温度分布形式大致相同，叶片前缘由于受到主流高温燃气的冲击，因而温度较高。靠近叶根区域的位置温度较低，原因是冷气从叶根进入冷却孔并与叶片进行换热，沿着流动方向冷却效率逐渐下降，因而叶片表面温度略有上升。叶片尾缘由于空间狭小且远离冷却孔，表面温度较高。与干燃气相比，湿燃气下的叶片表面温度较高，其叶片压力面前缘和尾缘温升较明显，吸力面前缘和中部的温升较明显。与干燃气相比，含湿量为150 g/kg的湿燃气叶片壁面温度有所升高，最低温度提高了20 K，最高温度提高了10 K。

图12 Tg,t=1 473 K时叶片壁面温度分布Fig.12 Wall temperature distribution of turbine blade at inlet temperature of 1 473 K

图13 Tg,t=1 773 K时叶片壁面温度分布Fig.13 Wall temperature distribution of turbine blade at inlet temperature of 1 773 K

图14给出了不同主流进口温度下，叶片出口与进口压比τ分别为0.82、0.75、0.68和0.615时，叶片中截面平均对流传热系数的相对值havg/havg0随含湿量的变化曲线。其中，havg0表示干燃气下的平均对流传热系数，havg/havg0表征湿燃气为工质时相对于干燃气的换热增强效果。由图14可知，叶片壁面传热系数随着湿燃气含湿量的增加而增大，且主流进口温度越高，湿燃气叶片壁面换热增强越显著。湿燃气含湿量从0 g/kg增加到150 g/kg，主流进口温度为1 473 K时叶片壁面传热系数提高了10%，主流进口温度为1 773 K时叶片壁面传热系数提高了14%。值得注意的是，相同含湿量下，湿燃气的换热增强效果随着压比的提高而增强。吸力面侧的换热增强效果相对于压力面侧对压比更敏感，并且与高主流进口温度下相比，低主流进口温度下的换热增强效果对压比也更敏感。相比压力面，吸力面上的流动较为复杂，存在由层流向湍流过渡的转捩区，因而对工质的变化更为敏感。温度越高，干燃气和水蒸气密度越低，但密度降幅也越低。因此，高主流进口温度下不同含湿量的湿燃气密度差异变小，因而该条件下的换热增强效果相对于低主流进口温度下更不敏感。

图14 不同主流进口温度下叶片中截面havg/havg0随含湿量的变化Fig.14 Variation of havg/havg0 with moist gas humidity rate at mid-span section under different inlet temperatures

2.3 湿空气含湿量对冷气侧流动的影响

图15给出了冷气进口温度为602 K和725 K时，在不同湿空气含湿量下冷却通道中截面处的速度V和雷诺数Re分布。由图15可知，10个冷却孔的平均速度和雷诺数有所差异，这是由于冷却孔进口处质量流量略微有所区别，以及冷却孔径不同。当冷气进口温度一定时，随着湿空气含湿量的增加，冷却孔中截面处的速度也增加。这是由于湿空气的密度随着含湿量的增加和温度的升高而减小，导致湿空气流速增大，冷却通道换热增强。

(a)Tc,t=602 K

2.4 湿空气含湿量对冷却效果的影响

图16给出了冷气进口温度为602 K和725 K时，不同湿空气含湿量下冷却孔中径处Nud的变化。由图16可知，随着湿空气含湿量的增加，不同冷却孔的Nud也随之增大，冷气侧换热增强，这与前文分析所得结论一致。

(a)Tc,t=602 K

图17给出了冷气进口温度为602 K和725 K时，干空气和含湿量为200 g/kg湿空气的叶片冷却效率η分布云图。图18给出了2种冷气进口温度下，叶片冷却效率的相对值ηavg/ηavg0随含湿量的变化曲线。其中，ηavg为湿空气下的叶片平均冷却效率，ηavg0为干空气下的叶片平均冷却效率。由图17可知，湿空气和干空气的冷却效率分布规律相似，即靠近叶根区域的η较大，原因是冷气从叶根进入冷却孔并与叶片进行换热，随着冷气温度的升高，冷却能力下降。叶片前缘由于受到主流高温燃气的冲击，η较低；叶片尾缘由于空间狭小且远离冷却孔，η同样较低。从图17可以看出，与干空气相比，采用含湿量为200 g/kg的湿空气时η提高了约8%。当冷气与主流的温比为0.409时，冷气温度较高时对工质的变化更敏感，725 K下湿空气的冷却效率略高于602 K下湿空气的冷却效率。

压力面

图18 叶片冷却效率相对值随湿空气含湿量的变化Fig.18 Variation of relative cooling efficiency with moist air humidity rate

图19给出了冷气进口温度为602 K和725 K时不同湿空气含湿量下叶片中截面壁面温度分布曲线。由图19可知，当冷气与主流的温比为0.409时，不同进口温度下湿空气对叶片表面温度的冷却效果不同。与干空气冷却相比，在叶片压力面和吸力面的中部区域，含湿量为200 g/kg的湿空气冷却的叶片壁面温度下降了20～30 K；在靠近叶片压力面和吸力面前缘和尾缘处，湿空气冷却的叶片壁面温度下降了10～15 K。整体上，含湿量为200 g/kg的湿空气冷却的叶片壁面温度下降了约20 K。图20给出了冷气进口温度为725 K时叶片中截面的温度分布云图。从图20可以看出，在湿空气含湿量为0 g/kg时，随着湿燃气含湿量的增加，叶片中截面温度亦有所升高，表现最为明显的是在叶片吸力面前缘和前4个冷却孔处。在叶片中间靠近尾缘处，高温区域范围亦有所扩大。在湿燃气含湿量为150 g/kg时，随着湿空气含湿量的增加，叶片中截面温度又有所下降，当湿空气含湿量从0 g/kg增加到200 g/kg时，叶片靠近前缘的区域温度下降了约50 K。

(a)Tc,t=602 K

图20 Tc,t=725 K时叶片中截面温度云图Fig.20 Wall temperature contours at mid-span section when Tc,t=725 K

2.5 传热关联式拟合

为了准确预测湿燃气传热系数的分布，对湿化燃气轮机叶片冷却结构的设计有着重要的作用。根据前文分析可知，透平叶片燃气侧的换热强度与湿燃气的流动状况、含湿量和进口温度密切相关。流动状况可用透平进口雷诺数Rein描述，含湿量和进口温度的影响可用透平进口普朗特数Pr描述。当Rein从5.3×104增加到9.0×104时，根据经验假设叶片表面换热强度努塞尔数Nu为：

(8)

式中：a、b和c为待求系数。

采用最小二乘法分别对压力面平均努塞尔数(Nul,avg)p和吸力面平均努塞尔数(Nul,avg)s进行拟合，确定的关联式为：

(9)

(10)

其中，5.3×104≤Rein≤9.0×104，0.757 9≤Pr≤0.832 0。

图21为平均努塞尔数的关联式计算值Nul,avge与数值计算值Nul,avgn的相关性分布。由图21可知，关联式计算值与数值计算值的误差在±10%以内，二者吻合良好。