2. 航空科学技术国家实验室, 北京 100191;
3. 中航工业第一飞机设计研究院, 陕西 西安 710089;
4. 上海飞机设计研究院 环控氧气部, 上海 201203
2. National Laboratory for Aeronautics and Astronautics, Beijing 100191, China;
3. AVIC The First Aircraft Institute, Xi'an 710089, China;
4. Environment Control System Department, Shanghai Aircraft Design and Research Institute of COMAC, Shanghai 201203, China
飞机结冰是危害飞行安全、导致飞行事故的重要因素，机翼和发动机短舱结冰对飞行安全的危害尤其严重。飞机结冰不仅会直接破坏飞机的气动外形，导致升力下降。对于发动机而言，结冰还可引起发动机喘振，甚至冰块脱落，损坏发动机叶片，严重时甚至会导致发动机熄火。可见对发动机短舱进行防除冰设计的意义非常重大。目前飞机机翼与发动机多采用热气防冰系统，需要从发动机引气，但从发动机引气的同时也会给发动机带来一定的代偿损失。因此必须在有限的发动机热气供给情况下保证防冰性能。
对飞机结冰的研究通常主要有三种：结冰飞行试验、结冰风洞试验、数值模拟。其中前两种试验方式需耗费较大的人力物力，且由于条件限制，无法覆盖所有的结冰条件。近年来，随着计算机硬件成本的不断降低，数值模拟方式在经济成本上较其他两者有了巨大的优势，所以数值模拟方法逐渐成为了研究飞机结冰/防冰问题的重要手段[1-6]。目前国内展开的防冰性能研究多是针对飞机机翼[7-10]，针对发动机短舱的较少[11-13]。本文针对某飞机发动机短舱在给定条件下的防冰性能进行数值分析计算，采用了欧拉法计算水滴撞击特性，利用距离反比插值法进行内外流场数据交互，得到了短舱表面温度和溢流水结果，并分析了发动机进气流量对表面温度的影响。
1 计算方法 1.1 三维防冰表面温度计算方法影响热气防冰表面温度的因素有外部空气流动与换热、腔内部热气流动与换热、蒙皮导热，因此表面温度是多场耦合的结果。计算防冰表面温度时需要考虑多个计算区域：外部空气-水滴两相流场区域、发动机短舱蒙皮外表面上的水膜、蒙皮、防冰腔内部，计算时需要对这4个计算区域进行耦合，最终得到收敛后的发动机短舱表面温度。具体流程包括：计算外部整机空气流场、在空气流场的基础上计算水滴场、计算防冰腔内部热气流场、对蒙皮内外热流进行耦合计算[8]。
1.2 外部流场及水滴场计算选取半整机模型划分网格，采用了混合型网格。大部分区域采用非结构化六面体网格绘制，仅在发动机附近流场区域使用结构化网格绘制并对唇口附近进行加密，以保证短舱附近计算结果的准确性与精确度。然后将划分好的网格导入Fluent软件中进行外部流场的计算。
空气流场计算时，由于飞行速度较大且发动机周围流场复杂，计算时采用Fluent内置的三维粘性可压流求解器。计算时使用如下的边界条件：外部流场边界使用远场边界条件，并将发动机入口段某截面处设为压力出口，用匹配流量的方式模拟发动机的抽吸效应对流场的影响。发动机外壁面采用恒壁温边界条件。计算收敛后得到空气速度和对壁面的对流换热热流量。
|
| 图 1 计算网格 Fig. 1 Computational grid |
水滴场使用Fluent中的用户自定义标量UDS(User Defined Scalars)进行计算。分别设置4个UDS：水滴x、y、z方向的速度、水滴相对容积分数，并建立水滴场的连续性和动量方程，使用Fluent自带求解器进行求解。计算时使用如下的边界条件：流场入口处认为水滴速度与空气速度一致并且水滴相对容积分数为1，认为水滴可以自由通过内、外流场出口；当水滴速度指向壁面时，认为水滴可以自由通过，反之指向流场时，说明此处为水滴遮蔽区，局部水收集系数β=0。
使用上述方法计算得到水滴场相应参数后，由定义可知欧拉法下β的计算公式为：
(1) 式中，u为水滴速度方向；n为发动机短舱壁面的单位法向矢量；V∞为自由流水滴速度；αn为水滴的相对容积分数，如下式：
(2) 式中，α∞是自由流中水滴的容积分数。
1.3 防冰腔内部流场计算选取计算的是发动机短舱笛形管射流防冰系统。为计算防冰表面温度分布，需要先对防冰腔内部的热气流动情况进行计算。使用的计算模型为整个发动机短舱环状防冰腔。射流孔均匀分布于整个笛型管上，热气从射流孔流出加热蒙皮，然后在腔内流动，并最终从排气孔流出至外部环境。网格分为蒙皮导热区域及内部流动区域，内部流场网格总数965万，其中固体区域77万，流体区域888万，全部使用结构化网格绘制。计算时使用如下的边界条件：笛形管射流孔出口为计算区域的入口，设为压力入口条件；排气孔为计算区域的出口，设为压力出口条件；防冰腔外壁面假设为均匀温度293.15K，作为表面温度迭代计算时的初始温度。
|
| 图 2 发动机短舱外部及内部构造 Fig. 2 External and internal structure of engine nacelle |
1.4 内外流场数据交互
由于是通过外部网格计算来获取短舱表面的局部水收集系数、热流、剪切力等数据，而计算防冰腔表面温度时则是使用内部网格计算，所以需要将外部网格计算时获得的结果数据导入至内部网格相应的网格点上。使用距离加权反比插值法[14]来进行插值计算，计算公式如下：
(3) 耦合计算蒙皮传热过程时需同时考虑蒙皮导热及内部流动换热两个方面，本文使用强固耦合方法[8]对其进行迭代求解，具体求解过程如下：
1) 使用CFD方法计算得到蒙皮外表面换热系数；
2) 对蒙皮外表面加载温度边界条件，使用UDF计算蒙皮表面水蒸发量、溢流水量及热载荷；
3) 引入松弛系数ω1，将松弛处理后的热载荷作为边界条件加载至蒙皮外表面，使用Fluent计算蒙皮导热及内部流动换热，并得到蒙皮外表面温度；
(4) 4) 引入松弛系数ω2，对此温度进行松弛处理，重新作为边界条件加载至蒙皮外表面；
(5) 5) 重复步骤1)~步骤4)，直至热载荷以及表面温度计算结果收敛。
本文计算发动机短舱防冰性能情况，不考虑结冰过程，此时蒙皮热载荷主要由如下4项构成：
1) 蒙皮与外部气流之间的对流换热；
2) 蒙皮表面水蒸发所带走的热量；
3) 外部流场中的水撞击至蒙皮表面所带来的能量；
4) 水流入/流出计算区域所带来/带走的能量；
蒙皮总热载荷计算公式如下：
式中:mimp、mevap、min、mout分别为撞击水质量流量、蒸发水质量流量、流进控制体的溢流水质量流量、流出控制体的溢流水质量流量，cp,w为水的比热容，ilv为汽化潜热，Tref为参考温度273.15K，Trec为恢复温度，定义如下：
(7) 式中:rc为恢复系数，Te为附面层外边界处温度。
在计算表面换热系数时采用CFD方法，并且在内、外网格数据交互时，使用距离加权反比插值算法来保证数据准确性。由于此方法在计算时是一个纯三维过程，相比以前采用附面层积分方法计算n个二维截面(外部流场的)的对流换热系数，然后插值到三维内部流场计算网格的外表面上，本文的方法更适合于三维问题。
2 算例结果与分析 2.1 计算状态说明使用上述发动机短舱热气防冰系统作为计算对象，进行三维内外耦合计算，获得了局部水收集系数、蒙皮外表面温度、表面溢流水等结果。
算例已知飞行气象条件：飞行马赫数0.36，液态水含量2.3g/m3，攻角7°，环境温度266K，压力75 265Pa，水滴直径20μm。防冰系统引气状态：热气喷口总压211 313Pa，总温493K。
2.2 计算结果以下各曲线图均以蒙皮外表面处的弧长作为x轴，相应的计算结果值作为y轴。s=0为几何驻点处，s为负数代表发动机短舱内侧表面，反之s为正数代表发动机短舱外侧表面。为更好地对计算结果进行展示，截取了5个切面分别进行对比分析，取法如图 3所示。
|
| 图 3 发动机短舱表面局部水收集系数云图 Fig. 3 Contour of the local water collection coefficient on engine nacelle |
插值计算结果中的局部水收集系数β值云图如图 3所示。从图 3中可知，撞击水量最大的区域仍位于几何驻点附近。
图 4、图 5为内外耦合计算收敛后的蒙皮外表面温度结果。从图中可以看出：
|
| 图 4 发动机短舱表面温度云图 Fig. 4 Temperature contour of engine nacelle surface |
|
| 图 5 发动机短舱表面温度曲线图 Fig. 5 Temperature curves of engine nacelle surface |
1) 由于笛形管射流孔均朝向短舱内侧表面，从图中可以明显看出，发动机短舱内侧表面温度较外侧高 。
2) 笛形管射流孔射流正对位置处的蒙皮表面明显偏高，但远离此区域后快速衰减，并达到周围平均温度。
3) 0°、90°、225°、270°截面位于一列单孔和一列双孔之间，90°和270°截面更靠近单孔，0°和225°截面更靠近双孔。所以0°、90°、225°、270°截面对应的温度分布都出现了三峰，且90°和270°截面中间峰更显著，而0°和225°截面对应的温度分布两侧峰较中间峰显著。180°截面位于两列双孔之间，所以此截面处的温度分布呈现了双峰分布形式。
4) 发动机短舱整体温度基本高于0℃，作为重点关注区域的内侧表面温度较高，且高于10℃，可满足防冰系统需求。
图 6为内外耦合计算收敛后的表面溢流水量结果。从图中可以看出：
|
| 图 6 发动机短舱表面溢流水流量云图(单位：g/s) Fig. 6 Contour of the runback water mass flux on engine nacelle surface (unit:g/s) |
1) 溢流水从驻点位置开始向短舱内外两侧流动，内外表面溢流水均未在防冰区内全部蒸发，会流向防冰区后部且可能形成冰脊；内表面溢流水量明显少于外表面。
2) 笛形管喷口正对的短舱周向位置，沿其弦向方向溢流水明显减少甚至全部蒸发。
3) 在270°截面附近溢流水分布异常，可能是此处短舱与机身的连接件影响了附近流场，从而导致空气与壁面之间剪切力方向发生改变，使得溢流水流动方向改变。
3 发动机进气流量对表面温度的影响液态水含量、环境温度、飞行攻角以及飞行速度等对短舱防冰表面温度的影响,与它们对飞机机翼的影响类似，本文不做赘述。这里主要考查发动机进气流量对短舱表面温度的影响。选取了如表 1状态点进行计算和对比分析，所有状态下防冰系统热气喷口总压205 214Pa，防冰系统热气喷口总温493K。
| Case | Temp/℃ | Mach | Engine Flow Rate/(kg·s-1) | α/(°) | LWC/(g·m-3) | MVD/μm |
| 1 | -2.8 | 0.21 | 106.64 | 7.34 | 2.41 | 20 |
| 2 | -2.8 | 0.21 | 209.11 | 7.34 | 2.41 | 20 |
发动机进气流量变化时，首先会影响空气流场，进而在此基础上影响水滴流场。进气道流量越大，发动机亦能更多地吸取周围的空气，从而使驻点向唇口外侧(+s方向)延伸[15]。从图 7可见，发动机进气流量增大时，驻点向发动机外侧表面偏移。从0°到180°截面，局部水收集系数数值相比小流量时的增大量越来越小，甚至在180°截面时局部水收集系数反而减小。这是两个影响因素作用的结果：一方面，进气流量增大会导致更多的水滴撞击；另一方面180°截面处的发动机短舱比较厚，且撞击驻点离几何驻点(s=0)最远，使得局部水收集系数减小。两者综合作用导致在180°截面附近，进气流量更大时局部水收集系数却减小。图 8是不同发动机引气流量下的温度结果。随着发动机进气流量增大，更多气流流入发动机进气道，
|
| 图 7 发动机二维截面处的局部水收集系数分布 Fig. 7 Distribution of local water collection coefficient of engine 2D cross section |
|
| 图 8 发动机二维截面处的表面温度分布 Fig. 8 Distribution of surface temperature of engine 2D cross section |
发动机内侧壁面对流换热系数增大，而外侧壁面对流换热系数却有所降低，使得发动机短舱表面温度大部分区域降低，但在外侧表面部分位置的温度反而略微升高。另外从图 9中所示的发动机短舱表面溢流水量图中,也可以明显地看出驻点的偏移。溢流水在驻点附近会出现谷值，然后沿气流方向流动;随着撞击水的加入，溢流水量逐渐增多直到撞击极限附近达到最大;然后随着蒸发的作用，溢流水量又逐渐减小。90°截面处溢流水分布曲线出现跳跃，原因是溢流水流动方向和此截面不在一个平面。由于撞击水量增大，表面温度降低，因此发动机表面溢流水总体呈上升趋势。状态2流量比较大，发动机内部抽吸空气能力增强，因此朝发动机内部的剪切力较大，导致某些截面处内表面的溢流水量较大。
|
| 图 9 发动机二维截面处的溢流水量分布 Fig. 9 Distribution of runback water mass flux of engine 2D cross section |
4 结 论
对三维发动机短舱的防冰性能进行了计算和研究。利用CFD方法计算表面换热系数,以及距离加权反比插值法对内外表面数据进行插值，更适用于三维防冰问题的计算。当其他飞行气象条件固定，只考虑发动机进气流量的影响时，随着进气流量的增大，发动机短舱蒙皮表面溢流水量增加，蒙皮表面温度有明显的下降趋势。本文的研究方法可为短舱防冰系统性能评估提供帮助。
| [1] | Messinger B L. Equilibrium temperature of an unheated icing surface as a function of air speed[J].Journal of the Aeronautical Sciences, 1953, 20(1):29–421. |
| [2] | Wright W B. User manual for the NASA Glenn ice accretion code LEWICE version 2.2.2[R]. NASA CR 2002-211793. |
| [3] | Saeed F, Morency F, Paraschivoiu I. Numerical simulation of a hot-air anti-icing system[R]. AIAA 2000-0630. |
| [4] | Honsek R, Habashi W G. FENSAP-ICE: Eulerian modeling of droplet impingement in the SLD regime of aircraft icing[R]. AIAA 2006-465. |
| [5] | Bu X Q, Lin G P, Yu J, et al. Numerical simulation of an airfoil electrothermal anti-icing system[J].Proc. IMechE Part G: J. Aerospace Engineering, 2013, 227(10):1608–1622. |
| [6] | Shen X B, Lin G P, Yu J, et al. Three-dimensional numerical simulation of ice accretion at the engine inlet[J].Journal of Aircraft, 2013, 50(2):635–642.DOI:10.2514/1.C031992. |
| [7] |
Chang S N, Yuan M M, Huo X H, et al. Investigations of the bleed air anti-icing system for an aircraft wing[J].Journal of Aerospace Power, 2008, 23(6):1141–1145. (in Chinese) 常士楠, 袁美名, 霍西恒, 等. 某型飞机机翼防冰系统计算分析[J]. 航空动力学报, 2008, 23(6) : 1141–1145. |
| [8] |
Bu X Q, Lin G P, Yu J. Three dimensional conjugate heat transfer simulation for the surface temperature of wing hot-air anti-icing system[J].Journal of Aerospace Power, 2009, 24(11):2495–2500. (in Chinese) 卜雪琴, 林贵平, 郁嘉. 三维内外热耦合计算热气防冰系统表面温度[J]. 航空动力学报, 2009, 24(11) : 2495–2500. |
| [9] |
Pan X Y. Performance analysis of micro-ejector anti-icing tunnel[D]. Nanjing University of Aeronautics and Astronautics, 2006. (in Chinese) 潘旭云. 某型飞机防冰系统机翼防冰腔性能研究[D]. 南京航空航天大学, 2006. |
| [10] |
Huo X H, Wang D W, Li G P, et al. Research of wing anti-ice system performance validation for civil aircraft[J].Civil Aircraft Design & Research, 2013, 111(4):13–16. (in Chinese) 霍西恒, 王大伟, 李革萍, 等. 某型飞机机翼防冰系统性能验证研究[J]. 民用飞机设计与研究, 2013, 111(4) : 13–16. |
| [11] |
Chang S N, Ai S X. Design and calculation for the anti-icing system of an aircraft engine let[J].Journal of Beijing University of Aeronautics and Astronautics, 2007, 33(6):649–652. (in Chinese) 常士楠, 艾素霄. 飞机发动机进气道防冰系统的设计计算[J]. 北京航空航天大学学报, 2007, 33(6) : 649–652. |
| [12] |
Zhu Y F, Fang Y F, Feng W C. Design and calculation of aircraft nacelle anti-icing system[J].Journal of Aerospace Power, 2012, 27(6):1326–1331. (in Chinese) 朱永峰, 方玉峰, 封文春. 某型飞机发动机短舱防冰系统设计计算[J]. 航空动力学报, 2012, 27(6) : 1326–1331. |
| [13] |
Chen J S. Performance analysis of anti-icer in airplane engine inlet[D]. Nanjing University of Aeronautics and Astronautics, 2008. (in Chinese) 陈景松. 发动机进气道前缘防冰腔性能研究[D]. 南京航空航天大学, 2008. |
| [14] |
Jin G D, Liu Y C, Niu W J. Comparison between inverse distance weighting method and Kriging[J].Journal of Changchun University of Technology, 2003, 24(3):53–57. (in Chinese) 靳国栋, 刘衍聪, 牛文杰. 距离加权反比插值法和克里金插值法的比较[J]. 长春工业大学学报, 2003, 24(3) : 53–57. |
| [15] |
Shen X B, Lin G P, Yang S H. 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. (in Chinese) 申晓斌, 林贵平, 杨胜华. 三维发动机进气道水滴撞击特性分析[J]. 北京航空航天大学学报, 2011, 37(1) : 1–5. |