CN113221271A - 数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法 - Google Patents

Abstract

本发明公开了一种数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法，方法中，建立三维叶片模型，使用有限元方法计算模型旋转状态下的各阶固有频率并提取该模型的广义刚度矩阵，通过叶端定时技术获取叶片实体旋转状态下的各阶固有频率。划分有限元模型网格区域，并设置初始参数。以有限元模型材料杨氏模量的刚度折减系数作为待解的损伤参数，分别估算、计算损伤参数的似然函数、组稀疏先验分布，进而计算其后验概率分布；使用EM算法估算贝叶斯概率框架隐含的超参数并最大化后验概率分布，多次迭代直至得到可使叶片模型固有频率与叶片实体固有频率接近且满足实际情况的损伤参数的数值解，实现旋转叶片裂纹损伤的定量识别。

Description

数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法

技术领域

本发明涉及机械故障诊断领域，尤其涉及一种数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法。

背景技术

旋转叶片是影响航空发动机性能和安全性的重要部件，旋转叶片在航空发动机运行过程中受高速、高温、高压等极端环境的影响，易产生裂纹损伤，进而诱发掉块、断裂等严重故障。航空发动机叶片裂纹的产生往往导致其附近区域刚度的下降及叶片模态信息(模态固有频率、模态振型等)的改变，准确识别上述模态信息，通过模态信息的变化实现对裂纹参数的定量识别，对于减少发动机运行维护成本、提高其运行安全性具有重要的意义。

本背景技术部分中公开的上述信息仅用于增强对本发明的理解，因此可能包含不构成本领域普通技术人员公知的现有技术的信息。

发明内容

针对现有技术中存在的问题，本发明的目的在于提供了一种数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法，使用有限元模型修正技术，结合组稀疏贝叶斯学习方法引入先验信息，通过最大化后验概率，从概率的角度计算有限元模型各单元刚度的下降信息，使用EM算法多次迭代计算，实现对叶片裂纹参数的定量识别。

为达到上述目的，本发明采用以下技术方案予以实现：

一种数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法包括以下步骤：

步骤S1：建立待识别的旋转叶片实体对应的无损伤三维叶片模型，使用有限元方法计算所述模型旋转状态下的各阶固有频率λ^e，测量并根据测得的无损伤旋转叶片实体的固有频率λ^m对所述模型进行修正并得到基准模型，修正准则为：

其中，D为旋转叶片的设计参数，其包括材料密度ρ、泊松比ν、各单元弹性模量E_i以及叶片形状参数，使用响应曲面设计方法构建旋转叶片的设计参数D与固有频率λ^e间的关系，计算使所述修正准则最小的设计参致，即得到基准模型；

步骤S2：划分网格区域并预设相关参数：将网格单元划分为独立的g个区域，并分别设置参数η，θ_i，γ_i，B_i的初始值，其中，参数η反映测量所述固有频率λ^m时产生的误差，所述误差服从多元高斯分布N(0，ηI)，I表示单位矩阵；参数θ_i为各单元刚度折减系数以反映因裂纹产生而导致该单元刚度下降的程度，设定第i个区域的刚度折减系数向量θⁱ服从多元高斯分布N(0，γ_iB_i)；参数γ_i用于控制组间裂纹稀疏程度；B_i用于控制组内刚度折减系数的关系，其使用托普利兹(Teoplitz)矩阵形式：

步骤S3：求解当前条件下的所述固有频率关于参数θ的灵敏度矩阵S：

其中，φ_j为旋转叶片模型第j阶质量归一化模态振型，

为广义刚度矩阵，上标T表示向量或矩阵转置；

步骤S4：基于贝叶斯理论构建损伤参数贝叶斯概率框架并求解损伤参数：

其中，Δλ＝λ^m-λ^e，θ^r表示基于本次迭代所计算的损伤参数值，c^-1为常数，损伤参数的似然函数p(Δλ|θ^r，η)和先验分布

分别为：

p(Δλ|θ^r，η)～N(Sθ^r，ηI)，

其中，∑₀为：

由此，损伤参数的后验概率密度函数服从高斯分布，

其中，均值向量μ＝∑₀S^T(ηI+S∑₀S^T)^-1Δλ，方差矩阵

使用遗传算法最大化

解得本次迭代的θ^r，或直接求取数值解θ^r←∑₀S^T(λI+S∑₀S^T)^-1Δλ，进而计算θ＝θ^r(1+θ^(-1))+θ^(-1)，其中，1表示元素全部为1的列向量，θ^(-1)表示上一次迭代计算得到的损伤参数；

步骤S5：评估迭代计算效果是否满足要求，判断标准如下：

若满足上述要求，则停止计算，输出参数θ，根据参数θ的索引及数值实现对裂纹的定位与量化，若不满足，则跳转至步骤S6继续进行迭代计算；

步骤S6：使用EM算法求解超参数

继而跳转至步骤S3开始下一轮迭代计算：

其中，

表示对应于θⁱ的∑中第i个主对角线分块矩阵，

表示对应于θⁱ的S矩阵的第i个子矩阵，

其中，

为对应于θⁱ的μ的第i个向量块，Tr表示矩阵的迹，即矩阵主对角线元素之和，

其中，

和

分别为矩阵

主对角线和次对角线元素的平均值。

所述的方法中，步骤S1包括以下步骤：

步骤S101：建立待识别旋转叶片实体形状相同的无损伤三维叶片模型，使用有限元方法计算其在工作转速范围内的各阶固有频率，

步骤S102、检查旋转叶片实体与叶端定时监测系统无故障后运行，使旋转叶片转速由0均匀上升至最高再降为0，停机后再次检查，若无故障，得到叶端定时系统输出信号，若有故障，则修正后重新测量直至无误，

步骤S103、使用离散傅里叶变换的方法处理监测系统输出信号，确定旋转叶片实体的固有频率λ^m，

步骤S104、使用响应曲面设计方法构建旋转叶片模型设计参数D与固有频率λ之间的显式函数关系λ^e＝f(D)，采用优化算法更新设计参数。

所述的方法中，步骤S2中，区域的划分将g个单元区域划分为裂纹常见区C₁和裂纹非常见区C₂，对于裂纹常见区C₁区域内对应参数γ_i的初始设置取值大于裂纹非常见区C₂区域内对应参数γ_i的初始设置取值，参数B_i中r不等于1。

所述的方法中，步骤S3中，每次迭代灵敏度矩阵S均重新计算，质量归一化模态振型φ_j、广义刚度矩阵

均考虑初始或上次迭代计算的参数θ的影响。

所述的方法中，广义刚度矩阵

与总体刚度矩阵同维度，本单元所含节点对应处数据与单元刚度矩阵相同，非本单元所含节点对应处数据置为0。

所述的方法中，步骤S4中，每次迭代Δλ均重新计算，计算固有频率λ^e考虑初始或上次迭代计算的参数θ的影响。

本发明相比于现有技术的优势及积极效果在于：

结合了有限元模型修正技术与组稀疏贝叶斯方法，能够利用叶端定时监测系统获取的旋转叶片各阶固有频率和部分先验知识，从概率学的角度精确定位并量化旋转叶片实体裂纹。相比于传统的监测方法，能够在不停机的情况下实现对旋转叶片裂纹的实时运行态监测，节约时间、降低成本，提高了航空发动机运行安全性。

附图说明

图1为本发明方法的流程示意图。

以下结合附图和实施例对本发明作进一步的解释。

具体实施方式

下面结合参照附图1更详细地对本发明的具体实施例进行阐述。

为了更好地理解，图1为一种数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法的流程示意图，如图1所示，一种数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法包括以下步骤：

建立与所研究旋转叶片实体对应的无损伤三维叶片模型，使用有限元方法获取该模型旋转状态下的各阶固有频率λ^e，根据测得的无损伤旋转叶片实体的固有频率λ^m对有限元模型进行修正并得到基准模型，修正准则为：

其中，D为旋转叶片的主要设计参数，包括材料的密度ρ、泊松比ν、各单元弹性模量E_i以及叶片形状参数等，使用响应曲面设计方法构建旋转叶片模型设计参数D与固有频率λ^e间的关系，计算合理范围内使上式最小的设计参数，即可得到基准模型；

S2：划分网格区域并预设相关参数：依据专家知识中叶片常见裂纹性态将网格单元划分为独立的g个区域，并分别设置参数η，θ_i，γ_i，B_i的初始值。其中，η反映测量系统识别固有频率λ^m时产生的误差，设定误差ε服从多元高斯分布N(0，ηI)，I表示单位矩阵；θ_i为各单元刚度折减系数，也是需要计算的损伤参数，反映因裂纹产生而导致该单元刚度下降的程度，设定第i个区域的刚度折减系数向量θⁱ服从多元高斯分布N(0，γ_iB_i)；γ_i可用于控制组间裂纹稀疏程度；B_i则控制组内刚度折减系数的关系，其形式可依据专家知识所给出的先验，无此先验的情况下可使用托普利兹(Teoplitz)矩阵形式：

S3：求解当前条件下的叶片模型固有频率关于损伤参数θ的灵敏度矩阵S：

其中，φ_j为旋转叶片有限元模型第j阶质量归一化模态振型，

为广义刚度矩阵(

需与总体刚度矩阵同维度，本单元所含节点对应处数据与单元刚度矩阵相同，非本单元所含节点对应处数据置0)，上标T表示向量或矩阵转置；

基于贝叶斯理论构建损伤参数贝叶斯概率框架并求解损伤参数：

其中，Δλ＝λ^m-λ^e，θ^r表示基于本次迭代所计算的损伤参数值，c^-1为一常数，不予考虑，损伤参数的似然函数p(Δλ|θ^r，η)和先验分布

分别为：

p(Δλ|θ^r，η)～N(Sθ^r，ηI)，

其中，∑₀为：

由此，损伤参数的后验概率密度函数仍任服从高斯分布，即

其中，均值向量μ＝∑₀S^T(ηI+S∑₀S^T)^-1Δλ，方差矩阵

使用遗传算法最大化

即可解得本次迭代的θ^r，亦可在条件允许的情况下直接求取数值解θ^r←∑₀S^T(λI+S∑₀S^T)^-1Δλ。进而计算θ＝θ^r(1+θ^(-1))+θ^(-1)，其中，1表示元素全部为1的列向量，θ^(-1)表示上一次迭代计算得到的损伤参数；

评估迭代计算效果是否满足要求，判断标准如下：

若满足上述要求，则停止计算，若不满足，则跳转至步骤S6继续进行迭代计算；

S6：使用EM算法求解超参数

继而跳转至步骤S3开始下一轮迭代计算：

其中，

表示对应于θⁱ的∑中第i个主对角线分块矩阵，

表示对应于θⁱ的S矩阵的第i个子矩阵。

其中，

为对应于θⁱ的μ的第i个向量块，Tr表示矩阵的迹，即矩阵主对角线元素之和(或矩阵所有特征值之和)。

所述的方法的优选实施方式中，其中，

和

分别为矩阵

主对角线和次对角线元素的平均值。

所述的方法的优选实施方式中，其中S1包括以下步骤：

所述的方法的优选实施方式中，S101、建立与所研究旋转叶片实体形状相同的无损伤三维叶片模型，使用有限元方法计算其在工作转速范围内的各阶固有频率，

所述的方法的优选实施方式中，S102、检查旋转叶片实体与叶端定时监测系统，无故障后运行，使旋转叶片转速由0均匀上升至最高再降为0，停机后再次检查，若无故障，即可得到叶端定时系统输出信号，若有故障，则修正后重新测量直至无误，

所述的方法的优选实施方式中，使用离散傅里叶变换的方法处理监测系统输出信号，确定旋转叶片实体的固有频率λ^m，

所述的方法的优选实施方式中，S104、使用响应曲面设计方法构建旋转叶片模型设计参数D与固有频率λ^e之间的显式函数关系λ^e＝f(D)，采用优化算法更新设计参数。

所述的方法的优选实施方式中，S2中，区域的划分可依据专家知识将g个单元区域划分为裂纹常见区C₁和裂纹非常见区C₂，对于C₁区域内对应γ_i的初始设置取较大值，C₂区域内对应γ_i的初始设置取较小值。B_i中r的设定需根据实际情况设定且≠1。

所述的方法的优洗实施方式中，S3中，需要注意的是：每次迭代灵敏度矩阵S均需重新计算，即质量归一化模态振型φ_j、广义刚度矩阵

均需考虑初始或上次迭代计算的θ的影响。

所述的方法的优选实施方式中，S4中，需要注意的是：每次迭代Δλ均需重新计算，即模型计算的固有频率λ^e需考虑初始或上次迭代计算的θ的影响。

在一个实施例中，叶端定时监测系统直接输出由非均匀欠采样信号转换后的采样信号。

基于叶端定时测量的旋转叶片裂纹组稀疏贝叶斯定量识别方法包括以下步骤：

S1中，分别通过叶端定时测量系统与有限元分析软件获取旋转叶片实体与模型的各阶固有频率，二者构造残差，进而使用基于响应面法的有限元模型修正方法更新叶片模型得到基准模型；

S2中，依据专家知识划分网格单元区域并设置初始参数；

S3中，计算固有频率关于损伤参数的灵敏度矩阵；

S4中，构建损伤参数贝叶斯概率框架，通过最大化后验概率或直接采用解析公式的方法获取损伤参数；

S5中，判断损伤参数是否满足需求，若满足则直接输出，实现裂纹损伤定位、量化；若不满足，则估算超参数(S6内容)后重新进行S3继续迭代；

所述方法中，S1中，建立旋转叶片三维模型，并对模型进行修正包括以下步骤：

S101、根据所需研究的旋转叶片形状，建立三维模型并分析。使用激光三维扫描仪对叶片实体进行扫描并在Geomagic studio中处理，获取到较为准确的旋转叶片三维模型，导入Ansys中处理并计算其在工作转速范围内的各阶固有频率λ^e。

S102、确定旋转叶片实体无损伤、叶端定时监测系统无故障，使旋转叶片转速由0均匀上升至最高再降为0，停机后再次检查，确保叶片无损伤、监测系统无故障，若任意部分不满足要求，则应更换或维修。

S103、最终得到由叶端定时系统输出的采样信号，对其作离散傅里叶变换，即可得到叶片实体的各阶固有频率λ^e。

S104、使用响应曲面设计方法构建旋转叶片模型设计参数D与固有频率λ^e之间的显式函数关系λ^e＝f(D)，设置如下修正准则：

其中，D为旋转叶片的主要设计参数，包括材料的密度ρ、泊松比ν、各单元弹性模量E_i以及叶片形状参数等，使用优化算法计算合理范围内使上式最小的设计参数，即可得到适用于后续操作的基准模型。

所述的方法中，S2中，划分网格区域并预设相关参数：

依据专家知识中叶片常见裂纹性态将S1中的有限元模型网格单元划分为相互独立的g个区域(包括裂纹常见区C₁和裂纹非常见区C₂)，并分别设置参数η，θ_i，γ_i，B_i的初始值。

其中，η反映测量系统识别固有频率λ^m时产生的误差，设定误差ε服从多元高斯分布N(0，ηI)，I表示单位矩阵；

θ_i为各单元刚度折减系数，也是需要计算的损伤参数，反映因裂纹产生而导致该单元刚度下降的程度，认为第i个单元的实际刚度矩阵为

其中K_i为由基准模型提取的单元刚度矩阵，后续操作包括计算灵敏度矩阵S及求解模型固有频率λ^e均需基于

设定第i个区域的刚度折减系数向量θⁱ服从多元高斯分布N(0，γ_iB_i)；

γ_i可用于控制组间裂纹稀疏程度，C₁区域内对应γ_i的初始设置取较大值，C₂区域内对应γ_i的初始设置取较小值；

B_i则控制组内刚度折减系数的关系，其形式可依据专家知识所给出的先验，无此先验的情况下可使用托普利兹(Teoplitz)矩阵形式：

其中，r的设定需根据实际情况设定且≠1。

所述的方法中，S3中，求解当前条件下的叶片模型固有频率关于损伤参数θ的灵敏度矩阵S，参照如下公式：

其中，φ_j为旋转叶片有限元模型第j阶质量归一化模态振型，

为广义刚度矩阵(

需与总体刚度矩阵同维度，本单元所含节点对应处数据与单元刚度矩阵相同，非本单元所含节点对应处数据置0)，上标T表示向量或矩阵转置。

所述的方法中，S4中，基于贝叶斯理论构建损伤参数贝叶斯概率框架并求解损伤参数：

其中，Δλ＝λ^m-λ^e，θ^r表示基于本次迭代所计算的损伤参数值，c^-1不予考虑，损伤参数的似然函数p(Δλ|θ^r，η)和先验分布

分别为：

p(Δλ|θ^r，η)～N(Sθ^r，ηI)，

其中，∑₀为：

由此，损伤参数的后验概率密度函数仍任服从高斯分布，即

其中，均值向量μ＝∑₀S^T(ηI+S∑₀S^T)^-1Δλ，方差矩阵

使用遗传算法最大化

即可解得本次迭代的θ^r，亦可在条件允许的情况下直接求取数值解θ^r←∑₀S^T(λI+S∑₀S^T)^-1Δλ。

进而计算θ＝θ^r(1+θ^(-1))+θ^(-1)，其中，1表示元素全部为1的列向量，θ^(-1)表示上一次迭代计算得到的损伤参数。

所述的方法中，S5中，评估迭代计算效果是否满足要求，判断标准如下：

若满足上述要求，则停止计算，输出裂纹损伤参数θ，根据θ的索引及数值即可实现对裂纹的定位与量化。若不满足，则跳转至步骤S6继续进行迭代计算，阈值的设置视具体要求而定。

所述的方法中，S6中，使用EM算法求解超参数

继而跳转至步骤S3开始下一轮迭代计算：

其中，

表示对应于θⁱ的∑中第i个主对角线分块矩阵，

表示对应于θⁱ的S矩阵的第i个子矩阵。

其中，

为对应于θⁱ的μ的第i个向量块，Tr表示矩阵的迹，即矩阵主对角线元素之和(或矩阵所有特征值之和)。

其中，

和

分别为矩阵

主对角线和次对角线元素的平均值。

虽然已经通过相关附图和具体实施方式对本发明内容进行了描述，但并不限制本领域的技术人员在本发明权利要求下做出多种形式的替换与修改，且均受到本发明的保护。

Claims (6)

1.一种数字孪生驱动的航空发动机旋转叶片裂纹定量识别方法，其特征在于，所述方法包括以下步骤：

步骤S1：建立待识别的旋转叶片实体对应的无损伤三维叶片模型，使用有限元方法计算所述模型旋转状态下的各阶固有频率λ^e，测量并根据测得的无损伤旋转叶片实体的固有频率λ^m对所述模型进行修正并得到基准模型，修正准则为：

其中，D为旋转叶片的设计参数，其包括材料密度ρ、泊松比v、各单元弹性模量E_i以及叶片形状参数，使用响应曲面设计方法构建旋转叶片的设计参数D与固有频率λ^e间的关系，计算使所述修正准则最小的设计参数，即得到基准模型；

步骤S2：划分网格区域并预设相关参数：将网格单元划分为独立的g个区域，并分别设置参数η，θ_i，γ_i，B_i的初始值，其中，参数η反映测量所述固有频率λ^m时产生的误差，所述误差服从多元高斯分布N(0，ηI)，I表示单位矩阵；参数θ_i为各单元刚度折减系数以反映因裂纹产生而导致该单元刚度下降的程度，设定第i个区域的刚度折减系数向量θⁱ服从多元高斯分布N(0，γ_iB_i)；参数γ_i用于控制组间裂纹稀疏程度；B_i用于控制组内刚度折减系数的关系，其使用托普利兹矩阵形式：

步骤S3：求解当前条件下的所述固有频率关于参数θ的灵敏度矩阵S：

其中，φ_j为旋转叶片模型第j阶质量归一化模态振型，

为广义刚度矩阵，上标T表示向量或矩阵转置；

步骤S4：基于贝叶斯理论构建损伤参数贝叶斯概率框架并求解损伤参数：

其中，Δλ＝λ^m-λ^e，θ^r表示基于本次迭代所计算的损伤参数值，c^-1为常数，损伤参数的似然函数p(Δλ|θ^r，η)和先验分布

分别为：

p(Δλ|θ^r，η)～N(Sθ^r，ηI)，

其中，∑₀为：

由此，损伤参数的后验概率密度函数服从高斯分布，

其中，均值向量μ＝∑₀S^T(ηI+S∑₀S^T)^-1Δλ，方差矩阵

使用遗传算法最大化

解得本次迭代的θ^r，或直接求取数值解θ^r←∑₀S^T(λI+S∑₀S^T)^-1Δλ，进而计算θ＝θ^r(1+θ^(-1))+θ^(-1)，其中，1表示元素全部为1的列向量，θ^(-1)表示上一次迭代计算得到的损伤参数；

步骤S5：评估迭代计算效果是否满足要求，判断标准如下：

若满足上述要求，则停止计算，输出参数θ，根据参数θ的索引及数值实现对裂纹的定位与量化，若不满足，则跳转至步骤S6继续进行迭代计算；

步骤S6：使用EM算法求解超参数

继而跳转至步骤S3开始下一轮迭代计算：

其中，

表示对应于θⁱ的∑中第i个主对角线分块矩阵，

表示对应于θⁱ的S矩阵的第i个子矩阵，

其中，

为对应于θⁱ的μ的第i个向量块，Tr表示矩阵的迹，即矩阵主对角线元素之和，

其中，

和

分别为矩阵

主对角线和次对角线元素的平均值。

2.根据权利要求1所述的方法，其特征在于，优选的，步骤S1包括以下步骤：

步骤S101：建立待识别旋转叶片实体形状相同的无损伤三维叶片模型，使用有限元方法计算其在工作转速范围内的各阶固有频率，

步骤S102、检查旋转叶片实体与叶端定时监测系统无故障后运行，使旋转叶片转速由0均匀上升至最高再降为0，停机后再次检查，若无故障，得到叶端定时系统输出信号，若有故障，则修正后重新测量直至无误，

步骤S103、使用离散傅里叶变换的方法处理监测系统输出信号，确定旋转叶片实体的固有频率λ^m，

步骤S104、使用响应曲面设计方法构建旋转叶片模型设计参数D与固有频率λ之间的显式函数关系λ^e＝f(D)，采用优化算法更新设计参数。

3.根据权利要求1所述的方法，其特征在于，步骤S2中，区域的划分将g个单元区域划分为裂纹常见区C₁和裂纹非常见区C₂，对于裂纹常见区C₁区域内对应参数γ_i的初始设置取值大于裂纹非常见区C₂区域内对应参数γ_i的初始设置取值，参数B_i中r不等于1。

4.根据权利要求1所述的方法，其特征在于，步骤S3中，每次迭代灵敏度矩阵S均重新计算，质量归一化模态振型φ_j、广义刚度矩阵

均考虑初始或上次迭代计算的参数θ的影响。

5.根据权利要求1所述的方法，其特征在于，广义刚度矩阵

与总体刚度矩阵同维度，本单元所含节点对应处数据与单元刚度矩阵相同，非本单元所含节点对应处数据置为0。

6.根据权利要求1所述的方法，其特征在于，步骤S4中，每次迭代Δλ均重新计算，计算固有频率λ^e考虑初始或上次迭代计算的参数θ的影响。

Images