结构动力学newmark法程序

1、用matlab编程实现Newmark - B法计算多自由度体系的动力响应姓名：王辉学号：12121190班级：硕1205班专业：桥梁与隧道工程2012年11月27日用matlab编程实现Newmark - B法计算多自由度体系的动力响应一、Newmark - B法的基本原理Newmark-法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非 线性的反应分析。Newmark-:法假定：ut . t =ut (1 - ：)ut 1ut j ：t(1-1)ut t =ut u ：t (-2 - )ut ut . J(1-2)式中，1和是按积分的精度和稳定性要求进行调整的参数。当 =0.5,

2、=0.25时, 为常平均加速度法，即假定从t到t+ t时刻的速度不变，取为常数-(ut +ut 锂)。研究表明，当 P> 0.5 » 0.25(0.5+)2 时，Newmark-B法是一种 2无条件稳定的格式。由式(2-141)和式(2-142 )可得到用ut ,t 及ut，ut，ut 表示的u.t，utj表达式，即有 1 1-1 ut .tF(ut . t -u J ut - (21)u t(1-3)pp .P .ut t t (ut t _ut) (1)ut (1 -厂):tut(1-4)考虑t+=t时刻的振动微分方程为：Mut,t Cut,t Kut,t 二Rt t(1-

3、5)将式(2-143)、式(2-144)代入(2-145)，得到关于ut+ -t的方程Kut 沁二Rt t(1-6)式中1K=K十M二C1 11R二Rt .t M(ut 二ut (-1)ut)P P . P .C( ut (1)ut ( -1) tut)求解式（2-146）可得ut . t，然后由式（2-143）和式（2-144）可解出ut ,t和ut由此，Newmark- 一：法的计算步骤如下：1. 初始计算：（1）形成刚度矩阵K、质量矩阵M和阻尼矩阵C；（2）给定初始值u0，u0和u0 ；(3)选择积分步长t、参数：、，并计算积分常数P.< 1 二 :t1-<2 ：t ：t -

4、 二云（一-2），-迸(1 )，：7 - -t ；（4）形成有效刚度矩阵K二K: °M ： 1【C；2. 对每个时间步的计算：（1）计算t+ t时刻的有效荷载：Ft * 二FtM（：。心：2ut ： 3ut）C（： 1ut ：,ut ： 5ut）（2）求解t+ t时刻的位移：'Kut,t =Ft,t.:（3）计算t+ t时刻的速度和加速度：u t倾=：0（ u t 讥一 ut ） - ： 2 u t - ： 3u tut t 二ut ： 6ut ： 7ut .tNewmark- 方法是一种无条件稳定的隐式积分格式， 时间步长氏的大小不影响解 的稳定性，氏的选择主要根据解的精度

5、确定。二、本文用Newmark - B法计算的基本问题四层框架结构在顶部受一个简谐荷载 F二F°sin（纟t）的作用，力的作用时间t1t1 =5s，计算响应的时间为100s，分2000步完成。阻尼矩阵由Rayleigh阻尼构造。具体数据如下图：mi=ikgUh=2kg77777777777/777777777/ms=3kg了/"/"/"m+Mkg#/"ki=800N.mk2=1600N.mk；=3200N.mki=6400N.ni图一：结构基本计算简图三、 计算Newmark - B法的源程序m=1,2,3,4;m=diag(m);k= 800

6、 -800 0 0;-800 2400 -1600 0;0 -1600 4800 -3200;0 0 -3200 8000;c=0.05*m+0.02*k;f0=100;t仁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*(beta/alfa-2); a6=dt*(1-beta); a7=dt*beta;d=zeros(4,nt);v=zeros(4,nt);a=zeros(4,nt);for

7、i=2:ntt=(i-1)*dt;if (t<t1), f=f0*sin(4*pi*t/t1);0;0;0; else f=0;0;0;0; end ke=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四、 计算结果截图最后程序分别计算出四个质点