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

1、.暨南大学研究生课程论文课程：结构动力学姓名：许可悦学号： 1634361002学院：力学与建筑工程学院专业：建筑与土木工程任课教师：李雪艳.基于 matlab 的四层框架结构动力响应与研究许可悦（暨南大学理工学院力学与土木工程学院，广州51063 ）摘要：本文用 matlab语言对四层建筑结构进行编程， 计算结构的自振频率、振型，分析该结构在自由振动和一般激励下的动力响应。采用了 newmark-法计算了在简谐正弦激励作用下结构的位移响应，并以此为初始条件结合瑞利阻尼矩阵计算了结构在简谐正弦荷载卸载后的结构自由振动的位移响应。关键词： matlab、newmark-法、瑞利阻尼矩the fo

2、ur layers of frame structure dynamic responsebased 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 编写程序, 往往可以达到事半功倍的效果 . 在图形处理方面,ma tlab 可以

6、给数据以二维、三维乃至四维的直观表现,并在图形色彩、视角、品性等方面具有较强的渲染和控制能力,使科技人员对大量原始数据的分析变得轻松和得心应手， 从根本上满足了科技人员对工程数学计算的要求,将科技人员及普通用户从繁重的数学运算中解放出来。.本文通过实例介绍了matlab 语言在结构动力学中的应用,通过结构的自振频率、振型以及动力响应在 matalb 中的实现 ,说明了 matlab 在结构动力学计算中的强大功能及其编程的便捷性 , 使科技人员真正地从繁杂的计算中解放出来。2 公式2.1 结构自振特性和特征值结构自振特性是指结构的振动频率和振型，计算经验指出，结构的阻尼对结构的频率和振型的影响很

7、小，所以求频率振型时可以不考虑阻尼的影响，此时系统的自由振动方程式如式（1）所示，即 k u m u 0（1）当系统做自由振动时，各质点做简谐振动，各节点的位移可表示为：u sin( t)（2）将（ 2）代入（ 1）式，并消去公因子得到 k 2 m 0（3）因此求解式（ 1）就是寻找式（ 3）的 2 值和非零向量 ，这种问题称为广义特征值问题，记=2，和 分别称为广义特征值和特征向量。式（3）可写成 k m 0（ 4）这是一个齐次的线性方程组，若要有的非零解，系数行列式必须等于零，即 k m 0（5）展开此式可得k11m 11k12m 12.kk 21m 21k 22m 22.k.k n1m

8、n1k n2m n2.k1nm2 nm.nn m1n2 nnn0如果弹性结构的总刚度矩阵k和总质量矩阵 m 的阶数都是 n，则上述行列式展开后为的n 次代数方程式，由此可求出n 个根，即 n 个广义特征值 i，从而求出结构的n个自振频率i=1,2,.,nii , (i 1,2,.,n) 。求得广义特征值 i后，就可利用式（ 4）算得对应的广义特征向量i ，它代表 n 个质点的振幅构成的振型。2.2 用 matlab对建筑结构自振频率、振型的分析如图所示四层刚架结构 ,各层质量分别为 m1 = 1kg, m2 = 2kg, m3 =3kg;m4=4kg各.层的侧移刚度分别为 k1 = 800n

9、/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);% 计算 mfor i=1:n;m(i,i)=m0(i);end% 计算 kk(n,n)=k0(n)

10、; 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+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.00001.00001.00001.0000-1.6

11、299-0.50000.66161.46822.1565-0.2500-0.06221.6557-1.01250.2500-0.38501.71002.3 结构动力响应求结构的动力响应，就要对公式（ 6）进行解答，可以用数值积分的方法对方程直接求解，即按时间增量 t 逐步求解运动微分方程，直至反应终了，这一方法称作逐步积分法。这里只讨论线性结构体系的问题，逐步积分法求解运动微分方程的基本思路是：（1）把连续的时间过程离散为t1，t2， .，tn 有限个点，对于运动微分方程 m u c u k u f （6）求出其的位移、速度和加速度在有限个时间离散点上的值。在每个时间间隔 t 内，假定位移、速

12、度和加速度符合某一简单的关系，而 t 的选择要求保证计算的稳定性与精确性。从这样的基本思路出发， 本文采用 newmark-法来求解结构动力响应。 newmark-法的计算步骤归纳如下：（1）基本数据准备和初始条件计算：1）选择时间长t 、参数和，并计算积分常数.a012 , a1, a21 , a311, a41, a5t (2),ttt22a6t (1 ), a7t2）确定运动的初始值 u 0、u 0 和 u 0 。（2）形成刚度矩阵 k，质量矩阵 m 和 阻尼矩阵 c?（3）形成等效刚度矩阵 k ，即? k k a0 m a1 c计算 ti+1 时刻的等效荷载? p i 1 m a0 u

13、 i a2 u i a3 u i c a1u i a4 u i a5 u i p i 1（5）求解 ti+1 时刻的位移，即? k u i 1 pi 1计算 ti+1 时刻的加速度和速度 u i 1a0 ( u i 1 u i ) a2 u i a3 u i u i 1 u i 1 a6 u i a7 u i 1循环第（ 4）至（ 6）计算步骤，可以得到线弹性体系在任一时刻的动力反应。2.4 结构在正弦荷载卸载后的自振响应newmark-法的基本原理newmark- 法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非线性的反应分析。newmark-法假定： u utt u t(1)

14、u tu u t u1) uttt (2tt t(1-1)t u t t t 2(1-2)式中，和 是按积分的精度和稳定性要求进行调整的参数。当=0.5，=0.25 时，为常平均加速度法，即假定从 t 到 t+1 ( u tu t t )t 时刻的速度不变，取为常数 2。研究表明，当 0.5， 0.25(0.5+ )2 时， newmark- 法是一种无条件稳定的格式。由式 (2-141)和式 (2-142)可得到用 u tt 及 u t ， u t ， u t 表示的 u tt ， u tt 表达式，即有. u utt12 ( u ttu t )1 u t( 11) u ttt2(1-3)t

15、tt( u ttu t ) (1) u t(1) t u t2(1-4)考虑 t+t 时刻的振动微分方程为： m u ttc u tt k u tt r t t(1-5)将式 (2-143)、式 (2-144) 代入 (2-145)，得到关于 ut+t 的方程 k u tt r tt(1-6)式中 k k 1 m c t 2t111 r r tt m (t 2 u tt u t( 21) u t ) c(u t(1) u t(1) t u t )t2求解式 (2-146)可得 u tt ，然后由式 (2-143)和式 (2-144)可解出 u tt 和 u t t 。由此， newmark-法

16、的计算步骤如下：1.初始计算：（1）形成刚度矩阵 k、质量矩阵 m 和阻尼矩阵 c；（2）给定初始值 u 0 ， u 0 和 u 0 ；（3）选择积分步长t、参数、，并计算积分常数11021t ，2t ，t，415t (2)，2， 6（4）形成有效刚度矩阵 k k 2.对每个时间步的计算：（1）计算 t+t 时刻的有效荷载：1312，t(1) ， 7t ；0 m 1 c； f t t f tt m ( 0 u t2 u t3 u t ) c(1u t4 u t5 u t )（2）求解 t+t 时刻的位移：k u tt f tt（3）计算 t+t 时刻的速度和加速度：. u tt0 ( u t

17、t u t )2 u t3 u t u tt u t6 u t7 u t tnewmark-方法是一种无条件稳定的隐式积分格式，时间步长t 的大小不影响解的稳定性，t 的选择主要根据解的精度确定。瑞利矩阵瑞利阻尼矩阵 c a0 m a 1 k 利用瑞利矩阵的正交性，质量矩阵和刚度矩阵的正交性我们可得到tttx i c x ia0 x i m x ia1 x i k x ic *ia0 m i*a1 k i*（其中 ci* 为广义阻尼矩阵）ci*2 ii m i*i1(a0 a1 i2 )2i（7）计算卸载后的位移响应结构的运动方程表达式为m y(t )c y(t )k y(t)0（ 8）设方程

18、的解为ny (t)x i d i (t )i1（9）将（ 9）代入（ 8）可得nnnm (x i di (t) )c (x i di (t)k (x i d i (t )0i 1i1i 1（10）左乘 x tj 可得到nnnttt0x j m x i d ix j c x i dix j k x i dii 1i 1i1（11）由瑞利阻尼矩阵、质量矩阵和刚度矩阵的正交性可以得到m *j d j (t)c*j d j(t)k *j d j (t )0( j1,2, n )（12）将 c*j 2 jj m *j 代入（ 12），两边同时除以 m *j可得d j ( t) 2 jj d j (t)j

19、 d j (t)0 （13）（13）的解为dj(t )aj ej j tsin( dj tj )（14）.其中j 12ajd j2 (0) ( j j d j (0) d j (0)2djjdj（ 15）tandj (0)djjj d j (0)d j (0)j（16）由简谐正弦荷载作用完毕时刻t=2s 的结构位移及速度条件作为结构的自振初始条件：ny (2)x i di (0)i 1（17）左乘x tj mx tj m y( 2)m*j d j(0)（18）tmy(2)x jd j (0)m *j（19）同理x tj m y(2)d j (0)m*j（20）3动力响应分析4tf = f 0

20、sin(t1)假设上图的四层框架结构在顶部受一个简谐荷载的作用，力的作用时间 =5s，计算响应的时间为 100s，分 2000 步完成。阻尼矩阵由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/al

21、fa/dt;a2=1/alfa/dt;a3=1/2/alfa-1;a4=beta/alfa-1;a5=dt/2*(beta/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

22、(:,i-1);%d(:,i)=inv(ke)*fe;d(:,i)=kefe;a(:,i)=a0*(d(:,i)-d(:,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)subp

23、lot(3,1,1),plot(d(2,:);title(2 质点的位移响应 )subplot(3,1,2),plot(v(2,:);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(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,: