%
%   Fuel_Temp_SH2.M   Solve for T(r) in a cylindrical fuel pin 
%               Analytically and via the Shooting Method for variable kf
%
%   Given differential eqn.:
%         T'' + (1/r + (1/kf)kf'T')T' + qqq/kf =  0
%   with  T'(0) = 0  and   T(Rf) = Ts 
%
%
%   File prepared by J. R. White, UMass-Lowell (Oct. 2003)
%

%
%   getting started
      clear all, close all,  nfig = 0;
%
%   set properties
      global qqq a1 a2
      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;  % coeffs for correlation of kf(T)
%
%   Analytical Solution (exact)
%
      re = linspace(0,Rf,30);   c = 1+a2;
      Te = (Ts^c + (c/a1)*(qqq/4)*(Rf*Rf - re.*re)).^(1/c);
%
%   Solution via Shooting Method (numerical approx)
%
%   set coefficients for BCs and domain for independent variable
      xo = 0;  xf = Rf;
      zbc = [0 1.0 0.0;      % left BC  --> T'(xo) = 0.0
             1 0 Ts];        % right BC --> T(xf) = Ts
%
%   solution via Shooting Method (numerical approx)  [use Ts as a rough guess for alpha]
      tol = 1e-6;  options = odeset('RelTol',tol);   % set tight tolerance for ODE soln
      [xs,zs] = bvp2sh('fuel_temp_sh2_eqns',[xo xf],zbc,options,Ts);
%
%   plot solution 
      nfig = nfig+1;   figure(nfig)
      plot(re,Te,'r-',xs,zs(:,1),'bo','LineWidth',2),grid
      title('Fuel Pin Radial Temperature Profile (Variable kf)')                  
      xlabel('Radial Distance (ft)'),ylabel('Temperature (F)')    
      legend('Analytical Solution','Shooting Method (SH2)')                    
%
%   end program
%
