%
%   PARALLEL_FLOW_1.M  Solve for the Flow Rates in Two Parallel Pipes
%                       (using the Linearized Iteration Method)
%
%   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 
%   the "linearized iteration" approach for solving mildly nonlinear problems.
%   In this method, the coefficient matrix is a function of the solution vector.  
%   However, with a given guess for x, A(x) can be determined, and the solution of
%   a linear system, A(x) * x = b, can be obtained.  This process is continued 
%   until convergence (or we hit the max number of iterations).  Note that a similar
%   solution is available using Matlab's fsolve command in file PARALLEL_FLOW_2.M.
%
%   Note that one often simply lets the next estimate of the solution be the 
%   current best guess, xold = xnew.  However, for this problem, the solution 
%   tends to oscillate around the real solution.  When this happens, we often 
%   use the next guess as the average of the previous two solution estimates,
%   xold = (xnew + xold)/2.  This is the approach taken here.  
%
%   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.
%
%   File prepared by J. R. White, UMass-Lowell (last update:  Nov. 2017)
%
      
      clear all;   close all;   nfig = 0;
%
%   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)
      Qold = [Qtot/2 Qtot/2]';      
%
%   start iteration loop
      itmax = 20;   it = 0;   tol = 1e-5;   emax = 1;
      while emax > tol && it <= itmax
        it = it+1;   
%   
%   set up A and b matrices given current guess and solve for new solution vector
        Re = (D./(nu*A)).*abs(Qold);                      % Reynolds number
        f = 0.25./(log10(eoD/3.7 + 5.74./Re.^0.9)).^2;    % friction factor
        K = (f.*L./D + sumK)./(2*g*A.*A);                 % resistance coefficients
        AA = [1  1;  K(1)*Qold(1)  -K(2)*Qold(2)];        % coefficient matrix
        bb = [Qtot 0]';                                   % RHS vector
        Qnew = AA\bb;                                     % find solution vector
%
%   calc & edit error (intermediate results)
        emax = max(abs(Qnew-Qold)./Qnew);
        fprintf(1,'\n Intermediate edit for Parallel_Flow_1 Problem \n');
        fprintf(1,'  it = %3d      max error = %8.2e \n',it,emax);
        fprintf(1,'      Qnew         Qold     \n');
        for j = 1:length(Qold)
          fprintf(1,' %10.4g   %10.4g    \n',Qnew(j),Qold(j));
        end
%
        Qold = (Qnew+Qold)/2;  % next guess is average of the two previous estimates
      end
%
%   print max relative error and iteration count
      fprintf(1,'\n Number of iterations to convergence = %3d\n',it);
      fprintf(1,' Max relative error at convergence =  %8.2e\n',emax);
      if it >= itmax 
        fprintf(1,' ***** WARNING  --  Hit max number of iterations!!!  *****\n');
      end
%
%   now compute a bunch of parameters to edit (for converged solution)
      Q = Qnew;                                       % simplify the variable names
      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_1:  Summary Results for Simple Parallel Piping System \n');
      fprintf(fid,'                             (using Linearized Iteration Method) \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