CN114398736A - 基于时变模型参数的滚动轴承剩余寿命预测方法 - Google Patents

基于时变模型参数的滚动轴承剩余寿命预测方法 Download PDF

Info

Publication number
CN114398736A
CN114398736A CN202210056078.XA CN202210056078A CN114398736A CN 114398736 A CN114398736 A CN 114398736A CN 202210056078 A CN202210056078 A CN 202210056078A CN 114398736 A CN114398736 A CN 114398736A
Authority
CN
China
Prior art keywords
bearing
time
value
degradation
particle
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.)
Pending
Application number
CN202210056078.XA
Other languages
English (en)
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.)
Zhengzhou University of Light Industry
Original Assignee
Zhengzhou University of Light Industry
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 Zhengzhou University of Light Industry filed Critical Zhengzhou University of Light Industry
Priority to CN202210056078.XA priority Critical patent/CN114398736A/zh
Publication of CN114398736A publication Critical patent/CN114398736A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/04—Ageing analysis or optimisation against ageing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明提供基于时变模型参数的滚动轴承剩余寿命预测方法，包括：1、轴承全寿命周期数据变换、预处理及特征提取；2、利用峭度K监测轴承故障的出现，确定轴承剩余寿命RUL预测的起点；3、利用相对均方根值RRMS作为退化模型输入，进行轴承寿命预测；4、建立时变参数退化模型；5、轴承剩余寿命预测。本发明方法不仅减轻了基于数据驱动和物理模型RUL预测方法的局限性，而且更加准确地预测轴承的未来退化趋势，从而使该模型得到更好的预测性能，可以用于各种变工况下轴承剩余寿命预测。

Description

基于时变模型参数的滚动轴承剩余寿命预测方法
技术领域
本发明涉及基于时变模型参数的滚动轴承剩余寿命预测方法，属于机械装备智能维护技术领域。
背景技术
滚动轴承作为旋转机械的重要零部件之一，一旦发生故障需要消耗大量人力、物力进行检修。通常轴承的工况复杂，工作环境恶劣，生产过程中轴承的性能和健康状态不可避免的会发生退化现象。若能对滚动轴承运行状态实时监测与剩余寿命(RUL)预测，并对故障提前预警将事后维修转变成视情维护，可有效降低机器的运维成本并且避免重大安全事故的发生，保证设备的安全运行。
目前寿命预测的方法可以分为基于物理模型和基于数据驱动的方法。数据驱动的方法是将历史数据和条件监控数据与机器学习技术相结合，例如人工神经网络，支持向量机，高斯过程回归，建立预测模型并利用训练模型来预测寿命。数据驱动方法的优点是，它们可以从可用的传感器数据中直接了解轴承的潜在退化趋势，因此，用户不需要知道轴承的确切故障机制。另一方面基于模型的方法是通过失效机理建立数学或者物理模型来描述轴承整体衰退趋势，然后利用统计估计技术，如卡曼滤波，粒子滤波，无迹粒子滤波等对轴承状态进行更新和估计并预测轴承的RUL。基于模型的方法主要优点是预测结果更直观，并且充分结合专家知识和来自机器的实时信息，可以有效的进行轴承的RUL。
但是上述方法主要针对恒定工况下的轴承的RUL预测，对于变转速工况下滚动轴承RUL预测并不适用，难以从非平稳时域振动信号中准确提取蕴含轴承退化信息的特征，另外，建立退化模型时通常需要明确的先验知识和广泛的经验数据，退化模型没有根据轴承实际退化状态实时更新，与轴承的实际退化趋势误差大，影响轴承RUL预测的精度。
发明内容
针对现有技术的不足，本发明的目的是提供基于时变模型参数的滚动轴承剩余寿命预测方法。
为了实现上述目的，本发明所采用的技术方案是：
一种基于时变模型参数的滚动轴承剩余寿命预测方法，包括如下步骤：
步骤1、轴承全寿命周期数据变换、预处理及特征提取；
将变转速工况下轴承全寿命数据进行角域转换，将原始非平稳时域信号的转速脉冲通过数值积分转换成转角时间曲线，并进行等角度重采样，将时域信号转换成角域平稳信号，然后对角域信号进行特征提取，得到滚动轴承的退化特征指标：峭度K、均方根RMS、相对均方根值RRMS；
步骤2.利用峭度K监测轴承故障的出现，进而确定轴承剩余寿命RUL预测的起点；
步骤3.利用相对均方根值RRMS作为退化模型输入，用其进行轴承寿命预测；
步骤4.建立时变参数退化模型
利用Paris定理以及有限元分析得到轴承时变物理参数ΔK与输入裂纹长度x的数值关系，将其传递至轴承的退化模型，并建立轴承退化过程中的状态转移方程和观测方程；
步骤5.轴承剩余寿命预测
利用无迹粒子滤波算法UPF算法对输入到退化模型中的RRMS进行状态更新和估计，进而预测每个时刻轴承的退化指标RRMS；然后利用退化指标RRMS预测值和其达到失效阈值的时间，进一步映射到滚动轴承的RUL。
所述步骤1的具体计算过程如下：
在短时间间隔内，可将旋转轴的运动看作均匀角加速运动，转角θ与时间t的关系如下:
θ(t)＝a0+a1t+a2t2 (1)
式中a0、a1、a2为多项式系数；
令三个连续键相时标为t1、t2、t3，则有：
Figure BDA0003476528720000021
令相邻两个键相脉冲的角度间隔为固定值2π，根据式(2)可以求出待定系数a0、a1、a2，代入式(1)反求出与转角变化相对应的时间：
Figure BDA0003476528720000022
对转角按等角度采样间隔Δθ进行离散化处理，式(3)可变为：
Figure BDA0003476528720000031
式中，Tn为第n个等角度采样点的时刻,Δθ为采样间隔；
根据等角度采样时刻Tn，对非平稳时域进行拉格朗日插值运算，求出非平稳时域在Tn时刻的幅值，进而得重采样的角域信号。
所述步骤2的具体过程如下：
首先计算轴承正常工作状态下峭度的均值μ与标准差σk，然后将峭度区间[μ-3σk,μ+3σk]定义为轴承正常工作区间，当轴承tm时刻的峭度km超出定义的轴承的正常区间，则判断轴承进入退化阶段，将tm时刻作为轴承进入退化阶段的起点，并开始对轴承进行剩余寿命预测。
所述步骤3相对均方根值RRMS的具体定义公式如下：
Figure BDA0003476528720000032
式中XRRMS(t)为t时刻信号的RRMS值；XRMS(t)为t时刻信号的RMS值；XRMS(T)为轴承寿命预测起点的RMS值。
所述步骤4的具体步骤为：
a.建立与试验轴承一致的有限元模型，设置裂纹的长度和深度信息，模拟疲劳裂纹的扩展，并依据试验条件施加约束和载荷；
b.由于普通单元形函数不能表征裂纹前沿应力和应变的奇异性，采用8节点等参奇异单元进行网格划分和仿真，奇异单元是一种节点位于四分之一处的二次单元，其形函数为：
Figure BDA0003476528720000033
Figure BDA0003476528720000034
式中，u，v为节点坐标，u2，u3为相应节点坐标，η为常数，L为单元边长；
c.得到一组特定扩展形式下的离散点的应力强度因子幅后，通过数值插值拟合求出ΔK与裂纹长度平方根值的数学关系；
d.利用Paris定理以及有限元分析得到时变物理参数ΔK与输入的裂纹长度x的数值关系，将其传递至轴承的退化模型；
Figure BDA0003476528720000041
xk＝xk-1+C(ΔK)m xk-1Δtk (4)
式中，C和m是与材料有关的常数，Δσ为一个常数，Δtk＝tk-tk-1，xk为当前时刻系统的状态值，tk为当前时刻，tk-1为上一时刻；
结合式(4)建立轴承退化过程中的状态转移方程和观测方程：
Figure BDA0003476528720000042
yk＝xk+vt
式中，α＝C(Δσ)m,β＝m/2；C和m是与材料有关的常数，Δσ为一个常数，Δtk＝tk-tk-1，wt,vt为系统随机噪声，xk为当前时刻系统的状态值,yk为当前系统状态的观测值，k为当前时刻，k-1为上一时刻。
所述步骤5的具体方法为：
首先利用峭度判断轴承剩余寿命预测的起点，然后将退化特征指标RRMS带入时变参数退化模型作为模型的输入，利用UPF算法对轴承状态进行更新和估计，进而得到每一时刻轴承RRMS的预测值，并通过其达到失效阈值的时间映射到轴承的剩余寿命。
所述无迹粒子滤波算法UPF采用无迹卡曼滤波算法UKF来构造粒子滤波的重要性采样分布，将最新观测信息融入重要性分布函数，加强了粒子的跟踪能力，提高目标跟踪的精度，有效解决粒子退化问题，其具体步骤如下：
(1)初始化，从初始分布采样N个粒子生成原始粒子集
Figure BDA0003476528720000043
且每个粒子对应的权值系数为：
Figure BDA0003476528720000044
其中，
Figure BDA0003476528720000045
为每个粒子初始权重；
(2)重要性分布，针对每个粒子，利用UKF算法对每个粒子的状态进行估计得到
Figure BDA0003476528720000046
根据高斯分布
Figure BDA0003476528720000047
得到粒子xi：
Figure BDA0003476528720000051
Figure BDA0003476528720000052
其中，
Figure BDA0003476528720000053
为每个粒子的状态值，
Figure BDA0003476528720000054
为每个粒子的状态分布；
(3)根据当前时刻t的测量值来更新各粒子权重，并对权重进行归一化处理：
Figure BDA0003476528720000055
Figure BDA0003476528720000056
其中，
Figure BDA0003476528720000057
为当前时刻的状态值，
Figure BDA0003476528720000058
为前一时刻的状态值，yt为当前时刻的观测值，
Figure BDA0003476528720000059
为当前时刻粒子权重，
Figure BDA00034765287200000510
为前一时刻粒子权重；
(4)重采样，为了增加有效粒子的数量，需要剔除权值较小的粒子，复制权值较大的粒子。其实现过程如下：
计算有效粒子的数量：
Figure BDA00034765287200000511
其中，Neff为有效粒子数量，
Figure BDA00034765287200000512
为当前时刻粒子权重；
如果Neff≤Nthr，Nthr表示阈值，则说明粒子已经严重退化，需要根据重要性权重进行重采样，将粒子映射为等权重的N个粒子；
(5)状态估计，利用更新后的粒子及其权值估计当前时刻t的状态。
技术方案中文词语的英文缩写：
剩余寿命预测：RUL、无迹粒子滤波：UPF、峭度：K、均方根：RMS、相对均方根：RRMS。
本发明有益效果：
本发明提出的一种基于时变退化模型参数的变转速滚动轴承剩余寿命预测方法，首先利用角域变换对原始时域信号进行处理并提取轴承实时退化特征，将变工况下非平稳的时域振动信号转换为平稳的角域信号，从而更加准确地提取各种具有代表性的退化信息。
本发明首先建立与轴承性能退化状态相一致的时变参数退化模型，使其更加符合轴承退化的趋势，并利用无迹粒子滤波算法(UPF)对轴承退化状态进行估计，并对其RUL进行预测。本发明基于Paris公式建立轴承退化模型过程中，模型的参数不是根据经验公式选取，而是通过建立与轴承退化模型一致的有限元模型，利用有限元分析结果对模型参数实时更新，使其更加符合轴承真实的退化趋势，从而提高轴承RUL的预测精度。
本发明方法不仅减轻了基于数据驱动和物理模型RUL预测方法的局限性，而且更加准确地预测轴承的未来退化趋势，从而使该模型得到更好的预测性能，可以用于各种变工况下轴承剩余寿命预测。
附图说明
图1为本发明滚动轴承剩余寿命预测方法的流程图。
图2为本发明应用例5个试验轴承全寿命周期健康指数图。
其中，a.试验轴承1；b.试验轴承2；c.试验轴承3；d.试验轴承4；e.试验轴承5；HI为：健康指标，K为峭度，RMS为均方根。
图3为本发明应用例5个试验轴承退化指标RRMS预测结果图。
其中，a.试验轴承1；b.试验轴承2；c.试验轴承3；d.试验轴承4；e.试验轴承5；Reallife为：真实值；Predictded life为：预测值。
图4为本发明应用例5个试验轴承剩余寿命预测结果图。
其中，a.试验轴承1；b.试验轴承2；c.试验轴承3；d.试验轴承4；e.试验轴承5；Reallife为实际寿命；Time-varying model为时变模型预测寿命；Failure threshold为失效阈值。
具体实施方式
以下结合实施例对本发明的具体实施方式作进一步详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
一种基于时变模型参数的滚动轴承剩余寿命预测方法，其步骤如下：
步骤1.轴承全寿命周期数据变换、预处理及特征提取；
将变转速工况下轴承全寿命数据进行角域转换，将原始非平稳时域信号的转速脉冲通过数值积分转换成转角时间曲线，并进行等角度重采样，将时域信号转换成角域平稳信号，得到一个完整的角域数据，然后对轴承每个时刻平稳的角域信号进行特征提取，从而得到滚动轴承的退化特征指标(峭度、均方根等)。
具体计算过程如下：
在短时间间隔内，旋转轴作均匀角加速运动，转角θ与时间t的关系如下:
θ(t)＝a0+a1t+a2t2 (1)
式中a0、a1、a2为多项式系数。
令三个连续键相时标为t1、t2、t3，则有：
Figure BDA0003476528720000071
令相邻两个键相脉冲的角度间隔为固定值2π，根据式(2)可以求出待定系数a0、a1、a2，代入式(1)反求出与转角变化相对应的时间：
Figure BDA0003476528720000072
对转角按等角度采样间隔Δθ进行离散化处理，式(3)可变为：
Figure BDA0003476528720000073
式中，Tn为第n个等角度采样点的时刻,Δθ为采样间隔。
根据等角度采样时刻Tn，对非平稳时域信号进行拉格朗日插值运算，求出非平稳时域信号在Tn时刻的幅值，进而得重采样的角域信号。
步骤2.利用峭度(K)监测轴承故障的出现，进而确定轴承剩余寿命(RUL)预测的起点。具体过程如下：
首先计算轴承正常工作状态下峭度的均值μ与标准差σk，然后将峭度区间[μ-3σk,μ+3σk]定义为轴承正常工作区间，当轴承tm时刻的峭度km超出定义的轴承的正常区间，则判断轴承进入退化阶段，将tm时刻作为轴承进入退化阶段的起点，并开始对轴承进行剩余寿命预测。
步骤3.利用相对均方根值(RRMS)作为退化模型输入，用其进行轴承寿命预测。
其中RRMS的具体定义公式如下：
Figure BDA0003476528720000081
式中XRRMS(t)为t时刻信号的RRMS值；XRMS(t)为t时刻信号的RMS值；XRMS(T)为轴承寿命预测起点的RMS值。
步骤4.建立时变参数退化模型
利用Paris定理以及有限元分析得到时变物理参数ΔK与输入x的数值关系，将其传递至轴承的退化模型，并建立轴承退化过程中的状态转移方程和观测方程。
所述建立时变参数退化模型的具体步骤为：
a.建立与试验轴承一致的有限元模型，设置裂纹的长度和深度信息，模拟疲劳裂纹的扩展，并依据试验条件施加约束和载荷；
b.由于普通单元形函数不能表征裂纹前沿应力和应变的奇异性，采用8节点等参奇异单元进行网格划分和仿真，奇异单元是一种节点位于四分之一处的二次单元，其形函数为：
Figure BDA0003476528720000082
Figure BDA0003476528720000083
式中，u，v为节点坐标，u2，u3为相应节点坐标，η为常数，L为单元边长。
c.得到一组特定扩展形式下的离散点的应力强度因子幅后，通过数值插值拟合求出ΔK与裂纹长度平方根值的数学关系；
d.利用Paris定理以及有限元分析得到时变物理参数ΔK与输入的裂纹长度x的数值关系，将其传递至轴承的退化模型。
Figure BDA0003476528720000084
xk＝xk-1+C(ΔK)m xk-1Δtk (4)
式中，C和m是与材料有关的常数,Δσ为一个常数，Δtk＝tk-tk-1，xk为当前时刻系统的状态值，k为当前时刻，k-1为上一时刻。
结合式(4)建立轴承退化过程中的状态转移方程和观测方程：
Figure BDA0003476528720000085
yk＝xk+vt
式中，α＝C(Δσ)m,β＝m/2。C和m是与材料有关的常数,Δσ为一个常数，Δtk＝tk-tk-1，wt,vt为系统随机噪声，xk为当前时刻系统的状态值,yk为当前系统状态的观测值，k为当前时刻，k-1为上一时刻。
步骤5.轴承剩余寿命预测
利用无迹粒子滤波算法(UPF)算法对输入到退化模型中的RRMS进行状态更新和估计，进而预测每个时刻轴承的退化指标RRMS。然后利用退化指标RRMS预测值和其达到失效阈值的时间，进一步映射到滚动轴承的RUL。
所述剩余寿命预测的方法为：
首先利用峭度判断轴承剩余寿命预测的起点，然后将退化特征指标RRMS带入时变参数退化模型作为模型的输入，利用UPF算法对轴承状态进行更新和估计，进而得到每一时刻轴承RRMS的预测值，并通过其达到失效阈值的时间映射到轴承的剩余寿命。
其中，无迹粒子滤波算法(UPF)采用无迹卡曼滤波算法(UKF)来构造粒子滤波的重要性采样分布，将最新观测信息融入重要性分布函数，加强了粒子的跟踪能力，提高目标跟踪的精度，有效解决粒子退化问题，其具体步骤如下。
(1)初始化，从初始分布采样N个粒子生成原始粒子集
Figure BDA0003476528720000091
且每个粒子对应的权值系数为：
Figure BDA0003476528720000092
其中，
Figure BDA0003476528720000093
为每个粒子初始权重。
(2)重要性分布，针对每个粒子，利用UKF算法对每个粒子的状态进行估计得到
Figure BDA0003476528720000094
根据高斯分布
Figure BDA0003476528720000095
得到粒子xi：
Figure BDA0003476528720000096
Figure BDA0003476528720000097
其中，
Figure BDA0003476528720000098
为每个粒子的状态值，
Figure BDA0003476528720000099
为每个粒子的状态分布。
(3)根据当前时刻t的测量值来更新各粒子权重，并对权重进行归一化处理：
Figure BDA0003476528720000101
Figure BDA0003476528720000102
其中，
Figure BDA0003476528720000103
为当前时刻的状态值，
Figure BDA0003476528720000104
为前一时刻的状态值，yt为当前时刻的观测值，
Figure BDA0003476528720000105
为当前时刻粒子权重，
Figure BDA0003476528720000106
为前一时刻粒子权重。
(4)重采样，为了增加有效粒子的数量，需要剔除权值较小的粒子，复制权值较大的粒子。其实现过程如下：
计算有效粒子的数量：
Figure BDA0003476528720000107
其中，Neff为有效粒子数量，
Figure BDA0003476528720000108
为当前时刻粒子权重。
如果Neff≤Nthr(Nthr表示阈值)，则说明粒子已经严重退化，需要根据重要性权重进行重采样，将粒子映射为等权重的N个粒子；
(5)状态估计，利用更新后的粒子及其权值估计当前时刻t的状态。
应用例：
一种基于时变退化模型参数的变转速滚动轴承剩余寿命预测方法，其流程如图1所示。
XJTU-SY轴承数据集由中国陕西西安交通大学设计科学与基础组件研究所和中国浙江长兴苏扬科技有限公司(SY)提供。XJTU-SY轴承数据集包含15个滚动元件轴承从运行到失效的完整数据，试验中设置采样频率为25.6kHz，采样间隔为1min，每次采样时长为1.28s。本次选取工况1(转速为2500r/min,径向力为12KN)下的5个型号为LDK UER204滚动轴承进行验证，对5个试验轴承全寿命周期数据以如下步骤进行处理：
步骤1：将原始非平稳时域信号的转速脉冲通过数值积分转换成转角时间曲线，并进行等角度重采样，将时域信号转换成角域平稳信号，然后对轴承每个时刻平稳的角域信号进行特征提取，从而得到滚动轴承的退化特征指标(峭度、均方根等)。提取特征的5个试验轴承全寿命周期健康指数如图2所示。
具体计算过程如下：
在短时间间隔内，旋转轴作均匀角加速运动，转角θ与时间t的关系如下:
θ(t)＝a0+a1t+a2t2 (1)
式中a0、a1、、a2为多项式系数。
令三个连续键相时标为t1、t2、t3，则有：
Figure BDA0003476528720000111
令相邻两个键相脉冲的角度间隔为固定值2π，根据式(2)可以求出待定系数a0、a1、a2，代入式(1)反求出与转角变化相对应的时间：
Figure BDA0003476528720000112
对转角按等角度采样间隔Δθ进行离散化处理，式(3)可变为：
Figure BDA0003476528720000113
式中，Tn为第n个等角度采样点的时刻,Δθ为采样间隔。
根据等角度采样时刻Tn，对非平稳时域信号进行拉格朗日插值运算，求出非平稳时域信号在Tn时刻的幅值，进而得重采样的角域信号。
步骤2：首先计算5个轴承正常工作状态下健康指标峭度的均值μ与标准差σk，然后将峭度区间[μ-3σk,μ+3σk]定义为轴承正常工作区间，当轴承tm时刻的峭度km超出定义的轴承的正常区间，则判断轴承进入退化阶段，将tm时刻作为为轴承进入退化阶段的起点，并对开始对轴承进行剩余寿命预测。
步骤3：利用相对均方根值(RRMS)作为退化模型输入，用其进行轴承寿命预测，其相对均方根值定义如下。
Figure BDA0003476528720000114
式中，XRRMS(t)为t时刻信号的RRMS值；XRMS(t)为t时刻信号的RMS值；XRMS(T)为轴承寿命预测起点的RMS值。
步骤4.建立时变参数退化模型：
a.建立与试验轴承一致的有限元模型，其中试验轴承为单列深沟球轴承，型号为SKF6205-2RSH，分别在轴承的内圈和外圈设置长度为1mm，深度为1mm的平面张开型裂纹。并根据实验条件将轴承内圈固定，外圈施加8500N的径向载荷，设置裂纹的长度和深度信息，模拟疲劳裂纹的扩展，并依据试验条件施加约束和载荷；
b.由于普通单元形函数不能表征裂纹前沿应力和应变的奇异性，采用8节点等参奇异单元进行网格划分和仿真，奇异单元是一种节点位于四分之一处的二次单元，其形函数为：
Figure BDA0003476528720000121
Figure BDA0003476528720000122
式中，u，v为节点坐标，u2，u3为相应节点坐标，η为常数，L为单元边长。
c.得到一组特定扩展形式下的离散点的应力强度因子幅后，通过数值插值拟合求出ΔK与裂纹长度平方根值的数学关系；
d.利用Paris定理以及有限元分析得到时变物理参数ΔK与输入x的数值关系，将其传递至轴承的退化模型，并建立轴承退化过程中的状态转移方程和观测方程，具体表达式如下。
Figure BDA0003476528720000123
xk＝xk-1+C(ΔK)m xk-1Δtk (4)
式中，C和m是与材料有关的常数,Δσ为一个常数，Δtk＝tk-tk-1。
结合式(4)建立轴承退化过程中的状态转移方程和观测方程：
Figure BDA0003476528720000124
yk＝xk+vt
式中，α＝C(Δσ)m,β＝m/2。C和m是与材料有关的常数,Δσ为一个常数，Δtk＝tk-tk-1，wt,vt为系统随机噪声。
步骤5.利用UPF算法对输入到退化模型中的RRMS进行状态更新和估计，进而预测每个时刻轴承的退化指标RRMS，5个试验轴承RRMS预测结果与实际值如图3所示，无迹粒子滤波算法(UPF)采用无迹卡曼滤波(UKF)来构造粒子滤波的重要性采样分布，将最新观测信息融入重要性分布函数，加强了粒子的跟踪能力，提高目标跟踪的精度，有效解决粒子退化问题，其具体步骤如下。
(1)初始化，令t＝0，从初始分布采样N个粒子生成原始粒子集
Figure BDA0003476528720000131
且每个粒子对应的权值系数为：
Figure BDA0003476528720000132
其中，w0 i为每个粒子初始权重。
(2)重要性分布，针对每个粒子，利用UKF算法对每个粒子的状态进行估计得到
Figure BDA0003476528720000133
根据高斯分布
Figure BDA0003476528720000134
得到粒子xi：
Figure BDA0003476528720000135
Figure BDA0003476528720000136
其中，
Figure BDA0003476528720000137
为每个粒子的状态值，
Figure BDA0003476528720000138
为每个粒子的状态分布。
(3)根据当前时刻t的测量值来更新各粒子权重，并对权重进行归一化处理：
Figure BDA0003476528720000139
Figure BDA00034765287200001310
其中，
Figure BDA00034765287200001311
为当前时刻的状态值，
Figure BDA00034765287200001312
为前一时刻的状态值，yt为当前时刻的观测值，
Figure BDA00034765287200001313
为当前时刻粒子权重，
Figure BDA00034765287200001314
为前一时刻粒子权重。
(4)重采样，为了增加有效粒子的数量，需要剔除权值较小的粒子，复制权值较大的粒子。其实现过程如下：
计算有效粒子的数量：
Figure BDA00034765287200001315
其中，Neff为有效粒子数量，
Figure BDA00034765287200001316
为当前时刻粒子权重。
如果Neff≤Nthr(Nthr表示阈值),则说明粒子已经严重退化，需要根据重要性权重进行重采样，将粒子映射为等权重的N个粒子；
(5)状态估计，利用更新后的粒子及其权值估计当前时刻t的状态。
步骤6.利用退化指标RRMS预测值和其达到失效阈值的时间，进一步映射到滚动轴承的RUL，5个试验轴承的寿命预测结果(Time-varying model)，实际值(Real life)，以及失效阈值(Failure threshold)如图4所示。

Claims (7)

1.一种基于时变模型参数的滚动轴承剩余寿命预测方法，其特征在于，包括如下步骤：
步骤1、轴承全寿命周期数据变换、预处理及特征提取；
将变转速工况下轴承全寿命数据进行角域转换，将原始非平稳时域信号的转速脉冲通过数值积分转换成转角时间曲线，并进行等角度重采样，将时域信号转换成角域平稳信号，然后对角域信号进行特征提取，得到滚动轴承的退化特征指标：峭度K、均方根RMS、相对均方根值RRMS；
步骤2.利用峭度K监测轴承故障的出现，进而确定轴承剩余寿命RUL预测的起点；
步骤3.利用相对均方根值RRMS作为退化模型输入，用其进行轴承寿命预测；
步骤4.建立时变参数退化模型
利用Paris定理以及有限元分析得到轴承时变物理参数ΔK与输入裂纹长度x的数值关系，将其传递至轴承的退化模型，并建立轴承退化过程中的状态转移方程和观测方程；
步骤5.轴承剩余寿命预测
利用无迹粒子滤波算法UPF算法对输入到退化模型中的RRMS进行状态更新和估计，进而预测每个时刻轴承的退化指标RRMS；然后利用退化指标RRMS预测值和其达到失效阈值的时间，进一步映射到滚动轴承的RUL。
2.如权利要求1所述的滚动轴承剩余寿命预测方法，其特征在于，所述步骤1的具体计算过程如下：
在短时间间隔内，可将旋转轴的运动看作均匀角加速运动，转角θ与时间t的关系如下:
θ(t)＝a0+a1t+a2t2 (1)
式中a0、a1、a2为多项式系数；
令三个连续键相时标为t1、t2、t3，则有：
Figure FDA0003476528710000011
令相邻两个键相脉冲的角度间隔为固定值2π，根据式(2)可以求出待定系数a0、a1、a2，代入式(1)反求出与转角变化相对应的时间：
Figure FDA0003476528710000021
对转角按等角度采样间隔Δθ进行离散化处理，式(3)可变为：
Figure FDA0003476528710000022
式中，Tn为第n个等角度采样点的时刻,Δθ为采样间隔；
根据等角度采样时刻Tn，对非平稳时域进行拉格朗日插值运算，求出非平稳时域在Tn时刻的幅值，进而得重采样的角域信号。
3.如权利要求1所述的滚动轴承剩余寿命预测方法，其特征在于，所述步骤2的具体过程如下：
首先计算轴承正常工作状态下峭度的均值μ与标准差σk，然后将峭度区间[μ-3σk,μ+3σk]定义为轴承正常工作区间，当轴承tm时刻的峭度km超出定义的轴承的正常区间，则判断轴承进入退化阶段，将tm时刻作为轴承进入退化阶段的起点，并开始对轴承进行剩余寿命预测。
4.如权利要求1所述的滚动轴承剩余寿命预测方法，其特征在于，所述步骤3相对均方根值RRMS的具体定义公式如下：
Figure FDA0003476528710000023
式中XRRMS(t)为t时刻信号的RRMS值；XRMS(t)为t时刻信号的RMS值；XRMS(T)为轴承寿命预测起点的RMS值。
5.如权利要求1所述的滚动轴承剩余寿命预测方法，其特征在于，所述步骤4的具体步骤为：
a.建立与试验轴承一致的有限元模型，设置裂纹的长度和深度信息，模拟疲劳裂纹的扩展，并依据试验条件施加约束和载荷；
b.由于普通单元形函数不能表征裂纹前沿应力和应变的奇异性，采用8节点等参奇异单元进行网格划分和仿真，奇异单元是一种节点位于四分之一处的二次单元，其形函数为：
Figure FDA0003476528710000031
Figure FDA0003476528710000032
式中，u，v为节点坐标，u2，u3为相应节点坐标，η为常数，L为单元边长；
c.得到一组特定扩展形式下的离散点的应力强度因子幅后，通过数值插值拟合求出ΔK与裂纹长度平方根值的数学关系；
d.利用Paris定理以及有限元分析得到时变物理参数ΔK与输入的裂纹长度x的数值关系，将其传递至轴承的退化模型；
Figure FDA0003476528710000033
xk＝xk-1+C(ΔK)mxk-1Δtk (4)
式中，C和m是与材料有关的常数，Δσ为一个常数，Δtk＝tk-tk-1，xk为当前时刻系统的状态值，tk为当前时刻，tk-1为上一时刻；
结合式(4)建立轴承退化过程中的状态转移方程和观测方程：
Figure FDA0003476528710000034
yk＝xk+vt
式中，α＝C(Δσ)m,β＝m/2；C和m是与材料有关的常数，Δσ为一个常数，Δtk＝tk-tk-1，wt,vt为系统随机噪声，xk为当前时刻系统的状态值,yk为当前系统状态的观测值，k为当前时刻，k-1为上一时刻。
6.如权利要求1所述的滚动轴承剩余寿命预测方法，其特征在于，所述步骤5的具体方法为：
首先利用峭度判断轴承剩余寿命预测的起点，然后将退化特征指标RRMS带入时变参数退化模型作为模型的输入，利用UPF算法对轴承状态进行更新和估计，进而得到每一时刻轴承RRMS的预测值，并通过其达到失效阈值的时间映射到轴承的剩余寿命。
7.如权利要求6所述的滚动轴承剩余寿命预测方法，其特征在于，所述无迹粒子滤波算法UPF采用无迹卡曼滤波算法UKF来构造粒子滤波的重要性采样分布，将最新观测信息融入重要性分布函数，加强了粒子的跟踪能力，提高目标跟踪的精度，有效解决粒子退化问题，其具体步骤如下：
(1)初始化，从初始分布采样N个粒子生成原始粒子集
Figure FDA0003476528710000041
且每个粒子对应的权值系数为：
Figure FDA0003476528710000042
其中，
Figure FDA0003476528710000043
为每个粒子初始权重；
(2)重要性分布，针对每个粒子，利用UKF算法对每个粒子的状态进行估计得到
Figure FDA0003476528710000044
根据高斯分布
Figure FDA0003476528710000045
得到粒子xi：
Figure FDA0003476528710000046
Figure FDA0003476528710000047
其中，
Figure FDA0003476528710000048
为每个粒子的状态值，
Figure FDA0003476528710000049
为每个粒子的状态分布；
(3)根据当前时刻t的测量值来更新各粒子权重，并对权重进行归一化处理：
Figure FDA00034765287100000410
Figure FDA00034765287100000411
其中，
Figure FDA00034765287100000412
为当前时刻的状态值，
Figure FDA00034765287100000413
为前一时刻的状态值，yt为当前时刻的观测值，
Figure FDA00034765287100000414
为当前时刻粒子权重，
Figure FDA00034765287100000415
为前一时刻粒子权重；
(4)重采样，为了增加有效粒子的数量，需要剔除权值较小的粒子，复制权值较大的粒子。其实现过程如下：
计算有效粒子的数量：
Figure FDA00034765287100000416
其中，Neff为有效粒子数量，
Figure FDA00034765287100000417
为当前时刻粒子权重；
如果Neff≤Nthr，Nthr表示阈值，则说明粒子已经严重退化，需要根据重要性权重进行重采样，将粒子映射为等权重的N个粒子；
(5)状态估计，利用更新后的粒子及其权值估计当前时刻t的状态。
CN202210056078.XA 2022-01-18 2022-01-18 基于时变模型参数的滚动轴承剩余寿命预测方法 Pending CN114398736A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210056078.XA CN114398736A (zh) 2022-01-18 2022-01-18 基于时变模型参数的滚动轴承剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210056078.XA CN114398736A (zh) 2022-01-18 2022-01-18 基于时变模型参数的滚动轴承剩余寿命预测方法

Publications (1)

Publication Number Publication Date
CN114398736A true CN114398736A (zh) 2022-04-26

Family

ID=81231121

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210056078.XA Pending CN114398736A (zh) 2022-01-18 2022-01-18 基于时变模型参数的滚动轴承剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN114398736A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116227366A (zh) * 2023-05-08 2023-06-06 浙江大学 两阶段电机绝缘寿命预测方法
CN118013244A (zh) * 2024-04-07 2024-05-10 中核武汉核电运行技术股份有限公司 一种轴承寿命预测方法及装置

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116227366A (zh) * 2023-05-08 2023-06-06 浙江大学 两阶段电机绝缘寿命预测方法
CN116227366B (zh) * 2023-05-08 2023-08-11 浙江大学 两阶段电机绝缘寿命预测方法
CN118013244A (zh) * 2024-04-07 2024-05-10 中核武汉核电运行技术股份有限公司 一种轴承寿命预测方法及装置

Similar Documents

Publication Publication Date Title
Zhao et al. Performance prediction using high-order differential mathematical morphology gradient spectrum entropy and extreme learning machine
CN109766583B (zh) 基于无标签、不均衡、初值不确定数据的航空发动机寿命预测方法
CN114398736A (zh) 基于时变模型参数的滚动轴承剩余寿命预测方法
CN106951695B (zh) 多工况下的机械设备剩余使用寿命计算方法及系统
CN113221280B (zh) 一种基于数字孪生的滚动轴承建模与模型更新方法及系统
Yan et al. A dynamic multi-scale Markov model based methodology for remaining life prediction
CN112507452A (zh) 航空发动机涡轮叶片可靠性数字孪生建模方法
CN113092115B (zh) 数模联合驱动的全寿命滚动轴承数字孪生模型构建方法
CN103983453A (zh) 一种航空发动机的执行机构和传感器故障诊断的区分方法
Cui et al. Comprehensive remaining useful life prediction for rolling element bearings based on time-varying particle filtering
Wang et al. Online bearing fault diagnosis using numerical simulation models and machine learning classifications
CN103793601A (zh) 基于异常搜索和组合预测的汽轮机组在线故障预警方法
CN110006552B (zh) 一种机组设备温度异常检测方法
CN112347571B (zh) 考虑模型和数据不确定性的滚动轴承剩余寿命预测方法
CN110703594A (zh) 一种航空发动机多变量孪生支持向量机的健康预测方法
CN115438726A (zh) 一种基于数字孪生技术的设备寿命与故障类型预测方法及系统
CN111881568A (zh) 一种提高风功率预测精度的方法及系统
Shao et al. Fault prognosis and diagnosis of an automotive rear axle gear using a RBF-BP neural network
Valeti et al. Remaining useful life estimation of wind turbine blades under variable wind speed conditions using particle filters
CN113344275B (zh) 一种基于lstm模型的浮式平台波浪爬升在线预报方法
CN114705432A (zh) 防爆电机轴承健康状态评估方法及系统
Kang et al. A dual-experience pool deep reinforcement learning method and its application in fault diagnosis of rolling bearing with unbalanced data
Wu et al. Fault diagnosis of pitch sensor bias for wind turbine based on the multi-innovation Kalman filter
CN114936494A (zh) 数据物理融合驱动的高温部件可靠性评定方法和系统
Niu et al. A hybrid bearing prognostic method with fault diagnosis and model fusion

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