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