风扇叶片飞失显式动力学仿真网格与时间步长研究

2021-12-15李宏新冯国全

航空发动机 2021年5期

徐雪，李宏新，冯国全，3

（1.中国航发沈阳发动机研究所，沈阳 110015；2.中国航空发动机集团有限公司，北京 100097；3.辽宁省航空发动机冲击力学重点试验室，沈阳 110015）

0 引言

风扇叶片飞失（Fan Blade Out，FBO）后发动机动力学行为及其影响的设计技术是民用航空发动机安全性领域最关键的技术壁垒之一，是各适航体系都要求的重要试验[1-3]。近年来，随着显式动力学计算软件逐步成熟，凭借其对撞击、碰摩和断裂等典型高度瞬态、高度非线性物理过程的有效模拟，成为FBO 设计和仿真工作的重要方法[4]。

在该领域，Nicolas 等[5]开展了FBO 过程仿真与整机试验测试结果的对比分析，提出风扇转子和静子的模型简化方法；Robin[6]进行了整机FBO仿真分析中风扇转子叶尖严重碰摩的瞬态动力学行为计算，描述了叶尖碰摩的微观物理过程；Sinha[7-8]系统地发展了考虑碰摩的转子动力学行为计算理论，并给出了典型算例的仿真分析结果以及与整机试验数据的对比；马辉等[9]开展了旋转叶片-机匣系统碰摩动力学研究，提出转子盘、叶片和机匣碰摩系统的动力学简化分析方法；刘璐璐等[10]开展了缩比的FBO 整机模拟试验与仿真研究，获得安装结构的外传振动特点。由于显式动力学的算法特点，其网格尺度和时间步长设置在分析中具有关键作用，不仅影响计算精度和资源的需求，还直接影响计算的收敛性和物理过程时序，因此网格无关性和时间步长无关性研究对于瞬态显式动力学分析结果具有决定性影响，但是在以往的研究中对于这方面的分析较少。

本文基于商用有限元显式动力学计算软件，开展了FBO 过程中风扇叶片的网格与时间步长无关性研究。

1 显式动力学网格与时间步长的理论关系

1.1 计算方法和主要影响因素

显式动力学分析的基本动力学方程[11]为

式中：[M](t)为系统质量矩阵；[C+CT+G](t)为阻尼矩阵，其中C为系统阻尼，CT为接触阻尼，G为陀螺力矩；[K+KT](t)为刚度矩阵，其中K为系统刚度，KT为接触刚度；{u}(t)为系统自由度向量；{F+FT}(t)为力向量，其中F为系统外力和惯性力，FT为接触相互作用力；∆t为时间增量步长；β、γ为Newmark 法中加速度矢量和速度矢量的时间步参数占比控制常数。

显式动力学仿真计算所需的计算资源和总运算时间取决于每个增量时间步中方程解算时间以及总的计算步数。方程解算时间一方面取决于矩阵和力向量的数据量大小（即计算模型网格的总自由度数），另一方面取决于上述矩阵中有多少数据需要参与到碰摩、损伤及穿透等复杂的迭代计算过程，例如采用罚函数法进行局部碰摩行为的迭代计算。这2 种因素都与结构的具体设计特点和物理过程有关，需要根据仿真问题进行具体分析。

总的计算步数由仿真目标时长和时间步长∆t决定，尤其是∆t与网格的尺度有紧密联系，而且与计算过程的精度和稳定程度息息相关，也需要根据仿真问题进行具体分析。β、γ只能间接地在小范围内调节时域积分计算的速度和效率。

由上述分析可知，针对给定的仿真结构和物理过程，能够直接影响计算速度和结果精度的只有网格划分和∆t2 种因素，而且二者之间又有紧密联系。同时由上述计算方法可知，在显式计算中没有前后时间步的迭代收敛计算，因此前序时间步中形成的计算偏差会随时间步的不断增加逐步积累放大。

1.2 网格与时间步长的理论关系

在∆t的选择上，由于瞬态显式动力学分析所关注的物理过程包括碰撞（时长约10-2s 数量级）、冲击（时长约10-3s 数量级）和爆炸（时长约10-4s 数量级）[13]等，由于其时间历程很短，可能与计算网格中应力波传递时间尺度相近（如图1 所示），因此必须考虑应力波的影响。

图1 物理过程结构

根据显式动力学计算理论[14]，积分步长必须小于临界时间步长，否则会使1 个计算时间步长内应力波传播通过多个网格，使结点参数变化剧烈，从而使计算过程不稳定。临界时间步长取决于计算网格中最小尺度单元的特征长度

材料中纵波的传播速度为

式中：∆tcr为临界时间步长；Lch为网格特征长度；c为材料中纵波的传播速度，即材料中的声速；E为材料压缩弹性模量；μ为材料泊松比；ρ为材料密度。

该特征长度在不同的网格类型和不同的计算程序中有多种设置方法，以单个六面体实体单元为例，Lch通常有最短边长、最短对角线以及其他计算方法

式中：lx、ly、lz分别为x、y、z方向的单元边长。

以六面体单元最短边长法计算临界时间步长为例：如果材料为TC4 钛合金[15]，在100 ℃时，E=110 GPa，μ=0.34，ρ=4440 kg/m3，算得声速为6175 m/s，网格最短边长为10 mm，则其开展瞬态显式动力学分析的临界时间步长为1.62×10-6s。

上述研究的临界时间步长是从计算稳定性角度给出的时间步长的理论上限，小于这个理论步长计算过程也可能顺利完成，但是其计算结果可能不符合真实的物理过程，因此还需要根据实际物理过程进一步计算分析。

2 FBO的物理过程及有限元建模

2.1 FBO的物理过程

FBO事件主要有以下几个典型的瞬态物理过程：

（1）飞失的风扇叶片与机匣和后续叶片的碰撞损伤过程；

（2）飞失后风扇转子不平衡量显著增大，低压转子在不平衡离心力作用下变形破坏过程；

（3）结构载荷不断增大，触发结构保险破坏失效[4]过程；

（4）剩余低压转子稳定在风车转速持续运转过程[16]。

前3 个物理过程不是孤立的，而是综合体现在FBO 后的很短时间内（转子旋转几周，总时长在10-2～10-1s 数量级）不平衡转子与机匣的碰摩、变形和快速减速过程，主要约束条件是适航要求中叶片包容和大振动载荷外传等条款[1-3]；第4个物理过程实际上也包含剩余叶片与机匣碰摩的物理过程，但是由于持续时间较长，一般采用隐式动力学方法，主要约束条件是适航要求中转子持续旋转条款[1-3]。

由于第3个物理过程与结构保险的具体结构形式息息相关，而且在仿真应用中多采用带有门限值的接触条件方法来模拟，并且第4 个物理过程关注的时间范围和采用的分析方法不在本文讨论的范围内，所以本文中只研究前2个物理过程。

2.2 模型简化、网格划分及边界条件设置

根据新一代双转子大涵道比涡扇发动机结构设计特点[17]和具体结构参数[18-19]，通过一系列由简到繁的模型简化试算分析，逐渐积累对于具体结构设计特点及物理过程的认识，建立了显式动力学分析简化模型，进行网格划分和参数设置，如图2 所示并见表1。

图2 仿真模型及网格

表1 风扇叶片参数及仿真目标参数

2.2.1 模型简化

为了聚焦叶片与风扇机匣（包括易磨耗涂层）的碰摩过程，排除前序叶片对于后续叶片的影响，采用仅保留1 片叶片的风扇转子，但是叶片的材料、根/尖半径、不同截面的弦长和厚度以及叶片的扭转角等主要参数与真实发动机的基本一致；为了准确模拟低压转子的动力学行为，低压转子的结构形式、尺度、相互连接的拓扑关系以及材料与真实发动机的基本一致；为了准确模拟叶尖与机匣碰摩所产生的接触刚度和阻尼，以及易磨耗涂层的剥落消耗过程，风扇机匣和易磨耗涂层的结构和材料与真实发动机的基本一致；为了降低仿真结构的复杂性，低压转子的3 主轴承都采用弹簧单元来模拟。

2.2.2 网格划分

为了更加准确地控制网格的尺度，并更好地模拟碰摩损伤过程，全部模型采用六面体网格；在FBO 物理过程中，风扇叶尖与风扇机匣的碰摩是接触速度最高、破坏作用最明显、对于结构系统影响也最大的细节过程，因此选择风扇叶片和易磨耗涂层2 个部件（碰摩对）作为主要关注部件，其网格尺度根据研究需要采用不同的尺度；其余部件网格出于简化计算的需要，采用较大的尺度，而且各计算项目保持一致。

2.2.3 边界条件设置

低压转子设置转动边界条件模拟转子转动行为；风扇机匣和转子轴承的静子端采用固定边界条件来约束整个转静子模型系统；转/静子之间事先不限定是否接触，由软件根据计算结果自行判断是否接触，接触后以恒定的摩擦系数方法确定接触面的边界条件。

2.3 测量参数及测点布置

针对前2 个物理过程，选择3 个有代表性的输出参数，并根据模型的实际情况布置参数测点，如图3所示。

图3 叶尖碰摩的测试点位置

为了反映FBO 状态下，离心力与叶片碰摩产生的附加接触刚度和阻尼的共同影响，在风扇盘前缘轴心位置设置了轴心轨迹测点，输出形式为以时间为X轴的空间曲线，衡量标准为不同网格尺度和时间步长的计算项目得到的空间曲线的一致性；为了反映叶片碰摩过程中叶片的动力学响应和损伤情况，在叶片中部沿叶高方向等距均布了5 个应力响应测点，应力取值采用米塞斯等效应力，形式为相同测试位置的应力时序变化，衡量标准为不同网格尺度和时间步长的计算项目得到相同位置应力随时间变化的趋势和量值的一致性；为了反映转子通过轴承外传的振动载荷，在1 支点轴承水平位置设置了振动加速度测点，形式为同一测试位置的加速度时序变化，衡量标准为不同网格尺度和时间步长的计算项目得到加速度时序的变化趋势和量值的一致性。

3 基于DOE 方法的参数研究方案及总体完成情况

采用网格尺度和时间步长为输入参数，通过大范围的初步试算，缩小参数范围，最终确定每个参数的4个取值，通过排列组合形成16种不同计算参数的仿真计算项目表。在相同的硬件和软件平台下开展了上述计算项目表中所有的仿真计算，完成情况、总计算耗时及系统节点数量见表2。

表2 不同参数的仿真计算项目及计算情况

4 仿真结果分析

4.1 网格无关性检验

在仿真分析中，时间步长为3×10-7s 和2×10-7s2 组不同网格尺度的仿真计算全部完成，分组对比分析如下。

4.1.1 网格尺度对转子轴心轨迹的影响

以时间步长同为3×10-7s的4个仿真计算为例，如图4 所示。从图中可见，网格最粗的计算结果A1-B3 线在第0.02 s 左右就与其他3 条曲线出现明显差别，而其他3条曲线在第0.05～0.06 s和第0.08～0.10 s也出现了较大偏差，从仿真视频来看，在第0.05～0.06 s 和第0.08～0.10 s 都出现了叶片与机匣的剧烈碰摩。

图4 时间步长为3×10-7 s的3维轨迹

以时间步长同为3×10-7s 的4 个仿真计算为例开展同样分析的结论基本一致，只是第0.08 s 之后的虚线变化时序与上述分析略有差别。

4.1.2 网格尺度对叶片应力的影响

以时间步长为2×10-7s 的4 个仿真计算结果为例，在不同的网格尺度下，叶身正中部测点的应力对比如图5 所示。从图中可见，网格尺度为40 mm 的计算结果曲线与其他3 条明显不一致，随着网格逐渐加密，叶身应力的差异逐渐缩小，尽管其总的趋势基本一致，但是在第0.08 s 之后具体的应力大小变化的时间顺序规律不一致。