Introduction

An isocurrent is a flow in which the cross-sectional shape, riverbed gradient, water depth, and flow rate do not change in terms of time and place, and in river relations, it is often analyzed using Manning's average flow velocity formula. The following formula is used for the analysis.

v=\cfrac{1}{n} R^{2/3} I^{1/2} \qquad R=\cfrac{A}{S} \qquad Q=A \cdot v

$ v $: Average flow velocity
$ n $: Roughness length
$ R $: Hydraulic radius
$ I $: Riverbed gradient
$ A $: Running water cross-sectional area
$ S $: Junbe
$ Q $: Flow rate

Since the cross-sectional area of running water and the edge of water are functions of the water depth, the calculation process is simple, assuming the water depth $ h $ (or water level elevation) first and calculating the flow down amount $ Q $ at that water depth. is there.

Shell script for execution

In the execution script, the following three programs are executed.

py_eng_uniflow.py: Inequal flow analysis program body
py_eng_uni_figHQ.py: Water level-flow rate curve drawing
py_eng_uni_figsec.py: Analysis cross section drawing

The text of the execution script is as follows.

`a_py_uni.txt`


python3 py_eng_uniflow.py inp_sec1.txt > out_PSHQ1.txt

python3 py_eng_uni_figHQ.py out_PSHQ1.txt fig_PSHQ1.png

python3 py_eng_uni_figsec.py inp_sec1.txt fig_sec1_PS.png n=0.04,i=0.001 79.4

Input data

1st row: 1st column roughness coefficient $ n $, 2nd column riverbed gradient $ I $ 2nd row and after: 1st column riverbed x coordinate (from left bank to right bank), y coordinate (elevation)

`inp_sec1.txt`

output data

The output data sequence is as follows. For isokinetic analysis only, it is sufficient to know the relationship between flow rate and altitude or water depth, but for the purpose of creating input conditions for unequal flow analysis, the cross-sectional area of flowing water, moist edge, and water surface width are also output. By the way, the water surface width is the information necessary for calculating the limit water depth.

1st row: Flow rate
Second row: Water level elevation
3rd row: Water depth
4th row: running water cross-sectional area
5th row: Junbe
6th row: Water surface width

`out_PSHQ1.txt`


         Q         WL          h          A          S        wsw
     0.000     61.500      0.000      0.000      0.000      0.000
     0.069     61.600      0.100      0.640     12.802     12.800
     0.436     61.700      0.200      2.560     25.603     25.600
     1.285     61.800      0.300      5.760     38.405     38.400
     2.768     61.900      0.400     10.240     51.206     51.200
     5.019     62.000      0.500     16.000     64.008     64.000
     8.760     62.100      0.600     22.426     64.563     64.513
..........      ..........

Isocurrent analysis program

In the program, the cross-sectional area $ A $ and the edge $ S $ are calculated. This is done by following the line segment representing the input terrain from the left bank (left side) with reference to the figure below, and calculating the area and terrain line length of the part surrounded by the specified elevation line (el) and terrain line. I try to do it.

`py_eng_uniflow.py`


import numpy as np
from math import *
import sys
import matplotlib.pyplot as plt

def INTERSEC(x1,y1,x2,y2,el):
    xx=1.0e10
    if 0.0<abs(y2-y1):
        a=(y2-y1)/(x2-x1)
        b=(x2*y1-x1*y2)/(x2-x1)
        xx=(el-b)/a
    return xx

# Main routine
param=sys.argv
fnameR=param[1]

# Input data
x=[]
y=[]
fin=open(fnameR,'r')
text=fin.readline()
text=text.strip()
text=text.split()
rn=float(text[0]) # Manning roughness coefficient
gi=float(text[1]) # Gradient of river bed
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    x=x+[float(text[0])]
    y=y+[float(text[1])]
fin.close()
nn=len(x)

print('{0:>10s} {1:>10s} {2:>10s} {3:>10s} {4:>10s} {5:>10s}'.format('Q','WL','h','A','S','wsw'))
ymin=min(y)
ymax=max(y)
dy=0.1
yy=np.arange(ymin,ymax+dy,dy)
for el in yy:
    S=0.0
    A=0.0
    xxl=0.0
    xxr=0.0
    for i in range(0,nn-1):
        x1=x[i]
        y1=y[i]
        x2=x[i+1]
        y2=y[i+1]
        xx=INTERSEC(x1,y1,x2,y2,el)
        dS=0.0
        dA=0.0
        if y1<el and y2<el:
            dS=sqrt((x2-x1)**2+(y2-y1)**2)
            dA=0.5*(2.0*el-y1-y2)*(x2-x1)
        if x1<=xx and xx<=x2:
            if y2<=el and el<=y1:
                dS=sqrt((x2-xx)**2+(y2-el)**2)
                dA=0.5*(x2-xx)*(el-y2)
                xxl=xx
            if y1<=el and el<=y2:
                dS=sqrt((xx-x1)**2+(el-y1)**2)
                dA=0.5*(xx-x1)*(el-y1)
                xxr=xx
        S=S+dS
        A=A+dA
    if 0.0<S:
        R=A/S
        v=1.0/rn*R**(2/3)*sqrt(gi)
        Q=A*v
    else:
        R=0.0
        v=0.0
        Q=0.0
    h=el-ymin
    wsw=xxr-xxl
    print('{0:10.3f} {1:10.3f} {2:10.3f} {3:10.3f} {4:10.3f} {5:10.3f}'.format(Q,el,h,A,S,wsw))

Water level-flow curve drawing program

This follows the standard two-dimensional graph creation procedure.

`py_eng_uni_figHQ.py`


import numpy as np
import matplotlib.pyplot as plt
import sys

def DINP(fnameR):
    data=np.loadtxt(fnameR,skiprows=1,usecols=(0,1))
    QQ=data[:,0]
    EL=data[:,1]
    return QQ,EL

# Parameter setting
param=sys.argv
fnameR=param[1]     # input filename
fnameF=param[2]     # output image file name

qmin=0.0;qmax=10000.0;dq=1000
emin=60.0;emax=90.0;de=2.0
# Input data
QQ,EL=DINP(fnameR)
# Plot
fig=plt.figure()
plt.plot(QQ,EL,color='black',lw=2.0,label='PS (n=0.04, i=0.001)')
plt.xlabel('Discharge (m$^3$/s)')
plt.ylabel('Elevation (EL.m)')
plt.xlim(qmin,qmax)
plt.ylim(emin,emax)
plt.xticks(np.arange(qmin,qmax+dq,dq))
plt.yticks(np.arange(emin,emax+de,de))
plt.grid()
plt.legend(shadow=True,loc='upper left',handlelength=3)
plt.savefig(fnameF,dpi=200,, bbox_inches="tight", pad_inches=0.2)
plt.show()
plt.close()

Output image

Calculated cross section drawing program

Although it is a matter of appearance, fill_between is used to fill the underwater surface with light blue.

`py_eng_uni_figsec.py`


import numpy as np
#from math import *
import sys
import matplotlib.pyplot as plt

# Main routine
param=sys.argv
fnameR=param[1]     # input filename
fnameF=param[2]     # output image file name
strt=param[3]       # graph title (strings)
fwl=float(param[4]) # target flood water level

# Input data
x=[]
y=[]
fin=open(fnameR,'r')
text=fin.readline()
text=text.strip()
text=text.split()
rn=float(text[0])
gi=float(text[1])
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    x=x+[float(text[0])]
    y=y+[float(text[1])]
fin.close()
xx=np.array(x)
yy=np.array(y)

# plot
xmin=-150
xmax=150.0
ymin=50.0
ymax=100
fig = plt.figure()
ax=fig.add_subplot(111)
ax.plot(xx,yy,color='black',lw=1.0,label='ground surface')
ax.fill_between(xx,yy,fwl,where=yy<=fwl,facecolor='cyan',alpha=0.3,interpolate=True)
ax.text(0,fwl,'10,000 years flood (EL.'+str(fwl)+')',fontsize=10,color='black',ha='center',va='top')
ax.set_title(strt)
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Elevation (EL.m)')
ax.set_xlim(xmin,xmax)
ax.set_ylim(ymin,ymax)
aspect=1.0*(ymax-ymin)/(xmax-xmin)*(ax.get_xlim()[1] - ax.get_xlim()[0]) / (ax.get_ylim()[1] - ax.get_ylim()[0])
ax.set_aspect(aspect)
plt.savefig(fnameF,dpi=200,, bbox_inches="tight", pad_inches=0.2)
plt.show()
plt.close()

Output image

that's all

Perform isocurrent analysis of open channels with Python and matplotlib

Introduction

Shell script for execution

a_py_uni.txt

Input data

inp_sec1.txt

output data

out_PSHQ1.txt

Isocurrent analysis program

py_eng_uniflow.py

Water level-flow curve drawing program

py_eng_uni_figHQ.py

Output image

Calculated cross section drawing program

py_eng_uni_figsec.py

Output image

`a_py_uni.txt`

`inp_sec1.txt`

`out_PSHQ1.txt`

`py_eng_uniflow.py`

`py_eng_uni_figHQ.py`

`py_eng_uni_figsec.py`