%
%   RECT1D_FIN_2.M   Introduction to FD Methods for Solving BVPs
%                Heat Transfer Analysis of a Rectangular Fin Arrangement
%
%   This file illustrate the finite difference (FD) method for solving BVPs.  
%   We solve the same rectangular fin heat transfer problem that was treated in 
%   program rect1d_fin_1.m.  Here, the focus is on comparing the FD solution
%   to the analytical solution for the temperature profile, T(x), along the
%   length of the fin.  The development of the FD equations programmed here is
%   contained within the Lecture Notes (see fd_intro.pdf).
%
%   File prepared by J. R. White, UMass-Lowell  (last update: Oct. 2017)
%

      clear all,  close all,  nfig = 0;
%
%   identify basic problem data
      w = 1;                            % unit width of fin (m)
      thk = 0.01;                       % fin thickness (m)
      L = 0.05;                         % fin length  (m)
      Tb = 200;                         % fin base temperature (C)
      Tinf = 30;                        % environment temperature (C)
      h = 500;                          % heat transfer coeff  (W/m^2-C)
      k = 200;                          % fin thermal conductivity (W/m-C)
%
%   compute some derived parameters 
      P = 2*w + 2*thk;                  % perimeter
      Ac = w*thk;                       % cross section area (conduction area)
      m = sqrt(h*P/(k*Ac));             % constant in derived equations (see notes)
      m2 = m*m;                         % constant in derived equations (see notes)
      hmk = h/(m*k);                    % constant in derived equations (see notes)
      bot = cosh(m*L) + hmk*sinh(m*L);  % constant in derived equations (see notes)
%
%       *************************
%       *  Analytical Solution  *
%       *************************
%
      Nx = 51;  xa = linspace(0,L,Nx)';  % independent spatial variable
      Ta = Tinf + (Tb - Tinf)*(cosh(m*(L-xa)) + hmk*sinh(m*(L-xa)))/bot;
%
%       *****************
%       *  FD Solution  *
%       *****************
%
%   Setup and solve matrix eqn for discrete temperatures
      N = input('Enter number of unknowns (N) for problem: ');
      dx = L/N;    dx2 = dx*dx;   xn = dx:dx:L;   xn = xn';
      A = zeros(N,N);     b = zeros(N,1);
%   node 1
      A(1,1) = -(2 + m2*dx2);  A(1,2) = 1;  b(1) = -m2*dx2*Tinf - Tb;
%   interior nodes
      for i = 2:N-1
        A(i,i-1) = 1;  A(i,i) = -(2 + m2*dx2);  A(i,i+1) = 1;  b(i) = -m2*dx2*Tinf;
      end
%   last node
      A(N,N-2) = 1;  A(N,N) = -(2*m2*dx2 + 1 + 2*h*dx/k);
      b(N) = -(2*h*dx/k + 2*m2*dx2)*Tinf;
%   solve
      Tn = A\b;   
%   add on the known base temperature to the solution vector
      xn = [0; xn];   Tn = [Tb; Tn];
%
%       *************************
%       *  Plot both solutions  *
%       *************************
%
      nfig = nfig+1;   figure(nfig)
      plot(xa*100,Ta,'r-',xn*100,Tn,'bo','LineWidth',2),grid
      title('RECT1D\_FIN\_2: T(x) Profile in Rectangular Fin')
      xlabel('Spatial Position (cm)')
      ylabel('Temperature (^oC)')
      legend('Analytical Solution',['FD Solution (N = ',num2str(N),')'])
%
%  end of problem
