结构动力学大作业.docx

1、结构动力学大作业结构动力学大作业班级： 学号：姓名： 1. Wilson-法原理简介图1-1 Wilson-法示意图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-2 W

2、ilson-法流程图2. Wilson-程序验算对线弹性条件下的Wilson-法进行MATLAB编程，源代码见附录。选取如下算例进行验证。对于一个单自由度的无阻尼结构，当其受到一个周期荷载时，其结构响应分为稳态解和瞬态解，由于没有阻尼的影响，其瞬态解并不会衰减，其理论表达式为： (2-1)式中，为位移响应，为激励，为刚度，为荷载频率与固有振动频率之比，为荷载频率，为结构固有频率。现令为1，为1，则为1，取为2/3。程序求得的解与解析解对比如图2-1所示（由于理论解与程序基本重合，所以将理论解乘以-1，方便比较）：a)位移b)速度c)加速度图2-1 Wilson-法结果验证2.1 t的影响上述算

3、例验证时选择的t非常小，因此看不出理论解与Wilson-法的求解区别，以下改变t的取值，探讨t对迭代的影响。图2-2 t对位移曲线的影响可以看出并不是t太大时计算结果很不准确，偏小，反映不出周期特征；当t合适时正好基本和理论解重合，也不是t越小越好越小时越能反映出一些细部特征，但这也不是很准确。2.2 的影响当1.37时，该算法是无条件稳定的算法，以下探讨对算法的影响。图2-3 对位移曲线的影响由上图可知随着值越大，位移的周期变大。3. 非线性问题求解由于实际结构并不一定为线性，其刚度会随着位移的的变化而改变，下图为求解非线性问题时的Wilson-法流程。此处要说明的是，刚度矩阵Ky(t)是与

4、位移相关的量，判断那时候速度的大小是为了确定其是否处于卸载段。具体可能得根据实际情况求解。图3-1 Wilson-法解非线性问题修改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极

5、值在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代表加速度if nargin9 theta=1.4;endL,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_

6、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=dt2/6;%计算拟刚度矩阵k0=k+a0*m+a1*c;%计算拟荷载for i=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=k0R; y_3(:,i+1)=a4*(y_th

7、eta-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代表输入

8、的荷载，c为阻尼矩阵，dt为时间间隔,m为质量矩阵%v0为初始的速度，y0为初始的位移.a_0为初始加速度%输出的矩阵y_1代表位移,y_2代表速度，y_3代表加速度if nargin0 k=60*(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=k0R; 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