在机械系统设计过程中, 动载荷是机械结构可靠性计算和疲劳分析的依据.但是, 很多结构内部的动载荷是很难直接测量的, 因此动载荷识别方法的研究具有重要意义.现有动载荷识别方法主要有两类:频域法[1-3]和时域法[4].频域法识别精度低, 对噪声敏感[5];时域法识别精度高, 抗噪性强.状态空间法[6]通常用于识别固定周期载荷;灵敏度法[7]将输入力近似为傅里叶级数, 用响应对载荷的灵敏度求解力; 零阶系统马尔可夫参数也被用于求解逆动力学问题[8-9]; 卡尔曼滤波法[10-11]和显式块反演算法[12]同样采用状态空间法进行载荷识别.状态空间法应用广泛, 然而作为显示方法, 它是条件稳定的, 其识别精度受采样频率和采样持续时长的限制.
本文基于Houbolt法提出了一种新的动载荷识别方法.由于Houbolt法本质上是隐式积分格式, 因此舍入误差与步长大小无关, 具有无条件稳定的优点, 同时本方法提高了识别精度.
1 基于Houbolt法的动载荷识别多自由度线性系统的运动微分方程为
(1) |
其中, M, C和K分别是系统质量矩阵、阻尼矩阵和刚度矩阵.
采用Houbolt法对式(1)进行逐步积分, 即利用两个向后差分公式表示在时刻t+Δt的速度和加速度为[13]
(2) |
(3) |
在t+Δt时刻, 将式(2)和式(3)代入方程(1), 可得
(4) |
其中,
由于Houbolt不是自起步的, 所以利用中心差分法对式(1)离散化处理, 可得
(5) |
(6) |
其中,
在式(4)中, 令t=2Δt, 可得
(7) |
将式(5)和式(6)以及初始条件代入式(7), 可得
(8) |
依次递推, 可得
(9) |
其中, 当1≤k≤3, 有
当4≤k≤n, 有
将式(9)表示成矩阵形式, 得
(10) |
其中:
在零初始条件下式(10)简化为
(11) |
考虑到可能存在的不适定性, 利用吉洪诺夫正则化方法计算动载荷F, 系统误差记为[14]
(12) |
为了找到e的最小值, 引入惩罚函数J为
(13) |
当F的一阶导数为零时, e达到最小值, 所以F的解为
(14) |
其中, λ为正则化系数, 由差异原则确定[15].
对于一个给定的系统, 基于Houbolt法的载荷识别方法可归纳为以下步骤:
1) 确定系统动力学参数M, C, K;
2) 由系统动力学参数计算矩阵A1, A2, A3, A4, B1, B2, B3, C1, C2, C3;
3) 分别计算式(9)中的矩阵Hk, k和Dk, k, 并组集得到式(10)中的矩阵H和D;
4) 选择采样频率, 由实验或仿真测得系统位移响应, 得到位移响应向量Y;
5) 根据式(14)计算待识别载荷向量F.
2 算例一线弹性悬臂梁, 长0.6 m, 宽0.06 m, 厚0.03 m.弹性模量E=3×1011 N/m, 密度ρ=2 800 kg/m3, 比例阻尼系数α=1.327, β=3.146×10-6, 泊松比μ=0.3.利用ANSYS建立有限元模型, 如图 1所示, 可求得系统的质量矩阵、刚度矩阵和阻尼矩阵.
在上述悬臂梁结构中第6号节点处施加激励F6(t)=400sin(3πt).选取采样频率为20 Hz, 采样时间为5 s, 由有限元仿真计算, 测得悬臂梁第5号、第9号节点处的横向位移响应时间历程, 作为输入数据.
根据施载点和响应点的位置提取矩阵H对应的行与列, 分别基于Houbolt法与传统的状态空间法计算得到动载荷F6的时间历程识别结果, 如图 2所示.
识别载荷与实际载荷的相对误差为
(15) |
其中, Fid和Freal分别是识别载荷和实际载荷.
上述工况中, 其他条件不变, 在第6号节点处施加激励F′6=400sin(3πt)+300sin(1.5πt), 同理计算可得动载荷F′6时域识别结果如图 4所示, 图 5为频域识别结果.
由上述识别结果分别计算时域上F6和F′6平均相对误差和频域上F′6对应频率成分的幅值相对误差, 结果如表 1和表 2所示.
由图 2~图 5及表 1~表 2的载荷识别结果对比可知, 本文提出的基于Houbolt法识别得到的载荷时间历程更加接近实际载荷, 识别相对误差比状态空间法更小.对于两种方法来说, 较大的相对误差均出现在波峰和波谷处, 大的波动会给载荷识别带来一定影响, 结果显示Houbolt法受影响较小, 而状态空间法识别得到的载荷此处相对误差较大.另一方面, 在频域上, 本文提出的方法识别得到的载荷对应频率成分的幅值大小也更接近实际载荷, 且没有多余的频率成分或频率丢失.而状态空间法识别得到的载荷在频域上含有多余频率成分, 且对应幅值大小与实际载荷相对误差较大.
2.2 多输入激励识别上述悬臂梁结构中, 其他条件不变, 在第6号和第9号节点处分别同时施加激励F6和F9, 选取采样频率为20 Hz, 采样时间为5 s, 测得悬臂梁第5号、第8号和第10号节点处的横向位移响应作为输入数据.根据施载点和响应点的位置提取矩阵H对应的行与列, 分别基于Houbolt法与传统的状态空间法计算F6和F9, 所得的时域识别结果如图 6和图 7所示, 频域识别结果如图 8和图 9所示.
(16) |
由上述F6和F9识别结果分别计算时域上平均相对误差如表 3所示, 以及频域上对应频率成分的幅值相对误差如表 4及表 5所示.
由图 6~图 9及表 3~表 5载荷识别结果分析可知, 本文提出的方法识别得到的第6号和第9号节点处载荷时间历程更加接近实际载荷, 相对误差更小, 且在波峰波谷处相对更加稳定, 而状态空间法识别得到的载荷有较大的波动.同时在频域上, 本文提出的方法识别得到的载荷对应频率成分的幅值大小更接近于实际载荷, 没有多余频率成分或漏频现象, 而状态空间法识别得到的载荷含有多余的不可忽略的频率成分, 且对应幅值大小与实际载荷相对误差较大.
2.3 多输入含噪声激励识别上节工况中, 其他条件不变, 加入10%随机噪声干扰, 计算第6号和第9号节点处载荷识别相对误差如图 10所示, 并由此计算得到相对误差如表 6所示.
由图 10及表 6载荷识别结果可知, 本文提出的方法在含噪声环境下识别相对误差总体上小于传统的状态空间法.因此本方法在载荷识别中受噪声干扰更小, 更加具有抗噪性和抗干扰能力.
3 结论1) 本文提出了基于Houbolt法的载荷识别方法, 推导出多自由度线性振动系统激励与系统动力学参数及响应递推表达式, 由正则化方法解得待识别动载荷.此方法克服了传统的状态空间法识别精度受采样频率和采样时长限制的缺点, 同时具有无条件稳定的优点.
2) 通过算例验证, 相比于传统的状态空间法, 本文提出的方法识别得到的动载荷在时域和频域上均更加接近实际动载荷, 识别相对误差均更小, 能够更加精确识别结构上的单输入及多输入激励载荷, 且克服了状态空间法识别载荷的增频现象.
3) 本方法在含一定随机噪声环境下, 能够更加准确地识别结构动载荷, 识别相对误差更小, 具有良好的抗噪性和抗干扰能力.
[1] |
Lage Y E, Maia N M M, Neves M M, et al.
Force identification using the concept of displacement transmissibility[J]. Journal of Sound and Vibration, 2013, 332(7): 1674–1686.
DOI:10.1016/j.jsv.2012.10.034 |
[2] |
Ma C, Hua H X, Xiang F.
The identification of external forces for a nonlinear vibration system in frequency domain[J]. Proceedings of the Institution of Mechanical Engineers, Part C:Journal of Mechanical Engineering Science, 2014, 228(9): 1531–1539.
DOI:10.1177/0954406213509085 |
[3] |
Obrien E J, Mcgetrick P J, Gonzalez A.
A drive-by inspection system via vehicle moving force identification[J]. Smart Structures and Systems, 2014, 13(5): 821–848.
DOI:10.12989/sss.2014.13.5.821 |
[4] |
Law S S, Zhu X Q.
Moving load identification problems and applications[J]. Current Opinion in Infectious Diseases, 2007, 20(1): 3–10.
DOI:10.1097/QCO.0b013e328012c5cd |
[5] |
Doyle J F.
Further developments in determining the dynamic contact law[J]. Experimental Mechanics, 1984, 24(4): 265–270.
DOI:10.1007/BF02323986 |
[6] |
Law S S, Fang Y L.
Moving force identification:optimal state estimation approach[J]. Journal of Sound and Vibration, 2001, 239(2): 233–254.
DOI:10.1006/jsvi.2000.3118 |
[7] |
Lu Z R, Law S S.
Force identification based on sensitivity in time domain[J]. Journal of Engineering Mechanics, 2006, 132(10): 1050–1056.
DOI:10.1061/(ASCE)0733-9399(2006)132:10(1050) |
[8] |
Kammer D.
Input force reconstruction using a time domain technique[J]. Journal of Vibration and Acoustics, 2013, 120(4): 868–874.
|
[9] |
Law S S, Bu J Q, Zhu X Q.
Time-varying wind load identification from structural responses[J]. Engineering Structures, 2005, 27(10): 1586–1598.
DOI:10.1016/j.engstruct.2005.05.007 |
[10] |
Liu J J, Ma C K, Kung I C, et al.
Input force estimation of a cantilever plate by using a system identification technique[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190(11): 1309–1322.
|
[11] |
Lei Y, Wu D T, Lin S Z.
Integration of decentralized structural control and the identification of unknown inputs for tall shear building models under unknown earthquake excitation[J]. Engineering Structures, 2013, 52(9): 306–316.
|
[12] |
Nordberg T P, Gustafsson I.
Dynamic regularization of input estimation problems by explicit block inversion[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(44/45/46/47): 5877–5890.
|
[13] |
张义民.
机械振动[M]. 北京: 清华大学出版社, 2007: 237-238.
( Zhang Yi-min. Mechanical vibration[M]. Beijing: Tsinghua University Press, 2007: 237-238. ) |
[14] |
Brezinski C, Redivo-Zaglia M, Rodriguez G, et al.
Multi-parameter regularization techniques for ill-conditioned linear systems[J]. Numerische Mathematik, 2003, 94(2): 203–228.
DOI:10.1007/s00211-002-0435-8 |
[15] |
Yang X J, Wang L.
A modified Tikhonov regularization method[J]. Journal of Computational and Applied Mathematics, 2015, 288: 180–192.
DOI:10.1016/j.cam.2015.04.011 |