JPH07140059A - 粒子挙動解析装置 - Google Patents

Abstract

(57)【要約】 【目的】 精度を確保した上で演算の高速化を可能とす
る粒子挙動解析装置を提供する。 【構成】 初回の演算時には、粒子以外に起因する電磁
場に基づいて各粒子の双極子モーメントが求められ、以
降の演算時には、前回の処理時に求められた他の粒子に
起因する電磁場と、粒子以外に起因する電磁場とに基づ
いて各粒子の双極子モーメントが求められる。

Description

【発明の詳細な説明】

【０００１】

【産業上の利用分野】この発明は粒子挙動解析装置に関
し、特に磁場や電場における粒子に加わる力を計算して
その挙動を解析する装置に関するものである。

【０００２】

【従来の技術】電子写真の現像器における課題（穂切り
の最適化、トナー飛散対策、クリーニング等）を解決す
るために現像剤の挙動をシミュレートするソフトウエア
が開発されている。

【０００３】そのシミュレートにおいては、以下の１）
〜３）の計算が繰返されている。 １） 現在の粒子の位置の計算 ２） 粒子に加わる力の計算 ３） 上記２）の力による微小時間経過後の粒子の位置
の計算

【０００４】

【発明が解決しようとする課題】上記２）の計算におい
て、粒子に加わる力は ステップ−１：粒子位置の場の状態を計算する ステップ−２：上記で計算された場の状態によって定ま
る量（磁気双極子モーメント等）と場の状態とからその
粒子に加わる力を求めるという２つのステップによって
計算される。

【０００５】しかし、上記の方法で計算すると計算量が
膨大となり、実質的な計算は不可能となる。

【０００６】一方、ＰａｒａｎｊｐｅとＥｌｒｏｄによ
る論文（１９８６）において、上記２）での“粒子に加
わる力”として、 （粒子に加わる力）＝（粒子以外に起因する場による
力）＋（各粒子の作る場による力）という力の分解が示
されている。ただし、この論文の内容によると上記の計
算は数式的には可能であるが、これだけでは以下の理由
でソフトウエアは作成できない。

【０００７】１） 数式で表せない一般的な場の扱いは
この文献では述べられていない。 ２） この論文の内容そのままでは粒子の双極子モーメ
ントと粒子位置の磁束密度は決定できない。

【０００８】さらに、Ｊａｐａｎ Ｈａｒｄｃｏｐｙ
’９２（松下電産）において、計算ステップごとに粒
子（ここではトナー）を含めた電場計算を行なう内容が
開示されている。この方法は精度および汎用性ともに高
いが、各計算ステップ毎に有限要素計算を行なうため計
算量が極端に多くなってしまう。

【０００９】この発明は上記のような課題を解決するた
めになされたもので、精度を確保した上で、演算の高速
化を可能とする粒子挙動解析装置を提供することを目的
とする。

【００１０】

【課題を解決するための手段】この発明に係る粒子挙動
解析装置は、有限な領域内での電磁場の影響下における
複数の粒子の挙動を解析する粒子挙動解析装置におい
て、各粒子の初期状態および電磁場の状態を入力するた
めの入力手段と、各粒子の状態を記憶する記憶手段と、
任意の時間間隔での各粒子の挙動を演算し、記憶手段に
記憶された各粒子の状態を更新する粒子挙動演算手段
と、記憶手段に記憶された各粒子の状態を出力するため
の出力手段とを備え、粒子挙動演算手段が、任意の時間
間隔での各粒子の位置における粒子以外に起因する第１
の電磁場を、入力手段から入力された電磁場の状態に基
づいて求める第１電磁場処理手段と、任意の時間間隔で
の各粒子の双極子モーメントを求める双極子モーメント
処理手段と、任意の時間間隔での各粒子の位置における
他の粒子に起因する第２の電磁場を、各粒子の双極子モ
ーメントに基づいて求める第２電磁場処理手段と、任意
の時間間隔での各粒子の状態を第１および第２の電磁場
および各粒子の双極子モーメントに基づいて求め、記憶
手段に記憶された各粒子の状態を更新する粒子状態処理
手段とを有し、双極子モーメント処理手段が、初回の処
理時には第２の電磁場を無視して第１の電磁場に基づい
て各粒子の双極子モーメントを求め、以降の処理時には
前回の処理時に求められた第２の電磁場と第１の電磁場
との両者に基づいて各粒子の双極子モーメントを求める
ように構成されているものである。

【００１１】

【作用】この発明においては、初回の処理時には第２の
電磁場を無視して第１の電磁場に基づいて各粒子の双極
子モーメントを求めるが、以後の処理においては、第１
および第２の電磁場に基づいて各粒子の双極子モーメン
トが求められる。

【００１２】

【実施例】図１はこの発明の一実施例による粒子挙動解
析装置のシステム構成を示すブロック図である。

【００１３】図を参照して、粒子挙動解析装置は、初期
条件入力装置１と、情報の一時記憶装置２（以下「メモ
リ」と称する）と、演算装置３と、計算結果出力装置４
とから構成される。

【００１４】初期条件入力装置１は、粒子の初期位置や
計算条件を入力するための手段であり、具体的には磁気
媒体に記録されたデータを読出す装置、カードおよびＣ
ＲＴ画面上で位置を指示する装置のいずれかが該当す
る。

【００１５】情報の一時記憶装置２は、入力された情報
または計算結果を一時的に格納する手段であり、具体的
には半導体記憶装置または磁気記憶装置のいずれかが該
当する。この装置には情報を格納する領域にアドレス
（番地）が割り振られており、アドレスを指定すること
でその領域に格納されているデータを読み出すことがで
きる。

【００１６】演算装置３は、一時記憶装置２に格納され
た情報を利用して、計算処理または比較判断処理を行な
うための演算手段であり、ＣＰＵを中心として構成され
る。

【００１７】計算結果出力装置４は、一時記憶装置２に
格納された情報を出力する手段であり、具体的にはＣＲ
Ｔ表示装置、プリンタ、ＸＹプロッタおよび磁気媒体へ
の記録装置のうち少なくとも１つが該当する。

【００１８】図２は図１の粒子挙動解析装置における粒
子の挙動計算の処理の流れを示したフローチャートであ
る。

【００１９】まず、ステップＳ１において、粒子の初期
配置や場の状態等が入力され、ステップＳ２において、
入力されたデータがメモリに格納される。そしてステッ
プＳ３において、メモリに格納されたデータに基づいて
粒子が場の状態によって受ける力が計算され、ステップ
Ｓ４において、粒子が受ける力に基づいて微小時間後の
粒子の位置と速度が計算される。なお粒子の位置と速度
とはそのつど必要に応じて出力することが可能である。

【００２０】ステップＳ５において、シミュレーション
の終了条件が満たされたか否かが判別される。すなわ
ち、所定の時間分のシミュレーションが行なわれたかま
たはシミュレーションの対象とすべき領域外に粒子が飛
び出したような場合が、シミュレーションの終了条件を
満たしたことになる。この終了条件が満たされていない
場合は、フローはステップＳ３に戻り、ステップＳ３お
よびステップＳ４の処理が繰返される。一方、シミュレ
ーションの終了条件が満たされた場合は、ステップＳ６
において、最終的に計算された粒子の位置と速度とを出
力して解析処理は終了する。

【００２１】図３は図２のステップＳ１においてなされ
る入力装置によるデータの入力手順を示したフローチャ
ートである。

【００２２】データの入力処理ルーチンに入るとステッ
プＳ１１において、粒子番号ｉが初期化され、ステップ
Ｓ１２において、粒子位置ベクトルＰ_os（ｉ）（以下
「Ｐ_os（ｉ）↓」として示す。他のベクトルの表示も同
様である）および粒子速度ベクトルＶ_el（ｉ）↓の初期
値が入力される。さらにステップＳ１３において、粒子
の半径、密度、透磁率、誘電率、磁気双極子モーメン
ト、電荷量等（これらを総称して「粒子物性Ｐ
_rop （ｉ）」と称する）の初期値を入力する。

【００２３】そしてステップＳ１４において、すべての
粒子について上記のデータが入力されたか否かが判別さ
れる。すべての粒子についての入力がされていない場合
は、ステップＳ１７において、粒子番号を１インクリメ
ントした後フローはステップＳ１２に戻り、すべての粒
子についてのデータの入力が終了するまでこの処理が繰
返される。

【００２４】一方、ステップＳ１４において、すべての
粒子についての入力が終了した場合は、ステップＳ１５
において、対象空間の初期状態（磁場、電場等）を、
「場のテーブル」として入力する。そしてステップＳ１
６において、シミュレーションの計算時間を「ｔ_end 」
として入力し、解析処理の時間間隔を表わす時間刻みを
「Δｔ」として入力して、データの入力処理は終了す
る。

【００２５】図４は図１の演算装置３の具体的構成を示
したブロック図である図を参照して、演算装置は場から
の力を計算する手段３１と、他の粒子による力を計算す
る手段３２と、位置・速度・粒子物性を更新する手段３
３と、終了条件の判定手段３４と、加算手段３５とから
構成されている。

【００２６】メモリ２から各種データＰ_os（ｉ）↓、Ｖ
_el（ｉ）↓、Ｐ_rop （ｉ）、場のテーブル、ｔ_end およ
びΔｔが入力される。場からの力を計算する手段３１
は、データＰ_os（ｉ）↓、Ｐ_rop （ｉ）および場のテー
ブルのデータによって計算を行ない、場からの力とし
て、Ｆ_field （ｉ）を他の粒子による力を計算する手段
３２と加算手段３５に出力する。他の粒子による力を計
算する手段３２は、入力されたデータＰ_os（ｉ）↓およ
びＰ_rop （ｉ）に基づいて場からの力を考慮しつつ計算
し、他の粒子による力としてＦ_ibead （ｉ）↓を加算手
段３５に出力する。加算手段３５において、場からの力
と他の粒子による力とが合成され、粒子が受ける力とし
てＦ（ｉ）↓が，粒子の位置等を更新する手段３３に出
力される。一方、更新する手段３３には、さらにデータ
Ｐ_os（ｉ）、Ｖ_el（ｉ）、Ｐ_rop （ｉ）およびΔｔが入
力され、これらのデータに基づいて微小時間経過後の粒
子の位置や速度等のデータを更新する。そして、更新さ
れたデータＰ_os（ｉ）↓およびＶ_el（ｉ）↓を、メモリ
２に出力する。

【００２７】一方、更新する手段３３から終了条件判定
手段３４には、データＰ_os（ｉ）↓、Ｖ_el（ｉ）↓、Ｐ
_rop （ｉ）およびΔｔが出力され、判定手段３４は、経
過時間のデータｔ_end に基づいて終了条件が満たされた
か否かを判定する。終了条件が満たされていないとき
は、判定手段３４は、更新されたデータＰ_os（ｉ）↓、
Ｖ_el（ｉ）↓およびＰ_rop （ｉ）を、計算する手段３１
および３２に戻し、上記と同様の計算を繰返す。

【００２８】図５〜図７は、図２のステップＳ２におい
て入力されたデータがメモリに格納された状態を示した
図である。

【００２９】図５を参照して、シミュレーション領域内
の粒子個数ｎ_B をメモリ上の任意の（１）アドレスＡ_nB
に格納する。以後、粒子個数ｎ_B が必要となるときに
は、アドレスＡ_nBの値を利用する。

【００３０】次に、対象となったｎ_B 個の粒子について
（２）粒子テーブルを作成する。粒子テーブルはメモリ
上の連続したアドレスに、粒子番号１〜ｎ_B の粒子につ
いて４つずつのアドレス情報、すなわち物性アドレス、
位置アドレス、計算量アドレスおよび近接粒子アドレス
を各４バイトずつのデータとして格納したものである。

【００３１】物性アドレスは粒子の半径、密度、比透磁
率、剛性係数およびダンピング係数を連続して格納して
いる（４）物性テーブルとしてのメモリ領域の先頭アド
レスの情報である。なお、剛性係数およびダンピング係
数は、粒子同士の衝突の際の力を計算するものであり、
本発明とは直接の関係を有さないので以後これらについ
ては言及しない。

【００３２】位置アドレスは、粒子の座標および粒子の
速度を連続して格納している（５）粒子位置テーブルと
してのメモリ領域の先頭アドレスの情報である。粒子位
置テーブルには、具体的にはｘ座標、ｙ座標、ｘ方向の
速度およびｙ方向の速度のデータが格納されている。な
お、物性テーブルおよび粒子位置テーブルの内部のアド
レスは連続したものでなくてもよい。

【００３３】なお、粒子位置テーブルはｘ，ｙの２次元
座標でその位置および速度を表しているが、三次元とし
ては、ｘ，ｙ，ｚ座標を用いることによって可能であ
り、他のデータは上記と同様のデータで処理することが
できる。

【００３４】なお、計算量アドレス、近接粒子のリスト
アドレスについては後述する。このように粒子の情報を
格納することで、たとえば３番目の粒子（粒子番号３）
の比透磁率は次のように参照できる。

【００３５】ここで１つのアドレス情報に４バイトのメ
モリ、１つの数値データに８バイトのメモリを使用する
とする。その場合、粒子テーブルの先頭アドレスはＡで
あるから、「３番目の粒子の物性アドレス」が格納され
ているアドレスはＡ＋３２となる。このアドレスのデー
タがたとえばＭ₂ であったなら、Ｍ₂ のアドレスから始
まるデータ列が３番目の粒子の物性データである。した
がって、３番目の粒子の比透磁率はＭ₂ ＋１６からの８
バイトのデータを参照すればよい。

【００３６】次に、場に関する有限要素法計算のための
データを図６に示すようにメモリに格納する。

【００３７】有限要素法はたとえば図８に示すように、
多角形の要素１〜５とその頂点である節点１〜５を利用
して計算を行なうことが一般的である。複雑に変化して
いる場の状態を計算するために１つ１つの要素とその要
素を構成する節点の特性値の間に成り立つ方程式を作成
し、それを連立させて全節点の特性値を求めることがで
きる。

【００３８】有限要素法の要素については多くの種類が
あるが、ここでは最も単純な三角形一次要素を想定して
説明する。また、特性値として磁気ポテンシャル値を使
用する。

【００３９】図６を参照して、全節点の数ｎ_nodeを
（１）アドレスＡｎ_nodeに格納する。続いて、任意のア
ドレスＡ_nodeからｎ_node個数のアドレスを収める連続し
た領域を（２）節点テーブルとして確保する。その領域
にはＡ_nodeアドレスから順番に１，２，３…の節点番号
に対応する節点位置テーブルへのアドレスが格納されて
いる。この（３）節点位置テーブルには節点テーブルに
示されたアドレスから順にｘ座標、ｙ座標、磁気ポテン
シャル値が格納されている。

【００４０】次に全要素の数ｎ_elemを（４）アドレスＡ
ｎ_elemに格納する。続いて、任意のアドレスＡ_elemから
ｎ_elem個数のアドレスを収める連続した領域を（５）要
素テーブルとして確保する。その領域はＡ_elemアドレス
から順番に１，２，３…の要素番号に対応する要素構成
テーブルへのアドレスが格納されている。この（６）要
素構成テーブルには、ある要素を構成する３つの節点の
番号が順に格納されている。この３つの節点番号はその
３点が反時計回りに結ばれるようになっている。

【００４１】このように節点と要素に関する情報を格納
することで、たとえば５番目の要素の２番目の節点のｘ
座標は以下のように読み出すことができる。

【００４２】まず、５番目の要素の情報を格納している
アドレスは、（５）要素テーブルからＡ_elem＋１６のア
ドレスと判明する。したがって、２番目の節点の番号
は、（Ａ_elem＋１６）＋４のアドレスの値を読出せばよ
い。このアドレスの値を基に（６）要素構成テーブルか
ら節点２の値を読出す。今この値が３であったとすれ
ば、（２）節点テーブルより、Ａ_node＋８のアドレスか
ら節点位置の座標が格納されていることがわかる。
（３）節点位置テーブルはｘ座標から始まっているか
ら、Ａ_node＋８のアドレスの内容を読出すと、その値は
５番目の要素の２番目の節点のｘ座標となっている。

【００４３】図７を参照して、アドレスＡ_t に経過時間
ｔ_end が格納され、アドレスＡ_stepに時間刻みΔｔが格
納される。

【００４４】以上の方法により、入力されたデータをメ
モリに格納し、任意のデータを読出すことができる。

【００４５】次に、図２のステップＳ３における粒子が
場の状態によって受ける力の計算について詳細に説明す
る。

【００４６】この実施例では、磁場中での磁気双極子の
運動について例を示す。粒子（磁気双極子）が磁場から
受ける力、すなわち並進方向の力Ｆ↓および回転方向の
力のモーメントＴ↓は、 Ｆ↓＝（ｍ↓・∇）Ｂ↓ …（１） Ｔ↓＝ｍ↓×Ｂ↓ …（２） ここで、 ｍ↓：粒子の磁気双極子モーメントベクトル Ｂ↓：粒子位置の磁束密度 である。

【００４７】したがって、ｍ↓、Ｂ↓およびＢ↓の偏微
分（∇Ｂ）が計算できれば、磁場から受ける力を計算す
ることができる。このうちＢ↓の偏微分値はＢ↓の差分
によって近似できるので、最低限必要な値はｍ↓とＢ↓
である。

【００４８】ｎ_B 個の粒子の内ｊ番目の粒子の磁気双極
子モーメントｍ_j ↓は、

【００４９】

【数１】

【００５０】ここで、 ａ：粒子半径 μ₀ ：真空の透磁率 μ：粒子の比透磁率 Ｂ_j ↓：粒子位置の磁束密度 として計算される。

【００５１】一方、粒子位置の磁束密度は、本来ならば
ある時点での粒子の位置、物性および粒子以外に起因す
る場の状態を併わせて計算されるものである。しかし、
粒子位置は各計算時間毎に変化するため、これらをすべ
て考慮するのは時間的な効率が非常に悪い。

【００５２】そこで、ｊ番目の粒子位置での磁束密度Ｂ
_j ↓が、

【００５３】

【数２】

【００５４】ここで、 Ｂ′_j ↓：粒子に起因しない磁束密度 Ｂ_kj↓：ｋ番目の粒子によるｊ番目の粒子位置での磁束
密度 という式で表せるものと仮定して計算の高速化を図る。
この仮定は粒子の磁気特性が線形である場合に正しくな
る。

【００５５】（４）式のＢ′_j ↓は粒子に起因しない磁
束密度であるが、図６で示される場のデータから次のよ
うに計算する。

【００５６】まず、ｊ番目の粒子の座標（ｘ，ｙ）を粒
子テーブルと粒子位置テーブルとから読出す。次に、粒
子座標を含む有限要素を探すために図９で示されたフロ
ーに基づいて処理を行なう。

【００５７】図９を参照して、まずステップＳ３１にお
いて要素番号ｎを初期化し、ステップＳ３２において、
要素番号ｎの要素を構成する３節点の番号をメモリから
取得する。次にステップＳ３３において、得られた３節
点の座標（ｘ₁ ，ｙ₁ ）、（ｘ₂ ，ｙ₂ ）および
（ｘ₃ ，ｙ₃ ）を節点番号に基づいてメモリから取得す
る。そしてステップＳ３４において、下記の３つの式に
ついて計算を行なう。

【００５８】Ａ＝（ｘ₁ −ｘ）（ｙ₂ −ｙ）−（ｘ−ｘ
₂ ）（ｙ−ｙ₁ ） Ｂ＝（ｘ₂ −ｘ₁ ）（ｙ−ｙ₁ ）−（ｘ₁ −ｘ）（ｙ₁
−ｙ₂ ） Ｃ＝（ｘ−ｘ₂ ）（ｙ₁ −ｙ₂ ）−（ｘ₂ −ｘ₁ ）（ｙ
₂ −ｙ） そしてステップＳ３５において、上記の３つの式のＡ，
Ｂ、Ｃの値がいずれも０以上であるか否かが判別され
る。これらの式の値がいずれも０以上である場合は、ス
テップＳ３６において、その粒子座標はその要素の中に
あると見なされて処理は終了する。

【００５９】一方、３つの式の値がいずれも０以上でな
い場合は、ステップＳ３７において、要素番号が１イン
クリメントされた後、ステップＳ３８において、要素番
号ｎが最後の要素番号ｎ_elemになるまでステップＳ３２
〜ステップＳ３５の処理が繰返される。このようにし
て、該当の粒子座標を含む有限要素が検出される。

【００６０】そして、粒子座標を含む要素が見つかった
場合、Ｂ′_j ↓のｘ成分Ｂ′_jxおよびｙ成分Ｂ′_jyを以
下のように計算する。

【００６１】

【数３】

【００６２】一方、（４）式のＢ_kj↓は以下の式によっ
て計算される。

【００６３】

【数４】

【００６４】ここで、 ｒ_kj↓：ｋ番目の粒子からｊ番目の粒子への位置ベクト
ル これはｋ番目の粒子の座標を（ｘ_k ，ｙ_k ）、ｊ番目の
粒子の座標を（ｘ_j ，ｙ_j ）とすると ｒ_kj↓＝（ｘ_j −ｘ_k ，ｙ_j −ｙ_k ）となる。

【００６５】これらの式を組合わせれば、すべての粒子
の磁気双極子モーメントとすべての粒子が受ける力が計
算できる。しかし、現実的には（３）式のｍ↓を求める
にはＢ↓が必要であり、（７）式のＢ_kj↓を求めるのに
ｍ↓が必要であることから、すべてのｍ↓、Ｂ↓を求め
るためには膨大な繰返しを伴なう収束計算が要求され
る。

【００６６】そこで、この発明の一実施例による粒子挙
動回析装置においては、ｍ↓を求める際、（３）式のＢ
↓として最初の１ステップの計算に限り、Ｂ′↓を使
い、第２ステップ以降は１計算ステップ前のＢ↓、すな
わち（４）式によって得られた値を用いることとする。
計算ステップの間隔が十分に小さいときは、その間の粒
子の移動量も小さいためこの手法を利用することができ
る。

【００６７】ここまでの計算において求められた各粒子
のｍ↓、Ｂ↓、Ｆ↓（粒子が受ける力）の各値は、メモ
リ上に図１０の形式で計算量テーブルとして格納されて
いる。この値は図５の粒子テーブルにある計算量アドレ
スを用いて読出すことができる。

【００６８】なお、（７）式によって求められる粒子に
起因する磁束密度であるが、この式から明らかなように
Ｂ_kj↓の大きさはｒ_kj↓が大きくなるにつれて急激に小
さくなる。そこで、Ｂ_kj↓を計算する粒子の組合わせを
粒子間距離の近いものだけに限定することで、さらに計
算の高速化が実現する。

【００６９】この目的で粒子間距離の小さい粒子対を各
計算ステップの最初に選び出し、リストしたものが図１
１である。

【００７０】図を参照して、図５の粒子テーブルにおけ
る粒子番号の各々に対応した近接粒子アドレス毎に近接
粒子として扱うべき粒子番号が格納されている。このよ
うなリストをメモリ上に置いて、この組合わせについて
のみＢ_kj↓を計算し、残りは０として処理すればよい。

【００７１】次に、図２のステップＳ４における微小時
間後の粒子の位置と速度の計算について説明する。

【００７２】先に述べた計算によって算出された各粒子
が受ける力に基づいて、粒子の位置と速度を次のように
更新する。ここでは前進オイラー法による例を示すが、
すべての粒子について次の処理を行なう。

【００７３】まず、加速度α_x をＦ_x ／ｍａｓｓ、α_y
をＦ_y ／ｍａｓｓとして計算する。なお、Ｆは先に述べ
た計算によって計算された粒子が受ける力であり、ｍａ
ｓｓはメモリの物性テーブルに格納された各粒子の物性
データ（密度と半径）から計算される質量である。すな
わち、 ｍａｓｓ＝４πａ³ ρ／３ ここで ａ：粒子の半径 ρ：粒子の密度 次に、粒子の速度Ｖ_x に“α_x ×（Δｔ）”、Ｖ_y に
“α_y ×（Δｔ）”を加えることによって速度が更新さ
れる。

【００７４】さらに、粒子の位置ｘに“Ｖ_x ×（Δｔ）
＋α_x ×（Δｔ）² ／２”、ｙに“Ｖ_y ×（Δｔ）＋α
_y ×（Δｔ）² ／２”を加えることによって粒子位置が
更新される。なお、ここでの速度Ｖ_x およびＶ_y は更新
される前の値である。

【００７５】更新された位置および速度はメモリの粒子
の位置テーブルに格納されると共に必要に応じて出力装
置によって外部に出力される。

【００７６】図１２は図２のステップＳ６の出力ルーチ
ンの具体的内容を示したフローチャートである。

【００７７】この処理に入るとまずステップＳ６１にお
いて粒子番号ｉが初期化され、ステップＳ６２において
粒子位置Ｐ_os（ｉ）↓および粒子速度Ｖ_el（ｉ）↓が出
力される。そしてステップＳ６３においてすべての粒子
についてその粒子位置と粒子速度が出力されたか否かが
判別される。すべての粒子について出力されていない場
合にはステップＳ６４において、粒子番号ｉが１インク
リメントされた後ステップＳ６２およびステップＳ６３
の処理が繰返される。そしてすべての粒子について出力
処理が終わった時点で出力処理は終了する。

【００７８】なお、上記実施例における時間刻み幅Δｔ
は初期値として与えた値を使うことも可能であるが、以
下のように計算の進行に伴って変化させていくこともで
きる。

【００７９】たとえば、「すべての粒子の１ステップ
（時間刻み幅）当たりの移動距離が最小粒子半径の１％
以下になるように規制する」場合は以下となる。

【００８０】１） 最も小さい粒子半径の１％を許容最
大移動距離Ｄｍａｘとする。 ２） 初回の計算ステップにおいては初期設定値を時間
刻み幅Δｔとして使用する。

【００８１】そして次回以降の計算ステップにおいて
は、 ３） 前回のステップにおいて計算された粒子の移動距
離Ｄ_-1（ｉ）は、

【００８２】

【数５】

【００８３】そして、この移動距離Ｄ_1-（ｉ）の最大値
ｍａｘ（Ｄ_-1（ｉ））と、前ステップの時間刻み幅Δｔ
_-1を使って、今回の仮の時間刻み幅Δｔを、 Δｔ＝Δｔ_-1・Ｄｍａｘ／ｍａｘ（Ｄ_-1（ｉ）） として決める。

【００８４】４） このΔｔを使って、現ステップの計
算を行なう。 ５） その結果、求められた粒子の移動距離Ｄ（ｉ）
は、

【００８５】

【数６】

【００８６】そしてこの移動距離Ｄ（ｉ）の最大値ｍａ
ｘ（Ｄ（ｉ））を計算する。 ６） ｍａｘ（Ｄ（ｉ））≦Ｄｍａｘであれば、仮の時
間刻み幅Δｔは適正な値であったと判断し、計算結果を
確定し、次のステップに進む。

【００８７】７） 一方、ｍａｘ（Ｄ（ｉ））＞Ｄｍａ
ｘであれば、 Δｔ′＝Δｔ・Ｄｍａｘ／ｍａｘ（Ｄ（ｉ）） で計算されるΔｔ′を新しい時間刻み幅Δｔとし、今回
のステップで得られた結果を破棄し、再計算を行なう。

【００８８】このような処理によって計算の時間刻みΔ
ｔを可変にすることによって、次のような効果が得られ
る。

【００８９】ａ．粒子の速度が大きくなった場合には、
時間刻み幅が小さくなることによって、誤差の拡大が防
止される。

【００９０】ｂ．粒子の速度が小さいときは、時間刻み
幅を大きくとることによって、高速計算が可能となる。

【００９１】すなわち、対象粒子集団の状態に合わせて
効率のよい計算を行なうことが可能となる。

【００９２】また、上記実施例においては、磁場につい
ての粒子の挙動を説明しているが、電場についても同様
に本発明は適用できる。図１３および図１４はこのよう
な観点で磁場と電場とを対応させた表である。図に示し
たように１／μ₀ がε₀ に変わっている以外はほぼ同じ
対応をしている。ところで磁場において「粒子に起因し
ない場」の源は「永久磁石」や「電磁石」であるが、電
場ではこれに対応してチャージャ線などの放電電極や、
感光体等のプレート電極等が挙げられる。

【００９３】

【発明の効果】この発明は以上説明したとおり、初回の
処理時には第２の電磁場を無視して第１の電磁場に基づ
いて各粒子の双極子モーメントを求めるが、以後の処理
においては第１および第２の電磁場に基づいて各粒子の
双極子モーメントを求めるので、解析精度を確保しつ
つ、高速な計算処理が可能となる。

【図面の簡単な説明】

【図１】この発明の一実施例による粒子挙動解析装置の
システム構成を示すブロック図である。

【図２】図１の装置における解析処理の流れを示したフ
ローチャートである。

【図３】図２のステップＳ１の入力ルーチンの具体的内
容を示したフローチャートである。

【図４】図１の演算装置３の具体的構成を示すブロック
図である。

【図５】図１の一時記憶装置２における粒子に関するデ
ータの格納状態を示した図である。

【図６】図１の一時記憶装置２における場に関するデー
タの格納状態を示した図である。

【図７】図１の一時記憶装置２における解析処理の経過
時間や時間刻み幅情報の格納状態を示した図である。

【図８】この発明の一実施例に用いる有限要素法の要素
と節点の一例を示した図である。

【図９】図２のステップＳ３の計算に関連して、対象と
する粒子を含む有限要素を探す処理を示すフローチャー
トである。

【図１０】図１の一時記憶装置２における計算量テーブ
ルの格納状態を示した図である。

【図１１】図１の一時記憶装置２における近接粒子リス
トテーブルのデータの格納状態を示した図である。

【図１２】図２のステップＳ６の出力ルーチンの具体的
内容を示したフローチャートである。

【図１３】本願発明の一実施例に対応して、磁場と電場
とを比較した表の一部である。

【図１４】この発明の一実施例に対応して、磁場と電場
とを比較した表の他の一部である。

【符号の説明】

１ 初期条件入力装置 ２ 情報の一時記憶装置 ３ 演算装置 ４ 計算結果出力装置 ３１ 場からの力を計算する手段 ３２ 他の粒子による力を計算する手段 ３３ 位置、速度、粒子物性を更新する手段 ３５ 加算手段 なお、図において同一符号は同一部分または相当部分を
示す。

───────────────────────────────────────────────────── フロントページの続き (51)Int.Cl.⁶ 識別記号 庁内整理番号 ＦＩ 技術表示箇所 // Ｇ０３Ｇ 15/08 ５０７ Ｋ (72)発明者 加藤 充 大阪市中央区安土町二丁目３番13号大阪国 際ビル ミノルタカメラ株式会社内 (72)発明者 相澤 龍彦 東京都世田谷区南烏山１−１−21−708 (72)発明者 田村 茂之 東京都町田市三輪緑山１−34−17 (72)発明者 岩井 隆 神奈川県逗子市池子２−19−64

Claims (1)

【特許請求の範囲】
1. 【請求項１】 有限な領域内での電磁場の影響下におけ
   る複数の粒子の挙動を解析する粒子挙動解析装置におい
   て、 各粒子の初期状態および電磁場の状態を入力するための
   入力手段と、 各粒子の状態を記憶する記憶手段と、 任意の時間間隔での各粒子の挙動を演算し、前記記憶手
   段に記憶された各粒子の状態を更新する粒子挙動演算手
   段と、 前記記憶手段に記憶された各粒子の状態を出力するため
   の出力手段とを備え、 前記粒子挙動演算手段が、 前記任意の時間間隔での各粒子の位置における粒子以外
   に起因する第１の電磁場を、前記入力手段から入力され
   た電磁場の状態に基づいて求める第１電磁場処理手段
   と、 前記任意の時間間隔での各粒子の双極子モーメントを求
   める双極子モーメント処理手段と、 前記任意の時間間隔での各粒子の位置における他の粒子
   に起因する第２の電磁場を、各粒子の双極子モーメント
   に基づいて求める第２電磁場処理手段と、 前記任意の時間間隔での各粒子の状態を前記第１および
   第２の電磁場および各粒子の双極子モーメントに基づい
   て求め、前記記憶手段に記憶された各粒子の状態を更新
   する粒子状態処理手段とを有し、 前記双極子モーメント処理手段が、初回の処理時には第
   ２の電磁場を無視して第１の電磁場に基づいて各粒子の
   双極子モーメントを求め、以降の処理時には前回の処理
   時に求められた第２の電磁場と第１の電磁場との両者に
   基づいて各粒子の双極子モーメントを求めるように構成
   されている、粒子挙動解析装置。