JP6015761B2 - 磁性体のシミュレーションプログラム、シミュレーション装置及びシミュレーション方法 - Google Patents

Description

本発明は、磁性体のシミュレーションプログラム、シミュレーション装置及びシミュレーション方法に関する。

モータや発電機といった磁性材料を用いた電気機器のシミュレーションは、コンピュータの性能向上と電磁界解析手法の進歩により、さまざまな場面で広く行なわれるようになってきている。電磁界解析の手法としては、差分法や有限要素法が一般的に用いられている。近年、ＣＯ₂削減や地球温暖化防止への取組みとして、電気機器の効率は非常に重要視されており、シミュレーションへの期待はさらに高まっている。

モータや変圧器といった低周波の交流電流が入力されるデバイスに最も用いられている材料には、電磁鋼板がある。電磁鋼板とは、電気エネルギーと磁気エネルギーとの変換効率を高くした板状の鋼である。電磁鋼板のヒステリシス性を表す際に問題になるのが磁気異方性である。電磁鋼板には２種類あり、特定の方向に対し磁気特性がよい方向性電磁鋼板とある程度どの方向に対しても磁気特性のよい無方向性電磁鋼板である。特に方向性電磁鋼板では磁界印加方向により磁化曲線が大きく異なる。

図１３は、方向性電磁鋼板の磁化曲線の一例を示すグラフである。このようなヒステリシス曲線を数値計算で取り扱う方法として、ストップモデルがある。この方法は数学的な観点から任意のヒステリシス特性を表現する方法として利用される（下記非特許文献１を参照。）。

また、ミクロな視点でモデル化する手法として、マイクロマグネティックスに基づいた技術が開示されている（下記特許文献１を参照。）。特許文献１の従来技術は磁壁移動を簡略化した形で取り扱う。

国際公開第２０１１／１１４４９２号

Ｔｅｔｓｕｊｉ Ｍａｔｓｕｏ， Ｍａｓａａｋｉ Ｓｈｉｍａｓａｋｉ， "Ｒｅｐｒｅｓｅｎｔａｔｉｏｎ Ｔｈｅｏｒｅｍｓ ｆｏｒ Ｓｔｏｐ ａｎｄ Ｐｌａｙ Ｍｏｄｅｌｓ Ｗｉｔｈ Ｉｎｐｕｔ−Ｄｅｐｅｎｄｅｎｔ Ｓｈａｐｅ Ｆｕｎｃｔｉｏｎｓ"， ＩＥＥＥ ＴＲＡＮＳＡＣＴＩＯＮＳ ＯＮ ＭＡＧＮＥＴＩＣＳ， ＶＯＬ． ４１， ＮＯ． ５， ＭＡＹ ２００５， Ｐａｇｅｓ １５４８−１５５１．

しかしながら、上述した非特許文献１の従来技術では、ヒステリシス曲線を数値計算する際に、ヒステロン分布を求める必要があり、ヒステロン分布を同定するために多くの測定をしなければならないという問題がある。また、測定される磁性体の磁気特性は材料のもつミクロな結晶構造などに起因するが、結晶構造が変わった場合の特性を計算することはできない。

また、上述した特許文献１の従来技術では、方向性電磁鋼板や無方向性電磁鋼板でみられる磁気特性の方向依存性を再現することはできない。したがって、非特許文献１の従来技術に、特許文献１の従来技術を適用しても、結晶構造がどのように変化すると特性がどのくらい変化するか把握するのが困難であるという問題がある。

１つの側面では、磁性体内部の方向依存性を高精度に再現することができる磁性体のシミュレーションプログラム、シミュレーション装置、およびシミュレーション方法を提供することを目的とする。

本発明の一側面によれば、磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成し、分割された領域の磁化ごとの磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択し、前記領域の磁化と、選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、前記領域内のピニングによるエネルギー障壁の高さに応じた磁化反転の反転角度を特定し、特定された反転角度により前記磁化を反転させる磁性体のシミュレーションプログラム、シミュレーション装置及びシミュレーション方法が提案される。

本発明の一側面によれば、磁性体内部の方向依存性を高精度に再現することができるという効果を奏する。

図１は、実施の形態１にかかるシミュレーションの一例を示す説明図である。 図２は、シミュレーションを実行するシミュレーション装置２００のハードウェア構成例を示すブロック図である。 図３は、図２に示したシミュレーション装置２００を用いたシステム構成例を示す説明図である。 図４は、磁性体ＤＢの記憶内容の一例を示す説明図である。 図５は、本実施の形態にかかるシミュレーション装置２００の機能的構成例を示すブロック図である。 図６は、フェライトに対する直流でのヒステリシス曲線、ならびに、１.０ＭＨｚでのヒステリシス曲線の実測データを示すグラフである。 図７は、本実施の形態でのシミュレーション結果を示すグラフである。 図８は、図５に示した反転処理部５００の機能的構成例を示すブロック図である。 図９は、実施の形態にかかるシミュレーション装置２００によるシミュレーション処理手順例を示すフローチャートである。 図１０は、図９に示した磁性体内の要素Ｃｇでのシミュレーション処理（ステップＳ９０３）の詳細な処理手順を示すフローチャート（その１）である。 図１１は、図９に示した磁性体内の要素Ｃｇでのシミュレーション処理（ステップＳ９０３）の詳細な処理手順を示すフローチャート（その２）である。 図１２は、図１０に示した反転処理（ステップＳ１００４）の詳細な処理手順例を示すフローチャートである。 図１３は、方向性電磁鋼板の磁化曲線の一例を示すグラフである。

以下に添付図面を参照して、この発明にかかる磁性体のシミュレーションプログラム、シミュレーション装置、およびシミュレーション方法の実施の形態を詳細に説明する。なお、本明細書では、図面および数式において太字で表現されているベクトルを、『［ ］』を用いて表記する。たとえば、ベクトルＸは、［Ｘ］とする。また、図面および数式においてベクトルの上に表記されているバーは、平均を表している。明細書中では、バーが表記されているベクトルを、『［ ａｖ］』を用いて表記する。たとえば、バーがあるベクトルＸは、［Ｘａｖ］とする。

（実施の形態１）
＜シミュレーションの一例＞
図１は、実施の形態１にかかるシミュレーションの一例を示す説明図である。磁性体モデル１００は、磁性体の一例である方向性電磁鋼板をモデル化した電子データである。方向性電磁鋼板は、たとえば、周波数ｆが１.０ＭＨｚ以上の交流磁界となる磁性体である。

（Ａ）は、方向性電磁鋼板の磁性体モデル１００を示す。（Ａ）において、ＲＤは加工時の圧延方向であり、ＴＤは圧延方向に直交する方向である。方向性電磁鋼板の磁性体は、圧延方向ＲＤと直交方向ＴＤとで特性が異なる。圧延方向ＲＤの透磁率は直交方向ＴＤに比べて高く、また、圧延方向ＲＤのヒステリシス損失は直交方向ＴＤに比べて少ない。本実施の形態では、圧延方向ＲＤと直交方向ＴＤとの特性の違いが再現できるようにシミュレーションが実行される。磁性体モデル１００は、Ｎ個（Ｎ≧２）の要素Ｃ１〜ＣＮに個々に割り当てられる。

（Ｂ）は、磁性体モデル１００のある要素Ｃｇの拡大図である。（Ａ）に示したように、磁性体モデル１００は３次元で表現されるが、（Ｂ）では便宜上２次元で表現する。要素Ｃｇ（１≦ｇ≦Ｎ）は、ｎ個（ｎ≧２）のメッシュｃ１〜ｃｎをもつ。ｎ個（ｎ≧２）のメッシュｃ１〜ｃｎは、たとえば、マイクロマグネティクスの磁化分布を表す。マイクロマグネティクスにおける磁性体内の磁気的エネルギーは、下記式（１）〜（７）の磁気異方性エネルギーＥ_ani、静磁エネルギーＥ_mag、交換相互エネルギーＥ_excおよびゼーマンエネルギーＥ_extで表現される。

ここで、Ｋ：磁気異方性係数、［ｋ］：磁化容易軸方向の単位ベクトル、［ｍ］：磁化方向の単位ベクトル、［Ｍ］：メッシュｃｉの磁化、Ｍ：磁化［Ｍ］の大きさ、Ｄ_ij：メッシュｃｉの幾何学的形状から決定される反磁界係数、［ｒ］：メッシュｃｉの位置ベクトル、［ｒ'］：メッシュｃｊの位置ベクトル、Ａ^*：スティフィネス・コンスタント、ａ：メッシュ間隔、［Ｈ_app］：外部印加磁界、ｎ：メッシュ数である。なお、式（５）のΣ［Ｍ_j］は、［Ｍ_i］のメッシュに隣接するメッシュの磁化の総和である。

λ₁₀₀とλ₁₁₁は磁歪定数である。磁歪定数は結晶軸によって異なる値をもち、λ₁₀₀は、（１００）方向の磁歪定数であり、λ₁₁₁は（１１１）方向の磁歪定数である。また、αおよびγは以下の式（８），（９）で与えられる。

ここで、［ａ_1,i］、［ａ_2,i］および［ａ_3,i］は、ｉ番目の磁化に対応する結晶容易軸であり、立方晶の場合は結晶方向を示すａ軸、ｂ軸、ｃ軸方向をもち、長さが１となるベクトルである。［σ］は応力の印加方向の単位ベクトルである。容易軸は方向性電磁鋼板では圧延方向ＲＤをｘ軸とする場合、［ａ_1,i］＝（１，０，０）の値をもつ。電磁鋼板の圧延方向ＲＤに対し直交方向ＴＤをｙ軸とすると、［ａ_2,i］＝（０，１／√２，１／√２）である。［ａ_3,i］は、［ａ_1,i］および［ａ_2,i］に直交する向きに設定される。

図１において、各メッシュｃｉ（１≦ｉ≦ｎ）内の矢印は、磁化［Ｍ_i］を示している。磁化［Ｍ_i］は微小時間Δτごとに変化する。各メッシュｃｉの磁化［Ｍ_i］の合計をメッシュ数ｎで割ることで、平均磁化［Ｍａｖ］が求められる。なお、メッシュｃｉにおける磁化方向の単位ベクトル［ｍ_i］は、［ｍ_i］＝［Ｍ_i］／Ｍ_iで表される。

また、磁性体の要素Ｃｇ内の磁気的な全エネルギーＥ_totは、下記式（１０）に示すように磁気異方性エネルギーＥ_ani、静磁エネルギーＥ_mag、交換エネルギーＥ_exc、ゼーマンエネルギーＥ_ext、および応力による磁気弾性エネルギーＥσ（［ｍ_i］）の和で表される。

これらのエネルギーは計算対象により異なり、上記エネルギーの組み合わせであり１つ以上の項をもつものとして計算される。ここで磁化状態を計算する際に磁化回転と磁化反転の２種類の過程を計算する。１つ目はランダウ−リフシッツ−ギルバート（Ｌａｎｄａｕ−Ｌｉｆｓｈｉｔｚ−Ｇｉｌｂｅｒｔ）方程式（以下、「ＬＬＧ方程式」）に伴う方法である。そこでは有効磁界に従って磁化の回転運動が計算される。

（Ｃ）は、磁化の回転運動の概念を示す。磁化回転運動ではそれぞれの磁化はＬＬＧ方程式に従って変化する。その場合、磁化は、エネルギーの小さい方向に回転するため、ｕ１に示すように、連続的にエネルギーの小さい磁化方向へと変化する。一方で磁化反転においては、ｕ２に示すように、磁化は、エネルギー障壁である磁壁のピニングＰを乗り越えて、１８０度または９０度に対応するエネルギーの位置へ不連続に変化する。本実施の形態では、シミュレーションでメッシュｃｉごとに（Ｃ）を再現することにより、（Ｄ）に示すような磁化回転または磁化反転後の状態を得る。

このように、（Ａ）のように巨視的にみれば磁性体の磁化の方向は圧延方向ＲＤとなるが、（Ｄ）のように微視的にみれば要素Ｃｇ内のメッシュｃｉの磁化［Ｍ_i］の方向は、一様ではない。これは、磁性体内部に存在する不純物などの存在を考慮してシミュレーションが実行されたからである。このようなシミュレーションをシミュレーション装置が実行することにより、磁性体内部において、磁化の方向が特定の結晶軸の方向に依存して反転するという方向依存性を、高精度に再現することができる。

＜磁化回転運動＞
ここで、（Ｃ）のｕ１に示す磁化回転運動の計算について説明する。磁化回転運動はエネルギーにより計算される有効磁界から計算される。有効磁界は式（１１）のように算出される。

式（１１）はメッシュｃｉ内の磁化［Ｍ_i］に対する磁気的エネルギーによる磁界である。有効磁界［Ｈ_i］は、次の時刻での磁化［Ｍ_i］の時間変化を求めるのに用いられる。式（１１）を磁気異方性Ｈ_k＝２Ｋ／Ｍによって規格化して、下記式（１２）を得る。

ここで、ｈ_m，ｈ_eは、下記式（１３），（１４）に示すように、それぞれ磁気異方性Ｈ_kで規格化された静磁界係数および交換相互作用係数である。また、ｈ_aは磁気異方性Ｈ_kで規格化された外部印加磁界Ｈ_app／Ｈ_kである。

また、磁性体の磁化の運動は、下記式（１５）に示すＬＬＧ方程式によって決定される。

ここで、αはダンピングコンスタントである。ダンピングコンスタントαは、減衰していく過程での速さを表す磁性体固有の定数である。式（１５）の右辺の第１項は歳差運動項であり、第２項はダンピング項である。フェライトの場合、たとえば、α＝０．１を用いた。τはＬＬＧ方程式の計算を実行する時間変数である。

ここでＬＬＧ方程式の分布計算方法を、下記式（１６）を用いて示す。式（１６）を計算するにあたり、ΔτはＬＬＧ方程式の計算時間ステップ幅である。ここではステップごとにステップ幅は一様とする。

［ｍ_i］^oldは、現在の時間ステップτでの［ｍ_i］であり、［ｍ_i］^newは、次の時間ステップ（τ＋Δτ）での［ｍ_i］である。

また、時間ステップ毎に以下の式（１７），（１８）が評価される。

ここで、［Ｍａｖ］^jは、要素Ｃｇ内の磁化［Ｍ_i］がいずれも収束したとき（時刻ｔ_j）での要素Ｃｇ内の平均磁化である。［Ｍａｖ］^j-1は、［Ｍａｖ］^jの直近で要素Ｃｇ内の磁化［Ｍ_i］がいずれも収束したときの平均磁化である。また、［Ｍａｖ］^j-2は、［Ｍａｖ］^j-1の直近で要素Ｃｇ内の磁化［Ｍ_i］がいずれも収束したときの平均磁化である。

上述した磁化の収束判定は、たとえば、下記式（１９）で実行される。

シミュレーション装置は、すべての［ｍ_i］に対して微小量ε_mより小さくなった場合に収束したと判定する。ここではε_mは１．ｅ−６を用いた。

＜磁化反転運動＞
つぎに、磁化反転運動について説明する。磁化反転過程は、エネルギーが小さくなる容易軸へと磁化の向きを急激に変えることに対応する。これは磁壁移動により微小領域における磁化が全体的に回転するのではなく、磁壁移動によりエネルギーの小さい磁区が増加することと対応する。しかし、一般的には磁壁移動は必ずしもスムーズに動かず、不純物等によるピニングと呼ばれる磁壁移動の抑制が生じる。磁壁のピニングがあるため、磁化がさらにエネルギーの低い状態へ遷移できない。磁化反転では１８０度および９０度に対応する不連続な磁化変化を考え、その際にエネルギーの大小関係から反転を計算する。１８０度の磁壁移動に相当する１８０度反転過程は、下記式（２０）で計算される。

ここで、ΔＥ_i ¹⁸⁰は、ピニングによるエネルギー障壁の高さに対応するパラメータである。電磁鋼板では主に立方晶による異方性をもつ。つまり、結晶軸であるａ軸，ｂ軸，ｃ軸方向に磁気異方性エネルギーが小さくなる。そのため、１８０度の反転だけでなく、９０度の反転も存在する。９０度反転過程は、下記式（２１）で計算される。

式（２１）において、［ａ_j］は容易軸ベクトルである。また、式（２１）のＥ_tot（［ａ_j］）は、［ｍ_i］の代わりに［ａ_j］をＥ_tot（）代入した場合の（１０）式を計算することにより計算される。磁化［Ｍ_i］の単位ベクトル［ｍ_i］は、結晶軸であるａ軸，ｂ軸，ｃ軸のいずれかの軸に平行になりやすいという性質がある。したがって、反転前の磁化［Ｍ_i］の単位ベクトル［ｍ_i］は、９０度反転によりエネルギーの低い結晶軸に反転させられる。たとえば、反転前の磁化［Ｍ_i］の単位ベクトル［ｍ_i］の方向がａ軸方向である場合、反転後の磁化［Ｍ_i］の単位ベクトル［ｍ_i］の方向は、ａ軸方向よりも磁気的エネルギーが低いｂ軸またはｃ軸方向のいずれかとなる。具体的には、反転後の磁化［Ｍ_i］の単位ベクトル［ｍ_i］の方向は、ｂ軸またはｃ軸方向のうち、より磁気的エネルギーが低い方向となる。

この性質を再現するため、本実施の形態では、９０度反転において式（２１）の不等式を満たした場合、反転前の磁化［Ｍ_i］の単位ベクトル［ｍ_i］は、反転前の磁化［Ｍ_i］の単位ベクトル［ｍ_i］よりも磁気的エネルギーの低い容易軸ベクトル［ａ_j］に更新される。反転前の磁化［Ｍ_i］の単位ベクトル［ｍ_i］と容易軸ベクトル［ａ_j］とは直交することになるため、簡単な計算により９０度反転が再現されることになる。

図１の（Ｃ）のｕ２に示したように、磁化反転では、式（２０）の１８０度反転と式（２１）の９０度反転の両方が起こり得る。この場合はエネルギー差が最も大きい状態へ遷移すると考えられるため、シミュレーション装置は、エネルギーが最も小さくなる状態に磁化を反転させることになる。

＜静磁界計算＞
このあと、いずれの磁化［Ｍ_i］も式（１９）の収束判定を満たした場合、要素Ｃｇでは、平均磁化［Ｍａｖ］を用いて静磁界計算が実行される。静磁界計算では、下記式（２２）を用いて、静磁界Ｈｓ_gが求められる。νは、透磁率μの逆数である。

有限要素法における静磁界計算は、電磁気学の基礎方程式であるマックスウェル方程式から導かれる上記式（１６）を用いて計算される。ここで、Ａは磁気ベクトルポテンシャル、Ｊ₀は計算の対象となる磁性体内を流れる電流である。

本実施の形態では、磁性体は高抵抗であるため、磁性体内に流れる電流はＪ₀＝０となる。また、外部電流が存在する場合は、外部電流の影響を考慮する必要があるため、Ｊ₀に外部電流の値を与える。そして、式（２２）の磁気ベクトルポテンシャルＡに要素Ｃｇの位置座標を与えることで、磁気ベクトルポテンシャルＡが求まる。磁気ベクトルポテンシャルＡは、磁束密度［Ｂ］＝μ［Ｈ］＋［Ｍ］とした場合、［Ｂ］＝ｒｏｔ（Ａ）で定義される。したがって、磁気ベクトルポテンシャルＡが求まれば磁束密度［Ｂ］も求まる。磁束密度［Ｂ］が求まれば、磁化［Ｍ］を与えることで静磁界［Ｈ］が求まる。

本例では、Ｍに［Ｍａｖ］を与えることで、静磁界計算により、静磁界［Ｈｓ_g］が求められる。静磁界［Ｈｓ_g］が磁化収束する都度、そのときの平均磁化［Ｍａｖ］と静磁界［Ｈｓ_g］との組が保存される。平均磁化［Ｍａｖ］と静磁界［Ｈｓ_g］は複数組得られるため、グラフにプロットすることで、ヒステリシス曲線ＨＬｇが作成される。そして、ヒステリシス曲線ＨＬｇ内の面積Ｓｇが要素Ｃｇでのヒステリシス損失となる。

このように、高抵抗磁性材料の磁性体であっても、共鳴現象による磁化変化を高精度に再現することができるため、高抵抗磁性材料の磁性体の共鳴現象によるヒステリシス損失を高精度に求めることができる。

＜コンピュータのハードウェア構成例＞
図２は、シミュレーションを実行するシミュレーション装置２００のハードウェア構成例を示すブロック図である。図２において、シミュレーション装置２００は、プロセッサ２０１、記憶装置２０２、入力装置２０３、出力装置２０４、および通信装置２０５が、バス２０６に接続されて構成されるコンピュータである。

プロセッサ２０１は、シミュレーション装置２００の全体の制御を司る。また、プロセッサ２０１は、記憶装置２０２に記憶されている各種プログラム（ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）や本実施の形態のプログラム）を実行することで、記憶装置２０２内のデータを読み出したり、実行結果となるデータを記憶装置２０２に書き込んだりする。

記憶装置２０２は、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、フラッシュメモリ、磁気ディスクドライブなどで構成され、プロセッサ２０１のワークエリアになったり、各種プログラム（ＯＳや本実施の形態のプログラム）や各種データ（各プログラムの実行により得られたデータを含む）を記憶したりする。

入力装置２０３は、キーボード、マウス、タッチパネルなどユーザの操作により、各種データの入力をおこなうインターフェースである。出力装置２０４は、プロセッサ２０１の指示により、データを出力するインターフェースである。出力装置２０４には、ディスプレイやプリンタが挙げられる。通信装置２０５は、ネットワークを介して外部からデータを受信したり、外部にデータを送信したりするインターフェースである。

＜システム構成例＞
図３は、図２に示したシミュレーション装置２００を用いたシステム構成例を示す説明図である。図３において、ネットワークＮＷは、サーバ３０１，３０２とクライアント３３１〜３３４とが通信可能なネットワークであり、たとえば、ＬＡＮ（Ｌｏｃａｌ Ａｒｅａ Ｎｅｔｗｏｒｋ）、ＷＡＮ（Ｗｉｄｅ Ａｒｅａ Ｎｅｔｗｏｒｋ）、インターネット、携帯電話網などで構成される。

サーバ３０２は、クラウド３２０を構成するサーバ群（サーバ３２１〜３２５）の管理サーバである。クライアント３３１〜３３４のうち、クライアント３３１はノート型パソコン、クライアント３３２はデスクトップ型パソコン、クライアント３３３は携帯電話機（スマートフォン、ＰＨＳ（Ｐｅｒｓｏｎａｌ Ｈａｎｄｙｐｈｏｎｅ Ｓｙｓｔｅｍ）でもよい）、クライアント３３４はタブレット型端末である。図３のサーバ３０１，３０２，３２１〜３２５、クライアント３３１〜３３４は、たとえば、図２に示したシミュレーション装置２００により実現される。なお、クライアント３３１〜３３４は、必ずしもネットワークＮＷに接続する必要はない。

＜磁性体ＤＢ（データベース）の記憶内容例＞
図４は、磁性体ＤＢの記憶内容の一例を示す説明図である。図４において、磁性体ＤＢ４００には、磁性体ごとに、磁性体ＩＤ、ダンピングコンスタントα、摩擦項係数β、慣性項係数γ、導電率σ、異方性磁界Ｈｋ、容易軸タイプＡＥが格納されている。

導電率σ、異方性磁界Ｈｋは、要素Ｃｇでの静磁界計算に用いられる。容易軸タイプＡＥとは、磁性体に適用される磁化容易軸（以下、単に「容易軸」）を示す。容易軸とは、磁化され易い結晶方位である。たとえば、鉄の場合、容易軸は６個存在する。容易軸タイプには、一軸異方性と立方異方性の２種類がある。一軸異方性の場合は、容易軸は、互いに反対方向に向いている２軸であり、立方異方性の場合は、容易軸は、ａ軸、ｂ軸、ｃ軸の３軸と各軸の反対方向の３軸とを合わせて６軸である。鉄の容易軸タイプＡＥは、立方異方性であるため、鉄の容易軸は６軸である。コンピュータにより磁性体ＩＤが指定されると、指定された磁性体ＩＤのレコードの各項の値が磁性体ＤＢ４００から読み出されることになる。磁性体ＤＢ４００は、具体的には、たとえば、図２に示した記憶装置２０２によりその機能を実現する。

＜シミュレーション装置の機能的構成例＞
図５は、本実施の形態にかかるシミュレーション装置２００の機能的構成例を示すブロック図である。シミュレーション装置２００は、磁性体ＤＢ４００と記憶領域５１１を有する。また、シミュレーション装置２００は、反転処理部５００と、平均磁化算出部５０１と、磁界算出部５０２と、磁化算出部５０３と、判断部５０４と、静磁界算出部５０５と、判定部５０６と、格納部５０７と、作成部５０８と、損失算出部５０９と、出力部５１０と、を有する。反転処理部５００〜出力部５１０は、具体的には、たとえば、図２に示した記憶装置２０２に記憶されたシミュレーションプログラムをプロセッサ２０１に実行させることによりその機能を実現する。また、記憶領域５１１は、図２に示した記憶装置２０２によりその機能を実現する。

反転処理部５００は、要素Ｃｇ内のメッシュｃ１〜ｃｎの磁化［Ｍ₁］〜［Ｍ_n］について反転処理を実行する。反転処理部５００の詳細については、図８で説明する。

平均磁化算出部５０１は、時刻τがΔτ変化する都度、要素Ｃｇ内のメッシュｃ１〜ｃｎの磁化［Ｍ₁］〜［Ｍ_n］の平均磁化［Ｍａｖ］を算出する。具体的には、たとえば、下記式（２３）により算出される。

磁界算出部５０２は、磁性体を構成する要素群のいずれかの要素から分割された領域群の領域ごとの磁化が更新されると、領域内の磁気的エネルギーで生成される有効磁界を、領域ごとに算出する。ここで、「領域内の磁気的エネルギーで生成される磁界」とは、メッシュｃｉ内の磁気的エネルギーで生成される磁界であり、上記式（１１）の右辺に相当する。

また、磁界算出部５０２は、更新された領域ごとの磁化の平均磁化の変化を阻止する方向に働く磁化変化速度と、磁化変化速度を一定に保つ方向に磁界を働かせる慣性と、を考慮して有効磁界を算出してもよい。この場合、上記式（１１）に替え、下記式（２４）が用いられる。

「更新された領域ごとの磁化の平均磁化の変化を阻止する方向に働く磁化変化速度」とは、平均磁化［Ｍａｖ］の変化を阻止する方向に働く磁化変化速度であり、式（２４）の右辺第３項の摩擦項に相当する。また、「磁化変化速度を一定に保つ方向に磁界を働かせる慣性」とは、式（２４）の摩擦項で表現された磁化変化速度を一定に保つ方向に磁界を働かせる慣性であり、式（２４）の右辺第２項の慣性項に相当する。なお、式（２４）において、右辺第１項は、上述したメッシュｃｉ内の磁気的エネルギーで生成される磁界であり、上記式（１１）の右辺に相当する。

磁界算出部５０２は、式（１１）または式（２４）により、時刻τがΔτ変化する都度、メッシュｃｉの有効磁界［Ｈ_i］を算出することになる。なお、式（２４）の右辺から慣性項を除いてもよい。有効磁界［Ｈ_i］は、次の時刻での磁化［Ｍ_i］を求めるのに用いられる。

このように、摩擦項を用いた有効磁界［Ｈ_i］を求めることで、メッシュｃｉにおいて、時刻（τ−Δτ）での磁化［Ｍ_i］から時刻τでの磁化［Ｍ_i］への変化量Δ［Ｍ_i］を、時刻（τ−Δτ）での磁化［Ｍ_i］と時刻τでの磁化［Ｍ_i］との摩擦を考慮して求めることができる。したがって、共鳴現象より経時的に変化する磁化［Ｍ_i］での磁化変化を阻止する挙動を再現することができる。

また、摩擦項とともに慣性項も加えた場合、摩擦項および慣性項を用いた有効磁界［Ｈ_i］を求めることで、メッシュｃｉにおいて、時刻（τ−Δτ）での磁化［Ｍ_i］から時刻τでの磁化［Ｍ_i］への変化量Δ［Ｍ_i］を、時刻（τ−Δτ）での磁化［Ｍ_i］と時刻τでの磁化［Ｍ_i］との摩擦および慣性を考慮して求めることができる。したがって、共鳴現象より経時的に変化する磁化［Ｍ_i］での磁化変化を阻止する挙動や磁化変化速度を一定に保つ挙動を再現することができる。

磁化算出部５０３は、領域ごとに算出された有効磁界と、領域ごとの磁化と、に基づいて、領域ごとに磁化変化量を求めることにより、領域ごとの磁化を算出する。磁化算出部５０３では、磁化算出に先立って磁化変化量を算出する。磁化変化量は、時刻がτからτ＋Δτに変化した場合の磁化［Ｍ_i］の単位ベクトル［ｍ_i］の変化量であり、具体的には、上記式（１６）の右辺第２項で表現される。式（１６）の単位ベクトルが［ｍ_i］^oldから［ｍ_i］^newに更新されると、更新後の［ｍ_i］^newに磁化［Ｍ_i］の大きさＭ_iを乗じることで、更新後の磁化［Ｍ_i］が算出される。算出された磁化［Ｍ_i］は、磁界算出部５０２に与えられ、つぎの時刻での有効磁界［Ｈ_i］が算出される。

判断部５０４は、磁化算出部５０３による領域ごとの変化前の磁化と変化後の磁化に基づいて、いずれかの要素での磁化収束を判断する。具体的には、たとえば、判断部５０４は、要素Ｃｇにおいて、上述した更新前の［ｍ_i］^oldと更新後の［ｍ_i］^newを用いて、上記式（１９）により判断する。

静磁界算出部５０５は、上記式（２２）を用いて、要素Ｃｇの静磁界Ｈｓ_gを算出する。具体的には、たとえば、静磁界算出部５０５は、判断部５０４で同一時刻τにおいていずれの磁化［Ｍ_i］も式（１９）の磁化収束条件を満たす場合に、要素Ｃｇの静磁界Ｈｓ_gを算出する。静磁界計算で用いられる透磁率μとしては、たとえば、真空の透磁率が用いられる。

判定部５０６は、静磁界が磁界収束条件を満たすか否かを判定する。具体的には、たとえば、判定部５０６は、今回の静磁界［Ｈｓ_g］^newと前回の静磁界［Ｈｓ_g］^oldとの差分ΔＨｓがしきい値εｈ以内であれば、磁界収束と判定する。

格納部５０７は、式（１９）の磁化収束条件を満たした時刻ｔｊでの平均磁化［Ｍ_i］と判定部５０６による磁界収束条件を満たした静磁界Ｈｓ_gとの組み合わせを記憶領域５１１に格納する。格納された平均磁化［Ｍ_i］と静磁界Ｈｓ_gとの組み合わせ群は、ヒステリシス曲線の元となる。

作成部５０８は、格納された平均磁化［Ｍ_i］と静磁界［Ｈｓ_g］との組み合わせ群を、横軸が磁界、縦軸が磁束密度のグラフにプロットすることで、ヒステリシス曲線を作成する。作成部５０８では、磁束密度［Ｂ］＝μ［Ｈ］＋［Ｍ］に平均磁化［Ｍａｖ］と静磁界［Ｈｓ_g］とを与えることで、時刻ｔｊごとの磁束密度［Ｂ］が得られる。これにより、ヒステリシス曲線を作成することができる。

損失算出部５０９は、作成部５０８で作成されたヒステリシス曲線内の面積を算出することでヒステリシス損失を算出する。式（２４）を適用した場合、式（２４）の摩擦項や慣性項を考慮したヒステリシス損失が得られる。

出力部５１０は、損失算出部５０９によって算出されたヒステリシス損失を出力する。出力部５１０は、ヒステリシス損失をディスプレイに表示してもよく、プリンタにより印刷出力してもよい。また、外部装置に送信したり、記憶装置２０２に格納してもよい。また、出力部５１０は、作成部５０８によって作成されたヒステリシス曲線も出力してもよい。ここで、ヒステリシス曲線の例について説明する。

図６は、フェライトに対する直流でのヒステリシス曲線、ならびに、１.０ＭＨｚでのヒステリシス曲線の実測データを示すグラフである。（ｂ）は（ａ）の拡大グラフである。ヒステリシス曲線のループ内の面積が磁性体内で消費された１サイクル当たりのエネルギー損失となることが知られている。

図７は、本実施の形態でのシミュレーション結果を示すグラフである。式（１１）を用いて計算したヒステリシス曲線が点線７０１で示されている。また、慣性項のみ０として計算したヒステリシス曲線が一点鎖線７０２で示されている。また、式（２４）により慣性項および摩擦項を含めて計算したヒステリシス曲線が実線７０３で示されている。実線７０３のヒステリシス曲線は楕円であり、実測でのヒステリシス曲線が再現されている。点線７０１および一点鎖線７０２のヒステリシス曲線では、楕円で再現できていない。図７では、周波数ｆ＝１.０ＭＨｚ、γ＝０．７５×１０^-11、β＝４．０×１０^-5とした。

＜反転処理部５００の機能的構成例＞
図８は、図５に示した反転処理部５００の機能的構成例を示すブロック図である。反転処理部５００は、生成部８０１と、選択部８０２と、特定部８０３と、反転判断部８０４と、反転部８０５と、を有する。

生成部８０１は、磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成する。具体的には、たとえば、生成部８０１は、磁性体ＤＢ４００を参照して、解析対象となる磁性体の容易軸タイプを特定する。そして、生成部８０１は、容易軸タイプに応じて、容易軸ベクトルを生成する。たとえば、容易軸タイプが、立方異方性である場合、ａ軸〜ｃ軸のそれぞれの正負の方向となる６個の単位ベクトルを容易軸ベクトルとして生成する。また、容易軸タイプが、１軸異方性の場合、たとえば、ａ軸の正負の方向となる２個の単位ベクトルを容易軸ベクトルとして生成する。より具体的には、生成部８０１は、圧延方向ＲＤが（１００）、その直交方向ＴＤが（０１１）となるように容易軸ベクトルを生成する。

選択部８０２は、分割された領域の磁化ごとの磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択する。具体的には、たとえば、選択部８０２は、分割された領域の磁化ごとに、上記式（１０）に示したような磁気的エネルギーＥ_totを算出する。選択部８０２は、式（１０）において、パラメータとして磁化ベクトル［ｍ_i］のかわりに、容易軸ベクトルを与える。これにより、分割された領域の磁化ごとに、磁気的エネルギーＥ_totが算出される。

そして、選択部８０２は、分割された領域の磁化ごとに算出された磁気的エネルギーＥ_totの中から最大ではない磁気的エネルギー、好ましくは最小の磁気的エネルギーを選択する。磁化反転過程では、エネルギーが小さくなる容易軸へと磁化の向きを急激に変える。したがって、選択部８０２は、磁気的エネルギーＥ_totの中から最大ではない磁気的エネルギーを選択することにより、選択された磁気的エネルギーＥ_totの算出元となる容易軸ベクトルへの変化を再現することができる。これにより、磁壁移動により微小領域における磁化の全体的な回転ではなく、磁壁移動によりエネルギーの小さい磁区の増加を再現することができる。

特定部８０３は、領域の磁化と、選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、領域内のピニングによるエネルギー障壁の高さに応じた磁化反転の反転角度を特定する。具体的には、たとえば、特定部８０３は、下記式（２５）を満たすか否かを判断する。特定の容易軸ベクトルとは、選択部８０２によって選択された磁気的エネルギーの算出に用いられた容易軸ベクトルである。ここでは、［ａ_ikmin］とする。

上記式（２５）は、磁化ベクトル［ｍ_i］と特定の容易軸ベクトル［ａ_ikmin］との内積である。すなわち、磁化ベクトル［ｍ_i］と特定の容易軸ベクトル［ａ_ikmin］とのなす角度が１３５度〜２２５度の間であれば、特定部８０３は、磁化ベクトル［ｍ_i］の反転角度を１８０度とし、そうでない場合は、磁化ベクトル［ｍ_i］の反転角度を９０度とする。なお、ここでは、磁化ベクトル［ｍ_i］と特定の容易軸ベクトル［ａ_ikmin］とのなす角度の範囲を１３５度〜２２５度としたが、９０度〜２７０度の範囲内であれば、任意に設定してもよい。

反転部８０５は、特定された反転角度により磁化を反転させる。具体的には、たとえば、反転角度が１８０度である場合、反転部８０５は、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を−［ｍ_i］に更新する。これにより、磁化［Ｍ_i］の単位ベクトル［ｍ_i］は１８０度反転させられた状態となる。また、反転角度が９０度である場合、反転部８０５は、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を特定の容易軸ベクトル［ａ_ikmin］に更新する。これにより、磁化［Ｍ_i］の単位ベクトル［ｍ_i］は９０度反転させられた状態となる。

また、特定部８０３は、磁化と、特定の容易軸ベクトルと、に基づいて、ピニングによるエネルギー障壁の高さに対応するパラメータΔＥを特定する。上述したように、磁化反転過程では、磁壁移動により微小領域における磁化が全体的に回転するのではなく、磁壁移動によりエネルギーの小さい磁区が増加するため、エネルギーが小さくなる容易軸へと磁化の向きが急激に変化する。しかし、一般的には磁壁移動は必ずしもスムーズに動かず、図１の（Ｃ）に示したように、不純物等によるピニングがおきる。

このため、特定部８０３では、１８０度または９０度の磁化反転に対応する不連続な磁化変化について、磁気的エネルギーの大小関係から反転を計算するためのパラメータを特定する。具体的には、１８０度反転の場合は、特定部８０３は、ピニングによるエネルギー障壁の高さに対応するパラメータΔＥをΔＥ_i ¹⁸⁰に設定する。９０度反転の場合は、特定部８０３は、ピニングによるエネルギー障壁の高さに対応するパラメータΔＥをΔＥ_i ⁹⁰に設定する。エネルギー障壁の高さに対応するパラメータは、反転角度ごとに固定の値を用いてもよいが、ここでは、ローレンツ分布を仮定して、乱数により磁化［Ｍ_i］ごとに異なる値とした。例えば、ΔＥ_i ¹⁸⁰は、１０±１５［Ｊ/ｍ³］となり、ΔＥ_i ⁹⁰は、４０±２０［Ｊ/ｍ³］となる。

反転判断部８０４は、磁化に関する領域内の磁気的エネルギーとピニングによるエネルギー障壁の高さに対応するパラメータとに基づいて、反転角度により磁化を反転させるか否かを判断する。具体的には、たとえば、反転判断部８０４は、１８０度反転の場合は、式（２０）に示す１８０度反転条件を満たすか否かを判断する。１８０度反転条件を満たした場合、反転部８０５は、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を−［ｍ_i］に更新する。これにより、磁化［Ｍ_i］の単位ベクトル［ｍ_i］は１８０度反転させられた状態となる。

同様に、反転判断部８０４は、９０度反転の場合は、式（２１）に示す９０度反転条件を満たすか否かを判断する。９０度反転条件を満たした場合、反転部８０５は、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を特定の容易軸ベクトル［ａ_ikmin］に更新する。式（２１）の９０度反転条件を満たした場合は、磁化［Ｍ_i］の単位ベクトル［ｍ_i］と特定の容易軸ベクトル［ａ_ikmin］とのなす角度が９０度であるため、磁化［Ｍ_i］の単位ベクトル［ｍ_i］は９０度反転させられた状態となる。いずれの場合も反転条件を満たさない場合には、反転部８０５は、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を更新しない。

＜シミュレーション処理手順例＞
図９は、実施の形態にかかるシミュレーション装置２００によるシミュレーション処理手順例を示すフローチャートである。まず、シミュレーション装置２００は、要素を特定する変数ｇをｇ＝１とし（ステップＳ９０１）、ｇ＞Ｎであるか否かを判断する（ステップＳ９０２）。Ｎは要素の総数である。ｇ＞Ｎでない場合（ステップＳ９０２：Ｎｏ）、シミュレーション装置２００は、磁性体内の要素Ｃｇでのシミュレーション処理を実行する（ステップＳ９０３）。

このあと、シミュレーション装置２００は、ｇをインクリメントし（ステップＳ９０４）、ステップＳ９０２に戻る。ステップＳ９０２において、ｇ＞Ｎである場合（ステップＳ９０２：Ｙｅｓ）、シミュレーションすべき要素Ｃｇがないため、シミュレーション装置２００は、出力部５１０による出力処理を実行して（ステップＳ９０５）、一連のシミュレーション処理を終了する。なお、図９では、要素Ｃｇごとに順番にシミュレーション処理（ステップＳ９０３）を実行することとしているが、並列に実行することとしてもよい。

図１０は、図９に示した磁性体内の要素Ｃｇでのシミュレーション処理（ステップＳ９０３）の詳細な処理手順を示すフローチャート（その１）である。まず、シミュレーション装置２００は、時刻変数ｊと更新時刻τをｊ＝０、τ＝０とする（ステップＳ１００１）。つぎに、シミュレーション装置２００は、時刻ｔｊでの外部印加磁界［Ｈ_app］を設定する（ステップＳ１００２）。外部印加磁界［Ｈ_app］は周波数ｆと時刻ｔｊで決まる磁界である。周波数ｆはあらかじめ与えられるものとする。外部印加磁界［Ｈ_app］は、上記式（６）で用いられる。

そして、シミュレーション装置２００は、磁化［Ｍ_i］の初期値を記憶装置２０２から読み込む（ステップＳ１００３）。磁化［Ｍ_i］の初期値は、たとえば、あらかじめ記憶装置２０２に記憶させておき、シミュレーション開始時に記憶装置２０２から読み込むものとする。

このあと、シミュレーション装置２００は、反転処理部５００により、磁化［Ｍ_i］の反転処理を実行する（ステップＳ１００４）。反転処理（ステップＳ１００４）の詳細については、図１２で説明する。

つぎに、シミュレーション装置２００は、平均磁化算出部５０１により、時刻τでの平均磁化［Ｍ_gａｖ］を算出する（ステップＳ１００５）。そして、シミュレーション装置２００は、磁界算出部５０２により、時刻τでの有効磁界［Ｈ_i］を算出する（ステップＳ１００６）。

このあと、シミュレーション装置２００は、時刻τを所定時間幅Δτ進め（ステップＳ１００７）、式（１５）のＬＬＧ方程式を用いて、磁化［Ｍ_i］の単位ベクトル［ｍ_i］の変化量を算出することで、磁化算出部５０３により、式（１６）を用いて磁化［Ｍ_i］を算出する（ステップＳ１００８）。そして、シミュレーション装置２００は、変化前後の磁化［Ｍ_i］を用いて、判断部５０４により磁化収束条件を満たすか否かを、メッシュｃｉごとに判断する（ステップＳ１００９）。

いずれか１つのメッシュｃｉでも磁化収束条件を満たさない場合（ステップＳ１００９：Ｎｏ）、ステップＳ１００５に戻って、ステップＳ１００８で更新後の磁化［Ｍ_i］により、平均磁化［Ｍａｖ］を算出することになる。一方、いずれのメッシュｃｉでも磁化収束条件を満たした場合（ステップＳ１００９：Ｙｅｓ）、図１１のステップＳ１１０１に移行する。

図１１は、図９に示した磁性体内の要素Ｃｇでのシミュレーション処理（ステップＳ９０３）の詳細な処理手順を示すフローチャート（その２）である。ステップＳ１００９：Ｙｅｓのあと、図１０において、シミュレーション装置２００は、［Ｍ］^jを、ステップＳ１００５で得られた最新の平均磁化［Ｍａｖ］に設定する（ステップＳ１１０１）。［Ｍ］^jは、式（１７）の摩擦項および式（１８）の慣性項で用いられる値である。後続のステップＳ１１０６でｊがインクリメントされるため、その後の有効磁界の算出（ステップＳ１００６）では、［Ｍ_g］^j-1となる。

このあと、シミュレーション装置２００は、静磁界算出部５０５により、式（２２）を用いて静磁界計算処理を実行する（ステップＳ１１０２）。そして、シミュレーション装置２００は、判定部５０６により、静磁界［Ｈｓ_g］が磁界収束条件を満たすか否かを判定する（ステップＳ１１０３）。満たさない場合（ステップＳ１１０３：Ｎｏ）、図１０のステップＳ１００４に戻る。このときの静磁界［Ｈｓ_g］と平均磁化［Ｍａｖ］の組み合わせは記憶領域５１１に格納されず、ヒステリシス曲線には反映されないことになる。

一方、満たした場合（ステップＳ１１０３：Ｙｅｓ）、シミュレーション装置２００は、時刻ｔｊでの静磁界［Ｈｓ_g］^jと平均磁化［Ｍａｖ］^jとして、記憶領域５１１に保持する（ステップＳ１１０４）。このあと、ｊ＞ｊｍａｘであるか否かを判断する（ステップＳ１１０５）。ｊｍａｘは変数ｊの最大値であり、ｔｊｍａｘがシミュレーション時間となる。ｊ＞ｊｍａｘでない場合（ステップＳ１１０５：Ｎｏ）、ｊをインクリメントして（ステップＳ１１０６）、図１０のステップＳ１００２に戻り、インクリメント後の時刻ｔｊで外部印加磁界［Ｈ_app］を再設定する。

また、ｊのインクリメントにより、ステップＳ１００１での最新の平均磁化［Ｍａｖ］^jは、平均磁化［Ｍａｖ］^j-1となり、平均磁化［Ｍａｖ］^j-1は平均磁化［Ｍａｖ］^j-2となる。これにより、摩擦項（式（１７））と慣性項（式（１８））の計算をおこなうことができることになる。一方、ｊ＞ｊｍａｘである場合（ステップＳ１１０５：Ｙｅｓ）、要素Ｃｇでのシミュレーション処理が終了するため、図９のステップＳ９０４に移行する。

図１２は、図１０に示した反転処理（ステップＳ１００４）の詳細な処理手順例を示すフローチャートである。まず、シミュレーション装置２００は、磁化番号ｉを＝１とし（ステップＳ１２０１）、ｉ＞ｉｍａｘであるか否かを判断する（ステップＳ１２０２）。ｉｍａｘは、磁化番号ｉの最大値である。ｉ＞ｉｍａｘでない場合（ステップＳ１２０２：Ｎｏ）、シミュレーション装置２００は、生成部８０１により、磁化［Ｍ_i］の容易軸ベクトル［ａ_ik］を生成する（ステップＳ１２０３）。ｋは、容易軸ベクトルの番号であり、一軸異方性の場合は、ｋ＝１，２となり、立方異方性の場合は、ｋ＝１，２，３，４，５，６となる。

そして、シミュレーション装置２００は、選択部８０２により、各ｋについて、磁気的エネルギーＥ_tot（［ａ_ik］）を算出する（ステップＳ１２０４）。シミュレーション装置２００は、選択部８０２により、最小の磁気的エネルギーＥ_tot（［ａ_ik］）を選択して、Ｅ_tot（［ａ_ikmin］）とする（ステップＳ１２０５）。容易軸ベクトル［ａ_ikmin］が、特定の容易軸ベクトルとなる。

つぎに、シミュレーション装置２００は、特定部８０３により、式（２５）を満たすか否かを判断する（ステップＳ１２０６）。満たす場合（ステップＳ１２０６：Ｙｅｓ）、シミュレーション装置２００は、特定部８０３により、反転角度θをθ＝１８０度とし、ピニングによるエネルギー障壁の高さに対応するパラメータΔＥをΔＥ_i ¹⁸⁰に設定する（ステップＳ１２０７）。そして、シミュレーション装置２００は、反転判断部８０４により、式（２０）に示す１８０度反転条件を満たすか否かを判断する（ステップＳ１２０８）。式（２０）に示す１８０度反転条件を満たす場合（ステップＳ１２０８：Ｙｅｓ）、シミュレーション装置２００は、反転部８０５により、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を−［ｍ_i］に更新して（ステップＳ１２０９）、ステップＳ１２１３に移行する。一方、式（２０）に示す１８０度反転条件を満たさない場合（ステップＳ１２０８：Ｎｏ）、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を更新せずに、ステップＳ１２１３に移行する。

また、ステップＳ１２０６において、式（２５）を満たさない場合（ステップＳ１２０６：Ｎｏ）、シミュレーション装置２００は、特定部８０３により、反転角度θをθ＝９０度とし、ピニングによるエネルギー障壁の高さに対応するパラメータΔＥをΔＥ_i ⁹⁰に設定する（ステップＳ１２１０）。そして、シミュレーション装置２００は、反転判断部８０４により、式（２１）に示す９０度反転条件を満たすか否かを判断する（ステップＳ１２１１）。式（２１）に示す９０度反転条件を満たす場合（ステップＳ１２１１：Ｙｅｓ）、シミュレーション装置２００は、反転部８０５により、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を特定の容易軸ベクトル［ａ_ikmin］に更新して（ステップＳ１２１２）、ステップＳ１２１３に移行する。一方、式（２１）に示す９０度反転条件を満たさない場合（ステップＳ１２１１：Ｎｏ）、磁化［Ｍ_i］の単位ベクトル［ｍ_i］を更新せずに、ステップＳ１２１３に移行する。

そして、シミュレーション装置２００は、ステップＳ１２１３において、ｉをインクリメントして（ステップＳ１２１３）、ステップＳ１２０２に戻る。ステップＳ１２０２において、ｉ＞ｉｍａｘである場合（ステップＳ１２０２：Ｙｅｓ）、反転処理（ステップＳ１００４）を終了し、ステップＳ１００５に移行する。

このように、本実施の形態によれば、方向性電磁鋼板や無方向性電磁鋼板でもみられる磁気特性の方向依存性を再現することができる。また、本実施の形態では、方向性電磁鋼板の圧延方向ＲＤだけでなく、直交方向ＴＤの磁気特性も再現することができる。また、式（２０）や式（２１）に示したように、簡単な演算により方向依存性を再現するため、方向依存性の再現に際して計算量の増加を抑制することができる。また、反転処理後の状態から磁化収束を判断したり、静磁界を計算したりするため、モータや変圧器などの電磁鋼板内で発生するヒステリシス損失を計算することができる。

また、本実施の形態によれば、摩擦を考慮して有効磁界［Ｈ_i］を求めることができるため、共鳴現象より経時的に変化する磁化［Ｍ_i］での磁化変化を阻止する挙動を再現することができる。したがって、その挙動に影響されたヒステリシス曲線が得られ、高周波での高抵抗磁性材料の磁性体についてヒステリシス損を高精度に求めることができる。

また、摩擦項とともに慣性項も加えた場合、摩擦および慣性を考慮して有効磁界［Ｈ_i］を求めることができるため、共鳴現象より経時的に変化する磁化［Ｍ_i］での磁化変化を阻止する挙動と、磁化変化速度を一定に保つ挙動と、を再現することができる。したがって、その挙動に影響されたヒステリシス曲線が得られ、高周波での高抵抗磁性材料の磁性体についてヒステリシス損失を高精度に求めることができる。

なお、上述した実施の形態では、シミュレーション装置２００において作成部５０８によりヒステリシス曲線を作成し、損失算出部５０９によりヒステリシス損失を算出することとしたが、作成部５０８および損失算出部５０９は、シミュレーション装置２００ではない他の装置に実行させてもよい。たとえば、シミュレーション装置２００では、格納部５０７により記憶領域５１１に、時刻ｔｊでの静磁界［Ｈｓ_g］^jと平均磁化［Ｍａｖ］^jとの組み合わせを格納して、作成部５０８および損失算出部５０９を有する他の装置に対し送信することとしてもよい。

上述した実施の形態に関し、さらに以下の付記を開示する。

（付記１）磁性体のシミュレーションプログラムにおいて、
コンピュータに、
前記磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成させ、
分割された領域の磁化ごとの磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択させ、
前記領域の磁化と、選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、前記領域内のピニングによるエネルギー障壁の高さに応じた磁化反転の反転角度を特定させ、
特定された反転角度により前記磁化を反転させる、
ことを特徴とする磁性体のシミュレーションプログラム。

（付記２）前記磁気的エネルギーの選択において、
前記コンピュータに、
前記磁気的エネルギー群の中から最小となる磁気的エネルギーを選択させることを特徴とする付記１記載の磁性体のシミュレーションプログラム。

（付記３）前記コンピュータに、
前記磁化と、前記特定の容易軸ベクトルと、に基づいて、前記ピニングによるエネルギー障壁の高さに対応するパラメータを特定させ、
前記磁化に関する前記領域内の磁気的エネルギーと、特定されたパラメータと、に基づいて、前記反転角度により前記磁化を反転させるか否かを判断させ、
前記磁化の反転において、
前記コンピュータに、
前記磁化を反転させると判断された場合、前記反転角度により前記磁化を反転させることを特徴とする付記１または２記載の磁性体のシミュレーションプログラム。

（付記４）前記磁化を反転させるか否かの判断において、
前記コンピュータに、
前記反転角度が１８０度である場合、前記磁化に関する前記領域内の磁気的エネルギーから前記磁化に対し方向が１８０度反転した磁化に関する前記領域内の磁気的エネルギーを引いた差分が、前記パラメータ以上であるか否かを判断させ、
前記磁化の反転において、
前記コンピュータに、
前記磁化を反転させると判断された場合、前記磁化を前記磁化に対し方向が１８０度反転した磁化に更新させることを特徴とする付記３記載の磁性体のシミュレーションプログラム。

（付記５）前記磁化を反転させるか否かの判断において、
前記コンピュータに、
前記反転角度が９０度である場合、前記磁化に関する前記領域内の磁気的エネルギーから前記特定の容易軸ベクトルに関する前記領域内の磁気的エネルギーを引いた差分が、前記パラメータ以上であるか否かを判断させ、
前記磁化の反転において、
前記コンピュータに、
前記磁化を反転させると判断された場合、前記磁化を前記特定の容易軸ベクトルに更新させることを特徴とする付記３記載の磁性体のシミュレーションプログラム。

（付記６）前記コンピュータに、
前記磁化の反転がおこなわれた後に、前記磁性体を構成する要素群のいずれかの要素から分割された各領域内の磁気的エネルギーから生成される磁界と、前記各領域の磁化を平均した平均磁化の変化を阻止する方向に働く磁化変化速度と、に基づく有効磁界を、領域ごとに算出させ、
前記領域ごとに算出された有効磁界と、前記領域ごとの磁化に基づいて、前記領域ごとに磁化変化量を求めて、前記領域ごとに前記変化後の磁化を算出させ、
前記領域ごとの変化前の磁化と前記変化後の磁化に基づいて、前記いずれかの要素における磁化が収束するかを判断させ、
前記いずれかの要素における磁化が収束すると判断された場合の前記領域ごとの平均磁化および前記平均磁化に基づく静磁界との組み合わせを記憶装置に格納させる、
ことを特徴とする付記１〜５のいずれか１項に記載の磁性体のシミュレーションプログラム。

（付記７）前記有効磁界の算出において、
前記コンピュータに、
前記領域ごとの磁化が変化する場合、前記領域内の磁気的エネルギーにより生成される磁界と、前記磁化変化速度と、前記磁化変化速度を一定に保つ方向に前記磁界を働かせる慣性と、に基づく有効磁界を、前記領域ごとに算出させることを特徴とする付記６記載の磁性体のシミュレーションプログラム。

（付記８）前記コンピュータに、
前記記憶装置に格納された前記平均磁化および前記静磁界の組み合わせ群から得られるヒステリシス曲線の面積からヒステリシス損失を算出させる、
ことを特徴とする付記６または７記載の磁性体のシミュレーションプログラム。

（付記９）磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成する生成部と、
分割された領域の磁化ごとの磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択する選択部と、
前記領域の磁化と、選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、前記領域内のピニングによるエネルギー障壁の高さに応じた磁化反転の反転角度を特定する特定部と、
前記特定部によって特定された反転角度により前記磁化を反転させる反転部と、
を有することを特徴とする磁性体のシミュレーション装置。

（付記１０）磁性体のシミュレーション方法において、
コンピュータが、
前記磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成し、
分割された領域の磁化ごとの磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択し、
前記領域の磁化と、選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、前記領域内のピニングによるエネルギー障壁の高さに応じた磁化反転の反転角度を特定し、
特定された反転角度により前記磁化を反転させる、
処理を実行することを特徴とする磁性体のシミュレーション方法。

２００ シミュレーション装置
２０２ 記憶装置
５００ 反転処理部
５０１ 平均磁化算出部
５０２ 磁界算出部
５０３ 磁化算出部
５０４ 判断部
５０５ 静磁界算出部
５０６ 判定部
５０７ 格納部
５０８ 作成部
５０９ 損失算出部
５１０ 出力部
８０１ 生成部
８０２ 選択部
８０３ 特定部
８０４ 反転判断部
８０５ 反転部

Claims (6)

1. 磁性体のシミュレーションプログラムにおいて、
   コンピュータに、
   前記磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成させ、
   分割された領域の磁化ごとに、生成された前記容易軸ベクトルそれぞれに磁化した場合の磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択させ、
   分割された領域の磁化ごとに、反転前の当該磁化と、前記容易軸ベクトルのうち選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、磁化反転の反転角度を特定させ、
   分割された領域の磁化ごとに、当該領域内のピニングによるエネルギー障壁の高さに応じて、当該磁化について特定された反転角度により当該磁化を反転させるか否かを判断させる、
   ことを特徴とする磁性体のシミュレーションプログラム。
2. 前記コンピュータに、
   分割された領域の磁化ごとに、当該磁化と、前記特定の容易軸ベクトルとの間の角度に基づいて、所定のパラメータ群の中から前記ピニングによるエネルギー障壁の高さに対応するパラメータを特定させ、
   前記磁化を反転させるか否かの判断において、
   前記コンピュータに、
   分割された領域の磁化ごとに、反転前後の当該磁化に関する前記領域内の磁気的エネルギーと、特定されたパラメータと、に基づいて、当該磁化について特定された前記反転角度により当該磁化を反転させるか否かを判断させ、
   前記コンピュータに、
   分割された領域の磁化ごとに、当該磁化を反転させると判断された場合、当該磁化について特定された前記反転角度により当該磁化を反転させることを特徴とする請求項１記載の磁性体のシミュレーションプログラム。
3. 前記磁化を反転させるか否かの判断において、
   前記コンピュータに、
   分割された領域の磁化ごとに、当該磁化について特定された前記反転角度が１８０度である場合、当該磁化に関する当該領域内の磁気的エネルギーから当該磁化に対し方向が１８０度反転した磁化に関する当該領域内の磁気的エネルギーを引いた差分が、特定された前記パラメータ以上であるか否かを判断させ、
   前記磁化の反転において、
   前記コンピュータに、
   分割された領域の磁化ごとに、当該磁化を反転させると判断された場合、当該磁化を当該磁化に対し方向が１８０度反転した磁化に更新させることを特徴とする請求項２記載の磁性体のシミュレーションプログラム。
4. 前記磁化を反転させるか否かの判断において、
   前記コンピュータに、
   分割された領域の磁化ごとに、当該磁化について特定された前記反転角度が９０度である場合、当該磁化に関する当該領域内の磁気的エネルギーから前記特定の容易軸ベクトルに関する当該領域内の磁気的エネルギーを引いた差分が、特定された前記パラメータ以上であるか否かを判断させ、
   前記磁化の反転において、
   前記コンピュータに、
   分割された領域の磁化ごとに、当該磁化を反転させると判断された場合、当該磁化を前記特定の容易軸ベクトルに更新させることを特徴とする請求項２記載の磁性体のシミュレーションプログラム。
5. 磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成する生成部と、
   分割された領域の磁化ごとに、生成された前記容易軸ベクトルそれぞれに磁化した場合の磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択する選択部と、
   分割された領域の磁化ごとに、反転前の当該磁化と、前記容易軸ベクトルのうち選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、磁化反転の反転角度を特定する特定部と、
   分割された領域の磁化ごとに、当該領域内のピニングによるエネルギー障壁の高さに応じて、前記特定部によって当該磁化について特定された反転角度により当該磁化を反転させるか否かを判断する反転部と、
   を有することを特徴とする磁性体のシミュレーション装置。
6. 磁性体のシミュレーション方法において、
   コンピュータが、
   前記磁性体を構成する要素群のいずれかの要素から分割された領域における容易軸ベクトルを生成し、
   分割された領域の磁化ごとに、生成された前記容易軸ベクトルそれぞれに磁化した場合の磁気的エネルギーを算出し、算出した磁気的エネルギー群の中から最大ではない磁気的エネルギーを選択し、
   分割された領域の磁化ごとに、反転前の当該磁化と、前記容易軸ベクトルのうち選択された磁気的エネルギーの場合の特定の容易軸ベクトルと、に基づいて、磁化反転の反転角度を特定し、
   分割された領域の磁化ごとに、当該領域内のピニングによるエネルギー障壁の高さに応じて、当該磁化について特定された反転角度により当該磁化を反転させるか否かを判断する、
   処理を実行することを特徴とする磁性体のシミュレーション方法。

Images