旋转部件复杂表面水滴撞击计算

申晓斌张志强林贵平穆作栋卜雪琴

(1.北京航空航天大学,北京 100191;2.人机工效与环境控制重点学科实验室,北京 100191)

旋转部件复杂表面水滴撞击计算

申晓斌1,2,*,张志强1,2,林贵平1,2,穆作栋1,2,卜雪琴1,2

(1.北京航空航天大学,北京 100191;2.人机工效与环境控制重点学科实验室,北京 100191)

为模拟旋转体结冰问题的过冷水滴运动及撞击过程,基于欧拉方法及单旋转坐标系模型,建立了三维旋转水滴运动模型,并提出了相应的数值求解方法。采用单旋转坐标系对空气及水滴两相流动过程进行处理,通过引入惯性力,将惯性系下的周期转动边界转换为定常流动边界;利用欧拉方法,使用单向耦合形式描述空气—水滴流场;在单旋转坐标系下,向控制方程内引入科里奥利加速度及牵连加速度,进行非惯性系下欧拉方程的修正,从而描述水滴运动过程;采用有限容积求解器对空气及水滴运动的控制方程组进行求解,通过引入源项定义单旋转坐标系下的惯性力,得到空气流场及水滴场的速度、体积分数分布,进而得到表面水滴撞击特性。采用上述方法对旋转帽罩与叶片模型进行算例分析,结果表明所建立的旋转水滴计算方法有效,对比静止状态表面的结果,旋转对帽罩的水滴撞击特性影响甚小,而对桨叶存在显著影响;由于帽罩具有中心对称的特性,因而旋转带来的切向速度变化对其水滴撞击特性影响不明显;桨叶表面水滴收集系数随旋转角速度增大而增大,同时收集系数在表面的分布会向迎风方向偏移,较大的角速度对应了更为显著的收集系数增幅与偏移现象。

飞机结冰;数值模拟;欧拉法;旋转部件;水滴撞击特性;单旋转坐标系模型

Keywords:aircraft icing;numerical simulation;Eularian method;rotating components; water collection efficiency;rotating reference frame

0 引 言

随着航空科学与技术的发展,以飞机结冰气象条件为典型代表的气象因素成为威胁飞行安全的重要因素。过冷水滴撞击导致的飞机结冰现象可能破坏升力部件气动外形,影响发动机进气效率,甚至引起飞行事故。飞机结冰过程及其防护方法研究涉及复杂的传热传质耦合问题,近年来得到了极大的重视,国家重点基础研究发展计划、国家安全重大基础研究计划、国家自然基金等都将其纳入了研究计划。本文针对飞机结冰的重要环节,即水滴撞击特性分析,开展计算方法研究。

水滴撞击特性决定了结冰的范围与质量,其计算是飞机结冰与防除冰仿真分析的基础。国内外开展了大量的理论分析与计算方法研究,发展了拉格朗日法与欧拉法两种主要的水滴撞击计算方法,开发了一系列的计算分析软件[1-4],并将其应用于各种复杂的飞机部件表面,取得了较好的效果。目前,固定部件表面的水滴撞击计算方法已经基本成熟,其准确性得到了学者们的共识。

对于旋转部件,特别是具有非中心对称特点的旋转表面(如螺旋桨、发动机旋转叶片等)而言,由于部件旋转引入的离心力等惯性力影响,其水滴撞击过程与固定表面存在显著差异,针对旋转部件的水滴撞击研究分析较少。加拿大的FENSAP-ICE软件[5-6]对旋转部件结冰问题的关注较早,分析了中心旋转体及复杂旋转表面的水滴撞击特性。国内,赵秋月、董威等[7]用拉格朗日法对发动机旋转整流帽罩的水滴撞击特性进行了分析。王治国等[8]应用CFX软件及其粒子输运模型计算了发动机旋转表面的水滴撞击特性。易贤等[9]采用多参考坐标系方法下的空气流场结合欧拉法下的水滴计算,得到了风力机三维旋转表面水滴收集率。吴孟龙、常士楠等[10]基于欧拉法模拟计算了旋转帽罩的水滴撞击特性。王健、胡娅萍等[11]对旋转整流罩的结冰生长开展了试验研究。

本文采用单旋转坐标系模型来研究复杂表面稳定旋转时的空气流动与水滴运动,将其抽象为旋转坐标系下的定常流动,用欧拉法来计算水滴撞击特性,并分析水滴撞击计算方法的有效性及旋转速度对计算结果的影响规律。

1 旋转流数学模型

固定部件绕流流动模型通常是建立在惯性参考坐标系中,将其扩展到旋转体流动通常有四种处理方法:单旋转坐标系模型、多参考系坐标模型、混合平面模型及滑移网格模型。其中单旋转坐标系模型将坐标固定在旋转部件上,坐标随部件以相同的速度转动,因此旋转边界相对于坐标系是静止的,非定常的旋转流动转化为定常流动,从而简化计算过程。本文采用单旋转坐标系模型,建立旋转坐标系下的空气流动与水滴运动模型,发展旋转部件复杂表面水滴撞击计算方法。

1.1 空气流动模型

考虑部件匀角速度定轴旋转,建立单旋转参考坐标系。当牵连运动为定轴转动时,动点的绝对加速度等于它的牵连加速度、相对加速度和科里奥利加速度的矢量和。因此,单旋转坐标系下定常空气流动的控制方法变为[12]:

式中,ρa为空气密度;p为压力;τr为基于相对速度的粘性应力;F为外力矢量;ω为旋转速度;r为控制体矢径;var为空气相对速度,由速度合成定理可得:

其中,va为空气绝对速度;ur为牵连速度,由下式计算:

此外,式(2)中的2ω×var为科里奥利加速度;ω×(ω ×r)为牵连加速度,即向心加速度。

1.2 水滴运动模型

在结冰气象条件下,过冷水滴随气流一起运动,迎风部件表面附近的水滴可能撞击到飞机壁面。由于水滴尺寸较小,空气中的液态水含量少,只考虑空气流动对水滴运动的影响,忽略水滴对气流的影响(单向耦合)。对于复杂表面的水滴撞击计算,欧拉法有显著的优势。欧拉法将水滴视为连续相,引入水滴容积分数来表示控制体内的水滴体积占总体积的大小,建立控制方程。用与空气流场同样的处理方法将水滴欧拉方程扩展到单旋转坐标系下,得:

式中,α是水滴容积分数;ρ为水滴的密度;vr为水滴相对速度,由速度合成定理可知:

其中,v为水滴绝对速度。

式(6)中的K为空气-水滴交换系数,由下式定义[13]:

其中,μa表示空气的动力粘度;dp为水滴直径;f是阻力函数,可由下式计算:

其中,阻力系数CD由下式计算[13]:

Re为相对雷诺数,由下式给出:

1.3 水滴收集系数计算模型

用欧拉法求解水滴控制方程后,可直接得到旋转部件表面的水滴撞击区域和各处的水滴撞击量,用于飞机结冰及防除冰计算分析。为便于分析比较,计算撞击表面水滴收集系数,其定义为表面上的实际水收集量与该表面上最大可能的水收集量之比。根据定义可推出旋转坐标系下的局部水滴收集系数β的表达式为[9]:

式中,n为局部壁面的单位法向矢量;V∞为远场的水滴速度;α∞为远场的水滴容积分数。

由式(12)与式(7)分析可知:旋转运动时的表面水滴收集系数不仅与远场空气速度有关系,还受旋转速度的影响;当旋转速度较大时,β值可能大于1[9],即旋转表面的水收集量比表面固定时最大可能的水收集量还要大。

2 计算条件与模型实现方法

2.1 几何模型与网格

以发动机进气道内部的旋转帽罩与叶片为研究对象,建立简化的几何模型,如图1所示。图中的外部圆柱轮廓表示发动机腔的内壁面,中间位置用60°的三角圆锥来表示发动机整流帽罩,用三个互为120°的椭圆柱来表示旋转叶片。为保证计算的收敛性,将发动机腔壁向前后分别进行拉长处理,前后分别为速度进口与压力出口。

将几何模型导入ICEM软件中进行网格划分,选用非结构方法,并在内外壁面进行边界层网格的设置。由于水滴撞击的区域一般在迎风表面,在整流帽罩与椭圆柱表面进行网格加密处理,最终的网格结果如图2所示。

图1 旋转部件几何模型Fig.1 Rotation part geometry model

图2 旋转部件表面网格Fig.2 Surface grid of rotation part

2.2 旋转模型实现方法

以通用流体计算软件FLUENT为平台,通过二次开发来实现旋转部件复杂表面的水滴撞击特性计算。由于空气流动与水滴运动单向耦合,将空气与水滴进行分步计算。

发动机整流帽罩与叶片绕中心轴以恒定速度转动,以转轴为基础建立单旋转坐标系。采用三维粘性不可压流求解器,选用k-ε湍流模型来迭代求解旋转坐标系下的守恒方程式(1)、式(2),得到旋转部件周围的空气绕流流场。需要注意的是,帽罩与叶片随旋转坐标系以相同的速度转动,其相对速度为0;而原本固定的发动机腔壁,在旋转坐标系下则是以相同的转速反向转动。

在空气流场的基础上,用欧拉法进行单旋转坐标系下的水滴流动。参考惯性系下的水滴求解方法,应用FLUENT软件提供的UDS(User Defined Scalar,用户自定义标量)输运方程来进行旋转水滴计算[13-14]。分析旋转下的欧拉方程,将水滴容积分数、水滴相对速度沿坐标轴的三个分量分别设为自定义标量。同时将旋转坐标下产生的牵连项与科里奥利项放入水滴动量方程的源项中,建立UDS框架下的水滴输运方程。通过有限容积法对输运方程求解,便可以得到旋转作用下的水滴场分布。需要注意的是,进行水滴撞击壁面后的吸水处理时,考虑旋转作用的影响,用水滴的相对速度来进行。

2.3 计算状态

首先用基本成熟的固定部件水滴撞击计算方法对几何模型进行静止时的计算。之后设一个非常小的旋转速度,应用建立的单旋转坐标系模型进行求解,与传统的固定坐标方法进行比较,验证有效性。然后,取其他两个不同旋转速度进行计算,分析旋转速度对水滴撞击特性的影响。

模型计算时,速度进口条件取30 m/s,出口压力取标准大气压,水滴直径取20μm,旋转计算时转速分别取0.01 rad/s、30 rad/s、60 rad/s。

3 计算结果与分析

3.1 旋转模型验证

图3所示为当叶片静止时,用传统固定坐标系模型[13,15]计算得到的局部水滴收集系数分布结果。从图3中可以看出,水滴主要撞击区域为帽罩的前缘点以及叶片的前缘驻点附近。

图4所示为旋转速度为0.01 rad/s时,用单旋转坐标系模型计算得到的帽罩与叶片表面速度分布结果。从图4中可以看到,旋转使壁面产生运动,离旋转轴越远速度越大。

图3 表面水滴收集系数分布云图(固定)Fig.3 Contours of collection efficiency(static)

图4 表面速度分布云图(ω=0.01 rad/s)Fig.4 Contours of velocity(ω=0.01 rad/s)

图5为通过旋转坐标下水滴撞击模型计算得到的转速为0.01 rad/s时的水滴收集系数分布结果。从图5中可以看到分布结果与静止时传统方法的结果类似。为更好地与传统固定参考系模型进行比较,取叶片的一个截面来进行显示,截面的位置为图4中上边叶片较亮线处。

图6为固定参考系模型计算得到的截面结果与旋转模型计算转速为0.01 rad/s时水滴收集系数结果的比较。从图6中可以看到水滴撞击范围与收集系数分布都非常一致,在一定程度上说明了所建立的单旋转坐标水滴撞击模型是有效的。

图5 表面水滴收集系数分布云图(ω=0.01 rad/s)Fig.5 Contours of collection efficiency(ω=0.01 rad/s)

图6 旋转计算与静止计算收集系数比较Fig.6 Comparison between collection efficient results

3.2 旋转速度的影响

图7 表面水滴收集系数分布云图(ω=30 rad/s)Fig.7 Contours of collection efficiency(ω=30 rad/s)

图8 表面水滴收集系数分布云图(ω=60 rad/s)Fig.8 Contours of collection efficiency(ω=60 rad/s)

图7与图8分别为旋转速度为30 rad/s与60 rad/s时(旋转速度方向为图中+x方向),表面水滴收集系数分布云图。从图中可以看到旋转作用下的最大局部水滴收集系数比静止时的结果大,且当转速为60 rad/s,水收集系数存在大于1的区域。图中叶片表面水滴撞击区域向右侧表面偏移,体现了旋转速度对水滴运动的影响,而帽罩由于中心对称,水滴收集系数变化较小。

着重分析旋转速度对叶片的影响,取上文所述截面进行比较显示,结果如图9所示。从图9中可以看到,截面上的局部水收集系数向-y轴方向移动,且旋转速度越大,区域移动越大。这是由于旋转作用使叶片表面的相对速度发生变化,从而改变驻点位置,且转速越快,合成速度产生的驻点变化越大。

图9 不同旋转速度水收集系数比较Fig.9 Comparison between collection efficient results at different rotational speeds

4 结 论

采用单旋转坐标系建立了旋转表面水滴撞击模型,并对发动机旋转帽罩与叶片进行了计算分析。结果表明:所建立的方法能有效地计算旋转部件复杂表面的水滴撞击特性;旋转作用使水滴能撞击区域向相对迎风的叶片侧表面移动,且旋转速度越大,水滴撞击区域移动越大。所建立模型的计算准确性还需要经过试验来验证。

[1]Gent R W.TRAJICE2-a combined water droplet trajectory and ice accretion prediction program for aerofoils[R].RAE TR 90054,1990.

[2]Hedde T,Guffond D.ONERA three-dimensional icing model [J].AIAA Journal,1995,33(6):1038-1045.

[3]Morency F,Beaugendre H,Habashi W G.FENSAP-ICE:a study of the effect of ice shapes on droplet impingement[R]. AIAA 2003-1223,2003.

[4]Bidwell C.Ice particle transport analysis with phase change for the E3 turbofan engine using LEWICE3D version 3.2[R]. AIAA 2012-3037,2012.

[5]Thomas R,Habashi W G.FENSAP-ICE simulation of icing on wind turbine blades,part 1:performance degradation[R]. AIAA 2013-0750,2013.

[6]David S,Habashi W G.FENSAP-ICE simulation of complex wind turbine icing events,and comparison to observed performance data[R].AIAA 2014-1399,2014.

[7]Zhao Q Y,Dong W,Zhu J J.Droplets impinging characteristic analysis of the rotating fairing of aero-engine[J].Gas Turbine Experiment and Research,2011,24(4):32-35.(in Chinese)赵秋月,董威,朱剑鋆.发动机旋转整流帽罩的水滴撞击特性分析[J].燃气涡轮试验与研究,2011,24(4):32-35.

[8]Wang Z G,Lou D C,Guo W.Numerical simulation of water droplet impingement characteristic on aero-engine rotating machinery[J].Gas Turbine Experiment and Research,2013, 26(01):35-39.(in Chinese)王治国,娄德仓,郭文.发动机旋转表面水滴撞击特性数值研究[J].燃气涡轮试验与研究,2013,26(01):35-39.

[9]Yi X,Wang K C,Ma H L,et al.3-D ementional simulation of droplet collection efficiency in large-scale wind turbine icing[J]. Acta Aerodynamica Sinica,2013,31(6):745-751.(in Chinese)易贤,王开春,马洪林,等.大型风力机结冰过程水滴收集率三维计算[J].空气动力学学报,2013,31(6):745-751.

[10]Wu M L,Chang S N,Leng M Y,et al.Simulation of droplet impingement characteristics of spinner based on Eulerian method[J].Journal of Beijing University of Aeronautics and Astronautics,2014,40(9):1263-1267.(in Chinese)吴孟龙,常士楠,冷梦尧,等.基于欧拉法模拟旋转帽罩水滴撞击特性[J].北京航空航天大学学报,2014,40(9):1263-1267.

[11]Wang J,Hu Y P,Ji H H,et al.Experiment of ice accretion and shedding on rotating spinner[J].Journal of Aerospace Power,2014,29(06):1352-1357.(in Chinese)王健,胡娅萍,吉洪湖,等.旋转整流罩积冰生长与脱落过程的实验[J].航空动力学报,2014,29(06):1352-1357.

[12]ANSYS Inc.ANSYS Fluent 13.0 UDF manual[M].New Hapmpshire:ANSYS Inc,2012.

[13]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,31(1):1-5.(in Chinese)申晓斌,林贵平,杨胜华.三维发动机进气道水滴撞击特性分析[J].北京航空航天大学学报,2011,31(1):1-5.

[14]ANSYS F.Fluent user’s manual[J].Software Release,2006, (6):449-456.

[15]Yang S H,Lin G P,Shen X B.An Eulerian method for water droplet impingement prediction and its implementations[C]// Proceeding of the 1st International Symposium on Aircraft Airworthiness.ISAA 2009 2-045,2009.

Droplet impingement calculation on complex suface of rotating part

Shen Xiaobin1,2,*,Zhang Zhiqiang1,2,Lin Guiping1,2,Mu Zuodong1,2,Bu Xueqin1,2
(1.Beihang University,Beijing100191,China; 2.Fundamental Science on Ergonomics and Environment Control Laboratory,Beijing100191,China)

A three-dimensional water impingement model for rotating surface has been established and solved to simulate the process of supercooled water droplet impinging on ice protection surfaces,in which an Eularian method and a single rotating coordinate system were applied.With the single rotating coordinate system,the unsteady cyclically rotating boundary condition was simplified to a steady boundary by introducing the inertial forces;with the Eularian method,the air flow and water droplets flow field were depicted in a one-way coupled form;the Eularian control equations were adjusted in the non-inertial reference system by introducing the Coriolis acceleration and convected acceleration;then a finite volume solver was applied to solve the governing equations of air flow and droplet flow,in which the inertial forces were defined as the source terms of the momentum equation,and the results of velocity and volumetric fraction of air flow and droplet flow were obtained as well as the droplet collection efficiency.With the mentioned method,a rotating model of cowling and blades was simulated for a stationary case and at different rotational speeds.The results show that the rotation has significant influences on the water collection efficiency of the blades.The value of collection efficiency increases with increasing rotational velocity since droplet velocity normal to the surface increases due to the rotation,and the distribution curve shifts to the windward direction;the effect on the cowling is slight,and a distinguishable change of the impingement characteristic has not been caused by the change of the tangential velocity due to the axis-symmetric shape of the cowling.

V211.3

A

10.7638/kqdlxxb-2015.0213

0258-1825(2016)06-0709-05

2015-12-21;

2015-12-25

中央高校基本科研业务费专项资金(YWF-14-HKXY-008)

申晓斌*(1985-),湖南怀化人,男,讲师,研究方向:飞机结冰与防除冰技术.E-mail:sxb762@163.com

申晓斌,张志强,林贵平,等.旋转部件复杂表面水滴撞击计算[J].空气动力学学报,2016,34(6):709-713.

10.7638/kqdlxxb-2015.0213 Shen X B,Zhang Z Q,Lin G P,et al.Droplet impingement calculation on complex suface of rotating part[J].Acta Aerodynamica Sinica,2016,34(6):709-713.