2. 北京航空航天大学 能源与动力工程学院, 北京 100083
高效的冷却设计对先进航空发动机涡轮等热端部件极为重要，而准确预测涡轮内部的流动现象以及温度分布则是进行高效冷却设计的前提，也是当今科学研究的热点问题之一[1]。
目前阶段，工程上预测涡轮部件的温度分布主要有3类方法：经验公式法、弱耦合方法、强耦合方法[2]。现代航空发动机涡轮叶片通常具有复杂的冷却结构，经验公式难以保证精度要求；弱耦合方法中流体域和固体域一般采用不同的计算程序，通过某种方式互相传递边界条件，2个计算域并没有紧密耦合起来，计算效率低且计算精度也不高[3]；与前2类方法相比，强耦合方法，即流/热耦合计算方法[4-5](conjugate heat transfer, 共轭传热)，计算效率最高且结果最为准确，在理论研究与工程实际中均有着非常广阔的应用前景。从20世纪90年代开始，国内外众多学者和研究机构对流/热耦合计算方法进行研究，并开发了许多有代表性的程序[6-10]。近些年来，随着计算机技术及计算流体力学技术的发展，CFD商业软件取得了相当大的进步，也很好的支持流/热耦合计算功能，如Numeca、CFX、Fluent等，许多研究人员利用商业软件对发动机涡轮部件的流动现象和机理进行流/热耦合数值模拟研究[11-13]。
另一方面，作为高速旋转的叶轮机械，发动机涡轮部件具有显著的周期性非定常特性，对流动和传热均有非常重要的影响[14]，典型如动静干涉(包括尾迹和势干涉)以及热斑迁移现象等, 在设计和研究中其周期性非定常特性不可忽视。国内外非常多的学者利用实验或者数值模拟方法进行了动静干涉、热斑迁移等方面的研究[15-18]。
需要注意的是，上文所提到的流/热耦合数值模拟方法的相关研究，即文献[6-13]，全部采用的是基于时间推进的定常计算方法；而动静干涉和热斑迁移的相关研究，即文献[5-18]，则全部是对单一流体域进行计算，即没有采用流/热耦合方法。从已经公开发表的文献来看，基于时间推进的非定常流/热耦合数值模拟方法的相关研究是比较少的，而针对于具有周期性特征的非定常流/热耦合数值模拟研究则更为少见。本文利用带周期性热斑的空心钝头叶片算例对该问题进行了研究，并在时域和频域范畴内对计算结果进行了较为深入的分析。
1 对流换热问题时间精度的校验本文采用自行开发的CFD程序HGFS进行相关计算，该程序的空间离散采用基于非结构化网格的有限体积法，时间推进为全隐式格式，对流通量采用迎风格式，黏性通量采用中心格式，程序的应用及验证可参见文献[19-21]。
无限大平板非定常对流换热算例参考了文献[22]中的例1.3，一块很大的厚度(2δ)为160 mm的碳钢板，均匀的初始温度为298 K，将其置于炉温为1 173 K的炉中加热，钢板在炉中的换热系数α=180 W/(m2·K)，导热系数λ=36 W/(m·K)，密度ρ=7 750 kg/m3，比热容c=470 J/(kg·K)。据此计算1 300 s之后和1 960 s之后平板内部温度分布。
温度分布理论解求解方法参见文献[22]，将t=1 300 s和1 960 s对应的结果列于表 1和表 2。其中x/δ的数值代表平板内部相对位置，等于1为平板表面，等于0为平板中心；θ/θm为不同位置过余温度与平板中心过余温度的比值；θ/θ0为不同位置的过余温度与初始过余温度的比值；T为公式方法计算得到的理论温度。
| x/δ | 1.0 | 0.9 | 0.8 | 0.6 | 0.4 | 0 |
| θ/θm | 0.83 | 0.86 | 0.89 | 0.94 | 0.96 | 1 |
| θ/θ0 | 0.398 | 0.413 | 0.427 | 0.451 | 0.460 | 0.48 |
| T/K | 824 | 812 | 799 | 778 | 770 | 753 |
| xδ | 1.0 | 0.9 | 0.8 | 0.6 | 0.4 | 0 |
| θ/θm | 0.83 | 0.86 | 0.89 | 0.94 | 0.96 | 1 |
| θ/θ0 | 0.249 | 0.258 | 0.267 | 0.282 | 0.288 | 0.30 |
| T/K | 955 | 947 | 939 | 926 | 921 | 910 |
建立二维计算网格，设定时间步长dt=0.5 s，调用HGFS程序的传热模块进行2阶时间精度的非定常计算。图 1a)~1b)为计算得到的平板内部横截面的温度云图。图 1c)~1d)为计算得到的沿钢板厚度方向的温度分布，实心符号为理论解，可以看出计算结果与理论解吻合得非常好。该算例表明，HGFS程序的传热求解模块具有比较好的时间精确性。
|
| 图 1 无限大平板对流换热数值模拟结果与理论解的对比 |
计算模型采用无限长的金属平板，厚0.1 m，宽1.0 m，上方有沿着宽度方向且速度为70 m/s的气流吹过，与平板进行对流换热。初始状态，平板和气流保持在相同的温度300 K。在某一瞬间将气流温度提高到800 K，而平板的下表面温度仍然保持为300 K。平板其他边界视为绝热壁面。研究在之后的一段时间内，气流和平板内部的温度分布随时间的变化规律。平板材料为钢，其导热系数λ=16.27 W/(m·K)，密度ρ=8 030 kg/m3，比热容c=502.48 J/(kg·K)。
二维计算网格如图 2所示，其中上方为流体域，下方为固体域。流体域向前延伸0.1 m作为引导区，为无黏壁面，而流体域和固体域的交界面为流/热耦合的黏性壁面。非定常计算的时间步长(dt)为2.0×10-5 s，气流从进口到出口需要计算约790个时间步。在进行非定常流热耦合计算中，流体域和固体域均采用2阶时间/空间离散精度。
|
| 图 2 平板对流换热算例的计算网格 |
图 3a)~3e)分别为物理时间从0.001~0.016 s的温度分布情况。图 3a)中t=0.001 s，对应的时间步为第50步，高温气体即将到达平板前缘；图 3e)中t=0.016 s，对应的时间步为第800步，高温气体已经到达出口。这一过程仅仅在毫秒量级，固体域温度几乎没有发生变化。
|
| 图 3 平板对流换热算例流体域各时刻温度分布(t=0~0.016 s) |
图 4a)~4f)为物理时间从2.0~60.0 s范围内，6个时刻的固体域温度分布情况。此时流体域已经充满了高温气体，边界层以外的区域温度几乎不变，因此仅对固体域x=0~0.2 m的部分进行放大展示，如图 3e)中的虚线方框所示。
|
| 图 4 平板对流换热算例固体域各时刻温度分布(t=2.0~60.0 s) |
图 4a)的物理时间t=2.0 s，对应的时间步(nstep)为第10万步，固体域温度传导距离仅为0.004 m，最高温度(T)为308.9 K；图 4f)中t=60.0 s，时间步为第300万步，温度传导距离约为0.07 m，仍小于平板厚度。实际上，若要温度传导到平板下方，计算需要约500万个时间步。
该算例比较直观地反映了非定常流/热耦合计算的时间推进过程，以及当前所存在的问题：
1) 在时域方法的框架内，流体域和固体域的物理时间必须是同步推进的，且时间步长的设定需要参考整体的最小时间尺度；
2) 流体域热对流的时间尺度通常远小于固体域热传导的时间尺度，从本算例看流体域时间尺度为10-3 s量级，而固体域时间尺度为秒的量级，时间尺度的差异导致固体域时间推进效率低下，需要用极多的计算步数来弥补。
3 带热斑的空心钝头叶片的数值模拟本算例的空心钝头叶片实际上为气冷涡轮叶片的简化模型，目的是研究在周期性热斑作用下流体域和固体域非定常温度的变化规律。
计算模型及网格如图 5所示，采用二维混合型网格，总共包括3个计算域，左侧部分为固定坐标系下的流体域，相当于静子，进口有热斑，但计算域内没有叶片；右侧线框部分为动叶流体域，相当于转子；U形部分为叶片固体域。
|
| 图 5 空心钝头叶片算例的计算网格 |
叶片长度为0.078 m，厚度为0.009 m，栅距为0.1 m。叶片表面为流、固交界面。进口总压Tt_inlet=101 325.0 Pa，出口静压poutlet=90 000 Pa，内壁面温度Tw=350 K且保持恒定，进口总温Tt_inlet为简单的矩形波分布其中非热斑处为400 K，热斑处为600 K，热斑位置固定。叶片沿着y轴负向移动，速度为0.1 m/s。叶片材料为金属钢。
3.1 固体域及前缘点的温度变化规律分析依据给定的栅距和叶片速度，设定计算物理时间步长dt=0.005 s。本算例计算了足够多的周期数(实际上多达数百个周期)，并将最后2个周期的数据进行处理和分析，如下文的图 6~图 10所示。
|
| 图 6 空心钝头叶片算例流体域各时刻温度分布 |
|
| 图 7 空心钝头叶片算例固体域各时刻温度分布 |
|
| 图 8 前缘点温度随时间的变化 |
|
| 图 9 叶片内部时均温度与定常温度的对比 |
|
| 图 10 截线上时均温度与定常温度的对比 |
图 6展示了流体域在一个周期内6个时刻的温度分布情况。在1/6~3/6周期的时间范围内，热斑位于2个叶片的中间；在5/6周期的时刻，热斑几乎恰好打在叶片前缘。
图 7展示了固体域在一个周期内6个时刻的温度分布情况，将叶片前缘进行了放大，并标出了各时刻前缘点的温度值。与图 6结合起来分析，可以观察到如下现象：1/6周期的时刻，热斑已经扫过并且远离叶片前缘，而此时叶片前缘的温度是最高的，为367.22 K；5/6周期的时刻，热斑几乎恰好打在叶片前缘，而此时叶片前缘的温度却是最低的，为366.14 K。这一现象说明，叶片固体域和流体域温度随时间的变化是不同步的，具体情况还需要进一步分析。
图 8为叶片前缘位置流体域、固体域温度随时间的变化曲线。注意，流体域内的监测点位于温度边界层(厚度约为0.2 mm)之外。从图 8可以得到如下结论：
1) 叶片前缘附近的流体域温度上升和下降的过程都非常迅速，而且高温区持续时间非常短，大约只有0.5个周期。这是因为流体域温度变化主要受到热斑的影响，而热斑宽度很小，通过时间较短。叶片固体域前缘点的温度变化则非常缓慢，上升段和下降段均持续了较长时间。
2) 在5/6~6/6周期的时间内，流体域温度达到最大值600 K，而固体域温度此时近似达到最小值366.1 K；而在6/6~1/6周期的时间内，流体域温度已经下降到最小值400 K，固体域温度才达到最大值367.15 K。在前缘位置，固体域温度的上升存在滞后效应，滞后时间约为1/6个周期。图 8所反映的温度变化规律与图 7和图 8是吻合的。
图 9为叶片内部时均温度与定常温度分布的对比，其中时均温度由非定常计算结果进行傅里叶变换得到，而定常温度由定常计算直接得到，可见二者之间有明显差别。图 9中标出了3条截线，截线上的温度分布曲线如图 10所示，其横坐标为叶片的相对厚度，可以看出定常温度比时均温度略高，且外壁面位置的差值最大，约为0.2 K。
3.2 沿叶片深度的温度变化规律分析如图 11所示，为研究叶片内部温度波动沿叶片深度的变化规律，在截线2上从外到内选择了7个均布的监测点。
|
| 图 11 叶片深度方向监测点的分布 |
非定常计算是以定常计算结果为初场进行的，图 12~16展示了监测点温度随物理时间的变化，其中横坐标的时间以计算周期的形式展示。图 12和图 13是叶片速度(Vblade)为0.1 m/s时的计算结果(即图 6~10对应的算例)，而图 14~16是叶片速度为0.5 m/s时的计算结果。当叶片速度为0.5 m/s，热斑扫过叶片的时间周期为0.2 s，一个周期内仍取200个时间积分点，对应的dt=0.001 s。
|
| 图 12 监测点1~7上温度随时间的变化(叶片速度0.1 m/s) |
|
| 图 13 监测点1, 2上温度随时间的变化(叶片速度0.1 m/s) |
|
| 图 14 监测点1~7上温度随时间的变化(叶片速度0.5 m/s) |
|
| 图 15 监测点1, 2上温度随时间的变化(叶片速度0.5 m/s, 周期0~60) |
|
| 图 16 监测点1, 2上温度随时间的变化(叶片速度0.5 m/s，周期120~180) |
图 12中展示了叶片速度为0.1 m/s时，所有监测点温度随时间的变化，可以看出各监测点近似呈现相同的波动规律，但随着叶片深度的增加，温度波动幅度逐渐减小，监测点1的波动幅度约为0.62 K，而50%深度的监测点4仅为0.01 K。
图 13展示的是监测点1, 2温度随时间的变化，实际上只是把图 12的纵坐标范围进行了调整。可以看出，非定常计算需要至少40个周期才能达到较为理想时间周期解；另一方面，监测点2的温度波动幅度为0.16 K，与监测点1(幅度为0.62 K)相比下降得非常迅速。
图 14展示了叶片速度为0.5 m/s时各监测点的温度变化，与图 12相比，可以看出各监测点的波动幅度均有所减小，这表明当扫过频率增加时，热斑对叶片温度分布的影响是逐渐减弱的。图 15和图 16分别为监测点1, 2在第0~60周期以及120~180周期内的温度随时间的变化曲线，可以看出在计算180个周期之后，温度波动近似满足时间周期性。
1) 在周期性热斑的作用下，叶片内部各位置的温度均会产生周期性波动，但沿着深度方向波动幅度迅速衰减，到达一定深度后波动现象可以忽略，即存在一定的“渗透深度”；
2) 随着热斑扫过频率的增加，叶片内部各位置的波动幅度均减小，例如叶片表面位置，当扫过频率为1 Hz时，监测点1温度波动幅度为0.62 K，而频率增加到5 Hz时，监测点1温度波动幅度减小为0.22 K；
3) 热斑扫过频率的增加使叶片温度变化的周期缩短，为满足时间精度，非定常计算的物理时间步长也要随之减小，但叶片导热速率并未变化，因此需要计算更多的周期才能达到时间周期解。当热斑扫过频率为1 Hz，需要计算约40个周期，共8 000个时间步；而热斑扫过频率为5 Hz，需要计算至少180个周期，共36 000个时间步。
4) 从计算量的角度进一步分析，本算例中热斑扫过频率为5 Hz时，所需计算步数已多达36 000步，而对于更高频率，如105 Hz(接近涡轮转子叶片频率)，粗略估算所需要的计算步数可达108以上，直接进行计算显然不切实际，过高的频率仍是限制基于时域的计算方法应用的重要因素。
图 17a)~b)分别展示了叶片速度为0.1 m/s和0.5 m/s时，叶片内部温度波动的一阶幅值与一阶谐波系数沿叶片深度方向的分布，其结果是通过对截线2上各点在一个周期内的温度波动序列进行傅里叶变换得到的。由于本算例中的温度波以一阶为主导，因此仅对一阶幅值和傅里叶系数进行了分析。图 17a)~17b)中实线为幅值，点划线为同相谐波系数A1，虚线为异相谐波系数B1，用双点线标出了温度波的渗透深度及对应的幅值，其中，渗透深度的定义为边界最大幅值的1.0%所对应的位置。
|
| 图 17 温度波动的一阶幅值与谐波系数沿着叶片深度方向的分布 |
从图 17a)~17b)中的计算结果可以看出，当叶片速度为0.1 m/s，即热斑扫过频率为1 Hz时，一阶最大幅值为0.343 K，渗透深度为叶片厚度的57.9%，约为5.21 mm；当叶片速度为0.5 m/s，即热斑扫过频率为5 Hz时，一阶最大幅值为0.173 K，渗透深度为叶片厚度的24.8%，约为2.23 mm。也可以解释为，当热斑的扫过频率增加到原来的5倍时，叶片内部温度波动的一阶幅值减小为原来的50.4%，而渗透深度减小为原来的42.8%。
图 17a)~17b)与图 12和图 14实际上分别是非定常温度波动在频域和时域范畴内的不同表现形式。时域范畴的结果可以直观反映出各位置温度随时间的变化规律，而频域范畴的结果则可以更直接的反映出幅值衰减规律和渗透深度等信息。
4 结论利用自编CFD程序平台对基于时域方法的非定常流/热耦合数值模拟问题进行研究，对平板对流换热过程和带有周期性热斑的空心钝头叶片进行数值模拟，得到的主要结论如下：
1) 流体域内主要传热方式为对流，温度扰动的传播速度直接取决于流动速度；而固体域内主要传热方式为导热，其速率与温度差和材料属性均相关。平板对流换热过程的数值模拟结果表明，流体域的时间尺度为10-3 s量级，而固体域的时间尺度为秒的量级，温度传导存在巨大时间尺度差异；
2) 空心钝头叶片算例的数值模拟结果表明，叶片固体域内的温度变化相对于流体域存在明显的滞后，即当高温热斑恰好打到叶片前缘时，叶片前缘点的瞬时温度几乎刚好为最低值，而当热斑扫过叶片之后，前缘温度才逐渐增加并达到最大，其滞后时间约为1/6个周期；
3) 热斑扫过频率的增加使非定常流/热耦合计算所需要的周期数大幅增加，当热斑扫过频率为1Hz时，需要计算约40个周期，共8 000个时间步；而热斑扫过频率为5 Hz时，需要计算36 000个时间步。而对于更高频率，如105 Hz，粗略估算所需要的计算步数可达108以上，计算量过于庞大，即过高的频率可能使基于时域的计算方法失效；
4) 从频域范畴对叶片内部温度波动规律进行分析，以基波为例，当热斑扫过频率为1 Hz时，最大幅值为0.343 K，渗透深度为叶片厚度的57.9%，约5.21 mm；当热斑扫过频率为5 Hz时，最大幅值为0.173 K，渗透深度为叶片厚度的24.8%，约2.23 mm。即热斑的扫过频率增加到原来的5倍时，幅值减小为原来的50.4%，而渗透深度减小为原来的42.8%。
| [1] |
林宏镇, 汪火光, 蒋章焰. 高性能航空发动机传热技术[M]. 北京: 国防工业出版社, 2005.
LIN Hongzhen, WANG Huoguang, JIANG Zhangyan. High Performance Aeroengine Heat Transfer Technology[M]. Beijing: National Defense Industry Press, 2005. (in Chinese) |
| [2] |
尹钊, 方祥军, 刘思永, 等. 超跨音内冷涡轮流热耦合数值模拟研究[J]. 中国电机工程学报, 2013, 33(14): 114-120.
YIN Zhao, FANG Xiangjun, LIU Siyong, et al. Conjugate Heat Transfer Analysis of Transonic Internally Cooled Turbine Blades[J]. Proceedings of the CSEE, 2013, 33(14): 114-120. (in Chinese) |
| [3] | YOSHIARA T, SASAKI D, NAKAHASHI K. Conjugate Heat Transfer Simulation of Cooled Turbine Blades Using Unstructured-Mesh CFD Solver[R]. AIAA-2011-0498 |
| [4] | BOHN D, HEUER T. Conjugate Flow And Heat Transfer Calculation of a High-Pressure Turbine Nozzle Guide Vane[R]. AIAA-2001-3304 |
| [5] | GUPTA J A. Application of Conjugate Heat Transfer(CHT) Methodology for Computation of Heat Transfer on a Turbine Blade[D]. Columbus: The Ohio State University, 2009 |
| [6] | BOHN D, REN J. Cooling Performance of the Steam-Cooled Vane in a Steam Turbine Cascade[R]. ASME-GT2005-68148 |
| [7] | KUSTERER K, HAGEDORN T, BOHN D, et al. Improvement of a Film-Cooled Blade by Application of the Conjugate Calculation Technique[J]. Journal of Turbomachinery, 2006, 128(3): 572-578. DOI:10.1115/1.2183314 |
| [8] | RIGBY D L, LEPICOVSKY J. Conjugate Heat Transfer Analysis of Internally Cooled Configurations[R]. ASME-GT2001-0405 |
| [9] | CROCE G A. Conjugate Heat Transfer Procedure for Gas Turbine Blades[J]. Annals of New York Academy of Sciences, 2001, 934: 273-280. |
| [10] |
冯国泰, 黄家骅, 李海滨, 等. 涡轮发动机三维多流场耦合数值仿真的数学模型[J]. 上海理工大学学报, 2001, 23(3): 189-192.
FENG Guotai, HUANG Jiaye, LI Haibin, et al. Mathematical Model of Three-Dimensional Multi-Field Coupled in Gas Turbine Engine[J]. Journal of University of Shanghai for Science and Technology, 2001, 23(3): 189-192. (in Chinese) DOI:10.3969/j.issn.1007-6735.2001.03.001 |
| [11] | YUTING J, QUN Z, PING D, et al. Numerical Simulation on Turbine Blade LeadingvEdge High-Efficiency Film Cooling by the Application of Water Mist[J]. Numerical Heat Transfer Part A, Applications, 2014, 66(12): 1341-1364. DOI:10.1080/10407782.2014.915690 |
| [12] | INSINNA M, GRIFFINI D, SALVADORI S, et al. Film Cooling Performance in a Transonic High-Pressure Vane:Decoupled Simulation and Conjugate Heat Transfer Analysis[J]. Energy Procedia, 2014, 45: 1126-1135. DOI:10.1016/j.egypro.2014.01.118 |
| [13] | HO K S, LIU J, URWILLER C, et al. Conjugate Heat Transfer Analysis of a Cooled Turbine Blade Using Frozen Rotor Approach[R]. ASME-GT2015-43699 |
| [14] | TUCKER P G. Computation of Unsteady Turbomachinery Flows:Part 1-Progress and Challenges[J]. Progress in Aerospace Sciences, 2011, 47(7): 522-545. DOI:10.1016/j.paerosci.2011.06.004 |
| [15] | FUNAZAKI K I, TETSUKA N, TANUMA T. Effects of Periodic Wake Passing upon Aerodynamic Loss of a Turbine Cascade Part Ⅰ: Measurements of Wake-Affected Cascade Loss by Use of a Pnuemetic Probe[R]. ASME-GT1999-093 |
| [16] | ARDEY S, FOTTNER L. Flow Field Measurements on a Large-Scale Turbine Cascade with Leading Edge Film Cooling by Two Rows of Holes[R]. ASME-GT1997-524 |
| [17] |
周莉, 张鑫, 蔡元虎.非定常尾迹宽度对气膜冷却效果的影响[R].中国电机工程学报, 2012, 31(29): 97-102 ZHOU Li, ZHANG Xin, CAI Yuanhu. Effect of Unsteady Wake Width on the Film-Cooling Effectiveness for a Gas Turbine Blade[J]. Proceedings of the CSEE, 2012, 31(29): 97-102(in Chinese) |
| [18] |
李虹杨, 郑赟. 动静干涉对涡轮转子叶片气膜冷却的影响[J]. 北京航空航天大学学报, 2016, 42(1): 139-146.
LI Hongyang, ZHENG Yun. Effect of Rotor-Stator Interaction on Film-Cooling of Turbine Blade[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(1): 139-146. (in Chinese) |
| [19] |
李虹杨, 郑赟. 尾迹对涡轮叶栅边界层转捩的影响[J]. 推进技术, 2017, 38(3): 532-538.
LI Hongyang, ZHENG Yun. Effect of Wake on Boundary Layer Transition of Turbine Cascade[J]. Journal of Propulsion Technology, 2017, 38(3): 532-538. (in Chinese) |
| [20] |
郑赟, 李虹杨, 刘大响. γ-Reθ转捩模型在高超声速下的应用及分析[J]. 推进技术, 2014, 35(3): 296-304.
ZHENG Yun, LI Hongyang, LIU Daxiang. Application and Analysis of γ-Reθ Transition Model in Hypersonic Flow[J]. Journal of Propulsion Technology, 2014, 35(3): 296-304. (in Chinese) |
| [21] | LI Hongyang, ZHENG Yun. Effect of Surface Roughness on Conjugate Heat Transfer of a Turbine Vane[R]. ASME-GT2016-56744 |
| [22] |
曹玉璋. 传热学[M]. 北京: 北京航空航天大学出版社, 2001: 40-44.
CAO Yuzhang. Heat Transfer[M]. Beijing: Beihang University Press, 2001: 40-44. (in Chinese) |
2. School of Energy and Power Engineering, Beijing University of Aeronautics and Astronautics, Beijing 100083, China
