含鸟撞变形叶片的压气机气动性能数值模拟

2022-07-05刘双丽王志强

航空发动机 2022年2期

刘双丽，陈伟，3，王志强，3，罗刚

（1.南京航空航天大学能源与动力学院，2.机械结构力学及控制国家重点实验室：南京 210016；3.辽宁省航空发动机冲击动力学重点实验室，沈阳 110015）

0 引言

航空发动机吸入飞鸟后，不但可能发生严重的机械损伤，同时极易导致整机性能恶化，推力降低甚至停车，危及飞行安全。为避免上述危险，各国适航管理机构颁布了适航规章，如美国FAR33部，欧洲的CSE800，中国的CCAR33部等，以法规形式规定了航空发动机吸入飞鸟后应具备的最低安全工作能力，如吸入中鸟后需要维持75%推力，吸入大型群鸟后需要维持50%推力等标志性要求。长期以来，围绕航空发动机吸鸟后风扇/压气机叶片的损伤、部件或整机的动态响应、适航符合性设计与验证等开展了大量研究，逐渐意识到发动机吸鸟后的结构损伤引发的气动性能变化，可能是导致发动机整机性能恶化和失去推力或工作能力的重要原因，因此针对变形受损叶片或受损部件的气动性能持续进行了相关研究。Imregun等采用稳态CFD分析方法进行了计算流体力学/有限元法（CFD/FEM）耦合计算，对人为假设的鸟撞击损伤叶片的气弹稳定性问题进行了研究；Kim等研究了70%转速下受到不同形式损伤的风扇转子叶片的气动稳定性；Bohari等通过在Rotor67转子叶片的前缘部分假设模拟鸟撞损伤变形，采用计算流体动力学方法获取了风扇在不同转速下的特性变化；Li等研究了Rotor67转子叶尖鸟撞损伤对旋转失速与喘振裕度的影响；Muir等对某涡扇发动机风扇在起飞阶段鸟撞的气动问题进行了研究；杨杰等基于SPH法鸟撞击发动机风扇数值模拟技术，应用NUMACA流体计算软件开展了受损叶片的流场数值计算；罗刚等建立了含有3联受损叶片的整级风扇CFD模型，对整级风扇气动性能变化进行了预估；陆嘉华等在平面叶栅试验验证的基础上，采用数值模拟方法研究了中鸟对典型风扇转子叶片造成损伤后导致的风扇气动性能变化。

综上所述，国外以高精度数值模拟方法为主，开展了较为深入的研究；国内开展了部分建模、数值模拟与平面叶栅试验验证工作，其他研究绝大部分仍停留在抗鸟撞结构分析与设计领域，对于发动机吸鸟后导致压气机转子部件的气动性能变化和整机性能降低过程缺乏了解，对发动机性能恶化的作用机理不够明确，在发动机吸鸟后气动性能变化的高精度虚拟分析手段亟需掌握。

本文针对含有鸟撞变形叶片风扇/压气机气动性能分析的需求，建立了NASA Rotor37转子在鸟撞前后的气动性能CFD分析模型，开展了数值模拟，采用试验数据进行了方法验证，并开展了对比分析与方法流程研究工作。

1 鸟撞变形叶片的几何建模

1.1 叶片模型选取

为研究压气机转子叶片受到鸟撞发生变形后其气动性能的变化，本文选择了NASA公开的Rotor37转子叶片作为研究对象，Rotor37转子是NASA Glenn Research Center于20世纪70年代设计的低展弦比的核心压气机跨声速第1级转子，其整级具有36片叶片，叶尖相对马赫数为1.48，转速为17188 r/min，该转子的性能参数在跨声速压气机中具有典型性，并有详细的试验数据，可作为数值模拟方法的有效验证依据。Rotor37转子单片叶片几何模型如图1所示。

图1 Rotor37转子单片叶片几何模型

1.2 叶片几何建模流程

采用基于迹线提取的半自动建模方法建立适用于气动性能分析的变形叶片几何输入模型。首先建立鸟撞叶片的有限元模型，在建模分网时排布叶片迹线的节点编号；随后在冲击动力学分析软件LS-DYNA中开展鸟撞叶片瞬态动力学分析，输出含有上述迹线节点坐标的文件，通过对节点编号进行编程排布，生成符合气动分析软件NUMECA适用的GEOMTURBO文件，导入NUMECA中生成几何模型。鸟撞变形叶片的几何建模流程如图2所示。

图2 鸟撞变形叶片的几何建模流程

1.3 鸟撞变形叶片的几何模型

建立了小鸟撞击含2联叶片的Rotor37转子叶片有限元模型，将Rotor37转子撞击前的2联叶片几何模型进行网格划分，选用Solid164实体单元，采用SPH方法进行鸟体建模，将鸟体设置为长径比为2的圆柱体模拟鸟，由于叶片尺寸较小，因此取鸟体质量为9 g。鸟体撞击位置设置在50%叶高，鸟体与叶片间采用点面侵彻接触方式，对叶片进行根部约束及旋转状态设置，鸟体以100 m/s的入射速度被100%转速旋转的叶片切割。在建模时，对叶片网格加密区域的迹线进行了节点编号，所建立的有限元模型如图3所示。

叶片材料采用PLASTIC_KINEMATIC模型，具体参数见表1。鸟体采用弹塑性水动力本构模型与多项式方程描述其大变形特征。鸟体材料参数见表2。

图3 NASA rotor37转子的鸟撞有限元模型

表1 叶片材料参数（TC4）

表2 鸟体材料参数

叶片受到鸟撞的数值模拟残余变形网格模型如图4所示。

图4 鸟撞叶片的数值模拟残余变形网格模型

从图中可见，模拟圆柱鸟体沿发动机进气方向撞击叶片，在撞击位置使叶片前缘发生了波浪形鼓包变形，其中受到鸟体撞击的叶片中间部位的残余变形幅值最大，对该模型的变形迹线节点进行编程排布，生成GEOMTURBO文件并导入NUMECA中建立含变形叶片的整级叶片几何模型，如图5所示。

图5 含变形叶片的整级叶片几何模型

从图中可见，模型包含了1片鸟撞后变形的叶片和35片完好叶片，构成了整级叶片气动分析用几何模型。

2 鸟撞变形叶片的气动性能数值模拟方法

2.1 含变形叶片的整级叶片流场计算模型的CFD网格划分

对整级叶片几何模型（图5）进行全3维的黏性流场计算用CFD网格划分，针对受损转子的网格划分采用NUMECA软件中的IGG-Autogrid5模块。叶片区域网格的拓扑结构采用O4H型，叶片上、下游网格为H型，叶尖径向间隙内网格为蝶形。第1层网格厚度为0.003 mm，径向设置了73个网格节点，叶尖径向间隙内设置了13个网格节点。该受损转子整级36片叶片通道的总网格量约为4166万。网格划分结果见表3并如图6所示。

表3 含变形叶片的Rotor37转子整级叶片的计算域网格参数

2.2 数值模拟方法

数值模拟和分析目的在于计算含变形叶片的压气机稳态气动性能，为保证计算的效率和稳定性，采用定常计算的方法，开展含变形叶片的全通道全3维流场数值模拟。受到鸟撞受损后的Rotor37转子整级全3维流场数值模拟采用NUMECA软件中的Fine/Turbo模块完成。采用时间推进法求解3维雷诺平均的NS方程，时间项采用4阶Runge-Kutta法迭代求解，空间离散采用2阶精度的中心差分格式进行，湍流模型选用S-A模型。Rotor37整级全3维流场定常数值模拟结果的可靠性采用文献[22]中的试验数据进行验证。

图6 受损Rotor37转子整级叶片流场计算网格

3 含鸟撞变形叶片的压气机气动性能数值模拟结果分析

3.1 Rotor37转子鸟撞前后气动特性对比

Rotor37转子受到鸟撞受损前后的气动性能对比如图7所示。图中黑色方点为文献[22]给出的试验测量得到的完好Rotor37转子特性，蓝线为计算得到的完好转子的气动特性，橙线为计算得到的受到鸟撞后转子的气动特性。

图7 Rotor37转子受到鸟撞受损前后气动性能对比

从图中可见，鸟撞导致叶片变形后，该转子的气动性能明显降低，在所有能稳定工作的流量状态下，受到鸟撞后转子的压比和效率都要明显小于鸟撞前的，相关参数变化结果见表4。

表4 鸟撞前后转子气动特征参数对比

从表中可见，数值模拟获得的转子受到鸟撞前的特性与试验结果比较接近，说明本文所采用的数值模拟方法是可靠的。而从鸟撞前后数值模拟获得的转子气动特性的对比可见，鸟撞后的转子气动性能较撞击前的虽然明显降低，但是鸟撞前后的特性变化规律一致，说明该数值模拟结果能够较合理地反映出含有变形叶片的转子气动性能影响的一般规律，说明本文针对受损叶片所采用的气动性能计算方法较为有效。

3.2 鸟撞后的Rotor37转子流场分析

3.2.1 堵塞和最大效率工况下的流场分析

在堵塞工况和最大效率工况下，Rotor37转子受到鸟撞后不同轴向截面相对马赫数分布分别如图8、9所示。图中的4个轴向截面分别位于转子上游、叶排中间、转子下游以及计算域出口。

图8 Rotor37转子受到鸟撞后不同轴向截面相对马赫数分布（堵塞工况）

图9 Rotor37转子受到鸟撞后不同轴向截面相对马赫数分布（最大效率工况）

从图8中可见，整体而言，在该状态下转子通道内的流动状态较好。在远离变形叶片的区域，流场结构较为一致，流动的周向不均匀性并不显著。但是在变形叶片区域，流场结构还是呈现出了明显不同。主要表现为在变形叶片的吸力面与相邻叶片压力面所形成的通道内存在大范围的低速区，其产生原因是叶片变形后，相应径向位置处的叶型几何进口角大幅增大，有的甚至超过90°（图6（d）的54%叶高位置），导致气流攻角也大幅增大，从而引发了受损叶片吸力面附面层的流动分离，堵塞了叶片通道，降低了整个转子的流动能力，进而使转子在堵塞工况下的流量减小，这与图7中的结果一致。进一步观察受损叶片通道中低速区的发展可见，在该状态下随着气流向下游的发展，在主流的掺混作用下，该低速区的影响逐渐减弱，在转子下游截面，其速度已得到明显提升。再观察计算域出口截面发现，出口截面的相对马赫数分布已比较均匀，上游低速区的影响已基本消失。

对比图8、9可见，在最大效率工况以及堵塞工况下，流场结构并没有明显变化，主要特点还是表现为在变形叶片附近的流动受到了叶片变形带来的影响；而在远离变形叶片的区域，各通道内流动的一致性较好。由于受损叶片的变形使来流攻角增大，导致气流在受损叶片吸力面附近形成了大范围的分离区。

在堵塞工况下不同S流面的相对马赫数分布如图10所示。结合图8、10可见，由于该受损叶片的变形位置主要在40%叶高以上，在叶根附近没有明显变形。从图10（a）中可见，在10%叶高位置，受损叶片通道的流动状态相对较好，未出现明显的附面层分离现象；而从图10（b）、（c）中可见，在54%叶高以及95%叶高处由于叶片变形严重，都出现了由附面层分离导致的低速流动区域。从图10中还可见，3个不同叶高截面上等熵马赫数分布由于受到受损叶片通道堵塞的影响，受损叶片前缘上游弓形激波的形态和位置与其它叶片也有明显区别，通道的堵塞作用类似于提高了激波后的反压，使得激波强度增大，激波后的马赫数出现了明显降低；该叶片弓形激波的位置相比于其它叶片前缘的激波更靠近上游；与受损叶片左侧相邻的约5个叶片通道（与转子旋转方向相反）的流动都受到了该叶片上游弓形激波的左支——外伸激波的影响，其上游的来流马赫数也相应偏低；由于受损叶片外伸激波减弱，在受损叶片左侧第6个叶片通道上游，激波后的马赫数迅速提高，形成了局部高马赫数区域；虽然受损叶片吸力面附面层的分离发生在通道上部，但是其形成的堵塞作用却是3维的，因为在10%叶高截面上也可以看到，其前缘的外伸激波强度较大，激波后的马赫数出现了明显降低，说明通道内附面层分离引起的堵塞作用是空间3维的。

图10 在堵塞工况下不同S1流面的相对马赫数分布

3.2.2 近失稳工况下的流场分析

Rotor37转子在受到鸟撞后在近失稳工况下不同轴向截面相对马赫数分布如图11所示。

Rotor37转子在受到鸟撞后，在近失稳工况下不同S流面相对马赫数分布、变形叶片区域和远离变形叶片区域相对马赫数为0.2的等值面分别如图12~14所示。

图11 Rotor37转子在受到鸟撞后近失稳工况下不同轴向截面相对马赫数分布

图12 Rotor37转子在受到鸟撞后在近失稳工况下变形叶片附近区域相对马赫数为0.2的等值面

与前2种工况相比，在近失稳工况下，Rotor37转子内的流动发生了较明显变化。主要表现为：在该工况下，变形叶片通道内的低速区范围明显增大，并且在其相邻通道内的近叶尖区域也出现了贯穿整个通道的低速区。从图12中的变形叶片区域的低马赫数等值面的位置和形状可见，2个通道的低速区在叶片下游连成了较大区域的低速团。从图13（c）中可见，出现这一现象的原因应该是，随着该转子的节流，转子下游的背压逐渐增大，转子上游的激波强度也逐渐增大，波后马赫数逐渐降低，攻角进一步增大，导致分离区也进一步增大。

图13 Rotor37转子在受到鸟撞后在近失稳工况下不同S1流面相对马赫数分布

图14 Rotor37转子在受到鸟撞后在近失稳工况下远离变形叶片区域相对马赫数为0.2的等值面

另外，从图11、13（c）以及图14中可见，在近失稳工况下，与变形叶片间隔4个叶片通道（与转子旋转方向相反）的叶尖区域也出现了大范围的低速区。该低速区占据的周向范围较大，约6个叶片通道。这可能是因为受到变形叶片上游外伸激波影响的区域来流马赫数较小，而该外伸激波影响开始消失的区域的流速又突然增大，导致这些位置处的叶片前缘激波强度增大，在叶尖泄漏流的共同影响下，这些叶片通道内的流动状态恶化，引起了流动的分离，形成了大范围的低速区。同样，这些区域的低速区也会引起流道的堵塞，其作用类似于受损叶片通道内的低速区，从而进一步导致在间隔了若干个叶片通道后，又形成了1个叶尖局部低速区。这些叶尖局部低速区与失速团的形式非常接近，其存在导致该转子很快进入失稳工况。此外，这些低速区的存在也导致总压损失和熵的增大，使得该转子的压比和效率明显降低，这与图7给出的该受损转子的性能变化情况吻合。

以上针对Rotor37转子受到鸟撞前后的气动性能以及3种典型工况下的流场结构分析表明，本文针对含鸟撞变形叶片的Rotor37转子所采用的全通道全3维黏性流场数值模拟方法可行，其结果合理，利用该方法可以较准确地计算出含变形叶片转子内部的流场细节，并获得其气动性能的变化特性。