% Conduction of heat in a slab with left and right faces both at zero % and slab initially at T0. 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); figure(1);clf figure(2);clf figure(3);clf figure(4);clf figure(5);clf 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:2:nterm sum=sum+((4/((i)*pi))*exp(-(((i)*pi)^2)*tau(1,k))*sn(i,:)); end framedata=[framedata;sum]; end framedata(1,:)=ones(1,101); figure(1);clf;% 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+.1) 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,.95,['\kappat/L^2 = ' num2str(0)]) text(.16,.8,num2str(0.02)) text(.4,.1,num2str(0.2)) text(1.04,.8,'T = 0') text(-.15,.8,'T = 0') set(gca,'Box','on') text(.1,0,'Press Enter to Continue') pause figure(2);clf;%Plot T(x,t) vs t for various x plot(tau,ones(1,101)) hold on for j=11:10:51 plot(tau,framedata(:,j)) hold on end text(.12,.1,['x/L = ' num2str(.1)]) text(.08,.3,num2str(0.2)) text(.05,.6,num2str(.3)) text(0.035,0.8, num2str(.4)) text(.03,.95,num2str(.5)) xlabel('Dimensionless Time, \kappat/L^2') ylabel('Dimensionless Temperature, T(x,t)/T_0') text(.1,.9,'Press Enter to Continue') pause figure(3);clf; %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 = 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,0,'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 on set(texthandl,'String','Press Enter to Continue'); pause figure(4);clf;%Color plot axis([-.2 1.2 -.2 1.2 ]); hold on framedata(:,1)=zeros(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=2:60 for j=1:101 zdata(j,:)=framedata(i,:); end Q=pcolor(xoverL,xoverL,zdata); pause(.005) end set(texthandl2,'String','Press Enter to Continue'); pause figure(5);clf;%3-D Plot of T(x,t) vs x and t [X,Y]=meshgrid(xoverL,tau); mesh(X,Y,framedata) colormap(cool) xlabel('Distance, x/L','Rotation',-34) ylabel('Time, \kappat/L^2','Rotation',12) zlabel('Dimensionless Temperature, T(x,t)/T_0') view(60,30) axis([0 1 0 .2 0 1])