%
%   RECT1D_FIN_3.M   Introduction to FD Methods for Solving Nonlinear BVPs
%                Heat Transfer Analysis of a Rectangular Fin Arrangement with
%              a Temperature-Dependent Heat Transfer Coefficient (Nonlinear Eqns.)
%
%   This file illustrates the finite difference (FD) method for solving nonlinear   
%   BVPs.  We solve the same rectangular fin heat transfer problem that was treated 
%   previously in program RECT1D_FIN_2.M, except now we allow the heat transfer 
%   coefficient to vary with temperature.  This makes the BVP nonlinear and, upon
%   discretization, the system of coupled algebraic equations that results also 
%   becomes nonlinear.  We will use the Linearized Iteration method (as discussed
%   in the Lecture Notes) to solve this nonlinear problem.
%
%   The basic iterative procedure is as follows:
%      1. guess the temperature profile (here let T(x) = constant for simplicity)
%      2. compute h(x) with the correlation given using the current estimate of T(x)
%      3. set up the matrix balance equations and solve for new T(x)
%      4. check convergence (go back to Step 2 if not converged)
%   Hopefully this scheme will converge rather quickly (remember that this one-point
%   iteration scheme is not guaranteed to converge).  Let's see...
%
%   The development of the FD equations are done in the Lecture Notes and the 
%   coefficient matrices for each new T(x) guess are constructed in RECT1D_FIN_3A.M.
%
%   File prepared by J. R. White, UMass-Lowell  (last update:  Nov. 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)
      k = 200;                     % fin thermal conductivity (W/m-C)
      a1 = 200; a2 = 2;            % constants in heat transfer coeff correlation
%
%   define # of points and set x-vector 
      N = input(' Input number of unknowns in problem (N): '); 
      dx = L/N;   x = dx:dx:L;   x = x';
%
%   set initial guess and start iterative loop
      Told = 0.5*(Tb+Tinf)*ones(size(x));      
      itmax = 20;  it = 0;  tol = 1e-5;  merr = 1;
      while merr > tol && it <= itmax
        h = a1 + a2*Told;                               % determine h(T) -> h(x)
%
%   compute temperature profile for given h(x)
        [A,b] = rect1d_fin_3a(h,Tb,Tinf,k,L,w,thk,N);   % determine coeff matrices
        Tnew = A\b;                                     % find solution vector
        it = it+1;                                      % increment iteration counter
        merr = max(abs(Tnew-Told)./Tnew);               % compute max relative change
%
%   print intermediate edit if the number of eqns is relatively small
        if N < 15
          fprintf('\n Intermediate edit for RECT1D_FIN_3.M \n');
          fprintf('  it = %3d      max error = %8.2e \n',it,merr);
          fprintf('     Tnew         Told        h used   \n');
          fprintf('      (C)          (C)       (W/m^2-C) \n');
          for j = 1:N
            fprintf('  %8.2f     %8.2f     %8.2f \n',Tnew(j),Told(j),h(j));
          end
        end
%
        Told = Tnew;  %    set Told (i.e. the guesss) for the next pass
%
      end   % end of iteration loop
%
%   print max relative error and iteration count
      fprintf(1,'\n\n Number of iterations to convergence = %3d\n',it);
      fprintf(1,' Max relative error at convergence =  %8.2e\n',merr);
%
%   add on the known base temperature to the solution vector [do same for h(x)]
      xp = [0; x];   Tp = [Tb; Tnew];   hb = a1 + a2*Tb;  hp = [hb; h];
%
%  plot temperature profile  
      nfig = nfig+1;   figure(nfig)
      plot(xp*100,Tp,'r--','LineWidth',2),grid
      title('T(x) Profile in Rectangular Fin (Nonlinear Case)')
      xlabel('Spatial Position (cm)'), ylabel('Temperature (^oC)')
%
%  plot heat transfer coeff profile  
      nfig = nfig+1;   figure(nfig)
      plot(xp*100,hp,'r--',xp*100,500*ones(size(hp)),'b-','LineWidth',2),grid
      title('h(x) Profile in Rectangular Fin (Nonlinear Case)')
      xlabel('Spatial Position (cm)'), ylabel('Heat Transfer Coeff (W/m^2-^oC)')
      legend('variable h', 'constant h')
%
%  end of problem
