"""
The tsnet.simulation.solver module contains methods to solver MOC
for different grid configurations, including:
1. inner_node
2. valve_node
3. pump_node
4. source_pump
5. valve_end
6. dead_end
7. rev_end
8. add_leakage
"""
from __future__ import print_function
import numpy as np
import warnings
[docs]def Reynold(V, D):
""" Calculate Reynold number
Parameters
----------
V : float
velocity
D : float
diameter
Returns
-------
Re : float
Reynold number
"""
nu = 1.004e-6 # kinematic viscosity [m^2/s]
Re = np.abs(V*D/nu)
return Re
[docs]def quasi_steady_friction_factor(Re, KD):
""" Update friction factor based on Reynold number
Parameters
----------
Re : float
velocity
KD : float
relative roughness height (K/D)
Returns
-------
f : float
quasi-steady friction factor
"""
a = -1.8*np.log10(6.9/Re + KD)
f = (1./a)**2.
return f
[docs]def unsteady_friction(Re, dVdt, dVdx, V, a, g):
""" Calculate unsteady friction
Parameters
----------
Re : float
velocity
dVdt : float
local instantaneous acceleration
dVdx : float
instantaneous convective acceleration
V : float
velocity
a : float
wave speed
g: float
gravitational acceleration
Returns
-------
Ju : float
unsteady friction factor
"""
# calculate Vardy's shear decay coefficient (C)
if Re< 2000: # laminar flow
C = 4.76e-3
else:
C = 7.41 / Re**(np.log10(14.3/Re**0.05))
# calculate Brunone's friction coefficient
k = np.sqrt(C)/2.
" TO DO: check the sign of unsteady friction"
Ju = k/g/2.* (dVdt + a* np.sign(V) * np.abs(dVdx))
return Ju
[docs]def cal_friction(friction, f, D, V, KD, dt, dVdt, dVdx, a, g ):
""" Calculate friction term
Parameters
----------
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
f : float
steady friction factor
D : float
pipe diameter
V : float
pipe flow velocity
KD : float
relative roughness height
dt : float
time step
dVdt : float
local instantaneous acceleration
dVdx : float
convective instantaneous acceleration
a : float
wave speed
g : float
gravitational accelerations
Returns
-------
float
total friction, including steady and unsteady
"""
tol = 1
if friction == 'steady':
Ju = 0
Js = f*dt/2./D*V*abs(V) #steady friction
else:
Re = Reynold(V, D)
if Re < tol:
Js = 0
else:
f = quasi_steady_friction_factor(Re, KD)
Js = f*dt/2./D*V*abs(V)
if friction == 'quasi-steady':
Ju = 0
elif friction == 'unsteady':
Ju = unsteady_friction(Re, dVdt, dVdx, V, a, g)
return Ju + Js
[docs]def cal_Cs( link1, link2, H1, V1, H2, V2, s1, s2, g, dt,
friction, dVdx1, dVdx2, dVdt1, dVdt2):
"""Calculate coefficients for MOC characteristic lines
Parameters
----------
link1 : object
Pipe object of C+ charateristics curve
link2 : object
Pipe object of C- charateristics curve
H1 : list
List of the head of C+ charateristics curve
V1 : list
List of the velocity of C+ charateristics curve
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
s1 : list
List of signs that represent the direction of the flow
in C+ charateristics curve
s2 : list
List of signs that represent the direction of the flow
in C- charateristics curve
dt : float
Time step
g : float
Gravity acceleration
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx1 : list
List of convective instantaneous acceleration on the
C+ characteristic curve
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt1 : list
List of local instantaneous acceleration on the
C+ characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
Returns
-------
A1: list
list of left adjacent pipe cross-section area
A2: list
list of right adjacent pipe cross-section area
C1: list
list of left adjacent pipe MOC coefficients
C2: list
list of right adjacent pipe MOC coefficients
"""
# property of left adjacent pipe
f1 = [link1[i].roughness for i in range(len(link1))] # unitless
D1 = [link1[i].diameter for i in range(len(link1))] # m
a1 = [link1[i].wavev for i in range(len(link1))] # m/s
A1 = [np.pi * D1[i]**2. / 4. for i in range(len(link1))] # m^2
C1 = np.zeros((len(link1),2), dtype=np.float64)
theta1 = [link1[i].theta for i in range((len(link1)))]
KD1 = [link1[i].roughness_height for i in range(len(link1))]
for i in range(len(link1)):
# J = f1[i]*dt/2./D1[i]*V1[i]*abs(V1[i])
J = cal_friction(friction, f1[i], D1[i], V1[i], KD1[i],
dt, dVdt1[i], dVdx1[i], a1[i], g )
C1[i,0] = s1[i]*V1[i] + g/a1[i]*H1[i] - s1[i]*J + g/a1[i]* dt *V1[i]*theta1[i]
C1[i,1] = g/a1[i]
# property of right adjacent pipe
f2 = [link2[i].roughness for i in range(len(link2))] # unitless
D2 = [link2[i].diameter for i in range(len(link2))] # m
a2 = [link2[i].wavev for i in range(len(link2))] # m/s
A2 = [np.pi * D2[i]**2. / 4. for i in range(len(link2))] # m^2
C2 = np.zeros((len(link2),2),dtype=np.float64)
theta2 = [link2[i].theta for i in range((len(link2)))]
KD2 = [link2[i].roughness_height for i in range(len(link2))]
for i in range(len(link2)):
# J = f2[i]*dt/2./D2[i]*V2[i]*abs(V2[i])
J = cal_friction(friction, f2[i], D2[i], V2[i], KD2[i],
dt, dVdt2[i], dVdx2[i], a2[i], g)
C2[i,0] = s2[i]*V2[i] + g/a2[i]*H2[i] - s2[i]* J + g/a2[i]* dt *V2[i]*theta2[i]
C2[i,1] = g/a2[i]
return A1, A2, C1, C2
[docs]def inner_node_unsteady(link, H0, V0, dt, g, dVdx, dVdt):
"""Inner boundary MOC using C+ and C- characteristic curve with unsteady friction
Parameters
----------
link : object
current pipe
H0 : list
head at previous time step
V0 : list
velocity at previous time step
dt : float
Time step
g : float
Gravity acceleration
dVdx : list
List of convective instantaneous acceleration
dVdt : list
List of local instantaneous acceleration
Returns
-------
HP : float
Head at current pipe inner nodes at current time
VP : float
Velocity at current pipe inner nodes at current time
"""
HP = np.zeros(len(H0))
VP = np.zeros(len(V0))
# property of current pipe
f = link.roughness # unitless
D = link.diameter # m
a = link.wavev # m/s
# A = np.pi * D**2. / 4. # m^2
theta = link.theta
KD = link.roughness_height
ga = g/a
tol = 1e-1
for i in range(1,len(H0)-1):
V1 = V0[i-1]; H1 = H0[i-1]
V2 = V0[i+1]; H2 = H0[i+1]
dVdx1 = dVdx[i-1] ; dVdx2 = dVdx[i]
dVdt1 = dVdt[i-1] ; dVdt2 = dVdt[i+1]
C = np.zeros((2,1), dtype=np.float64)
Re = Reynold(V1, D)
if Re <tol:
Js = 0
else:
f = quasi_steady_friction_factor(Re, KD)
Js = f*dt/2./D*V1*abs(V1)
Ju = unsteady_friction(Re, dVdt1, dVdx1, V1, a, g)
J1 = Js +Ju
C[0,0] = V1 + ga*H1 - J1 + ga* dt *V1*theta
Re = Reynold(V2, D)
if Re < tol:
Js = 0
else:
f = quasi_steady_friction_factor(Re, KD)
Js = f*dt/2./D*V2*abs(V2)
Ju = unsteady_friction(Re, dVdt2, dVdx2, V2, a, g)
J2 = Js +Ju
C[1,0] = -V2+ ga*H2 + J2 + ga* dt *V2*theta
HP[i] = (C[0,0] + C[1,0]) / 2./ga
VP[i] = np.float64(-C[1,0]+ ga*HP[i])
return HP[1:-1], VP[1:-1]
[docs]def inner_node_quasisteady(link, H0, V0, dt, g):
"""Inner boundary MOC using C+ and C- characteristic curve with unsteady friction
Parameters
----------
link : object
current pipe
H0 : list
head at previous time step
V0 : list
velocity at previous time step
dt : float
Time step
g : float
Gravity acceleration
dVdx : list
List of convective instantaneous acceleration
dVdt : list
List of local instantaneous acceleration
Returns
-------
HP : float
Head at current pipe inner nodes at current time
VP : float
Velocity at current pipe inner nodes at current time
"""
HP = np.zeros(len(H0))
VP = np.zeros(len(V0))
# property of current pipe
D = link.diameter # m
a = link.wavev # m/s
theta = link.theta
KD = link.roughness_height
ga = g/a
tol = 1e-1
for i in range(1,len(H0)-1):
V1 = V0[i-1]; H1 = H0[i-1]
V2 = V0[i+1]; H2 = H0[i+1]
C = np.zeros((2,1), dtype=np.float64)
Re = Reynold(V1, D)
if Re < tol:
J1 = 0
else:
f = quasi_steady_friction_factor(Re, KD)
J1 = f*dt/2./D*V1*abs(V1)
Re = Reynold(V2, D)
if Re < tol:
J2 = 0
else:
f = quasi_steady_friction_factor(Re, KD)
J2 = f*dt/2./D*V2*abs(V2)
C[0,0] = V1 + ga*H1 - J1 + ga* dt *V1*theta
C[1,0] = -V2+ ga*H2 + J2 + ga* dt *V2*theta
HP[i] = (C[0,0] + C[1,0]) / 2./ ga
VP[i] = np.float64(-C[1,0]+ ga*HP[i])
return HP[1:-1], VP[1:-1]
[docs]def inner_node_steady(link, H0, V0, dt, g):
"""Inner boundary MOC using C+ and C- characteristic curve with unsteady friction
Parameters
----------
link : object
current pipe
H0 : list
head at previous time step
V0 : list
velocity at previous time step
dt : float
Time step
g : float
Gravity acceleration
Returns
-------
HP : float
Head at current pipe inner nodes at current time
VP : float
Velocity at current pipe inner nodes at current time
"""
HP = np.zeros(len(H0))
VP = np.zeros(len(V0))
# property of current pipe
f = link.roughness # unitless
D = link.diameter # m
a = link.wavev # m/s
# A = np.pi * D**2. / 4. # m^2
theta = link.theta
ga = g/a
for i in range(1,len(H0)-1):
V1 = V0[i-1]; H1 = H0[i-1]
V2 = V0[i+1]; H2 = H0[i+1]
C = np.zeros((2,2), dtype=np.float64)
J1 = f*dt/2./D*V1*abs(V1)
C[0,0] = V1 + ga*H1 - J1 + ga* dt *V1*theta
C[0,1] = ga
J2 = f*dt/2./D*V2*abs(V2)
C[1,0] = -V2+ ga*H2 + J2 + ga* dt *V2*theta
C[1,1] = ga
HP[i] = (C[0,0] + C[1,0]) / (C[0,1] + C[1,1])
VP[i] = np.float64(-C[1,0]+ C[1,1]*HP[i])
return HP[1:-1], VP[1:-1]
[docs]def valve_node(KL_inv, link1, link2, H1, V1, H2, V2, dt, g, nn, s1, s2,
friction, dVdx1, dVdx2, dVdt1, dVdt2):
"""Inline valve node MOC calculation
Parameters
----------
KL_inv : int
Inverse of the valve loss coefficient at current time
link1 : object
Pipe object of C+ charateristics curve
link2 : object
Pipe object of C- charateristics curve
H1 : list
List of the head of C+ charateristics curve
V1 : list
List of the velocity of C+ charateristics curve
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
dt : float
Time step
g : float
Gravity acceleration
nn : int
The index of the calculation node
s1 : list
List of signs that represent the direction of the flow
in C+ charateristics curve
s2 : list
List of signs that represent the direction of the flow
in C- charateristics curve
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx1 : list
List of convective instantaneous acceleration on the
C+ characteristic curve
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt1 : list
List of local instantaneous acceleration on the
C+ characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
try :
list(link1)
except:
link1 = [link1]
V1 = [V1] ; H1 = [H1]
dVdx1 = [dVdx1]; dVdt1 = [dVdt1]
try :
list(link2)
except:
link2 = [link2]
V2 = [V2] ; H2 = [H2]
dVdx2 = [dVdx2]; dVdt2 = [dVdt2]
A1, A2, C1, C2 = cal_Cs(link1, link2, H1, V1, H2, V2, s1, s2, g, dt,
friction, dVdx1, dVdx2, dVdt1, dVdt2)
# parameters of the quadratic polynomial
aq = 1
bq = 2*g*KL_inv* (A2[0]/A1[0]/C1[0,1] + 1/C2[0,1])
cq = 2*g*KL_inv* (C2[0,0]/C2[0,1] - C1[0,0]/C1[0,1])
# solve the quadratic equation
delta = bq**2 - 4*aq*cq
if delta >= 0:
VP = (-bq + np.sqrt(delta))/(2*aq)
elif delta > -1.0e-7 and delta <0 :
VP = (-bq)/(2*aq)
else:
VP = (-bq)/(2*aq)
warnings.warn('Error: The quadratic equation has no real solution (valve)')
if VP >=0 : # positive flow
if nn == 0: # pipe start
VP = VP
HP = (C2[0,0] + VP) / C2[0,1]
else: # pipe end
VP = VP*A2[0]/A1[0]
HP = (C1[0,0] - VP) / C1[0,1]
else : # reverse flow
# reconstruct the quadratic equation
# parameters of the quadratic polynomial
aq = 1
bq = 2*g*KL_inv* (-A1[0]/A2[0]/C2[0,1]-1/C1[0,1])
cq = 2*g*KL_inv* (-C2[0,0]/C2[0,1]+C1[0,0]/C1[0,1])
# solve the quadratic equation
delta = bq**2 - 4*aq*cq
if delta >= 0:
VP = (-bq - np.sqrt(delta))/(2*aq)
elif delta > -1.0e-7 and delta <0 :
VP = (-bq)/(2*aq)
else:
VP = (-bq)/(2*aq)
warnings.warn('Error: The quadratic equation has no real solution (valve)')
if nn == 0: # pipe start
VP = VP*A1[0]/A2[0]
HP = (C2[0,0] + VP ) / C2[0,1]
else: # pipe end
VP = VP
HP = (C1[0,0] - VP) / C1[0,1]
return HP, VP
[docs]def pump_node(pumpc,link1, link2, H1, V1, H2, V2, dt, g, nn, s1, s2,
friction, dVdx1, dVdx2, dVdt1, dVdt2):
""" Inline pump node MOC calculation
Parameters
----------
pumpc : list
Parameters (a, b,c) to define pump characteristic cure,
so that
.. math:: h_p = a*Q**2 + b*Q + c
link1 : object
Pipe object of C+ charateristics curve
link2 : object
Pipe object of C- charateristics curve
H1 : list
List of the head of C+ charateristics curve
V1 : list
List of the velocity of C+ charateristics curve
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
dt : float
Time step
g : float
Gravity acceleration
nn : int
The index of the calculation node
s1 : list
List of signs that represent the direction of the flow
in C+ charateristics curve
s2 : list
List of signs that represent the direction of the flow
in C- charateristics curve
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx1 : list
List of convective instantaneous acceleration on the
C+ characteristic curve
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt1 : list
List of local instantaneous acceleration on the
C+ characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
try :
list(link1)
except:
link1 = [link1]
V1 = [V1] ; H1 = [H1]
dVdx1 = [dVdx1]; dVdt1 = [dVdt1]
try :
list(link2)
except:
link2 = [link2]
V2 = [V2] ; H2 = [H2]
dVdx2 = [dVdx2]; dVdt2 = [dVdt2]
A1, A2, C1, C2 = cal_Cs( link1, link2, H1, V1, H2, V2, s1, s2, g, dt,
friction, dVdx1, dVdx2, dVdt1, dVdt2)
# pump power function
ap, bp, cp = pumpc[0]
ap = ap * A1[0]**2.
bp = bp * A1[0]
# parameters of the quadratic polynomial
aq = 1
bq = 1/ap * (bp - 1/C1[0,1] - A1[0]/C2[0,1]/A2[0])
cq = 1/ap * (-C2[0,0]/C2[0,1] + C1[0,0]/C1[0,1] + cp)
# solve the quadratic equation
delta = bq**2. - 4.*aq*cq
if delta >= 0:
VP = (-bq + np.sqrt(delta))/(2*aq)
elif delta > -1.0e-7 and delta <0 :
VP = (-bq)/(2*aq)
else:
VP = (-bq)/(2*aq)
warnings.warn('Error: The quadratic equation has no real solution (pump)')
hp = ap*VP**2. + bp*VP + cp # head gain
if VP > 0 and hp >=0 : # positive flow & positive head gain
if nn == 0: # pipe start
VP = VP*A1[0]/A2[0]
HP = (C2[0,0] + VP ) / C2[0,1]
else: # pipe end
VP = VP
HP = (C1[0,0] - VP) / C1[0,1]
elif VP<0 :
warnings.warn( "Reverse flow stopped by check valve!")
VP = 0
if nn == 0: # pipe start
HP = (C2[0,0] + VP ) / C2[0,1]
else :
HP = (C1[0,0] - VP) / C1[0,1]
# hp = cp
# # suction or discharge side?
# if pumpc[1] == "s": # suction side
# if nn == 0: # pipe start
# HP = (C2[0,0] + VP ) / C2[0,1]
# else :
# HP = (C1[0,0] - VP) / C1[0,1]
# else: #discharge
# if nn == 0: # pipe start
# HP = (C1[0,0] - VP) / C1[0,1] + hp
# else :
# HP = (C2[0,0] + VP ) / C2[0,1] + hp
else: # positive flow and negative head gain
warnings.warn( "Negative head gain activates by-pass!")
hp = 0
# suction or discharge side?
if pumpc[1] == "s": # suction side
if nn == 0: # pipe start
HP = (C2[0,0] + VP ) / C2[0,1]
else :
HP = (C1[0,0] - VP) / C1[0,1]
else:
if nn == 0: # pipe start
HP = (C1[0,0] - VP) / C1[0,1] + hp
else :
HP = (C2[0,0] + VP ) / C2[0,1] +hp
return HP, VP
[docs]def source_pump(pump, link2, H2, V2, dt, g, s2,
friction, dVdx2, dVdt2):
"""Source Pump boundary MOC calculation
Parameters
----------
pump : list
pump[0]: elevation of the reservoir/tank
pump[1]: Parameters (a, b,c) to define pump characteristic cure,
so that
.. math:: h_p = a*Q**2 + b*Q + c
link2 : object
Pipe object of C- charateristics curve
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
dt : float
Time step
g : float
Gravity acceleration
s2 : list
List of signs that represent the direction of the flow
in C- charateristics curve
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
pumpc = pump[1]
Hsump = pump[0]
try :
list(link2)
except:
link2 = [link2]
V2 = [V2] ; H2 = [H2]
dVdx2 = [dVdx2]; dVdt2 = [dVdt2]
_, A2, _, C2 = cal_Cs( link2, link2, H2, V2, H2, V2, s2, s2, g, dt,
friction, dVdx2, dVdx2, dVdt2, dVdt2)
# pump power function
ap, bp, cp = pumpc
ap = ap * A2[0]**2.
bp = bp * A2[0]
# parameters of the quadratic polynomial
aq = ap * C2[0,1]**2.
bq = bp*C2[0,1] - 2.*ap*C2[0,0]*C2[0,1] - 1
cq = ap*C2[0,0]**2. - bp*C2[0,0] + Hsump + cp
# solve the quadratic equation
delta = bq**2. - 4.*aq*cq
if delta >= 0:
HP = (-bq - np.sqrt(delta))/(2*aq)
elif delta > -1.0e-7 and delta <0 :
HP = (-bq)/(2*aq)
else:
HP = (-bq)/(2*aq)
warnings.warn('The quadratic equation has no real solution (pump)')
if HP > Hsump:
VP = np.float64(-C2[0,0] + C2[0,1]*HP)
else :
HP = Hsump
VP = np.float64(-C2[0,0] + C2[0,1]*HP)
if VP <= 0 : # positive flow
warnings.warn( "Reverse flow stopped by check valve!")
VP = 0
HP = (C2[0,0] + VP ) / C2[0,1]
return HP, VP
[docs]def valve_end(H1, V1, V, nn, a, g, f, D, dt,
KD, friction, dVdx2, dVdt2):
""" End Valve boundary MOC calculation
Parameters
----------
H1 : float
Head of the C+ charateristics curve
V1 : float
Velocity of the C+ charateristics curve
V : float
Velocity at the valve end at current time
nn : int
The index of the calculation node
a : float
Wave speed at the valve end
g : float
Gravity acceleration
f : float
friction factor of the current pipe
D : float
diameter of the current pipe
dt : float
Time step
KD : float
relative roughness height
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
J = cal_friction(friction, f, D, V, KD, dt, dVdt2, dVdx2, a, g )
if nn == 0 :
# HP = H1 + a/g*(V - V1) + a/g*f*dt/(2.*D)*V1*abs(V1)
HP = H1 + a/g*(V - V1) + a/g*J
VP = V
else :
HP = H1 - a/g*(V - V1) - a/g*J
VP = V
return HP,VP
[docs]def dead_end(linkp, H1, V1, elev, nn, a, g, f, D, dt,
KD, friction, dVdx1, dVdt1):
"""Dead end boundary MOC calculation with pressure dependant demand
Parameters
----------
linkp : object
Current pipe
H1 : float
Head of the C+ charateristics curve
V1 : float
Velocity of the C+ charateristics curve
elev : float
Elevation at the dead end node
nn : int
The index of the calculation node
a : float
Wave speed at the valve end
g : float
Gravity acceleration
f : float
friction factor of the current pipe
D : float
diameter of the current pipe
dt : float
Time step
KD : float
relative roughness height
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx1 : list
List of convective instantaneous acceleration on the
C+ characteristic curve
dVdt1 : list
List of local instantaneous acceleration on the
C+ characteristic curve
"""
A = np.pi/4. * linkp.diameter**2.
J = cal_friction(friction, f, D, V1, KD, dt, dVdt1, dVdx1, a, g )
if nn == 0: # dead end is the start node of a pipe
k = linkp.start_node.demand_coeff + linkp.start_node.emitter_coeff
aq = 1
bq = -a/g*k/A
# cq = a/g *V1 - a/g*f*dt/(2.*D)*V1*abs(V1) - H1 - g/a*dt*V1*linkp.theta + elev
cq = a/g *V1 - a/g*J - H1 - g/a*dt*V1*linkp.theta + elev
# solve the quadratic equation
delta = bq**2. - 4.*aq*cq
if delta >= 0:
HP = (-bq - np.sqrt(delta))/(2*aq)
HP = HP**2. + elev
elif delta > -1.0e-7 and delta <0 :
HP = (-bq)/(2*aq)
HP = HP**2. +elev
else:
HP = (-bq)/(2*aq)
HP = HP**2. +elev
warnings.warn("""The quadratic equation has no real solution (dead end).
The results might not be accurate.""")
VP = V1 - g/a*H1 - f*dt/(2.*D)*V1*abs(V1) + g/a*HP - g/a*dt*V1*linkp.theta
else : # dead end is the end node of a pipe
k = linkp.end_node.demand_coeff + linkp.end_node.emitter_coeff
aq = 1
bq = a/g*k/A
# cq = -a/g *V1 + a/g*f*dt/(2.*D)*V1*abs(V1) - H1 - g/a*dt*V1*linkp.theta + elev
cq = -a/g *V1 + a/g*J - H1 - g/a*dt*V1*linkp.theta + elev
# solve the quadratic equation
delta = bq**2. - 4.*aq*cq
if delta >= 0:
HP = (-bq + np.sqrt(delta))/(2*aq)
HP = HP**2. + elev
elif delta > -1.0e-7 and delta <0 :
HP = (-bq)/(2*aq)
HP = HP**2. + elev
else:
HP = (-bq)/(2*aq)
HP = HP**2. + elev
warnings.warn("The quadratic equation has no real solution (dead end).\
The results might not be accurate.")
VP = V1 + g/a *H1 - f*dt/(2.*D)*V1*abs(V1) - g/a*HP + g/a*dt*V1*linkp.theta
return HP,VP
[docs]def rev_end( H2, V2, H, nn, a, g, f, D, dt,
KD, friction, dVdx2, dVdt2):
"""Reservoir/ Tank boundary MOC calculation
Parameters
----------
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
H : float
Head of the reservoir/tank
nn : int
The index of the calculation node
a : float
Wave speed at the valve end
g : float
Gravity acceleration
f : float
friction factor of the current pipe
D : float
diameter of the current pipe
dt : float
Time step
KD : float
relative roughness height
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
J = cal_friction(friction, f, D, V2, KD, dt, dVdt2, dVdx2, a, g )
if nn == 0 :
VP = V2 + g/a*(H - H2) - J
HP = H
else:
VP = V2 - g/a*(H - H2) - J
HP = H
return HP, VP
[docs]def add_leakage(emitter_coef, block_per, link1, link2, elev,
H1, V1, H2, V2, dt, g, nn, s1, s2,
friction, dVdx1=0, dVdx2=0, dVdt1=0, dVdt2=0):
r"""Leakage Node MOC calculation
Parameters
----------
emitter_coef : float
float, optional
Required if leak_loc is defined
The leakage coefficient of the leakage
.. math:: Q_leak = leak_A [ m^3/s/(m H20)^(1/2)] * \sqrt(H)
link1 : object
Pipe object of C+ charateristics curve
link2 : object
Pipe object of C- charateristics curve
H1 : list
List of the head of C+ charateristics curve
V1 : list
List of the velocity of C+ charateristics curve
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
dt : float
Time step
g : float
Gravity acceleration
nn : int
The index of the calculation node
s1 : list
List of signs that represent the direction of the flow
in C+ charateristics curve
s2 : list
List of signs that represent the direction of the flow
in C- charateristics curve
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx1 : list
List of convective instantaneous acceleration on the
C+ characteristic curve
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt1 : list
List of local instantaneous acceleration on the
C+ characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
emitter_coef = emitter_coef # m^3/s//(m H2O)^(1/2)
try :
list(link1)
except:
link1 = [link1]
V1 = [V1] ; H1 = [H1]
dVdx1 = [dVdx1]; dVdt1 = [dVdt1]
try :
list(link2)
except:
link2 = [link2]
V2 = [V2] ; H2 = [H2]
dVdx2 = [dVdx2]; dVdt2 = [dVdt2]
A1, A2, C1, C2 = cal_Cs(link1, link2, H1, V1, H2, V2, s1, s2, g, dt,
friction, dVdx1, dVdx2, dVdt1, dVdt2)
a = np.dot(C1[:,0], A1) + np.dot(C2[:,0],A2)
b = np.dot(C1[:,1], A1) + (1-block_per)*np.dot(C2[:,1],A2)
# parameters of the quadratic polynomial
# a1 = b**2.
# b1 = -(2.*a*b +emitter_coef**2)
# c1 = a**2.
a1 = b**2
b1 = -2*a*b - emitter_coef**2.
c1 = a**2 + emitter_coef**2.*elev
# solve the quadratic equation
delta = b1**2 - 4*a1*c1
if delta >= 0:
HP = (-b1 - np.sqrt(delta))/(2*a1)
elif delta > -1.0e-7 and delta <0 :
HP = (-b1)/(2*a1)
else:
HP = (-b1)/(2*a1)
warnings.warn('Error: The quadratic equation has no real solution (leakage)')
if nn == 0: # pipe start
VP = np.float64(-C2[:,0]+ C2[:,1]*HP)
else: # pipe end
VP = np.float64(C1[:,0] - C1[:,1]*HP)
return HP, VP
[docs]def surge_tank(tank, link1, link2, H1, V1, H2, V2, dt, g, nn, s1, s2,
friction, dVdx1, dVdx2, dVdt1, dVdt2):
"""Surge tank node MOC calculation
Parameters
----------
tank : int
tank shape parameters
[As, z, Qs]
As : cross-sectional area of the surge tank
z : water level in the surge tank at previous time step
Qs : water flow into the tank at last time step
link1 : object
Pipe object of C+ charateristics curve
link2 : object
Pipe object of C- charateristics curve
H1 : list
List of the head of C+ charateristics curve
V1 : list
List of the velocity of C+ charateristics curve
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
dt : float
Time step
g : float
Gravity acceleration
nn : int
The index of the calculation node
s1 : list
List of signs that represent the direction of the flow
in C+ charateristics curve
s2 : list
List of signs that represent the direction of the flow
in C- charateristics curve
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx1 : list
List of convective instantaneous acceleration on the
C+ characteristic curve
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt1 : list
List of local instantaneous acceleration on the
C+ characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
try :
list(link1)
except:
link1 = [link1]
V1 = [V1] ; H1 = [H1]
dVdx1 = [dVdx1]; dVdt1 = [dVdt1]
try :
list(link2)
except:
link2 = [link2]
V2 = [V2] ; H2 = [H2]
dVdx2 = [dVdx2]; dVdt2 = [dVdt2]
A1, A2, C1, C2 = cal_Cs(link1, link2, H1, V1, H2, V2, s1, s2, g, dt,
friction, dVdx1, dVdx2, dVdt1, dVdt2)
As, z, Qs = tank
at = 2.* As/dt
HP = ((np.dot(C1[:,0], A1) + np.dot(C2[:,0],A2) + at*z + Qs ) /
(np.dot(C1[:,1], A1) + np.dot(C2[:,1],A2) + at))
VP2 = -C2[:,0]+ C2[:,1]*HP
VP1 = C1[:,0] - C1[:,1]*HP
QPs = (np.sum(np.array(VP1)*np.array(A1)) -
np.sum(np.array(VP2)*np.array(A2)))
if nn == 0: # pipe start
VP = np.float64(VP2)
else: # pipe end
VP = np.float64(VP1)
return HP, VP, QPs
[docs]def air_chamber(tank, link1, link2, H1, V1, H2, V2, dt, g, nn, s1, s2,
friction, dVdx1, dVdx2, dVdt1, dVdt2):
"""Surge tank node MOC calculation
Parameters
----------
tank : int
tank shape parameters
[As, ht, C, z, Qs]
As : cross-sectional area of the surge tank
ht : tank height
C : air constant
z : water level in the surge tank at previous time step
Qs : water flow into the tank at last time step
link1 : object
Pipe object of C+ charateristics curve
link2 : object
Pipe object of C- charateristics curve
H1 : list
List of the head of C+ charateristics curve
V1 : list
List of the velocity of C+ charateristics curve
H2 : list
List of the head of C- charateristics curve
V2 : list
List of the velocity of C- charateristics curve
dt : float
Time step
g : float
Gravity acceleration
nn : int
The index of the calculation node
s1 : list
List of signs that represent the direction of the flow
in C+ charateristics curve
s2 : list
List of signs that represent the direction of the flow
in C- charateristics curve
friction : str
friction model, e.g., 'steady', 'quasi-steady', 'unsteady',
by default 'steady'
dVdx1 : list
List of convective instantaneous acceleration on the
C+ characteristic curve
dVdx2 : list
List of convective instantaneous acceleration on the
C- characteristic curve
dVdt1 : list
List of local instantaneous acceleration on the
C+ characteristic curve
dVdt2 : list
List of local instantaneous acceleration on the
C- characteristic curve
"""
try :
list(link1)
except:
link1 = [link1]
V1 = [V1] ; H1 = [H1]
dVdx1 = [dVdx1]; dVdt1 = [dVdt1]
try :
list(link2)
except:
link2 = [link2]
V2 = [V2] ; H2 = [H2]
dVdx2 = [dVdx2]; dVdt2 = [dVdt2]
A1, A2, C1, C2 = cal_Cs(link1, link2, H1, V1, H2, V2, s1, s2, g, dt,
friction, dVdx1, dVdx2, dVdt1, dVdt2)
# parameters
Hb = 10.3 # barometric pressure head
m = 1.2
As, ht, C, z, Qs = tank # tank properties and results at last time step
at = 2.* As/dt
Va = (ht-z)*As # air volume at last time step
Cor = 0
a = np.dot(C1[:,0], A1) + np.dot(C2[:,0],A2)
b = np.dot(C1[:,1], A1) + np.dot(C2[:,1],A2)
def tank_flow(QPs, Qs, a, b, As, ht, C, z, at, Va, Cor, m, Hb):
return (((a-QPs)/b + Hb - z - (Qs+QPs)/at - Cor*QPs*np.abs(QPs))
* (Va- (Qs+QPs)*As/at)**m - C)
def tank_flow_prime(QPs, Qs, a, b, As, ht, C, z, at, Va, Cor, m, Hb):
p1 = (-m*As/at * (Va- (Qs+QPs)*As/at)**(m-1)*
((a-QPs)/b + Hb - z - (Qs+QPs)/at - Cor*QPs*np.abs(QPs)))
p2 = (-1/b -1/at - Cor*2.*QPs*np.sign(QPs))* (Va- (Qs+QPs)*As/at)**m
return p1+p2
# solve nonlinear equation for tank flow at this time step
from scipy import optimize
QPs = optimize.newton(
tank_flow, Qs, fprime=tank_flow_prime,
args=(Qs, a, b, As, ht, C, z, at, Va, Cor, m, Hb),
tol=1e-10)
zp = z + (Qs+QPs)/at
HP = (a - QPs)/b
VP2 = -C2[:,0] + C2[:,1]*HP
VP1 = C1[:,0] - C1[:,1]*HP
if nn == 0: # pipe start
VP = np.float64(VP2)
else: # pipe end
VP = np.float64(VP1)
return HP, VP, QPs, zp