複雑な砂粒子形状の個別要素モデル化手法の提案

松島・竿本：GIT2001 北関東地盤工学研究発表会
複雑な砂粒子形状の個別要素モデル化手法の提案
Grain-Shape Modeling of Real Sands and Rocks for Discrete Element Method
筑波大学機能工学系
松島亘志
筑波大学大学院
竿本英貴
1. はじめに
近年のコンピュータ性能の目覚しい向上により、地盤材料・紛体材料などの粒状体の、集合体としての力学挙
動を、構成粒子個々の運動からシミュレートする数値解析法が注目されてきている。地盤工学においては、接触
点にばね・ダッシュポット・スライダーを仮定した DEM( Discrete Element Method, 個別要素法)が多く用いられ
ており、特に連続体として解析する場合の構成則の形やその材料定数を定めるための要素試験については、実際
の境界条件をそのまま模擬した「数値要素試験」が実現しつつある。この数値要素試験によれば、実実験で計測
困難な諸量(特に粒子運動や接触点力などの微視的な物理量)が計測可能であり、また境界条件等の自由度も高い
ため、実験と相補的に用いることで、極めて多くの新知見を得ることが期待されている。
しかしながら、現在までの解析例のほとんどは球形(２次元解析では円形)要素を用いたものであり、それ以外
も２次元では多角形や楕円、３次元では楕円体を用いたものがある程度で、実際の砂粒子の形状を忠実に再現し
ようと試みた例はない。その原因としては、①実際の砂粒子の形状が非常に複雑であり、モデル化が困難である
こと、②砂粒子の形状を精密に計測した研究例が少ないこと(特に３次元形状計測例は皆無)、③まだコンピュー
タ性能が不足しているという一般認識があること、などが挙げられる。一方で、粒子形状がマクロな力学挙動に
大きく影響を及ぼすことは、古くから実験などで確かめられており、現存する多くの実実験結果と直接比較する
ためには、粒子形状の適切なモデル化が必要であることは確かである。
そこで本研究では、複雑な砂粒子形状を任意の精度でモデル化し、それを DEM 解析に用いる手法を提案する。
DEM で自然粒子の複雑な形状を表現する手法はいくつか考えられるが、ここでは①３次元への拡張性、②多く
の汎用ソフトへの適用性、を考慮して球形要素(2D では円形要素)の剛接続によって表現する方法を考える。
[0] Start
[1] Grain data input
grain surface data of a sand grain
解くべき問題は「なるべく少ない数の円(３次元の
[2] Calculation of volume (area), gravity center, etc. of the grain
場合は球)で、所用の精度を有する粒子を表現する [3] Input the calculation conditions
number of elements and their initial configuration
ための円(球)の位置および半径を定める」ことで
time increment, spring constant, damping coefficient
ある。基本的には、粒子表面のうち、曲率の大き
[4] Iteration loop start
な部分は小さな円(球)で表現し、曲率の小さな部
[5] Loop for each grain-surface data point
[6] Detection of the element to describe the surface point
分は大きな円(球で)表現するのが効率的であるこ
[7] Calculation of force from the distance
とは察しがつくが、ここでの「曲率」はある範囲
[8] Return [5] up to the end of data point
[9] Solve the equation of motion for each element
を持った表面の平均的な曲率の意味となり、その
and find the next configuration
「範囲」も変数となる（すなわち、ある円(球)が
[10] Calculation of accuracy
粒子表面のどの領域をカバーするかということも、[11] If accuracy is not enough, go to [4]
円(球)の相対位置や半径で変わってくる）ため、 [12] End
2. 複雑な粒子形状のモデル化手法
問題はそれほど単純ではない。
図1
球要素による砂の最適形状表現アルゴリズム
2-1. 計算アルゴリズム
本研究では次のような動的反復計算アルゴリズムを提案する。すなわち、粒子形状を表現するための球(円)
の半径及び中心位置が、運動方程式に従って時間変化し、結果的に最も精度の良い解に到達するように、
-1-
松島・竿本：GIT2001 北関東地盤工学研究発表会
球(円)に働く仮想力を与える。この仮想力は、球要素と粒子表面の変位に比例し、粒子表面に引き寄せられ
る方向に働くものとする。仮想力は、ある表面のデータ点と、その表面を表現する球要素の表面の間にの
み働く。そして、粒子の運動に対して、減衰を考慮することで、球要素は振動しながら収束解に近づいて
ゆく。図 1 は、その計算のフローチャートを示している。
2-2. ２次元計算
以降は、図示による理解の容易な２次元問題について具体的な計算例を説明する。図 2 は、豊浦砂の砂粒
子形状の例である。これを複数の円要素の剛接によって、なるべく効率的にモデル化することが目標であ
る。まず、１個の円要素で計算した結果を図 3 に示す。はじめ、小さな半径を持っていた円要素は、砂の
表面から力を受けることで膨張し、またその中心位置も移動する。それらは若干振動しながら平衡点へと
移動し、10 step 以降は、ほとんど変化がなくなる。収束の仕方はばね定数と粘性係数によって変化するが、
半径の変化については、余り振動しない方が良い収束を得られるようである。図 4 は、random に発生させ
た異なる初期位置(ただし粒子輪郭内)から始めた場合の計算結果であるが、全ての場合で同じ収束解が得ら
れている。
このときのモデル化の精度を表現するのに、次のような形状誤差指標を用いる。
N
err =
∑ d k − rk
k =1
N Rav
ここに、 N は輪郭データ数、 Rav は粒子面積と等価な面積と有する円の半径、 d k 、 rk は、 k 番目のデータ
を表現する球(円)要素の中心から輪郭データまでの距離、および要素半径である。１要素の円による近似で
は、誤差は 9.04%となっている。これは、粒子半径に対して、平均として 9%程度の半径誤差を含んでいる
ことを意味する。
図 5 は、要素２つによる表現の例である。10 個の異なる初期位置の組み合わせに対して、異なるいくつか
の収束解が得られた。このように、複数の要素について計算する場合は、初期位置によって異なった解に
収束するので、random に初期位置を変化させたいくつかのケースについて計算を行い、最も精度の良い収
束解を選択することになる。ここでは、10 回の計算で、最も精度の良かった２例について示しているが、
最良解の誤差は 5.48%であった。
要素数を多くすると、新たな問題が生じる。図 6 は 3 要素による計算結果であるが、ある場合は、一つの
要素が、別の要素の内部に入り込んでしまい、機能していない状態が発生する。これは、アルゴリズムか
ら見れば当然で、輪郭点からの引力はその輪郭点を表現する粒子のみにしか働かないために、一旦このよ
うな状態になると、内部の粒子には全く力が働かないからである。この問題は、粒子の発生の仕方を工夫
することで、かなりの程度回避することができる。図 7 は 10 要素による形状近似の例であるが、初期位置
を輪郭に沿って図のように配置することで、全ての要素が形状表現に寄与するような結果となっている。
更に、付加的なアルゴリズムとして、機能していない粒子は、その時点で最も精度の悪い部分に再配置す
ることとし、より高精度な解を得ることを目指している。
要素数が多くなると、収束に多くの時間を要するようになるという問題もある。これは、収束計算の途中
で、要素が移動することにより、それらに働く力の系が刻々変化するためである。しかし、精度自体は単
調減少し、ある収束解に向かう。この収束性の改善は今後の課題である。
2-3. 3 次元計算
図 8 は３次元粒子データを基に、10 要素による形状近似を行った例である。３次元の場合、10 要素では充
分な近似ができていないが、要素数を増やすことにより、精度の向上が可能である。なにより、全く同様
-2-
400
400
350
350
k=100, cr=0.2, cx=cy=0.5
350
300
300
250
250
250
200
150
150
250
300
350
400
450
100
200
500
step=1
initial
250
300
350
400
450
100
200
500
250
300
x
豊浦砂粒子の例
350
400
450
500
x
図 3 円要素１つによる形状近似計算
k=100, cr=0.2, cx=cy=0.5
k=100, cr=0.2, cx=cy=0.5
400
The resultant circle is independent
from its initial position.
case 4 & 7: err=0.0738730
case 5, 8, 9 & 10: err=0.054805
300
300
250
250
250
150
y
300
y
350
200
k=100, cr=0.2, cx=cy=0.5
400
350
y
350
200
150
x (pixel)
400
step=5
y
200
図2
step=10
step=20
step=3
300
100
200
k=100, cr=0.2, cx=cy=0.5
400
step=2
y
y (pixel)
松島・竿本：GIT2001 北関東地盤工学研究発表会
200
200
150
150
err=0.0904271
100
200
250
300
350
400
450
500
100
200
250
300
350
400
450
500
100
200
250
x
x
図 4 様々な初期位置からの収束計算結果
300
350
400
450
500
x
図 5 ２要素による形状近似
のアルゴリズムで、３次元のモデル化を行えることが本手法の特徴である。
3. おわりに
より高精度な個別要素法解析のための、粒子形状モデル化手法の提案を行った。現在、豊浦砂の２次元形
状データを基にした粒子モデルの作成、およびそれを用いた個別要素法解析を行っている最中である。基
本的には円形要素の個別要素法であるから、解析において困難な障害が発生することは考えにくい。ただ
し、一つの粒子中で、円要素がいくつも重なっているために、同じ食い込み量でも反力が２倍、３倍にな
る部分が存在するという問題は予想され、ばね定数や時間刻みの設定に注意が必要である。
例えば１つの粒子を 10 個の円要素でモデル化した場合、1000 個の粒子の解析を行うには、
実質的には 10000
要素の解析の問題となる。接触判定の合理化を行えば、計算時間はほぼ要素数に比例するから、1000 個の
円形粒子の計算に比べて約 10 倍の計算量となる。要素試験目的であれば、この程度の負荷は充分実現可能
であると考えられる。
参考文献
[1] Muhlhaus, H.-B., Sakaguchi, H. and Moresi, L.: Particle in Cell and Discrete Element Models for Granular
Materials, Computer Methods and Advances in Geomechanics, Desai et al. eds., Balkema, Vol.1, pp.511-518. 2001.
[2] Matsushima, T. Konagai, K: Grain-shape effect on Peak Strength of Granular Materials, Computer Methods and
Advances in Geomechanics, Proc. 10ACMAG, Desai et al. eds., Vol. 1, pp.361-366. 2001.1.
-3-
松島・竿本：GIT2001 北関東地盤工学研究発表会
k=100, cr=0.2, cx=cy=0.5
400
k=100, cr=0.2, cx=cy=0.5
400
case 10, err=0.0404938
case 3, err=0.0739351
case 1, 4-9, err=0.0474879
350
300
300
300
250
250
250
y
y
350
y
350
200
200
200
150
150
150
100
200
250
300
350
400
450
100
200
500
k=100, cr=0.2, cx=cy=0.5
400
250
300
350
x
400
450
100
200
500
250
300
x
x
図6
３要素による形状近似の例
400
400
Initial state
350
豊浦砂の粒子輪郭の例
10個の円要素による輪郭表現
350
300
300
y
y
250
250
200
200
150
150
100
200
250
300
350
400
450
500
100
200
x
図7
250
300
350
x
400
450
500
10 要素による形状近似の例（左：初期配置、右最終状態）
6
4
-2
Y Axis
2
0
-4
-6
420 14
1086 is
12
14
Z
16
Ax
18
20
X Axis
22
24
26
28
14
12
Z Axis
10
8
6
4
2
14
0
16
18
20
22
X Axis
図8
24
26
64
20
-228
３次元形状近似の例（上からおよび側面からの観察結果）
-4-
350
400
450
500

Download Report