航空发动机吞冰损伤一体化数值模拟平台搭建与验证

沈安庆 ，徐惊雷 ，夏全忠 ，夹福年 ，姚艳玲

（1.南京航空航天大学能源与动力学院，南京 210016；2.中国航发四川燃气涡轮研究院，四川绵阳 621010）

0 引言

航空发动机进气系统吞冰可能造成发动机结构失效、性能降低，严重影响飞行安全[1]。发动机整机吸冰试验是适航取证必须通过的验证项目[2]，代价高昂，因此对发动机风扇叶片冰撞击问题开展显式动力学数值仿真方法研究十分必要。

国外20世纪70年代就开始对飞行器冰撞击问题开展研究。1973年，Hayduk 等[3-4]对厚度为1.0～1.6 mm 的2024-T3 铝合金薄板试验件进行冰撞击试验，以冰块的直径、速度和靶板厚度为变量，使用有限差分程序对薄板的冰块冲击进行了计算，靶板变形计算结果与试验较吻合；1991年，Reddy 等[5]建立了风扇叶片的有限元模型，将冰冲击载荷等效为分布在叶片前缘节点上的作用力，利用模态叠加法对冰撞击风扇叶片进行了瞬态分析，研究了冲击位置、叶片转速等对叶片变形的影响；2000年，Kim 等[6-7]首次基于通用显式动力学有限元软件DYNA-3D 建立了冰的有限元模型，为人工冰块建立了带失效准则的弹塑性本构材料模型，在速度为60～190 m/s 条件下开展了冰撞刚性靶板以及冰撞复合材料平板的试验，标定了数值仿真模型，但是该模型未考虑应变率对冰力学性能的影响；2005年，Anghileri 等[8]对采用拉格朗日法（Lagrange method），任意拉格朗日欧拉法（Arbitrary Lagrange Euler，ALE）和光滑粒子流体动力学（ Smooth Particle Hydrodynamic，SPH）法3 种数值方法对冰块进行了建模，并对冰块撞击进气道前缘进行了数值仿真，表明SPH 数值仿真模型计算速度最快，计算结果与试验吻合最好，而且能够反映冲击后冰块的破碎行为。自2000年以来，中国学者对飞行器适航方面的问题越来越重视，很多高校都开展了结冰与冰撞击方面的工作。汪洋等[9-10]利用分离式Hopkinson 压杆对温度为-25、-10 ℃的冰开展应变率为500～2000 s-1的动态压缩试验，研究了结冰温度、应变率对冰块动态力学性能的影响；孟卓等[11]比较了拉格朗日方法、ALE 方法和SPH 方法对冰块进行数值模拟的吻合程度，并模拟了冰块以185 m/s 的速度撞击飞机发动机进气道的过程；李静等[12]利用SPH 方法建立冰块模型，利用Carney 等[13]提出的冰材料参数，对厚度为0.9～1.6 mm 的2024 铝合金薄板进行了多个随机分布冰块冲击靶板的数值仿真研究，表明冰块冲击损伤具有叠加效应。目前已开展的研究多针对蒙皮、雷达整流罩等飞机机身部件的冰撞击问题，对航空发动机风扇叶片冰撞击问题涉及较少，并且缺乏工程化解决此类问题的专业软件系统。

本文通过对冰姿态流场以及冰撞击叶片动力学进行仿真研究，建立了吞冰流场模型和叶片损伤模型，并基于LS-DYNA、Fluent 的二次开发，形成了航空发动机叶片吞冰损伤的快速分析软件系统。

1 冰姿态流场分析

1.1 重叠网格和6自由度流场计算理论方法

本文基于冰块姿态流场计算和SPH 碰撞计算理论开展发动机冰块吞入过程及撞击特性的数值仿真分析研究。冰块吸入过程采用N-S 方程和刚体6 自由度运动方程耦合求解[14]。由于冰块和流场计算域的尺度差异较大，为了节约计算成本，采用重叠网格技术实现。重叠网格技术是CFD 计算刚体运动中比较成熟的一种算法，和网格变形及重构方法相比，能够实现边界层加密、六面体、混合网格等多种类型的计算，具有更好的稳健性和精度，运动部件重叠网格构建方法如图1所示。

图1 运动部件重叠网格构建方法

为验证冰块6 自由度流场计算精度，采用文献[15]的标准投弹模型验证条件：马赫数为1.2，静温为216.65 K，静压为20646 Pa，飞行高度为11600 m。基于重叠网格技术可以精确地模拟弹投放过程，计算和试验所得到的姿态高度吻合。弹6 自由度重叠网格姿态和试验数据对比如图2所示。

图2 弹6自由度重叠网格姿态和试验数据对比

在发动机吸入冰块的运动过程中，在从中间位置释放冰块的风洞试验结果中，发现冰块在运动一段距离后总会发生偏移而打到壁面某个位置。

本文建立了验证的冰块6 自由度运动模型，冰块在风洞中不同时刻的运动轨迹如图3所示。

图3 冰块在风洞中不同时刻的运动轨迹

冰块在风洞中的运动轨迹和常规预想的有较大差异，原因在于冰块的形状不规则，导致其在旋转运动过程中各表面压力不均衡，从而形成马格努斯效应，即冰块往一侧倾斜运动，而不是对称运动，这种现象在实际风洞试验中也存在。为减少计算的量和难度，发动机吸入模型采用了动量源项的方法进行模拟。基于Fluent中动量源项加入虚拟风扇模型，真实模拟轴流风扇效应。

风扇动量源项的轴向分量

径向分量（风扇动量源项的径向分量）

剪切分量

式中：ΔP(Q)为给定流量的压升；h为轴向的叶片厚度，即风扇叶片对应半径叶型弦长；Wfan为风扇功率；Ωop为风扇角速度；Vφ为局部剪切速度，即风扇对应半径位置处当地切向速度分量；ρ为流体密度；r为风扇半径，Rh为风扇轮毂半径，Rip风扇等效半径，Rt为风扇叶尖半径；c1为与几何相关的常数。

可以通过建立真实的发动机多级风扇模型验证上述方法。2 种模型模拟的唇口气流速度分布对比如图4 所示。从图中可见，简化动量源项模型和实际吸气模型的流场非常接近。

图4 2种模型模拟的唇口气流速度分布对比

在实际的发动机进气模型中，风扇导致的进气旋流效应在靠近唇口处不太明显，可忽略旋流效应导致的冰块在唇口附近的运动偏转。

1.2 吞冰流场仿真建模过程

吞冰过程的流场建模基于Fluent软件平台完成，为了精细模拟冰块的运动轨迹、欧拉角姿态及其撞击叶片瞬间的速度、位置，为后续SPH 碰撞计算提供初始边界输入，评估冰块撞壁及吸入情况。

对风洞中冰块轨迹进行分析，搭建了风洞稳压室、风洞管道、发动机进口风扇、冰块几何模型，同时建立重叠网格计算模型，网格类型采用四面体混合网格，数量为35 万。在冰块表面划分贴体边界层网格，采用中等网格密度。冰块在翻转过程中会有较大的分离流动产生，采用剪切应力输运（Shear Stress Transport，SST）湍流模型可有效模拟中度分离流动，数值仿真分析对象总模型及冰块边界层切面网格如图5所示。

图5 数值仿真分析对象总模型及冰块边界层切面网格

风扇模型基于实际参数给定，根据风扇设计部门给定的P-Q 曲线（即风扇流量-压降曲线）给定压降。基于Fluent 软件进行了冰块撞壁自动停止计算和反弹算法开发，程序会根据最小距离自动判断撞壁位置并停止计算，同时可以考虑冰块的碰撞反弹效应。

吸冰过程流场的计算需考虑冰块不同尺寸、不同放置方式、与风扇入口不同距离等情况。当投放位置确定，冰块尺寸小于0.3 倍原始尺寸（冰块原始尺寸为：长0.18D，宽0.10D，厚0.01D，D为风洞直径）时，才能顺利进入叶片区域，否则冰块会碰壁；冰块竖向放置时，1倍与0.5倍冰块均会碰壁。对确定投放位置，冰块与风扇入口不同距离条件下，2种情况下冰块均碰壁。

2 冰撞击叶片动力学分析

2.1 SPH建模方法

SPH 是核函数近似的配点型无网格方法，将连续的介质离散为质点，质量、速度、温度和压力等物理场分布在质点上，通过对质点动力学方程的求解获得质点的运动轨迹以及质点之间的相互作用，得到所有质点组成系统的力学行为。

粒子的近似函数为

式中：f（x）为任意空间上的物理量；h为光滑长度，随着时间和空间变化；W为核函数，其定义为

式中：θ（x）为辅助函数；W（x，h）为具有中心峰值的函数。

通常核函数采用立方B 样条函数，此时辅助函数θ（x）的形式为

式中：C为归一化常量，与空间的维数有关。

冰撞击问题中，SPH 方法采用点-面接触模拟流体与固体的相互作用，规避了ALE 方法流固耦合界面问题。采用SPH 方法建立冰块数值仿真模型能准确模拟冰块的破碎过程及冰碎片与叶片的相互作用。

2.2 冰本构模型及材料参数标定

冰是一种具有脆性失效的线弹性材料。冰的力学特性具有应变率相关性，随着应变率的增加，冰的抗压强度随之增大。冰的抗压强度远远大于抗拉强度，在压缩载荷下，冰破碎后仍存在剩余强度。冰块材料模型使用LS-DYN 软件的MAT155（关键字MAT_PLASTICITY_COMP-RESSION_

TENSION_EOS）模型，该模型考虑了弹塑性材料的拉压非对称性。状态方程采用关键字EOS_TABULATED_COMPACTION 定义，冰块受到的压缩载荷为

式中：p为冰块受到的压缩载荷；C(εv)为基于体积应变的冷压力；γT(εv)E为取决与体积应变和温度的热压力，其中E为热比内能，T为温度，γ为Gruneisen系数。

本文对直径为25.3 mm 的冰块以192 m/s 的速度撞击铝合金靶板工况进行计算，靶板边界处设置位移约束。冰块撞击靶板不同时刻破碎过程和飞溅分布情况如图6、7 所示，铝板中心位置处最大位移为11.98 mm 铝板位移变形从四周边界向平板中心不断增加。

图6 SPH冰块及铝板模型变形过程

图7 铝合金靶板中心位置处位移时间历程曲线

针对不同速度冰撞击工况进行仿真模型标定，得到适用于航空发动机风扇叶片冰撞击数值仿真的冰块本构参数和仿真建模方法。冰块撞击铝合金靶板的数值模拟结果与文献[3]的试验结果对比见表1。从表中可见，数值模拟结果与试验结果的误差在8%左右，结果基本吻合。通过试验验证了数值模拟结果的准确性，可用显式动力学方法进行叶片抗冰块冲击分析。

表1 冰块撞击铝合金靶板数值模拟结果与文献[3]中试验结果对比

2.3 冰撞击发动机风扇叶片数值仿真分析

对直径为25.3 mm 的冰块以192 m/s 的速度撞击铝合金叶片的情况进行计算，在冰块撞击风扇叶片的过程中，自第0.1～0.4 ms 时刻冰块破碎过程、叶片动态变形的SPH 模型模拟结果如图8 所示，在第0.6 ms时刻叶片根部等效应力达到487 MPa，叶片等效应力分布如图9所示。

图8 冰块撞击叶片过程SPH模型模拟结果

图9 叶片等效应力分布

冲击位置为75%叶高处，叶轮安装孔处设置位移约束，在本构模型仿真过程中，冰块的初始撞击角度为0°，在流动轨迹+撞击的一体化仿真计算中，冰块初始撞击角度通过流场计算结果传递给结构计算。冰块采用SPH 模型模拟，冰的材料本构和建模方法与前述靶板冰撞击情况相同。

叶片3 个测点处的塑性应变历程曲线如图10 所示。从图中可见，在第0.6 ms后塑性应变增加缓慢，可通过冰撞击叶片试验获得叶片冲击区域的动态应变曲线校核数值仿真结果。

图10 叶片塑性应变历程曲线

3 航空发动机风扇叶片冰撞击评估系统开发

吞冰软件系统基于LS-DYNA、Fluent 二次开发，集成冰块运动、撞击过程的仿真流程和仿真规范。3维虚拟仿真系统可快速实现发动机吞冰过程的高精度模拟，实现冰块运动姿态、冰块撞击叶片破碎过程及轨迹预测的一体化仿真流程，集成系统底层逻辑原理如图11 所示。软件可快速实现几何建模及导入、网格智能划分、物理模型及边界条件设定、仿真自动求解、报告自动生成等一系列复杂仿真流程的智能实现，规范化流程大幅度降低了复杂问题的仿真门槛，提高了技术人员进行仿真计算的效率。

图11 集成系统底层逻辑原理

冰姿态流场计算模块主要包含流场几何模型创建、网格生成、冰块材料设置、边界条件定义、求解参数定义、求解、结果输出等功能。流场计算界面如图12所示，图中给出了流场非定常求解参数变化历程。

图12 流场计算界面

在“冰块”面板中选择冰块几何形状，包括“球形”和“方形”2 种。点击“冰块尺寸设置”按钮，定义冰块初始重心坐标位置及尺寸。在“叶片”面板中需定义叶轮厚度、叶根半径、整流锥角等。完成上述参数设置后，点击“流体域几何生成”按钮，系统将自动完成冰姿态流场几何模型的创建。

在网格参数设置面板中，对进气道和冰块网格参数（包含面网格单元尺寸、边界层首层、边界层增长率、边界层层数）进行定义，点击“网格生成查看”按钮，系统将自动调用网格处理软件对几何模型进行网格划分。

在边界条件设置中，定义仿真试验环境压力、入口边界条件、叶片边界条件及彻体力加速度分量等。入口边界条件需定义入口初始流量、最终流量、流量变化时间、入口温度等；叶片边界条件需定义叶片转速、剪切拐点、压升（降）等。

计算过程为先进行稳态计算，为后续瞬态计算提供稳态初场。在计算设置面板中，定义初始稳态迭代步，瞬态计算的库朗特数、时间步长、时间步数、子循环迭代等。用户可定义自动保存时间步、计算并行线程数。点击“冰块运动求解计算”按钮，系统将自动调用Fluent对冰姿态流场进行仿真求解计算。

冰块撞击叶片计算模块界面如图13 所示，主要包含参数设置、SPH 模型建立、叶片网格划分、冰块撞击叶片求解、仿真报告生成等功能。

在初始条件设置中，用户可以选择从流场的计算结果继承冰块速度，也可以自定义冰块各方向的速度分量。上述参数设置后，点击“SPH 模型生成”按钮，系统将自动创建冰块SPH模型，点击“叶片网格生成”按钮，系统自动生成叶片结构网格，点击“冰块撞击叶片计算”按钮，系统将自动调用LS-DYNA对冰块撞击叶片过程进行求解，在图形窗口中显示冰块撞击叶片的破碎过程；求解完成后，点击“仿真报告生成”按钮，系统将自动调用Word 生成仿真报告。

4 结论

（1）本文研究了发动机进气吞冰过程中冰体6 自由度运动姿态和轨迹的流场模拟方法，验证了冰体本构模型参数，准确地模拟了冰体在撞击叶片过程中的损伤效应。

（2）基于吞冰流场计算和冰撞击过程计算一体化耦合过程建立了发动机叶片吞冰损伤的快速分析方法和流程。

（3）基于该方法可以有效地解决吞冰损伤复杂过程有限元模型的自动生成问题，极大地提高了分析效率。

（4）完成了航空发动机风扇叶片吞冰损伤一体化数值模拟平台搭建，通过多工况数值分析验证了数值模拟平台计算能力的稳定性和预测分析能力。