3924.sagews

"""
 Newton Rapson method solution of the example of 3 reservouirs
|_____| H1
       \  |_______| H2
        \/
        /HD
       /
|______|H3
"""
import numpy as np
import matplotlib.pyplot as p

no = np.array([1,2,3]) # pipe no
L = np.array([3,1,5]) # L [km]
D = np.array([20,20,25]) # [cm]
C = np.array([120,110,125]) # Hazen Williams
H = np.array([100,80,20]) # heads, m 
n = 1.852

# print L, D, C, H, n

R = L*1.526e7/((C**n)*(D**4.87)) # resistances


dQ = np.array([np.Inf, np.Inf]) 

Q = np.array([100, 100])# Q(3) = Q(1) + Q(2);

while np.max(np.abs(dQ)) > 1e-2: # repeat until the corrections are small
# 	f(Q1,Q2) = R1*Q1^n - R2*Q2^n - (H1-H2) = 0
	f = R[0]*Q[0]**n - R[1]*Q[1]**n - (H[0] - H[1])

# g(Q1,Q2) = R1*Q1^n + R3*(Q1 + Q2)^n - (H1-H3) = 0
	g = R[0]*Q[0]**n + R[2]*(Q[0]+Q[1])**n - (H[0] - H[2])

# corrections: dQ1, dQ2 are solved using the equations:
# n*R(1)*Q(1)^(n-1)*dQ1 - n*R2*Q(2)^(n-1)*dQ2 = -f;
# n*(R(1)*Q(1)^(n-1) + R(3)*(Q(1)+Q(2))^(n-1))*dQ1 + n*R(3)*(Q(1)+Q(2))^(n-1)*dQ2 = -g;
	LHS = np.array([[n*R[0]*Q[0]**(n-1), -n*R[1]*Q[1]**(n-1)]  , 
	[n*(R[0]*Q[0]**(n-1) + R[2]*(Q[0]+Q[1])**(n-1)),  n*R[2]*(Q[0]+Q[1])**(n-1)]])
        
	RHS = np.array([-f,-g])
	
	# print LHS, RHS

	dQ = np.linalg.solve(LHS,RHS)
	print dQ

	Q = Q + dQ
	# print Q


Q.resize((3))
Q[2] = Q[0] + Q[1] # flowing out of the node

print "discharges:"
print Q

HD = H[0] - R[0]*Q[0]**n # head in the node 
print "head at D %3.2f m" % HD
# or
HD = H[1] - R[1]*Q[1]**n
print "or %3.2f m" % HD
#  or
HD = H[2] + R[2]*Q[2]**n
print "or %3.2f m" % HD
print HD

[ 51.65948884 32.77065612] [-7.73483889 -3.36967404] [-0.17167735 -0.02473909] [ -8.13805980e-05 1.46059453e-05] discharges: [ 143.75289121 129.37625759 273.12914881] head at D 70.49 m or 70.49 m or 70.49 m 70.4932101892