过冷大水滴运动和撞击特性数值模拟

摘 要：为了更准确地模拟过冷大水滴在结冰部件中的撞击特性，基于LEWICE模型和Trujillo模型的过冷大小滴撞击模型，以NACA23012翼型为研究对象，通过Fluent软件的二次开发功能，考虑大水滴运动中变形以及撞击中反弹/飞溅的影响，对过冷大水滴的运动和撞击特性进行数值模拟，并将计算结果与NASA的试验数据进行对比。结果表明：当过冷水滴的MVD值较大时，相比于LEWICE模型和Trujillo模型，结合两个模型计算得到的翼型表面局部水收集系数的计算结果和撞击范围都更贴合试验数据，绝大部分计算结果和试验数据的误差在10%以内；MVD值越大，局部水收集系数值的计算结果和试验值越贴近。

关键词：过冷大水滴；数值模拟；变形；反弹；飞溅

中图分类号：TP391.9 文献标志码：B 文章编号：1671-5276（2024）04-0140-06

Numerical Simulation of Motion and Impingement Characteristics of Supercooled Large Droplets

LI Zhonghui， LI Shiming， LIU Lei， YU Dongyin

（College of Energy and Power Engineering， Nanjing University of Aeronautics and Astronautics， Nanjing 210016， China）

Abstract：In order to more accurately simulate the impingement characteristics of supercooled large droplets on frozen parts， based on LEWICE model，Trujillo model and two SLD impingement models， taking NACA23012 airfoil as the research object and with the secondary development function of Fluent software， conducts the numerical simulation on the motion and impingement characteristics of SLD in consideration of the influence of deformation during movement and rebound/splash during impingement. The calculated results are compared with the experimental data of NASA. The results show that when the MVD value of SLD is large， compared with LEWICE model and Trujillo model， the calculation results of local collection efficiency of airfoil surface and impingement range calculated by combining the two models are more consistent with the test data， and the error of most of the calculation results is less than 10% compared with the test data; the larger the MVD is， the closer the calculation result of the local collection efficiency value is to the test value.

Keywords：supercooled large droplets; numerical simulation; deformation; bounce; splash

0 引言

过冷大水滴（supercooled large droplets， SLD）是指云层中平均容积直径（median volume diameter，MVD）大于50μm的过冷水滴。飞机和发动机在含有过冷大水滴的环境中飞行时，会发生结冰现象。过冷大水滴在运动和撞击过程中，会发生变形、破碎、反弹、飞溅等复杂现象，所产生的二次小水滴可能会再次撞击^[1]。过冷大水滴的这些特性使其撞击范围和结冰范围较大，也容易产生溢流水，并导致在防护区域外形成脊状冰，产生安全隐患。水滴的运动和撞击特性是结冰和防冰预测的前提，而常规尺寸水滴运动和撞击特性的模拟方法通常不需要考虑水滴的破碎、飞溅和反弹等现象，但并不适用于大水滴。因此对过冷大水滴运动和撞击特性的研究是十分有必要的。

国内外研究学者对此开展了一系列理论、数值和试验研究。为了描述过冷大水滴运动和撞击过程中出现的变形、破碎、反弹和飞溅等现象，国外学者结合试验数据建立了相关的经验或半经验模型^[2]。针对过冷大水滴运动过程中发生的变形，1933年LIU等^[3]、1978年CLIFT等^[4]分别提出了两种过冷大水滴阻力系数的模型。针对过冷大水滴撞击壁面过程中发生的反弹、飞溅现象，WRIGHT、MACKLIN等^[5-6]提出了相应的计算模型，分别被LEWICE软件^[7]、FENSAP-ICE软件^[8]采用，TRUJILLO等^[9]也提出了类似的飞溅模型。2010年，IULIANO等^[10]采用Clift的阻力模型和不同的反弹/飞溅模型计算了SLD在翼型上的撞击特性。结果表明，采用Trujillo模型预测MVD为111μm的SLD撞击到NACA23012翼型表面的水收集系数峰值比试验值低约20%，采用LEWICE模型预测MVD为94μm的SLD撞击到MS（1）-0317翼型表面的水收集系数峰值高于试验值约30%。2021年，王龙等^[11]采用水平集法从微观角度研究水滴撞击过程，通过数值模拟研究半径1.3mm的过冷大水滴撞击特性，发现大水滴运动过程中存在回流。2004年，PAPADAKIS等^[12]在NASA冰风洞中获得了MVD分别为20μm、52μm、111μm、154μm和236μm的过冷水滴在NACA23012翼型表面的撞击特性试验数据。

如上所述，国内外学者计算得到的过冷大水滴撞击特性结果和试验数据相比，精度不够高。为了更准确地模拟过冷大水滴在结冰部件表面的撞击特性，本文采用考虑过冷大水滴运动过程中变形的Clift模型，结合考虑过冷大水滴撞击过程中反弹/飞溅的LEWICE模型和Trujillo模型，以NACA23012翼型为研究对象，通过Fluent软件的二次开发功能，对过冷大水滴的运动和撞击特性进行数值模拟，并将计算结果与NASA的试验数据进行对比分析。

1 数学模型和计算方法

1.1 过冷大水滴运动过程中的变形模型

本文采用欧拉法计算水滴流场。根据守恒定律可得水滴的连续方程和动量方程^[13]分别为：

式中：ρ_d为水滴密度；α为水滴的容积分数；u_d为水滴的速度矢量；u_a为空气的速度矢量；μ_a为空气的动力黏度；d为将水滴视为规则球体时的水滴直径；Re为以空气和水滴之间的相对速度作为特征速度的雷诺数，表达式为

式中ρ_a为空气密度。大水滴由于尺寸较大，其表面张力不容易克服所受的空气剪切力，因此在运动过程中会发生变形，水滴形状会由球状向类似扁圆盘状过渡。采用水滴韦伯数表征水滴惯性力和表面张力之间的相对大小，其表达式为

式中σ_d为水滴的表面张力。变形后的水滴受到的空气阻力发生改变，本文采用Clift提出的变形模型模拟大水滴形状改变所导致的阻力系数变化。该模型取球体阻力系数C_D，sphere和扁圆盘阻力系数C_D，disk的加权平均值来描述大水滴变形后的阻力系数C_D，其表达式为

式中e为插值系数，表达式为

式（5）中球体阻力系数C_D，sphere为

式（5）中扁圆盘阻力系数C_D，disk为

1.2 过冷大水滴撞击过程中的反弹/飞溅模型

过冷大水滴在撞击壁面后，会发生反弹和飞溅，其引起的水滴质量损失较大。反弹是指水滴在撞击到壁面后以一定角度离开壁面；飞溅是指水滴在撞击壁面后部分水滴停留在壁面，部分飞离壁面。本文使用LEWICE软件中SLD反弹/飞溅模型和TRUJILLO的SLD飞溅模型，模拟过冷大水滴撞击壁面反弹/飞溅现象。

1）LEWICE的模型

LEWICE软件应用了WRIGHT等提出的同时考虑了水滴的反弹和飞溅质量损失的大水滴撞击模型。模型中采用了MACKLIN等提出的与R_e，s数和O_h数有关的飞溅参数K，其表达式为

式中O_h数为

R_e，s数的表达式为

式中u_d，n为水滴速度在壁面上的法向分量。定义耦合水滴撞击能和入射频率的水滴撞击参数K_L，作为判断发生反弹、飞溅现象的衡量参数，其表达式为

式中：θ为水滴速度和撞击表面切向的夹角；f^*为无量纲水滴频率，定义为

由于反弹/飞溅效应引起的过冷大水滴撞击质量损失用无量纲质量损失系数f_m，lewice表示为

2）Trujillo模型

TRUJILLO等提出的大水滴撞击模型考虑了水滴飞溅引起的质量损失，满足以下条件时水滴飞溅发生^[14]：

式中：K_Ctr为判断飞溅发生的临界参数；R为撞击表面粗糙度相对水滴直径的比值，表达式为

式中R_s为撞击表面的最大粗糙度，其可由试验结果^[15]确定。该模型的无量纲质量损失系数f_m，trujillo为

1.3 计算方法

本文通过Fluent软件的二次开发功能，结合上述数学模型，对过冷大水滴的运动和撞击特性进行数值计算，包含以下3个模块。

1）空气流场计算。首先计算部件外空气流场，由于水滴含量很少，忽略其对空气的影响，使用Fluent软件的单相流求解器求解空气相的连续方程、动量方程以及相应的湍流方程，获得部件外空气流场的速度、压力等分布。

2）过冷水滴流场计算。在上一步计算好的空气流场基础上，加载编译用于求解水滴流场的UDF文件。该UDF文件中将过冷水滴的连续方程和3个坐标方向的动量方程视为4个标量输运方程，即使用Fluent软件的用户自定义标量UDS（user defined scalars）功能定义4个输运方程（UDS0—UDS3），UDS0—UDS3分别表示水滴容积分数和水滴在x、y、z 3个方向的速度分量，其中水滴动量方程的阻力系数应用Clift的SLD变形模型。以此对一定边界条件约束下水滴运动的连续方程和动量方程进行离散，迭代求解得到水滴流场的容积分数、水滴速度等的分布。

3）水滴撞击特性计算。在空气和过冷水滴两相流场结果的基础上，使用Fluent软件的用户自定义内存UDM（user defined memories）功能，依次求解得到质量损失系数f_in，lewice和f_m，trujillo，再结合上述两个SLD撞击模型求解得到最终的局部水收集系数β。

2 计算模型

2.1 计算域和边界条件

本文以NACA23012翼型为研究对象，其弦长为0.914 4m，展向高度为0.914 4m。建立如图1所示的翼型外部空气-过冷大水滴两相流计算域。该计算域的进口边界为速度进口，出口边界为压力出口，压力为1个大气压。本文在来流速度、温度和液态水含量均相同的情况下计算不同MVD的过冷大水滴运动和撞击特性，具体计算边界条件参数如表1所示。

2.2 网格划分及独立性验证

采用ANSYS ICEM商业软件对计算域进行网格划分，全域划分六面体结构化网格，并对翼型壁面网格进行加密处理，如图2所示。由于空气流场计算采用k-ε湍流模型并结合Realizable壁面函数，划分翼型壁面附近第一层网格尺度为0.2mm，对应的y⁺约为30。

为了消除网格疏密对计算结果的影响，分别采取单元总数分别为12万、37万、76万、101万和129万的5套网格进行两相流场计算，对比不同计算网格下图2所示的监测点处水滴速度，以此进行网格独立性验证。图3为监测位置处水滴速度随网格数量变化的结果。由图3可知，当网格总数大于76万时，监测点的水滴速度不再随网格数量变化，因此本文最终选用总数101万的网格进行计算。

3 计算结果及分析

为描述方便，在计算域展向的中部垂直于展向截取一个中截面如图4所示，因该翼型沿展向的各截面均相同，故分析该中截面上的流场计算结果即可。

3.1 运动特性分析

根据本文的计算方法，计算空气流场时忽略过冷水滴的影响，故6个算例的空气流场计算结果相同。计算域中截面上的绝对压力分布和空气速度分布分别如图5和图6所示。由图可知，翼型驻点处的压力最大，约为106kPa，沿着翼型上下表面，压力先降低后上升，但并不对称，且在相同的流向坐标处，上表面的压力较小。

虽然6个算例的空气流场完全一致，但加入过冷大水滴后，因水滴的直径不同，则空气对水滴的气动力不同，导致水滴的速度分布和容积分数分布不同。6个算例中截面上的水滴速度分布如图7所示（本刊黑白印刷，相关疑问咨询作者）。

图7（a）所示为MVD值为52μm时的水滴速度分布，驻点附近速度较低，沿流向往下游，速度先上升后下降，但相同位置处水滴速度值略小于空气速度值。与空气速度分布相比不同之处在于，翼型上部水滴速度值较大的红色区域呈扁长状，沿流向延伸至翼型的中后部。由图7（a）—图7（f）可见，随着MVD值增大，翼型外部水滴速度值较高的红色区域减小，为水滴的最大速度值也降低。这是因为随着过冷大水滴尺寸增大，水滴惯性增大，更不容易被空气带动，其流动方向趋向于直线。

不同MVD值下翼型附近水滴体积分数分布如图8所示。由图8（a）可见，MVD值为52μm时，在稍远离翼型的流场区域，水滴体积分数值与进口处基本相同。在紧贴翼型壁面处，沿着流向从驻点向下游，水滴的体积分数减小。越靠近翼型壁面，水滴的体积分数值越小，直至体积分数为0的遮蔽区。对比图8（a）—图8（f）可见，随着MVD值增大，水滴体积分数大的红色区域减小，这是因为随着水滴直径增大，水滴惯性增大，更多的水滴撞击到翼型壁面并被收集。由于翼型前部的遮挡，导致翼型中后部附近过冷水滴无法到达，形成体积分数几乎为0的遮蔽区。

3.2 过冷大水滴撞击特性分析

取计算域中截面与翼型壁面相交的轮廓线上的局部水收集系数，并与Papadakis的试验数据进行对比，如图9所示。图中横坐标y/c表示翼型表面从驻点处开始的弧长坐标y与弦长c之比，y/cgt;0为翼型上表面，y/clt;0为翼型下表面。由图9可知，水滴撞击集中在翼型前缘，翼型驻点处局部水收集系数最大，从驻点向后，局部水收集系数降低，在翼型前缘下游，局部水收集系数迅速降低为0。从图9（a）—图9（f）可以看出，随着MVD值增大，水滴的撞击范围略微增大。这也是由于水滴直径增大使得其惯性增大所致。

由图9（a）可见，MVD值为52μm时，翼型下表面和驻点附近水收集系数的计算值比试验数据略高，两者存在少许偏差，其中驻点处相差约0.1。从驻点附近往下游，翼型上表面水收集系数计算值和试验值的重合度高。由图9（b）—图9（f）可知，本文的计算结果与试验数据总体符合较好，在翼型表面的绝大部分范围内，计算值和试验值之间的差值均低于10%。

4 结语

为了更准确地模拟过冷大水滴在结冰部件表面的撞击特性，以NACA23012翼型为研究对象，通过Fluent软件的二次开发功能，采用Clift过冷大水滴变形模型，并基于LEWICE模型和Trujillo模型的两个过冷大水滴反弹、飞溅模型，对过冷大水滴的运动和撞击特性进行了数值模拟。计算得到了流场的空气速度、压力分布，水滴的速度、体积分数分布以及翼型表面局部水收集系数分布，并将局部水收集系数分布结果与试验数据进行对比，验证了本文计算模型和计算方法的合理性。本文主要研究结论如下。

1）空气-过冷大水滴两相流在翼型驻点附近速度较低，沿流向往下游，速度先上升后下降。相同位置处，过冷水滴速度值略低于空气速度值。翼型上部水滴速度值较大的红色区域呈扁长状，沿流向延伸至翼型的中后部。

2）其他条件相同的情况下，随MVD值增大，翼型外部过冷水滴高速度值区域范围减小，水滴速度最大值降低，水滴流线越来越趋近于直线。

3）翼型附近外水滴体积分数沿流向从驻点向下游逐渐减小，越靠近翼型壁面，水滴的体积分数值越小，直至为0。随MVD值增大，水滴体积分数较高的区域减小。

4）过冷大水滴撞击集中在翼型前缘，翼型驻点附近局部水收集系数最大，从驻点向后，局部水收集系数降低，在翼型前缘下游局部水收集系数迅速降低为0。

5）当过冷大水滴的MVD值较大时，与LEWICE模型和Trujillo模型相比，结合两个模型计算所得的翼型表面局部水收集系数计算结果和撞击范围与试验数据符合较好，MVD值越大，计算值与试验值越吻合，由此验证了采用本文的计算模型和计算方法在计算过冷大水滴运动与撞击特性时具有适用性和有效性。

参考文献：

[1] 林贵平，卜雪琴，申晓斌. 飞机结冰与防冰技术[M]. 北京：北京航空航天大学出版社，2016.

[2] WRIGHT W，POTAPCZUK M. Semi-empirical modelling of SLD physics[C]//42nd AIAA Aerospace Sciences Meeting and Exhibit. Reno，Nevada. Reston，Virginia：AIAA，2004：412.

[3] LIU A B，MATHER D，REITZ R D. Modeling the effects of drop drag and breakup on fuel sprays[C]//SAE Technical Paper Series. 400 Commonwealth Drive，Warrendale，PA，United States：SAE International，1993：930072.

[4] CLIFT R，GRACE J R，WEBER M E. Bubbles，drops，and particles[M]. New York：Academic Press，1978.

[5] WRIGHT W，POTAPCZUK M，LEVINSON L. Comparison of LEWICE and GlennICE in the SLD regime[C]//46th AIAA Aerospace Sciences Meeting and Exhibit. Reno，Nevada. Reston，Virginia：AIAA，2008：439.

[6] MACKLIN W C，METAXAS G J. Splashing of drops on liquid layers[J]. Journal of Applied Physics，1976，47（9）：3963-3970.

[7] WRIGHT W. User’s Manual for LEWICE Version 3.2[M]. [S.I.]：Nasal/Cr， 2008.

[8] HONSEK R，HABASHI W G. FENSAP-ICE：Eulerian modeling of droplet impingement in the SLD regime of aircraft icing[C]//44th AIAA Aerospace Sciences Meeting and Exhibit. Reno，Nevada. Reston，Virginia：AIAA，2006：465.

[9] TRUJILLO M F，MATHEWS W S，LEE C F，et al. Modelling and experiment of impingement and atomization of a liquid spray on a wall[J]. International Journal of Engine Research，2000，1（1）：87-105.

[10] IULIANO E，MINGIONE G，PETROSINO F，et al. Eulerian modeling of SLD physics towards more realistic aircraft icing simulation[C]//AIAA Atmospheric and Space Environments Conference. Toronto，Ontario，Canada. Reston，Virginia：AIAA，2010：7676.

[11] 王龙，汪宁，来永斌，等. 大水滴撞击翼型后的破碎过程研究[J]. 应用力学学报，2021，38（3）：1071-1078.

[12] PAPADAKIS M，RACHMAN A，WONG S C，et al. Water impingement experiments on a NACA 23012 airfoil with simulated glaze ice shapes[C]//42nd AIAA Aerospace Sciences Meeting and Exhibit. Reno，Nevada. Reston，Virginia：AIAA，2004：565.

[13] BOURGAULT Y，HABASHI W，DOMPIERRE J，et al. An Eulerian approach to supercooled droplets impingement calculations[C]//35th Aerospace Sciences Meeting and Exhibit. Reno，NV. Reston，Virginia：AIAA，1997：176.

[14] YARIN A L，WEISS D A. Impact of drops on solid surfaces：self-similar capillary waves，and splashing as a new type of kinematic discontinuity[J]. Journal of Fluid Mechanics，1995，283：141-173.

[15] ANDERSON D，HENTSCHEL D，RUFF G. Measurement and correlation of ice accretion roughness[C]//36th AIAA Aerospace Sciences Meeting and Exhibit. Reno，NV. Reston，Virginia：AIAA，1998：486.

收稿日期：2023-01-30