, 黄清云, 张义民东北大学 机械工程与自动化学院, 辽宁 沈阳 110819
收稿日期：2017-05-08
基金项目：国家重点基础研究发展计划项目(2014CB046303)。
作者简介：赵春雨(1963-)，男，辽宁黑山人，东北大学教授，博士生导师。
摘要：利用Houbolt逐步积分和向后差分法, 推导出多自由度线性振动系统的响应与系统动力学参数及激励的关系递推表达式, 进而提出了一种无条件稳定的动载荷识别方法.以一线弹性悬臂梁为例, 利用单频和双频单输入激励、双频双输入激励响应的仿真结果, 对其激励载荷进行了识别.将识别载荷的时域误差和频谱成分与传统的状态空间法识别结果进行了对比, 验证了本载荷识别方法的时域和频域相对误差均更小, 且无增频现象, 提高了识别精度, 且具有良好的抗噪性.
关键词：Houbolt法振动系统无条件稳定性载荷识别相对误差
Dynamic Load Identification Based on the Houbolt Method
ZHAO Chun-yu
, HUANG Qing-yun, ZHANG Yi-minSchool of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Corresponding author: ZHAO Chun-yu, E-mail: chyzhao@mail.neu.edu.cn
Abstract: The formulae of the dynamic response related to the excitation and dynamic parameters for a multiple-degree-of-freedom linear vibration system were deduced by virtue of the Houbolt integration and backward difference methods. Then a new method for load identification of unconditional stability was proposed. Taking a linear elastic cantilever beam as an example, the exciting loads were identified using the simulation results of the single inputs of single frequency and dual frequencies, and the dual inputs of dual frequencies. Compared with the conventional state space method, the relative errors and frequency spectrum of the loads identified by the proposed method were smaller in both time and frequency domains, and there was no frequency increasing. This proposed method has improved identification accuracy and has a good anti-noise performance.
Key words: Houbolt methodvibration systemunconditional stabilityload identificationrelative error
在机械系统设计过程中, 动载荷是机械结构可靠性计算和疲劳分析的依据.但是, 很多结构内部的动载荷是很难直接测量的, 因此动载荷识别方法的研究具有重要意义.现有动载荷识别方法主要有两类:频域法[1-3]和时域法[4].频域法识别精度低, 对噪声敏感[5]；时域法识别精度高, 抗噪性强.状态空间法[6]通常用于识别固定周期载荷；灵敏度法[7]将输入力近似为傅里叶级数, 用响应对载荷的灵敏度求解力; 零阶系统马尔可夫参数也被用于求解逆动力学问题[8-9]; 卡尔曼滤波法[10-11]和显式块反演算法[12]同样采用状态空间法进行载荷识别.状态空间法应用广泛, 然而作为显示方法, 它是条件稳定的, 其识别精度受采样频率和采样持续时长的限制.
本文基于Houbolt法提出了一种新的动载荷识别方法.由于Houbolt法本质上是隐式积分格式, 因此舍入误差与步长大小无关, 具有无条件稳定的优点, 同时本方法提高了识别精度.
1 基于Houbolt法的动载荷识别多自由度线性系统的运动微分方程为
 | (1) |
采用Houbolt法对式(1)进行逐步积分, 即利用两个向后差分公式表示在时刻t+Δt的速度和加速度为[13]
 | (2) |
 | (3) |
 | (4) |
 |
 | (5) |
 | (6) |
 |
 | (7) |
 | (8) |
 | (9) |
 |
 |
 |
 | (10) |
 |
 | (11) |
 | (12) |
 | (13) |
 | (14) |
对于一个给定的系统, 基于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所示, 可求得系统的质量矩阵、刚度矩阵和阻尼矩阵.
图 1(Fig. 1)
 | 图 1 悬臂梁单输入激励模型Fig.1 Cantilever beam model with single-input excitation |
2.1 单输入激励识别在上述悬臂梁结构中第6号节点处施加激励F6(t)=400sin(3πt).选取采样频率为20 Hz, 采样时间为5 s, 由有限元仿真计算, 测得悬臂梁第5号、第9号节点处的横向位移响应时间历程, 作为输入数据.
根据施载点和响应点的位置提取矩阵H对应的行与列, 分别基于Houbolt法与传统的状态空间法计算得到动载荷F6的时间历程识别结果, 如图 2所示.
图 2(Fig. 2)
 | 图 2 F6第6号节点激励时域识别结果Fig.2 Identification result in the time domain of excitation at No. 6 of F6 |
识别载荷与实际载荷的相对误差为
 | (15) |
由图 2的识别结果计算识别相对误差, 如图 3所示.
图 3(Fig. 3)
 | 图 3 F6第6号节点激励识别相对误差Fig.3 Relative error of excitation identification at No. 6 of F6 |
上述工况中, 其他条件不变, 在第6号节点处施加激励F′6=400sin(3πt)+300sin(1.5πt), 同理计算可得动载荷F′6时域识别结果如图 4所示, 图 5为频域识别结果.
图 4(Fig. 4)
 | 图 4 F′6第6号节点激励时域识别结果Fig.4 Identification result in the time domain of excitation at No. 6 of F′6 |
图 5(Fig. 5)
 | 图 5 F′6第6号节点激励频域识别结果Fig.5 Identification result in the frequency domain of excitation at No. 6 of F′6 (a)—本文方法辨识频谱；(b)—实际载荷频谱；(c)—状态空间法辨识频谱. |
由上述识别结果分别计算时域上F6和F′6平均相对误差和频域上F′6对应频率成分的幅值相对误差, 结果如表 1和表 2所示.
表 1(Table 1)
 表 1 时域识别结果平均相对误差Table 1 Average relative error of identification result in the time domain
| 表 1 时域识别结果平均相对误差 Table 1 Average relative error of identification result in the time domain |
表 2(Table 2)
表 2 F′6频率成分的幅值相对误差Table 2 Relative error of amplitude for the frequency components of F′6
| 表 2 F′6频率成分的幅值相对误差 Table 2 Relative error of amplitude for the frequency components of F′6 |
由图 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) |
 | 图 6 第6号节点激励时域识别结果Fig.6 Identification result in the time domain of excitation at No. 6 |
图 7(Fig. 7)
 | 图 7 第9号节点激励时域识别结果Fig.7 Identification result in the time domain of excitation at No. 9 |
图 8(Fig. 8)
 | 图 8 第6号节点激励频域识别结果Fig.8 Identification result in the frequency domain of excitation at No. 6 (a)—本文方法辨识频谱；(b)—实际载荷频谱；(c)—状态空间法辨识频谱. |
图 9(Fig. 9)
 | 图 9 第9号节点激励频域识别结果Fig.9 Identification result in the frequency domain of excitation at No. 9 (a)—本文方法辨识频谱；(b)—实际载荷频谱；(c)—状态空间法辨识频谱. |
由上述F6和F9识别结果分别计算时域上平均相对误差如表 3所示, 以及频域上对应频率成分的幅值相对误差如表 4及表 5所示.
表 3(Table 3)
表 3 多输入激励下时域识别结果相对误差Table 3 Relative error of identification result in time domain of multiple-input excitation
| 表 3 多输入激励下时域识别结果相对误差 Table 3 Relative error of identification result in time domain of multiple-input excitation |
表 4(Table 4)
表 4 F6频率成分的幅值相对误差Table 4 Relative error of amplitude for the frequency components of F6
| 表 4 F6频率成分的幅值相对误差 Table 4 Relative error of amplitude for the frequency components of F6 |
表 5(Table 5)
表 5 F9频率成分的幅值相对误差Table 5 Amplitude′s relative error of the frequency components of F9
| 表 5 F9频率成分的幅值相对误差 Table 5 Amplitude′s relative error of the frequency components of F9 |
由图 6~图 9及表 3~表 5载荷识别结果分析可知, 本文提出的方法识别得到的第6号和第9号节点处载荷时间历程更加接近实际载荷, 相对误差更小, 且在波峰波谷处相对更加稳定, 而状态空间法识别得到的载荷有较大的波动.同时在频域上, 本文提出的方法识别得到的载荷对应频率成分的幅值大小更接近于实际载荷, 没有多余频率成分或漏频现象, 而状态空间法识别得到的载荷含有多余的不可忽略的频率成分, 且对应幅值大小与实际载荷相对误差较大.
2.3 多输入含噪声激励识别上节工况中, 其他条件不变, 加入10%随机噪声干扰, 计算第6号和第9号节点处载荷识别相对误差如图 10所示, 并由此计算得到相对误差如表 6所示.
图 10(Fig. 10)
 | 图 10 多输入含噪声激励识别相对误差Fig.10 Relative error of multiple-input excitation identification with noise (a)—F6；(b)—F9. |
表 6(Table 6)
表 6 多输入含噪声激励下时域识别结果相对误差Table 6 Relative error of identification result in the time domain of multiple-input excitation with noise
| 表 6 多输入含噪声激励下时域识别结果相对误差 Table 6 Relative error of identification result in the time domain of multiple-input excitation with noise |
由图 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 |
