基于磁共振测速的复合冷却涡轮叶片流动分析

2023-08-31赖佑奎马海腾刘一粟欧阳华1

航空学报 2023年14期

赖佑奎，马海腾，刘一粟，欧阳华1，，*

1．上海交通大学中英国际低碳学院，上海 200240

2．上海交通大学机械与动力工程学院，上海 200240

近年来，为追求航空发动机更高的热效率和推重比，高压涡轮叶片的进口温度不断提升［1］，已远超镍基高温合金等单晶叶片材料的熔点［2］。目前主要通过耐高温材料［3］、热障涂层［4］、高效冷却技术来解决这一挑战。其中，涡轮叶片的冷却通常采用气膜冷却［5］、前缘旋流冲击冷却［3］、中弦带肋蛇形通道、尾缘带扰流柱通道、劈缝冷却等组合而成的复合冷却结构［6］。为了获得更高的冷却效率，对涡轮叶片复合冷却结构下的流动进行实验研究和数值仿真具有重要意义。

在实验方面，涡轮叶片的外部流场或简单冷却结构下的流场测量通常是基于光学技术，如粒子图像测速（Particle Image Velocimetry，PIV）［7-8］、多普勒激光测速（Laser Doppler Velocimetry， LDV）［9］等。LDV 多应用于空间单点的速度测量［10］，PIV 用于平面速度的测量［11-12］。对于复合冷却涡轮叶片，由于LDV 和PIV 结构间的相互遮挡而难以对内部的三维流场进行完整的测量。而源于医疗领域的磁共振测速（Magnetic Resonance Velocimetry， MRV）技术为解决这一问题提供了新思路。

斯坦福大学的Elkins 等［13］、Iaccarino 和Elkins［14］通过PIV、LDV 对MRV 在涡轮叶片简单冷却结构流动测量中的准确性进行了验证，并肯定了MRV 的测量效果。Benson 等［15］通过MRV阐明了前缘射流冲击冷却复杂流动的关键特征，同时识别了滞止区和再循环区。Saglam 等［16］的MRV 实验数据表明，冲击孔内流动受通道横流的影响，为数值模拟的验证提供了良好的基础。Elkins 等［17］利用MRV 测量了一个蛇形通道模型的平均质量流量，结果与超声波流量计测量的误差在7%以内。Bruschewski［18］、Wüstenhagen［19］等对MRV 数据进行了空间可视化处理，得到了蛇形通道的三维流线，并通过MRV 数据验证了不同的湍流模型。Baek 等［20］对一种蛇形通道进行了MRV 实验和数值分析，详细阐述了蛇形通道中的二次流和流动不均匀性。Yapa［21］对基于VKI 叶型的带扰流柱和半劈缝的叶片内外流掺混流动进行了MRV 测量，通过分析发现冷流产生的纵向涡在尾迹扩散中起主导作用。Siekman等［22］使用MRV 研究了一种带扰流柱的尾缘结构，与CFD 结果对比发现，CFD 在第一排扰流柱的最大速度预测值较MRV 数据更大。还有很多学者用MRV 测量了带肋三角形通道［23］、格子通道［24］的速度场。关于完整复合冷却结构的研究较少，仅Williams 等［25］测量了基于NACA-0012翼型的带有层板冷却叶片的气膜冷却孔上的速度场。目前，医用磁共振成像（MRI）已经实现脑部、心脏位置血液的瞬态流动显示［26-27］，但现阶段MRV 技术在高压涡轮叶片内部三维的测量还限于稳态流场。综上所述，大多学者的研究局限于1、2 种简单冷却结构的测量，而接近真实叶片的复合冷却的稳态测量研究还很少。

数值模拟方面的相关研究很多，这里重点介绍尾缘扰流柱、半劈缝等结构下流动的相关研究。白琰［28］通过Realizablek-ε模型对梯形通道侧向出流扰流柱阵列的流动进行了详细的分析，并通过气体实验验证了出口通道的质量流量。Bianchini 等［29］用SST（Shear Stress Transfer）模型模拟了错排扰流柱对流动的影响，与PIV 结果吻合较好。吴伟龙等［30］通过k-ωSST模型讨论了尾缘通道中不同大小的扰流柱在旋转状况下对叶片的内部流动的影响，并通过实验数据对比验证了数值方法可靠性。孔星傲等［31］提出了直线式、曲线式2 种倾斜劈缝结构，并通过标准k-ε模型的数值模拟研究了内部流场，得出倾斜劈缝可减小冷气转折角、抑制旋涡产生，使流动更加平缓，从而减小流阻。张丽等［32］对叶片尾缘扰流柱通道进行了计算，结果与流量分配实验吻合较好，另外还发现扰流柱区的流动方向基本为弦向，但也存在着一定的径向流动。综上所述，简单冷却结构下的流动通常可使用PIV 或LDV 进行验证，但复合冷却涡轮叶片的内部流动通常只能使用CFD 进行研究或局部验证。而MRV 为完整流场的实验验证提供了一种新途径。

对于上述问题，本文针对具有复合冷却结构的高压涡轮叶片（美国GE 公司的E3发动机高压涡轮叶片），开展基于MRV 技术的流动测量实验，获取涡轮叶片内部的三维速度数据。同时，通过MRV 实验数据对数值仿真结果进行校验，并着重对尾缘的流动情况进行分析，研究全高扰流柱、半高扰流柱对尾缘出流的影响。

1 MRV 原理

MRV 即相位对比磁共振成像（PC-MRI），是利用射频脉冲使氢质子发生核磁共振，进而得到磁共振图像。其中，三维空间位置的速度信号主要是通过3 个方向的相位编码实现，即利用氢质子在梯度磁场中的进动产生相位差［33-34］。

氢质子自旋角速度ω(r，t)的定义为

式中：r为位置，为时间t的函数，m；γ为氢质子的旋磁比，常取42.5 MHz/T；B0为主磁场场强，ΔB0为局部不均匀磁场的场强，T；G为附加的梯度磁场， T/m。

将式（1）在某时刻t0至下一时刻TE内进行积分，得到

在某一时刻t0，施加一次双极梯度磁场后，可得到相位Φ(1)，然后在下一时刻tE施加一次反向的双极梯度磁场后得到相位Φ(2)。这时，位移未改变的质子相位差为0；而位移改变的质子相位差不为0，且该处质子的相位差为

定义编码速度为

当ΔM1很小时，可将质子自旋的速度转化为3 个方向上的真实速度

式中：vencX、vencY、vencZ分别为X、Y、Z方向的编码速度，m/s；vX、vY、vZ分别为X、Y、Z方向的真实速度，m/s；ΔΦ为相位差，rad；ΔΦ/π 为实验得到的磁共振图像数据。

2 实验设置

2.1 实验系统与装置

本次实验在上海东软医疗科技有限公司的磁共振研发实验室完成。图1 为实验系统，主要包括东软医疗1.5 T 核磁共振仪、测试段、供水系统、数字温度计、浮子流量计等。图2 为测试段，由非磁性材料制作，复合冷却涡轮叶片模型是采用GE-E3叶型［35］，叶片弦长为187.7 mm，叶高213.2 mm，采用SLA 技术3D 打印而成（打印精度为0.1 mm）。其内部冷却结构如图3 所示，包括前缘冲击冷却、气膜孔、带肋蛇形通道、尾缘带全高、半高扰流柱通道、气膜孔和半劈缝等结构［36］。

图1 MRV 实验系统［36］Fig.1 MRV experiment setup［36］

图2 涡轮叶片模型和测试段Fig.2 Turbine blade model and test section

图3 涡轮叶片内部冷却结构Fig.3 Internal cooling structures of turbine blade

2.2 实验工况

根据美国航空航天局（NASA）试验报告［35］中第11 个状态点的试验数据和几何参数，计算第一级动叶的前缘和尾缘通道的冷气质量流量和雷诺数。尾缘通道雷诺数范围一般为1 000～60 000［6，37-40］，本次实验为保证得到较好的测量效果，取较低的雷诺数（2 718）进行，同时采用含0.015 mol/L 的硫酸铜水溶液作为本次实验工质。表1 列出了NASA 试验工况与本次实验工况对比。

表1 实验工况Table 1 Test conditions

基于水力直径的雷诺数定义为

式中：ρ、μ分别为密度、截面平均速度、动力学黏度；DH为截面水力直径，DH=4A/C，A、C分别为通道截面积和湿周长。

2.3 实验方法

测试段中心与主磁场线圈中心重合，并将接收信号的相控线圈包裹在测试段上，以得到更好的测量效果。考虑到温度对雷诺数的影响较大（22～25 ℃时，每提升1 ℃，雷诺数增加约2.3%），为减小温度变化对雷诺数的影响，将温度控制在22 ℃左右。核磁共振仪的扫描参数如表2 所示［36］。在流动情况下扫描3 组，在不流动的情况下扫描一组。其中，为避免出现混叠［41］，编码速度主要由CFD 结果的全局速度最大值进行初步确定。另外，通过对MRV 装置进行静态水模测试，得到的平均速度误差在2.62 cm/s 以内。

表2 MRV 测量参数［36］Table 2 MRV measurement parameters［36］

2.4 实验数据处理

本次实验得到的是稳态流场数据，磁共振图像数据（Dicom 文件格式）主要包括幅值图和左右RL、前后AP、头足FH 这3 个方向的相位图。在对流动情况下的3 组数据取平均后，需将其转换为3 个方向的坐标值及对应的速度值(X，Y，Z，vX，vY，vZ)。横断位的单切片层的坐标映射公式及速度转换公式为

式中：XX、XY、XZ为绕X方向旋转的余弦值；YX、YY、YZ为绕Y方向旋转的余弦值；Δi、Δj分别为X、Y方向分辨率；SX、SY、SZ为扫描中心点离坐标系原点的距离；i，j分别为单层图片的纵向和横向的像素点坐标；v(i，j)为该像素点的速度，m/s；D(i，j)为该像素点的值。其中，余弦值XXYZ、YXYZ，距离SXYZ及分辨率Δi、Δj由东软医疗提供。另外，冠状位计算公式与之类似。

速度幅值V的计算公式为

涡量ωY的计算公式为

其中，质量流量̇的计算公式为

不确定度σv的估计公式［17，34］为

式中：SNR 为信噪比，噪声值为背景区域的平均幅值，信号值为感兴趣区域的平均幅值。横断位和冠状位数据的SNR 分别为33.3、20.1［36］。

3 数值仿真设置

3.1 计算域和边界条件

基于商用软件ANSYS CFX 的k-ωSST 湍流模型进行雷诺平均数值模拟（RANS）计算，采用高精度的对流格式。提取实验叶片模型的内部流道，作为涡轮叶片整体三维模型Case0 的计算域，具体如图4 所示。设置水流进口流量ṁ为0.287 kg/s（前缘通道ṁLE、蛇形通道ṁMC、尾缘通道ṁTE均已标出），前缘LE1～LE21 出口静压为102 140 Pa，蛇形通道MC1～MC3 出口静压为102 267 Pa，尾缘TE1～TE10 出口静压为102 290 Pa（由于本次实验未测量出口测压管水头，为保证基于各个通道水力直径的雷诺数与实验一致，让前缘通道Z/D=17、蛇形通道Z/D=18、尾缘Z/D=14 这3 个通道处的流量与用MRV 数据计算的流量值保持一致）。温度为22 ℃，壁面均为绝热无滑移壁面。进口湍流度为5%。稳态计算的收敛标准为通量残差的均方根值（RMS）小于10-4。流量进口和压力出口位置已在图4 中标出，虚线箭头表示水流在前缘、蛇形通道、尾缘3 个通道中的流动方向，其中，扰流柱直径D为6.3 mm。

图4 Case0 计算域和边界条件Fig.4 Computational domain and boundary conditionsof Case0

图5 为Case0 计算域简化后的尾缘计算域Case1，尾缘处边界条件同Case0。简化的原因是前缘和蛇形通道的压力出口不利于后续CFD的工况设置（需保证尾缘通道的进口流量为0.058 kg/s），同时还能减小尾缘部分各工况的计算量。

图5 Case1 计算域Fig.5 Computational domain of Case1

为讨论全高、半高扰流柱对尾缘出流的影响，对4 组不同扰流柱配置进行了数值模拟研究，如图6 所示。图中，Case1 为基础构型，作为对照组；Case2 仅有高扰流柱排；Case3 仅有半高扰流柱排；Case4 无扰流柱。

图6 不同扰流柱结构设置Fig.6 Different pin fins settings

3.2 网格划分

采用商业软件ANSYS ICEM 对Case0 进行网格划分，流体区域的主体部分采用四面体网格、近壁区域采用棱柱层网格，并对气膜孔、冲击孔、扰流肋、扰流柱、半劈缝处进行局部加密。近壁面区域设有10 层边界层网格，首层网格高度为0.002 mm，网格增长率为1.2，保证y+＜1，前缘气膜孔、尾缘半劈缝处的网格细节见图7。

图7 Case0 网格细节Fig.7 Mesh detail of Case0

为验证网格独立性，对Case0 的近壁面网格进行加密，得到网格节点数为963 万、1 298 万、1 633 万的3组网格。计算得到尾缘TE9处的质量流量，见图8，当网格数量从1 298 万升至1 633 万后，尾缘TE9 处的流量差值＜0.5%。故采用网格节点数为1 298 万的网格进行计算，Case1（基础构型）的网格节点数为200 万。

图8 网格无关性验证Fig.8 Gird-independence validation

4 结果分析与讨论

4.1 全流域流动情况

根据式（12），计算得到冠状位数据的不确定度σv为1.2 cm/s，约占入口速度的3.7%［34］。

图9 为叶片内部流场的速度分布。在前缘位置，水流经过前缘通道到达冲击孔处，随着Z向距离增加，流速逐渐减小，并在叶顶Z/D=45～50处，几乎发生停滞。而在冲击孔处，流速显著增加，并喷射到涡轮叶片前缘内壁上，形成冲击冷却，之后经过气膜孔喷射到叶片外部。

图9 MRV 与CFD（Case0）速度云图对比Fig.9 Comparison of the velocity contour between MRV and CFD （Case0）

在蛇形通道处，水流经过直通道后到达Z/D=45 处的180°拐角（Z/D=50 处的2 个气膜孔有着全局最大的射流速度），先在拐角前进行了加速，然后出现了流动分离，水流直接冲击到第2 个直通道壁面上，由于Z向的速度分量未消失，使得较高速度的水流一直贴着通道壁流动。在Z/D=10 处的180°拐角，流动现象与Z/D=45 处的比较类似，但分离区域的形状不相似，主要由于拐角处流道的几何形状不同。流过Z/D=10 处拐角后，水流并未贴着第3 个直通道壁流出，可能是受到了V 形肋产生的二次流的扰动。

在尾缘处，水流以一定角度流入，并在尾缘入口形成了一个分离区域，后被Z/D=17，19 处的扰流柱分成了3 部分，靠近蛇形通道的部分贴着通道壁流到了叶顶处，中间部分则流到了半劈缝处，靠近气膜孔的一部分几乎全部从气膜孔和Z/D=32 处的半劈缝流出。

图9（a）、图9（b）对比表明，CFD 结果与MRV 数据在前缘、蛇形通道、尾缘处的流场分布基本吻合，可在此基础上开展尾缘部分的流动模拟，并进一步分析全高、半高扰流柱对尾缘出流的影响。

4.2 尾缘流动分析

图10 为尾缘部分的几何结构。中截面定义为图10（a）中Y1=0 处的截面。其中，全高、半高错排扰流柱的直径D均为6.3 mm，沿X1方向分布的全高扰流柱高度分别为H1/D、H2/D、H3/D分别为1.16、0.89、0.6，半高扰流柱分别为0.58、0.45、0.3，X1、Z1方向的间距直径比SX/D均为2.06，TE1～TE6 为直径3.25 mm 的气膜孔，TE7～TE9 处是高3.25 mm、宽20 mm 的半劈缝，TE10 处是直径3 mm、射流角度为30°的气膜孔。

图11 给出了尾缘中截面的速度云图，通过MRV 结果和Case1 结果的对比表明尾缘的整体流场较为吻合。气膜孔TE1、TE3、TE5 处的流动情况较为相似，扰流柱的存在使得气膜孔入口有较大的速度；在气膜孔TE2、TE4、TE6 处，受扰流柱后低速区的影响，气膜孔入口处的速度略低于其前一个气膜孔；在半劈缝TE7 处（Z1/D=18.49）存在明显的流动分离区域。

图12 为尾缘中截面的涡量云图。在气膜孔TE1～TE6 处，无宏观上的漩涡，一正一负的涡量是由壁面的黏性应力引起。而在半劈缝（Z1/D=18.0 处）的正涡量值明显大于Z1/D=21.0 处的负涡量值，与Z1/D=25.0，30.0 处情况不太一样，再结合图13 可发现，Z1/D=18.0 处的涡量分布现象应该是由宏观漩涡和黏性应力共同产生的。从定性的角度来说，MRV 数据与CFD 结果在尾缘中截面处的涡量分布趋势吻合较好。

图12 尾缘中截面（Y1/D=0）涡量分布Fig.12 Vorticity distribution of middle cross-section of TE （Y1/D=0）

图13 半劈缝处流线分布（Z1/D=18.49）Fig.13 Streamline distribution at slot （Z1/D=18.49）

图14 为MRV 数据与CFD 定量对比，MRV数据与Case1 在Y1/D=0，Z1/D=8.3 处的速度分布趋势较为接近，但在扰流柱后的低速区即X1/D=0.5，4.0，5.5 处的分布存在较大差距。

图14 Y1/D=0，Z1/D=8.3 处线速度分布Fig.14 Velocity distribution of line at Y1/D=0， and Z1/D=8.3

图15 为Y1/D=1.39，Z1/D=5.49 处流线分布，在扰流柱后存在2 个反向的漩涡，表明此处发生了明显的圆柱扰流现象。可判断出此处一正一负的涡量分布是漩涡引起的微团旋转运动。

图15 扰流柱处流线分布（X1/D=1.39，Z1/D=6.48）Fig.15 Streamline distribution at pin fin （X1/D=1.39，and Z1/D=6.48）

图16 展示了半劈缝处的速度分布。在TE7、TE8 处，由于扰流柱的影响，存在2 个速度较高的区域，呈现出不均匀的分布。而在TE9 处，速度分布非常均匀，且CFD 的结果与MRV 数据较为吻合。从TE7～TE9 的变化趋势上看，半劈缝截面的速度分布沿展向愈发均匀。

图16 半劈缝处速度云图（X1/D=7.06）Fig.16 Velocity contour of slots （X1/D=7.06）

图17 给出了气膜孔TE1～TE6、TE10 和半劈缝TE7～TE9 的质量流量。MRV 数据在气膜孔TE1～TE6，半劈缝TE7、TE8 处的流量略小于CFD 结果，而在半劈缝TE9、气膜孔TE10 处大于CFD 结果。尾缘整体出流趋势较为一致：靠近叶根的气膜孔的流量最大，沿展向逐渐减小；靠近叶根的半劈缝的流量最大，沿叶展向逐渐增大。通过MRV 数据计算得到的TE1～TE10 处的质量流量总和为55.44 g/s，比尾缘通道Z/D=14 处的质量流量小4.4%，局部出流位置存在较大的误差。引起误差主要原因是：MRV 数据图片分辨率为0.5 mm×0.5 mm，而气膜孔的直径和半劈缝的高为3.25 mm，像素点较少，边缘位置识别不清，使计算得到的流量存在一定误差。