不同燃气输运特性计算方法对涡轮平面叶栅对流换热的影响

薛 钰,刘景源

(南昌航空大学飞行器工程学院,南昌 330063)

以燃料与氧化剂化学放热产生的热能转变为机械功为动力的推进系统一直是研究的热点。研究性能高、质量轻、尺寸小的航空发动机一直是世界各国追求的目标。提高涡轮前燃气温度是达到上述目标的重要因素。而提高涡轮前燃气温度的贡献中,约70%是来自换热和冷却技术的发展[1]。因此航空发动机推进系统内的对流换热问题是研究的重点之一。

文献[1]针对燃气涡轮换热,系统地论述了换热分析及冷却等技术及研究进展。文献[2]对小型航空发动机的换热问题进行了综述,指出了燃烧室、涡轮叶片内外表面等独特的换热问题。文献[3]则对一种自激Helmholtz 脉冲燃烧室的弯尾管的换热问题进行了实验及数值模拟研究,为此类发动机燃烧室设计提供技术支撑。文献[4]应用实验及理论分析研究了燃烧室内在流动方向加装波纹管以提高涡轮发动机燃烧室在高温燃气热冲击下的使用期限。由于在波纹管壁面开孔导致了燃气和冷却空气气膜掺混,因此文献[4]研究掺混状态下冷却效果。关于涡轮叶片对流换热问题,气膜冷却[5-6]及多场耦合[7]是研究的热点。

上述研究的换热高温气体绝大多数采用高温空气代替高温燃气,但是燃烧室内产生的高温气体为多组分的燃气,无论高温燃气与航空发动机部件壁面的对流换热还是冷却空气和高温燃气的掺混,均涉及多组分气体的输运过程。多组分的气体动量输运导致了黏性的产生,而温度的差别导致换热。高温燃气和高温空气输运特性不同以及输运特性计算方法的不同均引起换热的差别。

为了计及多组分燃气输运特性的影响,需要计算燃气的输运系数。而燃气的输运系数由多组分燃气中每种气体的质量分数和每种气体的输运特性计算方法决定。每种组分的质量分数是燃烧室燃料成分、燃烧室出口温度等的函数。并且组分的输运特性为燃烧室出口温度的函数。另外,组分的输运特性的计算方法也有多种。因此研究燃气的输运特性计算方法对涡轮换热的影响,需要考虑燃烧室出口温度不同导致的燃气组分质量分数、每种组分不同输运特性计算方法等对对流换热的影响。发动机不同工作状态下,燃烧室喷入的燃油量不同，燃烧室出口温度、燃气质量分数亦不同。综上,上述不同因素的组合数巨大,应提出一种快速计算不同燃气输运特性计算方法对涡轮换热影响的方法,为工程应用提供参考。

随着数值模拟及计算机技术的不断进步,应进一步深入研究涡轮叶片对流换热问题。因此，在计及燃气多组分换热的基础上,提出一种不同燃气组分输运特性计算方法对对流换热影响的快速算法,分析了燃气各组分的不同输运特性计算方法之间的不同组合对多组分燃气对流换热强度的影响,并以平板及涡轮平面叶栅为研究对象,应用数值模拟及理论分析方法,研究了高温燃气动量输运及热量输运计算方法之间不同组合的差异对对流换热强度的影响范围。提出的多组分燃气输运特性快速计算方法可为燃气涡轮的对流换热问题提供支撑。

1 燃气输运特性计算

所用燃料的化学方程式为肖保国等[8]对RP-3航空燃油的化学动力学模型研究简化后的C9.74H20.52。假定燃料在燃烧内与空气完全燃烧,并且燃烧后燃气组分仅为O2、N2、CO2及H2O等4种。以下仅对上述4种组分进行计算及分析。

1.1 燃气输运特性计算方法

单组分气体分子黏性系数和热导率计算公式有多种,但常用的有简单的幂次律及Sutherland公式[9]。

分子黏性系数μ的幂次律及Sutherland公式为

(2)

热导率k的幂次律及Sutherland公式为

(4)

式中:μ0、k0分别为常温度T0下的分子黏性系数及热导率；S、n为常数,皆取自文献[9]。

由于燃气为多组分混合气体,可应用Wilke[10]建议的公式,计算燃气的分子黏性系数μmix和热导率kmix。

(6)

式中：μi、ki分别为各组分的分子黏性系数和热导率,由式(1)～式(4)计算可得；xi为燃气各组分的摩尔分数，可由文献[11]中的方法计算给出；φij为配分函数；ns为混合气体组分数目。

配分函数的表达式[10]为

(7)

式(7)中：μi、Mi分别为第i个组分的分子黏性系数及摩尔质量。

只需把式(7)中μi的换成ki即为热导率k的配分函数。

所用的第3种燃气不同组分的分子黏性系数及热导率的计算方法由文献[12]给出。

最后一种计算方法由国际工质物性计算软件REFPROP得到[13]。

1.2 燃气输运特性计算结果

图1所示为不同计算方法时燃气4种组分的分子黏性系数μ随温度变化曲线。其中μN为文献[12]给出的值,μP为式(1) 的计算值,μS为式(2) 的计算值,μR是软件REFPROP的计算值(其中只给出了2 000 K的H2O的值)[13]。由图1可见,燃气4种组分的μ皆随温度的升高单调增加,且不同方法给出的各组分μ的差别也随温度的升高而逐渐增大。

图1 不同计算方法时燃气各组分μ随温度的变化Fig.1 Different methods of calculating viscosity coefficient of gas components at different outlet temperatures

图2所示为不同方法给出的燃气4种组分的热导率k随温度变化曲线。其中各热导率k的含义与图1类似。图中燃气4种组分的k皆随温度的升高单调增加,且不同方法给出的k的差别也随温度的升高而逐渐变大。

图2 不同方法时燃气各组分k随温度的变化Fig.2 Different methods for calculating thermal conductivity of gas components at different temperatures

综上,不同计算方法得到的μ和k不同,且之间的差别随来流温度的升高而升高。k不同计算方法之间的差别比μ的大。

1.3 燃气各组分的普朗特数

对于燃气多组分的对流换热问题,动量的运输及热量的传递对对流换热的影响可以用无量纲的Reynolds数及普朗特数表征。

混合燃气的普朗特数为

Pr=μmixcPmix/kmix

(8)

式(8)中：燃气定压比热cPmix由燃气各组分的cPi和质量分数xi加权计算；各组分的cPi则用温度的四次多项式给出。

图3所示为不同计算方法给出的燃气各组分的Pr随来流温度变化曲线。其中PrN为文献[12]中各组分的Pr；PrR为软件REFPROP计算的值(只给到2 000 K的H2O的值)；而μP-kP为式(1)和式(3)计算μ及k后得到的(下标P及S分别为用幂次律及Sutherland公式给出的μ及k)。由于Pr数是μ和k的函数,而不同计算方法得到的μ和k不同,因此不同计算方法组合后得出的Pr数有一定的差别。

图3 不同计算方法燃气各组分Pr随来流温度的变化Fig.3 Different methods for calculating Pr of gas components at different temperatures

1.4 不同燃气输运特性计算方法对对流换热影响的快速算法

虽然燃气中每种组分的输运特性计算有多种组合,由于燃气当地努赛尔数主要由其雷诺数及普朗特数决定(组分扩散系数影响很小),因而给出不同组合下燃气普朗特数随温度变化所有曲线中的上下边界曲线后,即可分析不同燃气输运特性计算方法的不同组合对对流换热的影响程度。

由式(5)～式(7)易证,μmix随μi的增大单调递增,kmix随ki的增大单调递增。另外,从物理上也易得,当每种燃气组分的质量分数(不同来流温度下组分的质量分数可由燃烧室热力计算给出,下同)及燃气温度均给定下,μmix随μi的增大单调递增,kmix随ki的增大单调递增。

由于μmix随μi的增大单调递增,kmix随ki的增大单调递增,因此,由式(8)燃气Pr的上边界由各组分μ的计算方法中最大值和各组分k的计算方法中最小值代入式(5)和式(6)计算得出；燃气Pr的下边界由各组分μ的计算方法中最小值和各组分k的计算方法中最大值代入式(5)和式(6)计算得出。图4为燃气Pr的上下边界随来流温度变化曲线。由图可见,燃气Pr上下边界之间的差别随着燃气温度的升高逐渐增大。

图4 燃气Pr的上下边界随来流温度变化曲线Fig.4 Maximum and minimum values of Pr of the gas at different temperatures

2 平板对流换热数值模拟

2.1 几何模型及边界条件

平板几何模型如图5所示。整个计算区域300 mm×60 mm,其中板长250 mm,前缘点到来流边界为50 mm。

图5 平板几何模型及边界条件Fig.5 The geometry and boundary conditions of flat plate

数值计算采用不可压N-S方程组及燃气多组分输运方程组。平板绕流的雷诺数为25 000,数值模拟假设层流流动。采用速度入口、压强出口及壁面温度为1 000 K的等温、无滑移壁边界条件。来流为多组分燃气,温度分别为1 200、2 200 K(不同来流温度下组分的质量分数可由燃烧室热力计算给出,下同)；平板前缘到来流边界为对称条件。

2.2 计算结果及分析

2.2.1 网格无关性验证

采用网格点数为9×104、17×104、34×104的3套网格,对来流温度为1 200 K时平板对流换热进行网格无关性验证。

图6(a)和图6(b)所示为不同网格数时出口温度分布曲线和当地努塞尔数随Rex的变化曲线。其中图6(b)中Nux为式(9)给出的值[14]。

图6 出口温度分布及当地Nux变化曲线Fig.6 Temperature and Nux distribution at the outlet

(9)

从图6可见，3套网格的出口温度分布曲线和当地努塞尔数随Rex变化曲线几乎重合,表明当网格数大于9×104时,计算结果不受网格数影响。因此以下均采用17×104网格进行计算。图6(b)中当地努塞尔数和式(9)吻合良好,证明采用的数值模拟方案正确。

2.2.2 不同来流温度下燃气Pr数的上下边界对对流换热的影响

图7给出了来流温度分别为1 200、2 200 K时燃气Pr数分别取上下边界时平板出口温度分布曲线。由图7可见,随着平板法向距离的增加,不同来流温度下的出口温度均单调升高。这是由于流体与平板距离越小,受到对流换热的影响越大。图7中燃气Pr的上下边界出口温度分布曲线差别较小,上边界温度曲线略大于下边界,入口温度为2 200 K时两者之间差别略大于入口温度1 200 K的。在固定来流雷诺数不变下,燃气Pr上下边界出口温度曲线差别随温度升高逐渐增大,来流2 200 K时两者最大差别达到2.6%。

图7 来流温度不同时燃气Pr上下边界出口温度分布Fig.7 Outlet temperature distributions of the gas at different inlet temperatures under the upper and lower boundaries of the Pr

图8分别给出了来流温度1 200 K和2 200 K时燃气Pr上下边界的出口速度分布。由于平板壁面黏性滞止,出口速度随平板法向距离的增加而增加。图8中燃气Pr上下边界出口速度分布曲线差别较大,上边界对应的速度比下边界大,而来流温度为2 200 K时两者之间差别比来流温度为1 200 K的大很多。这表明在平板来流雷诺数不变时,燃气Pr上下边界对应的出口速度差别随来流温度的升高而增大,到2 200 K时两者的最大差别达17.9%。

图8 来流温度不同时燃气Pr上下边界出口速度分布Fig.8 Outlet velocity distributions of the gas under the upper and lower boundaries of the Pr of at inlet different temperatures

图9分别给出了来流温度分别为1 200、2 200 K时燃气Pr上下边界的当地Nux随平板当地雷诺数Rex变化曲线。由式(9),当固定Rex时,Nux随着Pr的增大而增大,但由于Pr变化相对Rex较小,所以Nux的变化并不大。图9中平板的Nux随着Rex的增加单调增加,Pr上边界对应的Nux高于下边界,入口温度为1 200 K时两者最大差别为2.7%,2 200 K时则达到4.0%。这表明在雷诺数相同时,燃气Pr上下边界之间Nux的差值随着温度的升高逐渐增大,其变化趋势与图4中Pr上下边界的变化类似。

图9 来流温度不同时燃气Pr上下边界对应Nux的变化Fig.9 Local Nux distributions of the gas under the upper and lower boundaries of the Pr at different inlet temperatures over the flat plate

3 涡轮叶栅对流换热数值模拟

以微小型涡轮发动机涡轮静叶1/2叶高处平面叶栅对流换热为例,应用数值模拟的方法,研究不同燃气输运特性计算方法对其对流换热的影响。

3.1 几何模型及边界条件

涡轮叶栅几何模型及边界条件如图10所示。其中流场入口及出口高度均为12.5 mm,叶栅弦长为14.5 mm。

图10 涡轮叶栅几何模型及边界条件Fig.10 The geometric model and boundary conditions of turbine cascade

数值模拟的控制方程为N-S方程组及燃气多组分输运方程组。入口流场为湍流,湍流模型为可实现性的二方程k-ε方程,湍流度为5%。采用压力基耦合求解器的SIMPLE算法求解数值模拟的控制方程组。叶栅模型采用速度入口、压力出口及1 000 K等温壁面边界条件。入口温度分别为1 200、1 500 K；入口速度根据雷诺数的定义及本节给定的数值10 000,根据给定的入口温度,由1.4节给出的不同温度下的燃气输运特性及燃气的热力参数求得；出口压强为1 atm；流场上下边界为周期性对称边界条件。

3.2 计算结果及分析

3.2.1 网格无关性验证

采用网格数为14×104、18×104、24×104的3套网格验证入口温度1 200 K时的网格无关性。

图11所示为不同网格数时出口温度分布曲线和当地Nux随叶栅Rex(原点为叶栅前缘点)的变化曲线。由图11(a)可见，3套网格曲线给出的数值模拟结果基本重合,而图11(b)和图11(c)中3套网格在叶背和叶盆上的Nux曲线也基本重合。由此,以下计算均采用18×104网格进行。

Fig.11 出口温度分布及当地Nux变化曲线Fig.11 Temperature and Nux distributions of the skin surface of turbine cascade at the cascade outlet

3.2.2 不同入口温度下燃气Pr上下边界对叶栅对流换热的影响

图12分别给出了入口温度为1 200、1 500 K时燃气Pr上下边界对应的出口温度分布曲线。由图12可见,出口温度中间低两边高。这是因为流体流过呈周期性分布的叶片时,会受到与叶片对流换热的影响,离中部叶栅越远,受到的影响就越小。另外,图12中燃气Pr上边界的出口温度低于下边界,入口温度为1 500 K时两者之间差别大于入口温度1 200 K的。这表明燃气Pr上下边界出口温度的差别随温度的升高而增大。

图12 来流温度不同时燃气Pr上下边界出口温度分布Fig.12 Outlet temperature distributions of the gas under the upper and lower boundaries of the Pr at different inlet temperatures over the turbine cascade

图13分别给出了入口温度为1 200、1 500 K时燃气Pr上下边界对应的出口速度分布曲线。由图13可见,出口速度曲线的变化趋势与图12的出口温度分布曲线类似。这是因为燃气组分的μ会随着温度的升高而增加,而来流雷诺数固定,则来流速度会随μ的增加而变大。另外,图13中燃气Pr上边界的出口速度大于下边界,入口温度为1 500 K时上下边界对应的速度之间差别大于入口温度1 200 K的。这表明燃气Pr上下边界出口速度的差别随温度的升高而增大。

Fig.13 来流温度不同时燃气Pr上下边界出口速度分布Fig.13 Outlet velocity distribution of the gas under the upper and lower boundaries of the Pr at different inlet temperatures

图14分别为入口温度为1 200、1 500 K时燃气Pr上下边界对应的叶栅当地Nux的变化曲线。从图14中可见,叶背及叶盆的Nux随叶栅当地Rex的增大而增大,但在叶片尾缘处受附近当地流场变化的影响而变小。图14中燃气Pr上边界对应的Nux大于下边界。入口温度为1 200 K时,两者在叶背及叶盆上最大差别分别为4.74%、4.58%；入口温度1 500 K时两者最大差别为9.15%及9.70%。这表明,燃气Pr不同对对流换热强度的影响较大,且随着温度的升高,燃气Pr上下边界对应的Nux的差别逐渐增大。

Fig.14 来流温度不同时叶背和叶盆上燃气Pr上下边界对应的当地Nux变化曲线Fig.14 Local Nux distributions of the gas under the upper and lower boundaries of the Pr over the Turbine cascade at different temperatures

4 结论

提出一种不同燃气组分输运特性计算方法对对流换热影响的快速算法。虽然多组分燃气的输运特性计算有多种组合,但由于当地努赛尔数由雷诺数及普朗特数决定,因而仅需得到不同组合下燃气普朗特数随温度变化所有曲线中的上下边界曲线,并以上下边界点对应的燃气输运特性的值作为输入条件,由给定的外形进行数值模拟即可评估不同输运特性计算方法的差别对对流换热影响程度。而后以平板和涡轮平面叶栅为例研究了不同来流温度下,不同燃气输运特性计算方法组合对对流换热的影响，得到如下结论。

(1)对平板对流换热,随着平板来流温度的升高(来流雷诺数不变),燃气普朗特数上下边界对应的当地努赛尔数的差别变大。当来流温度从1 200 K升高到2 200 K时,当地努赛尔数差别从2.7%增大到4.0%。

(2)对涡轮平面叶栅对流换热,燃气普朗特数上下边界对应的当地努赛尔数之间的差别随温度的升高而增大(来流雷诺数不变)。当来流温度从1 200 K升高到1 500 K时,当地努赛尔数差别从4.58%增大到9.70%。

(3)在雷诺数相同的条件下,燃气组分输运特性计算方法的不同组合对对流换热强度的影响较大。多组分燃气在燃烧室出口温度不同时各组分的质量分数也不同,并且单组分气体输运特性有多种计算方法,导致在研究热机及推进系统的对流换热问题中的燃气输运特性和Pr时,有众多的计算方法的组合。因此在进行对流换热计算时,可快速给出不同燃气输运特性计算方法所有组合的努赛尔数的变化区间的左右边界,为工程应用提供参考。