%
%   PIPE_PHead  Type II Problem -- Solve for Flow Rate given the Pressure Head 
%                               (from Class Example)
%
%   This file computes the pressure head associated with a vector of flow rates.  
%   Then a plot of Pressure Head versus Flow Rate is made.  With the desired 
%   pressure given in the problem statement, one can determine the actual flow 
%   rate 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.  Minor losses need to be included 
%   explicitly.  The total head loss is found from Darcy's Eqn.  Finally, the 
%   energy equation is then used to compute the pressure head for several flow 
%   rates.  
%
%   File prepared by J. R. White, UMass-Lowell (last update: April 2017)
%
      
      clear all;   close all;   nfig = 0;
%
%   set the geometry parameters (from Class Example)
      g = 32.2;                       % gravitational constant (ft/s^2)
      D =  2/12;                      % pipe diameter (ft)
      A =  pi*D*D/4;                  % flow area (ft^2)
      L = 32;                         % channel length (ft)
      e = 0.0005;                     % surface roughness (galvanized iron)(ft)
      eoD = e/D;                      % relative roughness
      ehead = 16-4;                   % elevation head (z2-z1)  (ft)
%
%   set fluid properties and volume flow rate (water at 70 F)
      rho = 1.94;                     % fluid density (slugs/ft^3)
      mu = 20.2e-6;                   % fluid viscosity (lbf-s/ft^2)
      Qgpm = 100:10:300;              % vector of volumetric flow rates (gpm)
      Q = Qgpm/448.83;                % vector of volumetric flow rates (ft^3/s)
      NQ = length(Q);                 % number of flow rates
%
%   treat minor losses (see Class Example)
      K = 0.5 + 2*0.9 + 0.19;         % total resistance coeff 
%   minor losses: sudden contraction + elbows + fully open globe value 
%
%   compute variety of parameters
      v = Q/A;                        % ave fluid velocity (ft/s)
      Re = rho*v*D/mu;                % Reynolds number
      f = zeros(size(Q));             % initialize and calc the friction factor
      for i = 1:NQ
        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/3.7 + 5.74/Re(i)^0.9))^2;  end
      end
      vhead = v.*v/(2*g);             % dynamic or velocity head (ft)
      hLmajor = f*(L/D).*vhead;       % major head loss (ft)
      hLminor = K*vhead;              % minor head loss (ft)
%
%   now evaluate the energy equation
      phead = ehead + vhead + hLmajor + hLminor; % desired elevation head
      P1 = 40;                        % pressure in tank (psi)
      desired_phead = P1*144/(rho*g); % desired pressure head
%
%   interpolate to find the desired Q
      Qfinal = interp1(phead,Q,desired_phead)*448.83;       
%
%   let's plot phead versus Qgpm
      nfig = nfig+1;   figure(nfig)
      plot(Qgpm,phead,'r-','LineWidth',2),grid,hold on
      r = axis;
      plot(r(1:2),[desired_phead desired_phead],'b-','LineWidth',2),hold off
      title('Pressure Head versus Flow Rate (Type II Class Example)')
      xlabel('Flow Rate (gal/min)')
      ylabel('Pressure Head (ft)')
      gtext(['Desired Q = ',num2str(Qfinal,'%10.3e'),' gpm'])
%
%   print results
      fid = 1;
      fprintf(fid,'   Results for Type II Class Example  \n');
      fprintf(fid,'     Geometry: \n');
      fprintf(fid,'       pipe dia (ft):           %10.3e \n',D);
      fprintf(fid,'       flow area (ft^2):        %10.3e \n',A);
      fprintf(fid,'       channel length (ft):     %10.3e \n',L);
      fprintf(fid,'       surface roughness (ft):  %10.3e \n',e);
      fprintf(fid,'       relative roughness :     %10.3e \n',eoD);
      fprintf(fid,'       elevations head (ft):    %10.3e \n',ehead);
      fprintf(fid,'     Fluid Properties: \n');
      fprintf(fid,'       density (slugs/ft^3):    %10.3e \n',rho);
      fprintf(fid,'       viscosity (lbf-s/ft^2):  %10.3e \n',mu);
      fprintf(fid,'    \n');
      fprintf(fid,'     Calculated Parameters: \n');
      fprintf(fid,'       flow rate    ave vel   Reynolds #  fric factor  vel head    hL major    hL minor   press head \n');
      fprintf(fid,'         (gpm)      (ft/s)                               (ft)        (ft)        (ft)        (ft) \n');
      for i = 1:NQ
        fprintf(fid,'      %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e \n', ...
                    Qgpm(i),v(i),Re(i),f(i),vhead(i),hLmajor(i),hLminor(i),phead(i));
      end
      fprintf(fid,'    \n');
      fprintf(fid,'     Interpolated Result for Tank Pressure of %4.1f psi:  Q = %10.3e gpm\n', ...
                  P1,Qfinal);
%
%   end of file