二维机翼混合相结冰数值模拟

卜雪琴，李皓，黄平，林贵平

北京航空航天大学 航空科学与工程学院，北京 100083

从现代飞机飞行开始，飞机结冰问题就逐渐引起人们的重视。早期的研究主要关注过冷水滴结冰，随着研究地不断深入，飞机结冰研究也逐渐由过冷水滴结冰扩展到冰晶结冰、冰晶和过冷水滴同时存在的混合相结冰[1-3]。Mason等[4]总结发现涡扇发动机的动力损失与冰晶密切相关。冰晶结冰的机理尚不明确，目前普遍认为冰晶引发结冰需要液态水的存在[5]。尽管纯冰晶条件下，冰晶撞击到机翼、尾翼等冷表面会弹开而无法结冰，但撞击到传感器等热表面或在发动机内部热环境下运动时，冰晶会(部分)融化成液态水，液态水使得冰晶黏附冻结，同时冰晶融化吸热使得热表面、热环境温度下降。传感器结冰会导致传感器精度下降。发动机内部核心部件冰晶结冰严重者可能引发发动机动力下降、喘振、熄火等危险，同时冰层脱落可能造成发动机的机械损伤。研究表明，冰晶结冰时经常伴随着总温传感器异常[4]。此外，在同时含有冰晶及过冷水滴的混合相条件下，由于液态水的存在，结冰不仅可能发生在机翼、尾翼等冷表面，造成飞机气动外形的改变，降低升力、增大阻力、影响飞机稳定性，而且可能发生在传感器、发动机内部压气机等热表面，使传感器指示异常，对发动机构成威胁。由于冰晶的存在，冰晶结冰传热传质模型和过冷水滴结冰模型有显著区别。另外，由于冰晶形状各异，可能是扁平盘状、柱状等非球形，其运动以及在热环境中运动的融化相变也是一个复杂过程。

冰晶/混合相结冰的试验研究需要先进的冰风洞作为支持，存在的难点包括冰晶的制备与运输、真实结冰气象条件的模拟等，具有成本高、周期长的特点[6-8]。而采用数值模拟能大大降低人力、物力等资源消耗，节约时间成本，因此也是冰晶/混合相结冰研究的重要手段。目前国外针对冰晶/混合相结冰数值模拟问题已经展开广泛研究，在已有结冰软件基础上进行了改进和补充，使其能够应用到冰晶结冰问题中来，已经取得了不错的进展[9-12]。而国内研究相对较少，Zhang等[13]使用FLUENT软件研究了混合相结冰问题。袁庆浩等[14]对航空发动机内部冰晶结冰文献进行了调研和总结，综述了冰晶结冰和过程水结冰的区别、冰晶结冰试验和计算方法等。姜飞飞等[15]计算分析了在涡扇发动机内涵道内运动过程中冰晶粒子的传热传质情况，得到了冰晶粒子的温度、半径、液态水质量分数的变化结果。

本文暂不考虑冰晶的非球形特性，认为是球形，旨在建立混合相条件下的结冰热力学模型和数值计算方法，并开展二维翼型混合相条件下的结冰模拟。基于Messinger模型进行扩展，考虑冰晶黏附效应，建立了混合相的结冰热力学模型。首先采用FLUENT软件计算空气流场，然后基于单向耦合假设求解欧拉法下的冰晶运动控制方程，获得冰晶撞击特性，最后分析结冰过程中的传热传质，通过FLUENT提供的用户自定义函数(UDF)编程求解混合相热力学模型，实现混合相结冰预测。将数值计算结果与文献试验数据进行对比，分析计算方法的有效性。

1 混合相结冰热力学模型

1.1 质量守恒

图1 控制体质量守恒Fig.1 Mass balance of control volume

建立质量守恒方程：

(1)

(2)

同时，还有如下关系：

(3)

(4)

(5)

式中，U∞,ic和U∞,d分别为远场的冰晶速度和水滴速度；β为收集系数；ηic为冰晶融化比，ηic=LWCic/IWC，LWCic表示部分融化冰晶中液态水的含量，IWC表示空气流场中的冰晶含量；A为控制体的底面积；LWCd为空气流场中过冷水滴的液态水含量。

(6)

式中：hc为对流换热系数；ρ为空气密度；cp为比热容；Pr为普朗特数；hm为对流传质系数；Sc为施密特数，物理意义为动量扩散与质量扩散之比，即

(7)

式中：v和μ分别表示空气的运动和动力黏度；D为扩散系数。

利用气体状态方程可以推导出对流传质的质量流量计算式[18]，当表面有溢流水时为蒸发传质，表面仅为冰时为升华传质：

(8)

式中：Le为路易斯数，Le=Sc/Pr；MW为相对分子质量；Ts为壁面温度；pv,s,sat为壁面平衡温度下的饱和蒸汽压；pT为总压；TT为总温；pe为边界层外边界的空气压力；rh为相对湿度，本文中取值为1。

1.2 能量守恒

建立能量守恒方程：

(9)

能量守恒方程式(9)中，每一项都有特定的计算式，首先有进入控制体各质量项的动能：

图2 控制体内能量守恒Fig.2 Energy balance of control volume

(10)

(11)

式中：Uimp,d和Uimp,ic分别为水滴和冰晶粒子撞击到结冰表面时候的速度。

控制体内涉及水的冻结、水的蒸发和冰的升华3种相变过程(水的蒸发和冰的升华取其一，表面有溢流水则为蒸发散热，表面只有冰时则为升华散热)，因此存在相变潜热项：

(12)

(13)

(14)

式中：Lf为结冰潜热，取值333 kJ/kg；Lev为蒸发潜热，取值2 500 kJ/kg；Lsub为升华潜热，取值2 833 kJ/kg。

对流换热引起的能量传递：

(15)

式中：Tinf为自由流温度。

(16)

(17)

(18)

(19)

(20)

式中：Tm为融化温度，即0 ℃；cp,w为水的比热容，设为常数4 174 J/kg·K；cp,ic为冰的比热容，设置为常数2 102 J/kg·K。

2 冰晶黏附模型

通常情况下认为过冷水滴撞击表面后会全部停留在表面参与结冰过程，而冰晶存在很大不同，试验表明冰晶撞击表面后是否黏附与液态水的存在有很大关系[19]。冰晶撞击表面可能的结果有：反弹、破碎和黏附[20-22]。Baumert等[23]的冰晶撞击NACA0012翼型表面的试验表明反弹后的粒子除驻点附近极少能再次撞击到壁面，对结冰过程影响不大，因此主要关注冰晶撞击后的黏附效应。

Trontin等[24]根据冰风洞试验结果，包括纯冰晶部分融化条件结冰试验以及同时含有过冷水滴和冰晶的混合相结冰试验，分析液态水含量对冰晶黏附撞击表面可能性的影响，建立起符合统计学规律和冰风洞试验的经验公式。Trontin等[24]指出根据加拿大NRC (National Research Council) 的试验结果，相同LWC/TWC时(TWC为冰水总含量)，部分融化的纯冰晶条件下的黏附系数要高于混合相条件下的黏附系数，因此引入了εs,ic和εs,d，分别代表只考虑融化冰晶且无液态过冷水滴的影响以及考虑有液态水滴和(部分融化)冰晶共同的影响，黏附系数ε取两者中较大者：

ε=max(εs,ic,εs,d)

(21)

式中：

(22)

εs,d=Kd(Yd+ηicYic)

(23)

其中：ηic为冰晶融化比，根据Currie等[25]的纯冰晶环境下的试验，这里取0<ηic<0.2；Yd为撞击到表面的过冷水滴质量占撞击到表面的总质量的比值；Yic是撞击到表面的冰晶质量占撞击到表面的总质量的比值；Yd和Yic两者之和为1；Kic和Kd是常数。

(24)

(25)

考虑到表面温度越低时，冰晶黏附越少，特别是温度较低的干态冰表面时，冰晶不会黏附。于是，黏附模型中的常数Kd和表面温度有关：

(26)

可以发现，式(23)中的Yd+ηicYic代表了撞击到表面的所有液态水(包括过冷水滴和融化的冰晶)占总质量的比值。因此当混合相中冰晶不发生融化时，ηic=0，则ε=max(εs,ic,εs,d)=εs,d=KdYd。

本文计算的混合相结冰算例中环境温度较低，认为冰晶在运动过程中不发生融化，LWCic=0，式(21)的黏附系数ε则简化为ε=KdYd。

考虑冰晶黏附效应后，撞击冰晶中实际黏附在结冰表面参与结冰过程的冰晶质量流量为

(27)

3 模型求解方法

3.1 整体计算方法

结冰热力学模型用C语言编写，并通过UDF挂靠到FLUENT中完成编译和执行。首先计算空气流场，在空气流场的基础上计算水滴和冰晶的流场，获得水滴和冰晶的收集系数，随后把计算得到的剪切力、对流换热系数、压强、收集系数等作为输入参数计算结冰热力学模型，获得最终的结冰量。计算流程如图3所示。

图3 结冰计算流程Fig.3 Icing calculation flowchart

3.2 结冰热力学模型求解方法

对于二维翼型，认为驻点处流入控制体的质量流量为0，驻点所在控制体内的水平均分为两部分分别向上下表面溢流；流入流出控制体的水之间存在约束关系，即前一个控制体流出水的质量流量等于后一个控制体流入水的质量流量：

(28)

(29)

混合相结冰热力学模型求解的具体方法如下：将结冰表面从驻点处分为上表面和下表面两部分，首先计算驻点所在的控制体，然后分别针对上表面和下表面按顺序逐个求解驻点下游的控制体，直到表面水溢流极限处，完成全部控制体的求解。对每一个控制体来说，结冰过程热力学模型的求解通过迭代过程来实现。定义流入流出控制体能量的净值Q*=进入控制体的热量-流出控制体的热量，即能量守恒方程左右相减的差值，通过迭代使得Q*为0即满足控制体能量守恒，完成能量守恒方程的求解。

(30)

式中：t是结冰时间；ρice是冰的密度；A是控制体的面积(二维情况下为长度)。

4 计算结果及对比验证

取Al-Khalil在Cox冰风洞做的混合相结冰试验作为验证算例[19]。该试验采用了NACA0012翼型，弦长0.914 4 m，水滴平均直径为20 μm，冰晶平均直径分别为150 μm和200 μm。结冰试验条件如表1所示。其中算例1和算例2环境温度较高，算例2比算例1增加了冰晶含量，其他参数完全相同；算例3和算例4环境温度较低，两者的液态水含量和冰晶含量不同，算例3的液态水含量少冰晶含量多，而算例4的液态水含量多冰晶含量少。

表1 结冰计算条件Table 1 Calculation conditions of icing

4.1 空气流场计算

空气流场的计算划分结构网格、翼型周围网格如图4所示。网格数量为24 477个，计算域半径10 m，保证远场边界条件的成立。选取了底层网格高度0.1 mm、0.05 mm、0.005 mm和0.001 mm进行网格无关性计算，最终确定底层网格为0.005 mm的网格能达到要求，此时壁面附近的y+≈1。空气为理想气体，采用转捩剪切应力输运(SST)湍流模型，假设表面等效沙粒高度为1 mm，使用SIMPLE算法，为了提高计算精度，除压力采用标准格式外，其他参数全部采用二阶迎风格式。

图4 计算网格Fig.4 Computation mesh

由于计算参数相近，4个算例的计算结果十分相近，算例1中翼型周围空气速度和温度分布如图5所示。由于攻角为0°，速度和温度分布均为上下对称，驻点处速度为0 m/s，温度最高，较远场处升高2 ℃，但仍低于0 ℃，这也说明冰晶在运动过程中不可能会发生融化，前面假设冰晶运动过程中不融化是合理的。

4.2 冰晶/水滴撞击特性计算

混合相中既有过冷水滴又有冰晶，本文有如下假设：① 水滴或冰晶粒子在空气中的含量小(小于0.7 g/m3)，各自以及混合相的容积分数(α=Vpriticle/Vmixture)在10-6或10-7量级，可认为水滴和冰晶互不影响，空气和水滴以及空气和冰晶之间作用均是单向耦合，忽略粒子对空气流场的影响；② 环境温度较低，忽略水滴、冰晶运动过程中的变形、蒸发、升华、冰晶融化、粒子碰撞融合等现象；③ 远场处粒子的初始速度和周围空气相同。

图5 空气速度和温度分布云图Fig.5 Contours of air velocity and temperature distribution

水滴粒子的撞击特性计算采用欧拉法，基于利用课题组在FLUENT基础上二次开发的用户自定义标量的方法，详见文献[26]，在此仅作简单描述。在空气流场计算收敛的基础上，自定义了包括水滴容积分数、水滴x、y方向速度的3个标量；阻力系数采用球形阻力系数，对水滴运动的质量、动量方程进行描述，利用FLUENT求解器进行求解，采用二阶迎风格式；为了不让水滴在撞击壁面处堆积使得计算发散，壁面采用自定义的排水边界。局部水滴收集系数计算式为

(31)

式中：αd和α∞,d分别是壁面当地和自由来流中水滴容积分数；ud是壁面当地水滴速度矢量；n是壁面单位法向量。

为了尽可能接近真实情况，水滴撞击特性的计算采用7种不同直径的水滴计算结果取加权平均：

(32)

式中：pd,i为第i种直径水滴所占的液态水含量比重；βd,i是相应水滴直径的局部水收集系数。水滴直径d分布如表2所示。

表2 水滴直径分布Table 2 Distribution of droplet diameters

算例1的计算结果如图6所示，图中s表示弧长，c表示弦长。其他算例计算结果相近。水滴直径显著影响撞击极限和局部收集系数大小。

图6 水滴局部收集系数Fig.6 Local collection efficiency of droplets

冰晶的撞击特性计算中，由于本文对比文献中的试验温度较低，故不考虑冰晶运动过程中的融化。由于文献试验中采用的冰晶粒子接近球形，故本文针对冰晶的撞击特性计算方法和水滴撞击特性计算的方法类似，冰晶运动中的阻力系数采用了Clift和Gauvin的球形粒子阻力系数模型[27]。冰晶粒子的局部撞击系数计算式为

(33)

算例2和算例3中冰晶的局部收集系数结果如图7所示，冰晶的直径越大，冰晶撞击翼型的范围越广。当量直径为200 μm和150 μm的冰晶收集系数峰值相差不大。算例3和算例4的飞行条件以及冰晶当量直径一样，其冰晶收集系数结果也一致。

图7 冰晶局部收集系数Fig.7 Local collection efficiency of ice crystals

4.3 黏附效应及结冰冰形计算

算例1不含冰晶，因此无冰晶黏附。算例2～算例4冰晶黏附系数和局部撞击系数的计算结果如图8所示。从黏附系数结果可以看出，仅算例2有冰晶黏附，算例3和算例4无冰晶黏附，说明算例3和算例4中冰晶撞击壁面后全部反弹，不参与结冰过程。

无冰晶的算例1表面温度和溢流水结果如图9(a)所示。算例2有冰晶黏附是因为此状态下环境温度相对较高，从图9(b)可看出算例2表面收集的水并没有在撞击处全部冻结，驻点附近(-0.03

图8 冰晶局部黏附系数Fig.8 Local adhesion efficiency of ice crystals

存在冰晶黏附，黏附系数小于0.2。算例3和算例4无冰晶黏附的原因是这两个状态下环境温度相对较低，从它们的表面温度和溢流水结果即图9(c) 和图9(d)中看到，撞击水在撞击处能够全部冻结而无溢流水，使得表面形成干态冷表面，此时黏附系数模型中Kd=0，没有形成冰晶黏附条件。换句话说，算例3和算例4中的冰晶量并没有对图10是本文冰形预测结果与文献[18]试验结果的对比。为了体现黏附模型的作用，本文预测结果展示了是否考虑黏附模型两种情况。不考虑冰晶黏附效应时，认为撞击到机翼表面的冰晶全部黏附并参与结冰过程，忽略了冰晶可能出现的反弹飞溅。从图中可看到对于不含冰晶的算例1，计算结果与试验结果相差不大，其余算例中无黏附模型的冰形计算结果均远大于试验结果。考虑黏附效应后，计算结果得到极大改善。这说明冰晶的黏附效应是影响混合相结冰的重要因素，数值模拟时必须考虑。

算例1和算例2因为表面存在溢流水，在溢流区表面温度为273.15 K(图9)，结冰为明冰。从图10(b)可以看到混合明冰结冰气象条件下所结冰偏向楔形状而不是槽状，这和Baumert等[23]混合相结冰冰形一般是楔形冰的结论一致。另外，算例1计算结果较试验结果偏小，而算例2偏大。分析原因如下，算例1与算例2的结冰条件不同之处在于算例2中增加了冰晶含量，从图10(a) 和图10(b)中可看到算例2的试验结冰量反而较算例1的略有减小，试验现象说明了混合相条件下除了冰晶黏附效应还存在冰晶侵蚀效应，这说明本文液态水结冰量预测值偏小，加上本文计算未考虑冰晶侵蚀效应，因此算例1冰形预测值偏小，而算例2预测值偏大。冰晶侵蚀效应机理还不明确，目前尚未有合适的模型。

图9 表面温度和溢流水结果Fig.9 Results of surface temperature and runback water

图10 冰形计算结果对比验证Fig.10 Comparison and validation of computational results of ice shape

算例3和算例4表面无溢流水，表面温度低于273.15 K(图9)，结冰为霜冰。从图10(c)和图10(d)可看出霜冰条件下的仿真结果与试验结果吻合较好。

5 结 论

本文针对混合相结冰问题，建立了混合相结冰热力学模型和冰晶黏附模型，利用FLUENT软件及其用户自定义函数开展了混合相结冰模拟，和试验结果对比说明了方法的有效性。主要结论如下：

1) 基于Messinger模型，增加冰晶有关质量和能量相，建立了适用于混合相结冰的结冰热力学模型。

2) 分别计算了NACA0012翼型在霜冰条件和明冰条件下的结冰情况，仿真结果与试验结果具有较好的一致性，说明了本文计算模型和方法对混合相结冰的适用性。混合明冰结冰气象条件下所结冰偏向楔形状而不是槽状。

3) 黏附效应是影响混合相结冰的关键因素，添加黏附效率经验公式能有效改善仿真计算的准确性。本文的黏附模型是经验公式，还有待从机理上进一步开展黏附模型的试验研究，以获得更精确的黏附模型。

4) 非球形冰晶在运动过程中的传热传质模型有待于增加到现有欧拉模型中，从而扩展对热环境下的非球形冰晶运动及黏附结冰的计算。

5) 冰晶的侵蚀效应是区别于过冷水滴结冰的重要方面，其理论模型有待深入研究。