(57)【要約】【課題】  少なくとも粒子同士の接触力を受ける粒子の挙動を計算する粒子挙動計算方法において、その計算を効率的にしかも正確に実行可能とすること。【解決手段】  （Ａ）に示すような対象領域を（Ｂ）に示すようにセルに分割し、セル内の粒子数が１以下となるまで分割を繰り返してツリーを作成する。このとき、セルの幅が規定値以下となったら分割を中止する、ツリーを利用して接触力の計算対象となる粒子を探索する、隣接するセルに含まれる粒子のみを接触力の計算対象とする、などの方法によって、接触力等の計算の効率化を図る。【選択図】  図５
2008134959.tif 000002
【特許請求の範囲】
【請求項１】
対象領域を小領域に順次分割し、分割後の上記小領域内の粒子数が１以下となるか分割後の上記小領域の各辺の長さが第１の距離よりも小さくなるかするまで、上記分割を繰り返す領域分割処理と、
上記対象領域及び上記各分割時における各小領域に対応する節点を有する木構造を作成する木構造作成処理と、
該木構造作成処理により作成された木構造に基づき、注目粒子から第２の距離以内にある粒子を探索する粒子探索処理と、
該粒子探索処理により探索された粒子と上記注目粒子との間に作用する接触力をそれぞれ計算する接触力計算処理と、
該接触力計算処理により計算された接触力に基づき、上記注目粒子の挙動を計算する挙動計算処理と、
を備えたことを特徴とする粒子挙動計算方法。
【請求項２】
対象領域を小領域に順次分割し、分割後の上記小領域内の粒子数が１以下となるまで上記分割を繰り返す領域分割処理と、
上記対象領域及び上記各分割時における各小領域に対応する節点を有し、かつ、少なくとも最下位の上記節点にはその節点に対応する上記小領域内の粒子の位置を表すデータを備えた木構造を作成する木構造作成処理と、
該木構造作成処理により作成された木構造に基づき、注目粒子が属する上記小領域に隣接する上記小領域を探索する小領域探索処理と、
該小領域探索処理により探索された小領域内の粒子と上記注目粒子との間に作用する接触力をそれぞれ計算する接触力計算処理と、
該接触力計算処理により計算された接触力に基づき、上記注目粒子の挙動を計算する挙動計算処理と、
を備えたことを特徴とする粒子挙動計算方法。
【請求項３】
少なくとも最下位の上記各節点が、その上位の上記節点に関する節点情報を備え、上記小領域探索処理では、その節点情報を用いて上記隣接する小領域が探索されることを特徴とする請求項２記載の粒子挙動計算方法。
【請求項４】
上記注目粒子から、第３の距離以内にある上記小領域に対しては、その小領域内の粒子毎に、上記注目粒子との間に作用する相互作用力を計算し、上記第３の距離より遠くにある上記小領域に対しては、上記第３の距離より遠くにある最上位の上記小領域内の全粒子を１つの粒子とみなして上記注目粒子との間に作用する相互作用力を計算する相互作用力計算処理を、
更に備え、
上記挙動計算処理では、上記相互作用力計算処理にて計算された相互作用力と上記接触力とに基づいて、上記注目粒子の挙動が計算されることを特徴とする請求項１〜３のいずれかに記載の粒子挙動計算方法。
【請求項５】
上記各節点が、少なくともその節点に対応する上記小領域内の粒子の重量の総和量または帯電量の総和量と上記小領域の重心位置とを表すデータを備えたことを特徴とする請求項４記載の粒子挙動計算方法。
【請求項６】
上記各節点が、その節点に対応する上記小領域内の粒子数及び粒子の位置を表すデータを備えたことを特徴とする請求項１〜５のいずれかに記載の粒子挙動計算方法。
【請求項７】
請求項１〜６のいずれかに記載の上記各処理を、コンピュータに実行させるためのプログラム。
【請求項８】
請求項７記載のプログラムが、コンピュータによって読み取り可能に記録されたことを特徴とする記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、少なくとも粒子同士の接触力を受ける粒子の挙動を計算する粒子挙動計算方法、並びに、その方法を実施するためのプログラム及び記録媒体に関する。
【背景技術】
【０００２】
従来より、トナーなどのように各種の力を受ける粒子の挙動を計算する粒子挙動計算方法が種々提案されている。この種の粒子挙動計算方法では、粒子同士の間に作用する力を全て計算すると、粒子数Ｎの２乗の組み合せに対して計算を行う必要があり、膨大な時間を要してしまう。そこで、対象領域に一定幅のメッシュを生成し、注目粒子の存在するセルまたはそのセルから所定距離内に存在するセルに含まれる粒子との相互作用のみを計算することが提案されている（例えば、特許文献１参照）。
【特許文献１】特開２００５−１２２３５４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００３】
ところが、対象領域に一定幅のメッシュを生成した場合、粒子が存在しない領域にもセルが形成されて膨大なセル数となるため、各粒子をセルに割り当てる計算の負荷が極めて大きくなる。また、粒子に作用する力を正確に計算するためには上記所定距離を大きくする必要があり、多数の粒子同士の間に作用する力を個々に計算する必要が生じて充分に計算を効率化することができない。
【０００４】
そこで、本発明は、少なくとも粒子同士の接触力を受ける粒子の挙動を計算する粒子挙動計算方法において、その計算を効率的にしかも正確に実行可能とすることを目的としてなされた。
【課題を解決するための手段】
【０００５】
上記目的を達するためになされた本発明の粒子挙動計算方法は、対象領域を小領域に順次分割し、分割後の上記小領域内の粒子数が１以下となるか分割後の上記小領域の各辺の長さが第１の距離よりも小さくなるかするまで、上記分割を繰り返す領域分割処理と、上記対象領域及び上記各分割時における各小領域に対応する節点を有する木構造を作成する木構造作成処理と、該木構造作成処理により作成された木構造に基づき、注目粒子から第２の距離以内にある粒子を探索する粒子探索処理と、該粒子探索処理により探索された粒子と上記注目粒子との間に作用する接触力をそれぞれ計算する接触力計算処理と、該接触力計算処理により計算された接触力に基づき、上記注目粒子の挙動を計算する挙動計算処理と、を備えたことを特徴としている。
【０００６】
このように構成された本発明の粒子挙動計算方法では、領域分割処理により、対象領域が小領域に順次分割され、分割後の上記小領域内の粒子数が１以下となるか分割後の上記小領域の各辺の長さが第１の距離よりも小さくなるかするまで、上記分割が繰り返される。すると、木構造作成処理により、上記対象領域及び上記各分割時における各小領域に対応する節点を有する木構造が作成される。このため、粒子の存在密度に応じて対象領域が分割され、極めて扱い易い木構造が作成される。
【０００７】
そして、このように作成された木構造に基づき、注目粒子から第２の距離以内にある粒子が粒子探索処理によって探索され、その探索された粒子と上記注目粒子との間に作用する接触力が接触力計算処理によって計算される。続いて、挙動計算処理では、接触力計算処理により計算された接触力に基づき、上記注目粒子の挙動が計算される。
【０００８】
このように、本発明では、対象領域の分割に対応した木構造に基づいて接触力の計算対象の粒子が探索されるので、接触力を効率的に計算することができる。しかも、上記対象領域の分割は、分割後の小領域の各辺の長さが第１の距離よりも小さくなるとそれ以上の分割がなされず、対象領域が必要以上に細かく分割されないので、上記計算を一層効率的にしかも正確に実行することができる。従って、本発明では、粒子の挙動を効率的にしかも正確に計算することができる。
【０００９】
また、上記目的を達するためになされた本発明の粒子挙動計算方法は、対象領域を小領域に順次分割し、分割後の上記小領域内の粒子数が１以下となるまで上記分割を繰り返す領域分割処理と、上記対象領域及び上記各分割時における各小領域に対応する節点を有し、かつ、少なくとも最下位の上記節点にはその節点に対応する上記小領域内の粒子の位置を表すデータを備えた木構造を作成する木構造作成処理と、該木構造作成処理により作成された木構造に基づき、注目粒子が属する上記小領域に隣接する上記小領域を探索する小領域探索処理と、該小領域探索処理により探索された小領域内の粒子と上記注目粒子との間に作用する接触力をそれぞれ計算する接触力計算処理と、該接触力計算処理により計算された接触力に基づき、上記注目粒子の挙動を計算する挙動計算処理と、を備えたことを特徴とするものであってもよい。
【００１０】
このように構成された本発明の粒子挙動計算方法では、領域分割処理により、対象領域が小領域に順次分割され、分割後の上記小領域内の粒子数が１以下となるまで、上記分割が繰り返される。すると、木構造作成処理により、上記対象領域及び上記各分割時における各小領域に対応する節点を有し、かつ、少なくとも最下位の上記節点にはその節点に対応する上記小領域内の粒子の位置を表すデータを備えた木構造が作成される。このため、粒子の存在密度に応じて対象領域が分割され、極めて扱い易い木構造が作成される。
【００１１】
そして、このように作成された木構造に基づき、注目粒子が属する小領域に隣接する小領域が小領域探索処理によって探索され、その探索された小領域内の粒子と上記注目粒子との間に作用する接触力が接触力計算処理によってそれぞれ計算される。続いて、挙動計算処理では、接触力計算処理により計算された接触力に基づき、上記注目粒子の挙動が計算される。
【００１２】
このように、本発明では、対象領域の分割に対応した木構造に基づいて、注目粒子が属する小領域に隣接する小領域が探索され、その小領域内の粒子が接触力の計算対象とされる。このため、本発明では、粒子間の接触力を効率的にしかも正確に計算することができ、延いては、粒子の挙動を効率的にしかも正確に計算することができる。
【００１３】
なお、上記各発明の粒子挙動計算方法は以下の構成に限定されるものではないが、後者の粒子挙動計算方法において、少なくとも最下位の上記各節点が、その上位の上記節点に関する節点情報を備え、上記小領域探索処理では、その節点情報を用いて上記隣接する小領域が探索されてもよい。
【００１４】
この場合、少なくとも最下位の節点が、その上位の節点に関する節点情報を備えているので、注目粒子が属する小領域から１つ上位の小領域へ遡って、その小領域の下位の小領域を容易に探索することができる。このような小領域は、注目粒子が属する小領域に隣接する小領域である。従って、この場合、小領域探索処理を一層効率化して、延いては、粒子の挙動を一層効率的にしかも正確に計算することができる。
【００１５】
また、上記各発明の粒子挙動計算方法において、上記注目粒子から、第３の距離以内にある上記小領域に対しては、その小領域内の粒子毎に、上記注目粒子との間に作用する相互作用力を計算し、上記第３の距離より遠くにある上記小領域に対しては、上記第３の距離より遠くにある最上位の上記小領域内の全粒子を１つの粒子とみなして上記注目粒子との間に作用する相互作用力を計算する相互作用力計算処理を、更に備え、上記挙動計算処理では、上記相互作用力計算処理にて計算された相互作用力と上記接触力とに基づいて、上記注目粒子の挙動が計算されてもよい。
【００１６】
この場合、接触力のみならず、注目粒子との間に作用する相互作用力にも基づいて注目粒子の挙動が計算されるので、注目粒子の挙動を一層正確に計算することができる。しかも、その相互作用力を計算する相互作用力計算処理では、第３の距離以内にある小領域に対しては、その小領域内の粒子毎に相互作用力が計算されるが、第３の距離より遠くにある小領域に対しては、その第３の距離より遠くにある最上位の小領域内の全粒子を１つの粒子とみなして相互作用力が計算される。従って、この場合、相互作用力の計算を一層効率的に実行することができ、延いては、粒子の挙動を一層効率的にしかも正確に計算することができる。
【００１７】
また、この場合、上記各節点が、少なくともその節点に対応する上記小領域内の粒子の重量の総和量または帯電量の総和量と上記小領域の重心位置とを表すデータを備えていてもよい。各節点がこのようなデータを備えている場合、第３の距離より遠くにある小領域内の全粒子を１つの粒子とみなして相互作用力を計算する際に、そのデータを利用することができる。従って、この場合、粒子の挙動を一層効率的に計算することができる。
【００１８】
また、上記各発明の粒子挙動計算方法において、上記各節点が、その節点に対応する上記小領域内の粒子数及び粒子の位置を表すデータを備えていてもよい。各節点がこのようなデータを備えている場合、そのデータを利用することにより、粒子の探索等を一層効率的に実行することができる。
【００１９】
また、本発明のプログラムは、上記いずれかの発明の粒子挙動計算方法の各処理を、コンピュータに実行させることを特徴としている。このため、本発明のプログラムをコンピュータに実行させれば、上記いずれかの発明の粒子挙動計算方法を容易に実行することができる。
【００２０】
また、本発明の記録媒体は、上記発明のプログラムが、コンピュータによって読み取り可能に記録されたことを特徴としている。このため、本発明の記録媒体に記録されたプログラムをコンピュータに実行させれば、上記いずれかの発明の粒子挙動計算方法を容易に実行することができる。
【発明を実施するための最良の形態】
【００２１】
［粒子挙動計算装置の構成］
以下に、本発明の実施の形態について、図面を参照して説明する。図１は、本発明が適用された粒子挙動計算方法を実行可能な粒子挙動計算装置としてのパーソナルコンピュータ１の構成を概略的に表す説明図である。本実施の形態の粒子挙動計算方法を実行するためのプログラムは、このパーソナルコンピュータ１のハードディスク装置（ＨＤＤ）３４に格納されている。パーソナルコンピュータ１は、この記録媒体の一例としてのハードディスク装置３４の他、ＣＰＵ３１，ＲＯＭ３２，ＲＡＭ３３を備えた本体３に、ディスプレイ５、キーボード７、マウス９等を接続した一般的な構成を有している。
【００２２】
［参考例としての粒子挙動計算方法］
先ず、本発明の実施の形態の粒子挙動計算方法を説明するに先立ち、その粒子挙動計算方法を理解するための参考となる参考例を説明する。本参考例では、パーソナルコンピュータ１に粒子挙動の計算が指示されると、ＣＰＵ３１はハードディスク装置３４に記録されたプログラムに基づき、次のような粒子挙動計算処理を実行する。
【００２３】
図２は、本参考例の粒子挙動計算処理のメインルーチンを表すフローチャートである。図２に示すように、この処理では、先ず、Ｓ１０（Ｓはステップを表す：以下同様）にて、キーボード７，マウス９等からの入力に基づいて粒子の初期配置が設定される。続いて、後述するツリー作成処理（Ｓ２０），力計算処理（Ｓ３０），位置計算処理（Ｓ４０）が順次実行され、Ｓ５０にて粒子挙動計算の終了時間であるか否かが判断される。終了時間でない場合は（Ｓ５０：Ｎ）、次のタイミングの粒子挙動を計算すべく、Ｓ５１にて時間が進められ、処理は前述のＳ２０へ移行する。こうして、時間を進めながらＳ２０〜Ｓ４０の処理がくり返し実行され、終了時間となると（Ｓ５０：Ｙ）、処理が終了する。
【００２４】
図３は、Ｓ２０のツリー作成処理を詳細に表すフローチャートである。図３に示すように、この処理では、先ず、Ｓ２１にて、対象領域全体が頂上セルとして設定され、続くＳ２２にて、下位の木作成処理が実行される。
【００２５】
図４は、この下位の木作成処理を詳細に表すフローチャートである。図４に示すように、この処理では、先ず、Ｓ２２０にて、注目セル（最初は頂上セル）に含まれる粒子数が１以下であるか否かが判断される。通常、頂上セル、すなわち対象領域全体には、複数の粒子が配置されているので、最初はＳ２２０にて否定判断され、処理はＳ２２１へ移行する。
【００２６】
Ｓ２２１では、注目セル内の粒子を小領域の一例としての子セル、すなわち１つ下位のセルに振り分ける処理がなされる。例えば、０≦ｘ≦１０，０≦ｙ≦１０で表れる平面状の対象領域に図５（Ａ）に例示するように粒子が配置された場合、図５（Ｂ）に太線で示すように４分割（対象領域が立体であれば８分割）された子セルが設定され、各子セルに粒子が振り分けられる。なお、各セルには、頂上セルを１として順次連続番号が設定される。
【００２７】
続くＳ２２２では、Ｓ２２１で作成された各子セルに対してその子セルのパラメータが算出される。ここで、パラメータとは、例えば、セル内の粒子数及び粒子の位置（座標）、上位のセルの番号（節点情報に相当）、セル内の粒子の重量の総和量、セル内の粒子の帯電量の総和量、セルの重心位置などである。
【００２８】
続くＳ２２３では、Ｓ２２２で設定された各子セルを注目セルとして、本図に示す下位の木作成処理が実行される。ここで、子セルの数は前述のように予め決められている。このため、Ｓ２２３は子セルの数だけ繰返し実行される。そこで、Ｓ２２３に続くＳ２２４では、Ｓ２２３の処理が子セルの数に対応する所定回繰り返されたか否かが判断され、まだの場合は（Ｓ２２４：Ｎ）、再びＳ２２３が実行される。そして、子セルの数だけＳ２２３が実行されると（Ｓ２２４：Ｙ）、処理はＳ３０の力計算処理（図２参照）へ移行する。
【００２９】
また、頂上セルよりも下位のセルを注目セルとした下位の木作成処理では、粒子数が１以下である場合が生じるが、この場合は（Ｓ２２０：Ｙ）、そのまま処理を終了して１つ上位の下位の木作成処理におけるＳ２２４へ移行する。頂上セルよりも下位のセルを注目セルとした下位の木作成処理でＳ２２３が、その注目セルの子セルの数だけ繰り返された場合も（Ｓ２２４：Ｙ）、同様に１つ上位の下位の木作成処理におけるＳ２２４へ移行する。
【００３０】
以上の処理により、Ｓ２０のツリー作成処理では、図５（Ｂ）に例示するように対象領域が順次分割され、その分割は各セルの粒子数が１以下となるまで繰り返される。そして、各セルには、頂上セルを「１」として順次増加する連続番号が、図５（Ｂ）に四角で囲んで示すように設定される。なお、図５（Ｂ）の例では、粒子数が０のセルには番号を付与していない。
【００３１】
このようなＳ２０のツリー作成処理によって、図６に例示するような木構造の一例としてのツリーが作成される。また、このツリーの各節点は、各セルに対応する番号と、上記各セルのパラメータとをデータとして備えている。なお、図６では、粒子数が０のセルに対応する節点（実際には節点としては形成されない）を丸にバツ印で、粒子数が１のセルに対応する節点を黒丸で、粒子数が２以上のセルに対応する節点を中間色の丸で、それぞれ表している。また、図６に例示するようなツリーの各枝（すなわち節点のデータ同士の対応関係）は、前述のように下位の木作成処理が１つ上位の下位の木作成処理に繰り上がる際に作成される。
【００３２】
次に、Ｓ３０の力計算処理では、このように作成されたツリーを用いて次のように各粒子に作用する力が計算される。図７はこの力計算処理を詳細に表すフローチャートである。図７に示すように、この処理では、上記各粒子のうちの１つを注目粒子として、その注目粒子に作用するクーロン力等の相互作用力を計算する相互作用力計算処理（Ｓ３１）、及び、その注目粒子に作用する粒子間の接触力を計算する接触力計算処理（Ｓ３２）が順次実行される。続くＳ３３では、それらの計算処理が全粒子に対して終了したか否かが判断され、終了していない場合は（Ｓ３３：Ｎ）、他の粒子を注目粒子として同様の相互作用力計算処理（Ｓ３１），接触力計算処理（Ｓ３２）が繰返し実行される。そして、全粒子に対して計算が終了すると（Ｓ３３：Ｙ）、処理はＳ４０の位置計算処理（図２参照）へ移行する。この位置計算処理（Ｓ４０）は、Ｓ３０にて計算された各粒子に作用する力を時間で２回積分することにより、Ｓ５１にて進められた時間における各粒子の位置を計算する周知の処理である。
【００３３】
次に、Ｓ３１の相互作用力計算処理、及び、Ｓ３２の接触力計算処理について、更に詳細に説明する。図８は、Ｓ３１の相互作用力計算処理を詳細に表すフローチャートである。図８に示すように、この処理では、先ず、Ｓ３１０にて、計算対象セルの重心と注目粒子との距離Ｄが計算される。なお、本処理は、頂上セルを計算対象セルとして開始される。
【００３４】
続くＳ３１１では、その計算対象セルが注目粒子から充分に遠いか否かが判断される。なお、この判断では、Ｓ３１０にて計算された距離Ｄが、計算対象セルの一辺の長さＬと以下の関係を満たす場合に、充分に遠いと判断される。後述のように、次式においてＬ／θが第３の距離の一例に相当する。
【００３５】
Ｌ／Ｄ＜θ      （但し、θは必要精度や系に応じて適宜設定されるパラメータ）
頂上セルは一辺の長さＬが極めて大きいため、一般的に最初はＳ３１１で否定判断されて、処理はＳ３１２へ移行する。Ｓ３１２では、計算対象セル内の粒子数が１個であるか否かが判断されるが、ここでも最初はＳ２２０と同様に否定判断されて、処理はＳ３１３へ移行する。Ｓ３１３では、計算対象セルの各子セルを新たな計算対象セルとして、本図に示す相互作用力計算処理がそれぞれ実行される。そして、各子セルに対する相互作用力計算処理がそれぞれ終了すると、処理は後述のＳ３１６へ移行する。
【００３６】
頂上セルよりも下位のセルを計算対象セルとした相互作用力計算処理では、その計算対象セルが注目粒子から充分遠いと判断される場合や（Ｓ３１１：Ｙ）、その計算対象セル内の粒子数が１個である場合があり（Ｓ３１２：Ｙ）、これらの場合は処理はＳ３１４へ移行する。
【００３７】
Ｓ３１４では、そのセルに対してＳ２２２で計算されているパラメータから、注目粒子に作用する相互作用力が計算される。例えば、相互作用力としてクーロン力のみを計算する場合、注目粒子の位置ベクトル及び帯電量をｘi，ｑi、計算対象セルの重心の位置ベクトル及び帯電量の総和をｘj，ｑjとすると、注目粒子に作用するクーロン力Ｆi は次式で表される。
【００３８】
【数１】
2008134959.tif 000003
上記処理では、例えば、図９に例示するように、計算対象セルＣの重心Ｇと注目粒子Ｐ０との距離Ｄが、計算対象セルＣの一辺の長さＬに対して充分に大きい場合は（Ｓ３１１：Ｙ）、計算対象セルＣに含まれる全ての粒子を１つの粒子とみなしてＳ３１４により相互作用力を計算することができる。また、計算対象セルの粒子数が１個の場合は（Ｓ３１２：Ｙ）、セルの重心の位置及び帯電量の総和は、そのセルに含まれる１個の粒子の位置及び帯電量と一致する。従って、この場合のＳ３１４の処理は、１対１の粒子間に作用する相互作用力を計算したのと同様の結果となる。
【００３９】
続くＳ３１５では、Ｓ３１４にて計算された相互作用力が、注目粒子に対して計算された他の相互作用力に加算され、Ｓ３１６にて、全ての子セルに対して上記計算が終了したか否かが判断される。そして、まだの場合は（Ｓ３１６：Ｎ）、未計算の子セルに対してＳ３１０〜Ｓ３１５の処理が実行される。一方、全ての子セルに対して上記計算が終了すると（Ｓ３１６：Ｙ）、処理は１つ上位のセルに対する相互作用計算処理のＳ３１６、または、Ｓ３２の接触力計算処理（図７参照）へ移行する。
【００４０】
次に、図１０は、Ｓ３２の接触力計算処理を詳細に表すフローチャートである。この処理では、注目粒子以外の粒子が１つ選択され、Ｓ３２０にてその粒子と注目粒子との間に作用する接触力が計算される。この接触力の計算は、一般に知られた計算方法で行えばよく、例えば、粉体工学会編「粉体シミュレーション入門」産業図書、１９９８年４月、などに詳しく開示されている。続くＳ３２１では、Ｓ３２０にて計算された接触力がその注目粒子に対して計算された他の接触力に加算され、Ｓ３２２にて、注目粒子に作用する接触力が全ての粒子に対して計算されたか否かが判断される。そして、まだの場合は（Ｓ３２２：Ｎ）、未計算の粒子に対してＳ３２０，Ｓ３２１の処理が実行され、注目粒子に作用する接触力が全ての粒子に対して計算されると（Ｓ３２２：Ｙ）、処理は前述のＳ３３へ移行する（図７参照）。
【００４１】
以上の処理により、各粒子に作用する接触力及び相互作用力に応じた、各粒子の詳細な挙動を計算することができる。また、上記処理では、粒子の存在密度に応じて対象領域が分割され、極めて扱い易いツリーを作成することができる。しかも、上記処理では、注目粒子から充分遠くにあるセルに対しては、そのセル内の全粒子を１つの粒子とみなして相互作用力が計算されるので、相互作用力の計算を一層効率的に実行することができる。また、前述のように予め計算されたセルのパラメータを利用することにより、このような計算も極めて効率的に実行することができる。
【００４２】
［第１の実施の形態としての粒子挙動計算方法］
次に、第１の実施の形態の粒子挙動計算方法を説明する。先ず、本実施の形態は、上記参考例に次のような改良を加えた方法である。すなわち、図１１に例示するように、粒子中の最大粒子径Ｒｍａｘの２倍よりもセルの一辺の長さＬが大きければ、全ての接触粒子を含むことができる。そこで、本実施の形態では、セルの一辺の長さＬ（以下、セルの幅ともいう）が第１の距離の一例としての規定値以下となった時点でセルの分割を中止し、接触力の計算においても所定の目標距離（第２の距離の一例で、具体的には例えば２×Ｒｍａｘ）より近い粒子のみを計算対象とする改良を行った。
【００４３】
以下、本実施の形態の上記参考例に対する改良点について説明する。本実施の形態では、メインルーチン及びツリー作成処理は図２，図３に示したものと同様で、ツリー作成処理のうちの下位の木作成処理（Ｓ２２）が次のように改良されている。すなわち、図１２に示すように、Ｓ２２０にて注目セルに含まれる粒子数が複数と判断された場合（Ｓ２２０：Ｎ）、新たに追加されたＳ２２０Ａにて、注目セルの幅が規定値より大きいか否かが判断される。そして、注目セルの幅が規定値よりも大きい場合は（Ｓ２２０Ａ：Ｙ）、前述のＳ２２１〜Ｓ２２４の処理が実行される。一方、注目セルの幅が規定値以下の場合は（Ｓ２２０Ａ：Ｎ）、注目セルに含まれる粒子数が１以下である場合と同様に（Ｓ２２０：Ｙ）、そのまま処理を終了して１つ上位の下位の木作成処理におけるＳ２２４、またはＳ３０の力計算処理（図２参照）へ移行する。
【００４４】
このように改良されたツリー作成処理では、セルの粒子数が１以下となるか（Ｓ２２０：Ｙ）、若しくはセルの幅（すなわち一辺の長さＬ）が規定値以下となるかするまで（Ｓ２２０Ａ：Ｙ）、セルの分割がなされる。このため、必要以上にセルを細かく分割することが抑制される。
【００４５】
図１３は、本実施の形態における力計算処理を表すフローチャートである。図１３に示すように、本実施の形態の力計算処理は、参考例の相互作用力計算処理とほぼ同様に構成されている。すなわち、頂上セルを計算対象セルとして処理が開始され、計算対象セルの重心と注目粒子との距離Ｄが計算される（Ｓ３１０）。なお、計算対象セルの重心ではなく、計算対象セルの中心と注目粒子との距離を計算しても差し支えない。その距離ＤがＬ／Ｄ＜θを満たす充分に遠い距離であるか否かが判断される点は（Ｓ３１１）、参考例の相互作用力計算処理と同様である。また、充分遠い場合（Ｓ３１１：Ｙ）、計算対象セルに含まれる全粒子を１つの粒子とみなして、セルのパラメータから相互作用力が計算され（Ｓ３１４）、既に計算されている相互作用力に加算されて（Ｓ３１５）、全ての子セルに計算が終了していなければ処理がＳ３１０へ戻る点（Ｓ３１６）も同様である。
【００４６】
一方、本実施の形態の力計算処理では、距離Ｄが充分に遠くない場合は（Ｓ３１１：Ｎ）、Ｓ３１７にて、計算対象セルに子セルがないか否かが判断される。そして、計算対象セルに子セルがある場合は（Ｓ３１７：Ｎ）、処理はＳ３１８へ移行し、計算対象セルの各子セルを新たな計算対象セルとして、本図に示す力計算処理がそれぞれ実行される。そして、各子セルに対する力計算処理がそれぞれ終了すると、処理はＳ３１６へ移行する。また、計算対象セルに子セルがない場合は（Ｓ３１７：Ｙ）、次に説明する粒子間力計算処理がＳ３１９にて実行された後、処理はＳ３１６へ移行する。こうして、全ての子セルに対してＳ３１０〜Ｓ３１９の処理が終了すると（Ｓ３１６：Ｙ）、処理は１つ上位のセルに対する力計算処理のＳ３１６、またはＳ４０の位置計算処理（図２参照）へ移行する。
【００４７】
図１４は、Ｓ３１９の粒子間力計算処理を詳細に表すフローチャートである。この処理では、その時点における計算対象セルに含まれる粒子が上記ツリーに基づいて１つ選択される。すなわち、計算対象セルに対応する節点にぶら下がる子セルの１つがツリーに基づいて選択される。そして、Ｓ３１９０にて、その粒子が注目粒子から上記目標距離より近いか否かが各粒子の位置（座標）のデータを用いて判断される。上記目標距離は、最大粒子径の２倍に設定されている。なお、目標距離は、実際の対象粒子の半径の和を求めてその値に設定しても差し支えない。目標距離より近い場合は（Ｓ３１９０：Ｙ）、Ｓ３１９１にて前述のＳ３２０と同様に接触力が計算された後、目標距離以遠の場合はそのまま（Ｓ３１９０：Ｎ）、処理はＳ３１９２へ移行する。
【００４８】
Ｓ３１９２では、選択された粒子から注目粒子に作用するクーロン力等の相互作用力が計算され、Ｓ３１９３にて、既に計算されている接触力，相互作用力に上記計算された接触力，相互作用力が加算される。続くＳ３１９４では、計算対象セルに含まれる全ての粒子に対して上記計算が終了したか否かが判断される。すなわち、上記計算対象セルに対応する節点にぶら下がる子セルについて全て上記計算が終了したか否かが判断される。そして、まだの場合は（Ｓ３１９４：Ｎ）、未計算の粒子に対してＳ３１９０〜３１９３の処理が実行される。そして、計算対象セルに含まれる全ての粒子に対して上記計算が終了すると（Ｓ３１９４：Ｙ）、処理は前述のＳ３１６へ移行する（図１３参照）。
【００４９】
本実施の形態では、上記参考例の効果に加えて、更に、対象領域の分割に対応したツリーに基づいて接触力の計算対象の粒子が選択可能であるので、接触力を効率的に計算することができる。しかも、セルの分割が必要以上に細かくなされないので、上記計算を一層効率的にしかも正確に実行することができる。従って、本実施の形態では、粒子の挙動を一層効率的にしかも正確に計算することができる。
【００５０】
なお、上記実施の形態においてＳ２２１の処理のうち領域を分割する処理が領域分割処理に、Ｓ２０のツリー作成処理における他の処理が木構造作成処理に、Ｓ３１９０の処理が粒子探索処理に、Ｓ３１９１の処理が接触力計算処理に、Ｓ４０の処理が挙動計算処理に、Ｓ３１４，Ｓ３１９２の処理が相互作用力計算処理に、それぞれ相当する。
【００５１】
［第２の実施の形態としての粒子挙動計算方法］
次に、第２の実施の形態としての粒子挙動計算方法を説明する。本実施の形態は、上記参考例に次のような改良を加えた方法である。すなわち、図１５（Ａ）に例示するように、注目粒子Ｐ０がセルＣ０に含まれている場合、接触力の計算は、そのセルＣ０に隣接するセルＣ１〜Ｃ６に含まれる粒子に対して行うだけでも充分な場合がある。そこで、本実施の形態では、以下に説明するように隣接するセル等の探索を行っている。
【００５２】
本実施の形態では、メインルーチン及びツリー作成処理は、下位の木作成処理も含めて図２〜図４に示したものと同様で、力計算処理（Ｓ３０）が次のように改良されている。図１６は、本実施の形態の力計算処理を表すフローチャートである。図１６に示すように、この処理では、先ずＳ３１にて、参考例と同様の相互作用力計算処理（図８参照）が実行され、続くＳ３５にて、接触セルリスト作成処理が実行される。
【００５３】
図１７は、Ｓ３５の接触セルリスト作成処理を詳細に表すフローチャートである。図１７に示すように、この処理では、先ずＳ３５０にて、注目粒子を含むセルの親セルの子セル、すなわち、１つ上位のセルの１つ下位のセルがリストに追加される。図１５（Ａ）のセルＣ０に対してはセルＣ１，Ｃ２，Ｃ３がこのようなセルに相当し、このような親セルの子セルは注目粒子を含むセルに必ず隣接している。また、このようなセルの探索は、上記計算されたセルのパラメータを利用して実行される。
【００５４】
続くＳ３５１では、上記親セルの子セル以外の、注目粒子を含むセルに近いセルの探索が実行される。図１８は、この処理の詳細を表すフローチャートである。図１８に示すように、この処理では、先ず、Ｓ３５１０にて、注目粒子の座標に基づき、近いセルとして探索すべきセルの座標範囲が計算される。なお、この処理では、注目粒子を含むセルの一辺の長さＬが粒子径に比べて充分に大きい場合は、図１５（Ａ）に網掛けをして例示するように隣接するセルの座標範囲を探索セル座標として計算すればよいが、図１５（Ａ）に円で例示するように、上記セルの一辺の長さＬに比べて粒子径が大きい場合は、図１５（Ｂ）に網掛けをして例示するように、上位のセルに遡って隣接する探索セル座標が計算される。
【００５５】
続くＳ３５１１では、ツリーのデータを用いて上記範囲のセルが探索され、Ｓ３５１２にて、当該探索されたセルがリストに追加されて処理は図１７のＳ３５２へ移行する。Ｓ３５２では、Ｓ３５０またはＳ３５１にて探索されたセルが接触セルリストに追加され、処理は図１６のＳ３６へ移行する。
【００５６】
図１９は、Ｓ３６の接触力計算処理を詳細に表すフローチャートである。この処理では、接触セルリストに登録されたセルが１つ選択され、更に、そのセルに含まれる粒子が１つ選択された上で、Ｓ３６０にて前述のＳ３２０と同様に接触力が計算される。続くＳ３６１では、Ｓ３６０にて計算された接触力が既に計算された他の接触力に加算され、Ｓ３６２にて、当該セルに含まれる全粒子に対して計算が終了したか否かが判断される。
【００５７】
まだの場合は（Ｓ３６２：Ｎ）、未計算の粒子に対してＳ３６０，Ｓ３６１の処理が実行され、当該セルに含まれる全粒子に対して上記計算が終了すると（Ｓ３６２：Ｙ）、処理はＳ３６３へ移行する。Ｓ３６３では、接触セルリストに登録された全てのセルに対して上記計算が終了したか否かが判断され、まだの場合は（Ｓ３６３：Ｎ）、残りのセルに対しても同様にＳ３６０〜Ｓ３６２の処理が実行される。そして、接触セルリストに登録された全てのセルに対して上記計算が終了すると（Ｓ３６３：Ｙ）、処理は参考例と同様のＳ３３（図１６参照）へ移行し、上記Ｓ３１〜Ｓ３６の処理が他の粒子を注目粒子として繰り返される。そして、全粒子に対して計算が終了すると（Ｓ３３：Ｙ）、処理はＳ４０の位置計算処理（図２参照）へ移行する。
【００５８】
本実施の形態では、上記参考例の効果に加えて、更に、対象領域の分割に対応したツリーに基づいて注目粒子が含まれるセルに隣接または近接するセルが探索され、そのセルに含まれる粒子のみが接触力の計算対象とされるので、接触力を効率的に計算することができる。しかも、上記セルの探索は、セルのパラメータを用いてなされるので、その探索を一層効率的に実行することができる。従って、本実施の形態では、粒子の挙動を一層効率的にしかも正確に計算することができる。
【００５９】
なお、上記実施の形態においてＳ２２１の処理のうち領域を分割する処理が領域分割処理に、Ｓ２０のツリー作成処理における他の処理が木構造作成処理に、Ｓ３５の処理が小領域探索処理に、Ｓ３６の処理が接触力計算処理に、Ｓ４０の処理が挙動計算処理に、Ｓ３１の処理が相互作用力計算処理に、それぞれ相当する。
【００６０】
次に、図２０は、第２の実施の形態（改良ツリー法）によって粒子挙動を計算した場合と、各粒子間に作用する力を総当りで計算した場合との計算時間を比較したグラフである。なお、本計算は、レーザプリンタにおけるトナーカートリッジ中のトナー粒子の挙動を表す系を想定して行った。また、本計算は、簡易的に接触力の計算のみで比較したが、接触力及び相互作用力の両方を計算する場合、その計算時間の差はもっと顕著となる。
【００６１】
図２０に示すように、第２の実施の形態では、総当りに比べて極めて効率的に粒子挙動を計算することができ、その効果は粒子数が多くなるほど顕著になった。すなわち、粒子数をＮとした場合、総当りの場合ではＮ2 のオーダーで計算時間が増加するが、第２の実施の形態ではＮｌｏｇＮのオーダーまで計算時間を圧縮することができた。
【００６２】
なお、本発明は上記実施の形態に何ら限定されるものではなく、本発明の要旨を逸脱しない範囲で種々の形態で実施することができる。例えば、第２の実施の形態でも、下位の木作成処理を第１の実施の形態と同様にして、第１実施の形態と第２の実施の形態とを組み合せてもよい。この場合、Ｓ３５０では注目粒子が含まれるセル（図１５（Ａ）の例ではＣ０）自身も親セルの子セルとしてリストに追加され、Ｓ３５１では、注目粒子が含まれるセルに隣接しているセルのみが探索される。
【００６３】
また、相互作用力としては、万有引力，ファンデルワールス力など、クーロン力以外の力も計算してもよい。更に、本発明の記録媒体は、上記実施の形態のようなハードディスク装置、或いはＲＯＭ，ＲＡＭ等の素子に限定されるものではなく、フレキシブルディスク，コンパクトディスク，インターネット上のウェブサーバ等であってもよい。
【図面の簡単な説明】
【００６４】
【図１】本発明が適用された粒子挙動計算方法を実行可能なパーソナルコンピュータの構成を概略的に表す説明図である。
【図２】参考例の粒子挙動計算処理のメインルーチンを表すフローチャートである。
【図３】その粒子挙動計算処理のツリー作成処理を表すフローチャートである。
【図４】そのツリー作成処理の下位の木作成処理を表すフローチャートである。
【図５】その処理による対象領域の分割を例示する説明図である。
【図６】上記処理によって作成されるツリーを例示する説明図である。
【図７】上記粒子挙動計算処理の力計算処理を表すフローチャートである。
【図８】その力計算処理の相互作用力計算処理を表すフローチャートである。
【図９】その相互作用力計算処理の考え方を例示する説明図である。
【図１０】上記力計算処理の接触力計算処理を表すフローチャートである。
【図１１】第１の実施の形態の考え方を例示する説明図である。
【図１２】第１の実施の形態の下位の木作成処理を表すフローチャートである。
【図１３】第１の実施の形態の力計算処理を表すフローチャートである。
【図１４】その力計算処理の粒子間力計算処理を表すフローチャートである。
【図１５】第２の実施の形態の考え方を例示する説明図である。
【図１６】第２の実施の形態の力計算処理を表すフローチャートである。
【図１７】その力計算処理の接触セルリスト作成処理を表すフローチャートである。
【図１８】その接触セルリスト作成処理の近いセルの探索処理を表すフローチャートである。
【図１９】上記力計算処理の接触力計算処理を表すフローチャートである。
【図２０】第２の実施の形態の効果を表すグラフである。
【符号の説明】
【００６５】
１…パーソナルコンピュータ        ３…本体            ５…ディスプレイ
７…キーボード                    ３１…ＣＰＵ        ３４…ハードディスク装置
Ｃ，Ｃ０，Ｃ１，Ｃ２，Ｃ３，Ｃ４，Ｃ５，Ｃ６…セル    Ｐ０…注目粒子
【図１】
2008134959.tif 000004
【図２】
2008134959.tif 000005
【図３】
2008134959.tif 000006
【図４】
2008134959.tif 000007
【図７】
2008134959.tif 000008
【図８】
2008134959.tif 000009
【図１０】
2008134959.tif 000010
【図１２】
2008134959.tif 000011
【図１３】
2008134959.tif 000012
【図１４】
2008134959.tif 000013
【図１６】
2008134959.tif 000014
【図１７】
2008134959.tif 000015
【図１８】
2008134959.tif 000016
【図１９】
2008134959.tif 000017
【図５】
2008134959.tif 000018
【図６】
2008134959.tif 000019
【図９】
2008134959.tif 000020
【図１１】
2008134959.tif 000021
【図１５】
2008134959.tif 000022
【図２０】
2008134959.tif 000023