%
%   CONDUCTING_ROD_1.M      A Root Finding Example
%      Steady State Temperature vs Current in an Bare Conducting Wire
%
%   This file solves for the steady state temperature in a circular conducting rod 
%   that is heated internally via ohmic heating and that loses energy by both 
%   convective and radiation heat transfer.  The analysis is done for a unit length 
%   of rod or wire.  A range of currents is used so that a plot of Tss vs I can be made.
%
%   The goal of this problem is to show how to solve nonlinear equations using the 
%   built-in fzero routine in Matlab.  It also demonstrates how, with a simple for loop, 
%   to do a parametric study for a single variable -- in this case, we use a range of  
%   values of the current, I, to show how the steady state temperature varies with I.  
%   For each value of I, a nonlinear energy balance equation must be solved (here we use 
%   fzero to do this) to determine the correct steady state temperature.
%
%   A table of results is also produced that shows the various components of the energy 
%   balance eqn for each value of I.
%
%   Note: When including radiation heat transfer, the computations need to be done 
%   using absolute temperatures.  Here, however, the final plot and table display 
%   temperature in C  --  the units conversion is done in the plot and fprintf 
%   commands.
%
%   The basic idea for this problem comes from Example 1.3 on pgs. 17-19
%   in the text, Fundamentals of Heat and Mass Transfer, by Incropera and DeWitt 
%   5th Ed (2002, John Wiley & Sons). 
%
%   File prepared by J. R. White, UMass-Lowell (last update:  Nov. 2017)
%
      
      clear all;   close all;   nfig = 0;
      global D emiss R_L sig h Tinf Tsur Isq
%
%   set parameters for the problem
      D = 0.001;            % wire diameter (m)
      emiss = 0.8;          % emissivity (dimensionless)
      R_L = 0.4;            % resistance per unit length (ohm/m)
      sig = 5.67e-8;        % Stefan-Boltzmann constant (W/m^2-K^4)
      h = 100;              % heat transfer coeff (W/m^2-K)
      Tinf = 25+273.15;     % bulk fluid temperature (K)
      Tsur = Tinf;          % temperature of surroundings for radiation HT (K)
%
      NI = 21;  I = linspace(0,10,NI); % set range for wire current (amps)
      Tss = zeros(size(I));            % allocate space for storage of SS temperatures
      Isq = I(1)*I(1);                 % I^2 for first value of current
      Tg = Tinf;                       % initial guess of Tss for first current
      Tss(1) = fzero(@conducting_rod_1a,Tg);  % find actual Tss for first value of I
%   Note: The above step was not really necessary for this case, since for I = 0
%         and Tsur = Tinf, Tss = Tinf (so we could simply set Tss(1) = Tinf.  However,
%         for the general case, this is not true...
%
%   find SS temperature for each current (implicit)
      for n = 2:NI
        Isq = I(n)*I(n);
        Tss(n) = fzero(@conducting_rod_1a,Tss(n-1));
      end
%
%   find SS temperature for each current assuming no radiation heat transfer (explicit)
       Tssa = Tinf + I.*I*R_L/(h*pi*D);  % this is often a poor approx (not so here)
%
%   plot Tss vs I
      nfig = nfig+1;   figure(nfig)
      plot(I,Tss-273.15,'r-',I,Tssa-273.15,'b--','LineWidth',2), grid on
      rr = axis;   rr(3) = 0;   axis(rr);
      title('Steady State Wire Temperature vs Current')
      xlabel('Current (amps)'),ylabel('Steady State Temperature (^oC)')
      legend('T_{ss} with radiation','T_{ss} with no radiation', ...
             'Location','NorthWest')
%
%   evaluate components of balance equation (note 'dot arithmetic' where appropriate)
      A = pi*D;                              % surface area per unit length
      gen = I.*I*R_L;                        % energy generation rate
      conv = h*A*(Tss - Tinf);               % loss rate by convection
      rad = emiss*sig*A*(Tss.^4 - Tsur^4);   % loss rate via radiation
      bal =  (gen - conv - rad);             % SS energy balance eqn 
%
%   write summary table of results 
      fprintf(1,'\n\n');
      fprintf(1,'           Summary Results from the CONDUCTING_ROD_1.M Program \n');
      fprintf(1,'\n');
      fprintf(1,'    Basic Design Parameters      \n');
      fprintf(1,'      Wire Diameter (m):                      %8.3f    \n',D);
      fprintf(1,'      Emissivity (dimensionless):             %8.2f    \n',emiss);
      fprintf(1,'      Resistance per Unit Length (ohm/m):     %8.2f    \n',R_L);
      fprintf(1,'      Stefan-Boltzmann Constant (W/m^2-K^4): %9.2e    \n',sig);
      fprintf(1,'      Heat Transfer Coeff (W/m^2-K):          %8.2f    \n',h);
      fprintf(1,'      Bulk Fluid Temperature (C):             %8.2f    \n',Tinf-273.15);
      fprintf(1,'      Temperature of Surroundings (C):        %8.2f    \n',Tsur-273.15);
      fprintf(1,'\n');
      fprintf(1,'    Components of Energy Balance vs. Current  \n')
      fprintf(1,'        I    Generation    Convection   Radiation     Balance      SS Temp \n');
      fprintf(1,'       (A)     (W/m)         (W/m)        (W/m)        (W/m)         (C)   \n');
      for n = 1:NI
        fprintf(1,'     %5.1f  %8.2f      %8.2f     %8.2f     %8.3f     %8.2f    \n', ...
                  I(n),gen(n),conv(n),rad(n),bal(n),Tss(n)-273.15);
      end
%
%   end of problem