wave act on high rise buildings

조회 수: 4 (최근 30일)
Latifah Hanum
Latifah Hanum 2020년 2월 25일
댓글: Latifah Hanum 2020년 2월 25일
something's wrong in the plot, they say the vectors aren't the same. anyone can help fix? attach full code
%Simulasi Model%
clc;clear all;close all;
%Masukkan Data Gempa%
load('Cosine.txt');
ag1=([Cosine]);
tt1=0:0.02:15;
ndata1=length(ag1)-1;
Pilihan=menu('Silahkan Pilih Menu Dibawah ini',...
'MDOF 7 Lantai Tanpa Peredam',...
'MDOF 7 Lantai Dengan Peredam',...
'MDOF 4 Lantai Tanpa Peredam',...
'MDOF 4 Lantai Dengan Peredam',...
'Exit Matlab')
switch Pilihan;
case (1)
clc;
M1=input('Masukkan Massa Lantai 1(N)M1=');
M2= input('Masukkan Massa Lantai 2(N)M2=');
M3=input('Masukkan Massa Lantai 3(N)M3=');
M4=input('Masukkan Massa Lantai 4(N)M4=');
M5=input('Masukkan Massa Lantai 1(N)M5=');
M6=input('Masukkan Massa Lantai 1(N)M6=');
M7=input('Masukkan Massa Lantai 1(N)M7=');
K1=input('Masukkan Nilai Kekakuan Struktur(N/m)=');
K2=10000;K3=K2;K4=K2;K5=K2;K6=K2;K7=K2;K8=K2;
%Studi Kasus 1
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Studi Kasus 2
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Studi Kasus 3
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Pilih Beban Gempa
gempa=menu('Silahkan pilih',...
'Beban Cosine');
switch gempa;
case(1)
ag=ag1;
tt=tt1;
ndata=ndata1;
end
M=[M1 0 0 0 0 0 0
0 M2 0 0 0 0 0
0 0 M3 0 0 0 0
0 0 0 M4 0 0 0
0 0 0 0 M5 0 0
0 0 0 0 0 M6 0
0 0 0 0 0 0 M7];
K=[K1+K2 -K2 0 0 0 0 0
-K2 K2+K3 -K3 0 0 0 0
0 -K3 K3+K4 -K4 0 0 0
0 0 -K4 K4+K5 -K5 0 0
0 0 0 -K5 K5+K6 -K6 0
0 0 0 0 -K6 K6+K7 -K7
0 0 0 0 0 -K7 K7+K8];
psi=10;
[V,D]=eig(inv(M)*K);
W=sqrt(D),%disp(W);
alpha=(0.01*psi*2*V*M*V*W)/(V*K*V),%disp(alpha);
dt=0.02;
t=0;
x=[0
0
0
0
0
0
0];
v=[0
0
0
0
0
0
0];
perpindahan1(:,1)=x;
kecepatan1(:,1)=v*10;
l=ones(7,1);
for i=1:ndata;
t1=t;
x1=x;
v1=v;
f1=inv(M)*((-1)*M*l*ag(i)-(K*x1));
t2=t+dt/2;
x2=x+v1*dt/2;
v2=v+f1*dt/2;
f2=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(K*x2));
t3=t+dt/2;
x3=x+v2*dt/2;
v3=v+f2*dt/2;
f3=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(K*x3));
t4=t+dt/2;
x4=x+v3*dt/2;
v4=v+f3*dt/2;
f4=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(K*x4));
t5=t+dt/2;
x5=x+v4*dt/2;
v5=v+f4*dt/2;
f5=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(K*x5));
t6=t+dt/2;
x6=x+v5*dt/2;
v6=v+f5*dt/2;
f6=inv(M)*((-1)*M*l*(ag(i+1))-(K*x6));
x=x+1/6*dt*(v1+2*v2+2*v3+2*v4+2*v5+v6);
v=v+1/6*dt*(f1+2*f2+2*f3+2*f4+2*f5+f6);
t=t+0.02;
perpindahan1(:,i+1)=x;
kecepatan1(:,i+1)=v*10;
end
%Nilai Maksimum Dari Respon%
disp(max(perpindahan1(7,:)));
disp(max(perpindahan1(6,:)));
disp(max(perpindahan1(5,:)));
disp(max(perpindahan1(4,:)));
disp(max(perpindahan1(3,:)));
disp(max(perpindahan1(2,:)));
disp(max(perpindahan1(1,:)));
disp(min(perpindahan1(7,:)));
disp(min(perpindahan1(6,:)));
disp(min(perpindahan1(5,:)));
disp(min(perpindahan1(4,:)));
disp(min(perpindahan1(3,:)));
disp(min(perpindahan1(2,:)));
disp(min(perpindahan1(1,:)));
disp(max(abs(kecepatan1(1,:))));
disp(max(abs(ag)));
%Plot Respon Struktur Dalam Grafik%
figure(1);han=plot(tt,perpindahan1),xlabel('Waktu(detik)'),ylabel('Displacement(meter)'),grid on,
legend('Lantai 1','Lantai 2','Lantai 3','Lantai 4','Lantai 5','Lantai 6','Lantai 7');
figure(2);subplot(1,1,1),plot(tt,kecepatan1),xlabel('Waktu(detik)'),ylabel('Kecepatan(meter/s)'),grid on;
case(2)
clc;
M1=input('Masukkan Massa Lantai 1(N)M1=');
M2= input('Masukkan Massa Lantai 2(N)M2=');
M3=input('Masukkan Massa Lantai 3(N)M3=');
M4=input('Masukkan Massa Lantai 4(N)M4=');
M5=input('Masukkan Massa Lantai 1(N)M5=');
M6=input('Masukkan Massa Lantai 1(N)M6=');
M7=input('Masukkan Massa Lantai 1(N)M7=');
K1=input('Masukkan Nilai Kekakuan Struktur(N/m)=');
K2=10000;K3=K2;K4=K2;K5=K2;K6=K2;K7=K2;K8=K2;
%Studi Kasus 1
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Studi Kasus 2
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Studi Kasus 3
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Pilih Beban Gempa
gempa=menu('Silahkan pilih',...
'Beban Cosine');
switch gempa;
case(1)
ag=ag1;
tt=tt1;
ndata=ndata1;
end
M=[M1 0 0 0 0 0 0
0 M2 0 0 0 0 0
0 0 M3 0 0 0 0
0 0 0 M4 0 0 0
0 0 0 0 M5 0 0
0 0 0 0 0 M6 0
0 0 0 0 0 0 M7];
K=[K1+K2 -K2 0 0 0 0 0
-K2 K2+K3 -K3 0 0 0 0
0 -K3 K3+K4 -K4 0 0 0
0 0 -K4 K4+K5 -K5 0 0
0 0 0 -K5 K5+K6 -K6 0
0 0 0 0 -K6 K6+K7 -K7
0 0 0 0 0 -K7 K7+K8];
psi=10;
[V,D]=eig(inv(M)*K);
W=sqrt(D);,%disp(W);
alpha=(0.01*psi*2*V*M*V*W)/(V*K*V);,%disp(alpha);
C=alpha*K;
dt=0.02;
t=0;
x=[0
0
0
0
0
0
0];
v=[0
0
0
0
0
0
0];
perpindahan2(:,1)=x;
kecepatan2(:,1)=v*10;
l=ones(7,1);
for i=1:ndata;
t1=t;
x1=x;
v1=v;
f1=inv(M)*((-1)*M*l*ag(i)-(C*v1)-K*x1);
t2=t+dt/2;
x2=x+v1*dt/2;
v2=v+f1*dt/2;
f2=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(C*v2)-K*x2);
t3=t+dt/2;
x3=x+v2*dt/2;
v3=v+f2*dt/2;
f3=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(C*v3)-K*x3);
t4=t+dt/2;
x4=x+v3*dt/2;
v4=v+f3*dt/2;
f4=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(C*v4)-K*x4);
t5=t+dt/2;
x5=x+v4*dt/2;
v5=v+f4*dt/2;
f5=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(C*v5)-K*x5);
t6=t+dt/2;
x6=x+v5*dt/2;
v6=v+f5*dt/2;
f6=inv(M)*((-1)*M*l*(ag(i+1))-(C*v6)-K*x6);
x=x+1/6*dt*(v1+2*v2+2*v3+2*v4+2*v5+v6);
v=v+1/6*dt*(f1+2*f2+2*f3+2*f4+2*f5+f6);
t=t+0.02;
perpindahan2(:,i+1)=x;
kecepatan2(:,i+1)=v*10;
end
%Nilai Maksimum Dari Respon%
disp(max(perpindahan2(7,:)));
disp(max(perpindahan2(6,:)));
disp(max(perpindahan2(5,:)));
disp(max(perpindahan2(4,:)));
disp(max(perpindahan2(3,:)));
disp(max(perpindahan2(2,:)));
disp(max(perpindahan2(1,:)));
disp(min(perpindahan2(7,:)));
disp(min(perpindahan2(6,:)));
disp(min(perpindahan2(5,:)));
disp(min(perpindahan2(4,:)));
disp(min(perpindahan2(3,:)));
disp(min(perpindahan2(2,:)));
disp(min(perpindahan2(1,:)));
disp(max(abs(kecepatan2(1,:))));
disp(max(abs(ag)));
%Plot Respon Struktur Dalam Grafik%
figure(1);han=plot(tt,perpindahan2),xlabel('Waktu(detik)'),ylabel('Displacement(meter)'),grid on,
legend('Lantai 1','Lantai 2','Lantai 3','Lantai 4','Lantai 5','Lantai 6','Lantai 7');
figure(2);subplot(1,1,1),plot(tt,kecepatan2),xlabel('Waktu(detik)'),ylabel('Kecepatan(meter/s)'),grid on;
case(3)
clc;
M1=input('Masukkan Massa Lantai 1(N)M1=');
M2= input('Masukkan Massa Lantai 2(N)M2=');
M3=input('Masukkan Massa Lantai 3(N)M3=');
M4=input('Masukkan Massa Lantai 4(N)M4=');
K1=input('Masukkan Nilai Kekakuan Struktur(N/m)=');
K2=10000;K3=K2;K4=K2;K5=K2;
%Studi Kasus 1
%M1=1000;M2=M1;M3=1;M4=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;
%Studi Kasus 2
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Studi Kasus 3
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Pilih Beban Gempa
gempa=menu('Silahkan pilih',...
'Beban Cosine');
switch gempa;
case(1)
ag=ag1;
tt=tt1;
ndata=ndata1;
end
M=[M1 0 0 0
0 M2 0 0
0 0 M3 0
0 0 0 M4];
K=[K1+K2 -K2 0 0
-K2 K2+K3 -K3 0
0 -K3 K3+K4 -K4
0 0 -K4 K4+K5];
psi=10;
[V,D]=eig(inv(M)*K);
W=sqrt(D);,%disp(W);
alpha=(0.01*psi*2*V*M*V*W)/(V*K*V);,%disp(alpha)
dt=0.02;
t=0;
x=[0
0
0
0];
v=[0
0
0
0];
perpindahan3(:,1)=x;
kecepatan3(:,1)=v*10;
l=ones(4,1);
for i=1:ndata;
t1=t;
x1=x;
v1=v;
f1=inv(M)*((-1)*M*l*ag(i)-K*x1);
t2=t+dt/2;
x2=x+v1*dt/2;
v2=v+f1*dt/2;
f2=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-K*x2);
t3=t+dt/2;
x3=x+v2*dt/2;
v3=v+f2*dt/2;
f3=inv(M)*((-1)*M*l*(ag(i+1))-K*x3);
x=x+1/6*dt*(v1+2*v2+v3);
v=v+1/6*dt*(f1+2*f2+f3);
t=t+0.02;
perpindahan3(:,i+1)=x;
kecepatan3(:,i+1)=v*10;
end
%Nilai Maksimum Dari Respon%
disp(max(perpindahan3(4,:)));
disp(max(perpindahan3(3,:)));
disp(max(perpindahan3(2,:)));
disp(max(perpindahan3(1,:)));
disp(min(perpindahan3(4,:)));
disp(min(perpindahan3(3,:)));
disp(min(perpindahan3(2,:)));
disp(min(perpindahan3(1,:)));
disp(max(abs(kecepatan3(1,:))));
disp(max(abs(ag)));
%Plot Respon Struktur Dalam Grafik%
figure(1);han=plot(tt,perpindahan3),xlabel('Waktu(detik)'),ylabel('Displacement(meter)'),grid on,
legend('Lantai 1','Lantai 2','Lantai 3','Lantai 4');
figure(2);subplot(1,1,1),plot(tt,kecepatan3),xlabel('Waktu(detik)'),ylabel('Kecepatan(meter/s)'),grid on;
case(4)
clc;
M1=input('Masukkan Massa Lantai 1(N)M1=');
M2= input('Masukkan Massa Lantai 2(N)M2=');
M3=input('Masukkan Massa Lantai 3(N)M3=');
M4=input('Masukkan Massa Lantai 4(N)M4=');
K1=input('Masukkan Nilai Kekakuan Struktur(N/m)=');
K2=10000;K3=K2;K4=K2;K5=K2;
%Studi Kasus 1
%M1=1000;M2=M1;M3=1;M4=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;
%Studi Kasus 2
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Studi Kasus 3
%M1=1000;M2=M1;M3=1;M4=M1;M5=M1;
%K=10000;
%K1=K;
%K2=K;K3=K;K4=K;K5=K;
%Pilih Beban Gempa
gempa=menu('Silahkan pilih',...
'Beban Cosine');
switch gempa;
case(1)
ag=ag1;
tt=tt1;
ndata=ndata1;
end
M=[M1 0 0 0
0 M2 0 0
0 0 M3 0
0 0 0 M4];
K=[K1+K2 -K2 0 0
-K2 K2+K3 -K3 0
0 -K3 K3+K4 -K4
0 0 -K4 K4+K5];
psi=10;
[V,D]=eig(inv(M)*K);
W=sqrt(D);%disp(W)
alpha=(0.01*psi*2*V*M*V*W)/(V*K*V);,%disp(alpha);
C=alpha*K;
dt=0.02;
t=0;
x=[0
0
0
0];
v=[0
0
0
0];
perpindahan4(:,1)=x;
kecepatan4(:,1)=v*10;
l=ones(4,1);
for i=1:ndata;
t1=t;
x1=x;
v1=v;
f1=inv(M)*((-1)*M*l*ag(i)-(C*v1)-K*x1);
t2=t+dt/2;
x2=x+v1*dt/2;
v2=v+f1*dt/2;
f2=inv(M)*((-1)*M*l*0.5*(ag(i)+ag(i+1))-(C*v2)-(K*x2));
t3=t+dt/2;
x3=x+v2*dt/2;
v3=v+f2*dt/2;
f3=inv(M)*((-1)*M*l*(ag(i+1))-(C*v3)-K*x3);
x=x+1/6*dt*(v1+2*v2+v3);
v=v+1/6*dt*(f1+2*f2+f3);
t=t+0.02;
perpindahan4(:,i+1)=x;
kecepatan4(:,i+1)=v*10;
end
%Nilai Maksimum Dari Respon%
disp(max(perpindahan4(4,:)));
disp(max(perpindahan4(3,:)));
disp(max(perpindahan4(2,:)));
disp(max(perpindahan4(1,:)));
disp(min(perpindahan4(4,:)));
disp(min(perpindahan4(3,:)));
disp(min(perpindahan4(2,:)));
disp(min(perpindahan4(1,:)));
disp(max(abs(kecepatan4(1,:))));
disp(max(abs(ag)));
%Plot Respon Struktur Dalam Grafik%
figure(1);han=plot(tt,perpindahan4),xlabel('Waktu(detik)'),ylabel('Displacement(meter)'),grid on,
legend('Lantai 1','Lantai 2','Lantai 3','Lantai 4');
figure(2);subplot(1,1,1),plot(tt,kecepatan4),xlabel('Waktu(detik)'),ylabel('Kecepatan(meter/s)'),grid on;
end
  댓글 수: 8
이전 댓글 6개 표시
darova
darova 2020년 2월 25일
Is it possible edit your post so that the code will be attached? I don't have time to scroll up and down all the time
Latifah Hanum
Latifah Hanum 2020년 2월 25일
Yes, I wanna try, thankyou

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

답변 (0개)

카테고리

Help Center 및 File Exchange에서 MATLAB에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by