结冰是影响飞行安全的一大隐患[1]。通常认为，结冰发生的原因是云层中含有温度低于冰点的亚稳态液态水，当飞机穿越由这些过冷水滴组成的云团 时，便会在迎风部件表面发生结冰现象[2]。研究表明，飞机结冰会破坏飞机的气动外形，影响稳定性和操纵性；严重的结冰现象极有可能会造成机毁人亡的重大飞行事故[3, 4]。因此，有效的防/除冰系统对于结冰环境下的飞行安全具有十分重要的意义。
防冰系统能保证部件表面温度维持在冰点以上。当云层中的过冷水滴撞击表面时，若防冰系统提供的热量不足以完全蒸发表面撞击水，未蒸发的液态水会在气动力作用下沿部件表面向后流动，形成溢流水。溢流水的存在及其在表面呈现的不同流动形态会直接影响部件表面防冰热流的分布情况，尤其在防冰不足的情况下，溢流水会对表面结冰区域和积冰冰形产生决定性影响。Messinger模型[5]作为结冰与防冰表面流动换热计算分析的基础模型，考虑部件表面的质量守恒和能量守恒，但并未引入溢流水流动的因素。此后的研究工作也大多围绕结冰预测及防冰分析而展开[6-13]。针对液膜流动及破裂问题，研究工作主要集中于降膜流动及破裂问题[14-17]，考虑飞机防冰表面的水膜破裂及溪流形态的文献还十分有限。Al-Khalil等人[18]曾对防冰表面的水膜流动进行数学建模，耦合了部件表面的能量平衡及溢流水的质量和能量平衡，并提出采用最小能量法预测水膜破裂。Silva等人[19-21]同样将热力学第一定律应用于翼型表面及溢流水的分析、建模，其重点在于对表面温度分布的预测；针对水膜破裂问题，沿用了Al-Khalil的最小能量法，但提出了湿润系数的不同定义。总体而言，已有模型的计算结果与实验数据相比还是存在一定的误差，同时在预测液膜破裂内在机理等方面仍存在一定的局限性。
本文主要关注溢流水在防冰表面的流动换热以及水膜破裂等问题。基于空气-水两层相互作用的质量、动量和能量守恒,建立防冰表面溢流水水膜流动换热及破裂的数学模型，计算分析了来流速度对表面溢流水分布、流态及表面换热的影响，获得了防冰表面连续水膜厚度、破裂处溪流厚度和宽度，以及连续水膜表面主要散热热流的分布情况，并分析了溢流水及其流态对表面换热的影响。
1 物理描述云层中的过冷水滴撞击飞机防冰部件表面，受表面加热作用蒸发，部分未蒸发的撞击水在气动力作用下形成溢流水，沿部件表面向后流动。溢流水在流动过程中会形成连续水膜，但通常情况下，这种流动是不稳定的。
在撞击区域内，由于水滴的收集量相对较大，部件表面通常完全湿润并形成连续流动的水膜。但在撞击区域之外，随着水膜不断向后流动，在表面张力、气动剪切力、壁面粘滞力和压差力等力的作用下以及蒸发、对流换热等传热传质过程的共同作用下，水膜在流动过程中会发生动量和能量损耗。在此过程中，水膜厚度会不断发生变化，而其流动状态也会变得不稳定，易发生破裂形成溪流。此时，部件表面将分为干、湿两个区域，如图 1 所示。
|
| 图 1 防冰部件表面溢流水流动形态示意图 Fig. 1 Runback water flow on anti-icing surface |
1.1 连续水膜流动换热
如图 2所示，部件表面形成连续水膜流动，流动过程中存在质量、动量和能量的交换。假设部件表面水膜流动及空气边界层均为不可压缩、定常、层流流动。由于水膜很薄，且空气中液态水含量(LWC)很小，因此仅考虑空气对水膜的单向耦合作用。相对空气剪切力的大小，可忽略水膜重力影响，而仅考虑空气剪切力及表面压力差对水膜流动的驱动效应。当防冰系统运行足够长时间后，部件表面温度始终保持在0℃以上，因此本文仅考虑无结冰情况下的稳态水膜流动。
|
| 图 2 防冰部件表面水膜流动示意图 Fig. 2 Water film and air boundary flow on anti-icing surface |
图 2中，mevap为水膜表面蒸发质量，mimp为水滴撞击质量，min和mout分别为流入和流出控制体积的质量，Qevap为蒸发散热，Qimp为水滴撞击动能，Qin和Qout分别为控制体积进出口焓值，Qc为表面对流换热，Qd为加热撞击水滴所损失的热能，Qw为壁面导热。
1.2 水膜破裂连续水膜在流动过程中易发生破裂，在表面形成若干条溪流，如图 3所示为假定的破裂处的溪流模型在y-z平面上的横截面形状。
|
| 图 3 溪流模型示意图 Fig. 3 Rivulet model |
图 3中，δ0为溪流厚度，W为溪流宽度，L表示单条溪流及其临界干表面的总宽度。假设溪流截面为圆的一部分，其半径为R，因此在气-液界面上溪流的几何外形可以表示为：
(1) 式中，θ为液相与气相之间形成的角度，即接触角。如图 3中所示，气-液-固三相接触点存在力平衡条件，即Young方程[22]：
其中，σsg、σsl、σlg分别表示固-气界面、固-液界面、液-气界面张力。
根据几何关系，溪流的厚度和宽度可以表示为：
(3) 
(4) 取控制体积内的水膜为研究对象，应同时满足质量守恒、动量定理及热力学第一定律，得到水膜流动的连续性方程、动量方程及能量方程：
(5) 
(6) 
(7) 其中，δ1和δ2分别表示水膜厚度和空气边界层厚度。p为机翼表面静压力；τ0为壁面粘滞应力，τair为气-液界面空气对水膜的剪切应力，剪切应力由τ=μ(du/dy)确定，μ为动力粘度；uδ1为气-液界面水膜速度；u∞,s为水滴撞击表面时沿s方向的速度；ρwater、λwater、Cp,water均为水的物性参数，分别表示水的密度、导热系数及定压比热容；hc为对流换热系数； Tw-s为气-液界面水膜表面温度；T∞为自由来流温度；Td为碰撞水滴的温度，可以视为与来流温度相等；Le为液态水的蒸发潜热；uwater和Twater分别为水膜边界层的速度分布和温度分布。根据边界条件可以表示为:
(8) 
(9) 其中，Twall为防冰部件表面温度；uδ1、δ1、τair、Tw-s均为关于s的函数。
考虑空气对水膜的单向耦合作用，在连续水膜流动换热模型中须引入空气边界层的动量方程。取空气边界层的控制体积为研究对象。与水膜控制体积的分析方法相类似，其动量方程为：
(10) 式(10)中，τ′air为气-液界面水膜对空气的粘滞力，与τair互为作用力与反作用力；ue为空气边界层靠近主流位置的速度；ρair为自由来流空气密度；uair为空气边界层的速度分布，同样可以根据边界条件求得：
(11) 在水膜破裂处，水膜与溪流之间存在质量守恒和能量守恒关系:
(12) 
(13) 其中，下标f和r分别指代水膜和溪流。
根据上述两个守恒方程,可以化简得到破裂点处无量纲控制方程：
(14) 其中，
(15) 
(16) 
(17) 
(18) 上式中δf为破裂处临界水膜厚度，δ+为对应的无量纲值,R+为无量纲溪流截面半径。
在给定破裂位置的条件下，根据式(14)可求得破裂位置溪流截面半径的无量纲参数，进而可以获得溪流的几何外形参数。
对于连续水膜流动换热及水膜破裂模型的详细推导过程可参考文献[23]。对于水膜破裂位置的确定，可以按照文献[18]所提出的最小能量法得出破裂处的临界液膜厚度。对于本文的算例而言，撞击区域内形成的连续水膜厚度均无法达到临界液膜厚度，可以认为在撞击极限处水膜将发生破裂。目前看来最小能量法得到的破裂点与实验观察到的水膜破裂位置还存在很大误差，还需要进一步开展研究工作。
2.3 防冰表面热流水膜破裂后，防冰部件表面被分为干表面、湿表面，其换热特性与连续水膜流动时的完全湿润表面以及无溢流水时的干表面有较大区别。若不考虑溢流水的存在，即撞击极限之外，完全为干表面，此时表面无蒸发热流且对流换热仅考虑部件表面与空气之间的热传递过程。若不考虑水膜破裂，即部件表面始终由连续水膜覆盖，完全为湿表面，此时仅考虑水膜表面的蒸发热流以及水膜与空气之间的对流换热。当考虑水膜破裂时，翼型表面被分为干/湿两个部分，此时仅湿表面存在蒸发热流，而对流换热热流必须综合考虑湿表面上溪流和空气之间的对流换热以及干表面上壁面与空气之间的对流换热。
因此，三种不同溢流水分布形态条件下，表面蒸发热流及对流换热热流可以分别根据以下公式进行计算：
(1) 完全湿润表面(即连续水膜表面)：
(19) 
(20) (2) 干/湿表面(即发生水膜破裂的表面)：
(21) 
(22) (3) 干表面(即无溢流水表面)：
(23) 
(24) 其中，F为水膜破裂处表面湿润系数，可表示为
(25) 求解表面水膜流动换热及水膜破裂模型的前提是部件表面空气流场及水滴撞击特性的计算结果。因此，先通过CFD软件及水滴撞击计算程序获得部件表面空气流场和表面水滴撞击特性。在Fluent中计算空气流场采用k-ε湍流模型；进出口边界分别采用速度入口和压力出口；翼型壁面采用无滑移边界条件。空气流场迭代收敛后，采用欧拉法计算部件表面局部水收集系数。外流场和局部水收集系数的计算结果将作为本文模型求解必要的边界条件。
水膜流动换热及破裂模型的控制方程均在Matlab中编程求解。首先迭代计算连续水膜控制方程；连续水膜模型迭代收敛后，可获得水膜表面各项热流值；其次，在给定破裂位置可对水膜破裂模型进行求解，并获得溪流的厚度和宽度，以及表面换热。整个计算过程如图 4所示。
|
| 图 4 计算过程 Fig. 4 Computational procedure |
4 计算结果及分析
Zhang K等人[24]利用DIP技术对NACA 0012翼型表面液态水的流动形态和水膜/溪流厚度进行了实验测量。实验中NACA 0012翼型弦长为0.101m，攻角0°。采用本文所述的数学模型针对该实验模型进行了连续水膜流动换热及水膜破裂的计算分析。在实验给定的低速流动条件下，连续水膜厚度及溪流厚度和宽度的计算结果在数量级和变化趋势上均与实验数据吻合较好。限于篇幅，本文不再详述，具体对比结果参见文献[23]。
防冰表面溢流水的流动受气动力和表面张力的影响最大，气流速度直接影响水膜及溪流所受气动力的大小。云层参数不变的条件下，防冰部件表面的流动状态由于飞行速度的影响会呈现不同的特性，从而影响表面换热过程。本文对不同来流速度下水膜流动及破裂情况进行了研究，分析了破裂前后防冰表面蒸发及对流换热的特点。 计算仍采用弦长为0.101m，攻角0°的NACA 0012翼型的几何模型，具体的计算状态参数如表 1所示。
| Twall/℃ | u∞/(m·s-1) | T∞/℃ | LWC/(g·m-3) | MVD/μm | |
| 状态1 | 4 | 40 | -5 | 1.0 | 20 |
| 状态2 | 4 | 60 | -5 | 1.0 | 20 |
| 状态3 | 4 | 80 | -5 | 1.0 | 20 |
| 状态4 | 4 | 100 | -5 | 1.0 | 20 |
4.1 水滴撞击特性
不同来流速度条件下翼型表面的局部水收集系数分布如图 5所示。从整体上看，随着来流速度的不断增加，翼型表面的局部水收集系数无论在峰值还是在撞击区域上都有所增加，表明来流速度高时将有更多的液态水到达翼型表面。
|
| 图 5 不同来流速度条件下局部水收集系数分布 Fig. 5 Local collection efficiency at different free stream velocities |
4.2 连续水膜厚度分布
如图 6所示为不同来流速度条件下的连续水膜厚度在翼型上表面的分布情况。从整体上看，在不同来流速度条件下，连续水膜沿翼型表面的分布均呈现逐渐增大的变化趋势。这是由于在撞击区域内，水滴撞击质量远大于表面蒸发质量，导致溢流水量逐渐变大，从而使得水膜在沿表面流动过程中逐渐增厚。此外，从图 6中还可以看出，随着来流速度的增大，表面同一位置处水膜厚度有明显的减小。这是由于前缘处水膜受气流法向作用力影响较大。当来流速度较大时，作用在水膜表面的法向气动力也较大，使得水膜变薄。
|
| 图 6 不同来流速度条件下连续水膜厚度分布 Fig. 6 Water film thickness at different free stream velocities |
由于本文将水膜破裂点设置在水滴撞击极限位置，由局部水收集系数的分布图(图 5)可以看出，来流速度越大，撞击区域越大，即水膜破裂位置随来流速度的增大而延后。从图 6中可以看出，当来流速度较大时，连续水膜在表面上的分布范围更大，破裂位置更加远离驻点，此时破裂处的连续水膜厚度相对较大。这一特性将对破裂后溪流形态产生重要影响。
4.3 连续水膜主要热流分布蒸发换热与对流换热是水膜表面主要的热量损失[25]。图 7和图 8分别为不同来流速度条件下水膜 表面蒸发热流与对流换热热流沿表面分布图。从图中可以看出，随着来流速度的增大，两项热流均有明显增加。这是由于当来流速度增加时，表面对流换热系数会随之增大，进而影响表面蒸发和对流换热热流。
|
| 图 7 不同来流速度条件下水膜表面蒸发散热热流量 Fig. 7 Evaporation heat fluxes on water film surface at different free stream velocities |
|
| 图 8 不同来流速度条件下水膜表面对流换热热流量 Fig. 8 Convection heat fluxes on water film surface at different free stream velocities |
4.4 水膜破裂及溪流形态
图 9和10所示为不同来流速度条件下，水膜破裂处(即水滴撞击极限处)溪流的厚度及宽度的分布情况。从图上可以看出，相比于连续水膜，在相同来流速度条件下，破裂处溪流的厚度和宽度都有所增加；同时，当来流速度不断增大时，溪流的厚度和宽度也都随之增大。正如上文所述，来流速度增大，表面溢流水量增加，在破裂处连续水膜厚度相对较大，因此溪流的厚度会有所增加。同样，溪流宽度也随之增大。但从计算结果的数值来看，几个不同来流速度条件下的溪流厚度和宽度差别并不是很大。
|
| 图 9 不同来流速度条件下破裂处溪流厚度分布 Fig. 9 Rivulet thickness at breakup point at different free stream velocities |
|
| 图 10 不同来流速度条件下破裂处溪流宽度分布 Fig. 10 Rivulet width at breakup point at different free stream velocities |
4.5 防冰溢流水表面换热特性
溢流水的流动形态会影响表面换热特性。图 11和12所示为不同来流速度条件下，三种不同溢流水流动形态表面蒸发和对流换热热流的分布情况。首先，在相同来流速度下，完全湿润表面(连续水膜表面)的蒸发热流总大于部分湿润表面(破裂处溪流干/湿表面)，这是由于部分湿润表面中的干表面无蒸发热流。其次，这两种状态下的蒸发热流随着来流速度的增加均呈现明显增大的变化趋势。这是由于来流速度较大时，表面传质系数相应增大，从而影响表面蒸发热流。
|
| 图 11 不同来流速度条件下表面蒸发热流量 Fig. 11 Evaporation heat fluxes on anti-icing surface at different free stream velocities |
|
| 图 12 不同来流速度条件下不同表面对流换热热流量 Fig. 12 Convection heat fluxes on anti-icing surface at different free stream velocities |
表面对流换热热流随来流速度增加的变化趋势与蒸发热流相同，而对于同一来流速度下三种不同溢流水形态条件下表面对流换热热流的区别并没有蒸发热流那么明显。从图 12还可以看出，相同速度下，完全干表面(无溢流水)的对流换热热流相对较大，部分湿润表面次之，完全湿润表面对流换热量最小。这说明，尽管表面溢流水形成的水膜/溪流厚度很小，仍会对对流换热表面的温度产生一定的影响，只是影响量有限。
5 结 论基于空气-水两层相互作用的质量、动量和能量守恒建立的防冰表面溢流水水膜流动换热及破裂的数学模型，研究分析来流速度对表面连续水膜厚度及主要散热项、破裂点溪流厚度与宽度的影响。同时，对三种不同溢流水流动形态下的表面蒸发与对流换热热流进行对比，分析溢流水对表面换热特性的影响。主要结论可以归纳为：
1) 撞击区域内，在撞击区域内的蒸发质量远小于撞击质量，因此溢流水质量会不断积累，使得连续水膜厚度沿表面的分布呈现不断增大的变化趋势；
2) 来流速度的增加使得表面撞击区域以及撞击水量增加，表现为连续水膜在表面的铺展面积及厚度均随之增加；
3) 来流速度对蒸发和对流换热热流的影响主要取决于换热系数，来流速度增大，表面换热剧烈，从而使得表面蒸发热流和对流换热热流增加；
4) 由于来流速度的增加，连续水膜在破裂处的厚度也会变大，使得破裂点溪流厚度和宽度增大；
5) 在撞击区域之外，翼型表面考虑溢流水的存在使得表面换热过程蒸发热流增加，而当溢流水流动形态由连续水膜破裂成溪流时，表面蒸发热流有所下降；对于表面对流换热而言，无溢流水流动的表面热流量略有增加，而完全湿润表面换热热流量最小。
致谢: 论文得到了上海交通大学机械与动力工程学院工程热物理研究所陈勇老师、雷桂林博士的支持和帮助，这里对他们的帮助表示感谢。| [1] | Melody J W. In-flight characterization of aircraft icing[D]. Illinois: University of Illinois at Urbana-Champaign, 2004. |
| [2] | Politovich M K. Aircraft icing[J].Encyclopedia of Atmospheric Sciences, 2003:68–75. |
| [3] |
Qiu X G, Han F H.
Aircraft anti-icing system[M]. Beijing: Aviation Industry Press, 1985 : 44 -58.
(in Chinese) 裘燮纲, 韩凤华. 飞机防冰系统[M]. 北京: 航空工业出版社, 1985 : 44 -58. |
| [4] | Addy H E. Ice accretions and icing effects for modern airfoils[R]. NASA TP 2000-210031. |
| [5] | Messinger B L. Equilibrium temperature of an unheated icing surface as a function of air speed[J].Journal of the Aeronautical Sciences, 1953, 20(1):29–42. |
| [6] | Mac A C. Numerical simulation of airfoil ice accretion[R]. AIAA 83-0112. |
| [7] |
Yi X, Zhu G L. Computation of glaze ice accretion on airfoil[J].Acta Aerodynamica Sinica, 2004, 22(4):490–493. (in Chinese) 易贤, 朱国林. 考虑传质传热效应的翼型积冰计算[J]. 空气动力学学报, 2004, 22(4) : 490–493. |
| [8] | Hedde T, Guffond D. ONERA three-dimensional icing model[J].AIAA Journal, 1995, 33(6):1038–1045. |
| [9] | Honsek R, Habashi W G. FENSAP-ICE: eulerian modeling of droplet impingement in the SLD regime of aircraft icing[R]. AIAA 2006-465. |
| [10] | Morency F, Tezok F, Paraschivoiu I. Anti-icing system simulation using CANICE[J].Journal of Aircraft, 1999, 36(3):999–1006. |
| [11] | Al-Khalil K M, Wright W B, Miller D R. Validation of NASA thermal ice protection computer codes: part 3-the validation of ANTICE[R]. AIAA 1997-0051. |
| [12] | Dong W, Zhu J, Zheng M, et al. Thermal analysis and testing of non-rotating cone leading edge hot air anti-icing system[J].Journal of Propulsion and Power, 2015, 31(1):896–902. |
| [13] | Dong W, Zhu J, Zhou Z X, et al. Heat transfer and temperature analysis of an aero-engine strut under icing condition[J].Journal of Aircraft, 2015, 52(1):216–225. |
| [14] | El-Genk M S, Saber H H. Minimum thickness of a flowing down liquid film on a vertical surface[J].International Journal of Heat and Mass Transfer, 2001, 44(15):2809–2825. |
| [15] | El-Genk M S, Saber H H. An investigation of the breakup of an evaporating liquid film, falling down a vertical, uniformly heated wall[J].Journal of Heat Transfer, 2002, 124(1):39–50. |
| [16] |
Zhou X Q, Chen P L. Performance of the water flow film along a vertical plate[J].Journal of Tongji University, 1995, 23(6):621–626. (in Chinese) 周孝清, 陈沛霖. 竖直平板上的水膜流动特征[J]. 同济大学学报:自然科学版, 1995, 23(6) : 621–626. |
| [17] |
Ye X M, Yan J G, Li C X, et al. Breakdown model of locally heated and subcooled films driving by interfacial shear and gravity[J].Journal of North China Electric Power University, 2010, 37(2):68–73. (in Chinese) 叶学明, 闫俊刚, 李春曦, 等. 界面切应力和重力驱动下受热过冷液膜的破断模型研究[J]. 华北电力大学学报, 2010, 37(2) : 68–73. |
| [18] | Al-Khalil K M. Numerical simulation of an aircraft anti-icing system incorporating a rivulet model for the runback water[D]. Toledo: The University of Toledo, 1991. |
| [19] | Silva G A L, Silvares O, Zerbini G J, et al. Numerical simulation of airfoil thermal anti-ice operation, part 1: mathematical modelling[J].Journal of Aircraft, 2007, 44(2):627–633. |
| [20] | Silva G A L, Silvares O, Zerbini G J, et al. Numerical simulation of airfoil thermal anti-ice operation, part 2: implementation and results[J].Journal of Aircraft, 2007, 44(2):634–641. |
| [21] | Silva G A L, Silvares O, Zerbini G J, et al. Differential boundary-layer ananlysis and runback water flow model applied to flow around airfoils with thermal anti-ice[R]. AIAA 2009-3967. |
| [22] | Thundat T, Oden P I, Warmack R J. Microcantilever sensors[J].Microscale Thermophysical Engineering, 1997, 1(3):185–199. |
| [23] | Dong W, Zheng M, Zhu J, et al. Calculation and analysis of runback water flow on anti-icing surface[R]. AIAA 2015-0538. |
| [24] | Zhang K, Blake J, Rothmayer A, et al. An experimental investigation on wind-driven rivulet/film flows over a NACA0012 airfoil by using digital image projection technique[R]. AIAA 2014-0741. |
| [25] | Zheng M, Dong W, Lei G L, et al. Study of the flow and heat transfer of water film on hot air anti-icing airfoil surface[R]. IHTC 15-8671. |