%
%   PARALLEL_FLOW_2.M  Solve for the Flow Rates in Two Parallel Pipes
%                       (using Matlab's fsolve command)
%
%   This file computes the flow rates in a relatively simple parallel piping system.
%   This problem is derived from Prob. 12.1 in the Fluid Mechanics text by Mott 
%   (5th Ed.).  The problem reduces to a coupled set of 2 nonlinear algebraic 
%   equations (for steady flow).  An iterative hand calculation was done in my 
%   Fluid Mechanics class (many years ago).
%   
%   This file simply automates the solution of the two nonlinear equations using 
%   Matlab's fsolve command.  Note that a similar solution is available using 
%   the Linearized Iteration Method in file PARALLEL_FLOW_1.M.
%
%   Also note that the Reynolds number must be positive.  Since Q might take on  
%   negative values (flow direction is opposite to the assumed direction) during  
%   the iterative process, we will use abs(Q) in the computation of Re.
%
%   The function file PARALLEL_FLOW_2A.M evaluates the two nonlinear equations for
%   any given flow distribution, Q, for use by the fsolve command.  The goal here is
%   to find Q such that F(Q) = 0, where Q is the vector of flow rates.
%
%   File prepared by J. R. White, UMass-Lowell (last update:  Nov. 2017)
%
      
      clear all;   close all;   nfig = 0;
      global D A L eoD sumK nu g Qtot
%
%   set the geometry parameters 
      g = 9.8;                               % gravitational constant (m/s^2)
      D =  [0.1023 0.0779]';                 % pipe diameter (4" and 3" Sch 40) (m)
      A =  pi*D.*D/4;                        % flow area (m^2)
      L = [60 60]';                          % pipe lengths (m)
      e = [4.6e-5 4.6e-5]';                  % surface roughness (m)
      eoD = e./D;                            % relative roughness
      fT = 0.25./(log10(eoD/3.7)).^2;        % Darcy friction factor (fully turbulent)
      sumK = [2*30*fT(1) (2*30+240)*fT(2)]'; % minor loss factors
%   minor losses: pipe 1 -> 2 standard elbows
%                 pipe 2 -> 2 standard elbows + valve
%
%   set fluid properties  (oil)
      gam = 8.80e+3;                  % fluid specific weight (N/m^3)
      nu = 4.8e-6;                    % fluid kinematic viscosity (m^2/s)
      Qtot = 0.06;                    % total flow rate in common pipe (m^3/s)
%
%   set reasonable initial guess for flow rates  (m^3/s)
      Qg = [Qtot/2 Qtot/2]';      
%
%   call fsolve to find solution to system of unconstrained nonlinear eqns
      Q = fsolve(@parallel_flow_2a,Qg);
%
%   now compute a bunch of parameters to edit (for converged solution)
      v = Q./A;                                       % average flow velocity
      Re = abs(v).*D/nu;                              % Reynolds number
      f = 0.25./(log10(eoD/3.7 + 5.74./Re.^0.9)).^2;  % Darcy friction factor
      K = (f.*(L./D) + sumK)./(2*g*A.*A);             % resistance coeff
      hL = K.*Q.^2;                                   % head loss
%
%   print summary results
      fid = 1;
      fprintf(fid,' \n');
      fprintf(fid,'     Parallel_Flow_2:  Summary Results for Simple Parallel Piping System \n');
      fprintf(fid,'                             (using Matlab''s fsolve command) \n');
      fprintf(fid,' \n');
      fprintf(fid,'       Geometry (Sch 40 steel pipe): \n');
      fprintf(fid,'     Pipe#   Diameter  Length   Flow Area  Surface Roughness  Rel Roughness  Minor Losses\n');
      fprintf(fid,'               (m)      (m)       (m^2)          (m)                             SumK \n');
      for n = 1:length(Q)
        fprintf(fid,' %7i     %7.4f %7.1f   %10.2e   %10.2e        %10.2e      %7.2f \n', ...
                    n,D(n),L(n),A(n),e(n),eoD(n),sumK(n));
      end
      fprintf(fid,' \n');
      fprintf(fid,'       Fluid Properties: \n');
      fprintf(fid,'         sp weight (N/m^3):            %10.3e \n',gam);
      fprintf(fid,'         kinematic viscosity ((m^2/s): %10.3e \n',nu);
      fprintf(fid,'         total flow rate (m^3/s):      %10.3f \n',Qtot);
      fprintf(fid,' \n');
      fprintf(fid,'       Calculated Parameters for each Pipe Segment: \n');
      fprintf(fid,'     Pipe#   Length   Flow Rate  Ave Vel  Reynolds#  FricFactor   Kcoeff   Head Loss \n');
      fprintf(fid,'               (m)     (m^3/s)    (m/s)                         (s^2/m^5)     (m) \n');
      for n = 1:length(Q)
        fprintf(fid,' %7i   %7.1f     %7.4f  %7.2f  %10.2e   %7.4f  %10.2e   %7.3f \n', ...
                    n,L(n),Q(n),v(n),Re(n),f(n),K(n),hL(n));
      end
      fprintf(fid,' \n');
      fprintf(fid,'       Total Calculated Flow Rate (m^3/s):              %7.4f \n',sum(Q));
      fprintf(fid,'       Calculated DeltaP across Branched Section (kPa): %7.2f \n',gam*hL(1)/1000);
%
%   end of file