(57)【要約】【課題】Ｃａｕｅｒ回路と有限要素法との連成解析による電磁界シミュレータの計算コストを低減する。【解決手段】解析モデル２００は、有限要素法によって複数の要素３１１，４１１に分割される。各要素に標準型Ｃａｕｅｒ回路が等価的に配置されるＣａｕｅｒ回路と有限要素法との連成解析のために、有限要素法に対応する第１の変数と、Ｃａｕｅｒ回路に対応する第２の変数とを含む連立一次方程式の定常解が、時刻ステップ毎の繰り返し計算によって求められる。【選択図】図６
2020154560.tif 000002
【特許請求の範囲】
【請求項１】
磁性素子を含む解析モデルを対象とする電磁界シミュレータであって、
有限要素法によって分割された前記解析モデルの各要素に対応して配置されたＣａｕｅｒ回路を連成することによって得られた、前記有限要素法に対応する第１の変数と、前記Ｃａｕｅｒ回路に対応する第２の変数とを含む第１の連立一次方程式について、時刻ステップ毎の繰り返し計算によって定常解を求めるための手段を備えた、電磁界シミュレータ。
【請求項２】
前記定常解を求めるための手段は、
各前記時刻ステップにおいて、前記第１の連立一次方程式を縮約することによって予め得られた、前記第２の変数を含まない第２の連立一次方程式から当該時刻ステップにおける前記第１の変数を導出する第１の演算を実行する手段と、
前記第１の演算によって導出された前記第１の変数が代入された前記第１の連立一次方程式から前記第２の変数を導出する第２の演算を実行する手段とを含む、請求項１記載の電磁界シミュレータ。
【請求項３】
前記Ｃａｕｅｒ回路を構成する複数の抵抗素子およびインダクタンス素子のうちのＣａｕｅｒ回路の初段のインダクタンス素子のインダクタンスは、プレイモデルに従って算出される、請求項１又は２に記載の電磁界シミュレータ。
【請求項４】
磁性素子を含む解析モデルを対象とする電磁界解析方法であって、
有限要素法によって分割された前記解析モデルの各要素に対応して配置されたＣａｕｅｒ回路を連成することによって得られた、前記有限要素法に対応する第１の変数と、前記Ｃａｕｅｒ回路に対応する第２の変数とに係る第１の連立一次方程式について、時刻ステップ毎の繰り返し計算によって定常解を求めるためのステップを備える、電磁界解析方法。
【請求項５】
前記定常解を求めるためのステップは、
各前記時刻ステップにおいて、前記第１の連立一次方程式を縮約することによって予め得られた、前記第２の変数を含まない第２の連立一次方程式から当該時刻ステップにおける前記第１の変数を導出する第１の演算を実行するステップと、
前記第１の演算によって導出された前記第１の変数が代入された前記第１の連立一次方程式から前記第２の変数を導出する第２の演算を実行するステップとを含む、請求項４記載の電磁界解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、電磁界シミュレータおよび電磁界解析方法に関し、より特定的には、磁性素子の周波数特性を解析するための電磁界シミュレータおよび電磁界解析方法に関する。
【背景技術】
【０００２】
磁性素子は、磁性材料からなるコアと、コアに励磁電流を供給するコイルとを有し、励磁されたコアは磁化曲線上で動作する。磁性素子の設計において、コンピュータを用いた電磁界解析が、研究開発や設計業務で広く用いられている。しかしながら、近年の昇圧回路等の大電流化に伴い、コアの磁気飽和を緩和するため磁路に磁気ギャップを設ける等、構造が複雑化することで、解析値と実測値とが大きく乖離する問題が生じている。また、商用周波数の正弦波などにリプルが重畳した電流波形など複数の周波数を含む複雑な電流波形で動作する磁性素子において、周波数の影響による磁気特性の変化を考慮した検討を行う場合には、任意の周波数で測定した磁気特性が必要になる。
【０００３】
例えば、国際公開第２０１０／０３８７９９号（特許文献１）には、直流重畳特性を解析するために、電流の増加に伴って変化する比透磁率についてマイナーループを考慮して計算する解析方法が記載されている。しかしながら、特許文献１の解析方法では、初磁化特性の他に大量のマイナーループ群をリプル幅に応じて測定する必要があるため、測定データが膨大になることが懸念される。
【０００４】
他の公知の解析方法として、プレイモデルと呼ばれる一定磁束密度ごとに測定した磁化特性のループ群を用いる手法が提案されている。しかしながら、複数の周波数もしくは周波数の変化に伴うインダクタンス値の変化を計算する手法は実用化されていないため、プレイモデルを用いて周波数の影響などを考慮する場合には、周波数毎のプレイモデルを測定する必要があるため、膨大なデータを測定することが必要となる。
【０００５】
あるいは、非特許文献１には、圧粉磁芯への適用を目的とした、Ｃａｕｅｒ回路を用いた磁気特性モデル化手法が提案されている。圧粉磁芯は、数十μｍ程度の微小な磁性粒子を絶縁被膜で覆い、圧縮成形することで作製される軟磁性材料であり、電磁鋼板に比べて渦電流損が低いことに加えて，加工が容易であるといった特長を有している。
【０００６】
圧粉磁芯の渦電流を考慮した磁界解析を行う際に、個々の磁性粒子の形状および大きさが同一であり、かつ、粒子同士が絶縁されていることを仮定した、いわゆる均質化に基づく解析が提案されている。しかしながら、実際の圧粉磁芯は磁性粒子では、形状および大きさは均一ではなく、成形時の圧力により絶縁の一部が破壊され、粒子同士が不均一に接触していることも少なくないため、粒子内部に流れる局所的な渦電流と、粒子を横断して流れる大域的な渦電流が存在する。このような理由から、単純な均質化法による広範な周波数範囲での磁気特性のモデル化には限界が存在する。
【０００７】
これに対して、非特許文献１では、直流ヒステリシス特性をプレイモデルで表現するとともに、渦電流による周波数依存性を標準Ｃａｕｅｒ回路（Ｃａｕｅｒ−Ｉ型回路）によって考慮することにより、広範な周波数範囲で高精度に測定結果を再現することが可能である。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】国際公開第２０１０／０３８７９９号
【非特許文献】
【０００９】
【非特許文献１】渡邉直也、高橋康人、藤原耕二、堺香代、島津英一郎「Ｃａｕｅｒ回路を用いた圧粉磁芯の磁気特性モデリングに関する基礎的検討」、電気学会静止器・回転機合同研究会資料、SA-17-65，RM-17-96，pp. 73-78（2017）
【発明の概要】
【発明が解決しようとする課題】
【００１０】
しかしながら、非特許文献１による磁気特性モデリングでは、渦電流による周波数依存性を考慮した標準Ｃａｕｅｒ回路を有限要素法に組み込むことで、Ｃａｕｅｒ回路と有限要素法との連成解析が実行されるため、通常の有限要素法で定義される未知変数に加えて、圧粉磁芯要素内でＣａｕｅｒ回路に対応する未知変数が増加するため、計算コストが大きく増加してしまうことが懸念される。
【００１１】
本発明は、このような問題点を解決するためになされたものであって、本発明の目的は、Ｃａｕｅｒ回路と有限要素法との連成解析による電磁界シミュレータおよび電磁界解析方法の計算コストを低減することである。
【課題を解決するための手段】
【００１２】
本発明のある局面では、磁性素子を含む解析モデルを対象とする電磁界シミュレータであって、有限要素法によって分割された解析モデルの各要素に対応して配置されたＣａｕｅｒ回路を連成することによって得られた、有限要素法に対応する第１の変数と、Ｃａｕｅｒ回路に対応する第２の変数とを含む第１の連立一次方程式について、時刻ステップ毎の繰り返し計算によって定常解を求めるための手段を備える。
【００１３】
本発明の他のある局面では、磁性素子を含む解析モデルを対象とする電磁界解析方法であって、有限要素法によって分割された解析モデルの各要素に対応して配置されたＣａｕｅｒ回路を連成することによって得られた、有限要素法に対応する第１の変数と、Ｃａｕｅｒ回路に対応する第２の変数とに係る第１の連立一次方程式について、時刻ステップ毎の繰り返し計算によって定常解を求めるためのステップを備える。
【発明の効果】
【００１４】
本発明によれば、Ｃａｕｅｒ回路と有限要素法との連成解析による電磁界シミュレータおよび電磁界解析方法の計算コストを低減することができる。
【図面の簡単な説明】
【００１５】
【図１】磁気特性測定の対象とした圧粉磁芯コアの諸元の一覧を示す図表である。
【図２】圧粉磁芯コアの周波数５０（Ｈｚ）でのヒステリシスループの測定結果を示すグラフである。
【図３】一定の最大磁束密度の下で周波数を変化させたときの圧粉磁芯のヒステリシスループの測定結果を示すグラフである。
【図４】図１〜図３で説明した圧粉磁芯にコイルを巻回して構成されたインダクタのインダクタンスの周波数特性の測定結果を示すグラフである。
【図５】本実施の形態に係る電磁界解析で用いられるＣａｕｅｒ−Ｉ型回路（標準Ｃａｕｅｒ回路）の回路図である。
【図６】本実施の形態に係る電磁界シミュレータの解析対象および解析モデルの第１の例を説明する概念図である。
【図７】本実施の形態に係る電磁界シミュレータによる電磁界解析の演算処理を説明するフローチャートである。
【図８】本実施の形態に係る電磁界シミュレータを実現するためのハードウェア構成の一例を説明するブロック図である。
【図９】ヒステリシスループの測定結果と本実施の形態に係るシミュレーション結果とを比較するグラフである。
【図１０】インダクタンスの周波数特性の測定結果と本実施の形態に係るシミュレーション結果とを比較するグラフである。
【図１１】本実施の形態に係る電磁界解析の所要時間を他の手法の電磁界解析と比較する図表である。
【図１２】本実施の形態に係る電磁界解析における時刻ステップあたりの収束計算結果を他の手法の電磁界解析と比較する図表である。
【図１３】本実施の形態に係る電磁界シミュレータの解析対象および解析モデルの第２の例を説明する概念図である。
【図１４】ＡＣリアクトルのシミュレーション結果における本実施の形態に係る電磁界解析の所要時間を他の手法の電磁界解析と比較する図表である。
【図１５】ＡＣリアクトルのシミュレーション結果における本実施の形態に係る電磁界解析における時刻ステップあたりの収束計算結果を他の手法の電磁界解析と比較する図表である。
【図１６】インダクタンスの周波数特性の本実施の形態に係るシミュレーション結果を示すグラフである。
【発明を実施するための形態】
【００１６】
以下、本発明の実施の形態について図面を参照して詳細に説明する。なお、以下では図中の同一または相当部分には同一符号を付してその説明は原則的に繰返さないものとする。
【００１７】
（圧粉磁芯の磁気特性測定）
上述のように、本実施の形態では、圧粉磁芯の圧縮加工によって形成されるコア（以下、圧粉磁芯コアとも称する）の磁気特性をシミュレーション可能な電磁界解析を指向している。このため、シミュレーション対象となる圧粉磁芯コアの磁気特性について、まず測定した。
【００１８】
図１は、磁気特性測定の対象とした圧粉磁芯コアの諸元の一覧を示す図表である。
図１に示されるように、内径１２．６（ｍｍ）、外径２０．２６（ｍｍ）、高さ６．４４（ｍｍ）の寸法を有し、密度が７．３５（ｇ／ｃｍ3）である環状の圧粉磁芯コアを対象として、磁気特性測定を実行した。
【００１９】
圧粉磁芯は、個々の粒子が絶縁されており渦電流の影響が非常に小さいため、周波数５０（Ｈｚ）で測定したヒステリシスループを直流ヒステリシス特性として扱うことができる。ヒステリシスループについては、最大磁束密度０．０１（Ｔ）〜０．４（Ｔ）まで、０．０１（Ｔ）刻みで、対称ヒステリシスループを測定した。磁気特性測定時の収束条件は、正弦波の波形率に対する誘導起電力波形の波形率が０．０５（％）以内、最大磁束密度の目標値に対する誤差が０．０５（％）以内、かつ、誘導起電力波形の歪率が２．０（％）以内とした。
【００２０】
図２には、図１に示した圧粉磁芯コアの周波数５０（Ｈｚ）でのヒステリシスループの測定結果が示される。図２中には、上記の測定点のうちの、最大磁束密度が０．１（Ｔ）、０．２（Ｔ）、０．３（Ｔ）、および、０．４（Ｔ）のそれぞれにおけるヒステリシスループ１０１〜１０４が示される。これらの対称ヒステリシスループからプレイモデルの同定を行うと、圧粉磁芯コアの直流ヒステリシス特性の表現に用いることができる。
【００２１】
直流ヒステリシス特性の次に、圧粉磁芯の磁気特性の渦電流による影響を明らかにするために、周波数を変化させてヒステリシスループを測定した。
【００２２】
図３には、図１に示した圧粉磁芯コアの周波数を変化させたときのヒステリシスループの測定結果が示される。ここでは、最大磁束密度を約５０（ｍＴ）としてヒステリシスループを測定するとともに、測定周波数は、発熱を考慮して、１（ｋＨｚ）から１（ＭＨｚ）までの範囲とした。
【００２３】
図３中には、測定周波数１（ｋＨｚ）、５０（ｋＨｚ）、２００（ｋＨｚ）、５００（ｋＨｚ）、８００（ｋＨｚ）、および、１（ＭＨｚ）におけるヒステリシスループ１１１〜１１６が示される。
【００２４】
図３での測定結果から、測定周波数が大きくなるにつれて、渦電流の影響によって損失が増加している、すなわち、ヒステリシスループが膨らむことが確認できる。さらに、最大磁束密度Ｂｍと、当該Ｂｍにおける磁界の強さＨｂから導出される透磁率μが、周波数の増加に伴って低下していることが理解される。
【００２５】
次に、図１〜図３で説明した圧粉磁芯コアに導線（コイル）を巻回することで構成されるインダクタでのインダクタンスの周波数特性について、ＬＣＲメータを用いて測定結果が、図４に示される。
【００２６】
図４には、コイルの巻数は１４として、周波数ｆを１（ｋＨｚ）から８（ＭＨｚ）まで変化させて、インダクタンスを測定した結果が示される。図４の縦軸は、各周波数でのインダクタンスの測定値Ｌを、周波数ｆ＝１（ｋＨｚ）におけるインダクタンスの測定値Ｌｄｃで除算することによって正規化されたインダクタンス（Ｌ／Ｌｄｃ）である。
【００２７】
図４から、ｆ＝１００（ｋＨｚ）を超えた付近から、高周波領域でインダクタンス値が低下していることが読み取れる。これは、図２に示すように、周波数の増大に伴って透磁率μが減少することに起因する。
【００２８】
（磁気特性の等価回路モデル）
次に、本実施の形態に係る電磁界解析で用いられる等価回路モデルについて説明する。
【００２９】
図５は、本実施の形態に係る電磁界解析で用いられるＣａｕｅｒ−Ｉ型回路（以下、単に、Ｃａｕｅｒ回路とも称する）の構成を示す回路図である。
【００３０】
図５に示されるように、Ｃａｕｅｒ回路は、抵抗素子およびインダクタが並列と直列接続を繰り返すことで構成されるはしご型回路（ラダー回路）である。図５では、抵抗素子Ｒ0〜ＲnおよびインダクタＬ0〜Ｌnによって構成されるＣａｕｅｒ回路が示される。最初のインダクタＬ0によって、直流ヒステリシス特性が表現するとともに、Ｌ0以降のインダクタＬ1〜Ｌnによって、渦電流による周波数特性を表現することができる。図５に示されるように、Ｃａｕｅｒ回路の端子間電圧は、磁束密度の時間的変化ｄＢ／ｄｔを示し、電流は、磁界Ｈを示す。
【００３１】
ここで、インダクタＬ0に対応する直流ヒステリシス特性は、プレイモデルによって算出することができる。また、インダクタＬ0以降の、インダクタＬ1〜Ｌnおよび抵抗素子Ｒ0〜Ｒnの回路素子値（抵抗値またはインダクタンス値）については、例えば、遺伝的アルゴリズムによる最適化手法によって、図２に示したヒステリシス特性の測定結果を再現可能な値に決定することができる。
【００３２】
次に、図５に示したＣａｕｅｒ回路の有限要素法への適用を説明する。
図５に示したＣａｕｅｒ回路において、節点電圧に対応する磁束密度Ｂ（ベクトル）を三次元直交座標に分解することで、下記の式（１）のように定義する。以下の各数式中において、ベクトルは太字で表記される。
【００３３】
【数１】
2020154560.tif 000003
【００３４】
式（１）において、下付添字iは、Ｃａｕｅｒ回路の段数を表す。また、Ｃａｕｅｒ回路における各節点の関係から、下記式（２），（３）が得られる。
【００３５】
【数２】
2020154560.tif 000004
【００３６】
【数３】
2020154560.tif 000005
【００３７】
式（２），（３）において、主磁束に対応するＨdc（ベクトル）は、磁化曲線またはプレイモデルを適用することによって求めることができる。さらに、静磁界の有限要素方程式にガラーキン法を適用すると，残差Ｇi（ベクトル）は、下記の式（４）によって示される。
【００３８】
【数４】
2020154560.tif 000006
【００３９】
式（４）中のＮk（ベクトル）は、辺要素基底関数である。式（４）中のＨ0に式（２）を代入することで、以下の式（５）を得ることができる。
【００４０】
【数５】
2020154560.tif 000007
【００４１】
同様に、Ｃａｕｅｒ回路の各段（ｉ）に対応する残差Ｇci（ベクトル）は、下記の式（６）で示される。
【００４２】
【数６】
2020154560.tif 000008
【００４３】
式（６）において、Ｂ0（ベクトル）は、主磁束に対応する磁束密度であり，有限要素法において定義された磁気ベクトルポテンシャルＡの回転により得られる。式（５），（６）にニュートンラフソン法を適用すると、下記の式（７）に示される線形化された方程式が得られる。
【００４４】
【数７】
2020154560.tif 000009
【００４５】
式（７）中のブロック行列は、以下の式（８）〜（１１）によって示される。
【００４６】
【数８】
2020154560.tif 000010
【００４７】
【数９】
2020154560.tif 000011
【００４８】
【数１０】
2020154560.tif 000012
【００４９】
【数１１】
2020154560.tif 000013
【００５０】
なお、ｅαは、α方向（α＝ｘ，ｙ，ｚ）の単位ベクトルであり、式（７）のブロック行列において、式（９）〜（１１）以外の成分は０である。ここで、式（７）を簡易的に、次の式（１２）のように表現する。
【００５１】
【数１２】
2020154560.tif 000014
【００５２】
式（１２）において、ＸFおよびＸCは，それぞれ有限要素法およびＣａｕｅｒ回路に対応する未知変数である。すなわち、変数ＸFは「第１の変数」の一実施例に対応し、変数ＸCは「第２の変数」の一実施例に対応し、式（１２）は、「第１の連立一次方程式」に対応する。Ｓｃｈｕｒ補元を用いることによって、Ｃａｕｅｒ回路に対応する未知変数ＸCを消去することにより、式（１２）を下記の式（１３）に変形することができる。
【００５３】
【数１３】
2020154560.tif 000015
【００５４】
通常、式（１３）内の行列（ＫFF−ＫFC・ＫCC-1・ＫCF）を陽的に求めることは，計算コストの観点から容易ではない。しかしながら、本実施の形態に係る磁気特性モデルにおいては、ＫCC内のＣａｕｅｒ回路に相当する係数行列は各要素内で体積の寄与を除いて同一であるので、逆行列ＫCC-1を容易に算出することができる。
【００５５】
ＫFCおよびＫCFも圧粉磁芯の要素内でのみ非零となる疎行列であるため、式（１２）を縮約した式（１３）の導出は低コストで実現することができる。式（１３）は、「第２の連立一次方程式」の一実施例に対応する。
【００５６】
また、行列（ＫFC・ＫCC-1・ＫCF）で計算されるＣａｕｅｒ回路の影響は圧粉磁芯コアの要素にのみ現れるため、行列（ＫFF−ＫFC・ＫCC-1・ＫCF）およびＫFFにおける非零要素位置は不変である。したがって、縮約した方程式の実装は容易であり、行列ベクトル積や前処理に要する計算コストは、通常のヒステリシス磁界解析に相当する、行列ＫFFのみの方程式での行列ベクトル積や前処理に要する計算コストとほぼ同等である。縮約後の式（１３）を解くことによって、有限要素法に対応する未知変数ＸFを導出した後、式（１４）により、Ｃａｕｅｒ回路に対応する未知変数ＸCを導出することができる。すなわち、式（１３）を解く演算は「第１の演算」に対応する。また、式（１４）では、式（１３）から導出された変数ＸFを式（１２）に代入することによって、変数ＸCが導出される。すなわち、式（１４）を解く演算は「第２の演算」に対応する。
【００５７】
【数１４】
2020154560.tif 000016
【００５８】
導出したＸFおよびＸCから、次の時刻ステップにおける、式（１３）の右辺項である残差ＧF，ＧCを計算することにより、非線形過渡磁界解析を時刻ステップ毎に、定常解が得られるまで繰り返すことができる。なお、式（１３），（１４）中において、行列ＫFFについては各時刻ステップで変化する一方で、行列ＫFC，ＫCF，ＫCCについては、式（９）〜（１１）から理解されるように定数である。したがって、これらのＫFC，ＫCF，ＫCCについては、過渡解析開始時に一度のみ計算すればよいので、計算コストが抑制される。
【００５９】
以上のように、Ｃａｕｅｒ回路による周波数依存磁気特性モデル化と、式（１３）に示す連立一次方程式の求解法を組み合わせることによって、幅広い周波数領域で高精度な磁界解析を、従来法と同程度の計算コストで実施することが可能である。
【００６０】
（解析例）
図６は、本実施の形態に係る電磁界シミュレータの解析対象および解析モデルの第１の例を説明する概念図である。
【００６１】
図６を参照して、リングコア３００は、例えば、圧粉磁芯の圧縮加工によって形成される。本実施の形態では、リングコア３００に図示しない導線（コイル）を巻回することで構成されるリングコアインダクタが解析対象とされる。
【００６２】
当該リングコアインダクタのうち、対称性から周方向に１／２００領域の上半分を切り出した全体の１／４００モデルを解析領域とすることによって、解析モデル２００が得られる。図６の例では、リングコア３００の解析領域３１０は、複数の要素３１１にメッシュ化され、リングコア３００に巻回されるコイル４００の解析領域４１０は、複数の要素４１１に分割される。解析領域４１０は、コイル４００のうちの、解析領域３１０に対応して巻回された部分領域に相当する。以下に例示するシミュレーションでは、メッシュ化された要素の合計数は１４００であり、そのうち、リングコア３００の解析領域３１０における分割数（すなわち、要素３１１の個数）は５５である。有限要素法と示したＣａｕｅｒ回路との連成解析は、各要素において、Ｃａｕｅｒ回路（図５）が配置されるのと等価である。
【００６３】
シミュレーションでは、線形反復法としてＢｉＣＧｓｔａｂ法を使用し、非線形反復法としては、直線探索を導入したニュートンラフソン法を使用した。ニュートンラフソン法では、磁束密度の大きさの修正量 |ΔＢ|が全要素において１０-4（Ｔ）以下になった時点で収束とし、ＢｉＣＧｓｔａｂ法の収束判定値は１０-4とした。さらに、周波数ｆは、１（ｋＨｚ）、２００（ｋＨｚ）、５００（ｋＨｚ）および１（ＭＨｚ）とし，磁化特性測定時（図２および図３）に得られた実測電流波形が、各要素への入力として与えられる。また、直流磁化特性としてプレイモデルを使用し、簡易ＴＰ−ＥＥＣ法を用いて、鎖交磁束から算出する電圧波形が、定常状態に至るまで繰り返し計算を行った。
【００６４】
当該シミュレーションのうちの、上述した、Ｃａｕｅｒ回路による電磁界解析の演算処理は、図７に示されるフローチャートに従って実行されることになる。
【００６５】
図７を参照して、シミュレーションが開始されると、ステップ（以下、単に「Ｓ」と表記する）１１０により、各要素での入力として、シミュレーションを行う周波数での磁化特性測定時に得られた実測電流波形（または、電圧波形）が与えられる。さらに、Ｓ１２０により、上述のように定数である、行列ＫFC，ＫCF，ＫCCが算出される。Ｓ１３０により、最初の時刻ステップにおける行列ＫFFが算出される。
【００６６】
続いて、Ｓ１４０では、行列ＫFC，ＫCF，ＫCCおよび行列ＫFFを用いて、上述の式（１３）に従って、有限要素法に対応する未知変数ＸFが導出される。Ｓ１５０では、Ｓ１４０で求められた未知変数ＸFを用いて、上述の式（１４）に従って、有限要素法に対応する未知変数Ｘｃが導出される。すなわち、Ｓ１４０では「第１の演算」が実行され、Ｓ１５０では「第２の演算」が実行される。
【００６７】
Ｓ１６０では、時刻ステップが進められて、Ｓ１７０において、現在の時刻ステップ（ｔ）での変数ＸF、および、変数Ｘｃを、上述の式（１２）に従って、次の時刻ステップ（ｔ＋Δｔ）での、ＧFおよびＧcが算出される。Ｓ１８０では、現在の時刻ステップ（ｔ）での算出値と、次の時刻ステップ（ｔ＋Δｔ）での算出値によって、上述の収束条件が満たされるか否かの収束判定が実行される。
【００６８】
収束条件の成立時（Ｓ１８０のＹＥＳ判定時）には、Ｓ２００により、最後に求められた算出値が定常解として出力される。収束条件の非成立時（Ｓ１８０のＮＯ判定時）には、Ｓ１９０により、上述のように各時刻ステップで変化する行列ＫFFについて更新した上で、処理がＳ１４０に戻されることにより、時刻ステップを逐次進行させた演算処理を繰り返し実行することができる。
【００６９】
図７に示された演算処理を、図８に例示したハードウェア構成の下で実行することで、本実施の形態に係る電磁界シミュレータを実現することが可能である。
【００７０】
図８を参照して、本実施の形態に係るシミュレータ６００は、コンピュータ等によって構成されて、ＣＰＵ（Central  Processing  Unit）６１０と、メモリ６２０と、入力装置６３０と、出力装置６４０と、外部機器との接続ポート６５０とを備える。ＣＰＵ６１０、メモリ６２０、入力装置６３０、出力装置６４０、および、接続ポート６５０は、内部バス等によって通信接続されており、相互にデータの入出力が可能である。
【００７１】
入力装置６３０は、キーボードおよびマウス等によって構成されて、ユーザからのコマンドおよびデータ等の入力を受ける。出力装置６４０は、液晶ディスプレイ等で構成されてユーザに視認される情報を出力する。接続ポート６５０は、ＵＳＢ（Universal  Serial  Bus)に代表される、外部機器および周辺機器のコネクタであり、外部機器および周辺機器との間でデータの入出力が可能である。
【００７２】
図７に示された演算処理を実行するためのプログラムを予めメモリ６２０に格納し、当該プログラムに従った演算処理をＣＰＵ６１０で実行することにより、本実施の形態に係る電磁界解析が実行される。
【００７３】
（シミュレーション結果１）
次に、本実施の形態に係る電磁界解析による、図１〜図３で説明した圧粉磁芯コアに導線（コイル）を巻回することで構成されるリングコアインダクタのシミュレーション結果について説明する。
【００７４】
図９には、ヒステリシスループの測定結果とシミュレーション結果とを比較するグラフが示される。
【００７５】
図９中には、シミュレーションを行った周波数１（ｋＨｚ）、２００（ｋＨｚ）、５００（ｋＨｚ）および１（ＭＨｚ）のそれぞれについて、図３に示したヒステリシスループの測定結果が実線で示される一方で、シミュレーションによって得られたヒステリシスループが点線で示される。
【００７６】
図９での点線および実線の比較から、本実施の形態に係る電磁界解析によって、低周波から高周波まで幅広い周波数範囲において，周波数に依存して変化するヒステリシス特性を精度良く表現できていることが理解される。
【００７７】
図１０には、インダクタンスの測定結果とシミュレーション結果とを比較するグラフが示される。
【００７８】
図１０中には、図４に示した、インダクタンスの周波数特性の測定結果が特性線１５１で示される。図１０の縦軸は、図４と同様の正規化されたインダクタンス（Ｌ／Ｌｄｃ）である。
【００７９】
さらに、本実施の形態に係る電磁界解析によるシミュレーション結果として得られたインダクタンス値に基づく特性線１５２が示される。シミュレーションでは、コイルの巻数および印加電流は、ＬＣＲメータを用いた、図４中のインダクタンス値と同様の値が設定されている。
【００８０】
また、比較のために、直流ヒステリシスのみを考慮して、渦電流を考慮しない通常の磁界解析、すなわち、Ｃａｕｅｒ回路を導入しない解析によるシミュレーション結果として得られたインダクタンス値に基づく特性線１５３が、さらに示されている。このような解析では、磁化特性の周波数依存性がモデル化されていないため、周波数変化に対して、インダクタンスは一定値となる。
【００８１】
一方で、本実施の形態では、Ｃａｕｅｒ回路に基づく磁気特性モデル化を行っているので、低周波数から高周波数までインダクタンスの周波数特性を高精度に推定できていることが理解される。なお、特性線１５２では、２（ＭＨｚ）を超える周波数帯において、インダクタンスの低下が緩やかになり、実測値（特性線１５１）と異なる傾向を示している。これは、上述したように、図３の測定結果を用いたＣａｕｅｒ回路の最適化計算において、１（ＭＨｚ）までの測定値しか使用していないことが原因であると考えられる。
【００８２】
さらに、式（１２）に代えて式（１３）および（１４）を用いることによる、連立一次方程式の縮約による効果を確認するために、シミュレーションに要した計算時間の比較を行った。具体的には、図６に示された解析モデル２００を共通の解析対象として、周波数ｆ＝１（ｋＨｚ）、２００（ｋＨｚ）、５００（ｋＨｚ）および１（ＭＨｚ）の各々において、渦電流を考慮しない通常の電磁界解析（ケースＡ）、Ｃａｕｅｒ回路を導入するが、本実施の形態のような連立一次方程式の縮約を行わない電磁界解析（ケースＢ）、および、本実施の形態に係る連立一次方程式の縮約を伴う電磁界解析（ケースＣ）を行った。
【００８３】
図１１には、これらケースＡ〜ケースＣでの計算所要時間を比較する図表が示される。
図１１を参照して、ケースＡでは、上述のように、磁化特性の周波数依存性がモデル化されていないため周波数が変化しても、計算所要時間は一定である。ケースＡおよびケースＣでの計算所要時間は、ケースＡでの計算所要時間に対する比で示される。
【００８４】
縮約を行っていないケースＢでは、周波数の増大に伴って計算時間が増加する一方で、縮約によって連立一次方程式の規模が縮小されたケースＣ（本実施の形態）では、シミュレーションを行った１（ｋＨｚ）〜１（ＭＨｚ）までの周波数範囲を通じて、ケースＢよりも計算時間が減少されている。
【００８５】
図１２には、ケースＡ〜ケースＣでの時刻ステップあたりでの収束計算結果を比較する図表が示される。具体的には、ニュートラフソン法（ＮＲ法）における時刻ステップあたりの収束条件が成立するまでの平均反復回数、ならびに、ＢｉＣＧｓｔａｂ法における時刻ステップあたりの、収束条件が成立するまでの平均反復回数および平均計算時間が比較される。
【００８６】
図１２から、非線形反復法であるニュートラフソン法では、ケースＡ〜ケースＣの間で反復回数はほとんど変化が見られないことが理解される。一方で、線形反復法であるＢｉＣＧｓｔａｂ法では、ケースＣ（本実施の形態）による平均反復回数が大幅に減少しており、この結果、図１１で説明した計算時間の短縮が達成されていることが理解される。
【００８７】
以上のように、本実施の形態に係る電磁界解析によれば、ヒステリシスのみを考慮した磁界解析（図１１および図１２のケースＡ）と同程度の計算コストで、磁気特性の周波数依存性を高精度に再現できている。特に、Ｃａｕｅｒ回路の影響が相対的に大きくなる高周波数領域において、連立一次方程式の縮約による計算コスト削減効果が大きい。
【００８８】
（シミュレーション結果２）
次に、本実施の形態に係る電磁界解析によるＡＣリアクトルのシミュレーション結果を説明する。
【００８９】
図１３は、本実施の形態に係る電磁界シミュレータの解析対象および解析モデルの第２の例を説明する概念図である。
【００９０】
図１３を参照して、ＡＣリアクトルは、８個の直方体ブロック５０１、および、４個の円柱ブロック５０２の集合体であるコア５００に対して、円柱部分に図示しない導線（コイル）を巻回することで構成されるＡＣリアクトルが解析対象とされる。各直方体ブロック５０１および各円柱ブロック５０２は、図６のリングコア３００と同様に、例えば、圧粉磁芯の圧縮加工によって形成される。各直方体ブロック５０１および各円柱ブロック５０２の密度は、リングコアと同程度である。
【００９１】
当該ＡＣリアクトルのうち、対称性から全体の１／８モデルを解析領域とすることによって、解析モデル２１０が得られる。
【００９２】
図１３の例では、コア５００の解析領域５１０は、複数の要素５１１にメッシュ化され、コア５００の円柱部分に巻回されるコイル４００の解析領域４１０は、複数の要素４１１に分割される。解析領域４１０は、コイル４００のうちの、解析領域５１０に対応して巻回された部分領域に相当する。又、ブロック間の接合面には、ギャップ要素４５０を設定することにより、接着剤に対応した微小な空隙を模擬した。以下に例示するシミュレーションでは、メッシュ化された要素の合計数は１９７８０であり、そのうち、コア５００の解析領域５１０における分割数（すなわち、要素５１１の個数）は５３６０である。
【００９３】
図１３に示した解析モデルに対して、本実施の形態に係る電磁界解析によって、上述のリングコアインダクタと同様のシミュレーションを行った。すなわち、図１３に示された解析モデルを共通の解析対象として、周波数ｆ＝１（ｋＨｚ）、２００（ｋＨｚ）、５００（ｋＨｚ）および１（ＭＨｚ）の各々において、上記と同様のケースＡ〜ケースＣでのシミュレーションを行った。
【００９４】
図１４には、図１１と同様の本実施の形態に係る電磁界解析の所要時間を他の手法の電磁界解析と比較する図表が示される。図１４においても、図１１と同様に、ケースＡでの周波数依存しない計算所要時間を基準値（１．０）として、縮約無しでＣａｕｅｒ回路をシミュレーションするケースＢ、および、縮約を伴うケースＣ（本実施の形態）による各周波数での計算所要時間を、上記基準時間に対する比で表している。
【００９５】
同様に、図１５には、図１２と同様の本実施の形態に係る電磁界解析における時刻ステップあたりの収束計算結果を他の手法の電磁界解析と比較する図表が示される。図１５においても、図１２と同様に、ニュートラフソン法における時刻ステップあたりの収束条件が成立するまでの平均反復回数、ならびに、ＢｉＣＧｓｔａｂ法における時刻ステップあたりの、収束条件が成立するまでの平均反復回数および平均計算時間が、ケースＡ〜ケースＣの間で比較される。
【００９６】
図１４においても、ケースＢでは、周波数が増大するにつれて、平均反復回数および計算時間が増大しているが、ケースＣでは、平均反復回数および計算時間の増大が抑制されていることが理解できる。
【００９７】
図６の解析モデル（リングコアインダクタ）では全要素数に対する圧粉磁芯コアの要素数が約４％（５５／１４００）であるのに対し、図１３の解析モデル（ＡＣリアクトル）では、全要素数に対して圧粉磁芯コアの要素数が２７％（５３６０／１９７８０）であり、圧粉磁芯コアの要素が占める割合が大きくなっている。このため、図１４および図１５では、図１１および図１２と比較して、ケースＢにおける、周波数の増大に対する平均反復回数および計算時間の増大が顕著になっている。
【００９８】
一方で、ケースＣにおける周波数の増大に対する平均反復回数および計算時間の増大は、図１４および図１５においても、図１１および図１２とほぼ同等である。したがって、解析モデル中の全要素数における圧粉磁芯コアの要素数が占める割合が大きいほど、本実施の形態に係る電磁界解析による計算コストの削減効果が高くなることが理解される。
【００９９】
また、図１６には、図１３に示した解析モデル（ＡＣリアクトル）を対象とした、本実施の形態に係る電磁界解析によるシミュレーション結果として得られたインダクタンス値が示される。
【０１００】
図１６においても、各周波数でのシミュレーション結果として得られたインダクタンス値Ｌが、周波数ｆ＝１（ｋＨｚ）でのシミュレーション結果として得られたインダクタンス値Ｌｃを基準として正規化された値（Ｌ／Ｌｄｃ）が縦軸に示される。
【０１０１】
図１６の周波数特性についても、図１０に示したリングコアでの周波数特性と同様に、１００（ｋＨｚ）超えた周波数からインダクタンスの低下が見られるが、ブロック間のギャップの影響でインダクタンスの低下は緩やかであり、定性的に妥当な周波数特性がシミュレーションから得られていることが理解される。
【０１０２】
本実施の形態に係る電磁界解析は、有限要素法におけるメッシュ作成を適切に行うことにより、任意の形状のコアに適用することが可能である。また、本実施の形態に係る電磁界解析は、上述のように、圧縮成形されたコアを解析対象とするときに大きな効果を享受できるものであるが、解析対象の材質は任意である。すなわち、他の磁性材料を用いて射出成形によって形成されたコア、又は、薄膜の磁性材料や板状の磁性材料を積層したコア等に対しても、本実施の形態に係る電磁界解析を適用することが可能である。
【０１０３】
以上のように本発明の実施の形態について説明を行ったが、上述の実施の形態を様々に変形することも可能である。また、本発明の範囲は上述の実施の形態に限定されるものではない。本発明の範囲は、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更を含むことが意図される。
【符号の説明】
【０１０４】
１０１〜１０４，１１１〜１１６  ヒステリシスループ、１５１〜１５３  特性線（インダクタンス周波数特性）、２００，２１０  解析モデル、３００  リングコア、３１０，４１０，５１０  解析領域、３１１，４１１，５１１  要素、４００  コイル、４５０  ギャップ要素、５００  コア、５０１  直方体ブロック、５０２  円柱ブロック、６００  シミュレータ、６２０  メモリ、６３０  入力装置、６４０  出力装置、６５０  接続ポート。
【図１】
2020154560.tif 000017
【図２】
2020154560.tif 000018
【図３】
2020154560.tif 000019
【図４】
2020154560.tif 000020
【図５】
2020154560.tif 000021
【図６】
2020154560.tif 000022
【図７】
2020154560.tif 000023
【図８】
2020154560.tif 000024
【図９】
2020154560.tif 000025
【図１０】
2020154560.tif 000026
【図１１】
2020154560.tif 000027
【図１２】
2020154560.tif 000028
【図１３】
2020154560.tif 000029
【図１４】
2020154560.tif 000030
【図１５】
2020154560.tif 000031
【図１６】
2020154560.tif 000032