CN115077919B - 一种适于航空发动机机载的整机振动评估方法 - Google Patents

一种适于航空发动机机载的整机振动评估方法 Download PDF

Info

Publication number
CN115077919B
CN115077919B CN202210595102.7A CN202210595102A CN115077919B CN 115077919 B CN115077919 B CN 115077919B CN 202210595102 A CN202210595102 A CN 202210595102A CN 115077919 B CN115077919 B CN 115077919B
Authority
CN
China
Prior art keywords
fault
vibration
force
engine
rotor
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210595102.7A
Other languages
English (en)
Other versions
CN115077919A (zh
Inventor
王俨剀
姚尚鹏
高术依
王彤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202210595102.7A priority Critical patent/CN115077919B/zh
Publication of CN115077919A publication Critical patent/CN115077919A/zh
Application granted granted Critical
Publication of CN115077919B publication Critical patent/CN115077919B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M15/00—Testing of engines
    • G01M15/04—Testing internal-combustion engines
    • G01M15/12—Testing internal-combustion engines by monitoring vibrations
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

一种适于航空发动机机载的整机振动评估方法，以降低虚警损失和误诊损失为目标。在航空发动机典型故障的智能诊断中，通过有限元模型进行故障响应仿真得到典型故障的振动响应；提取典型故障的时域以及频域特征作为bp神经网络输入参数。在模型训练过程中提出了考虑虚警损失和误诊损失的损失函数，以该损失函数最小化作为训练目标。并考虑到误诊情况下发动机存在严重实际损失风险，因此引入风险系数εi使得在模型参数迭代优化过程中更多的降低误诊损失。最后采用仿真数据进行模型准确性是否满足诊断要求的验证。本发明可针对不同结构的发动机转子进行建模仿真以及诊断模型训练，具有较为广泛的适用性。

Description

一种适于航空发动机机载的整机振动评估方法
技术领域
本发明涉及航空发动机智能诊断技术领域，具体是一种基于动力学模型的航空发动机的振动状态评估方法。
技术背景
在航空发动机振动状态评估领域，bp神经网络、支持向量机、卷积神经网络等智能诊断方法已有应用先例。其智能化体现在由事先训练完成的诊断模型自行计算得到振动状态评估结果，全程无需技术人员参与。智能诊断的关键在于数学模型与实际问题本质之间的逼近程度。对于机载智能诊断系统而言，确保飞行安全则是重中之重，这就需要诊断过程中将虚警损失与误诊损失降至最低。
在发动机振动状态评估中，整机振动数据能较好的反映故障本质。引起发动机异常振动的因素可达上百种，但对于特定型号的发动机，其典型振动故障种类有限且故障特征可描述。这便为利用振动数据进行智能诊断提供了可行性。
在以往公开的方法中，多数采用模式识别以及最大隶属度原则确定最终诊断结果。
在公开号为CN105758645A的发明创造中公开了一种基于概率神经网络的发动机振动状态评估系统，包括进行历史数据的预处理、特征选择和提取等流程，通过训练样本，得到期望的诊断模型，利用该模型进行诊断。该方法在模型训练过程中并未区分虚警与误诊两种错误诊断情况。
当训练样本特征不明显或存在干扰信号时，易造成虚警或误诊。当两类故障特征相近时，如不平衡故障与叶片掉块故障频域特征均表现为1倍频占优。但叶片掉块故障危害程度远大于不平衡故障，若将叶片掉块故障误诊为不平衡故障，则发动机存在重大安全隐患，这对于机载系统是不可接受的。
中国发明专利为ZL201710230772.8中公开了一种航空发动机结构类故障的智能诊断方法，采用数据分析方法对故障样本进行识别得到疑似故障，再依据故障因子决策表对疑似故障进行筛选得到主要疑似故障，最后采用模式识别确定故障类别。并利用否定检验降低了由干扰信号带来的虚警损失。该方法考虑到了干扰信号带来的虚警损失，并从动力学角度对结果进行了检验，缩小了模型学习范围，降低了虚警损失。该方法考虑到了干扰信号带来的虚警损失影响，但在模型训练过程中还是未能区分由于诊断结果错误带来的误诊损失与虚警损失影响。
大多数智能诊断方法以精确区分故障类别为首要目标，在训练样本选取时，直接选取实验器数据或是发动机实测数据进行特征提取，而后进行模型训练。这样做可以保证诊断模型尽可能贴近真实情况，但由于发动机实测数据存在大量噪声干扰，在模型训练过程中，这些干扰会极大的影响训练效果。致使诊断模型未能准确分辨出故障特征。究其原因就在于未能从动力学角度出发，建立故障模型，得到特征显著的训练数据，在模式识别过程中，未考虑到虚警损失和误诊损失对诊断结果的影响。为达到精确分类的目标，在训练过程中，模型复杂度会不断上升，结果导致模型对训练样本有很好的准确度，而对测试样本准确度大打折扣，也就难以在工程层面应用。
对于地面机械振动状态评估来说，准确率是追求的终极目标；而对于机载的实时状态评估系统来说，既要考虑诊断的准确性，也要考虑虚警损失。这二者实际上代表两类不同的识别损失。
发明内容
为克服现有技术中存在的虚警损失和误诊损失高、训练数据质量不一、难以应用于机载的问题，本发明提出了一种适于航空发动机机载的整机振动评估方法。
本发明的具体过程是：
步骤1，建立单转子有限元模型：
所述建立转子有限元模型的具体做法如下：
第一步。建立子坐标系。包括OXYZ为固定坐标系和Oxyz为旋转坐标系。
转子系统在变形状态下，其任意横截面相对于固定坐标系OXYZ的位置用V,W,B,Γ表示，V表示Y方向的位移，W表示Z方向的位移，B表示绕Y轴的转角，Γ表示绕Z轴的转角。
第二步，推导刚性盘的运动微分方程。
得固定坐标系OXYZ中刚性盘的运动微分方程为：
式中：是刚性盘的质量矩阵及惯性矩阵；Gd是刚性盘的陀螺效应矩阵；Qd是刚性盘的外力向量；d是盘单元上标；q为位移矢量；/>为q一次求导项；/>为q二次求导项。
第三步，确定等截面弹性轴段单元运动微分方程。
所述等截面弹性轴段单元在固定坐标系中的运动微分方程如下：
式中：是等截面轴段单元的质量矩阵以及惯性矩阵；Ge是梁单元的陀螺效应矩阵；Ke是梁单元的刚度矩阵；Qe是梁单元的外力向量；q为位移矢量；/>为q一次求导项；/>为q二次求导项。
得到等截面弹性轴段单元任意截面横向平动位移与横向转角位移的表达式后，通过形函数可将任意微元段坐标用单元坐标qe表示，这样微元段的势能dPe和动能dTe同样可用单元坐标qe表示：
式中：ΘΓ-Ψ′V为YOX平面的剪切变形；-ΘB-Ψ′W为ZOX平面的剪切变形；ρl为单位长度的质量；Id为单位长度的直径转动惯量；Ip为单位长度的直径转动惯量；Φ为转角。
将式(6)、式(7)带入到式(12)、式(13)，并沿元素全长积分得到轴段单元动能与势能表达式：
式中：
利用拉格朗日方程：
所述等截面弹性轴段单元运动微分方程的推导过程是：
在转子系统有限元法中，通常采用考虑了剪切变形的Timoshenko梁单元对弹性轴单元进行建模。每个Timoshenko梁单元有两个节点，如前所述每个节点有4个自由度，即每个梁单元有8个自由度。梁单元的8个自由度组成的广义位移向量为：
qe＝[q1 q2 q3 q4 q5 q6 q7 q8]T (4)
式中：q1,q5为轴元素两端在Y方向的位移；q2,q6为轴元素两端在Z方向的位移；q3,q7为轴元素两端绕Y轴的转角；q4,q8为轴元素两端绕Z轴的转角。
设轴段单元长度为l，轴向距离s处截面的横向位移为(V,W,B,Γ)，位移和转角的关系可近似表示为：
使用轴元素两端节点广义坐标表示该轴元素中任意微元段的位移和转角表示为：
式中：ψ1,ψ2,ψ3,ψ4是平动位移插值函数；θ1,θ2,θ3,θ4是转角位移插值函数。
所述式6与7中，形函数矩阵为：
所述式8与9中，
横向平动位移插值函数和横向转角位移插值函数的表达式为：
式中：
其中：ζ为微元段相对位置；为剪切变形系数；As为有效抗剪面积；A为横截面面积；E为弹性模量；I为截面惯性矩；G为剪切模量；D，d，l分别为元素外径、内径、长度；ν为泊松比。
第四步，计算轴承运动微分方程。
在考虑线性刚度和阻尼时，轴承的运动方程为
式中，Cb为轴承阻尼矩阵；Kb为轴承刚度矩阵；Qb_ex为轴承处外力；上标b表示轴承元素；q为位移矢量；为q一次求导项；/>为q二次求导项。
第五步，组装运动方程。
将盘元素、轴元素和轴承元素的运动方程组装成系统运动方程时，需要将各同类项的系数矩阵相加。
组装后的转子系统稳态运动方程为：
式中，M为转子系统的质量矩阵，C为转子系统的阻尼矩阵，G为转子系统的陀螺效应矩阵，K为转子系统的刚度矩阵，qs为转子系统的位移向量，为qs一次求导项；/>为qs二次求导项。Qs为转子系统的外力向量。
至此，建立了单转子动力学模型，然后将该模型在matlab环境中编程。即可进行下一步的动力学响应计算。
步骤2，求解步骤1所建立的单转子动力学模型的典型故障的动力学响应：
所述单转子动力学模型的典型故障为不平衡故障、不对中故障、转静碰摩故障和叶片掉块故障。
第一步，确定转子系统的不平衡量大小与位置、不平衡故障力表达式、以及不平衡时的发动机转速信息。将故障力施加于转子模型模拟盘位置。利用matlab编程求解单转子动力学模型的振动响应得到n1组不平衡故障时域响应。所述时域响应为不平衡故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n1组不平衡故障情况下发动机振动的频域响应，所述频域响应为转子系统不平衡故障情况下的振动幅值随频率的变化。
在步骤一第5步所组装的运动方程中，模拟盘上的不平衡质量产生的离心力是不平衡振动中周期性激励的主要来源。设转子不平衡质量为m，偏心距离为e，不平衡质量初始相位为θ，若转子角速度为ω，转子转动的时间为t，则不平衡故障水平方向沿x轴的离心力Fx1以及不平衡故障竖直方向沿y轴的离心力Fy1表示为：
Fx1＝meω2cos(ωt-θ)，Fy1＝meω2sin(ωt-θ) (26)
第二步，确定不对中量大小与位置、不对中故障力表达式、不对中时的发动机转速信息，将故障力施加于转子模型联轴器位置。利用matlab编程求解单转子动力学模型的振动响应得到n2组不对中故障时域响应，所述时域响应为不对中故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n2组不对中故障情况下发动机振动的频域响应，所述频域响应为转子系统不对中故障情况下的振动幅值随频率的变化。
在工程实际中，不对中故障经常以综合不对中的情况出现，所以在仿真时以综合不对中情况下的受力情况作为故障力。给出不对中故障力表达式：
mi联轴器质量，△E综合不对中，ω转子旋转速度，t时间，不对中相位，m0初始不平衡量，/>为初始不平衡相位，e偏心距离，Fx2水平方向的故障力，Fy2垂直方向的故障力。
第三步，确定机匣刚度与侵入量、转静碰摩故障力表达式、转静碰摩时的发动机转速信息，将故障力施加于转子模型模拟盘位置。力。利用matlab编程求解单转子动力学模型的振动响应得到n3组转静碰摩故障时域响应，所述时域响应为转静碰摩故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n3组转静碰摩故障情况下发动机振动的频域响应，所述频域响应为转子系统转静碰摩故障情况下的振动幅值随频率的变化。
设发生碰摩时，转子和静子之间发生弹性形变，用法向正压力和切向摩擦力来表示该碰摩时的碰摩力。碰摩时产生的法向弹性力为FN、，碰摩时产生的切向摩擦力为FT，δ为转静间隙，rd为运行时转子与机匣之间的相对位移，为碰撞点法向和水平方向夹角。
设碰摩发生时弹性形变刚度为Kr，摩擦系数为μ。则法向弹性力FN、和切向摩擦力FT分别表示为：
式(28)中，所述x为转盘质心水平方向的振动位移，y为转盘质心竖直方向的振动位移。当rd≤δ时，未发生碰摩，FN、FT均为0，当rd＞δ时，发生转静碰摩，得到转静碰摩故障的时域以及频域响应。
将所述FN、FT分解到x、y方向，可得：FN、FT沿x方向的合力Fx3以及FN、FT沿y方向的合力Fy3
式中，Fx3转静碰摩水平方向故障力，Fy3转静碰摩数值方向故障力，FN法向弹性力，FT为切向摩擦力，转静碰摩故障力相位。
第四步，确定掉块质量大小与位置、叶片掉块故障力表达式、叶片掉块时的发动机转速信息，将故障力施加于转子模型模拟盘位置。模拟n4组掉块质量。将n4组掉块质量代入式(30)得到n4组叶片掉块故障力，令系统外力矩阵等于叶片掉块故障力。利用matlab编程求解单转子动力学模型的振动响应得到n4组叶片掉块故障时域响应。所述时域响应为叶片掉块故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n4组叶片掉块故障情况下发动机振动的频域响应，所述频域响应为转子系统叶片掉块故障情况下的振动幅值随频率的变化。
将叶片掉块视为突加不平衡。其故障力表示为：
式中，Fx4为故障力在水平方向沿x轴的投影，Fy4为故障力在竖直方向沿y轴的投影，t为时间，ts为掉块发生时刻，ms为掉块质量，m0为初始不平衡质量，m1为不平衡质量，为叶片掉块相位，/>为初始不平衡相位
第五步，求解得n5组无故障情况下发动机振动的时域响应。所述时域响应为转子系统在不平衡量低于5g.cm情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n5组无故障情况下发动机振动的频域响应，所述频域响应为转子系统在不平衡量低于5g.cm情况下的振动幅值随频率的变化。
通过以上步骤得到发动机典型故障的整机振动时域响应数据和频域响应数据。
将所述的n1组不平衡故障的时域及频域响应、n2组不对中故障的时域及频域响应、n3组转静碰摩故障的时域及频域响应、n4组叶片掉块故障的时域及频域响应和n5组无故障情况下发动机振动的时域及频域响应，作为bp神经网络训练数据样本。
步骤3，故障响应数据分析处理：
所述故障响应数据分析处理的具体过程是：
第一步，对不平衡、不对中、转静碰摩和叶片掉块四种故障仿真数据进行傅里叶变换，得到各自的故障特征频率。
第二步，对不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障这五种情况下的动力学时域响以及频域响应进行特征提取。所述特征提取为，计算时域响应的方差及均值，以及，记录频域响应的1/2X，1X，2X，3X的幅值。根据第一步中得到的特征频率，取其0.5倍旋转基频、1倍旋转基频、2倍旋转基频、3倍旋转基频，共计四个特征频率的幅值。
对步骤2中得到的各组振动时域及频域响应数据提取：频域内的故障特征频率的0.5倍频率的幅值故障特征频率的1倍频率的幅值A1x、故障特征频率的2倍频率的幅值A2x、故障特征频率的3倍频率的幅值A3x以及时域内的均值Amean和方差s2共计六个特征参数。
第三步，对得到的n1、n2、n3、n4和n5组数据的六个特征参数进行归一化处理。
所述n1、n2、n3、n4和n5组数据之和为W组。
所述归一化处理如下：
记所述n1、n2、n3、n4和n5组数据的1/2x特征频率的幅值分别为归一化处理后所得的n1、n2、n3、n4和n5组数据为：
记n1、n2、n3、n4和n5组数据的1x特征频率的幅值分别为A1xi，i＝1，2，…W。
归一化处理后所得的n1、n2、n3、n4和n5组数据为：
记n1、n2、n3、n4和n5组数据的2x特征频率的幅值分别为A2xi，i＝1，2，…W。
归一化处理后所得的n1、n2、n3、n4和n5组数据为：
记n1、n2、n3、n4和n5组数据的3x特征频率的幅值分别为A3xi，i＝1，2，…W。
归一化处理后所得的W组数据为：
记n1、n2、n3、n4和n5组数据的时域振动幅值的均值Ameani，i＝1，2，…W。归一化处理后所得的100组数据为：
记n1、n2、n3、n4和n5组数据的时域振动幅值的方差归一化处理后所得的W组数据为：
步骤4，建立振动状态评估模型：
利用仿真数据进行bp神经网络模型训练，得到振动状态评估模型。
所述利用仿真数据进行bp神经网络模型训练的具体过程是：
第一步，建立故障样本：对于步骤3中得到的不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障这五种情况下的归一化后的n1、n2、n3、n4和n5组特征参数数据，依据“不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障”五种类型均分为5组。随机抽取每组的3/4个样本作为训练样本，剩余的1/4个作为测试样本。
建立故障样本时的打标签处理为：无故障的故障类型为正常，标记为1。不平衡类型与不对中类型均为异常，标记为2。转静碰摩类型与叶片掉块类型均为故障，标记为3。
第二步，给定损失函数表达式。
式(33)中，x为训练次数，Loss(x)为训练误差，εi为风险系数。y为故障标签所标识的值，取值范围为1，2，3。y*为模型计算结果；当y＝1时，发动机实际状态为正常。若此时的诊断结果为y＝2或y＝3，则称所述情况为虚警。当y＝2时，发动机实际状态为异常，当y＝3时，发动机实际状态为故障。当y＝2以及y＝3时，若此时的诊断结果为y＝1，则称所述情况为误诊。
εi为风险系数，振动状态评估中误诊带来的风险大于虚警情况下的风险，因此在损失函数计算中要考虑到误诊带来的风险占比要高于虚警带来的风险。
当实际发动机状态为1正常时，若bp神经网络模型计算结果为2异常或3故障，此时发动机实际状态为健康，而模型计算结果为异常或故障，此时诊断结果为虚警。此时εi取值为1。
当实际发动机状态为2异常时，若bp神经网络模型计算的结果为1正常，此时诊断结果为误诊且发动机存在故障风险，风险系数εi取值为2。若模型计算结果为故障，此时诊断结果为误诊且发动机实际故障风险低于误诊结果状态下的故障风险，风险系数εi取值为1。
当实际发动机状态为3故障时，若bp神经网络模型计算的结果为1正常，此时诊断结果为误诊且发动机存在严重故障风险，风险系数εi取值为8。若bp神经网络模型计算的结果为2异常，此时诊断结果为误诊且发动机实际故障风险高于误诊结果状态下的故障风险，风险系数εi取值为4。
辨识的四种故障在等级V故障状态中的隶属度值。将不平衡与不对中定义为异常状态，将转静碰摩与叶片掉块定义为故障状态，据此给出风险系数εi的取值。
所述风险系数εi的取值为：
当正常识别为故障或异常时，εi的取值为1；
当异常识别为故障时，εi的取值为1；
当异常识别为正常时，εi的取值为2；
当故障识别为异常时，εi的取值为4；
当故障识别为正常时，εi的取值为8。
第三步，计算bp神经网络模型的损失，
利用梯度下降法计算所述bp神经网络模型的损失时的具体过程是：
Ⅰ循环计算从i＝0直至i＝训练数据的长度：
ⅰ计算第i个训练数据的权重ω和偏差b相对于损失函数的梯度于是得到每一个训练数据的权重和偏差的梯度值。
ⅱ计算所有训练数据权重ω的梯度的总和。
ⅲ计算所有训练数据偏差b的梯度的总和。
Ⅱ更新各样本的权重值和偏差值
ⅰ使用上面第(2)、(3)步所得到的结果，计算所有样本的权重和偏差的梯度的平均值。
ⅱ通过公式(34)，更新每个样本的权重值和偏差值。
式中：α为学习率，ωi为第i次计算得到的权重值，bi为第i次计算得到的偏差值，L为损失函数。
重复所述循环计算从i＝0直至i＝训练数据的长度和更新各样本的权重值和偏差值的过程，直至损失函数收敛不变。
所述损失函数为：由式{31}运算所得的训练误差Loss(x)值随训练次数x值的变化函数。
第四步，输出模型参数，所述模型参数为bp神经网络中的权值。bp神经网络模型训练完成，所述的训练完成了的bp神经网络模型为振动状态评估模型。
至此，完成适对合航空发动机机载的整机振动状态评估。
本发明以降低虚警损失和误诊损失为目标。用于航空发动机典型故障的智能诊断，首先建立转子动力学模型，通过有限元模型进行故障响应仿真得到典型故障的振动响应。然后通过数据预处理提取典型故障的时域以及频域特征作为bp神经网络输入参数。在模型训练过程中提出了考虑虚警损失和误诊损失的损失函数，以该损失函数最小化作为训练目标。并考虑到误诊情况下发动机存在严重实际损失风险，因此引入风险系数εi使得在模型参数迭代优化过程中更多的降低误诊损失。最后采用仿真数据进行模型准确性是否满足诊断要求的验证。
训练过程与现有的神经网络训练过程区别在于，
Ⅰ损失函数的构造不同。现有的损失函数是单纯从数学角度出发定义的，而本发明中采用的损失函数有实际物理意义，表现为虚警损失和误诊损失。
Ⅱ训练样本选取不同。现有智能诊断模型训练数据多数直接采用实验器数据，由于实验器数据响应较为复杂，存在诸多噪声干扰，这些干扰会导致模型难以快速定位到故障特征。本发明建立转子有限元模型，利用特征明显的仿真数据进行神经网络模型训练。
本发明可针对不同结构的发动机转子进行建模仿真以及诊断模型训练，具有较为广泛的适用性。
与现有技术相比较，本发明的有益效果为：
1、本发明利用特征显著的仿真数据进行模型训练，降低了直接采用实验器数据时存在的噪声干扰问题的影响。采用特征显著的仿真数据使模型快速学习到典型故障的故障特征。
2、本发明提出的考虑虚警损失和误诊损失的损失函数可以最大限度地保证发动机安全，常规损失函数对这两者并未做区分，但应用到实际发动机时，“误诊”代表发动机存在故障，但诊断结果并未准确展现出来，此时发动机本身处于故障或异常状态，若将故障状态误诊为正常状态，其潜在危害不言而喻。这对于机载振动状态评估系统是不可接受的。
3、本发明在损失函数计算过程中引入的风险系数εi可以使得模型内部参数优化过程中更加关注减低误诊损失。“虚警”代表将正常识别为故障，此时发动机实际状态是正常，发动机不存在故障风险；而“误诊”代表将故障识别为正常，此时发动机实际状态为故障，存在严重故障风险，因此需要风险系数将二者区分。
4.本发明所述流程中不存在复杂的数据优化算法或大规模的神经网络，而是从故障本质即动力学特征出发，进行故障特征仿真，bp神经网络训练以及最终振动状态评估。除了所示出的不平衡，不对中，转静碰摩，叶片掉块四种典型故障外，还可用于其它能够刻画出动力学模型的典型故障，如齿轮故障，轴承故障，支座松动，盘腔积液。
附图说明
图1是本发明技术方案示意图。
图2是转子模型示意图。
图3是转子坐标系。
图4是轴元素示意图。
图5是故障时域响应以及频域响应图；其中，图5a是时域响应图。图5b是频域响应图。
图6是损失函数训练结果图。
图7是振动状态评估模型验证结果。
图8是本发明的流程图。
图中：1—刚性盘、2—联轴器、3—等截面弹性轴段、4—轴承结构、5—虚警损失、6—误诊损失、7—实际振动类别、8—振动状态评估模型计算得到的振动类别。
具体实施方式
本实施例是一种航空发动机机载的整机振动评估方法，其具体过程是：
步骤1，建立单转子有限元模型：
建立转子动力学模型，以保证模型辨识故障状态的准确性。通过转子模型进行典型故障的仿真，得到特征显著且无噪声干扰的高质量转子振动响应数据。将其作为神经网络训练数据有助于快速达到训练目标。
根据廖明夫等人编著的《航空发动机转子动力学》p208页所述有限元法建立转子模型的流程，建立故障仿真单转子发动机有限元模型。如图2所示。
所述建立单转子发动机有限元模型的具体过程是：
第一步。建立如图3所示的转子坐标系。
所述转子坐标系包括OXYZ固定坐标系和Oxyz旋转坐标系。其中，所述OXYZ固定坐标系与Oxyz旋转坐标系中的x轴与X轴重合，为轴承4中心的连线。Oxyz坐标系通过OXYZ坐标系绕着x轴旋转ωt角度得到，其中ω为转子旋转角速度。
转子系统在变形状态下，其任意横截面相对于固定坐标系OXYZ的位置用V,W,B,Γ表示，V为Y方向的位移，W为Z方向的位移，B为绕Y轴的转角，Γ为绕Z轴的转角。
第二步，推导刚性盘单元运动微分方程。
刚性盘1的动能表示为：
式中：Td是刚性盘的动能；是刚性盘盘心Y方向的平动位移求导；/>是刚性盘盘心Z方向的平动位移求导；/>是刚性盘的绕Y轴的转角求导；/>是刚性盘绕Z轴的转角求导；md是刚性盘质量；Id是刚性盘绕直径转动惯量；Ip是刚性盘极转动惯量；Ω是刚性盘自转角速度。
将刚性盘1的动能表达式(1)代入拉格朗日方程：
可求得固定坐标系OXYZ中刚性盘的运动微分方程为：
式中：分别是刚性盘的质量矩阵以及惯性矩阵；Gd是刚性盘的陀螺效应矩阵；Qd是刚性盘的外力向量；d是盘单元上标；q为位移矢量；/>为q一次求导项；/>为q二次求导项。
其中，刚性盘的质量矩阵为对称矩阵，陀螺效应矩阵为反对称矩阵，他们的表达式分别为：
/>
第三步，推导等截面弹性轴段单元运动微分方程。
在转子系统有限元法中，通常采用考虑了剪切变形的Timoshenko梁单元对等截面弹性轴3的单元进行建模。每个Timoshenko梁单元有两个节点，如前所述每个节点有4个自由度，即每个梁单元有8个自由度。梁单元的8个自由度组成的广义位移向量为：
qe＝[q1 q2 q3 q4 q5 q6 q7 q8]T (4)
式中：q1,q5为轴元素两端在Y方向的位移；q2,q6为轴元素两端在Z方向的位移；q3,q7为轴元素两端绕Y轴的转角；q4,q8为轴元素两端绕Z轴的转角。
如图4所示，设等截面弹性轴段单元长度为l，轴向距离s处截面的横向位移为(V,W,B,Γ)，位移和转角的关系可近似表示为：
使用等截面弹性轴段单元两端节点广义坐标表示该轴元素中任意微元段的位移和转角可表示为：
式中：ψ1,ψ2,ψ3,ψ4是平动位移插值函数；θ1,θ2,θ3,θ4是转角位移插值函数。
形函数矩阵为：
横向平动位移插值函数和横向转角位移插值函数的表达式为：
式中：
其中：ζ为微元段相对位置；为剪切变形系数；As为有效抗剪面积；A为横截面面积；E为弹性模量；I为截面惯性矩；G为剪切模量；D，d，l分别为元素外径、内径、长度；ν为泊松比。
得到等截面弹性轴段单元任意截面横向平动位移与横向转角位移的表达式后，通过形函数可将任意微元段坐标用单元坐标qe表示，这样微元段的势能dPe和动能dTe同样可用单元坐标qe表示：
式中：ΘΓ-Ψ′V为YOX平面的剪切变形；-ΘB-Ψ′W为ZOX平面的剪切变形；ρl为单位长度的质量；Id为单位长度的直径转动惯量；Ip为单位长度的极转动惯量；Φ为转角；为Φ求导所得。
将式(6)、式(7)带入到式(12)、式(13)，并沿元素全长积分得到等截面弹性轴段单元动能与势能表达式：
式中：
利用拉格朗日方程：
得到等截面弹性轴段单元在固定坐标系中的运动微分方程如下：
式中：是梁单元的质量矩阵；Ge是梁单元的陀螺效应矩阵；Ke是梁单元的刚度矩阵；Qe是梁单元的外力向量；q为位移矢量；/>为q一次求导项；/>为q二次求导项。
等截面弹性轴段单元系数矩阵表达式为：
质量矩阵为：
式中: />
式中：陀螺效应矩阵为：
式中：G1＝36，刚度矩阵为：
式中：KB1＝12，KB4＝6l。
第四步，推导轴承运动微分方程。
在考虑线性刚度和阻尼时，轴承4的运动方程为
式中，Cb为轴承阻尼矩阵；Kb为轴承刚度矩阵；Qb_ex为轴承处外力；上标b表示轴承元素。
第五步，组装运动方程。
将刚性盘元素、等截面弹性轴段元素和轴承元素的运动方程组装成系统运动方程时，需要将个同类项的系数矩阵相加。
对于等截面弹性轴段元素，由于相邻两个元素共用一个节点，故需要将共用的部分对应的系数矩阵叠加，而且它们之间相互作用的剪力相互抵消，方程右边只有外力作用；对于刚性盘元素，只需要把相应节点的刚性盘元素的各系数矩阵叠加到等截面弹性轴段轴元素系数矩阵中即可；对于轴承元素，其方程右边为轴承所受的外力，为了处理方便，组装时可以将其移到系统运动方程的左边，并用刚度和阻尼的表达式代替。这样只需把相应节点的刚度和阻尼系数矩阵叠加到等截面弹性轴段元素系数矩阵对应位置即可。
组装后的转子系统稳态运动方程为：
式中，M为转子系统的质量矩阵，C为转子系统的阻尼矩阵，G为转子系统的陀螺效应矩阵，K为转子系统的刚度矩阵，qs为转子系统的位移向量，为q一次求导项；/>为q二次求导项，Qs为转子系统的外力向量。
至此，建立了单转子动力学模型，将该单转子动力学模型在matlab环境中编程。即可进行下一步的动力学响应计算。
建立的转子模型见图2。
步骤2，求解步骤1所建立的单转子动力学模型的典型故障的动力学响应：
所述整机振动的典型故障为不平衡故障、不对中故障、转静碰摩故障和叶片掉块故障。
第一步，确定转子系统的不平衡量大小与位置、不平衡故障力表达式、以及不平衡时的发动机转速信息。将故障力施加于转子模型模拟盘，即附图2中的1位置。模拟n1组不平衡质量：5g·cm，5.1g·cm，5.2g·cm，5.3g·cm，5.4g·cm，5.5g·cm，5.6g·cm，5.7g·cm，5.8g·cm，5.9g·cm，6.0g·cm，6.1g·cm，6.2g·cm，6.3g·cm，6.4g·cm，6.5g·cm，6.6g·cm，6.7g·cm，6.8g·cm，6.9g·cm，7.0g·cm。将其代入式(26)得到n1组不平衡故障力，令系统外力矩阵等于不平衡故障力。利用matlab编程求解单转子动力学模型的振动响应得到n1组不平衡故障时域响应。所述时域响应为不平衡故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n1组不平衡故障情况下发动机振动的频域响应，所述频域响应为转子系统不平衡故障情况下的振动幅值随频率的变化。本实施例中，所述n1＝20。
在步骤一第5步所组装的运动方程中，模拟盘上的不平衡质量产生的离心力是不平衡振动中周期性激励的主要来源。设转子不平衡质量为m，偏心距离为e，不平衡质量初始相位为θ，若转子角速度为ω，转子转动的时间为t，则不平衡故障水平方向沿x轴的离心力Fx1以及不平衡故障竖直方向沿y轴的离心力Fy1表示为：
Fx1＝meω2cos(ωt-θ)，Fy1＝meω2sin(ωt-θ) (26)
第二步，确定不对中量大小与位置、不对中故障力表达式、不对中时的发动机转速信息，将故障力施加于转子模型联轴器2的位置。
模拟n2组综合不对中量为：
0.0019m，0.0020m，0.0021m，0.0022m，0.0023m，0.0024m，0.0025m，0.0026m，0.0027m，0.0028m，0.0029m，0.0030m，0.0031m，0.0032m，0.0033m，0.0034m，0.0035m，0.0036m，0.0037m，0.0038m，将其代入式(27)得到n2组不对中故障力，令系统外力矩阵等于故障力。利用matlab编程求解单转子动力学模型的振动响应得到n2组不对中故障时域响应，所述时域响应为不对中故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n2组不对中故障情况下发动机振动的频域响应，所述频域响应为转子系统不对中故障情况下的振动幅值随频率的变化。本实施例中，所述n2＝20。
在工程实际中，不对中故障经常以综合不对中的情况出现，所以在仿真时以综合不对中情况下的受力情况作为故障力。给出不对中故障力表达式：
mi为联轴器质量，△E为综合不对中量，ω为转子角速度，t为时间，为不对中相位，m0为初始不平衡量，/>为初始不平衡相位，e为偏心距离，Fx2为水平方向的不对中故障力，Fy2为垂直方向的不对中故障力。
第三步，确定机匣刚度与侵入量、转静碰摩故障力表达式、转静碰摩时的发动机转速信息，将故障力施加于转子模型模拟盘位置。模拟n3组侵入量：180mm，181mm，182mm，183mm，184mm，185mm，186mm，187mm，188mm，189mm，190mm，191mm，192mm，193mm，194mm，195mm，196mm，197mm，198mm，199mm。将其代入式(28)得到n3组转静碰摩故障力，令系统外力矩阵等于转静碰摩故障力。利用matlab编程求解单转子动力学模型的振动响应得到n3组转静碰摩故障时域响应，所述时域响应为转静碰摩故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n3组转静碰摩故障情况下发动机振动的频域响应，所述频域响应为转子系统转静碰摩故障情况下的振动幅值随频率的变化。本实施例中，所述n3＝20。
设发生碰摩时，转子和静子之间发生弹性形变，用法向正压力和切向摩擦力来表示该碰摩时的碰摩力。碰摩时产生的法向弹性力为FN、，碰摩时产生的切向摩擦力为FT，δ为转静间隙，rd为运行时转子与机匣之间的相对位移。
设碰摩发生时弹性形变刚度为Kr，摩擦系数为μ。则法向弹性力FN、和切向摩擦力FT分别表示为：
式(28)中，所述x为转盘质心水平方向的振动位移，y为转盘质心竖直方向的振动位移。当rd≤δ时，未发生碰摩，FN、FT均为0，当rd＞δ时，发生转静碰摩，得到转静碰摩故障的时域以及频域响应。/>
将所述FN、FT分解到x、y方向，可得：FN、FT沿x方向的合力Fx3以及FN、FT沿y方向的合力Fy3
式中，Fx3转静碰摩水平方向故障力，Fy3转静碰摩数值方向故障力，FN法向弹性力，FT为切向摩擦力，为转静碰摩故障力相位。
第四步，确定掉块质量大小与位置、叶片掉块故障力表达式、叶片掉块时的发动机转速信息，将故障力施加于转子模型模拟盘位置。模拟n4组掉块质量：20(g.cm)，21(g.cm)，22(g.cm)，23(g.cm)，24(g.cm)，25(g.cm)，26(g.cm)，27(g.cm)，28(g.cm)，29(g.cm)，30(g.cm)，31(g.cm)，32(g.cm)，33(g.cm)，34(g.cm)，35(g.cm)，36(g.cm)，37(g.cm)，38(g.cm)，39(g.cm)。将其代入式(30)得到n4组叶片掉块故障力，令系统外力矩阵等于叶片掉块故障力。利用matlab编程求解单转子动力学模型的振动响应得到n4组叶片掉块故障时域响应。所述时域响应为叶片掉块故障情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n4组叶片掉块故障情况下发动机振动的频域响应，所述频域响应为转子系统叶片掉块故障情况下的振动幅值随频率的变化。本实施例中，所述n4＝20。
由于叶片掉块是瞬间发生的，因此将叶片掉块视为是突加不平衡。其故障力表示为：
式中，Fx4为故障力在水平方向沿x轴的投影，Fy4为故障力在竖直方向沿y轴的投影，t为时间，ts为掉块发生时刻，ms为掉块质量，m0为初始不平衡质量，m1为不平衡质量，为叶片掉块相位，/>为初始不平衡相位
第五步，确定初始不平衡质量大小与位置、无故障情况下外力表达式、无故障时的发动机转速信息，将无故障情况下外力施加于转子模型模拟盘位置。模拟n5组初始不平衡质量：1(g.cm)，1.1(g.cm)，1.2(g.cm)，1.3(g.cm)，1.4(g.cm)，1.5(g.cm)，1.6(g.cm)，1.7(g.cm)，1.8(g.cm)，1.9(g.cm)，2(g.cm)，2.1(g.cm)，2.2(g.cm)，2.3(g.cm)，2.4(g.cm)，2.5(g.cm)，2.6(g.cm)，2.7(g.cm)，2.8(g.cm)，2.9(g.cm)。将其代入式(30)得到n5组叶片掉块故障力，令系统外力矩阵等于叶片掉块故障力。利用matlab编程求解单转子动力学模型的振动响应求解得n5组无故障情况下发动机振动的时域响应。所述时域响应为转子系统在不平衡量低于5g.cm情况下的振动幅值随时间的变化。对所述时域响应进行傅里叶变换，得到n5组无故障情况下发动机振动的频域响应，所述频域响应为转子系统在不平衡量低于5g.cm情况下的振动幅值随频率的变化。本实施例中，所述n5＝20。
式中，Fx5为无故障情况下外力在水平方向沿x轴的投影，Fy5为无故障情况下外力在竖直方向沿y轴的投影，m0为初始不平衡质量，为初始不平衡相位。/>
通过以上步骤得到发动机典型故障的整机振动时域响应数据和频域响应数据。即：上述所得的n1组不平衡故障的时域及频域响应、n2组不对中故障的时域及频域响应、n3组转静碰摩故障的时域及频域响应、n4组叶片掉块故障的时域及频域响应和n5组无故障情况下发动机振动的时域及频域响应，共计W组时域及频域响应数据，其中W＝n1+n2+n3+n4+n5。将其作为bp神经网络训练数据样本。故障响应如附图5所示。本实施例中，W＝100
步骤3，故障响应数据分析处理：
所述故障响应数据分析处理的具体过程是：
第一步，对不平衡、不对中、转静碰摩和叶片掉块四种故障仿真数据进行傅里叶变换，得到各自的故障特征频率。
表1四种典型故障特征频率
故障名称 不平衡 不对中 转静碰摩 叶片掉块
特征频率 1X 1X，2X nX，1/nX 1X
所述故障特征频率，是在频谱图中，当某一故障出现时，频域响应图上会出现与故障类型相对应的频率谱线，则这个故障对应的频率就叫做该故障的故障特征频率。
第二步，对不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障这五种情况下的的动力学时域响以及频域响应进行特征提取。所述特征提取为，计算时域响应的方差及均值，以及频域响应的1/2X，1X，2X，3X的幅值。根据第一步中得到的特征频率，取其0.5倍旋转基频、1倍旋转基频、2倍旋转基频、3倍旋转基频，共计四个特征频率的幅值。
对步骤2中100组振动时域及频域响应数据提取：频域内的故障特征频率的0.5倍频率的幅值故障特征频率的1倍频率的幅值A1x、故障特征频率的2倍频率的幅值A2x、故障特征频率的3倍频率的幅值A3x以及时域内的均值Amean和方差s2这六个特征参数。
第三步，对得到的N组数据的第二步中所述的六个特征参数进行归一化处理。
所述归一化处理如下：
记W组数据的1/2x特征频率的幅值分别为归一化处理后所得的W组数据为：
记W组数据的1x特征频率的幅值分别为A1xi，i＝1，2，…W。归一化处理后所得的W组数据为：
记W组数据的2x特征频率的幅值分别为A2xi，i＝1，2，…W。归一化处理后所得的W组数据为：
记W组数据的3x特征频率的幅值分别为A3xi，i＝1，2，…W。归一化处理后所得的W组数据为：
记W组数据的时域振动幅值的均值Ameani，i＝1，2，…W。归一化处理后所得的W组数据为：
记W组数据的时域振动幅值的方差s2i，i＝1，2，…W。归一化处理后所得的W组数据为：
步骤4，建立初步振动状态评估模型：
利用仿真数据进行bp神经网络模型训练，得到初步振动状态评估模型。
所述利用仿真数据进行bp神经网络模型训练的具体过程是：
第一步，建立故障样本:对于步骤3中得到的不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障这五种情况下的归一化后的W组特征参数数据，依据“不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障”这五种类型，分为5组。随机抽取每组的3/4个样本作为训练样本，其余剩下的1/4作为测试样本。然后进行打标签处理，所述打标签处理为：无故障的故障类型为正常，标记为1。不平衡与不对中的故障类型为异常，标记为2。转静碰摩叶片掉块的故障类型为故障，标记为3。如表2所示。
表2故障数据样本
第二步，给定损失函数表达式。
式(33)中，x为训练次数，Loss(x)为训练误差，εi为风险系数。y为故障标签所标识的值，取值范围为1，2，3。y*为模型计算结果；当y＝1时，发动机实际状态为正常。若此时的诊断结果为y＝2或y＝3，则称所述情况为虚警损失5。当y＝2时，发动机实际状态为异常，当y＝3时，发动机实际状态为故障。当y＝2以及y＝3时，若此时的诊断结果为y＝1，则称所述情况为误诊损失6。
εi为风险系数，振动状态评估中误诊带来的风险大于虚警情况下的风险，因此在损失函数计算中要考虑到误诊带来的风险占比要高于虚警带来的风险。
当实际发动机状态为1正常时，若bp神经网络模型计算结果为2异常或3故障，此时发动机实际状态为健康，而模型计算结果为异常或故障，此时诊断结果为虚警。此时εi取值为1。
当实际发动机状态为2异常时，若bp神经网络模型计算的结果为1正常，此时诊断结果为误诊且发动机存在故障风险，风险系数εi取值为2。若模型计算结果为故障，此时诊断结果为误诊且发动机实际故障风险低于误诊结果状态下的故障风险，风险系数εi取值为1。
当实际发动机状态为3故障时，若bp神经网络模型计算的结果为1正常，此时诊断结果为误诊且发动机存在严重故障风险，风险系数εi取值为8。若bp神经网络模型计算的结果为2异常，此时诊断结果为误诊且发动机实际故障风险高于误诊结果状态下的故障风险，风险系数εi取值为4。
根据王俨剀等人编著的《航空发动机振动状态评估》P418页给出的不同损失风险故障对5个健康状态等级的隶属度值一表。可以得到本发明要辨识的四种故障在等级V故障状态中的隶属度值。本发明中将不平衡与不对中定义为异常状态，将转静碰摩与叶片掉块定义为故障状态正是基于此表。基于此数据给出风险系数εi的取值。如表3所示：
表3风险系数取值
第三步，计算bp神经网络模型的损失，利用梯度下降法：
Ⅰ循环计算从i＝0直至i＝训练数据的长度：
i计算第i个训练数据的权重ω和偏差b相对于损失函数的梯度于是得到每一个训练数据的权重和偏差的梯度值。
ii计算所有训练数据权重ω的梯度的总和。
iii计算所有训练数据偏差b的梯度的总和。
Ⅱ.更新权值与偏差：
i使用I中第ii、iii步所得到的结果，计算所有样本的权重和偏差的梯度的平均值。
ii通过公式(34)，更新每个样本的权重值和偏差值
式中：α为学习率，ωi为第i次计算得到的权重值，bi为第i次计算得到的偏差值，L为损失函数。
重复所述循环计算上面的过程，直至损失函数收敛不变。
损失函数曲线如附图6所示。所述损失函数为：由式(33)运算所得的训练误差Loss(x)值随训练次数x值的变化函数。
第四步，输出模型参数，所述模型参数为bp神经网络中的权值。bp神经网络模型训练完成，所述的训练完成了的bp神经网络模型为振动状态评估模型。
至此，一种适合航空发动机机载的整机振动状态评估方法完成。
为了验证振动状态评估效果，进行如下流程。
所述振动状态评估效果验证，即校验所得到的振动状态评估模型的准确性
根据步骤二故障响应仿真流程，针对正常情况设置不同的初始不平衡量大小共计10组；针对不平衡故障设置不同的不平衡量大小共计10组；针对不对中故障设置不同的不对中量大小共计10组；针对转静碰摩故障设置不同的侵入量大小共计10组；针对叶片掉块故障设置不同的掉块质量大小共计10组；共计50组故障响应仿真数据作为验证数据。按照步骤3中所述数据处理方法进行数据预处理后代入步骤4得到的振动状态评估模型，得到振动状态评估结果。评估结果如附图7所示。附图7中分别是实际振动类别7和振动状态评估模型计算得到的振动类别8。

Claims (7)

1.一种适合航空发动机机载的整机振动评估方法，其特征在于，
步骤1，建立单转子有限元模型：
所述建立转子有限元模型的具体做法如下：
第一步；建立子坐标系；包括OXYZ为固定坐标系和Oxyz为旋转坐标系；
转子系统在变形状态下，其任意横截面相对于固定坐标系OXYZ的位置用V,W,B,Γ表示，V表示Y方向的位移，W表示Z方向的位移，B表示绕Y轴的转角，Γ表示绕Z轴的转角；
第二步，推导盘单元运动微分方程；
得固定坐标系OXYZ中刚性盘的运动微分方程为：
式中：是刚性盘的质量矩阵及惯性矩阵；Gd是刚性盘的陀螺效应矩阵；Qd是刚性盘的外力向量；d是盘单元上标；q为位移矢量；/>为q一次求导项；/>为q二次求导项；
第三步，确定等截面弹性轴段单元运动微分方程；
所述等截面弹性轴段在固定坐标系中的运动微分方程如下：
式中：是等截面单元的质量矩阵；Ge是梁单元的陀螺效应矩阵；Ke是梁单元的刚度矩阵；Qe是梁单元的外力向量；q为位移矢量；/>为q一次求导项；/>为q二次求导项
得到轴段单元任意截面横向平动位移与横向转角位移的表达式后，通过形函数可将任意微元段坐标用单元坐标qe表示，这样微元段的势能dPe和动能dTe同样可用单元坐标qe表示：
式中：ΘΓ-Ψ′V为YOX平面的剪切变形；-ΘB-Ψ′W为ZOX平面的剪切变形；ρl为单位长度的质量；Id为单位长度的直径转动惯量；Ip为单位长度的直径转动惯量；Φ为转角；
第四步，计算普通轴承运动微分方程；
在考虑线性刚度和阻尼时，轴承的运动方程为
式中，Cb为轴承阻尼矩阵；Kb为轴承刚度矩阵；Qb_ex为轴承处外力；上标b表示普通轴承元素；
第五步，组装运动方程；
将盘元素、轴元素和轴承元素的运动方程组装成系统运动方程时，需要将各同类项的系数矩阵相加；
组装后的转子系统稳态运动方程为：
式中，M为转子系统的质量矩阵，C为转子系统的阻尼矩阵，G为转子系统的陀螺效应矩阵，K为转子系统的刚度矩阵，qs为转子系统的位移向量，Qs为转子系统的外力向量；
至此，建立了单转子动力学模型，然后将该模型在matlab环境中编程；即可进行下一步的动力学响应计算；
步骤2，求解步骤1所建立的单转子模型的的典型故障的动力学响应：
所述整机振动的典型故障为不平衡故障、不对中故障、转静碰摩故障和叶片掉块故障；
第一步，确定转子系统的不平衡量大小与位置、不平衡故障力表达式、以及不平衡时的发动机转速信息；将故障力施加于转子模型模拟盘位置；利用matlab编程求解单转子动力学模型的振动响应得到20组不平衡故障时域响应；所述时域响应为不平衡故障情况下的振动幅值随时间的变化；对所述时域响应进行傅里叶变换，得到n1组不平衡故障情况下发动机振动的频域响应，所述频域响应为转子系统不平衡故障情况下的振动幅值随频率的变化；
在步骤一第5步所组装的运动方程中，模拟盘上的不平衡质量产生的离心力是不平衡振动中周期性激励的主要来源；设转子不平衡质量为m，偏心距离为e，不平衡质量初始相位为θ，若转子角速度为ω，转子转动的时间为t，则不平衡故障水平方向沿x轴的离心力Fx1以及不平衡故障竖直方向沿y轴的离心力Fy1表示为：
Fx1＝meω2cos(ωt-θ)，Fy1＝meω2sin(ωt-θ) (26)
第二步，确定不对中量大小与位置、不对中故障力表达式、不对中时的发动机转速信息，将故障力施加于转子模型联轴器位置；利用matlab编程求解单转子动力学模型的振动响应得到n2组不对中故障时域响应，所述时域响应为不对中故障情况下的振动幅值随时间的变化；对所述时域响应进行傅里叶变换，得到20组不对中故障情况下发动机振动的频域响应，所述频域响应为转子系统不对中故障情况下的振动幅值随频率的变化；
在工程实际中，不对中故障经常以综合不对中的情况出现，所以在仿真时以综合不对中情况下的受力情况作为故障力；给出不对中故障力表达式：
mi联轴器质量，△E综合不对中，ω转子旋转速度，t时间，不对中相位，m0初始不平衡量，/>为初始不平衡相位，e偏心距离，Fx2水平方向的故障力，Fy2垂直方向的故障力；
第三步，确定机匣刚度与侵入量、转静碰摩故障力表达式、转静碰摩时的发动机转速信息，将故障力施加于转子模型模拟盘位置；力；利用matlab编程求解单转子动力学模型的振动响应得到n3组转静碰摩故障时域响应，所述时域响应为转静碰摩故障情况下的振动幅值随时间的变化；对所述时域响应进行傅里叶变换，得到n3组转静碰摩故障情况下发动机振动的频域响应，所述频域响应为转子系统转静碰摩故障情况下的振动幅值随频率的变化；
设发生碰摩时，转子和静子之间发生弹性形变，用法向正压力和切向摩擦力来表示该碰摩时的碰摩力；碰摩时产生的法向弹性力为FN、，碰摩时产生的切向摩擦力为FT，δ为转静间隙，rd为运行时转子与机匣之间的相对位移，为碰撞点法向和水平方向夹角；
设碰摩发生时弹性形变刚度为Kr，摩擦系数为μ；则法向弹性力FN、和切向摩擦力FT分别表示为：
式(28)中，所述x为转盘质心水平方向的振动位移，y为转盘质心竖直方向的振动位移；当rd≤δ时，未发生碰摩，FN、FT均为0，当rd＞δ时，发生转静碰摩，得到转静碰摩故障的时域以及频域响应；
将所述FN、FT分解到x、y方向，可得：FN、FT沿x方向的合力Fx3以及FN、FT沿y方向的合力Fy3
式中，Fx3转静碰摩水平方向故障力，Fy3转静碰摩数值方向故障力，FN法向弹性力，FT为切向摩擦力，转静碰摩故障力相位；
第四步，确定掉块质量大小与位置、叶片掉块故障力表达式、叶片掉块时的发动机转速信息，将故障力施加于转子模型模拟盘位置；模拟n4组掉块质量；将n4组掉块质量代入式(30)得到n4组叶片掉块故障力，令系统外力矩阵等于叶片掉块故障力；利用matlab编程求解单转子动力学模型的振动响应得到n4组叶片掉块故障时域响应；所述时域响应为叶片掉块故障情况下的振动幅值随时间的变化；对所述时域响应进行傅里叶变换，得到n4组叶片掉块故障情况下发动机振动的频域响应，所述频域响应为转子系统叶片掉块故障情况下的振动幅值随频率的变化；
将叶片掉块视为突加不平衡；其故障力表示为：
式中，Fx4为故障力在水平方向沿x轴的投影，Fy4为故障力在竖直方向沿y轴的投影，t为时间，ts为掉块发生时刻，ms为掉块质量，m0为初始不平衡质量，m1为不平衡质量，为叶片掉块相位，/>为初始不平衡相位
第五步，求解得n5组无故障情况下发动机振动的时域响应；所述时域响应为转子系统在不平衡量低于5g.cm情况下的振动幅值随时间的变化；对所述时域响应进行傅里叶变换，得到n5组无故障情况下发动机振动的频域响应，所述频域响应为转子系统在不平衡量低于5g.cm情况下的振动幅值随频率的变化；
通过以上步骤得到发动机典型故障的整机振动时域响应数据和频域响应数据；
步骤3，故障响应数据分析处理：
所述故障响应数据分析处理的具体过程是：
第一步，对不平衡、不对中、转静碰摩和叶片掉块四种故障仿真数据进行傅里叶变换，得到各自的故障特征频率；
第二步，对不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障这五种情况下的的动力学时域响以及频域响应进行特征提取；所述特征提取为，计算时域响应的方差及均值，以及，记录频域响应的1/2X,1X,2X,3X的幅值；根据第一步中得到的特征频率，取其0.5倍旋转基频、1倍旋转基频、2倍旋转基频、3倍旋转基频，共计四个特征频率的幅值；
对步骤2中得到的各组振动时域及频域响应数据提取：频域内的故障特征频率的0.5倍频率的幅值故障特征频率的1倍频率的幅值A1x、故障特征频率的2倍频率的幅值A2x、故障特征频率的3倍频率的幅值A3x以及时域内的均值Amean和方差s2共计六个特征参数；
第三步，对得到的n1、n2、n3、n4和n5组数据的六个特征参数进行归一化处理；
所述n1、n2、n3、n4和n5组数据之和为W组；
所述归一化处理如下：
记所述n1、n2、n3、n4和n5组数据的1/2x特征频率的幅值分别为
归一化处理后所得的n1、n2、n3、n4和n5组数据为：
记n1、n2、n3、n4和n5组数据的1x特征频率的幅值分别为A1xi，i＝1,2,…W；归一化处理后所得的n1、n2、n3、n4和n5组数据为：
记n1、n2、n3、n4和n5组数据的2x特征频率的幅值分别为A2xi，i＝1,2,…W；归一化处理后所得的n1、n2、n3、n4和n5组数据为：
记n1、n2、n3、n4和n5组数据的3x特征频率的幅值分别为A3xi，i＝1,2,…W；归一化处理后所得的W组数据为：
记n1、n2、n3、n4和n5组数据的时域振动幅值的均值Ameani，i＝1,2,…W；归一化处理后所得的W组数据为：
记n1、n2、n3、n4和n5组数据的时域振动幅值的方差s2i，i＝1,2,…W；归一化处理后所得的W组数据为：
步骤4，建立初步振动状态评估模型：
利用仿真数据进行bp神经网络模型训练，得到初步振动状态评估模型；
所述利用仿真数据进行bp神经网络模型训练的具体过程是：
第一步，建立故障样本：对于步骤3中得到的不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障这五种情况下的归一化后的n1、n2、n3、n4和n5组特征参数数据，依据“不平衡故障、不对中故障、转静碰摩故障、叶片掉块故障和无故障”五种类型均分为5组；每组的3/4个样本随机抽取15个作为训练样本，剩余的1/4个作为测试样本；
建立故障样本时的打标签处理为：无故障的故障类型为正常，标记为1；不平衡类型与不对中类型均为异常，标记为2；转静碰摩类型与叶片掉块类型均为故障，标记为3；
第二步，给定损失函数表达式；
式(33)中，x为训练次数，Loss(x)为训练误差，εi为风险系数；y为故障标签所标识的值，取值范围为1，2，3；y*为模型计算结果；当y＝1时，发动机实际状态为正常；若此时的诊断结果为y＝2或y＝3,则称所述情况为虚警；当y＝2时，发动机实际状态为异常,当y＝3时，发动机实际状态为故障；当y＝2以及y＝3时，若此时的诊断结果为y＝1，则称所述情况为误诊；
εi为风险系数，振动状态评估中误诊带来的风险大于虚警情况下的风险，因此在损失函数计算中要考虑到误诊带来的风险占比要高于虚警带来的风险；
当实际发动机状态为1正常时，若bp神经网络模型计算结果为2异常或3故障，此时发动机实际状态为健康，而模型计算结果为异常或故障，此时诊断结果为虚警；
此时εi取值为1；
当实际发动机状态为2异常时，若bp神经网络模型计算的结果为1正常，此时诊断结果为误诊且发动机存在故障风险，风险系数εi取值为2；若模型计算结果为故障，此时诊断结果为误诊且发动机实际故障风险低于误诊结果状态下的故障风险，风险系数εi取值为1；
当实际发动机状态为3故障时，若bp神经网络模型计算的结果为1正常，此时诊断结果为误诊且发动机存在严重故障风险，风险系数εi取值为8；若bp神经网络模型计算的结果为2异常，此时诊断结果为误诊且发动机实际故障风险高于误诊结果状态下的故障风险，风险系数εi取值为4；
辨识的四种故障在等级V故障状态中的隶属度值；将不平衡与不对中定义为异常状态，将转静碰摩与叶片掉块定义为故障状态，据此给出风险系数εi的取值；
第三步，计算bp神经网络模型的损失，
第四步，输出模型参数，所述模型参数为bp神经网络中的权值；bp神经网络模型训练完成，所述的训练完成了的bp神经网络模型为振动状态评估模型；
至此，一种适合航空发动机机载的整机振动状态评估方法完成。
2.如权利要求书1所述适合航空发动机机载的整机振动评估方法，其特征在于，
所述等截面弹性轴段单元运动微分方程的推导过程是：
在转子系统有限元法中，通常采用考虑了剪切变形的Timoshenko梁单元对弹性轴单元进行建模；每个Timoshenko梁单元有两个节点，如前所述每个节点有4个自由度，即每个梁单元有8个自由度；梁单元的8个自由度组成的广义位移向量为：
qe＝[q1 q2 q3 q4 q5 q6 q7 q8]T (4)
式中：q1,q5为轴元素两端在Y方向的位移；q2,q6为轴元素两端在Z方向的位移；q3,q7为轴元素两端绕Y轴的转角；q4,q8为轴元素两端绕Z轴的转角；
设轴段单元长度为l，轴向距离s处截面的横向位移为(V,W,B,Γ)，位移和转角的关系可近似表示为：
3.如权利要求书2所述适合航空发动机机载的整机振动评估方法，其特征在于，使用轴元素两端节点广义坐标表示该轴元素中任意微元段的位移和转角表示为：
式中：ψ1,ψ2,ψ3,ψ4是平动位移插值函数；θ1,θ2,θ3,θ4是转角位移插值函数；
所述式(6)与(7)中，形函数矩阵为：
所述式(8)与(9)中，
横向平动位移插值函数和横向转角位移插值函数的表达式为：
式中：
其中：ξ为微元段相对位置；为剪切变形系数；As为有效抗剪面积；A为横截面面积；E为弹性模量；I为截面惯性矩；G为剪切模量；D，d，l分别为元素外径、内径、长度；ν为泊松比。
4.如权利要求1所述适合航空发动机机载的整机振动评估方法，其特征在于，将式(6)、式(7)带入到式(12)、式(13)，并沿元素全长积分得到轴段单元动能与势能表达式：
式中：
利用拉格朗日方程：
5.如权利要求1所述适合航空发动机机载的整机振动评估方法，其特征在于，将所述的n1组不平衡故障的时域及频域响应、n2组不对中故障的时域及频域响应、n3组转静碰摩故障的时域及频域响应、n4组叶片掉块故障的时域及频域响应和n5组无故障情况下发动机振动的时域及频域响应，作为bp神经网络训练数据样本。
6.如权利要求1所述适合航空发动机机载的整机振动评估方法，其特征在于，利用梯度下降法计算所述bp神经网络模型的损失时的具体过程是：
Ⅰ循环计算从i＝0直至i＝训练数据的长度：
ⅰ计算第i个训练数据的权重ω和偏差b相对于损失函数的梯度于是得到每一个训练数据的权重和偏差的梯度值；
ⅱ计算所有训练数据权重ω的梯度的总和；
ⅲ计算所有训练数据偏差b的梯度的总和；
Ⅱ更新各样本的权重值和偏差值
ⅰ使用上面第(2)、(3)步所得到的结果，计算所有样本的权重和偏差的梯度的平均值；
ⅱ通过公式34，更新每个样本的权重值和偏差值；
式中：α为学习率，ωi为第i次计算得到的权重值，bi为第i次计算得到的偏差值；L为损失函数；重复所述循环计算从i＝0直至i＝训练数据的长度和更新各样本的权重值和偏差值的过程，直至损失函数收敛不变；
所述损失函数为：由式{31}运算所得的训练误差Loss(x)值随训练次数x值的变化函数。
7.如权利要求1所述适合航空发动机机载的整机振动评估方法，其特征在于，步骤4中所述风险系数εi的取值为：
当正常识别为故障或异常时，εi的取值为1；
当异常识别为故障时，εi的取值为1；
当异常识别为正常时，εi的取值为2；
当故障识别为异常时，εi的取值为4；
当故障识别为正常时，εi的取值为8。
CN202210595102.7A 2022-05-28 2022-05-28 一种适于航空发动机机载的整机振动评估方法 Active CN115077919B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210595102.7A CN115077919B (zh) 2022-05-28 2022-05-28 一种适于航空发动机机载的整机振动评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210595102.7A CN115077919B (zh) 2022-05-28 2022-05-28 一种适于航空发动机机载的整机振动评估方法

Publications (2)

Publication Number Publication Date
CN115077919A CN115077919A (zh) 2022-09-20
CN115077919B true CN115077919B (zh) 2024-04-12

Family

ID=83248565

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210595102.7A Active CN115077919B (zh) 2022-05-28 2022-05-28 一种适于航空发动机机载的整机振动评估方法

Country Status (1)

Country Link
CN (1) CN115077919B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115683687B (zh) * 2023-01-03 2023-04-18 成都大汇物联科技有限公司 一种水电机械设备动静碰磨故障诊断方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2941049A1 (fr) * 2009-01-13 2010-07-16 Snecma Procede et systeme de surveillance de phenomenes vibratoires survenant dans un moteur a turbine a gaz d'aeronef en fonctionnement
CN103149029A (zh) * 2013-01-16 2013-06-12 南京航空航天大学 利用倒频谱识别航空发动机转静碰摩部位的方法
CN103471854A (zh) * 2013-09-26 2013-12-25 沈阳黎明航空发动机(集团)有限责任公司 一种航空发动机整机振动特性分析方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2941049A1 (fr) * 2009-01-13 2010-07-16 Snecma Procede et systeme de surveillance de phenomenes vibratoires survenant dans un moteur a turbine a gaz d'aeronef en fonctionnement
CN103149029A (zh) * 2013-01-16 2013-06-12 南京航空航天大学 利用倒频谱识别航空发动机转静碰摩部位的方法
CN103471854A (zh) * 2013-09-26 2013-12-25 沈阳黎明航空发动机(集团)有限责任公司 一种航空发动机整机振动特性分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于参数识别的航空发动机转子故障诊断与定位方法;李亚伟;荆建平;张永强;牛超阳;;噪声与振动控制;20180818(第04期);181-185 *

Also Published As

Publication number Publication date
CN115077919A (zh) 2022-09-20

Similar Documents

Publication Publication Date Title
Kwan et al. A novel approach to fault diagnostics and prognostics
CN115077919B (zh) 一种适于航空发动机机载的整机振动评估方法
Shrivastava et al. Estimation of single plane unbalance parameters of a rotor-bearing system using Kalman filtering based force estimation technique
CN109211546A (zh) 基于降噪自动编码器及增量学习的旋转机械故障诊断方法
Gradzki et al. Method of shaft crack detection based on squared gain of vibration amplitude
CN111413404A (zh) 基于叶尖定时和支持向量机原理的叶片裂纹在线测量方法
Benrahmoune et al. Detection and modeling vibrational behavior of a gas turbine based on dynamic neural networks approach
Nadai et al. Equipment failure prediction based on neural network analysis incorporating maintainers inspection findings
CN111426459A (zh) 基于叶尖定时和朴素贝叶斯算法的叶片裂纹在线测量方法
IL184970A (en) Method and system for diagnosing an aircraft from measurements performed on the aircraft
Outa et al. Prognosis and fail detection in a dynamic rotor using artificial immunological system
Li et al. Transmissibility function-based fault diagnosis methods for beam-like engineering structures: a review of theory and properties
Al-Haddad et al. Investigation of frequency-domain-based vibration signal analysis for UAV unbalance fault classification
Ugechi et al. Condition-based diagnostic approach for predicting the maintenance requirements of machinery
CN116754134A (zh) 基于试验与仿真数据融合的转子不平衡状态识别方法
Gohari et al. Unbalance rotor parameters detection based on artificial neural network
Roemer et al. Machine health monitoring and life management using finite-element-based neural networks
CN109139499A (zh) 一种变频水泵状态监测与故障诊断系统及方法
Aditiya et al. Fault diagnosis system of rotating machines using Hidden Markov Model (HMM)
Musthofa et al. Vibration analysis for the classification of damage motor PT Petrokimia Gresik using fast fourier transform and neural network
Chandrasekhar et al. Impact of material uncertainty on delamination detection in composite plate structures using modal curvatures and fuzzy logic
Han et al. Hybrid model based identification of local rubbing fault in rotor systems
Souza et al. Applying Mahalanobis-Taguchi method to detect faults in rotating machinery
Samy et al. Subsonic tests of a flush air data sensing system applied to a fixed-wing micro air vehicle
Kumhar et al. Dynamic-balance monitoring scheme for industrial fans using neural-network based machine learning approach

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant