%
%   PIPE_INSULATION_1.M    Function evaluation and plotting in Matlab 
%             A Simple Heat Transfer Analysis: Circular Pipe with Insulation
%
%   This file evaluates and plots the temperature profiles in a simple annular
%   region in the range ri < r < ro, where ri and ro are the inside and outside 
%   radii of a solid annulus.  The inside boundary has a fixed temperature, Ti.
%   The outside surface is exposed to a convective environment, with heat 
%   transfer coefficient h and bulk fluid temperature Tinf.  
%
%   A normalized temperature is defined (see notes) so that the resulting 
%   expression highlights only three design parameters  --  k, h, and ro  --  and 
%   a series of cases is analyzed to determine the effect of each of these 
%   parameters on the temperature profile and on the overall heat transfer from 
%   this particular system (relative to the reference case).
%
%   The goals of this exercise are:
%     1. to become comfortable with function evaluation and plotting in Matlab,
%     2. to evaluate and edit some important quantitative characteristics of 
%        a particular system for comparison purposes (ie. sensitivity studies), and
%     3. to hopefully learn a little bit about conduction and convection heat 
%        transfer in a simple 1-D system.
%
%   The basic idea for this problem was obtained from the discussion in Chapter 2
%   of the Heat Transfer text by J. P. Holman (5th Ed.)  (1981, McGraw-Hill).
%
%   File prepared by J. R. White, UMass-Lowell  (last update: August 2017)
%

      clear all,  close all,  nfig = 0;
%%
%   identify basic problem data (reference data)
      kref = 0.1;                % thermal conductivity of insulation (W/m-C)
      href = 5;                  % heat transfer coeff for reference cooling (W/m^2-C)
      ri = 0.01;                 % inside radius of insulation (m)
      roref = 0.01635;           % outside radius of insulation (m) (about 1/4 inch)
%%
%   compute reference normalized temperature profile (use vector arithmetic)
      Nr = 50;  r = linspace(ri,roref,Nr);
      uref = 1 - log(r/ri)./(log(roref/ri) + kref/(href*roref));
%
%   now compute the normalized profile for three separate perturbation cases
%   case 1: high thermal conductivity
      k1 = 1.0;     % thermal conductivity of (relatively poor) insulation (W/m-C)
      u1 = 1 - log(r/ri)./(log(roref/ri) + k1/(href*roref));
%   case 2: higher heat transfer coeff
      h2 = 50;      % heat transfer coeff for high cooling rate (W/m^2-C)
      u2 = 1 - log(r/ri)./(log(roref/ri) + kref/(h2*roref));
%   case 3: more insulation
      ro3 = 0.0608; % outside radius of insulation (m) (about 2 inches)
      r3 = linspace(ri,ro3,Nr);
      u3 = 1 - log(r3/ri)./(log(ro3/ri) + kref/(href*ro3));
%%
%   now let's plot these comparisons 
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(r*100,uref,'r-',r*100,u1,'g--',r*100,u2,'b:','LineWidth',2),grid
      title('Normalized Temperature Profiles for Pipe Insulation Problem')
      ylabel('Normalized Temperature')
      v1 = axis;  v1(2) = 1.65;  axis(v1)
      legend('Ref. Case','Higher k','Higher h','Location','SouthWest')
%
      subplot(2,1,2),plot(r*100,uref,'r-',r3*100,u3,'g--','LineWidth',2),grid
      xlabel('Radial Location (cm)'),ylabel('Normalized Temperature')
      legend('Ref. Case','More Insulation')
      v2 = axis;  v2(2) = 6;  axis(v2)
%%
%   just for fun, let's compute the heat transfer for each of these cases
%   Note: we will compute q two different ways for each case just to show
%         that they are equivalent...
%               qcond = conduction    and   qconv = convection
%   Since there is no internal heat source qcond is not a function of r (see notes).
%   Also, the right boundary condition forces qcond = qconv.
%   Let's see if this all works as expected...
%
%   ref case
      qcondref = kref/(log(roref/ri) + kref/(href*roref)); 
      qconvref = href*roref*uref(Nr);
%   case 1
      qcond1 = k1/(log(roref/ri) + k1/(href*roref)); 
      qconv1 = href*roref*u1(Nr);
%   case 2
      qcond2 = kref/(log(roref/ri) + kref/(h2*roref)); 
      qconv2 = h2*roref*u2(Nr);
%   case 3
      qcond3 = kref/(log(ro3/ri) + kref/(href*ro3)); 
      qconv3 = href*ro3*u3(Nr);
%   compute and edit desired ratios
      ratio1 = [qcond1/qcondref qconv1/qconvref];
      ratio2 = [qcond2/qcondref qconv2/qconvref];
      ratio3 = [qcond3/qcondref qconv3/qconvref];
      fprintf(1,'\n       Summary Results for the Pipe Insulation Problem \n\n');
      fprintf(1,'  Case #1 (higher k)       : qcond1/qcondref = %6.3f  and  qconv1/qconvref = %6.3f \n',ratio1);
      fprintf(1,'  Case #2 (higher h)       : qcond2/qcondref = %6.3f  and  qconv2/qconvref = %6.3f \n',ratio2);
      fprintf(1,'  Case #3 (more insulation): qcond3/qcondref = %6.3f  and  qconv3/qconvref = %6.3f \n',ratio3);
%
%  end of problem
