% canttipforceanimation.m Last modified 7/29/2008 % Cantilever beam animation via the steady state solution to the % Bernoulli-Euler beam equation. Beam is forced at the tip by a force % of amplitude Fo.Input is the ratio of driving freq. to the first % natural freq. % 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); fratio=-1 while fratio<=0|fratio==1 string1=' The ratio of excitation freq. to first natural freq. = '; fratio=input(string1); if fratio==1 disp(' ') disp(' This is the resonant case choose another freq. ratio'); disp(' ') end end xoverL=linspace(0,1,201); betaL=sqrt(fratio)*1.8751041; betax=betaL*xoverL; cL=cos(betaL); sL=sin(betaL); chL=cosh(betaL); shL=sinh(betaL); c1=cos(betax); s1=sin(betax); ch1=cosh(betax); sh1=sinh(betax); delta=-2*(1+chL*cL); A=(sL+shL)*(c1-ch1)-(cL+chL)*(s1-sh1); amplitude=(3/((betaL^3)*delta))*A; lim=max(abs(amplitude)); figure(1);clf; plot(xoverL,abs(amplitude)) hold on plot([0 1],[0 0]); axis([0 1 -.2*lim 1.2*lim]); text(.8,1.1*lim,['f/f_1 = ',num2str(fratio)]) xlabel('Dimensionless distance from the fixed end, x/L') ylabel('Dimensionless Amplitude, y(x)/y_s_t') text(.1,.6*lim,'y_s_t = F_0L^3/3EI') text(.1,1.1*lim,'Press Enter to Continue') amplitude=[0 amplitude]; hold off xoverL=[-.07 xoverL]; pause figure(2);clf; set(gca,'Box','on') axis([-.3 1.3 -2*lim 2*lim]); xlabel('Distance, x/L') ylabel('Displacement, y(x,t)/y_s_t') hold on texthandle2=text(-.1,-1.6*lim,['f/f_1 = ',num2str(fratio)]); text(.6,-1.6*lim,'y_s_t = F_0L^3/3EI') ys=amplitude*0; plot(xoverL,ys,'k','LineWidth',[1.5]); hold on ydata=[ys]; texthandle2=text(0,1.5*lim,'Press Enter to Animate one Frame at a Time'); xlabel('Distance , x/L'); pause t=linspace(0,8*pi,257); sine=sin(t); for k=1:2:65 ys=amplitude*sine(1,k); plot(xoverL, ys,'k','LineWidth',[2.5]); hold on set(texthandle2,'String', 'Press Enter to Animate One Frame at a Time'); pause end hold off set(texthandle2,'String','Press Enter to Continue'); pause figure(3):clf; xpatch=[0 0 -.15 -.08 -.12 -.07 -.1 0]; ypatch=3*lim*[-.2 .2 .14 .05 -.01 -.14 -.18 -.2]; axis([-.3 1.3 -2*lim 2*lim]); hold on set(gca,'Box','on'); hold on P=patch(xpatch,ypatch,'r'); hold on xlabel('Distance, x/L'); ylabel('Displacement, y(x,t)/y_s_t'); hold on x1=[-.3 1.3];y1=[0 0]; plot(x1,y1); hold on texthandle2=text(-.1,-1.6*lim,['f/f_1 = ',num2str(fratio)]); text(.6,-1.6*lim,'y_s_t = F_0L^3/3EI') ys=amplitude*0; L=plot(xoverL,ys,'k','EraseMode','xor','LineWidth',[3.5]); ydata=[ys(1,2:202)]; texthandl=text(0,1.5*lim,'Press Enter to Animate'); xlabel('Distance , x/L'); ylabel('Displacement, y(x,t)/y_s_t'); pause set(texthandl,'String',' '); t=linspace(0,8*pi,257); sine=sin(t); for k=2:257 ys=amplitude*sine(1,k);; set(L,'Ydata',ys); pause(.05) ydata=[ydata ;ys(1,2:202)]; end hold off set(texthandl,'String','Press Enter to Continue'); pause figure(4);clf; [X,T]=meshgrid(xoverL(1,2:202),t); mesh(X,T,ydata) axis([0 1 0 30 -1.2*lim 1.2*lim]); text(1,5,-lim,'Dimensionless Time, \omegat','Rotation',-10) text(1,30,-lim,'Dimensionless Distance, x/L','Rotation',32) zlabel('Displacement, y(x,t)/y_s_t') text(0,4, lim,'Press Enter to Continue','Rotation',-10) view(120,30)