% torsionbar.m Last modified 8/18/2008 % An elastic torsion bar fixed at the left end and free at the right end % with initial deformation proportional to the distance from the fixed end % and no initial velocity. Rotational displacement at the free end is pi/2 % radians. Solution via 40 terms ofthe Fourier series solution to the % wave equation. % 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); xoverL=linspace(0,1,101); omega1t=linspace(0,4*pi,401); cs=[]; nterm=40;%40 modes omegaoveromega1=[1]; for i=1:nterm fac=((2*i)-1)*pi/2; sine(i,:)=sin(fac*xoverL); B(1,i)=2*(((1/fac)^2)*sin(fac)); end for i=2:40 omegaoveromega1=[omegaoveromega1 2*i-1]; end for i=1:nterm cs1=cos(omegaoveromega1(1,i)*omega1t); cs=[cs;cs1]; end framedata=[]; for k=1:401 sum=zeros(1,101); for i=1:nterm sum=sum+B(1,i)*cs(i,k)*sine(i,:); end framedata=[framedata;(pi/2)*sum]; end figure(1);clf; xlabel('Dimensionless Distance, x/L') hold on ylabel('Torsional Displacement, \theta(x,t)') S=plot(xoverL,framedata(1,:),'k','EraseMode','xor'); text(-.1,1.1,'Press Enter to Continue') axis([-.2 1.2 -2 2]) box on pause for k=11:10:101 set(S,'Ydata',framedata(k,:)); axis([-.2 1.2 -2 2]) text(-.1,1.1,'Press Enter to Continue') pause end hold off figure(2);clf; axis([0 4*pi -2 2]) hold on text(5.2 ,1.57,'x/L=') box on for j=1:10:101 plot(omega1t,framedata(:,j)) hold on text(6,framedata(1,j),num2str(xoverL(1,j))) hold on end xlabel('Dimensionless Time, \omega_1t') ylabel('Torsional Displacement, \theta(x,t)') text(4,-1,'Press Enter to Continue') pause hold off %now do animation figure(3);clf; axis([-.2 1.2 -2 2]) box on xlabel('Dimensionless Distance, x/L') hold on ylabel('Torsional Displacement, \theta(x,t)') hold on L=plot(xoverL, zeros(1,101),'k','EraseMode','xor'); hold off texthandl1=text(-.1,1.1,'Press Enter to Set Initial Condition'); 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 hold off figure(4);clf; R=.2;R1=.16; phi=(pi/15); xcir=[];ycir=[]; for i=1:61 gamma=2*pi*i/60; xcir=[xcir R*cos(gamma)]; ycir=[ycir R*sin(gamma)]; end xellipse=xcir*sin(phi); yellipse=ycir; xfront=xellipse(1,16:46); xback=[xellipse(1,46:61) xellipse(1,2:16)]; yfront=yellipse(1,16:46); yback =[yellipse(1,46:61) yellipse(1,2:16)]; axis([-.2 1.2 -1 1]); ypatch=[-.5 [-.2 fliplr(yfront) .2] .5 .3 -.1 -.35 -.5]; xpatch=[0 [0 fliplr(xfront) 0] 0 -.15 -.08 -.17 0]; patch(xpatch, ypatch,'r'); hold on x1=linspace(0.1, 1,10); for i=1:9 plot(xfront+x1(1,i),yfront,'-'); hold on plot(xback+x1(1,i),yback,'--') hold on end plot(xfront+1,yfront,'-'); hold on plot(xback+1,yback,'-') hold on xbar1=[0 1];xbar2=[1 0]; ybar1=[.2 .2];ybar2=[ -.2 -.2]; hold on box on plot(xbar1,ybar1);plot(xbar2,ybar2) R1sp=R1*sin(phi);R1cp=R1*cos(phi); for i=1:10 xtic=x1(1,i);ytic=R1; tic(1,i)=plot([xtic xtic],[ytic,-ytic],'EraseMode','xor','LineWidth',[2]); end %pause xtic1=zeros(401,10); ytic1=zeros(401,10); for i=1:10; for j=1:401; xtic1(j,i)=R1sp*sin(framedata(j,10*i+1)); ytic1(j,i)=R1*cos(framedata(j,10*i+1)); end end t1=text(.3,.8,'Press Enter to Set Initial Condition'); xlabel('Distance Along the Bar, x/L') pause for i=1:10 set(tic(1,i), 'Xdata', x1(1,i)*ones(1,2)+[xtic1(1,i) -xtic1(1,i)],'Ydata',[ytic1(1,i),-ytic1(1,i)]); end set(t1,'String','Press Enter to Animate '); pause set(t1,'String',' '); for k=2:401 for i=1:10 set(tic(1,i), 'Xdata', x1(1,i)*ones(1,2)+[xtic1(k,i) -xtic1(k,i)],'Ydata',[ytic1(k,i),-ytic1(k,i)]); end pause(.04) end set(t1,'String', 'Press Enter to Continue') pause hold off