JP7264841B2 - 磁場解析方法、磁場解析装置、及びプログラム - Google Patents

Description

本発明は、磁場解析方法、磁場解析装置、及びプログラムに関する。

古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法に基づくシミュレーションが知られている（例えば、特許文献１参照）。分子動力学法のような粒子法は静的な現象だけでなく流れなどの動的な現象をも取り扱えるので、主に静的な現象を解析対象とする有限要素法などに代わるシミュレーション手法として注目されている。例えば、各粒子に磁気モーメントを付与し、各粒子の球対称性に基づく厳密解を利用して磁気的な物理量を演算することで、比較的精度の高いシミュレーション結果を高速に得ることのできる磁気ビーズ法が提案されている（例えば、特許文献２参照）。

上述の磁気ビーズ法では、対象物を複数個の要素に分割し、各要素を球形状の粒子とみなして球対称性に基づく厳密解を利用して磁場を演算している。このとき、対象物を複数の粒子で充填させようとすると粒子間に隙間が生じるため、その隙間や球形状の粒子に起因する誤差が生じ得る。計算精度を高めるために、対象物を立方体などの多面体要素で隙間なく分割し、かつ計算時間の短縮化を図ることが可能な手法が提案されている（例えば、特許文献３参照）。

特開２００９－３７３３４号公報 特開２０１５―１１１４０１号公報 特開２０１７―１９４８８４号公報

磁気ビーズ法を用いた磁場解析において、一般的に反復法が用いられる。磁気ビーズ法に用いられる反復法では、磁場解析モデルの各粒子に磁化を付与し、各粒子内の評価点における磁場を最適解に収束させる運動方程式を数値的に反復計算して解く。この反復計算の計算時間を短縮することが望まれている。

本発明の目的は、反復法を用いて磁場を求めるための反復計算の計算時間を短縮することが可能な磁場解析方法、磁場解析装置、及びプログラムを提供することである。

本発明の一観点によると、
複数の多面体要素に分割された解析モデルの、前記複数の多面体要素の各々の内部の評価点の磁場を反復法により計算する磁場解析方法であって、
反復計算の１つのタイムステップが、
前記複数の評価点の各々の現時点の磁場を表す第１磁場ベクトルに基づいて、１タイムステップ経過後の前記複数の評価点の各々の磁場を表す第２磁場ベクトルを求める手順と、
前記第１磁場ベクトルから前記第２磁場ベクトルまでの磁場の変化を表すベクトルである磁場変化ベクトルを求める手順と、
現時点のタイムステップで求めた前記磁場変化ベクトルと、直前のタイムステップで求めた前記磁場変化ベクトルとのなす角度が９０°未満のときに、現時点のタイムステップで求めた前記磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを、前記第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる前記第１磁場ベクトルとして設定する手順と
を含む磁場解析方法が提供される。

本発明の他の観点によると、
複数の多面体要素に分割された解析モデルの、前記複数の多面体要素の各々の内部の評価点の磁場を反復法により計算する磁場演算部を有する磁場解析装置であって、
前記磁場演算部は、反復計算の１つのタイムステップにおいて、
前記複数の評価点の各々の現時点の磁場を表す第１磁場ベクトルに基づいて、１タイムステップ経過後の前記複数の評価点の各々の磁場を表す第２磁場ベクトルを求める手順と、
前記第１磁場ベクトルから前記第２磁場ベクトルまでの磁場の変化を表すベクトルである磁場変化ベクトルを求める手順と、
現時点のタイムステップで求めた前記磁場変化ベクトルと、直前のタイムステップで求めた前記磁場変化ベクトルとのなす角度が９０°未満のときに、現時点のタイムステップで求めた前記磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを、前記第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる前記第１磁場ベクトルとして設定する手順と
を実行する磁場解析装置が提供される。

本発明のさらに他の観点によると、
複数の多面体要素に分割された解析モデルの、前記複数の多面体要素の各々の内部の評価点の磁場を反復法により計算する手順をコンピュータに実行させるプログラムであって、
反復計算の１つのタイムステップにおいて、
前記複数の評価点の各々の現時点の磁場を表す第１磁場ベクトルに基づいて、１タイムステップ経過後の前記複数の評価点の各々の磁場を表す第２磁場ベクトルを求める手順と、
前記第１磁場ベクトルから前記第２磁場ベクトルまでの磁場の変化を表すベクトルである磁場変化ベクトルを求める手順と、
現時点のタイムステップで求めた前記磁場変化ベクトルと、直前のタイムステップで求めた前記磁場変化ベクトルとのなす角度が９０°未満のときに、現時点のタイムステップで求めた前記磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを、前記第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる前記第１磁場ベクトルとして設定する手順と
をコンピュータに実行させるプログラムが提供される。

現時点のタイムステップで求めた磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを、第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる第１磁場ベクトルとして設定することにより、反復計算の計算時間を短縮することが可能である。

図１は、本実施例による磁場解析装置のブロック図である。 図２Ａは、仮想空間内に配置された解析対象である磁性体を複数のボクセルに分割した状態の模式図であり、図２Ｂは、１つのボクセル内に配置された複数の評価点を示す模式図である。 図３は、実施例による磁場解析方法を示すフローチャートである。 図４は、ｉ番目のボクセルとｊ番目のボクセルとの関係を示す模式図である。 図５は、反復計算の１タイムステップで求められる種々の磁場ベクトルを磁場ベクトル空間内に示した模式図である。 図６は、反復計算の１タイムステップで求められる種々の磁場ベクトルを磁場ベクトル空間内に示した模式図である。 図７は、解析モデルの斜視図である。 図８Ａ及び図８Ｂは、ある１つのボクセルの磁場の仮想時間軸における変化を示すグラフである。 図９は、図３のステップＳ１５でマルチタイムステップ法を適用しないで磁場解析を行った結果の暫定磁場ベクトルと第２磁場ベクトルとの仮想時間軸における変化を示すグラフである。 図１０は、実施例による方法、ＦＩＲＥ法を適用しマルチステップ法を適用しない方法、及びＦＩＲＥ法とマルチステップ法とのいずれも適用しない方法で、磁場解析を行った結果を示すグラフである。

図１は、本実施例による磁場解析装置３０のブロック図である。本実施例による磁場解析装置３０は、処理装置２０、記憶装置２５、入力装置２８、及び出力装置２９を含む。処理装置２０は、解析情報取得部２１、ボクセル分割部２２、磁場演算部２３、及び出力制御部２４を含む。

図１に示す各ブロックは、ハードウェア的には、コンピュータの中央処理ユニット（ＣＰＵ）をはじめとする素子や機械装置で実現することができ、ソフトウェア的にはコンピュータプログラム等によって実現することができる。図１では、ハードウェア及びソフトウェアの連携によって実現される機能ブロックが示されている。従って、これらの機能ブロックは、ハードウェア及びソフトウェアの組み合わせによって、種々の態様で実現することが可能である。

処理装置２０は入力装置２８及び出力装置２９と接続される。入力装置２８は、処理装置２０で実行される処理に関係するユーザからのコマンド及びデータの入力を受ける。入力装置２８として、例えばユーザが操作を行うことにより入力を行うキーボードやマウス、インターネット等のネットワークを介して入力を行う通信装置、ＣＤ、ＤＶＤ等の記録媒体から入力を行う読取装置等を用いることができる。出力装置２９は、処理装置２０からの制御により、解析結果を出力する。出力装置２９として、例えば画像を表示する表示装置、インターネット等のネットワークを介してデータを送信する通信装置、ＣＤ、ＤＶＤ等の記録媒体にデータを記録する書込装置等を用いることができる。

解析情報取得部２１は、入力装置２８を介して磁場解析情報を取得する。磁場解析情報には、磁場解析に必要な種々の情報が含まれる。例えば、仮想空間内に定義される解析対象物の形状や物性値を定義する情報、解析対象物をボクセルに分割するための情報、外部磁場を定義する情報、解析にあたり定義しておくべきパラメータの値等が含まれる。

図２Ａ及び図２Ｂを参照して、ボクセル分割部２２の機能について説明する。
図２Ａは、仮想空間内に配置された解析対象である磁性体を複数のボクセル４０に分割した状態の模式図であり、図２Ｂは、１つのボクセル４０内に配置された複数の評価点４１を示す模式図である。

ボクセル分割部２２は、図２Ａに示すように、磁場解析情報に基づき仮想空間内に配置された磁性体を複数のボクセル４０に分割する。図２Ａでは、ボクセル４０の分布を二次元的に表しているが、実際には磁性体は三次元の形状を有し、複数のボクセル４０は三次元空間内に配置される。また、図２Ａでは、１つのボクセル４０を正方形で表しているが、三次元に拡張するとボクセル４０は正六面体になる。なお、ボクセル４０を正六面体以外の多面体としてもよい。例えば、磁性体をボロノイ分割することによりボクセル４０を配置してもよい。ボクセル分割部２２は、ボクセル４０の各々に、１つのボクセルを特定するためのボクセル識別子（ボクセルＩＤ）を付与する。ボクセル識別子として、１から始まる通し番号ｉを用いる。

ボクセル分割部２２は、さらに、図２Ｂに示すようにボクセル４０の各々の内部に複数の評価点４１を配置する。評価点４１は、例えばボクセル４０の重心から各面に下した垂線の中点に配置される。１つのボクセル４０の内部に配置される評価点４１の個数は、多面体であるボクセル４０の表面を構成する平面の数に等しい。例えば、正六面体のボクセル４０の内部には、６個の評価点４１が配置される。１つのボクセル４０内に複数の評価点４１を配置するのは、磁場解析精度を高く保つためである。

ボクセル分割部２２は、ボクセル４０の各々の内部に配置された複数の評価点４１に、評価点識別子（評価点ＩＤ）を付与する。評価点識別子として、１から始まる通し番号ｐを用いる。複数の評価点４１のうち１つの評価点４１は、ボクセル識別子と評価点識別子との対によって特定される。

磁場演算部２３は、複数の評価点４１の各々に磁化を付与し、各評価点４１に発生する磁場を計算する。図２Ｂに示すようにｉ番目のボクセル４０（ｉ）の内部のｐ番目の評価点４１（ｉ，ｐ）に付与された磁化ベクトルをＭ_ｉｐと表記する。また、評価点４１（ｉ，ｐ）に対応する表面の外向きの単位法線ベクトルをｎ_ｉｐと表記する。磁場演算部２３が行う手順については、後に図３～図６を参照して説明する。

出力制御部２４は、解析結果、例えば記憶装置２５に記憶されている評価点４１ごとの磁場の計算結果を出力装置２９に出力する。

図３は、実施例による磁場解析方法を示すフローチャートである。
まず、解析情報取得部２１（図１）が解析条件を取得し、ボクセル分割部２２（図１）が解析対象の磁性体を複数のボクセル４０（図２Ａ）に分割する（ステップＳ１１）。さらに、ボクセル分割部２２は、複数のボクセル４０の各々の内部に複数の評価点４１（図２Ｂ）を配置する（ステップＳ１１）。

次に、磁場演算部２３（図１）が、評価点４１の各々について、ステップＳ１２からステップＳ１６までの１タイムステップ分の磁場の演算を、すべての評価点４１について演算が終了するまで繰り返す（ステップＳ１７、Ｓ１８）。すべての評価点４１についての１タイムステップ分の磁場の演算を、磁場の計算結果が収束するまで繰り返す（ステップＳ１９、Ｓ２０）。磁場の計算結果が収束すると、解析を終了し、解析結果を出力装置２９（図１）に出力する（ステップＳ２１）。

ステップＳ１２からステップＳ１６までの処理について説明する前に、反復法によって磁場を計算する方法について説明する。

反復法では、磁場を最適解に収束させる運動方程式を数値的に解く。この運動方程式は、以下の式で表される。

ここで、Ｈ_ｉｐは、ｉ番目のボクセル４０（ｉ）内のｐ番目の評価点４１（ｉ，ｐ）（図２Ｂ）に付与されている磁場ベクトルである。α_ｉｐは、評価点４１（ｉ，ｐ）の磁場を時間発展させるときに用いる仮想質量である。λ_ｉはラグランジュの未定乗数である。γは人工的に付与された減衰定数である。

磁場ベクトルＨ（ｒ_ｉｐ）は、磁性体内の磁化、及び外部磁場によって評価点４１（ｉ，ｐ）の位置に発生する磁場である。図４を参照して、式（１）の右辺の磁場ベクトルＨ（ｒ_ｉｐ）について説明する。

図４は、ｉ番目のボクセル４０（ｉ）とｊ番目のボクセル４０（ｊ）との関係を示す模式図である。ｊ番目のボクセル４０（ｊ）のｑ番目の評価点４１（ｊ、ｑ）に磁化Ｍ_ｊｑが付与されている。磁性体の磁化率をχと表記すると、磁化Ｍ_ｊｑは磁場Ｈ_ｊｑにより以下の式で表される。

評価点４１（ｉ，ｐ）の位置ベクトルをｒ_ｉｐと表記し、評価点４１（ｊ、ｑ）の位置ベクトルをｒ_ｊｑと表記する。評価点４１（ｊ、ｑ）の磁化Ｍ_ｊｑによって評価点４１（ｉ，ｐ）の位置に磁場が発生する。磁性体内のすべての評価点４１に付与された磁化によって評価点４１（ｉ，ｐ）に発生する磁場Ｈ（ｒ_ｉｐ）は、以下の式で表される。

ここで、ΔＳ_ｊｑ及びｎ_ｊｑは、それぞれボクセル４０（ｊ）の評価点４１（ｊ，ｑ）に対応する表面の面積及び外側を向く単位法線ベクトルである。Ｈ_{ｅｘｔ，ｉｐ}は、評価点４１（ｉ，ｐ）の位置における外部磁場を表す磁場ベクトルである。ｑに関するシグマは、ｊ番目のボクセル４０（ｊ）内のすべての評価点４１（ｊ、ｑ）についての和をとることを意味する。ｊに関するシグマは、すべてのボクセル４０（ｊ）についての和をとることを意味する。

式（１）の右辺のＣ_ｉ（Ｍ_ｐ）は、ｉ番目のボクセル４０（ｉ）の内部における磁化ベクトルの発散の体積積分を、ガウスの発散定理により面積分に置き換えたものであり、以下の式で表される。

ここで、式（４）の右辺のｐに関するシグマは、ｉ番目のボクセル４０（ｉ）内のすべての評価点４１について和をとることを意味する。

次に、式（１）の物理的な意味について説明する。
式（１）の左辺は、磁場Ｈの変化の加速度を表す加速度項に相当する。右辺の第１項はバネ力に相当し、第２項は外力に相当し、第４項は減衰項に相当する。右辺第３項は、運動方程式の束縛条件を表す束縛項である。式（１）に示した運動方程式には、その右辺第３項によって、ボクセル４０の各々の内部において磁化の発散がゼロであるという束縛条件が課されている。磁場の計算結果が収束した状態では、加速度項、減衰項、及び束縛項がゼロになる。この状態で、評価点４１（ｉ，ｐ）の磁場Ｈ_ｉｐは、磁性体内の評価点に付与された磁化により評価点４１（ｉ，ｐ）の位置に発生する磁場と等しくなる。

反復法においては、式（１）の運動方程式を、ある仮想的な時間刻み幅で数値的に反復計算し、評価点４１ごとに磁場Ｈ_ｉｐの時間変化を求める。このとき、磁場Ｈ_ｉｐの計算結果は、時間の経過に従って、その最適値に収束する。束縛条件付きの運動方程式を解く手法として、ＳＨＡＫＥ法（Ｊ．Ｍ．ティッセン「計算物理学」）が知られている。ＳＨＡＫＥ法では、まず束縛条件を課さずに運動方程式を解き、その後、束縛条件を付加する。以下に説明する実施例による方法では、ＳＨＡＫＥ法を適用する。

次に、図３に示したステップＳ１２からステップＳ１６までの手順について説明する。この手順の説明において、必要に応じて図５及び図６を参照する。
図５及び図６は、反復計算の１タイムステップで求められる種々の磁場ベクトルを磁場ベクトル空間内に示した模式図である。

まず、式（１）の運動方程式の束縛条件を課さず、ｎ番目のタイムステップが終了した時点（現時点）の評価点４１（ｉ，ｐ）における第１磁場ベクトルＨ１_ｉｐ（ｎ）、及び現時点で各評価点４１に付与されている磁化ベクトルＭ_ｊｑに基づいて、評価点４１（ｉ，ｐ）に発生する磁場を１タイムステップ分だけ時間発展させた磁場ベクトルである暫定磁場ベクトルＨＴ_ｉｐ（ｎ＋１）を計算する。以下、暫定磁場ベクトルＨＴ_ｉｐ（ｎ＋１）の計算方法について説明する。

束縛条件を課さない運動方程式（１）を、蛙跳び法で離散化すると、以下の式が得られる。

ここで、δｔは仮想時間の刻み幅であり、Ｆ_ｉｐ（ｎ）は以下の式で定義される。

式（５）～（７）を用いて、第１磁場ベクトルＨ１_ｉｐ（ｎ）を１タイムステップ分だけ時間発展させた暫定磁場ベクトルＨＴ_ｉｐ（ｎ＋１）を求めることができる。

暫定磁場ベクトルＨＴ_ｉｐ（ｎ＋１）を求めた後、運動方程式（１）の束縛条件が満たされるまで暫定磁場ベクトルＨＴ_ｉｐ（ｎ＋１）を更新して束縛条件を満たす磁場（第２磁場ベクトルＨ２_ｉｐ（ｎ＋１））を求める。

第２磁場ベクトルＨ２_ｉｐ（ｎ＋１）は、以下の式を、計算結果が収束するまで反復計算することにより求めることができる。

ここで、式（８）の右辺の暫定磁場ベクトルＨＴ_ｉｐ（ｎ＋１）は反復計算による更新前のベクトルを表し、左辺の暫定磁場ベクトルＨＴ_ｉｐ（ｎ＋１）は反復計算による更新後のベクトルを表す。式（９）の右辺の磁化ベクトルＭＴ_ｑ（ｎ＋１）は、暫定磁場ベクトルに磁化率χを乗じたものである。

式（８）を、暫定磁場ベクトルＨＴ_ｉｐが収束するまで反復計算して得られた暫定磁場ベクトルＨＴ_ｉｐの収束値が、第２磁場ベクトルＨ２_ｉｐ（ｎ＋１）に相当する。

第２磁場ベクトルＨ２_ｉｐ（ｎ＋１）が求まると、磁場演算部２３（図１）は、第２磁場ベクトルＨ２_ｉｐ（ｎ＋１）から第１磁場ベクトルＨ１_ｉｐ（ｎ）までの磁場の変化の速度を計算する（ステップＳ１４）。磁場の変化の速度は、以下の式で表される。

続いて、マルチタイムステップ法を適用して、次のタイムステップにおける第１磁場ベクトルＨ１_ｉｐ（ｎ＋１）を計算する（ステップＳ１５）。図５及び図６を参照して、ステップＳ１５における計算について説明する。

まず、直近に求めた磁場の変化の速度ベクトルと、前回のタイムステップで求めた磁場の変化の速度ベクトルとの内積Ｑを計算する。内積Ｑは以下の式で表される。

２つの速度ベクトルのなす角度が９０°未満の場合、内積Ｑは正の値をとる。図５は、内積Ｑが正の場合を示している。２つの速度ベクトルのなす角度が９０°より大きい場合、内積Ｑは負の値をとる。図６は、内積Ｑが負の場合を示している。２つの速度ベクトルのなす角度が９０°に等しい場合、内積Ｑはゼロになる。

磁場変化の速度ベクトルに、仮想時間の刻み幅δｔを乗じたベクトルは、第１磁場ベクトルＨ１_ｉｐ（ｎ）から第２磁場ベクトルＨ２_ｉｐ（ｎ＋１）までのベクトルの変化量（以下、磁場変化ベクトルという。）に相当する。式（１１）に代えて、現時点の磁場変化ベクトルと、直前のタイムステップで求めた磁場変化ベクトルとの内積を計算してもよい。磁場変化の２つの速度ベクトルのなす角度は、２つの磁場変化ベクトルのなす角度と等しい。

以下の式に基づいて、次のタイムステップで用いる第１磁場ベクトルＨ１_ｉｐ（ｎ＋１）を計算する。

ここで、ｆはマルチステップ法で適用される伸長短縮係数である。

図５に示すように、現時点の磁場変化の速度ベクトルと直前のタイムステップで求めた磁場変化の速度ベクトルとのなす角度が９０°未満のとき、伸長短縮係数ｆを１より大きくする。このとき、第１磁場ベクトルＨ１_ｉｐ（ｎ）に、現時点のタイムステップで求めた磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを加えたベクトルを、次のタイムステップにおいて用いる第１磁場ベクトルＨ１_ｉｐ（ｎ＋１）として設定する。

図６に示すように、現時点の磁場変化の速度ベクトルと直前のタイムステップで求めた磁場変化の速度ベクトルとのなす角度が９０°より大きいとき、伸長短縮係数ｆを１より小さくする。このとき、第１磁場ベクトルＨ１_ｉｐ（ｎ）に、現時点のタイムステップで求めた磁場変化ベクトルの向きを変えず、大きさを小さくしたベクトルを加えたベクトルを、次のタイムステップにおいて用いる第１磁場ベクトルＨ１_ｉｐ（ｎ＋１）として設定する。

次のタイムステップで用いる第１磁場ベクトルＨ１_ｉｐ（ｎ＋１）が求まったら、計算時間の短縮化を図るためのＦＩＲＥ法(fast inertial relaxation engine)を適用する（ステップＳ１６）。ＦＩＲＥ法に関しては、Erik Bitzek, Pekka Koskinen, Franz Gahler, Michael Moseler, and Peter Gumbsch, “Structural Relaxation Made Simple”, Physical Review Letters, 97, 170201 (2006)に説明されている。

次に、ＦＩＲＥ法について簡単に説明する。磁場の変化の速度ベクトルと加速度ベクトルとの内積が正またはゼロのとき、仮想質量α_ｉｐを現時点の値より小さくする。逆に、磁場の速度ベクトルと加速度ベクトルとの内積が負のとき、仮想質量α_ｉｐを現時点の値より大きくする。以下、仮想質量を変化させることの物理的な意味について説明する。

磁場の速度ベクトルと加速度ベクトルとの内積が正またはゼロの場合は、速度ベクトルと加速度ベクトルとのなす角度が９０°以下である。このとき、磁場の計算結果は収束値に近づいていると考えられる。この場合、仮想質量α_ｉｐを小さくすると、磁場の変化量が大きくなり、磁場の計算結果がより速く収束値に近づく。磁場の速度ベクトルと加速度ベクトルとの内積が負の場合は、速度ベクトルと加速度ベクトルとのなす角度が９０°より大きい。このとき、磁場の計算結果は収束値から遠ざかっていると考えられる。この場合、仮想質量α_ｉｐを大きくすると、磁場の変化量が小さくなり、磁場の計算結果が収束値から遠ざかる速さが遅くなる。

次に、上記実施例の優れた効果について説明する。
上記実施例では、図３に示したステップＳ１６でＦＩＲＥ法を適用しているため、磁場の計算結果の収束を速めることができる。その結果、計算時間を短縮することができる。

図５に示したように、現時点のタイムステップで求めた磁場変化の速度ベクトルと、直前のタイムステップで求めた磁場変化の速度ベクトルとのなす角度が９０°未満である場合、磁場の計算結果は収束値に近づいていると考えられる。このとき、式（１２）の伸長短縮係数ｆを１より大きくすることにより、磁場の計算結果をより速く収束値に近付けることができる。その結果、計算時間を短縮することができる。

図６に示したように、現時点のタイムステップで求めた磁場変化の速度ベクトルと、直前のタイムステップで求めた磁場変化の速度ベクトルとのなす角度が９０°より大きい場合、磁場の計算結果が収束値に近づいているか否か不明である。このとき、式（１２）の伸長短縮係数ｆを１未満とすることにより、磁場の計算結果が収束値から大きく離れてしまうことが回避される。その結果、計算時間を短縮することができる。

次に、図７～図１０を参照して、上記実施例の優れた効果を確認するために行った磁場解析結果について説明する。

図７は、磁場解析の対象となる解析モデルの斜視図である。解析モデルとなる磁性体の形状は、一辺の長さが１０ｍｍの立方体とした。この磁性体は、磁化が磁場に比例する線形磁性材料で形成されていることとした。解析モデルの１つの稜に平行で、空間的に均一な外部磁場Ｈ_ｅｘｔを解析モデルに与えた。外部磁場Ｅ_ｅｘｔの大きさは０．１Ｔとした。

まず、図３のステップＳ１５のマルチタイムステップ法を適用しないで磁場解析を行った。この磁場解析では、図５及び図６に示した伸長短縮係数ｆが常に１である。すなわち、第２磁場ベクトルＨ２_ｉｐ（ｎ＋１）が、次のタイムステップで使用される第１磁場ベクトルＨ１_ｉｐ（ｎ＋１）として設定される。

図８Ａ及び図８Ｂは、ある１つのボクセル４０（図２Ａ）の磁場の仮想時間軸における変化を示すグラフである。横軸はタイムステップ数を表し、縦軸は磁化の大きさを単位「Ａ／ｍ」で表す。図８Ｂは、図８Ａの横軸の原点近傍を拡大したものである。比較のために、図３のステップＳ１６においてＦＩＲＥ法を適用しないで解析を行った結果を破線で示す。

ＦＩＲＥ法を適用した場合には、ＦＩＲＥ法を適用しない場合に比べて、磁化の大きさがより速く収束値に近づいていることがわかる。この解析結果から、ＦＩＲＥ法を適用することにより計算時間を短縮することができるという優れた効果が得られることが確認された。

ただし、図８Ａに示したグラフから、解析初期段階における振動が収まった後、磁化の計算結果が収束するまで、磁化の計算結果の変化が緩やかであり、収束までに膨大なタイムステップが費やされていることが読み取れる。以下、図９を参照して、振動が収まった後に、磁化の計算結果が収束するまでに膨大なタイムステップが費やされる原因について説明する。

図９は、図３のステップＳ１５でマルチタイムステップ法を適用しないで磁場解析を行った結果の暫定磁場ベクトルＨＴ_ｉｐと第２磁場ベクトルＨ２_ｉｐとの仮想時間軸における変化を示すグラフである。マルチタイムステップ法を適用していないため、第２磁場ベクトルＨ２_ｉｐが、次のタイムステップにおける第１磁場ベクトルＨ１_ｉｐとして設定される。

磁化の大きさが急激に変化する解析の初期段階においては、外部磁場や他のボクセル４０からの寄与が支配的であるため、ステップＳ１２（図３）で求まる暫定磁場ベクトルＨＴ_ｉｐが第１磁場ベクトルＨ１_ｉｐから大きく変化する。この変化に対して、暫定磁場ベクトルＨＴ_ｉｐから第２磁場ベクトルＨ２_ｉｐまでの変化は相対的に小さい。このため、暫定磁場ベクトルＨＴ_ｉｐから第２磁場ベクトルＨ２_ｉｐまでの変化量が、磁場の収束の速さに与える影響は小さい。

これに対して、磁場の変化が緩やかになると、第１磁場ベクトルＨ１_ｉｐから暫定磁場ベクトルＨＴ_ｉｐまでの変化量と、暫定磁場ベクトルＨＴ_ｉｐから第２磁場ベクトルＨ２_ｉｐまでの変化量とがほぼ同等の大きさになる。このため、第１磁場ベクトルＨ１_ｉｐ（ｎ－１）から大きく変化した暫定磁場ベクトルＨＴ_ｉｐ（ｎ）は、運動方程式（１）の束縛条件を課して磁場を計算することにより、第２磁場ベクトルＨ２_ｉｐ（ｎ）が元の第１磁場ベクトルＨ１_ｉｐ（ｎ－１）の近傍まで戻される。その結果、１タイムステップでの磁場の変化量が極僅かになっている。

ステップＳ１６（図３）においてＦＩＲＥ法を適用することによって仮想質量を小さくすると、第１磁場ベクトルＨ１_ｉｐ（ｎ－１）から暫定磁場ベクトルＨＴ_ｉｐ（ｎ）までの磁場の変化量が大きくなる。この変化の方向は、磁場の変化が緩やかになった時間帯では、本来の変化の方向からずれている。このため、磁場の変化が緩やかになった時間帯では、ＦＩＲＥ法を適用する効果は低い。

上記実施例では、磁場の変化が緩やかになった時間帯で、図５に示すように伸長短縮係数ｆを１より大きくすることによって、磁場変化ベクトルを伸長させて、次のタイムステップで用いる第１磁場ベクトルＨ１_ｉｐ（ｎ＋１）を求めている。このため、磁場の計算結果をより速く収束値に近付けることができる。

図１０は、ＦＩＲＥ法とマルチタイムステップ法との両方を適用する実施例による方法、ＦＩＲＥ法を適用しマルチステップ法を適用しない方法、及びＦＩＲＥ法とマルチステップ法とのいずれも適用しない方法で、磁場解析を行った結果を示すグラフである。横軸はタイムステップ数を表し、縦軸は１つのボクセル４０の磁化の大きさを単位「Ａ／ｍ」で表す。現時点の磁場変化ベクトルと、直前のタイムステップにおける磁場変化ベクトルとのなす角度が９０°未満のとき、伸長短縮係数ｆを１．１とした。また、２つの磁場変化ベクトルのなす角度が９０°以上のとき、伸長短縮係数ｆを０．９とした。

ＦＩＲＥ法を適用することにより、磁場の計算結果がより速く収束値に近づいていることがわかる。さらに、マルチタイムステップ法を適用することにより、磁場の計算結果がより速く収束値に近づいていることがわかる。図１０に示した解析結果により、計算時間を短縮することができるという本実施例の優れた効果が確認された。

次に、伸長短縮係数ｆの値について説明する。伸長短縮係数ｆを大きくし過ぎると、計算が不安定になり、計算結果が発散しやすくなる。反復計算を安定させるために、伸長短縮係数ｆを２未満とすることが好ましい。現時点の磁化変化ベクトルと、直前のタイムステップにおける磁場変化ベクトルとのなす角度が９０°の場合、伸長短縮係数ｆは１または１未満にするとよい。

次に、上記実施例の変形例について説明する。
上記実施例ではステップＳ１６（図３）でＦＩＲＥ法を適用したが、必ずしもＦＩＲＥ法を適用する必要はない。マルチタイムステップ法のみを適用することによって、磁場の計算結果の変化が緩やかになった時間帯において、磁場の計算結果をより速く収束値に近付けることができる。

式（１２）の伸長短縮係数ｆの値は、ユーザが入力装置２８（図１）から入力するようにするとよい。現時点の磁場変化ベクトルと直前のタイムステップにおける磁場変化ベクトルとのなす角度が９０°以上の場合、伸長短縮係数ｆの値を１に設定してもよい。

上記実施例では、ボクセル４０（図２Ｂ）の重心から各表面に下した垂線の中点に評価点４１を配置している。評価点の位置は、ボクセル内のその他の位置としてもよい。例えば、ボクセル４０の重心から各表面に下した垂線の中点より、表面に近い位置に評価点４１を配置してもよい。

また、上記実施例では、解析対象である磁性体を複数の正六面体のボクセル４０に分割しているが、その他の形状の複数の多面体要素に分割してもよい。また、分割された複数の多面体要素の面数は、すべての多面体要素の間で同一でなくてもよい。評価点４１の個数は、その評価点４１が含まれる多面体要素の面数と同一でなくてもよく、面数より少なくてもよいし、面数より多くてもよい。

本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

２０ 処理装置
２１ 解析情報取得部
２２ ボクセル分割部
２３ 磁場演算部
２４ 出力制御部
２５ 記憶装置
２８ 入力装置
２９ 出力装置
３０ 磁場解析装置
４０ ボクセル
４１ 評価点

Claims (7)

1. 複数の多面体要素に分割された解析モデルの、前記複数の多面体要素の各々の内部の評価点の磁場を反復法により計算する磁場解析方法であって、
   反復計算の１つのタイムステップが、
   前記複数の評価点の各々の現時点の磁場を表す第１磁場ベクトルに基づいて、１タイムステップ経過後の前記複数の評価点の各々の磁場を表す第２磁場ベクトルを求める手順と、
   前記第１磁場ベクトルから前記第２磁場ベクトルまでの磁場の変化を表すベクトルである磁場変化ベクトルを求める手順と、
   現時点のタイムステップで求めた前記磁場変化ベクトルと、直前のタイムステップで求めた前記磁場変化ベクトルとのなす角度が９０°未満のときに、現時点のタイムステップで求めた前記磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを、前記第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる前記第１磁場ベクトルとして設定する手順と
   を含む磁場解析方法。
2. 現時点のタイムステップで求めた前記磁場変化ベクトルと、直前のタイムステップで求めた前記磁場変化ベクトルとのなす角度が９０°より大きいときに、現時点のタイムステップで求めた前記磁場変化ベクトルの向きを変えず、大きさを小さくしたベクトルを、前記第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる前記第１磁場ベクトルとして設定する手順を含む請求項１に記載の磁場解析方法。
3. 前記第２磁場ベクトルを求める手順において、前記複数の評価点に付与された磁場を最適解に収束させる運動方程式を数値的に解き、
   前記第２磁場ベクトルを求めた後、次のタイムステップを実行する前に、前記運動方程式の質量に相当するパラメータを変化させるＦＩＲＥ法を適用する請求項１または２に記載の磁場解析方法。
4. 前記運動方程式は、前記複数の多面体要素の各々の内部の磁化ベクトルの発散がゼロであるという束縛条件を表す束縛項を含んでおり、
   前記第２磁場ベクトルを求める手順において、
   現時点の前記第１磁場ベクトルを用いて、前記束縛条件を課さない条件で前記運動方程式を数値的に解いて、暫定磁場ベクトルを求め、
   前記束縛条件が満たされるまで、前記暫定磁場ベクトルを反復法により更新して前記第２磁場ベクトルを求める請求項３に記載の磁場解析方法。
5. さらに、前記磁場変化ベクトルを大きくする情報を取得する手順を含み、
   次のタイムステップにおいて用いる前記第１磁場ベクトルを設定する手順において、取得された情報に基づいて前記磁場変化ベクトルの大きさを大きくする請求項１乃至４のいずれか１項に記載の磁場解析方法。
6. 複数の多面体要素に分割された解析モデルの、前記複数の多面体要素の各々の内部の評価点の磁場を反復法により計算する磁場演算部を有する磁場解析装置であって、
   前記磁場演算部は、反復計算の１つのタイムステップにおいて、
   前記複数の評価点の各々の現時点の磁場を表す第１磁場ベクトルに基づいて、１タイムステップ経過後の前記複数の評価点の各々の磁場を表す第２磁場ベクトルを求める手順と、
   前記第１磁場ベクトルから前記第２磁場ベクトルまでの磁場の変化を表すベクトルである磁場変化ベクトルを求める手順と、
   現時点のタイムステップで求めた前記磁場変化ベクトルと、直前のタイムステップで求めた前記磁場変化ベクトルとのなす角度が９０°未満のときに、現時点のタイムステップで求めた前記磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを、前記第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる前記第１磁場ベクトルとして設定する手順と
   を実行する磁場解析装置。
7. 複数の多面体要素に分割された解析モデルの、前記複数の多面体要素の各々の内部の評価点の磁場を反復法により計算する手順をコンピュータに実行させるプログラムであって、
   反復計算の１つのタイムステップにおいて、
   前記複数の評価点の各々の現時点の磁場を表す第１磁場ベクトルに基づいて、１タイムステップ経過後の前記複数の評価点の各々の磁場を表す第２磁場ベクトルを求める手順と、
   前記第１磁場ベクトルから前記第２磁場ベクトルまでの磁場の変化を表すベクトルである磁場変化ベクトルを求める手順と、
   現時点のタイムステップで求めた前記磁場変化ベクトルと、直前のタイムステップで求めた前記磁場変化ベクトルとのなす角度が９０°未満のときに、現時点のタイムステップで求めた前記磁場変化ベクトルの向きを変えず、大きさを大きくしたベクトルを、前記第１磁場ベクトルに加えたベクトルを、次のタイムステップにおいて用いる前記第１磁場ベクトルとして設定する手順と
   をコンピュータに実行させるプログラム。
   

Images