Ungleiche Durchflussberechnung mit Python: Berechnung der Wasseroberflächenform von Wasserkanälen (Sandteichen), bei denen sich die Wasserkanalbreite und die Höhe des Wasserkanalbodens ändern

Verfolgen Sie die Wasseroberflächenform des Sandteichs durch Berechnung der ungleichen Strömung. Obwohl die Querschnittsform rechteckig ist, wird die Wasseroberflächenform berechnet, wenn sich die Wasserstraßenbreite und die Höhe des Wasserstraßenbodens ändern. Da es sich um eine normale Durchflussberechnung handelt, wird der Wasserstand von stromabwärts nach stromaufwärts verfolgt.

Die Grundgleichung des zu lösenden ungleichen Flusses lautet wie folgt.

\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 $ Durchflussrate (konstanter Wert)
$ x_2, z_2, A_2, R_2, n_2, h_2 $ Querschnittsquerschnittseigenschaften (Abstand, Kanalbetthöhe, Fließwasserquerschnittsfläche, hydraulischer Radius, Besatzungsrauheit Koeffizient, Wassertiefe)
$ x_1, z_1, A_1, R_1, n_1, h_1 $ Querschnittseigenschaften (Abstand, Kanalbodenhöhe, Fließwasserquerschnittsfläche, hydraulischer Radius, Rauheit der Besatzung) Koeffizient, Wassertiefe)

Index 1 bedeutet die stromabwärtige Seite und Index 2 bedeutet die stromaufwärtige Seite. Da dieser Fall normal ist, sind alle stromabwärtigen Bedingungen (Index 1) bekannt, und die stromaufwärtige Wassertiefe $ h_2 $ wird nacheinander berechnet. Die Querschnittseigenschaften werden als Funktion des Abstandes "x" in der Funktion "def sec (x, h)" definiert. Da es viele Querschnittsänderungen gibt, ist die Definition der Querschnittseigenschaften (`def sec (x, h)`) lang, aber der Kernteil der Berechnung (Dichotomie) ist sehr einfach. def sec(x,h)Ist die Entfernung (x) Und Wassertiefe (h) Und die Querschnittsspezifikationen (z, a, r, n) Ist eine Funktion, die ausgegeben wird. Wenn dies geändert wird, kann sie auch zur Berechnung beliebiger Querschnitte verwendet werden.

In diesem Fall ist die Querschnittsform ein einfaches Rechteck, aber da die Wasserflussquerschnittsfläche $ A $ und der Wasserflussradius $ R $ Funktionen der Wassertiefe $ h $ sind, wie unten gezeigt, wird die Wassertiefe stromaufwärts durch die Dichotomie berechnet. Ich suche.

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

Hier ist $ b $ die Kanalbreite und $ h $ die Wassertiefe. Die beiden Anfangswerte in der Dichotomie liegen in der Funktion `def Bisection (h1, x1, x2, q)`, ha = 3.0``` und hb = 7.0``` Es ist angegeben als. Diese Werte müssen entsprechend dem zu lösenden Problem entsprechend umgeschrieben werden.

Der vollständige Programmtext wird unten angezeigt.

# 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()        

Das Berechnungsergebnis ist wie folgt. Ich war überrascht, dass der Wasserstand unerwartet gesunken ist! Es lohnt sich zu berechnen.

       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

das ist alles

Recommended Posts

Ungleiche Durchflussberechnung mit Python: Berechnung der Wasseroberflächenform von Wasserkanälen (Sandteichen), bei denen sich die Wasserkanalbreite und die Höhe des Wasserkanalbodens ändern