%
%   PIPE_FRICTION_1.M   Fluid Flow Application Illustrating Pipe Friction   
%
%   The goal of this problem is to compute and plot the exit velocity associated 
%   with a tank and pipe system.  The problem is highly nonlinear because of the 
%   dependence of the friction factor on the flow velocity.
% 
%   The equation of interest is:
%          V^2 = 2g(h+L)/(1 + fL/D)  
%   where f is the friction factor (given by the Swamee-Jain correlation)
%
%   Here we will use three different pipe sizes and let L vary to determine
%   V vs L.   We expect V to increase with increasing L since the overall head,
%   (h + L), is increasing.  However, with larger L, the friction in the pipe 
%   tends to decrease the flow velocity V.  Thus, there is a trade-off here, 
%   and the goal of this problem is to study this phenomenon.  Let's see if we 
%   can discover this relationship via graphical analysis.  Also, we want to see 
%   how this changes with the pipe diameter, D, where increasing D should 
%   reduce the effect of friction.  Let's give it a try...
%
%   File last modified by J. R. White, UMass-Lowell (last update: Nov. 2017)
%

      clear all, close all, nfig = 0;
      global DD h LL g visc rough
%
%   initial data for problem
      D = [0.0518 0.0874 0.1342];          % different pipe diameters (ft) 
      h = 5;                               % height of fluid in tank (ft)
      g = 32.2;                            % acceleration constant (ft/s^2)
      visc = 1.21e-5;                      % kinematic viscosity of 60 F water (ft^2/s)
      rough = 1.5e-4;                      % steel pipe roughness (ft)
      NL = 50;  L = linspace(1/12,20,NL)'; % setup dependent variable (pipe length)
%
%   let's loop over the number of pipe sizes
      ND = length(D);   V = zeros(NL,ND);   
      for n = 1:ND
        DD = D(n);
%   now use fzero to solve for V for each value of L and D (uses previous V as guess)
        Vo = sqrt(2*g*h);         % initial guess for V for L = 0 ft
        LL = L(1);  V(1,n) = fzero(@pipe_friction_1a,Vo);
        for i = 2:NL
          LL = L(i);  V(i,n) = fzero(@pipe_friction_1a,V(i-1,n));  
        end
      end
%
%   now plot the results 
      nfig = nfig+1;   figure(nfig)
      plot(L,V(:,1),'r-',L,V(:,2),'b--',L,V(:,3),'g:','LineWidth',2)
      title('Pipe\_Friction\_1:  Exit Velocity vs Pipe Length')
      xlabel('Pipe Length (ft)'),ylabel('Exit Flow Velocity (ft/s)'), grid
      legend('Sch 40 1/2" pipe','Sch 40 1" pipe','Sch 40 1-1/2" pipe')
%
%   just for fun, let's calc Re and f for these cases
      Re = zeros(NL,ND);   f = zeros(NL,ND);
      for n = 1:ND
        Re(:,n) =  V(:,n)*D(n)/visc;                             % Reynolds number
        RR = rough/D(n);                                         % relative roughness
        f(:,n) = 0.25./(log10(RR/3.7 + 5.74./Re(:,n).^0.9)).^2;  % friction factor
      end
%
%   now plot these results 
      nfig = nfig+1;   figure(nfig)
      plot(L,Re(:,1),'r-',L,Re(:,2),'b--',L,Re(:,3),'g:','LineWidth',2)
      title('Pipe\_Friction\_1:  Reynolds Number vs Pipe Length')
      xlabel('Pipe Length (ft)'),ylabel('Reynolds Number'), grid
      legend('Sch 40 1/2" pipe','Sch 40 1" pipe','Sch 40 1-1/2" pipe')
%
      nfig = nfig+1;   figure(nfig)
      plot(L,f(:,1),'r-',L,f(:,2),'b--',L,f(:,3),'g:','LineWidth',2)
      title('Pipe\_Friction\_1:  Friction Factor vs Pipe Length')
      xlabel('Pipe Length (ft)'),ylabel('Friction Factor'), grid
      legend('Sch 40 1/2" pipe','Sch 40 1" pipe','Sch 40 1-1/2" pipe')
%
%   end of simulation