用matlab編程實現(xiàn)法計算多自由度體系的動力響應

1、用matlab編程實現(xiàn)Newmark - B法計算多自由度體系的動力響應姓名: 學號: 班級: 專業(yè):用matlab編程實現(xiàn)Newmark - B法計算多自由度體系的動力響應、Newmark - B法的基本原理Newmark-法是一種逐步積分的方法，避免了任何疊加的應用，能很好的適應非 線性的反應分析。Newmark-:法假定：ut t 二ut 心 - ：)ut ：ut A t(1-1)ut ,t =ut uf '4 (-2 - )ut ut. t ：t2(1-2)式中，1和是按積分的精度和穩(wěn)定性要求進行調整的參數(shù)。當上0.5, =0.25時,為常平均加速度法，即假定從t到t+ t時刻

2、的速度不變，取為常數(shù) -(ut +ut植)。研究表明，當 P> 0.5 丫0.25(0.5+)2 時，Newmark法是一種2無條件穩(wěn)定的格式。由式(2-141)和式(2-142)可得到用ut t及ut，ut，ut表示的ut t，ut t 表達式，即有 1 1-1 ut 二tF(ut ：莫 ut)ut -_1)u) t(1-3)pp .P .ut t 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的方程

3、Kut t 二Rt t(1-6)式中1 PK=K十M石C1 1 1R二Rt t M(ut 石ut (廠-1)uJP P . P A .C( tut (1)ut (21) tut)求解式（2-146）可得ut . .t，然后由式（2-143）和式（2-144）可解出ut . t：和ut. ?由此，Newmark- 一：法的計算步驟如下：1. 初始計算：（1）形成剛度矩陣K、質量矩陣M和阻尼矩陣C；(2)給定初始值u0，u0和u0 ；(3)選擇積分步長t、，并計算積分常數(shù)-<2 : ：t52(2)，二：t(1. )，JS 7 八 t ；（4）形成有效剛度矩陣K二K *0M C；2. 對每個時

4、間步的計算：（1）計算t+ 1時刻的有效荷載：Ft t=FtM（：out Mt ：3ut）C（： 1ut ：,ut ： 5ut）（2）求解t+ t時刻的位移：Kut-Ft ,t-:（3）計算t+ t時刻的速度和加速度：ut t = ： 0（ut t -ut） - ： 2ut - ： 3utut t 二ut ： 6ut ： 7ut tNewmark- 方法是一種無條件穩(wěn)定的隱式積分格式，時間步長氏的大小不影響解的穩(wěn)定性，氏的選擇主要根據(jù)解的精度確定。二、 本文用Newmark - B法計算的基本問題四層框架結構在頂部受一個簡諧荷載 F二F°sin（4 I）的作用，力的作用時間t1J=5

5、s，計算響應的時間為100s，分2000步完成。阻尼矩陣由Rayleigh阻尼構造。具體數(shù)據(jù)如下圖：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 -800 0 0;-800 2400 -1600 0;0 -1600 4800 -3200; 0 0 -3200 8000;c=0.05*m+0.02

6、*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 i=2:nt t=(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四、計算結果截圖最后程序分別計算出四個質點的位