Track the water surface shape of the sand basin by unequal flow calculation. Although the cross-sectional shape is rectangular, the water surface shape is calculated when the channel width and channel floor elevation change. Since it is a normal flow calculation, the water level is tracked from downstream to upstream.
The basic equation of inequality to be solved is as follows.
\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 $ | Flow rate (constant value) |
| $ x_2, z_2, A_2, R_2, n_2, h_2 $ | Upstream cross-sectional characteristics (distance, channel floor elevation, running water cross-sectional area, hydraulic radius, roughness of manning) Coefficient, water depth) |
| $ x_1, z_1, A_1, R_1, n_1, h_1 $ | Downstream cross-sectional characteristics (distance, channel floor elevation, running water cross-sectional area, hydraulic radius, roughness of manning) Coefficient, water depth) |
Subscript 1 means the downstream side, and subscript 2 means the upstream side.
Since this case is normal, all downstream conditions (subscript 1) are known, and the upstream water depth $ h_2 $ is calculated sequentially.
The cross-sectional characteristics are defined as a function of the distance `x``` in the function def sec (x, h) . Since there are many cross-section changes, the definition of cross-section characteristics ( def sec (x, h) ``) is long, but the core part of the calculation (dichotomy) is very simple.
def sec(x,h)Is the distance (x) And water depth (h), And the section specifications (z, a, r, n) Is output, so if this is changed, it can be applied to the calculation of arbitrary cross sections.
In this case, the cross-sectional shape is a simple rectangle, but since the water flow cross-sectional area $ A $ and the hydraulic radius $ R $ are functions of the water depth $ h $ as shown below, the upstream water depth is calculated by the dichotomy. I'm looking for.
\begin{equation*}
A=b\cdot h \qquad R=\cfrac{b\cdot h}{b+2\cdot h}
\end{equation*}
Here, $ b $ is the channel width and $ h $ is the water depth.
The two initial values in the dichotomy are in the function def bisection (h1, x1, x2, q)` ``, ha = 3.0 and` `hb = 7.0. It is specified as. These values need to be rewritten appropriately according to the problem to be solved.
The full program text is shown below.
# 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()
The calculation result is as follows. I was surprised that the water level dropped unexpectedly! Worth to calculate.
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
that's all