教科書の「不等流計算」は理解できるが、実務で汎用される「準二次元不等流計算」がなかなか理解し難い方に向けて、記事を書きました。
多分相当わかりやすいと思いますので、河川行政に関わる方やコンサルの方に読んで頂きたいです。
全4回の第2回目です。
本記事はGitHub、nbviewer、Colab、Marp(スライド形式)でも公開しています。
※Colabの解説記事はこちら
※Marpの解説記事はこちら
はじめに
- 下図に示す実河川ような横断面形状は一般断面や自然河道断面などと呼ばれる。前回は矩形断面を対象として不等流計算を行ったが今回は一般断面を対象とする。
- このような実務者向けの計算テクニックは水理公式集、河川砂防技術基準などに一部記述はあるのものの参考書籍が少ない。
import matplotlib.pyplot as plt
L=[0, 5,93,100,200,206,294,300]
Z=[6, 3.5, 3.5, 0, 0, 3, 3, 6]
plt.plot(L,Z, c='k')
plt.show()
基礎式
一般断面の不等流計算の基礎式は次式を用いる。
ここに、:流量、:水位、:河積、:コントロールボリュームに作用する力とする。
前項で示した矩形断面の式と比較して以下の変更を行っている。
- エネルギー保存則をベースとした式形としたが、一般断面の計算では摩擦損失以外の抵抗を考慮するため、運動量保存則をベースとした式形とした。それにより、右辺の生成項の式形が変更されている。
- (圧力項+重力項)をに書き換えている。
後者について補足すると、
例えば、下図のように、WL1、WL2と2つの水位が与えらた場合にそれぞれは定義できるだろうか。
それぞれの定義は難しい上には水位によって変わることが理解できる。これが、一般断面の大きな特徴であり、水位によって理論的な(物理的に意味を持つ)河床高が変わる。
の2変数を同時に解くことはできないため、としてのみを変数として計算を行なう。
import matplotlib.pyplot as plt
x=[0, 5,93,100,200,206,294,300]
y=[6, 3.5, 3.5, 0, 0, 3, 3, 6]
plt.plot(x,y,c='k')
plt.axhline(y=5,c='b',label='WL1')
plt.axhline(y=3.3,c='b',ls='dashed',label='WL2')
plt.legend()
plt.show()
マニング則についても以下のとおり若干の変更を加える。
ここに、:流量、:水位、:河積、:流速、:径深、:潤辺、:通水能、:マニングの粗度係数とする。
水位、潤辺、粗度係数は、水位の関数となる。式の簡略化のため、これらを合わせた通水能を用いることがある。
運動方程式の生成項が摩擦損失のみの場合、マニング則を用いると以下のとおりになる。
限界水深、等流水深の定義
一般断面の基礎式を用いて、限界水深、等流水深を矩形断面と同様に定義で設定することは難しいため、便宜的に以下のとおりに設定する。
等流水深
全ての損失による水頭の勾配が河床勾配と釣り合う状態を等流と定義してその水深を等流水深とする。
なお、一般断面では水深を用いないため、正確には等流時の水位(以降、等流水位と定義)となる。
生成項が摩擦損失のみの場合は、
となり、これを満足する水位を反復法などにより求めれば良い。
限界水深
フルード数が1となる水位を限界流時の水位(以降、限界水位と定義)とする。
フルード数は、
となるが、平方根の中に水深が含まれるため、そのままでは計算できない。
そのため、水深の代替として径深を使用することが多い。その他にはを用いることもある。
参考:FORUM8ソフトウェア:等流・不等流の計算・3DCAD Ver.9 Q&A
離散化
離散化は次式となる。
なお、:上流側、:下流側とする。
常流の場合、下流から逐次計算を行なうため未知数はのみとなる。
一般断面の場合、離散式の未知数による微分が難しいため、ニュートン法が使いづらい。そのため、二分法を使用することをおすすめする(もちろんニュートン法を使っても良い)。ただし、通常の二分法では安定的に計算することが難しいため、少し工夫が必要である。この点については次回詳述する。
数値計算方法2:断面特性(河積、潤辺)の計算方法
河道断面の座標を次図のように定義する。
河積,潤辺,径深は次式で計算する。
ここでは、水際位置と座標定義点が同一であるが、実計算ではその限りではないため、水面と河道断面の交点を求めて計算する必要がある。詳細は後述するプログラムを参照いただきたい。
数値計算方法3:(断面内の流速が一定と仮定した)合成粗度係数の計算方法
ここまでは河道断面内で粗度係数が一定であることを前提としたが、上図のとおり、断面内で粗度係数が異なる場合について考える。
各潤辺の粗度の影響が及ぶ範囲を下図のように考えて、各領域の流速がが全て等しいと仮定する(Einsteinの方法)。その場合、断面平均流速となる。断面全体および各潤辺領域でマニング則が成立すると仮定すると、次式が導出される。
出典:椿東一郎 水理学1 pp.147
この方法は平均流速公式1a、は合成粗度係数と呼ばれる。
当然ではあるが水位によって合成粗度係数は変わる。
導出方法
断面全体の流速が各領域の流速と一致すると仮定する。(Einsteinの方法)
マニング則を用いると次式が導かれる。
のため、
となり、合成粗度が導かれる。
断面特性を計算するプログラム例
ある河川断面の座標が次図のとおりとする。
import matplotlib.pyplot as plt
L=[0, 5,93,100,200,206,294,300]
Z=[6, 3.5, 3.5, 0, 0, 3, 3, 6]
plt.plot(L,Z, c='k')
plt.show()
また、粗度係数は次のとおりとする。(個数は測点間で定義するためi-1個)
n=[0.041, 0.041,0.030,0.030,0.030,0.040,0.040]
断面特性を計算するプログラムは次のように記述できる。
def H2ABSKn(l, z, n, H):
A, B, S, SN = float(0), float(0), float(0), float(0)
for i in range(1, len(l)):
dx = l[i] - l[i-1]
dy = z[i] - z[i-1]
hb, hf = H - z[i-1], H - z[i]
if hb <= float(0) :
if hf > float(0) :
dx_dh = dx / (hf - hb)
B += hf * dx_dh
A += 0.5 * hf * hf * dx_dh
Sp = hf * np.sqrt( dx_dh * dx_dh + 1.0)
S += Sp
SN += Sp * n[i-1]**1.5
elif hf <= float(0) :
if hb > float(0) :
dx_dh = dx / (hf - hb)
B -= hb * dx_dh
A -= 0.5 * hb * hb * dx_dh
Sp = hb * np.sqrt(dx_dh * dx_dh + 1.0)
S += Sp
SN += Sp * n[i-1]**1.5
else :
B += dx
A += 0.5 * dx * (hf + hb)
Sp = np.sqrt(dx**2 + dy**2)
S += Sp
SN += Sp * n[i-1]**1.5
if S <= float(0):
nd = float(0)
K = float(0)
else:
nd = (SN/S)**(2.0/3.0)
K = A**(5.0/3.0)/nd/S**(2.0/3.0)
return A, B, S, K, nd
- 水位が5.0の場合、各種断面諸元は次のとおりとなる。
import numpy as np
A, B, S, K, nd = H2ABSKn(np.array(L), np.array(Z), np.array(n), 5.0)
print('河積:{}\n水面幅:{}\n潤辺:{}\n通水能:{}\n合成粗度係数:{}'.format(A, B, S, K, nd))
河積:858.0
水面幅:296.0
潤辺:298.3606797749979
通水能:47342.84520415623
合成粗度係数:0.03664910724429057
不等流計算のプログラム例
Coming soon.
参考文献
リンク
リンク
リンク
関連記事
computational-sediment-hyd.hatenablog.jp