0 引言

飞机及发动机结冰是危及飞行安全的严重问题，与结冰有关的飞行事故每年都有发生。当飞机穿越含有过冷水滴的云层时，过冷水滴撞击到飞机机翼以及发动机的进气部件表面，极易发生冻结。飞机及发动机部件结冰会改变部件的气动性能，对飞机造成极大的安全隐患。融化的冰如果被发动机吸入，会导致发动机叶片的机械损伤，引起发动机熄火停车，造成致命后果。

飞机结冰和防冰过程涉及了复杂的多相流动换热过程，其撞击水滴冻结过程的数学模型至今仍不完善。有关飞机结冰和防冰过程中大水滴撞击、水膜的流动与换热、结冰引起的表面粗糙度对流动换热的影响等研究工作在近年来得到了研究人员的广泛关注^[1-6]。过冷水滴撞击特性是飞机结冰和防冰分析的基础。通常，直径小于50 μm的过冷水滴不考虑其在运动及撞击过程中的变形，在撞击机翼表面后不发生飞溅脱离现象。对于小水滴在机翼表面的撞击可以较为准确地采用拉格朗日方法和欧拉方法做出预测^[7-8]。而直径大于50 μm的过冷水滴通常称为过冷大水滴(Supercooled Large Droplets, SLD)。大水滴在撞击飞机迎风表面后，会发生飞溅、反弹、二次撞击等多种情况。大水滴的破碎和飞溅会使得撞击水滴的直径和部件表面的水收集系数发生变化，从而影响部件表面生成的冰型。对于飞机结冰分析中大水滴的撞击特性，目前主要还是采用实验数据总结的经验模型进行计算分析^[9-10]。由于SLD撞击特性的复杂性，目前还没有一个经验模型在各种结冰条件下都可以给出满意的撞击特性预测结果。因而对SLD撞击过程的机理研究还是非常迫切且必要的。

高速摄像技术的出现为大液滴撞击表面的实验研究提供了强有力的支持。近年来许多研究人员利用高速摄像技术和液滴生成装置研究了液滴撞击形态的变化。Thoroddsen等^[11]、Yarin等^[12]通过实验方法拍摄了液滴撞击壁面后的铺展、快速溅射、冠状溅射、收缩溅射、部分回弹、完全反弹等多种结果。Rioboo等^[13]着重研究了壁面浸润性及表面张力对液滴撞击壁面铺展的影响。Pan等^[14]研究了不同属性的液滴撞击不同粗糙度的固体表面的飞溅情况，通过高速摄像装置拍摄了大韦伯数(We数)液滴撞击壁面发生放射状飞溅和日冕皇冠状飞溅，发现当韦伯数足够大时即使在光滑壁面上也能发生液滴飞溅。Tsai等^[15]实验研究了微纳米超疏水表面对液滴撞击的影响。液滴撞击的实验研究也具有很大的局限性，通过高速摄像只能得到不同时刻的液滴撞击的形态图片，液滴撞击过程中的速度分布、压力分布、换热情况等均不能测出。同时由于高速液滴撞击实验条件实现困难，目前大量进行的液滴撞击试验的初始速度一般小于40m/s。

理论研究方面，Naber和Reitz^[16]在类比射流喷射的基础上考虑了液滴破裂、液滴碰撞、液滴融合以及空气涡流，提出了著名的N-R模型，认为韦伯数决定了液滴碰壁后形态变化。Rioboo等^[17]定义了液滴撞击壁面后的几种不同状态：沉积、回弹、部分回弹、迅速飞溅、冠状飞溅、收缩破裂。Cossali等^[18-19]、Pan和Law^[20]、Pan等^[21]通过大量实验及数据分析确定了液滴撞击湿表面不同形态结果的临界韦伯数。数值模拟方面，Gunjal等^[22]使用VOF (Volume of Fluid)方法模拟了液滴撞击壁面的过程，并通过实验验证了VOF模型模拟液滴撞击壁面结果的可行性。Kim等^[23]使用VOF模型数值模拟了不同速度、流变参数、表面张力条件下的液滴撞击壁面，结果表明液滴的铺展形态取决于雷诺数和韦伯数，液滴的收缩形态决定于毛细数、Bingham-Capillary数，并预测了液滴撞击壁面的最大铺展半径。Bussmann等^[24]使用三维模型研究了液滴撞击不对称表面过程中的动态壁面接触角的影响，在模拟液滴撞击倾斜壁面时使用了实验测得的液滴铺展动态壁面接触角作为壁面边界条件，与实验结果较为符合，并在此基础上提出了一种新的给定动态接触角的方法。梁超^[25]使用三维模型和VOF方法研究了液滴低速撞击等温壁面过程中的壁面接触角、液滴初始速度、表面张力系数、液滴动力粘度等的影响及液滴撞击热壁面的换热特性，此外还使用二维轴对称模型模拟了不同液滴初始速度、表面张力系数、液滴动力粘度、初始液膜厚度等条件下的液滴撞击液膜过程。贾小娟^[26]使用三维模型完成了双液滴垂直及斜向撞击液膜的数值研究。Burtnett^[27]使用二维轴对称模型模拟了50 μm液滴撞击微结构表面，对于干表面和渗透表面得到了不同的液滴撞击结果。Tan等^[28]采用大水滴破碎、飞溅、反弹模型进行了结冰翼型表面的水滴撞击特性计算，并与试验结果进行了对比验证。

针对飞机结冰所涉及的SLD直径范围内的水滴高速撞击壁面和液膜的研究较少。本文借助VOF方法开展了大水滴撞击壁面过程的动态数值模拟研究，分析了大水滴撞击过程中形态变化、飞溅现象以及撞击特性，数值模拟结果加深了对飞机表面过冷大水滴撞击过程物理机理的理解。

1 计算模型 1.1 计算模型

结冰云层中SLD直径分布大致在50 μm到500 μm的范围内。大水滴撞击飞机表面后会发生飞溅、反弹、二次撞击、铺展-收缩-振荡等多种情况。当飞机表面聚集的水滴足够多时，水滴会在飞机表面形成一层液膜，随后还会发生大水滴撞击薄液膜的情况。

数值模拟中采用二维轴对称模型，计算区域为矩形区域，中间为对称轴，底面为壁面边界条件。对于大水滴高速撞击壁面的数值模拟，为了能够细致地计算出水滴的形态变化，要求计算区域的网格足够密集，使用二维轴对称模型能有效地减小网格数量，缩短计算时间。水滴撞击表面的几何模型如图 1所示。

在数值模拟中采用的二维轴对称模型，考虑了一些流体流动过程中的三维效应，能在一定程度上代替三维网格。在水滴低速撞击的很多文献中均使用二维轴对称模型模拟三维的液滴撞击壁面和液滴撞击液膜。相对于三维模型，二维轴对称模型认为流动过程中只存在沿径向的梯度，沿周向均匀分布，因此二维轴对称模型不能够完整捕捉液滴撞击壁面过程中出现的沿周向不均匀、不连续的现象。但使用二维轴对称网格模拟液滴撞击壁面，计算得到的液滴的铺展系数、液膜厚度等数据可信度高，同时能够捕捉到液滴铺展和收缩过程中的二次液滴飞溅、液膜断裂等现象。

1.2 VOF模型

大水滴高速撞击壁面的数值模拟，是一种水-空气的两相流动分析。数值追踪两相流界面的方法包括Level-Set方法、粒子标记方法、VOF方法等。其中VOF方法相对于其他的界面追踪方法模型简单，追踪界面精确，能够考虑界面融合和分离现象^[29]。

Hirt和Nichols^[30]提出的VOF方法，计算中通过引入流体体积分数α来标记界面两边不同相的流体，动态追踪不同相的区域，采用界面重构技术来确定相界面。若假设计算域中共有n相流体，对于任意计算网格的第k相流体来说，其体积分数可能存在如下三种情况，即：

(1)

如果α_k等于0，说明在这个网格内不含有第k相流体；如果等于1则说明在这个网格内只存在第k相流体；而如果α_k处于0和1之间，表明这个网格内部同时存在多相流体。由于VOF模型假定多相流体不能混合，可以认为该网格处于多相流体界面或流体界面附近。根据体积分数的定义，在同一个网格内部，各相流体的体积分数和应该为1，即：

(2)

VOF方法在每一个网格内采用体平均来计算流动方程并通过流体体积分数捕捉相的界面，图 2为VOF界面追踪示意图。

在本文讨论范围内，所涉及的流体相只包含空气相和水滴相两相，在每个计算网格中，第k相(k=1，2)的连续性方程可表示如下：

(3)

在VOF模型中，认为在相同位置处的不同相流体具有同样的速度，因此对于空气相和水滴相采用同一套动量方程进行求解：

(4)

其中，P为压力，F_SF为由表面张力引起的动量源项。式(4)中的物性参数ρ和μ采用体积加权平均法获得：

(5)

(6)

1.3 表面张力模型

表面张力采用连续表面张力CSF(Continue Surface Force)模型，CSF模型将VOF计算中附加的表面张力处理为动量方程的源项。相界面两侧的内外压差等于表面曲率和与表面张力系数之积：

(7)

式中，p₁和p₂分别为界面两侧的流体压力，R₁和R₂为两个主曲率半径。

在水滴撞击的过程中表面张力以及壁面浸润性起着重要作用，正是由于表面张力的存在，水滴在撞击壁面达到最大铺展直径后会发生收缩反弹、收缩振荡或者液膜在表面张力作用下拉裂等现象。由于水滴撞击过程中的动态接触角难以测得，数值模拟中使用静态接触角代替动态接触角。

水滴撞击壁面是一个瞬态不可压缩过程，数值模拟采用PISO算法。时间步长和空间步长的选取对于水滴撞击特性的数值模拟非常重要，合适的空间步长和时间步长能够用较少的时间得到精确的结果。空间步长的选取通过网格无关性验证来确定。时间步长通过克朗数确定：

(8)

其中：dt为时间步长；dx为空间步长；Cr可以理解为内扩散相在单位时间内前进的距离与空间步长的比值，这里Cr取为0.25确定初始时间步长。

2 计算方法模型验证

结冰云层过冷大水滴直径为几十到几百微米，水滴撞击速度很大，要保证计算过程中合适的Cr，需要的网格步长很小。如果选择的网格较小，虽然能够更精确、细致地计算出大水滴撞击壁面后的形态变化，但是整个网格的数量将急剧增大，同时需要设定较小的时间步长，这样整个模拟的时间将大大增长；如果选择的网格步长过大，则不能精确地计算水滴撞击壁面后的形态变化，计算结果偏差较大，相的界面不清晰。

大水滴高速撞击的试验数据很少，因此这里以水滴速度较低的试验数据进行了计算方法的验证。选取试验条件水滴直径2.6mm，水滴初始速度0.5m/s，壁面接触角为160°作为对比算例。图 3为不同网格尺度下不同时刻水滴撞击壁面后的相图。图 4为撞击区域网格尺度分别为0.02mm、0.01mm、0.0075mm条件下水滴在壁面上的铺展系数(d/D)随无量纲时间(t · v/D)的变化曲线，其中d为水滴的铺展直径，D为水滴的初始直径。通过对比发现，相同时刻0.01mm与0.0075mm网格的水相形态和水滴铺展系数基本相同，而与0.02mm网格步长下的水滴形态差异相对较大。在本算例中，0.01mm是一个合适的网格步长，既能保证计算结果的精度，也能保证计算的速度，后面的模拟计算以此网格尺度所对应的无量纲网格尺度作为网格划分的参考依据。

图 5为数值模拟得到的相同时刻水滴撞击壁面后的相图与文献[31]相同时刻实验拍摄图片的对比图。通过对比发现，数值模拟得到的结果与文献中实验得到的结果基本相同，验证了本文计算方法的可靠性。

3 计算结果及分析

数值模拟中大水滴直径D选取500 μm、300 μm、100 μm和50 μm，水滴初始速度v₀选取30m/s、40m/s、50m/s、60m/s、70m/s。以铝合金材料表面模拟飞机部件表面，壁面接触角为75°。

3.1 大水滴直径对撞击特性的影响

研究在高速情况下不同直径的水滴撞击壁面特性时，保持水滴的初始速度为60m/s。观察不同直径下水滴撞击壁面后的形态变化，分析水滴撞击壁面后的铺展系数及铺展时间。

以250网格/直径的网格模拟500 μm水滴撞击壁面的动态过程，确定了水滴的最大铺展半径大约为2mm，因此取模型计算区域为3mm×3mm的对称区域。由于使用该网格步长/直径模拟时，铺展液膜最薄处仅有2层网格的厚度。为了更精确地反映高速水滴撞击壁面过程，对计算区域近壁面处进行了局部加密处理，使用400网格/直径并采用渐密网格。贴近壁面的第一层网格高为0.1 μm，网格高度方向尺度增长率为1.01。

图 6为水滴以不同直径撞击壁面的铺展半径随时间变化曲线。图 7为水滴以不同直径撞击壁面的铺展系数随无量纲时间变化曲线。图 8为水滴直径与最大铺展半径及最大铺展系数关系曲线。图 9为水滴直径与达到最大铺展对应时间的关系曲线。分析图 6~图 9可以看出：大水滴高速撞击壁面时，随着大水滴直径的增大，同一时刻水滴的铺展半径和铺展系数都增大，水滴达到最大铺展对应的时间增大，对应的无量纲时间增大。

数值计算分析发现：大水滴以高速撞击壁面时，水滴在壁面上的铺展系数较大，中心液膜的厚度极薄，一般小于10 μm；水滴直径小于300 μm时，液膜的收缩阶段会发生液膜破裂；随着水滴直径的增大，在最大铺展时刻对应的液膜厚度增大，水滴铺展阶段伴随有部分水滴飞溅，飞溅的水滴直径小于50 μm。当水滴初始直径为500 μm时，水滴在壁面上快速铺展，在铺展过程中有发生飞溅现象。由于撞击的初始速度高，水滴铺展系数大，达到最大铺展系数时刻液膜的厚度只有7 μm。在表面张力的作用下，整个薄液膜逐渐收缩，但在收缩阶段并没有发生液膜拉断的情况。随着水滴初始直径的缩小，水滴撞击铺展后的液膜厚度逐渐减小。在直径为300 μm条件下，铺展过程伴随飞溅现象，由于达到最大铺展系数时刻液膜较薄，在表面张力的作用下，收缩阶段发生中心液膜拉断现象；由于被拉断的中心液膜面积较小，表面张力的回复作用不明显，水滴形态变化所储存的能量较少，表面张力在液膜收缩为中心小水滴过程中做功较少，中心液膜收缩形成的中心水滴具备的动能较低，中心水滴没有脱离壁面发生反弹，而是在壁面上不断振荡。在水滴直径为100 μm和50 μm条件下，由于铺展的液膜厚度太薄，表面张力的作用更加剧烈，整个液膜被拉断为多个部分，同时在表面张力的强烈收缩作用下，破裂液膜收缩形成的水滴伴随有飞溅和反弹脱离壁面的情况。由上述分析可知，在水滴撞击壁面达到最大铺展后，水滴的形态变化达到最大，同时水滴的表面张力达到最大，薄液膜开始收缩。如果水滴铺展形成的液膜厚度达到临界值，在较大的表面张力作用下，液膜会被拉断形成外围的液环和中心的液膜；在极端情况下，整个液膜会形成多个外层液环甚至整个液膜被完全拉断为不规则的多个部分(完全破碎)。表 1给出了不同直径大水滴撞击壁面的液膜厚度及飞溅、破裂情况。

3.2 大水滴速度对撞击特性的影响

研究大水滴以不同速度撞击壁面的特性时，保持水滴的直径为300 μm，大水滴以初始速度分别为30m/s、40m/s、50m/s、60m/s、70m/s。观察不同初始速度的大水滴撞击壁面后的形态变化，分析不同时刻水滴的铺展半径及飞溅现象。

图 10为不同速度下大水滴撞击壁面后的局部放大相图。图 11为不同速度的大水滴撞击壁面后的铺展半径随时间的变化曲线。图 12为不同速度的大水滴撞击壁面后的铺展系数随无量纲时间的变化曲线。图 13为大水滴撞击速度与最大铺展半径及最大铺展系数关系曲线，图 14为大水滴撞击速度与达到最大铺展时间和达到最大铺展对应无量纲时间的关系曲线。从图 11~图 14可以看出：大水滴高速撞击壁面时，随着大水滴初始撞击速度的增大，同一时刻水滴的铺展半径和铺展系数同时增大，水滴的最大铺展半径也增大，水滴达到最大铺展对应的时间缩小，尽管水滴达到最大铺展对应的无量纲时间相差不大，但达到最大铺展时刻对应的液膜厚度减小。

表 2给出了不同速度的大水滴撞击壁面后的液膜厚度及飞溅情况。初始速度为30m/s的水滴具备的能量相对较低，撞击壁面形成的液膜铺展半径较小，在液膜收缩阶段不发生液膜断裂；初始速度为40m/s和50m/s的水滴具备的能量较高，液膜的铺展半径更大同时更不稳定，在液膜的收缩阶段，边缘液环与中心薄液膜连接处断裂，初始速度为50m/s的初始水滴形成的液膜更是断裂为三部分；初始速度为60m/s和70m/s的水滴能量更高，在薄液膜的中心部位发生了液膜断裂情况，中心液膜断裂并聚集为一个小水滴不断振荡。这是因为：在表面张力的作用下，水滴撞击壁面的收缩阶段，水滴形成的液膜有向液膜中心收拢的趋势。由于液膜边缘处表面张力的作用最为强烈，液膜边缘部分的水相率先向内运动，而液膜中心部位的水相运动趋势不明显，中心部位的液膜厚度基本不发生变化。边缘液膜运动的结果是液膜的铺展半径变小，同时在液膜边缘处水相堆积形成凸状液环。

4 结论

使用VOF方法数值模拟了直径50~500 μm范围内的大水滴以高速撞击壁面的撞击特性。通过改变大水滴的直径和大水滴的初始速度，计算分析了不同条件下大水滴形态变化、撞击结果以及铺展系数的变化。

1) 直径50~500 μm的大水滴以高速撞击壁面时，水滴在壁面上的最大铺展系数较大，达到最大铺展后的液膜厚度界于1~10 μm之间。由于液膜的厚度很薄，在表面张力的作用下，可能发生液膜被拉断的行为。

2) 随着水滴直径以及撞击速度的增大，水滴的铺展速度、最大铺展半径、最大铺展系数均增大。水滴直径增大时，水滴达到最大铺展所用的时间变长，所用的无量纲时间呈微弱增大趋势；同一直径水滴初始速度增大时，水滴达到最大铺展所用的时间减少，所用的无量纲时间变化相对较小。

3) 低速情况下水滴撞击光滑壁面不易发生飞溅现象，但在水滴高速撞击光滑壁面的大韦伯数情况下，数值模拟发现了明显的飞溅现象。

致谢: 数值计算分析过程中得到了上海交通大学机械与动力工程学院工程热物理研究所陈勇老师的支持和帮助，这里对其帮助表示感谢。