水工学論文集,第52巻,2008年2月

LES乱流モデルを用いた水平流動する単一粒径・二粒径混合粒子サーマルの数値解析
NUMERICAL SIMULATIONS OF UNIFORM OR MIXTURE SIZED PARTICLETHERMALS ON HORIZONTAL BED USING LES TURBULENCE MODEL

重枝 未玲1・秋山 壽一郎2・杉山 誉3
Mirei SHIGE-EDA, Juichiro AKIYAMA and Takashi SUGIYAMA
1正会員 博士(工) 九州工業大学准教授 工学部建設社会工学科(〒804-8550 北九州市戸畑区仙水町1-1)
2フェロー会員 Ph.D. 九州工業大学教授 工学部建設社会工学科(同上)3学生会員 九州工業大学大学院

工学研究科建設社会工学専攻(同上)Numerical simulations of uniform or mixture sized particle thermals on horizontal bed were performed
using 3D numerical model. The thermal consists of fine and coarse sized particles and water. Thenumerical model is constructed with Simplified Marker and Cell (SMAC) method, Monotone Upstreamcentered Schemes for Conservation Laws (MUSCL) technique, Smagorinsky turbulence model andtreatment of several diameter sizes of a particle. The simulated results are compared with experimental
data of maximum height, propagation speed and amount of particles deposition of particle thermal. It isfound that this model can reproduce behavior of mixture sized particle thermals on horizontal bed.
Key Words : Particle thermal, fine and coarse particle, numerical models, SMAC method, MUSCLmethod, Smagorinsky model

１． はじめに

粒子サーマル現象は，水門の開放による濁水の清水域への広がり，土砂直投工による濁りの発生など，水質汚濁現象と密接に関係した固-液混相乱流現象である．このため，水工水理学をはじめとする様々な分野で活発に研究がなされている1)，2)，3)．
底面に沿って流動する粒子サーマル現象は，土砂直投工などでみられる現象である．その挙動は，初期条件(初期総有効重力，粒子粒径，濁水塊の形状)，底面の状態(勾配，底面粗度)，水域の条件(水域水深，密度，流れ)の影響を強く受ける．また，底面での粒子の巻き上げやサーマルを構成する浮遊粒子の底面への堆積により，駆動力である総有効重力の増加／減少が生じ，サーマルの加速あるいは減速が発生する．これらの要素を適切に取り扱うことが，巻き上げ・堆積を伴う粒子サーマルの挙動を予測するために不可欠である．
このような粒子サーマルの挙動の予測を目的に，様々な数値モデルが構築・提案されている4)~10)．粒子の堆積メカニズムを考慮した層平均モデル4)，5)，6)，巻き上げ・堆積を考慮した層平均モデル7)，鉛直2次元モデルや3次元モデル8)，9)，10)がある．また，粒子サーマルは高濃度の固-液混相乱流現象であるため，個々の粒子をLagrange的に取り扱う手法8)なども提案されている．しかし，実スケール海域のような大きなスケールを対象とする場合，個々の粒子の挙動を追跡することは現在のコンピュータの能力では難しく，現状では1流体モデル9)，10)としてEuler的に取り扱われるのが通常である．
著者らは，これまでに単一粒径粒子で構成される粒子サーマルが水平流動する現象に，3次元モデルを適用し，実験結果に基づき数値モデルの予測精度について検証してきた11)，12)，13)．しかしながら，実現象では粒子群は粒度分布を有しており，その適用には粒度分布が粒子サーマルの挙動に及ぼす影響を適切に取り扱えるモデルが必要となる．
本研究では，以上のような背景を踏まえ，著者らが開発 し た SMAC(Simplified MAC method) 法 14) と MUSCL(Monotone Upstream-centered Schemes for ConservationLaws)法15)に基づくモデル13)に，LES乱流モデルを導入するとともに，混合粒径粒子の取り扱いを考慮した3次元数値モデルへと発展させた．さらに，単一粒径・二粒径混合粒子で構成される水平底面を流動する粒子サーマルの実験結果に基づき同モデルを検証した．

ここでは，Cs=0.15を用いている．また式(7)中の| S | と，２．モデルの概要

GS成分歪速度テンソル S lm は式(8)，(9)のように定義される．

個体粒子群については，その粒径および沈降速度は小さく，粒子の速度と流体の鉛直速度は沈降速度だけ異なると仮定する1流体モデルとして取り扱う．

(S lm =(1) 基礎方程式
基礎方程式は，連続の式，ブシネスク近似を施した運動方程式および粒子の沈降速度を考慮した体積濃度の移流拡散方程式である．乱流モデルとしてLESを用い，フィルタリング操作を施すことにより，格子スケール以上の成分と格子スケール以下の成分とに分離する．
フィルター操作を施した基礎方程式は，それぞれ式(1)，(2)，(3)で表される．
∂U l=0(1)

∂U l
∂U l
∂ 2U l ∂ u l′u m′
1 ∂p+Um=−+ν−+ εg l∂t∂x m
ρ 0 ∂x l∂x m2∂x m(2)∂ c i ∂ (U l + W si )c i
∂ 2 c i ∂ u l′c i′+= Dm−∂t∂xl∂ x l2∂xl(3)∂xl

nni =1i =1

C = ∑ c i ； ε = ∑ ci ⋅ s i(4)ここに，t=時間，l ，m=1，2，3であり，それぞれ流下，横断，鉛直方向に対応する添え字，U=l方向の流速，u’=l方向の流速のsub-gridスケール成分，p=圧力，C=粒子の体積濃度，ci=各粒径に対応する粒子の体積濃度，c’=体積濃度のsub-gridスケール成分，ε =相対密度差(=(ρ-ρa)/ρa)，ρ =粒子群の密度，s=粒子の水中比重，ρa=周囲水の密度，g=重力加速度，ν=水の動粘性係数，Dm=物質拡散係数(=ν/Sc)，Sc=シュミット数，i=粒子粒径に対する添え字で総和規約は適用しない．また，式(3)第2項中のWsは粒子の沈降速度でありRubeyの実験式から算出される．式(2)，(3)中の ul′u ′m ， u l′c l′ は，式(5)，(6)のように表せる．
 ∂U ∂U m  2
 − kδ lm
ul′u m′ = ν t  l +∂xl  3
 ∂xm

− ul′ci′ =ν t ∂C(5)(6)

S c ∂xl

ここに，νt=sub-gridスケールでの渦動粘性係数，k=乱れエネルギー，δlm=クロネッカーデルタ関数である．
渦動粘性係数νtには修正Smagorinskyの渦動粘性モデル11)を用いた．渦動粘性係数νtは式(7)のように表される．
ν t = (C s ∆ )  S −2

2g l ∂∆ρ 
ρS c ∂xl (7)ここに，Δ =フィルター幅，Cs=Smagorinsky定数であり，)(8)

1  ∂U l ∂U m +∂xl 2  ∂xm(9)| S |= 2S lm S lm

1/2(2) 数値解析手法
基礎方程式の離散化は，SMAC法15)に基づき行う．解析手法の詳細は，参考文献13)を参照されたい．以下に本計算手順の概要を示す．
a) まず，時間tでの流速ul，圧力p，全粒子の体積濃度Cを用いて，運動方程式から流速の予測子を求める．
同時に，各粒子粒径に対応する体積濃度ciの移流拡散方程式から，新しい時間ステップの粒子の体積濃度ciを求める．
b) 圧力の時間変化に対応したスカラーポテンシャルφのポアソン方程式をSOR法で解き，圧力と流速を更新し，新しい時間ステップの圧力と流速を求める．
移流項の離散化にはMUSCL法16)を，圧力項，粘性項および連続の式については中心差分により離散化を行った．計算格子にはスタガード格子を用いた．
(3) 境界条件
側壁，底面および水表面での境界条件はそれぞれ次式のように与えた．
・ 側壁と底面：ul=0(non-slip条件)，∂φ/∂n=0，∂ci/∂n=0
・ 水表面：∂u3/∂x3=∂u2/∂x3=0，u3=0(slip条件)，φ=0，∂ci/∂n=0
ここに，n=水表面，側壁および底面に垂直な方向を表す．
浮遊粒子の底面への堆積は，沈降粒子濃度フラック
スDを式(3)の境界条件として与え，これを底面から流出させることで表現した．なお，ここでは堆積した粒子の巻き上げはないとして取り扱っている．
底面近傍の沈降粒子濃度フラックスを表現するにあたり，基準点濃度や粒子の沈降速度の取り方などについては様々な提案がなされているものの，開水路浮遊砂流に対してさえも未だ明確な見解を得るに至っていない．また，本研究で取り扱っているようなサスペンションタイプの密度流については，十分な検証がなされないまま開水路浮遊砂流に関する知見を準用しているのが現状である．そこで，ここでは式(10) に示す沈降粒子濃度フラックスDを用いた．
D = αW s C b(10)ここに，α=底面への堆積率を規定するパラメーター，表-1 実験条件
粒径各粒径の 沈降速度 水深
Vf(m/s) h(m)d(m)
割合0.000044
1
Case-A
0.00155
0.00010900.000044
0.5
Case-B
0.00458 0.250.000109
0.50.000044
0
Case-C0.00886
0.000109
1Case

初期総有効重力 初期体積423W0(cm /s )A0(cm )

1754753610

図-2 撮影領域と堆積粒子の測定点

図-1 定義図

Ws= 粒子の自由沈降速度およびCb=底面近傍での粒子濃度である．沈降粒子濃度フラックスDに関する係数αは，浮遊粒子の沈積を伴う3次元濁水サージに関するデータが存在しない．ここでは，式(10) において，底面に最も近い計算セルの重心での粒子濃度をCbとして採用し，Ws としては中央粒径粒子の自由沈降速度を用いた．底面への堆積率を規定するパラメーターαについては，秋山ら12)が，LES 数値モデルにより底面衝突後の直投微細粒子群の挙動とその堆積プロセスに関する解析で用いたα=2.0を採用した．

３．実験の概要
実験は，直投濁水塊の底面衝突後の挙動を対象としている．実験装置は，長さ2.0m，幅2.0m，深さ0.6mの水槽であり，水槽の底面中央に長さ0.19m，幅0.19m，高さ0.10mの流入装置が設置されている．
濁水はガラスビーズでモデル化し，表-1のように初期総有効重力W0(cm4/s2)，初期体積A0(cm3)，周囲水水深h(m)はそれぞれ一定で，Case-A~Cのように混合比を変化させた．ここに，ε0=初期相対密度差，ρ0=濁水の密度，ρa=周囲水の密度およびg=重力加速度である．なお，粒子は，粒径d=0.000044mのものを細粒子，d=0.000109mのものを粗粒子と定義している．底面に設置した流入装置に，初期総有効重力W0の濁水を注入し，十分に攪拌した後，ゲートを瞬間的に開放することで3次元粒子サーマルを発生させた．濁水塊の流動特性量の定義は，
図-1に示す通りである．
測定項目は，各測定点での層厚H(m)，フロント位置および堆積粒子量Wd(g)である．層厚Hについては，図2に示すA～J 断面で，水表面上からスリット光を照射することで粒子サーマルを可視化し，デジタルビデオカメラで撮影した画像を解析することで求めた．その最大値を最大層厚Hmとした．フロント位置については，水槽上方から撮影した画像から特定した．フロント移
動速度Ufについてはフロント位置より算定した．堆積粒子量Wdは，図-2中に示す赤点を中心とした各0.2×0.2mの領域で堆積粒子を採集し，その乾燥重量より算定した．いずれの測定も3回行い実験値の信頼性を高めた．なお，いずれの実験も水槽の中心である原点を通る軸に対して概ね左右対称であったので，撮影領域は図-2に示す実験水槽の4分の1とした．

４．結果と考察
まず，単一粒径粒子で構成される粒子サーマルの実験結果に対して，次に二粒径混合粒子で構成される粒子サーマルの実験結果に対してモデルの検証を行う．
(1) 単一粒径粒子サーマルの実験結果に対する検証
図-3は，単一粒径粒子のCase-A，Case-Cの投下条件(A0，W0)で無次元化した無次元最大層厚Hm*(=Hm/(A01/3))，
無次元フロント移動速度Uf*(=Uf/(W01/2/A01/3))とx軸方向の無次元流下距離xf*(=xf/(A01/3))またはr軸方向の無次元流下距離rf*(=rf/(A0)1/3)との関係および最大堆積粒子量Wdmaxで無次元化した堆積粒子量Wd*の実験および解析結果の比較を示したものである．いずれの結果についても，参考までに乱流モデルを導入していない，単一粒径粒子を対象とした旧モデル13)の解析結果も併せて示している．
これらの図から，細粒子を用いたCase-Aでは，フロント移動速度，層厚，堆積粒子量のいずれについても解析結果は実験結果を再現できること，粗粒子を用いたCase-Cでは層厚が若干大きく評価されているものの，実験結果を概ね再現していること，などが確認できる．
また，本モデルは乱流モデルを組み込んでいない旧モデル13)と比較すると，特にフロント移動速度において精度が改善しており，本研究で導入した乱流モデルの有効性が認められる．これまでに，著者らは乱流モデルを組み込んでいないモデルを2次元や軸対象の粒子サーCase-A Exp x軸
Case-A 新モデル x軸
Case-A 旧モデル x軸

Case-C Exp x軸
Case-C 新モデル x軸
Case-C 旧モデル x軸

Case-A Exp r軸
Case-A 新モデル r軸
Case-A 旧モデル r軸

0.4

0.4

0.3
Hm*
0.2

0.3
Hm*
0.2

0.1

0.1

00

0.4
Uf*
0.3

0.4
Uf*
0.3

0.2

0.2

0.1

0.1

00

123 *
xf*，r
f4

05(a) 無次元最大層厚Hm*および無次元フロント移動速度
Uf*と無次元流下距離xf*との関係
0.0
0.5
1.0
Wd*(Wdmax=25.30g)

5.0y*

0.0
1.0
0.5
Wd*(Wdmax=26.82g)

0.8
0.6
0.4
0.2

2.5012Case-C Exp r軸
Case-C 新モデル r軸
Case-C 旧モデル r軸

3xf*，rf*456(c) 無次元最大層厚Hm*および無次元フロント移動速度
Uf*と無次元流下距離xf*との関係
0.0
0.5
1.0
Wd*(Wdmax=59.98g)

5.0

1.0
0.0
0.5
Wd*(Wdmax=75.24g)

0.8

y*

0.8

0.6

0.8
0.6

2.5

0.4

0.4
0.2

0.4

0.6

0.2

0.2

0

00

2.5x*

5.0

0

2.5x*

5.0(b) 無次元堆積粒子量Wd* (左：実験結果，右：計算結果)

02.5x*

5.002.5x*5.0(d) 無次元堆積粒子量Wd* (左：実験結果，右：計算結果)

図-3 単一粒径粒子サーマルの実験値と解析結果との比較 (左：Case-A，右：Case-C)

マルに適用し良い結果を得てきた．しかし，今回対象とした粒子サーマルは，3次元的に広がるものであり，2次元や軸対称の粒子サーマルに比べ境界により流れが規定されるものではない．そのため，流れの自由度がより大きくなり，計算格子で解像することが難しくなったため，LES乱流モデルを導入することで精度が改善したと考えられる．
以上の結果から，本モデルは，単一粒径粒子で構成される粒子サーマルの挙動を予測可能であることが確認された．
Case-B 0s(2) 二粒径混合粒子サーマルの実験結果に対する検証
図-4は，二粒径混合粒子で構成される濁水塊の広がりの実験および解析結果を，図-5は濁水塊の相対密度差および粗粒子成分と細粒子成分の相対密度差と流速ベクトルの解析結果を示したものである．
図-4より，実験と解析のいずれも広がりが軸対称であること，両結果を比較すると若干ではあるが解析値の広がりの方が早く，その形状もよりひし形に近いこCase-B 2s

Case-B 6s

図-4 流況写真と粒子サーマルの広がりの解析結果との比較(上：流況写真, 下：解析結果)と，などがわかる．図-5より，濁水塊が周囲水を連行しながら，サーマルを形成する様子，粗粒子成分の相対密度差εは，細粒子成分に比べ粒子の沈降堆積が早いため，早い時間で小さくなる様子，一方，細粒子成分の相対密度差εは周囲水と混合するため，粗粒子成分に(cm)20 Case-B 混合
15105020 Case-B 粗粒子成分
15105020 Case-B 細粒子成分
15105020010302s

40504s01020304050 (cm)

図-5 濁水塊の相対密度差コンターと流速ベクトル比べ拡散する様子，粗粒子成分のフロント移動速度が細粒子成分に比べ小さい様子などがわかる．また， 全体および細・粗粒子成分を比較すると，フロントが流下するにつれて，細粒子成分の相対密度差と全体の相対密度差は一致することが確認できる．そのため，この地点での濁水の挙動は細粒子成分により支配されていると考えられる．このように本モデルは，定性的ではあるが，粒径分布を有する水平面上の濁水塊の挙動を再現できることが確認できる．
図-6は，濁水塊の各特性量の無次元量(最大層厚Hm*，フロント移動速度Uf*)とx軸方向の流下距離xf*またはr軸方向の流下距離rf*との関係および無次元堆積粒子量Wd*について実験値と解析結果の比較を行ったものである．
なお堆積粒子量については，粗粒子成分Wdf*，細粒子成分Wds*についても併せて示している．
最大層厚Hm*とフロント移動速度Uf*については，解析値は流下に伴う減少を良好に再現できることがわかる．
図-3の結果からわかるように，単一粒径粒子では最大層厚Hm*とフロント移動速度Uf*の減少の割合は粗粒子成分が大きい．しかしながら，混合粒子では粗粒子に見られた急激な減少は認められない．これは，図-5で示したように，流下に伴い濁水を構成する粒子は細粒子成分になるためであると考えられる．このように本モデルは，粒径の違いが特性量に及ぼす影響を捉えていることがわかる．堆積粒子量Wd*については，解析結果は，粗粒子成分の堆積量を若干大きめに評価している
ものの，流入点直下の著しい堆積や流下距離の増加に伴う堆積量の減少など，実験結果を良好に再現できることがわかる．また，混合粒子群の粗粒子成分の堆積粒子量Wdf*と細粒子成分の堆積粒子量Wds*に着目すると，粗粒子成分については無次元堆積量0.2のコンターラインに，細粒子成分については無次元堆積量0.4のコンターラインに，若干のずれはあるものの実験結果を概ね再現できることがわかる．これについては，堆積粒子量の算定に用いた係数αは流れの状態によって変化すると予想されるが本解析では一定で取り扱っていること，本モデルでは堆積粒子の巻き上げを考慮していないこと，などが理由として考えられる．この点については今後検討したいと考えている．以上から，本モデ
ルは，定量的にも粒径分布を有する水平面上の濁水塊の挙動を再現できることが確認できる．

５．おわりに
本研究では，粒子サーマルを対象に，混合粒径粒子の取り扱いを考慮した3次元乱流数値モデルを構築し，単一粒径・二粒径混合粒子で構成される水平底面を流動する粒子サーマルの実験結果に基づき同モデルの検証を行った．その結果，本モデルが粒子サーマルの特性量である最大層厚，フロント移動速度および無次元堆積粒子量を十分な精度で再現可能であることがわかった．
ただし，本モデルは一流体モデルであるために，より粒径の大きな粒子を対象とした場合には予測精度が低下することが予想される．その場合には，粒子間の相互作用や固液間の相互作用などをモデル化する必要があると考えられるため，本モデルの適用限界について検討しておく必要がある．今後は，この点について検討するとともに，堆積粒子量の算定に用いた係数αについて検討し，堆積粒子の巻き上げのメカニズムや一般曲線座標などを導入することによって，より現実的な状況に適用可能なモデルに発展させたいと考えている．
謝辞：本研究を遂行するに当り，本学学部4年生の和田智之君には，データ整理等で多大な協力を得た．ここに記して感謝の意を表します．

Case-B Exp x軸
Case-B Sim x軸

参考文献

Case-B Exp r軸
Case-B Sim r軸

1)

0.4

Simpson, J. E.: Gravity Currents: In the environment and thelaboratory, Ellis Horwood, 1987.

0.3
H m*0.22)

Simpson, J. E.: Gravity currents in the laboratory, atmosphere andocean, Annual Review of Fluid Mechanics, Vol.14, pp.213-234,
1982.

0.1

3)

Annual Review of Fluid Mechanics, Vol.15, pp.47-76, 1983.

04)
0.4
Uf*Bonnecaze, R. T., Huppert, H. E. and Lister, J. R.：Particledriven gravity currents, Journal of Fluid Mechanics, Vol. 250,pp.339–369, 1993.

0.35)

Bonnecaze, R. T., Hallworth, M. A., Huppert, H. E. and Lister, J.
R.: Axisymmetric particle-driven gravity currents, Journal of0.2

Fluid Mechanics, Vol. 294, pp.93–121, 1995.

0.1
0Hopfinger, E. J.: Snow avalanche motion and related phenomena,6)01

23xf*，rf*

45Jha, A. K., Akiyama, J. and Ura, M.: Simulation of gravity
currents by flux-difference splitting scheme, Annual Journal of

6

Hydraulic Engineering, JSCE, Vol. 46, pp. 1037–1042, 2002.

(a) 無次元最大層厚Hm*および無次元フロント移動速度Uf*と無次元流下距離xf*との関係

7)

Bradford, S. and Katopodes, N. D.: Hydrodynamics of turbid
underflows. I:Formulation and numerical analysis, Journal of
Hydraulic Engineering, ASCE, Vol. 125, No. 10, pp. 1006–1015,

0.0
0.5
1.0
Wd*(Wdmax=46.10g)

5.0y*

0.0
1.0
0.5
Wd*(Wdmax=51.53g)

0.8

0.8
0.6
0.4
0.2

0.6

2.5

0.40.2

1999.
8)

小田一紀，重松孝昌，大西信幸，井上雅夫：改良型
DEMAC法による粒子群の沈降・分散挙動の解析，海岸工学論文集，第39巻，pp.971-975，1992．

9)

蒋勤，国栖広志，渡辺晃：土砂投入に伴う汚濁の拡散過程に関する数値シミュレーション，海岸工学論文集，第002.5x*

5.0

02.5x*

5.0

44巻，pp.1031-1035，1997．
10) 玉井昌宏，村岡浩爾：土砂直投工により発生する濁りの(b) 無次元堆積粒子量Wd* (左：実験結果，右：計算結果)

挙動特性に関する数値解析，水工学論文集，第42巻，

0.0
1.0
0.5
Wdf*(Wdfmax=37.38g)

11) 秋山壽一郎，応新亜，浦勝，重枝未玲：有限な水深を有pp.541-546，1998．
0.0
0.5
1.0
Wdf*(Wdfmax=30.61g)

5.0

y*

0.8
0.4

0.8
0.6
0.4

0.2

0.2

0.6

2.5

する静水中に直投された微細粒子群の数値シミュレー
ション，土木学会論文集，No.663/II-53，pp.43-53，2000．
12) Singh, J., Akiyama, J. and Shige-eda, M.: Numerical simulationof asymmetric gravity current induced by direct dumping ofsediments, Annual Journal of Hydraulic Engineering, Vol.49,

002.5x*

5.0

0

2.5x*

5.0(c) 無次元堆積粒子量Wdf* (左：実験結果，右：計算結果)

pp.1399-1404, 2005.
13) 秋山壽一郎，重枝未玲，前多陽子：3次元モデルによる水平面上の粒子サーマルの流動・堆積シミュレーション，
水工学論文集，第48巻，pp.1165-1170，2004．

y1.0
0.00.5
Wds*(Wdsmax=14.15g)

0.00.5
1.0
Wds*(Wdsmax=15.49g)

5.0
*

0.8

0.80.6

2.5

0.4
0.2x*

5.002.5x*

Computational Physics, Vol.6, pp.322-325, 1970.
15) Leer, van B.: Toward the ultimate conservative difference

0.4

scheme. 4 A new approach to numerical convection, Journal of

02.5

for incompressible fluid flow calculations, Journal of0.6

0.2014) Amsden, A. A. and Harlow, F. H.: A simplified MAC technique

5.0

Computational Physics, Vol.23, pp.276-299, 1977.
16) Rubey, W. N.: Settling velocities of gravel, sand and silt particles,
Amer. Jour. Sci., Vol.25, pp.325-338，1933．

(d) 無次元堆積粒子量Wds* (左：実験結果，右：計算結果)
図-6 混合粒径粒子サーマルの実験値と解析結果との比較（2007.9.30受付）