%
%   PIPE_HeadAdded  Type III Problem -- Solve for Pipe Diameter to achieve
%                a desired flow withthe given pump power (from Class Example)
%
%   This file computes the required head added to achieve a desired flow rate as
%   a function of the pipe diameter.  Then a plot of Head Added versus D is made.  
%   With the Pump Power (ie. actual head added) given in the problem statement, one 
%   can determine the smallest diameter pipe that can achieve this in the given system.
%  
%   All the geometry information and fluid properties must be identified.  
%   Depending on flow regime, either the laminar or turbulent flow relationship 
%   is used to find the friction factor.  There are no Minor losses for this problem.
%   The total head loss is found from Darcy's Eqn.  Finally, the energy equation is 
%   then used to compute the head added for the various pipe diameters.
%
%   File prepared by J. R. White, UMass-Lowell (last update:  April 2017)
%
      
      clear all;   close all;   nfig = 0;
%
%   set the problem parameters (from Class Example)
      g = 9.81;                       % gravitational constant (m/s^2)
      L = 500;                        % channel length (m)
      rho = 1000;                     % fluid density (kg/m^3)
      mu = 1.15e-3;                   % fluid viscosity (N-s/m^2)
      e = 0.00002;                    % surface roughness (flexible tubing)(m)
      Q = 0.04;                       % volumetric flow rates (m^3/s)
      actPA = 40e3;                   % actual power added by pump (W)
      acthA = actPA/(rho*g*Q);        % actual head added by pump (m)
%
%   set pipe diameters and related quantities 
      D =  0.075:0.005:0.15;          % set range on pipe diameters (m)
      ND = length(D);                 % number of pipe diameters
      A =  pi*D.*D/4;                 % flow area (m^2)
      eoD = e./D;                     % relative roughness
%
%   compute variety of parameters
      v = Q./A;                       % ave fluid velocity (m/s)
      Re = rho*v.*D/mu;               % Reynolds number
      f = zeros(size(D));             % initialize and calc the friction factor
      for i = 1:ND
        if Re(i) <= 2300,  f(i) = 64/Re(i);  end
        if Re(i) > 2300 && Re(i) < 4000
          disp('  *** Warning: Flow in Transition Region ***');  f(i) = 0.05;  end
        if Re(i) >= 4000, f(i) = 0.25/(log10(eoD(i)/3.7 + 5.74/Re(i)^0.9))^2;  end
      end
      vhead = v.*v/(2*g);             % dynamic or velocity head (m)
      hLmajor = f.*(L./D).*vhead;     % major head loss (m)
      hLminor = 0;                    % minor head loss (m)  (none for this problem)
%
%   now evaluate the energy equation
      hA = vhead + hLmajor + hLminor; % required head added to drive flow
%
%   interpolate to find the smallest D
      Dmin = interp1(hA,D,acthA);       
%
%   let's plot hA versus D
      nfig = nfig+1;   figure(nfig)
      plot(D*100,hA,'r-','LineWidth',2),grid,hold on
      r = axis;
      plot(r(1:2),[acthA acthA],'b-','LineWidth',2),hold off
      title('Head Added versus Pipe Diameter (Type III Class Example)')
      xlabel('Diameter (cm)')
      ylabel('Required Head Added (m)')
      gtext(['Smallest allowed D = ',num2str(Dmin*100,'%6.2f'),' cm'])
%
%   print results    
      fid = 1;
      fprintf(fid,'  \n');
      fprintf(fid,'   Calculated Parameters: \n');
      fprintf(fid,'     diameter     ave vel   Reynolds #  fric factor   vel head    hL major    hL minor       hA \n');
      fprintf(fid,'       (m)         (m/s)                                (m)         (m)         (m)          (m) \n');
      for i = 1:ND
        fprintf(fid,'    %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e \n', ...
                    D(i),v(i),Re(i),f(i),vhead(i),hLmajor(i),hLminor,hA(i));
      end
      fprintf(fid,'  \n');
      fprintf(fid,'   Interpolated Result for Pump Power of %4.1f kW:  Dmin = %10.3e cm\n',actPA/1000,Dmin*100);
%
%   end of file