(57)【特許請求の範囲】
【請求項１】
読み書き可能なメモリを有する情報処理装置において、密集状態にある複数の粒子の運動を計算する方法であって、
情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを前記メモリに読み込む入力ステップと、
情報処理装置が、前記複数の粒子それぞれについて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、
情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、
情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、
前記接触判定ステップにおいて、粒子ｉの頂点ｉ−ｖ０が他の粒子ｊにめり込んでいると判定された場合に、
前記計算ステップでは、
粒子ｊの複数の面のそれぞれについて、頂点ｉ−ｖ０の当該面に対するめり込み量と、頂点ｉ−ｖ０を端点にもつ粒子ｉのすべての辺それぞれの当該面に対するめり込み量と、の総和を算出し、
めり込み量の総和が最小となる粒子ｊの面を、めり込み基準面と定め、
前記めり込み基準面に対するめり込み量に基づいて、粒子ｉに作用する接触力を算出する
ことを特徴とする粒子挙動解析方法。
【請求項２】
前記計算ステップでは、粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０とのあいだの距離を、頂点ｉ−ｖ０の面ｊ−ｆ０に対するめり込み量とする
ことを特徴とする請求項１に記載の粒子挙動解析方法。
【請求項３】
前記計算ステップでは、粒子ｉの頂点ｉ−ｖ０を端点にもつ辺ｉ−ｅ０が粒子ｊの面ｊ−ｆ０とは異なる面ｊ−ｆ１に交差する点と、面ｊ−ｆ０とのあいだの距離を、辺ｉ−ｅ０の面ｊ−ｆ０に対するめり込み量とする
ことを特徴とする請求項２に記載の粒子挙動解析方法。
【請求項４】
頂点ｉ−ｖ０を端点にもつ辺の他方の端点である頂点ｉ−ｖ１も粒子ｊにめり込んでいる場合に、前記計算ステップでは、頂点ｉ−ｖ１のめり込み量も算出し、めり込み量の総和に加算する
ことを特徴とする請求項１〜３のうちいずれか１項に記載の粒子挙動解析方法。
【請求項５】
密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、
複数の粒子それぞれの形状を多面体により定義するデータを前記メモリに読み込む読込部と、
前記複数の粒子それぞれについて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、
前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、
前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、
前記接触判定部において、粒子ｉの頂点ｉ−ｖ０が他の粒子ｊにめり込んでいると判定された場合に、
前記計算部では、
粒子ｊの複数の面のそれぞれについて、頂点ｉ−ｖ０の当該面に対するめり込み量と、頂点ｉ−ｖ０を端点にもつ粒子ｉのすべての辺それぞれの当該面に対するめり込み量と
、の総和を算出し、
めり込み量の総和が最小となる粒子ｊの面を、めり込み基準面と定め、
前記めり込み基準面に対するめり込み量に基づいて、粒子ｉに作用する接触力を算出する
ことを特徴とする粒子挙動解析装置。
【請求項６】
前記計算部では、粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０とのあいだの距離を、頂点ｉ−ｖ０の面ｊ−ｆ０に対するめり込み量とする
ことを特徴とする請求項５に記載の粒子挙動解析装置。
【請求項７】
前記計算部では、粒子ｉの頂点ｉ−ｖ０を端点にもつ辺ｉ−ｅ０が粒子ｊの面ｊ−ｆ０とは異なる面ｊ−ｆ１に交差する点と、面ｊ−ｆ０とのあいだの距離を、辺ｉ−ｅ０の面ｊ−ｆ０に対するめり込み量とする
ことを特徴とする請求項６に記載の粒子挙動解析装置。
【請求項８】
頂点ｉ−ｖ０を端点にもつ辺の他方の端点である頂点ｉ−ｖ１も粒子ｊにめり込んでいる場合に、前記計算部では、頂点ｉ−ｖ１のめり込み量も算出し、めり込み量の総和に加算する
ことを特徴とする請求項５〜７のうちいずれか１項に記載の粒子挙動解析装置。
【請求項９】
請求項１〜４のうちいずれか１項に記載の粒子挙動解析方法の各ステップを情報処理装置に実行させることを特徴とする解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、密集状態にある多数の粒子の挙動を計算する方法に関するものであり、更には個々の粒子の形状を考慮して挙動を計算するものである。
【背景技術】
【０００２】
粉体や粒体などの粒子を取り扱う分野では、個々の粒子の挙動を知ることが重要である。従来これらの分野では、実験と観察から粒子の挙動を把握することが多かった。しかし近年、計算機内でその運動を再現させ、その挙動を把握する、いわゆるシミュレーション技術の適用が活発に行われている。正確な粒子挙動のシミュレーションを行うには、粒子のもつ形状、機械的特性等の情報をできるだけ詳細に計算モデルに反映させ、それらの因子を考慮して解を求めることが重要である。
【０００３】
一方、このような粒子を扱う装置の１つに、プリンタ、複写機、ファクシミリ等の電子写真技術を用いた画像形成装置がある。これらの装置では、直径数ミクロンの粉体粒子からなるトナーを用いて、帯電、露光、現像、転写、定着、クリーニングという６つのプロセスを経て、紙上に画像を形成する。このような画像形成装置では、トナー粒子の形状が、作成される画像に大きな影響を与えることがわかっている。そこで画像形成装置の開発においても、シミュレーション技術を用いて、密集状態にある個々のトナー粒子の挙動を求める試みが活発に行われ、近年の傾向として特にトナー粒子の形状を考慮することが求められている。
【０００４】
密集状態にある粒子の挙動をシミュレーションする方法として個別要素法（ＤＥＭ）がある。個別要素法の技術は非特許文献１に詳細に記載されている。本文献によれば、個々の粒子を球形で表現し、それらの接触に際してはそのめり込み量からヘルツの式とミンドリンの式を用いることで粒子の持つ機械特性を考慮した接触力を求める。そしてニュートンの運動方程式（並進方向）とオイラーの運動方程式（回転方向）を用いて、その接触力をもとに個々の粒子の時間進展に伴う運動を求める。
【０００５】
個別要素法の特徴は個々の粒子を要素と捉え、それらが接触する部分にバネ、ダンパー、スライダを仮定して、その接触力を求め、運動方程式を解くところにある。そこで、粒子の形状を球形に限定することなく考慮してその挙動を計算するものとして、非特許文献２、非特許文献３がある。これらの文献では粒子の形状を多面体で表現し離散要素法を用いてその挙動を計算している。そして、粒子の接触の検出に関して、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)
【非特許文献４】酒井幸市、OpenGLで作る力学アニメーション入門、森北出版(2005)
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明の目的は、密集状態にある多数の粒子の挙動を、その形状も考慮して、安定に計算可能にすることにある。また本発明の他の目的は、上記従来技術で用いているCommon-Planeを使用することなく、粒子同士で作用する接触力を適切に計算することにある。
【課題を解決するための手段】
【０００９】
本発明の第一態様は、読み書き可能なメモリを有する情報処理装置において、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを前記メモリに読み込む入力ステップと、情報処理装置が、前記複数の粒子それぞれについて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子ｉの頂点ｉ−ｖ０が他の粒子ｊにめり込んでいると判定された場合に、前記計算ステップでは、粒子ｊの複数の面のそれぞれについて、頂点ｉ−ｖ０の当該面に対するめり込み量と、頂点ｉ−ｖ０を端点にもつ粒子ｉのすべての辺それぞれの当該面に対するめり込み量と、の総和を算出し、めり込み量の総和が最小となる粒子ｊの面を、めり込み基準面と定め、前記めり込み基準面に対するめり込み量に基づいて、粒子ｉに作用する接触力を算出することを特徴とする粒子挙動解析方法を提供するものである。
【００１０】
本発明の第二態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを前記メモリに読み込む読込部と、前記複数の粒子それぞれについて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部において、粒子ｉの頂点ｉ−ｖ０が他の粒子ｊにめり込んでいると判定された場合に、前記計算部では、粒子ｊの複数の面のそれぞれについて、頂点ｉ−ｖ０の当該面に対するめり込み量と、頂点ｉ−ｖ０を端点にもつ粒子ｉのすべての辺それぞれの当該面に対するめり込み量と、の総和を算出し、めり込み量の総和が最小となる粒子ｊの面を、めり込み基準面と定め、前記めり込み基準面に対するめり込み量に基づいて、粒子ｉに作用する接触力を算出することを特徴とする粒子挙動解析装置を提供するものである。
【００１１】
本発明の第三態様は、上記した本発明に係る粒子挙動解析方法の各ステップを情報処理装置に実行させることを特徴とする解析プログラムを提供するものである。
【発明の効果】
【００１２】
本発明によれば、密集状態にある多数の粒子の挙動を、その形状も考慮して、安定に計
算可能にすることができる。
【図面の簡単な説明】
【００１３】
【図１】粒子挙動解析装置（情報処理装置）のハードウエア構成を示すブロック図。
【図２】粒子挙動解析装置の機能構成を示すブロック図。
【図３】接触パターンと各パターンにおける粒子めり込み量を示す図。
【図４】誤っためり込み基準面を選択したときの不具合を説明する図。
【図５】本発明におけるめり込み基準面を選択する方法を説明する図。
【図６】本発明を用いて選択されためり込み基準面を説明する図。
【図７】めり込み頂点辺の他端の頂点もめり込んでいる場合の処理を説明する図。
【図８】粒子受力計算の流れを示すフローチャート。
【図９】実施例としての計算条件を示す図。
【図１０】実施例の実験結果を示す図。
【発明を実施するための形態】
【００１４】
本発明の実施形態は、密集状態にある多数の粒子の挙動をコンピュータ・シミュレーションにより計算するための粒子挙動解析方法及び装置を提供するものである。本実施形態に係る粒子挙動解析方法の特徴の一つは、多数の多角形から構成される多面体（ポリゴンモデル）により物体（粒子、壁など）の形状を定義した点である。これにより任意の粒子形状や壁形状を表現できるようになり、多様なケースのシミュレーションを精度良く行うことが可能になる。本方法のもう一つの特徴は、粒子が他の物体から受ける接触力を計算するにあたり、当該粒子のめり込み量やめり込み方向などを算出する基準となる面（めり込み基準面）を、適切かつ安定したシミュレーション結果が得られるように設定する点にある。以下、本実施形態の具体的な構成例及び特徴について図に基づき説明する。
【００１５】
（解析装置のハードウエア構成）
図１は、本実施形態における粒子挙動解析装置が動作する計算装置を示すブロック図である。計算装置本体（情報処理装置又はコンピュータとも呼ぶ）２０は、各種判断及び処理を行う中央処理装置（ＣＰＵ）２１と、処理プログラム及び固定データを格納するＲＯＭ２２と、処理データを記憶するＲＡＭ２３と、入出力回路（Ｉ／Ｏ）２４とを有している。また必要に応じて、計算装置本体２０には、Ｉ／Ｏ２４を介して記憶装置、入力装置、表示装置、他のコンピュータなどが接続される。この計算装置本体２０には、入力データ３０がＩ／Ｏ２４を介して入力される。計算結果は、Ｉ／Ｏ２４を介して出力データ３１として出力される。
【００１６】
（解析装置の機能）
図２は、粒子挙動解析を行う解析プログラムによって実現される機能の構成を示すブロック図である。解析プログラムは、ＲＯＭ２２や記憶装置などのコンピュータ読み取り可能な記録媒体に記憶されており、必要に応じて、ＣＰＵ２１によって読み込まれ実行されるソフトウエアである。この解析プログラムは、主な機能（モジュール）として、入力データ読込部（取得部）Ｂ１１０、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、粒子運動計算部Ｂ１４０、計算結果出力部Ｂ１５０を備えている。また粒子受力計算部Ｂ１３０は、めり込み基準面設定部Ｂ１３１とめり込み量算出部Ｂ１３２を有している。全体制御部Ｂ１００は、これらのモジュールの処理の順番やデータの引渡しなどの全体制御を担っている。解析プログラムは、読み書き可能なメモリであるＲＡＭ２３を、計算で用いる各種のデータを一時的に格納するワークメモリとして利用する。
【００１７】
入力データ読込部Ｂ１１０は、粒子形状データ、壁形状データ、及び各種パラメータを読み込んでＲＡＭ２３に格納する。粒子形状データとは、粒子の表面を複数の多角形（ポリゴン）で表現したものである。ここで多角形としてはその頂点が一平面上にあることが
望ましく、３角形を用いるのが最も単純である。また粒子は凸形状（凸多面体）に限定すると、他の粒子や壁との接触状態の判定が簡便になる。壁形状データも同様であり、壁の表面を複数の多角形で表現したものである。粒子形状データ、壁形状データともに、多角形の頂点座標と、各多角形を構成する頂点の番号とで構成されるデータである。各種パラメータには、粒子と壁の機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）、粒子に働く力（重力、付着力、電磁力、空気の粘性特性等）、及び計算条件としての計算刻み時間、計算終了時刻がある。なお、壁とは、粒子が移動可能な空間を形成する（規定する）境界面である。
【００１８】
接触判定部Ｂ１２０は、粒子の位置データに基づいて、粒子が周囲の粒子または壁と接触しているかどうか、またどのように接触しているかを判定する。図３を用いて接触判定の内容を示す。図３の上段は４種類の接触パターンを３次元的に示したものであり、下段は各パターンにおける２つの粒子の接触によるめり込みの様子を示すものである。図３において、実線は６面体粒子ｉ、破線は６面体粒子ｊを示す。ｉ−ｖ０は粒子ｉを構成する頂点、ｉ−ｅ０は粒子ｉを構成する辺、ｊ−ｅ０とｊ−ｅ１は粒子ｊを構成する辺、ｊ−ｆ０とｊ−ｆ１は粒子ｊを構成する面、星印は粒子ｉと粒子ｊの接触点（以下、めり込み点と記す）を示す。また図３の上段における矢印は粒子ｉが粒子ｊに衝突する方向を示し、また下段における黒丸は粒子ｊの辺を示す。粒子が凸形状である場合、粒子ｉと粒子ｊの接触は、図３に示す以下の（ａ）〜（ｄ）の４つのパターンのいずれかとなる。
【００１９】
（ａ）頂点のめり込み：
粒子ｉを構成する頂点が粒子ｊの１つの面にめり込む場合である。粒子ｉの頂点が粒子ｊの内部にあれば、接触していると判定する。図３（ａ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込んでいる。このパターンでは、頂点ｉ−ｖ０が、粒子ｉの粒子ｊに対するめり込み点である。
【００２０】
（ｂ）辺のめり込み（ケース１）：
上記（ａ）に加えて、めり込んだ頂点を端点に持つ粒子ｉの辺が粒子ｊの辺にもめり込む場合である。（ａ）のめり込んだ頂点を端点に持つ粒子ｉの辺（めりこみ頂点辺と呼ぶ）において、粒子ｊ構成面との交点数が１であり、かつその交点を持つ面が（ａ）におけるめり込み面でなければ、これに該当する。図３（ｂ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込む（上記（ａ）に該当）とともに、めり込み頂点辺ｉ−ｅ０が粒子ｊの面ｊ−ｆ１と交差しており、辺ｊ−ｅ０にもめり込んでいる。このパターンでは、頂点ｉ−ｖ０に加え、めり込み頂点辺ｉ−ｅ０上にもう一つめり込み点を考える。
【００２１】
（ｃ）辺のめり込み（ケース２）：
粒子ｉの辺が粒子ｊの１つの辺だけにめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、その交点を持つ２つの面が共有辺を持ち、かつその共有辺が粒子ｉの辺に近接していれば、これに該当する。図３（ｃ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にめり込んでいる。このパターンでは、辺ｊ−ｅ０に近接する辺ｉ−ｅ０上の点を、めり込み点として考える。
【００２２】
（ｄ）辺のめり込み（ケース３）：
粒子ｉの辺が粒子ｊの２つの辺にめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、（ｃ）ではない場合がこれに該当する。図３（ｄ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０とｊ−ｅ１にめり込んでいる。このパターンでは、辺ｉ−ｅ０上に２つのめり込み点を考える。
【００２３】
粒子受力計算部Ｂ１３０は、各粒子が受ける力を計算するものである。粒子は他の粒子
または壁と接触している場合に、そのめり込み量に応じて接触力を受ける。図３の４つのパターン（ａ）〜（ｄ）に対するそれぞれのめり込み量の定義を図３中の下段の図を用いて説明する。
【００２４】
（ａ）頂点のめり込み：粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０との距離をめり込み量とする。
（ｂ）辺のめり込み（ケース１）：粒子ｉのめり込み頂点辺ｉ−ｅ０が粒子ｊの面ｊ−ｆ１と交差する点（めり込み点）と、（ａ）におけるめり込み頂点ｉ−ｖ０がめり込んでいる面ｊ−ｆ０との距離をめり込み量とする。
（ｃ）辺のめり込み（ケース２）：粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｄ）辺のめり込み（ケース３）：粒子ｉの辺ｉ−ｅ０に対する、粒子ｊの辺ｊ−ｅ０、ｊ−ｅ１との距離をそれぞれの辺とのめり込み量とする。
【００２５】
以上の４パターンにおいて、めり込み量を求めた各粒子表面の２点の方向を接触面法線方向と定める。この接触面法線方向は、粒子ｉの粒子ｊに対するめり込み点（接触点）におけるめり込み方向（又は侵入方向）ということもできる。そして、上記で求めためり込み量に対して、ヘルツの式を用いて接触面法線方向に作用する接触力（バネ力）を求める。また各めり込み点（接触点）に対して、接触が始まってからの接触面（上記接触面法線方向に垂直な面）に平行な方向（接触面内方向）の移動量を求め、ミンドリンの式を用いて粒子ｉに働く接触面内方向に作用する接触力（バネ力）を求める。また接触面法線方向の運動速度（めり込み速度）と面内方向の各運動速度（ずり速度）から、各方向に作用する粘性力（ダンパー力）を求める。そしてそれら接触力と粘性力の和に、粒子に働くその他の力（重力、付着力、電磁力、空気の粘性抵抗等）を加えて、各粒子が受ける並進力と回転力を求める。なお、本実施形態で用いるめり込み量から接触力を求める過程は、上記非特許文献１に記載されている手法と同じであるため、詳しい説明は割愛する。
【００２６】
粒子運動計算部Ｂ１４０は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。すなわち、粒子運動計算部Ｂ１４０は、並進力からニュートンの運動方程式を用いて並進方向の加速度を求め、更に速度と位置を求める。また粒子運動計算部Ｂ１４０は、回転力からオイラーの運動方程式を用いて角加速度を求め、更に回転速度と回転位置を求める。本処理は、非特許文献１で行われている手法と同じであるため、詳しい説明は割愛する。
【００２７】
計算結果出力部Ｂ１５０は、得られた各粒子の位置と速度の結果を出力データ３１として出力するものである。
【００２８】
なお、図２のめり込み基準面設定部Ｂ１３１とめり込み量算出部Ｂ１３２は、図３の（ａ），（ｂ）の場合に対して、粒子受力を算出する元となるめり込み量を求める部分である。これらは本実施形態に固有の部分であり、その内容を以下に詳述する。
【００２９】
まず、図３の（ａ），（ｂ）の場合に対して、粒子ｉの頂点ｉ−ｖ０がめり込んでいる面（図では面ｊ−ｆ０）を、頂点ｉ−ｖ０とめり込み頂点辺ｉ−ｅ０のめり込み量を算出する基準となる面という意味から、”めり込み基準面”と呼ぶことにする。めり込み基準面設定部Ｂ１３１では、粒子ｊの複数の構成面の中からこのめり込み基準面を選択する。本めり込み基準面は、粒子ｉの頂点ｉ−ｖ０に最も近接する粒子ｊの面を選択するのが最も簡易である。しかしそうすると、２粒子が接触した状態で運動したときに、計算が発散することがある。このことを図４を用いて説明する。
【００３０】
図４は粒子ｉの頂点ｉ−ｖ０が粒子ｊにめり込んだ状態で相対運動を行ったときの状態
を模式的に示すものであり、図の表記は図３と同じである。粒子ｉが粒子ｊに対して紙面左方向に移動しており、図４（ａ）は移動前、図４（ｂ）は移動後の状態を示す。移動前においては、粒子ｉの頂点ｉ−ｖ０に最近接な粒子ｊの面は面ｊ−ｆ０であることから、面ｊ−ｆ０がめり込み基準面ということになる。従って、粒子ｉの頂点ｉ−ｖ０とめり込み頂点辺ｉ−ｅ０がめり込んでいることとなり、それぞれのめり込み量は図４（ａ）に示した、δｖ０（ａ）とδｅ０（ａ）ということになる。次に移動後についてみると、粒子ｉの頂点ｉ−ｖ０に最近接な粒子ｊの面は面ｊ−ｆ１であることから、面ｊ−ｆ１がめり込み基準面ということになる。従って、粒子ｉの頂点ｉ−ｖ０とめり込み頂点辺ｉ−ｅ１がめり込んでいることとなり、それぞれのめり込み量は図４（ｂ）に示した、δｖ０（ｂ）とδｅ１（ｂ）ということになる。ここで移動前と後ではめり込み点の位置とめり込み量が大きく違っており、特に移動後のめり込み量δｅ１（ｂ）が極めて大きいため、極めて大きな接触力が発生する。その結果、計算が発散する。
【００３１】
そこでこのような問題を解決するため、本実施形態では粒子ｉの頂点ｉ−ｖ０のめり込み量と、そのめり込み頂点ｉ−ｖ０を端点に持つ粒子ｉのすべての辺（めり込み頂点辺）のめり込み量の総和が最小となる粒子ｊの面をめり込み基準面と定める。すなわち、次式を満たす粒子ｊの面をめり込み基準面と定める。
【数１】
0005901417.tif 000002
ここで、δｖ０は粒子ｉのめり込み頂点ｉ−ｖ０のめり込み量であり、δｅｋはそのめり込み頂点ｉ−ｖ０を端点に持つ粒子ｉのｋ番目のめり込み頂点辺のめり込み量である。なおｎはめり込み頂点辺の数である。
【００３２】
このことを図５を用いて説明する。図５は２つの粒子が接触している状態を示すものであり、図の表記は図３と同じである。図５（ａ）は斜視図、図５（ｂ）、（ｃ）はめり込み基準面を選択する過程を説明する図である。図５（ｂ）、（ｃ）において、２粒子は同じ状態にある。まず粒子ｊの面ｊ−ｆ０をめり込み基準面としたとき、粒子ｉの頂点ｉ−ｖ０とめり込み頂点辺ｉ−ｅ０がめり込んでいることとなり、それぞれのめり込み量は図５（ｂ）におけるδｖ０（ｂ）とδｅ０（ｂ）となる。次に粒子ｊの面ｊ−ｆ１をめり込み基準面としたとき、粒子ｉの頂点ｉ−ｖ０とめり込み頂点辺ｉ−ｅ１がめり込んでいることとなり、それぞれのめり込み量は、図５（ｃ）におけるδｖ０（ｃ）とδｅ１（ｃ）となる。そこで各めり込み基準面におけるめり込み量の総和である（δｖ０（ｂ）＋δｅ０（ｂ））と（δｖ０（ｃ）＋δｅ１（ｃ））を求め、両者を比較する。本図においては、（δｖ０（ｂ）＋δｅ０（ｂ））＜（δｖ０（ｃ）＋δｅ１（ｃ））であるから、粒子ｊの面ｊ−ｆ０をめり込み基準面として採用する。
【００３３】
本方法を採用したときの、図４（ｂ）に代わる図を図６に示す。めり込み基準面としては粒子ｊの面ｊ−ｆ０が採用され、粒子ｉの頂点ｉ−ｖ０とめり込み頂点辺ｉ−ｅ０がめり込んでいることとなる。その結果それぞれのめり込み量はδｖ０とδｅ０となり、図４（ａ）に示した粒子移動前に対してめり込み量は連続的に変化することとなり、安定した挙動解を求めることができる。すなわち、式１を満たす粒子ｊの面をめり込み基準面として採用することにより、各めり込み点におけるめり込み量が時間経過に対して連続的に変化するようになり、各めり込み点に作用する接触力も時間的に連続となる。その結果、接触力の時間的不連続性に起因する一部粒子の異常挙動に端を発した計算の発散等の問題がなくなり、安定して解を得られるようになる。なお、図５と図６に示した例では、めり込んでいる頂点辺は１つしかなかったが、３次元的な形状をした実際のモデルでは複数個存在するケースが発生する。従って式１における括弧中の第２項のように、複数のめり込み
頂点辺に対する考慮が必要となる。
【００３４】
ここで、めり込み頂点辺のもう一方の端点もめり込んでいる場合の処理について説明を追加しておく。図７は図４と同様に粒子ｉの頂点ｉ−ｖ０が粒子ｊにめり込んでいる状態を示すものである。めり込み頂点辺としてｉ−ｅ０とｉ−ｅ１の２つがあるが、そのうちｉ−ｅ１では頂点ｉ−ｖ０と反対側に位置する他方の頂点ｉ−ｖ１が粒子ｊの内部にあり、粒子ｊのどの面とも交差しない。このように他端側の頂点もめり込んでいるめり込み頂点辺に対しては、その他端側頂点のめり込み基準面に対するめり込み量も算出し、めり込み量の総和（式１）に加算するとよい。例えば、図７においてめり込み基準検討面ｊ−ｆ０に対するめり込み量は、図に示すδｖ０、δｅ０、δｅ１の３つということになる。このようにすることで、めり込み頂点辺ｉ−ｅｋの頂点ｉ−ｖ０と反対側の端点である頂点もめり込んでいる場合にも、正しいめり込み基準面が選択されるようになる。その結果、時間経過に対する各めり込み点に作用する接触力の連続性が完全に確実なものとなり、如何なる接触ケースに対しても、安定して解を得られるようになる。
【００３５】
めり込み量算出部Ｂ１３２はめり込み基準面設定部Ｂ１３１で選択しためり込み基準面に対して、めり込み頂点とめり込み頂点辺のめり込み量を算出するものであり、めり込み基準面設定部Ｂ１３１で求めためり込み量δｖ０、δｅ０等を改めて計算するものである。
【００３６】
次に粒子受力計算部Ｂ１３０とめり込み基準面設定部Ｂ１３１の処理フローを図８を用いて説明する。
【００３７】
図８（ａ）は、粒子受力計算部Ｂ１３０の処理フローを示すものである。まず粒子ｉの頂点ｉ−ｖ０を設定し（Ｓ１００）、その頂点が粒子ｊにめり込んでいるかどうか（粒子ｊの内部にあるかどうか）を判定する（Ｓ１１０）。もし粒子ｊにめり込んでいるならば、めり込み基準面設定部Ｂ１３１によってめり込み基準面を選択する（Ｓ１２０）。そして選択しためり込み基準面を元にめり込み量算出部Ｂ１３２によって頂点めり込み量を求め、接触力に換算して粒子ｉに作用する接触力として加算する（Ｓ１３０）。まためり込み頂点辺に対してもめり込み量算出部Ｂ１３２によってめり込み量を求め、接触力に換算して粒子ｉに作用する接触力として加算する（Ｓ１４０）。以上の処理を粒子ｉの全頂点に対して行い、粒子ｉの頂点と頂点辺に働く接触力が求まる。
【００３８】
図８（ｂ）は、めり込み基準面設定部Ｂ１３１の処理フローを示すものであり、図８（ａ）におけるＳ１２０の処理に該当するものである。まずめり込み基準面であるかどうかを検討する粒子ｊの面ｊ−ｆａを設定し（Ｓ２００）、めり込み量の総和Δを非常に大きな値に設定しておく（Ｓ２１０）。次に面ｊ−ｆａに対してめり込み頂点とめり込み頂点辺によるめり込み量の総和Δａの計算を行い（Ｓ２２０）、それがΔよりも小さいならば（Ｓ２３０）、その面をその時点での最有力めり込み基準面ａｍｉｎとして登録する（Ｓ２４０）。上記処理を粒子ｊのすべての面に対して行い、最終的に設定されたａｍｉｎがめり込み基準面である。
【００３９】
＜実施例＞
以上の粒子挙動計算方法を用いて計算を実施した例を、図９と図１０を用いて説明する。図９は解析モデルと解析条件を示すものである。図９（ａ）は解析モデルである。解析モデルは電子写真技術を用いた画像形成装置におけるクリーニングプロセスである。０．１ｍ／ｓｅｃで矢印の方向に移動する感光体ドラムに対してクリーニングブレードを当接させておき、感光体ドラムに乗って運ばれてきたトナーがどのようにクリーニングブレードで堰き止められるかを検討する。クリーニングブレードと感光体ドラムのなす角度βは３０°とした。計算条件は図９（ｂ）の通りである。またトナー粒子の形状として球形と
非球形の２種類を、それぞれ図９（ｃ）に示すような多面体を用いて表現し、両者の運動を比較した。
【００４０】
図１０に計算結果を示す。上段の図はトナー粒子がクリーニングブレードで堰き止められている様子を示し、図中実線は代表的な粒子の軌跡を示す。球形トナーの場合は堰き止められた部分でトナーが循環対流する。一方非球形トナーの場合は、堰き止められた状態で動かないことがわかる。下段の図はそのときの各トナー粒子が受ける圧縮応力の主軸の方向を線で示したものである。下段の図は上段の図のトナーが堰き止められている部分を若干拡大して表示している。球形トナーの場合は循環対流するドラム進行方向先端部分で、トナー粒子は図の上下方向の圧縮応力を受けている。一方非球形トナーの場合は、ドラム面に対して斜め上方向に圧縮応力を受けていることがわかる。従って球形トナーの場合はドラム面に対して垂直方向に、非球形トナーの場合は斜め上方向に、それぞれ堰き止められた部分でクリーニングブレードを押し上げようとする力が働くことになる。この力はトナー粒子がクリーニングブレードに潜り込もうとする力に相当するものであり、トナーのブレードすり抜け力に相当するものである。本結果からトナー粒子がブレードに与える力として、球形トナーの方が非球形トナーよりもすり抜け力として作用し易く、すり抜け易いことがわかる。本結果は実際の現象と符合するものであり、本発明が実際の現象を説明するのに非常に有効であることを示すものである。なお本実施例での計算モデルで、堰き止められたトナーは非常に大きな圧縮応力を受ける。そしてそのことは、粒子挙動計算としてロバスト性を要求される計算であることを意味する。本計算結果は本発明の効果である計算の安定性を土台に得られたものであり、本発明を使用することで、非常に厳しい計算条件下においても、安定に粒子挙動計算ができることを示すものである。なお上記実施例が、本発明による粒子挙動計算方法の、粒子を扱う分野への広汎な適用可能性を示唆していることはいうまでもない。
【符号の説明】
【００４１】
２０：計算装置本体（情報処理装置）、２１：ＣＰＵ、２２：ＲＯＭ、２３：ＲＡＭ
【図１】
0005901417.tif 000003
【図２】
0005901417.tif 000004
【図３】
0005901417.tif 000005
【図４】
0005901417.tif 000006
【図５】
0005901417.tif 000007
【図６】
0005901417.tif 000008
【図７】
0005901417.tif 000009
【図８】
0005901417.tif 000010
【図９】
0005901417.tif 000011
【図１０】
0005901417.tif 000012