(57)【要約】      （修正有）【課題】磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる磁場解析装置を提供する。【解決手段】磁場解析装置１は、解析条件と磁場の運動方程式を記憶する記憶手段と、前記解析条件から、前記磁性体を構成する粒子の磁場の運動方程式の係数および変数の初期値および定数を演算して、磁場の運動方程式の解を演算し、磁場の運動方程式の解に基づき、前記空間における任意の点の磁場を演算する演算手段とを有している。磁場解析装置１は、磁場の運動方程式と前記入力情報に基づき、演算手段を用いて空間における任意の磁場を演算する。【選択図】図１
2010160782.tif 000002
【特許請求の範囲】
【請求項１】
磁性体が配置された空間における任意の点の磁場を情報処理装置を使って解析する磁場解析装置であって、
解析条件と磁場の運動方程式を記憶する記憶手段と、
前記解析条件から、前記磁性体を構成する粒子の磁場の運動方程式の係数、変数の初期値および定数を演算して、前記磁場の運動方程式の解を演算し、該解に基づき、前記空間における任意の点の磁場を演算する演算手段と、
を有することを特徴とする磁場解析装置。
【請求項２】
前記解析条件は、
前記磁性体の形状と、
前記磁性体の磁化曲線を表す関数と、
を含むことを特徴とする請求項１記載の磁場解析装置。
【請求項３】
前記磁場の運動方程式の変数は、前記磁性体を構成する粒子の磁場ベクトルであることを特徴とする請求項１記載の磁場解析装置。
【請求項４】
前記磁場の運動方程式の係数は、
前記磁性体を構成する粒子の位置ベクトルと、
前記磁性体の磁化曲線を表す関数と、
から演算されることを特徴とする請求項１記載の磁場解析装置。
【請求項５】
前記磁場の運動方程式は、
磁場の項を含む粒子のラグランジアンを磁場ベクトルと磁場ベクトルの時間微分を正準変数としてラグランジュの運動方程式に代入することで求められたことを特徴とする請求項１〜４のいずれかに記載の磁場解析装置。
【請求項６】
前記空間には導体がさらに配置されており、
前記解析条件は導体の形状と電流を有することを特徴とする請求項１〜５のいずれかに記載の磁場解析装置。
【請求項７】
前記粒子のラグランジアンは、前記磁性体を構成する粒子をＮ個としたとき以下の式（Ａ）〜（Ｃ）で表される関数であることを特徴とする請求項５記載の磁場解析装置。
【数１】
2010160782.tif 000024
α：仮想質量
（太字の）ｒ：位置ベクトル
（太字の）Ｈ：磁場ベクトル
（太字の傍点付き）Ｈ：磁場ベクトルの時間微分
（太字の）Ｍ：磁化ベクトル
（太字の）ｎ：粒子の位置ベクトルを重心とする多面体要素境界面の法線ベクトル
（太字の）Ｈｅｘｔ：外部からの印加磁場ベクトル
ΔＳ：粒子の位置ベクトルを重心とする多面体要素境界面の面積
ΔＶ：粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積
χ：磁気感受率
μ０：真空の透磁率
λ：ラグランジュの未定定数
π：円周率
ｓ：粒子の位置ベクトルを重心とする多面体要素境界面数
添え字ｉｐ:粒子の位置ベクトルを重心とする多面体要素内のp点での物理量
添え字ｊｑ:粒子の位置ベクトルを重心とする多面体要素内のｑ点での物理量
ｍ：質量
ｖ：速度
φ（ｒｉ−ｒｊ）：ｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギ
添え字ｉ、ｊ：ｉ番目の粒子とｊ番目の粒子の物理量
【請求項８】
前記磁場の運動方程式は、以下の式（Ｄ）で表される式であることを特徴とする請求項７記載の磁場解析装置。
【数２】
2010160782.tif 000025
【請求項９】
前記磁場の運動方程式は、以下の式（Ｅ）、（Ｆ）、（Ｇ）で表される離散化された式と、以下の式（Ｈ）、（Ｉ）で表される束縛力を加えた式とを有することを特徴とする請求項７記載の磁場解析装置。
【数３】
2010160782.tif 000026
γ：減衰定数
δｔ：仮想的な時間刻み
添え字ｎ：ｎδｔにおける物理量
添え字ｎ−１／２：（ｎ−１／２）δｔにおける物理量
添え字ｎ＋１／２：（ｎ＋１／２）δｔにおける物理量
添え字ｎ＋１：（ｎ＋１）δｔにおける物理量
【請求項１０】
前記磁性体はモータのロータおよびステータであることを特徴とする請求項１〜９のいずれかに記載の磁場解析装置。
【請求項１１】
前記導体はモータのコイルであることを特徴とする請求項６に記載の磁場解析装置。
【請求項１２】
前記磁性体はＰＭモータ（Permanent Magnet Motor）のロータおよびステータであることを特徴とする請求項１〜９のいずれかに記載の磁場解析装置。
【請求項１３】
前記導体はＰＭモータ（Permanent Magnet Motor）のコイルであることを特徴とする請求項６に記載の磁場解析装置。
【請求項１４】
前記磁性体は磁気シールドを構成する磁性体であることを特徴とする請求項１〜９のいずれかに記載の磁場解析装置。
【請求項１５】
コンピュータを請求項１から１４のいずれかに記載の磁場解析装置として機能させるためのプログラム。
【請求項１６】
磁性体が配置された空間における任意の点の磁場を情報処理装置を使って解析する磁場解析方法であって、
解析条件と磁場の運動方程式を記憶する工程（ａ）と、
前記解析条件から、前記磁性体を構成する粒子の磁場の運動方程式の係数および変数の初期値および定数を演算して、磁場の運動方程式の解を演算し、
磁場の運動方程式の解に基づき、前記空間における任意の点の磁場を演算する工程（ｂ）と、
を有することを特徴とする磁場解析方法。
【請求項１７】
前記解析条件は、
前記磁性体の形状と
前記磁性体の磁化曲線を表す関数と、
を有することを特徴とする請求項１６記載の磁場解析方法。
【請求項１８】
前記磁場の運動方程式の変数は、前記磁性体を構成する粒子の磁場ベクトルであることを特徴とする請求項１６記載の磁場解析方法。
【請求項１９】
前記磁場の運動方程式の係数は、
前記磁性体を構成する粒子の位置ベクトルと、
前記磁性体の磁化曲線を表す関数と、
から演算されることを特徴とする請求項１６記載の磁場解析方法。
【請求項２０】
前記運動方程式は、
磁場の項を含む粒子のラグランジアンを磁場ベクトルと磁場ベクトルの時間微分を正準変数としてラグランジュの運動方程式に代入することで求められたことを特徴とする請求項１６〜１９のいずれかに記載の磁場解析方法。
【請求項２１】
前記空間には導体がさらに配置されており、
前記解析条件は導体の形状と電流を有することを特徴とする請求項１６〜２０のいずれかに記載の磁場解析方法。
【請求項２２】
前記粒子のラグランジアンは、磁性体を構成する粒子をＮ個としたとき、以下の式（Ａ）〜（Ｃ）で表される関数であることを特徴とする請求項２０記載の磁場解析方法。
【数４】
2010160782.tif 000027
α：仮想質量
（太字の）ｒ：位置ベクトル
（太字の）Ｈ：磁場ベクトル
（太字の傍点付き）Ｈ：磁場ベクトルの時間微分
（太字の）Ｍ：磁化ベクトル
（太字の）ｎ：粒子の位置ベクトルを重心とする多面体要素境界面の法線ベクトル
（太字の）Ｈｅｘｔ：外部からの印加磁場ベクトル
ΔＳ：粒子の位置ベクトルを重心とする多面体要素境界面の面積
ΔＶ：粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積
χ：磁気感受率
μ０：真空の透磁率
λ：ラグランジュの未定定数
π：円周率
ｓ：粒子の位置ベクトルを重心とする多面体要素境界面数
添え字ｉｐ:粒子の位置ベクトルを重心とする多面体要素内のｐ点での物理量
添え字ｊｑ:粒子の位置ベクトルを重心とする多面体要素内のｑ点での物理量、
ｍ：質量
ｖ：速度
φ（ｒｉ−ｒｊ）：ｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギ
添え字ｉ、ｊ：ｉ番目の粒子とｊ番目の粒子の物理量
【請求項２３】
前記磁場の運動方程式は、以下の式（Ｄ）で表される式であることを特徴とする請求項２２記載の磁場解析方法。
【数５】
2010160782.tif 000028
【請求項２４】
前記磁場の運動方程式は、以下の式（Ｅ）、（Ｆ）、（Ｇ）で表される離散化された式と、以下の式（Ｈ）、（Ｉ）で表される束縛力を加えた式とを有することを特徴とする請求項２３記載の磁場解析方法。
【数６】
2010160782.tif 000029
γ：減衰定数
δｔ：仮想的な時間刻み
添え字ｎ：ｎδｔにおける物理量、
添え字ｎ−１／２：（ｎ−１／２）δｔにおける物理量
添え字ｎ＋１／２：（ｎ＋１／２）δｔにおける物理量
添え字ｎ＋１：（ｎ＋１）δｔにおける物理量
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、磁場解析装置および磁場解析方法に関する。
【背景技術】
【０００２】
近年、モータ等の磁気エネルギーを駆動源とする電気機器は、高精度化、高効率化が要求されると同時に、開発スピードのアップが求められている。
【０００３】
このような要求を満たすためには電気機器の磁場解析による検証は不可欠となっている。
【０００４】
磁場解析の手法としては従来、有限要素法が広く用いられてきた。
【０００５】
例えば特許文献１では、空気中に磁気コアが存在する場合の磁場解析例として、有限要素法を用いた解析方法が記載されている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開平１０−１８８０４２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、有限要素法によって磁場解析を行う場合は、対象とする電気機器を含んだ空間全域をメッシュ分割する必要がある。
【０００８】
そのため、大規模、複雑形状の電気機器の解析や、積層鉄心等の高アスペクト比問題を取り扱う場合、解析に膨大な時間と解析機器のメモリを要し、十分な精度で効率よく磁場解析を行うのが困難であるという問題があった。
【０００９】
本発明は、このような問題に鑑みてなされたもので、その目的は、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる磁場解析装置を提供することにある。
【課題を解決するための手段】
【００１０】
前述した目的を達成するために、本発明者は鋭意検討の結果、情報処理装置を使って磁性体を構成する粒子の磁場の運動方程式の解を演算し、磁場の運動方程式の解に基づいて磁性体が配置された空間の任意の場所の磁場を演算することにより、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができることを見出し、本発明をするに至った。
【００１１】
即ち、第１の発明は、磁性体が配置された空間における任意の点の磁場を情報処理装置を使って解析する磁場解析装置であって、解析条件と磁場の運動方程式を記憶する記憶手段と、前記解析条件から、前記磁性体を構成する粒子の磁場の運動方程式の係数、変数の初期値および定数を演算して、磁場の運動方程式の解を演算し、その解に基づき、前記空間における任意の点の磁場を演算する演算手段と、を有することを特徴とする磁場解析装置である。
【００１２】
第２の発明は、コンピュータを第１の発明に記載の磁場解析装置として機能させるためのプログラムである。
【００１３】
第３の発明は、磁性体が配置された空間における任意の点の磁場を情報処理装置を使って解析する磁場解析方法であって、解析条件と磁場の運動方程式を記憶する工程（ａ）と、前記解析条件から、前記磁性体を構成する粒子の磁場の運動方程式の係数、変数の初期値および定数を演算して、磁場の運動方程式の解を演算し、磁場の運動方程式の解に基づき、前記空間における任意の点の磁場を演算する工程（ｂ）と、を有することを特徴とする磁場解析方法である。
【発明の効果】
【００１４】
本発明によれば、磁性体を有する系に対して、解析対象を含む空間全域をメッシュ分割する必要なく、十分な精度で効率よく磁場解析を行うことができる磁場解析装置を提供することができる。
【００１５】
また、磁性体を構成する粒子のラグランジアンから導出された磁場の運動方程式の解を演算すれば、外部からの入力が変化しない場合、系の全エネルギーが保存される。
【図面の簡単な説明】
【００１６】
【図１】磁場解析装置１のハードウェア構成を示す図である。
【図２】記憶装置５を示す図である。
【図３】ＳＰＭモータ３１の構成の概略を示す図である。
【図４】図３のステータティース４３の１つの周囲の拡大斜視図である。
【図５】磁場解析装置１を用いた磁場の解析の手順を示すフローチャートである。
【図６】図６（ａ）はコイル４５の拡大斜視図であって、図６（ｂ）はコイル４５を、ローカル導体に分割した例を示す図である。
【図７】直方体導体４５ａとローカル座標の関係を示す図である。
【図８】円弧状柱状導体４５ｃとローカル座標の関係を示す図である。
【図９】磁性体の磁場計算に用いる要素の説明図である。
【発明を実施するための形態】
【００１７】
以下、図面に基づいて本発明に好適な実施形態を詳細に説明する。
【００１８】
まず、図１を参照して、本実施形態に係る磁場解析装置１のハードウェア構成を説明する。
【００１９】
なお、図１のハードウェア構成は例示であり、これに限定されないのは当然である。
【００２０】
図１に示すように、磁場解析装置１は制御部３、記憶装置５、メディア入出力部６、入力部７、表示部９、プリンタポート１１等がバス１３を介して互いに接続されている。
【００２１】
制御部３は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等で構成され、記憶部としての記憶装置５に格納されたプログラムに従って、バス１３を介して接続された各装置を駆動制御する。
【００２２】
図２に示すように、記憶装置５には、磁場解析装置１の各構成部分を駆動制御するための制御プログラム１５、本発明を実施するための磁場解析プログラム１７が格納されている。
【００２３】
磁場解析プログラム１７は、解析条件を有する情報である入力情報２１と、入力情報２１に基づき、磁場の運動方程式に基づいて磁場を演算する演算プログラム１９とを有している。
【００２４】
メディア入出力部６は、フロッピー（登録商標）ディスク、ＣＤ、ＤＶＤ等のメディアとの間で情報の入出力を行う装置である。
【００２５】
入力部７は、キーボード、マウス等の入力装置であり、表示部９はディスプレイ等の表示機器である。
【００２６】
プリンタポート１１には出力装置としてのプリンタ１２等が接続される。
【００２７】
次に、磁場解析装置１を用いた磁場の解析の手順について図３〜９を参照して説明する。
【００２８】
ここではＰＭモータ（Permanent Magnet Motor）の一種である、ＳＰＭモータ３１（Surface Permanent Magnet Motor）の磁場解析を例にして説明する。
【００２９】
まず、ＳＰＭモータ３１の構成の概略を図３および図４を参照して説明する。
【００３０】
図３に示すように、ＳＰＭモータ３１は回転子（移動子）であるロータ３３と固定子であるステータ３５を有している。
【００３１】
ロータ３３は鉄等の磁性体である円柱状のロータコア３７を有し、ロータコア３７の表面には永久磁石３９が設けられている。
【００３２】
ロータコア３７の軸中心には棒状のロータシャフト４１が設けられている。
【００３３】
ステータ３５は磁性体である歯状のステータティース４３とステータティース４３の外側に設けられた円筒状の磁性体であるコアバック４４、コアバック４４の外側に設けられた円筒状のフレーム４６から構成されている。
【００３４】
図３および図４に示すように、ステータティース４３には、金属等の導電体であるコイル４５が巻きつけられている。
【００３５】
なお、実際のＳＰＭモータ３１ではコイル４５は仕様に応じたターン数でステータティース４３に巻きつけられて束となっているが、本実施形態では、図３および図４に描かれているように、コイル一本一本をモデル化せず、コイルの束を一つの導体として扱う。
【００３６】
このような構造のＳＰＭモータ３１は、永久磁石３９の磁場、およびコイル４５に電流を流すことにより発生する磁場によって、ロータ３３、ステータ３５が磁化する。磁性体の磁気エネルギーの偏差によりＳＰＭモータ３１は駆動する。
【００３７】
そのため、ＳＰＭモータ３１の磁場解析を行うためにはコイル４５、永久磁石３９がロータ３３、ステータ３５を構成する磁性体上に作る磁場ベクトルを計算し、これら磁性体の磁化現象を解析する必要がある。
【００３８】
次に、解析の手順について図５〜図９を参照して説明する。
【００３９】
なお、以下の手順においてはＳＰＭモータ３１を、複数に要素分割して要素ごとの粒子の集合体とし、剛体モデルとして扱っているが、本発明これに限定されることはなく、粒子を用いずに要素ごとの剛体モデルとして扱ってもよい。
【００４０】
また、以下の手順において、磁性体とは、ロータ３３、ステータ３５を構成する磁性体と永久磁石３９を指し示し、磁化曲線を表す関数によってこれらは区別される。
【００４１】
まず、磁場解析装置１の制御部３は磁場解析プログラム１７を起動し、解析したいＳＰＭモータ３１の解析条件としての三次元構造（形状、座標点）、質量密度、磁性体の磁化曲線を表す関数、導体の電流密度ベクトルを記憶装置５の入力情報２１として記憶する（図５のステップ１０１）。
【００４２】
これらの物理量は例えばメディア入出力部６を介してＣＤ−ＲＯＭ等の記録媒体から読込んだものであってもよい。
【００４３】
また、ＳＰＭモータ３１の三次元構造の情報とは例えばＣＡＤ等のデータである。
【００４４】
さらに、あらかじめ上記物理量が入力情報２１として記憶されている場合は、上記ステップは不要である。
【００４５】
電流密度ベクトルは、コイル４５の作る磁場ベクトルを計算する際に必要になる。
【００４６】
さらに、磁性体の磁化曲線を表す関数は磁化ベクトル（太字の）Ｍを計算する際に必要になる。
【００４７】
以上が図５のステップ１０１の詳細である。
【００４８】
次に、磁場解析装置１の制御部３は、入力情報２１の有する三次元構造の情報から、磁性体をＮ個の粒子（粒子の位置ベクトルを重心とする多面体要素（本実施例では立方体要素）に分割し、粒子の位置ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０２）。
【００４９】
ここで、Ｎは任意の整数であり、粒子の数Ｎおよび位置ベクトルは入力情報２１の有する三次元構造およびあらかじめ記憶装置５に記憶されている多面体要素の形状により計算される。
【００５０】
以上が図５のステップ１０２の詳細である。
【００５１】
次に、制御部３は、入力情報２１の有する三次元構造の情報から、導体（図６（ａ）に示すコイル４５）を図６（ｂ）に示すようにローカル導体（直方体導体４５ａ、４５ｂと円弧状柱状導体４５ｃ、４５ｄ）に分割し、それぞれの導体が作る磁場ベクトルを計算するための係数を計算し、記憶装置５に記憶する（図５のステップ１０３）。
【００５２】
ここで、図５のステップ１０３について、具体的に説明する。
【００５３】
まず、ローカル導体が直方体導体４５ａの場合について説明する。
【００５４】
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。
【００５５】
ローカル導体が直方体導体４５ａの場合は、制御部３は、図７に示すように、ローカル座標系（ｘｓ,ｙｓ,ｚｓ）を適用する。
【００５６】
このローカル座標系においてはローカル導体（直方体導体４５ａ）の重心を原点Ｏｓとし、直方体導体４５ａの寸法はｘｓ方向に２ａ、ｙｓ方向に２ｂ、ｚｓ方向に２ｃの長さを持つものとする。
【００５７】
また原点Ｏｓに粒子は位置するものとする。
【００５８】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体（直方体導体４５ａ）の寸法であるａ，ｂ，ｃと粒子位置ベクトルを計算し記憶装置５に記憶する。
【００５９】
次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
【００６０】
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。
【００６１】
ローカル導体が円弧状柱状導体４５ｃの場合は、制御部３は、図８に示すようにローカル座標系（ｘｃ,ｙｃ,ｚｃ）を適用する。
【００６２】
このローカル座標系においては原点Ｏｃは円弧の中心軸上に存在し、かつ円弧状柱状導体４５ｃの高さ方向（図８のｚｃ方向）に対して円弧状柱状導体４５ｃが対称となる点に存在するものとする。
【００６３】
また、ｘｃ,ｙｃ,ｚｃは、ｘｃ−ｙｃ平面でみると、＋ｘｃ軸を基点とし、円弧状柱状導体４５ｃの円弧が＋ｚｃ軸からみて反時計回りになるようして決定する。
【００６４】
円弧状柱状導体４５ｃの内径と外径の平均値をＲｃとし、径方向の厚さを２ｒａ、ｚｃ方向の高さを２ｚｂとする。
【００６５】
電流は＋ｘｃ軸を基点として、＋ｚｃから見て反時計回りの方向への角度をθとし、電流はこの方向に一様な電流密度ｊで流れているものとする。
【００６６】
粒子は円筒座標系で（Ｒｃ、θ／２、０）に位置するものとする。
【００６７】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体（円弧状柱状導体４５ｃ）の寸法であるｒａ、ｚｂ、θ、Ｒｃおよび粒子位置ベクトルを計算し記憶装置５に記憶する。
【００６８】
さらに、あらかじめ上記物理量が入力情報２１として記憶されている場合は、上記ステップは不要である。
【００６９】
以上が図５のステップ１０３の詳細である。
【００７０】
次に、磁場解析装置１の制御部３は、このステップまでに計算され記憶装置５に記憶されている磁性体を構成するＮ個の粒子の位置ベクトルと、あらかじめ記憶装置５に記憶されている多面体要素の形状により、磁性体を構成する粒子の磁場の運動方程式の係数を計算し記憶装置５に記憶する（図５のステップ１０４）。
【００７１】
ここで、図５のステップ１０４について、具体的に説明する。
【００７２】
分割した多面体要素（本実施形態では立方体要素）を２次元表示すると図９に示す形状となる。
【００７３】
ここで、粒子の位置ベクトルを（太字の）ｒｇ、粒子の位置ベクトルを重心とする多面体要素の要素境界面の中点をｑ点とし、粒子の位置ベクトルとｑ点との中間点にｐ点を定義する。
【００７４】
制御部３は、ステップ１０２において計算し記憶装置５に記憶している磁性体を構成する粒子（粒子を重心とする多面体要素（本実施形態では立方体要素））の位置ベクトルを読み込む。
【００７５】
なお、ステップ１１３において粒子の位置ベクトルが更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【００７６】
制御部３は、粒子の位置ベクトルとあらかじめ記憶装置５に記憶されている多面体要素の形状から、磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内のｐ点とｑ点の位置ベクトルと要素境界面への法線ベクトル（太字の）ｎを計算し、記憶装置５に記憶する。
【００７７】
次に、制御部３は、磁性体を構成する粒子を重心とする多面体要素の境界面積ΔＳ、および粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積ΔＶを計算し、記憶装置５に記憶する。
【００７８】
以上が図５のステップ１０４の詳細である。
【００７９】
次に、磁場解析装置１の制御部３は、ステップ１０３で記憶装置５に記憶されているコイルの寸法およびこのステップまでに記憶装置５に記憶されている電流密度ベクトルを用いて、ビオ・サバールの法則を積分することにより得られる解析解により、通電されたコイル４５が、磁性体を構成する粒子上に作る磁場ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０５）。
【００８０】
ここで、図５のステップ１０５について、具体的に説明する。
【００８１】
ここでは、任意の位置ベクトルが図４における磁性体としてのステータティース４３のある点Ｐであると仮定した場合に、通電されたコイル４５がＰ点の位置ベクトル上に作る磁場ベクトルを計算する手順を例に説明する。
【００８２】
まず、ローカル導体が直方体導体４５ａの場合について説明する。
【００８３】
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。
【００８４】
ローカル導体が直方体導体４５ａの場合は、図７に示すようにローカル座標系（ｘｓ,ｙｓ,ｚｓ）を適用する。
【００８５】
このローカル座標系はステップ１０３で説明している。
【００８６】
次に、Ｐ点の位置ベクトルをローカル座標系（ｘｓ,ｙｓ,ｚｓ）に変換する。
【００８７】
通電された直方体導体がＰ点に作る磁場ベクトルは、以下に示す式（１）〜（３）で記載される。
【００８８】
【数１】
2010160782.tif 000003
【００８９】
ここで、（太字の）ｒｐｓはローカル座標系（ｘｓ,ｙｓ,ｚｓ）でのＰ点の位置ベクトルであり、ｘｐｓ,ｙｐｓ,ｚｐｓはｘｓ,ｙｓ,ｚｓ方向の値である。 
πは円周率である。
【００９０】
Ｈｘｓ、Ｈｙｓ、Ｈｚｓはローカル座標系（ｘｓ,ｙｓ,ｚｓ）における磁場ベクトルの各成分である。 
ｊは電流密度である。
【００９１】
また、ｘｉ,ｙｊ,ｚｋはｘｓ,ｙｓ,ｚｓ方向の積分の上限、下限を表しており、式（４）に示す関係が成立する。
【００９２】
【数２】
2010160782.tif 000004
【００９３】
ここで、ａ，ｂ，ｃは直方体導体の寸法である。
【００９４】
次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
【００９５】
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。
【００９６】
ローカル導体が円弧状柱状導体４５ｃの場合は、図８に示すようにローカル座標系（ｘｃ,ｙｃ,ｚｃ）を適用する。
【００９７】
このローカル座標系はステップ１０３において説明している。
【００９８】
ローカル座標系（ｘｃ,ｙｃ,ｚｃ）に変換後のＰ点の位置ベクトルを以下の式（５）に示すように円筒座標系（太字の）ｒｐｃ＝（Ｒｐｃ、φｐｃ、Ｚｐｃ）に変換する。
【００９９】
通電された円弧状柱状導体４５ｃがＰ点に作る磁場ベクトルは、以下の式（６）〜（１１）で表される。
【０１００】
【数３】
2010160782.tif 000005
【０１０１】
ここで、Ｈｒｃ、Ｈｔｃ、Ｈｚｃは円筒座標系での磁場ベクトルの各成分である。 
ｊは電流密度である。
【０１０２】
ｒａ、θ、ｚｂは円弧状柱状導体４５ｃの寸法であり、Ｒｃは円弧の内径と外径の平均値である。
【０１０３】
ｓｇｎはＺｋの符号であり、Ｒｊ、Ｚｋは積分の上限、下限を表しており、式（１２）に示す関係が成立する。
【０１０４】
【数４】
2010160782.tif 000006
【０１０５】
以上が任意の点Ｐにコイルが作る磁場ベクトルの計算手順である。
【０１０６】
磁場解析装置１の制御部３は、ステップ１０４で既に計算され記憶装置５に記憶されている磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内部のｐ点の位置ベクトルを読み込む。
【０１０７】
次に、制御部３は、ｐ点の位置ベクトルをローカル座標系（ｘｓ,ｙｓ,ｚｓ）に変換し、ステップ１０３で計算され記憶装置５に記憶されている直方体導体４５ａの寸法と、ステップ１０１より入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込み、式（１）〜（４）に基づいて磁場ベクトルを計算する。
【０１０８】
なお、電流密度ベクトルがステップ１１６で更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【０１０９】
次に、制御部３は、式（１）〜（４）で計算された磁場ベクトルをグローバル座標系（ｘ,ｙ,ｚ）に変換し、変換後の磁場ベクトルを記憶装置５に記憶する。
【０１１０】
これにより、通電された直方体導体４５ａが、ｐ点に作る磁場ベクトルが求められる。
【０１１１】
さらに制御部３は、ｐ点の位置ベクトルをローカル座標系（ｘｃ,ｙｃ,ｚｃ）に変換し、さらに円筒座標系（太字の）ｒｐｃ＝（Ｒｐｃ、φｐｃ、Ｚｐｃ）に変換する。
【０１１２】
次に制御部３は、ステップ１０３で計算し記憶装置５に記憶している円弧状柱状導体４５ｃの寸法ｒａ、θ、ｚｂおよびＲｃを読み込み、また、制御部３はステップ１０１で入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込む。
【０１１３】
なお、電流密度ベクトルがステップ１１６で更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【０１１４】
制御部３は式（６）〜（１２）に基づいて磁場ベクトルを計算する。
【０１１５】
次に、計算された磁場ベクトルを直交座標系に変換し、さらにグローバル座標系（ｘ,ｙ,ｚ）に変換し記憶装置５に記憶する。
【０１１６】
これにより、通電された円弧状柱状導体４５ｃがｐ点に作る磁場ベクトルが求められる。
【０１１７】
以上が、図５のステップ１０５の詳細である。
【０１１８】
次に、磁場解析装置１の制御部３は、演算プログラム１９を用いて磁性体を構成する粒子の磁場の運動方程式から束縛を考慮せずに仮想時間刻みδｔ後の磁場ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０６）。
【０１１９】
ここで、図５のステップ１０６について、具体的に説明する。
【０１２０】
磁性体を構成するＮ個の粒子のラグランジアンを式（１３）〜式（１５）で表される形とする。
【０１２１】
【数５】
2010160782.tif 000007
【０１２２】
ここで、式（１３）において、αは仮想質量、太字のｒは位置ベクトル、太字のＨは磁場ベクトル、太字の傍点付きＨは磁場ベクトルの時間微分、太字のＭは磁化ベクトル、太字のｎは粒子の位置ベクトルを重心とする多面体要素境界面の法線ベクトル、太字のＨｅｘｔは外部からの印加磁場ベクトル、ΔＳは粒子の位置ベクトルを重心とする多面体要素の要素境界面の面積、ΔＶは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、χは磁気感受率、μ０は真空の透磁率、λはラグランジュの未定定数、πは円周率、ｓは粒子の位置ベクトルを重心とする多面体要素の要素境界面数である。
【０１２３】
また、各物理量の添え字ｉｐはｉ番目の粒子の位置ベクトルを重心とする多面体要素内でのｐ点の物理量、添え字ｊｑはｊ番目の粒子の位置ベクトルを重心とする多面体要素内でのｑ点の物理量を示す。
【０１２４】
ｐ点、ｑ点はステップ１０４において説明している。
【０１２５】
式（１３）において、ｍは粒子の質量、ｖは速度、φ（ｒｉ−ｒｊ）はｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギである。添え字ｉ、ｊはそれぞれ、ｉ、ｊ番目の粒子の物理量を示す。
【０１２６】
次に、正準変数を（太字の）Ｈｉｐ、（太字の傍点付き）Ｈｉｐとし、式（１３）〜式（１５）で示されるラグランジアンをラグランジュの運動方程式に代入すると、磁場の運動方程式は式（１６）のように記載できる。
【０１２７】
【数６】
2010160782.tif 000008
【０１２８】
ここで、式（１６）の右辺第２項は、ラグランジュの未定定数を通して束縛（磁化ベクトルの発散は０）を課している。
【０１２９】
式（１６）の右辺第３項は外部からの印加磁場ベクトルが変化したときにすばやく追従させるための減衰項であり、γは減衰定数である。
【０１３０】
式（１６）の右辺第２項に示される束縛を含んだ運動方程式を解くにあたり、本実施形態では、一般化された束縛の導入法であるＳＨＡＫＥ法を採用する。
【０１３１】
束縛を考慮せずに蛙跳び法により式（１６）を離散化すると以下の式（１７）、式（１８）、（１９）になる。
【０１３２】
【数７】
2010160782.tif 000009
【０１３３】
ここで、δｔは磁化現象の収束計算を行う上で用いる仮想時間刻みである。
【０１３４】
添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ−１／２は（ｎ−１／２）δｔにおける物理量、ｎ＋１／２は（ｎ＋１／２）δｔにおける物理量、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１３５】
磁場解析装置１の制御部３は、ステップ１０４において既に計算され記憶装置５に記憶されている粒子の位置ベクトルを重心とする多面体要素内のｐ点、ｑ点の位置ベクトルを読み込む。
【０１３６】
次に、制御部３は、ステップ１０５で既に計算され記憶装置５に記憶されているコイルが磁性体を構成する粒子を重心とする多面体要素内部のｐ点に作る磁場ベクトルを外部からの印加磁場ベクトルとして読み込む。
【０１３７】
さらに、制御部３は、ステップ１０４で既に計算され、記憶装置５に記憶されている、磁性体を構成する粒子の要素境界面積、法線ベクトルおよび粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
【０１３８】
また、制御部３はあらかじめ記憶装置５に記憶されている減衰定数、仮想質量、仮想時間刻みを読み込む。
【０１３９】
次に、制御部３はあらかじめ記憶装置５に記憶されている磁場ベクトル、磁場ベクトルの時間微分の初期値を読み込む。
【０１４０】
なお、制御部３は後述するステップ１０７において磁場ベクトル、磁場ベクトルの時間微分が更新されている場合は、その値を読み込む。
【０１４１】
なお、磁化ベクトルは、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に磁場ベクトルを代入することで計算される。
【０１４２】
制御部３は、式（１７）、式（１８）、（１９）を、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内の全てのｐ点に対して計算し、計算された磁場ベクトルを記憶装置５に記憶する。
【０１４３】
以上が、図５のステップ１０６の詳細である。
【０１４４】
次に、磁場解析装置１の制御部３はステップ１０６で計算し記憶装置５に記憶されている、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のｐ点の磁場ベクトルに束縛力を加え、計算された磁場ベクトルを、記憶装置５に記憶する（図５のステップ１０７）。
【０１４５】
ここで、図５のステップ１０７について、具体的に説明する。
【０１４６】
磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のｐ点の磁場ベクトルに以下に示す式（２０）、式（２１）に従って束縛力を加える。
【０１４７】
【数８】
2010160782.tif 000010
【０１４８】
ここで、太字のＨは磁場ベクトル、αは仮想質量、δｔは仮想時間刻み、γは減衰定数、Ｎは磁性体を構成する粒子数、ｓは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量、ｉｑは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｑ点の物理量、添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１４９】
磁場解析装置１の制御部３は、演算プログラム１９を用いて、ステップ１０６で計算され記憶装置５に記憶されている磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルを読み込む。
【０１５０】
制御部３は、あらかじめ記憶装置５に記憶されている仮想質量、仮想時間刻み、減衰定数を読み込む。
【０１５１】
制御部３は、磁性体を構成する粒子各々に対して式（２０）、式（２１）に基づく計算を行い、計算された磁場ベクトルを記憶装置５に記憶する。
【０１５２】
以上が、ステップ図５の１０７の詳細である。
【０１５３】
次に、磁場解析装置１の制御部３はステップ１０７で求めた磁場ベクトルが束縛条件を満たしているかを判断し、満たしていれば次のステップに進み、満たしていなければステップ１０７に戻る（図５のステップ１０８）。
【０１５４】
具体的には制御部３は、ステップ１０７で計算し記憶装置５に記憶されている磁性体を構成する粒子各々の磁場ベクトルから計算される磁化ベクトルを用いて式（２２）に基づく計算を行う。
【０１５５】
【数９】
2010160782.tif 000011
【０１５６】
ここで、添え字のｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１５７】
ｅｒｒｉは磁性体を構成する粒子の内、ｉ番目の粒子の束縛条件に対する誤差値である。
【０１５８】
磁化ベクトル（太字の）Ｍはステップ１０７において計算され記憶装置５に記憶されている磁場ベクトルを制御部３が読み込み、ステップ１０１より入力情報２１に記憶されている磁化曲線を表す関数に代入することで計算される。
【０１５９】
法線ベクトル（太字の）ｎには、ステップ１０４において既に計算され記憶装置５に記憶されているものを制御部３は読み込み代入する。
【０１６０】
制御部３は、すべての粒子に対して、誤差の値が式（２３）を満たさなければステップ１０７に戻る。
【０１６１】
【数１０】
2010160782.tif 000012
【０１６２】
式（２３）においてＡ１は任意の誤差判別値であり、あらかじめ記憶装置５に任意の値が記憶されている。
【０１６３】
以上が図５のステップ１０８の詳細である。
【０１６４】
次に、磁場解析装置１の制御部３は、磁性体の磁化現象が定常状態に到達したかを判断し、条件を満たしていれば次のステップに進む。（図５のステップ１０９）
ここで、図５のステップ１０９について、具体的に説明する。
【０１６５】
磁性体を構成する粒子が定常状態に到達したかは、以下の式（２４）により判断される。
【０１６６】
【数１１】
2010160782.tif 000013
【０１６７】
ここで、太字のＨは磁場ベクトル、μ0は真空の透磁率、Ｎは磁性体を構成する粒子の粒子数、ｓは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点での物理量である。
【０１６８】
さらに、添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ―１は（ｎ―１）δｔ、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１６９】
また、Ａ２は磁性体の磁場ベクトルが定常状態に到達したかを判断するための任意の誤差判定値である。
【０１７０】
磁場解析装置１の制御部３は、あらかじめ記憶装置５に記憶されているＡ２、μ０を読み込む。
【０１７１】
また、制御部３は、ステップ１０７で計算され記憶装置５に記憶されている、磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルを読み込む。
【０１７２】
次に、制御部３は式（２４）を計算し、式（２４）を満たしていれば次のステップに進み、満たしていなければステップ１０６に戻る。
【０１７３】
以上が図５のステップ１０９の詳細である。
【０１７４】
次に、磁場解析装置１の制御部３は、磁性体を構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁化ベクトル、磁束密度ベクトルと、コイルを構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁束密度ベクトルを計算し、記憶装置５に記憶する（図５のステップ１１０）。
【０１７５】
以下にステップ１１０を具体的に説明する。
【０１７６】
（磁性体を構成する粒子）
磁性体を構成する粒子の位置ベクトル上の磁場ベクトルは以下の式（２５）で表される。
【０１７７】
【数１２】
2010160782.tif 000014
【０１７８】
ここで、（太字の）ｅｘ＝（１,０,０）、（太字の）ｅｙ＝（０,１,０）、（太字の）ｅｚ＝（０,０,１）である。
【０１７９】
太字のＨは磁場ベクトル、太字のｎは法線ベクトルであり、添え字ｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を表し、添え字ｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量である。
【０１８０】
磁性体を構成する粒子の磁化ベクトルは、磁場ベクトルを磁化曲線を表す関数に代入することで求まる。
【０１８１】
磁性体を構成する粒子の位置ベクトル上での磁束密度ベクトルは、式（２６）で表される。
【０１８２】
【数１３】
2010160782.tif 000015
【０１８３】
ここで、太字のＢは磁束密度ベクトル、太字のＨは磁場ベクトル、太字のＭは磁化ベクトルであり、添え字のｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を指し示す。
【０１８４】
μ０は真空の透磁率である。
【０１８５】
磁場解析装置１の制御部３は、ステップ１１０までで計算されている、磁性体を構成する粒子を重心とする多面体要素内のｐ点での磁場ベクトルを記憶装置５から読み込む。
【０１８６】
次に、制御部３は、ステップ１０４で既に計算され記憶装置５に記憶されている法線ベクトルを読み込み、式（２５）に基づいて磁性体を構成する粒子の位置ベクトル上での磁場ベクトルを計算し記憶装置５に記憶する。
【０１８７】
さらに、制御部３は、このステップで計算された磁場ベクトルを、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に代入し、磁性体を構成する粒子の磁化ベクトルを計算し記憶装置５に記憶する。
【０１８８】
次に制御部３は、このステップで計算された磁場ベクトルと磁化ベクトルを式（２６）に代入し、磁性体を構成する粒子の位置ベクトル上の磁束密度ベクトルを計算し記憶装置５に記憶する。
【０１８９】
なお、真空の透磁率はあらかじめ記憶装置５に記憶されているものを制御部３が読み込む。
【０１９０】
（コイルを構成する粒子）
コイルを構成する粒子の位置ベクトル上の磁場ベクトルは式（２７）、式（２８）で記載される。
【０１９１】
【数１４】
2010160782.tif 000016
【０１９２】
ここで、太字のｒは位置ベクトル、太字のＨは磁場ベクトル、Ｎは磁性体を構成する粒子の数、ｓおよびΔＳは磁性体を構成する粒子を重心とする多面体要素の境界面数と境界面積、太字のＭは磁化ベクトル、太字のｎは法線ベクトルであり、添え字のｉはコイルを構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｊは磁性体を構成する粒子の内、ｊ番目の粒子の位置ベクトル上での物理量であり、ｊｑは磁性体を構成する粒子の内、ｊ番目の粒子の位置ベクトルを重心とする多面体要素のｑ点での物理量を指し示す。
【０１９３】
式（２７）の右辺第一項は、磁性体を構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルを記述する項であり、式（２８）で表される。
【０１９４】
式（２７）の右辺第二項は、コイルを構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルであり式（１）〜（１２）で表される。
【０１９５】
磁場解析装置１の制御部３は、ステップ１０３において計算し記憶装置５に記憶されているコイルを構成する粒子の位置ベクトルを記憶装置５より読み込む。
【０１９６】
なお、ステップ１１３においてコイルを構成する粒子の位置ベクトルが更新されている場合は、制御部３はその値を読み込む。
【０１９７】
次に、制御部３はステップ１０２より計算され記憶装置５に記憶されている磁性体を構成する粒子の数と多面体要素の面数、および粒子の位置ベクトルを重心とする多面体要素内のｑ点の位置ベクトルと多面体要素境界面積を読み込む。
【０１９８】
また、制御部３はこのステップで計算され記憶装置５に記憶されている、磁性体を構成する粒子の磁化ベクトルを読み込む。
【０１９９】
さらに、制御部３はステップ１０３において計算され記憶装置５に記憶されているコイルの寸法およびステップ１０１で入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込む。
【０２００】
なお、ステップ１１６において電流密度ベクトルが更新されている場合は、制御部３はその値を読み込む。
【０２０１】
その後、制御部３は式（２７）、式（２８）および式（１）〜（１２）に従い、コイルを構成する粒子の位置ベクトル上での磁場ベクトルを計算し、記憶装置５に記憶する。
【０２０２】
制御部３は、コイルを構成する粒子の位置ベクトル上での磁束密度ベクトル（太字の）Ｂｉを、すでに計算されているコイルを構成する粒子の位置ベクトル上の磁場ベクトルに真空の透磁率を掛けて計算し記憶装置５に記憶する。
【０２０３】
真空の透磁率は、あらかじめ記憶装置５に記憶されている。
【０２０４】
以上が、ステップ１１０の詳細である。
【０２０５】
このように、制御部３は、磁場の運動方程式を解くことにより磁性体の磁化現象を解析する。
【０２０６】
そのため、有限要素法のように空間全域をメッシュ分割する必要がなく、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる。
【０２０７】
また、運動方程式を解くため、行列を扱わず、計算に必要なメモリ量は粒子数に比例する。
【０２０８】
次に、磁場解析装置１の制御部３はＳＰＭモータ３１を構成する全ての粒子に働く力ベクトルを計算し、記憶装置５に記憶する（図５のステップ１１１）。
【０２０９】
以下に図５のステップ１１１を具体的に説明する。
【０２１０】
（磁性体を構成する粒子に働く力ベクトル）
正準変数を太字のｒi、太字の傍点付きｒiとし、式（１３）のラグランジアンをラグランジュの運動方程式に代入する。
【０２１１】
すると、磁性体を構成する粒子に働く力は式（２９）のように記載される。
【０２１２】
【数１５】
2010160782.tif 000017
【０２１３】
ここで、μ０は真空の透磁率、太字のＨは磁場ベクトル、太字のＭは磁化ベクトル、太字のｒは位置ベクトル、ΔＶは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、Ｎは磁性体を構成する粒子数であり、添え字のｉ、ｊは磁性体を構成する粒子の内、ｉ番目、ｊ番目の粒子の物理量を指し示し、添え字のｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量を指し示す。
【０２１４】
また、φ（ｒｉ−ｒｊ）は、ｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギである。
【０２１５】
本実施形態ではＳＰＭモータを剛体（粒子間の相対距離が不変）とする。
【０２１６】
剛体モデルとすると、磁性体を構成する粒子に働く力は式（３０）のように記載される。
【０２１７】
【数１６】
2010160782.tif 000018
【０２１８】
磁場解析装置１の制御部３は、ステップ１１０までにおいて計算され記憶装置５に記憶されている磁性体を構成する粒子を重心とする多面体要素内のｐ点での磁場ベクトルを読み込む。
【０２１９】
制御部３は、ｐ点での磁場ベクトルを、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に代入することによりｐ点での磁化ベクトルを計算し記憶装置５に記憶する。
【０２２０】
また、制御部３は、ステップ１０４で計算され記憶装置５に記憶されている磁性体を構成する粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
【０２２１】
次に、制御部３は、磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルの偏微分を計算し、その値に基づいて、式（３０）より力ベクトルを計算し記憶装置５に記憶する。
【０２２２】
（コイルを構成する粒子に働く力ベクトル）
コイルを構成する粒子のうち、ｉ番目の粒子の位置ベクトルを太字のｒｉ、その位置ベクトルでの電流の単位ベクトルを太字のｔｉ、磁束密度ベクトルを太字のＢｉ、電流の流れる方向のコイルの長さをＬｉとすると、コイルを構成する粒子の内、ｉ番目の粒子に働く力ベクトルは以下の式（３１）で記載される。
【０２２３】
【数１７】
2010160782.tif 000019
【０２２４】
磁場解析装置１の制御部３は、ステップ１０３において計算され記憶装置５に記憶されているコイルを構成する粒子の位置ベクトルを読み込む。
【０２２５】
なお、ステップ１１３においてコイルを構成する粒子の位置ベクトルが更新されている場合は制御部３はその値を読み込む。
【０２２６】
また、制御部３はステップ１０３において計算され記憶装置５に記憶されているコイルの寸法を読み込む。
【０２２７】
次に、制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されているコイルの電流密度ベクトルを読み込む。
【０２２８】
なお、制御部３は、ステップ１１６においてコイルの電流密度ベクトルが更新されている場合はその値を用いる。
【０２２９】
制御部３は、式（３１）に基づいて、コイルを構成する粒子に働く力ベクトルを計算し記憶装置５に記憶する。
【０２３０】
以上が、図５のステップ１１１の詳細である。
【０２３１】
次に、磁場解析装置１の制御部３は、ステップ１１１で計算した可動部（移動子）を構成する粒子に働く力ベクトルから、可動部（移動子）の回転量、角速度、並進量、並進速度を計算し、記憶装置５に記憶する（図５のステップ１１２）。
【０２３２】
以下に図５のステップ１１２を具体的に説明する。
【０２３３】
剛体の並進運動の運動方程式は式（３２）で表される。
【０２３４】
【数１８】
2010160782.tif 000020
【０２３５】
ここで、ｍｔｏｔは可動部（移動子）の質量、太字のｒｃは可動部（移動子）において任意に定めた代表点の位置ベクトル、太字のＦは可動部（移動子）を構成する粒子に働く力ベクトル、Ｎは可動部（移動子）を構成する粒子数である。
【０２３６】
次に、任意の固定点の周りの回転運動の運動方程式は式（３３−１）、式（３３−２）、式（３３−３）で表される。
【０２３７】
【数１９】
2010160782.tif 000021
【０２３８】
ここで、Ｉ１、Ｉ２、Ｉ３は剛体の主慣性モーメント、ω１、ω２、ω３は慣性主軸座標系での角速度、Ｎ１、Ｎ２、Ｎ３は慣性主軸座標系でのトルクである。
【０２３９】
本実施例のＳＰＭモータ３１のロータ３１はロータシャフト４１を通じてある回転軸に固定されている。
【０２４０】
固定された回転軸をｚ軸とすると、ロータ３１の運動は以下のように計算される。
【０２４１】
ロータ３１はロータシャフト４１を通じてｚ軸に固定されているため、ｚ軸周りの回転運動の運動方程式を解くことで、ロータ３１の運動は計算される。
【０２４２】
ｚ軸周りのロータ３１の回転運動の運動方程式は式（３４）で表される。
【０２４３】
【数２０】
2010160782.tif 000022
【０２４４】
ここで、Ｉｚはロータ３１の慣性モーメント、ωｚはｚ方向の角速度、Ｎｚはｚ方向のトルクである。
【０２４５】
式（３４）を蛙飛び法で離散化すると、以下の式（３５）、式（３６）を得る。
【０２４６】
【数２１】
2010160782.tif 000023
【０２４７】
ここで、θｚはｚ軸周りの回転量であり、ｄｔは微小時間刻み幅である。
【０２４８】
また、添え字のｎ−１／２、ｎ、ｎ＋１／２、ｎ＋１はそれぞれ（ｎ―１／２）ｄｔ、ｎｄｔ、（ｎ＋１／２）ｄｔ、（ｎ＋１）ｄｔ秒での物理量に対応している。
【０２４９】
時間刻み幅は記憶装置５にあらかじめ任意の値が記憶されているものを制御部３が読み込む。
【０２５０】
回転量、角速度の初期値は記憶装置５にあらかじめ任意の値が記憶装置５に記憶されている。
【０２５１】
なお、すでにステップ１１２においてロータ３１の回転量、角速度が更新されている場合はその値を用いる。
【０２５２】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されているロータ３１の三次元構造と質量密度により慣性モーメントを計算する。
【０２５３】
さらに、制御部３は、ロータ３１を構成する粒子の位置ベクトルと、ステップ１１１で計算し記憶装置５に記憶されているロータ３１を構成する粒子に働く力ベクトルに基づいて、式（３５）、式（３６）を計算し、ｄｔ秒後のロータ３１の回転量と角速度を計算し記憶装置５に記憶する。
【０２５４】
以上が、図５のステップ１１２の詳細である。
【０２５５】
次に、磁場解析装置１の制御部３はステップ１１２で計算した微小時間後の回転量や重心の移動量を基に、移動子を構成している粒子の位置ベクトルを計算し、記憶する（図５のステップ１１３）。
【０２５６】
次に、磁場解析装置１の制御部３は必要に応じてＳＰＭモータ３１を構成する粒子の位置ベクトル、力ベクトル、磁場ベクトル、磁束密度ベクトル、磁化ベクトルをプリンタポート１１を介してプリンタ１２より出力する（図５のステップ１１４）。
【０２５７】
次に、磁場解析装置１の制御部３は所定の終了条件（時間、移動量等）を満たしているかを判断し、満たしている場合は解析を終了し、満たしていない場合はステップ１１６に進む（図５のステップ１１５）。
【０２５８】
次に、磁場解析装置１の制御部３は必要に応じてコイル４５に流れる電流密度ベクトルを更新し、ステップ１０４に戻る（図５のステップ１１６）。
【０２５９】
以上が本実施形態に係る磁場解析の手順である。
【０２６０】
このように、本実施形態によれば、磁場解析装置１が粒子のラグラジアンから導出される運動方程式を解くことで磁性体の磁化現象を計算する。
【０２６１】
そのため、有限要素法のように空間全域をメッシュ分割する必要がなく、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる。
【０２６２】
また、行列を扱わないので、計算に必要なメモリは粒子（要素）数に比例する。非対称密行列を扱う計算に比べてメモリを大幅に減少させることができる。
【０２６３】
また、本実施形態によれば、磁性体を構成する粒子のラグランジアンから導出された磁場の運動方程式の解を演算すると、外部からの入力が変化しない場合、系の全エネルギーが保存される。
【産業上の利用可能性】
【０２６４】
上記した実施形態では、本発明を磁性体を構成する粒子のラグランジアンから導出される運動方程式を解くことで磁性体の磁化現象を計算した場合について説明したが、本発明は、何等、これに限定されることなく、運動方程式であればラグランジアンから導出されたもの以外のものであってもよい。
【０２６５】
また、上記した実施形態では本発明をＳＰＭモータ３１の磁場解析に使用したが、本発明は、何等、これに限定されることなく、磁性体が配置された空間における任意の磁場を解析するものであれば、例えばリニアモータ等のＳＰＭモータ３１以外のモータ、あるいは磁性体で空間を囲む磁気シールドにおいて、磁性体でシールドされた空間内の磁場解析に用いても良い。
【符号の説明】
【０２６６】
１…………磁場解析装置
３…………制御部
５…………記憶装置
７…………入力部
９…………表示部
１１………プリンタポート
１２………プリンタ
１３………バス
３１………ＳＰＭモータ
３３………ロータ
３５………ステータ
３７………ロータコア
３９………永久磁石
４１………ロータシャフト
４３………ステータティース
４４………コアバック
４５………コイル
４６………フレーム
【図１】
2010160782.tif 000030
【図２】
2010160782.tif 000031
【図３】
2010160782.tif 000032
【図４】
2010160782.tif 000033
【図５】
2010160782.tif 000034
【図６】
2010160782.tif 000035
【図７】
2010160782.tif 000036
【図８】
2010160782.tif 000037
【図９】
2010160782.tif 000038