この節の作者: 清水康行 <yasu@i-ric.org>
本テキストへの感想, 質問, 要望, バグ報告などはこちらへお願いします.
1次元非定常流(幅が一定の場合)¶
基礎式¶
水路幅一定, 長方形断面水路における非定常流れの連続式および運動方程式を以下で表します.
ここで,
(72) 式の左辺は保存形と呼ばれる形で表現されており,
水面形が不連続の場合にも適用可能な形式となっています.
水面形が連続な場合, 即ち,
と展開可能となり, (71) 式の連続式を用いることにより, (72) 式は次式となります. ただし, 拡散項は近似的な表現です.
以下の説明では, 運動方程式としてこの非保存形の (74) 式を用います. 本来は, 非保存形では不連続な流れは扱えないはずなのですが, CIP法を用いるとこの式でも不連続な流れに ある程度対応可能なことが知られています. おそらく,CIP法で使用している3次式の補完関数が, 関数の急変点にもある程度まで対応可能であることによるものと考えられます.
CIP法の適用の準備(分離解法)¶
議論をさらに単純化するために, 河床せん断力(摩擦)が無視でき, 粘性も無視できる 場合を考えます. 基礎式は,
ということになります. 式を見て分かると思いますが今までは (52) 式のように未知
変数が
と分離します. ただし,
とすることにします. 以上を整理しますと下表のようになります.
Phase名 |
使用する式 |
|
---|---|---|
[A] Non-advection Phase |
||
[B] Advection Phase |
Non-Advection Phase¶
(77) 式で
これより,
一方, (75) 式で
を満たす必要があります. したがって,
(81) 式と (82) 式を連立して得られる
CIP法で流体の計算を行う場合は一般に

Figure 24 : スタッカード格子¶
この計算点の配置を考慮して (81) 式を差分表示すると,
となり(
なる単位幅流量
および,
より

Figure 25 : Non-Advection Phase における
この繰り返し計算を通して得られる
Non-Advection Phaseにおける
Advection Phase¶
Advextion Phaseにおいては (78) 式に直接CIP法を適用することが出来ます.
即ち, (62) 式で
ただし,
従って, Advection Phaseにおける
となります. Advection Phaseにおける
および
で求めることができます.
境界条件¶
本計算法では
という方法がありこれは自由境界条件ということになります. なお,水位(水深)については, 計算条件に応じて適宜選択する必要があります.詳細は省略しますが,常流においては 下流端の水位(水深)境界条件を,射流においては上流端の水位(水深)境界条件を与えることが出来ます. なお, 当然のことですが境界条件は Non-Advection PhaseとAdvection Phaseの両方に適用する必要があります.
まとめ¶
以上の1次元自由水面を有する流れの計算をまとめると以下の フローチャートとなります.

Figure 26 : 1次元自由水面を有する流れの計算フローチャート¶
Pythonコード¶
上記1次元非定常流れの計算コードの例を下記に示します.
https://github.com/YasuShimizu/1D-Free-Water-Surface
このGitHubのフォルダにある拡張子がpyのファイルが計算用のファイル群ですが,このうちmain.pyが計算コードの本体で, この中の下記部分で,ymlファイル名をコメントアウトしてある部分を逐次有効にすることにより 異なる条件の水面計計算が可能となります.
with open('config.yml','r', encoding='utf-8') as yml:
#with open('config_jump.yml','r', encoding='utf-8') as yml: #跳水
#with open('config_trans.yml','r', encoding='utf-8') as yml: #遷移流
#with open('config_back_water.yml','r', encoding='utf-8') as yml: #堰上げ背水
#with open('config_drop_down.yml','r', encoding='utf-8') as yml: #低下背水
#with open('config_bump.yml','r', encoding='utf-8') as yml: #突起のある流れ
config = yaml.load(yml)
ymlファイルは計算条件を記述したファイルで,上記GitHubサイトに含まれてます. 一例として[config_trans.yml]を下記に示します.コメントにある項目を適宜修正することにより 条件を変更しながら計算を行うことが出来ます.なお,計算法の説明では簡単のために摩擦を無視した場合で 説明しましたが,計算コードでは運動方程式に摩擦項を追加してあります.ymlファイルで粗度係数をゼロとすれば 摩擦を無視した流れの計算が出来ます.
xl : 200 # 水路長(m)
nx: 51 # x方向格子数
j_channel: 2 # 水路の状態 1=一定勾配 2=勾配変化点あり
# j_channel=1の場合
slope: 0.002 # 全体の勾配
# j_channel=2の場合
x_slope : 130 # 勾配変化点
slope1 : 0.001 # 上流側勾配
slope2 : 0.01 # 下流側勾配
xb1 : 4 # 突起始り位置の上流端からの距離(m)
xb2 : 5 # 突起頂点の上流店からの距離(m)
xb3 : 6 # 突起終了点の上流店からの距離(m)
dbed : 0.00 #突起高さ(m)
qp : 1 #単位幅流量(m**2/2)
g : 9.8 #重力加速度
snm : 0.025 # 粗度係数
alh : 0.6 # 水位計算緩和係数
lmax : 50 # 水位計算最大繰り返し回数
errmax : 0.0001 # 打ち切り誤差
hmin : 0.01 # 最小水深(m)
j_upstm : 2 # 上流端の条件 1=Wall, 2=Free, 3=Fixed Value
j_dwstm : 3 # 下流端の条件
alpha_up: 1.0 # 上流端の水深(等流水深に対する割合)
alpha_dw: 1.0 # 下流端の水深(等流水深に対する割合)
etime : 400 # 計算終了時間(sec)
dt : 0.05 # 計算タイムステップ(sec)
tuk : 10 # アウトプットインターバル(sec)
なお,[config_trans.yml]中の xb1, xb2, xb3 および dbed は下図のような定義とします.

Figure 27 : 河床突起の設定パラメーター¶
計算例¶
堰上げ背水¶
[config_back_water.yml]を使用した計算例を下記に示します.計算結果は,上から水位(Water),河床形状(Bed),流速(Velocity),
単位幅流量(Discharge,フルード数(Fr)の縦断形です.
水位縦断には以下の (93) 式による等流水深

Figure 28 : 堰上げ背水の計算例¶
Figure 28 の計算では,下流端で等水深より20%深い水深に堰上げられた水位の影響が上流へ伝わる様子が計算されています.
低下背水¶
[config_drop_down.yml]を用いた計算結果を Figure 29 に示します.この例では,下流端で等流水深の70%に保たれた下流端水深を 境界条件とる低下背水の様子が計算されています.

Figure 29 : 低下背水の計算例¶
跳水¶
[config_jump.yml]を用いた計算結果を Figure 30 に示します.この例では,上流側が急こう配,下流側が緩勾配になっている水路 における跳水の様子が計算されています.

Figure 30 : 低下背水の計算例¶
遷移流¶
[config_trans.yml]を用いた計算結果を Figure 31 に示します.この例では,上流側が緩勾配,下流側が急こう配になっている水路 で常流から射流に遷移する流れの様子が計算されています.

Figure 31 : 遷移流の計算例¶
水路中央に突起がある流れ¶
[config_bump.yml]を用いた計算結果を Figure 32 に示します.この例では,水路中央部に突起がある流れの様子が計算されています. この例は常流なので,突起部分で水位は低下している様子が計算されています.(水理学で習いましたよね?覚えていますか?)

Figure 32 : 水路中央に突起がある流れ(常流)の計算例¶