JP6618417B2 - 解析装置および解析方法 - Google Patents

Description

本発明は、粒子系を解析する解析装置および解析方法に関する。

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

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

古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法に基づくシミュレーションが知られている（例えば、特許文献２参照）。粒子法は静的な現象だけでなく流れなどの動的な現象をも取り扱えるので、主に静的な現象を解析対象とする上述の有限要素法などに代わるシミュレーション手法として注目されている。例えば、各粒子に磁気モーメントを付与し、各粒子の球対称性に基づく厳密解を利用して磁気的な物理量を演算することで、比較的精度の高いシミュレーション結果を高速に得ることのできる磁気ビーズ法が提案されている（例えば、特許文献３参照）。

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

Erik Bitzek他、「Structural Relaxation Made Simple」、Physical Review Letters、２００６年１０月２７日、９７ L.D.Landau、E.M.Lifshitz、L.P.Litaevskii著、「Electrodynamics of Continuous Media 2nd ed」、１９８４年１月１日、２０５ページ 桑原敏彦、武田毅著、「三次元ラプラス問題のための解析的積分を用いた境界要素法」、電気学会論文誌Ａ １０６巻１１号、１９８６年、３４ページ

上述の磁気ビーズ法では、対象物を複数個の要素に分割し、各要素を球形状の粒子とみなして球対称性に基づく厳密解を利用して演算している。このとき、対象物を複数の粒子で充填させようとすると粒子間に隙間が生じるため、その隙間や球形状の粒子に起因する誤差が生じうる。計算精度を高めるためには、対象物を立方体などのボクセルで隙間なく分割し、ボクセルの形状に応じた厳密解に基づく面積分を実行する必要があるが、面積分を高精度で計算しようとすると計算時間が膨大となる。

本発明はこうした課題に鑑みてなされたものであり、その目的は、磁気ビーズ法による演算処理の精度向上および高速化を両立させた解析技術を提供することにある。

上記課題を解決するために、本発明のある態様の解析装置は、仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する磁気モーメント付与部と、磁気モーメントが付与される基準粒子が作る磁場の数値解を基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、基準粒子に付与される磁気モーメントのベクトル量と磁場の数値解のベクトル量との対応関係を複数の格子点ごとに保持する基準データ保持部と、基準データ保持部に保持される対応関係を用いて、磁気モーメント付与部によって磁気モーメントが付与される各粒子が作る磁場を演算する磁場演算部と、磁場演算部における演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、を備える。

本発明の別の態様は、解析方法である。この方法は、仮想空間内に定義される粒子系の粒子に磁気モーメントを付与し、磁気モーメントが付与される基準粒子が作る磁場の数値解を基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、基準粒子に付与される磁気モーメントのベクトル量と磁場の数値解のベクトル量との対応関係を複数の格子点ごとに保持し、あらかじめ算出して保持される対応関係を用いて、磁気モーメントが付与される各粒子が作る磁場を演算し、各粒子が作る磁場の演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算する。

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

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

実施の形態に係る解析装置の機能および構成を示すブロック図である。 図１の粒子データ保持部の一例を示すデータ構造図である。 図１の解析装置における一連の処理の一例を示すフローチャートである。 ｎ＝１次までの展開に対応する計算結果を示す図である。 ｎ＝２次までの展開に対応する計算結果を示す図である。 ｎ＝３次までの展開に対応する計算結果を示す図である。 ｎ＝４次までの展開に対応する計算結果を示す図である。 球面調和関数の次数と計算値との関係を示す図である。 ヒステリシスカーブを示す図である。 ビーズにより構成される導体球の解析モデルを示す図である。 導体球表面での時間変化する磁界の計算値を示すグラフである。 計算値の位相誤差と振幅誤差を示すグラフである。 導体球表面での磁界、電流および電磁力の計算値を示すグラフである。 基準粒子の周りに設定される近傍領域および遠方領域を模式的に示す図である。 近接する複数の格子点の磁場から観測点の磁場を算出する方法を模式的に示す図である。 近傍領域の任意の座標にて磁場を演算するために必要となる格子点を模式的に示す図である。 図１７（ａ）、（ｂ）、（ｃ）は、立方体表面の面積分を実行する方法の一例を模式的に示す図である。 微小三角形の面積分を実行する方法を模式的に示す図である。 実施の形態に係る解析装置の機能および構成を示すブロック図である。 図２０の解析装置における一連の処理の一例を示すフローチャートである。 図２０のＳ３６の処理の詳細を示すフローチャートである。 面積分による磁場の厳密解と双極子近似による磁場の数値解を比較したグラフである。 図２３（ａ）、（ｂ）は、本実施の形態に係る計算結果を示す図である。 本実施の形態に係る計算結果を示すグラフである。

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

（第１の実施の形態）
古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（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) は
を考えれば、
と成る。即ちｍ_ｖを十分小さく取れば、解くべき本来の方程式
の解に一致する。

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

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

精度検証：球体の一様磁化
図１０は、ビーズにより構成される導体球の解析モデルを示す図である。精度検証では、ｘ軸方向にＨ_{ｅｘｔ，ｘ}＝Ｈ_Ｏｓｉｎ（２πｆｔ）となる空間一様な外部磁界を印加した。ここで、Ｈ_Ｏ＝７．９５８×１０^４［Ａ／ｍ］、ｆ＝２００［Ｈｚ］、導体球の半径Ａ＝１０［ｍｍ］、電気電導率σ＝５９×１０^６［Ｓ／ｍ］とした。導体球を構成する粒子（原子）数は、６０９９個とし、ｆｃｃ構造に配置した。

図１１は、導体球表面での時間変化する磁界の計算値を示すグラフであり、導体表面近傍に位置する粒子の重心位置（ｒｘ＝９．８９［ｍｍ］、ｒｙ＝ｒｚ＝０）上の磁界Ｈｘの時間変化を示す。なお、本図では、シミュレーションによる計算結果とともに、外部磁界および厳密解による値も示している。誘導磁界の影響により、導体内部の磁界が外部磁界に対して位相が遅れて変化する様子を再現することができており、計算結果と厳密解の傾向は良く一致することが示された。

図１２は、計算値の位相誤差と振幅誤差を示すグラフであり、ｘ軸上に位置する粒子の重心位置における磁界の位相誤差ε_{ｐｈａｓｅ}と振幅誤差ε_{ａｍｐｌｉｔｕｄｅ}を示す。位相誤差ε_{ｐｈａｓｅ}および振幅誤差ε_{ａｍｐｌｉｔｕｄｅ}は、
と定義した。ここで、θ_{ｅｘａｃｔ}は厳密解の位相、θ_ＭＢＭは解析結果の位相である。また、Ｈ_{ｅｘａcｔ}は厳密解により得られた磁界、Ｈ_ＭＢＭは解析結果の磁界である。なお位相誤差は、Ｈ_{ｅｘａcｔ}が最大となる位相θ_{ｅｘａｃｔ}と、Ｈ_ＭＢＭが最大となる位相θ_ＭＢＭとを比較した。図１２より、ｘ軸上における解析結果と厳密解の磁界の位相誤差は最大で５．８６％、振幅誤差は６．３８％となり、厳密解と比較して精度の高い解析結果が得られることが示された。

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

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

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

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

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

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

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

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

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

誘導磁界が作る渦電流の計算
導体内部の誘導磁場は、Ｑ−２１より、
となる。渦電流は
なお、上式Ｒ−５２−２は、Ｒ−１２において、球面調和関数の展開次数ｎ＝０とした場合と同じである。

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

（第２の実施の形態）
上述の第１の実施の形態では、Ｑ−２１の右辺第１項に示される誘導磁化と、右辺第２項に示される磁気モーメントの相互作用とにより表される誘導磁場の式（Ｑ−２１）を用いて動磁場解析を行う場合を示した。本実施の形態では、Ｑ−２１の式を用いる代わりに、ベクトルポテンシャルの微分方程式から導出されるベクトルポテンシャルの厳密解を利用して誘導磁場を定式化する。以下、第１の実施の形態との相違点を中心に述べる。

ベクトルポテンシャルＡについて、以下の微分方程式が成立する。
上記微分方程式の解は、
となる。ここで、
として、各粒子に分割すれば、各粒子の位置ｒ_ｉでのベクトルポテンシャルは、
となる。ここで、ｖ_ｊは各粒子の体積である。なお、ｊの和に際し、ベクトル
が粒子表面を横切る場合には和を取らない。また、粒子径ａが満たすべき条件は、
である。

このとき、各対象物を粒子に分割することにより、ベクトルポテンシャルＡの並進対称性を失うおそれがある。そこで、本発明者は、下記Ｓ−４に示されるゲージ変換を導入することにより、ベクトルポテンシャルＡの並進対称性を回復できることを見いだした。なお、このゲージ変換を本明細書において「ビーズゲージ」ともいう。
ここで、ベクトルｇ_ｉは、計算対象とする粒子群（粒子数：Ｎ_ｉ）に関する局所重心ベクトルであり、
で表される。したがって、ベクトルｇ_ｉは、粒子群の重心位置を表している。また、Ｓ−４に示すビーズゲージは、磁束密度Ｂ_ｊと局所重心ベクトルｇ_ｉの外積を用いて記述されると言える。

上記のビーズゲージについて、磁束密度ＢとベクトルポテンシャルＡの関係式
を求めれば、
となる。したがって、粒子数Ｎ_ｉが十分に大きければ、並進対称性が満たされることがわかる。

ここで、Ｓ−４に示すゲージ変換から、Ｓ−３の右辺に含まれるベクトルポテンシャルの時間微分
を磁束密度Ｂの時間微分を用いて以下のように表すことができる。
さらに、ベクトルｄ_ｊｉを、
とすると、ベクトルポテンシャルの時間微分は、
と表される。したがって、Ｓ−３にＳ−９を代入して積分を実行すれば、ベクトルポテンシャルＡの厳密解を磁束密度Ｂの時間微分を用いて記述できる。

ベクトルポテンシャルＡの厳密解を導出するため、まず、磁束密度Ｂのｚ成分の時間微分である、
のみが作用する場合の積分を実行する。これにより、ベクトルポテンシャルＡは、
となる。磁束密度Ｂのｘ成分およびｙ成分の時間微分である、
の寄与も加えれば、
が得られる。

解析対象とする導体内部を記述する誘導磁場Ｈは、
の関係式より、
となる。ここで、
と置換すると、誘導磁場Ｈの式として、
が得られる。

磁場演算部１２１は、粒子系により生成される誘導磁場を、Ｓ−１３に示す誘導磁場の式を用いて演算する。誘導磁場の時間変化値を得るための計算は、上述の第１の実施の形態で示したものと同様の方法を用いればよい。磁場演算部１２１は、得られた誘導磁場から電流密度を計算してもよい。また、力演算部１２２は、得られた磁場および電流の値から、各粒子に作用する電磁力を計算してもよい。

本実施の形態において、Ｓ−５に示す局所重心ベクトルｇ_ｉは、互いに同電位となる粒子群に対して導入され、電気的に互いに絶縁される複数の粒子群に対しては、それぞれ別の局所重心ベクトルが導入される。例えば、解析対象とする粒子系において、互いが電気的に絶縁される複数の部分を有する場合、それぞれの部分を構成する粒子群に対して重心位置の異なる局所重心ベクトルが適用される。また、解析対象とする物体が時間経過とともに変形して複数の物体に分裂し、分裂した部分同士が互いに絶縁した状態となる場合には、それぞれの部分を構成する粒子群に対して重心位置の異なる局所重心ベクトルを適用する。

したがって、解析対象の粒子系に、電気的に互いが絶縁される第１粒子群および第２粒子群を含む複数の粒子群が配置される場合、磁場演算部１２１は、それぞれの粒子群に対応する誘導磁場の式を用いて解析をする。例えば、第１粒子群における誘導磁場は、第１粒子群の重心位置を表す第１局所重心ベクトルを用いたゲージ変換により導出される第１の誘導磁場の式を用いて演算される。また、第２粒子群における誘導磁場は、第２粒子群の重心位置を表す第２局所重心ベクトルを用いたゲージ変換により導出される第２の誘導磁場の式を用いて演算される。より具体的には、Ｓ−１３に示す誘導磁場の式において、右辺第１項に含まれるｄ_ｉｊがそれぞれの粒子群での計算において異なる値をとることとなる。

本実施の形態に係る誘導磁場の式（Ｓ−１３）を用いた解析手法について精度検証を行った。第１の実施の形態と同様に、図１０に示す導体球の解析モデルを用い、ｘ軸方向にＨ_{ｅｘｔ，ｘ}＝Ｈ_Ｏｓｉｎ（２πｆｔ）となる空間一様な外部磁界を印加した。計算条件は、Ｈ_Ｏ＝７．９４８×１０^４［Ａ／ｍ］、ｆ＝１００［Ｈｚ］、導体球の半径Ａ＝１０［ｍｍ］、電気電導率σ＝５９×１０^６［Ｓ／ｍ］とした。導体球を構成する粒子（原子）数は、６０９９個とし、ｆｃｃ構造に配置した。

図１３は、導体球表面での磁界、電流密度および電磁力の計算値を示すグラフである。本図は、導体表面近傍に位置する粒子の重心位置における、（ａ）磁界Ｈｘ、（ｂ）電流密度ｊｚ、（ｃ）電磁力Ｆｙの時間変化を示す。本図では、厳密解の値を実線で、計算値を破線で示している。図示されるように、本実施の形態に係る解析手法を用いると、計算結果と厳密解がよく一致することがわかった。

（第３の実施の形態）
上述の第１の実施の形態では、対象物を粒子に分割し、各粒子の球対称性に基づく厳密解を利用することで、高速かつ精度の優れた演算結果を得ている。しかしながら、連続体である対象物を粒子で分割すると粒子間に隙間が生じるため、その隙間に起因する計算誤差がどうしても生じてしまう。特に、粒子間の距離が近い場合の粒子近傍の磁場を計算する場合に誤差の影響が大きくなる。計算精度を高めるためには、対象物を立方体やボロノイ多面体といったボクセルに分割し、これらのボクセルにより隙間なく対象物を分割する必要がある。この場合、ボクセルの形状に応じた厳密解に基づく面積分を実行する必要があるため計算時間が膨大となる。

本実施の形態では、粒子近傍の磁場を演算する場合にボクセル形状に応じた厳密解を利用し、近似による誤差がほぼ無視できる粒子遠方の磁場を演算する場合に双極子近似を利用する。このとき、ボクセル形状に応じた厳密解の計算を都度実行するのではなく、あらかじめボクセルに付与される磁気モーメントのベクトル量とそのボクセルが作る磁場のベクトル量との対応関係を厳密解に基づいて算出しておき、その対応関係を参照することで各粒子の近傍における磁場を演算する。したがって、本実施の形態によれば、負荷の高い面積分の都度演算を回避しつつ、あらかじめ算出した厳密解に基づく対応関係を用いて磁場を高精度かつ高速に演算することができる。以下、本実施の形態について、上述の実施の形態との相違点を中心に説明する。

仮想空間内に定義される観測点ｒ_ｉの磁場は、各微小体積ｄｖ’の磁化Ｍ（ｒ’）が作る磁場を全空間にわたって積分することにより得られ、下記式（Ｔ−１）により表される。式（Ｔ−１）は、体積分の項と面積分の項に分けることによって、式（Ｔ−２）に変形できる。なお、Ｖは領域体積、Ｓは表面積、ｎ_Ｓは表面の法線ベクトルである。

磁気ビーズ法では、対象空間をＮ個の粒子（ビーズ）に分割し、各ビーズ内部の磁化ベクトルＭ_ｊが一定であるとし、式（Ｔ−１）を双極子近似して得られる下記式（Ｔ−３）により座標ｒ_ｉの磁場ベクトルＨ（ｒ_ｉ）を計算する。ここで、ｒ_ｉｊ＝ｒ_ｉ−ｒ_ｊ，ｎ_ｉｊ＝ｒ_ｉｊ／｜ｒ_ｉｊ｜であり、ΔＶ_ｊはビーズ体積であり、γ_ｊはビーズの充填率である。

一方、磁気モーメント法では、対象をＮ個の要素（ボクセル）に分割し、要素内部の磁化ベクトルＭ_ｊは一定であるとし、式（Ｔ−２）を変形して得られる下記式（Ｔ−４）により座標ｒ_ｉの磁場ベクトルＨ（ｒ_ｉ）を計算する。ここで、ｄｓ_ｊは各要素の微小表面を意味し、ボクセルが正六面体（立方体）であれば、６つの表面について面積分を実行する必要がある。

磁気ビーズ法と磁気モーメント法を比較すると、計算精度については磁気モーメント法が優れており、計算速度については磁気ビーズ法が優れている。磁気ビーズ法は、連続体をＮ個の球形状の粒子に分割し、各粒子が作る磁場を双極子近似して求めているため、各粒子近傍の磁場の計算精度が低い。その一方で、各粒子から離れた遠方領域の磁場については双極子近似による誤差が無視できるために計算精度が高い。磁気モーメント法は、各要素間に隙間が生じないように連続体を多面体要素に分割するため、磁気ビーズ法において生じるような近似誤差がなく、計算精度が高い。

磁気モーメント法では、上記式（Ｔ−４）に示される面積分を高精度で実行する必要があり、数値積分法や厳密解での計算が必須である。このとき、多面体要素が有する面数ｎ（立方体であればｎ＝６）の回数だけ面積分を実行する必要があり、積分回数が多くなることによって計算時間が膨大となる。一方、磁気ビーズ法では、ビーズの球形状に基づく厳密解を利用すれば、上記式（Ｔ−３）に示されるように面積分が含まれないため、計算時間が少なくて済む。

そこで、本実施の形態では、計算精度と計算速度を両立させるために、双極子近似が成立する遠方領域内の磁場の計算では、磁気ビーズ法に基づく式（Ｔ−３）を利用し、双極子近似が成立しない近傍領域内の磁場の計算では、磁気モーメント法に基づく式（Ｔ−４）を利用する。さらに、式（Ｔ−４）に基づく演算処理を高速化させるために、面積分の演算結果をあらかじめ算出し、磁化ベクトルＭ_ｊと、磁化ベクトルＭ_ｊがその周りに作る磁場Ｈ_ｊとの対応関係を定める磁場テーブルを用意する。シミュレーション実行時の演算処理では、面積分を都度実行する代わりにあらかじめ用意した磁場テーブルを参照して磁場を演算する。

図１４は、基準粒子Ｖ_ｊの周りに設定される近傍領域Ｃ_Ａおよび遠方領域Ｃ_Ｂを模式的に示す図である。図１４は、対象物が立方体形状の要素（ボクセル）で分割される場合を示しており、各要素の頂点に格子点Ｇ_ｋが設定される。各立方体要素の一辺の長さは、粒子の半径ａの２倍（２ａ）である。基準粒子（基準要素）Ｖ_ｊを中心として直交座標系（ξ，η，ζ）が設定され、立方体の中心に各要素の座標ｒ_ｊが設定される。基準粒子Ｖ_ｊの周りには半径ａ_ｓの仮想的な球が定められ、この球の内部が近傍領域Ｃ_Ａとなり、外側が遠方領域Ｃ_Ｂとなる。近傍領域Ｃ_Ａを規定する球の半径ａ_ｓは、双極子近似が成立する距離に応じて決定され、例えば、粒子径ａの４倍以上、好ましくは粒子径ａの６倍以上の値が半径ａ_ｓとして設定される。図示する例では、球の半径ａ_ｓを粒子径ａの６倍（つまり、ａ_ｓ＝６ａ）としている。

図１４において、基準粒子Ｖ_ｊが周りに作る磁場のうち、遠方領域Ｃ_Ｂの任意の座標ｒ_ｉにおける磁場Ｈ_Ｂ（ｒ_ｉ）は、式（Ｔ−３）に基づいて算出される。一方、基準粒子Ｖ_ｊが周りに作る磁場のうち、近傍領域Ｃ_Ａの内部の任意の座標ｒ_ｉにおける磁場Ｈ_Ａ（ｒ_ｉ）は、座標ｒ_ｉに近接する複数の格子点Ｇ_１〜Ｇ_４（実際には、後述する図１５に示す８個の格子点）の磁場を平均化することにより算出される。各格子点の磁場は、式（Ｔ−４）に基づいて算出されるものであり、格子点ごとにあらかじめ算出される磁場テーブルの値を用いて演算される。

図１５は、近接する複数の格子点の磁場Ｈ_１〜Ｈ_８から観測点ｒ_ｉの磁場Ｈ_Ａを算出する方法を模式的に示す図である。近傍領域Ｃ_Ａ内の磁場Ｈ_Ａは、磁場Ｈ_Ａの演算対象となる座標（ξ，η，ζ）を囲む８個の格子点における磁場Ｈ_ｋ（ｋ＝１〜８）を平均化することにより算出され、下記式（Ｔ−５）で表される。なお、Ｎ_ｋ（ξ，η，ζ）は、式（Ｔ−６）で表される。

ここで、基準粒子Ｖ_ｊが各格子点Ｇ_ｋ（位置ｒ_ｉ）に作る磁場Ｈ_ｋは、式（Ｔ−４）から下記式（Ｔ−７）で表すことができる。

式（Ｔ−７）で表される磁場Ｈ_ｋは、面積分を利用した厳密解である。ここで、式（Ｔ−７）を変形し、磁場Ｈ_ｋと磁化ベクトルＭ_ｊの対応関係を定めるテンソルＴ_ｋを用いて下記式（Ｔ−８）のように表すことができれば、任意の磁化ベクトルＭ_ｊを有する粒子が各格子点Ｇ_ｋに作る磁場Ｈ_ｋを容易に算出できるようになる。式（Ｔ−７）から、テンソルＴ_ｋの各成分は、下記式（Ｔ−９）で表すことができる。なお、ｎ_ｘ，ｎ_ｙ，ｎ_ｚはそれぞれｘ方向、ｙ方向、ｚ方向の単位ベクトルであり、ｘ_ｉｊ，ｙ_ｉｊ，ｚ_ｉｊはそれぞれ、基準粒子Ｖ_ｊから観測点（格子点）までの距離ｒ_ｉｊのｘ成分、ｙ成分、ｚ成分である。

そこで、本実施の形態では、式（Ｔ−９）で表されるテンソルＴ_ｋの各成分の数値解をあらかじめ算出し、その数値解と磁化ベクトルＭ_ｊの積から各格子点の磁場Ｈ_ｋを算出し、さらに、近傍領域Ｃ_Ａ内の磁場Ｈ_Ａを算出できるようにする。つまり、下記式（Ｔ−１０）を用いて近傍領域Ｃ_Ａ内の磁場Ｈ_Ａを演算する。あらかじめテンソルＴ_ｋの各成分の数値を求めておけば、式（Ｔ−１０）の計算に面積分が含まれないため、面積分を都度実行する場合よりも非常に短い時間で精度の高い計算結果を得ることができる。

図１６は、近傍領域Ｃ_Ａの任意の座標ｒ_ｉにて磁場Ｈ_Ａを演算するために必要となる格子点Ｇ_ｋを模式的に示す図である。仮想的な半径ａ_ｓの球の内部の全ての座標ｒ_ｉで磁場Ｈ_Ａを計算できるようにするためには、近傍領域Ｃ_Ａを囲うことのできるボクセル群の範囲内（太線の範囲内）の全ての格子点Ｇ_ｋについて、磁場Ｈ_Ａと磁化ベクトルＭ_ｊの対応関係を定めるテンソルＴ_ｋを算出しておけばよい。例えば、近傍領域Ｃ_Ａを定める範囲がビーズ径ａの６倍であれば、事前計算が必要となる格子点の数は５００個未満である。この数は、仮想空間内に定義される粒子系に含まれる全粒子数よりも少ない。そのため、全粒子が作る磁場を面積分を用いて１回計算するよりも短い時間で、磁気的な対応関係を定めるテンソルＴ_ｋを導出することができる。つまり、磁気テーブルを算出するための演算負荷はそれほど高くなく、全計算時間に占める磁気テーブルの算出時間の割合は小さい。

つづいて、磁気的な対応関係を定めるテンソルＴ_ｋの具体的な計算方法の一例を示す。テンソルＴ_ｋの各成分の数値解の演算方法は、式（Ｔ−９）の計算が実行できれば特に問わないが、例えば、非特許文献３に開示される方法を用いることができる。図１７（ａ）、（ｂ）、（ｃ）は、立方体表面の面積分を実行する方法の一例を模式的に示す図である。まず、面積分の対象となる六つの面の一つを取り出し、取り出した正方形を二つに分割した微小な三角形を考える。このようにして六つの面を１２個の微小な三角形に分割し、微小三角形ごとに面積分を実行する。

図１８は、微小三角形の面積分を実行する方法を模式的に示す図である。図１８に示すローカル座標系（ｘ’，ｙ’，ｚ’）において、面積分の対象となる微小三角形△１２３がｘ’ｙ’平面内に配置され、微小三角形△１２３により作られる磁場ΔＨ_ｋの座標ｒ_ｉがｚ’軸上に配置されるように座標変換を施す。このとき、微小三角形△１２３により作られる磁場ΔＨ_ｋについて、ΔＨ_ｋ＝ΔＴ_ｋＭ_ｊと表すと、微小三角形Δ１２３に起因するテンソルΔＴ_ｋは、下記式（Ｔ−１１）で表される。式（Ｔ−１１）に含まれる各変数は、式（Ｔ−１２）で定義される。

このようにして微小三角形Δ１２３に起因して得られるテンソルΔＴ_ｋを基準粒子Ｖ_ｊの全ての面について演算して足し合わせることで、磁気的な対応関係を定めるテンソルＴ_ｋを導出できる。磁気的な対応関係を定めるテンソルＴ_ｋは、近傍領域Ｃ_Ａの任意の座標ｒ_ｉで磁場Ｈ_Ａを演算するために必要となる全て格子点Ｇ_ｋについて計算され、その数値解が磁気テーブルとして保持される。

図１９は、実施の形態に係る解析装置２００の機能および構成を示すブロック図である。本実施の形態に係る解析装置２００は、基準データ保持部２３０をさらに備える点で図１に示す解析装置１００と相違する。基準データ保持部２３０は、基準粒子に付与される磁気モーメントのベクトル量と、基準粒子が周りに作る磁場の厳密解のベクトル量との対応関係を保持する。つまり、基準データ保持部２３０は、磁気的な対応関係を定めるテンソルＴ_ｋの数値解が格子点Ｇ_ｋごとに定められた磁気テーブルを保持する。

磁場演算部１２１は、粒子系取得部１０８および磁気モーメント付与部１１０によって定義された粒子系に応じて、上述の磁気テーブルを導出するための演算処理を実行する。磁場演算部１２１は、粒子系の定義に応じて近傍領域Ｃ_Ａを設定し、磁気的な対応関係を定めるテンソルＴ_ｋの算出が必要となる複数の格子点Ｇ_ｋを特定する。磁場演算部１２１は、原点に配置した基準粒子Ｖ_ｊが各格子点Ｇ_ｋの位置に作る磁場Ｈ_ｋに関連する演算を実行し、基準粒子Ｖ_ｊと各格子点Ｇ_ｋの配置関係を示す値とともにテンソルＴ_ｋの数値解を基準データ保持部２３０に記憶させる。

磁場演算部１２１は、粒子ｊが座標ｒ_ｉに作る磁場Ｈ_ｊ（ｒ_ｉ）を演算する場合、座標ｒ_ｉが粒子ｊの近傍領域内であれば磁気テーブルを参照して磁場を演算し、座標ｒ_ｉが粒子ｊの遠方領域（近傍領域外）であれば双極子近似による式を用いて磁場を演算する。磁場演算部１２１は、粒子ｊと座標ｒ_ｉの距離が所定値を超える場合、つまり、座標ｒ_ｉから大きく離れた位置にある粒子ｊについては、磁場の演算対象とせずに無視をしてもよい。例えば、遠方領域のうち、近傍領域に近い中間領域については、双極子近似による磁場演算を実行し、中間領域よりも遠い領域については磁場演算を実行しなくてもよい。中間領域の範囲は、例えば、粒子径ａの１０倍程度である。

磁場演算部１２１は、磁気テーブルを参照して磁場を演算する場合、座標ｒ_ｉの位置を図１６に示す直交座標系（ξ，η，ζ）に変換して磁場を算出し、その後、元の座標系に戻すことによって座標ｒ_ｉの磁場を演算する。このとき、粒子系の座標系と磁気テーブルの直交座標系（ξ，η，ζ）との間の座標変換は平行移動のみで済むため、座標変換に関する処理負荷を小さくできる。変形例においては、座標変換時に拡大縮小や回転などの座標変換処理を施してもよい。

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

図２１は、図２０のＳ３６の処理の詳細を示すフローチャートである。磁場演算部１２１は、変数ｉに初期値ｉ＝１を設定する（Ｓ５０）。磁場演算部１２１は、基準データ保持部２３０に保持される磁気テーブルを参照して近傍粒子が座標ｒ_ｉに作る磁場を算出し（Ｓ５２）、双極子近似式を用いて遠方粒子が座標ｒ_ｉに作る磁場を算出する（Ｓ５４）。磁場演算部１２１は、算出した磁場を足し合わせ、座標ｒ_ｉの粒子の磁化ベクトルＭ_ｉ（磁気モーメント）をＭ_ｉ＝χＨ_ｉの関係式を用いて演算する（Ｓ５６）。変数ｉに１を加え（Ｓ５８）、変数ｉが粒子数Ｎ以下であれば（Ｓ６０のＮ）、Ｓ５２〜Ｓ５８のステップを繰り返し、変数ｉが粒子数Ｎより大きければ（Ｓ６０のＹ）、Ｓ３６の処理を終了する。

図２２は、面積分による磁場の厳密解Ｂ１と双極子近似による磁場の数値解Ｂ２を比較したグラフである。図２２において、横軸はソースとなる粒子ｊと観測点ｉの距離ｒ_ｉｊを粒子径ａで規格化した値であり、縦軸はｘ軸方向の磁束密度Ｂｘの値である。図示されるように、粒子ｊと観測点ｉの距離ｒ_ｉｊが２倍となる場合に厳密解Ｂ１（実線）と近似解Ｂ２（破線）の差が大きくなり、距離ｒ_ｉｊが大きくなるにつれて差が小さくなり、距離ｒ_ｉｊが粒子径ａの６倍以上になると差がほとんど生じないことがわかる。このことから、近傍領域Ｃ_Ａとして粒子径ａの６倍程度の範囲を設定すればよく、６倍を超える範囲については双極子近似による誤差がほとんど無視できると言える。

図２３（ａ）、（ｂ）は、本実施の形態に係る計算結果を示す図である。本図は、ｘ方向に１０ｍｍ、ｙ方向およびｚ方向に５ｍｍが設定された直方体形状の磁性体の磁化現象の計算例を示す。外部磁場を一定の０．１［Ｔ］とし、磁化特性は線形（χ＝９９９）とした。図２３（ａ）は、比較例であり、全ての粒子に対して上記式（Ｔ−４）に基づく面積分を実行して磁場を演算した結果を示す。図２３（ｂ）は、本実施の形態に係る磁気テーブルを用いて磁場を演算した結果を示す。両者を比較したところ、計算結果がよく一致することが分かった。また、図２３（ａ）の例では１ステップあたりの計算時間が２．１５秒であるのに対し、図２３（ｂ）の例では１ステップあたりの計算時間が０．０３秒であった。したがって、本実施の形態に係る手法により、高い計算精度を維持したまま約７０倍の計算速度の高速化を実現できた。

図２４は、本実施の形態に係る計算結果を示し、上述の図２３に示す計算例のｙ軸上の磁束密度分布の絶対値｜Ｂ｜を示す。実線で示す磁束密度分布Ｂ３は、図２３（ｂ）の実施の形態に係る計算例に対応し、破線で磁束密度分布Ｂ４は、図２３（ａ）の比較例に係る計算例に対応する。図示されるように、両者の計算結果がよく一致することが分かる。

本実施の形態によれば、計算過程においてボクセル形状が変化しないという磁気ビーズ法の特徴をうまく利用して高い計算精度を維持しつつ計算速度の高速化を実現できる。有限要素法では、節点位置の変化によりボクセル形状が変化し、ボクセル表面の面要素の形状が変化するため、変化した面要素の形状に応じた面積分を実行する必要がある。一方、磁気ビーズ法では、ボクセルの形状が計算過程にわたって変わらないため、共通する基準粒子を用いた面積分の結果を流用できる。このように本実施の形態では、磁気ビーズ法の特徴を利用し、共通する基準粒子に基づく磁気テーブルをあらかじめ用意することで計算速度を向上させることができる。

なお変形例においては、粒子系にボクセル形状の異なる複数種類の粒子を定義してもよい。この場合、ボクセル形状の異なる複数種類の基準粒子ごとにあらかじめ磁気テーブルを用意しておき、任意の観測点ｒ_ｉの位置に磁場を作る粒子ｊのボクセル形状に応じて、異なる種類の磁気テーブルを参照して磁場を算出すればよい。これにより、複数種類のボクセル形状を利用する場合であっても上述の実施の形態と同様の効果を奏することができる。ボクセル形状は、多面体であればどのような形状であってもよく、正四面体、正八面体、正十二面体などであってもよい。また、対象空間をボロノイ分割することにより得られるボロノイ多面体であってもよい。

上述の実施の形態では、近傍粒子と遠方粒子とで使用する数式を切り替える場合を示した。変形例においては、遠方に存在する粒子（例えば、粒子径ａの１０倍より遠い粒子）については磁場演算をせず、磁場演算となる粒子について全て磁気テーブルを参照して演算してもよい。つまり、粒子径ａが１０倍となる範囲について磁気的な対応関係Ｔ_ｋを算出し、対応関係Ｔ_ｋに基づいて磁場演算を実行してもよい。

上述の実施の形態では、磁気テーブルにテンソルＴ_ｋが保持される場合を示したが、式（Ｔ−１０）に示される計算の実行に必要な形式であれば、その他の種類の数値解が磁気テーブルとして保持されてもよい。例えば、変形例においては、Ｈ_Ａ（ξ，η，ζ）＝Ｓ_ｋ（ξ，η，ζ）Ｍ_ｊと記述されるテンソルＳ_ｋの数値解を磁気テーブルとして保持してもよい。

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

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

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

Claims (9)

1. 仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する磁気モーメント付与部と、
   磁気モーメントが付与される基準粒子が作る磁場の数値解を前記基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、前記基準粒子に付与される磁気モーメントのベクトル量と前記磁場の数値解のベクトル量との対応関係を前記複数の格子点ごとに保持する基準データ保持部と、
   前記基準データ保持部に保持される前記対応関係を用いて、前記磁気モーメント付与部によって磁気モーメントが付与される各粒子が作る磁場を演算する磁場演算部と、
   前記磁場演算部における演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、を備えることを特徴とする解析装置。
2. 前記基準データ保持部は、前記基準粒子に付与される磁気モーメントのベクトル量と前記基準粒子が作る磁場の厳密解のベクトル量との対応関係を前記複数の格子点ごとに保持することを特徴とする請求項１に記載の解析装置。
3. 前記基準データ保持部は、前記基準粒子の近傍に設定される複数の格子点であって、前記磁気モーメント付与部によって磁気モーメントが付与される粒子数より少ない複数の格子点における前記対応関係を保持することを特徴とする請求項１または２に記載の解析装置。
4. 前記磁場演算部は、前記磁気モーメント付与部によって磁気モーメントが付与される各粒子が作る磁場のうち、各粒子の近傍領域内の磁場を前記基準データ保持部に保持される前記対応関係を用いて演算し、各粒子の近傍領域外の磁場を各粒子の磁気モーメントを双極子近似して演算することを特徴とする請求項１から３のいずれか一項に記載の解析装置。
5. 前記近傍領域は、各粒子の中心からの距離が粒子径の２倍より大きく、１０倍より小さい範囲となるように設定されることを特徴とする請求項４に記載の解析装置。
6. 前記磁場演算部は、各粒子の前記近傍領域内の観測点の磁場を演算する場合、前記観測点に近接する複数の格子点の磁場を前記対応関係から算出し、前記観測点に近接する複数の格子点について算出した磁場から前記観測点の磁場を演算することを特徴とする請求項４または５に記載の解析装置。
7. 前記磁気モーメント付与部は、粒子系に定義されるボクセル形状の異なる複数種類の粒子に磁気モーメントを付与し、
   前記基準データ保持部は、ボクセル形状の異なる複数種類の基準粒子ごとにあらかじめ算出した複数種類の前記対応関係を保持し、
   前記磁場演算部は、前記磁気モーメント付与部によって磁気モーメントが付与される各粒子のボクセル形状に応じた種類の前記対応関係を用いて、各粒子が作る磁場を演算することを特徴とする請求項１から６のいずれか一項に記載の解析装置。
8. 仮想空間内に定義される粒子系の粒子に磁気モーメントを付与し、
   磁気モーメントが付与される基準粒子が作る磁場の数値解を前記基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、前記基準粒子に付与される磁気モーメントのベクトル量と前記磁場の数値解のベクトル量との対応関係を前記複数の格子点ごとに保持し、
   あらかじめ算出して保持される前記対応関係を用いて、磁気モーメントが付与される各粒子が作る磁場を演算し、
   各粒子が作る磁場の演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算することを特徴とする解析方法。
9. 仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する機能と、
   磁気モーメントが付与される基準粒子が作る磁場の数値解を前記基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、前記基準粒子に付与される磁気モーメントのベクトル量と前記磁場の数値解のベクトル量との対応関係を前記複数の格子点ごとに保持する機能と、
   あらかじめ算出して保持される前記対応関係を用いて、磁気モーメントが付与される各粒子が作る磁場を演算する機能と、
   各粒子が作る磁場の演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算する機能と、をコンピュータに実現させることを特徴とするコンピュータプログラム。