論文名：河床波上の掃流砂運動及び乱流構造の解明とその河川数値解析モデルへの応用に関する研究
Title : Study on Interaction between Turbulence Structures and Bedload Transportand Its Application to River Morph-dynamic Models
木村 一郎 1
Ichiro KIMURA
1

北海道大学 工学研究院 准教授

要旨

デューンなどの河床波は，流れの抵抗則に影響を及ぼし，その形成過程の予測は治水上極めて重要である．河床波の発生に伴い，河床変動を支配する掃流
砂運動の形態や，河床上乱流構造が変化することが知られているが，その詳細な特性は十分解明されているとは言い難い．特に，デューンクレスト背後に生じる剥離渦のような乱流構造と，粒子の掃流機構については，未解明な点が多い．そこで，本研究では河床波上の掃流砂運動と乱流構造を実験的，数値解析的に検討し，その基本的メカニズムの解明を試みるものである．特に，次の３つの点に着目して検討を実施いした．
① 砂堆の鉛直 2 次元数値解析による再現計算を対象に，平衡流砂量式の適用性に関して２，３の考察を実施し，適切なモデルを同定する．
② DEM-URANS カップリングによるラグランジュ的アプローチにより，デューンの形成過程を数値解析的に検討する．
③ 洪水時のような非定常流下でのデューンの遷移過程について，水理実験により検討を行う．
《キーワード：河床波；砂堆；DEM；URANS；数値解析》

-1-

１．緒

言

デューンなどの河床波は，流れの抵抗則に影響を及ぼし，その形成過程の予測は治水上極めて重要である．河床波の発生に伴い，河床変動を支配する掃流砂運動の形態や，河床上乱流構造が変化する
ことが知られているが，その詳細な特性は未だ解明されていない．本研究では河床波上の掃流砂運動と乱流構造を実験的，数値解析的に検討し，その基本的メカニズムの解明を試みるものである．まず，第２章では，砂堆の鉛直 2 次元数値解析による再現計算を対象に，平衡流砂量式の適用性に関して２，３の考察を実施し，適切なモデルを道程する．次に，第３章では DEM-URANS カップリングによるラグランジュ的アプローチにより，デューンの形成過程を数値解析的に検討する．最後に，第４章では
洪水時のデューンの遷移過程について，水理実験により検討を行う．

２．平衡流砂量式による dune 河床の再現計算
２．１

はじめに

河床に見られる dune は，流れと流砂の相互作用によって発生することが知られている．従来，その形成機構に関する多くの研究が行われてきた．特に，古くより盛んに行われた線形安定解析 1),2) ,3),4) ,5) ,6)の成果によって，dune の発生機構が理論的に明らかにされてきた．中でも，流砂の河床勾配の影響に着目した研究は数多く，例えば Fredsøe2)や山口・泉 4)は，いわゆる平衡流砂量式に河床勾配が考慮された式を用いて卓越波長に対する河床勾配の影響の重要性を明らかにしている．これに対して，中川・
辻本 3)は流砂の遅れによる非平衡性に着目している．彼らは確率モデルで表される非平衡流砂量式を用い，流砂の非平衡性が dune の発生に重要であることを明らかにしている．また，山口・泉 6)は河床勾配の影響が考慮された流砂量式および確率モデルをそれぞれ用いて比較し，卓越波長に対する河床勾配の影響と非平衡性の影響が類似していることを示唆している．
理論的研究と並行して，近年，数値解析による dune の再現計算が行われるようになった．音田・細田 7)や Giri and Shimizu8)は，流砂量式に中川・辻本 3)の用いた確率モデルを適用することによって dune
河床を再現した．さらに Giri ら 9)は非定常流量下での計算を行い，確率モデルにおける砂粒の step
length を無次元掃流力の関数で表すことで dune―平坦床の遷移過程を再現した．確率モデルによるこれらの再現計算では，砂粒の step length による流砂の非平衡性に起因した dune 形成過程が合理的に再現されている．また，別の観点から関根 10)は土粒子の saltation モデルによって dune 形成過程の再現を試みるとともに，砂粒の step length についても議論している．
しかしながら，これらのモデルには河床勾配の影響が考慮されていないため，これらの再現計算では前述のような流砂の河床勾配の影響が検討されていない．また，確率モデルにおける粒子の step
length は，関根 10)による研究の他にもいくつかの研究によって明らかにされつつある 3),11),12)ものの未表‑2.1

Guy ら 17)による実験データ（1 ケースを抜粋）

実験条件

実験
No.

平均粒径
(mm)

Run35

0.93

実験結果

平均勾配

流量
(m2/s)

0.0013

0.089

河床形状

波長
(cm)

波高
(cm)

水深
(cm)

Fr

有効
掃流力

Dune

88.4

1.8

16.2

0.44

0.085

-2-

だ不確定な部分の多いパラメータである．そのため，上述の再現計算では dune 形成に対する流砂の非平衡性の影響が詳細には検討されていない．
本研究では，dune 形成に対する流砂の河床勾配の影響を明らかにするために，河床勾配の影響が考慮された平衡流砂量式を用いた dune の再現計算を行う．また，確率モデルによる非平衡流砂量式を用いた再現計算と比較することによって，dune の形成過程や平衡形状に対する河床勾配の影響および非平衡性の影響について検討する．
表‑2.2

計算 No.

流砂量式

Case‑a

芦田・道上

‑b

２．２

計算条件

Kovacs & Parker

‑c
‑d

計算条件

〃
確率モデル

‑e

〃

‑f

〃

鉛直二次元移動床モデル(1) 支配方程式
本研究では，Giri and Shimizu8)や Giri ら 9)の構築した鉛直二次元流れにおける移動床モデルを適用す
る．流れの支配方程式は次の連続式および運動方程式で表される．

∂u ∂v
=0
+
∂x ∂y

(2-1)

(

)

(

)

(2-2)

∂v
∂v
∂v
1 ∂p ∂∂+− v' u ' +− v' v' − g+u +v = −∂t∂x∂yρ ∂y ∂y
∂y

(2-3)

∂
∂u
∂u
1 ∂p ∂
∂u
+u
+v
=−
+
− u'u' +− u ' v'∂t∂x∂y
ρ ∂x ∂x
∂y

(

)

(

)

ここで， および はそれぞれ鉛直および水平方向の座標， および はそれぞれ鉛直および水平方向の流速， は圧力である．また，，

およびはレイノルズテンソルである．上式は水面と底面の境界面に応じた境界適合座標系に座標変換される．水面では運動学的条件によって水面変動を計算している．底面では no slip 条件となっている．無次元掃流力の算定には次式の対数則を適用し
ている．
(2-4)
ここで，

および

はそれぞれ底面付近の格子点における流速および河床からの距離， は底面摩擦

速度， はカルマン定数であり，である．ここで は粗度高さであり，粒径の 2.5 倍の値を-3-用いる．無次元掃流力 は次式で表される．
(2-5)
ここで

は砂の水中比重， は砂粒の粒径である．乱流モデルは木村・細田

13)による非線形 k-εモデルを適用する．
Giri ら 9)は，確率過程モデルによる非平衡流砂量式を用いている．確率モデルによる河床変動は次式で表される．

(2-6)
ここで，
数（

は河床高さ， は砂の空隙率，

および

）
，

および

および

はそれぞれ砂粒の二次元および三次元の形状係はそれぞれ砂の pick-up rate および deposit rate であり，それぞれ次式で表される．
(2-7)

(2-8)
ここで， は次のような step length の確率密度関数である．

(2-9)
ここで，

は砂粒の平均 step length である．Giri ら 9)は掃流砂に加えて浮遊砂による河床変動量を加味しているが，ここでは掃流砂による変動のみを検討する．
(2) 河床勾配の影響を考慮した平衡流砂量式
本研究では，流砂量式に河床勾配の影響が考慮された平衡流砂量式を適用し，dune の再現計算を試みる．Kovacs & Parker14)は芦田・道上 15)の平衡流砂量式に対して物理的に河床勾配の影響を組み込んだ流砂量式を導いている．本研究では，次式で表される芦田・道上式および Kovacs & Parker 式を適用する．
芦田・道上式;(2-10)
Kovacs & Parker 式;(2-11)
ここで，であり， は流砂量である．また， は粒子の無次元限界掃流力であり，

ここでは岩垣の式 16)により求めている．両式中の係数 は，河床近傍の流速とせん断力の関係を表す-4-

図‑2.1(a) 芦田・道上式による再現結果(Case‑a,

図‑2.1(b) Kovacs and Parker 式による再現結果).

(Case‑b,).

係数であり，本研究では式(4)に対応して次式で表される．

(2-12)
式(10)および(11)中のはクーロン動摩擦係数であり，Kovacs & Parker14)に従いものとする．ここで は砂の安息角である．式(11)ではで表されるの値が小さいほど河床勾配の影響が大きく評価されることになる．また，この河床勾配の項を無視すると，Kovacs & Parker 式(11)は芦田・道上
式(10)に一致する．
(3) 計算条件
本研究では，Guy ら 17)による実験の再現を行った． Guy らの実験データのうち再現計算を行った
実験ケースについて表‑2.1 に示す．また，計算は Case‑a~f まで行い，各ケースに用いた流砂量式およびの値または平均 step length の値を表‑1.2 に示す．計算はいずれも平坦床から開始し，180 分後まで行った．計算初期の平坦床には最大振幅 0.01mm の乱数によって微小な乱れを与えた．上下流端は周期境界条件としている．流下方向の格子間隔は 1cm とし，水深方向の格子間隔は Giri & Shimizu8)に従い底面付近ほど指数的に小さくなるように設定し，最も底面近傍の格子間隔は水深の 0.01 倍とした．
a）河床勾配の影響を表す
芦田・道上式では
たように
から
および

は

の値についての値に関係なく河床勾配の影響はゼロである．Kovacs & Parker 式では，前述しの値が小さいほど河床勾配の影響が大きく現れる．砂の安息角が程度であること

程度の値となることが知られており，本研究では表‑2.2 に示したようにを適用している．

-5-

図‑2.1(d) 確率モデルによる再現結果

図‑2.1(c) Kovacs and Parker 式による再現結果(Case‑c,(Case‑d,)).

).

Kovacs & Parker 式では，流れ方向の河床勾配が安息角になると流砂量が無限大になるため河床勾配
が安息角にごく近い値まで大きくなると何らかの数値的な処理が必要となるが，本研究における計算
内ではそこまで河床勾配が増大することはなかった．
b）砂粒の平均 step length について
中川・辻本 3)は，平坦床における砂粒の平均 step length は 40~250d 程度であるとしている．Giri ら
9)の再現計算でも同様の値が適用されている．これに対して Nelson ら 11)は，最近の実験結果より，平均 step length が掃流力（河床せん断力）と粒子の限界掃流力の関係によって表わされることを示している．彼らの実験によると，掃流力が限界掃流力の 2 倍の時は 15d となり，3 倍の時は 30d，4 倍の時は 40d〜50d となることが示されている．
表‑2.1 に示した実験砂では，岩垣の式 16)によって限界掃流力を求めると 0.034 となり，実験における有効掃流力はその 3 倍程度となっている．Nelson らの実験結果に従うと平均 step length は 30d 程度の値になると考えられる．本研究では，表‑2.2 に示すように妥当だと考えられる平均 step length の値
を Case‑e に，それよりも小さい値および大きい値をそれぞれ Case‑d および f に設定している．
ここで用いる確率過程モデルでは河床勾配の影響を考慮していないが，河床勾配が安息角を超えた
場合は斜面崩壊の計算によって安息角となるようにしている．確率過程モデルを用いた Case‑d f で
は安息角を

２．３

とした．

結果および考察(1) 平衡流砂量式による再現結果
図‑2.1(a)，(b)および(c)にそれぞれ Case‑a，b および c の計算結果を示す．図にはそれぞれ計算開始から 30 分，60 分，90 分，135 分および 180 分後の水面および河床形状を示した．図‑2.1(a)を見て-6-

図‑2.1(e) 確率モデルによる再現結果(Case1‑e,

図‑2.1(f) 確率モデルによる再現結果(Case1‑f,).

).

わかるように，芦田・道上式を用いた場合には，計算開始直後に現れた細かな擾乱が減衰することも発達することもなく，単に下流側へ移動するだけで dune 形状の再現はできない．これに対して図‑2.1(b)および(c)に示したように，Kovacs and Parker 式を用いた場合には擾乱の発達とともに波長が徐々に長くなりながら dune 形状が形成されている．このことから，平衡流砂量式を用いた場合には河床勾配の影響を考慮しないと dune 形状は再現できないことがわかる．
平坦床において Kovacs and Parker 式は芦田・道上式と一致することからも理解できるように，Case‑b
および Case‑c のごく初期においても図‑2.1(a)に見られるような細かい擾乱が現われている．図‑2.1(b)および(c)では，時間の経過とともにこれらの細かな擾乱や比較的波長の短い擾乱が消滅していくのがわかる．このことから，河床勾配の影響には波長の短い擾乱を減衰させる効果のあることがわかる．この効果によって，比較的長い波長の dune が発達し得るものと考えられる．
また，図‑2.1(b)と(c)を比較すると，の値が小さい(c)の方がより波長の長い dune が現われている．これより河床勾配の影響が大きい（の値が小さい）と dune の波長が長くなることがわかる．

これは，河床勾配の影響によって波長の短い擾乱が減衰するためと考えられる．
図‑2.1(a)~(c)を見ると，dune の発達していない(a)の水深が最も大きい．これは，芦田・道上式では河床勾配の影響が考慮されておらず局所的に勾配が安息角を超えて大きくなり得るために擾乱が小
さくても大きな形状抵抗が生じていると考えられる．また，(c)より(b)の水深が大きいのも，Kovacs
and Parker 式において設定した安息角の違いによって(c)より(b)の方が局所的な勾配が大きくなり，
形状抵抗が大きくなっていると考えられる．
(2) 確率モデルによる再現結果
図‑2.1(d)，(e)および(f)にそれぞれ Case‑d，e および f の計算結果を示す．図‑2.1(d)を見るとわかるように，平均 step length が極端に小さい Case‑d の場合は，図‑2.1(a)の芦田・道上式を用いた場-7-表‑2.3

計算結果

計算結果

実験結果

計算 No.

波長
(cm)

波高
(cm)

水深
(cm)

Fr

Case‑a

---

---

22.9

0.26

‑b

50

2.9

15.9

0.45

‑c

100

2.3

14.8

0.50

‑d

---

---

(19.5)(0.33)

‑e

100

4.4

20.3

0.31

‑f

(75)

(2.3)

(15.6)(0.46)

波長
(cm)

波高
(cm)

水深
(cm)

Fr

88.4

1.8

16.2

0.44

合と同様に細かい擾乱が卓越して発生することになる．これに対して図‑2.1(e)を見てわかるように，
平均 step length に適当な値を用いた Case‑e では比較的波長の長い擾乱がゆっくりと発生し，この擾乱
が徐々に発達して dune 形状が現われている．また，平均 step length を Case‑e よりも長く 40d とすると，図‑2.1(f)に見られるように比較的波長の長い擾乱が時間をかけてゆっくり発生するものの，計算時間内に平衡形状には至らなかった．このことから，砂粒の step length に適当な値を用いなければ dune
形状は再現できないことがわかる．
(３) 河床勾配の影響が考慮された平衡流砂量式および確率モデルによる再現結果の比較
表‑2.3 に Case‑a~f の計算結果より得られた dune の平均波長，平均波高，平均水深およびフルード
数 Fr を示す．ここで表に示した値は計算の最終時刻（）における値であり，この時点で平衡状態に達していないケースは参考値として括弧をつけて示した．表より波長については，Kovacs
and Parker 式および確率モデルのどちらの場合でもや平均 step length に適当な値を与えると実験結果をある程度良好に再現できることがわかる．
しかしながら図‑2.1(a)~(f)に示した河床形状をよく見ると，発達した dune 形状に見られる特徴は，
Kovacs and Parker 式を用いた場合と確率モデルを用いた場合とで大きく異なることがわかる．Kovacs
and Parker 式を用いた場合には，dune のクレスト部がほぼ平坦な形状であるのが特徴であり，波高も比較的低くなっている．これに対して確率モデルを用いた場合には，クレストの先端が尖っているのが特徴であり，波高も比較的高くなっている．実験で再現された dune の詳細な形状について不明なため，どちらの特徴が実際に現れているかは検討できないが，表‑2.3 を見てわかるように水深（またはFr）も大きく異なることから，形状に見られる特徴の違いに対応して流れの形状抵抗にも違いが現わ
れていることがわかる．Kovacs and Parker 式を用いた場合には形状抵抗が小さいために実験結果より水深が低く（Fr は高く），確率モデルを用いた場合には形状抵抗が大きいために実験結果より水深が高く（Fr は低く）計算されていることがわかる．
(４) 流砂の河床勾配の影響および非平衡性の影響
a）発生初期および発達過程に見られる波長の変化
Kovacs and Parker 式を用いた場合（図‑2.1(b)~(c)）は，初期に現れた擾乱のうち波長の短い擾乱
が重力の影響を受けて減衰することによって，比較的長い波長の dune が発達していく．これに対して確率モデルを用いた場合（図‑2.1(d)~(f)）は，step length の長さに応じた波長の擾乱が発生し，この擾乱が成長して dune 形状へと発達していく．このとき step length が長いほど発生する波長が長くなる
ため，結果として，dune の波長に対する河床勾配の影響と非平衡性の影響が類似することになる．

-8-b）平衡形状
前節で述べたように Kovacs and Parker 式を用いた場合の平衡形状は，dune のクレスト部がほぼ平坦な形状となり，形状抵抗が確率モデルを用いた場合よりも小さくなっている．これは，河床勾配の影
響に dune の発達を抑制する効果があるためと考えられる．擾乱がある程度発達して局所勾配の影響が強くなると，それ以降の発達が抑制され dune の平衡形状が決定される．ここで Kovacs and Parker 式を用いた場合の形状抵抗が実験よりも小さくなる要因の一つとして，局所勾配の影響が実際よりも大きく評価されているために dune の発達が実際よりも強く抑制されていることが考えられる．山口・泉
6)

は，Kovacs and Parker 式に圧力勾配の影響を加味した線形安定解析および弱非線形安定解析を行い，dune 河床上では圧力勾配の影響で局所勾配の影響が低減されることを示唆している．このことからも実際には dune 背後の剥離が勾配の影響を低減させている可能性がある．
一方，前述のように確率モデルを用いた場合は波高の大きい dune が発達し，形状抵抗が実験結果よりも大きくなる．これは，確率モデルで表される非平衡性には dune 形状を大きく発達させる効果があるためと考えられる．本研究では Giri ら 9)に従い平均 step length の値を流下方向に一様に与えている．
しかし，関根 10)は dune の発達過程において step length は空間的にも時間的にも変動することを示唆しており，本来ならば形状の変化に伴う空間的および時間的な step length の変化を考慮すべきであろう．
前出のように局所的な勾配の影響に dune の発達を抑える効果があるとすれば，本研究で用いた確率モ
デルにはそのような発達を抑制する効果が欠けている可能性がある．また，その逆で Kovacs and Parker
式を用いた場合には，確率モデルで表されるような流砂固有の非平衡性の効果が欠けているために，
クレスト部が平坦な形状しか現われず形状抵抗が小さくなっているとも考えられる．
２．４．結論本研究では，dune の形成過程に対する河床勾配の影響を明らかにするために，河床勾配の影響が考慮された平衡流砂量式による dune の再現計算を行った．従来の再現計算で用いられた確率モデルによる再現計算と比較した結果，dune の波長や平衡形状に対する河床勾配の影響および非平衡性の影響について次のことが明らかとなった．
・河床勾配の影響が考慮された平衡流砂量式を用いても比較的良好に dune 形状が再現される．
・現れる dune の波長に対する河床勾配の影響と非平衡性の影響は類似している．
・流砂量式によって平衡形状の特徴が大きく異なる．これは，河床勾配の影響には dune の発達を抑制する効果があるのに対して，非平衡性の影響には dune 形状を大きく発達させる効果があるためであると考えられる．
２．５

参考文献

1) Engelund, F.: Instability of erodible beds, J. Fluid Mech., Vol.42, pp.225-244, 1970.
2) Fredsøe, J.: On the development of dunes in erodible channels, J. Fluid Mech., Vol.64, pp.1-16, 1974.
3) 中川博次，辻本哲郎: 砂礫の運動に伴う移動床砂面の擾乱発生過程，土木学会論文報告集，No.291,
pp.53-62, 1979.
4) 山口里実, 泉 典洋: デューン−平坦床遷移過程に見られる亜臨界分岐現象，土木学会論文集，
No.740/II-64, pp.75-94, 2003.
-9-5) Colombini, M.:

Revisiting the linear theory of sand dune formation, J. Fluid Mech., Vol.502, pp.1-16,

doi:10.1017/ S0022112003007201, 2004.
6) 山口里実, 泉 典洋: デューンの不安定現象に対する流砂の非平衡性の影響，水工学論文集, Vol.51,
pp.1015-1020, 2007.
7) 音田慎一郎, 細田 尚: 水深積分モデルによる小規模河床波の発生・発達過程と流れの抵抗則の数値解
析，水工学論文集, Vol.48, pp.973-978, 2004.
8) Giri, S., Shimizu, Y.:

Numerical computation of sand dune migration with free surface flow, Water Resources

Research, Vol.42, w10422, doi:10.1029/ 2005WR004588, 2006.
9) Giri, S., Yamaguchi, S., Shimizu, Y., Nelson, J.: Simulating temporal response of bedform characteristics tovarying flow, RCEM 2007, pp.939-947, 2007.
10)関根正人: 土砂粒子の運動の解析を基礎とした河床波の形成過程シミュレーションの試み，土木学会
論文集，No.691/II-57, pp.85-92, 2001.
11)関根正人，吉川秀夫: 掃流砂の停止機構に関する研究，土木学会論文集，No.399/II-10, pp.105-112, 1988.
12)Nelson, J., Shimizu, Y., Giri, S., Shreve, R., McLean, S., Logan, B., Kinzel, P.: Bedform response to Flow
Variability, MARID 2008, keynote, 2008.
13)Kimura, I., Hosoda, T., : A nonlinear k-ε model with realizability for prediction of flows around bluff bodies,Int. J. Num. Meth. Fluids. Vol.42, pp.813-837, 2003.
14)Kovacs, A., Parker, G.: A new vectorial bedload formulation and its application to the time evolution ofstraight river channels, J. Fluid Mech., Vol.267, pp.153-183, 1994.
15)芦田和男，道上正規: 移動床流れの抵抗と流砂量に関する基礎的研究，土木学会論文報告集，No.208,
pp.59-69, 1972.
16)岩垣雄一: 限界掃流力の流体力学的研究，土木学会論文集，No.41, pp. 1-21, 1956.
17)Guy, H. P., Simons, D. B., Richardson, E. V.: Summary of alluvial channel data from flume experiments,
1956~1961, Geological Survey Professional Paper, 462-I, 1966.

- 10 -

３．ＤＥＭ−ＵＲＡＮＳカップリングアプローチによる砂堆の発生過程の再現計算
３．１

はじめに

砂河床移動水路においては河床材料特性と水理条件によって時々刻々と河床形態が姿を変える．河
床波は水理的抵抗として作用し，水位の上昇やそれに伴う河川の氾濫を引き起こす要因と成りえる．
抵抗変化機構の解明と予測手法の確立は，上記の防災的目的にとどまらず，流送土砂量の推定精度の向上という観点からも重要であり，これまでに数多くの研究が行われてきた．これらの河床形態に関する研究は流砂をフラックスとして捉える Euler 的視点に基づくものと，粒子個々の運動に着目する
Lagrange 的視点に基づくそれに大別される．Nakagawa and Tujimoto 1)は Einstein 型確率モデルを基に非平衡流砂量式を提案している．この式は解析的モデルゆえに簡潔な記述が導入されており，粒子間の相互作用や砂粒子の輸送距離に対する河床形状の影響は考慮されていない．Euler 的な再現モデルとしては Giri・Shimizu ら 2)が Kimura ら 3)の非線形ｋ−εモデルと Nakagawa and Tujimoto による流砂式を用いた数値計算によって Dune の形成から平坦床への遷移，及び再形成といった非定常の現象再現に成功している．Lagrange 的な再現に関しても，後藤 4)らは個別要素法を用いて微小河床凹凸発生機構の再現に成功している．その際，河床凹凸が微小であることから乱れの生じない一様流を仮定して計算を実行しているが，実河川の河床では砂粒子が移動することによって流れに変化が生じ，それが粒子に影響を及ぼすというフィードバックが発生している．
本研究の目的は以上の現象を踏まえ，流れを二次元鉛直 URANS モデル，粒子の挙動を個別要素法
によって計算し，Dune の発生・発達過程の再現を試みることである．また，形成された Dune 形状や水位・流速分布を既知の知見と比較することによって DEM-URANS カップリングによる Dune 形成数
値模擬の精度検討も同時に行う．
３．２

計算手法

（１）固液混相流モデルの概要
固液二相流モデルには Euler 型・Lagrange 型が存在する．液相とのアナロジーを前提とする方法を
Euler-Euler カップリング，固相粒子の追跡を行う方法を Euler-Lagrange カップリングと呼ぶ．本研究の 目 的 は個別 要 素 法を用 い て Dune と 流 れ の関 係 を 再現す る こ とであ る が 故に， 本 研 究で はEuler-Lagrange 型のモデルを採用する．
（２）乱流場の解析モデルの概要
掃流砂の挙動を再現するには乱流場を精度よく再現する手法を欠かすことが出来ない．このため，
本研究では Giri・Shimizu ら 2)と同様のフレームワークで計算を行い，時々刻々と変化する河床上の流れを計算した．基礎となるのは Kimura and Hosoda 1)の非線形ｋ−ε鉛直二次元乱流解析モデルである．
前述のように Giri・Shimizu ら 2)はこのモデルを用いて，デューンのクレスト部背後に生じる剥離現象などを再現すると共に，Nakagawa and Tujimoto1)による非平衡流砂モデルを用いたシミュレーションにより Dune の発生・発達過程や，流量変化にともなう遷移河床への移行過程などを良好に再現してみせた．流体の基礎方程式を以下に示す．

∂u ∂v
=0
+
∂x ∂y

(3-1)
- 11 -

(

)

(

)

(3-2)

∂v
∂v
∂v
1 ∂p ∂∂+u +v = −+− v' u ' +− v' v' − g∂t∂x∂yρ ∂y ∂y
∂y

(3-3)

∂
∂u
∂u
∂u
1 ∂p ∂+u+v=−+− u'u' +− u ' v'∂t∂x∂y
ρ ∂x ∂x
∂y

(

)

(

)

ここに x：主流方向座標，y：鉛直方向座標，u：主流方向座標，v：鉛直方向座標，p：圧力，− u ' u ',−u ' v',−v' v'：
レイノルズ応力の各成分である．これらの式をΣ座標に変換して用いる．圧力項に関しては静水圧と動水圧に分離して次のように表す．
H

p = p0 + p' = g ∫ρdz + p' = pg (H − y + p')

(3-4)

y

p’：動水分圧，H：水位であり，動水位については SOR 法による収束計算によりステップ毎に求める．
水位の変化については，次の運動学的条件を用いて計算する．

∂H∂H= v , H = yb + h
+u
∂x
∂t

(3-5)

yb：河床高，h：水深である．レイノルズ応力は非線形ｋ−εモデルを用いて次のように評価する．

− u 'i u ' j = νt S ij −S1ij =

3k2
1kδij − νt ∑ Cβ Sβij − Sβααδij 33
ε β=1 ∂U i ∂U j∂U r ∂U r
1  ∂U ∂U j ∂U r ∂U i 
S 3ij =+
S 2ij =  r∂x j ∂xr 
2  ∂xi ∂xr∂x r ∂x r∂xi ∂x j(3-6)

ここに，νt：渦動粘性係数，δij：クロネッカーのデルタ，ｋ：乱れエネルギー，ε：乱れエネルギ
ーの拡散率，Cβ：モデル係数である．渦動粘性係数は次式で与える．

νt = Cμ

k2
0. 3


,
, Cμ = min  0.09,
2 
1 + 0.09Φ εΦ = max[S ,Ω] , S =S ij =

k 1
k 1
S ij S ij , Ω =ΩijΩijε 2
ε 2∂U i ∂U j∂U i ∂U j−+, Ω ij =
∂xi
∂xi
∂x j
∂x j

(3-7)

ｋおよびεについては次の輸送方程式より求める．

∂k∂k∂k ∂  vt ∂k  ∂  vt ∂k 
+ 
 + Ph −ε
+u
+v
= 
∂t
∂x
∂y ∂x σk ∂x  ∂y σk ∂y - 12 -(3-8)

∂ε ∂ε ∂ε ∂  vt ∂ε ∂  vt ∂εεε2 + 
 + C1ε Ph − C2ε
+u
+v
= 
k
k
∂t
∂x
∂y ∂x σk ∂x  ∂y σk ∂y 

(3-9)

モデル定数として，σk =1.0，σ3=1.3，C1ε=1.44，C2ε=1.92 を用いた．生成項 Ph は次式で与える．

Ph = −u i ' u j '

∂u i
∂x j

(3-10)

壁面の境界条件は壁関数で与え，底面摩擦については粗面の対数則で評価する．このときの粗度高さ
ks は粒径の 2.5 倍とする．水面付近の乱れエネルギー減衰を表す為に次の関数を渦動粘性係数に乗ず
る．




(h − y )εs 
f s = 1 − exp − B
,
3


k s2



B = 10(3-11)

なお，移流項は CIP 法を用いて計算を実行する．
（３）粒子の運動方程式
砂粒子番号 i の併進の運動方程式は，

σA3d 3pi

1= CDρA2 d 3pi u − u pi (u − u pi )
2
dt

du pi

 du du pi −
du pi dτ 
3 du
3  du2 ν  dτ + 6ρA2 d pi
+ρA3d pi
+ CMρA3d pi  −
dτ
πT∫0
dt
dt
dt
T −τ


T

(3-12)+ A3d 3pi (σ−ρ)g + FLM + FLS + FpINTi

であり，また回転運動の方程式は

I

dωpi
dt

= TpINTi(3-13)と表記できる．ρ：流体密度，σ：粒子密度，CM：付加質量係数，CD：抗力係数，A2：砂粒子の２次元形状係数(=4/π)，A3：砂粒子の 3 次元形状係数(=6/π)，d：粒径，u：周囲流体流速，up：砂粒子速度，ν：動粘性係数，τ：底面せん断力，g：重力加速度である．運動方程式の左辺第一項は砂粒子に作用する慣性力，右辺第一項は砂粒子に作用する流体抗力，右辺第二項は周囲流体の加速に伴う力，
右辺第三項は付加質量の加速に必要な力，右辺第四項（Basset 項）は砂粒子と周囲流体の非定常な相対運動履歴に関する力，右辺第五項は重力と浮力で，右辺第六・七項（Magnus 力，Saffman 力）は粒子に作用する揚力である．
これらの運動方程式において Basset 項，Magnus 力，Saffman 力は他の項と比較してオーダーが小さく，近似的には無視できると報告されている．現時点ではこれらの項を考慮していないが，研究の進展によってはこれらの項の導入も検討する．
（４）個別要素法の概要
粒子濃度が増加し，個々の粒子が多数の接点で近傍の粒子と接触する状態に至ると古典的剛体球モ- 13 -

デルは全く機能せず，粒子間の重なりを許容するモデルが必要となる．個別要素法は粒子の重なりの程度に応じて粒子間反発力を与える事によって，粒子間の作用力を決定するモデルである．個々の粒子間作用力は粒子間の接点に弾性スプリング，粘性ダッシュポッド，スライダーという機械的要素を導入してモデル化され，これらの機械的要素は接平面内とその法線方向に配置される．粒子 i，j の接触判定は以下の式を用いて行われる．

Lij =di + d j(3-14)

2

ここに，Lij：粒子 i，j の中心間距離であり，Lij =(x − x ) + (y − y )
2

i

2

j

i

(3-15)

jと書ける．接触が確認された粒子 i の時間Δt 間における変位および角変位であるΔxi，Δyi，ΔΦi を接線方向に設定したローカル座標系(ξ，η)に変換することにより，粒子間の法線・接線方向の相対変位を求める．

0   di / 2 Δξij Δxi −Δx j   0
Δη  = [TGL ] Δy −Δy  + ΔΦ ΔΦ  d / 2
ij 
j
i
j j



 i

(3-16)

ここに，TGL:グローバル座標・ローカル座標間の座標変換行列

 cosαij
− sinαij[TGL ] = 

sinαij = −sinαij 
cosαij 

yi − y j
Lij; cosαij = −(3-17)

xi − x j
Lijであり，作用力 Fξ,Fηを以下のように表現する．

Fξ(t ) = en (t ) + d n (t )Fη(t ) = es (t ) + d s (t )en (t ) = en (t −Δt ) + k nΔξij

es (t ) = es (t − Δt ) + k sΔηij

d n (t ) = cnΔξijd s (t ) = cs

Δt

(3-18)Δηij
Δt

ここに，kn，ks：法線・接線方向バネ係数，cn，cs:：法線・接線方向粘性定数，μ：静止摩擦係数である．これらの係数の決定方法は後藤ら 5)の研究を参考にされたい．法線方向に配置された引張りに抵抗しないジョイントは圧縮を正として，

Fξ(t ) = Fη(t ) = 0( en < 0)(3-19)- 14 -

と表現され，接線方向に活動するジョイントは，Fη(t ) = μ・sign[en (t ), es (t )] ( es > μen )(3-20)と表される，ここに Sign(a,b)は a の絶対値に b の符号を乗じたものと定義する．このローカル座標で導かれた粒子間作用力を以下の式でグローバル座標に変換する．

 FpINTxi di−1  Fξ F
 = −∑ [TGL ]  F  , T pINTi = − ∑ [Fη]2 jj pINTyi 
ij  η  ij
ij

(3-21)

上記の方程式により粒子 j が粒子 i に加える力を導き，粒子 i に加えられるこれらの力を合算することにより粒子の運動を記述する．

0.1

y (m)

0.05
0
-0.05
0

0.1

0.2

0.3

0.4

x (m)

図‑3.1 packing 後の河床
（５）河床の概要
数値模擬を正確に実行する際には実河床を精度よく再現することが欠かせない．関根ら 6)の実験によると，河床面上に露出している河床粒子中心位置は平均河床まわりに分散 d/3 の正規分布に従って存在することが確認されている．河床形成には様々な手法が存在するが，本研究では個別要素法のコードを直接採用することができる重力法の適用が妥当と判断した．重力法とは計算領域に粒子をラン
ダムに発生させ，粒子に重力を作用させることにより自然に河床を形成させる手法である．今ら 7)は
DEM を用いた重力法により分散 d/3 の正規分布を持つ河床を形成することが可能であると報告している．図‑3.1 は packing 後の河床の様子を表したものであるが，微小な起伏の存在を確認できる．
（６）計算領域
鉛直方向 0.21m，主流方向 0.4m の計算領域を設け，x 方向 40，y 方向 88 のグリッドで区切り計算を実行している．初期河床勾配は 1/1000 であり，粒径 0.002m(2mm)の粒子 2000 個を DEM の計算対象とする．Dune の発生条件としてはフルード数に纏わる条件 Fr<0.8 と砂粒レイノルズ数 Re*>20 が良く知られているが，初期条件では上記要件を満すべくτ*=0.15 なる水理条件を設定した．τ*は以下に示す
無次元掃流力（u*：摩擦速度，τ0：底面せん断力）である．

τ* ≡

τ0
u*2
, u* ≡
ρ(σ /ρ− 1) gd(3-22)- 15 -

図‑3.2 Dune の発生過程（７）Euler-Lagrange coupling
Lagrange 的粒子計算と Euler 的水理計算のカップリングは各グリッドにおける間隙率によって河床を定義することにより実行される．即ち，2.6 で述べたグリッドで仕切られた区域を占める粒子の割合
が 70%を越えると，そのグリッドを河床と認識し，その河床分布を液相の計算に反映させる．ところ
で，実際の砂河床では底面流速が 0 となる河床の上部に砂粒子が小跳躍を繰り返す saltation 層が存在するが，上述の手法では saltation 層を定量的に表現することが出来ない．そこで，以下に示す芦田・
道上の式を用いて掃流砂量と無次元掃流力の関係付けを試みる．
- 16 -

qB * ≡ τ 
τ = 17τ*2 1 − *c 1 − *c τ* (σ/ρ− 1)gd 3 τ* 
3

qB(3-23)

ここに，qB:掃流砂量，qB*:無次元掃流砂量，τ*c：無次元限界掃流力である．粒子の間隙率により河床を定義した後に，掃流砂量と無次元掃流力の関係が適切となるように河床高さを調節することで数値模擬と実河川のマクロな整合性を確保した．

３．３

結果と考察

（１）Dune の発生過程
Dune 発生過程を時系列に表示したものを図‑3.2 に示す．左上から右下に向かって 0.1sec おきの河床高さを表す．高さを表す縦軸は一枡 0.02m であり，横軸は流下方向を表す．図‑3.1 で示したように河床は初期条件において d/3 程度の分散を持っている．0sec から 31sec に至る河床では，後藤ら 4)が粒子配列の偏りをもって説明した微小な起伏が存在するのみで周期的な波形は認められない．32sec に幾分スケールの大きい河床波が偶発的に形成されると，この凸部が時間経過と共に自己増幅し，33.5sec 以降には 0.02m 程度の波高を持つに至る．この河床波はこれ以上成長することなく，主流方向に移動することとなる．本研究の平均水深は約 0.1m であるので，再現された河床波の波高は水深の約 1/5 となる．一般的に Dune の鉛直スケール波高は平均水深の 1/3~1/4 とされているので，今回再現された Dune スケールは実物より幾分小さい．波長に関しては，水深の約 5 倍程度の値をとるとされているが，本研究の再現結果では 0.4m の波長を示している．これは平均水深の 4 倍程度なので，波高と同様に若干規模が小さいことが伺える．

0.3

0.2

河床高 ・水 位(m)

河床

0.2

水位

0.1

流下方向流速

0.15

0.05
0.1
0
0

0.04

0.08

0.12

0.16

0.2

0.24

0.28

0.32

0.36

-0.05

0.4

流 下方 向流速 (m/s)

0.25

0.15

0.05

0

流下方向距離(m)図‑3.3 Dune 上の水位・流速・河床高分布（２）Dune 上の河床流速・水位分布
Dune 周囲の水位と流速について関して図‑3.3 のような結果を得た．計算領域における最小河床高- 17 -

さを 0 としてグラフを作成している．流下方向流速とは河床における x 方向流速成分であり，本研究で用いた計算コードを固定平坦床に適用すると河床の流下方向流速は 0.07(m/s)程度の値をとる．一般
的に Dune の上流勾配部では流速が増大し，頂上から下流勾配部では流速の値が減衰，または負になることが知られている．数値模擬の結果を見ると，河床が順勾配の際には流速が平坦河床の 2 倍程度に増幅し，逆勾配へ変化すると共に平坦河床流速の 1/2 程度へと減衰していることが認められる．ゆえに，DEM-URANS カップリングは良好な再現結果を示しているといえるだろう．数値模擬では実現
象において存在するであろう負の流速が存在しないが，これは Dune の波高が実現象よりも小さく流れに十分な剥離が発生していない為と考えられる．Dune 上の水面形状に関して，河床波と水面が逆位相を示すことが知られている．本研究の数値模擬の結果では x=0.35（河床波頂上付近）から周期境界条件を跨いで x=0.12 地点まで上述の理論が成り立つことが確認される．x=0.12~0.32 地点では水位の低下が確認されないが，これは計算領域に大きな Dune 頂上(x=0.35)と小さな Dune 頂上(x=0.11)が混在
し，両 Dune 中間区域において河床勾配の正負が明確に区別されない為と考えられる．
（３）DEM-URANS における Dune 発生機構の考察
上記の結果を踏まえて，DEM-URANS カップリングにおける Dune 発生機構について考察を加える．
程度の差こそあれ，河床には粒子の配列偏差に起因する起伏が常に存在する．この起伏が偶発的に規模の大きなものとなると凸部背後に流速の減衰が生ずる．粒子速度と比較して流体速度が相対的に遅い場合には粒子に進行方向と逆向きの抗力が働き，Dune 下流勾配部に存在する粒子の進行を妨げる．
この抵抗により凸部移動が抑制され，そこに掃流砂が打ち付けられることにより微小凹凸が自己増幅すると考えられる．Dune の成長と共に上流勾配部で流体の加速が誘発されるが，この加速により河床
波頂上部の粒子が甚だしい流速に曝され，進行方向に大きな抗力が発生する．この抗力により河床波の成長は抑制され，Dune の形状が安定すると考えられる．
（４）本研究を踏まえての改善点
一般に Dune の頂上背部には安息角程度の勾配が現れるとされているが，図‑3.3 の Dune 頂上背部
勾配は安息角と比較してやや緩やかな傾斜を示している．これは DEM の砂粒子を球体とみなす仮定に対して，実際の砂粒子は不規則な表面形状をもち，球体よりも噛み合わせが卓越していることが原因と考えられる．加えて，湿潤砂粒子は粘性を帯びることにより，単純な球の衝突とは異なる挙動を示す．ゆえに，今後は仮想粒子間の作用力を粗面湿潤砂粒子の挙動に見合うように改善し，Dune 頂上
背後の傾斜角増大を図る．
３．４．おわりに本研究では，DEM-URANS カップリングアプローチを用いて Dune 発生過程の再現と水理的状況について検討を行った．これまでにも Euler 的な Dune の発生や Lagrange 的な河床微小凹凸発生機構を再現した研究は行われているが，本研究では DEM と URANS を組み合わせて Dune の発生の再現を行うという独自の試みを行っている．
Lagrange 的な Dune の再現を試みることによってクレスト部背後の粒子の挙動や遷移状態における粒子輸送，混合粒系 Dune 発生過程における sorting など，Euler 的な手法では表現することの出来ない現象を対象とした数値模擬を行うことが出来る可能性がある．Dune の再現性を十分に高めた後に，上述した現象に取り組む予定である．

- 18 -

３．５

参考文献

1) Nakagawa, H. and Tsujimoto, T.: Sand bed instability due to bed load motion, Proc. ASCE 106, 2029-2051,
HY12, 1980
2) Giri, S. and Simizu, Y.: Numerical computation of sand dune migration with free surface flow, Water
Resources Research, Vol.42, w10422,doi:10.1029/2005 WR004588, 2006
3) Kimura, I. and Hosoda, T.: A non-linear k-εmodel with reliability for prediction of flows around bluff
bodies, International Journal for Numerical Methods in Fluids, Int. J. Numer. Meth. Fluids 2003; 42:813-837(DOI:10.1002/fld.540).
4) 後藤仁志・原田英治・酒井哲郎
集

第４５巻

水工学論文

2001

5) 後藤仁志・原田英治・酒井哲郎
化

移動床における砂粒子間干渉による砂面攪乱の発生

個別要素法に基づく移動床シミュレーションのモデル定数の最適

土木学会論文集 No.691/Ⅱ-57 2001

6) 関根正人・吉川秀夫

掃流砂の停止機構に関する研究

7) 今翔平・清水康行・木村一郎・山口里美
プ長に関する研究

土木学会論文集

第 399 号/Ⅱ-10 1988

DEM-RANS カップリングによる河床波上掃流砂ステッ

応用力学論文集 Vol.12 2000.

- 19 -

４．非定常流れにおける砂堆の消滅過程に関する実験的研究
４．１

はじめに

小規模河床形態は流れの抵抗の支配的要因であり，河川水位に大きな影響を与える．そのため河川
工学上及び防災上，小規模河床形態の形成・消失機構を明らかにすることは重要である．特に河床波の形成・消失といった遷移現象は洪水時に現れるため，洪水時に見られるような流量の非定常性が河床形状の変化に対してどのような影響を及ぼすのかを明らかにする必要がある．さらに，近年では地球規模気候変動による降雨パターンの変化が指摘されている．これは流量ハイドログラフ特性の変化をもたらし，これまで蓄積してきた河床抵抗則の知見の適合性が下がる懸念もある．このため，河川
防災上，非定常流量下における河床抵抗及び水位の予測は重要課題となっている．
河床波は流れや流砂，河床形状の相互作用によって発生する現象であることが知られている．しかし，河川のように時間的な流量変化に対する河床形状の応答に関してはいくつかの実験が行われている 1) ものの，河床波の形成・消失といった遷移現象に関しては定量的にはほとんど解明されていない．
なかでも，遷移過程において同流量時であっても流量増加過程と減少過程では異なる河床形状が現れる現象が報告されている 2, 3)．流量増加過程では河床波が形成されるため抵抗および水位は高い値を示すのに対して，減少過程では河床波が形成されないため抵抗および水位は低い値を示すといった現象である．
このような非定常流量下における河床形状の遷移現象を予測することは容易ではなく，これまでにも非定常流れの河床波の変形過程についていくつかの実験が行われてきた 1, 2, 4, 5)．これらの実験により，流量変化に対する河床形状変化の時間的な遅れの影響がある程度明らかになっているものの，未だ定量的には明らかになっていない．近年，移動床数値モデルによって非定常流量下における河床
dune の遷移現象を再現し，流れの抵抗の変動過程を定量的に予測しようという試みがなされている（音田・細田 6)，Giri ら 7, 8)）．このような移動床数値モデルによる抵抗予測の精度向上という観点からも，より現実的な条件下における河床形態の遷移現象の再現実験および詳細なデータの蓄積は重要である．
dune から平坦遷移過程に着目した実験は，澤井 4)，山口ら 5)，外山ら 2) が行なっており，河床波の変形過程や水理量が再現されている．しかし，澤井の行った実験では流量の代わりに水路勾配を時間的に変化させ，擬似的な非定常流れで河床波の変形過程を再現している．図‑4.2 実験に用いた非定常流量ハイドログラフ（実績）山口らや外山らの実験では水深に対する水路幅が小さく，流量ハイドロ
グラフが階段状であり，側壁抵抗や急激な流量変化による攪乱が河床形態に及ぼす影響を無視できな
いと考えられる．そこで本研究では，平均河床勾配（および水面勾配）一定および水路幅水深比 B/h =5 − 6 の条件において，急激な流量変化による撹乱の影響を排除するために滑らかな流量ハイドログ
ラフの非定常流量下における移動床実験を行った．また，岸・黒木 9) は河床形態に対する抵抗則を径深・粒径比 R/d によって分類しており，遷移条件が R/d によって異なることを説明している．彼らによると，R/d > 500 の領域では遷移条件や遷移過程における抵抗値がある程度明確なのに対して R/d <500 では実験値のばらつきが大きいというように，R/d = 500 を境に遷移過程における傾向が異なることも示唆されている．本研究でも R/d が実験スケールを表わす重要なパラメータの一つと考えており，山口らや外山らの実験では R/d が 100 程度と比較的小さい値であったのに対して，本研究では- 20 -

R/d > 500 の条件における遷移過程の再現実験を試みた．

図‑4.1 実験水路の概要

図‑4.2 実験に用いた流量のハイドログラフ
表‑4.1 実験における測定項目- 21 -

４．２

実験の概要

図‐4.1 に実験に用いた水路の概要を示す．長さ 25m，幅 100cm のガラス製側壁を有する水路を使用し，発生条件の水理量に合わせるため，表面を塗装した木製側壁で幅 70cm に仕切り，測定部とした．河床勾配は 0.0015 で設定し，河床材料は平均粒径 0.18mm のほぼ均一粒径を用いた．水路の上下流端にはそれぞれ長さ 2.4m 及び 2m，高さ 30cm の固定床を設置した．さらに上下流端の給水部，放流部それぞれ 1m を除いた，水路中央部の 18.6m を移動床測定区間とした．給砂は河床の平均勾配が変化しないよう，固定床と移動床の境界に注意を払って上流端側から行なった．図‐4.2 に実験で用いた流量ハイドログラフを示す．ハイドログラフは渡邊ら 3) により示されている次式で与えられる波
形とした．

(4-1)
ここで，D0：初期水深 h0 で無次元化した水深(D0 =h/h0)，t ：ハイドログラフの継続時間 T で無次元化した時間（t = t/T），α, β, γ, δ は定数.α, β, γ, δ は，北海道の河川の実績洪水波形（昭和 56 年
8 月豊平川など）に適合するように設定し，α=0.5，β=0.125，γ=1.2，δ=1.1 とした．300 分，60 分
流量ハイドログラフともにこの値を用いた．

図‑4.3 実験水路の河床形態区分
表‐4.1 に本実験での測定項目を示す．予備実験の結果から，初期平坦河床に流量 40l/s を通水し，dune が安定的に形成された時点を初期条件として，ピーク流量を 60l/s とする，300 分及び 60 分の2 種類のハイドログラフにより dune の消失過程及び河床抵抗の変動過程を測定した．流量ハイドロ
グラフの形成には 1 分毎に流量調節が可能な装置を用い，流量の急激な変化による攪乱の河床形態への影響を抑えた．また，河床形状の時系列変化過程を解析するため，上流側から 9m− 19m の区間で
- 22 -

デジタルカメラを用いて 1 分毎に測定区間の垂直方向，水平方向から画像撮影を行なった．9m 及び19m 地点に超音波水位計を横断方向にそれぞれ 3 台設置し，水位の測定を行なった．基本的に，結果諸量の算出にはこの水位を使用した．実験期間中の水面勾配を一定に保つ目的で，実験水路上流端から 5m，9m，19m，23.5m の位置にゲージを設置し，1 分毎の水位を目視により補完測定した．
４．３

実験結果および考察（１）河床形態区分
図‐4.3 に実験条件の最小流量 40l/s，中間流量 50l/s，最大流量 60l/s における河床形態区分を示す．
縦軸は径深粒径比 R/d，横軸は無次元掃流力τ*である．流量の設定については，河床材料の粒径 d =0.18mm に対して，芦田・道上及び Garde-Raju の河床形態区分で dune 発生域となり，かつ実験水路の最大給水能力 180l/s の範囲内で R/d >500 を満たす組み合わせとして，河床勾配 ib=0.0015，流量 Q= 40 − 60l/s を設定した．この条件で平坦河床に実際に通水したところ，60l/s では平坦河床のままであり，50l/s では ripple が卓越し，40l/s では明確な dune の安定形成が見られた．このため，dune の安定形成と非定常流量ハイドログラフの初期値を 40l/s とし，平坦河床となる 60l/s を最大流量として設定した．
実験条件は，芦田・道上の河床形態区分では Lower Regime，GardeRaju の河床形態の区分 10) では
Transition 領域に属した．また，実験条件下のフルード数は概ね Fr = 0.55 − 0.7 で常流であった．粒子
レイノルズ数は Re = 6 前後であり，ripple の発生領域であるが，本実験の条件では ripple の波長は54mm〜270mm（粒径の 300〜1500 倍）となり，本実験で観測された dune の波長 1.2m〜1.5m とは明らかに異なるため ripple と dune の区別は明確であった．

図‑4.4 dune 形成過程における河床抵抗係数の時系列変化(2) dune 形成の再現性について
dune 形成の再現性については，定常流れによる dune 形成の予備実験を 2 回行ない，河床抵抗係数の時系列変化により検証した．この目的は，外観上は dune が時間とともに発達して定常状態に達し
安定していく様子が観測されたものの，定性的な評価にとどまるため，河床抵抗係数という具体的な- 23 -

数値を用いて再現性を定量的に評価，確認するためであった．河床抵抗係数 Cf は以下の式により求められる 11)．

(4-2)
(4-3)
ただし，u ：摩擦速度，u：平均流速，R：径深，ie：エネルギー勾配で，ここではエネルギー勾配の代わりに水面勾配を用いて算出した．図‐4.4 に河床抵抗係数 Cf の時系列変化を示す．Cf の算出に
は超音波水位計で計測した水位から求めた水面勾配を用いた．超音波水位計は水面付近のノイズ的な波を観測し，算出された水面勾配及び河床抵抗係数 Cf にばらつきが見られたことから，Cf の移動平均により評価した．図中の線は Cf の移動平均線を表しており，移動平均 MA は下式により定義され
る．

(4-4)
ただし，Pm, Pm−1, Pm−n+1 は直近 n 個の平均値であり，ここでは n = 10 とした．予備実験 1 回目，2 回
目ともに平坦河床から dune が発達していくに従い Cf が上昇した．途中，両ケースで河床抵抗係数
Cf の上昇過程に時間差があるものの，両ケースとも 180 分を過ぎた辺りで Cf の上昇が緩やかに止まり，180 − 210 分にかけて Cf = 0.007 − 0.009 付近で収束する挙動を示した．このことから，通水開始
後 180 − 210 分程度の時間で dune の発達が定常状態に達し安定形成されていると判断し，通水開始
後 210 分（3 時間 30 分）を今回の dune 消失過程に関する水理実験の非定常ハイドログラフの開始時点とした．

図‑4.5 最大波高の時系列変化（300 分ケース）- 24 -

図‑4.6 最大波高の時系列変化（60 分ケース）

図‑4.7

dune の移動状況（t=55 から 70min まで 3 分間隔）

- 25 -

図‑4.8 dune が消失し平坦になった状況（t=150min）(3) 河床形態の時系列変化
dune が消失していく過程を，dune の最大波高で評価した．最大波高は，各時点で最も高い頂部と最も低い底部の差とした．
図‐4.5 及び図‐4.6 に初期条件の dune が消失していく過程における最大波高の変遷を示す．図には流量波形の上昇期，ピーク期，減少期を示した．初期条件とした dune は ripple との混在が見られ，
dune 自体の波高は明瞭に測定できなかった．300 分ケースでは，非定常開始後 30 分頃までは dune とripple がやや混在し，40 分過ぎには ripple が消えて明確な dune が観測されるようになった．60 分
ケースでは，非定常開始後 15 分頃までは dune と ripple がやや混在し，20 分過ぎには ripple が消えて明確な dune が観測されるようになった．両ケースとも，流量増加とともに ripple が消失して dune
形状が明確に観測される過程で，dune 波高が増大する現象が見られた．
300 分ハイドログラフのケースでは，流量ピーク期直前から波高の増大が止まり，その後時間をかけて徐々に波高が低下して最終的に平坦床へと遷移していく様子が観測された．平坦床へと遷移する
間も，水路上流端では新たな dune が発生し続けていた．水路上の dune が流下していく過程で波高が低下すると同時に，上流端で新たに発生する dune の波高も低下していくというように，水路上の duneの波高自体が時間とともに徐々に低下して，120 分付近で dune が完全に消失し水路全体が平坦床へと遷移する過程が観測された．この間，波高の低下が観測されたのに対して，波長については大きな変化は無かった．

図‑4.9 非定常流れにおける水位〜流量曲線（300 分）- 26 -

図‑4.10 非定常流れにおける水位〜流量曲線（60 分）
60 分ハイドログラフのケースでは，300 分ケースと同様に dune の流下にしたがって波高が低くなり，上流端で発生する dune の波高自体も時間とともに徐々に低下していく過程が観測された．しかし，dune の完全な消失と平坦河床への遷移までには至らず，図‐4.6 に示されるように，非定常開始
30 分を過ぎた頃から低くなった dune 波高が再び上昇し始めて明確な dune が再形成されるようになった．この間，dune 波高は低下と上昇を示したのに対して，300 分ケースと同様，波長に大きな変化は見られなかった．図‐4.7 及び図‐4.8 に 300 分流量ハイドログラフ実験時の dune の移動状況と平坦河床になった状態を示す．dune は時間とともに波高は低くなったが，波長についてはほぼ一定値を保ちながら下流に移動している様子が確認された．なお，この写真は dune の側面部の断面形状を表している．実験で形成された dune は横断方向には直線形状を取らず，曲線で角度がついていたため，側面形状のみから dune の全体形状を計測するのは困難であった．
(4) 水位流量曲線について
非定常流れにおける河床抵抗値の多価性，河床形態の追従遅れを検証するため，水位〜流量曲線による評価を行なった．図‐4.9 及び図‐4.10 に 300 分，60 分非定常ハイドログラフにおける水位〜流量曲線を示す．図中の矢印は非定常ハイドログラフの時間に沿った動きを示す．300 分ハイドログラ
フは，流量上昇とともに水位が上昇し，ピーク流量 60l/s で dune が消失した時点で水位が低下した．
その後，流量減少の期間は水位の低減が継続した．60l/s から 58l/s にかけて水位が一時的に上昇しているが，目視では明確な河床波は確認されなかった．300 分の時間内では, 同流量に対する水位が 2 つある二価性を示している．これに対して，60 分ハイドログラフでは，dune の消失過程で水位が低下したが，流量減少の期間に dune が再形成されており，それに対応して再び水位が上昇するループを描いた．図‐4.9 に示すように，300 分ケースの最後の部分では水深が徐々に上昇しており，より長い
流量ハイドログラフを用いた場合は河床抵抗値が非定常開始時の初期値に戻るループを描く可能性もある．図‐4.9 に見られる二価性が非定常ハイドログラフの継続時間の不足，すなわち河床の追従遅れにより生じるものか，水理学的な本質解であるか，今後さらに長いハイドログラフで検証を行なう必要がある．また，300 分ハイドログラフで水位が急低下した直後に急上昇をすることについては，実- 27 -

験技術上 dune が消失した際の急激な水位低下に下流端の水位調節操作が対応し切れず，時間差を生じて水位が戻ったように見える可能性もあるが，具体的な機構は不明である．今後，この水位変動の機構解明に向けてさらに検証していきたい．

図‑4.11 流速係数の時系列変化(5) 流速係数の時系列変化による河床形状の評価
非定常流量ハイドログラフによる河床形状及び河床抵抗の変化について，流速係数による評価を試みた．流速係数は以下により定義される 12)．

(4-5)
ただし， ：流速係数，u：平均流速，u ：摩擦速度，R：径深，ks：相当粗度である．
図‐4.11 に流速係数の時系列変化を示す．図は木下 13) にならい，横軸には径深の代わりに平均水深を用いた．図中の矢印はハイドログラフの時系列を示している．300 分ハイドログラフでは，dune が発達している初期には流速係数は相当粗度 5.4mm （粒径 0.18mm の 30 倍相当）付近であり，dune の消失過程とともに流速係数が上昇し，粒径の 1 倍相当まで上昇した後，再び低下し，相当粗度 0.9mm（粒径の 5 倍相当）付近から低下が緩やかになった．60 分ハイドログラフについては，流速係数が- 28 -

上昇したものの，dune が完全には消失しなかったことから，相当粗度 0.9mm （粒径の 5 倍相当）までには到達せず，再び低下を始めた．60 分非定常ハイドログラフでは相当粗度は粒径の 5 倍以上の範囲で推移していた．
４．４

本章の結論

本研究では，dune が遷移し消失する過程とそれに対応した河床抵抗の変動過程に着目した非定常流
量波形による移動床実験を行った．径深粒径比 R/d > 500 で，側壁の影響を極力軽減するため水路幅
水深比が B/h = 5−6 の浅水条件で，かつ流量の急激な変化による攪乱が河床に与える影響を軽減する
ため，流量を微調節できる装置を用いて，非定常流れにおける dune と遷移河床間に見られる現象を再現した．本研究の結果，以下のことが明らかになった．
・ dune 波高の減少・上昇の時系列変化に対応して河床抵抗係数 Cf も同じトレンドで減少・上昇の変化をした．
・ 300 分流量ハイドログラフでは dune が完全に消滅しその後平坦河床を保ったが，60 分流量ハイド
ログラフでは dune が完全に消失せず，不完全な消失状態から再び dune が形成され始めた．
・ 300 分流量ハイドログラフでは流量全体に対して水位の二価性が見られたが，60 分流量ハイドログ
ラフでは水位が戻るループを描いた．
・ 300 分流量ハイドログラフでは流速係数は粒径の 30 倍相当から 1 倍相当まで低下し，再び上昇して 5 倍相当付近で緩やかになった．
・ 60 分流量ハイドログラフでは，流速係数は粒径の 5 倍以上までしか上昇せず，dune の再形成に伴って上昇に転じた．
４．５

参考文献

1) 芦田和男・澤井健二：河床波の変形過程に関する実験的研究, 京都大学防災研究所年報, 第 23 号
B-2, pp457-473,1980.
2) 外山明宏・清水康行・山口里実・木村一郎：非定常流場における dune の形成消滅過程に関する実験とシミュレーション，水工学論文集, 第 53 巻, pp727-732, 2009.
3) 渡邊康玄・Tubino M. and Zelezzi G.：掃流砂により形成される交互砂州の非定常流下での挙動，北海道開発土木研究所月報，No.576, pp4-12, 2001.
4) 澤井健二：水理量の時間変化に伴う河床波の変形過程に関する研究，第 31 回水理講演会論文集，647-652, 1987
5) 山口里実・泉典洋・五十嵐章：デューンの遷移過程に関する実験，水工学論文集，第 47 巻, pp613-618,
2003.
6) 音田慎一郎・細田尚：水深積分モデルによる小規模河床波の発生・発達過程と流れの抵抗則の数値解析，水工学論文集，第 48 巻, pp973-978, 2004.
7) Giri, S., Shimizu , Y.: Numerical computation of sand dune migration with free surface flow, Water
Resources Research, Vol.42, w10422, doi:10.1029/2005WR004588, 2006.
8) Giri, S., Yamaguchi, S., Shimizu, Y., Nelson, J.: Simulating temporal response of bedform characteristics tovarying flows, River, Coastal and Estuarine Morphodynamics., RCEM2007, pp939-947, 2007.
- 29 -

9) 岸力・黒木幹男：移動床流における河床形状と流体抵抗(I)，北海道大学工学部研究報告，第 67 号,
pp1-23, 1973.
10) Garde, R. J., Ranga Raju, K. G.: Regime Criteria for Alluvial Streams, Journal of the Hydraulics Division,
ASCE,No.HY-6, Nov., 1963.
11) 関根正人：移動床流れの水理学，共立出版，2005.
12) 吉川秀夫：流砂の水理学，丸善株式会社，1985.
13) 木下良作：石狩川河道変遷調査，科学技術庁資源局資料，36 号，1961.

- 30 -