"情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法。粒子iの辺及び粒子jの面を抽出するステップは、前記第1の平面に垂直で且つ粒子jの中心側を向く単位ベクトルn1を決定し、単位ベクトルn1と粒子iの各頂点の位置ベクトルPiの内積n1Piを計算し、単位ベクトルn1と粒子jの各頂点の位置ベクトルPjの内積n1Pjを計算し、粒子iの全ての頂点についての内積n1Piのうち、最大となるものをD1max、粒子jの全ての頂点についての内積n1Pjのうち、最小となるものをD1minとした場合に、内積n1PiがD1min以上D1max以下となる粒子iの頂点、及び、内積n1PjがD1min以上D1max以下となる粒子jの頂点からなる集合Aを定義し、集合Aに含まれる頂点を含む粒子iの辺及び粒子jの面を抽出するステップであることを特徴とする請求項1に記載の粒子挙動解析方法。情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側であって、且つ、前記第1の平面とは異なる角度で前記直線に交わり且つ粒子jに接する第3の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側であって、且つ、前記第3の平面に平行で且つ粒子iに接する第4の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法。粒子iの辺及び粒子jの面を抽出するステップは、前記第1の平面に垂直で且つ粒子jの中心側を向く単位ベクトルn1を決定し、単位ベクトルn1と粒子iの各頂点の位置ベクトルPiの内積n1Piを計算し、単位ベクトルn1と粒子jの各頂点の位置ベクトルPjの内積n1Pjを計算し、粒子iの全ての頂点についての内積n1Piのうち、最大となるものをD1max、粒子jの全ての頂点についての内積n1Pjのうち、最小となるものをD1minとした場合に、内積n1PiがD1min以上D1max以下となる粒子iの頂点、及び、内積n1PjがD1min以上D1max以下となる粒子jの頂点からなる集合Aを定義し、前記第3の平面に垂直で且つ粒子jの中心側を向く単位ベクトルn2を決定し、単位ベクトルn2と粒子iの各頂点の位置ベクトルPiの内積n2Piを計算し、単位ベクトルn2と粒子jの各頂点の位置ベクトルPjの内積n2Pjを計算し、粒子iの全ての頂点についての内積n2Piのうち、最大となるものをD2max、粒子jの全ての頂点についての内積n2Pjのうち、最小となるものをD2minとした場合に、内積n2PiがD2min以上D2max以下となる粒子iの頂点、及び、内積n2PjがD2min以上D2max以下となる粒子jの頂点からなる集合Bを定義し、集合Aと集合Bの積集合に含まれる頂点を含む粒子iの辺及び粒子jの面と、集合Aに含まれる頂点と集合Bに含まれる頂点の両方を含む粒子iの辺及び粒子jの面とを抽出するステップであることを特徴とする請求項3に記載の粒子挙動解析方法。前記単位ベクトルn1を、粒子iの中心と粒子jの中心を結ぶ直線に平行となるように決定することを特徴とする請求項2又は4に記載の粒子挙動解析方法。前記単位ベクトルn1を、粒子iの中心と粒子jの中心を結ぶ直線に平行となるように決定し、前記単位ベクトルn2を、D2max-D2minがD1max-D1minよりも小さくなるように決定することを特徴とする請求項4に記載の粒子挙動解析方法。密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置。前記接触判定部は、前記第1の平面に垂直で且つ粒子jの中心側を向く単位ベクトルn1を決定し、単位ベクトルn1と粒子iの各頂点の位置ベクトルPiの内積n1Piを計算し、単位ベクトルn1と粒子jの各頂点の位置ベクトルPjの内積n1Pjを計算し、粒子iの全ての頂点についての内積n1Piのうち、最大となるものをD1max、粒子jの全ての頂点についての内積n1Pjのうち、最小となるものをD1minとした場合に、内積n1PiがD1min以上D1max以下となる粒子iの頂点、及び、内積n1PjがD1min以上D1max以下となる粒子jの頂点からなる集合Aを定義し、集合Aに含まれる頂点を含む粒子iの辺及び粒子jの面を抽出することを特徴とする請求項7に記載の粒子挙動解析装置。密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側であって、且つ、前記第1の平面とは異なる角度で前記直線に交わり且つ粒子jに接する第3の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側であって、且つ、前記第3の平面に平行で且つ粒子iに接する第4の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置。前記接触判定部は、前記第1の平面に垂直で且つ粒子jの中心側を向く単位ベクトルn1を決定し、単位ベクトルn1と粒子iの各頂点の位置ベクトルPiの内積n1Piを計算し、単位ベクトルn1と粒子jの各頂点の位置ベクトルPjの内積n1Pjを計算し、粒子iの全ての頂点についての内積n1Piのうち、最大となるものをD1max、粒子jの全ての頂点についての内積n1Pjのうち、最小となるものをD1minとした場合に、内積n1PiがD1min以上D1max以下となる粒子iの頂点、及び、内積n1PjがD1min以上D1max以下となる粒子jの頂点からなる集合Aを定義し、前記第3の平面に垂直で且つ粒子jの中心側を向く単位ベクトルn2を決定し、単位ベクトルn2と粒子iの各頂点の位置ベクトルPiの内積n2Piを計算し、単位ベクトルn2と粒子jの各頂点の位置ベクトルPjの内積n2Pjを計算し、粒子iの全ての頂点についての内積n2Piのうち、最大となるものをD2max、粒子jの全ての頂点についての内積n2Pjのうち、最小となるものをD2minとした場合に、内積n2PiがD2min以上D2max以下となる粒子iの頂点、及び、内積n2PjがD2min以上D2max以下となる粒子jの頂点からなる集合Bを定義し、集合Aと集合Bの積集合に含まれる頂点を含む粒子iの辺及び粒子jの面と、集合Aに含まれる頂点と集合Bに含まれる頂点の両方を含む粒子iの辺及び粒子jの面とを抽出することを特徴とする請求項9に記載の粒子挙動解析装置。前記接触判定部は、前記単位ベクトルn1を、粒子iの中心と粒子jの中心を結ぶ直線に平行となるように決定することを特徴とする請求項8又は10に記載の粒子挙動解析装置。前記接触判定部は、前記単位ベクトルn1を、粒子iの中心と粒子jの中心を結ぶ直線に平行となるように決定し、前記単位ベクトルn2を、D2max-D2minがD1max-D1minよりも小さくなるように決定することを特徴とする請求項10に記載の粒子挙動解析装置。請求項1~6のうちいずれか1項に記載の粒子挙動解析方法の各ステップを情報処理装置に実行させることを特徴とする解析プログラム。本発明の目的は、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することにある。また本発明の目的は、密集状態にある多数の粒子の挙動を計算するに際し、粒子同士の接触判定を高速に行うことにある。上記課題を解決するために、本発明は以下の構成を採用する。本発明の第一態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。本発明の第二態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側であって、且つ、前記第1の平面とは異なる角度で前記直線に交わり且つ粒子jに接する第3の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側であって、且つ、前記第3の平面に平行で且つ粒子iに接する第4の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。本発明の第三態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。本発明の第四態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子iが粒子jに接触しているか否かを判定する場合に、粒子iの中心と粒子jの中心を結ぶ直線に交わり且つ粒子jに接する第1の平面よりも粒子jの中心側であって、且つ、前記第1の平面とは異なる角度で前記直線に交わり且つ粒子jに接する第3の平面よりも粒子jの中心側の空間に位置する粒子iの辺、及び、前記第1の平面に平行で且つ粒子iに接する第2の平面よりも粒子iの中心側であって、且つ、前記第3の平面に平行で且つ粒子iに接する第4の平面よりも粒子iの中心側の空間に位置する粒子jの面、を抽出し、前記抽出した粒子iの辺と粒子jの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。本発明の第五態様は、上述した粒子挙動解析方法の各ステップを情報処理装置コンピュータに実行させる解析プログラムである。本発明によれば、粒子同士の接触判定に要する時間を短縮でき、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することができる。本発明の実施形態は、密集状態にある多数の粒子の挙動をコンピュータシミュレーションにより計算するための粒子挙動解析方法及び装置を提供するものである。本実施形態に係る粒子挙動解析方法の特徴の一つは、多数の多角形から構成される多面体ポリゴンモデルにより各粒子の形状を定義した点である。これにより任意の粒子形状を表現できるようになり、多様なケースのシミュレーションを精度良く行うことが可能になる。本方法のもう一つの特徴は、粒子同士の接触判定を行う際に、全ての頂点、辺、面について網羅的に接触判定計算を行うのではなく、予め接触可能性のある辺及び面を特定し、それらのあいだについてのみ接触判定を行う点である。これにより、接触判定に要する時間を短縮でき、粒子の挙動計算の高速化を図ることができる。以下、本実施形態の具体的な構成及び特徴について図に基づき詳しく説明する。粒子挙動計算のモデル図1は、本実施形態における粒子挙動計算の力学モデルを説明するものであり、異形粒子非球形粒子が他の物体と接触した際に働く接触力を説明する図である。このモデルでは、粒子が他の物体と接触した際に発生するめり込みを、その接触面に対して法線方向と接線方向に分解して考える。粒子には、それぞれのめり込み量に対してスプリングで押し返す力が働く。また各スプリングと並列にダンパーを設定する。これはめり込み速度に比例して働く力であり、運動を減衰させるものである。またすべりを考慮するためのスライダを設定している。このように、スプリング、ダンパー、スライダを設定して衝突を考えることは従来の個別要素法と同じである。解析装置のハードウエア構成図2は、本実施形態における解析装置が動作する計算装置を示すブロック図である。情報処理装置コンピュータ10は、各種判断及び処理を行う中央処理装置CPU11と、処理プログラム及び固定データを格納するROM12と、処理データを記憶するRAM13と、入出力回路I/O14とから構成されている。また必要に応じて、情報処理装置10には、I/O14を介して記憶装置、入力装置、表示装置、他のコンピュータなどが接続される。I/O14は、これら外部装置と情報処理装置10との間でデータをやり取りするための回路である。この情報処理装置10には、入力データ15がI/O14を介して入力される。計算結果は、I/O14を介して出力データ16として出力される。解析装置の機能図3は、粒子挙動解析を行う解析プログラムによって実現される機能の構成を示すブロック図である。解析プログラムは、ROM12や記憶装置などのコンピュータ読み取り可能な記録媒体に記憶されており、必要に応じて、CPU11によって読み込まれ実行されるソフトウエアである。この解析プログラムは、主な機能モジュールとして、入力データ読込部取得部B110、接触判定部B120、粒子受力計算部B130、粒子運動計算部B140、計算結果出力部B150を備えている。全体制御部B100は、これらのモジュールの処理の順番やデータの引渡しなどの全体制御を担っている。入力データ読込部B110は、粒子形状データ、壁形状データ、及び各種パラメータを読み込んでRAM13に格納する。粒子形状データとは、粒子の表面を複数の多角形ポリゴンで表現したものである。ここで多角形としてはその頂点が一平面上にあることが望ましく、3角形を用いるのが最も単純である。また粒子は凸形状凸多面体に限定すると、他の粒子や壁との接触状態の判定が簡便になる。壁形状データも同様であり、壁の表面を複数の多角形で表現したものである。粒子形状データ、壁形状データともに、多角形の頂点座標と、各多角形を構成する頂点の番号とで構成されるデータである。各種パラメータには、粒子と壁の機械特性ヤング率、ポアソン比、減衰係数、摩擦係数、比重、粒子に働く力重力、付着力、電磁力、空気の粘性特性等、及び計算条件としての計算刻み時間Δt、計算終了時刻がある。なお、壁とは、粒子が移動可能な空間を形成する規定する境界面である。接触判定部B120は、粒子の位置情報各頂点の座標、粒子の中心座標に基づいて、粒子が他の物体周囲の粒子または壁と接触しているかどうか、またどのように接触しているかを判定する。粒子受力計算部B130は、各粒子が他の物体から受ける力を計算する機能である。粒子は他の粒子または壁と接触している場合に、そのめり込み量に応じて接触抗力を受ける。粒子運動計算部B140は、粒子受力計算部B130で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。計算結果出力部B150は、得られた各粒子の位置と速度の結果を出力データ16として出力するものである。以下に、接触判定部B120、粒子受力計算部B130、及び粒子運動計算部B140の詳しい内容と計算方法を示す。接触判定部B120まずは、接触判定部B120における粒子同士の接触判定の内容を図4を用いて説明する。図4の上段は4種類の接触パターンを3次元的に示したものであり、下段は各パターンにおける2つの粒子の接触によるめり込みの様子を示すものである。図4において、実線は6面体粒子i、破線は6面体粒子jを示す。i-v0は粒子iを構成する多角形の頂点、i-e0は粒子iを構成する多角形の辺、j-e0とj-e1は粒子jを構成する多角形の辺、j-f0は粒子jを構成する多角形の面、星印は粒子iと粒子jの接触点を示す。また図4の上段における矢印は粒子iが粒子jに衝突する方向を示し、また下段における黒丸は粒子jの辺を示す。粒子が凸形状である場合、粒子iと粒子jの接触は、図4に示す次のa~dの4つのパターンのいずれかとなる。なお、図4では、粒子同士の接触判定を例に挙げているが、粒子と壁の接触判定も同じように考えることができる。a頂点のめり込み:粒子iを構成する多角形の頂点が粒子jの1つの面にめり込む場合である。粒子iの頂点が粒子jの内部にあれば、接触していると判定する。図4aにおいては、粒子iの頂点i-v0が粒子jの面j-f0にめり込んでいる。このパターンでは、頂点i-v0を、粒子iの粒子jに対する接触点として考える。b辺のめり込みケース1:上記aに加えて、めり込んだ頂点を端点に持つ粒子iの辺が粒子jの辺にもめり込む場合である。aのめり込んだ頂点を端点に持つ粒子iの辺において、粒子j構成面との交点数が1であり、かつその交点を持つ面がaにおけるめり込み面でなければ、これに該当する。図4bにおいては、粒子iの頂点i-v0が粒子jの面j-f0にめり込む上記aに該当とともに、粒子iの辺i-e0が粒子jの辺j-e0にもめり込んでいる。このパターンでは、頂点i-v0に加え、辺j-e0に最も近い辺i-e0上の点も、粒子iの粒子jに対する接触点として考える。c辺のめり込みケース2:粒子iの辺が粒子jの1つの辺だけにめり込む場合である。粒子iの辺i-e0において、粒子j構成面との交点数が2であり、その交点を持つ2つの面が共有辺を持ち、かつその共有辺が粒子iの辺に近接していれば、これに該当する。図4cにおいては、粒子iの辺i-e0が粒子jの辺j-e0にめり込んでいる。このパターンでは、辺j-e0に最も近い辺i-e0上の点を、粒子iの粒子jに対する接触点として考える。d辺のめり込みケース3:粒子iの辺が粒子jの2つの辺にめり込む場合である。粒子iの辺i-e0において、粒子j構成面との交点数が2であり、cではない場合がこれに該当する。図4dにおいては、粒子iの辺i-e0が粒子jの辺j-e0とj-e1にめり込んでいる。このパターンでは、辺j-e0に最も近い辺i-e0上の点、及び、辺j-e1に最も近い辺i-e0上の点を、粒子iの粒子jに対する接触点として考える。上記の接触判定を行うには、図4aの場合には、粒子iの頂点と粒子jの面の位置関係、図4b~dにおいては、粒子iの辺と粒子jの面の交差を調べる必要がある。辺と面のめり込み量の計算を行うための接触判定の仕方を6面体の粒子を例に、図5を用いて説明する。図5の20は粒子jを構成する1つの多角形の面を表し、この法線ベクトルをnとする。21は粒子iを構成する辺の1つi-e0を表している。この辺21の両端の頂点をそれぞれP0px0,py0,pz0,P1px1,py1,pz1とする。また、頂点j-v0,j-v1,j-v2,j-v3を含む平面無限と辺21の交差点を点Qとする。このとき、点Qが多角形面20上にあるか否かの判定を行う。まず、粒子iの辺21が粒子jの多角形面20の4つの頂点j-v0,j-v1,j-v2,j-v3を含む平面無限を横切っているかどうかを調べる。粒子jの多角形面20上の一つの任意点A0と粒子iの辺21の頂点P0,P1をそれぞれ結び、ベクトルP0A0,P1A0とする。ベクトルP0A0は、点P0から点A0に向かうベクトルであり、ベクトルP1A0は、点P1から点A0に向かうベクトルである。このとき、2つのベクトルP0A0,P1A0と法線ベクトルnの内積L0=P0A0nとL1=P1A0nが同じ符号ならば、粒子iの辺21が粒子jの多角形面20の頂点を含む平面無限を横切っていないことがわかる。逆に、内積L0とL1の符号が異なる場合は、粒子iの辺21が粒子jの多角形面20の頂点を含む平面無限を横切っていることがわかる。次に、粒子iの辺21が粒子jの多角形面20の頂点を含む平面無限を横切っている場合、粒子jの多角形面20の頂点を含む平面無限と粒子iの辺21との交点Qの座標qx,qy,qzを下記式で求める。ここで、rは内分比であり、r=L0/L0-L1で求まる。次に、交点Qと頂点j-v0,j-v1,j-v2,j-v3を結ぶ4つのベクトルQv0,Qv1,Qv2,Qv3について、隣のベクトルとの外積を下記式を用いて求め、ベクトルV0,V1,V2,V3を得る。ベクトルV0は、ベクトルQv0とQv1の両方に垂直なベクトル、すなわち多角形面20に対し垂直なベクトルとなり、その向きはベクトルQv0とQv1のなす角が0~πかπ~2πかで反転する。他のベクトルV1~V3も同様の性質をもつ。交点Qが多角形面20の内部にある場合は、Qv0とQv1、Qv1とQv2、Qv2とQv3、Qv3とQv0のいずれのなす角も0~πとなるから、ベクトルV0~V3の向きは一致する。一方、交点Qが多角形面20の外部にある場合は、なす角が0~πとなる組と、π~2πとなる組とがでるため、ベクトルV0~V3の向きは一致しなくなる。このような性質を利用することで、交点Qが多角形面20の内側に存在するか否か、つまり粒子iの辺21と粒子jの多角形面20とが交差するか否かを判定することができる。ここでは、下記式のように内積M0~M3を計算し、それらの符号を調べることで判定を行う。M0~M3の符号が全て同じ場合は、粒子jの多角形面20と粒子iの辺21は交差しており、一つでも符号が逆のものがあれば交差していないことがわかる。上記のような接触判定では、1つの四角形面に対して1つの辺が交差していることを判定するのに、43個の加減算、47個の乗算、5回の条件演算を要する。そのため、粒子の多面体の面数が多いほど、また粒子の数が多いほど、膨大な計算時間がかかってしまう。そこで、本実施形態では接触判定の前に接触可能頂点を絞り込み、その接触可能頂点を含む辺、面のみを接触の可能性がある辺、面として接触判定を行う。以下、接触の可能性がある辺及び面を抽出する方法として、2つの方法を説明する。第1の抽出方法図6は、辺及び面の第1の抽出方法について説明する図である。粒子i、jは3次元形状であるが、簡単のため2次元で表示している。第1の抽出方法は、粒子iの粒子jに対する接触を判定する際に、粒子iの中心と粒子jの中心を結ぶ直線に交わる2つの平面61,62を設定し、平面61,62で挟まれた空間に位置する粒子iの辺及び粒子jの面のみを判定対象として抽出するものである。ここで、平面61第1の平面と平面62第2の平面はそれぞれ粒子j,iに接し、且つ、互いに平行になるように設定される。言い換えると、第1の抽出方法は、平面61よりも粒子jの中心側の空間に位置する粒子iの辺と、平面62よりも粒子iの中心側の空間に位置する粒子jの面とを、判定対象として抽出する方法である。なお、辺又は面が空間に位置するとは、辺又は面の少なくとも一部分が空間内に存在すること空間と交わりを有することを意味する。第1の抽出方法を実現する具体的なアルゴリズムとして、本実施形態の接触判定部B120では次のような計算を行う。粒子iの中心から粒子jの中心に向かう中心結線ベクトルを考え、原点を通り中心結線ベクトルに垂直な面を基準面Sとする。粒子iの各頂点と基準面Sとの距離D1を計算し、最大値をD1maxとする。粒子jの各頂点と基準面Sとの距離を計算し、最小値をD1minとする。基準面Sとの距離がD1min以上D1max以下となる粒子i、jの各頂点を接触可能頂点として選ぶ。図6で塗りつぶし黒丸で示されている頂点が接触可能頂点である。そして、接触可能頂点を含む辺及び面が判定対象として抽出される。図6で実線で示されている辺が抽出された辺である。D1max、D1minの算出方法をより詳しく説明する。中心結線ベクトルと向きが同じで原点を起点とする単位ベクトルをn1=n1x,n1y,n1zとし、粒子iの各頂点の位置ベクトルをPik=xik,yik,zikとする。n1のことを中心結線単位ベクトルとも称す。ここで、kは粒子iの各構成頂点の番号である。基準面Sと粒子iの各頂点との距離D1は、下記式のように、単位ベクトルn1と頂点の位置ベクトルPikの内積n1Pikで算出される。同様に、粒子jの各頂点の位置ベクトルをPjk=xjk,yjk,zjkとすると、基準面Sと粒子jの各頂点との距離D1は下記式により算出される。粒子i及び粒子jの全ての頂点について距離D1を求めた後、D1max及びD1minは下記式で算出される。なお、説明を簡単にするためD1を距離と表現したが、実際にはD1は内積n1Pik、内積n1Pjkのようにベクトルの内積で定義されるものである。よって、n1とPik又はPjkの位置関係によっては、各頂点に対するD1の一部またはすべてが負になることもあるが、その場合でも上記式によりD1max、D1minを求めることができる。また、ここでは簡単のためn1の起点を原点としたが、ベクトルの平行移動に対して内積の大小は変わらないため、起点をどの位置にとってもよい。また、ここでは単位ベクトルn1を中心結線ベクトルに平行にとったが、ベクトルn1の向きはこれに限られない。ベクトルn1が中心結線ベクトルに対して直交していなければつまり、平面61,62が粒子i,jの中心を結んだ直線に交わっていれば、同様の処理が可能である。以上述べた方法により接触可能頂点を選択し、判定対象とする粒子iの辺と粒子jの面を絞り込むことで、接触判定を行う辺、面の組み合わせを減らし、計算時間を短縮することができる。第2の抽出方法上述した第1の抽出方法では、接触判定を行う範囲を平面61と62で挟まれた空間に限定している。しかし、図6で示されるように、平面61と62のあいだに、実際には他の粒子と接触していない頂点が多く含まれる場合もある。そこで、第2の抽出方法として、接触可能頂点をさらに絞り込む方法について説明する。第2の抽出方法は、図7に示すように、平面61,62とは異なる角度で平面71,72を設定し、4つの平面61,62,71,72で挟まれた空間に位置する粒子iの辺及び粒子jの面のみを判定対象として抽出するものである。ここで、平面71第3の平面と平面72第4の平面はそれぞれ粒子j,iに接し、且つ、互いに平行になるように設定される。言い換えると、第2の抽出方法は、平面61及び平面71よりも粒子jの中心側の空間に位置する粒子iの辺と、平面62及び平面72よりも粒子iの中心側の空間に位置する粒子jの面とを、判定対象として抽出する方法である。第2の抽出方法を実現する具体的なアルゴリズムを説明する。図7は2つのベクトルによる接触可能頂点の絞り込み方法を示す図である。図7左上は上記の中心結線単位ベクトルn1によって絞り込まれた接触可能頂点を示している。塗りつぶし黒丸で示されているのが接触可能頂点である。図7左下は、重なり最小単位ベクトルn2によって絞り込まれた接触可能頂点を表している。重なり最小単位ベクトルとは、粒子iと粒子jの重なりを最小とする向きの単位ベクトルをいい、以降の説明では単に重なり最小ベクトルとも称す。n2による接触可能頂点の絞り込み方法はn1による接触可能頂点の絞り込み方法と同じである。つまり、ベクトルn2と粒子iの各頂点座標の位置ベクトルPikの内積n2Pikの最大値D2maxと、ベクトルn2と粒子jの各頂点座標の位置ベクトルPjkの内積n2Pjkの最小値D2minを求める。そして、内積がD2min以上D2max以下となる頂点を接触可能頂点として抽出する。図7右は中心結線ベクトルによる接触可能頂点集合をA、重なり最小ベクトルによる接触可能頂点集合をBとしたときに、共通部分である積集合A∩Bで選択される頂点を表している。図7において、集合Aに含まれる頂点数は10、集合A内の頂点を構成頂点とする辺数は12、集合Bに含まれる頂点数は8、辺数は10である。一方、積集合A∩Bに含まれる頂点数は4、辺数は6となり、1つのベクトルのみの場合より判定対象を絞り込むことができる。よって、積集合A∩Bに含まれる頂点を最終的な接触可能頂点と定義する。積集合A∩Bが空集合でも接触する場合がある。図8を用いて説明する。四角で表されている頂点が中心結線ベクトルにより選択された頂点集合A、三角で表されているのが重なり最小ベクトルにより選択された頂点集合Bである。塗りつぶしの円で表されているのが積集合A∩Bに含まれる接触可能頂点である。粒子jは接触可能頂点を持たないが、図8に示すように粒子iと接触している。このような場合を含めるため、判定対象とする辺、面は次のように決定する。積集合A∩Bに含まれる頂点接触可能頂点を持つ辺、面。集合Aに含まれる頂点と集合Bに含まれる頂点の両方を持つ辺、面。以上の2パターンによって、全ての接触の可能性がある辺と面を選択することができる。本方法により、接触判定を行う辺、面の組み合わせを大幅に減らすことができ、計算時間のより一層の短縮が可能となる。本方法で用いる重なり最小ベクトルn2の算出方法について説明する。図9は重なり最小ベクトルを求めるフローチャートである。図9に示す処理は、接触判定部B120が実行するものである。まず、重なり最小ベクトルの初期値を設定するステップS1。初期値は中心結線単位ベクトルと平行な単位ベクトルとする。粒子iの中心位置ベクトルをCi、粒子jの中心位置ベクトルをCjとすると、重なり最小ベクトルn2の初期値は下記式で表される。以降の処理で、ベクトルn2の向きを微小に変動させながら、粒子iと粒子jの重なりが最小となるものを見つける。その処理で用いる微小量δの初期値δmax=0.1程度をステップS1で設定する。次に、2粒子の重なり量L=D2max-D2minを算出するステップS2。次に、接触の可能性があるか否かを判定するステップS3。Lが0以下の場合、2粒子は非接触であるので、接触判定から除外されるステップS4。Lが0より大きい場合は接触の可能性があるため、重なり最小ベクトルn2を更新探索するためのループに入る。重なり最小ベクトルn2に直交し、たがいに直交するベクトルu、vを適当に定めるステップS5。そして、試行する重なり最小ベクトルn21~n24を下記式のように設定するステップS6。上記式の1行目で決定される試行重なり最小ベクトルn21は、図9中に示されるように、直交ベクトルuの方向へ微小量δで決まる量だけ回転させたものとなる。このようにして決定された各試行重なり最小ベクトルn2kに対して2粒子の重なり量Lnewkを新たに算出するk=1,2,3,4。算出された4つのLnewkの最小値とLを比較するステップS7。min{Lnewk}がLより小さければ、重なり最小ベクトルn2をn2kに更新する。また、LもLnewkに更新するステップS8。重なり量Lを最小とするループにおいて、min{Lnewk}がLより大きくなった場合はステップS9、微小量δをδ/2と小さくした後、試行重なりベクトルの設定から繰り返すステップS10。上記のループを繰り返して、δがδmin=0.001程度より小さくなったときのn2が最終的な重なり最小ベクトルとなるステップS11。このような処理により、粒子iとjの重なり量Lを最も小さくするベクトルn2を得ることができる。なお、図9のステップS5で示したように、ループを繰り返すたびにu,vを面内で回転させると、局所解への陥りを抑制し、重なり最小ベクトルの最適解を得やすいという効果が得られる。図9で示した方法は重なり最小ベクトルn2を決定するための一手法にすぎない。他の公知の最適化手法を利用して、重なり量Lを最小とするベクトルn2を求めることも可能である。また、少なくともベクトルn1とn2の向きが異なっていれば、第1の抽出方法よりも判定対象を絞り込む効果は得られるため、ベクトルn1に対して任意の回転を行ったものをベクトルn2として用いるだけでもよい。接触判定処理接触判定部B120の処理フローを図10で説明する。まず、接触判定を行う2粒子i、jを設定するステップS100。次に2粒子i、jの中心座標から中心結線単位ベクトルn1を求めるステップS110。n1を基に2粒子i、jの重なり量Lを算出するステップS120。重なり量Lを基に接触可能性を判定しステップS130、L>0でなければ、2粒子は非接触として接触判定を終了するステップS135。次に、2粒子i、jの位置座標情報から重なり最小単位ベクトルn2を求めるステップS140。ステップS140の処理内容は、図9で詳しく説明したとおりである。ここで、ステップS110で設定したn1を用いて接触可能頂点集合Aを設定しステップS150、ステップS140で設定したn2を用いて接触可能頂点集合Bを設定するステップS160。次に、接触可能頂点集合A∩Bを設定するステップS170。こうして設定した接触可能頂点を含む粒子iの辺と粒子jの面を抽出するステップS180。以上で、接触判定に用いる判定対象が絞り込まれる。次に、粒子iが粒子jに接触しているか否かの判定処理を行う。まず、接触判定を行う粒子iの辺を設定するステップS190。接触判定を行う粒子jの面を設定するステップS200。ここで設定した粒子iの辺と粒子jの面との交差の有無を判定するステップS210。ここで、粒子jの面に粒子iの辺が交差していれば、交差回数のカウントアップを行うステップS220。この操作をステップS180で抽出された粒子jの面の全てに対して行うステップS230。次に、粒子iの辺と粒子jの面の交差回数が2回かどうかの判定を行うステップS240。交差回数が2回ならば、めり込み量を算出し、メモリに記憶しておくステップS250。本操作をS180で抽出された粒子iの辺の全てに対して行うステップS260。なお、図10では、説明を簡単にするため、粒子iのある辺が粒子jの2つの面と交差した場合に、接触有りと判定し、粒子iとjのあいだのめり込み量を算出する処理例を示した。しかし実際の処理では、図4に示したように、粒子iと粒子jの接触を4つのパターンで判定する。前述したようにパターンごとに判定条件が異なるが、いずれのパターンにおいても粒子iの辺と粒子jの面との交差判定が基本となるので、図10のステップS180で抽出した辺及び面を利用すればよい。各パターンa~dにおけるめり込み量の定義について、図4中の下段の図を用いて説明する。a頂点のめり込み:粒子iの頂点i-v0と粒子jの面j-f0との距離をめり込み量とする。b辺のめり込みケース1:粒子iの辺i-e0と粒子jの辺j-e0との距離をめり込み量とする。c辺のめり込みケース2:粒子iの辺i-e0と粒子jの辺j-e0との距離をめり込み量とする。d辺のめり込みケース3:粒子iの辺i-e0に対する、粒子jの辺j-e0、j-e1との距離をそれぞれの辺とのめり込み量とする。粒子受力計算部B130以上の4パターンにおいて、めり込み量距離の計算に用いた各粒子表面の2点を結ぶ方向を接触面法線方向と定める。この接触面法線方向は、接触点の粒子jの面又は辺に対するめり込み方向ということもできる。そして、接触判定部B120で求めためり込み量に基づいて、ヘルツの式を用いて、粒子iに働く接触面法線方向の接触力バネ力を求める。また各接触点に対して、接触が始まってからの接触面上記接触面法線方向に垂直な面内方向の移動量を求め、ミンドリンの式を用いて粒子iに働く接触面内方向の接触力バネ力を求める。また接触面法線方向の運動速度めり込み速度と面内方向の各運動速度ずり速度から、各方向に作用する粘性力ダンパー力を求める。そしてそれら接触力と粘性力の和に、粒子iに働くその他の力重力、付着力、電磁力、空気の粘性抵抗等を加えて、各粒子iが受ける並進力と回転力を求める。なお、本めり込み量から接触力を求める過程は、上記非特許文献1に記載されている手法と同じであるため、詳しい説明は割愛する。粒子運動計算部B140次に、粒子運動計算部B140における各粒子の刻み時間後の速度と位置の算出方法を説明する。この処理は、粒子受力計算部B130で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。すなわち、粒子運動計算部B140は、並進力からニュートンの運動方程式を用いて並進方向の加速度を求め、更に速度と位置を求める。また粒子運動計算部B140は、回転力からオイラーの運動方程式を用いて角加速度を求め、更に回転速度と回転位置を求める。本処理は、非特許文献1に記載されている手法と同じであるため、詳しい説明は割愛する。実施例以上の接触判定方法を用いて粒子挙動計算を実施した例を、図11~図12を用いて説明する。図11aは解析モデルを示す。解析モデルは電子写真技術を用いた画像形成装置におけるクリーニングプロセスである。0.1m/secで矢印の方向に移動する感光体ドラムに対してクリーニングブレードをβ=30°で当接させておき、感光体ドラム上に乗って運ばれてきたトナーの運動を見る。各トナーの形状を粒子形状データとして与え、感光体ドラム及びクリーニングブレードの形状を壁形状データとして与える。図11bは、解析条件として与えた、トナー、感光体ドラム、クリーニングブレードそれぞれの機械特性ヤング率、ポアソン比、減衰係数、摩擦係数、比重を示している。この条件で、刻み時間Δtを10nsecに設定し、約2百万ステップの計算を実行した。図12は、a球形粒子とb異形粒子の計算結果を比較したものである。球形粒子の計算は従来の計算方法で実施している。異形粒子は4面体の角をとった形状であり、その表面を8つの3角形と6つの4角形の合計12の多角形で構成している。なお本計算においては、トナーがブレードに堰き止められた際の挙動をみるため、強制的にトナーがすり抜けないようにしている。球形粒子に対して異形粒子は回転しにくいことから、全体的に動き辛く、トナー粒子の密集度が上がる傾向にある。図12の計算結果では、異形粒子の密集度が高くなっている様子が再現されている。このトナーがブレードに堰き止められた際の挙動は、クリーニング不良であるトナーのすり抜けとの関係があり、流動的である程すり抜けやすいといわれている。このことから、製品設計を行う上で、トナーの形状によるクリーニング不良の傾向を粒子挙動計算で予測することが可能となる。また、図13は面数の異なる異形粒子の自由落下の計算を示した図である。aは、頂点数12、面数14、辺数24の場合の異形粒子で、bは頂点数98、面数108、辺数204の異形粒子の場合である。aとbの場合において、256個の異形粒子の自由落下の計算を刻み時間Δt=10nsecで8百万ステップ実行し、本発明を使用しない場合と使用した場合での計算時間の比較を行った。aでは、従来の方法では計算に19時間35分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込み第1の抽出方法を用いた場合では、8時間09分となり、高速化率は2.4倍となる。bでは、従来の方法では計算に524時間32分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込みを用いると22時間24分となり、高速化率は23.4倍となった。このことから、粒子の面の数が多くなればなるほど、本発明を用いることによる計算時間の短縮効果が大きくなることがわかる。図14は中心結線ベクトルと重なり最小ベクトルを用いた高速化の効果を示す図である。計算には頂点数198、面数392、辺数588の異形粒子を用い、図13と同様の自由落下計算を300000ステップ、Δt=10nsecで行った。接触判定に中心結線ベクトルのみを用いた場合第1の抽出方法と中心結線ベクトルと重なり最小ベクトルの両方を用いた場合第2の抽出方法で計算時間を比較した。中心結線ベクトルのみを用いた場合は4726秒、中心結線ベクトルと重なり最小ベクトルの両方を用いた場合は2009秒となり、高速化率は2.35倍となった。このことから、中心結線ベクトルと重なり最小ベクトルの両方を用いた接触判定により、さらなる計算時間の短縮が実現できていることが分かる。"