2. 固液混相乱流モデル
(1)固液混合系流れ
固液混合系流れは,計算格子(グリッドスケール:GS)
で平滑化した連続式およびNavier-Stokes式によって記述 される.粒子と流体が共存する計算格子では,計算格子 における各相の占有率で重み付けした物性値を用いる.
また,粒子内部での固液混合系の粘性係数およびせん断 歪みはゼロの取扱いとする.以下に基礎方程式を示すが,
乱流モデルとして,0方程式モデルのLES(Smagorinsky モデル)を採用した.
………(1)
………(2)
………(3)
………(4)
………(5)
………(6)
ここに,固液混合系流れに対して,ui:xi方向の流体速 度,ρ:密度,p:圧力,kSGS:サブグリッドスケール
(SGS)乱れエネルギー,gi:重力加速度,µ:粘性係数,
νt:SGS乱流粘性係数,Dij:歪み速度テンソル,Cs:
Smagorinsky定数,∆:計算格子幅,φp:各計算格子に含
まれる粒子の体積占有率である.なお,本研究では,直 交座標系xi(i=1,2)を考え,式中の下付き添え字p,lは それぞれ,粒子および流体が占有する領域における物性
シートフロー漂砂における鉛直分級過程の高解像度計算
Numerical Simulation for Vertical Sorting in Sheetflow by Particle-Laden Turbulence Flow Model
原田英治
1・後藤仁志
2Eiji HARADA and Hitoshi GOTOH
Sheetflow is the key to predict a beach deformation because massive sediment is transported by a high shear stress. To understand the mechanism of sheetflow dynamics, numerical simulations and experiments have been performed, however, sheetflow is a highly-concentrated particulate flow internal structure of which, so it is difficult to be investigated experimentally. In the present study, numerical simulation of a sheetflow by using of high resolution model the solid-liquid two-phase turbulence flow is performed, and the internal structure of sheetflow layer has been investigated in detail by a high-resolution solution in 1/10 times of particle diameter.
1. はじめに
シートフロー漂砂は多量の土砂輸送を伴うため,海浜 変形を論じる上で予測精度の鍵となる現象と言える.こ れまで,水理実験および数値シミュレーションからシー トフロー漂砂機構の検討が進められ,一定の成果が挙げ られてはいるが,シートフロー漂砂が高濃度混相流であ るため,その内部構造の理解は必ずしも十分ではない.
混相流計測技術には未だ制約が多いことを考えると,内 部構造の詳細な理解には,数値シミュレーションによる 検討が有効であることに疑念は無い.数値シミュレーシ ョンによって,シートフロー漂砂の内部構造を詳細に検 討するには,i)固相-固相,ii)固相-液相,iii)固液混 相乱流の物理素過程に対する適切な評価が要求される が,現在の研究水準では,これらi),ii),iii)の要件を 十分に満足するモデルは発展途上にある.また,計算機 性能の制約から,特に,粒子周りの流れ場が粗視化され て取り扱われる場合が多く,粒子スケール以下の乱れ構 造を含めてシートフロー漂砂を検討した研究は見当たら ない.本研究では,粒子スケールの1/10程度の細かい流 体計算格子を用いる高解像度計算によって,シートフロ ー漂砂における鉛直分級過程を計算力学的観点より検討 し,シートフロー漂砂の内部構造に詳細にアプローチす る.また,シートフロー漂砂が高粒子レイノルズ数の現 象であることから,流体の高解像度計算にはサブグリッ ドスケールの乱れのモデル化が不可欠となる.ここでは,
原田ら(2008)の高解像度型固液混相流計算手法の枠組 みに,LES(Smagorinskyモデル)(Smagorinsky, 1963)
を導入し,LES導入の効果がシートフロー漂砂における 鉛直分級過程に与える影響について,詳細に検討した.
1 正会員 博(工) 京都大学准教授 工学研究科都市環境工学 専攻
2 正会員 博(工) 京都大学教授 工学研究科都市環境工学専攻
値や物理量を示す.また,Smagorinsky定数Csの与え方 には議論の余地があるが,ここでは,Cs=0.1の一定値を 与えた.
(2)粒子モデル
粒子間接触力を評価しつつ多数粒子運動を追跡するた め,個々の粒子運動には個別要素法を適用する.並進及 び回転の運動方程式は,流体力および粒子間相互作用力 に起因する駆動力を用いて,以下の様に記述される.
………(7)
………(8)
………(9)
ここに,mp:粒子質量,vpi:xi方向の粒子の移動速度,
fpi:粒子を含む計算格子における粒子に作用する流体力,
Fflow i:粒子に作用する流体力,Fpint i:粒子間相互作用力,
Ip:粒子の慣性テンソル,ωp:粒子の角速度,Tflow i:流 体力に起因するトルク,Tpint i:粒子間力に起因するトル クである.なお,Fflow iは粒子を含む計算格子におけるfpi
の総和によって算定し,Tflow iは粒子の境界を含む計算格 子に対するfpiとブロック構成要素の中心座標から計算格 子中心座標までの相対位置ベクトルr(式(10)参照)
からモーメントを算定し評価した.
要素間相互作用力ベクトルFpint iは,接触粒子間の法線
(n)および接線(s)方向に配置された弾性スプリング および粘性ダッシュポットによって評価される(図-1参 照).また,非粘着性材料を対象として,法線(n)方向 には引っ張りに抵抗しないジョイントを,接線(s)方 向には一定の力が作用すると滑動するジョイントをそれ ぞれ配置した(後藤ら,1995).なお,弾性スプリング
(k)および粘性ダッシュポット(c)は,石材を対象と した物性値を用いて,ヘルツの弾性接触理論を準用して 設定した.なおこのモデル定数群を用いた場合,均一粒 径を対象とした数値移動床に一方向流を作用させたシミ ュレーションの結果,底面せん断力と流砂量の関係は,
乱流モデル導入の有無によらず活発な粒子移動を呈する 高Shields数の条件においては概ね一致することが確認さ れた.
(3)数値解法
先ずCIP-CUP法(Yabeら,1991)にならい固液混合系 流れを求める.次に,粒子要素に作用する流体力を式
(9)を用いて計算する.続いて,個別要素法による粒子 運動計算を行う.なお,個別要素法の安定した計算には,
計算時間刻みを小さくする必要がある.そのため,固液 混合系流れの計算刻み1ステップに対して,これと同等 の計算刻みとするために複数回数の個別要素法の計算ル ープを実施することになる.この複数回の個別要素法の 計算では,流体力fpiは計算効率を考慮して一定値を与え ている.最後に,粒子の移動計算後に得られた位置およ び速度ベクトルに応じて,各計算格子での粒子占有率を 修正し,粒子を含むセルにおける流れ場については,
………(10)
によって修正する.そして,この境界条件の下で上記に 示した一連のプロセスを,所定の計算ステップまで繰り 返す.
図-1 スプリング-ダッシュポットシステム
図-2 濃度重心の時系列
3. シートフロー漂砂の鉛直分級過程
(1)計算条件
水中に配置した3粒径混合状態の数値移動床を対象計 算領域とする(図-3のinitial condition参照).粒子の比重
は2.48であり,表-1に示す粒子をランダムに配置した.
重力の作用下で粒子のパッキング完了後,周期T=6.0s,
平均粒径に対するShields数Ψ=0.4(Komarら,1974)の 振動流を,計算領域の水平(x軸)方向に圧力勾配を付加 することで作用させ,対象とする水理条件がシートフロ ー漂砂の発現領域であることは別途確認している.また,
本研究で実施した計算条件における掃流粒子の平均ステ ップレングスが概ね計算領域幅程度となるようにx軸方 向の計算領域幅を決定した.また,流体計算で用いる正 方計算格子幅∆と基準粒径dとの関係は,∆/d=1/8とした.
上記同条件をSmagorinskyモデルを導入していないケ ース(以後,「LES無し」)とSmagorinskyモデルを導入し たケース(以後,「LES有り」)について,シートフロー 条件での混合粒径粒子の鉛直分級過程のシミュレーショ ンを実施した.
(2)濃度重心
各粒径の濃度重心の時系列を図-2に示す.ただし,パ ッキング完了後の平均移動床面から大粒子径の2倍以上 の高度に抜け出した粒子は浮遊に遷移したと見なし,除 外して算定した.LESの有無によらず双方の結果に,大 粒子の濃度重心の時間上昇過程が示されており,鉛直分 級の発生が理解できる.また,双方のケースの同時刻 t/T=0.6付近において,小粒子の移動開始が示されている ことも,双方に類似した傾向として確認されるが,振動 流作用後しばらくは,乱流の影響が顕著ではない状態で あ る こ と が う か が え る . 大 粒 子 が 上 昇 し 始 め る 時 刻 t/T=1.5についても双方に大きな違いはないが,その後の 大・中粒子の挙動に関して,「LES有り」は,大きな変動 を伴い「LES無し」と顕著な違いを示すことから,乱れ による流れ構造の変化が分級過程に影響を与えているこ とが推察される.
(3)瞬間画像
図-3に「LES有り」の鉛直分級過程の代表的なスナッ プショットを流速ベクトルと併せて示す.移動床表層付 近(-2.5<y/d<5.0)には,振動流の作用によって複雑な渦 図-3 鉛直分級過程のスナップショットおよび流速ベクトル分布「LES有り」
構造の発生が確認される.また,同時に粒子運動が活発 化し,鉛直分級が発達する様子が見て取れる.詳細に見 ると,時刻t/T=0.75では,0.0<x/d<5.0の移動床表層に示 されるような大規模渦によって表層を被覆していた小・
中粒子が水流中に巻き上げられている様子が示されてい る.また,この様な表層粒子の運動は近接する移動床表 層直下の粒子運動の活発化にも影響を及ぼしていること が,粒子配列の変化から理解される.さらに,この粒子 配列の変化によって生じた粒子間の隙間に流体が入り込 み,流体と粒子の運動量交換の促進によって交換層内の 粒子運動が発達するものと考えられる.時刻t/T=1.25で は,初期粒子配置では,y/d<0.0に埋没していた大粒子3 および5が移動床表層へ上昇しており,鉛直分級の発達 が明瞭に確認される.また,y/d>7.5には小粒子の存在が 確認されるが,シミュレーション結果を時間を追って確 認すると,これらの小粒子は計算終了まで概ねy/d>7.5の 速い流速領域中を移動していた(図-4の数密度分布を参 照).すなわち,この時刻では流体運動に追随し易い一 部の小粒子が掃流から浮遊へと遷移し,掃流と浮遊の粒 子 移 動 形 態 が 混 在 す る 状 況 に あ る と 言 え る . 時 刻
t/T=1.75以後にも,移動床表層付近に発達する渦によっ
て攪拌される粒子運動が継続して見て取れ,交換層内で の活発な流体と粒子の運動量交換がうかがえる.
(4)粒子数密度分布
粒子数密度分布の発達過程を図-4に示す.LES導入の 効果を明確にするために,「LES有り」および「LES無し」
の計算結果を示す.LESの有無に因らずy/d=0.0の移動床 表層付近における大粒子数密度の増加が認められる.ま た,LESの導入によってシートフロー層上部に存在する 浮遊粒子の存在が時刻t/T=1.25以降に顕在化することが 見て取れ,LES導入による乱れの効果が理解できる.
(5)乱れ構造
交換層付近での乱れ構造を検討するために,グリッ ドスケールの歪み速度テンソルの大きさのスペクトル を (a) 貯 留 層 (- 5 . 0 <y / d≤- 2 . 5),(b) 交 換 層 下 部
(-2.5<y/d≤0.0),(c)交換層上部(0.0<y/d≤2.5)そして
(d)低濃度粒子層(2.5<y/d≤5.0)に分けて図-5に示す.
なお,t/T=1.0以降の各層での平均歪み速度テンソルの 大きさの時系列を用いてスペクトル解析を実施した.貯 留層では,粒子が密に存在するため,流体の存在率が低 図-4 粒子数密度分布の時間発達過程
構の検討を続けたい.また,この種の数値シミュレーシ
ョンの3次元計算には大きな計算負荷が伴う.そのため,
対象とする問題に対して要求される計算精度に応じて,
2次元あるいは3次元計算の選択は必要であると考える.
その意味でも,2次元計算で得られる解の再現性につい ての検討は重要であり,今後,2次元計算結果と実験結 果との定量比較を実施するとともに,Smagorinsky定数 Csの値の与え方についても検討していきたい.
謝辞:本研究に助成頂いた財団法人中部電力基礎技術研 究所に感謝申し上げる.また,本研究の遂行にあたって,
京都大学大学院修士課程 鶴田修己君の熱心な協力を得た ことを記して,謝意を表する.
参 考 文 献
後藤仁志・酒井哲郎(1995):表層せん断を受ける砂層の動的 挙動の数値解析,土木学会論文集,No. 521/II-32,pp.
101-112.
原田英治・後藤仁志(2008):高解像度固液混相流モデルを用 いた水中投入ブロック群沈降・堆積過程の数値シミュレ ーション,海岸工学論文集,第55巻,pp. 961-965.
Komar, P.D. and Miller, M.C.(1974): The initiation of oscillatory ripple marks and the development of plane-bed at high stresses under waves, Jour. Sedimentary Petrology, Vol.45, No.3, pp.697-703.
Smagorinsky, J.(1963): General circulation experiments with the primitive equations, Mon. Weath. Rev., Vol.91, No.3, pp.99- 164.
Yabe,T. and Wang, P.Y.(1991): Unified Numerical Procedure for Compressible and Incompressible Fluid,J.Phys.Soc.Japan,Vol.
60,No. 7, pp. 2105-2108.
いことから,乱れエネルギー生成に起因する歪み速度分 布が他の層と比較して小さいレベルにあると考えられ る.一方,交換層では,流体と粒子が混在し,活発な運 動の交換が促進されているため,貯留層と比較して,低 周波から高周波にわたり乱れエネルギーの生成も活発で あると考えられる.また,交換層より上方の低濃度粒子 層では,交換層のスペクトルよりも抑えられた形状が見 られる.これは,粒子混在率が交換層と比較して低下し ているため,流体と粒子界面での速度歪みが発生し難い ことが原因であると考えられる.以上のように,粒子混 入による乱れ場への影響は有意であり,適正な予測には 乱流モデル導入が必要である.
4. おわりに
混合粒径粒子のシートフロー漂砂における鉛直分級過 程を高解像度型の固液混相流モデルにLES(Smagorinsky モデル)を導入した固液混相乱流モデルを用いて,計算 力学的観点から詳細に検討した.LES導入によって粒子 運動が活発化することが,粒径別の濃度重心の時系列や,
粒子数密度分布の時間発達過程から確認された.また,
各時刻の流速場の瞬間像より,移動床表層での大規模渦 の発生とそれによる粒子運動の活発化が示され,鉛直分 級過程での交換層内部の流体・粒子相互作用構造の一端 が明らかになった.さらに,移動床表層付近の乱れ構造 についても歪み速度テンソルのスペクトル解析から,交 換層での粒子混入による乱れ場の促進を示し,乱流モデ ル導入の意義を示した.今後は,モデルの3次元化を喫 緊の課題として,計算力学的観点での詳細な鉛直分級機
図-5 歪み速度テンソルの構造