范玉川^1,2, 黄清云³, 鲁艳², 赵春雨¹
1. 东北大学 机械工程与自动化学院，辽宁 沈阳 110819;
2. 潍柴动力股份有限公司，山东 潍坊 261001;
3. 香港城市大学 工学院，香港 999077
收稿日期：2018-12-10
基金项目：国家自然科学基金资助项目(51775094)。
作者简介：范玉川(1988-)，男，河南辉县人，东北大学博士研究生;
赵春雨(1963-), 男, 辽宁黑山人，东北大学教授，博士生导师。

摘要：基于Newmark-β数值仿真方法，对于刚度变化的单自由度非线性体系，采用修正的Newton-Raphson迭代方法最小化由切线刚度代替变化刚度代入的误差，推导出非线性体系在已知外部激励、体系特性下的动力响应迭代求解过程，并反向推导出在已知动力响应、体系特性下动载荷的反求迭代求解过程.通过算例分析验证了应用该修正迭代方法进行非线性体系的载荷识别是可行的，克服了无迭代方法的误差累积缺点.
关键词：Newmark-β法Newton-Raphson迭代非线性载荷识别误差累积
Dynamic Load Identification of a Nonlinear System Based on Newmark-β Method
FAN Yu-chuan^1,2, HUANG Qing-yun³, LU Yan², ZHAO Chun-yu¹
1. School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China;
2. Weichai Power Company Limited, Weifang 261001, China;
3. College of Engineering, City University of Hong Kong, Hong Kong 999077, China
Corresponding author: ZHAO Chun-yu, E-mail: chyzhao@mail.neu.edu.cn
Abstract: Based on the Newmark-β numerical simulation method, a modified Newton-Raphson iteration method is used to minimize the substitution error of tangent stiffness for variable stiffness for a single-degree-of-freedom nonlinear system. Then, the iterative solution process of dynamic response of the nonlinear system under the known external excitation and system characteristics is deduced, and the inverse iterative solution process of dynamic load under the known dynamic response and system characteristics is deduced. A numerical example shows that the modified iteration method is feasible for load identification of nonlinear systems and overcomes the shortcoming of error accumulation without the iteration method.
Key words: Newmark-βNewton-Raphson iterationnonlinearityload identificationerror accumulation
动载荷辨识在机械动力结构可靠性分析、故障诊断和健康监测等许多工程研究中发挥着重要作用.当前，动态荷载识别主要包括两种方法：频域法^[1-3]和时域法^[4].频域动载荷识别的研究起步较早，理论也较为成熟.频域法根据传递函数矩阵与系统响应谱的关系确定动力谱，或在模态坐标变换后计算频域内模态力的动力特性^[5].时域法是基于载荷与响应之间的复卷积关系进行反分析，直接在时间域中检索动态载荷的时间历程.时域方法不需要傅立叶变换，结果直观，近年来对时域方法的研究也有很大的发展.Fan等^[6]基于Wilson-θ数值算法提出了一种改进的动载荷识别算法，取得了良好的识别效果.Liu等^[7]提出了一种新的时域动态伽辽金法(TDM)，利用形状函数逼近动态载荷、核函数响应和实测结构响应三个参数，然后以形状函数为加权函数建立了前向模型TDGM.Qiao等^[8]提出了动载荷识别的稀疏算法，对冲击载荷有较好的识别效果.徐婧等^[9]基于Newmark-β数值方法提出了一种动载荷识别迭代算法，取得了较高的识别精度.然而，少有学者对非线性体系的动载荷识别理论进行研究.
本文基于Newmark-β数值方法，采用修正的Newton-Raphson迭代方法使由切线刚度代替变化刚度代入的误差最小化，推导出非线性体系动载荷反求迭代过程，消除了无迭代方法的误差及其逐步累积.
1 基于Newmark-β法的动力响应求解由Newmark-β时间步进法^[10]，在时间间隔[t_i, t_i+1]内基于下面的公式:

(1a)

(1b)

式中：β=1/4；γ=1/2.
为了方便地扩展到非线性体系，使用增量

(2)

(3)

式(1)可以重新写为

(4a)

(4b)

由式(4b)可得

(5)

在式(4a)中代入式(5)，可得

(6)

对任一非弹性系统，在i时刻满足运动方程

(7)

式中：m和c分别是体系的质量和阻尼；p_i是外力；(f_s)_i是i时刻的抗力；u_i，和分别是体系的位移、速度和加速度.同时在i+1时刻满足

(8)

对于非线性体系，由式(7)和式(8)之差给出增量平衡方程：

(9)

在非线性体系中，增量抗力为

(10)

其中，(k_i)_T为切线刚度，去掉下标T，然后代入式(9)中可得

(11)

与线性体系相比，非线性反应分析要求出每个时间步开始时的切线刚度k_i.
将式(5)和式(6)代入式(11)，可得

(12)

其中：

(13)

(14)

在刚度不变或者变化不大的线性体系中，

(15)

由式(16)可得增量位移：

(16)

又由式(6)和式(5)可得Δ和Δ，则由式(2)可得u_i+1, 和.
然而，在每个时间步内用切线刚度来代替变化的刚度，会产生一定误差.如图 1中的力-变形关系所示，在时间步开始时刻的位移为图示的a点，在a点使用切线刚度，从时刻i到时刻i+1的数值积分导出位移u_i+1，标识为点b.如果能沿着正确的曲线，那么结果会是b′点对应的位移.这个偏差经过一系列时间步的累积，会引入非常大的误差.这些误差可用一种修正的Newton-Raphson迭代方法使其最小化.
图 1(Fig. 1)

图 1 非线性体系无迭代求解Fig.1 Non iterative solution of nonlinear systems

图 2详细叙述了该迭代方法.第一个迭代步是将式(12)应用于上述过程，有

(17)

图 2(Fig. 2)

图 2 非线性体系一个时间步内的迭代Fig.2 Iteration in a time step for nonlinear systems

确定Δu⁽¹⁾(相应于图 1中的b点)，作为最终Δu(相应于图 1中的b′点)的第一次近似值.与Δu⁽¹⁾相关的真实力是Δf⁽¹⁾，它比小，定义残余力为ΔR⁽²⁾=-Δf⁽¹⁾.这个残余的力所引起的附加位移Δu⁽²⁾由下式确定：

(18)

使用这个附加位移，以寻找残余力的新值，继续这个过程直到达到收敛为止.这个从i到i+1的时间步内的迭代过程称为修正的Newton-Raphson法，如图 3所示.
图 3(Fig. 3)

图 3 修正的Newton-Raphson迭代流程图Fig.3 Flow chart for modified Newton-Raphson iteration

其中，在静力分析中，Δf^(j)=f_s^(j)-f_s^(j-1)，剩余的项来自于体系的动力学，如式(14)中含有c和m的项所反映的，可将其表达为.
在l次迭代后，当增量位移Δu^(l)与当前求出的值相比变得足够小，即

(9)

迭代过程结束，于是从i到i+1时间步内的位移增量为

(20)

这是一个Δu_i的相对精确值，代替无迭代获得的值，而无迭代的值与一次迭代后获得的Δu⁽¹⁾相同.
已知Δu_i后，开始时间步进求解过程，如图 4所示.
图 4(Fig. 4)

图 4 非线性系统Newmark-β法求解流程图Fig.4 Flow chart for solving the nonlinear system based on Newmark-β method

原始的Newton-Raphson法比上面描述的修正迭代过程收敛得更快，代价是附加的计算.在每次迭代中，用切线刚度k_i^(j)代替k_i，可获得对收敛性的改进，用较少的迭代次数即可收敛.然而，它需要在每次迭代时求切线刚度，涉及附加的计算，对于多自由度体系这种附加的计算可能是大量的.因此采用修正的Newton-Raphson迭代方法能减小计算量.
2 基于Newmark-β法的载荷识别根据上述动力反应求解过程，由外力时间历程和体系特性参数，可求得非线性体系的位移、速度和加速度时间历程.这是正向求解过程.
若已知非线性体系特性参数和体系位移时间历程，现推导非线性体系外力时间历程的反求过程：
选定时间步长Δt，由体系的特性m，k_i和c，算法参数γ和β，根据已知体系位移时间历程u_i由式(2)可得其增量形式Δu_i.
由式(6)和式(5)进一步计算可得Δ和Δ.
由式(14)可得.
现选定一无限小量ε，确定一增量位移Δu⁽¹⁾, 使得与Δu相比足够小，即满足式(19).
由第l次迭代的增量位移Δu⁽¹⁾确定第l次迭代的残余力为

(21)

设第l-1次迭代的增量位移Δu^(l-1)，则由增量位移确定迭代前的位移为

(22)

(23)

由u_i+1^(l-1)和u_i+1^(l-2)可分别确定f_s^(l-1)和f_s^(l-2)，则根据修正的Newton-Raphson迭代过程计算可得Δf^(l-1)为

(24)

由图 3中的ΔR^(j+1)=ΔR^(j)-Δf^(j)可得

(25)

由图 3中的Δx^(j)=ΔR^(j)可得

(26)

将式(21)，式(24)和式(26)代入式(25)可得

(27)

上式中，只有一个未知数Δu^(l-1)，因此可解得其值，再由式(26)进一步求得ΔR^(l-1).
以此类推，一直迭代下去.依次可求得Δu^(l-2)，Δu^(l-3)，Δu^(l-4), …及ΔR^(l-2)，ΔR^(l-3)，ΔR^(l-4)…
直至经过m次迭代，所有增量位移的和无限接近Δu_i时迭代结束，即满足

(28)

此时解得的ΔR^(l-m)，即为采用迭代法求解得到的Δ.
由图 4中的可得

(29)

则由式(3)可得

(30)

由上述迭代过程，可依次求得体系外载荷时间历程.
3 算例现有一单自由度系统，各个系统参数为：质量为0.25 kg, 阻尼0.15.并且系统均为零初始条件.由于该系统具有材料非线性，其恢复力与变形关系如图 5所示.
图 5(Fig. 5)

图 5 恢复力-变形关系Fig.5 Relation of restoring force and deformation

其非线性弹性力与变形关系式可以表示为

(31)

已知该系统在施加一个动载荷
作用下前1.5 s的系统响应x(t)，现采用上述方法识别该动载荷.
选取时间步长Δt为0.1 s，选取Newmark-β法中的计算参数β=1/4, γ=1/2，并进行初始计算：

(32)

代入非线性系统动载荷识别迭代算法中，可依次求得前1.5 s的载荷时间历程，同时与采用无迭代方法求解得到的载荷识别结果进行比较，如图 6所示.
图 6(Fig. 6)

图 6 非线性系统载荷时间历程识别结果比较Fig.6 Comparison of load time history identification results for nonlinear systems

由图 6结果分析可知，在刚度发生变化的非线性系统中，基于Newmark-β逐步积分法，采用修正的Newton-Raphson方法迭代求解所得动载荷识别结果比无迭代方法求解得到的动载荷识别结果更加接近实际作用载荷，消除了无迭代方法载荷识别的误差累积，具有更好的识别效果和更小的识别误差.
4 结论1) 本文基于Newmark-β数值仿真方法，在刚度变化的非线性体系中，利用修正的Newton-Raphson迭代方法使由切线刚度代替变化刚度带来的误差最小化，推导出在已知外部激励、体系特性下的动力响应迭代求解过程，在减小计算量的同时，增加了数值求解精度.
2) 基于动力反应正向求解过程，反向推导出在已知动力响应、体系特性下动载荷的反求迭代求解过程.通过算例验证，该方法具有可行性，相比于无迭代的方法，本文基于Newmark-β法采用的修正Newton-Raphson迭代方法得到的识别结果更加接近实际动载荷，克服了无迭代方法的误差累积缺点.
参考文献

[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, Xiao 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]	O'Brien 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]	Ewins D J. Modal testing:theory, practice and application[M]. 2nd ed. Baldock: Research Studies Press, 2000.
[6]	Fan Y C, Zhao C Y, Yu H Y. Research on dynamic load identification based on explicit Wilson-θ and improved regularization algorithm[J].Shock and Vibration, 2019, 2019: 8756546.
[7]	Liu J, Meng X, Jiang X, et al. Time-domain Galerkin method for dynamic load identification[J].International Journal for Numerical Methods in Engineering, 2016, 105(8): 620–640.DOI:10.1002/nme.4991
[8]	Qiao B J, Zhang X, Wang C, et al. Sparse regularization for force identification using dictionaries[J].Journal of Sound and Vibration, 2016, 368: 71–86.DOI:10.1016/j.jsv.2016.01.030
[9]	徐菁, 张方, 姜金辉, 等. 运用数值迭代的动载荷识别算法[J].振动工程学报, 2014(5): 702–707. ( Xu Jing, Zhang Fang, Jiang Jin-hui, et al. Dynamic load identification algorithm based on numerical iteration[J].Journal of Vibration Engineering, 2014(5): 702–707.DOI:10.3969/j.issn.1004-4523.2014.05.008)
[10]	Chopra A K, Naeim F. Dynamics of structures—theory and applications to earthquake engineering[J].Earthquake Spectra, 2001, 17(3): 549–550.DOI:10.1193/1.1586188