CN107577902B - 一种基于ukf的飞机疲劳结构剩余寿命预测方法 - Google Patents

一种基于ukf的飞机疲劳结构剩余寿命预测方法 Download PDF

Info

Publication number
CN107577902B
CN107577902B CN201710995498.3A CN201710995498A CN107577902B CN 107577902 B CN107577902 B CN 107577902B CN 201710995498 A CN201710995498 A CN 201710995498A CN 107577902 B CN107577902 B CN 107577902B
Authority
CN
China
Prior art keywords
state parameter
fatigue
follows
equation
crack
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
CN201710995498.3A
Other languages
English (en)
Other versions
CN107577902A (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.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
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 Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201710995498.3A priority Critical patent/CN107577902B/zh
Publication of CN107577902A publication Critical patent/CN107577902A/zh
Application granted granted Critical
Publication of CN107577902B publication Critical patent/CN107577902B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

一种基于UKF的飞机疲劳结构剩余寿命预测方法,本发明涉及基于UKF的飞机疲劳结构剩余寿命预测方法。本发明为了解决现有方法飞机疲劳结构剩余寿命低的缺点。本发明包括:步骤一:基于Paris疲劳裂纹扩展公式,建立状态空间评估模型;步骤二:对步骤一建立的状态空间评估模型利用无迹卡尔曼滤波算法进行滤波,得到准确的状态参数向量xk;步骤三:利用步骤二得到的准确的状态参数向量xk,进行结构的裂纹扩展剩余寿命预测。通过对比实验可知,本发明算法的预测结果优于EKF算法,且预测得到的RUL绝对相对误差小于10%。本发明应用于飞机疲劳结构剩余寿命预测领域。

Description

一种基于UKF的飞机疲劳结构剩余寿命预测方法
技术领域
本发明涉及基于UKF的飞机疲劳结构剩余寿命预测。
背景技术
由于多结构设计和飞行过程中承受随机动态载荷,疲劳损伤是飞机结构的主要失效方式之一。由于材料的微观结构变化能够对结构的疲劳性能造成很大的影响,而且材料的微观结构很难在制造过程中得到控制,即使采用相同的材料制造相同的飞机结构,制造出的构件也会表现出不同的疲劳性能。特别是起飞和降落过程中的增压/减压循环载荷易使机身壁板产生疲劳裂纹。据统计,因交变载荷造成的疲劳断裂事故占到机械结构失效总数的95%(高镇同,熊峻江.疲劳可靠性[M].北京航空航天大学出版社,2000)。而且在达到疲劳寿命时,结构在断裂失效前没有变形等明显先兆,具有极大的危险性。
事实上,更多的疲劳裂纹并未造成事故,但却由于需要对飞机结构进行维修严重地影响了正常的飞行训练,削弱了部队战斗力(刘芳,赵建印,宋贵宝.任务准备阶段机群战备完好率评估模型[J].哈尔滨工业大学学报,2008,40(3):488-491)。因为结构件查出疲劳裂纹而使大批飞机停飞的事时有发生,大大降低了飞机的战备完好率,影响装备可靠性。为避免疲劳裂纹可能产生的突然失效而需要对飞机结构进行复杂的检查、探伤、维修,也会使得飞机的维护费用、日常运行费用大大增加,导致飞机的经济性变差(WANG Y W,GOGU C,BINAUD N,et al.A cost driven predictive maintenance policy for structuralairframe maintenance[J].中国航空学报(英文版),2017,30(3):1242-1257.)。随着传感器技术的不断发展,在飞机的设计上越来越多的传感器被安装来采集关键部位的疲劳信息(白生宝,肖迎春,刘马宝,等.智能涂层传感器监测裂纹的工程适用性[J].无损检测,2015,37(1):42-44.PARK C,PETERS K.Optimization of embedded sensor placement forstructural health monitoring of composite airframes[J].AIAA Journal,2012,50(11):2536-2545.STASZEWSKI WJ.Fatigue crack detection using smart sensortechnologies[J].Fatigue&Fracture of Engineering Materials&Structures,2010,31(8):609-610.),支持更为科学的视情维修。视情维修(周炳海,陶红玉,綦法群.带随机突变的两阶段退化系统视情维修建模[J].哈尔滨工业大学学报,2016,48(1):87-93.UN J Z,ZUO H F,WANG W B,et al.Prognostics uncertainty reduction by fusing on-linemonitoring data based on a state-space-based degradation model[J].MechanicalSystems and Signal Processing,2014,45:396-407.XIA T B,XI L F,ZHOU X J,etal.Condition-based maintenance for intelligent monitored series system withindependent machine failure modes.International Journal of ProductionResearch,2013,51(15):4585-4596.)由于能够充分利用采集的实时状态信息来减少不必要的定期停机维修,被广泛的应用于飞机维修计划的制定。由于视情维修计划的制定需要参照待维修结构的预期失效时间,因此结构的剩余寿命(remaining useful life,RUL)预测对于基于视情维修的计划制定十分的重要。
疲劳结构寿命的预测建立在对裂纹扩展趋势的估计基础上。疲劳裂纹扩展规律Paris公式被广泛的应用于估计金属结构的疲劳裂纹扩展趋势。Paris公式中的两个疲劳性能参数(C和m)名义值一般由结构的疲劳试验确定。由于材料的微观结构变化能够对结构的疲劳性能造成很大的影响,而且材料的微观结构很难在制造过程中得到控制,即使采用相同的材料制造相同的飞机结构,制造出的构件也会表现出不同的疲劳性能。考虑到一个机队包含多架飞机,一架飞机包含数百个结构形式相同的飞机壁板,疲劳性能参数(C和m)的随机性会影响机队每架飞机每个壁板的疲劳裂纹扩展,最终导致壁板之间呈现出不同的疲劳失效时间。由于飞机壁板疲劳失效时间的确定对于机队飞行任务和维修计划的制定十分重要,因此壁板疲劳性能参数的分散性和疲劳失效时间的确定需要得到充分的考虑。
发明内容
本发明的目的是为了解决现有方法飞机疲劳结构剩余寿命低的缺点,而提出一种基于UKF的飞机疲劳结构剩余寿命预测方法。
一种基于UKF的飞机疲劳结构剩余寿命预测方法包括以下步骤:
步骤一:基于Paris(帕里斯)疲劳裂纹扩展公式,建立状态空间评估模型:
(1)系统状态参数转移模型为:
Figure BDA0001442415040000021
其中xk为状态参数向量,ak为结构裂纹,Ck和mk为结构的材料性能参数,g(·)为裂纹扩展方程,wa,k-1、wC,k-1、wm,k-1分别为ak、Ck、mk的系统过程噪声,f(·)非线性系统状态转移方程,Wk-1为系统噪声向量,k为时间或载荷周期;
(2)系统状态参数观测方程为:
Figure BDA0001442415040000022
其中zk为观测向量,za,k为结构裂纹的观测值,zgrow_a,k为结构裂纹增量的观测值,vk为结构裂纹的观测噪声,h'(xk)为观测方程,Vk为系统的观测噪声向量;
由于方程中状态参数a,C和m存在分散性,因此需要对其进行评估(滤波)。
步骤二:对步骤一建立的状态空间评估模型利用无迹卡尔曼滤波算法进行滤波,得到准确的状态参数向量xk
步骤三:利用步骤二得到的准确的状态参数向量xk,进行结构的裂纹扩展剩余寿命预测。
飞机疲劳结构为壁板、耳片、机翼的盒段、发动机的叶片、起落架等。
本发明的有益效果为:
针对飞机壁板疲劳性能参数的分散性问题,本发明基于Paris疲劳裂纹扩展公式建立了离散的状态参数评估模型,并应用无迹卡尔曼滤波算法(UKF)对疲劳性能参数和裂纹的扩展趋势进行了估计。最后依据估计的裂纹扩展趋势,对飞机壁板的剩余寿命(RUL)进行预测。
通过对比实验可知,UKF算法(本发明算法)的预测结果优于EKF算法,且预测得到的RUL绝对相对误差小于10%。由此可见UKF算法能够很好的处理飞机结构疲劳裂纹剩余寿命预测问题。
附图说明
图1为机身壁板图;
图2为预制裂纹图;
图3为5个试样的裂纹扩展曲线图;
图4为试件1的疲劳性能参数拟合结果图;
图5为试件2的疲劳性能参数拟合结果图;
图6为试件3的疲劳性能参数拟合结果图;
图7为试件4的疲劳性能参数拟合结果图;
图8为试件5的疲劳性能参数拟合结果图;
图9为试件1的仿真结果图;
图10为试件2的仿真结果图;
图11为试件3的仿真结果图;
图12为试件4的仿真结果图;
图13为试件5的仿真结果图。
具体实施方式
具体实施方式一:一种基于UKF的飞机疲劳结构剩余寿命预测方法包括以下步骤:
步骤一:基于Paris(帕里斯)疲劳裂纹扩展公式,建立状态空间评估模型:
(1)系统状态参数转移模型为:
Figure BDA0001442415040000041
其中xk为状态参数向量,ak为结构裂纹,Ck和mk为结构的材料性能参数,g(·)为裂纹扩展方程,wa,k-1、wC,k-1、wm,k-1分别为ak、Ck、mk的系统过程噪声,f(·)非线性系统状态转移方程,Wk-1为系统噪声向量,k为时间或载荷周期;所述系统状态参数为结构裂纹和材料性能参数;
(2)系统状态参数观测方程为:
Figure BDA0001442415040000042
其中zk为观测向量,za,k为结构裂纹的观测值,zgrow_a,k为结构裂纹增量的观测值,vk为结构裂纹的观测噪声,h'(xk)为观测方程,Vk为系统的观测噪声向量;
由于方程中状态参数a,C和m存在分散性,因此需要对其进行评估(滤波)。
步骤二:对步骤一建立的状态空间评估模型利用无迹卡尔曼滤波算法进行滤波,得到准确的状态参数向量xk
步骤三:利用步骤二得到的准确的状态参数向量xk,进行结构的裂纹扩展剩余寿命预测。
具体实施方式二:本实施方式与具体实施方式一不同的是:所述步骤一中建立状态空间评估模型的具体过程为:
由于工况以及裂纹位置的不同,许多学者提出了针对其特定问题的疲劳裂纹扩展模型。通过对不同种类的军用飞机进行机身结构疲劳试验,Molent L等[4]得出了Paris疲劳裂纹扩展公式能够充分处理飞机典型结构的疲劳裂纹扩展问题。因此,选取Paris疲劳裂纹扩展公式来建立适应于UKF算法的疲劳性能状态参数评估模型。传统的Paris公式如下:
Figure BDA0001442415040000051
式中a表示裂纹长度,N表示应力循环次数,da/dN表示裂纹扩展速率,C和m是材料系数,即疲劳性能参数;ΔK表示应力强度因子幅,ΔK与压差p,机身半径r和壁板厚度t,存在如下关系:
Figure BDA0001442415040000052
由式(1)可知,裂纹的扩展是一个连续的累积过程。随着测试技术和传感器技术(如无损检测和智能涂层)的发展,机身壁板在单位循环周期结束的裂纹长度可以通过传感器得到,因此裂纹的扩展可以转化为一个离散的累积过程。通过Euler(欧拉)方法,取dN=1,式(1)转换为如下的离散递归形式:
Figure BDA0001442415040000053
其中pk-1为第k-1个时刻或者循环周期的压差,ak-1为第k-1个时刻或者循环周期的裂纹长度;
由于压差在不同的飞行周期中存在波动,因此将压差作为一个随机变量表示如下:
Figure BDA0001442415040000054
式中Δpk表示均值压差
Figure BDA0001442415040000055
的波动,且其服从
Figure BDA0001442415040000056
Figure BDA0001442415040000057
为Δp的方差,pk的差异用来反映不同飞行周期中飞机巡航高度的变化。根据式(4),式(3)进一步表达为:
Figure BDA0001442415040000058
由于压差的变化通常十分小,以
Figure BDA00014424150400000516
中心,对式(5)进行一阶泰勒展开,得到:
Figure BDA0001442415040000059
式中,
Figure BDA00014424150400000510
为一阶偏导,表示如下:
Figure BDA00014424150400000511
Figure BDA00014424150400000512
当做系统过程噪声,
Figure BDA00014424150400000513
为给定的常量,式(6)进一步表示如下:
Figure BDA00014424150400000514
Figure BDA00014424150400000515
式中wa,k-1为系统过程噪声,且服从wa,k-1~N(0,Qa,k),Qa,k为wa,k-1的方差,Qa,k表示如下:
Figure BDA0001442415040000061
裂纹长度是通过传感器测量得到的,由于测量环境和传感器自身误差等因素,测量过程中难免引入误差,裂纹长度测量模型表示如下:
zk=h(ak)+vk (11)
式中h表示测量函数,取为恒等函数;vk表示裂纹长度测量误差,且其服从vk~N(0,Rk),Rk为vk的方差;
由式(3)可知,疲劳性能参数C和m的分散性,会造成裂纹长度呈现出一定的分散性。因此将裂纹长度a作为另一个需要评估的状态参数,则待评估的状态参数向量表示为x=[a C m]T,式(8)进一步转换为系统状态参数转移模型;
wC,k-1服从wC,k-1~N(0,QC,k-1),QC,k-1为wC,k-1的方差,wm,k-1服从wm,k-1~N(0,Qm,k-1),QC,k-1为wC,k-1的方差,QC,k-1和Qm,k-1分别由材料的疲劳试验确定。系统过程噪声方差以矩阵的形式表示如下:
Figure BDA0001442415040000062
由式(3)可知,疲劳性能参数C和m对裂纹扩展增量有很大的影响,因此选取从第k-1到第k步的裂纹扩展增量ak-ak-1作为附加的测量变量,由式(11)得到系统的观测方程;
由于裂纹的增量计算误差可以被裂纹长度测量误差所包含,因此取裂纹扩展增量计算误差为0,系统的测量误差以矩阵的形式表示为:
Figure BDA0001442415040000063
其中Rk-1为裂纹长度的观测误差方差。
其它步骤及参数与具体实施方式一相同。
具体实施方式三:本实施方式与具体实施方式一或二不同的是:所述步骤二中对步骤一建立的状态空间评估模型利用无迹卡尔曼滤波算法进行滤波,得到准确的状态参数向量xk的具体过称为:
基于逼近非线性系统模型的概率密度函数比获得非线性系统的解析模型要简单的思想,通过无迹变换,UKF算法利用一组确定性采样的Sigma点来逼近系统状态参数分布,即将Sigma点带入非线性模型中,得到对应的非线性模型数值点集,通过得到的点集求得变换后的均值和方差。假设n维状态向量x的均值和方差分别为:
Figure BDA00014424150400000711
和P,根据无迹变换,在
Figure BDA00014424150400000712
的周围可以构造一组Sigma点{χ02,...,χ2n}和对应的权重{ω02,...,ω2n}如下:
Figure BDA0001442415040000071
Figure BDA0001442415040000072
其中,λ=α2(n+κ)-n;
Figure BDA0001442415040000073
Figure BDA0001442415040000074
分别表示Sigma点χi的均值权重和方差权重;参数α确定了
Figure BDA0001442415040000075
周围Sigma点集的分布,用于调节采样得到的Sigma点和
Figure BDA0001442415040000076
的距离,且0≤α≤1;参数β包含了先验分布的高阶矩信息,对于正态分布β=2为最优值;参数κ为第二尺度参数,用于保证(n+κ)Px为半正定矩阵,当n≥3时κ=0,当n<3时κ=n-3。
对于上述的疲劳性能状态参数评估问题,在第k-1步的Sigma点χi的产生过程具体描述如下。首先,矩阵
Figure BDA0001442415040000077
可以表示为:
Figure BDA0001442415040000078
由式(16)和(17),产生的Sigma点如下:
Figure BDA0001442415040000079
步骤二一:进行初始化:k=0;对状态参数向量进行初始化,设置状态参数均值和方差分别为:
Figure BDA00014424150400000710
和P0
步骤二二:计算Sigma点集和对应的权重,根据式(16),(17),(18)和(19)得到Sigma点集{χ0,k-11,k-1,...,χ6,k-1}及对应的权重{ω0,k-11,k-1,...,ω6,k-1};
步骤二三:进行时间预测;
步骤二三一:利用非线性系统转移状态方程传递Sigma点:根据系统状态参数转移模型,得到转移后的点集为:χ'i,k-1=f(χi,k-1)(i=0,...,6);
步骤二三二:进行状态参数的预测:
Figure BDA0001442415040000081
xk|k-1为预测的状态参数;所述状态参数为ak、Ck和mk
步骤二三三:进行状态参数方差矩阵计算:
Figure BDA0001442415040000082
步骤二三四:进行sigma采样点观测值计算:γi,k-1=h(χi',k-1)(i=0,...,6);
步骤二三五:进行测量值预测:
Figure BDA0001442415040000083
步骤二三六:进行测量值方差计算:
Figure BDA0001442415040000084
步骤二四:进行测量更新;
步骤二四一:系统状态参数与测量值的协方差计算:
Figure BDA0001442415040000085
步骤二四二:进行卡尔曼增益矩阵计算:
Figure BDA0001442415040000086
步骤二四三:进行状态参数向量和方差更新:xk=xk|k-1+K(zk-zk|k-1),其中zk为在第k步对应时刻由传感器获得的裂纹长度;Px,k=Px,k|k-1-KkPz,kKk T;Kk T是卡尔曼增益矩阵Kk的转置;
步骤二五:取k=1,重新执行步骤二二至步骤二四,直至k的取值达到设定阈值为止,得到准确的状态参数向量xk
其它步骤及参数与具体实施方式一或二相同。
具体实施方式四:本实施方式与具体实施方式一至三之一不同的是:所述步骤三中利用步骤二得到的准确的状态参数向量xk,进行结构的裂纹扩展剩余寿命预测的具体过程为:
利用步骤二得到的准确的状态参数向量xk=[ak,Ck,mk]T(第k步对应时刻的疲劳性能参数的评估值xk=[ak,Ck,mk]T),根据式(3)建立的疲劳裂纹扩展离散递归模型,得到k+l时刻的疲劳裂纹长度ak+l,l为裂纹长度扩展到大于临界裂纹的最小剩余步骤对应的时刻,l>0;
Figure BDA0001442415040000091
满足式(15)的最小l定义为k时刻预测的疲劳裂纹剩余寿命(RUL):
ak+l≥ac (15)
式中ac表示结构临界裂纹长度。
其它步骤及参数与具体实施方式一至三之一相同。
采用以下实施例验证本发明的有益效果:
实施例一:
(a)飞机结构疲劳试验
选取飞机的一个典型结构—机身壁板来验证提出的方法,在相同的试验环境下对5个相同的壁板试样进行了疲劳试验。如图1所示为飞机的壁板,壁板的两个位置装有疲劳裂纹监测传感器。
试样的材料为在飞机制造中广泛应用的2024-T351铝合金。试验设备为Inston8801电液伺服疲劳试验机。试验采用试验机配置的da/dN软件记录裂纹长度及载荷循环次数。试验采用正弦波对试件进行加载,载荷频率为290Hz,最大载荷为4.5kN,应力比为0.2。如图2所示,为了控制裂纹的扩展方向,在紧固空的一侧预制了一个3mm的初始裂纹。如图3所示为5个试样的疲劳裂纹扩展曲线,图4-图8为5个试样的疲劳性能参数拟合结果。
(b)疲劳性能参数随机性分析
对式(1)两边取对数,得到
Figure BDA0001442415040000092
显然,log(da/dN)和logΔK之间呈线性关系。通过将试验所得的5个试样的a~N曲线采用7点递增多项式法进行数据处理,并将得到的数据运用最小二乘法进行线性拟合,可求得式(21)中的疲劳性能参数C和m。图4-图8分别为5个试样的疲劳性能参数拟合结果,具体的C和m值如表1所示。
表1 5个试样的疲劳性能参数C和m
Figure BDA0001442415040000101
由表1可知,在相同的试验条件下,5个试样的两个疲劳性能参数表现出不同程度的分散性。现有的大量试验数据统计研究发现(WU W F,Ni C C.Probabilistic models offatigue crack propagation and their experimental verification[J].Probabilistic Engineering Mechanics,2004,19(3):247-257.),logC和m服从正态分布,为UKF算法的应用提供了条件。通过对5个试样的数据进行统计分析,可以得到logC~N(-32.6893,0.67642)和m~N(3.9454,0.14022)。
(c)基于UKF的剩余寿命预测结果分析
用于UKF算法仿真实验的参数设置如表2所示。依据对疲劳性能参数数据的统计分析,分别选取3个状态参数的分布均值作为其初始状态。
表2仿真实验参数设置
Figure BDA0001442415040000102
基于表2的实验参数设置,选取5个试件的前30000个载荷循环的疲劳裂纹数据,利用UKF算法对疲劳性能状态参数进行估计,根据估计得到的状态参数结果并对未来的10000个载荷循环下的结构裂纹长度进行预测。为了更好的分析UKF算法在处理问题上的收敛性,使得到的仿真结果更可信,针对5个试件分别进行了30次仿真,并选取EKF算法进行了对比实验分析。如图9-图13所示,为5个试件的仿真结果。UKF算法为本发明算法。
分别利用30次仿真的均值及其与真实值的绝对相对误差来评估UKF算法的收敛性。均值和绝对相对误差的计算如式(21)和(22)。
Figure BDA0001442415040000111
Figure BDA0001442415040000112
其中,
Figure BDA0001442415040000113
表示在第i次仿真实验,第k个载荷循环时利用算法得到的裂纹长度估计值或者预测值。ak表示在第k个载荷循环时试件的真实裂纹长度。
Figure BDA0001442415040000114
为在第k个载荷循环时,30次仿真实验获得的均值。AREk为在第k个载荷循环时,30次仿真结果的绝对相对误差。
为了分析30次仿真实验结果的分散性,分别计算了30次仿真的结果与真实值的方差,以及30次仿真的结果与其均值的方差,计算公式如下:
Figure BDA0001442415040000115
Figure BDA0001442415040000116
其中,Vk为30次仿真的结果与真实值的方差。
Figure BDA0001442415040000117
为30次仿真的结果与其均值的方差。
分别选取试件1和2在不同载荷循环下(k=5000,10000,15000,20000,25000,30000,32000,34000,36000,38000,40000)得到的4个评估参数(
Figure BDA0001442415040000118
AREk,Vk
Figure BDA0001442415040000119
)来分析UKF和EKF算法在飞机结构疲劳裂纹预测问题上的表现。
表3试件1和2的4个评估参数计算结果
Figure BDA00014424150400001110
Figure BDA0001442415040000121
根据
Figure BDA0001442415040000122
和AREk反映了算法在处理问题上的收敛性,由表3可知,UKF和EKF算法在参数评估阶段都表现出了很好的收敛性,在预测阶段两个算法的收敛性都有所下降,但绝对相对误差都小于10%,且UKF算法在收敛性上整体都优于EKF算法。根据Vk反映了30次仿真的实验值与真实值的距离,由表3可知,UKF算法在参数评估和预测阶段的分散性都十分的小,且UKF算法在分散性的表现上整体都优于EKF算法。根据
Figure BDA0001442415040000123
反映了30次仿真的实验值与其均值的距离,由表3可知,UKF算法在参数评估和预测阶段的分散性都十分的小,而EKF算法在预测阶段分散性显著增大。由以上分析可知,UKF算法可以较为准确的预测出结构疲劳裂纹的真实尺寸。基于以上的分析结果,利用UKF和EKF算法对5个试件的剩余寿命进行了30次仿真预测,结果如表4所示。
表4 5个试件的剩余寿命预测结果
Figure BDA0001442415040000131
由表4可知,UKF算法在试件1,2,3和5上的预测结果都要优于EKF算法,且预测得到的RUL绝对相对误差都小于10%。由此可见UKF算法能够很好的处理飞机结构疲劳裂纹剩余寿命预测问题。
本发明还可有其它多种实施例,在不背离本发明精神及其实质的情况下,本领域技术人员当可根据本发明作出各种相应的改变和变形,但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (3)

1.一种基于UKF的飞机疲劳结构剩余寿命预测方法,其特征在于:所述基于UKF的飞机疲劳结构剩余寿命预测方法的具体过程为:
步骤一:基于Paris疲劳裂纹扩展公式,建立状态空间评估模型:
(1)系统状态参数转移模型为:
Figure FDA0002611266830000011
其中xk为状态参数向量,ak为结构裂纹,Ck和mk为结构的材料性能参数,g(·)为裂纹扩展方程,wa,k-1、wC,k-1、wm,k-1分别为ak、Ck、mk的系统过程噪声,f(·)非线性系统状态转移方程,Wk-1为系统噪声向量,k为时间或载荷周期;
(2)系统状态参数观测方程为:
Figure FDA0002611266830000012
其中zk为观测向量,za,k为结构裂纹的观测值,zgrow_a,k为结构裂纹增量的观测值,vk为结构裂纹的观测噪声,h'(xk)为观测方程,Vk为系统的观测噪声向量;
步骤二:对步骤一建立的状态空间评估模型利用无迹卡尔曼滤波算法进行滤波,得到准确的状态参数向量xk;具体过程为:
步骤二一:进行初始化:k=0;对状态参数向量进行初始化,设置状态参数均值和方差分别为:
Figure FDA0002611266830000013
和P0
步骤二二:计算Sigma点集和对应的权重,得到Sigma点集{χ0,k-11,k-1,...,χ6,k-1}及对应的权重{ω0,k-11,k-1,...,ω6,k-1};
步骤二三:进行时间预测;
步骤二三一:利用非线性系统转移状态方程传递Sigma点:根据系统状态参数转移模型,得到转移后的点集为:χ′i,k-1=f(χi,k-1)(i=0,...,6);
步骤二三二:进行状态参数的预测:
Figure FDA0002611266830000014
xk|k-1为预测的状态参数;所述状态参数为ak、Ck和mk
步骤二三三:进行状态参数方差矩阵计算:
Figure FDA0002611266830000021
步骤二三四:进行sigma采样点观测值计算:γi,k-1=h(χ′i,k-1)(i=0,...,6);
步骤二三五:进行测量值预测:
Figure FDA0002611266830000022
步骤二三六:进行测量值方差计算:
Figure FDA0002611266830000023
步骤二四:进行测量更新;
步骤二四一:系统状态参数与测量值的协方差计算:
Figure FDA0002611266830000024
步骤二四二:进行卡尔曼增益矩阵计算:
Figure FDA0002611266830000025
步骤二四三:进行状态参数向量和方差更新:xk=xk|k-1+K(zk-zk|k-1),其中zk为在第k步对应时刻由传感器获得的裂纹长度;Px,k=Px,k|k-1-KkPz,kKk T;Kk T是卡尔曼增益矩阵Kk的转置;
步骤二五:取k=1,重新执行步骤二二至步骤二四,直至k的取值达到设定阈值为止,得到准确的状态参数向量xk
步骤三:利用步骤二得到的准确的状态参数向量xk,进行结构的裂纹扩展剩余寿命预测。
2.根据权利要求1所述的一种基于UKF的飞机疲劳结构剩余寿命预测方法,其特征在于:所述步骤一中建立状态空间评估模型的具体过程为:
Paris公式如下:
Figure FDA0002611266830000026
式中a表示裂纹长度,N表示应力循环次数,da/dN表示裂纹扩展速率,C和m是材料系数;△K表示应力强度因子幅,△K与压差p,机身半径r和壁板厚度t,存在如下关系:
Figure FDA0002611266830000031
通过Eule方法,取dN=1,式(1)转换为如下的离散递归形式:
Figure FDA0002611266830000032
其中pk-1为第k-1个时刻或者循环周期的压差,ak-1为第k-1个时刻或者循环周期的裂纹长度;
将压差作为一个随机变量表示如下:
Figure FDA0002611266830000033
式中△pk表示均值压差
Figure FDA0002611266830000034
的波动,且其服从
Figure FDA0002611266830000035
Figure FDA0002611266830000036
为△p的方差,根据式(4),式(3)进一步表达为:
Figure FDA0002611266830000037
Figure FDA0002611266830000038
中心,对式(5)进行一阶泰勒展开,得到:
Figure FDA0002611266830000039
式中,
Figure FDA00026112668300000310
为一阶偏导,表示如下:
Figure FDA00026112668300000311
Figure FDA00026112668300000312
当做系统过程噪声,式(6)进一步表示如下:
Figure FDA00026112668300000313
Figure FDA00026112668300000314
式中wa,k-1服从wa,k-1~N(0,Qa,k),Qa,k为wa,k-1的方差,Qa,k表示如下:
Figure FDA00026112668300000315
裂纹长度测量模型表示如下:
zk=h(ak)+vk (11)
式中h表示测量函数,取为恒等函数;vk表示裂纹长度测量误差,且其服从vk~N(0,Rk),Rk为vk的方差;
将裂纹长度a作为另一个需要评估的状态参数,则待评估的状态参数向量表示为x=[aC m]T,式(8)进一步转换为系统状态参数转移模型;
wC,k-1服从wC,k-1~N(0,QC,k-1),QC,k-1为wC,k-1的方差,wm,k-1服从wm,k-1~N(0,Qm,k-1),QC,k-1为wC,k-1的方差,系统过程噪声方差以矩阵的形式表示如下:
Figure FDA0002611266830000041
选取从第k-1到第k步的裂纹扩展增量ak-ak-1作为附加的测量变量,由式(11)得到系统的观测方程;
取裂纹扩展增量计算误差为0,系统的测量误差以矩阵的形式表示为:
Figure FDA0002611266830000042
其中Rk-1为裂纹长度的观测误差方差。
3.根据权利要求2所述的一种基于UKF的飞机疲劳结构剩余寿命预测方法,其特征在于:所述步骤三中利用步骤二得到的准确的状态参数向量xk,进行结构的裂纹扩展剩余寿命预测的具体过程为:
利用步骤二得到的准确的状态参数向量,根据式(3)建立的疲劳裂纹扩展离散递归模型,得到k+l时刻的疲劳裂纹长度ak+l,l为裂纹长度扩展到大于临界裂纹的最小剩余步骤对应的时刻,l>0;
Figure FDA0002611266830000043
满足式(15)的最小l定义为k时刻预测的疲劳裂纹剩余寿命:
ak+l≥ac (15)
式中ac表示结构临界裂纹长度。
CN201710995498.3A 2017-10-23 2017-10-23 一种基于ukf的飞机疲劳结构剩余寿命预测方法 Active CN107577902B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710995498.3A CN107577902B (zh) 2017-10-23 2017-10-23 一种基于ukf的飞机疲劳结构剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710995498.3A CN107577902B (zh) 2017-10-23 2017-10-23 一种基于ukf的飞机疲劳结构剩余寿命预测方法

Publications (2)

Publication Number Publication Date
CN107577902A CN107577902A (zh) 2018-01-12
CN107577902B true CN107577902B (zh) 2020-11-13

Family

ID=61036897

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710995498.3A Active CN107577902B (zh) 2017-10-23 2017-10-23 一种基于ukf的飞机疲劳结构剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN107577902B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112380705B (zh) * 2020-11-16 2024-02-06 江南大学 基于非线性预测滤波算法的金属疲劳裂纹扩展预测方法
CN112651599A (zh) * 2020-12-03 2021-04-13 中国航空工业集团公司沈阳飞机设计研究所 一种飞机结构维修成本评估方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103870662A (zh) * 2014-04-01 2014-06-18 青岛科技大学 一种预测储罐剩余寿命的方法
CN105740625A (zh) * 2016-01-31 2016-07-06 太原科技大学 一种齿轮的实时剩余寿命预测方法
CN106886663A (zh) * 2017-03-29 2017-06-23 北京理工大学 齿轮弯曲疲劳寿命预测方法及装置

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160246287A1 (en) * 2014-03-13 2016-08-25 Rolls-Royce Corporation Probabilistic evaluation of turbomachinery design to predict high cycle fatigue failure

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103870662A (zh) * 2014-04-01 2014-06-18 青岛科技大学 一种预测储罐剩余寿命的方法
CN105740625A (zh) * 2016-01-31 2016-07-06 太原科技大学 一种齿轮的实时剩余寿命预测方法
CN106886663A (zh) * 2017-03-29 2017-06-23 北京理工大学 齿轮弯曲疲劳寿命预测方法及装置

Also Published As

Publication number Publication date
CN107577902A (zh) 2018-01-12

Similar Documents

Publication Publication Date Title
Coppe et al. Uncertainty reduction of damage growth properties using structural health monitoring
CN107133400B (zh) 一种飞机结构疲劳可靠度贝叶斯组合预测方法
Wang et al. Determination of Paris' law constants and crack length evolution via Extended and Unscented Kalman filter: An application to aircraft fuselage panels
CN107730014B (zh) 一种基于cbm的机队维修决策方法
Zeng et al. Flutter prediction for flight/wind-tunnel flutter test under atmospheric turbulence excitation
CN107577902B (zh) 一种基于ukf的飞机疲劳结构剩余寿命预测方法
Tibaduiza Burgos Design and validation of a structural health monitoring system for aeronautical structures.
Vidanović et al. Validation of the CFD code used for determination of aerodynamic characteristics of nonstandard AGARD-B calibration model
ul Hassan et al. Comparison of different life distribution schemes for prediction of crack propagation in an aircraft wing
DeCarlo et al. Bayesian calibration of aerothermal models for hypersonic air vehicles
Boller Structural health management of ageing aircraft and other infrastructure
Coppe et al. Least squares-filtered Bayesian updating for remaining useful life estimation
CN115408755B (zh) 一种考虑时变效应的组合梁桥动力疲劳可靠性评估方法
Aeschliman et al. A proposed methodology for computational fluid dynamics code verification, calibration, and validation
Xu et al. Determination of impact events on a plate-like composite structure
Niculescu et al. Integrated System for Monitoring Aircraft Structural Condition by Using the Strain Gauge Marks Method
Wu et al. Uncertainty analysis of flutter predictions with focus on the AGARD 445.6 wing
Coates et al. An inverse method for selection of Fourier coefficients for flight load identification
Kim et al. Bayesian approach for fatigue life prediction from field data
Coppe et al. Reducing uncertainty in damage growth properties by structural health monitoring
Gallagher et al. Developing normalized crack growth curves for tracking damage in aircraft
Ovchinnikov et al. Authenticity of the equivalent vibration tests
Shen et al. Overview of vibrational-based nondestructive evaluation techniques
Vecchio et al. Numerical Evaluation of Damage Distribution over a Slat Track Using Flight Test Data
Park et al. Modeling the effect of structural tests on uncertainty in estimated failure stress

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