(57)【要約】【課題】寸法共鳴現象が発生する磁性材料の磁場解析の精度を向上させる。【解決手段】記憶部１１は、コアの形状を分割することで生成された複数の要素領域を含む形状モデル１３と、コアに用いる磁性材料の物理的特性を示す特性情報１４と、コアに巻かれるコイルに流れる電流の時間変化を示すコイル電流情報１５とを記憶する。処理部１２は、特性情報１４と第１の時刻におけるコイルの電流密度１６ａとを用いて、第１の時刻における磁場に関する指標値１６ｂ，１６ｃ，１６ｄを算出する。処理部１２は、指標値１６ｂ，１６ｃ，１６ｄを用いて要素領域の電荷密度１７を算出する。処理部１２は、特性情報１４と第２の時刻におけるコイルの電流密度１８ａと電荷密度１７とを用いて、第２の時刻における磁場に関する指標値１８ｂ，１８ｃ，１８ｄを算出する。【選択図】図１
2018136753.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータに、
コアの形状を分割することで生成された複数の要素領域を含む形状モデルと、前記コアに用いる磁性材料の物理的特性を示す特性情報と、前記コアに巻かれるコイルに流れる電流の時間変化を示すコイル電流情報とを取得し、
前記コイル電流情報に基づいて第１の時刻における前記コイルの第１の電流密度を特定し、前記特性情報と前記第１の電流密度とを用いて、前記複数の要素領域によって特定される前記形状モデル上の複数の位置に対応させて前記第１の時刻における磁場に関する複数の第１の指標値を算出し、
前記複数の第１の指標値を用いて、前記複数の要素領域それぞれの電荷密度を算出し、
前記コイル電流情報に基づいて前記第１の時刻より後の第２の時刻における前記コイルの第２の電流密度を特定し、前記特性情報と前記第２の電流密度と前記複数の要素領域それぞれの前記電荷密度とを用いて、前記複数の位置に対応させて前記第２の時刻における磁場に関する複数の第２の指標値を算出する、
処理を実行させる磁性材料シミュレーションプログラム。
【請求項２】
前記電荷密度の算出では、前記複数の第１の指標値を用いて、前記第１の時刻における前記複数の要素領域それぞれの第３の電流密度を算出し、前記第３の電流密度を用いて、前記複数の要素領域それぞれの前記電荷密度を算出する、
請求項１記載の磁性材料シミュレーションプログラム。
【請求項３】
前記コンピュータに更に、
前記複数の第２の指標値と前記電荷密度とを用いて、前記第２の時刻における前記複数の要素領域それぞれの第４の電流密度を算出し、前記第４の電流密度と前記電荷密度とを用いて、前記複数の要素領域それぞれの他の電荷密度を算出する、
処理を実行させる請求項１記載の磁性材料シミュレーションプログラム。
【請求項４】
前記特性情報は、前記磁性材料に含まれる複数の結晶粒それぞれの厚さと前記複数の結晶粒の間に存在する粒界層の厚さとの間の関係を示すサイズ情報を含む、
請求項１記載の磁性材料シミュレーションプログラム。
【請求項５】
前記複数の第１の指標値および前記複数の第２の指標値は、ベクトルポテンシャルおよびスカラポテンシャルの少なくとも一方を含む、
請求項１記載の磁性材料シミュレーションプログラム。
【請求項６】
コンピュータが実行する磁性材料シミュレーション方法であって、
コアの形状を分割することで生成された複数の要素領域を含む形状モデルと、前記コアに用いる磁性材料の物理的特性を示す特性情報と、前記コアに巻かれるコイルに流れる電流の時間変化を示すコイル電流情報とを取得し、
前記コイル電流情報に基づいて第１の時刻における前記コイルの第１の電流密度を特定し、前記特性情報と前記第１の電流密度とを用いて、前記複数の要素領域によって特定される前記形状モデル上の複数の位置に対応させて前記第１の時刻における磁場に関する複数の第１の指標値を算出し、
前記複数の第１の指標値を用いて、前記複数の要素領域それぞれの電荷密度を算出し、
前記コイル電流情報に基づいて前記第１の時刻より後の第２の時刻における前記コイルの第２の電流密度を特定し、前記特性情報と前記第２の電流密度と前記複数の要素領域それぞれの前記電荷密度とを用いて、前記複数の位置に対応させて前記第２の時刻における磁場に関する複数の第２の指標値を算出する、
磁性材料シミュレーション方法。
【請求項７】
コアの形状を分割することで生成された複数の要素領域を含む形状モデルと、前記コアに用いる磁性材料の物理的特性を示す特性情報と、前記コアに巻かれるコイルに流れる電流の時間変化を示すコイル電流情報とを記憶する記憶部と、
前記コイル電流情報に基づいて第１の時刻における前記コイルの第１の電流密度を特定し、前記特性情報と前記第１の電流密度とを用いて、前記複数の要素領域によって特定される前記形状モデル上の複数の位置に対応させて前記第１の時刻における磁場に関する複数の第１の指標値を算出し、前記複数の第１の指標値を用いて、前記複数の要素領域それぞれの電荷密度を算出し、前記コイル電流情報に基づいて前記第１の時刻より後の第２の時刻における前記コイルの第２の電流密度を特定し、前記特性情報と前記第２の電流密度と前記複数の要素領域それぞれの前記電荷密度とを用いて、前記複数の位置に対応させて前記第２の時刻における磁場に関する複数の第２の指標値を算出する処理部と、
を有する磁性材料シミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は磁性材料シミュレーションプログラム、磁性材料シミュレーション方法および磁性材料シミュレーション装置に関する。
【背景技術】
【０００２】
コンピュータの演算能力の向上に伴い、コンピュータを利用して様々な物理現象のシミュレーションが行われるようになっている。物理現象のシミュレーションの１つとして、磁性材料を用いたコアと当該コアの周りに巻かれたコイルとを含む電気部品について、コイルに電流を流したときに発生する磁場を解析する磁場解析が挙げられる。磁場解析のシミュレーションでは、数値解析法として有限要素法を用いることがある。
【０００３】
有限要素法では、解析対象の形状を「要素」もしくは「メッシュ」と言われる小領域に分割したモデルを生成し、要素の中心・要素の頂点である節点・要素の辺などの位置に対して、未知の物理値を示す変数を割り当てる。そして、物理現象を示す基礎方程式とモデル上に割り当てた複数の変数から連立方程式（通常、大規模かつ疎な係数行列をもつ連立一次方程式）を生成し、連立方程式を解くことで変数の近似解を求める。
【０００４】
なお、有限要素法を用いた磁気シールド解析方法が提案されている。提案の磁気シールド解析方法では、第１ステップと第２ステップと第３ステップとを含む解析ルーチンを実行する。第１ステップでは、解析対象の形状に含まれる複数のサブ領域それぞれに対して磁界状態を設定する。第２ステップでは、第１ステップで設定された磁界状態と積分方程式を用いて境界条件を算出する。第３ステップでは、第２ステップで算出された境界条件と有限要素法を用いて各サブ領域の磁界状態を更新する。上記の解析ルーチンを、時間軸に沿って時刻を変化させながら繰り返し実行する。また、渦電流の影響を考慮して上記の３つのステップを再実行し、各サブ領域の磁界状態を更新する。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００６−３５１６２２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
ところで、マンガン亜鉛（ＭｎＺｎ）コアなど多結晶質の磁性材料を用いたコアでは、コアに巻かれたコイルに高周波電流を流すと寸法共鳴現象が発生することがある。寸法共鳴現象は、磁性材料に含まれる複数の結晶粒の間に存在する隙間がコンデンサのように振る舞うことで、特定の周波数で磁性材料の透磁率が急激に増大する物理現象である。コアの透磁率の虚部が増大するとコア内の損失が大きくなるため、磁性材料を用いたコアの磁場解析では寸法共鳴現象を考慮することが好ましい。しかし、従来の有限要素法の磁場解析では、上記の磁性材料の物理的性質がモデルや基礎方程式に組み込まれておらず、シミュレーションによって寸法共鳴現象を再現することが難しかったという問題がある。
【０００７】
１つの側面では、本発明は、寸法共鳴現象が発生する磁性材料の磁場解析の精度を向上できる磁性材料シミュレーションプログラム、磁性材料シミュレーション方法および磁性材料シミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【０００８】
１つの態様では、コンピュータに以下の処理を実行させる磁性材料シミュレーションプログラムが提供される。コアの形状を分割することで生成された複数の要素領域を含む形状モデルと、コアに用いる磁性材料の物理的特性を示す特性情報と、コアに巻かれるコイルに流れる電流の時間変化を示すコイル電流情報とを取得する。コイル電流情報に基づいて第１の時刻におけるコイルの第１の電流密度を特定し、特性情報と第１の電流密度とを用いて、複数の要素領域によって特定される形状モデル上の複数の位置に対応させて第１の時刻における磁場に関する複数の第１の指標値を算出する。複数の第１の指標値を用いて、複数の要素領域それぞれの電荷密度を算出する。コイル電流情報に基づいて第１の時刻より後の第２の時刻におけるコイルの第２の電流密度を特定し、特性情報と第２の電流密度と複数の要素領域それぞれの電荷密度とを用いて、複数の位置に対応させて第２の時刻における磁場に関する複数の第２の指標値を算出する。
【０００９】
また、１つの態様では、コンピュータが実行する磁性材料シミュレーション方法が提供される。また、１つの態様では、磁性材料シミュレーション装置が提供される。
【発明の効果】
【００１０】
１つの側面では、寸法共鳴現象が発生する磁性材料の磁場解析の精度が向上する。
【図面の簡単な説明】
【００１１】
【図１】第１の実施の形態の磁性材料シミュレーション装置の例を示す図である。
【図２】シミュレーション装置のハードウェア例を示すブロック図である。
【図３】トロイダルコアとコイルの形状例を示す図である。
【図４】コイルの周波数とトロイダルコアの比透磁率の関係例を示すグラフである。
【図５】磁性材料の微細構造の例を示す図である。
【図６】磁性材料の微細構造と等価な電気回路の例を示す図である。
【図７】トロイダルコアとコイルの有限要素モデルの例を示す図である。
【図８】有限要素モデルに対する変数の割り当て例を示す図である。
【図９】シミュレーション装置の機能例を示すブロック図である。
【図１０】節点テーブルと要素テーブルの例を示す図である。
【図１１】パラメータテーブルの例を示す図である。
【図１２】結果テーブルの例を示す図である。
【図１３】シミュレーションの手順例を示すフローチャートである。
【図１４】寸法共鳴現象を考慮しないシミュレーション例を示すグラフである。
【図１５】寸法共鳴現象を考慮したシミュレーション例を示すグラフである。
【図１６】磁束密度可視化画像の例を示す図である。
【図１７】電流密度可視化画像の例を示す図である。
【発明を実施するための形態】
【００１２】
以下、本実施の形態を図面を参照して説明する。
［第１の実施の形態］
第１の実施の形態を説明する。
【００１３】
図１は、第１の実施の形態の磁性材料シミュレーション装置の例を示す図である。
第１の実施の形態の磁性材料シミュレーション装置１０は、磁性材料を用いたコアと当該コアの周りに巻いたコイルとを含む電気部品の磁場解析を、有限要素法を用いたコンピュータシミュレーションとして実行する。磁性材料シミュレーション装置１０は、ユーザが操作するクライアント装置でもよいしサーバ装置でもよい。
【００１４】
磁性材料シミュレーション装置１０は、記憶部１１および処理部１２を有する。記憶部１１は、ＲＡＭ（Random Access Memory）などの揮発性の記憶装置、または、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性の記憶装置である。処理部１２は、例えば、ＣＰＵ（Central Processing Unit）やＤＳＰ（Digital Signal Processor）などのプロセッサである。ただし、処理部１２は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの特定用途の電子回路を含んでもよい。プロセッサは、ＲＡＭなどのメモリに記憶されたプログラムを実行する。例えば、磁性材料シミュレーションプログラムが実行される。なお、複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」と言うことがある。
【００１５】
記憶部１１は、シミュレーションに用いる形状モデル１３、特性情報１４およびコイル電流情報１５を記憶する。形状モデル１３、特性情報１４およびコイル電流情報１５は、ユーザが作成してもよいしモデリングソフトウェアなどのソフトウェアが生成してもよい。また、形状モデル１３、特性情報１４およびコイル電流情報１５は、磁性材料シミュレーション装置１０で作成してもよいし他の装置から取得してもよい。
【００１６】
形状モデル１３は、コアの形状を分割することで生成された複数の要素領域を含む有限要素モデルである。コアの例として、ドーナツ型のトロイダルコアが挙げられる。形状モデル１３は、コアに巻かれるコイルの形状に対応する要素領域や、コアおよびコイルの周辺の空間に対応する要素領域などを更に含んでもよい。複数の要素領域は、例えば、三角形・四角形・四面体・六面体など、同一種類の二次元または三次元の領域である。
【００１７】
特性情報１４は、コアに用いる磁性材料の物理的特性を示す。特性情報１４は、例えば、磁性材料の導電率や誘電率を含む。また、特性情報１４は、例えば、磁性材料に含まれる複数の結晶粒それぞれの厚さとそれら複数の結晶粒の間に存在する粒界層の厚さとの間の関係を示すサイズ情報を含む。コイル電流情報１５は、コアに巻かれるコイルに流れる電流の時間変化を示す。電流の時間変化は、振幅や周波数によって連続的に表現されてもよいし、各時刻の電流密度として離散的に表現されてもよい。
【００１８】
処理部１２は、コイル電流情報１５に基づいて、シミュレーション上の時刻＃１におけるコイルを流れる電流の電流密度１６ａを特定する。コイル電流情報１５が振幅や周波数によって表現されている場合、処理部１２は、時刻＃１における電流密度１６ａを計算により求める。コイル電流情報１５が離散的に表現されている場合、処理部１２は、複数の時刻の電流密度の中から時刻＃１の電流密度１６ａを選択する。
【００１９】
処理部１２は、特性情報１４と電流密度１６ａとを用いて、形状モデル１３上の複数の位置に対応させて時刻＃１の磁場に関する複数の指標値を算出する。形状モデル１３上の複数の位置は、要素領域の中心・要素領域を囲む節点・要素領域を囲む辺など、形状モデル１３に含まれる複数の要素領域に基づいて特定される位置である。例えば、処理部１２は、要素領域を囲む辺ｅ１に対して指標値１６ｂを算出し、辺ｅ２に対して指標値１６ｃを算出し、辺ｅ３に対して指標値１６ｄを算出する。磁場に関する指標値は、シミュレーションによって算出したい物理現象を示す数値であり、例えば、ベクトルポテンシャルなどである。例えば、処理部１２は、所定の計算式＃１と特性情報１４と電流密度１６ａから、形状モデル１３上の複数の位置に対して割り当てた複数の変数を含む連立方程式を生成し、連立方程式を解くことで当該複数の変数の解として複数の指標値を算出する。
【００２０】
次に、処理部１２は、指標値１６ｂ，１６ｃ，１６ｄを含む時刻＃１の複数の指標値を用いて、形状モデル１３に含まれる複数の要素領域それぞれに蓄えられた電荷の電荷密度を算出する。例えば、処理部１２は、ある要素領域の電荷密度１７を、当該要素領域を囲む辺ｅ１，ｅ２，ｅ３に対応する指標値１６ｂ，１６ｃ，１６ｄを用いて算出する。電荷密度は、例えば、上記の計算式＃１とは異なる所定の計算式＃２に従って算出される。このとき、処理部１２は、時刻＃１の複数の指標値を用いて、形状モデル１３に含まれる複数の要素領域それぞれに流れる電流の電流密度を算出し、その電流密度から複数の要素領域それぞれの電荷密度を算出してもよい。
【００２１】
次に、処理部１２は、コイル電流情報１５に基づいて、時刻＃１より後であるシミュレーション上の時刻＃２におけるコイルを流れる電流の電流密度１８ａを特定する。処理部１２は、特性情報１４と電流密度１８ａと複数の要素領域それぞれの電荷密度とを用いて、形状モデル１３上の複数の位置に対応させて時刻＃２の磁場に関する複数の指標値を算出する。例えば、処理部１２は、辺ｅ１に対して指標値１８ｂを算出し、辺ｅ２に対して指標値１８ｃを算出し、辺ｅ３に対して指標値１８ｄを算出する。
【００２２】
指標値１８ｂ，１８ｃ，１８ｄの種類は指標値１６ｂ，１６ｃ，１６ｄと同様であり、例えば、ベクトルポテンシャルなどである。ここで使用される計算式は、基本的に上記の計算式＃１と同様である。ただし、時刻＃１のシミュレーションの際に算出された電荷密度１７が計算式の中で用いられる。例えば、処理部１２は、特性情報１４と電流密度１８ａと複数の要素領域それぞれの電荷密度から、上記の複数の変数を含む連立方程式を生成し、連立方程式を解くことで当該複数の変数の解として複数の指標値を算出する。
【００２３】
第１の実施の形態の磁性材料シミュレーション装置１０によれば、有限要素法を用いた磁性材料シミュレーションにおいて、時刻＃１におけるシミュレーション結果から各要素領域に蓄えられた電荷の電荷密度が算出される。そして、算出された電荷密度が、時刻＃１の後（例えば、時刻＃１の直後）の時刻＃２におけるシミュレーションに用いられる。これにより、磁性材料に含まれる複数の結晶粒の間に存在する隙間がコンデンサのように振る舞うことで、特定の周波数で透磁率が急激に増大するという寸法共鳴現象を再現することができる。よって、磁性材料シミュレーションの精度が向上する。
【００２４】
［第２の実施の形態］
次に、第２の実施の形態を説明する。
第２の実施の形態のシミュレーション装置１００は、ＭｎＺｎトロイダルコアとコイルとを含む電気部品の磁場解析を、コンピュータシミュレーションとして実行する。
【００２５】
図２は、シミュレーション装置のハードウェア例を示すブロック図である。
シミュレーション装置１００は、ＣＰＵ１０１、ＲＡＭ１０２、ＨＤＤ１０３、画像信号処理部１０４、入力信号処理部１０５、媒体リーダ１０６および通信インタフェース１０７を有する。上記ユニットはバスに接続される。
【００２６】
ＣＰＵ１０１は、プログラムの命令を実行する演算回路を含むプロセッサである。ＣＰＵ１０１は、ＨＤＤ１０３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ１０２にロードし、プログラムを実行する。なお、ＣＰＵ１０１は複数のプロセッサコアを備えてもよく、シミュレーション装置１００は複数のプロセッサを備えてもよく、以下の処理を複数のプロセッサまたはプロセッサコアを用いて並列に実行してもよい。また、複数のプロセッサの集合を「マルチプロセッサ」または「プロセッサ」と言うことがある。
【００２７】
ＲＡＭ１０２は、ＣＰＵ１０１が実行するプログラムやＣＰＵ１０１が演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。なお、シミュレーション装置１００は、ＲＡＭ以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。
【００２８】
ＨＤＤ１０３は、ＯＳ（Operating System）やアプリケーションソフトウェアなどのソフトウェアのプログラム、および、データを記憶する不揮発性の記憶装置である。なお、シミュレーション装置１００は、フラッシュメモリやＳＳＤ（Solid State Drive）などの他の種類の記憶装置を備えてもよく、複数の不揮発性の記憶装置を備えてもよい。
【００２９】
画像信号処理部１０４は、ＣＰＵ１０１からの命令に従って、シミュレーション装置１００に接続されたディスプレイ３１に画像を出力する。ディスプレイ３１としては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、プラズマディスプレイ、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなど、任意の種類のディスプレイを用いることができる。
【００３０】
入力信号処理部１０５は、シミュレーション装置１００に接続された入力デバイス３２から入力信号を取得し、ＣＰＵ１０１に出力する。入力デバイス３２としては、マウスやタッチパネルやタッチパッドやトラックボールなどのポインティングデバイス、キーボード、リモートコントローラ、ボタンスイッチなどを用いることができる。また、シミュレーション装置１００に複数の種類の入力デバイスが接続されていてもよい。
【００３１】
媒体リーダ１０６は、記録媒体３３に記録されたプログラムやデータを読み取る読み取り装置である。記録媒体３３として、例えば、磁気ディスク、光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。磁気ディスクには、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤが含まれる。光ディスクには、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）が含まれる。
【００３２】
媒体リーダ１０６は、例えば、記録媒体３３から読み取ったプログラムやデータを、ＲＡＭ１０２やＨＤＤ１０３などの他の記録媒体にコピーする。読み取られたプログラムは、例えば、ＣＰＵ１０１によって実行される。なお、記録媒体３３は可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体３３やＨＤＤ１０３を、コンピュータ読み取り可能な記録媒体と言うことがある。
【００３３】
通信インタフェース１０７は、ネットワーク３４に接続され、ネットワーク３４を介して他の装置と通信を行うインタフェースである。通信インタフェース１０７は、スイッチなどの通信装置とケーブルで接続される有線通信インタフェースでもよいし、基地局と無線リンクで接続される無線通信インタフェースでもよい。
【００３４】
次に、解析対象の電気部品の物理的性質について説明する。
図３は、トロイダルコアとコイルの形状例を示す図である。
解析対象の電気部品は、トロイダルコア４１とコイル４２を含む。トロイダルコア４１には、多結晶質の磁性材料であるＭｎＺｎ焼結磁石が用いられている。焼結磁石は、例えば、磁性粉末を高温で焼き固めることで作られる。トロイダルコア４１は、中心に穴の空いたドーナツ型の形状をもつ。コイル４２は、トロイダルコア４１に巻かれた電線である。コイル４２は、トロイダルコア４１の穴を通るように螺旋状に巻かれている。コイル４２は、例えば、トロイダルコア４１の表面に接している。コイル４２には、外部から交流電流が供給される。コイル４２を流れる電流は、高周波電流であることがある。
【００３５】
ここで、トロイダルコア４１の比透磁率（真空の透磁率に対するトロイダルコア４１の透磁率の比）は、コイル４２に流れる電流の周波数に応じて変化する。特に、比較的高い特定の周波数の電流がコイル４２に流れると、トロイダルコア４１では、比透磁率が急激に増大するという寸法共鳴現象が発生する。寸法共鳴現象は、磁性材料に含まれる複数の結晶粒の間に存在する隙間（粒界層）がコンデンサのように振る舞うことで発生する。
【００３６】
図４は、コイルの周波数とトロイダルコアの比透磁率の関係例を示すグラフである。
グラフ５１は、コイル４２を流れる電流の周波数とトロイダルコア４１の比透磁率との間の関係を示す。グラフ５１に示すように、トロイダルコア４１の比透磁率の実部および虚部はそれぞれ、比較的高い特定の周波数においてピークをもつ。例えば、比透磁率の実部は、約８０００ｋＨｚにおいて約１５００の最大値をとる。また、例えば、比透磁率の虚部は、約９０００ｋＨｚにおいて約１３００の最大値をとる。
【００３７】
トロイダルコア４１の比透磁率の虚部が大きくなると、トロイダルコア４１の単位体積当たりの損失は増大する。そのため、電気部品の設計者は、コイル４２の周波数とトロイダルコア４１の損失との間の関係を把握しておくことが好ましい。シミュレーション装置１００は、このような寸法共鳴現象を反映した磁場解析を行う。
【００３８】
図５は、磁性材料の微細構造の例を示す図である。
トロイダルコア４１に用いる磁性材料では、焼結の過程で、結晶粒６１を含む複数の結晶粒とそれら複数の結晶粒の間に存在する粒界層６２とが形成される。結晶粒６１は、比較的高い導電率、すなわち、比較的低い抵抗率を有する。粒界層６２は、比較的低い導電率、すなわち、比較的高い抵抗率を有する。導電率の高い結晶粒６１を導電率の低い粒界層６２で囲むことで、渦電流の発生が抑制される。ただし、特定の高周波電流がコイル４２に流れると、結晶粒６１に渦電流が発生しやすくなる。
【００３９】
磁性材料の中で、サイズが均一である複数の結晶粒が規則的に並んでいると仮定する。磁性材料をモデル化するにあたり、１つの結晶粒とその周囲に存在する粒界層を１つの単位領域として抽出する。ここでは説明を簡単にするため、二次元の単位領域を考える。図５に示すように、結晶粒６１と結晶粒６１の三辺を囲む粒界層６２とが１つの単位領域として抽出される。ある方向の粒界層の厚さは、隣接する結晶粒の間の距離に相当する。
【００４０】
このような単位領域について、電位差Ｖ1，Ｖ2、電流密度ｉ1，ｉ2、導電率σ1，σ2および距離Ｌ，Ｄが定義される。電位差Ｖ1は、結晶粒６１の一辺とその対辺との間の電位の差である。電位差Ｖ2は、粒界層６２の一辺とその対辺（結晶粒６１の一辺と隣の結晶粒の辺）との間の電位の差である。電流密度ｉ1は、結晶粒６１の中の電流密度である。電流密度ｉ2は、粒界層６２の中の電流密度である。導電率σ1は、結晶粒６１の導電率である。導電率σ2は、粒界層６２の導電率である。距離Ｌは、結晶粒６１の一辺とその対辺との間の距離、すなわち、結晶粒６１の厚さである。距離Ｄは、粒界層６２の一辺とその対辺との間の距離、すなわち、粒界層６２の厚さである。
【００４１】
この単位領域の電気的性質は、これと等価な電気回路によって表現できる。
図６は、磁性材料の微細構造と等価な電気回路の例を示す図である。
上記の磁性材料の単位領域では、粒界層６２が電荷を蓄えるコンデンサとして振る舞うことがある。そこで、上記の単位領域の電気的性質は、図６の電気回路と等価となる。この電気回路は、抵抗回路６３，６４およびコンデンサ６５を含む。抵抗回路６３に対して、抵抗回路６４とコンデンサ６５が並列に接続されている。抵抗回路６３は結晶粒６１に対応し、抵抗回路６４およびコンデンサ６５は粒界層６２に対応する。
【００４２】
すなわち、抵抗回路６３の一端と他端との間の電位差はＶ1である。抵抗回路６４の一端と他端の間の電位差はＶ2である。更に、図６の電気回路について、抵抗Ｒ1，Ｒ2、電流Ｉ1，Ｉ2，ＩC、電荷ＱおよびキャパシタンスＣが定義される。抵抗Ｒ1は抵抗回路６３の抵抗、すなわち、結晶粒６１の抵抗である（単位はΩ（オーム））。抵抗Ｒ2は抵抗回路６４の抵抗、すなわち、粒界層６２の抵抗である。電流Ｉ1は抵抗回路６３を流れる電流である。電流Ｉ2は抵抗回路６４を流れる電流である。電流ＩCはコンデンサ６５を流れる電流でありＩ1−Ｉ2に相当する。電荷Ｑはコンデンサ６５に蓄えられる電荷、すなわち、粒界層６２の電荷である（単位はＣ（クーロン））。キャパシタンスＣはコンデンサ６５のキャパシタンス（静電容量）、すなわち、粒界層６２のキャパシタンスである（単位はＦ（ファラド））。また、電荷密度ｑが定義される。電荷密度ｑはコンデンサ６５の電荷密度、すなわち、粒界層６２の電荷密度である（単位はＣ／ｍ2）。
【００４３】
上記の電気回路モデルに基づいて、抵抗Ｒ1，Ｒ2、キャパシタンスＣおよび電荷Ｑは数式（１）〜（４）のように定義することができる。数式（３）においてεrは粒界層６２の比誘電率であり、ε0は真空の誘電率である。
【００４４】
【数１】
2018136753.tif 000003
【００４５】
【数２】
2018136753.tif 000004
【００４６】
【数３】
2018136753.tif 000005
【００４７】
【数４】
2018136753.tif 000006
【００４８】
結晶粒６１の電位差Ｖ1と粒界層６２の電位差Ｖ2の合計である合計電位Ｖは、上記の数式（１），（２）から、数式（５）のように算出される。Ｅを単位領域の平均電界とすると、距離Ｄは距離Ｌより十分に小さいため、合計電位Ｖは数式（６）のように近似される。数式（５），（６）から、平均電界Ｅは数式（７）のように算出される。数式（７）において寸法比ｒは、結晶粒６１の距離Ｌに対する粒界層６２の距離Ｄの比である。
【００４９】
【数５】
2018136753.tif 000007
【００５０】
【数６】
2018136753.tif 000008
【００５１】
【数７】
2018136753.tif 000009
【００５２】
また、前述した電流Ｉ1，Ｉ2，ＩCの間の関係を示す電流保存の式として数式（８）が成立し、数式（８）と数式（４）から電流密度保存の式として数式（９）が成立する。
【００５３】
【数８】
2018136753.tif 000010
【００５４】
【数９】
2018136753.tif 000011
【００５５】
また、抵抗回路６４の電位差とコンデンサ６５の電位差が共にＶ2になるという関係から数式（１０）が成立する。この数式（１０）と数式（２）〜（４）から、電流密度ｉ2と電荷密度ｑとの間に数式（１１）が成立する。数式（７）の電流密度ｉ2に数式（１１）を代入すると、平均電界Ｅは数式（１２）のように算出される。数式（１２）を変形すると、電流密度ｉ1は数式（１３）のように算出される。
【００５６】
【数１０】
2018136753.tif 000012
【００５７】
【数１１】
2018136753.tif 000013
【００５８】
【数１２】
2018136753.tif 000014
【００５９】
【数１３】
2018136753.tif 000015
【００６０】
粒界層６２に蓄えられる電荷の電荷密度ｑは、結晶粒６１の電流密度ｉ1と粒界層６２の電流密度ｉ2の差を時間積分することで算出できる。よって、数式（１１）を利用して、電荷密度ｑは数式（１４）のように算出される。このようにして、結晶粒６１の電流密度ｉ1と粒界層６２の電荷密度ｑとの間の関係が定義される。
【００６１】
【数１４】
2018136753.tif 000016
【００６２】
次に、連立方程式を生成するための磁場解析の基礎方程式について説明する。粒界層６２のキャパシタンス効果を考慮しない場合、磁場解析の基礎方程式として数式（１５），（１６）を用いることができる。数式（１５）は磁場の支配方程式であり、数式（１６）は電流密度保存の式である。数式（１５），（１６）においてσは導電率、Ａはベクトルポテンシャル、φはスカラポテンシャルである。また、数式（１５）においてμは透磁率、Ｊ0はコイル４２の電流密度である。なお、以下では数式においてベクトル値をとる変数にベクトル記号を付して表記することがある。数式（１５），（１６）のベクトルポテンシャルＡと数式（１５）の電流密度Ｊ0はベクトル値をとる。これに対し、数式（１５），（１６）のスカラポテンシャルφはスカラ値をとる。ただし、本文の説明では、ある変数がスカラ値をとるかベクトル値をとるかを明記しないことがある。
【００６３】
【数１５】
2018136753.tif 000017
【００６４】
【数１６】
2018136753.tif 000018
【００６５】
単位領域の渦電流密度ｉは、数式（１５），（１６）によって決まるベクトルポテンシャルＡとスカラポテンシャルφを用いて数式（１７）のように算出される。渦電流密度ｉはベクトル値をとる。この渦電流密度ｉは、単位領域の中心に定義され、平均電界Ｅに応じて発生する電流密度であると定義される。すなわち、ｉ＝σＥが成立する。
【００６６】
【数１７】
2018136753.tif 000019
【００６７】
次に、粒界層６２のキャパシタンス効果を考慮して上記の基礎方程式を変形する。図５の微細構造のもとでは、結晶粒６１が粒界層６２よりも十分に大きく、結晶粒６１の電流密度ｉ1が単位領域全体の電流密度の大部分を占めている。このため、数式（１３）に現れる結晶粒６１の電流密度ｉ1および導電率σ1を、単位領域を代表する電流密度および導電率とみなし、単位領域の中心に定義されていると考えることができる。そこで、単位領域の渦電流密度ｉは、数式（１８）のように算出される電流密度ｉ1と近似される。電流密度ｉ1はベクトルポテンシャルＡとスカラポテンシャルφと電荷密度ｑに依存する。数式（１７）と異なり、数式（１８）には電荷密度ｑの項が追加されている。電流密度ｉ1および電荷密度ｑは、単位領域の中心の位置に定義される。
【００６８】
【数１８】
2018136753.tif 000020
【００６９】
電流密度ｉ1および導電率σ1を反映させると、数式（１５）の磁場の支配方程式は数式（１９）のように変形できる。数式（１５）と異なり、数式（１９）には電荷密度ｑの項が追加されている。また、数式（１６）の電流密度保存の式は数式（２０）のように変形できる。数式（１６）と異なり、数式（２０）には電荷密度ｑの項が追加されている。
【００７０】
【数１９】
2018136753.tif 000021
【００７１】
【数２０】
2018136753.tif 000022
【００７２】
前述の数式（１４）から、それぞれベクトル値をとる電流密度ｉ1と電荷密度ｑについても数式（２１）が成立する。数式（２１）の時間積分を時間方向の数値積分として表すと数式（２２）が成立する。数式（２２）においてｑnはシミュレーション内の時刻ｎにおける電荷密度であり、ｑn+1はシミュレーション内の時刻ｎ＋１（時刻ｎの次の時刻）における電荷密度である。数式（２２）においてΔｔは、シミュレーション内の時間経過の最小単位であり、時刻ｎと時刻ｎ＋１の差である。数式（２２）のΔｔは、コイル４２の電荷密度Ｊ0の変化を観測する時間周期としての意味をもつ。数式（２２）を変形すると数式（２３）が成立する。時刻ｎ＋１の電荷密度ｑは、時刻ｎ＋１の電流密度ｉ1と直前の時刻である時刻ｎの電荷密度ｑとに依存する。
【００７３】
【数２１】
2018136753.tif 000023
【００７４】
【数２２】
2018136753.tif 000024
【００７５】
【数２３】
2018136753.tif 000025
【００７６】
シミュレーション装置１００は、以上の数式（１８）〜（２０），（２３）を、粒界層６２のキャパシタンス効果を考慮した磁場解析の基礎方程式として使用できる。すなわち、時刻ｎ＋１の電荷密度Ｊ0と時刻ｎの電荷密度ｑから、数式（１９），（２０）に従って、時刻ｎ＋１のベクトルポテンシャルＡと時刻ｎ＋１のスカラポテンシャルφが算出される。次に、時刻ｎ＋１のベクトルポテンシャルＡと時刻ｎ＋１のスカラポテンシャルφと時刻ｎの電荷密度ｑから、数式（１８）に従って、時刻ｎ＋１の電流密度ｉ1が算出される。そして、時刻ｎ＋１の電流密度ｉ1と時刻ｎの電荷密度ｑから、数式（２３）に従って、時刻ｎ＋１の電荷密度ｑが算出される。以上を時間軸に沿って複数の離散時刻について繰り返し実行することで、キャパシタンス効果を考慮した磁場解析が可能となる。
【００７７】
次に、トロイダルコア４１およびコイル４２を含む電気部品の形状を表すモデルについて説明する。また、このモデル上の複数の位置に変数を割り当て、上記の基礎方程式に基づいてこれら複数の変数の解を求める連立方程式を生成する方法について説明する。
【００７８】
図７は、トロイダルコアとコイルの有限要素モデルの例を示す図である。
トロイダルコア４１およびコイル４２を含む電気部品の形状は、モデル７０によって表現される。モデル７０は細分化された領域である要素の集合である。モデル７０は、トロイダルコア４１の領域を示す要素とコイル４２の領域を示す要素とを含む。モデル７０は、トロイダルコア４１およびコイル４２の周辺の空気領域を示す要素を含んでもよい。
【００７９】
各要素は複数の節点および複数の辺によって囲まれた閉領域である。各要素の形状は、例えば、三角形・四角形・四面体・六面体などである。三角形の１つの要素は、３つの節点および３つの辺によって形成される。四面体の１つの要素は、４つの節点および６つの辺によって形成される。同じ節点が２以上の要素によって共有されることがある。また、同じ辺が２以上の要素によって共有されることがある。
【００８０】
図８は、有限要素モデルに対する変数の割り当て例を示す図である。
ここでは説明を簡単にするため、二次元の三角形の要素７１を考える。要素７１は、節点７２，７３，７４および辺７５，７６，７７によって形成されている。辺７５は節点７２と節点７３を結ぶ線分である。辺７６は節点７３と節点７４を結ぶ線分である。辺７７は節点７４と節点７２を結ぶ線分である。
【００８１】
磁場解析に用いるモデル７０が与えられると、シミュレーション装置１００は、モデル７０に含まれる複数の辺それぞれに対して未知のベクトルポテンシャルＡを割り当てる。また、シミュレーション装置１００は、モデル７０に含まれる複数の節点それぞれに対して未知のスカラポテンシャルφを割り当てる。未知のベクトルポテンシャルＡおよび未知のスカラポテンシャルφは、連立方程式の変数に相当する。ここで割り当てるベクトルポテンシャルＡとスカラポテンシャルφは、それぞれスカラ値をとる。
【００８２】
シミュレーションを通じて、ベクトルポテンシャルＡおよびスカラポテンシャルφの近似解が算出される。ベクトル値をとる電流密度ｉ1およびベクトル値をとる電荷密度ｑは、各要素の中心の位置に割り当てられる。前述のように、電流密度ｉ1および電荷密度ｑはベクトルポテンシャルＡおよびスカラポテンシャルφに依存する。ある要素の電流密度ｉ1および電荷密度ｑは、その要素を形成する節点に割り当てられたスカラポテンシャルφとその要素を形成する辺に割り当てられたベクトルポテンシャルＡから算出される。
【００８３】
例えば、節点７２に対して未知のスカラポテンシャルφ1が割り当てられる。節点７３に対して未知のスカラポテンシャルφ2が割り当てられる。節点７４に対して未知のスカラポテンシャルφ3が割り当てられる。また、辺７５に対して未知のベクトルポテンシャルＡ1が割り当てられる。辺７６に対して未知のベクトルポテンシャルＡ2が割り当てられる。辺７７に対して未知のベクトルポテンシャルＡ3が割り当てられる。また、要素７１の中心に対して電流密度ｉ1および電荷密度ｑが割り当てられる。その電流密度ｉ1および電荷密度ｑはＡ1，Ａ2，Ａ3およびφ1，φ2，φ3に依存する。
【００８４】
次に、連立方程式の生成方法を説明する。ある要素のベクトルポテンシャルは数式（２４）のように算出される。数式（２４）において、左辺のＡは要素の中心に定義されるベクトルポテンシャルでありベクトル値をとる。ｎeは１つの要素を形成する辺の数である。要素が二次元の三角形である場合はｎe＝３であり、要素が三次元の四面体である場合はｎe＝６である。Ａiはｉ番目の辺に割り当てられた未知のベクトルポテンシャルでありスカラ値をとる。Ｎeiはｉ番目の辺に対する補間関数でありベクトル的性質をもつ。各要素のベクトルポテンシャルＡは、ｎe個のＡiとｎe個のＮeiから算出される。
【００８５】
【数２４】
2018136753.tif 000026
【００８６】
また、ある要素のスカラポテンシャルは数式（２５）のように算出される。数式（２５）において、左辺のφは要素の中心に定義されるスカラポテンシャルでありスカラ値をとる。ｎnは１つの要素を形成する節点の数である。要素が二次元の三角形である場合はｎn＝３であり、要素が三次元の四面体である場合はｎn＝４である。φiはｉ番目の節点に割り当てられた未知のスカラポテンシャルでありスカラ値をとる。Ｎniはｉ番目の節点に対する補間関数でありスカラ的性質をもつ。各要素のスカラポテンシャルφは、ｎn個のφiとｎn個のＮniから算出される。
【００８７】
【数２５】
2018136753.tif 000027
【００８８】
ここで、基礎方程式の１つである数式（１９）に対して、ベクトルポテンシャルＡに数式（２４）を代入しスカラポテンシャルφに数式（２５）を代入する。そして、数式（１９）の両辺に数式（２６）によって定義される重み関数ｗを乗じて空間積分を行う。数式（２６）においてＮeiは数式（２４）と同様の補間関数であり、ｗeiは辺の重みを表すスカラ値である。上記の計算によって数式（２７）が算出される。数式（２７）からは、ベクトルポテンシャルＡiの数、すなわち、モデル７０に含まれる辺の数に相当する方程式を含む連立方程式が生成される。
【００８９】
【数２６】
2018136753.tif 000028
【００９０】
【数２７】
2018136753.tif 000029
【００９１】
数式（２７）の左辺の第１項は、行列を用いて数式（２８）のように展開することができる。｛ｗe｝は１×辺数のサイズの行列（行ベクトル）、［Ｎe・Ｎe］は辺数×辺数のサイズの行列、｛Ａn+1−Ａn｝は辺数×１のサイズの行列（列ベクトル）である。数式（２７）の左辺の第２項は、行列を用いて数式（２９）のように展開することができる。［Ｎe・∇Ｎn］は辺数×節点数のサイズの行列、｛φn+1｝は節点数×１のサイズの行列である。数式（２７）の左辺の第３項は、行列を用いて数式（３０）のように展開することができる。［Ｎe・ｑ］は辺数×１のサイズの行列である。数式（２７）の左辺の第４項は、行列を用いて数式（３１）のように展開することができる。［∇×Ｎe・∇×Ｎe］は辺数×辺数のサイズの行列、｛Ａn+1｝は辺数×１のサイズの行列である。数式（２７）の右辺は、行列を用いて数式（３２）のように展開することができる。［Ｎe・Ｊ0］は辺数×１のサイズの行列である。
【００９２】
【数２８】
2018136753.tif 000030
【００９３】
【数２９】
2018136753.tif 000031
【００９４】
【数３０】
2018136753.tif 000032
【００９５】
【数３１】
2018136753.tif 000033
【００９６】
【数３２】
2018136753.tif 000034
【００９７】
数式（２７）〜（３２）をまとめると数式（３３）が得られる。数式（３３）からは、モデル７０に含まれる辺の数に相当する方程式を含む連立方程式が生成される。
【００９８】
【数３３】
2018136753.tif 000035
【００９９】
同様に、基礎方程式の１つである数式（２０）に対して、ベクトルポテンシャルＡに数式（２４）を代入しスカラポテンシャルφに数式（２５）を代入する。そして、数式（２０）の両辺に数式（３４）によって定義される重み関数Ｎを乗じて空間積分を行う。数式（３４）においてＮniは数式（２５）と同様の補間関数であり、ｗniは節点の重みを表すスカラ値である。
【０１００】
上記の計算によって算出される数式の左辺の第１項は、行列を用いて数式（３５）のように展開することができる。｛ｗn｝は１×節点数のサイズの行列、［∇Ｎn・Ｎe］は節点数×辺数のサイズの行列、｛Ａn+1−Ａn｝は辺数×１のサイズの行列である。算出される数式の左辺の第２項は、行列を用いて数式（３６）のように展開することができる。［∇Ｎn・∇Ｎn］は節点数×節点数のサイズの行列、｛φn+1｝は節点数×１のサイズの行列である。算出される数式の左辺の第３項は、行列を用いて数式（３７）のように展開することができる。［∇Ｎn・ｑ］は節点数×１のサイズの行列である。
【０１０１】
【数３４】
2018136753.tif 000036
【０１０２】
【数３５】
2018136753.tif 000037
【０１０３】
【数３６】
2018136753.tif 000038
【０１０４】
【数３７】
2018136753.tif 000039
【０１０５】
数式（３５）〜（３７）をまとめると数式（３８）が得られる。数式（３８）からは、モデル７０に含まれる節点の数に相当する方程式を含む連立方程式が生成される。
【０１０６】
【数３８】
2018136753.tif 000040
【０１０７】
上記の数式（３５），（３８）を結合してまとめると、最終的な連立方程式として数式（４５）が得られる。数式（４５）は、係数行列Ｇと解ベクトルＸの積が右辺ベクトルＦに等しいことを表している。係数行列Ｇは、４つの部分係数行列ＧAA，ＧAφ，ＧφA，Ｇφφを結合したものである。部分係数行列ＧAAは、数式（３９）によって定義される辺数×辺数のサイズの行列である。部分係数行列ＧAφは、数式（４０）によって定義される辺数×節点数のサイズの行列である。部分係数行列ＧφAは、数式（４１）によって定義される節点数×辺数のサイズの行列である。部分係数行列Ｇφφは、数式（４２）によって定義される節点数×節点数のサイズの行列である。
【０１０８】
係数行列Ｇの中でＧAAは左上に配置され、ＧAφは右上に配置され、ＧφAは左下に配置され、Ｇφφは右下に配置される。係数行列Ｇは、一辺の長さが辺数＋節点数である正方行列である。解ベクトルＸは、モデル７０の辺に割り当てた未知のベクトルポテンシャルＡとモデル７０の節点に割り当てた未知のスカラポテンシャルφを並べた列ベクトルである。解ベクトルＸは、連立方程式の変数を並べた変数ベクトルに相当する。解ベクトルＸは、行数が辺数＋節点数である列ベクトルである。
【０１０９】
右辺ベクトルＦは、２つの部分ベクトルＦA，Ｆφを結合したものである。部分ベクトルＦAは、数式（４３）によって定義される辺数×１のサイズの列ベクトルである。部分ベクトルＦAは、１つ前の時刻のベクトルポテンシャルＡと電荷密度ｑを用いて算出される。部分ベクトルＦφは、数式（４４）によって定義される節点数×１のサイズの列ベクトルである。部分ベクトルＦφは、１つ前の時刻のベクトルポテンシャルＡと電荷密度ｑを用いて算出される。右辺ベクトルＦの中でＦAは上に配置されＦφは下に配置される。右辺ベクトルＦは、行数が辺数＋節点数である列ベクトルである。
【０１１０】
【数３９】
2018136753.tif 000041
【０１１１】
【数４０】
2018136753.tif 000042
【０１１２】
【数４１】
2018136753.tif 000043
【０１１３】
【数４２】
2018136753.tif 000044
【０１１４】
【数４３】
2018136753.tif 000045
【０１１５】
【数４４】
2018136753.tif 000046
【０１１６】
【数４５】
2018136753.tif 000047
【０１１７】
このように、基礎方程式である数式（１９），（２０）に従って、モデル７０に含まれる辺の数と節点の数の合計、すなわち、未知のベクトルポテンシャルの数と未知のスカラポテンシャルの数の合計に相当する方程式を含む連立方程式が生成される。数式（４５）が示す連立方程式の解は、共役勾配（ＣＧ：Conjugate Gradient）法や不完全コレスキー分解付き共役勾配（ＩＣＣＧ：Incomplete Cholesky CG）法などの反復法で算出できる。
【０１１８】
電流密度ｉ1を示す数式（１８）は、離散時刻のシミュレーションでは数式（４６）のように変形することができる。時刻ｎ＋１のベクトルポテンシャルＡn+1、時刻ｎ＋１のスカラポテンシャルφn+1、時刻ｎのベクトルポテンシャルＡnおよび時刻ｎの電荷密度ｑnから、各要素の時刻ｎ＋１の電流密度ｉ1n+1が数式（４６）に従って算出される。そして、時刻ｎ＋１の電流密度ｉ1n+1および時刻ｎの電荷密度ｑnから、各要素の時刻ｎ＋１の電荷密度ｑn+1が前述の数式（２３）に従って算出される。
【０１１９】
【数４６】
2018136753.tif 000048
【０１２０】
次に、シミュレーション装置１００の機能および処理手順について説明する。
図９は、シミュレーション装置の機能例を示すブロック図である。
シミュレーション装置１００は、モデル記憶部１１１、パラメータ記憶部１１２、結果記憶部１１３、方程式生成部１２１、解算出部１２２および結果表示部１２３を有する。モデル記憶部１１１、パラメータ記憶部１１２および結果記憶部１１３は、例えば、ＲＡＭ１０２またはＨＤＤ１０３に確保した記憶領域を用いて実装される。方程式生成部１２１、解算出部１２２および結果表示部１２３は、例えば、ＣＰＵ１０１が実行するプログラムモジュールを用いて実装される。
【０１２１】
モデル記憶部１１１は、モデル７０に含まれる複数の要素を示すモデルデータを記憶する。モデルデータは、ユーザの操作に基づいて生成されてもよいしモデリングソフトウェアなどのソフトウェアによって自動生成されてもよい。また、モデルデータは、シミュレーション装置１００が生成してもよいし他の装置が生成してもよい。
【０１２２】
パラメータ記憶部１１２は、シミュレーションに用いる各種のパラメータを記憶する。パラメータは、例えば、ユーザによって指定される。パラメータには、コイル４２を流れる電流の周波数などの稼働条件を示すパラメータが含まれる。また、パラメータには、磁性材料の導電率などの物理的特性を示す物性値が含まれる。また、パラメータには、シミュレーションの時間幅などのシミュレーション条件を示すパラメータが含まれる。
【０１２３】
結果記憶部１１３は、シミュレーション結果を記憶する。シミュレーション結果には、ベクトルポテンシャル、スカラポテンシャル、電流密度および電荷密度が含まれる。
方程式生成部１２１は、モデル記憶部１１１に記憶されたモデルデータに基づいて、モデル７０上の複数の位置に変数を割り当てる。第２の実施の形態では、方程式生成部１２１は、要素の辺にベクトルポテンシャルの変数を割り当て、要素の節点にスカラポテンシャルの変数を割り当てる。そして、方程式生成部１２１は、パラメータ記憶部１１２に記憶されたパラメータと結果記憶部１１３に記憶された１つ前のステップのシミュレーション結果を用いて、変数の解を求める連立方程式を生成する。方程式生成部１２１は、パラメータに従って、時間軸に沿って連立方程式の生成を繰り返す。
【０１２４】
解算出部１２２は、方程式生成部１２１が生成した連立方程式を反復法によって解き、変数の近似解を求める。第２の実施の形態では、解算出部１２２は、要素の辺のベクトルポテンシャルの近似解と要素の節点のスカラポテンシャルの近似解を求める。また、解算出部１２２は、最新のベクトルポテンシャルおよびスカラポテンシャルと結果記憶部１１３に記憶された１つ前のステップのシミュレーション結果を用いて、各要素の電流密度と電荷密度を算出する。解算出部１２２は、最新のステップのベクトルポテンシャル、スカラポテンシャル、電流密度および電荷密度を、結果記憶部１１３に保存する。
【０１２５】
結果表示部１２３は、方程式生成部１２１および解算出部１２２による反復処理が終了すると、結果記憶部１１３に記憶されたシミュレーション結果をディスプレイ３１に表示させる。このとき、結果表示部１２３は、シミュレーション結果に含まれる数値をディスプレイ３１に表示させてもよい。また、結果表示部１２３は、シミュレーション結果から変換された可視化画像をディスプレイ３１に表示させてもよい。また、結果表示部１２３は、ディスプレイ３１以外の出力装置にシミュレーション結果を出力してもよい。
【０１２６】
図１０は、節点テーブルと要素テーブルの例を示す図である。
節点テーブル１１４は、モデル記憶部１１１に記憶されている。節点テーブル１１４は、節点番号、Ｘ座標およびＹ座標の項目を含む。節点番号は、節点の識別番号である。Ｘ座標およびＹ座標は、直交座標系における節点の位置を示す座標である。なお、モデル７０が三次元モデルである場合、節点テーブル１１４に更にＺ座標の項目が含まれる。
【０１２７】
要素テーブル１１５は、モデル記憶部１１１に記憶されている。要素テーブル１１５は、要素番号、節点１、節点２および節点３の項目を含む。要素番号は、要素の識別番号である。節点１、節点２および節点３は、要素を形成する節点の節点番号を示す。要素を形成する辺は、要素を形成する２つの節点の組によって識別することができる。なお、要素が四面体である場合、要素テーブル１１５に更に節点４の項目が含まれる。すなわち、１つの要素番号に対して、要素を形成する節点の数だけ節点番号が対応付けられる。
【０１２８】
図１１は、パラメータテーブルの例を示す図である。
パラメータテーブル１１６は、パラメータ記憶部１１２に記憶されている。パラメータテーブル１１６は、パラメータ名とその値を対応付ける。
【０１２９】
パラメータには、総ステップ数および時間幅Δｔが含まれる。総ステップ数は、シミュレーション上の離散的な時刻の数、すなわち、連立方程式を生成して解を求める処理の繰り返し回数を表す。時間幅Δｔは、あるステップと次のステップの間のシミュレーション上の時間差、すなわち、経過時間の最小単位を表す。
【０１３０】
また、パラメータには、コイル電流およびコイル電流周波数が含まれる。コイル電流は、コイル４２を流れる電流の振幅を表す。コイル電流周波数は、コイル４２を流れる電流の周波数を表す。コイル４２を流れる電流が正弦的に変化すると仮定して、コイル電流およびコイル電流周波数から、各時刻におけるコイル４２の電流密度Ｊ0を算出できる。ただし、パラメータテーブル１１６に、各時刻の電流密度Ｊ0が記録されていてもよい。
【０１３１】
また、パラメータには、複数の材料それぞれの物性値が含まれる。材料１は、トロイダルコア４１およびコイル４２の周囲に存在する空気である。空気に関するパラメータには透磁率μが含まれる。材料２はコイル４２の材料である。コイル４２に関するパラメータには、透磁率μおよび導電率σが含まれる。材料３は、トロイダルコア４１に用いる磁性材料である。磁性材料に関するパラメータには、透磁率μ、導電率σ1、比誘電率εrおよび寸法比ｒが含まれる。ここでは３つの材料のパラメータを記載しているが、パラメータテーブル１１６には解析対象に応じて１以上の材料のパラメータが記載される。
【０１３２】
図１２は、結果テーブルの例を示す図である。
結果テーブル１１７は、結果記憶部１１３に記憶されている。結果テーブル１１７は、ステップ番号、ベクトルポテンシャル、スカラポテンシャル、電流密度および電荷密度の項目を含む。ステップ番号は、シミュレーションの各ステップ（各時刻）を識別する番号であり、０以上かつ総ステップ数以下の整数である。
【０１３３】
ベクトルポテンシャルの項目には、あるステップで算出された複数の辺のベクトルポテンシャルが列挙される。ただし、ステップ０における各辺のベクトルポテンシャルは０に初期化される。スカラポテンシャルの項目には、あるステップで算出された複数の節点のスカラポテンシャルが列挙される。ただし、ステップ０におけるスカラポテンシャルは算出しなくてよい。電流密度の項目には、あるステップで算出された複数の要素の電流密度が列挙される。ただし、ステップ０における電流密度は算出しなくてよい。電荷密度の項目には、あるステップで算出された複数の要素の電荷密度が列挙される。ただし、ステップ０における各要素の電荷密度は０に初期化される。
【０１３４】
図１３は、シミュレーションの手順例を示すフローチャートである。
（Ｓ１０）方程式生成部１２１は、モデル記憶部１１１からモデルデータを読み出す。また、方程式生成部１２１は、パラメータ記憶部１１２からパラメータを読み出す。
【０１３５】
（Ｓ１１）方程式生成部１２１は、モデルデータが示す複数の辺それぞれにベクトルポテンシャル変数（未知のベクトルポテンシャルＡを示す変数）を割り当てる。また、方程式生成部１２１は、モデルデータが示す複数の節点それぞれにスカラポテンシャル変数（未知のスカラポテンシャルφを示す変数）を割り当てる。
【０１３６】
（Ｓ１２）方程式生成部１２１は、ステップ番号ｎを０に初期化する。また、方程式生成部１２１は、パラメータが示す総ステップ数（ステップ数の上限）を特定する。
（Ｓ１３）解算出部１２２は、ステップＳ１１で割り当てられたベクトルポテンシャル変数のステップ０における値（ベクトルポテンシャルＡ0）を０に初期化する。また、解算出部１２２は、モデルデータが示す複数の要素それぞれのステップ０における電荷密度ｑ0を０に初期化する。解算出部１２２は、結果記憶部１１３に記憶された結果テーブル１１７に、ベクトルポテンシャルＡ0と電荷密度ｑ0を登録する。
【０１３７】
（Ｓ１４）方程式生成部１２１は、ステップＳ１０で読み出したパラメータに基づいて、ステップｎ＋１におけるコイル４２の電流密度Ｊ0を算出する。
（Ｓ１５）方程式生成部１２１は、数式（４５）に従って、ステップｎ＋１に対応する連立方程式を示す係数行列Ｇと右辺ベクトルＦを生成する。係数行列Ｇの生成にあたっては、ステップＳ１０で読み出したパラメータが用いられる。右辺ベクトルＦの生成にあたっては、ステップＳ１４で算出した電流密度Ｊ0と、ステップＳ１０で読み出したパラメータが用いられる。更に、右辺ベクトルＦの生成にあたっては、結果テーブル１１７に記憶されたステップｎのベクトルポテンシャルＡnと電荷密度ｑnが用いられる。
【０１３８】
（Ｓ１６）解算出部１２２は、ステップＳ１５で生成された連立方程式の解をＣＧ法やＩＣＣＧ法などの反復法によって算出する。これにより、ステップｎ＋１のベクトルポテンシャルＡn+1とスカラポテンシャルφn+1が算出される。解算出部１２２は、結果テーブル１１７にベクトルポテンシャルＡn+1とスカラポテンシャルφn+1を登録する。
【０１３９】
（Ｓ１７）解算出部１２２は、数式（４６）に従って、ステップｎ＋１における各要素の電流密度ｉ1n+1を算出する。電流密度ｉ1n+1の算出にあたっては、読み出したパラメータ、ステップｎ＋１のベクトルポテンシャルＡn+1、ステップｎのベクトルポテンシャルＡn、ステップｎ＋１のスカラポテンシャルφn+1およびステップｎの電荷密度ｑnが用いられる。解算出部１２２は、結果テーブル１１７に電流密度ｉ1n+1を登録する。
【０１４０】
（Ｓ１８）解算出部１２２は、数式（２３）に従って、ステップｎ＋１における各要素の電荷密度ｑn+1を算出する。電荷密度ｑn+1の算出にあたっては、読み出したパラメータ、ステップｎ＋１の電流密度ｉ1n+1およびステップｎの電荷密度ｑnが用いられる。解算出部１２２は、結果テーブル１１７に電荷密度ｑn+1を登録する。
【０１４１】
（Ｓ１９）方程式生成部１２１は、現在のステップ番号ｎが総ステップ数未満であるか判断する。ステップ番号ｎが総ステップ数未満である場合にはステップＳ２０に処理が進み、ステップ番号ｎが総ステップ数に達した場合にはシミュレーションが終了する。
【０１４２】
（Ｓ２０）方程式生成部１２１は、ステップ番号ｎに１を加えステップＳ１４に進む。
次に、シミュレーション結果の例を説明する。
図１４は、寸法共鳴現象を考慮しないシミュレーション例を示すグラフである。
【０１４３】
粒界層６２のキャパシタンス効果を考慮せず寸法共鳴現象が再現されないシミュレーションを行った場合、グラフ５２のようなシミュレーション結果が得られる。キャパシタンス効果を考慮しないシミュレーションは、前述の数式（１５）〜（１７）を基礎方程式として使用したシミュレーションである。グラフ５２は、前述のグラフ５１と同様、コイル４２を流れる電流の周波数とトロイダルコア４１の比透磁率との間の関係を示している。グラフ５１が示すトロイダルコア４１の実際の比透磁率と異なり、グラフ５２は比透磁率の実部の増大や虚部の増大を表現しておらず寸法共鳴現象を再現していない。
【０１４４】
図１５は、寸法共鳴現象を考慮したシミュレーション例を示すグラフである。
粒界層６２のキャパシタンス効果を考慮して寸法共鳴現象が再現されるシミュレーションを行った場合、グラフ５３のようなシミュレーション結果が得られる。キャパシタンス効果を考慮したシミュレーションは、前述の数式（１８）〜（２０），（２３）を基礎方程式として使用したシミュレーションである。グラフ５３は、前述のグラフ５１，５２と同様、コイル４２を流れる電流の周波数とトロイダルコア４１の比透磁率との間の関係を示している。グラフ５１が示す実際の比透磁率と同様、グラフ５３は比透磁率の実部の増大や虚部の増大を表現しており寸法共鳴現象を再現している。また、グラフ５３において比透磁率の実部が最大となるピーク周波数や虚部が最大となるピーク周波数は、グラフ５１が示す実際のピーク周波数と十分に近似している。
【０１４５】
図１６は、磁束密度可視化画像の例を示す図である。
前述のように、結果表示部１２３は、結果テーブル１１７に登録された結果データを可視化画像に変換してディスプレイ３１に表示することが可能である。結果表示部１２３が生成可能な可視化画像の１つとして、図１６のような磁束密度可視化画像５４が挙げられる。磁束密度可視化画像５４は、トロイダルコア４１の内部の磁束密度の大きさおよび方向を矢印を用いて表現している。磁束密度はベクトルポテンシャルＡから算出できる。
【０１４６】
図１７は、電流密度可視化画像の例を示す図である。
結果表示部１２３が生成可能な可視化画像の他の１つとして、図１７のような電流密度可視化画像５５が挙げられる。電流密度可視化画像５５は、トロイダルコア４１の内部の渦電流密度ｉの大きさおよび方向を矢印を用いて表現している。
【０１４７】
第２の実施の形態のシミュレーション装置１００によれば、有限要素法を用いた磁性材料シミュレーションにおいて、粒界層６２のキャパシタンス効果が各要素に蓄えられる電荷の電荷密度ｑとして表現され、磁場に関する基礎方程式に電荷密度ｑが追加される。時刻ｎ＋１におけるベクトルポテンシャルＡn+1およびスカラポテンシャルφn+1が、時刻ｎのベクトルポテンシャルＡnおよび電荷密度ｑnから算出される。時刻ｎ＋１における各要素の電流密度ｉ1n（渦電流密度）が、時刻ｎ＋１のベクトルポテンシャルＡn+1およびスカラポテンシャルφn+1と時刻ｎのベクトルポテンシャルＡnおよび電荷密度ｑnから算出される。時刻ｎ＋１における各要素の電荷密度ｑn+1が、時刻ｎ＋１の電流密度ｉ1nと時刻ｎの電荷密度ｑnから算出され、時刻ｎ＋２の計算に使用される。
【０１４８】
これにより、多結晶質の磁性材料で発生する寸法共鳴現象を再現することができ、磁性材料シミュレーションの精度を向上させることができる。
【符号の説明】
【０１４９】
１０  磁性材料シミュレーション装置
１１  記憶部
１２  処理部
１３  形状モデル
１４  特性情報
１５  コイル電流情報
１６ａ，１８ａ  電流密度
１６ｂ，１６ｃ，１６ｄ，１８ｂ，１８ｃ，１８ｄ  指標値
１７  電荷密度
【図１】
2018136753.tif 000049
【図２】
2018136753.tif 000050
【図３】
2018136753.tif 000051
【図４】
2018136753.tif 000052
【図５】
2018136753.tif 000053
【図６】
2018136753.tif 000054
【図７】
2018136753.tif 000055
【図８】
2018136753.tif 000056
【図９】
2018136753.tif 000057
【図１０】
2018136753.tif 000058
【図１１】
2018136753.tif 000059
【図１２】
2018136753.tif 000060
【図１３】
2018136753.tif 000061
【図１４】
2018136753.tif 000062
【図１５】
2018136753.tif 000063
【図１６】
2018136753.tif 000064
【図１７】
2018136753.tif 000065