不等流計算により沈砂池の水面形を追跡する。
断面形状は矩形であるが、水路幅・水路床標高が変化する場合の水面形計算となる。
常流計算であるため、下流から上流に向かって水位を追跡する。

解くべき不等流の基本方程式は以下の通り。

\begin{equation*}
\left(\cfrac{Q^2}{2 g A_2{}^2}+h_2+z_2\right) 
- \left(\cfrac{Q^2}{2 g A_1{}^2}+h_1+z_1\right)
= \cfrac{1}{2}\left(\cfrac{n_2{}^2 Q^2}{R_1{}^{4/3} A_1{}^2} + \cfrac{n_1{}^2 Q^2}{R_2{}^{4/3} A_2{}^2}\right)\cdot(x_2-x_1)
\end{equation*}

$Q$	流量（一定値）
$x_2, z_2, A_2, R_2, n_2, h_2$	上流側断面特性（距離程、水路床標高、流水断面積、動水半径、マニングの粗度係数、水深）
$x_1, z_1, A_1, R_1, n_1, h_1$	下流側断面特性（距離程、水路床標高、流水断面積、動水半径、マニングの粗度係数、水深）

添字1 は下流側を意味し、添字2 は上流側を意味する。
このケースは常流であるため、下流側条件（添字1）はすべて既知であり、上流側の水深 $h_2$ を逐次求めていく。
断面特性は、関数 def sec(x,h) の中で距離程 x の関数として定めている。
断面変化箇所が多いので、断面特性の定義箇所（def sec(x,h)）が長いが、計算の核の部分（二分法）は非常に単純である。
def sec(x,h) は、距離程（x）と水深（h）を入力し、断面諸元（z, A, R, n）を出力する関数なので、ここを変えれば任意断面の計算にも応用することができる。

この事例では断面形は単純な矩形であるが、通水断面積 $A$ および動水半径 $R$ は以下に示すとおり、水深 $h$ の関数であるため、二分法により上流側水深を求めている。

\begin{equation*}
A=b\cdot h \qquad R=\cfrac{b\cdot h}{b+2\cdot h}
\end{equation*}

ここに $b$ は水路幅、$h$ は水深。
二分法における２つの初期値は、関数 def bisection(h1,x1,x2,q) の中で、ha=3.0 および hb=7.0 として指定している。これらの値は、解くべき問題に応じて、適切に書き換える必要がある。

プログラム全文を以下に示す。

# Non-Uniform Flow Analysis (Subcritical flow)
import numpy as np


def sec(x, h):
    # definition of section property
    # x  : distance
    # h  : water depth
    # zz : invert level
    # aa : secion area
    # rr : hydraulic radius
    # nn : Manning's roughness coefficient
    n0=0.014 # roughness coefficient (normal value)
    zz,aa,rr,nn=0,0,0,0
    if 0.0 <= x < 11.0:
        ds=11
        nn=n0
        z1=562.2; b1=4.0
        z2=560.5; b2=26.0
        zz=z1-(z1-z2)/ds*x
        bb=b1+(b2-b1)/ds*x
        aa=bb*h
        rr=bb*h/(bb+2*h)
    if 11.0 <= x < 19.0:
        ds=8.0
        nn=n0
        z1=560.5; b1=12
        z2=560.5; b2=12
        zz=z1-(z1-z2)/ds*(x-11)
        bb=b1+(b2-b1)/ds*(x-11)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 19.0 <= x < 47.0:
        ds=29.0
        nn=n0
        z1=560.5     ; b1=12.0
        z2=z1+ds*0.02; b2=12.0
        zz=z1-(z1-z2)/ds*(x-19)
        bb=b1+(b2-b1)/ds*(x-19)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 47.0 <= x < 62.0:
        ds=15.0
        nn=n0
        z1=561.06    ; b1=12.0
        z2=z1+ds*0.02; b2=6.0
        zz=z1-(z1-z2)/ds*(x-47)
        bb=b1+(b2-b1)/ds*(x-47)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 62.0 <= x < 69.0:
        ds=7.0
        nn=n0
        z1=561.36    ; b1=6.0
        z2=z1+ds*0.02; b2=6.0
        zz=z1-(z1-z2)/ds*(x-62)
        bb=b1+(b2-b1)/ds*(x-62)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 69.0 <= x < 69.5:
        ds=0.5
        nn=n0
        z1=561.5; b1=6.0
        z2=562.0; b2=6.0
        zz=z1-(z1-z2)/ds*(x-69)
        bb=b1+(b2-b1)/ds*(x-69)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 69.5 <= x < 73.5:
        ds=4.0
        nn=n0
        z1=562.0; b1=6.0
        z2=562.0; b2=6.0
        zz=z1-(z1-z2)/ds*(x-69.5)
        bb=b1+(b2-b1)/ds*(x-69.5)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 73.5 <= x < 74.0:
        ds=0.5
        nn=n0
        z1=562.0; b1=6.0
        z2=561.5; b2=6.5
        zz=z1-(z1-z2)/ds*(x-73.5)
        bb=b1+(b2-b1)/ds*(x-73.5)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 74.0 <= x < 77.0:
        ds=3.0
        nn=n0
        z1=561.5; b1=6.5
        z2=561.5; b2=6.5
        zz=z1-(z1-z2)/ds*(x-74.0)
        bb=b1+(b2-b1)/ds*(x-74.0)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 77.0 <= x < 96.625:
        ds=19.625
        nn=n0
        z1=561.5; b1=6.5
        z2=563.0; b2=6.5
        zz=z1-(z1-z2)/ds*(x-77.0)
        bb=b1+(b2-b1)/ds*(x-77.0)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 96.625 <= x <= 102.125:
        ds=5.5
        nn=n0
        z1=563.0; b1=6.5
        z2=563.0; b2=6.5
        zz=z1-(z1-z2)/ds*(x-96.625)
        bb=b1+(b2-b1)/ds*(x-96.625)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    return zz,aa,rr,nn


def func(h2,h1,x1,x2,q):
    g=9.8
    z1,a1,r1,n1=sec(x1,h1)
    z2,a2,r2,n2=sec(x2,h2)
    e2=q**2/2/g/a2**2+h2+z2
    e1=q**2/2/g/a1**2+h1+z1
    e3=0.5*(n1**2*q**2/r1**(4/3)/a1**2 + n2**2*q**2/r2**(4/3)/a2**2)*(x2-x1)
    return e2-e1-e3


def bisection(h1,x1,x2,q):
    ha=3.0 # lower initial value for bisection method
    hb=7.0 # upper initial value for bisection method
    for k in range(100):
        hi=0.5*(ha+hb)
        fa=func(ha,h1,x1,x2,q)
        fb=func(hb,h1,x1,x2,q)
        fi=func(hi,h1,x1,x2,q)
        if fa*fi<0: hb=hi
        if fb*fi<0: ha=hi
        #print(fa,fi,fb)
        if np.abs(hb-ha)<1e-10: break
    return hi


def main():
    g=9.8    # gravity acceleration
    q=42.0   # discharge
    # starting point (sub-critical flow)
    h1=3.9 # water depth at starting point
    x1=0.0   # origin of distance coordinate
    z1,a1,r1,n1=sec(x1,h1) # section property
    v1=q/a1  # flow velocity
    print('{0:>8s}{1:>8s}{2:>8s}{3:>8s}{4:>10s}{5:>8s}'.format('x','z','h','z+h','Energy','v'))
    print('{0:8.3f}{1:8.3f}{2:8.3f}{3:8.3f}{4:10.5f}{5:8.3f}'.format(x1,z1,h1,z1+h1,z1+h1+v1**2/2/g,v1))

    # calculation point
    #xp=np.arange(1,103,1)
    xp=np.array([11,19,47,62,69,69.5,73.5,74,77,96.625,102.125])

    # water level calculation to upstream direction
    for x2 in xp:
        h2=bisection(h1,x1,x2,q)
        z2,a2,r2,n2=sec(x2,h2)
        v2=q/a2
        print('{0:8.3f}{1:8.3f}{2:8.3f}{3:8.3f}{4:10.5f}{5:8.3f}'.format(x2,z2,h2,z2+h2,z2+h2+v2**2/2/g,v2))
        x1=x2 # distance
        h1=h2 # water depth


#==============
# Execution
#==============
if __name__ == '__main__': main()

計算結果は以下の通り。
意外に水面低下が大きいのにびっくり！
計算する価値がある。

       x       z       h     z+h    Energy       v
   0.000 562.200   3.900 566.100 566.46982   2.692
  11.000 560.500   5.971 566.471 566.47522   0.293
  19.000 560.500   5.971 566.471 566.47523   0.293
  47.000 561.060   5.410 566.470 566.47528   0.323
  62.000 561.360   5.091 566.451 566.47541   0.687
  69.000 561.500   4.950 566.450 566.47553   0.707
  69.500 562.000   4.444 566.444 566.47554   0.788
  73.500 562.000   4.444 566.444 566.47562   0.788
  74.000 561.500   4.954 566.454 566.47563   0.652
  77.000 561.500   4.954 566.454 566.47567   0.652
  96.625 563.000   3.431 566.431 566.47615   0.942
 102.125 563.000   3.431 566.431 566.47634   0.942

以　上