JP2013183551A

JP2013183551A - 解析装置

Info

Publication number: JP2013183551A
Application number: JP2012046390A
Authority: JP
Inventors: Daiji Ichijima; 大路市嶋; Shuji Miyazaki; 修司宮崎
Original assignee: Sumitomo Heavy Industries Ltd
Current assignee: Sumitomo Heavy Industries Ltd
Priority date: 2012-03-02
Filing date: 2012-03-02
Publication date: 2013-09-12
Anticipated expiration: 2032-03-02
Also published as: EP2637117A3; JP5839473B2; EP2637117A2; US9322886B2; US20130231879A1

Abstract

【課題】シミュレーションにおいて磁気的な現象を好適に扱える解析技術を提供する。
【解決手段】解析装置１００は、仮想空間内に定義される粒子系の情報を取得する粒子系取得部１０８と、粒子系の粒子に磁気モーメントを付与する磁気モーメント付与部１１０と、磁気モーメント付与部１１０によって磁気モーメントが付与された粒子を含む粒子系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部１２０と、数値演算部１２０における演算結果を使用して、粒子系が生成する磁場を演算する磁場演算部１３２と、を備える。
【選択図】図１

Description

本発明は、粒子系を解析する解析装置に関する。

近年、コンピュータの計算能力の向上に伴い、モータなどの電気機器の設計開発の現場では磁場解析を取り入れたシミュレーションがよく使用されるようになっている。シミュレーションを使用すると、実際にプロトタイプを製作しなくてもある程度の評価が可能となるので、設計開発のスピードが向上しうる。

例えば特許文献１には、磁場解析を実行する演算処理装置を具えるモータ解析装置が記載されている。演算処理装置は、ユーザの操作に基づく外部指令に応じて、有限要素法による静磁場解析やマクスウェルの応力法によるトルク計算を実行する。有限要素法による静磁場解析のためにメッシュ分割が行われる。メッシュ分割は、コア領域及びハウジング領域、さらには外部空気層領域についても行われる。
磁場解析の他の手法は差分法および磁気モーメント法を含む。

特開平１１−１４６６８８号公報特開２００９−３７３３４号公報

Erik Bitzek他、「Structural Relaxation Made Simple」、Physical Review Letters、２００６年１０月２７日、９７

これらの手法では解析対象をメッシュ分割する。しかしながら、解析対象が流体の場合固定メッシュが使われ、移動境界を伴うと解析が困難になる。また、解析対象が磁性体や塑性体の大変形を伴う場合もメッシュ分割は困難である。特に、モータなどの回転体を含む磁場解析でメッシュ分割を行おうとする場合、メッシュが重合（一重なら大変形）するため、解析は困難であった。

本発明はこうした課題に鑑みてなされたものであり、その目的は、シミュレーションにおいて磁気的な現象を好適に扱える解析技術の提供にある。

本発明のある態様は、解析装置に関する。この解析装置は、仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する磁気モーメント付与部と、磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部と、数値演算部における演算結果を使用して、粒子系に関連する磁気的な物理量を演算する磁気的物理量演算部と、を備える。

この態様によると、粒子に付与された磁気モーメントを使用して磁気的な物理量を導出できる。

なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。

本発明によれば、シミュレーションにおいて磁気的な現象を好適に扱うことができる。

実施の形態に係る解析装置の機能および構成を示すブロック図である。図１の粒子データ保持部の一例を示すデータ構造図である。図１の解析装置における一連の処理の一例を示すフローチャートである。ｎ＝１次までの展開に対応する計算結果を示す図である。ｎ＝２次までの展開に対応する計算結果を示す図である。ｎ＝３次までの展開に対応する計算結果を示す図である。ｎ＝４次までの展開に対応する計算結果を示す図である。球面調和関数の次数と計算値との関係を示す図である。ヒステリシスカーブを示す図である。永久磁石（１．０Ｔ）による磁化の計算結果を示す図である。永久磁石（１．０Ｔ）による磁化の計算結果を示す図である。永久磁石（１．０Ｔ）による磁化の計算結果を示す図である。

以下、各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。

古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）や、量子分子動力学法（Quantum Molecular Dynamics Method）や、ＭＤ法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法（Renormalized Molecular Dynamics、以下ＲＭＤ法と称す）に基づくシミュレーションが知られている（例えば、特許文献２参照）。正確にはＭＤ法やＲＭＤ法は運動論的手法（物理量の算出には統計力学を使う）であり、粒子法は、連続体を記述する微分方程式を離散化する手法であり、別ものであるが、ここではＭＤ法やＲＭＤ法も粒子法と呼ぶ。

粒子法は静的な現象だけでなく流れなどの動的な現象をも取り扱えるので、主に静的な現象を解析対象とする上述の有限要素法などに代わるシミュレーション手法として注目されている。

粒子法には、連続体を粒子で離散化することにより解析対象の粒子系を得る、という微分的な見方がある。例えば、粒子法において流体を扱う場合、ナビエストークス（Navier-Stokes）方程式を粒子で離散化することが多い。

一方、粒子法の別の見方として、多くの粒子を集めて連続体を形成する、という積分的な見方もある。これは例えば、小さな鉄の粒を集めて固めて大きな鉄球を形成するという見方である。

一般に、多くの磁気モーメントが存在する空間内のある点の磁場を求める場合、重ね合わせの原理により、各磁気モーメントがその点に作る磁場を磁気モーメントに亘って足し合わせる。本発明者は、このような磁気モーメントの集まりから磁場を求める手法と、積分的な見方をした場合の粒子法と、の親和性を独自に見い出し、粒子法における粒子に磁気モーメントを付与することに想到した。これにより、対流や大変形の解析に強いという粒子法の利点を維持しつつ、粒子法の適用範囲を磁場解析にまで広げることが可能となる。

図１は、実施の形態に係る解析装置１００の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。

本実施の形態ではＲＭＤ法に倣って粒子系を解析する場合について説明するが、繰り込みを行わないＭＤ法やＤＥＭ（Distinct Element Method）やＳＰＨ（Smoothed Particle Hydrodynamics）やＭＰＳ（Moving Particle Semi-implicit）などの他の粒子法に倣って粒子系を解析する場合にも、本実施の形態に係る技術的思想を適用できることは本明細書に触れた当業者には明らかである。

解析装置１００は入力装置１０２およびディスプレイ１０４と接続される。入力装置１０２は、解析装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。

解析装置１００は、粒子系取得部１０８と、磁気モーメント付与部１１０と、数値演算部１２０と、磁場演算部１３２と、表示制御部１１８と、粒子データ保持部１１４と、を備える。

粒子系取得部１０８は、入力装置１０２を介してユーザから取得する入力情報に基づき、１、２または３次元の仮想空間内に定義されるＮ（Ｎは自然数）個の粒子からなる粒子系のデータを取得する。粒子系はＲＭＤ法を使用して繰り込まれた粒子系である。

粒子系取得部１０８は、入力情報に基づき仮想空間内にＮ個の粒子を配置し、配置されたそれぞれの粒子に速度を付与する。粒子系取得部１０８は、配置された粒子を特定する粒子ＩＤと、その粒子の位置と、その粒子の速度と、を対応付けて粒子データ保持部１１４に登録する。

磁気モーメント付与部１１０は、入力装置１０２を介してユーザから取得する入力情報に基づき、粒子系取得部１０８によって取得された粒子系の粒子に磁気モーメントを付与する。例えば、磁気モーメント付与部１１０は、ディスプレイ１０４を介してユーザに、粒子系の粒子の磁気モーメントの入力を求める。磁気モーメント付与部１１０は、入力された磁気モーメントを粒子ＩＤと対応付けて粒子データ保持部１１４に登録する。

以下では粒子系の粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。

数値演算部１２０は、粒子データ保持部１１４によって保持されるデータが表す粒子系の各粒子の運動を支配する支配方程式を数値的に演算する。特に数値演算部１２０は、離散化された粒子の運動方程式にしたがった繰り返し演算を行う。粒子の運動方程式は磁気モーメントに依存する項を有する。
数値演算部１２０は、力演算部１２２と、粒子状態演算部１２４と、状態更新部１２６と、終了条件判定部１２８と、を含む。

力演算部１２２は粒子データ保持部１１４によって保持される粒子系のデータを参照し、粒子系の各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。粒子に働く力は、磁気モーメント間の相互作用に基づく力を含む。力演算部１２２は、粒子系のｉ番目（１≦ｉ≦Ｎ）の粒子について、そのｉ番目の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。

力演算部１２２は、各近接粒子について、その近接粒子とｉ番目の粒子との間のポテンシャルエネルギ関数およびその近接粒子とｉ番目の粒子との距離に基づいて、その近接粒子がｉ番目の粒子に及ぼす力を演算する。特に力演算部１２２は、その近接粒子とｉ番目の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。力演算部１２２は、近接粒子がｉ番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、ｉ番目の粒子に働く力を算出する。

粒子状態演算部１２４は粒子データ保持部１１４に保持される粒子系のデータを参照し、粒子系の各粒子について、離散化された粒子の運動方程式に力演算部１２２によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。

粒子状態演算部１２４は、力演算部１２２によって演算された力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき所定の微小な時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって演算された力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。

粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。

状態更新部１２６は、粒子データ保持部１１４に保持される粒子系の各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。

終了条件判定部１２８は、数値演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、外部から終了の指示を受け付けたことや、粒子系が定常状態に達したことである。終了条件判定部１２８は、終了条件が満たされる場合、数値演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を力演算部１２２に戻す。すると力演算部１２２は、状態更新部１２６によって更新された粒子の位置で再び力を演算する。

磁場演算部１３２は、数値演算部１２０における演算が終了すると、粒子データ保持部１１４によって保持されるデータが表す粒子系が生成する磁場を演算する。磁場演算部１３２は、粒子データ保持部１１４によって保持される各粒子の磁気モーメントに基づいて磁場を演算する。磁場は、粒子系に関連する磁気的な物理量である。磁場演算部１３２は磁場の代わりにまたは磁場に加えて磁束密度や磁化を演算してもよい。
また、磁場演算部１３２は、数値演算部１２０における繰り返し演算の各ステップにおいて、粒子の位置に生成される磁場を演算してもよい。

表示制御部１１８は、粒子データ保持部１１４に保持されるデータが表す粒子系の各粒子の位置、速度、磁気モーメントに基づき、ディスプレイ１０４に粒子系の時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。

図２は、粒子データ保持部１１４の一例を示すデータ構造図である。粒子データ保持部１１４は、粒子ＩＤと、粒子の位置と、粒子の速度と、粒子の磁気モーメントと、を対応付けて保持する。

上述の実施の形態において、保持部の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。

以上の構成による解析装置１００の動作を説明する。
図３は、解析装置１００における一連の処理の一例を示すフローチャートである。粒子系取得部１０８は、ＲＭＤ法に倣って繰り込まれた粒子系を取得する（Ｓ２０２）。磁気モーメント付与部１１０は、取得された粒子系の粒子に磁気モーメントを付与する（Ｓ２０４）。力演算部１２２は、粒子間の距離から粒子に働く力を演算する（Ｓ２０６）。粒子状態演算部１２４は、演算された力を含む粒子の運動方程式から粒子の速度と位置を演算する（Ｓ２０８）。状態更新部１２６は、粒子データ保持部１１４に保持される粒子の位置、速度を演算された位置、速度で更新する（Ｓ２１０）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２１２）。終了条件が満たされない場合（Ｓ２１２のＮ）、処理はＳ２０６に戻される。終了条件が満たされる場合（Ｓ２１２のＹ）、磁場演算部１３２は粒子系によって生成される磁場を演算する（Ｓ２１４）。

粒子法は一般に、流体や分離現象や切断加工などの動的な解析対象をより好適に扱える。
そして本実施の形態に係る解析装置１００によると、粒子法による解析と磁気モーメントに基づく磁場解析とを連成することができる。したがって、動的な解析対象を、磁場などの磁気的な性質も含めてより精度高くかつ短時間で解析することが可能となる。

例えば解析対象が砂鉄の場合、従来の有限要素法などのメッシュベースの解析手法ではメッシュを好適に設定することができないので、うまくシミュレーションすることが困難である。そこで本実施の形態に係る解析装置１００によると、砂鉄の磁性を考慮に入れた形で砂鉄を粒子法で好適に解析することができる。

また例えば解析対象がモータの場合、固定体に対して回転体が回転するので、やはりメッシュベースの解析手法は不向きであり静解析に止まる。そこで、本実施の形態に係る解析装置１００によると、回転体の回転を粒子法により高精度にシミュレーションしつつ、回転体と固定体との磁気的相互作用を動的に扱うことができる。

以下、解析装置１００で使用される解析手法の原理について説明する。

磁場はマクロなベクトルポテンシャル：
から以下の様に求まる。スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにある１個のスピンが作る磁場はスピンＳ_ｊｚ（磁気モーメントｍ_ｊｚ）を発生させる円電流を
とし、
ここで、
ｎ＝０のみを残すと一様に磁化：
した球の作る磁場と同等になる。
なお、真球のバルク材に本手法を適用する場合、ｎ＝０のみ採用すると、一様な磁化が得られず、磁気モーメント法と同様な結果になる。高次の項を加える事で改善できる。ＦＥＭは一様になる。

デカルト座標系に変換（
を使って、)すれば
ここで、ｓｉｎθ_ｊ＝０のときは、Ｂ_θ＝０から、Ｂ_ｘ＝Ｂ_ｙ＝０であることに注意する。磁気モーメントｍ_ｘ、ｍ_ｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場
も同様に求まる。

における磁場
、Ｎ個のスピンからの寄与、は
ここで、
は球内の平均を意味する。Ｐ−１３とＰ−１６を連立することで, 磁気モーメントｍ_ｐが求まる。磁性体外部に生じる磁場は、Ｐ−１３であり、磁性体内部に生じる磁場は、
Ｐ−１４は球体表面に生じる磁化による反磁場が考慮され、
を満足する。

Ｐ−１３はマトリクスが優対角でない(ｊ＝ｐが除外されている) から収束は遅い。そこで、以下の運動方程式から解を得る。
ここで、ｍ_ｖ＝５×１０^−２ｄｔ^２は仮想質量、γ＝ｍ_ｖ／（１０ｄｔ）は、減衰係数である。

ＦＩＲＥ（非特許文献１参照）を追加すればさらに高速化が可能。以下にＦＩＲＥのコードを示す。なお、粒子数８０２、残差＜１０^−８における計算で１０．０５[ｓ］である。ＦＩＲＥ無しでは、残差が１０^−５以下にならなかった。

以下にLegendre 関数と倍Legendre 関数の特別な場合を列挙する。（ｘ：＝ｃｏｓθ_ｊ）

以下に、ｘ軸方向にＢ_ｏｘ＝０．１［Ｔ］の一様な外場においた磁性球（帯磁率が９９９）の磁化
を計算した結果を示す。球を構成する粒子（原子）数は４００９個とし、ｆｃｃ構造に配置した。厳密解は
である。図４、図５、図６、図７に多重極展開の１次〜４次までの結果を示す。共に中心を通過する断面を表示した。多重極展開の次数を上げると, 厳密解と良く合う事が判る。次数が低いと収束が悪い。また計算時間の大半は球面調和関数の呼び出しである。

図４は、ｎ＝１次までの展開に対応する計算結果を示す図である。厳密解μ_ｏＭ_ｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３４６２３１[Ｔ］である。磁化は外場Ｂ_ｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。残差１ｅ−１３以下で、計算時間は約８０分であった。

図５は、ｎ＝２次までの展開に対応する計算結果を示す図である。厳密解μ_ｏＭ_ｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３２１８９１[Ｔ]である。磁化は外場Ｂ_ｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。残差１ｅ−８以下で、計算時間は約３時間であった。残差は１ｅ−９以下に落ちなかった。

図６は、ｎ＝３次までの展開に対応する計算結果を示す図である。厳密解μ_ｏＭ_ｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３１２４４４[Ｔ]である。磁化は外場Ｂ_ｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる.

図７は、ｎ＝４次までの展開に対応する計算結果を示す図である。厳密解μ_ｏＭ_ｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３１２２２２[Ｔ]である。磁化は外場Ｂ_ｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。約５時間後において、残差は１ｅ−９以下に落ちなかった。

図８は、球面調和関数の次数と計算値との関係を示す図である。次数を上げると厳密解に近づくことが分かる。

解析装置１００で使用される解析手法の原理を、次のように説明することもできる。

Ｍ_ｐは原子核の質量、ｍ_ｅは電子の質量である。ｅは原子核と電子の電荷(の絶対値)である。最終的に電荷は電流に置き換わり、電子の質量は原子核の質量と和して原子の質量Ｍとなるから、陽には現れない。δＶは磁性球の体積、ｍ_ｖは仮想質量、λ_ｉはラグランジェの未定常数である。厳密に拘束を満足させるには、ＳＨＡＫＥ法に代表される収束計算が必要であるが、本手法ではλ_ｉ＝１として、正解の周りの微小振動を誤差として扱う。

スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにあるスピンＳ_ｊｚ（磁気モーメントｍ_ｊｚ）を発生させる円電流は
であるから、磁場はマクロなベクトルポテンシャル：
から以下の様に求まる。
ここで、
ｎ＝０のみを残すと一様に磁化：
した球の作る磁場と同等になる。

デカルト座標系に変換(
を使って、）すれば
ここで、ｓｉｎθ_ｊ＝０のときは、Ｈ_θ＝０から、Ｈ_ｘ＝Ｈ_ｙ＝０であることに注意する。磁気モーメントｍ_ｘ、ｍ_ｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場
も同様に求まる。

におけるモーメントが感じる外磁場
(Ｎ個のスピンからの寄与) は
誘導磁場
は、
から求まる。
なお、繰り込みに際し、電気伝導度σは
今はδ＝２としている。

Ｑ−２１右辺のLagrange 微分
は、粒子法であるから
で置き換えることが出来る。ただし、磁化の並進移動は粒子の移動により考慮されるが、磁化ベクトルの回転は別に考慮する必要がある(ランダウ=リフシッツ、電磁気学、５１章)。

磁性体外部に生じる磁場は、Ｑ−２６であり、磁性体内部に生じる磁場は、

も、非線形であることに注意する。全外場
と磁化
との関係は、
である。Ｑ−２６とＱ−２７を連立することで、磁気モーメント
が求まる。

非線形磁化とヒステリシス
Ｑ−２７とＱ−２４から磁場
を消去する：
Ｑ−２８を解くことで
が得られる（図９にはＢ＝ｆ（Ｈ）の場合を示す）。

図９は、ヒステリシスカーブを示す図である。ヒステリシスカーブとＢ＋２μ_ｏＨ＝３Ｂ_ｏの交点からＭを求める。今
を永久磁化とするＢ−Ｈカーブ：
が与えられたとする。
と連立し
を得る。一般のＢ＝ｆ（Ｈ）に対しては、Newton-Laphson法で解く。

磁場の粒子化法
Ｑ−１９はマトリクスが優対角でない（ｊ＝ｐが除外されている）から収束は遅い。そこで、以下の運動方程式から解を得る。
ここで、ｍ_ｖは仮想質量である。推奨値は１００×２．５×１０^−１７ｄｔである。

もし静磁場解析に限れば、グリーン関数Ｇ（ｒ_ｐｊ）をγδ_ｐｊ／μ_ｏで置き換える。このとき、ｍ_ｖ＝５×１０^１ｄｔ^２、γ＝ｍ_ｖ／（１００ｄｔ）とする。ｍ_ｖを十分小さく取れば、正解の周りの微小な減衰振動が起こり、誤差を小さくできる。

模型を考える：
減衰項を持つ強制振動の解(Landau and Lifshitz,Mechanics,Sec.26) は
を考えれば、
と成る。即ちｍ_ｖを十分小さく取れば、解くべき本来の方程式
の解に一致する。

離散化
蛙飛び法（速度はｎとｎ＋１の中間点における評価）に従い離散化すれば、収束回数をｓと書き、
対角要素のみ未来を使えば、陽的解法に出来、約３０％程度効率改善になる。計算の大半は特殊関数の呼び出しである。

ＦＩＲＥによる高速化：静磁場解析にのみ有効
ＦＩＲＥ（非特許文献１参照）を追加すればさらに高速化が可能。以下にＦＩＲＥのコードを示す。なお、粒子数８０２、残差＜１０^−８における計算で１０．０５[ｓ]である。ＦＩＲＥ無しでは、残差が１０^−５以下にならなかった。

精度検証：球体の一様磁化
以下に、ｘ軸方向にＢ_ｏｘ＝０．１[Ｔ]の一様な外場においた磁性球（帯磁率が９９９）の磁化
を計算した結果を示す。電気伝導度はゼロとした。球を構成する粒子(原子) 数は４００９個とし、ｆｃｃ構造に配置した。厳密解は
である。この場合、図４、図５、図６、図７、図８に示される計算結果と同様な計算結果が得られた。

さらに、ルービング(looping) パターンを確認した結果を図１０、図１１、図１２に示す。従来の磁気モーメント法等に見られる磁場の「巻き」は見られない。そして、粒子の配列においても、規則正しい配列（ｆｃｃ、ｂｃｃ、ｓｃ）ならば、結果は同じである事が観れる。

図１０は、永久磁石（１．０Ｔ）による磁化の計算結果を示す図である。永久磁石はｂｃｃ、磁性体はｆｃｃ格子上に配列した。
図１１は、永久磁石（１．０Ｔ）による磁化の計算結果を示す図である。永久磁石はｆｃｃ、磁性体はｂｃｃ格子上に配列した。
図１２は、永久磁石（１．０Ｔ）による磁化の計算結果を示す図である。永久磁石はｆｃｃ、磁性体はｓｃ格子上に配列した。原子半径をｆｃｃ格子のものを使ったので磁性体内部の磁場が小さくでた。ｓｃ格子を充填する原子半径にすれば正解を得る。

Legendre関数の級数表示
以下にLegendre関数と倍Legendre関数の特別な場合を列挙する。（ｘ：＝ｃｏｓθ_ｊ）

粒子に働く力の導出
粒子（電子と核）と、粒子磁場間の相互作用を記述するLagrangeanは
である。ベクトルポテンシャル
と電磁場の関係は、
である。
を使い、以下の運動方程式を得る。
原子Ｍ＝Ｍ_ｐ＋ｍ_ｅに対する運動方程式は、
となる。右辺：Lorentz力を平均電流密度
、電気モーメントを
で置き換えると
Maxwellの方程式とＱ−２４から
最終的に原子に働く力は、電気分極を無視し、
となる。

粒子に働く力は、直接Maxwellの応力テンソルからも求まる。
球に局在した電流が作るモーメント
に作用する力は、最低の次数では
である。従い、ポテンシャルエネルギーは、電気分極によるポテンシャルエネルギーを合わせて、
となるが、渦電流や真電流を包含しない。

結局、
を計算すれば磁性球に働く磁力が求まる。

磁化による電流
と渦電流
を完全に切り分けることは恐らく出来ないが、それぞれ、

微係数の計算
以下に必要な微分係数を書き出す。
等からｒ_ｐｊ＞ａのみ列挙すると（なお、
はｓｉｎθに比例するから、θ＝０で発散は無い）、

磁気モーメントが作る磁化電流の計算
ｚ軸に平行に磁化したとき
今一度、ｚ軸に磁化したときの磁場を書く
磁化電流を
と書けば、磁化により生じる磁力
は
である。先ず、
を求める。

ｘ軸に平行に磁化したとき
ｘ軸に磁化したときの磁場を書く：
磁化電流
を求める。

ｙ軸に平行に磁化したとき
ｙ軸に磁化したときの磁場を書く：
磁化電流
を求める。

誘導磁界が作る渦電流の計算
誘導磁場は、Ｓ−３から
ここでＧ（ｒ）はグリーン関数(伝播関数) である。渦電流は

粒子に働く力の計算
ｍ_ｘ、ｍ_ｙ、ｍ_ｚと渦電流による寄与を加えれば、全電流が得られる。
座標
に粒子に働く力は
最後に
の微分から出てくる項を付加する必要がある。これは磁場の誤差
から生じる見かけの力であり、運動と磁場の全エネルギーを保存させる。

以上、実施の形態に係る解析装置１００の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。

実施の形態では、数値演算部１２０において粒子の位置と速度の両方を演算する場合について説明したが、これに限られない。例えば、数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法もあり、本実施の形態に係る技術的思想をそのような手法に適用してもよい。

１００解析装置、１０２入力装置、１０４ディスプレイ、１０８粒子系取得部、１１０磁気モーメント付与部、１１４粒子データ保持部、１１８表示制御部、１２０数値演算部、１２２力演算部、１２４粒子状態演算部、１２６状態更新部、１２８終了条件判定部、１３２磁場演算部。

Claims

仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する磁気モーメント付与部と、
前記磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部と、
前記数値演算部における演算結果を使用して、粒子系に関連する磁気的な物理量を演算する磁気的物理量演算部と、を備えることを特徴とする解析装置。
粒子系は繰り込み群分子動力学法を使用して繰り込まれた粒子系であることを特徴とする請求項１に記載の解析装置。
前記数値演算部における支配方程式は磁気モーメントに依存する項を有することを特徴とする請求項１または２に記載の解析装置。
前記磁気的物理量演算部は、粒子系が生成する磁場または磁束密度を演算することを特徴とする請求項１から３のいずれかに記載の解析装置。
仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する機能と、
磁気モーメントが付与された粒子を含む粒子系の各粒子の運動を支配する支配方程式を数値的に演算する機能と、
演算結果を使用して、粒子系に関連する磁気的な物理量を演算する機能と、をコンピュータに実現させることを特徴とするコンピュータプログラム。