import numpy as np
import pandas as pd
class section(object):
def __init__(self, zb, B, manning, distance=np.nan):
self.distance = distance
self.zb = zb
self.B = B
self.manning = manning
def H2ABS(self, H):
dep = (H - self.zb)
B = self.B
A = dep*B
S = 2.0 * dep + B
return A, B, S
def A2HBS(self, A, Hini=float(-9999)):
B = self.B
dep = A/B
H = dep + self.zb
S = 2.0 * dep + B
return H, B, S
def calIeAlphaBetaRcUsubABS(self, Q, H):
A, B, S = self.H2ABS(H)
K = A**(5/3)/S**(2/3)/self.manning
Ie = Q**2/K**2
Usub = Q/A
Alpha = float(1.0)
Beta = float(1.0)
Rc = A/S
return Ie, Alpha, Beta, Rc, Usub, A, B, S
def calH0ABS(self, Q, ib):
dh = float(0.5)
delta = lambda f : np.abs(ib - f)/f
zb = self.zb
H = zb + float(0.01)
arr = self.calIeAlphaBetaRcUsubABS(Q, H)
ie = arr[0]
while delta(ie) > 10**(-8):
if ib < ie:
H += dh
else :
dh *= float(0.5)
H -= dh
arr = self.calIeAlphaBetaRcUsubABS(Q, H)
ie = arr[0]
A, B, S = arr[-3], arr[-2], arr[-1]
return H, A, B, S
def UnSteadyflow(sections, A, Q, H, Abound, Qbound, dt):
g = float(9.8)
imax = len(A)
Anew, Qnew, Hnew = np.zeros(imax), np.zeros(imax), np.zeros(imax)
ie = np.zeros(imax)
Beta = np.zeros(imax)
# continuous equation
for i in range(1, imax-1) :
dx = 0.5*(sections[i-1].distance - sections[i+1].distance)
Anew[i] = A[i] - dt * ( Q[i] - Q[i-1] ) / dx
Anew[imax-1] = Abound
Anew[0] = (Anew[1] - A[1]) + A[0]
for i in range(imax) :
s = sections[i]
Hnew[i], _, _ = s.A2HBS(Anew[i], H[i])
arr = s.calIeAlphaBetaRcUsubABS(Q[i], H[i])
ie[i] = arr[0]
Beta[i] = arr[2]
# moumentum equation
for i in range(1, imax-1):
ic, im, ip = i, i-1, i+1
dxp = sections[ic].distance - sections[ip].distance
dxm = sections[im].distance - sections[ic].distance
dxc = 0.5*(sections[im].distance - sections[ip].distance)
Cr1 = 0.5*( Q[ic]/A[ic] + Q[ip]/A[ip] )*dt/dxp
Cr2 = 0.5*( Q[ic]/A[ic] + Q[im]/A[im] )*dt/dxm
dHdx1 = ( Hnew[ip] - Hnew[ic] ) / dxp
dHdx2 = ( Hnew[ic] - Hnew[im] ) / dxm
dHdx = (float(1.0) - Cr1) * dHdx1 + Cr2 * dHdx2
Qnew[ic] = Q[ic] - dt * ( Beta[ic]*Q[ic]**2/A[ic] - Beta[im]*Q[im]**2/A[im] ) / dxc \
- dt * g * Anew[ic] * dHdx \
- dt * g * A[ic] * ie[ic]
Qnew[imax-1] = Qnew[imax-2]
Qnew[0] = Qbound
return Anew, Qnew, Hnew
def kinematicwave(sections, A, H, Qbound, dt):
g = float(9.8)
imax = len(A)
Anew, Qnew, Hnew = np.zeros(imax), np.zeros(imax), np.zeros(imax)
# manning
Qnew[0] = Qbound
for i in range(1, imax) :
s = sections[i]
sb = sections[i-1]
dx = sb.distance - s.distance
ib = (sb.zb - s.zb)/dx
Ap, Bp, Sp = s.H2ABS(H[i])
Qnew[i] = (Ap/Sp)**(2/3) * ib**0.5 * Ap/s.manning
Anew[i] = A[i] - dt * ( Qnew[i] - Qnew[i-1] ) / dx
Anew[0] = (Anew[1] - A[1]) + A[0]
for i in range(imax) :
s = sections[i]
Hnew[i], _, _ = s.A2HBS(Anew[i], H[i])
return Anew, Qnew, Hnew
df = pd.read_csv('201910shiroyamaDam.csv', index_col='date', parse_dates=True)
dfp = df['2019/10/12 8:00':]
dfp = dfp.resample('10S').interpolate()
Qbound = dfp['Qout'].values
Qini = Qbound[0]
ib =float(1/1000)
x = np.linspace(0, 20000,101, dtype=np.float64, endpoint=True)
zb = x * ib
zb = zb[::-1]
x = x[::-1]
n = float(0.03)
B = float(100)
river = []
for dis, zbp in zip(x,zb):
s = section(zbp, B, n, dis)
river.append(s)
A = np.zeros_like(river)
H = np.zeros_like(river)
Q = np.full_like(river, Qini)
for i, s in enumerate(river):
H[i], A[i], _, _ = s.calH0ABS(Q[i], ib)
dt = float(10)
Qub = Qini
n = len(Qbound)
i = len(river)
Amatd = np.zeros((n,i))
Qmatd = np.zeros((n,i))
Hmatd = np.zeros((n,i))
Amatd[0] = A
Qmatd[0] = Q
Hmatd[0] = H
for nn in range(1, n):
_, Abound, _, _ = s.calH0ABS(Q[-1], ib)
A, Q, H = UnSteadyflow(river, A, Q, H, Abound, Qbound[nn], dt)
Amatd[nn] = A
Qmatd[nn] = Q
Hmatd[nn] = H
# set initial condition
A = np.zeros_like(river)
H = np.zeros_like(river)
Q = np.full_like(river, Qini)
for i, s in enumerate(river):
H[i], A[i], _, _ = s.calH0ABS(Q[i], ib)
# simulation
dt = float(10)
Qub = Qini
n = len(Qbound)
i = len(river)
Amatk = np.zeros((n,i))
Qmatk = np.zeros((n,i))
Hmatk = np.zeros((n,i))
Amatk[0] = A
Qmatk[0] = Q
Hmatk[0] = H
for nn in range(1, n):
A, Q, H = kinematicwave(river, A, H, Qbound[nn], dt)
Amatk[nn] = A
Qmatk[nn] = Q
Hmatk[nn] = H
import hvplot.pandas
dfoutd = pd.DataFrame({'上流から0km': Qmatd[:,0],
'上流から10km': Qmatd[:,50],
'上流から20km': Qmatd[:,100]}
,index = dfp.index)
dfoutk = pd.DataFrame({'上流から0km': Qmatk[:,0],
'上流から10km': Qmatk[:,50],
'上流から20km': Qmatk[:,100]}
,index = dfp.index)
g = dfoutd.hvplot(ylabel='discharge[m3/s]',title='Dynamic wave', width=600, grid=True) + dfoutk.hvplot(ylabel='discharge[m3/s]',title='Kinematic wave', width=600, grid=True)
g = g.cols(1)
g
hvplot.save(g,'fig.html')