么 虹，宁义君
（中国航空工业空气动力研究院 高速高雷诺数气动力航空科技重点实验室，沈阳110034）
摘 要：为了掌握发动机进气流量对水滴撞击特性的影响规律，采用欧拉方法建立水滴运动方程，基于分区多块结构网格体系求解水滴运动轨迹，提出了复杂外形飞机水滴撞击特性的计算方法.基于该方法对发动机进气道唇口进行水滴撞击特性分析，研究了不同发动机进气流量对进气道唇口水滴撞击特性的影响.结果表明，发动机进气流量越大，水滴收集系数最大值越大.因此，航空发动机吊舱进气道前缘的防除冰系统设计必须考虑不同飞行状态下进气流量的影响.
关 键 词：发动机；吊舱；欧拉法；过冷水滴；撞击特性；防除冰系统；结构网格；数值模拟
当飞机在含有过冷水滴的云层中飞行时，飞机迎风表面诸如机翼前缘、发动机进气口和雷达整流罩等部件上会出现结冰现象，而结冰现象的产生又有很大的随机性、突然性.飞机结冰是导致飞行安全事故的主要隐患，危害性巨大.飞机发动机吊舱进气道前缘（唇口）结冰不仅直接改变进气道前缘气动外形，减少发动机进气流量，降低发动机推力，且当进气道前缘冰层脱落时，冰块将随着气流进入发动机内部从而损伤叶片.发动机唇口结冰是飞机结冰中最危险的情况之一，会对飞行安全造成极大威胁，因此，发动机进气道的防除冰系统设计对飞机安全飞行而言极其重要.
当进行飞机防除冰系统的防护面积确定和结冰外形数值预测时，首先需要求解水滴运动轨迹，获得结冰防护表面的水滴撞击特性.目前水滴运动方程的建立主要有两种方法，一种是拉格朗日法，即根据牛顿第二定律得到单个水滴的运动方程，由该方法获得的运动方程形式简单，便于求解，适于二维简单翼型外形的计算；另一种是欧拉方法，在欧拉坐标系下建立水滴运动控制方程，通过求解控制方程得到流场中水滴的分布，进而获得水滴收集系数.在三维或复杂外形的水滴收集系数计算方面，欧拉方法相比拉格朗日法具有明显优势，近年来得到了较大发展.
飞机表面水滴撞击特性对于防除冰系统设计及其性能分析而言非常重要.国外学者基于二维翼型模型、轴对称模型和三维飞机部件模型等开展了大量的数值模拟方法和试验测量方法等方面的研究工作.Langmuir首次描述了水滴运动无量纲控制方程，并计算了圆柱和圆球的水滴撞击特性 ［1］.基于Langmuir简单几何模型的水滴计算工作，Gelder首次在风洞内测量了某航空发动机进气道的水滴收集系数 ［2］.20世纪50年代，Brun基于不可压缩流动理论计算了航空发动机进气道水滴撞击特性 ［3］.从早期发表的文献上看，只有很少的水滴计算研究工作是基于三维发动机进气道开展的.Kim基于可压缩势流理论计算了三维发动机进气道的水滴运动轨迹 ［4］；美国航空航天局开发了专业飞机结冰计算分析软件LEWICE系列 ［5－6］；Bidwell采用LEW ICE3D飞机结冰计算软件计算了轴对称进气道和波音 737-300飞机发动机进气道的水滴撞击特性和结冰冰形 ［7］.
国内针对飞机结冰数值预测方面的研究始于90年代.根据一些经验公式进行过冷水滴撞击、结冰情况与防除冰问题的估算，并在2000年左右开始采用数值模拟手段研究飞机结冰问题.目前针对二维翼型等简单模型水滴撞击特性的研究已较为成熟，但对于飞机进气道等复杂外形结冰问题的研究还较少，且一般将进气道简化为二维模型进行计算分析.顾海君 ［8］采用二维欧拉坐标系下的空气-过冷水滴两相流控制方程计算研究了某型发动机进气道唇口的水滴撞击特性；杨倩等 ［9］采用拉格朗日法模拟了某型直升机进气道前缘处的水滴撞击特性；申晓斌等 ［10］将欧拉方法与计算流体力学软件 FLUENT相结合实现了可压缩流场中三维复杂外形表面水滴收集系数的计算，并成功应用于三维发动机进气道水滴撞击特性的计算和研究.
本文利用中航工业空气动力研究院自主研发的结构网格并行 CFD计算平台 ENSMB ［11］并基于欧拉法建立三维水滴运动方程，在分区多块结构网格体系下求解水滴运动轨迹，实现了三维复杂外形飞机表面可压缩流场水滴撞击特性的计算，并建立了航空发动机进气道水滴撞击特性分析软件系统.基于该软件系统，本文对发动机进气道唇口开展三维水滴撞击特性计算分析，研究了不同发动机进气流量对进气道唇口水滴收集系数的影响.
运输类飞机结冰适航条例规定云层中过冷水滴平均有效直径一般低于50μm.通常采用水滴体积因子α描述云中过冷水滴体积，α的量级一般约为10 －6.可以认为极小水滴的存在和水滴体积因子不影响空气流场的流动过程，可以采用如下假设建立水滴运动控制方程：
1）与空间流场相比云中水滴总体积很小，不影响空气流场的流动；
2）云中水滴尺寸在微米量级，尺寸足够小，且均匀分布，可以将水滴作为具有等效直径的球体，且该球体在运动过程中不变形、不破碎、不发生碰撞或凝聚，且与物面撞击时不发生飞溅等现象；
3）水滴在流场中的运动过程不发生热交换和蒸发现象，且水滴温度、密度等物性参数不变；
4）流场中的湍流脉动不影响水滴运动；
5）水滴仅受到空气阻力、重力和空气浮力等外力作用.
基于以上假设，空气和水滴控制方程无耦合关系，即空气流动控制方程与水滴运动控制方程不发生耦合.空气流场可以单独求解，之后在空气流场数值计算结果基础上求解水滴运动方程.
水滴运动控制方程包括连续性方程和动量方程，其无量纲形式分别为
式中：ρ为液态水密度；ρ a为空气密度；u为水滴速度矢量；u a为空气速度矢量；g为重力加速度；K为惯性参数；F r表征流体惯性力和重力的相对大小，记作Froude（弗劳德）数.
动量方程右边第一项为水滴所受到的空气阻力，其与水滴相对速度、阻力系数C D和水滴雷诺数成正比.水滴雷诺数表达式为
式中：d为水滴中值直径；U a，∞为自由来流速度；μ a为空气分子粘性系数.
惯性参数K表达式为
式中：L
∞为物体的特征尺寸；ρ
d为水滴密度；f为阻力函数，且
动量方程右边第二项代表浮力和重力，且与当地 Froude数成正比，且 Froude数表达式为
2.1 积分形式控制方程
将式（1）、（2）改写成积分形式，即
式中：Ω为控制体体积；S为控制体表面积.式（6）中相应变量表达式为
式中，v、w、v a和w a为相应参量.将式（6）在进行空间流场计算的结构网格上进行离散求解.
2.2 边界条件
水滴的自由来流边界条件设置为水滴速度与自由来流空气速度相同，且水滴体积因子设置为自由来流液态水含量与水滴密度的比值.对于物面边界条件而言，若水滴速度方向指向物面内，则意味着水滴与物面发生撞击，同时水滴速度取为物面第一层网格的流场值；若水滴速度方向指向物面外，则意味着水滴由物面内向物面外流动，这种现象在真实情况下不可能发生，此时将水滴体积因子取为一个接近于0的较小量.
结冰表面局部水滴收集系数指结冰表面的法向水滴流量.应用欧拉方法中的水滴控制方程计算水滴收集系数，且计算公式为
式中，n为物面法向向量.在实际计算中可以使用网格中距结冰表面第一内点的水滴速度代替 u.
首先进行发动机进气道三维通气模型和进排气模型的对比计算，从而分析发动机进气道进出口边界设置对唇口水滴撞击特性的影响.进气道三维通气模型网格示意图如图1所示.
图1 进气道三维通气模型网格示意图
Fig.1 Schematic grids in 3D ventilation model for air inlet
在通气模型计算中速度为85 m／s，迎角为0°，海拔高度为0 km.通过流场计算结果可以得到流过进气道的空气质量流量为124 kg／s.因而设置进排气模型进口边界条件的进气流量亦为124 kg／s，并与通气模型的计算结果进行对比.
为了方便进行数据对比，沿发动机进气道周向设置角度 θ，示意图如图2所示.由图2可见，从发动机正面看，12点钟位置为0°，按顺时针方向依次为0°、45°、90°、135°和180°.由于进气道模型左右对称，对比分析时只截取0°、90°和180°位置处的数据即可.
图2 周向角度坐标设置示意图
Fig.2 Schematic position setting of circumference angle
发动机通气模型和进排气模型在不同角度位置截面处的水滴收集系数分布如图3所示.由图3可见，两种模型的水滴收集系数分布相差很小，表明进排气模型进出口边界位置的设置是合理的，其对唇口水滴撞击特性的影响很小.
为了分析不同进气流量对发动机唇口水滴撞击特性的影响，需要设置相应的模型参数，具体结果如表1所示.
当迎角为0°、水滴中值直径为20μm时，各进气流量状态下进排气模型的水滴收集系数分布如图4～6所示.图4～6中横轴为发动机吊舱进气道前缘（唇口）弧长，其0点表示前缘点，正方向指向进气道外侧方向，负方向指向进气道内侧方向.
由图4～6可见，进气流量越大，驻点越向进气道外侧移动，水滴撞击极限也越向外侧移动.同时进气流量越大代表进气道流速越大，水滴收集系数最大值也越大.因此，进排气对进气道唇口的水滴撞击特性影响很大，在确定进气道结冰防护范围时必须考虑进排气的影响.
当迎角为0°、水滴中值直径为30μm时，各进气流量状态下进排气模型的水滴收集系数分布如图7～9所示.观察图7～9可以得到与水滴中值直径为20μm时相同的结论.
图3 进气道进出口边界位置的设置对水滴收集系数的影响
Fig.3 Effect of boundary position setting of inlet and outlet of air inlet on water droplet collection coefficient
表1 三维进排气模型参数
Tab.1 Parameters for 3D air suction and exhaust model
本文采用欧拉方法建立水滴运动方程，基于分区多块结构网格体系并采用有限体积方法求解水滴运动轨迹，实现了三维复杂外形飞机表面可压缩流场水滴撞击特性的计算.本文对发动机进气道唇口进行三维水滴撞击特性计算分析，研究了不同发动机进气流量对进气道唇口水滴撞击特性的影响.结果表明，发动机进气流量的大小对进气道前缘的水滴撞击特性影响很大.发动机进气流量越大，水滴收集系数最大值越大，且水滴撞击极限点越向进气道外侧移动，因此，在航空发动机吊舱进气道前缘的防除冰系统设计中必须考虑不同飞行状态下进气流量的影响.
图4 0°截面处水滴收集系数分布（d＝20μm）
Fig.4 Distribution of water droplet collection coefficient at 0°section（d＝20μm）
图5 90°截面处水滴收集系数分布（d＝20μm）
Fig.5 Distribution of water droplet collection coefficient at 90°section（d＝20μm）
图6 180°截面处水滴收集系数分布（d＝20μm）
Fig.6 Distribution of water d rop let collection coefficient at 180°section（d＝20μm）
图7 0°截面处水滴收集系数分布（d＝30μm）
Fig.7 Distribution of water droplet collection coefficient at 0°section（d＝30μm）
图8 90°截面处水滴收集系数分布（d＝30μm）
Fig.8 Distribution of water drop let collection coefficient at 90°section（d＝30μm）
图9 180°截面处水滴收集系数分布（d＝30μm）
Fig.9 Distribution of water drop let collection coefficient at 180°section（d＝30μm）
参考文献：（
References：）：
［1］Langmuir I.Mathematical investigation of water droplet trajectories［J］.Journal of Atmospheric and Oceanic Technology，1946（5）：160－170.
［2］Gelder T F.Droplet impingement and ingestion by supersonic nose inlet in subsonic tunnel conditions［R］. NACA Technical Notes，1958.
［3］Brun R J.Cloud-droplet ingestion in engine inlets with inlet velocity ratios of1.0 and 0.7［R］.NACA Technical Notes，1956.
［4］Kim JJ.Computational particle trajectory analysis on a 3-dimensional engine inlet［R］.AIAA Paper，1985.
［5］WrightW B.Validation results for LEW ICE 3.0［R］.NASA CR，2005.
［6］WrightW B.Users manual for the improved NASA Lew is ice accretion code LEW ICE 1.6［R］.NASA CR，1995.
［7］Bidwell C S.Collection efficiency and ice accretion calculations for a Boeing 737-300 inlet［R］.SAE International，1996.
［8］顾海君.发动机进气道唇口水滴撞击特性的数值模拟研究［D］.南京：南京航空航天大学，2005.（GU Hai-jun.Numerical simulation research of water droplet impingement permormace of engine inlet［D］.Nanjing：Nanjing University of Aeronautics and Astronautics，2005.）
［9］杨倩，常士楠，袁修干.发动机进气道水滴撞击特性分析［J］.北京航空航天大学学报，2002，28（3）：362－365.（YANG Qian，CHANG Shi-nan，YUAN Xiu-gan. Analysis on droplet trajectories of an engine inlet［J］. Journal of Beijing University of Aeronautics and Astronautics，2002，28（3）：362－365.）
［10］申晓斌，林贵平，杨胜华.三维发动机进气道水滴撞击特性分析［J］.北京航空航天大学学报，2011，37（1）：1－5.（SHEN Xiao-bin，LIN Gui-ping，YANG Sheng-hua. Analysis on three dimensional water droplets impingement characteristics of engine inlet［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37（1）：1－5.）
［11］李艳亮，董军，杜希奇.零质量射流对翼型增升影响的数值模拟研究［J］.航空计算技术，2008，38（5）：41－45.（LIYan-liang，DONG Jun，DU Xi-qi.Computation of zero-mass synthetic jets for improving airfoil lift［J］. Aeronautical Computing Technique，2008，38（5）：41－45.）
（责任编辑：尹淑英 英文审校：尹淑英）
Effect of air inflow of engine on impingement characteristics of water drop lets at leading edge
YAO Hong，NING Yi-jun
（Key Laboratory of Aviation Science and Technology on Aerodynamics of High Speed and High Reynolds Number，AVIC Aerodynamics Research Institute，Shenyang 110034，China）
Abstract：In order to understand the effect of air inflow of engine on the impingement characteristics of water droplets，the motion equations for the water droplets were established with Eulerian method.Based on the multi-block structured grid system in subregion，the motion trajectory of water droplets was solved.A calculation method for the impingement characteristics of water droplets on the airplane with complex configuration was proposed.Based on the proposed method，the analysis on the impingement characteristics of water droplets for the air inlet lip of engine was carried out.In addition，the effect of different air inflow of engine on the impingement characteristics of water droplets at the air inlet lip was investigated.The results show that the bigger the air inflow of engine is，the bigger the maximum value of water droplet collection coefficient is.Therefore，the effect of air inflow at different flying states must be considered in the design for the anti-icing system at the leading edge of nacelle air inlet of aero-engine.
Key words：engine；nacelle；Eulerian method；supercooled water droplet；impingement characteristic；anti-icing system；structural grid；numerical simulation
中图分类号：TM 303
文献标志码：A
文章编号：1000－1646（2016）06－0645－06
doi：10.7688／j.issn.1000－1646.2016.06.09
收稿日期：2016－07－10.
基金项目：中航工业集团公司创新基金资助项目（2013A62601R）.
作者简介：幺 虹（1963－），女，辽宁本溪人，高级工程师，主要从事飞行器非定常空气动力学等方面的研究.
网络出版地址：11－07 12∶32在中国知网优先数字出版.
网络出版地址：http：∥www.cnki.net／kcms／detail／21.1189.T. 20161107.1232.018.htm l