• ES1524
    了解作者
  • matlab
    开发工具
  • 9KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 2
    下载次数
  • 2020-04-06 18:37
    上传日期
能自动计算单自由度在简谐振动具有初始位移和初始速度的响应程序
jianxie.rar
  • jianxie
  • M638.m
    1.1KB
  • danjianxiezhendong.m
    773B
  • guangyiafa.m
    6.9KB
  • danshuxianggezhen.m
    459B
  • newmark_ode.m
    707B
  • gezhen.m
    7.2KB
  • pingmiangangjia.m
    6.6KB
内容介绍
%%%%———结构动力学多自由度体系计算,包括隔震支座,采用Newmark法为基础,————— %————考虑Bouc-Wen模型进行计算—————————————————— %(考虑了隔震支座刚度修改的情况) %%%——————————————————————————————————————— clear; %清除环境中存在的变量 clc; %清屏 tic; %计算程序运行时间 %—————————————初始量设定—————————————————————— %————质量矩阵—————————————— m=[4.5 3.456 3.456 3.456 3.456 3.456 3.456 3.456 3.456]*1e5; %各层的集中质量单位为kg M=diag(m); %质量矩阵 %————————刚度矩阵—————————— kkk=[1.805 3.404 3.257 2.849 2.686 2.43 2.073 1.687 1.366]*1e8; %各层层间刚度,单位为N/m k1=[kkk(2) kkk(3) kkk(4) kkk(5) kkk(6) kkk(7) kkk(8) kkk(9)]; %层间刚度,单位为N/m,为形成刚度矩阵用 k2=[kkk(1)+kkk(2),kkk(2)+kkk(3),kkk(3)+kkk(4),kkk(4)+kkk(5),kkk(5)+kkk(6),kkk(6)+kkk(7),kkk(7)+kkk(8),kkk(8)+kkk(9),kkk(9)]; k3=[kkk(2) kkk(3) kkk(4) kkk(5) kkk(6) kkk(7) kkk(8) kkk(9)]; K=diag(k2)-diag(k1,1)-diag(k3,-1); %初始刚度矩阵(如何形成初始刚度矩阵?公式怎么来的?) %——————阻尼矩阵—————————————— c=[26.17 490 467 410 386 348 298 243 196]*1e3; %各层阻尼系数,单位是N*sec/m c1=[c(2),c(3),c(4),c(5),c(6),c(7),c(8),c(9)]; %为形成阻尼矩阵所用 c2=[c(1)+c(2),c(2)+c(3),c(3)+c(4),c(4)+c(5),c(5)+c(6),c(6)+c(7),c(7)+c(8),c(8)+c(9),c(9)]; c3=[c(2),c(3),c(4),c(5),c(6),c(7),c(8),c(9)]; C=diag(c2)-diag(c1,1)-diag(c3,-1); %阻尼矩阵(如何形成阻尼刚度矩阵?公式怎么来的?) %————————————考虑隔震支座计算中涉及的参数取值———————————— mb=4.5e5; %隔震支座的质量,单位为 kg kb=1.805e8; %隔震支座的初始刚度,单位为 N/m cb=26.17e3; %隔震支座的阻尼,单位为 N*sec/m alphaB=0.6; DyB=0.04; %单位换算成m AB=1.0; betaB=0.5; gamaB=0.5; nB=3; %———————————Newmark法计算多自由度结构反应——————————————— %——————————————基本参数设定———————————————————— ag=load('E:\qiao\tqqdaima\Earthquake_1000.txt'); %导入加速度数据 acceleration0=ag(1:10000); %换算成m/s^2 (1:10000) dt=0.001; %计算时间步长,同地震波数据时间步长一致 t=0:dt:(length(acceleration0)-1)*dt; %计算时长 gama=0.5; %定义gama值 beta=0.25; %定义beta值 am=0.5; %定义am值 af=0.5; %定义af值 L=ones(9,1); %影响系数取值 p=-M*L*acceleration0'; %地震力作用 %———————Newmark法——————————————————————— %——————参数计算—————— Gaa=M/(beta*dt)+gama*C/beta; % Gbb=M/(2*beta)+dt*C*(gama/(2*beta)-1);%newmark方法里面的二个参数,等效荷载增量(加速度的量的参数),书本6.2.28 {振动力学基础与matlab应用} %————————初始量赋值—————————— displacement(:,1)=zeros(9,1); %初始相对位移向量 velocity(:,1)=zeros(9,1); %初始相对速度向量 acceleration(:,1)=zeros(9,1); %初始相对加速度向量 %——————隔震支座的初始量赋值—————————— RestoringForce(1)=0; % The stiffness restoring force of the lead-core rubber-bearing xb(1)=0; %位移 dxb(1)=0; %速度,位移对时间的导数 vb(1)=0; %跟时间有关的参数 dvb(1)=0; %vb对时间的导数 %——————做循环计算————————————————————————————— N=length(acceleration0)-1; for i=1:N %—————————————————————————————————————— Ka=K;%(1-af)*K; GKK=K+gama*C/(beta*dt)+M/(beta*dt^2);%newmark方法里面的参数,等效刚度矩阵,书本6.2.28 {振动力学基础与matlab应用} %—————————————————————————————————————— dp(:,i)=(p(:,i+1)-p(:,i))+Gaa*velocity(:,i)+Gbb*acceleration(:,i);%书本6.2.28 {振动力学基础与matlab应用}等效刚度矩阵 dq(:,i)=GKK\dp(:,i);%书本6.2.28 {振动力学基础与matlab应用} 6.2.27 dq1(:,i)=gama*dq(:,i)/(beta*dt)-gama*velocity(:,i)/beta+dt*acceleration(:,i)*(1-gama/(2*beta));%书本6.2.28 {振动力学基础与matlab应用}6.2.25 dq2(:,i)=dq(:,i)/(beta*dt^2)-velocity(:,i)/(beta*dt)-acceleration(:,i)/(2*beta);%书本6.2.28 {振动力学基础与matlab应用}6.2.26 displacement(:,i+1)=displacement(:,i)+dq(:,i); velocity(:,i+1)=velocity(:,i)+dq1(:,i); acceleration(:,i+1)=acceleration(:,i)+dq2(:,i); %%————————隔震支座刚度矩阵的修正计算————————————————— xb(i+1)=displacement(1,i+1); %隔震支座的相对位移(为什么?没看懂) dxb(i+1)=velocity(1,i+1); %隔震支座的相对速度 vb(i+1)=vb(i)+dvb(i)*dt;%这个表示什么?上面一个表示速度dxb,vb表示哪一个量? dvb(i+1)=(AB*dxb(i)-betaB*abs(dxb(i))*abs(vb(i))^(nB-1)*vb(i)-gamaB*dxb(i)*abs(vb(i)^nB))/DyB; dvbdvx=(AB-betaB*abs(dxb(i+1))/dxb(i+1)*(abs(vb(i+1)))^(nB-1)*vb(i+1)-gamaB*(abs(vb(i+1))^nB))/DyB; kbb=alphaB*kb+(1-alphaB)*kb*DyB*dvbdvx; K(1,1)=kkk(2)+kbb; %修正后的刚度,并进入下一次的循环计算中 %%—————————————————————————————————————— RestoringForce(i+1)=alphaB*kb*xb(i+1)+(1-alphaB)*kb*DyB*vb(i+1); end %——————————————计算各层的绝对加速度———————————————— for i=1:length(m) absacceleration(i,:)=acceleration(i,:)+acceleration0'; %各层绝对加速度 end %——————————————反应量计算结果绘图————————————————— %——————滞回力曲线—————————————— figure plot(xb,RestoringForce); xlabel('位移tb/m'),ylabel('滞回力RestoringForce/N'); title('Bouc-Wen模型隔震支座刚度滞回力曲线'); set(gca,'FontSize',20); set(get(gca,'XLabel'),'FontSize',20); set(get(gca,'YLabel'),'FontSize',20); set(get(gca,'title'),'FontSize',20); grid on %——————————隔震支座的位移和绝对加速度—————————————————— figure plot(t,xb); xlabel('时间t/s'),ylabel('相对位移xb/m'); title('隔震支座相对位移'); set(gca,'FontSize',20); set(get(gca,'XLabel'),'FontSize',20); set(get(gca,'YLabel'),'FontSize',20); set(get(gca,'title'),'FontSize',20); grid on figure plot(t,absacceleration(1,:)); xlabel('时间t/s'),ylabel('绝对加速度m/s^2'); title('隔震支座相对位移'); set(gca,'FontSize',20); set(get(gca,'XLabel'),'FontSize',20); set(get(gca,'YLabel'),'FontSize',20); set(get(gca,'title'),'FontSize',20); grid on %————————一层的绝对加速度和层间相对位移————————————————— figure plot(t,displacement(2,:)-displacement(1,:)); xlabel('时间t/s'),ylabel('层间相对位移/m'); title('一层层间相对位移'); set(gca,'FontSize',20); set(get(gca,'XLabel'),'FontSize',20); set(get(gca,'YLabel'),'FontSize',20); set(get(gca,'title'),'FontSize',20); grid on figure plot(t,absacceleration(2,:)); xlabel('时间t/s'),ylabel('绝对加速度m/s^2'); title('一层绝对加速度'); set(gca,'FontSize',20); set(get(gca,'XLabel'),'FontSize',20); set(get(gca,'YLabel'),'FontSize',20); set(get(gca,'title'),'FontSize',20); grid on %——————————第八层绝对加速度和层间相对位移———————————————— figure plot(t,displacement(9,:)-displacement(8,:)); xlabel('时间t/s'),ylabel('层间相对位移/m'); title('八层层间相对位移'); set(gca,'FontSize',20); set(get(gca,'XLabel'),'FontSize',20); set(get(gca,'YLabel'),'FontSize',20); set(get(gca,'title'),'FontSize',20); grid on figure plot(t,absacceleration(9,:)); xlabel('时间t/s'),ylabel('绝对加速度m/s^2'); title('八层绝对加速度'); set(gca,'FontSize',20); set(get(gca,'XLabel'),'FontSize',20); set(get(gca,'YLabel'),'FontSize',20); set(get(gca,'title'),'FontSize',20); grid on %———————————————————————————————————————— toc;
评论
    相关推荐