(57)【特許請求の範囲】
【請求項１】
仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する磁気モーメント付与部と、
磁気モーメントが付与される基準粒子が作る磁場の数値解を前記基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、前記基準粒子に付与される磁気モーメントのベクトル量と前記磁場の数値解のベクトル量との対応関係を前記複数の格子点ごとに保持する基準データ保持部と、
前記基準データ保持部に保持される前記対応関係を用いて、前記磁気モーメント付与部によって磁気モーメントが付与される各粒子が作る磁場を演算する磁場演算部と、
前記磁場演算部における演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、を備えることを特徴とする解析装置。
【請求項２】
前記基準データ保持部は、前記基準粒子に付与される磁気モーメントのベクトル量と前記基準粒子が作る磁場の厳密解のベクトル量との対応関係を前記複数の格子点ごとに保持することを特徴とする請求項１に記載の解析装置。
【請求項３】
前記基準データ保持部は、前記基準粒子の近傍に設定される複数の格子点であって、前記磁気モーメント付与部によって磁気モーメントが付与される粒子数より少ない複数の格子点における前記対応関係を保持することを特徴とする請求項１または２に記載の解析装置。
【請求項４】
前記磁場演算部は、前記磁気モーメント付与部によって磁気モーメントが付与される各粒子が作る磁場のうち、各粒子の近傍領域内の磁場を前記基準データ保持部に保持される前記対応関係を用いて演算し、各粒子の近傍領域外の磁場を各粒子の磁気モーメントを双極子近似して演算することを特徴とする請求項１から３のいずれか一項に記載の解析装置。
【請求項５】
前記近傍領域は、各粒子の中心からの距離が粒子径の２倍より大きく、１０倍より小さい範囲となるように設定されることを特徴とする請求項４に記載の解析装置。
【請求項６】
前記磁場演算部は、各粒子の前記近傍領域内の観測点の磁場を演算する場合、前記観測点に近接する複数の格子点の磁場を前記対応関係から算出し、前記観測点に近接する複数の格子点について算出した磁場から前記観測点の磁場を演算することを特徴とする請求項４または５に記載の解析装置。
【請求項７】
前記磁気モーメント付与部は、粒子系に定義されるボクセル形状の異なる複数種類の粒子に磁気モーメントを付与し、
前記基準データ保持部は、ボクセル形状の異なる複数種類の基準粒子ごとにあらかじめ算出した複数種類の前記対応関係を保持し、
前記磁場演算部は、前記磁気モーメント付与部によって磁気モーメントが付与される各粒子のボクセル形状に応じた種類の前記対応関係を用いて、各粒子が作る磁場を演算することを特徴とする請求項１から６のいずれか一項に記載の解析装置。
【請求項８】
仮想空間内に定義される粒子系の粒子に磁気モーメントを付与し、
磁気モーメントが付与される基準粒子が作る磁場の数値解を前記基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、前記基準粒子に付与される磁気モーメントのベクトル量と前記磁場の数値解のベクトル量との対応関係を前記複数の格子点ごとに保持し、
あらかじめ算出して保持される前記対応関係を用いて、磁気モーメントが付与される各粒子が作る磁場を演算し、
各粒子が作る磁場の演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算することを特徴とする解析方法。
【請求項９】
仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する機能と、
磁気モーメントが付与される基準粒子が作る磁場の数値解を前記基準粒子の周りに設定される複数の格子点ごとにあらかじめ算出し、前記基準粒子に付与される磁気モーメントのベクトル量と前記磁場の数値解のベクトル量との対応関係を前記複数の格子点ごとに保持する機能と、
あらかじめ算出して保持される前記対応関係を用いて、磁気モーメントが付与される各粒子が作る磁場を演算する機能と、
各粒子が作る磁場の演算結果を用いて、各粒子の運動を支配する支配方程式を数値的に演算する機能と、をコンピュータに実現させることを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子系を解析する解析装置および解析方法に関する。
【背景技術】
【０００２】
近年、コンピュータの計算能力の向上に伴い、モータなどの電気機器の設計開発の現場では磁場解析を取り入れたシミュレーションがよく使用される。シミュレーションを使用すると、実際にプロトタイプを製作しなくてもある程度の評価が可能となるので、設計開発のスピードが向上しうる。
【０００３】
例えば特許文献１には、磁場解析を実行する演算処理装置を具えるモータ解析装置が記載されている。演算処理装置は、ユーザの操作に基づく外部指令に応じて、有限要素法による静磁場解析やマクスウェルの応力法によるトルク計算を実行する。有限要素法による静磁場解析のためにメッシュ分割が行われる。メッシュ分割は、コア領域及びハウジング領域、さらには外部空気層領域についても行われる。磁場解析の他の手法は差分法および磁気モーメント法を含む。
【０００４】
古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法に基づくシミュレーションが知られている（例えば、特許文献２参照）。粒子法は静的な現象だけでなく流れなどの動的な現象をも取り扱えるので、主に静的な現象を解析対象とする上述の有限要素法などに代わるシミュレーション手法として注目されている。例えば、各粒子に磁気モーメントを付与し、各粒子の球対称性に基づく厳密解を利用して磁気的な物理量を演算することで、比較的精度の高いシミュレーション結果を高速に得ることのできる磁気ビーズ法が提案されている（例えば、特許文献３参照）。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開平１１−１４６６８８号公報
【特許文献２】特開２００９−３７３３４号公報
【特許文献３】特開２０１５−１１１４０１号公報
【非特許文献】
【０００６】
【非特許文献１】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）の値から力を算出する。力演算部１２２は、近接粒子がｉ番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、ｉ番目の粒子に働く力を算出する。
【００３３】
粒子状態演算部１２４は粒子データ保持部１１４に保持される粒子系のデータを参照し、粒子系の各粒子について、離散化された粒子の運動方程式に力演算部１２２によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。
【００３４】
粒子状態演算部１２４は、力演算部１２２によって演算された力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき所定の微小な時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって演算された力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。
【００３５】
粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。
【００３６】
状態更新部１２６は、粒子データ保持部１１４に保持される粒子系の各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。
【００３７】
終了条件判定部１２８は、数値演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、外部から終了の指示を受け付けたことや、粒子系が定常状態に達したことである。終了条件判定部１２８は、終了条件が満たされる場合、数値演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を力演算部１２２に戻す。すると力演算部１２２は、状態更新部１２６によって更新された粒子の位置、速度で再び力を演算する。
【００３８】
表示制御部１１８は、粒子データ保持部１１４に保持されるデータが表す粒子系の各粒子の位置、速度、磁気モーメントに基づき、ディスプレイ１０４に粒子系の時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。
【００３９】
図２は、粒子データ保持部１１４の一例を示すデータ構造図である。粒子データ保持部１１４は、粒子ＩＤと、粒子の位置と、粒子の速度と、粒子の磁気モーメントと、を対応付けて保持する。
【００４０】
上述の実施の形態において、保持部の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。
【００４１】
以上の構成による解析装置１００の動作を説明する。
図３は、解析装置１００における一連の処理の一例を示すフローチャートである。粒子系取得部１０８は、ＲＭＤ法に倣って繰り込まれた粒子系を取得する（Ｓ１２）。磁気モーメント付与部１１０は、取得された粒子系の粒子に磁気モーメントを付与する（Ｓ１４）。磁場演算部１２１は、粒子系によって生成される磁場を演算する（Ｓ１６）。力演算部１２２は、粒子間の距離と各粒子の磁気モーメントから粒子に働く力を演算する（Ｓ１８）。粒子状態演算部１２４は、演算された力を含む粒子の運動方程式から粒子の速度と位置を演算する（Ｓ２０）。状態更新部１２６は、粒子データ保持部１１４に保持される粒子の位置、速度を演算された位置、速度で更新する（Ｓ２２）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２４）。終了条件が満たされない場合（Ｓ２４のＮ）、処理はＳ１６〜Ｓ２２のステップを繰り返す。終了条件が満たされる場合（Ｓ２４のＹ）、表示制御部１１８は、演算結果を出力する（Ｓ２６）。
【００４２】
粒子法は一般に、流体や分離現象や切断加工などの動的な解析対象をより好適に扱える。
そして本実施の形態に係る解析装置１００によると、粒子法による解析と磁気モーメントに基づく磁場解析とを連成することができる。したがって、動的な解析対象を、磁場などの磁気的な性質も含めてより精度高くかつ短時間で解析することが可能となる。
【００４３】
例えば解析対象が砂鉄の場合、従来の有限要素法などのメッシュベースの解析手法ではメッシュを好適に設定することができないので、うまくシミュレーションすることが困難である。そこで本実施の形態に係る解析装置１００によると、砂鉄の磁性を考慮に入れた形で砂鉄を粒子法で好適に解析することができる。
【００４４】
また例えば解析対象がモータの場合、固定体に対して回転体が回転するので、やはりメッシュベースの解析手法は不向きであり静解析に止まる。そこで、本実施の形態に係る解析装置１００によると、回転体の回転を粒子法により高精度にシミュレーションしつつ、回転体と固定体との磁気的相互作用を動的に扱うことができる。さらに、誘導モータのように電磁誘導現象を含む解析対象であっても、電磁誘導現象を考慮に入れた誘導磁化と、誘導磁化に基づいた磁気モーメントを付与することにより、磁場の時間変化を加味した磁気的相互作用を扱うことができる。
【００４５】
以下、解析装置１００で使用される解析手法の原理について説明する。
【００４６】
磁場はマクロなベクトルポテンシャル：
【数１】
0006618417.tif 000002
から以下の様に求まる。スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにある１個のスピンが作る磁場はスピンＳｊｚ（磁気モーメントｍｊｚ）を発生させる円電流を
【数２】
0006618417.tif 000003
とし、
【数３】
0006618417.tif 000004
ここで、
【数４】
0006618417.tif 000005
ｎ＝０のみを残すと一様に磁化：
【数５】
0006618417.tif 000006
した球の作る磁場と同等になる。
なお、真球のバルク材に本手法を適用する場合、ｎ＝０のみ採用すると、一様な磁化が得られず、磁気モーメント法と同様な結果になる。高次の項を加える事で改善できる。ＦＥＭは一様になる。
【００４７】
デカルト座標系に変換（
【数６】
0006618417.tif 000007
を使って、)すれば
【数７】
0006618417.tif 000008
ここで、ｓｉｎθｊ＝０のときは、Ｂθ＝０から、Ｂｘ＝Ｂｙ＝０であることに注意する。磁気モーメントｍｘ、ｍｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場
【数８】
0006618417.tif 000009
も同様に求まる。
【数９】
0006618417.tif 000010
【００４８】
【数１０】
0006618417.tif 000011
における磁場
【数１１】
0006618417.tif 000012
（Ｎ個のスピンからの寄与）は
【数１２】
0006618417.tif 000013
ここで、
【数１３】
0006618417.tif 000014
は全外場を意味する。Ｐ−１４とＰ−１５を連立することで, 磁気モーメントｍｐ が求まる。磁性体外部に生じる磁場は、Ｐ−１５であり、磁性体内部に生じる磁場は、
【数１４】
0006618417.tif 000015
Ｐ−１４は球体表面に生じる磁化による反磁場が考慮され、
【数１５】
0006618417.tif 000016
を満足する。
【００４９】
Ｐ−１３はマトリクスが優対角でない(ｊ＝ｐが除外されている) から収束は遅い。そこで、以下の運動方程式から解を得る。
【数１６】
0006618417.tif 000017
ここで、ｍｖ＝５×１０−２ｄｔ２は仮想質量、γ＝ｍｖ／（１０ｄｔ）は、減衰係数である。
【００５０】
ＦＩＲＥ（非特許文献１参照）を追加すればさらに高速化が可能。以下にＦＩＲＥのコードを示す。なお、粒子数８０２、残差＜１０−８における計算で１０．０５[ｓ］である。ＦＩＲＥ無しでは、残差が１０−５以下にならなかった。
【数１７】
0006618417.tif 000018
【００５１】
以下にLegendre 関数と倍Legendre 関数の特別な場合を列挙する。（ｘ：＝ｃｏｓθｊ）
【数１８】
0006618417.tif 000019
【００５２】
以下に、ｘ軸方向にＢｏｘ＝０．１［Ｔ］の一様な外場においた磁性球（帯磁率が９９９）の磁化
【数１９】
0006618417.tif 000020
を計算した結果を示す。球を構成する粒子（原子）数は４００９個とし、ｆｃｃ構造に配置した。厳密解は
【数２０】
0006618417.tif 000021
である。図４、図５、図６、図７に多重極展開の１次〜４次までの結果を示す。共に中心を通過する断面を表示した。多重極展開の次数を上げると, 厳密解と良く合う事が判る。次数が低いと収束が悪い。また計算時間の大半は球面調和関数の呼び出しである。
【００５３】
図４は、ｎ＝１次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３４６２３１[Ｔ］である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。残差１ｅ−１３以下で、計算時間は約８０分であった。
【００５４】
図５は、ｎ＝２次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３２１８９１[Ｔ]である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。残差１ｅ−８以下で、計算時間は約３ 時間であった。残差は１ｅ−９以下に落ちなかった。
【００５５】
図６は、ｎ＝３次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３１２４４４[Ｔ]である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる.
【００５６】
図７は、ｎ＝４次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３１２２２２[Ｔ]である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。約５時間後において、残差は１ｅ−９以下に落ちなかった。
【００５７】
図８は、球面調和関数の次数と計算値との関係を示す図である。次数を上げると厳密解に近づくことが分かる。
【００５８】
解析装置１００で使用される解析手法の原理を、次のように説明することもできる。
【００５９】
【数２１】
0006618417.tif 000022
Ｍｐは原子核の質量、ｍｅは電子の質量である。ｅは原子核と電子の電荷(の絶対値)である。最終的に電荷は電流に置き換わり、電子の質量は原子核の質量と和して原子の質量Ｍとなるから、陽には現れない。δＶは磁性球の体積、ｍｖは仮想質量、λｉはラグランジェの未定常数である。厳密に拘束を満足させるには、ＳＨＡＫＥ法に代表される収束計算が必要であるが、本手法ではλｉ＝１として、正解の周りの微小振動を誤差として扱う。
【００６０】
スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにあるスピンＳｊｚ（磁気モーメントｍｊｚ）を発生させる円電流は
【数２２】
0006618417.tif 000023
であるから、磁場はマクロなベクトルポテンシャル：
【数２３】
0006618417.tif 000024
から以下の様に求まる。
【数２４】
0006618417.tif 000025
ここで、
【数２５】
0006618417.tif 000026
ｎ＝０のみを残すと一様に磁化：
【数２６】
0006618417.tif 000027
した球の作る磁場と同等になる。
【００６１】
デカルト座標系に変換(
【数２７】
0006618417.tif 000028
を使って、）すれば
【数２８】
0006618417.tif 000029
ここで、ｓｉｎθｊ＝０のときは、Ｈθ＝０から、Ｈｘ＝Ｈｙ＝０であることに注意する。磁気モーメントｍｘ、ｍｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場
【数２９】
0006618417.tif 000030
も同様に求まる。
【数３０】
0006618417.tif 000031
【００６２】
【数３１】
0006618417.tif 000032
におけるモーメントが感じる外磁場
【数３２】
0006618417.tif 000033
(Ｎ個のスピンからの寄与) は
【数３３】
0006618417.tif 000034
誘導磁場
【数３４】
0006618417.tif 000035
は、ｐ点が導体内部である場合、
【数３５】
0006618417.tif 000036
ｐ点が導体外部である場合、
【数３６】
0006618417.tif 000037
Ｃ１、Ｃ２は、例えば、非特許文献２に示される式において
【数３７】
0006618417.tif 000038
となる極限と比較したときの補償係数であり、Ｃ１＝１／３、Ｃ２＝３／５の値をとる。
【００６３】
ここで、Ｑ−２１の右辺第１項（またはＱ−２１−２）は、時間変化する外部磁場により各粒子に誘起される誘導磁化Ｍａｎｔを表す項である。また、Ｑ−２１の右辺第２項およびＱ−２２の右辺は、Ｑ−２１−４式で表される各粒子の誘導磁化に基づいた磁気モーメントｍａｎｔ同士の相互作用により得られる磁場を表す項である。誘導磁場をこのように表すことで、時間変化する外部磁場による影響を考慮した動磁場解析を実現することができる。
【００６４】
なお、繰り込みに際し、電気伝導度σは
【数３８】
0006618417.tif 000039
であり、今はδ＝２としている。
【００６５】
Ｑ−２１−２右辺のLagrange 微分
【数３９】
0006618417.tif 000040
は、粒子法であるから
【数４０】
0006618417.tif 000041
で置き換えることが出来る。ただし、磁化の並進移動は粒子の移動により考慮されるが、磁化ベクトルの回転は別に考慮する必要がある(ランダウ=リフシッツ、電磁気学、§５１)。
【００６６】
磁性体内部に生じる磁場は、
【数４１】
0006618417.tif 000042
であり、全外場ＨＯ（ｘｐ）と、磁化Ｍ（Ｈ）との関係は、
【００６７】
【数４２】
0006618417.tif 000043
である。Ｑ−２６とＱ−２７を連立することで、磁気モーメント
【数４３】
0006618417.tif 000044
が求まる。また、磁性体外部の磁場は、
【数４４】
0006618417.tif 000045
となる。
【００６８】
非線形磁化とヒステリシス
Ｑ−２７とＱ−２４から磁場
【数４５】
0006618417.tif 000046
を消去する：
【数４６】
0006618417.tif 000047
Ｑ−２８を解くことで
【数４７】
0006618417.tif 000048
が得られる（図９にはＢ＝ｆ（Ｈ）の場合を示す）。
【００６９】
図９は、ヒステリシスカーブを示す図である。ヒステリシスカーブとＢ＋２μｏＨ＝３Ｂｏの交点からＭを求める。今
【数４８】
0006618417.tif 000049
を永久磁化とするＢ−Ｈカーブ：
【数４９】
0006618417.tif 000050
が与えられたとする。
【数５０】
0006618417.tif 000051
と連立し
【数５１】
0006618417.tif 000052
を得る。一般のＢ＝ｆ（Ｈ）に対しては、Newton-Laphson法で解く。
【００７０】
磁場の粒子化法
Ｑ−１９はマトリクスが優対角でない（ｊ＝ｐが除外されている）から収束は遅い。そこで、以下の運動方程式から解を得る。
【数５２】
0006618417.tif 000053
ここで、ＨＯ（ｘｐ）は、Ｑ−２６により表される。また、ｍｖは仮想質量である。推奨値は２５０×ｄｔ×ｄｔである。ｍｖを十分小さく取れば、正解の周りの微小な減衰振動が起こり、誤差を小さくできる。
【００７１】
模型を考える：
【数５３】
0006618417.tif 000054
減衰項を持つ強制振動の解(Landau and Lifshitz,Mechanics,Sec.26) は
【数５４】
0006618417.tif 000055
【数５５】
0006618417.tif 000056
を考えれば、
【数５６】
0006618417.tif 000057
と成る。即ちｍｖを十分小さく取れば、解くべき本来の方程式
【数５７】
0006618417.tif 000058
の解に一致する。
【００７２】
離散化
蛙飛び法（速度はｎとｎ＋１の中間点における評価）に従い離散化すれば、収束回数をｓと書き、
【数５８】
0006618417.tif 000059
ここで、誘導磁場Ｈｉｎｄ（ｘｉ）の離散化において、
【数５９】
0006618417.tif 000060
ととる。このとき、対角要素のみ未来を使えば、陽的解法に出来、約３０％程度効率改善になる。計算の大半は特殊関数の呼び出しである。
【００７３】
ＦＩＲＥによる高速化：静磁場解析にのみ有効
ＦＩＲＥ（非特許文献１参照）を追加すればさらに高速化が可能。以下にＦＩＲＥのコードを示す。なお、粒子数８０２、残差＜１０−８における計算で１０．０５[ｓ]である。ＦＩＲＥ無しでは、残差が１０−５以下にならなかった。
【数６０】
0006618417.tif 000061
【００７４】
精度検証：球体の一様磁化
図１０は、ビーズにより構成される導体球の解析モデルを示す図である。精度検証では、ｘ軸方向にＨｅｘｔ，ｘ＝ＨＯｓｉｎ（２πｆｔ）となる空間一様な外部磁界を印加した。ここで、ＨＯ＝７．９５８×１０４［Ａ／ｍ］、ｆ＝２００［Ｈｚ］、導体球の半径Ａ＝１０［ｍｍ］、電気電導率σ＝５９×１０６［Ｓ／ｍ］とした。導体球を構成する粒子（原子）数は、６０９９個とし、ｆｃｃ構造に配置した。
【００７５】
図１１は、導体球表面での時間変化する磁界の計算値を示すグラフであり、導体表面近傍に位置する粒子の重心位置（ｒｘ＝９．８９［ｍｍ］、ｒｙ＝ｒｚ＝０）上の磁界Ｈｘの時間変化を示す。なお、本図では、シミュレーションによる計算結果とともに、外部磁界および厳密解による値も示している。誘導磁界の影響により、導体内部の磁界が外部磁界に対して位相が遅れて変化する様子を再現することができており、計算結果と厳密解の傾向は良く一致することが示された。
【００７６】
図１２は、計算値の位相誤差と振幅誤差を示すグラフであり、ｘ軸上に位置する粒子の重心位置における磁界の位相誤差εｐｈａｓｅと振幅誤差εａｍｐｌｉｔｕｄｅを示す。位相誤差εｐｈａｓｅおよび振幅誤差εａｍｐｌｉｔｕｄｅは、
【数６１】
0006618417.tif 000062
と定義した。ここで、θｅｘａｃｔは厳密解の位相、θＭＢＭは解析結果の位相である。また、Ｈｅｘａcｔは厳密解により得られた磁界、ＨＭＢＭは解析結果の磁界である。なお位相誤差は、Ｈｅｘａcｔが最大となる位相θｅｘａｃｔと、ＨＭＢＭが最大となる位相θＭＢＭとを比較した。図１２より、ｘ軸上における解析結果と厳密解の磁界の位相誤差は最大で５．８６％、振幅誤差は６．３８％となり、厳密解と比較して精度の高い解析結果が得られることが示された。
【００７７】
Legendre関数の級数表示
以下にLegendre関数と倍Legendre関数の特別な場合を列挙する。（ｘ：＝ｃｏｓθｊ）
【数６２】
0006618417.tif 000063
【００７８】
粒子に働く力の導出
粒子（電子と核）と、粒子磁場間の相互作用を記述するLagrangeanは
【数６３】
0006618417.tif 000064
である。ベクトルポテンシャル
【数６４】
0006618417.tif 000065
と電磁場の関係は、
【数６５】
0006618417.tif 000066
である。
【数６６】
0006618417.tif 000067
を使い、以下の運動方程式を得る。
【数６７】
0006618417.tif 000068
原子Ｍ＝Ｍｐ＋ｍｅに対する運動方程式は、
【数６８】
0006618417.tif 000069
となる。右辺：Lorentz力を平均電流密度
【数６９】
0006618417.tif 000070
、電気モーメントを
【数７０】
0006618417.tif 000071
で置き換えると
【数７１】
0006618417.tif 000072
Maxwellの方程式とＱ−２４から
【数７２】
0006618417.tif 000073
最終的に原子に働く力は、電気分極を無視し、
【数７３】
0006618417.tif 000074
となる。
【００７９】
粒子に働く力は、直接Maxwellの応力テンソルからも求まる。
【数７４】
0006618417.tif 000075
球に局在した電流が作るモーメント
【数７５】
0006618417.tif 000076
に作用する力は、最低の次数では
【数７６】
0006618417.tif 000077
である。従い、ポテンシャルエネルギーは、電気分極によるポテンシャルエネルギーを合わせて、
【数７７】
0006618417.tif 000078
となるが、渦電流や真電流を包含しない。
【００８０】
【数７８】
0006618417.tif 000079
結局、
【数７９】
0006618417.tif 000080
を計算すれば磁性球に働く磁力が求まる。
【００８１】
磁化による電流
【数８０】
0006618417.tif 000081
と渦電流
【数８１】
0006618417.tif 000082
を完全に切り分けることは恐らく出来ないが、それぞれ、
【数８２】
0006618417.tif 000083
【００８２】
微係数の計算
以下に必要な微分係数を書き出す。
【数８３】
0006618417.tif 000084
等からｒｐｊ＞ａのみ列挙すると（なお、
【数８４】
0006618417.tif 000085
はｓｉｎθに比例するから、θ＝０で発散は無い）、
【数８５】
0006618417.tif 000086
【数８６】
0006618417.tif 000087
【００８３】
磁気モーメントが作る磁化電流の計算
ｚ軸に平行に磁化したとき
今一度、ｚ軸に磁化したときの磁場を書く
【数８７】
0006618417.tif 000088
磁化電流を
【数８８】
0006618417.tif 000089
と書けば、磁化により生じる磁力
【数８９】
0006618417.tif 000090
は
【数９０】
0006618417.tif 000091
である。先ず、
【数９１】
0006618417.tif 000092
を求める。
【数９２】
0006618417.tif 000093
【００８４】
ｘ軸に平行に磁化したとき
ｘ軸に磁化したときの磁場を書く：
【数９３】
0006618417.tif 000094
【数９４】
0006618417.tif 000095
磁化電流
【数９５】
0006618417.tif 000096
を求める。
【数９６】
0006618417.tif 000097
【００８５】
ｙ軸に平行に磁化したとき
ｙ軸に磁化したときの磁場を書く：
【数９７】
0006618417.tif 000098
磁化電流
【数９８】
0006618417.tif 000099
を求める。
【数９９】
0006618417.tif 000100
【００８６】
誘導磁界が作る渦電流の計算
導体内部の誘導磁場は、Ｑ−２１より、
【数１００】
0006618417.tif 000101
となる。渦電流は
【数１０１】
0006618417.tif 000102
なお、上式Ｒ−５２−２は、Ｒ−１２において、球面調和関数の展開次数ｎ＝０とした場合と同じである。
【００８７】
粒子に働く力の計算
ｍｘ、ｍｙ、ｍｚと渦電流による寄与を加えれば、全電流が得られる。
【数１０２】
0006618417.tif 000103
座標
【数１０３】
0006618417.tif 000104
に粒子に働く力は
【数１０４】
0006618417.tif 000105
最後に
【数１０５】
0006618417.tif 000106
の微分から出てくる項を付加する必要がある。これは磁場の誤差
【数１０６】
0006618417.tif 000107
から生じる見かけの力であり、運動と磁場の全エネルギーを保存させる。
【数１０７】
0006618417.tif 000108
【００８８】
（第２の実施の形態）
上述の第１の実施の形態では、Ｑ−２１の右辺第１項に示される誘導磁化と、右辺第２項に示される磁気モーメントの相互作用とにより表される誘導磁場の式（Ｑ−２１）を用いて動磁場解析を行う場合を示した。本実施の形態では、Ｑ−２１の式を用いる代わりに、ベクトルポテンシャルの微分方程式から導出されるベクトルポテンシャルの厳密解を利用して誘導磁場を定式化する。以下、第１の実施の形態との相違点を中心に述べる。
【００８９】
ベクトルポテンシャルＡについて、以下の微分方程式が成立する。
【数１０８】
0006618417.tif 000109
上記微分方程式の解は、
【数１０９】
0006618417.tif 000110
となる。ここで、
【数１１０】
0006618417.tif 000111
として、各粒子に分割すれば、各粒子の位置ｒｉでのベクトルポテンシャルは、
【数１１１】
0006618417.tif 000112
となる。ここで、ｖｊは各粒子の体積である。なお、ｊの和に際し、ベクトル
【数１１２】
0006618417.tif 000113
が粒子表面を横切る場合には和を取らない。また、粒子径ａが満たすべき条件は、
【数１１３】
0006618417.tif 000114
である。
【００９０】
このとき、各対象物を粒子に分割することにより、ベクトルポテンシャルＡの並進対称性を失うおそれがある。そこで、本発明者は、下記Ｓ−４に示されるゲージ変換を導入することにより、ベクトルポテンシャルＡの並進対称性を回復できることを見いだした。なお、このゲージ変換を本明細書において「ビーズゲージ」ともいう。
【数１１４】
0006618417.tif 000115
ここで、ベクトルｇｉは、計算対象とする粒子群（粒子数：Ｎｉ）に関する局所重心ベクトルであり、
【数１１５】
0006618417.tif 000116
で表される。したがって、ベクトルｇｉは、粒子群の重心位置を表している。また、Ｓ−４に示すビーズゲージは、磁束密度Ｂｊと局所重心ベクトルｇｉの外積を用いて記述されると言える。
【００９１】
上記のビーズゲージについて、磁束密度ＢとベクトルポテンシャルＡの関係式
【数１１６】
0006618417.tif 000117
を求めれば、
【数１１７】
0006618417.tif 000118
となる。したがって、粒子数Ｎｉが十分に大きければ、並進対称性が満たされることがわかる。
【００９２】
ここで、Ｓ−４に示すゲージ変換から、Ｓ−３の右辺に含まれるベクトルポテンシャルの時間微分
【数１１８】
0006618417.tif 000119
を磁束密度Ｂの時間微分を用いて以下のように表すことができる。
【数１１９】
0006618417.tif 000120
さらに、ベクトルｄｊｉを、
【数１２０】
0006618417.tif 000121
とすると、ベクトルポテンシャルの時間微分は、
【数１２１】
0006618417.tif 000122
と表される。したがって、Ｓ−３にＳ−９を代入して積分を実行すれば、ベクトルポテンシャルＡの厳密解を磁束密度Ｂの時間微分を用いて記述できる。
【００９３】
ベクトルポテンシャルＡの厳密解を導出するため、まず、磁束密度Ｂのｚ成分の時間微分である、
【数１２２】
0006618417.tif 000123
のみが作用する場合の積分を実行する。これにより、ベクトルポテンシャルＡは、
【数１２３】
0006618417.tif 000124
となる。磁束密度Ｂのｘ成分およびｙ成分の時間微分である、
【数１２４】
0006618417.tif 000125
の寄与も加えれば、
【数１２５】
0006618417.tif 000126
が得られる。
【００９４】
解析対象とする導体内部を記述する誘導磁場Ｈは、
【数１２６】
0006618417.tif 000127
の関係式より、
【数１２７】
0006618417.tif 000128
となる。ここで、
【数１２８】
0006618417.tif 000129
と置換すると、誘導磁場Ｈの式として、
【数１２９】
0006618417.tif 000130
が得られる。
【００９５】
磁場演算部１２１は、粒子系により生成される誘導磁場を、Ｓ−１３に示す誘導磁場の式を用いて演算する。誘導磁場の時間変化値を得るための計算は、上述の第１の実施の形態で示したものと同様の方法を用いればよい。磁場演算部１２１は、得られた誘導磁場から電流密度を計算してもよい。また、力演算部１２２は、得られた磁場および電流の値から、各粒子に作用する電磁力を計算してもよい。
【００９６】
本実施の形態において、Ｓ−５に示す局所重心ベクトルｇｉは、互いに同電位となる粒子群に対して導入され、電気的に互いに絶縁される複数の粒子群に対しては、それぞれ別の局所重心ベクトルが導入される。例えば、解析対象とする粒子系において、互いが電気的に絶縁される複数の部分を有する場合、それぞれの部分を構成する粒子群に対して重心位置の異なる局所重心ベクトルが適用される。また、解析対象とする物体が時間経過とともに変形して複数の物体に分裂し、分裂した部分同士が互いに絶縁した状態となる場合には、それぞれの部分を構成する粒子群に対して重心位置の異なる局所重心ベクトルを適用する。
【００９７】
したがって、解析対象の粒子系に、電気的に互いが絶縁される第１粒子群および第２粒子群を含む複数の粒子群が配置される場合、磁場演算部１２１は、それぞれの粒子群に対応する誘導磁場の式を用いて解析をする。例えば、第１粒子群における誘導磁場は、第１粒子群の重心位置を表す第１局所重心ベクトルを用いたゲージ変換により導出される第１の誘導磁場の式を用いて演算される。また、第２粒子群における誘導磁場は、第２粒子群の重心位置を表す第２局所重心ベクトルを用いたゲージ変換により導出される第２の誘導磁場の式を用いて演算される。より具体的には、Ｓ−１３に示す誘導磁場の式において、右辺第１項に含まれるｄｉｊがそれぞれの粒子群での計算において異なる値をとることとなる。
【００９８】
本実施の形態に係る誘導磁場の式（Ｓ−１３）を用いた解析手法について精度検証を行った。第１の実施の形態と同様に、図１０に示す導体球の解析モデルを用い、ｘ軸方向にＨｅｘｔ，ｘ＝ＨＯｓｉｎ（２πｆｔ）となる空間一様な外部磁界を印加した。計算条件は、ＨＯ＝７．９４８×１０４［Ａ／ｍ］、ｆ＝１００［Ｈｚ］、導体球の半径Ａ＝１０［ｍｍ］、電気電導率σ＝５９×１０６［Ｓ／ｍ］とした。導体球を構成する粒子（原子）数は、６０９９個とし、ｆｃｃ構造に配置した。
【００９９】
図１３は、導体球表面での磁界、電流密度および電磁力の計算値を示すグラフである。本図は、導体表面近傍に位置する粒子の重心位置における、（ａ）磁界Ｈｘ、（ｂ）電流密度ｊｚ、（ｃ）電磁力Ｆｙの時間変化を示す。本図では、厳密解の値を実線で、計算値を破線で示している。図示されるように、本実施の形態に係る解析手法を用いると、計算結果と厳密解がよく一致することがわかった。
【０１００】
（第３の実施の形態）
上述の第１の実施の形態では、対象物を粒子に分割し、各粒子の球対称性に基づく厳密解を利用することで、高速かつ精度の優れた演算結果を得ている。しかしながら、連続体である対象物を粒子で分割すると粒子間に隙間が生じるため、その隙間に起因する計算誤差がどうしても生じてしまう。特に、粒子間の距離が近い場合の粒子近傍の磁場を計算する場合に誤差の影響が大きくなる。計算精度を高めるためには、対象物を立方体やボロノイ多面体といったボクセルに分割し、これらのボクセルにより隙間なく対象物を分割する必要がある。この場合、ボクセルの形状に応じた厳密解に基づく面積分を実行する必要があるため計算時間が膨大となる。
【０１０１】
本実施の形態では、粒子近傍の磁場を演算する場合にボクセル形状に応じた厳密解を利用し、近似による誤差がほぼ無視できる粒子遠方の磁場を演算する場合に双極子近似を利用する。このとき、ボクセル形状に応じた厳密解の計算を都度実行するのではなく、あらかじめボクセルに付与される磁気モーメントのベクトル量とそのボクセルが作る磁場のベクトル量との対応関係を厳密解に基づいて算出しておき、その対応関係を参照することで各粒子の近傍における磁場を演算する。したがって、本実施の形態によれば、負荷の高い面積分の都度演算を回避しつつ、あらかじめ算出した厳密解に基づく対応関係を用いて磁場を高精度かつ高速に演算することができる。以下、本実施の形態について、上述の実施の形態との相違点を中心に説明する。
【０１０２】
仮想空間内に定義される観測点ｒｉの磁場は、各微小体積ｄｖ’の磁化Ｍ（ｒ’）が作る磁場を全空間にわたって積分することにより得られ、下記式（Ｔ−１）により表される。式（Ｔ−１）は、体積分の項と面積分の項に分けることによって、式（Ｔ−２）に変形できる。なお、Ｖは領域体積、Ｓは表面積、ｎＳは表面の法線ベクトルである。
【数１３０】
0006618417.tif 000131
【０１０３】
磁気ビーズ法では、対象空間をＮ個の粒子（ビーズ）に分割し、各ビーズ内部の磁化ベクトルＭｊが一定であるとし、式（Ｔ−１）を双極子近似して得られる下記式（Ｔ−３）により座標ｒｉの磁場ベクトルＨ（ｒｉ）を計算する。ここで、ｒｉｊ＝ｒｉ−ｒｊ，ｎｉｊ＝ｒｉｊ／｜ｒｉｊ｜であり、ΔＶｊはビーズ体積であり、γｊはビーズの充填率である。
【数１３１】
0006618417.tif 000132
【０１０４】
一方、磁気モーメント法では、対象をＮ個の要素（ボクセル）に分割し、要素内部の磁化ベクトルＭｊは一定であるとし、式（Ｔ−２）を変形して得られる下記式（Ｔ−４）により座標ｒｉの磁場ベクトルＨ（ｒｉ）を計算する。ここで、ｄｓｊは各要素の微小表面を意味し、ボクセルが正六面体（立方体）であれば、６つの表面について面積分を実行する必要がある。
【数１３２】
0006618417.tif 000133
【０１０５】
磁気ビーズ法と磁気モーメント法を比較すると、計算精度については磁気モーメント法が優れており、計算速度については磁気ビーズ法が優れている。磁気ビーズ法は、連続体をＮ個の球形状の粒子に分割し、各粒子が作る磁場を双極子近似して求めているため、各粒子近傍の磁場の計算精度が低い。その一方で、各粒子から離れた遠方領域の磁場については双極子近似による誤差が無視できるために計算精度が高い。磁気モーメント法は、各要素間に隙間が生じないように連続体を多面体要素に分割するため、磁気ビーズ法において生じるような近似誤差がなく、計算精度が高い。
【０１０６】
磁気モーメント法では、上記式（Ｔ−４）に示される面積分を高精度で実行する必要があり、数値積分法や厳密解での計算が必須である。このとき、多面体要素が有する面数ｎ（立方体であればｎ＝６）の回数だけ面積分を実行する必要があり、積分回数が多くなることによって計算時間が膨大となる。一方、磁気ビーズ法では、ビーズの球形状に基づく厳密解を利用すれば、上記式（Ｔ−３）に示されるように面積分が含まれないため、計算時間が少なくて済む。
【０１０７】
そこで、本実施の形態では、計算精度と計算速度を両立させるために、双極子近似が成立する遠方領域内の磁場の計算では、磁気ビーズ法に基づく式（Ｔ−３）を利用し、双極子近似が成立しない近傍領域内の磁場の計算では、磁気モーメント法に基づく式（Ｔ−４）を利用する。さらに、式（Ｔ−４）に基づく演算処理を高速化させるために、面積分の演算結果をあらかじめ算出し、磁化ベクトルＭｊと、磁化ベクトルＭｊがその周りに作る磁場Ｈｊとの対応関係を定める磁場テーブルを用意する。シミュレーション実行時の演算処理では、面積分を都度実行する代わりにあらかじめ用意した磁場テーブルを参照して磁場を演算する。
【０１０８】
図１４は、基準粒子Ｖｊの周りに設定される近傍領域ＣＡおよび遠方領域ＣＢを模式的に示す図である。図１４は、対象物が立方体形状の要素（ボクセル）で分割される場合を示しており、各要素の頂点に格子点Ｇｋが設定される。各立方体要素の一辺の長さは、粒子の半径ａの２倍（２ａ）である。基準粒子（基準要素）Ｖｊを中心として直交座標系（ξ，η，ζ）が設定され、立方体の中心に各要素の座標ｒｊが設定される。基準粒子Ｖｊの周りには半径ａｓの仮想的な球が定められ、この球の内部が近傍領域ＣＡとなり、外側が遠方領域ＣＢとなる。近傍領域ＣＡを規定する球の半径ａｓは、双極子近似が成立する距離に応じて決定され、例えば、粒子径ａの４倍以上、好ましくは粒子径ａの６倍以上の値が半径ａｓとして設定される。図示する例では、球の半径ａｓを粒子径ａの６倍（つまり、ａｓ＝６ａ）としている。
【０１０９】
図１４において、基準粒子Ｖｊが周りに作る磁場のうち、遠方領域ＣＢの任意の座標ｒｉにおける磁場ＨＢ（ｒｉ）は、式（Ｔ−３）に基づいて算出される。一方、基準粒子Ｖｊが周りに作る磁場のうち、近傍領域ＣＡの内部の任意の座標ｒｉにおける磁場ＨＡ（ｒｉ）は、座標ｒｉに近接する複数の格子点Ｇ１〜Ｇ４（実際には、後述する図１５に示す８個の格子点）の磁場を平均化することにより算出される。各格子点の磁場は、式（Ｔ−４）に基づいて算出されるものであり、格子点ごとにあらかじめ算出される磁場テーブルの値を用いて演算される。
【０１１０】
図１５は、近接する複数の格子点の磁場Ｈ１〜Ｈ８から観測点ｒｉの磁場ＨＡを算出する方法を模式的に示す図である。近傍領域ＣＡ内の磁場ＨＡは、磁場ＨＡの演算対象となる座標（ξ，η，ζ）を囲む８個の格子点における磁場Ｈｋ（ｋ＝１〜８）を平均化することにより算出され、下記式（Ｔ−５）で表される。なお、Ｎｋ（ξ，η，ζ）は、式（Ｔ−６）で表される。
【数１３３】
0006618417.tif 000134
【数１３４】
0006618417.tif 000135
【０１１１】
ここで、基準粒子Ｖｊが各格子点Ｇｋ（位置ｒｉ）に作る磁場Ｈｋは、式（Ｔ−４）から下記式（Ｔ−７）で表すことができる。
【数１３５】
0006618417.tif 000136
【０１１２】
式（Ｔ−７）で表される磁場Ｈｋは、面積分を利用した厳密解である。ここで、式（Ｔ−７）を変形し、磁場Ｈｋと磁化ベクトルＭｊの対応関係を定めるテンソルＴｋを用いて下記式（Ｔ−８）のように表すことができれば、任意の磁化ベクトルＭｊを有する粒子が各格子点Ｇｋに作る磁場Ｈｋを容易に算出できるようになる。式（Ｔ−７）から、テンソルＴｋの各成分は、下記式（Ｔ−９）で表すことができる。なお、ｎｘ，ｎｙ，ｎｚはそれぞれｘ方向、ｙ方向、ｚ方向の単位ベクトルであり、ｘｉｊ，ｙｉｊ，ｚｉｊはそれぞれ、基準粒子Ｖｊから観測点（格子点）までの距離ｒｉｊのｘ成分、ｙ成分、ｚ成分である。
【数１３６】
0006618417.tif 000137
【数１３７】
0006618417.tif 000138
【０１１３】
そこで、本実施の形態では、式（Ｔ−９）で表されるテンソルＴｋの各成分の数値解をあらかじめ算出し、その数値解と磁化ベクトルＭｊの積から各格子点の磁場Ｈｋを算出し、さらに、近傍領域ＣＡ内の磁場ＨＡを算出できるようにする。つまり、下記式（Ｔ−１０）を用いて近傍領域ＣＡ内の磁場ＨＡを演算する。あらかじめテンソルＴｋの各成分の数値を求めておけば、式（Ｔ−１０）の計算に面積分が含まれないため、面積分を都度実行する場合よりも非常に短い時間で精度の高い計算結果を得ることができる。
【数１３８】
0006618417.tif 000139
【０１１４】
図１６は、近傍領域ＣＡの任意の座標ｒｉにて磁場ＨＡを演算するために必要となる格子点Ｇｋを模式的に示す図である。仮想的な半径ａｓの球の内部の全ての座標ｒｉで磁場ＨＡを計算できるようにするためには、近傍領域ＣＡを囲うことのできるボクセル群の範囲内（太線の範囲内）の全ての格子点Ｇｋについて、磁場ＨＡと磁化ベクトルＭｊの対応関係を定めるテンソルＴｋを算出しておけばよい。例えば、近傍領域ＣＡを定める範囲がビーズ径ａの６倍であれば、事前計算が必要となる格子点の数は５００個未満である。この数は、仮想空間内に定義される粒子系に含まれる全粒子数よりも少ない。そのため、全粒子が作る磁場を面積分を用いて１回計算するよりも短い時間で、磁気的な対応関係を定めるテンソルＴｋを導出することができる。つまり、磁気テーブルを算出するための演算負荷はそれほど高くなく、全計算時間に占める磁気テーブルの算出時間の割合は小さい。
【０１１５】
つづいて、磁気的な対応関係を定めるテンソルＴｋの具体的な計算方法の一例を示す。テンソルＴｋの各成分の数値解の演算方法は、式（Ｔ−９）の計算が実行できれば特に問わないが、例えば、非特許文献３に開示される方法を用いることができる。図１７（ａ）、（ｂ）、（ｃ）は、立方体表面の面積分を実行する方法の一例を模式的に示す図である。まず、面積分の対象となる六つの面の一つを取り出し、取り出した正方形を二つに分割した微小な三角形を考える。このようにして六つの面を１２個の微小な三角形に分割し、微小三角形ごとに面積分を実行する。
【０１１６】
図１８は、微小三角形の面積分を実行する方法を模式的に示す図である。図１８に示すローカル座標系（ｘ’，ｙ’，ｚ’）において、面積分の対象となる微小三角形△１２３がｘ’ｙ’平面内に配置され、微小三角形△１２３により作られる磁場ΔＨｋの座標ｒｉがｚ’軸上に配置されるように座標変換を施す。このとき、微小三角形△１２３により作られる磁場ΔＨｋについて、ΔＨｋ＝ΔＴｋＭｊと表すと、微小三角形Δ１２３に起因するテンソルΔＴｋは、下記式（Ｔ−１１）で表される。式（Ｔ−１１）に含まれる各変数は、式（Ｔ−１２）で定義される。
【０１１７】
【数１３９】
0006618417.tif 000140
【数１４０】
0006618417.tif 000141
【０１１８】
このようにして微小三角形Δ１２３に起因して得られるテンソルΔＴｋを基準粒子Ｖｊの全ての面について演算して足し合わせることで、磁気的な対応関係を定めるテンソルＴｋを導出できる。磁気的な対応関係を定めるテンソルＴｋは、近傍領域ＣＡの任意の座標ｒｉで磁場ＨＡを演算するために必要となる全て格子点Ｇｋについて計算され、その数値解が磁気テーブルとして保持される。
【０１１９】
図１９は、実施の形態に係る解析装置２００の機能および構成を示すブロック図である。本実施の形態に係る解析装置２００は、基準データ保持部２３０をさらに備える点で図１に示す解析装置１００と相違する。基準データ保持部２３０は、基準粒子に付与される磁気モーメントのベクトル量と、基準粒子が周りに作る磁場の厳密解のベクトル量との対応関係を保持する。つまり、基準データ保持部２３０は、磁気的な対応関係を定めるテンソルＴｋの数値解が格子点Ｇｋごとに定められた磁気テーブルを保持する。
【０１２０】
磁場演算部１２１は、粒子系取得部１０８および磁気モーメント付与部１１０によって定義された粒子系に応じて、上述の磁気テーブルを導出するための演算処理を実行する。磁場演算部１２１は、粒子系の定義に応じて近傍領域ＣＡを設定し、磁気的な対応関係を定めるテンソルＴｋの算出が必要となる複数の格子点Ｇｋを特定する。磁場演算部１２１は、原点に配置した基準粒子Ｖｊが各格子点Ｇｋの位置に作る磁場Ｈｋに関連する演算を実行し、基準粒子Ｖｊと各格子点Ｇｋの配置関係を示す値とともにテンソルＴｋの数値解を基準データ保持部２３０に記憶させる。
【０１２１】
磁場演算部１２１は、粒子ｊが座標ｒｉに作る磁場Ｈｊ（ｒｉ）を演算する場合、座標ｒｉが粒子ｊの近傍領域内であれば磁気テーブルを参照して磁場を演算し、座標ｒｉが粒子ｊの遠方領域（近傍領域外）であれば双極子近似による式を用いて磁場を演算する。磁場演算部１２１は、粒子ｊと座標ｒｉの距離が所定値を超える場合、つまり、座標ｒｉから大きく離れた位置にある粒子ｊについては、磁場の演算対象とせずに無視をしてもよい。例えば、遠方領域のうち、近傍領域に近い中間領域については、双極子近似による磁場演算を実行し、中間領域よりも遠い領域については磁場演算を実行しなくてもよい。中間領域の範囲は、例えば、粒子径ａの１０倍程度である。
【０１２２】
磁場演算部１２１は、磁気テーブルを参照して磁場を演算する場合、座標ｒｉの位置を図１６に示す直交座標系（ξ，η，ζ）に変換して磁場を算出し、その後、元の座標系に戻すことによって座標ｒｉの磁場を演算する。このとき、粒子系の座標系と磁気テーブルの直交座標系（ξ，η，ζ）との間の座標変換は平行移動のみで済むため、座標変換に関する処理負荷を小さくできる。変形例においては、座標変換時に拡大縮小や回転などの座標変換処理を施してもよい。
【０１２３】
図２０は、解析装置２００における一連の処理の一例を示すフローチャートである。粒子系取得部１０８は、ＲＭＤ法に倣って繰り込まれた粒子系を取得する（Ｓ３０）。磁気モーメント付与部１１０は、取得された粒子系の粒子に磁気モーメントを付与する（Ｓ３２）。磁場演算部１２１は、定義された粒子系に応じて磁気テーブルを作成し、基準データ保持部２３０に保持させる（Ｓ３４）。磁場演算部１２１は、粒子系によって生成される磁場を演算する（Ｓ３６）。力演算部１２２は、粒子間の距離と各粒子の磁気モーメントから粒子に働く力を演算する（Ｓ３８）。粒子状態演算部１２４は、演算された力を含む粒子の運動方程式から粒子の速度と位置を演算する（Ｓ４０）。状態更新部１２６は、粒子データ保持部１１４に保持される粒子の位置、速度を演算された位置、速度で更新する（Ｓ４２）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ４４）。終了条件が満たされない場合（Ｓ４４のＮ）、Ｓ３６〜Ｓ４２のステップを繰り返す。終了条件が満たされる場合（Ｓ４４のＹ）、表示制御部１１８は、演算結果を出力する（Ｓ４６）。
【０１２４】
図２１は、図２０のＳ３６の処理の詳細を示すフローチャートである。磁場演算部１２１は、変数ｉに初期値ｉ＝１を設定する（Ｓ５０）。磁場演算部１２１は、基準データ保持部２３０に保持される磁気テーブルを参照して近傍粒子が座標ｒｉに作る磁場を算出し（Ｓ５２）、双極子近似式を用いて遠方粒子が座標ｒｉに作る磁場を算出する（Ｓ５４）。磁場演算部１２１は、算出した磁場を足し合わせ、座標ｒｉの粒子の磁化ベクトルＭｉ（磁気モーメント）をＭｉ＝χＨｉの関係式を用いて演算する（Ｓ５６）。変数ｉに１を加え（Ｓ５８）、変数ｉが粒子数Ｎ以下であれば（Ｓ６０のＮ）、Ｓ５２〜Ｓ５８のステップを繰り返し、変数ｉが粒子数Ｎより大きければ（Ｓ６０のＹ）、Ｓ３６の処理を終了する。
【０１２５】
図２２は、面積分による磁場の厳密解Ｂ１と双極子近似による磁場の数値解Ｂ２を比較したグラフである。図２２において、横軸はソースとなる粒子ｊと観測点ｉの距離ｒｉｊを粒子径ａで規格化した値であり、縦軸はｘ軸方向の磁束密度Ｂｘの値である。図示されるように、粒子ｊと観測点ｉの距離ｒｉｊが２倍となる場合に厳密解Ｂ１（実線）と近似解Ｂ２（破線）の差が大きくなり、距離ｒｉｊが大きくなるにつれて差が小さくなり、距離ｒｉｊが粒子径ａの６倍以上になると差がほとんど生じないことがわかる。このことから、近傍領域ＣＡとして粒子径ａの６倍程度の範囲を設定すればよく、６倍を超える範囲については双極子近似による誤差がほとんど無視できると言える。
【０１２６】
図２３（ａ）、（ｂ）は、本実施の形態に係る計算結果を示す図である。本図は、ｘ方向に１０ｍｍ、ｙ方向およびｚ方向に５ｍｍが設定された直方体形状の磁性体の磁化現象の計算例を示す。外部磁場を一定の０．１［Ｔ］とし、磁化特性は線形（χ＝９９９）とした。図２３（ａ）は、比較例であり、全ての粒子に対して上記式（Ｔ−４）に基づく面積分を実行して磁場を演算した結果を示す。図２３（ｂ）は、本実施の形態に係る磁気テーブルを用いて磁場を演算した結果を示す。両者を比較したところ、計算結果がよく一致することが分かった。また、図２３（ａ）の例では１ステップあたりの計算時間が２．１５秒であるのに対し、図２３（ｂ）の例では１ステップあたりの計算時間が０．０３秒であった。したがって、本実施の形態に係る手法により、高い計算精度を維持したまま約７０倍の計算速度の高速化を実現できた。
【０１２７】
図２４は、本実施の形態に係る計算結果を示し、上述の図２３に示す計算例のｙ軸上の磁束密度分布の絶対値｜Ｂ｜を示す。実線で示す磁束密度分布Ｂ３は、図２３（ｂ）の実施の形態に係る計算例に対応し、破線で磁束密度分布Ｂ４は、図２３（ａ）の比較例に係る計算例に対応する。図示されるように、両者の計算結果がよく一致することが分かる。
【０１２８】
本実施の形態によれば、計算過程においてボクセル形状が変化しないという磁気ビーズ法の特徴をうまく利用して高い計算精度を維持しつつ計算速度の高速化を実現できる。有限要素法では、節点位置の変化によりボクセル形状が変化し、ボクセル表面の面要素の形状が変化するため、変化した面要素の形状に応じた面積分を実行する必要がある。一方、磁気ビーズ法では、ボクセルの形状が計算過程にわたって変わらないため、共通する基準粒子を用いた面積分の結果を流用できる。このように本実施の形態では、磁気ビーズ法の特徴を利用し、共通する基準粒子に基づく磁気テーブルをあらかじめ用意することで計算速度を向上させることができる。
【０１２９】
なお変形例においては、粒子系にボクセル形状の異なる複数種類の粒子を定義してもよい。この場合、ボクセル形状の異なる複数種類の基準粒子ごとにあらかじめ磁気テーブルを用意しておき、任意の観測点ｒｉの位置に磁場を作る粒子ｊのボクセル形状に応じて、異なる種類の磁気テーブルを参照して磁場を算出すればよい。これにより、複数種類のボクセル形状を利用する場合であっても上述の実施の形態と同様の効果を奏することができる。ボクセル形状は、多面体であればどのような形状であってもよく、正四面体、正八面体、正十二面体などであってもよい。また、対象空間をボロノイ分割することにより得られるボロノイ多面体であってもよい。
【０１３０】
上述の実施の形態では、近傍粒子と遠方粒子とで使用する数式を切り替える場合を示した。変形例においては、遠方に存在する粒子（例えば、粒子径ａの１０倍より遠い粒子）については磁場演算をせず、磁場演算となる粒子について全て磁気テーブルを参照して演算してもよい。つまり、粒子径ａが１０倍となる範囲について磁気的な対応関係Ｔｋを算出し、対応関係Ｔｋに基づいて磁場演算を実行してもよい。
【０１３１】
上述の実施の形態では、磁気テーブルにテンソルＴｋが保持される場合を示したが、式（Ｔ−１０）に示される計算の実行に必要な形式であれば、その他の種類の数値解が磁気テーブルとして保持されてもよい。例えば、変形例においては、ＨＡ（ξ，η，ζ）＝Ｓｋ（ξ，η，ζ）Ｍｊと記述されるテンソルＳｋの数値解を磁気テーブルとして保持してもよい。
【０１３２】
以上、実施の形態に係る解析装置の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。
【０１３３】
実施の形態では、数値演算部１２０において粒子の位置と速度の両方を演算する場合について説明したが、これに限られない。例えば、数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法もあり、本実施の形態に係る技術的思想をそのような手法に適用してもよい。
【符号の説明】
【０１３４】
２００  解析装置、  １０２  入力装置、  １０４  ディスプレイ、  １０８  粒子系取得部、  １１０  磁気モーメント付与部、  １１４  粒子データ保持部、  １１８  表示制御部、  １２０  数値演算部、  １２１  磁場演算部、  １２２  力演算部、  １２４  粒子状態演算部、  １２６  状態更新部、  １２８  終了条件判定部、２３０  基準データ保持部。
【図１】
0006618417.tif 000142
【図２】
0006618417.tif 000143
【図３】
0006618417.tif 000144
【図４】
0006618417.tif 000145
【図５】
0006618417.tif 000146
【図６】
0006618417.tif 000147
【図７】
0006618417.tif 000148
【図８】
0006618417.tif 000149
【図９】
0006618417.tif 000150
【図１０】
0006618417.tif 000151
【図１１】
0006618417.tif 000152
【図１２】
0006618417.tif 000153
【図１３】
0006618417.tif 000154
【図１４】
0006618417.tif 000155
【図１５】
0006618417.tif 000156
【図１６】
0006618417.tif 000157
【図１７】
0006618417.tif 000158
【図１８】
0006618417.tif 000159
【図１９】
0006618417.tif 000160
【図２０】
0006618417.tif 000161
【図２１】
0006618417.tif 000162
【図２２】
0006618417.tif 000163
【図２３】
0006618417.tif 000164
【図２４】
0006618417.tif 000165