哈工大结构动力学作业-威尔逊-θ法

结构动力学大作业（威尔逊- 法）

姓名：

学号：

班级：

专业：

威尔逊-θ法原理及应用

【摘要】在求解单自由度体系振动方程时我们用了常加速度法及线加速度法等数值分析方法。在多自由度体系中，也有类似求解方法，即中心差分法及威尔逊-θ法。实际上后两种方法也能求解单自由度体系振动方程。对于数值方法，有三个重要要求：收敛性、稳定性及精度。本文推导了威尔逊-θ法的公式，并利用MATLAB 编程来研究单自由度体系的动力特性。 【关键词】威尔逊-θ法 冲击荷载 阻尼比

【正文】威尔逊-θ法可以很方便的求解任意荷载作用下单自由度体系振动问题。实际上，当 1.37θ>时，威尔逊-θ法是无条件收敛的。 一、威尔逊-θ法的原理

威尔逊-θ法是线性加速度法的一种拓展（当1θ=时，两者相同），其基本思路和实现方法是求出在时间段[],t t t θ+∆时刻的运动，其中1θ≥，然后通过内插得到

i t t +∆时刻的运动（见图 1.1）。

图 1.1

1、公式推导

推导由t 时刻的状态求t t θ+∆时刻的状态的递推公式：

{}{}{}{})(t t t t t y

y t

y y -∆+=∆++θτθτ

对τ积分

{}{}{}{}{})(22

t t t t t t y

y t y y y

-∆++=∆++θτθττ

{}{}{}{}{}{})(623

2

t t t t t t t y

y t y y y y -∆+++=∆++θτ

θτττ

t ∆=θτ

{}{}{}{}{})(2

1

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

22)(3∆---∆=∆+∆+θθθθ

[]{}[]{}[]{}{}P y k y C y

m =++ []{}[]{}[]{}{}t t t t t t t t P y k y C y

m ∆+∆+∆+∆+=++θθθθ

[]{}

{}t

t t

t R y k ∆+∆+=θθ

[][][][]

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 ∆++∆++∆+∆+-+=∆+θθθθθ

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('位移图形')