%
%   Fuel_Temp_Exact.M   Plot the analytical T(r) in a cylindrical fuel 
%                         pin with constant and variable kf
%
%   This file simply compares the analytical solution for the variable 
%   thermal conductivity case versus the constant kf case.  The heat 
%   conduction equations was solved in class for the two cases and 
%   the resultant temperature profiles are simply plotted below.
%
%   File prepared by J. R. White, UMass-Lowell (Oct. 2004)
%

%
%   getting started
      clear all, close all,  nfig = 0;
      global Ts qqq Rf a1 a2
%
%   set properties
      qqq = 1.35e+7;            % average volumetric heat source (BTU/hr-ft^3)
      kfo = 1.1;                % first choice for kf for UO2 (BTU/hr-ft-F)
      Rf = 0.0175;              % fuel pin radius (ft)
      Ts = 650;                 % surface temperature (F)
      a1 = 94.4;   a2 = -0.556; % constants for kf = a1*T^a2
      re = linspace(0,Rf,30);   % set r-vector
%
%   Analytical Solution for constant kf (initial guess)
      Tco = Ts + (qqq/(4*kfo))*(Rf*Rf - re.*re);     
%
%   Analytical Solution for variable kf  
      c = 1+a2;  Tv = (Ts^c + (c/a1)*(qqq/4)*(Rf*Rf - re.*re)).^(1/c);
%
%   Analytical Solution for constant kf (better guess -- use Tave)
      Tint = quadl('fuel_temp_exact_int',0,Rf);
      V = pi*Rf*Rf;   Tave = Tint/V;
      kfb = a1*Tave.^a2;;      % better choice for kf
      Tcb = Ts + (qqq/(4*kfb))*(Rf*Rf - re.*re);     
%
%   plot kf(T) 
      T = linspace(500,4000,50);  kf = a1*T.^a2;
      nfig = nfig+1;   figure(nfig)
      plot([500 4000],[kfo kfo],'b--',[500 4000],[kfb kfb],'g-.', ...
            T,kf,'r-','LineWidth',2),grid
      title('UO2 Thermal Conductivity versus Temperature')                  
      xlabel('Temperature (F)'),ylabel('Thermal Conductivity (BTU/hr-ft-F)') 
      legend(['Initial kf = ',num2str(kfo,'%4.2f'),'  BTU/hr-ft-F'], ...
             ['kf(Tave) = ',num2str(kfb,'%4.2f'),'  BTU/hr-ft-F'],'Variable kf')                    
%
%   plot temperature profiles 
      nfig = nfig+1;   figure(nfig)
      plot(re,Tco,'b--',re,Tcb,'g-.',re,Tv,'r-','LineWidth',2),grid
      title('Fuel Pin Radial Temperature Profiles')                  
      xlabel('Radial Distance (ft)'),ylabel('Temperature (F)')    
      legend(['Initial kf = ',num2str(kfo,'%4.2f'),' BTU/hr-ft-F'], ...
             ['kf(Tave) = ',num2str(kfb,'%4.2f'),' BTU/hr-ft-F'], ...
             ['Variable kf -- Tave = ',num2str(Tave,'%5.1f'),' F'])  
%
%   end main program