基于欧拉法模拟旋转帽罩水滴撞击特性

吴孟龙 常士楠 冷梦尧 王 超

(北京航空航天大学 航空科学与工程学院，北京100191)

飞机在含有大量过冷水滴的云层中飞行时，其迎风面会发生结冰现象.发动机进口旋转帽罩处于飞机的迎风面，也会发生结冰现象.旋转帽罩结冰会降低发动机入口气流品质，使发动机性能降低.并且，旋转帽罩表面结冰脱落，被吸入发动机内部，会导致发动机损毁，造成重大飞行事故［1］.因此十分有必要对旋转帽罩进行结冰研究.飞机结冰数值模拟研究一般分为空气流场求解、水滴轨迹及撞击特性求解、结冰表面热平衡分析以及结冰冰形计算4个部分.本文主要关注于旋转帽罩水滴轨迹及撞击特性分析.

过冷水滴在空气流场中的运动属于气粒两相流动.按照对水滴颗粒的不同处理方式，主要有两种气粒两相流动研究方法［2］:

1)把气体当作连续介质，而将颗粒视为离散体系，即拉格朗日法;

2)把气体与颗粒都看成共同存在且相互渗透的连续介质，把颗粒作为拟流体，即欧拉法.

目前，国内外对于静止部件的水滴撞击特性的求解既有采用拉格朗日法［3］，也有采用欧拉法［4-6］.而对于旋转部件的水滴撞击特性，国内外研究较少，并且多数采用拉格朗日法［7-8］.

本文在欧拉法求解静止部件水滴撞击特性的基础上，提出一种旋转帽罩水滴撞击特性数值模拟方法，实现了欧拉法求解旋转帽罩表面的水滴撞击特性.通过对不同转速旋转帽罩水滴撞击特性进行数值模拟研究，发现旋转帽罩转速越大，其表面水滴撞击极限越小.但是，在典型飞行条件与气象条件下，离心力对水滴运动轨迹的影响远小于惯性力对其的影响，旋转帽罩转速对水滴撞击特性影响较小.

1 水滴相数学模型

基于欧拉法在静止部件水滴撞击特性的成功应用［9-12］，本文提出了一种旋转帽罩水滴撞击特性求解方法.本文借助于Fluent软件，采用S-A湍流模型实现旋转帽罩空气流场的求解.而将水滴视为拟流体，采用欧拉法，借助于Fluent UDS模块，实现旋转帽罩水滴流场的求解.

结冰条件下，由于水滴的体积分数一般在10-6量级，可认为空气和水滴之间是单向耦合的，即空气流场影响水滴的运动，而水滴的运动对空气流场没有影响［9，13］.这使得空气相和水滴相可以分离求解:先得到空气流场，再计算水滴流场.并且对水滴进行一些合理的假设:①水滴是球形的，无变形或破裂;②水滴之间无碰撞、聚合，撞击到壁面后无飞溅;③水滴和周围空气无质量和热量的传递;④只考虑作用在水滴上的空气阻力.

1.1 控制方程

基于以上假设，稳态水滴相控制方程为

式(1)右边项为耗散函数，主要作用为消除数值计算过程中产生的奇异值，使计算可以顺利进行并得到合理的数值，在1.3节中将会详细阐述.式中，α是水滴体积分数;ρd和ud分别为水滴的密度和速度;ua为空气速度;K为水滴惯性系数，其表达式为

其中μa为空气的动力黏性系数.基于假设①认为水滴为球形，因此可应用球形水滴的阻力系数函数f，其表达式为

其中，水滴的阻力系数CD为

相对雷诺数Red计算公式为

1.2 边界条件

旋转帽罩空气-水滴流场数值模拟中，空气相与水滴相的边界条件处理方式不同.其中，空气相入口采用速度入口边界条件，出口采用出流边界条件，壁面采用动壁面中的旋转、无滑移壁面边界条件.

水滴相的入口采用速度入口边界条件，速度大小与空气速度保持一致;出口则采用出流边界条件.水滴相的壁面边界条件需要进行特殊处理.这是由于水滴撞到壁面后，会积聚到壁面上，从而从流场中离开.因此，在水滴撞击区域，如图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

1.3 耗散函数

由于水滴相控制方程属于双曲型方程［10］，数值迭代过程会出现局部数值奇异，因此常采用一些耗散函数，将局部奇异值耗散，使计算能够得到收敛.本文采用的耗散函数［9］为

其中，αp为当前网格的水滴体积分数;αNi为周围网格的水滴体积分数;a为经验常数，a取值过大，会使得计算结果误差增大，a取值过小，则会使计算结果发散.本文根据文献［9］中a的取值，并通过与Fensap进行算例计算结果比较，取a=0.001.本文通过编写UDF函数，改变材料项中的耗散值，从而将耗散函数添加到水滴相控制方程.

1.4 局部水滴收集系数计算

局部水滴收集系数β为当地表面微元实际水滴收集率与最大可能水滴收集率的比值，其计算公式为

式中，n为壁面单元的单位法向量;V∞为自由来流的水滴速度;α为壁面单元的水滴体积分数;α0为入口处单元的水滴体积分数.

2 计算模型与计算条件

2.1 网格划分

由于椭圆型帽罩较易结冰［14］，因此本文采用椭圆型旋转帽罩作为计算数模，其中长轴半径为0.09 m，短轴半径为0.06 m.在网格划分过程中，对旋转帽罩尾部进行适当的延伸.在旋转帽罩前缘、尾部以及周围划分10倍远场网格.所划网格如图2所示.

图2 旋转帽罩计算域网格Fig.2 Computational mesh of the spinner

2.2 计算条件

本文关注重点在于旋转速度对于旋转帽罩水滴撞击特性的影响，因此，取不同旋转速度的飞行及气象条件作为计算条件，分析研究旋转速度对旋转帽罩水滴撞击特性的影响.参照典型民用航空发动机旋转帽罩转速，选取旋转帽罩转速与其他飞行条件如表1所示.

表1 飞行及气象条件Table 1 Flight and weather conditions

3 计算结果与分析

利用上述计算方法对表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，速度78 m/s，转速5000 r/min水滴轨迹图Fig.5 Case 3，velocity=78 m/s，rotate speed=5000 r/min，the scheme of droplet trajectory

3.2 来流速度对水滴撞击特性影响

选取来流速度不同，其他条件相同的状态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.3 水滴平均容积直径对水滴撞击特性影响

选取水滴平均容积直径不同，其他条件相同的状态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、状态9进行比较.由图8可以看出，随着旋转帽罩转速的增大，驻点处局部水滴收集系数略有减小，并且水滴撞击极限也有所减小.

图8 状态8和9，y=0截面局部水滴收集系数分布曲线Fig.8 Case 8 and 9，cross-section y=0，local water droplets collect coefficient distribution curves

4 结论

1)欧拉法求解旋转帽罩的空气-水滴两相流流场，继而得到旋转帽罩的水滴撞击特性是可行的.

2)在飞行条件与气象条件保持不变的情况下，随着旋转帽罩转速增加，水滴撞击极限稍有减小.

3)在典型飞行条件与气象条件下，离心力对水滴运动轨迹的影响远小于惯性力对其的影响.旋转帽罩转速对水滴的运动轨迹影响较小，水滴撞击特性变化不大.

References)

［1］裘燮纲，韩凤华.飞机防冰系统［M］.北京:航空专业教材编审组，1985:50-53 Qiu 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-298 Zhou Lixing.Theory and numerical modeling of turbulent gasparticle 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)

［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)

［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)

［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

［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)

［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)

［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

［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

［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)