结构动力学大作业.docx
《结构动力学大作业.docx》由会员分享，可在线阅读，更多相关《结构动力学大作业.docx（9页珍藏版）》请在冰豆网上搜索。
结构动力学大作业
结构动力学大作业
班级：
学号：
姓名：
1.Wilson-θ法原理简介
图1-1Wilson-θ法示意图
Wilson-θ法是基于对加速度a的插值近似得到的，图1-1为Wilson-θ法的原理示意图。
推导由t时刻的状态求t+△t时刻的状态的递推公式：
（1-1）
对τ积分可得速度与位移的表达式如下：
（1-2）
（1-3）
其中τ=θt，由式（1-2）、（1-3）可以解出：
（1-4）
（1-5）
将式（1-4）、（1-5）带入运动方程：
（1-6）
（1-7）
注意到此时的式子为{
}和上一个时刻
、
、
以及t+θ△t时刻的荷载
相关，可以运用迭代的思想来求解，下图给出线弹性条件下Wilson-θ法的流程图：
图1-2Wilson-θ法流程图
2.Wilson-θ程序验算
对线弹性条件下的Wilson-θ法进行MATLAB编程，源代码见附录。
选取如下算例进行验证。
对于一个单自由度的无阻尼结构，当其受到一个周期荷载时，其结构响应分为稳态解和瞬态解，由于没有阻尼的影响，其瞬态解并不会衰减，其理论表达式为：
（2-1）
式中，
为位移响应，
为激励，
为刚度，
为荷载频率与固有振动频率之比，
为荷载频率，
为结构固有频率。
现令
为1，
为1，则
为1，
取为2/3。
程序求得的解与解析解对比如图2-1所示（由于理论解与程序基本重合，所以将理论解乘以-1，方便比较）：
a）位移
b）速度
c）加速度
图2-1Wilson-θ法结果验证
2.1△t的影响
上述算例验证时选择的△t非常小，因此看不出理论解与Wilson-θ法的求解区别，以下改变△t的取值，探讨△t对迭代的影响。
图2-2△t对位移曲线的影响
可以看出并不是△t太大时计算结果很不准确，偏小，反映不出周期特征；当△t合适时正好基本和理论解重合，也不是△t越小越好越小时越能反映出一些细部特征，但这也不是很准确。
2.2θ的影响
当θ>1.37时，该算法是无条件稳定的算法，以下探讨θ对算法的影响。
图2-3θ对位移曲线的影响
由上图可知随着θ值越大，位移的周期变大。
3.非线性问题求解
由于实际结构并不一定为线性，其刚度会随着位移的的变化而改变，下图为求解非线性问题时的Wilson-θ法流程。
此处要说明的是，刚度矩阵[K]y（t）是与位移相关的量，判断那时候速度的大小是为了确定其是否处于卸载段。
具体可能得根据实际情况求解。
图3-1Wilson-θ法解非线性问题
修改MATLAB程序，并用该程序来计算如下例题：
对该问题采用Wilson-θ法非线性方式计算，采用△t=0.1s和△t=0.05s两种方式，计算位移、速度和加速度曲线如下图所示：
a）位移
b）速度
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
ify_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
升级会员 