% conductionsphere.m Last modified 10/23/2008 % Conduction of heat in a sphere of radius R % with exterior face (r=R) at temperature of T0 % and initial interior at zero temperature. % The diffusivity is kappa=k/(rho*c) % This is the well-known "Cooking a Baked Potato Problem" problem posed when % investigating how long it will take to cook a potato starting at % temperature T0 and with a constant external temperature TH. % 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); roverR=linspace(0,1,101);roverR(1,1)=.0001; nterm=20;% number of terms in the Fourier series solution. for i=1:nterm phi(i,:)=sin(i*pi*roverR); end tau=linspace(0,.25,201); T1=(ones(1,101)); F=ones(1,101)./roverR; framedata=[]; for k=1:201%time loop sum=zeros(1,101); for i=1:nterm sum=sum+(((-1)^i)/i)*phi(i,:)*exp(-(i*i*pi*pi*tau(1,k))); end sum=sum*(2/pi).*F; 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 axis([-.2 1.2 -.2 1.2]) hold on 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]/[T_H-T_0]') for k=21:20:201% time loop plot(roverR,framedata(k,:),'k') hold on end text(.7,.05,['\kappat/R^2 = ' num2str(0)]) text(.65,.2,num2str(0.025)) text(.6,.3,num2str(0.05)) text(0.3,.9,num2str(0.25)) set(gca,'Box','on') text(.1,1.13,'Press Enter to Continue') pause hold off figure(2);clf;%Plot T(x,t) vs t for various x axis([0 .25 -.2 1.2]); hold on plot(tau,ones(1,201)) hold on for j=1:10:101% space loop plot(tau,framedata(:,j)) hold on end text(.02,.95,['r/R = ' num2str(1)]) text(.04,.87,['r/R = ' num2str(0.9)]) text(.125,.35,['r/R = ' num2str(0)]) xlabel('Dimensionless Time, \kappat/R^2') ylabel('Dimensionless Temperature, [T(r,t)-T_0]/[T_H-T_0]') text(.02,1.1,'Press Enter to Continue') set(gca,'Box','on') pause hold off 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]/[T_H-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 k=2:201% time loop set(L,'Ydata',framedata(k,:)) pause(.05) end hold on plot(roverR,framedata(1,:),'k') hold on 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,.25,101); plotdata=[]; plotdata(1,:)=framedata(1,:); for i=2:101 plotdata=[plotdata ;framedata((2*i)-1,:)]; end [X,Y]=meshgrid(roverR,tau1); mesh(X,Y,plotdata,'Edgecolor',[0 0 0]) newmapdata=colormap; newmapdata=ones(size(newmapdata)); colormap=newmapdata; %colormap(cool) text(.2,0,-.2,'Dimensionless Radius, r/R','Rotation',10) text(0,.17,-.25,'Time, \kappat/R^2','Rotation',-34) zlabel('[T(r,t) - T_0]/[T_H - T_0]') view(330,30) axis([0 1 0 .25 0 1]) hold off text(.05,.25,.93,'Press Enter to Quit','Rotation',10) pause