混合距離モデルを用いた河床デューンの弱非線形安定解析
Weakly Nonlinear Stability Analysis of Fluvial Dunesby the Use of the Mixing Length Turbulent Model
泉 典洋Norihiro IZUMI北海道大学 大学院工学研究科 教授

要

旨

多重尺度摂動展開を用いた増幅率展開法を使って，河床デューンの弱非線形安定解析を行った．これまで開水路乱流については，Engelund1)や Fredsøe2)の用いていた定剪断層近似を用いていたが，本解析では，現実の開水路乱流をより良好に表現できると言われている混合距離モデルを用いた．
支配方程式を摂動展開することによって得られる一連の摂動方程式を解くに当たっては，Chebyshev 多項式展開を適用したスペクトル法を用いた．
解析の結果によれば，抵抗係数（C-1）が大きい領域において，デューン−平坦床遷移過程は亜臨界分岐によって特徴付けられることが明らかとなった．
これによって，デューン−平坦床遷移時に見られるヒステリシス現象の一部が亜臨界分岐によって生じていることが示唆された．
《キーワード：デューン；線形安定解析；弱非線形安定解析；亜臨界分岐；混合距離モデル》

‑ 1 ‑

1.

はじめにフルード数がある程度小さい時，河床には水深の数割程度の波高を持ったデューンと呼ばれる河床波が形成される．デューンで覆われた河床（デューン河床）はフルード数の増加とともに平坦床へと遷移し，平坦床はフルード数の減少とともにデューン河床へと遷移する．ところがこの遷移時に，デューン河床から平坦床へと遷移する際の臨界フルード数の方が，平坦床からデューン河床へ遷移する際の臨界フルード数より大きい，いわゆるヒステリシス現象が現れることが知られている．
著者らはこれまで Engelund1) および Fredsøe2) の提案した定剪断層近似 (constant stress layer approximation) を用いてデューンの弱非線形安定解析を行い，平坦床–デューン遷移時の分岐形態が亜臨界分岐であることが，ヒステリシスの現れる原因の一つである可能性を示している3),4) ．本研究は，定剪断層近似より実際の乱流をよりよく表現できる混合距離モデルを用いた Colombini5) の線形安定解析をベースにして，河床デューンの弱非線形安定解析を行ったものである．

2.

定式化
開水路内の乱流は，レイノルズ平均を取った次の Navier-Stokes 方程式によって表わすことができる．

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

ここで x および y はそれぞれ流下方向および水深方向の座標，U および V はそれぞれ x および y 方向の流速，S および P はそれぞれ平均河床勾配および圧力，T i j (i, j = x, y) はレイノルズ応力テンソルである．上式中では流れの時間変化は河床形状の時間変化に比較して十分速いと仮定し，非定常項を無視している．非定常の効果は，後述する河床形状の時間変化を表わす式でのみ考慮する．また上式ではすでに次のような無次元化が行われている．

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

ここで U ∗f 0 および D∗0 はそれぞれ平坦床基準状態における底面摩擦速度および水深である．混合距離モデルを用いるとレイノルズ応力テンソルは次のように表すことができる．

∂U,T xx = 2νT∂x∂VT yy = 2νT,∂y∂U,
νT = l∂y2T xy = νTl = κ (y − Z)∂U ∂V+∂y∂x

 D − y 1/2!
(5a–c)(5d, e)D

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

∂P∂ψ ∂2 ψ∂ψ ∂2 ψ∂∂ψ=−−+1+2νT2∂y ∂x∂y ∂x ∂y∂x∂x∂x∂y"!#22∂ ψ ∂ ψ∂
νT−+∂y∂y2 ∂x22!

(7)

y, VR(x)+D(x)FlowB(x)R(x)x, U

図–1 流れと座標系の概念図．

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

ξ = x,

η=y − R(x)D(x)(9a, b)

ここで R は対数速度分布において流速がゼロとなる高さ（以降，基準高さと呼ぶ．図–1 参照），(),x はx に関する微分を表わしている．上式の変数変換を用いると無次元混合距離 l は次のように表わされる．

1/2
RR − Zl = κD η +1− −ηDD(10)

水面および底面における境界条件は次のようになる．

u · ens = 0atη=1(11)ens · T · ens = 0atη=1(12)

ets · T · ens = 0atη=1(13)u · enb = 0atη=0(14)u · etb = 0atη=0(15)

ここで u は流速ベクトル (= (u, v))，ens および ets は水面に対するそれぞれ法線および接線方向の単位ベクトル，enb および etb は底面に対するそれぞれ法線および接線方向の単位ベクトルである．上式中の水面位置の座標では R/D を微小として無視している．また T は応力テンソルであり次式で表される．
"#T xy−P + T xx(16)T=T xy−P + T yy3.

基本解
安定解析の基本状態は平坦床等流状態である．基本状態では各変数は次のように表わすことができる．

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

(17)

このとき支配方程式は次のように単純化される．

1+
νT 0 = l02

dT xy0= 0,dη

dU0,dηT xy0 = νT 0dU0dη

l0 = κ (η + R0 ) (1 − R0 − η)1/2

3(18a, b)(18c, d)

ds/6
mds/30y = R0y=0ds

図–2 砂の粒径と原点 (y = 0) および基準高さ (y = R0 ) の関係．

ここで ()0 は基本状態における解を表わしている．上式を解くにあたっては次の境界条件を用いる．

U = 0,T xy0 = 1 − R0at η = 0(19, 20)式 (18)–(20) より次の対数分布則が得られる．

η + R01U0 = lnκR0!

上式を η = 0 から 1 まで積分すると次の抵抗則（抵抗係数 C ）が得られる．
"!
#
∗Ua0
1
1 + R0−1(1 + R0 ) ln
C = ∗ =−1Uf0 κ
R0

(21)

(22)

∗ は基本状態における水深平均流速である．
上式中の Ua0
ここで R0 がどの程度の物理的大きさを持っているか調べておこう．対数分布則を変数変換以前の変数を用いて書き表すと次のようになる．
!
!
y1
30y
1+ 8.5 = ln
(23)
U0 = lnκks
κ
md s

ここで k s は D∗0 で無次元化した等価粗度高さ，m は k s /d s ，d s は D∗0 で無次元化した砂の粒径である．
y = md s /30 の時，U0 = 0 であるから R0 = md s /30 となる．等価粗度高さは砂の粒径の 1∼3 倍 (m = 1 ∼ 3)と言われていることから，R0 は d s の 1/30 から 1/10 程度の値の範囲にあることがわかる．通常の混合距離モデルでは，原点近傍に存する特異性のために流速が急激に減少し −∞ に発散してしまうため，実際の流速分布を適切に表すことはできない．ところが，上記の結果によれば，砂面上における基準高さ
R0 は砂の粒径に比して十分小さくなるため，原点を砂粒子の最上点より d s /6 だけ下（体積平均を取った場合の河床面）に取れば，基準高さ y = R0 は砂粒子の最上点 y = d s /6 よりもさらに下に位置することになる（図–2 参照）．したがって混合距離モデルが持つ原点近傍の特異性は無視できて，砂面上の流れは混合距離モデルによって十分表現できることが期待できる．

4. 掃流層モデル
混合距離モデルを用いたデューンの線形安定解析では底面の極近傍で河床形状と流速の間の位相差が大きく変化するため，流砂量を見積もる剪断力の位置によって安定性は大きく変化する．このことに気付いた Colombini5) は，独自の掃流層モデルを導入することによって，混合距離モデルを用いたデューンの線形安定解析結果を大幅に改善した．モデルの概要は以下の通りである．
河床上における砂の連続式は次のように表される．

∂B ∂Φ
+=0
∂t
∂x(24)

1/2 )，q∗ は掃流砂量，R
ここで B は掃流層上面の高さ（図–1,3 参照），Φ は無次元掃流砂量 (= q∗s /(R s gd∗3ss )s∗ )，λ は空隙率である．Colombini5)(1−λ)]tは砂粒子の水中比重 (= 1.65)，t は時間 (= [d∗s (R s gd∗s )1/2 ]/[D∗2pp04

y(a)

y(b)

FlowFlowB
RRx, τx, ττcτr

τb τr

図–3 掃流層の定義5) ．

は，掃流砂量の算出に Meyer-Peter & Müller 公式を採用し，掃流層上面における剪断力を用いて次のように表されるとした．
Φ = 8 (θb − θc )3/2
(25)
ここで流砂の非平衡性は無視されている．θb および θc はそれぞれ掃流層上面におけるシールズ剪断力および限界シールズ剪断力であり，次のように表される．

τ∗bDS=
τb
ρR s gd∗s
R s d∗s!
∂B
θchθc = θch − µ S −, µ=∂xtan Ψθb =(26)(27)

ここで θch は平坦床に対応する θc ，Ψ は摩擦角である．さらに Colombini5) は従来の実験結果（図–4）から掃流層厚さ hb を次のように表わした．
!0.55
τr − τchb = lb d s , lb = 1 + 1.3(28a, b)τc
ここで τr および τc はそれぞれ基準高さ (η = R) における無次元の剪断力および限界剪断力である．掃流層厚さ hb と掃流層上面の高さ B の間には次の関係がある．

B = hb + d s /6 = (lb + 1/6) d s(29)

果たして Colombini5) の掃流層モデルが物理的な見地から本当に妥当であるかどうかについては議論の余地があるが，後述するように線形安定解析の結果と実験結果との一致が大幅に改善されているのは事実である．ここでは最初のステップとして Colombini5) のモデルをそのまま用いて弱非線形解析へと拡張することとする．

5.

線形安定解析
先に求めた基本解に対し摂動を与える．そのとき各変数は次のように展開することができる．

(ψ, P, D, Z, R, B) = (ψ0 , P0 , 1, 0, R0 , B0 )
+ A ψ̂1 , P̂1 , D̂1 , R̂1 , R̂1 , R̂1(30)

ここで A は摂動の振幅を表すパラメータであり，線形安定解析のスキームでは無限小と考える．上式を支配方程式 (7) および (8) に代入して A のオーダーで整理すると，O(A) において次式が得られる．

L̂ψ ψ̂1 + L̂D D̂1 + L̂R R̂1 = 0(31)∂P̂1+ R̂ψ ψ̂1 + R̂R R̂1 + R̂D D̂1 = 0∂ξ

(32)

5

lb

τr/τc

図–4 サルテーション高さに関する従来の実験結果5) と提案式の比較．実線は式 (28b) を表している．

また境界条件 (11)–(15) （式 (13) は常に成立するため不要）および Exner 方程式 (24) からは次式が得られる．
∂ψ̂1 /∂ξ = 0, P̂1 = 0 at η = 1(33, 34)∂ψ̂1 /∂ξ = 0,∂ψ̂1 /∂η = 0 at

η=0Êψ ψ̂(ηB ) + ÊR R̂1 + ÊD D̂1 = 0(35, 36)(37)

ここで ηB = (B − R)/D である．また L̂φ および R̂φ ，Êφ (φ = ψ, D, R) は線形演算子であり，スペースの関係から具体的な形については省略する．
摂動がフーリエ級数で表されるとして，ノーマルモード解析を行う．単一波数に注目し摂動を次のように表す．
ψ̂1 , P̂1 , D̂1 , R̂1 = (ψ1 , P1 , D1 , R1 ) exp i (αξ − Ωt)(38)
ここで α および Ω はそれぞれ摂動の波数および複素周波数を表している．上式を式 (31)–(37) に代入すると次式が得られる．

Lψ (η) ψ1 (η) + LD (η) D1 + LR (η) R1 = 0ψ(39)

DR(40)

ψ1 (1) = 0,P1 (1) = 0(41, 42)ψ1 (0) = 0,Dψ1 (0) = 0(43, 44)iαP1 (η) + R (η) ψ1 (η) + R (η) D1 + R (η) R1 = 0

ここで D = d/dη である．また式 (37) と，式 (40) および (42) からそれぞれ次式が得られる．

Eψ ψ1 (ηB ) + ER R1 + ED D1 = 0
ψ

R(45)

D(46)

an T n (ζ)(47)

R (1) ψ1 (1) + R (1) R1 + R (1) D1 = 0
ψ1 を Chebyshev 多項式展開を用いて次のように表す．
ψ1 =NXn=0

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

ζ j = cos( jπ/N),( j = 1, · · · , N − 2)6(49)-100-1-101.5
2-10F

1-100-1

10.8-10.6

0.5

0.40.200.5

1.01.5α
図–5 増幅率 Im[Ω] のコンタ．C −1 = 20, µ = 0.1, m = 2.5

それを境界条件と合わせると次の線形代数方程式系が得られる．
La = 0(50)a = [a0 , a1 , · · · , aN , D1 , R1 ] Ěψ T 0 (ζB ) · · · Ěψ T N (ζB ) ĚD ĚR 


···
T N (−1)00 
 T 0 (−1)
 ĎT (−1) · · · ĎT (−1)
00 0N
T 0 (1)···
T N (1)00 

···Řψ T N (1)ŘD ŘR 
L =  Řψ T 0 (1) Ľψ T 0 (ζ1 ) · · · Ľψ T N (ζ1 ) ĽD ĽR 
······· 
······· 
Ľψ T 0 (ζN−2 ) · · · Ľψ T N (ζN−2 ) ĽD ĽR(51)

ここで(52)式 (52) 中の ˇ は η を ζ に変数変換した線形演算子を表している．式 (50) が自明でない解を持つためには可解条件 |L| = 0 が満足される必要がある．その条件から複素周波数 Ω は次のような関数形で求められる．
Ω = f (α, F; C, µ, m)(53)
ここで求められた Ω の虚部が摂動の増幅率に相当する．
図–5 に摂動の増幅率 Im[Ω] の α–F 平面上におけるコンタを示した．フルード数 F が小さい領域と大きい領域に摂動の増幅率 Im[Ω] が正となる領域が現れているが，前者がデューン，後者がアンチデューンの領域である．デューンの領域についてはこれまでの解析結果と大きな差はないが，アンチデューンの領域は大きく変化し，波数が 0.5–1.0 の範囲にのみ存在している．Colombini5) は Runge-Kutta 積分によるシューティング法を用いて同様の線形安定解析を行い，この図とほぼ同じ図を得ている．Colombini5)によれば，デューンおよびアンチデューンに関する従来の実験結果を図中にプロットすると，いずれもIm[Ω] > 0 の領域の Im[Ω] がピークを取る付近に現れ理論と実験結果の一致は大幅に改善されるという．

6. 弱非線形安定解析(1) 摂動展開
臨界フルード数の極近傍（わずかに小さい領域）を考え，次式で定義される微小パラメータ  を導入する．
!1/2
Fc − F
=(54)
Fc7ここで Fc は臨界フルード数である．また多重尺度法を用いて，速い時間 T 0 と遅い時間 T 1 を導入する．
そのとき次の関係が成り立つ．

T 0 = t,∂∂∂=+ 2∂t ∂T 0
∂T 1

T 1 =  2 t,

(55a, b, c)

これらに合わせて各変数を次のように展開する．

(ψ, D, R) = (ψ0 , 1, R0 ) +3X

 i (ψi , Di , Ri )

(56)

i=1

上式を式 (7) および (8)，(11)–(15) に代入して  のオーダーで整理する．基本擾乱として O() における解を次のような形で与える．
(ψ1 , D1 , R1 ) = (ψ11 , D11 , R11 ) E(57)

ここで E = exp i (αξ − Ωt) である．すると O( 2 ) および O( 3 ) における解はそれぞれ次のような形を有していることが予想される．

(ψ2 , D2 , R2 ) = A2 (ψ22 , D22 , R22 ) E 2 + c.c.
+ AĀ (ψ20 , D20 , R20 ) + (ψ00 , 0, 0)
3(58)

3

(ψ3 , D3 , R3 ) = A (ψ33 , D33 , R33 ) E + c.c.
+ (ψ31 , D31 , R31 ) E + c.c.

(59)

式 (56)–(59) を式 (6)，(7)，(11)–(15) に代入して整理すれば  の各オーダーにおいて次のような結果が得られる．
O():ψL1 ψ11 + L1D D11 + LR1 R11 = 0(60)

ψE1 ψ11 (ηB ) + E1D D11 + ER1 R11 = 0(61)ψ1 (0) = 0(62)D1 ψ11 (0) = 0ψDR1 (1)ψ11 (1) + R1 (1)D1 + R1R (1)R1 = 0(63)ψ1 (1) = 0(65)L2 ψ22 + L2D D22 + LR2 R22 = N22(66)ψE2 ψ22 (ηB ) + E2D D22 + ER2 R22 = S22(67)(64)

O( 2 ):ψ

ψ1 (0) = 0D1 ψ11 (0) = 0
ψ
DR1 (1)ψ11 (1) + R1 (1)D1 + R1R (1)R1 = P22
ψ1 (1) = 0
ψ

(68)
(69)
(70)
(71)

L0 ψ20 + L0D D20 + LR0 R20 = N20(72)ψE0 ψ20 (ηB ) + E0D D20 + ER0 R20 = S20(73)ψ20 (0) = 0Dψ20 (0) = 0ψDR0 (1)ψ20 (1) + R0 (1)D20 + R0R (1)R20 = P20
ψ1 (1) = 08(74)
(75)
(76)
(77)

表–1 臨界フルード数と Landau 定数．

C −1µm

αFc

Ω

λ0

λ1

212121

0.1
0.1
0.1

1.7
1.7
1.7

0.217
0.227
0.237

0.82540.82560.8254

12.6
11.9
13.2

11.2 + 51.2i
11.6 + 54.2i
11.9 + 57.1i−2040 + 57500i−3770 + 58600i−5290 + 59100i

2222220.1
0.1
0.1

1.71.7
1.7

0.238
0.248
0.258

0.84910.84920.8491

26.7
28.1
29.4

30.4 + 139i
31.5 + 147i
32.5 + 155i

5210 + 173000i
1890 + 176000i−999 + 177000i

21
21
21

0.05
0.05
0.05

1.7
1.7
1.7

0.219
0.229
0.239

0.82680.82700.8269

12.2
12.9
13.5

11.7 + 53.0i
12.1 + 56.1i
12.5 + 59.0i−1540 + 59900i−3300 + 61000i−4830 + 61700i

2121
21

0.1
0.1
0.1

2.5
2.5
2.5

0.258
0.268
0.278

0.8496
0.8496
0.849533.6
35.2
36.841.3 + 177i
42.8 + 187i
44.2 + 196i

11500 + 203000i
8240 + 206000i
5360 + 208000i

O( 3 ):
ψ

L1 ψ31 + L1D D31 + LR1 R31 = N31(78)ψE1 ψ31 (ηB ) + E1D D31 + ER1 R31 = S31(79)ψ31 (0) = 0
Dψ31 (0) = 0
ψR1 (1)ψ31 (1) + R1D (1)D31 + R1R (1)R31 = P31
ψ1 (1) = 0
φ

φ(80)
(81)
(82)
(83)

φ

ここで Ln および Rn ，En (φ = ψ, D, R; n = 0, 1, 2) は，それぞれ Lφ および Rφ ，Eφ において Im[Ω] = 0 お
よび (F, k) = (Fc , nk) とした線形演算子，Ni j および Si j ，Pi j ((i, j) = (1, 1), (2, 2), (2, 0), (3, 1)) はより低次の解からなる非同次項を表わしている．

(2) 数値解法
線形解析と同様に Chebyshev 多項式展開を用いたスペクトル法を用いて解く．ψi j を次のように展開する．
NX(i j)
ψi j =
an T n (ζ)
(84)
n=0

上式を式 (60)–(83) に代入すると，一連の線形代数方程式系が得られる．
O() では次式が得られる．
L1 a11 = 0
(i j)

(i j)

(85)

(i j)

ここで ai j = (a0 , a1 , · · · , aN , Di j , Ri j )T ，Ln は L において Im[Ω] = 0 および (F, α) = (Fc , nα) とした行列を表わしている．線形安定解析の結果から判るように |L1 | = 0 であり，上式は自明でない解を持つ．
O( 2 ) では次式が得られる．

ここでL2 a22 = N22(86)L0 a20 = N20 ,(87)Ni j = Ši j , 0, 0, 0, P̌i j , Ňi j (ζ1 ), · · · , Ňi j (ζN−2 ) T(88)9A
λ1 < 0 ߩ႐ว⿥⥃⇇ಽጘλ1 > 0 ߩ႐ว੝⥃⇇ಽጘFcF

λ0 > 0 ߩ㗔ၞ

λ0 < 0 ߩ㗔ၞ

図–6 超臨界分岐と亜臨界分岐．

O( 3 ) では次式が得られる．

N31 = 

L1 a31 = N31(2)A2 A∗ Š(1)31 + AŠ31 − R11000(89)dA
dT 1(2)A2 A∗ P̌(1)31 + AP̌31(1)(2)(ζ1 ) + AŇ31(ζ1 )A2 A∗ Ň31··(1)(2)
2∗A A Ň31 (ζN−2 ) + AŇ31(ζN−2 )













(90)

ここで |L1 | = 0 であるから，式 (89) が解を持つためには次の可解条件が満たされなければならない．

dA
dT 1

Ě1 T 0 (ζB )ψ···Ě1 T N (ζB )ψĚ1DT 0 (−1)···

T N (−1)

00Ď1 T 0 (−1)ψ···

Ď1 T N (−1)ψ00T 0 (1)ψŘ1 T 0 (−1)ψĽ1 T 0 (ζ1 )···

T N (1)
ψŘ1 T N (−1)ψĽ1 T N (ζ1 )00Ř1D(2)A2 A∗ P̌(1)31 (−1) + AP̌31 (−1)(2)(1)
A2 A∗ Ň31(ζ1 ) + AŇ31(ζ1 )······

Ľ1D(2)A2 A∗ Š(1)31 + AŠ31 − R11=0(91)··············ψψ(1)(2)
Ľ1 T 0 (ζN−2 ) · · · Ľ1 T N (ζN−2 ) Ľ1D A2 A∗ Ň31(ζN−2 ) + AŇ31(ζN−2 )上式から次の Landau 方程式が得られる．

dA= λ0 A + λ1 |A|2 A
dT 1(92)√
ここで λ0 は線形増幅率，λ1 は Landau 定数である．式 (92) の平衡解（平衡振幅）は −λ0 /λ1 と表される．λ1 < 0 の時，λ0 > 0 の領域では実数の平衡解が存在する超臨界分岐 (supercritical bifurcation) となるが，λ1 > 0 の時，λ0 < 0 の領域にしか実数の平衡解が存在しない亜臨界分岐 (subcritical bifurcation) となる（図–6 参照）．

7.

結果と考察

表–1 に線形解析より得られた臨界フルード数および弱非線形解析より得られた Landau 定数を示す．
表中 1–3 行目に C −1 = 21，µ = 0.1，m = 1 の場合について示した．このとき臨界フルード数は 0.8256であり，それに対応する波数（臨界波数）は 0.227 である．表より臨界波数周辺ではいずれも λ0 > 0 お

10よび Re(λ1 ) < 0 となっていることが分かる．これはデューン–平坦床遷移が超臨界分岐であることを意味している．C −1 のみを 22 に増加させた場合の結果を表中 4–6 行目に示した．C −1 が大きくなると，臨界波数およびそれより小さな波数領域において λ0 > 0 および Re(λ1 ) > 0 となっており，デューン–平坦床遷移は亜臨界分岐となる．
次に，局所勾配の影響を表すパラメータ µ を小さくした場合と，相当粗度高さ k s と粒径 d s の比 m を大きくした場合の結果を表中 7–12 行目に示した．µ が小さくなると Landau 定数の実部は負の値を取り，m が大きくなると正の値を取る．局所勾配の影響が小さくなると超臨界分岐となり，粒径に対する相当粗度が大きくなると亜臨界分岐となることがわかる

8.

結論

開水路乱流を良好に記述できると言われている混合距離仮説を用いて，河床デューンの非線形安定解析を行った．その結果，デューン–平坦床遷移の分岐形態について次のようなことが明らかとなった．

• 摩擦係数 C −1 が大きい領域で亜臨界分岐となる．
• 局所勾配の影響が小さいと超臨界分岐となり，粒径に対する相当粗度が大きくなると亜臨界分岐になる．
参考文献
1) Engelund, F.: Instability of erodible beds, J. Fluid Mech., 42, pp. 225–244, 1970.
2) Fredsøe, J.: On the development of dunes in erodible channels, J. Fluid Mech., 64, pp. 1–16, 1974.
3) 山口里実, 泉 典洋: デューン–平坦床遷移過程にみられる亜臨界分岐現象, 土木学会論文集 No. 740/II–64, pp. 75–94, 2003.
4) 泉 典洋, 山口里実: デューン–平坦床遷移再考, 土木学会論文集 B, Vol. 62, No. 4, pp. 360–375, 2006.
5) Colombini, M.: Revisiting the linear theory of sand dune formation, J. Fluid Mech., 502, pp. 1–16, 2004.

11