基于MATLAB的四層框架結(jié)構(gòu)動(dòng)力響應(yīng)與研究

1、暨南大學(xué)研究生課程論文課程：結(jié)構(gòu)動(dòng)力學(xué)姓名：許可悅學(xué)號(hào)：1634361002學(xué)院：力學(xué)與建筑工程學(xué)院專業(yè)：建筑與土木工程任課教師：李雪艷基于MATLAB的四層框架結(jié)構(gòu)動(dòng)力響應(yīng)與研究許可悅（暨南大學(xué)理工學(xué)院力學(xué)與土木工程學(xué)院，廣州 51063）摘要：本文用MATLAB語言對(duì)四層建筑結(jié)構(gòu)進(jìn)行編程，計(jì)算結(jié)構(gòu)的自振頻率、振型，分析該結(jié)構(gòu)在自由振動(dòng)和一般激勵(lì)下的動(dòng)力響應(yīng)。采用了Newmark-法計(jì)算了在簡諧正弦激勵(lì)作用下結(jié)構(gòu)的位移響應(yīng)，并以此為初始條件結(jié)合瑞利阻尼矩陣計(jì)算了結(jié)構(gòu)在簡諧正弦荷載卸載后的結(jié)構(gòu)自由振動(dòng)的位移響應(yīng)。關(guān)鍵詞：MATLAB、Newmark-法、瑞利阻尼矩The four layers

2、 of frame structure dynamic response based on MATLAB and researchXu Keyue (Jinan university institute of mechanics and civil engineering department, Guangzhou)Abstract：This paper uses MATLAB language to program the the four layers of frame structure , calculates the self-vibration frequency and vibr

3、ation mode of the structure, and analyzes the dynamic response of the structure under free vibration and general excitation. Adopted the Newmark - beta method to calculate the displacement of the structure under the action of a harmonic sine excitation response, and the initial conditions in combina

4、tion with the Rayleigh damping matrix to calculate the structure in the structure of harmonic sine load after unloading free vibration displacement response.Key words：MATLAB; Newmark-method;Rayleigh orthogonal damping1引言在社會(huì)發(fā)展的今天，很多科技人員都會(huì)遇到數(shù)值分析計(jì)算機(jī)應(yīng)用等問題，一些傳統(tǒng)的高級(jí)程序語言如 FORTRAN 等雖然能在一定程度上減輕計(jì)算量 ,但它們要求應(yīng)用人員要

5、具有較強(qiáng)的編程能力和對(duì)算法有深入的研究. 另外 ,在運(yùn)用這些高級(jí)程序語言進(jìn)行計(jì)算結(jié)果的可視化分析及圖形處理方面 ,對(duì)非計(jì)算機(jī)專業(yè)的普通用戶來說 ,存在著很大的難度. MATLAB 正是在這一應(yīng)用要求背景下產(chǎn)生的數(shù)學(xué)類科技應(yīng)用軟件。MATLAB 是是以矩陣計(jì)算為基礎(chǔ)的程序設(shè)計(jì)語言 , MATLAB 具有功能豐富和完備的數(shù)學(xué)函數(shù)庫及工具箱 ,大量繁雜的數(shù)學(xué)運(yùn)算和分析可通過調(diào)用 MATLAB 函數(shù)直接求解 ,大大提高了編程效率 ,其程序編譯和執(zhí)行速度遠(yuǎn)遠(yuǎn)超過了傳統(tǒng)的 FORTRAN 語言 ,因而用 MA T2LAB 編寫程序 , 往往可以達(dá)到事半功倍的效果. 在圖形處理方面 ,MA TLAB 可以給

6、數(shù)據(jù)以二維、三維乃至四維的直觀表現(xiàn) ,并在圖形色彩、視角、品性等方面具有較強(qiáng)的渲染和控制能力 ,使科技人員對(duì)大量原始數(shù)據(jù)的分析變得輕松和得心應(yīng)手，從根本上滿足了科技人員對(duì)工程數(shù)學(xué)計(jì)算的要求 ,將科技人員及普通用戶從繁重的數(shù)學(xué)運(yùn)算中解放出來。本文通過實(shí)例介紹了 MATLAB 語言在結(jié)構(gòu)動(dòng)力學(xué)中的應(yīng)用 ,通過結(jié)構(gòu)的自振頻率、振型以及動(dòng)力響應(yīng)在 MATALB 中的實(shí)現(xiàn) ,說明了 MATLAB 在結(jié)構(gòu)動(dòng)力學(xué)計(jì)算中的強(qiáng)大功能及其編程的便捷性 ,使科技人員真正地從繁雜的計(jì)算中解放出來。 2公式2.1結(jié)構(gòu)自振特性和特征值 結(jié)構(gòu)自振特性是指結(jié)構(gòu)的振動(dòng)頻率和振型，計(jì)算經(jīng)驗(yàn)指出，結(jié)構(gòu)的阻尼對(duì)結(jié)構(gòu)的頻率和振型的影響

7、很小，所以求頻率振型時(shí)可以不考慮阻尼的影響，此時(shí)系統(tǒng)的自由振動(dòng)方程式如式（1）所示，即 （1） 當(dāng)系統(tǒng)做自由振動(dòng)時(shí)，各質(zhì)點(diǎn)做簡諧振動(dòng)，各節(jié)點(diǎn)的位移可表示為： （2） 將（2）代入（1）式，并消去公因子得到 （3）因此求解式（1）就是尋找式（3）的2 值和非零向量，這種問題稱為廣義特征值問題，記=2，和分別稱為廣義特征值和特征向量。式（3）可寫成 （4）這是一個(gè)齊次的線性方程組，若要有的非零解，系數(shù)行列式必須等于零，即 （5）展開此式可得 如果彈性結(jié)構(gòu)的總剛度矩陣K和總質(zhì)量矩陣M的階數(shù)都是n，則上述行列式展開后為的n次代數(shù)方程式，由此可求出n個(gè)根，即n個(gè)廣義特征值i，i=1,2,.,n，從而求出

8、結(jié)構(gòu)的n個(gè)自振頻率。求得廣義特征值i后，就可利用式（4）算得對(duì)應(yīng)的廣義特征向量i，它代表n個(gè)質(zhì)點(diǎn)的振幅構(gòu)成的振型。2.2用MATLAB對(duì)建筑結(jié)構(gòu)自振頻率、振型的分析如圖所示四層剛架結(jié)構(gòu) ,各層質(zhì)量分別為 m1 = 1kg, m2 = 2kg, m3 =3kg;m4=4kg.各層的側(cè)移剛度分別為 k1 = 800N /m , k2 = 1600N /m , k3 =3200N /m ,k4=6400N /m.求剛架的固有頻率和振型. 用matlab語言編程：% four_layerclc;clear;% k0每段的剛度k0(1)=800;k0(2)=1600;k0(3)=3200;k0(4)=6

9、400;% m每段的質(zhì)量m0(1)=1;m0(2)=2;m0(3)=3;m0(4)=4;% 層數(shù)n=4;% 定義m為質(zhì)量矩陣，k為總剛度矩陣m=zeros(n,n);k=zeros(n,n);% 計(jì)算mfor i=1:n; m(i,i)=m0(i);end% 計(jì)算kk(n,n)=k0(n);for i=1:n-1; k(i,i)=k0(i+1)+k0(i);endfor i=1:n-1; k(i,i+1)=-k0(i+1); k(i+1,i)=-k0(i+1);endmn=mk; %mn=inv(m)*k;% 求特征值w2=eig(mn);% 求角頻率w=sqrt(w2);% 頻率f=w/(2

10、*pi);% 周期T=1/f;for i=1:n; L=k-w2(i)*m; L00=L(2:n,2:n); L01=L(2:n,1); X=-inv(L00)*L01; xa(:,i)=X;endx1=ones(1,n);x=x1,xax = 1.0000 1.0000 1.0000 1.0000 -1.6299 -0.5000 0.6616 1.4682 2.1565 -0.2500 -0.0622 1.6557 -1.0125 0.2500 -0.3850 1.71002.3結(jié)構(gòu)動(dòng)力響應(yīng)求結(jié)構(gòu)的動(dòng)力響應(yīng)，就要對(duì)公式（6）進(jìn)行解答，可以用數(shù)值積分的方法對(duì)方程直接求解，即按時(shí)間增量t逐步求解

11、運(yùn)動(dòng)微分方程，直至反應(yīng)終了，這一方法稱作逐步積分法。這里只討論線性結(jié)構(gòu)體系的問題，逐步積分法求解運(yùn)動(dòng)微分方程的基本思路是：（1）把連續(xù)的時(shí)間過程離散為t1，t2，.，tn有限個(gè)點(diǎn)，對(duì)于運(yùn)動(dòng)微分方程 （6）求出其的位移、速度和加速度在有限個(gè)時(shí)間離散點(diǎn)上的值。在每個(gè)時(shí)間間隔t內(nèi)，假定位移、速度和加速度符合某一簡單的關(guān)系，而t的選擇要求保證計(jì)算的穩(wěn)定性與精確性。 從這樣的基本思路出發(fā)，本文采用Newmark-法來求解結(jié)構(gòu)動(dòng)力響應(yīng)。Newmark-法的計(jì)算步驟歸納如下：（1）基本數(shù)據(jù)準(zhǔn)備和初始條件計(jì)算： 1）選擇時(shí)間長t、參數(shù)和，并計(jì)算積分常數(shù) 2）確定運(yùn)動(dòng)的初始值。（2）形成剛度矩陣K，質(zhì)量矩陣M和

12、阻尼矩陣C（3）形成等效剛度矩陣，即計(jì)算ti+1時(shí)刻的等效荷載（5）求解ti+1時(shí)刻的位移，即計(jì)算ti+1時(shí)刻的加速度和速度循環(huán)第（4）至（6）計(jì)算步驟，可以得到線彈性體系在任一時(shí)刻的動(dòng)力反應(yīng)。2.4結(jié)構(gòu)在正弦荷載卸載后的自振響應(yīng) Newmark-b法的基本原理Newmark-b法是一種逐步積分的方法，避免了任何疊加的應(yīng)用，能很好的適應(yīng)非線性的反應(yīng)分析。Newmark-b法假定： (1-1) (1-2)式中，b和g是按積分的精度和穩(wěn)定性要求進(jìn)行調(diào)整的參數(shù)。當(dāng)b=0.5，g=0.25時(shí)，為常平均加速度法，即假定從t到t+Dt時(shí)刻的速度不變，取為常數(shù)。研究表明，當(dāng)b0.5， g0.25(0.5+b

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

14、算t+Dt時(shí)刻的有效荷載：（2）求解t+Dt時(shí)刻的位移： （3）計(jì)算t+Dt時(shí)刻的速度和加速度：Newmark-b方法是一種無條件穩(wěn)定的隱式積分格式，時(shí)間步長Dt的大小不影響解的穩(wěn)定性，Dt的選擇主要根據(jù)解的精度確定。瑞利矩陣瑞利阻尼矩陣?yán)萌鹄仃嚨恼恍裕|(zhì)量矩陣和剛度矩陣的正交性我們可得到 （其中為廣義阻尼矩陣） （7）計(jì)算卸載后的位移響應(yīng) 結(jié)構(gòu)的運(yùn)動(dòng)方程表達(dá)式為 （8）設(shè)方程的解為 （9）將（9）代入（8）可得 （10）左乘可得到（11）由瑞利阻尼矩陣、質(zhì)量矩陣和剛度矩陣的正交性可以得到 （12）將代入（12），兩邊同時(shí)除以可得（13）（13）的解為 （14） 其中 （15） （16）

15、由簡諧正弦荷載作用完畢時(shí)刻t=2s的結(jié)構(gòu)位移及速度條件作為結(jié)構(gòu)的自振初始條件： （17）左乘 （18） （19）同理 （20）3 動(dòng)力響應(yīng)分析假設(shè)上圖的四層框架結(jié)構(gòu)在頂部受一個(gè)簡諧荷載的作用，力的作用時(shí)間=5s，計(jì)算響應(yīng)的時(shí)間為100s，分2000步完成。阻尼矩陣由Rayleigh阻尼構(gòu)造。用matlab語言編程：clc;clear;% 質(zhì)量矩陣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*k;f0=100;t1=5;n

16、t=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 (tt1) 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*

17、(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; d(:,i)=kefe; 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% 質(zhì)點(diǎn)1figure(1)subplot(3,1,1),plot(d(1,:);title(1質(zhì)點(diǎn)的位移響應(yīng))subplot(3,1,2),plot(v(1,:);title(1質(zhì)點(diǎn)的速度響應(yīng)

18、)subplot(3,1,3),plot(a(1,:);title(1質(zhì)點(diǎn)的加速度響應(yīng))% 質(zhì)點(diǎn)2figure(2)subplot(3,1,1),plot(d(2,:);title(2質(zhì)點(diǎn)的位移響應(yīng))subplot(3,1,2),plot(v(2,:);title(2質(zhì)點(diǎn)的速度響應(yīng))subplot(3,1,3),plot(a(2,:);title(2質(zhì)點(diǎn)的加速度響應(yīng))% 質(zhì)點(diǎn)3figure(3)subplot(3,1,1),plot(d(3,:);title(3質(zhì)點(diǎn)的位移響應(yīng))subplot(3,1,2),plot(v(3,:);title(3質(zhì)點(diǎn)的速度響應(yīng))subplot(3,1,3),pl

19、ot(a(3,:);title(3質(zhì)點(diǎn)的加速度響應(yīng))% 質(zhì)點(diǎn)4figure(4)subplot(3,1,1),plot(d(4,:);title(4質(zhì)點(diǎn)的位移響應(yīng))subplot(3,1,2),plot(v(4,:);title(4質(zhì)點(diǎn)的速度響應(yīng))subplot(3,1,3),plot(a(4,:);title(4質(zhì)點(diǎn)的加速度響應(yīng))% 4個(gè)質(zhì)點(diǎn)的位移響應(yīng)figure(5)plot(d(1,:),b);hold on;plot(d(2,:),r);hold on;plot(d(3,:),k);hold on;plot(d(4,:),g);title(各個(gè)質(zhì)點(diǎn)的位移響應(yīng))legend(質(zhì)點(diǎn)1,質(zhì)點(diǎn)2,質(zhì)點(diǎn)3,質(zhì)點(diǎn)4);% 4個(gè)質(zhì)點(diǎn)的速度響應(yīng)figure(6)plot(v(1,:),b);h