结构动力学大作业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 t
t t t y y y y t
τθτ
θ++∆=+-∆ (1-1)
对τ积分可得速度与位移的表达式如下：
{}{}{}{}{}2
()2t t t t t t y
y y y y
t
τθττθ++∆=++-∆ (1-2)
{}{}{}{}{}{}2
3
()26t t t t t t t y y y y y y
t
τ
θτττθ++∆=+++-∆ (1-3)
其中τ=θt ，由式(1-2)、(1-3)可以解出：
{}{}{}{}{}266
()2()t t t t
t t t y y y y y t t
θθθθ+∆+∆=---∆∆
(1-4)
{}{}{}{}{}3()22
t t t t t t t t
y
y 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 t
m 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编程，源代码见附录。

选取如下算例进行验证。

对于一个单自由度的无阻尼结构，当其受到一个周期荷载时，其结构响应分为稳态解和瞬态解，由于没有阻尼的影响，其瞬态解并不会衰减，其理论表达式为：
021
()()(sin sin )1p x t t t k ωβωβ
=
-- (2-1)
式中，()x t 为位移响应，0p 为激励，k 为刚度，β为荷载频率与固有振动频率之比，ω为荷载频率，ω为结构固有频率。

现令0p 为1，k 为1，则ω为1，ω取为2/3。

程序求得的解与解析解对比如图2-1所示（由于理论解与程序基本重合，所以将理论解乘以-1，方便比较）：
位移y
时间t a )位移
速度v
时间t
b )速度
加速度a
时间t
c )加速度
图2-1Wilson -θ法结果验证
2.1 △t 的影响
上述算例验证时选择的△t 非常小，因此看不出理论解与Wilson -θ法的求解区别，以下改变△t 的取值，探讨△t 对迭代的影响。

位移y
时间t
图2-2 △t 对位移曲线的影响
可以看出并不是△t 太大时计算结果很不准确，偏小，
反映不出周期特征；当△t 合适时正好基本和理论解重合，也不是△t 越小越好越小时越能反映出一些细部特征，但这也不是很准确。

2.2θ的影响
当θ>1.37时，该算法是无条件稳定的算法，以下探讨θ对算法的影响。

位移y
时间t
图2-3θ对位移曲线的影响
由上图可知随着θ值越大，位移的周期变大。

3. 非线性问题求解
由于实际结构并不一定为线性，其刚度会随着位移的的变化而改变，下图为求解非线性问题时的Wilson -θ法流程。

此处要说明的是，刚度矩阵[K ]y(t)是与位移相关的量，判断那时候速度的大小是为了确定其是否处于卸载段。

具体可能得根据实际情况求解。

图3-1Wilson-θ法解非线性问题修改MATLAB程序，并用该程序来计算如下例题：
对该问题采用Wilson-θ法非线性方式计算，采用△t=0.1s和△t=0.05s两种方式，计算位移、速度和加速度曲线如下图所示：
位
移y
时间t
a)位移
速
度
u
时间t
b)速度
加
速
度a
时间t
c)加速度
图3-2 非线性分析结果
由上图可知，结构在0.6s时达到位移极值，在△t=0.1s和△t=0.05s算得的值分别为0.096m和0.108m，速度极值在0.9s取到分别为-0.468和-0.580，加速度极值在△t=0.1s时为0.7s时取到，为-2.127，在△t=0.05s极值在0.75s时取到，为-2.9。

4.附录
Wilson-θ法源程序
function [y_1,y_2,y_3]=wilson_theta(p,m,c,k,dt,v0,y0,a_0,theta)
%p代表输入的荷载，c为阻尼矩阵，dt为时间间隔,m为质量矩阵,k为刚度矩阵
%v0为初始的速度，y0为初始的位移.a_0为初始加速度
%输出的矩阵y_1代表位移,y_2代表速度，y_3代表加速度
ifnargin<9
theta=1.4;
end
[L,r]=size(p);
y_1=NaN(L,r);y_2=NaN(L,r);y_3=NaN(L,r);
y_1(:,1)=y0;y_2(:,1)=v0;y_3(:,1)=a_0;
%计算积分常数
a0=6/((theta*dt)^2);a1=3/theta/dt;a2=2*a1;
a3=theta*dt/2;a4=a0/theta;a5=-a2/theta;
a6=1-3/theta;a7=dt/2;a8=dt^2/6;
%计算拟刚度矩阵
k0=k+a0*m+a1*c;
%计算拟荷载
fori=1:r-1
R=p(:,i)+theta*(p(:,i+1)-p(:,i))+m*(a0*y_1(:,i)+a2*y_2(:,i)+2*y_3(:,i))+c*(a1*y_1(:,i)+2*y_2(:,i)+a3 *y_3(:,i));
y_theta=k0\R;
y_3(:,i+1)=a4*(y_theta-y_1(:,i))+a5*y_2(:,i)+a6*y_3(:,i);
y_2(:,i+1)=y_2(:,i)+a7*(y_3(:,i)+y_3(:,i+1));
y_1(:,i+1)=y_1(:,i)+y_2(:,i)*dt+a8*(y_3(:,i+1)+2*y_3(:,i));
end
上述代码只适合分析线弹性结构，对于非线性结构，编写起来比较繁琐，针对不同的情况可能需要具体处理，所以本文只给出了针对本文例题的代码。

与上述代码不同的是以下代码增加了一个判断的语句。

function [y_1,y_2,y_3]=wilson_theta2(p,m,c,dt,v0,y0,a_0,theta)
%p代表输入的荷载，c为阻尼矩阵，dt为时间间隔,m为质量矩阵
%v0为初始的速度，y0为初始的位移.a_0为初始加速度
%输出的矩阵y_1代表位移,y_2代表速度，y_3代表加速度
ifnargin<9
theta=1.4;
end
[L,r]=size(p);
y_1=NaN(L,r);y_2=NaN(L,r);y_3=NaN(L,r);
y_1(:,1)=y0;y_2(:,1)=v0;y_3(:,1)=a_0;
%计算积分常数
a0=6/((theta*dt)^2);a1=3/theta/dt;a2=2*a1;
a3=theta*dt/2;a4=a0/theta;a5=-a2/theta;
a6=1-3/theta;a7=dt/2;a8=dt^2/6;
fori=1:r-1
if y_2(:,i)>0
k=60*(y_1(:,i)<=0.05)+3/y_1(:,i)*(y_1(:,i)>0.05);
else
k=60;
end
%计算拟刚度矩阵
k0=k+a0*m+a1*c;
%计算拟荷载
R=p(:,i)+theta*(p(:,i+1)-p(:,i))+m*(a0*y_1(:,i)+a2*y_2(:,i)+2*y_3(:,i))+c*(a1*y_1(:,i)+2*y_2(:,i)+a3 *y_3(:,i));
y_theta=k0\R;
y_3(:,i+1)=a4*(y_theta-y_1(:,i))+a5*y_2(:,i)+a6*y_3(:,i);
y_2(:,i+1)=y_2(:,i)+a7*(y_3(:,i)+y_3(:,i+1));
y_1(:,i+1)=y_1(:,i)+y_2(:,i)*dt+a8*(y_3(:,i+1)+2*y_3(:,i));
end。