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

2、ur layers of frame structure dynamic responsebased 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

3、and vibration 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

4、 combination 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 引言在社會發(fā)展的今天， 很多科技人員都會遇到數值分析計算機應用等問題， 一些傳統(tǒng)的高級程序語言如 FORTRAN等雖然能在一定程度上減輕計算量 ,

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

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

7、小，所以求頻率振型時可以不考慮阻尼的影響，此時系統(tǒng)的自由振動方程式如式（1）所示，即 K u M u 0（1）當系統(tǒng)做自由振動時，各質點做簡諧振動，各節(jié)點的位移可表示為：u sin( t)（2）將（ 2）代入（ 1）式，并消去公因子得到 K 2 M 0（3）因此求解式（ 1）就是尋找式（ 3）的 2 值和非零向量 ，這種問題稱為廣義特征值問題，記=2，和 分別稱為廣義特征值和特征向量。式（3）可寫成 K M 0（ 4）這是一個齊次的線性方程組，若要有的非零解，系數行列式必須等于零，即 K M 0（5）展開此式可得K11M 11K12M 12.KK 21M 21K 22M 22.K.K n1M

8、n1K n2M n2.K1nM2 nM.nn M1n2 nnn0如果彈性結構的總剛度矩陣K和總質量矩陣 M 的階數都是 n，則上述行列式展開后為的n 次代數方程式，由此可求出n 個根，即 n 個廣義特征值 i，從而求出結構的n個自振頻率i=1,2,.,nii , (i 1,2,.,n) 。求得廣義特征值 i后，就可利用式（ 4）算得對應的廣義特征向量i ，它代表 n 個質點的振幅構成的振型。2.2 用 MATLAB對建筑結構自振頻率、振型的分析如圖所示四層剛架結構 ,各層質量分別為 m1 = 1kg, m2 = 2kg, m3 =3kg;m4=4kg各.層的側移剛度分別為 k1 = 800N

9、/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)=6400;% m每段的質量 m0(1)=1;m0(2)=2;m0(3)=3;m0(4)=4;% 層數n=4;% 定義 m 為質量矩陣， k 為總剛度矩陣 m=zeros(n,n);k=zeros(n,n);% 計算 mfor i=1:n;m(i,i)=m0(i);end% 計算 kk(n,n)=k0(n)

10、; 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*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.00001.00001.00001.0000-1.6

11、299-0.50000.66161.46822.1565-0.2500-0.06221.6557-1.01250.2500-0.38501.71002.3 結構動力響應求結構的動力響應，就要對公式（ 6）進行解答，可以用數值積分的方法對方程直接求解，即按時間增量 t 逐步求解運動微分方程，直至反應終了，這一方法稱作逐步積分法。這里只討論線性結構體系的問題，逐步積分法求解運動微分方程的基本思路是：（1）把連續(xù)的時間過程離散為t1，t2， .，tn 有限個點，對于運動微分方程 M u C u K u F （6）求出其的位移、速度和加速度在有限個時間離散點上的值。在每個時間間隔 t 內，假定位移、速

12、度和加速度符合某一簡單的關系，而 t 的選擇要求保證計算的穩(wěn)定性與精確性。從這樣的基本思路出發(fā)， 本文采用 Newmark-法來求解結構動力響應。 Newmark-法的計算步驟歸納如下：（1）基本數據準備和初始條件計算：1）選擇時間長t 、參數和，并計算積分常數.a012 , a1, a21 , a311, a41, a5t (2),ttt22a6t (1 ), a7t2）確定運動的初始值 u 0、u 0 和 u 0 。（2）形成剛度矩陣 K，質量矩陣 M 和 阻尼矩陣 C?（3）形成等效剛度矩陣 K ，即? K K a0 M a1 C計算 ti+1 時刻的等效荷載? P i 1 M a0 u

13、 i a2 u i a3 u i C a1u i a4 u i a5 u i P i 1（5）求解 ti+1 時刻的位移，即? K u i 1 Pi 1計算 ti+1 時刻的加速度和速度 u i 1a0 ( u i 1 u i ) a2 u i a3 u i u i 1 u i 1 a6 u i a7 u i 1循環(huán)第（ 4）至（ 6）計算步驟，可以得到線彈性體系在任一時刻的動力反應。2.4 結構在正弦荷載卸載后的自振響應Newmark-法的基本原理Newmark- 法是一種逐步積分的方法，避免了任何疊加的應用，能很好的適應非線性的反應分析。Newmark-法假定： u utt u t(1)

14、u tu u t u1) uttt (2tt t(1-1)t u t t t 2(1-2)式中，和 是按積分的精度和穩(wěn)定性要求進行調整的參數。當=0.5，=0.25 時，為常平均加速度法，即假定從 t 到 t+1 ( u tu t t )t 時刻的速度不變，取為常數 2。研究表明，當 0.5， 0.25(0.5+ )2 時， Newmark- 法是一種無條件穩(wěn)定的格式。由式 (2-141)和式 (2-142)可得到用 u tt 及 u t ， u t ， u t 表示的 u tt ， u tt 表達式，即有. u utt12 ( u ttu t )1 u t( 11) u ttt2(1-3)t

15、tt( u ttu t ) (1) u t(1) t u t2(1-4)考慮 t+t 時刻的振動微分方程為： M u ttC u tt K u tt R t t(1-5)將式 (2-143)、式 (2-144) 代入 (2-145)，得到關于 ut+t 的方程 K u tt R tt(1-6)式中 K K 1 M C t 2t111 R R tt M (t 2 u tt u t( 21) u t ) C(u t(1) u t(1) t u t )t2求解式 (2-146)可得 u tt ，然后由式 (2-143)和式 (2-144)可解出 u tt 和 u t t 。由此， Newmark-法

16、的計算步驟如下：1.初始計算：（1）形成剛度矩陣 K、質量矩陣 M 和阻尼矩陣 C；（2）給定初始值 u 0 ， u 0 和 u 0 ；（3）選擇積分步長t、參數、，并計算積分常數11021t ，2t ，t，415t (2)，2， 6（4）形成有效剛度矩陣 K K 2.對每個時間步的計算：（1）計算 t+t 時刻的有效荷載：1312，t(1) ， 7t ；0 M 1 C； F t t F tt M ( 0 u t2 u t3 u t ) C(1u t4 u t5 u t )（2）求解 t+t 時刻的位移：K u tt F tt（3）計算 t+t 時刻的速度和加速度：. u tt0 ( u t

17、t u t )2 u t3 u t u tt u t6 u t7 u t tNewmark-方法是一種無條件穩(wěn)定的隱式積分格式，時間步長t 的大小不影響解的穩(wěn)定性，t 的選擇主要根據解的精度確定。瑞利矩陣瑞利阻尼矩陣 C a0 M a 1 K 利用瑞利矩陣的正交性，質量矩陣和剛度矩陣的正交性我們可得到TTTX i c X ia0 X i m X ia1 X i k X ic *ia0 M i*a1 K i*（其中 ci* 為廣義阻尼矩陣）ci*2 ii M i*i1(a0 a1 i2 )2i（7）計算卸載后的位移響應結構的運動方程表達式為m y(t )c y(t )k y(t)0（ 8）設方程

18、的解為ny (t)X i D i (t )i1（9）將（ 9）代入（ 8）可得nnnm (X i Di (t) )c (X i Di (t)k (X i D i (t )0i 1i1i 1（10）左乘 X Tj 可得到nnnTTT0X j m X i D iX j c X i DiX j k X i Dii 1i 1i1（11）由瑞利阻尼矩陣、質量矩陣和剛度矩陣的正交性可以得到M *j D j (t)C*j D j(t)K *j D j (t )0( j1,2, N )（12）將 C*j 2 jj M *j 代入（ 12），兩邊同時除以 M *j可得D j ( t) 2 jj D j (t)j

19、 D j (t)0 （13）（13）的解為Dj(t )Aj ej j tsin( Dj tj )（14）.其中j 12AjD j2 (0) ( j j D j (0) D j (0)2DjjDj（ 15）tanDj (0)Djjj D j (0)D j (0)j（16）由簡諧正弦荷載作用完畢時刻t=2s 的結構位移及速度條件作為結構的自振初始條件：ny (2)X i Di (0)i 1（17）左乘X Tj mX Tj m y( 2)M*j D j(0)（18）Tmy(2)X jD j (0)M *j（19）同理X Tj m y(2)D j (0)M*j（20）3動力響應分析4tF = F 0

20、sin(t1)假設上圖的四層框架結構在頂部受一個簡諧荷載的作用，力的作用時間 =5s，計算響應的時間為 100s，分 2000 步完成。阻尼矩陣由Rayleigh 阻尼構造。用 matlab 語言編程：clc;clear;% 質量矩陣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;nt=2000;dt=0.01;.alfa=0.25;beta=0.5;a0=1/alfa/dt/dt;a1=beta/al

21、fa/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:ntt=(i-1)*dt;if (tt1)f=f0*sin(4*pi*t/t1);0;0;0;elsef=0;0;0;0;endke=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

22、(:,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% 質點 1figure(1)subplot(3,1,1),plot(d(1,:);title(1 質點的位移響應 ) subplot(3,1,2),plot(v(1,:);title(1 質點的速度響應 ) subplot(3,1,3),plot(a(1,:);title(1 質點的加速度響應 )% 質點 2figure(2)subp

23、lot(3,1,1),plot(d(2,:);title(2 質點的位移響應 )subplot(3,1,2),plot(v(2,:);title(2 質點的速度響應 )subplot(3,1,3),plot(a(2,:);title(2 質點的加速度響應 )% 質點 3figure(3).subplot(3,1,1),plot(d(3,:);title(3 質點的位移響應 )subplot(3,1,2),plot(v(3,:);title(3 質點的速度響應 )subplot(3,1,3),plot(a(3,:);title(3 質點的加速度響應 )% 質點 4figure(4)subplot(3,1,1),plot(d(4,:);title(4 質點的位移響應 ) subplot(3,1,2),plot(v(4,:);title(4 質點的速度響應 ) subplot(3,1,3),plot(a(4,:);title(4 質點的加速度響應 )% 4個質點的位移響應figure(5)plot(d(1,:),b);hold on;plot(d(2,:),r);hold on;plot(d(3,: