% torsionbar2.m Last modified 8/18/2008 % An elastic torsion bar fixed at the left end and driven by a suddenly % applied torque T at the right end and no initial twist or velocity. % Solution via 40 terms of the Fourier series solution to the 1-d % 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)=((-1)^(i+1))/((2*i-1)^2); end for i=2:40 omegaoveromega1=[omegaoveromega1 2*i-1]; end for i=1:nterm term=ones(1,401)-cos(omegaoveromega1(1,i)*omega1t); cs=[cs;term]; 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;sum]; end framedata=(8/(pi^2))*framedata; figure(1);clf; xlabel('Dimensionless Distance, x/L') hold on ylabel('Torsional Displacement, JG \theta(x,t)/TL') S=plot(xoverL,framedata(1,:),'k','EraseMode','xor'); text(-.1,2,'Press Enter to Continue') axis([-.2 1.2 -.2 2.2]) hold on box on pause for k=11:10:201 set(S,'Ydata',framedata(k,:)); text(-.1,2,'Press Enter to Continue') pause end hold off figure(2);clf; axis([0 4*pi -.2 2.2]) hold on text(8 ,2,'x/L=') box on for j=1:10:101 plot(omega1t,framedata(:,j)) hold on text(9.4,framedata(300,j),num2str(xoverL(1,j))) hold on end xlabel('Dimensionless Time, \omega_1t') ylabel('Torsional Displacement, JG \theta(x,t)/TL') text(2,2.05,'Press Enter to Continue') hold off pause %now do animation figure(3);clf; axis([-.2 1.2 -.2 2.2]) box on xlabel('Dimensionless Distance, x/L') hold on ylabel('Torsional Displacement, JG \theta(x,t)/TL') hold on L=plot(xoverL, zeros(1,101),'k','EraseMode','xor'); hold off texthandl1=text(-.1,2,'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; framedata=(3*pi/8)*framedata; 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 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 xarrow=[2*fliplr(xfront) 0 .03 0 0 .03]+ones(1,36); yarrow=[2*fliplr(yfront) .4 .4 .38 .42 .4]; plot(xarrow,yarrow); text(1,.58,'T'); hold on 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