結(jié)構(gòu)動力學(xué)newmark法程序

1、用matlab編程實(shí)現(xiàn)Newmark - B法計(jì)算多自由度體系的動力響應(yīng)姓名：王輝學(xué)號：12121190班級：碩1205班專業(yè)：橋梁與隧道工程2012年11月27日用matlab編程實(shí)現(xiàn)Newmark - B法計(jì)算多自由度體系的動力響應(yīng)一、Newmark - B法的基本原理Newmark-法是一種逐步積分的方法，避免了任何疊加的應(yīng)用，能很好的適應(yīng)非 線性的反應(yīng)分析。Newmark-:法假定：ut . t =ut (1 - ：)ut 1ut j ：t(1-1)ut t =ut u ：t (-2 - )ut ut . J(1-2)式中，1和是按積分的精度和穩(wěn)定性要求進(jìn)行調(diào)整的參數(shù)。當(dāng) =0.5,

2、=0.25時(shí), 為常平均加速度法，即假定從t到t+ t時(shí)刻的速度不變，取為常數(shù)-(ut +ut 鋰)。研究表明，當(dāng) P> 0.5 » 0.25(0.5+)2 時(shí)，Newmark-B法是一種 2無條件穩(wěn)定的格式。由式(2-141)和式(2-142 )可得到用ut ,t 及ut，ut，ut 表示的u.t，utj表達(dá)式，即有 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時(shí)刻的振動微分方程為：Mut,t Cut,t Kut,t 二Rt t(1-

3、5)將式(2-143)、式(2-144)代入(2-145)，得到關(guān)于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- 一：法的計(jì)算步驟如下：1. 初始計(jì)算：（1）形成剛度矩陣K、質(zhì)量矩陣M和阻尼矩陣C；（2）給定初始值u0，u0和u0 ；(3)選擇積分步長t、參數(shù)：、，并計(jì)算積分常數(shù)P.< 1 二 :t1-<2 ：t ：t -

4、 二云（一-2），-迸(1 )，：7 - -t ；（4）形成有效剛度矩陣K二K: °M ： 1【C；2. 對每個(gè)時(shí)間步的計(jì)算：（1）計(jì)算t+ t時(shí)刻的有效荷載：Ft * 二FtM（：。心：2ut ： 3ut）C（： 1ut ：,ut ： 5ut）（2）求解t+ t時(shí)刻的位移：'Kut,t =Ft,t.:（3）計(jì)算t+ t時(shí)刻的速度和加速度：u t傾=：0（ u t 譏一 ut ） - ： 2 u t - ： 3u tut t 二ut ： 6ut ： 7ut .tNewmark- 方法是一種無條件穩(wěn)定的隱式積分格式， 時(shí)間步長氏的大小不影響解 的穩(wěn)定性，氏的選擇主要根據(jù)解的精度

5、確定。二、本文用Newmark - B法計(jì)算的基本問題四層框架結(jié)構(gòu)在頂部受一個(gè)簡諧荷載 F二F°sin（纟t）的作用，力的作用時(shí)間t1t1 =5s，計(jì)算響應(yīng)的時(shí)間為100s，分2000步完成。阻尼矩陣由Rayleigh阻尼構(gòu)造。具體數(shù)據(jù)如下圖：mi=ikgUh=2kg77777777777/777777777/ms=3kg了/"/"/"m+Mkg#/"ki=800N.mk2=1600N.mk；=3200N.mki=6400N.ni圖一：結(jié)構(gòu)基本計(jì)算簡圖三、 計(jì)算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四、 計(jì)算結(jié)果截圖最后程序分別計(jì)算出四個(gè)質(zhì)點(diǎn)