%
%   CYLFINFD.M   Heat Transfer in a Cylindrical Fin (Finite Difference Solution)
%
%   This file solves the cylindrical fin heat transfer problem using the finite
%   difference method. 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 2nd order ODE is first converted into a set of N simultaneous
%   algebraic equations (one for each node in the physical system) and then 
%   these are solved using Matlab's standard equation solver.  The mesh interval
%   is constant for this example.
%
%   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('cylfinfd.out','w');
      fprintf(fid,'\n *** CYLFINFD.OUT ***  Data and Results from CYLFINFD.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);
%
%   loop over different N if desired
      u0 = 1;
      N = input('Input number of unknowns/nodes for problem (zero to quit) = ');
      while N ~= 0
        clear A B u x xx zz zzp
        A = zeros(N,N);  B = zeros(N,1);
%
        dx = (b-a)/N;   dx2 = dx*dx;   x = linspace(a+dx,b,N);
        fn = (b-dx/4)/(b-dx/2);
%
%   setup matrix eqns. (treat boundary terms as special cases)
        for n = 2:N-1
          A(n,n-1) = 1-dx/(2*x(n));   A(n,n) = -(2+alf2*dx2);
          A(n,n+1) = 1+dx/(2*x(n));   B(n) = 0;
        end
%   left boundary (fixed temp)
        A(1,1) = -(2+alf2*dx2);   A(1,2) = 1+dx/(2*x(1));
        B(1) = -(1-dx/(2*x(1)))*u0;
%   right boundary (zero heat flux)
        A(N,N-1) = 1;     A(N,N) = -(1+alf2*dx2*fn/2);   B(N) = 0;
%
%   solve system of eqns
        u = A\B;
%
%   add left boundary point to solution for plotting
        xx = [a x];      zz = [u0 u'];
%
%   estimate temperature gradient for plotting
        zzp(1) = (zz(2)-zz(1))/dx;   zzp(N+1) = 0;
        for n = 2:N,   zzp(n) = (zz(n+1)-zz(n-1))/(2*dx);   end
%
%   calc heat transferred and fin efficiency (using cond. & conv. estimates)
        Qactual1 = -k*(2*pi*a*thk)*(Tw-Tinf)*zzp(1);   dr = dx*rs;
        Qactual2 = h*2*pi*((rw+dr/2)^2-rw^2)*(Tw-Tinf)*zz(1) +    ...
                   h*2*pi*(rs^2-(rs-dr/2)^2)*(Tw-Tinf)*zz(N+1);
        for n = 2:N
          r2 = (xx(n)+dx/2)*rs;   r1 = (xx(n)-dx/2)*rs;
          area = 2*pi*(r2^2-r1^2);
          Qactual2 = Qactual2+h*area*(Tw-Tinf)*zz(n);
        end
%
%   edit key parameters
        Ttip = Tinf + (Tw-Tinf)*zz(N+1);    fineff = Qactual2/Qideal;
        fprintf(fid,'\n \nFINITE DIFFERENCE METHOD RESULTS (N = %3i) \n',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 (conduction)\n',Qactual1);
        fprintf(fid,'Qactual (BTU/hr) =  %8.3f (convection)\n',Qactual2);
        fprintf(fid,'Fin Eff =           %8.3f \n',fineff);
%
%   plot normalized profiles
        nfig = nfig+1;   figure(nfig)
        subplot(2,1,1),plot(xx,zz,'LineWidth',2)
        title(['CylFinFD:  Normalized Temp Profile for Cylindrical Fin (FD Method N = ', ...
               num2str(N),')'])
        grid,ylabel('temperature')
        subplot(2,1,2),plot(xx,zzp,'LineWidth',2)
        title(['CylFinFD:  Normalized Temp Gradient for Cylindrical Fin (FD Method N = ', ...
               num2str(N),')'])
        grid,xlabel('normalized distance'),ylabel('temp gradient')
%
%   loop back for a different N if desired
        N = input('Enter another value of N if desired (zero to quit)? ');
      end
%
%   close output file
      fclose(fid);
%
%   end simulation
