clc
clear all
fprintf('-----------------STRUCTURAL DYNAMIC-------------------------------\n')
fprintf('-----------------Response Spectrum Analysis------------------------\n')
N=6; %number of lumped Masses/storey's
H=3.2; %height of each storey
M=[373710 0 0 0 0 0;
0 373710 0 0 0 0;
0 0 373710 0 0 0;
0 0 0 373710 0 0;
0 0 0 0 269040 0;
0 0 0 0 0 26510]; %Mass Matrix
K=[1152600000 -576300000 0 0 0 0;
-576300000 1152600000 -576300000 0 0 0;
0 -576300000 1152600000 -576300000 0 0;
0 0 -576300000 1152600000 -576300000 0;
0 0 0 -113200000 689500000 -113200000;
0 0 0 0 -113200000 113200000]; %stiffness matrix
%Finding Eigen Vector and Eigen Value
fprintf('\nEigen Value and Eigen Vector :');
[Modeshape,omega]=eig(K,M);
fprintf('\nEigen Vector\n');
disp(Modeshape);
fprintf('\nEigen Value\n');
disp(omega);
m=diag(M);
%Modal Expansion of effective EQ forces (Chopra pg.no.551)
for j=1:N
for i=1
num(i,j)=m(i)*(Modeshape(i,j));
den(i,j)=m(i)*Modeshape(i,j)*Modeshape(i,j);
end
for i=2:N
num(i,j)= num(i-1,j)+m(i)*(Modeshape(i,j));
den(i,j)=den(i-1,j)+m(i)*Modeshape(i,j)*Modeshape(i,j);
end
mod_exp(j)=num(N,j)/den(N,j);
end
%Modal Mass calculation from IS1893
for i=1:1:N
modal_mass(i)=((num(N,i)*9.81)^2/(9.81*9.81*den(N,i)))
end
%Calculation of Time Peritod of each Mode
for i=1:1:N
Timemode(i)=(2*pi)/sqrt(omega(i,i))
end
importance_fact=input('\nImportance factor Of Building:')
response_red_fact=input('\nEnter Response Reduction Factor R:')
zonefact=input('\nEnter Zone Factor Z:')
selection=input('\nSelect any one method for calculation of base shear\n 1)IS Code Method \n2)Elecentro Ground motion')
if selection==1
%Obtaining Sa/g Values for Diffrent terrains using IS1893
fprintf('\nSelect Type of Terrain:\n');
fprintf('\n1)Rocky\n2)Medium Soils\n3)Soft Soils')
sel_terrain=input('\nEnter Type of Terrain 1,2 or 3')
%According to IS1893
if sel_terrain==1
for i=1:N
if Timemode(i)>=0 && Timemode(i)<=0.10
Sabyg(i)=1+15*Timemode(i);
elseif Timemode(i)>=0.10 && Timemode(i)<=0.40
Sabyg(i)=2.5;
elseif Timemode(i)>=0.40 && Timemode(i)<=4.0
Sabyg(i)=1/Timemode(i);
end
end
elseif sel_terrain==2
for i=1:N
if Timemode(i)>=0 && Timemode(i)<=0.10
Sabyg(i)=1+15*Timemode(i);
elseif Timemode(i)>=0.10 && Timemode(i)<=0.40
Sabyg(i)=2.5;
elseif Timemode(i)>=0.40 && Timemode(i)<=4.0
Sabyg(i)=1.36/Timemode(i);
end
end
else
for i=1:N
if Timemode(i)>=0 && Timemode(i)<=0.10
Sabyg(i)=1+15*Timemode(i);
elseif Timemode(i)>=0.10 && Timemode(i)<=0.40
Sabyg(i)=2.5;
elseif T(i)>=0.40 && Timemode(i)<=4.0
Sabyg(i)=1.67/Timemode(i);
end
end
end
else
%Calculation of Sabyg Starts - using elcentro motion
K = 45*10^4;
%T = 31.18;
Z = 0.05; % Case 1 ( 5% Damping)
dt = 0.02;
matno=1;
maxpga=0.515;
for lkn=1:1:1
load elecentro.txt;
a1=elecentro;
ug1_pga=max(a1);
a1=a1.*(maxpga/ug1_pga);
ts=0.02;
fprintf('1')
T=0.01:0.01:4; % building time period
for i=1:1:length(T)
U(2) = 0; % SD
V(2) = 0; % SV
V1(2) = 0; % PSV
A1(2) = 0; % PSA
W(i) = ((2*pi)./T(i));
m(i) = ((K)./(W(i).^2));
P = m(i).*(a1);
c(i) = (2*Z.*m(i).*(W(i)));
K1(i) = (m(i)./(dt^2))+(c(i)./(2*dt));
a(i) = (m(i)./(dt^2))-(c(i)./(2*dt));
b(i) = K-(2*(m(i))./(dt^2));
A(2) = ((P(1)-(c(i).*V(2))-(K.*U(2)))./m(i)); % SA
U(1) = U(2)-(dt.*V(2))+(A(2).*(dt^2)/2);
t=0:0.02:31.18;
for j=2:1:length(t)-1
P7(j) = P(j)-(a(i).*U(j-1))-(b(i).*U(j));
U(j+1) = P7(j)./K1(i); % SD
V(j+1) = ((U(j+1)-U(j-1))./(2*dt)); % SV
V1(j+1)= U(j+1).*W(i); % PSV
A(j+1) = (((U(j+1)-(2*U(j))+U(j-1)))./(dt^2)); % SA
A1(j+1)= V(j+1).*(W(i)); % PSA
end
%display(U);
% Maximum Values:
Umax(i)=max(abs(U));
Vmax(i)=max(abs(V));
V1max(i)=Umax(i).*W(i);
Amax(i)=max(abs(A));
A1max(i)=V1max(i).*W(i);
Amaxall(:,matno)=mean(Amax(i));
matno=matno+1;
end
end
meanamaxall=Amaxall;
stddevumaxall=std(Amaxall,0,1);
plotamaxall=meanamaxall+stddevumaxall;
T=0.01:0.01:4;
plot(T,plotamaxall);
title('Design Spectrum Diagram');
xlabel('Time period')
ylabel('Spectral Acceleration in m')
hold off
for j=1:1:N
Sabyg(j)=plotamaxall((round(Timemode(j).*100)));
end
%Calculation of Sa/g Ends
end
%Horizontal Seismic Coefficient Ah
for i=1:1:N
Ah(i)=zonefact*importance_fact*Sabyg(i)/(2*response_red_fact);
end
%Design Lateral Force at each floor in each mode
for k=1:N
for i=1:N
Qik(i,k)=(Ah(k)*Modeshape(i,k)*mod_exp(k)*modal_mass(k)*9.81)*1*10^-4
end
end
%Storey shear forces in each mode
for j=1:N
for i=1:N
if i==1
storey_shear(i,j)=Qik(1,j)+Qik(2,j)+Qik(3,j)+Qik(4,j)+Qik(5,j)+Qik(6,j)
end
if i==2
storey_shear(i,j)=Qik(2,j)+Qik(3,j)+Qik(4,j)+Qik(5,j)+Qik(6,j)
end
if i==3
storey_shear(i,j)=Qik(3,j)+Qik(4,j)+Qik(5,j)+Qik(6,j)
end
if i==4
storey_shear(i,j)=Qik(4,j)+Qik(5,j)+Qik(6,j)
end
if i==5
storey_shear(i,j)=Qik(5,j)+Qik(6,j)
end
if i==6
storey_shear(i,j)=Qik(6,j)
end
end
end
%Using SRSS Method to find Peak Modal Responses of each storey
for i=1:N
peak_storey_shear(i)=sqrt(storey_shear(i,1)^2+storey_shear(i,2)^2+storey_shear(i,3)^2+storey_shear(i,4)^2+storey_shear(i,5)^2)
end
%Lateral Forces at each storey due to all modes considered
for i=1:N
if i~=N%For non terrace floors
force(i)=peak_storey_shear(i)-peak_storey_shear(i+1)
else%For terrace/roof floor
force(i)=peak_storey_shear(i)
end
end
%calculation for base shear
for i=1:N
if i==1%For ground floor
base_shear=force(1)
else
base_shear=base_shear+force(i)
end
end
for i=1:N
if i==1
base_moment=force(1)*H
elseif i>1
base_moment=base_moment+force(i)*H*i
end
end
tot_mass=sum(diag(M));
base_reaction=tot_mass/950
display(base_shear)
display(base_moment)
display(base_reaction)
% % % output:
%base_shear =
%
% 116.1745
%
%
% base_moment =
%
% 1.0834e+003
%
%
% base_reaction =
%
% 1.8846e+003
