newmark法程序法計(jì)算多自由度體系的動(dòng)力響應(yīng)

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

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

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

4、F=F. sin()的作用，力的作用時(shí)間fl =5s,計(jì)算響應(yīng)的時(shí)間為100s,分2000步完成。阻尼矩陣由Rayleigh阻尼構(gòu)造。具體數(shù)據(jù)如下圖：圖一:結(jié)構(gòu)基本計(jì)算簡圖三、計(jì)算Newmark - p法的源程序clcclearm=l,2,3,4;m=diag(m);%計(jì)算質(zhì)量矩陣k=800-80000;-8002400-16000;0；00-；%計(jì)算剛度矩陣c=2*m* sqrt(k/m);tl=5;%力的作用時(shí)間nt=2000;%分 2000 步完成dt=;%時(shí)間步長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、力作用時(shí)間內(nèi)對(duì)結(jié)構(gòu)進(jìn)行加載elsef=0:0;0;0;%力作用時(shí)間外結(jié)構(gòu)不受力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); %/+/時(shí)刻的有效荷載d(:,i)=inv(ke)*fe;%求解 t+t 時(shí)刻的位移a(:,i)=aO*(d(:,i)-d(:,i-l)-a2*v(:,i-1 )-a3*a(:,i-1);%計(jì)算 t+t 時(shí)刻的加速度v(:,i)=v(:,i-1 )+a6*a(:,i-1 )+a7*a(:,i

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