・関 暢晴

応用力学論文集Vol.8 (2005年8月） 土木学会

複雑形状を有する粒状体材料の画像データに基づく個別要素モデリング

Image-based distinct Element Modeling of Granular Material with Irregular Shapes

山田貴博

・関 暢晴

・白坂 一

・松井和己

・塩見忠彦

Takahiro YAMADA, Nobuharu SEKI, Hajime SHIRASAKA, Kazumi MATSUI and Tadahiko SHIOMI

∗正会員 学博 横浜国立大学 大学院環境情報研究院（〒240-8501横浜市保土ケ谷区常盤台79-7）

∗∗工修 日揮株式会社（〒220-6001横浜市西区みなとみらい2-3-1）

∗∗∗工修 富士写真フイルム株式会社（〒106-8620東京都港区西麻布二丁目26番30号）

†正会員 工博 横浜国立大学 大学院環境情報研究院（〒240-8501横浜市保土ケ谷区常盤台79-7）

††正会員 Ph.D.（株）竹中工務店 技術研究所（〒270-1395千葉県印西市大塚1-5-1）

Mechanical behaviors of granular materials, which are developed from the interaction of grains, are very complicated to understand. In such a field, the computational mechanics is a powerful tool to investigate various features by the simulation in the grain level. To model the behavior of granular materials, geometrical structures should be represented precisely since stability of the material is secured by contact between grains. In this paper, we present an image-based modeling based on the distinct element method for granular materials. The 3-D image data of a gravel specimen obtained by CT scan are utilized to represent a real geometrical structure and individual grains are separated in the geometrical model by an image processing technique. The distinct element method that deals with grains with irregular shape obtained by the image data is also presented.

Key Words: granular material, image-based modeling, distinct element method, irregular shape

1. はじめに

工学のさまざまな分野において，粉体や粒状体等の 個々の要素の集合体に対する力学的挙動を理解するこ とが必要となる．このような粒状体等においては，多 くの粒子，要素が集まったマクロ的な(平均的な)力学 挙動は複雑な非線形性を示す．したがって，このよう な現象を理解するためには，個々の粒子，要素のレベ ルの微視的な現象をとらえることが有用であり，コン ピュータを用いた数値シミュレーションはこのような 場面において有効である．特に地盤工学の分野におい ては，礫のような粒状体の力学特性が重要であり，や はりその非線形挙動は粒状体の微視構造における現象 に起因するものとして考える必要がある．また，礫な どは自然物によって構成される粒状体であり，その粒 子形状と配置からなる複雑な３次元の幾何学構造を有 することから，力学挙動の数値シミュレーションを行 う際には，幾何学的構造を正確に表現することが重要 であると考えられる．

このような中でHollister等¹⁾は，従来の手法では数 値解析モデルの作成が困難であるような非常に複雑な 幾何学的構造をもつ対象に対して，CTスキャナなど の非破壊検査装置から得たディジタル画像を直接に用 いたモデル作成の概念（イメージベーストモデリング）

を提案している．そこで，本研究では粒状体材料に対

するイメージベーストモデリングとして，粒状体材料 の３次元画像データの取得と画像データに基づく幾何 学モデルの作成およびそのモデルに基づく非線形挙動 の数値シミュレーションを考える．

筆者等²⁾はコンクリート材料を複雑な３次元構造を 持つ２相複合材料と捉え，イメージベーストモデリン グを適用し，異種材料界面における破壊を考慮した材 料の非線形挙動の数値シミュレーション手法を提案し ている．この研究では，複合材料内部の界面における 剥離，滑りのモデル化について検討しているが，基材 で固定された複合材料内部における界面の破壊を，直 交格子を用いた有限要素法において表現できるものに 留まっている．礫のような粒状体は，接触状態により 形態が維持され，変形による非線形性を呈する際には，

剥離，滑りによる大きな変位が現れる．したがって，筆 者らの提案する直交格子を用いた有限要素法に基づく コンクリート材料に対するイメージベーストモデリン グを粒状体の問題に適用することは困難である．

複雑な形状を有する粒子間の接触，滑りを考慮する 力学的なモデル化の手法としては，接触を考慮した有 限要素法が考えられる．粒状体における個々の構成要 素を分離し，それぞれを有限要素分割できればこのよ うな手法を適用することも可能である．しかしながら，

多数の粒子に対してそれぞれを有限要素法で分割し，粒 子間の接触を取り扱うのは計算時間の面から現実的に

は困難である．また，粒子の内部の応力状態は，粒状体 材料の挙動に対して一般に重要ではないことから，粒 子内部の応力を評価する有限要素法による計算は効率 的であるとはいえない．

これに対して，多数の粒子における接触，滑りを考慮 することで，粒状体などの挙動をとらえる手法として 古くから個別要素法が用いられている．個別要素法は，

球形粒子に対する手法として出発したものであるが近 年では，複雑な形状を有する粒子を取り扱うことが可 能な個別要素法がいくつか提案されている3),4),5),6),7)． そこで本研究では，これらの研究をふまえ，画像から 得られる詳細な幾何学的構造を表現可能な複雑形状を 有する粒子に対する個別要素法に基づく計算手法を提 案する．

以上のように本研究は，粒状体材料に対する個別要 素法に基づくイメージベーストモデリング手法を提案 するものである．

2. 礫の３次元画像

本研究では，粒状体材料として礫を取り上げ，その 画像に基づくモデル化を考える．用いた礫のモデル供 試体は，１辺が5cmの立方体容器に礫を詰め，エポキ シ樹脂により固定することにより作成したものである．

この供試体に対して，幅方向，奥行き方向が355画素 となるX線CTスキャンを高さ方向に連続的に行い，

高さ方向200画素となる３次元画像データを作成した

（図−1参照）．得られた３次元画像の画素の間隔は，

幅方向，奥行き方向が0.14mm，高さ方向が0.25mmで あり，また，画像はCT値により256階調となってい る．この画像データ全体から数値計算のためのモデル を作成したとすると，非常に大規模な計算となってし まうことから，本研究では，第１段階として，この画 像データの１部である図−2に示す50×50×50画素 を用いて幾何学モデルを作成し，数値計算を行う．

3. 幾何学モデルの生成

本研究では，まず上述の礫の画像に基づき，個々の 粒子を分離した幾何学モデルを生成することを考える．

３次元画像から非構造格子である有限要素分割を生 成する研究については，血管流解析や骨の応力解析等 のバイオメカニクスの分野で多く行われている．しか しながら，このような分野における手法は，基本的に １つの連続的な領域を生成することが目的であるため，

本研究で取り扱う粒状体のように，個別の粒子構成要 素が分離されたモデルの生成は考えられてこなかった．

また，複雑な形状を有する粒状体では初期配置も現象 を左右する重要な材料固有の幾何学的構造であり，画 像データから実際の配置の状態を取得することが必要 であると考える．以下では，粒状体の画像データから

5cm 355 pixel 5cm

355 pixel 5cm

200 pixel

図−1 礫供試体とCT画像

図−2 礫の３次元画像

個々の粒子を分離し，それらの形状と配置に関する幾 何学モデルを同時に生成する手法について述べる．

3.1 等値面によるモデリング

３次元画像から幾何学モデルを生成する場合には，濃 淡ディジタル画像から等濃度面を抽出し境界とみなす 手法が広く用いられている⁸⁾．図−3は，CT値による 濃淡画像の等値面を表したものである．CT値は，対応 する画素位置をサンプリング点として，その位置の平 均化された密度に対応する情報であることから，その 濃淡による等値面は画素間隔以上の解像度をもつ幾何 学的情報と考えられる．その結果，図−3の等値面は，

滑らかな曲面を構成したものとなっている．

しかしながら，粒状体材料のように接触した粒子が 多数存在する画像を対象とした場合，生成された等値 面においては接触した粒子は同一の粒子として認識さ れる．特に粒子の接触面における剥離，滑りを考慮し た数値シミュレーションを考える場合には，モデル化

においてあらかじめ粒子を分離し，接触面を同定して おくことが必要となる．したがって，本研究では上述 の濃淡画像の等値面に粒子の分離情報を付加すること により，幾何学的な精度を低下させることなく，粒子 を分離した幾何学モデルを生成することを考える．

3.2 画像処理による粒子の分離

２次元画像処理においては，二値化処理，距離変換，

ラベル付け，エッジ検出，等高線処理等の粒子の分離 あるいは領域分割を行う手法が用いられている．本研 究で取り扱う礫の３次元画像について，２次元画像処 理による粒子の分離を試みたが，礫の画像は非常に複 雑な３次元構造を表しているものであり，３次元構造 内の粒子の分離は成功しなかった．

そこで本研究では，三次元画像処理によって粒子間 に対して境界面を作成し，連結する粉状体の分離及び 個別認識手法を考える．２次元および３次元画像に対 して，多数の個体を認識し，分離を行うことのできる 手法としてWatershed分割法⁹⁾が知られており，さま ざまな分野の画像処理で用いられている．そこで本研 究においても，Watershed法を粒子の分離に適用する ことを考える．

Watershed法を適用するにあたり，画像に対して前

処理としての変換を行う．具体的には，まずしきい値処 理を行い，画像を２値化する．得られた２値化画像か ら領域部分と境界部分を濃度値の勾配として表現した 画像を生成する骨格化処理を三次元ユークリッド距離 変換を用いて行う．さらに，距離変換処理によって生じ る濃度値の変動を抑制すべく，３次元平滑化処理を行 う．以上のように変換された画像に対して，Watershed 法を用いる．

Watershed法ではマーカと呼ばれる領域の核を濃度

勾配の極小値をとる部分に設定し，隣接画素へ勾配に 沿って広げていく．異なるマーカを拡大させた結果と して接触する部分が現れ，これを異なる領域の境界面 として領域分割をおこなう．このとき，濃度勾配の大 きい画像のエッジ部分に境界を作成されることとなる．

以上の処理で得られたラベル付けされた画素に対する 断面画像を図−4に示す．

一方，Watershed法では，領域の分割がしばしば過 分割になることが知られている．これは，Watershed 法では，濃度の勾配を領域の分割の指標としているた め，粒子の形状にくびれのような部分があった場合に は，粒子が分離されたものとして認識されることとなっ てしまうためである．図−4においても，本来は同一 の粒子と考える領域ものが複数の領域に分割されて認 識されている部分を見いだすことができる．したがっ て，Watershed法における過分割を何らかの手法で抑 制する必要性がある．

本研究では先に述べた過分割の抑制に対して，礫の ような粒状体において各粒子の接触面は粒子の表面全

図−3 ３次元画像から得られる等値面

体と比較すれば小さいという対象が有する物理的な先 験情報を利用する．過分割された領域同士は，大きな接 触面を持っているものと考え，そのような領域を再結 合することにより，過分割を除去する．具体的には，分 割された領域間の接触面の大きさを評価し，領域の接 触面の法線方向の投影面積と比べて十分大きい場合に は過分割されたものと判定し，領域の再結合を行った．

いま，画素の大きさが幾何学モデルの基本的な精度 を定めていることから，粒子を表現する画素数のある 程度少ないものについては，１つの粒子としてモデル 化するために必要な精度を有していないものと考えら れる．したがって，本研究では，２０画素程度以上で 表される大きさ以上の粒子のみを考慮し，それより小 さい粒子を除去した．

図−4の画像に対して再結合処理を行った結果を図−

5に示す．この結果から，本研究で提案する礫のための 再結合処理は，過分割された領域を適切に結合してい るが，部分的には不自然に分離された粒子も残ってい る，しかしながら，かみ合ったような形で物理的に不 自然に分離された粒子は，後述する個別要素法による 計算においては，実際は一体となって運動することと なる．したがって，ここで用いた手法によって，完全 に物理的な実体と対応したものとは言えないが数値計 算上は妥当な領域分割を与えていると考える．

3.3 粒子毎モデルの生成

Watershed法を用いると，各画素にラベル付けがさ

れた図−6のような画像が得られる．この画像を元に，

各画素は空隙，粒子間境界および粒子領域に分類する ことができる．この分類に基づき濃淡画像において領 域分割を行うため，各粒子に対してその粒子の領域の みを抽出するため，分類された画素とその境界画素を １，それ以外を０とする図−7のような３次元ウイン ドウを作成する．このウインドウにおける値と，あら かじめ生成すべき等値面のレベルを０に調整した濃淡

画像の画素値の積をとることによって，ラベル付けら れた各粒子毎に分割した濃淡画像が生成できる．

得られた濃淡画像に対して等値面生成を行うことに より，粒子毎の表面パッチとしての図−8のような幾 何学モデルが作成される．このとき，粒子間の境界の 画素を０とし，０を等値面の生成レベルとすることに より，隣接粒子に対してそれぞれの表面パッチは境界 が接する形でモデル化されることとなる．

以上より得られた全体の幾何学モデルを図−９に示 す．なお，ここまでの画像処理にはMatlab¹⁰⁾を用いた．

図−4 Watershed法による粒子の分離

図−5 再結合された粒子の領域

5 5 5 5 5 5 5 5 5

5

1 1 1 1 1

8 8

8 8 8 8 8 8 8 8 8 8 8 8

8 8 8 8 8 8 8 8

8 8

8 8 8 8 8 8 1

5

0 0 0 0 0 0

0 0 0 0

0 0

0 0 0

0 0

0=境界 1=空隙

5,8=ラベル付けられた領域

図−6 Watershed法により分割された画素の例

4. 複雑形状を有する要素を用いた個別要素 法

4.1 形状表現と接触判定

上述の幾何学モデルにおいては，三角形パッチによ り表面が表された複雑な形状を有する粒子の集合とし て粒状体材料がモデル化される．そこで，本研究では 表面パッチとして表されている幾何学モデルに対して 個別要素法を適用する．

複雑形状に対する個別要素法として，間隙ばねで 結んだ複数の球形個別要素によって表現する手法が提案

図−7 領域分割のためのウインドウ

図−8 各粒子の幾何学モデル

図−9 全体の幾何学モデル

されている^3),4)．また，Munjizaらは多面体などに対し て有限要素法による接触解析で用いられている計算法 を用いた手法を提案している^5),6)．しかしながら，前節 で得られたような複雑な形状を有する粒子を多数考え た場合の接触判定や接触力の評価を直接行うことを想 定すると，これらのアプローチでは非常に大きな計算 量を要することから現実的ではないと考える．

一方，Matsushimaらは比較的少ない球要素の集合 体として複雑形状を近似的に表す手法⁷⁾を提案してい る．本研究では画像データから得られた粒子の初期配 置を考慮した数値計算を考えていることから，初期配 置における釣り合いを満足させるための微調整などを 行うため直接的に球要素の配置を設定する手法が必要 と考える．

そこで本研究では，図−10のように等値面として生 成された三角形パッチの頂点が中心なる一定の大きさ の球状要素を配置し，この球状要素によって接触判定 と接触力の評価を行う手法を採用する．このとき，１ つの粒子は構成する球状要素を剛に結合した剛体と考 える．剛体としての粒子の運動は，球状要素から計算 される接触力の合力を外力とした運動方程式によって 記述されるものとする．このように，複雑形状を球要 素で表すことで，接触判定を高速に行うことが可能と なり，多数の複雑形状を有する粒子から構成される粒 状体材料の数値計算が可能となる．実際に，３次元画像 に対して生成された個別要素モデルを図−11に示す．

本研究で用いている手法は，有限要素法における接 触判定と接触力の評価に対して提案されたピンボール アルゴリズム¹¹⁾と類似な手法といえる．ピンボールア ルゴリズムによる接触力評価では，滑らかな面におけ る接触状態において不適切な摩擦力が評価されること が指摘されている⁶⁾．しかしながら本研究で取り扱って いる礫のような粒状体においては，広い面で接触が起 きることは稀であると考えられることから，粒子の表 面に配置した球要素による接触判定，摩擦力評価を採 用した．また，本研究で用いたモデルでは，ピンボー ルアルゴリズムで見られるような重なりの少ない表面 における球要素の配置¹¹⁾⁶⁾ではなく，重なりを多くし たことで滑らかな面における接触力の問題を緩和して いる．具体的には，画素の大きさの1.5 倍程度の半径 を持つ球要素が画素の大きさの0.5〜1.0倍の間隔で表 面に配置されたものとなっている．

なお，球要素の接触判定に対しては，空間領域をセ ルに分割し接触探索領域を対象要素が含まれるセルと その近傍のセルに含まれるものに限定することで効率 化をはかる近傍セル法を用いた．

本研究で提案した粒状体の個別要素モデル作成の流 れを図−12に示す．

図−10 表面パッチと個別要素モデル

図−11 複雑形状に対する個別要素モデル

4.2 接触力評価

接触した球要素間に対しては，通常の個別要素法と 同様な相互作用を考える．具体的には，接触面の法線 方向すなわち接触した２つの要素の中心を結ぶ直線の 方向とそれに対して垂直となる接触面の接線方向に接 触バネとダッシュポットを考える．さらに，接触面の 接線方向に対してはモール・クーロンの破壊基準に従 う摩擦スライダを考える．なお，接触バネのバネ定数 はヘルツの定理を用いて粒子のヤング率と重心から球 要素中心までの距離から決定し，ダッシュポットの減 衰定数は臨界減衰となる定数を与えた．

4.3 数値計算法

上述のように粒子は同一半径の球要素が表面に配置 され，それらが剛に結合されたものとして表されてい る．本研究では，粒子における球要素の配置を記述す るため，粒子毎に粒子の重心を原点とする局所座標系 を設定し，局所座標系において球要素の中心座標を与 える．粒子の質量，重心，重心周りの回転慣性は，画 像処理の段階で２値化により得られるボクセル領域か ら計算する．

運動方程式に対しては，局所座標系において記述さ れる並進３成分および回転３成分の６つの方程式を中 心差分法により離散化する．運動方程式の外力として

は球状要素から計算される接触力の合力を与える．全 体座標系における粒子の位置と姿勢は，局所座標系上 で得られた速度と角速度を全体座標系に座標変換し，積 分することにより得られる．なお，角速度を積分する 際には，四元数を用いることにより，任意の回転量が 精度の低下無く評価できる手法を用いた．

5. 数値計算例

数値計算は，画像から得られた粒子の個別要素モデ ルに対して，６面をすべて剛体壁とした直方体の計算 領域を設定し，剛体壁の変位および剛体壁に対する合 力によって境界条件を与えた．このとき，粒子とは壁 面に対しては粒子間と同様な接触条件を設定した．数 値計算において用いた物性値を表１に示す．また用い たモデルの粒子数は 34 ，表面には位置された球要素 の総数は6552である．

5.1 初期配置の作成

従来の個別要素法を用いた計算においては，初期配 置は粒子の詰め込みを与えられた領域に対して行う計 算によるものが多かった．一方，本研究では，詰め込ま れた状態の画像から粒子の形状データと初期配置を得 るため詰め込みの計算は行わない．しかしながら，画像 データから完全に釣り合った接触状態にある粒子の初 期配置を作成することは困難である．また，計算対象と なる粒子群に対して設定される境界条件も画像データ の元となった粒状体におけるものとは異なることから，

得られた個別要素による計算モデルは釣り合い状態に はない．そこで，本研究では画像データにおける配置に 近い釣り合い状態にある初期配置を得るため，最も簡 便な方法として落下法により初期配置の計算を行った．

落下法による初期配置計算では，まず画像から得ら れた粒子配置から粒子間で接触が起こらない程度に鉛 直方向にずらした粒子配置を生成する．その状態を初 期条件として，剛体壁で囲まれた計算領域において重 力による落下の計算を行う．図−13に落下法による初 期配置計算における初期条件と計算された釣り合い状 態にある初期配置を示す．ここでは，より短時間の計 算で釣り合い状態に収束させるため，粒子群を上下に ２グループに分け，２段階の落下計算を行った．

図−9の元となった画像データにおける配置と図−

13の得られた初期配置を比較すると，用いている画像 が本来の釣り合ったすべての粒子に対するデータでは なくその一部を切り取ったものであるため，特に境界 近傍に位置する比較的小さな粒子において大きく変位 し，元の画像とは異なった位置となっているものが見 られるが，おおむね画像に対応した初期配置が得られ ている．したがって，本研究で提案する手法では，初 期配置を含めた粒状体材料の計算モデルが得られてい ると考える．

図−12 モデル生成の流れ

表１ 物性値

ヤング率 37.80 GPa

ポアソン比 0.29 要素間の摩擦係数 0.30 壁面との摩擦係数 0.25

単位体積重量 2.65 g/cm³

図−13 落下法による初期配置の作成

5.2 3軸圧縮問題

計算手法に対する基本的な妥当性を検討するため，

図−14のような３軸圧縮問題の計算を行う．境界条件 としては，下面と２つの直交する側面を拘束する．側 圧に関しては，拘束していない側面に対して一定の力 を載荷し，水平方向の平均的な軸応力が一定である条 件を課す．さらに載荷は変形制御とし，上面を一定速 度で平行に移動させることで圧縮力を課す．

図−15に計算された粒状体の変形状態，図−１６に 応力とひずみの関係を示す．この結果においては，接 触した要素の滑りが荷重の増加に伴って発生し，段階 的に変形が進む様子が観測されている．

図−1４ 3軸圧縮問題

t = 0.0 s

t = 1.5 s

t = 3.0 s

図−15 3軸圧縮問題における粒状体の変形

5.3 せん断問題

粒状体材料に特有な挙動について検討を行うため，

図−17のようなせん断問題の計算を行う．載荷は変形 制御とし，直方体の計算領域に対してせん断変形を与 えるよう，側面を傾斜させる．一方，下面は固定し，上 面は平均軸応力を０とした条件として，自由とする．

図−18に粒状体の変形状態を示す．粒子の回転によ り，接触と分離が繰り返される複雑な挙動が見られる．

また，図−19に解析領域に対して評価した体積ひずみ とせん断ひずみの関係を示す．この図で見られる特性 は，実験で観測されるせん断変形により体積ひずみが 引き起こされるダイレイタンシ特性に定性的に対応し たものであると考えられる．

以上の計算はPentium 4 3.6GHzのPCを用い，２ 時間程度の計算時間を要した．

6. 結び

粒状体材料の３次元画像データに基づき，粒子の複 雑な形状と配置を考慮した幾何学モデルを構成する手 法と，得られた幾何学モデルに従った個別要素法に基づ く粒状体材料のシミュレーション手法を提案した．提案 した手法は，その初期配置を考慮できる手法であるこ とから，土質材料のように初期配置が材料の力学挙動 に大きく影響を与える問題に対して従来の個別要素法 に基づく手法では得ることが困難な知見を数値シミュ レーションで得ることが可能になると考える．

本論文で行った数値計算は，研究の第１段階とし て，画像の一部分のみを用いた少ない粒子のモデルにと

0 20 40 60 80

0 10 20 30 40 50

stress (MPa)

strain (%)

図−16 3軸圧縮問題における応力とひずみの関係

図−17 せん断問題

どまっている．また，接触バネの定数などパラメータ 設定の検討についても十分とは言えず，得られた計算 結果を定量的に考察する段階には至っていない．今後 は，より多くの粒子群による大規模な問題で提案した 手法の数値特性を把握するとともに，粒状体材料の定 量的特性評価に対する本手法の妥当性と有効性を検討 することが必要である．

t = 0.0 s

t = 1.5 s

t = 3.0 s

図−18 せん断問題における粒状体の変形

-1.0 -0.5 0 0.5 1.0

0 5 10 15 20 25 30

volumetric strain (%)

shearing strain (%)

図−19 せん断問題におけるせん断ひずみと 体積ひずみの関係

参考文献

1) Hollister, S.J. and Kikuchi, N.: Homogenization theory and digital imaging a basis for studying the mechanics and design principles of bone tis- sue, Biotech. Bioeng, Vol. 43, pp.586–596, 1994.

2) 永井学志,山田貴博,和田章: コンクリート材料の 圧縮破壊のイメージベース有限要素解析, コンク リート工学年次論文報告集, Vol. 21, pp.157–162, 1999.

3) Meguro, K. , Iwashita, K. and Hakuno, M. : Fracture Analyses of Media Composed of Irreg- ularly Shaped Resions by the Extended Distinct Element Method, Structural Eng. / Earthquake Eng., JSCE, Vol.8, pp.131–142, 1991.

4) 栗田哲史, 中瀬仁, 安中正, 嶋田昌義, 藤谷昌弘:

ケーソン式護岸の模型振動実験に対する個別要素 法シミュレーション,土木学会論文集, 525/I-33号, pp.299-308, 1995.

5) Munjiza, A. , Latham, J. P. and John, N. W. M.:

3D dynamics of discrete element systems com- prising irregular discrete elements – integration solution for finite rotations in 3D, Int. J. Numer.

Meths. Engng.,Vol. 56, pp.35–55, 2003.

6) Munjiza, A.: The Combined Finite-Discrete Ele- ment Method, John Wiley Sons Inc, 2004.

7) Matsushima, T. and Saomoto, H.: Discrete element modeling for irregularly-shaped sand grains, Proc. NUMGE2002: Numerical Methods in Geotechnical Engineering, Mestat ed., pp.239–

246, 2002.ls

8) Lorensen, W. E. and Cline. H. E. : Marching Cubes: A High Resolution 3D Surface Construc- tion Algorithm, Computer Graphics, 21, pp.163–

169, 1987.

9) Vincent, L. and Soille, P.: Watersheds in digi- tal spaces : an efficient algorithm based on im- mersion simulation, IEEE transactions on pat- tern analysis and machine intelligence, Vol. 13, pp.583–598, 1991.

10) http://www.mathworks.com/

11) Belytschko,T. and Neal, MO, Contact-impact by the Pinball Algorithm with Penalty and La- grangian Methods, Int. J. Num. Meths Engng., Vol. 31, pp.547-572, 1991.

（2005年4月15日受付）