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