第 35 回数値流体力学シンポジウム
E10-1

APM 法による石礫型土石流に関する数値実験
Numerical experiments on stony-type debris flow by APM method○ 加藤宏季, 中央大, 東京都文京区春日 1-13-27, E-mail: a17.jnt5@g.chuo-u.ac.jp
福岡捷二, 中央大,同上, E-mail: sfuku@tamacc.chuo-u.ac.jp
Three-dimensional numerical experiments were conducted using CFD-DEM with a one-fluid model for the solid-liquidboundary. The results of the numerical experiments showed that bank erosion advanced at the period of the passage ofthe debris flow front. In addition, the advance speed of the debris flow in the front with a large concentration of debris
became small, while the velocity of the debris flow in the tail part with a small concentration became large. Furthermore,
the maximum bank erosion was just after the debris flow front had passed. These results are consistent with previous
findings on debris flows. The water pressure increased with the passage of the debris flow front, but the hydrodynamic
force acting on particles did not change.

１． 序論
石礫型土石流は，連続体としての運動をする水と離散的な運動をする石礫や流木などが一体となって流れ下る現象である．流下中の石礫粒子に作用する力や水の流れ場の様子を詳細に精度よく計測することは難しいため，数値解法による現象解明が期待される．本来であれば，石礫粒子は離散体モデル，水は連続体モデルで考えられ，これらの相互作用を考慮した形でモデル化されることが望ましい．しかし，既往研究 1)では土石流を水と土砂の混合流体とみなし，連続体の枠組みでの Euler-Euler 型のモデル化が主流であった．この考え方は，Euler-Lagrange 型のモデル化に比べて計算負荷が小さく，小さい粒径の土砂も考慮でき，現地スケールでの解析が可能であるため実用的であるものの，本質的な現象解明とは言い難い．
近年では，物体より小さな計算格子を用いた大規模な計算により，流れ場を Euler 的に詳細に解き，物体の運動を個別要素法に基づき Lagrange 的に解析する固液混相流場の解析法 2)3)が開発されている．Euler-Lagrange 型の解析で重要な固液境界の取り扱いにおいては，直接強制型埋め込み境界法や，一流体モデルがある．直接強制型埋め込み境界法 4)は，境界セルの流速を補完し，NavierStokes 方程式に boundary body force と呼ばれる外力項を与える手法である．各計算格子内に含まれる固相の体積割合と固体の向きを流れの解析に考慮できるという利点を持つ一方，本研究で対象とする土石流のように，固体と液体の境界形状が時空間的に変化する場においては，流速を統一的に補完することは容易ではない．
それに対して，一流体モデルは流れの運動に及ぼす固体の影響を統一的に考慮可能である．牛島ら 2)は，一流体モデルを用いた 3 次
元固気液多相場での解析モデル(3D MICS)を提案している．Fukuda
et al.3)は，固液境界の取り扱い法として一流体モデルを用いた任意形状の粒子群の解析を可能としたモデルを構築した．
福田ら 5)は，現地スケールへ適応するために各プロセスに複数の計算領域ブロックを割り当てる動的負荷分散手法を構築し，2014 年 8 月に広島の安佐南区八木で発生した土石流を対象とした数値解析を行い，土砂や流木の流動機構を説明している．
本研究では，福田らの APM 法を用いて数値実験を行い，石礫型土石流流下時における圧力場と粒子に働く流体力との関係を考察した．
２． 数値計算法
本研究では，流れを Euler 的に，石礫を Lagrange 的に解析する
3)．流れの解析においては，固液混相流の一流体モデルを用いて粒子の固相部分を密度の異なる液体として考え，非圧縮性流れとしCopyright © 2021 by JSFM1て解析する．普遍性が期待できる格子サイズより小さい流れ場の乱れについては，Smagorinsky モデルによりモデル化し，LES として解いている．石礫粒子の運動については，接触力は DEM によりモデル化し，流体力は流れ場の運動方程式の圧力項と拡散項を体積積分することによって，抗力係数や揚力係数を用いずに評価する．
(1) 流れ場の解析
流れの基礎方程式を以下に示す．
𝜕𝑢𝑖(1)= 0 (𝑖 = 1,2,3)
𝜕𝑥𝑖
𝜕𝑢𝑖
𝜕𝑢𝑖
1 𝜕𝑃
𝜕
(2)+ 𝑢𝑗= 𝑔𝑖 −+{2(𝜈 + 𝜈𝑡 )𝑆𝑖𝑗 }𝜕𝑡𝜕𝑥𝑗
𝜌 𝜕𝑥𝑖 𝜕𝑥𝑗
1 𝜕𝑢𝑖 𝜕𝑢𝑗(3)𝑆𝑖𝑗 = (+)
2 𝜕𝑥𝑗 𝜕𝑥𝑖𝜇(4)𝜈=𝜌𝜈𝑡 = (𝐶𝑠 Δ)2 √2𝑆𝑖𝑗 𝑆𝑖𝑗(5)

ここに，𝑢𝑖 は流体計算格子内で質量平均した流速，𝑃は圧力と SGS応力の等方成分の和，𝜌は流体計算格子内で体積平均した密度，𝑔𝑖は重力加速度，𝜇は流体計算格子内で体積平均した粘性係数，𝐶𝑠 はSmagorinsky 定数，∆は計算格子幅である．下付きの添え字𝑖, 𝑗はそ
れぞれ各方向成分を示す．
流れの解析における固体境界の取り扱いについては，固液混相の一流体モデルを用いる．すなわち，各計算格子内に含まれる各相同士の間には界面が存在し，本来であれば各相が混ざり合うことはないが，本研究では複数の相が混ざり合った混合流体として扱う．この時，物性値𝜙(密度，粘性係数)および流速は，それぞれ次のように体積平均，質量平均とした．
(6)
𝜙 = 𝛼𝜙𝑠 + (1 − 𝛼)𝜙𝑓 , 𝜙𝑓 = 𝑓𝜙𝑙 + (1 − 𝑓)𝜙𝑔
(1𝛼ρ𝑠 𝑢𝑠𝑖 + − 𝛼)𝜌𝑓 𝑢𝑓𝑖
(7)
𝑢𝑖 =
𝜌
ここで，𝛼：流体計算格子内の気相を除く領域での固相の体積割合，𝑓：流体計算格子内の液相(粒子も含む)の体積割合，下付きの添え字𝑠, 𝑓, 𝑔はそれぞれ固相，液相，気相を示し，下付きの添え字𝑙は粒子も流体とみなす一流体モデル上の液相を示す．𝛼と𝑓の関係をFig.1 に示す．連続式，運動方程式はスタガード格子を用い SMAC法を用いて解き，圧力の収束計算には SOR 法を用いた．
固相の体積割合𝛼の算定精度は，混合流体の流速𝑢𝑖 の評価に大きく影響を及ぼす．本研究では，固相の体積割合𝛼はサブセル法 2)を用いて評価する．サブセル法では，Fig.2 に示すように，計算格子

第 35 回数値流体力学シンポジウム
E10-1particle

ΔcellFig.1 relationship between 𝛼 and 𝑓よりも小さい体積割合算定用のサブセルを用いて体積割合𝛼を評価する．サブセルの中心が固体表面よりも内側に位置するときにそのサブセルを固相と判定する．本研究では，サブセルは一様な大きさの立方体とし，その大きさはΔ = 2Δ𝑠𝑢𝑏 とした．
また，自由水面の変化は VOF 法により解析している．
(2) 粒子運動の解析
基礎方程式を以下に示す．
(6)
𝑀𝑟̈𝐺 = 𝑀𝑔 + 𝐹𝑓 + 𝐹𝑐𝜔̇ 𝑟 = 𝐼𝑟−1 {𝑅−1 (𝑁𝑓 + 𝑁𝑐 ) − 𝜔𝑟 × 𝐼𝑟 𝜔𝑟 }𝑛𝐹𝑓,𝑖= ∫ {−𝛺𝑠𝜕𝑝𝑛
𝜕𝑛+ 𝜌𝑛{2(𝜈 𝑛 + 𝜈𝑡𝑛 )𝑆𝑖𝑗}} 𝑑𝛺𝜕𝑥𝑖𝜕𝑥𝑗𝑛𝑁𝑓,𝑖= ∫ 𝜀𝑖𝑗𝑘 𝑟𝑓,𝑗 {−𝛺𝑠𝜕𝑝𝑛𝜕𝑛 }}{2(𝜈 𝑛 + 𝜈𝑡𝑛 )𝑆𝑘𝑙+ 𝜌𝑛𝑑𝛺𝜕𝑥𝑘𝜕𝑥𝑙

(7)
(8)

solid phase
ΔΔsubΔsubsub cell

Fig.2 sub cell method(9)

ここで，𝑀：質量，𝑔：重力加速度，𝑟𝐺 ：重心位置，𝐹：粒子が受ける力，𝑁：粒子が受けるトルク，𝜔：角速度，𝐼：慣性テンソル，
𝑅：ローカル座標系からグローバル座標系への変換行列，𝜀𝑖𝑗𝑘 ：エディトンのイプシロン，𝑟：作用点の重心からの距離である．下付
き添え字𝑓は流体の成分を，c は固体の成分を，𝑟は石礫固定座標系をあらわす．
３． 数値実験条件
本研究で用いた数値移動床水路は,水路長 60 m，水路幅 8 m の直線水路である．水路勾配は，土石流の発生条件 1)を満たすように，20°とした．座標軸は，流下方向に x 軸，横断方向に y 軸，z 軸方向のことを垂直方向と呼ぶ)をとる．
垂直上向きに z 軸(以下，渓床・渓岸の材料として，Fig.3 に示す最大粒径が 1.0m，最小粒径が 0.1m の，10 種類の粒径からなる粒度分布を用いた．すべての粒子は形状と密度が等しい．本研究では，非球形粒子を用いており，粒径とは非球形粒子の体積と同等の体積を持つ球の直径のことをいう．
水路中に粒子群をランダムに投下し，渓岸の傾斜角は渓床から約 30°となるように初期堆積縦横断河道をつくった．数値実験水路の初期状態を Fig.4 に示す．
実験では，渓谷の上流端で，一定流量 10.0[m3 ⁄s]を与え続けた．
通水中に土砂の供給は行っていない．下流端では，圧力をゲージ圧 0 とし，土砂は自由に流出できるとしている．土石流の発生を促すために初期条件として，渓床表層より上層部分の流速に対数分布則を与えた．
本解析で用いたパラメータを Table.1 に示す．福田らの検討 6)を参考に，パラメータを決定した．
Copyright © 2021 by JSFM2

Table-1 Parameters used in simulation∆𝑥, ∆𝑦, ∆𝑧: Fluid-computational
cell size∆𝑡: Time step for fluidsimulation∆𝑡′: Time step for particle
simulation
𝜌𝑤 : Density of water𝜌𝑠 : Density of particles
𝜇𝑤 : Viscosity coefficient ofwater𝐶𝑠 : Smagorinsky constant
𝐸 : Elastic modulus
𝑝𝑜𝑠: Poisson’s ratio
𝜙𝑝 : Internal friction angle
𝑏: repulsion factor
Number of fluid-computational
cells
Number of particles

0.025

m

1.0 × 10−4s

2.0 × 10−6s

1,000kg⁄m3

2,650kg⁄m3

8.9 × 10−4Pa ∙s0.173
5.0 × 1010

Pa0.3326.6
0.7
153,600,000
66,619°

第 35 回数値流体力学シンポジウム
E10-1

8.0 m

valley bankpercentage [%]60.0 mvalley bed

diameter [m]Fig.3 particle size distribution

Fig.4 initial state

front
015x [m](a) before passing of the debris flow(a) before passing of
the debris flow

60
front
015x [m](b) during passing of the debris flow front(b) during passing of thedebris flow front

60
0front

15x [m]60(c) after passing of the debris flow front(c) after passing of thedebris flow front
Fig.5 cross sections

Fig.6 schematic diagram of flow conditions

Copyright © 2021 by JSFM3

第 35 回数値流体力学シンポジウム
E10-1
z [m]

z [m]
velocity[m⁄s]

volume fraction[−]

frontfront
x [m]

x [m]

Fig.7 contour of flow velocity

Fig.8 contour of volume fraction

waterpressure[N⁄m2 ]
z [m]

x [m]
(a) initial statez [m]debris flow frontx [m](b) during passing of the debris flow front
z [m]

x [m]
(c) after passing the debris flow front
Fig.9 contour of water pressure
Copyright © 2021 by JSFM4４．数値実験結果
Fig.5 には，上流端から 20m 地点での土石流フロント通過前，通過中，通過後の横断面形を示す．Fig.6 はそのときの流況の模式図を示す．Fig.5 より，渓岸侵食は土石流フロント通過中よりは，土石流フロント通過後付近の時間帯に大きく起こっていることがわ
かる．Fig7，Fig.8 には，𝑡 = 20秒のときの水路中央断面での流速コンター，粒子の体積割合をそれぞれ示す．Fig.7，Fig.8 より，土石流フロントでは石礫粒子の体積割合が高く，巨礫が集中している．さらに，流速は小さくなっている．一方，後続部では体積割合が小さく，流れている粒子の大きさは相対的に小さい．そして，流速は大きくなっている．これらの傾向は，Lyu7)らの実験結果や，石礫型土石流の特徴 1)と一致している．
Fig.9 に，土石流フロント通過前，通過中，通過後における水路中央断面での圧力コンター図を示す．土石流フロント通過により，増加している様子がわかる．土石流による侵食現象の要因として，
現地 8)や模型実験 9)10)での計測結果から，堆積層の間隙水圧が急上昇することによって液状化のような状態となることが言われている．本研究で用いた計算手法(APM 法)では，個別の粒子に働く流体力が求められる．そこで，間隙水圧の上昇と実際に粒子に働く流体力との関係性について検討する．また，侵食メカニズムの解明に向けた検討であるため，水路床付近の動かない粒子ではなく，侵食を受けた層の粒子群を対象に検討を行う．Fig.10 には，Fig.11で示す粒子群(x=16~22m かつ z=2.2~2.8m に位置する渓床粒子)に働く無次元力を，各時間ごとに平均化した値の時間変化を示す．
本研究では，垂直方向に働く力𝐹𝑧 を粒子の水中重量の z 方向成分((𝜌𝑠 − 𝜌𝑤 )Vgcos 𝜃 )で無次元化を行う．𝜃は水路勾配である． 土石流による侵食現象の要因として，現地 8)や模型実験 9)10)での計測結果から，堆積層の間隙水圧が急上昇することによって液状化のような状態となることが言われているが，Fig.10 より，侵食可能層の粒子が受ける流体力は時間的に大きく変化していないことがわかる．
Johnson et al.10)の実験では，66%の砂(粒径：0.0625-2mm)，33%の礫(粒径：2-32mm)，1%の泥(粒径：0.0625mm 以下)で構成される粒子群を用いている．本研究で用いた粒子群の方が粗いため，本研究では空隙が大きく間隙水が流れやすい場であると考えられる．すなわち，土石流フロントが到達することで非排水条件での圧縮荷重が働くことによる間隙水圧の上昇というのは見られず，流体力が大きく変化するほどの圧力上昇とはならなかったと考えられる．

第 35 回数値流体力学シンポジウム
E10-1

contact forcehydrodynamic force

weightdimensionless force [-]total force

time [s]Fig.11 the averaged dimensionless force of a particle

Fig.10 featured particles

５．まとめ
本研究では，APM 法を用いて石礫型土石流を対象とした数値実験を行った．計算結果は，既往の研究と同様の傾向を示した．また，土石流フロントが到達することにより圧力は大きくなったが，粒子に働く流体力は大きくならなかった．これは，用いた粒子群が粗く，流体力に影響するほどの圧力変化が起こらなかったためであると考えられる．
また，本研究で使用したモデルでは渓岸粒子に働くサクションの効果を考慮していない．実現象を考えると地下水面より高い位置の不飽和土にはサクションが働くことにより安息角以上の傾斜角になったとしても崩壊せずに安定している．このサクションの効果を取り込むことを今後の課題とする．

Switzerland”, Journal of Geophysical Research, 116(2011), F01002.
(9) Iverson, R. M., Logan, M., LaHusen, R. G., and Berti, M.,” Theperfect debris flow? Aggregated results from 28 large ‐ scaleexperiments”, Journal of Geophysical Research, 115(2010), F03005.
(10) Johnson, C. G., Kokelaar, B. P., Iverson, R.M., Logan, M. LaHusen,
R. G., and Gray, J. M. N. T.,” Grain-size segregation and levee
formation in geophysical mass flows”, Journal of GeophysicalResearch, 117(2012), F01032.

参考文献(1) 高橋保，"土石流の発生と流動に関する研究"，京都防災研究所年報，第 20 号 B-2，pp.405-435，1977.
(2) 牛島省，山田修三，藤岡奨，禰津家久，"3 次元自由水面流れによる物体輸送の数値解法（3D MICS）の提案と適用性の検討"，土木学会論文集 B, Vol.62, 1 (2006), pp.100-110.
(3) Fukuoka, S. and Fukuda, T., “Effects of sizes and shapes of gravel
particles on sediment transports and bed variations in a numericalmovable-bed channel”, Adv. Water Resour., 72 (2014), pp.84-96.
(4) Fadlun, A., Verzicco, R., Orlandi, P. and Mohd-Yusof,J.,” Combined
immersed-boundary finite-difference methods for three-dimensional
complex flow simulations”, Journal of Computational Physics,Vol.10(2000), pp.35-60.
(5) Fukuda, T. and Fukuoka, S., “Interface-resolved large eddysimulation of a field debris flow with large stones and woods”, Adv.
Water Resour., Vol.129 (2017), pp.297-310.
(6) 福田朝生，福岡捷二，"土石流水路実験結果を用いた数値移動
床水路による高濃度粒子群と水流の力学的相互作用の検討"，
土木学会論文集 B(水工学), Vol.72, 1 (2016), pp.859-854.
(7) Lyu, L., Wang, Z., Cui, P. and Xu, M.,” The role of bank erosion on
the initiation and motion of gully debris flows”, Geomorphology.,
Vol.285(2017), pp.137-151(8) Berger, C., McArdell, B. W. and Schlunegger, F.,” Direct
measurement of channel erosion by debris flows, Illgraben,
Copyright © 2021 by JSFM5