结构动力学大作业（威尔逊- 法）
：
学号：
班级：
专业：
威尔逊-θ法原理及应用
【摘要】在求解单自由度体系振动方程时我们用了常加速度法及线加速度法等数值分析方法。
在多自由度体系中，也有类似求解方法，即中心差分法及威尔逊-θ法。
实际上后两种方法也能求解单自由度体系振动方程。
对于数值方法，有三个重要要求：收敛性、稳定性及精度。
本文推导了威尔逊-θ法的公式，并利用MATLAB 编程来研究单自由度体系的动力特性。
【关键词】威尔逊-θ法 冲击荷载 阻尼比
【正文】威尔逊-θ法可以很方便的求解任意荷载作用下单自由度体系振动问题。
实际上，当 1.37θ>时，威尔逊-θ法是无条件收敛的。
一、威尔逊-θ法的原理
威尔逊-θ法是线性加速度法的一种拓展（当1θ=时，两者相同），其基本思路和实现方法是求出在时间段[],t t t θ+∆时刻的运动，其中1θ≥，然后通过插得到i t t +∆时刻的运动（见图 1.1）。
图 1.1
1、公式推导
推导由t 时刻的状态求t t θ+∆时刻的状态的递推公式：
对τ积分
{}{}{}{}{}{})(623
2
t t t t t t t y
y t y y y y &&&&&&&-∆+++=∆++θτ
θτττ
{}{}{}{}{})2(6)(2t 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-∆--∆=∆+∆+θθθθ
[]{}{}
{}[]{}{}{}[]{}{}{})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 &&&&&&∆++∆++∆+∆+-+=∆+θθθθθ
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>Lp
P(Lp+1:Lu)=0; end
elseif isnumeric(Pt)%荷载为散点 if Lu<=Lp
P=Pt(1:Lu);
else
P(1:Lp)=Pt;
P(Lp+1:Lu)=0;
end
end
y=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-1
kk=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);
end
plot(uds,y,'r'),xlabel('时间t'),ylabel('位移y'),title('位移图形')
二、利用威尔逊-θ法求冲击荷载下的结构反应
1、矩形脉冲
研究不同时长脉冲作用下，体系振动位移。
取单自由度刚度为1N/m，质量为1/(4*pi^2)kg，频率为2*pi1s-，周期为1s，阻尼c=0，荷载为1N，积分步长为0.1，θ=1.42，初位移为0，初速度为0时的质点位移时间图如下：
图2.1 1/41/4d n t s T ==
图2.2 1/21/2d n t s T ==
图2.3 1d n t s T ==
图2.4 1.25 1.25d n t s T ==
图2.5 1.5 1.5d n t s T ==
由图形可看出：当/2d n t T >时，最大位移发生在荷载离开前； 当/2d n t T <时，最大位移发生在荷载离开后； 当/2d n t T =时，最大位移发生在荷载离开时。
特别的，当d n t T =时，d t t =后没有位移。
2、其他脉冲
图2.6 负斜率直线
图2.7 正斜率直线
图2.8 二次抛物线
图2.9 5次抛物线
三、利用威尔逊-θ法求不同阻尼下结构自振反应
本体系度刚度为1N/m ，质量为1/(4*pi^2)kg ，频率为2*pi 1s -。
故其临界阻尼
21/r c m ω==pi 。
分别取结构阻尼c 为：0.05/pi ，0.1/pi ，1/pi ，1.5/pi ，进行计算。
计算结果见下图：
图3.1 0.05ξ=
ξ=
图3.2 0.1
ξ=图3.3 1
图3.4 1.5ξ=
由图形对比可知，当1ξ<时，阻尼越大，结构运动衰减越快，此时结构处于小阻尼状态；当1ξ=体系不能振动，此时的c 为临界阻尼；当1ξ>时，为超阻尼系统，不发生自由振动。
四、利用威尔逊-θ法例
已知：如下图，W=438.18kN;k=40181.1kN/m
求：体系位移反应。
图4.1 结构 图4.2 荷载
由MATLAB 计算出的结果见图 4.3。
图4.3 位移曲线
手算结果：
kN 1088.433⨯==g W m
s /119.301088.4310437=⨯⨯=
ω s
208.02==ωπ
T
T t <<=05.01
3310.02543010 5.375102
S =⨯⨯⨯=⨯ 3
335.37510 4.0571043.881030.19
S A m m ω-⨯===⨯⨯⨯ 手算结果与电算结果近似相等，说明当冲击荷载作用时间远小于结构自振周期是，可近似认为最大位移S A m ω
=（S 为荷载与坐标轴所围成的面积）。