clear all clc % 3 Story Building (prob 8-11 in Hurty and Rubenstein) except that the % the masses are numbered from story 1 to story 3, 1,2,3 upward % Stiffness and mass matricies in kips,feet and sec. % The base motion is ug(t)=A1*sin(omega*t) ft. This is not the steady % -state solution as given by the phasor (frequency response) method. % The structure will never reach a steady state sinusoidal response as long % as the damping is zero. ndof=3; % number of degrees of freedom disp('The mass matrix is') m=[1.5 0 0;0 1 0;0 0 1.5]*38.4 disp('The stiffness matrix is') k=[126 -58 0;-58 100 -42;0 -42 42]*1000 disp('The damping matrix is') c=[0 0 0;0 0 0;0 0 0] %c=.01*k D=inv(m)*k; [PHI1,lambda]=eig(D); omegan1=sqrt(diag(lambda))'; [omegan1, index]=sort((omegan1')); disp('The radian undamped natural frequencies are') omegan=omegan1' Phi=zeros(ndof,ndof); for i=1:ndof PHI(:,i)=PHI1(:,index(i,1)); end disp('The modal matrix is') PHI disp('The Hertzian undamped natural frequencies are') fn=omegan/(2*pi) disp('The undamped natural periods are') Tn=ones(1,3)./fn% natural periods in sec. M=PHI'*m*PHI; K=PHI'*k*PHI; C=PHI'*c*PHI; for i=1:ndof zeta(1,i)=C(i,i)/(2*sqrt(K(i,i)*M(i,i))); omegad(1,i)=omegan(1,i)*sqrt(1-zeta(1,i)^2); end zeta omegad dt=.01; disp('The ground motion is ug(t)=A1*sin(omega*t)'); A1=input(' Ground motion amplitude is A1 = '); omega=input(' Ground motion radian frequency is omega = '); disp(' ') disp('Press Enter to Continue') pause t=[0:dt:400*dt]; ugdd=-A1*(omega^2)*sin(omega*t);%ground acceleration ug=A1*sin(omega*t);%ground motion ugd=A1*omega*cos(omega*t);%ground acceleration minv=inv(m); p=-[m(1,1);m(2,2);m(3,3)]*ugdd; pp1=-[m(1,1);m(2,2);m(3,3)]*[ugdd(1,2:401) A1*(omega^2)*sin(omega*402*dt)]; P=PHI'*p;% PP1=PHI'*pp1; %ut=u+ug uzero=zeros(3,1); Q=zeros(3,401); udotzero=-A1*omega*[1;1;1]; Q(:,1)=uzero; qzero=inv(PHI)*uzero;qdotzero=inv(PHI)*udotzero; for i=1:ndof x=[qzero(i,1);qdotzero(i,1)]; z=zeta(1,i); wn=omegan(1,i); wD=omegad(1,i); Co=cos(wD*dt);Si=sin(wD*dt); E=exp(-z*wn*dt);H=z/sqrt(1-z*z);J=sqrt(1-z*z); A=E*(Co+H*Si);B=E*Si/wD;Ap=-E*((wn/J)*Si);Bp=E*(Co-H*Si); phi=[A B;Ap Bp];%fill the transition matrix C=(1/K(i,i))*(((2*z)/(wn*dt))+E*((((1-2*z*z)/(wD*dt))-H)*Si-(1+((2*z)/(wn*dt)))*Co)); D=(1/K(i,i))*(1-((2*z)/(wn*dt))+E*(((2*z*z-1)/(wD*dt))*Si+((2*z/(wn*dt)))*Co)); Cp=(1/K(i,i))*(-(1/dt)+E*(((wn/J)+(H/dt))*Si+(1/dt)*Co)); Dp=(1/(K(i,i)*dt))*(1-E*(H*Si+Co)); gamma=[C D;Cp Dp];%fill the input distribution matrix for k1=2:400 x=phi*x+gamma*[P(i,k1-1);PP1(i,k1-1)]; Q(i,k1)=x(1,1); end end U=PHI*Q; Ut=U+[ug;ug;ug]; figure(1);clf plot(t,ug) xlabel('Time t,sec.') ylabel('Base Motion, u_g, ft.') text(200*dt,.08, 'Press Enter to Continue') pause figure(2);clf plot(t,ug,t,Ut(1,:),t,Ut(2,:),t,Ut(3,:)) xlabel('Time t,sec.') ylabel('Displacement, ft.') legend('u_g(t)','ut_1(t)','ut_2(t)','ut_3(t)') Umax=max(max(abs(Ut(1:3,:)))); text(t(1,190),-.8*Umax,['Max. Total Response = ' num2str(Umax)]); text(t(1,190),-.9*Umax, 'Press Enter to Continue') pause U1=U*5/Umax;%scale relative displacement for animation ug1=ug*5/Umax;% scaled ground motion box3x=[ 0 30 30 0 0];box3y=[34 34 36 36 34]; box2x=[ 0 30 30 0 0];box2y=[22 22 24 24 22]; box1x=[ 0 30 30 0 0];box1y=[10 10 12 12 10]; col1x=zeros(1,7);col1y=[0 1 3 5 7 9 10]; col2x=col1x+30*ones(1,7);col2y=col1y; col3x=col1x;col3y=col1y+12*ones(1,7); col4x=col2x;col4y=col3y; col5x=col1x;col5y=col1y+24*ones(1,7); col6x=col2x;col6y=col5y; base1x=[-5 35 35 -5 -5];base1y=[0 0 -4 -4 0]; damp1x=[3 3 8];damp1y=[0 5 5]; damp2x=[22 8 8 22];damp2y=[6 6 4 4]; damp3x=[15 15];damp3y=[6 4]; damp4x=[27 27 15];damp4y=[10 5 5]; arrowx=[0 6 4.5 4.5 6];arrowy=[0 0 1 -1 0]; figure(3);clf% animate the displacement excited motion axis square;hold on axis off;hold on axis([-10 40 -10 40]);hold on Base1=patch(base1x,base1y,'r','EdgeColor','k','EraseMode','xor');hold on B1=patch(box1x,box1y,'r','EraseMode','xor');hold on B2=patch(box2x,box2y,'r','EraseMode','xor');hold on B3=patch(box3x,box3y,'r','EraseMode','xor');hold on C1=plot(col1x,col1y,'k','EraseMode','xor');hold on C2=plot(col2x,col2y,'k','EraseMode','xor');hold on C3=plot(col3x,col3y,'k','EraseMode','xor');hold on C4=plot(col4x,col4y,'k','EraseMode','xor');hold on C5=plot(col5x,col5y,'k','EraseMode','xor');hold on C6=plot(col6x,col6y,'k','EraseMode','xor');hold on D1=plot(damp1x,damp1y,'k');hold on D2=plot(damp2x,damp2y,'k');hold on D3=plot(damp3x,damp3y,'k');hold on D4=plot(damp4x,damp4y,'k');hold on D21=plot(damp1x,damp1y+[ 1 1 1 ]*12,'k');hold on D22=plot(damp2x,damp2y+[ 1 1 1 1]*12,'k');hold on D23=plot(damp3x,damp3y+[ 1 1 ]*12,'k');hold on D24=plot(damp4x,damp4y+[ 1 1 1 ]*12,'k');hold on D31=plot(damp1x,damp1y+[ 1 1 1]*24,'k');hold on D32=plot(damp2x,damp2y+[ 1 1 1 1]*24,'k');hold on D33=plot(damp3x,damp3y+[ 1 1 ]*24,'k');hold on D34=plot(damp4x,damp4y+[ 1 1 1]*24,'k');hold on Arrow1=plot(arrowx+31*ones(1,5),arrowy+10*ones(1,5)); Arrow2=plot(arrowx+31*ones(1,5),arrowy+22*ones(1,5)); Arrow3=plot(arrowx+31*ones(1,5),arrowy+34*ones(1,5)); Arrow4=plot(arrowx+32*ones(1,5),arrowy+ones(1,5)); Tx1=text(15,14,'m_1');Tx2=text(15,26,'m_2');Tx3=text(15,38,'m_3'); Tx4=text(31,5,'k_1/2');Tx5=text(31,17,'k_2/2');Tx6=text(31,29,'k_3/2'); Tx8=text(32,8,'u_1');Tx9=text(32,20,'u_2');Tx10=text(32,32,'u_3'); Tx11=text(15,7.5,'c_1');Tx12=text(15,19.5,'c_2');Tx13=text(15,31.5,'c_3'); Txt14=text(38,0,'u_g(t)'); Tx7=text(-5,-7,'Press Enter to Set Initial Displacements'); pause clf axis square;hold on axis off;hold on axis([-10 40 -10 40]);hold on i=1; e=[0 0 .15 .5 .85 1 1]; deltax1=U1(1,1); deltax2=U1(2,1)-U1(1,1); deltax3=U1(3,1)-U1(2,1); Base1=patch(base1x,base1y,'r','EdgeColor','k','EraseMode','xor');hold on B1=patch(box1x+U1(1,i)*ones(1,5),box1y,'r','EraseMode','xor');hold on B2=patch(box2x+U1(2,i)*ones(1,5),box2y,'r','EraseMode','xor');hold on B3=patch(box3x+U1(3,i)*ones(1,5),box3y,'r','EraseMode','xor');hold on C1=plot(col1x+(deltax2)*e+(ug1(1,1)+deltax1)*ones(1,7),col1y,'k','EraseMode','xor');hold on C2=plot(col2x+(deltax2)*e+(ug1(1,1)+deltax1)*ones(1,7),col2y,'k','EraseMode','xor');hold on C3=plot(col3x+deltax2*e+(ug1(1,1)+deltax2)*ones(1,7),col3y,'k','EraseMode','xor');hold on C4=plot(col4x+deltax2*e+(ug1(1,1)+deltax2)*ones(1,7),col4y,'k','EraseMode','xor');hold on C5=plot(col5x+deltax3*e+(ug1(1,1)+deltax3)*ones(1,7),col5y,'k','EraseMode','xor');hold on C6=plot(col6x+deltax3*e+(ug1(1,1)+deltax3)*ones(1,7),col6y,'k','EraseMode','xor');hold on D1=plot(damp1x+ug1(1,i)*[1 1 1],damp1y,'k');hold on D2=plot(damp2x+ug1(1,i)*[1 1 1 1],damp2y,'k');hold on D3=plot(damp3x+[1 1]*U1(1,i),damp3y,'k');hold on D4=plot(damp4x+[1 1 1 ]*U1(1,i),damp4y,'k');hold on D21=plot(damp1x+[1 1 1]*U1(1,i),damp1y+[ 1 1 1 ]*12,'k');hold on D22=plot(damp2x+[1 1 1 1]*U1(1,i),damp2y+[ 1 1 1 1]*12,'k');hold on D23=plot(damp3x+[1 1]*U1(2,i),damp3y+[ 1 1 ]*12,'k');hold on D24=plot(damp4x+[1 1 1]*U1(2,i),damp4y+[ 1 1 1 ]*12,'k');hold on D31=plot(damp1x+[1 1 1]*U1(2,i),damp1y+[ 1 1 1]*24,'k');hold on D32=plot(damp2x+[1 1 1 1]*U1(2,i),damp2y+[ 1 1 1 1]*24,'k');hold on D33=plot(damp3x+[1 1]*U1(3,i),damp3y+[ 1 1 ]*24,'k');hold on D34=plot(damp4x+[1 1 1]*U1(3,1),damp4y+[ 1 1 1]*24,'k');hold on Tx7=text(-7,-7,['Press Enter to Animate Max. Total Resp. = ' num2str(Umax) ' ft.']); Tx13=text(-7,-9,['t = ' num2str(0) ' sec.']); pause % now animate for i=1:401 deltax1=U1(1,i); deltax2=U1(2,i)-U1(1,i); deltax3=U1(3,i)-U1(2,i); set(Base1,'Xdata',base1x+ug1(1,i)*[1 1 1 1 1]); set(B1,'Xdata',box1x+(U1(1,i)+ug1(1,i))*ones(1,5));hold on set(B2,'Xdata',box2x+(U1(2,i)+ug1(1,i))*ones(1,5));hold on set(B3,'Xdata',box3x+(U1(3,i)+ug1(1,i))*ones(1,5));hold on set(C1,'Xdata',col1x+deltax1*e+ug1(1,i)*ones(1,7));hold on set(C2,'Xdata',col2x+(deltax1)*e+(ug1(1,i))*ones(1,7));hold on set(C3,'Xdata',col3x+deltax2*e+(ug1(1,i)+U1(1,i))*ones(1,7));hold on set(C4,'Xdata',col4x+deltax2*e+(ug1(1,i)+U1(1,i))*ones(1,7));hold on set(C5,'Xdata',col5x+deltax3*e+(ug1(1,i)+U1(2,i))*ones(1,7));hold on set(C6,'Xdata',col6x+deltax3*e+(ug1(1,i)+U1(2,i))*ones(1,7));hold on set(D1,'Xdata',damp1x+ug1(1,i)*ones(1,3));hold on set(D2,'Xdata',damp2x+ug1(1,i)*ones(1,4));hold on set(D3,'Xdata',damp3x+(ug1(1,i)+U1(1,i))*ones(1,2)); set(D4,'Xdata',damp4x+(ug1(1,i)+U1(1,i))*ones(1,3)); set(D21,'Xdata',damp1x+(ug1(1,i)+U1(1,i))*ones(1,3)); set(D22,'Xdata',damp2x+(ug1(1,i)+U1(1,i))*ones(1,4)); set(D23,'Xdata',damp3x+(ug1(1,i)+U1(2,i))*ones(1,2)); set(D24,'Xdata',damp4x+(ug1(1,i)+U1(2,i))*ones(1,3)); set(D31,'Xdata',damp1x+(ug1(1,i)+U1(2,i))*ones(1,3)); set(D32,'Xdata',damp2x+(ug1(1,i)+U1(2,i))*ones(1,4)); set(D33,'Xdata',damp3x+(ug1(1,i)+U1(3,i))*ones(1,2)); set(D34,'Xdata',damp4x+(ug1(1,i)+U1(3,i))*ones(1,3)); set(Tx13,'String',['t = ' num2str(t(1,i)) ' sec.']); pause(.04) end