基于MATLAB的四层框架结构动力响应与研究

暨南大学研究生课程论文课程：结构动力学姓名：许可悦学号：1634361002学院：力学与建筑工程学院专业：建筑与土木工程任课教师：李雪艳基于MATLAB的四层框架结构动力响应与研究许可悦（暨南大学理工学院力学与土木工程学院，广州51063）摘要：本文用MATLAB语言对四层建筑结构进行编程，计算结构的自振频率、振型，分析该结构在自由振动和一般激励下的动力响应。采用了Newmark-β法计算了在简谐正弦激励作用下结构的位移响应，并以此为初始条件结合瑞利阻尼矩阵计算了结构在简谐正弦荷载卸载后的结构自由振动的位移响应。关键词：MATLAB、Newmark-β法、瑞利阻尼矩ThefourlayersofframestructuredynamicresponsebasedonMATLABandresearchXuKeyue(Jinanuniversityinstituteofmechanicsandcivilengineeringdepartment,Guangzhou)Abstract：ThispaperusesMATLABlanguagetoprogramthethefourlayersofframestructure,calculatestheself-vibrationfrequencyandvibrationmodeofthestructure,andanalyzesthedynamicresponseofthestructureunderfreevibrationandgeneralexcitation.AdoptedtheNewmark-betamethodtocalculatethedisplacementofthestructureundertheactionofaharmonicsineexcitationresponse,andtheinitialconditionsincombinationwiththeRayleighdampingmatrixtocalculatethestructureinthestructureofharmonicsineloadafterunloadingfreevibrationdisplacementresponse.Keywords：MATLAB;Newmark-βmethod;Rayleighorthogonaldamping1引言在社会发展的今天，很多科技人员都会遇到数值分析计算机应用等问题，一些传统的高级程序语言如FORTRAN等虽然能在一定程度上减轻计算量,但它们要求应用人员要具有较强的编程能力和对算法有深入的研究.另外,在运用这些高级程序语言进行计算结果的可视化分析及图形处理方面,对非计算机专业的普通用户来说,存在着很大的难度.MATLAB正是在这一应用要求背景下产生的数学类科技应用软件。MATLAB是是以矩阵计算为基础的程序设计语言,MATLAB具有功能丰富和完备的数学函数库及工具箱,大量繁杂的数学运算和分析可通过调用MATLAB函数直接求解,大大提高了编程效率,其程序编译和执行速度远远超过了传统的FORTRAN语言,因而用MAT2LAB编写程序,往往可以达到事半功倍的效果.在图形处理方面,MATLAB可以给数据以二维、三维乃至四维的直观表现,并在图形色彩、视角、品性等方面具有较强的渲染和控制能力,使科技人员对大量原始数据的分析变得轻松和得心应手，从根本上满足了科技人员对工程数学计算的要求,将科技人员及普通用户从繁重的数学运算中解放出来。本文通过实例介绍了MATLAB语言在结构动力学中的应用,通过结构的自振频率、振型以及动力响应在MATALB中的实现,说明了MATLAB在结构动力学计算中的强大功能及其编程的便捷性,使科技人员真正地从繁杂的计算中解放出来。2公式2.1结构自振特性和特征值结构自振特性是指结构的振动频率和振型，计算经验指出，结构的阻尼对结构的频率和振型的影响很小，所以求频率振型时可以不考虑阻尼的影响，此时系统的自由振动方程式如式（1）所示，即0}]{[}]{[uMuK（1）当系统做自由振动时，各质点做简谐振动，各节点的位移可表示为：)sin(}{}{tu（2）将（2）代入（1）式，并消去公因子得到0][][2MK（3）因此求解式（1）就是寻找式（3）的ω2值和非零向量{φ}，这种问题称为广义特征值问题，记λ=ω2，λ和{φ}分别称为广义特征值和特征向量。式（3）可写成0][][MK（4）这是一个齐次的线性方程组，若要有{φ}的非零解，系数行列式必须等于零，即0][][MK（5）展开此式可得0.....................221122222221211112121111nnnnnnnnnnnnMKMKMKMKMKMKMKMKMK如果弹性结构的总刚度矩阵[K]和总质量矩阵[M]的阶数都是n，则上述行列式展开后为λ的n次代数方程式，由此可求出n个根，即n个广义特征值λi，i=1,2,...,n，从而求出结构的n个自振频率),...,2,1(,niii。求得广义特征值λi后，就可利用式（4）算得对应的广义特征向量{φi}，它代表n个质点的振幅构成的振型。2.2用MATLAB对建筑结构自振频率、振型的分析如图所示四层刚架结构,各层质量分别为m1=1kg,m2=2kg,m3=3kg;m4=4kg.各层的侧移刚度分别为k1=800N/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);%%计算mfori=1:n;m(i,i)=m0(i);end%%计算kk(n,n)=k0(n);fori=1:n-1;k(i,i)=k0(i+1)+k0(i);endfori=1:n-1;k(i,i+1)=-k0(i+1);k(i+1,i)=-k0(i+1);endmn=m\k;%mn=inv(m)*k;%%求特征值w2=eig(mn);%%求角频率w=sqrt(w2);%%频率f=w/(2*pi);%%周期T=1/f;fori=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',xa']'x=1.00001.00001.00001.0000-1.6299-0.50000.66161.46822.1565-0.2500-0.06221.6557-1.01250.2500-0.38501.71002.3结构动力响应求结构的动力响应，就要对公式（6）进行解答，可以用数值积分的方法对方程直接求解，即按时间增量Δt逐步求解运动微分方程，直至反应终了，这一方法称作逐步积分法。这里只讨论线性结构体系的问题，逐步积分法求解运动微分方程的基本思路是：（1）把连续的时间过程离散为t1，t2，...，tn有限个点，对于运动微分方程][}]{[}]{[}]{[FuKuCuM（6）求出其的位移、速度和加速度在有限个时间离散点上的值。在每个时间间隔Δt内，假定位移、速度和加速度符合某一简单的关系，而Δt的选择要求保证计算的稳定性与精确性。从这样的基本思路出发，本文采用Newmark-β法来求解结构动力响应。Newmark-β法的计算步骤归纳如下：（1）基本数据准备和初始条件计算：1）选择时间长Δt、参数β和γ，并计算积分常数),2(2,1,121,1,,15432120taaatatatatata76),1(2）确定运动的初始值000}{}{}{uuu和、。（2）形成刚度矩阵[K]，质量矩阵[M]和]阻尼矩阵[C]（3）形成等效刚度矩阵]ˆ[K，即][][][]ˆ[10CaMaKK计算ti+1时刻的等效荷载]}{}{}{][[]}{}{}{][[}{}ˆ{54132011iiiiiiiiuauauaCuauauaMPP（5）求解ti+1时刻的位移，即11]ˆ[}]{ˆ[iiPuK计算ti+1时刻的加速度和速度iiiiiuauauuau}{}{)}{}({}{3210117611}{}{}{}{iiiiuauauu循环第（4）至（6）计算步骤，可以得到线弹性体系在任一时刻的动力反应。2.4结构在正弦荷载卸载后的自振响应Newmark-Newmark-析。Newmark-tuuuutttttt]}{}){1[(}{}{(1-1)2]}{}){21[(}{}{}{tuutuuutttttt(1-2)=0.5=0.25时，为常平均加速度法，即假定从t到t+t时刻的速度不变，取为常数)}{}({21tttuu0.5，0.25(0.5+)2时，Newmark-由式(2-141)和式(2-142)可得到用ttu}{及tu}{，tu}{，tu}{表示的ttu}{，ttu}{表达式，即有tttttttuutuutu}){121(}{1)}{}({1}{2(1-3)tttttttutuuutu}{)21(}){1()}{}({}{(1-4)考虑t+t时刻的振动微分方程为：ttttttttRuKuCuM}{}]{[}]{[}]{[(1-5)将式(2-143)、式(2-144)代入(2-145)，得到关于ut+t的方程ttttRuK}{}]{[(1-6)式中][][1][][2CtMtKK)}{)12(}){1(}{]([)}){121(}{1}{1]([}{}{2ttttttttutuutCuututMRR求解式(2-146)可得ttu}{，然后由式(2-143)和式(2-144)可解出ttu}{和ttu}{。由此，Newmark-1.初始计算：（1）形成刚度矩阵[K]、质量矩阵[M]和阻尼矩阵[C]；（2）给定初始值0}{u，0}{u和0}{u；（3t201t，t1，t12，1213，14，)2(25t，)1(6t，t7；（4）形成有效刚度矩阵][][][][10CMKK；2.对每个时间步的计算：（1）计算t+t时刻的有效荷载：)}{}{}{]([)}{}{}{]([}{}{541320ttttttttttuuuCuuuMFF（2）求解t+t时刻的位移：ttttFuK}{}{（3）计算t+t时刻的速度和加速度：tttttttuuuuu}{}{)}{}({}{320ttttttuuuu}{}{}{}{76Newmark-t的大小不影响解的稳定性，t的选择主要根据解的精度确定。瑞利矩阵瑞利阻尼矩阵][][][10KaMaC利用瑞利矩阵的正交性，质量矩阵和刚度矩阵的正交性我们可得到iTiiTiiTiXkXaXmXaXcX10*1*0*iiiKaMac（