clear all clc disp(' This script evaluates the forced response of a reinforced-concrete') disp(' chimney of Example 15.1 in Dynamics of Structuresby A.K. Chopra.') disp(' Problem is solved by the Chopra interpolation of excitation method.') disp('The parametric matricies are in kips,feet and sec. ') disp(' ') disp(' Press Enter to Continue') pause ndof=5; % number of degrees of freedom disp('The mass matrix is') m=[1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 .5]*208.6 disp('The stiffness matrix is') k=[18.83 -11.9 4.773 -1.193 .1989;-11.9 14.65 -10.71 4.177 -.6961]; k=[k;4.773 -10.71 14.06 -9.514 2.586;-1.193 4.177 -9.514 9.878 -3.646]; k=[k;.1989 -.6961 2.586 -3.646 1.608];% k=k*3.165e4 alpha=[0 0]; alpha(1,1)=0;%%% c=alpha(1,1)*m+alpha(1,2)*k D=inv(m)*k; [PHI1,lambda]=eig(D); omegan1=[]; for i=1:ndof omegan1=[omegan1 sqrt(lambda(i,i))]; end [omegan1, index]=sort((omegan1')); disp('The radian undamped natural frequencies are') omegan=omegan1' Phi=zeros(ndof,ndof); % reorder eigenvectors for i=1:ndof PHI(:,i)=PHI1(:,index(i,1)); 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 K=PHI'*k*PHI; disp('The Hertzian undamped natural frequencies are') fn=omegan/(2*pi) disp('The undamped natural periods are') Tn=ones(1,ndof)./fn% natural periods in sec. C=PHI'*c*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 24 48 72 96 120]; plot(100*[zeros(5,1) PHI'],z1) axis([-30 30 -20 130]); legend('i=1','i=2','i=3','i=4','i=5') patch(basex1,basey1,'r');hold on text(-25, 110,'Mode Shapes') text(-25,-18,'Press Enter to Continue'); pause dt=.02; minv=inv(m); t=[0:dt:600*dt]; % step forcing function at chimney tip p1=1000;% force applied to 5th floor, kips p=[0;0;0;0;p1]*ones(1,601); figure(2);clf plot([-.5 0 t],[0 0 p1*ones(1,601)]) axis([-.5 t(1,201) -100 1100]); xlabel('Time t,sec.') ylabel('Chimney Top Force, f(t), kips') text(1,.8*max(max(p)), 'Press Enter to Continue') pause P=PHI'*p; PP1=[P(:,2:601) P(:,601)]; Q=zeros(ndof,201); uzero=[0 0 0 0 0]'; udotzero=[0 0 0 0 0]';% initial velocities U=[uzero]; 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 Q(i,1)=qzero(i,1); x=[qzero(i,1);qdotzero(i,1)]; for k1=2:601 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(1,:),t,U(2,:),t,U(3,:),t,U(4,:),t,U(5,:)) xlabel('Time t,sec.') ylabel('Displacement, ft.') legend('u_1(t)','u_2(t)','u_3(t)','u_4(t)','u_5(t)') Umax=max(max(abs(U(1:5,:)))); text(.4,.8*Umax, 'Press Enter to Continue') pause % Draw frame and prepare to animate U=U*5/Umax;%scale responses for animation xcoord=[-10 -9 -8 -7 -6 -5]; xcoord=[xcoord fliplr(-xcoord) -10]; ycoord=[0 24 48 72 96 120]; ycoord=[ycoord fliplr(ycoord) 0]; basex=[-20 20 18 14 7 -3 -8 -20]; basey=[0 0 -11 -7 -15 -9 -12 0]; figure(4);clf axis square axis off axis([-30 30 -20 130]); patch(basex,basey,'r');hold on Chimney=patch(xcoord,ycoord,'y'); txt1=text(-15,127,'Press Enter to Continue'); pause i=1; set(Chimney,'Xdata',xcoord+[0 U(:,1)' fliplr(U(:,1)') 0 0]); patch(basex,basey,'r');hold on txt2=text(-25,-15,['t = ' num2str(t(1,1)) ' sec.']); txt3=text(15,-15,['Max. Disp. = ' num2str(Umax) 'ft.']); set(txt1,'String','Press Enter to Amimate '); pause for i=1:601 set(Chimney,'Xdata',xcoord+[0 U(:,i)' fliplr(U(:,i)') 0 0]); set(txt2,'String',['t = ' num2str(t(1,i)) ' sec.']); pause(.01) end