《newmark法程序法计算多自由度体系的动力响应.doc》由会员分享，可在线阅读，更多相关《newmark法程序法计算多自由度体系的动力响应.doc（14页珍藏版）》请在三一文库上搜索。
1、用matlab编程实现Newmark - B法计算多自由度体系的动力响应用matlab编程实现Newmark - B法计算多自由度体系的动力响应一、Newmark - B法的基本原理Newmark- 法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应 非线性的反应分析。Newmark-：法假定：ut t 二ut 心一 ：)utut 订 t(1-1) ut 沁二ut u ：t ( -)ut ut (1-2)2式中，1和是按积分的精度和稳定性要求进行调整的参数。当 亠0.5, =0.25时, 为常平均加速度法，即假定从t到t+ t时刻的速度不变，取为常数11 22(ut +U仏)。研究表明
2、，当 P> 0.5 丫0.25(0.5+)时，Newmark-B法是一种 无条件稳定的格式。由式(2-141)和式(2-142)可得到用ut t及ut，ut，ut表示的ut t，ut t表达式，即有. 1 1.1.ut(ut .t -ut) -一tut 一(厂-1)ut(1-3)ut t：二(ut t -ut) (1 -)ut(1 - 厂)：tut(1-4)考虑t+=t时刻的振动微分方程为:Mut t Cut . Kut t 二Rt t(1-5)将式(2-143)、式(2-144)代入(2-145),得到关于ut+ t的方程Ku t t 二Rt t(1-6)式中1 -K仆M pC 1 1
3、-1 R -Rt t MK-Mtuh (2_-1)ut)P P . P A -C(叭(-1)ut (厂-0 ：tut)求解式(2-146)可得ut t，然后由式(2-143)和式(2-144)可解出ut t和ut超。 由此，Newmark-：法的计算步骤如下：1.初始计算：(1)(2)形成刚度矩阵K、质量矩阵M和阻尼矩阵C； 给定初始值uo，(3)选择积分步长t、Uo 和Uo ；参数：、，并计算积分常数：'0十丄-1,3 2 ，(4)形成有效刚度矩阵(2)，: 6=：t(1 一 ：)，： 7 二：t ；K“K ： °M-Cl ；2.对每个时间步的计算:：4>5(1) 计
4、算t+ t时刻的有效荷载：Ft 亠二 Ft亠M (-：%ut ”：2ut 心3口t)C(： iut ： 4ut ： 5ut)(2) 求解t+t时刻的位移：Kut t -Ft ,t.:(3) 计算t+ t时刻的速度和加速度：u tdjt： = ： 0 ( ut：yt：- ut ) - ： 2 ut - ： 3u tut t 二ut ： 6ut ： 7utNewmark- 方法是一种无条件稳定的隐式积分格式，时间步长耳的大小不影响解的稳定性，氏的选择主要根据解的精度确定。二、 本文用Newmark - B法计算的基本问题四层框架结构在顶部受一个简谐荷载 F二F°sin(4 -)的作用，力
5、的作用时间t1J=5s，计算响应的时间为100s，分2000步完成。阻尼矩阵由Rayleigh阻尼构造。 具体数据如下图：mi=ikgm：=2kgITb=3kgm+=4kgki=800N.mk2=1600N.mk3=3200N.mk6400N.m图一：结构基本计算简图三、计算Newmark - B法的源程序clcclearm=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=2*m*0.05*sqrt(k/m);t仁5;%力的作用时间nt=20
6、00;%分 2000 步完成dt=0.01;%时间步长alfa=0.25;% =0.25beta=0.5;% =0.5a0=1/alfa/dt/dt;% ：0a仁beta/alfa/dt;a2=1/alfa/dt;a3=1/2/alfa-1;a4=beta/alfa-1;p%4 一 -1a5=dt/2*(beta/alfa-2);a6=dt*(1-beta);a7=dt*beta;d=zeros( 4,n t);v=zeros( 4,n t);a=zeros( 4,n t);for i=2: nt:t% ： 5 巧(一-2)% ： 6 = ；t(1 - )%_ -氏%初位移为0%初速度为0%
7、初加速度为0t=(i-1)*dt;if (t<t1),f= 300*sin(6*pi*t)-50*cos(3*pi*t);0;0;0; elsef=0;0;0;0;%力作用时间外结构不受力end%力作用时间内对结构进行加载ke=k+aO*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);% t+飢时刻的有效荷载 %求解t+ t时刻的位移%计算t+ t时刻的加速度 %计算t+ t时刻的速度d(:,i)=i nv(ke)*fe;a(:,i)=a0*(
8、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);endT=0:dt:19.99;close all figure plot(T,d) title('? + ? e p? o? x u ?') legend('?eo?,'?e?t','?xlabel(' ±e?(s)') ylabel('? o ?(m)') grid figure plot(T,a)title('? + ? e p? o?
9、u?ex u ?') legend('?eo?,'?e?t','?xlabel(' ±e?(s)')ylabel('?o? u? e(m/sA2)') grid%figure plot(T,v) title(' ?-? e p? u?ex u ?') legend('?eo?,'?e?t','?xlabel(' ±e?(s)') ylabel('? u? e(m/sA2)') gridey'：?许 ey',&
10、#39;?%离散系统dt为采样周期19.99为终端时间%控制窗口数量%绘制位移函数图像%添加标题为各质点位移总图 e?'）%添加图例的标注%对x轴进行标注为时间（s）%对y轴进行标注为位移（m）%显示画图中的个网线%绘制加速度函数图像%添加标题为各质点加速度总图e?'）%对x轴进行标注为时间（s）2%对y轴进行标注为加速度（m/s ）%绘制速度函数图像%添加标题为各质点速度总图pey'，'?回？'）%对x轴进行标注为时间（s）%对y轴进行标注为加速度（m/s2）四、计算结果截图7最后程序分别计算出四个质点的位移、速度、加速度响应现将部分截图如下：1、位移响应:图二：1质点的位移响应8图三：4质点的位移响应2、速度响应图四：1质点的速度响应图五：4质点的速度响应3、加速度响应图六：1质点的加速度响应图七：4质点的加速度响应