MATLAB在结构动力学中的应用.doc

MATLAB在结构动力学中的应用第17卷第1期2007年3月湖南工程学院JournalofHunanInstituteofEngineeringVo1.17.No.1M8r.20o7MATt.&B在结构动力学中的应用周后志,冷辉平(湖南科技大学土木工程学院,湖南湘潭411201)摘要:通过两个实例介绍了MATLAB语言在结构动力学中的应用,通过结构的自振频率,振型以及动力响应在MATALB中的实现,说明了MATLAB在结构动力学计算中的强大功能及其编程的便捷性,使科技人员真正地从繁杂的计算中解放出来.关键词:MATLAB;自振特性;动力响应中图分类号:TU311.3文献标识码:A文章编号:1671119X(2007)010091040引言在科学技术飞速发展的当今社会,计算机的应用能力已成为评价科技人员综合能力的一项重要内容.在进行科学研究与工程技术应用的过程中,科技人员往往会遇到大量繁重的数学运算和数值分析,一些传统的高级程序语言如FORTRAN等虽然能在一定程度上减轻计算量,但它们要求应用人员要具有较强的编程能力和对算法有深入的研究.另外,在运用这些高级程序语言进行计算结果的可视化分析及图形处理方面,对非计算机专业的普通用户来说,存在着很大的难度.MATLAB正是在这一应用要求背景下产生的数学类科技应用软件.MATLAB是Matrix和Laboratory前三个字母的缩写,是以矩阵计算为基础的程序设计语言,其指令格式与教科书中的数学表达式非常相近,用MAT-LAB编写程序犹如在便笺上列写公式和求解,因而被称为便笺式的编程语言uJ.同时,MATLAB具有功能丰富和完备的数学函数库及工具箱,大量繁杂的数学运算和分析可通过调用MATLAB函数直接求解,大大提高了编程效率,其程序编译和执行速度远远超过了传统的FORTRAN语言,因而用MAT-LAB编写程序,往往可以达到事半功倍的效果.在图形处理方面,MATLAB可以给数据以二维,三维乃至四维的直观表现,并在图形色彩,视角,品性等方面具有较强的渲染和控制能力,使科技人员对大量原始数据的分析变得轻松和得心应手J,从根本上满足了科技人员对工程数学计算的要求,将科技人员及普通用户从繁重的数学运算中解放出来.下面我们将演示运用MATLAB程序语言求解结构动力学中结构的自振频域,振型以及结构动力响应的例子,为此我们先回顾一下有关结构的自振特性以及动力响应的知识.1结构的自振特性和特征值J结构的自振特性是指结构的振动频率和振型,求结构的自振频率和振型也称对结构进行模态分析,是结构动力学计算的主要内容之一.计算经验指出,结构的阻尼对结构的频率和振型的影响很小,所以求频率振型时可以不考虑阻尼的影响.此时系统的自由振动方程如式(1)所示,即K+菇=O(1)当系统做自由振动时,各质点做简谐振动,各节点的位移可表示为=coswt(2)将(2)代入(1)式,并消去公凶子COSWt得到(K一)=O或K币=0)2币(3)因此求解式(1)就是寻找式(3)的tO值和非零向量,这种问题称为广义特征值问题.记A=,A和分别称为广义特征值和特征向量.式(3)可写成(K一)=O(4)收稿日期:2oo607一o6作者简介:周后志(1975一),男,硕士研究生,研究方向:工程结构损伤诊断与识别.湖南工程学院2007年这是一个齐次的线性方程组,若要有西的非零解,系数行列式必须等于零,即llAlO(5)展开此式可得.一A.-一A.一A.一A2一AK也一M也一A一AK一M=0如果弹性结构的总体刚度矩阵和总体质量矩阵的阶数都是T/,则上述行列式展开后为A的/7,次代数方程式,由此可求出T/,个根,即T/,个广义特征值A,=1,2,T/,从而求出结构的T/,个自振频率=,(=1,2,).求得广义特征值A后,就可利用式(4)算得对应的广义特征向量,它代表T/,个质点的振幅构成的振型.在弹性结构的振动分析中,结构的自由度总数T/,往往很大,因此无法直接从上述代数方程中求解广义特征值A(i=1,2,),而大多采用变换法和迭代法,如熟知的雅可比法,正迭代和逆迭代法等,采用这些方法计算的过程比较繁琐,而应用IVIATIdB程序语言则可使的求解结构的特征值的数学计算简单的多.2用MATLAB对结构自振频率,振型的分析用MATLAB对结构的自振频率,振型分析计算的便捷性用实例来证明:例题1如图1所示三层刚架结构,各层的楼面质量分别为ml=180,m2=270,m3=360;各层的侧移刚度分别为kl=98lVIN/m,k2=196lN/m,k3=294lVIN/m,求刚架的固有频率和振型.mfl80tKl=98MNm2=270tK2=l96MNm3=360tk3=294MN图1三层刚架运用matl语言编程,很容易就可得到:%求结构的自振特性(自振频率和振型)clear;%刚度,质量输入,编号从下至上为1,2,3ko(1)=98000000;ko(2)=196000000;ko(3)=2940o00o0:m0(1)=180000;m0(2)=270000;m0(3)=3600oo:n=3;%层数fori=1:n:ka(i,1)=k0(i);%组成刚度数组;m(i,i)=n(i);%组成质量矩阵;end%形成总刚度矩阵k(n,n)=ka(n);fori=1:n一1:k(i,i)=ka(i+1)+ka(i);endfori=1:n一1:k(i,i+1)=-k(i+1);k(i+1,i)=一ka(i+1);endmn=inv(m)=lck;w2=eig(nm);%求特征值W=sqrt(w2)%角频率f=w/(2=lc3.1415926)%频率T=1./f%周期%求振型fori=1:n:L=kw2(i)m;L00=L(2:n,2:n);LO1=L(2:n,1);X=一inv(too)=lcL01:xa(:,i)=X;endxl=ones(1,n);x=xl,)(a】运行程序得:x=1.00001.00001.0Oo0O.6485一O.6O662.54190.30180.67902.4396此程序可以求得任意多层刚架结构的固有频率和振型.3结构动力响应求结构的动力响应,在数学上就是要求出运动方程式(6)的解答.式(6)是一个二阶常系数微分方第1期周后志等:MATLAB在结构动力学中的应用93程组,可以用数值积分的方法对方程直接求解,即按时间增量t逐步求解运动微分方程,直至反应终了,这一方法称作逐步积分法.逐步积分法既可用于求解线性结构体系问题在整个动力反应过程中KM,C矩阵保持不变的问题;也可用于求解非线性结构体系的问题一KM,C矩阵随动力反应的过程而变化的问题.这里只讨论线性结构体系的问题.逐步积分法求解运动微分方程的基本思路是:(1)把连续的时间过程离散为t,t,t有限个点,对于运动微分方程M+c+Kx-F(6)只要求它们在上述每个时间离散点上得到满足,也就是说,最终求解得到的只是位移,速度和加速度在有限个时间离散点上的值,而不是连续函数x,.(2)在每个时间间隔At内,假定位移,速度和加速度符合某一简单的关系.而At的选择,要求保证计算的稳定性与精确性.从这样的基本思路出发,这里采用Wilson一0法来求解结构的动力响应.Wilson一0法的计算步骤归纳如下:(1)初始计算1)形成总体刚度矩阵和总体质量矩阵,而阻尼的影响一般不考虑.2)形成初始值戈(0),(0),并计算(0).3)选取时间步长At和(一般0=1.4),计算积分常数,得631OAt口2,A?a0口3丁,a2.3AtAt口5一6一7丁口84)形成有效刚度矩阵芷=K+口.+c(2)对于每一时间步长,循环计算1)计算HAt时刻的有效载荷向量(+OAt)=F(t)+0(F(+OAt)一F()+(a0()+a()+2()+C(aX()+2()+a3()2)求解HOAt时刻的位移向量(HOAt)佥(+OAt):分(+OAt)3)计算HAt时刻的加速度,速度,位移向量(+)=a(X(+OAt)一x(t)+a()+a()(+)=()+a,(1(+)+互()X(t+)=X(t)+()+口8(+)+2()4用MATLAB对结构动力响应的分析例题2有一简单的结构系统,不考虑阻尼影响,求此结构在03.5S的位移响应.已知结构初始状态是静止的,受突加不变荷载作用.(见图2)m1=2K1=2F=Im2=lk2=正图2二层刚架运用MATLAB语言编程得(其动力响应图见图3):结构动力位移响应时间图3结构动力响应%求结构的动力响应(wilson一0计算结构动力响应)theta=1.4:%l1eta=0delmt=O.28;%delta=t/n<O.1T(T:结构系统湖南工程学院2007年自由振动最小周期,由上面的程序可得T=2.8s)%计算初始加速度值Y(O),静止:v(O)=0,x(O)=O:运动方程my+cv+kx=P(t)y:加速度,v:速度,x:位移,P(t)=F;xO=0;O;vO=0;O;yO=0;10;aO=6/(thetadeltat);al=3/(thetadeltat);a2=2al;a3=thetadeltat/2;a4=aO/theta;a5=-a2/them;a6=1-3/theta;a7=deltat/2;a8=deltat2/6;K=k+a0m;%有效刚度矩阵fori=1:14F(1:2,i)=0;10;Y(1:2,i)=0;O;V(1:2,i)=0;O;x(1:2,i)=0;O;t=(i一1)deltat;fl=F(1:2,i)+theta0;0+m(a0】c0+a2v0+2),o)+C(a1.】c0+2v0+a3),o);%n=f(t+thetadeltat)xl=inv(K)fl;%xl=x(t+thetadeltat)y2=a4(xlx0)+a5,ro+a6y0;%y2=Y(t+deltat)v2=v0+a7(y2+),o);%v2=v(t+deltat)x2=x0+deltatv0+a8(y2+2),o);%x2=x(t+deltat)Y(1,i)=y2(1,1);Y(2,i)=y2(2,1);V(1,i)=v2(1,1);V(2,i)=v2(2,1);X(1,i)=x2(1,1);X(2,i)=x2(2,1);F(1,i)=0;F(2,i)=10;yO=y2;xO=x2;vO=v2;endt=0:deltat:3.64;X1=X(1,:)X2=X(2,:)plot(t,X1,窘:t,X2,的;axis(051515);title(结构动力位移响应);xlabel(时间);ylabel(幅度);运行程序得:X1=Columns1through8O.0o6o0.05250.19600.48960.95161.54252.16232.6702Columns9through142.92262.81822.33401.5415X2=Columns1through80.36631.33932.63943.92354.87935.30935.17814.6064Columns9through143.81823.06052.52332.28625结束语从以上两个程序及其运行图我们可以清楚地看到,运用MATLAB编程简单,大大提高了结构动力学中振动问题的求解效率,而且计算结果可用图像清楚直观地表达,效果良好.用MATLAB编写程序,其程序编译和执行速度远远超过了传统高级语言,可以达到事半功倍的效果.使科技人员对大量原始数据的分析变得轻松自如和得心应手,真正将科技人员从繁重的数学运算中解放了出来.参考文献1石博强.MATLAB数学计算范例教程M.北京:中国铁道出版社,2004.2云舟工作室.MATLAB教学建模基础教程M.北京:人民邮电出版社,2001.3(美)克夫,彭津.结构动力学M.王光远,等译.北京:科学出版社,1983.4侯新录结构分析中的有限元法与程序设计一用vj,alC实现M.j匕京:中国建捌业出J6I土,2m4.ApplicationofMATLABinStructureDynamicsZHOUHouzhi,LENGHuiping(SchoolofCivilEngineering,HunanScienceandTechnologyUniversity,Xiangtan411201,China)Abstract:ThispaperintroducesMATLABapplicationinthestructuredynamicsthroughtwoexamples,Withtherealizatio