JP4508582B2

JP4508582B2 - ヒステリシス磁界解析法及びシステム

Info

Publication number: JP4508582B2
Application number: JP2003313124A
Authority: JP
Inventors: 健治宮田; 燦李; 芳英和田山; 弘隆竹島
Original assignee: Hitachi Ltd
Current assignee: Hitachi Ltd
Priority date: 2003-09-04
Filing date: 2003-09-04
Publication date: 2010-07-21
Anticipated expiration: 2023-09-04
Also published as: JP2005083764A

Description

本発明は、磁気回路のヒステリシス特性（以下、ヒステリシス磁界という。）の解析法及び装置に係り、具体的には、ヒステリシス磁界を高精度かつ高速に求める方法及びシステムに関する。

電気機器等の磁気回路を構成する磁性体は、必ずヒステリシス特性を有する。磁気回路のヒステリシス特性は、電気機器の動特性や損失に影響を及ぼすことから、設計等に際してヒステリシス磁界を解析することは重要である。

従来、ヒステリシス磁界の解析方法として知られている代表的な手法は、プライザッハモデルである（例えば、非特許文献1）。この方法は、磁気回路を構成する磁性体を多数の磁気双極子の集合体としてモデル化し、ミクロ的観点から解析する方法である。

Magnetic Hysteresis：E. Della Torre著、IEEE Press、New York、1999

モータ等の回転機や磁気共鳴撮像（ＭＲＩ）等の磁気回路においては、電磁鋼鈑、純鉄等の磁性体を使用するが、従来、設計段階でそれらの磁気回路のヒステリシス特性の影響は考慮されていなかった。

これらのヒステリシス磁界の解析に用いる従来のプライザッハモデルは、ヒステリシス特性の一面だけをモデル化するに止まっていて、種々のヒステリシス特性を正確に反映したモデルではない。例えば、磁界Ｈの変動幅が同一な場合、プライザッハモデルによれば、ヒステリシスのメジャーループ内のどの部分においても同一のマイナーループを形成することになる。しかし、実際の磁性体においては、そのような結果は観察されない。このように、今までのモデルはヒステリシス特性を正確に表現できないため、磁界解析にヒステリシス特性の影響を取り入れることができなかった。

一方、実際的な面であるが、磁気回路が複雑になるにつれて、人手による計算は殆ど不可能になり、計算機で計算しなければならないが、磁気回路の形状が複雑になればなるほど、計算機に計算させる計算要素の数は幾何級数的に拡散する。例えば、ＭＲＩのような形状が複雑な磁気回路においては、数百万の要素を必要とする。このような場合、解析には膨大なメモリと計算時間を要する。そこに、ヒステリシス特性を考慮すると、さらに膨大なメモリと計算時間を必要とし、計算が不可能なことになる。例えば、プライザッハモデルにおいては、そのモデルが最小限に必要とする二次元三角形の積分領域の境界条件を各要素が必ず記憶していなければならない。磁性体の要素数が数十万であっても、その境界条件を記憶するためのメモリを用意しなければならないことは、ＭＲＩのような巨大系では磁性体のために要素数に匹敵するメモリを用意しなければならないことになる。
このようなメモリと計算時間の制約を克服するためには、メモリを必要としないで、かつ高速の計算アルゴリズムが必要となる。

本発明は、磁気回路のヒステリシス特性を高精度かつ高速に求めることを課題とする。

上記課題を解決するために、本発明は、ヒステリシス曲線に関する現象論を直接反映させた関数論で表現することを基本とする。具体的には、解析対象の磁性体について実測したメジャーループのヒステリシスデータに基づいてヒステリシス関数を演算し、そのヒステリシス関数に基づいてヒステリシス磁界解析することを特徴とする。

すなわち、本発明のヒステリシス磁界解析法及びシステムは、解析対象の磁性体についてメジャーループのヒステリシスデータを予め計測し、該計測されたヒステリシスデータに基づいてマイナーループを含む前記磁性体のヒステリシス関数を求め、該求めたヒステリシス関数に基づいてヒステリシス磁界を解析することを特徴とする。

この場合において、前記ヒステリシスデータは、磁界強度と磁束密度又は磁界と磁化の空間データを用いることができる。また、前記ヒステリシス関数は、磁界の増加に対応する上昇曲線関数と、磁界の減少に対応する下降曲線関数とで記述され、前記上昇曲線関数と前記下降曲線関数は前記メジャーループの飽和点に対応する2つの特異点で接続されるものとすることができる。

また、ヒステリシス関数を求める際に、予め設定された少なくとも１つのパラメータ関数を用いることができる。また、ヒステリシス関数の演算は、ヒステリシス軌道を過去にさかのぼって最初に出会う分岐点と、現在のヒステリシス軌道が目指す分岐点の２個の分岐点間を結ぶ曲線関数でヒステリシス軌道を記述することができる。

ヒステリシス磁界解析は、有限要素法を用いて解析することができる。この場合において、コイル電流等の磁界発生ソースの変動に対して、前の演算周期で求めた透磁率分布を用いて近似解析し、ヒステリシス曲線の上昇か下降かを要素ごとに判別し、その結果を用いて前記ヒステリシス関数の演算により用意された前記ヒステリシス関数により微分透磁率を求め、求めた透磁率を用いて磁束密度の変化量を求めることができる。また、微分透磁率を求めることに代えて、磁性材料の磁化変化量を求め、求めた磁化変化率を用いて磁束密度の変化量を求めることができる。

また、本発明のヒステリシス磁界解析システムは、予め計測された解析対象の磁性体についてメジャーループのヒステリシスデータを入力し、該入力されたヒステリシスデータに基づいてマイナーループを含む前記磁性体のヒステリシス関数を演算により求め、該求めたヒステリシス関数に基づいてヒステリシス磁界を解析する処理をコンピュータに実行させるプログラムにより実現できる。また、そのプログラムを記録媒体に格納することにより、任意のコンピュータにより本発明のヒステリシス磁界解析システムを実現できる。

本発明によれば、磁気回路のヒステリシス特性を高精度かつ高速に求めることができる。

以下、本発明の実施の形態について説明する。図1は、本発明のヒステリシス磁界解析法の一実施例の手順を示すフローチャートである。図2は、実施例１のヒステリシス磁界解析法を適用してなる一実施の形態の装置である。図3は、本発明のヒステリシス磁界解析法の実施例２の手順を示すフローチャートである。
（基本原理１の説明）
ここで、本発明の実施例のヒステリシス磁界解析法を説明する前に、本発明のヒステリシス磁界解析法の原理について説明する。本発明は、磁界Ｈと磁化Ｍの空間においてメジャーループ内部の任意のヒステリシス曲線（マイナーループ）を記述するヒステリシス関数をＭ(Ｈ)としたとき、Ｍ(Ｈ)＝Ｍ_１＋（Ｍ_２−Ｍ_１）ｆ（Ｈ；Ｈ_１、Ｈ_２）とする。ここに、（Ｈ_１、Ｍ_１）(Ｈ_２、Ｍ_２)は、過去の履歴から決定される（Ｈ、Ｍ）空間上の両側の２つの分岐点である。また、関数ｆ（Ｈ；Ｈ_１、Ｈ_２）はＨ＝Ｈ_１のとき０で、Ｈ＝Ｈ_２のとき１となるように規格化された関数であり、実測されたメジャーループデータから同定できる。このヒステリシス関数Ｍ(Ｈ)を用いて、微分透磁率を順次更新しながら磁界解析することにより、順次（Ｈ、Ｍ）点を解析して高速高精度で消費メモリの少ないヒステリシス磁界解析を実現する。

すなわち、図４に示すメジャーループの計測データが与えられることを前提に、図５に示すマイナーループを作成することを特徴とする。以下に、このマイナーループの作成手順を説明する。ヒステリシス曲線を記述するヒステリシス関数の基本式は次式の数１で表せる。

ここで、Ｍｉはヒステリシスループにおける出発点の磁化を示し、Ｍｆは終着点の磁化を示す。Ｈ_ｉとＨ_ｆとの間には、Ｍ_ｉ＝ｆ_±（Ｈ_ｉ）、Ｍ_ｆ＝ｆ_±（Ｈ_ｆ）の関係がある。数１に基づくヒステリシスモデルは、ヒステリシスループを形成する両端の２つの分岐点(Ｈ_ｉ、Ｍ_ｉ)、(Ｈ_ｆ、Ｍ_ｆ)を確定点として固定し、その間をつなぐ曲線関数の形をしている。ヒステリシス軌道は、途中で新たな分岐点を形成して、軌道を分岐させた場合、分岐点を新たなペアに変更して、数１によってその後のヒステリシス軌道を記述することになる。これが、本発明におけるヒステリシスモデルの基本的な考え方である。数１を、別の言葉で表現すると、ヒステリシス軌道を過去にさかのぼって最初に出会う分岐点と、現在のヒステリシス軌道が目指す分岐点の計２個の分岐点間を結ぶ曲線関数でヒステリシス軌道を記述するということである。

ところで、図５に示したマイナーループは閉曲線で示しているが、実際のマイナーループは一般的に完全な閉曲線を形成するわけではない。動作点(Ｈ_ｉ、Ｍ_ｉ)から出発したヒステリシス軌道は下降曲線を描き、点(Ｈ_ｆ、Ｍ_ｆ)で反転して上昇曲線を描くが、完全には元の点(Ｈ_ｉ、Ｍ_ｉ) には戻らない。しかし、本発明では、マイナーループが閉ループを形成するか否かに関わらず、基本関数である数１でヒステリシス曲線を記述するものとする。

次に、ヒステリシス軌跡を描く手順について、図６を用いて説明する。図６（ａ）に示すように、最初は磁性体が消磁されている原点（0、0）から出発するとする。Ｈ＝０から磁界を増加させると、初期磁化曲線に沿って飽和磁化点(Ｈ_Ｓ、Ｍ_Ｓ)に向かう。このヒステリシス曲線を求めるために数１を用いる。原点（0、0）から出発する初期磁化曲線ｆ_０→Ｈｓ(Ｈ)は、次の数２で表せる。

ここで、磁界Ｈの増加率を図６(ｂ)に示す磁界Ｈ＝Ｈ_１で反転させると、この点はヒステリシス軌道における分岐点となる。この分岐点における磁化Ｍ_１の値は、数２から次の数３のように与えられる。

分岐点(Ｈ_１、Ｍ_１)で反転したＭ−Ｈカーブは、マイナス側の飽和点（−Ｈｓ、 −Ｍｓ）に向かう。その関数ｆ_{Ｈ１→Ｈｓ}(Ｈ)を数１を用いて表すと次の数４になる。

数４は、図６(ｂ)に破線で示す点(Ｈ_１、Ｍ_１)から飽和点（−Ｈｓ、 −Ｍｓ）に向かう曲線の関数を示している。この曲線上の磁界Ｈ_２で磁界の減少を止めた場合、Ｍ_２は数４から次の数５により求められる。

この点（Ｈ_２、Ｍ_２）から磁界の増加率を反転させて磁界が増加に転じると、点(Ｈ_１、Ｍ_１)に向かうマイナーループを作る曲線ｆ_{Ｈ２→Ｈ１}(Ｈ)を同様な方法で求めることが可能であり、数６のようになる。

このように、初期磁化曲線上の任意の点からマイナーループを作ることが可能である。例えば、図７に、上昇と下降の初期磁化曲線からなるマイナーループを作成した例を示す。このような方法を用いるとメジャーループの上昇及び下降曲線上に生じるマイナーループも生成することが可能である。

ところで、非特許文献１には、鉄と４―７９Ｍｏパーマロイの各種条件のマイナーループに関するR. M. Bozorthらによる実験結果が示されている。これによれば、いろいろな動作点におけるマイナーループの形状、傾き、あるいは入れ子状の多重構造をもつマイナーループの振る舞い等について測定されており、マイナーループの特色を良く表している。同様なマイナーループを本発明法のヒステリシス関数で表現したものを図８に示す。図８から判るように、本発明法を用いたマイナーループは、非特許文献１に記載された実測のマイナーループが有する特徴をよく再現している。すなわち、磁化Ｍの絶対値が小さいところのマイナーループは膨らみが比較的大きく、磁化Ｍの絶対値が大きくなって磁気飽和点に近づくに従って、マイナーループは扁平な形状になり、２つの分岐点を結ぶ直線の勾配も次第に小さくなる。

図９は、本ヒステリシスモデルを用いて、磁界Ｈの変動幅を同一にして異なる出発点から作成したマイナーループを示している。この場合、標準的な方法として古くから利用されているプライザッハモデルでは同じ大きさのマイナーループが作成されるが、本発明法では、同図に示すように、実際の現象を反映した、異なった大きさのマイナーループを作成することが可能である。
（基本原理２の説明）
上において、メジャーループの１本のデータのみを用いてマイナーループを作成する本発明法を説明した。ところが、物質によっては、上述の方法だけで表現可能なものだけではなく、ある物質では、実際のマイナーループと比べて誤差が大きい場合もある。また、高精度な計算を要求される場合がある。このような、要求を満たすためには上述の方法を拡張した次に述べるパラメータ関数を導入する必要がある。

その方法を、簡単に説明する。図１０にメジャーループとその内側に形成されたマイナーループ１がマイナーループ２を含む２つのマイナーループがあるとする。ここで、本発明のヒステリシス関数モデルでは、簡単のため、次に定義するマイナーループモデルを採用している。

マイナーループモデルの定義：図１０のような多重マイナーループにおいて、内側マイナーループ２が外側マイナーループ１の外に出るためには、外側マイナーループ１を形成する２個の分岐点のどちらかを通過しなくてはならない。また、その分岐点に達した時点でそのマイナーループ１の内部でたどったヒステリシス軌道に関する記憶（履歴）はすべて失われる。

このように定義されたマイナーループモデルは、内側のマイナーループ軌道が目指す点を、その外側のマイナーループがもつ２個の分岐点に固定している。また、あらゆるマイナーループは閉ループを形成するものとする。メジャーループは、一番大きなマイナーループとみなすことができ、メジャーループの特異点は、両側の２つの飽和点である。このモデルは、メジャーループの飽和点に達すると、それまでのヒステリシス軌道の記憶が消失するという現象を、全てのマイナーループに対しても拡張したモデルになっている。つまり、マイナーループの両側の特異点は、メジャーループの飽和点と同じ意味合いを持っている。こういう意味で、本発明のマイナーループモデルでは、マイナーループの２個の分岐点は数学的な特異点としての意味合いをもっている。

次に、さらに次の関数モデルを設定する。

（関数モデルの定義）
１つのループが与えられた場合、そのループの中のあらゆる曲線関数は、ループを表現する曲線関数ならびにループ両端の２個の特異点を用いて表現可能である。

ここで、図１０に示すループ０、１、２の上昇曲線の関数を、それぞれ、次のように設定する。

ｆ₊（Ｈ）：メジャーループ（ループ０）の上昇曲線関数
ｆ_１,+（Ｈ）：マイナーループ（ループ１）の上昇曲線関数
ｆ_２,+（Ｈ）：マイナーループ（ループ２）の上昇曲線関数
この関数モデルにおいて、ｆ₊（Ｈ）を用いて、ｆ_１,+（Ｈ）とｆ_２,+（Ｈ）を作成でき、ｆ_１,+（Ｈ）を用いてｆ_２,+（Ｈ）を作成できる。

これらの関数モデルに示した２個の特異点間の関係にはまだ自由度が残っており、この自由度について考察する。一般的にｆ₊（Ｈ）、ｆ_１,+（Ｈ）ならびにｆ_２,+（Ｈ）は、同一の微分関数を持たない。そこで、ループ０、１、２間に関して微分関数を同一にするある種の関数変換を考える。ここで、ループｉ
( ｉ＝０、１、２ )からループｊ ( ｊ＝０、１、２；ｊ＞ｉ
)へのヒステリシス軌道に関する上昇曲線の変換を、Ｔ［ｆ(Ｈ)；ｉ、ｊ］とおくと、ｆ_２,+（Ｈ）は次の数７のように変換される。

また、ｆ_２,+（Ｈ）はループ１の上昇曲線関数ｆ_１,+（Ｈ）を用いて次の数８のように変換できる。

同様にして、ｆ_１,+（Ｈ）は次の数９のように表現できる。

ここで、数７、数８、数９を用いると、次の数１０の関係式が導かれる。

数１０の関係式は、変換Ｔが群をなすことを意味する。次に、変換Ｔに関する群の例を５つ示す。
（１）Translation（移動）
関数ｆ(Ｈ)をＨ軸方向にηだけ移動する変換群Ｔ：ｆ(Ｈ)→ｆ(Ｈ-η) に関して、数１１が成り立つと、数１２が成り立つ。

（２）Dilation（引き伸ばし）
関数ｆ(Ｈ)をＨ軸方向にα倍引き伸ばす変換群Ｔ：ｆ(Ｈ)→ｆ(Ｈ/α) に関して、数１３が成り立つとすると、数１４が成り立つ。

（３）Power‐f（歪ませる）
関数ｆ(Ｈ)をｆ(Ｈ)γに歪ませる変換群Ｔ：ｆ(Ｈ)→ｆ(Ｈ)γに関して、数１５を定義すると、数１６となる。

（４）Power‐H（歪ませる）
関数ｆ(Ｈ)をｆ(Ｈλ)に歪ませる変換群Ｔ：ｆ(Ｈ)→ｆ(Ｈλ)に関して、数１７を定義すると、数１８となる。

（５）Rotation（回転）
関数ｆ(Ｈ)を角度θ回転する変換群Ｔ：Ｈ_ｍ＋ｊｆ_ｍ(Ｈ)→(Ｈ_ｍ＋ｊｆ_ｍ(Ｈ))ｅ^ｊθ に関して、数１９を定義する。

ここで、Ｈ_ｍはＨをその最大値で正規化した無次元数で、ｆ_ｍ(Ｈ)はｆ(Ｈ)をその最大値で正規化した無次元数である。数１９の定義式を用いると、数２０の関係になる。

以上、示した５種の変換群は、いずれも数１０を満たしている。ここで、これら変換群を複数重ね合わせた合成変換群も数１０を満たす必要がある。そのためには、いずれの変換群も互いに可換な群でなければならない。ところが、Ｍ=ｆ(Ｈ)として、全体座標系（Ｈ、Ｍ）の空間で単純に上記変換を施した場合、上記変換群は互いに非可換な群になってしまう。互いに可換な群にするためには、まず、変換を施す前の状態における座標系（Ｈ、Ｍ）を局所座標系としてこの関数系に固定する。translationの変換後も、rotation、dilation、power‐H、power‐fはいずれもこの局所座標系のもとに実施すればよい。この点に気をつければ、変換の順番は任意でよい。このため、dilation、power‐H、power‐fを先に実施して、つぎにrotation、最後にtranslationを施せば、確実に合成可換群を作成できる。

上記の変換群に関わる５個のパラメータα、γ、λ、θ、ηを内側マイナーループの出発点における磁化Ｍ_ｉの関数として数値解析的に求めることを考える。Ｍ_ｉから始まるマイナーループに関するヒステリシス関数Ｆ(Ｈ、Ｈ_ｉ、Ｍ_ｉ、Ｈ_ｆ、Ｍ_ｆ、
α、γ、λ、θ、η：ｆ_＋)と実測データｆ(Ｈ_ｉ、Ｍ_ｉ)との距離Ｄisを次の数２１で定義する。

数２１を用いて、距離Ｄisが最小になるように、各Ｍ_ｉにおいて、α、γ、λ、θ、ηの数値を決定する。こうして、５つのパラメータ関数α(Ｍ_ｉ)、γ(Ｍ_ｉ)、λ(Ｍ_ｉ)、θ(Ｍ_ｉ)、η(Ｍ_ｉ)が求まる。このようにして求めたパラメータ関数に基づいて、各マイナーループのＭｉに対応して、参照関数であるメジャーループ関数を変化処理し、その関数をマイナーループを記述する基本式に用いる。

なお、ここでは、例として５個のパラメータ関数を示したが、有限要素法では磁気ベクトルポテンシャルＡを未知数にして、磁束密度Ｂが求まることがしばしばである。この場合、上記のヒステリシスモデルを使うためには、Ｈ‐Ｂ空間ではなく、Ｂ‐Ｈ空間でヒステリシスが表現されなければならない。すなわち、Ｍ＝ｆ(Ｈ)の関係からＨ＝ｆ^−１(Ｍ)の関係に変換して用いなければならない。あるいは、はじめから、Ｈ＝ｇ(Ｍ)となるような関係式を導いておく必要がある。

次に、本発明のヒステリシスモデルを有限要素法に組み込む方法について説明する。有限要素法で磁界解析する場合の一般的な方法として、磁気ベクトルポテンシャルＡ、電気スカラポテンシャルφを用いたＡ−φ法で記述すると、変位電流の影響が無視できる準静的な場の基礎方程式は、数２２となる。

ここで、μは透磁率、σは導電率、Ｊ_０はコイルの電流密度ベクトル、Ｍ_０は磁石の残留磁束密度ベクトルである。数２２の右辺の項が場を励磁するためのソース項である。

数２２において、磁石の残留磁束密度ベクトルＭ_０は一定であるから、コイルの電流密度ベクトルＪ_０が、Ｊ_０＋ΔＪ_０に微小変化した場合をソース項の変化として考える。ヒステリシス磁界解析では、磁界Ｈ、磁化Ｍの(Ｈ、Ｍ)空間において、上昇曲線に沿うか、あるいは下降曲線に沿うかで、振る舞いが大きく異なる。そのため、Ｊ_０が微小変化してＪ_０＋ΔＪ_０になった場合、(Ｈ、Ｍ)空間で上昇曲線側か下降曲線側かを早めに判別しておくことが重要である。

そこで、数２２を演算周期の時刻ｔ＝ｎΔｔ（ｎ＝１、２、３、・・・)における方程式とみなし、次の時刻ｔ＝（ｎ＋１）Δｔにおける方程式として、数２３を考える。数２３におけるＡ’、φ’、μ’は、それぞれ数２４のとおりである。また、ΔＡ、Δφ、Δμは、それぞれ時間がΔｔ経過したことに伴うＡ、φ、μの微小変化量を表す。

まず、近似的にμ’＝μとおいて、線形磁界解析を実行する。時刻ｔ＝（ｎ＋１）Δｔにおける磁束密度の変化量ΔＢ_ｎ＋１をΔＢ_ｎ＋１＝ｒｏｔ(ΔＡ)により求める。次に、要素分割した各要素について上昇曲線か下降曲線かを判別する。この判定法は、ヒステリシスを示す全領域（全要素）をＩＤ＝１としておき、磁束密度の大きさ｜Ｂ｜が小さくなって、ある時点で磁束密度ベクトルが反対方向に向いた要素では、ＩＤ＝−１とする。また、同じような現象がおきたら、ＩＤ＝１にする。つまり、磁束密度ベクトルの向きが変わるたびにＩＤの符号を反転する。そして、ΔＢ・ＩＤ＞０のとき上昇、ΔＢ・ＩＤ＜０のとき下降と判定する。具体的には、ｔ＝ｎΔｔからｔ＝（ｎ＋１）Δｔの時間経過に伴う磁束密度の変動量ΔＢ_ｎ＋１を求め、そのときのＩＤとの積の正負によって（Ｈ、Ｍ）空間におけるヒステリシス軌跡が上昇方向に向かうのか、それとも下降方向に向かうのかを判定する。

上昇曲線か下降曲線かの判別の後、非線形磁界解析に移行する。ヒステリシス解析において、上昇曲線あるいは下降曲線をそのままたどっていく場合には、通常のニュートン・ラプソン（Newton‐Raphson）法で問題なく解ける。しかし、分岐点を形成してヒステリシス軌跡が反転した場合、取り扱いが厄介になる。分岐点の前後でＢＨカーブの勾配が特異的に変化するため、ニュートン・ラプソン法をそのまま適用するには問題があるので、ここでは単純反復法で解を求める。その代わり、磁気ベクトルポテンシャルＡそのものに関する数２３を直接解くのではなく、ソースの変化量ΔＪ_０に対する場の変化量ΔＡを求め、その結果として、Ａ＋ΔＡを求める。

まず、磁界ベクトルの変化量ΔＨ、磁化ベクトルの変化量ΔＭに対して、数２５の関係にあるから、微分磁化率χｄが近似的に確定する。

図１０には、上昇曲線をたどってきた時刻ｔ＝ｎΔｔの動作点(Ｈ、Ｍ)と時刻ｔ＝（ｎ＋１)Δｔにおける動作点(Ｈ、Ｍ)を上昇曲線上に沿う場合と反転して下降曲線上に沿う場合の２ケースについて示してある。この図が示すように、上昇曲線側か下降曲線側かで、微分磁化率χｄは大きく異なる。

ここで、図１１に、本発明におけるヒステリシスモデルを有限要素法にて解析する場合に有用になる磁束密度Bを独立変数とした（B、H）空間でのヒステリシス曲線を示す。まず、磁束密度Ｂの変化量をより正確に求めるための基礎式を導いておく。磁化をＭとすると、数２６の関係から、数２７が成立する。つまり、上述で求めたヒステリシス関数は、Ｍを表現する式になっている。そこで、ＨとＢ、Ｍの関係を示す数２６により、本発明におけるヒステリシスモデルを有限要素法に結合する。

数２７に数２５を代入すると、数２８が得られる。

さらに、数２５より、数２９が得られる。

ここで、数２６にrotationを施して、数３０を得る。

数３０にファラデー（Faraday）の式である数３１を代入すると、数３２が得られる。

数３２の変化分をとると、数３３となる。

数３３に数２９を代入して式を整理すると、数３４が得られる。

ここで、μｄは微分透磁率であり、微分磁化率χｄと数３５の関係にある。

なお、微分透磁率μｄや微分磁化率χｄは、求めようとしている時刻ｔ＝(ｎ＋1)Δｔにおける値ではなく、時刻ｔ＝(ｎ＋１／２)Δｔ近傍における値である点に注意しなければならない。ヒステリシスのない通常の非線形磁界解析において、透磁率μが時刻ｔ＝(ｎ＋1)Δｔにおける値になるように、ニュートン・ラプソン法でセルフコンシステントに解を求める場合と事情が異なる。数２５に示した磁界の変化量ΔＨと磁化の変化量ΔＭの関係式に示した勾配μ_０χｄは、図１２に示す三角形の斜辺の勾配に相当する。

また、別の解析法として、数３４に代えて、数３３を変形した次式の数３６を用いて解析することもできる。

ここで、磁気ベクトルポテンシャルＡの変化量ΔＡに関する基本式である数３４は、左辺第２項にＡの時間に関する偏微分項を含んでいる。これは次のように取り扱う。ニューマーク（Newmark）のβ法を用いると、次式の数３７になる。

数３７において、βは0≦β≦１なる固定定数であり、安定した時間応答解析のために１／２≦β≦１なるβを用いる。ここで、ΔＡ(ｔ)＝０であり、ΔＡ(ｔ＋Δt)はこれまでの式の展開と同様に単にΔＡとおき、数３４に代入すると、次式の数３８が得られる。

数３７を用いて磁気ベクトルポテンシャルの変化量ΔＡを求め、磁束密度変化量ΔＢを更新する。新たなΔＢに対して、微分透磁率μｄを更新し、再度、数３７を解く。このような操作を、ΔＢが収束するまで繰り返す。ΔＢが収束後、透磁率分布が求まり、時刻ｔ＝(ｎ＋1)Δｔでの解析は、基本的に終了する。

ところで、毎回、場の変化量のみを求めるとなると、数値誤差が積算する恐れがあるため、求めた透磁率分布μ’＝μ＋ΔμとＡ’=Ａ+ΔＡが数２３を満足しているかどうかをチェックし、必要ならば、Ａ’=Ａ＋ΔＡならびにμ’＝μ＋Δμを補正する。解は真の解の近傍にあるため、この補正のための計算は短時間で済む。

次に、上述した本発明のヒステリシス磁界解析法の具体的な一実施例を図１に示したフローチャートを用いて説明する。磁界解析手法としては、一般的な有限要素法を適用する。図２は、本発明のヒステリシス磁界解析法を適用してなる一実施例のヒステリシス磁界解析システムの構成図である。本システムは、図示のように、計算機３０、表示装置３１、記憶媒体３２から構成される。なお、記憶媒体３２は計算機３０の内部に設けることができる。

図１のヒステリシス磁界解析法は、数３４を基本解析式とするフローチャートである。本実施例による演算を実行する前に、磁性体と空間と含む磁気回路について解析対象の領域を設定し、その領域を有限個の微小な要素に分割する。この分割法は、周知のように、磁束密度の分布が大きく変化する領域は細かく分割し、あまり変化がないところは粗めに分割する。分割された要素の全てについて図１の演算処理を行う。本実施例では、演算周期の時間刻み幅をΔｔとし、図１は第（ｎ＋１）番目の時刻ｔ＝（ｎ＋１）Δｔにおける磁界解析の流れを示している。ただし、ここで、以下のことを念頭においておく必要がある。すなわち、磁性体の透磁率は一般的に磁界強度によって変化するということである。ヒステリシスを伴う磁界解析を実施する場合、磁性体の透磁率分布はその後のヒステリシスの振舞いに影響を与える。また、ヒステリシス磁界解析では、磁束密度の変化量はヒステリシス軌跡における分岐発生の有無の判定に重要である。このため、各時間ステップにおいて、磁性体の透磁率分布ならびに磁束密度の変化量を計算機３０の記憶媒体３２に格納しておくことが重要である。

図１の第１ステップ１１において、第ｎ番目の時刻ｔ＝ｎΔｔで計算機３０の記憶媒体３２に格納していた透磁率データを呼び出し、そのデータ並びに数２２（又は、数２３）を用いて、時刻ｔ＝（ｎ＋１）Δｔにおける透磁率固定の線形磁界解析を実施する。

第２ステップ１２において、ヒステリシス軌跡が上昇方向に向かうのか、それとも下降方向に向かうのかを判定する。この判定法は、ヒステリシスを示す全領域（全要素）をＩＤ＝１としておき、磁束密度の大きさ｜Ｂ｜が小さくなって、ある時点で磁束密度ベクトルが反対方向に向いた要素では、ＩＤ＝−１とする。また、同じような現象がおきたら、ＩＤ＝１にする。つまり、磁束密度ベクトルの向きが変わるたびにＩＤの符号を反転する。そして、ΔＢ・ＩＤ＞０のとき上昇、ΔＢ・ＩＤ＜０のとき下降と判定する。具体的には、ｔ＝ｎΔｔからｔ＝（ｎ＋１）Δｔの時間経過に伴う磁束密度の変動量ΔＢ_ｎ＋１を求め、そのときのＩＤとの積の正負によって（Ｈ、Ｍ）空間におけるヒステリシス軌跡が上昇方向に向かうのか、それとも下降方向に向かうのかを判定する。

第３ステップ１３において、上昇方向か下降方向が確定しているので、用いるヒステリシス関数（数１）を確定できる。つまり、上昇の場合は数６の式で表されるヒステリシス関数を用いる。また、下降の場合は、数４の式で表されるヒステリシス関数を用いる。なお、それらのヒステリシス関数において、（Ｈｉ、Ｍｉ）、（Ｈｆ、Ｍｆ）は、与えられたメインループの計測データ等の状況に応じて変わる。確定したヒステリシス関数と第２ステップで近似的に求めた磁束密度の変動量ΔＢ_ｎ＋１から数２５、数３５に基づき、微分透磁率μｄを求め、ヒステリシスを呈する磁性材料の各要素に割り当てる。次いで、数３４を用いて、磁気ベクトル変化量ΔＡを求め、磁束密度変化量ΔＢを求める。

第４ステップ１４において、磁束密度変化量ΔＢが収束したか否か判定する。磁束密度変化量ΔＢが収束していなければ、第４ステップ１４において、その磁束密度変化量ΔＢと用意したヒステリシス関数（数４、数６）から微分透磁率μを更新して、再度、第３ステップ１３で磁束密度変化量ΔＢを求める。

このようにして、第３ステップ１３および第４ステップ１４の操作を磁束密度変化量ΔＢが収束するまで反復する。磁束密度変化量ΔＢが収束した場合は、第５ステップ１５において、ｔ＝（ｎ＋１）Δｔにおける磁性体の磁束密度変化量ΔＢ_ｎ＋１ならびに各要素の新たな透磁率μ＝Ｂ／Ｈを計算機３０の記憶媒体３２に格納し、次の時刻ｔ＝(ｎ＋２) Δｔの解析に移行する。

図３に、本発明のヒステリシス磁界解析法の具体的な他の実施例のフローチャートを示す。本実施例のヒステリシス磁界解析法は、図２に示したヒステリシス磁界解析システムに適用できることは言うまでもない。本実施例は、数３６を基本解析式とするヒステリシス磁界解析法であり、時間刻み幅を実施例１と同じくΔｔとして、第（ｎ＋１）番目の時刻ｔ＝（ｎ＋１）Δｔにおける磁界解析工程を順を追って説明する。

第１ステップ２１において、計算機３０の記憶媒体３２に格納していた第ｎ番目の時刻ｔ＝ｎΔｔにおける透磁率データを呼び出し、そのデータならびに数２３を用いて、時刻ｔ＝（ｎ＋１）Δｔにおける透磁率固定の線形磁界解析を実施する。

第２ステップ２２において、ｔ＝ｎΔｔからｔ＝（ｎ＋１）Δｔの時間経過に伴う磁束密度の変動量ΔＢ_ｎ＋１を求め、図１の第２ステップ１２と同様に、ヒステリシス軌跡が上昇方向に向かうのか、それとも下降方向に向かうのかを判定する。

第３ステップ２３において、上昇方向か下降方向が確定しているので、図１の第３ステップと同様に、用いるヒステリシス関数（数１、数４、数６）を確定できる。確定したヒステリシス関数と第２ステップで近似的に求めた磁束密度の変動量ΔＢから磁化変化量ΔＭを求め、ヒステリシスを呈する磁性材料の各要素に割り当てる。次いで、数３６を用いて、磁気ベクトル変化量ΔＡを求め、磁束密度変化量ΔＢを求める。

第４ステップ２４において、磁束密度変化量ΔＢが収束したかどうかを判定し、収束していなければ、第６ステップ２５において、磁束密度変化量ΔＢと用意したヒステリシス関数から磁化変化量ΔＭを更新する。そして、第３ステップ２３に戻って、更新した磁化変化量ΔＭをヒステリシスを呈する磁性材料の各要素に割り当て、再度、数３６を用いて、磁気ベクトル変化量ΔＡを求め、磁束密度変化量ΔＢを求める。以上の過程を磁束密度変化量ΔＢが収束するまで繰り返す。

第４ステップ２４において磁束密度変化量ΔＢが収束したら、第７ステップ２６において、ｔ＝（ｎ＋１）Δｔにおける磁性体の磁束密度変化量ΔＢ_ｎ＋１ならびに各要素の新たな透磁率μ＝Ｂ／Ｈを計算機３０の記憶媒体３２に格納し、次の時刻ｔ＝(ｎ＋２) Δｔでの解析に移行する。

ここで、本発明をＭＲＩのヒステリシス磁界解析に適用した実施例を、図１３〜図１５を参照して説明する。図１３は、永久磁石を用いた開放型MRIの断面図を示す。ノン実施例の開放型ＭＲＩは、円盤状の永久磁石を上下に対抗して配置して構成されている。ＭＲＩは、主にヨーク１、永久磁石２、ポールピース３、半径及び垂直軸の原点（ｒ=０、ｚ=０）を含む有限領域に、傾斜磁場を形成するための傾斜磁場コイル４で構成されている。ＭＲＩでは均一静磁場に傾斜磁場を重畳させてｒ=0、ｚ=0の原点を含む有限領域の映像を作り出す。ここで、均一静磁場は永久磁石２ならびにヨーク１、ポールピース３で形成される磁気回路を通して形成される。この均一静磁場の安定性が、映像の分解能と直接関係するため最も重要な要素である。

図１５に示すように、均一静磁場を作るためにポールピース３の表面には幾つかの突起部５を持つような形になる。有限要素法で磁界解析する場合、磁束密度の飽和領域あるいはそれに近いところを表示するように入力すると、磁界解析により、例えば図１３の突起部５が特定される。特定した場所においてユーザが、例えば図１３のように領域６、領域７と名前を付け、マーキングする。入力手段から傾斜磁場コイル４にかかる電流の時間変化、例えば図１４のように入力し、ヒステリシス磁界解析を実行するように入力して、有限要素法による磁界解析を実行する。そして、解析結果を（Ｈ、Ｂ）空間におけるヒステリシス曲線で表示するように入力すると、図１５のような結果が得られる。

図１５において、領域６では、永久磁石２による磁界Ｈ０と傾斜磁場コイル４による磁界ＨＧＣが同じ向きに働くとすれば、反対側の領域７では、永久磁石２による磁界Ｈ０と傾斜磁場コイル４による磁界ＨＧＣが反対向きに働く。このため、ヒステリシス曲線は、領域６では、（Ｈ、Ｂ）空間の上側に小さなマイナーループを形成するが、領域７では下側に比較的大きなマイナーループを形成する。このようにヒステリシスループは、均一磁場の生成において障害となる。このため、本実施例によれば、ヒステリシス曲線を表示することにより、均一磁場を乱す場所を特定することができるという効果がある。また、ヒステリシスによる影響と渦電流による影響を解析上分離でき、ＭＲＩの映像に与える影響を原因別に探求できるという効果もある。

以上説明したように、本発明に関るヒステリシス磁界解析法によれば、ヒステリシス現象をほぼ忠実に再現できるため、従来法に比べて高精度な解析が可能になり、加えて磁気双極子集合モデルを使わずにマクロ関数で表現するため、計算機コアのメモリ量が少なくてすみ、かつ高速にヒステリシス解析ができるという効果がある。

本発明の実施例１のヒステリシス磁界解析法の手順を示すフローチャートである。本発明のヒステリシス磁界解析法を適用してなる一実施形態のヒステリシス磁界解析システムの構成図である。本発明の実施例２のヒステリシス磁界解析法の手順を示すフローチャートである。本発明の原理を説明するためのヒステリシス特性のメジャーループを示す図である。本発明の原理を説明するためのヒステリシス特性のマイナーループを示す図である。本発明により求めるマイナーループの例を説明する図である。本発明に係るヒステリシスモデルを用いて作成した上昇と下降の初期磁化曲線からなるマイナーループを示す図である。本発明に係るヒステリシスモデルを用いて作成した初期磁化曲線からのマイナーループを示す図である。本発明に係るヒステリシスモデルを用いて作成した磁界変動幅が同一で出発点が異なったマイナーループを示す図である。本発明のヒステリシス磁界解析法を説明するためのマイナーループを示す図である。本発明のヒステリシス磁界解析法を有限要素法にて実現する場合に有用になる磁束密度Bを独立変数とした（B、H）空間でのヒステリシス曲線を示す図である。ヒステリシス磁界解析法を説明する図である。本発明の解析対象の一実施例の永久磁石型MRIの縦断面図を示す図である。図１３の永久磁石型MRIの傾斜磁場コイルに流すコイル電流の時間変化を示す図である。図１３の永久磁石型MRIの解析結果の一例を示す図である。

符号の説明

３０計算機
３１表示装置
３２記憶媒体

Claims

解析対象の磁性体についてメジャーループのヒステリシスデータを予め計測し、該計測されたヒステリシスデータに基づいてマイナーループを含む前記磁性体のヒステリシス関数を演算により求め、該求めたヒステリシス関数に基づいてヒステリシス磁界を解析する方法であって、
前記ヒステリシスデータは、磁界強度と磁束密度又は磁界強度と磁化の空間データであり、
前記ヒステリシス関数は、磁界強度を横軸、磁束密度又は磁化を縦軸として表したヒステリシス軌道において、磁界強度の増加に対応する上昇曲線関数と、磁界強度の減少に対応する下降曲線関数とで記述され、前記上昇曲線関数と前記下降曲線関数は前記メジャーループの飽和点に対応する２つの特異点を共有するとともに、ヒステリシス軌道を過去にさかのぼって最初に出会う分岐点と、現在のヒステリシス軌道が目指す分岐点の２個の分岐点間を結ぶ曲線関数で前記ヒステリシス軌道が記述されてなるヒステリシス磁界解析法。
解析対象の磁性体についてメジャーループのヒステリシスデータを予め計測し、該計測されたヒステリシスデータと予め設定された少なくとも１つのパラメータ関数とに基づいてマイナーループを含む前記磁性体のヒステリシス関数を求め、該求めたヒステリシス関数に基づいてヒステリシス磁界を解析する方法であって、
前記ヒステリシスデータは、磁界強度と磁束密度又は磁界強度と磁化の空間データであり、
前記ヒステリシス関数は、磁界強度を横軸とし磁束密度又は磁化を縦軸として表したヒステリシス軌道において、磁界強度の増加に対応する上昇曲線関数と、磁界強度の減少に対応する下降曲線関数とで記述され、前記上昇曲線関数と前記下降曲線関数は前記メジャーループの飽和点に対応する２つの特異点を共有するとともに、ヒステリシス軌道を過去にさかのぼって最初に出会う分岐点と、現在のヒステリシス軌道が目指す分岐点の２個の分岐点間を結ぶ曲線関数でヒステリシス軌道が記述されてなるヒステリシス磁界解析法。
前記ヒステリシス磁界解析は、複数の要素に分割された解析空間を有限要素法を用いて解析するものとし、磁界発生コイル電流の変動に対して、要素ごとに求めた透磁率分布を用いて近似解析し、ヒステリシス曲線の上昇曲線側か下降曲線側かを要素ごとに判別し、その判別結果を用いて前記ヒステリシス関数の演算により求められた前記ヒステリシス関数により微分透磁率を求め、求めた透磁率を用いて磁束密度の変化量を求めることを特徴とする請求項1又は２に記載のヒステリシス磁界解析法。
前記ヒステリシス磁界解析は、複数の要素に分割された解析空間を有限要素法を用いて解析するものとし、磁界発生コイル電流の変動に対して、要素ごとに求めた透磁率分布を用いて近似解析し、ヒステリシス曲線の上昇曲線側か下降曲線側かを要素ごとに判別し、その判別結果を用いて前記ヒステリシス関数の演算により求められた前記ヒステリシス関数から磁性材料の磁化変化量を求め、求めた磁化変化量を用いて磁束密度の変化量を求めることを特徴とする請求項1又は２に記載のヒステリシス磁界解析法。
解析対象の磁性体について予め計測したメジャーループのヒステリシスデータを入力する手段と、該手段により入力されたヒステリシスデータに基づいてマイナーループを含む前記磁性体のヒステリシス関数を演算する手段と、該手段により求められたヒステリシス関数に基づいてヒステリシス磁界を解析する手段とを備えてなり、
前記ヒステリシスデータは、磁界強度と磁束密度又は磁界強度と磁化の空間データであり、
前記ヒステリシス関数は、磁界強度の増加に対応する上昇曲線関数と、磁界強度の減少に対応する下降曲線関数とで記述され、前記上昇曲線関数と前記下降曲線関数は前記メジャーループの飽和点に対応する２つの特異点を共有するとともに、ヒステリシス軌道を過去にさかのぼって最初に出会う分岐点と、現在のヒステリシス軌道が目指す分岐点の２個の分岐点間を結ぶ曲線関数でヒステリシス軌道が記述されてなるヒステリシス磁界解析システム。
前記ヒステリシスデータを入力する手段は、前記ヒステリシスデータとともに、少なくとも１つのパラメータ関数を入力することを特徴とする請求項５に記載のヒステリシス磁界解析システム。
前記ヒステリシス磁界解析手段は、複数の要素に分割された解析空間を有限要素法を用いて解析するものとし、磁界発生コイル電流の変動に対して、要素ごとに求めた透磁率分布を用いて近似解析し、ヒステリシス曲線の上昇曲線側か下降曲線側かを要素ごとに判別し、その判別結果を用いて前記ヒステリシス関数演算手段により用意された前記ヒステリシス関数により微分透磁率を求め、求めた透磁率を用いて磁束密度の変化量を求めることを特徴とする請求項５又は６に記載のヒステリシス磁界解析システム。
前記ヒステリシス磁界解析手段は、複数の要素に分割された解析空間を有限要素法を用いて解析するものとし、磁界発生コイル電流の変動に対して、要素ごとに求めた透磁率分布を用いて近似解析し、ヒステリシス曲線の上昇曲線側か下降曲線側かを要素ごとに判別し、その判別結果を用いて前記ヒステリシス関数演算手段により用意された前記ヒステリシス関数から磁性材料の磁化変化量を求め、求めた磁化変化量を用いて磁束密度の変化量を求めることを特徴とする請求項５又は６に記載のヒステリシス磁界解析システム。