MATLAB在结构动力学中的应用

第17卷第1期湖南工程学院学报2007年3月Journal of Hunan Institute of EngineeriMATLAB在结构动力学中的应用周后志,冷辉平(湖南科技大学土木工程学院,湖南湘潭411201)摘要:通过两个实例介绍了 MATLAB语言在结构动力学中的应用,通过结构的自振频率、振型以及动力响应在 MATALB中的实现,说明了 MATLAB在结构动力学计算中的强大功能及其编程的便捷性,使科技人员真正地从繁杂的计算中解放出来关键词: MATLAB;自振特性;动力响应中图分类号:TU311.3文献标识码:A文章编号:1671-119X(2007)01-0091-04原始数据的分析变得轻松和得心应手(),从根本上0引言满足了科技人员对工程数学计算的要求,将科技人员及普通用户从繁重的数学运算中解放出来在科学技术飞速发展的当今社会,计算机的应下面我们将演示运用 MATLAB程序语言求解用能力已成为评价科技人员综合能力的一项重要内结构动力学中结构的自振频域、振型以及结构动力容在进行科学研究与工程技术应用的过程中,科技响应的例子,为此我们先回顾一下有关结构的自振人员往往会遇到大量繁重的数学运算和数值分析,特性以及动力响应的知识一些传统的高级程序语言如 FORTRAN等虽然能在1结构的自振特性和特征值定程度上减轻计算量,但它们要求应用人员要具有较强的编程能力和对算法有深人的研究另外,在结构的自振特性是指结构的振动频率和振型,运用这些髙级程序语言进行计算结果的可视化分析求结构的自振频率和振型也称对结构进行模态分及图形处理方面,对非计算机专业的普通用户来说,析,是结构动力学计算的主要内容之一计算经验指存在着很大的难度 MATLAB正是在这一应用要求出,结构的阻尼对结构的频率和振型的影响很小,所背景下产生的数学类科技应用软件以求频率振型时可以不考虑阻尼的影响此时系统MATLAB是 Matrix和 Laboratoty 1前三个字母的的自由振动方程如式(1)所示,即缩写,是以矩阵计算为基础的程序设计语言,其指令[K]{x}+[M]{x}=0格式与教科书中的数学表达式非常相近,用MAT当系统做自由振动时,各质点做简谐振动,各节LAB编写程序犹如在便笺上列写公式和求解,因而点的位移可表示为被称为“便笔式”的编程语言.同时, MATLAB具Ixi=lIe(2)有功能丰富和完备的数学函数库及工具箱,大量繁将(2)代入(1)式,并消去公因子 cost得到杂的数学运算和分析可通过调用 MATLAB函数直([K]-a2[M])|φ=0接求解,大大提高了编程效率,其程序编译和执行速或[K]{d}=a2M]{d}(3)度远远超过了传统的 FORTRAN语言,因而用MAT因此求解式(1)就是寻找式(3)的a2值和非零LAB编写程序,往往可以达到事半功倍的效果在向量{d},这种问题称为广义特征值问题记A=图形处理方面 MATLAB可以给数据以二维、三维乃a2,A和{d}分别称为广义特征值和特征向量至四维的直观表现,并在图形色彩、视角、品性等方面具有较强的渲染和控制能力,使科技人员对大量中国煤化工(4)CNMHG收稿日期:2006-07-06作者简介:周后志(1975-),男,硕士研究生,研究方向:工程结构损伤诊断与识别湖南工程学院学报2007年这是一个齐次的线性方程组,若要有{}的非运用maab语言编程,很容易就可得到零解,系数行列式必须等于零,即%求结构的自振特性(自振频率和振型)I[K]-A[M]=0clear展开此式可得%刚度、质量输入,编号从下至上为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;如果弹性结构的总体刚度矩阵[K和总体质量矩阵n=3;%层数M]的阶数都是n,则上述行列式展开后为A的nfor i=l: n;次代数方程式,由此可求出n个根,即n个广义特征ka(i,1)=k0(i);%组成刚度数组;值A;,=1,2,…,n,从而求出结构的n个自振频率m(i,i)=mo(i);%组成质量矩阵;,(i=1,2,…n)%形成总刚度矩阵求得广义特征值λ2后,就可利用式(4)算得对k(n, n)=ka(n)应的广义特征向量{中},它代表n个质点的振幅构for i=lan-1:成的振型在弹性结构的振动分析中,结构的自由度总数k(i, i)=ka(i+1)+ka(i)n往往很大,因此无法直接从上述代数方程中求解广义特征值λ(i=1,2,…,n),而大多采用变换法和迭代法,如熟知的雅可比法、正迭代和逆迭代法k(i,i+1)=-ka(i+1);k(i+1)=-ka(i+1)等,采用这些方法计算的过程比较繁琐,而应用MATLAB程序语言则可使的求解结构的特征值的数mn=inv(m)*kw2=eig(mn);%求特征值w=sqt(w2)%角学计算简单的多频率2用 MATLAB对结构自振频率振型的分析f=w/(2*3.1415926)%频率T=1./%周期用 MATLAB对结构的自振频率、振型分析计算%求振型for i=l:n;的便捷性用实例来证明例题1如图1所示三层刚架结构,各层的楼面L=k-w2(i)*m;D00=L(2:n,2:n);质量分别为ml=180t,m2=270t,m3=360t;各层01=L(2:n,1)的侧移刚度分别为kl=98MN/m,k2=19MNmX=-inv(L00)* LOl;xa(:, i)=Xk3=294MN/m,求刚架的固有频率和振型XK1-98MN运行程序得:x=1.00001.00001.00000.6485-0.60662.5419K2=196MN0.30180.67902.4396此程序可以求得任意多层刚架结构的固有频率和振型k3=294MN中国煤化工CNMHG求结构的动力响应,在数学上就是要求出运动图1三层刚架方程式(6)的解答式(6)是一个二阶常系数微分方第1期周后志等 MATLAB在结构动力学中的应用程组,可以用数值积分的方法对方程直接求解,即按[]{x(x+△)}=P(t+6△)时间增量Mt逐步求解运动微分方程,直至反应终3)计算t+△t时刻的加速度速度、位移向量了,这一方法称作逐步积分法逐步积分法既可用于x(t+△t)}=a4({X(t+0△)}-{x(l)})+求解线性结构体系问题——在整个动力反应过程中a5{x(t)}+a6{x(t)}K],[M],[C]矩阵保持不变的问题;也可用于求{x(t+△t)}={x(t)}+an(I(t+△)}+{x解非线性结构体系的问题-[K],[M],[C]矩阵随动力反应的过程而变化的问题这里只讨论线性结()})构体系的问题.逐步积分法求解运动微分方程的基X(t+△t)}={X(t)}+△t{x(t)}+a3({x(t本思路是:△t)}+2{x(t)})(1)把连续的时间过程离散为t1,t2,…,t有限个点,对于运动微分方程4用 MATLAB对结构动力响应的分析[M]{x}+[C]{x}+[K]{x}=[F]只要求它们在上述每个时间离散点上得到满例题2有一简单的结构系统,不考虑阻尼影足,也就是说最终求解得到的只是位移速度和加响求此结构在0-358的位移响应已知结构初速度在有限个时间离散点上的值而不是连续函数始状态是静止的,受突加不变荷载作用(见图2)X}、{x}、{x}(2)在每个时间间隔△t内,假定位移、速度和KI-2加速度符合某一简单的关系而△t的选择,要求保证计算的稳定性与精确性从这样的基本思路出发,这里采用Wlon-0k2=4法来求解结构的动力响应. Wilson-6法的计算步骤归纳如下(1)初始计算图2二层刚架1)形成总体刚度矩阵和总体质量矩阵,而阻尼运用 MATLAB语言编程得(其动力响应图见的影响一般不考虑图3):2)形成初始值{x(0)},{x(0)},并计算{x结构动力位移响应(0)3)选取时间步长△t和θ(一般θ=14),计算积分常数,得6(6△t)2q≈424)形成有效刚度矩阵[k]=[K]+a0[M]+a1[C](2)对于每一时间步长,循环计算时间1)计算t+△t时刻的有效载荷向量图3结构动力响应](t+0△t)}={F(1)}+({F(t+0△)}-F()1)+[M](ax()|+a1x()}+21x力响中国煤化工∞0n-计算结构动(t)})+[C](a1{X(t)}+2{x(t)}+a3{x(t)})CNMH2)求解t+6△t时刻的位移向量|X(t+0△t)deltat=0.28;% delta=t/n<0.1T(T:结构系统湖南工程学院学报2007年自由振动最小周期,由上面的程序可得T=28s) ylabel(‘幅度”);%计算初始加速度值y(0),静止:v(0)=0,x运行程序得(0)=0:运动方程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结束语f(t+theta deltat)从以上两个程序及其运行图我们可以清楚地看xl =inv (K)'fl; %xl =x(t +theta deltat)到,运用 MATLAB编程简单,大大提高了结构动力y2=a4·(x1-xo)+a5°v0+a6y0;%学中振动问题的求解效率,而且计算结果可用图像y2=y(t+deltat)清楚直观地表达,效果良好.用 MATLAB编写程序v2=v0+a7(2 +y0): %v2 =v(t+deltat其程序编译和执行速度远远超过了传统高级语言,x2=x0+ deltat yO+a8·(y2+2°y0);可以达到事半功倍的效果使科技人员对大量原始%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,)=x2(2,1);[1]石博强 MATLAB数学计算范例教程[M].北京:中国F(1,i)=0;F(2,i)=10;铁道出版社,2004.y0=y2;x0=x2;v0=v2;[2]云舟工作室 MATLAB教学建模基础教程[M].北京:人民邮电出版社,2001t=0: deltat: 3. 64[3](美)克夫,彭津结构动力学[M]王光远,等译北京:科学出版社,1983X1=X(1,)X2=X(2,)plot(t,Ⅺ1,g,t,2,);[4]侯新录结构分析中的有限元法与程序设计一用vsulaxis([05-1515])C+·实现[M]北京:中国建材工业出版社,204tde(‘结构动力位移响应); xlabel(时间);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-中国煤化工ful function and its convenience in the structure dyhelp people be frefrom the complicated calculationCNMHGKey words: MATLAB; free-vibration characteristic; dynamic response