(57)【要約】【課題】  本発明は個別要素法に基づいた粉体挙動計算方法であって、特に粒子間相互作用を複数考慮する必要のある、例えば電子写真における二成分現像方式における現像剤挙動のように、粒子間接触だけでなく、静電的及び磁気的相互作用を考慮しなければならない粉体挙動計算を高速に実施可能な粒体挙動計算方法及び該方法を実行するためのプログラムを格納したコンピュータ読み取り可能な記録媒体を提供することを目的とする。【解決手段】  本発明の粒子挙動計算方法によれば、粒子間距離に依存した複数の相互作用毎の有効距離を定めた後、粒子間距離計算を実施する必要のある粒子を選択するための小領域分割処理における小領域の大きさを相互作用の有効距離を考慮して適切に設定することができ、不要な距離計算の少ない、効率の良い粒子挙動計算を実施できる。【選択図】        図６
2005122354.tif 000002
【特許請求の範囲】
【請求項１】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、
相互作用の有効距離の最大値rmaxと最小値rminに応じて、前記小領域の大きさを定める小領域サイズ設定工程と
を有することを特徴とする粒子挙動計算方法。
【請求項２】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
前記計算領域における粒子の密度を算出する粒子密度算出工程と、
粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、
相互作用の有効距離の最大値rmaxと最小値rminと前記粒子の密度に応じて、前記小領域の大きさを定める小領域サイズ設定工程と
を有することを特徴とする粒子挙動計算方法。
【請求項３】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、
前記計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、
前記部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminに応じて、前記小領域の大きさを定める小領域サイズ設定工程と
を有することを特徴とする粒子挙動計算方法。
【請求項４】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
前記計算領域における粒子の密度を算出する粒子密度算出工程と、
粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、
前記計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、
前記部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminと前記粒子の密度に応じて、前記小領域の大きさを定める小領域サイズ設定工程と
を有することを特徴とする粒子挙動計算方法。
【請求項５】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、
前記計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、
前記部分計算領域毎の粒子の密度を算出する粒子密度算出工程と、
前記部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminと前記粒子の密度に応じて、前記小領域の大きさを定める小領域サイズ設定工程と
を有することを特徴とする粒子挙動計算方法。
【請求項６】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
前記計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、
前記部分計算領域毎に、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、
前記部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminに応じて、前記小領域の大きさを定める小領域サイズ設定工程と
を有することを特徴とする粒子挙動計算方法。
【請求項７】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
前記計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、
前記部分計算領域毎に、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、
前記部分計算領域毎の粒子の密度を算出する粒子密度算出工程と、
前記部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminと前記粒子の密度に応じて、前記小領域の大きさを定める小領域サイズ設定工程と
を有することを特徴とする粒子挙動計算方法。
【請求項８】
各粒子の物質間に、粒子間の距離に依存した少なくとも２種類以上の相互作用力を考慮し、粒子の挙動を計算する計算領域を多数の小領域に分割し、分割した小領域内に含まれる粒子の挙動を計算する粒体挙動計算方法において、
前記計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、
前記各部分計算領域に含まれる粒子の個数を計数する粒子計数工程と、
前記各部分計算領域に含まれる粒子の個数に応じて境界線もしくは境界面を修正して再度複数の前記部分計算領域に分割する計算領域再分割工程と
を有することを特徴とする粒体挙動計算方法。
【請求項９】
前記小領域は、２次元の場合正方形、もしくは３次元の場合立方体であって、その１辺の長さＬを、相互作用の有効距離の最大値rmaxと最小値rminの関係がrmin≦Ｌ≦rmaxとなるように設定する請求項１〜８のいずれかに記載の粒子挙動計算方法。
【請求項１０】
前記小領域サイズ設定工程による前記部分計算領域毎の前記小領域の大きさを設定する演算手段を前記部分計算領域毎に割り当てる工程を含む請求項３〜８のいずれかに記載の粒子挙動計算方法。
【請求項１１】
割り当てられた前記各演算手段による演算処理を並行して実行させる工程を含む請求項１０記載の粒子挙動計算方法。
【請求項１２】
前記計算領域再分割工程は、前記各部分計算領域に含まれる粒子の個数又は前記部分計算領域毎の各挙動計算に要する計算時間に応じて境界線もしくは境界面を修正して再度複数の前記部分計算領域に分割する請求項８記載の粒子挙動計算方法。
【請求項１３】
前記計算領域再分割工程は、前記各部分計算領域に含まれる粒子の個数及び前記部分計算領域毎の各挙動計算に要する計算時間に応じて境界線もしくは境界面を修正して再度複数の前記部分計算領域に分割する請求項８記載の粒子挙動計算方法。
【請求項１４】
請求項１〜１３のいずれかに記載の粒子挙動計算方法を実行するためのプログラムを格納したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は粒子挙動計算方法に関し、詳細には粉体間に働く複数の相互作用を考慮する必要のある粉体挙動、例えば電子写真式画像形成プロセスにおける現像プロセスでの二成分現像剤の挙動などの計算機シミュレーションに関する。
【背景技術】
【０００２】
近年、粉体や粒体等の粒子の挙動シミュレーションについては、個別要素法あるいは離散粒子法と呼ばれる方法が普及されており、特にCundallらによって提案された個別要素法に基礎をおいた粒子挙動シミュレータを利用した粉体挙動解析が実施されるようになってきた。これは、コンピュータ性能が飛躍的に進歩し、従来のスーパーコンピュータと同等あるいはそれ以上の計算性能を有するコンピュータが、大量に、しかも安価に利用できるようになってきたことが背景として考えられる。
【０００３】
一方、個別要素法に基づいた粒子挙動計算アルゴリズムでは、粒子数が多くなると、計算量が膨大になり、いくら計算機の性能が向上したとはいっても、実際の系と同等の粒子数での計算を実行することは困難な場合が多い。粉体挙動解析方法では、壁または角との接触判定の実施回数を軽減し、かつ角との接触での挙動を適切に処理する方法が重要である。また、相互作用は粒子間接触のみならず、静電気力や磁気相互作用を考慮しなければならない。
【０００４】
そこで、従来の粒子挙動シミュレーション方法として、特許文献１には、容器内粒子の混合攪拌に伴って粒子間又は粒子−容器壁面間に作用する付着力を、粒子表面及び容器壁面の表面粗さを用いて算出する方法が提案されている。また、特許文献１における計算対象として電子写真におけるキャリアとトナーからなる２成分現像方式の現像剤における帯電を考慮している。また、特許文献２には、接触力だけでなく磁気相互作用を考慮し、少なくとも２つの磁気搬送ローラ間の磁性粉体の受け渡しの挙動を、粒子要素法と磁場解析とを組み合わせて数値解析することによって磁性粉体の搬送量と磁性粉体の攪拌率を精度良く把握できる磁性粉体の挙動予測方法が提案されている。
【特許文献１】特開２００２−１０９４４５号公報
【特許文献２】特開２００３−１９５６４２号公報
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかしながら、特許文献１は、計算速度の向上を考慮したものではなく、更に粒径の大きくことなる粒子（キャリアとトナー）の混合状態のシミュレーションであることに対して考慮されていない。また、特許文献２によれば、接触力だけでなく磁気相互作用を考慮した粉体挙動シミュレータに関わる計算方法については公知のものを利用しており、計算速度を向上させるような機能については考慮されていない。
【０００６】
本発明はこれらの問題点を解決するためのものであり、個別要素法に基づいた粉体挙動計算方法であって、特に粒子間相互作用を複数考慮する必要のある、例えば電子写真における二成分現像方式における現像剤挙動のように、粒子間接触だけでなく、静電的及び磁気的相互作用を考慮しなければならない粉体挙動計算を高速に実施可能な粒体挙動計算方法及び該方法を実行するためのプログラムを格納したコンピュータ読み取り可能な記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００７】
前記問題点を解決するために、本発明の粒子挙動計算方法によれば、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、相互作用の有効距離の最大値rmaxと最小値rminに応じて、小領域の大きさを定める小領域サイズ設定工程とを有している。よって、不要な距離計算の少ない、効率の良い粒子挙動計算を実施できる。
【０００８】
次に、本発明の粒子挙動計算方法によれば、計算領域における粒子の密度を算出する粒子密度算出工程と、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、相互作用の有効距離の最大値rmaxと最小値rminと粒子の密度に応じて、小領域の大きさを定める小領域サイズ設定工程とを有している。よって、不要な距離計算の少ない、効率の良い粒子挙動計算方法を実施できる。
【０００９】
また、本発明の粒子挙動計算方法によれば、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、記計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminに応じて、小領域の大きさを定める小領域サイズ設定工程とを有している。よって、計算対象の不均一性に依存しない、効率の良い粒子挙動計算を実施できる。
【００１０】
更に、本発明の粒子挙動計算方法によれば、計算領域における粒子の密度を算出する粒子密度算出工程と、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminと前記粒子の密度に応じて、小領域の大きさを定める小領域サイズ設定工程とを有している。よって、計算対象の不均一性に依存しない、効率の良い粒子挙動計算を実施できる。
【００１１】
また、本発明の粒子挙動計算方法によれば、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、部分計算領域毎の粒子の密度を算出する粒子密度算出工程と、部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminと粒子の密度に応じて、小領域の大きさを定める小領域サイズ設定工程とを有している。よって、計算対象の不均一性に依存しない、効率の良い粒子挙動計算を実施できる。
【００１２】
更に、本発明の粒子挙動計算方法によれば、計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、部分計算領域毎に、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminに応じて、小領域の大きさを定める小領域サイズ設定工程とを有している。よって、計算対象の不均一性に依存しない、効率の良い粒子挙動計算を実施できる。
【００１３】
また、本発明の粒子挙動計算方法によれば、計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、部分計算領域毎に、粒子間に作用する相互作用毎の有効距離を定める有効距離設定工程と、部分計算領域毎の粒子の密度を算出する粒子密度算出工程と、部分計算領域毎に、相互作用の有効距離の最大値rmaxと最小値rminと粒子の密度に応じて、小領域の大きさを定める小領域サイズ設定工程とを有している。よって、計算対象における粒子の分布状態、計算対象の不均一性に依存しない、効率の良い粒子挙動計算を実施できる。
【００１４】
更に、本発明の粒子挙動計算方法によれば、計算領域の全体を複数の部分計算領域に分割する境界線あるいは境界面を定める計算領域分割工程と、各部分計算領域に含まれる粒子の個数を計数する粒子計数工程と、各部分計算領域に含まれる粒子の個数に応じて境界線もしくは境界面を修正して再度複数の部分計算領域に分割する計算領域再分割工程とを有している。よって、各部分計算領域毎の粒子の個数を均一になるように計算領域全体を分割することができ、より効率の良い粒子挙動計算を実施できる。
【００１５】
また、小領域は、２次元の場合正方形、もしくは３次元の場合立方体であるので、各小領域に含まれる粒子のチェックが容易となる。更に、小領域の１辺の長さＬを、相互作用の有効距離の最大値rmaxと最小値rminの関係がrmin≦Ｌ≦rmaxとなるように適切に設定することにより、不要な距離計算の少ない、効率の良い粒子挙動計算を実施できる。
【００１６】
更に、小領域サイズ設定工程による部分計算領域毎の小領域の大きさを設定する演算手段を部分計算領域毎に割り当てる工程を含むことにより、また割り当てられた各演算手段による演算処理を並行して実行させることにより、より効率の良い粒子挙動計算を実施できる。
【００１７】
また、計算領域再分割工程は、各部分計算領域に含まれる粒子の個数及び／又は部分計算領域毎の各挙動計算に要する計算時間に応じて境界線もしくは境界面を修正して再度複数の部分計算領域に分割することにより、各部分計算領域の計算時間を同時に実行させることができ、より効率良く粒子挙動計算を実施することができる。
【００１８】
更に、別の発明として、上記記載の粒子挙動計算方法を実行するためのプログラムを格納したコンピュータ読み取り可能な記録媒体に特徴がある。よって、既存のシステムを変えることなく、粒子挙動計算システムを汎用的に構築することができる。
【発明の効果】
【００１９】
本発明の粒子挙動計算方法によれば、粒子間距離に依存した複数の相互作用毎の有効距離を定めた後、粒子間距離計算を実施する必要のある粒子を選択するための小領域分割処理における小領域の大きさを相互作用の有効距離を考慮して適切に設定することができ、不要な距離計算の少ない、効率の良い粒子挙動計算を実施できる。
【発明を実施するための最良の形態】
【００２０】
本発明の粒子挙動計算方法は、粒子間に作用する相互作用毎の有効距離を定め、相互作用の有効距離の最大値rmaxと最小値rminに応じて、粒子の挙動を計算する計算領域を多数に分割した小領域の大きさを定める。
【実施例】
【００２１】
はじめに、本発明が適用する個別要素法に基づく粒体挙動計算の原理について説明すると、この個別要素法とは粉体を形成する各粒子に作用する力を全て求め、その作用力による各粒子の運動方程式を逐次解くことで、粒体挙動を計算するものである。
【００２２】
粒子に作用する力としては、重力のような外力の他、壁と粒子間の接触及び粒子間接触による接触力が考えられる。その力の大きさは、例えば図１に示すような接触モデルに従って計算される。図１の（ａ）は、法線方向に作用する接触力のモデルで、バネ１１とダッシュポット１２により構成され、作用力の大きさは、法線方向の速度と粒子の接触状態により決定される。
【００２３】
各粒子の接触状態を計算するためには、各粒子間の距離計算を実施する必要がある。全ての粒子間で距離計算を行うと、粒子数ｎの場合、ｎ（ｎ−１）／２回の距離計算が必要となり、粒子数ｎが大きくなると、その計算量は膨大なものとなってしまう。
【００２４】
このような不具合を解消するために、従来からいろいろな方法が提案され、実施されている。その中で、一般的でかつ効率のよい手法として、次のような方法がある。その方法とは、計算対象領域を小領域に分割し、接触判定（距離計算）を実施する粒子を、着目する粒子のある小領域とその近傍の小領域内に存在する粒子に制限する方法である。この方法を図２に従って説明すると、図２において、着目粒子２１を含む小領域２４及びその周辺の８セルを含む領域２５に含まれる粒子２２を距離計算の対象粒子とする。それ以外の粒子２３は、着目粒子２１と接触する可能性はないので、距離計算を実施する必要はない。図２は、計算領域を２次元で考えているので、格子間隔Ｌを対象粒子の直径としたＸ１〜Ｘ５、Ｙ１〜Ｙ４で形成するメッシュ分割したセルを小領域としている。このような手法は、近傍セル法、グリッド法などと呼ばれる方法である。
【００２５】
一方、電子写真における現像剤挙動の計算を実施する場合、上記のような粒子間接触による作用力のほかに、静電的な作用力や磁気的な作用力を考慮する必要がある。図３では、電荷量Ｑ１、Ｑ２をそれぞれ持った２つの粒子３１、３２が距離ｒの位置にある場合を示している。２つの粒子３１，３２は接触していないが、下記の式で示すようなクーロン力が作用する。
【００２６】
｜Ｆ１｜＝｜Ｆ２｜＝Ｑ１・Ｑ２／４πε０ｒ２
【００２７】
従って、静電的な作用力など、接触力以外の粒子間相互作用を考慮した粒子挙動計算を実施する場合には、上述した近傍セル法のように、メッシュ間隔Ｌを粒子間隔程度にする方法では不具合が生じる。
【００２８】
図４は、例えば上記式の粒子間距離ｒに対して、有効距離を定め、その値をメッシュ間隔Ｌとした場合である。着目粒子４１を含むセル４２とその近傍８セルの領域４３内の粒子のみ距離計算の対象とし、それ以外の粒子４４等は、距離計算を行う必要はない。これは、領域４３の外部にある粒子４４と着目粒子４１との距離は、常に静電力の有効距離より大きくなるためである。
【００２９】
ここで、静電的相互作用の有効距離を粒子直径ｄの８倍とすることを考える。このとき、近傍セル法においてメッシュ間隔ＬをＬ＝ｄとした場合、接触相互作用のみ考慮する場合には、９セルを探索すれば十分であるが、静電相互作用（有効距離８ｄとして）を考慮する場合には、２次元の場合で２８９セル（１７×１７）、３次元の場合だと４、９１３セル（１７×１７×１７）に含まれる粒子との距離計算を実施する必要があり、その探索が必要なセル数の増加がメモリを圧迫し、計算時間を増大させる原因となってしまう場合がある。
【００３０】
一方、Ｌ＝８ｄと静電的相互作用の有効距離に一致させた場合では、探索するセル数は、図４に示すように、２次元の場合で９セル、３次元の場合でも２７セルでよい。この場合では、セル探索が負担にはならないが、距離計算対象粒子の探索領域面積Ｓは、２次元の場合で、Ｓ＝３Ｌ×３Ｌ＝２４ｄ×２４ｄとなる。Ｌ＝ｄの場合では、Ｓ＝１７Ｌ×１７Ｌ＝１７ｄ×１７ｄであり、探索領域が大きくなるため、無駄な距離計算が多く発生してしまう。特に、粒子密度が高い場合には不要な距離計算が増加し、これが計算時間を増大させることになる。全体的な粒子密度が低くとも、粒子が局在する場合なども考えられ、メッシュ間隔Ｌの最適設定はその計算対象である粒子の存在状態に強く依存するが、少なくとも、粒径ｄ、静電的相互作用の有効距離を粒径の８倍と設定する場合で、小領域を等間隔のメッシュ分割によるセルとすれば、メッシュ間隔Ｌはｄ≦Ｌ≦８ｄの範囲に最適値があることになる。よって、上記のように粒子間距離に依存した複数の相互作用が存在する場合の小領域の大きさ（上記ではメッシュ間隔Ｌ）の設定方法を相互作用の有効距離の最小値と最大値（上記では、ｄと８ｄ）に応じて定める。また、小領域を等間隔のメッシュにより形成する場合において、そのメッシュ間隔Ｌを相互作用の有効距離の最小値rminと最大値rmaxとの中間値としてｄ≦Ｌ≦８ｄとなるように定める。
【００３１】
ここで、従来の近傍セル法による接触による相互作用のみを考慮した粒子挙動計算処理の処理フローを図５に示す。同図において、先ず計算に必要な各種物理パラメータや粒子の初期配置などの計算条件を読み込み（ステップＳ１０１）、次にメッシュを生成し（ステップＳ１０２）、セル内に存在する粒子をチェックする（ステップＳ１０３）。着目粒子との接触判定の候補として、着目している粒子の存在するセルの近傍セル内の粒子を選定し（ステップＳ１０４）、各粒子に作用する接触力を算出・積算する（ステップＳ１０５）。ステップ１０６では、全粒子の計算の終了をチェックし、終了していない場合（ステップＳ１０６；ＮＯ）、ステップＳ１０４に戻り、終了していれば（ステップＳ１０６；ＹＥＳ）、各粒子の運動方程式を解き、加速度、速度及び変位を計算し（ステップＳ１０７）、それぞれの粒子位置を更新する（ステップＳ１０８）。ステップＳ１０９では、最終時間ステップであるかどうかをチェックし、最終時間ステップに達していれば（ステップＳ１０９；ＹＥＳ）計算を終了し、それ以外は（ステップＳ１０９；ＮＯ）、時間ステップを更新し（ステップＳ１１０）、セル内の粒子チェックのステップＳ１０３に戻る。
【００３２】
図６は本発明の第１の実施例に係る粒子挙動計算処理を示すフローチャートである。同図において、小領域は、等間隔のメッシュ分割により、２次元の場合では正方形、３次元の場合では、立法体とした。同図において、先ず計算に必要な各種物理パラメータや粒子の初期配置などの計算条件を読み込み（ステップＳ２０１）、次に相互作用の有効距離を設定し（ステップＳ２０２）、メッシュサイズＬを設定する（ステップＳ２０３）。次に、メッシュを生成し（ステップＳ２０４）、セル内に存在する粒子をチェックする（ステップＳ２０５）。着目粒子との接触判定の候補として、着目している粒子の存在するセルの近傍セル内の粒子を選定し（ステップＳ２０６）、各粒子に作用する接触力を算出・積算する（ステップＳ２０７）。ステップ２０８では、全粒子の計算の終了をチェックし、終了していない場合（ステップＳ２０８；ＮＯ）、ステップＳ２０６に戻り、終了していれば（ステップＳ２０８；ＹＥＳ）、各粒子の運動方程式を解き、加速度、速度及び変位を計算し（ステップＳ２０９）、それぞれの粒子位置を更新する（ステップＳ２１０）。ステップＳ２１１では、最終時間ステップであるかどうかをチェックし、最終時間ステップに達していれば（ステップＳ２１１；ＹＥＳ）計算を終了し、それ以外は（ステップＳ２１１；ＮＯ）、時間ステップを更新し（ステップＳ２１２）、セル内の粒子チェックのステップＳ２０５に戻る。
【００３３】
このように、図５に示した従来の処理と異なるステップは、ステップＳ２０２における相互作用の有効距離を設定することと、ステップＳ２０３でのメッシュサイズＬを設定することである。このように、相互作用の有効距離を設定したのちに、メッシュサイズＬを設定することで、相互作用有効距離を考慮してメッシュサイズＬを設定することが可能となっている。また、相互作用の有効距離は、接触では計算対象である粒子の最大粒子径であり、静電相互作用のように、数式上は無限遠方まで作用する場合には、必要に応じてカットオフ距離を設定して、それを有効距離とする方法などが最も単純な方法として考えられる。この有効距離の設定方法については、計算対象あるいは解析の目的などに応じて適切に定めるべきものであり、本発明においては、有効距離の設定方法を制限するものではない。また、図５のステップＳ１０５での接触力算出が、図６では相互作用力算出（ステップＳ２０７）となっており、接触力に限らず、複数の相互作用の計算を全て実施している。なお、計算対象領域を２次元の場合には正方形に、３次元の場合には立方体の小領域に分割することにより、粒子が存在する小領域を容易に見出すことができ、かつ相互作用の有効距離と小領域の大きさを定める値として正方形もしくは立法体の一辺の長さをメッシュサイズＬとして明確に関連付けることができる。
【００３４】
図７は本発明の第２の実施例に係る粒子挙動計算処理を示すフローチャートである。図６の処理フローとほぼ同じであるが、ステップＳ３０２のような粒子密度算出処理が追加されている。すなわち、ステップＳ３０５でのメッシュサイズＬを設定する前に、ステップＳ３０２での粒子密度算出処理と、ステップＳ３０３での相互作用の有効距離設定処理を実施するため、粒子の密集度並びに相互作用の有効距離の最小値と最大値を考慮してメッシュサイズＬを定めることができるような計算フローとなっているので、より最適なメッシュサイズＬが得られる。なお、本実施例では、粒子密度を算出しているが、粒子の密集度を表すものであれば、粒子密度の算出以外の値を用いてもよい。
【００３５】
計算対象として、磁場中に置かれた磁性粒子、例えば電子写真における２成分現像剤などを考える場合、磁気的な相互作用によって、図８に示すように磁石８１を配置している近傍に粒子が局在するような状況が想定される。このような場合、磁石８の近傍とそれ以外とでは、メッシュサイズを変更した法がより効果的であることがある。そのような場合において、領域を分割した処理が効果的である。
【００３６】
なお、計算対象領域を２次元の場合には正方形に、３次元の場合には立方体の小領域に分割することにより、粒子が存在する小領域を容易に見出すことができ、かつ相互作用の有効距離と小領域の大きさを定める値として正方形もしくは立法体の一辺の長さをメッシュサイズＬとして明確に関連付けることができる。
【００３７】
図９は本発明の第３の実施例に係る粒子挙動計算処理を示すフローチャートである。本実施例は、上記の第１、第２の実施例における計算領域内でその粒子分布が局在している場合などに効果的な処理が可能なように計算領域を複数に分割して処理する方法を提供するものである。同図において、先ず計算に必要な各種物理パラメータや粒子の初期配置などの計算条件を読み込み（ステップＳ４０１）、次に相互作用の有効距離を設定する（ステップＳ４０２）。次に、計算領域を２つの領域に分割する境界線（２次元の場合）あるいは境界面（３次元の場合）を定め、計算領域Ａと計算領域Ｂとに分割する（ステップＳ４０３）。引き続いて、分割した計算領域毎に、後述する粒子挙動の計算を実施する(ステップＳ４０４、Ｓ４０５)。そして、時間ステップをチェックし（ステップＳ４０６）、最終時間ステップであれば（ステップＳ４０６；ＹＥＳ）計算を終了し、それ以外であれば（ステップＳ４０６；ＮＯ）時間ステップを更新し（ステップＳ４０７）、最初の粒子挙動計算のステップＳ４０４に戻る。
【００３８】
ここで、図９での分割した領域での粒子挙動計算処理を図１０のフローチャートに従って説明すると、先ず部分計算領域内でのメッシュサイズＬを設定し（ステップＳ５０１）、それに応じたメッシュを生成する（ステップＳ５０２）。これによって、領域を分割しているのでそれぞれに応じた適切なメッシュサイズを設定可能となる。引き続き、セル内の粒子をチェックし（ステップＳ５０３）着目粒子との距離算出候補を選定し（ステップＳ５０４）、相互作用を算出・積算する（ステップＳ５０５）。次に、部分計算領域内における全粒子の計算終了の確認を実施し（ステップＳ５０６）、全粒子の計算が終了していなければ（ステップＳ５０６；ＮＯ）、ステップＳ５０４に戻る。全粒子の作用力計算が終了していれば（ステップＳ５０６；ＹＥＳ）、その作用力から各粒子の運動方程式を解き、加速度・速度及び変位を求め（ステップＳ５０７）、部分計算領域内の各粒子の位置を更新する（ステップＳ５０８）。
【００３９】
よって、第３の実施例によれば、粒子の局在など、局所的な状態の変化による影響を回避し、分割された部分領域毎に最適なメッシュサイズＬを定めることが可能となる。なお、粒子の局在等がある場合に限るものではなく、粒子分布が均一な場合であってもよい。また、計算対象領域を２次元の場合には正方形に、３次元の場合には立方体の小領域に分割することにより、粒子が存在する小領域を容易に見出すことができ、かつ相互作用の有効距離と小領域の大きさを定める値として正方形もしくは立法体の一辺の長さをメッシュサイズＬとして明確に関連付けることができる。
【００４０】
図１１は本発明の第４の実施例に係る粒子挙動計算処理を示すフローチャートである。図１１は、図９とほぼ同じ処理フローであるが、粒子密度算出（ステップＳ６０２）が追加されている。ステップＳ６０５，Ｓ６０６における、分割された領域での粒子挙動計算の処理は図９の第３の実施例と同じく、図１０の処理と同じである。本実施例において、ステップＳ６０２で部分計算領域内メッシュサイズＬ設定に先立って粒子密度が算出されているので、粒子の密集度に応じてメッシュサイズＬを適切に設定することができる。他のステップの動作については、第３の実施例として示した図９と同様なので、ここでは説明を省略する。
【００４１】
なお、計算対象領域を２次元の場合には正方形に、３次元の場合には立方体の小領域に分割することにより、粒子が存在する小領域を容易に見出すことができ、かつ相互作用の有効距離と小領域の大きさを定める値として正方形もしくは立法体の一辺の長さをメッシュサイズＬとして明確に関連付けることができる。
【００４２】
次に、図１２は分割した領域での別の粒子挙動計算処理を示すフローチャートである。同図に示す分割した領域での粒子挙動計算処理では、ステップＳ７０１で部分計算領域内の粒子密度を算出して、メッシュサイズＬを設定された相互作用の有効距離と部分計算領域内の粒子の密集度に応じて設定することができる。これによって、より最適なメッシュサイズＬを定めることができる。なお、それ以降の処理（ステップＳ７０２〜Ｓ７０９）は、図１０の処理（ステップＳ５０１〜Ｓ５０８）と同一であるので説明を省略する。
【００４３】
なお、計算対象領域を２次元の場合には正方形に、３次元の場合には立方体の小領域に分割することにより、粒子が存在する小領域を容易に見出すことができ、かつ相互作用の有効距離と小領域の大きさを定める値として正方形もしくは立法体の一辺の長さをメッシュサイズＬとして明確に関連付けることができる。
【００４４】
また、図１３は分割した領域での更に別の粒子挙動計算処理を示すフローチャートである。同図に示す分割した領域での粒子挙動計算処理では、ステップＳ８０１で部分計算領域内における相互作用の有効距離を設定する。これによって分割された領域毎の状態（粒子の配置や磁場の強弱や電荷の有無等）に応じた相互作用有効距離が設定でき、計算領域全体の状態が均質でない場合でも適切な相互作用距離が設定できる。次に、部分計算領域内のメッシュサイズＬを定め（ステップＳ８０２）メッシュを生成する（ステップＳ８０３）。メッシュ生成前に部分領域に適した相互作用の有効距離を定め、それに応じたメッシュサイズＬが定められているので、部分領域の粒子挙動計算に最適なメッシュが生成されることになる。続いて、セル内の粒子をチェックする（ステップＳ８０４）から着目粒子と距離計算対象粒子選定を設定し（ステップＳ８０５）、各粒子間距離に応じた相互作用力（接触している場合、接触力も含む）を算出し、各粒子毎の作用力を積算する（ステップＳ８０６）。次に、部分計算領域内全ての粒子計算が完了していることを確認する（ステップＳ８０７）。すなわち、部分計算領域内の全ての粒子計算が完了していれば（ステップＳ８０７；ＹＥＳ）、各粒子の運動方程式を解き、加速度、速度及び変位を求め（ステップＳ８０８）、そうでない場合には（ステップＳ８０７；ＹＥＳ）、ステップＳ８０５の着目粒子と距離計算対象粒子選定に戻る。ステップＳ８０８で変位が求めた後、ステップＳ８０９にて部分領域内の各粒子の位置を更新する。
【００４５】
なお、計算対象領域を２次元の場合には正方形に、３次元の場合には立方体の小領域に分割することにより、粒子が存在する小領域を容易に見出すことができ、かつ相互作用の有効距離と小領域の大きさを定める値として正方形もしくは立法体の一辺の長さをメッシュサイズＬとして明確に関連付けることができる。
【００４６】
更に、図１４は更に別の分割した領域での粒子挙動計算処理を示すフローチャートである。同図において、先ず部分計算領域内における粒子密度を算出し（ステップＳ９０１）、次に相互作用の有効距離を設定する（ステップＳ９０２）。その後に部分計算領域内のメッシュサイズＬを定め（ステップＳ９０３）、メッシュを生成している（ステップＳ９０４）。このように、メッシュ生成前に部分領域に適した相互作用の有効距離を部分計算領域内の粒子密度や部分計算領域の状況に応じて適切に定められた相互作用の有効距離からメッシュサイズＬを定めているので、部分領域の粒子挙動計算に最適なメッシュが生成されることになる。続いて、セル内の粒子をチェックし（ステップＳ９０５）、着目粒子と距離計算対象粒子選定を設定し（ステップＳ９０６）、各粒子間距離に応じた相互作用力（接触している場合、接触力も含む）を算出し、各粒子毎の作用力を積算する（ステップＳ９０７）。次に、部分計算領域内全ての粒子計算が完了していることを確認する（ステップＳ９０８）。すなわち、部分計算領域内の全ての粒子計算が完了していれば（ステップＳ９０８；ＹＥＳ）、各粒子の運動方程式を解き、加速度、速度及び変位を求め（ステップＳ９０９）、そうでない場合には（ステップＳ９０８；ＮＯ）、ステップＳ９０５の着目粒子と距離計算対象粒子選定に戻る。そして、ステップＳ９０９で変位が求めた後、ステップＳ９１０にて部分領域内の各粒子の位置を更新する。
【００４７】
なお、計算対象領域を２次元の場合には正方形に、３次元の場合には立方体の小領域に分割することにより、粒子が存在する小領域を容易に見出すことができ、かつ相互作用の有効距離と小領域の大きさを定める値として正方形もしくは立法体の一辺の長さをメッシュサイズＬとして明確に関連付けることができる。
【００４８】
図１５は本発明の第５の実施例に係る粒子挙動計算処理を示すフローチャートである。同図において、先ず計算に必要な各種物理パラメータや粒子の初期配置などの計算条件を読み込み（ステップＳ１００１）、次に粒子密度の算出（ステップＳ１００２）、相互作用の有効距離設定（ステップＳ１００３）を実施した後、計算領域を分割する（ステップＳ１００４）。ここでは領域Ａと領域Ｂの２つの部分計算領域に分割している。分割後、分割した領域毎に含まれる粒子数をカウントし（ステップＳ１００５，Ｓ１００６）、それぞれのカウント数ＮＡ，ＮＢの差ＮＡ−ＮＢを所定値と比較してその差が所定値以下であれば（ステップＳ１００７；ＹＥＳ）、部分領域内毎での粒子挙動計算をそれぞれ行い（ステップＳ１００８，Ｓ１００９）、カウント差が所定値以下でない場合（ステップＳ１００７；ＮＯ）、ステップＳ１００４での計算領域分割処理に戻り、カウント差が少なくなるように分割計算領域に再度分割する。なお、部分領域内の粒子挙動計算は、例えば図１０に示した処理フローによって実施される。そして、粒子計算終了後、時間ステップを確認し（ステップＳ１０１０）、時間ステップが最終時間ステップまで到達していれば（ステップＳ１０１０；ＹＥＳ）、処理を終了し、それ以外は、時間ステップを更新し（ステップＳ１０１１）、最初の部分領域内粒子計算処理のステップＳ１００８に戻る。このような処理によって、領域Ａと領域Ｂに含まれる粒子個数をほぼ均等に保つことができ、領域毎の差異を少なくすることができる。
【００４９】
図１６は本発明の第６の実施例に係る粒子挙動計算処理を示すフローチャートである。同図において、先ず計算に必要な各種物理パラメータや粒子の初期配置などの計算条件を読み込み（ステップＳ１１０１）、次に粒子密度の算出（ステップＳ１１０２）、相互作用の有効距離設定（ステップＳ１１０３）を実施した後、計算領域を分割する（ステップＳ１１０４）。ここでは領域Ａと領域Ｂの２つの部分計算領域に分割している。分割後、分割した領域毎に含まれる粒子数をカウントし（ステップＳ１１０５，Ｓ１１０６）、それぞれのカウント数ＮＡ，ＮＢの差ＮＡ−ＮＢを所定値と比較してその差が所定値以下であれば（ステップＳ１１０７；ＹＥＳ）、部分領域内毎での粒子挙動計算をそれぞれ行い（ステップＳ１１０８，Ｓ１１０９）、カウント差が所定値以下でない場合（ステップＳ１１０７；ＮＯ）、ステップＳ１１０４での計算領域分割処理に戻る。なお、部分領域内の粒子挙動計算は、例えば図１０に示した処理フローによって実施される。そして、粒子計算終了後、時間ステップを確認し（ステップＳ１１１０）、時間ステップが最終時間ステップまで到達していれば（ステップＳ１１１０；ＹＥＳ）、処理を終了し、それ以外は、時間ステップを更新し（ステップＳ１１１１）、ステップＳ１１０４での計算領域分割処理に戻る。このような処理によって、領域Ａと領域Ｂに含まれる粒子個数を各時間ステップでほぼ均等に保つことができ、計算過程での粒子の移動の影響を受けることなく、部分領域毎の差異を少なくすることができる。
【００５０】
図１７は本発明の第７の実施例に係る粒子挙動計算処理を示すフローチャートである。同図において、先ず計算に必要な各種物理パラメータや粒子の初期配置などの計算条件を読み込み（ステップＳ１２０１）、次に粒子密度の算出（ステップＳ１２０２）、相互作用の有効距離設定（ステップＳ１２０３）を実施した後、計算領域を分割する（ステップＳ１２０４）。分割後、分割した領域毎に粒子挙動計算手段をそれぞれの部分領域に割り当てる（ステップＳ１２０５）。そして、部分計算領域の数だけの粒子挙動計算手段を備えている場合には、各粒子挙動計算手段によりそれぞれの部分領域での粒子挙動計算を並行して実施する（ステップＳ１２０６）。なお、粒子挙動計算手段が１つしかない場合には、図１６の処理フローのように領域Ａ内の粒子挙動計算のステップＳ１１０８の後に、領域Ｂ内の粒子挙動計算のステップＳ１１０９を実施する逐次処理となるが、複数の粒子挙動計算手段が利用可能な場合には、領域Ａ内と領域Ｂ内の粒子挙動計算を同時に実施することができる。このように並列処理が実施できることで、全体の粒子挙動計算時間を短縮することが可能である。このように、粒子計算処理を並行して実施するためには、粒子計算処理手段を複数個備えた装置が必要となるが、複数のＣＰＵを備えたパーソナルコンピュータなども普及しており、容易に実現可能である。そして、部分計算領域内のそれぞれの粒子挙動計算が終了したのち、時間ステップを確認し（ステップＳ１２０７）、時間ステップが最終時間ステップに到達していれば（ステップＳ１２０７；ＹＥＳ）、処理を終了し、それ以外は（ステップＳ１２０７；ＮＯ）、時間ステップを更新し（ステップＳ１２０８）、ステップＳ１２０６での粒子挙動計算へと戻る。なお、領域内粒子挙動計算は、例えば、例えば図１０に示した処理フローによって実施される。
【００５１】
図１８は本発明の第８の実施例に係る粒子挙動計算処理を示すフローチャートである。同図において、先ず計算に必要な各種物理パラメータや粒子の初期配置などの計算条件を読み込み（ステップＳ１３０１）、次に粒子密度の算出（ステップＳ１３０２）、相互作用の有効距離設定（ステップＳ１３０３）を実施した後、計算領域を分割する（ステップＳ１３０３）。分割後、分割した領域毎に粒子挙動計算を並行して実施する（ステップＳ１３０６）が、その計算時間を計測するためのステップ（計測の開始と終了）を粒子挙動計算ステップの前後に設けてある（ステップＳ１３０５，Ｓ１３０６）。部分領域における粒子挙動計算が全て完了し、それぞれの計算時間計測が終了した後（ステップＳ１３０７）、時間ステップを確認し（ステップ１３０８）、最終時間ステップに到達している場合（ステップＳ１３０８；ＹＥＳ）処理を終了し、そうでない場合には（ステップＳ１３０８；ＮＯ）、時間ステップを更新する（ステップＳ１３０９）へ進み、部分領域における計算時間を比較する（ステップＳ１３１０）。部分領域毎の計算時間の差が所定値以下であれば（ステップＳ１３１０；ＹＥＳ）、そのままステップＳ１３０５〜Ｓ１３０７の各領域毎の粒子挙動計算を行うが、部分領域毎の計算時間の差が所定値以上である場合（ステップＳ１３１０；ＮＯ）、ステップＳ１３０４の領域分割へと戻り、計算時間の差を小さくするような修正を施す。このように、各部分領域における挙動計算を並行して実施する場合に、それぞれの計算時間を均一化する処理を実施することで、全体の計算時間を短縮することが可能となる。
【００５２】
次に、図１９は本発明の粒子挙動計算方法の各実施例を実行するプログラムを起動するための具体的な装置の構成を示すブロック図である。つまり、同図は上記実施例における粒子挙動計算方法によるソフトウェアを実行するマイクロプロセッサ等から構築されるハードウェアを示すものである。同図において、粒子挙動計算システムはインターフェース（以下Ｉ／Ｆと略す）５１、ＣＰＵ５２、ＲＯＭ５３、ＲＡＭ５４、表示装置５５、ハードディスク５６、キーボード５７及びＣＤ−ＲＯＭドライブ５８を含んで構成されている。また、汎用の処理装置を用意し、ＣＤ−ＲＯＭなどの読取可能な記録媒体５９には、本発明の粒子挙動計算方法を実行するプログラムが記憶されている。更に、Ｉ／Ｆ５１を介して外部装置から制御信号が入力され、キーボード５７によって操作者による指令又は自動的に本発明のプログラムが起動される。そして、ＣＰＵ５２は当該プログラムに従って上述の粒子挙動計算方法に伴う計算処理を施し、その処理結果をＲＡＭ５４やハードディスク５６等の記憶装置に格納し、必要により表示装置５５などに出力する。以上のように、本発明の粒子挙動計算方法を実行するプログラムが記録した記録媒体を用いることにより、既存のシステムを変えることなく、粒子挙動計算システムを汎用的に構築することができる。
【００５３】
なお、本発明は上記実施例に限定されるものではなく、特許請求の範囲内の記載であれば多種の変形や置換可能であることは言うまでもない。
【図面の簡単な説明】
【００５４】
【図１】接触モデルを示す概略図である。
【図２】小領域に分割した計算対象領域で接触判定を実施する粒子を、着目する粒子のある小領域とその近傍の小領域内に存在する粒子に制限する様子を示す概略図である。
【図３】異なる電荷量を持った２つの粒子が距離ｒの位置にある場合の静電的な作用力や磁気的な作用力を示す概略図である。
【図４】メッシュ間隔Ｌを有効距離とした場合、小領域に分割した計算対象領域で接触判定を実施する粒子を、着目する粒子のある小領域とその近傍の小領域内に存在する粒子に制限する様子を示す概略図である。
【図５】従来の近傍セル法による接触による相互作用のみを考慮した粒子挙動計算処理を示すフローチャートである。
【図６】本発明の第１の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図７】本発明の第２の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図８】磁石８１を配置している近傍に粒子が局在するような状況を示す概略図である。
【図９】本発明の第３の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図１０】分割した領域での粒子挙動計算処理を示すフローチャートである。
【図１１】本発明の第４の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図１２】分割した領域での別の粒子挙動計算処理を示すフローチャートである。
【図１３】分割した領域での更に別の粒子挙動計算処理を示すフローチャートである。
【図１４】分割した領域での更に別の粒子挙動計算処理を示すフローチャートである。
【図１５】本発明の第５の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図１６】本発明の第６の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図１７】本発明の第７の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図１８】本発明の第８の実施例に係る粒子挙動計算処理を示すフローチャートである。
【図１９】本発明の粒子挙動計算方法の各実施例を実行するプログラムを起動するための具体的な装置の構成を示すブロック図である。
【符号の説明】
【００５５】
１１；バネ、１２；ダッシュポット、１３；スライダ、２１，４１；着目粒子、
２２，２３，３１，３２，４４；粒子、２４；小領域、２５，４３；領域、
４２；セル。
【図１】
2005122354.tif 000003
【図２】
2005122354.tif 000004
【図３】
2005122354.tif 000005
【図４】
2005122354.tif 000006
【図５】
2005122354.tif 000007
【図６】
2005122354.tif 000008
【図７】
2005122354.tif 000009
【図８】
2005122354.tif 000010
【図９】
2005122354.tif 000011
【図１０】
2005122354.tif 000012
【図１１】
2005122354.tif 000013
【図１２】
2005122354.tif 000014
【図１３】
2005122354.tif 000015
【図１４】
2005122354.tif 000016
【図１５】
2005122354.tif 000017
【図１６】
2005122354.tif 000018
【図１７】
2005122354.tif 000019
【図１８】
2005122354.tif 000020
【図１９】
2005122354.tif 000021