第17卷第1期湖南工程學(xué)院學(xué)報(bào)2007年3月Journal of Hunan Institute of EngineeriMATLAB在結(jié)構(gòu)動(dòng)力學(xué)中的應(yīng)用周后志,冷輝平(湖南科技大學(xué)土木工程學(xué)院,湖南湘潭411201)摘要:通過(guò)兩個(gè)實(shí)例介紹了 MATLAB語(yǔ)言在結(jié)構(gòu)動(dòng)力學(xué)中的應(yīng)用,通過(guò)結(jié)構(gòu)的自振頻率、振型以及動(dòng)力響應(yīng)在 MATALB中的實(shí)現(xiàn),說(shuō)明了 MATLAB在結(jié)構(gòu)動(dòng)力學(xué)計(jì)算中的強(qiáng)大功能及其編程的便捷性,使科技人員真正地從繁雜的計(jì)算中解放出來(lái)關(guān)鍵詞: MATLAB;自振特性;動(dòng)力響應(yīng)中圖分類號(hào):TU311.3文獻(xiàn)標(biāo)識(shí)碼:A文章編號(hào):1671-119X(2007)01-0091-04原始數(shù)據(jù)的分析變得輕松和得心應(yīng)手(),從根本上0引言滿足了科技人員對(duì)工程數(shù)學(xué)計(jì)算的要求,將科技人員及普通用戶從繁重的數(shù)學(xué)運(yùn)算中解放出來(lái)在科學(xué)技術(shù)飛速發(fā)展的當(dāng)今社會(huì),計(jì)算機(jī)的應(yīng)下面我們將演示運(yùn)用 MATLAB程序語(yǔ)言求解用能力已成為評(píng)價(jià)科技人員綜合能力的一項(xiàng)重要內(nèi)結(jié)構(gòu)動(dòng)力學(xué)中結(jié)構(gòu)的自振頻域、振型以及結(jié)構(gòu)動(dòng)力容在進(jìn)行科學(xué)研究與工程技術(shù)應(yīng)用的過(guò)程中,科技響應(yīng)的例子,為此我們先回顧一下有關(guān)結(jié)構(gòu)的自振人員往往會(huì)遇到大量繁重的數(shù)學(xué)運(yùn)算和數(shù)值分析,特性以及動(dòng)力響應(yīng)的知識(shí)一些傳統(tǒng)的高級(jí)程序語(yǔ)言如 FORTRAN等雖然能在1結(jié)構(gòu)的自振特性和特征值定程度上減輕計(jì)算量,但它們要求應(yīng)用人員要具有較強(qiáng)的編程能力和對(duì)算法有深人的研究另外,在結(jié)構(gòu)的自振特性是指結(jié)構(gòu)的振動(dòng)頻率和振型,運(yùn)用這些髙級(jí)程序語(yǔ)言進(jìn)行計(jì)算結(jié)果的可視化分析求結(jié)構(gòu)的自振頻率和振型也稱對(duì)結(jié)構(gòu)進(jìn)行模態(tài)分及圖形處理方面,對(duì)非計(jì)算機(jī)專業(yè)的普通用戶來(lái)說(shuō),析,是結(jié)構(gòu)動(dòng)力學(xué)計(jì)算的主要內(nèi)容之一計(jì)算經(jīng)驗(yàn)指存在著很大的難度 MATLAB正是在這一應(yīng)用要求出,結(jié)構(gòu)的阻尼對(duì)結(jié)構(gòu)的頻率和振型的影響很小,所背景下產(chǎn)生的數(shù)學(xué)類科技應(yīng)用軟件以求頻率振型時(shí)可以不考慮阻尼的影響此時(shí)系統(tǒng)MATLAB是 Matrix和 Laboratoty 1前三個(gè)字母的的自由振動(dòng)方程如式(1)所示,即縮寫,是以矩陣計(jì)算為基礎(chǔ)的程序設(shè)計(jì)語(yǔ)言,其指令[K]{x}+[M]{x}=0格式與教科書中的數(shù)學(xué)表達(dá)式非常相近,用MAT當(dāng)系統(tǒng)做自由振動(dòng)時(shí),各質(zhì)點(diǎn)做簡(jiǎn)諧振動(dòng),各節(jié)LAB編寫程序猶如在便箋上列寫公式和求解,因而點(diǎn)的位移可表示為被稱為“便筆式”的編程語(yǔ)言.同時(shí), MATLAB具Ixi=lIe(2)有功能豐富和完備的數(shù)學(xué)函數(shù)庫(kù)及工具箱,大量繁將(2)代入(1)式,并消去公因子 cost得到雜的數(shù)學(xué)運(yùn)算和分析可通過(guò)調(diào)用 MATLAB函數(shù)直([K]-a2[M])|φ=0接求解,大大提高了編程效率,其程序編譯和執(zhí)行速或[K]osqwwom0=a2M]osqwwom0(3)度遠(yuǎn)遠(yuǎn)超過(guò)了傳統(tǒng)的 FORTRAN語(yǔ)言,因而用MAT因此求解式(1)就是尋找式(3)的a2值和非零LAB編寫程序,往往可以達(dá)到事半功倍的效果在向量osqwwom0,這種問(wèn)題稱為廣義特征值問(wèn)題記A=圖形處理方面 MATLAB可以給數(shù)據(jù)以二維、三維乃a2,A和osqwwom0分別稱為廣義特征值和特征向量至四維的直觀表現(xiàn),并在圖形色彩、視角、品性等方面具有較強(qiáng)的渲染和控制能力,使科技人員對(duì)大量中國(guó)煤化工(4)CNMHG收稿日期:2006-07-06作者簡(jiǎn)介:周后志(1975-),男,碩士研究生,研究方向:工程結(jié)構(gòu)損傷診斷與識(shí)別湖南工程學(xué)院學(xué)報(bào)2007年這是一個(gè)齊次的線性方程組,若要有{}的非運(yùn)用maab語(yǔ)言編程,很容易就可得到零解,系數(shù)行列式必須等于零,即%求結(jié)構(gòu)的自振特性(自振頻率和振型)I[K]-A[M]=0clear展開此式可得%剛度、質(zhì)量輸入,編號(hào)從下至上為1、23Ku-AMu K1-AMr.. Kin-AM,k0(1)=980000002)=196000000(3)K2-AM2I Kx -AM2..K2, -AM2=0=29400000m0(1)=18000m0(2)=2700;m0(3)=Kn- AM,, K-AMa…Km-AM。360000;如果彈性結(jié)構(gòu)的總體剛度矩陣[K和總體質(zhì)量矩陣n=3;%層數(shù)M]的階數(shù)都是n,則上述行列式展開后為A的nfor i=l: n;次代數(shù)方程式,由此可求出n個(gè)根,即n個(gè)廣義特征ka(i,1)=k0(i);%組成剛度數(shù)組;值A(chǔ);,=1,2,…,n,從而求出結(jié)構(gòu)的n個(gè)自振頻率m(i,i)=mo(i);%組成質(zhì)量矩陣;,(i=1,2,…n)%形成總剛度矩陣求得廣義特征值λ2后,就可利用式(4)算得對(duì)k(n, n)=ka(n)應(yīng)的廣義特征向量{中},它代表n個(gè)質(zhì)點(diǎn)的振幅構(gòu)for i=lan-1:成的振型在彈性結(jié)構(gòu)的振動(dòng)分析中,結(jié)構(gòu)的自由度總數(shù)k(i, i)=ka(i+1)+ka(i)n往往很大,因此無(wú)法直接從上述代數(shù)方程中求解廣義特征值λ(i=1,2,…,n),而大多采用變換法和迭代法,如熟知的雅可比法、正迭代和逆迭代法k(i,i+1)=-ka(i+1);k(i+1)=-ka(i+1)等,采用這些方法計(jì)算的過(guò)程比較繁瑣,而應(yīng)用MATLAB程序語(yǔ)言則可使的求解結(jié)構(gòu)的特征值的數(shù)mn=inv(m)*kw2=eig(mn);%求特征值w=sqt(w2)%角學(xué)計(jì)算簡(jiǎn)單的多頻率2用 MATLAB對(duì)結(jié)構(gòu)自振頻率振型的分析f=w/(2*3.1415926)%頻率T=1./%周期用 MATLAB對(duì)結(jié)構(gòu)的自振頻率、振型分析計(jì)算%求振型for i=l:n;的便捷性用實(shí)例來(lái)證明例題1如圖1所示三層剛架結(jié)構(gòu),各層的樓面L=k-w2(i)*m;D00=L(2:n,2:n);質(zhì)量分別為ml=180t,m2=270t,m3=360t;各層01=L(2:n,1)的側(cè)移剛度分別為kl=98MN/m,k2=19MNmX=-inv(L00)* LOl;xa(:, i)=Xk3=294MN/m,求剛架的固有頻率和振型XK1-98MN運(yùn)行程序得:x=1.00001.00001.00000.6485-0.60662.5419K2=196MN0.30180.67902.4396此程序可以求得任意多層剛架結(jié)構(gòu)的固有頻率和振型k3=294MN中國(guó)煤化工CNMHG求結(jié)構(gòu)的動(dòng)力響應(yīng),在數(shù)學(xué)上就是要求出運(yùn)動(dòng)圖1三層剛架方程式(6)的解答式(6)是一個(gè)二階常系數(shù)微分方第1期周后志等 MATLAB在結(jié)構(gòu)動(dòng)力學(xué)中的應(yīng)用程組,可以用數(shù)值積分的方法對(duì)方程直接求解,即按[]{x(x+△)}=P(t+6△)時(shí)間增量Mt逐步求解運(yùn)動(dòng)微分方程,直至反應(yīng)終3)計(jì)算t+△t時(shí)刻的加速度速度、位移向量了,這一方法稱作逐步積分法逐步積分法既可用于x(t+△t)}=a4({X(t+0△)}-{x(l)})+求解線性結(jié)構(gòu)體系問(wèn)題——在整個(gè)動(dòng)力反應(yīng)過(guò)程中a5{x(t)}+a6{x(t)}K],[M],[C]矩陣保持不變的問(wèn)題;也可用于求{x(t+△t)}={x(t)}+an(I(t+△)}+{x解非線性結(jié)構(gòu)體系的問(wèn)題-[K],[M],[C]矩陣隨動(dòng)力反應(yīng)的過(guò)程而變化的問(wèn)題這里只討論線性結(jié)()})構(gòu)體系的問(wèn)題.逐步積分法求解運(yùn)動(dòng)微分方程的基X(t+△t)}={X(t)}+△t{x(t)}+a3({x(t本思路是:△t)}+2{x(t)})(1)把連續(xù)的時(shí)間過(guò)程離散為t1,t2,…,t有限個(gè)點(diǎn),對(duì)于運(yùn)動(dòng)微分方程4用 MATLAB對(duì)結(jié)構(gòu)動(dòng)力響應(yīng)的分析[M]{x}+[C]{x}+[K]{x}=[F]只要求它們?cè)谏鲜雒總€(gè)時(shí)間離散點(diǎn)上得到滿例題2有一簡(jiǎn)單的結(jié)構(gòu)系統(tǒng),不考慮阻尼影足,也就是說(shuō)最終求解得到的只是位移速度和加響求此結(jié)構(gòu)在0-358的位移響應(yīng)已知結(jié)構(gòu)初速度在有限個(gè)時(shí)間離散點(diǎn)上的值而不是連續(xù)函數(shù)始狀態(tài)是靜止的,受突加不變荷載作用(見圖2)X}、{x}、{x}(2)在每個(gè)時(shí)間間隔△t內(nèi),假定位移、速度和KI-2加速度符合某一簡(jiǎn)單的關(guān)系而△t的選擇,要求保證計(jì)算的穩(wěn)定性與精確性從這樣的基本思路出發(fā),這里采用Wlon-0k2=4法來(lái)求解結(jié)構(gòu)的動(dòng)力響應(yīng). Wilson-6法的計(jì)算步驟歸納如下(1)初始計(jì)算圖2二層剛架1)形成總體剛度矩陣和總體質(zhì)量矩陣,而阻尼運(yùn)用 MATLAB語(yǔ)言編程得(其動(dòng)力響應(yīng)圖見的影響一般不考慮圖3):2)形成初始值{x(0)},{x(0)},并計(jì)算{x結(jié)構(gòu)動(dòng)力位移響應(yīng)(0)3)選取時(shí)間步長(zhǎng)△t和θ(一般θ=14),計(jì)算積分常數(shù),得6(6△t)2q≈424)形成有效剛度矩陣[k]=[K]+a0[M]+a1[C](2)對(duì)于每一時(shí)間步長(zhǎng),循環(huán)計(jì)算時(shí)間1)計(jì)算t+△t時(shí)刻的有效載荷向量圖3結(jié)構(gòu)動(dòng)力響應(yīng)](t+0△t)}={F(1)}+({F(t+0△)}-F()1)+[M](ax()|+a1x()}+21x力響中國(guó)煤化工∞0n-計(jì)算結(jié)構(gòu)動(dòng)(t)})+[C](a1{X(t)}+2{x(t)}+a3{x(t)})CNMH2)求解t+6△t時(shí)刻的位移向量|X(t+0△t)deltat=0.28;% delta=t/n<0.1T(T:結(jié)構(gòu)系統(tǒng)湖南工程學(xué)院學(xué)報(bào)2007年自由振動(dòng)最小周期,由上面的程序可得T=28s) ylabel(‘幅度”);%計(jì)算初始加速度值y(0),靜止:v(0)=0,x運(yùn)行程序得(0)=0:運(yùn)動(dòng)方程my+ew+kx=p(t)y:加速度,vX1=速度,x:位移,p(t)=F;xO=[0;0];v0=[0;0];Columns 1 through 8y0=[0;10];0.00600.05250.19600.48960.95161.5425a0=6/(theta'deltat)2; al =3/( theta 'deltat); 2. 1623 2. 67022 =2 al: a3 =theta'deltat/2houa4=a0theta; a5 =-a2/theta; a6=1-3/theta;2.92262.81822.33401.5415a7 =deltat/2; a8= deltat 2/6K=k+a0·m;%有效剛度矩陣Columns 1 through 8for i=1:140.36631.33932.63943.92354.87935.3093F(1:2,i)=[0;10];Y(1:2,i)=[0;0];v(1:5.17814.60642,i)=[0;0];X(1:2,i)=[0;0]Columns 9 through 14t=(i-1)· deltat3.81823.06052.52332.2862f1=F(1:2,i)+ theta[0;0]+m(a0x0+a20+2°y0)+c'(alx0+20+a3y0);%f5結(jié)束語(yǔ)f(t+theta deltat)從以上兩個(gè)程序及其運(yùn)行圖我們可以清楚地看xl =inv (K)'fl; %xl =x(t +theta deltat)到,運(yùn)用 MATLAB編程簡(jiǎn)單,大大提高了結(jié)構(gòu)動(dòng)力y2=a4·(x1-xo)+a5°v0+a6y0;%學(xué)中振動(dòng)問(wèn)題的求解效率,而且計(jì)算結(jié)果可用圖像y2=y(t+deltat)清楚直觀地表達(dá),效果良好.用 MATLAB編寫程序v2=v0+a7(2 +y0): %v2 =v(t+deltat其程序編譯和執(zhí)行速度遠(yuǎn)遠(yuǎn)超過(guò)了傳統(tǒng)高級(jí)語(yǔ)言,x2=x0+ deltat yO+a8·(y2+2°y0);可以達(dá)到事半功倍的效果使科技人員對(duì)大量原始%x2=x(t+ deltat)數(shù)據(jù)的分析變得輕松自如和得心應(yīng)手,真正將科技Y(1,i)=y2(1,1);Y(2,i)=y2(2,1);人員從繁重的數(shù)學(xué)運(yùn)算中解放了出來(lái)v(1,i)=v2(1,1);V(2,i)=v2(2,1);參考文獻(xiàn)X(1,i)=x2(1,1);X(2,)=x2(2,1);[1]石博強(qiáng) MATLAB數(shù)學(xué)計(jì)算范例教程[M].北京:中國(guó)F(1,i)=0;F(2,i)=10;鐵道出版社,2004.y0=y2;x0=x2;v0=v2;[2]云舟工作室 MATLAB教學(xué)建?；A(chǔ)教程[M].北京:人民郵電出版社,2001t=0: deltat: 3. 64[3](美)克夫,彭津結(jié)構(gòu)動(dòng)力學(xué)[M]王光遠(yuǎn),等譯北京:科學(xué)出版社,1983X1=X(1,)X2=X(2,)plot(t,Ⅺ1,g,t,2,);[4]侯新錄結(jié)構(gòu)分析中的有限元法與程序設(shè)計(jì)一用vsulaxis([05-1515])C+·實(shí)現(xiàn)[M]北京:中國(guó)建材工業(yè)出版社,204tde(‘結(jié)構(gòu)動(dòng)力位移響應(yīng)); xlabel(時(shí)間);pplication of MATLAB in Structure DynamicsZHOU Hou-zhi, LENG Hui -pingSchool of Civil Engineering, Hunan Science and Technology University, Xiangtan 411201, China)Abstract: This paper introduces MATLAB application in the structure dynamics through two examples, With therealization of the natural vibration frequency, vibration mode andtik rB, matLab power-中國(guó)煤化工ful function and its convenience in the structure dyhelp people be frefrom the complicated calculationCNMHGKey words: MATLAB; free-vibration characteristic; dynamic response