特許 6249912 解析装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6249912

(24)【登録日】2017年12月1日

(45)【発行日】2017年12月20日

(54)【発明の名称】解析装置

(51)【国際特許分類】

   G06F 19/00 20110101AFI20171211BHJP

【ＦＩ】

   G06F19/00 110

【請求項の数】9

【全頁数】40

(21)【出願番号】特願2014-183427(P2014-183427)

(22)【出願日】2014年9月9日

(65)【公開番号】特開2015-111401(P2015-111401A)

(43)【公開日】2015年6月18日

【審査請求日】2016年12月15日

(31)【優先権主張番号】特願2013-228374(P2013-228374)

(32)【優先日】2013年11月1日

(33)【優先権主張国】JP

(73)【特許権者】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105924

【弁理士】

【氏名又は名称】森下 賢樹

(74)【代理人】

【識別番号】100109047

【弁理士】

【氏名又は名称】村田 雄祐

(74)【代理人】

【識別番号】100109081

【弁理士】

【氏名又は名称】三木 友由

(74)【代理人】

【識別番号】100116274

【弁理士】

【氏名又は名称】富所 輝観夫

(72)【発明者】

【氏名】宮崎 修司

(72)【発明者】

【氏名】市嶋 大路

【審査官】 宮地 匡人

(56)【参考文献】

【文献】 特開２０１３−１８３５５１（ＪＰ，Ａ）

【文献】 特開２０１１−２５４６７８（ＪＰ，Ａ）

【文献】 特開２００２−３６６８７８（ＪＰ，Ａ）

【文献】 宮崎修司 他，磁場の運動方程式による磁場解析，電気学会研究会資料 静止器回転機合同研究会（ＳＡ・ＲＭ），２０１０年 ９月２９日，ＳＡ−１０−９２（ＲＭ−１０−１０１），ｐ．１−５

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １９／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、
前記磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、
前記磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、
を備え、
前記磁場演算部は、外部磁場の時間変動により各粒子に誘起される誘導磁化と、前記誘導磁化に基づいた磁気モーメント同士の相互作用により得られる磁場を用いて、誘導磁場を数値的に演算することを特徴とする解析装置。

【請求項2】

前記磁場演算部は、前記誘導磁化を表す項と、前記誘導磁化に基づいた磁気モーメントを有する粒子間の相互作用を表す磁場の補正項の双方を用いて、前記誘導磁場を演算することを特徴とする請求項１に記載の解析装置。

【請求項3】

仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、
前記磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、
前記磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、
を備え、
前記磁場演算部は、解析対象とする粒子群の重心位置を表す局所重心ベクトルを用いたゲージ変換により導出される誘導磁場の式を用いて、前記粒子群における誘導磁場を数値的に演算することを特徴とする解析装置。

【請求項4】

前記磁場演算部は、前記粒子系のベクトルポテンシャルに対して、前記局所重心ベクトルと前記粒子系の磁束密度の外積を用いて記述されるゲージに変換することにより導出される誘導磁場の式を用いて、前記粒子群における誘導磁場を数値的に演算することを特徴とする請求項３に記載の解析装置。

【請求項5】

前記粒子系には、電気的に互いが絶縁される第１粒子群および第２粒子群を含む複数の粒子群が配置可能とされており、
前記磁場演算部は、前記第１粒子群の重心位置を表す第１局所重心ベクトルを用いたゲージ変換により導出される第１の誘導磁場の式を用いて前記第１粒子群における誘導磁場を演算し、前記第２粒子群の重心位置を表す第２局所重心ベクトルを用いたゲージ変換により導出される第２の誘導磁場の式を用いて前記第２粒子群における誘導磁場を演算することを特徴とする請求項３または４に記載の解析装置。

【請求項6】

粒子系は繰り込み群分子動力学法を使用して繰り込まれた粒子系であることを特徴とする請求項１から５のいずれか一項に記載の解析装置。

【請求項7】

前記粒子状態演算部における支配方程式は磁気モーメントに依存する項を有することを特徴とする請求項１から６のいずれか一項に記載の解析装置。

【請求項8】

仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する機能と、
磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する機能と、
磁気的な物理量の演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する機能と、
をコンピュータに実現させ、
前記磁気的な物理量を演算する機能は、外部磁場の時間変動により各粒子に誘起される誘導磁化と、前記誘導磁化に基づいた磁気モーメント同士の相互作用により得られる磁場を用いて、誘導磁場を数値的に演算することを特徴とするコンピュータプログラム。

【請求項9】

仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する機能と、
磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する機能と、
磁気的な物理量の演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する機能と、
をコンピュータに実現させ、
前記磁気的な物理量を演算する機能は、解析対象とする粒子群の重心位置を表す局所重心ベクトルを用いたゲージ変換により導出される誘導磁場の式を用いて、前記粒子群における誘導磁場を数値的に演算することを特徴とするコンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

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

【0003】

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

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特開平１１−１４６６８８号公報

【特許文献2】特開２００９−３７３３４号公報

【非特許文献】

【0005】

【非特許文献1】Erik Bitzek他、「Structural Relaxation Made Simple」、Physical Review Letters、２００６年１０月２７日、９７

【非特許文献2】L.D.Landau、E.M.Lifshitz、L.P.Litaevskii著、「Electrodynamics of Continuous Media 2nd ed」、１９８４年１月１日、２０５ページ

【発明の概要】

【発明が解決しようとする課題】

【0006】

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

【0007】

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

【課題を解決するための手段】

【0008】

本発明のある態様の解析装置は、仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、を備える。磁場演算部は、外部磁場の時間変動により各粒子に誘起される誘導磁化と、誘導磁化に基づいた磁気モーメント同士の相互作用により得られる磁場を用いて、誘導磁場を数値的に演算する。

【0009】

本発明の別の態様もまた、解析装置である。この装置は、仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、を備える。磁場演算部は、解析対象とする粒子群の重心位置を表す局所重心ベクトルを用いたゲージ変換により導出される誘導磁場の式を用いて、粒子群における誘導磁場を数値的に演算する。

【0010】

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

【発明の効果】

【0011】

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

【図面の簡単な説明】

【0012】

【図1】実施の形態に係る解析装置の機能および構成を示すブロック図である。

【図2】図１の粒子データ保持部の一例を示すデータ構造図である。

【図3】図１の解析装置における一連の処理の一例を示すフローチャートである。

【図4】ｎ＝１次までの展開に対応する計算結果を示す図である。

【図5】ｎ＝２次までの展開に対応する計算結果を示す図である。

【図6】ｎ＝３次までの展開に対応する計算結果を示す図である。

【図7】ｎ＝４次までの展開に対応する計算結果を示す図である。

【図8】球面調和関数の次数と計算値との関係を示す図である。

【図9】ヒステリシスカーブを示す図である。

【図10】ビーズにより構成される導体球の解析モデルを示す図である。

【図11】導体球表面での時間変化する磁界の計算値を示すグラフである。

【図12】計算値の位相誤差と振幅誤差を示すグラフである。

【図13】導体球表面での磁界、電流および磁界の計算値を示すグラフである。

【発明を実施するための形態】

【0013】

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

【0014】

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

【0015】

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

【0016】

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

【0017】

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

【0018】

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

【0019】

さらに、本発明者は、粒子法における各粒子に付与される磁気モーメント同士の相互作用項に加えて各粒子に誘起される誘導磁化を考慮に入れることで、粒子法の適用範囲を電磁誘導現象を含む動磁場解析にまで拡張できることに想到した。これにより、誘導モータにおける固定子と回転子の相互作用といった、大変形を伴うとともに電磁誘導の考慮が必要な物理現象に対して、粒子法による動磁場解析を実現することができる。

【0020】

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

【0021】

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

【0022】

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

【0023】

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

【0024】

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

【0025】

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

【0026】

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

【0027】

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

【0028】

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

【0029】

磁場演算部１２１は、粒子データ保持部１１４によって保持されるデータが表す粒子系が生成する磁場を演算する。磁場演算部１２１は、粒子データ保持部１１４によって保持される各粒子の磁気モーメントに基づいて磁場を演算する。磁場は、粒子系に関連する磁気的な物理量である。磁場演算部１２１は磁場の代わりにまたは磁場に加えて磁束密度や磁化を演算してもよい。

【0030】

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

【0031】

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

【0032】

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

【0033】

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

【0034】

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

【0035】

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

【0036】

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

【0037】

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

【0038】

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

【0039】

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

【0040】

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

【0041】

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

【0042】

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

【0043】

また例えば解析対象がモータの場合、固定体に対して回転体が回転するので、やはりメッシュベースの解析手法は不向きであり静解析に止まる。そこで、本実施の形態に係る解析装置１００によると、回転体の回転を粒子法により高精度にシミュレーションしつつ、回転体と固定体との磁気的相互作用を動的に扱うことができる。さらに、誘導モータのように電磁誘導現象を含む解析対象であっても、電磁誘導現象を考慮に入れた誘導磁化と、誘導磁化に基づいた磁気モーメントを付与することにより、磁場の時間変化を加味した磁気的相互作用を扱うことができる。

【0044】

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

【0045】

磁場はマクロなベクトルポテンシャル：

【数1】

から以下の様に求まる。スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにある１個のスピンが作る磁場はスピンＳ_ｊｚ（磁気モーメントｍ_ｊｚ）を発生させる円電流を

【数2】

とし、

【数3】

ここで、

【数4】

ｎ＝０のみを残すと一様に磁化：

【数5】

した球の作る磁場と同等になる。
なお、真球のバルク材に本手法を適用する場合、ｎ＝０のみ採用すると、一様な磁化が得られず、磁気モーメント法と同様な結果になる。高次の項を加える事で改善できる。ＦＥＭは一様になる。

【0046】

デカルト座標系に変換（

【数6】

を使って、)すれば

【数7】

ここで、ｓｉｎθ_ｊ＝０のときは、Ｂ_θ＝０から、Ｂ_ｘ＝Ｂ_ｙ＝０であることに注意する。磁気モーメントｍ_ｘ、ｍ_ｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場

【数8】

も同様に求まる。

【数9】

【0047】

【数10】

における磁場

【数11】

（Ｎ個のスピンからの寄与）は

【数12】

ここで、

【数13】

は全外場を意味する。Ｐ−１４とＰ−１５を連立することで, 磁気モーメントｍ_ｐ が求まる。磁性体外部に生じる磁場は、Ｐ−１５であり、磁性体内部に生じる磁場は、

【数14】

Ｐ−１４は球体表面に生じる磁化による反磁場が考慮され、

【数15】

を満足する。

【0048】

Ｐ−１３はマトリクスが優対角でない(ｊ＝ｐが除外されている) から収束は遅い。そこで、以下の運動方程式から解を得る。

【数16】

ここで、ｍ_ｖ＝５×１０^−２ｄｔ^２は仮想質量、γ＝ｍ_ｖ／（１０ｄｔ）は、減衰係数である。

【0049】

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

【数17】

【0050】

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

【数18】

【0051】

以下に、ｘ軸方向にＢ_ｏｘ＝０．１［Ｔ］の一様な外場においた磁性球（帯磁率が９９９）の磁化

【数19】

を計算した結果を示す。球を構成する粒子（原子）数は４００９個とし、ｆｃｃ構造に配置した。厳密解は

【数20】

である。図４、図５、図６、図７に多重極展開の１次〜４次までの結果を示す。共に中心を通過する断面を表示した。多重極展開の次数を上げると, 厳密解と良く合う事が判る。次数が低いと収束が悪い。また計算時間の大半は球面調和関数の呼び出しである。

【0052】

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

【0053】

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

【0054】

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

【0055】

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

【0056】

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

【0057】

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

【0058】

【数21】

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

【0059】

スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにあるスピンＳ_ｊｚ（磁気モーメントｍ_ｊｚ）を発生させる円電流は

【数22】

であるから、磁場はマクロなベクトルポテンシャル：

【数23】

から以下の様に求まる。

【数24】

ここで、

【数25】

ｎ＝０のみを残すと一様に磁化：

【数26】

した球の作る磁場と同等になる。

【0060】

デカルト座標系に変換(

【数27】

を使って、）すれば

【数28】

ここで、ｓｉｎθ_ｊ＝０のときは、Ｈ_θ＝０から、Ｈ_ｘ＝Ｈ_ｙ＝０であることに注意する。磁気モーメントｍ_ｘ、ｍ_ｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場

【数29】

も同様に求まる。

【数30】

【0061】

【数31】

におけるモーメントが感じる外磁場

【数32】

(Ｎ個のスピンからの寄与) は

【数33】

誘導磁場

【数34】

は、ｐ点が導体内部である場合、

【数35】

ｐ点が導体外部である場合、

【数36】

Ｃ_１、Ｃ_２は、例えば、非特許文献２に示される式において

【数37】

となる極限と比較したときの補償係数であり、Ｃ_１＝１／３、Ｃ_２＝３／５の値をとる。

【0062】

ここで、Ｑ−２１の右辺第１項（またはＱ−２１−２）は、時間変化する外部磁場により各粒子に誘起される誘導磁化Ｍ_ａｎｔを表す項である。また、Ｑ−２１の右辺第２項およびＱ−２２の右辺は、Ｑ−２１−４式で表される各粒子の誘導磁化に基づいた磁気モーメントｍ_ａｎｔ同士の相互作用により得られる磁場を表す項である。誘導磁場をこのように表すことで、時間変化する外部磁場による影響を考慮した動磁場解析を実現することができる。

【0063】

なお、繰り込みに際し、電気伝導度σは

【数38】

であり、今はδ＝２としている。

【0064】

Ｑ−２１−２右辺のLagrange 微分

【数39】

は、粒子法であるから

【数40】

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

【0065】

磁性体内部に生じる磁場は、

【数41】

であり、全外場Ｈ_Ｏ（ｘ_ｐ）と、磁化Ｍ（Ｈ）との関係は、

【0066】

【数42】

である。Ｑ−２６とＱ−２７を連立することで、磁気モーメント

【数43】

が求まる。また、磁性体外部の磁場は、

【数44】

となる。

【0067】

非線形磁化とヒステリシス
Ｑ−２７とＱ−２４から磁場

【数45】

を消去する：

【数46】

Ｑ−２８を解くことで

【数47】

が得られる（図９にはＢ＝ｆ（Ｈ）の場合を示す）。

【0068】

図９は、ヒステリシスカーブを示す図である。ヒステリシスカーブとＢ＋２μ_ｏＨ＝３Ｂ_ｏの交点からＭを求める。今

【数48】

を永久磁化とするＢ−Ｈカーブ：

【数49】

が与えられたとする。

【数50】

と連立し

【数51】

を得る。一般のＢ＝ｆ（Ｈ）に対しては、Newton-Laphson法で解く。

【0069】

磁場の粒子化法
Ｑ−１９はマトリクスが優対角でない（ｊ＝ｐが除外されている）から収束は遅い。そこで、以下の運動方程式から解を得る。

【数52】

ここで、Ｈ_Ｏ（ｘ_ｐ）は、Ｑ−２６により表される。また、ｍ_ｖは仮想質量である。推奨値は２５０×ｄｔ×ｄｔである。ｍ_ｖを十分小さく取れば、正解の周りの微小な減衰振動が起こり、誤差を小さくできる。

【0070】

模型を考える：

【数53】

減衰項を持つ強制振動の解(Landau and Lifshitz,Mechanics,Sec.26) は

【数54】

【数55】

を考えれば、

【数56】

と成る。即ちｍ_ｖを十分小さく取れば、解くべき本来の方程式

【数57】

の解に一致する。

【0071】

離散化
蛙飛び法（速度はｎとｎ＋１の中間点における評価）に従い離散化すれば、収束回数をｓと書き、

【数58】

ここで、誘導磁場Ｈ_ｉｎｄ（ｘ_ｉ）の離散化において、

【数59】

ととる。このとき、対角要素のみ未来を使えば、陽的解法に出来、約３０％程度効率改善になる。計算の大半は特殊関数の呼び出しである。

【0072】

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

【数60】

【0073】

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

【0074】

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

【0075】

図１２は、計算値の位相誤差と振幅誤差を示すグラフであり、ｘ軸上に位置する粒子の重心位置における磁界の位相誤差ε_{ｐｈａｓｅ}と振幅誤差ε_{ａｍｐｌｉｔｕｄｅ}を示す。位相誤差ε_{ｐｈａｓｅ}および振幅誤差ε_{ａｍｐｌｉｔｕｄｅ}は、

【数61】

と定義した。ここで、θ_{ｅｘａｃｔ}は厳密解の位相、θ_ＭＢＭは解析結果の位相である。また、Ｈ_{ｅｘａcｔ}は厳密解により得られた磁界、Ｈ_ＭＢＭは解析結果の磁界である。なお位相誤差は、Ｈ_{ｅｘａcｔ}が最大となる位相θ_{ｅｘａｃｔ}と、Ｈ_ＭＢＭが最大となる位相θ_ＭＢＭとを比較した。図１２より、ｘ軸上における解析結果と厳密解の磁界の位相誤差は最大で５．８６％、振幅誤差は６．３８％となり、厳密解と比較して精度の高い解析結果が得られることが示された。

【0076】

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

【数62】

【0077】

粒子に働く力の導出
粒子（電子と核）と、粒子磁場間の相互作用を記述するLagrangeanは

【数63】

である。ベクトルポテンシャル

【数64】

と電磁場の関係は、

【数65】

である。

【数66】

を使い、以下の運動方程式を得る。

【数67】

原子Ｍ＝Ｍ_ｐ＋ｍ_ｅに対する運動方程式は、

【数68】

となる。右辺：Lorentz力を平均電流密度

【数69】

、電気モーメントを

【数70】

で置き換えると

【数71】

Maxwellの方程式とＱ−２４から

【数72】

最終的に原子に働く力は、電気分極を無視し、

【数73】

となる。

【0078】

粒子に働く力は、直接Maxwellの応力テンソルからも求まる。

【数74】

球に局在した電流が作るモーメント

【数75】

に作用する力は、最低の次数では

【数76】

である。従い、ポテンシャルエネルギーは、電気分極によるポテンシャルエネルギーを合わせて、

【数77】

となるが、渦電流や真電流を包含しない。

【0079】

【数78】

結局、

【数79】

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

【0080】

磁化による電流

【数80】

と渦電流

【数81】

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

【数82】

【0081】

微係数の計算
以下に必要な微分係数を書き出す。

【数83】

等からｒ_ｐｊ＞ａのみ列挙すると（なお、

【数84】

はｓｉｎθに比例するから、θ＝０で発散は無い）、

【数85】

【数86】

【0082】

磁気モーメントが作る磁化電流の計算
ｚ軸に平行に磁化したとき
今一度、ｚ軸に磁化したときの磁場を書く

【数87】

磁化電流を

【数88】

と書けば、磁化により生じる磁力

【数89】

は

【数90】

である。先ず、

【数91】

を求める。

【数92】

【0083】

ｘ軸に平行に磁化したとき
ｘ軸に磁化したときの磁場を書く：

【数93】

【数94】

磁化電流

【数95】

を求める。

【数96】

【0084】

ｙ軸に平行に磁化したとき
ｙ軸に磁化したときの磁場を書く：

【数97】

磁化電流

【数98】

を求める。

【数99】

【0085】

誘導磁界が作る渦電流の計算
導体内部の誘導磁場は、Ｑ−２１より、

【数100】

となる。渦電流は

【数101】

なお、上式Ｒ−５２−２は、Ｒ−１２において、球面調和関数の展開次数ｎ＝０とした場合と同じである。

【0086】

粒子に働く力の計算
ｍ_ｘ、ｍ_ｙ、ｍ_ｚと渦電流による寄与を加えれば、全電流が得られる。

【数102】

座標

【数103】

に粒子に働く力は

【数104】

最後に

【数105】

の微分から出てくる項を付加する必要がある。これは磁場の誤差

【数106】

から生じる見かけの力であり、運動と磁場の全エネルギーを保存させる。

【数107】

【0087】

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

【0088】

ベクトルポテンシャルＡについて、以下の微分方程式が成立する。

【数108】

上記微分方程式の解は、

【数109】

となる。ここで、

【数110】

として、各粒子に分割すれば、各粒子の位置ｒ_ｉでのベクトルポテンシャルは、

【数111】

となる。ここで、ｖ_ｊは各粒子の体積である。なお、ｊの和に際し、ベクトル

【数112】

が粒子表面を横切る場合には和を取らない。また、粒子径ａが満たすべき条件は、

【数113】

である。

【0089】

このとき、各対象物を粒子に分割することにより、ベクトルポテンシャルＡの並進対称性を失うおそれがある。そこで、本発明者は、下記Ｓ−４に示されるゲージ変換を導入することにより、ベクトルポテンシャルＡの並進対称性を回復できることを見いだした。なお、このゲージ変換を本明細書において「ビーズゲージ」ともいう。

【数114】

ここで、ベクトルｇ_ｉは、計算対象とする粒子群（粒子数：Ｎ_ｉ）に関する局所重心ベクトルであり、

【数115】

で表される。したがって、ベクトルｇ_ｉは、粒子群の重心位置を表している。また、Ｓ−４に示すビーズゲージは、磁束密度Ｂ_ｊと局所重心ベクトルｇ_ｉの外積を用いて記述されると言える。

【0090】

上記のビーズゲージについて、磁束密度ＢとベクトルポテンシャルＡの関係式

【数116】

を求めれば、

【数117】

となる。したがって、粒子数Ｎ_ｉが十分に大きければ、並進対称性が満たされることがわかる。

【0091】

ここで、Ｓ−４に示すゲージ変換から、Ｓ−３の右辺に含まれるベクトルポテンシャルの時間微分

【数118】

を磁束密度Ｂの時間微分を用いて以下のように表すことができる。

【数119】

さらに、ベクトルｄ_ｊｉを、

【数120】

とすると、ベクトルポテンシャルの時間微分は、

【数121】

と表される。したがって、Ｓ−３にＳ−９を代入して積分を実行すれば、ベクトルポテンシャルＡの厳密解を磁束密度Ｂの時間微分を用いて記述できる。

【0092】

ベクトルポテンシャルＡの厳密解を導出するため、まず、磁束密度Ｂのｚ成分の時間微分である、

【数122】

のみが作用する場合の積分を実行する。これにより、ベクトルポテンシャルＡは、

【数123】

となる。磁束密度Ｂのｘ成分およびｙ成分の時間微分である、

【数124】

の寄与も加えれば、

【数125】

が得られる。

【0093】

解析対象とする導体内部を記述する誘導磁場Ｈは、

【数126】

の関係式より、

【数127】

となる。ここで、

【数128】

と置換すると、誘導磁場Ｈの式として、

【数129】

が得られる。

【0094】

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

【0095】

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

【0096】

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

【0097】

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

【0098】

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

【0099】

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

【0100】

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

【符号の説明】

【0101】

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

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】