Python, matplotlibで開水路の等流解析を行う

はじめに

等流とは，時間的・場所的に断面形状，河床勾配，水深，流量が変化しない流れであり，河川関係ではManningの平均流速公式を用いて解析される場合が多い．
解析に用いるのは，以下の式である．

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

• $v$ : 平均流速
• $n$ : 粗度係数
• $R$ : 動水半径
• $I$ : 河床勾配
• $A$ : 流水断面積
• $S$ : 潤辺
• $Q$ : 流量

流水断面積および潤辺は，水深の関数となるので，最初に水深$h$（もしくは水位標高）を仮定し，その水深での流下量$Q$を計算する方法が，計算過程が単純である．

実行用シェルスクリプト

実行用スクリプトでは，以下の３つのプログラムを実行している．

• py_eng_uniflow.py ： 不等流解析プログラム本体
• py_eng_uni_figHQ.py ： 水位ー流量曲線描画
• py_eng_uni_figsec.py ： 解析断面描画

実行用スクリプト本文は以下の通り．

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

入力データ

1行目：1列目粗度係数$n$，2列目河床勾配$I$
2行目以降：1列目河床x座標（左岸から右岸へ），y座標（標高）

inp_sec1.txt

0.04 0.001
-138 100
-114  90
 -70  80
 -45  70
 -32  62
   0  61.5
  32  62
  60  70
  80  80
  98  84
 120  84

出力データ

出力データ並びは以下の通り．
等流解析のみであれば，流量と標高もしくは水深の関係がわかれば良いが，不等流解析の入力条件作成の目的も兼ね，流水断面積，潤辺，水面幅も出力させている．ちなみに水面幅は限界水深算定に必要な情報である．

• 1列目：流量
• 2列目：水位標高
• 3列目：水深
• 4列目：流水断面積
• 5列目：潤辺
• 6列目：水面幅

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
..........      ..........

等流解析プログラム

プログラム内で，断面積$A$と潤辺$S$の計算を行っている．
これは下図を参照して，左岸（左側）から，入力された地形を表す線分を追っていき，指定した標高ライン（el）と地形線に囲まれた部分の面積および地形線長さを計算するようにしている．

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

水位ー流量曲線描画プログラム

これは，標準的な２次元グラフ作成の手順に則っている．

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

出力画像

計算断面描画プログラム

見かけの問題であるが，fill_betweenを用いて水面下を水色で塗りつぶしている．

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

出力画像

以上