粒子挙動解析方法、粒子挙動解析装置、及び解析プログラム
【課題】密集状態にある多数の粒子の挙動を、各粒子の機械特性に加えて、その形状も考慮して計算可能にする技術を提供する。
【解決手段】粒子ｉと物体ｊの形状を多面体により定義する。粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する。算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出し、算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する。
【解決手段】粒子ｉと物体ｊの形状を多面体により定義する。粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する。算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出し、算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、密集状態にある多数の粒子の挙動を計算する方法に関するものであり、更には粒子形状を考慮した粒子挙動計算に関するものである。
【背景技術】
【０００２】
粉体や粒体などの粒子を取り扱う分野では、個々の粒子の挙動を知ることが重要である。従来これらの分野では、実験と観察から粒子の挙動を把握することが多かった。しかし近年、計算機内でその運動を再現させ、その挙動を把握する、いわゆるシミュレーション技術の適用が活発に行われている。正確な粒子挙動のシミュレーションを行うには、粒子のもつ形状、機械的特性等の情報をできるだけ詳細に計算モデルに反映させ、それらの因子を考慮して解を求めることが重要である。また、土木分野における砕石、砂利等の比較的大きな物体も「粒子」と考えることにより、上記シミュレーション技術をこれらの物体の挙動を計算する場合へも適用可能である。
【０００３】
一方、このような粒子を扱う装置の１つに、プリンタ、複写機、ファクシミリ等の電子写真技術を用いた画像形成装置がある。これらの装置では、直径数ミクロンの粉体粒子からなるトナーを用いて、帯電、露光、現像、転写、定着、クリーニングという６つのプロセスを経て、紙上に画像を形成する。そしてこのような画像形成装置の開発においても、シミュレーション技術を用いて個々のトナー粒子の挙動を求め、上記６つのプロセスの設計を効率化する試みが活発に行われている。
【０００４】
粒子が密集した状態での各粒子の挙動をシミュレーションする代表的な方法として個別要素法（ＤＥＭ）がある。個別要素法の技術は非特許文献１に詳細に記載されているが、個々の粒子を球形で表現し、それらの接触に際してはそのめり込み量からヘルツの式とミンドリンの式を用いることで粒子の持つ機械特性を考慮した接触抗力を求める。そしてニュートンの運動方程式（並進方向）とオイラーの運動方程式（回転方向）を用いて、その接触抗力をもとに個々の粒子の時間進展に伴う運動を求めるものである。そして本個別要素法を上記画像形成装置に適用したものとして非特許文献２がある。本文献ではトナー粒子を球として表現し、その挙動をシミュレーションしている。
【０００５】
一方、上記画像形成過程において、トナー粒子の形状が、作成される画像に大きな影響を与えることがわかっている。そこでこのトナー形状がプロセス性能に及ぼす影響をみたものとして非特許文献３がある。複数の球を連結して１つの粒子を表現することによって異形形状（非球形）をしたトナーを表現し（クラスタ粒子）、球形トナーと楕円球トナーによるクリーニングプロセスにおけるトナーすり抜け現象の違いをシミュレーションしている。
【０００６】
このような個別要素法を用いた画像形成装置内での粒子の挙動計算での問題点の１つに、粒子形状の模擬範囲の狭さがある。最適なトナーの形状をシミュレーションで模索するには、様々な粒子形状を扱えるようにする必要がある。そしてそれは球を連結させて表現できるものだけでは不十分である。
【０００７】
一方、異形形状をした物体の衝突前後の運動を計算する方法として非特許文献４がある。これは物体の形状を考慮して接触判定を行い、運動量保存則を用いることで衝突を考慮して、物体の運動を計算するものである。物体の形状を考慮することができるが、衝突に関しては２物体間に限定されており、密集状態にある多数の粒子の運動を計算することはできない。また密集状態にある多数の異形粒子からなる構造物の構造計算を行う手法とし
て非特許文献５がある。しかしこれは、一般的な構造計算が連続体であることに対して、不連続体に対応させたものであり、個々の粒子の運動を計算するものではない。
【０００８】
また、異形形状をした物体の表面を複数の多角形（ポリゴン）で表現している先行例として特許文献１と特許文献２がある。特許文献１は物体や部品の接触や相互干渉を判定する方法として、幾何形状を用いたポリゴンモデルの衝突判定の方法が記載されている。特許文献２は空間あるいは媒質、及び空間あるいは媒質内に散りばめられる物質の投影画像を作成するのに、物質パターンを１つ以上の３角形状のポリゴンにより構成している。しかしこれらはコンピュータグラフィックスやＣＡＤの分野のために開発されたものであり、複数の物体が密集状態にあるときの個々の物体の運動を計算するものではない。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開平１１−３２８４４５号公報
【特許文献２】特開２００１−０４３３９２号公報
【非特許文献】
【００１０】
【非特許文献１】粉体工学会編、粉体シミュレーション入門(1998)
【非特許文献２】仲野、磁性１成分現像システムにおいてトナー電荷がトナーチェーンに及ぼす影響、日本画像学会、Japan Hardcopy 2003 Fall Meeting、p41-44(2003)
【非特許文献３】中山信行、電子写真クリーニングプロセスのシミュレーション、第18回「電磁力関連のダイナミックス」シンポジウム、p165-170(2006)
【非特許文献４】酒井幸市、OpenGLで作る力学アニメーション入門、森北出版(2005)
【非特許文献５】河野昭子他、離散体モデルを用いたバラスト軌道の挙動シミュレーション、http://www.rtri.or.jp/infoce/rrr/2008/10/200810_04.pdf
【発明の概要】
【発明が解決しようとする課題】
【００１１】
本発明の目的は、密集状態にある多数の粒子の挙動を、各粒子の機械特性に加えて、その形状も考慮して計算可能にする技術を提供することにある。
【課題を解決するための手段】
【００１２】
本発明の第一態様は、情報処理装置において、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得ステップと、粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する接触判定ステップと、前記接触判定ステップで算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算ステップと、前記粒子受力計算ステップで算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算ステップと、前記粒子運動計算ステップで算出された粒子の運動を出力する出力ステップと、を実行する粒子挙動解析方法である。
【００１３】
本発明の第二態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得部と、粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出
された各接触点のめり込み方向及びめり込み量を算出する接触判定部と、前記接触判定部で算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算部と、前記粒子受力計算部で算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算部と、前記粒子運動計算部で算出された粒子の運動を出力する計算結果出力部と、を有する粒子挙動解析装置である。
【００１４】
本発明の第三態様は、上述した粒子挙動解析方法の各ステップを情報処理装置に実行させる解析プログラムである。
【発明の効果】
【００１５】
本発明によれば、密集状態にある多数の粒子の挙動を、各粒子の機械特性に加えて、その形状も考慮して計算することができる。
【図面の簡単な説明】
【００１６】
【図１】粒子挙動計算の力学モデルを説明する図。
【図２】解析装置のハードウエア構成を示すブロック図。
【図３】解析装置の機能構成を示すブロック図。
【図４】接触パターンと各パターンにおける粒子めり込み量を示す図。
【図５】粒子の物体への大きなめり込みを示す図。
【図６】解析装置の処理の流れを示すフローチャート。
【図７】実施例としてのトナーの挙動解析の解析モデルと解析条件を示す図。
【図８】実施例と従来例の計算結果の比較を示す図。
【図９】トナー形状がブレードのすり抜けに与える影響を示す図。
【図１０】トナー形状がブレードに及ぼす影響を計算した例を示す図。
【発明を実施するための形態】
【００１７】
以下に本発明の好適な実施形態について説明する。本実施形態は、密集状態にある多数の粒子の挙動をコンピュータ・シミュレーションにより計算するための粒子挙動解析方法及び装置を提供するものである。本実施形態に係る粒子挙動解析方法の特徴の一つは、複数の多角形から構成される多面体（ポリゴンモデル）により各粒子の表面形状を定義した点である。これにより任意の粒子形状を表現できるようになり、従来方法に比べて、多様なケースのシミュレーションが可能になる。本方法のもう一つの特徴は、多面体同士の接触現象を、頂点の面へのめり込みと、辺の辺へのめり込みの２種類で捉えることで、いくつかの接触点による点接触に単純化した点である。具体的には、ある粒子と他の物体（他の粒子を含む）との接触判定において、物体の面にめり込む粒子の頂点、及び、物体の辺にめり込む粒子の辺上の点をそれぞれ接触点として抽出し、その接触点にかかる接触抗力を計算する。これにより、多面体同士の接触抗力の計算を簡易に且つ少ない計算コストで精度よく行うことができる。本方法のその他の特徴は以下の実施形態の説明で述べる。
【００１８】
（粒子挙動計算のモデル）
図１は、本実施形態における粒子挙動計算の力学モデルを説明するものであり、異形粒子（非球形粒子）が他の物体と接触した際に働く接触力を説明する図である。このモデルでは、粒子が他の物体と接触した際に発生するめり込みを、その接触面に対して法線方向と接線方向に分解して考える。粒子は、それぞれのめり込み量に対してスプリングで押し返される。また各スプリングと並列にダンパーを設定する。これはめり込み速度に比例して働く力であり、運動を減衰させるものである。またすべりを考慮するためのスライダを設定している。このように、スプリング、ダンパー、スライダを設定して衝突を考えることは従来の個別要素法と同じである。
【００１９】
（解析装置のハードウエア構成）
図２は、本実施形態における解析装置のハードウエア構成を示すブロック図である。情報処理装置（コンピュータ）２０は、各種判断及び処理を行う中央処理装置（ＣＰＵ）２１と、処理プログラム及び固定データを格納するＲＯＭ２２と、処理データを記憶するＲＡＭ２３と、入出力回路（Ｉ／Ｏ）２４とを備える。また必要に応じて、情報処理装置２０には、Ｉ／Ｏ２４を介して記憶装置、入力装置、表示装置、他のコンピュータなどが接続される。Ｉ／Ｏ２４は、これら外部装置と情報処理装置２０との間でデータをやり取りするための回路である。この情報処理装置２０には、入力データ３０がＩ／Ｏ２４を介して入力される。計算結果は、Ｉ／Ｏ２４を介して出力データ３１として出力される。
【００２０】
（解析装置の機能）
図３は、粒子挙動解析を行う解析プログラムの機能構成を示すブロック図である。この解析プログラムは、ＲＯＭ２２や記憶装置などのコンピュータ読み取り可能な記録媒体に記憶されており、必要に応じて、ＣＰＵ２１によって読み込まれ実行されるソフトウエアである。この解析プログラムは、主な機能（モジュール）として、入力データ読込部（取得部）Ｂ１１０、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、粒子運動計算部Ｂ１４０、計算結果出力部Ｂ１５０を備えている。全体制御部Ｂ１００は、これらのモジュールの処理の順番やデータの引渡しなどの全体制御を担っている。
【００２１】
入力データ読込部Ｂ１１０は、粒子形状データ、壁形状データ、及び各種パラメータを読み込んでＲＡＭ２３に格納する。粒子形状データとは、粒子の表面を複数の多角形（ポリゴン）で表現したものである。ここで多角形としてはその頂点が一平面上にあることが望ましく、３角形を用いるのが最も単純である。また粒子は凸形状に限定すると、他粒子との接触状態の判定が簡便になる。壁形状データも同様であり、壁の表面を複数の多角形で表現したものである。粒子形状データ、壁形状データともに、多角形の頂点座標と、各多角形を構成する頂点の番号とで構成されるデータである。各種パラメータには、粒子と壁の機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）、粒子に働く力（重力、付着力、電磁力、空気の粘性特性等）、及び計算条件としての計算刻み時間（Δｔ）、計算終了時刻がある。壁とは、粒子が移動可能な空間を形成する（規定する）境界面である。
【００２２】
接触判定部Ｂ１２０は、粒子の位置データ（各頂点の座標）に基づいて、粒子が他の物体（周囲の粒子または壁）と接触しているかどうか、またどのように接触しているかを判定する。図４を用いて接触判定の内容を示す。図４は（ａ）〜（ｄ）の４種類の接触パターンを示している。図４の上段は各パターンを３次元的に示したものであり、下段は各パターンにおける２つの粒子の接触によるめり込みの様子を模式的に示すものである。図４において、実線は６面体粒子ｉ、破線は６面体粒子ｊを示す。ｉ−ｖ０は粒子ｉを構成する多角形の頂点、ｉ−ｅ０は粒子ｉを構成する多角形の辺、ｊ−ｅ０とｊ−ｅ１は粒子ｊを構成する多角形の辺、ｊ−ｆ０は粒子ｊを構成する多角形の面、星印は粒子ｉと粒子ｊの接触点を示す。また図４の上段における矢印は粒子ｉが粒子ｊに衝突する方向を示し、また下段における黒丸は粒子ｊの辺を示す。粒子が凸形状である場合、粒子ｉと粒子ｊの接触は、図４に示す次の（ａ）〜（ｄ）の４つのパターンのいずれかとなる。なお、図４では、粒子同士の接触判定を例に挙げているが、粒子と壁の接触判定も全く同じように考えることができる。
【００２３】
（ａ）頂点のめり込み：
粒子ｉを構成する多角形の頂点が粒子ｊの１つの面にめり込む場合である。粒子ｉの頂点が粒子ｊの内部にあれば、接触していると判定する。図４（ａ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込んでいる。このパターンでは、頂点ｉ−ｖ０を、粒子ｉの粒子ｊに対する接触点として考える。
【００２４】
（ｂ）辺のめり込み（ケース１）：
上記（ａ）に加えて、めり込んだ頂点を端点に持つ粒子ｉの辺が粒子ｊの辺にもめり込む場合である。（ａ）のめり込んだ頂点を端点に持つ粒子ｉの辺において、粒子ｊ構成面との交点数が１であり、かつその交点を持つ面が（ａ）におけるめり込み面でなければ、これに該当する。図４（ｂ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込む（上記（ａ）に該当）とともに、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にもめり込んでいる。このパターンでは、頂点ｉ−ｖ０に加え、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点も、粒子ｉの粒子ｊに対する接触点として考える。
【００２５】
（ｃ）辺のめり込み（ケース２）：
粒子ｉの辺が粒子ｊの１つの辺のみにめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、交点を持つ２つの面のそれぞれについて、交点に最近接な辺を抽出し、抽出した辺が同じであれば、これに該当する。図４（ｃ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にめり込んでいる。このパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２６】
（ｄ）辺のめり込み（ケース３）：
粒子ｉの辺が粒子ｊの２つの辺にめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、（ｃ）ではない場合がこれに該当する。図４（ｄ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０とｊ−ｅ１にめり込んでいる。このパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点、及び、辺ｊ−ｅ１に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２７】
粒子受力計算部Ｂ１３０は、各粒子が受ける力を計算するものである。粒子は他の粒子または壁と接触している場合に、そのめり込み量に応じて接触抗力を受ける。図４の４つのパターンに対する各接触点のめり込み量を図４中の下段の図を用いて説明する。
（ａ）頂点のめり込み： 粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０との距離をめり込み量とする。
（ｂ）辺のめり込み（ケース１）： 粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。（頂点ｉ−ｖ０のめり込み量は（ａ）と同じである。）
（ｃ）辺のめり込み（ケース２）： 粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｄ）辺のめり込み（ケース３）： 粒子ｉの辺ｉ−ｅ０に対する、粒子ｊの辺ｊ−ｅ０、ｊ−ｅ１との距離を各接触点のめり込み量とする。
【００２８】
なお図４（ｂ）辺のめり込み（ケース１）と図４（ｄ）辺のめり込み（ケース３）の場合について、接触点とめり込み量を以下の（ｂ´）及び（ｄ´）のように定義してもよい。
（ｂ´）辺のめり込み（ケース１）： めり込み頂点ｉ−ｖ０に加え、この頂点ｉ−ｖ０を端点にもつ辺ｉ−ｅ０が粒子ｊ構成面と交差する交点も接触点とする。そして、粒子ｊ構成面のうち、本接触点とめり込み頂点ｉ−ｖ０に最も近接する面（言い換えれば、本接触点からの距離とめり込み頂点ｉ−ｖ０からの距離の合計が最も小さくなる面）をめり込み基準面とする。図４では、面ｊ−ｆ０がめり込み基準面となる。そして接触点とめり込み基準面との距離をめり込み量とする。
（ｄ´）辺のめり込み（ケース３）： めり込み辺ｉ−ｅ０が粒子ｊ構成面と交差する２つの交点を接触点とする。そして、粒子ｊ構成面のうち、本２つの接触点に最も近接する面（言い換えれば、各接触点からの距離の合計が最も小さくなる面）をめり込み基準面とする。図４では、面ｊ−ｆ０がめり込み基準面となる。そして各接触点とめり込み基準
面との距離を、各接触点のめり込み量とする。
【００２９】
以上の４パターンにおいて、めり込み量（距離）の計算に用いた各粒子表面の２点を結ぶ方向を接触面法線方向と定める。この接触面法線方向は、接触点の（粒子ｊの面又は辺に対する）めり込み方向ということもできる。そして、各接触点に対して、以上で求めためり込み量から、式（１）を用いて粒子ｉに働く接触面法線方向（めり込み方向）の接触抗力を求める。ここで式（１）はヘルツの接触力の式であり、非特許文献１に記載されているものである。式中のｋｎは弾性ばね係数で、粒子のヤング率、ポアソン比、粒子の半径から算出されるものである。なお異形粒子における粒子の半径は、その粒子の外接球の半径で代用する。
【００３０】
また上記の各接触点に対して、接触が始まってからの接触面（上記接触面法線方向に垂直な面）内方向の移動量（ずり量）を求める。そして式（２）を用いて粒子ｉに働く接触面内方向の接触抗力を求める。ここで式（２）はミンドリンの接触力の式であり、非特許文献１に記載されているものである。式中のｋｔは弾性ばね係数で、粒子のヤング率、ポアソン比、粒子の半径から算出されるものである。なお異形粒子における粒子の半径としては、その外接球の半径で代用する。
【００３１】
【数１】
【００３２】
また式（３）、式（４）を用いて粘性力を求める。これは非特許文献１に記載されているもので、めり込み速度（接触面法線方向速度）とずり速度（接触面内方向速度）は、接触点における粒子の速度から算出する。そして式（５）、式（６）のように、各粒子が受ける接触抗力の和に、粒子に働くその他の力（重力、付着力、電磁力、空気の粘性抵抗等）を加えて、各粒子が受ける並進力Ｆと回転力Ｍを求める。
【００３３】
【数２】
【００３４】
【数３】
【００３５】
一般に異形粒子の場合、接触は点、線、面のいずれでも発生し得る。例えば図４の（ｂ）、（ｄ）のパターンでは実際の接触は線接触である。また（ａ）のパターンが粒子ｊの１つの面に３つ以上同時に発生すれば面接触となる。本実施形態では、それらの現象を頂点の面へのめり込み、辺の辺へのめり込みの２種類で捉えることで、全ての接触パターンを１つ又は２つの接触点による点接触に単純化している。これにより、粒子の受ける力を離散化して考慮することができるため、ヘルツの式及びミンドリンの式の適用が可能となり、接触抗力の計算を簡易に且つ少ない計算コストで精度よく行うことができる。また、接触抗力を離散化したことで、後述する粒子の運動計算において、ニュートンの運動方程式及びオイラーの運動方程式が簡易に適用できるようになる。
【００３６】
なお、粒子受力計算部Ｂ１３０のめり込み量の算出において、粒子の複数の頂点が他の１つの物体にめり込む場合には、各頂点（接触点）のめり込み量を小さくする補正を行うことが好ましい。図５を用いてそのことを具体的に説明する。図５では、多面体からなる球形の粒子が物体Ａにめり込んでおり、そのめり込み量が大きいことから粒子の複数の頂点が物体Ａにめり込んでいる。このような場合、粒子は複数の頂点にて矢印で示した方向の接触抗力を受けることとなるため、上記モデルで計算されためり込み量δをそのまま用
いると接触抗力の計算値が過大になる。そこで本実施形態では、式（７）に示すように、補正後のめり込み量δ´の合計値（Σδ´）が補正前のめり込み量δの内の最大値δｍａｘと等しくなるように、各頂点のめり込み量δを補正する。そして、補正後のめり込み量δ´を各頂点の接触抗力の算出に用いる。図５のように複数の頂点がめり込む現象は、粒子表面が微細に多角形分割され、めり込み量が大きい場合などに発生する。本補正によれば、このような場合にも妥当な接触抗力が計算され、粒子挙動計算を正確に行うことができるようになる。なお、式（７）は各頂点のめり込み量δに一定の係数を乗じる補正を行うものであるが、補正式はこれに限らず、めり込み量δや頂点の位置等に応じて補正係数の値を変えてもよいし、非線形な補正を行ってもよい。
【００３７】
【数４】
【００３８】
粒子運動計算部Ｂ１４０は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。式（８）はニュートンの運動方程式である。粒子運動計算部Ｂ１４０は、式（８）から並進加速度を求め、その結果より式（９）、式（１０）から速度と位置を求める。また式（１１）はオイラーの運動方程式である。粒子運動計算部Ｂ１４０は、式（１１）から角加速度を求め、その結果より式（１２）、式（１３）から回転速度と回転位置を求める。なおここで行う処理は、非特許文献１で行われている従来技術と同じである。
【００３９】
【数５】
【００４０】
【数６】
【００４１】
計算結果出力部Ｂ１５０は、得られた各粒子の位置と速度の結果を出力データ３１として出力するものである。
【００４２】
（解析処理のフロー）
次に粒子挙動計算を行う処理の流れについて、図６のフローチャートを用いて説明する。
【００４３】
まず入力データ読込部Ｂ１１０が、入力データを読み込む（Ｓ１００）。これにより、粒子形状データ、壁形状データ、及び各種パラメータが取得される。そして、全体制御部Ｂ１００が、接触の判定を行う粒子ｉと物体ｊをそれぞれ設定する（Ｓ１１０，Ｓ１２０）。物体ｊは、粒子ｉ以外の粒子または壁である。次に接触判定部Ｂ１２０が粒子ｉと物体ｊの接触判定を実施する（Ｓ１３０）。ここで両者が接触していると判定されたならば、粒子受力計算部Ｂ１３０が接触抗力を求め、粒子受力に加算する（Ｓ１４０）。Ｓ１３０とＳ１４０の処理をすべての粒子と物体の組み合わせに対して実施する（Ｓ１５０，Ｓ１６０）。次に求まった各粒子の受力を元に、粒子運動計算部Ｂ１４０が、刻み時間後の粒子の速度と位置を求める（Ｓ１７０）。Ｓ１１０〜Ｓ１７０の処理をユーザが指定した時間まで反復した後（Ｓ１８０，Ｓ１９０）、計算結果出力部Ｂ１５０が、計算結果をファイル出力する（Ｓ２００）。計算結果出力部Ｂ１５０は、計算結果を表示装置に出力したりプリンタに出力したりすることもできる。以上の処理を行うことによって、時間進展に伴う粒子の挙動を計算によって求めることができる。
【００４４】
（実施例）
以上の粒子挙動計算方法を用いて計算を実施した例を、図７〜図１０を用いて説明する。
図７（ａ）は解析モデルを示す。解析モデルは電子写真技術を用いた画像形成装置におけるクリーニングプロセスである。0.1m/secで矢印の方向に移動する感光体ドラムに対してクリーニングブレードをβ＝３０度で当接させておき、感光体ドラム上に乗って運ばれてきたトナーの運動を見る。各トナーの形状を粒子形状データとして与え、感光体ドラム及びクリーニングブレードの形状を壁形状データとして与える。このときのトナーと感光体ドラムの摩擦係数をμＴＰ、トナーとクリーニングブレードの摩擦係数をμＴＢとしている。トナー径は約６μｍに設定した。もしトナーがブレードをすり抜けるとクリーニング不良ということになる。図７（ｂ）は、解析条件として与えた、トナー、感光体ドラム、クリーニングブレードそれぞれの機械特性（ヤング率、ポアソン比、減衰係数、摩擦係
数、比重）を示している。この条件で、刻み時間Δｔを10nsecに設定し、約２百万ステップの計算を実行した。
【００４５】
図８は、（ａ）球形粒子と（ｂ）異形粒子の計算結果を比較したものである。球形粒子の計算は従来の計算方法で実施している。異形粒子は４面体の角をとった形状であり、その表面を８つの３角形と６つの４角形の合計１２の多角形で構成している。なお本計算においては、トナーがブレードに堰き止められた際の挙動をみるため、強制的にトナーがすり抜けないようにしている。球形粒子に対して異形粒子は回転しにくいことから、全体的に動き辛く、トナー粒子の密集度が上がる傾向にある。図８の計算結果では、異形粒子の密集度が高くなっている様子が再現されている。
【００４６】
図９は、感光体ドラムとクリーニングブレードの摩擦係数がトナーのすり抜けに与える影響が、トナーの形状によってどのように変化するかを見たものである。図９の上段に示したように、トナーを、（ａ）正２０面体の球、（ｂ）正１２８０面体の球、（ｃ）扁平球にした場合について計算した。なおここで扁平球とは球を碁石のように押し潰した形状である。計算は感光体ドラム上にトナーを１つ設定して、トナーがすり抜けるかどうかを調べた。図中、○印は、トナーが回転状態で堰き止められる、△はトナーが非回転状態で堰き止められる、×はブレードをすり抜けることを示す。（ａ）の結果から感光体ドラムとブレードの摩擦係数μＴＰ、μＴＢが高くなるとトナーがすり抜けることがわかる。そして（ｂ）の結果から、トナーが完全球に近づくと、感光体とブレードの摩擦係数が小さくてもすり抜けが発生するようになることがわかる。これは球形度が上がることで、トナーが回転し易くなり、それによってすり抜け易くなるものである。一方、（ｃ）に示すように、トナーを扁平にすると、すり抜け難くなることがわかる。これはトナーが回転し難くなることによるものである。（ａ）と（ｃ）を比較すると、すり抜けない部分の分布が異なっていることがわかる。これは（ａ）のすり抜けはトナーが回転しながらすり抜けるのに対して、（ｃ）のすり抜けは回転しないですり抜けることによるものである。以上のように、本方法を使用することで、様々なトナー形状がブレードのすり抜けに与える影響を再現でき、製品設計に極めて有用な知見を得ることが可能となった。
【００４７】
図１０は図７の系におけるβを７６°、粒子数を４５００個としたときのトナーの挙動計算結果である。トナー粒子間の付着力（トナー凝集力）として１ｎＮ、トナーと感光体ドラム、クリーニングブレード間の付着力として１ｎＮを与えている。図７（ａ）の紙面奥行き方向の解析領域をトナー３個分にとり、その境界条件は周期境界とした。図１０（ａ）は球形トナー、図１０（ｂ）は異形トナーの結果を示す。それぞれ左図はブレードの感光体ドラム当接部上流でトナーが堰き止められてクリーニングされる様子、右図はブレードが受ける力（全粒子から受ける力の積分値）の時間変化を示す。異形形状としては図８と同じく４面体の角をとった形状とした。また右図中、Ｆｘ、Ｆｙはそれぞれ左図の水平右方向、垂直上方向の力の成分を示し、｜Ｆｘｙ｜は絶対値である。左図から球形トナーの場合は堰き止められたトナーが上方に堆積していくこと、異形トナーの場合は塊状となって溜まることがわかる。この違いはトナーの回転のし易さから発現するものであり、実際の現象を観察しても同様な違いが見られる。また右図から球形粒子の方がブレード受力は若干大きく、特に垂直方向の力Ｆｙにその差が顕著であることがわかる。これもトナー粒子の回転のし易さに起因するものであり、回転し易いと粒子にブレードと感光体ドラムの間に潜り込もうとする力が働くことによるものである。また異形トナーの場合、ブレード受力に時間的な脈動が見られる。これは溜まったトナーの塊形状が時間をおいて変化することによるもので、本塊形状の周期的な変動は実際の系でも観察されるものである。このように本方法を用いることで、トナーの形状因子がクリーニング特性に与える影響のメカニズムを解析することが可能となり、トナー形状、ブレード材質、構成等の設計に大きく寄与できる。
【００４８】
以上、本発明の実施例を電子写真技術を用いた画像形成装置におけるクリーニングプロセスに適用した例について説明した。ここでの説明から明らかなように本発明による粒子挙動計算方法は、粒子を扱う分野に幅広く適用できるものである。更には、土木分野における砕石、砂利等の比較的大きな物体も「粒子」と考えることにより、本発明をこれらの物体の挙動の計算へも適用可能である。
【００４９】
なお本発明ではめり込み量に対する接触抗力の計算に、異形形状粒子にも拘わらず、球形粒子時に成立するヘルツの式とミンドリンの式を使用している。従って個々の接触点での接触抗力値の誤差は大きいと思われる。しかしその誤差は粒子の異形性が挙動に与える影響に比べて十分に小さく、粒子の挙動傾向を見るには本方法で十分であることは、これらの計算結果が、観察結果とよく一致することからわかっている。
【符号の説明】
【００５０】
２０：情報処理装置、Ｂ１１０：入力データ読込部、Ｂ１２０：接触判定部、Ｂ１３０：粒子受力計算部、Ｂ１４０：粒子運動計算部、Ｂ１５０：計算結果出力部
【技術分野】
【０００１】
本発明は、密集状態にある多数の粒子の挙動を計算する方法に関するものであり、更には粒子形状を考慮した粒子挙動計算に関するものである。
【背景技術】
【０００２】
粉体や粒体などの粒子を取り扱う分野では、個々の粒子の挙動を知ることが重要である。従来これらの分野では、実験と観察から粒子の挙動を把握することが多かった。しかし近年、計算機内でその運動を再現させ、その挙動を把握する、いわゆるシミュレーション技術の適用が活発に行われている。正確な粒子挙動のシミュレーションを行うには、粒子のもつ形状、機械的特性等の情報をできるだけ詳細に計算モデルに反映させ、それらの因子を考慮して解を求めることが重要である。また、土木分野における砕石、砂利等の比較的大きな物体も「粒子」と考えることにより、上記シミュレーション技術をこれらの物体の挙動を計算する場合へも適用可能である。
【０００３】
一方、このような粒子を扱う装置の１つに、プリンタ、複写機、ファクシミリ等の電子写真技術を用いた画像形成装置がある。これらの装置では、直径数ミクロンの粉体粒子からなるトナーを用いて、帯電、露光、現像、転写、定着、クリーニングという６つのプロセスを経て、紙上に画像を形成する。そしてこのような画像形成装置の開発においても、シミュレーション技術を用いて個々のトナー粒子の挙動を求め、上記６つのプロセスの設計を効率化する試みが活発に行われている。
【０００４】
粒子が密集した状態での各粒子の挙動をシミュレーションする代表的な方法として個別要素法（ＤＥＭ）がある。個別要素法の技術は非特許文献１に詳細に記載されているが、個々の粒子を球形で表現し、それらの接触に際してはそのめり込み量からヘルツの式とミンドリンの式を用いることで粒子の持つ機械特性を考慮した接触抗力を求める。そしてニュートンの運動方程式（並進方向）とオイラーの運動方程式（回転方向）を用いて、その接触抗力をもとに個々の粒子の時間進展に伴う運動を求めるものである。そして本個別要素法を上記画像形成装置に適用したものとして非特許文献２がある。本文献ではトナー粒子を球として表現し、その挙動をシミュレーションしている。
【０００５】
一方、上記画像形成過程において、トナー粒子の形状が、作成される画像に大きな影響を与えることがわかっている。そこでこのトナー形状がプロセス性能に及ぼす影響をみたものとして非特許文献３がある。複数の球を連結して１つの粒子を表現することによって異形形状（非球形）をしたトナーを表現し（クラスタ粒子）、球形トナーと楕円球トナーによるクリーニングプロセスにおけるトナーすり抜け現象の違いをシミュレーションしている。
【０００６】
このような個別要素法を用いた画像形成装置内での粒子の挙動計算での問題点の１つに、粒子形状の模擬範囲の狭さがある。最適なトナーの形状をシミュレーションで模索するには、様々な粒子形状を扱えるようにする必要がある。そしてそれは球を連結させて表現できるものだけでは不十分である。
【０００７】
一方、異形形状をした物体の衝突前後の運動を計算する方法として非特許文献４がある。これは物体の形状を考慮して接触判定を行い、運動量保存則を用いることで衝突を考慮して、物体の運動を計算するものである。物体の形状を考慮することができるが、衝突に関しては２物体間に限定されており、密集状態にある多数の粒子の運動を計算することはできない。また密集状態にある多数の異形粒子からなる構造物の構造計算を行う手法とし
て非特許文献５がある。しかしこれは、一般的な構造計算が連続体であることに対して、不連続体に対応させたものであり、個々の粒子の運動を計算するものではない。
【０００８】
また、異形形状をした物体の表面を複数の多角形（ポリゴン）で表現している先行例として特許文献１と特許文献２がある。特許文献１は物体や部品の接触や相互干渉を判定する方法として、幾何形状を用いたポリゴンモデルの衝突判定の方法が記載されている。特許文献２は空間あるいは媒質、及び空間あるいは媒質内に散りばめられる物質の投影画像を作成するのに、物質パターンを１つ以上の３角形状のポリゴンにより構成している。しかしこれらはコンピュータグラフィックスやＣＡＤの分野のために開発されたものであり、複数の物体が密集状態にあるときの個々の物体の運動を計算するものではない。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開平１１−３２８４４５号公報
【特許文献２】特開２００１−０４３３９２号公報
【非特許文献】
【００１０】
【非特許文献１】粉体工学会編、粉体シミュレーション入門(1998)
【非特許文献２】仲野、磁性１成分現像システムにおいてトナー電荷がトナーチェーンに及ぼす影響、日本画像学会、Japan Hardcopy 2003 Fall Meeting、p41-44(2003)
【非特許文献３】中山信行、電子写真クリーニングプロセスのシミュレーション、第18回「電磁力関連のダイナミックス」シンポジウム、p165-170(2006)
【非特許文献４】酒井幸市、OpenGLで作る力学アニメーション入門、森北出版(2005)
【非特許文献５】河野昭子他、離散体モデルを用いたバラスト軌道の挙動シミュレーション、http://www.rtri.or.jp/infoce/rrr/2008/10/200810_04.pdf
【発明の概要】
【発明が解決しようとする課題】
【００１１】
本発明の目的は、密集状態にある多数の粒子の挙動を、各粒子の機械特性に加えて、その形状も考慮して計算可能にする技術を提供することにある。
【課題を解決するための手段】
【００１２】
本発明の第一態様は、情報処理装置において、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得ステップと、粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する接触判定ステップと、前記接触判定ステップで算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算ステップと、前記粒子受力計算ステップで算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算ステップと、前記粒子運動計算ステップで算出された粒子の運動を出力する出力ステップと、を実行する粒子挙動解析方法である。
【００１３】
本発明の第二態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得部と、粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出
された各接触点のめり込み方向及びめり込み量を算出する接触判定部と、前記接触判定部で算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算部と、前記粒子受力計算部で算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算部と、前記粒子運動計算部で算出された粒子の運動を出力する計算結果出力部と、を有する粒子挙動解析装置である。
【００１４】
本発明の第三態様は、上述した粒子挙動解析方法の各ステップを情報処理装置に実行させる解析プログラムである。
【発明の効果】
【００１５】
本発明によれば、密集状態にある多数の粒子の挙動を、各粒子の機械特性に加えて、その形状も考慮して計算することができる。
【図面の簡単な説明】
【００１６】
【図１】粒子挙動計算の力学モデルを説明する図。
【図２】解析装置のハードウエア構成を示すブロック図。
【図３】解析装置の機能構成を示すブロック図。
【図４】接触パターンと各パターンにおける粒子めり込み量を示す図。
【図５】粒子の物体への大きなめり込みを示す図。
【図６】解析装置の処理の流れを示すフローチャート。
【図７】実施例としてのトナーの挙動解析の解析モデルと解析条件を示す図。
【図８】実施例と従来例の計算結果の比較を示す図。
【図９】トナー形状がブレードのすり抜けに与える影響を示す図。
【図１０】トナー形状がブレードに及ぼす影響を計算した例を示す図。
【発明を実施するための形態】
【００１７】
以下に本発明の好適な実施形態について説明する。本実施形態は、密集状態にある多数の粒子の挙動をコンピュータ・シミュレーションにより計算するための粒子挙動解析方法及び装置を提供するものである。本実施形態に係る粒子挙動解析方法の特徴の一つは、複数の多角形から構成される多面体（ポリゴンモデル）により各粒子の表面形状を定義した点である。これにより任意の粒子形状を表現できるようになり、従来方法に比べて、多様なケースのシミュレーションが可能になる。本方法のもう一つの特徴は、多面体同士の接触現象を、頂点の面へのめり込みと、辺の辺へのめり込みの２種類で捉えることで、いくつかの接触点による点接触に単純化した点である。具体的には、ある粒子と他の物体（他の粒子を含む）との接触判定において、物体の面にめり込む粒子の頂点、及び、物体の辺にめり込む粒子の辺上の点をそれぞれ接触点として抽出し、その接触点にかかる接触抗力を計算する。これにより、多面体同士の接触抗力の計算を簡易に且つ少ない計算コストで精度よく行うことができる。本方法のその他の特徴は以下の実施形態の説明で述べる。
【００１８】
（粒子挙動計算のモデル）
図１は、本実施形態における粒子挙動計算の力学モデルを説明するものであり、異形粒子（非球形粒子）が他の物体と接触した際に働く接触力を説明する図である。このモデルでは、粒子が他の物体と接触した際に発生するめり込みを、その接触面に対して法線方向と接線方向に分解して考える。粒子は、それぞれのめり込み量に対してスプリングで押し返される。また各スプリングと並列にダンパーを設定する。これはめり込み速度に比例して働く力であり、運動を減衰させるものである。またすべりを考慮するためのスライダを設定している。このように、スプリング、ダンパー、スライダを設定して衝突を考えることは従来の個別要素法と同じである。
【００１９】
（解析装置のハードウエア構成）
図２は、本実施形態における解析装置のハードウエア構成を示すブロック図である。情報処理装置（コンピュータ）２０は、各種判断及び処理を行う中央処理装置（ＣＰＵ）２１と、処理プログラム及び固定データを格納するＲＯＭ２２と、処理データを記憶するＲＡＭ２３と、入出力回路（Ｉ／Ｏ）２４とを備える。また必要に応じて、情報処理装置２０には、Ｉ／Ｏ２４を介して記憶装置、入力装置、表示装置、他のコンピュータなどが接続される。Ｉ／Ｏ２４は、これら外部装置と情報処理装置２０との間でデータをやり取りするための回路である。この情報処理装置２０には、入力データ３０がＩ／Ｏ２４を介して入力される。計算結果は、Ｉ／Ｏ２４を介して出力データ３１として出力される。
【００２０】
（解析装置の機能）
図３は、粒子挙動解析を行う解析プログラムの機能構成を示すブロック図である。この解析プログラムは、ＲＯＭ２２や記憶装置などのコンピュータ読み取り可能な記録媒体に記憶されており、必要に応じて、ＣＰＵ２１によって読み込まれ実行されるソフトウエアである。この解析プログラムは、主な機能（モジュール）として、入力データ読込部（取得部）Ｂ１１０、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、粒子運動計算部Ｂ１４０、計算結果出力部Ｂ１５０を備えている。全体制御部Ｂ１００は、これらのモジュールの処理の順番やデータの引渡しなどの全体制御を担っている。
【００２１】
入力データ読込部Ｂ１１０は、粒子形状データ、壁形状データ、及び各種パラメータを読み込んでＲＡＭ２３に格納する。粒子形状データとは、粒子の表面を複数の多角形（ポリゴン）で表現したものである。ここで多角形としてはその頂点が一平面上にあることが望ましく、３角形を用いるのが最も単純である。また粒子は凸形状に限定すると、他粒子との接触状態の判定が簡便になる。壁形状データも同様であり、壁の表面を複数の多角形で表現したものである。粒子形状データ、壁形状データともに、多角形の頂点座標と、各多角形を構成する頂点の番号とで構成されるデータである。各種パラメータには、粒子と壁の機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）、粒子に働く力（重力、付着力、電磁力、空気の粘性特性等）、及び計算条件としての計算刻み時間（Δｔ）、計算終了時刻がある。壁とは、粒子が移動可能な空間を形成する（規定する）境界面である。
【００２２】
接触判定部Ｂ１２０は、粒子の位置データ（各頂点の座標）に基づいて、粒子が他の物体（周囲の粒子または壁）と接触しているかどうか、またどのように接触しているかを判定する。図４を用いて接触判定の内容を示す。図４は（ａ）〜（ｄ）の４種類の接触パターンを示している。図４の上段は各パターンを３次元的に示したものであり、下段は各パターンにおける２つの粒子の接触によるめり込みの様子を模式的に示すものである。図４において、実線は６面体粒子ｉ、破線は６面体粒子ｊを示す。ｉ−ｖ０は粒子ｉを構成する多角形の頂点、ｉ−ｅ０は粒子ｉを構成する多角形の辺、ｊ−ｅ０とｊ−ｅ１は粒子ｊを構成する多角形の辺、ｊ−ｆ０は粒子ｊを構成する多角形の面、星印は粒子ｉと粒子ｊの接触点を示す。また図４の上段における矢印は粒子ｉが粒子ｊに衝突する方向を示し、また下段における黒丸は粒子ｊの辺を示す。粒子が凸形状である場合、粒子ｉと粒子ｊの接触は、図４に示す次の（ａ）〜（ｄ）の４つのパターンのいずれかとなる。なお、図４では、粒子同士の接触判定を例に挙げているが、粒子と壁の接触判定も全く同じように考えることができる。
【００２３】
（ａ）頂点のめり込み：
粒子ｉを構成する多角形の頂点が粒子ｊの１つの面にめり込む場合である。粒子ｉの頂点が粒子ｊの内部にあれば、接触していると判定する。図４（ａ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込んでいる。このパターンでは、頂点ｉ−ｖ０を、粒子ｉの粒子ｊに対する接触点として考える。
【００２４】
（ｂ）辺のめり込み（ケース１）：
上記（ａ）に加えて、めり込んだ頂点を端点に持つ粒子ｉの辺が粒子ｊの辺にもめり込む場合である。（ａ）のめり込んだ頂点を端点に持つ粒子ｉの辺において、粒子ｊ構成面との交点数が１であり、かつその交点を持つ面が（ａ）におけるめり込み面でなければ、これに該当する。図４（ｂ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込む（上記（ａ）に該当）とともに、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にもめり込んでいる。このパターンでは、頂点ｉ−ｖ０に加え、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点も、粒子ｉの粒子ｊに対する接触点として考える。
【００２５】
（ｃ）辺のめり込み（ケース２）：
粒子ｉの辺が粒子ｊの１つの辺のみにめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、交点を持つ２つの面のそれぞれについて、交点に最近接な辺を抽出し、抽出した辺が同じであれば、これに該当する。図４（ｃ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にめり込んでいる。このパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２６】
（ｄ）辺のめり込み（ケース３）：
粒子ｉの辺が粒子ｊの２つの辺にめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、（ｃ）ではない場合がこれに該当する。図４（ｄ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０とｊ−ｅ１にめり込んでいる。このパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点、及び、辺ｊ−ｅ１に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
【００２７】
粒子受力計算部Ｂ１３０は、各粒子が受ける力を計算するものである。粒子は他の粒子または壁と接触している場合に、そのめり込み量に応じて接触抗力を受ける。図４の４つのパターンに対する各接触点のめり込み量を図４中の下段の図を用いて説明する。
（ａ）頂点のめり込み： 粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０との距離をめり込み量とする。
（ｂ）辺のめり込み（ケース１）： 粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。（頂点ｉ−ｖ０のめり込み量は（ａ）と同じである。）
（ｃ）辺のめり込み（ケース２）： 粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｄ）辺のめり込み（ケース３）： 粒子ｉの辺ｉ−ｅ０に対する、粒子ｊの辺ｊ−ｅ０、ｊ−ｅ１との距離を各接触点のめり込み量とする。
【００２８】
なお図４（ｂ）辺のめり込み（ケース１）と図４（ｄ）辺のめり込み（ケース３）の場合について、接触点とめり込み量を以下の（ｂ´）及び（ｄ´）のように定義してもよい。
（ｂ´）辺のめり込み（ケース１）： めり込み頂点ｉ−ｖ０に加え、この頂点ｉ−ｖ０を端点にもつ辺ｉ−ｅ０が粒子ｊ構成面と交差する交点も接触点とする。そして、粒子ｊ構成面のうち、本接触点とめり込み頂点ｉ−ｖ０に最も近接する面（言い換えれば、本接触点からの距離とめり込み頂点ｉ−ｖ０からの距離の合計が最も小さくなる面）をめり込み基準面とする。図４では、面ｊ−ｆ０がめり込み基準面となる。そして接触点とめり込み基準面との距離をめり込み量とする。
（ｄ´）辺のめり込み（ケース３）： めり込み辺ｉ−ｅ０が粒子ｊ構成面と交差する２つの交点を接触点とする。そして、粒子ｊ構成面のうち、本２つの接触点に最も近接する面（言い換えれば、各接触点からの距離の合計が最も小さくなる面）をめり込み基準面とする。図４では、面ｊ−ｆ０がめり込み基準面となる。そして各接触点とめり込み基準
面との距離を、各接触点のめり込み量とする。
【００２９】
以上の４パターンにおいて、めり込み量（距離）の計算に用いた各粒子表面の２点を結ぶ方向を接触面法線方向と定める。この接触面法線方向は、接触点の（粒子ｊの面又は辺に対する）めり込み方向ということもできる。そして、各接触点に対して、以上で求めためり込み量から、式（１）を用いて粒子ｉに働く接触面法線方向（めり込み方向）の接触抗力を求める。ここで式（１）はヘルツの接触力の式であり、非特許文献１に記載されているものである。式中のｋｎは弾性ばね係数で、粒子のヤング率、ポアソン比、粒子の半径から算出されるものである。なお異形粒子における粒子の半径は、その粒子の外接球の半径で代用する。
【００３０】
また上記の各接触点に対して、接触が始まってからの接触面（上記接触面法線方向に垂直な面）内方向の移動量（ずり量）を求める。そして式（２）を用いて粒子ｉに働く接触面内方向の接触抗力を求める。ここで式（２）はミンドリンの接触力の式であり、非特許文献１に記載されているものである。式中のｋｔは弾性ばね係数で、粒子のヤング率、ポアソン比、粒子の半径から算出されるものである。なお異形粒子における粒子の半径としては、その外接球の半径で代用する。
【００３１】
【数１】
【００３２】
また式（３）、式（４）を用いて粘性力を求める。これは非特許文献１に記載されているもので、めり込み速度（接触面法線方向速度）とずり速度（接触面内方向速度）は、接触点における粒子の速度から算出する。そして式（５）、式（６）のように、各粒子が受ける接触抗力の和に、粒子に働くその他の力（重力、付着力、電磁力、空気の粘性抵抗等）を加えて、各粒子が受ける並進力Ｆと回転力Ｍを求める。
【００３３】
【数２】
【００３４】
【数３】
【００３５】
一般に異形粒子の場合、接触は点、線、面のいずれでも発生し得る。例えば図４の（ｂ）、（ｄ）のパターンでは実際の接触は線接触である。また（ａ）のパターンが粒子ｊの１つの面に３つ以上同時に発生すれば面接触となる。本実施形態では、それらの現象を頂点の面へのめり込み、辺の辺へのめり込みの２種類で捉えることで、全ての接触パターンを１つ又は２つの接触点による点接触に単純化している。これにより、粒子の受ける力を離散化して考慮することができるため、ヘルツの式及びミンドリンの式の適用が可能となり、接触抗力の計算を簡易に且つ少ない計算コストで精度よく行うことができる。また、接触抗力を離散化したことで、後述する粒子の運動計算において、ニュートンの運動方程式及びオイラーの運動方程式が簡易に適用できるようになる。
【００３６】
なお、粒子受力計算部Ｂ１３０のめり込み量の算出において、粒子の複数の頂点が他の１つの物体にめり込む場合には、各頂点（接触点）のめり込み量を小さくする補正を行うことが好ましい。図５を用いてそのことを具体的に説明する。図５では、多面体からなる球形の粒子が物体Ａにめり込んでおり、そのめり込み量が大きいことから粒子の複数の頂点が物体Ａにめり込んでいる。このような場合、粒子は複数の頂点にて矢印で示した方向の接触抗力を受けることとなるため、上記モデルで計算されためり込み量δをそのまま用
いると接触抗力の計算値が過大になる。そこで本実施形態では、式（７）に示すように、補正後のめり込み量δ´の合計値（Σδ´）が補正前のめり込み量δの内の最大値δｍａｘと等しくなるように、各頂点のめり込み量δを補正する。そして、補正後のめり込み量δ´を各頂点の接触抗力の算出に用いる。図５のように複数の頂点がめり込む現象は、粒子表面が微細に多角形分割され、めり込み量が大きい場合などに発生する。本補正によれば、このような場合にも妥当な接触抗力が計算され、粒子挙動計算を正確に行うことができるようになる。なお、式（７）は各頂点のめり込み量δに一定の係数を乗じる補正を行うものであるが、補正式はこれに限らず、めり込み量δや頂点の位置等に応じて補正係数の値を変えてもよいし、非線形な補正を行ってもよい。
【００３７】
【数４】
【００３８】
粒子運動計算部Ｂ１４０は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。式（８）はニュートンの運動方程式である。粒子運動計算部Ｂ１４０は、式（８）から並進加速度を求め、その結果より式（９）、式（１０）から速度と位置を求める。また式（１１）はオイラーの運動方程式である。粒子運動計算部Ｂ１４０は、式（１１）から角加速度を求め、その結果より式（１２）、式（１３）から回転速度と回転位置を求める。なおここで行う処理は、非特許文献１で行われている従来技術と同じである。
【００３９】
【数５】
【００４０】
【数６】
【００４１】
計算結果出力部Ｂ１５０は、得られた各粒子の位置と速度の結果を出力データ３１として出力するものである。
【００４２】
（解析処理のフロー）
次に粒子挙動計算を行う処理の流れについて、図６のフローチャートを用いて説明する。
【００４３】
まず入力データ読込部Ｂ１１０が、入力データを読み込む（Ｓ１００）。これにより、粒子形状データ、壁形状データ、及び各種パラメータが取得される。そして、全体制御部Ｂ１００が、接触の判定を行う粒子ｉと物体ｊをそれぞれ設定する（Ｓ１１０，Ｓ１２０）。物体ｊは、粒子ｉ以外の粒子または壁である。次に接触判定部Ｂ１２０が粒子ｉと物体ｊの接触判定を実施する（Ｓ１３０）。ここで両者が接触していると判定されたならば、粒子受力計算部Ｂ１３０が接触抗力を求め、粒子受力に加算する（Ｓ１４０）。Ｓ１３０とＳ１４０の処理をすべての粒子と物体の組み合わせに対して実施する（Ｓ１５０，Ｓ１６０）。次に求まった各粒子の受力を元に、粒子運動計算部Ｂ１４０が、刻み時間後の粒子の速度と位置を求める（Ｓ１７０）。Ｓ１１０〜Ｓ１７０の処理をユーザが指定した時間まで反復した後（Ｓ１８０，Ｓ１９０）、計算結果出力部Ｂ１５０が、計算結果をファイル出力する（Ｓ２００）。計算結果出力部Ｂ１５０は、計算結果を表示装置に出力したりプリンタに出力したりすることもできる。以上の処理を行うことによって、時間進展に伴う粒子の挙動を計算によって求めることができる。
【００４４】
（実施例）
以上の粒子挙動計算方法を用いて計算を実施した例を、図７〜図１０を用いて説明する。
図７（ａ）は解析モデルを示す。解析モデルは電子写真技術を用いた画像形成装置におけるクリーニングプロセスである。0.1m/secで矢印の方向に移動する感光体ドラムに対してクリーニングブレードをβ＝３０度で当接させておき、感光体ドラム上に乗って運ばれてきたトナーの運動を見る。各トナーの形状を粒子形状データとして与え、感光体ドラム及びクリーニングブレードの形状を壁形状データとして与える。このときのトナーと感光体ドラムの摩擦係数をμＴＰ、トナーとクリーニングブレードの摩擦係数をμＴＢとしている。トナー径は約６μｍに設定した。もしトナーがブレードをすり抜けるとクリーニング不良ということになる。図７（ｂ）は、解析条件として与えた、トナー、感光体ドラム、クリーニングブレードそれぞれの機械特性（ヤング率、ポアソン比、減衰係数、摩擦係
数、比重）を示している。この条件で、刻み時間Δｔを10nsecに設定し、約２百万ステップの計算を実行した。
【００４５】
図８は、（ａ）球形粒子と（ｂ）異形粒子の計算結果を比較したものである。球形粒子の計算は従来の計算方法で実施している。異形粒子は４面体の角をとった形状であり、その表面を８つの３角形と６つの４角形の合計１２の多角形で構成している。なお本計算においては、トナーがブレードに堰き止められた際の挙動をみるため、強制的にトナーがすり抜けないようにしている。球形粒子に対して異形粒子は回転しにくいことから、全体的に動き辛く、トナー粒子の密集度が上がる傾向にある。図８の計算結果では、異形粒子の密集度が高くなっている様子が再現されている。
【００４６】
図９は、感光体ドラムとクリーニングブレードの摩擦係数がトナーのすり抜けに与える影響が、トナーの形状によってどのように変化するかを見たものである。図９の上段に示したように、トナーを、（ａ）正２０面体の球、（ｂ）正１２８０面体の球、（ｃ）扁平球にした場合について計算した。なおここで扁平球とは球を碁石のように押し潰した形状である。計算は感光体ドラム上にトナーを１つ設定して、トナーがすり抜けるかどうかを調べた。図中、○印は、トナーが回転状態で堰き止められる、△はトナーが非回転状態で堰き止められる、×はブレードをすり抜けることを示す。（ａ）の結果から感光体ドラムとブレードの摩擦係数μＴＰ、μＴＢが高くなるとトナーがすり抜けることがわかる。そして（ｂ）の結果から、トナーが完全球に近づくと、感光体とブレードの摩擦係数が小さくてもすり抜けが発生するようになることがわかる。これは球形度が上がることで、トナーが回転し易くなり、それによってすり抜け易くなるものである。一方、（ｃ）に示すように、トナーを扁平にすると、すり抜け難くなることがわかる。これはトナーが回転し難くなることによるものである。（ａ）と（ｃ）を比較すると、すり抜けない部分の分布が異なっていることがわかる。これは（ａ）のすり抜けはトナーが回転しながらすり抜けるのに対して、（ｃ）のすり抜けは回転しないですり抜けることによるものである。以上のように、本方法を使用することで、様々なトナー形状がブレードのすり抜けに与える影響を再現でき、製品設計に極めて有用な知見を得ることが可能となった。
【００４７】
図１０は図７の系におけるβを７６°、粒子数を４５００個としたときのトナーの挙動計算結果である。トナー粒子間の付着力（トナー凝集力）として１ｎＮ、トナーと感光体ドラム、クリーニングブレード間の付着力として１ｎＮを与えている。図７（ａ）の紙面奥行き方向の解析領域をトナー３個分にとり、その境界条件は周期境界とした。図１０（ａ）は球形トナー、図１０（ｂ）は異形トナーの結果を示す。それぞれ左図はブレードの感光体ドラム当接部上流でトナーが堰き止められてクリーニングされる様子、右図はブレードが受ける力（全粒子から受ける力の積分値）の時間変化を示す。異形形状としては図８と同じく４面体の角をとった形状とした。また右図中、Ｆｘ、Ｆｙはそれぞれ左図の水平右方向、垂直上方向の力の成分を示し、｜Ｆｘｙ｜は絶対値である。左図から球形トナーの場合は堰き止められたトナーが上方に堆積していくこと、異形トナーの場合は塊状となって溜まることがわかる。この違いはトナーの回転のし易さから発現するものであり、実際の現象を観察しても同様な違いが見られる。また右図から球形粒子の方がブレード受力は若干大きく、特に垂直方向の力Ｆｙにその差が顕著であることがわかる。これもトナー粒子の回転のし易さに起因するものであり、回転し易いと粒子にブレードと感光体ドラムの間に潜り込もうとする力が働くことによるものである。また異形トナーの場合、ブレード受力に時間的な脈動が見られる。これは溜まったトナーの塊形状が時間をおいて変化することによるもので、本塊形状の周期的な変動は実際の系でも観察されるものである。このように本方法を用いることで、トナーの形状因子がクリーニング特性に与える影響のメカニズムを解析することが可能となり、トナー形状、ブレード材質、構成等の設計に大きく寄与できる。
【００４８】
以上、本発明の実施例を電子写真技術を用いた画像形成装置におけるクリーニングプロセスに適用した例について説明した。ここでの説明から明らかなように本発明による粒子挙動計算方法は、粒子を扱う分野に幅広く適用できるものである。更には、土木分野における砕石、砂利等の比較的大きな物体も「粒子」と考えることにより、本発明をこれらの物体の挙動の計算へも適用可能である。
【００４９】
なお本発明ではめり込み量に対する接触抗力の計算に、異形形状粒子にも拘わらず、球形粒子時に成立するヘルツの式とミンドリンの式を使用している。従って個々の接触点での接触抗力値の誤差は大きいと思われる。しかしその誤差は粒子の異形性が挙動に与える影響に比べて十分に小さく、粒子の挙動傾向を見るには本方法で十分であることは、これらの計算結果が、観察結果とよく一致することからわかっている。
【符号の説明】
【００５０】
２０：情報処理装置、Ｂ１１０：入力データ読込部、Ｂ１２０：接触判定部、Ｂ１３０：粒子受力計算部、Ｂ１４０：粒子運動計算部、Ｂ１５０：計算結果出力部
【特許請求の範囲】
【請求項１】
情報処理装置において、密集状態にある複数の粒子の運動を計算する方法であって、
情報処理装置が、
粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得ステップと、
粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する接触判定ステップと、
前記接触判定ステップで算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算ステップと、
前記粒子受力計算ステップで算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算ステップと、
前記粒子運動計算ステップで算出された粒子の運動を出力する出力ステップと、
を実行することを特徴とする粒子挙動解析方法。
【請求項２】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、前記辺ｊ−ｅ１に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項１に記載の粒子挙動解析方法。
【請求項３】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ１との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項２に記載の粒子挙動解析方法。
【請求項４】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０が物体ｊの面と交差する点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０が物体ｊの２つの面それぞれと交差する２つの点、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項１に記載の粒子挙動解析方法。
【請求項５】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項４に記載の粒子挙動解析方法。
【請求項６】
前記粒子受力計算ステップにおいて、
前記接触点のめり込み量から、ヘルツの式を用いて、前記めり込み方向の接触抗力を算出することを特徴とする請求項１〜５のうちいずれか１項に記載の粒子挙動解析方法。
【請求項７】
前記粒子受力計算ステップにおいて、
前記めり込み方向に垂直な面内での前記接触点の移動量を算出し、算出された移動量から、ミンドリンの式を用いて、前記めり込み方向に垂直な方向の接触抗力を算出することを特徴とする請求項１〜６のうちいずれか１項に記載の粒子挙動解析方法。
【請求項８】
前記粒子受力計算ステップにおいて、
粒子ｉの複数の頂点が物体ｊにめり込んでいる場合に、前記接触判定ステップで算出された前記複数の頂点それぞれのめり込み量を小さくする補正を行い、補正されためり込み量から前記複数の頂点それぞれにかかる接触抗力を算出することを特徴とする請求項１〜７のうちいずれか１項に記載の粒子挙動解析方法。
【請求項９】
密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、
粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得部と、
粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する接触判定部と、
前記接触判定部で算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算部と、
前記粒子受力計算部で算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算部と、
前記粒子運動計算部で算出された粒子の運動を出力する計算結果出力部と、
を有することを特徴とする粒子挙動解析装置。
【請求項１０】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、前記辺ｊ−ｅ１に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項９に記載の粒子挙動解析装置。
【請求項１１】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ１との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項１０に記載の粒子挙動解析装置。
【請求項１２】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０が物体ｊの面と交差する点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０が物体ｊの２つの面それぞれと交差する２つの点、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項９に記載の粒子挙動解析装置。
【請求項１３】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項１２に記載の粒子挙動解析装置。
【請求項１４】
前記粒子受力計算部は、
前記接触点のめり込み量から、ヘルツの式を用いて、前記めり込み方向の接触抗力を算出することを特徴とする請求項９〜１３のうちいずれか１項に記載の粒子挙動解析装置。
【請求項１５】
前記粒子受力計算部は、
前記めり込み方向に垂直な面内での前記接触点の移動量を算出し、算出された移動量から、ミンドリンの式を用いて、前記めり込み方向に垂直な方向の接触抗力を算出することを特徴とする請求項９〜１４のうちいずれか１項に記載の粒子挙動解析装置。
【請求項１６】
前記粒子受力計算部は、
粒子ｉの複数の頂点が物体ｊにめり込んでいる場合に、前記接触判定部で算出された前記複数の頂点それぞれのめり込み量を小さくする補正を行い、補正されためり込み量から前記複数の頂点それぞれにかかる接触抗力を算出することを特徴とする請求項９〜１５のうちいずれか１項に記載の粒子挙動解析装置。
【請求項１７】
請求項１〜８のうちいずれか１項に記載の粒子挙動解析方法の各ステップを情報処理装置に実行させることを特徴とする解析プログラム。
【請求項１】
情報処理装置において、密集状態にある複数の粒子の運動を計算する方法であって、
情報処理装置が、
粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得ステップと、
粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する接触判定ステップと、
前記接触判定ステップで算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算ステップと、
前記粒子受力計算ステップで算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算ステップと、
前記粒子運動計算ステップで算出された粒子の運動を出力する出力ステップと、
を実行することを特徴とする粒子挙動解析方法。
【請求項２】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、前記辺ｊ−ｅ１に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項１に記載の粒子挙動解析方法。
【請求項３】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ１との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項２に記載の粒子挙動解析方法。
【請求項４】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０が物体ｊの面と交差する点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０が物体ｊの２つの面それぞれと交差する２つの点、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項１に記載の粒子挙動解析方法。
【請求項５】
前記接触判定ステップにおいて、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項４に記載の粒子挙動解析方法。
【請求項６】
前記粒子受力計算ステップにおいて、
前記接触点のめり込み量から、ヘルツの式を用いて、前記めり込み方向の接触抗力を算出することを特徴とする請求項１〜５のうちいずれか１項に記載の粒子挙動解析方法。
【請求項７】
前記粒子受力計算ステップにおいて、
前記めり込み方向に垂直な面内での前記接触点の移動量を算出し、算出された移動量から、ミンドリンの式を用いて、前記めり込み方向に垂直な方向の接触抗力を算出することを特徴とする請求項１〜６のうちいずれか１項に記載の粒子挙動解析方法。
【請求項８】
前記粒子受力計算ステップにおいて、
粒子ｉの複数の頂点が物体ｊにめり込んでいる場合に、前記接触判定ステップで算出された前記複数の頂点それぞれのめり込み量を小さくする補正を行い、補正されためり込み量から前記複数の頂点それぞれにかかる接触抗力を算出することを特徴とする請求項１〜７のうちいずれか１項に記載の粒子挙動解析方法。
【請求項９】
密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、
粒子ｉの形状、及び、粒子ｉが接触し得る他の物体ｊの形状をそれぞれ複数の多角形から構成される多面体により定義するデータを取得する取得部と、
粒子ｉの物体ｊに対する接触判定を行って、物体ｊの面にめり込む粒子ｉの頂点、及び、物体ｊの辺にめり込む粒子ｉの辺上の点を、粒子ｉの物体ｊに対する接触点として抽出するとともに、抽出された各接触点のめり込み方向及びめり込み量を算出する接触判定部と、
前記接触判定部で算出されためり込み量から粒子ｉの各接触点にかかる接触抗力を算出する粒子受力計算部と、
前記粒子受力計算部で算出された接触抗力から、ニュートンの運動方程式とオイラーの運動方程式を用いて、粒子ｉの時間進展に伴う運動を算出する粒子運動計算部と、
前記粒子運動計算部で算出された粒子の運動を出力する計算結果出力部と、
を有することを特徴とする粒子挙動解析装置。
【請求項１０】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点と、前記辺ｊ−ｅ１に最も近い前記辺ｉ−ｅ０上の点と、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項９に記載の粒子挙動解析装置。
【請求項１１】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離、及び、前記辺ｉ−ｅ０と前記辺ｊ−ｅ１との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項１０に記載の粒子挙動解析装置。
【請求項１２】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０を接触点として抽出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と、前記辺ｊ−ｅ０が物体ｊの面と交差する点と、を接触点として抽出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｊ−ｅ０に最も近い前記辺ｉ−ｅ０上の点を接触点として抽出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、前記辺ｉ−ｅ０が物体ｊの２つの面それぞれと交差する２つの点、を接触点として抽出し、
（ｅ）その他の場合に、粒子ｉの物体ｊに対する接触はないと判定する
ことを特徴とする請求項９に記載の粒子挙動解析装置。
【請求項１３】
前記接触判定部は、
（ａ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込んでいる場合に、前記頂点ｉ−ｖ０と前記面ｊ−ｆ０との距離を、前記接触点のめり込み量として算出し、
（ｂ）粒子ｉの頂点ｉ−ｖ０が物体ｊの１つの面ｊ−ｆ０にめり込み、かつ、前記頂点ｉ−ｖ０を端点に持つ粒子ｉの辺ｉ−ｅ０が物体ｊの辺ｊ−ｅ０にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出し、
（ｃ）粒子ｉの辺ｉ−ｅ０が物体ｊの１つの辺ｊ−ｅ０のみにめり込んでいる場合に、前記辺ｉ−ｅ０と前記辺ｊ−ｅ０との距離を、前記接触点のめり込み量として算出し、
（ｄ）粒子ｉの辺ｉ−ｅ０が物体ｊの２つの辺ｊ−ｅ０及びｊ−ｅ１にめり込んでいる場合に、物体ｊの面のうち前記各接触点からの距離の合計が最も小さくなる面をめり込み基準面とし、前記各接触点と前記めり込み基準面との距離を、前記各接触点のめり込み量として算出する
ことを特徴とする請求項１２に記載の粒子挙動解析装置。
【請求項１４】
前記粒子受力計算部は、
前記接触点のめり込み量から、ヘルツの式を用いて、前記めり込み方向の接触抗力を算出することを特徴とする請求項９〜１３のうちいずれか１項に記載の粒子挙動解析装置。
【請求項１５】
前記粒子受力計算部は、
前記めり込み方向に垂直な面内での前記接触点の移動量を算出し、算出された移動量から、ミンドリンの式を用いて、前記めり込み方向に垂直な方向の接触抗力を算出することを特徴とする請求項９〜１４のうちいずれか１項に記載の粒子挙動解析装置。
【請求項１６】
前記粒子受力計算部は、
粒子ｉの複数の頂点が物体ｊにめり込んでいる場合に、前記接触判定部で算出された前記複数の頂点それぞれのめり込み量を小さくする補正を行い、補正されためり込み量から前記複数の頂点それぞれにかかる接触抗力を算出することを特徴とする請求項９〜１５のうちいずれか１項に記載の粒子挙動解析装置。
【請求項１７】
請求項１〜８のうちいずれか１項に記載の粒子挙動解析方法の各ステップを情報処理装置に実行させることを特徴とする解析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１１−１７０８３２（Ｐ２０１１−１７０８３２Ａ）
【公開日】平成２３年９月１日（２０１１．９．１）
【国際特許分類】
【出願番号】特願２０１０−１８６３５３（Ｐ２０１０−１８６３５３）
【出願日】平成２２年８月２３日（２０１０．８．２３）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【公開日】平成２３年９月１日（２０１１．９．１）
【国際特許分類】
【出願日】平成２２年８月２３日（２０１０．８．２３）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
[ Back to top ]