k − ε乱流モデルを用いた中層密度流の解析
水工学研究室指導教官

増戸洋幸

細山田 得三福嶋

1. はじめにSurface

洪水時に河川流域から生産される多量の微細砂

ρ1（200µm）を含む濁水が成層した貯水池に流入すると，密ρ1＋Δρbm

度躍層に到達した濁水はそれ自身の密度分布と貯水池水流入流体

uihの密度分布の相対的な関係により，流入濁水の一部またはすべてが躍層界面に沿って水平に向きを変えて，楔状祐介

ρ1＋Δρ0zに貫入する．これを中層密度流と呼ぶ．中層密度流は，濁水長期化の主たる原因として挙げられ，下流水域の環x

図− 1 中層密度流の模式図

境問題や堆砂問題と密接に関係している．中層密度流の流動特性を明確にすることにより，このような問題に対○連続式して効果的な対策を講じることが期待できる．

2. 研究の目的
Kao（1977）は理論的に中層密度流の先端移動速度を求めているが，非粘性流体の仮定に基づいているため，実現象を精度良く再現できてはいない．そこで本研究
では，乱流運動エネルギー k と分子粘性逸散率 ε を未∂ρ ∂ρu ∂ρw++=0∂t∂x∂z(1)○ x 方向の運動量方程式∂p∂(ρu) ∂(ρu2 ) ∂(ρuw)++=−∂t∂z ¶ ∂xµ ∂x∂ 2 (ρu) ∂ 2 (ρu)∂∂+ν+−(ρu02 ) −(ρu0 w0 )∂x2
∂z 2∂x∂z(2)知数とする k − ε 乱流モデルを用いて中層密度流の数値解析を行った．それぞれの輸送方程式の離散化手法と○ z 方向の運動量方程式

これは移流拡散流束に対する独自の離散化手法，圧力方∂p∂(ρw) ∂(ρuw) ∂(ρw2 )++= −ρg −∂t µ∂x∂z ¶∂z∂ 2 (ρw) ∂ 2 (ρw)∂∂
+ν+−(ρu0 w0 ) −(ρw02 )
∂x2
∂z 2
∂x∂z

程式の導入により，基礎変数を扱っているのにも関わら(3)しては S.V.Patanker（1985）の SIMPLE（Semi-ImplicitMethod for Pressure-Lincked Equations）法を採用した．

ず，その取り扱いが容易であり実用的である．
数値計算プログラムの妥当性を検証するために福嶋（1981）の行った二次元貯水池密度流の実験結果と比較を行った．その後，流入流体と下層流体の密度を変化させた中層密度流貫入実験の実験結果との比較を行い，計算モデルの妥当性を検証した．次いで，中層密度流の流動特性の解明のため濃度コンターと流速ベクトルを求めた．

3. 解析モデル
図− 1 のように，上層が ρ1 下層が ρ1 + ∆ρ0 の密度をもつ二層界面にこれらの中間の密度 ρ1 + ∆ρbm をもつ流体が進入する場合を考える．流れは鉛直二次元であり，十分に発達した乱流であるものとする．密度は塩分濃度c のみの関数と仮定する．数値計算に用いた基礎方程式を以下に示す．

○拡散方程式µ 2¶∂ c∂2c∂c ∂(uc) ∂(wc)++=D+ 2
∂t∂x
∂z
∂x2∂z
∂ 0 0∂ 0 0− cu −cw∂x∂z○乱流運動エネルギー k の方程式µ¶∂k ∂(uk) ∂(wk)∂νt ∂k++=∂t∂x∂z∂x σk ∂xµ¶µ ¶2∂
νt ∂k∂uνt ∂c++ νt+ Rg−ε∂z σk ∂z∂zσt ∂z(4)(5)○分子粘性逸散率 ε の方程式µ¶µ¶∂ε ∂(uε) ∂(wε)∂νt ∂ε∂νt ∂ε++=+∂t∂x∂z∂x σε ∂x∂z σε ∂z( µ ¶)2ε∂uνt ∂cε2+c1ε
νt+ (1 + c3ε )Rg− c2εk∂zσt ∂zk(6)Surfacez2690cm
PUMPHEAD
TANK

P90cm

z1600cmz(j)zRESERVOR

ρ1
ρ1+ρ0xtanθ=0.180cm

Bottomzx(i)※ () is the number of gridsx

図− 2 福嶋（1981）の実験装置
HEAD
TANK

図− 4 福嶋（1981）の計算格子図
Surface

400cm
Pz2

PUMPz(j)ρ1
ρ1+ρ0zRESERVORx60cm

z1
Bottomzx(i)※ () is the number of gridsx

図− 3 中層密度流貫入実験の実験装置

図− 5 中層密度流貫入実験の計算格子図

ここで νt は渦動粘性係数であり，次のように定義される．

15cm，高さ 60cm のアクリル製水槽である．密度流流出k2
νt = cµε部には流出口が高さ 4.0cm ×幅 13.6cm，ゲート開閉式(7)のアクリル製の箱が用いられた．上層流体は淡水を使用し，下層流体は一様密度に調整した塩水が用いられた．

u は流下方向の流速，w は鉛直方向の流速，c は塩分のウラニン色素で黄緑色に着色した塩水を一定ヘッドタン

体積濃度，g は重力加速度，ρ は計算領域中の密度，Rクに導き，流量制御バルブを通してゲート開閉式の箱には相対密度差，D は拡散係数，σt は乱流シュミット数，貯められ，ゲートを引き上げることにより密度流が発生

σk ，σε ，σµ ，σ1ε ，σ2ε ，σ3ε は k − ε 乱流モデルの数値された．実験は，下層流体の塩分濃度を固定し流入流体定数である．の塩分濃度を変化させるパターン（Run1-1，1-2，1-3）と，流入流体の塩分濃度を固定し下層流体の塩分濃度を4. 室内実験について変化させるパターン（Run2-1，2-2，2-3）について行わ4-1. 福嶋（1981）の実験概要れた．

実験装置を図− 2 に示す．実験水路は長さ 700cm，幅水槽と平行に 2 台のデジタルビデオカメラが設置され，20cm，水深 80cm，底勾配 1/10 の片面ガラス張り水槽

1 台目のカメラで流出口から x = 70cm の位置まで，2 台である．上層流体は淡水を使用し，下層流体は一様密度目のカメラで x = 70cm から x = 160cm の位置まで撮影に調整した塩水を用い，淡塩水の密度成層場が形成されされた．撮影された画像を画像解析することにより，流た．フルオレイセンナトリウムで黄緑色に着色した塩水

入水の流動形態，厚さ，先端形状，先端移動速度の場所を一定ヘッドタンクに導き，流量制御バルブを通して上的，時間的変化が求められた．

流端の堰より流入させる．堰の近傍では激しい混合が生じ流れが不安定なため，測定区間は堰より 90cm 離れた

5. 数値計算法について地点より下流とされた．

5-1. 福嶋（1981）との比較計算概要

水槽と平行に 35mm カメラを数台設置し，着色した塩

計算領域，計算格子を図− 4 に示す．計算領域は x 軸方

水の先端が測定区間に達したときを T =0sec とし，10 秒

向 410cm，z 軸方向 80cm とし，格子間隔は ∆x = 1.0cm，間隔で撮影された．同時に 8mm カメラにより塩水の先∆z = 0.5cm とした．底面は 1/10 勾配の階段状とした．

端を追跡した．このようにして得られた 35mm フィルム

計算時間は時間ステップを ∆t = 0.1sec として，140secと 8mm フィルムから流入水の流動形態，厚さ，先端形まで行った．流入部境界条件は各々の実験条件に従い与状，先端移動速度の場所的，時間的変化が求められた．

えた．以上の条件のもとで数値計算を行い，中層密度流4-2. 中層密度流貫入実験の実験概要

先端の移動速度，先端部厚さ，中層密度流厚の実験値と実験装置を図− 3 に示す．実験水路は長さ 400cm，幅比較をした．

5-2. 中層密度流貫入実験との比較計算概要
6向 170cm，z 軸方向 40cm とし，格子間隔は ∆x = 0.5cm，5∆z = 0.5cm とした．計算時間は時間ステップを ∆t =0.1sec として，40sec まで行った．流入部境界条件は各々の実験条件に従い与えた．以上の条件のもとで数値計算を行い，中層密度流先端の移動速度，先端部厚さの実験移動速度(cm/sec)

計算領域，計算格子を図− 5 に示す．計算領域は x 軸方計算値
実験値

4321値と比較をした．また，得られた数値計算結果から濃度0コンター図と流速ベクトル図を作成した．

050100150200x(cm)250300

図− 6 先端移動速度の流下距離変化

6. 結果と考察
6-1. 福嶋（1981）との比較
図− 6 は中層密度流先端の移動速度の流下距離変化で25ある．実験値を点で，計算値が線として示している．中であるという実験結果を解析結果がよく再現している．

15

図− 7 は中層密度流の先端部厚さの流下距離変化であ厚さ(cm)

層密度流先端の移動速度は，距離的に変化せずほぼ一定

20る．先端部厚さは，流下初期（x < 100cm）において計算計算値
実験値10値が大きくなる結果となった．実験の場合には x = 0cm5において放物線状の速度分布を持っていると考えられる00が，計算では流入部条件として速度を断面内で一様とし50て与えているため，流入開始直後に実験値より広く拡散

100150200x(cm)250300

図− 7 先端部厚さの流下距離変化してしまっていると考えられる．x > 100cm は，実験結果と計算結果とでよく一致しているといえる．
35厚の流下方向変化である．中層密度流厚は，先端部に膨30らみをもち急激に閉塞した後徐々に広がるという実験結25果と同様の傾向は得られたが，先端部より上流において計算値が実験値より大きくなる結果となった．実験では厚さ(cm)

図− 8 は先端位置が x = 280cm における中層密度流計算値
実験値

201510密度流の厚さを染料によって可視化された塩水を画像解5析することにより求めているが，数値計算では塩分濃度0を基準にし求めているため，実験で可視できている領域050100よりも濃度の領域の方が広かったと考えられる．

150x(cm)

200250300

図− 8 中層密度流厚の流下方向変化（x = 280cm）

図− 9 は計算開始から 94.5 秒後の中層密度流先端付近における流速ベクトルと濃度コンター図を重ねて描いたものある．図− 9 の流速ベクトル，濃度コンター図を見

80

行方向上向きの流速ベクトルが発生している．また，中60

層密度流先端上部を中心とした渦状の流れが形成されてz(cm)

ると，中層密度流前方で内部孤立波の発生へと繋がる進

2.0cm/sec

0.00.20.40.60.81.040

いる．中層密度流先端部の形状は上下対称形となってお20

り，実験結果をよく説明している．

0150

200

250x(cm)

300

350
Time=94.5sec

図− 9 流速ベクトルと濃度コンター図

6-2. 中層密度流貫入実験との比較

10

計算値
実験値

下距離変化である．実験値を点で，計算値が線として示している．計算結果は実験結果をよく再現している．また，その他の実験条件でも計算結果と実験結果はよく一致した．

移動速度(cm/sec)

図− 10 は Run1-3 の中層密度流先端の移動速度の流864

図− 11 は Run1-3 の中層密度流の先端部厚さの流下2距離変化である．計算結果は実験結果をよく再現してい00る．また，流入流体の濃度を変えた場合，厚さの違いは50100150x(cm)ほとんどみられなかった．下層流体の濃度を変えた場合

図− 10 先端移動速度の流下距離変化

は，濃度 2%のとき厚さが最も厚く，3%と 4%の場合は（Run1-3，中層:3.0%−下層:4.0%）

顕著な差がみられなかった．
8

図− 12 は Run1-1 の先端位置が x = 130cm における数値計算結果の流速ベクトルと濃度コンターを重ねて描実験画像である．図− 12 を見ると上部は緩やかにふくらみ，下部は滑らかな上下非対称の形状をなしている．

6厚さ(cm)いたものである．図− 13 は同じ位置における Run1-1 の計算値
実験値

75432

この数値計算結果は染料で可視化した図− 13 の実験画1像とよく一致する．また流速ベクトルは先端部で巻き上00げ，上部で進行方向と逆向きの流れ，先端後方で巻き込50みの流れが形成されており，全体として循環流が形成さ（Run1-3，中層:3.0%−下層:4.0%）

図− 14 は Run1-3 の先端位置が x = 130cm における数値計算結果の流速ベクトルと濃度コンターを重ねて描30いたものである．先端部は Run1-1 の場合より上部のふくらみがなだらかになり，下部によりふくらみをもつ形

5.0cm/s ec

0.00.51.025z(cm)る循環流はみられるものの，流速ベクトルは小さい．ま150

図− 11 先端部厚さの流下距離変化

れている．

状をなしている．Run1-1 でみられたような上層におけ100x(cm)た，下層において進行方向とは逆向きの流れが発生して2015いる．
10100以上より，流入流体の密度が増加すると先端部の形状は下部にふくらみをもちはじめ，上層の流速ベクトルに

140
1.0‑4.0Time=26.0sec

30る．先端部の形状については，流入流体の密度が増加すると鉛直方向の圧力が増加するため上部のふくらみはな

中層密度流に連行される周囲水量が上層では減少し下層

130（Run1-1，中層:1.0%−下層:4.0%）は逆向きの流速ベクトルが発生し大きくなることがわか

25z(cm)

流速ベクトルについては，鉛直方向の圧力の増加により120x(cm)

図− 12 流速ベクトルと濃度コンター図

よる循環流は徐々に弱くなり，下層において進行方向とくなり下部にふくらみが形成されると考えられる．また

110

2015では増加するため，それを補償する流れに変化が生じたと考えられる．

10100

110

120x(cm)

130

図− 13 実験画像
（Run1-1，中層:1.0%−下層:4.0%）140図− 15 は Run2-1 の先端位置が x = 130cm における数値計算結果の流速ベクトルと濃度コンターを重ねて描305.0cm/s ec

0.00.51.0

いたものである．先端部は下部に大きなふくらみをなし
25

果と同様に上層において循環流が発生しており，下層においても僅かではあるが循環流が形成されている．

z(cm)ており，実験結果をよく説明する．また他の条件での結2015

図− 16 は Run2-3 の先端位置が x = 130cm における数値計算結果の流速ベクトルと濃度コンターを重ねて描

10100いたものである．Run2-3 の数値計算結果は Run2-1 の110120x(cm)1301403.0‑4.0

場合より下部のふくらみがなくなっている．上層におい図− 14 流速ベクトルと濃度コンター図て発生している循環流は Run2-1 の場合より大きくなり，（Run1-3，中層:3.0%−下層:4.0%）

Time=23.0sec

下層における流速ベクトルは小さくなっている．
以上より，下層流体の密度が増加すると先端部の形状

305.0cm/s ec

0.00.51.0は下部のふくらみが徐々になくなり，上層の循環流は強25ることがわかる．先端部の形状は，下層流体の密度が増加することにより流入流体が下層流体を押しのけることz(cm)くなり，下層における逆向きの流速ベクトルは小さくな

2015ができないため下部のふくらみが減少したと考えられる．
流速ベクトルについては，中層密度流が押し上げられた

10100

ことにより下層では連行される水量が減少し，逆に上層

110120x(cm)1301401.5‑2.0Time=29.5sec

図− 15 流速ベクトルと濃度コンター図では連行される水量が増加したため流れが変化したと考（Run2-1，中層:1.5%−下層:2.0%）

えられる．
理論的検討として成層と乱流との比をあらわすパラメーターであるリチャードソン数を流入流体の密度を考

305.0cm/s ec

0.00.51.0

慮できるよう式 (8) のように拡張し，計算結果の整理を
25

Ri0 =

(ρ2 − ρbm )gh02
ρ1 Uih0(8)ここで，ρ1 は上層密度，ρ2 は下層密度，ρbm は流入流体の密度，h0 は流動層厚さ（= 40cm），Uih0 は流入流体z(cm)行った．

201510100110120x(cm)1301401.5‑4.0の初期速度（= 3.0cm/sec）である．図− 17 は中層密度

図− 16 流速ベクトルと濃度コンター図

流の上流側の層厚をリチャードソン数に対して整理した（Run2-3，中層:1.5%−下層:4.0%）

ものである．横軸はリチャードソン数，縦軸は上流側のTime=24.0sec

100

層厚 δih を流動層厚さ h0 で無次元化している．リチャーなる結果となった．リチャードソン数が大きくなると密度成層の乱れが小さくなるので，中層密度流の層厚が小δih/h0

ドソン数が大きくなるにしたがって無次元層厚が小さくさくなると考えられる．この結果から，実際の貯水池においても温度躍層と流入濁水の関係からリチャードソン数を求めることにより，中層密度流の層厚が推定可能であると考えられる．

10‑1100101102
Richardson No. Ri0103

図− 17 中層密度流の層厚とリチャードソン数の関係7. 結論
本研究では濁水長期化の主たる原因である中層密度流について，その流動特性を解明するため，乱流特性量の解析に k − ε 乱流モデルを，離散化手法に SIMPLE 法を用いて数値計算を行った．そして数値計算結果と実験結果との比較検討を行った．その結果を以下にまとめる．

参考文献

1) Timothy W. Kao：Density Currents and Their Applications，Journal of the Hydraulics Div.，ASCE，Vol.103，No.HY.5，pp.543〜555，1977．
2) パタンカー, S. V.（水谷幸夫・香月正司共訳）：コンピュータによる熱移動と流れの数値解析, 北出版,

• 先端移動速度と先端部厚さについて，数値計算結果

1985.

は実験結果とよく一致した．また，流入流体と下層流体の密度の変化による傾向の違いがみられた．

3) 福嶋祐介：成層化した貯水池に流入する濁水の流動機構に関する基礎的研究，東京大学学位論文，1981.

• 流動特性に関して，流入流体の密度が増加すると先端部の形状は下部にふくらみをもちはじめ，上層の流速ベクトルによる循環流は徐々に弱くなり，下層において進行方向とは逆向きの流速ベクトルが発生し大きくなることがわかった．また，下層流体の密度が増加する場合，先端部の形状は下部のふくらみが徐々になくなり，上層の循環流は強くなり，下層における逆向きの流速ベクトルは小さくなることが示された．

• リチャードソン数が大きくなるにしたがって層厚が小さくなる結果となった．これにより，実際の貯水池においても温度躍層と流入濁水の関係からリチャードソン数を求めることにより，中層密度流の層厚推定が可能であると考えられる．
以上の結果より，本研究で作成した数値計算モデルによる中層密度流解析の妥当性が示された．しかしながらいくつかの残された問題がある．以下にそれらの問題を列挙する．

1. 本研究では流入濁質の第一近似として塩水を対象としているため，物質の沈降性は考慮していない．しかしながら自然界において，流入濁質は沈降速度が無視できない程度の粒子も当然含んでおり，それらは下層密度流により躍層界面まで輸送され，中層密度流によりさらに貯水池内へ拡散される．したがって，沈降性濁質物質の挙動を考慮するモデルの導入が必要である．

2. 実際の貯水池内では，貯水池への潜り込みから下層密度流への遷移，そして躍層界面における分岐を経て中層密度流となる．実際の貯水池内での中層密度流の流動特性を解明するために，これら複合的な現象を考慮するモデルの構築が不可欠である．