JP5915157B2 - 磁性体特性解析プログラム、磁性体特性解析装置、及び磁性体特性解析方法 - Google Patents

Description

本発明は、磁性体特性解析プログラム、磁性体特性解析装置、及び磁性体特性解析方法に関する。

従来、磁性体を用いた機器（モータ、トランス等）の特性や損失を解析するために、コンピュータを用いた磁性体の特性解析シミュレーションが行われている。係るシミュレーションにおいては、有限要素法による磁場解析技術が活用されている。有限要素法による磁場解析技術を用いると、解析対象を小さなメッシュに分割することで複雑な形状を扱うことができる。

また、磁性体の特性解析に用いられる他の手法として、マイクロマグネティクスが知られている。マイクロマグネティクスを用いると、コンピュータが高精度に磁性体の磁区構造を再現することができる。

磁性体の特性解析に関連し、解析対象物から発生する電流ベクトルによる磁界を算出するための第一の方程式を有限要素法および境界積分法により解析可能な形式で作成する磁界磁性体特性解析装置についての発明が開示されている（例えば、特許文献１参照）。特許文献１に記載の磁界磁性体特性解析装置は、解析対象物から発生する磁化ベクトルによる磁界を算出するための第二の方程式を有限要素法および境界積分法により解析可能な形式で作成する。そして、特許文献１に記載の磁界磁性体特性解析装置は、第一の方程式と第二の方程式とを用いて、第一の磁界および第二の磁界それぞれを算出し、第一と第二の磁界を合算した合算磁界を解析対象物における磁界の解析結果とする。

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

しかしながら、上記特許文献に記載の装置の如く、有限要素法と他の手法を用いて磁性体の特性解析を行う場合、有限要素法に基づく演算の処理負荷が高く、全体として処理速度が低下するという課題が生じる。

一つの側面では、本発明は、磁性体の特性の解析スピードの向上を図ることを目的とする。

上記目的を達成するための一態様は、
磁性体特性解析プログラムにおいて、
磁性体特性解析装置に、
解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行わせ、
前記有限要素法による磁場解析により算出される磁界を固定値として前記各要素に作用する有効磁界を算出させ、
前記算出した有効磁界を用いてＬＬＧ（Landau Lifshitz Gilbert）方程式を時間積分して前記各要素内で磁化ベクトルを算出させ、
前記算出した磁化ベクトルを平均して前記各要素についての平均磁化を算出させることを特徴とする磁性体特性解析プログラムである。

一実施態様によれば、磁性体の特性の解析スピードの向上を図ることができる。

本発明の一実施例に係る磁性体特性解析装置１のハードウエア構成例である。 本発明の第１実施例に係る磁性体特性解析装置１の機能構成例である。 有限要素法においてメッシュへの分割が行われる様子を示す図である。 複数節点で構成される三角形要素と四角形要素を示す図である。 一次元の補間関数Ｎを一次元座標上で表現した図である。 立方体を構成する複数の磁化ベクトルを示す図である。 磁化ベクトルの平均値と外部磁界の関係を示す図である。 第１演算処理部４０と第２演算処理部５０により実行される処理を概念的に示す図である。 制御部６０により実行される処理の流れを示すフローチャートである。 本実施例の磁性体特性解析装置１の解析結果を表示装置２４により表示した画面例である。 参考手法において実行される処理を概念的に示す図である。 参考手法を実行するための処理の流れを示すフローチャートである。 本実施例の手法による解析結果と、参考手法による解析結果を比較した図である。 本実施例の手法と参考手法におけるループ（Ｂ）の実行回数を比較した図である。 本実施例の手法と参考手法の処理時間を比較した図である。 解析対象となった複数の磁性体モデルを示す図である。 磁性体特性解析装置１に与えられた各磁性体モデルの寸法を示す図である。 異なる形状を有する複数の磁性体モデルについて、本実施例の手法により特性解析シミュレーションを行った結果を示す図である。

以下、本発明を実施するための形態について、添付図面を参照しながら実施例を挙げて説明する。

以下、図面を参照し、本発明の一実施例に係る磁性体特性解析プログラム、磁性体特性解析装置、及び磁性体特性解析方法について説明する。

［構成］
図１は、本発明の一実施例に係る磁性体特性解析装置１のハードウエア構成例である。磁性体特性解析装置１は、例えば、ＣＰＵ（Central Processing Unit）１０と、ドライブ装置１２と、補助記憶装置１６と、メモリ装置１８と、インターフェース装置２０と、入力装置２２と、表示装置２４と、を備える。これらの構成要素は、バスやシリアル回線等を介して接続されている。

ＣＰＵ１０は、例えば、プログラムカウンタや命令デコーダ、各種演算器、ＬＳＵ（Load Store Unit）、汎用レジスタ等を有する演算処理装置としてのプロセッサである。

ドライブ装置１２は、記憶媒体１４からプログラムやデータを読み込み可能な装置である。プログラムを記録した記憶媒体１４がドライブ装置１２に装着されると、プログラムが記憶媒体１４からドライブ装置１２を介して補助記憶装置１６にインストールされる。記憶媒体１４は、例えば、ＣＤ（Compact Disc）、ＤＶＤ（Digital Versatile Disc）、ＵＳＢ（Universal Serial Bus）メモリ等の可搬型の記憶媒体である。また、補助記憶装置１６は、例えば、ＨＤＤ（Hard Disk Drive）やフラッシュメモリである。

プログラムのインストールは、上記のように記憶媒体１４を用いる他、インターフェース装置２０がネットワークを介して他のコンピュータよりダウンロードし、補助記憶装置１６にインストールすることによって行うこともできる。ネットワークは、インターネット、ＬＡＮ（Local Area Network）、無線ネットワーク等である。また、プログラムは、磁性体特性解析装置１の出荷時に、予め補助記憶装置１６やＲＯＭ（Read Only Memory）等に格納されていてもよい。

このようにしてインストール又は予め格納されたプログラムをＣＰＵ１０が実行することにより、図１に示す態様の情報処理装置を、本実施例の磁性体特性解析装置１として機能させることができる。

メモリ装置１８は、例えば、ＲＡＭ（Random Access Memory）、ＥＥＰＲＯＭ（Electrically Erasable and Programmable Read Only Memory）やフラッシュメモリ等である。インターフェース装置２０は、上記ネットワークとの接続等を制御する。

入力装置２２は、例えば、キーボード、マウス、ボタン、タッチパッド、タッチパネル、マイク等である。また、表示装置２４は、例えば、ＬＣＤ（Liquid Crystal Display）やＣＲＴ（Cathode Ray Tube）等の表示装置である。磁性体特性解析装置１は、表示装置２４の他、プリンタ、スピーカ等の他の種類の出力装置を備えてもよい。

図２は、本発明の第１実施例に係る磁性体特性解析装置１の機能構成例である。磁性体特性解析装置１は、設定入力受付部３０と、第１演算処理部４０と、第２演算処理部５０と、制御部６０と、を備える。第２演算処理部５０は、更に、有効磁界算出５１と、磁化ベクトル算出５２と、平均磁化算出部５３と、を有する。これらの機能ブロックは、補助記憶装置１６等に格納されたプログラム・ソフトウエアをＣＰＵ１０が実行することにより機能する。なお、これらの機能ブロックが明確に分離したプログラムによって実現される必要はなく、サブルーチンや関数として他のプログラムによって呼び出されるものであってもよい。また、機能ブロックの一部が、ＬＳＩ（Large Scale Integrated circuit）、ＩＣ（Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）等のハードウエア手段であっても構わない。

設定入力受付部３０は、各種パラメータ設定、及び解析対象の物性情報等の入力を受け付ける。各種パラメータには、磁場の非定常計算ステップＮ、ＬＬＧループ回数Ｍ、収束判定閾値ε、自由度Ｐ、外部磁界の時間当たりの変動量等が含まれる。また、解析対象の物性情報には、基本物性値、材料の性質、異方性の設定等が含まれる。これらの設定入力は、例えばインターフェース装置２０や入力装置２２に対して行われ、入力されたデータは補助記憶装置１６やメモリ装置１８に格納される。

第１演算処理部４０は、有限要素法による磁場解析を行う。また、第２演算処理部５０は、マイクロマグネティックスに基づく磁性体の特性解析を行う。以下、有限要素法とマイクロマグネティックスについて説明する。

［有限要素法とマイクロマグネティックス］
（有限要素法の基本概念）
有限要素法は、解析対象を有限個の小領域（要素）であるメッシュに分割し、解を求める手法である。流体、構造、磁場等の物理状態を支配する微分方程式は一般に支配方程式と呼ばれ、空気や固体・液体等の物質中に連続的に定義される物理量（連続場）が未知変数となる。有限要素法は、解析的に解くことが難しい微分方程式の近似解を数値的に得る方法の一つである。有限要素法では、本来は無限の自由度を持つ支配方程式の定義領域Ωを有限個の小領域（要素）にメッシュに分割することにより、解を求めることが可能になる。図３は、有限要素法においてメッシュへの分割が行われる様子を示す図である。

要素とは、解析対象となる領域を分割した最小単位の領域であり、複数の節点によって構成される。図４は、複数節点で構成される三角形要素と四角形要素を示す図である。このように、有限要素法では、モデル全体を小さな要素で構成して解を求めるため、複雑な形状の解析対象領域を扱うことができる。

有限要素法による磁場解析手法として、ガラーキン法（または重み付き残差法）と称される手法が用いられる。ガラーキン法とは、支配方程式の両辺に任意形状の重み関数を乗じ、空間積分することで結果的に自由度を低減（離散化）する手法である。

式（１）は、物理変数φの空間分布を表現する支配方程式である。ここで、φは空間Ｘに依存する。有限要素法では、支配方程式が非線形であっても、解を求める過程では式（２）で示される線形方程式として扱うことができる。式（２）中、Ｆはφに対する空間微分の演算子の関数、Ｆ₀はφに依存しない定数である。

Ｎ個の節点を持つ要素内の解φは、補間関数Ｎを用いて式（３）のように近似することができる。式（３）中、φ_iは、位置ｘ_iにおける場（フィールド）の値、ｎは節点の数である。また、重み関数ｗ（ｘ）は、補間関数Ｎを用いて式（４）のように定義される。ここで、位置ｘ_iにおける補間関数Ｎの値は、式（５）で示される性質を有する。図５は、一次元の補間関数Ｎを一次元座標上で表現した図である。

次に、補間関数Ｎを用いて、二次元空間で定義された式を有限要素法で解く手法について説明する。まず、式（２）の両辺に任意の重み関数ｗを乗じ、面積で積分すると式（６）となる。そして、重み関数ｗと求める解φを、補間関数Ｎを用いて表現すると、一要素当たりの式は式（７）となる。

このように空間積分を行うことによって、空間の自由度が消滅し、結果的に未知数φに関する方程式となる。そして、任意の重み関数の数ｎだけφ_jに関する方程式が得られるので、この連立方程式を解くことで解φを求めることが可能となる。

（有限要素法の磁場解析への適用）
空気中に磁化ベクトル→Ｍが存在する場合の磁場の方程式は、式（８）で表される。以下、数式等においてアルファベット上に矢印が付されたベクトル表記を、明細書文中では「→Ａ」等と表記する。また、＜＞は平均を表すものとする。式（８）中、→Ａはベクトルポテンシャル、ｖ₀は真空の透磁率の逆数、→Ｊ₀は励磁電流、→Ｍは磁化ベクトルである。

式（８）の両辺に重み関数ｗを乗算し、一つの要素について空間積分を行うと、式（９）が得られる。

更に、式（９）の左辺第二項について部分積分を行うと、式（１０）が得られる。

なお、式（１０）の右辺第二項を導出する過程で関係式（１１）を用いている。また、磁化ベクトル→Ｍに起因する項は、部分積分によって式（１２）に書き換えられる。

式（９）の空間積分によって得られる係数ｃ_ij、ｄ_ij、ｅ_i、ｆ_i、ｇ_iを、次式（１３）〜（１７）で表す。

そして、式（９）を任意の値の重みｗ_iに対する恒等式とみなすと、重みｗ_iの係数はゼロとなり、連立方程式（１８）が得られる。

式（１８）は、一つの要素における未知数Ａ_jに関する連立方程式である。全ての要素に対して未知数Ａ_jに関する連立方程式を合成することで、定義領域全体の連立方程式が構成できる。ここで、要素単位の連立方程式である式（１８）のインデックスｉ、ｊをそれぞれ、領域全体のインデックスＩ、Ｊに書き換え、時刻ｎ、ｎ＋１に対して陰解法で記述すると式（１９）となる。式（１９）における未知数Ａ_j ⁿ⁺¹を左辺に移行すると式（２０）が得られる。

そして、式（２１）のように、未知数Ａ_j ⁿ⁺¹に関する連立方程式を解くことにより、時刻ｎ＋１における解Ａ_j ⁿ⁺¹を求めることができる。なお、解析領域の境界から磁束が流入又は流出しない一般的な場合には、式（２１）における→Ｂがゼロとなり、式（２２）が成立する。

こうして解として得られるベクトルポテンシャルＡ_j ⁿ⁺¹により、式（２３）のように磁束密度Ｂⁿ⁺¹が得られる。ここで、Ｎは一つの要素内に定義される補間関数であり、Ｂⁿ⁺¹は対象となる要素における値である。

そして、式（２３）の磁束密度Ｂⁿ⁺¹を用いて、式（２４）の磁界Ｈ_j ⁿ⁺¹が得られる。

（マイクロマグネティクスをベースとする磁気特性の表現方法）
本実施例の第２演算処理部５０が有する磁化ベクトル算出部５２は、マイクロマグネティクスにおけるＬＬＧ方程式を積分し、ＬＬＧ方程式が所定値未満となる定常状態の磁化ベクトル→Ｍを求める。ＬＬＧ方程式は、式（２５）で表される。ＬＬＧ方程式は、磁性体の磁性体の磁区構造やそれに基づく磁気特性を解析するために用いられる。式（２５）中、→Ｍ_iは磁化ベクトル、→Ｈ_effは磁化ベクトルに作用する有効磁界、γはジャイロ磁気定数（Hz・m／A）、αはダンピングコンスタントである。

磁化ベクトルに作用する有効磁界は、式（２６）で表される。式（２６）中、→Ｈ_aniは結晶磁気異方性磁界、→Ｈ_excは磁気交換結合磁界、→Ｈ_externalは外部磁界、→Ｈはその他の磁界（磁界による静磁界、励磁電流、渦電流によるもの）である。

ＬＬＧ方程式は、複数の部分構造を有する磁性体モデルの磁気特性や磁区構造を解析するために活用されてきた。ＬＬＧ方程式を解くことで得られる個々の部分構造の磁化ベクトルは、それぞれ空間上の位置に配置され、磁化ベクトルの空間分布による磁界（静磁界、交換結合磁界）も磁気特性に寄与する。しかし、磁化ベクトルの空間分布に依存する有効磁界（静磁界、交換結合磁界）の影響を減らしても、磁性体の特性を近似的に再現することが可能である。従って、詳細な幾何形状を無視した適当なサイズの直方体モデルに対してＬＬＧ方程式を適用して磁化ベクトルを算出し、磁化ベクトルの平均値を求めることによって、磁性体の特性を近似的に表現することができる。

より具体的には、有限要素法によりメッシュに分割された要素を更に分割した内部構造について磁化ベクトルを求め、内部構造毎の磁化ベクトルの平均値を、当該要素の磁化ベクトル平均値とする。

図６は、立方体を構成する部分構造のそれぞれについて算出される磁化ベクトルを示す図である。図中、矢印が磁化ベクトルの向きを示している。

磁化ベクトルの数（自由度）Ｐを、図６で示すようにＰ＝ｌ×ｍ×ｎとすると、磁化ベクトル平均値（平均磁化）は式（２７）で表わされる。

図７は、磁化ベクトルの平均値と外部磁界の関係を示す図である。図示するように、磁化ベクトルの平均値を用いた結果は、磁気ヒステリシスを正確に再現することができる。

ところで、ＬＬＧ方程式を磁性体に適用した場合、扱える磁性体の大きさは数百ミクロンが限界である。従って、本実施例の磁性体特性解析装置１では、有限要素法を用いてメッシュへの分割を行った各要素に対してＬＬＧ方程式を適用する。これによって、解析対象の形状やサイズに拘わらず、ＬＬＧ方程式を用いた磁気特性の解析を行うことができる。

［具体的動作］
図８は、第１演算処理部４０と第２演算処理部５０により実行される処理を概念的に示す図である。また、図９は、制御部６０により実行される処理の流れを示すフローチャートである。

図９に示すように、まず、制御部６０は、パラメータの初期化を行う（Ｓ１００）。パラメータの初期化には、磁場の非定常計算ステップＮ及びＬＬＧループ回数Ｍの指定、収束判定閾値εの設定、磁場の方程式の時間ステップｎの初期化等が含まれる。制御部６０は、設定入力受付部３０が受け付けた値、或いはデフォルト値として補助記憶装置１６等に格納された値を読み込み、メモリ装置１８上に展開する。

次に、制御部６０は、物性情報を設定し、外部磁界→Ｈ_externalの初期値を設定する（Ｓ１０２）。物性情報や外部磁界→Ｈ_externalは、上記パラメータと同様に、設定入力受付部３０が受け付けた値、或いはデフォルト値として補助記憶装置１６等に格納された値であり、メモリ装置１８上に展開される。

次に、制御部６０は、第１演算処理部４０に有限要素法におけるメッシュへの分割を行わせ、複数の要素を設定させる（Ｓ１０４）。第１演算処理部４０は、解析対象となる磁性体を含む空間を、複数の要素であるメッシュに分割する。

次に、制御部６０は、第２演算処理部５０の磁化ベクトル算出部５２がＬＬＧ方程式の計算に用いる自由度（各要素内で、磁化ベクトル→Ｍ_iを計算する内部構造の数）Ｐを設定する（Ｓ１０６）。

次に、制御部６０は、分割した各要素内の内部構造について、初期の磁化ベクトルを例えばランダムに設定する（Ｓ１０８）。そして、制御部６０は、設定した磁化ベクトルに基づいて各要素の平均磁化＜→Ｍ＞を算出する（Ｓ１１０）。なお、Ｓ１０６〜Ｓ１１０の処理は、第２演算処理部５０に行わせてもよい。

このように処理状態の設定が終了すると、制御部６０は、第１演算処理部４０に有限要素法の磁場解析を行わせる（Ｓ１１２）。

第１演算処理部４０は、まず、Ｓ１１０又は後述するＳ１２０において算出された平均磁化＜→Ｍ＞を図８における式（２８）に代入し、ベクトルポテンシャル→Ａ_j ⁿ⁺¹を算出する。次に、第１演算処理部４０は、ベクトルポテンシャル→Ａ_j ⁿ⁺¹に基づき式（２９）により磁束密度→Ｂⁿ⁺¹を算出し、式（３０）により磁界→Ｈⁿ⁺¹を算出する。ここで、磁界→Ｈⁿ⁺¹は、磁束密度→Ｂⁿ⁺¹を真空透磁率μ₀で除算して得られる成分→Ｈ₁ ⁿ⁺¹と、平均磁化＜Ｍ＞の符号を逆にした値を真空透磁率μ₀で除算して得られる成分→Ｈ₂ ⁿ⁺¹に分解することができる。第１演算処理部４０は、これらのうち磁束密度→Ｂⁿ⁺¹を真空透磁率μ₀で除算して得られる成分→Ｈ₁ ⁿ⁺¹を第２演算処理部５０の有効磁界算出部５１に出力する。

第１演算処理部４０による有限要素法の磁場解析は、分割された各要素を含む全空間について行われる。従って、磁場解析の結果には、ある要素における平均磁化＜→Ｍ＞に起因して他の要素に作用する磁界等が全て反映されることになる。この結果、ＬＬＧ方程式の時間ステップ方程式による磁性体の特性解析が、一定以上の大きさを有する磁性体デバイスの特性解析まで拡張される。

第１演算処理部４０による有限要素法の磁場解析が完了すると、制御部６０は、第２演算処理部５０に以下の処理を実行させる。第２演算処理部５０の有効磁界算出部５１は、図８における式（３１）により有効磁界を算出する（Ｓ１１６）。

ここで、有効磁界算出部５１は、式（３１）における磁界の項→Ｈ₁ ⁿ⁺¹については、Ｓ１１２において第１演算処理部４０により算出された値を固定値として用いるが、他の成分＜→Ｍ＞／μ₀については、Ｓ１２０で更新した値を再帰的に反映させる。

次に、第２演算処理部５０の磁化ベクトル算出部５２は、各要素の内部構造について、図８における式（３２）に示すＬＬＧ方程式を設定し、ＬＬＧ方程式を時間積分して磁化ベクトル→Ｍ_iを算出する（Ｓ１１８）。

そして、第２演算処理部５０の平均磁化算出部５３は、各要素について、図８における式（３３）により、磁化ベクトル→Ｍ_iの平均値である平均磁化＜→Ｍ＞を求めて平均磁化＜→Ｍ＞を更新し、ＬＬＧ方程式の時間ステップｋを１増加させる（Ｓ１２０）。

次に、第２演算処理部５０は、ＬＬＧ方程式の時間ステップｋがＬＬＧループ回数Ｍに等しくなったか否かを判定する（Ｓ１２２）。ＬＬＧループ回数Ｍは、例えば数十〜数百程度の値に設定される。

ＬＬＧ方程式の時間ステップｋがＬＬＧループ回数Ｍに等しくない場合、第２演算処理部５０は、Ｓ１１６に戻り処理を実行する（ループ（Ａ））。

ＬＬＧ方程式の時間ステップｋがＬＬＧループ回数Ｍに等しくなると、第２演算処理部５０の処理が一時的に終了し、制御部６０が平均磁化＜→Ｍ＞の残差δ＜Ｍ＞を算出する（Ｓ１２４）。残差δ＜Ｍ＞は、Ｓ１１２〜Ｓ１２０の処理において最終的に算出された平均磁化＜→Ｍ＞から、前回実行されたＳ１１２〜Ｓ１２０の処理において最終的に算出された平均磁化＜→Ｍ＞を差し引いた値に基づく指標値である。残差δ＜Ｍ＞は、例えば分割された各要素について算出された平均磁化＜→Ｍ＞の平均値＜＜→Ｍ＞＞の差分である。また、残差δ＜Ｍ＞は、これに限らず、各要素についての平均磁化＜→Ｍ＞の差分の集合であっても構わない。

残差δ＜Ｍ＞は、例えば式（３４）又は式（３５）により算出される。式（３５）の方が高精度に残差δ＜Ｍ＞が算出されるが、処理負荷が大きいため、磁性体特性解析装置１は、解析に求められる精度に応じていずれか任意の式を採用してもよいし、ユーザ設定等によって切り換えてもよい。

残差δ＜Ｍ＞を算出すると、制御部６０は、残差δ＜Ｍ＞が収束判定閾値ε未満であるか否かを判定する（Ｓ１２６）。本ステップの判定は、残差δ＜Ｍ＞が平均磁化＜→Ｍ＞の平均値＜＜→Ｍ＞＞の差分である場合は、当該平均値＜＜→Ｍ＞＞の差分が収束判定閾値ε未満であるか否かを判定するものとなる。また、残差δ＜Ｍ＞が各要素についての平均磁化＜→Ｍ＞の差分の集合である場合は、例えば差分の全てが収束判定閾値ε未満であるかを判定するものとなる。

残差δ＜Ｍ＞が収束判定閾値ε以上である場合、Ｓ１１２に戻り、第１演算処理部４０が有限要素法の磁場解析を行う（ループ（Ｂ））。

残差δ＜Ｍ＞が収束判定閾値ε未満である場合、現在の外部磁界→Ｈ_externalに関して平均磁化＜→Ｍ＞が十分に収束したと判断することができる。従って、磁性体特性解析装置１は、外部磁界→Ｈ_externalと、平均磁化＜→Ｍ＞の規格化した値（例えば各要素について平均した値）を、解析結果の一部として出力する（Ｓ１２８）。係る解析結果は、例えばメモリ装置１８に格納される。

次に、制御部６０は、磁場の方程式の時間ステップｎが磁場の非定常計算ステップＮに等しくなったか否かを判定する（Ｓ１３０）。

磁場の方程式の時間ステップｎが磁場の非定常計算ステップＮに等しくない場合、制御部６０は、Ａⁿ⁺¹＝Ａⁿ、ｎ＝ｎ＋１とし、時間ステップを進める（Ｓ１３２）。この際に、制御部６０は、予め定められた外部磁界の時間当たりの変動量等に基づき、外部磁界→Ｈ_externalを一定量変更する。そして、Ｓ１１２に戻り、第１演算処理部４０が有限要素法の磁場解析を行う（ループ（Ｃ））。

ループ（Ｃ）に係る処理を繰り返し実行することによって、外部磁界→Ｈ_externalと、当該外部磁界が与えられた際の平均磁化＜→Ｍ＞の収束値を規格化した値との組み合わせであるデータが複数個出力されることになる。係る複数個のデータが、本実施例の磁性体特性解析装置１の解析結果の一つとなる。図１０は、本実施例の磁性体特性解析装置１の解析結果を表示装置２４により表示した画面例である。図示するように、磁性体特性解析装置１の解析結果は、例えば外部磁界と平均磁化の関係を示すＢＨカーブを示す。

磁場の方程式の時間ステップｎが磁場の非定常計算ステップＮに等しくなると、制御部６０は、本フローチャートの処理を終了する。

なお、Ｓ１３２における外部磁界→Ｈ_externalの変更程度は、例えば磁場の方程式の時間ステップｎが磁場の非定常計算ステップＮに等しくなる時点で一巡するように、予め設定されている。

本実施例の磁性体特性解析装置１は、上記説明したループ（Ａ）、（Ｂ）の処理によって、平均磁化＜→Ｍ＞の収束を早め、処理速度を向上させることができる。第１演算処理部４０がベクトルポテンシャルＡ_j ⁿ⁺¹に基づき算出する磁界の成分Ｈ₁ ⁿ⁺¹は、磁化分布、励磁電流、渦電流等に起因する成分であり、平均磁化＜→Ｍ＞が変化する時間と比較すると、緩やかに変化する。ベクトルポテンシャルＡ_j ⁿ⁺¹が変化する時間スケールは、例えばマイクロ秒オーダーであるのに対し、平均磁化＜→Ｍ＞が変化する時間スケールは、例えばナノ秒オーダーである。このため、平均磁化＜→Ｍ＞が変化する時間スケールは、ベクトルポテンシャルＡ_j ⁿ⁺¹が変化する時間スケールの１／１０００以下となり得る。

この結果、磁界の成分Ｈ₁ ⁿ⁺¹は、磁界の成分Ｈ₂ ⁿ⁺¹に比して緩やかに時間変化し、磁界の成分Ｈ₂ ⁿ⁺¹の変化に対して静止しているものとして扱うことができる。これに対し、磁界の成分Ｈ₂ ⁿ⁺¹は、ＬＬＧ方程式の時間発展に関するループ（Ａ）の中で平均磁化＜→Ｍ＞と共に変化する磁界であり、解析結果となり得るＢＨカーブの形に大きく影響する。

更に、有限要素法による磁場解析では、分割された各要素に関する連立方程式を解く必要があるため、Ｓ１１２の磁場解析は、Ｓ１１８のＬＬＧ方程式の時間積分に比して、多くの処理時間を要するものである。

従って、ループ（Ａ）において、磁界の成分Ｈ₁ ⁿ⁺¹を固定しつつ、平均磁化＜→Ｍ＞に比例する磁界の成分Ｈ₂ ⁿ⁺¹を変化させることにより、ループ（Ｂ）における平均磁化＜→Ｍ＞の収束を早めることができる。

また、Ｓ１１２の処理を一回実行した後、ループ（Ａ）の処理を複数回実行することにより、例えばこれらを一回ずつ交互に実行する場合に比して、多くの処理時間を要する有限要素法による磁場解析の実行頻度を低下させ、高速且つ高精度に磁性体の特性を解析することができる。

［参考手法との比較］
ここで、有限要素法による磁場解析と、ＬＬＧ方程式の時間積分とを一回ずつ交互に実行する手法（以下、参考手法と称する）との比較について説明する。なお、本実施例の処理結果は、与えられた外部磁界→Ｈ_externalに対して平均磁化＜→Ｍ＞が収束した際の平均磁化＜→Ｍ＞に関して、参考手法による処理結果と一致する筈である。

図１１は、参考手法において実行される処理を概念的に示す図である。また、図１２は、参考手法を実行するための処理の流れを示すフローチャートである。

参考手法では、図１１に示すように、有限要素法による磁場解析の結果である磁界→Ｈⁿ⁺¹は、成分→Ｈ₁ ⁿ⁺¹、→Ｈ₂ ⁿ⁺¹に分解されず、磁界→Ｈⁿ⁺¹全体として有効磁界の算出に用いられる。ＬＬＧ方程式の算出・積分、及び平均磁化の算出に関しては、図８に示す本実施例の手法と同様である。

また、図１２に示す参考の処理の流れは、図９に示す本実施例の処理の流れと比較すると、Ｓ１１４、Ｓ１２２の処理が省略され、ループ（Ａ）が無くなっている。このため、前述した理由により、平均磁化＜→Ｍ＞が収束するまでのループ（Ｂ）の実行回数が、本実施例の処理に比して多くなってしまう。

本出願の出願人は、図８、９に示す本実施例の手法と、図１１、１２に示す参考手法の双方により磁性体の特性解析シミュレーションを行った。係るシミュレーションは、ＬＬＧループ回数Ｍを２００として実行された。図１３は、本実施例の手法による解析結果と、参考手法による解析結果を比較した図である。図示するように、本実施例の手法によれば、参考手法と同様の結果を得られることが確かめられた。

更に、上記磁性体の特性解析シミュレーションにより、本実施例の手法を用いることにより、参考手法と比較して処理時間が短縮されることが確かめられた。図１４、１５は、係る結果を示す図である。

図１４は、本実施例の手法と参考手法におけるループ（Ｂ）の実行回数を比較した図である。図示するように、本実施例の手法では、参考手法に比して平均磁化＜→Ｍ＞の収束が早まるため、ループ（Ｂ）の実行回数は約１／４０となった。なお、ＬＬＧループ回数Ｍを２００に設定しているため、図９におけるＳ１１６〜Ｓ１２０の実行回数は、本実施例の手法の方が５倍（２００／４０）増加している。

図１５は、本実施例の手法と参考手法の処理時間を比較した図である。図示するように本実施例の手法では、参考手法に比して処理時間が１／７程度に短縮された。

この処理時間の短縮効果は、以下のような理由によると考えられる。有限要素法による磁場解析に要する処理時間を１とすると、図９のＳ１１６〜Ｓ１２０の処理に要する処理時間は、０．０２〜０．０２５程度である。従って、参考手法によりループ（Ｂ）を一回実行する時間を（１＋０．０２＝１．０２）とすると、本実施例の手法においてループ（Ｂ）を一回実行する時間は（１＋０．０２×２００＝５）となる。更に、本実施例の手法はループ（Ｂ）の実行回数が参考手法の約１／４０であるため、本実施例の手法におけるループ（Ｃ）の一回の実行時間は、（５／４０＝０．１２５）となる。これを参考手法におけるループ（Ｃ）の一回の実行時間（１．０２）で除算すると、（０．１２５／１．０２＝０．１２２（約１／８））となる。ループ（Ｃ）の実行回数は、外部磁界→Ｈ_externalの変更幅等によって決定される任意事項であるため、ループ（Ｃ）の一回の実行時間比を、処理時間比とみなすことができる。

以上より、本実施例の手法によれば、参考手法と同様の解析結果が得られると共に、処理時間を短縮することができることが分かる。

［解析結果例］
また、本出願の出願人は、異なる形状を有する複数の磁性体モデルについて、本実施例の手法により特性解析シミュレーションを行った。図１６は、解析対象となった複数の磁性体モデル（タイプＡ、タイプＢ、タイプＣ、タイプＤ、タイプＥ、タイプＦ）を示す図である。図示するように、解析対象となった複数の磁性体モデルは、それぞれ厚さが異なる直方体形状の磁性体モデルである。最も厚いタイプＡの磁性体モデルは、形状によるＢＨカーブへの影響が最も大きく、最も薄いタイプＢの磁性体モデルは、厚さ０を想定しているため、ＢＨカーブへの影響が最も小さくなる。図１７は、磁性体特性解析装置１に与えられた各磁性体モデルの寸法を示す図である。なお、本シミュレーションにおいて、外部磁界→Ｈ_externalは、図１６におけるＸ軸方向に与えられた。

そして、図１８は、異なる形状を有する複数の磁性体モデルについて、本実施例の手法により特性解析シミュレーションを行った結果を示す図である。図示するように、磁性体寸法が最も厚いタイプＡのＢＨカーブの傾きが最も小さく、厚さが薄くなるに従って傾きが増加し、厚さ０を想定したタイプＧの結果に近づくことが確かめられた。この結果、本実施例の手法により、ＬＬＧ方程式を用いた磁性体モデルに対して、有限要素法による磁場解析によって得られる磁性体の形状効果を正しく考慮できていることが分かる。

［まとめ］
以上説明した本実施例の磁性体特性解析プログラム、磁性体特性解析装置、及び磁性体特性解析方法によれば、磁界の成分Ｈ₁ ⁿ⁺¹を固定しつつ、磁界の成分Ｈ₂ ⁿ⁺¹を変化させてＬＬＧ方程式の積分を複数回行うため、平均磁化＜→Ｍ＞の収束を早めることができる。また、Ｓ１１２の処理を一回実行した後、ループ（Ａ）の処理を複数回実行するため、上記参考手法に比して、多くの処理時間を要する有限要素法による磁場解析の実行頻度を低下させ、高速且つ高精度に磁性体の特性を解析することができる。

以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。

例えば、上記実施例では、ＬＬＧループ回数Ｍを固定値としたが、ループ（Ｂ）を繰り返し実行する回数（ＬＬＧ方程式の時間ステップｋ）が増加するのに応じて、ＬＬＧループ回数Ｍが徐々に小さくなるように変更してもよい。平均磁化＜→Ｍ＞の収束が近づくと、平均磁化＜→Ｍ＞の変動幅が徐々に小さくなるからである。

以上の説明に関し、さらに以下の項を開示する。
（付記１）
磁性体特性解析プログラムにおいて、
磁性体特性解析装置に、
解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行わせ、
前記有限要素法による磁場解析により算出される磁界を固定値として前記各要素に作用する有効磁界を算出させ、
前記算出した有効磁界を用いてＬＬＧ（Landau Lifshitz Gilbert）方程式を時間積分して前記各要素内で磁化ベクトルを算出させ、
前記算出した磁化ベクトルを平均して前記各要素についての平均磁化を算出させることを特徴とする磁性体特性解析プログラム。
（付記２）
前記磁性体特性解析プログラムは、
前記磁性体特性解析装置に、
前記算出した磁化ベクトルを平均して前記各要素についての平均磁化の算出を、前記平均磁化に依存する磁界の変化を前記有効磁界に再帰的に反映させて複数回行わせることを特徴とする付記１記載の磁性体特性解析プログラム。
（付記３）
前記磁性体特性解析プログラムは、
前記磁性体特性解析装置に、
有限要素法による磁場解析を行う処理と、
前記有効磁界を算出させ、前記磁化ベクトルを算出させ、前記各要素についての平均磁化を算出させる処理と、を交互に繰り返し実行させることを特徴とする付記１又は２記載の磁性体特性解析プログラム。
（付記４）
前記磁性体特性解析プログラムは更に、
前記磁性体特性解析装置に、
設定値である外部磁界を反映させて前記有効磁界を算出させ、
前記複数回行う処理を実行させた後、該複数回行う処理において算出された平均磁化と、前回の前記複数回行う処理において算出された平均磁化との変化を示す指標値が所定値未満であるか否かを判定させ、該指標値が所定値未満である場合に前記外部磁界を変更させることを特徴とする、
付記１ないし３のいずれか１項に記載の磁性体特性解析プログラム。
（付記５）
前記磁性体特性解析プログラムは更に、
前記指標値が所定値未満となった際に前記複数回行う処理により算出された平均磁化に応じた値と、該指標値が所定値未満となった際に設定されている外部磁界と、の組み合わせを解析結果として前記磁性体特性解析装置に出力させることを特徴とする、
付記４記載の磁性体特性解析プログラム。
（付記６）
前記磁性体特性解析プログラムは更に、
前記磁性体特性解析装置に、
前記複数回行う処理における実行回数を、前記外部磁界が不変である期間の経過に応じて徐々に減少させることを特徴とする、
付記３記載の磁性体特性解析プログラム。
（付記７）
解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行う磁場解析部と、
前記有限要素法による磁場解析により算出される磁界を固定値として前記各要素に作用する有効磁界を算出する有効磁界算出部と、
前記算出した有効磁界を用いてＬＬＧ（Landau Lifshitz Gilbert）方程式を時間積分して前記各要素内で磁化ベクトルを算出する磁化ベクトル算出部と、
前記算出した磁化ベクトルを平均して前記各要素についての平均磁化を算出する平均磁化算出部を有することを特徴とする磁性体特性解析装置。
（付記８）
前記平均磁化算出部は、前記平均磁化の算出を、前記平均磁化に依存する磁界の変化を前記有効磁界に再帰的に反映させて複数回行うことを特徴とする付記７記載の磁性体特性解析装置。
（付記９）
前記磁場解析部による処理と、
前記有効磁界算出部による処理、前記磁化ベクトル算出部による処理、及び前記平均磁化算出部による処理と、を交互に繰り返し実行することを特徴とする付記７又は８記載の磁性体特性解析装置。
（付記１０）
前記有効磁界算出部は、設定値である外部磁界を反映させて前記有効磁界を算出し、
前記平均磁化算出部が平均磁化を算出した後、平均磁化の今回値と前回値との変化を示す指標値が所定値未満であるか否かを判定し、該指標値が所定値未満である場合に前記外部磁界を変更することを特徴とする、
付記７ないし９のいずれか１項に記載の磁性体特性解析装置。
（付記１１）
前記指標値が所定値未満となった際に前記平均磁化算出部により算出された平均磁化に応じた値と、該指標値が所定値未満となった際に設定されている外部磁界と、の組み合わせを解析結果として出力することを特徴とする、
付記１０記載の磁性体特性解析装置。
（付記１２）
前記平均磁化算出部における実行回数を、前記外部磁界が不変である期間の経過に応じて徐々に減少させることを特徴とする、
付記９記載の磁性体特性解析装置。
（付記１３）
磁性体特性解析方法において、
磁性体特性解析装置が、
解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行い、
前記有限要素法による磁場解析により算出される磁界を固定値として前記各要素に作用する有効磁界を算出し、
前記算出した有効磁界を用いてＬＬＧ（Landau Lifshitz Gilbert）方程式を時間積分して前記各要素内で磁化ベクトルを算出し、
前記算出した磁化ベクトルを平均して前記各要素についての平均磁化を算出することを特徴とする磁性体特性解析方法。
（付記１４）
前記磁性体特性解析方法は、
前記磁性体特性解析装置が、
前記算出した磁化ベクトルを平均して前記各要素についての平均磁化の算出を、前記平均磁化に依存する磁界の変化を前記有効磁界に再帰的に反映させて複数回行うことを特徴とする付記１３記載の磁性体特性解析方法。
（付記１５）
前記磁性体特性解析方法は、
前記磁性体特性解析装置が、
有限要素法による磁場解析を行う処理と、
前記有効磁界を算出し、前記磁化ベクトルを算出し、前記各要素についての平均磁化を算出する処理と、を交互に繰り返し実行することを特徴とする付記１３又は１４記載の磁性体特性解析方法。
（付記１６）
前記磁性体特性解析方法は更に、
前記磁性体特性解析装置が、
設定値である外部磁界を反映させて前記有効磁界を算出し、
前記複数回行う処理を実行させた後、該複数回行う処理において算出された平均磁化と、前回の前記複数回行う処理において算出された平均磁化との変化を示す指標値が所定値未満であるか否かを判定させ、該指標値が所定値未満である場合に前記外部磁界を変更させることを特徴とする、
付記１３ないし１５のいずれか１項に記載の磁性体特性解析方法。
（付記１７）
前記磁性体特性解析方法は更に、
前記指標値が所定値未満となった際に前記複数回行う処理により算出された平均磁化に応じた値と、該指標値が所定値未満となった際に設定されている外部磁界と、の組み合わせを解析結果として前記磁性体特性解析装置が出力することを特徴とする、
付記１６記載の磁性体特性解析方法。
（付記１８）
前記磁性体特性解析方法は更に、
前記磁性体特性解析装置が、
前記複数回行う処理における実行回数を、前記外部磁界が不変である期間の経過に応じて徐々に減少させることを特徴とする、
付記１５記載の磁性体特性解析方法。
（付記１９）
磁性体特性解析プログラムにおいて、
磁性体特性解析装置に、
解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行う処理と、
前記各要素について、前記有限要素法による磁場解析において算出される磁界を用いて得られるＬＬＧ方程式を時間積分して前記平均磁化を算出する処理を行う処理とを、
前記平均磁化を算出する処理を前記有限要素法による磁場解析を行う処理よりも多く実行させることを特徴とする磁性体特性解析プログラム。

１ 磁性体特性解析装置
１０ ＣＰＵ
１２ ドライブ装置
１４ 記憶媒体
１６ 補助記憶装置
１８ メモリ装置
２０ インターフェース装置
２２ 入力装置
２４ 表示装置
３０ 設定入力受付部
４０ 第１演算処理部
５０ 第２演算処理部
５１ 有効磁界算出部
５２ 磁化ベクトル算出部
５３ 平均磁化算出部
６０ 制御部

Claims (7)

1. 磁性体特性解析プログラムにおいて、
   磁性体特性解析装置に、
   解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行わせ、
   前記有限要素法による磁場解析により算出される磁界のうち磁束密度で決まる第１成分を固定値として前記各要素に作用する有効磁界を算出させ、
   前記算出した有効磁界を用いてＬＬＧ（Landau Lifshitz Gilbert）方程式を時間積分して前記各要素内で磁化ベクトルを算出させ、
   前記算出した磁化ベクトルを平均して前記各要素についての平均磁化を算出させ、
   前記固定値は、前記平均磁化に依存する磁界の変化を前記有効磁界に再帰的に反映させて前記各要素についての平均磁化を算出する再帰処理おいて固定する値であることを特徴とする磁性体特性解析プログラム。
2. 前記磁性体特性解析プログラムは、
   前記磁性体特性解析装置に、
   前記有限要素法による磁場解析を行う処理と、
   前記有効磁界を算出させ、前記磁化ベクトルを算出させ、前記各要素についての平均磁化を算出させる前記再帰処理と、を交互に繰り返し実行させることを特徴とする請求項１記載の磁性体特性解析プログラム。
3. 前記磁性体特性解析プログラムは更に、
   前記磁性体特性解析装置に、
   設定値である外部磁界を反映させて前記有効磁界を算出させ、
   前記再帰処理を実行させた後、該再帰処理において算出された平均磁化と、前回の再帰処理において算出された平均磁化との変化を示す指標値が所定値未満であるか否かを判定させ、該指標値が所定値未満である場合に前記外部磁界を変更させることを特徴とする、
   請求項１または２に記載の磁性体特性解析プログラム。
4. 前記磁性体特性解析プログラムは更に、
   前記指標値が所定値未満となった際に前記再帰処理により算出された平均磁化に応じた値と、該指標値が所定値未満となった際に設定されている外部磁界と、の組み合わせを解析結果として前記磁性体特性解析装置に出力させることを特徴とする、
   請求項３記載の磁性体特性解析プログラム。
5. 前記磁性体特性解析プログラムは更に、
   前記磁性体特性解析装置に、
   設定値である外部磁界を反映させて前記有効磁界を算出させ、
   前記交互に繰り返し実行する実行回数を、前記外部磁界が不変である期間の経過に応じて徐々に減少させることを特徴とする、
   請求項２記載の磁性体特性解析プログラム。
6. 解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行う磁場解析部と、
   前記有限要素法による磁場解析により算出される磁界のうち磁束密度で決まる第１成分を固定値として前記各要素に作用する有効磁界を算出する有効磁界算出部と、
   前記算出した有効磁界を用いてＬＬＧ（Landau Lifshitz Gilbert）方程式を時間積分して前記各要素内で磁化ベクトルを算出する磁化ベクトル算出部と、
   前記算出した磁化ベクトルを平均して前記各要素についての平均磁化を算出する平均磁化算出部と、
   を有し、
   前記平均磁化算出部は、前記平均磁化に依存する磁界の変化を前記有効磁界に再帰的に反映させて前記各要素についての平均磁化を算出する再帰処理おいて、前記固定値を用いることを特徴とする
   磁性体特性解析装置。
7. 磁性体特性解析方法において、
   磁性体特性解析装置が、
   解析対象が分割された各要素について与えられる平均磁化を用いて、有限要素法による磁場解析を行い、
   前記有限要素法による磁場解析により算出される磁界のうち磁束密度で決まる第１成分を固定値として前記各要素に作用する有効磁界を算出し、
   前記算出した有効磁界を用いてＬＬＧ（Landau Lifshitz Gilbert）方程式を時間積分して前記各要素内で磁化ベクトルを算出し、
   前記算出した磁化ベクトルを平均して前記各要素についての平均磁化を算出し、
   前記平均磁化の算出は、前記平均磁化に依存する磁界の変化を前記有効磁界に再帰的に反映させて前記各要素についての平均磁化を算出する再帰処理おいて、前記固定値を用いることを特徴とする磁性体特性解析方法。