【信号处理】地面震动作用下多层建筑的线性分析附matlab代码

最新推荐文章于 2024-05-25 23:53:59 发布

天天Matlab代码科研顾问

最新推荐文章于 2024-05-25 23:53:59 发布

阅读量12

点赞数

文章标签：信号处理 matlab 算法人工智能前端

本文链接：https://blog.csdn.net/m0_57702748/article/details/140031780

版权

1 简介

地面震动作用下多层建筑的线性分析附matlab代码

2 完整代码


          
          
           
           clear; clc; close all;
           
           
%Input
           
           
N=4;                   %Number of stories
           
           
Ks=1000;               %Interstory stiffness
           
           
Ms=5;                  %Floor mass
           
           
fs=100;                %Sampling frequency
           
           
ag=randn(1,10000);     %Ground acceleration
           
           
%----------------------------
           
           
%Processing
           
           
M=Ms*diag(ones(N,1));                  %Mass matrix
           
           
K=MultiStory_Stiffness(Ks,N);          %Stiffness matrix
           
           
C=0.01*M+0.01*K;                       %Damping matrix
           
           
f=-M*ones(N,1).*ag;                    %Effective force
           
           
Result=MDOF_simulation(M,C,K,f,fs);    %MDOF solver
           
           
%----------------------------
           
           
%Plot Response
           
           
t=[0:1/fs:(length(ag)-1)/fs];
           
           
figure;
           
           
for i=1:1:N
           
           
subplot(N,3,3*i-2)
           
           
plot(t,Result.Acceleration(i,:)); xlabel('Time (sec)'); ylabel(strcat('ACC',num2str(i)));
           
           
subplot(N,3,3*i-1)
           
           
plot(t,Result.Velocity(i,:)); xlabel('Time (sec)'); ylabel(strcat('VEL',num2str(i)));
           
           
subplot(N,3,3*i)
           
           
plot(t,Result.Displacement(i,:)); xlabel('Time (sec)'); ylabel(strcat('DSP',num2str(i)));
           
           
end
           
           
%---------------------------
           
           
%Plot Modeshapes
           
           
figure;
           
           
for i=1:1:N
           
           
plot([0 ; Result.Parameters.ModeShape(:,i)],0:N,'*-');
           
           
hold on
           
           
end
           
           
hold off
           
           
xlabel('Amplitude');
           
           
ylabel('Floor');
           
           
grid on
           
           
daspect([1 1 1]);
           
           
title('Modeshapes');


          
          
           
           function Result=MDOF_simulation(M,C,K,f,fs)
           
           
n=size(f,1);
           
           
dt=1/fs; %sampling rate
           
           
[Vectors, Values]=eig(K,M);
           
           
Freq=sqrt(diag(Values))/(2*pi); % undamped natural frequency
           
           
steps=size(f,2);
           
           
Mn=diag(Vectors'*M*Vectors); % uncoupled mass
           
           
Cn=diag(Vectors'*C*Vectors); % uncoupled damping
           
           
Kn=diag(Vectors'*K*Vectors); % uncoupled stifness
           
           
wn=sqrt(diag(Values));
           
           
zeta=Cn./(2*sqrt(Mn.*Kn));   % damping ratio
           
           
wd=wn.*sqrt(1-zeta.^2);
           
           
fn=Vectors'*f; % generalized input force matrix
           
           
t=[0:dt:dt*steps-dt];
           
           
%forced vibration
           
           
for i=1:1:n
           
           
    h(i,:)=(1/(Mn(i)*wd(i))).*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t); %transfer function of displacement
           
           
    hd(i,:)=(1/(Mn(i)*wd(i))).*(-zeta(i).*wn(i).*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t)+wd(i).*exp(-zeta(i)*wn(i)*t).*cos(wd(i)*t)); %transfer function of velocity
           
           
    hdd(i,:)=(1/(Mn(i)*wd(i))).*((zeta(i).*wn(i))^2.*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t)-zeta(i).*wn(i).*wd(i).*exp(-zeta(i)*wn(i)*t).*cos(wd(i)*t)-wd(i).*((zeta(i).*wn(i)).*exp(-zeta(i)*wn(i)*t).*cos(wd(i)*t))-wd(i)^2.*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t)); %transfer function of acceleration
           
           
    qq=conv(fn(i,:),h(i,:))*dt;
           
           
    qqd=conv(fn(i,:),hd(i,:))*dt;
           
           
    qqdd=conv(fn(i,:),hdd(i,:))*dt;
           
           
    q(i,:)=qq(1:steps); % modal displacement
           
           
    qd(i,:)=qqd(1:steps); % modal velocity
           
           
    qdd(i,:)=qqdd(1:steps); % modal acceleration
           
           
end
           
           
x=Vectors*q; %displacement
           
           
v=Vectors*qd; %vecloity
           
           
a=Vectors*qdd; %vecloity
           
           
% Free vibration
           
           
xi=zeros(n,1);  % displacement initial condition
           
           
vi=zeros(n,1);  % velocity initial condition
           
           
xno=Vectors'*M*xi./Mn;
           
           
vno=Vectors'*M*vi./Mn;
           
           
for i=1:1:n
           
           
AA=(vno(i)+xno(i).*zeta(i).*wn(i))./wd(i);
           
           
BB=xno(i);
           
           
qf(i,:)=exp(-zeta(i)*wn(i)*t).*(AA.*sin(wd(i)*t)+BB.*cos(wd(i)*t));
           
           
qdf(i,:)=wd(i)*exp(-zeta(i)*wn(i)*t).*(AA.*cos(wd(i)*t)-BB.*sin(wd(i)*t))-zeta(i).*wn(i).*exp(-zeta(i)*wn(i)*t).*(AA.*sin(wd(i)*t)+BB.*cos(wd(i)*t));
           
           
qddf(i,:)=wd(i)^2*exp(-zeta(i)*wn(i)*t).*(-AA.*sin(wd(i)*t)-BB.*cos(wd(i)*t))-2*zeta(i).*wn(i).*wd(i).*exp(-zeta(i)*wn(i)*t).*(AA.*cos(wd(i)*t)-BB.*sin(wd(i)*t))+zeta(i)^2.*wn(i)^2*exp(-zeta(i)*wn(i)*t).*(-AA.*sin(wd(i)*t)-BB.*cos(wd(i)*t));
           
           
end
           
           
x=x+Vectors*qf;
           
           
v=v+Vectors*qdf;
           
           
a=a+Vectors*qddf;
           
           
Result.Displacement=x;
           
           
Result.Velocity=v;
           
           
Result.Acceleration=a;
           
           
Result.Parameters.Freq=Freq;
           
           
Result.Parameters.DampRatio=zeta*100;
           
           
Result.Parameters.ModeShape=Vectors;
           
           
end


          
          
           
           function K=MultiStory_Stiffness(Ks,N)
           
           
%Input
           
           
%Ks: Interstory stiffness of columns
           
           
%N: Number of storys
           
           
%Output
           
           
%K: Stiffness matrix
           
           
KK=Ks*ones(N+1,1);
           
           
KK(N+1)=0;
           
           
K=zeros(N,N);
           
           
for i=1:1:N
           
           
for j=1:1:N              
           
           
if i==j 
           
           
K(i,j)=KK(i)+KK(i+1);
           
           
end            
           
           
if i==j+1
           
           
K(i,j)=-KK(i);
           
           
end      
           
           
if j==i+1
           
           
K(i,j)=-KK(j);
           
           
end        
           
           
end
           
           
end
           
           
end