结构动力学多自由度线性体系Wilson-θ法程序编写

结构动力学大作业班级：学号：姓名：目录1. Wilson-θ法原理简介 (2)2. Wilson-θ程序验算 (3)2.1△t的影响 (4)2.2 θ的影响 (5)3. 非线性问题求解 (5)4. 附录 (8)Wilson-θ法源程序 (8)1. Wilson -θ法原理简介图1-1Wilson-θ法示意图Wilson-θ法是基于对加速度a 的插值近似得到的，图1-1为Wilson-θ法的原理示意图。

推导由t 时刻的状态求t +△t 时刻的状态的递推公式：{}{}{}{}()t tt t t y y y y tτθτθ++∆=+-∆ (1-1)对τ积分可得速度与位移的表达式如下：{}{}{}{}{}2()2t t t t t t yy y y ytτθττθ++∆=++-∆ (1-2){}{}{}{}{}{}23()26t t t t t t t y y y y y ytτθτττθ++∆=+++-∆ (1-3)其中τ=θt ，由式(1-2)、(1-3)可以解出：{}{}{}{}{}266()2()t t t tt t t y y y y y t tθθθθ+∆+∆=---∆∆(1-4){}{}{}{}{}3()22t t t t t t t tyy y y y t θθθθ+∆+∆∆=---∆(1-5)将式(1-4)、(1-5)带入运动方程：[]{}[]{}[]{}{}m y C y k y P ++=(1-6)[]{}[]{}[]{}{}t t t t t t t tm y C y k y P θθθθ+∆+∆+∆+∆++= (1-7)注意到此时的式子为{{}t t y θ+∆}和上一个时刻{}t y 、{}t y、{}t y 以及t +θ△t 时刻的荷载{}t t P θ+∆相关，可以运用迭代的思想来求解，下图给出线弹性条件下Wilson -θ法的流程图：图1-2Wilson-θ法流程图2.Wilson-θ程序验算对线弹性条件下的Wilson-θ法进行MATLAB编程，源代码见附录。

第三章离散化结构动力方程的解法(2013.4.24）§3.1 绪言对于一个实际结构，由有限元法离散化处理后，应用瞬时最小势能原理可导出动力方程[]{}[]{}[]{}{}++=（3.1）M u C u K u F(t)这里，{}u、{}u、{}u及{}F t分别表示加速度、速度、位移及所()作用的外力矢量，他们都是与时间有关的。

从数学的角度来看，式（3.1）是一个常系数的二阶线性常微分方程组，对于它的求解原则上并无困难。

但是，由于[]M、[]C 和[]K的阶数非常高，使得式（3.1）的求解必须花费很大的代价，便促使人们去寻求一些效率高的近似计算方法。

目前，用于求解式（3.1）的方法，大致可分为两大类。

一是坐标变换法，它是对结构动力方程式（3.1），在求解之前，进行模态坐标变换，实际上就是一种Ritz变换，即把原物理空间的动力方程变换到模态空间中去求解。

现在，普遍使用的方法是模态（振型）迭加法，即用结构的前q阶实际主模态集（主振型阵）构成坐标变换阵进行变换。

通过这一变换，实现降阶，求较好的近似解，而且，还用解除耦合的办法，简化方程的计算。

还有一种所谓假设模态法，即是用一组假设模态，构成模态坐标变换阵进行变换，获得一组降阶的而不解耦的模态基坐标方程。

显然，这种方法的计算精度，取决于所假设的模态。

用Ritz矢量法求解的近似模态作为假设模态，可得到满足要求的精度。

二是直接积分法，它是对式（3.1）在求解之前，不进行坐标变换，直接进行数值积分计算。

这种方法的特点是对时域进行离散，将式（3.1）分为各离散时刻的方程，然后，将该时刻的加速度和速度用相邻时刻的各位移线性组合而成，于是，式（3.1）就化为一个由位移组成的该离散时刻上的响应值，通常又称为逐步积分法。

线性代数方程组的解法与静力时刻的位移来线性组合，就导致了各种不同的方法。

主要有中央差分法，Houbolt 方法，Wilson -θ法和Newmark 方法等。

结构动力学大作业(威尔逊- 法）姓名：学号:班级：专业：威尔逊—θ法原理及应用【摘要】在求解单自由度体系振动方程时我们用了常加速度法及线加速度法等数值分析方法。

在多自由度体系中，也有类似求解方法,即中心差分法及威尔逊—θ法。

实际上后两种方法也能求解单自由度体系振动方程。

对于数值方法，有三个重要要求：收敛性、稳定性及精度。

本文推导了威尔逊-θ法的公式，并利用MATLAB 编程来研究单自由度体系的动力特性。

【关键词】威尔逊—θ法 冲击荷载 阻尼比【正文】威尔逊-θ法可以很方便的求解任意荷载作用下单自由度体系振动问题。

实际上,当 1.37θ>时，威尔逊—θ法是无条件收敛的. 一、威尔逊—θ法的原理威尔逊-θ法是线性加速度法的一种拓展（当1θ=时,两者相同），其基本思路和实现方法是求出在时间段[],t t t θ+∆时刻的运动，其中1θ≥,然后通过内插得到i t t +∆时刻的运动（见图 1。

1)。

图 1。

11、公式推导推导由t 时刻的状态求t t θ+∆时刻的状态的递推公式:{}{}{}{})(t t t t t yy t y y -∆+=∆++θτθτ对τ积分{}{}{}{}{})(22t t t t t t yy t y y y-∆++=∆++θτθττ{}{}{}{}{}{})(6232t t t t t t t yy t y y y y -∆+++=∆++θτθτττt ∆=θτ{}{}{}{}{})(21t t t t t t t yy t y t y y -∆+∆+=∆+∆+θθθθ{}{}{}{}{})2(6)(2t t t t t tt yy t y t y y +∆+∆+=∆+∆+θθθθ {}{}{}{}{}t t t t t t t y y t y y t y26)()(62-∆--∆=∆+∆+θθθθ{}{}{}{}{}t t t t t t t yty y y t y22)(3∆---∆=∆+∆+θθθθ[]{}[]{}[]{}{}P y k y C ym =++ []{}[]{}[]{}{}t t t t t t t t P y k y C y m ∆+∆+∆+∆+=++θθθθ[]{}{}t t tt R y k ∆+∆+=θθ[][][][]c tm t k k ∆+∆+=θθ3)(62[]{}{}{}[]{}{}{}[]{}{}{})223()26)(6()(2t tt t t t t tt ty ty y t c y y t y t m P P P R ∆++∆++∆+∆+-+=∆+θθθθθ2、MA TLAB 源程序: clc;clear;K=input （’请输入结构刚度k （N/m ）')； M=input （'请输入质量（kg ）')；C=input （’请输入阻尼（N ＊s/m ）'）; t=sym （’t ’);%产生符号对象t Pt=input(’请输入荷载）;Tp=input （’请输入荷载加载时长（s)')； Tu=input （'请输入需要计算的时间长度(s ） ’）； dt=input （'请输入积分步长（s)'）; Sita=input('请输入θ’）;uds=0：dt:Tu;％确定各积分步时刻pds=0：dt：Tp;Lu=length（uds）;Lp=length（pds）;if isa(Pt,'sym'）%荷载为函数P=subs(Pt,t，uds）; %将荷载在各时间步离散if Lu〉LpP（Lp+1:Lu)=0;endelseif isnumeric（Pt）%荷载为散点if Lu〈=LpP=Pt(1:Lu)；elseP(1:Lp)=Pt;P(Lp+1：Lu）=0；endendy=zeros（1,Lu)；％给位移矩阵分配空间y1=zeros(1,Lu);％给速度矩阵分配空间y2=zeros(1，Lu)；%给加速度矩阵分配空间pp=zeros(1,Lu-1);％给广义力矩阵分配空间yy=zeros(1，Lu-1)；％给y（t+theta*t）矩阵分配FF=zeros（1，Lu);%给内力矩阵分配空间y(1)=input('请输入初位移(m)’）；y1（1）=input（’请输入初速度(m/s)');％——-—-——-———--———--初始计算-—-—------———————--——--——y2(1)=(P(1）—C＊y1(1)-K*y（1)）/M；％初始加速度FF（1）=P(1）-M＊y2(1）；l=6/（Sita*dt）^2；q=3/(Sita*dt);r=6/（Sita＊dt);s=Sita＊dt/2；for z=1:Lu—1kk=K+l＊M+q*C;pp(z)=P（z)+Sita＊（P（z+1）—P(z））+(l＊y（z）+r＊y1（z）+2＊y2（z))＊M+（q*y(z）+2＊y1（z）+s＊y2（z)）*C；yy（z)=pp（z）/kk;y2（z+1）=l/Sita＊（yy（z)—y（z)）-l＊dt＊y1(z)+（1-3/Sita）*y2（z）；y1（z+1）=y1（z）+dt/2*（y2（z+1)+y2(zp));y(z+1）=y（z）+y1(z)*dt+dt＊dt/6*(y2(z+1)+2＊y2(z））；FF(z+1）=P(z+1）—M*y2(z+1）;endplot （uds ，y ，’r ’）,xlabel('时间 t ’）,ylabel('位移 y ’），title （'位移图形’） 二、利用威尔逊-θ法求冲击荷载下的结构反应1、矩形脉冲研究不同时长脉冲作用下，体系振动位移。

中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析作者：于津津贾慧敏宋敏来源：《教育周报·教育论坛》2018年第05期摘要：在动荷载作用下的物体位移、速度和加速度的计算中，中心差分法、纽马克法和威尔逊-θ法三类方法都是可取的，为结构动力学的理论研究提供了参考。

但三类方法与精确值之间均存在一定的误差，本文基于这一问题进行研究和计算，通过图表展示这三类方法与精确值之间的关系。

关键词：结构动力学;中心差分法;纽马克法;威尔逊-θ法一、引言：结构动响应的数值计算问题，主要针对多自由或者连续体经过空间散离后建立的二阶常微分方程组形式的运动控制方程：[M] {¨x}+[ C] {﹒x}+[ K ]{x}=Q ;;（1）为了探究三种方法相较于精确解的误差，用如下具体问题进行具体分析。

如图1所示，该体系在冲击荷载 p（t）=[0 10]T 作用下，求该体系的位移反应表达式，质量单位Kg，弹簧k单位N/cm。

另：自由振动的周期T1=4.45，T2=2.8，使用中心差分法计算，取时间步长Δt=0.1，T2=0.28，并假定X0=0;V0=0试计算这个系统在前12个时间步长的反应。

取δ=0.25，γ=0.5，用纽马克法计算该系统的动力反应。

取θ=1.4，用威尔逊-θ法计算该系统的反应。

二、计算方法简介1、精确解计算根据精确解的计算公式可得：X1（t）=1-5/3×cos（2^0.5×t）+2/3×cos（5^0.5×t）X2（t）;=3-5/3×cos（2^0.5×t）-4/3×cos（5^0.5×t）速度的计算公式为位移的导数，加速度的公式为速度的导数。

（下同）2、中心差分法用位移向前一步的差分表示的速度后一步的差分表示的速度的平均来确定当前时刻的速度，得到以当前时刻t为中心的前后时刻位移的差分表示的速度，即：若：x=x0-1/（2×a1）×d x0+1/（2×a0）×d2x0; ;x1（t）=x0（1）;x2（t）=x0（2）;3、纽马克法当在t时刻的响应{x}t，{﹒x}t，{¨x}t，已知时，要求下一时刻t+Δt的响应值{x} t+Δt，{﹒x} t+Δt，{¨x} t+Δt，令在待求时刻动力学方程成立，即：{﹒x} t+Δt={﹒x}t+Δt（1-γ）{¨x} t +γΔt{¨x} t+Δt ;;（2）{x} t+Δt={x}t+{﹒x}tΔt+（0.5-δ）{¨x} tΔt^2+δ{¨x} t+ΔtΔt^2 ;（3）β，γ为按积分精度和稳定性要求而确定的参数，由式3可解得：1/{¨X}t+Δt;=βΔt 2（{X}t+Δt -{x}t）-βΔt ×1/{﹒x}t-（2β-1） ×1/{¨x}t ;;（4）将（4）带入（2）得：{﹒x}t+Δt =γ/βΔt 2×（{x}t+Δt -{x}t）+（1 –γ/β）{﹒x}t +（1 -1/2β）t{¨x}t ;;（5）{x}t +Δt 可由t +Δt 时刻的运动方程求得，即：[M]{¨X}t+Δt +[C]{¨X}t +Δt +[K]{X}t +Δt =[F] t +Δt ;;（6）將式（4）、式（5）代入式（6），可求得求得{X}t+Δt，后求{﹒X}t +Δt 和{¨X}t +Δt。

结构动力响应数值计算方法对比分析作者：李涵来源：《青年生活》2019年第21期摘要：中心差分法、纽马克法、威尔逊-法是结构动力学中常用的三种方法，为了系统的比较其优缺性，本文针对一个双自由度的体系，首先根据已知条件计算出振动微分方程，运用Matlab计算出可求出12个步长内相应的位移值，即精确解。

然后分别运用中心差分法，纽马克法，威尔逊-法求出其近似解;最后通过三种方法的近似解与精确解相对比，进而分析出三种计算方法的优缺性，为结构动力计算提供依据。

关键词：动力计算、中心差分法、纽马克法、威尔逊-法1、动力体系概况2、精确解推导针对该双自由度体系，理论推导出系统的位移表达式，通过代入各时刻周期得出位移在各时刻的具体数值，即位移精确解。

对位移方程求一阶导数得出速度方程，求二阶导数求出加速度方程。

代入各时刻的周期值，通过Matlab计算得出位移、速度、加速度的数值如下：3、三种数值计算方法3.1、中心差分法中心差分法是基于用有限差分代替位移对时间的求导，对位移一阶求导得到速度，对位移二阶求导得加速度。

通过Matlab计算出前4个步长所对应的位移响应。

3.2、纽马克法纽马克-β法是一种将线性加速度方法普遍化的方法。

通过Matlab计算出前4个步长所对应的位移响应。

3.3威尔逊-法通过Matlab计算出前4个步长所对应的位移响应。

4、近似解与精确解对比分析从上述结构的位移、速度、加速度可以看出，三种方法都能大致表示该体系大体运动趋势，并且误差较小。

其中，在描述物体位移时，中心差分法较后两种方法更为精确。

然而在描述速度和加速度时，中心差分法表现出了较大的误差，而纽马克和威尔逊法则能更详尽的表征物体速度和加速度。

5、结论中心差分法、纽马克法和威尔逊-θ法均是结构动力计算中的常用方法。

本文针对具体的计算实例，分别计算出三种方法的动力响应结果，并与精确解进行对比。

经过分析，中心差分法能更精确的表示物体位移响应，而纽马克和威尔逊法在表征物体速度和加速度方面相较于中心差分法更为精确，三种方法，各有其优缺点，应视具体情况采用相应的计算方法。

收稿日期:２０１８－０４－１９基金项目:国家自然科学基金资助项目(５１７７５０９４).作者简介:范玉川(１９８８－)ꎬ男ꎬ河南新乡人ꎬ东北大学博士研究生ꎻ赵春雨(１９６３－)ꎬ男ꎬ辽宁黑山人ꎬ东北大学教授ꎬ博士生导师ꎻ张义民(１９５８－)ꎬ男ꎬ吉林长春人ꎬ东北大学"长江学者奖励计划"特聘教授ꎬ博士生导师.第４０卷第５期２０１９年５月东北大学学报(自然科学版)ＪｏｕｒｎａｌｏｆＮｏｒｔｈｅａｓｔｅｒｎＵｎｉｖｅｒｓｉｔｙ(ＮａｔｕｒａｌＳｃｉｅｎｃｅ)Ｖｏｌ.４０ꎬＮｏ.５Ｍａｙ２０１９㊀ｄｏｉ:１０.１２０６８/ｊ.ｉｓｓｎ.１００５－３０２６.２０１９.０５.０１３基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别研究范玉川１ꎬ赵春雨１ꎬ鲁㊀艳２ꎬ张义民１(１ 东北大学机械工程与自动化学院ꎬ辽宁沈阳㊀１１０８１９ꎻ２ 郑州信大先进技术研究院ꎬ河南郑州㊀４５０００１)摘㊀㊀㊀要:推导出多自由度动力学方程的Ｗｉｌｓｏｎ－θ数值算法显式表达形式ꎬ进而提出了一种显式Ｗｉｌｓｏｎ－θ的动载荷识别算法.该算法避免了Ｗｉｌｓｏｎ－θ算法的隐式迭代形式的迭代误差ꎬ在拥有显式算法特性的同时具备隐式算法的特性.当θ取合适的值时ꎬ该算法是无条件稳定的.通过悬臂梁的算例和实验对算法的识别效果进行了验证ꎬ并与传统的状态空间法的识别结果进行了对比.结果表明:该算法不仅能够对矩形载荷㊁谐波载荷和随机载荷进行准确地识别ꎬ并且比状态空间法的识别精度更高.关㊀键㊀词:Ｗｉｌｓｏｎ－θ法ꎻ显式表达ꎻ载荷识别ꎻ无条件稳定ꎻ状态空间法中图分类号:Ｏ３２６㊀㊀㊀文献标志码:Ａ㊀㊀㊀文章编号:１００５－３０２６(２０１９)０５－０６７３－０５ＲｅｓｅａｒｃｈｏｎＤｙｎａｍｉｃＬｏａｄＩｄｅｎｔｉｆｉｃａｔｉｏｎＢａｓｅｄｏｎＥｘｐｌｉｃｉｔＷｉｌｓｏｎ￣θＭｅｔｈｏｄＦＡＮＹｕ￣ｃｈｕａｎ１ꎬＺＨＡＯＣｈｕｎ￣ｙｕ１ꎬＬＵＹａｎ２ꎬＺＨＡＮＧＹｉ￣ｍｉｎ１(１ ＳｃｈｏｏｌｏｆＭｅｃｈａｎｉｃａｌ＆ＡｕｔｏｍａｔｉｏｎꎬＮｏｒｔｈｅａｓｔｅｒｎＵｎｉｖｅｒｓｉｔｙꎬＳｈｅｎｙａｎｇ１１０８１９ꎬＣｈｉｎａꎻ２ ＺｈｅｎｇｚｈｏｕＸｉｎｄａＩｎｓｔｉｔｕｔｅｏｆＡｄｖａｎｃｅｄＴｅｃｈｎｏｌｏｇｙꎬＺｈｅｎｇｚｈｏｕ４５０００１ꎬＣｈｉｎａ.Ｃｏｒｒｅｓｐｏｎｄｉｎｇａｕｔｈｏｒ:ＺＨＡＯＣｈｕｎ￣ｙｕꎬＥ￣ｍａｉｌ:ｃｈｙｚｈａｏ＠ｍａｉｌ.ｎｅｕ.ｅｄｕ.ｃｎ)Ａｂｓｔｒａｃｔ:ＴｈｅｅｘｐｌｉｃｉｔｅｘｐｒｅｓｓｉｏｎｏｆＷｉｌｓｏｎ￣θｎｕｍｅｒｉｃａｌａｌｇｏｒｉｔｈｍｆｏｒｍｕｌｔｉ￣ｄｏｆｓ(ｄｅｇｒｅｅｏｆｆｒｅｅｄｏｍｓ)ｄｙｎａｍｉｃｅｑｕａｔｉｏｎｉｓｄｅｒｉｖｅｄꎬａｓｗｅｌｌａｓａｎｅｘｐｌｉｃｉｔＷｉｌｓｏｎ￣θｄｙｎａｍｉｃｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎａｌｇｏｒｉｔｈｍｉｓｐｒｏｐｏｓｅｄꎬａｖｏｉｄｉｎｇｔｈｅｉｔｅｒａｔｉｏｎｅｒｒｏｒｗｈｉｌｅｋｅｅｐｉｎｇｔｈｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｔｈｅｉｍｐｌｉｃｉｔｉｔｅｒａｔｉｏｎａｌｇｏｒｉｔｈｍ.Ｔｈｅａｌｇｏｒｉｔｈｍｉｓｕｎｃｏｎｄｉｔｉｏｎａｌｌｙｓｔａｂｌｅｗｈｅｎａｐｐｌｙｉｎｇａｐｐｒｏｐｒｉａｔｅθｖａｌｕｅ.Ｔｈｅｒｅｃｏｇｎｉｔｉｏｎｅｆｆｅｃｔｏｆｔｈｅａｌｇｏｒｉｔｈｍｉｓｖｅｒｉｆｉｅｄｂｙａｎｅｘａｍｐｌｅａｎｄａｎｅｘｐｅｒｉｍｅｎｔｏｆａｃａｎｔｉｌｅｖｅｒｂｅａｍꎬａｎｄｔｈｅｒｅｓｕｌｔｓｗｅｒｅｃｏｍｐａｒｅｄｗｉｔｈｔｈｏｓｅｆｒｏｍｔｈｅｔｒａｄｉｔｉｏｎａｌｓｔａｔｅｓｐａｃｅｍｅｔｈｏｄ.Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅａｌｇｏｒｉｔｈｍｎｏｔｏｎｌｙｃａｎａｃｃｕｒａｔｅｌｙｉｄｅｎｔｉｆｙｔｈｅｒｅｃｔａｎｇｕｌａｒｌｏａｄꎬｔｈｅｈａｒｍｏｎｉｃｌｏａｄａｎｄｔｈｅｒａｎｄｏｍｌｏａｄꎬｂｕｔａｌｓｏｈａｓｈｉｇｈｅｒｒｅｃｏｇｎｉｔｉｏｎａｃｃｕｒａｃｙｔｈａｎｓｔａｔｅｓｐａｃｅｍｅｔｈｏｄ.Ｋｅｙｗｏｒｄｓ:Ｗｉｌｓｏｎ￣θꎻｅｘｐｌｉｃｉｔｆｏｒｍｕｌａｒꎻｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎꎻｕｎｃｏｎｄｉｔｉｏｎａｌｌｙｓｔａｂｌｅꎻｓｔａｔｅ￣ｓｐａｃｅｍｅｔｈｏｄ㊀㊀在机械系统设计过程中ꎬ动载荷是机械结构进行疲劳分析以及进行可靠性计算的基本依据.但是ꎬ在动态系统中ꎬ有些力是很难直接测量的ꎬ特别是结构系统内部各部件之间的相互作用力ꎬ在难以直接测量的情况下ꎬ就需要通过逆动力学的分析技术来得到这些力ꎬ因而开展动载荷识别技术的研究是很有必要的.动载荷识别方法主要分为频域法[１－２]和时域法[３]两大类.由于时域法能够识别各种类型的载荷ꎬ识别精度高ꎬ并且其识别结果具有明确的物理意义ꎬ因而ꎬ时域法越来越受到专家学者们的青睐.Ｌｉｕ等[４]推导了一种基于Ｎｅｗｍａｒｋ－β法的动载荷识别方法ꎬ将传统的隐式Ｎｅｗｍａｒｋ－β算法转化为Ａｘ＝ｂ方程解的显式形式ꎬ与隐式的Ｎｅｗｍａｒｋ－β法具有相同的特点ꎬ并与状态空间法进行了对比ꎬ表明该算法具有明显的优势.Ｌｉ等[５]提出了一种基于二阶泰勒级数展开的时域动态结构载荷识别方法ꎬ该算法将响应表示为一㊀㊀种泰勒级数的逼近形式ꎬ推导出一系列公式ꎬ并建立了系统响应㊁系统特性和输入激励相结合的显式离散方程.Ｌｉｕ等[６]提出了一种新的时域动态Ｇａｌｅｒｋｉｎ算法ꎬ将形状函数作为加权函数ꎬ建立了前向模型ＴＤＧＭꎬ与传统的格林函数法相比ꎬＴＤＧＭ能有效地克服噪声的影响ꎬ提高动态载荷识别的精度.在求逆运算中ꎬ载荷识别结果通常对结构模型中的响应和误差测量中的噪声非常敏感[７].张方等[８－９]对复杂结构的载荷识别进行了研究ꎬ推导了一种基于广义正交多项式特征技术的动载荷识别模型ꎬ可以在一定精度范围内通过有限的测量点信息对无限未知量的分布动载荷进行识别.Ａｌｌｅｎ等[１０]采用一种ＳＷＡＴ方法ꎬ可以同时识别出冲击型载荷和稳态动态载荷.本文将对Ｗｉｌｓｏｎ－θ法进行变换ꎬ推导其显式表达形式ꎬ提出一种基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别算法ꎬ并通过算例和实验验证了该算法的有效性.１㊀动载荷识别算法线性阻尼结构中ꎬ多自由度结构的动力学方程可以表示为Ｍ㊆ｘｔ＋Ｃ̇ｘｔ＋Ｋｘｔ＝Ｐｔ.(１)其中:ＭꎬＣ和Ｋ分别表示质量矩阵㊁阻尼矩阵以及刚度矩阵ꎻＰｔ是作用在结构上的外加载荷ꎻ㊆ｘｔꎬ̇ｘｔ和ｘｔ分别表示加速度响应㊁速度响应以及位移响应.本文假定阻尼为瑞利阻尼:Ｃ＝α１Ｍ＋α２Ｋ.(２)其中ꎬα１和α２表示阻尼系数.１ １㊀Ｗｉｌｓｏｎ－θ法的显式表达Ｗｉｌｓｏｎ－θ法假定在[ｔꎬｔ＋θΔｔ](θȡ１)的时间间隔内ꎬ加速度呈线性变化ꎬ如图１所示.令δ为自ｔ时刻开始的时间变量ꎬ适用于０ɤδɤθΔｔꎬ由线性加速度的假设可知ꎬ在适用范围内的加速度为㊆ｘｔ＋δ＝㊆ｘｔ＋δθΔｔ(㊆ｘｔ＋θΔｔ－㊆ｘｔ).(３)积分后可得̇ｘｔ＋δ＝̇ｘｔ＋㊆ｘｔδ＋δ２２θΔｔ(㊆ｘｔ＋θΔｔ－㊆ｘｔ)ꎬ(４)ｘｔ＋δ＝ｘｔ＋̇ｘｔδ＋１２㊆ｘｔδ２＋δ３６θΔｔ(㊆ｘｔ＋θΔｔ－㊆ｘｔ).(５)若δ＝θΔｔꎬ由式(４)和式(５)可得ｔ＋θΔｔ瞬时的速度和位移:̇ｘｔ＋θΔｔ＝̇ｘｔ＋θΔｔ２(㊆ｘｔ＋θΔｔ＋㊆ｘｔ)ꎬ(６)ｘｔ＋θΔｔ＝ｘｔ＋θΔｔ̇ｘｔ＋θ２Δｔ２６(㊆ｘｔ＋θΔｔ＋２㊆ｘｔ).(７)图１㊀Ｗｉｌｓｏｎ－θ法模型Ｆｉｇ １㊀ＭｏｄｅｌｏｆＷｉｌｓｏｎ￣θ联立式(６)和式(７)ꎬ可以用ｔ＋θΔｔ时刻的位移表示ｔ＋θΔｔ时刻的加速度和速度ꎬ即㊆ｘｔ＋θΔｔ＝６θ２Δｔ２(ｘｔ＋θΔｔ－ｘｔ)－６θΔｔ̇ｘｔ－２㊆ｘｔꎬ(８)̇ｘｔ＋θΔｔ＝３θΔｔ(ｘｔ＋θΔｔ－ｘｔ)－２̇ｘｔ－θΔｔ２㊆ｘｔ.(９)则ｔ＋θΔｔ时刻的动力方程可以表示为Ｍ㊆ｘｔ＋θΔｔ＋Ｃ̇ｘｔ＋θΔｔ＋Ｋｘｔ＋θΔｔ＝Ｐｔ＋θΔｔ.(１０)式中:Ｐｔ＋θΔｔ＝Ｐｔ＋θ(Ｐｔ＋Δｔ－Ｐｔ).(１１)将式(８)和式(９)以及式(１１)代入式(１０)ꎬ即得关于ｘｔ＋θΔｔ的求解方程为ｘｔ＋θΔｔ＝＾Ｋ－１＾Ｐｔ＋θΔｔ.(１２)式中:＾Ｋ＝Ｋ＋３θΔｔＣ＋６θ２Δｔ２Ｍꎬ(１３)＾Ｐｔ＋θΔｔ＝Ｐｔ＋θ(Ｐｔ＋Δｔ－Ｐｔ)＋Ｍ(６θ２Δｔ２ｘｔ＋６θΔｔ̇ｘｔ＋２㊆ｘｔ)＋Ｃ(３θΔｔｘｔ＋２̇ｘｔ＋θΔｔ２㊆ｘｔ).(１４)求解方程式(１２)ꎬ得ｘｔ＋θΔｔ.再把ｘｔ＋θΔｔ代入式(８)就可获得㊆ｘｔ＋θΔｔ.在式(３)中取δ＝Δｔꎬ同时将式(８)代入ꎬ可得㊆ｘｔ＋Δｔ＝６θ３Δｔ２(ｘｔ＋θΔｔ－ｘｔ)－６θ２Δｔ̇ｘｔ＋(１－３θ)㊆ｘｔ.(１５)取δ＝Δｔꎬ将式(３)分别代入式(４)和式(５)ꎬ有̇ｘｔ＋Δｔ＝̇ｘｔ＋Δｔ２(㊆ｘｔ＋Δｔ＋㊆ｘｔ)ꎬ(１６)ｘｔ＋Δｔ＝ｘｔ＋Δｔ̇ｘｔ＋Δｔ２６(㊆ｘｔ＋Δｔ＋２㊆ｘｔ).(１７)由式(１２)~式(１４)可得４７６东北大学学报(自然科学版)㊀㊀㊀第４０卷㊀㊀ｘｔ＋θΔｔ＝＾Ｋ－１(１－θ)Ｐｔ＋θ＾Ｋ－１Ｐｔ＋Δｔ＋(６θ２Δｔ２Ｍ＾Ｋ－１＋３θΔｔＣ＾Ｋ－１)ｘｔ＋(６θΔｔＭ＾Ｋ－１＋２Ｃ＾Ｋ－１)̇ｘｔ＋(２Ｍ＾Ｋ－１＋θΔｔ２Ｃ＾Ｋ－１)㊆ｘｔ.(１８)将式(１８)代入式(１５)ꎬ已知恒等式:Ｉ＝＾Ｋ－１＾Ｋꎬ可得㊆ｘｔ＋Δｔ＝Ｃ０Ｐｔ＋Ｃ１Ｐｔ＋Δｔ＋Ｃｄｘｔ＋Ｃｖ̇ｘｔ＋Ｃａ㊆ｘｔ.(１９)其中:Ｃ０＝６θ３Δｔ２＾Ｋ－１(１－θ)ꎻＣ１＝６θ２Δｔ２＾Ｋ－１ꎻＣｄ＝－６θ３Δｔ２＾Ｋ－１ＫꎻＣｖ＝６θ２Δｔ(６θ２Δｔ２Ｍ＾Ｋ－１＋２θΔｔＣ＾Ｋ－１－１)ꎻＣａ＝１２θ３Δｔ２Ｍ＾Ｋ－１＋３θ２ΔｔＣ＾Ｋ－１＋１－３θ.将式(１９)代入式(１６)可得̇ｘｔ＋Δｔ＝Ｂ０Ｐｔ＋Ｂ１Ｐｔ＋Δｔ＋Ｂｄｘｔ＋Ｂｖ̇ｘｔ＋Ｂａ㊆ｘｔ.(２０)其中:Ｂ０＝Δｔ２Ｃ０ꎻＢ１＝Δｔ２Ｃ１ꎻＢｄ＝Δｔ２ＣｄꎻＢｖ＝Δｔ２Ｃｖ＋１ꎻＢａ＝Δｔ２(Ｃａ＋１).将式(１９)代入式(１７)可得ｘｔ＋Δｔ＝Ａ０Ｐｔ＋Ａ１Ｐｔ＋Δｔ＋Ａｄｘｔ＋Ａｖ̇ｘｔ＋Ａａ㊆ｘｔ.(２１)其中:Ａ０＝Δｔ２６Ｃ０ꎻＡ１＝Δｔ２６Ｃ１ꎻＡｄ＝Δｔ２６Ｃｄ＋１ꎻＡｖ＝(Δｔ６Ｃｖ＋１)ΔｔꎻＡａ＝Δｔ２６(Ｃａ＋２).在此用ｘｉꎬ̇ｘｉꎬ㊆ｘｉ和Ｐｉ表示第ｉ时刻式(１９)~式(２１)中的位移㊁速度㊁加速度和激励ꎬ并将其表示为矩阵形式:ｘｉ＋１̇ｘｉ＋１㊆ｘｉ＋１éëêêêùûúúú＝Ａ０㊀Ａ１Ｂ０㊀Ｂ１Ｃ０㊀Ｃ１éëêêêùûúúúＰｉＰｉ＋１éëêêùûúú＋Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀Ｃａéëêêêùûúúúｘｉ̇ｘｉ㊆ｘｉéëêêêùûúúú.(２２)在第ｉ时刻的响应可以表示为ｘｉ̇ｘｉ㊆ｘｉéëêêêùûúúú＝ðｉ－１ｊ＝０ＡｄＡｖＡａＢｄＢｖＢａＣｄＣｖＣａéëêêêùûúúúｊＡ０Ａ１Ｂ０Ｂ１Ｃ０Ｃ１éëêêêùûúúúＰｉ－ｊ－１Ｐｉ－ｊ[]＋ＡｄＡｖＡａＢｄＢｖＢａＣｄＣｖＣａéëêêêùûúúúｉｘ０̇ｘ０㊆ｘ０éëêêêùûúúú.(２３)式中ꎬ两个指数ｉ和ｊ分别表示相应矩阵的幂.式(２３)是一种新型的Ｗｉｌｓｏｎ－θ法的显式表达ꎬ每一个时间步中位移㊁速度㊁加速度响应可被同时求解出来ꎬ而通常的方法是利用每个时间步的迭代算法来计算ꎬ显然这种显式的算法更有优势.１ ２㊀基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别令:ｙｉ＝ｘｉ̇ｘｉ㊆ｘｉéëêêêùûúúú－Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀Ｃａéëêêêùûúúúｉｘ０̇ｘ０㊆ｘ０éëêêêùûúúú.则式(２３)可以改写为ｙｉ＝ðｉ－１ｊ＝０Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀ＣａéëêêêùûúúúｊＡ０㊀Ａ１Ｂ０㊀Ｂ１Ｃ０㊀Ｃ１éëêêêùûúúúＰｉ－ｊ－１Ｐｉ－ｊ[].(２４)令:Ｈｋ＝Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀Ｃａéëêêêùûúúú(ｋ－１)Ａ０㊀Ａ１Ｂ０㊀Ｂ１Ｃ０㊀Ｃ１éëêêêùûúúúꎬ(２５)Ｒｉ＝Ｐｉ－１Ｐｉéëêêùûúú.(２６)式(２９)可以写成从时间段１到ｎｔ上的矩阵的卷积形式:Ｙ＝ＨＦ.(２７)其中:Ｙ＝[ｙＴ１ꎬｙＴ２ꎬ ꎬｙＴｎｔ]ＴꎻＨ＝Ｈ１０ ０Ｈ２Ｈ１ ０⋮⋮⋮ＨｎｔＨｎｔ－１ Ｈ１éëêêêêêùûúúúúúꎻＦ＝Ｒ１Ｒ２⋮Ｒｎｔéëêêêêêùûúúúúú.式(２７)可以改写成载荷识别的表达式Ｆ＝Ｈ－１Ｙ.(２８)对于一个给定的系统ꎬＨ是常数ꎬＹ可以从系统测量的响应中得到ꎬ考虑到式(２８)可能存在不适定性ꎬ可以通过Ｔｉｋｈｏｎｏｖ正则化方法来确定Ｆꎬ在阻尼最小二乘意义下对目标函数进行优化.Ｊ(Ｆꎬλ)＝ Ｙ－ＨＲ ２＋λ Ｆ ２.(２９)其中ꎬλ是正则化参数ꎬ它的值可以通过Ｌ曲线法来确定.基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别方法可以归纳为以下步骤:１)建立系统的有限元模型ꎬ确定系统的质量矩阵Ｍ㊁刚度矩阵Ｋ和阻尼矩阵Ｃꎻ２)选取时间步Δｔ和Ｗｉｌｓｏｎ－θ法的参数θꎬ参数θ的值可以取为１ ４ꎬ从而保证算法的无条件稳定(在Ｗｉｌｓｏｎ－θ法中ꎬ只要θ的值大于５７６第５期㊀㊀㊀范玉川等:基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别研究㊀㊀１ ３７ꎬ那么该算法就无条件稳定ꎬ但是θ取的过大ꎬ截断误差增大ꎬ精度会下降).３)为递推式(２２)计算矩阵Ａ１ꎬＡ２ꎬＡｄꎬＡｖꎬＡａꎬＢ１ꎬＢ２ꎬＢｄꎬＢｖꎬＢａꎬＣ１ꎬＣ２ꎬＣｄꎬＣｖ和Ｃａꎻ４)计算式(２５)中的矩阵Ｈｋ(ｋ＝０ꎬ ꎬｎｔ－１)和式(２８)中的组合矩阵Ｈꎻ５)通过实验或仿真计算获得系统的响应数据㊆ｘꎬ̇ｘ和ｘꎻ６)利用式(２９)中的Ｔｉｋｈｏｎｏｖ正则化方法来确定Ｆ.２㊀仿真算例及抗噪性能研究如图２所示ꎬ一矩形截面悬臂梁长０ ６４ｍꎬ截面尺寸为:宽０ ０５６ｍ㊁高０ ００８ｍꎬ弹性模量Ｅ＝２００ＧＰａꎬ材料密度为７８４０ｋｇ/ｍ３ꎬ该悬臂梁被划分为１８个单元ꎬ共１９个节点.图２㊀悬臂梁结构图Ｆｉｇ ２㊀Ｓｔｒｕｃｔｕｒｅｏｆｃａｎｔｉｌｅｖｅｒｂｅａｍ在第１０个节点上施加载荷ꎬ为了尽可能全面地验证本文载荷识别算法的准确性ꎬ分别施加方波载荷㊁谐波载荷和随机载荷ꎬ考虑到噪声的影响ꎬ在响应数据中加入５％的随机噪声ꎬ选取第１２个节点的响应信息进行载荷识别ꎬ时间间隔Δｔ＝０ ００１ｓ.选取第１２个节点上的数据进行载荷识别计算ꎬ本文算法与状态空间法进行对比ꎬ如图３所示ꎬｐ２为方波载荷㊁ｐ３为谐波载荷㊁ｐ４为随机载荷.两种方法对于这三种不同类型的载荷都能够比较准确的识别ꎬ但是状态空间法的识别效果明显要差一些ꎬ而基于显式Ｗｉｌｓｏｎ－θ法的载荷识别方法识别得很好ꎬ对实际载荷的还原度较高.为便于量化分析基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别算法的识别精度ꎬ可引入式(３０)计算误差:Ｅｒｒｏｒ＝Ｆｉｄ－ＦｒｅａｌＦｒｅａｌˑ１００％.(３０)式中:Ｅｒｒｏｒ表示相对识别误差ꎻＦｉｄ表示识别得到的外激励ꎻＦｒｅａｌ表示真实的外激励.运用不同的方法ꎬ在不同的采样频率和采样时间下研究动载荷识别的误差ꎬ结果列于表１中.图３㊀载荷识别结果Ｆｉｇ ３㊀Ｒｅｓｕｌｔｓｏｆｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎ(ａ) ｐ２ꎻ(ｂ) ｐ３ꎻ(ｃ) ｐ４.表１㊀两种载荷识别方法的识别误差对比Ｔａｂｌｅ１㊀Ｃｏｍｐａｒｉｓｏｎｏｆｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｅｒｒｏｒｓｗｉｔｈｄｉｆｆｅｒｅｎｔｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｍｅｔｈｏｄｓ％载荷Ｗｉｌｓｏｎ－θ法状态空间法ｐ２３ ２１５ ０４Ｐ３０ ６１３ ８４Ｐ４３ ３０５ １６㊀㊀数据表明ꎬ本文提出的基于显式Ｗｉｌｓｏｎ－θ法的载荷识别方法的识别精度明显高于基于状态空间法的载荷识别方法的识别精度ꎬ并且两种算法都显示谐波载荷的识别精度要远远好于方波载荷和随机载荷的识别精度.６７６东北大学学报(自然科学版)㊀㊀㊀第４０卷㊀㊀㊀㊀图４显示了随着噪声的增大载荷识别精度的变化情况ꎬ两种方法的识别误差都是随着噪声的增大而增大ꎬ但本文算法随噪声的增大识别误差变化较小ꎬ表明本文算法抗噪性能更好㊁识别误差更小.图４㊀载荷识别误差随噪声变化情况Ｆｉｇ ４㊀Ｖａｒｉａｔｉｏｎｏｆｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｅｒｒｏｒｗｉｔｈｎｏｉｓｅ３㊀实验验证实验用的悬臂梁模型参数与仿真悬臂梁模型参数一致ꎬ选用ＹＥ６２５１振动力学实验系统.激励载荷为谐波载荷ꎬ激励频率设置为１０Ｈｚꎬ算法中采用有限元法进行建模ꎬ代入悬臂梁的相关参数ꎬ选用测量得到的第１１个节点上的响应数据进行载荷识别ꎬ识别结果如图５所示.由图可知ꎬ该算法对实验数据的识别较准确ꎬ识别的载荷曲线与实际的载荷曲线基本吻合.图５㊀实验识别结果Ｆｉｇ ５㊀Ｒｅｓｕｌｔｓｏｆｅｘｐｅｒｉｍｅｎｔ由于实际工程中响应数据的测量位置往往是有限的ꎬ不能够随意选取ꎬ因而在本实验中分别选取第４㊁第７㊁第１０㊁第１３和第１６个节点的测量数据进行载荷的识别ꎬ同样利用式(３０)进行识别误差的计算ꎬ计算结果分别为５ ６７％ꎬ５ ７９％ꎬ５ ３８％ꎬ５ ８６％ꎬ５ ９８％.从结果可以看到ꎬ选取不同的测量点ꎬ识别得到的载荷与实际载荷的误差没有太大差别ꎬ也就是说测量点的选取对识别结果影响不大ꎻ该算法能够适应实际工程的应用环境.４㊀结㊀㊀论１)利用隐式Ｗｉｌｓｏｎ－θ法ꎬ推导出了它的一种显式表达形式ꎬ进而提出了一种基于显式Ｗｉｌｓｏｎ－θ法的载荷识别算法ꎬ该算法既具有显式算法的优势ꎬ同时具有隐式算法的特点.２)通过算例验证了本文算法在识别矩形波载荷㊁谐波载荷和随机载荷方面的识别精度都高于状态空间法的识别精度ꎬ并且在算例中加入了５％的随机噪声.实验结果表明:算法不仅能够对矩形载荷㊁谐波载荷和随机载荷准确识别ꎬ并且识别精度比状态空间法的识别精度更高.参考文献:[１]㊀ＪｉｅＨꎬＺｈａｎｇＸ.Ａｎｏｐｔｉｍｉｚａｔｉｏｎｍｅｔｈｏｄｏｆｌｏａｄｉｄｅｎ￣ｔｉｆｉｃａｔｉｏｎｉｎｆｒｅｑｕｅｎｃｙｄｏｍａｉｎ[Ｊ].Ｎｏｉｓｅ＆ＶｉｂｒａｔｉｏｎＣｏｎｔｒｏｌꎬ２００９ꎬ２９(６):３４－３６[２]㊀ＨｅＺＣꎬＬｉｎＸＹꎬＬｉＥ.Ａｎｏｖｅｌｍｅｔｈｏｄｆｏｒｌｏａｄｂｏｕｎｄｓｉｄｅｎｔｉｆｉｃａｔｉｏｎｆｏｒｕｎｃｅｒｔａｉｎｓｔｒｕｃｔｕｒｅｓｉｎｆｒｅｑｕｅｎｃｙｄｏｍａｉｎ[Ｊ].ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＭｅｔｈｏｄｓꎬ２０１７(３):１８５００５１.[３]㊀ＬａｗＳＳꎬＣｈａｎＴＨＴꎬＺｈｕＱＸ.Ｒｅｇｕｌａｒｉｚａｔｉｏｎｉｎｍｏｖｉｎｇｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＥｎｇｉｎｅｅｒｉｎｇＭｅｃｈａｎｉｃｓＡＳＣＥꎬ２００１ꎬ１２７(２):１３６－１４８.[４]㊀ＬｉｕＫꎬＬａｗＳＳꎬＺｈｕＸＱꎬｅｔａｌ.Ｅｘｐｌｉｃｉｔｆｏｒｍｏｆａｎｉｍｐｌｙ￣ＣＩＴｍｅｔｈｏｄｆｏｒｉｎｖｅｒｓｅｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＳｏｕｎｄａｎｄＶｉｂｒａｔｉｏｎꎬ２０１４ꎬ３３:７３０－７４４.[５]㊀ＬｉＸＷꎬＤｅｎｇＺＭ.Ｉｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｄｙｎａｍｉｃｌｏａｄｓｂａｓｅｄｏｎｓｅｃｏｎｄ￣ｏｒｄｅｒｔａｙｌｏｒ￣ｓｅｒｉｅｓｅｘｐａｎｓｉｏｎｍｅｔｈｏｄ[Ｊ].ＳｈｏｃｋａｎｄＶｉｂｒａｔｉｏｎꎬ２０１６(２０１６):１－９.[６]㊀ＬｉｕＪꎬＭｅｎｇＸꎬＪｉａｎｇＣꎬｅｔａｌ.Ｔｉｍｅ￣ｄｏｍａｉｎＧａｌｅｒｋｉｎｍｅｔｈｏｄｆｏｒｄｙｎａｍｉｃｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｆｏｒＭｕｍｅｒｉｃａｌＭｅｔｈｏｄｓｉｎＥｎｇｉｎｅｅｒｉｎｇꎬ２０１６ꎬ１０５:６２０－６４０.[７]㊀ＱｉａｎＢꎬＺｈａｎｇＸꎬＷａｎｇＣꎬｅｔａｌ.Ｓｐａｒｓｅｒｅｇｕｌａｒｉｚａｔｉｏｎｆｏｒｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎｕｓｉｎｇｄｉｃｔｉｏｎａｒｉｅｓ[Ｊ].ＪｏｕｒｎａｌｏｆＳｏｕｎｄａｎｄＶｉｂｒａｔｉｏｎꎬ２０１６ꎬ３６８:７１－８６.[８]㊀张方ꎬ秦远田ꎬ邓吉宏.复杂分布动载荷识别技术研究[Ｊ].振动工程学报ꎬ２００６ꎬ１９(１):８１－８５.(ＺｈａｎｇＦａｎｇꎬＱｉｎＹｕａｎ￣ｔｉａｎꎬＤｅｎｇＪｉ￣ｈｏｎｇ.Ｒｅｓｅａｒｃｈｏｆｉｄｅｎｔｉｆｉｃａｔｉｏｎｔｅｃｈｎｏｌｏｇｙｏｆｄｙｎａｍｉｃｌｏａｄｄｉｓｔｒｉｂｕｔｅｄｏｎｔｈｅｓｔｒｕｃｔｕｒｅ[Ｊ].ＪｏｕｒｎａｌｏｆＶｉｂｒａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇꎬ２００６ꎬ１９(１):８１－８５.)[９]㊀徐菁ꎬ张方ꎬ姜金辉ꎬ等.运用数值迭代的动载荷识别算法[Ｊ]ꎬ振动工程学报ꎬ２０１４ꎬ２７(１５):７０２－７０７.(ＸｕＪｉｎｇꎬＺｈａｎｇＦａｎｇꎬＪｉａｎｇＪｉｎ￣ｈｕｉꎬｅｔａｌ.Ａｎａｌｇｏｒｉｔｈｍｏｆｄｙｎａｍｉｃｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｂａｓｅｄｏｎｎｕｍｅｒｉｃａｌｉｔｅｒａｔｉｏｎꎬ[Ｊ].ＪｏｕｒｎａｌｏｆＶｉｂｒａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇꎬ２０１４ꎬ２７(１５):７０２－７０７.)[１０]ＡｌｌｅｎＭＳꎬＣａｒｎｅＴＧ.Ｄｅｌａｙｅｄｍｕｌｔｉ￣ｓｔｅｐｉｎｖｅｒｓｅｓｔｒｕｃｔｕｒａｌｆｉｌｔｅｒｆｏｒｒｏｂｕｓｔｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＭｅｃｈａｎｉｃａｌＳｙｓｔｅｍｓａｎｄＳｉｇｎａｌＰｒｏｃｅｓｓｉｎｇꎬ２００８ꎬ２２(５):１０３６－１０５４.７７６第５期㊀㊀㊀范玉川等:基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别研究。