JP2011254678A - 解析装置および解析方法 - Google Patents

Abstract

【課題】磁場・機構・弾性・制御を連成した解析を実現する。
【解決手段】解析装置１は、物体を複数の粒子からなる粒子系として記述した上で当該物体の振る舞いを解析する。解析装置１は、粒子間の相互作用による相互作用ポテンシャルエネルギを演算する相互作用演算部６４と、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより各粒子に作用する磁場を演算する磁場演算部６６と、相互作用演算部６４によって演算された相互作用ポテンシャルエネルギと磁場演算部６６によって演算された磁場とに基づき各粒子の運動を演算し、各粒子の位置を更新する機構・弾性演算部６８と、を備える。
【選択図】図１

Description

本発明は、解析装置および解析方法に関する。

近年、モータなどの電気機器の設計開発の現場では、コンピュータの計算能力の向上に伴い、磁場解析を取り入れたシミュレーションがよく使用されるようになっている。シミュレーションを使用すると、実際にプロトタイプを製作しなくてもある程度の評価が可能となるので、設計開発のスピードが向上しうる。

例えば特許文献１には、磁場解析を実行する演算処理装置を具えるモータ解析装置が記載されている。演算処理装置は、ユーザの操作に基づく外部指令に応じて、有限要素法による静磁場解析やマクスウェルの応力法によるトルク計算を実行する。有限要素法による静磁場解析のためにメッシュ分割が行われる。メッシュ分割は、コア領域及びハウジング領域、さらには外部空気層領域についても行われる。

また、繰り込み群分子動力学を使用したシミュレーションの手法が知られている（特許文献２参照）。

特開平１１−１４６６８８号公報 特開２００６−２８５８６６号公報

上記特許文献１のように磁場解析に有限要素法を使用する場合、空間全域をメッシュ分割する必要がある。このため、機構との連成を行う場合、複雑に動くような系に対しては、ステップ毎にリメッシュを実施しなくてはならない。また、有限要素法では接触問題を取り扱うことは困難である。さらに、磁気エネルギに起因する変形、この変形に伴う磁気回路の変化、この変化による変形量の変動といった一連の現象を計算することも困難である。

本発明はこうした課題に鑑みてなされたものであり、その目的は、高精度なシミュレーションを実現する解析技術の提供にある。

本発明のある態様は解析装置に関する。この解析装置は、物体を複数の粒子からなる粒子系として記述した上で当該物体の振る舞いを解析する解析装置であって、粒子間の相互作用による相互作用ポテンシャルエネルギを演算する相互作用演算部と、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより各粒子に作用する磁場を演算する磁場演算部と、相互作用演算部によって演算された相互作用ポテンシャルエネルギと磁場演算部によって演算された磁場とに基づき各粒子の運動を演算し、各粒子の位置、速度を更新する機構・弾性演算部と、を備える。

この態様によると、相互作用演算と磁場演算と機構・弾性演算とを連成することができる。

なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。

本発明によれば、高精度なシミュレーションを実現できる。

実施の形態に係る解析装置の機能および構成を示すブロック図である。 ＳＰＭモータの構成の概略を示す図である。 ＳＰＭモータのステータティースの１つの周囲の拡大斜視図である。 図１の解析装置における一連の処理を示すフローチャートである。 コイルの説明図である。 直方体導体とローカル座標との関係を示す図である。 円弧状柱状導体とローカル座標との関係を示す図である。 図４に示される処理の一部の詳細を示すフローチャートである。 磁性体の磁場計算に用いる要素の説明図である。

以下、本発明を好適な実施の形態をもとに図面を参照しながら説明する。各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付する場合がある。また、適宜重複した説明は省略する。

実施の形態に係る解析装置では、モータなどの現実の物体に対して、磁場・機構・弾性・制御の４つの連成計算によるシミュレーションが行われる。特に、機構・弾性に関しては分子動力学、磁場に関しては磁場の運動方程式を用いて連成計算が実施される。これにより、従来技術では検討が困難であった現象の詳細な分析が可能となる。

図１は、実施の形態に係る解析装置１の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。

解析装置１は、解析対象の物体を複数の粒子からなる粒子系として記述した上で当該物体の振る舞いを解析する。ここでの粒子系としては、特に繰り込まれた粒子系が使用される。繰り込みについては、特許文献２に詳しい。
解析装置１は、制御部３と、記憶装置５と、メディア入出力部６と、入力部７と、表示部９と、プリンタポート１１と、を備え、これらの部材はバス１３を介して互いに接続されている。

制御部３は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等で構成され、記憶手段としての記憶装置５に格納されたプログラムに従って、バス１３を介して接続された各装置を駆動制御する。

記憶装置５は、初期条件を有する情報である入力情報を有している。メディア入出力部６は、フロッピー（登録商標）ディスク、ＣＤ、ＤＶＤ等のメディアとの間で情報の入出力を行う装置である。入力部７は、キーボード、マウス等の入力装置であり、表示部９はディスプレイ等の表示機器である。プリンタポート１１には出力装置としてのプリンタ１２等が接続される。

制御部３は、粒子モデル生成部６０と、繰り込み部６２と、相互作用演算部６４と、磁場演算部６６と、機構・弾性演算部６８と、制御パラメータ更新部７０と、終了判定部７２と、を含む。各部の詳細は後述する。

以下では、解析装置１による解析対象の物体としてモータ、特に永久磁石モータ（Permanent Magnet Motor）の一種であるＳＰＭモータ３１（Surface Permanent Magnet Motor）を採用する場合を説明する。

まず、ＳＰＭモータ３１の構成の概略を図２および図３を参照して説明する。
図２に示すように、ＳＰＭモータ３１は回転子（移動子）であるロータ３３と固定子であるステータ３５を有している。ロータ３３は鉄等の磁性体である円柱状のロータコア３７を有し、ロータコア３７の表面には永久磁石３９が設けられている。ロータコア３７の軸中心には棒状のロータシャフト４１が設けられている。

ステータ３５は磁性体である歯状のステータティース４３とステータティース４３の外側に設けられた円筒状の磁性体であるコアバック４４、コアバック４４の外側に設けられた円筒状のフレーム４６から構成されている。

図２および図３に示すように、ステータティース４３には、金属等の導電体であるコイル４５が巻きつけられている。なお、実際のＳＰＭモータ３１ではコイル４５は仕様に応じたターン数でステータティース４３に巻きつけられて束となっているが、本実施の形態では、図２および図３に描かれているように、コイル一本一本をモデル化せず、コイルの束を１つの導体として扱う。

このような構造のＳＰＭモータ３１は、永久磁石３９の磁場、およびコイル４５に電流を流すことにより発生する磁場によって、ロータ３３、ステータ３５が磁化する。磁性体の磁気エネルギの偏差によりＳＰＭモータ３１は駆動する。

そのため、ＳＰＭモータ３１の解析を行うためにはコイル４５、永久磁石３９がロータ３３、ステータ３５を構成する磁性体上に作る磁場ベクトルを計算し、これら磁性体の磁化現象を解析する必要がある。

図４は、実施の形態に係る解析装置１における一連の処理を示すフローチャートである。以下の手順において、磁性体とは、ロータ３３、ステータ３５を構成する磁性体と永久磁石３９を指し示し、磁化曲線を表す関数によってこれらは区別される。導体とは、例えばコイル４５である。

（ステップＳ２０２）
解析装置１は、ＳＰＭモータ３１の形状に関する情報や材料の特徴（材料定数など）を含む初期条件を取得し、記憶装置５に入力情報として記憶する。この初期条件は例えばメディア入出力部６を介してＣＤ−ＲＯＭ等の記録媒体から読み込んだものであってもよい。さらに、あらかじめ上記初期条件が入力情報として記憶されている場合は、ステップＳ２０２は不要である。

初期条件は、ＳＰＭモータ３１の三次元構造（形状、座標点）、質量密度、磁性体の磁化曲線を表す関数、導体の電流密度ベクトル、を含む。ＳＰＭモータ３１の三次元構造の情報とは例えばＣＡＤ等のデータである。電流密度ベクトルは、コイル４５の作る磁場ベクトルを計算する際に必要になる。磁性体の磁化曲線を表す関数は磁化ベクトルを計算する際に必要になる。

（ステップＳ２０４）
粒子モデル生成部６０は、記憶装置５が記憶する初期条件からＳＰＭモータ３１の粒子モデルを生成する。粒子モデル生成部６０は、初期条件に含まれる三次元構造の情報から、ＳＰＭモータ３１をＮ個の粒子に分割し（Ｎは２以上の整数）、各粒子の位置ベクトルを計算し、記憶装置５に記憶する。Ｎ個の粒子は粒子系Ｓを構成する。粒子は、原子、分子単位であってもよいし、粒子の位置ベクトルを重心とする多面体要素であってもよい。

粒子の数Ｎおよび位置ベクトルは初期条件に含まれる三次元構造およびあらかじめ記憶装置５に記憶されている多面体要素の形状により計算される。以下ではＳＰＭモータ３１の磁性体に対応する粒子は、粒子の位置ベクトルを重心とする立方体要素であるとする。

また、導体であるコイル４５に対応する粒子は、以下のように決定する。
図５は、コイル４５の説明図である。図５（ａ）は、コイル４５の拡大斜視図であり、図５（ｂ）は、コイル４５をローカル導体に分割した例を示す図である。図６は、直方体導体とローカル座標との関係を示す図である。図７は、円弧状柱状導体とローカル座標との関係を示す図である。

粒子モデル生成部６０は、初期条件に含まれる三次元構造の情報から、図５（ａ）に示されるコイル４５を図５（ｂ）に示すようにローカル導体（直方体導体４５ａ、４５ｂと円弧状柱状導体４５ｃ、４５ｄ）に分割し、それぞれの導体が作る磁場ベクトルを計算するための係数を計算し、記憶装置５に記憶する。

まず、ローカル導体が直方体導体４５ａの場合について説明する。
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。

ローカル導体が直方体導体４５ａの場合は、粒子モデル生成部６０は、図６に示すように、ローカル座標系（ｘ_ｓ,ｙ_ｓ,ｚ_ｓ）を適用する。このローカル座標系においてはローカル導体（直方体導体４５ａ）の重心を原点Ｏ_ｓとし、直方体導体４５ａの寸法はｘ_ｓ方向に２ａ、ｙ_ｓ方向に２ｂ、ｚ_ｓ方向に２ｃの長さを持つものとする。また原点Ｏ_ｓに粒子は位置するものとする。

粒子モデル生成部６０は、初期条件に含まれるコイル４５の三次元構造を読み込み、ローカル導体（直方体導体４５ａ）の寸法であるａ、ｂ、ｃと粒子位置ベクトルとを計算し記憶装置５に記憶する。

次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。

ローカル導体が円弧状柱状導体４５ｃの場合は、粒子モデル生成部６０は、図７に示すようにローカル座標系（ｘ_ｃ,ｙ_ｃ,ｚ_ｃ）を適用する。このローカル座標系においては原点Ｏ_ｃは円弧の中心軸上に存在し、かつ円弧状柱状導体４５ｃの高さ方向（図７のz_ｃ方向）に対して円弧状柱状導体４５ｃが対称となる点に存在するものとする。

また、ｘ_ｃ,ｙ_ｃ,ｚ_ｃは、ｘ_ｃ−ｙ_ｃ平面でみると、＋ｘ_ｃ軸を基点とし、円弧状柱状導体４５ｃの円弧が＋ｚ_ｃ軸からみて反時計回りになるようして決定する。円弧状柱状導体４５ｃの内径と外径の平均値をＲ_ｃとし、径方向の厚さを２ｒ_ａ、ｚ_ｃ方向の高さを２ｚ_ｂとする。

電流は＋ｘ_ｃ軸を基点として、＋ｚ_ｃから見て反時計回りの方向への角度をθとし、電流はこの方向に一様な電流密度ｊで流れているものとする。粒子は円筒座標系で（Ｒ_ｃ、θ／２、０）に位置するものとする。

粒子モデル生成部６０は、初期条件に含まれるコイル４５の三次元構造を読み込み、ローカル導体（円弧状柱状導体４５ｃ）の寸法であるｒ_ａ、ｚ_ｂ、θ、Ｒ_ｃおよび粒子位置ベクトルを計算し記憶装置５に記憶する。
なお、粒子モデル生成部６０は磁性体、導体の別なく同じように粒子を設定してもよい。

（ステップＳ２０６）
繰り込み部６２は、特許文献２で説明されている繰り込み群分子動力学の手法を使用して粒子系Ｓを繰り込み処理し、繰り込まれた粒子系Ｓ’を生成する。以下では、第１の繰り込み因子をα、第２の繰り込み因子をγ＝０、第３の繰り込み因子δ＝２、空間の次元数ｄ＝３とする。特許文献２によると、粒子系Ｓのパラメータと繰り込まれた粒子系Ｓ’のパラメータとの間には、
の関係が成り立つ。また、本発明者による独自の考察（後述）によると、
の関係が成り立つ。

（ステップＳ２０８）
相互作用演算部６４は、繰り込まれた粒子系Ｓ’について、粒子間の相互作用による相互作用ポテンシャルエネルギφ’を演算する。相互作用ポテンシャルエネルギφ’の演算については特許文献２に詳しい。特に本実施の形態では、相互作用ポテンシャルエネルギφ’は弾性を考慮した形を有する。

（ステップＳ２１０）
磁場演算部６６は、繰り込まれた粒子系Ｓ’について、各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより各粒子に作用する磁場を演算する。詳細は後述する。

（ステップＳ２１２）
機構・弾性演算部６８は、相互作用演算部６４によって演算された相互作用ポテンシャルエネルギφ’と磁場演算部６６によって演算された磁場とに基づき各粒子の運動を演算し、各粒子の位置、速度を更新する。詳細は後述する。

（ステップＳ２１４）
制御部３は、ユーザからの出力指示の有無を確認する。
（ステップＳ２１６）
制御部３は、ユーザからの出力指示がある場合（ステップＳ２１４のＹ）、繰り込まれた粒子系Ｓ’にステップＳ２０６に対応するリスケーリングを施す。
（ステップＳ２１８）
制御部３は、リスケーリングの結果得られる位置ベクトル、力ベクトル、磁場ベクトル、磁束密度ベクトル、磁化ベクトルなどをプリンタポート１１を介してプリンタ１２より出力する。

（ステップＳ２２０）
終了判定部７２は、ユーザからの出力指示がない場合（ステップＳ２１４のＮ）、所定の終了条件（時間、移動量等）を満たしているかを判断する。終了判定部７２は、終了条件が満たされている場合は（ステップＳ２２０のＹ）、解析を終了する。
制御部３は、終了条件が満たされていない場合は（ステップＳ２２０のＮ）、処理をステップＳ２０８に戻す。すなわち、制御部３は、更新された各粒子の位置、速度を使用して、再度相互作用ポテンシャルエネルギφ’や磁場を演算する。

なお、制御パラメータ更新部７０は、終了条件が満たされていない場合であって（ステップＳ２２０のＮ）ユーザからの求めがある場合は、繰り込まれた粒子系Ｓ’の外部から与えられる制御パラメータであって、各粒子に作用する磁場に影響を及ぼす制御パラメータ、例えば電流密度ベクトルを更新する。磁場演算部６６は更新された電流密度ベクトルを使用して磁場を演算する。

ステップＳ２１０およびステップＳ２１２について以下に詳述する。
図８は、図４のステップＳ２１０およびステップＳ２１２の詳細を示すフローチャートである。

（ステップＳ１０２）
磁場演算部６６は、繰り込まれた粒子系Ｓ’におけるコイルの寸法および電流密度ベクトルを用いて、ビオ・サバールの法則を積分することにより得られる解析解により、通電されたコイルが磁性体に対応する粒子の位置ベクトル上に作る磁場ベクトルを計算し、記憶装置５に記憶する。

図９は、磁性体の磁場計算に用いる要素の説明図である。繰り込まれた粒子系Ｓ’において、分割した立方体要素を２次元表示すると図９に示す形状となる。
ここで、粒子の位置ベクトルを（太字の）ｒ_ｇ’、粒子の位置ベクトルを重心とする立方体要素の要素境界面の中点をｑ’点とし、粒子の位置ベクトルとｑ’点との中間点にｐ’点を定義する。

図８に戻る。
まず、ローカル導体が直方体導体４５ａの場合について説明する。
なお、ローカル導体が直方体導体４５ｂの場合は、ローカル導体が直方体導体４５ａの場合と同様であるため、説明を省略する。

ローカル導体が直方体導体４５ａの場合は、図６に対応して繰り込まれた粒子系Ｓ’で定義されるローカル座標系（ｘ_ｓ’,ｙ_ｓ’,ｚ_ｓ’）を適用する。
次に、任意の点であるＰ点の位置ベクトルをローカル座標系（ｘ_ｓ’,ｙ_ｓ’,ｚ_ｓ）に変換する。

通電された直方体導体がＰ点に作る磁場ベクトルは、以下に示す式（１）〜（３）で記載される。

ここで、（太字の）ｒ_ｐｓ’はローカル座標系（ｘ_ｓ’,ｙ_ｓ’,ｚ_ｓ’）でのＰ点の位置ベクトルであり、ｘ_ｐｓ’,ｙ_ｐｓ’,ｚ_ｐｓ’はｘ_ｓ’,ｙ_ｓ’,ｚ_ｓ’方向の値である。πは円周率である。Ｈ_ｘ’ｓ’、Ｈ_ｙ’ｓ’、Ｈ_ｚ’ｓ’はローカル座標系（ｘ_ｓ’,ｙ_ｓ’,ｚ_ｓ’）における磁場ベクトルの各成分である。ｊ’は電流密度である。

また、ｘ_ｉ’,ｙ_ｊ’,ｚ_ｋ’はｘ_ｓ’,ｙ_ｓ’,ｚ_ｓ’方向の積分の上限、下限を表しており、式（４）に示す関係が成立する。

ここで、ａ’，ｂ’，ｃ’は繰り込み後の直方体導体の寸法である。

次に、ローカル導体が円弧状柱状導体４５ｃの場合について説明する。
なお、ローカル導体が円弧状柱状導体４５ｄの場合は、ローカル導体が円弧状柱状導体４５ｃの場合と同様であるため、説明を省略する。

ローカル導体が円弧状柱状導体４５ｃの場合は、図７に対応して繰り込まれた粒子系Ｓ’で定義されるローカル座標系（ｘ_ｃ’,ｙ_ｃ’,ｚ_ｃ’）を適用する。
ローカル座標系（ｘ_ｃ’,ｙ_ｃ’,ｚ_ｃ’）に変換後のＰ点の位置ベクトルを以下の式（５）に示すように円筒座標系（太字の）ｒ_ｐｃ’＝（Ｒ_ｐｃ’、φ_ｐｃ’、Ｚ_ｐｃ’）に変換する。
通電された円弧状柱状導体４５ｃがＰ点に作る磁場ベクトルは、以下の式（６）〜（１１）で表される。

ここで、Ｈ_ｒｃ’、Ｈ_ｔｃ’、Ｈ_ｚｃ’は円筒座標系での磁場ベクトルの各成分である。ｊ’は電流密度である。ｒ_ａ’、θ’、ｚ_ｂ’は円弧状柱状導体４５ｃの寸法であり、Ｒ_ｃ’は円弧の内径と外径の平均値である。ｓｇｎはＺ_ｋ’の符号であり、Ｒ_ｊ’、Ｚ_ｋ’は積分の上限、下限を表しており、式（１２）に示す関係が成立する。

以上が任意の点Ｐにコイル４５が作る磁場ベクトルの計算手順である。
磁場演算部６６は、繰り込まれた粒子系Ｓ’における、磁性体に対応するすべての粒子の位置ベクトルを重心とする立方体要素内部のｐ’点の位置ベクトルを読み込む。

次に、磁場演算部６６は、図９のｐ’点の位置ベクトルをローカル座標系（ｘ_ｓ’,ｙ_ｓ’,ｚ_ｓ’）に変換し、直方体導体４５ａの繰り込まれた粒子系Ｓ’での寸法と、繰り込まれた粒子系Ｓ’での電流密度ベクトルとを読み込み、式（１）〜（４）に基づいて磁場ベクトルを計算する。

次に、磁場演算部６６は、式（１）〜（４）で計算された磁場ベクトルをグローバル座標系（ｘ’,ｙ’,ｚ’）に変換し、変換後の磁場ベクトルを記憶装置５に記憶する。
これにより、通電された直方体導体４５ａが、ｐ’点に作る磁場ベクトルが求められる。

さらに磁場演算部６６は、ｐ’点の位置ベクトルをローカル座標系（ｘ_ｃ’,ｙ_ｃ’,ｚ_ｃ’）に変換し、さらに円筒座標系（太字の）ｒ_ｐｃ’＝（Ｒ_ｐｃ’、φ_ｐｃ’、Ｚ_ｐｃ’）に変換する。

次に磁場演算部６６は、繰り込まれた粒子系Ｓ’での円弧状柱状導体４５ｃの寸法ｒ_ａ’、θ’、ｚ_ｂ’およびＲ_ｃ’を読み込み、また、磁場演算部６６は繰り込まれた粒子系Ｓ’での電流密度ベクトルを読み込む。

磁場演算部６６は式（６）〜（１２）に基づいて磁場ベクトルを計算する。
次に、磁場演算部６６は、計算された磁場ベクトルを直交座標系に変換し、さらにグローバル座標系（ｘ’,ｙ’,ｚ’）に変換し記憶装置５に記憶する。
これにより、通電された円弧状柱状導体４５ｃがｐ’点に作る磁場ベクトルが求められる。

（ステップＳ１０４）
磁場演算部６６は、繰り込まれた粒子系Ｓ’において、永久磁石３９が磁性体に対応する粒子の位置ベクトル上に作る磁場ベクトルを計算し、記憶装置５に記憶する。

（ステップＳ１０６）
磁場演算部６６は、ステップＳ１０２で計算した磁場ベクトルとステップＳ１０４で計算した磁場ベクトルとの和を計算する。

（ステップＳ１０８）
磁場演算部６６は、繰り込まれた粒子系Ｓ’において、磁性体に対応する粒子の位置ベクトルと、繰り込まれた粒子系Ｓ’における立方体要素の形状により、磁性体に対応する粒子の磁場の運動方程式の係数を計算し記憶装置５に記憶する。

磁場演算部６６は、繰り込まれた粒子系Ｓ’における磁性体に対応する粒子（粒子を重心とする立方体要素）の位置ベクトルを読み込む。
なお、ステップＳ２１２おいて粒子の位置ベクトルが更新され記憶装置５に記憶されている場合は、磁場演算部６６はその値を読み込む。

磁場演算部６６は、粒子の位置ベクトルと繰り込まれた粒子系Ｓ’における立方体要素の形状から、磁性体に対応するすべての粒子の位置ベクトルを重心とする立方体要素内のｐ’点とｑ’点の位置ベクトルと要素境界面への法線ベクトル（太字の）ｎ’を計算し、記憶装置５に記憶する。

次に、磁場演算部６６は、磁性体に対応する粒子を重心とする立方体要素の境界面積ΔＳ’、および粒子の位置ベクトルを重心とする立方体要素の境界面の頂点と粒子位置を頂点とする立方体要素体積ΔＶ’を計算し、記憶装置５に記憶する。

次に、磁場演算部６６は、磁性体に対応する粒子の磁場の運動方程式から拘束条件を考慮せずに仮想時間刻みδｔ’後の磁場ベクトルを計算し、記憶装置５に記憶する。

磁性体に対応する（Ｎ’−ｗ）個の粒子のラグランジアンを式（１３）〜式（１５）で表される形とする。ｗは自然数であり、繰り込まれた粒子系Ｓ’においてコイル４５に対応する粒子の数である。

ここで、式（１３）において、α’は仮想質量、太字のｒ’は位置ベクトル、太字のＨ’は磁場ベクトル、太字の傍点付きＨ’は磁場ベクトルの時間微分、太字のＭ’は磁化ベクトル、太字のｎ’は粒子の位置ベクトルを重心とする立方体要素境界面の法線ベクトル、太字のＨｅxｔ’は外部からの印加磁場ベクトル、ΔＳ’は粒子の位置ベクトルを重心とする立方体要素の要素境界面の面積、ΔＶ’は粒子の位置ベクトルを重心とする立方体要素境界面の頂点と粒子位置を頂点とする立方体要素体積、χ’は磁気感受率、μ₀’は真空の透磁率、λ’はラグランジュの未定定数、πは円周率、ｓ’は粒子の位置ベクトルを重心とする立方体要素の要素境界面数である。

また、各物理量の添え字ｉｐ’はｉ番目の粒子の位置ベクトルを重心とする立方体要素内でのｐ’点の物理量、添え字ｊｑ’はｊ番目の粒子の位置ベクトルを重心とする立方体要素内でのｑ’点の物理量を示す。

式（１３）において、ｍ’は粒子の質量、ｖ’は速度、φ’（ｒ_ｉ’−ｒ_ｊ’）はｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギである。添え字ｉ、ｊはそれぞれ、ｉ、ｊ番目の粒子の物理量を示す。

次に、正準変数を（太字の）Ｈ_ｉｐ’’、（太字の傍点付き）Ｈ_ｉｐ’’とし、式（１３）〜式（１５）で示されるラグランジアンをラグランジュの運動方程式に代入すると、繰り込まれた粒子系Ｓ’における磁場の運動方程式は式（１６）のように記載できる。

ここで、式（１６）の右辺第２項は、ラグランジュの未定定数を通して拘束条件（磁化ベクトルの発散は０）を課している。
式（１６）の右辺第３項は外部からの印加磁場ベクトルが変化したときにすばやく追従させるための減衰項であり、γ’は減衰定数である。

式（１６）の右辺第２項に示される拘束条件を含んだ運動方程式を解くにあたり、本実施形態では、一般化された拘束条件の導入法であるＳＨＡＫＥ法を採用する。
拘束条件を考慮せずに蛙跳び法により式（１６）を離散化すると以下の式（１７）、式（１８）、（１９）になる。

ここで、δｔ’は磁化現象の収束計算を行う上で用いる仮想時間刻みである。
添え字のｎは任意の整数であり、ｎδｔ’における物理量、ｎ−１／２は（ｎ−１／２）δｔ’における物理量、ｎ＋１／２は（ｎ＋１／２）δｔ’における物理量、ｎ＋１は（ｎ＋１）δｔ’における物理量に対応している。

磁場演算部６６は、繰り込まれた粒子系Ｓ’における粒子の位置ベクトルを重心とする立方体要素内のｐ’点、ｑ’点の位置ベクトルを読み込む。

次に、磁場演算部６６は、ステップＳ１０６で既に計算され記憶装置５に記憶されているコイル４５および永久磁石３９が磁性体に対応する粒子を重心とする立方体要素内部のｐ’点に作る磁場ベクトルを外部からの印加磁場ベクトルとして読み込む。
さらに、磁場演算部６６は、既に計算され記憶装置５に記憶されている、磁性体に対応する粒子の要素境界面積、法線ベクトルおよび粒子の位置ベクトルを重心とする立方体要素境界面の頂点と粒子位置を頂点とする立方体要素体積を読み込む。
また、磁場演算部６６はあらかじめ記憶装置５に記憶されている減衰定数、仮想質量、仮想時間刻みを読み込む。

次に、磁場演算部６６はあらかじめ記憶装置５に記憶されている磁場ベクトル、磁場ベクトルの時間微分の初期値を読み込む。
なお、磁場演算部６６は後述するステップＳ１１０において磁場ベクトル、磁場ベクトルの時間微分が更新されている場合は、その値を読み込む。
なお、磁化ベクトルは、初期条件に含まれる磁化曲線を表す関数と磁場ベクトルとを使用して計算される。

磁場演算部６６は、式（１７）、式（１８）、式（１９）を、磁性体に対応する粒子の位置ベクトルを重心とする立方体要素内の全てのｐ’点に対して計算し、計算された磁場ベクトルを記憶装置５に記憶する。

（ステップＳ１１０）
磁場演算部６６は、ステップＳ１０８で計算し記憶装置５に記憶されている、磁性体に対応する粒子の位置ベクトルを重心とする立方体要素内のｐ’点の磁場ベクトルに拘束条件を加え、計算された磁場ベクトルを、記憶装置５に記憶する。

磁性体に対応する粒子の位置ベクトルを重心とする立方体要素内のｐ’点の磁場ベクトルに以下に示す式（２０）、式（２１）に従って拘束条件を加える。

ここで、太字のＨ’は磁場ベクトル、α’は仮想質量、δｔ’は仮想時間刻み、γ’は減衰定数、ｓ’は磁性体に対応する粒子の位置ベクトルを重心とする立方体要素の面数であり、添え字のｉは磁性体に対応する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｉｐ’は磁性体に対応する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする立方体要素内のｐ’点の物理量、ｉｑ’は磁性体に対応する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする立方体要素内のｑ’点の物理量、添え字のｎは任意の整数であり、ｎδｔ’における物理量、ｎ＋１は（ｎ＋１）δｔ’おける物理量に対応している。

磁場演算部６６は、ステップＳ１０８で計算され記憶装置５に記憶されている磁性体に対応する粒子を重心とする立方体要素内のｐ’点の磁場ベクトルを読み込む。
磁場演算部６６は、あらかじめ記憶装置５に記憶されている仮想質量、仮想時間刻み、減衰定数を読み込む。
磁場演算部６６は、磁性体に対応する粒子各々に対して式（２０）、式（２１）に基づく計算を行い、計算された磁場ベクトルを記憶装置５に記憶する。

（ステップＳ１１２）
次に、磁場演算部６６はステップＳ１１０で求めた磁場ベクトルが拘束条件を満たしているかを判断し、満たしていれば次のステップに進み、満たしていなければステップＳ１１０に戻る。

具体的には磁場演算部６６は、ステップＳ１１０で計算し記憶装置５に記憶されている磁性体に対応する粒子各々の磁場ベクトルから計算される磁化ベクトルを用いて式（２２）に基づく計算を行う。

ここで、添え字のｎ＋１は（ｎ＋１）δｔ’おける物理量に対応している。
ｅｒｒ_ｉ’は磁性体に対応する粒子の内、ｉ番目の粒子の拘束条件に対する誤差値である。

磁化ベクトル（太字の）Ｍ’は、ステップＳ１１０において計算され記憶装置５に記憶されている磁場ベクトルを磁場演算部６６が読み込み、これと初期条件に含まれる磁化曲線を表す関数とを使用することで計算される。

法線ベクトル（太字の）ｎ’には、既に計算され記憶装置５に記憶されているものを磁場演算部６６は読み込み代入する。

磁場演算部６６は、すべての粒子に対して、誤差の値が式（２３）を満たさなければステップＳ１１０に戻る。

式（２３）においてＡ_１’は任意の誤差判別値であり、あらかじめ記憶装置５に任意の値が記憶されている。

（ステップＳ１１４）
次に、磁場演算部６６は、磁性体の磁化現象が定常状態に到達したかを判断し、条件を満たしていれば次のステップに進む。
磁性体に対応する粒子が定常状態に到達したかは、以下の式（２４）により判断される。

ここで、太字のＨ’は磁場ベクトル、μ₀’は真空の透磁率、Ｎ’−ｗは磁性体に対応する粒子の粒子数、ｓ’は磁性体に対応する粒子の位置ベクトルを重心とする立方体要素の面数であり、添え字のｉｐ’は磁性体を構成する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする立方体要素内のｐ’点での物理量である。
さらに、添え字のｎは任意の整数であり、ｎδｔ’における物理量、ｎ―１は（ｎ―１）δｔ’、ｎ＋１は（ｎ＋１）δｔ’おける物理量に対応している。
また、Ａ_２’は磁性体の磁場ベクトルが定常状態に到達したかを判断するための任意の誤差判定値である。

磁場演算部６６は、あらかじめ記憶装置５に記憶されているＡ_２’、μ₀’を読み込む。
また、磁場演算部６６は、ステップＳ１１０で計算され記憶装置５に記憶されている、磁性体に対応する粒子を重心とする立方体要素内のｐ’点の磁場ベクトルを読み込む。

次に、磁場演算部６６は式（２４）を計算し、式（２４）を満たしていれば次のステップに進み、満たしていなければステップＳ１０８に戻る。
ステップＳ２１０は、ステップＳ１０２，ステップＳ１０４、ステップＳ１０６、ステップＳ１０８、ステップＳ１１０、ステップＳ１１２、ステップＳ１１４、を含む。

（ステップＳ１１６）
機構・弾性演算部６８は、ＳＰＭモータ３１に対応するすべての粒子の位置ベクトル上の磁場ベクトル、磁化ベクトル、磁束密度ベクトルを計算し、記憶装置５に記憶する。すなわち、機構・弾性演算部６８は、繰り込まれた粒子系Ｓ’において、磁性体に対応するすべての粒子の位置ベクトル上の磁場ベクトル、磁化ベクトル、磁束密度ベクトルと、コイル４５に対応するすべての粒子の位置ベクトル上の磁場ベクトル、磁束密度ベクトルを計算し、記憶装置５に記憶する。

（磁性体に対応する粒子）
磁性体に対応する粒子の位置ベクトル上の磁場ベクトルは以下の式（２５）で表される。

ここで、（太字の）ｅ_x’’＝（１,０,０）、（太字の）ｅ_ｙ’’＝（０,１,０）、（太字の）ｅ_ｚ’’＝（０,０,１）である。
太字のＨ’は磁場ベクトル、太字のｎ’は法線ベクトルであり、添え字ｉは磁性体に対応する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を表し、添え字ｉｐ’は磁性体に対応する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする立方体要素内のｐ’点の物理量である。
磁性体に対応する粒子の磁化ベクトルは、磁場ベクトルと磁化曲線を表す関数とを使用して求めることができる。

磁性体に対応する粒子の位置ベクトル上での磁束密度ベクトルは、式（２６）で表される。

ここで、太字のＢ’は磁束密度ベクトル、太字のＨ’は磁場ベクトル、太字のＭ’は磁化ベクトルであり、添え字のｉは磁性体に対応する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量を指し示す。μ₀’は真空の透磁率である。

機構・弾性演算部６８は、既に計算されている、磁性体に対応する粒子を重心とする立方体要素内のｐ’点での磁場ベクトルを記憶装置５から読み込む。

次に、機構・弾性演算部６８は、既に計算され記憶装置５に記憶されている法線ベクトルを読み込み、式（２５）に基づいて磁性体に対応する粒子の位置ベクトル上での磁場ベクトルを計算し記憶装置５に記憶する。

さらに、機構・弾性演算部６８は、このステップで計算された磁場ベクトルと初期条件に含まれる磁化曲線を表す関数とを使用して、磁性体に対応する粒子の磁化ベクトルを計算し記憶装置５に記憶する。

次に機構・弾性演算部６８は、このステップで計算された磁場ベクトルと磁化ベクトルを式（２６）に代入し、磁性体に対応する粒子の位置ベクトル上の磁束密度ベクトルを計算し記憶装置５に記憶する。
なお、真空の透磁率はあらかじめ記憶装置５に記憶されているものを制御部３が読み込む。

（コイルに対応する粒子）
コイル４５に対応する粒子の位置ベクトル上の磁場ベクトルは式（２７）、式（２８）で記載される。

ここで、太字のｒ’は位置ベクトル、太字のＨ’は磁場ベクトル、Ｎ’−ｗは磁性体に対応する粒子の数、ｓ’およびΔＳ’は磁性体に対応する粒子を重心とする立方体要素の境界面数と境界面積、太字のＭ’は磁化ベクトル、太字のｎ’は法線ベクトルであり、添え字のｉはコイル４５に対応する粒子の内、ｉ番目の粒子の位置ベクトル上での物理量、ｊは磁性体に対応する粒子の内、ｊ番目の粒子の位置ベクトル上での物理量であり、ｊｑ’は磁性体に対応する粒子の内、ｊ番目の粒子の位置ベクトルを重心とする立方体要素のｑ’点での物理量を指し示す。

式（２７）の右辺第一項は、磁性体に対応する粒子がコイル４５に対応する粒子の位置ベクトル上に作る磁場ベクトルを記述する項であり、式（２８）で表される。
式（２７）の右辺第二項は、コイル４５に対応する粒子がコイル４５に対応する粒子の位置ベクトル上に作る磁場ベクトルであり式（１）〜（１２）で表される。

機構・弾性演算部６８は、繰り込まれた粒子系Ｓ’におけるコイル４５に対応する粒子の位置ベクトルを記憶装置５より読み込む。
なお、ステップＳ１２０においてコイル４５に対応する粒子の位置ベクトルが更新されている場合は、機構・弾性演算部６８はその値を読み込む。

機構・弾性演算部６８は計算され記憶装置５に記憶されている磁性体に対応する粒子の数と立方体要素の面数、および粒子の位置ベクトルを重心とする立方体要素内のｑ’点の位置ベクトルと立方体要素境界面積を読み込む。
また、機構・弾性演算部６８はこのステップで計算され記憶装置５に記憶されている、磁性体に対応する粒子の磁化ベクトルを読み込む。
さらに、機構・弾性演算部６８は、繰り込まれた粒子系Ｓ’におけるコイル４５の寸法および電流密度ベクトルを読み込む。

その後、機構・弾性演算部６８は式（２７）、式（２８）および式（１）〜（１２）に従い、コイル４５に対応する粒子の位置ベクトル上での磁場ベクトルを計算し、記憶装置５に記憶する。

機構・弾性演算部６８は、コイル４５に対応する粒子の位置ベクトル上での磁束密度ベクトル（太字の）Ｂ_ｉ’を、すでに計算されているコイル４５に対応する粒子の位置ベクトル上の磁場ベクトルに真空の透磁率を掛けて計算し記憶装置５に記憶する。なお、真空の透磁率は、あらかじめ記憶装置５に記憶されている。

このように、制御部３は、磁場の運動方程式を解くことにより磁性体の磁化現象を解析する。そのため、有限要素法のように空間全域をメッシュ分割する必要がなく、磁性体を有する系に対して十分な精度で効率よく磁場解析を行うことができる。
また、運動方程式を解くため、行列を扱わず、計算に必要なメモリ量は粒子数に比例する。

（ステップＳ１１８）
機構・弾性演算部６８はＳＰＭモータ３１に対応する全ての粒子に働く力ベクトルを計算し、記憶装置５に記憶する。

（磁性体に対応する粒子に働く力ベクトル）
正準変数を太字のｒ_ｉ’、太字の傍点付きｒ_ｉ’とし、式（１３）のラグランジアンをラグランジュの運動方程式に代入する。
すると、磁性体に対応する粒子に働く力ベクトルは式（２９）のように記載される。

ここで、μ₀’は真空の透磁率、太字のＨ’は磁場ベクトル、太字のＭ’は磁化ベクトル、太字のｒ’は位置ベクトル、ΔＶ’は粒子の位置ベクトルを重心とする立方体要素境界面の頂点と粒子位置を頂点とする立方体要素体積、Ｎ’−ｗは磁性体に対応する粒子数であり、添え字のｉ、ｊは磁性体に対応する粒子の内、ｉ番目、ｊ番目の粒子の物理量を指し示し、添え字のｉｐ’は磁性体に対応する粒子の内、ｉ番目の粒子の位置ベクトルを重心とする立方体要素内のｐ’点の物理量を指し示す。
また、φ’（ｒ_ｉ’−ｒ_ｊ’）はｉ番目の粒子とｊ番目の粒子の相互作用ポテンシャルエネルギであり、ステップＳ２０８において求められている。

式（２９）の右辺第１項は、相互作用演算部６４によって演算された相互作用ポテンシャルエネルギφ’から導かれる粒子に働く力である。特に相互作用ポテンシャルエネルギφ’が弾性を考慮して決定されている場合、この力は弾性力に相当する。式（２９）の右辺第２項は、磁場演算部６６によって演算された磁場から導かれる粒子に働く磁気力である。式（２９）は、磁性体に対応する粒子には弾性力と磁気力とを足し合わせた力が働くことを示す。

（コイルに対応する粒子に働く力ベクトル）
コイル４５に対応する粒子のうち、ｉ番目の粒子の位置ベクトルを太字のｒ_ｉ’、その位置ベクトルでの電流の単位ベクトルを太字のｔ_ｉ’、磁束密度ベクトルを太字のＢ_ｉ’、電流の流れる方向のコイルの長さをＬ_ｉ’とすると、コイル４５に対応する粒子の内、ｉ番目の粒子に働く力ベクトルは以下の式（３０）で記載される。

（ステップＳ１２０）
機構・弾性演算部６８は、繰り込まれた粒子系Ｓ’における粒子の運動方程式を数値的に解くことにより、各粒子の位置ベクトル、速度ベクトルを更新する。この運動方程式で使用する力ベクトルとしては、ステップＳ１１８で求められた力ベクトルを使用する。
ステップＳ２１２は、ステップＳ１１６、ステップＳ１１８、ステップＳ１２０、を含む。

本実施の形態に係る解析装置１によると、相互作用演算部６４における弾性解析、磁場演算部６６における磁場解析、機構・弾性演算部６８における機構解析、制御パラメータ更新部７０における制御解析、が連成されている。したがって、解析対象の物体に対する高精度なシミュレーションが可能となる。

特に、磁場・機構・制御の連成解析に繰り込み群分子動力学の手法を導入して弾性をさらに連成させたので、シミュレーションにおいて解析対象の弾性的な振る舞いに基づく磁場の変化なども考慮に入れることができる。

また、本実施の形態に係る解析装置１では、分子動力学法計算により、接触問題なども容易に解くことが可能である。さらに、１つのラグランジアンを基に計算するため、系のエネルギを保存した形での計算が可能となる。また、磁場の運動方程式による磁場計算を行うことで空間メッシュは不要となる。

また、本実施の形態に係る解析装置１では、磁気エネルギに起因する変形、それに伴う磁気回路の変化、磁気回路変化による変形量の変化といった一連の現象を高精度に計算することが可能となり、実際の現象により近いシミュレーションが可能となる。

以上、実施の形態に係る解析装置１の構成と動作について説明した。この実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。

上記では、磁性体に対応する粒子のラグランジアンから導出される磁場の運動方程式を解くことで磁性体の磁化現象を計算した場合について説明したが、本発明は、何等、これに限定されることなく、運動方程式であればラグランジアンから導出されたもの以外のものであってもよい。

また、上記では本実施の形態をＳＰＭモータ３１の解析に使用したが、本発明は、何等、これに限定されることなく、磁性体が配置された空間における任意の点の磁場を解析するものであれば、例えばリニアモータ等のＳＰＭモータ３１以外のモータ、あるいは磁性体で空間を囲む磁気シールドにおいて、磁性体でシールドされた空間内の解析に用いても良い。

以下、粒子系Ｓのパラメータと繰り込まれた粒子系Ｓ’のパラメータとの間の関係について、本発明者が独自に行った考察を説明する。
第１の繰り込み因子をα、第２の繰り込み因子をγ、第3の繰り込み因子をδ、空間の次元数をdとした場合、特許文献２によると、
となる。

（磁場の繰り込み）
N個の原子上にスピン（核磁気モーメント）μ_ｉ（ｉ＝１、２、…、Ｎ）があるバルクを考える。バルク内部の核スピンμ_ｉが遠方に作る磁気誘導は、
と表される。これを粗視化すると、
となる。

以下に示されるように磁気モーメントを繰り込む。
繰り込まれた磁束密度は、
と表される。

磁束密度と同様に磁場を繰り込むと、
であり、繰り込まれた磁化は、
である。粗視化された磁気モーメントｍ＝Σμ_ｉに働く力は、
であり、以下に示されるようにスケールされる。

粗視化された電荷の流れｊが作る磁束密度は、
である。粗視化された磁気モーメントｍ_ｉと電流密度ｊの作る磁束密度との相互作用は、
であるから、電流密度は以下に示されるようにスケールされる。

以上を纏めると、
となる。

１ 解析装置、 ３ 制御部、 ５ 記憶装置、 ６ メディア入出力部、 ７ 入力部、 ９ 表示部、 １１ プリンタポート、 １２ プリンタ、 １３ バス、 ６０ 粒子モデル生成部、 ６２ 繰り込み部、 ６４ 相互作用演算部、 ６６ 磁場演算部、 ６８ 機構・弾性演算部、 ７０ 制御パラメータ更新部、 ７２ 終了判定部。

Claims (6)

1. 物体を複数の粒子からなる粒子系として記述した上で当該物体の振る舞いを解析する解析装置であって、
   粒子間の相互作用による相互作用ポテンシャルエネルギを演算する相互作用演算部と、
   各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより各粒子に作用する磁場を演算する磁場演算部と、
   前記相互作用演算部によって演算された相互作用ポテンシャルエネルギと前記磁場演算部によって演算された磁場とに基づき各粒子の運動を演算し、各粒子の位置を更新する機構・弾性演算部と、を備えることを特徴とする解析装置。
2. 前記機構・弾性演算部は、前記相互作用演算部によって演算された相互作用ポテンシャルエネルギから導かれる粒子に働く第１の力と、前記磁場演算部によって演算された磁場から導かれる当該粒子に働く第２の力と、を足し合わせて当該粒子に働く合計の力とすることを特徴とする請求項１に記載の解析装置。
3. 前記粒子系の外部から与えられる制御パラメータであって各粒子に作用する磁場に影響を及ぼす制御パラメータを更新する制御更新部をさらに備えることを特徴とする請求項１または２に記載の解析装置。
4. 前記粒子系は繰り込まれた粒子系であることを特徴とする請求項１から３のいずれかに記載の解析装置。
5. 物体を複数の粒子からなる粒子系として記述した上で当該物体の振る舞いを解析する解析方法であって、
   粒子間の相互作用による相互作用ポテンシャルエネルギを演算するステップと、
   各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより各粒子に作用する磁場を演算するステップと、
   演算された相互作用ポテンシャルエネルギと演算された磁場とに基づき各粒子の運動を演算し、各粒子の位置を更新するステップと、を含むことを特徴とする解析方法。
6. 物体を複数の粒子からなる粒子系として記述した上で当該物体の振る舞いを解析する機能をコンピュータに実現させるコンピュータプログラムであって、
   粒子間の相互作用による相互作用ポテンシャルエネルギを演算する機能と、
   各粒子に作用する磁場が満たすべき関係を粒子ごとに運動方程式の形式で記述した磁場の運動方程式を数値的に解くことにより各粒子に作用する磁場を演算する機能と、
   演算された相互作用ポテンシャルエネルギと演算された磁場とに基づき各粒子の運動を演算し、各粒子の位置を更新する機能と、を前記コンピュータに実現させることを特徴とするコンピュータプログラム。