%
%   Fuel_Temp_FD1.M   Solve for T(r) in a cylindrical fuel pin (kf = const)
%                      Analytically and via the Finite Difference Method
%
%   Given differential eqn.:   T'' + (1/r)T' + qqq/kf =  0   
%                        with  T'(0) = 0  and   T(Rf) = Ts 
%                                                                             
%   This demo uses the finite difference technique as the solution method                 
%   for a particular ODE.  The mesh interval is constant.  This follows
%   the development done in class...
%  
%   File prepared by J. R. White, UMass-Lowell (Sept. 2003).
%

%
%   getting started
      clear all, close all,   nfig = 0;
%
%   set properties
      qqq = 1.35e+7;     % average volumetric heat source (BTU/hr-ft^3)
      kf = 1.1;          % thermal conductivity of UO2 (BTU/hr-ft-F)
      Rf = 0.0175;       % fuel pin radius (ft)
      Ts = 650;          % surface temperature (F)
%
%   Analytical Solution (exact)
%
      re = linspace(0,Rf,30);
      Te = Ts + (qqq/(4*kf))*(Rf*Rf - re.*re);
%
%   Solution via Finite Difference Method (numerical approx)
%
%   set domain limits. define # of points, and set r-vector 
      ro = 0;   rf = Rf;   
      N = input('Input number of unknowns in problem (N): '); 
      dr = (rf-ro)/N;   h = dr;   h2 = h*h;    r = ro:h:rf-h;
%
%   matrix eqns. --  interior points
      a = zeros(N,N);   b = zeros(N,1);
      for i = 2:N-1                                                          
        a(i,i-1) = 1 - h/(2*r(i));    
        a(i,i) = -2;                                                     
        a(i,i+1) = 1 + h/(2*r(i));     
        b(i) = -qqq*h2/kf;     
      end                                                                       
%
%   matrix eqns.  --  left boundary  (Case A: Use Forward Approx of order h^2)
      a(1,1) = 2;   a(1,2) = -5;   a(1,3) = 4;   a(1,4) = -1;   
      b(1) = -qqq*h2/(2*kf);    
%
%   matrix eqns.  --  right boundary                                                            
      a(N,N-1) = 1 - h/(2*r(N));     a(N,N) = -2;                                                  
      b(N) = -qqq*h2/kf - (1 + h/(2*r(N)))*Ts;
%
%   solve system of eqns                                                      
      T = a\b;
%
%   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 (N = ',num2str(N),')'])                    
%
%   end of demo
