-
0 引言
-
1994年美国发生ATR-72飞机事故，2006年我国安徽省发生军用飞机失事，据调查发现，事故均为机翼上形成的冰脊导致的。形成冰脊的原因一般是在结冰气象条件下，机翼前缘热气防冰能力不足以完全蒸发撞击在前缘上的水滴，未蒸发的部分在空气吹拂下溢流至防冰区域外冻结成冰，即溢流结冰，并逐渐累积形成冰脊。从20世纪90年代开始，随着大家对冰脊的认识逐渐加深，发现冰脊会严重破坏飞机的气动外形，使绕过机翼的气流性质变化，对飞机性能造成较大影响，学者开始对冰脊的危害展开深入研究。Bragg[1]等研究者通过一系列试验和计算，模拟冰脊形状，开始研究冰脊对飞机气动性能的影响，分析了不同冰脊条件下几种翼型的气动特性，并经过进一步研究分析，指出在特定的结冰条件下，位于翼型的主升力区会有脊状冰生长聚集，导致翼型的最大升力系数下降。Broeren[2]等通过试验研究对比了多个不同冰形对翼型的气动特性影响，结果表明在翼型的最大升力系数、失速迎角方面，脊状冰比流向冰和角状冰的危害更大。李焱鑫[3]等人对生成了溢流冰脊的超临界翼型进行了数值计算分析，结合雷诺应力模型，研究了冰脊对超临界翼型的气动性能影响，结果表明溢流冰脊造成流场提前发生气动分离，气动性能严重下降。
-
除了研究冰脊的危害，还有不少学者研究了冰脊的形成机理及影响因素。对于过冷大水滴条件下冰脊的形成，Miller[4]等通过实验进行过仔细的研究，得到了冰脊的位置和高度受不同攻角、空速、温度、撞击水滴大小、襟翼位置的影响结果。Tan[5]等通过风洞试验对过冷大水滴的机理进行了研究，包括碰撞分析、飞溅和分裂特性，为后续研究冰脊的形成机理及其对气动特性的影响奠定了基础。Honsek[6]等通过数值模拟计算，采用欧拉法研究了过冷大水滴的碰撞特性，模拟计算得到的水滴收集率与试验结果一致。易贤[7]采用CFD方法，先获得空气流场分布，再采用水滴项控制方程的有限体积求解方法，获得水滴收集率，并在NACA0012翼型上进行计算并得到了验证。顾洪宇[8]进行了机翼热气防冰的数值模拟，对内外流场及蒙皮进行了热质耦合传热计算，得到模拟溢流冰脊，并研究了不同结冰环境温度和结冰时间下形成的冰脊对机翼气动特性的影响。
-
现役主流商用飞机多采用热气防冰系统来降低结冰所带来的危害，即通过在机翼前缘的缝翼内部装备有热气防冰系统对机翼防护区进行加热，避免机翼关键部位严重结冰。对于热气防冰构型，在外界结冰环境一定时，供气参数是影响防冰系统效果的主要参数。现有文献中对冰脊的研究更多的是关注影响冰脊形成的外部环境因素以及冰脊形成后对气动性能的影响，而对供气参数这个内部环境因素的影响研究较少。针对热气防冰构型的机翼，本文基于流固耦合传热思想，采用欧拉法进行水滴撞击计算，并将计算结果与文献的计算结果进行对比，验证本文算法的正确性。在得到内外流场、水滴撞击计算结果后，进行内外流场及蒙皮热质耦合传热计算，研究了供气温度、热防护区域对机翼表面溢流结冰特性的影响规律，并分析了产生的溢流冰对翼面气动性能的影响，对后续相关研究提供参考和借鉴。
-
1 数值模拟方法
-
本文采用的热气防冰系统的模型视图如图1所示，模型弦长为1 m，沿展向宽度为20 mm。
-
图1 热气防冰系统模型（单位:mm）
-
图1中参数含义如下:1 蒙皮与外流场交界面; 2 蒙皮; 3 蒙皮与内流场交界面; 4 内流入口; 5 供气笛形管; 6 内流出口; H1 热防护区域; H2 蒙皮厚度（1 mm）; H3 内流入口与前缘距离（25 mm）; D 笛形管直径（50 mm）; d1 内流入口直径（3 mm）; d2 内流出口直径（7 mm）。根据文献[9]，笛形管喷口位置对热气防冰结果影响甚微，所以为了便于计算，直接将喷口定义在-x方向。该模型并非最优，但是具有一定代表性。
-
机翼热气防冰数值模拟计算的主要工作包括网格划分、流场计算、水滴撞击特性计算、冰积聚和后流水计算、蒙皮传热耦合计算等几个部分。
-
1.1 计算网格生成
-
机翼热气防冰性能计算所需的计算网格包括:机翼外流场计算网格、防冰腔内流场计算网格和蒙皮网格，如图2所示。
-
图2 网格分布
-
外流场和固体蒙皮为结构网格，由六面体结构单元构成，单元数分别为963 500和200 000，内流场网格为非结构网格，由四面体和棱柱单元构成，单元数为1 826 223。
-
1.2 流场计算
-
自然界流体流动是有章可循的，任何流动系统都是遵守着某种规律在运行，而研究这种规律的流场控制方程，实则是在三大基本守恒定律，即质量守恒定律、动量守恒定律和能量守恒定律的基础上构建的微分方程组。
-
1）质量守恒方程:
-
式中，ρ表示密度（kg/m3），u表示速度（m/s），t为时间（s），下标air表示空气解; 根据文献[10]，空气的可压缩性会对水滴撞击特性产生影响，但是在马赫数不超过0.6时，这种影响较小，因此本文忽略空气的可压缩性，假定密度不随时间发生变化，质量守恒方程变为:
-
2）动量守恒方程:
-
描述黏性流体流动的动量守恒方程也称为N-S方程，数学表达式为:
-
式中，σ表示应力张量:
-
式中，p表示静压，μ为动态黏度。对于无黏流动，动态黏度设置为0，即得到Euler方程。对于层流，黏度由Sutherland公式决定:
-
式中，T表示Kelvin单位的温度，带下标∞的变量T∞=288 K，μ∞=17.9×10-6。
-
3）能量守恒方程:
-
式中，H为总焓，γ为空气比热比（空气为1.4），κ为热撞到系数，计算方式与动态粘度类似:
-
再加入气体状态方程，使得整个方程组封闭:
-
式中，R表示气体状态常数。
-
1.3 水滴撞击特性计算
-
采用欧拉两相流方法模拟水滴撞击现象，物理模型使用纯水滴模型，通常选择单一水滴尺寸分布。由于水滴直径较小，欧拉法的两相流控制方程通常对水滴进行一些必要的假设[11]:1）水滴均匀分布，保持为球形且不变形、不破裂; 2）水滴之间不碰撞、不凝聚，撞到壁面后也不飞溅; 3）水滴和周围空气之间不发生热交换; 4）只考虑空气对水滴的阻力影响。基于以上假设，控制方程如下:
-
1）连续方程:
-
2）动量方程:
-
式中，αw为水滴在控制容积中所占的体积比，计算式为:αw=Vw/V，其中Vw是水滴体积，V是控制容积。
-
1.4 冰积聚和后流水
-
固体表面的薄液体在空气剪切力与重力合力作用后向后流动，基于表面温度，部分薄层液体会结冰（冰积聚）、蒸发或升华，求解冰积聚一般就是在物面求解两偏微分方程组。第一个方程为质量守恒方程:
-
其中，右端三项分别为水滴撞击导致的质量转移（薄层的源项）、蒸发导致的质量转移、冰积聚导致的质量转移（薄层沉积）。
-
第二个偏微分方程为能量守恒方程:
-
其中，右端前三项分别为过冷水滴撞击导致的热转移、蒸发导致的热转移、冰积聚导致的热转移。右端最后两项为辐射和对流导致的热转移。
-
系数ρf，cf，cs，σ，ks，Levap，Lsubl，Lfusion为流体和固体的物理属性系数; 参考状态T∞，V∞，LWC为气流和水滴参数; τa，wall为当地物面剪切应力，Qh为对流热通量; β为水滴收集系数，Vd为水滴撞击速度; mevap为蒸发质量通量; ms为给定时间段内的冰积聚总质量。
-
薄层厚度hf、空气/水薄层/物面交界的平衡温度Tf和瞬时冰的积聚质量mice，由兼容性关系来封闭方程系统，作出如下约束:
-
上述不等式保证了平衡温度在冰点0℃以下时，模型预测没有液态水存在，平衡温度在冰点以上时没有冰形成。
-
1.5 蒙皮传热耦合计算
-
通过多次调用外流场、水滴撞击特性和内流场的计算结果进行耦合传热计算，最终得到飞机防护区蒙皮内/外表面温度分布以及防护区外的表面溢流结冰特性，数值模拟流程图如图3所示。
-
2 验证算例
-
为了验证计算模型和方法的正确性，本文以NACA0012翼型（弦长1 m）为研究对象，进行了水滴撞击特性计算和溢流结冰计算，并将计算结果和公开发表的文献中的数值计算结果进行比较。
-
图3 机翼蒙皮耦合传热模拟流程
-
算例1:水滴撞击特性计算。飞行与气象状态参数:远场速度U∞=139 m/s; 液态水含量LWC∞=0.55 g/m3; 水滴粒径MVD=16 μm; 攻角α=5°，环境压力101 325 Pa，环境温度300 K。
-
计算结果和文献[12]中拉格朗日法的计算结果进行比较。图4为机翼前缘附近压力云图和马赫数云图，图5是液态水含量云图，蓝色的代表水滴遮蔽区，从图中可以看到在壁面附近空气绕过壁面流动，而水滴则在机翼前缘发生撞击，由于攻角的存在，驻点往下表面偏移。由于机翼的几何形状，导致一部分水滴撞击壁面，一部分水滴掠过机翼，并在机翼上下表面形成水滴遮蔽区，在此区域内几乎没有水的存在。由于攻角的存在，上表面水滴遮蔽区范围增大，下表面水滴遮蔽区范围减小。
-
图6是局部水收集系数结果，并且和文献[12]进行对比，结果吻合较好。从算例的计算结果和文献值的比较可知，本文对水滴撞击特性的计算模型及方法正确，特别是局部水收集系数的极值，和文献值很接近，撞击极限存在些许误差，但误差较小。
-
图4 空气压力云图及MA云图
-
图5 液态水含量云图
-
图6 局部水收集系数与文献[12]的比较
-
算例2:溢流结冰计算。飞行与气象状态参数:远场速度U∞=67.05 m/s; 液态水含量LWC∞=1 g/m3; 水滴粒径MVD=20 μm; 攻角α=0°，环境压力101 325 Pa，环境温度263.15 K，供气温度335.4 K，笛形管喷射速度Vin=367.255 7 m/s; 结冰时间分别为360 s（case a）和2 000 s（case b）。
-
溢流结冰计算结果如图7所示，表1是带冰翼型的升阻力数据，与文献[8]中的结果相近。
-
图7 溢流结冰计算结果
-
3 算例分析
-
本章以NACA0012翼型为例，在验证算例的结冰计算条件基础上，调整供气温度值，重新进行防冰腔内流场计算，得到稳定状态下的温度分布和速度分布，然后与外流场和蒙皮进行耦合传热计算，得到稳定状态下蒙皮表面温度分布，再进行表面溢流结冰计算，得到溢流冰形，观察溢流冰的位置及高度，通过对比来研究供气温度的影响。再更改热气防冰区域大小，重新进行内、外流场及蒙皮网格划分，并计算出相同供气温度下内流场温度分布、蒙皮耦合传热平衡温度分布，最后进行表面溢流结冰计算，得到溢流冰形，观察溢流冰的位置及高度，通过对比来研究热气防冰区域的影响。
-
3.1 供气温度对溢流冰形的影响
-
选取5组供气温度（310 K、320 K、330 K、340 K、350 K），先计算出供气温度为330 K时内流场及蒙皮耦合传热计算结果，得到溢流冰形数据，再通过降低供气温度（320 K、310 K）和升高供气温度（340 K、350 K），分别计算得到内流场及蒙皮耦合传热计算结果，得到溢流冰形数据，通过对比来研究供气温度的影响。
-
热防护区域H1=0.1 m，飞行与气象状态参数:远场速度U∞=67.05 m/s; 液态水含量LWC=1 g/m3; 水滴粒径MVD=20 μm; 攻角α=0°，环境压力101 325 Pa，环境温度263.15 K，结冰时间t=2 000 s。
-
1）供气温度降低时的影响
-
防冰腔内供入热气（T=330 K），进行内流场计算，得到温度分布如图8（a）所示; 然后进行内外流场及蒙皮传热耦合计算，得到蒙皮表面平衡温度，如图8（b）所示，在此供气温度条件下，热防护区蒙皮表面最低温度288.6 K，超过273.15 K，不会发生结冰; 再进行蒙皮表面溢流结冰计算，结果如图8（c）所示，在热防护区后方，机翼上下表面均有溢流冰产生，表明此温度无法将撞击水完全蒸发，而是产生了溢流，并在热防护区外冻结形成脊冰。进一步降低供气温度至320 K、310 K，得到溢流结冰计算结果如图9所示，结果表明溢流冰形位置更加靠近热防护区，而且冰形高度越高。
-
图8 供气温度330 K时计算结果
-
图9 溢流结冰计算结果（对比330 K、320 K、310 K）
-
2）供气温度升高时的影响
-
将供气温度升高至340 K、350 K，计算结果如图10所示，随着供气温度的升高，溢流结冰位置更加靠近热防护区，冰形高度也越低，结冰区域减小，进一步验证了供气温度升高导致溢流水量减少，进而对溢流冰形高度及冻结区域范围的影响。
-
图10 溢流结冰计算结果（对比330 K、340 K、350 K）
-
继续提高供气温度至358.5 K，内流场计算结果、蒙皮耦合传热计算结果及溢流结冰计算结果如图11（a）所示，机翼上表面只形成了非常少量溢流冰，而机翼下表面形成的溢流冰明显多于机翼上表面; 提高供气温度至359 K，计算结果如图11（b）所示，机翼上表面无溢流冰产生，表明在此供气温度下，机翼上表面撞击水在热防护区域内被完全蒸发，这是由于在重力作用下，下表面溢流水量比上表面多，在此供气温度下，机翼下表面的溢流水大部分被蒸发，有少量流到了热防护区外并冻结成冰; 进一步提高供气温度至360.5 K，计算结果如图11（c）所示，此时机翼上下表面溢流冰均不存在，表面供气温度达到360.5 K时，机翼上下表面撞击水在热防护区域内均被完全蒸发。
-
图11 溢流结冰计算结果
-
通过以上模拟计算结果可知，供气温度会直接影响溢流结冰能否形成以及冰形的高度和位置，并且即使机翼上表面撞击水被完全蒸发，机翼下表面的溢流水在重力作用下仍会流至热防护区外冻结成冰，对飞行安全仍然存在不良影响，所以在设计机翼热气防冰系统时，必须充分考虑各种因素的影响。
-
3.2 热防护区域尺寸对溢流结冰的影响
-
选取H1=0.08 m，0.09 m，0.11 m，0.12 m，分别模拟计算得到对应的溢流冰形，如图12所示。
-
图12 溢流结冰计算结果对比
-
机翼上下表面产生溢流结冰后翼型的气动性能参数（升力、阻力）计算结果如表2所示。
-
对比发现，随着热防护区域增大，机翼上下表面产生溢流结冰的高度逐渐减小，对升阻力影响逐渐降低。但是，热防护区域越大，为保证整个防冰区域蒙皮温度超过冰点，所需要的供气流量就要越高，对发动机所产生的损耗也就越大，而且机翼内部需要安装大量系统件及留有足够大的油箱尺寸，所以热防护区域并不是越大越好，而是需要综合考虑发动机性能及机翼内部空间布局要求，选定一个合适的热气防冰尺寸，既满足防冰要求，又不至于对发动机要求过高，同时机翼内部又能留出合适的空间。
-
4 结论
-
本文通过机翼流场计算、水滴撞击特性计算、冰积聚和后流水计算、蒙皮传热耦合计算、溢流结冰计算等，得到不同状态下的溢流冰形数据，分析了供气温度、热防护区域对翼面溢流冰形成的影响，计算结果表明:供气温度直接正向影响热交换后蒙皮表面温度，在供气温度低于330 K时，供气温度越低，撞击水蒸发量越少，溢流水量越多，而且溢流水温度低，冻结位置靠近热防护区，冰形高度及结冰区域范围明显增加; 供气温度高于330 K时，随着供气温度增加，溢流水量减少，与外流场热交换加快，在距离热防护区较近位置处开始冻结，溢流冰高度及结冰区域范围明显降低。在飞行状态、结冰气象条件及供气温度不变的情况下，热防护区域范围对溢流结冰结果也会产生影响，热防护区域增加时，参与热交换的蒙皮范围增加，导致蒙皮表面平衡温度降低，但是仍然超过273.15 K（0℃），会提高撞击水蒸发量，虽然溢流冰冻结位置更加靠近热防护区，结冰区域范围增加，但是相对整个翼型而言，冻结位置远离了驻点，而且冰形高度降低，通过对比升力和阻力发现，对气动特性影响减弱。
-
参考文献
-
[1] BRAGG M B.Aircraft aerodynamic effects due to large droplet ice accretions[C].34th Aerospace Sciences Meeting and Exhibit,January 15-18,1996.[S.l.:s.n.],1996.
-
[2] BROEREN A P,BRAGG M B.Effect of airfoil geometry on performance with simulated inter-cycle ice accretions[J].Journal of Aircraft,2005,42(1):121.
-
[3] 李焱鑫,张辰,刘洪,等.大粒径过冷水溢流结冰的翼型气动影响分析[J].空气动力学学报,2014,32(3):376-382.
-
[4] MILLER D R,ADDY JR H E,IDE R F.A study of large droplet ice accretions in the NASA-Lewis IRT at near-freezing conditions[C].34th Aerospace Sciences Meeting and Exhibit,January 15-18,1996.[S.l.:s.n.],1996.
-
[5] TAN S C,PAPADAKIS M,MILLER D,et al.Experimental study of large droplet splashing and breakup [C].45th AIAA Aerospace Sciences Meeting and Exhibit,January 08-11,2007.[S.l.:s.n.],2007.
-
[6] HONSEKR,HABASHI W G,AUBE M S.Eulerian modeling of in-flight icing due to super-cooled large droplets [J].Journal of Aircraft,2008,45(4):1290-1296.
-
[7] 易贤.飞机积冰的数值计算与积冰试验相似准则研究[D].绵阳:中国空气动力研究与发展中心,2007.
-
[8] 顾洪宇,桑为民,庞润,等.机翼热气防冰及冰脊形成数值模拟[J].气体物理,2019,4(4):41-49.
-
[9] 卜雪琴,林贵平,郁嘉.机翼电加热防冰表面内外传热的耦合计算[J].航空动力学报,2010,25(7):1491-1496.
-
[10] 陈维建.飞机机翼结冰的数值模拟研究[D].南京:南京航空航天大学,2007.
-
[11] 吴孟龙,常士楠,冷梦尧,等.基于欧拉法模拟旋转帽罩水滴撞击特性[J].北京航空航天大学学报,2014,40(9):1263-1267.
-
[12] WIROGO S,SRIRAMBHATLA S.An Eulerian method to calculate the collection efficiency on two and three dimensional bodies[C].41st Aerospace Sciences Meeting and Exhibit,January 06-09,2003.[S.l.:s.n.],2003.
-
摘要
机翼热气防冰数值模拟中,根据N-S方程求解流场,用Euler法获得水滴撞击特性,通过蒙皮导热将求解得到的内外流场进行耦合传热并达到稳定后,开始模拟结冰,来进行机翼热气防冰及形成溢流结冰的数值计算。计算结果表明,热空气防冰数值模拟是可行和合理的。采用数值模拟方法对机翼热空气防冰过程进行了模拟,得到了由于引气温度不足或机翼热空气保护面积不同而导致的不同溢冰高度和位置。分析了供气温度、防冰区域对翼面溢流冰形成的影响。结果表明:供气温度直接正向影响热交换后蒙皮表面温度。供气温度越低,溢流冰形高度越高,对气动特性影响也越大;热防护区域范围对溢流结冰结果也会产生影响,热防护区域越大,冻结位置距离驻点越远,而且冰形高度越低,对气动特性影响也越弱。
Abstract
In the numerical simulation process of wing hot air anti-icing, the flow field is solved according to the Navier Stokes equation, and the water droplet impact characteristics are solved by the Euler method, and the solved internal and external flow fields are coupled with heat transfer through the skin heat conduction to achieve stability, and then the simulation of icing is started to carry out the numerical calculation of wing hot air anti-icing and overflow ice formation. The calculation results show that the numerical simulation of hot air anti-icing is feasible and reasonable. In this paper, the numerical simulation method was used to simulate the anti-icing process of wing hot air, and the different ice overflow height and locations caused by insufficient bleeding air temperature or different wing hot air protection area were obtained. The influence of the bleeding air temperature and the anti-icing area on the formation of wing overflow ice was analyzed. The results show that the bleeding air temperature directly and positively affects the skin surface temperature after heat exchange. The lower the bleeding air temperature is, the higher the overflow ice is, and the greater the influence on the lift is. The scope of the hot air protection area will also affect the overflow icing results. The larger the hot air protection area is, the farther the overflow ice is, and the lower the ice height is, the weaker the impact on the lift is.
