JP2005207900A - 磁界解析法，磁界解析プログラム及び磁界解析プログラムを記録した記録媒体 - Google Patents

Abstract

【課題】
コイル電流の変化に伴う非定常磁界解析を高速処理する解析法を提供すること。
【解決手段】
本発明に関る磁界解析法では、磁気ベクトルポテンシャルＡならびに電気スカラーポテンシャルφそのものを解くのではなく、それらの変化量を逐次解いていく。また、微分透磁率と変動場を交互に修正しながら計算して収束解を得る。さらに、微分透磁率を用いた変動場に関する方程式を有限要素法で解く。それにより、コイル電流の変化に伴う非定常磁界解析を高速処理することが可能となる。
【選択図】図１

Description

本発明は、磁性材料の非線形磁界解析の方法，磁界解析プログラム及び磁界解析プログラムを記録した記録媒体に関するものである。

従来の代表的な非線形磁界解析法として、有限要素法による解析法があり、ＩＣＣＧ法による反復解法ならびに透磁率を逐次修正するNewton−Raphson法を併用している。例えば、中田高義，高橋則雄共著：「電気工学の有限要素法」森北出版（１９８６）（以下、非特許文献１と呼ぶ）に示されている。

中田高義，高橋則雄共著：「電気工学の有限要素法」森北出版 （１９８６）

モータ，回転機，医療用ＭＲＩ等では、その磁気回路において電磁鋼鈑，純鉄等の磁性体を使い、その磁性体は透磁率が磁束密度の大きさに依存するといった非線形性を有している。このような非線形磁界を解く場合、代表的な解法として、磁気ベクトルポテンシャルＡならびに電気スカラーポテンシャルφを用いるＡ−φ法では、変位電流の影響が無視できる準静的な場の基礎方程式である

をＩＣＣＧ法（不完全コレスキー型前処理付き共役勾配法）等による反復解法ならびに透磁率を逐次修正するNewton−Raphson 法を併用して解いている。ここに、μは透磁率、σは導電率、Ｊ₀はコイルの電流密度ベクトル、Ｍ₀は磁石の残留磁束密度ベクトルである。なお、磁石に関しては、便宜上、Ｂ＝μＨ＋Ｍ₀ とした。また、磁気ベクトルポテンシャルＡのみを用いるＡ法では、

の方程式を解く。しかし、いずれの方法でも、大規模解析になればなるほど、ＩＣＣＧ法等による反復解法における反復回数や、Newton−Raphson 法における反復回数が増える。非定常解析においては、前時間ステップで得られた解をＩＣＣＧ法等による反復解法における新たな初期値として用いても、大幅な解析の高速化を図ることはできない。

本発明は、磁界解析の高速化を図ることを目的とする。

本発明の一つの特徴によれば、コイル電流等の磁界発生ソースの変動量に伴う変動場を解くことにより磁界解析を行う。

本発明のその他の特徴は、以下の発明を実施するための最良の形態の欄で詳細に説明する。

本発明に関する磁界解析法によれば、磁性材料を含む非線形磁界を高速に解析できるという効果がある。

以下、本発明の実施例を説明する。

本発明に関る磁界解析法では、磁気ベクトルポテンシャルＡならびに電気スカラーポテンシャルφそのものを解くのではなく、それらの変化量を逐次解いていく。これに関する基礎方程式は、まだ一般的に周知されていないため、ここでＡ−φ法をベースに導いておく。コイル電流密度ベクトルの微小変化量をΔＪ₀として、磁界ベクトルの変化量ΔＨ，磁化ベクトルの変化量ΔＭに対して、

なる関係式における微分磁化率χ_d を導入する。磁束密度Ｂ，磁界Ｈ，磁化Ｍの関係式として、

より、各項の微小変化量に関する式は、

となる。（数４）を用いると、

となる。数５にベクトル微分演算子である回転（rotation）を施して、

とし、これにFaradayの式

を代入すると、

を得る。この式の変化分をとると、

これに、数７を代入して式を整理すると、

を得る。ここに、μ_dは微分透磁率で、微分磁化率χ_dと

の関係にある。微分透磁率とは、磁性材料の磁界Ｈと磁束密度Ｂの関係式であるＢ＝
Ｂ(Ｈ)において、

で定義される物理量であり、μ＝Ｂ／Ｈで定義される通常の透磁率とは異なる。

なお、ここで、微分透磁率μ_d は、求めようとしている時刻ｔ＝ｎΔｔにおける値ではなく、前時刻ｔ＝(ｎ−１)Δｔとの中間の時刻であるｔ＝(ｎ−１／２)Δｔあるいはその近傍における値である点に注意しなければならない。通常の非線形磁界解析において、透磁率μが時刻ｔ＝ｎΔｔでの値になるようにNewton−Raphson 法でセルフコンシステントに解を求める場合と事情が異なる。

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

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

となる。ここに、

は、前の時刻におけるΔＡの時間に関する偏微分値である。

これに加えて、数１６にベクトル微分演算子である発散（divergence）を施した式

が基本式となる。

数１６，数１８を用いて磁気ベクトルポテンシャルの変化量ΔＡ並び電気スカラーポテンシャルの変化量△φを求め、磁束密度変化量ΔＢを更新する。新たなΔＢに対して、微分透磁率μ_d を更新し、再度、数１６，数１８を解く。このような操作を、ΔＢが収束するまで繰り返す。ΔＢが収束後、時刻ｔ＝ｎΔｔでの解析はこれで終了する。

なお、Ａ法における基礎方程式は、

であり、数１６に相当する時間軸方向に関する差分式は、

である。

数１６，数１８あるいは数２０を例えば有限要素法で離散化すると、行列方程式

を得る。ここで、左辺のＰは離散化して得られた行列で、未知変数ベクトルｘは、未知変数△Ａ，△φをひとつのベクトル形式にまとめたものである。また、右辺のｂは数１６，数２０の右辺の項で形成されたソース項を表すベクトルである。この行列方程式を反復法で解く場合に、未知変数ベクトルｘの初期値として、前の時刻で得られた解ｘ₀ を選ぶと、解の収束が速まる。さらに、前時間ステップで求めた解ｘ₀ そのものを用いるよりも、それに次式で示す修正因子ｇ

ここに、ｂ：前時間ステップの行列方程式の右辺ベクトル
ｂ′：現時間ステップの行列方程式の右辺ベクトル
をかけたｇｘ₀ を初期値にすると、求める解にさらに近づくため、ＩＣＣＧ法における反復回数は１／１０程度に激減し、計算時間は、１／４程度に短くなる。ただし、この単一の修正因子ｇによる高速化は、電流変動する全てのコイルが相対的に同じ電流変化を示すときにのみ有効である。電流変動する全てのコイルに異なる電流波形の電流が流れる場合は、それぞれのコイルに対して個別の修正因子ｇを用意し、場の変化量に対するそれぞれの寄与を個別に求め、それらの線形和を初期値とすれば、解析は高速化できる。具体的に式で表現すると、コイルｋ（ｋ＝１，２，…Ｎ）を流れる電流変動によるソース変動量をｂ_kとおいて、Ｐχ_k＝ｂ_kなる解ｘ_kを上と同様に求め、その総和であるΣｘ_k を初期値にして、数２１をＩＣＣＧ法等の反復解法で求めれば、解ｘを高速に求めることができる。

以下、図を用いて本発明の第一実施例における磁界解析法の流れを説明する。用いる磁界解析法として、一般的な有限要素法をとりあげる。

図１に解析システムの一例の図を示す。本解析システムは、計算機１，表示装置２，記憶媒体３から構成される。ここでは、記憶媒体３を明示するために、計算機１の外に出しているが、計算機１内部に記憶媒体３を設置しても良い。

図２は本発明の一実施形態として、数１６，数１８あるいは数２０を基本解析式とする磁界解析法の流れを示す。時間刻み幅をΔｔとして、図２は第ｎ番目の時刻ｔ＝ｎΔｔにおける磁界解析の流れを示している。以下、図２を用いて、順を追って解析の流れを説明する。

第１ステップ１１において、初期磁化曲線から求めた微分透磁率μ_d を磁性材料の各要素に割り当て、行列方程式数２１に関する行列Ｐを作成する。ここで、数５で定義されている微分透磁率μ_dは、数値計算上は

で計算することになるが、ここで、ΔＨあるいはΔＢは、適当な微小量をとり、用意した初期磁化曲線からΔＨを与えた場合はΔＢを、ΔＢを与えた場合はΔＨを求める。こうして、数２３より仮の微分透磁率μ_d を求める。ここで、作成する行列Ｐは、すべて作成しなおすのではなく、前時間ステップｔ＝(ｎ−１)Δｔにおいて作成した行列に対して、磁性材料の要素に関連する成分のみについてのみ修正すればよい。次に、第２ステップ１２において、ソース項の変化量に相当する行列方程式数２１の右辺ベクトルｂを作成する。第３ステップ１３において、ＩＣＣＧ法等の反復解法を用いて、行列方程式数２１の解を求める。この際、反復解法の初期値として、前時間ステップの解ｘ₀ あるいはそれに数
２２に示したｇをかけたｇｘ₀ を用いると、高速に収束解が求まる。第４ステップ１４において、第３ステップで求めた解から、磁束密度変化量ΔＢを求める。第一回目は、第５ステップ１５を直接実行する。すなわち、磁束密度変化量ΔＢと初期磁化曲線からΔＨを求め、数２３より微分透磁率μ_d を更新して、再度、第１ステップ１１に戻り、第４ステップ１４までの計算を実行する。求めた磁束密度変化量ΔＢが前の計算で求めたものに基準範囲内になければ、すなわちΔＢが収束していなければ、第５ステップ１５において、磁束密度変化量ΔＢと初期磁化曲線から、微分透磁率μ_d を更新して、第１ステップ１１に戻り、第４ステップ１４までの計算を実行する。磁束密度変化量ΔＢが収束すれば、第６ステップ１６において、ｔ＝ｎΔｔにおける磁性体における解ｘおよび数２１の右辺ベクトルｂを計算機の記憶媒体３に記憶して、次の時刻ｔ＝(ｎ−１)Δｔにおける解析に移行する。

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

となる。ここに、μは透磁率、σは導電率、Ｊ₀はコイルの電流密度ベクトル、Ｍ₀は磁石の残留磁束密度ベクトルである。式（１）の右辺の項が場を励磁するためのソース項である。

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

そこで、数２４を時刻ｔ＝ｎΔｔ（ｎ＝１，２，３，…）における方程式とみなし、次の時刻ｔ＝(ｎ＋１)Δｔにおける方程式として、

を考える。ここに、

で、ΔＡ，Δφ，Δμは、それぞれ時間がΔｔ経過したことに伴うＡ，φ，μの微小変化量を表す。

まず、近似的にμ′＝μとおいて、線形磁界解析を実行する。時刻ｔ＝(ｎ＋１)Δｔにおける磁束密度の変化量ΔＢ_n+1をΔＢ_n+1＝rot(ΔＡ) により求め、要素分割した各要素に対して、内積ΔＢ_n・ΔＢ_n+1の正負によって、上昇曲線か下降曲線かを判別する。具体的には、ΔＢ_n・ΔＢ_n+1＞０のとき分岐点を形成せずにそのままヒステリシス曲線をたどり、ΔＢ_n・ΔＢ_n+1＜０のとき分岐点を形成し、ヒステリシス曲線上を反転して別の曲線上をたどる。言い換えると、上昇曲線上に動作点があった場合、次の時点でΔＢ_n・
ΔＢ_n+1＞０のとき、そのまま上昇曲線上をたどり、ΔＢ_n・ΔＢ_n+1＜０ のとき、下降曲線に移行する。下降曲線上に動作点があった場合、次の時点でΔＢ_n・ΔＢ_n+1＞０のとき、そのまま下降曲線上をたどり、ΔＢ_n・ΔＢ_n+1＜０のとき、上昇曲線に移行する。

上昇曲線か下降曲線かの判別のあと、非線形磁界解析に移行する。ヒステリシス解析において、上昇曲線あるいは下降曲線をそのままたどっていく場合には、通常のNewton−
Raphson 法で問題なく解ける。しかし、分岐点を形成してヒステリシス軌跡が反転した場合、取り扱いが厄介になる。分岐点の前後でＢＨカーブの勾配が特異的に変化するため、Newton−Raphson 法をそのまま適用するには問題がある。ＢＨカーブの勾配が特異的に変化しないように分岐点を越えた反対側の領域までなめらかな仮想的なＢＨカーブを形成しておけば問題ないが、その方法については、今後比較検討することにして、ここでは取りあえず、Newton−Raphson 法を使わずに単純反復法で解を求めることにした。その代わり、磁気ベクトルポテンシャルそのものに関する数２６を直接解くのではなく、ソースの変化量ΔＪ_n に対する場の変化量ΔＡを求めることにし、その結果として、Ａ＋ΔＡを求めることにした。

磁界ベクトルの変化量ΔＨ，磁化ベクトルの変化量ΔＭに対して、

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

ここで、磁束密度Ｂの変化量をより正確に求めるための基礎式を導いておく。まず、磁化をＭとして、

より、

となる。数２７より、

となるので、再び数２７より、

となる。数２８にrotationを施して、

とし、これにFaradayの式

を代入すると、

を得る。この式の変化分をとると、

これに、数３１を代入して式を整理すると、

を得る。ここに、μは微分透磁率で、微分磁化率χと

の関係にある。数３６は数２４の変化分に関する式に相当するが、透磁率μが微分透磁率μに置き換わることに注意しなければならない。

数３６は、ヒステリシスのない非線形磁界解析の基本式である数２４に対して、透磁率μを微分透磁率μに置き換え、ソース項のコイル電流密度Ｊ₀を変化量ΔＪ₀に置き換えた形になっており、既存の磁界解析プログラムにヒステリシス解析機能を容易に追加することができる。

なお、ここで、微分透磁率μや微分磁化率χは、求めようとしている時刻ｔ＝
(ｎ＋１)Δｔにおける値ではなく、時刻ｔ＝(ｎ＋１／２)Δｔ近傍における値である点に注意しなければならない。ヒステリシスのない通常の非線形磁界解析において、透磁率μが時刻ｔ＝(ｎ＋１)Δｔでの値になるようにNewton−Raphson 法でセルフコンシステントに解を求める場合と事情が異なる。数２７に示した磁界の変化量ΔＨと磁化の変化量ΔＭの関係式に示した勾配μ₀χ は、図３に示す三角形の斜辺の勾配に相当する。このため、微分透磁率や微分磁化率という表現よりも、より数値解析の特性を反映した「差分透磁率」「差分磁化率」という表現の方が誤解を与えなくてよい。しかし、このような名称は一般的ではないため、あえてここでは、微分透磁率や微分磁化率という表現のままにしておく。

上昇曲線か下降曲線かを判別する際に求めたΔＢは、数において近似的にμ′＝μとおいて線形磁界解析で求めたもので解析精度が低いので、改めて、数３６を解くことにより、ΔＢ＝rot(ΔＡ) をより精度の高いものに更新できる。このようにして、順次、微分透磁率μとΔＢを反復更新することにより、時刻ｔ＝(ｎ＋１)Δｔにおける（Ｈ，Ｂ）空間での動作点を各要素について求めることができる。

磁気ベクトルポテンシャルの変化量ΔＡに関する基本数３６において、時間軸方向に差分をとる。Newmarkのβ法を用いると、

となる。ここに、βは０＜β＜１なる固定定数であり、安定した時間応答解析のために１／２＜β＜１なるβを用いる。

ここで、ΔＡ(ｔ)＝０であり、ΔＡ(ｔ＋Δｔ)はこれまでの式の展開と同様に単にΔＡとおき、∂ΔＡ(ｔ)／∂ｔをΔＡ₀と略記して、数３６に代入すると、

となる。数３９を用いて磁気ベクトルポテンシャルの変化量ΔＡを求め、磁束密度変化量ΔＢを更新する。新たなΔＢに対して、微分透磁率μを更新し、再度、数３９を解く。

このような操作を、ΔＢが収束するまで繰り返す。ΔＢが収束後、透磁率分布が求まり、時刻ｔ＝(ｎ＋１)Δｔでの解析は、基本的にはこれで終了する。ただし、毎回、場の変化量のみを求めるとなると、数値誤差が積算する恐れがあるため、求めた透磁率分布μ′＝μ＋ΔμとＡ′＝Ａ＋ΔＡが数２５を満足しているかどうかをチェックし、必要ならば、Ａ′＝Ａ＋ΔＡならびにμ′＝μ＋Δμを補正する。解は真の解の近傍にあるため、この補正のための計算は短時間で済む。

以上の解析の流れをフローチャートにして示すと、図４のようになる。

本磁界解析では、Ａ−φ法を用いているので、図４にも示したように、磁束密度変化量ΔＢを求めた後、ＢＨ関数を用いて、磁界変化量ΔＨを求め、微分透磁率（差分透磁率）μ＝ΔＢ／ΔＨを導出しなければならない。このため、新たに導入するヒステリシスモデルでは、磁界Ｈを磁束密度Ｂの関数として表現しなければならない。

このヒステリシスモデルは、まず、Ｈの関数形としてのＢを次のような形に表す。

ここに、ｆ(Ｈ)はヒステリシス軌跡を記述するＨに関する関数である。また、Ｈ_i，Ｂ_iはヒステリシス曲線の出発点に相当する分岐点で、Ｈ_f，Ｂ_fはヒステリシス曲線が目指す最終点に相当する分岐点である。

これらの分岐点を図５に例示しておく。図５には、出発点と最終点をなす分岐点をそれぞれ３個示してある（添字０，１，２で区別）。ヒステリシス軌跡上の矢印で示すように、いずれも上昇曲線側のものを示しているが、下降曲線側の場合は、それぞれ、出発点と最終点をなす分岐点の位置が互いに入れ換わる。数４０を変形すると、

を得る。ｆ^-1(…)は、関数ｆ(Ｈ)の逆関数であり、関数ｆ(Ｈ)は１価関数であるため、数値的に容易に求めることができる。あるいは、見方を変えて、ＨをＢの関数形として直接次のような形に表してもよい。

数４１を用いた本ヒステリシスモデルによるヒステリシス軌跡の計算例を図６〜図９に示す。図６，図８に示した磁束密度の時間変化に対するヒステリシス軌跡を、それぞれ図７，図９に示す。これらの図で、磁束密度Ｂおよび磁界Ｈは、それぞれ最大値Ｂｍａｘ，Ｈｍａｘで規格化した値で表示してある。図６，図７の例は、振動しながら次第に磁束密度が増大する場合の例であり、図８，図９の例は、逆に振動しながら次第に磁束密度が減少する場合で、いわゆる消磁の例である。ここで、図７，図９では、メジャーループは破線で示してあり、実線はヒステリシス軌跡を示しており、図中の矢印は、それに沿ってヒステリシスの軌跡を描くことを示している。

本発明の実施形態を実現する解析システムを示す図。 本発明の一実施例である磁界解析法の流れを示す図。 （Ｈ，Ｍ）空間における変化分。 ヒステリシス解析フローチャート。 ヒステリシス軌跡における分岐点。 磁束密度の時間変化（例１）。 （Ｂ，Ｈ）空間におけるヒステリシス解析結果（例１）。 磁束密度の時間変化（例２）。 （Ｂ，Ｈ）空間におけるヒステリシス解析結果（例２）。

符号の説明

１…第１および第２の実施形態における計算機、２…第１および第２の実施形態における表示装置、３…第１および第２の実施形態における記憶媒体、１１…本実施例におけるステップ１、１２…本実施例におけるステップ２、１３…本実施例におけるステップ３、１４…本実施例におけるステップ４、１５…本実施例におけるステップ５、１６…本実施例におけるステップ６。

Claims (8)

1. コイル電流等の磁界発生ソースの変動量に伴う変動場を解くことを特徴とする磁界解析法。
2. 微分透磁率と変動場を交互に修正しながら計算して収束解を得ることを特徴とする磁界解析法。
3. 微分透磁率を用いた変動場に関する方程式を有限要素法で解くことを特徴とする磁界解析法。
4. 請求項第１から３項において、行列方程式を反復解法で解く場合に、変動場を表現する未知数ベクトルの初期値として、前時間ステップでの解を用いることを特徴とする磁界解析法。
5. 請求項第１から３項において、行列方程式を反復解法で解く場合に、前時間ステップでのコイル電流等の磁界発生ソースの変動量に関する前記行列方程式の右辺ベクトルをｂとし、現時間ステップでのコイル電流等の磁界発生ソースの変動量に関する前記行列方程式の右辺ベクトルをｂ′として、変動場を表現する未知数ベクトルの初期値として、前時間ステップでの解に因子ｇ＝ｂ・ｂ′／｜ｂ｜｜ｂ′｜をかけたものを用いることを特徴とする磁界解析法。
6. コンピュータに請求項第１項から請求項第５項までのいずれかの処理を実行させることを特徴とするプログラム。
7. 請求項第６項記載のプログラムを記録したコンピュータ読み取り可能な記録媒体。
8. メッシュ数を変化させて磁界解析を行うことを特徴とする磁界解析法。
   

Images