砂尘对压气机叶片的侵蚀及性能影响

2023-11-27杨天南郑培英

航空发动机 2023年5期

杨天南，张轲，郑培英，张海

（1.海装沈阳局驻沈阳地区第二军事代表室，沈阳 110043；2.中国航发沈阳发动机研究所，沈阳 110015；3.哈尔滨工程大学动力与能源工程学院，哈尔滨 150001）

0 引言

飞机在起飞和降落的过程中会面临多砂环境，此时由风扬起的尘土、砂粒等异物颗粒会随空气被吸入航空发动机内，这些异物会在气流曳力的作用下不断加速，对发动机部件产生磨损，导致其性能降低甚至引起发动机故障，降低其运行的可靠性和使用寿命[1]。在实际应用中，由于砂尘侵蚀和沉积造成的事故和损失并不在少数。运输机或直升机在起飞、低速滑行或着陆时产生的漩涡和低速推力反向射流，会将砂粒和灰尘带入发动机，而恶劣的工作环境，如沙尘暴、火山喷发产生的火山灰、沙漠地带、浓烟浓雾等，会使砂粒冲蚀和沉积的影响更加明显，甚至影响结构的完整性[2]。1982 年，一架波音747 客机在飞行过程中驶入因火山喷发产生的火山灰云后，发生了引擎熄火的严重事故[3]；2000 年，火山灰云同样造成了飞机涡轮叶片前缘侵蚀[4]；同时，恶劣的砂尘环境导致军用直升机发动机耗油率和维修成本大幅提高，发动机实际使用寿命不足设计值的1/8，使用环境是造成直升机在沙漠地带受损严重的主要原因之一[5-6]。

Brun 等[7]的研究表明，直径在10 μm 以上的砂粒就会对压气机叶片造成侵蚀；傅国如[8]的研究表明，直径超过5 mm 的砂粒可将压气机第1 级叶片的前缘打出裂口；李钊等[9-11]在对压气机叶片断裂或失效的故障分析中发现，部分叶片会由于砂粒的点蚀或侵蚀作用在压气机叶片表面形成疲劳点，降低叶片疲劳强度，从而在点蚀处萌生裂纹、扩展直至叶片断裂。由此可见，砂粒粒径与其侵蚀作用有着明显关联，发动机吞砂防护是一项复杂的工程，就需要详细分析发动机中砂粒运动规律及砂粒磨蚀原理。针对砂尘颗粒对叶片的撞击，Balen等[12]利用1台单级轴流压缩机开展了试验研究，吞砂的过程持续了605 s，共计吸入25 kg砂尘，表明吞砂后压气机载荷系数减小3.5%，效率降低4%；Ghenaiet 等[13]对轴流风机进行了9 h 的连续试验，发现吸入砂粒9 h 后，压气机的绝热效率、总压比和喘振裕度分别降低了7.1%、9.5%和55%；Bammert 等[14]利用叶栅风洞对叶片进行了吞砂试验，发现砂粒的吸入会导致压气机工作点前移。

砂粒在发动机中的运动规律复杂，伴随着气流与砂粒相互作用和砂粒与壁面的碰撞作用，是人们无法从理论出发得出1 套准确可行的数学模型。目前，两相流数值模拟技术已逐渐成熟，两相流3 维计算已在工程上得到广泛应用。本文以Rotor 37 为仿真对象，在不同颗粒质量流量、不同粒径、不同转速下进行数值仿真，分析颗粒在不同吞砂条件下的动力学特性，并对吞砂过程中动叶的性能变化进行分析。

1 数值方法

航空发动机压缩部件吞砂是一种稀疏颗粒气固两相流动，砂粒在气流中运动受到气动力、惯性力等力的作用运动，遇到壁面发生碰撞反弹。在数值模拟中，一般将砂粒视为离散相，将气流视为连续相，建立各相的运动方程及相与相间的耦合作用是数值模拟的关键。颗粒流动数值模拟方法主要有2 种，欧拉-欧拉法和欧拉朗格朗日法，2 种方法在颗粒流模拟中都应用广泛，各有优缺点，2 种方法的选择主要是基于两相流动状态和模拟需求，下面对2 种方法的异同和选取依据进行详细说明。

1.1 欧拉-欧拉法和欧拉-拉格朗日法

欧拉-欧拉法和欧拉-拉格朗日法都是基于经典的连续介质力学方法，在计算不同类型的颗粒流动方面有着各自的特点和优势。欧拉-欧拉法是将离散的颗粒相和连续的流体相都看作是连续介质，同时在欧拉坐标系中考虑离散相和连续相的运动，对每一相都建立连续方程、动量方程和能量方程，来考察空间各点上的两相物理量及其变化；欧拉-拉格朗日法是将离散相和连续相分别用不同的方程来控制，连续相的介质的运动由经典的Navier-Stokes方程控制，捕捉每一个空间节点上连续相的物理量变化，而离散相的运动则由独立的动量方程控制，对离散相颗粒在空间中的运动轨迹进行追踪，一般将连续相计算方法称为欧拉法，将对离散相进行实时追踪的计算方法称为拉格朗日法。欧拉法和拉格朗日法分别如图1、2所示。

图1 欧拉法：利用空间固定点跟踪元素

图2 拉格朗日法：对粒子进行实时追踪

在发动机压缩部件吞砂过程中，砂粒受到气动力、惯性力等力的作用发生运动，遇到壁面发生碰撞反弹，离散相与连续相间存在动量和能量交换。《军用直升机防砂尘要求》（GJB 1171）[15]和《航空涡喷涡扇发动机吞砂试验要求》（GJB 2026）[16]规定，发动机吞砂试验砂粒最大浓度为53 mg/m3，质量分数约为0.00004%（砂粒/空气），换算砂粒的体积分数则更小。各地的砂尘粒度也不一致[17]，因此在计算中需要酌情考虑。

采用不同方法来计算气流连续相和颗粒离散相的欧拉-拉格朗日法可以较好地模拟出砂粒运动特性以及两相间的耦合作用，较低的砂粒体积分数也满足欧拉-拉格朗日方法的使用要求，综合考虑，本文采用欧拉-拉格朗日法来计算压缩部件吞砂过程中的气固两相流动，其中通过动量和能量源来建立连续相与离散相之间的耦合作用。

1.2 欧拉-拉格朗日法控制方程

1.2.1 连续相控制方程

（1）质量守恒方程

式中：ρ为空气工质的密度；t为时间；u→为混合工质连续相的速度。

（2）动量守恒方程

式中：p为作用在控制体上的表面压力为控制体所受到的体积力，包括重力、惯性力及砂粒与连续相间的相互作用力等力为连续相流体粘性剪切应力张量。

（3）能量守恒方程

式中：ht为连续相总焓；λ为砂粒与气流间的热传导率；Sh为热量源，为砂粒与连续相气流间传递热量，离散项与连续相间保持能量平衡。

1.2.2 离散相控制方程

发动机吞砂过程中，砂粒作为稀疏的离散相分布在气流中，与气流间存在一定的滑移速度，砂粒的运动会受到多种作用力控制，包括自身重力、虚拟质量力，以及受气流作用的气动曳力和压力梯度力；在旋转运动过程中，砂粒还受到离心力和科氏力的作用。离散相颗粒的运动与受力如图3 所示。另外，砂粒通过对流换热与气流间产生热量交换。

图3 离散相颗粒运动与受力

（1）砂粒运动方程

其中重力、虚拟质量力和压力梯度力与气动曳力和惯性力相比是小量（约相差三个数量级）。忽略重力、虚拟质量力和压力梯度力，砂粒主要受气动曳力和惯性力作用，在拉格朗日坐标系下，砂粒动量方程为

其中，气动曳力

式中：CD为阻力系数，up为球形颗粒在静止流体中的运动速度，Ap和dp分别为球形颗粒面积和直径，ρˉ为流体密度。

（2）砂粒热量传递方程

式中：Ap为砂粒表面积；Cs为砂粒定压比热容；hc为气流与砂粒颗粒尖的对流换热系数；T和Tp分别为气流和砂粒表面的温度。

砂粒在运动过程中受到各种气流作用力影响，作用力也会反作用于气流，砂粒与气流间也存在一定热量交换。欧拉-拉格朗日法以控制体为单元，将离散相砂粒产生的动量和能量变化量作为源项引入连续相控制方程，通过迭代计算保证两相之间作用源的平衡，从而建立离散相和连续相间的耦合作用。

1.3 砂粒侵蚀模型

砂粒对壁面的侵蚀程度一般用侵蚀速率来衡量，侵蚀速率定义为壁面材料在单位时间上损失的质量（单位：kg/m2/s），在数值计算中，通过计算每个颗粒对壁面的累计损伤来计算侵蚀速率

式中：Af为壁面单元面积̇冲击壁面颗粒质量流量；er为侵蚀率，侵蚀速率与颗粒的运动状态有关。

关于侵蚀建模的研究多偏向于基于试验的工程应用数值解，通过理论分析和试验结果归纳出不同颗粒状态下的侵蚀率

式中：β0为侵蚀率最大时碰撞角度，当碰撞角度为90°时，a取值4；K1、K2、K3为与颗粒类型和碰撞壁面材料有关的系数，通过试验测试获得，石英颗粒撞击不同材料（铝合金、不锈钢、钛合金）时，K1、K2、K3取值见表1。

表1 不同材料K1、K2、K3取值

2 模型与边界条件

采用计算流体力学与离散元法（Computational Fluid Dynamics-Discrete Element Method，CFD-DEM）耦合方法模拟压气机在不同吞砂条件下各部件的磨蚀情况，总体思路是使用NUMECA AutoGrid 对Rotor37 进行网格划分，并进行网格无关性验证，随后根据Rotor37 的试验数据进行数值模拟的校核，选取合适的边界条件进行耦合仿真，最后通过改变吞砂条件，综合比较不同吞砂条件下的侵蚀变化。

Rotor 37 是典型的跨声速轴流压气机模型之一，其具体设计参数见表2。Rotor 37 网格划分形式及拓扑结构如图4 所示。通过无关性验证，本文对于Rotor 37 单通道的网格数选用80 万，整周共有36个叶片，则整周的计算网格数为2880万。

表2 Rotor 37部分设计参数

图4 Rotor 37网格划分形式及拓扑结构

对于单相流场计算，使用FLUENT 进行仿真，选用理想气体作为流体工质，湍流模型选用标准k-e湍流模型，近壁面函数选用标准壁面函数，进口边界条件选用压力进口，进口总压给定101325 Pa，进口总温给定288.15 K，湍流粘度比给定3%；出口边界条件选用压力出口，静压按径向平衡压力分布给定9000 Pa，出口温度给定300 K，轮毂、叶片和机匣给定光滑壁面，选用多重坐标系（Moving Reference Frames，MRF）法定义压气机的旋转，其中轮毂和叶片根据坐标及流道方向，给定转速值为-17188 r/min，机匣设为静止，求解方法选用Couple 算法，残差量级给定10-5。稳态计算收敛后，获得该工况下的流量、效率、压比等参数。

由于CFD-DEM 耦合仿真为瞬态计算，与稳态计算不同的是，瞬态计算的求解速度和准确性会受到初场的影响，合理的初场会提高瞬态计算过程中的收敛速度和稳定性，因此为了节省计算时间并保证瞬态计算结果的准确，以近设计点工况的稳态仿真结果作为瞬态仿真的初场。沙砾粒径选取为平均150 μm，质量比为1%，颗粒速度为170 m/s。

3 数值结果分析

3.1 颗粒轨迹及分布

砂粒的运动情况直接决定着撞击发生的位置，进而影响到磨损及侵蚀出现的位置及程度。因此首先要对颗粒的运动情况展开分析。整周内颗粒轨迹如图5 所示。在入口处，绝大部分的颗粒的运动方向和速度与气流保持一致，随气流进入压气机后撞击叶片，轨迹发生偏移，偏移方向与叶片的转动方向一致，由于在碰撞过程中颗粒会受到来自叶片的作用力，颗粒在改变运动方向的同时，移动速度也会略有提升。

图5 整周内颗粒轨迹

在某稳定运行时刻（t=0.005 s）的颗粒分布如图6所示，其中以颗粒速度作为颗粒的着色变量。从图中可见，颗粒在尚未与压气机叶片接触前，其在空气中是均匀分布的，当颗粒随气流流过叶片前缘后，其速度及运动方向均发生改变。与此同时，在前一片相邻叶片的叶背处，颗粒的密度降低，从图中还可见明显的颗粒空隙，产生这种现象的原因主要是颗粒与叶片之间的速度方向不同。颗粒的初始速度是沿压气机轴线的进气方向，而叶片则是围绕轴线做旋转运动，颗粒和叶片之间存在相对运动，由此颗粒与叶片便产生了速度夹角，发生颗粒与叶背处的分离，向下游相邻叶片的叶盆处移动聚集，使得叶背处颗粒密度降低。颗粒在叶片分流的作用下，从入口处的均匀分布变成了相对集中的颗粒束，在后续向叶片出口移动的过程中，由于气流尾迹的扰动作用，颗粒束之间再次出现掺混，粒子分布趋于平均。

图6 在第0.005 s时刻的颗粒分布

不同叶高处的颗粒质量浓度如图7 所示。颗粒质量浓度数值越大，说明此时颗粒聚集越明显，而着色区域越大，说明颗粒的分布越广泛。相较于叶片的吸力面，压力面处的颗粒体积分数更大，颗粒分布范围更广。在叶片尾缘到出口位置，随着叶高逐渐增大，这部分流场中的颗粒体积分数逐渐增加，颗粒分布范围也在扩大，相比于在进口段各截面的均匀分布，颗粒逐渐远离压气机叶根和轮毂部分，向叶顶和机匣处移动，这是由于压气机作为旋转机械，流场内会存在对颗粒的离心力，使得颗粒轨迹发生偏移，随着气流向叶顶和机匣处移动。

图7 不同叶高处的颗粒质量浓度

3.2 磨损情况

颗粒的直径除了对侵蚀量产生影响，也会改变粒子轨迹，使得侵蚀位置发生变化，为了探究粒径对磨损的影响和变化规律，本节在100%转速、砂粒质量流量分数为1%的近设计点工况下进行不同粒径的吞砂仿真模拟，选用的粒子直径包括100、150和300 μm。砂粒撞击下的磨损量如图8 所示。从图中可见粒径对磨损分布及磨损量的影响。叶片整体的磨损量随颗粒直径的增大而增大，小颗粒产生的磨损相对均匀，随着颗粒直径的增加，叶片出现了类似于“斑块状”的磨损区域。随着粒子直径的增加，同等吞砂质量流量下，颗粒的数量会减少，使得气流中所含颗粒的密度减小，出现侵蚀分布不均的现象。

图8 砂粒撞击下的磨损量

不同砂粒质量分数对叶片产生的侵蚀分布如图9 所示，根据云图显示，吞砂量的改变基本不会影响侵蚀区域的形态，只会改变磨损量的大小，在相同粒径下，磨损量随颗粒质量分数的增加而增加。

图9 不同砂粒质量分数对叶片产生的侵蚀分布

3.3 性能影响

本节主要通过改变吞砂条件，探究压气机的性能变化规律。颗粒质量分数为1%时，不同吞砂粒径下的压比特性如图10 所示。将仿真流量按设计点流量20.19 kg/s 进行归一化处理，首先是流量-压比特性线，从图中可见，相对于纯空气工况的特性线，吞砂时各背压对应的流量都在减小，压比特性线整体向左下方移动，在砂粒质量分数相同时，粒径越小，对特性线的影响越明显，向左下方移动的幅度越大。与纯空气工况相比，吞砂会使压气机的增压能力减弱，对于直径为20、150、300 μm 的吞砂工况，其增压范围分别缩小了3.83%、2.13%和1.01%，也说明了小粒径颗粒对压气机压比的影响更为明显。

图10 不同吞砂粒径下的压比特性

不同粒径下的压气机的温比、效率特性线如图11 所示，效率受到压比和温比的共同影响，且温比的影响更为明显。从温比图中可见，在同一背压下，吞砂直径为20 μm 时，其温比在小流量时较纯空气工况略微降低，在大流量时略微增加，温比的变化使得20 μm 工况下的效率较纯空气工况先增大后减小，而对于150、300 μm 工况，同背压下的温比较纯空气工况略微增加，因此使得吞砂工况的效率降低。对于20 μm的小粒径，相比效率，其对压气机流量的影响更为明显，使得效率特性线左移明显，而如150、300 μm 的大颗粒，其对效率的影响比小直径颗粒更为显著，此时效率特性线整体向左下方移动。

图11 不同粒径下的压气机的温比、效率特性线

在90%转速下的压比和效率的特性线如图12所示。模拟工况为1%吞砂量，粒径分别为20、150 和300 μm，流量系数通过各点的仿真流量与效率峰值点对应的流量归一化处理得到。在90%转速下，吞砂依旧造成了性能恶化，不论是压比特性线还是效率特性线，都存在不同程度的左移。在压比图中，小粒径（20 μm）的吞砂工况特性线左移幅度最大，即颗粒质量分数相同时，颗粒直径越小，流量和压比下降越明显。在效率图中，不同吞砂工况的效率特性线均不同程度地向左下方偏移，同样地，20 μm 工况产生的偏移程度最为明显，在小流量的工况时，20 μm 工况的效率线略有左移，这一点与设计转速时相同。