function cantvib3 % Last modified 7/15/2008 % This function illustrates the free vibration of a cantilever beam with % an initial deflection curve that is both quadratic and cubic in the % spatial variable. This amounts to both an initial transverse load and % tip bending moment released at t=0.The solution is accomplished by % discretizing the fourth derivative of the Bernoulli-Euler beam equation % a sixteen nodes. This script requires another function appended to the % end of the main script called beamrhs which calculates the righthand % side of the differential equation for ODE45. % This m-file was written at the University of Wyoming in % the Electrical and Computer Engineering Department and % is to be distributed without cost. clc set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); n=16;% number of nodes M=[7 -4 1 zeros(1,n-3)]; M=[M;-4 6 -4 1 zeros(1,n-4)]; for i=3:n-2 M=[M;zeros(1,i-3) 1 -4 6 -4 1 zeros(1,n-i-2)]; end M=[M;zeros(1,n-4) 1 -4 5 -2]; M=[M;zeros(1,n-3) 2 -4 2]; M=((n^4)/12.3624)*M; global A; xoverl=[0:1/n:1]; xoverl1=xoverl(:,2:n+1); w0=(.667*(xoverl1.^2)+.333*(xoverl1.^3)); lim1=max(w0); x0=[w0 zeros(1,n)]; A=[zeros(n,n) eye(n,n);-M zeros(n,n)]; tspan=linspace(0,4*pi,601); omega1t=linspace(0,4*pi,601); [t,x]=ode45(@beamrhs,tspan,x0); w=x(:,1:n); sz=size(t); w=[zeros(sz(1,1),1) w]; figure(1);clf;% Plot y(x,t) vs x/L for various times h1=axes('YDir','reverse'); axis([-.2 1.2 -1.2*lim1 1.2*lim1]) hold on text(0,-1,'y(x,0)=y_0[0.667(x/L)^2+0.333(x/L)^3]') hold on xp=[ 0 0 -.15 -.13 -.16 -.14 -.15 0]; yp=[-.37 .37 .2 .05 -.15 -.2 -.3 -.37]; patch(xp,yp,'r'); hold on plot([-.1 0],[0 0],'k','LineWidth',[3.5]) hold on plot([0 1.1],[0 0]) hold on text(1.1,-.05,'x/L') hold on for i=1:10:150 plot(xoverl,w(i,:),'k','LineWidth',[3.5]); hold on end xlabel('Dimensionless Distance, x/L') ylabel('Dimensionless Displacement, y(x,t)/y_0') box on hold off text(.2,1,'Press Enter to Continue'); pause figure(2);clf;%plot y(xi,t) vs omega1t for the nodes xi h1=axes('YDir','reverse'); axis([0 4*pi -1.2*lim1 1.2*lim1]); hold on for i=2:2:n plot(omega1t,x(:,i)) hold on plot([2*pi,8],[x(301,i) .1+(i/n)],'k') hold on text(8.1,.1+(i/n),['x/L=' num2str(xoverl(1,i+1))]); end xlabel('Dimensionless Time, \omega_1t') ylabel('Displacement, y(x,t)/y_0') text(1,1,'Press Enter to Continue'); %box on hold off pause [T,X]=meshgrid(omega1t,xoverl); figure(3);clf; axis([0 4*pi 0 1 -1.2*lim1 1.2*lim1]); hold on mesh(T,X,w') colormap([0 0 0]); view(60,30) xlabel('Dimensionless Time, \omega_1t','rotation',-31) ylabel('Dimensionless Location, x/L','rotation',12) zlabel('Displacement, y(x,t)/y_0') grid on hold off pause figure(4);clf;%Animate y(x,omega1t) h1=axes('YDir','reverse'); axis([-.2 1.2 -1.2*lim1 1.2*lim1]) hold on xp=[ 0 0 -.15 -.13 -.16 -.14 -.15 0]; yp=[-.37 .37 .2 .05 -.15 -.2 -.3 -.37]; patch(xp,yp,'r'); hold on plot([-.1 0],[0 0],'k','LineWidth',[3.5]) hold on plot([0 1.1],[0 0]) L=plot(xoverl,zeros(1,n+1),'k','EraseMode','xor','LineWidth',[3.5]); hold on text(1.1,-.05,'x/L') hold on box on text(0,-1,'y(x,0)=y_0[0.667(x/L)^2+0.333(x/L)^3]') thandl=text(0.2,1,'Press Enter to Set Initial Condition'); xlabel('Dimensionless Distance, x/L') ylabel('Dimensionless Displacement, y(x,t)/y_0') pause for i=1:sz(1,1) set(L,'YData',w(i,:)); if i==1 set(thandl,'String','Press Enter to Animate'); pause end drawnow; pause(.01) end hold off set(thandl,'String','Press Enter to Continue'); pause function xdot=beamrhs(t,x) global A; xdot=A*x;