%
%   Fuel_Temp_FD2.M   Solve for T(r) in a cylindrical fuel pin 
%                 Analytically and via the FD Method  (variable kf)
%                                                                             
%   This demo uses the finite difference technique as the solution method                 
%   for a particular heat transfer problem  -->  determine T(r) in a 
%   cylindrical fuel pin in a nuclear reactor when the thermal conductivity
%   of the UO2 fuel is temperature dependent.  The mesh interval is 
%   constant.  This follows the development done in class (see notes).
%
%   For the variable thermal conductivity case, an iterative scheme is required 
%   to solve the nonlinear system.  The basic scheme is as follows:
%      1. guess the temperature profile (here let T(r) = constant for simplicity)
%      2. compute K(r) with K(T) correlation using current estimate of T(r)
%      3. setup matrix balance equations and solve for new T(r)
%      4. check convergence (go back to Step 2 if not converged)
%   This scheme should converge rather quickly since the linear dependence of K(T)
%   only introduces a mild nonlinearity (over a reasonable temperature range). 
%
%   Function Fuel_Temp_FD2_mx.m is used to set up the matrix equations.
%  
%   File prepared by J. R. White, UMass-Lowell (Oct. 2003).
%

%
%   getting started
      clear all, close all,  nfig = 0;
      global  qqq Ts Ks h h2
%
%   set properties
      qqq = 1.35e+7;           % average volumetric heat source (BTU/hr-ft^3)
      Rf = 0.0175;             % fuel pin radius (ft)
      Ts = 650;                % surface temperature (F)
      a1 = 94.4;  a2 = -0.556; % paramters for kf UO2 (kf in BTU/hr-ft-F)
%
%   Analytical Solution for variable kf  
      ro = 0;   rf = Rf;   
      re = linspace(ro,rf,30);
      c = 1+a2;  Te = (Ts^c + (c/a1)*(qqq/4)*(Rf*Rf - re.*re)).^(1/c);
%
%   Solution via Finite Difference Method (numerical approx)
%
%   set domain limits. define # of points, and set r-vector 
      N = input('Input number of unknowns in problem (N): '); 
      dr = (rf-ro)/N;   h = dr;   h2 = h*h;   r = ro:h:rf-h;   r = r';
%
      Told = Ts*ones(length(r),1);      % initial guess for temp profile
      K = a1*Told.^a2;                  % initial guess for K(T)->K(r) 
      Ks = a1*Ts^a2;                    % kf at surface temp 
%
      itmax = 10;   itend = itmax;   it = 1;   tol = 1e-5;
      while it <= itmax
%
%   compute temperature profile for given K profile
        [AA,BB] = fuel_temp_fd2_mx(r,K);  % determine coeff matrices
        Tnew = AA\BB;                     % find solution vector
%
%   check convergence
        emax = max(abs(Tnew-Told)./Tnew);
%
        fprintf(1,'\n Intermediate edit for Variable kf Case\n')
        fprintf(1,'  it = %3d      max error = %8.2e \n',it,emax)
        fprintf(1,'     Tnew         Told         K used\n')
        for j = 1:N
          fprintf(1,'  %8.2f     %8.2f     %8.2f\n',Tnew(j),Told(j),K(j))
        end
%
        if emax < tol;
          itend = it;   it = itmax+1; 
        else 
          it = it+1;    K = a1*Tnew.^a2;   Told = Tnew;
        end
      end
%
%   print max relative error and iteration count
      fprintf(1,'\n\n Number of iterations to convergence = %3d\n',itend)
      fprintf(1,' Max relative error at convergence =  %8.2e\n',emax)
%
%   save K(T) -> K(r) and T(r) profiles
      KT = K;   T = Tnew;
%
%   add any necessary boundary values to solution for plotting                             
      rp = [r; rf];   Tp = [T; Ts]; 
%                                                                             
%   plot results 
      nfig = nfig+1;  figure(nfig)                                                             
      plot(re,Te,'r-',rp,Tp,'bo','LineWidth',2),grid 
      title('Fuel Pin Radial Temperature Profile')                  
      xlabel('Radial Distance (ft)'),ylabel('Temperature (F)')    
      legend('Analytical Solution',['FD Method (FD2) (N = ',num2str(N),')'])                    
%
%   end of demo
%
