3921.sagews

"""
Demostration of the Hardy-Cross relaxation solution of the pipeline network
 (1)  (2)   (3)    (7)
----o-----o-----o------
     \(4) |(5) /(6)
      \   |   /
       \--o--/
       
       
"""
import numpy as n
import pylab as p

L = n.array([2,2,3,3,2,3,2],dtype='f') # km
D = n.array([25,25,20,20,15,20,25],dtype='f') # cm
C = n.array([120,100,110,120,130,120,130],dtype='f') # Hazen Williams 

# resistance per unit length as a function
r = lambda l,d,c: 1.526e7/(c**1.852*d**4.87)*l
# head loss as a function:
hf = lambda R,Q: R*n.sign(Q)*n.abs(Q)**1.852

R = r(L,D,C)# resistances

# define branches - each row contains numbers of pipes:
branch = n.array([[2,5,4],[3,6,5]])-1 # Python counts from zero, not 1

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


# initial guess that 
Q = n.array([500,250,250,-250,0,-250,500],dtype='f') # m**3/hr

dQ = 500.0

# main loop
while abs(dQ) > 0.5:
	for i in n.arange(rows):
		y = hf(R,Q) 
		
		yq = n.abs(1.852*y/n.abs(Q))
		yq[n.isnan(yq)] = 0.0
		
		# for the first branch
		
		sumyq = n.sum(yq[branch[i,:]])
		sumy = n.sum(y[branch[i,:]])
		
		dQ = -1*sumy/sumyq
		print("dQ = %f" % 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)

Warning: invalid value encountered in divide dQ = 70.284711 dQ = -18.083459 dQ = -6.477074 dQ = 0.765176 dQ = -0.338565 dQ = 0.067968 Discharges [m^3/hr]: [ 500. 313.46908569 232.74967957 -186.53092957 46.21875763 -267.25033569 500. ] Equivalent length = 7.44 km