粒子シミュレーション装置及び粒子シミュレーション方法
【課題】メモリへの書込みの競合を防ぐことができ、並列プロセッサを用いて離散要素法の演算効率を向上させる。
【解決手段】作業空間５１内の複数の粒子が他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置１０において、複数の粒子のそれぞれには粒子番号が付されており、作業空間５１が複数のセルに分割され、各セルにはセル番号が付されており、複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得し、取得されたセル番号に基づいて、複数の粒子の粒子番号を付け替える粒子番号変更部１４を備える。
【解決手段】作業空間５１内の複数の粒子が他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置１０において、複数の粒子のそれぞれには粒子番号が付されており、作業空間５１が複数のセルに分割され、各セルにはセル番号が付されており、複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得し、取得されたセル番号に基づいて、複数の粒子の粒子番号を付け替える粒子番号変更部１４を備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子シミュレーション装置及び粒子シミュレーション方法に関する。
【背景技術】
【０００２】
土砂崩れや地盤液状化などの自然災害予測や粉粒体現象の解析を高精度に行うためには、数値シミュレーションを用いて実験による観察が難しい微視的な知見を得ることが重要である。しかし多くの場合、固体粒子群は数えきれないほどの多数の粒子が多点で同時接触をしながら運動している。したがって、固体粒子群の運動を正確に予測するためには、物理に基づき、可能な限り多数の粒子を扱い、ときには長時間スケールを扱えるシミュレーションが必要となる。
【０００３】
固体粒子群の運動を扱うシミュレーション手法として、Ｖｏｉｇｔモデルを用いて粘弾性体としての粒子運動を再現するＤＥＭ（ＤｉｓｃｒｅｔｅＥｌｅｍｅｎｔ Ｍｅｔｈｏｄ：離散要素法）が知られている（例えば特許文献１）。ＤＥＭは、多体衝突および摩擦や粒子の回転運動が扱えるため、高濃度で多数の粒子を扱う粉体産業や自然科学の分野で広く用いられてきた。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００７−１７２１６９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかし、ＤＥＭ で扱える粒子数は現実よりも遥かに少なく、実現象を正確に再現するためにＤＥＭ に対する高速化アルゴリズムの開発が求められている。
【０００６】
これまでＤＥＭは、接触相手粒子の探索方法の効率化や並列化およびベクトル化を行うことで高速化されてきた。しかし、ＤＥＭにとって並列化やベクトル化は容易ではない。並列化ではノード毎に粒子数が異なる上にノード間で情報交換すべき粒子が不定であるために並列化率が向上し難い。また、ベクトル化ではメモリへの書き込み競合を回避する処理が多く必要となるためにその処理で計算時間を浪費する。そのため、ＣＰＵ（Central Processing Unit）がムーアの法則を破り桁違いに高性能化されない限り、ＤＥＭの飛躍的な高速化は期待できない。
【０００７】
一方、最近ではＣＰＵと比べて価格に対する演算性能の比が格段に大きいＧＰＵ（Graphics Processing Unit）を用いた高速化が盛んに行われている。しかし、粒子に働く力の全てを考慮している完全なＤＥＭに対してＧＰＵ を用いて高速化した例は無く、さらにＧＰＵは並列計算を行うため、メモリ書き込み競合を避けるべく計算量を増大させているためにＧＰＵの性能を十分に生かせていない。
【０００８】
本発明は、上記課題を解決するためになされたものであり、メモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる粒子シミュレーション装置及び粒子シミュレーション方法を提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記課題を解決するため、本発明に係る粒子シミュレーション装置は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、複数の粒子のそれぞれには粒子番号が付されており、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得手段と、セル番号取得手段により取得されたセル番号に基づいて、複数の粒子の前記粒子番号を付け替える粒子番号変更手段と、粒子番号変更手段により付け替えられた粒子番号に基づき、複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択手段と、近傍粒子選択手段により選択された近傍粒子と前記１つの粒子との位置関係に基づいて、１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算手段と、接触力計算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段とを備えることを特徴とする。
【００１０】
同様に、上記課題を解決するため、本発明に係る粒子シミュレーション方法は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、複数の粒子のそれぞれには粒子番号が付されており、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得ステップと、セル番号取得ステップにおいて取得されたセル番号に基づいて、複数の粒子の前記粒子番号を付け替える粒子番号変更ステップと、粒子番号変更ステップにおいて付け替えられた粒子番号に基づき、複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択ステップと、近傍粒子選択ステップにおいて選択された近傍粒子と前記１つの粒子との位置関係に基づいて、１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算ステップと、接触力計算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップとを備えることを特徴とする。
【００１１】
このような粒子シミュレーション装置及び粒子シミュレーション方法によれば、作業空間の複数の粒子をセルに割り当てるときに、各粒子が配置されるセルのセル番号が取得され、このセル番号に基づいて粒子の粒子番号が付け替えられるため、粒子をセルに割り当てる際にメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００１２】
上記課題を解決するため、本発明に係る粒子シミュレーション装置は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算手段と、接触力計算手段により計算された粒子間の接触力情報を配列に格納する接触力格納手段と、接触力格納手段により粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算手段と、総和演算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段と、を備えることを特徴とする。
【００１３】
同様に、上記課題を解決するため、本発明に係る粒子シミュレーション方法は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算ステップと、接触力計算ステップにおいて計算された粒子間の接触力情報を配列に格納する接触力格納ステップと、接触力格納ステップにおいて粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算ステップと、総和演算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップと、を備えることを特徴とする。
【００１４】
このような粒子シミュレーション装置及び粒子シミュレーション方法によれば、算出された粒子間の接触力情報が配列に格納され、この配列を参照して粒子ごとの接触力が総和演算で得られるため、接触力の総和演算時に発生するメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【発明の効果】
【００１５】
本発明に係る粒子シミュレーション装置及び粒子シミュレーション方法によれば、メモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【図面の簡単な説明】
【００１６】
【図１】本発明の第１実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図２】本実施形態で粒子モデルに適用するＶｏｉｇｔモデルを示す概略図である。
【図３】本実施形態における作業領域を示す図である。
【図４】粒子番号変更の流れを示す図である。
【図５】作業領域における粒子番号を変更する処理と、各セルに粒子番号の最大値及び最小値を記憶する処理を示す図である。
【図６】近傍粒子表の構成の一例を示す図である。
【図７】接触候補リストの構成の一例を示す図である。
【図８】現在の接触候補リストのペア粒子について１ステップ前にペアであったか調べる処理を示す図である。
【図９】接触力参照表の構成の一例を示す図である。
【図１０】接触力参照表の構成の一例を示す図である。
【図１１】本発明の第１実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【図１２】本発明の第２実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図１３】近傍粒子表の構成の一例を示す図である。
【図１４】本発明の第２実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【図１５】本発明の第３実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図１６】粒子をセルに格納する流れを示す図である。
【図１７】本発明の第３実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【発明を実施するための形態】
【００１７】
以下、図面を参照しながら本発明の実施形態を詳細に説明する。なお、図面の説明において同一又は同等の要素には同一の符号を付し、重複する説明を省略する。
【００１８】
（第1実施形態）
図１は、本発明の第１実施形態に係る粒子シミュレーション装置１０の機能ブロック図である。図１に示すように、粒子シミュレーション装置１０は、粒子情報保持部１１、粒子情報取得部１２、接触候補リスト更新判定部１３、粒子番号変更部（セル番号取得手段、粒子番号変更手段）１４、近傍粒子表作成部（近傍粒子選択手段）１５、接触候補リスト作成部（接触候補粒子選択手段）１６、接触力参照表作成部１７、接触判定部（接触判定手段、接触力演算手段、接触力計算手段、接触力格納手段）１８、接触力計算部（接触力演算手段、総和演算手段）１９、粒子情報更新部（粒子情報更新手段）２０を備えている。
【００１９】
粒子シミュレーション装置１０は、物理的には、ＧＰＵ（Graphics Processing Unit）、主記憶装置であるＲＡＭ及びＲＯＭ、ハードディスク装置等の補助記憶装置、入力デバイスである入力キー等の入力装置、ディスプレイ等の出力装置、ネットワークカード等のデータ送受信デバイスである通信モジュールなどを有するコンピュータとして構成されている。図１において説明した粒子シミュレーション装置１０の各機能は、ＧＰＵ、ＲＡＭ等のハードウェア上に所定のコンピュータソフトウェアを読み込ませることにより、ＧＰＵの制御のもとで入力装置、出力装置、通信モジュールを動作させるとともに、ＲＡＭや補助記憶装置におけるデータの読み出し及び書き込みを行うことで実現される。
【００２０】
まず本実施形態で適用する粒子モデルについて説明する。粒子間接触力の法線方向および接線方向成分Ｆｎ，Ｆｔは、図２に示すＶｏｉｇｔモデル７１を用いて以下の式により表される。
【数１】
ここでｘとｖは相対変位および相対速度ベクトル、μtは滑り摩擦係数、ｍｉｎは絶対値が小さい方の値を表す関数、Ｋは線形バネ７２における弾性係数、ηはダッシュポット７３の粘性減衰係数であり、反発係数eと次式の関連がある。
【数２】
【００２１】
したがって、ＤＥＭで扱う並進および重心を中心とする回転の運動方程式は以下の式で表される。
【数３】
ここで、ｕｐとωは粒子の速度および角速度、ｍｐとＩｐは粒子の質量および慣性モーメント、Ｒは粒子中心から接触点へ向かう位置ベクトルである。また、（５）式の右辺第二項は転がり摩擦力を表し、μｒは転がり摩擦係数、ｂは接触面幅である。（１）、（２）式のＶｏｉｇｔモデル７１は、後述する接触判定部１８で粒子間接触力を算出するのに用いられ、また（４）、（５）式の運動方程式は、粒子情報更新部２０で粒子の座標を算出するのに用いられる。
【００２２】
また、本実施形態において粒子が運動する領域である作業領域５１は、図３に示すように、一辺の大きさがDmax（1．0＋α）の立方体セルで区分される。ここで、Dmaxは最大粒子径、αは更新頻度を調整するパラメータであり例えばα＝０．２である。また、作業領域５１内の各セルにはセル番号cell_idが付されている。セル番号の付け方は後述する。
【００２３】
粒子情報保持部１１は、作業空間５１内の複数の粒子のそれぞれについて粒子情報を保持している。粒子情報は、座標、並進速度、回転速度、粒子半径を含む。座標とは、図３に示す三次元の作業空間５１における三次元座標（ｐｏｓ．ｘ，ｐｏｓ．ｙ，ｐｏｓ．ｚ）である。粒子情報は、後述する粒子情報更新部２０により更新される。
【００２４】
粒子情報取得部１２は、粒子情報保持部１１から粒子情報を取得する。粒子情報取得部１２は、粒子シミュレーションの開始時に粒子情報保持部１１から予め設定されている初期状態の粒子情報を取得し、シミュレーション動作中には、粒子情報保持部１１が更新される度に新たな粒子情報を粒子情報保持部１１から取得する。
【００２５】
接触候補リスト更新判定部１３は、後述する接触候補リスト５３を更新するか否かを判定する。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。また、粒子の積算移動距離は、例えば、粒子情報取得部１２が新しい粒子情報を粒子情報保持部１１から取得する度に、１ステップ前の粒子情報との座標データの差分からこの直近の１ステップにおける粒子の移動距離を算出し、これを粒子ごとに積算して求める。なお、粒子移動距離の積算値は、接触候補リスト５３の更新が行われるたびに初期化される。
【００２６】
粒子番号変更部１４は、粒子を所属するセル番号順に並べ、その順に粒子番号を付け直す。セル番号は図３に示すように一次元であり、例えば次式のように表される。
cell_id＝i.x＋id.x×i.y＋id.x×id.y×i.z （６）
i.m = int(pos.m／cellsize.m)
ここでcell_idはセル番号であり、i．mはｍ方向のセル位置であり、id．mはｍ方向のセル数である。また、pos．mはｍ軸方向の座標を表し、cellsize．mはm軸方向のセルサイズを表し、int()は実数を整数に変換する関数である。
【００２７】
図４に示すように、粒子に自身が所属するセル番号を記憶させる（図４（ａ））。次に、セル番号が昇順になるように粒子を並べる（図４（ｂ））。そして、セル番号順に粒子番号を付け直す（図４（ｃ））。
【００２８】
図４の流れを式で表すと以下のとおりである。
Pcell[p_id]=cell_id
B_pid[p_id]=p_id
Step1
Pcell={5,13,0,21,5,3,6} B_pid={0,1,2,3,4,5,6}
Step2
Pcell={0,3,5,5,6,13,21} B_pid={2,5,0,4,6,1,3}
Step3
old_p_id=B_pid[p_id] Val[p_id]=Val[old_p_id]
ここで、p_idは粒子番号であり、Valは粒子座標、並進速度、回転速度、半径である。Step1のPcellの中身を周知のソート手法（例えばバイトニックソート）を用いて昇順にし、それと同じようにB_pidの中身も移動させる。この時点でBの中身はソート前の粒子番号が記憶されている。Step3でB_pidを使い、ソート前の粒子情報をソート後の粒子情報に書き換える。
【００２９】
この粒子番号を変更する処理により、例えば図５（ａ）に示すように作業空間５１に配置された粒子の粒子番号が、図５（ｂ）に示すように、セル番号に沿った粒子番号に変更される。
【００３０】
近傍粒子表作成部１５は、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２を作成する。近傍粒子表作成部１５は、着目粒子の所属するセルおよびその隣接セルに存在する粒子のうち、自身より大きな粒子番号の粒子で粒子間距離がある値以下の粒子を近傍粒子とする。以下に詳細を説明する。
【００３１】
近傍粒子表作成部１５は、セル内の粒子の番号のうち、最小値Pnum_in_cellminと最大値Pnum_in_cellmaxをセルに記憶させる（図５（ｃ））。番号cell_idのセルに隣接するセルの番号nei_cell_idは次式で表される。
nei_cell_id[n]= cell_id + con[n] (n=0〜26) （７）
ここで、conはid．xとid．yによって決まる定数である。セル番号で粒子番号をソートしてあるため、各隣接セルにはPnum_in_cellmin[nei_cell_id]からPnum_in_cellmax[nei_cell_id]までの番号の粒子が連続して存在することがわかる。これらの粒子を近傍粒子と呼ぶ。
【００３２】
各粒子ごとに、図６に示すような近傍粒子を格納する箱nei_pn[i][box_id]をn個用意し、近傍粒子のうち粒子中心間距離がdis_lim以下で、自身より大きい番号（i=p_id）の近傍粒子を番号（box_id）が小さい方の箱から入れて行き、一方、小さい番号（i）の近傍粒子は番号（box_id）が大きい方の箱から順に入れる。dis_lim = (ri + rj)(1．0+α)であり、ri,rjは近傍粒子ペアi,jの粒子半径、αは接触候補の更新回数を調整するパラメータである。図６に示す一連の箱の集合を、作業空間５１内の全ての粒子ｉについて用意したものを、本実施形態では近傍粒子表５２という。
【００３３】
さらに、粒子番号ｉが自身より大きい粒子および小さい粒子の数n_jgi， n_jliをそれぞれ全粒子に対して求める。図６の例では、例えば５番の粒子（ｉ＝５）の近傍粒子の粒子番号が１，３，４，７，９番であったとすると、
nei_pn[5][0]=7
nei_pn[5][1]=9
nei_pn[5][n-3]=4
nei_pn[5][n-2]=3
nei_pn[5][n-1]=1
n_jgi[5]=2
n_jli[5]=3
となる。
【００３４】
接触候補リスト作成部１６は、まず、図７に示すように、n_jgiのプレフィックス和s_jgiを次式から求める。
【数４】
【００３５】
そして、各粒子iに対して、box_id を 0 〜 n_jgi [i ] − 1まで変化させて、次式を用いて、候補リスト番号Lnumおよび接触候補相手の粒子番号jを求め、接触候補リスト５３list_i，list_jを作る。
Lnum = s_jgi [i - 1] + box_id
j = nei_pn[i][box_id]
list_i[Lnum] = i
list_j[Lnum] = j
【００３６】
例えば、図７の例では、接触候補リスト５３の内容は、以下の表１のようになる。
【表１】
つまり、接触候補リスト５３は、作業空間５１の全粒子について、接触している可能性がある粒子ペアｌｉｓｔ＿ｉ，ｌｉｓｔ＿ｊに番号Ｌｎｕｍを付加する。
【００３７】
次に、１ステップ前の接触候補ペアBlist_i，Blist_jを探索し、現在の接触候補ペアlist_i, list_jが１ステップ前にもペアであったか調べる。ペアであったら、図２に示した粒子モデルにおける法線ベクトルと接線方向のバネによる力（バネの伸縮量）を新しい候補リスト番号の配列に格納する。
【００３８】
粒子番号i=list_i[Lnum], j=list_j[Lnum]の接触候補ペアが有るとき、このペア粒子の１ステップ前の粒子番号はBi =Bpn[i]， Bj = Bpn[j]である。Bj＞BiのときはBs_jgi[Bi−1]〜Bs_jgi[Bi]−1の範囲の接触候補Blist_j[BLnum]を調べBjと一致するか調べる。一致したときはそのリスト番号BLnumを用いて１ステップ前の法線ベクトルと接線方向の接触力BSpring[BLnum]を現在のリスト番号配列Spring[Lnum]に代入する。逆にBi＞Bjのとき、Bjに対して粒子番号が大きい粒子が接触候補として格納されているリスト番号BLnumの範囲はBs_jig[Bj−1]〜Bs_jgi[Bj]−1であるから、この範囲でBlist_j[BLnum]がBiと一致するか調べ、一致するときはSpring[Lnum]に−BSpring[BLnum]を代入する。
【００３９】
この処理について、図８を参照して、現在の接触候補リスト番号Lnum = 6におけるｉ＝５，ｊ＝９番のペア粒子について１ステップ前にペアであったか調べる処理を説明する。
【００４０】
まず、ｉ＝５，ｊ＝９の１ステップ前の粒子番号はBpn[]関数を使ってBi =４=Bpn[5], Bj=０=Bpn[9]と求まる。ここでBi＞Bjであるので、Bjに対する接触候補の粒子番号はリスト番号BLnumがBs_jgi[Bj - 1]〜Bs_jgi[Bj] - 1 = 0~0の範囲のBlist_j[bl]に格納されている。そこで、このBLnumの範囲でBlist_j[bl] がbiと一致するか調べると、BLnum = 0のときにBlist_j[0] = 4でありBi = 4と一致する。
【００４１】
したがって、現在のリスト番号Lnum = 6の候補ペアは１ステップ前もペアであった事がわかり、そのリスト番号はBLnum = 0である。そしてBLnum =0の法線ベクトルと接線方向の接触力の逆符号の値をLnum = 6の配列に代入する（Spring [6]=−BSpring [0]）ことで、１ステップ前の接触状態を引き継ぐ事ができる。
【００４２】
BSpring [0]は4番の粒子から0番の粒子に対する法線ベクトルおよび接触力である。Spring[6]は9番の粒子から5番の粒子に対する法線ベクトルおよび接触力である。したがって、 Spring [6]=−BSpring [0]である。逆にBj > Biの条件であればSpring []=BSpring []となる
【００４３】
接触力参照表作成部１７は、box_idから対応する候補リスト番号Lnumを算出する関数Ref[i][box_id]を作る。本実施形態では、この関数Ref[i][box_id]を接触力参照表５４とする。接触力参照表５４は、図９に示すように、接触候補リスト５３のi粒子がｊ粒子から受ける接触力Fc[Lnum]を任意の一粒子ｉ（図９ではi＝５）ごとに参照先（Lnum）をまとめたものである。
【００４４】
例として5番と9番の接触候補ペア粒子について説明する(図１０)。 i=5番の粒子はk=box_id=1の箱（nei_pn）を調べて接触候補粒子の番号がj=9でありその候補リスト番号がLnum=6であることを知っている。したがって、 Ref[5][1]=6すなわちRef[i][k]=Lnum である。
【００４５】
一方、ｊ＝９番の粒子に対するRef関数は、9番の箱（nei_pn）の後ろの領域（N−n_jli−1<box_id<N）にi=5番の粒子が格納されているはずなので、その箱を調べてi=5番の粒子が入ってる箱番号kを求める。ここではk = N‐1であり、 Ref[9][3]=6すなわちRef[j][n_jgi[j]−k+N−1] =Lnumである。
【００４６】
接触判定部１８は、接触候補リスト５３を参照して、２つの粒子間の接触判定をすると共に、接触していると判定した場合には粒子間の接触力Fijr、Fijtを算出する。接触候補リスト５３を走査してi 粒子とj 粒子の粒子間距離を計算し、接触している場合は、上述したＶｏｉｇｔモデルを用いて接触力の並進成分Fn +Ft と回転成分Ri ×Ft を計算し、それぞれを接触候補ペアリスト番号Lnumを要素としてもつ配列Fijt[Lnum]，Fijr[Lnum]に記憶する。
【００４７】
接触判定および接触力は接触候補ペアリスト番号をスレッド化して計算する。GPU のグローバルメモリはアクセスパターンがスレッドの順である場合に最も良いアクセス効率が得られる。しかし、ペア粒子の番号はリスト番号に対してランダムであるため、粒子情報が格納されているグローバルメモリをランダムにアクセスしメモリアクセス速度が大きく低下する。一方、テクスチャメモリはキャッシュ機能が有るため、ランダムアクセスによるメモリアクセス速度の低下を抑えることができる。そこで、粒子情報を格納しているグローバルメモリ領域をテクスチャメモリから参照できるように設定した。さらに粒子には所属するセル番号順になるように番号を付けているために接触候補ペア同士の粒子番号が近く、加えて接触候補ペアリストは粒子番号i が小さい順にリスト内に並ぶように作られている。そのため、接触候補ペアリスト番号でスレッド化することでキャッシュヒット率が高くなり、粒子間の接触判定および接触力計算のルーチンが高速化される。
【００４８】
接触力計算部１９は、接触力参照表５４を用いて各粒子に働く接触力Fr、Ftを計算する。
【００４９】
i番の粒子に働く全接触力(並進成分)は接触力参照表５４を用いて、0〜n_jgi−1番の箱（box_id）はそのまま接触力を参照し、n_jgi〜n_jgi+n_jli−1番の箱（box_id）は接触力を逆符号にして参照し、それらの総和で次式のように表わされる。なぜなら、0〜num_g−1番の箱はp_id番の粒子が受ける接触力を参照するのに対して、 n_jgi〜n_jgi+n_jli−1番の箱はi番の粒子が相手粒子に与える接触力を参照するため、i番の粒子が受ける力はその逆符号の値になる。
【数５】
【００５０】
i番の粒子に働く全接触力（回転成分）は接触力参照表５４を用いて、0〜n_jgi−1番の箱（box_id）はそのまま接触力を参照し、n_jgi〜n_jgi+n_jli−1番の箱（box_id）は接触力にβ を掛けてから足し合わせ、次式のように表される。ここで、β はj粒子 の中心から接触点までの距離をi粒子 の中心から接触点までの距離で割った値である。
【数６】
（１０）
【００５１】
コンピュータを用いた演算においては、具体的には粒子番号iをスレッド化し、box_idでループを回してFtall, Frallを求める。ここでFijへのメモリアクセスはランダムアクセスになるのでテクスチャメモリを使用する。上述のようにテクスチャメモリはキャッシュヒット率は高いので、グローバルメモリを直接アクセスするよりは高速に処理できる。なお接触力参照表５４に関しては配列Ref[A][B]のBの要素サイズを16の倍数になるように配列宣言することで連続アクセスが可能となる。
【００５２】
粒子情報更新部２０は、接触力計算部１９により算出された各粒子の接触力と、各粒子の位置、速度に応じて粒子情報を更新する。具体的には、並進の運動方程式をleap-flog法を用いて解くことで新しい時刻の並進速度ならびに粒子座標を求める。さらに、回転の運動方程式も同様にleap-flog 法を用いて解き、ここで得られる回転速度を予測値とする。この予測される回転速度に対して転がり摩擦力を計算し、接触力による回転力と転がり摩擦力を考慮した回転の運動方程式を解くことにより新しい時刻の回転速度を得る。
【００５３】
次に、図１１を参照して、粒子シミュレーション装置１０において実行される処理について説明すると共に、本実施形態に係る粒子シミュレーション方法について説明する。図１１は本実施形態の粒子シミュレーション装置１０において実行される処理を示すフローチャートである。粒子情報取得部１２により粒子情報が粒子情報保持部１１から取得される（Ｓ１０１）。接触候補リスト更新判定部１３により、接触候補リストの更新が必要か否かが判定される（Ｓ１０２）。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。更新が必要と判定された場合にはステップＳ１０３に移行し、更新が不要と判定された場合にはステップＳ１０７に移行する。
【００５４】
次に、粒子番号変更部１４により粒子番号が変更される（Ｓ１０３：セル番号取得ステップ、粒子番号変更ステップ）。粒子番号変更部１４は、粒子を所属するセル番号順に並べ、その順に粒子番号を付け直す。続いて、近傍粒子表作成部１５により、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２が作成される（Ｓ１０４：近傍粒子選択ステップ）。
【００５５】
次に、接触候補リスト作成部１６により、作業空間５１の全粒子について、接触している可能性がある粒子ペアの情報を含む接触候補リスト５３が作成される（Ｓ１０５：接触候補粒子選択ステップ）。接触力参照表作成部１７により、接触力参照表５４が作成される（Ｓ１０６）。
【００５６】
次に、接触判定部１８により、接触候補リスト５３を参照して、２つの粒子間の接触判定が行われ（Ｓ１０７：接触判定ステップ）、接触していると判定された場合には粒子間の接触力Fijr、Fijtが算出され、配列に格納される（Ｓ１０８：接触力演算ステップ、接触力計算ステップ、接触力格納ステップ）。
【００５７】
そして、接触力計算部１９により、接触力参照表５４を参照して、粒子と接触する他の粒子が選択され、（９），（１０）式を用いて、対応する配列Fijr、Fijtから接触力が取得され総和演算されて、各粒子に働く接触力Fr、Ftが計算される（Ｓ１０９：接触力演算ステップ、総和演算ステップ）。粒子情報更新部２０により、ステップＳ１０９において算出された接触力、現在の粒子情報に基づいて粒子情報が更新される（Ｓ１１０：粒子情報更新ステップ）。
【００５８】
そして、粒子シミュレーションが終了したか否かが判定される（Ｓ１１１）。終了判定は、例えば所定時間が経過した場合にシミュレーションが終了したものと判定される。粒子シミュレーションが終了していないと判定された場合にはステップＳ１０１に戻り、次のステップに入る。粒子シミュレーションが終了したと判定された場合には処理を終了する。
【００５９】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、作業空間５１の複数の粒子をセルに割り当てるときに、粒子番号変更部１４が、各粒子が配置されるセルのセル番号を取得し、このセル番号に基づいて粒子の粒子番号を付け替えるため、粒子をセルに割り当てる際にメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００６０】
また、接触判定部１８が算出した粒子間の接触力情報を配列に格納し、接触力計算部１９が、この配列を参照して粒子ごとの接触力を総和演算で得るため、接触力の総和演算時に発生するメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００６１】
（第２実施形態）
次に、本発明の第２実施形態について説明する。図１２は、本発明の第２実施形態に係る粒子シミュレーション装置３０の機能ブロック図である。図１２に示すように、本実施形態に係る粒子シミュレーション装置３０と第１実施形態に係る粒子シミュレーション装置１０との相違点は、（１）接触力参照表作成部１７を備えず、接触力参照表を作成しない点、（２）接触力計算部３３が、接触力参照表ではなく接触候補リストを参照して各粒子の接触力を総和演算する点、（３）近傍粒子表作成部３１が、粒子番号に関係なく近傍粒子表を作成する点である。図１２に示す本実施形態の粒子シミュレーション装置３０に含まれるほかの構成要素は、第１実施形態と同様の機能を有するものなので説明を省略する。
【００６２】
近傍粒子表作成部３１は、各粒子ごとに、図１３に示すような近傍粒子を格納する箱nei_pn[i][box_id]をn個用意し、近傍粒子のうち粒子中心間距離がdis_lim以下の近傍粒子を番号が小さい方の箱から入れる。図１３に示す一連の箱の集合を、作業空間５１内の全ての粒子ｉについて用意したものを、本実施形態では近傍粒子表６２という。本実施形態では、近傍粒子の粒子番号ｊの大小に関係なく図１３では左詰で近傍粒子表６２を作成する。さらに、箱に入れた粒子数n_jiをカウントする。
【００６３】
図１３の例では、例えば５番の粒子（ｉ＝５）の近傍粒子の粒子番号が１，３，４，７，９番であったとすると、
nei_pn[5][0]=7
nei_pn[5][1]=9
nei_pn[5][2]=4
nei_pn[5][3]=3
nei_pn[5][4]=1
n_ji[5]=5
となる。
【００６４】
接触候補リスト作成部３２は、まず、n_jiのプレフィックス和s_jiを次式から求める。
【数７】
【００６５】
そして、各粒子iに対して、box_id を 0 〜 n_ji [i ]−1まで変化させて、次式を用いて、候補リスト番号Lnumおよび接触候補相手の粒子番号jを求め、接触候補リスト６３list_i，list_jを作る。
Lnum = s_ji [i - 1] + box_id
j = nei_pn[i][box_id]
list_i[Lnum] = i
list_j[Lnum] = j
【００６６】
接触力計算部３３は、接触候補リスト６３を用いて各粒子に働く接触力Fr、Ftを計算する。
【００６７】
i番の粒子に働く全接触力(並進成分)は接触候補リストlist_jを用いて、0〜n_ji[i]−1番のリストを検索してj粒子を呼び出し、j粒子との接触力Fij[box_id]を計算しそのまま接触力を足し合わせる（（１２）式）。iに対してスレッド化し、作用反作用を使わないのでメモリ同時書き込みは起きない。
【数８】
【００６８】
i番の粒子に働く全接触力（回転成分）についても並進成分と同様に、（１３）式を用いて算出する。
【数９】
【００６９】
図１４は本実施形態の粒子シミュレーション装置３０において実行される処理を示すフローチャートである。図１４に示すように、本実施形態に係るフローチャートと、図１１に示した第１実施形態に係るフローチャートとの相違点は、（１）接触力参照表作成ステップ（Ｓ１０６）を含まない点と、（２）接触力計算部３３により、接触候補リスト６３を参照して、粒子と接触する他の粒子が選択され、（１２），（１３）式を用いて、対応する配列Fijr、Fijtから接触力が取得され総和演算されて、各粒子に働く接触力Fr、Ftが計算される（Ｓ２０９）点である。
【００７０】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、作業空間５１の複数の粒子をセルに割り当てるときに、粒子番号変更部１４が、各粒子が配置されるセルのセル番号を取得し、このセル番号に基づいて粒子の粒子番号を付け替えるため、粒子をセルに割り当てる際にメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００７１】
（第３実施形態）
次に、本発明の第３実施形態について説明する。図１５は、本発明の第３実施形態に係る粒子シミュレーション装置４０の機能ブロック図である。図１５に示すように、本実施形態に係る粒子シミュレーション装置４０と第１実施形態に係る粒子シミュレーション装置１０との相違点は、粒子番号変更部１４の代わりに粒子格納部４１を備える点である。つまり、本実施形態では、粒子番号の変更は行われない。図１５に示す本実施形態の粒子シミュレーション装置４０に含まれるほかの構成要素は、第１実施形態と同様の機能を有するものなので説明を省略する。
【００７２】
粒子格納部４１は、作業空間５１内の全ての粒子をいずれかのセルに格納する。各粒子は、自身が所属するセルのセル番号と、このセルに格納された順位とを記憶する。
【００７３】
例えば、図１６に示すように、一つのセルに最大で粒子が４個格納される場合を考える。まず、全粒子（ｉ＝１，５，６，８）をセルに格納しようとするが、一度には一粒子しかセルに格納することができない。この例の場合、まず５番粒子だけがセルに格納される。つぎに、再び全粒子をセルに入れようとする。この時、すでにセルに入っている粒子を調べて、その粒子以外をセルに入れようとする。図１６の例の場合、５番以外の粒子（ｉ＝１，６，８）を格納しようとして８番の粒子だけが格納される。この操作を4回繰り返し、全粒子（ｉ＝１，５，６，８）をセルに格納する。
【００７４】
図１７は本実施形態の粒子シミュレーション装置４０において実行される処理を示すフローチャートである。図１７に示すように、本実施形態に係るフローチャートと、図１１に示した第１実施形態に係るフローチャートとの相違点は、粒子番号変更ステップ（Ｓ１０３）の代わりに、粒子格納部４１により、作業空間５１内の全ての粒子がいずれかのセルに格納される（Ｓ３０３）点である。
【００７５】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、接触判定部１８が算出した粒子間の接触力情報を配列に格納し、接触力計算部１９が、この配列を参照して粒子ごとの接触力を総和演算で得るため、接触力の総和演算時に発生するメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【符号の説明】
【００７６】
１０，３０，４０…粒子シミュレーション装置、１１…粒子情報保持部、１２…粒子情報取得部、１３…接触候補リスト更新判定部、１４…粒子番号変更部（セル番号取得手段、粒子番号変更手段）、１５，３１…近傍粒子表作成部（近傍粒子選択手段）、１６，３２…接触候補リスト作成部（接触候補粒子選択手段）、１７…接触力参照表作成部、１８…接触判定部（接触判定手段、接触力演算手段、接触力計算手段、接触力格納手段）、１９，３３…接触力計算部（接触力演算手段、総和演算手段）、２０…粒子情報更新部（粒子情報更新手段）、４１…粒子格納部、５１…作業領域、５２，６２…近傍粒子表、５３，６３…接触候補リスト、５４…接触力参照表。
【技術分野】
【０００１】
本発明は、粒子シミュレーション装置及び粒子シミュレーション方法に関する。
【背景技術】
【０００２】
土砂崩れや地盤液状化などの自然災害予測や粉粒体現象の解析を高精度に行うためには、数値シミュレーションを用いて実験による観察が難しい微視的な知見を得ることが重要である。しかし多くの場合、固体粒子群は数えきれないほどの多数の粒子が多点で同時接触をしながら運動している。したがって、固体粒子群の運動を正確に予測するためには、物理に基づき、可能な限り多数の粒子を扱い、ときには長時間スケールを扱えるシミュレーションが必要となる。
【０００３】
固体粒子群の運動を扱うシミュレーション手法として、Ｖｏｉｇｔモデルを用いて粘弾性体としての粒子運動を再現するＤＥＭ（ＤｉｓｃｒｅｔｅＥｌｅｍｅｎｔ Ｍｅｔｈｏｄ：離散要素法）が知られている（例えば特許文献１）。ＤＥＭは、多体衝突および摩擦や粒子の回転運動が扱えるため、高濃度で多数の粒子を扱う粉体産業や自然科学の分野で広く用いられてきた。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００７−１７２１６９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかし、ＤＥＭ で扱える粒子数は現実よりも遥かに少なく、実現象を正確に再現するためにＤＥＭ に対する高速化アルゴリズムの開発が求められている。
【０００６】
これまでＤＥＭは、接触相手粒子の探索方法の効率化や並列化およびベクトル化を行うことで高速化されてきた。しかし、ＤＥＭにとって並列化やベクトル化は容易ではない。並列化ではノード毎に粒子数が異なる上にノード間で情報交換すべき粒子が不定であるために並列化率が向上し難い。また、ベクトル化ではメモリへの書き込み競合を回避する処理が多く必要となるためにその処理で計算時間を浪費する。そのため、ＣＰＵ（Central Processing Unit）がムーアの法則を破り桁違いに高性能化されない限り、ＤＥＭの飛躍的な高速化は期待できない。
【０００７】
一方、最近ではＣＰＵと比べて価格に対する演算性能の比が格段に大きいＧＰＵ（Graphics Processing Unit）を用いた高速化が盛んに行われている。しかし、粒子に働く力の全てを考慮している完全なＤＥＭに対してＧＰＵ を用いて高速化した例は無く、さらにＧＰＵは並列計算を行うため、メモリ書き込み競合を避けるべく計算量を増大させているためにＧＰＵの性能を十分に生かせていない。
【０００８】
本発明は、上記課題を解決するためになされたものであり、メモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる粒子シミュレーション装置及び粒子シミュレーション方法を提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記課題を解決するため、本発明に係る粒子シミュレーション装置は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、複数の粒子のそれぞれには粒子番号が付されており、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得手段と、セル番号取得手段により取得されたセル番号に基づいて、複数の粒子の前記粒子番号を付け替える粒子番号変更手段と、粒子番号変更手段により付け替えられた粒子番号に基づき、複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択手段と、近傍粒子選択手段により選択された近傍粒子と前記１つの粒子との位置関係に基づいて、１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算手段と、接触力計算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段とを備えることを特徴とする。
【００１０】
同様に、上記課題を解決するため、本発明に係る粒子シミュレーション方法は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、複数の粒子のそれぞれには粒子番号が付されており、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得ステップと、セル番号取得ステップにおいて取得されたセル番号に基づいて、複数の粒子の前記粒子番号を付け替える粒子番号変更ステップと、粒子番号変更ステップにおいて付け替えられた粒子番号に基づき、複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択ステップと、近傍粒子選択ステップにおいて選択された近傍粒子と前記１つの粒子との位置関係に基づいて、１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算ステップと、接触力計算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップとを備えることを特徴とする。
【００１１】
このような粒子シミュレーション装置及び粒子シミュレーション方法によれば、作業空間の複数の粒子をセルに割り当てるときに、各粒子が配置されるセルのセル番号が取得され、このセル番号に基づいて粒子の粒子番号が付け替えられるため、粒子をセルに割り当てる際にメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００１２】
上記課題を解決するため、本発明に係る粒子シミュレーション装置は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算手段と、接触力計算手段により計算された粒子間の接触力情報を配列に格納する接触力格納手段と、接触力格納手段により粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算手段と、総和演算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段と、を備えることを特徴とする。
【００１３】
同様に、上記課題を解決するため、本発明に係る粒子シミュレーション方法は、作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算ステップと、接触力計算ステップにおいて計算された粒子間の接触力情報を配列に格納する接触力格納ステップと、接触力格納ステップにおいて粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算ステップと、総和演算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップと、を備えることを特徴とする。
【００１４】
このような粒子シミュレーション装置及び粒子シミュレーション方法によれば、算出された粒子間の接触力情報が配列に格納され、この配列を参照して粒子ごとの接触力が総和演算で得られるため、接触力の総和演算時に発生するメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【発明の効果】
【００１５】
本発明に係る粒子シミュレーション装置及び粒子シミュレーション方法によれば、メモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【図面の簡単な説明】
【００１６】
【図１】本発明の第１実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図２】本実施形態で粒子モデルに適用するＶｏｉｇｔモデルを示す概略図である。
【図３】本実施形態における作業領域を示す図である。
【図４】粒子番号変更の流れを示す図である。
【図５】作業領域における粒子番号を変更する処理と、各セルに粒子番号の最大値及び最小値を記憶する処理を示す図である。
【図６】近傍粒子表の構成の一例を示す図である。
【図７】接触候補リストの構成の一例を示す図である。
【図８】現在の接触候補リストのペア粒子について１ステップ前にペアであったか調べる処理を示す図である。
【図９】接触力参照表の構成の一例を示す図である。
【図１０】接触力参照表の構成の一例を示す図である。
【図１１】本発明の第１実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【図１２】本発明の第２実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図１３】近傍粒子表の構成の一例を示す図である。
【図１４】本発明の第２実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【図１５】本発明の第３実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図１６】粒子をセルに格納する流れを示す図である。
【図１７】本発明の第３実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【発明を実施するための形態】
【００１７】
以下、図面を参照しながら本発明の実施形態を詳細に説明する。なお、図面の説明において同一又は同等の要素には同一の符号を付し、重複する説明を省略する。
【００１８】
（第1実施形態）
図１は、本発明の第１実施形態に係る粒子シミュレーション装置１０の機能ブロック図である。図１に示すように、粒子シミュレーション装置１０は、粒子情報保持部１１、粒子情報取得部１２、接触候補リスト更新判定部１３、粒子番号変更部（セル番号取得手段、粒子番号変更手段）１４、近傍粒子表作成部（近傍粒子選択手段）１５、接触候補リスト作成部（接触候補粒子選択手段）１６、接触力参照表作成部１７、接触判定部（接触判定手段、接触力演算手段、接触力計算手段、接触力格納手段）１８、接触力計算部（接触力演算手段、総和演算手段）１９、粒子情報更新部（粒子情報更新手段）２０を備えている。
【００１９】
粒子シミュレーション装置１０は、物理的には、ＧＰＵ（Graphics Processing Unit）、主記憶装置であるＲＡＭ及びＲＯＭ、ハードディスク装置等の補助記憶装置、入力デバイスである入力キー等の入力装置、ディスプレイ等の出力装置、ネットワークカード等のデータ送受信デバイスである通信モジュールなどを有するコンピュータとして構成されている。図１において説明した粒子シミュレーション装置１０の各機能は、ＧＰＵ、ＲＡＭ等のハードウェア上に所定のコンピュータソフトウェアを読み込ませることにより、ＧＰＵの制御のもとで入力装置、出力装置、通信モジュールを動作させるとともに、ＲＡＭや補助記憶装置におけるデータの読み出し及び書き込みを行うことで実現される。
【００２０】
まず本実施形態で適用する粒子モデルについて説明する。粒子間接触力の法線方向および接線方向成分Ｆｎ，Ｆｔは、図２に示すＶｏｉｇｔモデル７１を用いて以下の式により表される。
【数１】
ここでｘとｖは相対変位および相対速度ベクトル、μtは滑り摩擦係数、ｍｉｎは絶対値が小さい方の値を表す関数、Ｋは線形バネ７２における弾性係数、ηはダッシュポット７３の粘性減衰係数であり、反発係数eと次式の関連がある。
【数２】
【００２１】
したがって、ＤＥＭで扱う並進および重心を中心とする回転の運動方程式は以下の式で表される。
【数３】
ここで、ｕｐとωは粒子の速度および角速度、ｍｐとＩｐは粒子の質量および慣性モーメント、Ｒは粒子中心から接触点へ向かう位置ベクトルである。また、（５）式の右辺第二項は転がり摩擦力を表し、μｒは転がり摩擦係数、ｂは接触面幅である。（１）、（２）式のＶｏｉｇｔモデル７１は、後述する接触判定部１８で粒子間接触力を算出するのに用いられ、また（４）、（５）式の運動方程式は、粒子情報更新部２０で粒子の座標を算出するのに用いられる。
【００２２】
また、本実施形態において粒子が運動する領域である作業領域５１は、図３に示すように、一辺の大きさがDmax（1．0＋α）の立方体セルで区分される。ここで、Dmaxは最大粒子径、αは更新頻度を調整するパラメータであり例えばα＝０．２である。また、作業領域５１内の各セルにはセル番号cell_idが付されている。セル番号の付け方は後述する。
【００２３】
粒子情報保持部１１は、作業空間５１内の複数の粒子のそれぞれについて粒子情報を保持している。粒子情報は、座標、並進速度、回転速度、粒子半径を含む。座標とは、図３に示す三次元の作業空間５１における三次元座標（ｐｏｓ．ｘ，ｐｏｓ．ｙ，ｐｏｓ．ｚ）である。粒子情報は、後述する粒子情報更新部２０により更新される。
【００２４】
粒子情報取得部１２は、粒子情報保持部１１から粒子情報を取得する。粒子情報取得部１２は、粒子シミュレーションの開始時に粒子情報保持部１１から予め設定されている初期状態の粒子情報を取得し、シミュレーション動作中には、粒子情報保持部１１が更新される度に新たな粒子情報を粒子情報保持部１１から取得する。
【００２５】
接触候補リスト更新判定部１３は、後述する接触候補リスト５３を更新するか否かを判定する。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。また、粒子の積算移動距離は、例えば、粒子情報取得部１２が新しい粒子情報を粒子情報保持部１１から取得する度に、１ステップ前の粒子情報との座標データの差分からこの直近の１ステップにおける粒子の移動距離を算出し、これを粒子ごとに積算して求める。なお、粒子移動距離の積算値は、接触候補リスト５３の更新が行われるたびに初期化される。
【００２６】
粒子番号変更部１４は、粒子を所属するセル番号順に並べ、その順に粒子番号を付け直す。セル番号は図３に示すように一次元であり、例えば次式のように表される。
cell_id＝i.x＋id.x×i.y＋id.x×id.y×i.z （６）
i.m = int(pos.m／cellsize.m)
ここでcell_idはセル番号であり、i．mはｍ方向のセル位置であり、id．mはｍ方向のセル数である。また、pos．mはｍ軸方向の座標を表し、cellsize．mはm軸方向のセルサイズを表し、int()は実数を整数に変換する関数である。
【００２７】
図４に示すように、粒子に自身が所属するセル番号を記憶させる（図４（ａ））。次に、セル番号が昇順になるように粒子を並べる（図４（ｂ））。そして、セル番号順に粒子番号を付け直す（図４（ｃ））。
【００２８】
図４の流れを式で表すと以下のとおりである。
Pcell[p_id]=cell_id
B_pid[p_id]=p_id
Step1
Pcell={5,13,0,21,5,3,6} B_pid={0,1,2,3,4,5,6}
Step2
Pcell={0,3,5,5,6,13,21} B_pid={2,5,0,4,6,1,3}
Step3
old_p_id=B_pid[p_id] Val[p_id]=Val[old_p_id]
ここで、p_idは粒子番号であり、Valは粒子座標、並進速度、回転速度、半径である。Step1のPcellの中身を周知のソート手法（例えばバイトニックソート）を用いて昇順にし、それと同じようにB_pidの中身も移動させる。この時点でBの中身はソート前の粒子番号が記憶されている。Step3でB_pidを使い、ソート前の粒子情報をソート後の粒子情報に書き換える。
【００２９】
この粒子番号を変更する処理により、例えば図５（ａ）に示すように作業空間５１に配置された粒子の粒子番号が、図５（ｂ）に示すように、セル番号に沿った粒子番号に変更される。
【００３０】
近傍粒子表作成部１５は、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２を作成する。近傍粒子表作成部１５は、着目粒子の所属するセルおよびその隣接セルに存在する粒子のうち、自身より大きな粒子番号の粒子で粒子間距離がある値以下の粒子を近傍粒子とする。以下に詳細を説明する。
【００３１】
近傍粒子表作成部１５は、セル内の粒子の番号のうち、最小値Pnum_in_cellminと最大値Pnum_in_cellmaxをセルに記憶させる（図５（ｃ））。番号cell_idのセルに隣接するセルの番号nei_cell_idは次式で表される。
nei_cell_id[n]= cell_id + con[n] (n=0〜26) （７）
ここで、conはid．xとid．yによって決まる定数である。セル番号で粒子番号をソートしてあるため、各隣接セルにはPnum_in_cellmin[nei_cell_id]からPnum_in_cellmax[nei_cell_id]までの番号の粒子が連続して存在することがわかる。これらの粒子を近傍粒子と呼ぶ。
【００３２】
各粒子ごとに、図６に示すような近傍粒子を格納する箱nei_pn[i][box_id]をn個用意し、近傍粒子のうち粒子中心間距離がdis_lim以下で、自身より大きい番号（i=p_id）の近傍粒子を番号（box_id）が小さい方の箱から入れて行き、一方、小さい番号（i）の近傍粒子は番号（box_id）が大きい方の箱から順に入れる。dis_lim = (ri + rj)(1．0+α)であり、ri,rjは近傍粒子ペアi,jの粒子半径、αは接触候補の更新回数を調整するパラメータである。図６に示す一連の箱の集合を、作業空間５１内の全ての粒子ｉについて用意したものを、本実施形態では近傍粒子表５２という。
【００３３】
さらに、粒子番号ｉが自身より大きい粒子および小さい粒子の数n_jgi， n_jliをそれぞれ全粒子に対して求める。図６の例では、例えば５番の粒子（ｉ＝５）の近傍粒子の粒子番号が１，３，４，７，９番であったとすると、
nei_pn[5][0]=7
nei_pn[5][1]=9
nei_pn[5][n-3]=4
nei_pn[5][n-2]=3
nei_pn[5][n-1]=1
n_jgi[5]=2
n_jli[5]=3
となる。
【００３４】
接触候補リスト作成部１６は、まず、図７に示すように、n_jgiのプレフィックス和s_jgiを次式から求める。
【数４】
【００３５】
そして、各粒子iに対して、box_id を 0 〜 n_jgi [i ] − 1まで変化させて、次式を用いて、候補リスト番号Lnumおよび接触候補相手の粒子番号jを求め、接触候補リスト５３list_i，list_jを作る。
Lnum = s_jgi [i - 1] + box_id
j = nei_pn[i][box_id]
list_i[Lnum] = i
list_j[Lnum] = j
【００３６】
例えば、図７の例では、接触候補リスト５３の内容は、以下の表１のようになる。
【表１】
つまり、接触候補リスト５３は、作業空間５１の全粒子について、接触している可能性がある粒子ペアｌｉｓｔ＿ｉ，ｌｉｓｔ＿ｊに番号Ｌｎｕｍを付加する。
【００３７】
次に、１ステップ前の接触候補ペアBlist_i，Blist_jを探索し、現在の接触候補ペアlist_i, list_jが１ステップ前にもペアであったか調べる。ペアであったら、図２に示した粒子モデルにおける法線ベクトルと接線方向のバネによる力（バネの伸縮量）を新しい候補リスト番号の配列に格納する。
【００３８】
粒子番号i=list_i[Lnum], j=list_j[Lnum]の接触候補ペアが有るとき、このペア粒子の１ステップ前の粒子番号はBi =Bpn[i]， Bj = Bpn[j]である。Bj＞BiのときはBs_jgi[Bi−1]〜Bs_jgi[Bi]−1の範囲の接触候補Blist_j[BLnum]を調べBjと一致するか調べる。一致したときはそのリスト番号BLnumを用いて１ステップ前の法線ベクトルと接線方向の接触力BSpring[BLnum]を現在のリスト番号配列Spring[Lnum]に代入する。逆にBi＞Bjのとき、Bjに対して粒子番号が大きい粒子が接触候補として格納されているリスト番号BLnumの範囲はBs_jig[Bj−1]〜Bs_jgi[Bj]−1であるから、この範囲でBlist_j[BLnum]がBiと一致するか調べ、一致するときはSpring[Lnum]に−BSpring[BLnum]を代入する。
【００３９】
この処理について、図８を参照して、現在の接触候補リスト番号Lnum = 6におけるｉ＝５，ｊ＝９番のペア粒子について１ステップ前にペアであったか調べる処理を説明する。
【００４０】
まず、ｉ＝５，ｊ＝９の１ステップ前の粒子番号はBpn[]関数を使ってBi =４=Bpn[5], Bj=０=Bpn[9]と求まる。ここでBi＞Bjであるので、Bjに対する接触候補の粒子番号はリスト番号BLnumがBs_jgi[Bj - 1]〜Bs_jgi[Bj] - 1 = 0~0の範囲のBlist_j[bl]に格納されている。そこで、このBLnumの範囲でBlist_j[bl] がbiと一致するか調べると、BLnum = 0のときにBlist_j[0] = 4でありBi = 4と一致する。
【００４１】
したがって、現在のリスト番号Lnum = 6の候補ペアは１ステップ前もペアであった事がわかり、そのリスト番号はBLnum = 0である。そしてBLnum =0の法線ベクトルと接線方向の接触力の逆符号の値をLnum = 6の配列に代入する（Spring [6]=−BSpring [0]）ことで、１ステップ前の接触状態を引き継ぐ事ができる。
【００４２】
BSpring [0]は4番の粒子から0番の粒子に対する法線ベクトルおよび接触力である。Spring[6]は9番の粒子から5番の粒子に対する法線ベクトルおよび接触力である。したがって、 Spring [6]=−BSpring [0]である。逆にBj > Biの条件であればSpring []=BSpring []となる
【００４３】
接触力参照表作成部１７は、box_idから対応する候補リスト番号Lnumを算出する関数Ref[i][box_id]を作る。本実施形態では、この関数Ref[i][box_id]を接触力参照表５４とする。接触力参照表５４は、図９に示すように、接触候補リスト５３のi粒子がｊ粒子から受ける接触力Fc[Lnum]を任意の一粒子ｉ（図９ではi＝５）ごとに参照先（Lnum）をまとめたものである。
【００４４】
例として5番と9番の接触候補ペア粒子について説明する(図１０)。 i=5番の粒子はk=box_id=1の箱（nei_pn）を調べて接触候補粒子の番号がj=9でありその候補リスト番号がLnum=6であることを知っている。したがって、 Ref[5][1]=6すなわちRef[i][k]=Lnum である。
【００４５】
一方、ｊ＝９番の粒子に対するRef関数は、9番の箱（nei_pn）の後ろの領域（N−n_jli−1<box_id<N）にi=5番の粒子が格納されているはずなので、その箱を調べてi=5番の粒子が入ってる箱番号kを求める。ここではk = N‐1であり、 Ref[9][3]=6すなわちRef[j][n_jgi[j]−k+N−1] =Lnumである。
【００４６】
接触判定部１８は、接触候補リスト５３を参照して、２つの粒子間の接触判定をすると共に、接触していると判定した場合には粒子間の接触力Fijr、Fijtを算出する。接触候補リスト５３を走査してi 粒子とj 粒子の粒子間距離を計算し、接触している場合は、上述したＶｏｉｇｔモデルを用いて接触力の並進成分Fn +Ft と回転成分Ri ×Ft を計算し、それぞれを接触候補ペアリスト番号Lnumを要素としてもつ配列Fijt[Lnum]，Fijr[Lnum]に記憶する。
【００４７】
接触判定および接触力は接触候補ペアリスト番号をスレッド化して計算する。GPU のグローバルメモリはアクセスパターンがスレッドの順である場合に最も良いアクセス効率が得られる。しかし、ペア粒子の番号はリスト番号に対してランダムであるため、粒子情報が格納されているグローバルメモリをランダムにアクセスしメモリアクセス速度が大きく低下する。一方、テクスチャメモリはキャッシュ機能が有るため、ランダムアクセスによるメモリアクセス速度の低下を抑えることができる。そこで、粒子情報を格納しているグローバルメモリ領域をテクスチャメモリから参照できるように設定した。さらに粒子には所属するセル番号順になるように番号を付けているために接触候補ペア同士の粒子番号が近く、加えて接触候補ペアリストは粒子番号i が小さい順にリスト内に並ぶように作られている。そのため、接触候補ペアリスト番号でスレッド化することでキャッシュヒット率が高くなり、粒子間の接触判定および接触力計算のルーチンが高速化される。
【００４８】
接触力計算部１９は、接触力参照表５４を用いて各粒子に働く接触力Fr、Ftを計算する。
【００４９】
i番の粒子に働く全接触力(並進成分)は接触力参照表５４を用いて、0〜n_jgi−1番の箱（box_id）はそのまま接触力を参照し、n_jgi〜n_jgi+n_jli−1番の箱（box_id）は接触力を逆符号にして参照し、それらの総和で次式のように表わされる。なぜなら、0〜num_g−1番の箱はp_id番の粒子が受ける接触力を参照するのに対して、 n_jgi〜n_jgi+n_jli−1番の箱はi番の粒子が相手粒子に与える接触力を参照するため、i番の粒子が受ける力はその逆符号の値になる。
【数５】
【００５０】
i番の粒子に働く全接触力（回転成分）は接触力参照表５４を用いて、0〜n_jgi−1番の箱（box_id）はそのまま接触力を参照し、n_jgi〜n_jgi+n_jli−1番の箱（box_id）は接触力にβ を掛けてから足し合わせ、次式のように表される。ここで、β はj粒子 の中心から接触点までの距離をi粒子 の中心から接触点までの距離で割った値である。
【数６】
（１０）
【００５１】
コンピュータを用いた演算においては、具体的には粒子番号iをスレッド化し、box_idでループを回してFtall, Frallを求める。ここでFijへのメモリアクセスはランダムアクセスになるのでテクスチャメモリを使用する。上述のようにテクスチャメモリはキャッシュヒット率は高いので、グローバルメモリを直接アクセスするよりは高速に処理できる。なお接触力参照表５４に関しては配列Ref[A][B]のBの要素サイズを16の倍数になるように配列宣言することで連続アクセスが可能となる。
【００５２】
粒子情報更新部２０は、接触力計算部１９により算出された各粒子の接触力と、各粒子の位置、速度に応じて粒子情報を更新する。具体的には、並進の運動方程式をleap-flog法を用いて解くことで新しい時刻の並進速度ならびに粒子座標を求める。さらに、回転の運動方程式も同様にleap-flog 法を用いて解き、ここで得られる回転速度を予測値とする。この予測される回転速度に対して転がり摩擦力を計算し、接触力による回転力と転がり摩擦力を考慮した回転の運動方程式を解くことにより新しい時刻の回転速度を得る。
【００５３】
次に、図１１を参照して、粒子シミュレーション装置１０において実行される処理について説明すると共に、本実施形態に係る粒子シミュレーション方法について説明する。図１１は本実施形態の粒子シミュレーション装置１０において実行される処理を示すフローチャートである。粒子情報取得部１２により粒子情報が粒子情報保持部１１から取得される（Ｓ１０１）。接触候補リスト更新判定部１３により、接触候補リストの更新が必要か否かが判定される（Ｓ１０２）。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。更新が必要と判定された場合にはステップＳ１０３に移行し、更新が不要と判定された場合にはステップＳ１０７に移行する。
【００５４】
次に、粒子番号変更部１４により粒子番号が変更される（Ｓ１０３：セル番号取得ステップ、粒子番号変更ステップ）。粒子番号変更部１４は、粒子を所属するセル番号順に並べ、その順に粒子番号を付け直す。続いて、近傍粒子表作成部１５により、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２が作成される（Ｓ１０４：近傍粒子選択ステップ）。
【００５５】
次に、接触候補リスト作成部１６により、作業空間５１の全粒子について、接触している可能性がある粒子ペアの情報を含む接触候補リスト５３が作成される（Ｓ１０５：接触候補粒子選択ステップ）。接触力参照表作成部１７により、接触力参照表５４が作成される（Ｓ１０６）。
【００５６】
次に、接触判定部１８により、接触候補リスト５３を参照して、２つの粒子間の接触判定が行われ（Ｓ１０７：接触判定ステップ）、接触していると判定された場合には粒子間の接触力Fijr、Fijtが算出され、配列に格納される（Ｓ１０８：接触力演算ステップ、接触力計算ステップ、接触力格納ステップ）。
【００５７】
そして、接触力計算部１９により、接触力参照表５４を参照して、粒子と接触する他の粒子が選択され、（９），（１０）式を用いて、対応する配列Fijr、Fijtから接触力が取得され総和演算されて、各粒子に働く接触力Fr、Ftが計算される（Ｓ１０９：接触力演算ステップ、総和演算ステップ）。粒子情報更新部２０により、ステップＳ１０９において算出された接触力、現在の粒子情報に基づいて粒子情報が更新される（Ｓ１１０：粒子情報更新ステップ）。
【００５８】
そして、粒子シミュレーションが終了したか否かが判定される（Ｓ１１１）。終了判定は、例えば所定時間が経過した場合にシミュレーションが終了したものと判定される。粒子シミュレーションが終了していないと判定された場合にはステップＳ１０１に戻り、次のステップに入る。粒子シミュレーションが終了したと判定された場合には処理を終了する。
【００５９】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、作業空間５１の複数の粒子をセルに割り当てるときに、粒子番号変更部１４が、各粒子が配置されるセルのセル番号を取得し、このセル番号に基づいて粒子の粒子番号を付け替えるため、粒子をセルに割り当てる際にメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００６０】
また、接触判定部１８が算出した粒子間の接触力情報を配列に格納し、接触力計算部１９が、この配列を参照して粒子ごとの接触力を総和演算で得るため、接触力の総和演算時に発生するメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００６１】
（第２実施形態）
次に、本発明の第２実施形態について説明する。図１２は、本発明の第２実施形態に係る粒子シミュレーション装置３０の機能ブロック図である。図１２に示すように、本実施形態に係る粒子シミュレーション装置３０と第１実施形態に係る粒子シミュレーション装置１０との相違点は、（１）接触力参照表作成部１７を備えず、接触力参照表を作成しない点、（２）接触力計算部３３が、接触力参照表ではなく接触候補リストを参照して各粒子の接触力を総和演算する点、（３）近傍粒子表作成部３１が、粒子番号に関係なく近傍粒子表を作成する点である。図１２に示す本実施形態の粒子シミュレーション装置３０に含まれるほかの構成要素は、第１実施形態と同様の機能を有するものなので説明を省略する。
【００６２】
近傍粒子表作成部３１は、各粒子ごとに、図１３に示すような近傍粒子を格納する箱nei_pn[i][box_id]をn個用意し、近傍粒子のうち粒子中心間距離がdis_lim以下の近傍粒子を番号が小さい方の箱から入れる。図１３に示す一連の箱の集合を、作業空間５１内の全ての粒子ｉについて用意したものを、本実施形態では近傍粒子表６２という。本実施形態では、近傍粒子の粒子番号ｊの大小に関係なく図１３では左詰で近傍粒子表６２を作成する。さらに、箱に入れた粒子数n_jiをカウントする。
【００６３】
図１３の例では、例えば５番の粒子（ｉ＝５）の近傍粒子の粒子番号が１，３，４，７，９番であったとすると、
nei_pn[5][0]=7
nei_pn[5][1]=9
nei_pn[5][2]=4
nei_pn[5][3]=3
nei_pn[5][4]=1
n_ji[5]=5
となる。
【００６４】
接触候補リスト作成部３２は、まず、n_jiのプレフィックス和s_jiを次式から求める。
【数７】
【００６５】
そして、各粒子iに対して、box_id を 0 〜 n_ji [i ]−1まで変化させて、次式を用いて、候補リスト番号Lnumおよび接触候補相手の粒子番号jを求め、接触候補リスト６３list_i，list_jを作る。
Lnum = s_ji [i - 1] + box_id
j = nei_pn[i][box_id]
list_i[Lnum] = i
list_j[Lnum] = j
【００６６】
接触力計算部３３は、接触候補リスト６３を用いて各粒子に働く接触力Fr、Ftを計算する。
【００６７】
i番の粒子に働く全接触力(並進成分)は接触候補リストlist_jを用いて、0〜n_ji[i]−1番のリストを検索してj粒子を呼び出し、j粒子との接触力Fij[box_id]を計算しそのまま接触力を足し合わせる（（１２）式）。iに対してスレッド化し、作用反作用を使わないのでメモリ同時書き込みは起きない。
【数８】
【００６８】
i番の粒子に働く全接触力（回転成分）についても並進成分と同様に、（１３）式を用いて算出する。
【数９】
【００６９】
図１４は本実施形態の粒子シミュレーション装置３０において実行される処理を示すフローチャートである。図１４に示すように、本実施形態に係るフローチャートと、図１１に示した第１実施形態に係るフローチャートとの相違点は、（１）接触力参照表作成ステップ（Ｓ１０６）を含まない点と、（２）接触力計算部３３により、接触候補リスト６３を参照して、粒子と接触する他の粒子が選択され、（１２），（１３）式を用いて、対応する配列Fijr、Fijtから接触力が取得され総和演算されて、各粒子に働く接触力Fr、Ftが計算される（Ｓ２０９）点である。
【００７０】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、作業空間５１の複数の粒子をセルに割り当てるときに、粒子番号変更部１４が、各粒子が配置されるセルのセル番号を取得し、このセル番号に基づいて粒子の粒子番号を付け替えるため、粒子をセルに割り当てる際にメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【００７１】
（第３実施形態）
次に、本発明の第３実施形態について説明する。図１５は、本発明の第３実施形態に係る粒子シミュレーション装置４０の機能ブロック図である。図１５に示すように、本実施形態に係る粒子シミュレーション装置４０と第１実施形態に係る粒子シミュレーション装置１０との相違点は、粒子番号変更部１４の代わりに粒子格納部４１を備える点である。つまり、本実施形態では、粒子番号の変更は行われない。図１５に示す本実施形態の粒子シミュレーション装置４０に含まれるほかの構成要素は、第１実施形態と同様の機能を有するものなので説明を省略する。
【００７２】
粒子格納部４１は、作業空間５１内の全ての粒子をいずれかのセルに格納する。各粒子は、自身が所属するセルのセル番号と、このセルに格納された順位とを記憶する。
【００７３】
例えば、図１６に示すように、一つのセルに最大で粒子が４個格納される場合を考える。まず、全粒子（ｉ＝１，５，６，８）をセルに格納しようとするが、一度には一粒子しかセルに格納することができない。この例の場合、まず５番粒子だけがセルに格納される。つぎに、再び全粒子をセルに入れようとする。この時、すでにセルに入っている粒子を調べて、その粒子以外をセルに入れようとする。図１６の例の場合、５番以外の粒子（ｉ＝１，６，８）を格納しようとして８番の粒子だけが格納される。この操作を4回繰り返し、全粒子（ｉ＝１，５，６，８）をセルに格納する。
【００７４】
図１７は本実施形態の粒子シミュレーション装置４０において実行される処理を示すフローチャートである。図１７に示すように、本実施形態に係るフローチャートと、図１１に示した第１実施形態に係るフローチャートとの相違点は、粒子番号変更ステップ（Ｓ１０３）の代わりに、粒子格納部４１により、作業空間５１内の全ての粒子がいずれかのセルに格納される（Ｓ３０３）点である。
【００７５】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、接触判定部１８が算出した粒子間の接触力情報を配列に格納し、接触力計算部１９が、この配列を参照して粒子ごとの接触力を総和演算で得るため、接触力の総和演算時に発生するメモリへの書込みの競合を防ぐことができ、ＧＰＵなどの並列プロセッサを用いてＤＥＭ（離散要素法）の演算効率を向上させることができる。
【符号の説明】
【００７６】
１０，３０，４０…粒子シミュレーション装置、１１…粒子情報保持部、１２…粒子情報取得部、１３…接触候補リスト更新判定部、１４…粒子番号変更部（セル番号取得手段、粒子番号変更手段）、１５，３１…近傍粒子表作成部（近傍粒子選択手段）、１６，３２…接触候補リスト作成部（接触候補粒子選択手段）、１７…接触力参照表作成部、１８…接触判定部（接触判定手段、接触力演算手段、接触力計算手段、接触力格納手段）、１９，３３…接触力計算部（接触力演算手段、総和演算手段）、２０…粒子情報更新部（粒子情報更新手段）、４１…粒子格納部、５１…作業領域、５２，６２…近傍粒子表、５３，６３…接触候補リスト、５４…接触力参照表。
【特許請求の範囲】
【請求項１】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、
前記複数の粒子のそれぞれには粒子番号が付されており、前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得手段と、
前記セル番号取得手段により取得されたセル番号に基づいて、前記複数の粒子の前記粒子番号を付け替える粒子番号変更手段と、
前記粒子番号変更手段により付け替えられた粒子番号に基づき、前記複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択手段と、
前記近傍粒子選択手段により選択された近傍粒子と前記１つの粒子との位置関係に基づいて、前記１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、
前記選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、
前記接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算手段と、
前記接触力計算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段と
を備えることを特徴とする粒子シミュレーション装置。
【請求項２】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、
前記複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、
前記選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、
前記接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算手段と、
前記接触力計算手段により計算された粒子間の接触力情報を配列に格納する接触力格納手段と、
前記接触力格納手段により粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算手段と、
前記総和演算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段と、
を備えることを特徴とする粒子シミュレーション装置。
【請求項３】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、
前記複数の粒子のそれぞれには粒子番号が付されており、前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得ステップと、
前記セル番号取得ステップにおいて取得されたセル番号に基づいて、前記複数の粒子の前記粒子番号を付け替える粒子番号変更ステップと、
前記粒子番号変更ステップにおいて付け替えられた粒子番号に基づき、前記複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択ステップと、
前記近傍粒子選択ステップにおいて選択された近傍粒子と前記１つの粒子との位置関係に基づいて、前記１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、
前記選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、
前記接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算ステップと、
前記接触力計算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップと
を備えることを特徴とする粒子シミュレーション方法。
【請求項４】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、
前記複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、
前記選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、
前記接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算ステップと、
前記接触力計算ステップにおいて計算された粒子間の接触力情報を配列に格納する接触力格納ステップと、
前記接触力格納ステップにおいて粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算ステップと、
前記総和演算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップと、
を備えることを特徴とする粒子シミュレーション方法。
【請求項１】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、
前記複数の粒子のそれぞれには粒子番号が付されており、前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得手段と、
前記セル番号取得手段により取得されたセル番号に基づいて、前記複数の粒子の前記粒子番号を付け替える粒子番号変更手段と、
前記粒子番号変更手段により付け替えられた粒子番号に基づき、前記複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択手段と、
前記近傍粒子選択手段により選択された近傍粒子と前記１つの粒子との位置関係に基づいて、前記１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、
前記選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、
前記接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算手段と、
前記接触力計算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段と
を備えることを特徴とする粒子シミュレーション装置。
【請求項２】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、
前記複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択手段と、
前記選択候補粒子選択手段により選択された接触候補粒子との間で接触判定を行う接触判定手段と、
前記接触判定手段により前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算手段と、
前記接触力計算手段により計算された粒子間の接触力情報を配列に格納する接触力格納手段と、
前記接触力格納手段により粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算手段と、
前記総和演算手段により計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新手段と、
を備えることを特徴とする粒子シミュレーション装置。
【請求項３】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、
前記複数の粒子のそれぞれには粒子番号が付されており、前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記複数の粒子のそれぞれについて、位置情報に基づいて各粒子が配置されるセルのセル番号を取得するセル番号取得ステップと、
前記セル番号取得ステップにおいて取得されたセル番号に基づいて、前記複数の粒子の前記粒子番号を付け替える粒子番号変更ステップと、
前記粒子番号変更ステップにおいて付け替えられた粒子番号に基づき、前記複数の粒子の１つの粒子の近傍に位置する近傍粒子を選択する近傍粒子選択ステップと、
前記近傍粒子選択ステップにおいて選択された近傍粒子と前記１つの粒子との位置関係に基づいて、前記１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、
前記選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、
前記接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算し、粒子ごとの接触力を総和演算する接触力演算ステップと、
前記接触力計算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップと
を備えることを特徴とする粒子シミュレーション方法。
【請求項４】
作業空間内の複数の粒子について他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション方法であって、
前記複数の粒子のそれぞれの位置関係に基づいて、前記複数の粒子の１つの粒子と接触している可能性の高い接触候補粒子を選択する接触候補粒子選択ステップと、
前記選択候補粒子選択ステップにおいて選択された接触候補粒子との間で接触判定を行う接触判定ステップと、
前記接触判定ステップにおいて前記１つの粒子と接触していると判定された粒子との間の接触力を計算する接触力計算ステップと、
前記接触力計算ステップにおいて計算された粒子間の接触力情報を配列に格納する接触力格納ステップと、
前記接触力格納ステップにおいて粒子間の接触力情報が格納された配列を参照して、粒子ごとの接触力を総和演算する総和演算ステップと、
前記総和演算ステップにおいて計算された粒子ごとの接触力に基づいて粒子の位置及び速度を含む粒子情報を更新する粒子情報更新ステップと、
を備えることを特徴とする粒子シミュレーション方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【公開番号】特開２０１０−２３８０３０（Ｐ２０１０−２３８０３０Ａ）
【公開日】平成２２年１０月２１日（２０１０．１０．２１）
【国際特許分類】
【出願番号】特願２００９−８６２７２（Ｐ２００９−８６２７２）
【出願日】平成２１年３月３１日（２００９．３．３１）
【出願人】（５０４１９４８７８）独立行政法人海洋研究開発機構 (110)
【公開日】平成２２年１０月２１日（２０１０．１０．２１）
【国際特許分類】
【出願日】平成２１年３月３１日（２００９．３．３１）
【出願人】（５０４１９４８７８）独立行政法人海洋研究開発機構 (110)
[ Back to top ]