结构动力学多自由度线性体系Wilson-θ法程序编写
结构动力学大作业2
结构动力学大作业班级：学号：姓名：目录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编程，源代码见附录。
wilson法和newmark法的理论过程
第三章离散化结构动力方程的解法(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、矩形脉冲研究不同时长脉冲作用下，体系振动位移。
威尔逊—θ法在matlab中的实现
^
^
k u(ti +θ∆t) = P(ti +θ∆t)
其中
^
k=k+
6
m+ 3 c
(θ∆t)2 θ∆t
（10） （11）
•
^
P(ti
+ θ∆t )
=
Pi
+θ
( Pi +1
−
Pi
)
+
6 [ (θ∆t ) 2
ui
+
6 θ∆t
•
ui
+
••
2 ui
]m
+
(3 θ∆t
ui
+
2 ui
+
θ∆t 2
••
ui
)c
（12）
（2）
•
τ 2 ••
τ 3 ••
••
u(ti +τ ) = u(ti ) +τ u(ti ) + 2 u(ti ) + 6θ∆t [u(ti +θ∆t) − u(ti )]
（3）
当 τ = θ∆t 时，由式（2）和式（3）得到
•
u(ti
+ θ∆t )
=
•
u(ti
)
+ θ∆t
••
u(ti
)
+
θ∆t 2
加速度 0.00E+00 -7.41E-04 -1.21E-03 -1.42E-03 -1.58E-03 -1.54E-03 -1.44E-03 -1.18E-03 -8.64E-04 -4.67E-04 -2.70E-05 3.80E-04 8.26E-04 1.10E-03 2.48E-03 7.78E-03 1.49E-02 8.72E-03 -1.65E-03 -6.50E-04 -1.54E-03 -2.86E-03 -4.51E-03 -6.70E-03 -8.43E-03 -3.74E-03 1.78E-03
中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析
中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析作者：于津津贾慧敏宋敏来源：《教育周报·教育论坛》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。
第十四章结构动力学多自由度
第 j 阶振型的广义质量 第 j 阶振型的广义刚度
j
K
* j
M
* j
32 / 67
§14-6 多自由度结构的自由振动
5、多自由度体系振型正交性应用举例
例题1 检验框架结构振型的正确性
m2
A1
1 1.618
A2
1 0.618
h2
m1
M
m
0
0 m
K
12EI l3
2 1
1
1
h1
A1TM A2 0
1 2 n
K 2M A 0
M
1
2
I
A
0
A1,A2, ,An
25 / 67
§14-6 多自由度结构的自由振动
注意 自振频率与主振型一一对应 振型只表明振动的形状，不能唯一确定其幅值 振型是多自由度特有的概念
3、振型的标准化 补充条件，使主振型用确定的幅值表示 方法1：规定振型中某元素为1，其它元素就是 相对于它的比值。通常选第一个元素或最大一 个元素，令其等于1。
k21 L kn1
k22 L kn2
L L L
k1n y1 FP1
k2n L knn
yM2 yn
FMP2 FPn
简写为 M Y&& Y& K Y FP
其中[M]为质量矩阵，在集中质点的结构中它是对角矩阵；[β] 为阻尼矩阵；[K]为刚度矩阵，根据反力互等定理，它是对称矩 阵。
一、刚度法
取任意质量mi列动力平衡方程：
FIi FRi Fei FPi 0 n
根据叠加原理有
Fei kij y j
jn1
FRi ij yj
j 1
结构动力学多自由度
▪ 振型方程：
(K i2M)ji 0 (i 1, 2, 3, n)
▪∵
K 2i M 0
▪ ∴ 第i 个振型方程中的n 个方程中只有n-1个是独立的！ ▪ ——无法得到j1i、 j2i、 … 、 jni 的确定值， ▪ 但可以确定各质点振幅之间的相对比值： ▪ —— 振型的幅值是任意的，但形状是惟一的。
一致质量矩阵：
L
pava m13v1 0 fI ( x)v( x)dx
L
0
m( x) 3( x)v3
L
1( x)v1dx
mij 0 EI ( x)i ( x) j ( x)dx
L
cij 0 c( x) i ( x) j ( x)dx
其中，c(x)表示分布的粘滞阻尼特性。
一致节点荷载
L
vˆ 表示体系的形状，不随时间变化。
v 2vˆ sin(t ) 2v 2mvˆ sin(t ) kvˆ sin(t ) 0
k 2m vˆ 0
k 2m vˆ 0
即： k 2m 0
上式的N个根，表述体系可能存在的N个振型的频率。
1
2
3
N
2)
2
)
y32
(t
)
jˆ
32
s
in
(
2
t
2
)
1
jˆ
2i
yi
(t
)
jˆ3i
s
in(i
t
i
)
jˆ ni
1
jˆ 21
jˆ 31
jˆ 32
1
jˆ 22
将N个振型中的每一振型形式，用F表示N个振型所组成的方阵。
11 12 13 1N
结构动力响应数值计算方法对比分析
结构动力响应数值计算方法对比分析作者：李涵来源：《青年生活》2019年第21期摘要：中心差分法、纽马克法、威尔逊-法是结构动力学中常用的三种方法，为了系统的比较其优缺性，本文针对一个双自由度的体系，首先根据已知条件计算出振动微分方程，运用Matlab计算出可求出12个步长内相应的位移值，即精确解。
然后分别运用中心差分法，纽马克法，威尔逊-法求出其近似解;最后通过三种方法的近似解与精确解相对比，进而分析出三种计算方法的优缺性，为结构动力计算提供依据。
关键词：动力计算、中心差分法、纽马克法、威尔逊-法1、动力体系概况2、精确解推导针对该双自由度体系，理论推导出系统的位移表达式，通过代入各时刻周期得出位移在各时刻的具体数值，即位移精确解。
对位移方程求一阶导数得出速度方程，求二阶导数求出加速度方程。
代入各时刻的周期值，通过Matlab计算得出位移、速度、加速度的数值如下：3、三种数值计算方法3.1、中心差分法中心差分法是基于用有限差分代替位移对时间的求导，对位移一阶求导得到速度，对位移二阶求导得加速度。
通过Matlab计算出前4个步长所对应的位移响应。
3.2、纽马克法纽马克-β法是一种将线性加速度方法普遍化的方法。
通过Matlab计算出前4个步长所对应的位移响应。
3.3威尔逊-法通过Matlab计算出前4个步长所对应的位移响应。
4、近似解与精确解对比分析从上述结构的位移、速度、加速度可以看出，三种方法都能大致表示该体系大体运动趋势，并且误差较小。
其中，在描述物体位移时，中心差分法较后两种方法更为精确。
然而在描述速度和加速度时，中心差分法表现出了较大的误差，而纽马克和威尔逊法则能更详尽的表征物体速度和加速度。
5、结论中心差分法、纽马克法和威尔逊-θ法均是结构动力计算中的常用方法。
本文针对具体的计算实例，分别计算出三种方法的动力响应结果，并与精确解进行对比。
经过分析，中心差分法能更精确的表示物体位移响应，而纽马克和威尔逊法在表征物体速度和加速度方面相较于中心差分法更为精确，三种方法，各有其优缺点，应视具体情况采用相应的计算方法。
子结构拟动力试验的CD-Wilson θ法研究
Ｒｅ ｓ ｅ ａ ｒ ｃ ｈ ｏ ｎ ＣＤ － Ｗ ｉ ｌ ｓ ｏ ｎ ０ ｍｅ ｔ ｈ ｏ ｄ ｉ ｎ ｓ ｕ ｂ ｓ ｔ ｒ ｕ ｃ ｔ ｕｒ ｅ ｐ ｓ ｅ ｕ ｄｏ ｄ ｙ ｎ ａ ｍｉ ｃ ｔ ｅ ｓ ｔ ＷＡ Ｎ Ｇ Ｘ ｉ ａ ｏ ｆ ｅ ｎ ｇ ， Ｃ Ａ Ｉ Ｘ ｉ ｎ ｊ ｉ ａ ｎ ｇ ， Ｔ Ｉ Ａ Ｎ Ｓ ｈ ｉ ｚ ｈ ｕ ’
Ａｂ ｓ ｔ ｒ ａ ｃ ｔ ： Ｔ ｈ ｅ ａ ｄ ｖ ａ ｎ ｔ ａ ｇ ｅ ｓ ｏ ｆ ｃ ｏ ｍｂ ｉ ｎ ａ ｔ ｉ ｏ ｎ ｏ ｆ ｎ ｕ ｍｅ ｉ ｒ ｃ ａ ｌ ｉ ｎ ｔ ｅ ｇ ｒ ａ ｔ ｉ ｏ ｎ ｍｅ ｔ ｈ ｏ ｄ ｃ ａ ｌ ｌ ｂ ｅ ｃ ｏ ｎ ｓ ｉ ｄ ｅ ｒ ｅ ｄ ａ ｓ ｔ ｈ ｅ ｅ ｘ ｐ ｌ ｉ ｃ ｉ ｔ ｎ ｕ ｍｅ ｉｃ ｒ ｌ ａ ｉ ｎ ｔ ｅ ｒａ ｇ ｔ ｉ ｏ ｎ ｍｅ ｔ ｈ ｄ ｏ ｗ ｉ ｈｏ ｔ ｕ ｔ ｉ ｔ ｅ ｒ ａ ｔ ｉ ｏ ｎ ａ ｎ ｄ ｔ ｈ ｅ ｉ ｍｐ ｉ ｌ ｃ ｉ ｔ ｎ ｕ ｍｅ ｒ ｉ ｃ ｌ ａ ｉ ｎ ｔ ｅ ｒａ ｇ ｔ ｉ ｏ ｎ ｍｅ ｔ ｈ ｄ ｏ ｗ ｉ ｔ ｈ ｕ ｎ ｃ ｏ ｎ ｄ ｉ ｔ ｉ ｏ ｎ ａ ｌ ｓ ｔ ａ ｂ ｉ ｉ ｌ ｔ ｙ ｉ ｎ ｔ ｈ ｅ ｓ ｕ ｂ ｓ ｔ ｒ ｕ ｃ ｔ ｕ ｒ ｅ ｐ ｓ ｅ ｕ ｄ ｄｙ ｏ ｎ ａ ｍｉ ｃ ｔ ｅ ｓ ｔ ． Ａｃ ｃ ｏ ｒ ｄ ｉ ｎ ｇ ｔ ｏ ｈｅ ｔ Ｃ Ｄ・ Ｎｅ ｗ ｍａ ｒｋ ｍｅ ｔ ｈ ｄ ｏ ｉｔ ｗ ｈ ｌ ｒｇ ａ ｅ ｕ ｐ ｐ ｅ ｒ ｆ ｅｑ ｒ ｕ ｅ ｎ ｃ ｅ ｔ ｈ ａ ｔ ｌ ｅ ａ ｄ ｓ ｔ ｏ ｔ ｈ ｅ ｔ ｉ ｍｅ ｓ ｔ ｅ ｐ ｂ ｅ ｃ ｏ ｍｅ ｓ ｍａ ｌ ｌ ｅ ｒ ， ｔ ｈ ｅ ＣＤ－ Ｗｉ ｌ ｓ ｏ ｎ ０ ｍｅ ｔ ｈ ｄ ｏ ｉ ｓ ｃ ｏ ｍｂ ｉ ｎ ｅ ｄ ｂ ｙ ｅ ｘ ｐ ｌ ｉ ｃ ｉ ｔ ｃ ｅ ｎ ｔ ｅ ｒ ｄ ｉ ｆ ｅ ｒ ｅ ｎ ｃ ｅ ｍｅ ｔ ｈ ｄ ｏ ｎｄ ａ ｉ ｍｐ ｌ ｉ ｃ ｉ ｔ Ｗｉ ｌ ｓ ｏ ｎ － ０ ｍｅ ｔ ｈ ｄ． ｏ Ａｆ ｔ ｅ ｒ ｃ ｌｃ ａ ｕ ｌ ａ ｔ ｉ ｏ ｎ ａ ｎ ｄ ｍｏ ｄｉ ｉ ｆ ｃ ａ ｔ ｉ ｏ ｎ ｔ ｈ ｅ Ｗｉ ｌ ｓ ｏ ｎ － ０ ｍｅ ｔ ｈ ｄ， ｏ
基于显式Wilsonθ法的动载荷识别研究
收稿日期:２０１８－０４－１９基金项目:国家自然科学基金资助项目(５１７７５０９４).作者简介:范玉川(１９８８－)ꎬ男ꎬ河南新乡人ꎬ东北大学博士研究生ꎻ赵春雨(１９６３－)ꎬ男ꎬ辽宁黑山人ꎬ东北大学教授ꎬ博士生导师ꎻ张义民(１９５８－)ꎬ男ꎬ吉林长春人ꎬ东北大学"长江学者奖励计划"特聘教授ꎬ博士生导师.第４０卷第５期２０１９年５月东北大学学报(自然科学版)ＪｏｕｒｎａｌｏｆＮｏｒｔｈｅａｓｔｅｒｎＵｎｉｖｅｒｓｉｔｙ(ＮａｔｕｒａｌＳｃｉｅｎｃｅ)Ｖｏｌ.４０ꎬＮｏ.５Ｍａｙ２０１９㊀ｄｏｉ:１０.１２０６８/ｊ.ｉｓｓｎ.１００５－３０２６.２０１９.０５.０１３基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别研究范玉川１ꎬ赵春雨１ꎬ鲁㊀艳２ꎬ张义民１(１ 东北大学机械工程与自动化学院ꎬ辽宁沈阳㊀１１０８１９ꎻ２ 郑州信大先进技术研究院ꎬ河南郑州㊀４５０００１)摘㊀㊀㊀要:推导出多自由度动力学方程的Ｗｉｌｓｏｎ－θ数值算法显式表达形式ꎬ进而提出了一种显式Ｗｉｌｓｏｎ－θ的动载荷识别算法.该算法避免了Ｗｉｌｓｏｎ－θ算法的隐式迭代形式的迭代误差ꎬ在拥有显式算法特性的同时具备隐式算法的特性.当θ取合适的值时ꎬ该算法是无条件稳定的.通过悬臂梁的算例和实验对算法的识别效果进行了验证ꎬ并与传统的状态空间法的识别结果进行了对比.结果表明:该算法不仅能够对矩形载荷㊁谐波载荷和随机载荷进行准确地识别ꎬ并且比状态空间法的识别精度更高.关㊀键㊀词:Ｗｉｌｓｏｎ－θ法ꎻ显式表达ꎻ载荷识别ꎻ无条件稳定ꎻ状态空间法中图分类号:Ｏ３２６㊀㊀㊀文献标志码:Ａ㊀㊀㊀文章编号:１００５－３０２６(２０１９)０５－０６７３－０５ＲｅｓｅａｒｃｈｏｎＤｙｎａｍｉｃＬｏａｄＩｄｅｎｔｉｆｉｃａｔｉｏｎＢａｓｅｄｏｎＥｘｐｌｉｃｉｔＷｉｌｓｏｎ￣θＭｅｔｈｏｄＦＡＮＹｕ￣ｃｈｕａｎ１ꎬＺＨＡＯＣｈｕｎ￣ｙｕ１ꎬＬＵＹａｎ２ꎬＺＨＡＮＧＹｉ￣ｍｉｎ１(１ ＳｃｈｏｏｌｏｆＭｅｃｈａｎｉｃａｌ＆ＡｕｔｏｍａｔｉｏｎꎬＮｏｒｔｈｅａｓｔｅｒｎＵｎｉｖｅｒｓｉｔｙꎬＳｈｅｎｙａｎｇ１１０８１９ꎬＣｈｉｎａꎻ２ ＺｈｅｎｇｚｈｏｕＸｉｎｄａＩｎｓｔｉｔｕｔｅｏｆＡｄｖａｎｃｅｄＴｅｃｈｎｏｌｏｇｙꎬＺｈｅｎｇｚｈｏｕ４５０００１ꎬＣｈｉｎａ.Ｃｏｒｒｅｓｐｏｎｄｉｎｇａｕｔｈｏｒ:ＺＨＡＯＣｈｕｎ￣ｙｕꎬＥ￣ｍａｉｌ:ｃｈｙｚｈａｏ＠ｍａｉｌ.ｎｅｕ.ｅｄｕ.ｃｎ)Ａｂｓｔｒａｃｔ:ＴｈｅｅｘｐｌｉｃｉｔｅｘｐｒｅｓｓｉｏｎｏｆＷｉｌｓｏｎ￣θｎｕｍｅｒｉｃａｌａｌｇｏｒｉｔｈｍｆｏｒｍｕｌｔｉ￣ｄｏｆｓ(ｄｅｇｒｅｅｏｆｆｒｅｅｄｏｍｓ)ｄｙｎａｍｉｃｅｑｕａｔｉｏｎｉｓｄｅｒｉｖｅｄꎬａｓｗｅｌｌａｓａｎｅｘｐｌｉｃｉｔＷｉｌｓｏｎ￣θｄｙｎａｍｉｃｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎａｌｇｏｒｉｔｈｍｉｓｐｒｏｐｏｓｅｄꎬａｖｏｉｄｉｎｇｔｈｅｉｔｅｒａｔｉｏｎｅｒｒｏｒｗｈｉｌｅｋｅｅｐｉｎｇｔｈｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｔｈｅｉｍｐｌｉｃｉｔｉｔｅｒａｔｉｏｎａｌｇｏｒｉｔｈｍ.Ｔｈｅａｌｇｏｒｉｔｈｍｉｓｕｎｃｏｎｄｉｔｉｏｎａｌｌｙｓｔａｂｌｅｗｈｅｎａｐｐｌｙｉｎｇａｐｐｒｏｐｒｉａｔｅθｖａｌｕｅ.Ｔｈｅｒｅｃｏｇｎｉｔｉｏｎｅｆｆｅｃｔｏｆｔｈｅａｌｇｏｒｉｔｈｍｉｓｖｅｒｉｆｉｅｄｂｙａｎｅｘａｍｐｌｅａｎｄａｎｅｘｐｅｒｉｍｅｎｔｏｆａｃａｎｔｉｌｅｖｅｒｂｅａｍꎬａｎｄｔｈｅｒｅｓｕｌｔｓｗｅｒｅｃｏｍｐａｒｅｄｗｉｔｈｔｈｏｓｅｆｒｏｍｔｈｅｔｒａｄｉｔｉｏｎａｌｓｔａｔｅｓｐａｃｅｍｅｔｈｏｄ.Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅａｌｇｏｒｉｔｈｍｎｏｔｏｎｌｙｃａｎａｃｃｕｒａｔｅｌｙｉｄｅｎｔｉｆｙｔｈｅｒｅｃｔａｎｇｕｌａｒｌｏａｄꎬｔｈｅｈａｒｍｏｎｉｃｌｏａｄａｎｄｔｈｅｒａｎｄｏｍｌｏａｄꎬｂｕｔａｌｓｏｈａｓｈｉｇｈｅｒｒｅｃｏｇｎｉｔｉｏｎａｃｃｕｒａｃｙｔｈａｎｓｔａｔｅｓｐａｃｅｍｅｔｈｏｄ.Ｋｅｙｗｏｒｄｓ:Ｗｉｌｓｏｎ￣θꎻｅｘｐｌｉｃｉｔｆｏｒｍｕｌａｒꎻｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎꎻｕｎｃｏｎｄｉｔｉｏｎａｌｌｙｓｔａｂｌｅꎻｓｔａｔｅ￣ｓｐａｃｅｍｅｔｈｏｄ㊀㊀在机械系统设计过程中ꎬ动载荷是机械结构进行疲劳分析以及进行可靠性计算的基本依据.但是ꎬ在动态系统中ꎬ有些力是很难直接测量的ꎬ特别是结构系统内部各部件之间的相互作用力ꎬ在难以直接测量的情况下ꎬ就需要通过逆动力学的分析技术来得到这些力ꎬ因而开展动载荷识别技术的研究是很有必要的.动载荷识别方法主要分为频域法[１－２]和时域法[３]两大类.由于时域法能够识别各种类型的载荷ꎬ识别精度高ꎬ并且其识别结果具有明确的物理意义ꎬ因而ꎬ时域法越来越受到专家学者们的青睐.Ｌｉｕ等[４]推导了一种基于Ｎｅｗｍａｒｋ－β法的动载荷识别方法ꎬ将传统的隐式Ｎｅｗｍａｒｋ－β算法转化为Ａｘ＝ｂ方程解的显式形式ꎬ与隐式的Ｎｅｗｍａｒｋ－β法具有相同的特点ꎬ并与状态空间法进行了对比ꎬ表明该算法具有明显的优势.Ｌｉ等[５]提出了一种基于二阶泰勒级数展开的时域动态结构载荷识别方法ꎬ该算法将响应表示为一㊀㊀种泰勒级数的逼近形式ꎬ推导出一系列公式ꎬ并建立了系统响应㊁系统特性和输入激励相结合的显式离散方程.Ｌｉｕ等[６]提出了一种新的时域动态Ｇａｌｅｒｋｉｎ算法ꎬ将形状函数作为加权函数ꎬ建立了前向模型ＴＤＧＭꎬ与传统的格林函数法相比ꎬＴＤＧＭ能有效地克服噪声的影响ꎬ提高动态载荷识别的精度.在求逆运算中ꎬ载荷识别结果通常对结构模型中的响应和误差测量中的噪声非常敏感[７].张方等[８－９]对复杂结构的载荷识别进行了研究ꎬ推导了一种基于广义正交多项式特征技术的动载荷识别模型ꎬ可以在一定精度范围内通过有限的测量点信息对无限未知量的分布动载荷进行识别.Ａｌｌｅｎ等[１０]采用一种ＳＷＡＴ方法ꎬ可以同时识别出冲击型载荷和稳态动态载荷.本文将对Ｗｉｌｓｏｎ－θ法进行变换ꎬ推导其显式表达形式ꎬ提出一种基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别算法ꎬ并通过算例和实验验证了该算法的有效性.１㊀动载荷识别算法线性阻尼结构中ꎬ多自由度结构的动力学方程可以表示为Ｍ㊆ｘｔ＋Ｃ̇ｘｔ＋Ｋｘｔ＝Ｐｔ.(１)其中:ＭꎬＣ和Ｋ分别表示质量矩阵㊁阻尼矩阵以及刚度矩阵ꎻＰｔ是作用在结构上的外加载荷ꎻ㊆ｘｔꎬ̇ｘｔ和ｘｔ分别表示加速度响应㊁速度响应以及位移响应.本文假定阻尼为瑞利阻尼:Ｃ＝α１Ｍ＋α２Ｋ.(２)其中ꎬα１和α２表示阻尼系数.１ １㊀Ｗｉｌｓｏｎ－θ法的显式表达Ｗｉｌｓｏｎ－θ法假定在[ｔꎬｔ＋θΔｔ](θȡ１)的时间间隔内ꎬ加速度呈线性变化ꎬ如图１所示.令δ为自ｔ时刻开始的时间变量ꎬ适用于０ɤδɤθΔｔꎬ由线性加速度的假设可知ꎬ在适用范围内的加速度为㊆ｘｔ＋δ＝㊆ｘｔ＋δθΔｔ(㊆ｘｔ＋θΔｔ－㊆ｘｔ).(３)积分后可得̇ｘｔ＋δ＝̇ｘｔ＋㊆ｘｔδ＋δ２２θΔｔ(㊆ｘｔ＋θΔｔ－㊆ｘｔ)ꎬ(４)ｘｔ＋δ＝ｘｔ＋̇ｘｔδ＋１２㊆ｘｔδ２＋δ３６θΔｔ(㊆ｘｔ＋θΔｔ－㊆ｘｔ).(５)若δ＝θΔｔꎬ由式(４)和式(５)可得ｔ＋θΔｔ瞬时的速度和位移:̇ｘｔ＋θΔｔ＝̇ｘｔ＋θΔｔ２(㊆ｘｔ＋θΔｔ＋㊆ｘｔ)ꎬ(６)ｘｔ＋θΔｔ＝ｘｔ＋θΔｔ̇ｘｔ＋θ２Δｔ２６(㊆ｘｔ＋θΔｔ＋２㊆ｘｔ).(７)图１㊀Ｗｉｌｓｏｎ－θ法模型Ｆｉｇ １㊀ＭｏｄｅｌｏｆＷｉｌｓｏｎ￣θ联立式(６)和式(７)ꎬ可以用ｔ＋θΔｔ时刻的位移表示ｔ＋θΔｔ时刻的加速度和速度ꎬ即㊆ｘｔ＋θΔｔ＝６θ２Δｔ２(ｘｔ＋θΔｔ－ｘｔ)－６θΔｔ̇ｘｔ－２㊆ｘｔꎬ(８)̇ｘｔ＋θΔｔ＝３θΔｔ(ｘｔ＋θΔｔ－ｘｔ)－２̇ｘｔ－θΔｔ２㊆ｘｔ.(９)则ｔ＋θΔｔ时刻的动力方程可以表示为Ｍ㊆ｘｔ＋θΔｔ＋Ｃ̇ｘｔ＋θΔｔ＋Ｋｘｔ＋θΔｔ＝Ｐｔ＋θΔｔ.(１０)式中:Ｐｔ＋θΔｔ＝Ｐｔ＋θ(Ｐｔ＋Δｔ－Ｐｔ).(１１)将式(８)和式(９)以及式(１１)代入式(１０)ꎬ即得关于ｘｔ＋θΔｔ的求解方程为ｘｔ＋θΔｔ＝＾Ｋ－１＾Ｐｔ＋θΔｔ.(１２)式中:＾Ｋ＝Ｋ＋３θΔｔＣ＋６θ２Δｔ２Ｍꎬ(１３)＾Ｐｔ＋θΔｔ＝Ｐｔ＋θ(Ｐｔ＋Δｔ－Ｐｔ)＋Ｍ(６θ２Δｔ２ｘｔ＋６θΔｔ̇ｘｔ＋２㊆ｘｔ)＋Ｃ(３θΔｔｘｔ＋２̇ｘｔ＋θΔｔ２㊆ｘｔ).(１４)求解方程式(１２)ꎬ得ｘｔ＋θΔｔ.再把ｘｔ＋θΔｔ代入式(８)就可获得㊆ｘｔ＋θΔｔ.在式(３)中取δ＝Δｔꎬ同时将式(８)代入ꎬ可得㊆ｘｔ＋Δｔ＝６θ３Δｔ２(ｘｔ＋θΔｔ－ｘｔ)－６θ２Δｔ̇ｘｔ＋(１－３θ)㊆ｘｔ.(１５)取δ＝Δｔꎬ将式(３)分别代入式(４)和式(５)ꎬ有̇ｘｔ＋Δｔ＝̇ｘｔ＋Δｔ２(㊆ｘｔ＋Δｔ＋㊆ｘｔ)ꎬ(１６)ｘｔ＋Δｔ＝ｘｔ＋Δｔ̇ｘｔ＋Δｔ２６(㊆ｘｔ＋Δｔ＋２㊆ｘｔ).(１７)由式(１２)~式(１４)可得４７６东北大学学报(自然科学版)㊀㊀㊀第４０卷㊀㊀ｘｔ＋θΔｔ＝＾Ｋ－１(１－θ)Ｐｔ＋θ＾Ｋ－１Ｐｔ＋Δｔ＋(６θ２Δｔ２Ｍ＾Ｋ－１＋３θΔｔＣ＾Ｋ－１)ｘｔ＋(６θΔｔＭ＾Ｋ－１＋２Ｃ＾Ｋ－１)̇ｘｔ＋(２Ｍ＾Ｋ－１＋θΔｔ２Ｃ＾Ｋ－１)㊆ｘｔ.(１８)将式(１８)代入式(１５)ꎬ已知恒等式:Ｉ＝＾Ｋ－１＾Ｋꎬ可得㊆ｘｔ＋Δｔ＝Ｃ０Ｐｔ＋Ｃ１Ｐｔ＋Δｔ＋Ｃｄｘｔ＋Ｃｖ̇ｘｔ＋Ｃａ㊆ｘｔ.(１９)其中:Ｃ０＝６θ３Δｔ２＾Ｋ－１(１－θ)ꎻＣ１＝６θ２Δｔ２＾Ｋ－１ꎻＣｄ＝－６θ３Δｔ２＾Ｋ－１ＫꎻＣｖ＝６θ２Δｔ(６θ２Δｔ２Ｍ＾Ｋ－１＋２θΔｔＣ＾Ｋ－１－１)ꎻＣａ＝１２θ３Δｔ２Ｍ＾Ｋ－１＋３θ２ΔｔＣ＾Ｋ－１＋１－３θ.将式(１９)代入式(１６)可得̇ｘｔ＋Δｔ＝Ｂ０Ｐｔ＋Ｂ１Ｐｔ＋Δｔ＋Ｂｄｘｔ＋Ｂｖ̇ｘｔ＋Ｂａ㊆ｘｔ.(２０)其中:Ｂ０＝Δｔ２Ｃ０ꎻＢ１＝Δｔ２Ｃ１ꎻＢｄ＝Δｔ２ＣｄꎻＢｖ＝Δｔ２Ｃｖ＋１ꎻＢａ＝Δｔ２(Ｃａ＋１).将式(１９)代入式(１７)可得ｘｔ＋Δｔ＝Ａ０Ｐｔ＋Ａ１Ｐｔ＋Δｔ＋Ａｄｘｔ＋Ａｖ̇ｘｔ＋Ａａ㊆ｘｔ.(２１)其中:Ａ０＝Δｔ２６Ｃ０ꎻＡ１＝Δｔ２６Ｃ１ꎻＡｄ＝Δｔ２６Ｃｄ＋１ꎻＡｖ＝(Δｔ６Ｃｖ＋１)ΔｔꎻＡａ＝Δｔ２６(Ｃａ＋２).在此用ｘｉꎬ̇ｘｉꎬ㊆ｘｉ和Ｐｉ表示第ｉ时刻式(１９)~式(２１)中的位移㊁速度㊁加速度和激励ꎬ并将其表示为矩阵形式:ｘｉ＋１̇ｘｉ＋１㊆ｘｉ＋１éëêêêùûúúú＝Ａ０㊀Ａ１Ｂ０㊀Ｂ１Ｃ０㊀Ｃ１éëêêêùûúúúＰｉＰｉ＋１éëêêùûúú＋Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀Ｃａéëêêêùûúúúｘｉ̇ｘｉ㊆ｘｉéëêêêùûúúú.(２２)在第ｉ时刻的响应可以表示为ｘｉ̇ｘｉ㊆ｘｉéëêêêùûúúú＝ðｉ－１ｊ＝０ＡｄＡｖＡａＢｄＢｖＢａＣｄＣｖＣａéëêêêùûúúúｊＡ０Ａ１Ｂ０Ｂ１Ｃ０Ｃ１éëêêêùûúúúＰｉ－ｊ－１Ｐｉ－ｊ[]＋ＡｄＡｖＡａＢｄＢｖＢａＣｄＣｖＣａéëêêêùûúúúｉｘ０̇ｘ０㊆ｘ０éëêêêùûúúú.(２３)式中ꎬ两个指数ｉ和ｊ分别表示相应矩阵的幂.式(２３)是一种新型的Ｗｉｌｓｏｎ－θ法的显式表达ꎬ每一个时间步中位移㊁速度㊁加速度响应可被同时求解出来ꎬ而通常的方法是利用每个时间步的迭代算法来计算ꎬ显然这种显式的算法更有优势.１ ２㊀基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别令:ｙｉ＝ｘｉ̇ｘｉ㊆ｘｉéëêêêùûúúú－Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀Ｃａéëêêêùûúúúｉｘ０̇ｘ０㊆ｘ０éëêêêùûúúú.则式(２３)可以改写为ｙｉ＝ðｉ－１ｊ＝０Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀ＣａéëêêêùûúúúｊＡ０㊀Ａ１Ｂ０㊀Ｂ１Ｃ０㊀Ｃ１éëêêêùûúúúＰｉ－ｊ－１Ｐｉ－ｊ[].(２４)令:Ｈｋ＝Ａｄ㊀Ａｖ㊀ＡａＢｄ㊀Ｂｖ㊀ＢａＣｄ㊀Ｃｖ㊀Ｃａéëêêêùûúúú(ｋ－１)Ａ０㊀Ａ１Ｂ０㊀Ｂ１Ｃ０㊀Ｃ１éëêêêùûúúúꎬ(２５)Ｒｉ＝Ｐｉ－１Ｐｉéëêêùûúú.(２６)式(２９)可以写成从时间段１到ｎｔ上的矩阵的卷积形式:Ｙ＝ＨＦ.(２７)其中:Ｙ＝[ｙＴ１ꎬｙＴ２ꎬ ꎬｙＴｎｔ]ＴꎻＨ＝Ｈ１０ ０Ｈ２Ｈ１ ０⋮⋮⋮ＨｎｔＨｎｔ－１ Ｈ１éëêêêêêùûúúúúúꎻＦ＝Ｒ１Ｒ２⋮Ｒｎｔéëêêêêêùûúúúúú.式(２７)可以改写成载荷识别的表达式Ｆ＝Ｈ－１Ｙ.(２８)对于一个给定的系统ꎬＨ是常数ꎬＹ可以从系统测量的响应中得到ꎬ考虑到式(２８)可能存在不适定性ꎬ可以通过Ｔｉｋｈｏｎｏｖ正则化方法来确定Ｆꎬ在阻尼最小二乘意义下对目标函数进行优化.Ｊ(Ｆꎬλ)＝ Ｙ－ＨＲ ２＋λ Ｆ ２.(２９)其中ꎬλ是正则化参数ꎬ它的值可以通过Ｌ曲线法来确定.基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别方法可以归纳为以下步骤:１)建立系统的有限元模型ꎬ确定系统的质量矩阵Ｍ㊁刚度矩阵Ｋ和阻尼矩阵Ｃꎻ２)选取时间步Δｔ和Ｗｉｌｓｏｎ－θ法的参数θꎬ参数θ的值可以取为１ ４ꎬ从而保证算法的无条件稳定(在Ｗｉｌｓｏｎ－θ法中ꎬ只要θ的值大于５７６第５期㊀㊀㊀范玉川等:基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别研究㊀㊀１ ３７ꎬ那么该算法就无条件稳定ꎬ但是θ取的过大ꎬ截断误差增大ꎬ精度会下降).３)为递推式(２２)计算矩阵Ａ１ꎬＡ２ꎬＡｄꎬＡｖꎬＡａꎬＢ１ꎬＢ２ꎬＢｄꎬＢｖꎬＢａꎬＣ１ꎬＣ２ꎬＣｄꎬＣｖ和Ｃａꎻ４)计算式(２５)中的矩阵Ｈｋ(ｋ＝０ꎬ ꎬｎｔ－１)和式(２８)中的组合矩阵Ｈꎻ５)通过实验或仿真计算获得系统的响应数据㊆ｘꎬ̇ｘ和ｘꎻ６)利用式(２９)中的Ｔｉｋｈｏｎｏｖ正则化方法来确定Ｆ.２㊀仿真算例及抗噪性能研究如图２所示ꎬ一矩形截面悬臂梁长０ ６４ｍꎬ截面尺寸为:宽０ ０５６ｍ㊁高０ ００８ｍꎬ弹性模量Ｅ＝２００ＧＰａꎬ材料密度为７８４０ｋｇ/ｍ３ꎬ该悬臂梁被划分为１８个单元ꎬ共１９个节点.图２㊀悬臂梁结构图Ｆｉｇ ２㊀Ｓｔｒｕｃｔｕｒｅｏｆｃａｎｔｉｌｅｖｅｒｂｅａｍ在第１０个节点上施加载荷ꎬ为了尽可能全面地验证本文载荷识别算法的准确性ꎬ分别施加方波载荷㊁谐波载荷和随机载荷ꎬ考虑到噪声的影响ꎬ在响应数据中加入５％的随机噪声ꎬ选取第１２个节点的响应信息进行载荷识别ꎬ时间间隔Δｔ＝０ ００１ｓ.选取第１２个节点上的数据进行载荷识别计算ꎬ本文算法与状态空间法进行对比ꎬ如图３所示ꎬｐ２为方波载荷㊁ｐ３为谐波载荷㊁ｐ４为随机载荷.两种方法对于这三种不同类型的载荷都能够比较准确的识别ꎬ但是状态空间法的识别效果明显要差一些ꎬ而基于显式Ｗｉｌｓｏｎ－θ法的载荷识别方法识别得很好ꎬ对实际载荷的还原度较高.为便于量化分析基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别算法的识别精度ꎬ可引入式(３０)计算误差:Ｅｒｒｏｒ＝Ｆｉｄ－ＦｒｅａｌＦｒｅａｌˑ１００％.(３０)式中:Ｅｒｒｏｒ表示相对识别误差ꎻＦｉｄ表示识别得到的外激励ꎻＦｒｅａｌ表示真实的外激励.运用不同的方法ꎬ在不同的采样频率和采样时间下研究动载荷识别的误差ꎬ结果列于表１中.图３㊀载荷识别结果Ｆｉｇ ３㊀Ｒｅｓｕｌｔｓｏｆｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎ(ａ) ｐ２ꎻ(ｂ) ｐ３ꎻ(ｃ) ｐ４.表１㊀两种载荷识别方法的识别误差对比Ｔａｂｌｅ１㊀Ｃｏｍｐａｒｉｓｏｎｏｆｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｅｒｒｏｒｓｗｉｔｈｄｉｆｆｅｒｅｎｔｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｍｅｔｈｏｄｓ％载荷Ｗｉｌｓｏｎ－θ法状态空间法ｐ２３ ２１５ ０４Ｐ３０ ６１３ ８４Ｐ４３ ３０５ １６㊀㊀数据表明ꎬ本文提出的基于显式Ｗｉｌｓｏｎ－θ法的载荷识别方法的识别精度明显高于基于状态空间法的载荷识别方法的识别精度ꎬ并且两种算法都显示谐波载荷的识别精度要远远好于方波载荷和随机载荷的识别精度.６７６东北大学学报(自然科学版)㊀㊀㊀第４０卷㊀㊀㊀㊀图４显示了随着噪声的增大载荷识别精度的变化情况ꎬ两种方法的识别误差都是随着噪声的增大而增大ꎬ但本文算法随噪声的增大识别误差变化较小ꎬ表明本文算法抗噪性能更好㊁识别误差更小.图４㊀载荷识别误差随噪声变化情况Ｆｉｇ ４㊀Ｖａｒｉａｔｉｏｎｏｆｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｅｒｒｏｒｗｉｔｈｎｏｉｓｅ３㊀实验验证实验用的悬臂梁模型参数与仿真悬臂梁模型参数一致ꎬ选用ＹＥ６２５１振动力学实验系统.激励载荷为谐波载荷ꎬ激励频率设置为１０Ｈｚꎬ算法中采用有限元法进行建模ꎬ代入悬臂梁的相关参数ꎬ选用测量得到的第１１个节点上的响应数据进行载荷识别ꎬ识别结果如图５所示.由图可知ꎬ该算法对实验数据的识别较准确ꎬ识别的载荷曲线与实际的载荷曲线基本吻合.图５㊀实验识别结果Ｆｉｇ ５㊀Ｒｅｓｕｌｔｓｏｆｅｘｐｅｒｉｍｅｎｔ由于实际工程中响应数据的测量位置往往是有限的ꎬ不能够随意选取ꎬ因而在本实验中分别选取第４㊁第７㊁第１０㊁第１３和第１６个节点的测量数据进行载荷的识别ꎬ同样利用式(３０)进行识别误差的计算ꎬ计算结果分别为５ ６７％ꎬ５ ７９％ꎬ５ ３８％ꎬ５ ８６％ꎬ５ ９８％.从结果可以看到ꎬ选取不同的测量点ꎬ识别得到的载荷与实际载荷的误差没有太大差别ꎬ也就是说测量点的选取对识别结果影响不大ꎻ该算法能够适应实际工程的应用环境.４㊀结㊀㊀论１)利用隐式Ｗｉｌｓｏｎ－θ法ꎬ推导出了它的一种显式表达形式ꎬ进而提出了一种基于显式Ｗｉｌｓｏｎ－θ法的载荷识别算法ꎬ该算法既具有显式算法的优势ꎬ同时具有隐式算法的特点.２)通过算例验证了本文算法在识别矩形波载荷㊁谐波载荷和随机载荷方面的识别精度都高于状态空间法的识别精度ꎬ并且在算例中加入了５％的随机噪声.实验结果表明:算法不仅能够对矩形载荷㊁谐波载荷和随机载荷准确识别ꎬ并且识别精度比状态空间法的识别精度更高.参考文献:[１]㊀ＪｉｅＨꎬＺｈａｎｇＸ.Ａｎｏｐｔｉｍｉｚａｔｉｏｎｍｅｔｈｏｄｏｆｌｏａｄｉｄｅｎ￣ｔｉｆｉｃａｔｉｏｎｉｎｆｒｅｑｕｅｎｃｙｄｏｍａｉｎ[Ｊ].Ｎｏｉｓｅ＆ＶｉｂｒａｔｉｏｎＣｏｎｔｒｏｌꎬ２００９ꎬ２９(６):３４－３６[２]㊀ＨｅＺＣꎬＬｉｎＸＹꎬＬｉＥ.Ａｎｏｖｅｌｍｅｔｈｏｄｆｏｒｌｏａｄｂｏｕｎｄｓｉｄｅｎｔｉｆｉｃａｔｉｏｎｆｏｒｕｎｃｅｒｔａｉｎｓｔｒｕｃｔｕｒｅｓｉｎｆｒｅｑｕｅｎｃｙｄｏｍａｉｎ[Ｊ].ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＭｅｔｈｏｄｓꎬ２０１７(３):１８５００５１.[３]㊀ＬａｗＳＳꎬＣｈａｎＴＨＴꎬＺｈｕＱＸ.Ｒｅｇｕｌａｒｉｚａｔｉｏｎｉｎｍｏｖｉｎｇｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＥｎｇｉｎｅｅｒｉｎｇＭｅｃｈａｎｉｃｓＡＳＣＥꎬ２００１ꎬ１２７(２):１３６－１４８.[４]㊀ＬｉｕＫꎬＬａｗＳＳꎬＺｈｕＸＱꎬｅｔａｌ.Ｅｘｐｌｉｃｉｔｆｏｒｍｏｆａｎｉｍｐｌｙ￣ＣＩＴｍｅｔｈｏｄｆｏｒｉｎｖｅｒｓｅｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＳｏｕｎｄａｎｄＶｉｂｒａｔｉｏｎꎬ２０１４ꎬ３３:７３０－７４４.[５]㊀ＬｉＸＷꎬＤｅｎｇＺＭ.Ｉｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｄｙｎａｍｉｃｌｏａｄｓｂａｓｅｄｏｎｓｅｃｏｎｄ￣ｏｒｄｅｒｔａｙｌｏｒ￣ｓｅｒｉｅｓｅｘｐａｎｓｉｏｎｍｅｔｈｏｄ[Ｊ].ＳｈｏｃｋａｎｄＶｉｂｒａｔｉｏｎꎬ２０１６(２０１６):１－９.[６]㊀ＬｉｕＪꎬＭｅｎｇＸꎬＪｉａｎｇＣꎬｅｔａｌ.Ｔｉｍｅ￣ｄｏｍａｉｎＧａｌｅｒｋｉｎｍｅｔｈｏｄｆｏｒｄｙｎａｍｉｃｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｆｏｒＭｕｍｅｒｉｃａｌＭｅｔｈｏｄｓｉｎＥｎｇｉｎｅｅｒｉｎｇꎬ２０１６ꎬ１０５:６２０－６４０.[７]㊀ＱｉａｎＢꎬＺｈａｎｇＸꎬＷａｎｇＣꎬｅｔａｌ.Ｓｐａｒｓｅｒｅｇｕｌａｒｉｚａｔｉｏｎｆｏｒｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎｕｓｉｎｇｄｉｃｔｉｏｎａｒｉｅｓ[Ｊ].ＪｏｕｒｎａｌｏｆＳｏｕｎｄａｎｄＶｉｂｒａｔｉｏｎꎬ２０１６ꎬ３６８:７１－８６.[８]㊀张方ꎬ秦远田ꎬ邓吉宏.复杂分布动载荷识别技术研究[Ｊ].振动工程学报ꎬ２００６ꎬ１９(１):８１－８５.(ＺｈａｎｇＦａｎｇꎬＱｉｎＹｕａｎ￣ｔｉａｎꎬＤｅｎｇＪｉ￣ｈｏｎｇ.Ｒｅｓｅａｒｃｈｏｆｉｄｅｎｔｉｆｉｃａｔｉｏｎｔｅｃｈｎｏｌｏｇｙｏｆｄｙｎａｍｉｃｌｏａｄｄｉｓｔｒｉｂｕｔｅｄｏｎｔｈｅｓｔｒｕｃｔｕｒｅ[Ｊ].ＪｏｕｒｎａｌｏｆＶｉｂｒａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇꎬ２００６ꎬ１９(１):８１－８５.)[９]㊀徐菁ꎬ张方ꎬ姜金辉ꎬ等.运用数值迭代的动载荷识别算法[Ｊ]ꎬ振动工程学报ꎬ２０１４ꎬ２７(１５):７０２－７０７.(ＸｕＪｉｎｇꎬＺｈａｎｇＦａｎｇꎬＪｉａｎｇＪｉｎ￣ｈｕｉꎬｅｔａｌ.Ａｎａｌｇｏｒｉｔｈｍｏｆｄｙｎａｍｉｃｌｏａｄｉｄｅｎｔｉｆｉｃａｔｉｏｎｂａｓｅｄｏｎｎｕｍｅｒｉｃａｌｉｔｅｒａｔｉｏｎꎬ[Ｊ].ＪｏｕｒｎａｌｏｆＶｉｂｒａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇꎬ２０１４ꎬ２７(１５):７０２－７０７.)[１０]ＡｌｌｅｎＭＳꎬＣａｒｎｅＴＧ.Ｄｅｌａｙｅｄｍｕｌｔｉ￣ｓｔｅｐｉｎｖｅｒｓｅｓｔｒｕｃｔｕｒａｌｆｉｌｔｅｒｆｏｒｒｏｂｕｓｔｆｏｒｃｅｉｄｅｎｔｉｆｉｃａｔｉｏｎ[Ｊ].ＭｅｃｈａｎｉｃａｌＳｙｓｔｅｍｓａｎｄＳｉｇｎａｌＰｒｏｃｅｓｓｉｎｇꎬ２００８ꎬ２２(５):１０３６－１０５４.７７６第５期㊀㊀㊀范玉川等:基于显式Ｗｉｌｓｏｎ－θ法的动载荷识别研究。
Wilson_方法在结构动力分析中的应用研究
第20卷 第3期2003年08月工 程 数 学 学 报JOURNA L OF E NGI NEERI NG MATHE MATICSV ol.20N o.3Aug.2003文章编号:100523085(2003)0320082205Wilson 2θ方法在结构动力分析中的应用研究Ξ马小燕1, 王一兵2(12西安电力高等专科学校数学组 22西安交通大学建筑工程与力学学院,西安710049)摘 要:介绍了逐步积分Wils on 2θ方法在砖混建筑结构动力分析中的应用。
通过对单质点结构剪切模型在正弦激励下的模拟计算,研究Wils on 2θ方法中各参数对该方法计算精度的影响。
提出Wils on 2θ方法适合砖混建筑结构动力分析的最优参数选择。
关键词:Wils on 2θ方法;逐步积分;剪切模型;精度;稳定性分类号:AMS (2000)70J35 中图分类号:T B115 文献标识码:A1 前言逐步积分方法在工程实际当中应用非常广泛。
许多工程振动问题无法得到其精确解。
例如在地震工程问题中,研究建筑结构在地震激励作用下结构的响应时,由于地震激励的复杂性(地震波的随机性,基础场地的个性等)无法得到其精确解。
这种情况就需要用数值方法求解,逐步积分是求解复杂振动问题最常用的数值方法。
研究逐步积分方法的求解精度,对工程实际应用是非常重要。
本文采用W ils on 2θ方法,对建筑结构动力分析模型(剪切模型),受正弦地震激励作用下,进行了模拟计算。
并与精确解进行比较,分析研究W ils on 2θ方法的计算精度和稳定性。
图1 剪切模型2 Wilson 2θ法2.1 数学模型和假设结构动力分析的单质点剪切模型如图1所示。
m 为质点质量,k 层间侧向刚度,当地面有水平加速度¨x g 时质点运动方程为(无阻尼)。
m ¨x +kx =-¨x g(1)式中x 表示质点相对地面的水平位移,¨x 表示质点相对地面的水平加速度,¨x g 为地面运动加速度(地震加速度)。
基于Newmark-β法的Matlab简单程序编写
基于Newmark-β法的Matlab简单程序编写李坤明;黄菲【摘要】我们知道,对于结构动力学问题的直接求解往往是比较困难的,而Newmark-β法的出现,很巧妙地将这种困难化解.Newmark-β法是一种时程分析的方法,它将动力学微分方程的求解问题转化为若干个代数方程逐步求解,是一种隐式积分法,可以得到足够精确的解.本文基于Newmark-β法的理论知识,用Matlab 对一个简单动力学问题的求解进行编程运算.【期刊名称】《四川建材》【年(卷),期】2018(044)009【总页数】3页(P65-66,90)【关键词】Newmark-β法;编程;结构力学【作者】李坤明;黄菲【作者单位】西华大学土木建筑与环境学院,四川成都 610039;西华大学土木建筑与环境学院,四川成都 610039【正文语种】中文【中图分类】TP3131 逐步积分法在结构动力学问题的分析中，对于承受任意动荷载的线性结构，我们可以采用Duhamel、频域分析等方法，这些方法都可以很方便地求解出所需要的结果。
但是上面的两种方法都用到了叠加原理，所以它们只适用于求解线性结构体系，同时也必须要求这种体系的特性在反应过程中不能发生变化[1]。
然而，在另一方面，我们的实际生活中有很多重要的结构动力学问题，整个反应体系并不能单纯地认为是线性变化的，比如，在足以引起强烈破坏的地震作用下，很多建筑物的反应，我们都必须考虑非线性反应的影响。
所以，我们还必须要发展对于非线性结构的动力响应的分析方法。
对于非线性体系的动力响应问题，发展出了最有效的分析方法—逐步积分法[1]。
在这种方法中，我们采用一系列的短时间增量Δt来计算反应，而通常为了方便起见，我们将Δt取为等间距的时间步长。
在每个时间间隔的起点和终点建立动力平衡方程，并以一个假设的反应机理为根据，近似地计算在时间间隔内体系的运动(通常忽略在时间间隔内产生的不平衡)，体系的非线性特性可以用每个时间间隔的起点所求得的当前变形状态的新特性来说明。
多自由度体系wilson-θ法程序编写
多自由度体系wilson-θ法程序编写多自由度体系Wilson-θ法是一种广泛应用于多体动力学和结构动力学领域的数值计算方法。
本文将介绍如何使用Python编程语言编写多自由度体系Wilson-θ法的程序。
一、引言多自由度体系Wilson-θ法是一种基于有限元法的数值计算方法，适用于求解多体动力学和结构动力学中的问题。
该方法通过将体系分解为一系列有限元子系统，并采用θ矩阵方法进行求解，能够有效地处理大规模的多自由度体系。
二、程序编写1. 导入必要的库和模块在编写程序之前，需要导入必要的库和模块，包括numpy、scipy 和matplotlib等。
这些库提供了必要的数学运算、数值分析和图形绘制等功能。
2. 定义体系结构和有限元节点首先需要定义多自由度体系的结构和有限元节点的位置。
可以使用网格划分工具将体系划分为有限元网格，并定义每个节点的位置和编号。
3. 构建有限元矩阵和求解器使用Wilson-θ法进行数值计算，需要构建有限元矩阵和求解器。
该矩阵可以采用三角矩阵的形式进行表示，并使用θ矩阵方法进行求解。
在程序中，需要实现矩阵的构建、求解器的初始化等操作。
4. 迭代求解体系响应使用构建好的矩阵和求解器，可以进行迭代求解多自由度体系的响应。
在每次迭代中，需要输入当前时刻的体系响应作为初值，并输出下一时刻的响应结果。
5. 结果可视化最后，可以使用matplotlib等库将求解得到的响应结果进行可视化。
可以将时间历程、振型、频率响应等结果进行绘制，以便更好地分析体系的动态特性。
三、示例代码以下是一个简单的示例代码，用于演示如何使用Python编程语言编写多自由度体系Wilson-θ法的程序。
代码中假设体系由3个自由度的弹簧-质量系统组成，采用三角形矩阵进行求解。
```pythonimport numpy as npfrom scipy.sparse import csc_matrix, dia_matriximport matplotlib.pyplot as plt# 定义体系结构和有限元节点nodes = np.array([[0], [0.5], [1]]) # 节点位置数组degrees_of_freedom = 3 # 自由度数量system_size = len(nodes) # 体系大小node_indices = np.arange(system_size) # 节点编号数组# 构建有限元矩阵和求解器theta_matrix = csc_matrix(dia_matrix(system_size - degrees_of_freedom, 0)) # θ矩阵mass_matrix = csc_matrix(np.diag([0.5, 0.5, 1])) # 质量矩阵solution = np.zeros((system_size, degrees_of_freedom)) # 初始响应数组forces = np.zeros((system_size, degrees_of_freedom)) # 输出力数组forces[:degrees_of_freedom] = np.zeros((system_size, degrees_of_freedom)) # 初始输出力数组为零向量solver = theta_matrix.dot(solution) +theta_matrix.dot(forces) + mass_matrix # 初始化求解器theta_vector = np.zeros(system_size) # θ向量用于控制有限元矩阵的构造和更新# 进行迭代求解体系响应for iteration in range(100): # 迭代次数限制为100次response = solver.dot(theta_vector) # 输入当前时刻的响应作为初值进行迭代求解下一时刻的响应结果输出为力向量output_forces在每个节点上作用在体系的上结果可与theta向量用于控制有限元矩阵的构造和更新为了演示程序的基本结构和流程以上给出了一个简单的示例代码其中包含的主要内容有定义体系结构和有限元节点构建有限元矩阵和求解器以及进行迭代求解体系响应结果可视化等当然在实际应用中可能还需要考虑更多的因素例如如何处理边界条件如何处理体系的非线性特性等等因此在实际应用中需要根据具体问题对程序进行适当的修改和优化以下是一些可能需要的注意事项和技巧：1. 选择合适的有限元网格划分工具和算法，以确保计算的精度和效率。
结构动力学数值积分方法-wilson
(10)
求解(10)式，得到位移 u t t 。然后采用(6)式计算得到加速度 a t t 。
a t t
6
t
2
6 u t t u t t v t 2a t
(11)
将加速度 a t t 代入到(1)式中，并令 t ，可以计算出加速度 a t t 。
其中，激振力可以用下式表达。
p t t p t p t t p t
(8)
(9)
由振动控制方程得到关于 u t t 表达式如下：
6m 3c k u t t 2 t t p t p t t p t 6 6 t 3 m u t v t 2 a t u t 2v t a t c 2 t 2 t t
2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5
0
1
2
3
4
5
6
7
8
9
3/3
西安交大 航天学院 仲继泽(hnskzjz@)
结构动力学数值积分方法(Wilson-θ) 及其 Matlab 实现
Xjtu-Zjz
Wilson-θ 方法是线性加速度方法的一种拓展。 假设加速度在时间段 t , t t 内线性变化，首先计算时刻 t t 的位移，然后通过内插得到 t t 时刻的位移 (u)、速度(v)及加速度(a)。其中 1 ，相关研究证明， 1.37 时，Wilson-θ 方法 是无条件稳定的。一般情况下，取 1.4 (至于为什么，本人没有做深入的研究)。 Wilson-θ 方法理论推导如下： 记 t 为时间段 t , t t 内的任一时刻，对于一个单自由度振动系统来说， 其初始时刻 t 的位移 u t 、速度 v t 及加速度 a t 为已知量，其 t 时刻的加速 度可以用 t 时刻和 t t 时刻的加速度线性表达，见下式。
f6_??????
Consistent-mass
M
C
e
ml 420
54 22l
156 13l
13l 4l 2
22l 3l 2
13l 22l 3l 2
4l 2
静力凝聚法
当采用集中质量阵，而结构体系的自由度又存在转角时，转动自由 度的惯性力为0，此时可以采用静力凝聚法，消去无质量的自由 度。若平动和转动自由度分别用下标t和θ区分，则采用集中质量 方法时，存在转角自由度体系的运动方程可写为如下形式，
T T (q1, q2 ,, qN , q1, q2 ,, qN )
V V (q1, q2 ,, qN )
Wnc Q1q1 Q2q2 QNqN q1, q2, …, qN为广义坐标，
t2
(T V )dt
t2 Wncdt 0
t1
t1
பைடு நூலகம்
Q1 , Q2, …, QN为非保守力，例如外力、阻尼力等。
下面通过算例来介绍如何应用Lagrange方程，从算例中可 以看到，用Lagrange运动方程建立的运动方程不限于线 性问题。
6.2 Lagrange运动方程
算例6.1 如图所示一复合摆，摆的杆长分别为l1和l2，摆 的 质 量 分 别 为 m1 和 m2 ， 忽 略 杆 的 分 布 质 量 ， 采 用 Lagrange方程建立体系无阻尼自由运动方程。
结构动力学
（2003秋）
结构动力学
第六章
多自由度体系的运动方程
第六章 多自由度体系的运动方程
以前各章讨论的对象均为单自由度体系，它的运 动仅需一个运动方程来描述，求解这个运动方 程，就可以得到单自由度体系的位移、速度和 加速度以及能量等。
工程中所涉及的结构一般都是多自由度的，例如 多层建筑结构、大跨桥梁结构、空间网架结构 等等。为合理反映振动过程中惯性力的影响， 需要采用更多的自由度描述结构体系的质量分 布并确定体系的变形。
结构动力学方程的数值解法研究
结构动力学方程的数值解法研究
韩爱红;钱晓军
【期刊名称】《低温建筑技术》
【年(卷),期】2015(037)008
【摘要】本文介绍了求解结构动力学方程的数值积分方法,主要包括:线性加速度法、Wilson-θ法、Newmark-β法、中心差分法、Houbolt法和精细积分法;论述了各种数值积分方法的基本原理、稳定性和适用范围;通过算例,指出了各种数值积分方
法的优缺点,证实了精细积分法在计算精度和在长步长、长时间内中保持稳定性的
优越性.
【总页数】3页(P81-83)
【作者】韩爱红;钱晓军
【作者单位】华北水利水电大学,郑州450008;华北水利水电大学,郑州450008【正文语种】中文
【中图分类】TU311.3
【相关文献】
1.反应堆时空动力学方程的解法研究
2.高速光驱粘弹性阻尼减振框架结构的动力学方程及有限元数值解
3.结构动力学方程的显式与隐式数值计算
4.斜压海洋动力学
的一种三维数值模式Ⅰ.动力学方程数值格式5.分数阶导数型粘弹性结构动力学方
程的数值算法
因版权原因，仅展示原文概要，查看原文内容请购买。
威尔逊法
P ( t )t i
P i(
t )
（9）
式（6） 、式（7）和式（9）代入式（8） ，得到计算 u(ti t ) 的方程
ˆ (t t ) P ˆ t( t ) ku i i
其中
ˆk k 6 3 m c 2 (t ) t
（10）
（11）
ˆ (t t ) P ( P P ) 6 u 6 u 2u m P i i i 1 i i (t ) 2 i t i 3 t ( ui 2ui ui )c t 2
(
3 t ui 2 ui ui ) C t 2
ˆ ˆ K u (ti t ) P(ti t ) 6 6 3 ui 1 3 2 (u (t i t ) u i ui ) 2 ui (1 ) ui t t t ui 1 ui (ui 1 ui ) 2 t 2 u u t u (ui 1 2 ui ) i 1 i i 6
（16）
其中等效刚度矩阵和等效载荷向量分别为
6 3 ˆ K K (t ) 2 M t C ˆ (t t ) P (P P
i i i 1
6 6 Pi ) u ui 2 ui M 2 i t (t )
（12）
又式（10）得到 u(ti t ) 代入式（6）求得 u(ti t ) ，再把 u(ti t ) 代 入式（1） ，并取 t ，得到
i 令式（2）和式（3）中的 1 ，并取 t ，可得到 t t 时刻的速度和 3 2
u( t u i t ) i1
6 t
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
多自由度线性体系Wilson -θ法程序编写
【摘要】本文主要介绍了通过使用Matlab 软件，Wilson-θ法编写多自由度线性
体系的程序的原理、流程图、具体算例以及使用注意事项。
通过该程序可以得到剪切型结构在任意函数荷载作用下各质点的位移函数。
【关键词】Matlab ；多自由度；Wilson-θ法
1.wilson-θ法原理
wilson-θ法中最主要的步骤就是推导由t 时刻的状态求t t ∆+时刻的状态的递推公
式，现推导如下：
对τ积分
解出
代入
整理，得
其中
本程序的核心就是对以上公式的循环使用。
{}{}{}{})(t t t t t y y t
y
y -∆+=∆++θτθτ
t ∆=θτ{}{}{}{}{})(22
t t t t t t y
y t y y y
-∆++=∆++θτθττ{}{}{}{}{}{})(623
2t t t t t t t y
y t
y y y y -∆+++=∆++θτ
θτττ{}{}{}{}{})(21
t t t t t t t y
y t y t y y -∆+∆+=∆+∆+θθθθ{}{}{}{}{})2(6
)(2
t t t t t t t y
y t y t y y +∆+∆+=∆+∆+θθθθ{}{}{}{}{}t t t t t t t y y t y y t y 26
)()(62
-∆--∆=∆+∆+θθθθ{}{}{}{}{}t t t t t t t y
t
y y y t y 2
2)(3∆---∆=∆+∆+θθθθ[]{}[]{}[]{}{}t t t t t t t t P y k y C y
m ∆+∆+∆+∆+=++θθθθ []{}[]{}[]{}{}P y k y C y
m =++ []{}[]
R y k t t =∆+θ[]
[][][]
c t
m t k k ∆+∆+
=θθ3
)(6
2
[]{}{}{}[]{}{}{}[]{}{}{})223()26)(6()(2t t t t t t t t t t y
t y y t c y y t y t m P P P R ∆++∆++∆+∆+-+=∆+θθθθθ{}{}{}{})
(t t t t t t P P P P -+=∆+∆+θθ
2.程序流程图
求出各常数值
For I=1 to n
[][][][]
c a m a k k 1
++=
3.具体应用算例
如图所示，两自由度框架结构，其中
初始静止，求各层位移。
将相应的数据输入到程序中，得出各层位移关于时间的图像。
图1为第一层，图2为第二层。
将所得数值解与精确解相比较，图中实线为数值解，虚线为精确解。
由两张图，我们可以看出数值解大致是与精确相近的，但是仍然有些许的不同，这可能是算法中仍然有缺陷，说明程序仍然有待改善。
kg
1020021==m m kN/m 30001=k kN/m 20002=k kN 217.0=P s
/145.11=
θ图1
4.程序使用注意事项
（1）本程序针对于剪切型刚架结构，对于其他结构无法使用。
（2）本程序中各质点的荷载必须是函数的形式（包括常数），即对于只有某些点的荷载无法使用，且荷载函数输入时，必须采用inline语句。
例如荷载为常数10，则输入inline（’10’）；如荷载函数为sin（at），则输入inline（’sin(a*t)’,’t’）。
（3）本程序主要针对无阻尼情况，若有阻尼，只需输入阻尼矩阵即可。
（4）θ的值应大于1.37，通常取1.4，优化值为1.420815。
（5）从第一层开始为m
1，m
2
……
【参考文献】
[1] 王焕定. 结构力学（第3版）[M]. 北京：高等教育出版社，2010.
[2] Anil K.Chopra. 结构动力学理论及其在地震工程中的应用（第2版）[M]. 北京：高等教育出版社，2007.
图2。