全转速系数矩阵降维重构的燃机不平衡量逆推方法

2020-12-01王辰左彦飞江志农胡明辉冯坤

航空学报 2020年11期

王辰，左彦飞,*，江志农，胡明辉，冯坤

1. 北京化工大学发动机健康监控及网络化教育部重点实验室，北京 100029 2. 中国航发动力所-北京化工大学航空发动机振动健康监控联合实验室，北京 100029

转子是燃气涡轮发动机的核心部件，其振动问题不可避免[1]。转子不平衡是导致燃气涡轮发动机整机振动过大的最主要原因之一，直接影响发动机的可靠性及使用寿命[2]。传统动平衡技术虽较为成熟，但需拆下转子分解后，各部件单独平衡再逐级装配并装入机匣，这样动平衡效果易受装配过程中不确定因素的影响，且效率低下，往往拆装多次进行动平衡，振动也无法达标；而整机动平衡无需拆装转子，可直接在试车台上完成，在保证动平衡效果的同时提高效率、降低成本，是燃气涡轮发动机动平衡技术的一个重要发展方向[3]。同时，燃气涡轮发动机由于其自身结构及工作原理等原因，只能在机匣外表面安装接触式振动传感器进行振动测试，与化工、能源等领域广泛使用的旋转机械测振手段不同，为整机动平衡带来了更为复杂的理论问题和实践困难。

整机动平衡的关键是通过机匣测点振动计算转子不平衡量。现有转子不平衡量的计算方法很多，如基于最小二乘影响系数法改进的不平衡量计算方法[4]、基于谐波平衡-频时转换的转子不平衡量计算方法[5]等。此外，不平衡振动相位测算是确定不平衡量的重要步骤，常用方法有基于互功率谱的不平衡振动相位检测方法[6]、基于最小二乘法基准信号的振动相位计算方法[7]等。上述研究为转子不平衡量的计算识别奠定了基础，不过多是通过直接测量转子本身振动响应而不是基于机匣外表面测振数据进行计算的。

对燃气涡轮发动机而言，振动测点大多只能布置在机匣外表面，无法直接测量转子自身的振动响应[8]，转子的不平衡振动要经过轴承、挤压油膜阻尼器和结构框架才能传递到机匣外表面，导致机匣表面测点在不同转速下的不平衡振动幅值、相位变化规律十分复杂，增大了转子不平衡量的计算难度[9]，如图1所示。

因此可见，深入分析研究整机的动力特性是必要的。对于转子支承系统，邓旺群等[10]建立了涡轴发动机转子有限元模型，计算了转子的临界转速，对转子的复杂振型进行了研究；Wagner和Helfrich[11]对燃气涡轮发动机转子在预应力作用下的动力学特性展开了相关研究，获得了转子复杂的动力学行为；Kuaia等[12]对燃气涡轮发动机转子在不同热条件下的动力特性进行了研究，掌握了温度对转子动力学行为复杂的影响规律。对于机匣系统，温登哲和陈予恕[9]指出了燃气涡轮发动机机匣动力学特性的复杂性，阐述了机匣对整机动力学特性的影响。马英群等[13]研究了多转子激励下的机匣振动响应特性，得出振动能量在机匣上的复杂传递特性。对于整机系统，Salles等[14]考虑了陀螺效应以及转子与机匣之间的摩擦，对燃气涡轮发动机的动力学特性进行了研究，为发动机结构优化设计提供了帮助；张大义等[15]采用整机三维模型对双转子涡扇发动机的固有振动特性进行了计算和评估，考察了转子系统总应变与共振裕度是否满足转子动力学设计要求；刘永泉等[16]针对高性能燃气涡轮发动机结构复杂性和高温高转速工况下动力学稳定性问题，明确了转子动力学特性设计分析与整机振动控制技术的关系；孙凯等[17]针对发动机工程研制过程中整机结构分析在模型规模、模型质量及分析效率等方面面临的一系列挑战提出模型简化、混合分网、模型修正和连接结构的建模方法；Zuo和Wang[18]对燃气涡轮发动机的3D有限元模型减缩方法进行了研究，提高了仿真分析的精度与效率。

图1 某双转子发动机振动传递及测点示意图Fig.1 Vibration transmission and measure points of a double-rotor engine

总之，目前对燃气涡轮发动机的动力学特性研究较全面，已经普遍认识到柔性转子、机匣及整机系统结构和动力特性十分复杂，且转子-机匣系统在特定转速下存在不同程度的相互耦合。这为燃气涡轮发动机的整机动平衡带来了更大困难，即单一转速下进行的动平衡无法保证全工作转速区间内的有效性，且由于振动的复杂性，振动测点及配平位置的优选组合难以确定，而目前对于上述困难应对方法的研究较少。虽然冯坤等[19]考虑了机匣的复杂动力特性，对轴承位置的不平衡响应矢量进行了逆推方法研究，但未考虑转子在工作转速区间内的复杂动力特性以及振动测点、配平叶盘的选取问题，且整机动平衡的成熟实施案例也颇为少见，仅在少数民用燃气涡轮发动机上有所应用[20]。

为此，针对如何由机匣表面测点振动响应逆推转子不平衡量的问题，通过机匣测点振动响应与转子不平衡量之间的对应关系，将转子各叶盘上分布的不平衡量集中等效到拟进行配平的叶盘上，提出一种基于机匣振动响应的转子等效不平衡量逆推方法，并以典型双转子发动机高保真模型进行数值仿真分析，验证该方法的有效性。

1 不平衡振动传递机理

首先，将燃气涡轮发动机的“转子-轴承-机匣”系统简化为转子-支承系统，则系统在不平衡激励下稳态响应方程为

(1)

式中：M为转子质量矩阵；z为转子的振动位移向量；C为系统阻尼矩阵；Ω为转子的转速；G为转子陀螺矩阵；K为转子刚度矩阵；q为转子的不平衡量向量，包含不平衡量的大小和相位；t为时间变量。

式(1)的解可表示为

z=A(Ω,q)eiΩ t+φ=Ω2qeiΩ t+φ[-MΩ2+

Ω(C+ΩG)+K]-1

(2)

式中：A为Ω和q的函数；φ为转子振动位移相位。

由于转子支承点即轴承的振动位移向量zb为转子振动位移向量z的子集，故zb与z中相关元素存在确定的函数关系f1；同理，轴承振动位移相位φb与φ也存在确定的函数关系f2，即

zb=A(Ω,q)eiΩ t+φbf1(M,C,K,G)=

A(Ω,q)f1eiΩ t+φb

(3)

φb=φf2(M,C,K,G)=φf2

(4)

然后，将转子-支承系统中的支承恢复为“轴承-机匣系统”，将转子支承点位移zb作为激励施加在轴承位置，求取机匣测点的振动响应。对于机匣，其动刚度Kcasing为

Kcasing=Kc-McΩ2+iCcΩ

(5)

式中：Kc为机匣刚度矩阵；Mc为机匣质量矩阵；Cc为机匣阻尼矩阵。

同理，对于轴承，忽略其质量，则动刚度Kbearing为

Kbearing=Kb-MbΩ2+iCbΩ=Kb+iCbΩ

(6)

式中：Kb为轴承刚度矩阵；Mb为轴承质量矩阵；Cb为轴承阻尼矩阵。

“机匣-轴承”系统可等效为机匣系统与轴承系统的串联，故系统的等效总刚度Kbc为

(7)

当在轴承上施加振动位移激励zb时，施加到机匣的激励力p可表示为

p=Kbczb

(8)

故可求得机匣的振动位移向量zc为

(9)

则机匣振动位移幅值Bc为

Bc=|A(Ω,q)f1|×

W(Ω)|q|

(10)

式中：W(Ω)为关于Ω的函数。

而机匣振动相位φc与轴承振动位移相位φb的相位差为

φc-φb=arctan[CbΩ(Kc-McΩ2)-KbCcΩ]/

[Kb(Kc-McΩ2+Kb)+CbΩ(CcΩ+CbΩ)]=

arctanV(Ω)

(11)

式中：V(Ω)为关于Ω的函数。

故机匣测点振动位移幅值及相位可表示为

(12)

由式(12)可知，在忽略非线性及不确定因素并将系统简化处理的情况下，机匣测点振动幅值及相位与不平衡激励之间的关系已较为复杂。不过，机匣测点振动幅值Bc、机匣测点振动相位φc均与转子转速Ω、不平衡量q直接相关。在特定转速Ω下，机匣测点的振动响应幅值与转子的不平衡量q线性相关；且当不平衡量相位确定时，转子的振动响应相位φ为确定值，则机匣测点的振动响应相位φc也为确定值。

因此可通过动力学仿真分析或实际测试，获得燃气涡轮发动机在不同稳定转速下，机匣测点的振动响应与转子不平衡量之间的函数关系，为基于机匣实测振动响应逆向求解不平衡量奠定基础。

2 等效不平衡量逆推方法

由于燃气涡轮发动机整机具有复杂的动力特性，机匣外表面测点振动幅值及相位随不平衡量分布、转速变化的规律也十分复杂，因此需要考虑转子不平衡位置、机匣外表面测点位置、发动机工作转速区间等因素对等效不平衡量逆推的影响，建立针对整机动平衡更为有效的不平衡量逆推方法。

2.1 单一转速等效不平衡量逆推方程组

设某型燃气涡轮发动机存在m个叶盘，机匣外表面存在n个振动测点。在任意j叶盘位置不平衡量qj作用下，系统的振动位移响应频率与不平衡激励频率(即转子转频)相同，则k测点在转速ω下的振动位移zk(ω)可表示如下：

zk(ω)=Bkeiφk

(13)

式中：k为正整数，且k≤n；Bk为k测点的振动位移幅值；φk为k测点的振动位移相位角。定义测点k振动位移与叶盘j不平衡量之间的不平衡响应系数Skj(ω)：

(14)

Skj(ω)包含两部分信息，分别是振动响应幅值对不平衡量大小的放大系数以及振动响应与不平衡量之间的相位差，故Skj(ω)描述的是测点k对j号叶盘位置不平衡量的敏感程度。通过有限元数值仿真或实际测试可获得燃气涡轮发动机特定转速下不平衡响应系数矩阵S(ω)：

(15)

矩阵S(ω)的行向量代表某测点对各叶盘不平衡量的响应系数，列向量则代表各测点对某叶盘不平衡量的响应系数。

在特定转速ω下，依据有限元数值仿真或实际测试结果得到机匣各测点的振动响应数据，求取各测点的振动工频幅值、相位组成的响应向量：

z(ω)=[z1(ω)z2(ω)…zk(ω)…zn(ω)]T

(16)

定义各级盘的等效不平衡量向量q为

q=[q1q2…qk…qm]T

(17)

将不平衡响应系数矩阵、等效不平衡量矩阵和响应矩阵构造成不平衡量逆推方程组，表达该转速下等效不平衡量与响应之间的关系：

S(ω)q=z(ω)

(18)

若忽略机匣连接结构的微弱非线性，不考虑轴承间隙、挤压油膜阻尼器等存在的非线性因素，式(18)为线性方程组。求解式(18)可得各盘的等效不平衡量：

q=S-1(ω)z(ω)

(19)

由于燃气涡轮发动机工作转速范围广，且在不同的工作转速下，在某一工作转速下进行的动平衡可能无法满足全转速范围内的动平衡要求。因此，将不平衡量逆推方程组从特定转速拓展到全转速范围，构造全转速等效不平衡量逆推方程组。

2.2 全转速等效不平衡量逆推方程组

采用h个稳定转速对全转速范围进行离散，求取各稳定转速下的不平衡响应系数矩阵。

将ω离散为ω1,ω2, …,ωh，求取各稳定转速下的不平衡响应系数矩阵S(ω)，可得全转速区间不平衡响应系数矩阵SS：

SS=[S(ω1)S(ω2) …S(ωh)]T

(20)

各测点在全转速区间的振动响应向量zz为

(21)

构造全转速等效不平衡量逆推方程组：

SS(n×h)×mqm×1=zz(n×h)×1

(22)

对等效不平衡量q进行求解，即

(23)

但需要注意的是，全转速区间不平衡响应系数矩阵SS具有n×h行、m列，即方程组由n×h个方程和m个未知数组成，这就意味着方程个数远大于未知数个数。

在数学求解上，为减少线性相关度很高的冗余方程个数，降低方程组病态程度，从而提高求解稳定性，需要对全转速等效不平衡量逆推方程组进行简化，筛选出稳定求解必要的有效方程。而对于燃机整机系统，数学上的简化与配平叶盘位置、振动测点及转速筛选本质相通。因此，选择合适的配平叶盘、振动测点及对应转速可大幅提高转子不平衡量逆推的准确性，提高整机动平衡成功的概率。

2.3 全转速等效不平衡量逆推方程组的降维重构

为提高全转速等效不平衡量逆推方程组的求解稳定性，并考虑到工程应用可行性，需要同时缩减逆推方程组的方程个数和未知数个数，利用全转速系数矩阵的物理意义，结合数学方法，从配平位置选取、振动测点选取及全转速优化3个方面，对全转速等效不平衡量逆推方程组进行降维及重新构造，得到降维重构的全转速等效不平衡量逆推方程组(下称降维重构逆推方程组)，步骤如图2 所示。

图2 逆推方程组降维及重构流程Fig.2 Flow chart of dimension reduction and re- construction of backstepping equations

2.3.1 筛选配平叶盘

利用高斯消元(Gauss-Jordan)法对全转速区间不平衡响应系数矩阵SS进行初等行变换，并调整主元所在列最大元素数值相对0的容差，以控制所求极大无关列向量的数量，即根据转子实际振动情况筛选出在全转速区间内能够描述转子主要振动形态的叶盘数量及位置，作为全转速区间内动平衡的配平叶盘。

2.3.2 筛选转速及振动测点

将全转速区间不平衡响应系数矩阵SS转置，并挑选出与配平叶盘对应的行向量，大幅缩减不平衡响应系数矩阵的规模，形成缩减转置不平衡响应系数矩阵SS_mp；利用高斯消元法对矩阵SS_mp进行初等行变换，并调整主元所在列最大元素数值相对0的容差以控制所求极大无关列向量的数量，即利用机匣测点在全转速区间内与被筛选叶盘之间的不平衡振动响应关系，筛选出与各配平叶盘对应的振动测点及转速作为全转速区间内动平衡的数据来源。

2.3.3 构造降维后的等效不平衡量逆推方程组

构造降维后的逆推方程组系数矩阵Sdr，建立等效不平衡量列向量qdr与振动数据列向量zdr，初步完成对降维后逆推方程组Sdrqdr=zdr的构建，形成降维重构逆推方程组，展开后的形式为

(24)

式中：下标o,p, …,q表示测点号，i,j,…,k表示叶盘号；qi为叶盘i的不平衡量；zo为测点o的振动数据。

同时要求方程组式(24)有解且必须存在非零解，即通过调整振动测点的筛选，使降维重构逆推方程组系数矩阵的秩与其增广矩阵的秩相等。

2.3.4 优化降维重构逆推方程组求解稳定性

通过调整配平叶盘的筛选确定出条件数最小的降维重构逆推方程组，即最有利于方程组求解稳定性的配平叶盘组合，以减弱由实际振动数据测量及采集误差等因素对逆推结果引起的扰动。

根据工程应用情况，考虑配平叶盘个数后，将筛选出的不平衡响应系数和实测振动数据代入降维重构逆推方程组进行求解，即得到将转子不平衡量等效分配至各配平叶盘上的不平衡量。

3 数值仿真验证

3.1 发动机有限元模型及边界条件

建立某型燃气涡轮发动机高保真有限元模型，在各安装节处施加约束，构成与实际情况类似的边界条件。在不显著影响动力学特性的前提下，简化高、低压转子与机匣之间的支承关系以及内外涵道的连接关系。利用文献[18]所提方法生成超单元减缩模型，以大幅提高计算效率，如图3和图4所示。

3.2 基于仿真分析的全转速区间不平衡响应系数矩阵SS计算

以在1级叶盘位置施加幅值为10 g·m、相位为0°、频率区间为6～175 Hz(对应的转速为300～10 500 r/min)、频率间隔为1 Hz的扫频不平衡激励为例，进行谐响应分析计算。振动测点分别选取在风扇机匣(截面1)水平方向(前水)与垂直方向(前垂)、中介机匣前法兰(截面2)水平方向(中前水)与垂直方向(中前垂)、中介机匣后法兰(截面3)水平方向(中后水)与垂直方向(中后垂)、涡轮后机匣(截面4)水平方向(后水)与垂

图3 某型燃气涡轮发动机高保真有限元模型Fig.3 High fidelity finite element model of gas turbine engine

图4 某型燃气涡轮发动机整机减缩模型Fig.4 Supercell model of gas turbine engine

直方向(后垂)、燃烧室后机匣(截面5)水平(燃水)，总共9个测点，依次命名为测点1～测点9，如图5 所示。

分析计算其在各级叶盘处承受单位不平衡激励作用时，系统在工作转速范围内(300～10 500 r/min)的振动响应特征。依次提取机匣振动测点的位移响应(幅值与相位)，绘制曲线如图6所示。

图6中曲线为1级叶盘在单位不平衡激励下各测点的振动幅值-转速曲线和相位-转速曲线。任意测点对应的两条曲线即为在1级叶盘激励位置到该测点之间系统的幅频特性曲线与相频特性曲线，直接反映整机系统在1级叶盘单位不平衡激励下的振动响应特性。可以看出，由于转子-机匣结构动力特性的影响，机匣测点振动响应变化规律十分复杂，体现为以下3点：

1) 同一测点振动响应随转速变化明显。以图6(a)为例，截面1的测点2在3 120 r/min及7 140 r/min 转速下的振动幅值及相位均有显著变化。

2) 在同一转速下，同一截面不同方向测点的振动响应存在差异。以图6(b)为例，在3 900 r/min 转速下，位于截面2的测点3和测点4在振动幅值及相位均存在明显差异。

3) 不同截面同一方向测点的振动响应存在差异。以图6(d)为例，同为水平方向的截面4测点7与截面5测点9在3 900 r/min转速下的振动幅值和相位均有明显差异。

图5 机匣振动测点分布Fig.5 Distribution of vibration measuring points of casing

此外，激励位置也对各测点的振动响应产生影响，在不同级叶盘施加单位不平衡激励，同一测点振动响应也存在一定差异。以测点1振动响应为例，1级叶盘施加单位不平衡激励与9级叶盘施加单位不平衡激励，在特定转速或一定转速区间内测点1振动幅值与相位存在明显差异，如图7 所示，在3 960 r/min位置，测点振动幅值差异明显。

图7 1级和9级叶盘施加单位不平衡激励测点1振动响应Fig.7 Vibration response of measuring point 1 excited by unit unbalance vector from 1st blisk and 9th blisk respectively

以上数值仿真结果验证了不平衡激励位置、测点及转速对燃机整机系统不平衡振动响应的影响，在进行整机动平衡时必须予以考虑。

对于式(20)的构建，可依次分析计算第1～9级叶盘分别施加单位不平衡激励时各截面测点振动响应。

以第1级叶盘施加单位不平衡激励机匣外表面测点振动响应为例，图6中每组幅频特性曲线与相频特性曲线分别对应全转速区间不平衡响应系数矩阵SS中第1列向量的幅值和相位。以此类推，在各级叶盘分别施加单位不平衡激励，获取叶盘激励位置到测点之间系统的幅频特性曲线与相频特性曲线，可得到全转速区间不平衡响应系数矩阵SS的第2～9列向量，从而得到完整的全转速区间不平衡响应系数矩阵SS。

3.3 基于仿真的等效不平衡量逆推准确性验证

为验证逆推方法的正确性，在该发动机模型转子的1～9级叶盘上同时随机施加不平衡激励，模拟转子上的随机不平衡量，如表1所示。

将转子转频区间设定为6～175 Hz(对应转速为300～10 500 r/min)，频率间隔为1 Hz，使用ANSYS谐响应法计算发动机模型不平衡激励振动响应。提取机匣各测点振动响应并绘制300～10 500 r/min转速范围内各测点振动的幅值-转速图，如图8所示。

表1 随机不平衡量Table 1 Random unbalance vector

图8 1～9级叶盘同时施加不平衡激励时机匣振动响应Fig.8 Vibration response of casing under unbalance excitation from 1st-9th blisks at the same time

将振动幅值、相位信息以复数形式表示，作为所提出方法的数据输入，配合全转速区间不平衡响应系数矩阵，可以得到全转速等效不平衡量逆推方程组：

SS·q=zz

(25)

考虑到该方法的实际可选配平面及可操作性，按照2.3节所提降维重构方法，经过优化选取5个叶盘作为配平平面，并经配平叶盘位置、机匣测点及转速的筛选后，确定出使用第1、2、3、4、7级叶盘进行动平衡，即将转子的不平衡量集中等效至这5级叶盘上；使用1、2、8测点的振动数据，以及对应的转速360、1 560、5 580 r/min。得到最终的降维重构逆推方程组：

(26)

(27)

式中：ω1、ω21和ω88为分别筛选出的转速，ω1=360 r/min，ω21=1 560 r/min，ω88=5 580 r/min；逆推方程组系数矩阵S中各元素Sqk(ωh)为q(1、2、8)测点、k(1、2、3、4、7)级叶盘在ωh(360、1 560、5 580 r/min)转速下的不平衡响应系数；qdr=[q1,q2,q3,q4,q7]T中各元素分别为等效至第1、2、3、4、7级叶盘的不平衡量；zdr=[z1(ω1),z8(ω1),z1(ω21),z2(ω88)]T中各元素分别为1、2、8测点在360、1 560、5 580 r/min转速下的振动幅值和相位。

求解式(27)，由于方程个数小于未知数个数，故通过响应系数矩阵伪逆求取式(27)的最小范数解[21]，得到列向量qdr的幅值和相位，即等效不平衡量，如表2所示。

在这5级叶盘上，与等效不平衡量等幅值、反相位施加配平矢量进行谐响应分析，得到配平后的机匣各测点振动响应幅值，并与配平前进行对比，如图9所示。可见，配平前机匣各测点振动幅值在10-6～10-3m数量级，而配平后机匣各测点振动幅值在10-12～10-7m数量级，使残余不平衡量幅值接近于0，减振效果十分显著，说明逆推出的不平衡量与原始不平衡量具有高度等效性，同时也证明了该逆推方法的有效性。