ポスト処理による鉄損算出

電気機器の損失の主成分である積層鉄心に発生する鉄損を電磁気解析で算定する方法として、現在様々な手法が考案されています。簡便な方法として、計算で得られた磁束密度よりポスト処理として算出する方法があります。この度EMSolutionでは、ポスト処理として近似的に鉄損算出する機能を追加しました。鉄損計算式は、文献（1）より、磁束密度の絶対値の最大値を用いる方法と、磁束密度波形から直接求める方法を採用しております。下記式中では直交座標系（x, y, z）としていますが、円筒座標系（r, θ, z）とすれば回転磁界も考慮できます。
文献（1）に基づき、鉄損W_iは（1）式に示されるように、渦電流損W_e、ヒステリシス損W_hに分離して算出します。さらに渦電流損W_eは、電磁鋼板の厚み方向（積層方向）に還流するW_e^と、鋼板面（積層面内方向）を流れる鋼板面（積層面内方向）を流れるW_e||に分けられます。

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　（1）

積層方向の渦電流損W_e^、ヒステリシス損W_hはポスト処理としてFEM解析により得られた磁束密度の履歴より計算されます。積層面内の渦電流損W_e||は、FEM解析より得られた渦電流J_eを用いて直接得られます。

算出法�@：磁束密度の最大値から求める方法
単一周波数f、振幅B_max（磁束密度の最大値）の交番磁界が印加された場合、鉄損は下記の鉄損推定式から計算されます。

　　　　　　　　　　（2）

ここで、DV_iは有限要素iの体積（m³）を表し、K_eは渦電流損係数（W/kg/T²/Hz²）、K_hはヒステリシス損係数（W/kg/T²/Hz）を表します。また、a、b、gは累乗を表す係数で、デフォルトではa = b = g = 2と設定されます。磁束密度波形が正弦波形を仮定しているため、振幅を磁束密度の最大値の絶対値として算出します。
積層面内の渦電流損W_e||は、FEM解析より次式で直接計算されます。

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　（3）
、 　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　（4）

ここで、Tは周期（s）、aは占積率、s_sは鉄部の導電率（S/m）を表します。これによりFEM解析では積層面内方向の渦電流のみを占積率を考慮して計算することになります。

算出法�A：磁束密度波形から直接求める方法
磁束密度波形が正弦波で無い場合、またマイナーループを含む場合には、（2）式の近似精度は高くないため、ここで示すような磁束密度波形から直接求める方法が有効な方法となります。このときの積層方向の渦電流損W_e^は、以下の式で計算されます。

　　　　　　　　　　　（5）

ここで、Dは鉄心の密度（kg/m³）、Nは一周期あたりの反復（サンプリング数）回数、Dtは時間刻み（s）、Bの上添え字のkは時間ステップを表します。（5）式で表される渦電流損には、積層方向である磁束密度のz成分に関する項が除かれています。これは、積層方向の磁場により誘起される積層面内を流れる渦電流損を表しており、FEM解析より（3）式で求められるためです。
また、ヒステリシス損W_hは以下の式で計算されます。

　　　　　　　　　　　　　　　（6）

ここで、NEは鉄心領域中の要素数、N_pxⁱ、N_pyⁱ、N_pzⁱはi 番目の有限要素における磁束密度の各成分の時間変化に対する極大・極小値の個数、B_mx^ij、B_my^ij、B_mz^ijはメジャーおよびマイナーループの振幅であり、Fig.1に示すような方法で、磁束密度波形の全極大・極小値より決定しています。その際、あるループにおける磁束密度の振幅の2乗が2回加算されるため、加算後に2で割っています。これは、メジャーループとマイナーループのヒステリシス損を同一のアルゴリズムで求めています。

なお、鉄損係数K_e、K_hは、計測された各周波数f、磁束密度Bにおける鉄損W_iから、代表磁束密度（例えばB_m=1T）でのf（Hz）-W_i/f（W/kg/Hz）グラフを作成し、それを線形近似した傾きをK_e、切片をK_hとして使用します。
また、積層鉄心を均質化法（PACKING）を適用して解析し、鉄損を算出する場合、「均質化法による積層鉄心の解析と積層鉄心鉄部の磁束密度の出力」で説明した積層鉄心鉄部の磁束密度を使用します。

次に具体例として、文献（1）の「鉄損計算法検証用モデルの解析および実験」の直流ブラシレスモータを使用して鉄損を算出したものを示します（Fig.2）。このモデルは、ステータ、ロータ共に無方向性電磁鋼板50A1300から成る積層鉄心で、コイル巻き線なしで強制的にロータを回転させたときのロータの永久磁石による鉄損を測定し、解析により得られた結果と比較、検討するための検証用モデルです。鉄損計算には、円筒座標の磁束密度成分を用いるため、B_rの一周期である機械角180度計算した結果を用いました。積層鉄心を導電率ゼロの鉄心として扱った二次元、三次元非線形静磁場解析結果と、積層鉄心を均質化法（PACKING）により近似した渦電流を含む非線形過渡解析結果を示します。ここで、BHカーブは、L+C方向（圧延方向とその垂直方向の平均値）を使用し、積層方向はz方向です。ただし、占積率は文献中に明記されていないため、98.5%と仮定しました。また、ロータの回転数は1500r/minで比較しました。
Fig.3に、ポスト処理により得られた三次元静磁場解析結果の積層方向の渦電流損とヒステリシス損分布を示します。Fig.4に、それぞれの鉄損算出結果を実験値と比較して示します。Table Iにはその鉄損の内訳を示します。算出法�@、�Aの結果と比べると、算出法�@の方が値は小さくなっています。これは、磁束密度波形の歪みとマイナーループよるものと思われます。また、Fig.5に均質化法を適用した渦電流を含む三次元非線形過渡解析結果である面内方向の渦電流損分布を示します。値としては小さいですが上側の端部に渦電流が集中しています。
なお、算出法�@では印加磁場は交番磁界を想定しています。そのため、ロータの磁束密度は位置に対して時間変化していないにも関わらず、磁束密度の絶対値の最大値を使用しているため、実際には発生していない鉄損を算出してしまいます。そのため、ステータ部のみの値を用いて比較しています。積層鉄心に交番磁界を印加する簡単なモデルで解析した例を、「均質化法による積層鉄心の解析と積層鉄心鉄部の磁束密度の出力」で示しています。そこでは算出法�@、�A共に良い一致を示しています。

今回使用した検証モデルの実機では、無方向性電磁鋼板50A1300をロータ・ステータ共に回積みして作成しています。さらに接着剤により積層を作成し、磁気特性を劣化させる要因である応力や焼きばめの影響を除いています。そのため、実機の測定結果と解析結果が非常に良く合うものと考えられます。ただし、上記の鉄損推定法�@、�A共に、特にヒステリシス損は、磁束密度のべき乗に比例するというスタインメッツの式を基に計算するものです。実際には、ヒステリシス損はヒステリシスループの面積より計算されるもので、それを考慮するような他の提案手法^（2）と比べると精度的に劣ると思われますが、ご承知の上お使い頂ければと思います。

＜参考文献＞
（1）　回転機の三次元電磁界解析高度化調査専門委員会：「回転機の電磁界解析高度化技術」、電気学会技術報告、第942号 （2004）
（2）　回転機の電磁界解析高精度モデリング技術調査専門委員会：「回転機の電磁界解析高精度モデリング技術」、電気学会技術報告、第1044号（2006）

鉄損算出はポスト処理として計算され、多数ステップ計算した非線形静磁場解析、過渡解析で行うことができます。
�@　ポスト処理として計算するので、メイン計算実行後、PRE_PROCESSING, MAKING_MATRICES, SOLVING_EQUATIONを0とし、POST_PROCESSINGのみ1と設定します。


�A　一周期ごとの鉄損を計算するために、EMSolution Handbook「9. 出力ステップ、フェーズ」のSTEP_INTERVALを一周期分のステップ数に設定します。


�B　「10. 入出力ファイル」のAVERAGEを1とします。これは通常、周期ごとの平均熱発熱量を計算するために使用されますが、鉄損計算の場合も周期ごとに出力するために使用します。


�C　鉄損分布の出力ファイルとして、iron_lossファイルが出力されます。「10. 入出力ファイル」のMAGNETIZATIONの次に新たにできたIRON_LOSSに、1（W/kg）または2（W/m³）を設定します。また、面内方向の渦電流損である発熱を出力したい場合は、HEAT=1と設定してください。


�D　一周期分の鉄損をoutputファイルに出力するために、「11. プリント出力：11.1 出力オプション」の3カラム目のMAGNETIC_ENERGYの後に新たにできたIRON_LOSS を1（磁束密度の最大値の絶対値による算出法�@）、または2（磁束密度波形による算出法�A）と設定します。また、HEAT=1と設定すれば、面内方向の渦電流損も出力することができます。


�E　算出法�@であるIRON_LOSS=1と設定した場合、次のカラムに周波数FREQUENCY（Hz）、係数であるALPHA、BETA、GAMMAを入力します。ALPHA、BETA、GAMMAは省略することもでき、その場合デフォルトとして全て2.0と設定されます。


�F　鉄損係数を設定するために、「16.1.1 体積要素特性」のANISOTROPYの後に新たにできたIRON_LOSS に1を設定します。鉄損計算に必要なパラメータとして、磁束密度の成分を評価する局所座標COORD_ID、密度MASS_DENSITY（kg/m³）、渦電流損係数KE（W/kg/T²/Hz²）、ヒステリシス損係数KH（W/kg/T²/Hz）を設定します。また、通常の計算で入力したPACKINGと導電率の異方性の入力例も合わせて示します。 ただし、PACKINGを使用して全体座標系で計算した結果を、円筒座標系にて鉄損算出および磁束密度を出力したい場合は、PACKINGのCOORD_IDには新たにz方向が積層方向となるように定義した局所円筒座標IDを設定し直す必要があります。


PACKINGを使用しない場合は、鉄損算出パラメータのCOORD_IDに、積層方向をz方向として定義した局所円筒座標IDを指定します。


�G　EMSolutionを実行すると、outputファイル中にプロパティごとに、積層方向の渦電流損、ヒステリシス損が出力されます。ポストファイルとしてiron_lossに積層方向の渦電流損、ヒステリシス損分布が出力されます。HEATオプションを使用した場合は、同じようにoutputファイルにプロパティごとの面内方向の渦電流損が、ポストファイルheatにその分布が出力されます。

＜出力例：三次元静磁場解析の鉄損算出結果＞

＜出力ファイル内容：FEMAP Neutral file形式の場合＞
iron_loss:
LOSS-elem-1　異常渦電流損 （=1:W/kg or =2:W/m³）
LOSS-elem-2　異常渦電流損 （W）
LOSS-elem-3　ヒステリシス損 （=1:W/kg or =2:W/m³）
LOSS-elem-4　ヒステリシス損 （W）

使用データ：input3D_static：三次元静磁場解析用
　　　　　 　　input3D_transient：三次元渦電流解析（PACKING）用
　　　　　 　　inputPost_ironLoss_Bdiff：算出法�A鉄損算出用。算出法�@に変更するには、IRON_LOSS=1とし、FREQUENCYのコメントアウトを外してください。
　＊ 例題で使用したメッシュファイルは形状秘匿のため公開することができません。ご了承下さい。
メッシュデータも含めた簡単なモデルは「均質化法による積層鉄心の解析と積層鉄心鉄部の磁束密度の出力」にありますので、合わせてご参照ください。