% conduction4.m Last modified 7/17/2008 %Conduction of heat in a slab with convection at both surfaces % solved by finite differeences in space and state variables in time, % using 21 nodes counting the faces of the slab. % 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); xoverL=[-.1 xoverL 1.1]; N=20;% 19 internal nodes h1Loverk=1; h2Loverk=1; A=zeros(21,21); B=zeros(21,1); B(1,1)=2*h1Loverk*N; for i=2:20 A(i,i)=-2*N*N; end for i=2:20 A(i,i+1)=N*N; A(i+1,i)=N*N; end A(1,1)=-2*N*N-2*N*h1Loverk; A(1,2)=2*N*N; A(2,1)=N*N; A(21,21)=-2*N*N-2*N*h2Loverk; A(21,20)=A(1,2); %A(20,21)=A(2,1); T=.001; [Ad,Bd]=c2d(A,B,T); x=zeros(21,1); framedata=x'; for k=1:2000 x=Ad*x+Bd; if fix(k/10)==k/10; framedata=[framedata;x']; end 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(.18,.2,num2str(0.2)) text(.22,.6,num2str(2)) text(1.04,.8,'T = 0') text(-.15,.8,'T = T_0') set(gca,'Box','on') text(.5,.7,['h_1L/k = ' num2str(h1Loverk)]); text(.5,.6,['h_2L/k = ' num2str(h2Loverk)]); text(-.05,.4,'h_1'); text(1.02,.4,'h_2'); 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=2:2:22 plot(tau,framedata(:,j)') hold on end text(1.2,.25,['x/L = ' num2str(1)]) text(1.2,.66,['x/L = ' num2str(0)]) %text(.082,.42,['x/L = ' num2str(0.2)]) %text(0.09,0.4,['x/L = ' num2str(.2)]) xlabel('Dimensionless Time, \kappat/L^2') ylabel('Dimensionless Temperature, T(x,t)/T_0') hold off text(.4,.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 = T_0') text(1.05,.8,'T = 0') text(-.05,.4,'h_1'); text(1.02,.4,'h_2'); 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 on text(.5,.7,['h_1L/k = ' num2str(h1Loverk)]); text(.5,.6,['h_2L/k = ' num2str(h2Loverk)]); 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,23); for j=1:21 zdata(j,:)=framedata(1,:); end Q=pcolor(xoverL,linspace(0,1,21),zdata); xlabel('Dimensionless Distance, x/L'); hold on ylabel('Dimension y'); texthandl2=text(0,1.1,'Press Enter to Continue'); pause set(texthandl2,'String',' '); y=linspace(0,1,21); for i=2:2:200 for j=1:21 zdata(j,:)=framedata(i,:); end Q=pcolor(xoverL,y,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) axis([-.1 1.1 0 2 0 1.1]); hold on text(0.05,2,1,'Press Enter to Continue','Rotation',-30); text(0,0,-.2,'Location, x/L','Rotation',-30) text(1.1,.4,-.2,'Dimensionless Time, \kappat/L^2','Rotation',12) zlabel('Dimensionless Temperature, T(x,t)/T_0') view(60,30)