% Water_Hammer_Streeter_example_13_9
% 19.05.2019


 % Water hammer example 13.9 from Streeter [1], pages 548 to 551 

 % Water hammer program, with original source code in Basic (IBM-PC) 
 % System is shown in Figure 13.23 
 %    |      |
 %    |      |
 %    |      |
 %    |      |
 %    |______=================================X   
 %   Reservoir            Pipe              Valve
 
  %                Figure 13.23
  
 % The system consists of a reservoir with an initial water level "Hres", 
 % a horizontal pipe of length "L", inside diameter "D", and friction 
 % factor "f" (constant value assumed). At the end of the pipe is located 
 % a valve  with a -Product of discharge coefficient and opening area-   
 % defined by the matrix "CVA". 			

 % Code Basic (IBM-PC) is shown  in  figure 13.24 
 % Results table in  figure 13.25 
 % The system of Figure 13.23 has an initial valve opening 
 % "Cd * A" = 0.06 m² (at time T = 0)                         '
 % At time intervals of CVD = 5 [s], the valve opening "Cd * A"  takes 
 % the values:  % 0.03, 0.01, 0.003, 0.001, 0.0005, 0.0002, 0 and then 
 % remains closed.                                 '
 % Determine system transients until a time Tmax = 40 [s], after the valve 
 % begins  to close.                     '
 % Example uses a pipe division of 4 segments (N = 4). 
 % Wave velocity (celerity) A = 1200 (m/s)          
 % Values for CVA are given at time intervals DCV = 5 [s]  
 
 % Calculated values are stored in the two dimensional matrices Q and H
 %                   T-range:  0  to  Tmax
 % Second dimension  I: Index of the node position I = 1  to I = NS
 % Some row numbers from Streeter code are maintained, to facilitate the
 % comparison.
 % Equation numbers of the type (13.6.8) corresponds to the Streeter
 % equations
 % Equation numbers of the type (E-15), corresponds to deducted equations
 % Deduction of equations and its boundary vales are shown in 
 % www.piping-rools.net in the file:
 % Water hammer. Method of characteristics. Equations, Boundary conditions 
 % and Initial conditions
 
 
clear variables;
close all;
clc;

% Constants		
g = 9.80665;    % [m/s²]		
		
% Counter		
num = 0;		

% ..................... Input data .................................
% 1700.  	
		
% 1.- Reservoir height		
Hres =	100;	% m
% 2.- Pipe length		
L =	4800;	    %m
% 3.- Pipe internal diameter		
D =	2;	        %m
% 4.- Friction factor		
f = 	0.022;	% -
% 5.- Celerity		
a =	1200;	    % m/s
% 6.- Number of sections		
N =	4;	        % -
% 7.- Time interval of CVA values		
DCV =	5;   	% s
% 8.- Maximum elapsed time		
Tmax =	48;	    % s   

% Matrix CVA. Product of discharge coefficient and opening area
% Inputs 9 to 19

CVA(1) =	0.06;      % At time = 0 s
CVA(2) = 	0.03;      % At time = 5 s
CVA(3) = 	0.01;      % At time = 10 s
CVA(4) = 	0.003;     % ...
CVA(5) = 	0.001;
CVA(6) = 	0.0005;
CVA(7) = 	0.0002;
CVA(8) =	0;
CVA(9) = 	0;
CVA(11) =	0;


%........................... Help variables ...........................
% 2600:                                                                                                   
% Area of pipe section                                                                                 
 AR = (pi / 4) * D ^ 2;                                                                                  
% Constant for Eq (13.6.9)                                                                              
B = a / (g * AR);                                                           % Eq. (13.6.8a)     
                                                                                                         
% Last node "NS", corresponds to the valve position                                                    
NS = N + 1 ;                                                                                            
                                                                                                        
% Segment lenght                                                                                        
Dx = L / N; 

% Time increment    DT = DX / a                                             % Eq. (13.6.2)    
DT = L / (a * N);

% .........Two dimensional matrices  H and Q with zero-values .........
% 2620
H = zeros(Tmax,NS);                   % H matrix Tmax x 5, with zeros
Q = zeros(Tmax,NS);                   % Q matrix Tmax x 5, with zeros

% ...................... Initial values ...............................		
% 2700:
% Initial  time "Ti" 		
Ti = 0 ;
% To use the time as index, since the index cannot be zero, use
T = Ti+1;                       % Thus, for the time = 0, T = 1

% Valve CV value at t = 0 (index 1), from input matrix CVA        		
CV0 = CVA(1);

% 2720
% Asignment to the matrix element H(1,1), the reservoir pressure value
% at time = 0 (index 1), at position node 1
H(1,1) = Hres;                                                              % Eq. (E-20)

% 2740
% Initial flow rate		
Q0 = (Hres / (f * L / (2 * g * D * AR ^ 2) + 1 / (2 * g * CV0 ^ 2)))^0.5 ;  % Eq. (E-15) 
Q(1,1) = Q0;
% at time = 0 (index 1) and node 1
 		
% Initial piezometric height at the valve  Hvalve_0   
% at time = 0 (index 1)  and node  N = NS  		
Hvalve_0 = Q0 ^ 2 / (2 * g * CV0 ^ 2);                                      % Eq (E-13)
		
% Resistance coefficient "R" per segment  		
R = (Hres - Hvalve_0) / (Q0 ^ 2 * N);                                       % Eq (E-21)


% ........... Initial values  of Q and H  .............................
% at time =0 (index 1)
% Initial values are calculated for all the nodes  ( 1  to  NS)
% 2900
for I = 1:NS   		
    Q(1,I) = Q0;  		                                                    % Eq (E-15)
    H(1,I) = Hres - (I - 1) * R * Q0 ^ 2;                                   % Eq (E-19)
% Q and H values checked
end


% ......................Start of new calculation....................
% 3900
% K value for interpolation of CV Integer.      
% Kmax: maximum K-value, corresponding to the time Tmax

Kmax = floor(Tmax / DCV);                    
% B = floor(A) rounds the elements of A to 
% the nearest integer less than or equal to A   
                                             
% 3901                                    
while T< Tmax	                            

 num = num +1;
 T = T + DT;
 
Tant = T-1; 
% where T is the time of new calculated term
% and Tant is the previous time ( corresponding to terms already known)
% K value for interpolation of CV Integer.

K = floor(T / DCV) + 1 ; 		                                            % Eq. (H-1)


% ..... CV value is interpolated in each group of DCV seconds .....
% It corresponds to the previous time Tant.
% 4100
CV = CVA(K) + (Tant - (K - 1) * DCV) * (CVA(K + 1) - CVA(K)) / DCV ;        % Eq. (H-2)


% ......... Boundary condition at the valve (node NS).............		
% 4200
% Flow rate and piezometric height at the valve (node NS).  
% Previous node: N												 	

% CP and BP values at (T,NS)
CP = H(Tant,N)+ B*Q(Tant,N);                                                % Eq. (G-5) Eq(13.6.11a
BP = B + R * abs(Q(Tant,N));                                                % Eq. (G-6) Eq(13.6.11b) 

% Q(T,NS) and H(T,NS)
Q(T,NS) = -g * CV^2 * BP + ((g * BP * CV^2)^2+ CV^2 *2*g *CP)^0.5 ;         % Eq. (G-4)
H(T,NS) = CP - BP * Q(T,NS) ;                                               % Eq. (G-3)


%......... Boundary condition at the reservoir .....................
% 4600                                                                     
% The value of H at the reservoir (node 1), at any time, 
% is already defined as                                   
H(T,1) = Hres;                                                              % Eq.(E-20)
% Flow rate at the reservoir at a time T
Q(T,1) = (H(Tant,1)-H(Tant,2)+ B * Q(Tant,2)) / (B + R * abs(Q(Tant,2)));   % Eq.(F-4)b


%........ Interior sections .......................................
% 4800
% Whewn calculating a term of a node I, the term to the link of it
% has and index "I-1" and it is denoted as  Ilink
% Similarly, the term of the right of it has an index "I+1" and it is 
% denoted as  Iright

for I = 2:N      %  I = 2:N
    Ilink = I-1;
    Iright = I+1;
    CP = H(Tant,Ilink) + B * Q(Tant,Ilink);                                 % Eq(13.6.11a) 
    BP = B + R * abs(Q(Tant,Ilink));                                        % Eq(13.6.11b)
    CM = H(Tant,Iright) - B * Q(Tant,Iright);                               % Eq(13.6.12a)
    BM = B + R * abs(Q(Tant,Iright));                                       % Eq(13.6.12a)
    Q(T,I) = (CP - CM) / (BP + BM);                                         % Eq(13.6.15) 
    H(T,I) = CP - BP * Q(T,I);                                              % Eq(13.6.16)
end 

end    % goto 3901 to continue the while loop
       % or the loop will exit when T > Tmax

 
       
% ...................  Printing of results  ...........................

% Matrix of the flows  Q in Command Window
% Q
% Matrix of the Piezometric heights in Command Window
% H
% Maximum element of the matrix H
 [max_H_value] = max(H(:));
% This value has to be compared with the value obtained when using the
% case of the total pressure developed by a sudden  shutoff.		
% (Joukovsky / Michaud equations)		
% P_Joukovsky = 198.6	m.w.c.
% In this Streeter example, Hmax = 163.0326 (= 82% of P_Joukovsky)
% Location of the maximum value of  H
 [row, col] = find(ismember(H, max(H(:))));
% Value of the element (10,5) ot the H matrix
  H(10,5);
  
% Linearly space vector. From 1 to Tmax, with a number of elements: Tmax
 T = linspace(1,Tmax,Tmax);
 % __________________________________________
 % Column with piezometric heights of node 5
 % H(T,5)
 %___________________________________________
 % Column with flow rate of node 3
 % Q(T,3)
 %___________________________________________
 % Matrix with piezometric heights of all nodes
 plot(T,H(T,1),'r', T,H(T,2),'m',T,H(T,3),'c', T,H(T,4),'b', T,H(T,5),'g')
 title('Piezometric height' );
 xlabel('t');			
 ylabel('H' );
 legend('I=1. Valve','I=2','I=3','I=4','I=5.Resrvoir' );
 %___________________________________________
 % Plot of flow rates of all nodes
 % plot(T,Q(T,1), T,Q(T,2), T,Q(T,3), T,Q(T,4), T,Q(T,5))
 %___________________________________________  
 
 % Program end
 
 % References
 % [1] Mecánica de fluidos
 % Victor L. Streeter and E. Benjamin Wylie
 % 8th edition, 3rd in spanish























