航空发动机旋转整流帽罩水滴撞击特性计算

Calculation of water droplet impingement characteristics on rotating fairing of aero-engine

摘要

航空发动机整流帽罩作为飞机发动机结冰的关键部件之一，其表面的水滴撞击特性直接影响后续的结冰状态。为了得出整流帽罩在不同工况下的水滴撞击特性，采用欧拉法和单旋转坐标系，建立了航空发动机旋转整流帽罩的三维水滴撞击模型。分别模拟了整流帽罩在静止及旋转工况下的水滴撞击特性。结果表明，在静止工况下，随着来流马赫数的增加，水滴收集系数和撞击区均显著增大。相比单一直径，当考虑水滴的Langmuir-D分布时，驻点处水滴收集系数减少，撞击区域反而增加，下游区域的收集系数也会略有增大。在旋转工况下，由于整流帽罩的流线型外形特征，转速对水滴撞击特性的影响可以忽略不计。

Abstract

As one of the critical components of aircraft engine icing，the water droplet impingement characteristics on the rotating fairing surface directly influences the subsequent icing state.To investigate the water droplet impingement behavior under different operating conditions，a three-dimensional water droplet impingement model for the rotating fairing of aero-engine was established by using the Euler method and single rotating coordinate system.The water droplet impingement characteristics of the rotating fairing under stationary and rotating conditions were simulated respectively.The results show that under stationary conditions，as the freestream Mach number increases，both the water droplet collection coefficient and impingement area​ exhibit significant growth； compared to single-diameter droplets，when considering the Langmuir-D distribution​​of water droplets，the collection coefficient at the stagnation point decreases，while the impingement zone increases and the downstream collection coefficient also slightly increases.Under rotating conditions，the​​rotational speed​​ has negligible effects on droplet impingement characteristics due to the streamlined aerodynamic profile of the rotating fairing.

Graphical abstract

关键词

旋转整流帽罩 / 水滴撞击 / 欧拉法 / 单旋转坐标系 / 水滴收集系数

Key words

rotating fairing / water droplet impingement / Euler method / single rotating coordinate system / water droplet collection coefficient

随着航空科技的发展，以飞机结冰为典型代表的相关问题，对飞行安全构成了严重威胁。因受到过冷水滴撞击，发动机的气动外形可能会出现结冰现象，使发动机进气效率受到影响，甚至可能造成飞行事故。涉及传热传质耦合复杂问题的飞机结冰过程及其防护措施，近年来备受关注。

国内外学者进行了大量的理论分析及方法研究，发展了基于拉格朗日方程和欧拉方程两种主要数值计算方法，并将其应用于相关计算分析软件系列程序中，在分析飞机机翼表面^［1-3］及发动机表面^［4］的水滴撞击特性方面取得了良好的进展。目前，用于计算固定表面水滴撞击特性的方法已经基本完善和成熟，其精度已经得到了许多实验的验证。

对于旋转部件，特别是轴对称旋转表面，例如旋转整流帽罩和航空发动机叶片，由于受切向速度引入的离心力影响，水滴撞击与固定表面存在明显差异。关于旋转部件水滴撞击的研究相对较少。Reid等^［5］、Switchenko等^［6］较早地关注了旋转部件的结冰问题，利用FENSAP-ICE软件分析了轴对称旋转表面和其他复杂旋转表面的水滴撞击特征。赵秋月等^［7］采用拉格朗日法对航空发动机旋转整流帽罩的水滴撞击特性进行了分析。王治国等^［8］利用CFX模型和粒子输运模型对航空发动机旋转面水滴撞击进行了计算。易贤等^［9］采用多参考坐标系计算了空气流场，并利用欧拉法对风力机三维旋转叶片的集水效率进行了计算。吴孟龙等^［10］基于欧拉法模拟了水滴在旋转整流帽罩上的撞击。王健等^［11］对旋转整流罩的积冰现象进行了实验研究。

本文针对飞机结冰问题的部分重要环节，即整流帽罩的水滴撞击特性，开展了模拟仿真计算。通过引入惯性力将非稳态周期性旋转边界条件简化为稳态边界条件。用欧拉法将空气流场和水滴流场描述为单向耦合形式。在非惯性参考系中引入科氏加速度以修正欧拉控制方程。使用有限体积求解器求解空气、水滴流动的控制方程，将惯性力定义为动量方程的源项，获得空气、水滴流动的速度、体积分数及水滴收集系数。采用上述方法，模拟计算并分析了航空发动机整流帽罩在静止和不同转速情况下的水滴撞击特性。

1 旋转流场的数学模型

惯性坐标系下通常建立固定部件的绕流流动模型。将非旋转条件扩展为旋转条件的方法有4种，分别是单旋转坐标模型、多参考坐标系模型、混合平面模型及滑动网格模型。其中，单旋转坐标模型是假定坐标固定于旋转元件上，坐标旋转时的切向速度与旋转元件的切向速度一致。因此，本文将旋转边界转换为静止参考坐标，将非定常流动转换为定常流动，从而简化模型的处理和计算过程。

1.1 空气流场模型

首先，考虑部件匀角速度定轴旋转，建立单旋转参考坐标系。当牵连运动为定轴运动时，绝对加速度等于牵连加速度、相对加速度和科氏加速度的总和。因此，在单旋转坐标系下的空气流场控制方程为^［11］

\nabla {\rho }_{\mathrm{a}}{\mathit{v}}_{\mathrm{a}\mathrm{r}}=\mathrm{0}（1）

\begin{array}{l}\nabla \left({\rho }_{\mathrm{a}}{\mathit{v}}_{\mathrm{a}\mathrm{r}}{\mathit{v}}_{\mathrm{a}\mathrm{r}}\right)+{\rho }_{\mathrm{a}}\left[\mathrm{2}\omega {\mathit{v}}_{\mathrm{a}\mathrm{r}}+\omega \left(\omega \times \mathit{r}\right)\right]\\ =-\nabla p+\nabla {\mathit{\tau }}_{\mathrm{r}}+\mathit{F}\end{array}（2）

式中：{\rho }_{\mathrm{a}}为空气密度；p为压力；{\mathit{\tau }}_{\mathrm{r}}为以相对速度为基础的黏性应力；\mathit{F}为外力矢量；\omega为转速；\mathit{r}为控制体积的失径；2\omegav_ar为科氏加速度；\omega（\omega \times \mathit{r}）为牵连加速度。{\mathit{v}}_{\mathrm{a}\mathrm{r}}可由式（3）计算得到。

{\mathit{v}}_{\mathrm{a}\mathrm{r}}={\mathit{v}}_{\mathrm{a}}-{\mathit{u}}_{\mathrm{r}}（3）

式中：{\mathit{\nu }}_{\mathrm{a}}为空气的绝对速度；{\mathit{u}}_{\mathrm{r}}为牵连速度，可由式（4）计算得到。

{\mathit{u}}_{\mathrm{r}}=\omega \mathit{r}（4）

1.2 水滴撞击模型

在结冰的天气条件下，过冷水滴随着气流一起运动，可能会撞击到飞机的迎风表面。由于水滴直径小、空气中水滴的含量少，因此忽略水滴对气流的影响，仅考虑空气流动对水滴运动的影响（单向耦合）。对于复杂表面，欧拉方法具有明显的优势，它将水滴视为连续相，通过引入水滴体积分数来表示控制体内的水滴体积占总体积的大小，从而建立控制方程。水滴的欧拉方程在单旋转坐标系中表示为^［12］

▽ραv_r=0（5）

▽（ραv_rv_r）+ρ［2\omegav_r+\omega (\omega \times \mathit{r})］

=ραK（ v_ar- v_r）（6）

式中：\alpha为水滴体积分数；ρ为水的密度； v_r为水滴的相对速度，通过式（7）计算。

v_r=v-u_r（7）

式中：\mathit{v}为水滴的绝对速度。式（6）中的K为空气-水交换系数，定义为^［13］

K=\frac{\mathrm{18}{\mu }_{\mathrm{a}}f}{\rho d_{\mathrm{p}}^{\mathrm{2}}}（8）

式中：{\mu }_{\mathrm{a}}为空气的动力黏度；d_{\mathrm{p}}^{}为水滴直径；f为阻力函数，通过式（9）计算。

f = \frac{C_{\mathrm{D}}Re}{\mathrm{24}}（9）

式中：C_{\mathrm{D}}为球体的阻力系数，可由式（10）计算得到。

C_D=\left\{\begin{array}{l}\mathrm{24}(\mathrm{1}+\mathrm{0.15}R{e}^{\mathrm{0.687}})/Re,Re\le \mathrm{1}\mathrm{000}\\ \mathrm{0.44},Re>\mathrm{1}\mathrm{000}\end{array}\right.（10）

Re为相对雷诺数，定义为

Re=\frac{\rho a|{\mathit{v}}_{ar}-{\mathit{v}}_r|d_p}{{\mu }_a}（11）

1.3 水滴收集系数模型

求解了水滴的欧拉方程之后，可直接得到旋转表面的水滴撞击区域及水滴撞击量，进一步计算出积冰和除冰热。将实际水滴收集值与局部表面可能的最大值之比定义为水滴的收集系数，使计算更为简单化。根据定义，旋转坐标下的水滴收集系数\beta可以描述为

\beta =\frac{\alpha }{{\alpha }_{\mathrm{\infty }}}·\frac{{\mathit{v}}_{\mathrm{r}}\mathit{n}}{V\mathrm{\infty }}（12）

式中：\mathit{n}为局部法向量；V\mathrm{\infty }为远场水滴速度；{\alpha }_{\mathrm{\infty }}为水滴在远场中的体积分数。

由式（7）和式（12）可知，旋转表面的水滴收集系数既与远场空气流速有关，也受旋转速度的影响。当转速足够大时，\beta可能大于1，说明旋转表面的收集系数可能大于固定表面的最大收集系数。

2 计算模型

2.1 几何和网格

建立航空发动机旋转整流帽罩的简化三维几何模型，如图1所示。图1中，半球体表面为整流帽罩。为了避免尾流对帽罩表面迎风部分的影响并增强计算的收敛性，帽罩采用如图1所示的圆锥体进行延伸。外表面描述远场边界条件，定义飞行、天气状况的相关数值。

本文采用ICEM软件对三维模型进行网格划分。由于六面体结构化网格在数值模拟中具有较强的优势，因此本文采用六面体结构化网格对流场进行拓扑划分。为了更准确地计算表面水滴撞击特性，对整流帽罩的表面网格进行局部加密。在下游和远场位置，适当增加网格尺寸以加快计算速度，网格细节如图2所示。

2.2 仿真模型计算方法

由于空气运动和水滴流动是单向耦合的，因此分两个单独步骤并利用Fluent软件完成计算。

发动机整流帽罩沿中心轴线以恒定速度旋转，建立以中心轴为基准的单旋转坐标系。采用三维可压缩流体求解器和k-\epsilon湍流模型来迭代求解方程（1）和方程（2），得到旋转部件周围的空气流场，由于整流帽罩与坐标系同步旋转，整流帽罩相对速度为零。实际的固定远场则表现为相反的转速。

在空气流场结果的基础上，采用更新的欧拉模型对水滴流动进行求解。采用Fluent的UDS（用户定义标量）输运方程进行了旋转条件下的水滴计算^［14-15］。水滴体积分数和相对速度矢量的3个投影标量被定义为4个UDSs。同时，将旋转引入的对流加速度和科氏加速度定义为动量方程的源项。通过求解4个UDSs输运方程，可以得到水滴的体积分数和速度，通过式（12）计算水滴收集系数。

2.3 计算条件

首先，计算在非旋转条件下不同来流速度的影响；其次，引入旋转边界条件进行对比。计算条件如表1所示。

所采用的欧拉模型只能求解单直径的水滴。但在真实的大气中，水滴的直径通常具有某种分布类型。此外，单个水滴直径不能反映由旋转引入的惯性影响。因此，本文采用L-D分布，其平均容积直径（mean volume cliameter，MVD）为11.5 \mathrm{\mu }\mathrm{m}，由7种不同直径的水滴组成。表2给出了直径分布和相应的体积分数。计算完7个直径的水滴收集系数后，最终的总收集系数可由式（13）计算。

\beta ={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{7}}}{p}_i{\beta }_i（13）

式中：p_i和{\beta }_i分别为各直径的体积分数和收集系数。

3 仿真结果与分析

3.1 非旋转情况下来流速度的影响

首先比较了不同马赫数下非旋转情况下的水滴体积分数和收集系数，如图3所示。由图3可以看出，随着马赫数的增加，下游的零体积分数区域变小，迎风面水滴撞击区域明显增大。这是因为随着马赫数的增加，尽管空气只能在整流帽罩表面绕流，但受惯性力的作用，水滴的运动将变得更加难以改变，使得水滴更有可能撞击到迎风面。利用式（12）计算出整流帽罩表面上的水滴收集系数，其云图如图3所示。可以看出，随着流速的增大，收集系数值和撞击区也随之增大。

为了进一步分析水滴收集系数，以对称面和整流帽罩表面的交线作为横坐标，不同马赫数下L-D分布和单直径水滴的收集系数变化曲线如图4所示。图4中x轴为沿曲线的长度，零值位于驻点，在零攻角的情况下，驻点即整流帽罩表面的中心点。图4中3条实心曲线表示由式（13）计算出的不同马赫数下的总收集系数值。为了比较满足L-D分布的液滴直径与具有相同MVD值的单个液滴直径之间的差异，图4中以离散点表示出了利用单一液滴直径计算获得的水滴收集系数。

由图4可以得出，随着马赫数的增加，相同位置处的水滴收集系数更大，驻点处的水滴收集系数最高。当马赫数为0.15时，驻点处水滴收集系数为0.3左右。当马赫数增加到0.6时，最大水滴收集系数几乎翻倍，约为0.55。水滴撞击量越大，融冰量越大，融冰系统热负荷越大。从图4还可以看出，由于水滴的惯性作用，撞击范围随速度的增大而变大。尽管差异不大，而且在下游位置的水滴收集系数值很小，但它仍然引入了更大的除冰区域。

水滴直径的L-D分布也与单一直径的计算结果存在很大的差异。从图4可以看出，L-D分布会导致高撞击区域特别是驻点处的收集系数值变小，但会扩大撞击区域并增加下游位置的收集系数值。因此，由直径分布所导致的水滴撞击特性的改点是显著的，在分析结冰或设计除冰系统时必须计算这些特性。

3.2 旋转的影响

本文计算了不同转速下空气流速大小和水滴体积分数，如图5所示。

旋转通过壁面剪切应力对空气流场的边界层产生影响，因此，不同的转速有着不同的壁面绝对速度。由于整流帽罩形状的特殊性，其边界层较薄，所以在不同的转速下，大部分空气流场都是相同的。由图5还可以看出，不同条件下水滴体积分数云图基本相同。这是因为旋转只在水滴到达边界层并非接近表面时才会影响水滴的运动。由于水滴的惯性，水滴受到剪切应力的影响较小，同时也没有足够的时间来改变运动。

不同转速下的水滴最大收集系数和碰撞极限如表3所示。最大收集系数为驻点对应的值，碰撞极限是沿曲线的长度。结果表明，随着转速的增大，最大收集系数降低，碰撞极限增大。但是差值非常小，可以忽略不计。因此，在大多数情况下，旋转对水滴撞击特性的影响可以忽略不计。

4 结论

本文利用单一旋转坐标系和欧拉方法，建立了水滴在旋转表面撞击特性的三维撞击模型。并利用该方法对航空发动机旋转整流帽罩进行了数值计算和分析。结果表明，随着马赫数的增加，水滴收集系数增加，撞击区域增大。水滴直径分布对撞击特性也有明显的影响，其结果比单一直径更准确全面。而整流帽罩的旋转对水滴撞击特性的影响则可以忽略不计。

参考文献

[1] Gent R W.TRAJICE2-a combined water droplet trajectory and ice accretion prediction program for aerofoils：RAE TR90054［R］.RAE TR90054，1990.

[2] Hedde T， Guffond D.ONER A three-dimensional-icing model［J］.AIAA Journal，1995，33（6）：1038-1045.

[3] Morency F， Beaugendre H， Habashi W.FENSAP-ICE：effect of ice shapes on 3D eulerian droplet impingement［C］//41st Aerospace Sciences Meeting and Exhibit.Reno：AIAA，2003：1223.

[4] Bidwell C.Ice particle transport analysis with phase change for the E3 turbofan engine using LEWICE3D［C］//4th AIAA Atmospheric and Space Environments Conference.New Orleans：AIAA，2012：2012-3037.

[5] Reid T， Baruzzi G， Ozcer I，et al.FENSAP-ICE simulation of icing on wind turbine blades，part 1：performance degradation［C］//51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.Grapevine ：AIAA，2013：750.

[6] Switchenko D， Habashi W， Reid T，et al.FENSAP-ICE simulation of complex wind turbine icing events，and comparison to observed performance data［C］//32nd ASME Wind Energy Symposium.National Harbor：AIAA，2014：1399.

[7] 赵秋月，董威，朱剑鋆.发动机旋转整流帽罩的水滴撞击特性分析［J］.燃气涡轮试验与研究，2011，24（4）：32-35.

[8] 王治国，娄德仓，郭文.发动机旋转表面水滴撞击特性数值研究［J］.燃气涡轮试验与研究，2013，26（1）：35-39.

[9] 易贤，王开春，马洪林，等.大型风力机结冰过程水滴收集率三维计算［J］.空气动力学学报，2013，31（6）：745-751.

[10] 吴孟龙，常士楠，冷梦尧，等.基于欧拉法模拟旋转帽罩水滴撞击特性［J］.北京航空航天大学学报，2014，40（9）：1263-1267.

[11] 王健，胡娅萍，吉洪湖，等.旋转整流罩积冰生长与脱落过程的实验［J］.航空动力学报，2014，29（6）：1352-1357.

[12] ANSYS Inc.ANSYS fluent 13.0 UDF manual［EB/OL］.（2010-08-01）［2024-06-10］.

[13] 申晓斌，林贵平，杨胜华.三维发动机进气道水滴撞击特性分析［J］.北京航空航天大学学报，2011，37（1）：1-5.

[14] 毛晓东，卜雪琴，赵国昌，等.基于UDS框架的水滴撞击特性数值计算方法［J］.沈阳航空航天大学学报，2016，33（1）：8-12.

[15] 杨胜华，林贵平，申晓斌.三维复杂表面水滴撞击特性计算［J］.航空动力学报，2010，25（2）：284-290.

基金资助

辽宁省“兴辽英才计划”项目(XLYC1802092)