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

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

2、ark-法、瑞利阻尼矩 The four layers of frame structure dynamic response based on MATLAB and research Xu 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 sel

3、f-vibration frequency 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

4、 initial conditions in 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 damping 1 引言 在社會發(fā)展的今天， 很多科技人員都會遇到數(shù)值分析計算機應(yīng)用等問題， 一些傳統(tǒng)的高級程序語

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

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

7、自振特性和特征值 結(jié)構(gòu)自振特性是指結(jié)構(gòu)的振動頻率和振型， 計算經(jīng)驗指出， 結(jié)構(gòu)的阻尼對結(jié)構(gòu)的頻率和振型的影響 很小，所以求頻率振型時可以不考慮阻尼的影響，此時系統(tǒng)的自由振動方程式如式（1）所示，即 Ku M u 0 1） 當(dāng)系統(tǒng)做自由振動時，各質(zhì)點做簡諧振動，各節(jié)點的位移可表示為： u sin( t ) 2） 將（ 2）代入（ 1）式，并消去公因子得到 2 K 2M 0 3） 因此求解式（ 1）就是尋找式（ 3）的2 值和非零向量 ，這種問題稱為廣義特征值問題，記 = 2，和 分別稱為廣義特征值和特征向量。式（ 3）可寫成 K M 0 4） 這是一個齊次的線性方程組，若要有 K M 的非零解，

8、 0 系數(shù)行列式必須等于零，即 5） 展開此式可得 K11 K21 M 11 M 21 K12 K22 M12 M 22 K1n K2n M 1n M 2n Kn1 M n1 Kn2 M n2 Knn M nn 如果彈性結(jié)構(gòu)的總剛度矩陣 K和總質(zhì)量矩陣 M 的階數(shù)都是 代數(shù)方程式，由此可求出 n，則上述行列式展開后為的 n 次 n 個根，即 n 個廣義特征值 i， i=1,2,.,n，從而求出結(jié)構(gòu)的 n 個自振頻率 i ,(i 1,2,.,n) 。 求得廣義特征值 i后，就可利用式（ 4）算得對應(yīng)的廣義特征向量 i，它代表 n 個質(zhì)點的振幅 構(gòu)成的振型 2.2 用 MATLAB對建筑結(jié)構(gòu)自振頻

9、率、振型的分析 如圖所示四層剛架結(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_layer clc; clear; % k0每段的剛度 k0(1)=800; k0(2)=1600; k0(3)=3200; k0(4)=6400; % m每段的質(zhì)量 m0(1)=1; m0(2)=2; m0(3)=3; m0(4)=4; % 層數(shù) n=4;

10、% 定義 m 為質(zhì)量矩陣， k 為總剛度矩陣 m=zeros(n,n); k=zeros(n,n); % 計算 m for i=1:n; m(i,i)=m0(i); end % 計算 k k(n,n)=k0(n); for i=1:n-1; k(i,i)=k0(i+1)+k0(i); end for i=1:n-1; k(i,i+1)=-k0(i+1); k(i+1,i)=-k0(i+1); end mn=mk; %mn=inv(m)*k; % 求特征值 w2=eig(mn); % 求角頻率 w=sqrt(w2); % 頻率 f=w/(2*pi); % 周期 T=1/f; for i=1:n;

11、 L=k-w2(i)*m; L00=L(2:n,2:n); L01=L(2:n,1); X=-inv(L00)*L01; xa(:,i)=X; end x1=ones(1,n); x=x1,xa 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.7100 2.3 結(jié)構(gòu)動力響應(yīng) 求結(jié)構(gòu)的動力響應(yīng)，就要對公式（ 6）進(jìn)行解答，可以用數(shù)值積分的方法對方程直接求解，即按 時間增量 t 逐步求解運動微分方程，直至反應(yīng)終了，這一方法

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

13、t,a3 1,a4 1,a5 2t (2), a6 t(1 ),a7t 2)確定運動的初始值 u 0、u 0和u 0。 (2) 形成剛度矩陣 K，質(zhì)量矩陣 M和阻尼矩陣 C (3) 形成等效剛度矩陣 K? ，即 K? K a0M a1C 計算 ti+1 時刻的等效荷載 P?i 1 Pi 1 Ma0ui a2ui a3ui Ca1ui a4ui a5ui ( 5)求解 ti+1 時刻的位移，即 K?ui1 P?i 1 計算 ti+1 時刻的加速度和速度 u i 1 a0(ui 1 ui ) a2ui a3u i u i 1 u i 1 a6u i a7 u i 1 循環(huán)第( 4)至( 6)計算步

14、驟，可以得到線彈性體系在任一時刻的動力反應(yīng)。 2.4 結(jié)構(gòu)在正弦荷載卸載后的自振響應(yīng) Newmark- 法的基本原理 Newmark- 法是一種逐步積分的方法，避免了任何疊加的應(yīng)用，能很好的適應(yīng)非線性的反應(yīng)分 析。 Newmark- 法假定： u t t ut (1 )u t ut t t (1-1) 12 u t t ut u t ( )ut u t t t2 2 (1-2) 式中， 和 是按積分的精度和穩(wěn)定性要求進(jìn)行調(diào)整的參數(shù)。當(dāng)=0.5， =0.25 時，為常平均 1 (u t u t t) 加速度法，即假定從 t 到 t+ t時刻的速度不變，取為常數(shù) 2 t t t 。研究表明，當(dāng) 0

15、.5， 0.25(0.5+ )2 時， Newmark- 法是一種無條件穩(wěn)定的格式。 由式(2-141)和式(2-142)可得到用 ut t及ut，ut，ut表示的ut t，ut t表達(dá)式，即有 u t t 1 1 1 2 (u t t ut) ut ( t2t 2 1) u t u t t t (u t t u t) (1)u t (1 2 ) tu 考慮 t+ t 時刻的振動微分方程為： M ut t Cut t Kut t R t t 將式(2-143)、式 (2-144) 代入(2-145)，得到關(guān)于 ut+ t 的方程 Kut t Rt t 式中 1 K K 2M C t 2t 11

16、1 R Rt t M( t2 uttut (2 1)u t ) C( ut ( 1)ut ( 1) tut ) t2 (1-3) (1-4) (1-5) (1-6) 求解式(2-146)可得 u t t ，然后由式 (2-143)和式(2-144)可解出ut 由此， Newmark- 法的計算步驟如下： 1.初始計算： (1)形成剛度矩陣 K、質(zhì)量矩陣 M和阻尼矩陣 C； t 和u t t 2) 3) 給定初始值 u0 ， u0和u0； 選擇積分步長 t、參數(shù) 、 ，并計算積分常數(shù) 1 t2 1 t， 1 32 2t (2) t(1 ) ， 7 t； K (4)形成有效剛度矩陣 K 2.對每個

17、時間步的計算： t 時刻的有效荷載： 0M 1C； 1) 計算 t+ Ft 2) t F t C( 求解 t+ t M ( 0ut 1ut t 時刻的位移： 4ut 2ut3ut ) 5ut) K ut t Ft t 3)計算 t+ t 時刻的速度和加速度： u t t 0(ut t u t)2ut3u t t 的大小不影響解的穩(wěn)定性， u t t ut6u t7u t t Newmark- 方法是一種無條件穩(wěn)定的隱式積分格式，時間步長 t 的選擇主要根據(jù)解的精度確定。 瑞利矩陣 瑞利阻尼矩陣 C a0M a1K 利用瑞利矩陣的正交性，質(zhì)量矩陣和剛度矩陣的正交性我們可得到 T T T X i

18、c X i a0 X i m X i a1 X i k X i ci a0 M ia1 K i 其中 ci 為廣義阻尼矩陣） ci2 i i M i i 2 i (a0 a1 i2) (7) 計算卸載后的位移響應(yīng) 結(jié)構(gòu)的運動方程表達(dá)式為 y(t) c y(t) y(t) 0 8） 設(shè)方程的解為 y(t) n 9） X i D i (t) i1 將（9） 代入（ 8）可得 n m( i 左乘 X T j 可得到 X iDi (t) 1 c( i1 X i Di (t) k( i1 X iDi (t) 0 （10） n T X Tj i1 T X Tj c X i Di i1 T X Tj k X

19、 i Di0 1 11) 由瑞利阻尼矩陣、質(zhì)量矩陣和剛度矩陣的正交性可以得到 M*jDj （t） C*jDj（t） K*jDj （t） 0 ( j 1,2 N)(12) 將C*j 2 j jM j 代入（ 12），兩邊同時除以 Mj 可得 Dj(t) 2 j jDj (t) jDj(t) 0 (13) 13) 的解為 t j) (14) Dj (t) Aje j jt sin( Djt 其中 tan j Aj Dj2(0) ( j j Dj(0) Dj(0)2 Dj 15) Dj (0) Dj Dj (0)j jDj(0) 16) 由簡諧正弦荷載作用完畢時刻 t=2s 的結(jié)構(gòu)位移及速度條件作為

20、結(jié)構(gòu)的自振初始條件： 左乘 X Tj m 同理 y(2) X i1 i Di (0) （17） X Tj m y(2) M*jDj (0) （18） Dj(0) X Tj m y(2) M *j （19） Dj(0) X Tj m y(2) * M *j （20） n 3 動力響應(yīng)分析 F =F0 sin( 4t t) 假設(shè)上圖的四層框架結(jié)構(gòu)在頂部受一個簡諧荷載 0 t1 的作用，力的作用時間 =5s，計算響 應(yīng)的時間為 100s，分 2000 步完成。阻尼矩陣由 Rayleigh阻尼構(gòu)造。 用 matlab 語言編程： clc; clear; % 質(zhì)量矩陣 m=1,2,3,4; m=diag

21、(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/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

22、=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*(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)=

23、v(:,i-1)+a6*a(:,i-1)+a7*a(:,i); end % 質(zhì)點 1 figure(1) subplot(3,1,1),plot(d(1,:);title(1 質(zhì)點的位移響應(yīng) ) subplot(3,1,2),plot(v(1,:);title(1 質(zhì)點的速度響應(yīng) ) subplot(3,1,3),plot(a(1,:);title(1 質(zhì)點的加速度響應(yīng) ) % 質(zhì)點 2 figure(2) subplot(3,1,1),plot(d(2,:);title(2 質(zhì)點的位移響應(yīng) ) subplot(3,1,2),plot(v(2,:);title(2 質(zhì)點的速度響應(yīng) ) subp

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