飞机在含有大量过冷水滴的云层中飞行时,其迎风面会发生结冰现象.发动机进口旋转帽罩处于飞机的迎风面,也会发生结冰现象.旋转帽罩结冰会降低发动机入口气流品质,使发动机性能降低.并且,旋转帽罩表面结冰脱落,被吸入发动机内部,会导致发动机损毁,造成重大飞行事故[1].因此十分有必要对旋转帽罩进行结冰研究.飞机结冰数值模拟研究一般分为空气流场求解、水滴轨迹及撞击特性求解、结冰表面热平衡分析以及结冰冰形计算4个部分.本文主要关注于旋转帽罩水滴轨迹及撞击特性分析.
过冷水滴在空气流场中的运动属于气粒两相流动.按照对水滴颗粒的不同处理方式,主要有两种气粒两相流动研究方法[2]:
1) 把气体当作连续介质,而将颗粒视为离散体系,即拉格朗日法;
2) 把气体与颗粒都看成共同存在且相互渗透的连续介质,把颗粒作为拟流体,即欧拉法.
目前,国内外对于静止部件的水滴撞击特性的求解既有采用拉格朗日法[3],也有采用欧拉法[4,5,6].而对于旋转部件的水滴撞击特性,国内外研究较少,并且多数采用拉格朗日法[7,8].
本文在欧拉法求解静止部件水滴撞击特性的基础上,提出一种旋转帽罩水滴撞击特性数值模拟方法,实现了欧拉法求解旋转帽罩表面的水滴撞击特性.通过对不同转速旋转帽罩水滴撞击特性进行数值模拟研究,发现旋转帽罩转速越大,其表面水滴撞击极限越小.但是,在典型飞行条件与气象条件下,离心力对水滴运动轨迹的影响远小于惯性力对其的影响,旋转帽罩转速对水滴撞击特性影响较小. 1 水滴相数学模型
基于欧拉法在静止部件水滴撞击特性的成功应用[9,10,11,12],本文提出了一种旋转帽罩水滴撞击特性求解方法.本文借助于Fluent软件,采用S-A湍流模型实现旋转帽罩空气流场的求解.而将水滴视为拟流体,采用欧拉法,借助于Fluent UDS模块,实现旋转帽罩水滴流场的求解.
结冰条件下,由于水滴的体积分数一般在10-6量级,可认为空气和水滴之间是单向耦合的,即空气流场影响水滴的运动,而水滴的运动对空气流场没有影响[9,13].这使得空气相和水滴相可以分离求解:先得到空气流场,再计算水滴流场.并且对水滴进行一些合理的假设:①水滴是球形的,无变形或破裂;②水滴之间无碰撞、聚合,撞击到壁面后无飞溅;③水滴和周围空气无质量和热量的传递;④只考虑作用在水滴上的空气阻力. 1.1 控制方程
基于以上假设,稳态水滴相控制方程为
旋转帽罩空气-水滴流场数值模拟中,空气相与水滴相的边界条件处理方式不同.其中,空气相入口采用速度入口边界条件,出口采用出流边界条件,壁面采用动壁面中的旋转、无滑移壁面边界条件.
水滴相的入口采用速度入口边界条件,速度大小与空气速度保持一致;出口则采用出流边界条件.水滴相的壁面边界条件需要进行特殊处理.这是由于水滴撞到壁面后,会积聚到壁面上,从而从流场中离开.因此,在水滴撞击区域,如图 1所示情况a和b,壁面微元单位法向量n与水滴速度ud点积n·ud≤0时,水滴的体积分数与水滴速度保持不变;在非撞击区域,如图 1所示情况c,n·ud>0时,水滴的体积分数取零,水滴速度取下游网格单元的速度值[9,10,12].
 |
| 图 1 水滴相壁面边界条件 Fig. 1 Wall boundary condition of droplet flowfield |
由于水滴相控制方程属于双曲型方程[10],数值迭代过程会出现局部数值奇异,因此常采用一些耗散函数,将局部奇异值耗散,使计算能够得到收敛.本文采用的耗散函数[9]为
局部水滴收集系数β为当地表面微元实际水滴收集率与最大可能水滴收集率的比值,其计算公式为
由于椭圆型帽罩较易结冰[14],因此本文采用椭圆型旋转帽罩作为计算数模,其中长轴半径为0.09m,短轴半径为0.06m.在网格划分过程中,对旋转帽罩尾部进行适当的延伸.在旋转帽罩前缘、尾部以及周围划分10倍远场网格.所划网格如图 2所示.
 |
| 图 2 旋转帽罩计算域网格 Fig. 2 Computational mesh of the spinner |
本文关注重点在于旋转速度对于旋转帽罩水滴撞击特性的影响,因此,取不同旋转速度的飞行及气象条件作为计算条件,分析研究旋转速度对旋转帽罩水滴撞击特性的影响.参照典型民用航空发动机旋转帽罩转速,选取旋转帽罩转速与其他飞行条件如表 1所示.
| 状态 | 旋转速度(r/min) | 压强/Pa | 飞行速度/(m/s) | 温度/K | MVD/μm | LWC/(g/m3) |
| 状态1 | 1000 | 89867 | 78 | 260 | 20 | 1 |
| 状态2 | 3000 | 89867 | 78 | 260 | 20 | 1 |
| 状态3 | 5000 | 89867 | 78 | 260 | 20 | 1 |
| 状态4 | 5000 | 89867 | 170 | 260 | 20 | 1 |
| 状态5 | 5000 | 89867 | 50 | 260 | 20 | 1 |
| 状态6 | 5000 | 89867 | 78 | 260 | 10 | 1 |
| 状态7 | 5000 | 89867 | 78 | 260 | 40 | 1 |
| 状态8 | 0 | 89867 | 50 | 260 | 10 | 1 |
| 状态9 | 5000 | 89867 | 50 | 260 | 10 | 1 |
| 注：MVD—水滴平均容积直径；LWC—液态水含量. | ||||||
利用上述计算方法对表 1中所列计算条件下的旋转帽罩水滴撞击特性进行了数值模拟计算.
水滴在旋转帽罩空气流场中运动,其运动轨迹受到空气的黏性阻力与离心力的影响.由于水滴具有惯性,因此会保持原有运动形式,撞击到旋转帽罩表面.而水滴受到的黏性阻力与离心力则会使水滴偏离原始轨迹,阻碍水滴撞击旋转帽罩.水滴的惯性、受到的黏性阻力与离心力跟水滴的平均容积直径、来流速度以及旋转帽罩旋转速度有关.由于旋转帽罩表面的水滴撞击特性取决于水滴的运动轨迹.因此,本文对不同旋转速度、来流速度以及水滴平均容积直径的计算结果进行了分析.
本文首先取状态3计算结果与成熟的商业结冰计算软件Fensap进行对比.沿旋转帽罩轴向做y=0截面,作局部水滴收集系数沿弦长变化曲线如图 3所示.
 |
| 图 3 状态3，y=0截面局部水滴收集系数分布曲线 Fig. 3 Case 3, cross-section y=0,local water droplets collect coefficient distribution curves |
由图 3可以看出,本文计算结果与Fensap计算结果十分吻合,因此可以证明本文计算方法是合理的.由于本文采用的耗散函数与Fensap所采用的耗散函数不同,导致帽罩后部存在一定差异. 3.1 旋转帽罩转速对水滴撞击特性影响
选取旋转速度不同,其他条件相同的状态1~3的局部水滴撞击特性计算结果进行分析.沿旋转帽罩轴向做y=0截面,作局部水滴收集系数沿弦长变化曲线如图 4所示.由图 4可以看出,各转速下,旋转帽罩局部水滴收集系数分布结果大致相同.随着转速的增加,水滴撞击极限逐渐减小.这与文献[7,15]采用拉格朗日法求解旋转帽罩水滴撞击特性所得结论相同.
 |
| 图 4 状态1~3,y=0截面局部水滴收集系数分布曲线 Fig. 4 Case 1~3,cross-section y=0,local water droplets collect coefficient distribution curves |
图 5所示为状态3的水滴运动轨迹,可以看出水滴偏离基准线,随气流发生旋转,但其偏转角度较小,即水滴受到旋转帽罩的旋转影响较小.
 |
| 图 5 状态3,速度78m/s,转速5000r/min水滴轨迹图 Fig. 5 Case 3,velocity=78m/s,rotate speed= 5000r/min,the scheme of droplet trajectory |
选取来流速度不同,其他条件相同的状态3~5的水滴撞击特性计算结果进行分析.当其他飞行条件与气象条件保持不变,来流速度减小时,水滴受到的惯性力减小.惯性力对水滴运动轨迹的影响减小,相应地黏性阻力与离心力对水滴的运动轨迹影响增大.因此,水滴更易偏离原始轨迹,并且撞击到旋转帽罩表面的趋势减弱,从而得到如图 6所示旋转帽罩表面局部水滴收集系数明显减小,水滴撞击极限范围也明显减小的结果.
 |
| 图 6 状态3~5,y=0截面局部水滴收集系数分布曲线 Fig. 6 Case 3~5,cross-section y=0,local water droplets collect coefficient distribution curves |
选取水滴平均容积直径不同,其他条件相同的状态3、状态6与状态7的水滴撞击特性计算结果进行分析.当其他飞行条件与气象条件保持不变,水滴平均容积直径减小时,水滴本身具有的惯性减小.所以,随着水滴平均容积直径的减小,惯性力对于水滴运动轨迹的影响会减小,相应地水滴受到的黏性阻力与离心力对水滴运动轨迹的影响会增大.因此,得到如图 7所示,随着水滴平均容积直径的减小,旋转帽罩局部水滴收集系数明显减小,水滴撞击极限范围也明显减小的结果.
 |
| 图 7 状态3,6,7,y=0截面局部水滴收集系数分布曲线 Fig. 7 Case 3,6,7,cross-section y=0,local water droplets collect coefficient distribution curves |
由以上计算结果与分析可知,当来流速度较小,水滴平均容积直径较小时,水滴运动轨迹受黏性阻力与离心力的影响较大.
 |
| 图 8 状态8和9,y=0截面局部水滴收集系数分布曲线 Fig. 8 Case 8 and 9,cross-section y=0,local water droplets collect coefficient distribution curves |
为了使旋转帽罩转速对水滴撞击特性的影响更加明显,应该尽量减小惯性力对水滴运动轨迹的影响,加大黏性阻力与离心力对水滴运动轨迹的影响.因此采用来流速度较小,水滴平均容积直径较小的状态8、状态9进行比较.由图 8可以看出,随着旋转帽罩转速的增大,驻点处局部水滴收集系数略有减小,并且水滴撞击极限也 有所减小. 4 结 论
1) 欧拉法求解旋转帽罩的空气-水滴两相流流场,继而得到旋转帽罩的水滴撞击特性是可行的.
2) 在飞行条件与气象条件保持不变的情况下,随着旋转帽罩转速增加,水滴撞击极限稍有减小.
3) 在典型飞行条件与气象条件下,离心力对水滴运动轨迹的影响远小于惯性力对其的影响.旋转帽罩转速对水滴的运动轨迹影响较小,水滴撞击特性变化不大.
| [1] | 裘燮纲,韩凤华.飞机防冰系统[M].北京:航空专业教材编审组,1985:50-53Qiu Xiegang,Han Fenghua.Aircraft anti-icing system[M].Beijing:Compilation and Examination Group of Aero Specialized Teaching Materials,1985:50-53(in Chinese) |
| [2] | 周力行.湍流气粒两相流动和燃烧的理论与数值模拟[M].北京:科学出版社,1994:293-298Zhou Lixing.Theory and numerical modeling of turbulent gas-particle flows and combustion[M].Beijing:Science Press,1994:293-298(in Chinese) |
| [3] | 杨倩,常士楠,袁修干.水滴撞击特性的数值计算方法研究[J].航空学报,2002,23(2):173-176 Yang Qian,Chang Shinan,Yuan Xiugan.Study on numerical method for determining the droplet trajectories[J].Acta Aeronuatica et Astronautica Sinca,2002,23(2):173-176(in Chinese) |
| Cited By in Cnki (69) | |
| [4] | 常士楠,苏新明,邱义芬.三维机翼结冰模拟[J].航空学报,2011,32(2):212-222 Chang Shinan,Su Xinming,Qiu Yifen.Ice accretion simulation on three dimensional wings[J].Acta Aeronautica et Astronautica Sinca,2011,32(2):212-222(in Chinese) |
| Cited By in Cnki (9) | |
| [5] | 申晓斌,林贵平,杨胜华.三维发动机进气道水滴撞击特性分析[J].北京航空航天大学学报,2011,37(1):1-5 Shen Xiaobin,Lin Guiping,Yang Shenghua.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) |
| Cited By in Cnki (3) | |
| [6] | Bourgault Y,Boutanios Z,Habashi W G.Three-dimensional eulerian approach to droplet impingement simulation using FENSAP-ICE,part 1:model, algorithm, and validation[J].Journal of Aircraft,2000,37(1):95-103 |
| Click to display the text | |
| [7] | 赵秋月,董威,朱剑鋆.发动机旋转整流帽罩的水滴撞击特性分析[J].燃气涡轮试验与研究,2011,24(4):32-35 Zhao Qiuyue,Dong Wei,Zhu Jianjun.Droplets impinging characteristic analysis of the rotating fairing of aero-engine[J].Gas Turbine Experiment and Research,2011,24(4):32-35(in Chinese) |
| Cited By in Cnki (2) | |
| [8] | Das K,Hamed A,Basu D.Ice shape prediction for turbofan rotating blades[R].AIAA 2006-0209,2006 |
| [9] | 杨胜华,林贵平,申晓斌.三维复杂表面水滴撞击特性计算[J].航空动力学报,2010,25(2):284-290 Yang Shenghua,Lin Guiping,Shen Xiaobin.Water droplet impingement prediction for three-dimensional complex surfaces[J].Journal of Aerospace Power,2010,25(2):284-290(in Chinese) |
| Cited By in Cnki (10) | |
| [10] | Slater S A,Young J B.The calculation of inertial particle transport in dilute gas-particle flows[J].International Journal of Multiphase Flow,2001,27(1):61-87 |
| Click to display the text | |
| [11] | Wirogo S,Srirambhatla S.An Eulerian method to calculate the collection efficiency on two and three dimensional bodies[R].AIAA 2003-1073,2003 |
| [12] | Tong X,Luke E A.Eulerian simulation of icing collection efficiency using a singularity diffusion model[R].AIAA 2005-1246,2005 |
| [13] | Cao Y,Ma C,Zhang Q,et al.Numerical simulation of ice accretions on an aircraft wing[J].Aerospace Science and Technology,2012,23(1):296-304 |
| Click to display the text | |
| [14] | Linke D A.Systems of commercial turbofan engines[M].Berlin:Springer,2008:179-183 |
| [15] | 李西园.旋转部件结冰缩比试验相似性研究[D].北京:北京航空航天大学,2010 Li Xiyuan.Similarity of icing scaling test for aircraft rotating components[D].Beijing:Beijing University of Aeronautics and Astronautics,2010(in Chinese) |