自然災害科学 J. JSNDS 41 特別号 95-110（2022 ）多面体要素の個別要素法による粒状体シ
ミュレーションへの転がり摩擦の導入に関する基礎的研究
野畑舞愛郎1・古川愛子1・清野純史1

Fundamental Study on Introduction of Rolling Friction
to Granular Material Simulation using Distinct
Element Method with Polyhedral Elements
Maaro Nobata1, Aiko Furukawa1 and Junji Kiyono1
Abstract
Circular elements are frequently used in the DEM simulations of granular materials. Since
circular elements are easy to rotate, rolling friction has been used as a method to suppress
element rotation. On the other hand, in DEM simulations using polyhedral elements, rolling
friction is rarely used because the shape of the element itself has the effect of suppressing
rotation. However, if the element shape is different from that of the actual granular material, therolling resistance is also different. Therefore, the analysis results are different from the actual
behavior. Therefore, in this study, rolling friction is newly introduced into the DEM using thepolyhedral elements. It was confirmed whether the proposed rolling friction works in the DEM
using the polyhedral elements and whether the introduction of rolling friction reduces thedifferences in simulation results using elements with different shapes. Specifically, simulations
were performed to determine the angle of repose for elements with two different shapes(decahedral and octahedral) with varying coefficients of sliding friction and rolling friction toinvestigate the effect of rolling friction on the angle of repose. The effects of sliding and rolling
friction on the analytical results were verified through failure simulations of a rock mass.
キーワード：多面体要素の個別要素法，転がり摩擦，要素形状，安息角，斜面崩壊
Key words: DEM with polyhedral elements, rolling friction, element shape, angle of repose, slope failure

1 ．はじめにの解析手法の 1 つである。個別要素法は連続体の1971年に Cundall1）によって提唱された個別要解析手法では表現の難しい粒状体の挙動2-6）や構素法（Distinct Element Method：DEM）は離散体造物の破壊後の挙動2,7,8）のシミュレーションがで

1

京都大学大学院工学研究科都市社会工学専攻
Department of Urban Management, Kyoto University

95

96

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

きる手法として有用性が認められてきた。
個別要素法を用いた粒状体シミュレーションの精度と落石シミュレーションの解析精度の間に関係のあることが明らかとされている22）。

ほとんどのケースでは円・球形要素が用いられて

東ら23）は，石造アーチ橋における中詰め材（砕いるが2,3,5,6），多面体要素を用いた個別要素法の例

石）の大きさが壁面の崩壊し易さに及ぼす検討のとしては，古川らによって提案された改良版個別

中で，改良版個別要素法を用いて振動台実験の再9）要素法（Refined DEM）
，BLOCK3D10），iDEM11），

現解析を行っている。なお，砕石は十面体要素で商用ソフトである 3DEC ，オープンソースのモデル化しており，改良版個別要素法では転がりBlaze-DEM13），Yade14）などがある。

摩擦を考慮していない。解析により，振動台実験

12）円・球形要素を用いた個別要素法は，要素間ので確認された“砕石が大きいほうが壁石の耐震性

接触判定のアルゴリズムが簡単であるため計算時が高い”という傾向は再現できたが，加振中の挙間が大幅に短縮されるというメリットがある。し動を完全に再現することはできなかった。その理

かし，実際の粒子形状は完全な円や球ではないた

由の 1 つとして，実際の砕石の形状と解析モデル

め，円・球形要素を用いた個別要素法では要素のの要素形状の違いがあると考えた。複雑な形状を

回転量が実挙動に比べて大きくなりすぎるという持つ砕石を十面体要素でモデル化したが，砕石同問題が生じる。そこで，円・球形要素を用いた個士の形状に起因するかみ合いを表現することがで別要素法における回転を抑制するための方法としきないため，転がり抵抗に差が生じたのではないて，転がり摩擦が提案された 。

かと考えた。解析結果を実験結果に近づける方法

15）円・球形要素を用いた個別要素法における転がとして，①実際の砕石の複雑な形状を忠実にモデり摩擦に関しては数多くの研究がなされてきル化する方法，②転がり摩擦の導入により形状にた15-21）。Wensrich ら16）は，個別要素法では精緻な起因する転がり抵抗の差を解消する方法，の 2 つモデル化が難しい複雑な形状を持つ粒子の挙動を，が考えられる。近年，複雑な形状が扱える解析

粒子形状をモデル化する代わりに円形要素に転がコードも存在4,14）するが，全ての砕石や岩石・岩り摩擦を導入することで簡易に再現できるかどう

塊の形状を 3D スキャナーで測量して形状を完璧か研究を行った。中瀬ら20）は，球形要素によってにモデル化することは現実的でない。そこで本研

地盤モデルの作成を行った場合，接線方向と法線

究では，多面体要素で形状を大まかに再現した上方向の接触力だけでは，粒子間の滑り摩擦係数をで，多面体要素に転がり摩擦を導入することで，どれだけ大きくしても地盤強度定数としての内部要素形状の違いによる解析結果の差異を埋めるこ

摩擦角は25°程度に留まってしまうことと，転がとができないかと考えた。

り摩擦の必要性を示した。

本研究は，多面体要素と実際の粒子形状の差を，一方，多面体要素を用いた個別要素法では，要多面体要素に転がり摩擦を導入することで解消で素形状自体に回転を抑制する効果があることもあきるか検討する。球形要素でなく多面体要素に転り，転がり摩擦を考慮した解析は行われていない。

がり摩擦を導入する理由は，要素形状が実際の粒しかし，要素形状と実際の粒子形状が対応してい子形状と大きく異なる場合に，転がり摩擦だけでない場合は，転がり抵抗も異なるため，解析結果差を埋めることができないためである。 例えばと実挙動に差異が生じると考えられる。

後述する十面体と八面体を用いた安息角のシミュ

なお，円・球形要素を用いた個別要素法においレーションでは，転がり摩擦係数の増加に対してて，転がり摩擦を用いずに，複数の円・球形要素

安息角の増加は頭打ちとなり，転がり摩擦だけでを結合して複雑な形状を模擬する解析も行われて再現できる安息角に上限があることがわかった。

いる。実際の粒子サイズに比べて径の小さい

また，球形に近い十面体の方が再現できる安息角

円・球形要素をオーバーラップを許容して多数結の上限が小さいこともわかった。十面体と八面体合するほど形状の再現精度が高まり，形状の再現は既に転がり抵抗を有しているため，転がり抵抗

16,22）

97

自然災害科学 J. JSNDS 41 特別号（2022 ）のない球形要素の場合，転がり摩擦だけで再現で要素に対して安息角を求めるシミュレーションをきる安息角の上限は更に小さくなると考えられる。

行い，安息角に及ぼす転がり摩擦の影響を比較検

つまり，転がり摩擦は万能ではなく，球形要素を証する。 4 章では，岩塊を用いた斜面崩壊シミュ

用いた場合，転がり摩擦だけで形状の差をカバー

レーションを行い，異なる形状の要素でも，転がすることは困難な場合があると考えられる。従っり摩擦により安息角が等しくなるよう調整するこて，岩石・岩塊等の粒子形状の大まかな特徴を捉とによって，同じ崩壊挙動の再現ができるか検証

えた多面体で近似した上で，それでも残る形状のを行った。また同じ要素形状で転がり摩擦を変化

差を転がり摩擦によって埋めることを考えている。

させた際の斜面上での挙動の変化についても検討

本研究で用いる改良版個別要素法では，ばね定数を行った。

は物性値から決定できるため，調整が必要なパラメータは転がり摩擦係数だけとなる。実際の粒子の安息角を再現するように決定することを考えて

2 ．改 良版個別要素法への転がり摩擦の
導入
2. 1

いる。
本研究は，改良版個別要素法に新たに転がり摩改良版個別要素法の概要（1）個別要素法との違い擦を導入すること，転がり摩擦の導入によって転個別要素法では，構造物を剛体要素の集合体とがり抵抗が向上したか検証すること，要素形状にしてモデル化する。要素同士が接触した際に，接よる解析結果の差異を転がり摩擦によって埋める

触要素間に仮想のばね・ダッシュポットを設置し，ことができるのか検証を行うことを目的とする。

要素同士の貫入距離に基づいた力を伝達すること

具体的には 2 種類の形状（十面体と八面体）の要で，要素間の相互作用を簡易にモデル化する。個素について滑り摩擦係数と転がり摩擦係数を様々別要素法では，時々刻々と変化する各要素間の接に変化させた安息角を求めるシミュレーションを

触問題を個別に取り扱いながら解析が行われる。

行い，安息角に及ぼす転がり摩擦の影響を調べた。

円形要素を用いた個別要素法では，要素間のば

そして，要素形状が異なっていても，転がり摩擦

ねのばね定数は Hertz の接触理論等に基づき材料を考慮することで安息角を一致させることができ物性値から決定される24）。しかし，多面体要素をるか検証した。また，十面体と八面体の要素を用用いた個別要素法には，物性値から要素間をつないて斜面崩壊挙動の数値解析を行った。安息角がぐばねのばね定数を理論的に決定できないという等しくなるように転がり摩擦を調整することで，問題があった。

異なる要素形状を用いても同じような斜面崩壊挙動を得ることができるか検証を行った。

本研究で用いる改良版個別要素法は，この問題を解決したもの9）である。改良点としては，要素

本論文の構成は以下の通りである。 2 章では，表面を複数のセグメントに分割し（図 1 ⒜），そ

改良版個別要素法の概要と転がり摩擦のモデル化

れぞれのセグメントの代表点にばねとダッシュ

手法について述べる。 3 章では， 2 種類の形状のポットを取り付けた（図 1 ⒝）という点にある。

⒜

要素表面のセグメント分割
図1⒝

ばね・ダッシュポットの設置

改良版個別要素法の概要

98

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

これにより，セグメント毎の力のつり合いに基づントの法線方向と接線方向それぞれに設置される。

きばね定数を物性値から理論的に決定することが減衰定数を h とし，単位接触面積あたりの減衰可能となった。なお，構造物を剛体要素の集合体としてモデル化する点は，従来の個別要素法と同じである。
なお，個別要素法を改良した手法は全て「改良版」個別要素法となるが，本研究で用いる「改良版個別要素法」の改良点は多面体の個別要素法のばね定数を 3 次元応力－ひずみ関係に基づき物性値から決定できるようにした点である。既往の文

係数は次のように表わされるとする。
cn＝2h maνekn

cs＝2h maνeks（3）

ここに，maνe は単位接触面積あたりの要素 A，Bの質量の和であり，次式で表される。
maνe＝ρAℓA＋ρBℓB（4）

ρA，ρB は要素 A，B の質量密度である。

献との連続性から本稿でも「改良版個別要素法」
2. 2

と呼ぶ。
（2）ばね定数の決定

Zhou らの転がり摩擦モデル

転がり摩擦に関しては，以前より多くの研究がばねは，要素表面に対して，法線方向（n）となされ様々なモデルが提案されてきた15-19, 21）。本

接線方向（s）の両方に取り付けられる。法線，接研究では，様々な個別要素法の研究で用いられて線方向の単位面積あたりのばね定数は次式で表さおり，他の転がり摩擦モデルと比較して簡単なれる。

Zhou ら18）のモデルを参考にする。このモデルは，
E

kn＝（1－ν 2）ℓks＝E
2（1＋ν）ℓ

法線方向の接触力に比例した抵抗モーメントを回（1）

転方向と逆方向に作用させるモデルを一般化したものである。

ここに，E は要素の弾性係数，ν はポアソン比，ℓは要素重心から表面までの距離である。式（ 1 ）Zhou らの円形要素を用いた DEM では，次式を用いて転がり抵抗モーメントが計算されている。

は，セグメント毎の力のつり合い式に 3 次元応力
―ひずみ関係を適用して導出したものである9）。

Tr＝－μr Rr |Fn|

2 つの要素 A，B が接触しているとする。要素
A，B の弾性係数を EA，EB，ポアソン比を νA，νB，

ωrel|ωrel|（5）

ここに，Tr は転がり抵抗モーメント（ベクトル），

重心から表面までの距離をℓA，ℓB で表す。本研

μr は転がり摩擦係数，Rr は有効転がり半径（要素 i

究では，式（ 1 ）で求めたばねが直列につながっと要素 j の半径が ri と rj のとき，Rr＝ri rj /(ri＋rj )），ていると想定し，要素間の単位面積あたりのばね

Fn は法線方向の接触力，ωrel は接触している 2 要

定数を次式で与えることとする。

素の相対角速度ベクトル，ωrel /|ωrel| は，相対角

kn＝

ks＝

1
ℓA
ℓB
＋
EA /（1－νA2） EB /（1－νB2）

速度ベクトル ωrel と平行な方向の単位ベクトルで
（2）

1
ℓA

ℓB
＋
EA /2（1＋νA ） EB /2（1＋νB ）（3）ダッシュポットの減衰係数の決定ある。
式（ 5 ）は円形要素を対象としたものであるた
め，本研究ではこれを多面体要素に拡張する。
2. 3

改良版個別要素法への転がり摩擦の導入式（ 5 ）では円形要素の相対角速度ベクトルの逆方向に転がり抵抗モーメントを与えているが，

ダッシュポットはばねと並列に設置され，要素3 次元の多面体要素を用いた改良版個別要素法に同士が衝突する際のエネルギー消散を表現するた

おいて，接触要素間の相対角速度ベクトルを求めめのものである。ばねとダッシュポットはセグメ

るには非常に複雑な計算が必要になってしまう。

99

自然災害科学 J. JSNDS 41 特別号（2022 ）要素 ie（時間 t-t）
重⼼ G

重⼼ G

𝑟𝑟⃗�

𝑟𝑟⃗���
セグメントの代表点（接触点）
要素 je （時間 t）

𝑟𝑟⃗���

相対変位ベクトルの回転成分の増分

重⼼ G

𝑢𝑢
�⃗� � 𝑟𝑟⃗𝑡𝑡 � 𝑟𝑟⃗𝑡𝑡�𝑡𝑡

相対速度ベクトルの回転成分𝑟𝑟⃗� 𝑟𝑟⃗���𝑢𝑢��⃗� �∆𝑡𝑡

要素 je から⾒た位置関係
図2

要素 ie
（時間 t）

相対速度ベクトルの回転成分

そこで，相対角速度ベクトルを用いる代わりに，した。要素 je が受ける転がり抵抗モーメントは，

図 2 に示すように，相対速度ベクトルの回転成分式（ 7 ）の rt を要素 je の重心から接触点に向かうベu̇θ を求めて，要素重心から接触点に向かうベク

クトルに変更することによって求めることとした。

トル rt と，相対速度ベクトルの回転成分方向の単位ベクトルの外積で求まるベクトルの逆方向に

Tr＝－μr |Fn|rt×

転がり抵抗モーメントを与えることとした。

u̇θ
|u̇θ|

（7）

改良版個別要素法では，要素 ie の各セグメン

式（ 7 ）では，要素重心から接触点までの距離が要

トの代表点（ばねを付けている点）と要素 je の接

素 ie と要素 je で異なるため，それぞれの要素が受

触判定を行っている。従って，接触点は要素 ie のける抵抗モーメントの大きさは異なる。 2 要素がセグメントの代表点となる。図 2 は，時間 t にお受ける抵抗モーメントの大きさが異なることにつける要素 ie と要素 je の位置関係と，時間 t －∆t にいて考察する。実際の 2 つの粒子が噛み合うとき，

おいて要素 je から見た要素 ie の位置関係を示す。

作用・反作用の法則により互いに逆向きで大きさ

時間 t－∆t において要素 ie の重心 G からセグメンの等しい力を及ぼし合っている。しかし， 2 粒子トの代表点に向かうベクトルを rt－∆t，時間 t におそれぞれの重心から噛み合い箇所までの距離が異いて要素 ie の重心 G からセグメントの代表点になるため，噛み合いによって生じるモーメントの向かうベクトルを rt とする。このとき，相対速大きさは必ずしも 2 粒子で同じではない。従って，度ベクトルの回転成分を式（ 7 ）により計算されるモーメントの大きさが 2u̇θ＝rt－rt－∆t∆t

要素で異なるのは致命的な問題ではないと考える。
（6）で表すこととする。そして，rt と相対速度ベクトルの回転成分方向の単位ベクトル u̇θ /|u̇θ| との外積の逆方向に転がり抵抗モーメントを与えた。

3 ．転がり摩擦が安息角に与える影響
3. 1

解析概要

本章では，転がり摩擦を導入した改良版個別要素法を用いて，粒状体の安息角シミュレーション

具体的には，式（ 7 ）に示す転がり抵抗モーメを行い，転がり摩擦が安息角に与える影響についントを要素 ie に与えることによって，要素 ie がて検討する。要素形状として，十面体と八面体の受ける転がり摩擦による抵抗力を表現することと2 種類の形状を用いた。改良版個別要素法は多面100

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

体要素を基本としており，完全な球形要素を用いであり，この縮尺を保った状態で基準辺の長さにて解析を行うことが難しいため，球形に近い転がより十面体，八面体のサイズの定義を行う。本研

りやすい要素形状として十面体要素を採用した。

究の解析では，斜面を転がる岩塊を想定し，基準

また，転がりにくい要素形状として八面体要素を

辺が 1 m の要素モデルを用いた。

用いた。様々な滑り摩擦係数と転がり摩擦係数の3. 3組み合わせに対して安息角を計算し，転がり摩擦

安息角を求めるための解析手順

本研究の解析手順について説明する。解析は 2の導入によって任意の安息角を再現できるかどう段階で実施した。

か検証した。

（1）第 1 段階：要素の充填
3. 2

要素モデル

安息角を求めるための解析に先立ち，十面体または八面体要素の集合体を充填させるための解析

本研究で作成した十面体と八面体の要素モデル

4 ⒜⒝に示すように，細長い直方を図 3 に示す。各辺の長さの比は図に示した通り基準辺 を実施した。図1.0

基準辺 1.0

0.5

0.3

0.4

0.3

0.3

0.4

0.3

⒜

1.00.5
0.4
1.0
0.4

十面体要素

基準辺 1.0
基準辺 1.0

0.5
1.00.5

⒝
図3

1.0

八面体要素

要素モデル（左から，立面図，上面図，側面図）

y
x

z

z
x

⒜

解析モデル（ 3 次元）
図4

⒝

y

xy，yz，xz 平面から見た岩塊要素
モデルの配置状況（十面体の例）⒞

充填後の要素岩塊要素の充填に用いる容器と岩塊の解析モデルと充填後の様子

101

自然災害科学 J. JSNDS 41 特別号（2022 ）体の容器の上方に要素を規則的に配置し，重力になお，図 4 ⒝のような規則的な配置以外に乱数

よって自由落下させて，容器内に充填させた。

を用いたランダム配置や，自由落下以外に振動を容器のサイズは幅方向（ x 方向）に7.5 m，奥行与えて充填させる方法など様々試したが， 4 方向

き方向（ y 方向）に7.5 m，高さ方向（ z 方向）にからの平均を取れば安息角の算定結果に大きな影

35 m である。この容器内に，図 4 ⒝に示すよう響はなかったことから，規則配置と自由落下を組に，基準辺 1 m の岩塊要素モデルを x，y 方向そみ合わせた簡単な方法を採用した。

れぞれに 7 個ずつ（ 1 m× 7 個＝ 7 m），z 方向に10個（ 1 m×10個＝10 m），計490個を容器底面か

3. 4

解析諸元

ら20～30 m の高さに配置した。岩塊要素を自由

本解析では，十面体要素で構成される岩塊と，落下させて充填する際に，落下した要素が散らば

八面体要素で構成される岩塊それぞれに対して，

り隙間を埋めるるように x，y 方向それぞれ0.5 m

滑り摩擦係数と転がり摩擦係数を様々に変えた数ずつ余裕を持たせた容器の大きさに設定している。

値解析を実施し，安息角を算出した。

具体的な解析手順について述べる。容器を全自滑り摩擦係数 μs は0.1，0.2，0.3，0.4，0.8の 5

由度固定し，容器の底面から高さ20～30 m の位

通り，転がり摩擦係数 μr は0.0，0.2，0.3，0.5，

置に岩塊の集合体を配置し，要素重心に重力加速0.8の 5 通りとした。なお，滑り摩擦係数と転がり

度を鉛直下向き（z 軸の負方向）に与え，自由落

摩擦係数は全ての要素間で同じ値を利用している。

下させる。岩塊が自由落下し，容器および他の岩要素間の衝突時のエネルギー消散を表すための塊要素と衝突し，岩塊の挙動が容器内で安定する

ダッシュポットの減衰定数には，臨界減衰（式までの10秒間の解析を行い，岩塊を充填させた。

（ 3 ）における減衰定数 h を1.0に設定）を採用した。

（2）第 2 段階：安息角の算出

砕石の安息角の実験を対象にした数値解析23）にお第 1 段階の解析後，充填された容器内の岩塊モ
デルの座標を読み取り，平らな要素の上に配置す

いて，臨界減衰を採用したときに実験結果を精度よく再現したことから臨界減衰を採用した。

る。再び，重力下に置くと，これまで岩塊を支え解析に使用した密度，ヤング率，ポアソン比をていた容器の支えがないため，図 5 のように斜面

表 1 に 示 す。 表 1 に 示 す 物 性 値 を， 式（ 1 ）～状に崩れる。崩れた岩塊モデルの安息角は，解析（ 4 ）に代入することで，ばね定数とダッシュポッ

結果の画像からオンライン分度器25）を用いて計測

トの減衰係数を決定した。容器の物性値はアクリした。測定場所によって測定結果に差が出ること

ル板を想定して一般的な値を採用した。

を防ぐため，x-z 平面と y-z 平面でそれぞれ両方

岩塊の密度，ポアソン比は既往の研究を参考に向から測定し 4 回の測定結果の平均値をそのモデ設定した26）。ヤング率は，実際のヤング率の1/100

ルの安息角としている。 4 方向からの測定結果のに設定した。改良版個別要素法では，ヤング率が差は，最小で 1 度，最大で 5 度であった。

大きいほど短い解析時間間隔が要求され，解析時間が長くなる。本研究では，解析時間短縮のために，実際より小さいヤング率を使用した。ヤング率（剛性）の値は，圧縮試験のように要素間に大きな圧縮力がかかり要素同士がオーバーラップする

表1

図5

解析結果と安息角の測定の様子解析に用いる物性値

要素

密度（t/m3）

ヤング率（N/m2）

ポアソン比

岩塊

3.2

2.0×108

0.18

アクリル板

1.19

3.0×1090.35

安息角（度）102

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

40に対する安息角の解析結果を示す。仮説の通り，

30

ヤング率は安息角にほとんど影響しないことを確認できたことから，本研究では表 1 に示す物性値20を採用した。

10

ヤング率2.0x10^8 (N/mm^2)

ヤング率2.0x10^10 (N/mm^2)

0

0

図6

0.2

0.4

0.6

3. 5

解析結果

例として，十面体要素で滑り摩擦係数が0.4の0.8

転がり摩擦係数

ときの，第 2 段階の解析結果を図 7 に示す。滑り岩塊のヤング率が安息角に及ぼす影響（十面体要素，滑り摩擦係数0.4）

摩擦係数が同じでも，転がり摩擦係数が増加する
につれて安息角が増加する傾向が見てとれる。つ
まり， 2 章で提案した方法によって転がり抵抗がことで変形が生じるような解析では結果に大きな増加していることを確認することができる。

影響を及ぼすが，要素同士の相対運動（滑りや回滑り摩擦係数を0.1，0.2，0.3，0.4，0.8の 5 通転）によって変形が生じるケースでは支配的でなりに対して，転がり摩擦係数を0.0，0.2，0.3，0.5，

いことが多い27）。本研究の解析も，後者のケース

0.8と 5 段階に増加させたときの安息角を読み取であると考え，実際より小さいヤング率を採用しりグラフにまとめたものを図 8 に示す。同図⒜は2た。仮説の確認のため，ヤング率が2.0×10 N/m

十面体要素，同図⒝は八面体要素の結果である。

と2.0×108 N/m2 のときの安息角の算定結果を比

図 8 より， 5 通りの滑り摩擦係数すべてにおい較した。図 6 に，岩塊を十面体要素でモデル化し，て，十面体要素でも八面体要素でも，転がり摩擦

滑り摩擦係数を0.4としたときの 2 通りのヤング率係数の増加により安息角が増加していることがわ

10⒜

転がり摩擦係数0.0⒝

転がり摩擦係数0.2⒟

転がり摩擦係数0.5⒠

転がり摩擦係数0.8

40

40

30

30

20
10
0

0

転がり摩擦係数0.3安息角を算出するための解析結果の例（十面体要素，滑り摩擦係数0.4）

安息角（度）

安息角（度）

図7⒞μs=0.1μs=0.2μs=0.4μs=0.80.2

0.4
転がり摩擦係数μs=0.3
0.6

0.8

20
10
0

0

μs=0.1μs=0.2μs=0.4μs=0.80.2⒝

0.4
転がり摩擦係数

⒜

十面体要素

図8

様々な滑り摩擦係数 μs と転がり摩擦係数に対する安息角の解析結果

八面体要素

μs=0.3

0.6

0.8

103

自然災害科学 J. JSNDS 41 特別号（2022 ）かる。すなわち， 2 章で提案した方法によって，改良版個別要素法でも転がり抵抗を表現できていることを確認することができた。
また，図 8 ⒜と⒝の比較より，滑り摩擦係数が等しくても要素形状が異なる場合，従来の転がり

4 ．斜面崩壊シミュレーション
4. 1

解析概要

3 章の解析結果から，岩塊の要素間に転がり摩擦を導入することによって，より広い範囲の安息角を表現できることが分かった。

摩擦を考慮しない改良版個別要素法（転がり摩擦

本章では，より実践的な問題として，斜面にお係数が0.0のケースに相当）では異なる安息角とける岩塊の崩壊挙動の解析を行った。要素形状やなっており，球形に近い十面体要素の方が安息角

滑り摩擦係数と転がり摩擦係数の値によって，斜が小さくなっている。しかし，転がり摩擦係数を面上の岩塊の挙動がどのように変化するのか検証

調整することによって，異なる要素形状でも同程する。また，滑り摩擦係数が等しく形状が異なる

度の安息角を再現できるケースのあることも確認

要素（十面体要素と八面体要素）を用いたときに，できる。

安息角が等しくなるように転がり摩擦係数を設定

さらに， 5 通りの滑り摩擦係数すべてにおいて，十面体要素でも八面体要素でも，転がり摩擦係数

することで，斜面崩壊挙動の解析結果が等しくなるかどうかを検証する。

が0.0から0.2に増加した時に安息角は最も大きく増加している。そして，転がり摩擦係数が大きく

4. 2

解析モデル

なるにつれ，転がり摩擦係数の増加に対する安息

本解析に用いる岩塊の要素は前節で用いたもの

角の増加割合は緩やかになることが確認された。

と同じである（図 3 ）。基準辺の長さも前節と同

つまり，滑り摩擦係数や転がり摩擦係数を大きく様に 1 m とした。

していっても再現できる安息角の大きさには限界

斜面モデルを図 9 ⒜に示す。斜面モデルは，水が存在することから，転がり摩擦係数だけで要素

平面と斜面の 2 つの平面で構成される。水平面は，サイズの違いを完全に埋めるには限界があること

幅40 m，長さ30 m である。斜面の幅も40 m であもわかった。

り，水平面と斜面のなす角は30°である。

以上の解析結果から，安息角の大きさという点においては，転がり摩擦の導入効果を確認することができ，滑り摩擦だけでは再現することのでき座標軸は，水平面の長さ方向に x 軸，幅方向にy 軸，高さ方向に z 軸とした。
岩塊は，x 軸方向に 6 列（ 6 m），y 軸方向に10なかった広範囲の安息角の再現を行うことができ列（10 m），z 軸方向に 3 段（ 3 m）の計180個の要た。しかし，転がり摩擦だけでは再現できる安息

素で構成した。xy 平面，yz 平面，xz 平面から見角に上限があることもわかった。球形に近い十面た岩塊要素の配置状況を図 9（b）に示す。なお，

体要素の方が八面体要素よりも再現できる安息角岩塊の個数は既往研究28）の粒径大のモデルを参考の上限が小さいことから，球形要素の場合は，転にした。岩塊を包含する 6 m×10 m× 3 m の直がり摩擦だけで再現できる安息角の上限は更に小方体を考えたとき，斜面上の水平面からの高さがさくなると推察される。このことから，粒子同士

24 m の位置（斜面に沿って長さ48 m の位置）に，の噛み合いが大きく安息角が大きい粒状体の再現

岩塊を包含する直方体の x 座標と z 座標が最小とに球形要素を用いた場合，転がり摩擦だけで所望なる辺がくるように岩塊の初期位置を設定した。

の回転抵抗を表現するのは困難な場合があると考図 9 に示す岩塊要素に重量加速度を与え，固定えられる。この様な場合は，粒子形状の大まかなされた斜面と平面上の岩塊の挙動を解析した。

特徴を捉えた多面体要素でモデル化した上で，そ岩塊の物性値は前節と同じ表 1 の値を用いた。

れでも残る形状の差を転がり摩擦によって埋める

斜面と水平面を表す要素も岩塊と同じ物性値としことが有用であると考えられる。

た。要素同士が衝突する際のエネルギー消散を表すためのダッシュポットの減衰定数には，臨界減

104

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

岩塊
3m

長さ48m
10m
斜面

長さ30m

6m

高さ24m

z軸
30°

x軸

y軸

水平面
幅40m
⒜

y

z

x
⒝

全体モデル

z
x

y

岩塊要素の配置状況（十面体の例，左から，xy 平面，yz 平面，xz 平面）
図9斜面の解析モデル衰（式（ 3 ）における減衰定数 h を1.0に設定）を

滑り摩擦係数は0.8と0.4の 2 通りと設定してい仮定した。実際の岩塊の数値解析の際は，実際のる。斜面の角度を30°（tan 30°＝0.577）としてい

跳ね返り係数を元に設定する必要がある。

る た め， 滑 り 摩 擦 係 数 が0.8で0.577よ り 大 き い
ケース 1 ～ 6 は転がりが卓越することを想定した

4. 3

解析ケース

ケ ー ス， 滑 り 摩 擦 係 数 が0.4で0.577よ り 小 さ い

解析ケースとして，表 2 に示す10通りのケース

ケース 7 ～10は滑りが卓越することを想定したを設定した。表中の安息角の範囲とは， 3 章の安ケースとしている。

息角の測定において， 4 方向からの測定結果の範囲を示したものである。
表2

斜面崩壊シミュレーションの解析ケース

ケース

要素形状

滑り摩擦係数

転がり摩擦係数

安息角の平均値（°）

安息角の範囲
（°）

ケース 1

八面体

0.8

0

28.75

27～31

ケース 2

十面体

0.8

0

21.75

21～22

ケース 3

八面体

0.8

0.2

30.5

28～33

ケース 4

八面体

0.8

0.5

32.25

31～33

ケース 5

十面体

0.8

0.2

26

25～28

ケース 6

十面体

0.8

0.5

30.5

29～33

ケース 7

八面体

0.4

0.5

31.5

30～33

ケース 8

十面体

0.4

0.8

31.5

30～34

ケース 9

八面体

0.4

0

24.5

23～25

ケース10

十面体

0.4

0.2

25

22～27

105

自然災害科学 J. JSNDS 41 特別号（2022 ）

4. 4

解析結果

大到達距離は，180個の岩塊のうち，全頂点の x転がりが卓越すると想定したケース 1 ～ 6 の結座標の最大値として求めた。重心の移動距離は，果を図10に，滑りが卓越すると想定したケース 7180個の岩塊の重心の x 座標の平均値から，初期～10の結果を図11に示す。いずれも，斜面上を転値（ 3 m）を差し引くことによって求めた。法尻がる全ての岩塊が静止した時の様子を示している。

までの到達時間とは，斜面と水平面が交差する地

解析結果の定量的な比較のために，岩塊の x 軸点に岩塊がはじめて到達した時間として求めた。

方向最大到達距離，岩塊重心の x 軸方向移動距離，

結果を表 3 に示す。

法尻までの到達時間の 3 つの指標を算出した。最

⒜

ケース 1

⒝

ケース 2

⒞

ケース 3

⒟

ケース 4

⒠

ケース 5

⒡

ケース 6

図10

転がりが卓越するケース 1 ～ 6 の斜面上で静止した岩塊の様子

106

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

⒜

ケース 7

⒝

ケース 8

⒞

ケース 9

⒟

ケース10

図11

滑りが卓越するケース 7 ～10の斜面上で静止した岩塊の様子
表3

各ケースの岩塊の斜面上での挙動

ケース

岩塊の x 軸方向最大到達距離（m）

岩塊重心の x 軸方向移動距離（m）

法尻への到着時間
（sec）

ケース 1
ケース 2

22.18

2.87

到達しない

44.15

17.05

9.5

ケース 3

13.38

1.64

到達しない

ケース 4

15.29

1.39

到達しない

ケース 5

41.71

6.94

13

ケース 6

15.77

3.13

到達しない

ケース 7

43.17

33.2

15

ケース 8

44.37

33.01

14.6

ケース 9

42.8

34.28

15.2

ケース10

43.82

32.92

14.2

（1）ケース 1 ～ 6 の解析結果

慮していないことから，要素形状が解析結果に表

ケース 1 ～ 6 では，滑り摩擦係数0.8が斜面のれている。球形に近く安息角が小さい十面体要素

角度の正接0.577よりも大きいため，岩塊は滑るの方が転がりやすく，最大到達距離も重心の移動

ことなく斜面上を転がった（図10，表 3 ）。

距離も大きくなっており，法尻にも到達している。

ケース 1 ， 2 は，滑り摩擦係数が0.8で等しく，
要素形状が異なるケースである。転がり摩擦を考

一方の八面体要素は転がりにくいため，水平面にたどり着くことなく斜面上で静止した。

107

自然災害科学 J. JSNDS 41 特別号（2022 ）ケース 1 ， 3 ， 4 は，滑り摩擦係数が0.8で等ス 6 のほうがわずかに岩塊の散らばりが大きく，

しく，要素形状も八面体で等しく，転がり摩擦係

到達距離も大きくなっているが，個別要素法の解数の影響を比較するためのケースである。これら

析結果のばらつきを考えると，両ケースはよく似のケースでは要素が水平面に辿り着くことなく解た傾向を示しており，要素形状が異なっても安息

析の早いタイミングで斜面上で静止した。表 2 に角が一致するように転がり摩擦係数を設定するこ

各ケースの安息角を示しているが，ケース 1 ， 3 ，とで，概ね近い結果が得られたと言える。

4 は安息角が斜面の角度（30°）より高いため静止

水平面まで到達したのは，ケース 2 とケース 5

したものと考えられる。ケース 1 の安息角の平均の 2 ケースであった。これらは安息角が斜面の角値は30°より小さいが， 4 方向から測定したときの度30度より小さいケースであった。

安息角の範囲は27～31°となっており，最大値の（2）ケース 7 ～10の解析結果

31°が斜面の角度を上回っていることから静止しケース 7 ～10では，滑り摩擦係数0.4が斜面のたものと考えた。また表 3 からわかるように，転角度の正接0.577よりも小さいため，岩塊は斜面

がり摩擦の最も小さいケース 1 が，最終ステップ

上を転がることなく滑動した。

における到達距離は最大となった。もともと転がケース 7 ～10は要素形状と安息角が異なるもの

りにくい形状を有する八面体要素ではあるが，転の，全てのケースで最大到達距離，重心の移動距

がり摩擦係数を考慮したケース 3 ， 4 ではケース

離，法尻までの到達時間に違いはほとんど見られ

1 よりもさらに転がり難くなり，斜面上の移動距なかった。ケース 7 ， 8 では安息角が31.5度と斜離が減少する様子を確認することができた。転が面の角度より大きいものの，要素群は斜面上を滑り摩擦係数が最大で安息角の平均値も大きいケーり落ちた。

ス 4 が，土塊重心の移動距離が最小となったが，これらの結果から，滑り摩擦係数0.4が斜面の最大到達距離についてはケース 3 の方が最大と正接0.577よりも小さい場合，転がりよりも滑りなった。これは，図10⒞⒟の比較からわかるようが支配的となるため，転がり摩擦が岩塊の挙動にに，ケース 4 では 1 個だけ遠くまで転がった要素

対して与える影響はほとんど見られなかった。

があるため最大到達距離が逆転したと考えられる。

（3）まとめケース 2 ， 5 ， 6 は，滑り摩擦係数が0.8で等以上の10ケースの解析結果から，岩塊の要素群

しく，要素形状も十面体で等しく，転がり摩擦係

が斜面上を滑るのか，それとも転がるのかは滑り数の影響を比較するためのケースである。転がり

摩擦係数によって決まり，斜面の一番下まで転が摩擦の小さいケース 2 ， 5 では，多くの要素が斜り続けるか，それとも途中で静止するかは，安息面上に留まってしまったが，水平面まで到達する角によって決まる傾向が確認された。

要素もあった。転がり摩擦が最大のケース 6 では，滑り摩擦係数が大きく転がりが卓越するケース

安息角が斜面の角度を上回っているため，すべて

では，転がり摩擦によって斜面崩壊挙動が異なり，の要素が斜面上に留まった。これらの 3 ケースを安息角が等しいケース 3 ， 6 では斜面上の岩塊挙比較すると，転がり摩擦係数が大きくなるほど最動が似た結果となった。

大到達距離，重心の移動距離ともに減少しており，一方で，滑り摩擦係数が小さく滑動が卓越する

法尻への到達時間も長くなるか到達しなくなってケースでは，要素の形状や転がり摩擦係数に関係いる。

なく岩塊の要素群が斜面状を滑り落ちるため，転最後に，要素形状は異なるが滑り摩擦係数と安がり摩擦の影響はほとんど見られなかった。

息角の等しいケース 3 ， 6 の比較を行う。両ケー

なお，安息角が等しくなるように転がり摩擦係

スともに安息角が斜面の角度より高くなっており，数を調整したケース 3 ， 6 において，斜面上での水平面まで到達することなく，解析開始直後の早挙動に差が生じたのは，多数の要素との接触状態

いタイミングで斜面上で静止した。十面体のケーが維持される平面上の安息角シミュレーションと，

108

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

斜面上で要素が広がりやすい斜面崩壊シミュレー

素法に転がり摩擦を導入し，その影響を安息角とションでは，要素間の接触状況が異なるからと考

斜面崩壊のシミュレーションによって明らかにし察している。本章では，岩塊要素の初期配置は図た。

9 ⒝に示した 1 通りとした。初期配置は前章の安まず，転がり摩擦を考慮した安息角のシミュ

息角のシミュレーション結果に大きな影響を及ぼ

レーションによって以下の知見が得られた。

さなかったことから，岩塊が斜面上を滑るかどう

1 ）十面体と八面体のいずれの要素でも，転がり

かや，ケース間の相対的な転がり易さの違いは初摩擦係数の増加により安息角が増加している

期配置を変えても同様の傾向が得られると推察しことが確認できた。提案手法によって，改良ているが，定量的な評価も含めて今後の検討課題

版個別要素法において転がり抵抗を表現できとしたい。

ていることを確認した。
2 ）転がり摩擦係数を調整することによって，異4. 5岩塊の崩落予測への適用に関する考察

実際の岩塊の形状は凹凸があり複雑である。岩石の崩落シミュレーションに関する既往研究より，なる要素形状でも同程度の安息角を再現できることがわかった。しかし，再現できる安息
角には上限があることを確認した。

形状の精度と解析の精度には関係があることが報

3 ）球形に近い十面体要素の方が再現できる安息

告されているが，多面体要素を用いたり，球要素

角の上限が小さいことから，円形要素に転がを多数結合しても岩塊の形状を完全に再現するこり摩擦を導入して再現できる安息角の上限はとは困難である。

更に小さいと推察した。転がり摩擦だけでは

本章では，十面体要素と八面体要素を用いて土塊を模擬し，土塊の斜面崩落シミュレーションを実施した。見かけの安息角が等しくなるように，任意の転がり抵抗を再現することはできず，要素形状も重要であることがわかった。
次に，転がり摩擦を考慮した斜面崩壊シミュ

十面体要素と八面体要素に転がり摩擦係数を設定

レーションによって，以下の知見が得られた。

することによって，斜面上の十面体と八面体の土4 ）岩塊の要素群が斜面上を滑るのか，それとも塊の挙動が概ね似た結果となった。

転がるのかは滑り摩擦係数によって決まり，

このことから，実際の岩塊の形状を完全に再現

斜面の一番下まで転がり続けるか，それともする代わりに，岩塊の安息角を再現できるように

途中で静止するかは，安息角によって決まる

転がり摩擦を設定することで，形状に誤差があっ

傾向が確認された。

ても，実際の岩塊挙動を予測できる可能性が示唆

5 ）滑り摩擦係数が大きく転がりが卓越するケーされた。ただし，実際の岩塊の安息角が必要で，スでは，転がり摩擦によって斜面崩壊挙動が現場試験をしたり，斜面勾配から決定する必要が異なり，要素形状が異なっても，転がり摩擦ある。

を調整して安息角を近づけることによって，今後は，模型実験により解析手法の妥当性を検斜面上の挙動を近づけることができた。

証する必要がある。岩塊の形状を大まかに再現す

6 ）滑り摩擦係数が小さく滑動が卓越するケース

る多面体要素を決定し，安息角を再現するようなでは，転がり摩擦の影響はほとんど見られな

転がり摩擦係数を試行錯誤的に決定することで，かった。

実際の岩塊と多面体要素の形状差を転がり摩擦に

以下に本研究の課題を挙げる。本研究では，仮よって埋め，崩落予測精度向上できるかどうか検

想の斜面モデルを作成して解析を行ったが，模型証する必要がある。

実験を行い，解析結果を実挙動と比較することで改良版個別要素法における転がり摩擦の有意性を

5 ．結論

本研究では，多面体要素を用いた改良版個別要検討したいと考えている。

自然災害科学 J. JSNDS 41 特別号（2022 ）

参考文献
1 ）P. A. Cundall：Formulation of a three-dimensional
distinct element model-Part I. A scheme to detect
and represent contacts in system composed ofmany polyhedral blocks，International Journal ofRock Mechanics and Mining Sciences &Geomechanics, Vol. 25, No. 3, pp.107-116, 1988.
2 ）伯野元彦：破壊のシミュレーション，拡張個別要素法で破壊を追う，森北出版株式会社，1997．
3 ）Matsushima, T., Saomoto, H., Tsubokawa, Y.,
Yamada, Y.: Grain Rotation versus Continuum
Rotation during Shear Deformation of Granular
Assembly, Soils & Foundations, Vol. 43, No. 4,
pp.95-106, 2003.
4 ）Xiang W., Zhen-Yu Y., Hao X., Dong S., Yun-Tian
F.: A spherical-harmonic-based approach to discreteelement modeling of 3D irregular particles,
International Journal for Numerical Methods inEngineering, Vol. 122, Issue 20, pp.5626-5655, 2021.
5 ）Beichuan Y., Richard A. Regueiro: Threedimensional discrete element method parallel
computation of Cauchy stress distribution overgranular materials, International Journal for
Numerical
andAnalytical
MethodsinGeomechanics, Vol. 43, Issue 5, pp.974-1004, 2019.
6 ）J. Fortin, O. Millet, G. de Saxce: Numericalsimulation of granular materials by an improved
discrete element method, International Journal
for Numerical Methods in Engineering, Vol. 62,
Issue 5, pp.639-663, 2004.
7 ）Daniele M., Matthew J. D., Andrea Penna, Distinct
element modelling of the in-plane cyclic responseof URM walls subjected to shear-compression,
Earthquake Engineering & Structural Dynamics,
Vol. 48, Issue 12, pp.1322-1344, 2019.
8 ）B. Chetouane, F. Dubois, M. Vinches, C. Bohatier:
NSCD discrete element method for modelling
masonry structures, International Journal for
Numerical Methods in Engineering, Vol. 64, Issue
1, pp.65-94, 2005.
9 ）Furukawa, A., J. Kiyono, K. Toki：Proposal of a
Numerical Simulation Method for Elastic, Failure
and Collapse Behaviors of Structures and its
Application to Seismic Response Analysis of
Masonry Walls, Journal of Disaster Research, Vol.
6, No. 1, 2011.

10910）Ghaboussi, Jamshid Ricardo Barbosa: Threedimensional discrete element method for granularmaterials. International Journal for Numerical andAnalytical Methods in Geomechanics, No. 14,
pp.451-472, 1990.
11）Seung Jae Lee, Youssef M.A. Hashash, iDEM: An
impulse-based discrete element method for fast
granular dynamics, International Jpurnal for
Numerical Methods in Engineering, Vol. 104, Issue2, pp.79-103, 2015.
12）Itasca Consulting Group, Inc., 3DEC 7.0
documenetation, 2019.
13）Nicolin Govender, Daniel N. Wilke, Schalk Kok:
Blaze-DEMGPU: Modular high performance DEM
framework for the GPU architecture, SoftwareX,
Vol. 5, pp.62-66, 2016.
14）V. Smilauer et al.: Yade Documentation 3rd ed.,2021. DOI:10.5281/zenodo.5705394 (最 終 閲 覧
日：2022年 7 月 1 日).
15）坂口秀・尾崎叡司・五十嵐徹：円形要素法を用いた DEM における回転の抑制に関する研究，
神戸大学農学部研究報告，Vol. 20（2），pp.239246，1993．
16）C. M. Wensrich, A. Katterfeld：Rolling friction as
technique for modelling particle shape in DEM,
Powder Technology, Vol. 217, pp.409-417, 2012.
17）Jun Ai, Jian-Fei Chen, J. Michael Rotter, Jin Y. Ooi:
Assessment of rolling resistance models in discretelement simulations, Powder Technology, Vol.
206, pp.269-282, 2011.
18）Y. C. Zhou, B. D. Wright, R. Y. Yang, B. H. Xu, A.
B. Yu：Rolling friction in the dynamic simulation of
sandpile formation, Physica A, Vol. 269, pp.536553, 1999.
19）M. J. Jiang, H.S. Yu, D. Harris: A novel discretemodel for granular material incorporating rolling
resistance, Computers and Geotechnics, Vol. 32,
pp.340-357, 2005.
20）中瀬仁・岩本哲也・曹国強・田部井和人・坂口秀・松島亘志：個別要素法簡易モデルによる地
震起因斜面崩落土の堆積範囲評価法の提案，土木学会論文集 A1（構造・地震工学），Vol. 73，
pp.I_694-I_703，2017．
21）福元豊・坂口秀・村上章：地盤材料の破壊基準を表現するためのシンプルな個別要素モデル ,土木学会論文集 A2（応用力学），Vol. 67，No. 1，
pp.105-112，2011．

110

野畑・古川・清野：多面体要素の個別要素法による粒状体シミュレーションへの転がり摩擦の導入に関する基礎的研究

22）橘一光・森口周二・寺田賢二郎・高瀬慎介・京
谷孝史・加藤準治：個別要素法を用いた落石シ
ミュレーションにおける形状精度と解析精度の定量的関連付け，土木学会論文集 A2（応用力学），
Vol. 70，No. 2（応用力学論文集，Vol. 17），I_519I_530，2014．
23）東祐輔・古川愛子・清野純史：中詰め材の大きさに着目した石橋の耐震性に関する基礎的研究，
自然災害学科学，Vol. 40，特集号，pp.67-80，2021．
24）Johnson, K L.: Contact mechanics, 2nd edn. Cambridge
University Press, 1985.
25）オンライン分度器，https://www.ginifab.com/
feeds/angle_measurement/online_protractor.
ja.php（最終閲覧日：2022年 7 月 1 日）．
26）星野寔・根岸正充：岩石の力学的性質に関する要

一実験，応用地質，11巻， 3 号，pp.91-96，1970．
27）松島亘志：個別要素法の基礎と応用 Ⅶ．実際に個別要素法プログラムを動かしてみよう，
http://www.kz.tsukuba.ac.jp/~tmatsu/dem/seg_
program/2010第 3 回 DEM セミナー松島資料.pdf（最終閲覧日：2022年 7 月 1 日）．
28）中瀬仁・曹国強・田部井和人・栃木均・松島亘
志：個別要素法による原子力発電所周辺の地震起因性斜面崩落挙動のモデル化と適用性，土木学会論文集 A1（構造・地震工学），Vol. 71，No.
4（ 地 震 工 学 論 文 集 第34巻 ），pp.I_476-I_492，
2015．
（投 稿 受 理：2022年 4 月 1 日
訂正稿受理：2022年 7 月 4 日）

旨

個別要素法を用いた粒状体シミュレーションでは，円形要素が用いられることが多い。円形
要素は回転し易いため，要素の回転を抑制する方法として転がり摩擦が用いられてきた。一方，多面体要素を用いた個別要素法では，要素形状自体に回転を抑制する効果があるため，転がり摩擦を導入した解析事例はない。しかし，要素形状が実際の粒状体の形状と異なる場合，転がり抵抗も異なる。そのため，解析結果も実挙動と異なるものとなってしまう。そこで本研究では，多面体要素を用いた個別要素法に新たに転がり摩擦を導入した。提案した転がり摩擦が個別要素法において機能するか，転がり摩擦の導入により形状の異なる要素を用いた解析結果の差異を低減できるかどうか検証を行った。具体的には，十面体と八面体の 2 種類の形状の要素について滑り摩擦係数と転がり摩擦係数を変化させた安息角を求めるシミュレーションを行い，転がり摩擦が安息角に及ぼす影響を調べた。また，岩塊の崩壊シミュレーションも行い，滑り
摩擦および転がり摩擦が解析結果に及ぼす影響を検証した。