Non-uniform flow analysis by Python: calculation of water surface frofile of settling basin
A program to calculate water surface profile of settling basin which has an open channel with varied width and varied invert level. Since the condition is subcritical flow, the calculation will be carried out from downstream section to upstream section.
Basic equation of non-uniform flow analysis is shown below.
Discharge (constant value) | |
Properties of upstream section (distance, invert level, flow section area, hydraulic radius, Manning's roughness coefficient, water depth) | |
Properties of downstream section (distance, invert level, flow section area, hydraulic radius, Manning's roughness coefficient, water depth) |
Where, subscripu '1' is for downstream section and subscript '2' is for upstream section.
The equations for calculations of section area and hydraulic radius for rectangular cross section are shown below.
where, means channel width, means water depth.
Since this case is a condition under subcritical flow, the conditions of downstream with subscript of '1' are all known values, and calculation is carried out to obtain the water depth of using the know conditions.
The section properties are defined as a function of x
in def sec(x,h)
.
Since the calculation model has many varied section, the part of definition of section properties (def sec(x,h)
) is long.
However, the core part for bi-section method is simple.
Two required initial values for bi-section method are defined as ha=3.0
and hb=7.0
in the function of def bisection(h1,x1,x2,q)
. These initial values should be changed appropriately, depending on the problem to be solved.
A function def sec(x,h)
has a fuction to output section properties (z, A, R, n
) from input data of distance (x
) and water depth (h
). When this function def sec(x,h)
is changed, this program can be applied to a non-uniform flow analysis of open channel with various cross section.
A program 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()
Calculated results are shown below. A water level drop is larger than expected.
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. Thank you.