%
%   CYLFINA.M   Heat Transfer in a Cylindrical Fin (Analytical Solution)
%
%   This file solves the cylindrical fin heat transfer problem using modified 
%   Bessel functions. The base problem is defined via the following equation:
%
%     x^2*u'' + x*u' - ALF2*x^2*u = 0    where   ALF2 = 2*h*rs^2/[k*thk]
%
%     with B.C.   u(a) = 1    and    u'(b) = 0
%      and        a = rw/rs   and    b = rs/rs = 1
%
%     where   u = normalized temp = [T(r) - Tinf]/[Tw - Tinf]
%             x = normalized distance = r/rs
%
%     with   rw, rs = inside and outside radius of fin, respectively
%            thk = thickness of fin
%     and    h, k, Tw, and Tinf are all given quantities (fixed)
%
%   From the normalized solution we can construct absolute profiles (if desired):
%       T(r)  =  Tinf + u(x)[Tw - Tinf]
%       T'(r) =  u'(x)[Tw - Tinf]/rs
%
%   The above development is given as part of the course notes in the Math
%   Methods course (10/24.539).
%
%   File prepared by J. R. White, UMass-Lowell (Aug. 2003)
%

%
%   getting started
      clear all, close all, nfig = 0;
%
%   basic data for the problem
      rw = 1/12;    rs = 1.5/12;    % inside and outside radius (ft)
      thk = .0625/12;               % thickness of fin (ft)
      Tw = 200;    Tinf = 70;       % inside wall and ambient temps (F)
      h = 20;                       % heat transfer coeff (BTU/hr-ft^2-F)
      k = 75;                       % thermal conductivity (BTU/hr-ft-F)
%
%   derived constants
      a = rw/rs;   b = rs/rs;
      alf2 = (2*h*rs*rs)/(k*thk);    alpha = sqrt(alf2);
      Qideal = 2*pi*h*(rs*rs-rw*rw)*(Tw-Tinf);
%
%   write base data to output file
      fid = fopen('cylfina.out','w');
      fprintf(fid,'\n *** CYLFINA.OUT ***  Data and Results from CYLFINA.M \n');
      fprintf(fid,'\n \nBASIC DATA FOR PROBLEM \n');
      fprintf(fid,'Inside and outside radius:  rw =  %6.2f ft \t rs =   %5.2f ft \n',rw,rs);
      fprintf(fid,'Thickness of fin:           thk = %6.2f ft \n',thk);
      fprintf(fid,'Inside/ambient temps:       Tw =  %6.2f F  \t Tinf = %5.2f F \n',Tw,Tinf);
      fprintf(fid,'Heat transfer coeff:        h =   %6.2f Btu/hr-ft^2-F \n',h);
      fprintf(fid,'Thermal conductivity:       k =   %6.2f Btu/hr-ft-F \n',k);
%
%   define solution domain
      x = linspace(a,b,50);       [nr,nc] = size(x);
%
%   evaluate constants in temperature equation 
      a1 = besselk(1,alpha*b);      a2 = besseli(1,alpha*b);
      b1 = besseli(0,alpha*a);      b2 = besselk(0,alpha*a);
      denom = a1*b1+a2*b2;
%
%   create normalized temp and gradient
      u = (a1*besseli(0,alpha*x) + a2*besselk(0,alpha*x))/denom;
      up = alpha*(a1*besseli(1,alpha*x) - a2*besselk(1,alpha*x))/denom;
%
%   calc heat transferred (conduction) & fin eff. (edit key parameters)
      Qactual = -k*(2*pi*a*thk)*(Tw-Tinf)*up(1);
      Ttip = Tinf + (Tw-Tinf)*u(nc);
      fineff = Qactual/Qideal;
      fprintf(fid,'\n \nANALYTICAL SOLUTION RESULTS \n');
      fprintf(fid,'Wall Temp (F) =     %8.3f \n',Tw);
      fprintf(fid,'Ambient Temp (F) =  %8.3f \n',Tinf);
      fprintf(fid,'Tip Temp (F) =      %8.3f \n',Ttip);
      fprintf(fid,'Qideal (BTU/hr) =   %8.3f \n',Qideal);
      fprintf(fid,'Qactual (BTU/hr) =  %8.3f \n',Qactual);
      fprintf(fid,'Fin Eff =           %8.3f \n',fineff);
%
%   plot normalized profiles
      nfig = nfig+1;    figure(nfig)
      subplot(2,1,1)
      plot(x,u,'LineWidth',2),grid
      title('CylFinA:  Normalized Temp Profile for Cylindrical Fin (Analytical)')
      ylabel('temperature')
      subplot(2,1,2)
      plot(x,up,'LineWidth',2),grid
      title('CylFinA:  Normalized Temp Gradient for Cylindrical Fin (Analytical)')
      xlabel('normalized distance'),ylabel('temp gradient')
%
%   close output file
      fclose(fid);
%
%   end simulation
