1、暨南大学研究生课程论文课程：结构动力学姓名：许可悦学号：1634361002学院：力学与建筑工程学院专业：建筑与土木工程任课教师：李雪艳基于 MATLAB 的四层框架结构动力响应与研究许可悦（暨南大学理工学院力学与土木工程学院，广州 51063）摘要：本文用 MATLAB 语言对 四层建筑结构进行编程，计算结构的自振频率、振型，分析该结构在自由振动和一般激励下的动力响应。采用了 Newmark- 法计算了在简谐正弦激励作用下结构的位移响应，并以此为初始条件结合瑞利阻尼矩阵计算了结构在简谐正弦荷载卸载后的结构自由振动的位移响应。关键词：MATLAB、Newmark- 法、瑞利阻尼矩The fou
2、r layers of frame structure dynamic response based on MATLAB and researchXu 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 self-vibration frequency
3、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 initial conditions in
4、 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 damping1 引言在社会发展的今天，很多科技人员都会遇到数值分析计算机应用等问题，一些传统的高级程序语言如 FORTRAN 等虽然能在一定程度上减轻计算量
5、,但它们要求应用人员要具有较强的编程能力和对算法有深入的研究. 另外 ,在运用这些高级程序语言进行计算结果的可视化分析及图形处理方面 ,对非计算机专业的普通用户来说 ,存在着很大的难度. MATLAB 正是在这一应用要求背景下产生的数学类科技应用软件。MATLAB 是是以矩阵计算为基础的程序设计语言 , MATLAB 具有功能丰富和完备的数学函数库及工具箱 ,大量繁杂的数学运算和分析可通过调用 MATLAB 函数直接求解 ,大大提高了编程效率 ,其程序编译和执行速度远远超过了传统的 FORTRAN 语言 ,因而用 MA T2LAB 编写程序 , 往往可以达到事半功倍的效果. 在图形处理方面 ,
6、MA TLAB 可以给数据以二维、三维乃至四维的直观表现 ,并在图形色彩、视角、品性等方面具有较强的渲染和控制能力 ,使科技人员对大量原始数据的分析变得轻松和得心应手，从根本上满足了科技人员对工程数学计算的要求 ,将科技人员及普通用户从繁重的数学运算中解放出来。本文通过实例介绍了 MATLAB 语言在结构动力学中的应用 ,通过结构的自振频率、振型以及动力响应在 MATALB 中的实现 ,说明了 MATLAB 在结构动力学计算中的强大功能及其编程的便捷性 ,使科技人员真正地从繁杂的计算中解放出来。2 公式2.1 结构自振特性和特征值结构自振特性是指结构的振动频率和振型，计算经验指出，结构的阻尼对
7、结构的频率和振型的影响很小，所以求频率振型时可以不考虑阻尼的影响，此时系统的自由振动方程式如式（1）所示，即（1）0uMK当系统做自由振动时，各质点做简谐振动，各节点的位移可表示为：（2）)sin(tu将（2）代入（1）式，并消去公因子得到（3）02MK因此求解式（1）就是寻找式（3）的 2 值和非零向量，这种问题称为广义特征值问题，记= 2， 和 分别称为广义特征值和特征向量。式（3）可写成（4）0MK这是一个齐次的线性方程组，若要有 的非零解，系数行列式必须等于零，即（5）0展开此式可得0.21 222 1121 nnn nMKMK如果弹性结构的总刚度矩阵K和总质量矩阵M的阶数都是 n，则
8、上述行列式展开后为 的 n次代数方程式，由此可求出 n 个根，即 n 个广义特征值 i，i=1,2,.,n，从而求出结构的 n 个自振频率 。),.21(,ii求得广义特征值 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.求刚架的固有频率和振型.用 mat
9、lab 语言编程：% 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);% 计算 mfor i=1:n;m(i,i)=m0(i);end% 计算 kk(n,n)=k0(n);for i=1:n-1;k(i,i)=k0(i+1)+k0(i);endfor i=1:n-1;k(i,i+1)=-k0(i+1);k(i
10、+1,i)=-k0(i+1);endmn=mk; %mn=inv(m)*k;% 求特征值w2=eig(mn);% 求角频率w=sqrt(w2);% 频率f=w/(2*pi);% 周期T=1/f;for i=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,xax =1.0000 1.0000 1.0000 1.0000-1.6299 -0.5000 0.6616 1.46822.1565 -0.2500 -0.0622 1.6557-1.0125 0.2500 -
11、0.3850 1.71002.3 结构动力响应求结构的动力响应，就要对公式（6）进行解答，可以用数值积分的方法对方程直接求解，即按时间增量 t 逐步求解运动微分方程，直至反应终了，这一方法称作逐步积分法。这里只讨论线性结构体系的问题，逐步积分法求解运动微分方程的基本思路是：（1）把连续的时间过程离散为 t1，t2 ，.，tn 有限个点，对于运动微分方程（6） FuKCuM求出其的位移、速度和加速度在有限个时间离散点上的值。在每个时间间隔 t 内，假定位移、速度和加速度符合某一简单的关系，而 t 的选择要求保证计算的稳定性与精确性。从这样的基本思路出发，本文采用 Newmark- 法来求解结构动
12、力响应。Newmark- 法的计算步骤归纳如下：（1）基本数据准备和初始条件计算：1）选择时间长 t、参数 和 ，并计算积分常数),2(,1,21,1 543220 taattattt76),(2）确定运动的初始值 。00u和、（2）形成刚度矩阵K，质量矩阵M和 阻尼矩阵C（3）形成等效刚度矩阵 ，即K10CaMK计算 ti+1 时刻的等效荷载 54132011 iiiiiiii uauaCuaP （5）求解 ti+1 时刻的位移，即11iiuK计算 ti+1 时刻的加速度和速度 iiiii uauau)(32101 17611 iiii 循环第（4）至（6）计算步骤，可以得到线弹性体系在任一
13、时刻的动力反应。2.4 结构在正弦荷载卸载后的自振响应Newmark- 法的基本原理Newmark- 法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非线性的反应分析。Newmark- 法假定：(1-1)tuutttt )1(1-2)2)2(tt tttt 式中， 和 是按积分的精度和稳定性要求进行调整的参数。当 =0.5， =0.25 时，为常平均加速度法，即假定从 t 到 t+ t 时刻的速度不变，取为常数 。研究表明，当)(21ttu 0.5， 0.25(0.5+ )2 时，Newmark- 法是一种无条件稳定的格式。由式(2-141)和式(2-142)可得到用 及 ， ， 表
14、示的 ， 表达式，即有tuttut tut(1-3)tttttuu )12(1)(12 (1-4)ttttt uu)()1()( 考虑 t+ t 时刻的振动微分方程为：(1-5)tttt RuKCuM 将式(2-143)、式(2-144) 代入(2-145)，得到关于 ut+ t 的方程(1-6)ttRuK式中 12CtMt )12()1(2 ttt tttt uuuCR求解式(2-146)可得 ，然后由式(2-143) 和式(2-144)可解出 和 。t tut由此，Newmark- 法的计算步骤如下：1.初始计算：（1）形成刚度矩阵K、质量矩阵M和阻尼矩阵C；（2）给定初始值 ， 和 ；0
15、u0u（3）选择积分步长 t、参数 、 ，并计算积分常数， ， ， ，201tt1t12123， ， ， ；4)(5t )1(6tt7（4）形成有效刚度矩阵 ；10CMK2.对每个时间步的计算：（1）计算 t+ t 时刻的有效荷载： )( )(541 320ttt ttt uuCF（2）求解 t+ t 时刻的位移： ttFuK（3）计算 t+ t 时刻的速度和加速度： ttttt uuu )( 320 tttt 76Newmark- 方法是一种无条件稳定的隐式积分格式，时间步长 t 的大小不影响解的稳定性， t 的选择主要根据解的精度确定。瑞利矩阵瑞利阻尼矩阵 10KaMC利用瑞利矩阵的正交性
16、，质量矩阵和刚度矩阵的正交性我们可得到iTiiTiiTi XkamXac10（其中 为广义阻尼矩阵）*10*iii KaMc*ic*2iiiMc（7）)(21210iiia计算卸载后的位移响应结构的运动方程表达式为（8）0)()()(tyktctym设方程的解为（9）niitDXty1)()(将（9）代入（8）可得（10）0)()()( 111 niiniinii tXktXctDXm左乘 可得到TjX（11）0111 ni iTjniiTjni iTj DXkXcDm由瑞利阻尼矩阵、质量矩阵和刚度矩阵的正交性可以得到（12）0)()()(*tKtCtMjjj ),21(Nj将 代入（12），
17、两边同时除以 可得*2jjjC*jM（13）0)()()(tDttDjjjj（13）的解为（14）)sin()(jDjtjjeAtj 其中（15）21jjDj 22)0()0DjjjjA（16）)0()(tanjjDj由简谐正弦荷载作用完毕时刻 t=2s 的结构位移及速度条件作为结构的自振初始条件：（17）niiDXy1)0()2(左乘 mXTj（18）)0()2(*jTjMymX（19）*)0(jTjD同理（20）*)2()0(jTjMymX3 动力响应分析假设上图的四层框架结构在顶部受一个简谐荷载 的作用，力的作用时间=5s ，计算响014=sin()tF应的时间为 100s，分 2000
18、 步完成。阻尼矩阵由 Rayleigh 阻尼构造。用 matlab 语言编程：clc;clear;% 质量矩阵m=1,2,3,4;m=diag(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/
19、alfa-2);a6=dt*(1-beta);a7=dt*beta;d=zeros(4,nt);v=zeros(4,nt);a=zeros(4,nt);for i=2:ntt=(i-1)*dt;if (tt1)f=f0*sin(4*pi*t/t1);0;0;0; elsef=0;0;0;0; endke=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(:
20、,i-1)-a2*v(:,i-1)-a3*a(:,i-1);v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);end% 质点 1figure(1)subplot(3,1,1),plot(d(1,:);title(1 质点的位移响应)subplot(3,1,2),plot(v(1,:);title(1 质点的速度响应)subplot(3,1,3),plot(a(1,:);title(1 质点的加速度响应 )% 质点 2figure(2)subplot(3,1,1),plot(d(2,:);title(2 质点的位移响应)subplot(3,1,2),plot(v(2,:
21、);title(2 质点的速度响应)subplot(3,1,3),plot(a(2,:);title(2 质点的加速度响应 )% 质点 3figure(3)subplot(3,1,1),plot(d(3,:);title(3 质点的位移响应)subplot(3,1,2),plot(v(3,:);title(3 质点的速度响应)subplot(3,1,3),plot(a(3,:);title(3 质点的加速度响应 )% 质点 4figure(4)subplot(3,1,1),plot(d(4,:);title(4 质点的位移响应)subplot(3,1,2),plot(v(4,:);title(
22、4 质点的速度响应)subplot(3,1,3),plot(a(4,:);title(4 质点的加速度响应 )% 4 个质点的位移响应figure(5)plot(d(1,:),b);hold on;plot(d(2,:),r);hold on;plot(d(3,:),k);hold on;plot(d(4,:),g);title(各个质点的位移响应)legend(质点 1,质点 2,质点 3,质点 4);% 4 个质点的速度响应figure(6)plot(v(1,:),b);hold on;plot(v(2,:),r);hold on;plot(v(3,:),k);hold on;plot(v
23、(4,:),g);title(各个质点的速度响应)legend(质点 1,质点 2,质点 3,质点 4);% 4 个质点的加速度响应figure(7)plot(a(1,:),b);hold on;plot(a(2,:),r);hold on;plot(a(3,:),k);hold on;plot(a(4,:),g);title(各个质点的加速度响应)legend(质点 1,质点 2,质点 3,质点 4);输出其位移、速度和加速度图像：4 结语MATLAB 是以矩阵计算为基础的程序设计软件，其指令格式与教科书中的数学表达式非常接近，大量繁杂的数学运算和分析可通过调用 MATLAB 函数直接求解，大大提高了编程的效率，可以达到事半功倍的效果。本文通过实例介绍了 MATLAB 语言在结构动力学中的应用 ,通过结构的自振频率、振型以及动力响应在 MATALB 中的实现 ,说明了 MATLAB 在结构动力学计算中的强大功能及其编程的便捷性 ,使科技人员真正地从繁杂的计算中解放出来 .参考文献1R.克拉夫.结构动力学M.高等教育出版社,2006,112刘晶波. 结构动力学M.机械工业出版社,2005,013陈立宇.基于 MATLAB 的高层框架结构动力响应与仿真研究.安徽建筑工业学院.20124周后志,冷辉平.MATLAB 在结构动力学中的应用.湖南科技大学土木工程学院.2007
