Ukaolv.zip

%%%%———结构动力学多自由度体系计算，包括隔震支座，采用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('C:\Users\admin\Desktop\结构动力学project\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); %————————初始量赋值—————————— 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); %—————————————————————————————————————— dp(:,i)=(p(:,i+1)-p(:,i))+Gaa*velocity(:,i)+Gbb*acceleration(:,i); dq(:,i)=GKK\dp(:,i); dq1(:,i)=gama*dq(:,i)/(beta*dt)-gama*velocity(:,i)/beta+dt*acceleration(:,i)*(1-gama/(2*beta)); dq2(:,i)=dq(:,i)/(beta*dt^2)-velocity(:,i)/(beta*dt)-acceleration(:,i)/(2*beta); 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; 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;