磁場解析装置、磁場解析方法、及び連成解析装置
【課題】高速かつ高精度に磁場解析を行う。
【解決手段】磁場解析装置は、物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより物体に作用する磁場を演算する磁場演算部を備える。磁場の運動方程式には調整可能なパラメタが導入されており、磁場演算部は、物体の磁化曲線の非線形性が数値解の収束性に与える影響を緩和するようパラメタを調整する。
【解決手段】磁場解析装置は、物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより物体に作用する磁場を演算する磁場演算部を備える。磁場の運動方程式には調整可能なパラメタが導入されており、磁場演算部は、物体の磁化曲線の非線形性が数値解の収束性に与える影響を緩和するようパラメタを調整する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、磁場解析装置および磁場解析方法に関する。
【背景技術】
【０００２】
例えば特許文献１には、磁場解析を実行する演算処理装置を具えるモータ解析装置が記載されている。演算処理装置は、ユーザの操作に基づく外部指令に応じて、有限要素法による静磁場解析やマクスウェルの応力法によるトルク計算を実行する。有限要素法による静磁場解析のためにメッシュ分割が行われる。メッシュ分割は、コア領域及びハウジング領域、さらには外部空気層領域についても行われる。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平１１−１４６６８８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかし、静磁場解析は時間的な変動を考慮しないから、駆動状態のモータの電磁気的挙動を必ずしも十分に反映することができない。その結果、解析の精度も不十分となりがちである。また、演算の所要時間は短いことが実用上好ましい。
【０００５】
そこで、本発明は、高速かつ高精度に解析結果を得ることのできる磁場解析装置及び磁場解析方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明の一態様によれば、物体を含む系における磁場を解析するための磁場解析装置が提供される。この装置は、物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより前記物体に作用する磁場を演算する磁場演算部を備える。前記磁場の運動方程式には調整可能なパラメタが導入されており、前記磁場演算部は、前記物体の磁化曲線の非線形性が数値解の収束性に与える影響を緩和するよう前記パラメタを調整する。
【０００７】
この態様によれば、運動方程式を数値的に解くので、有限要素法に基づく周知の磁場解析法とは異なり逆行列の演算が不要である。このため、少ない計算資源で精度よく磁場解析を実行することができる。また、磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整可能パラメタが調整されるので、磁気飽和に伴う数値解の収束遅れを低減することができる。
【０００８】
本発明の他の態様は、磁場解析方法である。この方法は、磁化曲線が非線形である材料を含む物体に作用する磁場が満たすべき関係を運動方程式の形式で記述した磁場の運動方程式を数値的に解く磁場演算ステップと、前記磁場の運動方程式に導入されている調整可能パラメタを、前記磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整する調整ステップと、備える。
【０００９】
本発明の一実施例は、この磁場解析方法をコンピュータに実行させるためのプログラムであってもよいし、このプログラムを記録したコンピュータ読み取り可能な記録媒体であってもよい。
【００１０】
本発明の他の態様は、連成解析装置である。この装置は、物体を構成する各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式に基づき磁場を演算する磁場解析演算部と、演算された磁場に基づき前記物体の運動を演算する機構解析演算部と、を備える。
【００１１】
本発明の他の態様は、第１の物理量と、該第１の物理量と非線形の関係をもつ第２の物理量とを有する物体の特性を解析するための解析装置である。この装置は、運動方程式の形式で記述された第１の物理量が満たすべき関係式を数値的に解く数値演算部を備え、前記関係式は、第２の物理量に依存する項と、第２の物理量に依存しない項と、を含み、第２の物理量に依存しない項の係数が第２の物理量に依存する項に連動するよう定義されている。
【発明の効果】
【００１２】
本発明によれば、高速かつ高精度に磁場解析を行うことができる。
【図面の簡単な説明】
【００１３】
【図１】本発明の一実施形態に係る磁場解析装置のハードウェア構成を示す図である。
【図２】本発明の一実施形態に係る磁場解析装置の記憶装置を示す図である。
【図３】本発明の一実施形態に係るＳＰＭモータの構成の概略を示す図である。
【図４】図３のＳＰＭモータのステータティースの１つの周囲の拡大斜視図である。
【図５】本発明の一実施形態に係る磁場解析装置を用いた磁場の解析の手順を示すフローチャートである。
【図６】図６（ａ）はコイルの拡大斜視図であって、図６（ｂ）はコイルをローカル導体に分割した例を示す図である。
【図７】本発明の一実施形態に係る直方体導体とローカル座標の関係を示す図である。
【図８】本発明の一実施形態に係る円弧状柱状導体とローカル座標の関係を示す図である。
【図９】本発明の一実施形態に係る磁性体の磁場計算に用いる要素の説明図である。
【図１０】磁性体（例えば鉄）の磁化曲線の一例を模式的に示す図である。
【図１１】本発明の一実施形態に係る磁場演算処理を説明するためのフローチャートである。
【図１２】本発明の一実施形態による磁場演算の例を示す図である。
【図１３】本発明の一実施形態による磁場演算の例を示す図である。
【発明を実施するための形態】
【００１４】
以下、図面に基づいて本発明に好適な実施形態を詳細に説明する。
【００１５】
まず、図１を参照して、本実施形態に係る磁場解析装置１のハードウェア構成を説明する。
【００１６】
なお、図１のハードウェア構成は例示であり、これに限定されないのは当然である。
【００１７】
図１に示すように、磁場解析装置１は制御部３、記憶装置５、メディア入出力部６、入力部７、表示部９、プリンタポート１１等がバス１３を介して互いに接続されている。
【００１８】
制御部３は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等で構成され、記憶手段としての記憶装置５に格納されたプログラムに従って、バス１３を介して接続された各装置を駆動制御する。
【００１９】
図２に示すように、記憶装置５には、磁場解析装置１の各構成部分を駆動制御するための制御プログラム１５、本発明を実施するための磁場解析プログラム１７が格納されている。
【００２０】
磁場解析プログラム１７は、解析条件を有する情報である入力情報２１と、入力情報２１と磁場の運動方程式に基づいて磁場を演算し力学系の運動方程式を演算する演算プログラム１９とを有している。
【００２１】
メディア入出力部６は、フロッピー（登録商標）ディスク、ＣＤ、ＤＶＤ等のメディアとの間で情報の入出力を行う装置である。
【００２２】
入力部７は、キーボード、マウス等の入力装置であり、表示部９はディスプレイ等の表示機器である。
【００２３】
プリンタポート１１には出力装置としてのプリンタ１２等が接続される。
【００２４】
次に、磁場解析装置１を用いた磁場の解析の手順について図３〜９を参照して説明する。
【００２５】
ここでは永久磁石モータ（Permanent Magnet Motor）の一種である、ＳＰＭモータ３１（Surface Permanent Magnet Motor）の磁場解析を例にして説明する。
【００２６】
まず、ＳＰＭモータ３１の構成の概略を図３および図４を参照して説明する。
【００２７】
図３に示すように、ＳＰＭモータ３１は回転子（移動子）であるロータ３３と固定子であるステータ３５を有している。
【００２８】
ロータ３３は鉄等の磁性体である円柱状のロータコア３７を有し、ロータコア３７の表面には永久磁石３９が設けられている。
【００２９】
ロータコア３７の軸中心には棒状のロータシャフト４１が設けられている。
【００３０】
ステータ３５は磁性体である歯状のステータティース４３とステータティース４３の外側に設けられた円筒状の磁性体であるコアバック４４、コアバック４４の外側に設けられた円筒状のフレーム４６から構成されている。
【００３１】
図３および図４に示すように、ステータティース４３には、金属等の導電体であるコイル４５が巻きつけられている。
【００３２】
なお、実際のＳＰＭモータ３１ではコイル４５は仕様に応じたターン数でステータティース４３に巻きつけられて束となっているが、本実施形態では、図３および図４に描かれているように、コイル一本一本をモデル化せず、コイルの束を一つの導体として扱う。
【００３３】
このような構造のＳＰＭモータ３１は、永久磁石３９の磁場、およびコイル４５に電流を流すことにより発生する磁場によって、ロータ３３、ステータ３５が磁化する。磁性体の磁気エネルギーの偏差によりＳＰＭモータ３１は駆動する。
【００３４】
そのため、ＳＰＭモータ３１の磁場解析を行うためにはコイル４５、永久磁石３９がロータ３３、ステータ３５を構成する磁性体上に作る磁場ベクトルを計算し、これら磁性体の磁化現象を解析する必要がある。
【００３５】
次に、解析の手順について図５〜図９を参照して説明する。
【００３６】
なお、以下の手順においてはＳＰＭモータ３１を、複数に要素分割して要素ごとの粒子の集合体とし、剛体モデルとして扱っているが、本発明はこれに限定されることはなく、粒子を用いずに要素ごとの剛体モデルとして扱ってもよい。
【００３７】
また、以下の手順において、磁性体とは、ロータ３３、ステータ３５を構成する磁性体と永久磁石３９を指し示し、磁化曲線を表す関数によってこれらは区別される。
【００３８】
まず、磁場解析装置１の制御部３は磁場解析プログラム１７を起動し、解析したいＳＰＭモータ３１の解析条件としての三次元構造（形状、座標点）、質量密度、磁性体の磁化曲線を表す関数、導体の電流密度ベクトルを記憶装置５の入力情報２１として記憶する（図５のステップ１０１）。
【００３９】
これらの物理量は例えばメディア入出力部６を介してＣＤ−ＲＯＭ等の記録媒体から読込んだものであってもよい。
【００４０】
また、ＳＰＭモータ３１の三次元構造の情報とは例えばＣＡＤ等のデータである。
【００４１】
さらに、あらかじめ上記物理量が入力情報２１として記憶されている場合は、上記ステップは不要である。
【００４２】
電流密度ベクトルは、コイル４５の作る磁場ベクトルを計算する際に必要になる。
【００４３】
さらに、磁性体の磁化曲線を表す関数は磁化ベクトル（太字の）Ｍを計算する際に必要になる。
【００４４】
以上が図５のステップ１０１の詳細である。
【００４５】
次に、磁場解析装置１の制御部３は、入力情報２１の有する三次元構造の情報から、磁性体をＮ個の粒子（粒子の位置ベクトルを重心とする多面体要素（本実施例では立方体要素））に分割し、粒子の位置ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０２）。
【００４６】
ここで、Ｎは任意の整数であり、粒子の数Ｎおよび位置ベクトルは入力情報２１の有する三次元構造およびあらかじめ記憶装置５に記憶されている多面体要素の形状により計算される。
【００４７】
以上が図５のステップ１０２の詳細である。
【００４８】
次に、制御部３は、入力情報２１の有する三次元構造の情報から、導体（図６（ａ）に示すコイル４５）を図６（ｂ）に示すようにローカル導体（直方体導体４５ａ、４５ｂと円弧状柱状導体４５ｃ、４５ｄ）に分割し、それぞれの導体が作る磁場ベクトルを計算するための係数を計算し、記憶装置５に記憶する（図５のステップ１０３）。
【００４９】
ここで、図５のステップ１０３について、具体的に説明する。
【００５０】
まず、ローカル導体が直方体導体４５ａの場合について説明する。
【００５１】
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。
【００５２】
ローカル導体が直方体導体４５ａの場合は、制御部３は、図７に示すように、ローカル座標系（ｘｓ,ｙｓ,ｚｓ）を適用する。
【００５３】
このローカル座標系においてはローカル導体（直方体導体４５ａ）の重心を原点Ｏｓとし、直方体導体４５ａの寸法はｘｓ方向に２ａ、ｙｓ方向に２ｂ、ｚｓ方向に２ｃの長さを持つものとする。
【００５４】
また原点Ｏｓに粒子は位置するものとする。
【００５５】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体（直方体導体４５ａ）の寸法であるａ，ｂ，ｃと粒子位置ベクトルを計算し記憶装置５に記憶する。
【００５６】
次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
【００５７】
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。
【００５８】
ローカル導体が円弧状柱状導体４５ｃの場合は、制御部３は、図８に示すようにローカル座標系（ｘｃ,ｙｃ,ｚｃ）を適用する。
【００５９】
このローカル座標系においては原点Ｏｃは円弧の中心軸上に存在し、かつ円弧状柱状導体４５ｃの高さ方向（図８のzｃ方向）に対して円弧状柱状導体４５ｃが対称となる点に存在するものとする。
【００６０】
また、ｘｃ,ｙｃ,ｚｃは、ｘｃ−ｙｃ平面でみると、＋ｘｃ軸を基点とし、円弧状柱状導体４５ｃの円弧が＋ｚｃ軸からみて反時計回りになるようして決定する。
【００６１】
円弧状柱状導体４５ｃの内径と外径の平均値をＲｃとし、径方向の厚さを２ｒａ、ｚｃ方向の高さを２ｚｂとする。
【００６２】
電流は＋ｘｃ軸を基点として、＋ｚｃから見て反時計回りの方向への角度をθとし、電流はこの方向に一様な電流密度ｊで流れているものとする。
【００６３】
粒子は円筒座標系で（Ｒｃ、θ／２、０）に位置するものとする。
【００６４】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体（円弧状柱状導体４５ｃ）の寸法であるｒａ、ｚｂ、θ、Ｒｃおよび粒子位置ベクトルを計算し記憶装置５に記憶する。
【００６５】
さらに、あらかじめ上記物理量が入力情報２１として記憶されている場合は、上記ステップは不要である。
【００６６】
以上が図５のステップ１０３の詳細である。
【００６７】
次に、磁場解析装置１の制御部３は、このステップまでに計算され記憶装置５に記憶されている磁性体を構成するＮ個の粒子の位置ベクトルと、あらかじめ記憶装置５に記憶されている多面体要素の形状により、磁性体を構成する粒子の磁場の運動方程式の係数を計算し記憶装置５に記憶する（図５のステップ１０４）。
【００６８】
ここで、図５のステップ１０４について、具体的に説明する。
【００６９】
分割した多面体要素（本実施形態では立方体要素）を２次元表示すると図９に示す形状となる。
【００７０】
ここで、粒子の位置ベクトルを（太字の）ｒｇ、粒子の位置ベクトルを重心とする多面体要素の要素境界面の中点をｑ点とし、粒子の位置ベクトルとｑ点との中間点にｐ点を定義する。
【００７１】
制御部３は、ステップ１０２において計算し記憶装置５に記憶している磁性体を構成する粒子（粒子を重心とする多面体要素（本実施形態では立方体要素））の位置ベクトルを読み込む。
【００７２】
なお、ステップ１１３において粒子の位置ベクトルが更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【００７３】
制御部３は、粒子の位置ベクトルとあらかじめ記憶装置５に記憶されている多面体要素の形状から、磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内のｐ点とｑ点の位置ベクトルと要素境界面への法線ベクトル（太字の）ｎを計算し、記憶装置５に記憶する。
【００７４】
次に、制御部３は、磁性体を構成する粒子を重心とする多面体要素の境界面積ΔＳ、および粒子の位置ベクトルを重心とする多面体要素の境界面の頂点と粒子位置を頂点とする多面体要素体積ΔＶを計算し、記憶装置５に記憶する。
【００７５】
以上が図５のステップ１０４の詳細である。
【００７６】
次に、磁場解析装置１の制御部３は、ステップ１０３で記憶装置５に記憶されているコイルの寸法およびこのステップまでに記憶装置５に記憶されている電流密度ベクトルを用いて、ビオ・サバールの法則を積分することにより得られる解析解により、通電されたコイル４５が、磁性体を構成する粒子上に作る磁場ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０５）。
【００７７】
ここで、図５のステップ１０５について、具体的に説明する。
【００７８】
ここでは、任意の位置ベクトルが図４における磁性体としてのステータティース４３のある点Ｐであると仮定した場合に、通電されたコイル４５がＰ点の位置ベクトル上に作る磁場ベクトルを計算する手順を例に説明する。
【００７９】
まず、ローカル導体が直方体導体４５ａの場合について説明する。
【００８０】
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。
【００８１】
ローカル導体が直方体導体４５ａの場合は、図７に示すようにローカル座標系（ｘｓ,ｙｓ,ｚｓ）を適用する。
【００８２】
このローカル座標系はステップ１０３で説明している。
【００８３】
次に、Ｐ点の位置ベクトルをローカル座標系（ｘｓ,ｙｓ,ｚｓ）に変換する。
【００８４】
通電された直方体導体がＰ点に作る磁場ベクトルは、以下に示す式（１）〜（３）で記載される。
【００８５】
【数１】
【００８６】
ここで、（太字の）ｒｐｓはローカル座標系（ｘｓ,ｙｓ,ｚｓ）でのＰ点の位置ベクトルであり、ｘｐｓ,ｙｐｓ,ｚｐｓはｘｓ,ｙｓ,ｚｓ方向の値である。
πは円周率である。
【００８７】
Ｈｘｓ、Ｈｙｓ、Ｈｚｓはローカル座標系（ｘｓ,ｙｓ,ｚｓ）における磁場ベクトルの各成分である。
ｊは電流密度である。
【００８８】
また、ｘｉ,ｙｊ,ｚｋはｘｓ,ｙｓ,ｚｓ方向の積分の上限、下限を表しており、式（４）に示す関係が成立する。
【００８９】
【数２】
【００９０】
ここで、ａ，ｂ，ｃは直方体導体の寸法である。
【００９１】
次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
【００９２】
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。
【００９３】
ローカル導体が円弧状柱状導体４５ｃの場合は、図８に示すようにローカル座標系（ｘｃ,ｙｃ,ｚｃ）を適用する。
【００９４】
このローカル座標系はステップ１０３において説明している。
【００９５】
ローカル座標系（ｘｃ,ｙｃ,ｚｃ）に変換後のＰ点の位置ベクトルを以下の式（５）に示すように円筒座標系（太字の）ｒｐｃ＝（Ｒｐｃ、φｐｃ、Ｚｐｃ）に変換する。
【００９６】
通電された円弧状柱状導体４５ｃがＰ点に作る磁場ベクトルは、以下の式（６）〜（１１）で表される。
【００９７】
【数３】
【００９８】
ここで、Ｈｒｃ、Ｈｔｃ、Ｈｚｃは円筒座標系での磁場ベクトルの各成分である。
ｊは電流密度である。
【００９９】
ｒａ、θ、ｚｂは円弧状柱状導体４５ｃの寸法であり、Ｒｃは円弧の内径と外径の平均値である。
【０１００】
ｓｇｎはＺｋの符号であり、Ｒｊ、Ｚｋは積分の上限、下限を表しており、式（１２）に示す関係が成立する。
【０１０１】
【数４】
【０１０２】
以上が任意の点Ｐにコイルが作る磁場ベクトルの計算手順である。
【０１０３】
磁場解析装置１の制御部３は、ステップ１０４で既に計算され記憶装置５に記憶されている磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内部のｐ点の位置ベクトルを読み込む。
【０１０４】
次に、制御部３は、ｐ点の位置ベクトルをローカル座標系（ｘｓ,ｙｓ,ｚｓ）に変換し、ステップ１０３で計算され記憶装置５に記憶されている直方体導体４５ａの寸法と、ステップ１０１より入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込み、式（１）〜（４）に基づいて磁場ベクトルを計算する。
【０１０５】
なお、電流密度ベクトルがステップ１１６で更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【０１０６】
次に、制御部３は、式（１）〜（４）で計算された磁場ベクトルをグローバル座標系（ｘ,ｙ,ｚ）に変換し、変換後の磁場ベクトルを記憶装置５に記憶する。
【０１０７】
これにより、通電された直方体導体４５ａが、ｐ点に作る磁場ベクトルが求められる。
【０１０８】
さらに制御部３は、ｐ点の位置ベクトルをローカル座標系（ｘｃ,ｙｃ,ｚｃ）に変換し、さらに円筒座標系（太字の）ｒｐｃ＝（Ｒｐｃ、φｐｃ、Ｚｐｃ）に変換する。
【０１０９】
次に制御部３は、ステップ１０３で計算し記憶装置５に記憶している円弧状柱状導体４５ｃの寸法ｒａ、θ、ｚｂおよびＲｃを読み込み、また、制御部３はステップ１０１で入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込む。
【０１１０】
なお、電流密度ベクトルがステップ１１６で更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【０１１１】
制御部３は式（６）〜（１２）に基づいて磁場ベクトルを計算する。
【０１１２】
次に、計算された磁場ベクトルを直交座標系に変換し、さらにグローバル座標系（ｘ,ｙ,ｚ）に変換し記憶装置５に記憶する。
【０１１３】
これにより、通電された円弧状柱状導体４５ｃがｐ点に作る磁場ベクトルが求められる。
【０１１４】
以上が、図５のステップ１０５の詳細である。
【０１１５】
次に、磁場解析装置１の制御部３は、演算プログラム１９を用いて磁性体を構成する粒子の磁場の運動方程式から束縛を考慮せずに仮想時間刻みδｔ後の磁場ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０６）。
【０１１６】
ここで、図５のステップ１０６について、具体的に説明する。
【０１１７】
磁性体を構成するＮ個の粒子のラグランジアンを式（１３）〜式（１５）で表される形とする。
【０１１８】
【数５】
【０１１９】
ここで、式（１３）において、αは仮想質量、太字のｒは位置ベクトル、太字のＨは磁場ベクトル、太字の傍点付きＨは磁場ベクトルの時間微分、太字のＭは磁化ベクトル、太字のｎは粒子の位置ベクトルを重心とする多面体要素境界面の法線ベクトル、太字のＨｅxｔは外部からの印加磁場ベクトル、ΔＳは粒子の位置ベクトルを重心とする多面体要素の要素境界面の面積、ΔＶは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、χは磁気感受率、μ0は真空の透磁率、λはラグランジュの未定定数、πは円周率、ｓは粒子の位置ベクトルを重心とする多面体要素の要素境界面数である。
【０１２０】
また、各物理量の添え字ｉｐはｉ番目の粒子の位置ベクトルを重心とする多面体要素内でのｐ点の物理量、添え字ｊｑはｊ番目の粒子の位置ベクトルを重心とする多面体要素内でのｑ点の物理量を示す。
【０１２１】
ｐ点、ｑ点はステップ１０４において説明している。
【０１２２】
式（１３）において、ｍは粒子の質量、ｖは速度、φ（ｒｉ−ｒｊ）はｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギである。添え字ｉ、ｊはそれぞれ、ｉ、ｊ番目の粒子の物理量を示す。
【０１２３】
次に、正準変数を（太字の）Ｈｉｐ、（太字の傍点付き）Ｈｉｐとし、式（１３）〜式（１５）で示されるラグランジアンをラグランジュの運動方程式に代入すると、磁場の運動方程式は式（１６）のように記載できる。
【０１２４】
【数６】
【０１２５】
ここで、式（１６）の右辺第２項は、ラグランジュの未定定数を通して束縛（磁化ベクトルの発散は０）を課している。
【０１２６】
式（１６）の右辺第３項は外部からの印加磁場ベクトルが変化したときにすばやく追従させるための減衰項であり、γは減衰定数である。
【０１２７】
式（１６）の右辺第２項に示される束縛を含んだ運動方程式を解くにあたり、本実施形態では、一般化された束縛の導入法であるＳＨＡＫＥ法を採用する。
【０１２８】
束縛を考慮せずに蛙跳び法により式（１６）を離散化すると以下の式（１７）、式（１８）、（１９）になる。
【０１２９】
【数７】
【０１３０】
ここで、δｔは磁化現象の収束計算を行う上で用いる仮想時間刻みである。
【０１３１】
添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ−１／２は（ｎ−１／２）δｔにおける物理量、ｎ＋１／２は（ｎ＋１／２）δｔにおける物理量、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１３２】
磁場解析装置１の制御部３は、ステップ１０４において既に計算され記憶装置５に記憶されている粒子の位置ベクトルを重心とする多面体要素内のｐ点、ｑ点の位置ベクトルを読み込む。
【０１３３】
次に、制御部３は、ステップ１０５で既に計算され記憶装置５に記憶されているコイルが磁性体を構成する粒子を重心とする多面体要素内部のｐ点に作る磁場ベクトルを外部からの印加磁場ベクトルとして読み込む。
【０１３４】
さらに、制御部３は、ステップ１０４で既に計算され、記憶装置５に記憶されている、磁性体を構成する粒子の要素境界面積、法線ベクトルおよび粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
【０１３５】
また、制御部３はあらかじめ記憶装置５に記憶されている減衰定数、仮想質量、仮想時間刻みを読み込む。
【０１３６】
次に、制御部３はあらかじめ記憶装置５に記憶されている磁場ベクトル、磁場ベクトルの時間微分の初期値を読み込む。
【０１３７】
なお、制御部３は後述するステップ１０７において磁場ベクトル、磁場ベクトルの時間微分が更新されている場合は、その値を読み込む。
【０１３８】
なお、磁化ベクトルは、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に磁場ベクトルを代入することで計算される。
【０１３９】
制御部３は、式（１７）、式（１８）、式（１９）を、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内の全てのｐ点に対して計算し、計算された磁場ベクトルを記憶装置５に記憶する。
【０１４０】
以上が、図５のステップ１０６の詳細である。
【０１４１】
次に、磁場解析装置１の制御部３はステップ１０６で計算し記憶装置５に記憶されている、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のｐ点の磁場ベクトルに束縛力を加え、計算された磁場ベクトルを、記憶装置５に記憶する（図５のステップ１０７）。
【０１４２】
ここで、図５のステップ１０７について、具体的に説明する。
【０１４３】
磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のｐ点の磁場ベクトルに以下に示す式（２０）、式（２１）に従って束縛力を加える。
【０１４４】
【数８】
【０１４５】
ここで、太字のＨは磁場ベクトル、αは仮想質量、δｔは仮想時間刻み、γは減衰定数、Ｎは磁性体を構成する粒子数、ｓは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量、ｉｑは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｑ点の物理量、添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１４６】
磁場解析装置１の制御部３は、演算プログラム１９を用いて、ステップ１０６で計算され記憶装置５に記憶されている磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルを読み込む。
【０１４７】
制御部３は、あらかじめ記憶装置５に記憶されている仮想質量、仮想時間刻み、減衰定数を読み込む。
【０１４８】
制御部３は、磁性体を構成する粒子各々に対して式（２０）、式（２１）に基づく計算を行い、計算された磁場ベクトルを記憶装置５に記憶する。
【０１４９】
以上が、図５のステップ１０７の詳細である。
【０１５０】
次に、磁場解析装置１の制御部３はステップ１０７で求めた磁場ベクトルが束縛条件を満たしているかを判断し、満たしていれば次のステップに進み、満たしていなければステップ１０７に戻る（図５のステップ１０８）。
【０１５１】
具体的には制御部３は、ステップ１０７で計算し記憶装置５に記憶されている磁性体を構成する粒子各々の磁場ベクトルから計算される磁化ベクトルを用いて式（２２）に基づく計算を行う。
【０１５２】
【数９】
【０１５３】
ここで、添え字のｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１５４】
ｅｒｒｉは磁性体を構成する粒子の内、ｉ番目の粒子の束縛条件に対する誤差値である。
【０１５５】
磁化ベクトル（太字の）Ｍはステップ１０７において計算され記憶装置５に記憶されている磁場ベクトルを制御部３が読み込み、ステップ１０１より入力情報２１に記憶されている磁化曲線を表す関数に代入することで計算される。
【０１５６】
法線ベクトル（太字の）ｎには、ステップ１０４において既に計算され記憶装置５に記憶されているものを制御部３は読み込み代入する。
【０１５７】
制御部３は、すべての粒子に対して、誤差の値が式（２３）を満たさなければステップ１０７に戻る。
【０１５８】
【数１０】
【０１５９】
式（２３）においてＡ１は任意の誤差判別値であり、あらかじめ記憶装置５に任意の値が記憶されている。
【０１６０】
以上が図５のステップ１０８の詳細である。
【０１６１】
次に、磁場解析装置１の制御部３は、磁性体の磁化現象が定常状態に到達したかを判断し、条件を満たしていれば次のステップに進む。（図５のステップ１０９）
ここで、図５のステップ１０９について、具体的に説明する。
【０１６２】
磁性体を構成する粒子が定常状態に到達したかは、以下の式（２４）により判断される。
【０１６３】
【数１１】
【０１６４】
ここで、太字のＨは磁場ベクトル、μ0は真空の透磁率、Ｎは磁性体を構成する粒子の粒子数、ｓは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点での物理量である。
【０１６５】
さらに、添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ―１は（ｎ―１）δｔ、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１６６】
また、Ａ２は磁性体の磁場ベクトルが定常状態に到達したかを判断するための任意の誤差判定値である。
【０１６７】
磁場解析装置１の制御部３は、あらかじめ記憶装置５に記憶されているＡ２、μ0を読み込む。
【０１６８】
また、制御部３は、ステップ１０７で計算され記憶装置５に記憶されている、磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルを読み込む。
【０１６９】
次に、制御部３は式（２４）を計算し、式（２４）を満たしていれば次のステップに進み、満たしていなければステップ１０６に戻る。
【０１７０】
以上が図５のステップ１０９の詳細である。
【０１７１】
次に、磁場解析装置１の制御部３は、磁性体を構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁化ベクトル、磁束密度ベクトルと、コイルを構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁束密度ベクトルを計算し、記憶装置５に記憶する（図５のステップ１１０）。
【０１７２】
以下にステップ１１０を具体的に説明する。
【０１７３】
（磁性体を構成する粒子）
磁性体を構成する粒子の位置ベクトル上の磁場ベクトルは以下の式（２５）で表される。
【０１７４】
【数１２】
【０１７５】
ここで、（太字の）ｅx＝（１,０,０）、（太字の）ｅｙ＝（０,１,０）、（太字の）ｅｚ＝（０,０,１）である。
【０１７６】
太字のＨは磁場ベクトル、太字のｎは法線ベクトルであり、添え字ｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を表し、添え字ｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量である。
【０１７７】
磁性体を構成する粒子の磁化ベクトルは、磁場ベクトルを磁化曲線を表す関数に代入することで求まる。
【０１７８】
磁性体を構成する粒子の位置ベクトル上での磁束密度ベクトルは、式（２６）で表される。
【０１７９】
【数１３】
【０１８０】
ここで、太字のＢは磁束密度ベクトル、太字のＨは磁場ベクトル、太字のＭは磁化ベクトルであり、添え字のｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を指し示す。
【０１８１】
μ0は真空の透磁率である。
【０１８２】
磁場解析装置１の制御部３は、ステップ１１０までで計算されている、磁性体を構成する粒子を重心とする多面体要素内のｐ点での磁場ベクトルを記憶装置５から読み込む。
【０１８３】
次に、制御部３は、ステップ１０４で既に計算され記憶装置５に記憶されている法線ベクトルを読み込み、式（２５）に基づいて磁性体を構成する粒子の位置ベクトル上での磁場ベクトルを計算し記憶装置５に記憶する。
【０１８４】
さらに、制御部３は、このステップで計算された磁場ベクトルを、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に代入し、磁性体を構成する粒子の磁化ベクトルを計算し記憶装置５に記憶する。
【０１８５】
次に制御部３は、このステップで計算された磁場ベクトルと磁化ベクトルを式（２６）に代入し、磁性体を構成する粒子の位置ベクトル上の磁束密度ベクトルを計算し記憶装置５に記憶する。
【０１８６】
なお、真空の透磁率はあらかじめ記憶装置５に記憶されているものを制御部３が読み込む。
【０１８７】
（コイルを構成する粒子）
コイルを構成する粒子の位置ベクトル上の磁場ベクトルは式（２７）、式（２８）で記載される。
【０１８８】
【数１４】
【０１８９】
ここで、太字のｒは位置ベクトル、太字のＨは磁場ベクトル、Ｎは磁性体を構成する粒子の数、ｓおよびΔＳは磁性体を構成する粒子を重心とする多面体要素の境界面数と境界面積、太字のＭは磁化ベクトル、太字のｎは法線ベクトルであり、添え字のｉはコイルを構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｊは磁性体を構成する粒子の内、ｊ番目の粒子の位置ベクトル上での物理量であり、ｊｑは磁性体を構成する粒子の内、ｊ番目の粒子の位置ベクトルを重心とする多面体要素のｑ点での物理量を指し示す。
【０１９０】
式（２７）の右辺第一項は、磁性体を構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルを記述する項であり、式（２８）で表される。
【０１９１】
式（２７）の右辺第二項は、コイルを構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルであり式（１）〜（１２）で表される。
【０１９２】
磁場解析装置１の制御部３は、ステップ１０３において計算し記憶装置５に記憶されているコイルを構成する粒子の位置ベクトルを記憶装置５より読み込む。
【０１９３】
なお、ステップ１１３においてコイルを構成する粒子の位置ベクトルが更新されている場合は、制御部３はその値を読み込む。
【０１９４】
次に、制御部３はステップ１０２より計算され記憶装置５に記憶されている磁性体を構成する粒子の数と多面体要素の面数、および粒子の位置ベクトルを重心とする多面体要素内のｑ点の位置ベクトルと多面体要素境界面積を読み込む。
【０１９５】
また、制御部３はこのステップで計算され記憶装置５に記憶されている、磁性体を構成する粒子の磁化ベクトルを読み込む。
【０１９６】
さらに、制御部３はステップ１０３において計算され記憶装置５に記憶されているコイルの寸法およびステップ１０１で入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込む。
【０１９７】
なお、ステップ１１６において電流密度ベクトルが更新されている場合は、制御部３はその値を読み込む。
【０１９８】
その後、制御部３は式（２７）、式（２８）および式（１）〜（１２）に従い、コイルを構成する粒子の位置ベクトル上での磁場ベクトルを計算し、記憶装置５に記憶する。
【０１９９】
制御部３は、コイルを構成する粒子の位置ベクトル上での磁束密度ベクトル（太字の）Ｂｉを、すでに計算されているコイルを構成する粒子の位置ベクトル上の磁場ベクトルに真空の透磁率を掛けて計算し記憶装置５に記憶する。
【０２００】
真空の透磁率は、あらかじめ記憶装置５に記憶されている。
【０２０１】
以上が、ステップ１１０の詳細である。
【０２０２】
このように、制御部３は、磁場の運動方程式を解くことにより磁性体の磁化現象を解析する。
【０２０３】
そのため、有限要素法のように空間全域をメッシュ分割する必要がなく、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる。
【０２０４】
また、運動方程式を解くため、行列を扱わず、計算に必要なメモリ量は粒子数に比例する。
【０２０５】
次に、磁場解析装置１の制御部３はＳＰＭモータ３１を構成する全ての粒子に働く力ベクトルを計算し、記憶装置５に記憶する（図５のステップ１１１）。
【０２０６】
以下に図５のステップ１１１を具体的に説明する。
【０２０７】
（磁性体を構成する粒子に働く力ベクトル）
正準変数を太字のｒｉ、太字の傍点付きｒｉとし、式（１３）のラグランジアンをラグラ
ンジュの運動方程式に代入する。
【０２０８】
すると、磁性体を構成する粒子に働く力は式（２９）のように記載される。
【０２０９】
【数１５】
【０２１０】
ここで、μ0は真空の透磁率、太字のＨは磁場ベクトル、太字のＭは磁化ベクトル、太字のｒは位置ベクトル、ΔＶは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、Ｎは磁性体を構成する粒子数であり、添え字のｉ、ｊは磁性体を構成する粒子の内、ｉ番目、ｊ番目の粒子の物理量を指し示し、添え字のｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量を指し示す。
【０２１１】
また、φ（ｒｉ−ｒｊ）はｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギである。
【０２１２】
本実施形態ではＳＰＭモータを剛体（粒子間の相対距離が不変）とする。
【０２１３】
剛体モデルとすると、磁性体を構成する粒子に働く力は式（３０）のように記載される。
【０２１４】
【数１６】
【０２１５】
磁場解析装置１の制御部３は、ステップ１１０までにおいて計算され記憶装置５に記憶されている磁性体を構成する粒子を重心とする多面体要素内のｐ点での磁場ベクトルを読み込む。
【０２１６】
制御部３は、ｐ点での磁場ベクトルを、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に代入することによりｐ点での磁化ベクトルを計算し記憶装置５に記憶する。
【０２１７】
また、制御部３は、ステップ１０４で計算され記憶装置５に記憶されている磁性体を構成する粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
【０２１８】
次に、制御部３は、磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルの偏微分を計算し、その値に基づいて、式（３０）より力ベクトルを計算し記憶装置５に記憶する。
【０２１９】
（コイルを構成する粒子に働く力ベクトル）
コイルを構成する粒子のうち、ｉ番目の粒子の位置ベクトルを太字のｒｉ、その位置ベクトルでの電流の単位ベクトルを太字のｔｉ、磁束密度ベクトルを太字のＢｉ、電流の流れる方向のコイルの長さをＬｉとすると、コイルを構成する粒子の内、ｉ番目の粒子に働く力ベクトルは以下の式（３１）で記載される。
【０２２０】
【数１７】
【０２２１】
磁場解析装置１の制御部３は、ステップ１０３において計算され記憶装置５に記憶されているコイルを構成する粒子の位置ベクトルを読み込む。
【０２２２】
なお、ステップ１１３においてコイルを構成する粒子の位置ベクトルが更新されている場合は制御部３はその値を読み込む。
【０２２３】
また、制御部３はステップ１０３において計算され記憶装置５に記憶されているコイルの寸法を読み込む。
【０２２４】
次に、制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されているコイルの電流密度ベクトルを読み込む。
【０２２５】
なお、制御部３は、ステップ１１６においてコイルの電流密度ベクトルが更新されている場合はその値を用いる。
【０２２６】
制御部３は、式（３１）に基づいて、コイルを構成する粒子に働く力ベクトルを計算し記憶装置５に記憶する。
【０２２７】
以上が、図５のステップ１１１の詳細である。
【０２２８】
次に、磁場解析装置１の制御部３は、ステップ１１１で計算した可動部（移動子）を構成する粒子に働く力ベクトルから、可動部（移動子）の回転量、角速度、並進量、並進速度を計算し、記憶装置５に記憶する（図５のステップ１１２）。
【０２２９】
以下に図５のステップ１１２を具体的に説明する。
【０２３０】
剛体の並進運動の運動方程式は式（３２）で表される。
【０２３１】
【数１８】
【０２３２】
ここで、ｍｔｏｔは可動部（移動子）の質量、太字のｒｃは可動部（移動子）において任意に定めた代表点の位置ベクトル、太字のＦは可動部（移動子）を構成する粒子に働く力ベクトル、Ｎは可動部（移動子）を構成する粒子数である。
【０２３３】
次に、任意の固定点の周りの回転運動の運動方程式は式（３３−１）、式（３３−２）、式（３３−３）で表される。
【０２３４】
【数１９】
【０２３５】
ここで、Ｉ１、Ｉ２、Ｉ３は剛体の主慣性モーメント、ω１、ω２、ω３は慣性主軸座標系での角速度、Ｎ１、Ｎ２、Ｎ３は慣性主軸座標系でのトルクである。
【０２３６】
本実施例のＳＰＭモータ３１のロータ３３はロータシャフト４１を通じてある回転軸に固定されている。
【０２３７】
固定された回転軸をｚ軸とすると、ロータ３３の運動は以下のように計算される。
【０２３８】
ロータ３３はロータシャフト４１を通じてｚ軸に固定されているため、ｚ軸周りの回転運動の運動方程式を解くことで、ロータ３３の運動は計算される。
【０２３９】
ｚ軸周りのロータ３３の回転運動の運動方程式は式（３４）で表される。
【０２４０】
【数２０】
【０２４１】
ここで、Ｉｚはロータ３３の慣性モーメント、ωｚはｚ方向の角速度、Ｎｚはｚ方向のトルクである。
【０２４２】
式（３４）を蛙飛び法で離散化すると、以下の式（３５）、式（３６）を得る。
【０２４３】
【数２１】
【０２４４】
ここで、θｚはｚ軸周りの回転量であり、ｄｔは微小時間刻み幅である。
【０２４５】
また、添え字のｎ−１／２、ｎ、ｎ＋１／２、ｎ＋１はそれぞれ（ｎ―１／２）ｄｔ、ｎｄｔ、（ｎ＋１／２）ｄｔ、（ｎ＋１）ｄｔ秒での物理量に対応している。
【０２４６】
時間刻み幅は記憶装置５にあらかじめ任意の値が記憶されているものを制御部３が読み込む。
【０２４７】
回転量、角速度の初期値は記憶装置５にあらかじめ任意の値が記憶装置５に記憶されている。
【０２４８】
なお、すでにステップ１１２においてロータ３３の回転量、角速度が更新されている場合はその値を用いる。
【０２４９】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されているロータ３３の三次元構造と質量密度により慣性モーメントを計算する。
【０２５０】
さらに、制御部３は、ロータ３３を構成する粒子の位置ベクトルと、ステップ１１１で計算し記憶装置５に記憶されているロータ３３を構成する粒子に働く力ベクトルに基づいて、式（３５）、式（３６）を計算し、ｄｔ秒後のロータ３３の回転量と角速度を計算し記憶装置５に記憶する。
【０２５１】
以上が、図５のステップ１１２の詳細である。
【０２５２】
次に、磁場解析装置１の制御部３はステップ１１２で計算した微小時間後の回転量や重心の移動量を基に、移動子を構成している粒子の位置ベクトルを計算し、記憶する（図５のステップ１１３）。
【０２５３】
次に、磁場解析装置１の制御部３は必要に応じてＳＰＭモータ３１を構成する粒子の位置ベクトル、力ベクトル、磁場ベクトル、磁束密度ベクトル、磁化ベクトルをプリンタポート１１を介してプリンタ１２より出力する（図５のステップ１１４）。
【０２５４】
次に、磁場解析装置１の制御部３は所定の終了条件（時間、移動量等）を満たしているかを判断し、満たしている場合は解析を終了し、満たしていない場合はステップ１１６に進む（図５のステップ１１５）。
【０２５５】
次に、磁場解析装置１の制御部３は必要に応じてコイル４５に流れる電流密度ベクトルを更新し、ステップ１０４に戻る（図５のステップ１１６）。
【０２５６】
以上が本実施形態に係る磁場解析の手順である。
【０２５７】
このように、本実施形態によれば、磁場解析装置１がＳＰＭモータ３１の磁場解析を元に、移動子であるロータ３３の移動量を求めている。
【０２５８】
そのため、本実施形態では磁場と運動の相互作用を演算でき、磁場解析と機構の連成が可能となる。
【０２５９】
また、本実施形態では、コイル４５に流れる電流密度ベクトルを更新しながら解析を行うことができる。
【０２６０】
そのため、コイルに通電する電流を変化させながら解析を行うことができ、従来よりもさらに精度良く磁場解析を行うことができる。
【０２６１】
上記した実施形態では、本発明を磁性体を構成する粒子のラグランジアンから導出される磁場の運動方程式を解くことで磁性体の磁化現象を計算した場合について説明したが、本発明は、何等、これに限定されることなく、運動方程式であればラグランジアンから導出されたもの以外のものであってもよい。
【０２６２】
また、上記した実施形態では本発明をＳＰＭモータ３１の磁場解析に使用したが、本発明は、何等、これに限定されることなく、磁性体が配置された空間における任意の点の磁場を解析するものであれば、例えばリニアモータ等のＳＰＭモータ３１以外のモータ、あるいは磁性体で空間を囲む磁気シールドにおいて、磁性体でシールドされた空間内の磁場解析に用いても良い。
【０２６３】
上述のように本発明の一実施形態によれば、物体に作用する磁場についての磁場解析と、その磁場により生じる物体の運動についての機構解析との連成解析をするための連成解析装置が提供される。すなわち、この連成解析装置は、磁場解析演算部と、機構解析演算部と、を備える。一実施例においては、制御部３は、磁場解析演算部と、機構解析演算部と、を備えてもよい。磁場解析演算部は、物体を構成する各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式に基づき磁場を演算する。機構解析演算部は、演算された磁場に基づき物体の運動を演算し、物体の各粒子の位置を更新する。磁場解析演算部はさらに、更新された各粒子の位置に基づいて磁場を再計算する。このようにして磁場解析と機構解析とを粒子ごとに繰り返すことにより、連成解析が可能である。なお以下では便宜上、磁場解析演算部を磁場演算部とも称し、機構解析演算部を位置演算部とも称する。
【０２６４】
ここで磁場演算部は、上述のように、物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより物体に作用する磁場を演算する。解析対象物体は磁化曲線が非線形である材料を含んでもよく、例えば磁性体を解析対象物体としてもよい。本願発明者は、いわゆる解析力学の手法を用いて、物体の各粒子に作用する磁場の支配方程式を運動方程式の形式で記述することに成功した。上述の一実施例においては、上記の式（１３）に定義されるラグランジアンをラグランジュの運動方程式に代入することにより、式（１６）のように各粒子に作用する磁場を運動方程式の形式で記述している。
【０２６５】
磁場演算部は、磁場の運動方程式を数値的に解くことにより、物体に作用する磁場を求める。具体的には例えば、磁場演算部はまず、離散化された磁場の運動方程式に磁場及びその時間微分の初期値を代入し、微小時間経過後の磁場及びその時間微分を求める。ここで、微小時間とは磁場の収束計算を行うために用いられる仮想的な時間であってもよい。つまりこの仮想時間は実際の経過時間とは異なる。磁場演算部は、演算された微小時間経過後の磁場が収束判定基準を満たすか否かを判定する。収束判定条件は、一実施例においては式（２４）で与えられる。磁場演算部は、収束判定基準が満たされるまで、離散化された磁場の運動方程式を繰り返し演算することにより磁場の数値解を求める。
【０２６６】
この方法によれば、有限要素法に基づく周知の解析法とは異なり逆行列の演算が不要であるので、少ない計算資源で磁場解析を実行することができるという点で好ましい。また、この方法は、高精度の解析結果を、有限要素法のように解析対象物体とその周囲の空間とを含む系全体をメッシュ分割することなく得られるという点でも好ましい。
【０２６７】
一実施例においては、上述の連成解析がモータの磁場解析に利用される。本願発明者の検証によると、モータの磁場解析に要求される高精度のレベルに収束判定条件を設定したときに、演算終了までの所要時間が比較的長くなる場合があることがわかった。この収束の遅れは、解析対象物体における磁気飽和に起因すると考えられる。すなわち、磁化曲線の非線形性が磁場の数値解の収束性に影響している。
【０２６８】
図１０は、磁性体（例えば鉄）の磁化曲線の一例を模式的に示す図である。図１０の横軸は磁場Ｈを示し、縦軸は磁化Ｍを示す。図１０に示されるように、磁場Ｈが小さいときには磁化Ｍは磁場Ｈに比例する。磁場Ｈが大きくなるにつれて、磁化Ｍと磁場Ｈとの比例定数（いわゆる磁気感受率）は徐々に小さくなっていく。磁場Ｈが小さい場合の比例定数は（例えば図中の点Ａ）、磁場Ｈが大きい場合（図中の点Ｂ）の比例定数の１０００倍程度の大きさとなり得る。このように、磁性体の磁化曲線は非線形である。
【０２６９】
磁場の運動方程式には、磁場及びその時間微分項だけではなく、磁化に依存する磁化依存項が含まれる。ある粒子についての磁場の運動方程式における磁化依存項は例えば、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項である。このような相互作用を表す項は、例えば式（１４）に示されるように、磁化Ｍを陽に含む。したがって、式（１６）に示される磁場の運動方程式は、収束計算を進めるにつれて非線形に変化する項を含むと言える。収束計算の磁場の初期値は通常小さい値（例えばゼロ）にとるから、式１４のＨ（ｒｉｐ）は、収束計算の当初に比べて完了時には１０００分の１程度の小さい値となり得る。方程式におけるこの大きな非線形性によって数値解の収束に遅れが生じると考えられる。
【０２７０】
磁場の運動方程式には、任意に設定することが許容される調整可能パラメタを導入することができる。例えば式（１６）に示される磁場の運動方程式においては、仮想質量α及び減衰定数γは調整可能パラメタである。パラメタの設定作業を簡単にすることを重視する場合には、これらの調整可能パラメタを一定値に設定してもよい。この場合、収束計算中に調整可能パラメタは不変である。一方、磁場演算部は、収束計算中に調整可能パラメタを変化させてもよい。
【０２７１】
磁場の運動方程式は一般に、束縛を考慮に入れなければ、２階の線形常微分方程式の形式をとる。この形式は、１自由度の振動の基礎方程式と共通である。よって、磁場の運動方程式の収束特性は、振動の基礎方程式の減衰特性に対応する。このことから、振動の方程式の各項の係数に応じて減衰特性が変化することに準じて、磁場の運動方程式の調整可能パラメタを変更することにより所望の収束特性に調整することが可能である。
【０２７２】
そこで、本発明の一実施形態においては、磁場演算部は、磁場の運動方程式に導入されている調整可能パラメタを、磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整してもよい。磁場演算部は、磁化依存項による数値解の収束性への影響を抑えるように調整可能パラメタを変化させてもよい。一実施例においては、磁場演算部は、収束判定基準を満たすか否かを判定する収束判定処理において収束判定基準が満たされていないと判定したときに、調整可能パラメタを変化させる係数調整処理を実行してもよい。上述の実施例においては、制御部３は、誤差の値が誤差判別値Ａ２を超えると判定した場合に（図５のステップ１０９のＮｏ）、仮想質量α及び減衰定数γを調整する係数調整ステップを実行してもよい。
【０２７３】
ある物体の粒子に作用する磁場Ｈの運動方程式は一般に、式（３７）で表される。ここで、係数ＣＭは磁場Ｈの２階時間微分項の係数を表し、係数ＣＤは磁場Ｈの１階時間微分項の係数を表し、係数ＣＫは磁場Ｈの項の係数を表す。係数ＣＭ、ＣＤ、ＣＫはいずれも解析者が調整することのできるパラメタである。一方、係数ＦＮＬは、磁場Ｈと磁化Ｍとの比例定数に依存する。Ｈ０は外部からの印加磁場を表す。
【０２７４】
【数２２】
【０２７５】
よって、調整可能なパラメタである係数ＣＭ、ＣＤを磁化曲線に基づいて調整することにより、磁化曲線の非線形性が数値解の収束性に与える影響を排除することができる。すなわち、磁場Ｈの時間微分項の係数ＣＭ及びＣＤを、式（３８）に示されるように磁化依存係数ＦＮＬを用いてＣＭ’ＦＮＬ及びＣＤ’ＦＮＬと定義する。磁化依存係数ＦＮＬの時間変化率が十分に小さければ、磁化曲線の非線形性が数値解の収束性に与える影響を事実上排除することができる。
【０２７６】
【数２３】
【０２７７】
式（３８）の磁場の項及びその時間微分項には、磁場Ｈに対し非線形に変化する共通パラメタＦＮＬが導入されている。つまり、共通の磁化依存係数ＦＮＬをもたない運動方程式（３９）の収束特性に、式（３８）に示される磁場の運動方程式の収束特性は実質的に等しくなる。
【０２７８】
式（３９）の磁場Ｈ及びその時間微分項の係数ＣＭ’、ＣＤ’、及びＣＫはいずれも磁化Ｍに依存せず、解析者が任意に設定することのできるパラメタである。よって、式（３８）において磁化依存係数ＦＮＬが毎回の収束計算ごとに変化しても、係数ＣＭ’、ＣＤ’、及びＣＫから定まる一定の収束特性を得ることができる。なおこの場合、所望の収束特性をもたらすよう係数ＣＭ’、ＣＤ’、及びＣＫを予め設定することが好ましい。
【０２７９】
なお本発明は、磁化曲線の非線形性の影響を磁場の運動方程式の収束特性から排除する実施例には限られない。基礎方程式の変数の係数がその変数に対して非線形に変化する場合における運動方程式型の数値解法にも適用可能である。例えばいわゆる大変形解析にも適用可能であろう。
【０２８０】
一実施例においては、第１の物理量と、該第１の物理量と非線形の関係をもつ第２の物理量とを有する物体の特性を解析するための解析装置が提供される。例えば上述の実施例では第１の物理量及び第２の物理量はそれぞれ磁場及び磁化に対応する。この解析装置は、運動方程式の形式で記述された第１の物理量が満たすべき関係式を数値的に解く数値演算部を備える。この関係式は、第２の物理量に依存する項と、第２の物理量に依存しない項と、を含み、第２の物理量に依存しない項の係数が第２の物理量に依存する項に連動するよう定義されている。一実施例においては、第１の物理量及びその時間微分項のそれぞれの係数には第２の物理量に連動する連動部分が共通に導入されていてもよい。
【０２８１】
また、一実施例においては、調整可能パラメタは、磁場の運動方程式を数値的に解く過程で導出される値を用いて定義されていることが好ましい。つまり、磁場演算部は、磁場演算処理において算出された値を用いて係数調整処理を実行してもよい。このようにすれば、磁場の運動方程式を解く過程で得られた計算済みの値を流用して調整可能パラメタを調整することができる。調整可能パラメタの設定のために新たな演算をする必要がなくなるので、計算コストを低減することができる。よって、収束演算の所要時間を短くすることもできる。
【０２８２】
そこで、一実施例においては、調整可能パラメタは、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項に連動して調整されてもよい。式（１６）を解くことにより磁場を求める上述の実施例においては例えば、仮想質量αｉｐを式（４０）及び式（４１）により可変量として定義する。係数ｋｉｐは磁場の運動方程式に存在するＨ（ｒｉｐ）の項を用いて定義されているから、調整可能パラメタの設定のために新たな演算をする必要がない。なお係数α０は解析者が任意に設定することのできるパラメタである。
【０２８３】
【数２４】
【０２８４】
この仮想質量αｉｐにより、磁化曲線の非線形性の影響が磁場の運動方程式の収束特性から排除されることを説明する。上述のようにＳＨＡＫＥ法を適用して式（１６）を解くことを考慮すると、式（１６）の右辺第２項は０に収束するから、式（４２）のように表される。
【０２８５】
【数２５】
【０２８６】
式（１４）及び式（４１）から、Ｈ（ｒｉｐ）は式（４３）で表される。
【０２８７】
【数２６】
【０２８８】
式（４０）及び式（４３）を式（４２）に代入することにより、磁場の運動方程式は式（４４）で表される。式（４４）の磁場の項及びその時間微分項には、磁場Ｈに対し非線形に変化する共通パラメタ（１−ｋｉｐ）が導入されている。よって、磁場の運動方程式（４４）の収束特性は、共通パラメタ（１−ｋｉｐ）を有しない式（４５）の収束特性に実質的に等しくなる。よって、磁気飽和に起因して係数ｋｉｐ（式（１６）で言えばＨ（ｒｉｐ））が収束計算中に変化しても、磁場の運動方程式（４４）の収束特性は変化しない。係数α０、γを予め設定することにより、所望の収束特性を実現することができる。
【０２８９】
【数２７】
【０２９０】
なお、上述のパラメタの更新は、磁場解析の対象物体が磁気感受率の異なる物質を含む場合にも有効であると考えられる。例えば物体が磁石と鉄とを含む場合にも有効である。物質が異なることによる粒子ごとの磁気感受率の違いは、磁気飽和が生じているか否かによる磁気感受率の違いと同様に上記の磁場の運動方程式において取り扱われているからである。
【０２９１】
図１１は、本発明の一実施形態に係る磁場演算処理を説明するためのフローチャートである。図１１に示されるように、磁場解析装置１の制御部３はまず、磁場の運動方程式を作成する（Ｓ１０）。一実施例においては、制御部３は、図５に示されるステップ１０１乃至１０５を実行することにより磁場の運動方程式（１６）を作成する。ここで、仮想質量αｉｐは式（４０）及び式（４１）により可変量として定義されている。
【０２９２】
次に制御部３は、磁場の運動方程式に基づいて磁場を演算する（Ｓ１１）。一実施例においては、制御部３は、図５に示されるステップ１０６乃至１０８を実行することにより、磁場の運動方程式（１６）に基づいて磁場を演算する。制御部３は、演算された磁場が収束判定基準を満たすか否かを判定する（Ｓ１２）。収束判定基準は例えば図５に示すステップ１０９と同様であり、制御部３は、誤差が誤差判別値を満たすか否かを判定する。
【０２９３】
制御部３は、収束判定基準が満たされていないと判定した場合には（Ｓ１２のＮｏ）、調整可能パラメタを更新して（Ｓ１３）、磁場の運動方程式に基づく磁場演算を再度行う（Ｓ１１）。制御部３は、収束判定基準が満たされるまで、調整可能パラメタの更新と磁場演算とを繰り返す。制御部３は、調整可能パラメタとして仮想質量αｉｐを更新する。すなわち、制御部３は係数ｋｉｐを更新する。
【０２９４】
制御部３は、収束判定基準が満たされたと判定した場合には（Ｓ１２のＹｅｓ）、演算された磁場に基づいて解析対象物体の運動を演算する（Ｓ１４）。一実施例においては、制御部３は、図５に示されるステップ１１０乃至１１４を実行する。次に制御部３は、所定の終了条件が満たされているか否かを判定し（Ｓ１５）、満たされていると判定した場合には解析を終了する（Ｓ１５のＹｅｓ）。制御部３は、終了条件が満たされていないと判定した場合には（Ｓ１５のＮｏ）、運動による解析対象物体の新たな位置に基づいて再度磁場を演算する。このようにして、制御部３は、物体に作用する磁場とそれによる物体の運動との連成解析を実行する。
【０２９５】
図１２及び図１３は、本発明の一実施形態による磁場演算の例を示す図である。図１２は仮想質量αｉｐを一定値に固定した場合の磁場Ｈと仮想時間δｔとの関係を示し、図１３は上述のように係数ｋｉｐを用いて仮想質量αｉｐを更新しながら収束計算をしたときの磁場Ｈと仮想時間δｔとの関係を示す。図１２及び図１３の磁場演算の例は、１次元に配列された１０個の粒子からなる物体に対して行ったものである。粒子１及び粒子１０が物体の両端である。図１２の演算結果に比べて図１３の演算結果は磁場が短時間で定常状態に達していることがわかる。仮想質量αｉｐを更新しながら収束計算を行うことにより、磁場の運動方程式の数値解の磁気飽和に伴う収束の遅れを低減することができる。
【符号の説明】
【０２９６】
１ 磁場解析装置、 ３ 制御部、 ５ 記憶装置、 ７ 入力部、 ９ 表示部、 １１ プリンタポート、 １２ プリンタ、 １３ バス、 ３１ ＳＰＭモータ、 ３３ ロータ、 ３５ ステータ、 ３７ ロータコア、 ３９ 永久磁石、 ４１ ロータシャフト、 ４３ ステータティース、 ４４ コアバック、 ４５ コイル、 ４６ フレーム。
【技術分野】
【０００１】
本発明は、磁場解析装置および磁場解析方法に関する。
【背景技術】
【０００２】
例えば特許文献１には、磁場解析を実行する演算処理装置を具えるモータ解析装置が記載されている。演算処理装置は、ユーザの操作に基づく外部指令に応じて、有限要素法による静磁場解析やマクスウェルの応力法によるトルク計算を実行する。有限要素法による静磁場解析のためにメッシュ分割が行われる。メッシュ分割は、コア領域及びハウジング領域、さらには外部空気層領域についても行われる。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平１１−１４６６８８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかし、静磁場解析は時間的な変動を考慮しないから、駆動状態のモータの電磁気的挙動を必ずしも十分に反映することができない。その結果、解析の精度も不十分となりがちである。また、演算の所要時間は短いことが実用上好ましい。
【０００５】
そこで、本発明は、高速かつ高精度に解析結果を得ることのできる磁場解析装置及び磁場解析方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明の一態様によれば、物体を含む系における磁場を解析するための磁場解析装置が提供される。この装置は、物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより前記物体に作用する磁場を演算する磁場演算部を備える。前記磁場の運動方程式には調整可能なパラメタが導入されており、前記磁場演算部は、前記物体の磁化曲線の非線形性が数値解の収束性に与える影響を緩和するよう前記パラメタを調整する。
【０００７】
この態様によれば、運動方程式を数値的に解くので、有限要素法に基づく周知の磁場解析法とは異なり逆行列の演算が不要である。このため、少ない計算資源で精度よく磁場解析を実行することができる。また、磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整可能パラメタが調整されるので、磁気飽和に伴う数値解の収束遅れを低減することができる。
【０００８】
本発明の他の態様は、磁場解析方法である。この方法は、磁化曲線が非線形である材料を含む物体に作用する磁場が満たすべき関係を運動方程式の形式で記述した磁場の運動方程式を数値的に解く磁場演算ステップと、前記磁場の運動方程式に導入されている調整可能パラメタを、前記磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整する調整ステップと、備える。
【０００９】
本発明の一実施例は、この磁場解析方法をコンピュータに実行させるためのプログラムであってもよいし、このプログラムを記録したコンピュータ読み取り可能な記録媒体であってもよい。
【００１０】
本発明の他の態様は、連成解析装置である。この装置は、物体を構成する各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式に基づき磁場を演算する磁場解析演算部と、演算された磁場に基づき前記物体の運動を演算する機構解析演算部と、を備える。
【００１１】
本発明の他の態様は、第１の物理量と、該第１の物理量と非線形の関係をもつ第２の物理量とを有する物体の特性を解析するための解析装置である。この装置は、運動方程式の形式で記述された第１の物理量が満たすべき関係式を数値的に解く数値演算部を備え、前記関係式は、第２の物理量に依存する項と、第２の物理量に依存しない項と、を含み、第２の物理量に依存しない項の係数が第２の物理量に依存する項に連動するよう定義されている。
【発明の効果】
【００１２】
本発明によれば、高速かつ高精度に磁場解析を行うことができる。
【図面の簡単な説明】
【００１３】
【図１】本発明の一実施形態に係る磁場解析装置のハードウェア構成を示す図である。
【図２】本発明の一実施形態に係る磁場解析装置の記憶装置を示す図である。
【図３】本発明の一実施形態に係るＳＰＭモータの構成の概略を示す図である。
【図４】図３のＳＰＭモータのステータティースの１つの周囲の拡大斜視図である。
【図５】本発明の一実施形態に係る磁場解析装置を用いた磁場の解析の手順を示すフローチャートである。
【図６】図６（ａ）はコイルの拡大斜視図であって、図６（ｂ）はコイルをローカル導体に分割した例を示す図である。
【図７】本発明の一実施形態に係る直方体導体とローカル座標の関係を示す図である。
【図８】本発明の一実施形態に係る円弧状柱状導体とローカル座標の関係を示す図である。
【図９】本発明の一実施形態に係る磁性体の磁場計算に用いる要素の説明図である。
【図１０】磁性体（例えば鉄）の磁化曲線の一例を模式的に示す図である。
【図１１】本発明の一実施形態に係る磁場演算処理を説明するためのフローチャートである。
【図１２】本発明の一実施形態による磁場演算の例を示す図である。
【図１３】本発明の一実施形態による磁場演算の例を示す図である。
【発明を実施するための形態】
【００１４】
以下、図面に基づいて本発明に好適な実施形態を詳細に説明する。
【００１５】
まず、図１を参照して、本実施形態に係る磁場解析装置１のハードウェア構成を説明する。
【００１６】
なお、図１のハードウェア構成は例示であり、これに限定されないのは当然である。
【００１７】
図１に示すように、磁場解析装置１は制御部３、記憶装置５、メディア入出力部６、入力部７、表示部９、プリンタポート１１等がバス１３を介して互いに接続されている。
【００１８】
制御部３は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等で構成され、記憶手段としての記憶装置５に格納されたプログラムに従って、バス１３を介して接続された各装置を駆動制御する。
【００１９】
図２に示すように、記憶装置５には、磁場解析装置１の各構成部分を駆動制御するための制御プログラム１５、本発明を実施するための磁場解析プログラム１７が格納されている。
【００２０】
磁場解析プログラム１７は、解析条件を有する情報である入力情報２１と、入力情報２１と磁場の運動方程式に基づいて磁場を演算し力学系の運動方程式を演算する演算プログラム１９とを有している。
【００２１】
メディア入出力部６は、フロッピー（登録商標）ディスク、ＣＤ、ＤＶＤ等のメディアとの間で情報の入出力を行う装置である。
【００２２】
入力部７は、キーボード、マウス等の入力装置であり、表示部９はディスプレイ等の表示機器である。
【００２３】
プリンタポート１１には出力装置としてのプリンタ１２等が接続される。
【００２４】
次に、磁場解析装置１を用いた磁場の解析の手順について図３〜９を参照して説明する。
【００２５】
ここでは永久磁石モータ（Permanent Magnet Motor）の一種である、ＳＰＭモータ３１（Surface Permanent Magnet Motor）の磁場解析を例にして説明する。
【００２６】
まず、ＳＰＭモータ３１の構成の概略を図３および図４を参照して説明する。
【００２７】
図３に示すように、ＳＰＭモータ３１は回転子（移動子）であるロータ３３と固定子であるステータ３５を有している。
【００２８】
ロータ３３は鉄等の磁性体である円柱状のロータコア３７を有し、ロータコア３７の表面には永久磁石３９が設けられている。
【００２９】
ロータコア３７の軸中心には棒状のロータシャフト４１が設けられている。
【００３０】
ステータ３５は磁性体である歯状のステータティース４３とステータティース４３の外側に設けられた円筒状の磁性体であるコアバック４４、コアバック４４の外側に設けられた円筒状のフレーム４６から構成されている。
【００３１】
図３および図４に示すように、ステータティース４３には、金属等の導電体であるコイル４５が巻きつけられている。
【００３２】
なお、実際のＳＰＭモータ３１ではコイル４５は仕様に応じたターン数でステータティース４３に巻きつけられて束となっているが、本実施形態では、図３および図４に描かれているように、コイル一本一本をモデル化せず、コイルの束を一つの導体として扱う。
【００３３】
このような構造のＳＰＭモータ３１は、永久磁石３９の磁場、およびコイル４５に電流を流すことにより発生する磁場によって、ロータ３３、ステータ３５が磁化する。磁性体の磁気エネルギーの偏差によりＳＰＭモータ３１は駆動する。
【００３４】
そのため、ＳＰＭモータ３１の磁場解析を行うためにはコイル４５、永久磁石３９がロータ３３、ステータ３５を構成する磁性体上に作る磁場ベクトルを計算し、これら磁性体の磁化現象を解析する必要がある。
【００３５】
次に、解析の手順について図５〜図９を参照して説明する。
【００３６】
なお、以下の手順においてはＳＰＭモータ３１を、複数に要素分割して要素ごとの粒子の集合体とし、剛体モデルとして扱っているが、本発明はこれに限定されることはなく、粒子を用いずに要素ごとの剛体モデルとして扱ってもよい。
【００３７】
また、以下の手順において、磁性体とは、ロータ３３、ステータ３５を構成する磁性体と永久磁石３９を指し示し、磁化曲線を表す関数によってこれらは区別される。
【００３８】
まず、磁場解析装置１の制御部３は磁場解析プログラム１７を起動し、解析したいＳＰＭモータ３１の解析条件としての三次元構造（形状、座標点）、質量密度、磁性体の磁化曲線を表す関数、導体の電流密度ベクトルを記憶装置５の入力情報２１として記憶する（図５のステップ１０１）。
【００３９】
これらの物理量は例えばメディア入出力部６を介してＣＤ−ＲＯＭ等の記録媒体から読込んだものであってもよい。
【００４０】
また、ＳＰＭモータ３１の三次元構造の情報とは例えばＣＡＤ等のデータである。
【００４１】
さらに、あらかじめ上記物理量が入力情報２１として記憶されている場合は、上記ステップは不要である。
【００４２】
電流密度ベクトルは、コイル４５の作る磁場ベクトルを計算する際に必要になる。
【００４３】
さらに、磁性体の磁化曲線を表す関数は磁化ベクトル（太字の）Ｍを計算する際に必要になる。
【００４４】
以上が図５のステップ１０１の詳細である。
【００４５】
次に、磁場解析装置１の制御部３は、入力情報２１の有する三次元構造の情報から、磁性体をＮ個の粒子（粒子の位置ベクトルを重心とする多面体要素（本実施例では立方体要素））に分割し、粒子の位置ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０２）。
【００４６】
ここで、Ｎは任意の整数であり、粒子の数Ｎおよび位置ベクトルは入力情報２１の有する三次元構造およびあらかじめ記憶装置５に記憶されている多面体要素の形状により計算される。
【００４７】
以上が図５のステップ１０２の詳細である。
【００４８】
次に、制御部３は、入力情報２１の有する三次元構造の情報から、導体（図６（ａ）に示すコイル４５）を図６（ｂ）に示すようにローカル導体（直方体導体４５ａ、４５ｂと円弧状柱状導体４５ｃ、４５ｄ）に分割し、それぞれの導体が作る磁場ベクトルを計算するための係数を計算し、記憶装置５に記憶する（図５のステップ１０３）。
【００４９】
ここで、図５のステップ１０３について、具体的に説明する。
【００５０】
まず、ローカル導体が直方体導体４５ａの場合について説明する。
【００５１】
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。
【００５２】
ローカル導体が直方体導体４５ａの場合は、制御部３は、図７に示すように、ローカル座標系（ｘｓ,ｙｓ,ｚｓ）を適用する。
【００５３】
このローカル座標系においてはローカル導体（直方体導体４５ａ）の重心を原点Ｏｓとし、直方体導体４５ａの寸法はｘｓ方向に２ａ、ｙｓ方向に２ｂ、ｚｓ方向に２ｃの長さを持つものとする。
【００５４】
また原点Ｏｓに粒子は位置するものとする。
【００５５】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体（直方体導体４５ａ）の寸法であるａ，ｂ，ｃと粒子位置ベクトルを計算し記憶装置５に記憶する。
【００５６】
次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
【００５７】
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。
【００５８】
ローカル導体が円弧状柱状導体４５ｃの場合は、制御部３は、図８に示すようにローカル座標系（ｘｃ,ｙｃ,ｚｃ）を適用する。
【００５９】
このローカル座標系においては原点Ｏｃは円弧の中心軸上に存在し、かつ円弧状柱状導体４５ｃの高さ方向（図８のzｃ方向）に対して円弧状柱状導体４５ｃが対称となる点に存在するものとする。
【００６０】
また、ｘｃ,ｙｃ,ｚｃは、ｘｃ−ｙｃ平面でみると、＋ｘｃ軸を基点とし、円弧状柱状導体４５ｃの円弧が＋ｚｃ軸からみて反時計回りになるようして決定する。
【００６１】
円弧状柱状導体４５ｃの内径と外径の平均値をＲｃとし、径方向の厚さを２ｒａ、ｚｃ方向の高さを２ｚｂとする。
【００６２】
電流は＋ｘｃ軸を基点として、＋ｚｃから見て反時計回りの方向への角度をθとし、電流はこの方向に一様な電流密度ｊで流れているものとする。
【００６３】
粒子は円筒座標系で（Ｒｃ、θ／２、０）に位置するものとする。
【００６４】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体（円弧状柱状導体４５ｃ）の寸法であるｒａ、ｚｂ、θ、Ｒｃおよび粒子位置ベクトルを計算し記憶装置５に記憶する。
【００６５】
さらに、あらかじめ上記物理量が入力情報２１として記憶されている場合は、上記ステップは不要である。
【００６６】
以上が図５のステップ１０３の詳細である。
【００６７】
次に、磁場解析装置１の制御部３は、このステップまでに計算され記憶装置５に記憶されている磁性体を構成するＮ個の粒子の位置ベクトルと、あらかじめ記憶装置５に記憶されている多面体要素の形状により、磁性体を構成する粒子の磁場の運動方程式の係数を計算し記憶装置５に記憶する（図５のステップ１０４）。
【００６８】
ここで、図５のステップ１０４について、具体的に説明する。
【００６９】
分割した多面体要素（本実施形態では立方体要素）を２次元表示すると図９に示す形状となる。
【００７０】
ここで、粒子の位置ベクトルを（太字の）ｒｇ、粒子の位置ベクトルを重心とする多面体要素の要素境界面の中点をｑ点とし、粒子の位置ベクトルとｑ点との中間点にｐ点を定義する。
【００７１】
制御部３は、ステップ１０２において計算し記憶装置５に記憶している磁性体を構成する粒子（粒子を重心とする多面体要素（本実施形態では立方体要素））の位置ベクトルを読み込む。
【００７２】
なお、ステップ１１３において粒子の位置ベクトルが更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【００７３】
制御部３は、粒子の位置ベクトルとあらかじめ記憶装置５に記憶されている多面体要素の形状から、磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内のｐ点とｑ点の位置ベクトルと要素境界面への法線ベクトル（太字の）ｎを計算し、記憶装置５に記憶する。
【００７４】
次に、制御部３は、磁性体を構成する粒子を重心とする多面体要素の境界面積ΔＳ、および粒子の位置ベクトルを重心とする多面体要素の境界面の頂点と粒子位置を頂点とする多面体要素体積ΔＶを計算し、記憶装置５に記憶する。
【００７５】
以上が図５のステップ１０４の詳細である。
【００７６】
次に、磁場解析装置１の制御部３は、ステップ１０３で記憶装置５に記憶されているコイルの寸法およびこのステップまでに記憶装置５に記憶されている電流密度ベクトルを用いて、ビオ・サバールの法則を積分することにより得られる解析解により、通電されたコイル４５が、磁性体を構成する粒子上に作る磁場ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０５）。
【００７７】
ここで、図５のステップ１０５について、具体的に説明する。
【００７８】
ここでは、任意の位置ベクトルが図４における磁性体としてのステータティース４３のある点Ｐであると仮定した場合に、通電されたコイル４５がＰ点の位置ベクトル上に作る磁場ベクトルを計算する手順を例に説明する。
【００７９】
まず、ローカル導体が直方体導体４５ａの場合について説明する。
【００８０】
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。
【００８１】
ローカル導体が直方体導体４５ａの場合は、図７に示すようにローカル座標系（ｘｓ,ｙｓ,ｚｓ）を適用する。
【００８２】
このローカル座標系はステップ１０３で説明している。
【００８３】
次に、Ｐ点の位置ベクトルをローカル座標系（ｘｓ,ｙｓ,ｚｓ）に変換する。
【００８４】
通電された直方体導体がＰ点に作る磁場ベクトルは、以下に示す式（１）〜（３）で記載される。
【００８５】
【数１】
【００８６】
ここで、（太字の）ｒｐｓはローカル座標系（ｘｓ,ｙｓ,ｚｓ）でのＰ点の位置ベクトルであり、ｘｐｓ,ｙｐｓ,ｚｐｓはｘｓ,ｙｓ,ｚｓ方向の値である。
πは円周率である。
【００８７】
Ｈｘｓ、Ｈｙｓ、Ｈｚｓはローカル座標系（ｘｓ,ｙｓ,ｚｓ）における磁場ベクトルの各成分である。
ｊは電流密度である。
【００８８】
また、ｘｉ,ｙｊ,ｚｋはｘｓ,ｙｓ,ｚｓ方向の積分の上限、下限を表しており、式（４）に示す関係が成立する。
【００８９】
【数２】
【００９０】
ここで、ａ，ｂ，ｃは直方体導体の寸法である。
【００９１】
次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
【００９２】
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。
【００９３】
ローカル導体が円弧状柱状導体４５ｃの場合は、図８に示すようにローカル座標系（ｘｃ,ｙｃ,ｚｃ）を適用する。
【００９４】
このローカル座標系はステップ１０３において説明している。
【００９５】
ローカル座標系（ｘｃ,ｙｃ,ｚｃ）に変換後のＰ点の位置ベクトルを以下の式（５）に示すように円筒座標系（太字の）ｒｐｃ＝（Ｒｐｃ、φｐｃ、Ｚｐｃ）に変換する。
【００９６】
通電された円弧状柱状導体４５ｃがＰ点に作る磁場ベクトルは、以下の式（６）〜（１１）で表される。
【００９７】
【数３】
【００９８】
ここで、Ｈｒｃ、Ｈｔｃ、Ｈｚｃは円筒座標系での磁場ベクトルの各成分である。
ｊは電流密度である。
【００９９】
ｒａ、θ、ｚｂは円弧状柱状導体４５ｃの寸法であり、Ｒｃは円弧の内径と外径の平均値である。
【０１００】
ｓｇｎはＺｋの符号であり、Ｒｊ、Ｚｋは積分の上限、下限を表しており、式（１２）に示す関係が成立する。
【０１０１】
【数４】
【０１０２】
以上が任意の点Ｐにコイルが作る磁場ベクトルの計算手順である。
【０１０３】
磁場解析装置１の制御部３は、ステップ１０４で既に計算され記憶装置５に記憶されている磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内部のｐ点の位置ベクトルを読み込む。
【０１０４】
次に、制御部３は、ｐ点の位置ベクトルをローカル座標系（ｘｓ,ｙｓ,ｚｓ）に変換し、ステップ１０３で計算され記憶装置５に記憶されている直方体導体４５ａの寸法と、ステップ１０１より入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込み、式（１）〜（４）に基づいて磁場ベクトルを計算する。
【０１０５】
なお、電流密度ベクトルがステップ１１６で更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【０１０６】
次に、制御部３は、式（１）〜（４）で計算された磁場ベクトルをグローバル座標系（ｘ,ｙ,ｚ）に変換し、変換後の磁場ベクトルを記憶装置５に記憶する。
【０１０７】
これにより、通電された直方体導体４５ａが、ｐ点に作る磁場ベクトルが求められる。
【０１０８】
さらに制御部３は、ｐ点の位置ベクトルをローカル座標系（ｘｃ,ｙｃ,ｚｃ）に変換し、さらに円筒座標系（太字の）ｒｐｃ＝（Ｒｐｃ、φｐｃ、Ｚｐｃ）に変換する。
【０１０９】
次に制御部３は、ステップ１０３で計算し記憶装置５に記憶している円弧状柱状導体４５ｃの寸法ｒａ、θ、ｚｂおよびＲｃを読み込み、また、制御部３はステップ１０１で入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込む。
【０１１０】
なお、電流密度ベクトルがステップ１１６で更新され記憶装置５に記憶されている場合は、制御部３はその値を読み込む。
【０１１１】
制御部３は式（６）〜（１２）に基づいて磁場ベクトルを計算する。
【０１１２】
次に、計算された磁場ベクトルを直交座標系に変換し、さらにグローバル座標系（ｘ,ｙ,ｚ）に変換し記憶装置５に記憶する。
【０１１３】
これにより、通電された円弧状柱状導体４５ｃがｐ点に作る磁場ベクトルが求められる。
【０１１４】
以上が、図５のステップ１０５の詳細である。
【０１１５】
次に、磁場解析装置１の制御部３は、演算プログラム１９を用いて磁性体を構成する粒子の磁場の運動方程式から束縛を考慮せずに仮想時間刻みδｔ後の磁場ベクトルを計算し、記憶装置５に記憶する（図５のステップ１０６）。
【０１１６】
ここで、図５のステップ１０６について、具体的に説明する。
【０１１７】
磁性体を構成するＮ個の粒子のラグランジアンを式（１３）〜式（１５）で表される形とする。
【０１１８】
【数５】
【０１１９】
ここで、式（１３）において、αは仮想質量、太字のｒは位置ベクトル、太字のＨは磁場ベクトル、太字の傍点付きＨは磁場ベクトルの時間微分、太字のＭは磁化ベクトル、太字のｎは粒子の位置ベクトルを重心とする多面体要素境界面の法線ベクトル、太字のＨｅxｔは外部からの印加磁場ベクトル、ΔＳは粒子の位置ベクトルを重心とする多面体要素の要素境界面の面積、ΔＶは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、χは磁気感受率、μ0は真空の透磁率、λはラグランジュの未定定数、πは円周率、ｓは粒子の位置ベクトルを重心とする多面体要素の要素境界面数である。
【０１２０】
また、各物理量の添え字ｉｐはｉ番目の粒子の位置ベクトルを重心とする多面体要素内でのｐ点の物理量、添え字ｊｑはｊ番目の粒子の位置ベクトルを重心とする多面体要素内でのｑ点の物理量を示す。
【０１２１】
ｐ点、ｑ点はステップ１０４において説明している。
【０１２２】
式（１３）において、ｍは粒子の質量、ｖは速度、φ（ｒｉ−ｒｊ）はｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギである。添え字ｉ、ｊはそれぞれ、ｉ、ｊ番目の粒子の物理量を示す。
【０１２３】
次に、正準変数を（太字の）Ｈｉｐ、（太字の傍点付き）Ｈｉｐとし、式（１３）〜式（１５）で示されるラグランジアンをラグランジュの運動方程式に代入すると、磁場の運動方程式は式（１６）のように記載できる。
【０１２４】
【数６】
【０１２５】
ここで、式（１６）の右辺第２項は、ラグランジュの未定定数を通して束縛（磁化ベクトルの発散は０）を課している。
【０１２６】
式（１６）の右辺第３項は外部からの印加磁場ベクトルが変化したときにすばやく追従させるための減衰項であり、γは減衰定数である。
【０１２７】
式（１６）の右辺第２項に示される束縛を含んだ運動方程式を解くにあたり、本実施形態では、一般化された束縛の導入法であるＳＨＡＫＥ法を採用する。
【０１２８】
束縛を考慮せずに蛙跳び法により式（１６）を離散化すると以下の式（１７）、式（１８）、（１９）になる。
【０１２９】
【数７】
【０１３０】
ここで、δｔは磁化現象の収束計算を行う上で用いる仮想時間刻みである。
【０１３１】
添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ−１／２は（ｎ−１／２）δｔにおける物理量、ｎ＋１／２は（ｎ＋１／２）δｔにおける物理量、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１３２】
磁場解析装置１の制御部３は、ステップ１０４において既に計算され記憶装置５に記憶されている粒子の位置ベクトルを重心とする多面体要素内のｐ点、ｑ点の位置ベクトルを読み込む。
【０１３３】
次に、制御部３は、ステップ１０５で既に計算され記憶装置５に記憶されているコイルが磁性体を構成する粒子を重心とする多面体要素内部のｐ点に作る磁場ベクトルを外部からの印加磁場ベクトルとして読み込む。
【０１３４】
さらに、制御部３は、ステップ１０４で既に計算され、記憶装置５に記憶されている、磁性体を構成する粒子の要素境界面積、法線ベクトルおよび粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
【０１３５】
また、制御部３はあらかじめ記憶装置５に記憶されている減衰定数、仮想質量、仮想時間刻みを読み込む。
【０１３６】
次に、制御部３はあらかじめ記憶装置５に記憶されている磁場ベクトル、磁場ベクトルの時間微分の初期値を読み込む。
【０１３７】
なお、制御部３は後述するステップ１０７において磁場ベクトル、磁場ベクトルの時間微分が更新されている場合は、その値を読み込む。
【０１３８】
なお、磁化ベクトルは、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に磁場ベクトルを代入することで計算される。
【０１３９】
制御部３は、式（１７）、式（１８）、式（１９）を、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内の全てのｐ点に対して計算し、計算された磁場ベクトルを記憶装置５に記憶する。
【０１４０】
以上が、図５のステップ１０６の詳細である。
【０１４１】
次に、磁場解析装置１の制御部３はステップ１０６で計算し記憶装置５に記憶されている、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のｐ点の磁場ベクトルに束縛力を加え、計算された磁場ベクトルを、記憶装置５に記憶する（図５のステップ１０７）。
【０１４２】
ここで、図５のステップ１０７について、具体的に説明する。
【０１４３】
磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のｐ点の磁場ベクトルに以下に示す式（２０）、式（２１）に従って束縛力を加える。
【０１４４】
【数８】
【０１４５】
ここで、太字のＨは磁場ベクトル、αは仮想質量、δｔは仮想時間刻み、γは減衰定数、Ｎは磁性体を構成する粒子数、ｓは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量、ｉｑは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｑ点の物理量、添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１４６】
磁場解析装置１の制御部３は、演算プログラム１９を用いて、ステップ１０６で計算され記憶装置５に記憶されている磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルを読み込む。
【０１４７】
制御部３は、あらかじめ記憶装置５に記憶されている仮想質量、仮想時間刻み、減衰定数を読み込む。
【０１４８】
制御部３は、磁性体を構成する粒子各々に対して式（２０）、式（２１）に基づく計算を行い、計算された磁場ベクトルを記憶装置５に記憶する。
【０１４９】
以上が、図５のステップ１０７の詳細である。
【０１５０】
次に、磁場解析装置１の制御部３はステップ１０７で求めた磁場ベクトルが束縛条件を満たしているかを判断し、満たしていれば次のステップに進み、満たしていなければステップ１０７に戻る（図５のステップ１０８）。
【０１５１】
具体的には制御部３は、ステップ１０７で計算し記憶装置５に記憶されている磁性体を構成する粒子各々の磁場ベクトルから計算される磁化ベクトルを用いて式（２２）に基づく計算を行う。
【０１５２】
【数９】
【０１５３】
ここで、添え字のｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１５４】
ｅｒｒｉは磁性体を構成する粒子の内、ｉ番目の粒子の束縛条件に対する誤差値である。
【０１５５】
磁化ベクトル（太字の）Ｍはステップ１０７において計算され記憶装置５に記憶されている磁場ベクトルを制御部３が読み込み、ステップ１０１より入力情報２１に記憶されている磁化曲線を表す関数に代入することで計算される。
【０１５６】
法線ベクトル（太字の）ｎには、ステップ１０４において既に計算され記憶装置５に記憶されているものを制御部３は読み込み代入する。
【０１５７】
制御部３は、すべての粒子に対して、誤差の値が式（２３）を満たさなければステップ１０７に戻る。
【０１５８】
【数１０】
【０１５９】
式（２３）においてＡ１は任意の誤差判別値であり、あらかじめ記憶装置５に任意の値が記憶されている。
【０１６０】
以上が図５のステップ１０８の詳細である。
【０１６１】
次に、磁場解析装置１の制御部３は、磁性体の磁化現象が定常状態に到達したかを判断し、条件を満たしていれば次のステップに進む。（図５のステップ１０９）
ここで、図５のステップ１０９について、具体的に説明する。
【０１６２】
磁性体を構成する粒子が定常状態に到達したかは、以下の式（２４）により判断される。
【０１６３】
【数１１】
【０１６４】
ここで、太字のＨは磁場ベクトル、μ0は真空の透磁率、Ｎは磁性体を構成する粒子の粒子数、ｓは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点での物理量である。
【０１６５】
さらに、添え字のｎは任意の整数であり、ｎδｔにおける物理量、ｎ―１は（ｎ―１）δｔ、ｎ＋１は（ｎ＋１）δｔおける物理量に対応している。
【０１６６】
また、Ａ２は磁性体の磁場ベクトルが定常状態に到達したかを判断するための任意の誤差判定値である。
【０１６７】
磁場解析装置１の制御部３は、あらかじめ記憶装置５に記憶されているＡ２、μ0を読み込む。
【０１６８】
また、制御部３は、ステップ１０７で計算され記憶装置５に記憶されている、磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルを読み込む。
【０１６９】
次に、制御部３は式（２４）を計算し、式（２４）を満たしていれば次のステップに進み、満たしていなければステップ１０６に戻る。
【０１７０】
以上が図５のステップ１０９の詳細である。
【０１７１】
次に、磁場解析装置１の制御部３は、磁性体を構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁化ベクトル、磁束密度ベクトルと、コイルを構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁束密度ベクトルを計算し、記憶装置５に記憶する（図５のステップ１１０）。
【０１７２】
以下にステップ１１０を具体的に説明する。
【０１７３】
（磁性体を構成する粒子）
磁性体を構成する粒子の位置ベクトル上の磁場ベクトルは以下の式（２５）で表される。
【０１７４】
【数１２】
【０１７５】
ここで、（太字の）ｅx＝（１,０,０）、（太字の）ｅｙ＝（０,１,０）、（太字の）ｅｚ＝（０,０,１）である。
【０１７６】
太字のＨは磁場ベクトル、太字のｎは法線ベクトルであり、添え字ｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を表し、添え字ｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量である。
【０１７７】
磁性体を構成する粒子の磁化ベクトルは、磁場ベクトルを磁化曲線を表す関数に代入することで求まる。
【０１７８】
磁性体を構成する粒子の位置ベクトル上での磁束密度ベクトルは、式（２６）で表される。
【０１７９】
【数１３】
【０１８０】
ここで、太字のＢは磁束密度ベクトル、太字のＨは磁場ベクトル、太字のＭは磁化ベクトルであり、添え字のｉは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を指し示す。
【０１８１】
μ0は真空の透磁率である。
【０１８２】
磁場解析装置１の制御部３は、ステップ１１０までで計算されている、磁性体を構成する粒子を重心とする多面体要素内のｐ点での磁場ベクトルを記憶装置５から読み込む。
【０１８３】
次に、制御部３は、ステップ１０４で既に計算され記憶装置５に記憶されている法線ベクトルを読み込み、式（２５）に基づいて磁性体を構成する粒子の位置ベクトル上での磁場ベクトルを計算し記憶装置５に記憶する。
【０１８４】
さらに、制御部３は、このステップで計算された磁場ベクトルを、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に代入し、磁性体を構成する粒子の磁化ベクトルを計算し記憶装置５に記憶する。
【０１８５】
次に制御部３は、このステップで計算された磁場ベクトルと磁化ベクトルを式（２６）に代入し、磁性体を構成する粒子の位置ベクトル上の磁束密度ベクトルを計算し記憶装置５に記憶する。
【０１８６】
なお、真空の透磁率はあらかじめ記憶装置５に記憶されているものを制御部３が読み込む。
【０１８７】
（コイルを構成する粒子）
コイルを構成する粒子の位置ベクトル上の磁場ベクトルは式（２７）、式（２８）で記載される。
【０１８８】
【数１４】
【０１８９】
ここで、太字のｒは位置ベクトル、太字のＨは磁場ベクトル、Ｎは磁性体を構成する粒子の数、ｓおよびΔＳは磁性体を構成する粒子を重心とする多面体要素の境界面数と境界面積、太字のＭは磁化ベクトル、太字のｎは法線ベクトルであり、添え字のｉはコイルを構成する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｊは磁性体を構成する粒子の内、ｊ番目の粒子の位置ベクトル上での物理量であり、ｊｑは磁性体を構成する粒子の内、ｊ番目の粒子の位置ベクトルを重心とする多面体要素のｑ点での物理量を指し示す。
【０１９０】
式（２７）の右辺第一項は、磁性体を構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルを記述する項であり、式（２８）で表される。
【０１９１】
式（２７）の右辺第二項は、コイルを構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルであり式（１）〜（１２）で表される。
【０１９２】
磁場解析装置１の制御部３は、ステップ１０３において計算し記憶装置５に記憶されているコイルを構成する粒子の位置ベクトルを記憶装置５より読み込む。
【０１９３】
なお、ステップ１１３においてコイルを構成する粒子の位置ベクトルが更新されている場合は、制御部３はその値を読み込む。
【０１９４】
次に、制御部３はステップ１０２より計算され記憶装置５に記憶されている磁性体を構成する粒子の数と多面体要素の面数、および粒子の位置ベクトルを重心とする多面体要素内のｑ点の位置ベクトルと多面体要素境界面積を読み込む。
【０１９５】
また、制御部３はこのステップで計算され記憶装置５に記憶されている、磁性体を構成する粒子の磁化ベクトルを読み込む。
【０１９６】
さらに、制御部３はステップ１０３において計算され記憶装置５に記憶されているコイルの寸法およびステップ１０１で入力情報２１として記憶装置５に記憶されている電流密度ベクトルを読み込む。
【０１９７】
なお、ステップ１１６において電流密度ベクトルが更新されている場合は、制御部３はその値を読み込む。
【０１９８】
その後、制御部３は式（２７）、式（２８）および式（１）〜（１２）に従い、コイルを構成する粒子の位置ベクトル上での磁場ベクトルを計算し、記憶装置５に記憶する。
【０１９９】
制御部３は、コイルを構成する粒子の位置ベクトル上での磁束密度ベクトル（太字の）Ｂｉを、すでに計算されているコイルを構成する粒子の位置ベクトル上の磁場ベクトルに真空の透磁率を掛けて計算し記憶装置５に記憶する。
【０２００】
真空の透磁率は、あらかじめ記憶装置５に記憶されている。
【０２０１】
以上が、ステップ１１０の詳細である。
【０２０２】
このように、制御部３は、磁場の運動方程式を解くことにより磁性体の磁化現象を解析する。
【０２０３】
そのため、有限要素法のように空間全域をメッシュ分割する必要がなく、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる。
【０２０４】
また、運動方程式を解くため、行列を扱わず、計算に必要なメモリ量は粒子数に比例する。
【０２０５】
次に、磁場解析装置１の制御部３はＳＰＭモータ３１を構成する全ての粒子に働く力ベクトルを計算し、記憶装置５に記憶する（図５のステップ１１１）。
【０２０６】
以下に図５のステップ１１１を具体的に説明する。
【０２０７】
（磁性体を構成する粒子に働く力ベクトル）
正準変数を太字のｒｉ、太字の傍点付きｒｉとし、式（１３）のラグランジアンをラグラ
ンジュの運動方程式に代入する。
【０２０８】
すると、磁性体を構成する粒子に働く力は式（２９）のように記載される。
【０２０９】
【数１５】
【０２１０】
ここで、μ0は真空の透磁率、太字のＨは磁場ベクトル、太字のＭは磁化ベクトル、太字のｒは位置ベクトル、ΔＶは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、Ｎは磁性体を構成する粒子数であり、添え字のｉ、ｊは磁性体を構成する粒子の内、ｉ番目、ｊ番目の粒子の物理量を指し示し、添え字のｉｐは磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする多面体要素内のｐ点の物理量を指し示す。
【０２１１】
また、φ（ｒｉ−ｒｊ）はｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギである。
【０２１２】
本実施形態ではＳＰＭモータを剛体（粒子間の相対距離が不変）とする。
【０２１３】
剛体モデルとすると、磁性体を構成する粒子に働く力は式（３０）のように記載される。
【０２１４】
【数１６】
【０２１５】
磁場解析装置１の制御部３は、ステップ１１０までにおいて計算され記憶装置５に記憶されている磁性体を構成する粒子を重心とする多面体要素内のｐ点での磁場ベクトルを読み込む。
【０２１６】
制御部３は、ｐ点での磁場ベクトルを、ステップ１０１より入力情報２１として記憶装置５に記憶されている磁化曲線を表す関数に代入することによりｐ点での磁化ベクトルを計算し記憶装置５に記憶する。
【０２１７】
また、制御部３は、ステップ１０４で計算され記憶装置５に記憶されている磁性体を構成する粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
【０２１８】
次に、制御部３は、磁性体を構成する粒子を重心とする多面体要素内のｐ点の磁場ベクトルの偏微分を計算し、その値に基づいて、式（３０）より力ベクトルを計算し記憶装置５に記憶する。
【０２１９】
（コイルを構成する粒子に働く力ベクトル）
コイルを構成する粒子のうち、ｉ番目の粒子の位置ベクトルを太字のｒｉ、その位置ベクトルでの電流の単位ベクトルを太字のｔｉ、磁束密度ベクトルを太字のＢｉ、電流の流れる方向のコイルの長さをＬｉとすると、コイルを構成する粒子の内、ｉ番目の粒子に働く力ベクトルは以下の式（３１）で記載される。
【０２２０】
【数１７】
【０２２１】
磁場解析装置１の制御部３は、ステップ１０３において計算され記憶装置５に記憶されているコイルを構成する粒子の位置ベクトルを読み込む。
【０２２２】
なお、ステップ１１３においてコイルを構成する粒子の位置ベクトルが更新されている場合は制御部３はその値を読み込む。
【０２２３】
また、制御部３はステップ１０３において計算され記憶装置５に記憶されているコイルの寸法を読み込む。
【０２２４】
次に、制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されているコイルの電流密度ベクトルを読み込む。
【０２２５】
なお、制御部３は、ステップ１１６においてコイルの電流密度ベクトルが更新されている場合はその値を用いる。
【０２２６】
制御部３は、式（３１）に基づいて、コイルを構成する粒子に働く力ベクトルを計算し記憶装置５に記憶する。
【０２２７】
以上が、図５のステップ１１１の詳細である。
【０２２８】
次に、磁場解析装置１の制御部３は、ステップ１１１で計算した可動部（移動子）を構成する粒子に働く力ベクトルから、可動部（移動子）の回転量、角速度、並進量、並進速度を計算し、記憶装置５に記憶する（図５のステップ１１２）。
【０２２９】
以下に図５のステップ１１２を具体的に説明する。
【０２３０】
剛体の並進運動の運動方程式は式（３２）で表される。
【０２３１】
【数１８】
【０２３２】
ここで、ｍｔｏｔは可動部（移動子）の質量、太字のｒｃは可動部（移動子）において任意に定めた代表点の位置ベクトル、太字のＦは可動部（移動子）を構成する粒子に働く力ベクトル、Ｎは可動部（移動子）を構成する粒子数である。
【０２３３】
次に、任意の固定点の周りの回転運動の運動方程式は式（３３−１）、式（３３−２）、式（３３−３）で表される。
【０２３４】
【数１９】
【０２３５】
ここで、Ｉ１、Ｉ２、Ｉ３は剛体の主慣性モーメント、ω１、ω２、ω３は慣性主軸座標系での角速度、Ｎ１、Ｎ２、Ｎ３は慣性主軸座標系でのトルクである。
【０２３６】
本実施例のＳＰＭモータ３１のロータ３３はロータシャフト４１を通じてある回転軸に固定されている。
【０２３７】
固定された回転軸をｚ軸とすると、ロータ３３の運動は以下のように計算される。
【０２３８】
ロータ３３はロータシャフト４１を通じてｚ軸に固定されているため、ｚ軸周りの回転運動の運動方程式を解くことで、ロータ３３の運動は計算される。
【０２３９】
ｚ軸周りのロータ３３の回転運動の運動方程式は式（３４）で表される。
【０２４０】
【数２０】
【０２４１】
ここで、Ｉｚはロータ３３の慣性モーメント、ωｚはｚ方向の角速度、Ｎｚはｚ方向のトルクである。
【０２４２】
式（３４）を蛙飛び法で離散化すると、以下の式（３５）、式（３６）を得る。
【０２４３】
【数２１】
【０２４４】
ここで、θｚはｚ軸周りの回転量であり、ｄｔは微小時間刻み幅である。
【０２４５】
また、添え字のｎ−１／２、ｎ、ｎ＋１／２、ｎ＋１はそれぞれ（ｎ―１／２）ｄｔ、ｎｄｔ、（ｎ＋１／２）ｄｔ、（ｎ＋１）ｄｔ秒での物理量に対応している。
【０２４６】
時間刻み幅は記憶装置５にあらかじめ任意の値が記憶されているものを制御部３が読み込む。
【０２４７】
回転量、角速度の初期値は記憶装置５にあらかじめ任意の値が記憶装置５に記憶されている。
【０２４８】
なお、すでにステップ１１２においてロータ３３の回転量、角速度が更新されている場合はその値を用いる。
【０２４９】
制御部３は、ステップ１０１より入力情報２１として記憶装置５に記憶されているロータ３３の三次元構造と質量密度により慣性モーメントを計算する。
【０２５０】
さらに、制御部３は、ロータ３３を構成する粒子の位置ベクトルと、ステップ１１１で計算し記憶装置５に記憶されているロータ３３を構成する粒子に働く力ベクトルに基づいて、式（３５）、式（３６）を計算し、ｄｔ秒後のロータ３３の回転量と角速度を計算し記憶装置５に記憶する。
【０２５１】
以上が、図５のステップ１１２の詳細である。
【０２５２】
次に、磁場解析装置１の制御部３はステップ１１２で計算した微小時間後の回転量や重心の移動量を基に、移動子を構成している粒子の位置ベクトルを計算し、記憶する（図５のステップ１１３）。
【０２５３】
次に、磁場解析装置１の制御部３は必要に応じてＳＰＭモータ３１を構成する粒子の位置ベクトル、力ベクトル、磁場ベクトル、磁束密度ベクトル、磁化ベクトルをプリンタポート１１を介してプリンタ１２より出力する（図５のステップ１１４）。
【０２５４】
次に、磁場解析装置１の制御部３は所定の終了条件（時間、移動量等）を満たしているかを判断し、満たしている場合は解析を終了し、満たしていない場合はステップ１１６に進む（図５のステップ１１５）。
【０２５５】
次に、磁場解析装置１の制御部３は必要に応じてコイル４５に流れる電流密度ベクトルを更新し、ステップ１０４に戻る（図５のステップ１１６）。
【０２５６】
以上が本実施形態に係る磁場解析の手順である。
【０２５７】
このように、本実施形態によれば、磁場解析装置１がＳＰＭモータ３１の磁場解析を元に、移動子であるロータ３３の移動量を求めている。
【０２５８】
そのため、本実施形態では磁場と運動の相互作用を演算でき、磁場解析と機構の連成が可能となる。
【０２５９】
また、本実施形態では、コイル４５に流れる電流密度ベクトルを更新しながら解析を行うことができる。
【０２６０】
そのため、コイルに通電する電流を変化させながら解析を行うことができ、従来よりもさらに精度良く磁場解析を行うことができる。
【０２６１】
上記した実施形態では、本発明を磁性体を構成する粒子のラグランジアンから導出される磁場の運動方程式を解くことで磁性体の磁化現象を計算した場合について説明したが、本発明は、何等、これに限定されることなく、運動方程式であればラグランジアンから導出されたもの以外のものであってもよい。
【０２６２】
また、上記した実施形態では本発明をＳＰＭモータ３１の磁場解析に使用したが、本発明は、何等、これに限定されることなく、磁性体が配置された空間における任意の点の磁場を解析するものであれば、例えばリニアモータ等のＳＰＭモータ３１以外のモータ、あるいは磁性体で空間を囲む磁気シールドにおいて、磁性体でシールドされた空間内の磁場解析に用いても良い。
【０２６３】
上述のように本発明の一実施形態によれば、物体に作用する磁場についての磁場解析と、その磁場により生じる物体の運動についての機構解析との連成解析をするための連成解析装置が提供される。すなわち、この連成解析装置は、磁場解析演算部と、機構解析演算部と、を備える。一実施例においては、制御部３は、磁場解析演算部と、機構解析演算部と、を備えてもよい。磁場解析演算部は、物体を構成する各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式に基づき磁場を演算する。機構解析演算部は、演算された磁場に基づき物体の運動を演算し、物体の各粒子の位置を更新する。磁場解析演算部はさらに、更新された各粒子の位置に基づいて磁場を再計算する。このようにして磁場解析と機構解析とを粒子ごとに繰り返すことにより、連成解析が可能である。なお以下では便宜上、磁場解析演算部を磁場演算部とも称し、機構解析演算部を位置演算部とも称する。
【０２６４】
ここで磁場演算部は、上述のように、物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより物体に作用する磁場を演算する。解析対象物体は磁化曲線が非線形である材料を含んでもよく、例えば磁性体を解析対象物体としてもよい。本願発明者は、いわゆる解析力学の手法を用いて、物体の各粒子に作用する磁場の支配方程式を運動方程式の形式で記述することに成功した。上述の一実施例においては、上記の式（１３）に定義されるラグランジアンをラグランジュの運動方程式に代入することにより、式（１６）のように各粒子に作用する磁場を運動方程式の形式で記述している。
【０２６５】
磁場演算部は、磁場の運動方程式を数値的に解くことにより、物体に作用する磁場を求める。具体的には例えば、磁場演算部はまず、離散化された磁場の運動方程式に磁場及びその時間微分の初期値を代入し、微小時間経過後の磁場及びその時間微分を求める。ここで、微小時間とは磁場の収束計算を行うために用いられる仮想的な時間であってもよい。つまりこの仮想時間は実際の経過時間とは異なる。磁場演算部は、演算された微小時間経過後の磁場が収束判定基準を満たすか否かを判定する。収束判定条件は、一実施例においては式（２４）で与えられる。磁場演算部は、収束判定基準が満たされるまで、離散化された磁場の運動方程式を繰り返し演算することにより磁場の数値解を求める。
【０２６６】
この方法によれば、有限要素法に基づく周知の解析法とは異なり逆行列の演算が不要であるので、少ない計算資源で磁場解析を実行することができるという点で好ましい。また、この方法は、高精度の解析結果を、有限要素法のように解析対象物体とその周囲の空間とを含む系全体をメッシュ分割することなく得られるという点でも好ましい。
【０２６７】
一実施例においては、上述の連成解析がモータの磁場解析に利用される。本願発明者の検証によると、モータの磁場解析に要求される高精度のレベルに収束判定条件を設定したときに、演算終了までの所要時間が比較的長くなる場合があることがわかった。この収束の遅れは、解析対象物体における磁気飽和に起因すると考えられる。すなわち、磁化曲線の非線形性が磁場の数値解の収束性に影響している。
【０２６８】
図１０は、磁性体（例えば鉄）の磁化曲線の一例を模式的に示す図である。図１０の横軸は磁場Ｈを示し、縦軸は磁化Ｍを示す。図１０に示されるように、磁場Ｈが小さいときには磁化Ｍは磁場Ｈに比例する。磁場Ｈが大きくなるにつれて、磁化Ｍと磁場Ｈとの比例定数（いわゆる磁気感受率）は徐々に小さくなっていく。磁場Ｈが小さい場合の比例定数は（例えば図中の点Ａ）、磁場Ｈが大きい場合（図中の点Ｂ）の比例定数の１０００倍程度の大きさとなり得る。このように、磁性体の磁化曲線は非線形である。
【０２６９】
磁場の運動方程式には、磁場及びその時間微分項だけではなく、磁化に依存する磁化依存項が含まれる。ある粒子についての磁場の運動方程式における磁化依存項は例えば、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項である。このような相互作用を表す項は、例えば式（１４）に示されるように、磁化Ｍを陽に含む。したがって、式（１６）に示される磁場の運動方程式は、収束計算を進めるにつれて非線形に変化する項を含むと言える。収束計算の磁場の初期値は通常小さい値（例えばゼロ）にとるから、式１４のＨ（ｒｉｐ）は、収束計算の当初に比べて完了時には１０００分の１程度の小さい値となり得る。方程式におけるこの大きな非線形性によって数値解の収束に遅れが生じると考えられる。
【０２７０】
磁場の運動方程式には、任意に設定することが許容される調整可能パラメタを導入することができる。例えば式（１６）に示される磁場の運動方程式においては、仮想質量α及び減衰定数γは調整可能パラメタである。パラメタの設定作業を簡単にすることを重視する場合には、これらの調整可能パラメタを一定値に設定してもよい。この場合、収束計算中に調整可能パラメタは不変である。一方、磁場演算部は、収束計算中に調整可能パラメタを変化させてもよい。
【０２７１】
磁場の運動方程式は一般に、束縛を考慮に入れなければ、２階の線形常微分方程式の形式をとる。この形式は、１自由度の振動の基礎方程式と共通である。よって、磁場の運動方程式の収束特性は、振動の基礎方程式の減衰特性に対応する。このことから、振動の方程式の各項の係数に応じて減衰特性が変化することに準じて、磁場の運動方程式の調整可能パラメタを変更することにより所望の収束特性に調整することが可能である。
【０２７２】
そこで、本発明の一実施形態においては、磁場演算部は、磁場の運動方程式に導入されている調整可能パラメタを、磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整してもよい。磁場演算部は、磁化依存項による数値解の収束性への影響を抑えるように調整可能パラメタを変化させてもよい。一実施例においては、磁場演算部は、収束判定基準を満たすか否かを判定する収束判定処理において収束判定基準が満たされていないと判定したときに、調整可能パラメタを変化させる係数調整処理を実行してもよい。上述の実施例においては、制御部３は、誤差の値が誤差判別値Ａ２を超えると判定した場合に（図５のステップ１０９のＮｏ）、仮想質量α及び減衰定数γを調整する係数調整ステップを実行してもよい。
【０２７３】
ある物体の粒子に作用する磁場Ｈの運動方程式は一般に、式（３７）で表される。ここで、係数ＣＭは磁場Ｈの２階時間微分項の係数を表し、係数ＣＤは磁場Ｈの１階時間微分項の係数を表し、係数ＣＫは磁場Ｈの項の係数を表す。係数ＣＭ、ＣＤ、ＣＫはいずれも解析者が調整することのできるパラメタである。一方、係数ＦＮＬは、磁場Ｈと磁化Ｍとの比例定数に依存する。Ｈ０は外部からの印加磁場を表す。
【０２７４】
【数２２】
【０２７５】
よって、調整可能なパラメタである係数ＣＭ、ＣＤを磁化曲線に基づいて調整することにより、磁化曲線の非線形性が数値解の収束性に与える影響を排除することができる。すなわち、磁場Ｈの時間微分項の係数ＣＭ及びＣＤを、式（３８）に示されるように磁化依存係数ＦＮＬを用いてＣＭ’ＦＮＬ及びＣＤ’ＦＮＬと定義する。磁化依存係数ＦＮＬの時間変化率が十分に小さければ、磁化曲線の非線形性が数値解の収束性に与える影響を事実上排除することができる。
【０２７６】
【数２３】
【０２７７】
式（３８）の磁場の項及びその時間微分項には、磁場Ｈに対し非線形に変化する共通パラメタＦＮＬが導入されている。つまり、共通の磁化依存係数ＦＮＬをもたない運動方程式（３９）の収束特性に、式（３８）に示される磁場の運動方程式の収束特性は実質的に等しくなる。
【０２７８】
式（３９）の磁場Ｈ及びその時間微分項の係数ＣＭ’、ＣＤ’、及びＣＫはいずれも磁化Ｍに依存せず、解析者が任意に設定することのできるパラメタである。よって、式（３８）において磁化依存係数ＦＮＬが毎回の収束計算ごとに変化しても、係数ＣＭ’、ＣＤ’、及びＣＫから定まる一定の収束特性を得ることができる。なおこの場合、所望の収束特性をもたらすよう係数ＣＭ’、ＣＤ’、及びＣＫを予め設定することが好ましい。
【０２７９】
なお本発明は、磁化曲線の非線形性の影響を磁場の運動方程式の収束特性から排除する実施例には限られない。基礎方程式の変数の係数がその変数に対して非線形に変化する場合における運動方程式型の数値解法にも適用可能である。例えばいわゆる大変形解析にも適用可能であろう。
【０２８０】
一実施例においては、第１の物理量と、該第１の物理量と非線形の関係をもつ第２の物理量とを有する物体の特性を解析するための解析装置が提供される。例えば上述の実施例では第１の物理量及び第２の物理量はそれぞれ磁場及び磁化に対応する。この解析装置は、運動方程式の形式で記述された第１の物理量が満たすべき関係式を数値的に解く数値演算部を備える。この関係式は、第２の物理量に依存する項と、第２の物理量に依存しない項と、を含み、第２の物理量に依存しない項の係数が第２の物理量に依存する項に連動するよう定義されている。一実施例においては、第１の物理量及びその時間微分項のそれぞれの係数には第２の物理量に連動する連動部分が共通に導入されていてもよい。
【０２８１】
また、一実施例においては、調整可能パラメタは、磁場の運動方程式を数値的に解く過程で導出される値を用いて定義されていることが好ましい。つまり、磁場演算部は、磁場演算処理において算出された値を用いて係数調整処理を実行してもよい。このようにすれば、磁場の運動方程式を解く過程で得られた計算済みの値を流用して調整可能パラメタを調整することができる。調整可能パラメタの設定のために新たな演算をする必要がなくなるので、計算コストを低減することができる。よって、収束演算の所要時間を短くすることもできる。
【０２８２】
そこで、一実施例においては、調整可能パラメタは、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項に連動して調整されてもよい。式（１６）を解くことにより磁場を求める上述の実施例においては例えば、仮想質量αｉｐを式（４０）及び式（４１）により可変量として定義する。係数ｋｉｐは磁場の運動方程式に存在するＨ（ｒｉｐ）の項を用いて定義されているから、調整可能パラメタの設定のために新たな演算をする必要がない。なお係数α０は解析者が任意に設定することのできるパラメタである。
【０２８３】
【数２４】
【０２８４】
この仮想質量αｉｐにより、磁化曲線の非線形性の影響が磁場の運動方程式の収束特性から排除されることを説明する。上述のようにＳＨＡＫＥ法を適用して式（１６）を解くことを考慮すると、式（１６）の右辺第２項は０に収束するから、式（４２）のように表される。
【０２８５】
【数２５】
【０２８６】
式（１４）及び式（４１）から、Ｈ（ｒｉｐ）は式（４３）で表される。
【０２８７】
【数２６】
【０２８８】
式（４０）及び式（４３）を式（４２）に代入することにより、磁場の運動方程式は式（４４）で表される。式（４４）の磁場の項及びその時間微分項には、磁場Ｈに対し非線形に変化する共通パラメタ（１−ｋｉｐ）が導入されている。よって、磁場の運動方程式（４４）の収束特性は、共通パラメタ（１−ｋｉｐ）を有しない式（４５）の収束特性に実質的に等しくなる。よって、磁気飽和に起因して係数ｋｉｐ（式（１６）で言えばＨ（ｒｉｐ））が収束計算中に変化しても、磁場の運動方程式（４４）の収束特性は変化しない。係数α０、γを予め設定することにより、所望の収束特性を実現することができる。
【０２８９】
【数２７】
【０２９０】
なお、上述のパラメタの更新は、磁場解析の対象物体が磁気感受率の異なる物質を含む場合にも有効であると考えられる。例えば物体が磁石と鉄とを含む場合にも有効である。物質が異なることによる粒子ごとの磁気感受率の違いは、磁気飽和が生じているか否かによる磁気感受率の違いと同様に上記の磁場の運動方程式において取り扱われているからである。
【０２９１】
図１１は、本発明の一実施形態に係る磁場演算処理を説明するためのフローチャートである。図１１に示されるように、磁場解析装置１の制御部３はまず、磁場の運動方程式を作成する（Ｓ１０）。一実施例においては、制御部３は、図５に示されるステップ１０１乃至１０５を実行することにより磁場の運動方程式（１６）を作成する。ここで、仮想質量αｉｐは式（４０）及び式（４１）により可変量として定義されている。
【０２９２】
次に制御部３は、磁場の運動方程式に基づいて磁場を演算する（Ｓ１１）。一実施例においては、制御部３は、図５に示されるステップ１０６乃至１０８を実行することにより、磁場の運動方程式（１６）に基づいて磁場を演算する。制御部３は、演算された磁場が収束判定基準を満たすか否かを判定する（Ｓ１２）。収束判定基準は例えば図５に示すステップ１０９と同様であり、制御部３は、誤差が誤差判別値を満たすか否かを判定する。
【０２９３】
制御部３は、収束判定基準が満たされていないと判定した場合には（Ｓ１２のＮｏ）、調整可能パラメタを更新して（Ｓ１３）、磁場の運動方程式に基づく磁場演算を再度行う（Ｓ１１）。制御部３は、収束判定基準が満たされるまで、調整可能パラメタの更新と磁場演算とを繰り返す。制御部３は、調整可能パラメタとして仮想質量αｉｐを更新する。すなわち、制御部３は係数ｋｉｐを更新する。
【０２９４】
制御部３は、収束判定基準が満たされたと判定した場合には（Ｓ１２のＹｅｓ）、演算された磁場に基づいて解析対象物体の運動を演算する（Ｓ１４）。一実施例においては、制御部３は、図５に示されるステップ１１０乃至１１４を実行する。次に制御部３は、所定の終了条件が満たされているか否かを判定し（Ｓ１５）、満たされていると判定した場合には解析を終了する（Ｓ１５のＹｅｓ）。制御部３は、終了条件が満たされていないと判定した場合には（Ｓ１５のＮｏ）、運動による解析対象物体の新たな位置に基づいて再度磁場を演算する。このようにして、制御部３は、物体に作用する磁場とそれによる物体の運動との連成解析を実行する。
【０２９５】
図１２及び図１３は、本発明の一実施形態による磁場演算の例を示す図である。図１２は仮想質量αｉｐを一定値に固定した場合の磁場Ｈと仮想時間δｔとの関係を示し、図１３は上述のように係数ｋｉｐを用いて仮想質量αｉｐを更新しながら収束計算をしたときの磁場Ｈと仮想時間δｔとの関係を示す。図１２及び図１３の磁場演算の例は、１次元に配列された１０個の粒子からなる物体に対して行ったものである。粒子１及び粒子１０が物体の両端である。図１２の演算結果に比べて図１３の演算結果は磁場が短時間で定常状態に達していることがわかる。仮想質量αｉｐを更新しながら収束計算を行うことにより、磁場の運動方程式の数値解の磁気飽和に伴う収束の遅れを低減することができる。
【符号の説明】
【０２９６】
１ 磁場解析装置、 ３ 制御部、 ５ 記憶装置、 ７ 入力部、 ９ 表示部、 １１ プリンタポート、 １２ プリンタ、 １３ バス、 ３１ ＳＰＭモータ、 ３３ ロータ、 ３５ ステータ、 ３７ ロータコア、 ３９ 永久磁石、 ４１ ロータシャフト、 ４３ ステータティース、 ４４ コアバック、 ４５ コイル、 ４６ フレーム。
【特許請求の範囲】
【請求項１】
物体を含む系における磁場を解析するための磁場解析装置であって、
物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより前記物体に作用する磁場を演算する磁場演算部を備え、
前記磁場の運動方程式には調整可能なパラメタが導入されており、
前記磁場演算部は、前記物体の磁化曲線の非線形性が数値解の収束性に与える影響を緩和するよう前記パラメタを調整することを特徴とする磁場解析装置。
【請求項２】
前記磁場の運動方程式は、磁場の項及びその時間微分項を含み、
前記パラメタは、前記時間微分項の係数を含み、当該係数は、前記磁化曲線に基づいて調整されることを特徴とする請求項１に記載の磁場解析装置。
【請求項３】
前記磁場の運動方程式は、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項を含み、前記パラメタは、当該項に連動して調整されることを特徴とする請求項１または２に記載の磁場解析装置。
【請求項４】
前記パラメタは、前記磁場の運動方程式を数値的に解く過程で導出される値を用いて定義されていることを特徴とする請求項１から３のいずれかに記載の磁場解析装置。
【請求項５】
演算された磁場に基づいて前記物体の各粒子の位置を更新する位置演算部をさらに備え、
前記磁場演算部は、更新された各粒子の位置に基づいて磁場を再計算することを特徴とする請求項１から４のいずれかに記載の磁場解析装置。
【請求項６】
磁化曲線が非線形である材料を含む物体に作用する磁場が満たすべき関係を運動方程式の形式で記述した磁場の運動方程式を数値的に解く磁場演算ステップと、
前記磁場の運動方程式に導入されている調整可能パラメタを、前記磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整する調整ステップと、備えることを特徴とする磁場解析方法。
【請求項７】
請求項６に記載の方法をコンピュータに実行させるためのプログラム。
【請求項８】
請求項７に記載のプログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項９】
物体を構成する各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式に基づき磁場を演算する磁場解析演算部と、
演算された磁場に基づき前記物体の運動を演算する機構解析演算部と、を備えることを特徴とする連成解析装置。
【請求項１０】
第１の物理量と、該第１の物理量と非線形の関係をもつ第２の物理量とを有する物体の特性を解析するための解析装置であって、
運動方程式の形式で記述された第１の物理量が満たすべき関係式を数値的に解く数値演算部を備え、
前記関係式は、第２の物理量に依存する項と、第２の物理量に依存しない項と、を含み、第２の物理量に依存しない項の係数が第２の物理量に依存する項に連動するよう定義されていることを特徴とする解析装置。
【請求項１】
物体を含む系における磁場を解析するための磁場解析装置であって、
物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより前記物体に作用する磁場を演算する磁場演算部を備え、
前記磁場の運動方程式には調整可能なパラメタが導入されており、
前記磁場演算部は、前記物体の磁化曲線の非線形性が数値解の収束性に与える影響を緩和するよう前記パラメタを調整することを特徴とする磁場解析装置。
【請求項２】
前記磁場の運動方程式は、磁場の項及びその時間微分項を含み、
前記パラメタは、前記時間微分項の係数を含み、当該係数は、前記磁化曲線に基づいて調整されることを特徴とする請求項１に記載の磁場解析装置。
【請求項３】
前記磁場の運動方程式は、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項を含み、前記パラメタは、当該項に連動して調整されることを特徴とする請求項１または２に記載の磁場解析装置。
【請求項４】
前記パラメタは、前記磁場の運動方程式を数値的に解く過程で導出される値を用いて定義されていることを特徴とする請求項１から３のいずれかに記載の磁場解析装置。
【請求項５】
演算された磁場に基づいて前記物体の各粒子の位置を更新する位置演算部をさらに備え、
前記磁場演算部は、更新された各粒子の位置に基づいて磁場を再計算することを特徴とする請求項１から４のいずれかに記載の磁場解析装置。
【請求項６】
磁化曲線が非線形である材料を含む物体に作用する磁場が満たすべき関係を運動方程式の形式で記述した磁場の運動方程式を数値的に解く磁場演算ステップと、
前記磁場の運動方程式に導入されている調整可能パラメタを、前記磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整する調整ステップと、備えることを特徴とする磁場解析方法。
【請求項７】
請求項６に記載の方法をコンピュータに実行させるためのプログラム。
【請求項８】
請求項７に記載のプログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項９】
物体を構成する各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式に基づき磁場を演算する磁場解析演算部と、
演算された磁場に基づき前記物体の運動を演算する機構解析演算部と、を備えることを特徴とする連成解析装置。
【請求項１０】
第１の物理量と、該第１の物理量と非線形の関係をもつ第２の物理量とを有する物体の特性を解析するための解析装置であって、
運動方程式の形式で記述された第１の物理量が満たすべき関係式を数値的に解く数値演算部を備え、
前記関係式は、第２の物理量に依存する項と、第２の物理量に依存しない項と、を含み、第２の物理量に依存しない項の係数が第２の物理量に依存する項に連動するよう定義されていることを特徴とする解析装置。
【図１】
【図２】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図３】
【図２】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図３】
【公開番号】特開２０１１−１８００９２（Ｐ２０１１−１８００９２Ａ）
【公開日】平成２３年９月１５日（２０１１．９．１５）
【国際特許分類】
【出願番号】特願２０１０−４７０４１（Ｐ２０１０−４７０４１）
【出願日】平成２２年３月３日（２０１０．３．３）
【出願人】（０００００２１０７）住友重機械工業株式会社 (2,241)
【Ｆターム（参考）】
【公開日】平成２３年９月１５日（２０１１．９．１５）
【国際特許分類】
【出願日】平成２２年３月３日（２０１０．３．３）
【出願人】（０００００２１０７）住友重機械工業株式会社 (2,241)
【Ｆターム（参考）】
[ Back to top ]