JP7244757B2 - 情報処理装置、粒子シミュレーション方法および粒子シミュレーションシステム - Google Patents

Description

本発明は情報処理装置、粒子シミュレーション方法および粒子シミュレーションシステムに関する。

物体の挙動を解析するシミュレーション方法の１つに、離散要素法（ＤＥＭ：Discrete Element Method）がある。離散要素法は、個別要素法（Distinct Element Method）と言うこともある。離散要素法は、解析対象の物体を円や球などの粒子の集合としてモデル化し、粒子間の接触および滑動を考慮して粒子の運動を追跡する。離散要素法は、土砂や薬品やプリンタトナーなどの粉体の挙動解析に利用されることがある。

離散要素法は、ある粒子が接触している他の粒子を検出し、他の粒子から受ける接触力を算出し、接触力に応じて粒子の所定時間後の位置を算出する。シミュレーション上の時刻を所定時間刻みで進めて、接触している他の粒子の検出と接触力の算出と位置の更新を繰り返すことで、時間軸に沿った粒子の位置の変化をシミュレートできる。

なお、解析対象領域を複数のセルに分割し、着目する粒子と接触しているか否か判定する相手粒子を、着目する粒子が属するセルと隣接する隣接セルに属する粒子に限定することで、離散要素法の計算量を削減する解析装置が提案されている。また、着目する粒子と所定時間内に接触する可能性がある粒子を示す近傍リストを生成し、各時刻に接触しているか否か判定する相手粒子を近傍リストの粒子に限定することで、離散要素法の計算量を削減するシミュレーションシステムが提案されている。また、粒子を複数のノードに割り振って離散要素法を並列実行するシミュレーションシステムが提案されている。

特開２０１０－７２３７９号公報 国際公開第２０１２／０３４１７６号 国際公開第２０１４／０９４４１０号

離散要素法では、ある粒子が他の粒子から受ける接触力を算出するにあたり、２つの粒子が接触している間の累積の変位量を使用することがある。例えば、接触力を法線方向と接線方向に分解して算出する場合に、接線方向の接触力が、接触開始から現時刻までの間に接線方向にずれた相対的な移動量の累積値に依存することがある。よって、シミュレーション上の時刻を進めるにあたって、継続して接触している粒子のペア毎に前時刻の変位量を現時刻に複写（コピー）する複写処理が発生する。

ここで、ある粒子から見て前時刻で接触しておりかつ現時刻でも接触している他の粒子を特定し、特定した他の粒子の変位量を複写する複写処理は、計算量が大きい。例えば、複写処理の単純な方法として、各粒子の現時刻の位置に基づいて、ある粒子と現時刻で接触している他の粒子を抽出し、抽出した他の粒子を、前時刻で接触していた粒子を示す接触データの中から検索する方法が考えられる。この方法では、前時刻の接触データを何度も走査することになり、比較処理の計算量が大きくなる。

そこで、複写処理の手順を変更して計算量を削減することが考えられる。しかし、粒子が配置され得る全体領域を分割して、複数のノードが異なる分割領域を担当する並列粒子シミュレーションでは、複写処理の手順の変更が容易ではないという問題がある。

並列粒子シミュレーションの場合、ある分割領域に属する粒子が、隣接する他の分割領域に属する他の粒子と境界を挟んで接触することがある。そのため、複数のノードは、境界付近の粒子のデータを相互に通知し合うことになる。ただし、各ノードが粒子の移動を効率的に追跡できるのは当該ノードが担当する分割領域であり、隣接する他の分割領域については粒子の移動を効率的に追跡できるわけではない。そのため、担当する分割領域に相手粒子が属する場合と隣接する他の分割領域に相手粒子が属する場合とを区別せずに、統一的な方法で複写処理を行おうとすると、計算量の削減が難しい。

１つの側面では、本発明は、粒子シミュレーションの計算量を削減できる情報処理装置、粒子シミュレーションシステムおよび粒子シミュレーション方法を提供することを目的とする。

１つの態様では、記憶部と処理部とを有する情報処理装置が提供される。記憶部は、第１の領域に属する第１の粒子に対して、第１の時刻において第１の粒子と接触している接触粒子と、第１の時刻までの接触中の間の接触粒子に対する第１の粒子の累積の変位量とを示す第１の接触データを記憶する。処理部は、第１の時刻より後の第２の時刻において第１の領域に属する粒子の位置を示す第１の位置データを算出し、第２の時刻において第１の領域と隣接する第２の領域に属する粒子の位置を示す第２の位置データを他の情報処理装置から取得する。処理部は、第１の接触データと第１の位置データの算出結果とに基づいて、第１の時刻において第１の粒子と接触しており第１の領域に属し、第２の時刻において第１の領域に属する第２の粒子を検出する。処理部は、第１の位置データに基づいて、第２の時刻において第１の粒子と第２の粒子とが接触しているか判定し、接触している場合、第１の接触データから第２の時刻に対応する第２の接触データに、第２の粒子の変位量を複写する。処理部は、第１の位置データと第２の位置データとに基づいて、第２の粒子以外の第１の領域または第２の領域に属する粒子であって、第２の時刻において第１の粒子と接触している第３の粒子を検出する。処理部は、第１の接触データから第３の粒子を検索し、第３の粒子が第１の接触データに登録されている場合、第１の接触データから第２の接触データに第３の粒子の変位量を複写する。

また、１つの態様では、コンピュータが実行する粒子シミュレーション方法が提供される。また、１つの態様では、粒子シミュレーションシステムが提供される。

１つの側面では、粒子シミュレーションの計算量を削減できる。

第１の実施の形態の情報処理システムの例を説明する図である。 第２の実施の形態の情報処理システムの例を示す図である。 計算ノードのハードウェア例を示すブロック図である。 離散要素法における粒子モデルの例を示す図である。 シミュレーション対象領域の分割例を示す図である。 周辺粒子の移動例を示す第１の図である。 周辺粒子の移動例を示す第２の図である。 接線方向の変位量の複写例を示す図である。 計算ノードの機能例を示すブロック図である。 結果テーブルの例を示す図である。 粒子テーブルおよび粒子番号対応テーブルの例を示す図である。 接触粒子リストの例を示す図である。 離散要素法の手順例を示すフローチャートである。 接触粒子リスト生成の手順例を示すフローチャートである。 接触粒子リスト生成の手順例を示すフローチャート（続き１）である。 接触粒子リスト生成の手順例を示すフローチャート（続き２）である。

以下、本実施の形態を図面を参照して説明する。
［第１の実施の形態］
第１の実施の形態を説明する。

図１は、第１の実施の形態の情報処理システムの例を説明する図である。
第１の実施の形態の情報処理システムは、離散要素法を用いて粉体の動きをシミュレートする粒子シミュレーションを行う。第１の実施の形態の情報処理システムは、情報処理装置１０，１０－１などの複数の情報処理装置を含む。第１の実施の形態の情報処理システムを、並列処理システム、シミュレーションシステム、シミュレータなどと言うこともできる。粒子シミュレーションの並列処理では、粒子が配置され得るシミュレーション上の空間が分割され、異なる情報処理装置に異なる領域が割り当てられる。情報処理装置１０，１０－１などの各情報処理装置は、割り当てられた領域（対象領域）に属する粒子の移動をシミュレートする。各情報処理装置を、コンピュータ、ノード、シミュレーション装置などと言うこともできる。複数の情報処理装置は、ネットワークに接続され、相互にデータを送信できる。後述するように、情報処理装置１０は領域１３を担当し、情報処理装置１０－１は領域１３－１を担当する。

情報処理装置１０は、記憶部１１および処理部１２を有する。記憶部１１は、ＲＡＭ（Random Access Memory）などの揮発性の半導体メモリでもよいし、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性ストレージでもよい。処理部１２は、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、ＤＳＰ（Digital Signal Processor）などのプロセッサである。ただし、処理部１２は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの特定用途の電子回路を含んでもよい。プロセッサは、ＲＡＭなどのメモリ（記憶部１１でもよい）に記憶されたプログラムを実行する。複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」と言うことがある。

情報処理装置１０は、離散要素法に従って、着目する粒子と現時刻で接触している接触粒子を判定し、着目する粒子が接触粒子から受ける接触力を算出し、接触力に基づいて着目する粒子の次時刻の位置を算出することを繰り返す。そこで、情報処理装置１０は、シミュレーション上の各時刻において、接触粒子を示す接触データを生成する。

接触力の算出では、着目する粒子と接触粒子とが接触している間における累積の変位量が使用される。変位量は、例えば、着目する粒子が接触粒子に対して接線方向に移動した相対的な移動量を、接触開始から累積したものである。そこで、情報処理装置１０は、シミュレーション上の各時刻の変位量を接触データに記録しておく。このとき、接触が継続している接触粒子については、前時刻の接触データから現時刻の接触データに変位量を複写し、複写した変位量に現時刻の増分を加えることになる。

記憶部１１は、接触データ１４（第１の接触データ）を記憶する。記憶部１１は、後述する位置データ１５（第１の位置データ）、位置データ１５－１（第２の位置データ）および接触データ１６（第２の接触データ）を更に記憶してもよい。

接触データ１４は、情報処理装置１０が担当する領域１３に属する粒子Ｊ１（第１の粒子）に対する接触データであって、シミュレーション上の時刻ｔ１（第１の時刻）の接触データである。接触データ１６は、時刻ｔ１より後のシミュレーション上の時刻ｔ２（第２の時刻）の粒子Ｊ１に対する接触データである。時刻ｔ２は、例えば、時刻ｔ１から所定の単位時間後の時刻である。接触データ１４は、時刻ｔ１において粒子Ｊ１と接触している接触粒子と、時刻ｔ１までの接触中における接触粒子に対する粒子Ｊ１の変位量とを示す。接触データ１６は、時刻ｔ２において粒子Ｊ１と接触している接触粒子と、時刻ｔ２までの接触中における接触粒子に対する粒子Ｊ１の変位量とを示す。例えば、接触データ１４，１６は、接触粒子の識別情報と変位量とを含む。複数の接触粒子が存在する場合、例えば、複数の識別情報と複数の変位量とが対応付けられて列挙される。

位置データ１５は、時刻ｔ２において情報処理装置１０が担当する領域１３に属する粒子の位置を示す。位置データ１５－１は、時刻ｔ２において情報処理装置１０－１が担当する領域１３－１に属する粒子の位置を示す。領域１３－１は領域１３と隣接する。位置データ１５－１は、領域１３－１に属する粒子のうち、領域１３と領域１３－１の間の境界から所定の距離以内にある一部の粒子のみを示してもよい。位置データ１５は情報処理装置１０で生成され、位置データ１５－１は情報処理装置１０－１で生成される。

処理部１２は、シミュレーション上の時刻を時刻ｔ１から時刻ｔ２に進める。その際、処理部１２は、位置データ１５を生成する。例えば、処理部１２は、接触データ１４などの時刻ｔ１の接触データに基づいて、領域１３に属する各粒子が受ける接触力を算出し、算出した接触力と時刻ｔ１における領域１３の各粒子の位置を示す位置データとに基づいて、時刻ｔ２における位置を算出する。また、処理部１２は、情報処理装置１０－１から受信された位置データ１５－１を取得する。位置データ１５－１は、位置データ１５と同様の方法により情報処理装置１０－１によって生成される。

時刻を時刻ｔ１から時刻ｔ２に進めることで、粒子が領域１３と領域１３－１の境界を跨いで移動することもある。時刻ｔ１において領域１３に属する粒子は、時刻ｔ２において領域１３に留まっていることもあるし、時刻ｔ２において領域１３－１に移動していることもある。また、時刻ｔ１において領域１３－１に属する粒子は、時刻ｔ２において領域１３－１に留まっていることもあるし、時刻ｔ２において領域１３に移動していることもある。各粒子は、移動先の領域を担当する情報処理装置により管理される。

次に、処理部１２は、時刻ｔ２の接触データ１６を生成する。その際、処理部１２は、接触データ１４と位置データ１５の算出結果とに基づいて、時刻ｔ１において粒子Ｊ１と接触しかつ領域１３に属し、時刻ｔ２において領域１３に属する第２の粒子を検出する。

例えば、処理部１２は、時刻ｔ１と時刻ｔ２とで各粒子に異なる識別情報を割り当てる。ただし、処理部１２は、位置データ１５を算出する際に、領域１３に属する粒子については、時刻ｔ１の識別情報と時刻ｔ２の識別情報との間の対応関係を示す対応データを生成しておく。情報処理装置１０が粒子の移動を追跡する領域は領域１３のみであるため、対応データには、時刻ｔ１と時刻ｔ２の両方で領域１３に属する粒子が登録される。対応データには、領域１３から領域１３－１に移動した粒子、領域１３－１から領域１３に移動した粒子、および、領域１３－１に留まる粒子は登録されない。その場合、処理部１２は、接触データ１４に含まれる時刻ｔ１の識別情報に対応する時刻ｔ２の識別情報を対応データから検索することで、上記の条件を満たす第２の粒子を検出することができる。

処理部１２は、位置データ１５に基づいて、時刻ｔ２において粒子Ｊ１と上記の第２の粒子とが接触しているか判定する。接触している場合、処理部１２は、接触データ１４から接触データ１６に第２の粒子の変位量を複写する。例えば、処理部１２は、時刻ｔ１において粒子Ｊ１と接触しかつ領域１３に属し、時刻ｔ２において領域１３に属する粒子Ｊ２を検出する。時刻ｔ２において粒子Ｊ１と粒子Ｊ２とは接触している。処理部１２は、接触データ１４から接触データ１６に、粒子Ｊ２に対応する変位量Ｄ２を複写する。

次に、処理部１２は、位置データ１５，１５－１に基づいて、上記の第２の粒子以外の領域１３，１３－１に属する粒子であって、時刻ｔ２において粒子Ｊ１と接触している第３の粒子を検出する。処理部１２は、接触データ１４から第３の粒子を検索する。接触データ１４に上記の第３の粒子が登録されている場合、処理部１２は、接触データ１４から接触データ１６に第３の粒子の変位量を複写する。

例えば、処理部１２は、時刻ｔ２において粒子Ｊ１と接触している粒子Ｊ３－１，Ｊ３－２，Ｊ３－３を検出する。粒子Ｊ３－１は、時刻ｔ１において領域１３に属し、時刻ｔ２において領域１３－１に移動した粒子である。粒子Ｊ３－２は、時刻ｔ１において領域１３－１に属し、時刻ｔ２において領域１３に移動した粒子である。粒子Ｊ３－３は、時刻ｔ１と時刻ｔ２の両方において領域１３－１に属する粒子である。処理部１２は、接触データ１４から粒子Ｊ３－１，Ｊ３－２，Ｊ３－３を検索する。粒子Ｊ３－１，Ｊ３－２，Ｊ３－３は、接触データ１４に登録されている。そこで、処理部１２は、接触データ１４から接触データ１６に、粒子Ｊ３－１，Ｊ３－２，Ｊ３－３に対応する変位量Ｄ３－１，Ｄ３－２，Ｄ３－３を複写する。なお、接触データ１４，１６では、領域１３に属する接触粒子と領域１３－１に属する接触粒子が区別して登録されてもよい。

第１の実施の形態の情報処理システムによれば、粒子シミュレーションにおいて、時刻ｔ２で粒子Ｊ１と接触している接触粒子とその変位量を示す接触データ１６が生成される。よって、時刻ｔ２に粒子Ｊ１が接触粒子から受ける接触力を効率的に算出することができ、時刻ｔ２より後の時刻ｔ３の粒子Ｊ１の位置を効率的に算出することができる。

また、接触データ１６を生成する際に、時刻ｔ１の接触データ１４を用いて、時刻ｔ１で粒子Ｊ１と接触しかつ領域１３に属し、時刻ｔ２で領域１３に留まっている第２の粒子が検出される。第２の粒子については、時刻ｔ２において粒子Ｊ１と接触していることを確認した上で、時刻ｔ１の接触データ１４から時刻ｔ２の接触データ１６に変位量が複写される。これにより、領域１３に留まっている粒子については接触データ１４を繰り返し検索しなくてよく、比較処理の計算量を削減できる。よって、接触データ１４から接触データ１６への変位量の複写を効率化できる。

更に、上記の第２の粒子以外に時刻ｔ２において粒子Ｊ１と接触している第３の粒子が検出される。第３の粒子には、領域１３から領域１３－１に移動した粒子、領域１３－１から領域１３に移動した粒子、および、領域１３－１に留まっている粒子が含まれる。第３の粒子については、接触データ１４が検索され、接触データ１４に登録されている場合、接触データ１４から接触データ１６に変位量が複写される。

このように、領域１３に留まる接触粒子とそれ以外の接触粒子とで、異なる複写手順が採用される。これは、情報処理装置１０の対象領域である領域１３に属する粒子については移動を効率的に追跡できるのに対し、隣接領域である領域１３－１に属する粒子については移動を効率的に追跡することが難しいためである。第１の実施の形態では、領域１３に留まる接触粒子とそれ以外の接触粒子とを分けることで、少なくとも前者について複写手順が効率化される。その結果、複写処理の計算量の削減が可能となる。

［第２の実施の形態］
次に、第２の実施の形態を説明する。
図２は、第２の実施の形態の情報処理システムの例を示す図である。

第２の実施の形態の情報処理システムは、離散要素法を用いて多数の粒子の動きをシミュレートする並列処理システムである。離散要素法では、粒子間の接触を考慮して各粒子の位置および速度の変化を追跡する。離散要素法では、各粒子が初期位置に配置されたステップ１から開始して、シミュレーション上の時刻を微細時間だけ進めることでステップを１つずつ進め、各粒子の位置および速度を時間軸に沿って更新する。離散要素法で取り扱う粒子は、二次元の円でもよいし三次元の球でもよい。ただし、第２の実施の形態では説明を簡単にするため、粒子が二次元の円である場合を想定する。

第２の実施の形態の情報処理システムは、計算ノード１００，１００－１，１００－２，１００－３などの複数の計算ノードおよび管理ノード２００を含む。複数の計算ノードおよび管理ノード２００は、ネットワーク３０に接続されている。ネットワーク３０は、例えば、ＬＡＮ（Local Area Network）などのデータ通信ネットワークである。なお、計算ノード１００は、第１の実施の形態の情報処理装置１０に対応する。計算ノード１００－１は、第１の実施の形態の情報処理装置１０－１に対応する。

複数の計算ノードは、各ステップにおける多数の粒子の位置および速度を分散して計算するサーバコンピュータである。粒子が配置され得るシミュレーション対象の全体領域が複数の分割領域に細分化され、１つの分割領域が１つのプロセスに割り当てられる。離散要素法に使用される各計算ノードは、１つの分割領域が割り当てられた１つのプロセスを実行する。よって、複数の計算ノードを用いて、複数の分割領域に対応する複数のプロセスが並列に実行される。第２の実施の形態では説明を簡単にするため、全体領域を４つの分割領域に細分化し、計算ノード１００，１００－１，１００－２，１００－３の４つの計算ノードを用いて４つのプロセスを並列に実行する場合を想定する。ただし、分割領域の個数を増やして、使用する計算ノードを増やすことも可能である。

計算ノード１００，１００－１，１００－２，１００－３はそれぞれ、担当する分割領域の中に位置する粒子について、ステップを１つ進めて位置および速度を更新する。その結果、１つ前のステップでは担当する分割領域に位置していた粒子が、担当する分割領域の外に移動してしまうことがある。そこで、計算ノード１００，１００－１，１００－２，１００－３は、ステップ毎に通信を行い、分割領域の境界を跨がって移動した粒子の粒子データを交換する。また、担当する分割領域の粒子のうち、隣接する分割領域との境界付近に位置する粒子は、隣接する分割領域の粒子と接触している可能性がある。そこで、計算ノード１００，１００－１，１００－２，１００－３は、ステップ毎に通信を行い、境界付近の粒子の粒子データを複写して相互に通知する。なお、第２の実施の形態では、計算ノード１００に着目して離散要素法の処理を説明する。

管理ノード２００は、複数の計算ノードを用いた離散要素法の並列処理を制御する。具体的には、管理ノード２００は、各粒子の初期位置などの入力データをユーザから受け付ける。管理ノード２００は、シミュレーション対象の全体領域を複数の分割領域に細分化し、各粒子の初期位置に応じて粒子データを複数の分割領域に分配し、複数の分割領域それぞれの入力データを生成する。管理ノード２００は、計算ノード１００，１００－１，１００－２，１００－３それぞれに、複数の分割領域のうちの１つの分割領域を割り当て、割り当てた分割領域の入力データを送信する。

全てのステップの実行が終了すると、管理ノード２００は、計算ノード１００，１００－１，１００－２，１００－３から結果データを収集する。結果データは、各ステップにおける粒子の位置を含む。計算ノード１００，１００－１，１００－２，１００－３はそれぞれ、担当する分割領域に粒子が存在する間における当該粒子の位置を含む結果データを管理ノード２００に送信することになる。管理ノード２００は、結果データをＨＤＤなどの不揮発性ストレージに保存してもよい。また、管理ノード２００は、結果データを表示装置に表示してもよい。例えば、管理ノード２００は、各粒子を点で表現し、各粒子の位置変化を点の移動で表現した動画像を表示してもよい。また、例えば、管理ノード２００は、結果データを他の情報処理装置に送信してもよい。

なお、計算ノード１００，１００－１，１００－２，１００－３および管理ノード２００を、独立のコンピュータではなく、多数のプロセッサを有する並列処理装置の中の異なるプロセッサとしてもよい。また、管理ノード２００の処理を、計算ノード１００，１００－１，１００－２，１００－３のうちの１つの計算ノードに実行させてもよい。

図３は、計算ノードのハードウェア例を示すブロック図である。
計算ノード１００は、ＣＰＵ１０１、ＲＡＭ１０２、ＨＤＤ１０３、画像インタフェース１０４、入力インタフェース１０５、媒体リーダ１０６および通信インタフェース１０７を有する。上記ユニットはバスに接続されている。ＣＰＵ１０１は、第１の実施の形態の処理部１２に対応する。ＲＡＭ１０２またはＨＤＤ１０３は、第１の実施の形態の記憶部１１に対応する。他の計算ノードや管理ノード２００も同様のハードウェアを有する。

ＣＰＵ１０１は、プログラムの命令を実行するプロセッサである。ＣＰＵ１０１は、ＨＤＤ１０３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ１０２にロードし、プログラムを実行する。なお、ＣＰＵ１０１は複数のプロセッサコアを備えてもよく、計算ノード１００は複数のプロセッサを備えてもよい。複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」と言うことがある。

ＲＡＭ１０２は、ＣＰＵ１０１が実行するプログラムやＣＰＵ１０１が演算に使用するデータを一時的に記憶する揮発性の半導体メモリである。なお、計算ノード１００は、ＲＡＭ以外の種類のメモリを備えてもよく、複数のメモリを備えてもよい。

ＨＤＤ１０３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、および、データを記憶する不揮発性ストレージである。なお、計算ノード１００は、フラッシュメモリやＳＳＤ（Solid State Drive）など他の種類のストレージを備えてもよく、複数のストレージを備えてもよい。

画像インタフェース１０４は、ＣＰＵ１０１からの命令に従って、計算ノード１００に接続された表示装置１１１に画像を出力する。表示装置１１１として、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイ、プロジェクタなど、任意の種類の表示装置を使用することができる。また、計算ノード１００に、プリンタなど表示装置１１１以外の出力デバイスが接続されてもよい。

入力インタフェース１０５は、計算ノード１００に接続された入力デバイス１１２から入力信号を受け付ける。入力デバイス１１２として、マウス、タッチパネル、タッチパッド、キーボードなど、任意の種類の入力デバイスを使用することができる。また、計算ノード１００に複数種類の入力デバイスが接続されてもよい。

媒体リーダ１０６は、記録媒体１１３に記録されたプログラムやデータを読み取る読み取り装置である。記録媒体１１３として、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤなどの磁気ディスク、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）などの光ディスク、半導体メモリなど、任意の種類の記録媒体を使用することができる。媒体リーダ１０６は、例えば、記録媒体１１３から読み取ったプログラムやデータを、ＲＡＭ１０２やＨＤＤ１０３などの他の記録媒体に複写する。読み取られたプログラムは、例えば、ＣＰＵ１０１によって実行される。なお、記録媒体１１３は可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体１１３やＨＤＤ１０３を、コンピュータ読み取り可能な記録媒体と言うことがある。

通信インタフェース１０７は、ネットワーク３０に接続され、ネットワーク３０を介して他の計算ノードや管理ノード２００と通信する。通信インタフェース１０７は、例えば、スイッチやルータなどの有線通信装置に接続される有線通信インタフェースである。

次に、離散要素法について説明する。
図４は、離散要素法における粒子モデルの例を示す図である。
離散要素法では、着目する粒子にかかる合力を算出し、着目する粒子の質量と合力から加速度を算出し、現在の速度と加速度から次の速度を算出し、現在の位置と次の速度から次の位置を算出する。着目する粒子にかかる合力は、重力などの外力と、着目する粒子と接触している１以上の他の粒子（接触粒子）から受ける接触力とを合成したものである。そのため、離散要素法では、粒子間の接触力を表現するモデルが使用される。

ここでは、粒子４１と粒子４２とが接触している場合を考える。粒子４１の位置は、粒子４１の中心の座標で表現される。粒子４２の位置は、粒子４２の中心の座標で表現される。粒子４１の半径と粒子４２の半径は同じでもよいし異なってもよい。シミュレーションを簡単にするため、全ての粒子の半径が同じであると仮定することもできる。粒子４１と粒子４２との間では、強く衝突した場合に食い込みが発生することがある。そこで、粒子４１の中心と粒子４２の中心との間の距離が、粒子４１の半径と粒子４２の半径の合計以下である場合、粒子４１と粒子４２が接触していると判定される。

粒子４１が粒子４２と接触している場合、粒子４１が粒子４２から受ける接触力は法線方向と接線方向に分解される。法線方向の接触力は、並列に配置されたダッシュポット４３とばね４５によりモデル化される。接線方向の接触力は、並列に配置されたダッシュポット４４とばね４６によりモデル化される。ダッシュポット４３，４４は、速度に比例する粘性抵抗力を、進行方向の逆方向に与える部品である。速度が大きいほど粘性抵抗力が大きくなり、速度が小さいほど粘性抵抗力が小さくなる。ばね４５，４６は、当初の長さからの変位量に比例する弾性抵抗力を、変位方向の逆方向に与える部品である。変位量が大きいほど弾性抵抗力が大きくなり、変位量が小さいほど弾性抵抗力が小さくなる。

法線方向の接触力は、ダッシュポット４３の粘性抵抗力とばね４５の弾性抵抗力の和である。よって、ステップｔ（ｔは１以上の整数）における法線方向の接触力Ｆ_Ｎ（ｔ）は数式（１）のように定義される。ｋ_Ｎはばね４５の弾性係数であり、ｃ_Ｎはダッシュポット４３の粘性減衰係数である。弾性係数ｋ_Ｎと粘性減衰係数ｃ_Ｎは入力データとして与えられる。弾性係数ｋ_Ｎと粘性減衰係数ｃ_Ｎは、全ての粒子について同じでもよいし粒子によって異なってもよい。Ｄ_Ｎ（ｔ）はステップｔの法線方向の変位量である。ΔＤ_Ｎ（ｔ）は、ステップｔの粒子４２に対する粒子４１の相対速度のうち法線方向の成分である。

変位量Ｄ_Ｎ（ｔ）は、粒子４１と粒子４２の間の食い込み量に相当する。よって、変位量Ｄ_Ｎ（ｔ）は、粒子４１の半径と粒子４２の半径の和から、粒子４１の中心と粒子４２の中心の間の距離を差し引くことで算出することができる。相対速度ΔＤ_Ｎ（ｔ）は、粒子４１の速度と粒子４２の速度から算出することができる。

接線方向の接触力は、ダッシュポット４４の粘性抵抗力とばね４６の弾性抵抗力の和である。よって、ステップｔにおける接線方向の接触力Ｆ_Ｔ（ｔ）は数式（２）のように定義される。ｋ_Ｔはばね４６の弾性係数であり、ｃ_Ｔはダッシュポット４４の粘性減衰係数である。弾性係数ｋ_Ｔと粘性減衰係数ｃ_Ｔは入力データとして与えられる。弾性係数ｋ_Ｔと粘性減衰係数ｃ_Ｔは、全ての粒子について同じでもよいし粒子によって異なってもよい。Ｄ_Ｔ（ｔ）はステップｔの接線方向の変位量である。ΔＤ_Ｔ（ｔ）は、ステップｔの粒子４２に対する粒子４１の相対速度のうち接線方向の成分である。

変位量Ｄ_Ｔ（ｔ）は、粒子４１と粒子４２が接触を開始してから接触状態のまま接線方向にずれた累積の移動量である。粒子４２に対する粒子４１の相対速度はステップ毎に変化し得るため、変位量Ｄ_Ｔ（ｔ）は、粒子４１と粒子４２が接触を開始してから離れるまでの各ステップにおいて単位ステップ時間Δｔ（１ステップに相当する時間幅）の間の接線方向の移動量を算出し、それら接線方向の移動量を累積することで算出される。相対速度ΔＤ_Ｔ（ｔ）は、粒子４１の速度と粒子４２の速度から算出することができる。

以上の法線方向の接触力と接線方向の接触力を合成することで、粒子４１が粒子４２から受ける接触力を算出することができる。粒子４１が更に他の粒子と接触している場合、粒子４１が他の粒子から受ける接触力も算出される。粒子４１が１以上の接触粒子から受ける接触力と重力などの外力とを合成することで、粒子４１にかかる合力が算出される。

前述のように、法線方向の変位量Ｄ_Ｎ（ｔ）は、ステップｔにおける粒子４１，４２の位置から算出することができる。法線方向の相対速度ΔＤ_Ｎ（ｔ）は、ステップｔにおける粒子４１，４２の速度から算出することができる。接線方向の相対速度ΔＤ_Ｔ（ｔ）は、ステップｔにおける粒子４１，４２の速度から算出することができる。よって、変位量Ｄ_Ｎ（ｔ）および相対速度ΔＤ_Ｎ（ｔ），ΔＤ_Ｔ（ｔ）は、ステップｔの粒子４１，４２の状態から算出され、過去のステップの粒子４１，４２の状態に依存しない。

これに対して、接線方向の変位量Ｄ_Ｔ（ｔ）は、ステップｔ－１における変位量Ｄ_Ｔ（ｔ－１）に、ステップｔ－１からステップｔの間の接線方向の移動量を加算することで算出される。よって、変位量Ｄ_Ｔ（ｔ）は、ステップｔにおける粒子４１，４２の状態に加えて、過去のステップの粒子４１，４２の状態に依存する。そこで、粒子４１が粒子４２と接触している間、接線方向の変位量を状態として保持することになる。このように、離散要素法では、着目する粒子について、着目する粒子と接触している粒子毎の接線方向の変位量を前ステップから現ステップに複写して使用することになる。

図５は、シミュレーション対象領域の分割例を示す図である。
シミュレーション対象の全体領域は、例えば、分割領域５１，５２，５３，５４に分割される。分割領域５１，５２，５３，５４はそれぞれ正方形の領域である。分割領域５１は全体領域の左上に相当し、分割領域５２は全体領域の右上に相当し、分割領域５３は全体領域の左下に相当し、分割領域５４は全体領域の右下に相当する。よって、分割領域５１は、右辺で分割領域５２と接し、下辺で分割領域５３と接し、右下点で分割領域５４と接する。分割領域５２は、左辺で分割領域５１と接し、左下点で分割領域５３と接し、下辺で分割領域５４と接する。分割領域５３は、上辺で分割領域５１と接し、右上点で分割領域５２と接し、右辺で分割領域５４と接する。分割領域５４は、左上点で分割領域５１と接し、上辺で分割領域５２と接し、左辺で分割領域５３と接する。

分割領域５１は計算ノード１００に割り当てられる。分割領域５２は計算ノード１００－１に割り当てられる。分割領域５３は計算ノード１００－２に割り当てられる。分割領域５４は計算ノード１００－３に割り当てられる。

分割領域５１は、他の分割領域と接する境界から所定の距離以内にある点の集合である境界領域５１ａ，５１ｂ，５１ｃを含む。境界領域５１ｃは、下辺から所定の距離以内にありかつ右辺から所定の距離以内にある点の集合である。境界領域５１ａは、下辺から所定の距離以内にある点の集合のうち境界領域５１ｃを除く部分である。境界領域５１ｂは、右辺から所定の距離以内にある点の集合のうち境界領域５１ｃを除く部分である。

分割領域５２は、他の分割領域と接する境界から所定の距離以内にある点の集合である境界領域５２ａ，５２ｂ，５２ｃを含む。境界領域５２ｃは、下辺から所定の距離以内にありかつ左辺から所定の距離以内にある点の集合である。境界領域５２ａは、下辺から所定の距離以内にある点の集合のうち境界領域５２ｃを除く部分である。境界領域５２ｂは、左辺から所定の距離以内にある点の集合のうち境界領域５２ｃを除く部分である。

分割領域５３は、他の分割領域と接する境界から所定の距離以内にある点の集合である境界領域５３ａ，５３ｂ，５３ｃを含む。境界領域５３ｃは、上辺から所定の距離以内にありかつ右辺から所定の距離以内にある点の集合である。境界領域５３ａは、上辺から所定の距離以内にある点の集合のうち境界領域５３ｃを除く部分である。境界領域５３ｂは、右辺から所定の距離以内にある点の集合のうち境界領域５３ｃを除く部分である。

分割領域５４は、他の分割領域と接する境界から所定の距離以内にある点の集合である境界領域５４ａ，５４ｂ，５４ｃを含む。境界領域５４ｃは、上辺から所定の距離以内にありかつ左辺から所定の距離以内にある点の集合である。境界領域５４ａは、上辺から所定の距離以内にある点の集合のうち境界領域５４ｃを除く部分である。境界領域５４ｂは、左辺から所定の距離以内にある点の集合のうち境界領域５４ｃを除く部分である。

前述のように、２つの分割領域の境界を跨がって、一方の分割領域に属する粒子が他方の分割領域に属する粒子と接触する可能性がある。そこで、計算ノード１００，１００－１，１００－２，１００－３は、境界領域に属する粒子の粒子データを相互に通知する。

具体的には、計算ノード１００は、境界領域５１ｂ，５１ｃの粒子データを計算ノード１００－１に送信し、境界領域５１ａ，５１ｃの粒子データを計算ノード１００－２に送信し、境界領域５１ｃの粒子データを計算ノード１００－３に送信する。計算ノード１００－１は、境界領域５２ｂ，５２ｃの粒子データを計算ノード１００に送信し、境界領域５２ｃの粒子データを計算ノード１００－２に送信し、境界領域５２ａ，５２ｃの粒子データを計算ノード１００－３に送信する。計算ノード１００－２は、境界領域５３ａ，５３ｃの粒子データを計算ノード１００に送信し、境界領域５３ｃの粒子データを計算ノード１００－１に送信し、境界領域５３ｂ，５３ｃの粒子データを計算ノード１００－３に送信する。計算ノード１００－３は、境界領域５４ｃの粒子データを計算ノード１００に送信し、境界領域５４ａ，５４ｃの粒子データを計算ノード１００－１に送信し、境界領域５４ｂ，５４ｃの粒子データを計算ノード１００－２に送信する。

これにより、分割領域５１を担当する計算ノード１００は、分割領域５１に属する粒子の粒子データに加えて、境界領域５２ｂ，５２ｃ，５３ａ，５３ｃ，５４ｃに属する粒子の粒子データをもつことになる。計算ノード１００は、分割領域５１に属する粒子それぞれについて、接触粒子を判定して位置および速度を計算する。計算ノード１００から見て、分割領域５１を「ｍａｉｎ領域」と言うことがあり、境界領域５２ｂ，５２ｃ，５３ａ，５３ｃ，５４ｃの集合を「ｓｋｉｒｔ領域」と言うことがある。

同様に、分割領域５２を担当する計算ノード１００－１は、分割領域５２に属する粒子の粒子データに加えて、境界領域５１ｂ，５１ｃ，５３ｃ，５４ａ，５４ｃに属する粒子の粒子データをもつことになる。計算ノード１００－１は、分割領域５２に属する粒子それぞれについて、接触粒子を判定して位置および速度を計算する。計算ノード１００－１から見て、分割領域５２を「ｍａｉｎ領域」と言うことがあり、境界領域５１ｂ，５１ｃ，５３ｃ，５４ａ，５４ｃの集合を「ｓｋｉｒｔ領域」と言うことがある。

分割領域５３を担当する計算ノード１００－２は、分割領域５３に属する粒子の粒子データに加えて、境界領域５１ａ，５１ｃ，５２ｃ，５４ｂ，５４ｃに属する粒子の粒子データをもつことになる。計算ノード１００－２は、分割領域５３に属する粒子それぞれについて、接触粒子を判定して位置および速度を計算する。計算ノード１００－２から見て、分割領域５３を「ｍａｉｎ領域」と言うことがあり、境界領域５１ａ，５１ｃ，５２ｃ，５４ｂ，５４ｃの集合を「ｓｋｉｒｔ領域」と言うことがある。

分割領域５４を担当する計算ノード１００－３は、分割領域５４に属する粒子の粒子データに加えて、境界領域５１ｃ，５２ａ，５２ｃ，５３ｂ，５３ｃに属する粒子の粒子データをもつことになる。計算ノード１００－３は、分割領域５４に属する粒子それぞれについて、接触粒子を判定して位置および速度を計算する。計算ノード１００－３から見て、分割領域５４を「ｍａｉｎ領域」と言うことがあり、境界領域５１ｃ，５２ａ，５２ｃ，５３ｂ，５３ｃの集合を「ｓｋｉｒｔ領域」と言うことがある。

また、前述のように、ステップを１つ進めると、ある分割領域に属する粒子が他の分割領域に移動することがある。粒子が属する分割領域は、その粒子の中心の座標によって判断される。そこで、計算ノード１００，１００－１，１００－２，１００－３は、担当する分割領域の外に移動した粒子の粒子データを相互に交換する。

具体的には、計算ノード１００は、分割領域５１の右辺から出た粒子の粒子データを計算ノード１００－１に送信し、下辺から出た粒子の粒子データを計算ノード１００－２に送信し、右下から出た粒子の粒子データを計算ノード１００－３に送信する。計算ノード１００－１は、分割領域５２の左辺から出た粒子の粒子データを計算ノード１００に送信し、左下から出た粒子の粒子データを計算ノード１００－２に送信し、下辺から出た粒子の粒子データを計算ノード１００－３に送信する。計算ノード１００－２は、分割領域５３の上辺から出た粒子の粒子データを計算ノード１００に送信し、右上から出た粒子の粒子データを計算ノード１００－１に送信し、右辺から出た粒子の粒子データを計算ノード１００－３に送信する。計算ノード１００－３は、分割領域５４の左上から出た粒子の粒子データを計算ノード１００に送信し、上辺から出た粒子の粒子データを計算ノード１００－１に送信し、左辺から出た粒子の粒子データを計算ノード１００－２に送信する。

図６は、周辺粒子の移動例を示す第１の図である。
ｍａｉｎ領域６１は、図５の分割領域５１に相当する。ｓｋｉｒｔ領域６２は、図５の境界領域５２ｂ，５２ｃ，５３ａ，５３ｃ，５４ｃの集合に相当する。ｍａｉｎ領域６１に属する１つの粒子に着目すると、着目する粒子の周辺に存在する周辺粒子は、ｍａｉｎ領域６１に属するかｓｋｉｒｔ領域６２に属するか、および、着目する粒子と接触しているか否かによって、４通りの状態に分けられる。そして、前ステップの状態と現ステップの状態の組み合わせによって、周辺粒子の移動は１６パターンに分けられる。ここでは、１６パターンのうちの４パターンについて図示している。

粒子６３は、ｍａｉｎ領域６１の中の着目する粒子である。粒子６３－１～６３－４は、粒子６３の周辺に存在する周辺粒子である。前ステップにおいて、粒子６３－１，６３－３は、ｍａｉｎ領域６１に属しており粒子６３と接触している。粒子６３－２，６３－４は、ｓｋｉｒｔ領域６２に属しており粒子６３と接触している。現ステップでは、粒子６３－１は、ｍａｉｎ領域６１に属しており粒子６３と接触している。粒子６３－２は、ｓｋｉｒｔ領域６２からｍａｉｎ領域６１に移動しており粒子６３と接触している。粒子６３－３は、ｍａｉｎ領域６１からｓｋｉｒｔ領域６２に移動しており粒子６３と接触している。粒子６３－４は、ｓｋｉｒｔ領域６２に属しており粒子６３と接触している。

このような粒子６３－１～６３－４は、前ステップおよび現ステップの両方で粒子６３と接触している。よって、粒子６３の合力を計算するにあたり、粒子６３－１～６３－４それぞれの接線方向の変位量を前ステップから現ステップに複写することになる。

図７は、周辺粒子の移動例を示す第２の図である。
ここでは、前述の１６パターンのうちの残りの１２パターンについて図示している。粒子６３－５～６３－１６は、粒子６３の周辺に存在する周辺粒子である。

前ステップにおいて、粒子６３－５，６３－１１は、ｍａｉｎ領域６１に属しており粒子６３と接触している。粒子６３－６，６３－７，６３－１２，６３－１３は、ｍａｉｎ領域６１に属しており粒子６３と接触していない。粒子６３－８，６３－１４は、ｓｋｉｒｔ領域６２に属しており粒子６３と接触している。粒子６３－９，６３－１０，６３－１５，６３－１６は、ｓｋｉｒｔ領域６２に属しており粒子６３と接触していない。

現ステップでは、粒子６３－５は、ｍａｉｎ領域６１に属しており粒子６３から離れている。粒子６３－６は、ｍａｉｎ領域６１に属しており粒子６３と接触している。粒子６３－７は、ｍａｉｎ領域６１に属しており粒子６３と接触していない。粒子６３－８は、ｓｋｉｒｔ領域６２からｍａｉｎ領域６１に移動しており粒子６３から離れている。

粒子６３－９は、ｓｋｉｒｔ領域６２からｍａｉｎ領域６１に移動しており粒子６３と接触している。粒子６３－１０は、ｓｋｉｒｔ領域６２からｍａｉｎ領域６１に移動しており粒子６３と接触していない。粒子６３－１１は、ｍａｉｎ領域６１からｓｋｉｒｔ領域６２に移動しており粒子６３から離れている。粒子６３－１２は、ｍａｉｎ領域６１からｓｋｉｒｔ領域６２に移動しており粒子６３と接触している。

粒子６３－１３は、ｍａｉｎ領域６１からｓｋｉｒｔ領域６２に移動しており粒子６３と接触していない。粒子６３－１４は、ｓｋｉｒｔ領域６２に属しており粒子６３から離れている。粒子６３－１５は、ｓｋｉｒｔ領域６２に属しており粒子６３と接触している。粒子６３－１６は、ｓｋｉｒｔ領域６２に属しており粒子６３と接触していない。

このような粒子６３－５～６３－１６は、前ステップおよび現ステップの少なくとも一方で粒子６３と接触していない。よって、粒子６３の合力を計算するにあたり、粒子６３－５～６３－１６の接線方向の変位量の複写は発生しない。

図８は、接線方向の変位量の複写例を示す図である。
ここでは、接線方向の変位量を複写する比較的単純な方法を考える。
ｍａｉｎ領域内の１つの粒子に対して、前ステップの接触粒子リスト７１と前ステップの変位量リスト７２が生成されている。接触粒子リスト７１は、前ステップにおいて当該粒子と接触している接触粒子の粒子番号を列挙したデータである。接触粒子は、ｍａｉｎ領域に属していることもあるしｓｋｉｒｔ領域に属していることもある。変位量リスト７２は、接触粒子リスト７１に登録された粒子番号と対応付けて、粒子番号が示す接触粒子の前ステップまでの接線方向の変位量を列挙したデータである。

ステップを１つ進める、すなわち、シミュレーション上の時刻を単位ステップ時間だけ進めると、上記と同じ粒子に対して、現ステップの接触粒子リスト７３と現ステップの変位量リスト７４が生成される。接触粒子リスト７３は、現ステップにおいて当該粒子と接触している接触粒子の粒子番号を列挙したデータである。変位量リスト７４は、接触粒子リスト７３に登録された粒子番号と対応付けて、粒子番号が示す接触粒子の現ステップまでの接線方向の変位量を列挙したデータである。

現ステップにおいて接触粒子リスト７３および変位量リスト７４の生成が完了すると、前ステップの接触粒子リスト７１および変位量リスト７２は削除してよい。ただし、接触粒子リスト７１と接触粒子リスト７３に同じ接触粒子が登録されている場合、その接触粒子に対応する変位量が変位量リスト７２から変位量リスト７４に複写される。これにより、前ステップまでの変位量が現ステップの開始時点で複写されることになる。

そして、複写された変位量に、前ステップから現ステップまでの間の接線方向の相対的な移動量（変位量の増分）を加えて、現ステップまでの変位量が算出される。接線方向の変位量の増分は、例えば、着目する粒子の現ステップの速度の接線方向の成分に単位ステップ時間を乗じたものから、接触粒子の現ステップの速度の接線方向の成分に単位ステップ時間を乗じたものを引くことで算出することができる。

接線方向の変位量を複写する１つの方法として、以下の方法が考えられる。まず、前ステップの接触粒子リスト７１とは関係なく、着目する粒子とｍａｉｎ領域またはｓｋｉｒｔ領域に属する各粒子との間で接触の有無を判定し、現ステップの接触粒子を検出する。これにより、現ステップの接触粒子リスト７３が生成される。

接触粒子リスト７３が生成されると、接触粒子リスト７３に登録された接触粒子毎に、接触粒子リスト７１を走査して同じ接触粒子を接触粒子リスト７１から検索する。同じ接触粒子が接触粒子リスト７１から検出されると、接触粒子リスト７１の当該接触粒子の位置に対応する変位量リスト７２の位置から変位量を読み出す。そして、接触粒子リスト７３の当該接触粒子の位置に対応する変位量リスト７４の位置に変位量を複写する。一方、同じ接触粒子が接触粒子リスト７１に存在しない場合、接触粒子リスト７３の当該接触粒子の位置に対応する変位量リスト７４の位置に変位量＝０を書き込む。

例えば、接触粒子リスト７３の１番目の要素が接触粒子リスト７１の１番目の要素と同じ接触粒子を示している場合、変位量リスト７２の１番目の要素が変位量リスト７４の１番目の要素として複写される。また、接触粒子リスト７３の２番目の要素が接触粒子リスト７１の３番目の要素と同じ接触粒子を示している場合、変位量リスト７２の３番目の要素が変位量リスト７４の２番目の要素として複写される。また、接触粒子リスト７３の３番目および４番目の要素と同じ接触粒子を示す要素が接触粒子リスト７１に存在しない場合、変位量リスト７４の３番目および４番目の要素が変位量＝０に設定される。

上記の変位量の複写方法では、接触粒子リスト７１，７３に登録され得る接触粒子の平均個数をＮとすると、接触粒子リスト７１，７３の間で粒子番号を比較する比較処理の計算量がＯ（Ｎ^２）になる。ステップを１つ進める毎にｍａｉｎ領域に属する各粒子に対してこのような比較処理を行うことは負荷が高い。そこで、第２の実施の形態の計算ノード１００，１００－１，１００－２，１００－３は、比較処理の計算量が削減されるように、現ステップの接触粒子リスト７３を生成するアルゴリズムを改善する。現ステップの接触粒子リスト７３の生成方法の詳細は後述する。

なお、接触粒子リスト７１，７３に含まれる粒子番号は、１つの計算ノードの中で一意なローカルの識別番号である。例えば、計算ノード１００は、計算ノード１００のｍａｉｎ領域に属する粒子およびｓｋｉｒｔ領域に属する粒子に対して、計算ノード１００－１，１００－２，１００－３とは独立に粒子番号を付与する。分割領域間で粒子の出入りがあると、同じ粒子に対して前ステップと現ステップで異なる粒子番号が付与されることがある。ただし、計算ノード１００は、ｍａｉｎ領域に留まっている粒子については、前ステップの粒子番号と現ステップの粒子番号の対応を把握している。また、第２の実施の形態では、粒子の識別番号として、粒子番号に加えてグローバルＩＤを併用する。グローバルＩＤは、全ての計算ノードを通じて一意なグローバルの識別番号である。

次に、計算ノード１００の機能について説明する。
図９は、計算ノードの機能例を示すブロック図である。
計算ノード１００は、モデル記憶部１２１、中間データ記憶部１２２、結果データ記憶部１２３、データ入出力部１２４、ステップ実行部１２５およびプロセス間通信部１２６を有する。モデル記憶部１２１、中間データ記憶部１２２および結果データ記憶部１２３は、例えば、ＲＡＭ１０２またはＨＤＤ１０３の記憶領域を用いて実現される。データ入出力部１２４、ステップ実行部１２５およびプロセス間通信部１２６は、例えば、ＣＰＵ１０１が実行するプログラムを用いて実現される。計算ノード１００－１，１００－２，１００－３も、計算ノード１００と同様のモジュールを有する。

モデル記憶部１２１は、計算ノード１００が担当する分割領域についてのモデルを示す入力データを記憶する。入力データには、分割領域の大きさ、分割領域の座標、壁の位置などの領域データが含まれる。また、入力データには、分割領域に初期配置されている粒子の位置が含まれる。また、入力データには、重力、粒子の半径、接触モデルの弾性係数、接触モデルの粘性減衰係数などの物理パラメータが含まれる。また、入力データには、単位ステップ時間、最大ステップ数などの制御パラメータが含まれる。

中間データ記憶部１２２は、シミュレーション途中で生成される中間データを記憶する。中間データには、ｍａｉｎ領域に属する粒子の位置や速度を示す粒子データと、ｓｋｉｒｔ領域に属する粒子の位置や速度を示す粒子データが含まれる。また、中間データには、前ステップの粒子番号と現ステップの粒子番号の対応情報が含まれる。また、中間データには、ｍａｉｎ領域に属する粒子毎の接触粒子リストや変位量リストが含まれる。

結果データ記憶部１２３は、計算ノード１００が担当する分割領域のシミュレーション結果を示す結果データを記憶する。結果データには、ステップ毎に、当該ステップで分割領域の中に存在していた粒子およびその粒子の位置を示す情報が含まれる。各粒子の位置を時系列に追跡することで、粉体の動きを可視化することが可能となる。

データ入出力部１２４は、並列処理を制御する管理ノード２００と通信する。データ入出力部１２４は、計算ノード１００が担当する分割領域の入力データを管理ノード２００から受信し、入力データをモデル記憶部１２１に格納する。また、データ入出力部１２４は、ステップ数が最大ステップ数に到達すると、結果データ記憶部１２３から結果データを読み出して管理ノード２００に送信する。管理ノード２００では、複数の分割領域の結果データが集約され、全体領域における粉体の動きが可視化される。

ステップ実行部１２５は、モデル記憶部１２１に記憶された入力データを用いて、担当する分割領域の中での粒子の動きを時間軸に沿って１ステップずつ計算する。ステップを１つ進めることは、シミュレーション上の時刻を所定の微細時間だけ進めることに相当する。１つのステップは、分割領域に属する粒子毎に、接触粒子を判定して接触粒子から受ける接触力を算出し、接触力と重力を合成した合力を算出し、現在の位置と速度と合力から次の時刻の位置と速度を算出することを含む。ステップ実行部１２５は、中間データ記憶部１２２に記憶された中間データを適宜更新する。また、ステップ実行部１２５は、各ステップにおける粒子の位置を結果データ記憶部１２３に書き込む。

プロセス間通信部１２６は、計算ノード１００のステップ実行部１２５で実行されるプロセスと計算ノード１００－１，１００－２，１００－３で実行されるプロセスとの間のプロセス間通信を実現する。プロセス間通信部１２６は、中間データ記憶部１２２から一部の中間データを読み出して計算ノード１００－１，１００－２，１００－３に送信する。また、プロセス間通信部１２６は、計算ノード１００－１，１００－２，１００－３から中間データを受信して中間データ記憶部１２２に格納する。

計算ノード１００－１，１００－２，１００－３に送信する中間データには、計算ノード１００が担当する分割領域の外に出た粒子の粒子データが含まれる。また、送信する中間データには、計算ノード１００－１，１００－２，１００－３にとってｓｋｉｒｔ領域に相当する境界領域にある粒子の粒子データが含まれる。計算ノード１００－１，１００－２，１００－３から受信する中間データには、計算ノード１００が担当する分割領域の中に入った粒子の粒子データが含まれる。また、受信する中間データには、計算ノード１００にとってｓｋｉｒｔ領域に相当する境界領域にある粒子の粒子データが含まれる。

図１０は、結果テーブルの例を示す図である。
結果テーブル１３１は、結果データ記憶部１２３に記憶される。結果テーブル１３１は、ステップ番号、グローバルＩＤおよび位置の項目を含む。ステップ番号の項目には、ステップを識別する１以上の整数が登録される。ステップ番号＝１は、シミュレーション開始時の初期状態に相当する。グローバルＩＤの項目には、粒子を識別する識別番号であって、全ての粒子を通じて一意なグローバルの識別番号が登録される。位置の項目には、粒子の位置を示す座標が登録される。なお、入力データには、粒子の初期配置の情報として、ステップ番号＝１に対応するグローバルＩＤと位置が含まれている。

図１１は、粒子テーブルおよび粒子番号対応テーブルの例を示す図である。
ｍａｉｎ粒子テーブル１３２は、中間データ記憶部１２２に記憶される。ｍａｉｎ粒子テーブル１３２は、ｍａｉｎ領域に属する粒子の最新の粒子データを示す。ｍａｉｎ粒子テーブル１３２は、グローバルＩＤ、位置および速度の項目を含む。グローバルＩＤの項目には、ｍａｉｎ領域に属する粒子に付与されたグローバルＩＤが登録される。位置の項目には、ｍａｉｎ領域に属する粒子の位置を示す座標が登録される。速度の項目には、ｍａｉｎ領域に属する粒子の速度ベクトルが登録される。シミュレーション対象領域が二次元である場合、座標や速度ベクトルは二次元ベクトルである。

ｍａｉｎ粒子テーブル１３２のレコードの順番を示すインデックスが、ｍａｉｎ領域に属する粒子の粒子番号として使用される。例えば、ｍａｉｎ粒子テーブル１３２の１番目のレコードが示す粒子は、粒子番号＝１をもつ。ｍａｉｎ領域の外に粒子が出た場合、その粒子に対応するレコードがｍａｉｎ粒子テーブル１３２から削除される。ｍａｉｎ領域の中に粒子が入った場合、その粒子に対応するレコードがｍａｉｎ粒子テーブル１３２に挿入される。よって、ステップを１つ進めると既存のレコードのインデックスが変わることがあり、結果として既存の粒子の粒子番号が変わることがある。

なお、ｍａｉｎ粒子テーブル１３２を多次元配列として形成することもできる。また、ｍａｉｎ粒子テーブル１３２を、グローバルＩＤ配列と位置配列と速度配列の３つの一次元配列として形成することもできる。インデックスは配列の添え字に相当する。ｍａｉｎ粒子テーブル１３２を３つの一次元配列として形成した場合、同じインデックスが指すグローバルＩＤと位置と速度は、同じ粒子についての情報として対応付けられる。

ｓｋｉｒｔ粒子テーブル１３３は、中間データ記憶部１２２に記憶される。ｓｋｉｒｔ粒子テーブル１３３は、ｓｋｉｒｔ領域に属する粒子の最新の粒子データを示す。ｓｋｉｒｔ粒子テーブル１３３は、グローバルＩＤ、位置および速度の項目を含む。グローバルＩＤの項目には、ｓｋｉｒｔ領域に属する粒子に付与されたグローバルＩＤが登録される。位置の項目には、ｓｋｉｒｔ領域に属する粒子の位置を示す座標が登録される。速度の項目には、ｓｋｉｒｔ領域に属する粒子の速度ベクトルが登録される。

ｓｋｉｒｔ粒子テーブル１３３のレコードの順番を示すインデックスに、所定のプレフィックスを付加したものが、ｓｋｉｒｔ領域に属する粒子の粒子番号として使用される。プレフィックスは、ｍａｉｎ領域に属する粒子の粒子番号と衝突しないように付与される。例えば、ｓｋｉｒｔ粒子テーブル１３３の１番目のレコードが示す粒子は、粒子番号＝９００００１をもつ。ｓｋｉｒｔ粒子テーブル１３３の粒子データは、計算ノード１００－１，１００－２，１００－３から受信したものである。

なお、ｓｋｉｒｔ粒子テーブル１３３を多次元配列として形成することもできる。また、ｓｋｉｒｔ粒子テーブル１３３を、グローバルＩＤ配列と位置配列と速度配列の３つの一次元配列として形成することもできる。インデックスは配列の添え字に相当する。ｓｋｉｒｔ粒子テーブル１３３を３つの一次元配列として形成した場合、同じインデックスが指すグローバルＩＤと位置と速度は、同じ粒子についての情報として対応付けられる。

ｍａｉｎ粒子テーブル１３２の位置および速度は、計算ノード１００によって更新され得る。また、ｍａｉｎ粒子テーブル１３２に登録されたレコードは、次のステップでも保持され得る。これに対して、ｓｋｉｒｔ粒子テーブル１３３の位置および速度は、計算ノード１００－１，１００－２，１００－３によって算出されるものであり、計算ノード１００によって更新されない。また、ｓｋｉｒｔ粒子テーブル１３３に登録されたレコードは、ステップ毎にリセットされ次のステップでは保持されない。計算ノード１００は、ステップ毎に、計算ノード１００－１，１００－２，１００－３から受信したｓｋｉｒｔ領域の粒子データを集約して、新規のｓｋｉｒｔ粒子テーブル１３３を生成すればよい。

粒子番号対応テーブル１３４は、中間データ記憶部１２２に記憶される。粒子番号対応テーブル１３４は、前ステップにおいてｍａｉｎ領域に属しており現ステップでもｍａｉｎ領域に留まっている粒子について、前ステップの粒子番号（旧粒子番号）と現ステップの粒子番号を対応付ける。粒子番号対応テーブル１３４は、旧粒子番号の一次元配列である。粒子番号対応テーブル１３４のインデックスは、現ステップの粒子番号に対応する。

１つのインデックスと当該インデックスが指す旧粒子番号は、同じ粒子についての現ステップの粒子番号と前ステップの粒子番号を表している。例えば、粒子番号対応テーブル１３４の１番目のレコードが旧粒子番号＝３を示している場合、前ステップの粒子番号＝３の粒子と現ステップの粒子番号＝２の粒子が同一である。旧粒子番号から現ステップの粒子番号を高速に検索できるように粒子番号対応テーブル１３４を形成してもよい。

なお、計算ノード１００は、ｍａｉｎ領域の粒子の出入りをｍａｉｎ粒子テーブル１３２に反映させる際に粒子番号対応テーブル１３４を生成する。計算ノード１００は、ｍａｉｎ粒子テーブル１３２のレコード削除やレコード挿入を行うことから、ｍａｉｎ領域に留まる粒子の粒子番号の変化を把握している。一方、計算ノード１００は、ｓｋｉｒｔ粒子テーブル１３３をステップ毎にリセットすることから、ｓｋｉｒｔ領域に留まる粒子やｍａｉｎ領域とｓｋｉｒｔ領域の間で移動した粒子の粒子番号の変化は把握しない。

図１２は、接触粒子リストの例を示す図である。
第２の実施の形態では、計算ノード１００は、ｍａｉｎ領域に属する粒子それぞれに対して、接触粒子リストとグローバルＩＤリストと変位量リストを生成する。現ステップの接触粒子リストとグローバルＩＤリストと変位量リストが完成するまで、一時的に前ステップの接触粒子リストとグローバルＩＤリストと変位量リストが併存する。また、第２の実施の形態では、計算ノード１００は、相手粒子である接触粒子を、ｍａｉｎ領域に属する接触粒子とｓｋｉｒｔ領域に属する接触粒子に分類する。そして、計算ノード１００は、接触粒子リストとグローバルＩＤリストと変位量リストをそれぞれ、ｍａｉｎ領域に関するリストとｓｋｉｒｔ領域に関するリストに分割する。

具体的には、中間データ記憶部１２２は、前ステップについて、ｍａｉｎ接触粒子リスト１４１、グローバルＩＤリスト１４２、変位量リスト１４３、ｓｋｉｒｔ接触粒子リスト１４４、グローバルＩＤリスト１４５および変位量リスト１４６を記憶する。また、中間データ記憶部１２２は、現ステップについて、ｍａｉｎ接触粒子リスト１５１、グローバルＩＤリスト１５２、変位量リスト１５３、ｓｋｉｒｔ接触粒子リスト１５４、グローバルＩＤリスト１５５および変位量リスト１５６を記憶する。

ｍａｉｎ接触粒子リスト１４１は、前ステップにおいて着目する粒子と接触している接触粒子のうち、ｍａｉｎ領域に属する接触粒子の粒子番号を列挙したデータである。グローバルＩＤリスト１４２は、ｍａｉｎ接触粒子リスト１４１に登録された接触粒子のグローバルＩＤを列挙したデータである。変位量リスト１４３は、着目する粒子とｍａｉｎ接触粒子リスト１４１に登録された接触粒子との間の前ステップまでの接線方向の変位量を列挙したデータである。ｍａｉｎ接触粒子リスト１４１とグローバルＩＤリスト１４２と変位量リスト１４３の長さは同じである。ｍａｉｎ接触粒子リスト１４１のｎ番目（ｎは１以上の整数）の粒子番号と、グローバルＩＤリスト１４２のｎ番目のグローバルＩＤと、変位量リスト１４３のｎ番目の変位量は、同一の接触粒子の情報である。

ｓｋｉｒｔ接触粒子リスト１４４は、前ステップにおいて着目する粒子と接触している接触粒子のうち、ｓｋｉｒｔ領域に属する接触粒子の粒子番号を列挙したデータである。グローバルＩＤリスト１４５は、ｓｋｉｒｔ接触粒子リスト１４４に登録された接触粒子のグローバルＩＤを列挙したデータである。変位量リスト１４６は、着目する粒子とｓｋｉｒｔ接触粒子リスト１４４に登録された接触粒子との間の前ステップまでの接線方向の変位量を列挙したデータである。ｓｋｉｒｔ接触粒子リスト１４４とグローバルＩＤリスト１４５と変位量リスト１４６の長さは同じである。ｓｋｉｒｔ接触粒子リスト１４４のｎ番目の粒子番号と、グローバルＩＤリスト１４５のｎ番目のグローバルＩＤと、変位量リスト１４６のｎ番目の変位量は、同一の接触粒子の情報である。

ｍａｉｎ接触粒子リスト１５１は、現ステップにおいて着目する粒子と接触している接触粒子のうち、ｍａｉｎ領域に属する接触粒子の粒子番号を列挙したデータである。グローバルＩＤリスト１５２は、ｍａｉｎ接触粒子リスト１５１に登録された接触粒子のグローバルＩＤを列挙したデータである。変位量リスト１５３は、着目する粒子とｍａｉｎ接触粒子リスト１５１に登録された接触粒子との間の現ステップまでの接線方向の変位量を列挙したデータである。ｍａｉｎ接触粒子リスト１５１とグローバルＩＤリスト１５２と変位量リスト１５３の長さは同じである。ｍａｉｎ接触粒子リスト１５１のｎ番目の粒子番号と、グローバルＩＤリスト１５２のｎ番目のグローバルＩＤと、変位量リスト１５３のｎ番目の変位量は、同一の接触粒子の情報である。

ｓｋｉｒｔ接触粒子リスト１５４は、現ステップにおいて着目する粒子と接触している接触粒子のうち、ｓｋｉｒｔ領域に属する接触粒子の粒子番号を列挙したデータである。グローバルＩＤリスト１５５は、ｓｋｉｒｔ接触粒子リスト１５４に登録された接触粒子のグローバルＩＤを列挙したデータである。変位量リスト１５６は、着目する粒子とｓｋｉｒｔ接触粒子リスト１５４に登録された接触粒子との間の現ステップまでの接線方向の変位量を列挙したデータである。ｓｋｉｒｔ接触粒子リスト１５４とグローバルＩＤリスト１５５と変位量リスト１５６の長さは同じである。ｓｋｉｒｔ接触粒子リスト１５４のｎ番目の粒子番号と、グローバルＩＤリスト１５５のｎ番目のグローバルＩＤと、変位量リスト１５６のｎ番目の変位量は、同一の接触粒子の情報である。

なお、ｍａｉｎ領域に属する粒子のうち、ｓｋｉｒｔ領域から近い粒子がもつｓｋｉｒｔ接触粒子リスト１４４，１５４、グローバルＩＤリスト１４５，１５５および変位量リスト１４６，１５６は、空でない（要素数が１以上である）ことがある。一方、ｓｋｉｒｔ領域から遠い粒子がもつｓｋｉｒｔ接触粒子リスト１４４，１５４、グローバルＩＤリスト１４５，１５５および変位量リスト１４６，１５６は、空である（要素数が０である）と期待される。よって、多くの粒子のｓｋｉｒｔ接触粒子リスト１４４，１５４、グローバルＩＤリスト１４５，１５５および変位量リスト１４６，１５６は空になる。

次に、計算ノード１００の処理手順について説明する。
図１３は、離散要素法の手順例を示すフローチャートである。
（Ｓ１０）データ入出力部１２４は、計算ノード１００が担当する分割領域の入力データを管理ノード２００から受信する。ステップ実行部１２５は、入力データに基づいてｍａｉｎ粒子テーブル１３２を初期化する。初期化されたｍａｉｎ粒子テーブル１３２には、入力データに含まれる各粒子のグローバルＩＤおよび位置が登録される。ｍａｉｎ粒子テーブル１３２の各粒子の速度は０に設定される。ただし、入力データに速度の初期値が含まれる場合、その初期値がｍａｉｎ粒子テーブル１３２に登録される。

（Ｓ１１）プロセス間通信部１２６は、ｍａｉｎ粒子テーブル１３２から、ｍａｉｎ領域に属する粒子のうち境界領域にある粒子の粒子データを検索する。境界領域は、ｍａｉｎ領域のうちｓｋｉｒｔ領域から所定の距離以内にある領域である。粒子が境界領域に属するか否かは、入力データに含まれる分割領域の範囲の座標とｍａｉｎ粒子テーブル１３２に含まれる位置の座標から判定できる。１つの粒子の粒子データは、ｍａｉｎ粒子テーブル１３２の１つのレコードに相当し、グローバルＩＤと位置と速度を含む。プロセス間通信部１２６は、ｍａｉｎ粒子テーブル１３２から検索された境界領域の粒子データを複写して計算ノード１００－１，１００－２，１００－３に送信する。

（Ｓ１２）プロセス間通信部１２６は、計算ノード１００－１，１００－２，１００－３から境界領域の粒子データを受信する。ここで受信する粒子データは、計算ノード１００にとってのｓｋｉｒｔ領域に属する粒子の粒子データである。プロセス間通信部１２６は、受信した粒子データをｓｋｉｒｔ粒子テーブル１３３に登録する。１つの粒子の粒子データは、ｓｋｉｒｔ粒子テーブル１３３の１つのレコードに相当し、グローバルＩＤと位置と速度を含む。なお、ｓｋｉｒｔ粒子テーブル１３３に前ステップの粒子データが残っている場合、前ステップの粒子データは破棄される。

（Ｓ１３）ステップ実行部１２５は、ｍａｉｎ領域に属する粒子それぞれに対応する接触粒子リストとグローバルＩＤリストと変位量リストを生成する。ここでは、ステップ実行部１２５は、現ステップのｍａｉｎ接触粒子リスト１５１、ｓｋｉｒｔ接触粒子リスト１５４、グローバルＩＤリスト１５２，１５５および変位量リスト１５３，１５６を生成する。ステップ実行部１２５は、前ステップのｍａｉｎ接触粒子リスト１４１、ｓｋｉｒｔ接触粒子リスト１４４、グローバルＩＤリスト１４２，１４５および変位量リスト１４３，１４６を、複写処理が完了した後に破棄する。複写処理では、前ステップの変位量リスト１４３，１４６から現ステップの変位量リスト１５３，１５６に変位量が複写されることがある。接触粒子リスト生成の詳細は後述する。

（Ｓ１４）ステップ実行部１２５は、ｍａｉｎ領域に属する粒子それぞれの変位量リスト１５３，１５６に含まれる変位量を更新する。変位量リスト１５３の更新では、ステップ実行部１２５は、ｍａｉｎ粒子テーブル１３２から接触粒子の速度を読み出し、接触粒子に対する着目する粒子の相対速度を算出し、相対速度の接線方向の成分を算出する。ステップ実行部１２５は、相対速度の接線方向の成分に基づいて、前ステップから現ステップまでの接線方向の変位量の増分を算出し、変位量リスト１５３の変位量に増分を加算する。変位量リスト１５６の更新では、同様に、ステップ実行部１２５は、ｓｋｉｒｔ粒子テーブル１３３から接触粒子の速度を読み出し、前ステップから現ステップまでの接線方向の変位量の増分を算出して変位量リスト１５６の変位量に加算する。なお、変位量リスト１５３，１５６の更新は、ステップＳ１３の中で行ってもよい。

（Ｓ１５）ステップ実行部１２５は、ｍａｉｎ領域に属する粒子それぞれにかかる合力を算出する。合力の算出では、ステップ実行部１２５は、ｍａｉｎ接触粒子リスト１５１に登録された接触番号の粒子データをｍａｉｎ粒子テーブル１３２から読み出し、ｓｋｉｒｔ接触粒子リスト１５４に登録された接触番号の粒子データをｓｋｉｒｔ粒子テーブル１３３から読み出す。ステップ実行部１２５は、接触粒子それぞれについて、着目する粒子の位置および速度と、接触粒子の位置および速度と、接触粒子の接線方向の変位量と、入力データに含まれる弾性係数および粘性減衰係数から、前述の接触モデルに従って接触力を算出する。ステップ実行部１２５は、着目する粒子が１以上の接触粒子から受ける接触力と重力などの外力とを合成して、着目する粒子にかかる合力を算出する。

（Ｓ１６）ステップ実行部１２５は、ｍａｉｎ領域に属する粒子それぞれについて、ステップＳ１５で算出した合力と粒子の質量から加速度を算出し、現在の位置と現在の速度と加速度から、単位ステップ時間後の位置と速度を算出する。ステップ実行部１２５は、ｍａｉｎ粒子テーブル１３２に登録された各粒子の位置および速度を更新する。また、ステップ実行部１２５は、結果テーブル１３１に各粒子の次ステップの位置を追記する。

（Ｓ１７）プロセス間通信部１２６は、ｍａｉｎ粒子テーブル１３２から、計算ノード１００が担当する分割領域（ｍａｉｎ領域）の外に出た粒子の粒子データを抽出する。粒子が分割領域の外に出たか否かは、入力データに含まれる分割領域の範囲の座標とｍａｉｎ粒子テーブル１３２に含まれる位置の座標から判定できる。プロセス間通信部１２６は、ｍａｉｎ粒子テーブル１３２から抽出した粒子データを計算ノード１００－１，１００－２，１００－３に送信する。送信する粒子データには、グローバルＩＤと位置と速度が含まれる。抽出した粒子データはｍａｉｎ粒子テーブル１３２から削除される。なお、ｍａｉｎ粒子テーブル１３２から１以上のレコードを削除した場合、隙間が生じないように残っているレコードをシフトして詰めるようにしてもよい。これにより、残っているレコードのインデックスが変更されることがある。

（Ｓ１８）プロセス間通信部１２６は、計算ノード１００－１，１００－２，１００－３から、計算ノード１００が担当する分割領域（ｍａｉｎ領域）の中に入った粒子の粒子データを受信する。受信する粒子データには、グローバルＩＤと位置と速度が含まれる。プロセス間通信部１２６は、受信した粒子データをｍａｉｎ粒子テーブル１３２に挿入する。例えば、新たなレコードがｍａｉｎ粒子テーブル１３２の末尾に追加される。

（Ｓ１９）プロセス間通信部１２６は、ステップＳ１７，Ｓ１８におけるｍａｉｎ粒子テーブル１３２のレコードの削除および挿入を監視し、削除されずに残っている既存レコードのインデックスの変化を追跡する。既存レコードは、現ステップでｍａｉｎ領域に属しており次ステップでもｍａｉｎ領域に留まっている粒子を示す。プロセス間通信部１２６は、現ステップの粒子番号と次ステップの粒子番号の対応を示す粒子番号対応テーブル１３４を生成する。古い粒子番号対応テーブル１３４は破棄してよい。

（Ｓ２０）ステップ実行部１２５は、ステップの繰り返しを終了するか判断する。例えば、ステップ実行部１２５は、ステップ数が入力データで指定された最大ステップ数に到達したか判断する。繰り返しを終了する場合はステップＳ２１に進み、繰り返しをまだ終了しない場合はステップＳ１１に戻る。

（Ｓ２１）データ入出力部１２４は、結果データ記憶部１２３から結果テーブル１３１を読み出し、結果テーブル１３１の結果データを管理ノード２００に送信する。
図１４は、接触粒子リスト生成の手順例を示すフローチャートである。

接触粒子リスト生成は、前述のステップＳ１３で実行される。接触粒子リスト生成は、ｍａｉｎ領域に属する粒子それぞれに対して実行される。
以下のステップＳ３０～Ｓ３７は、前ステップではｍａｉｎ領域で接触しており現ステップでもｍａｉｎ領域で接触している接触粒子（ｍａｉｎ－ｍａｉｎ粒子と言うことがある）の変位量の複写を含む。この接触粒子は前述の粒子６３－１に相当する。

（Ｓ３０）ステップ実行部１２５は、前ステップのｍａｉｎ接触粒子リスト１４１に登録された旧粒子番号の中から、旧粒子番号を１つ選択する。
（Ｓ３１）ステップ実行部１２５は、粒子番号対応テーブル１３４から、ステップＳ３０で選択した旧粒子番号に対応する現ステップの粒子番号を検索する。

（Ｓ３２）ステップ実行部１２５は、ステップＳ３１において該当する粒子番号が粒子番号対応テーブル１３４から検索されたか判断する。該当する粒子番号がある場合はステップＳ３３に進み、該当する粒子番号が無い場合はステップＳ３７に進む。

（Ｓ３３）ステップ実行部１２５は、ｍａｉｎ粒子テーブル１３２から、ステップＳ３１で検索された粒子番号に対応する位置を読み出す。ステップ実行部１２５は、着目する粒子の位置と粒子番号が示す相手粒子の位置とに基づいて、接触の有無を判定する。着目する粒子と相手粒子の間の距離が、２つの粒子の半径の和以下であれば接触していると判定され、２つの粒子の半径の和より大きければ接触していないと判定される。なお、着目する粒子の位置は、ｍａｉｎ粒子テーブル１３２から予め読み出しておく。

（Ｓ３４）ステップ実行部１２５は、２つの粒子が接触しているか判断する。接触している場合はステップＳ３５に進み、接触していない場合はステップＳ３７に進む。
（Ｓ３５）ステップ実行部１２５は、現ステップのｍａｉｎ接触粒子リスト１５１の末尾に、相手粒子の粒子番号を追加する。また、ステップ実行部１２５は、現ステップのグローバルＩＤリスト１５２の末尾に、相手粒子のグローバルＩＤを追加する。グローバルＩＤは、ステップＳ３０で選択した旧粒子番号に対応するものを前ステップのグローバルＩＤリスト１４２から読み出してもよいし、現ステップの粒子番号に対応するものをｍａｉｎ粒子テーブル１３２から読み出すようにしてもよい。

（Ｓ３６）ステップ実行部１２５は、前ステップの変位量リスト１４３から、ステップＳ３０で選択した旧粒子番号に対応する変位量を読み出す。変位量リスト１４３の所望の変位量には、ｍａｉｎ接触粒子リスト１４１の旧粒子番号と同じインデックスを用いてアクセスすることが可能である。ステップ実行部１２５は、読み出した変位量を、現ステップの変位量リスト１５３の末尾に追加する。

（Ｓ３７）ステップ実行部１２５は、ステップＳ３０においてｍａｉｎ接触粒子リスト１４１の全ての旧粒子番号を選択したか判断する。全ての旧粒子番号を選択した場合はステップＳ３８に進み、未選択の旧粒子番号がある場合はステップＳ３０に戻る。

図１５は、接触粒子リスト生成の手順例を示すフローチャート（続き１）である。
以下のステップＳ３８～Ｓ４６は、前ステップではｓｋｉｒｔ領域で接触しており現ステップではｍａｉｎ領域で接触している接触粒子（ｓｋｉｒｔ－ｍａｉｎ粒子と言うことがある）の変位量の複写を含む。この接触粒子は前述の粒子６３－２に相当する。

（Ｓ３８）ステップ実行部１２５は、ｍａｉｎ粒子テーブル１３２から、着目する粒子以外の粒子番号であって、ステップＳ３０～Ｓ３７の対象となっていない粒子番号を１つ選択する。ステップＳ３０～Ｓ３７の対象となったか否かは、フラグによって管理するようにしてもよい。ここで選択する粒子番号が示す相手粒子は、ｍａｉｎ領域に属する粒子のうち接触の有無をまだ判定していないものである。

（Ｓ３９）ステップ実行部１２５は、ｍａｉｎ粒子テーブル１３２から、ステップＳ３８で選択した粒子番号に対応する位置を読み出す。ステップ実行部１２５は、着目する粒子の位置と粒子番号が示す相手粒子の位置とに基づいて、接触の有無を判定する。

（Ｓ４０）ステップ実行部１２５は、２つの粒子が接触しているか判断する。接触している場合はステップＳ４１に進み、接触していない場合はステップＳ４６に進む。
（Ｓ４１）ステップ実行部１２５は、現ステップのｍａｉｎ接触粒子リスト１５１の末尾に、相手粒子の粒子番号を追加する。また、ステップ実行部１２５は、ｍａｉｎ粒子テーブル１３２から、相手粒子の粒子番号に対応するグローバルＩＤを読み出し、現ステップのグローバルＩＤリスト１５２の末尾に、相手粒子のグローバルＩＤを追加する。

（Ｓ４２）ステップ実行部１２５は、前ステップのグローバルＩＤリスト１４５を先頭から末尾に向かって走査し、相手粒子のグローバルＩＤを検索する。グローバルＩＤリスト１４５は、前ステップでｓｋｉｒｔ領域に属する粒子を示している。

（Ｓ４３）ステップ実行部１２５は、ステップＳ４２において該当するグローバルＩＤが検出されたか判断する。該当するグローバルＩＤがある場合はステップＳ４５に進み、該当するグローバルＩＤが無い場合はステップＳ４４に進む。

（Ｓ４４）ステップ実行部１２５は、現ステップの変位量リスト１５３の末尾に、変位量が０であることを示す数値を追加する。そして、ステップＳ４６に進む。
（Ｓ４５）ステップ実行部１２５は、前ステップの変位量リスト１４６から、ステップＳ４２で検出されたグローバルＩＤに対応する変位量を読み出す。変位量リスト１４６の所望の変位量には、グローバルＩＤリスト１４５のグローバルＩＤと同じインデックスを用いてアクセスすることが可能である。ステップ実行部１２５は、読み出した変位量を、現ステップの変位量リスト１５３の末尾に追加する。

（Ｓ４６）ステップ実行部１２５は、ステップＳ３８においてｍａｉｎ粒子テーブル１３２の全ての粒子番号を選択したか判断する。全ての粒子番号を選択した場合はステップＳ４７に進み、未選択の粒子番号がある場合はステップＳ３８に戻る。

図１６は、接触粒子リスト生成の手順例を示すフローチャート（続き２）である。
以下のステップＳ４７～Ｓ５７は、前ステップではｍａｉｎ領域で接触しており現ステップではｓｋｉｒｔ領域で接触している接触粒子（ｍａｉｎ－ｓｋｉｒｔ粒子と言うことがある）の変位量の複写を含む。この接触粒子は前述の粒子６３－３に相当する。また、ステップＳ４７～Ｓ５７は、前ステップではｓｋｉｒｔ領域で接触しており現ステップでもｓｋｉｒｔ領域で接触している接触粒子（ｓｋｉｒｔ－ｓｋｉｒｔ粒子と言うことがある）の変位量の複写を含む。この接触粒子は前述の粒子６３－４に相当する。

（Ｓ４７）ステップ実行部１２５は、ｓｋｉｒｔ粒子テーブル１３３を参照して、現ステップでｓｋｉｒｔ領域に属する粒子の粒子番号を１つ選択する。粒子番号は、ｓｋｉｒｔ粒子テーブル１３３のインデックスに所定のプレフィックスを付加したものである。

（Ｓ４８）ステップ実行部１２５は、ｓｋｉｒｔ粒子テーブル１３３から、ステップＳ４７で選択した粒子番号に対応する位置を読み出す。ステップ実行部１２５は、着目する粒子の位置と粒子番号が示す相手粒子の位置とに基づいて、接触の有無を判定する。

（Ｓ４９）ステップ実行部１２５は、２つの粒子が接触しているか判断する。接触している場合はステップＳ５０に進み、接触していない場合はステップＳ５７に進む。
（Ｓ５０）ステップ実行部１２５は、現ステップのｓｋｉｒｔ接触粒子リスト１５４の末尾に、相手粒子の粒子番号を追加する。また、ステップ実行部１２５は、ｓｋｉｒｔ粒子テーブル１３３から、相手粒子の粒子番号に対応するグローバルＩＤを読み出し、現ステップのグローバルＩＤリスト１５５の末尾に、相手粒子のグローバルＩＤを追加する。

（Ｓ５１）ステップ実行部１２５は、前ステップのグローバルＩＤリスト１４２を先頭から末尾に向かって走査し、相手粒子のグローバルＩＤを検索する。グローバルＩＤリスト１４２は、前ステップでｍａｉｎ領域に属する粒子を示している。

（Ｓ５２）ステップ実行部１２５は、ステップＳ５１において該当するグローバルＩＤが検出されたか判断する。該当するグローバルＩＤがある場合はステップＳ５６に進み、該当するグローバルＩＤが無い場合はステップＳ５３に進む。

（Ｓ５３）ステップ実行部１２５は、前ステップのグローバルＩＤリスト１４５を先頭から末尾に向かって走査し、相手粒子のグローバルＩＤを検索する。グローバルＩＤリスト１４５は、前ステップでｓｋｉｒｔ領域に属する粒子を示している。

（Ｓ５４）ステップ実行部１２５は、ステップＳ５３において該当するグローバルＩＤが検出されたか判断する。該当するグローバルＩＤがある場合はステップＳ５６に進み、該当するグローバルＩＤが無い場合はステップＳ５５に進む。

（Ｓ５５）ステップ実行部１２５は、現ステップの変位量リスト１５６の末尾に、変位量が０であることを示す数値を追加する。そして、ステップＳ５７に進む。
（Ｓ５６）ステップ実行部１２５は、前ステップの変位量リスト１４３から、ステップＳ５１で検出されたグローバルＩＤに対応する変位量を読み出す。変位量リスト１４３の所望の変位量には、グローバルＩＤリスト１４２のグローバルＩＤと同じインデックスを用いてアクセスすることが可能である。または、ステップ実行部１２５は、前ステップの変位量リスト１４６から、ステップＳ５３で検出されたグローバルＩＤに対応する変位量を読み出す。変位量リスト１４６の所望の変位量には、グローバルＩＤリスト１４５のグローバルＩＤと同じインデックスを用いてアクセスすることが可能である。ステップ実行部１２５は、読み出した変位量を、現ステップの変位量リスト１５６の末尾に追加する。

（Ｓ５７）ステップ実行部１２５は、ステップＳ４７においてｓｋｉｒｔ粒子テーブル１３３の全ての粒子番号を選択したか判断する。全ての粒子番号を選択した場合は接触粒子リスト生成が終了し、未選択の粒子番号がある場合はステップＳ４７に戻る。

第２の実施の形態の情報処理システムによれば、離散要素法の粒子シミュレーションにおいて、ｍａｉｎ領域に属する粒子それぞれに対して、接触粒子を列挙した接触粒子リストと接触粒子の接線方向の変位量を列挙した変位量リストが生成される。よって、ｍａｉｎ領域に属する粒子が接触粒子から受ける接触力を効率的に算出することができる。また、接触粒子リストが、ｍａｉｎ領域に属する接触粒子を示すｍａｉｎ接触粒子リストとｓｋｉｒｔ領域に属する接触粒子を示すｓｋｉｒｔ接触粒子リストに分割される。よって、接触粒子が存在する領域に応じて接触粒子の情報を効率的に管理することができる。

また、現ステップのｍａｉｎ接触粒子リストを生成する際には、前ステップのｍａｉｎ接触粒子リストに登録されている相手粒子について優先的に接触の有無が判定される。その後に、前ステップのｍａｉｎ接触粒子リストに登録されていない相手粒子について接触の有無が判定される。よって、ｍａｉｎ－ｍａｉｎ粒子について、前ステップの変位量リストから現ステップの変位量リストへの変位量の複写を効率化できる。

例えば、ｍａｉｎ粒子テーブルから現ステップのｍａｉｎ領域の接触粒子を検出し、その接触粒子を前ステップのｍａｉｎ接触粒子リストから探すようにすると、前ステップのｍａｉｎ接触粒子リストを繰り返し走査することになる。ｍａｉｎ接触粒子リストの平均長をＮとすると、その計算量はＯ（Ｎ^２）になる。これに対して、第２の実施の形態では、前ステップのｍａｉｎ接触粒子リストを１回走査することで、変位量を複写することになるｍａｉｎ－ｍａｉｎ粒子を特定でき、複写する変位量を前ステップの変位量リストから効率的に読み出すことができる。よって、前ステップのｍａｉｎ接触粒子リストの走査回数を減らすことができ、その計算量をＯ（Ｎ）に抑制することができる。

また、ｍａｉｎ領域に属する粒子それぞれに対して、接触粒子リストおよび変位量リストに加えて、接触粒子のグローバルＩＤを列挙したグローバルＩＤリストが生成される。よって、ｓｋｉｒｔ－ｍａｉｎ粒子、ｍａｉｎ－ｓｋｉｒｔ粒子およびｓｋｉｒｔ－ｓｋｉｒｔ粒子について、前ステップのｍａｉｎ接触粒子リストや前ステップのｓｋｉｒｔ接触粒子リストから所望の接触粒子を効率的に検索することができる。これにより、各粒子の粒子番号は前ステップと現ステップとで異なる可能性があるところ、ｍａｉｎ－ｍａｉｎ粒子以外の接触粒子についても変位量を効率的に複写することができる。

１０，１０－１ 情報処理装置
１１ 記憶部
１２ 処理部
１３，１３－１ 領域
１４，１６ 接触データ
１５，１５－１ 位置データ

Claims (9)

1. 第１の領域に属する第１の粒子に対して、第１の時刻において前記第１の粒子と接触している接触粒子と、前記第１の時刻までの接触中の間の接触粒子に対する前記第１の粒子の累積の変位量とを示す第１の接触データを記憶する記憶部と、
   前記第１の時刻より後の第２の時刻において前記第１の領域に属する粒子の位置を示す第１の位置データを算出し、前記第２の時刻において前記第１の領域と隣接する第２の領域に属する粒子の位置を示す第２の位置データを他の情報処理装置から取得し、
   前記第１の接触データと前記第１の位置データの算出結果とに基づいて、前記第１の時刻において前記第１の粒子と接触しており前記第１の領域に属し、前記第２の時刻において前記第１の領域に属する第２の粒子を検出し、
   前記第１の位置データに基づいて、前記第２の時刻において前記第１の粒子と前記第２の粒子とが接触しているか判定し、接触している場合、前記第１の接触データから前記第２の時刻に対応する第２の接触データに、前記第２の粒子の前記変位量を複写し、
   前記第１の位置データと前記第２の位置データとに基づいて、前記第２の粒子以外の前記第１の領域または前記第２の領域に属する粒子であって、前記第２の時刻において前記第１の粒子と接触している第３の粒子を検出し、
   前記第１の接触データから前記第３の粒子を検索し、前記第３の粒子が前記第１の接触データに登録されている場合、前記第１の接触データから前記第２の接触データに前記第３の粒子の前記変位量を複写する処理部と、
   を有する情報処理装置。
2. 前記第１の接触データでは、前記第１の時刻において各粒子に割り当てられた第１の識別情報を粒子の識別に使用し、前記第１の位置データでは、前記第２の時刻において各粒子に割り当てられた第２の識別情報を粒子の識別に使用し、
   前記第１の位置データの算出では、前記第１の時刻において前記第１の領域に属しかつ前記第２の時刻において前記第１の領域に属する粒子に対して、前記第１の識別情報と前記第２の識別情報との間の対応関係を示す対応データを生成し、
   前記第２の粒子の検出では、前記第１の接触データに含まれる前記第１の識別情報に対応する前記第２の識別情報を、前記対応データから検索する、
   請求項１記載の情報処理装置。
3. 前記第１の接触データでは、各粒子に割り当てられており前記第１の時刻と前記第２の時刻との間で変化しない第３の識別情報を粒子の識別に更に使用し、前記第１の位置データおよび前記第２の位置データは、前記第３の識別情報を含んでおり、
   前記第３の粒子の前記変位量の複写では、前記第３の識別情報を用いて前記第１の接触データから前記第３の粒子を検索する、
   請求項２記載の情報処理装置。
4. 前記第１の接触データは、前記第１の領域に属する接触粒子に対応する第１の対象領域データと、前記第２の領域に属する接触粒子に対応する第１の隣接領域データとを含み、
   前記第２の接触データは、前記第１の領域に属する接触粒子に対応する第２の対象領域データと、前記第２の領域に属する接触粒子に対応する第２の隣接領域データとを含み、
   前記第２の粒子の検出では、前記第１の対象領域データと前記第１の位置データの算出結果とに基づいて、前記第２の粒子を検出し、
   前記第２の粒子の前記変位量の複写では、前記第１の対象領域データから前記第２の対象領域データに前記第２の粒子の前記変位量を複写する、
   請求項１記載の情報処理装置。
5. 前記第１の接触データは、前記第１の領域に属する接触粒子に対応する第１の対象領域データと、前記第２の領域に属する接触粒子に対応する第１の隣接領域データとを含み、
   前記第２の接触データは、前記第１の領域に属する接触粒子に対応する第２の対象領域データと、前記第２の領域に属する接触粒子に対応する第２の隣接領域データとを含み、
   前記第２の時刻において前記第３の粒子が前記第１の領域に属する場合、前記第３の粒子の前記変位量の複写では、前記第１の隣接領域データから前記第３の粒子を検索し、前記第３の粒子が前記第１の隣接領域データに登録されている場合、前記第１の隣接領域データから前記第２の対象領域データに前記第３の粒子の前記変位量を複写する、
   請求項１記載の情報処理装置。
6. 前記第１の接触データは、前記第１の領域に属する接触粒子に対応する第１の対象領域データと、前記第２の領域に属する接触粒子に対応する第１の隣接領域データとを含み、
   前記第２の接触データは、前記第１の領域に属する接触粒子に対応する第２の対象領域データと、前記第２の領域に属する接触粒子に対応する第２の隣接領域データとを含み、
   前記第２の時刻において前記第３の粒子が前記第２の領域に属する場合、前記第３の粒子の前記変位量の複写では、前記第３の粒子が前記第１の対象領域データおよび前記第１の隣接領域データのうちの一方のデータに登録されている場合、前記一方のデータから前記第２の隣接領域データに前記第３の粒子の前記変位量を複写する、
   請求項１記載の情報処理装置。
7. 前記変位量は、前記第１の時刻までの接触中の間に前記第１の粒子が前記接触粒子に対して接線方向に移動した相対的な移動量の累積値である、
   請求項１記載の情報処理装置。
8. コンピュータが、
   第１の領域に属する第１の粒子に対して、第１の時刻において前記第１の粒子と接触している接触粒子と、前記第１の時刻までの接触中の間の接触粒子に対する前記第１の粒子の累積の変位量とを示す第１の接触データを取得し、
   前記第１の時刻より後の第２の時刻において前記第１の領域に属する粒子の位置を示す第１の位置データを算出し、前記第２の時刻において前記第１の領域と隣接する第２の領域に属する粒子の位置を示す第２の位置データを他のコンピュータから受信し、
   前記第１の接触データと前記第１の位置データの算出結果とに基づいて、前記第１の時刻において前記第１の粒子と接触しており前記第１の領域に属し、前記第２の時刻において前記第１の領域に属する第２の粒子を検出し、
   前記第１の位置データに基づいて、前記第２の時刻において前記第１の粒子と前記第２の粒子とが接触しているか判定し、接触している場合、前記第１の接触データから前記第２の時刻に対応する第２の接触データに、前記第２の粒子の前記変位量を複写し、
   前記第１の位置データと前記第２の位置データとに基づいて、前記第２の粒子以外の前記第１の領域または前記第２の領域に属する粒子であって、前記第２の時刻において前記第１の粒子と接触している第３の粒子を検出し、
   前記第１の接触データから前記第３の粒子を検索し、前記第３の粒子が前記第１の接触データに登録されている場合、前記第１の接触データから前記第２の接触データに前記第３の粒子の前記変位量を複写する、
   粒子シミュレーション方法。
9. 第１の領域に属する第１の粒子に対して、第１の時刻において前記第１の粒子と接触している接触粒子と、前記第１の時刻までの接触中の間の接触粒子に対する前記第１の粒子の累積の変位量とを示す第１の接触データを記憶する第１の情報処理装置と、
   前記第１の時刻より後の第２の時刻において前記第１の領域と隣接する第２の領域に属する粒子の位置を示す第２の位置データを、前記第１の情報処理装置に送信する第２の情報処理装置と、
   を有し、前記第１の情報処理装置は、
   前記第２の時刻において前記第１の領域に属する粒子の位置を示す第１の位置データを算出し、
   前記第１の接触データと前記第１の位置データの算出結果とに基づいて、前記第１の時刻において前記第１の粒子と接触しており前記第１の領域に属し、前記第２の時刻において前記第１の領域に属する第２の粒子を検出し、
   前記第１の位置データに基づいて、前記第２の時刻において前記第１の粒子と前記第２の粒子とが接触しているか判定し、接触している場合、前記第１の接触データから前記第２の時刻に対応する第２の接触データに、前記第２の粒子の前記変位量を複写し、
   前記第１の位置データと前記第２の位置データとに基づいて、前記第２の粒子以外の前記第１の領域または前記第２の領域に属する粒子であって、前記第２の時刻において前記第１の粒子と接触している第３の粒子を検出し、
   前記第１の接触データから前記第３の粒子を検索し、前記第３の粒子が前記第１の接触データに登録されている場合、前記第１の接触データから前記第２の接触データに前記第３の粒子の前記変位量を複写する、
   粒子シミュレーションシステム。

Images