newmark法程序法計算多自由度體系的動力響應

1、用matlab編程實現(xiàn)法計算多自由度體系的動力響應用matlab編程實現(xiàn)法計算多自由度體系的動力響應一、法的基本原理Newmark-b法是一種逐步積分的方法，避免了任何疊加的應用，能很好的適應非線性的反應分析。Newmark-b法假定： (1-1) (1-2)式中，b和g是按積分的精度和穩(wěn)定性要求進行調(diào)整的參數(shù)。當b=0.5，g=0.25時，為常平均加速度法，即假定從t到t+Dt時刻的速度不變，取為常數(shù)。研究表明，當b0.5， g0.25(0.5+b)2時，Newmark-b法是一種無條件穩(wěn)定的格式。由式(2-141)和式(2-142)可得到用及，表示的，表達式，即有 (1-3) (1-4)考

2、慮t+Dt時刻的振動微分方程為： (1-5)將式(2-143)、式(2-144) 代入(2-145)，得到關(guān)于ut+Dt的方程 (1-6)式中求解式(2-146)可得，然后由式(2-143)和式(2-144)可解出和。由此，Newmark-b法的計算步驟如下：1.初始計算：（1）形成剛度矩陣K、質(zhì)量矩陣M和阻尼矩陣C；（2）給定初始值， 和；（3）選擇積分步長Dt、參數(shù)b、g，并計算積分常數(shù)，；（4）形成有效剛度矩陣；2.對每個時間步的計算：（1）計算t+Dt時刻的有效荷載：（2）求解t+Dt時刻的位移： （3）計算t+Dt時刻的速度和加速度：Newmark-b方法是一種無條件穩(wěn)定的隱式積分格

3、式，時間步長Dt的大小不影響解的穩(wěn)定性，Dt的選擇主要根據(jù)解的精度確定。二、 本文用法計算的基本問題四層框架結(jié)構(gòu)在頂部受一個簡諧荷載的作用，力的作用時間=5s，計算響應的時間為100s，分2000步完成。阻尼矩陣由Rayleigh阻尼構(gòu)造。具體數(shù)據(jù)如下圖：圖一：結(jié)構(gòu)基本計算簡圖三、 計算法的源程序clcclearm=1,2,3,4; m=diag(m); %計算質(zhì)量矩陣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) ; t1=5; %力的作用時間

4、 nt=2000; %分2000步完成 dt=0.01; %時間步長 alfa=0.25; %g=0.25 beta=0.5; %b=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); %初位移為0 v=zeros(4,nt); % 初速度為0 a=zeros(4,nt); % 初加速度為0 for i=2:nt t=(i-1)*

5、dt; if (t<t1), f=300*sin(6*pi*t)-50*cos(3*pi*t);0;0;0; %力作用時間內(nèi)對結(jié)構(gòu)進行加載 else f=0;0;0;0; %力作用時間外結(jié)構(gòu)不受力 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); % t+Dt時刻的有效荷載 d(:,i)=inv(ke)*fe; %求解t+Dt時刻的位移 a(:,i)=a0*(d(:,i)-d(:,i-1)-a2*v(:,i-1)-

6、a3*a(:,i-1); %計算t+Dt時刻的加速度 v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i); %計算t+Dt時刻的速度 endT=0:dt:19.99; %離散系統(tǒng)dt為采樣周期 19.99為終端時間close allfigure %控制窗口數(shù)量plot(T,d) %繪制位移函數(shù)圖像title('¸÷ÖʵãλÒÆ×Üͼ') %添加標題為各質(zhì)點位移總圖 legend('ÖÊ

7、;µãÒ»','Öʵã¶þ','ÖʵãÈý','ÖʵãËÄ') %添加圖例的標注xlabel('ʱ¼ä(s)') %對x軸進行標注為時間（s）ylabel('λÒÆ(m)') %對y軸進行標注為位移（m

8、）grid %顯示畫圖中的個網(wǎng)線figure plot(T,a) %繪制加速度函數(shù)圖像title('¸÷Öʵã¼ÓËÙ¶È×Üͼ') %添加標題為各質(zhì)點加速度總圖legend('ÖʵãÒ»','Öʵã¶þ','Öʵã

9、;Èý','ÖʵãËÄ')xlabel('ʱ¼ä(s)') %對x軸進行標注為時間（s）ylabel('¼ÓËÙ¶È(m/s2)') %對y軸進行標注為加速度（m/s2）grid%figure plot(T,v) %繪制速度函數(shù)圖像title('¸÷ÖʵãËÙ¶

10、;È×Üͼ') %添加標題為各質(zhì)點速度總圖legend('ÖʵãÒ»','Öʵã¶þ','ÖʵãÈý','ÖʵãËÄ')xlabel('ʱ¼ä(s)') %對x軸進行標注為時間（s）ylabel('ËÙ¶È(m/s2)&