1、newmark法程序法計算 多自由度體系的動力響應Last updatedon 也0 afternoon of January 3, 2021用matlab編程實現Newmark - p法計算多自由度體系的動力響應用matlab編程實現Newmark0法計算多自由度體系的動力響應、Newmark - p法的基本原理Ncwtmrk-法是一種逐步積分的方法，避免了任何疊加的應用，能很好的適應非線性的反應分析。Newmark-法假定：“+=Wt +(1 0)“L +0可“4(11) =“+0 + (*-川耳 +7汁川(1-2) 式中，和是按積分的精度和穩定性要求進行調整的參數。當二，二時，為常平均加

2、速度 法，即假定從/到屮時刻的速度不變，取為常數+(+何心)。研究表明.當2, 寸)，時，Newmark-法是一種無條件穩定的格式。由式(2-141)和式(2J42)可得到用何心及“.表示的何3，何“表達 式，即有嘰廠V(也-町嗆冋-&-呵g心唸(何2一冋)+(詳)何+(】一欽阿(T考慮r+r時刻的振動微分方程為:M 耐+K” =砒(15)將式(2143)、式(2144)代入(2145),得到關于曲的方程斤何=刃z(l6)式中 求解式(2-146)可得嘰 然后由式(2-143)和式(2-144河解出和u)/+J/ o由此，Newmark-法的計算步驟如下：1 初始計算:(1) 形成剛度矩陣K、

3、質量矩陣M和阻尼矩陣C； 給定初始值,“()和可。；(3) 選擇積分步長/、參數、，并計算積分常數1_ P11,a = W 4=莎廠莎，勺=-1,務=半(-2), a6 = At( - /?), an=pt ;/ 2 /(4) 形成有效剛度矩陣R = K + %M + C；2. 對每個時間步的計算：(1) 計算/+/時刻的有效荷載：(2) 求解/+/時刻的位移：(3) 計算/+/時刻的速度和加速度：Newmark-方法是一種無條件穩定的隱式積分格式，時間步長t的大小不影響解的穩定 性，f的選擇主要根據解的精度確定。二、本文用Newmark - p法計算的基本問題四層框架結構在頂部受一個簡諧荷載

4、F=F. sin()的作用，力的作用時間fl =5s,計算響應的時間為100s,分2000步完成。阻尼矩陣由Rayleigh阻尼構造。具體數據如下圖：圖一:結構基本計算簡圖三、計算Newmark - p法的源程序clcclearm=l,2,3,4;m=diag(m);%計算質量矩陣k=800-80000;-8002400-16000;0；00-；%計算剛度矩陣c=2*m* sqrt(k/m);tl=5;%力的作用時間nt=2000;%分 2000 步完成dt=;%時間步長alfa=;%=beta=;%=aO=l/cilfa/dt/dt;%c _ 1= a 1 =beta/alfa/dt;% a

5、x = -yAta2=l/alfa/dt;%Zn = yAta3=l/2/alfa-l;%a = - 12 丫a4=beta/alfa - l;%a4 =1/a5=dt/2*(beta/alfa-2);% a5 =-(-2)2 /a6=dt*( 1 -beta);% a6 =4/(1 一 0)a7=dt*beta;%7 =04 d=zeros(4,nt);%初位移為 0 v=zeros(4,nt);%初速度為 0 a=zeros(4,nt);%初加速度為 0fori=2:ntt=(i-I)*dt;if(ttl),f=300*sin(6*pi*t)-50*cos(3*pi*t);0;0;0;%

6、力作用時間內對結構進行加載elsef=0:0;0;0;%力作用時間外結構不受力endke=k+aO*m+a 1水c;%有效剛度矩陣fe=f+m*(aO*d(:J-l)+a2*v(:,i-l)+a3*a(:J-1 )+c*(a 1 *d(:,i-1 )+a4* v(:,i-1 )+a5 *a(: J-1); %/+/時刻的有效荷載d(:,i)=inv(ke)*fe;%求解 t+t 時刻的位移a(:,i)=aO*(d(:,i)-d(:,i-l)-a2*v(:,i-1 )-a3*a(:,i-1);%計算 t+t 時刻的加速度v(:,i)=v(:,i-1 )+a6*a(:,i-1 )+a7*a(:,i

7、)；%計算 t+t 時刻的速度endT=O:dt:;%離散系統dt為采樣周期為終端時間closeallfigure%控制窗口數量plot(T,d)%繪制位移函數圖像titleC-epoxiii)%添加標題為各質點位移總圖legend(,epd,ept,epey；e|j,)% 添加圖例的標注xlabel(*e?(s)%對x軸進行標注為時間(s)ylabelCd(m)%對y軸進行標注為位移(m)grid%顯示畫圖中的個網線figureplot(T,a)%繪制加速度函數圖像titleCenouexuf)%添加標題為各質點加速度總圖 legend(,epde|jte|jey,e|j,)xlabel(e?(s)%對x軸進行標注為時間(s)ylabelC6ue(m/sA2)%對y軸進行標注為加速度(m/s2)grid%figureplot(T,v)%繪制速度函數圖像titleC-epuexuf)%添加標題為各質點速度總圖legend(epd,ente|jey,e|j,)xlabel(e?(s)%對x軸進行標注為時間(s)ylabel(ue(m/sA2)%對y軸進行標注為加速度(m/s2)grid四、計算結果截圖最后