(57)【特許請求の範囲】
【請求項１】
離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを求めるために用いる粒子データ演算用プログラムであって、
座標データによって特定される複数の粒子を含む領域を複数のセルに分割する領域分割ステップと、
前記領域分割ステップで得た前記複数のセルの中の１つのセルを基準にしたときに、他のセルが、基準にした前記１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように、前記他のセルのアドレスとして、基準にした前記１つのセルから離れたセル数分に対応する値を、基準にした前記１つのセルのアドレスに増減させたアドレスを付し、前記複数の粒子の座標データに基いて前記複数の粒子の中心が位置する複数の前記セルを特定し、前記複数の粒子の中心が位置する前記複数のセルを特定するデータをセルデータメモリに記憶させるマッピングステップと、
前記複数のセル内の各粒子に関して、前記複数のセルに付された前記アドレスに基いて、着目した１つの前記粒子が位置する前記セルと隣接する他の複数のセル内に他の粒子の中心がある場合には、前記着目した１つの粒子と前記他の粒子とを接触候補粒子ペアであると決定する接触候補粒子ペア決定ステップと、
前記接触候補粒子ペア決定ステップで決定された前記接触候補粒子ペアについてのデータを接触候補粒子ペアデータメモリに記憶させる接触候補粒子ペア記憶ステップと、
前記接触候補粒子ペアの粒子どうしが接触しているか否かを判定し、前記着目した１つの粒子の単位時間経過後の位置に基いて、前記着目した１つの粒子と接触している他の粒子及び前記着目した１つの粒子と接触する可能性の高い他の粒子との間の運動に関する物理量を演算する物理量演算ステップと、
前記物理量演算ステップの演算結果に基いて、前記複数の粒子の前記座標データを更新する座標データ更新ステップと、
前記物理量演算ステップで演算した結果から前記接触候補粒子ペアを構成する複数の前記粒子の中で、積算変位量が所定量以上になったものがあったか否かを判定する変位量判定ステップと、
前記変位量判定ステップで前記所定量以上の積算変位量を判定したときには、前記マッピングステップへ戻ってそれ以後の前記各ステップを実行する接触候補粒子ペアデータ更新ステップと、
前記変位量判定ステップで前記所定量以上の積算変位量を判定しないときには所定の時間間隔で前記物理量演算ステップ、前記座標データ更新ステップ及び変位量判定ステップを繰り返す繰り返しステップとを前記コンピュータに実行させるための粒子データ演算用プログラム。
【請求項２】
離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを求めるために用いる粒子データ演算用プログラムであって、
座標データによって特定される複数の粒子を含む領域を複数のセルに分割するセル分割ステップと、
前記領域分割ステップで得た前記複数のセルの中の１つのセルを基準にしたときに、他のセルが、基準にした前記１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように、前記他のセルのアドレスとして、基準にした前記１つのセルから離れたセル数分に対応する値を、基準にした前記１つのセルのアドレスに増減させたアドレスを付し、前記複数の粒子の座標データに基いて前記複数の粒子の中心が位置する複数の前記セルを特定し、前記複数の粒子の中心が位置する前記複数のセルを特定するデータをセルデータメモリに記憶させるマッピングステップと、
前記複数のセル内の各粒子に関して、前記複数のセルに付された前記アドレスに基いて、着目した１つの前記粒子が位置する前記セルと隣接する他の複数のセル内に他の粒子の中心がある場合には、前記着目した１つの粒子と前記他の粒子とを接触候補粒子ペアであると決定する接触候補粒子ペア決定ステップと、
前記接触候補粒子ペア決定ステップで決定された前記接触候補粒子ペアについてのデータを接触候補粒子ペアデータメモリに記憶させる接触候補粒子ペア記憶ステップと、
前記接触候補粒子ペアの粒子どうしが接触しているか否かを判定し、前記着目した１つの粒子の単位時間経過後の位置に基いて、前記着目した１つの粒子と接触している他の粒子及び前記着目した１つの粒子と接触する可能性の高い他の粒子との間の運動に関する物理量を演算する物理量演算ステップと、
前記物理量演算ステップの演算結果に基いて、前記複数の粒子の前記座標データを更新する座標データ更新ステップと、
前記演算ステップで演算した結果から前記接触候補粒子ペアを構成する複数の前記粒子の中で積算変位量が１／４セル以上になったものがあったか否かを判定する変位量判定ステップと、
前記変位量判定ステップで１／４セル以上の積算変位量を判定したときには、再度前記マッピングステップ及び前記接触候補決定ステップを実行して更新を行い、その後従前は接触候補粒子ペアであったが更新後に接触候補粒子ペアではなくなった接触候補粒子ペアのデータを前記接触候補粒子ペアデータメモリから削除し、従前接触候補粒子ペアではなかったが更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアのデータを前記接触候補粒子ペアデータメモリに追加し、残りの従前の接触候補粒子ペアのデータについてはそのまま保存を継続する接触候補粒子ペアデータ更新ステップと、
前記変位量判定ステップで１／４セル以上の積算変位量を判定しないときには所定の時間間隔で前記物理量演算ステップ、前記座標データ更新ステップ及び変位量判定ステップを繰り返す繰り返しステップとを前記コンピュータに実行させるための粒子データ演算用プログラム。
【請求項３】
前記マッピングステップでは、前記複数のセルが２次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０）とし、他のセルのアドレスを（±ｐ，±ｑ）で表現し、前記複数のセルが３次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０，０）とし、他のセルのアドレスを（±ｐ，±ｑ，±ｒ）で表現し、ただしｐ，ｑ，ｒは、基準とした前記１つのセルから各次元方向に数えてプラス方向またはマイナス方向に何番目のセルであるかを示す整数であり、
前記接触候補粒子ペア決定ステップでは、着目した１つの前記粒子が位置する前記セルのアドレスの数値と、他の粒子が位置する他のセルの前記アドレスの数値の差分の絶対値がすべて１以下になるセル内に位置する場合に、着目した前記１つの粒子と前記他の粒子とを前記接触候補粒子ペアと決定する請求項１または２に記載の粒子データ演算用プログラム。
【請求項４】
前記マッピングステップでは、前記複数のセルが２次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０）とし、他のセルのアドレスを（±ｐ，±ｑ）で表現し、前記複数のセルが３次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０，０）とし、他のセルのアドレスを（±ｐ，±ｑ，±ｒ）で表し、ただしｐ，ｑ，ｒは、基準とした前記１つのセルから各次元方向に数えてプラス方向またはマイナス方向に何番目のセルであるかを示す整数であり、
前記接触候補決定ステップでは、着目した１つの前記粒子が位置する前記セルのアドレスの数値と、他の粒子が位置する他のセルの前記アドレスの数値の差分の絶対値がすべて１以下になるセル内に位置する場合に、着目した前記１つの粒子と前記他の粒子とを前記接触候補粒子ペアと決定し、
前記接触候補粒子ペアデータ更新ステップでは、従前接触候補粒子ペアを構成していた２つの粒子が位置する２つの前記セルについてのそれぞれの前記アドレスの数値の差分の絶対値が１より大きいか否かを調べ、前記絶対値が１より大きいものを前記接触候補粒子ペアではなくなった接触候補粒子ペアと決定し、新たに前記接触候補決定ステップで決定した新たな接触候補粒子ペアについて、該新たな接触候補粒子ペアを構成する２つの前記粒子が位置していた従前の２つの前記セルについてのそれぞれの前記アドレスの数値の差分の絶対値が１より大きいか否かを調べ、前記絶対値が１より大きいものを前記更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアと決定することを特徴とする請求項１または２に記載の粒子データ演算用プログラム。
【請求項５】
前記接触候補粒子ペアデータメモリは、各接触候補粒子ペアについてのデータが次の接触候補粒子ペアについてのデータへの参照を持つリンクリスト構造によりデータを保存している請求項２または４に記載の粒子データ演算用プログラム。
【請求項６】
離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを求める粒子データ演算方法であって、
座標データによって特定される複数の粒子を含む領域を複数のセルに分割するセル分割ステップと、
前記領域分割ステップで得た前記複数のセルの中の１つのセルを基準にしたときに、他のセルが、基準にした前記１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように、前記他のセルのアドレスとして、基準にした前記１つのセルから離れたセル数分に対応する値を、基準にした前記１つのセルのアドレスに増減させたアドレスを付し、前記複数の粒子の座標データに基いて前記複数の粒子の中心が位置する複数の前記セルを特定し、前記複数の粒子の中心が位置する前記複数のセルを特定するデータをセルデータメモリに記憶させるマッピングステップと、
前記複数のセル内の各粒子に関して、前記複数のセルに付された前記アドレスに基いて、着目した１つの前記粒子が位置する前記セルと隣接する他の複数のセル内に他の粒子の中心がある場合には、前記着目した１つの粒子と前記他の粒子とを接触候補粒子ペアであると決定する接触候補粒子ペア決定ステップと、
前記接触候補粒子ペア決定ステップで決定された前記接触候補粒子ペアについてのデータを接触候補粒子ペアデータメモリに記憶させる接触候補粒子ペア記憶ステップと、
前記接触候補粒子ペアの粒子どうしが接触しているか否かを判定し、前記着目した１つの粒子の単位時間経過後の位置に基いて、前記着目した１つの粒子と接触している他の粒子及び前記着目した１つの粒子と接触する可能性の高い他の粒子との間の運動に関する物理量を演算する物理量演算ステップと、
前記物理量演算ステップの演算結果に基いて、前記複数の粒子の前記座標データを更新する座標データ更新ステップと、
前記演算ステップで演算した結果から前記接触候補粒子ペアを構成する複数の前記粒子の中で積算変位量が１／４セル以上になったものがあったか否かを判定する変位量判定ステップと、
前記変位量判定ステップで１／４セル以上の積算変位量を判定したときには、再度前記マッピングステップ及び前記接触候補決定ステップを実行して更新を行い、その後従前は接触候補粒子ペアであったが更新後に接触候補粒子ペアではなくなった接触候補粒子ペアのデータを前記接触候補粒子ペアデータメモリから削除し、従前接触候補粒子ペアではなかったが更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアのデータを前記接触候補粒子ペアデータメモリに追加し、残りの従前の接触候補粒子ペアのデータについてはそのまま保存を継続する接触候補粒子ペアデータ更新ステップと、
前記変位量判定ステップで１／４セル以上の積算変位量を判定しないときには所定の時間間隔で前記物理量演算ステップ、前記座標データ更新ステップ及び変位量判定ステップを繰り返す繰り返しステップとからなる粒子データ演算方法。
【請求項７】
前記マッピングステップでは、前記複数のセルが２次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０）とし、他のセルのアドレスを（±ｐ，±ｑ）で表現し、前記複数のセルが３次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０，０）とし、他のセルのアドレスを（±ｐ，±ｑ，±ｒ）で表現し、ただしｐ，ｑ，ｒは、基準とした前記１つのセルから各次元方向に数えてプラス方向またはマイナス方向に何番目のセルであるかを示す整数であり、
前記接触候補粒子ペア決定ステップでは、着目した１つの前記粒子が位置する前記セルのアドレスの数値と、他の粒子が位置する他のセルの前記アドレスの数値の差分の絶対値がすべて１以下になるセル内に位置する場合に、着目した前記１つの粒子と前記他の粒子とを前記接触候補粒子ペアと決定する請求項６に記載の粒子データ演算方法。
【請求項８】
前記マッピングステップでは、前記複数のセルが２次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０）とし、他のセルのアドレスを（±ｐ，±ｑ）で表現し、前記複数のセルが３次元で配置されている場合には、基準となる前記１つのセルのアドレスを（０，０，０）とし、他のセルのアドレスを（±ｐ，±ｑ，±ｒ）で表し、ただしｐ，ｑ，ｒは、基準とした前記１つのセルから各次元方向に数えてプラス方向またはマイナス方向に何番目のセルであるかを示す整数であり、
前記接触候補決定ステップでは、着目した１つの前記粒子が位置する前記セルのアドレスの数値と、他の粒子が位置する他のセルの前記アドレスの数値の差分の絶対値がすべて１以下になるセル内に位置する場合に、着目した前記１つの粒子と前記他の粒子とを前記接触候補粒子ペアと決定し、
前記接触候補粒子ペアデータ更新ステップでは、従前接触候補粒子ペアを構成していた２つの粒子が位置する２つの前記セルについてのそれぞれの前記アドレスの数値の差分の絶対値が１より大きいか否かを調べ、前記絶対値が１より大きいものを前記接触候補粒子ペアではなくなった接触候補粒子ペアと決定し、新たに前記接触候補決定ステップで決定した新たな接触候補粒子ペアについて、該新たな接触候補粒子ペアを構成する２つの前記粒子が従前位置した２つの前記セルについてのそれぞれの前記アドレスの数値の差分の絶対値が１より大きいか否かを調べ、前記絶対値が１より大きいものを前記更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアと決定することを特徴とする請求項６に記載の粒子データ演算方法。
【請求項９】
離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを演算により求める粒子データ演算装置であって、
座標データメモリに記憶された座標データによって特定される複数の粒子を含む領域を複数のセルに分割する領域分割手段と、
前記領域分割手段で得た前記複数のセルの中の１つのセルを基準にしたときに、他のセルが、基準にした前記１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように、前記他のセルのアドレスとして、基準にした前記１つのセルから離れたセル数分に対応する値を、基準にした前記１つのセルのアドレスに増減させたアドレスを付し、前記複数の粒子の座標データに基いて前記複数の粒子の中心が位置する複数の前記セルを特定し、前記複数の粒子の中心が位置する前記複数のセルを特定するデータをセルデータメモリに記憶させるマッピング手段と、
前記複数のセル内の各粒子に関して、前記複数のセルに付された前記アドレスに基いて、着目した１つの前記粒子が位置する前記セルと隣接する他の複数のセル内に他の粒子の中心がある場合には、前記着目した１つの粒子と前記他の粒子とを接触候補粒子ペアであると決定する接触候補粒子ペア決定手段と、
前記接触候補粒子ペア決定手段で決定した前記接触候補粒子ペアについてのデータを記憶する接触候補粒子ペアデータメモリと、
前記接触候補粒子ペアの粒子どうしが接触しているか否かを判定し、前記着目した１つの粒子の単位時間経過後の位置に基いて、前記着目した１つの粒子と接触している他の粒子及び前記着目した１つの粒子と接触する可能性の高い他の粒子との間の運動に関する物理量を演算する物理量演算手段と、
前記物理量演算手段の演算結果を記憶する物理量データメモリと、
前記物理量演算手段の演算結果に基いて、前記複数の粒子の前記座標データを更新する座標データ更新手段と、
前記演算量演算手段で演算した結果から前記接触候補粒子ペアを構成する複数の前記粒子の中で積算変位量が１／４セル以上になったものがあったか否かを判定する変位量判定手段と、
前記変位量判定手段が、１／４セル以上の積算変位量を判定したときには、再度前記マッピング手段及び前記接触候補粒子ペア決定手段を動作させて更新を行い、その後従前は接触候補粒子ペアであったが更新後に接触候補粒子ペアではなくなった接触候補粒子ペアのデータを前記接触候補粒子ペアデータメモリから削除し、従前接触候補粒子ペアではなかったが更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアのデータを前記接触候補粒子ペアデータメモリに追加し、残りの従前の接触候補粒子ペアのデータについてはそのまま保存を継続する接触候補粒子ペアデータ更新手段と、
前記変位量判定手段が１／４セル以上の積算変位量を判定しないときには、所定の時間間隔で前記物理量演算手段、前記座標データ更新手段及び変位量判定手段が繰り返し動作する粒子データ演算装置。
【請求項１０】
前記接触候補粒子ペアデータメモリは、各接触候補粒子ペアについてのデータが次の接触候補粒子ペアについてのデータへの参照を持つリンクリスト構造によりデータを保存している請求項９に記載の粒子データ演算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを演算により求めるために使用する粒子データ演算用プログラム、このプログラムを利用して動作する粒子データ演算装置及びこのプログラムを実行する際に使用される粒子データ演算方法に関するものである。
【背景技術】
【０００２】
「計算工学」Ｖｏｌ．８，Ｎｏ．３  ２００３の２８頁乃至３３頁に掲載された「ＤＥＭ最近の話題について」と題する解説（非特許文献１）に示されるように、崖崩れにおける地中の粒子（岩）の動きをシミュレーションする場合や、管を地中に埋設した場合における地中内における管の周囲の粒子（砂利石）の移動状況をシミュレーションする場合などに、離散要素法が用いられている。この離散要素法では、各粒子についての運動方程式を、短い時間ステップΔｔの間隔で数値積分しながら運動に関する物理量を更新する。すなわち離散要素法（ＤＥＭ）の処理は、大きく分けると、１）接触判定、２）接触力計算、３）数値積分による運動に関する物理量の更新、の３つから成り立っている。そして離散要素法では、１）の接触判定の結果、接触しているとみなされた粒子のペアについてのみ２）、３）の計算を進める。そのため、１）のプロセスを高速化すれば、全体の演算速度を速くすることができる。各粒子の接触・分離は、各ステップごとに逐次変化するため、各粒子の接触判定そのものは、各ステップ毎に必要であると考えられている。
【非特許文献１】「計算工学」Ｖｏｌ．８，Ｎｏ．３  ２００３の２８頁乃至３３頁
【発明の開示】
【発明が解決しようとする課題】
【０００３】
しかしながら各ステップ毎に、各粒子の接触判定を行うために、全ての粒子間の距離を演算することは、コンピュータにとって大きな負担となり、結果として演算速度を大幅に遅くすることになる。しかしながら、Δｔが小さな値であれば、１ステップ内における粒子の動きは小さい。そのため着目した粒子の近傍には存在しているが、着目した粒子とは接触していない近傍粒子と他の近傍粒子との組み合わせ、即ち接触候補粒子ペアの組み合わせ（Ｐａｉｒ  ｃａｎｄｉ）は、殆ど変わらないと考えられる。そこで１ステップ内における粒子の動きが小さいことを根拠として、ある数ステップの間、接触判定を行わずに、従前の接触候補粒子ペアをそのまま利用することが考えられている。実際に、一部の離散要素法を用いたシミュレーション技術では、高速化のために、このようなアルゴリズムが組み込まれている。しかしながらこのような接触判定の一部省略は、接触していない粒子間の力を計算したり、本当は接触している粒子間の接触力を計算しないことになるため、シミュレーション結果の信頼性を大幅に低下させる問題がある。
【０００４】
本発明の目的は、接触判定の演算速度を速めることができて、しかもシミュレーションの精度の低下を阻止することができる粒子データ演算用プログラム並びに粒子データ演算装置及び方法を提供することにある。
【０００５】
本発明の他の目的は、接触候補粒子ペアの更新のタイミングを合理的に決定して、上記目的を達成する粒子データ演算用プログラム並びに粒子データ演算装置及び方法を提供することにある。
【０００６】
本発明の更に他の目的は、メモリの書き換え量を減らして更に演算の高速化を図ることができる粒子データ演算用プログラム並びに粒子データ演算装置及び方法を提供することにある。
【課題を解決するための手段】
【０００７】
本発明は、離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを求めるために用いる粒子データ演算用プログラムを対象とする。本発明のプログラムでは、以下の領域分割ステップと、マッピングステップと、接触候補粒子ペア決定ステップと、接触候補粒子ペア記憶ステップと、接触候補粒子ペア決定ステップと、接触候補粒子ペア記憶ステップと、物理量演算ステップと、座標データ更新ステップと、変位量判定ステップと、接触候補粒子ペアデータ更新ステップと、繰り返しステップとをコンピュータに実行させる。
【０００８】
領域分割ステップでは、座標データによって特定される複数の粒子を含む、領域を複数のセルに分割する。ここで座標データで粒子を特定することとは、領域に対して設定した複数の座標軸の数値で粒子の位置を特定することである。したがって粒子の座標データは、粒子の中心の正確な位置データを提供するものである。領域を分割して得た複数のセルは、領域を細分化した領域単位であって、２次元で配置されていても、３次元で配置されていてもよい。扱い易いセルの形状は、２次元であれば正方形であり、３次元であれば正立方体または正六面体であるが、他の形状であってもよいのは勿論である。そして粒子とは、シミュレーションによって運動をシミュレートする対象物である。なお設定するセルの大きさは、対象とする粒子が内部に収納可能な大きさとするのが好ましい。セルが正方形で粒子が円であるとすると、セルの一辺の長さは、粒子の直径よりも大きくするのが好ましい。
【０００９】
マッピングステップでは、複数のセルについて、複数のセルの中のある１つのセルを基準にしたときに他のセルが前述の１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように、他のセルのアドレスとして、基準にした１つのセルから離れたセル数分に対応する値を、基準にした１つのセルのアドレスに増減させたアドレスを付して、セルデータメモリを作る。次に、前述の粒子の座標データに基いて複数の粒子の中心が位置する複数のセルを特定し、粒子の中心が位置するセルを特定するデータをセルデータメモリに記憶させる。
【００１０】
次に、接触候補粒子ペア決定ステップでは、複数のセル内の各粒子に関して、複数のセルに付されたアドレスに基いて、着目した１つの粒子の中心が位置するセルに対して、隣接する他の複数のセル内に他の粒子の中心がある場合には、着目した１つの粒子と前記他の粒子を接触候補粒子ペアと決定する。ここで「隣接する」とは、着目した１つのセルと他のセルとの間に、更に別のセルが存在しない状態を言う。したがってセルの形状が２次元の正方形の場合には、着目した１つのセルの周囲に位置する８個のセルが「隣接する他の複数のセル」となる。この接触候補粒子ペア決定ステップは、複数のセル内の全ての粒子に関して実施される。したがって粒子の数が多くなれば、かなりの数の粒子候補ペアが決定されることになる。そして接触候補粒子ペア記憶ステップでは、接触候補粒子ペア決定ステップで決定された接触候補粒子ペアについてのデータを接触候補粒子ペアデータメモリに記憶させる。
【００１１】
なお前述のマッピングステップでは、例えば、複数のセルが２次元で配置されている場合には、基準となる１つのセルのアドレスを（０，０）とし、他のセルのアドレスを（±ｐ，±ｑ）で表現することができる。また複数のセルが３次元で配置されている場合には、基準となる１つのセルのアドレスを（０，０，０）とし、他のセルのアドレスを（±ｐ，±ｑ，±ｒ）で表現することができる。ただしｐ，ｑ，ｒは、基準とした１つのセルから各次元方向（一次元方向、二次元方向、三次元方向）に数えてプラス方向またはマイナス方向に何番目のセルであるかを示す整数である。このようにアドレスを決めた場合、接触候補粒子ペア決定ステップでは、着目した１つの粒子が位置するセルのアドレスの数値と、他の粒子が位置する他のセルのアドレスの数値の差分の絶対値が、すべて１以下になるセル内に位置する場合に、着目した１つの粒子と他の粒子とを接触候補粒子ペアと決定することができる。このようにすると、アドレスの数値を用いて演算により簡単に接触候補粒子ペアを決定することができる。
【００１２】
物理量演算ステップでは、接触候補粒子ペアの粒子どうしが接触しているか否かを判定する。この判定は、粒子の座標データに基いて演算により行う。予め粒子の大きさが判っているので、粒子の座標データから、接触候補粒子ペアの２つの粒子が接触しているか否かは、演算で判定できる。そして１つの粒子の単位時間経過後の位置に基いて、この１つの粒子と接触している他の粒子及びこの一つの粒子と接触する可能性の高い他の粒子との間の接触力計算を行い、この接触力の総和の計算値（数値積分）に基いて、粒子の運動に関する物理量（位置や速度）を演算して物理量メモリに記憶する。この演算は、離散要素法（ＤＥＭ）の処理における、１）接触判定、２）接触力計算、３）数値積分による運動に関する物理量の更新に相当するものである。本発明では、接触候補粒子ペアについてのみ、１）接触判定、２）接触力計算、３）数値積分による運動に関する物理量の更新が行われる。そのためすべての粒子間でこれらの処理を行う場合と比べて、演算速度は速くなる。
【００１３】
座標データ更新ステップでは、物理量演算ステップの演算結果（特に粒子の位置）に基いて、各粒子の座標データを更新する。更新された粒子の座標データは、座標データメモリに記憶される。なお接触候補粒子ペア以外の粒子については、物理量の更新がないため、従前の座標データがそのまま使用される。新たに接触候補粒子ペアを構成することになった粒子についても、この更新が行われる。
【００１４】
次に、変位量判定ステップでは、物理量演算ステップで演算した結果から、接触候補粒子ペアを構成する複数の粒子の中で、積算変位量が所定量以上（具体的には、１／４セル以上）になったものがあったか否かを判定する。ここで積算変位量とは、接触候補粒子ペアデータ更新ステップで更新動作をするまでの間、物理量演算ステップで演算した変位量を積算した値を意味する。そして、変位量判定ステップで所定量以上（具体的には１／４セル以上）の積算変位量を判定したときには、接触候補粒子ペアデータ更新ステップでは、マッピングステップへ戻ってそれ以後の各ステップを実行する。本発明は、接触候補粒子ペアを構成する粒子のいずれかにおいて、所定量以上（具体的には１／４セル以上）の積算変位量が発生すると、接触候補粒子ペアの中で前のステップでは接触していなかった接触候補粒子ペアにおいて接触状態が発生する可能性が高くなることを基礎とする。本発明では、このような状況になったときに、マッピングステップに戻って、接触候補粒子ペアのペアデータの更新を行う。本発明のように、接触候補粒子ペアの更新タイミングを定めると、接触候補粒子ペアについての接触状態の発生可能性が高まった状況において、接触候補粒子ペアのペアデータの更新をすることになるため、接触判定の一部省略を行って演算速度を速めても、シミュレーションの精度が低下することはない。なお変位量判定ステップで所定量以上（具体的には１／４セル以上）の積算変位量を判定しないときには、所定の時間間隔で物理量演算ステップ、座標データ更新ステップ及び変位量判定ステップを繰り返す、繰り返しステップをコンピュータに実行させる。すなわち本発明では、新たな接触の発生の可能性がない状況下では、従前の接触候補粒子ペアについて継続して物理量演算を行う。したがって接触候補粒子ペアの更新回数を減らしても、演算精度が低下することはない。
【００１５】
接触候補粒子ペアのペアデータの更新に際しては、具体的に、以下のようにすると、メモリの書き換え速度を速めることができて、結果として演算速度を更に速くすることができる。すなわち変位量判定ステップで、所定量以上の積算変位量を判定したときには、再度マッピングステップ及び接触候補決定ステップを実行して、次の接触候補粒子ペアデータ更新ステップを実行する。接触候補粒子ペアデータ更新ステップでは、従前は接触候補粒子ペアであったが更新後に接触候補粒子ペアではなくなった接触候補粒子ペアのデータを、接触候補粒子ペアデータメモリから削除する。そして従前は接触候補粒子ペアではなかったが更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアのデータを接触候補粒子ペアデータメモリに追加する。残りの従前の接触候補粒子ペアのデータについてはそのまま保存を継続する。このようにすると、データの削除と追加という操作だけで、ペアデータの更新を行うことができるので、データの書き換え作業が容易になり、しかも演算速度が速くなる。このような接触候補粒子ペアデータメモリは、各接触候補粒子ペアについてのデータが次の接触候補粒子ペアについてのデータへの参照を持つリンクリスト構造によりデータを保存しているのが好ましい。
【００１６】
なお接触候補粒子ペアデータ更新ステップでは、例えば具体的に以下のようにして削除する接触候補粒子ペアと追加する接触候補粒子ペアとを決定することができる。まず従前接触候補粒子ペアを構成していた２つの粒子が位置する２つのセルについてのそれぞれのアドレスの数値の差分の絶対値が１より大きいか否かを調べる。そして絶対値が１より大きいものを接触候補粒子ペアではなくなった接触候補粒子ペアと決定する。また新たに接触候補決定ステップで決定した新たな接触候補粒子ペアについて、この新たな接触候補粒子ペアを構成する２つの粒子が従前位置していた２つのセルについてのそれぞれのアドレスの数値の差分の絶対値が１より大きいか否かを調べる。そして、絶対値が１より大きいものを更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアと決定する。
【００１７】
本発明の方法を実施して、離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを演算により求める粒子データ演算装置は、セル分割手段と、マッピング手段と、接触候補粒子ペア決定手段と、接触候補粒子ペア記憶手段と、物理量演算手段と、座標データ更新手段と、変位量判定手段と、接触候補粒子ペアデータ更新手段とを有する。セル分割手段は、座標データメモリに記憶された座標データによって特定される複数の粒子を含む領域を複数のセルに分割する。マッピング手段は、領域分割手段で得た複数のセルの中の１つのセルを基準にしたときに他のセルが基準にした１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように、他のセルのアドレスとして、基準にした１つのセルから離れたセル数分に対応する値を、基準にした１つのセルのアドレスに増減させたアドレスを付し、複数の粒子の座標データに基いて複数の粒子の中心が位置する複数のセルを特定し、複数の粒子の中心が位置する複数のセルを特定するデータをセルデータメモリに記憶させる。接触候補粒子ペア決定手段は、複数のセル内の各粒子に関して、複数のセルに付されたアドレスに基いて、着目した１つの粒子が位置するセルと隣接する他の複数のセル内に他の粒子の中心がある場合には、着目した１つの粒子と前記他の粒子とを接触候補粒子ペアであると決定する。接触候補粒子ペアデータメモリは、接触候補粒子ペア決定ステップで決定された接触候補粒子ペアについてのデータを記憶する。物理量演算手段は、接触候補粒子ペアの粒子どうしが接触しているか否かを判定し、着目した１つの粒子の単位時間経過後の位置に基いて、着目した１つの粒子と接触している他の粒子及び着目した１つの粒子と接触する可能性の高い他の粒子との間の運動に関する物理量を演算する。物理量演算手段の演算結果は、物理量メモリに記憶される。座標データ更新手段は、物理量演算手段の演算結果に基いて、複数の粒子の座標データを更新する。変位量判定手段は、演算量演算手段で演算した結果から接触候補粒子ペアを構成する複数の粒子の中で積算変位量が１／４セル以上になったものがあったか否かを判定する。そして接触候補粒子ペアデータ更新手段は、変位量判定手段が、１／４セル以上の積算変位量を判定したときには、再度マッピング手段及び接触候補決定手段を動作させて更新を行い、その後従前は接触候補粒子ペアであったが更新後に接触候補粒子ペアではなくなった接触候補粒子ペアのデータを接触候補粒子ペアデータメモリから削除し、従前接触候補粒子ペアではなかったが更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアのデータを接触候補粒子ペアデータメモリに追加し、残りの従前の接触候補粒子ペアのデータについてはそのまま保存を継続する。そして変位量判定手段が１／４セル以上の積算変位量を判定しないときには、所定の時間間隔で物理量演算手段、座標データ更新手段及び変位量判定手段が繰り返し動作する。
【発明の効果】
【００１８】
本発明によれば、接触判定そのものを間引くのではなく、接触候補粒子ペアの更新を合理的に間引くことにより、結果として接触判定の一部省略を行って演算速度を速めることができて、しかもシミュレーションの精度の低下を阻止することができる。また本発明によれば、接触候補粒子ペアの更新のタイミングを合理的に決定することができるので、メモリの書き換え量を減らして更に演算の高速化を図ることができる。
【発明を実施するための最良の形態】
【００１９】
以下図面を参照して、本発明の粒子データ演算方法を実施するために用いる本発明の粒子データ演算用プログラムの実施の形態の一例を詳細に説明する。
【００２０】
離散要素法（ＤＥＭ）の処理は、大きく分けると、１）粒子の接触判定、２）粒子間の接触力計算、３）数値積分による粒子の運動に関する物理量の更新、の３つから成り立っている。そして離散要素法（ＤＥＭ）では、１）の接触判定の結果、接触しているとみなされた粒子のペアについてのみ２）、３）の計算に進めるため、１）の粒子の接触判定のプロセスが律速段階となる。すなわち、上記１）乃至３）のステップのうちで、演算速度を最も遅くするステップが１）のステップであり、そのために１）のステップの速度で全体の反応速度が支配されることになる。
【００２１】
本発明は、この点に鑑みてなされたものである。粒子間の接触判定は、基本的に２つの粒子間の距離計算に必要な各座標成分差の２乗演算とその平方根演算を伴うため、比較的計算機には負荷の重い計算である。したがって、できるだけ、接触判定の回数は少ないにこしたことはない。つまり、「接触していない」と判定されるケースが少ないアルゴリズムが「効率の良い接触判定ルーチン」となる。そこで、接触判定の前に、接触候補粒子ペアの絞込みが必要となる。
【００２２】
図１は、コンピュータにおいて実行される本発明の粒子データ演算用プログラムの一実施の形態のアルゴリズムをメモリと一緒に示したフローチャートであり、図２は図１のフローチャートの一部のステップのアルゴリズムをより詳細に示したフローチャートである。この実施の形態のプログラムでは、領域分割ステップと、マッピングステップと、接触候補粒子ペア決定ステップと、接触候補粒子ペア記憶ステップと、物理量演算ステップと、座標データ更新ステップと、変位量判定ステップと、接触候補粒子ペアデータ更新ステップと、繰り返しステップとをコンピュータに実行させる。
【００２３】
本実施の形態のプログラムでは、接触候補粒子ペア（Ｐａｉｒ  ｃａｎｄｉ）の絞込みを、周知の近傍セル法を用いて行う。近傍セル法を用いるためには、粒子が存在する領域を複数のセルに分割する必要がある。粒子が存在する領域とは、例えば、崖崩れをシミュレーションする場合に離散要素法を適用するとすれば、崖の所定領域がこの粒子が存在する領域となる。そして地中に埋まっている石または岩が粒子となる。そこで最初に領域分割ステップＳＴ１では、シミュレーションを行うために予め定めた座標データによって特定される複数の粒子を含む領域を複数のセルに分割する。座標データは、座標データメモリＭ１に記憶されている。ここで座標データは、対象とする領域に対して設定した二次元または三次元の座標に基いて、シミュレーションで用いる複数の粒子の位置を特定するデータである。
【００２４】
シミュレーションを行う場合、領域を分割する複数のセルは、２次元で配置されていても、３次元で配置されていてもよい。扱い易いセルの形状は、２次元であれば正方形であり、３次元であれば正立方体である。以下の説明では、説明を簡単にして理解を容易にするために、複数のセルが２次元状態で配置されている場合を例にして説明をする。図３には、２次元に配置した複数のセルの一例を示してある。図３において円の図形が粒子であり、円の図形中のｉの文字は、粒子に付された名前である。なお設定するセルの大きさは、対象とする粒子が内部に収納可能な大きさとするのが好ましい。セルが正方形で粒子が円であるとすると、セルの一辺の長さは、粒子の直径よりも大きくするのが好ましい。
【００２５】
次にマッピングステップＳＴ２では、複数のセルについてある１つのセルを基準にしたときに他のセルが該１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように、他のセルのアドレスとして、基準にした１つのセルから離れたセル数分に対応する値を、基準にした１つのセルのアドレスに増減させたアドレスを付す。そして、前述の粒子の座標データに基いて、粒子の中心が位置するセルを特定し、粒子の中心が位置するセルのアドレスをセルデータメモリＭ２に記憶させる。実際的には、アドレスの他にそのセル内に中心が存在する粒子を特定することができる名称を付して、セルデータメモリＭ２に記憶させている。
【００２６】
図３に示す例では、粒子ｉが存在するセルにＣikの名称を付している。以後、着目する１つの粒子を着目粒子と言う。着目粒子の番号をｉとし、粒子ｉが所属するセル（ｃｅｌｌ）のｋ方向（２次元の場合ｋ＝±ｐ，±ｑ、３次元の場合ｋ＝±ｐ，±ｑ，±ｒ）の番号をＣikとする。ここでは基準となる１つのセルのアドレスを（０，０）として、他のセルのアドレスを（±ｐ，±ｑ）で表現する。また複数のセルが３次元で配置されている場合には、基準となる１つのセルのアドレスを（０，０，０）とし、他のセルのアドレスを（±ｐ，±ｑ，±ｒ）で表現することができる。ただしｐ，ｑ，ｒは、基準とした１つのセルから各次元方向（一次元方向、二次元方向、三次元方向）に数えてプラス方向またはマイナス方向に何番目のセルであるかを示す整数である。
【００２７】
次に、接触候補粒子ペア決定ステップＳＴ３では、複数のセル内の各粒子に関して、セルデータメモリＭ２に記憶されているアドレスに基いて、着目した１つの着目粒子が位置するセルと隣接する他の複数のセル内に他の粒子の中心がある場合には、この１つの着目粒子とこの他の粒子を接触候補粒子ペアとして決定する。ここで「隣接する」とは、１つのセルと他のセルとの間に、更に別のセルが存在しない状態を言う。したがってセルの形状が２次元の正方形の場合には、１つのセルの周囲に位置する８個のセルが「隣接する他の複数のセル」となる。この接触候補粒子ペア決定ステップＳＴ３は、複数のセル内の全ての粒子に関して実施される。したがって粒子の数が多くなれば、多数の数の粒子候補ペアが決定されることになる。そして接触候補粒子ペア記憶ステップでは、接触候補粒子ペア決定ステップで決定された接触候補粒子ペアについてのデータを接触候補粒子ペアデータメモリＭ３に記憶させる。
【００２８】
接触候補粒子ペア決定ステップＳＴ３では、１つの着目粒子が位置するセルのアドレスの数値と、他のセルのアドレスの数値の差分の絶対値が、すべて１以下になるセル内に他の粒子の中心が位置する場合には、その着目粒子と他の粒子とを接触候補粒子ペアとして決定する。このようにすると、アドレスの数値を用いて演算により簡単に接触候補粒子ペアを決定することができる。なお接触候補粒子ペアの決定方法及び更新のタイミングについては、後に詳しく説明する。
【００２９】
物理量演算ステップ（ＳＴ４乃至ＳＴ８）では、まず接触候補粒子ペアの粒子どうしが接触しているか否かを判定する（ＳＴ４及びＳＴ５）。この判定ステップＳＴ５では、粒子の座標データに基いて演算により接触の判定の有無を行う。そこでステップＳＴ４では、接触候補粒子ペアデータメモリＭ４から接触候補粒子ペアのデータを入手し、また座標データメモリＭ１から接触候補粒子ペアを構成する粒子の座標データを入手する。そして判定ステップＳＴ５では、予め設定した粒子の大きさが判っているので、各粒子の座標データに基いて、接触候補粒子ペアに関して、２つの粒子の中心間の距離を演算する。演算した接触候補粒子ペアの２つの粒子の中心間の距離から、接触候補粒子ペアの２つの粒子が接触しているか否かを判定する。
【００３０】
次に、接触候補粒子ペアを構成する１つの粒子（着目粒子）の単位時間経過後の位置に基いて、この着目粒子と接触している他の粒子及びこの一つの粒子と接触する可能性の高い他の粒子（以下近傍粒子とも言う）との間の接触力計算を行う（ステップＳＴ６）。着目粒子の単位時間経過後の位置に関しては、例えば、最初は粒子に作用する重力を粒子に加わる外力として、粒子が位置を変える（変位）ものとする。各粒子に重力が作用しているものとして、着目粒子に作用する力（これを接触力という）を計算する。すなわち着目粒子に加わる重力ベクトル、１以上の近傍粒子から加わる力のベクトルが、その着目力粒子に加わる接触力となる。接触力の計算に当たっては、物理量データメモリＭ４に記憶されている各粒子についての従前の接触力データ（接触力総和）も当然にして用いる。そしてこれらの接触力の合成ベクトルを、接触力総和として計算する（ステップＳＴ７）。接触力総和の計算結果は、物理量データメモリＭ４に記憶される。そして接触力総和計算の結果は、時間積分される（ステップＳＴ８）。時間積分により得た数値積分により、着目粒子と近傍粒子との運動に関する物理量（具体的には位置と速度）を演算して（ステップＳＴ９）、その結果を物理量データメモリＭ４に記憶する。
【００３１】
物理量データメモリＭ４に記憶している従前の物理量と今回の演算結果の物理量が相違する場合には、物理量の更新を行う（ステップＳＴ９）。そして物理量の更新に基いて座標データメモリＭ１に記憶されている各粒子の座標データを更新する（座標データ更新ステップ）。すなわち演算により求めた物理量によって、座標データメモリＭ１に記憶されている粒子の座標データの変更を行う。なお接触候補粒子ペア以外の粒子については、物理量の更新がないため、従前の座標データがそのまま保持されて使用される。
【００３２】
ここまでの演算は、離散要素法（ＤＥＭ）の処理における、１）接触判定、２）接触力計算、３）数値積分による運動に関する物理量の更新に相当するものである。本例では、接触候補粒子ペアを構成する粒子についてのみ、これらの演算を行って、１）接触判定、２）接触力計算、３）数値積分による運動に関する物理量の更新を行う。このように接触候補粒子ペアについてのみ演算を行えば、領域内のすべての粒子間でこれらの演算を行う場合と比べて、演算速度は大幅に速くなる。演算を行う単位時間は、シミュレーションの対象によって任意に定めることになる。
【００３３】
すべての接触候補粒子ペアを構成する粒子についての演算が終了したことを判定すると（ステップＳＴ１０）、変位量判定ステップＳＴ１１へと進む。この変位量判定ステップＳＴ１１では、前述の物理量演算ステップ（ＳＴ４〜ＳＴ９）で演算した結果から、接触候補粒子ペアを構成する複数の粒子の中で積算変位量が所定量以上の変位、具体的には１／４セル以上の変位になったものがあったか否かを判定する。ここで「積算変位量が所定量以上の変位」は、新たに接触する粒子が発生する可能性が実質的に無い変位量を基準にして定める。すなわち粒子の直径と一つのセルの輪郭を構成する一辺の寸法とを考慮し、どの程度粒子が変位すると、新たに接触する粒子ペアが発生するかを検討して定めることになる。一つのセルには最も大きな粒子を一つ含む必要性から、セルの一辺は、最大粒子の直径より僅かに大きければよく、例えば、本実施の形態では、粒子の直径ｒと一つのセルの輪郭を構成する一辺の寸法Ｌとの比（ｒ：Ｌ）が、１：１．０１の関係にある場合を考える。なおこの比は、本実施の形態におけるものであって、この比の関係にのみ本発明が限定されるものではない。
【００３４】
そしてセルが２次元の正方形であれば、この一つの辺は、正方形の一辺とし、セルが正六面体であれば、この一つの辺は正六面体の一面の一辺とする。このように仮定した場合、「所定量以上」とは「１／４セル以上」とするのが好ましい。
【００３５】
本実施の形態では、変位量判定ステップ（ＳＴ１１）で１／４セル以上の変位を判定したときには、接触候補粒子ペアデータ更新ステップＳＴ１２へと進む。また変位量判定ステップで１／４セル以上の変位を判定しないときには、ステップＳＴ１３へと進み、所定の時間間隔で物理量演算ステップ、座標データ更新ステップ（ＳＴ９）及び変位量判定ステップ（ＳＴ１０，ＳＴ１１）を繰り返す、繰り返しステップ（ＳＴ４〜ＳＴ１１，ＳＴ１３）をコンピュータに実行させる。所定の時間間隔は、ステップＳＴ１３で更新される。後に詳しく説明するように、接触候補粒子ペアの粒子の変位量が、１／４セルより少ない場合には、新たに接触する粒子が発生する可能性が実質的に無いということを前提にしている。したがって新たな接触の発生の可能性がない状況下では、従前の接触候補粒子ペアについて継続して物理量演算を行う。このように接触候補粒子ペアの更新回数を減らしても、演算精度が低下することはない。
【００３６】
積算変位量が、１／４セル以上になったことをステップＳＴ１１で判定すると、接触候補粒子ペアの中で前のステップでは接触していなかった接触候補粒子ペアにおいて接触状態が発生する可能性が高くなる。そこでこのような状況になったときに、ステップＳＴ１２でセルデータ及び接触候補粒子ペアのペアデータの更新を行う。具体的には、マッピングステップＳＴ２に戻って、接触候補粒子ペアのペアデータの更新を行う。このように、接触候補粒子ペアの更新タイミングを定めると、接触候補粒子ペアについての接触状態の発生可能性が高まった状況において、接触候補粒子ペアのペアデータの更新をすることになるため、接触判定の一部省略を行って演算速度を速めても、シミュレーションの精度が低下することはない。
【００３７】
図２は、更新の際におけるマッピングステップ及び接触候補決定ステップを実行するプログラムのルーチンのアルゴリズムの詳細を示すフローチャートである。このアルゴリズムでは、接触候補粒子ペアのペアデータの更新に際して、以下のようにすることにより、メモリの書き換え速度を速める。すなわち前述の変位量判定ステップで１／４セル以上の変位を判定したときには、再度マッピングステップ及び接触候補決定ステップを実行して、次の接触候補粒子ペアデータ更新ステップを実行する。接触候補粒子ペアデータ更新ステップでは、従前は接触候補粒子ペアであったが更新後に接触候補粒子ペアではなくなった接触候補粒子ペアのデータを、接触候補粒子ペアデータメモリから削除する。そして従前は接触候補粒子ペアではなかったが更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアのデータを接触候補粒子ペアデータメモリに追加する。残りの従前の接触候補粒子ペアのデータについてはそのまま保存を継続する。このようにすると、データの削除と追加という操作だけで、ペアデータの更新を行うことができるので、データの書き換え作業が容易になり、しかも演算速度が速くなる。本実施の形態では、接触候補粒子ペアデータメモリＭ３は、各接触候補粒子ペアについてのデータが次の接触候補粒子ペアについてのデータへの参照を持つリンクリスト構造によりデータを保存している。
【００３８】
なお接触候補粒子ペアデータ更新ステップでは、例えば具体的に以下のようにして削除する接触候補粒子ペアと追加する接触候補粒子ペアとを決定する。ステップＳＴ２１では、更新された座標データに基いて粒子のセルマッピングを実行する。このステップＳＴ２１は、図１のステップ２と同じである。ステップＳＴ２１のセルマッピングでは、粒子が存在するセルに対して名称をつける。本実施の形態では、便宜的に各セルに対して「Ｃik」の名前または番号をつけている。すなわちｉは着目粒子の名前であり、ｋは着目粒子ｉが存在するセルのアドレスである。ここで粒子ｉが所属するセル（ｃｅｌｌ）のｋ方向のアドレスは、２次元の場合ｋ＝±ｐ，±ｑ、３次元の場合ｋ＝±ｐ，±ｑ，±ｒとなる。ここでは基準となる１つのセルのアドレスを（０，０）として、他のセルのアドレスを（±ｐ，±ｑ）で表現する。また複数のセルが３次元で配置されている場合には、基準となる１つのセルのアドレスを（０，０，０）とし、他のセルのアドレスを（±ｐ，±ｑ，±ｒ）で表現する。ただしｐ，ｑ，ｒは、基準とした１つのセルから各次元方向（一次元方向、二次元方向、三次元方向）に数えてプラス方向またはマイナス方向に何番目のセルであるかを示す整数である。図４に示すように、着目粒子ｉが含まれるセルを中心として±ｐ，±ｑ（以下の説明では、便宜的に２次元でアドレスを表現する）の方向に±１の領域（図中のハッチングを施した部分）を、着目粒子ｉの近傍セルと定義する。このとき、近傍セル内に存在する粒子を着目粒子ｉの近傍粒子ｊとし、粒子ｊが所属するセルの±ｐ，±ｑの方向のセルの名称をＣjkとする。
【００３９】
本実施の形態では、以下の考え方で、前述の接触候補粒子ペアの更新タイミングを決定している。図１におけるステップＳＴ４からステップＳＴ１３までの複数の計算ステップがｎ回繰り返された後の着目粒子ｉの周辺には、以下の３つのような状況の発生が想定される。
【００４０】
１．着目粒子ｉ、と近傍粒子ｊは継続して近傍関係にある（継続）
２．着目粒子ｉ、と近傍粒子ｊは互いに離れてしまい、もはや近傍関係にない（解約）。
【００４１】
３．ｎ計算ステップ前（ｎは正の整数）には着目粒子ｉとは近傍関係になかった粒子ｊ´が着目粒子ｉに接近し、着目粒子ｉと粒子ｊ´が新たに近傍関係となり新しい接触候補粒子ペアができる（新規）。
【００４２】
状況１は、ｎ計算ステップ前の接触候補粒子ペアデータを継続して利用することで、その中に見出すことができる。
【００４３】
また状況２も、ｎ計算ステップ前の接触候補粒子ペアデータの中に見出すことができる。着目粒子ｉと近傍粒子ｊがもはや近傍関係にないことが確認された時点で、その組のデータだけをメモリーから消去すればよい。
【００４４】
さらに状況３の関係は、ｎ計算ステップ前の接触候補粒子ペアデータの中に見出すことはできない。従ってこの状況が発生するタイミングで接触候補粒子ペアの更新が必要である。逆に言うと、すべての着目粒子に対してこの状況３が発生しない限り、接触候補粒子ペアは更新する必要がない。したがって、接触候補粒子ペアの更新タイミングは、この状況３が発生するタイミングと設定すればよい。しかしながら、状況３は、更新前の接触候補粒子ペアには見出せない関係であるため、更新前の近傍関係以外の基準で決めなければならない。そこで、図２に示したプログラムのアルゴリズムでは、「状況３が発生するタイミング」を「状況３が発生するであろうタイミング」に置き換えて、接触候補粒子ペアの更新基準を決めている。
【００４５】
ある時刻ｔにセルＣikに存在した着目粒子ｉが、Δｔ秒間に図５のように移動したとき、移動後に着目粒子が存在するセルをＣik´とする。そして着目粒子ｉの移動に伴うセルの移動量ベクトルをｄikとすると、移動量ベクトルｄikは次のように表される。
【００４６】
ｄik ＝Ｃik´−Ｃik
同様に着目粒子ｉの近傍粒子ｊのΔｔ秒間の移動についても考え、図６のようにｊの移動後に存在するセルをＣjk´とする。前述と同様に、近傍粒子ｊの移動に伴うセルの移動量ベクトルをｄjkとすると、この移動量ベクトルｄjkは次のように表される。
【００４７】
ｄjk＝Ｃjk´−Ｃjk
ここで、全ての粒子に対して、あるｎ計算ステップ間の±ｐ，±ｑ方向の移動量がセル幅の半分より小さくなるようなｎを選ぶ。その理由を次に示す。すなわちあるｋ方向（例えばｐ方向）のセルのアドレスをＣik＝ｍ−１、Ｃjk＝ｍ＋１のように表すとき、粒子ｉ，ｊが図７のような位置関係にある場合は、粒子ｉ，ｊは接触候補粒子ペアではない。ここでＣik＝ｍ−１のセルとＣik＝ｍ＋１のセル内にある粒子ｉ，ｊの状況は、図７に示す状況もあれば、図８（Ａ）及び（Ｂ）に示す状況もある。明らかに図８（Ｂ）のケースと比較して図８（Ａ）のケースは、接触候補粒子ペアに近い状態にある。しかしながら図８（Ａ）の場合は、やはり接触候補粒子ペアではない。
【００４８】
接触候補粒子ペアは、実際に接触している粒子のペアと接触に近い状況にある粒子のペアを含んでいるわけであるが、その必要条件として、接触するペアは必ず接触候補粒子ペアに含まれていなければならない。これを逆に言うと、本当は接触しているのに、そのペアが接触候補粒子ペアにないという場合があってはならない。したがって、接触候補粒子ペアの更新タイミングに具備されるべき条件としては、接触候補粒子ペアが更新されるまでに新しい接触は一つもないことが要求される。
【００４９】
前述のように、図８（Ａ）に示す状態でも、粒子ｉ，ｊは接触候補粒子ペアではない。この粒子ｉ，ｊが最短経路で接触するためには、粒子ｉ，ｊが、それぞれお互いに向かってセル幅の半分よりも多く移動する場合である。したがって、接触候補粒子ペアの条件に最も近い２つの接触候補粒子ペアではない粒子が、接触しないための条件は、すべの粒子の移動量または変位量がセル幅の半分より小さいことである。そこで、本実施の形態では、接触候補粒子ペアの更新タイミングを次のように決める。
【００５０】
＊セル幅は想定する最大粒子の直径よりも大きいこと。
【００５１】
＊すべての粒子の中で移動量のｋ（ｐ，ｑまたはｐ，ｑ，ｒ）方向成分を調べ、その絶対値の最大値がセル幅の１／４より大きくなったとき、接触候補粒子ペアを更新する。
【００５２】
接触候補粒子ペアの更新タイミングがきて、前述のステップＳＴ２１での再度のセルマッピングが行われると、ステップＳＴ２２，ＳＴ２３及びＳＴ２５で、まず従前接触候補粒子ペアを構成していた２つの粒子が位置する２つのセルについてのそれぞれのアドレスの数値の差分の絶対値が１より大きいか否かを調べる。ステップＳＴ２２では、前回まで用いていた接触候補粒子ペア（ｉ，ｊ）について更新されたセルデータＣik´とＣjk´とを比較する。そしてステップＳＴ２３では、以下の動作が行われる。
【００５３】
まず、時刻Ｔ＝ｔ＋Δｔにおいて
｜Ｃik´−Ｃjk´｜≦１
ならばｉ，ｊは近傍関係にある。これは、図９に示した変位後のＣjk´の隣接するセルの範囲（ハッチングで示したセル領域）に粒子ｊが存在することを意味する。また、
｜Ｃik´−Ｃjk´｜＞１
ならばｉ，ｊは近傍関係にない。すなわち図９に示した変位後のＣjk´の隣接するセルの範囲（ハッチングで示したセル領域）に粒子ｊが存在しないことを意味する。
【００５４】
時刻Ｔ＝ｔ＋Δｔにおいて近傍関係にあるｉ，ｊのペアのうち、時刻Ｔ＝ｔにおいても近傍関係であったものは、接触候補粒子ペアデータを継続して用いることができるが（継続接触候補データ）、時刻Ｔ＝ｔでは近傍関係でなかったものについては、あらたに接触候補粒子ペアデータとして登録する必要がある（新規接触候補データ）。また、逆に、時刻Ｔ＝ｔにおいて接触候補粒子ペアであったものが、時刻Ｔ＝ｔ＋Δｔで接触候補粒子ペアでなくなった場合、データを消去する必要が生じる（接触候補データの解約）。
【００５５】
図１０に示すように、着目粒子ｉの移動に伴って、粒子ｉの近傍セルも移動する。例えば、上記の場合、粒子ｉが右斜め上のセル内に移動すると、粒子ｉを中心にした領域Ｒ１も領域Ｒ２へと移動することになる。
【００５６】
図１１に示すように、もともと着目粒子ｉの近傍セル内（領域Ｒ１内）にあった近傍粒子ｊが、着目粒子ｉの移動後に、新しい近傍セル内（Ｒ２領域内）に存在しない場合、着目粒子ｉと近傍粒子ｊの接触候補粒子ペアの関係は解約される。この関係は、
｜Ｃik−Ｃjk｜≦１  ａｎｄ  ｜Ｃik´−Ｃjk´｜＞１
が満たされるときに成り立ち、この着目粒子ｉと近傍粒子ｊのペアデータはメモリーから完全に消去する（ステップＳＴ２７）。
【００５７】
全ての接触候補粒子ペアについて、上記のチェックを完了したらステップＳＴ２８へと進む。図２のステップＳＴ２８乃至ＳＴ３７では、チェック後の接触候補粒子ペアについて新たなペアデータを作成しセルデータの入れ替えを実行する。ステップＳＴ２８では、更新されたセルデータから近傍関係にある接触候補粒子ペア（ｉ，ｊ）を抽出する。そしてステップＳＴ２９では、近傍関係にあるペア（ｉ，ｊ）について、更新前のセルデータを比較する。すなわち更新後における接触候補粒子ペアの更新前の関係を比較する。この比較を行うことにより、更新前も接触候補関係にあったのか（ステップＳＴ３０〜ＳＴ３２）、更新後に接触候補関係になったとのか（ステップＳＴ３４〜ＳＴ３６）を判断する。
【００５８】
図１２に示すように、もともと着目粒子ｉの近傍セル外（領域Ｒ１外）にあった近傍粒子ｊが、着目粒子ｉの移動後に、新しい近傍セル内（領域Ｒ２内）に入ってきた場合、新規の着目粒子ｉと近傍粒子ｊの接触候補粒子ペアの関係が発生する。この関係は、
｜Ｃik−Ｃjk｜＞１  ａｎｄ  ｜Ｃik´−Ｃjk´｜≦１
が満たされるときに成り立つ（ステップＳＴ３４及びＳＴ３５）。そして、この着目粒子ｉと近傍粒子ｊのペアデータをメモリー上に新たに作る（ステップＳＴ３６）。
【００５９】
また図１３に示すように、もともと着目粒子ｉの近傍セル内（領域Ｒ１内）にあった近傍粒子ｊが、着目粒子ｉの移動後も、新しい近傍セル内（Ｒ２領域内）に存在する場合には、着目粒子ｉと近傍粒子ｊの接触候補粒子ペアの関係は維持される。
【００６０】
この関係は、
｜Ｃik−Ｃjk｜≦１  ａｎｄ  ｜Ｃik´−Ｃjk´｜≦１
が満たされるときに成り立つ（ステップＳＴ３０及びＳＴ３１）。そして、この着目粒子ｉと近傍粒子ｊのペアデータはそのまま継続して利用する（ステップＳＴ３２）。
【００６１】
上記の動作をまとめると、まず従前接触候補粒子ペアを構成していた２つの粒子が位置する２つのセルについてのそれぞれのアドレスの数値の差分の絶対値が１より大きいか否かを調べて、接触候補か否かを調べることになる。すなわち２次元であれば、着目粒子ｉが存在するセルのアドレス（ｐi，ｑi)と近傍粒子ｊのアドレス（ｐj，ｑj)のそれぞれの数値の差分の絶対値｜ｐi−ｐｊ｜及び｜ｑi−ｑｊ｜をとり、そのいずれもが１より大きいものを接触候補粒子ペアではなくなった接触候補粒子ペアと決定する。また新たに接触候補決定ステップで決定した新たな接触候補粒子ペアについて、この新たな接触候補粒子ペアを構成する２つの粒子が従前に位置した２つのセルについてのそれぞれのアドレス（ｐi，ｑi)及び（ｐj，ｑj)の数値の差分の絶対値｜ｐi−ｐｊ｜及び｜ｑi−ｑｊ｜がいずれも１より大きいか否かを調べる。そして、絶対値が１より大きかったものを更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアと決定する。
【００６２】
接触候補粒子ペアデータそのものの中には、接触、非接触、接触時の情報を記述しておき、これをすべの接触候補粒子ペアについてメモリーに保存しておくと、データ量は膨大になる。そしてこの膨大な量のデータについての、継続、新規、解約が接触候補粒子ペアの更新時に発生する。そのため本実施の形態では、図１４に示すように、データ構造としては、挿入、削除が容易な接触候補粒子ペアのデータの保存に、各接触候補粒子ペアについてのデータが次の接触候補粒子ペアについてのデータへの参照を持つリンクリスト構造によるデータの保存方法を採用している。この方法では、各接触候補粒子ペアデータをそれぞれカードに記入しているものと考える。そして各カードは、各接触候補粒子ペアについてのデータが次の接触候補粒子ペアについてのデータへの参照を持つリンクリスト構造で作成されているものと考える。前述の接触候補粒子ペアの解除または解約を実行する際には、解除または解約の対象とする接触候補粒子ペアのデータを記録したカードを除く（削除する）ものと考える。また前述の新たに発生した接触候補粒子ペアについては、その接触候補粒子ペアのデータを記録したカードを作成して、新規に追加するものと考える。さらに継続する接触候補粒子ペアのデータについてのカードは、そのまま維持するものと考える。このようなリンクリスト構造により、データを保存すると、データの取り回し時間を短縮し、トータルした計算時間を大幅に縮めることができる。
【００６３】
上記のステップＳＴ２８乃至ＳＴ３３のチェックは、全ての接触候補粒子ペアについて行う。そして全ての接触候補粒子ペアについてのチェックが完了した段階で、セルデータの入れ替えを行う（ステップＳＴ３７）。
【００６４】
上記が、本実施の形態で採用している接触候補粒子ペアの更新方法である。結局、
ｄik ＝Ｃik´−Ｃik
ｄjk ＝Ｃjk´−Ｃjk
を求めるための各粒子のセル情報Ｃik（アドレス）をセルデータメモリに記憶させておくことで、接触候補粒子ペアの更新は非常に簡単に行うことができる。
【００６５】
本実施の形態の方法では、接触判定そのものを間引くのではなく、接触候補粒子ペアの更新を合理的に間引いて、しかも接触候補データの取り回しを効率よく行うことで、接触判定に必要な計算時間とメモリーの節約を行う。したがって本発明のプログラムでは、接触候補粒子ペア更新のタイミングを合理的に決めることで、計算精度を保ちつつ、接触候補粒子ペアの更新を間引くことで、計算時間を大幅に縮めることが可能になった。
【００６６】
上記実施の形態では、二次元で各接触候補粒子ペアの更新を説明したが、三次元であっても、上記説明が本質的に変わるところはない。
【００６７】
図１５は、本発明の方法を実施して、離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを演算により求める粒子データ演算装置の一例の構成を示すブロック図である。この装置は、セル分割手段１と、マッピング手段２と、接触候補粒子ペア決定手段３と、物理量演算手段４と、座標データ更新手段６と、変位量判定手段５と、接触候補粒子ペアデータ更新手段７と、座標データメモリＭ１と、セルデータメモリＭ２と、接触候補粒子ペアデータメモリＭ３と、物理量メモリＭ４とを有する。セル分割手段１は、座標データメモリＭ１に記憶された座標データによって特定される複数の粒子を含む領域を複数のセルに分割する。マッピング手段２は、領域分割手段１で得た複数のセルの中の１つのセルを基準にしたときに、他のセルが、基準にした１つのセルから２次元方向または３次元方向に何セル分離れているかが判るように複数のセルのそれぞれにアドレスを付す。また複数の粒子の座標データに基いて複数の粒子の中心が位置する複数のセルを特定する。そして複数の粒子の中心が位置する複数のセルを特定するデータをセルデータメモリＭ２に記憶させる。
【００６８】
接触候補粒子ペア決定手段３は、複数のセル内の各粒子に関して、複数のセルに付されたアドレスに基いて、着目した１つの粒子が位置するセルと隣接する他の複数のセル内に他の粒子の中心がある場合には、着目した１つの粒子と他の粒子とを接触候補粒子ペアであると決定する。接触候補粒子ペアデータメモリＭ３は、接触候補粒子ペア決定手段３で決定された接触候補粒子ペアについてのデータを記憶する。物理量演算手段４は、接触候補粒子ペアの粒子どうしが接触しているか否かを判定し、着目した１つの粒子の単位時間経過後の位置に基いて、着目した１つの粒子と接触している他の粒子及び着目した１つの粒子と接触する可能性の高い他の粒子との間の運動に関する物理量を演算する。物理量データメモリは、物理量演算手段４の演算結果を記憶する。
【００６９】
座標データ更新手段６は、物理量演算手段４の演算結果に基いて、複数の粒子の座標データを更新する。変位量判定手段５は、演算量演算手段４で演算した結果から接触候補粒子ペアを構成する複数の粒子の中で積算変位量が１／４セル以上になったものがあったか否かを判定する。そして接触候補粒子ペアデータ更新手段７は、変位量判定手段５が、１／４セル以上の積算変位量を判定したときには、再度マッピング手段２及び接触候補粒子ペア決定手段を動作させて更新を行う。具体的には、後従前は接触候補粒子ペアであったが更新後に接触候補粒子ペアではなくなった接触候補粒子ペアのデータを接触候補粒子ペアデータメモリＭ３から削除する。また従前接触候補粒子ペアではなかったが更新後に新たに接触候補粒子ペアとなった接触候補粒子ペアのデータを接触候補粒子ペアデータメモリＭ３に追加する。そして、残りの従前の接触候補粒子ペアのデータについてはそのまま保存を継続する。そして変位量判定手段５が１／４セル以上の積算変位量を判定しないときには、所定の時間間隔で物理量演算手段４、座標データ更新手段６及び変位量判定手段５は繰り返し動作をする。
【００７０】
上記の装置の主要部をコンピュータを用いて実現する場合には、図１及び図２に示したアルゴリズムのプログラムをコンピュータで実行すればよい。その場合、セル分割手段１は、図１のステップＳＴ１で実現される。またマッピング手段２は、図１のステップＳＴ２で実現される。さらに接触候補粒子ペア決定手段３は、図１のステップＳＴ３により実現される。さらに物理量演算手段４は、図１のステップＳＴ４〜ＳＴ１０によって実現され、座標データ更新手段６はステップＳＴ９によって実現され、変位量判定手段５はステップＳＴ１１によって実現され、接触候補粒子ペアデータ更新手段７はステップＳＴ１２によって実現される。
【図面の簡単な説明】
【００７１】
【図１】コンピュータにおいて実行される本発明の粒子データ演算用プログラムの一実施の形態のアルゴリズムをメモリと一緒に示したフローチャートである。
【図２】図１のフローチャートの一部のステップのアルゴリズムをより詳細に示したフローチャートである。
【図３】セルへの名称付けを説明するために用いる図である。
【図４】着目粒子が存在するセルに対する近傍セルを示す図である。
【図５】着目粒子が移動した場合のセルの名称変更を示す図である。
【図６】着目粒子と近傍粒子との関係を説明するための図である。
【図７】２つの粒子が接触候補粒子ペアになる状況の前提を説明するために用いる図である。
【図８】（Ａ）及び（Ｂ）は、セル内における粒子の位置の例を示す図である。
【図９】着目粒子が異なるセルに移動することにより、着目粒子の隣接するセルの範囲も移動することを説明するために用いる図である。
【図１０】着目粒子が移動したときの状況を説明するために用いる図である。
【図１１】着目粒子から近傍粒子が離れて、接触候補粒子ペアにならない状態を説明するために用いる図である。
【図１２】着目粒子に対して近傍粒子ではなかった粒子が、近傍粒子になる状態を説明するために用いる図である。
【図１３】着目粒子が移動しても、従前の近傍粒子が、再度近傍粒子になる状態を説明するために用いる図である。
【図１４】リンクリスト構造によりデータを保存する場合を概念的に説明するために用いる図である。
【図１５】本発明の方法を実施して、離散要素法に基いてコンピュータを用いてシミュレーションを行う際に使用する粒子データを演算により求める粒子データ演算装置の一例の構成を示すブロック図である。
【符号の説明】
【００７２】
Ｍ１  座標データメモリ
Ｍ２  セルデータメモリ
Ｍ３  接触後方ペアデータメモリ
Ｍ４  物理量データメモリ
１  セル分割手段
２  マッピング手段
３  接触候補粒子ペア決定手段
４  物理量演算手段
５  変位量判定手段
６  座標データ更新手段
７  接触候補粒子ペアデータ更新手段
【図１】
0004766383.tif 000002
【図２】
0004766383.tif 000003
【図３】
0004766383.tif 000004
【図４】
0004766383.tif 000005
【図５】
0004766383.tif 000006
【図６】
0004766383.tif 000007
【図７】
0004766383.tif 000008
【図８】
0004766383.tif 000009
【図９】
0004766383.tif 000010
【図１０】
0004766383.tif 000011
【図１１】
0004766383.tif 000012
【図１２】
0004766383.tif 000013
【図１３】
0004766383.tif 000014
【図１４】
0004766383.tif 000015
【図１５】
0004766383.tif 000016