基于蒙特卡洛方法的冰晶撞击特性计算_马乙楗.pdf

第 55 卷第 2 期2023 年 4 月Vol.55 No.2Apr.2023南 京 航 空 航 天 大 学 学 报Journal of Nanjing University of Aeronautics&Astronautics基于蒙特卡洛方法的冰晶撞击特性计算马乙楗1,2，柴得林1，王强1,2，易贤1,2，余放3（1.中国空气动力研究与发展中心结冰与防除冰重点实验室,绵阳 621000；2.中国空气动力研究与发展中心空气动力学国家重点实验室,绵阳 621000；3.中国航发商用航空发动机有限责任公司,上海 200241）摘要：航空发动机摄入冰晶易发生内部结冰、造成发动机推力损失，严重时引发航空事故。冰晶结冰数值模拟是研究冰晶结冰的重要手段，而撞击收集系数计算是冰晶结冰数值模拟的关键环节。基于 NNWICE 平台开展冰晶结冰前期的阻力系数、运动轨迹和撞击收集系数的数值模拟与计算研究，实现了基于蒙特卡洛方法的冰晶撞击收集系数的准确计算；在此基础上发展了拉格朗日框架下冰晶粒子运动、传热传质耦合数值计算方法；通过完成二维 NACA0012翼型的相关算例，对比分析了不同阻力模型和冰晶相变对撞击收集系数的影响。结果表明阻力模型对冰晶撞击收集系数影响较小，冰晶相变过程会对冰晶运动轨迹及撞击收集系数产生显著影响，为进一步发展冰晶结冰数值模拟计算方法奠定了基础。关键词：冰晶结冰；撞击收集系数；蒙特卡洛方法；冰晶相变；数值模拟中图分类号：V211.41 文献标志码：A 文章编号：1005-2615（2023）02-0291-11Calculation of Ice Crystal Impact Characteristics Using Monte Carlo MethodMA Yijian1，2，CHAI Delin1，WANG Qiang1，2，YI Xian1，2，YU Fang3（1.Key Laboratory of Icing and Anti-icing，China Aerodynamics Research and Development Center，Mianyang 621000，China；2.State Key Laboratory of Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，China；3.China National Aero Engine Corporation，Shanghai 200241，China）Abstract:The intake of ice crystals by an aero-engine is prone to internal icing，causing loss of engine thrust and，in serious cases，aviation accidents.The numerical simulation of ice crystal icing is an important tool to study ice crystal icing，and the calculation of the collection coefficient is a key part of the numerical simulation of ice crystal icing.Based on the NNWICE platform，the numerical simulation and calculation of the drag coefficient，motion trajectory and collection coefficient of ice crystal icing in the early stage of ice crystal icing is carried out，and the accurate calculation of ice crystal collection coefficient based on Monte Carlo method is achieved.Therefore，the coupled numerical calculation method of ice crystal particle motion and heat and mass transfer under Lagrangian framework is developed；by completing the relevant calculation cases of two-dimensional NACA0012 airfoil，the effects of different drag models and ice crystal phase transition on the impact collection coefficient are compared and analyzed.The results show that the drag model has a little effect on the ice crystal collection coefficient，while the ice crystal phase transition will have a significant effect on the ice crystal motion trajectory and collection coefficient，which lays the foundation for further development of the numerical simulation calculation method of ice crystal icing.Key words:ice crystal icing；collection coefficiency；Monte Carlo method；ice crystal phase transition；numerical simulationDOI：10.16356/j.1005-2615.2023.02.015基金项目：国家自然科学基金重点基金（12132019）；国家重大科技专项（J2019-III-0010-0054）；国家自然科学面上基金（12172372）。收稿日期：2022-05-26；修订日期：2022-07-05通信作者：易贤，男，研究员，博士生导师，E-mail：yixian_。引用格式：马乙楗，柴得林，王强，等.基于蒙特卡洛方法的冰晶撞击特性计算 J.南京航空航天大学学报，2023，55（2）：291-301.MA Yijian，CHAI Delin，WANG Qiang，et al.Calculation of ice crystal impact characteristics using monte carlo methods J.Journal of Nanjing University of Aeronautics&Astronautics，2023，55（2）：291-301.第 55 卷南 京 航 空 航 天 大 学 学 报飞机结冰一直是影响飞行安全的重大问题之一。随着航空科技的发展，研究者发现在 6 km 以上含有极少过冷水的高空，航空发动机会发生失速、喘振，燃烧室熄火、停车以及叶片损伤等现象。这些现象是由该高度下云层所含冰晶进入推进系统的核心流动通道并发生结冰造成的1-2。由于模拟某一航空发动机在一定工况下冰晶结冰的试验设备异常复杂，不仅需要消耗大量人力、物力、财力，而且试验结果也仅适用于某些特定工况，没有普适性，因此有必要开展冰晶结冰数值计算方法研究3。国外经过多年的发展已开发出较为完善的结冰 软 件，如 美 国 的 GlennICE、LEWICE，法 国 的ONERA3D 和加拿大的 FEASAP 等，上述软件开发者主要通过改进拓展过冷水滴结冰的计算方法发展了相应软件对冰晶结冰情况的数值模拟能力4-6。而国内由于冰晶结冰研究起步晚，尚缺乏相应的结冰软件，仅有部分学者开展了对冰晶结冰问题的有限研究，如文献 3，7 系统总结了航空发动机中冰晶结冰的研究现状；文献 3，8 利用 Flu-ent 软件结合自定义函数（UDF）进行了冰晶结冰模拟；文献 9 对大型结冰风洞中冰晶的热/力平衡特性数值进行了研究。冰晶结冰数值研究主要包括冰晶运动、撞击等动力学数值研究以及冰晶融化、冻结等热力学数值研究。其中冰晶轨迹以及撞击收集系数计算是冰晶结冰计算研究的重点，也是本文的主要研究内容。冰晶在运动过程中，一方面受到空气阻力的作用，另一方面与空气发生传热传质，而阻力和传热传质又与冰晶自身形状密切相关，这些因素共同影响冰晶的运动轨迹及撞击收集系 数。针 对 阻 力 计 算，Pitter、Haider、Ganser 和Hlzer 等10-13相继提出了颗粒在空气中流动的不同阻力模型。针对传热传质，Iuliano 等14、Trontin等6，15和 Norde 等16均对自由流中的非球形冰晶颗粒传热传质进行了研究，发展并不断完善冰晶的传热传质模型。针对冰晶形状，Zhang等17假设冰晶为球形，冰晶在运动过程中形状、密度等参数不发生改变；郭向东等9假设冰晶为六棱柱和六角板状，研究了冰晶形状对其运动过程中换热和轨迹的影响；Nilamdeen 等18假设冰晶为六角平板状，并提出冰晶在流场的流动特性类似于圆盘状颗粒。此外，冰晶在空气中的含量较高且粒径较大，是否考虑冰晶与空气的相互影响是冰晶运动轨迹计算中亟待解决的问题。Ros 等19研究了尺寸为 50、100、200 m的非球状冰晶运动时，仅采用空气对冰晶影响的单向耦合的方法；Nilamdeen 等20认为虽然冰晶在空气中含量高，但满足单向耦合计算条件。目前，国内学者仅对冰晶运动及撞击收集系数进行了初步研究3，7，尚未形成包括冰晶形状、阻力模型、冰晶传热传质等过程的系统的撞击收集系数求解方法。本文依托 NNWICE 平台21，在拉格朗日框架下发展了针对典型非球状冰晶颗粒运动轨迹的数值计算方法。采用可克服粒子运动轨迹交叉、鲁棒性良好的蒙特卡洛撞击收集系数计算方法实现了冰晶局部撞击收集系数的计算。最后通过二维NACA0012 翼型算例验证了上述数值计算方法和3 种冰晶阻力模型的有效性和准确性。并在此基础上发展了考虑冰晶运动过程中传热传质现象的单向耦合数值计算方法，指出了冰晶相变对撞击极限及撞击收集系数的影响规律。1 计算方法本文基于 NNWICE平台采用现有拉格朗日框架下的蒙特卡洛方法，引入非球状冰晶粒子的传热传质模型、冰晶阻力模型，发展拉格朗日框架下冰晶粒子运动、传热传质耦合计算方法。为简化冰晶结 冰 计 算 问 题，针 对 上 述 模 型 提 出 以 下 假设15-16，22：（1）冰晶颗粒在运动过程中不发生旋转，以动力学最稳定状态运动；（2）冰晶粒子在运动过程中不发生碰撞、聚合；（3）冰晶发生融化后融化的水膜附着在冰晶粒子表面，在运动过程中水膜不与粒子分离也不发生飞溅；（4）冰晶粒子内部温度保持恒定，粒子表面仅存在对流换热过程，忽略辐射换热；（5）冰晶运动过程中，液态冰晶水的密度和固态冰晶水的密度保持恒定；（6）冰晶达到融化温度后开始融化，并且在融化过程中冰水混合物温度保持为融化温度不变。冰晶撞击收集系数计算主要分为以下 5 个步骤：网格划分、流场求解、设置云雾场条件、冰晶运动轨迹计算和冰晶撞击收集系数计算，具体求解流程如图 1所示。图 1 冰晶撞击收集系数计算流程图Fig.1 Calculation flowchart of ice crystal collection coeffi-ciency292第 2 期马乙楗，等：基于蒙特卡洛方法的冰晶撞击特性计算1.1冰晶运动模型1.1.