以下、図面に基づいて本発明に好適な実施形態を詳細に説明する。
まず、図1を参照して、本実施形態に係る磁場解析装置1のハードウェア構成を説明する。
なお、図1のハードウェア構成は例示であり、これに限定されないのは当然である。
図1に示すように、磁場解析装置1は制御部3、記憶装置5、メディア入出力部6、入力部7、表示部9、プリンタポート11等がバス13を介して互いに接続されている。
制御部3は、CPU(Central Processing Unit)、ROM(Read Only Memory)、RAM(Random Access Memory)等で構成され、記憶手段としての記憶装置5に格納されたプログラムに従って、バス13を介して接続された各装置を駆動制御する。
図2に示すように、記憶装置5には、磁場解析装置1の各構成部分を駆動制御するための制御プログラム15、本発明を実施するための磁場解析プログラム17が格納されている。
磁場解析プログラム17は、解析条件を有する情報である入力情報21と、入力情報21と磁場の運動方程式に基づいて磁場を演算し力学系の運動方程式を演算する演算プログラム19とを有している。
メディア入出力部6は、フロッピー(登録商標)ディスク、CD、DVD等のメディアとの間で情報の入出力を行う装置である。
入力部7は、キーボード、マウス等の入力装置であり、表示部9はディスプレイ等の表示機器である。
プリンタポート11には出力装置としてのプリンタ12等が接続される。
次に、磁場解析装置1を用いた磁場の解析の手順について図3〜9を参照して説明する。
ここでは永久磁石モータ(Permanent Magnet Motor)の一種である、SPMモータ31(Surface Permanent Magnet Motor)の磁場解析を例にして説明する。
まず、SPMモータ31の構成の概略を図3および図4を参照して説明する。
図3に示すように、SPMモータ31は回転子(移動子)であるロータ33と固定子であるステータ35を有している。
ロータ33は鉄等の磁性体である円柱状のロータコア37を有し、ロータコア37の表面には永久磁石39が設けられている。
ロータコア37の軸中心には棒状のロータシャフト41が設けられている。
ステータ35は磁性体である歯状のステータティース43とステータティース43の外側に設けられた円筒状の磁性体であるコアバック44、コアバック44の外側に設けられた円筒状のフレーム46から構成されている。
図3および図4に示すように、ステータティース43には、金属等の導電体であるコイル45が巻きつけられている。
なお、実際のSPMモータ31ではコイル45は仕様に応じたターン数でステータティース43に巻きつけられて束となっているが、本実施形態では、図3および図4に描かれているように、コイル一本一本をモデル化せず、コイルの束を一つの導体として扱う。
このような構造のSPMモータ31は、永久磁石39の磁場、およびコイル45に電流を流すことにより発生する磁場によって、ロータ33、ステータ35が磁化する。磁性体の磁気エネルギーの偏差によりSPMモータ31は駆動する。
そのため、SPMモータ31の磁場解析を行うためにはコイル45、永久磁石39がロータ33、ステータ35を構成する磁性体上に作る磁場ベクトルを計算し、これら磁性体の磁化現象を解析する必要がある。
次に、解析の手順について図5〜図9を参照して説明する。
なお、以下の手順においてはSPMモータ31を、複数に要素分割して要素ごとの粒子の集合体とし、剛体モデルとして扱っているが、本発明はこれに限定されることはなく、粒子を用いずに要素ごとの剛体モデルとして扱ってもよい。
また、以下の手順において、磁性体とは、ロータ33、ステータ35を構成する磁性体と永久磁石39を指し示し、磁化曲線を表す関数によってこれらは区別される。
まず、磁場解析装置1の制御部3は磁場解析プログラム17を起動し、解析したいSPMモータ31の解析条件としての三次元構造(形状、座標点)、質量密度、磁性体の磁化曲線を表す関数、導体の電流密度ベクトルを記憶装置5の入力情報21として記憶する(図5のステップ101)。
これらの物理量は例えばメディア入出力部6を介してCD−ROM等の記録媒体から読込んだものであってもよい。
また、SPMモータ31の三次元構造の情報とは例えばCAD等のデータである。
さらに、あらかじめ上記物理量が入力情報21として記憶されている場合は、上記ステップは不要である。
電流密度ベクトルは、コイル45の作る磁場ベクトルを計算する際に必要になる。
さらに、磁性体の磁化曲線を表す関数は磁化ベクトル(太字の)Mを計算する際に必要になる。
以上が図5のステップ101の詳細である。
次に、磁場解析装置1の制御部3は、入力情報21の有する三次元構造の情報から、磁性体をN個の粒子(粒子の位置ベクトルを重心とする多面体要素(本実施例では立方体要素))に分割し、粒子の位置ベクトルを計算し、記憶装置5に記憶する(図5のステップ102)。
ここで、Nは任意の整数であり、粒子の数Nおよび位置ベクトルは入力情報21の有する三次元構造およびあらかじめ記憶装置5に記憶されている多面体要素の形状により計算される。
以上が図5のステップ102の詳細である。
次に、制御部3は、入力情報21の有する三次元構造の情報から、導体(図6(a)に示すコイル45)を図6(b)に示すようにローカル導体(直方体導体45a、45bと円弧状柱状導体45c、45d)に分割し、それぞれの導体が作る磁場ベクトルを計算するための係数を計算し、記憶装置5に記憶する(図5のステップ103)。
ここで、図5のステップ103について、具体的に説明する。
まず、ローカル導体が直方体導体45aの場合について説明する。
なお、ローカル導体が直方体導体45bの場合は、ローカル導体が直方体導体45aの場合と同様であるため、説明を省略する。
ローカル導体が直方体導体45aの場合は、制御部3は、図7に示すように、ローカル座標系(xs,ys,zs)を適用する。
このローカル座標系においてはローカル導体(直方体導体45a)の重心を原点Osとし、直方体導体45aの寸法はxs方向に2a、ys方向に2b、zs方向に2cの長さを持つものとする。
また原点Osに粒子は位置するものとする。
制御部3は、ステップ101より入力情報21として記憶装置5に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体(直方体導体45a)の寸法であるa,b,cと粒子位置ベクトルを計算し記憶装置5に記憶する。
次に、ローカル導体が円弧状柱状導体45cの場合について説明する。
なお、ローカル導体が円弧状柱状導体45dの場合は、ローカル導体が円弧状柱状導体45cの場合と同様であるため、説明を省略する。
ローカル導体が円弧状柱状導体45cの場合は、制御部3は、図8に示すようにローカル座標系(xc,yc,zc)を適用する。
このローカル座標系においては原点Ocは円弧の中心軸上に存在し、かつ円弧状柱状導体45cの高さ方向(図8のzc方向)に対して円弧状柱状導体45cが対称となる点に存在するものとする。
また、xc,yc,zcは、xc−yc平面でみると、+xc軸を基点とし、円弧状柱状導体45cの円弧が+zc軸からみて反時計回りになるようして決定する。
円弧状柱状導体45cの内径と外径の平均値をRcとし、径方向の厚さを2ra、zc方向の高さを2zbとする。
電流は+xc軸を基点として、+zcから見て反時計回りの方向への角度をθとし、電流はこの方向に一様な電流密度jで流れているものとする。
粒子は円筒座標系で(Rc、θ/2、0)に位置するものとする。
制御部3は、ステップ101より入力情報21として記憶装置5に記憶されている導体の三次元構造と電流密度ベクトルを読み込み、ローカル導体(円弧状柱状導体45c)の寸法であるra、zb、θ、Rcおよび粒子位置ベクトルを計算し記憶装置5に記憶する。
さらに、あらかじめ上記物理量が入力情報21として記憶されている場合は、上記ステップは不要である。
以上が図5のステップ103の詳細である。
次に、磁場解析装置1の制御部3は、このステップまでに計算され記憶装置5に記憶されている磁性体を構成するN個の粒子の位置ベクトルと、あらかじめ記憶装置5に記憶されている多面体要素の形状により、磁性体を構成する粒子の磁場の運動方程式の係数を計算し記憶装置5に記憶する(図5のステップ104)。
ここで、図5のステップ104について、具体的に説明する。
分割した多面体要素(本実施形態では立方体要素)を2次元表示すると図9に示す形状となる。
ここで、粒子の位置ベクトルを(太字の)rg、粒子の位置ベクトルを重心とする多面体要素の要素境界面の中点をq点とし、粒子の位置ベクトルとq点との中間点にp点を定義する。
制御部3は、ステップ102において計算し記憶装置5に記憶している磁性体を構成する粒子(粒子を重心とする多面体要素(本実施形態では立方体要素))の位置ベクトルを読み込む。
なお、ステップ113において粒子の位置ベクトルが更新され記憶装置5に記憶されている場合は、制御部3はその値を読み込む。
制御部3は、粒子の位置ベクトルとあらかじめ記憶装置5に記憶されている多面体要素の形状から、磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内のp点とq点の位置ベクトルと要素境界面への法線ベクトル(太字の)nを計算し、記憶装置5に記憶する。
次に、制御部3は、磁性体を構成する粒子を重心とする多面体要素の境界面積ΔS、および粒子の位置ベクトルを重心とする多面体要素の境界面の頂点と粒子位置を頂点とする多面体要素体積ΔVを計算し、記憶装置5に記憶する。
以上が図5のステップ104の詳細である。
次に、磁場解析装置1の制御部3は、ステップ103で記憶装置5に記憶されているコイルの寸法およびこのステップまでに記憶装置5に記憶されている電流密度ベクトルを用いて、ビオ・サバールの法則を積分することにより得られる解析解により、通電されたコイル45が、磁性体を構成する粒子上に作る磁場ベクトルを計算し、記憶装置5に記憶する(図5のステップ105)。
ここで、図5のステップ105について、具体的に説明する。
ここでは、任意の位置ベクトルが図4における磁性体としてのステータティース43のある点Pであると仮定した場合に、通電されたコイル45がP点の位置ベクトル上に作る磁場ベクトルを計算する手順を例に説明する。
まず、ローカル導体が直方体導体45aの場合について説明する。
なお、ローカル導体が直方体導体45bの場合は、ローカル導体が直方体導体45aの場合と同様であるため、説明を省略する。
ローカル導体が直方体導体45aの場合は、図7に示すようにローカル座標系(xs,ys,zs)を適用する。
このローカル座標系はステップ103で説明している。
次に、P点の位置ベクトルをローカル座標系(xs,ys,zs)に変換する。
通電された直方体導体がP点に作る磁場ベクトルは、以下に示す式(1)〜(3)で記載される。
ここで、(太字の)rpsはローカル座標系(xs,ys,zs)でのP点の位置ベクトルであり、xps,yps,zpsはxs,ys,zs方向の値である。
πは円周率である。
Hxs、Hys、Hzsはローカル座標系(xs,ys,zs)における磁場ベクトルの各成分である。
jは電流密度である。
また、xi,yj,zkはxs,ys,zs方向の積分の上限、下限を表しており、式(4)に示す関係が成立する。
ここで、a,b,cは直方体導体の寸法である。
次に、ローカル導体が円弧状柱状導体45cの場合について説明する。
なお、ローカル導体が円弧状柱状導体45dの場合は、ローカル導体が円弧状柱状導体45cの場合と同様であるため、説明を省略する。
ローカル導体が円弧状柱状導体45cの場合は、図8に示すようにローカル座標系(xc,yc,zc)を適用する。
このローカル座標系はステップ103において説明している。
ローカル座標系(xc,yc,zc)に変換後のP点の位置ベクトルを以下の式(5)に示すように円筒座標系(太字の)rpc=(Rpc、φpc、Zpc)に変換する。
通電された円弧状柱状導体45cがP点に作る磁場ベクトルは、以下の式(6)〜(11)で表される。
ここで、Hrc、Htc、Hzcは円筒座標系での磁場ベクトルの各成分である。
jは電流密度である。
ra、θ、zbは円弧状柱状導体45cの寸法であり、Rcは円弧の内径と外径の平均値である。
sgnはZkの符号であり、Rj、Zkは積分の上限、下限を表しており、式(12)に示す関係が成立する。
以上が任意の点Pにコイルが作る磁場ベクトルの計算手順である。
磁場解析装置1の制御部3は、ステップ104で既に計算され記憶装置5に記憶されている磁性体を構成するすべての粒子の位置ベクトルを重心とする多面体要素内部のp点の位置ベクトルを読み込む。
次に、制御部3は、p点の位置ベクトルをローカル座標系(xs,ys,zs)に変換し、ステップ103で計算され記憶装置5に記憶されている直方体導体45aの寸法と、ステップ101より入力情報21として記憶装置5に記憶されている電流密度ベクトルを読み込み、式(1)〜(4)に基づいて磁場ベクトルを計算する。
なお、電流密度ベクトルがステップ116で更新され記憶装置5に記憶されている場合は、制御部3はその値を読み込む。
次に、制御部3は、式(1)〜(4)で計算された磁場ベクトルをグローバル座標系(x,y,z)に変換し、変換後の磁場ベクトルを記憶装置5に記憶する。
これにより、通電された直方体導体45aが、p点に作る磁場ベクトルが求められる。
さらに制御部3は、p点の位置ベクトルをローカル座標系(xc,yc,zc)に変換し、さらに円筒座標系(太字の)rpc=(Rpc、φpc、Zpc)に変換する。
次に制御部3は、ステップ103で計算し記憶装置5に記憶している円弧状柱状導体45cの寸法ra、θ、zbおよびRcを読み込み、また、制御部3はステップ101で入力情報21として記憶装置5に記憶されている電流密度ベクトルを読み込む。
なお、電流密度ベクトルがステップ116で更新され記憶装置5に記憶されている場合は、制御部3はその値を読み込む。
制御部3は式(6)〜(12)に基づいて磁場ベクトルを計算する。
次に、計算された磁場ベクトルを直交座標系に変換し、さらにグローバル座標系(x,y,z)に変換し記憶装置5に記憶する。
これにより、通電された円弧状柱状導体45cがp点に作る磁場ベクトルが求められる。
以上が、図5のステップ105の詳細である。
次に、磁場解析装置1の制御部3は、演算プログラム19を用いて磁性体を構成する粒子の磁場の運動方程式から束縛を考慮せずに仮想時間刻みδt後の磁場ベクトルを計算し、記憶装置5に記憶する(図5のステップ106)。
ここで、図5のステップ106について、具体的に説明する。
磁性体を構成するN個の粒子のラグランジアンを式(13)〜式(15)で表される形とする。
ここで、式(13)において、αは仮想質量、太字のrは位置ベクトル、太字のHは磁場ベクトル、太字の傍点付きHは磁場ベクトルの時間微分、太字のMは磁化ベクトル、太字のnは粒子の位置ベクトルを重心とする多面体要素境界面の法線ベクトル、太字のHextは外部からの印加磁場ベクトル、ΔSは粒子の位置ベクトルを重心とする多面体要素の要素境界面の面積、ΔVは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、χは磁気感受率、μ0は真空の透磁率、λはラグランジュの未定定数、πは円周率、sは粒子の位置ベクトルを重心とする多面体要素の要素境界面数である。
また、各物理量の添え字ipはi番目の粒子の位置ベクトルを重心とする多面体要素内でのp点の物理量、添え字jqはj番目の粒子の位置ベクトルを重心とする多面体要素内でのq点の物理量を示す。
p点、q点はステップ104において説明している。
式(13)において、mは粒子の質量、vは速度、φ(ri−rj)はi番目の粒子とj番目の粒子の相互作用ポテンシャルエネルギである。添え字i、jはそれぞれ、i、j番目の粒子の物理量を示す。
次に、正準変数を(太字の)Hip、(太字の傍点付き)Hipとし、式(13)〜式(15)で示されるラグランジアンをラグランジュの運動方程式に代入すると、磁場の運動方程式は式(16)のように記載できる。
ここで、式(16)の右辺第2項は、ラグランジュの未定定数を通して束縛(磁化ベクトルの発散は0)を課している。
式(16)の右辺第3項は外部からの印加磁場ベクトルが変化したときにすばやく追従させるための減衰項であり、γは減衰定数である。
式(16)の右辺第2項に示される束縛を含んだ運動方程式を解くにあたり、本実施形態では、一般化された束縛の導入法であるSHAKE法を採用する。
束縛を考慮せずに蛙跳び法により式(16)を離散化すると以下の式(17)、式(18)、(19)になる。
ここで、δtは磁化現象の収束計算を行う上で用いる仮想時間刻みである。
添え字のnは任意の整数であり、nδtにおける物理量、n−1/2は(n−1/2)δtにおける物理量、n+1/2は(n+1/2)δtにおける物理量、n+1は(n+1)δtおける物理量に対応している。
磁場解析装置1の制御部3は、ステップ104において既に計算され記憶装置5に記憶されている粒子の位置ベクトルを重心とする多面体要素内のp点、q点の位置ベクトルを読み込む。
次に、制御部3は、ステップ105で既に計算され記憶装置5に記憶されているコイルが磁性体を構成する粒子を重心とする多面体要素内部のp点に作る磁場ベクトルを外部からの印加磁場ベクトルとして読み込む。
さらに、制御部3は、ステップ104で既に計算され、記憶装置5に記憶されている、磁性体を構成する粒子の要素境界面積、法線ベクトルおよび粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
また、制御部3はあらかじめ記憶装置5に記憶されている減衰定数、仮想質量、仮想時間刻みを読み込む。
次に、制御部3はあらかじめ記憶装置5に記憶されている磁場ベクトル、磁場ベクトルの時間微分の初期値を読み込む。
なお、制御部3は後述するステップ107において磁場ベクトル、磁場ベクトルの時間微分が更新されている場合は、その値を読み込む。
なお、磁化ベクトルは、ステップ101より入力情報21として記憶装置5に記憶されている磁化曲線を表す関数に磁場ベクトルを代入することで計算される。
制御部3は、式(17)、式(18)、式(19)を、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内の全てのp点に対して計算し、計算された磁場ベクトルを記憶装置5に記憶する。
以上が、図5のステップ106の詳細である。
次に、磁場解析装置1の制御部3はステップ106で計算し記憶装置5に記憶されている、磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のp点の磁場ベクトルに束縛力を加え、計算された磁場ベクトルを、記憶装置5に記憶する(図5のステップ107)。
ここで、図5のステップ107について、具体的に説明する。
磁性体を構成する粒子の位置ベクトルを重心とする多面体要素内のp点の磁場ベクトルに以下に示す式(20)、式(21)に従って束縛力を加える。
ここで、太字のHは磁場ベクトル、αは仮想質量、δtは仮想時間刻み、γは減衰定数、Nは磁性体を構成する粒子数、sは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のiは磁性体を構成する粒子の内、i番目の粒子の位置ベクトル上での物理量、ipは磁性体を構成する粒子の内、i番目の粒子の位置ベクトルを重心とする多面体要素内のp点の物理量、iqは磁性体を構成する粒子の内、i番目の粒子の位置ベクトルを重心とする多面体要素内のq点の物理量、添え字のnは任意の整数であり、nδtにおける物理量、n+1は(n+1)δtおける物理量に対応している。
磁場解析装置1の制御部3は、演算プログラム19を用いて、ステップ106で計算され記憶装置5に記憶されている磁性体を構成する粒子を重心とする多面体要素内のp点の磁場ベクトルを読み込む。
制御部3は、あらかじめ記憶装置5に記憶されている仮想質量、仮想時間刻み、減衰定数を読み込む。
制御部3は、磁性体を構成する粒子各々に対して式(20)、式(21)に基づく計算を行い、計算された磁場ベクトルを記憶装置5に記憶する。
以上が、図5のステップ107の詳細である。
次に、磁場解析装置1の制御部3はステップ107で求めた磁場ベクトルが束縛条件を満たしているかを判断し、満たしていれば次のステップに進み、満たしていなければステップ107に戻る(図5のステップ108)。
具体的には制御部3は、ステップ107で計算し記憶装置5に記憶されている磁性体を構成する粒子各々の磁場ベクトルから計算される磁化ベクトルを用いて式(22)に基づく計算を行う。
ここで、添え字のn+1は(n+1)δtおける物理量に対応している。
erriは磁性体を構成する粒子の内、i番目の粒子の束縛条件に対する誤差値である。
磁化ベクトル(太字の)Mはステップ107において計算され記憶装置5に記憶されている磁場ベクトルを制御部3が読み込み、ステップ101より入力情報21に記憶されている磁化曲線を表す関数に代入することで計算される。
法線ベクトル(太字の)nには、ステップ104において既に計算され記憶装置5に記憶されているものを制御部3は読み込み代入する。
制御部3は、すべての粒子に対して、誤差の値が式(23)を満たさなければステップ107に戻る。
式(23)においてA1は任意の誤差判別値であり、あらかじめ記憶装置5に任意の値が記憶されている。
以上が図5のステップ108の詳細である。
次に、磁場解析装置1の制御部3は、磁性体の磁化現象が定常状態に到達したかを判断し、条件を満たしていれば次のステップに進む。(図5のステップ109)
ここで、図5のステップ109について、具体的に説明する。
磁性体を構成する粒子が定常状態に到達したかは、以下の式(24)により判断される。
ここで、太字のHは磁場ベクトル、μ0は真空の透磁率、Nは磁性体を構成する粒子の粒子数、sは磁性体を構成する粒子の位置ベクトルを重心とする多面体要素の面数であり、添え字のipは磁性体を構成する粒子の内、i番目の粒子の位置ベクトルを重心とする多面体要素内のp点での物理量である。
さらに、添え字のnは任意の整数であり、nδtにおける物理量、n―1は(n―1)δt、n+1は(n+1)δtおける物理量に対応している。
また、A2は磁性体の磁場ベクトルが定常状態に到達したかを判断するための任意の誤差判定値である。
磁場解析装置1の制御部3は、あらかじめ記憶装置5に記憶されているA2、μ0を読み込む。
また、制御部3は、ステップ107で計算され記憶装置5に記憶されている、磁性体を構成する粒子を重心とする多面体要素内のp点の磁場ベクトルを読み込む。
次に、制御部3は式(24)を計算し、式(24)を満たしていれば次のステップに進み、満たしていなければステップ106に戻る。
以上が図5のステップ109の詳細である。
次に、磁場解析装置1の制御部3は、磁性体を構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁化ベクトル、磁束密度ベクトルと、コイルを構成するすべての粒子の位置ベクトル上の磁場ベクトル、磁束密度ベクトルを計算し、記憶装置5に記憶する(図5のステップ110)。
以下にステップ110を具体的に説明する。
(磁性体を構成する粒子)
磁性体を構成する粒子の位置ベクトル上の磁場ベクトルは以下の式(25)で表される。
ここで、(太字の)ex=(1,0,0)、(太字の)ey=(0,1,0)、(太字の)ez=(0,0,1)である。
太字のHは磁場ベクトル、太字のnは法線ベクトルであり、添え字iは磁性体を構成する粒子の内、i番目の粒子の位置ベクトル上での物理量を表し、添え字ipは磁性体を構成する粒子の内、i番目の粒子の位置ベクトルを重心とする多面体要素内のp点の物理量である。
磁性体を構成する粒子の磁化ベクトルは、磁場ベクトルを磁化曲線を表す関数に代入することで求まる。
磁性体を構成する粒子の位置ベクトル上での磁束密度ベクトルは、式(26)で表される。
ここで、太字のBは磁束密度ベクトル、太字のHは磁場ベクトル、太字のMは磁化ベクトルであり、添え字のiは磁性体を構成する粒子の内、i番目の粒子の位置ベクトル上での物理量を指し示す。
μ0は真空の透磁率である。
磁場解析装置1の制御部3は、ステップ110までで計算されている、磁性体を構成する粒子を重心とする多面体要素内のp点での磁場ベクトルを記憶装置5から読み込む。
次に、制御部3は、ステップ104で既に計算され記憶装置5に記憶されている法線ベクトルを読み込み、式(25)に基づいて磁性体を構成する粒子の位置ベクトル上での磁場ベクトルを計算し記憶装置5に記憶する。
さらに、制御部3は、このステップで計算された磁場ベクトルを、ステップ101より入力情報21として記憶装置5に記憶されている磁化曲線を表す関数に代入し、磁性体を構成する粒子の磁化ベクトルを計算し記憶装置5に記憶する。
次に制御部3は、このステップで計算された磁場ベクトルと磁化ベクトルを式(26)に代入し、磁性体を構成する粒子の位置ベクトル上の磁束密度ベクトルを計算し記憶装置5に記憶する。
なお、真空の透磁率はあらかじめ記憶装置5に記憶されているものを制御部3が読み込む。
(コイルを構成する粒子)
コイルを構成する粒子の位置ベクトル上の磁場ベクトルは式(27)、式(28)で記載される。
ここで、太字のrは位置ベクトル、太字のHは磁場ベクトル、Nは磁性体を構成する粒子の数、sおよびΔSは磁性体を構成する粒子を重心とする多面体要素の境界面数と境界面積、太字のMは磁化ベクトル、太字のnは法線ベクトルであり、添え字のiはコイルを構成する粒子の内、i番目の粒子の位置ベクトル上での物理量、jは磁性体を構成する粒子の内、j番目の粒子の位置ベクトル上での物理量であり、jqは磁性体を構成する粒子の内、j番目の粒子の位置ベクトルを重心とする多面体要素のq点での物理量を指し示す。
式(27)の右辺第一項は、磁性体を構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルを記述する項であり、式(28)で表される。
式(27)の右辺第二項は、コイルを構成する粒子がコイルを構成する粒子の位置ベクトル上に作る磁場ベクトルであり式(1)〜(12)で表される。
磁場解析装置1の制御部3は、ステップ103において計算し記憶装置5に記憶されているコイルを構成する粒子の位置ベクトルを記憶装置5より読み込む。
なお、ステップ113においてコイルを構成する粒子の位置ベクトルが更新されている場合は、制御部3はその値を読み込む。
次に、制御部3はステップ102より計算され記憶装置5に記憶されている磁性体を構成する粒子の数と多面体要素の面数、および粒子の位置ベクトルを重心とする多面体要素内のq点の位置ベクトルと多面体要素境界面積を読み込む。
また、制御部3はこのステップで計算され記憶装置5に記憶されている、磁性体を構成する粒子の磁化ベクトルを読み込む。
さらに、制御部3はステップ103において計算され記憶装置5に記憶されているコイルの寸法およびステップ101で入力情報21として記憶装置5に記憶されている電流密度ベクトルを読み込む。
なお、ステップ116において電流密度ベクトルが更新されている場合は、制御部3はその値を読み込む。
その後、制御部3は式(27)、式(28)および式(1)〜(12)に従い、コイルを構成する粒子の位置ベクトル上での磁場ベクトルを計算し、記憶装置5に記憶する。
制御部3は、コイルを構成する粒子の位置ベクトル上での磁束密度ベクトル(太字の)Biを、すでに計算されているコイルを構成する粒子の位置ベクトル上の磁場ベクトルに真空の透磁率を掛けて計算し記憶装置5に記憶する。
真空の透磁率は、あらかじめ記憶装置5に記憶されている。
以上が、ステップ110の詳細である。
このように、制御部3は、磁場の運動方程式を解くことにより磁性体の磁化現象を解析する。
そのため、有限要素法のように空間全域をメッシュ分割する必要がなく、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる。
また、運動方程式を解くため、行列を扱わず、計算に必要なメモリ量は粒子数に比例する。
次に、磁場解析装置1の制御部3はSPMモータ31を構成する全ての粒子に働く力ベクトルを計算し、記憶装置5に記憶する(図5のステップ111)。
以下に図5のステップ111を具体的に説明する。
(磁性体を構成する粒子に働く力ベクトル)
正準変数を太字のri、太字の傍点付きriとし、式(13)のラグランジアンをラグラ
ンジュの運動方程式に代入する。
すると、磁性体を構成する粒子に働く力は式(29)のように記載される。
ここで、μ0は真空の透磁率、太字のHは磁場ベクトル、太字のMは磁化ベクトル、太字のrは位置ベクトル、ΔVは粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積、Nは磁性体を構成する粒子数であり、添え字のi、jは磁性体を構成する粒子の内、i番目、j番目の粒子の物理量を指し示し、添え字のipは磁性体を構成する粒子の内、i番目の粒子の位置ベクトルを重心とする多面体要素内のp点の物理量を指し示す。
また、φ(ri−rj)はi番目の粒子とj番目の粒子の相互作用ポテンシャルエネルギである。
本実施形態ではSPMモータを剛体(粒子間の相対距離が不変)とする。
剛体モデルとすると、磁性体を構成する粒子に働く力は式(30)のように記載される。
磁場解析装置1の制御部3は、ステップ110までにおいて計算され記憶装置5に記憶されている磁性体を構成する粒子を重心とする多面体要素内のp点での磁場ベクトルを読み込む。
制御部3は、p点での磁場ベクトルを、ステップ101より入力情報21として記憶装置5に記憶されている磁化曲線を表す関数に代入することによりp点での磁化ベクトルを計算し記憶装置5に記憶する。
また、制御部3は、ステップ104で計算され記憶装置5に記憶されている磁性体を構成する粒子の位置ベクトルを重心とする多面体要素境界面の頂点と粒子位置を頂点とする多面体要素体積を読み込む。
次に、制御部3は、磁性体を構成する粒子を重心とする多面体要素内のp点の磁場ベクトルの偏微分を計算し、その値に基づいて、式(30)より力ベクトルを計算し記憶装置5に記憶する。
(コイルを構成する粒子に働く力ベクトル)
コイルを構成する粒子のうち、i番目の粒子の位置ベクトルを太字のri、その位置ベクトルでの電流の単位ベクトルを太字のti、磁束密度ベクトルを太字のBi、電流の流れる方向のコイルの長さをLiとすると、コイルを構成する粒子の内、i番目の粒子に働く力ベクトルは以下の式(31)で記載される。
磁場解析装置1の制御部3は、ステップ103において計算され記憶装置5に記憶されているコイルを構成する粒子の位置ベクトルを読み込む。
なお、ステップ113においてコイルを構成する粒子の位置ベクトルが更新されている場合は制御部3はその値を読み込む。
また、制御部3はステップ103において計算され記憶装置5に記憶されているコイルの寸法を読み込む。
次に、制御部3は、ステップ101より入力情報21として記憶装置5に記憶されているコイルの電流密度ベクトルを読み込む。
なお、制御部3は、ステップ116においてコイルの電流密度ベクトルが更新されている場合はその値を用いる。
制御部3は、式(31)に基づいて、コイルを構成する粒子に働く力ベクトルを計算し記憶装置5に記憶する。
以上が、図5のステップ111の詳細である。
次に、磁場解析装置1の制御部3は、ステップ111で計算した可動部(移動子)を構成する粒子に働く力ベクトルから、可動部(移動子)の回転量、角速度、並進量、並進速度を計算し、記憶装置5に記憶する(図5のステップ112)。
以下に図5のステップ112を具体的に説明する。
剛体の並進運動の運動方程式は式(32)で表される。
ここで、mtotは可動部(移動子)の質量、太字のrcは可動部(移動子)において任意に定めた代表点の位置ベクトル、太字のFは可動部(移動子)を構成する粒子に働く力ベクトル、Nは可動部(移動子)を構成する粒子数である。
次に、任意の固定点の周りの回転運動の運動方程式は式(33−1)、式(33−2)、式(33−3)で表される。
ここで、I1、I2、I3は剛体の主慣性モーメント、ω1、ω2、ω3は慣性主軸座標系での角速度、N1、N2、N3は慣性主軸座標系でのトルクである。
本実施例のSPMモータ31のロータ33はロータシャフト41を通じてある回転軸に固定されている。
固定された回転軸をz軸とすると、ロータ33の運動は以下のように計算される。
ロータ33はロータシャフト41を通じてz軸に固定されているため、z軸周りの回転運動の運動方程式を解くことで、ロータ33の運動は計算される。
z軸周りのロータ33の回転運動の運動方程式は式(34)で表される。
ここで、Izはロータ33の慣性モーメント、ωzはz方向の角速度、Nzはz方向のトルクである。
式(34)を蛙飛び法で離散化すると、以下の式(35)、式(36)を得る。
ここで、θzはz軸周りの回転量であり、dtは微小時間刻み幅である。
また、添え字のn−1/2、n、n+1/2、n+1はそれぞれ(n―1/2)dt、ndt、(n+1/2)dt、(n+1)dt秒での物理量に対応している。
時間刻み幅は記憶装置5にあらかじめ任意の値が記憶されているものを制御部3が読み込む。
回転量、角速度の初期値は記憶装置5にあらかじめ任意の値が記憶装置5に記憶されている。
なお、すでにステップ112においてロータ33の回転量、角速度が更新されている場合はその値を用いる。
制御部3は、ステップ101より入力情報21として記憶装置5に記憶されているロータ33の三次元構造と質量密度により慣性モーメントを計算する。
さらに、制御部3は、ロータ33を構成する粒子の位置ベクトルと、ステップ111で計算し記憶装置5に記憶されているロータ33を構成する粒子に働く力ベクトルに基づいて、式(35)、式(36)を計算し、dt秒後のロータ33の回転量と角速度を計算し記憶装置5に記憶する。
以上が、図5のステップ112の詳細である。
次に、磁場解析装置1の制御部3はステップ112で計算した微小時間後の回転量や重心の移動量を基に、移動子を構成している粒子の位置ベクトルを計算し、記憶する(図5のステップ113)。
次に、磁場解析装置1の制御部3は必要に応じてSPMモータ31を構成する粒子の位置ベクトル、力ベクトル、磁場ベクトル、磁束密度ベクトル、磁化ベクトルをプリンタポート11を介してプリンタ12より出力する(図5のステップ114)。
次に、磁場解析装置1の制御部3は所定の終了条件(時間、移動量等)を満たしているかを判断し、満たしている場合は解析を終了し、満たしていない場合はステップ116に進む(図5のステップ115)。
次に、磁場解析装置1の制御部3は必要に応じてコイル45に流れる電流密度ベクトルを更新し、ステップ104に戻る(図5のステップ116)。
以上が本実施形態に係る磁場解析の手順である。
このように、本実施形態によれば、磁場解析装置1がSPMモータ31の磁場解析を元に、移動子であるロータ33の移動量を求めている。
そのため、本実施形態では磁場と運動の相互作用を演算でき、磁場解析と機構の連成が可能となる。
また、本実施形態では、コイル45に流れる電流密度ベクトルを更新しながら解析を行うことができる。
そのため、コイルに通電する電流を変化させながら解析を行うことができ、従来よりもさらに精度良く磁場解析を行うことができる。
上記した実施形態では、本発明を磁性体を構成する粒子のラグランジアンから導出される磁場の運動方程式を解くことで磁性体の磁化現象を計算した場合について説明したが、本発明は、何等、これに限定されることなく、運動方程式であればラグランジアンから導出されたもの以外のものであってもよい。
また、上記した実施形態では本発明をSPMモータ31の磁場解析に使用したが、本発明は、何等、これに限定されることなく、磁性体が配置された空間における任意の点の磁場を解析するものであれば、例えばリニアモータ等のSPMモータ31以外のモータ、あるいは磁性体で空間を囲む磁気シールドにおいて、磁性体でシールドされた空間内の磁場解析に用いても良い。
上述のように本発明の一実施形態によれば、物体に作用する磁場についての磁場解析と、その磁場により生じる物体の運動についての機構解析との連成解析をするための連成解析装置が提供される。すなわち、この連成解析装置は、磁場解析演算部と、機構解析演算部と、を備える。一実施例においては、制御部3は、磁場解析演算部と、機構解析演算部と、を備えてもよい。磁場解析演算部は、物体を構成する各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式に基づき磁場を演算する。機構解析演算部は、演算された磁場に基づき物体の運動を演算し、物体の各粒子の位置を更新する。磁場解析演算部はさらに、更新された各粒子の位置に基づいて磁場を再計算する。このようにして磁場解析と機構解析とを粒子ごとに繰り返すことにより、連成解析が可能である。なお以下では便宜上、磁場解析演算部を磁場演算部とも称し、機構解析演算部を位置演算部とも称する。
ここで磁場演算部は、上述のように、物体を多数の粒子に分割し、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより物体に作用する磁場を演算する。解析対象物体は磁化曲線が非線形である材料を含んでもよく、例えば磁性体を解析対象物体としてもよい。本願発明者は、いわゆる解析力学の手法を用いて、物体の各粒子に作用する磁場の支配方程式を運動方程式の形式で記述することに成功した。上述の一実施例においては、上記の式(13)に定義されるラグランジアンをラグランジュの運動方程式に代入することにより、式(16)のように各粒子に作用する磁場を運動方程式の形式で記述している。
磁場演算部は、磁場の運動方程式を数値的に解くことにより、物体に作用する磁場を求める。具体的には例えば、磁場演算部はまず、離散化された磁場の運動方程式に磁場及びその時間微分の初期値を代入し、微小時間経過後の磁場及びその時間微分を求める。ここで、微小時間とは磁場の収束計算を行うために用いられる仮想的な時間であってもよい。つまりこの仮想時間は実際の経過時間とは異なる。磁場演算部は、演算された微小時間経過後の磁場が収束判定基準を満たすか否かを判定する。収束判定条件は、一実施例においては式(24)で与えられる。磁場演算部は、収束判定基準が満たされるまで、離散化された磁場の運動方程式を繰り返し演算することにより磁場の数値解を求める。
この方法によれば、有限要素法に基づく周知の解析法とは異なり逆行列の演算が不要であるので、少ない計算資源で磁場解析を実行することができるという点で好ましい。また、この方法は、高精度の解析結果を、有限要素法のように解析対象物体とその周囲の空間とを含む系全体をメッシュ分割することなく得られるという点でも好ましい。
一実施例においては、上述の連成解析がモータの磁場解析に利用される。本願発明者の検証によると、モータの磁場解析に要求される高精度のレベルに収束判定条件を設定したときに、演算終了までの所要時間が比較的長くなる場合があることがわかった。この収束の遅れは、解析対象物体における磁気飽和に起因すると考えられる。すなわち、磁化曲線の非線形性が磁場の数値解の収束性に影響している。
図10は、磁性体(例えば鉄)の磁化曲線の一例を模式的に示す図である。図10の横軸は磁場Hを示し、縦軸は磁化Mを示す。図10に示されるように、磁場Hが小さいときには磁化Mは磁場Hに比例する。磁場Hが大きくなるにつれて、磁化Mと磁場Hとの比例定数(いわゆる磁気感受率)は徐々に小さくなっていく。磁場Hが小さい場合の比例定数は(例えば図中の点A)、磁場Hが大きい場合(図中の点B)の比例定数の1000倍程度の大きさとなり得る。このように、磁性体の磁化曲線は非線形である。
磁場の運動方程式には、磁場及びその時間微分項だけではなく、磁化に依存する磁化依存項が含まれる。ある粒子についての磁場の運動方程式における磁化依存項は例えば、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項である。このような相互作用を表す項は、例えば式(14)に示されるように、磁化Mを陽に含む。したがって、式(16)に示される磁場の運動方程式は、収束計算を進めるにつれて非線形に変化する項を含むと言える。収束計算の磁場の初期値は通常小さい値(例えばゼロ)にとるから、式14のH(rip)は、収束計算の当初に比べて完了時には1000分の1程度の小さい値となり得る。方程式におけるこの大きな非線形性によって数値解の収束に遅れが生じると考えられる。
磁場の運動方程式には、任意に設定することが許容される調整可能パラメタを導入することができる。例えば式(16)に示される磁場の運動方程式においては、仮想質量α及び減衰定数γは調整可能パラメタである。パラメタの設定作業を簡単にすることを重視する場合には、これらの調整可能パラメタを一定値に設定してもよい。この場合、収束計算中に調整可能パラメタは不変である。一方、磁場演算部は、収束計算中に調整可能パラメタを変化させてもよい。
磁場の運動方程式は一般に、束縛を考慮に入れなければ、2階の線形常微分方程式の形式をとる。この形式は、1自由度の振動の基礎方程式と共通である。よって、磁場の運動方程式の収束特性は、振動の基礎方程式の減衰特性に対応する。このことから、振動の方程式の各項の係数に応じて減衰特性が変化することに準じて、磁場の運動方程式の調整可能パラメタを変更することにより所望の収束特性に調整することが可能である。
そこで、本発明の一実施形態においては、磁場演算部は、磁場の運動方程式に導入されている調整可能パラメタを、磁化曲線の非線形性が数値解の収束性に与える影響を緩和するように調整してもよい。磁場演算部は、磁化依存項による数値解の収束性への影響を抑えるように調整可能パラメタを変化させてもよい。一実施例においては、磁場演算部は、収束判定基準を満たすか否かを判定する収束判定処理において収束判定基準が満たされていないと判定したときに、調整可能パラメタを変化させる係数調整処理を実行してもよい。上述の実施例においては、制御部3は、誤差の値が誤差判別値A2を超えると判定した場合に(図5のステップ109のNo)、仮想質量α及び減衰定数γを調整する係数調整ステップを実行してもよい。
ある物体の粒子に作用する磁場Hの運動方程式は一般に、式(37)で表される。ここで、係数CMは磁場Hの2階時間微分項の係数を表し、係数CDは磁場Hの1階時間微分項の係数を表し、係数CKは磁場Hの項の係数を表す。係数CM、CD、CKはいずれも解析者が調整することのできるパラメタである。一方、係数FNLは、磁場Hと磁化Mとの比例定数に依存する。H0は外部からの印加磁場を表す。
よって、調整可能なパラメタである係数CM、CDを磁化曲線に基づいて調整することにより、磁化曲線の非線形性が数値解の収束性に与える影響を排除することができる。すなわち、磁場Hの時間微分項の係数CM及びCDを、式(38)に示されるように磁化依存係数FNLを用いてCM’FNL及びCD’FNLと定義する。磁化依存係数FNLの時間変化率が十分に小さければ、磁化曲線の非線形性が数値解の収束性に与える影響を事実上排除することができる。
式(38)の磁場の項及びその時間微分項には、磁場Hに対し非線形に変化する共通パラメタFNLが導入されている。つまり、共通の磁化依存係数FNLをもたない運動方程式(39)の収束特性に、式(38)に示される磁場の運動方程式の収束特性は実質的に等しくなる。
式(39)の磁場H及びその時間微分項の係数CM’、CD’、及びCKはいずれも磁化Mに依存せず、解析者が任意に設定することのできるパラメタである。よって、式(38)において磁化依存係数FNLが毎回の収束計算ごとに変化しても、係数CM’、CD’、及びCKから定まる一定の収束特性を得ることができる。なおこの場合、所望の収束特性をもたらすよう係数CM’、CD’、及びCKを予め設定することが好ましい。
なお本発明は、磁化曲線の非線形性の影響を磁場の運動方程式の収束特性から排除する実施例には限られない。基礎方程式の変数の係数がその変数に対して非線形に変化する場合における運動方程式型の数値解法にも適用可能である。例えばいわゆる大変形解析にも適用可能であろう。
一実施例においては、第1の物理量と、該第1の物理量と非線形の関係をもつ第2の物理量とを有する物体の特性を解析するための解析装置が提供される。例えば上述の実施例では第1の物理量及び第2の物理量はそれぞれ磁場及び磁化に対応する。この解析装置は、運動方程式の形式で記述された第1の物理量が満たすべき関係式を数値的に解く数値演算部を備える。この関係式は、第2の物理量に依存する項と、第2の物理量に依存しない項と、を含み、第2の物理量に依存しない項の係数が第2の物理量に依存する項に連動するよう定義されている。一実施例においては、第1の物理量及びその時間微分項のそれぞれの係数には第2の物理量に連動する連動部分が共通に導入されていてもよい。
また、一実施例においては、調整可能パラメタは、磁場の運動方程式を数値的に解く過程で導出される値を用いて定義されていることが好ましい。つまり、磁場演算部は、磁場演算処理において算出された値を用いて係数調整処理を実行してもよい。このようにすれば、磁場の運動方程式を解く過程で得られた計算済みの値を流用して調整可能パラメタを調整することができる。調整可能パラメタの設定のために新たな演算をする必要がなくなるので、計算コストを低減することができる。よって、収束演算の所要時間を短くすることもできる。
そこで、一実施例においては、調整可能パラメタは、他の粒子の磁化に起因して当該粒子に作用する磁場を表す項に連動して調整されてもよい。式(16)を解くことにより磁場を求める上述の実施例においては例えば、仮想質量αipを式(40)及び式(41)により可変量として定義する。係数kipは磁場の運動方程式に存在するH(rip)の項を用いて定義されているから、調整可能パラメタの設定のために新たな演算をする必要がない。なお係数α0は解析者が任意に設定することのできるパラメタである。
この仮想質量αipにより、磁化曲線の非線形性の影響が磁場の運動方程式の収束特性から排除されることを説明する。上述のようにSHAKE法を適用して式(16)を解くことを考慮すると、式(16)の右辺第2項は0に収束するから、式(42)のように表される。
式(14)及び式(41)から、H(rip)は式(43)で表される。
式(40)及び式(43)を式(42)に代入することにより、磁場の運動方程式は式(44)で表される。式(44)の磁場の項及びその時間微分項には、磁場Hに対し非線形に変化する共通パラメタ(1−kip)が導入されている。よって、磁場の運動方程式(44)の収束特性は、共通パラメタ(1−kip)を有しない式(45)の収束特性に実質的に等しくなる。よって、磁気飽和に起因して係数kip(式(16)で言えばH(rip))が収束計算中に変化しても、磁場の運動方程式(44)の収束特性は変化しない。係数α0、γを予め設定することにより、所望の収束特性を実現することができる。
なお、上述のパラメタの更新は、磁場解析の対象物体が磁気感受率の異なる物質を含む場合にも有効であると考えられる。例えば物体が磁石と鉄とを含む場合にも有効である。物質が異なることによる粒子ごとの磁気感受率の違いは、磁気飽和が生じているか否かによる磁気感受率の違いと同様に上記の磁場の運動方程式において取り扱われているからである。
図11は、本発明の一実施形態に係る磁場演算処理を説明するためのフローチャートである。図11に示されるように、磁場解析装置1の制御部3はまず、磁場の運動方程式を作成する(S10)。一実施例においては、制御部3は、図5に示されるステップ101乃至105を実行することにより磁場の運動方程式(16)を作成する。ここで、仮想質量αipは式(40)及び式(41)により可変量として定義されている。
次に制御部3は、磁場の運動方程式に基づいて磁場を演算する(S11)。一実施例においては、制御部3は、図5に示されるステップ106乃至108を実行することにより、磁場の運動方程式(16)に基づいて磁場を演算する。制御部3は、演算された磁場が収束判定基準を満たすか否かを判定する(S12)。収束判定基準は例えば図5に示すステップ109と同様であり、制御部3は、誤差が誤差判別値を満たすか否かを判定する。
制御部3は、収束判定基準が満たされていないと判定した場合には(S12のNo)、調整可能パラメタを更新して(S13)、磁場の運動方程式に基づく磁場演算を再度行う(S11)。制御部3は、収束判定基準が満たされるまで、調整可能パラメタの更新と磁場演算とを繰り返す。制御部3は、調整可能パラメタとして仮想質量αipを更新する。すなわち、制御部3は係数kipを更新する。
制御部3は、収束判定基準が満たされたと判定した場合には(S12のYes)、演算された磁場に基づいて解析対象物体の運動を演算する(S14)。一実施例においては、制御部3は、図5に示されるステップ110乃至114を実行する。次に制御部3は、所定の終了条件が満たされているか否かを判定し(S15)、満たされていると判定した場合には解析を終了する(S15のYes)。制御部3は、終了条件が満たされていないと判定した場合には(S15のNo)、運動による解析対象物体の新たな位置に基づいて再度磁場を演算する。このようにして、制御部3は、物体に作用する磁場とそれによる物体の運動との連成解析を実行する。
図12及び図13は、本発明の一実施形態による磁場演算の例を示す図である。図12は仮想質量αipを一定値に固定した場合の磁場Hと仮想時間δtとの関係を示し、図13は上述のように係数kipを用いて仮想質量αipを更新しながら収束計算をしたときの磁場Hと仮想時間δtとの関係を示す。図12及び図13の磁場演算の例は、1次元に配列された10個の粒子からなる物体に対して行ったものである。粒子1及び粒子10が物体の両端である。図12の演算結果に比べて図13の演算結果は磁場が短時間で定常状態に達していることがわかる。仮想質量αipを更新しながら収束計算を行うことにより、磁場の運動方程式の数値解の磁気飽和に伴う収束の遅れを低減することができる。