0 引言

飞机在低温高湿的气象条件下飞行时，云层中的过冷水滴撞击表面，会发生积冰现象，随即引起一系列安全问题^[1-3]。根据FAA飞行安全部的统计结果，1990—1999年期间，由气象原因引起的飞行事故共3 230起，其中因结冰引起的事故高达388起（占12%）^[4-5]。2003—2008年期间，又有380起与结冰有关的事故报告。由此可见，机翼结冰是影响飞行安全的重大隐患。

空气中的过冷水滴撞击到机翼表面后会在机翼表面形成积冰，实验表明：积冰表面并不光滑，而是由许多微小半球形颗粒覆盖的。半球形颗粒在空间或时间上并不是常量，当过冷液滴以冰点以下的温度撞击机翼表面时，液滴并不会全部冻结，而是液滴内部冻结，外部还是保持液态，随着时间积累液态部分逐渐增多，在气流的作用下与邻近水珠聚结变大。当液态部分增加到一定程度后，空气剪切力大于水的表面张力时，液态部分停止变大。这时，水珠变为流动的水珠或汇集成细流，具体过程如图1所示^[6]。

Fortin等^[7]在Hansman^[8]的工作基础上提出了机翼表面粗糙度的反馈和水脱落计算分析模型。该模型定义了粗糙元素为水珠和冰上的液态水层有规则的接触角，粗糙度高度等于水珠在移动之前达到的最大高度，水珠的高度则是通过分析水珠的生长和作用在水珠上的力来计算的。图2展示了一个具有高度和接触角特征的球形水珠的二维视图。

水珠的增长和聚结过程侧面反映了表面粗糙度随时间的演化过程。表面粗糙度会影响边界层气流分布以及机翼表面热传导，从而反过来影响增生积冰层的形状和大小。因此，粗糙度和积冰之间有一个随时间变化的双向耦合过程^[9]。

Anderson^[10]将粗糙度尺寸与冻结系数和积冰参数关联起来，发现冻结系数越大粗糙度尺寸越小，但其提出的相关性是基于有限的数据，普适性较低。尼古拉兹^[11]在研究粗糙壁面湍流边界层时曾在水管内壁粘贴上粒度均匀的砂砾，然后测量水管的沿程损失和流速的关系，首次提出了粗糙度函数与粗糙高度的经验关系式。后来Connelly等^[12]研究发现，粗糙度函数不仅取决于粗糙高度，同时也取决于粗糙元的形状。并且提出了在粗糙高度K远大于摩擦长度时，可采用“等效砂砾粗糙高度”来模拟粗糙度对流动的影响。但这种方法的主要问题是在定义冰粗糙度的k_s值方面，还没有达成共识。通过k_s预测冰的吸积对于飞机的设计和适航取证非常重要。CANICE^[13]、LEWICE^[14-15]及FENSAP^[16]软件中广泛使用经验关联式求解k_s，以确定转捩位置。但目前大多数k_s经验关联式仅基于有限的实验数据得出，所以其使用范围十分有限。

实际结冰情况下机翼表面粗糙度是水珠、气流、温度等多因素耦合作用的结果，很难用确切的公式描述其变化规律。对不同k_s值进行气动及传热敏感性分析，探究其对积冰的影响规律，为机翼防除冰设计提供一定的理论指导，则具备一定的现实意义。

1 数值模型

本文采用欧拉法求解空气—水滴两相流，空气和水滴均被视为连续相。空气相由质量、动量和能量守恒方程控制，具体控制方程如下：

式中：t代表时间；C_p代表定压比热容；下标v表示气相；T_S，p，ρ_α，$\overrightarrow{U}$，u，v，w分别表示温度、压力、密度、合速度、沿x方向速度、沿y方向速度、沿z方向速度；λ表示流体导热系数；S_U代表动量源项；S_T代表能量源项。

引入水滴体积分数的概念，求解水滴相的连续性方程和动量方程，进而得到水滴撞击特性，具体控制方程如下：

式中：α表示水滴的体积分数；u_d为水滴的速度矢量；u_a为空气的速度矢量；ρ_d为水滴密度；G_d为水滴的重力；K为空气—水滴动量交换系数。

2 几何模型及验证

2.1 几何模型

计算模型选用特征长度为533 mm的NACA0012翼型，具体结构参数如图3所示。

为保证计算精度，机翼边界层处采用结构化网格，且保证机翼壁面第一层网格的y+值小于1，边界层网格数量不少于10层，并将网格高度增长率控制在1.2以内，网格细节如图4所示。

2.2 模型验证

选取AIWT风洞的实验工况进行模拟，具体边界条件如表1所示，实验翼型采用弦长为533 mm的NACA0012翼型。

将模拟得到的冰形（k_s=0.9 mm）与实验冰形进行对比，图5为对比结果。两者整体轮廓相似，均是机翼下缘冰层较长，前缘冰层较厚，只是在波动细节处略有差异。

3 结果与讨论

3.1 传热敏感性分析

图6为光滑翼型（k_s=0）与不同k_s值粗糙翼型的热流密度对比图。不同k_s值热流密度曲线沿翼型的变化趋势相似，均在机翼前缘驻点处达最大值，而后逐渐减小趋于稳定。在k_s∈[0，0.9]mm区间，热流密度随着k_s值增加而急剧变大，峰值最大相差约为5倍。这是由于k_s值增加，会使机翼附近气流在x方向速度减小（如图7所示），y方向速度近乎不变（如图8所示），转捩位置发生前移，增大了空气与机翼的摩擦，从而增强对流换热效应所致。

k_s值也会对模拟冰形产生影响。为了考虑积冰时间对最终冰形的影响，本文探究了27 s积冰时间和360 s积冰时间冰形之间的差异结果。如图9所示，积冰时间为27 s时冰形之间差异较小，随着积冰时间增加，冰形之间的差异会愈加明显。由图10可知，k_s值对冰形高度以及撞击极限都有较大影响。在k_s值较小时，冰层可延伸到机翼前缘的后方；随着k_s值增加，撞击极限会减小。这是因为低k_s值对应的对流换热系数较小，机翼上方的溢流水也更多，在相同结冰环境下，其在弦向的撞击极限也会更长。因此，针对不同的结冰环境进行k_s值敏感性分析，探究k_s值与不同冰形的关联度，对工程模拟具有较强的实际意义。

3.2 气动敏感性分析

将积冰时间27 s和360 s时不同k_s值模拟得到的冰形进行气动模拟，探究k_s值对机翼气动性能的影响。

给出了27 s积冰时间下，k_s值为0 mm和0.9 mm时模拟冰形对机翼气动特性的影响结果。由图11（a）可知，两者升力系数均呈现先上升后下降的趋势。这是由于飞机攻角增加到一定程度后，机翼上方发生气流分离所致，此时对应的攻角即为失速攻角。k_s值为0 mm时，机翼失速攻角在15°左右，升力系数最大可达1.42。当k_s值增加到0.9 mm时，机翼的失速攻角会减小到13°，最大升力系数则减小到1.28左右。由图11（b）可知，两者阻力系数均随着攻角变大而不断上升，k_s值为0 mm时阻力系数最大为0.136，当k_s值增大到0.9 mm时阻力系数峰值可达0.18。图12给出了不同k_s值冰形所对应的升阻比变化曲线，该曲线象征机翼的升力增益效果。由图12可知，两者升阻比均呈现先上升后下降的趋势，k_s值增加会显著降低机翼升力增益效果，两者相差最大约为14%。

图13比较了在k_s值相同时，积冰时间对机翼气动特性的影响。随着积冰时间增加，升力系数峰值减小至0.78，失速攻角减小到8°，阻力系数峰值上升到0.27。这是由于随着积冰时间增加，冰层逐渐增厚，对机翼上方气流的扰动愈加显著所致。由图14可知，随着积冰时间增加，升阻比会进一步下降，对机翼升力增益效果将进一步降低。

4 结论

本文采用数值模拟的方法，对k_s值进行了敏感性分析，探究了其对机翼传热特性的影响，并评估了不同k_s值对应预测冰形的气动性能，得到如下结论：

1）随着k_s值变大，热流密度会急剧上升，不同k_s值对应热流密度峰值最大相差约为5倍。

2）在积冰时间较短时，k_s值对预测冰形影响较小；但随着积冰时间增加，k_s值对预测冰形的影响将愈加显著。在k_s值较小时，冰层可延伸到机翼前缘的后方；随着k_s值增加，撞击极限会逐渐减小。

3）在积冰时间相同时，随着k_s值增加，预测冰形的升力系数及失速攻角减小，阻力系数变大，升阻比也会降低。

4）在k_s值相同时，随着积冰时间增加，冰层加厚，翼型升力系数及失速攻角会进一步减小，对机翼升力增益效果显著降低。

参考文献