% convboundaries.m Last modified 7/18/2008 % This module solves the conduction in a finite slab with convective % boundaries. The temperature at the left boundary layer is suddenly elevated % T0 degrees above the original equilibrium temperature. The Biot Number at % both boundaries is assumed the same (unity). % 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); % Below ids an alternative Biot number and the associated eigenvalues %Bi=.1; %lambdaL=[.44352 3.203994 6.314854 9.445951 12.58227 15.72069 18.86016]; %lambdaL=[lambdaL 22.00024 25.1407 28.2814 31.42229 34.56331 37.70442]; Bi=1;% Bi=Biot Number=hL/k lambdaL=[1.306542 3.673195 6.58462 9.631684 12.72324 15.83411]; lambdaL=[lambdaL 18.95497 22.08166 25.21203 28.34486 31.47944 34.61528]; xoverL=linspace(0,1,201); delta=.005;%delta(x/L), integration increment for i=1:12 z=lambdaL(1,i)*xoverL; phi(i,:)=(Bi/lambdaL(1,i))*sin(z)+cos(z); end %Check Orthogonality of Eigenfunctions for i=1:12 for j=1:12 I(i,j)=(sum(phi(i,:).*phi(j,:))-(phi(i,1)*phi(j,1)+phi(i,201)*phi(j,201))/2)*delta; end end den=Bi+2; A=(Bi+1)/den; B=-Bi/den; f=A*ones(1,201)+B*(xoverL);% Steady State Temperature in Slab sum1=zeros(1,201); for i=1:12 b(1,i)=(1/I(i,i))*(sum(-f.*phi(i,:))-(-f(1,1)*phi(i,1)-f(1,201)*phi(i,201))/2)*delta; sum1=sum1+b(1,i)*phi(i,:); end T=sum1+f; tempdata=[T]; T(1,1)=0; tau=linspace(0,2,201);% tau=kappa*t/(L^2), Dimensionless Time for t=2:201 sum2=zeros(1,201); for i=1:12 sum2=sum2+b(1,i)*exp(-(lambdaL(1,i)^2)*tau(1,t))*phi(i,:); end T=sum2+f; tempdata=[tempdata;T]; end figure(1);clf% Plot T vs time for 11 values of x/L axis([0 2 -.1 1.1]); hold on for i=1:20:201 plot(tau,tempdata(:,i)); hold on end box on ylabel('Dimensionless Temperature Change, T(x,t)/T_0') xlabel('Dimensionless Time, \kappat/L^2') text(1.7,.28,'x/L = 0'); text(1.7,.7,'x/L = 1'); text(1.7,.1,'\Deltax/L = 0.1'); hold off text(1, 1,'Press Enter to Continue'); pause xoverL1=[-.1 -.1 xoverL 1.1]; tempdata=[ones(201,1) ones(201,1) tempdata zeros(201,1)]; tempdata(1,:)=zeros(1,204); tempdata(1,1)=1; figure(2);clf% Plot T versus x/L for various times patch([0 0 1 1 0],[-.1 1.1 1.1 -.1 -.1],'c') hold on axis([-.2 1.2 -.1 1.1]) hold on for t=1:10:101 plot(xoverL1,tempdata(t,:)) hold on end text(.6,.9,['Biot No. = ' num2str(Bi)]) xlabel('Dimensionless distance, x/L') ylabel('Dimensionless Temperature Change, T(x,t)/T_0') text(.2,.025,'\kappat/L^2 = 0') text(.2,.1,'0.1') text(.2,.53,'1') box on hold off text(.1, 1,'Press Enter to Continue') pause figure(3);clf% Animate T vs x/L for closely spaced values of time axis([-.2 1.2 -.1 1.1]) hold on patch([0 0 1 1 0],[-.1 1.1 1.1 -.1 -.1],'c') hold on L=plot(xoverL1,tempdata(1,:),'EraseMode','xor'); ylabel('Dimensionless Temperature Change, T(x,t)/T_0') xlabel('Dimensionless Distance, x/L') text(.6,.9,['Biot No. = ' num2str(Bi)]) plot([-.2 1.2],[-.005 -.005]) hold on box on T1=text(.1,1,'Press Enter to Animate'); pause set(T1,'String',' ') for t=1:201 set(L,'Ydata',tempdata(t,:)); pause(.1) end box on hold off set(T1,'String','Press Enter to Continue'); pause figure(4);clf% 3-D plot of T vs x/L and t [X,Y]=meshgrid(xoverL1,tau); mesh(X,Y,tempdata); colormap(cool); axis([-.1 1.1 0 2 0 1.1]); hold on text(0,0,-.2,'Location, x/L','Rotation',-32) text(1.1,.4,-.2,'Dimensionless Time, \kappat/L^2','Rotation',12) zlabel('Dimensionless Temperature, T(x,t)/T_0') view(60,30) hold off text(0.05,2,.9,'Press Enter to Continue','Rotation',-32); pause