粒子挙動解析方法、粒子挙動解析装置、及び解析プログラム
【課題】密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算する。
【解決手段】密集状態にある多数の粒子の挙動解析において、粒子ｉと粒子ｊの間の接触を判定する場合、粒子ｉから粒子ｊに向かう向きの単位ベクトルと粒子ｉ、ｊの重なりを最小とする向きの単位ベクトルを定義し、両単位ベクトルと２粒子の各頂点の位置ベクトルの内積から接触可能頂点を絞り込むことで、接触判定を行う辺、面の組み合わせを大幅に減らす。
【解決手段】密集状態にある多数の粒子の挙動解析において、粒子ｉと粒子ｊの間の接触を判定する場合、粒子ｉから粒子ｊに向かう向きの単位ベクトルと粒子ｉ、ｊの重なりを最小とする向きの単位ベクトルを定義し、両単位ベクトルと２粒子の各頂点の位置ベクトルの内積から接触可能頂点を絞り込むことで、接触判定を行う辺、面の組み合わせを大幅に減らす。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、密集状態にある多数の粒子について、その形状を考慮して挙動を解析する技術に関するものであり、特に挙動解析計算における粒子同士の接触判定に関する。
【背景技術】
【０００２】
粉体や粒体などの粒子を取り扱う分野では、個々の粒子の挙動を知ることが重要である。従来これらの分野では、実験と観察から粒子の挙動を把握することが多かった。しかし近年、計算機内でその運動を再現させ、その挙動を把握する、いわゆるシミュレーション技術の適用が活発に行われている。正確な粒子挙動のシミュレーションを行うには、粒子のもつ形状、機械的特性等の情報をできるだけ詳細に計算モデルに反映させ、それらの因子を考慮して解を求めることが重要である。また、土木分野における砕石、砂利等の比較的大きな物体も「粒子」と考えることにより、上記シミュレーション技術をこれらの物体の挙動を計算する場合へも適用可能である。
【０００３】
一方、このような粒子を扱う装置の１つに、プリンタ、複写機、ファクシミリ等の電子写真技術を用いた画像形成装置がある。これらの装置では、直径数ミクロンの粉体粒子からなるトナーを用いて、帯電、露光、現像、転写、定着、クリーニングという６つのプロセスを経て、紙上に画像を形成する。このような画像形成装置では、トナー粒子の形状が、作成される画像に大きな影響を与えることがわかっている。そこでこのような画像形成装置の開発においても、シミュレーション技術を用いて、密集状態にある個々のトナー粒子の挙動を求める試みが活発に行われ、近年の傾向として特にトナー粒子の形状を考慮することが求められている。
【０００４】
粒子が密集した状態での各粒子の挙動をシミュレーションする代表的な方法として個別要素法（ＤＥＭ）がある。個別要素法の技術は非特許文献１に詳細に記載されているが、個々の粒子を球形で表現し、それらの接触に際してはそのめり込み量からヘルツの式とミンドリンの式を用いることで粒子の持つ機械特性を考慮した接触抗力を求める。そしてニュートンの運動方程式（並進方向）とオイラーの運動方程式（回転方向）を用いて、その接触抗力をもとに個々の粒子の時間進展に伴う運動を求めるものである。
【０００５】
個別要素法の特徴は個々の粒子を要素と捉え、それらが接触する部分にバネ、ダンパー、スライダを仮定して、その接触力を求め、運動方程式を解くところにある。この考え方を球形以外の形状をもつ粒子の挙動計算に適用したものとして、非特許文献２、非特許文献３がある。これらの文献では粒子の形状を多面体で表現し上述の個別要素法を用いてその挙動を計算している。そして、粒子の接触の検出に関して、Common-planeと呼ばれる面を想定して、接触の有無、接触力の作用点、接触力の作用方向、及び接触力の大きさを決定する。しかしながら、同文献の方法ではCommon-planeと粒子の接触断面積を求めることが必要となるが、断面積を求めるには粒子の各頂点、辺、面とCommon-planeとの接触状態を調べなければならないため、長い計算時間を要する。また本発明者らが同文献の方法を実際に試行しようとしたところ、同文献は接触力を作用させる位置についての情報開示が乏しく、再現することが困難であることが分かった。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】粉体工学会編、粉体シミュレーション入門(1998)
【非特許文献２】P.A. Cundall: Formulation of a three-dimensional distinct element model - Part I. a scheme to detect and represent contacts in a system composed of many polyhedral blocks, Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 25(3), pp.107-116 (1988)
【非特許文献３】R. Hart, P.A. Cundall, and J. Lemos: Formulation of a three-dimensional distinct element model - Part II. mechanical calculations for motion and interaction of a system composed of many polyhedral blocks, Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 25(3), pp.117-125 (1988)
【発明の概要】
【発明が解決しようとする課題】
【０００７】
本発明の目的は、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することにある。また本発明の目的は、密集状態にある多数の粒子の挙動を計算するに際し、粒子同士の接触判定を高速に行うことにある。
【課題を解決するための手段】
【０００８】
上記課題を解決するために、本発明は以下の構成を採用する。
本発明の第一態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。
【０００９】
本発明の第二態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。
【００１０】
本発明の第三態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計
算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。
【００１１】
本発明の第四態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。
【００１２】
本発明の第五態様は、上述した粒子挙動解析方法の各ステップを情報処理装置（コンピュータ）に実行させる解析プログラムである。
【発明の効果】
【００１３】
本発明によれば、粒子同士の接触判定に要する時間を短縮でき、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することができる。
【図面の簡単な説明】
【００１４】
【図１】粒子挙動計算の力学モデルを説明する図。
【図２】解析装置のハードウエア構成を示す図。
【図３】解析装置の機能構成を示す図。
【図４】接触パターンと各パターンにおける粒子めり込み量を示す図。
【図５】辺と面の交差の有無の判定方法を示す図。
【図６】接触可能頂点の絞り込み方法を示す図。
【図７】２ベクトルによる接触可能頂点の絞り込み方法を示す図。
【図８】接触可能頂点を持たない粒子との接触を示す図。
【図９】重なり最小ベクトルを求める方法を示す図。
【図１０】接触判定部の処理フローを示す図。
【図１１】実施例の計算条件を示す図。
【図１２】実施例の計算結果を示す図。
【図１３】中心結線ベクトルを用いた高速化の効果を示す図。
【図１４】中心結線ベクトルと重なり最小ベクトルを用いた高速化の効果を示す図。
【発明を実施するための形態】
【００１５】
本発明の実施形態は、密集状態にある多数の粒子の挙動をコンピュータ・シミュレーションにより計算するための粒子挙動解析方法及び装置を提供するものである。本実施形態に係る粒子挙動解析方法の特徴の一つは、多数の多角形から構成される多面体（ポリゴンモデル）により各粒子の形状を定義した点である。これにより任意の粒子形状を表現でき
るようになり、多様なケースのシミュレーションを精度良く行うことが可能になる。本方法のもう一つの特徴は、粒子同士の接触判定を行う際に、全ての頂点、辺、面について網羅的に接触判定計算を行うのではなく、予め接触可能性のある辺及び面を特定し、それらのあいだについてのみ接触判定を行う点である。これにより、接触判定に要する時間を短縮でき、粒子の挙動計算の高速化を図ることができる。以下、本実施形態の具体的な構成及び特徴について図に基づき詳しく説明する。
【００１６】
（粒子挙動計算のモデル）
図１は、本実施形態における粒子挙動計算の力学モデルを説明するものであり、異形粒子（非球形粒子）が他の物体と接触した際に働く接触力を説明する図である。このモデルでは、粒子が他の物体と接触した際に発生するめり込みを、その接触面に対して法線方向と接線方向に分解して考える。粒子には、それぞれのめり込み量に対してスプリングで押し返す力が働く。また各スプリングと並列にダンパーを設定する。これはめり込み速度に比例して働く力であり、運動を減衰させるものである。またすべりを考慮するためのスライダを設定している。このように、スプリング、ダンパー、スライダを設定して衝突を考えることは従来の個別要素法と同じである。
【００１７】
（解析装置のハードウエア構成）
図２は、本実施形態における解析装置が動作する計算装置を示すブロック図である。情報処理装置（コンピュータ）１０は、各種判断及び処理を行う中央処理装置（ＣＰＵ）１１と、処理プログラム及び固定データを格納するＲＯＭ１２と、処理データを記憶するＲＡＭ１３と、入出力回路（Ｉ／Ｏ）１４とから構成されている。また必要に応じて、情報処理装置１０には、Ｉ／Ｏ１４を介して記憶装置、入力装置、表示装置、他のコンピュータなどが接続される。Ｉ／Ｏ１４は、これら外部装置と情報処理装置１０との間でデータをやり取りするための回路である。この情報処理装置１０には、入力データ１５がＩ／Ｏ１４を介して入力される。計算結果は、Ｉ／Ｏ１４を介して出力データ１６として出力される。
【００１８】
（解析装置の機能）
図３は、粒子挙動解析を行う解析プログラムによって実現される機能の構成を示すブロック図である。解析プログラムは、ＲＯＭ１２や記憶装置などのコンピュータ読み取り可能な記録媒体に記憶されており、必要に応じて、ＣＰＵ１１によって読み込まれ実行されるソフトウエアである。この解析プログラムは、主な機能（モジュール）として、入力データ読込部（取得部）Ｂ１１０、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、粒子運動計算部Ｂ１４０、計算結果出力部Ｂ１５０を備えている。全体制御部Ｂ１００は、これらのモジュールの処理の順番やデータの引渡しなどの全体制御を担っている。
【００１９】
入力データ読込部Ｂ１１０は、粒子形状データ、壁形状データ、及び各種パラメータを読み込んでＲＡＭ１３に格納する。粒子形状データとは、粒子の表面を複数の多角形（ポリゴン）で表現したものである。ここで多角形としてはその頂点が一平面上にあることが望ましく、３角形を用いるのが最も単純である。また粒子は凸形状（凸多面体）に限定すると、他の粒子や壁との接触状態の判定が簡便になる。壁形状データも同様であり、壁の表面を複数の多角形で表現したものである。粒子形状データ、壁形状データともに、多角形の頂点座標と、各多角形を構成する頂点の番号とで構成されるデータである。各種パラメータには、粒子と壁の機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）、粒子に働く力（重力、付着力、電磁力、空気の粘性特性等）、及び計算条件としての計算刻み時間（Δｔ）、計算終了時刻がある。なお、壁とは、粒子が移動可能な空間を形成する（規定する）境界面である。
【００２０】
接触判定部Ｂ１２０は、粒子の位置情報（各頂点の座標、粒子の中心座標）に基づいて
、粒子が他の物体（周囲の粒子または壁）と接触しているかどうか、またどのように接触しているかを判定する。粒子受力計算部Ｂ１３０は、各粒子が他の物体から受ける力を計算する機能である。粒子は他の粒子または壁と接触している場合に、そのめり込み量に応じて接触抗力を受ける。粒子運動計算部Ｂ１４０は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。計算結果出力部Ｂ１５０は、得られた各粒子の位置と速度の結果を出力データ１６として出力するものである。
【００２１】
以下に、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、及び粒子運動計算部Ｂ１４０の詳しい内容と計算方法を示す。
【００２２】
（接触判定部Ｂ１２０）
まずは、接触判定部Ｂ１２０における粒子同士の接触判定の内容を図４を用いて説明する。図４の上段は４種類の接触パターンを３次元的に示したものであり、下段は各パターンにおける２つの粒子の接触によるめり込みの様子を示すものである。図４において、実線は６面体粒子ｉ、破線は６面体粒子ｊを示す。ｉ−ｖ０は粒子ｉを構成する多角形の頂点、ｉ−ｅ０は粒子ｉを構成する多角形の辺、ｊ−ｅ０とｊ−ｅ１は粒子ｊを構成する多角形の辺、ｊ−ｆ０は粒子ｊを構成する多角形の面、星印は粒子ｉと粒子ｊの接触点を示す。また図４の上段における矢印は粒子ｉが粒子ｊに衝突する方向を示し、また下段における黒丸は粒子ｊの辺を示す。粒子が凸形状である場合、粒子ｉと粒子ｊの接触は、図４に示す次の（ａ）〜（ｄ）の４つのパターンのいずれかとなる。なお、図４では、粒子同士の接触判定を例に挙げているが、粒子と壁の接触判定も同じように考えることができる。
【００２３】
（ａ）頂点のめり込み：
粒子ｉを構成する多角形の頂点が粒子ｊの１つの面にめり込む場合である。粒子ｉの頂点が粒子ｊの内部にあれば、接触していると判定する。図４（ａ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込んでいる。このパターンでは、頂点ｉ−ｖ０を、粒子ｉの粒子ｊに対する接触点として考える。
【００２４】
（ｂ）辺のめり込み（ケース１）：
上記（ａ）に加えて、めり込んだ頂点を端点に持つ粒子ｉの辺が粒子ｊの辺にもめり込む場合である。（ａ）のめり込んだ頂点を端点に持つ粒子ｉの辺において、粒子ｊ構成面との交点数が１であり、かつその交点を持つ面が（ａ）におけるめり込み面でなければ、これに該当する。図４（ｂ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込む（上記（ａ）に該当）とともに、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にもめり込んでいる。このパターンでは、頂点ｉ−ｖ０に加え、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点も、粒子ｉの粒子ｊに対する接触点として考える。
【００２５】
（ｃ）辺のめり込み（ケース２）：
粒子ｉの辺が粒子ｊの１つの辺だけにめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、その交点を持つ２つの面が共有辺を持ち、かつその共有辺が粒子ｉの辺に近接していれば、これに該当する。図４（ｃ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にめり込んでいる。このパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２６】
（ｄ）辺のめり込み（ケース３）：
粒子ｉの辺が粒子ｊの２つの辺にめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、（ｃ）ではない場合がこれに該当する。図４（ｄ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０とｊ−ｅ１にめり込んでいる。こ
のパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点、及び、辺ｊ−ｅ１に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２７】
上記の接触判定を行うには、図４（ａ）の場合には、粒子ｉの頂点と粒子ｊの面の位置関係、図４（ｂ）〜（ｄ）においては、粒子ｉの辺と粒子ｊの面の交差を調べる必要がある。辺と面のめり込み量の計算を行うための接触判定の仕方を６面体の粒子を例に、図５を用いて説明する。
【００２８】
図５の２０は粒子ｊを構成する１つの多角形の面を表し、この法線ベクトルをｎとする。２１は粒子ｉを構成する辺の１つｉ−ｅ０を表している。この辺２１の両端の頂点をそれぞれＰ０（ｐｘ０，ｐｙ０，ｐｚ０），Ｐ１（ｐｘ１，ｐｙ１，ｐｚ１）とする。また、頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を含む平面（無限）と辺２１の交差点を点Ｑとする。このとき、点Ｑが多角形面２０上にあるか否かの判定を行う。
【００２９】
まず、粒子ｉの辺２１が粒子ｊの多角形面２０の４つの頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を含む平面（無限）を横切っているかどうかを調べる。
【００３０】
粒子ｊの多角形面２０上の一つの任意点Ａ０と粒子ｉの辺２１の頂点Ｐ０，Ｐ１をそれぞれ結び、ベクトルＰ０Ａ０，Ｐ１Ａ０とする。ベクトルＰ０Ａ０は、点Ｐ０から点Ａ０に向かうベクトルであり、ベクトルＰ１Ａ０は、点Ｐ１から点Ａ０に向かうベクトルである。
【００３１】
このとき、２つのベクトルＰ０Ａ０，Ｐ１Ａ０と法線ベクトルｎの内積Ｌ０＝Ｐ０Ａ０・ｎとＬ１＝Ｐ１Ａ０・ｎが同じ符号ならば、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っていないことがわかる。逆に、内積Ｌ０とＬ１の符号が異なる場合は、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っていることがわかる。
【００３２】
次に、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っている場合、粒子ｊの多角形面２０の頂点を含む平面（無限）と粒子ｉの辺２１との交点Ｑの座標（ｑｘ，ｑｙ，ｑｚ）を下記式で求める。
【数１】
ここで、ｒは内分比であり、ｒ＝Ｌ０／（Ｌ０−Ｌ１）で求まる。
【００３３】
次に、交点Ｑと頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を結ぶ４つのベクトルＱｖ０，Ｑｖ１，Ｑｖ２，Ｑｖ３について、隣のベクトルとの外積を下記式を用いて求め、ベクトルＶ０，Ｖ１，Ｖ２，Ｖ３を得る。
【数２】
【００３４】
ベクトルＶ０は、ベクトルＱｖ０とＱｖ１の両方に垂直なベクトル、すなわち多角形面２０に対し垂直なベクトルとなり、その向きはベクトルＱｖ０とＱｖ１のなす角が０〜πかπ〜２πかで反転する。他のベクトルＶ１〜Ｖ３も同様の性質をもつ。交点Ｑが多角形面２０の内部にある場合は、Ｑｖ０とＱｖ１、Ｑｖ１とＱｖ２、Ｑｖ２とＱｖ３、Ｑｖ３とＱｖ０のいずれのなす角も０〜πとなるから、ベクトルＶ０〜Ｖ３の向きは一致する。一方、交点Ｑが多角形面２０の外部にある場合は、なす角が０〜πとなる組と、π〜２πとなる組とがでるため、ベクトルＶ０〜Ｖ３の向きは一致しなくなる。このような性質を利用することで、交点Ｑが多角形面２０の内側に存在するか否か、つまり粒子ｉの辺２１と粒子ｊの多角形面２０とが交差するか否かを判定することができる。
【００３５】
ここでは、下記式のように内積Ｍ０〜Ｍ３を計算し、それらの符号を調べることで判定を行う。
【数３】
Ｍ０〜Ｍ３の符号が全て同じ場合は、粒子ｊの多角形面２０と粒子ｉの辺２１は交差しており、一つでも符号が逆のものがあれば交差していないことがわかる。
【００３６】
上記のような接触判定では、１つの四角形面に対して１つの辺が交差していることを判定するのに、４３個の加減算、４７個の乗算、５回の条件演算を要する。そのため、粒子の多面体の面数が多いほど、また粒子の数が多いほど、膨大な計算時間がかかってしまう。そこで、本実施形態では接触判定の前に接触可能頂点を絞り込み、その接触可能頂点を含む辺、面のみを接触の可能性がある辺、面として接触判定を行う。以下、接触の可能性がある辺及び面を抽出する方法として、２つの方法を説明する。
【００３７】
（第１の抽出方法）
図６は、辺及び面の第１の抽出方法について説明する図である。粒子ｉ、ｊは３次元形状であるが、簡単のため２次元で表示している。
【００３８】
第１の抽出方法は、粒子ｉの粒子ｊに対する接触を判定する際に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わる２つの平面６１，６２を設定し、平面６１，６２で挟まれた空間に位置する粒子ｉの辺及び粒子ｊの面のみを判定対象として抽出するものである。ここで、平面６１（第１の平面）と平面６２（第２の平面）はそれぞれ粒子ｊ，ｉに接し、且つ、互いに平行になるように設定される。言い換えると、第１の抽出方法は、平面６１よりも粒子ｊの中心側の空間に位置する粒子ｉの辺と、平面６２よりも粒子ｉの中心側の空間に位置する粒子ｊの面とを、判定対象として抽出する方法である。なお、辺又は面が「空間に位置する」とは、辺又は面の少なくとも一部分が空間内に存在すること（空間と交わりを有すること）を意味する。
【００３９】
第１の抽出方法を実現する具体的なアルゴリズムとして、本実施形態の接触判定部Ｂ１２０では次のような計算を行う。粒子ｉの中心から粒子ｊの中心に向かう中心結線ベクトルを考え、原点を通り中心結線ベクトルに垂直な面を基準面Ｓとする。粒子ｉの各頂点と基準面Ｓとの距離Ｄ１を計算し、最大値をＤ１ｍａｘとする。粒子ｊの各頂点と基準面Ｓとの距離を計算し、最小値をＤ１ｍｉｎとする。基準面Ｓとの距離がＤ１ｍｉｎ以上Ｄ１
ｍａｘ以下となる粒子ｉ、ｊの各頂点を接触可能頂点として選ぶ。図６で塗りつぶし（黒丸）で示されている頂点が接触可能頂点である。そして、接触可能頂点を含む辺及び面が判定対象として抽出される。図６で実線で示されている辺が抽出された辺である。
【００４０】
Ｄ１ｍａｘ、Ｄ１ｍｉｎの算出方法をより詳しく説明する。中心結線ベクトルと向きが同じで原点を起点とする単位ベクトルをｎ１＝（ｎ１ｘ，ｎ１ｙ，ｎ１ｚ）とし、粒子ｉの各頂点の位置ベクトルをＰｉｋ＝（ｘｉｋ，ｙｉｋ，ｚｉｋ）とする。ｎ１のことを中心結線単位ベクトルとも称す。ここで、ｋは粒子ｉの各構成頂点の番号である。基準面Ｓと粒子ｉの各頂点との距離Ｄ１は、下記式のように、単位ベクトルｎ１と頂点の位置ベクトルＰｉｋの内積ｎ１・Ｐｉｋで算出される。
【数４】
同様に、粒子ｊの各頂点の位置ベクトルをＰｊｋ＝（ｘｊｋ，ｙｊｋ，ｚｊｋ）とすると、基準面Ｓと粒子ｊの各頂点との距離Ｄ１は下記式により算出される。
【数５】
【００４１】
粒子ｉ及び粒子ｊの全ての頂点について距離Ｄ１を求めた後、Ｄ１ｍａｘ及びＤ１ｍｉｎは下記式で算出される。
【数６】
【００４２】
なお、説明を簡単にするためＤ１を距離と表現したが、実際にはＤ１は内積ｎ１・Ｐｉｋ、内積ｎ１・Ｐｊｋのようにベクトルの内積で定義されるものである。よって、ｎ１とＰｉｋ又はＰｊｋの位置関係によっては、各頂点に対するＤ１の一部またはすべてが負になることもあるが、その場合でも上記式によりＤ１ｍａｘ、Ｄ１ｍｉｎを求めることができる。また、ここでは簡単のためｎ１の起点を原点としたが、ベクトルの平行移動に対して内積の大小は変わらないため、起点をどの位置にとってもよい。また、ここでは単位ベクトルｎ１を中心結線ベクトルに平行にとったが、ベクトルｎ１の向きはこれに限られない。ベクトルｎ１が中心結線ベクトルに対して直交していなければ（つまり、平面６１，６２が粒子ｉ，ｊの中心を結んだ直線に交わっていれば）、同様の処理が可能である。
【００４３】
以上述べた方法により接触可能頂点を選択し、判定対象とする粒子ｉの辺と粒子ｊの面を絞り込むことで、接触判定を行う辺、面の組み合わせを減らし、計算時間を短縮することができる。
【００４４】
（第２の抽出方法）
上述した第１の抽出方法では、接触判定を行う範囲を平面６１と６２で挟まれた空間に限定している。しかし、図６で示されるように、平面６１と６２のあいだに、実際には他の粒子と接触していない頂点が多く含まれる場合もある。そこで、第２の抽出方法として、接触可能頂点をさらに絞り込む方法について説明する。
【００４５】
第２の抽出方法は、図７に示すように、平面６１，６２とは異なる角度で平面７１，７２を設定し、４つの平面６１，６２，７１，７２で挟まれた空間に位置する粒子ｉの辺及
び粒子ｊの面のみを判定対象として抽出するものである。ここで、平面７１（第３の平面）と平面７２（第４の平面）はそれぞれ粒子ｊ，ｉに接し、且つ、互いに平行になるように設定される。言い換えると、第２の抽出方法は、平面６１及び平面７１よりも粒子ｊの中心側の空間に位置する粒子ｉの辺と、平面６２及び平面７２よりも粒子ｉの中心側の空間に位置する粒子ｊの面とを、判定対象として抽出する方法である。
【００４６】
第２の抽出方法を実現する具体的なアルゴリズムを説明する。図７は２つのベクトルによる接触可能頂点の絞り込み方法を示す図である。図７左上は上記の中心結線単位ベクトルｎ１によって絞り込まれた接触可能頂点を示している。塗りつぶし（黒丸）で示されているのが接触可能頂点である。
【００４７】
図７左下は、重なり最小単位ベクトルｎ２によって絞り込まれた接触可能頂点を表している。「重なり最小単位ベクトル」とは、粒子ｉと粒子ｊの重なりを最小とする向きの単位ベクトルをいい、以降の説明では単に「重なり最小ベクトル」とも称す。ｎ２による接触可能頂点の絞り込み方法はｎ１による接触可能頂点の絞り込み方法と同じである。つまり、ベクトルｎ２と粒子ｉの各頂点座標の位置ベクトルＰｉｋの内積ｎ２・Ｐｉｋの最大値Ｄ２ｍａｘと、ベクトルｎ２と粒子ｊの各頂点座標の位置ベクトルＰｊｋの内積ｎ２・Ｐｊｋの最小値Ｄ２ｍｉｎを求める。そして、内積がＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる頂点を接触可能頂点として抽出する。
【００４８】
図７右は中心結線ベクトルによる接触可能頂点集合をＡ、重なり最小ベクトルによる接触可能頂点集合をＢとしたときに、共通部分である積集合Ａ∩Ｂで選択される頂点を表している。図７において、集合Ａに含まれる頂点数は１０、集合Ａ内の頂点を構成頂点とする辺数は１２、集合Ｂに含まれる頂点数は８、辺数は１０である。一方、積集合Ａ∩Ｂに含まれる頂点数は４、辺数は６となり、１つのベクトルのみの場合より判定対象を絞り込むことができる。よって、積集合Ａ∩Ｂに含まれる頂点を最終的な接触可能頂点と定義する。
【００４９】
積集合Ａ∩Ｂが空集合でも接触する場合がある。図８を用いて説明する。四角で表されている頂点が中心結線ベクトルにより選択された頂点（集合Ａ）、三角で表されているのが重なり最小ベクトルにより選択された頂点（集合Ｂ）である。塗りつぶしの円で表されているのが積集合Ａ∩Ｂに含まれる接触可能頂点である。粒子ｊは接触可能頂点を持たないが、図８に示すように粒子ｉと接触している。このような場合を含めるため、判定対象とする辺、面は次のように決定する。
・積集合Ａ∩Ｂに含まれる頂点（接触可能頂点）を持つ辺、面。
・集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を持つ辺、面。
【００５０】
以上の２パターンによって、全ての接触の可能性がある辺と面を選択することができる。本方法により、接触判定を行う辺、面の組み合わせを大幅に減らすことができ、計算時間のより一層の短縮が可能となる。
【００５１】
本方法で用いる重なり最小ベクトルｎ２の算出方法について説明する。図９は重なり最小ベクトルを求めるフローチャートである。図９に示す処理は、接触判定部Ｂ１２０が実行するものである。
【００５２】
まず、重なり最小ベクトルの初期値を設定する（ステップＳ１）。初期値は中心結線単位ベクトルと平行な単位ベクトルとする。粒子ｉの中心位置ベクトルをＣｉ、粒子ｊの中心位置ベクトルをＣｊとすると、重なり最小ベクトルｎ２の初期値は下記式で表される。
【数７】
以降の処理で、ベクトルｎ２の向きを微小に変動させながら、粒子ｉと粒子ｊの重なりが最小となるものを見つける。その処理で用いる微小量δの初期値δｍａｘ（＝０．１程度）をステップＳ１で設定する。
【００５３】
次に、２粒子の重なり量Ｌ（＝Ｄ２ｍａｘ−Ｄ２ｍｉｎ）を算出する（ステップＳ２）。
次に、接触の可能性があるか否かを判定する（ステップＳ３）。Ｌが０以下の場合、２粒子は非接触であるので、接触判定から除外される（ステップＳ４）。Ｌが０より大きい場合は接触の可能性があるため、重なり最小ベクトルｎ２を更新（探索）するためのループに入る。
【００５４】
重なり最小ベクトルｎ２に直交し、たがいに直交するベクトルｕ、ｖを適当に定める（ステップＳ５）。そして、試行する重なり最小ベクトルｎ２（１）〜ｎ２（４）を下記式のように設定する（ステップＳ６）。
【数８】
上記式の１行目で決定される試行重なり最小ベクトルｎ２（１）は、図９中に示されるように、直交ベクトルｕの方向へ微小量δで決まる量だけ回転させたものとなる。
【００５５】
このようにして決定された各試行重なり最小ベクトルｎ２（ｋ）に対して２粒子の重なり量Ｌｎｅｗ（ｋ）を新たに算出する（ｋ＝１，２，３，４）。算出された４つのＬｎｅｗ（ｋ）の最小値とＬを比較する（ステップＳ７）。ｍｉｎ｛Ｌｎｅｗ（ｋ）｝がＬより小さければ、重なり最小ベクトルｎ２をｎ２（ｋ）に更新する。また、ＬもＬｎｅｗ（ｋ）に更新する（ステップＳ８）。重なり量Ｌを最小とするループにおいて、ｍｉｎ｛Ｌｎｅｗ（ｋ）｝がＬより大きくなった場合は（ステップＳ９）、微小量δをδ／２と小さくした後、試行重なりベクトルの設定から繰り返す（ステップＳ１０）。
【００５６】
上記のループを繰り返して、δがδｍｉｎ（＝０．００１程度）より小さくなったときのｎ２が最終的な重なり最小ベクトルとなる（ステップＳ１１）。このような処理により、粒子ｉとｊの重なり量Ｌを最も小さくするベクトルｎ２を得ることができる。
【００５７】
なお、図９のステップＳ５で示したように、ループを繰り返すたびにｕ，ｖを面内で回転させると、局所解への陥りを抑制し、重なり最小ベクトルの最適解を得やすいという効果が得られる。
【００５８】
図９で示した方法は重なり最小ベクトルｎ２を決定するための一手法にすぎない。他の公知の最適化手法を利用して、重なり量Ｌを最小とするベクトルｎ２を求めることも可能である。また、少なくともベクトルｎ１とｎ２の向きが異なっていれば、第１の抽出方法よりも判定対象を絞り込む効果は得られるため、ベクトルｎ１に対して任意の回転を行ったものをベクトルｎ２として用いるだけでもよい。
【００５９】
（接触判定処理）
接触判定部Ｂ１２０の処理フローを図１０で説明する。
まず、接触判定を行う２粒子ｉ、ｊを設定する（ステップＳ１００）。次に２粒子ｉ、ｊの中心座標から中心結線単位ベクトルｎ１を求める（ステップＳ１１０）。ｎ１を基に２粒子ｉ、ｊの重なり量Ｌを算出する（ステップＳ１２０）。重なり量Ｌを基に接触可能性を判定し（ステップＳ１３０）、Ｌ＞０でなければ、２粒子は非接触として接触判定を終了する（ステップＳ１３５）。次に、２粒子ｉ、ｊの位置座標情報から重なり最小単位ベクトルｎ２を求める（ステップＳ１４０）。ステップＳ１４０の処理内容は、図９で詳しく説明したとおりである。
【００６０】
ここで、ステップＳ１１０で設定したｎ１を用いて接触可能頂点集合Ａを設定し（ステップＳ１５０）、ステップＳ１４０で設定したｎ２を用いて接触可能頂点集合Ｂを設定する（ステップＳ１６０）。次に、接触可能頂点集合Ａ∩Ｂを設定する（ステップＳ１７０）。こうして設定した接触可能頂点を含む粒子ｉの辺と粒子ｊの面を抽出する（ステップＳ１８０）。以上で、接触判定に用いる判定対象が絞り込まれる。
【００６１】
次に、粒子ｉが粒子ｊに接触しているか否かの判定処理を行う。まず、接触判定を行う粒子ｉの辺を設定する（ステップＳ１９０）。接触判定を行う粒子ｊの面を設定する（ステップＳ２００）。ここで設定した粒子ｉの辺と粒子ｊの面との交差の有無を判定する（ステップＳ２１０）。ここで、粒子ｊの面に粒子ｉの辺が交差していれば、交差回数のカウントアップを行う（ステップＳ２２０）。この操作をステップＳ１８０で抽出された粒子ｊの面の全てに対して行う（ステップＳ２３０）。次に、粒子ｉの辺と粒子ｊの面の交差回数が２回かどうかの判定を行う（ステップＳ２４０）。交差回数が２回ならば、めり込み量を算出し、メモリに記憶しておく（ステップＳ２５０）。本操作をＳ１８０で抽出された粒子ｉの辺の全てに対して行う（ステップＳ２６０）。
【００６２】
なお、図１０では、説明を簡単にするため、粒子ｉのある辺が粒子ｊの２つの面と交差した場合に、接触有りと判定し、粒子ｉとｊのあいだのめり込み量を算出する処理例を示した。しかし実際の処理では、図４に示したように、粒子ｉと粒子ｊの接触を４つのパターンで判定する。前述したようにパターンごとに判定条件が異なるが、いずれのパターンにおいても粒子ｉの辺と粒子ｊの面との交差判定が基本となるので、図１０のステップＳ１８０で抽出した辺及び面を利用すればよい。
【００６３】
各パターン（ａ）〜（ｄ）におけるめり込み量の定義について、図４中の下段の図を用いて説明する。
（ａ）頂点のめり込み：粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０との距離をめり込み量とする。
（ｂ）辺のめり込み（ケース１）：粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｃ）辺のめり込み（ケース２）：粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｄ）辺のめり込み（ケース３）：粒子ｉの辺ｉ−ｅ０に対する、粒子ｊの辺ｊ−ｅ０、ｊ−ｅ１との距離をそれぞれの辺とのめり込み量とする。
【００６４】
（粒子受力計算部Ｂ１３０）
以上の４パターンにおいて、めり込み量（距離）の計算に用いた各粒子表面の２点を結ぶ方向を接触面法線方向と定める。この接触面法線方向は、接触点の（粒子ｊの面又は辺に対する）めり込み方向ということもできる。そして、接触判定部Ｂ１２０で求めためり込み量に基づいて、ヘルツの式を用いて、粒子ｉに働く接触面法線方向の接触力（バネ力）を求める。また各接触点に対して、接触が始まってからの接触面（上記接触面法線方向に垂直な面）内方向の移動量を求め、ミンドリンの式を用いて粒子ｉに働く接触面内方向の接触力（バネ力）を求める。また接触面法線方向の運動速度（めり込み速度）と面内方向の各運動速度（ずり速度）から、各方向に作用する粘性力（ダンパー力）を求める。そしてそれら接触力と粘性力の和に、粒子ｉに働くその他の力（重力、付着力、電磁力、空気の粘性抵抗等）を加えて、各粒子ｉが受ける並進力と回転力を求める。なお、本めり込み量から接触力を求める過程は、上記非特許文献１に記載されている手法と同じであるため、詳しい説明は割愛する。
【００６５】
（粒子運動計算部Ｂ１４０）
次に、粒子運動計算部Ｂ１４０における各粒子の刻み時間後の速度と位置の算出方法を説明する。この処理は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。すなわち、粒子運動計算部Ｂ１４０は、並進力からニュートンの運動方程式を用いて並進方向の加速度を求め、更に速度と位置を求める。また粒子運動計算部Ｂ１４０は、回転力からオイラーの運動方程式を用いて角加速度を求め、更に回転速度と回転位置を求める。本処理は、非特許文献１に記載されている手法と同じであるため、詳しい説明は割愛する。
【００６６】
（実施例）
以上の接触判定方法を用いて粒子挙動計算を実施した例を、図１１〜図１２を用いて説明する。図１１（ａ）は解析モデルを示す。解析モデルは電子写真技術を用いた画像形成装置におけるクリーニングプロセスである。０．１ｍ／ｓｅｃで矢印の方向に移動する感光体ドラムに対してクリーニングブレードをβ＝３０°で当接させておき、感光体ドラム上に乗って運ばれてきたトナーの運動を見る。各トナーの形状を粒子形状データとして与え、感光体ドラム及びクリーニングブレードの形状を壁形状データとして与える。
【００６７】
図１１（ｂ）は、解析条件として与えた、トナー、感光体ドラム、クリーニングブレードそれぞれの機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）を示している。この条件で、刻み時間Δｔを１０ｎｓｅｃに設定し、約２百万ステップの計算を実行した。
【００６８】
図１２は、（ａ）球形粒子と（ｂ）異形粒子の計算結果を比較したものである。球形粒子の計算は従来の計算方法で実施している。異形粒子は４面体の角をとった形状であり、その表面を８つの３角形と６つの４角形の合計１２の多角形で構成している。なお本計算においては、トナーがブレードに堰き止められた際の挙動をみるため、強制的にトナーがすり抜けないようにしている。球形粒子に対して異形粒子は回転しにくいことから、全体的に動き辛く、トナー粒子の密集度が上がる傾向にある。図１２の計算結果では、異形粒子の密集度が高くなっている様子が再現されている。
【００６９】
このトナーがブレードに堰き止められた際の挙動は、クリーニング不良であるトナーのすり抜けとの関係があり、流動的である程すり抜けやすいといわれている。このことから、製品設計を行う上で、トナーの形状によるクリーニング不良の傾向を粒子挙動計算で予測することが可能となる。
【００７０】
また、図１３は面数の異なる異形粒子の自由落下の計算を示した図である。（ａ）は、
頂点数１２、面数１４、辺数２４の場合の異形粒子で、（ｂ）は頂点数９８、面数１０８、辺数２０４の異形粒子の場合である。（ａ）と（ｂ）の場合において、２５６個の異形粒子の自由落下の計算を刻み時間Δｔ＝１０ｎｓｅｃで８百万ステップ実行し、本発明を使用しない場合と使用した場合での計算時間の比較を行った。（ａ）では、従来の方法では計算に１９時間３５分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込み（第１の抽出方法）を用いた場合では、８時間０９分となり、高速化率は２．４倍となる。（ｂ）では、従来の方法では計算に５２４時間３２分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込みを用いると２２時間２４分となり、高速化率は２３．４倍となった。このことから、粒子の面の数が多くなればなるほど、本発明を用いることによる計算時間の短縮効果が大きくなることがわかる。
【００７１】
図１４は中心結線ベクトルと重なり最小ベクトルを用いた高速化の効果を示す図である。計算には頂点数１９８、面数３９２、辺数５８８の異形粒子を用い、図１３と同様の自由落下計算を３０００００ステップ、Δｔ＝１０ｎｓｅｃで行った。接触判定に中心結線ベクトルのみを用いた場合（第１の抽出方法）と中心結線ベクトルと重なり最小ベクトルの両方を用いた場合（第２の抽出方法）で計算時間を比較した。中心結線ベクトルのみを用いた場合は４７２６秒、中心結線ベクトルと重なり最小ベクトルの両方を用いた場合は２００９秒となり、高速化率は２．３５倍となった。このことから、中心結線ベクトルと重なり最小ベクトルの両方を用いた接触判定により、さらなる計算時間の短縮が実現できていることが分かる。
【符号の説明】
【００７２】
１０：情報処理装置
Ｂ１００：全体制御部
Ｂ１１０：入力データ読込部
Ｂ１２０：接触判定部
Ｂ１３０：粒子受力計算部
Ｂ１４０：粒子運動計算部
Ｂ１５０：計算結果出力部
【技術分野】
【０００１】
本発明は、密集状態にある多数の粒子について、その形状を考慮して挙動を解析する技術に関するものであり、特に挙動解析計算における粒子同士の接触判定に関する。
【背景技術】
【０００２】
粉体や粒体などの粒子を取り扱う分野では、個々の粒子の挙動を知ることが重要である。従来これらの分野では、実験と観察から粒子の挙動を把握することが多かった。しかし近年、計算機内でその運動を再現させ、その挙動を把握する、いわゆるシミュレーション技術の適用が活発に行われている。正確な粒子挙動のシミュレーションを行うには、粒子のもつ形状、機械的特性等の情報をできるだけ詳細に計算モデルに反映させ、それらの因子を考慮して解を求めることが重要である。また、土木分野における砕石、砂利等の比較的大きな物体も「粒子」と考えることにより、上記シミュレーション技術をこれらの物体の挙動を計算する場合へも適用可能である。
【０００３】
一方、このような粒子を扱う装置の１つに、プリンタ、複写機、ファクシミリ等の電子写真技術を用いた画像形成装置がある。これらの装置では、直径数ミクロンの粉体粒子からなるトナーを用いて、帯電、露光、現像、転写、定着、クリーニングという６つのプロセスを経て、紙上に画像を形成する。このような画像形成装置では、トナー粒子の形状が、作成される画像に大きな影響を与えることがわかっている。そこでこのような画像形成装置の開発においても、シミュレーション技術を用いて、密集状態にある個々のトナー粒子の挙動を求める試みが活発に行われ、近年の傾向として特にトナー粒子の形状を考慮することが求められている。
【０００４】
粒子が密集した状態での各粒子の挙動をシミュレーションする代表的な方法として個別要素法（ＤＥＭ）がある。個別要素法の技術は非特許文献１に詳細に記載されているが、個々の粒子を球形で表現し、それらの接触に際してはそのめり込み量からヘルツの式とミンドリンの式を用いることで粒子の持つ機械特性を考慮した接触抗力を求める。そしてニュートンの運動方程式（並進方向）とオイラーの運動方程式（回転方向）を用いて、その接触抗力をもとに個々の粒子の時間進展に伴う運動を求めるものである。
【０００５】
個別要素法の特徴は個々の粒子を要素と捉え、それらが接触する部分にバネ、ダンパー、スライダを仮定して、その接触力を求め、運動方程式を解くところにある。この考え方を球形以外の形状をもつ粒子の挙動計算に適用したものとして、非特許文献２、非特許文献３がある。これらの文献では粒子の形状を多面体で表現し上述の個別要素法を用いてその挙動を計算している。そして、粒子の接触の検出に関して、Common-planeと呼ばれる面を想定して、接触の有無、接触力の作用点、接触力の作用方向、及び接触力の大きさを決定する。しかしながら、同文献の方法ではCommon-planeと粒子の接触断面積を求めることが必要となるが、断面積を求めるには粒子の各頂点、辺、面とCommon-planeとの接触状態を調べなければならないため、長い計算時間を要する。また本発明者らが同文献の方法を実際に試行しようとしたところ、同文献は接触力を作用させる位置についての情報開示が乏しく、再現することが困難であることが分かった。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】粉体工学会編、粉体シミュレーション入門(1998)
【非特許文献２】P.A. Cundall: Formulation of a three-dimensional distinct element model - Part I. a scheme to detect and represent contacts in a system composed of many polyhedral blocks, Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 25(3), pp.107-116 (1988)
【非特許文献３】R. Hart, P.A. Cundall, and J. Lemos: Formulation of a three-dimensional distinct element model - Part II. mechanical calculations for motion and interaction of a system composed of many polyhedral blocks, Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 25(3), pp.117-125 (1988)
【発明の概要】
【発明が解決しようとする課題】
【０００７】
本発明の目的は、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することにある。また本発明の目的は、密集状態にある多数の粒子の挙動を計算するに際し、粒子同士の接触判定を高速に行うことにある。
【課題を解決するための手段】
【０００８】
上記課題を解決するために、本発明は以下の構成を採用する。
本発明の第一態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。
【０００９】
本発明の第二態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。
【００１０】
本発明の第三態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計
算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。
【００１１】
本発明の第四態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。
【００１２】
本発明の第五態様は、上述した粒子挙動解析方法の各ステップを情報処理装置（コンピュータ）に実行させる解析プログラムである。
【発明の効果】
【００１３】
本発明によれば、粒子同士の接触判定に要する時間を短縮でき、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することができる。
【図面の簡単な説明】
【００１４】
【図１】粒子挙動計算の力学モデルを説明する図。
【図２】解析装置のハードウエア構成を示す図。
【図３】解析装置の機能構成を示す図。
【図４】接触パターンと各パターンにおける粒子めり込み量を示す図。
【図５】辺と面の交差の有無の判定方法を示す図。
【図６】接触可能頂点の絞り込み方法を示す図。
【図７】２ベクトルによる接触可能頂点の絞り込み方法を示す図。
【図８】接触可能頂点を持たない粒子との接触を示す図。
【図９】重なり最小ベクトルを求める方法を示す図。
【図１０】接触判定部の処理フローを示す図。
【図１１】実施例の計算条件を示す図。
【図１２】実施例の計算結果を示す図。
【図１３】中心結線ベクトルを用いた高速化の効果を示す図。
【図１４】中心結線ベクトルと重なり最小ベクトルを用いた高速化の効果を示す図。
【発明を実施するための形態】
【００１５】
本発明の実施形態は、密集状態にある多数の粒子の挙動をコンピュータ・シミュレーションにより計算するための粒子挙動解析方法及び装置を提供するものである。本実施形態に係る粒子挙動解析方法の特徴の一つは、多数の多角形から構成される多面体（ポリゴンモデル）により各粒子の形状を定義した点である。これにより任意の粒子形状を表現でき
るようになり、多様なケースのシミュレーションを精度良く行うことが可能になる。本方法のもう一つの特徴は、粒子同士の接触判定を行う際に、全ての頂点、辺、面について網羅的に接触判定計算を行うのではなく、予め接触可能性のある辺及び面を特定し、それらのあいだについてのみ接触判定を行う点である。これにより、接触判定に要する時間を短縮でき、粒子の挙動計算の高速化を図ることができる。以下、本実施形態の具体的な構成及び特徴について図に基づき詳しく説明する。
【００１６】
（粒子挙動計算のモデル）
図１は、本実施形態における粒子挙動計算の力学モデルを説明するものであり、異形粒子（非球形粒子）が他の物体と接触した際に働く接触力を説明する図である。このモデルでは、粒子が他の物体と接触した際に発生するめり込みを、その接触面に対して法線方向と接線方向に分解して考える。粒子には、それぞれのめり込み量に対してスプリングで押し返す力が働く。また各スプリングと並列にダンパーを設定する。これはめり込み速度に比例して働く力であり、運動を減衰させるものである。またすべりを考慮するためのスライダを設定している。このように、スプリング、ダンパー、スライダを設定して衝突を考えることは従来の個別要素法と同じである。
【００１７】
（解析装置のハードウエア構成）
図２は、本実施形態における解析装置が動作する計算装置を示すブロック図である。情報処理装置（コンピュータ）１０は、各種判断及び処理を行う中央処理装置（ＣＰＵ）１１と、処理プログラム及び固定データを格納するＲＯＭ１２と、処理データを記憶するＲＡＭ１３と、入出力回路（Ｉ／Ｏ）１４とから構成されている。また必要に応じて、情報処理装置１０には、Ｉ／Ｏ１４を介して記憶装置、入力装置、表示装置、他のコンピュータなどが接続される。Ｉ／Ｏ１４は、これら外部装置と情報処理装置１０との間でデータをやり取りするための回路である。この情報処理装置１０には、入力データ１５がＩ／Ｏ１４を介して入力される。計算結果は、Ｉ／Ｏ１４を介して出力データ１６として出力される。
【００１８】
（解析装置の機能）
図３は、粒子挙動解析を行う解析プログラムによって実現される機能の構成を示すブロック図である。解析プログラムは、ＲＯＭ１２や記憶装置などのコンピュータ読み取り可能な記録媒体に記憶されており、必要に応じて、ＣＰＵ１１によって読み込まれ実行されるソフトウエアである。この解析プログラムは、主な機能（モジュール）として、入力データ読込部（取得部）Ｂ１１０、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、粒子運動計算部Ｂ１４０、計算結果出力部Ｂ１５０を備えている。全体制御部Ｂ１００は、これらのモジュールの処理の順番やデータの引渡しなどの全体制御を担っている。
【００１９】
入力データ読込部Ｂ１１０は、粒子形状データ、壁形状データ、及び各種パラメータを読み込んでＲＡＭ１３に格納する。粒子形状データとは、粒子の表面を複数の多角形（ポリゴン）で表現したものである。ここで多角形としてはその頂点が一平面上にあることが望ましく、３角形を用いるのが最も単純である。また粒子は凸形状（凸多面体）に限定すると、他の粒子や壁との接触状態の判定が簡便になる。壁形状データも同様であり、壁の表面を複数の多角形で表現したものである。粒子形状データ、壁形状データともに、多角形の頂点座標と、各多角形を構成する頂点の番号とで構成されるデータである。各種パラメータには、粒子と壁の機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）、粒子に働く力（重力、付着力、電磁力、空気の粘性特性等）、及び計算条件としての計算刻み時間（Δｔ）、計算終了時刻がある。なお、壁とは、粒子が移動可能な空間を形成する（規定する）境界面である。
【００２０】
接触判定部Ｂ１２０は、粒子の位置情報（各頂点の座標、粒子の中心座標）に基づいて
、粒子が他の物体（周囲の粒子または壁）と接触しているかどうか、またどのように接触しているかを判定する。粒子受力計算部Ｂ１３０は、各粒子が他の物体から受ける力を計算する機能である。粒子は他の粒子または壁と接触している場合に、そのめり込み量に応じて接触抗力を受ける。粒子運動計算部Ｂ１４０は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。計算結果出力部Ｂ１５０は、得られた各粒子の位置と速度の結果を出力データ１６として出力するものである。
【００２１】
以下に、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、及び粒子運動計算部Ｂ１４０の詳しい内容と計算方法を示す。
【００２２】
（接触判定部Ｂ１２０）
まずは、接触判定部Ｂ１２０における粒子同士の接触判定の内容を図４を用いて説明する。図４の上段は４種類の接触パターンを３次元的に示したものであり、下段は各パターンにおける２つの粒子の接触によるめり込みの様子を示すものである。図４において、実線は６面体粒子ｉ、破線は６面体粒子ｊを示す。ｉ−ｖ０は粒子ｉを構成する多角形の頂点、ｉ−ｅ０は粒子ｉを構成する多角形の辺、ｊ−ｅ０とｊ−ｅ１は粒子ｊを構成する多角形の辺、ｊ−ｆ０は粒子ｊを構成する多角形の面、星印は粒子ｉと粒子ｊの接触点を示す。また図４の上段における矢印は粒子ｉが粒子ｊに衝突する方向を示し、また下段における黒丸は粒子ｊの辺を示す。粒子が凸形状である場合、粒子ｉと粒子ｊの接触は、図４に示す次の（ａ）〜（ｄ）の４つのパターンのいずれかとなる。なお、図４では、粒子同士の接触判定を例に挙げているが、粒子と壁の接触判定も同じように考えることができる。
【００２３】
（ａ）頂点のめり込み：
粒子ｉを構成する多角形の頂点が粒子ｊの１つの面にめり込む場合である。粒子ｉの頂点が粒子ｊの内部にあれば、接触していると判定する。図４（ａ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込んでいる。このパターンでは、頂点ｉ−ｖ０を、粒子ｉの粒子ｊに対する接触点として考える。
【００２４】
（ｂ）辺のめり込み（ケース１）：
上記（ａ）に加えて、めり込んだ頂点を端点に持つ粒子ｉの辺が粒子ｊの辺にもめり込む場合である。（ａ）のめり込んだ頂点を端点に持つ粒子ｉの辺において、粒子ｊ構成面との交点数が１であり、かつその交点を持つ面が（ａ）におけるめり込み面でなければ、これに該当する。図４（ｂ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込む（上記（ａ）に該当）とともに、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にもめり込んでいる。このパターンでは、頂点ｉ−ｖ０に加え、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点も、粒子ｉの粒子ｊに対する接触点として考える。
【００２５】
（ｃ）辺のめり込み（ケース２）：
粒子ｉの辺が粒子ｊの１つの辺だけにめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、その交点を持つ２つの面が共有辺を持ち、かつその共有辺が粒子ｉの辺に近接していれば、これに該当する。図４（ｃ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にめり込んでいる。このパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２６】
（ｄ）辺のめり込み（ケース３）：
粒子ｉの辺が粒子ｊの２つの辺にめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、（ｃ）ではない場合がこれに該当する。図４（ｄ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０とｊ−ｅ１にめり込んでいる。こ
のパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点、及び、辺ｊ−ｅ１に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２７】
上記の接触判定を行うには、図４（ａ）の場合には、粒子ｉの頂点と粒子ｊの面の位置関係、図４（ｂ）〜（ｄ）においては、粒子ｉの辺と粒子ｊの面の交差を調べる必要がある。辺と面のめり込み量の計算を行うための接触判定の仕方を６面体の粒子を例に、図５を用いて説明する。
【００２８】
図５の２０は粒子ｊを構成する１つの多角形の面を表し、この法線ベクトルをｎとする。２１は粒子ｉを構成する辺の１つｉ−ｅ０を表している。この辺２１の両端の頂点をそれぞれＰ０（ｐｘ０，ｐｙ０，ｐｚ０），Ｐ１（ｐｘ１，ｐｙ１，ｐｚ１）とする。また、頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を含む平面（無限）と辺２１の交差点を点Ｑとする。このとき、点Ｑが多角形面２０上にあるか否かの判定を行う。
【００２９】
まず、粒子ｉの辺２１が粒子ｊの多角形面２０の４つの頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を含む平面（無限）を横切っているかどうかを調べる。
【００３０】
粒子ｊの多角形面２０上の一つの任意点Ａ０と粒子ｉの辺２１の頂点Ｐ０，Ｐ１をそれぞれ結び、ベクトルＰ０Ａ０，Ｐ１Ａ０とする。ベクトルＰ０Ａ０は、点Ｐ０から点Ａ０に向かうベクトルであり、ベクトルＰ１Ａ０は、点Ｐ１から点Ａ０に向かうベクトルである。
【００３１】
このとき、２つのベクトルＰ０Ａ０，Ｐ１Ａ０と法線ベクトルｎの内積Ｌ０＝Ｐ０Ａ０・ｎとＬ１＝Ｐ１Ａ０・ｎが同じ符号ならば、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っていないことがわかる。逆に、内積Ｌ０とＬ１の符号が異なる場合は、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っていることがわかる。
【００３２】
次に、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っている場合、粒子ｊの多角形面２０の頂点を含む平面（無限）と粒子ｉの辺２１との交点Ｑの座標（ｑｘ，ｑｙ，ｑｚ）を下記式で求める。
【数１】
ここで、ｒは内分比であり、ｒ＝Ｌ０／（Ｌ０−Ｌ１）で求まる。
【００３３】
次に、交点Ｑと頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を結ぶ４つのベクトルＱｖ０，Ｑｖ１，Ｑｖ２，Ｑｖ３について、隣のベクトルとの外積を下記式を用いて求め、ベクトルＶ０，Ｖ１，Ｖ２，Ｖ３を得る。
【数２】
【００３４】
ベクトルＶ０は、ベクトルＱｖ０とＱｖ１の両方に垂直なベクトル、すなわち多角形面２０に対し垂直なベクトルとなり、その向きはベクトルＱｖ０とＱｖ１のなす角が０〜πかπ〜２πかで反転する。他のベクトルＶ１〜Ｖ３も同様の性質をもつ。交点Ｑが多角形面２０の内部にある場合は、Ｑｖ０とＱｖ１、Ｑｖ１とＱｖ２、Ｑｖ２とＱｖ３、Ｑｖ３とＱｖ０のいずれのなす角も０〜πとなるから、ベクトルＶ０〜Ｖ３の向きは一致する。一方、交点Ｑが多角形面２０の外部にある場合は、なす角が０〜πとなる組と、π〜２πとなる組とがでるため、ベクトルＶ０〜Ｖ３の向きは一致しなくなる。このような性質を利用することで、交点Ｑが多角形面２０の内側に存在するか否か、つまり粒子ｉの辺２１と粒子ｊの多角形面２０とが交差するか否かを判定することができる。
【００３５】
ここでは、下記式のように内積Ｍ０〜Ｍ３を計算し、それらの符号を調べることで判定を行う。
【数３】
Ｍ０〜Ｍ３の符号が全て同じ場合は、粒子ｊの多角形面２０と粒子ｉの辺２１は交差しており、一つでも符号が逆のものがあれば交差していないことがわかる。
【００３６】
上記のような接触判定では、１つの四角形面に対して１つの辺が交差していることを判定するのに、４３個の加減算、４７個の乗算、５回の条件演算を要する。そのため、粒子の多面体の面数が多いほど、また粒子の数が多いほど、膨大な計算時間がかかってしまう。そこで、本実施形態では接触判定の前に接触可能頂点を絞り込み、その接触可能頂点を含む辺、面のみを接触の可能性がある辺、面として接触判定を行う。以下、接触の可能性がある辺及び面を抽出する方法として、２つの方法を説明する。
【００３７】
（第１の抽出方法）
図６は、辺及び面の第１の抽出方法について説明する図である。粒子ｉ、ｊは３次元形状であるが、簡単のため２次元で表示している。
【００３８】
第１の抽出方法は、粒子ｉの粒子ｊに対する接触を判定する際に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わる２つの平面６１，６２を設定し、平面６１，６２で挟まれた空間に位置する粒子ｉの辺及び粒子ｊの面のみを判定対象として抽出するものである。ここで、平面６１（第１の平面）と平面６２（第２の平面）はそれぞれ粒子ｊ，ｉに接し、且つ、互いに平行になるように設定される。言い換えると、第１の抽出方法は、平面６１よりも粒子ｊの中心側の空間に位置する粒子ｉの辺と、平面６２よりも粒子ｉの中心側の空間に位置する粒子ｊの面とを、判定対象として抽出する方法である。なお、辺又は面が「空間に位置する」とは、辺又は面の少なくとも一部分が空間内に存在すること（空間と交わりを有すること）を意味する。
【００３９】
第１の抽出方法を実現する具体的なアルゴリズムとして、本実施形態の接触判定部Ｂ１２０では次のような計算を行う。粒子ｉの中心から粒子ｊの中心に向かう中心結線ベクトルを考え、原点を通り中心結線ベクトルに垂直な面を基準面Ｓとする。粒子ｉの各頂点と基準面Ｓとの距離Ｄ１を計算し、最大値をＤ１ｍａｘとする。粒子ｊの各頂点と基準面Ｓとの距離を計算し、最小値をＤ１ｍｉｎとする。基準面Ｓとの距離がＤ１ｍｉｎ以上Ｄ１
ｍａｘ以下となる粒子ｉ、ｊの各頂点を接触可能頂点として選ぶ。図６で塗りつぶし（黒丸）で示されている頂点が接触可能頂点である。そして、接触可能頂点を含む辺及び面が判定対象として抽出される。図６で実線で示されている辺が抽出された辺である。
【００４０】
Ｄ１ｍａｘ、Ｄ１ｍｉｎの算出方法をより詳しく説明する。中心結線ベクトルと向きが同じで原点を起点とする単位ベクトルをｎ１＝（ｎ１ｘ，ｎ１ｙ，ｎ１ｚ）とし、粒子ｉの各頂点の位置ベクトルをＰｉｋ＝（ｘｉｋ，ｙｉｋ，ｚｉｋ）とする。ｎ１のことを中心結線単位ベクトルとも称す。ここで、ｋは粒子ｉの各構成頂点の番号である。基準面Ｓと粒子ｉの各頂点との距離Ｄ１は、下記式のように、単位ベクトルｎ１と頂点の位置ベクトルＰｉｋの内積ｎ１・Ｐｉｋで算出される。
【数４】
同様に、粒子ｊの各頂点の位置ベクトルをＰｊｋ＝（ｘｊｋ，ｙｊｋ，ｚｊｋ）とすると、基準面Ｓと粒子ｊの各頂点との距離Ｄ１は下記式により算出される。
【数５】
【００４１】
粒子ｉ及び粒子ｊの全ての頂点について距離Ｄ１を求めた後、Ｄ１ｍａｘ及びＤ１ｍｉｎは下記式で算出される。
【数６】
【００４２】
なお、説明を簡単にするためＤ１を距離と表現したが、実際にはＤ１は内積ｎ１・Ｐｉｋ、内積ｎ１・Ｐｊｋのようにベクトルの内積で定義されるものである。よって、ｎ１とＰｉｋ又はＰｊｋの位置関係によっては、各頂点に対するＤ１の一部またはすべてが負になることもあるが、その場合でも上記式によりＤ１ｍａｘ、Ｄ１ｍｉｎを求めることができる。また、ここでは簡単のためｎ１の起点を原点としたが、ベクトルの平行移動に対して内積の大小は変わらないため、起点をどの位置にとってもよい。また、ここでは単位ベクトルｎ１を中心結線ベクトルに平行にとったが、ベクトルｎ１の向きはこれに限られない。ベクトルｎ１が中心結線ベクトルに対して直交していなければ（つまり、平面６１，６２が粒子ｉ，ｊの中心を結んだ直線に交わっていれば）、同様の処理が可能である。
【００４３】
以上述べた方法により接触可能頂点を選択し、判定対象とする粒子ｉの辺と粒子ｊの面を絞り込むことで、接触判定を行う辺、面の組み合わせを減らし、計算時間を短縮することができる。
【００４４】
（第２の抽出方法）
上述した第１の抽出方法では、接触判定を行う範囲を平面６１と６２で挟まれた空間に限定している。しかし、図６で示されるように、平面６１と６２のあいだに、実際には他の粒子と接触していない頂点が多く含まれる場合もある。そこで、第２の抽出方法として、接触可能頂点をさらに絞り込む方法について説明する。
【００４５】
第２の抽出方法は、図７に示すように、平面６１，６２とは異なる角度で平面７１，７２を設定し、４つの平面６１，６２，７１，７２で挟まれた空間に位置する粒子ｉの辺及
び粒子ｊの面のみを判定対象として抽出するものである。ここで、平面７１（第３の平面）と平面７２（第４の平面）はそれぞれ粒子ｊ，ｉに接し、且つ、互いに平行になるように設定される。言い換えると、第２の抽出方法は、平面６１及び平面７１よりも粒子ｊの中心側の空間に位置する粒子ｉの辺と、平面６２及び平面７２よりも粒子ｉの中心側の空間に位置する粒子ｊの面とを、判定対象として抽出する方法である。
【００４６】
第２の抽出方法を実現する具体的なアルゴリズムを説明する。図７は２つのベクトルによる接触可能頂点の絞り込み方法を示す図である。図７左上は上記の中心結線単位ベクトルｎ１によって絞り込まれた接触可能頂点を示している。塗りつぶし（黒丸）で示されているのが接触可能頂点である。
【００４７】
図７左下は、重なり最小単位ベクトルｎ２によって絞り込まれた接触可能頂点を表している。「重なり最小単位ベクトル」とは、粒子ｉと粒子ｊの重なりを最小とする向きの単位ベクトルをいい、以降の説明では単に「重なり最小ベクトル」とも称す。ｎ２による接触可能頂点の絞り込み方法はｎ１による接触可能頂点の絞り込み方法と同じである。つまり、ベクトルｎ２と粒子ｉの各頂点座標の位置ベクトルＰｉｋの内積ｎ２・Ｐｉｋの最大値Ｄ２ｍａｘと、ベクトルｎ２と粒子ｊの各頂点座標の位置ベクトルＰｊｋの内積ｎ２・Ｐｊｋの最小値Ｄ２ｍｉｎを求める。そして、内積がＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる頂点を接触可能頂点として抽出する。
【００４８】
図７右は中心結線ベクトルによる接触可能頂点集合をＡ、重なり最小ベクトルによる接触可能頂点集合をＢとしたときに、共通部分である積集合Ａ∩Ｂで選択される頂点を表している。図７において、集合Ａに含まれる頂点数は１０、集合Ａ内の頂点を構成頂点とする辺数は１２、集合Ｂに含まれる頂点数は８、辺数は１０である。一方、積集合Ａ∩Ｂに含まれる頂点数は４、辺数は６となり、１つのベクトルのみの場合より判定対象を絞り込むことができる。よって、積集合Ａ∩Ｂに含まれる頂点を最終的な接触可能頂点と定義する。
【００４９】
積集合Ａ∩Ｂが空集合でも接触する場合がある。図８を用いて説明する。四角で表されている頂点が中心結線ベクトルにより選択された頂点（集合Ａ）、三角で表されているのが重なり最小ベクトルにより選択された頂点（集合Ｂ）である。塗りつぶしの円で表されているのが積集合Ａ∩Ｂに含まれる接触可能頂点である。粒子ｊは接触可能頂点を持たないが、図８に示すように粒子ｉと接触している。このような場合を含めるため、判定対象とする辺、面は次のように決定する。
・積集合Ａ∩Ｂに含まれる頂点（接触可能頂点）を持つ辺、面。
・集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を持つ辺、面。
【００５０】
以上の２パターンによって、全ての接触の可能性がある辺と面を選択することができる。本方法により、接触判定を行う辺、面の組み合わせを大幅に減らすことができ、計算時間のより一層の短縮が可能となる。
【００５１】
本方法で用いる重なり最小ベクトルｎ２の算出方法について説明する。図９は重なり最小ベクトルを求めるフローチャートである。図９に示す処理は、接触判定部Ｂ１２０が実行するものである。
【００５２】
まず、重なり最小ベクトルの初期値を設定する（ステップＳ１）。初期値は中心結線単位ベクトルと平行な単位ベクトルとする。粒子ｉの中心位置ベクトルをＣｉ、粒子ｊの中心位置ベクトルをＣｊとすると、重なり最小ベクトルｎ２の初期値は下記式で表される。
【数７】
以降の処理で、ベクトルｎ２の向きを微小に変動させながら、粒子ｉと粒子ｊの重なりが最小となるものを見つける。その処理で用いる微小量δの初期値δｍａｘ（＝０．１程度）をステップＳ１で設定する。
【００５３】
次に、２粒子の重なり量Ｌ（＝Ｄ２ｍａｘ−Ｄ２ｍｉｎ）を算出する（ステップＳ２）。
次に、接触の可能性があるか否かを判定する（ステップＳ３）。Ｌが０以下の場合、２粒子は非接触であるので、接触判定から除外される（ステップＳ４）。Ｌが０より大きい場合は接触の可能性があるため、重なり最小ベクトルｎ２を更新（探索）するためのループに入る。
【００５４】
重なり最小ベクトルｎ２に直交し、たがいに直交するベクトルｕ、ｖを適当に定める（ステップＳ５）。そして、試行する重なり最小ベクトルｎ２（１）〜ｎ２（４）を下記式のように設定する（ステップＳ６）。
【数８】
上記式の１行目で決定される試行重なり最小ベクトルｎ２（１）は、図９中に示されるように、直交ベクトルｕの方向へ微小量δで決まる量だけ回転させたものとなる。
【００５５】
このようにして決定された各試行重なり最小ベクトルｎ２（ｋ）に対して２粒子の重なり量Ｌｎｅｗ（ｋ）を新たに算出する（ｋ＝１，２，３，４）。算出された４つのＬｎｅｗ（ｋ）の最小値とＬを比較する（ステップＳ７）。ｍｉｎ｛Ｌｎｅｗ（ｋ）｝がＬより小さければ、重なり最小ベクトルｎ２をｎ２（ｋ）に更新する。また、ＬもＬｎｅｗ（ｋ）に更新する（ステップＳ８）。重なり量Ｌを最小とするループにおいて、ｍｉｎ｛Ｌｎｅｗ（ｋ）｝がＬより大きくなった場合は（ステップＳ９）、微小量δをδ／２と小さくした後、試行重なりベクトルの設定から繰り返す（ステップＳ１０）。
【００５６】
上記のループを繰り返して、δがδｍｉｎ（＝０．００１程度）より小さくなったときのｎ２が最終的な重なり最小ベクトルとなる（ステップＳ１１）。このような処理により、粒子ｉとｊの重なり量Ｌを最も小さくするベクトルｎ２を得ることができる。
【００５７】
なお、図９のステップＳ５で示したように、ループを繰り返すたびにｕ，ｖを面内で回転させると、局所解への陥りを抑制し、重なり最小ベクトルの最適解を得やすいという効果が得られる。
【００５８】
図９で示した方法は重なり最小ベクトルｎ２を決定するための一手法にすぎない。他の公知の最適化手法を利用して、重なり量Ｌを最小とするベクトルｎ２を求めることも可能である。また、少なくともベクトルｎ１とｎ２の向きが異なっていれば、第１の抽出方法よりも判定対象を絞り込む効果は得られるため、ベクトルｎ１に対して任意の回転を行ったものをベクトルｎ２として用いるだけでもよい。
【００５９】
（接触判定処理）
接触判定部Ｂ１２０の処理フローを図１０で説明する。
まず、接触判定を行う２粒子ｉ、ｊを設定する（ステップＳ１００）。次に２粒子ｉ、ｊの中心座標から中心結線単位ベクトルｎ１を求める（ステップＳ１１０）。ｎ１を基に２粒子ｉ、ｊの重なり量Ｌを算出する（ステップＳ１２０）。重なり量Ｌを基に接触可能性を判定し（ステップＳ１３０）、Ｌ＞０でなければ、２粒子は非接触として接触判定を終了する（ステップＳ１３５）。次に、２粒子ｉ、ｊの位置座標情報から重なり最小単位ベクトルｎ２を求める（ステップＳ１４０）。ステップＳ１４０の処理内容は、図９で詳しく説明したとおりである。
【００６０】
ここで、ステップＳ１１０で設定したｎ１を用いて接触可能頂点集合Ａを設定し（ステップＳ１５０）、ステップＳ１４０で設定したｎ２を用いて接触可能頂点集合Ｂを設定する（ステップＳ１６０）。次に、接触可能頂点集合Ａ∩Ｂを設定する（ステップＳ１７０）。こうして設定した接触可能頂点を含む粒子ｉの辺と粒子ｊの面を抽出する（ステップＳ１８０）。以上で、接触判定に用いる判定対象が絞り込まれる。
【００６１】
次に、粒子ｉが粒子ｊに接触しているか否かの判定処理を行う。まず、接触判定を行う粒子ｉの辺を設定する（ステップＳ１９０）。接触判定を行う粒子ｊの面を設定する（ステップＳ２００）。ここで設定した粒子ｉの辺と粒子ｊの面との交差の有無を判定する（ステップＳ２１０）。ここで、粒子ｊの面に粒子ｉの辺が交差していれば、交差回数のカウントアップを行う（ステップＳ２２０）。この操作をステップＳ１８０で抽出された粒子ｊの面の全てに対して行う（ステップＳ２３０）。次に、粒子ｉの辺と粒子ｊの面の交差回数が２回かどうかの判定を行う（ステップＳ２４０）。交差回数が２回ならば、めり込み量を算出し、メモリに記憶しておく（ステップＳ２５０）。本操作をＳ１８０で抽出された粒子ｉの辺の全てに対して行う（ステップＳ２６０）。
【００６２】
なお、図１０では、説明を簡単にするため、粒子ｉのある辺が粒子ｊの２つの面と交差した場合に、接触有りと判定し、粒子ｉとｊのあいだのめり込み量を算出する処理例を示した。しかし実際の処理では、図４に示したように、粒子ｉと粒子ｊの接触を４つのパターンで判定する。前述したようにパターンごとに判定条件が異なるが、いずれのパターンにおいても粒子ｉの辺と粒子ｊの面との交差判定が基本となるので、図１０のステップＳ１８０で抽出した辺及び面を利用すればよい。
【００６３】
各パターン（ａ）〜（ｄ）におけるめり込み量の定義について、図４中の下段の図を用いて説明する。
（ａ）頂点のめり込み：粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０との距離をめり込み量とする。
（ｂ）辺のめり込み（ケース１）：粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｃ）辺のめり込み（ケース２）：粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｄ）辺のめり込み（ケース３）：粒子ｉの辺ｉ−ｅ０に対する、粒子ｊの辺ｊ−ｅ０、ｊ−ｅ１との距離をそれぞれの辺とのめり込み量とする。
【００６４】
（粒子受力計算部Ｂ１３０）
以上の４パターンにおいて、めり込み量（距離）の計算に用いた各粒子表面の２点を結ぶ方向を接触面法線方向と定める。この接触面法線方向は、接触点の（粒子ｊの面又は辺に対する）めり込み方向ということもできる。そして、接触判定部Ｂ１２０で求めためり込み量に基づいて、ヘルツの式を用いて、粒子ｉに働く接触面法線方向の接触力（バネ力）を求める。また各接触点に対して、接触が始まってからの接触面（上記接触面法線方向に垂直な面）内方向の移動量を求め、ミンドリンの式を用いて粒子ｉに働く接触面内方向の接触力（バネ力）を求める。また接触面法線方向の運動速度（めり込み速度）と面内方向の各運動速度（ずり速度）から、各方向に作用する粘性力（ダンパー力）を求める。そしてそれら接触力と粘性力の和に、粒子ｉに働くその他の力（重力、付着力、電磁力、空気の粘性抵抗等）を加えて、各粒子ｉが受ける並進力と回転力を求める。なお、本めり込み量から接触力を求める過程は、上記非特許文献１に記載されている手法と同じであるため、詳しい説明は割愛する。
【００６５】
（粒子運動計算部Ｂ１４０）
次に、粒子運動計算部Ｂ１４０における各粒子の刻み時間後の速度と位置の算出方法を説明する。この処理は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。すなわち、粒子運動計算部Ｂ１４０は、並進力からニュートンの運動方程式を用いて並進方向の加速度を求め、更に速度と位置を求める。また粒子運動計算部Ｂ１４０は、回転力からオイラーの運動方程式を用いて角加速度を求め、更に回転速度と回転位置を求める。本処理は、非特許文献１に記載されている手法と同じであるため、詳しい説明は割愛する。
【００６６】
（実施例）
以上の接触判定方法を用いて粒子挙動計算を実施した例を、図１１〜図１２を用いて説明する。図１１（ａ）は解析モデルを示す。解析モデルは電子写真技術を用いた画像形成装置におけるクリーニングプロセスである。０．１ｍ／ｓｅｃで矢印の方向に移動する感光体ドラムに対してクリーニングブレードをβ＝３０°で当接させておき、感光体ドラム上に乗って運ばれてきたトナーの運動を見る。各トナーの形状を粒子形状データとして与え、感光体ドラム及びクリーニングブレードの形状を壁形状データとして与える。
【００６７】
図１１（ｂ）は、解析条件として与えた、トナー、感光体ドラム、クリーニングブレードそれぞれの機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）を示している。この条件で、刻み時間Δｔを１０ｎｓｅｃに設定し、約２百万ステップの計算を実行した。
【００６８】
図１２は、（ａ）球形粒子と（ｂ）異形粒子の計算結果を比較したものである。球形粒子の計算は従来の計算方法で実施している。異形粒子は４面体の角をとった形状であり、その表面を８つの３角形と６つの４角形の合計１２の多角形で構成している。なお本計算においては、トナーがブレードに堰き止められた際の挙動をみるため、強制的にトナーがすり抜けないようにしている。球形粒子に対して異形粒子は回転しにくいことから、全体的に動き辛く、トナー粒子の密集度が上がる傾向にある。図１２の計算結果では、異形粒子の密集度が高くなっている様子が再現されている。
【００６９】
このトナーがブレードに堰き止められた際の挙動は、クリーニング不良であるトナーのすり抜けとの関係があり、流動的である程すり抜けやすいといわれている。このことから、製品設計を行う上で、トナーの形状によるクリーニング不良の傾向を粒子挙動計算で予測することが可能となる。
【００７０】
また、図１３は面数の異なる異形粒子の自由落下の計算を示した図である。（ａ）は、
頂点数１２、面数１４、辺数２４の場合の異形粒子で、（ｂ）は頂点数９８、面数１０８、辺数２０４の異形粒子の場合である。（ａ）と（ｂ）の場合において、２５６個の異形粒子の自由落下の計算を刻み時間Δｔ＝１０ｎｓｅｃで８百万ステップ実行し、本発明を使用しない場合と使用した場合での計算時間の比較を行った。（ａ）では、従来の方法では計算に１９時間３５分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込み（第１の抽出方法）を用いた場合では、８時間０９分となり、高速化率は２．４倍となる。（ｂ）では、従来の方法では計算に５２４時間３２分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込みを用いると２２時間２４分となり、高速化率は２３．４倍となった。このことから、粒子の面の数が多くなればなるほど、本発明を用いることによる計算時間の短縮効果が大きくなることがわかる。
【００７１】
図１４は中心結線ベクトルと重なり最小ベクトルを用いた高速化の効果を示す図である。計算には頂点数１９８、面数３９２、辺数５８８の異形粒子を用い、図１３と同様の自由落下計算を３０００００ステップ、Δｔ＝１０ｎｓｅｃで行った。接触判定に中心結線ベクトルのみを用いた場合（第１の抽出方法）と中心結線ベクトルと重なり最小ベクトルの両方を用いた場合（第２の抽出方法）で計算時間を比較した。中心結線ベクトルのみを用いた場合は４７２６秒、中心結線ベクトルと重なり最小ベクトルの両方を用いた場合は２００９秒となり、高速化率は２．３５倍となった。このことから、中心結線ベクトルと重なり最小ベクトルの両方を用いた接触判定により、さらなる計算時間の短縮が実現できていることが分かる。
【符号の説明】
【００７２】
１０：情報処理装置
Ｂ１００：全体制御部
Ｂ１１０：入力データ読込部
Ｂ１２０：接触判定部
Ｂ１３０：粒子受力計算部
Ｂ１４０：粒子運動計算部
Ｂ１５０：計算結果出力部
【特許請求の範囲】
【請求項１】
情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、
情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、
情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、
情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、
情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、
前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析方法。
【請求項２】
粒子ｉの辺及び粒子ｊの面を抽出するステップは、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出するステップである
ことを特徴とする請求項１に記載の粒子挙動解析方法。
【請求項３】
情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、
情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、
情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、
情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、
情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、
前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析方法。
【請求項４】
粒子ｉの辺及び粒子ｊの面を抽出するステップは、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
前記第３の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ２を決定し、
単位ベクトルｎ２と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ２・Ｐｉを計算し、
単位ベクトルｎ２と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ２・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ２・Ｐｉのうち、最大となるものをＤ２ｍａｘ、粒子ｊの全ての頂点についての内積ｎ２・Ｐｊのうち、最小となるものをＤ２ｍｉｎとした場合に、内積ｎ２・ＰｉがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ２・ＰｊがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｊの頂点からなる集合Ｂを定義し、
集合Ａと集合Ｂの積集合に含まれる頂点を含む粒子ｉの辺及び粒子ｊの面と、集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を含む粒子ｉの辺及び粒子ｊの面とを抽出するステップである
ことを特徴とする請求項３に記載の粒子挙動解析方法。
【請求項５】
前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定する
ことを特徴とする請求項２又は４に記載の粒子挙動解析方法。
【請求項６】
前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定し、
前記単位ベクトルｎ２を、Ｄ２ｍａｘ−Ｄ２ｍｉｎがＤ１ｍａｘ−Ｄ１ｍｉｎよりも小さくなるように決定する
ことを特徴とする請求項４に記載の粒子挙動解析方法。
【請求項７】
密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、
複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、
前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、
前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、
前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、
前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析装置。
【請求項８】
前記接触判定部は、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出する
ことを特徴とする請求項７に記載の粒子挙動解析装置。
【請求項９】
密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、
複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、
前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、
前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、
前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、
前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析装置。
【請求項１０】
前記接触判定部は、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
前記第３の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ２を決定し、
単位ベクトルｎ２と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ２・Ｐｉを計算し、
単位ベクトルｎ２と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ２・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ２・Ｐｉのうち、最大となるものをＤ２ｍａｘ、粒子ｊの全ての頂点についての内積ｎ２・Ｐｊのうち、最小となるものをＤ２ｍｉｎとした場合に、内積ｎ２・ＰｉがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ２・ＰｊがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｊの頂点からなる集合Ｂを定義し、
集合Ａと集合Ｂの積集合に含まれる頂点を含む粒子ｉの辺及び粒子ｊの面と、集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を含む粒子ｉの辺及び粒子ｊの面とを抽出する
ことを特徴とする請求項９に記載の粒子挙動解析装置。
【請求項１１】
前記接触判定部は、前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定する
ことを特徴とする請求項８又は１０に記載の粒子挙動解析装置。
【請求項１２】
前記接触判定部は、
前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定し、
前記単位ベクトルｎ２を、Ｄ２ｍａｘ−Ｄ２ｍｉｎがＤ１ｍａｘ−Ｄ１ｍｉｎよりも小さくなるように決定する
ことを特徴とする請求項１０に記載の粒子挙動解析装置。
【請求項１３】
請求項１〜６のうちいずれか１項に記載の粒子挙動解析方法の各ステップを情報処理装置に実行させることを特徴とする解析プログラム。
【請求項１】
情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、
情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、
情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、
情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、
情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、
前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析方法。
【請求項２】
粒子ｉの辺及び粒子ｊの面を抽出するステップは、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出するステップである
ことを特徴とする請求項１に記載の粒子挙動解析方法。
【請求項３】
情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、
情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、
情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、
情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、
情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、
前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析方法。
【請求項４】
粒子ｉの辺及び粒子ｊの面を抽出するステップは、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
前記第３の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ２を決定し、
単位ベクトルｎ２と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ２・Ｐｉを計算し、
単位ベクトルｎ２と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ２・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ２・Ｐｉのうち、最大となるものをＤ２ｍａｘ、粒子ｊの全ての頂点についての内積ｎ２・Ｐｊのうち、最小となるものをＤ２ｍｉｎとした場合に、内積ｎ２・ＰｉがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ２・ＰｊがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｊの頂点からなる集合Ｂを定義し、
集合Ａと集合Ｂの積集合に含まれる頂点を含む粒子ｉの辺及び粒子ｊの面と、集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を含む粒子ｉの辺及び粒子ｊの面とを抽出するステップである
ことを特徴とする請求項３に記載の粒子挙動解析方法。
【請求項５】
前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定する
ことを特徴とする請求項２又は４に記載の粒子挙動解析方法。
【請求項６】
前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定し、
前記単位ベクトルｎ２を、Ｄ２ｍａｘ−Ｄ２ｍｉｎがＤ１ｍａｘ−Ｄ１ｍｉｎよりも小さくなるように決定する
ことを特徴とする請求項４に記載の粒子挙動解析方法。
【請求項７】
密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、
複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、
前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、
前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、
前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、
前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析装置。
【請求項８】
前記接触判定部は、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出する
ことを特徴とする請求項７に記載の粒子挙動解析装置。
【請求項９】
密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、
複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、
前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、
前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、
前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、
前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析装置。
【請求項１０】
前記接触判定部は、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
前記第３の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ２を決定し、
単位ベクトルｎ２と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ２・Ｐｉを計算し、
単位ベクトルｎ２と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ２・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ２・Ｐｉのうち、最大となるものをＤ２ｍａｘ、粒子ｊの全ての頂点についての内積ｎ２・Ｐｊのうち、最小となるものをＤ２ｍｉｎとした場合に、内積ｎ２・ＰｉがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ２・ＰｊがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｊの頂点からなる集合Ｂを定義し、
集合Ａと集合Ｂの積集合に含まれる頂点を含む粒子ｉの辺及び粒子ｊの面と、集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を含む粒子ｉの辺及び粒子ｊの面とを抽出する
ことを特徴とする請求項９に記載の粒子挙動解析装置。
【請求項１１】
前記接触判定部は、前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定する
ことを特徴とする請求項８又は１０に記載の粒子挙動解析装置。
【請求項１２】
前記接触判定部は、
前記単位ベクトルｎ１を、粒子ｉの中心と粒子ｊの中心を結ぶ直線に平行となるように決定し、
前記単位ベクトルｎ２を、Ｄ２ｍａｘ−Ｄ２ｍｉｎがＤ１ｍａｘ−Ｄ１ｍｉｎよりも小さくなるように決定する
ことを特徴とする請求項１０に記載の粒子挙動解析装置。
【請求項１３】
請求項１〜６のうちいずれか１項に記載の粒子挙動解析方法の各ステップを情報処理装置に実行させることを特徴とする解析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【公開番号】特開２０１３−８９１００（Ｐ２０１３−８９１００Ａ）
【公開日】平成２５年５月１３日（２０１３．５．１３）
【国際特許分類】
【出願番号】特願２０１１−２３０３０４（Ｐ２０１１−２３０３０４）
【出願日】平成２３年１０月２０日（２０１１．１０．２０）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【公開日】平成２５年５月１３日（２０１３．５．１３）
【国際特許分類】
【出願日】平成２３年１０月２０日（２０１１．１０．２０）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
[ Back to top ]