两机叶片丨复旦大学：基于GRU神经网络的航空发动机叶片颤振快速预测方法

颤振是由非定常空气动力与弹性结构体相互作用产生的自激振动失稳现象。随着现代航空技术的发展，航空发动机不断追求高增压比使得结构的气动负荷加重，同时追求低耗油率使得涡扇发动机涵道比增大，因此风扇直径及其叶尖速度也随之增大，进而导致颤振问题变得更加突出和复杂。颤振不仅可能导致航空发动机性能的急剧下降，还可能引发结构的疲劳破坏，甚至造成灾难性的飞行事故。因此，准确预测叶片的颤振，可以有效降低设计和测试中的颤振风险，保证航空发动机的可靠性。

当前，颤振预测研究已取得许多进展。计算流体动力学(Computational Fluid Dynamics,CFD)已成为研究发动机部件颤振机理的常用手段，计算结构动力学(Computational Structural Dynamics,CSD)的发展也为颤振分析提供了强大的数值分析工具。CFD方法可以计算不同飞行条件下的气动力特性，CSD则可以分析结构在不同气动载荷作用下的响应。多位学者对不同耦合形式的CFD/CSD计算方法进行了研究，但若想获得较高的计算精度，仍然需要付出很大的计算代价。

影响系数法是航空发动机叶片颤振分析中一种常用的气动力建模方法，该方法结合气动弹性控制方程可以构建特征值问题来预测结构的颤振特性。Bendiksen等就采用影响系数法分析了叶栅的弯曲与扭转振型的气动力，然后采用结构气动弹性特征值法分析了弯曲与扭转耦合的颤振边界。Clark扩展了Bendiksen等的特征值方法，采用条形理论将多个叶片的展向截面纳入分析，并研究了开式转子的频率合并颤振问题，他们发现质量比、频率分离对颤振的临界转速有着显著影响。质量比更高的叶片需要更大的旋转速度才会经历耦合模态颤振,并且增加模态之间的频率分离也会增加临界转速。Ananth等采用Moore-Greitzer模型(MG模型)对经典的激波流动进行了建模，采用线性化Whitehead理论确定叶栅的气动载荷，最后通过将二维结构控制方程转化为特征值问题来检验叶栅的气动弹性稳定性。他们的研究提供了一种准稳态预测方法，可用于分析经典激波流动条件下的耦合弯曲-扭转颤振。此外，Sun等还提出了一种扩展影响系数法，并结合刚体单元转子模型，得到了考虑柔性支承的齿轮副的动态啮合刚度。

从能量角度也可对结构的颤振特性进行评估。基于能量平衡原理，通过比较结构振动过程中的能量输入和耗散，可以判断颤振是否发生。能量法的物理概念清晰，当能量从非定常流转移到叶片时，即一个周期内非定常流对叶片的净功为正，这意味着将发生颤振，反之，非定常流对叶片的净功为负时，意味着叶片稳定。能量法通常可以提供更快的计算速度，但是需要注意气动压力相位相对于叶片运动相位的关系，这对于准确预测叶片表面的非定常气动压力非常重要。能量法已被诸多学者广泛应用到发动机叶片的颤振分析中，其能快速处理单模态和多模态颤振，及弱阻尼和强阻尼等复杂情况，适用于工程应用中降低计算成本，进行快速预测和稳定性评估。

模态振型是影响颤振稳定性的重要因素，不仅叶片表面非定常压力分布会其影响，而且模态振型的改变会导致最小气动阻尼相对于节径的变化。模态分解通过将结构的复杂振动行为分解为若干个简单的模态，从而简化颤振分析的过程，是结构动力学中的一种基本分析方法。Vogt等研究了典型低压涡轮转子叶片的模态振型灵敏度，对一组三个正交模态进行了测试，发现最稳定的模态是具有轴向到弦向特征的弯曲模态，而扭转模态主导的模态灵敏度较高。此外，可以将叶片的第一阶弯曲模态(1F)分解为弯曲分量和扭转分量，并以两种分量模态为基础进行颤振分析。Isomura通过对跨声速风机近颤振工况的数值模拟，首次研究了风机叶片模态扭转分量对颤振的影响，发现1F模态扭转分量在与弯曲分量同相时会增加叶片阻尼，反相时则会增大叶片的激励。模态振型分解能够在复杂流动下揭示结构的动态特性，识别出关键模态分量，为颤振预测和控制提供依据。

神经网络方法利用数值仿真获得的数据，能够建立在特定条件下预测广义气动力的数学模型,如连续时间递归神经网络(continuous-time recurrent neural network,CTRNN)、卷积神经网络(convolutional neural network,CNN）、循环神经网络(recurrent neural network,RNN)等方法。门递归单元神经网络模型（gated recurrent unit，GRU）是一种有效的循环神经网络变体，它在处理时间序列数据方面表现更为出色。RNN在长序列数据上训练时，容易出现梯度消失与梯度爆炸的问题，导致训练困难。GRU则通过引入更新门和重置门来控制信息的流动，从而解决了这一问题，并且与长短记忆网络（long short-term memory，LSTM）模型相比，GRU模型减少了门结构的数量，直接通过隐藏状态传递信息，使得计算效率大大提升。在时间序列预测的场景下，GRU的表现也优于LSTM。GRU还可以与其他模型组成混合模型，进一步提高预测精度。对于本文气动力模型，希望通过单个叶片振动的广义位移时间历史来推断作用在各叶片上的瞬时广义气动力。该模型的输入是单个叶片的广义位移时间序列，因此在建模过程中无需通过类似CNN的结构对输入数据的空间特征进行提取和表征。另外，由于叶片瞬时气动力主要依赖于临近时间窗口内的叶片振动位移历史，因此我们也没有在模型中引入类似大语言模型中的长文本注意力机制，而只选择了单纯的时间序列神经网络模型。考虑到GRU网络相对LSTM或RNN网络的计算效率优势，最终选择GRU网络作为本文的气动力模型。

本研究提出了一种基于模态振型分解并结合GRU神经网络模型计算气动功的颤振快速预测方法。首先，通过模态振型分解，将叶片的模态振型分解为纯弯曲和纯扭转两个基本模态，计算在这两种基本模态组合下的气动力，用于GRU模型进行时间序列预测，以得到不同弯扭比、不同叶间相位角和不同振动频率下的中心叶片表面气动力。利用该气动力模型，可以对具有不同弯扭比的特征模态的叶片进行快速颤振预测，从而加快发动机叶片的颤振设计。

1.1 气动阻尼的计算

对于给定的发动机叶片，模态力可以定义为：

式中：F_m,i是0号中心叶片的振动对i号叶片施加的模态力；i代表叶片编号，本文共计算了7个叶片，取值范围为-3～3，其中i=0代表中心叶片（如图1所示）；p为叶片表面压力；Φ为模态振型矢量；ds为叶片表面面元；n为叶片表面单位法向量。

图 1 叶片编号

通过离散傅里叶变换计算模态力的一次谐波：

由此，将本算例中0号中心叶片上的气动阻尼定义为：

这里不考虑叶片的结构阻尼，因此，当气动阻尼ζ(σ)大于零时叶片不发生颤振，当气动阻尼ζ(σ)小于零时叶片发生颤振。

1.2 模态振型分解

本研究利用商用有限元分析软件对NASA Rotor67转子叶片模型进行了模态分析。叶片材料假定为钛合金TC4，假设转子盘具有较强的刚度，1F模态的模态振型和固有频率不随节径变化。转子根部采用固支边界条件，考虑了离心力的预应力影响，有限元网格和计算的1F模态振型如图2所示。计算得到的1F模态振型的固有频率为533Hz。

之后，将1F模态振型分解为两个正交模态分量。首先，根据叶片表面上的有限元节点排布在不同展向截面位置提取弯曲模态分量。在展向的第i个截面上，通过对前后缘位移矢量Φ_1F,i取算术平均值来获得弯曲模态分量Φ_plunge,i，如式(7)所示：

得到弯曲模态分量后，扭转模态分量Φ_twist,i定义为：

图2 结构分析网格和模态振型

分解后的弯曲模态Φ_plunge和扭转模态Φ_twist如图3所示。

图3 1F模态分解出的模态振型

1.3 GRU神经网络

GRU是一种时序神经网络，具有仅使用少量参数就可以捕捉复杂时序数据特征的能力。Hinton等通过研究发现，只要增加神经网络模型中的层数，将上一层的学习结果当作下一层的训练数据，就可以有效提高神经网络的学习能力，因此构建GRU模型时，通常会选择至少两层隐藏层，而隐藏层维度则需要根据数据特征选择适当的维度。在本研究中，预测模型由两个相同的GRU层和一个全连接层组成，其结构如图4所示。本文通过预实验测试了不同的参数组合，最终发现设置GRU层的隐藏层特征维度为32、循环层数为2可以有效避免过拟合，同时能够确保模型有效捕捉到输入数据中的时序依赖关系。然后，通过一个全连接层将隐藏层的输出转化为所需的14维输出维度。在模型训练中使用的批处理大小为32，学习率为0.0005，这两个参数的设置是经过一系列预实验和参数调整后的结果，目的是为了保证模型在训练过程中能够稳定收敛，同时保持较好的计算效率。

图4 GRU网络结构图

本文中的数值模拟计算域由如图1所示的7个叶片通道组成，中心叶片被施加小振幅振动，其他叶片静止。采用时域数值模拟方法计算中间叶片振动引起的非定常流场，从而获得各叶片上的非定常气动力的时间序列。考虑到中心叶片非定常压力扰动沿周向衰减很快，因此在计算中只考虑了中心叶片和其两侧相邻的各3个叶片，共7个叶片，叶片编号为-3,- 2,⋯,+2,+3 。

在GRU模型的训练中，为了生成覆盖结构固有频率的某一频率范围内的非定常气动力时间序列样本，在时域数值模拟中，中心叶片的强迫振动形式设定为变频率运动。中心叶片的运动时间函数定义为:

2.1 定常计算

采用商用软件CFX对叶片周围流场进行数值模拟计算。NASA Rotor67有22个叶片，是两级跨声速压气机的第一个转子，设计点流量为33.25kg/s，压力比为1.63，转速为16043r/min。为了获得叶片在静止状态下的平均流场，首先在设计转速下对NASA Rotor67进行单通道的定常流场计算。计算的进口边界给定总压、总温度和气流角，出口边界采用径向平衡方程，湍流模型采用标准壁函数的k−ε模型，计算网格节点数为46万，如图5所示。

图5 NASA Rotor67的计算网格

为了验证CFX数值仿真计算的可靠性，本文通过改变背压对单通道叶片进行不同工况下的定常计算，并将计算结果与Strazisar等的实验结果进行对比验证。图6展示了通道压比和效率曲线的对比，图中红点为Strazisar等的实验结果，蓝色曲线为本文计算结果，可以看出计算结果与实验数据是相符的。

2.2 非定常计算

在非定常计算中，本文先由式(9)生成一个从480Hz变化到580Hz的变频运动函数,最大位移为0.001m,采样间隔为3.1222×10⁻⁶s，采样点个数为13000个。弯曲模态位移和扭转模态位移之间的相位差设定为ϕ=30°，两种模态位移的时间函数曲线如图7所示。

在非定常计算中，中心叶片以式(10)的形式振动：

其中，P(t)与T(t)分别表示弯曲模态和扭转模态位移时间函数，并且定义弯扭比Q=a:b,a表示弯曲模态权重系数，b表示扭转模态权重系数。对于叶片的固有1F模态，Q=1:1。计算得到的7个叶片上的弯曲模态力和扭转模态力时间序列如图8所示，从图中可以看出，弯曲模态力相比扭转模态力数值更大，且同一叶片上具有更大的振动幅度。

图6 通道压比和效率的计算结果与实验结果对比

图7 弯曲模态和扭转模态位移的时间函数曲线

图8 模态力数据快照

2.3 GRU模型的训练

GRU模型的输入量为振动叶片弯模态与扭模态的广义位移时序函数，输出为周围7个叶片弯曲模态和扭转模态的广义气动力。本文中GRU模型的训练流程如图9所示。为了对GRU模型进行训练，将在非定常数值模拟中获得的叶片变频振动下的气动载荷时间序列作为训练数据放入GRU模型中进行训练，同时通过非定常数值模拟生成另一个测试集，该测试集由中心叶片以1F模态固有频率533Hz振动时计算获得。训练的损失函数曲线如图10所示。可以看到，在训练开始以后，训练集和测试集的损失函数都迅速下降，但是在50个迭代步之后，训练集和测试集误差都基本不再减小，故训练停止。训练后的GRU模型对于训练集中变频振动工况下各叶片弯曲和扭转模态力的预测和计算数据的对比如图11所示。

图9 GRU模型训练流程

图10 GRU模型的损失函数曲线

图11 变频激励模态力训练数据与预测数据对比

用广义力时间序列的均方根误差(root-mean-square error,RMSE)来衡量GRU模型预测的准确度，其计算公式如式(11)所示：

其中，M表示一个周期内的样本数，y_k与y_pk分别表示模态力仿真值和模态力预测值。e_RMSE反映了广义力仿真值和模型预测值之间的逐点偏离程度，因此其中同时包含了时间序列的幅值和相位误差信息。同时，由于不同叶片上不同模态的广义力大小差别很大,本文计算了一个周期内每个叶片的e_RMSE与对应模态力峰峰值(S_pp)之间的比值，作为模态力预测的相对误差，其结果如表1所示。

表1 训练集预测值误差分析

通过表1和图11可以看到，GRU模型对于0号振动叶片及附近叶片上的弯曲和扭转气动力的预测都十分准确，相对误差值均较小。而对于较远处的叶片，由于其周围流场受到0号振动叶片的扰动较小，非定常气动载荷的幅度也相对更小，因此这些叶片上的气动载荷的GRU模型预测结果相对误差较大。

对于作为测试集的533Hz叶片简谐振动工况，其GRU模型的预测结果和数值模拟结果的对比如图12所示，相应的相对误差计算如表2所示。

图12 533Hz激励模态力仿真数据与预测数据对比

表2 测试集预测值误差分析

从图12和表2可以看出，GRU模型对于测试集简谐振动工况下各叶片上的弯曲和扭转气动力的预测整体上也比较准确。尽管模型对于较远处叶片（+3叶片或-3叶片）上非定常气动力的预测相对误差还是略大，但是由于较远处叶片上的气动载荷的绝对值很小，因此其相对误差对于叶片颤振分析中气动阻尼计算的影响相对较小，这一点本文将在关于气动阻尼的计算部分进行详细考察。

3.1 变叶间相位角对气动阻尼的影响

在完成了GRU模型的训练之后,本文将模型应用于叶片的快速颤振分析。首先考查叶片在固有1F模态、固有频率533Hz下的颤振稳定性。设定Q=1:1不变,从−180°到180°遍历不同的叶间相位角(IBPA)，根据式(1-6)描述的过程计算气动阻尼随IBPA的变化情况，结果如图13所示。

图13 Q=1:3、固有频率为533Hz时不同IBPA下的

气动阻尼变化趋势

可以看出，中心叶片的气动阻尼在叶片几乎同相位振动（IBPA为0°附近）时最小，最容易发生颤振，而在叶片反相位（IBPA为180°附近）时最大，最不容易发生颤振。这与实际物理情况相符，因为在叶片同相位振动时，叶片对流场的扰动会相互叠加放大，更容易导致流固耦合系统的失稳，而反相位振动时叶片对流场的扰动会相互抵消，从而减小系统的扰动放大，抑制颤振。另外，还可以观察到，对于NASA Rotor67叶片，在固有频率533Hz下，1F模态在各叶间相位角情况下均不会发生颤振，这也与NASA Rotor67叶片的实际情况相符。

3.2 变频率对气动阻尼的影响

在发动机叶片的颤振设计中，考虑到发动机内流环境包含的复杂流动情况，往往需要考查叶片在一定频率范围内的颤振特性。因此，可以首先在给定频率下，遍历IBPA计算叶片气动阻尼，然后取出各IBPA下最小的气动阻尼（因为该点是颤振风险最大的点）作为该频率下获得的叶片气动阻尼。这里设定Q=1:1（即叶片固有1F模态），将激励频率以10Hz为间隔，从480Hz步进到580Hz，在每个频率下计算各叶间相位角下的最小气动阻尼，计算结果如图14所示。

图14 Q=1:3时各叶间相位角下的最小气动阻尼

随激励频率变化趋势

可以看到,对于NASA Rotor67叶片的1F固有模态，随着激励频率的增大，叶片的最小气动阻尼逐渐增大，呈现出气动弹性趋于稳定的趋势。

3.3 变弯扭比的气动阻尼计算验证

为了验证本文得到的GRU模型对于不同弯扭比下叶片的气动阻尼预测的准确性，选取Q=1:3、受迫振动激励频率为520Hz的工况进行分析。首先进行非定常广义气动力的数值仿真计算，同时应用GRU模型对相同工况下的气动力时间序列进行推断，然后应用式(1~6)对不同IBPA下的气动阻尼进行计算，对比结果如图15所示。

图15 Q=1:3、激励频率为520Hz时不同IBPA下的气动阻尼对比

由图15可以看出，GRU模型预测的气动力计算得到的气动阻尼与数值模拟获得的结果相比误差很小，由此验证了本文提出的GRU模型可以用于叶片气动阻尼的快速推断，且能较准确地评估叶片是否发生颤振。

为了在更大工况范围内考察GRU模型在气动阻尼计算中的准确性，本文在520~580Hz范围内改变叶片简谐振动频率，分别计算了Q=1:3和Q=3:1时不同IBPA下的叶片气动阻尼。同时，为了对气动阻尼的计算误差进行量化分析，本文计算了各工况下不同IBPA下的气动阻尼均方误差(DRMSE)，并且计算了相对于各工况下气动阻尼变化范围（气动阻尼最大值与最小值之差，△D）的相对误差，结果如表3所示。

表3 气动阻尼的误差分析

从表3可以看出,在叶片固有频率(533Hz)附近，GRU模型均能对叶片气动阻尼进行较准确的预测，但随着简谐振动频率的增加，计算得到的气动阻尼误差不断增大。考虑到在本文的模型训练中，样本气动载荷是由在480~580Hz范围内变化的变频叶片振动激励的，因此若要在高频率下获得更准确的气动阻尼估算，需要进一步增大训练样本的扫频范围。

3.4 变弯扭比对气动阻尼的影响

应用本文提出的方法分析叶片弯扭比对气动阻尼的影响，将弯扭比分别取为Q=1:3和Q=3:1，计算500~560Hz激励频率范围内叶片在不同IBPA下的最小气动阻尼。两种弯扭比下的气动阻尼对比如图16所示。

从图16可以看出，与Q=3:1的叶片相比，Q=1:3时中心叶片的最小气动阻尼值始终更小，即叶片的特征模态中扭转分量的占比越大，叶片越容易发生颤振，这与Cumpsty等在研究中发现的弯曲和扭转模态对叶片颤振的影响一致。

图16 不同Q下最小气动阻尼随激励频率的变化对比

本文基于GRU神经网络模型及模态振型分解方法提出了一种有效的快速颤振预测方法，并以NASA Rotor67模型为例对方法进行了验证。具体结论可归纳为以下几点：

1) 在GRU神经网络模型的训练中，通过引入覆盖一定频率范围的变频强迫振动生成的训练样本，可以使模型在这个频率范围内都适用。

2) 通过在生成训练样本时引入弯扭模态之间的初始相位差，可以使训练样本同时覆盖不同弯扭比的情况。在本文的验证算例中,GRU模型在弯扭比Q=1:3和Q=3:1时，均可以较准确地对叶片的气动阻尼进行快速计算。

3) 通过本文提出的模型，可以对弯扭模态对叶片颤振的影响进行快速有效的评估，避免了大量的非定常流场数值计算，这大大提高颤振分析和颤振设计的效率。

在未来工作中，本文将考虑将更多颤振影响因素纳入目前的模型中，从而在更广的参数空间中构建颤振分析快速计算模型。同时，也考虑将模型应用于不同几何构型、不同工况下的发动机叶片颤振分析。此外，鉴于颤振问题的相似性，该模型也可以进一步推广到机翼颤振等外流颤振问题中。