% conduction.m Last updated 7/15/2008 % One dimensional conduction of heat in a slab with left face at T0 and % right face at zero and slab initially 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. clc disp(' This script displays the tempreature in a slab with initial') disp(' temperature zero, temperature T0 at x=0 and temperature zero') disp(' at x=L. The diffusivity is kappa=k/(rho*c).') clear all set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); xoverL=linspace(0,1,101); sn=[]; nterm=40; for i=1:nterm sn(i,:)=sin(i*pi*xoverL); end tau=linspace(0,.2,101); T1=(ones(1,101)-xoverL); framedata=[]; for k=1:101 sum=zeros(1,101); for i=1:nterm sum=sum-((2/((i)*pi))*exp(-(((i)*pi)^2)*tau(1,k))*sn(i,:)); end T=T1+sum; framedata=[framedata;T]; end framedata(1,:)=zeros(1,101); figure(1);clf reset % Plot T(x,t) vs x for various t f1=[ 0 0.08 .2 .31 .43 .51 .58 .64 .76 .88 1]; xpatch1=[f1 ones(1,11)-f1 f1(1,1)]; y1=[-.1 -.04 -.16 -.12 -.05 -.14 -.07 -.11 -.16 -.12 -.08]; ypatch1=[y1 1.1*ones(1,11)+(y1) y1(1,1)]; patch(xpatch1,ypatch1,'c') hold on plot(xoverL,framedata(1,:),'k') hold on xlabel('Dimensionless Distance, x/L') hold on ylabel('Dimensionless Temperature, T(x,t)/T_0') for k=11:10:101 plot(xoverL,framedata(k,:),'k') hold on axis([-.2 1.2 -.2 1.2]) end text(.05,.05,['\kappat/L^2 = ' num2str(0)]) text(.16,.2,num2str(0.02)) text(.4,.54,num2str(0.2)) text(1.04,.8,'T = 0') text(-.15,.8,'T = T_0') set(gca,'Box','on') hold off text(.1,1.13,'Press Enter to Continue') pause figure(2);clf reset%Plot T(x,t) vs t for various x plot(tau,ones(1,101)) hold on for j=11:10:101 plot(tau,framedata(:,j)) hold on end text(.02,.97,['x/L = ' num2str(0)]) text(.03,.8,['x/L = ' num2str(0.1)]) text(.1,-.08,['x/L = ' num2str(1)]) text(0.04,0.6,['x/L = ' num2str(.2)]) xlabel('Dimensionless Time, \kappat/L^2') ylabel('Dimensionless Temperature, T(x,t)/T_0') text(.02,1.1,'Press Enter to Continue') pause figure(3);clf reset%now do animation xpatch=[0 0 .2 .3 .5 .6 .85 1 1 .82 .66 .43 .2 0]; ypatch=[-.1 1.1 1.15 1.05 1.13 1.1 1.14 1.1 -.1 -.03 -.14 -.1 -.16 -.1]; patch(xpatch,ypatch,'c') axis([-.2 1.2 -.2 1.2]) hold on text(-.15,.8,'T = T_0') text(1.05,.8,'T = 0') box on xlabel('Dimensionless Distance, x/L') hold on ylabel('Dimensionless Temperature, T(x,t)/T_0') hold on L=plot(xoverL, framedata(1,:),'k','EraseMode','xor'); hold on texthandl=text(.4,1,'Press Enter to Animate'); pause set(texthandl,'String',' '); for i=2:101 set(L,'Ydata',framedata(i,:)) pause(.05) end hold on plot(xoverL,framedata(1,:),'k') hold off set(texthandl,'String','Press Enter to Continue'); pause figure(4);clf reset%Color plot axis([-.2 1.2 -.2 1.2 ]); hold on framedata(:,1)=ones(101,1); colormap(cool); zdata=zeros(101,101); for j=1:101 zdata(j,:)=framedata(1,:); end Q=pcolor(xoverL,xoverL,zdata); xlabel('Dimensionless Distance, x/L'); hold on ylabel('Dimension y'); texthandl2=text(0,1.1,'Press Enter to Continue'); pause set(texthandl2,'String',' '); for i=3:2:61 for j=1:101 zdata(j,:)=framedata(i,:); end Q=pcolor(xoverL,xoverL,zdata); pause(.01) end hold off set(texthandl2,'String','Press Enter to Continue'); pause figure(5);clf reset%3-D Plot of T(x,t) vs x and t tau1=linspace(0,.1,51); [X,Y]=meshgrid(xoverL,tau1); mesh(X,Y,framedata(1:51,:),'Edgecolor',[0 0 0]) newmapdata=colormap; newmapdata=ones(size(newmapdata)); colormap=newmapdata; %mesh(X,Y,framedata(1:51,:)) %colormap(cool) text(.2,0,-.2,'Distance, x/L','Rotation',-34) text(1,.05,-.2,'Time, \kappat/L^2','Rotation',12) zlabel('Dimensionless Temperature, T(x,t)/T_0') view(60,30) axis([0 1 0 .1 0 1])