"""
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]
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]
|