先に紹介した不等流計算プログラム(常流計算)では、下流端協会での水位を指定する必要があり、それに等流水深を用いています。すなわち等流水深を知る必要があるわけです。
一応、非線形方程式を解くことになりますが、このような計算をサッとできるのが、Pythonのいいところだと思います。
矩形断面の等流水深計算のための基本式は以下の通り。
\begin{gather}
Q=A\cdot v \\
v=\cfrac{1}{n}\cdot R^{2/3}\cdot i^{1/2} \\
R=\cfrac{b\cdot h}{b+2\cdot h}
\end{gather}
$Q$ | (既知) 流量 |
$b$ | (既知) 水路幅 |
$n$ | (既知) マニングの粗度係数 |
$i$ | (既知) 水路床勾配 |
$h$ | (未知) 等流水深 |
前述式を $f=0$ の形に変形した以下の方程式を、$h$ について解きます。
\begin{equation}
f=Q-\cfrac{b\cdot h}{n}\cdot \left(\cfrac{b\cdot h}{b+2\cdot h}\right)^{2/3}\cdot i^{1/2}=0
\end{equation}
今回のプログラムでは、scipy.optimize.brentq
を用いて解いています。
Brent法に与える解の初期値は、h1=0, h2=10
として与えています。
限界水深 $h_c$ の計算はおまけです。
プログラムは以下に示すとおり。
# normal depth and critical depth of rectangular cross section
import numpy as np
from scipy import optimize
def cal_hc(q,b):
# critical depth
g=9.8
hc=(q**2/g/b**2)**(1/3)
return hc
def func(h,q,b,n,i):
f=q-b*h/n*(b*h/(b+2*h))**(2/3)*i**(1/2)
return f
def main():
q=42.0 # discharge
b=4.0 # channel width
n=0.014 # Manning's roughness coefficient
i=0.001 # invert gradient
h1=0.0
h2=10.0
hh=optimize.brentq(func,h1,h2,args=(q,b,n,i))
print('hn=',hh) # normal depth
hc=cal_hc(q,b)
print('hc=',hc) # critical depth
#==============
# Execution
#==============
if __name__ == '__main__': main()
計算結果は以下の通り。
hn= 3.866645305835682
hc= 2.2407023732785825
That's all.
Thank you.