第4章
雪崩シミュレーションを用いた被害範囲推定

那須雪崩の流下状況推定を目的とした流動シミュレーション

Simulations for estimating flow behaviors of Nasu Avalanche
伊藤陽一 1)，森口 周二 2)，小田 憲一 3)，西村 浩一 4)
Y. Ito1), S. Moriguchi2), K. Oda3), K. Nishimura4)
1)

1)

防災科学技術研究所 雪氷防災研究センター
Snow and Ice Research Center, National Research Institute for Earth Science and Disaster
Resilience2)
東北大学 災害科学国際研究所
2)International Research Institute of Disaster Science, Tohoku University3)
日本大学理工学部 土木工学科
3)
Department of Civil Engineering, College of Science and Technology, Nihon University4)
名古屋大学大学院 環境学研究科
4)
Graduate School of Environmental Studies, Nagoya University

Abstract: In this study, we reproduced Nasu Avalanche using two types of simulation tools to estimate flow
path and flow velocity. The flow path was firstly estimated by using TITAN2D, and the simulated flow area
was compared to the observed data reported in the field survey. Then, three-dimensional stabilized finite
element method was applied to the simulation to check detailed flow behaviors. According to the obtained
results, similar tendency of flow path was seen between results of two simulations, and maximum flow
velocity was estimated to be 11~12 m/s.
Key words: avalanche, flow area, flow velocity, computational fluid dynamics

1

はじめに雪崩発生後に行われた現地調査では，斜面下方の樹林帯で雪崩の流動に起因すると思われる積雪層構造の乱れなどが確認され，雪崩の流下範囲をある程度絞りこむことができた．しかし，斜面上方で雪崩が発生・流下したと思われる痕跡は視認できなかったため，発生箇所は天狗岩付近の斜面であろうとの推測にとどまった．このような背景のもと，雪崩運動シミュレー
ションによって雪崩の流下経路や流下速度を推定する試みを行った．本報告では，2 つの連続体モデルを用い，それぞれの特徴を生かし，現地調査で得られた流下推定範囲をどの程度再現できるかをまとめる。
2
簡便な連続体モデル（TITAN2D）を用いた雪崩の流下状況の推定
2.1 連続体モデル（TITAN2D）
雪崩の流下状況（流下経路や速度など）を推定するために，まず初めに簡便な雪崩運動モデルを用いて雪崩発生箇所を推定する試みを行った．使用したモデルは，土石流や地すべりなどの乾燥粒状体の流れを対象に開発された TITAN2D1）で，インターネット上で無料公開されている最新版（v4.0.0）を使用した．TITAN2D は，粒状体の流れを非圧縮性の連続体の流れとみなし，さらに流れの厚さが幅・長さに対して十分小さく，鉛直方向の速度勾配を無視できると仮定（浅水流近似）することにより，3 次元の流れを準 3 次元な 2 次元流れに変換して計算を行い流れの厚さや速度の時空間分布を求めることが可能である．なお，雪崩は流下時に斜面上の積雪を内部に取り込み，質量が流下距離にともない増加する場合があることが知られているが 2），TITAN2D では質量変化は考慮せず一定としている．
TITAN2D で用いられている質量保存式・運動量保存式は以下のとおりである．斜面を流下する流れに対し，斜面と平行に x 座標と y 座標（XY 平面）を，斜面と鉛直に z 座標をおいた場合に，厚さ方向に平均化された質量保存式は式(1)で表すことができる．

76̅̅

0(1)

ここで， は流れの厚さ， は時間， ̅ および ̅ はそれぞれ x および y 方向の速度（厚さ方向の平均値）である．また，同様に厚さ方向に平均化された運動量保存式（x 方向）は式(2)で表される．
̅̅̅ ̅̅

0.51̅tansgn̅

(2)sinここで， および は重力加速度の x および z 成分を，は底面摩擦角，は内部摩擦角である．式(2)の右辺の第一項は駆動力としての重力項である．第二項は底面におけるクーロン摩擦抵抗の効果を表し， は底面の曲率半径で無限大に近づくほど底面が平坦で曲率の効果が無視できるように示されている．第三項は流れ内部の速度勾配にともなう内部摩擦の効果を表し，は流れが流下方向や横断方向に広がる方向（ ）か縮む方向（ ）かで決定される土は ̅ ⁄ の正負によって正なら 1 を，負なら−1 を返す関数である．
圧係数，sgn̅ ⁄実際の TITAN2D を用いた計算は，雪崩発生後に無人航空機による空撮で得られた 1 m グリッドの雪面標高データ上に楕円柱を設定し，この楕円柱が崩壊し連続体として斜面上を流下するものとして行われる．雪崩の発生箇所や始動体積は，雪面標高データの座標や楕円の径および厚さを任意に設定することで決定される．その他，雪崩の流動状態（流れの拡がりや到達距離）に関連する底面摩擦角，内部摩擦角も入力パラメータとして設定することで計算を行った．計算時間は斜面の規模や地形の複雑さによって大きく変化するが，本解析では長くても 1～2 分程度で計算が終了することが多かった．
2.2 TITAN2D を用いた流下経路の推定結果
雪崩の発生箇所は不明であるが，現場の様子から天狗岩直下の斜面付近と推定されていたため（図－1），初めにこの推定結果を参考に雪崩流下状況を試算した．このときの初期条件は以下のとおりである．
位置：図－1 の推定発生区とほぼ同様

始動積雪：長径 55 m，短径 48 m，厚さ 0.6 m の楕円柱（体積 1244 m3）

入力パラメータ：底面摩擦角 30°，内部摩擦角 25°
なお，始動積雪の厚さは，雪崩発生時の推定積雪深に加えて雪の取り込みで雪崩の質量が増加することを見込んでやや多めに設定している．
雪崩の初期位置と，流下後 20 秒経過した後の流下範囲を図－2 に示す．図－1 に比べると，本流（埋没地点のある沢）に加えて南側の沢に流れが集中し，現地調査で雪崩の痕跡が確認された北側の沢にはほとんど流入しないという計算結果が得られた．このことから，実際の雪崩発生箇所は現地調査で推定された位置よりもさらに北側（埋没地点から天狗岩に向かって右方向）であった可能性が考えられる．そこで，雪崩の初期体積や標高などは変更せずに，初期位置のみ北方向に 5 m 程度ずつ変化させながら計算を行い，現地調査で得られた流下推定範囲と一致するような雪崩発生箇所を求めることにした．何例か試行錯誤しながら計算を行った結果，北方向に約 15 m 初期位置をずらした場合に，現地調査と比較的よく一致する結果が得られた（図－3）．なお，このときの雪崩の平均速度は最大で 12 m/s 程度と計算された（図－4）．

77

北天狗岩発生区
（推定）

埋没地点

図－1 現地調査（2017 年 4 月）後の雪崩の流下経路の推定図
天狗岩

埋没地点

図－2 雪崩の初期位置（左）と始動から 20 秒後（右）の流れ厚さの計算結果．

78

図－3 雪崩の初期位置を約 15 m 北方向にずらした場合（左）の始動から 20 秒後（右）の流れ厚さの計算結果

図－4 雪崩の平均速度の計算結果（図－3 の計算時）
3 3 次元安定化有限要素法を用いた雪崩の流下状況の推定
3.1 解析手法と支配方程式
本解析では，非構造四面体要素を用いた SUPG/PSPG 法 3)に基づく 3 次元安定化有限要素法を用いて雪崩の流下挙動を表現する．雪崩は非圧縮性の Bingham 流体と仮定し，下記の運動方程式と連続式を解くことにより流動挙動が表現される．
g0(3)

(4)

ここで， は密度，uは速度， は圧力， は等価粘性係数である．等価粘性係数は，せん断ひずみ速度 の関数として表現され，クーロンの破壊規準を導入した Bingham 流体モデルの割線粘性係数（疑似的に Newton 流体としてみた場合の等価な粘性係数）として下記のように求められる．

79

(5)
ここで， は最小粘性係数， は内部摩擦角， は粘着力である．なお，せん断ひずみ速度 は，せん断ひずみ速度テンソル を用いて下記のようにスカラー値で表現する．
12(6)

自由表面の表現には VOF（Volume of Fluid）法 4)を用いており，下記の密度関数 の移流方程式を解くことで自由表面の位置の時間発展を解いている．
0(7)

本手法は，もともと土砂流動を対象として開発された手法 5)であるが，これまでの研究で雪崩に対する適用性も確認されている 6),7)．また，完全な 3 次元条件での計算となるため，計算コストが高いという欠点があるが，地形形状や植生，および防護工の影響なども含めて解析が可能であり，適用範囲が広いのが特徴である．なお，本手法においても，TAITAN2D と同様に斜面上の積雪を取り込む効果は考慮しておらず，質量変化は一定としている．
3.2 解析モデルと解析パラメータTITAN2D の結果を踏まえて，図 5 に示す 2 ケースについて解析を実施した．Case 1 は頂上部の大岩（天狗の鼻）の直下に発生区を設定した条件であり，Case 2 は Case 1 の発生区の位置から北側に発生区を 15m 移動させた条件である．なお，図中には，現地調査の結果 8)から推定された雪崩の到達範囲および，救出活動地点の位置も含めている．解析メッシュは，図 6 に示すように地形に沿って高さ 10m の範囲で空間を分割することで作成している．地表面に近いメッシュほど細かくなっており，地表面付近では四面体メッシュの 1 辺の長さが約 15cm の解像度でモデル化している．
地形形状については，TAITAN2D と同様に，雪崩発生後に無人航空機による空撮で得られた 1 m グリッドの雪面標高データから作成しているが，計測結果の生データは植生の影響を含んでおり，部分的に不自然な凹凸が存在した．この凹凸は，雪崩の流動挙動や計算安定性にも影響を与える可能性があったため，平滑化処理を施し，図 6 のように滑らかな地形形状に変換している．なお，この処理のために，ラプラシアンフィルターという処理方法を用いたが，元の地形特性を変化させないように，平滑化の度合いを試行錯誤によって決定している．
本解析に用いた解析パラメータは表－1 に示すとおりであり，密度については現地の調査報
告結果 8)，粘着力と最小粘性係数については既往の研究 7)を参考に決定した．また，内部摩擦角については，別途実施した解析結果に基づいて決定している．なお，計算にはスーパーコンピューターを用いて並列数約 2500 程度の規模で行っている．

図－5 解析モデル（左：Case 1，右：Case 2）80図－6 解析メッシュ

図－7 地形形状の平滑化処理（左：平滑化処理前，右：平滑化処理後）密度(kg/m3)200

表－1 解析パラメータ
内部摩擦角(deg) 粘着力(N/m2)
15.0
100

最小粘性係数(Pa/s)
1.00

3.3 解析結果
得られた解析結果に基づいて可視化した流動挙動を図－8 と図－9 に示す．図－8 が Case 1の結果，図－9 が Case 2 の結果である．Case 1 の結果では，流下後の雪崩の到達範囲が現地調査から推定された到達範囲から外れた南側の沢に流れる結果となっており，実現象とは異なる流動挙動になっていると考えられる．これに対して Case 2 では，Case 1 よりも実際の到達範囲に近い結果となっている．本解析の結果は，現地調査で推定されている被災箇所と完全には一致していないが，全体的な流下方向として Case 2 の結果が実現象に近いと考えられる．そのため，発生区は Case 2 で設定した位置付近に存在した可能性が高いと考えられる．
一方で，Case 2 の結果から最大流速分布を求めたところ， 最大で 11m/s 程度となる箇所が確認された．本研究で用いた手法では，底面（斜面表面）の境界条件が non-slip（底面で速度が 0）条件となっており，底面付近のメッシュサイズが 10cm 程度であるため，底面摩擦が実現象よりも過大評価されている可能性が高い．また，実際の雪崩が表層雪崩であることを考えても，やはり今回の解析結果は流動速度を過小評価している可能性がある．このため，実際の流速の最大値は 11m/s 以上であった可能性が高い．

81

図－8 解析結果（Case 1）

図－9 解析結果（Case 2）

82

4

まとめ
本報告では，2 つの連続体モデルを用いた雪崩発生個所の推定や流下状況の推定についてまとめた．TAITAN2D を用いた解析では，雪崩の発生個所を現地調査で推定された範囲よりも北側に 15m の位置とすることで流下予測範囲に到達する結果が得られた．また，3 次元安定化有限要素法を用いた詳細な解析においても同様の結果となることが分かった．流下速度はTAITAN2D の結果では最大 12m/s，3 次元安定化有限要素法を用いた結果では最大 11m/s 以上であると推定された．以上のことから，両解析によって得られる結果はおおむね一致しており，雪崩の流下範囲を推定する際は TAITAN2D を使用した大まかな解析を採用し，より詳細に流下状況を推定する際は 3 次元安定化有限要素法を用いた解析によって行う方法が雪崩のソフト対策に有効である可能性を示唆することができた．
謝辞
3 次元安定化有限要素法を用いた数値解析に関して，八戸工業大学の高瀬慎介講師より多くの助言を頂いた．記して感謝の意を表す．

引用文献
1) Pitman E.B., C.C. Nichita, A.K. Patra, A.C. Bauer, M.F. Sheridan, M. Bursik: Computing granularavalanches and landslides. Physics of Fluids, 15(12), 3638−3646, 2003.
2) Sovilla, B., F. Sommavilla, A. Tomaselli: Measurements of mass balance in dense snow avalanche
events. Annals of Glaciology, 32, 230−236, 2001.
3) Brooks,A.N.,Hughes,T.J.R.:streamline-upwind/Petrov-Galerkin formulations for convection
dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer
Methods in Applied MEchanics and Engineering, Vol.32.1,pp.199-259,1982.
4) C.W. Hirt and B.D. Nichols: “Volume of fluid method for the dynamics of free boundaries”, J.
Comp. Phys., 39,pp.201- 225,1981.
5) Moriguchi S., Yashima A., Sawada K., Uzuoka R. and Ito M., Numerical simulation of flow failure ofGeomaterials based on fluid dynamics, Soils and Foundations, Vol.45, No.2, pp.155-166, 2005.
6) Oda K., Moriguchi S., Kamiishi I., Yashima A., Sawada K. and Sato A., Simulation of snow avalanche
model test using computational fluid dynamics, Annals of Glaciology, pp.57-64, Vol.52, No.58, pp.57-64,
August 2011.
7) 山口裕矢, 高瀬慎介, 森口周二, 寺田賢二郎, 小田憲一, 上石勲: 非ニュートン流体モデルを用いた雪崩の 3 次元非構造有限要素解析, 日本計算工学会論文集, Vol.2017, p.20170011,2017.
8) 上石勲, 山口悟, 根本征樹, 佐藤研吾, 伊藤陽一, 中村一樹, 内山庄一郎:災害調査那須町雪崩災害第 2 回調査（速報）, http://www.bosai.go.jp/saigai/2017/pdf/20170410_01.pdf(2018.2.25
閲覧).

83