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