3907.sagews

"""
Implementation of the Hardy-Cross relaxation solution of the pipeline network, 
given in the following example:
     (2)
A------------B
|          / |
|         /  |
|        /   |
|(1) (3)/ (4)|
|      /     |
|     /      |
|    /       |
|   /        |
|  /         |
| /    (5)   |
C------------D

QA = 500 # m**3/h
QB = 0   # m**3/h
QC = 0   # m**3/h
QD = -500 # m**3/h

R1 = 0.000298
R2 = 0.000939
R3 = 0.00695
R4 = 0.000349
R5 = 0.000298

"""

# import all the numerical library, replaces Matlab
from numpy import * 


# Initial conditions:
QA = 500 # m**3/h
QB = 0   # m**3/h
QC = 0   # m**3/h
QD = -500 # m**3/h

# array is like a vector in Matlab
L = array([2,3,3,2,3],dtype='f') # Length in km, 'f' means 'float'
D = array([25,20,20,15,20],dtype='f') # Diameter, cm
C = array([100,110,120,130,120],dtype='f') # Hazen Williams friction coefficient
 

# resistance as a single line function that receives C, D, L and returns R
r = lambda L,D,C: L*1.526e7/(C**1.852*D**4.87)

R = r(L,D,C)
print "Resistances: ", R

# Multiline definition of a function
# starts with def name of the function and inputs in parentheses
# next line is with indentation: 4 spaces or a TAB
# ends with the "return" and the list of outputs
# see the following example

#def resistance(L,d,c):
#    r = 1.526e7/(c**1.852*d**4.87)*L
#    return r           


# head loss as a function
# note that one cannot raise a negative value to the power of 1.852

hf = lambda R,Q: R*sign(Q)*power(abs(Q),1.852)



# define branches - each row contains numbers of pipes:
branch = array([[2,3,1],[3,4,5]])-1 # Python counts from zero, not 1, the first pipe will be (0)

rows,cols = branch.shape # rows = num of branches, cols = pipes in each


# initial guess that 
Q = array([-250, 250.0, 0.0, 250, -250.0]) # m**3/hr

dQ = 1.0

# main loop
while abs(dQ) > 0.5:
    # arange(N) gives a vector from 0 to N-1, i.e. arange(3) = 0,1,2
	for i in arange(rows):
                # estimate the losses in each pipe
		y = hf(R,Q) 
		

		yq = abs(1.852*y/abs(Q))
                
                # Remove NaNs for the exceptional cases
		yq[isnan(yq)] = 0.0
		
		# for the first branch
		# Sum is a sum :)
		sumyq = sum(yq[branch[i,:]])
		sumy = sum(y[branch[i,:]])
		
                # Estimate the correction dQ for the next step:
		dQ = -1*sumy/sumyq

		print("dQ = %f" % dQ)
		
                # += is equal to Q = Q + dQ
		Q[branch[i,:]] += dQ
		

print("Discharges [m^3/hr]:\n")
print Q

# we're looking for equivalent pipe solution with:
# Deq = 25 # cm
# Ceq = 120 # 

# y = hf(R[1],Q[1])+hf(R[2],Q[2])

# Leq = y/(r(1,25,120)*Q[0]**1.852)

# Leq = Leq + L[0] + L[6]*(120/C[6])**1.852

# print("Equivalent length = %3.2f km" % Leq)

Resistances: [ 0.00093889 0.00349946 0.00297863 0.0069502 0.00297863] Warning: invalid value encountered in divide dQ = -77.877620 dQ = -44.395536 dQ = 14.902768 dQ = -2.315951 dQ = 0.672938 dQ = -0.092987 Discharges [m^3/hr]: [-312.30191449 187.69808551 -109.10638868 203.19552581 -296.80447419]