% conduction3.m Last modified 7/17/08 % Conduction of heat in a slab with left face at T0 and right face at zero % and slab initially at zero temperature. The problem is solved % approximately by finite differences in space and state variables in time. % 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); xoverL=linspace(0,1,21); N=19;% 19 internal nodes A=zeros(19,19); B=zeros(19,1); B(1,1)=400; for i=1:19 A(i,i)=-800; end for i=1:18 A(i,i+1)=400; A(i+1,i)=400; end T=.001; [Ad,Bd]=c2d(A,B,T); x=zeros(19,1); framedata=x'; for k=1:200 x=Ad*x+Bd; framedata=[framedata;x']; end framedata=[ones(201,1) framedata zeros(201,1)]; framedata(1,1)=0; 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) 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:20:201 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;%Plot T(x,t) vs t for various x tau=linspace(0,.2,201); plot(tau,ones(1,201)) hold on for j=3:2:21 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') hold off text(.02,1.1,'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 = 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:201 set(L,'Ydata',framedata(i,:)) pause(.02) end hold on plot(xoverL,framedata(1,:),'k') hold off set(texthandl,'String','Press Enter to Continue'); pause figure(4);clf;%Color plot axis([-.2 1.2 -.2 1.2 ]); hold on framedata(:,1)=ones(201,1); colormap(cool); zdata=zeros(21,21); for j=1:21 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:2:200 for j=1:21 zdata(j,:)=framedata(i,:); end Q=pcolor(xoverL,xoverL,zdata); pause(.005) end hold off 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',-32) ylabel('Time, \kappat/L^2','Rotation',12) zlabel('Dimensionless Temperature, T(x,t)/T_0') text(.2,.2,.9,'Press Enter to Continue','Rotation',-32) view(60,30) axis([0 1 0 .2 0 1])