% conductionspheresinusoid.m Last modified 7/31/2008 % Conduction of heat in a sphere of radius R. % The exterior face (r=R) has a sinusoidal variation in temperature % T0*cos(omega*t). The steady-state sinusoidal temperature interior % to the sphere is illustrated. % The diffusivity is kappa=k/(rho*c) and the input quantity is the % dimensionless frequency (R^2)*omega/kappa. % 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); freqpar=-1; while freqpar<=0 freqpar=input('The dimensionless frequency is (R^2)*omega/kappa = '); end alpha=sqrt(freqpar/2)*(1+1i); roverR=linspace(0,1,101); roverR(1,1)=.0001; omegat=linspace(0,6*pi,301); T1=(ones(1,101)); F=(ones(1,101)./roverR)/(exp(alpha)-exp(-alpha)); G=F.*(exp(alpha*roverR)-exp(-alpha*roverR)); phase=angle(G); mag=abs(G); framedata=[]; for t=1:301 framedata(t,:)=mag.*cos(omegat(1,t)*ones(1,101)+phase); end figure(1);clf;% Plot T(r,t) vs r for various t axis([0 1 -1.2 1.2]) hold on xlabel('Dimensionless Radius, r/R') hold on ylabel('Dimensionless Temperature, T(r,t)/T_0') for t=1:10:101% time loop plot(roverR,framedata(t,:),'k') hold on if t<=21 xline=[.9-.01*(t-1) .8-.01*(t-1)]; yline=[framedata(t,91-t+1) framedata(t,91-t+1)+.1]; plot(xline,yline,'k') text1=num2str(.02*(t-1)); text(.65-.01*(t-1),framedata(t,91-t+1)+.1,['\omegat = ' text1 '\pi']) end end text(.2, -1.1,['Dimensionless Frequency = R^2\omega/\kappa = ' num2str(freqpar)]); set(gca,'Box','on') text(.2,1.1,'Press Enter to Continue') hold off pause figure(2);clf;%Plot T(x,t) vs t for various x axis([0 6*pi -1.2 1.2]); hold on for j=1:10:101% space loop plot(omegat,framedata(:,j)) hold on end for j=101:-10:71 plot([omegat(1,201) 15.8], [framedata(201,j) framedata(201,j)+.1],'k') text(15.9,framedata(201,j)+.1,['r/R = ' num2str(roverR(1,j))]) end text(1,-1.1,['Dimensionless Frequency = R^2\omega/\kappa = ' num2str(freqpar)]) xlabel('Dimensionless Time, \omegat') ylabel('Dimensionless Temperature, T(r,t)/T_0') text(1,1.1,'Press Enter to Continue') set(gca,'Box','on') hold off pause figure(3);clf;%now do animation axis([0 1 -1.2 1.2]) hold on 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'); text(.2,-1.1,['Dimensionless Frequency = R^2\omega/\kappa = ' num2str(freqpar)]) texthandl=text(.2,1.1,'Press Enter to Animate'); pause set(texthandl,'String',' '); for k=2:301% time loop set(L,'Ydata',framedata(k,:)) pause(.05) end 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,omegat(1,1:101)); mesh(X,Y,framedata(1:101,:)) colormap(cool) text(.2,0,-1.5,'Dimensionless Radius, r/R','Rotation',12) text(0,6,-1.45,'Dimensionless Time, \omegat','Rotation',-33) zlabel('Dimensionless Temperature, T(r,t)/T_0') text(1,5.5,.9,'Press Enter to Quit','Rotation',-33) view(330,30) axis([0 1 0 2*pi -1.1 1.1]) pause