clear all clc disp(' This script evaluated the forced response of a three story steel') disp(' frame from Problem 8.11 in Dymanics of Structures by ') disp(' Hurty and Rubenstein, Prentice-Hall, 1964') disp(' ') disp(' The matricies are in kips,feet and sec. ') disp(' Solved by using he interpolation of exciation technique') disp(' Press Enter to Continue') pause 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 alpha=[0 0]; % Damping matrix is a linear combination of the mass and stiffness matrices c=alpha(1,1)*m+alpha(1,2)*k A=inv(m)*k; [PHI1,lambda]=eig(A); omegan=sqrt(diag(lambda))'; [omegan, index]=sort((omegan)); disp(' The radian undamped natural frequencies are') omegan PHI=zeros(ndof,ndof); % reorder eigenvectors for i=1:ndof PHI(:,i)=PHI1(:,index(1,i)); end % normalize eigenvectrs with respect to tne mass matrix for i=1:ndof fact=PHI(:,i)'*m*PHI(:,i); PHI(:,i)=PHI(:,i)/sqrt(fact); end disp(' The ModalMatrix 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. C=PHI'*c*PHI; K=PHI'*k*PHI; z=zeros(1,ndof);omegad=zeros(1,ndof); for i=1:ndof z(1,i)=C(i,i)/(2*omegan(1,i)); omegad(1,i)=omegan(1,i)*sqrt(1-z(1,i)^2); end disp(' The radian damped natural frequencies are') omegad disp(' The modal damping ratios are') zeta=z disp(' Press Enter to Continue') pause figure(1);clf axis square;hold on; axis off; hold on basex1=0.5*[-20 20 18 14 7 -3 -8 -20]; basey1=[0 0 -11 -7 -15 -9 -12 0]; z1=[0 20 40 60]; plot(100*[zeros(3,1) PHI'],z1) axis([-30 30 -20 80]); legend('i=1','i=2','i=3') plot([0 0],[0 65],'k'); text(-25,65,'Mode Shapes'); patch(basex1,basey1,'r');hold on text(-25, 110,'Mode Shapes') text(-25,-18,'Press Enter to Continue'); pause uzero=zeros(3,1); %uzero=uzero';% initial displacements dt=.01; minv=inv(m); t=[0:dt:200*dt]; % bomb blast forcing function p2=6000*exp(-10*t)-3000*exp(-5*t);% force applied to 2nd and 3rd floors %p1=0.5*(p2+[p2(1,2:201) -.0001]);% use average forcing function over the interval p=[p2;p2;p2/2]; pp1=[p(:,2:201) [0;0;0]]; figure(2);clf plot(t,p2) xlabel('Time t,sec.') ylabel('Second Floor Force, p_1(t)') text(1,.8*max(max(p)), 'Press Enter to Continue') pause P=PHI'*p; PP1=PHI'*pp1; %uzero=[0 0 0]';% initial displacements U=[uzero]; Q=zeros(ndof,201); udotzero=[0 0 0]';% initial velocities qzero=PHI'*m*uzero; qdotzero=PHI'*m*udotzero; for i=1:ndof wn=omegan(1,i); wD=omegad(1,i); z=zeta(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 Si; Q(i,1)=qzero(i,1); x=[qzero(i,1);qdotzero(i,1)]; for k1=2:201 x=phi*x+gamma*[P(i,k1-1);PP1(i,k1-1)]; Q(i,k1)=x(1,1); end end U=PHI*Q; figure(3);clf plot(t,U) xlabel('Time t,sec.') ylabel('Displacement, ft.') legend('u_1(t)','u_2(t)','u_3(t)') Umax=max(max(abs(U(1:3,:)))); text(.4,.8*Umax, 'Press Enter to Continue') pause % Draw frame and prepare to animate U=U*5/Umax;%scale responses for animation 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; 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]; basex=[-5 35 29 22 14 10 4 -5]; basey=[0 0 -4 -2 -3 -1 -5 0]; arrowx=[0 6 4.5 4.5 6];arrowy=[0 0 .7 -.7 0]; figure(4);clf axis square axis off axis([-10 40 -10 40]); patch(basex,basey,'r');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);hold on D2=plot(damp2x,damp2y);hold on D3=plot(damp3x,damp3y);hold on D4=plot(damp4x,damp4y);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(1.4*arrowx-8.4*ones(1,5),arrowy+12*ones(1,5),'k'); Arrow5=plot(1.4*arrowx-8.4*ones(1,5),arrowy+24*ones(1,5),'k'); Arrow6=plot(1.4*arrowx-8.4*ones(1,5),arrowy+36*ones(1,5),'k'); Tx1=text(15,13.5,'m_1');Tx2=text(15,25.5,'m_2');Tx3=text(15,37.5,'m_3'); Tx4=text(31,4,'k_1/2');Tx5=text(31,16,'k_2/2');Tx6=text(31,28,'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'); Tx14=text(-9,15,'p_1(t)');Tx15=text(-9,27,'p_2(t)');Tx14=text(-9,39,'p_3(t)'); Tx7=text(-5,-7,'Press Enter to Set Initial Displacements'); pause clf axis square axis off axis([-10 40 -10 40]); patch(basex,basey,'r') hold on i=1; e=[0 0 .15 .5 .85 1 1]; deltax1=U(1,1); deltax2=U(2,1)-U(1,1); deltax3=U(3,1)-U(2,1); patch(basex,basey,'r');hold on B1=patch(box1x+ones(1,5)*U(1,1),box1y,'r','EraseMode','xor');hold on B2=patch(box2x+ones(1,5)*U(2,1),box2y,'r','EraseMode','xor');hold on B3=patch(box3x+ones(1,5)*U(3,1),box3y,'r','EraseMode','xor');hold on C1=plot(col1x+deltax1*e,col1y,'k','EraseMode','xor');hold on C2=plot(col2x+deltax1*e,col2y,'k','EraseMode','xor');hold on C3=plot(col3x+deltax2*e+U(1,i)*ones(1,7),col3y,'k','EraseMode','xor');hold on C4=plot(col4x+deltax2*e+U(1,i)*ones(1,7),col4y,'k','EraseMode','xor');hold on C5=plot(col5x+deltax3*e+U(2,i)*ones(1,7),col5y,'k','EraseMode','xor');hold on C6=plot(col6x+deltax3*e+U(2,i)*ones(1,7),col6y,'k','EraseMode','xor');hold on D1=plot(damp1x,damp1y,'k');hold on D2=plot(damp2x,damp2y,'k');hold on D3=plot(damp3x+U(1,1)*[1 1],damp3y,'k');hold on D4=plot(damp4x+U(1,1*[1 1 1]),damp4y,'k');hold on D21=plot(damp1x+U(1,1)*[1 1 1],damp1y+[ 1 1 1 ]*12,'k');hold on D22=plot(damp2x+U(1,1)*[1 1 1 1],damp2y+[ 1 1 1 1]*12,'k');hold on D23=plot(damp3x+U(2,1)*[1 1],damp3y+[ 1 1 ]*12,'k');hold on D24=plot(damp4x+U(2,1)*[1 1 1],damp4y+[ 1 1 1 ]*12,'k');hold on D31=plot(damp1x+U(2,1)*[1 1 1],damp1y+[ 1 1 1]*24,'k');hold on D32=plot(damp2x+U(2,1)*[1 1 1 1],damp2y+[ 1 1 1 1]*24,'k');hold on D33=plot(damp3x+U(3,1)*[1 1],damp3y+[ 1 1 ]*24,'k');hold on D34=plot(damp4x+U(3,1)*[1 1 1],damp4y+[ 1 1 1]*24,'k');hold on Tx7=text(-5 ,-7,['Press Enter to Animate, Max. Response = ' num2str(Umax) ' ft.']); Tx8=text(-5,-9,['t = ' num2str(t(1,i)) ' sec.']); pause for i=1:201 deltax1=U(1,i); deltax2=U(2,i)-U(1,i); deltax3=U(3,i)-U(2,i); set(B1,'Xdata',box1x+U(1,i)*ones(1,5));hold on set(B2,'Xdata',box2x+U(2,i)*ones(1,5));hold on set(B3,'Xdata',box3x+U(3,i)*ones(1,5));hold on set(C1,'Xdata',col1x+deltax1*e);hold on set(C2,'Xdata',col2x+deltax1*e);hold on set(C3,'Xdata',col3x+deltax2*e+U(1,i)*ones(1,7));hold on set(C4,'Xdata',col4x+deltax2*e+U(1,i)*ones(1,7));hold on set(C5,'Xdata',col5x+deltax3*e+U(2,i)*ones(1,7));hold on set(C6,'Xdata',col6x+deltax3*e+U(2,i)*ones(1,7));hold on set(D1,'Xdata',damp1x) set(D2,'Xdata',damp2x) set(D3,'Xdata',damp3x+U(1,i)*ones(1,2)); set(D4,'Xdata',damp4x+U(1,i)*ones(1,3)); set(D21,'Xdata',damp1x+U(2,i)*ones(1,3)); set(D22,'Xdata',damp2x+U(2,i)*ones(1,4)) set(D23,'Xdata',damp3x+U(2,i)*ones(1,2)); set(D24,'Xdata',damp4x+U(2,i)*ones(1,3)); set(D31,'Xdata',damp1x+U(2,i)*ones(1,3)); set(D32,'Xdata',damp2x+U(2,i)*ones(1,4)); set(D33,'Xdata',damp3x+U(3,i)*ones(1,2)); set(D34,'Xdata',damp4x+U(3,i)*ones(1,3)); set(Tx8,'String',['t = ' num2str(t(1,i)) ' sec.']); pause(.04) end