粒子シミュレーション装置及び粒子シミュレーション方法
【課題】粒子径分布を有する粉粒体のＤＥＭ計算のメモリ使用量を抑制する。
【解決手段】粒子群の位置・速度情報を含む粒子情報を保持する粒子情報保持部１１と、粒子の位置に応じた順序で特定する粒子番号を各粒子に付与する粒子番号変更部１４と、着目粒子と接触している可能性がある他の粒子との接触候補ペアを選択する接触候補リスト作成部１６と、接触候補ペアの粒子間の接触力を粒子情報に基づき算出して接触力テーブルに格納する接触判定部１８と、粒子ごとに、当該粒子より粒子径の大きい粒子との接触力を接触力参照表５４を用いて接触力テーブルから抽出し、粒子径の小さい粒子との接触力を、積算接触候補数s_jgi[i]を用いて接触力テーブルの格納位置を特定して抽出し、接触力を総和演算する接触力計算部１９と、粒子ごとの接触力に基づいて粒子情報を更新する粒子情報更新部２０と、を備える。
【解決手段】粒子群の位置・速度情報を含む粒子情報を保持する粒子情報保持部１１と、粒子の位置に応じた順序で特定する粒子番号を各粒子に付与する粒子番号変更部１４と、着目粒子と接触している可能性がある他の粒子との接触候補ペアを選択する接触候補リスト作成部１６と、接触候補ペアの粒子間の接触力を粒子情報に基づき算出して接触力テーブルに格納する接触判定部１８と、粒子ごとに、当該粒子より粒子径の大きい粒子との接触力を接触力参照表５４を用いて接触力テーブルから抽出し、粒子径の小さい粒子との接触力を、積算接触候補数s_jgi[i]を用いて接触力テーブルの格納位置を特定して抽出し、接触力を総和演算する接触力計算部１９と、粒子ごとの接触力に基づいて粒子情報を更新する粒子情報更新部２０と、を備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子シミュレーション装置及び粒子シミュレーション方法に関する。
【背景技術】
【０００２】
土木工学の研究分野から端を発したＤＥＭ（Discrete Element Method：離散要素法）は、その後粉体工学へ応用され、昨今では、粒子性物質のみならず連続体として扱いが困難な非常に多種多様の問題に対する数値解析手法として益々需要は広がってきた。一方、離散的な粒子運動を数値的に追跡する手法には、ＭＤ（Molecular Dynamics）、ＬＧＡ（Lattice Gas Automaton）、ＳＰＨ（Smoothed Particle Hydrodynamics）などＤＥＭ以外にも数多く開発されてきている。しかし、数ある数値粒子解析手法の中で、ＤＥＭの応用範囲が広い理由は、粒子の大きさや形状という個々の粒子の幾何情報を考慮しながら粒子間相互作用力を求めニュートンの運動方程式を解くという計算を非常に高速に行えることである。勿論、ＤＥＭ以外にもＤＤＡ（Discontinuous Deformation Analysis）なども、大きさと形状を考慮することが可能な手法であるが、もともと計算負荷が高いので、実用的にあまり多数の粒子の扱うことができない。
【０００３】
ところで、砂や小麦粉のように我々が目にする実際の粒状体や粉体（以下、両者をまとめて「粉粒体」と呼ぶ）は、篩い分けによる粒度調整など特別な処理を施さない限り、一般的には様々な大きさの粒子が混在した状態で存在している（以下、これを「粒子径分布」がある状態と呼ぶ）。実はこの粒子径分布は、粉粒体の均質性や変形や流動といった力学的な性質を支配する大きな要因であり、粉粒体を扱う産業では粒径分離効率や粒度偏析による粉体層の不均質化および混合効率、輸送・充填の制御にはいつも問題となる。したがって、粒子径分布が粒子群全体の振る舞いに与える影響を正しく理解することは、粉粒体を扱うためには必要不可欠である。
【０００４】
ＤＥＭは、前述したように粒子の大きさを考慮できる手法であるから、粒子径分布がもたらすマイクロメカニックスをＤＥＭで調べたくなる気持ちはごく自然である。しかし、粒子径が１／ｎになると粒子１個の体積が(１／ｎ)３になり、もとのサイズの粒子１個分の体積を充填するためにｎ３個の粒子が必要となる（厳密には空隙が生じるため、ｎ３個よりは少なくなるがオーダーは変わらない）ため、実は、計算機のメモリ量や計算時間の制約上、あまり大きな粒子径分布は扱えない。例えば自然土には数センチの砂利から数ミクロンの粘土粒子まで粒子径は約５桁程度に渡って分布している。すると、全てのサイズで同程度の体積を占めている場合、最大粒径の粒子数を１個にしても、全体では(１０５)３個オーダーの粒子数となってしまう。
【０００５】
また、ＤＥＭで粒子径分布を扱う場合には、上述の粒子数の問題に加えて、粒子径の違いによる接触点数のアンバランスも問題となる。例えば、最大粒子と最小粒子の粒径比をｒ：１とした時、最大粒子の表面を最小粒子が埋め尽くす場合と、最小粒子の表面を最大粒子が埋め尽くす場合を考えると、表面を覆う粒子の個数はそれぞれ４(ｒ＋１)２／１２および４(ｒ＋１)２／ｒ２となる。すなわち、接触点数のアンバランスは最大でｒ２：１のオーダーになる。個々の粒子の接触点数のアンバランスは、粒子ごとの接触力の計算回数と接触の総和計算量のアンバランスにつながり、並列計算には大きな障壁となる。つまり、ＤＥＭでは、均一粒径の場合と違って粒径比が大きくなると、単純に粒径の幾何学的な情報だけからでも粒子数が増加し接触点数のアンバランスが生じることが見込まれ、結果としてかなり計算負荷が高くなることが分かる。このことに加えて、実際のＤＥＭのプログラム内部の処理においても、粒子径分布が存在すると、接触判定の対象となる接触候補の絞り込み検索と接触候補ペアを格納するための処理が、均一粒径の場合と比較して計算負荷が著しく高くなってしまう。
【０００６】
この粒子径分布がある問題に対して、接触候補ペアの抽出と格納に関連する処理を高速化するためのＤＥＭプログラミングの手法は幾つか提案されてきている（例えば特許文献１）。
【０００７】
一方、最近ではＣＰＵと比べて価格に対する演算性能の比が格段に大きいＧＰＵ（Graphics Processing Unit）を用いた数値シミュレーションの高速化が盛んに行われている。これに伴い、粒子系シミュレーションについても、高価なスーパーコンピュータを用いたベクトル化処理よりも、ベクトル演算と同等の超並列演算が可能で安価なＧＰＵに対する高速化アルゴリズムが提案されてきている（例えば特許文献２）。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２００７−２８６５１４号公報
【特許文献２】特開２００９−６９９３０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかし、特許文献１など従来のＤＥＭ高速化手法は、まずセルに粒子を格納してからセルサイズやセル探索範囲を最適化する手法であって、並列化やベクトル化による大規模高速化をことさら意識したものではない。具体的な例を挙げると、粒子番号のループで粒子をセルに格納する処理をベクトル化のようなＳＩＭＤ（Single Instruction Multiple Data）型並列処理で行うと、２つ以上の粒子が同一セルに格納されるときにはメモリ書き込み競合が生じるために非常に都合が悪い。このメモリ書き込み競合問題は、単一粒径成分の場合でも一つのセルに複数の粒子が入ってしまう場合に発生するが、粒子径分布がある場合は一つのセルに入る粒子数自体に大きな分布ができるので、並列化やベクトル化によるプログラムの高速化の際には、さらに厄介な大きな問題となる。
【００１０】
また、特許文献２などのＧＰＵ高速化アルゴリズムでは、接触候補粒子ペアの検出処理においては、ＣＰＵベースのプログラムと同様にセルに粒子を格納してからセルを探索する方法が用いられている。しかし上述したように、ＧＰＵでも複数粒子が同一セルに格納される場合には、やはりメモリ書き込み競合が生じるため、この問題を回避するために、幾つかの対策が提案されている。一つは、粒子の格納と確認操作を複数回繰り返す方法で、他方は、コンパイラ付属の排他処理関数（Atomic関数）を用いる方法である。また、粒子間接触力の計算ルーチンにおいてもメモリ書き込み競合を防ぐために、粒子ｉ‐ｊ間の接触力計算を粒子ｉと粒子ｊそれぞれに対して同じ計算を二度行うことで問題を回避する工夫がなされている。
【００１１】
しかし、前述のように、粒子径分布があると、セルによって格納される粒子数に大きな偏りが生じやすく、また大粒子に対する接触粒子数は小粒子に対するそれと比べて非常に多いために、ロードバランスを大きく損ないやすい。例えば、これらの方法を用いた三次元ＤＥＭで二成分系粒子群に対するＧＰＵ計算を行うと、二成分の粒子径比が１：４になるだけで単一成分の場合と比べて計算時間は３倍近くにもなることが報告されている。
【００１２】
そして、粒子径分布を考慮した粉粒体のＤＥＭシミュレーションを行う場合、粒子径を均一とした同じ粒子数の場合に比べて、接触力を考慮すべき近傍粒子の数が急激に増大するので、ＤＥＭ計算を実現するために確保しなければならないメモリ使用量が増大するという問題があった。
【００１３】
ＧＰＵ演算は、ＣＰＵ演算のようにメモリをハードウェア的に増設することができないため、あくまでグラフィックカードに搭載された限られたメモリ容量内でしか処理できない。勿論、ＧＰＵメモリの不足分を転送通信によってＣＰＵメモリに頼る方法も考えられるが、ＧＰＵメモリとＣＰＵメモリ間の転送速度は非常に遅いため、せっかくのＧＰＵの高い計算効率が全く生かされず推奨できない。このように、ＧＰＵをＤＥＭ計算に使用する場合、現状のＧＰＵのメモリ容量が最大でも４ＧＢ程度であるので、このメモリ使用量の問題が、ＧＰＵを粒子シミュレーションに適用することの制約となっている。そして、この問題は、ＧＰＵだけでなく、ＧＰＵと同様の共有メモリ型のＳＩＭＤ（Single Instruction Multiple Data）型並列処理を行うハードウェア（例えばベクトルプロセッサ）（以下、「共有メモリ型並列計算機」という）にも該当する問題である。
【００１４】
本発明は、上記課題を解決するためになされたものであり、粒子径分布を有する粉粒体のＤＥＭ（離散要素法）計算を共有メモリ型並列計算機により実行する際に、メモリ使用量を抑制することができる粒子シミュレーション装置及び粒子シミュレーション方法を提供することを目的とする。
【課題を解決するための手段】
【００１５】
上記課題を解決するため、本発明に係る粒子シミュレーション装置は、作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力する入力手段と、粒子情報を保持する粒子情報保持手段と、粒子のそれぞれについて、位置情報に応じた順序で特定する特定情報を付与する付与手段と、粒子の中から前記特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択手段と、選択候補選択手段により選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の位置情報及び速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出手段と、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を前記接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成手段と、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出手段と、粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を接触力参照表の参照情報を用いて接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、積算接触候補数を用いて接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算手段と、接触力総和演算手段により計算された粒子ごとの接触力に基づいて、粒子の新たな位置情報及び速度情報を算出する算出手段と、算出手段により算出された新たな位置情報及び速度情報を用いて、粒子情報保持手段に保持される粒子情報を更新する粒子情報更新手段と、を備えることを特徴とする。
【００１６】
同様に、上記課題を解決するため、本発明に係る粒子シミュレーション方法は、作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置における粒子シミュレーション方法であって、粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力し、この粒子情報を粒子情報保持手段に保持する入力ステップと、粒子のそれぞれについて、位置情報に応じた順序で特定する特定情報を付与する付与ステップと、粒子の中から特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択ステップと、接触候補選択ステップにおいて選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の位置情報及び速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出ステップと、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成ステップと、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出ステップと、粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を接触力参照表の参照情報を用いて接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、積算接触候補数を用いて接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算ステップと、接触力総和演算ステップにおいて計算された粒子ごとの接触力に基づいて、粒子の新たな位置情報及び速度情報を算出する算出ステップと、算出ステップにおいて算出された新たな位置情報及び速度情報を用いて、粒子情報保持手段に保持される粒子情報を更新する粒子情報更新ステップと、を含むことを特徴とする。
【００１７】
このような粒子シミュレーション装置及び粒子シミュレーション方法によれば、作業空間内の粒子のそれぞれについて、当該粒子の近傍に位置する、この粒子と接触する可能性の高い他の粒子のうち当該粒子の粒子径より小さい粒子については、粒子の数を示す積算接触候補数のみを保持しておき、この積算接触候補数を利用して接触力テーブルの格納位置を特定して当該粒子と他の粒子との接触力を抽出する。これにより、接触力参照表には、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子に関する情報を含める必要がなく、当該粒子の粒子径より大きい粒子に関する情報のみを保持すればよくなる。粒子径分布を有する粒子群をシミュレーション対象とする場合、自分より粒子径の小さい粒子まで接触を考慮すると、同時に接触可能な粒子の数が増大するため、予め接触力参照表のために確保すべきメモリ領域を多くしなければならない。これに対し、本発明の上記構成により、自分より粒子径の大きい粒子との接触のみを考慮すればよく、接触を考慮する粒子数を減らすことができ、接触力参照表を作成するために予め割り当てるメモリ領域を少なくすることが可能となり、この結果、メモリ使用量を抑制することができる。このようにメモリ使用量が抑制されることにより、ＧＰＵなどの共有メモリ型並列計算機を利用して粒子シミュレーションを実行することが可能となる。
【００１８】
また、粒子シミュレーション装置では、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、粒子情報には、位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、付与手段が、粒子のそれぞれについて、セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することが好適である。
【００１９】
同様に、粒子シミュレーション方法では、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、粒子情報には、位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、付与ステップが、前記粒子のそれぞれについて、セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することが好適である。
【００２０】
この構成により、各粒子にセル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報が付与されるので、各粒子がセル番号順、さらには同一セル内では粒子径順にソーティングされる。このため、接触力参照表を作成する際に、特定情報に基づき粒子径が小さい粒子を容易に認識することができ、粒子径が小さい粒子については接触力参照表にいれるか否かの判定を行う必要がなくなる。この結果、接触力参照表の作成において計算効率を向上させることができる。
【００２１】
また、粒子シミュレーション装置では、接触候補選択手段が、粒子郡の中から特定情報の昇順で着目粒子を選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との組と、この組に昇順で付される接触候補番号とからなる接触候補ペア情報を、作業空間内の全粒子について含む接触候補リストを作成し、接触力算出手段が、接触候補リストに含まれる接触候補ペア情報のそれぞれについて、該接触候補ペア情報に関する２つの粒子の間の接触力を、当該粒子の前記位置情報及び速度情報に基づき算出し、算出した接触力と接触候補番号とを関連付けた接触力情報を接触力テーブルに格納し、接触力参照表作成手段が、粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の大きい他の粒子との接触力を前記接触力テーブルから参照するための参照情報として、該他の粒子及び当該粒子に関する接触候補番号を、該他の粒子及び当該粒子と関連付けて保持する接触力参照表を作成し、積算接触候補数算出手段が、粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の小さい他の粒子の数を示す接触候補数を前記特定情報順で積算した値である積算接触候補数s_jgi[i]（ｉは当該粒子の特定情報を示す）を算出し、接触力総和演算手段が、粒子のそれぞれについて、接触力参照表の参照情報を用いて接触力を接触力テーブルから抽出し、また、積算接触候補数s_jgi[i]を用いて接触候補番号がs_jgi[i-1]〜s_jgi[i]‐1の接触力を接触力テーブルから抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算することが好適である。
【発明の効果】
【００２２】
本発明に係る粒子シミュレーション装置及び粒子シミュレーション方法によれば、粒子径分布を有する粉粒体のＤＥＭ（離散要素法）計算を共有メモリ型並列計算機により実行する際に、メモリ使用量を抑制することができる。
【図面の簡単な説明】
【００２３】
【図１】本発明の一実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図２】本実施形態で粒子モデルに適用するＶｏｉｇｔモデルを示す概略図である。
【図３】本実施形態における作業領域を示す図である。
【図４】粒子番号変更の流れを示す図である。
【図５】作業領域における粒子番号を変更する処理と、各セルに粒子番号の最大値及び最小値を記憶する処理を示す図である。
【図６】近傍粒子表の構成の一例を示す図である。
【図７】接触候補リストの構成の一例を示す図である。
【図８】現在の接触候補リストのペア粒子について１ステップ前にペアであったか調べる処理を示す図である。
【図９】接触力参照表の構成の一例を示す図である。
【図１０】接触力参照表及び積算接触候補数と、接触力ボックスとの対応関係を示す図である。
【図１１】本実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【図１２】粒子モデルの最大粒子径と最小粒子径との比を変化させた場合の本実施形態及び比較例のメモリ使用量を示し図である。
【発明を実施するための形態】
【００２４】
以下、図面を参照しながら本発明の実施形態を詳細に説明する。なお、図面の説明において同一又は同等の要素には同一の符号を付し、重複する説明を省略する。
【００２５】
図１は、本発明の一実施形態に係る粒子シミュレーション装置１０の機能ブロック図である。図１に示すように、粒子シミュレーション装置１０は、粒子情報保持部（入力手段、粒子情報保持手段）１１、粒子情報取得部１２、接触候補リスト更新判定部１３、粒子番号変更部（付与手段）１４、近傍粒子表作成部１５、接触候補リスト作成部（接触候補粒子選択手段、積算接触候補数算出手段）１６、接触力参照表作成部（接触力参照表作成手段）１７、接触判定部（接触力算出手段）１８、接触力計算部（接触力総和演算手段）１９、粒子情報更新部（算出手段、粒子情報更新手段）２０を備えている。
【００２６】
粒子シミュレーション装置１０は、物理的には、ＧＰＵ（Graphics Processing Unit）、主記憶装置であるＲＡＭ（Random AccessMemory）及びＲＯＭ（ReadOnly Memory）、ハードディスク装置等の補助記憶装置、入力デバイスである入力キー等の入力装置、ディスプレイ等の出力装置、ネットワークカード等のデータ送受信デバイスである通信モジュールなどを有するコンピュータとして構成されている。図１を参照して以下に説明する粒子シミュレーション装置１０の各機能は、ＧＰＵ、ＲＡＭ等のハードウェア上に所定のコンピュータソフトウェアを読み込ませることにより、ＧＰＵの制御のもとで入力装置、出力装置、通信モジュールを動作させるとともに、ＲＡＭや補助記憶装置におけるデータの読み出し及び書き込みを行うことで実現される。
【００２７】
ここで、まず本実施形態で適用する粒子モデルについて説明する。本実施形態の粒子シミュレーション装置１０が対象とするのは、粒子径分布を有する粉粒体群のＤＥＭシミュレーションである。つまり、様々な大きさの粒子径をもつ粒子が混在した状態の粒子モデルが適用される。
【００２８】
そして、このような粒子モデルにおいて、各粒子間に発生する粒子間接触力の法線方向および接線方向成分Ｆｎ，Ｆｔは、図２に示すＶｏｉｇｔモデル７１を用いて以下の式により表される。
【数１】
ここでｘとｖは相対変位および相対速度ベクトル、Ｋは線形バネ７２における弾性係数、ηはダッシュポット７３の粘性減衰係数であり、それぞれの添字ｎ，ｔにより法線方向成分及び接線方向成分を表す。μtは滑り摩擦係数、ｍｉｎは絶対値が小さい方の値を表す関数である。粘性減衰係数ηは、反発係数eと次式の関連がある。
【数２】
ここでｍｉ，ｍｊは、それぞれ粒子ｉ，ｊの質量を表す。添字ｑは、法線方向成分の場合はｑ＝ｎとなり、接線方向成分の場合にはｑ＝ｔとなる。
【００２９】
したがって、ＤＥＭで扱う並進および重心を中心とする回転の運動方程式は以下の式で表される。
【数３】
ここで、ｕｐとωｐは粒子の速度および角速度、ｍｐとＩｐは粒子の質量および慣性モーメント、Ｒは粒子中心から接触点へ向かう位置ベクトルであり、添字ｐは、粒子ｉの場合はｐ＝iとなり、粒子ｊの場合はｐ＝ｊとなる。また、（５）式の右辺第二項は転がり摩擦力を表し、μｒは転がり摩擦係数、ｂは接触面幅である。（１）、（２）式のＶｏｉｇｔモデル７１は、後述する接触判定部１８で粒子間接触力を算出するのに用いられ、また（４）、（５）式の運動方程式は、粒子情報更新部２０で粒子の座標を算出するのに用いられる。
【００３０】
また、本実施形態において粒子が運動する領域である作業領域５１は、図３に示すように、一辺の大きさがcsizeの立方体セルで区分される。セルサイズcsizeは、パラメータφ (0≦φ≦ 1) を用いて以下の式で与える。
csize＝2.0{(1.0−φ)rmax＋φrmin)}(1.0＋α) （６）
ここで、rmaxは作業領域に配置される粒子モデルの最大粒子径であり、rminは最小粒子径であり、αは更新頻度を調整するパラメータであり例えばα＝０．２である。また、作業領域５１内の各セルにはセル番号cell_idが付されている。セル番号は、図３に示すように一次元であり、例えば次式のように表される。
cell_id＝i.x＋id.x×i.y＋id.x×id.y×i.z （７）
i.m＝{0, 1, 2,…, id.m-1} （m＝x, y, z）
ここでcell_idはセル番号であり、i.mはｍ方向のセル位置であり、id.mはｍ方向のセル数である。つまり、粒子シミュレーション装置１０は、シミュレーション条件として粒子の最大粒子径rmax、最小粒子径rmin、パラメータφ，αが設定されると、作業領域５１のセルサイズcsizeを自動的に算出して設定し、さらに画定したセルにセル番号を割り当てることができる。
【００３１】
粒子情報保持部１１は、作業空間５１内の複数の粒子のそれぞれについて粒子情報を保持している。粒子情報は、座標（位置情報）、並進速度（速度情報）、回転速度（速度情報）、粒子半径、質量、粒子番号、セル番号を含む。座標とは、図３に示す三次元の作業空間５１における三次元座標（pos.x，pos.y，pos.z）である。粒子番号は各粒子を識別するための番号である。セル番号は当該粒子が所属するセルのセル番号である。なお、シミュレーションを開始する際には、粒子情報保持部１１は、シミュレーションに適用する粒子径分布に応じて各粒子の粒子半径、質量を設定し、座標、速度、粒子番号を各粒子にランダムに割り当てる。セル番号には、粒子ごとにランダムに割り当てられた座標に基づき特定される、この粒子が属するセルの番号を付与する。そして、このように生成した各情報を初期状態の粒子情報として入力し保持する。また、シミュレーションの実行中には、後述する粒子情報更新部２０により粒子情報の各情報が更新されると、粒子情報保持部１１は、保持している粒子情報を、粒子情報更新部２０により更新された新しいものに置き換える。
【００３２】
粒子情報取得部１２は、粒子情報保持部１１から粒子情報を取得する。粒子情報取得部１２は、粒子シミュレーションの開始時に粒子情報保持部１１から予め設定されている初期状態の粒子情報を取得し、シミュレーション動作中には、粒子情報保持部１１が更新される度に新たな粒子情報を粒子情報保持部１１から取得する。粒子情報取得部１２は、取得した粒子情報を接触候補リスト更新判定部１３に送信する。
【００３３】
接触候補リスト更新判定部１３は、後述する接触候補リスト５３を更新するか否かを判定する。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。また、粒子の積算移動距離は、例えば、粒子情報取得部１２が新しい粒子情報を粒子情報保持部１１から取得する度に、１ステップ前の粒子情報との座標データの差分からこの直近の１ステップにおける粒子の移動距離を算出し、これを粒子ごとに積算して求める。接触候補リスト更新判定部１３は、１ステップ前の粒子情報と、パラメータαと、粒子の積算移動距離とを記憶しており、粒子情報取得部１２から新たな粒子情報を受信すると、積算移動距離を更新し、この更新された積算移動距離がｒαより大きい粒子の有無を判定する。そして、積算移動距離が所定値以上の粒子が有る場合に、接触候補リストを更新するものと判定し、新たな粒子情報を粒子番号変更部１４に送信する。なお、粒子移動距離の積算値は、接触候補リスト５３の更新が行われるたびに初期化される。また、更新された積算移動距離がｒαより大きい粒子が無い場合には、接触候補リストを更新しないものと判定し、粒子番号変更部１４には粒子情報を送信せず、接触判定部１８に新たな粒子情報を送信する。
【００３４】
粒子番号変更部１４は、作業領域５１内の各粒子のそれぞれについて、セル番号に応じた順序で特定し（すなわち位置情報に応じた順序で特定し）、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した粒子情報（特定情報）を付与する。言い換えると、粒子番号変更部１４は、作業領域５１内の各粒子が所属するセル及びその粒子径に基づき、各粒子の粒子番号を変更して、粒子番号をソーティングする。
【００３５】
粒子番号変更部１４は、接触候補リスト更新判定部１３から粒子情報を受信すると、図４に示すように、各粒子の粒子番号に、この粒子が所属するセルのセル番号と、この粒子の粒子径とを関連付けて1組の情報セットとし、これらの情報セットを粒子番号順に記憶する（図４（ａ））。次に、これらの情報セットをセル番号が昇順になるように並べ換え（図４（ｂ））、セル番号順に粒子番号を付け直す（図４（ｃ））。続いて、同一のセル番号の粒子群ごとに粒子径を比較し、これらの粒子群の中で粒子径の大きい順に並べ換え（図４（ｄ））、粒子径順に粒子番号を付け直す（図４（ｅ））。
【００３６】
この粒子番号を変更する処理により、例えば図５（ａ）に示すように作業空間５１に配置された粒子の粒子番号が、図５（ｂ）に示すように、セル番号及び粒子径に沿った粒子番号に変更される。粒子番号変更部１４は、粒子番号の変更が反映された粒子情報を近傍粒子表作成部１５に送信する。
【００３７】
近傍粒子表作成部１５は、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２を作成する。近傍粒子表作成部１５は、着目粒子の所属するセル、およびこのセルに隣接するセルや隣接セルのさらに外側に隣接する近傍セルに存在する粒子のうち、自身より粒子径が大きい粒子、或いは粒子径が同じ場合は初期配置時の粒子番号が自身より大きい粒子で、粒子間距離がある値以下の粒子を「近傍粒子」とする。以下に詳細を説明する。
【００３８】
近傍粒子表作成部１５は、近傍粒子表５２を作成するためのメモリ領域を予め確保している。確保するメモリ領域は、各粒子と接触する可能性がある近傍粒子の個数を考慮して、各粒子の近傍粒子表ごとに同一の大きさである。
【００３９】
近傍粒子表作成部１５は、粒子番号変更部１４から粒子情報を受信すると、同じセル番号をもつ粒子の中で、セル内の粒子の番号のうち、最小値Pnum_in_cellminと最大値Pnum_in_cellmaxをセルに記憶させる（図５（ｃ））。図５（ｃ）では、各セルの左上に最大粒子番号、右下に最小粒子番号が記憶されている。つまり、セルcell_idごとに、このセルに含まれる粒子の番号の最小値Pnum_in_cellmin[cell_id]及び最大値Pnum_in_cellmax[cell_id]の情報を格納する配列のみを保持しておく。これにより、各セル番号cell_idのセルには、Pnum_in_cellmin[cell_id]番からPnum_in_cellmax[cell_id]番までの粒子が粒子径の大きい順に所属していることがわかる。
【００４０】
番号cell_idのセルに隣接するセルの番号nei_cell_idは次式で表される。
nei_cell_id[n]= cell_id + con[n] （８）
ここで、conはid.xとid.yによって決まる定数である。セル番号で粒子番号をソートしてあるため、各隣接セルにはPnum_in_cellmin[nei_cell_id]からPnum_in_cellmax[nei_cell_id]までの番号の粒子が連続して存在することがわかる。
【００４１】
次に、近傍粒子表作成部１５は、各粒子ごとに、図６に示すような近傍粒子を格納する箱（配列）nei_pn[i][box_id]を用意し、粒子ｉの周囲の探索範囲（定義は後述する）に含まれる粒子のうち
（ａ）粒子ｉとの粒子間距離がdis_limよりも小さく、
（ｂ）粒子ｉの径よりも大きく、
（ｃ）粒子ｉの径と同じ場合は、初期配置時の粒子番号（i=p_id）が粒子ｉよりも大きい
という条件を満たす粒子を番号（box_id）が小さい方の箱から詰めて格納させる。ここで、dis_lim = (rｉ + rｊ)(1．0+α)であり、rｉ，rｊは近傍粒子ペアｉ，ｊの粒子半径である。αは接触候補の更新回数を調整するパラメータであり、例えばα＝０．２である。rｉ，rｊは粒子番号変更部１４から受信した粒子情報に含まれ、αは近傍粒子表作成部１５が予め保持している値であり、近傍粒子表作成部１５は、これらのデータを用いてdis_limを算出する。
【００４２】
上記条件（ａ）〜（ｃ）を満たす粒子を、粒子iに対する「近傍粒子」と呼ぶ。図６に示す一連の箱の集合を、作業空間５１内の全ての粒子ｉについて用意したものを、本実施形態では近傍粒子表５２という。近傍粒子表５２は、粒子i及びこの粒子iの近傍粒子の数box_idに関する二次元配列nei_pn[i][box_id]である。
【００４３】
この方法によって、着目粒子iについて、この粒子i以上の粒子径を持った粒子ｊのみが近傍粒子表５２に記憶されるため、近傍粒子表５２に必要なメモリ使用量の増加を抑制することができる。
【００４４】
ここで、近傍粒子表作成部１５は、粒子ｉの近傍粒子を決定するために、粒子ｉから(rmax＋ri)(1.0＋α)の距離にある粒子を探索する必要がある。粒子ｉが属するセルに隣接する２６個のセルのほかに、粒子モデルの最大粒子の半径ｒmaxに応じてさらに外側のセルまで探索範囲を広げることとなる。
【００４５】
近傍粒子表作成部１５が探索する近傍セルの範囲は、上記（６）式で算出されたセルサイズcsizeと、作業領域内の最大粒子の粒子径を示す最大粒子半径ｒmaxと、着目粒子半径ｒｉから次式によって求められる。
nei_cell＝2・ceil{2.0(rmax＋ri)(1.0＋α)／csize}−1 （９）
ceilは小数点以下を切り上げる関数であり、近傍粒子表作成部１５に予め保持されている。nei_cellは粒子ｉが所属するセルからｘ、ｙ、ｚ方向それぞれに対して＋−方向に探索するセル数を表す。
【００４６】
これにより、セルサイズが小さくなるにつれて探索するセル数が急激に増えて計算負荷が増加する懸念が生じる。しかし、隣接する２６個のセルより外側のセルに対しては、粒子径でソートした粒子の並び順を利用して近傍粒子の探索を効率化することができる。つまり、各セルにはPnum_in_cellmin [cell id]番からPnum_in_cellmax [cell id]番までの粒子が粒子径の大きい順に所属しており、各セル内の粒子が近傍粒子であるか否かを調べる際も粒子径の大きい順に調べて行くことになる。この時、粒子iが所属するセルとその近傍セルの位置関係から上記の（ａ）の条件を満たさない粒子径が求まるのでPnum_in_cellmax [cell id]番までの全ての粒子間距離を調べる必要はなく、（ａ）の条件を満たさない粒子径の粒子番号になったら探索を中止することができる。ここで、粒子iに対して(a)の条件を満たさない粒子径rj とは以下の式を満足するもののことである。
2 ceil{2.0(rj + ri )(1.0 + α)／csize}−1 ＜ dis cell （１０）
dis cell は粒子iの所属するセルから粒子ｊの所属するセルまでのセル数を表し、ｘ、ｙ、ｚ方向のセル数で最小のものを表す。
【００４７】
さらに、近傍粒子表作成部１５は、粒子ｉに対する近傍粒子数n_jli[i]と，上記条件（ａ）を満たすが（ｂ）及び（ｃ）を満たしていない粒子数n_jgi[i]を、それぞれ一次元配列に記憶させておく。
【００４８】
図６の例では、例えば４番の粒子（ｉ＝４）に対してセルを検索した結果５，６，７，８番の粒子が（ａ）の条件を満たし，そのうちの近傍粒子が７，８番で、５，６番の粒子は（ｂ），（ｃ）の条件を満たさない場合、
nei_pn[4][0]=7
nei_pn[4][1]=8
n_jgi[4]=2
n_jli[4]=2
となる。
【００４９】
このように、近傍粒子表作成部１５は、粒子番号変更部１４から粒子情報を受信すると、着目粒子ｉごとに、探索範囲nei_cellのセルを画定し、この探索範囲内に含まれる他の粒子のうち上記（ａ）〜（ｃ）の条件を満たす粒子の粒子番号を近傍粒子表５２の配列nei_pn[i][box_id]に格納する。さらに、近傍粒子数を配列n_jli[i]に格納し、上記条件（ａ）を満たすが（ｂ）及び（ｃ）を満たしていない粒子数、すなわち探索範囲に存在するものの粒子径が着目粒子より小さく、近傍粒子表５２に格納されなかった粒子の数を配列n_jgi[i]に格納する。そして、これらnei_pn[i][box_id] 、n_jli[i]、n_jgi[i]の情報を、粒子情報と共に接触候補リスト作成部１６に送信する。
【００５０】
本実施形態では、近傍粒子表５２には、着目粒子の粒子径以上の粒子のみを近傍粒子として保持している。複数の粒子径をもつ粒子群からなる多成分系では大粒子に対して接触候補になる小粒子の数が膨大になる可能性が有る。つまり、着目粒子より粒子径の小さい粒子とは同時に接触可能な数が非常に多くなる。このため、着目粒子より粒子径の小さい粒子まで近傍粒子に含むとすると、予め近傍粒子表のために用意する配列の数も最大値を考慮して大きくせざるを得ない。その場合、近傍粒子表５２の配列に必要なメモリ量が膨大になり、ＧＰＵのグローバルメモリを食い潰す危険が有る。本実施形態では、自分より大きい粒子しか近傍粒子表５２に格納しないので、近傍粒子表５２のために用意される配列のメモリ使用量を抑制することができる。
【００５１】
接触候補リスト作成部１６は、作業空間５１の全粒子の中から粒子番号（特定情報）の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する。言い換えると、接触候補リスト作成部１６は、作業空間５１の全粒子について、接触している可能性がある粒子ペアからなる接触候補リスト５３を作成する。接触候補リスト５３とは、具体的には、接触している可能性がある２つの粒子ｉ，ｊの粒子番号を格納する一組の一次元配列list_i[Lnum]、list_j[Lnum]である。ここでLnumは接触粒子ペアごとに付される接触候補番号である。
【００５２】
接触候補リスト作成部１６は、近傍粒子表作成部１５から情報（近傍粒子表５２nei_pn[i][box_id] 、近傍粒子数n_jli[i] 、n_jgi[i]、及び粒子情報）を受信すると、まず、図７に示すように、n_jgiのプレフィックス和s_jgiを次式から求める。
【数４】
【００５３】
このプレフィックス和s_jgi[i]は、粒子iの近傍粒子表５２に格納されなかった粒子数、すなわち粒子ｉの周囲の所定の探索範囲（上記（９）式で画定される近傍のセル）に存在する粒子のうち粒子ｉの粒子径より小さい粒子の数n_jgi[i]（接触候補数）の、粒子番号ｉまでの積算値をあらわす一次元配列であって、本実施形態では「積算接触候補数」とよぶ。言い換えると、積算接触候補数s_jgi[i]は、接触候補数n_jgi[i]を粒子番号ｉ順で積算した値である。
【００５４】
接触候補リスト作成部１６は、次に、粒子番号の昇順で着目粒子ｉを順次選択し、各粒子ｉに対して、上記（９）式で画定された粒子ｉの周囲の探索範囲のセルを再度探索して、粒子ｉの周囲の探索範囲に含まれる粒子のうち
ａ）粒子ｉとの粒子間距離がdis_limよりも小さく
ｂ）粒子ｉの径よりも小さく、
ｃ）粒子ｉの径と同じ場合は、初期配置時の粒子番号が粒子ｉよりも小さい
という条件を満たす粒子ｊを抽出する。そして、これらの粒子ｉ，ｊのペアに接触候補番号Lnumを昇順で付与し、それぞれの粒子の粒子番号を１次元配列list_i，list_jに格納して、接触候補リスト５３を作成する。
【００５５】
この条件ａ）〜ｃ）を満たすこととは、上述の近傍粒子表に格納する近傍粒子を抽出する条件（ａ）を満たし（ｂ），（ｃ）を満たさないのと同等である。つまり近傍粒子表５２に格納されない粒子ペア（例えばｉ＝４の着目粒子の場合にはｉ＝５，６の粒子とのペア）が接触候補ペアである。
【００５６】
例えば、図７の例では、接触候補リスト５３の内容は、以下の表１のようになる。
【表１】
つまり、接触候補リスト５３は、作業空間５１の全粒子について、接触している可能性がある粒子ペアlist_i，list_jに番号Lnumを付加した接触候補ペア情報を含むものである。
【００５７】
接線方向の接触力を求める時に必要となるバネの伸縮量は積算値で与えられるため、前ステップで接触していて現ステップでも接触している場合は接触候補ペアリストを更新する時に接線方向バネの伸縮量を引き継ぐ必要が有る。そこで、次に、１ステップ前の接触候補ペアBlist_i，Blist_jを探索し、現在の接触候補ペアlist_i, list_jが１ステップ前にもペアであったか調べる。ペアであったら、図２に示した粒子モデルにおける法線ベクトルと接線方向のバネによる力（バネの伸縮量）を新しい候補リスト番号の配列に格納する。
【００５８】
ここで、１ステップ前及び現在の各パラメータ及び配列の表記は以下のとおりである。
【表２】
【００５９】
粒子番号i=list_i[Lnum], j=list_j[Lnum]の接触候補ペアが有るとき、このペア粒子の１ステップ前の粒子番号はBi=Bpn[i]，Bj=Bpn[j]である。Bs_jgi[Bi−1]〜Bs_jgi[Bi]−1の範囲の接触候補Blist_j[BLnum]を調べBjと一致するか調べる。一致したときはそのリスト番号BLnumを用いて１ステップ前の法線ベクトルと接線方向の接触力BSpring[BLnum]を現在のリスト番号配列Spring[Lnum]に代入する。
【００６０】
この処理について、図８を参照して、現在の接触候補リスト番号Lnum = 2におけるｉ＝4，ｊ＝6番のペア粒子について１ステップ前にペアであったか調べる処理を説明する。
【００６１】
まず、ｉ＝4，ｊ＝6の１ステップ前の粒子番号はBpn[]関数を使ってBi=2=Bpn[4]，Bj=3=Bpn[6]と求まる。ここでBiとBjとの関係は、ｉとｊの関係と同様に、接触候補リスト５３の条件ａ）〜ｃ）を必ず満たしている。なぜなら、条件ａ）〜ｃ）はソートによって変化しない粒子径と初期配置時の粒子番号を使用しているためである。したがってBiに対する接触候補の粒子番号はリスト番号BLnumがBs_jgi[Bi−1]〜Bs_jgi[Bi]−1（この例では１〜１）の範囲のBlist_j[BLnum]に格納されている。そこで、このBLnumの範囲でBlist_j[BLnum] がBjと一致するか調べると、BLnum = 1のときにBlist_j[1]= 3でありBj = 3と一致する。
【００６２】
したがって，現在のリスト番号Lnum=2の候補ペアは1ステップ前もペアであった事がわかり、そのリスト番号はBLnum=1である。そしてBLnum=1の法線ベクトルと接線方向の接触力の値をLnum=2の配列に代入する（Spring[2]=BSpring [1]）ことで、1ステップ前の接触状態を引き継ぐ事ができる。
【００６３】
このように、接触候補リスト作成部１６は、近傍粒子表作成部１５から各種情報報（近傍粒子表５２nei_pn[i][box_id] 、近傍粒子数n_jli[i] 、n_jgi[i]、及び粒子情報）を受信すると、着目粒子ｉごとに、上記探索範囲内に含まれる他の粒子のうち上記ａ）〜ｃ）の条件を満たす粒子ｊとのペアに番号Lnumを付与し、これらの粒子ペアｉ，ｊの粒子番号をそれぞれ接触候補リスト５３の配列list_i[Lnum]，list_j[Lnum]に格納する。また、着目粒子ｉの周辺セルに存在し、粒子ｉの粒子径より小さい粒子の数n_jgi[i]の積算値である積算接触候補数s_jgi[i]を算出する。さらに、現在の接触候補ペアlist_i[Lnum]，list_j[Lnum]のうち、１ステップ前も接触候補だったペアを抽出し、これらのペアの１ステップ前の接触力をBSpring[BLnum]から引き継ぎ、Spring[Lnum]に格納する。そして、導出した接触候補リスト５３list_i[Lnum]，list_j[Lnum]と、近傍粒子表５２nei_pn[i][box_id]とを接触力参照表作成部１７に送信し、接触候補リスト５３list_i[Lnum]，list_j[Lnum]、積算接触候補数s_jgi[i]、１ステップ前から引き継いだ接触候補ペアの接触力Spring[Lnum]の情報を、粒子情報と共に接触判定部１８に送信する。
【００６４】
接触力参照表作成部１７は、box_idから対応する候補リスト番号Lnumを算出する関数Ref[i][box_id]を作成する。本実施形態では、この関数をあらわす二次元配列Ref[i][box_id]を接触力参照表５４と呼ぶ。接触力参照表５４は、図９に示すように、接触候補リスト５３のi粒子がｊ粒子から受ける接触力Fij[Lnum]を任意の一粒子ｉ（図９ではi＝４）ごとに参照先（Lnum）（参照情報）をまとめたものである。特に本実施形態では、接触力参照表５４には、着目粒子ｉに関する近傍粒子表５２に格納されている近傍粒子、すなわち着目粒子ｉより粒子径の大きい粒子との接触力の参照先（Lnum）のみが格納されている。ここで、接触力の参照先（Lnum）とは、後述する粒子間接触力の並進成分及び回転成分を保持する２種類の一次元配列Fijt[Lnum]，Fijr[Lnum]及び後述する（１３）式で用いるβ[Lnum]からなる接触力Fij[Lnum]から構成される接触力ボックスから、該当する接触力を取得するために用いる接触候補ペアリスト番号（接触候補番号）Lnumである。
【００６５】
接触力参照表作成部１７は、接触力参照表５４を作成するためのメモリ領域を予め確保している。確保するメモリ領域は、各粒子と接触する可能性がある近傍粒子の個数を考慮して、各粒子の接触力参照表５４ごとに同一の大きさである。
【００６６】
接触力参照表作成部１７は、具体的には、接触候補リスト作成部１６より接触候補リスト５３list_i[Lnum]，list_j[Lnum]及び近傍粒子表５２nei_pn[i][box_id]を受信すると、接触候補リスト５３のLnum番目の粒子ペアである粒子ｉ（＝list_i[Lnum]）及び粒子ｊ（＝list_j[Lnum]）について、粒子ｊの近傍粒子表５２を探索して粒子ｉが記憶されている近傍粒子表番号ｋ（＝box_id）を取得する。上述したように、接触候補ペアの粒子ｊは、粒子ｉより粒子径が小さいので、粒子ｉは粒子ｊの近傍粒子として粒子ｊの近傍粒子表５２に記憶されているからである。そして、粒子ｊに関する接触力参照表５４のｋ番目の配列Ref[j][k]に、粒子ｉ，ｊの接触候補リスト番号Lnumを記憶させる。
【００６７】
なお、上述の接触候補リスト作成部１６による接触候補リスト５３の作成処理、及び接触力参照表作成部１７による接触力参照表５４の作成処理は、以下のアルゴリズムで実施することができる。以下のアルゴリズム中の「接触候補ペア条件==true」とは、接触候補リスト５３の条件ａ）〜ｃ）を満たすことを意味するものである。
【数５】
【００６８】
本実施形態では、接触力参照表５４には、着目粒子の粒子径以上の粒子のみを保持している。このため、上述した近傍粒子表５２と同じ理由により、接触力参照表５４のために用意される配列のメモリ使用量を抑制することができる。
【００６９】
接触力参照表作成部１７は、導出した接触力参照表５４Ref[i][box_id]を接触力計算部１９に送信する。
【００７０】
接触判定部１８は、接触候補リスト５３を参照して、全てのペア粒子について２つの粒子間の接触判定をすると共に、接触していると判定した場合には粒子間の接触力Fijr、Fijtを算出する。
【００７１】
具体的には、接触判定部１８は、接触候補リスト作成部１６から接触候補リスト５３list_i[Lnum]，list_j[Lnum]、積算接触候補数s_jgi[i]、１ステップ前から引き継いだ接触候補ペアの接触力Spring[Lnum]、及び粒子情報を受信するか、または、接触候補リスト更新部１３から粒子情報を受信すると、接触候補リスト５３を走査して接触候補ペアの粒子ｉ，ｊの接触判定を行う。接触判定は、例えば、粒子ｉと粒子ｊとの粒子間距離を計算し、粒子間距離が両粒子の粒子径の和ｒｉ＋ｒｊ以下となる場合に両粒子が接触しているものと判定する。
【００７２】
接触判定において、両粒子が接触していると判定した場合には、上述したＶｏｉｇｔモデルの（１），（２）式を用いて接触力Ｆｎ，Ｆｔを計算する。なお、上記（２）式において現時点のバネによる接線方向接触力はKtxt=spring[]＋ΔFtとして計算する。spring[]は、上述のように１ステップ前のバネによる接線方向接触力であり、ΔFtは、１ステップ前から現時点までに増加（もしくは減少）したバネによる接線方向接触力である。ΔFtは粒子間の相対速度に１ステップの時間Δｔと接線方向成分のバネの弾性係数Ｋｔを乗じて算出される。法線方向成分及び接線方向成分のバネ７２の弾性係数Ｋｎ，Ｋｔと、粒子間の反発係数ｅは、シミュレーション前に設定される所与の値であり、ダッシュポット７３の粘性減衰係数μｎ，μｔは、所与のＫｎ，Ｋｔ，ｅと、粒子情報に含まれる質量ｍｉ，ｍｊを用いて（３）式から導出する。粒子間の相対変位ｘｎ，ｘｔ及び相対速度ｖｎ，ｖｔは、粒子情報に含まれる粒子座標（位置情報）、並進・回転速度（速度情報）に基づき算出する。Ｆｎ，Ｆｔは、このように取得した各変数及びパラメータを用いて、（１）式及び（２）式に示す演算処理を行って算出する。そして、算出したＦｎ，Ｆｔを用いて、接触力の並進成分Fn +Ft と回転成分Ri ×Ft を計算し、これらの成分を、接触候補ペアリスト番号（接触候補番号）Lnumを要素としてもつ配列Fijt[Lnum]，Fijr[Lnum]（上述の接触力ボックスに含まれる配列）に記憶する。
【００７３】
また、接触判定において、接触していないと判定された場合には、接触力は０なので、接触候補ペアリスト番号Lnumを要素としてもつ配列Fijt[Lnum]，Fijr[Lnum]（上述の接触力ボックス）には０を格納する。
【００７４】
この接触判定部１８による接触力算出処理は、以下のアルゴリズムで表現することができる。
【数６】
【００７５】
このように、接触判定部１８は、接触候補リスト５３が作成／更新される毎に、接触候補ペアリスト番号Lnumに対応付けるように接触力ボックスを作成／更新する。また、接触候補リストが更新されない場合でも、粒子情報が更新された場合には、既存の接触候補リストを用いて、接触力ボックスを作成／更新する。
【００７６】
接触判定および接触力は接触候補ペアリスト番号をスレッド化して計算する。ＧＰＵのグローバルメモリはアクセスパターンがスレッドの順である場合に最も良いアクセス効率が得られる。しかし、ペア粒子の番号はリスト番号に対してランダムであるため、粒子情報が格納されているグローバルメモリをランダムにアクセスしメモリアクセス速度が大きく低下する。一方、テクスチャメモリはキャッシュ機能が有るため、ランダムアクセスによるメモリアクセス速度の低下を抑えることができる。そこで、粒子情報を格納しているグローバルメモリ領域をテクスチャメモリから参照できるように設定した。さらに粒子には所属するセル番号順になるように番号を付けているために接触候補ペア同士の粒子番号が近く、加えて接触候補ペアリストは粒子番号ｉが小さい順にリスト内に並ぶように作られている。そのため、接触候補ペアリスト番号でスレッド化することでキャッシュヒット率が高くなり、粒子間の接触判定および接触力計算のルーチンが高速化される。
【００７７】
接触判定部１８は、作成／更新した接触力ボックスの情報を、積算接触候補数s_jgi[i]及び粒子情報と共に接触力計算部１９に送信する。
【００７８】
接触力計算部１９は、接触力参照表５４と積算接触候補数s_jgiを用いて各粒子ｉに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]を計算する。接触力参照表５４を接触力参照表作成部１７から受信し、積算接触候補数を接触判定部１８から受信すると、まず、着目粒子ｉごとに、接触力参照表５４及び積算接触候補数s_jgiを用いて、この着目粒子と接触している粒子との接触力を、接触力ボックスから抽出する。具体的には、着目粒子ｉより粒子径が大きい粒子については、接触力参照表５４に保持されるボックス番号（Lnum）（参照情報）に基づき、このボックス番号に対応する接触力Fij[Lnum]を接触力ボックスから抽出する。着目粒子ｉより粒子径が小さい粒子については、積算接触候補数s_jgi[]を用いて、接触力テーブルの格納位置を特定して接触力を抽出する。より詳細には、s_jgi[i-1]〜s_jgi[i]−1番のボックス番号（接触候補番号）に対応する接触力Fij[s_jgi[i-1]]〜Fij[s_jgi[i]−1]を接触力ボックスから抽出する。ここで接触力ボックスに格納されている接触力Fij[]とは、具体的には並進成分Fijt[Lnum]、回転成分Fijr[Lnum]、及び後述する（１３）式で用いるβ[Lnum]の３つの１次元配列を一組として保持するものである。
【００７９】
この接触力を接触力ボックスから抽出する処理を図１０を参照して具体的に説明する。図１０は、これまでの図６，９と同様に、着目粒子ｉ＝４の例である。接触力参照表５４には、着目粒子（ｉ＝４）より粒子径が大きい粒子（ｉ＝７，８）との接触力の参照先情報が含まれている。ｉ＝４の着目粒子とｉ＝７の粒子との接触力を格納する接触力ボックスのボックス番号（＝接触候補リスト番号）Lnumは４であり、ｉ＝８の粒子との接触力を格納する接触力ボックスのボックス番号Lnumは７である。接触力計算部１９は、接触力参照表５４に含まれるこのボックス番号を用いて、接触力ボックスから、ｉ＝７の粒子とのとの接触力Fij[４]，ｉ＝８の粒子との接触力Fij[７]を抽出する。
【００８０】
また、図１０に示すように、着目粒子ｉ＝４の積算接触候補数s_jgi[i]は、s_jgi[4]＝３であり、着目粒子の粒子番号より１つ前の粒子（図１０ではi＝３）の積算接触候補数s_jgi[i-1]は、s_jgi[3]＝３である。したがって、着目粒子ｉ＝４より粒子径が小さい粒子（この例ではi＝５，６）については、s_jgi[i-1]〜s_jgi[i]−1番、すなわち１番及び２番のボックス番号に接触力Fij[１]，Fij[２]を接触力ボックスから抽出する。
【００８１】
次に、着目粒子iごとに、接触力ボックスから抽出した接触力Fij[]（＝Fijt[]、Fijr[]、β[]）を用いて、粒子iに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]を計算する。
【００８２】
ｉ番の粒子に働く全接触力（並進成分）Ftall[i]は、接触候補リスト番号がs_jgi[i-1]〜s_jgi[i]−1番の接触力はそのままの符号で足し合わせる。さらに、接触力参照表５４（Ref）を用いて参照表番号0〜n_jli−1番の表から参照される接触力は逆符号にして足し合わせる。そして次式のように表される。
【数７】
【００８３】
ｉ番の粒子に働く全接触力（回転成分）Frall[i]も同様に足し合わせる。接触候補リスト番号がs_jgi[i-1]〜s_jgi[i]−1番の接触力はそのまま参照し、接触力参照表５４（Ref）を用いて参照表番号0〜n_jli−1番の表から参照される接触力は逆符号にしてからβを掛けてから足し合わせ、次式のように表される。ここで、βはｊ粒子の中心から接触点までの距離をｉ粒子の中心から接触点までの距離で割った値である。
【数８】
【００８４】
ここで、s_jgi[i-1]〜s_jgi[i]−1 番の接触力は粒子ｉが粒子ｊから受ける接触力であるのに対して、0〜n_jli−1番の接触力参照表５４で参照される接触力は粒子ｉが粒子ｊに与える接触力である。そのため， 0〜n_jli−1番の接触力参照表５４に対して、接触力の並進成分については符号を反対にしてから足し合わせ、接触力の回転成分についてはβ[box id]を掛けてから足し合わせる。この様に、粒子ｉが粒子ｊから受ける力を接触力参照表５４Ref[i][box id]を使わずにs_jgi[i]を用いて直接、接触力の入った配列Fijt[]，Fijr[]を参照することで、接触力参照表５４のメモリ使用量を増加させずに、大粒子に対して膨大な数の接触候補ペアとなる小粒子から受ける接触力を足し合わせることが可能になった。
【００８５】
コンピュータを用いた演算においては、具体的には粒子番号ｉをスレッド化し、box_idでループを回してFtall，Frallを求める。ここでFijへのメモリアクセスはランダムアクセスになるのでテクスチャメモリを使用する。上述のようにテクスチャメモリはキャッシュヒット率が高いので、グローバルメモリを直接アクセスするよりは高速に処理できる。なお接触力参照表５４に関しては配列Ref[A][B]のBの要素サイズを16の倍数になるように配列宣言することで連続アクセスが可能となる。
【００８６】
接触力計算部１９は、上述のように算出した、各粒子ｉに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]を、粒子情報と共に粒子情報更新部２０に送信する。
【００８７】
粒子情報更新部２０は、接触力計算部１９から各粒子ｉに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]と、粒子情報とを受信すると、各粒子の位置、速度、接触力計算部１９から受信した接触力に基づき、粒子の座標（位置情報）と、並進速度及び回転速度（速度情報）を算出し、これらの算出した座標、並進速度及び回転速度を用いて、粒子情報保持部１１に保持されている粒子情報を更新する。具体的には、並進の運動方程式（上記（４）式）をleap-flog法を用いて解くことで新しい時刻の並進速度ならびに粒子座標を求める。さらに、回転の運動方程式（上記（５）式）も同様にleap-flog 法を用いて解き、ここで得られる回転速度を予測値とする。この予測される回転速度に対して転がり摩擦力を計算し、接触力による回転力と転がり摩擦力を考慮した回転の運動方程式を解くことにより新しい時刻の回転速度を得る。
【００８８】
粒子情報更新部２０は、算出した各粒子の位置座標、並進速度、回転速度を用いて、粒子情報の座標、並進速度、回転速度を更新し、更新した粒子情報を粒子情報保持部１１に送信して粒子情報保持部１１の粒子情報を最新のものに更新する。
【００８９】
次に、図１１を参照して、粒子シミュレーション装置１０において実行される処理について説明すると共に、本実施形態に係る粒子シミュレーション方法について説明する。図１１は本実施形態の粒子シミュレーション装置１０において実行される処理を示すフローチャートである。粒子情報保持部１１に入力され保持されている粒子情報が、粒子情報取得部１２により粒子情報保持部１１から取得される（Ｓ１０１：入力ステップ）。接触候補リスト更新判定部１３により、接触候補リストの更新が必要か否かが判定される（Ｓ１０２）。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。更新が必要と判定された場合にはステップＳ１０３に移行し、更新が不要と判定された場合にはステップＳ１０８に移行する。
【００９０】
次に、粒子番号変更部１４により粒子番号が変更される（Ｓ１０３：付与ステップ）。粒子番号変更部１４は、粒子を所属するセル番号順に並べ、その順に粒子番号を付け直し、さらに、同じセル番号の粒子について粒子径順に並べ直し、その順に粒子番号を付け直す。続いて、近傍粒子表作成部１５により、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２が作成される（Ｓ１０４）。近傍粒子表５２に近傍粒子として格納する条件は、上述の条件（ａ）〜（ｃ）である。
【００９１】
次に、接触候補リスト作成部１６により、近傍粒子表５２の条件（ａ）を満たすが（ｂ）及び（ｃ）を満たしていない粒子数n_jgi[i] のプレフィックス和である積算接触候補数s_jgi[i]が算出され（Ｓ１０５：積算接触候補数算出ステップ）、次いで、作業空間５１の全粒子について、接触している可能性がある粒子ペアの情報を含む接触候補リスト５３が作成される（Ｓ１０６：接触候補選択ステップ）。接触候補リスト５３に格納する条件は、上述の条件ａ）〜ｃ）である。接触力参照表作成部１７により、接触力参照表５４が作成される（Ｓ１０７：接触力参照表作成ステップ）。
【００９２】
次に、接触判定部１８により、接触候補リスト５３を参照して、２つの粒子間の接触判定が行われ（Ｓ１０８）、接触していると判定された場合には粒子間の接触力Fijr、Fijtが（１），（２）式を用いて算出され、接触力ボックスの配列に格納される（Ｓ１０９：接触力算出ステップ）。また、粒子が接触していないと判定された場合には接触力Fijr、Fijtを０として、接触力ボックスの配列に格納される。
【００９３】
接触力計算部１９により、接触力参照表５４及び積算接触候補数s_jgi[i]を参照して、着目粒子と接触する他の粒子が選択され、接触力ボックスから対応する配列Fijr、Fijtから接触力が取得される。そして（１２），（１３）式を用いて接触力が総和演算されて、各粒子に働く接触力Fr、Ftが計算される（Ｓ１１０：接触力総和演算ステップ）。粒子情報更新部２０により、ステップＳ１０９において算出された接触力、現在の粒子情報に基づいて、各粒子の座標（位置情報）と、並進速度及び回転速度（速度情報）とが算出され、これらの算出した座標、並進速度及び回転速度を用いて、粒子情報保持部１１に保持されている粒子情報が更新される（Ｓ１１１：粒子情報更新ステップ）。
【００９４】
そして、粒子シミュレーションが終了したか否かが判定される（Ｓ１１２）。終了判定は、例えば所定時間が経過した場合にシミュレーションが終了したものと判定される。粒子シミュレーションが終了していないと判定された場合にはステップＳ１０１に戻り、次のステップに入る。粒子シミュレーションが終了したと判定された場合には処理を終了する。
【００９５】
次に、図１２を参照して、本実施形態の粒子シミュレーション装置１０によるメモリ使用量について説明する。図１２には、粒子モデルの最大粒子径と最小粒子径との比を変化させた場合の本実施形態（実線）及び比較例（破線）のメモリ使用量を示している。なお、粒子モデルの粒子数は１，０００，０００個とし、最大粒子径の粒子（最大粒子）と最小粒子径の粒子（最小粒子）とが混在する二成分系とした。比較例は、近傍粒子表５２及び接触力参照表５４に着目粒子ｉの粒子径より小さい粒子の情報も格納する点が本実施形態と異なるものである。図１２の横軸は、「最大粒子径／最小粒子径」を表し、縦軸は「メモリ資料量（MB）」を表している。
【００９６】
比較例の場合、近傍粒子表５２nei_pn[i][box_id]及び接触力参照表５４Ref[i][box_id]に必要なメモリ使用量は、１データに４バイト必要として、粒子数×最大接触候補粒子数×４byteとなる。ここで、最大接触候補粒子数は、近傍粒子表５２及び接触力参照表５４に格納する近傍粒子の最大個数であって、最大粒子の表面を最小粒子が埋め尽くす状態を考えればよく、この場合、
最大接触候補粒子数＝４π(ｒmax+ｒmin)２／πｒmin２ （１４）
となる。ここで、ｒｍａｘは最大粒子の粒子径であり、ｒｍｉｎは最小粒子の粒子径である。つまり、比較例の場合、図１２に示すように、最小粒子径に対して最大粒子径が大きくなるほど、つまり最大粒子径／最小粒子径が大きくなるほど、最大接触候補粒子数が増大し、この結果、メモリ使用量が増大する。例えば、粒子径比１：１０の条件（最大粒子径／最小粒子径＝１０）では、一粒子当たりの接触候補粒子数は最大で５００個程度に増加するため、必要なメモリ量で２ＧＢにもなる。
【００９７】
一方、本実施形態の場合、着目粒子より粒子径の小さい粒子は近傍粒子として考慮されず、近傍粒子表５２及び接触力参照表５４に格納されないので、最大接触候補粒子数を算出する上記（１４）式は成り立たず、本実施形態の最大接触候補粒子数は、最大と最小の粒子径に依存しない。本実施形態では、最大接触候補粒子数は、上述した近傍粒子の判定条件（ａ）において粒子間距離の閾値であるdis_lim = (rｉ + rｊ)(1.0+α)に含まれるパラメータαのみに依存し、粒子径分布には依存せず一定となる。α＝０．２の場合は、最大接触候補粒子数は如何なる粒子径分布条件においても１６も有れば十分であるので、ここでは最大接触候補粒子数を１６としてメモリ使用量を計算した。その結果、図１２に示すように、本実施形態のメモリ使用量は、粒子数×最大接触候補粒子数×４byte ＝1000000×16×４＝64MBの一定値となった。
【００９８】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、作業空間５１内の粒子のそれぞれについて、当該粒子ｉの近傍に位置する、この粒子ｉと接触する可能性の高い他の粒子のうち当該粒子ｉの粒子径より小さい粒子については、粒子の数を示す積算接触候補数s_jgi[i]のみを保持しておく。そして、接触判定部１８が、この積算接触候補数s_jgi[i]を利用して接触力テーブルの格納位置を特定して当該粒子と他の粒子との接触力を抽出する。これにより、接触力参照表５４には、当該粒子ｉの近傍に位置する他の粒子のうち当該粒子ｉの粒子径より小さい粒子に関する情報を含める必要がなく、当該粒子の粒子径より大きい粒子に関する情報のみを保持すればよくなる。粒子径分布を有する粒子群をシミュレーション対象とする場合、自分より粒子径の小さい粒子まで接触を考慮すると、同時に接触可能な粒子の数が増大するため、予め接触力参照表５４のために確保すべきメモリ領域を多くしなければならない。これに対し、上記構成により、自分より粒子径の大きい粒子との接触のみを考慮すればよく、接触を考慮する粒子数を減らすことができ、接触力参照表５４を作成するために予め割り当てるメモリ領域を少なくすることが可能となり、この結果、メモリ使用量を抑制することができる。このようにメモリ使用量が抑制されることにより、ＧＰＵなどの共有メモリ型並列計算機を利用して粒子シミュレーションを実行することが可能となる。
【００９９】
また、粒子番号変更部１４が、作業空間５１内の全粒子について、セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した粒子番号を付与し、各粒子をセル番号順、さらには同一セル内では粒子径順にソーティングする。このため、接触力参照表を作成する際に、粒子番号に基づき粒子径が小さい粒子を容易に認識することができ、粒子径が小さい粒子については接触力参照表にいれるか否かの判定を行う必要がなくなる。この結果、接触力参照表の作成において計算効率を向上させることができる。
【０１００】
以上、本発明について好適な実施形態を挙げて説明したが、本発明は上記実施形態に限定されるものではない。例えば、上記実施形態では、粒子番号変更部１４は、（１）粒子を所属するセル番号順に並べ、さらに、（２）同じセル番号の粒子について粒子径順に並べ直す、というソーティングを行っているが、（１）粒子を所属するセル番号順に並べる処理のみを行うよう構成してもよい。
【０１０１】
また、上記実施形態では、本発明の適用対象としてＧＰＵを挙げたが、ＧＰＵと同様の共有メモリ型のＳＩＭＤ（Single Instruction Multiple Data）型並列処理を行うハードウェア（例えばベクトルプロセッサ）、すなわち共有メモリ型並列計算機にも適用可能である。
【符号の説明】
【０１０２】
１０…粒子シミュレーション装置、１１…粒子情報保持部（入力手段、粒子情報保持手段）、１２…粒子情報取得部、１３…接触候補リスト更新判定部、１４…粒子番号変更部（付与手段）、１５…近傍粒子表作成部、１６…接触候補リスト作成部（接触候補選択手段、積算接触候補数算出手段）、１７…接触力参照表作成部（接触力参照表作成手段）、１８…接触判定部（接触力算出手段）、１９…接触力計算部（接触力総和演算手段）、２０…粒子情報更新部（算出手段、粒子情報更新手段）、５１…作業領域、５２…近傍粒子表、５３…接触候補リスト、５４…接触力参照表。
【技術分野】
【０００１】
本発明は、粒子シミュレーション装置及び粒子シミュレーション方法に関する。
【背景技術】
【０００２】
土木工学の研究分野から端を発したＤＥＭ（Discrete Element Method：離散要素法）は、その後粉体工学へ応用され、昨今では、粒子性物質のみならず連続体として扱いが困難な非常に多種多様の問題に対する数値解析手法として益々需要は広がってきた。一方、離散的な粒子運動を数値的に追跡する手法には、ＭＤ（Molecular Dynamics）、ＬＧＡ（Lattice Gas Automaton）、ＳＰＨ（Smoothed Particle Hydrodynamics）などＤＥＭ以外にも数多く開発されてきている。しかし、数ある数値粒子解析手法の中で、ＤＥＭの応用範囲が広い理由は、粒子の大きさや形状という個々の粒子の幾何情報を考慮しながら粒子間相互作用力を求めニュートンの運動方程式を解くという計算を非常に高速に行えることである。勿論、ＤＥＭ以外にもＤＤＡ（Discontinuous Deformation Analysis）なども、大きさと形状を考慮することが可能な手法であるが、もともと計算負荷が高いので、実用的にあまり多数の粒子の扱うことができない。
【０００３】
ところで、砂や小麦粉のように我々が目にする実際の粒状体や粉体（以下、両者をまとめて「粉粒体」と呼ぶ）は、篩い分けによる粒度調整など特別な処理を施さない限り、一般的には様々な大きさの粒子が混在した状態で存在している（以下、これを「粒子径分布」がある状態と呼ぶ）。実はこの粒子径分布は、粉粒体の均質性や変形や流動といった力学的な性質を支配する大きな要因であり、粉粒体を扱う産業では粒径分離効率や粒度偏析による粉体層の不均質化および混合効率、輸送・充填の制御にはいつも問題となる。したがって、粒子径分布が粒子群全体の振る舞いに与える影響を正しく理解することは、粉粒体を扱うためには必要不可欠である。
【０００４】
ＤＥＭは、前述したように粒子の大きさを考慮できる手法であるから、粒子径分布がもたらすマイクロメカニックスをＤＥＭで調べたくなる気持ちはごく自然である。しかし、粒子径が１／ｎになると粒子１個の体積が(１／ｎ)３になり、もとのサイズの粒子１個分の体積を充填するためにｎ３個の粒子が必要となる（厳密には空隙が生じるため、ｎ３個よりは少なくなるがオーダーは変わらない）ため、実は、計算機のメモリ量や計算時間の制約上、あまり大きな粒子径分布は扱えない。例えば自然土には数センチの砂利から数ミクロンの粘土粒子まで粒子径は約５桁程度に渡って分布している。すると、全てのサイズで同程度の体積を占めている場合、最大粒径の粒子数を１個にしても、全体では(１０５)３個オーダーの粒子数となってしまう。
【０００５】
また、ＤＥＭで粒子径分布を扱う場合には、上述の粒子数の問題に加えて、粒子径の違いによる接触点数のアンバランスも問題となる。例えば、最大粒子と最小粒子の粒径比をｒ：１とした時、最大粒子の表面を最小粒子が埋め尽くす場合と、最小粒子の表面を最大粒子が埋め尽くす場合を考えると、表面を覆う粒子の個数はそれぞれ４(ｒ＋１)２／１２および４(ｒ＋１)２／ｒ２となる。すなわち、接触点数のアンバランスは最大でｒ２：１のオーダーになる。個々の粒子の接触点数のアンバランスは、粒子ごとの接触力の計算回数と接触の総和計算量のアンバランスにつながり、並列計算には大きな障壁となる。つまり、ＤＥＭでは、均一粒径の場合と違って粒径比が大きくなると、単純に粒径の幾何学的な情報だけからでも粒子数が増加し接触点数のアンバランスが生じることが見込まれ、結果としてかなり計算負荷が高くなることが分かる。このことに加えて、実際のＤＥＭのプログラム内部の処理においても、粒子径分布が存在すると、接触判定の対象となる接触候補の絞り込み検索と接触候補ペアを格納するための処理が、均一粒径の場合と比較して計算負荷が著しく高くなってしまう。
【０００６】
この粒子径分布がある問題に対して、接触候補ペアの抽出と格納に関連する処理を高速化するためのＤＥＭプログラミングの手法は幾つか提案されてきている（例えば特許文献１）。
【０００７】
一方、最近ではＣＰＵと比べて価格に対する演算性能の比が格段に大きいＧＰＵ（Graphics Processing Unit）を用いた数値シミュレーションの高速化が盛んに行われている。これに伴い、粒子系シミュレーションについても、高価なスーパーコンピュータを用いたベクトル化処理よりも、ベクトル演算と同等の超並列演算が可能で安価なＧＰＵに対する高速化アルゴリズムが提案されてきている（例えば特許文献２）。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２００７−２８６５１４号公報
【特許文献２】特開２００９−６９９３０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかし、特許文献１など従来のＤＥＭ高速化手法は、まずセルに粒子を格納してからセルサイズやセル探索範囲を最適化する手法であって、並列化やベクトル化による大規模高速化をことさら意識したものではない。具体的な例を挙げると、粒子番号のループで粒子をセルに格納する処理をベクトル化のようなＳＩＭＤ（Single Instruction Multiple Data）型並列処理で行うと、２つ以上の粒子が同一セルに格納されるときにはメモリ書き込み競合が生じるために非常に都合が悪い。このメモリ書き込み競合問題は、単一粒径成分の場合でも一つのセルに複数の粒子が入ってしまう場合に発生するが、粒子径分布がある場合は一つのセルに入る粒子数自体に大きな分布ができるので、並列化やベクトル化によるプログラムの高速化の際には、さらに厄介な大きな問題となる。
【００１０】
また、特許文献２などのＧＰＵ高速化アルゴリズムでは、接触候補粒子ペアの検出処理においては、ＣＰＵベースのプログラムと同様にセルに粒子を格納してからセルを探索する方法が用いられている。しかし上述したように、ＧＰＵでも複数粒子が同一セルに格納される場合には、やはりメモリ書き込み競合が生じるため、この問題を回避するために、幾つかの対策が提案されている。一つは、粒子の格納と確認操作を複数回繰り返す方法で、他方は、コンパイラ付属の排他処理関数（Atomic関数）を用いる方法である。また、粒子間接触力の計算ルーチンにおいてもメモリ書き込み競合を防ぐために、粒子ｉ‐ｊ間の接触力計算を粒子ｉと粒子ｊそれぞれに対して同じ計算を二度行うことで問題を回避する工夫がなされている。
【００１１】
しかし、前述のように、粒子径分布があると、セルによって格納される粒子数に大きな偏りが生じやすく、また大粒子に対する接触粒子数は小粒子に対するそれと比べて非常に多いために、ロードバランスを大きく損ないやすい。例えば、これらの方法を用いた三次元ＤＥＭで二成分系粒子群に対するＧＰＵ計算を行うと、二成分の粒子径比が１：４になるだけで単一成分の場合と比べて計算時間は３倍近くにもなることが報告されている。
【００１２】
そして、粒子径分布を考慮した粉粒体のＤＥＭシミュレーションを行う場合、粒子径を均一とした同じ粒子数の場合に比べて、接触力を考慮すべき近傍粒子の数が急激に増大するので、ＤＥＭ計算を実現するために確保しなければならないメモリ使用量が増大するという問題があった。
【００１３】
ＧＰＵ演算は、ＣＰＵ演算のようにメモリをハードウェア的に増設することができないため、あくまでグラフィックカードに搭載された限られたメモリ容量内でしか処理できない。勿論、ＧＰＵメモリの不足分を転送通信によってＣＰＵメモリに頼る方法も考えられるが、ＧＰＵメモリとＣＰＵメモリ間の転送速度は非常に遅いため、せっかくのＧＰＵの高い計算効率が全く生かされず推奨できない。このように、ＧＰＵをＤＥＭ計算に使用する場合、現状のＧＰＵのメモリ容量が最大でも４ＧＢ程度であるので、このメモリ使用量の問題が、ＧＰＵを粒子シミュレーションに適用することの制約となっている。そして、この問題は、ＧＰＵだけでなく、ＧＰＵと同様の共有メモリ型のＳＩＭＤ（Single Instruction Multiple Data）型並列処理を行うハードウェア（例えばベクトルプロセッサ）（以下、「共有メモリ型並列計算機」という）にも該当する問題である。
【００１４】
本発明は、上記課題を解決するためになされたものであり、粒子径分布を有する粉粒体のＤＥＭ（離散要素法）計算を共有メモリ型並列計算機により実行する際に、メモリ使用量を抑制することができる粒子シミュレーション装置及び粒子シミュレーション方法を提供することを目的とする。
【課題を解決するための手段】
【００１５】
上記課題を解決するため、本発明に係る粒子シミュレーション装置は、作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力する入力手段と、粒子情報を保持する粒子情報保持手段と、粒子のそれぞれについて、位置情報に応じた順序で特定する特定情報を付与する付与手段と、粒子の中から前記特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択手段と、選択候補選択手段により選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の位置情報及び速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出手段と、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を前記接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成手段と、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出手段と、粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を接触力参照表の参照情報を用いて接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、積算接触候補数を用いて接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算手段と、接触力総和演算手段により計算された粒子ごとの接触力に基づいて、粒子の新たな位置情報及び速度情報を算出する算出手段と、算出手段により算出された新たな位置情報及び速度情報を用いて、粒子情報保持手段に保持される粒子情報を更新する粒子情報更新手段と、を備えることを特徴とする。
【００１６】
同様に、上記課題を解決するため、本発明に係る粒子シミュレーション方法は、作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置における粒子シミュレーション方法であって、粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力し、この粒子情報を粒子情報保持手段に保持する入力ステップと、粒子のそれぞれについて、位置情報に応じた順序で特定する特定情報を付与する付与ステップと、粒子の中から特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択ステップと、接触候補選択ステップにおいて選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の位置情報及び速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出ステップと、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成ステップと、粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出ステップと、粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を接触力参照表の参照情報を用いて接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、積算接触候補数を用いて接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算ステップと、接触力総和演算ステップにおいて計算された粒子ごとの接触力に基づいて、粒子の新たな位置情報及び速度情報を算出する算出ステップと、算出ステップにおいて算出された新たな位置情報及び速度情報を用いて、粒子情報保持手段に保持される粒子情報を更新する粒子情報更新ステップと、を含むことを特徴とする。
【００１７】
このような粒子シミュレーション装置及び粒子シミュレーション方法によれば、作業空間内の粒子のそれぞれについて、当該粒子の近傍に位置する、この粒子と接触する可能性の高い他の粒子のうち当該粒子の粒子径より小さい粒子については、粒子の数を示す積算接触候補数のみを保持しておき、この積算接触候補数を利用して接触力テーブルの格納位置を特定して当該粒子と他の粒子との接触力を抽出する。これにより、接触力参照表には、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子に関する情報を含める必要がなく、当該粒子の粒子径より大きい粒子に関する情報のみを保持すればよくなる。粒子径分布を有する粒子群をシミュレーション対象とする場合、自分より粒子径の小さい粒子まで接触を考慮すると、同時に接触可能な粒子の数が増大するため、予め接触力参照表のために確保すべきメモリ領域を多くしなければならない。これに対し、本発明の上記構成により、自分より粒子径の大きい粒子との接触のみを考慮すればよく、接触を考慮する粒子数を減らすことができ、接触力参照表を作成するために予め割り当てるメモリ領域を少なくすることが可能となり、この結果、メモリ使用量を抑制することができる。このようにメモリ使用量が抑制されることにより、ＧＰＵなどの共有メモリ型並列計算機を利用して粒子シミュレーションを実行することが可能となる。
【００１８】
また、粒子シミュレーション装置では、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、粒子情報には、位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、付与手段が、粒子のそれぞれについて、セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することが好適である。
【００１９】
同様に、粒子シミュレーション方法では、作業空間が複数のセルに分割され、各セルにはセル番号が付されており、粒子情報には、位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、付与ステップが、前記粒子のそれぞれについて、セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することが好適である。
【００２０】
この構成により、各粒子にセル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報が付与されるので、各粒子がセル番号順、さらには同一セル内では粒子径順にソーティングされる。このため、接触力参照表を作成する際に、特定情報に基づき粒子径が小さい粒子を容易に認識することができ、粒子径が小さい粒子については接触力参照表にいれるか否かの判定を行う必要がなくなる。この結果、接触力参照表の作成において計算効率を向上させることができる。
【００２１】
また、粒子シミュレーション装置では、接触候補選択手段が、粒子郡の中から特定情報の昇順で着目粒子を選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との組と、この組に昇順で付される接触候補番号とからなる接触候補ペア情報を、作業空間内の全粒子について含む接触候補リストを作成し、接触力算出手段が、接触候補リストに含まれる接触候補ペア情報のそれぞれについて、該接触候補ペア情報に関する２つの粒子の間の接触力を、当該粒子の前記位置情報及び速度情報に基づき算出し、算出した接触力と接触候補番号とを関連付けた接触力情報を接触力テーブルに格納し、接触力参照表作成手段が、粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の大きい他の粒子との接触力を前記接触力テーブルから参照するための参照情報として、該他の粒子及び当該粒子に関する接触候補番号を、該他の粒子及び当該粒子と関連付けて保持する接触力参照表を作成し、積算接触候補数算出手段が、粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の小さい他の粒子の数を示す接触候補数を前記特定情報順で積算した値である積算接触候補数s_jgi[i]（ｉは当該粒子の特定情報を示す）を算出し、接触力総和演算手段が、粒子のそれぞれについて、接触力参照表の参照情報を用いて接触力を接触力テーブルから抽出し、また、積算接触候補数s_jgi[i]を用いて接触候補番号がs_jgi[i-1]〜s_jgi[i]‐1の接触力を接触力テーブルから抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算することが好適である。
【発明の効果】
【００２２】
本発明に係る粒子シミュレーション装置及び粒子シミュレーション方法によれば、粒子径分布を有する粉粒体のＤＥＭ（離散要素法）計算を共有メモリ型並列計算機により実行する際に、メモリ使用量を抑制することができる。
【図面の簡単な説明】
【００２３】
【図１】本発明の一実施形態に係る粒子シミュレーション装置の機能ブロック図である。
【図２】本実施形態で粒子モデルに適用するＶｏｉｇｔモデルを示す概略図である。
【図３】本実施形態における作業領域を示す図である。
【図４】粒子番号変更の流れを示す図である。
【図５】作業領域における粒子番号を変更する処理と、各セルに粒子番号の最大値及び最小値を記憶する処理を示す図である。
【図６】近傍粒子表の構成の一例を示す図である。
【図７】接触候補リストの構成の一例を示す図である。
【図８】現在の接触候補リストのペア粒子について１ステップ前にペアであったか調べる処理を示す図である。
【図９】接触力参照表の構成の一例を示す図である。
【図１０】接触力参照表及び積算接触候補数と、接触力ボックスとの対応関係を示す図である。
【図１１】本実施形態の粒子シミュレーション装置において実行される処理を示すフローチャートである。
【図１２】粒子モデルの最大粒子径と最小粒子径との比を変化させた場合の本実施形態及び比較例のメモリ使用量を示し図である。
【発明を実施するための形態】
【００２４】
以下、図面を参照しながら本発明の実施形態を詳細に説明する。なお、図面の説明において同一又は同等の要素には同一の符号を付し、重複する説明を省略する。
【００２５】
図１は、本発明の一実施形態に係る粒子シミュレーション装置１０の機能ブロック図である。図１に示すように、粒子シミュレーション装置１０は、粒子情報保持部（入力手段、粒子情報保持手段）１１、粒子情報取得部１２、接触候補リスト更新判定部１３、粒子番号変更部（付与手段）１４、近傍粒子表作成部１５、接触候補リスト作成部（接触候補粒子選択手段、積算接触候補数算出手段）１６、接触力参照表作成部（接触力参照表作成手段）１７、接触判定部（接触力算出手段）１８、接触力計算部（接触力総和演算手段）１９、粒子情報更新部（算出手段、粒子情報更新手段）２０を備えている。
【００２６】
粒子シミュレーション装置１０は、物理的には、ＧＰＵ（Graphics Processing Unit）、主記憶装置であるＲＡＭ（Random AccessMemory）及びＲＯＭ（ReadOnly Memory）、ハードディスク装置等の補助記憶装置、入力デバイスである入力キー等の入力装置、ディスプレイ等の出力装置、ネットワークカード等のデータ送受信デバイスである通信モジュールなどを有するコンピュータとして構成されている。図１を参照して以下に説明する粒子シミュレーション装置１０の各機能は、ＧＰＵ、ＲＡＭ等のハードウェア上に所定のコンピュータソフトウェアを読み込ませることにより、ＧＰＵの制御のもとで入力装置、出力装置、通信モジュールを動作させるとともに、ＲＡＭや補助記憶装置におけるデータの読み出し及び書き込みを行うことで実現される。
【００２７】
ここで、まず本実施形態で適用する粒子モデルについて説明する。本実施形態の粒子シミュレーション装置１０が対象とするのは、粒子径分布を有する粉粒体群のＤＥＭシミュレーションである。つまり、様々な大きさの粒子径をもつ粒子が混在した状態の粒子モデルが適用される。
【００２８】
そして、このような粒子モデルにおいて、各粒子間に発生する粒子間接触力の法線方向および接線方向成分Ｆｎ，Ｆｔは、図２に示すＶｏｉｇｔモデル７１を用いて以下の式により表される。
【数１】
ここでｘとｖは相対変位および相対速度ベクトル、Ｋは線形バネ７２における弾性係数、ηはダッシュポット７３の粘性減衰係数であり、それぞれの添字ｎ，ｔにより法線方向成分及び接線方向成分を表す。μtは滑り摩擦係数、ｍｉｎは絶対値が小さい方の値を表す関数である。粘性減衰係数ηは、反発係数eと次式の関連がある。
【数２】
ここでｍｉ，ｍｊは、それぞれ粒子ｉ，ｊの質量を表す。添字ｑは、法線方向成分の場合はｑ＝ｎとなり、接線方向成分の場合にはｑ＝ｔとなる。
【００２９】
したがって、ＤＥＭで扱う並進および重心を中心とする回転の運動方程式は以下の式で表される。
【数３】
ここで、ｕｐとωｐは粒子の速度および角速度、ｍｐとＩｐは粒子の質量および慣性モーメント、Ｒは粒子中心から接触点へ向かう位置ベクトルであり、添字ｐは、粒子ｉの場合はｐ＝iとなり、粒子ｊの場合はｐ＝ｊとなる。また、（５）式の右辺第二項は転がり摩擦力を表し、μｒは転がり摩擦係数、ｂは接触面幅である。（１）、（２）式のＶｏｉｇｔモデル７１は、後述する接触判定部１８で粒子間接触力を算出するのに用いられ、また（４）、（５）式の運動方程式は、粒子情報更新部２０で粒子の座標を算出するのに用いられる。
【００３０】
また、本実施形態において粒子が運動する領域である作業領域５１は、図３に示すように、一辺の大きさがcsizeの立方体セルで区分される。セルサイズcsizeは、パラメータφ (0≦φ≦ 1) を用いて以下の式で与える。
csize＝2.0{(1.0−φ)rmax＋φrmin)}(1.0＋α) （６）
ここで、rmaxは作業領域に配置される粒子モデルの最大粒子径であり、rminは最小粒子径であり、αは更新頻度を調整するパラメータであり例えばα＝０．２である。また、作業領域５１内の各セルにはセル番号cell_idが付されている。セル番号は、図３に示すように一次元であり、例えば次式のように表される。
cell_id＝i.x＋id.x×i.y＋id.x×id.y×i.z （７）
i.m＝{0, 1, 2,…, id.m-1} （m＝x, y, z）
ここでcell_idはセル番号であり、i.mはｍ方向のセル位置であり、id.mはｍ方向のセル数である。つまり、粒子シミュレーション装置１０は、シミュレーション条件として粒子の最大粒子径rmax、最小粒子径rmin、パラメータφ，αが設定されると、作業領域５１のセルサイズcsizeを自動的に算出して設定し、さらに画定したセルにセル番号を割り当てることができる。
【００３１】
粒子情報保持部１１は、作業空間５１内の複数の粒子のそれぞれについて粒子情報を保持している。粒子情報は、座標（位置情報）、並進速度（速度情報）、回転速度（速度情報）、粒子半径、質量、粒子番号、セル番号を含む。座標とは、図３に示す三次元の作業空間５１における三次元座標（pos.x，pos.y，pos.z）である。粒子番号は各粒子を識別するための番号である。セル番号は当該粒子が所属するセルのセル番号である。なお、シミュレーションを開始する際には、粒子情報保持部１１は、シミュレーションに適用する粒子径分布に応じて各粒子の粒子半径、質量を設定し、座標、速度、粒子番号を各粒子にランダムに割り当てる。セル番号には、粒子ごとにランダムに割り当てられた座標に基づき特定される、この粒子が属するセルの番号を付与する。そして、このように生成した各情報を初期状態の粒子情報として入力し保持する。また、シミュレーションの実行中には、後述する粒子情報更新部２０により粒子情報の各情報が更新されると、粒子情報保持部１１は、保持している粒子情報を、粒子情報更新部２０により更新された新しいものに置き換える。
【００３２】
粒子情報取得部１２は、粒子情報保持部１１から粒子情報を取得する。粒子情報取得部１２は、粒子シミュレーションの開始時に粒子情報保持部１１から予め設定されている初期状態の粒子情報を取得し、シミュレーション動作中には、粒子情報保持部１１が更新される度に新たな粒子情報を粒子情報保持部１１から取得する。粒子情報取得部１２は、取得した粒子情報を接触候補リスト更新判定部１３に送信する。
【００３３】
接触候補リスト更新判定部１３は、後述する接触候補リスト５３を更新するか否かを判定する。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。また、粒子の積算移動距離は、例えば、粒子情報取得部１２が新しい粒子情報を粒子情報保持部１１から取得する度に、１ステップ前の粒子情報との座標データの差分からこの直近の１ステップにおける粒子の移動距離を算出し、これを粒子ごとに積算して求める。接触候補リスト更新判定部１３は、１ステップ前の粒子情報と、パラメータαと、粒子の積算移動距離とを記憶しており、粒子情報取得部１２から新たな粒子情報を受信すると、積算移動距離を更新し、この更新された積算移動距離がｒαより大きい粒子の有無を判定する。そして、積算移動距離が所定値以上の粒子が有る場合に、接触候補リストを更新するものと判定し、新たな粒子情報を粒子番号変更部１４に送信する。なお、粒子移動距離の積算値は、接触候補リスト５３の更新が行われるたびに初期化される。また、更新された積算移動距離がｒαより大きい粒子が無い場合には、接触候補リストを更新しないものと判定し、粒子番号変更部１４には粒子情報を送信せず、接触判定部１８に新たな粒子情報を送信する。
【００３４】
粒子番号変更部１４は、作業領域５１内の各粒子のそれぞれについて、セル番号に応じた順序で特定し（すなわち位置情報に応じた順序で特定し）、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した粒子情報（特定情報）を付与する。言い換えると、粒子番号変更部１４は、作業領域５１内の各粒子が所属するセル及びその粒子径に基づき、各粒子の粒子番号を変更して、粒子番号をソーティングする。
【００３５】
粒子番号変更部１４は、接触候補リスト更新判定部１３から粒子情報を受信すると、図４に示すように、各粒子の粒子番号に、この粒子が所属するセルのセル番号と、この粒子の粒子径とを関連付けて1組の情報セットとし、これらの情報セットを粒子番号順に記憶する（図４（ａ））。次に、これらの情報セットをセル番号が昇順になるように並べ換え（図４（ｂ））、セル番号順に粒子番号を付け直す（図４（ｃ））。続いて、同一のセル番号の粒子群ごとに粒子径を比較し、これらの粒子群の中で粒子径の大きい順に並べ換え（図４（ｄ））、粒子径順に粒子番号を付け直す（図４（ｅ））。
【００３６】
この粒子番号を変更する処理により、例えば図５（ａ）に示すように作業空間５１に配置された粒子の粒子番号が、図５（ｂ）に示すように、セル番号及び粒子径に沿った粒子番号に変更される。粒子番号変更部１４は、粒子番号の変更が反映された粒子情報を近傍粒子表作成部１５に送信する。
【００３７】
近傍粒子表作成部１５は、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２を作成する。近傍粒子表作成部１５は、着目粒子の所属するセル、およびこのセルに隣接するセルや隣接セルのさらに外側に隣接する近傍セルに存在する粒子のうち、自身より粒子径が大きい粒子、或いは粒子径が同じ場合は初期配置時の粒子番号が自身より大きい粒子で、粒子間距離がある値以下の粒子を「近傍粒子」とする。以下に詳細を説明する。
【００３８】
近傍粒子表作成部１５は、近傍粒子表５２を作成するためのメモリ領域を予め確保している。確保するメモリ領域は、各粒子と接触する可能性がある近傍粒子の個数を考慮して、各粒子の近傍粒子表ごとに同一の大きさである。
【００３９】
近傍粒子表作成部１５は、粒子番号変更部１４から粒子情報を受信すると、同じセル番号をもつ粒子の中で、セル内の粒子の番号のうち、最小値Pnum_in_cellminと最大値Pnum_in_cellmaxをセルに記憶させる（図５（ｃ））。図５（ｃ）では、各セルの左上に最大粒子番号、右下に最小粒子番号が記憶されている。つまり、セルcell_idごとに、このセルに含まれる粒子の番号の最小値Pnum_in_cellmin[cell_id]及び最大値Pnum_in_cellmax[cell_id]の情報を格納する配列のみを保持しておく。これにより、各セル番号cell_idのセルには、Pnum_in_cellmin[cell_id]番からPnum_in_cellmax[cell_id]番までの粒子が粒子径の大きい順に所属していることがわかる。
【００４０】
番号cell_idのセルに隣接するセルの番号nei_cell_idは次式で表される。
nei_cell_id[n]= cell_id + con[n] （８）
ここで、conはid.xとid.yによって決まる定数である。セル番号で粒子番号をソートしてあるため、各隣接セルにはPnum_in_cellmin[nei_cell_id]からPnum_in_cellmax[nei_cell_id]までの番号の粒子が連続して存在することがわかる。
【００４１】
次に、近傍粒子表作成部１５は、各粒子ごとに、図６に示すような近傍粒子を格納する箱（配列）nei_pn[i][box_id]を用意し、粒子ｉの周囲の探索範囲（定義は後述する）に含まれる粒子のうち
（ａ）粒子ｉとの粒子間距離がdis_limよりも小さく、
（ｂ）粒子ｉの径よりも大きく、
（ｃ）粒子ｉの径と同じ場合は、初期配置時の粒子番号（i=p_id）が粒子ｉよりも大きい
という条件を満たす粒子を番号（box_id）が小さい方の箱から詰めて格納させる。ここで、dis_lim = (rｉ + rｊ)(1．0+α)であり、rｉ，rｊは近傍粒子ペアｉ，ｊの粒子半径である。αは接触候補の更新回数を調整するパラメータであり、例えばα＝０．２である。rｉ，rｊは粒子番号変更部１４から受信した粒子情報に含まれ、αは近傍粒子表作成部１５が予め保持している値であり、近傍粒子表作成部１５は、これらのデータを用いてdis_limを算出する。
【００４２】
上記条件（ａ）〜（ｃ）を満たす粒子を、粒子iに対する「近傍粒子」と呼ぶ。図６に示す一連の箱の集合を、作業空間５１内の全ての粒子ｉについて用意したものを、本実施形態では近傍粒子表５２という。近傍粒子表５２は、粒子i及びこの粒子iの近傍粒子の数box_idに関する二次元配列nei_pn[i][box_id]である。
【００４３】
この方法によって、着目粒子iについて、この粒子i以上の粒子径を持った粒子ｊのみが近傍粒子表５２に記憶されるため、近傍粒子表５２に必要なメモリ使用量の増加を抑制することができる。
【００４４】
ここで、近傍粒子表作成部１５は、粒子ｉの近傍粒子を決定するために、粒子ｉから(rmax＋ri)(1.0＋α)の距離にある粒子を探索する必要がある。粒子ｉが属するセルに隣接する２６個のセルのほかに、粒子モデルの最大粒子の半径ｒmaxに応じてさらに外側のセルまで探索範囲を広げることとなる。
【００４５】
近傍粒子表作成部１５が探索する近傍セルの範囲は、上記（６）式で算出されたセルサイズcsizeと、作業領域内の最大粒子の粒子径を示す最大粒子半径ｒmaxと、着目粒子半径ｒｉから次式によって求められる。
nei_cell＝2・ceil{2.0(rmax＋ri)(1.0＋α)／csize}−1 （９）
ceilは小数点以下を切り上げる関数であり、近傍粒子表作成部１５に予め保持されている。nei_cellは粒子ｉが所属するセルからｘ、ｙ、ｚ方向それぞれに対して＋−方向に探索するセル数を表す。
【００４６】
これにより、セルサイズが小さくなるにつれて探索するセル数が急激に増えて計算負荷が増加する懸念が生じる。しかし、隣接する２６個のセルより外側のセルに対しては、粒子径でソートした粒子の並び順を利用して近傍粒子の探索を効率化することができる。つまり、各セルにはPnum_in_cellmin [cell id]番からPnum_in_cellmax [cell id]番までの粒子が粒子径の大きい順に所属しており、各セル内の粒子が近傍粒子であるか否かを調べる際も粒子径の大きい順に調べて行くことになる。この時、粒子iが所属するセルとその近傍セルの位置関係から上記の（ａ）の条件を満たさない粒子径が求まるのでPnum_in_cellmax [cell id]番までの全ての粒子間距離を調べる必要はなく、（ａ）の条件を満たさない粒子径の粒子番号になったら探索を中止することができる。ここで、粒子iに対して(a)の条件を満たさない粒子径rj とは以下の式を満足するもののことである。
2 ceil{2.0(rj + ri )(1.0 + α)／csize}−1 ＜ dis cell （１０）
dis cell は粒子iの所属するセルから粒子ｊの所属するセルまでのセル数を表し、ｘ、ｙ、ｚ方向のセル数で最小のものを表す。
【００４７】
さらに、近傍粒子表作成部１５は、粒子ｉに対する近傍粒子数n_jli[i]と，上記条件（ａ）を満たすが（ｂ）及び（ｃ）を満たしていない粒子数n_jgi[i]を、それぞれ一次元配列に記憶させておく。
【００４８】
図６の例では、例えば４番の粒子（ｉ＝４）に対してセルを検索した結果５，６，７，８番の粒子が（ａ）の条件を満たし，そのうちの近傍粒子が７，８番で、５，６番の粒子は（ｂ），（ｃ）の条件を満たさない場合、
nei_pn[4][0]=7
nei_pn[4][1]=8
n_jgi[4]=2
n_jli[4]=2
となる。
【００４９】
このように、近傍粒子表作成部１５は、粒子番号変更部１４から粒子情報を受信すると、着目粒子ｉごとに、探索範囲nei_cellのセルを画定し、この探索範囲内に含まれる他の粒子のうち上記（ａ）〜（ｃ）の条件を満たす粒子の粒子番号を近傍粒子表５２の配列nei_pn[i][box_id]に格納する。さらに、近傍粒子数を配列n_jli[i]に格納し、上記条件（ａ）を満たすが（ｂ）及び（ｃ）を満たしていない粒子数、すなわち探索範囲に存在するものの粒子径が着目粒子より小さく、近傍粒子表５２に格納されなかった粒子の数を配列n_jgi[i]に格納する。そして、これらnei_pn[i][box_id] 、n_jli[i]、n_jgi[i]の情報を、粒子情報と共に接触候補リスト作成部１６に送信する。
【００５０】
本実施形態では、近傍粒子表５２には、着目粒子の粒子径以上の粒子のみを近傍粒子として保持している。複数の粒子径をもつ粒子群からなる多成分系では大粒子に対して接触候補になる小粒子の数が膨大になる可能性が有る。つまり、着目粒子より粒子径の小さい粒子とは同時に接触可能な数が非常に多くなる。このため、着目粒子より粒子径の小さい粒子まで近傍粒子に含むとすると、予め近傍粒子表のために用意する配列の数も最大値を考慮して大きくせざるを得ない。その場合、近傍粒子表５２の配列に必要なメモリ量が膨大になり、ＧＰＵのグローバルメモリを食い潰す危険が有る。本実施形態では、自分より大きい粒子しか近傍粒子表５２に格納しないので、近傍粒子表５２のために用意される配列のメモリ使用量を抑制することができる。
【００５１】
接触候補リスト作成部１６は、作業空間５１の全粒子の中から粒子番号（特定情報）の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する。言い換えると、接触候補リスト作成部１６は、作業空間５１の全粒子について、接触している可能性がある粒子ペアからなる接触候補リスト５３を作成する。接触候補リスト５３とは、具体的には、接触している可能性がある２つの粒子ｉ，ｊの粒子番号を格納する一組の一次元配列list_i[Lnum]、list_j[Lnum]である。ここでLnumは接触粒子ペアごとに付される接触候補番号である。
【００５２】
接触候補リスト作成部１６は、近傍粒子表作成部１５から情報（近傍粒子表５２nei_pn[i][box_id] 、近傍粒子数n_jli[i] 、n_jgi[i]、及び粒子情報）を受信すると、まず、図７に示すように、n_jgiのプレフィックス和s_jgiを次式から求める。
【数４】
【００５３】
このプレフィックス和s_jgi[i]は、粒子iの近傍粒子表５２に格納されなかった粒子数、すなわち粒子ｉの周囲の所定の探索範囲（上記（９）式で画定される近傍のセル）に存在する粒子のうち粒子ｉの粒子径より小さい粒子の数n_jgi[i]（接触候補数）の、粒子番号ｉまでの積算値をあらわす一次元配列であって、本実施形態では「積算接触候補数」とよぶ。言い換えると、積算接触候補数s_jgi[i]は、接触候補数n_jgi[i]を粒子番号ｉ順で積算した値である。
【００５４】
接触候補リスト作成部１６は、次に、粒子番号の昇順で着目粒子ｉを順次選択し、各粒子ｉに対して、上記（９）式で画定された粒子ｉの周囲の探索範囲のセルを再度探索して、粒子ｉの周囲の探索範囲に含まれる粒子のうち
ａ）粒子ｉとの粒子間距離がdis_limよりも小さく
ｂ）粒子ｉの径よりも小さく、
ｃ）粒子ｉの径と同じ場合は、初期配置時の粒子番号が粒子ｉよりも小さい
という条件を満たす粒子ｊを抽出する。そして、これらの粒子ｉ，ｊのペアに接触候補番号Lnumを昇順で付与し、それぞれの粒子の粒子番号を１次元配列list_i，list_jに格納して、接触候補リスト５３を作成する。
【００５５】
この条件ａ）〜ｃ）を満たすこととは、上述の近傍粒子表に格納する近傍粒子を抽出する条件（ａ）を満たし（ｂ），（ｃ）を満たさないのと同等である。つまり近傍粒子表５２に格納されない粒子ペア（例えばｉ＝４の着目粒子の場合にはｉ＝５，６の粒子とのペア）が接触候補ペアである。
【００５６】
例えば、図７の例では、接触候補リスト５３の内容は、以下の表１のようになる。
【表１】
つまり、接触候補リスト５３は、作業空間５１の全粒子について、接触している可能性がある粒子ペアlist_i，list_jに番号Lnumを付加した接触候補ペア情報を含むものである。
【００５７】
接線方向の接触力を求める時に必要となるバネの伸縮量は積算値で与えられるため、前ステップで接触していて現ステップでも接触している場合は接触候補ペアリストを更新する時に接線方向バネの伸縮量を引き継ぐ必要が有る。そこで、次に、１ステップ前の接触候補ペアBlist_i，Blist_jを探索し、現在の接触候補ペアlist_i, list_jが１ステップ前にもペアであったか調べる。ペアであったら、図２に示した粒子モデルにおける法線ベクトルと接線方向のバネによる力（バネの伸縮量）を新しい候補リスト番号の配列に格納する。
【００５８】
ここで、１ステップ前及び現在の各パラメータ及び配列の表記は以下のとおりである。
【表２】
【００５９】
粒子番号i=list_i[Lnum], j=list_j[Lnum]の接触候補ペアが有るとき、このペア粒子の１ステップ前の粒子番号はBi=Bpn[i]，Bj=Bpn[j]である。Bs_jgi[Bi−1]〜Bs_jgi[Bi]−1の範囲の接触候補Blist_j[BLnum]を調べBjと一致するか調べる。一致したときはそのリスト番号BLnumを用いて１ステップ前の法線ベクトルと接線方向の接触力BSpring[BLnum]を現在のリスト番号配列Spring[Lnum]に代入する。
【００６０】
この処理について、図８を参照して、現在の接触候補リスト番号Lnum = 2におけるｉ＝4，ｊ＝6番のペア粒子について１ステップ前にペアであったか調べる処理を説明する。
【００６１】
まず、ｉ＝4，ｊ＝6の１ステップ前の粒子番号はBpn[]関数を使ってBi=2=Bpn[4]，Bj=3=Bpn[6]と求まる。ここでBiとBjとの関係は、ｉとｊの関係と同様に、接触候補リスト５３の条件ａ）〜ｃ）を必ず満たしている。なぜなら、条件ａ）〜ｃ）はソートによって変化しない粒子径と初期配置時の粒子番号を使用しているためである。したがってBiに対する接触候補の粒子番号はリスト番号BLnumがBs_jgi[Bi−1]〜Bs_jgi[Bi]−1（この例では１〜１）の範囲のBlist_j[BLnum]に格納されている。そこで、このBLnumの範囲でBlist_j[BLnum] がBjと一致するか調べると、BLnum = 1のときにBlist_j[1]= 3でありBj = 3と一致する。
【００６２】
したがって，現在のリスト番号Lnum=2の候補ペアは1ステップ前もペアであった事がわかり、そのリスト番号はBLnum=1である。そしてBLnum=1の法線ベクトルと接線方向の接触力の値をLnum=2の配列に代入する（Spring[2]=BSpring [1]）ことで、1ステップ前の接触状態を引き継ぐ事ができる。
【００６３】
このように、接触候補リスト作成部１６は、近傍粒子表作成部１５から各種情報報（近傍粒子表５２nei_pn[i][box_id] 、近傍粒子数n_jli[i] 、n_jgi[i]、及び粒子情報）を受信すると、着目粒子ｉごとに、上記探索範囲内に含まれる他の粒子のうち上記ａ）〜ｃ）の条件を満たす粒子ｊとのペアに番号Lnumを付与し、これらの粒子ペアｉ，ｊの粒子番号をそれぞれ接触候補リスト５３の配列list_i[Lnum]，list_j[Lnum]に格納する。また、着目粒子ｉの周辺セルに存在し、粒子ｉの粒子径より小さい粒子の数n_jgi[i]の積算値である積算接触候補数s_jgi[i]を算出する。さらに、現在の接触候補ペアlist_i[Lnum]，list_j[Lnum]のうち、１ステップ前も接触候補だったペアを抽出し、これらのペアの１ステップ前の接触力をBSpring[BLnum]から引き継ぎ、Spring[Lnum]に格納する。そして、導出した接触候補リスト５３list_i[Lnum]，list_j[Lnum]と、近傍粒子表５２nei_pn[i][box_id]とを接触力参照表作成部１７に送信し、接触候補リスト５３list_i[Lnum]，list_j[Lnum]、積算接触候補数s_jgi[i]、１ステップ前から引き継いだ接触候補ペアの接触力Spring[Lnum]の情報を、粒子情報と共に接触判定部１８に送信する。
【００６４】
接触力参照表作成部１７は、box_idから対応する候補リスト番号Lnumを算出する関数Ref[i][box_id]を作成する。本実施形態では、この関数をあらわす二次元配列Ref[i][box_id]を接触力参照表５４と呼ぶ。接触力参照表５４は、図９に示すように、接触候補リスト５３のi粒子がｊ粒子から受ける接触力Fij[Lnum]を任意の一粒子ｉ（図９ではi＝４）ごとに参照先（Lnum）（参照情報）をまとめたものである。特に本実施形態では、接触力参照表５４には、着目粒子ｉに関する近傍粒子表５２に格納されている近傍粒子、すなわち着目粒子ｉより粒子径の大きい粒子との接触力の参照先（Lnum）のみが格納されている。ここで、接触力の参照先（Lnum）とは、後述する粒子間接触力の並進成分及び回転成分を保持する２種類の一次元配列Fijt[Lnum]，Fijr[Lnum]及び後述する（１３）式で用いるβ[Lnum]からなる接触力Fij[Lnum]から構成される接触力ボックスから、該当する接触力を取得するために用いる接触候補ペアリスト番号（接触候補番号）Lnumである。
【００６５】
接触力参照表作成部１７は、接触力参照表５４を作成するためのメモリ領域を予め確保している。確保するメモリ領域は、各粒子と接触する可能性がある近傍粒子の個数を考慮して、各粒子の接触力参照表５４ごとに同一の大きさである。
【００６６】
接触力参照表作成部１７は、具体的には、接触候補リスト作成部１６より接触候補リスト５３list_i[Lnum]，list_j[Lnum]及び近傍粒子表５２nei_pn[i][box_id]を受信すると、接触候補リスト５３のLnum番目の粒子ペアである粒子ｉ（＝list_i[Lnum]）及び粒子ｊ（＝list_j[Lnum]）について、粒子ｊの近傍粒子表５２を探索して粒子ｉが記憶されている近傍粒子表番号ｋ（＝box_id）を取得する。上述したように、接触候補ペアの粒子ｊは、粒子ｉより粒子径が小さいので、粒子ｉは粒子ｊの近傍粒子として粒子ｊの近傍粒子表５２に記憶されているからである。そして、粒子ｊに関する接触力参照表５４のｋ番目の配列Ref[j][k]に、粒子ｉ，ｊの接触候補リスト番号Lnumを記憶させる。
【００６７】
なお、上述の接触候補リスト作成部１６による接触候補リスト５３の作成処理、及び接触力参照表作成部１７による接触力参照表５４の作成処理は、以下のアルゴリズムで実施することができる。以下のアルゴリズム中の「接触候補ペア条件==true」とは、接触候補リスト５３の条件ａ）〜ｃ）を満たすことを意味するものである。
【数５】
【００６８】
本実施形態では、接触力参照表５４には、着目粒子の粒子径以上の粒子のみを保持している。このため、上述した近傍粒子表５２と同じ理由により、接触力参照表５４のために用意される配列のメモリ使用量を抑制することができる。
【００６９】
接触力参照表作成部１７は、導出した接触力参照表５４Ref[i][box_id]を接触力計算部１９に送信する。
【００７０】
接触判定部１８は、接触候補リスト５３を参照して、全てのペア粒子について２つの粒子間の接触判定をすると共に、接触していると判定した場合には粒子間の接触力Fijr、Fijtを算出する。
【００７１】
具体的には、接触判定部１８は、接触候補リスト作成部１６から接触候補リスト５３list_i[Lnum]，list_j[Lnum]、積算接触候補数s_jgi[i]、１ステップ前から引き継いだ接触候補ペアの接触力Spring[Lnum]、及び粒子情報を受信するか、または、接触候補リスト更新部１３から粒子情報を受信すると、接触候補リスト５３を走査して接触候補ペアの粒子ｉ，ｊの接触判定を行う。接触判定は、例えば、粒子ｉと粒子ｊとの粒子間距離を計算し、粒子間距離が両粒子の粒子径の和ｒｉ＋ｒｊ以下となる場合に両粒子が接触しているものと判定する。
【００７２】
接触判定において、両粒子が接触していると判定した場合には、上述したＶｏｉｇｔモデルの（１），（２）式を用いて接触力Ｆｎ，Ｆｔを計算する。なお、上記（２）式において現時点のバネによる接線方向接触力はKtxt=spring[]＋ΔFtとして計算する。spring[]は、上述のように１ステップ前のバネによる接線方向接触力であり、ΔFtは、１ステップ前から現時点までに増加（もしくは減少）したバネによる接線方向接触力である。ΔFtは粒子間の相対速度に１ステップの時間Δｔと接線方向成分のバネの弾性係数Ｋｔを乗じて算出される。法線方向成分及び接線方向成分のバネ７２の弾性係数Ｋｎ，Ｋｔと、粒子間の反発係数ｅは、シミュレーション前に設定される所与の値であり、ダッシュポット７３の粘性減衰係数μｎ，μｔは、所与のＫｎ，Ｋｔ，ｅと、粒子情報に含まれる質量ｍｉ，ｍｊを用いて（３）式から導出する。粒子間の相対変位ｘｎ，ｘｔ及び相対速度ｖｎ，ｖｔは、粒子情報に含まれる粒子座標（位置情報）、並進・回転速度（速度情報）に基づき算出する。Ｆｎ，Ｆｔは、このように取得した各変数及びパラメータを用いて、（１）式及び（２）式に示す演算処理を行って算出する。そして、算出したＦｎ，Ｆｔを用いて、接触力の並進成分Fn +Ft と回転成分Ri ×Ft を計算し、これらの成分を、接触候補ペアリスト番号（接触候補番号）Lnumを要素としてもつ配列Fijt[Lnum]，Fijr[Lnum]（上述の接触力ボックスに含まれる配列）に記憶する。
【００７３】
また、接触判定において、接触していないと判定された場合には、接触力は０なので、接触候補ペアリスト番号Lnumを要素としてもつ配列Fijt[Lnum]，Fijr[Lnum]（上述の接触力ボックス）には０を格納する。
【００７４】
この接触判定部１８による接触力算出処理は、以下のアルゴリズムで表現することができる。
【数６】
【００７５】
このように、接触判定部１８は、接触候補リスト５３が作成／更新される毎に、接触候補ペアリスト番号Lnumに対応付けるように接触力ボックスを作成／更新する。また、接触候補リストが更新されない場合でも、粒子情報が更新された場合には、既存の接触候補リストを用いて、接触力ボックスを作成／更新する。
【００７６】
接触判定および接触力は接触候補ペアリスト番号をスレッド化して計算する。ＧＰＵのグローバルメモリはアクセスパターンがスレッドの順である場合に最も良いアクセス効率が得られる。しかし、ペア粒子の番号はリスト番号に対してランダムであるため、粒子情報が格納されているグローバルメモリをランダムにアクセスしメモリアクセス速度が大きく低下する。一方、テクスチャメモリはキャッシュ機能が有るため、ランダムアクセスによるメモリアクセス速度の低下を抑えることができる。そこで、粒子情報を格納しているグローバルメモリ領域をテクスチャメモリから参照できるように設定した。さらに粒子には所属するセル番号順になるように番号を付けているために接触候補ペア同士の粒子番号が近く、加えて接触候補ペアリストは粒子番号ｉが小さい順にリスト内に並ぶように作られている。そのため、接触候補ペアリスト番号でスレッド化することでキャッシュヒット率が高くなり、粒子間の接触判定および接触力計算のルーチンが高速化される。
【００７７】
接触判定部１８は、作成／更新した接触力ボックスの情報を、積算接触候補数s_jgi[i]及び粒子情報と共に接触力計算部１９に送信する。
【００７８】
接触力計算部１９は、接触力参照表５４と積算接触候補数s_jgiを用いて各粒子ｉに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]を計算する。接触力参照表５４を接触力参照表作成部１７から受信し、積算接触候補数を接触判定部１８から受信すると、まず、着目粒子ｉごとに、接触力参照表５４及び積算接触候補数s_jgiを用いて、この着目粒子と接触している粒子との接触力を、接触力ボックスから抽出する。具体的には、着目粒子ｉより粒子径が大きい粒子については、接触力参照表５４に保持されるボックス番号（Lnum）（参照情報）に基づき、このボックス番号に対応する接触力Fij[Lnum]を接触力ボックスから抽出する。着目粒子ｉより粒子径が小さい粒子については、積算接触候補数s_jgi[]を用いて、接触力テーブルの格納位置を特定して接触力を抽出する。より詳細には、s_jgi[i-1]〜s_jgi[i]−1番のボックス番号（接触候補番号）に対応する接触力Fij[s_jgi[i-1]]〜Fij[s_jgi[i]−1]を接触力ボックスから抽出する。ここで接触力ボックスに格納されている接触力Fij[]とは、具体的には並進成分Fijt[Lnum]、回転成分Fijr[Lnum]、及び後述する（１３）式で用いるβ[Lnum]の３つの１次元配列を一組として保持するものである。
【００７９】
この接触力を接触力ボックスから抽出する処理を図１０を参照して具体的に説明する。図１０は、これまでの図６，９と同様に、着目粒子ｉ＝４の例である。接触力参照表５４には、着目粒子（ｉ＝４）より粒子径が大きい粒子（ｉ＝７，８）との接触力の参照先情報が含まれている。ｉ＝４の着目粒子とｉ＝７の粒子との接触力を格納する接触力ボックスのボックス番号（＝接触候補リスト番号）Lnumは４であり、ｉ＝８の粒子との接触力を格納する接触力ボックスのボックス番号Lnumは７である。接触力計算部１９は、接触力参照表５４に含まれるこのボックス番号を用いて、接触力ボックスから、ｉ＝７の粒子とのとの接触力Fij[４]，ｉ＝８の粒子との接触力Fij[７]を抽出する。
【００８０】
また、図１０に示すように、着目粒子ｉ＝４の積算接触候補数s_jgi[i]は、s_jgi[4]＝３であり、着目粒子の粒子番号より１つ前の粒子（図１０ではi＝３）の積算接触候補数s_jgi[i-1]は、s_jgi[3]＝３である。したがって、着目粒子ｉ＝４より粒子径が小さい粒子（この例ではi＝５，６）については、s_jgi[i-1]〜s_jgi[i]−1番、すなわち１番及び２番のボックス番号に接触力Fij[１]，Fij[２]を接触力ボックスから抽出する。
【００８１】
次に、着目粒子iごとに、接触力ボックスから抽出した接触力Fij[]（＝Fijt[]、Fijr[]、β[]）を用いて、粒子iに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]を計算する。
【００８２】
ｉ番の粒子に働く全接触力（並進成分）Ftall[i]は、接触候補リスト番号がs_jgi[i-1]〜s_jgi[i]−1番の接触力はそのままの符号で足し合わせる。さらに、接触力参照表５４（Ref）を用いて参照表番号0〜n_jli−1番の表から参照される接触力は逆符号にして足し合わせる。そして次式のように表される。
【数７】
【００８３】
ｉ番の粒子に働く全接触力（回転成分）Frall[i]も同様に足し合わせる。接触候補リスト番号がs_jgi[i-1]〜s_jgi[i]−1番の接触力はそのまま参照し、接触力参照表５４（Ref）を用いて参照表番号0〜n_jli−1番の表から参照される接触力は逆符号にしてからβを掛けてから足し合わせ、次式のように表される。ここで、βはｊ粒子の中心から接触点までの距離をｉ粒子の中心から接触点までの距離で割った値である。
【数８】
【００８４】
ここで、s_jgi[i-1]〜s_jgi[i]−1 番の接触力は粒子ｉが粒子ｊから受ける接触力であるのに対して、0〜n_jli−1番の接触力参照表５４で参照される接触力は粒子ｉが粒子ｊに与える接触力である。そのため， 0〜n_jli−1番の接触力参照表５４に対して、接触力の並進成分については符号を反対にしてから足し合わせ、接触力の回転成分についてはβ[box id]を掛けてから足し合わせる。この様に、粒子ｉが粒子ｊから受ける力を接触力参照表５４Ref[i][box id]を使わずにs_jgi[i]を用いて直接、接触力の入った配列Fijt[]，Fijr[]を参照することで、接触力参照表５４のメモリ使用量を増加させずに、大粒子に対して膨大な数の接触候補ペアとなる小粒子から受ける接触力を足し合わせることが可能になった。
【００８５】
コンピュータを用いた演算においては、具体的には粒子番号ｉをスレッド化し、box_idでループを回してFtall，Frallを求める。ここでFijへのメモリアクセスはランダムアクセスになるのでテクスチャメモリを使用する。上述のようにテクスチャメモリはキャッシュヒット率が高いので、グローバルメモリを直接アクセスするよりは高速に処理できる。なお接触力参照表５４に関しては配列Ref[A][B]のBの要素サイズを16の倍数になるように配列宣言することで連続アクセスが可能となる。
【００８６】
接触力計算部１９は、上述のように算出した、各粒子ｉに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]を、粒子情報と共に粒子情報更新部２０に送信する。
【００８７】
粒子情報更新部２０は、接触力計算部１９から各粒子ｉに働く全接触力の並進成分Ftall[i]、及び回転成分Frall[i]と、粒子情報とを受信すると、各粒子の位置、速度、接触力計算部１９から受信した接触力に基づき、粒子の座標（位置情報）と、並進速度及び回転速度（速度情報）を算出し、これらの算出した座標、並進速度及び回転速度を用いて、粒子情報保持部１１に保持されている粒子情報を更新する。具体的には、並進の運動方程式（上記（４）式）をleap-flog法を用いて解くことで新しい時刻の並進速度ならびに粒子座標を求める。さらに、回転の運動方程式（上記（５）式）も同様にleap-flog 法を用いて解き、ここで得られる回転速度を予測値とする。この予測される回転速度に対して転がり摩擦力を計算し、接触力による回転力と転がり摩擦力を考慮した回転の運動方程式を解くことにより新しい時刻の回転速度を得る。
【００８８】
粒子情報更新部２０は、算出した各粒子の位置座標、並進速度、回転速度を用いて、粒子情報の座標、並進速度、回転速度を更新し、更新した粒子情報を粒子情報保持部１１に送信して粒子情報保持部１１の粒子情報を最新のものに更新する。
【００８９】
次に、図１１を参照して、粒子シミュレーション装置１０において実行される処理について説明すると共に、本実施形態に係る粒子シミュレーション方法について説明する。図１１は本実施形態の粒子シミュレーション装置１０において実行される処理を示すフローチャートである。粒子情報保持部１１に入力され保持されている粒子情報が、粒子情報取得部１２により粒子情報保持部１１から取得される（Ｓ１０１：入力ステップ）。接触候補リスト更新判定部１３により、接触候補リストの更新が必要か否かが判定される（Ｓ１０２）。更新判定は、例えば、粒子半径をｒとして、粒子の積算移動距離がｒαよりも大きい粒子が存在したときに更新する。ここでαは更新頻度を調整するパラメータであり、例えばα＝０．２である。更新が必要と判定された場合にはステップＳ１０３に移行し、更新が不要と判定された場合にはステップＳ１０８に移行する。
【００９０】
次に、粒子番号変更部１４により粒子番号が変更される（Ｓ１０３：付与ステップ）。粒子番号変更部１４は、粒子を所属するセル番号順に並べ、その順に粒子番号を付け直し、さらに、同じセル番号の粒子について粒子径順に並べ直し、その順に粒子番号を付け直す。続いて、近傍粒子表作成部１５により、作業空間５１内の個々の粒子について、この粒子の近傍に位置する近傍粒子の情報を含む近傍粒子表５２が作成される（Ｓ１０４）。近傍粒子表５２に近傍粒子として格納する条件は、上述の条件（ａ）〜（ｃ）である。
【００９１】
次に、接触候補リスト作成部１６により、近傍粒子表５２の条件（ａ）を満たすが（ｂ）及び（ｃ）を満たしていない粒子数n_jgi[i] のプレフィックス和である積算接触候補数s_jgi[i]が算出され（Ｓ１０５：積算接触候補数算出ステップ）、次いで、作業空間５１の全粒子について、接触している可能性がある粒子ペアの情報を含む接触候補リスト５３が作成される（Ｓ１０６：接触候補選択ステップ）。接触候補リスト５３に格納する条件は、上述の条件ａ）〜ｃ）である。接触力参照表作成部１７により、接触力参照表５４が作成される（Ｓ１０７：接触力参照表作成ステップ）。
【００９２】
次に、接触判定部１８により、接触候補リスト５３を参照して、２つの粒子間の接触判定が行われ（Ｓ１０８）、接触していると判定された場合には粒子間の接触力Fijr、Fijtが（１），（２）式を用いて算出され、接触力ボックスの配列に格納される（Ｓ１０９：接触力算出ステップ）。また、粒子が接触していないと判定された場合には接触力Fijr、Fijtを０として、接触力ボックスの配列に格納される。
【００９３】
接触力計算部１９により、接触力参照表５４及び積算接触候補数s_jgi[i]を参照して、着目粒子と接触する他の粒子が選択され、接触力ボックスから対応する配列Fijr、Fijtから接触力が取得される。そして（１２），（１３）式を用いて接触力が総和演算されて、各粒子に働く接触力Fr、Ftが計算される（Ｓ１１０：接触力総和演算ステップ）。粒子情報更新部２０により、ステップＳ１０９において算出された接触力、現在の粒子情報に基づいて、各粒子の座標（位置情報）と、並進速度及び回転速度（速度情報）とが算出され、これらの算出した座標、並進速度及び回転速度を用いて、粒子情報保持部１１に保持されている粒子情報が更新される（Ｓ１１１：粒子情報更新ステップ）。
【００９４】
そして、粒子シミュレーションが終了したか否かが判定される（Ｓ１１２）。終了判定は、例えば所定時間が経過した場合にシミュレーションが終了したものと判定される。粒子シミュレーションが終了していないと判定された場合にはステップＳ１０１に戻り、次のステップに入る。粒子シミュレーションが終了したと判定された場合には処理を終了する。
【００９５】
次に、図１２を参照して、本実施形態の粒子シミュレーション装置１０によるメモリ使用量について説明する。図１２には、粒子モデルの最大粒子径と最小粒子径との比を変化させた場合の本実施形態（実線）及び比較例（破線）のメモリ使用量を示している。なお、粒子モデルの粒子数は１，０００，０００個とし、最大粒子径の粒子（最大粒子）と最小粒子径の粒子（最小粒子）とが混在する二成分系とした。比較例は、近傍粒子表５２及び接触力参照表５４に着目粒子ｉの粒子径より小さい粒子の情報も格納する点が本実施形態と異なるものである。図１２の横軸は、「最大粒子径／最小粒子径」を表し、縦軸は「メモリ資料量（MB）」を表している。
【００９６】
比較例の場合、近傍粒子表５２nei_pn[i][box_id]及び接触力参照表５４Ref[i][box_id]に必要なメモリ使用量は、１データに４バイト必要として、粒子数×最大接触候補粒子数×４byteとなる。ここで、最大接触候補粒子数は、近傍粒子表５２及び接触力参照表５４に格納する近傍粒子の最大個数であって、最大粒子の表面を最小粒子が埋め尽くす状態を考えればよく、この場合、
最大接触候補粒子数＝４π(ｒmax+ｒmin)２／πｒmin２ （１４）
となる。ここで、ｒｍａｘは最大粒子の粒子径であり、ｒｍｉｎは最小粒子の粒子径である。つまり、比較例の場合、図１２に示すように、最小粒子径に対して最大粒子径が大きくなるほど、つまり最大粒子径／最小粒子径が大きくなるほど、最大接触候補粒子数が増大し、この結果、メモリ使用量が増大する。例えば、粒子径比１：１０の条件（最大粒子径／最小粒子径＝１０）では、一粒子当たりの接触候補粒子数は最大で５００個程度に増加するため、必要なメモリ量で２ＧＢにもなる。
【００９７】
一方、本実施形態の場合、着目粒子より粒子径の小さい粒子は近傍粒子として考慮されず、近傍粒子表５２及び接触力参照表５４に格納されないので、最大接触候補粒子数を算出する上記（１４）式は成り立たず、本実施形態の最大接触候補粒子数は、最大と最小の粒子径に依存しない。本実施形態では、最大接触候補粒子数は、上述した近傍粒子の判定条件（ａ）において粒子間距離の閾値であるdis_lim = (rｉ + rｊ)(1.0+α)に含まれるパラメータαのみに依存し、粒子径分布には依存せず一定となる。α＝０．２の場合は、最大接触候補粒子数は如何なる粒子径分布条件においても１６も有れば十分であるので、ここでは最大接触候補粒子数を１６としてメモリ使用量を計算した。その結果、図１２に示すように、本実施形態のメモリ使用量は、粒子数×最大接触候補粒子数×４byte ＝1000000×16×４＝64MBの一定値となった。
【００９８】
本実施形態に係る粒子シミュレーション装置１０及び粒子シミュレーション方法によれば、作業空間５１内の粒子のそれぞれについて、当該粒子ｉの近傍に位置する、この粒子ｉと接触する可能性の高い他の粒子のうち当該粒子ｉの粒子径より小さい粒子については、粒子の数を示す積算接触候補数s_jgi[i]のみを保持しておく。そして、接触判定部１８が、この積算接触候補数s_jgi[i]を利用して接触力テーブルの格納位置を特定して当該粒子と他の粒子との接触力を抽出する。これにより、接触力参照表５４には、当該粒子ｉの近傍に位置する他の粒子のうち当該粒子ｉの粒子径より小さい粒子に関する情報を含める必要がなく、当該粒子の粒子径より大きい粒子に関する情報のみを保持すればよくなる。粒子径分布を有する粒子群をシミュレーション対象とする場合、自分より粒子径の小さい粒子まで接触を考慮すると、同時に接触可能な粒子の数が増大するため、予め接触力参照表５４のために確保すべきメモリ領域を多くしなければならない。これに対し、上記構成により、自分より粒子径の大きい粒子との接触のみを考慮すればよく、接触を考慮する粒子数を減らすことができ、接触力参照表５４を作成するために予め割り当てるメモリ領域を少なくすることが可能となり、この結果、メモリ使用量を抑制することができる。このようにメモリ使用量が抑制されることにより、ＧＰＵなどの共有メモリ型並列計算機を利用して粒子シミュレーションを実行することが可能となる。
【００９９】
また、粒子番号変更部１４が、作業空間５１内の全粒子について、セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した粒子番号を付与し、各粒子をセル番号順、さらには同一セル内では粒子径順にソーティングする。このため、接触力参照表を作成する際に、粒子番号に基づき粒子径が小さい粒子を容易に認識することができ、粒子径が小さい粒子については接触力参照表にいれるか否かの判定を行う必要がなくなる。この結果、接触力参照表の作成において計算効率を向上させることができる。
【０１００】
以上、本発明について好適な実施形態を挙げて説明したが、本発明は上記実施形態に限定されるものではない。例えば、上記実施形態では、粒子番号変更部１４は、（１）粒子を所属するセル番号順に並べ、さらに、（２）同じセル番号の粒子について粒子径順に並べ直す、というソーティングを行っているが、（１）粒子を所属するセル番号順に並べる処理のみを行うよう構成してもよい。
【０１０１】
また、上記実施形態では、本発明の適用対象としてＧＰＵを挙げたが、ＧＰＵと同様の共有メモリ型のＳＩＭＤ（Single Instruction Multiple Data）型並列処理を行うハードウェア（例えばベクトルプロセッサ）、すなわち共有メモリ型並列計算機にも適用可能である。
【符号の説明】
【０１０２】
１０…粒子シミュレーション装置、１１…粒子情報保持部（入力手段、粒子情報保持手段）、１２…粒子情報取得部、１３…接触候補リスト更新判定部、１４…粒子番号変更部（付与手段）、１５…近傍粒子表作成部、１６…接触候補リスト作成部（接触候補選択手段、積算接触候補数算出手段）、１７…接触力参照表作成部（接触力参照表作成手段）、１８…接触判定部（接触力算出手段）、１９…接触力計算部（接触力総和演算手段）、２０…粒子情報更新部（算出手段、粒子情報更新手段）、５１…作業領域、５２…近傍粒子表、５３…接触候補リスト、５４…接触力参照表。
【特許請求の範囲】
【請求項１】
作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、
前記粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力する入力手段と、
前記粒子情報を保持する粒子情報保持手段と、
前記粒子のそれぞれについて、前記位置情報に応じた順序で特定する特定情報を付与する付与手段と、
前記粒子の中から前記特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択手段と、
前記選択候補選択手段により選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の前記位置情報及び前記速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出手段と、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を前記接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成手段と、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出手段と、
前記粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を前記接触力参照表の参照情報を用いて前記接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、前記積算接触候補数を用いて前記接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算手段と、
前記接触力総和演算手段により計算された粒子ごとの接触力に基づいて、前記粒子の新たな位置情報及び速度情報を算出する算出手段と、
前記算出手段により算出された新たな位置情報及び速度情報を用いて、前記粒子情報保持手段に保持される粒子情報を更新する粒子情報更新手段と、
を備えることを特徴とする粒子シミュレーション装置。
【請求項２】
前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記粒子情報には、前記位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、
前記付与手段が、前記粒子のそれぞれについて、前記セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することを特徴とする、請求項１に記載の粒子シミュレーション装置。
【請求項３】
前記接触候補選択手段が、前記粒子郡の中から前記特定情報の昇順で着目粒子を選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との組と、この組に昇順で付される接触候補番号とからなる接触候補ペア情報を、前記作業空間内の全粒子について含む接触候補リストを作成し、
前記接触力算出手段が、前記接触候補リストに含まれる前記接触候補ペア情報のそれぞれについて、該接触候補ペア情報に関する２つの粒子の間の接触力を、当該粒子の前記位置情報及び前記速度情報に基づき算出し、算出した接触力と前記接触候補番号とを関連付けた接触力情報を接触力テーブルに格納し、
前記接触力参照表作成手段が、前記粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の大きい他の粒子との接触力を前記接触力テーブルから参照するための参照情報として、該他の粒子及び当該粒子に関する前記接触候補番号を、該他の粒子及び当該粒子と関連付けて保持する接触力参照表を作成し、
前記積算接触候補数算出手段が、前記粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の小さい他の粒子の数を示す接触候補数を前記特定情報順で積算した値である積算接触候補数s_jgi[i]（ｉは当該粒子の特定情報を示す）を算出し、
前記接触力総和演算手段が、前記粒子のそれぞれについて、前記接触力参照表の前記参照情報を用いて接触力を前記接触力テーブルから抽出し、また、前記積算接触候補数s_jgi[i]を用いて前記接触候補番号がs_jgi[i-1]〜s_jgi[i]‐1の接触力を前記接触力テーブルから抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する、
ことを特徴とする請求項２に記載の粒子シミュレーション装置。
【請求項４】
作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置における粒子シミュレーション方法であって、
前記粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力し、この粒子情報を粒子情報保持手段に保持する入力ステップと、
前記粒子のそれぞれについて、前記位置情報に応じた順序で特定する特定情報を付与する付与ステップと、
前記粒子の中から前記特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択ステップと、
前記接触候補選択ステップにおいて選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の前記位置情報及び前記速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出ステップと、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を前記接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成ステップと、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出ステップと、
前記粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を前記接触力参照表の参照情報を用いて前記接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、前記積算接触候補数を用いて前記接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算ステップと、
前記接触力総和演算ステップにおいて計算された粒子ごとの接触力に基づいて、前記粒子の新たな位置情報及び速度情報を算出する算出ステップと、
前記算出ステップにおいて算出された新たな位置情報及び速度情報を用いて、前記粒子情報保持手段に保持される粒子情報を更新する粒子情報更新ステップと、
を含むことを特徴とする粒子シミュレーション方法。
【請求項５】
前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記粒子情報には、前記位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、
前記付与ステップが、前記粒子のそれぞれについて、前記セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することを特徴とする、請求項４に記載の粒子シミュレーション方法。
【請求項１】
作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、
前記粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力する入力手段と、
前記粒子情報を保持する粒子情報保持手段と、
前記粒子のそれぞれについて、前記位置情報に応じた順序で特定する特定情報を付与する付与手段と、
前記粒子の中から前記特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択手段と、
前記選択候補選択手段により選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の前記位置情報及び前記速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出手段と、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を前記接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成手段と、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出手段と、
前記粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を前記接触力参照表の参照情報を用いて前記接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、前記積算接触候補数を用いて前記接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算手段と、
前記接触力総和演算手段により計算された粒子ごとの接触力に基づいて、前記粒子の新たな位置情報及び速度情報を算出する算出手段と、
前記算出手段により算出された新たな位置情報及び速度情報を用いて、前記粒子情報保持手段に保持される粒子情報を更新する粒子情報更新手段と、
を備えることを特徴とする粒子シミュレーション装置。
【請求項２】
前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記粒子情報には、前記位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、
前記付与手段が、前記粒子のそれぞれについて、前記セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することを特徴とする、請求項１に記載の粒子シミュレーション装置。
【請求項３】
前記接触候補選択手段が、前記粒子郡の中から前記特定情報の昇順で着目粒子を選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との組と、この組に昇順で付される接触候補番号とからなる接触候補ペア情報を、前記作業空間内の全粒子について含む接触候補リストを作成し、
前記接触力算出手段が、前記接触候補リストに含まれる前記接触候補ペア情報のそれぞれについて、該接触候補ペア情報に関する２つの粒子の間の接触力を、当該粒子の前記位置情報及び前記速度情報に基づき算出し、算出した接触力と前記接触候補番号とを関連付けた接触力情報を接触力テーブルに格納し、
前記接触力参照表作成手段が、前記粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の大きい他の粒子との接触力を前記接触力テーブルから参照するための参照情報として、該他の粒子及び当該粒子に関する前記接触候補番号を、該他の粒子及び当該粒子と関連付けて保持する接触力参照表を作成し、
前記積算接触候補数算出手段が、前記粒子のそれぞれについて、当該粒子と接触している可能性のある、当該粒子より粒子径の小さい他の粒子の数を示す接触候補数を前記特定情報順で積算した値である積算接触候補数s_jgi[i]（ｉは当該粒子の特定情報を示す）を算出し、
前記接触力総和演算手段が、前記粒子のそれぞれについて、前記接触力参照表の前記参照情報を用いて接触力を前記接触力テーブルから抽出し、また、前記積算接触候補数s_jgi[i]を用いて前記接触候補番号がs_jgi[i-1]〜s_jgi[i]‐1の接触力を前記接触力テーブルから抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する、
ことを特徴とする請求項２に記載の粒子シミュレーション装置。
【請求項４】
作業空間内で粒子径分布を有する粒子群について、他の粒子との接触力に基づく位置・速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置における粒子シミュレーション方法であって、
前記粒子群の各粒子ごとの位置情報及び速度情報を含む粒子情報を入力し、この粒子情報を粒子情報保持手段に保持する入力ステップと、
前記粒子のそれぞれについて、前記位置情報に応じた順序で特定する特定情報を付与する付与ステップと、
前記粒子の中から前記特定情報の昇順で着目粒子を順次選択し、該着目粒子と接触している可能性のある、当該着目粒子より粒子径の小さい他の粒子との接触候補ペアを選択する接触候補選択ステップと、
前記接触候補選択ステップにおいて選択された接触候補ペアのそれぞれについて、該接触候補ペアの２つの粒子の間の接触力を、これらの粒子の前記位置情報及び前記速度情報に基づき算出し、算出した接触力の情報をこの接触力に対応する接触候補ペアの順番に応じて接触力テーブルに格納する接触力算出ステップと、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より大きい粒子との接触力を前記接触力テーブルから参照するための参照情報を含む接触力参照表を作成する接触力参照表作成ステップと、
前記粒子のそれぞれについて、当該粒子の近傍に位置する他の粒子のうち当該粒子の粒子径より小さい粒子の数を示す接触候補数の積算値である積算接触候補数を算出する積算接触候補数算出ステップと、
前記粒子のそれぞれについて、当該粒子より粒子径の大きい粒子との接触力を前記接触力参照表の参照情報を用いて前記接触力テーブルから抽出し、当該粒子より粒子径の小さい粒子との接触力を、前記積算接触候補数を用いて前記接触力テーブルの格納位置を特定して抽出し、これらの抽出した接触力を用いて粒子ごとの接触力を総和演算する接触力総和演算ステップと、
前記接触力総和演算ステップにおいて計算された粒子ごとの接触力に基づいて、前記粒子の新たな位置情報及び速度情報を算出する算出ステップと、
前記算出ステップにおいて算出された新たな位置情報及び速度情報を用いて、前記粒子情報保持手段に保持される粒子情報を更新する粒子情報更新ステップと、
を含むことを特徴とする粒子シミュレーション方法。
【請求項５】
前記作業空間が複数のセルに分割され、各セルにはセル番号が付されており、
前記粒子情報には、前記位置情報に基づき特定される、当該粒子が配置されるセルのセル番号が含まれており、
前記付与ステップが、前記粒子のそれぞれについて、前記セル番号に応じた順序で特定し、さらに同一のセルに配置される粒子間で、粒子径に応じた順序で特定した特定情報を付与することを特徴とする、請求項４に記載の粒子シミュレーション方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【公開番号】特開２０１１−１４５９４３（Ｐ２０１１−１４５９４３Ａ）
【公開日】平成２３年７月２８日（２０１１．７．２８）
【国際特許分類】
【出願番号】特願２０１０−７１４２（Ｐ２０１０−７１４２）
【出願日】平成２２年１月１５日（２０１０．１．１５）
【出願人】（５０４１９４８７８）独立行政法人海洋研究開発機構 (110)
【公開日】平成２３年７月２８日（２０１１．７．２８）
【国際特許分類】
【出願日】平成２２年１月１５日（２０１０．１．１５）
【出願人】（５０４１９４８７８）独立行政法人海洋研究開発機構 (110)
[ Back to top ]