翅片阵列对翼面的冲击换热影响数值模拟研究

张天伦， 王克辰， 张 栩， 周文武

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

当飞机穿过含有过冷水滴的低温云层时，机翼表面尤其是前缘会形成积冰，破坏翼型的气动外形，引起飞行性能下降甚至造成机毁人亡.Petty 等[1]总结得出：1975—1988年有超过800起由于飞机机翼结冰造成的安全事故.因此，高效的机翼防除冰技术对飞行安全至关重要.在众多的防除冰方法中，热气膜防除冰是一种常见有效的防冰手段，从发动机引入的高温空气先经调压阀进入笛形管(热气供应腔)，而后通过喷射孔冲击到机翼内表面从而达到防除冰目的.整个过程中，气膜防除冰效果与机翼冲击换热的性能紧密相关，因此提升曲面冲击换热性能一直是相关领域的研究重点.

国内外学者针对平板冲击换热开展了大量的研究.前期研究表明：射流孔的几何形状、射流孔距靶面的距离、靶面的几何结构对射流的流动和换热都有显著的影响.Jambunathan等[2]总结了不同雷诺数、不同喷嘴形状的平板局部和区域平均传热的实验结果和经验关联式.Lytle等[3]利用激光-多普勒测速和壁面压力测量研究了平板射流冲击的流动结构.Narayanan等[4]通过实验研究了3种流态冲击平板时的流场、表面压力和换热率，结果表明换热的二次峰位置与近壁的流动脉动速度变化相对应.Loureiro等[5]考虑了11个不同的实验数据集，总结了经验表达式，其中包括一个分段努塞尔数表达式.Siddique等[6]通过数值模拟确定了不同热流输入边界条件下的局部努塞尔数.除了平板流动和换热机理的研究，部分学者也探索了在平板表面添加粗糙单元(如三角形的翅片[7-9]、带矩形槽道[8-9]、针肋[10-12]、微型肋加糙[13]、球形突起[14])对射流冲击换热的影响.

除了平板模型，半圆形曲面冲击也是重要的研究对象.在改变圆形曲面的几何结构方面，Singh[15]利用ANSYS Fluent 17.2对射流冲击复合凹坑和凸凹表面进行了数值模拟，研究了凹面上凹陷和凸起的交替位置对冲击传热的强化作用.不仅如此，Singh等[16]还建立了一个基于人工神经网络的替代模型，以确定半圆形凸出物的最佳位置，最大限度地提高凹形表面的传热.Gau等[17]通过实验研究了表面曲率对凹、凸表面冲击冷却流动和换热过程的影响.Lyu等[18]测试了4个凹面(包括一个半圆柱凹面和3个具有不同宽深比的抛物线凹面)的射流冲击换热实验结果，探究了曲率对凹表面单排人字形喷嘴射流冲击换热的影响.此外，刘锡晨等[19]进行了圆形和冠齿形射流口的射流实验研究，实验结果表明冠齿形喷嘴可以增大半圆形凹面的努塞尔数.

针对曲面冲击在翼型防除冰的应用研究，Hannat等[20]通过ANSYS CFX数值模拟了三维防冰腔内部的流动与换热，并与测试案例进行了结果对比.张靖周等[21]采用数值模拟的方法研究了射流孔直径、射流孔间距等参数对凹腔表面温度分布的影响.郭之强等[22]研究了光滑表面与具有长方体突起结构对内凹面的流动和换热特性的影响，模拟结果表明凸起结构可以更加均匀地发散壁面射流.归晓烨[23]计算曲面槽道内部强制对流，运用 Gnielinski 公式对机翼防冰系统槽道进行了分析，得到了槽道截面长度和宽度变化对换热系数的影响.关于冲击弯曲表面的数值模拟湍流模型的使用，Seyedein等[24-25]通过研究认为二方程的湍流模型只能模拟出部分的射流流动特性，改进k-ε和k-ωSST模型可以较好地拟合实际的射流冲击曲面的换热和流动特性.

虽然，针对射流冲击机理以及结防冰应用开展了较多研究，然而之前研究大多以平板或半圆形曲面为主，这与真实的翼型曲面相差较大.因此，本文基于真实翼型结构，在几乎不改动翼型设计的基础上，数值探讨如何在内壁面添加微小翅片阵列实现曲面冲击换热性能提高.以平板和NACA0012对称翼型为研究模型，首先探索在平板上添加导流翅片阵列的流动换热性能，而后推广到曲面翼型.相关研究成果能为进一步提高机翼防除冰性能提供方法借鉴和设计参考.

1 射流冲击流动结构

平板冲击射流是一种被广泛研究的流动现象，冲击到平板后从驻点向外扩散的流动通常可以分为3个区域(见图1)：第1个区域为自由射流区，在自由射流区段，来自周围的流体被卷入射流，从而射流速度u0不断降低，其垂直速度方向的射流面积逐渐增大；第2个区域为冲击区，在该区域，垂直平板方向速度快速下降并偏转，静压则快速上升，冲击区在平板上的尺寸大约为两个喷嘴直径大小；第3个区域为壁面射流区，射流完全转向，速度平行壁面，随着流体的发展，逐渐出现速度边界层，壁面射流表现出比平行流动更高的传热水平，这是由于壁面射流和周围空气之间的剪切所产生的湍流，在换热表面被输送到边界层引起的.

图1 射流流动分区Fig.1 Division of jet flow

在冲击区和壁面射流区的过渡段处增加翅片阵列，起到翅片前扰流、翅片段整流的作用，一方面增大驻点附近的湍流水平，提高了驻点区域的换热；另一方面增大远离驻点区域的换热面积和流动速度，提高了远离驻点区域的换热.最终实现整体区域换热水平的提升.

2 平板射流冲击

2.1 平板模型

图2为平板计算域模型及网格划分情况.平板冲击模型和带有的翅片阵列结构如图2(a)和2(b)所示，考虑机翼防冰的笛形腔内的喷嘴长度与喷嘴直径相当，因此喷嘴长度和喷嘴直径d均为4 mm，喷嘴出口距靶面的高度为5d，平板厚度为0.3d.为保证出口流动对冲击区域无影响，整体的计算域尺寸需足够大，设置为60 mm×60 mm.靶面上的几何结构除了基准的无翅片以外，还包括8片、12片直翅片和12片弯翅片3种，其中直翅片的长度、高度和厚度分别为d、0.5d和 0.25d，离驻点的距离为1.25d.弯翅形状为1/4圆弧，圆弧半径为喷嘴直径大小.

图2 平板计算域模型与网格Fig.2 Flat computing domain model and grid

2.2 数值计算

2.2.1计算域与边界条件设置 平板冲击的计算域同模型，采用ANSYS ICEM进行非结构化网格划分，包括流体域和固体域，流体域与固体域的交界面节点一一对应.从图2(c)和2(d)中可以看出，在射流区域、翅片附近进行了网格加密，在靶面的流体一侧增加了壁面棱柱层，棱柱层的第1层厚度为0.005 mm，对应的无量纲壁面距离y+<1，网格质量在0.3以上.生成网格后在Fluent中进行边界条件的设置和湍流模型的选择，为了保证入口雷诺数与真实防冰系统射流雷诺数(Re≈10 000)匹配，入口速度按照不同的雷诺数(Re=5 500，8 600，10 000，15 000)分别设置为u0=23.70，36.518，42.61，63.91 m/s，温度为335 K.流体域的四周为压力出口，平板固体区域的下底面为恒热流密度面，热流密度为 -2 000 W/m2，四周面为绝热壁面，设置参考Dobbertean[9]的模拟过程.空气的物性参数求解使用的是理想气体状态方程.流体域计算所使用的湍流模型为Realizablek-ε，Enchanted Wall Function模型，该模型作为改进的k-ε算法，能够模拟射流撞击，特别是圆形射流，并且为求解靶面近壁处的流动换热，采用增强壁面函数型的近壁面处理.固体域的温度分布求解为边界条件已知的控制体导热微分方程的求解，材料设置为铝.采用二阶迎风格式的Coupled算法，对整个计算域进行了稳态的数值模拟.

2.2.2CFD验证 共设计3套不同数目的网格进行网格独立性的验证，网格数分别为120万、225万和322万.网格独立性验证的结果表明，当网格数大于225万时，壁面上的温度随着网格加密的变化在0.3 K以内，综合考虑计算精度和模拟效率，最终选择网格数目225万的网格.

为了确保网格和湍流模型的准确性，本文将数值模拟得到的驻点位置处的的局部努塞尔数大小与Garimella等[26]通过实验归纳的实验关联式计算值进行了比较.局部努塞尔数的定义式如下：

(1)

式中：h0为局部的对流换热系数；λ为空气的导热系数；qint为交界面处的热流密度；Tint和Tj分别为交界面固体的温度和射流的入口温度.

Garimella等[26]归纳的实验关联式如下：

Nu0=0.462Re0.585Pr0.4(H/d)0.024

(2)

式中：Pr为空气的普朗特数；H为垂直方向上离靶面的高度.式(2)的适用范围为Re=4 000～23 000，H/d=1～5.

本文数值模拟得到4种雷诺数下平板驻点位置处的Nu0与实验关联式所得的结果对比如表1所示.表中：Er为本文模拟所得Nu0与实验关联式所得Nu0之间的相对误差.由表可见，Er在10%以内，其中低雷诺数(Re=5 500)下数值模拟所得Nu0大于实验关联式计算的Nu0，高雷诺数下模拟所得Nu0小于实验关联式的Nu0.这是由于随着雷诺数的增大， 流动再循环和混合更强，下游的热空气回流到驻点附近，导致Nu0的数值随着表面热流密度的变化有轻微的改变[27]，而本文在流固交界面处的热流密度与文献[26]中实验的热流密度不同.

表1 不同雷诺数下驻点位置的努塞尔数Tab.1 Nusselt numbers of stagnation point at different Reynolds numbers

2.2.3结果分析 为了更加综合地评估增加翅片之后的换热效能，流固交界面处的局部努塞尔数是不适用的，参考Dobbertean[9]对有厚度平板的努塞尔数的定义，定义综合的努塞尔数：

(3)

式中：h为综合定义的对流换热系数；qcon为固体外壁面的热流密度，该值在平板各处为一个恒定值；To为固体外壁面的温度.该式从研究对象来看将整个平板的换热进行了考虑，包括外表面的热量流失、交界面的对流换热、交界面到外表面的热传导，因此是一个综合的换热效能.

针对新定义的Nu，首先对比了4种结构在Re=5 500 时的结果，如图3所示.图中：x、y分别为以射流孔圆心为原点，X和Y方向下的距离.由图可见，加8片、12片直翅片相比于不加翅片均可以提升换热能力，并且8片的换热效果最好.这是因为在12片直翅片的结构中，射流在向四周扩散时发生了堵塞，12片相较于8片周向面积减小了17.1%.加12片弯翅片相对于无翅片的换热能力反而减弱，这是由于弯翅片不仅存在直翅片的堵塞问题，而且翅片是弯曲的，气流在每个翅片的吸力面后缘均会出现明显的分离，使得气流更难向整个平板扩展出去，翅片换热面积增大所增加的换热小于因为射流无法扩散而导致的热量堆积，所以换热能力反而减弱了.

针对图3分析得到的结果，后续对不同雷诺数换热性能的分析，仅使用无翅片和增强换热效果最优的8片直翅片的结果进行比较，如图4所示.从图4(a)、4(b)和4(c)、4(d)中最小和最大两个雷诺数工况可以看出，安装直翅片相较于无翅片的平板，Nu都是有提高的，这也意味着恒热流外壁面的温度有提升.

图3 4种结构的努塞尔数分布云图Fig.3 Nusselt number distribution of 4 kinds of structures

图4 不同雷诺数下的努塞尔数分布云图Fig.4 Nusselt number distribution at different Reynolds numbers

从(0，0)驻点位置处沿与X轴正方向22.5°方向取一条线上的Nu，得到一系列点线图，即不同雷诺数下的Nu分布图，如图5所示.图中：r为沿22.5°线上的点到靶面中心的距离；η为相对于无翅片时，有翅片Nu的提升比例.由图可以直观地看出换热效能的提升效果，随着雷诺数的增大，Nu随之增大.相同雷诺数下，在径向距离5d以内，有翅片的Nu均高于无翅片平板的Nu.在低雷诺数(Re=5 500)下，η最大，平均为5.5%左右，其余3个雷诺数下的η在1.3%左右.这是由于在低雷诺数下换热效果由翅片面积带来的提升与翅片增强局部对流换热系数带来的提升都比较明显，而随着雷诺数的增大，后者带来的提升逐渐占据主导地位，所以整体提升效果没有低雷诺数时那么明显.此外，由于翅片的存在导致向外扩散的周向面积减小，从而向外扩散的射流流速增快，换热增强，径向距离在翅片之外的区域换热也有提升.

图5 综合努塞尔数对比Fig.5 Comprehensive Nusselt number comparison

3 翼型曲面射流冲击

3.1 翼型计算模型与网格

热气膜防除冰冲击的壁面是机翼前缘的内表面，在平板上确认增加翅片对于换热有提升，将该结构应用到翼型曲面上进行验证.在翼型曲面模型中，整体计算域的宽度、射流孔距驻点的高度以及翼型曲面板的厚度保持不变，在曲面模型中研究了翅片长度对换热的影响，在原先1倍喷嘴直径长度的翅片基础上，增加了2倍喷嘴直径长度的翅片模型.整体网格如图6所示，图中：S为沿着外壁面上与叶片型线相切的方向.同样在近壁面和翅片附近进行了加密，边界条件和湍流模型的设置与平板的模型相同.

图6 曲面计算域模型与网格Fig.6 Curved surface computing domain model and grid

3.2 翼型曲面冲击换热结果

图7为无翅片、8片1倍喷嘴直径长度翅片和2倍喷嘴直径长度翅片的曲面板外壁面在Re=5 500 工况下的综合努塞尔数分布，这里将翼型曲面按照图6所示的型线S方向进行展平.图中：s为翼型外表面上的点沿型线S方向距翼型前缘的距离；z为以射流孔圆心为原点，Z方向的距离.增加翅片后红色的高温区域明显增多，并且2倍的翅片的效果要优于1倍的.因此，后续不同雷诺数下的流动换热的分析在无翅片和2倍喷嘴直径长度翅片的模型之间展开.

图7 3种模型的综合努塞尔数Fig.7 Comprehensive Nusselt number distribution of three models

翼型曲面板模型同样研究了4种雷诺数下的换热性能.图8对比了不同雷诺数下的无翅片和有翅片的Nu分布，中间两个雷诺数类似，因此选择最大、最小两个雷诺数的结果进行展示.在4种雷诺数的工况下，加翅片阵列模型的Nu在各个位置均要高于无翅片模型的Nu.在高雷诺数Re=15 000的时候，换热效果的提升最为明显，其现象与平板的结果相反.

图8 不同雷诺数下的努塞尔数分布云图(曲面)Fig.8 Nusselt number distribution at different Reynolds numbers (curved surface)

图9显示了从(s,z)=(0, 0)驻点位置处延展的两条曲线上的综合努塞尔数及其有2倍喷嘴直径长度翅片时相对无翅片的提升比例η.从图中可以看出，随着雷诺数的增大，综合努塞尔数相应提升，并且带翅片的换热始终要优于不带翅片的.整体的提升比例在4%～10%左右，比平板的提升更为明显.

结合图8和图9观察到的现象，在翼型曲面上，随着雷诺数的增大，综合努塞尔数变化规律与平板不同，平板在小雷诺数下提升效果更明显，而翼型曲面在高雷诺数下效果更好；此外，整体提升效果也不同，平板提升比例明显小于翼型曲面.造成这种现象的原因有：① 曲面的曲率改变了射流的流场结构，从而对对流换热的效果产生了影响；② 在翼型模型中采用的是2倍喷嘴直径长度的翅片，而在平板上为1倍喷嘴直径长度的翅片，翅片的长度从增大换热面积上有不同的提升效果.

图9 曲面综合努塞尔数对比Fig.9 Curved surface comprehensive Nusselt number comparison

3.3 翼型曲面冲击流场

翅片阵列对流场的影响对于平板和曲面的模型是类似的，均可以提高冲击壁面附近的湍流度.湍流水平的提升可以导致传热率的增大，为了进一步说明换热得到改善的内在机理，对Re=8 600 的工况进行了瞬态的数值仿真计算.图10给出了距离曲面板内壁面0.5 mm处的近壁曲面上的湍动能强度k以及X、Y、Z方向上的速度u、v、w的均方根误差(RMSE)，RMSE越大，代表速度的波动越大，即在数值模拟计算中速度的脉动量越大.翅片的存在使得k增大了1倍多；并且主要改变了Y、Z两个方向上的速度脉动量，这两个方向也正是平行于曲面的两个方向.说明翅片的结构在驻点附近对冲击射流的流场起到了扰流和增强湍流度的作用.

图11为翅片方向上的时均流线图.图中：|V| 为流场合速度的大小.由放大的流线图可见，射流在到达翅片时与翅片发生了撞击，速度方向发生了改变，翅片前缘同时也是四周扩散流动的驻点位置，存在高水平的换热.不过整体流场的结构没有太大的变化，结合图10的信息，可知翅片影响的是每一个时刻的流场，特别是平行于翼型曲面方向上的速度脉动量，而翅片对整体时均流场的影响比较小.

图10 贴近内壁面的速度均方根误差及湍动能(Re=8 600)Fig.10 RMSE of velocity and turbulent kinetic energy close to inner wall (Re=8 600)

图11 截面流线图Fig.11 Cross-section streamline diagram

4 结论

为了提高防除冰目标曲面的射流冲击换热，研究了翅片结构对平板和翼型曲面板上射流冲击换热的影响，对比不同雷诺数下的Nu，研究结果表明：

(1) 直翅片阵列的结构在平板可以提高射流冲击的换热强度，8片的提升效果最好，12片次之，8片1倍喷嘴直径长度直翅片相较于无翅片可以提升1%～6%.

(2) 翅片的结构在翼型板上同样可以提高射流冲击的换热强度，8片2倍喷嘴直径长度长直翅片在曲面板上的提升为4%～10%，换热水平的提升可以提高防除冰效果.虽然增加翅片在平板和翼型曲面上均可以提升换热强度，不过，由于曲面曲率的影响，射流在曲面上的流动传热和平板是有区别的，后续的研究将针对真实的翼型进行探索.

(3) 翅片结构不仅可以增大驻点以及周边区域的湍动能强度从而提升换热性能(湍动能的提升达到了1倍左右)，也可以从增加换热面积的角度进一步提高换热.

(4) 低雷诺数下，翅片通过改变流动结构增加的换热量与增大换热面积增加的换热量都是比较明显的；高雷诺数下，则以改变流动结构的提升为主.