% clampedclampedbeam.m Last modified 7/28/2008 % Animate the free vibration of a uniform Bernoulli-Euler Clamped-Clamped % beam with initial conditions composed of quadratic and % cubic factors in x/L i.e. y(x,0)=2(x/L)^2-(x/L)^3-4(x/L)^4+3(x/L)^5 % The initial velocity is zero. % 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); betaL=[4.7300408 7.8532046 10.9956078 14.1371655 17.2787596]; betaLsq=betaL.*betaL; alpha=[.9825022158 1.000777311 .9999664501 1.000001450 .9999999373]; xoverL=linspace(0,1,101); omega1t=linspace(0,4*pi,401); sn=[]; % initial condition function y0=2*(xoverL.^2)-1*(xoverL.^3)-4*(xoverL.^4)+3*(xoverL.^5); y0max=max(abs(y0)); nterm=5; b=zeros(1,5); for i=1:nterm z=betaL(1,i)*xoverL; phi(i,:)=cosh(z)-cos(z)-alpha(1,i)*(sinh(z)-sin(z)); cs(i,:)=cos(betaLsq(1,i)*omega1t/betaLsq(1,1)); b(1,i)=.01*sum(phi(i,:).*y0); end xoverL=[-.06 xoverL 1.06]; framedata=[]; for k=1:401 sum=zeros(1,101); for i=1:nterm sum=sum+b(1,i)*cs(i,k)*phi(i,:); end framedata=[framedata;sum]; end lim=1; framedata=[zeros(401,1) framedata zeros(401,1)]/y0max; figure(1);clf; h1=axes('YDir','reverse'); axis([-.2 1.2 -1.2*lim 1.2*lim]) xpatch1=[0 0 -.1 -.08 -.09 0]; ypatch1=[.2 -.2 -.15 -.05 .05 .2]*lim/.4; xpatch2=ones(1,6)-xpatch1; xlabel('Dimensionless Distance, x/L') ylabel('Dimensionless Displacement, y(x,t)/y_0_m_a_x') patch(xpatch1,ypatch1,'r'); patch(xpatch2,ypatch1,'r'); hold on S=plot(xoverL,framedata(1,:),'k','LineWidth',[3.5]); hold on box on text(0,1.05*lim,'y(x,0)=2(x/L)^2-(x/L)^3-4(x/L)^4+3(x/L)^5'); texthandl=text(.1 ,-1.1*lim,'Press Enter to Animate One Frame at a Time'); pause for k=11:10:201 plot(xoverL,framedata(k,:),'k','LineWidth',[3.5]); set(texthandl,'String','Press Enter to Continue'); pause end set(gca,'Box','on') figure(2);clf; h2=axes('YDir','reverse'); hold on axis([0 4*pi -1.2*lim 1.2*lim]); plot(omega1t,framedata(:,2)); hold on for j=12:10:102 plot(omega1t,framedata(:,j)); hold on end for j=12:10:92 plot([pi 5.5],[framedata(101,j) -j*.01]); text(5.55,-j*.01,['x/L=' num2str(xoverL(1,j))]); end xlabel('Dimensionless Time, \omega_1t') ylabel('Dimensionless Displacement, y(x,t)') text(1,-1.1*lim,'Press Enter to Continue') hold off pause %now do animation figure(3);clf; h3=axes('YDir','reverse'); axis([-.2 1.2 -1.2*lim 1.2*lim]) xlabel('Dimensionless Distance, x/L') hold on ylabel('Dimensionless Displacement, y(x,t).y_0_m_a_x') hold on patch(xpatch1,ypatch1,'r') patch(xpatch2,ypatch1,'r') L=plot(xoverL, zeros(1,103),'k','EraseMode','xor','LineWidth',[3.5]); plot([xoverL(1,1) xoverL(1,2)],[0 0],'k','Linewidth',[3.5]); hold on box on texthandl1=text(.1,-1.1*lim,'Press Enter to Set Initial Condition'); text(0,1.05*lim,'y(x,0)=2(x/L)^2-(x/L)^3-4(x/L)^4+3(x/L)^5'); pause set(L,'Ydata',framedata(1,:)) set(texthandl1,'String','Press Enter to Animate'); pause set(texthandl1,'String',' '); for i=2:401 set(L,'Ydata',framedata(i,:)) pause(.04) end hold on set(texthandl1,'String','Press Enter to Continue'); pause figure(4);clf; h4=axes('ZDir','reverse'); hold on axis([0 1 0 4*pi -1.05*lim 1.05*lim]) hold on [X,Y]=meshgrid(xoverL(1,2:102),omega1t); mesh(X,Y,framedata(:,2:102)) xlabel('Distance, x/L','Rotation',-31) ylabel('Dimensionless Time, \omega_1t','Rotation',11) zlabel('Dimensionless Temperature, y(x,t)/y_0_m_a_x') view(60,30) grid on text(0,4*pi,-.9,'Press Enter to End','Rotation',-31) pause