% cantbeamimpulse.m Last modified 7/28/2008 %Solves and animates the problem of transient motion of a uniform cantilever %forced by an impulse of intensity Io at location x = a. %This m-file was written at the University of Wyoming in the Electrical %and Computer Engineering Department and is to be distributed without %cost. clear all set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); aoverL=-1; while ((aoverL<=0)|(aoverL>1)) aoverL=input('Input the location of the impulse, a/L = '); end phiofL=[2 -2 2 -2 2]; boverL=1-aoverL; pt=(100*aoverL)+1; betaL=[1.8751041 4.69409113 7.85475743 10.99554074 14.13716839]; gamma=(betaL/betaL(1,1)).^2; betaL4=betaL.^4; alpha=[.7340995 1.0184664 .9992245 1.000033553 .9999985501]; betan=betaL4/betaL4(1,1); x1=linspace(0,1,101); phiofx=[]; omega1t=linspace(0,8*pi,1001); for i=1:5 arg1=betaL(1,i)*x1; phiofx=[phiofx;cosh(arg1)-cos(arg1)-alpha(1,i)*(sinh(arg1)-sin(arg1))]; end phiofa=phiofx(:,pt)'; Y=[]; for t=1:1001 sum=zeros(1,101); for i=1:5 sum=sum+(phiofa(1,i)*phiofx(i,:)/gamma(1,i))*sin(gamma(1,i)*omega1t(1,t)); end Y=[Y;sum]; end Y=[zeros(1,101);Y]; Y=[zeros(1002,1) Y]; figure(1);clf;% Plot y(x,t) at x/L=0.2,0.4,0.6,0.8,1 h3=axes('YDir','reverse'); hold on axis([0 8*pi -5 5]); hold on for j=22:20:102 plot(omega1t,Y(2:1002,j)) hold on end for j=22:20:102 plot([omega1t(1,431) 13],[Y(432,j) -(2+0.5*((j-2)/20))],'k'); hold on text(13.2,-(2+0.5*((j-2)/20)),['x/L = ' num2str(x1(1,j-1))]); hold on end box on xlabel('Dimensionless Time, \omega_1t') ylabel('Beam Displacement, y(x,t)\muL\omega_1/I_0') text(3,4,['a/L = ' num2str(aoverL)]) hold off text(2,-4,'Press Enter to Continue'); pause figure(2);clf; h3=axes('YDir','reverse'); hold on axis([-.2 1.3 -5 5]); hold on box on xpatch1=[0 0 -.1 -.08 -.09 0]; ypatch1=[.2 -.2 -.15 -.05 .05 .2]*6; xlabel('Dimensionless Distance, x/L') ylabel('Beam Displacement, y(x,t)\muL\omega_1/I_0') patch(xpatch1,ypatch1,'r') hold on plot([-.05 x1], Y(1,:),'k','LineWidth',[3.5]); hold on for i=25:24:253 plot([-.05 x1], Y(i,:),'k','LineWidth',[3.5]); hold on text(0,-4,'Press Enter to Continue'); pause end text(0,4,'y(x,0)=0'); text(0,4.5,['Impulse Location a/L = ' num2str(aoverL)]) hold off text(0,-4,'Press Enter to Continue'); pause figure(3);clf;% plot motion of point where load is applied h3=axes('YDir','reverse'); hold on axis([0 8*pi -5 5]); hold on plot(omega1t,Y(2:1002,pt+1)); xlabel('Dimensionless Time, \omega_1t') ylabel('Forced Point Displacement, y(a,t)\muL\omega_1/I_0') text(2,-4,'Press Enter to Continue') text(2,4,['a/L = ' num2str(aoverL)]) box on pause hold off; figure(4);clf;% Animate the beam h3=axes('YDir','reverse'); hold on axis([-.2 1.3 -5 5]) hold on xpow=[1.1 1.02 1.05 .95 .88 .92 .9 1.01 1.04 1.03 1.1]; ypow=-10*[.1 .14 .17 .13 .15 .07 .04 .06 .05 .08 .1]; xpatch1=[0 0 -.1 -.08 -.09 0]; ypatch1=[.2 -.2 -.15 -.05 .05 .2]*6; xlabel('Dimensionless Distance, x/L') ylabel('Beam Displacement, y(x,t)\muL\omega_1/I_0') patch(xpatch1,ypatch1,'r') L=plot([-.05 x1], Y(1,:),'k','EraseMode','xor','LineWidth',[3.5]); box on texthandl1=text(0,-4,'Press Enter to Set Initial Condition'); text(0,4,'y(x,0)=0'); text(0,4.5,['Impulse Location a/L = ' num2str(aoverL)]) pause xpow=xpow+(-1+x1(1,pt))*ones(1,11); POW=patch(xpow,ypow,'m'); set(texthandl1,'String','Press Enter to Animate'); pause set(POW,'Visible','off') set(texthandl1,'String',' '); for i=3:1:1001 set(L,'Ydata',Y(i,:)) pause(.005) end set(texthandl1,'String','Press Enter to Continue'); pause hold off;