文章快速检索  
  高级检索
变截面悬臂梁结构动载荷辨识方法
吴肖1, 曾捷1, 胡子康1, 李明2, 胡锡涛1     
1. 南京航空航天大学 机械结构力学及控制国家重点实验室, 南京 210016;
2. 中国航空综合技术研究所 装备服务产品部, 北京 100028
摘要: 在航空航天领域,作用在结构上动载荷的确定对结构健康监测是非常必要和重要的。为此,本文以类似机翼结构的变截面悬臂梁结构为研究对象,提出了一种基于光纤光栅传感器与卡尔曼滤波器的动载荷识别方法。首先,根据变截面梁单元形式,推导出变截面梁的质量矩阵与刚度矩阵,建立动力学运动方程。然后,以光纤光栅传感器测得的应变信息作为观测信号,通过卡尔曼滤波器生成的增益矩阵、新息序列矩阵以及协方差矩阵,得到灵敏度矩阵和估计力的增益矩阵。在此基础上,利用广义回归模型及其最小二乘算法,估算出动载荷大小、判断出动载荷激励位置。借助数值仿真与实验手段,分别验证了该方法对于单点正弦激励、方波激励、锯齿波激励以及多点同时激励等工况下的动载荷识别效果。结果表明,本文所提算法具有较好的动载荷识别效果和噪声抑制能力,能够为未来风洞试验和真实飞行试验环境中诸如大展弦比机翼表面气动压力等载荷实时辨识、气动外形自适应控制以及结构健康监测提供技术支撑。
关键词: 变截面悬臂梁结构    光纤光栅传感器    卡尔曼滤波器    动载荷识别    递归最小二乘法    
Dynamic load identification method for variable cross-section cantilever structure
WU Xiao1, ZENG Jie1, HU Zikang1, LI Ming2, HU Xitao1     
1. State Key Laboratory of Mechanical Structure Mechanics and Control, Nanjing University of Aeronautics and Astronautics, Nanjing 210016;
2. Equipment Service Product Department, China Aviation Integrated Technology Research Institute, Beijing 100028
Abstract: In the field of aerospace, it is both necessary and important to determine the dynamic load on the structure for structural health monitoring. In this paper, a method of dynamic load identification based on FBG sensors and the Kalman filter is proposed. According to the element form of variable cross-section beams, the mass matrix and stiffness matrix of variable cross-section beams are firstly derived, and the dynamic equation of motion established. Then, taking the strain information measured by FBG sensors as the observation signal, the gain matrix, innovation sequence matrix and covariance matrix generated by the Kalman filter are used to obtain the sensitivity matrix and the gain matrix of estimation force. On this basis, the generalized regression model and its least square algorithm are adopted to estimate the dynamic load and determine the excitation position of the dynamic load. By means of numerical simulation and experiments, the dynamic load identification effect of this method under single point sine excitation, square wave excitation, saw tooth wave excitation and multi-point simultaneous excitation is verified. Results show that the algorithm proposed has good dynamic load identification effect and noise suppression ability. It can also provide technical support for future wind tunnel tests and real flight test environments such as real-time load identification and aerodynamic shape adaptive control, and structural health monitoring such as high aspect ratio wing surface aerodynamic pressure.
Keywords: variable section cantilever structure    fiber Bragg grating (FBG) sensor    Kalman filter    dynamic load identification    recursive least square method    

载荷识别主要用于计算恶劣服役环境下结构主承力构件所受载荷或者无法通过传感器获取的外部激励[1]。动载荷识别一般分为直接测量法和间接识别法2种[2],直接测量法是通过在飞行器表面设置大量测压孔感知整个表面的压力分布。此方式特点是测量精度高、计算量小,但缺点是打孔给飞行器结构引入损伤,一定程度上影响了飞行器结构力学性能。

间接识别算法是在结构表面布置传感器网络直接采集载荷响应信息来反求载荷,此种方法不需破坏结构、操作较简单。国内外许多学者对此进行了大量深入的研究,Hong等提出了一种倒谱域的机械动态载荷谱识别方法,基于高阶谱的谐波恢复,建立了ARMA(Auto Regvessive Moving Avevage Model)模型,确定格林函数的系数,求出动载荷,但只是进行了仿真验证,缺乏工程验证[3]。Maes等基于杆件的模态特性,假设了Timoshenko梁理论,在考虑杆的转动惯量和传感器质量的基础上,提出了一种估算组合结构杆件轴向力的方法,但不适用于变截面等复杂结构的动载荷识[4]。Ma等提出了一种基于卡尔曼滤波和递推最小二乘算法的在线递推动载荷方法。此算法对高斯白噪声具有很好的抑制作用,通过正弦、三角脉冲、矩形脉冲、一系列脉冲和随机脉冲仿真,验证了该方法对由噪声测量值估计梁结构外部激励具有良好的性能,但未实现载荷作用位置辨识[5]。间接识别法属于数学上的逆问题,因而存在矩阵病态导致的识别精度不高和计算量偏大的问题[6]。目前关于间接识别法的研究中,除了矩阵病态问题,还有传感器的安装和布局等难点[7-8]

由于光纤传感技术能够实现针对应变、温度、振动等多种物理量的直接测量,具有便于分布式组网、灵敏度高、芯径细、重量轻、柔韧性好、适于强电磁干扰以及腐蚀性危险环境使用等独特优点[9]。基于此,本文研究了一种基于光纤光栅传感器与卡尔曼滤波器的变截面悬臂梁结构动载荷识别算法。该方法基于自适应卡尔曼滤波器反馈修正的思想进行载荷识别,一定程度上解决了间接识别法中矩阵求逆所导致的病态问题,给出了载荷所用位置判断依据,并通过仿真与实验验证了此算法具有较好的抗噪能力、较高的载荷辨识精度和工程适用性。

1 光纤光栅应变传感原理

光纤光栅传感器(Fiber Bragg Grating,FBG)基于光纤周期或纤芯折射率的变化,当光纤Bragg光栅感受到外界温度、应变变化时,会引起光栅栅区折射率的变化,从而使得其对应中心波长发生改变[10]。光纤光栅传感原理,如图 1所示。

图 1 光纤光栅传感器原理 Fig. 1 Principle of fiber Bragg grating sensor

由光纤光栅传感器原理可知,假定环境温度恒定,此时FBG仅受轴向应变因素影响,当中心波长变化量为ΔλB时,其应变表达式为

$ \varepsilon = \frac{1}{{1 - {P_{\rm{e}}}}} \cdot \frac{{\Delta {\lambda _{\rm{B}}}}}{{{\lambda _{\rm{B}}}}} $ (1)
 

式中:Pe为光纤光栅有效弹光系数,取为常数0.22; λB为FBG反射光谱中心波长; ΔλB为中心波长偏移量[11]

2 动载荷识别原理 2.1 变截面梁有限元分析

针对如图 2所示长度为l、厚度为h、上底边宽度为b0、下底边宽度为b1的锥形变截面梁单元形式,其宽度b(x)、截面积A(x)和惯性矩I(x)可以表示为[12]

图 2 锥形变截面梁单元示意图 Fig. 2 Schematic diagram of tapered variable section beam element
$ \left\{ {\begin{array}{*{20}{l}} {b(x) = {b_0}\left( {1 - \frac{x}{l}} \right) + {b_1}\frac{x}{l}}\\ {A(x) = h\left[ {{b_0}\left( {1 - \frac{x}{l}} \right) + {b_1}\frac{x}{l}} \right]}\\ {I(x) = \frac{{{h^3}}}{{12}}\left[ {{b_0}\left( {1 - \frac{x}{l}} \right) + {b_1}\frac{x}{l}} \right]} \end{array}} \right. $ (2)
 

根据结构力学的有限元分析可知,等截面梁单元刚度矩阵Ke可以表示为

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{K}}^e} = \int\limits_0^l {{\mathit{\boldsymbol{B}}^{\rm{T}}}} EI\mathit{\boldsymbol{B}}{\rm{d}}x = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{EI}}{{{l^3}}}\left[ {\begin{array}{*{20}{c}} {12}&{6l}&{ - 12}&{6l}\\ {6l}&{4{l^2}}&{ - 6l}&{2{l^2}}\\ { - 12}&{ - 6l}&{12}&{ - 6l}\\ {6l}&{2{l^2}}&{ - 6l}&{4{l^2}} \end{array}} \right]} \end{array} $ (3)
 

对于等截面梁单元的一致质量矩阵Me可以表示为

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{M}}^e} = \int_0^l {{\mathit{\boldsymbol{N}}^{\rm{T}}}} \rho A\mathit{\boldsymbol{N}}{\rm{d}}x = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\rho Al}}{{420}}\left[ {\begin{array}{*{20}{c}} {156}&{22l}&{54}&{ - 13l}\\ {22l}&{4{l^2}}&{13l}&{ - 3{l^2}}\\ {54}&{13l}&{156}&{ - 22l}\\ { - 13l}&{ - 3{l^2}}&{ - 22l}&{4{l^2}} \end{array}} \right]} \end{array} $ (4)
 

式中:EI为梁单元抗弯刚度;A为梁的截面面积;b为梁截面宽度[13]N为梁单元的形函数矩阵;ρ为材料密度;B为几何矩阵。将式(2)代入式(3)、式(4),可以得到变截面梁单元刚度矩阵与质量矩阵。

2.2 单元节点划分

通过有限元法将梁离散成n个单元,在每个单元内选取2个应变测量点,共计2n个应变测量点。由此,可通过形函数分别建立起每个单元内所测应变与单元节点的位移、转角之间的映射关系:

$ \left[ {\begin{array}{*{20}{c}} {{\varepsilon _1}}\\ {{\varepsilon _2}}\\ \vdots \\ {{\varepsilon _{2n - 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{B}}(1,{\xi _1})}&{}&{}&{}&{}&{}\\ {\mathit{\boldsymbol{B}}(2,{\xi _2})}&{}&{}&{}&{}&{}\\ {}&{\mathit{\boldsymbol{B}}(3,{\xi _3})}&{}&{}&{}&{}\\ {}&{\mathit{\boldsymbol{B}}(4,{\xi _4})}&{}&{}&{}&{}\\ {}&{}&\mathit{\boldsymbol{O}}&{}&{}&{}\\ {}&{}&{}&{\mathit{\boldsymbol{B}}(i,{\xi _i})}&{}&{}\\ {}&{}&{}&{}&\mathit{\boldsymbol{O}}&{}\\ {}&{}&{}&{}&{}&{\mathit{\boldsymbol{B}}(2n - 1,{\xi _{2n - 1}})}\\ {}&{}&{}&{}&{}&{\mathit{\boldsymbol{B}}(2n,{\xi _{2n}})} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _1}}\\ {{\theta _1}}\\ {{\omega _2}}\\ {{\theta _2}}\\ \vdots \\ {{\omega _{n - 1}}}\\ {{\theta _{n - 1}}}\\ {{\omega _n}}\\ {{\theta _n}} \end{array}} \right] $ (5)
 

式中:ε1~ε2nn个单元对应的2n个应变值,实验中可由光纤光栅传感器测得;ω1~ωnn个单元节点的位移;θ1nn个单元节点的转角;B(i, ξi)=[-6+12ξi, l(-4+6ξi), 6-12ξi, l(-2+6ξi)]×y/l2, ξ为光纤光栅传感器粘贴位置,y=h/2为梁中性轴距离其上表面的距离。

2.3 建立状态方程

离散化后的梁结构为二阶振动系统,其2n个自由度系统的振动方程为

$ \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}^{\prime \prime }}(t) + \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{X}}^\prime }(t) + \mathit{\boldsymbol{KX}}(t) = \mathit{\boldsymbol{F}}(t) $ (6)
 

式中:MCK分别为2n×2n的质量矩阵、阻尼矩阵和刚度矩阵;F(t)为2n×1的载荷向量;X(t)为系统的位移;X'(t)为系统的速度;X″(t)为系统的加速度。

对于阻尼矩阵在实模态空间内不能解耦的情况,引入状态坐标变换[14]:

$ \mathit{\boldsymbol{Y}}(t) = [\mathit{\boldsymbol{X}}(t)\quad {\mathit{\boldsymbol{X}}^\prime }(t)] $ (7)
 

将式(7)代入式(6),可把上述矩阵形式的微分方程转换为状态空间方程,表示为

$ {\mathit{\boldsymbol{Y}}^\prime }(t) = \left[ {\begin{array}{*{20}{c}} {{{\bf{0}}_{2n \times 2n}}}&{{\mathit{\boldsymbol{I}}_{2n \times 2n}}}\\ { - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{K}}}&{ - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{C}}} \end{array}} \right]\mathit{\boldsymbol{Y}}(t) + \left[ {\begin{array}{*{20}{c}} {{{\bf{0}}_{2n \times 2n}}}\\ {{\mathit{\boldsymbol{M}}^{ - 1}}} \end{array}} \right]\mathit{\boldsymbol{F}}(t) $ (8)
 

$ \mathit{\boldsymbol{D}} = \left[ \begin{gathered} {{\bf{0}}_{2n \times 2n}}\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{I}}_{2n \times 2n}} \hfill \\ - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{K}}\;\;\;\;\; - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{C}} \hfill \\ \end{gathered} \right], \mathit{\boldsymbol{G}} = \left[ \begin{gathered} {{\bf{0}}_{2n \times 2n}} \hfill \\ {\mathit{\boldsymbol{M}}^{ - 1}} \hfill \\ \end{gathered} \right]$,则状态方程可以表示为

$ {\mathit{\boldsymbol{Y}}^\prime }(t) = \mathit{\boldsymbol{DY}}(t) + \mathit{\boldsymbol{GF}}(t) $ (9)
 

式中:F(t)=[F1  F2  …  F2n]T; D为状态矩阵,由梁的质量矩阵、刚度矩阵和阻尼矩阵决定。系统观测方程可表示为

$ \mathit{\boldsymbol{Z}}(t) = \mathit{\boldsymbol{HY}}(t) $ (10)
 

式中:H为观测矩阵;Z(t)为观测序列,即结构的应变响应[15]

$ \mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}_1}}&{}&{}&{}&{}\\ {{\mathit{\boldsymbol{B}}_2}}&{}&{}&{}&{}\\ {}&{{\mathit{\boldsymbol{B}}_3}}&{}&{}&{}\\ {}&{{\mathit{\boldsymbol{B}}_4}}&{}&{}&{}\\ {}&{}&\mathit{\boldsymbol{O}}&{}&{}\\ {}&{}&{}&\mathit{\boldsymbol{O}}&{}\\ {}&{}&{}&{}&{{\mathit{\boldsymbol{B}}_{2n - 1}}}\\ {}&{}&{}&{}&{{\mathit{\boldsymbol{B}}_{2n}}} \end{array}} \right] $ (11)
 

采样时间间隔设置为ΔT,将式(9)和式(10)离散化,得到线性离散系统的状态方程和观测方程:

$ {\mathit{\boldsymbol{Y}}(k + 1) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} Y}}(k) + \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}(\mathit{\boldsymbol{F}}(k) + \mathit{\boldsymbol{w}}(k))} $ (12)
 
$ {\mathit{\boldsymbol{Z}}(k) = \mathit{\boldsymbol{HY}}(k) + \mathit{\boldsymbol{v}}(k)} $ (13)
 
$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = {\rm{exp}}(\mathit{\boldsymbol{D}}*\Delta T)} $ (14)
 
$ {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = \int_{(k - 1)\Delta T}^{k\Delta T} {{\rm{exp}}} [\mathit{\boldsymbol{D}}(k\Delta T - \tau )]\mathit{\boldsymbol{G}}{\rm{d}}\tau } $ (15)
 

式中:Y(k)为状态向量;Φ为状态转移矩阵(描述动态系统从时间k到时间k+1状态的转移);Ψ为驱动矩阵;w(k)和v(k)分别为互不相干的系统白噪声和观测白噪声,其系统特性为[16]:数学期望(系统白噪声)为0,即E[w(k)]=0;系统白噪声自适应函数为E[w(k)wT(k)]=Qkdkl;数学期望(观测白噪声)为0,即E[v(k)]=0;观测白噪声自适应函数为E[v(k)vT(k)]=Rkdkl, dkl为自相关函数;QR为协方差矩阵, Q=QwI2n×2nR=RvI2n×2n, QwRvw(k)白噪声给定的协方差初始值。

2.4 卡尔曼滤波器原理

卡尔曼滤波器可分为2个部分:时间更新预测和测量更新修正[17-18]

时间更新:将光纤光栅传感器所测应变响应作为输入变量进行时间方程更新。

$ \mathit{\boldsymbol{\bar Y}}(k/(k - 1)) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} \bar Y}}(k - 1/(k - 1)) $ (16)
 

当前时刻的状态预计:利用上一时刻的最优解Y(k-1/(k-1))与状态转移矩阵Φ,得到当前时刻的状态估计值Y(k/(k-1))。

$ \mathit{\boldsymbol{P}}(k/(k - 1)) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} P}}(k - 1/(k - 1)){\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}} + \mathit{\boldsymbol{ \boldsymbol{\varPsi} Q}}{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}} $ (17)
 

根据上一状态的协方差P(k-1/(k-1))和系统过程中的协方差Q,求出当前状态的协方差预测值P(k/(k-1))。

测量更新:

$ \mathit{\boldsymbol{S}}(k) = \mathit{\boldsymbol{HP}}(k/(k - 1)){\mathit{\boldsymbol{H}}^{\rm{T}}} + \mathit{\boldsymbol{R}} $ (18)
 

利用式(17)得到当前时刻的协方差预测值P(k/(k-1)),再结合观测矩阵H,求得新息协方差S(k)。

$ {\mathit{\boldsymbol{K}}_a}(k) = \mathit{\boldsymbol{P}}(k/(k - 1)){\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{S}}^{ - 1}}(k) $ (19)
 

利用新息协方差S(k)以及观测矩阵H,得到卡尔曼增益矩阵Ka(k)。

$ {\mathit{\boldsymbol{P}}(k/k) = [\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_a}(k)\mathit{\boldsymbol{H}}]\mathit{\boldsymbol{P}}(k/(k - 1))} $ (20)
 
$ {\mathit{\boldsymbol{\bar Z}}(k) = \mathit{\boldsymbol{Z}}(k) - \mathit{\boldsymbol{H\bar Y}}(k/(k - 1))} $ (21)
 

利用观测方程Z(k)、观测矩阵H以及当前时间的状态估计值Y(k/(k-1)),求得当前时间的新息序列Z(k)。

$ \mathit{\boldsymbol{\bar Y}}(k/k) = \mathit{\boldsymbol{\bar Y}}(k/(k - 1)) + {\mathit{\boldsymbol{K}}_a}(k)\mathit{\boldsymbol{\bar Z}}(k) $ (22)
 

利用当前时间的状态估计值Y(k/(k-1))、卡尔曼增益矩阵Ka(k)以及当前时间的新息序列Z(k),求得当前时间的最优解Y(k/k)[19]

2.5 载荷识别器原理

基于卡尔曼滤波器的载荷估计器是以光纤光栅传感器测得的应变响应作为观测信号,通过卡尔曼滤波器生成的增益矩阵、新息序列以及协方差矩阵,得到灵敏度矩阵和估计力的增益矩阵。在此基础上,利用广义回归模型及其最小二乘算法估计载荷特征。

基于卡尔曼滤波器的动载荷识别器方程如下所示[20]

$ {{\mathit{\boldsymbol{B}}_s}(k) = \mathit{\boldsymbol{H}}[\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{M}}_{\rm{s}}}(k - 1) + \mathit{\boldsymbol{I}}]\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} $ (23)
 
$ {{\mathit{\boldsymbol{M}}_{\rm{s}}}(k) = [\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_a}(k)\mathit{\boldsymbol{H}}][\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{M}}_{\rm{s}}}(k - 1) + \mathit{\boldsymbol{I}}]} $ (24)
 
$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_{\rm{b}}}(k) = {\gamma ^{ - 1}}{\mathit{\boldsymbol{P}}_{\rm{b}}}(k - 1)\mathit{\boldsymbol{B}}_{\rm{s}}^{\rm{T}}(k)[{\mathit{\boldsymbol{B}}_{\rm{s}}}(k){\gamma ^{ - 1}} \cdot }\\ {{\mathit{\boldsymbol{P}}_{\rm{b}}}(k - 1)\mathit{\boldsymbol{B}}_{\rm{s}}^{\rm{T}}(k) + \mathit{\boldsymbol{S}}(k){]^{ - 1}}} \end{array} $ (25)
 
$ {{\mathit{\boldsymbol{P}}_{\rm{b}}}(k) = [\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_{\rm{b}}}(k){\mathit{\boldsymbol{B}}_{\rm{s}}}(k)]{\gamma ^{ - 1}}{\mathit{\boldsymbol{P}}_{\rm{b}}}(k - 1)} $ (26)
 
$ {\mathit{\boldsymbol{\hat F}}(k) = \mathit{\boldsymbol{\hat F}}(k - 1) + {\mathit{\boldsymbol{K}}_{\rm{b}}}(k)[\mathit{\boldsymbol{\bar Z}}(k) - {\mathit{\boldsymbol{B}}_{\rm{s}}}(k)\mathit{\boldsymbol{\hat F}}(k - 1)]} $ (27)
 

式中:Bs(k)和Ms(k)代表灵敏度矩阵;Kb(k)是用来估计力的增益矩阵;$ \mathit{\boldsymbol{\hat F}}\left( k \right)$为此时刻的激振力的估计值;γ是衰减因子用来调节自适应能力。

基于光纤光栅传感器和卡尔曼滤波器的动载荷识别流程,如图 3所示。

图 3 基于光纤光栅传感器与卡尔曼滤波器的动载荷识别流程 Fig. 3 Dynamic load identification process based on FBG sensor and Kalman filter
3 有限元数值仿真 3.1 变截面悬臂梁结构建模

利用Ansys Workbench对锥形梁单元形式的变截面梁结构进行有限元数值仿真,材料密度为2 770 kg/m3,弹性模量为71 GPa,泊松比为0.3,上底边长为100 mm,下底边长250 mm,模型长度为1 200 mm,厚度为6 mm。

将变截面悬臂梁结构模型平均划分为3个单元,从左到右依次定义相应3个单元节点:A~ C。在每个单元内的1/3处和2/3处设置应变监测提取点,则共有6个应变提取点,从固支端开始依次将其定义为1~6。单元节点划分与应变提取点位置,如图 4所示,其中矩形标记为单元内应变提取点,圆形标记为单元节点。

图 4 单元节点划分与测点位置图 Fig. 4 Unit node division and measuring point location

为验证此载荷辨识算法的抗噪能力,在提取应变响应后加上σ=1×10-9Qw=1×10-15的高斯白噪声作为系统的观测序列代入相关算法。设置衰减因子γ=0.65,采样频率为500 Hz。

假设变截面悬臂梁结构系统阻尼为比例阻尼,阻尼矩阵是质量矩阵和刚度矩阵的线性组合,可以表示为

$ \mathit{\boldsymbol{C}} = \alpha \mathit{\boldsymbol{M}} + \beta \mathit{\boldsymbol{K}} $ (28)
 

式中:α为质量比例常数,β为刚度比例常数,设置α=0.01,β=0.02。

为进一步描述动载荷识别效果,量化动载荷识别误差,采用平均绝对误差MAE和均方根误差RMSE来描述动载荷识别误差。

$ {{\rm{MAE}} = \frac{1}{N} \cdot \sum\limits_{i = 1}^N | {\mathit{\boldsymbol{F}}_i} - {{\mathit{\boldsymbol{\hat F}}}_i}|} $ (29)
 
$ {{\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^N {{{({\mathit{\boldsymbol{F}}_i} - {{\mathit{\boldsymbol{\hat F}}}_i})}^2}} }}{{N - 1}}} } $ (30)
 

式中:Fi为在梁结构上施加的模拟真实激励;${\mathit{\boldsymbol{\hat F}}_i} $为算法识别的载荷;N为样本总数。

3.2 有限元仿真结果验证 3.2.1 单点加载下动载荷辨识效果

1) 动载荷大小识别

依次在单元节点ABC施加正弦激励:F=5sin(3πt),识别效果如图 5所示。

图 5 载荷识别效果 Fig. 5 Load identification effect

根据图 5可知,单点载荷作用下,单元节点识别载荷效果良好,载荷识别曲线与仿真施加载荷曲线基本吻合,各单元节点载荷识别误差相差不大。单元节点A处的动载荷识别精度较好,单元节点B和单元节点C处的载荷识别曲线在峰值处略微偏小,这主要是因为悬臂梁模型在外加载荷激励过程中,距离固支端越远的测点,提取的应变响应受到外部激励扰动越明显,使得载荷辨识精度相对于固支端有所降低。

单点激励下各单元节点载荷识别误差,如表 1所示。

表 1 单点激励下各单元节点载荷识别误差 Table 1 Load identification errors of each element node under single point excitation
载荷施加位置 MAE/N RMSE/N
单元节点A 0.106 2 0.276 5
单元节点B 0.199 9 0.408 6
单元节点C 0.248 0 0.579 4

2) 动载荷激励位置识别

根据式(5)函数关系,本文提出通过提取变截面模型每个单元内2个位置的应变响应来计算相应单元节点处的位移和转角作为载荷识别算法的输入变量,而输出值是每个单元节点对应的载荷与弯矩。因此,算法识别出来的载荷与其相应单元节点呈一一对应关系,也即在载荷激励位置,将会存在显著的载荷响应,据此可以判断出动载荷施加位置。

以在单元节点B处施加激励载荷为例,单元节点B识别的有效载荷和单元节点A、单元节点C识别到的扰动载荷响应曲线,如图 6所示。可以发现,在单元节点B(距离固支端0.8 m)处识别的载荷显著大于单元节点A(距离固支端0.4 m)、单元节点C(距离固支端1.2 m),且单元节点A、单元节点C处识别的载荷波动很小,无明显变化规律,因此可以判断载荷施加位置为单元节点B处,这与模拟仿真所加载荷位置相符。

图 6 各单元节点载荷识别曲线 Fig. 6 Load identification curves of each element node
3.2.2 三单元节点同时加载下动载荷辨识效果

为进一步验证此算法在复合工况下动载荷识别效果,选择在单元节点A、单元节点B和单元节点C处同时施加3个正弦激励,具体加载情况如表 2所示。

表 2 三单元节点同时加载工况 Table 2 Three-unit joint simultaneous loading conditions
加载单元节点 加载大小
第1单元节点A FA=-20sin(4πt)
第2单元节点B FB=30sin(6πt)
第3单元节点C FC=10sin(2πt)

3个单元节点同时加载条件下动载荷识别效果,如图 7所示。从图 7可以看出,3个单元节点同时加载工况下,单元节点ABC处对应的载荷识别曲线与有限元仿真结果基本吻合,仅在峰值处均与存在一定偏差。单元节点C处识别的载荷受其他激励载荷影响,识别曲线出现一定的波动。相对于单点激励工况,由于加载幅值的增加,3个单元节点的载荷识别误差略有增大。

图 7 3个单元节点同时加载下动载荷识别效果 Fig. 7 Dynamic load identification effect under simultaneous loading of three element nodes

3个单元节点同时加载条件下动载荷识别误差,如表 3所示。

表 3 3个单元节点同时加载下动载荷识别误差 Table 3 Identification errors of dynamic load under simultaneous loading of three element nodes
载荷施加位置 MAE/N RMSE/N
单元节点A 0.208 4 0.302 1
单元节点B 0.351 4 0.479 5
单元节点C 0.408 0 0.550 3
4 实验验证 4.1 变截面悬臂梁结构动载荷识别实验系统

变截面悬臂梁结构单元节点划分与传感器布局示意图,如图 8所示。变截面悬臂梁结构分布式光纤动态载荷识别实验系统,如图 9所示。

图 8 单元节点划分与传感器布局 Fig. 8 Unit node division and sensor layout
图 9 分布式光纤动载荷监测与识别实验系统 Fig. 9 Distributed optical fiber dynamic load monitoring and identification experimental system

实验平台主要由光纤光栅传感器、光纤光栅解调仪、计算机、33500B型信号发生器、HEAS-50型功率放大器、力显示器、CZL-1012型力传感器、NI USB-6341数据采集卡、HEV-50型激振器等组成。

4.2 单点加载工况下动载荷辨识效果

在单元节点C和单元节点B进行单点激励,通过算法识别出来的动载荷与力传感器直接测得的真实激励进行数据分析与对比,激励工况及动载荷识别效果如下。

4.2.1 正弦激励

单元节点B正弦激励(7 Hz)动载荷识别效果,如图 10所示。由图 10(b)可以判别载荷激励位置为单元节点B

图 10 单元节点B正弦激励(7 Hz)动载荷识别效果 Fig. 10 Dynamic load identification effect of unit node B with sinusoidal excitation (7 Hz)
4.2.2 方波激励

单元节点C处方波激励(1 Hz)动载荷识别效果,如图 11所示。由图 11(b)可以判别载荷激励位置为单元节点C

图 11 单元节点C方波激励(1 Hz)动载荷识别效果 Fig. 11 Identification effect of dynamic load excited by C square wave (1 Hz) of unit node
4.2.3 锯齿波激励

单元节点C处锯齿波激励(1 Hz)动载荷识别效果,如图 12所示。由图 12(b)可以判别载荷激励位置为单元节点C

图 12 单元节点C锯齿波激励(1 Hz)动载荷识别效果 Fig. 12 Identification effect of dynamic load excited by sawtooth wave (1 Hz) of element node C

根据图 10~图 12可知,在单点加载工况下,载荷识别曲线与真实激励曲线吻合度较好,只是在最大幅值处出现一定偏差。由于在施加方波、锯齿波时激振器会产生一定的顿挫,在顿挫瞬间,激励本身出现振荡,使得载荷识别曲线相较于实际施加载荷曲线存在一定波动。

通过各单元节点载荷识别曲线,根据每个单元节点对应的识别载荷可以判断出真实激励单元节点。

单点激励下,不同频率与波形激励载荷对应的动载荷辨识误差,如表 4所示。

表 4 不同频率与波形激励载荷对应的载荷辨识误差 Table 4 Load identification errors corresponding to different frequencies and waveform excitation loads
波形 载荷施加位置 MAE/N RMSE/N
正弦(7 Hz) 单元节点B 0.486 5 0.542 3
方波(1 Hz) 单元节点C 0.490 5 0.635 1
锯齿波(1 Hz) 单元节点C 0.504 4 0.653 4
4.3 两单元节点同时加载工况下动载荷辨识效果

在单元节点B(激励频率2 Hz)和单元节点C(激励频率4 Hz)同时施加外部激励荷载,载荷识别效果如图 13所示。

图 13 单元节点B、C同时激励下动载荷识别效果 Fig. 13 Identification effect of dynamic load under simul taneous excitation of B and C nodes

根据图 13所示,在两点同时激励的复合动载荷加载工况中,载荷辨识曲线与实际加载曲线轮廓基本吻合,算法能够较好地识别出单元节点B和单元节点C处对应的外部激励载荷特征,具有良好的识别效果。根据图 14可以判断出单元节点载荷激励位置,分别为单元节点B和单元节点C

图 14 单元节点载荷激励位置识别 Fig. 14 Location identification of load excitation of element node

动载荷识别结果中呈现2个激励频率波形相互叠加的情况,这是因为在单元节点B和单元节点C处同时激励的过程中,2个激振器通过变截面悬臂梁结构相互作用,激励本身出现耦合所致。这表明该算法针对复合载荷同时加载工况,具有较好的动载荷识别能力。两单元节点同时激励实验误差,如表 5所示。

表 5 两单元节点同时激励对应的载荷辨识误差 Table 5 Load identification errors of two element nodes simultaneously excited
组合情况 载荷施加位置 MAE/N RMSE/N
正弦(2 Hz) 单元节点B 0.446 1 0.636 2
正弦(4 Hz) 单元节点C 0.429 1 0.610 5

综合上述分析,动载荷识别结果也存在一定偏差,分析其原因可能包括:一是系统结构离散的单元数较少;二是光纤光栅传感器与变截面悬臂梁结构上表面未均匀胶结或粘贴方向有偏差,引起应变测量误差;三是光纤光栅传感器中心波长偏移量受结构振动扰动影响较大,尤其是在最大振幅时识别精度有所下降。

5 结论

本文针对飞行器结构载荷监测需求,推导了变截面梁模型质量矩阵、刚度矩阵等参数信息,提出了基于分布式光纤传感器与卡尔曼滤波器的变截面悬臂梁动载荷辨识方法。

1) 基于自适应卡尔曼滤波器反馈修正的思想进行载荷识别,一定程度上解决了常规间接识别法中矩阵求逆所导致的病态问题。

2) 借助ANSYS有限元分析软件,开展针对结构在不同载荷激励状态下的数值仿真,验证了该算法对于多种不同载荷形式辨识的可行性以及良好的噪声抑制能力。

3) 基于形函数建立的变截面悬臂梁单元所测应变与单元节点的位移、转角之间映射关系,提出了动载荷激励位置判别方法。

4) 构建了基于分布式光纤传感器的变截面悬臂梁应变监测与动载荷识别系统,试验结果表明,该算法对于复合载荷工况具有较好的识别能力和工程适用性。

5) 本文所述方法采用全分布式光纤传感网络,相较于电学式测量方法,具有抗电磁干扰、质量轻不改变被测结构动力学特征、适用于狭小空间机载环境的非视觉测量等特点,能够为风洞和飞行试验中飞行器载荷辨识提供技术支撑。

参考文献
[1] JIANG R, JIANG R J. Identification of dynamic load and vehicle parameters based on bridge dynamic responses[J]. Journal of Geophysical Research Atmospheres, 2003, 83(B7): 3535-3538.
Click to display the text
[2] PECORA R, CONCILIO A, DIMINO I, et al. Structural design of an adaptive wing trailing edge for enhanced cruise performance[C]//24th AIAA/AHS Adaptive Structures Conference.Reston: AIAA, 2016.
Click to display the text
[3] HONG C H, QIAO S Y, WU M. Simulating study of dynamic load spectra identification method of machinery in cepstrum domain[J]. Journal of China University of Mining and Technology(English Edition), 2006, 16(1): 22-24.
Click to display the text
[4] MAES K, PEETERS J, REYNDERS E, et al. Identification of axial forces in beam members by local vibration measurements[J]. Journal of Sound and Vibration, 2013, V332: 5417-5432.
Click to display the text
[5] MA C K, CHANG J M, LIN D C. Input forces estimation of beam structures by an inverse method[J]. Journal of Sound and Vibration, 2003, 259(2): 387-407.
Click to display the text
[6] LIU J, MENG X, ZHANG D, et al. An efficient method to reduce ill-posedness for structural dynamic load identification[J]. Mechanical Systems and Signal Processing, 2017, 95: 273-285.
Click to display the text
[7] 杨智春, 贾有. 动载荷的识别方法[J]. 力学进展, 2015, 45(1): 29-54.
YANG Z C, JIA Y. Identification method of dynamic load[J]. Progress in Mechanics, 2015, 45(1): 29-54. (in Chinese)
Cited By in Cnki (47) | Click to display the text
[8] 兑红娜, 王勇军, 董江, 等. 基于飞行参数的飞机结构载荷最优回归模型[J]. 航空学报, 2018, 39(11): 80-89.
DUI H N, WANG Y J, DONG J, et al. Optimal regression model of aircraft structural load based on flight parameters[J]. Acta Aeronautics et Astronauctics Sinica, 2018, 39(11): 80-89. (in Chinese)
Cited By in Cnki (5) | Click to display the text
[9] 刘铁根, 王双, 江俊峰, 等. 航空航天光纤传感技术研究进展[J]. 仪器仪表学报, 2014(8): 1681-1692.
LIU T G, WANG S, JIANG J F, et al. Research progress of aerospace optical fiber sensing technology[J]. Journal of Instrumentation, 2014(8): 1681-1692. (in Chinese)
Cited By in Cnki (137) | Click to display the text
[10] WEI F F, LIU D J, MALLIK A K, et al. Magnetic field sensor based on a Tri-Microfiber coupler ring in magnetic fluid and a fiber bragg grating[J]. Sensors (Basel, Switzerland), 2019, 19(23): 5100.
Click to display the text
[11] XU B, HUANG J, XU X F, et al. Ultrasensitive no gas sensor based on the graphene oxide-coated long-period fiber grating[J]. ACS Applied Materials and Interfaces, 2019, 11(43): 40868-40874.
Click to display the text
[12] 王晓臣, 蒲军平. 变截面梁有限元分析[J]. 浙江工业大学学报, 2008, 36(3): 312-314.
WANG X C, PU J P. Finite element analysis of variable cross section beam[J]. Journal of Zhejiang University of Technology, 2008, 36(3): 312-314. (in Chinese)
Cited By in Cnki (53) | Click to display the text
[13] 胡海岩. 机械振动基础[M]. 北京: 北京航空航天大学出版社, 2018.
HU H Y. Foundation of mechanical vibration[M]. Beijing: Beijing University of Aeronautics and Astronautics Press, 2018. (in Chinese)
[14] 马志贵, 施文龙. 基于MATLAB的变截面梁单元有限元分析[J]. 兰州工业学院学报, 2015, 22(5): 37-39.
MA Z G, SHI W L. Finite element analysis of variable section beam element based on MATLAB[J]. Journal of Lanzhou Institute of Technology, 2015, 22(5): 37-39. (in Chinese)
Cited By in Cnki (4) | Click to display the text
[15] 宋雪刚.基于光纤光栅传感器的动载荷识别技术研究[D].南京: 南京航空航天大学, 2018.
SONG X G. Research on dynamic load identification technology based on FBG sensor[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2018(in Chinese).
Cited By in Cnki | Click to display the text
[16] GUO J, HUANG W, WILLIAMS B M. Adaptive Kalman filter approach for stochastic short-term traffic flow rate prediction and uncertainty quantification[J]. Transportation Research Part C:Emerging Technologies, 2014, 43: 50-64.
Click to display the text
[17] SUN F, TANG L J. Cubature Kalman filter-Kalman filter algorithm[J]. Control & Decision, 2012, 27(10): 1561-1565.
Click to display the text
[18] 侯康.基于卡尔曼滤波器的桥梁损伤快速诊断方法研究[D].哈尔滨: 哈尔滨工业大学, 2016.
HOU K. Research on fast diagnosis method of bridge damage based on Kalman filter[D]. Harbin: Harbin University of Technology, 2016(in Chinese).
Cited By in Cnki (2) | Click to display the text
[19] AZIZ K, EMRE K. Optimizing a Kalman filter with an evolutionary algorithm for nonlinear quadrotor attitude dynamics[J]. Journal of Computational Science, 2020, 39: 101051.
Click to display the text
[20] 宋雪刚, 刘鹏, 程竹明, 等. 基于光纤光栅传感器和卡尔曼滤波器的载荷识别算法[J]. 光学学报, 2018, 38(3): 2-3.
SONG X G, LIU P, CHENG Z M, et al. Load identification algorithm based on FBG sensor and Kalman filter[J]. Acta Optica Sinica, 2018, 38(3): 2-3. (in Chinese)
Cited By in Cnki (5) | Click to display the text
http://dx.doi.org/10.7527/S1000-6893.2020.23806
中国航空学会和北京航空航天大学主办。
0

文章信息

吴肖, 曾捷, 胡子康, 李明, 胡锡涛
WU Xiao, ZENG Jie, HU Zikang, LI Ming, HU Xitao
变截面悬臂梁结构动载荷辨识方法
Dynamic load identification method for variable cross-section cantilever structure
航空学报, 2020, 41(9): 223806.
Acta Aeronautica et Astronautica Sinica, 2020, 41(9): 223806.
http://dx.doi.org/10.7527/S1000-6893.2020.23806

文章历史

收稿日期: 2020-01-07
退修日期: 2020-02-06
录用日期: 2020-04-08
网络出版时间: 2020-04-13 09:42

相关文章

工作空间