% conductioncyl.m Last Modified 7/18/2008 % Conduction of heat in an infinite cylinder of radius R % with exterior face (r=R) at T0 and and initial interior at zero % temperature. The diffusivity is kappa=k/(rho*c) % 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); % eigenvalues--zeros of the J0 Bessel function lambdaR=[2.404826 5.520078 8.653728 11.791534 14.9309177 18.0761064]; lambdaR=[lambdaR 21.211637 24.352472 27.493479 30.634606]; roverR=linspace(0,1,101); nterm=10; for i=1:nterm J0(i,:)=besselj(0,lambdaR(1,i)*roverR); J1(1,i)=besselj(1,lambdaR(1,i)); end tau=linspace(0,.4,201); T1=(ones(1,101)); framedata=[]; for k=1:201 sum=zeros(1,101); for i=1:nterm sum=sum+(2/(lambdaR(1,i)*J1(1,i)))*exp(-((lambdaR(1,i)^2)*tau(1,k)))*J0(i,:); end T=T1-sum; framedata=[framedata;T]; end framedata(1,:)=zeros(1,101); framedata(1,101)=1; figure(1);clf;% Plot T(x,t) vs x for various t plot([0 0],[-.2,1.2]); hold on plot([1 1],[-.2,1.2]); hold on plot(roverR,framedata(1,:),'k') hold on xlabel('Dimensionless Radius, r/R') hold on ylabel('Dimensionless Temperature, T(r,t)/T_0') for k=21:20:201 plot(roverR,framedata(k,:),'k') hold on axis([-.2 1.2 -.2 1.2]) end text(.7,.05,['\kappat/R^2 = ' num2str(0)]) text(.6,.2,num2str(0.04)) text(.5,.3,num2str(0.08)) text(0.2,.9,num2str(0.4)) set(gca,'Box','on') hold off text(.1,1.13,'Press Enter to Continue') pause figure(2);clf;%Plot T(x,t) vs t for various x axis([0 .3 -.2 1.2]); hold on plot(tau,ones(1,201)) hold on for j=11:10:101 plot(tau,framedata(:,j)) hold on end box on text(.02,.95,['r/R = ' num2str(1)]) text(.04,.87,['r/R = ' num2str(0.9)]) text(.16,.35,['r/R = ' num2str(0)]) xlabel('Dimensionless Time, \kappat/R^2') ylabel('Dimensionless Temperature, T(r,t)/T_0') hold off text(.02,1.1,'Press Enter to Continue') pause figure(3);clf;%now do animation axis([-.2 1.2 -.2 1.2]) hold on plot([0 0],[-.2,1.2]); plot([1 1],[-.2,1.2]); box on xlabel('Dimensionless Radius, r/R') hold on ylabel('Dimensionless Temperature, T(r,t)/T_0') hold on L=plot(roverR, framedata(1,:),'k','EraseMode','xor'); hold on texthandl=text(.4,1.1,'Press Enter to Animate'); pause set(texthandl,'String',' '); for i=2:201 set(L,'Ydata',framedata(i,:)) pause(.05) end hold on plot(roverR,framedata(1,:),'k') hold off set(texthandl,'String','Press Enter to Continue'); pause figure(4);clf;%3-D Plot of T(r,t) vs x and t tau1=linspace(0,.2,101); [X,Y]=meshgrid(roverR,tau1); mesh(X,Y,framedata(1:101,:)) colormap(cool) text(.2,0,-.2,'Dimensionless Radius, r/R','Rotation',10) text(0,.17,-.25,'Time, \kappat/R^2','Rotation',-34) zlabel('Dimensionless Temperature, T(r,t)/T_0') view(330,30) axis([0 1 0 .2 0 1]) text(.05,.2,.9,'Press Enter to Continue','Rotation',10) pause