过冷大水滴变形与破碎的影响因素

2019-01-18李维浩易贤李伟斌王应宇

航空学报 2018年12期

李维浩，易贤，李伟斌，王应宇

中国空气动力研究与发展中心空气动力学国家重点实验室，绵阳 621000

飞机结冰是威胁飞行安全的主要因素之一[1-2]。长期以来，研究者一直认为飞机结冰是由平均容积直径小于50 μm的过冷小水滴所致，因此相关结冰研究也主要针对过冷小水滴开展[3]。1994年ATR-72空难之后[4]，过冷大水滴(Supercooled Large Droplet, SLD)结冰问题开始引起关注，关于SLD结冰的研究热潮一直持续至今。研究表明，流场中的SLD在与物面相向运动的过程中，存在明显的变形，并在一定条件下会发生破碎，这种水滴的变形与破碎现象会对结冰产生直接影响。因此，在进行过冷小水滴结冰研究时关于水滴在运动过程中保持球体、不会变形和破碎的假设，并不适用于SLD结冰研究[5]。

针对SLD的变形与破碎问题，国内外做了大量的实验研究[6-8]，并提出了多个描述水滴变形与破碎的模型，这些实验都是在特定流场中进行的，缺乏普适性，所以将这些经验模型用于模拟水滴在机翼周围复杂流场的变形与破碎缺乏准确性[9]。针对这一不足，2005年Tan和Papadakis[10]提出用泰勒类比破碎(Taylor Analogy Breakup, TAB)模型[11]来描述水滴的变形与破碎过程，大大提高了SLD结冰数值计算的精度。虽然TAB模型只能模拟一种振荡模式，而水滴破碎通常由多种振荡模式共同作用导致，但目前还没有可靠的破碎模型来替代[12]，因此TAB模型依然被广泛应用于SLD结冰计算中。目前，国内外关于SLD变形与破碎数值计算方面的研究主要集中在变形与破碎对水滴收集系数和冰形的影响上[13-15]，缺少关于变形与破碎影响因素的系统研究。鉴于此，本文采用数值计算的手段，基于水滴运动计算的拉格朗日法和TAB模型，系统研究了水滴直径、流场速度和绕流物体尺寸对水滴破碎与变形及其效应的影响规律，在此基础上，提出了一个表征水滴变形与破碎及相应效应的综合影响参数。

1 数值计算方法

1.1 空气流场

空气流场通过求解低速黏性流场的时均Navier-Stokes方程获得，采用SIMPLE系列算法求解[16]，控制方程的通用形式为

(1)

式中：t为时间；φ为输运变量；ρa为空气密度；ua为空气速度；Γφ为扩散系数；qφ为源项。φ、Γφ和qφ取不同的值，可代表流场的连续性方程、动量方程以及湍动能、湍流耗散率等标量的输运方程。式(1)中各项的物理意义、具体表达式和方程的求解方法，可详见文献[16]。

1.2 水滴轨迹

获得空气流场之后，本文采用拉格朗日法[17-18]进行水滴轨迹计算。根据牛顿第二定律，水滴运动方程可以写为

(2)

式中：xd为水滴位移；ρd为水滴密度；Vd为水滴体积；g为重力加速度；Ad为水滴的迎风面积；CD为阻力系数；ua为当地流场速度；ud为水滴速度。引入相对雷诺数：

(3)

式中：Deq为水滴当量直径；μa为空气黏性系数。则式(2)可改写为

(4)

以水滴的初始位置和流场的速度分布作为定解条件，水滴运动方程的求解属于一阶常微分方程的初值问题，本文采用一阶欧拉法进行数值求解。

1.3 水滴变形与破碎

采用TAB模型[10-11]模拟水滴变形与破碎。TAB模型是假定水滴振荡运动类似于弹簧-质量-阻尼系统的受迫运动(如图1所示)。

图1 泰勒类比破碎模型Fig.1 Taylor analogy breakup model

描述弹簧-质量-阻尼系统的受迫运动的数学方程可以写为

(5)

式中：位移x是位于初始球体赤道处的点的位移；外力F是施加在水滴上的气动力；恢复系数k是由水滴表面张力引起的；阻尼系数f为水滴黏度。引入系数Cb、CF、Ck和Cd，且定义

Y=x/(Cbr)

(6)

(7)

(8)

(9)

式中：Y为水滴的无量纲变形量；r为水滴半径。则式(5)可写为

(10)

其中：

(11)

式(10)采用如下差分格式进行离散求解：

(12)

式中：

(13)

当位移超过初始球体半径的一半，即Y>1时，认为水滴破碎发生[10]，破碎后产生的子水滴直径与母水滴直径的关系表示为

(14)

其中：D为母水滴直径；D32为子水滴粒径分布的Sauter平均直径[19]，定义为对应尺寸水滴的体积与表面积之比等于所有水滴的总体积与总表面积之比；K为水滴的变形和振动能量与基频总能的比值，计算中K=10/3。子水滴直径服从χ2分布[11]，即

(15)

式中：Γ(x)为伽马函数；n为该分布的自由度，可取任意大于0的整数，子水滴直径为自变量。由Sauter平均直径的定义可求得n与D32的关系为

n=D32+2

(16)

子水滴的数量可由质量守恒求得，且认为子水滴是无变形、无振荡的，并且具有与母水滴相同的速度。

1.4 水滴阻力

常见的飞机飞行速度范围以及云层中经常出现的水滴尺寸范围，都超出了使用Stokes公式计算水滴黏性阻力的理论允许范围，因此不能使用Stokes阻力计算公式，而使用低雷诺数下球体的严格阻力系数为[20]

(17)

水滴变形会引起阻力增加，阻力系数进行以下修正[10]：

CD=CD,S(1+2.632Y)

(18)

1.5 数值计算流程

数值计算的流程可以总结为：① 输入水滴初始的位置、直径等信息，读取之前已经获得的流场信息；② 判断水滴所处网格位置，得到当地流场速度，使用一阶欧拉法迭代求解水滴轨迹，并且每迭代一步都判断水滴是否撞击到物体表面或者飞出流场；③ 使用TAB模型计算水滴的无量纲变形量Y，当Y≥1时发生破碎，计算出破碎产生子水滴的平均直径，采用相同方法计算子水滴的轨迹，并计算子水滴的粒径分布和位置分布。

计算水滴轨迹及变形与破碎的具体流程如图2 所示。

图2 数值计算流程Fig.2 Workflow of numerical calculation

2 SLD变形与破碎特性计算分析

2.1 计算条件

SLD在机翼附近变形与破碎和水滴直径以及水滴和空气的相对速度有关，而相对速度主要决定于初始来流速度和物体的尺寸外形。对于NACA四位数字翼型，前缘半径rL≈1.10c2/b，c为翼型厚度，b为弦长。因此，本文选取弦长均为1 m的NACA0008、NACA0012、NACA0015 3种厚度不同的对称翼型，来流速度为90、120、150 m/s，水滴直径为50、100、150、200、300 μm作为计算条件。

2.2 典型SLD流场区域划分

图3是直径300 μm的水滴在速度150 m/s的流场下，在NACA0015翼型附近发生破碎前后的典型轨迹，图3(b)是破碎点附近(图3(a)的虚线圈)的放大图。水滴破碎的原因是受到空气剪切作用，当水滴表面张力小于其受到的空气剪切力时，发生破碎。当水滴离翼型较远时，运动速度与流场一致，与空气之间没有相互剪切力，而当水滴与翼型距离减小时，与空气的速度差增大，在剪切力作用下，水滴先是变形，进而破碎成若干子水滴继续前进，并碰撞在翼型的不同区域和位置。

从上面分析可知，包含SLD的典型流场可以划分成以下几个区域：① 水滴未撞击区域，② 水滴撞击但未破碎区域；③ 水滴待破碎区域；④ 水滴破碎区域。具体分布如图4所示，图中蓝色区域的水滴不会撞击在翼型表面，灰色区域水滴在撞击翼型表面之前仅发生变形不会发生破碎，绿色区域水滴即将发生破碎，褐色区域是水滴破碎后存在子水滴的区域。在水滴待破碎区域，水滴变形量随时间和空间推移而增加，越靠近翼型表面变形量越大，在两区域交界处发生破碎。

图3 水滴破碎前后的典型轨迹Fig.3 Typical trajectories of droplet before and after breakup

图4 SLD结冰流场典型区域划分Fig.4 Typical zoning of SLD icing flow field

2.3 SLD破碎区域面积

表1～表3分别给出了3种翼型流场在不同速度和水滴直径条件下的破碎区域面积。可以看到，在所有计算状态中，只有几个速度高、水滴直径大的状态存在水滴破碎区域。水滴破碎区域面积受到流场速度、水滴直径以及翼型前缘半径的共同影响，速度越高、直径越大、前缘半径越大，水滴破碎区域面积越大。

表1 NACA0008翼型水滴破碎面积Table 1 Breakup area of NACA0008 airfoil mm2

表2 NACA0012翼型水滴破碎区域面积Table 2 Breakup area of NACA0012 airfoil mm2

表3 NACA0015翼型水滴破碎区域面积Table 3 Breakup area of NACA0015 airfoil mm2

2.4 水滴收集特性

图5是有/无变形与破碎模型极限轨迹的对比(速度v=150 m/s，水滴直径d=300 μm，NACA0015翼型)。图中可以看出添加变形与破碎模型后撞击极限变小，这是由于水滴变形引起阻力系数增大，水滴的跟随性增强，轨迹的偏折程度会变大。图6给出的是有/无变形与破碎模型的水滴收集系数β对比，横坐标是碰撞点离驻点的物面距离s与弦长b的比值。其中图6(a)显示的是直径50 μm的SLD在速度150 m/s的流场中NACA0015翼型上的水滴收集系数，虽然该流场没有水滴破碎，但由于SLD在运动中会发生变形，变形导致水滴阻力系数增大，水滴运动轨迹与空气流线一致性更高，从而导致局部水收集系数和撞击极限都变小。图6(b)显示的是在相同流场下直径300 μm的SLD的水滴收集系数，该流场包含有较大的水滴破碎区域，考虑变形与破碎之后，局部水收集系数和撞击极限变化程度比图6(a)状态更大。

从图6可以看出，不同的结冰条件，变形与破碎对水滴收集特性影响不一样。图7(a)给出的是对于NACA0012翼型，不同速度条件下变形与破碎对水滴总收集系数影响的曲线；图7(b)显示的是速度为120 m/s时，不同翼型条件下变形与破碎对水滴总收集系数影响的曲线，图中的纵坐标表示的是考虑变形与破碎之后水滴总收集系数变化的比例。可以看出水滴直径、流场速度、机翼前缘半径等与变形和破碎影响水滴总收集系数的程度呈正相关关系。速度为90 m/s时，影响程度较小，总收集系数变化不超过4%，速度增加至150 m/s之后，最大变化量可达到10%以上。

图5 极限轨迹对比Fig.5 Comparison of limiting trajectories

图6 局部水滴收集系数对比Fig.6 Comparison of local droplet collection coefficients

2.5 子水滴分布

图8给出的是破碎后产生的子水滴的粒径分布，母水滴是在同一流场状态下(v=150 m/s,d=150 μm, NACA0012)，在远场的不同纵坐标处释放的，纵坐标y分别为0和0.05。可以看出，水滴破碎后，会产生大量的子水滴，其分布类似正态分布，在同一流场状态下，在不同位置释放的水滴破碎后产生子水滴，其粒径分布十分近似。

图7 总收集系数变化的比例Fig.7 Changing ratio of total collection coefficients

图8 不同位置释放水滴的子水滴粒径分布Fig.8 Diameter distribution of sub-droplet dropping at different positions releasing droplet

图9给出的是3种不同状态下，在远场的y=0处释放的直径为300 μm的水滴簇破碎后子水滴的粒径分布。可见，当水滴直径相同时，在不同翼型和不同流场速度条件下得到的子水滴粒径分布是十分接近的。所以，水滴发生破碎之后子水滴的粒径分布主要与SLD的初始直径有关，流场速度、翼型、释放位置等因素影响十分小。

图10是SLD破碎后产生子水滴撞击物面的位置分布，水滴直径为300 μm，流场速度为150 m/s，释放点位于NACA0015翼型前某位置。子水滴的位置用碰撞点距离未考虑破碎时的撞击点(已考虑变形)的物面距离s来表示。水滴破碎后，子水滴在空气动力作用下，撞击位置比不考虑破碎的撞击位置要远离驻点，并散落在一定区域，图10 中水滴散落的范围约为整个撞击区域的0.63%。可见水滴破碎会改变撞击在机翼表面的水滴大小及分布，并导致物面局部液态水的再分配，进而影响液态水在机翼表面的溢流、传热等物理过程。

图9 不同流场的子水滴粒径分布Fig.9 Diameter distribution of sub-droplet in different flow fields

图10 SLD破碎产生子水滴撞击位置分布Fig.10 Impingement location distribution of sub-droplet after SLD breakup

3 水滴变形与破碎的综合影响参数

通过前文的分析可知，水滴变形与破碎对撞击特性的影响程度随流速、直径以及翼型前缘半径增大而增大。为了综合考虑各个影响因素，定义水滴变形与破碎的综合影响参数X为

X=v2drL

(19)

式中：rL为机翼前缘半径。以X为横坐标，子水滴在翼型表面的覆盖范围为纵坐标，拟合成直线，如图11所示。线性回归的相关系数R为0.910 3，说明X与子水滴覆盖范围有较强的线性相关关系。拟合直线的横截距为0.040 4，未发生破碎状态的X都小于或稍微大于0.040 4，与拟合直线得到结果吻合。

为了进一步验证X作为表征SLD变形与破碎效应参数的适用性，对X与水滴收集系数的关系进行了研究。图12为X与总收集系数差值的关系图。图中纵坐标为考虑变形与破碎和不考虑变形与破碎对应的总收集系数的差值。可以看出总收集系数的变化量与X呈现较强的线性相关性，初步说明X是适合衡量SLD变形与破碎效应的参数。当X=0.040 4时，考虑变形与破碎和不考虑变形与破碎的总收集系数差值为0.042，比对应状态总收集系数小一个数量级。水滴直径为200 μm、流场速度为90 m/s、翼型为NACA0015的状态下X=0.040 1，最为接近拟合直线求得的水滴破碎阈值，该状态下总收集系数变化率约为4.56%。说明X<0.040 4时，水滴变形与破碎对SLD的撞击特性影响较小。