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

暨南大学研究生课程论文课程：结构动力学姓名：许可悦学号：1634361002学院：力学与建筑工程学院专业：建筑与土木工程任课教师：李雪艳基于MATLAB的四层框架结构动力响应与研究许可悦（暨南大学理工学院力学与土木工程学院，广州 51063）摘要：本文用MATLAB语言对四层建筑结构进行编程，计算结构的自振频率、振型，分析该结构在自由振动和一般激励下的动力响应。采用了Newmark-法计算了在简谐正弦激励作用下结构的位移响应，并以此为初始条件结合瑞利阻尼矩阵计算了结构在简谐正弦荷载卸载后的结构自由振动的位移响应。关键词：MATLAB、Newmark-法、瑞利阻尼矩The four 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 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 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 等虽然能在一定程度上减轻计算量 ,但它们要求应用人员要具有较强的编程能力和对算法有深入的研究. 另外 ,在运用这些高级程序语言进行计算结果的可视化分析及图形处理方面 ,对非计算机专业的普通用户来说 ,存在着很大的难度. MATLAB 正是在这一应用要求背景下产生的数学类科技应用软件。MATLAB 是是以矩阵计算为基础的程序设计语言 , MATLAB 具有功能丰富和完备的数学函数库及工具箱 ,大量繁杂的数学运算和分析可通过调用 MATLAB 函数直接求解 ,大大提高了编程效率 ,其程序编译和执行速度远远超过了传统的 FORTRAN 语言 ,因而用 MA T2LAB 编写程序 , 往往可以达到事半功倍的效果. 在图形处理方面 ,MA TLAB 可以给数据以二维、三维乃至四维的直观表现 ,并在图形色彩、视角、品性等方面具有较强的渲染和控制能力 ,使科技人员对大量原始数据的分析变得轻松和得心应手，从根本上满足了科技人员对工程数学计算的要求 ,将科技人员及普通用户从繁重的数学运算中解放出来。本文通过实例介绍了 MATLAB 语言在结构动力学中的应用 ,通过结构的自振频率、振型以及动力响应在 MATALB 中的实现 ,说明了 MATLAB 在结构动力学计算中的强大功能及其编程的便捷性 ,使科技人员真正地从繁杂的计算中解放出来。 2公式2.1结构自振特性和特征值 结构自振特性是指结构的振动频率和振型，计算经验指出，结构的阻尼对结构的频率和振型的影响很小，所以求频率振型时可以不考虑阻尼的影响，此时系统的自由振动方程式如式（1）所示，即 （1） 当系统做自由振动时，各质点做简谐振动，各节点的位移可表示为： （2） 将（2）代入（1）式，并消去公因子得到 （3）因此求解式（1）就是寻找式（3）的2 值和非零向量，这种问题称为广义特征值问题，记=2，和分别称为广义特征值和特征向量。式（3）可写成 （4）这是一个齐次的线性方程组，若要有的非零解，系数行列式必须等于零，即 （5）展开此式可得 如果弹性结构的总刚度矩阵K和总质量矩阵M的阶数都是n，则上述行列式展开后为的n次代数方程式，由此可求出n个根，即n个广义特征值i，i=1,2,.,n，从而求出结构的n个自振频率。求得广义特征值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);% 计算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+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.4682 2.1565 -0.2500 -0.0622 1.6557 -1.0125 0.2500 -0.3850 1.71002.3结构动力响应求结构的动力响应，就要对公式（6）进行解答，可以用数值积分的方法对方程直接求解，即按时间增量t逐步求解运动微分方程，直至反应终了，这一方法称作逐步积分法。这里只讨论线性结构体系的问题，逐步积分法求解运动微分方程的基本思路是：（1）把连续的时间过程离散为t1，t2，.，tn有限个点，对于运动微分方程 （6）求出其的位移、速度和加速度在有限个时间离散点上的值。在每个时间间隔t内，假定位移、速度和加速度符合某一简单的关系，而t的选择要求保证计算的稳定性与精确性。 从这样的基本思路出发，本文采用Newmark-法来求解结构动力响应。Newmark-法的计算步骤归纳如下：（1）基本数据准备和初始条件计算： 1）选择时间长t、参数和，并计算积分常数 2）确定运动的初始值。（2）形成刚度矩阵K，质量矩阵M和阻尼矩阵C（3）形成等效刚度矩阵，即计算ti+1时刻的等效荷载（5）求解ti+1时刻的位移，即计算ti+1时刻的加速度和速度循环第（4）至（6）计算步骤，可以得到线弹性体系在任一时刻的动力反应。2.4结构在正弦荷载卸载后的自振响应 Newmark-b法的基本原理Newmark-b法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非线性的反应分析。Newmark-b法假定： (1-1) (1-2)式中，b和g是按积分的精度和稳定性要求进行调整的参数。当b=0.5，g=0.25时，为常平均加速度法，即假定从t到t+Dt时刻的速度不变，取为常数。研究表明，当b0.5， g0.25(0.5+b)2时，Newmark-b法是一种无条件稳定的格式。由式(2-141)和式(2-142)可得到用及，表示的，表达式，即有 (1-3) (1-4)考虑t+Dt时刻的振动微分方程为： (1-5)将式(2-143)、式(2-144) 代入(2-145)，得到关于ut+Dt的方程 (1-6)式中求解式(2-146)可得，然后由式(2-143)和式(2-144)可解出和。由此，Newmark-b法的计算步骤如下：1.初始计算：（1）形成刚度矩阵K、质量矩阵M和阻尼矩阵C；（2）给定初始值， 和；（3）选择积分步长Dt、参数b、g，并计算积分常数，；（4）形成有效刚度矩阵；2.对每个时间步的计算：（1）计算t+Dt时刻的有效荷载：（2）求解t+Dt时刻的位移： （3）计算t+Dt时刻的速度和加速度：Newmark-b方法是一种无条件稳定的隐式积分格式，时间步长Dt的大小不影响解的稳定性，Dt的选择主要根据解的精度确定。瑞利矩阵瑞利阻尼矩阵利用瑞利矩阵的正交性，质量矩阵和刚度矩阵的正交性我们可得到 （其中为广义阻尼矩阵） （7）计算卸载后的位移响应 结构的运动方程表达式为 （8）设方程的解为 （9）将（9）代入（8）可得 （10）左乘可得到（11）由瑞利阻尼矩阵、质量矩阵和刚度矩阵的正交性可以得到 （12）将代入（12），两边同时除以可得（13）（13）的解为 （14） 其中 （15） （16）由简谐正弦荷载作用完毕时刻t=2s的结构位移及速度条件作为结构的自振初始条件： （17）左乘 （18） （19）同理 （20）3 动力响应分析假设上图的四层框架结构在顶部受一个简谐荷载的作用，力的作用时间=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/alfa/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:nt t=(i-1)*dt; if (tt1) f=f0*sin(4*pi*t/t1);0;0;0; else f=0;0;0;0; end ke=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(:,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,:);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,:),k);hold on;plot(d(4,:),g);title(各个质点的位移响应)legend(质点1,质点2,质点3,质点4);% 4个质点的速度响应figure(6)plot(v(1,:),b);h