用matlab编程实现法计算多自由度体系的动力响应

1、用matlab编程实现newmark - b法计算多自由度体系的动力响应姓名:学号: 班级:专业:用matlab编程实现newmark - b法计算多自由度体系的动力响应一、newmark - b法的基本原理newmark-法是一种逐步积分的方法，避免了任何叠加的应用,能很好的适应非 线性的反应分析。newmark-法假定:ut t ut (1)ut ut t t(1-1)1 2z/ zut t ut u t (-)ut ut t t(1-2)2式中，和是按积分的精度和稳定性要求进行调整的参数。当=0.5, =0.25 时，为常平均加速度法，即假定从t到t+ t时刻的速度不变，取为常数12 -

2、(ut ut t)。研究表明，当 0.5 ,0.2 5(0.5+)时，newmark-2法是一种无条件稳定的格式。由式(2-141)和式(2-142)可得到用ut t 及ut,ut,ut 表示的ut t,ut t表达式，即有111ut t r(ut t ut) ut (- 1)ut(1-3)t2t 2ut t 一 (ut t ut) (1 -)ut (1 一)tut(1-4)t2考虑t + t时刻的振动微分方程为：mut t cut t kut t rt t(1-5)将式(2-143)、式(2-144)代入(2-145),得到关于ut+ t的方程kut t rt t(1-6)式中1-k k 中

3、m c1,、1,、1r rt t m(下ut -ut (2- 1)ut)c( ut ( 1)ut (- 1) tut)可解出ut t和求解式(2-146)可得ut t,然后由式(2-143)和式(2-144) ut t。由此，newmark-法的计算步骤如下：1 .初始计算：(1)形成刚度矩阵k、质量矩阵m和阻尼矩阵q ;(2)给定初始值u。，u0和u。;(3)选择积分步长t、参数,并计算积分常数3d，2)t(1),7(4)形成有效刚度矩阵k k0m 1c;2 .对每个时间步的计算：(1)计算t+ t时刻的有效荷载：ft t ft t m( 0ut2ut3ut)c( 1ut4ut5ut)(2)

4、求解t+ t时刻的位移： kut t ft t(3)计算t+ t时刻的速度和加速度：ut t 0(ut t ut)2ut3utut t ut6ut7ut tnewmark-方法是一种无条件稳定的隐式积分格式，时间步长 t的大小不影响解的稳定性，t的选择主要根据解的精度确定。本文用newmark - b法计算的基本问题四层框架结构在顶部受一个简谐荷载 f=f0sin(42)的作用，力的作用时间t1t1二5s,计算响应的时间为100s,分2000步完成。阻尼矩阵由rayleigh阻尼构造具体数据如下图：nii=ik2mi=2kgki=800n.mm；=3kg匕二1600n.mttzzzzzzzzz

5、zzzzzzzzzzk3=3200n.mk产6400n.ni图一：结构基本计算简图三、 计算newmark - b法的源程序m=1,2,3,4;m=diag(m);k= 800 -800 0 0;-800 2400 -1600 0;0 -1600 4800 -3200;0 0 -3200 8000;c=0.05*m+0.02*k;f0=100;t1=5;nt=2000;dt=0.01;alfa=0.25;beta=0.5;a0=1/alfa/dt/dt; a1=beta/alfa/dt; a2=1/alfa/dt;a3=1/2/alfa-1; a4=beta/alfa-1;a5=dt/2*(b

6、eta/alfa-2);a6=dt*(1-beta);a7=dt*beta;d=zeros(4,nt);v=zeros(4,nt);a=zeros(4,nt);for i=2:ntt=(i-1)*dt;if (tt1), f=f0*sin(4*pi*t/t1);0;0;0; else f=0;0;0;0; endke=k+a0*m+a1*c;fe=f+m*(a0*d(:,i-1)+a2*v(:,i-1)+a3*a(:,i-1)+c*(a1*d(:,i-1)+a4*v(:,i-1)+a5*a(:,i-1);d(:,i)=inv(ke)*fe;a(:,i)=a0*(d(:,i)-d(:,i-1)-a2*v(:,i-1)-a3*a(:,i-1);v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);end计算结果截图最后程序分别计算出四个质点的位移、速度、加速度响应。现将部分截图如下：1、位移响应:图二：1质点的位移响应图三：4质点的位移响应2、速度响应图四：1质点的速度响应图