% conductionslabsinusoid.m Last modified 7/31,2008 % Conduction of heat in a slab of thickness L. % The left face (x=0) has a sinusoidal variation in temperature % T0*cos(omega*t). The face at x=L is insulated. The steady-state % sinusoidal temperature interior of the slab is illustrated. % The diffusivity is kappa=k/(rho*c) and the input quantity is the % dimensionless frequency (L^2)*omega/kappa. % 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); freqpar=-1; while freqpar<=0 freqpar=input('The dimensionless frequency is (L^2)*omega/kappa = '); end alpha=sqrt(freqpar/2)*(1+1i); xoverL=linspace(0,1,101); omegat=linspace(0,6*pi,301); T1=(ones(1,101)); F=(exp(-alpha*(ones(1,101)-xoverL))+exp(alpha*(ones(1,101)-xoverL))); G=F/(exp(-alpha)+exp(alpha)); phase=angle(G); mag=abs(G); framedata=[]; for t=1:301 framedata(t,:)=mag.*cos(omegat(1,t)*ones(1,101)+phase); end figure(1);clf;% Plot T(x,t) vs x for various t axis([0 1 -1.2 1.2]) hold on xlabel('Dimensionless Distance, x/L') hold on ylabel('Dimensionless Temperature, T(x,t)/T_0') for t=1:10:101% time loop plot(xoverL,framedata(t,:),'k') hold on if t<=21 xline=[.9-.01*(t-1) .8-.01*(t-1)]; yline=[framedata(t,91-t+1) framedata(t,91-t+1)+.1]; plot(xline,yline,'k') text1=num2str(.02*(t-1)); text(.65-.01*(t-1),framedata(t,91-t+1)+.1,['\omegat = ' text1 '\pi']) end end text(.2, -1.1,['Dimensionless Frequency = L^2\omega/\kappa = ' num2str(freqpar)]); set(gca,'Box','on') text(.2,1.1,'Press Enter to Continue') hold off pause figure(2);clf;%Plot T(x,t) vs t for various x axis([0 6*pi -1.2 1.2]); hold on for j=1:10:101% space loop plot(omegat,framedata(:,j)) hold on end for j=101:-10:71 plot([omegat(1,201) 15.8], [framedata(201,j) framedata(201,j)+.1],'k') text(15.9,framedata(201,j)+.1,['r/R = ' num2str(xoverL(1,j))]) end text(1,-1.1,['Dimensionless Frequency = L^2\omega/\kappa = ' num2str(freqpar)]) xlabel('Dimensionless Time, \omegat') ylabel('Dimensionless Temperature, T(x,t)/T_0') text(1,1.1,'Press Enter to Continue') set(gca,'Box','on') hold off pause figure(3);clf;%now do animation axis([0 1 -1.2 1.2]) hold on 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'); text(.2,-1.1,['Dimensionless Frequency = L^2\omega/\kappa = ' num2str(freqpar)]) texthandl=text(.2,1.1,'Press Enter to Animate'); pause set(texthandl,'String',' '); for k=2:301% time loop set(L,'Ydata',framedata(k,:)) pause(.05) end set(texthandl,'String','Press Enter to Continue'); pause figure(4);clf;%3-D Plot of T(x,t) vs x and t [X,Y]=meshgrid(xoverL,omegat(1,1:101)); mesh(X,Y,framedata(1:101,:)) colormap(cool) text(.2,0,-1.5,'Dimensionless Distance, x/L','Rotation',12) text(0,6,-1.45,'Dimensionless Time, \omegat','Rotation',-33) zlabel('Dimensionless Temperature, T(x,t)/T_0') view(330,30) axis([0 1 0 2*pi -1.1 1.1]) text(.2,2*pi,.9,'Press Enter to Continue','Rotation',12) pause