掃流砂の非平衡性を考慮した河床デューンの線形安定解析
Linear Stability Analysis of Fluvial Dunes
with the Use of Non-Equilibrium Bed Transport Formula
泉 典洋 1・佐藤博重 2
Norihiro IZUMI and Hiroshige SATO
12北海道大学大学院 工学研究科 教授
北海道大学大学院 工学研究科 博士課程前期

要

旨

洪水時，lower regime で形成された河床デューンは，流量がある値を超えると消滅し，平坦床が現れたりあるいはアンチデューンが現れたりする．そして洪水後，流量が減少するにつれて，河床は再びデューンに覆われることになる．
デューンの形成は河道の抵抗を大きく増加させることから，流れに対する河道の抵抗を精度良く予測するためには，河床の安定性とそれによるデューン，アンチデューンの形成条件を明らかにすることが重要となる．Colombini1)は，開水路乱流のモデルとして混合距離モデルを用いた線形安定解析を行い，その解析結果が実験結果を良好に再現できることを見出している．
本論文は，Colombini1)が用いた混合距離モデルに加えて，非平衡性を加味した掃流砂量式を用いた線形安定解析を行い，流砂の非平衡性が河床の安定性とデューンおよびアンチデューンの形成条件に及ぼす影響を調べたものである．
解析によれば，lower regime では，平坦床はわずかに安定化し，デューン形成領域が縮小する傾向にあるのに対して，upper regime では，平坦床は不安定化し，アンチデューン形成領域が拡張する傾向にあることが明らかとなった．
《キーワード：線形安定解析；河床デューン；アンチデューン；非平衡性；掃流砂》

‑ 1 ‑

y, V
R(x)+D(x)
Flow
B(x)
R(x)

x, U

図–1 The coordinate system.

1.

はじめに比較的小さい流量の下では，河床はデューンと呼ばれる小規模河床波に覆われる．洪水などのために流量が増加するとこのデューンは消滅し，平坦床あるいはアンチデューンと呼ばれる別の河床波が現れる．洪水後流量が減少するにつれて，河床は再びデューンに覆われることになる．
河床波の形成については，河道抵抗推定の重要性から半世紀以上にわたって数多くの研究が行われてきた．適当な流れの方程式に土砂の連続式を組み合わせ，様々なスケールの河床形態について線形安定解析の概念を上手く適用することで，その形成プロセスの説明に成功した研究者も多い．デューンやアンチデューン，リップルの解析に，流体の安定性に関する研究の中で発展してきた手法が用いられてきた2),5),4),6),3),7) .
Colombini1) は，Engelund2) によって提案された定剪断層モデルよりも開水路乱流をより精密に表現できると言われている混合距離モデルを用いたデューンの線形安定解析を提案している．この論文では，Colombini のモデルに，掃流砂が有する非平衡性の影響を加味したモデルを用いて，線形安定解析を行ったものである．

2.

定式化

(1) 無次元支配方程式
Reynolds 平均を取った無次元の Navier-Stokes 方程式は次のように書ける．
UU

∂U
∂P
∂T xx ∂T xy∂U+V=−+1++∂x∂y∂x∂x∂y∂T xy ∂T yy∂V∂V∂P+V=−+ S −1 ++∂x∂y∂y∂x∂y∂U ∂V+=0∂x
∂y(1)
(2)
(3)

ここで x および y はそれぞれ流れ方向および横断方向の座標，U および V はそれぞれ流速の流れ方向および横断方向成分，S は平均河床勾配，P は圧力，T i j (i, j = x, y) は二次元の Reynolds 応力テンソルである．上式では，流れの時間変化が河床形状の時間変化に比べて十分早いという準定常の仮定が用いられており，非定常項は既に落としてある．また上式の変数は既に平坦床等流状態における摩擦速度 U ∗f 0 ，
水深 D∗0 ，水の密度 ρ を用いて次のように無次元化されている．

(U ∗ , V ∗ ) = U ∗f 0 (U, V)(4a)(x∗ , y∗ ) = D∗0 (x, y)P∗ , T i∗j = ρU ∗2f 0 P, T i j(4b)(4c)

混合距離モデルを用いると，Reynolds 応力テンソルは次のように表すことができる．

T xx = 2νT2∂U
∂x(5a)

T xy = νT

∂U ∂V+∂y∂x!
(5b)∂V(5c)∂y∂U(5d)
νT = l2∂y D − y 1/2
l = κ (y − Z)(5e)D
ここで νT は U ∗f 0 および D∗0 で無次元化した渦動粘性係数，l および Z ，D はそれぞれ D∗0 で無次元化した混合距離および河床高さ，水深，κ はカルマン定数 (= 0.4) である．
ここで次のように定義される流れ関数を導入する．
!
∂ψ ∂ψ,−(6)(U, V) =∂y ∂x
T yy = 2νTすると式 (1) は次のように書き表される．

!
∂ψ ∂2 ψ∂ψ ∂2 ψ∂∂2 ψ∂P−=−+1+2νT∂y ∂x∂y ∂x ∂y2∂x
∂x
∂x∂y"
!#
2∂∂ ψ ∂2 ψ+νT−∂y∂y2 ∂x2
ψ を用いて書き直した式 (2) と式 (7) から P を消去すると次式が得られる．
!
∂ψ ∂∇2 ψ ∂ψ ∂∇2 ψ∂2∂2 ψ
νT−−4∂y ∂x∂x ∂y∂x∂y∂x∂y!"!#222∂∂ ψ ∂2 ψ∂νT=0+−−∂x2 ∂y2∂y2 ∂x2
水面および底面における境界条件の適用を容易にするために次のような変数変換を導入する．
!
y − R (x)(ξ, η) = x,D (x)

(7)

(8)

(9)

ここで R は対数分布則において流速がゼロとなる基準高さである．上の変換を用いると，無次元の混合距離は次のように書くことができる．
1/2
R − Z
Rl = κD η +
1− −η
(10)
DD

(2) 境界条件
水面および底面では次の境界条件が満足されなければならない．
u · ens = 0 at
ens · T · ens = 0

η=1(11a)atη=1(11b)

ets · T · ens = 0 at

η=1(11c)u · enb = 0 at

η=0(11d)u · etb = 0 at

η=0(11e)

ここで u (= (u, v)) は流速ベクトル，ets および ens はそれぞれ水面に対する単位接線ベクトルおよび単位法線ベクトル，etb および enb はそれぞれ底面に対する単位接線ベクトルおよび単位法線ベクトル，T は次式で表される応力テンソルである．
"#−P + T xx
T xy
T=(12)T xy−P + T yy

3yFlowBRxττcτb τr

図–2 Definition of the bedload layer.

3.

基準状態
線形安定解析の基準となる状態は平坦床等流状態である．そのときの変数を次のように表す．

(U, V, D, Z, R) = (U0 , 0, 1, 0, R0 )

(13)

すると支配方程式は次のように書き直される．

1+dT xy0=0dη

(14a)

dU0dη

(14b)

dU0dη

(14c)

T xy0 = νT 0
νT 0 = l02l0 = κ (η + R0 ) (1 − R0 − η)1/2(14d)ここで添え字 0 は基準状態を表す．式 (14a)–(14d) を適当な境界条件の下で解くと，次の対数流速分布が得られる．
!
1
η + R0U0 = ln(15)κR0
ここで R0 は基準状態での基準高さを表しており，河床材料の粒径と次のような関係がある．

R0 =ds6(16)

ここで d s は無次元化された砂の粒径である (= d∗s /D∗0 )．

4.

非平衡掃流砂量
土砂の運動は次式で表される．

ρs

∗!
∂q∗s
∗ ∂q s+ u p ∗ = τb − τq∂t∗∂x(17)

ここで ρ s は砂粒子の密度，u∗p は砂粒子の速度，τb および τq は掃流層上面および下面における無次元の剪断応力である（図–3 参照）．掃流層上面の剪断応力 τb が流れによって決定されるのに対して，掃流層下面の剪断応力 τq は砂の動きによって決まることに注意する必要がある．
砂粒子の速度 u∗p は掃流層内の平均速度で近似し，次のように表す．

u∗p =q∗s(1 − λ p )h∗b4(18)図–3 Definition of shear stresses at the upper and lower surfaces of the bedload layer.

ここで h∗b は掃流層厚さ，λ p は空隙率である．空隙率は一般に 0.4 程度であることが知られている．ところが，活発な運動をしている掃流層内では λ p はそれよりかなり大きくなることが予想される．
砂粒子の加速減速の時間スケールは河床形状変化の時間スケールに比べて十分早いことから，式 (17)の時間微分項についても無視する．
無次元化した式 (17) は次のように書ける．

βΦ ∂Φ= θb − θqhb ∂x(19)1/2 ), R は水中比重 (= 1.65), θ および θ はそれぞれ掃流層
ここで Φ は無次元の掃流砂量 (= q∗s /(R s gd∗3sbqs )上面および下面の Shields 応力 ((θb , θq ) = (τb , τq )/ρR s gd∗s ) である．
無次元パラメータ β は次のように書ける．

!
d∗s 2 (R s + 1)d2sρs=β=ρ(1 − λ p ) D∗0
1 − λp(20)

前述したように，掃流層下面の Shields 応力 θq は砂粒子の運動によって決定される．もし平衡状態にあれば θq は θb と釣り合っており，砂粒子は加速も減速もしない．したがって平衡状態では，通常の掃流砂量式で θb の代わりに θq を代入した式で表される．平衡掃流砂量式として，ここでは次の Meyer-Peter& Müller 公式を用いる．
3/2Φ = 8 θq − θc(21a)∂Bθc = θch − µ S −∂xµ=!

θchtan φ(21b)(21c)

ここで θc は限界 Shields 応力，θch は平坦床に対応する θc ， B は掃流層上面の高さ，φ は摩擦角である．
これらの式を，より一般的な非平衡状態に拡張する．
Colombini1) によれば，掃流層厚さ hb は次のように表される．

h b = lb d s
τr − τcτclb = 1 + 1.3(22a)!0.55(22b)

ここで τr および τc はそれぞれ基準高さにおける剪断応力および限界剪断応力である．
河床形状の時間変化は次の Exner 方程式で表される．

∂B ∂Φ
+=0∂t
∂x5(23)5.

線形安定解析

(1) 流れの方程式の線形化
全ての変数を次のように展開する．
(ψ, P, D, Z, R, B) = (ψ0 , P0 , 1, 0, R0 , B0 )
+ A ψ̂1 , P̂1 , D̂1 , R̂1 , R̂1 , R̂1(24)

ここで A は擾乱の振幅であり，十分小さいと仮定する．上式を支配方程式 (7) および (8) に代入する．
O(A) では次式が得られる．
L̂ψ ψ̂1 + L̂D D̂1 + L̂R R̂1 = 0(25)∂P̂1+ Ŵψ ψ̂1 + Ŵ D D̂1 + ŴR R̂1 = 0∂ξ(26)

ここで L̂φ および Ŵφ は線形演算子であるが，その具体的な形はスペースの都合上省略する．境界条件(11a)–(11e) から次式が得られる．
∂ψ̂1P̂1 == 0 at η = 1(27)∂ξψ1 =∂ψ̂1= 0 at∂η

η=0擾乱は次のように表されると仮定する．
ψ̂1 , P̂1 , D̂1 , R̂1 = (ψ1 , P1 , D1 , R1 ) exp[i(αξ − Ωt)]

(28)

(29)

ここで α および Ω はそれぞれ擾乱の波数および複素増幅率である．上式を式 (25)–(28) に代入すると次式が得られる．
Lψ (η)ψ1 (η) + LD (η)D1 + LR (η)R1 = 0(30)Wψ (1)ψ1 (1) + W D (1)D1 + WR (1)R1 = 0(31)ψ1 (1) = 0(32)ψ1 (0) = 0(33)Dψ1 (0) = 0(34)

ここで D = ∂/∂η である．
摂動された流関数 ψ1 を Chebyshev 多項式を用いて次のように展開する．

ψ1 =NXan T n (ζ)

(35)

n=0

ここで T n は n 次の Chebyshev 多項式，ζ は領域 [-1, 1] で定義される Chebyshev 多項式の変数である．計算の精度を上げるため，次の変数変換を用いた．
( )
ln (η + R0 ) /R0
ζ=2−1(36)ln [(1 + R0 ) /R0 ]これらの式を支配方程式に代入し，全ての式を次の Gauss-Labatte 点で評価する．

ζi = cos ( jπ/N) ,( j = 1, 2, · · · , N − 2)(37)

すると次の形の線形系が得られる．

La = MR1(38a)a = [a0 , a1 , · · · , aN , D1 ](38b)

ここで60 
 ĎT 0 (−1) · · · ĎT N (−1)
 T (−1)···
T N (−1)0 0 Ľψ T (1)···
Ľψ T N (1)ĽD 0
···
···L = 
···
 Ľψ T (ζ ) · · · Ľψ T (ζ ) ĽD 0N−2
NN−2
 W̌ψ T (1) · · · W̌ψ T (1) W̌ D 0N
T 0 (1)···
T N (1)0hi
M = 0, 0, ĽR , · · · , ĽR , W̌R , 0(38c)(38d)

また ˇ は η を ζ に変数変換した線形演算子を表している．
式 (38) は容易に解くことができて，Chebyshev 多項式の係数 an (n = 0, · · · , N) および D1 を求めることができる．すると摂動された流関数 ψ1 と水深 D1 が R1 を因数として含む形で得られる．

(2) 非平衡層流砂量式の線形化
無次元掃流砂量 Φ および Shields 応力は次のように展開される．
Φ, θb , θq = (Φ0 , θb0 , θb0 )+ A Φ1 , θb1 , θq1 exp i(αξ − Ωt)(39)

上式を用いれば式 (19) と (21) は次のように書き直される．

Φ0 ∂Φ1= θb1 − θq1
hb ∂x(40)Φ1 = 12 (θb0 − θc )1/2 θq1(41)!
12 (θb0 − θc )1/2 hb
Φ1 =θb1
hb + 12iαβ (θb0 − θc )1/2 Φ0(42)β

上式から次式が得られる．

摂動された Shields 応力 θb1 は前の節で得られた結果を用いて見積もることができる．上式を線形化された Exner 方程式 (23) に代入すると，次の一般形で表される Ω が得られる．

Ω = f (α, F; C, µ, m)(43)Ω の虚部が擾乱の成長率を表している．

6. 結果と考察
図–4 は α–F 平面における擾乱の成長率 Im[Ω] の中立曲線を示したものである．太い実線が Colombini1)の解析結果であり，他の線が本解析の結果である．
Froude 数が大きい領域と小さい領域に二つの不安定境域が現れるが，前者がアンチデューンに，後者がデューンに相当している．空隙率 λ p が小さいと掃流砂の非平衡性の影響はほとんどないが，λ p が大きくなるにつれて無視できないくらい大きな影響を持つことがわかる．空隙率によって非平衡性の影響が決まるのは式 (20) および (42) からも明らかである．
アンチデューンの領域は非平衡性の影響によって空隙率が大きくなると小波数領域と Froude 数の小さい領域に広がることがわかる．一方，デューンの領域は比較的波数が大きく，Froude 数が大きい領域において小さくなることがわかる．さらにそれ以上ではデューンが消滅する限界 Froude はわずかに減少することがわかる．非平衡性の影響は平坦床基準状態を，デューン領域で安定化させ，アンチデューン領域で不安定化させる傾向にあることがわかる．

7Colombini’sanalysis (2004)

λp =0.4

1.5

0.70.9F
1.00.500.5

ǩ

1.01.5

図–4 Instability diagrams obtained by the present model (C −1 = 20, µ = 0.1).

7.

結論

混合距離モデルと非平衡掃流砂量式を用いて，河床デューンの線形安定解析を行った．Colombini1) の解析結果と比較すると，非平衡性の影響は掃流層内の空隙率が大きくなるとき無視し得ないほど大きなものとなる．また非平衡性の影響によってアンチデューン領域は大きくなり，デューンの領域は小さくなる．また臨界 Froude 数は若干小さくなることがわかった．
掃流層内の空隙率の正確な値については不明であるため，ここではパラメータとして与えている．非平衡性の影響をより正確に評価するためには空隙率 λ p のより詳しい研究が必要である．今後の課題となろう．
参考文献
1) Colombini, M. (2004). Revisiting the linear theory of sand dunes formation, J. Fluid Mech. 502, 1–16.
2) Engelund, F. (1970). Instability of erodible beds, J. Fluid Mech. 42, 225–244.
3) Engelund, F. & Fredsoe J. (1982). Sediment ripples and dunes. Annu. Rev. Fluid Mech. 14, 13–37.
4) Fredsoe, J. (1974). On the development of dunes in erodible channels. J. Fluid Mech. 64, 1–16.
5) Hayashi, T. (1970). Formation of dunes and antidunes in open channels. J. Hydraul. Div. ASCE 96, 357–366.
6) Richards, K. J. (1980). The formation of ripples and dunes on an erodible bed. J. Fluid Mech. 99, 597–618.
7) Sumer, B. M. & Bakioglu, M. (1984). On the formation of ripples on an erodible bed. J. Fluid Mech. 144, 177–190.

8