from numpy import matrix, multiply, power, abs, transpose, sign, sum
from numpy.linalg import inv
L_pipe=matrix([3,2,3,2,3,2,3]) #[km]
D_pipe=matrix([30,20,20,20,20,20,20]) #[cm]
CHW_pipe=matrix([100,100,100,100,100,100,100]) # Friction coefficient of Hazen Williams
_=multiply(power(D_pipe,4.87),power(CHW_pipe,1.852))
r_pipe =1.526e7*L_pipe/_
# for each pipe the initial guess
Q0_pipe =matrix([110,30,110,110,80,60,140]) #[m^3/hr]
# for each node, the inputs/outputs
Qd = matrix([-220,0,0,0,20,200]).transpose() #[m^3/hr]
H1_assumed = 100.0 #[m]
Q_pipe=Q0_pipe
H=0.0
OLD_H=1.0
N_iter = 100
# start the loop
while sum(abs(OLD_H-H))> 0.01:
# for each pipe, get the 'm'
M = -1.0/multiply(r_pipe,power(abs(Q_pipe),0.852))
# construct the LHS
LHS=matrix([[-(M[0,0]+M[0,3]),M[0,0],0,M[0,3],0,0],
[M[0,0],-(M[0,0]+M[0,1]+M[0,4]),M[0,1],0,M[0,4],0],
[0,M[0,1],-(M[0,1]+M[0,2]+M[0,6]),M[0,2],0,M[0,6]],
[M[0,3],0,M[0,2],-(M[0,2]+M[0,3]),0,0],
[0,M[0,4],0,0,-(M[0,4]+M[0,5]),M[0,5]],
[0,0,M[0,6],0,M[0,5],-(M[0,5]+M[0,6])]])
LHS[0,0]=10**10
# right hand side vector = outputs
RHS = -Qd
RHS[0,0]=H1_assumed*10**5
OLD_H=H
# the solution is to invert the matrix and solve the equation
H=inv(LHS)*RHS
hloss = matrix([H[0,0]-H[1,0], H[1,0]-H[2,0], H[3,0]-H[2,0] , H[0,0]-H[3,0] , H[1,0]-H[4,0], H[4,0]-H[5,0] , H[2,0]-H[5,0]])
_ = power(abs(hloss)/r_pipe,1.0/1.852)
Q_pipe=multiply(sign(abs(hloss)),_)
print "Heads:\n", H
print ""
print "Discharges: \n", Q_pipe.transpose()
|
|
Heads:
[[ 9.99978003e-04]
[ -7.23216343e+00]
[ -1.24834360e+01]
[ -4.99277441e+00]
[ -2.99969110e+01]
[ -4.02213640e+01]]
Discharges:
[[ 162.85041749]
[ 58.71301737]
[ 57.14050793]
[ 57.14050793]
[ 104.1368128 ]
[ 84.13654684]
[ 115.86135711]]
Heads:
[[ 9.99978003e-04]
[ -7.23216343e+00]
[ -1.24834360e+01]
[ -4.99277441e+00]
[ -2.99969110e+01]
[ -4.02213640e+01]]
Discharges:
[[ 162.85041749]
[ 58.71301737]
[ 57.14050793]
[ 57.14050793]
[ 104.1368128 ]
[ 84.13654684]
[ 115.86135711]]
|