% conductionspheresource.m Last modified 7/31/2008 % Conduction of heat in a sphere of radius R % with exterior face (r=R) at temperature of 0 % and initial interior at zero temperature with a % heat source Q0 at the origin. Uses 40 terms of the Fourier Series. % 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); roverR=linspace(0.2,1,81); nterm=40;% number of terms in the Fourier series solution. for i=1:nterm phi(i,:)=sin(i*pi*(ones(1,81)-roverR)); end tau=linspace(0,.1,201); framedata=[]; for k=1:201%time loop sum=zeros(1,81); for i=1:nterm sum=sum+((-1)^(i+1))*(1/(i*pi))*phi(i,:)*exp(-(i*i*pi*pi*tau(1,k))); end sum=(ones(1,81)-roverR)-2*sum; T=sum./roverR; framedata=[framedata;T]; end Tinfinity=(ones(1,81)-roverR)./roverR; framedata(1,:)=zeros(1,81); %framedata(1,81)=0; figure(1);clf;% Plot T(r,t) vs r for various t axis([-.2 1.2 0 4]) hold on plot([0 0],[0,4]); hold on plot([1 1],[0,4]); hold on plot(roverR,framedata(1,:),'k') hold on xlabel('Dimensionless Radius, r/R') ylabel('Temperature,4\pikRT(r,t)/Q_0') for k=21:20:201% time loop plot(roverR,framedata(k,:),'k') hold on end plot(roverR,Tinfinity,'r','LineWidth',[1.5]); text(.36,2,['\kappat/R^2 = Infinity']); text(.02,.1,['\kappat/R^2 = ' num2str(0)]) text(.12,.5,num2str(0.01)) text(.38,1.2,num2str(0.1)) set(gca,'Box','on') text(.4,3.5,'Press Enter to Continue') hold off pause figure(2);clf;%Plot T(x,t) vs t for various x axis([0 tau(1,201) 0 4]); hold on hold on for j=1:10:81% space loop plot(tau,framedata(:,j)) hold on end text(.018,2.3,['r/R = ' num2str(0.2)]) text(.035,1.2,['r/R = ' num2str(0.3)]) %text(.125,.35,['r/R = ' num2str(0)]) xlabel('Dimensionless Time, \kappat/R^2') ylabel('Temperature,4\pikRT(r,t)/Q_0') text(.02,3.5,'Press Enter to Continue') set(gca,'Box','on') hold off pause figure(3);clf;%now do animation axis([-.2 1.2 -.2 4]) hold on plot([0 0],[-.2,4]); plot([1 1],[-.2,4]); box on xlabel('Dimensionless Radius, r/R') ylabel('Temperature,4\pikRT(r,t)/Q_0') L=plot(roverR, framedata(1,:),'k','EraseMode','xor'); hold on texthandl=text(.4,3.5,'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,Tinfinity,'r','LineWidth',[1.5]); text(.36,2,['\kappat/R^2 = Infinity']); 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(.10,0,-.8,'Location, r/R','Rotation',-32) text(1,.05,-.7,'Time, \kappat/R^2','Rotation',10) zlabel('Temperature,4\pikRT(r,t)/Q_0') view(60,30) axis([0 1 0 .2 0 3.5]) text(.1,.2,3.1,'Press Enter to Continue','Rotation',-32) pause