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