La forme de la surface de l'eau de l'étang de sable est suivie par un calcul de débit inégal. Bien que la forme de la section transversale soit rectangulaire, la forme de la surface de l'eau est calculée lorsque la largeur de la voie navigable et l'élévation du plancher de la voie navigable changent. Puisqu'il s'agit d'un calcul de débit normal, le niveau d'eau est suivi de l'aval vers l'amont.
L'équation de base d'un écoulement inégal à résoudre est la suivante.
\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 $ | Débit (valeur constante) |
| $ x_2, z_2, A_2, R_2, n_2, h_2 $ | Caractéristiques de la section amont (distance, élévation du plancher du chenal, aire de la section de l'eau courante, rayon hydraulique, rugosité de l'équipage Coefficient, profondeur de l'eau) |
| $ x_1, z_1, A_1, R_1, n_1, h_1 $ | Caractéristiques de la section en aval (distance, hauteur du plancher du chenal, surface de la section de l'eau courante, rayon hydraulique, rugosité de l'équipage Coefficient, profondeur de l'eau) |
L'indice 1 signifie le côté aval et l'indice 2 signifie le côté amont.
Comme ce cas est normal, toutes les conditions en aval (indice 1) sont connues et la profondeur d'eau en amont $ h_2 $ est calculée séquentiellement.
Les caractéristiques de la section sont définies en fonction de la distance x '' dans la fonction def sec (x, h) ''.
Puisqu'il existe de nombreux changements transversaux, la définition des caractéristiques transversales (`` def sec (x, h) '') est longue, mais la partie centrale du calcul (dichotomie) est très simple.
def sec(x,h)Est-ce que la distance (x) Et la profondeur de l'eau (h), Et les spécifications de section transversale (z, a, r, n) Est une fonction qui produit, donc si cela est changé, elle peut également être appliquée au calcul de sections transversales arbitraires.
Dans ce cas, la forme de la section transversale est un simple rectangle, mais comme la section transversale d'écoulement d'eau $ A $ et le rayon d'écoulement d'eau $ R $ sont des fonctions de la profondeur d'eau $ h $ comme indiqué ci-dessous, la profondeur d'eau en amont est calculée par la dichotomie. Je cherche.
\begin{equation*}
A=b\cdot h \qquad R=\cfrac{b\cdot h}{b+2\cdot h}
\end{equation*}
Ici, $ b $ est la largeur du chenal et $ h $ est la profondeur de l'eau.
Les deux valeurs initiales de la dichotomie sont dans la fonction def bisection (h1, x1, x2, q) , `ha = 3.0``` and `hb = 7.0``` Il est spécifié comme. Ces valeurs doivent être réécrites de manière appropriée en fonction du problème à résoudre.
Le texte complet du programme est présenté ci-dessous.
# 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()
Le résultat du calcul est le suivant. J'ai été surpris que le niveau d'eau ait chuté de manière inattendue! Vaut la peine de calculer.
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
c'est tout