はじめに
等流とは,時間的・場所的に断面形状,河床勾配,水深,流量が変化しない流れであり,河川関係では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$を計算する方法が,計算過程が単純である.
実行用シェルスクリプト
実行用スクリプトでは,以下の3つのプログラムを実行している.
- py_eng_uniflow.py : 不等流解析プログラム本体
- py_eng_uni_figHQ.py : 水位ー流量曲線描画
- py_eng_uni_figsec.py : 解析断面描画
実行用スクリプト本文は以下の通り.
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座標(標高)
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列目:水面幅
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)と地形線に囲まれた部分の面積および地形線長さを計算するようにしている.
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))
水位ー流量曲線描画プログラム
これは,標準的な2次元グラフ作成の手順に則っている.
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を用いて水面下を水色で塗りつぶしている.
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()
出力画像
以上