この節の作者: 清水康行 <yasu@i-ric.org>

1次元非定常流(幅が変化する場合)

基礎式

前章では流下方向に川幅が一定の場合の1次元不定流の計算法について述べましたが， ここでは流下方向川幅が変化する場合について述べます． 幅が一定の場合，流量は単位幅流量 \(q\) で扱いましたが，ここでは全流量 \(Q\) を用います． 一般断面の1次元非定常流れの連続式および運動方程式を以下で表します.

(94)\[{{\partial A}\over{\partial t}}+ {{\partial Q}\over{\partial x}}= 0\]

(95)\[{{1}\over{g A}}{{\partial Q}\over{\partial t}} +{{1}\over{g A}}{{\partial (A u^2)}\over{\partial x}} +{{\partial H}\over{\partial x}} + I_f =0\]

ここで, \(t\) は時間, \(x\) は流下方向距離, \(A\) は断面積, \(Q\) は流量, \(g\) は重力加速度, \(u\) は断面平均流速, \(H\) は水位, \(I_f\) はエネルギー勾配 です. (94) 式は上下流の流量の差が断面積の時間変化に等しいという連続の式 です. (95) 式の第1項以外は定常流の運動量の式です. 第1項は運動量の 時間変化を表したものです.

上記の式は任意の断面形に適用可能な式ですが, 最終的には1次元の河床変動計算に応用しますので以降, 断面は矩形断面と して扱うことにします. 即ち,

(96)\[A=Bh, \ \ Q=Au=Bhu\]

の関係が成立するものとします. ただし, \(B\) は水路幅, \(h\) は水深で, \(H=h+z_b\), \(z_b\) は河床高です.

(96) 式の関係を用い, \(h\) は \(x\) と \(t\) の関数, \(B\) は \(x\) のみの関数 ( \(B\) は時間的には変化しない)であることを考慮すると (94) , (95) 式は以下のように変形されます.

(97)\[B{{\partial h}\over{\partial t}} + {{\partial (uBh)}\over{\partial x}} = 0\]

(98)\[{{\partial (uBh)}\over{\partial t}} + {{\partial (u^2 Bh)}\over{\partial x}}= -gBh{{\partial H}\over{\partial x}} -gBhI_f\]

(98) 式の左辺は保存形と呼ばれる形で表現されており, 水面形が不連続の場合にも適用可能な形式となっています. 水面形が連続な場合, 即ち, \(\displaystyle{{\partial h}\over{\partial x}}\) や \(\displaystyle{{\partial u}\over{\partial x}}\) が 存在する場合には (98) 式の左辺は,

(99)\[Bu{{\partial h}\over{\partial t}}+Bh{{\partial u}\over{\partial t}} + u{{\partial(Buh)}\over{\partial x}}+Buh{{\partial u}\over{\partial x}}\]

と展開可能となり, (97) 式の連続式を用いることにより, (98) 式は次式となります.

(100)\[{{\partial u}\over{\partial t}}+u{{\partial u}\over{\partial x}} = -g{{\partial H}\over{\partial x}} -gI_f\]

エネルギー勾配 \(I_f\) はマニング則を用い, 広矩形断面を仮定して,

(101)\[u={1 \over n_m} h^{2/3} I_f^{1/2}\]

の関係を用いると, (100) 式は以下のように変形できます. ただし \(n_m\) はマニングの粗度係数です.

(102)\[{{\partial u}\over{\partial t}}+u{{\partial u}\over{\partial x}} = -g{{\partial H}\over{\partial x}}-{{g n_m^2 u^2}\over{h^{4/3}}}\]

以下の説明では, 運動方程式としてこの非保存形の (102) 式を用います. 本来は, 非保存形では不連続な流れは扱えないはずなのですが, CIP法を用いると (102) 式でも不連続な流れにある程度対応可能なことは前の章でも説明したとおりです.

CIP法の適用の準備(分離解法)

基礎式を再記すると,

(103)\[{{\partial h}\over{\partial t}}+ {1 \over B}{{\partial (uBh)}\over{\partial x}} = 0\]

(104)\[{{\partial u}\over{\partial t}} + u{{\partial u}\over{\partial x}} = -g{{\partial (h + z_b)}\over{\partial x}} -{{g n_m^2 u^2}\over{h^{4/3}}}\]

ということになります. (103) 式で \(h\) を, (104) 式で \(u\) を計算します. ここで, 分離解法を (104) 式に適用します. 即ち,

(105)\[{{\partial u}\over{\partial t}}= -g \displaystyle{\left( {{\partial h}\over{\partial x}} + {{\partial z_b}\over{\partial x}} \right) } -{{g n_m^2 u^2}\over{h^{4/3}}}\]

(106)\[{{\partial u}\over{\partial t}}+u{{\partial u}\over{\partial x}}=0\]

と分離します. ただし, \(H=h+z_b\) の関係を用いてあります. (104) 式で, \(u^n \rightarrow u^{n+1}\) を計算するのを, 中間に \(\widetilde{u}\) を設け, (105) 式で \(u^n \rightarrow \widetilde{u}\) を, (106) 式で \(\widetilde{u} \rightarrow u^{n+1}\) を計算します. なお (105) 式には \(h\) が含まれているため, (105) 式を計算する際に, (103) 式を同時 に満たす必要があります. これを \(h^n \rightarrow \widetilde{h}\) と 表しておきます. なお (106) 式のほうには \(h\) が含まれませんので, (103) 式と関連付けようがありませんので, \(h^{n+1}\) は (103) 式と (105) 式を連立して出てくる \(\widetilde{h}\) をそのまま用いて,

(107)\[h^{n+1}=\widetilde{h}\]

とすることにします. 以上を整理しますと以下の表のようになります．

Table 4 1次元非定常流れの分離解法

Phase名	\(u\) と \(h\) の更新	使用する式
[A] Non-advection Phase	\(u^n \rightarrow \widetilde{u}, h^n \rightarrow \widetilde{h}\)	(105) 式, (103) 式
[B] Advection Phase	\(\widetilde{u} \rightarrow u^{n+1}, \widetilde{h} \rightarrow h^{n+1}\)	(106) 式, (107) 式

Non-Advection Phase

(105) 式で \(u^n \rightarrow \widetilde{u}\) を計算して 同時に \(h^n \rightarrow \widetilde{h}\) を行うので, (105) 式を以下のように表します.

(108)\[{{\widetilde{u}-u^n}\over{\Delta t}}=-g \displaystyle{ \left( {{\partial \widetilde{h}}\over{\partial x}} + {{\partial z_b}\over{\partial x}} \right) } -{{g n_m^2 \widetilde{u}^2}\over{\widetilde{h}^{4/3}}}\]

これより,

(109)\[\widetilde{u}= u^n - g \Delta t {{\partial \widetilde{h}}\over{\partial x}} - g \Delta t {{\partial z_b}\over{\partial x}} - g \Delta t {{n_m^2 \widetilde{u}^2}\over{\widetilde{h}^{4/3}}}\]

一方, (103) 式で \(h^n \rightarrow \widetilde{h}\) を求めるので,

(110)\[{{\widetilde{h}-h^n}\over{\Delta t}} + {1 \over B} {{\partial}\over{\partial x}} (\widetilde{u} \widetilde{h} B ) =0\]

を満たす必要があります. したがって, (109) 式と (110) 式を連立して得られる \(\widetilde{h}\) および \(\widetilde{u}\) を求める必要があります.

CIP法で流体の計算を行う場合は \(u\) と \(h\) の計算点を交互に配置する スタッカード格子が用いられます. ここでは， Figure 33 に示すような計算点 の配置を採用することにします.

Figure 33 : スタッカード格子

Figure 33 の計算点の配置を考慮して (109) 式を差分表示すると,

(111)\[\widetilde{u}(i) = u^n(i)-g\Delta t {{\widetilde{h}(i+1) - \widetilde{h}(i)}\over{\Delta x}} -g\Delta t {{z_b(i+1) - z_b(i)}\over{\Delta x}} - g \Delta t {{n_m^2 \widetilde{u}(i)^2}\over{\widetilde{h_p}^{4/3}}}\]

ただし, \(\widetilde{h_p}=[\widetilde{h}(i)+\widetilde{h}(i+1)]/2\) , となり, これを (109) 式左辺第2項に代入するために,

(112)\[\widetilde{Q}=\widetilde{u} \widetilde{h} B\]

なる流量 \(Q\) を定義し(\(Q\) 計算点は \(u\) と同じとする), 格子点の配置を考慮して,

(113)\[\widetilde{Q}(i)= \widetilde{u}(i) {{\widetilde{h}(i+1)+\widetilde{h}(i)}\over{2}} B_u(i)\]

ただし, \(B_u\) は \(u\) の計算点における B の値であり，

(114)\[B(i)=\cfrac{B_u(i)+B_u(i-1)}{2}\]

で表される. (113) 式の \(Q\) を用いた (110) 式の差分式

(115)\[\widetilde{h}(i) = h^n(i) + {1 \over B(i)} \left[ \widetilde{Q}(i)-\widetilde{Q}(i-1) \right] {{\Delta t}\over{\Delta x}}\]

より \(\widetilde{h}\) を求めることが出来ます. ただし, \(bsa_7\) 式には右辺の \(\widetilde{Q}`にも :math:\)widetilde{h}` が陰な形式で含まれるため, 繰り返し計算が必要になります. この手順は Figure 34 に示すとおりです.

Figure 34 : Non-Advection Phase における \(\widetilde{h}(i)\) および \(\widetilde{u}(i)\) の計算手順

Non-Advection Phaseにおける \(\widetilde{u}(i)\) が求まった段階で, Non-Advection Phaseにおける \(\displaystyle{\widetilde{{\partial u}\over{\partial x}}}\) を求めておく必要があります. これは, (61) 式で \(f\) を \(u\) と置いた次式で求めることが出来ます.

(116)\[\widetilde{{\partial u}\over{\partial x}}(i)= {{\partial u}\over{\partial x}}^n(i) +{{1}\over{2\Delta t}} \left[ \widetilde{u}(i+1)-u^n(i+1)-\widetilde{u}(i-1)+u^n(i-1) \right]\]

Advection Phase

Advextion Phaseにおいては (106) 式に直接CIP法を適用することが出来ます. これについては何度も説明しましたので結果のみ書きますと, 次式で求めることが出来ます.

(117)\[{u(i)}^{n+1}=U(X)= \left[(aX+b)X+ \widetilde{{\partial u}\over{\partial x}}(i) \right] X +\widetilde{u}(i)\]

ただし, \(\widetilde{u}(i)\) および, \(\displaystyle{\widetilde{{\partial u}\over{\partial x}}(i)}\) はそれぞれNon-Advection Phaseにおける \(u\) および \(u\) の空間微分量であり, それぞれ, (111) 式および (115) 式で求まった値です. また, \(U(X)\) は Non-Advection Phaseにおける \(u\) の空間分布関数であり, 次式で与えられます.

(118)\[U(X)= \left[(aX+b)X+ \widetilde{{\partial u}\over{\partial x}}(i) \right] X +\widetilde{u}(i)\]

従って, Advection Phaseにおける \(u\) は \(U(X)\) を用いて,

(119)\[u(i)^{n+1}=U(X)\]

となります. Advection Phaseにおける \(u\) の空間微分量 も同様に,

(120)\[\widehat{{\partial u}\over{\partial x}}(i) = 3aX^2 + 2bX + \widetilde{{\partial u}\over{\partial x}}(i)\]

および

(121)\[{{\partial u}\over{\partial x}}^{n+1}(i) = \widehat{{\partial u}\over{\partial x}}(i) -{{u(i+1)-u(i-1))}\over{2\Delta x}} \widehat{{\partial u}\over{\partial x}}(i) \Delta t\]

で求めることができます.

まとめ

以上，流下方向に川幅が変化する1次元自由水面流れの計算をまとめると Figure 35 のフローチャートとなります.

Figure 35 : 幅が変化する1次元自由水面を有する流れの計算フローチャート

Pythonコード

上記水路幅が縦断方向に変化する1次元非定常流れの計算コードの例を下記に示します．

https://github.com/YasuShimizu/1D-Free-Water-Surface-with-Width-Variation

このGitHubのフォルダにある拡張子がpyのファイルが計算用のファイル群ですが，このうちmain.pyが計算コードの本体と なります．計算条件は[config.yml]を下記で記述します．コメントにある項目を適宜修正することにより 条件を変更しながら計算を行うことが出来ます．

xl : 200 # 水路長(m)
nx: 51  # x方向格子数
j_channel: 1 # 水路の状態 1=一定勾配 2=勾配変化点あり
# j_channel=1の場合
slope: 0.002 # 全体の勾配
# j_channel=2の場合
x_slope : 100 # 勾配変化点
slope1 : 0.001 # 上流側勾配
slope2 : 0.01 # 下流側勾配
xb1 : 90 # 突起始り地点の上流端からの距離（m)
xb2 : 100 # 突起頂点の上流端からの距離(m)
xb3 : 110 # 突起終了地点の上流からの距離(m)
dbed : 0. #突起高さ(m)
xw1 : 80 # 川幅変化点1の位置(m)
xw2 : 100 # 川幅変化点2の位置(m)
xw3 : 120 # 川幅変化点3の位置(m)
w0: 1 # 川幅1(m)(上流端)
w1: 1 # 川幅2(m)
w2: 0.75 # 川幅3(m)
w3: 1 # 川幅4(m)
w4: 1 # 川幅5(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: 0.7 # 下流端の水深(等流水深に対する割合)
etime : 400 # 計算終了時間(sec)
dt : 0.05 # 計算タイムステップ
tuk : 10 # アウトプットインターバル

川幅の変化は，xw1, xw2, xw3, w0, w1, w2, w3,w4 で与えることとし， これらのパラメーターの説明は下図のとおりとなります．

Figure 36 : 川幅の設定に関するパラメータ

計算例

中央部に狭窄区間を含む直線水路の水面計の計算例

幅1mで中央部に幅が0.75mの狭窄部を含む水路の計算例を下図にしまします． この例では，下流端を等流水深より下げた，低下背水としてあります． 計算結果の図には，以下の式で表される等流水位，限界水位も並記して あります．

(122)\[h_0=\left(\cfrac{Q n}{B \sqrt{i}}\right)^{3/5}, \ \ h_c=\left(\cfrac{Q^2}{B^2 g}\right)^{1/3}\]

Figure 37 : 水路の途中に狭窄区間を含む流れの低下背水の計算例