OppenFOA
AMによよるDyynamicc SGSモ
モデルをを用いいた
傾
傾斜壁面面粒子サ
サーマ
マルの数
数値解析
析衞藤 俊彦1・菊菊地卓郎郎2・大澤 範一3

1正正会員

8）
長岡工工業高等専門学学校准教授 環環境都市工学科科（〒940-8532新潟県長岡市市西片貝町888
E-maill:eto@nagaoka--ct.ac.jp
2正会員
員 福島工業高等専門学校准准教授 建設 環境工学科（〒970-8034福島県いわき市平平上荒川字長尾30）
ma-nct.ac.jp
E-mail:kikkuchi@fukushim
3正会員 東京
京建設コンサル
ルタント河川計
計画部（〒170-0004東京都豊豊島区北大塚1丁丁目15番6号）
E-mail:ohhsawa-n@token
ncon.co.jp

海底渓
渓谷や大規模なな湖底で発生すする泥水流は傾
傾斜壁面粒子ササーマルの流動動現象であり，，沿岸域の漂砂砂や渓谷の地地形変化，海岸岸浸食と密接にに関わっているる．本研究でははこの流動現象
象の予測を目的的として，解析析モデ
ルの構構築を行った．オープンソーースの三次元流流動解析ツール
ルOpenFOAMをを基に，これをを粒子サーマルの解を乱流計算に用いて析モデ
デルへと発展ささせた．この解解析モデルに， LESモデルとしてDynamiic SGSモデルをった．その結果，傾解析をを行い，既存のの実験結果と比比較することににより，解析モモデルの妥当性性の検証を行って，実験結果を良好
斜壁面
面粒子サーマルルの流下速度，最大厚さ，浮浮遊粒子総量のの流下方向距離離変化についてに再現現することがでできた．

mal, numerica
al analysis, LEES, Dynamic SGS
SKeey Words : OppenFOAM, incclined wall susspension therm1. はじめにに中に砂が浮遊した泥水が流入するとき，その均一な水中密度差によりり泥水流が形成成される．泥泥水流は海底渓
渓谷や大きな湖底でで発生し，沿岸岸域の砂や濁度物質を深部部に運man et搬する機構のの一部を担っていると考ええられる（Inmal.1)，Parker2)，Parker et al.3)）．そのため），海岸での漂漂砂や渓谷の地形変変化，海岸浸食
食と密接に関わる流動現象
象であり，その流動動機構を適切にに表現するシミュレーション手法の開発は，重要な研究テテーマである．
斜面上流かからの土砂の供供給が連続的である場合，定常なプルーム状状の泥水流が形形成され，上流からの土砂砂の供給が十分でなない場合，非定定常なサーマル状の泥水流流が形成される．後後者は傾斜壁面面粒子サーマルの流動現象象である．一定な傾傾斜角θをもつつ斜面上の傾斜
斜壁面粒子ササーマ
ルを模式的にに図-1に示す．上流から流入入してきた粒粒子を混合した流体体が，全体として非定常なフロント部をを形成しながら流下下する．
傾斜壁面粒粒子サーマルはは粒子の沈降性により，乱乱流強度が小さくななると浮遊状態態を維持できず，粒子が沈沈降しサーマルからら離脱する．ままた底面での粒粒子の巻き上上げがある場合は，底面からサーーマルに粒子が供給されるる．そのため負の浮浮力の総量は流流動に伴い増減し，サーママルの

図-1 傾斜壁面粒子ササーマルの模式
式図

度は減速またはは加速されるここととなる．傾傾斜壁面粒子子速度
サーーマルの流動をを予測するためめには，粒子の沈降や巻きき上げげの取り扱いがが非常に重要ととなる．
傾斜壁面粒子サ傾サーマルについいては，固体粒粒子の沈降をを考慮慮したサーマルル理論に基づく
く数値モデルによる解析がが4), 5), 6), 7)8)
行わわれている（福福嶋ら）．また衞藤ら
ら はk-ε乱流流モデデルを用いた鉛鉛直二次元モデデルを構築し，実験結果と比較較することによよりモデルの妥妥当性を検証した．さらにに9)
衞藤
藤ら はこのモデルを粒子のの底面からの巻巻き上げを考考慮ししたモデルへとと発展させ，流流動方向へ加速速する鉛直二二次元
元サーマルの流流動特性が明ららかにされた．しかしなががら自自然界で発生すする泥水流は，複雑な実地形形を流下し横横断方方向にも拡散すする為，その流流動機構を適切切に推定し施施

デルタである．νSGS はSGS渦動粘性係数であり，次式で設設計などに資するためには，三次元解析が必要不可欠である．
表す．
流体の三次元流動解析ツールとして，OpenFOAM10)が(6)
 SGS  Cs  2 Sある．OpenFOAMはサーマルのような混相流を含む様々
Δはフィルター幅，|S|は次式で表す．
な流動現象に応じたソルバを有する流体解析ツールであり，オープンソースであるため，容易に改良を加えるこ(7)S  2Sij Sijともできる．また多くの乱流モデルが利用でき，非構造格子を用いた複雑な計算領域に対する解析も可能である． 式(6)中に含まれるCsはSmagorinskyモデル定数である．こ加えて解析の際の前処理，後処理のためのアプリケーシのモデル定数は，Germano et al.11)のDynamic SGSモデルにョンが備わっており，並列計算にも対応するなど，非常に完成度が高い．
そこで本研究ではOpenFOAMを基に，傾斜壁面粒子サーマルの特徴である粒子の沈降，堆積を考慮できるように基礎方程式，境界条件を改良し，粒子サーマルの解析モデルへと発展させた．この解析モデルに， LESモデル
としてDynamic SGSモデルを乱流計算に用いて三次元流動解析を行い，傾斜壁面粒子サーマルの既存の実験結果と比較することにより，その妥当性の検証を行った．

基づき新たにテストフィルタを導入し，Lilly12)の最小二乗法を用いて動的に次式で計算される．

C1 Lij M ij
22 M ij M ij

2(8)式(8)右辺の〈 〉は空間平均操作を表す．またLij ，Mij は次式で与えられる．

~~~
Lij  ui u j  ui u j~~~
M ij   S Sij  S S

2. 解析モデルの概要

C  Cs

2ij

(9)~



(10)

(1) 基礎方程式
本研究では傾斜壁面粒子サーマルの解析に対して，OpenFOAMのソルバの中からtwoLiquidMixingFoamソルバ

チルダ付きはテストフィルタ操作を表す．γはテストフ
ィルタ幅とグリッドフィルタ幅の比でありGermanoの推奨値である2とした．

を採用した．これは非圧縮，混合性の2流体を対象としたソルバであり，基礎方程式は連続式，Navier-Stokes方程式，拡散方程式である．乱流モデルとしてLESモデルを用い，フィルタ操作を施した方程式を以下に示す．
ui(1)
0xi
ui ui u jp(2)
   ij  2S ij   g i(2) 解析手法
基礎方程式は有限体積法に基づき離散化し，PISO法13)による非定常計算を行う．移流項の離散化には二次精度のTVDスキームを採用し，制限関数にはsuperbeeを用いた．

tx jxi

xi

c ui cc   SGS c  Dm 2 txi
xi  S c xi xi2(3)ここでバー付きの変数はグリッドフィルタ操作された変数を示す． t は時間，添え字 i，j =1，2，3はそれぞれ流下方向，流下方向と直交する方向，横断方向成分，u は流速，c は初期粒子濃度で基準化された粒子濃度， ̅ は(3) 境界条件
境界条件は，壁面において流速はno-slip条件とし，粒子濃度は勾配0，乱流特性量は壁関数で与えた．水面においては，水面と直交方向の流速が0となる条件を与え，粒子濃度は勾配0の条件を与えた．

(4) 流動中の粒子の沈降，堆積の考慮
本研究で用いるtwoLiquidMixingFoamソルバは泥水流の圧力，ρは流体密度，g は重力加速度，νは動粘性係数，
Dmは拡散係数，Sc は乱流シュミット数を示す． τij はSGS流動機構に大きな影響をもたらす粒子の沈降を考慮できない．そのため，式(3)の拡散方程式に粒子の沈降速度応力，Sij はひずみ速度テンソルのGS成分であり，以下の式で表される．

 ij  2 SGS Sij  k ij

23(4)1  u u j Sij   i 2  x j xi (5)ここで k は乱流運動エネルギー，δ はクロネッカーのの項を加えた新たな式を導入し，解析で用いた．この式を以下に示す．
c  ui  Wsi c 2c   SGS c  Dm 2 xi  Sc xi txi

xi

(11)

ここでWs は粒子の沈降速度であり，ストークスの式より求める．
粒子濃度の底面への堆積は，重枝ら14)の水平流動する

粒子サーマルルの研究で用いいた境界条件を採用した．これ
は式(11)におおいて，沈降粒粒子濃度フラックスDを底底面から流出させるることで，粒子子の底面への堆積を表すももので．
ある．沈降粒粒子濃度フラッックスDは次式式で示される．

D  WsCb(12)αは底面へのの堆積率を規定定するパラメーターであり，重Cb は底面近傍枝らの研究とと同様にα=2..0を用いた．C
傍での粒子濃度であある．また解析
析では，堆積積した粒子の再再浮上については無無いものとする
る．

3. 数値解析析の概要(1) 実験概要
要解析モデル
ルの妥当性の検証のため，福嶋ら7)が行行った
硫酸バリウム
ム粒子を用いたた傾斜壁面粒
粒子サーマルのの室内
図-2 福嶋ら7)のの実験装置
実験結果を用用いた．この実実験は淡水が満たされた傾傾斜角30°）（上図：5°，10° 下図：35°，10°，30°の傾斜水水路に，それぞぞれ初期濃度
度を持の硫酸バリウムム粒子混合水を水路上流端端よりつ比重4.37の表-1 福嶋ら7)7)の実験条件
1L流入させるることで，水水路内に粒子ササーマルを発発生させその流動機機構を観察したたものである．図-2に実験
験装置
角 初期粒子濃濃度 (%) 初期相
相対密度差 R
傾斜角Case図を示す．水水路幅はどちちらの水路ともも15cmである
る．粒
110.00314
0子混合水は初初め，上流端のの流入ゲートを持つボッククスに230.0100
5350.0196貯められていいる．実験はは初期濃度1，3，5%でそれれぞれ
410.003540行われた．実実験条件を表--1に示す．表の初期相対密密度差530.0122
10Rは，粒子濃濃度に基づき設設定された粒子子重量から，実験
650.0219
後ボックス内内に取り残されれた粒子重量量を差し引いたた値，
710.003530
つまりサーママルとして流下下したバリウム重量より計計算さ830.0126
30れた値である
る．実験ではササーマルの流下速度，最大大厚さ
950.0235及び浮遊粒子子総量の流下方方向距離変化が測定されてている．
浮遊粒子総量は上流端ののボックスおおよび流下方方向に
表-2 解析
析条件20cm毎の区間間に沈降したたバリウム重量量を測定し，これを初期バリウウム重量から順順次差し引くことにより計計算さ
初期粒
粒子 粒径 沈沈降速度
傾
相対密度差 R0
初期相
Case 傾斜角%)(μm)(%
濃度(m/s)
れたものであある．実験で用用いられた硫酸バリウム粒粒子の11120.0.0002252
0.0003426152
粒径分布は22.0～50.8μmでである．
23453512830
150.0.0012129
0.0.0013897
0.0.0003514

0.0012144176
0.0023421147
0.0003973555(2) 解析条件
件
本研究は三三次元の領域をを対象として解解析を行ったた．座
1053300.0.0013897
0.0014814455
標軸は水路のの斜面方向をxx軸，それと直直交方向にy軸軸とし， 65300.0.0013897
0.0025699701
17200.0.0006229
0.0004153283
横断方向をzz軸とした．解解析領域の上流流端に硫酸ババリウ
8303340.0.0017772
0.0015608161
ム粒子混合水水の領域を与ええ，それ以外外は淡水領域ととし，
95370.0.0020972
0.0027610789
解析開始と同同時に上流端のの粒子混合水が流動を始めめる条件とした．
本解析モデルは
本は混合材料の粒粒径分布を考慮慮できないた
た
解析領域のの格子間隔は，，x軸，z軸は1cm刻みとし，，y軸め，代表粒径によよる単一粒径をを仮定し解析を実施する．
は傾斜角5°，10°では底底面から10cmままで0.5cm，傾傾斜角実験
験で用いられたたバリウム粒子子の代表粒径は不明である
る
30°は底面かから20cmまで00.5cm刻みとしし，それ以降はは1cmが，予備計算におおいて浮遊粒子子総量の流下方方向距離変化
化刻みとした．時間刻みははクーラン数がが0.5を越えなない範が実実験結果と一致致するような粒粒径を求め，これを解析でで囲で自動調整整させた．

図-3 流下速
速度，最大厚さ，浮遊粒子総量量の流下距離変変化（傾斜斜角5°）とした．その結結果，初期粒粒子濃度が大ききいほの設定粒径とど解析の代表表粒径も大きくくなる傾向が見られた．ここれは実験においててバリウム粒子子同士の凝集によって見かかけの粒径が大きくくなったためと考えられ，初期粒子濃度度が大きい条件ほどど凝集が顕著となり，実質質的な代表粒径径が大きくなったもものと考えられれる．
また解析ででの初期相対密密度差R0は，ボックス内部部の流体に設定するる必要がある．解析開始後，ボックス内内部に沈降する粒子子重量を考慮し，ボックスから流出したた粒子重量が実験とと解析で一致すするように予備備計算により検討し，初期相対対密度差R0を求求めた．
解析の条件件として，予備備計算により得られた代表
表粒径とその粒径ででの沈降速度おおよびボックス内部に設定定した
初期相対密度度差R0を表-2にに示す．

図-44 流下速度，最最大厚さ，浮遊遊粒子総量の流下下距離変化（傾斜角10°）4. 解析結果(1) 解析結果と実実験結果との比比較
図-3は傾斜角5図°の水路におおけるサーマル
ルの流下速度
度，最大大厚さ，浮遊粒粒子総量の流下下距離変化の比比較グラフでである
る．最大厚さははサーマル先端端部において実実験結果とのの一致致に配慮した上上で，各時間間における最大大粒子濃度のの0.3倍倍程度を境界部部の判定値として評価した
た．Case1，2のの解析
析結果は，実験験結果を概ねね再現している
る．Case3の解
解析結
結果は，流下速速度がやや小ささく評価されているもののの，減速速のタイミンググや浮遊粒子総総量の縦断的な変化においいて，実験結果の傾傾向は表せていいる．
0°の水路にお図-4は傾斜角10図
おけるサーマルルの流下速度
度，最大大厚さ，浮遊粒粒子総量の流下下距離変化の比比較グラフでで

図-6 粒子濃度ココンター，流速ベベクトル図（流流下距離100cm，5，下図：Case8）
上図：Case5）

図-5 流下下速度，最大厚ささ，浮遊粒子総総量の流下距離変変化（傾斜斜角30°）
図-7 サーマル
ル粒子濃度コンンターの三次元形状図ある．すべててのCaseにおいいて，解析結果
果は実験結果
果を非（流下距距離100cm，上図図：Case5，下図図：Case8）常によく再現現している．
図-5は傾斜斜角30°の水路路におけるサーマルの流下下速度，遊粒子総量が実実験結果と一致致するように粒粒子の代表粒
粒
浮遊最大厚さ，浮浮遊粒子総量のの流下距離変化の比較グラ
ラフで径をを解析においてて適切に与えるることにより，本解析モデデある．流下速速度の解析結果
果は概ね実験結果の傾向をを表せ
ルでで傾斜壁面粒子子サーマルの流流動を予測可能能であることているもののの，流下速度をを小さく評価している部分分も確が示示された．
Case7の流動初認できる．最最大厚さの解解析結果は，C初期の区間を除き，すべてのCasseで概ねよい一一致となってている．
(2) 粒子サーマル
ル内部の流動機機構，形状の考察
Case7の実験結結果においてて，流動初期のの厚さはCase8，9に8，解析結果の00cmにおける
図-6はCase5，8
図の流下距離10る比べ顕著に小小さい値となっっていることから，流入ゲゲート
水路
路横断方向中央央での粒子濃度度コンターと流流速ベクトル
ル
開放時のゲーート内のバリウウム粒子混合合水の厚さに問問題が図でである．サーマ
マルフロント部部の上部を中心とした循環環あったのでははないかと考えられる．ままたCase7～9ににおけ流がが形成されておおり，後方で濃濃度塊が取り残残されている
る
る浮遊粒子総総量の一致は非非常に良い．
Caseの方が周囲のががわかる．またた勾配が急なC
囲水をより多多以上の結果
果より，傾斜壁
壁面粒子サーマルの実験にに対し，く連連行し形状が発発達しており，一般的な傾斜斜壁面サーマ
マルの特徴（秋山ら15））をよく表している．
図-7はCase5，8，解析結果の流下距離100cmにおけるサーマル粒子濃度コンターの三次元形状図である．サーマルの表面でc = 0.01となるように表示している．サーマ

of Fluid Mech., Vol.8, pp.275-310, 1976.
Parker, G.: Conditions for the ignition of catastro-phicallyerosive turbidity currents, Marine Geology, Vol.46,
pp.306-327, 1982.
3) Parker, G., M. Garcia, Y. Fukushima and W. Yu: Experiments on turbidity currents over an erodible bed, J. of Hyルの先端部は丸みを帯びた形状をしており，上部や後方
draulic Research, Vol.25, No.1, pp.123-147, 1987.
では複雑な形状となっているのが確認できる．また勾配4) 福嶋祐介，金子幸弘: 固体粒子を浮遊する傾斜サーマの違いにより後方に巻き上げられる粒子の挙動や，フロル減速域，水工学論文集，第 41 巻，pp.537-542, 1997.
ント部の発達状況が大きく異なっていることがわかる．
5) 福嶋祐介，萩原達司: 浮遊傾斜サーマルの減速域の流動特性，水工学論文集，第 42 巻，pp.517-522, 1998.
6) 福嶋祐介，萩原達司，坂本充男: 粒子浮遊サーマルの5. 結論
実験による煙型雪崩流動モデルの適用性の検討，雪氷，第 60 巻 6 号，pp.453-462, 1998.
泥水流の流動形態の一つである傾斜壁面粒子サーマル
7) 福嶋祐介，萩原達司，坂本充男: 固体粒子浮遊による傾斜サーマルの流動特性，水工学論文集，第 43 巻，を対象として，オープンソースの解析ツールである
pp.929-934, 1999.
OpenFOAMを基に解析を行った．粒子の沈降，堆積を考8) 衞藤俊彦，福嶋祐介: 乱流モデルによる非定常泥水流
慮できるように基礎方程式，境界条件に変更を加え，粒の流動解析，海岸工学論文集，第 49 巻，pp.446-450,子サーマルの解析モデルへと発展させた．この解析モデ
2002.
9) 衞藤俊彦，福嶋祐介: 加速を伴う泥水サーマルの流動
ルに，LESモデルとしてDynamic SGSモデルを乱流計算解析，水工学論文集，第 47 巻，pp.1171-1176, 2003.に用いて三次元流動解析を行い，傾斜壁面粒子サーマル
10) OpenFOAM – The Open Source Computational Fluidの既存の実験結果と比較した．その結果，解析条件におDynamics(CFD) Toolbox, http://www.openfoam.com, 参
いて粒子の代表粒径を適切に与えることにより傾斜壁面照 2015-03.
11) Germano, M., Piomelli, U., Moin, P. and Cabot, W. H.: A
粒子サーマルの流動を予測可能であることが示された．
dynamic subgrid-scale eddy viscosity model, Phys. Fluid,
また解析結果より，粒子サーマル内部の流動機構や三次Vol.A3, pp.1760-1765, 1991.
元形状図を描き，サーマル内部の流動や形状の発達過程
12) Lilly, D. K.: A proposed modification of the Germanoを示すことができた．
subgrid-scale closure method, Phys. Fluid, Vol.A4,
pp.633-635, 1992.
自然界で発生するような泥水流の流動予測ができる解
13) Issa, R. I.: Solution of the implicitly discretized fluid flow
析モデルへ発展させていくための今後の課題として，本
equations by operator-splitting, J. Comp. Phys, Vol.62,
解析モデルに走路上からの粒子の巻き上げの計算を導入
pp.40-65, 1986.
することが挙げられる．またサーマル内部の流動や三次
14) 重枝未玲，秋山壽一郎，杉山誉: LES 乱流モデルを用いた水平流動する単一粒径・二粒径混合粒子サーマ
元形状についての定量的な検証が必要であるとともに，ルの数値解析，水工学論文集，第 52 巻，pp.1249複雑な地形のような，より実現象に近い条件を対象に，1254, 2008.
本解析モデルの妥当性を検証していくことが必要である． 15) 秋山壽一郎，浦勝，早川史郎: 保存性 2 次元傾斜サー
マルと傾斜プルーム先端部の流動特性の比較，水工学論文集，第 40 巻，pp.543-548, 1996.
参考文献
1)

2)

Inman, D.L., C.E. Nordstrom and R.E. Flick: Currents insubmarine canyons: an air-sea-land interaction, Annual Rev.

(2016. 3. 16 受付)

NUMERICAL ANALYSIS OF INCLINED WALL SUSPENSION THERMALSUSING Dynamic SGS MODEL IN OpenFOAMToshihiko ETO, Takuro KIKUCHI and Norihito OHSAWA
The inclined wall suspension thermals are analyzed using OpenFOAM, which is one of the threedimensional computational fluid daynamics tools. Dynamic SGS model is used for analyzing. The analysis model consists Finite Volume Method (FVM), TVD-superbee scheme and treatment of the particle
deposition. The variations of the travel speed, the maximum height and amount of total suspended particles of a thermal are obtained from numerical results and compared with the experimental results. Theagreement of numerical results and experimental ones is fairly well and it is found that this model can reproduce behavior of suspension thermals on inclined wall.