(57)【要約】【課題】シミュレーションにおいて磁気的な現象を好適に扱える解析技術を提供する。【解決手段】解析装置１００は、仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部１１０と、磁気モーメント付与部１１０によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部１２１と、磁場演算部１２１における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部１２４と、を備える。磁場演算部１２１は、外部磁場の時間変動により各粒子に誘起される誘導磁化と、前記誘導磁化に基づいた磁気モーメント同士の相互作用により得られる磁場を用いて、誘導磁場を数値的に演算する。【選択図】図１
2015111401.tif 000002
【特許請求の範囲】
【請求項１】
仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、
前記磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、
前記磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、
を備え、
前記磁場演算部は、外部磁場の時間変動により各粒子に誘起される誘導磁化と、前記誘導磁化に基づいた磁気モーメント同士の相互作用により得られる磁場を用いて、誘導磁場を数値的に演算することを特徴とする解析装置。
【請求項２】
前記磁場演算部は、前記誘導磁化を表す項と、前記誘導磁化に基づいた磁気モーメントを有する粒子間の相互作用を表す磁場の補正項の双方を用いて、前記誘導磁場を演算することを特徴とする請求項１に記載の解析装置。
【請求項３】
仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、
前記磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、
前記磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、
を備え、
前記磁場演算部は、解析対象とする粒子群の重心位置を表す局所重心ベクトルを用いたゲージ変換により導出される誘導磁場の式を用いて、前記粒子群における誘導磁場を数値的に演算することを特徴とする解析装置。
【請求項４】
前記磁場演算部は、前記粒子系のベクトルポテンシャルに対して、前記局所重心ベクトルと前記粒子系の磁束密度の外積を用いて記述されるゲージに変換することにより導出される誘導磁場の式を用いて、前記粒子群における誘導磁場を数値的に演算することを特徴とする請求項３に記載の解析装置。
【請求項５】
前記粒子系には、電気的に互いが絶縁される第１粒子群および第２粒子群を含む複数の粒子群が配置可能とされており、
前記磁場演算部は、前記第１粒子群の重心位置を表す第１局所重心ベクトルを用いたゲージ変換により導出される第１の誘導磁場の式を用いて前記第１粒子群における誘導磁場を演算し、前記第２粒子群の重心位置を表す第２局所重心ベクトルを用いたゲージ変換により導出される第２の誘導磁場の式を用いて前記第２粒子群における誘導磁場を演算することを特徴とする請求項３または４に記載の解析装置。
【請求項６】
粒子系は繰り込み群分子動力学法を使用して繰り込まれた粒子系であることを特徴とする請求項１から５のいずれか一項に記載の解析装置。
【請求項７】
前記粒子状態演算部における支配方程式は磁気モーメントに依存する項を有することを特徴とする請求項１から６のいずれか一項に記載の解析装置。
【請求項８】
仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する機能と、
磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する機能と、
磁気的な物理量の演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する機能と、
をコンピュータに実現させ、
前記磁気的な物理量を演算する機能は、外部磁場の時間変動により各粒子に誘起される誘導磁化と、前記誘導磁化に基づいた磁気モーメント同士の相互作用により得られる磁場を用いて、誘導磁場を数値的に演算することを特徴とするコンピュータプログラム。
【請求項９】
仮想空間内に定義される粒子系の粒子に磁気モーメントを付与する機能と、
磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する機能と、
磁気的な物理量の演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する機能と、
をコンピュータに実現させ、
前記磁気的な物理量を演算する機能は、解析対象とする粒子群の重心位置を表す局所重心ベクトルを用いたゲージ変換により導出される誘導磁場の式を用いて、前記粒子群における誘導磁場を数値的に演算することを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子系を解析する解析装置に関する。
【背景技術】
【０００２】
近年、コンピュータの計算能力の向上に伴い、モータなどの電気機器の設計開発の現場では磁場解析を取り入れたシミュレーションがよく使用される。シミュレーションを使用すると、実際にプロトタイプを製作しなくてもある程度の評価が可能となるので、設計開発のスピードが向上しうる。
【０００３】
例えば特許文献１には、磁場解析を実行する演算処理装置を具えるモータ解析装置が記載されている。演算処理装置は、ユーザの操作に基づく外部指令に応じて、有限要素法による静磁場解析やマクスウェルの応力法によるトルク計算を実行する。有限要素法による静磁場解析のためにメッシュ分割が行われる。メッシュ分割は、コア領域及びハウジング領域、さらには外部空気層領域についても行われる。磁場解析の他の手法は差分法および磁気モーメント法を含む。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開平１１−１４６６８８号公報
【特許文献２】特開２００９−３７３３４号公報
【非特許文献】
【０００５】
【非特許文献１】Erik Bitzek他、「Structural Relaxation Made Simple」、Physical Review Letters、２００６年１０月２７日、９７
【非特許文献２】L.D.Landau、E.M.Lifshitz、L.P.Litaevskii著、「Electrodynamics of Continuous Media 2nd ed」、１９８４年１月１日、２０５ページ
【発明の概要】
【発明が解決しようとする課題】
【０００６】
これらの手法では解析対象をメッシュ分割する。しかしながら、解析対象が流体の場合固定メッシュが使われ、移動境界を伴うと解析が困難になる。また、解析対象が磁性体や塑性体の大変形を伴う場合もメッシュ分割は困難である。特に、モータなどの回転体を含む磁場解析でメッシュ分割を行おうとする場合、メッシュが重合（一重なら大変形）するため、解析は困難であった。
【０００７】
本発明はこうした課題に鑑みてなされたものであり、その目的は、シミュレーションにおいて磁気的な現象を好適に扱える解析技術の提供にある。
【課題を解決するための手段】
【０００８】
本発明のある態様の解析装置は、仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、を備える。磁場演算部は、外部磁場の時間変動により各粒子に誘起される誘導磁化と、誘導磁化に基づいた磁気モーメント同士の相互作用により得られる磁場を用いて、誘導磁場を数値的に演算する。
【０００９】
本発明の別の態様もまた、解析装置である。この装置は、仮想空間内に定義される粒子系に磁気モーメントを付与する磁気モーメント付与部と、磁気モーメント付与部によって磁気モーメントが付与された粒子を含む粒子系に関連する磁気的な物理量を演算する磁場演算部と、磁場演算部における演算結果を使用して、各粒子の運動を支配する支配方程式を数値的に演算する粒子状態演算部と、を備える。磁場演算部は、解析対象とする粒子群の重心位置を表す局所重心ベクトルを用いたゲージ変換により導出される誘導磁場の式を用いて、粒子群における誘導磁場を数値的に演算する。
【００１０】
なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。
【発明の効果】
【００１１】
本発明によれば、シミュレーションにおいて磁気的な現象を好適に扱うことができる。
【図面の簡単な説明】
【００１２】
【図１】実施の形態に係る解析装置の機能および構成を示すブロック図である。
【図２】図１の粒子データ保持部の一例を示すデータ構造図である。
【図３】図１の解析装置における一連の処理の一例を示すフローチャートである。
【図４】ｎ＝１次までの展開に対応する計算結果を示す図である。
【図５】ｎ＝２次までの展開に対応する計算結果を示す図である。
【図６】ｎ＝３次までの展開に対応する計算結果を示す図である。
【図７】ｎ＝４次までの展開に対応する計算結果を示す図である。
【図８】球面調和関数の次数と計算値との関係を示す図である。
【図９】ヒステリシスカーブを示す図である。
【図１０】ビーズにより構成される導体球の解析モデルを示す図である。
【図１１】導体球表面での時間変化する磁界の計算値を示すグラフである。
【図１２】計算値の位相誤差と振幅誤差を示すグラフである。
【図１３】導体球表面での磁界、電流および磁界の計算値を示すグラフである。
【発明を実施するための形態】
【００１３】
以下、各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。
【００１４】
（第１の実施の形態）
古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）や、量子分子動力学法（Quantum Molecular Dynamics Method）や、ＭＤ法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法（Renormalized Molecular Dynamics、以下ＲＭＤ法と称す）に基づくシミュレーションが知られている（例えば、特許文献２参照）。正確にはＭＤ法やＲＭＤ法は運動論的手法（物理量の算出には統計力学を使う）であり、粒子法は、連続体を記述する微分方程式を離散化する手法であり、別ものであるが、ここではＭＤ法やＲＭＤ法も粒子法と呼ぶ。
【００１５】
粒子法は静的な現象だけでなく流れなどの動的な現象をも取り扱えるので、主に静的な現象を解析対象とする上述の有限要素法などに代わるシミュレーション手法として注目されている。
【００１６】
粒子法には、連続体を粒子で離散化することにより解析対象の粒子系を得る、という微分的な見方がある。例えば、粒子法において流体を扱う場合、ナビエストークス（Navier-Stokes）方程式を粒子で離散化することが多い。
【００１７】
一方、粒子法の別の見方として、多くの粒子を集めて連続体を形成する、という積分的な見方もある。これは例えば、小さな鉄の粒を集めて固めて大きな鉄球を形成するという見方である。
【００１８】
一般に、多くの磁気モーメントが存在する空間内のある点の磁場を求める場合、重ね合わせの原理により、各磁気モーメントがその点に作る磁場を磁気モーメントに亘って足し合わせる。本発明者は、このような磁気モーメントの集まりから磁場を求める手法と、積分的な見方をした場合の粒子法と、の親和性を独自に見い出し、粒子法における粒子に磁気モーメントを付与することに想到した。これにより、対流や大変形の解析に強いという粒子法の利点を維持しつつ、粒子法の適用範囲を磁場解析にまで広げることが可能となる。
【００１９】
さらに、本発明者は、粒子法における各粒子に付与される磁気モーメント同士の相互作用項に加えて各粒子に誘起される誘導磁化を考慮に入れることで、粒子法の適用範囲を電磁誘導現象を含む動磁場解析にまで拡張できることに想到した。これにより、誘導モータにおける固定子と回転子の相互作用といった、大変形を伴うとともに電磁誘導の考慮が必要な物理現象に対して、粒子法による動磁場解析を実現することができる。
【００２０】
図１は、実施の形態に係る解析装置１００の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。
【００２１】
本実施の形態ではＲＭＤ法に倣って粒子系を解析する場合について説明するが、繰り込みを行わないＭＤ法やＤＥＭ（Distinct Element Method）やＳＰＨ（Smoothed Particle Hydrodynamics）やＭＰＳ（Moving Particle Semi-implicit）などの他の粒子法に倣って粒子系を解析する場合にも、本実施の形態に係る技術的思想を適用できることは本明細書に触れた当業者には明らかである。
【００２２】
解析装置１００は入力装置１０２およびディスプレイ１０４と接続される。入力装置１０２は、解析装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。
【００２３】
解析装置１００は、粒子系取得部１０８と、磁気モーメント付与部１１０と、数値演算部１２０と、表示制御部１１８と、粒子データ保持部１１４と、を備える。
【００２４】
粒子系取得部１０８は、入力装置１０２を介してユーザから取得する入力情報に基づき、１、２または３次元の仮想空間内に定義されるＮ（Ｎは自然数）個の粒子からなる粒子系のデータを取得する。粒子系はＲＭＤ法を使用して繰り込まれた粒子系である。
【００２５】
粒子系取得部１０８は、入力情報に基づき仮想空間内にＮ個の粒子を配置し、配置されたそれぞれの粒子に速度を付与する。粒子系取得部１０８は、配置された粒子を特定する粒子ＩＤと、その粒子の位置と、その粒子の速度と、を対応付けて粒子データ保持部１１４に登録する。
【００２６】
磁気モーメント付与部１１０は、入力装置１０２を介してユーザから取得する入力情報に基づき、粒子系取得部１０８によって取得された粒子系の粒子に磁気モーメントを付与する。例えば、磁気モーメント付与部１１０は、ディスプレイ１０４を介してユーザに、粒子系の粒子の磁気モーメントの入力を求める。磁気モーメント付与部１１０は、入力された磁気モーメントを粒子ＩＤと対応付けて粒子データ保持部１１４に登録する。
【００２７】
以下では粒子系の粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。
【００２８】
数値演算部１２０は、粒子データ保持部１１４によって保持されるデータが表す粒子系の各粒子の運動を支配する支配方程式を数値的に演算する。特に数値演算部１２０は、離散化された粒子の運動方程式にしたがった繰り返し演算を行う。粒子の運動方程式は磁気モーメントに依存する項を有する。
数値演算部１２０は、磁場演算部１２１と、力演算部１２２と、粒子状態演算部１２４と、状態更新部１２６と、終了条件判定部１２８と、を含む。
【００２９】
磁場演算部１２１は、粒子データ保持部１１４によって保持されるデータが表す粒子系が生成する磁場を演算する。磁場演算部１２１は、粒子データ保持部１１４によって保持される各粒子の磁気モーメントに基づいて磁場を演算する。磁場は、粒子系に関連する磁気的な物理量である。磁場演算部１２１は磁場の代わりにまたは磁場に加えて磁束密度や磁化を演算してもよい。
【００３０】
力演算部１２２は粒子データ保持部１１４によって保持される粒子系のデータを参照し、粒子系の各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。粒子に働く力は、磁気モーメント間の相互作用に基づく力を含む。力演算部１２２は、粒子系のｉ番目（１≦ｉ≦Ｎ）の粒子について、そのｉ番目の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。
【００３１】
力演算部１２２は、各近接粒子について、その近接粒子とｉ番目の粒子との間のポテンシャルエネルギ関数およびその近接粒子とｉ番目の粒子との距離に基づいて、その近接粒子がｉ番目の粒子に及ぼす力を演算する。特に力演算部１２２は、その近接粒子とｉ番目の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。力演算部１２２は、近接粒子がｉ番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、ｉ番目の粒子に働く力を算出する。
【００３２】
粒子状態演算部１２４は粒子データ保持部１１４に保持される粒子系のデータを参照し、粒子系の各粒子について、離散化された粒子の運動方程式に力演算部１２２によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。
【００３３】
粒子状態演算部１２４は、力演算部１２２によって演算された力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき所定の微小な時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって演算された力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。
【００３４】
粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、粒子系のｉ番目の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。
【００３５】
状態更新部１２６は、粒子データ保持部１１４に保持される粒子系の各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。
【００３６】
終了条件判定部１２８は、数値演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、外部から終了の指示を受け付けたことや、粒子系が定常状態に達したことである。終了条件判定部１２８は、終了条件が満たされる場合、数値演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を力演算部１２２に戻す。すると力演算部１２２は、状態更新部１２６によって更新された粒子の位置、速度で再び力を演算する。
【００３７】
表示制御部１１８は、粒子データ保持部１１４に保持されるデータが表す粒子系の各粒子の位置、速度、磁気モーメントに基づき、ディスプレイ１０４に粒子系の時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。
【００３８】
図２は、粒子データ保持部１１４の一例を示すデータ構造図である。粒子データ保持部１１４は、粒子ＩＤと、粒子の位置と、粒子の速度と、粒子の磁気モーメントと、を対応付けて保持する。
【００３９】
上述の実施の形態において、保持部の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。
【００４０】
以上の構成による解析装置１００の動作を説明する。
図３は、解析装置１００における一連の処理の一例を示すフローチャートである。粒子系取得部１０８は、ＲＭＤ法に倣って繰り込まれた粒子系を取得する（Ｓ１２）。磁気モーメント付与部１１０は、取得された粒子系の粒子に磁気モーメントを付与する（Ｓ１４）。磁場演算部１２１は、粒子系によって生成される磁場を演算する（Ｓ１６）。力演算部１２２は、粒子間の距離と各粒子の磁気モーメントから粒子に働く力を演算する（Ｓ１８）。粒子状態演算部１２４は、演算された力を含む粒子の運動方程式から粒子の速度と位置を演算する（Ｓ２０）。状態更新部１２６は、粒子データ保持部１１４に保持される粒子の位置、速度を演算された位置、速度で更新する（Ｓ２２）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２４）。終了条件が満たされない場合（Ｓ２４のＮ）、処理はＳ１６〜Ｓ２２のステップを繰り返す。終了条件が満たされる場合（Ｓ２４のＹ）、表示制御部１１８は、演算結果を出力する（Ｓ２６）。
【００４１】
粒子法は一般に、流体や分離現象や切断加工などの動的な解析対象をより好適に扱える。
そして本実施の形態に係る解析装置１００によると、粒子法による解析と磁気モーメントに基づく磁場解析とを連成することができる。したがって、動的な解析対象を、磁場などの磁気的な性質も含めてより精度高くかつ短時間で解析することが可能となる。
【００４２】
例えば解析対象が砂鉄の場合、従来の有限要素法などのメッシュベースの解析手法ではメッシュを好適に設定することができないので、うまくシミュレーションすることが困難である。そこで本実施の形態に係る解析装置１００によると、砂鉄の磁性を考慮に入れた形で砂鉄を粒子法で好適に解析することができる。
【００４３】
また例えば解析対象がモータの場合、固定体に対して回転体が回転するので、やはりメッシュベースの解析手法は不向きであり静解析に止まる。そこで、本実施の形態に係る解析装置１００によると、回転体の回転を粒子法により高精度にシミュレーションしつつ、回転体と固定体との磁気的相互作用を動的に扱うことができる。さらに、誘導モータのように電磁誘導現象を含む解析対象であっても、電磁誘導現象を考慮に入れた誘導磁化と、誘導磁化に基づいた磁気モーメントを付与することにより、磁場の時間変化を加味した磁気的相互作用を扱うことができる。
【００４４】
以下、解析装置１００で使用される解析手法の原理について説明する。
【００４５】
磁場はマクロなベクトルポテンシャル：
【数１】
2015111401.tif 000003
から以下の様に求まる。スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにある１個のスピンが作る磁場はスピンＳｊｚ（磁気モーメントｍｊｚ）を発生させる円電流を
【数２】
2015111401.tif 000004
とし、
【数３】
2015111401.tif 000005
ここで、
【数４】
2015111401.tif 000006
ｎ＝０のみを残すと一様に磁化：
【数５】
2015111401.tif 000007
した球の作る磁場と同等になる。
なお、真球のバルク材に本手法を適用する場合、ｎ＝０のみ採用すると、一様な磁化が得られず、磁気モーメント法と同様な結果になる。高次の項を加える事で改善できる。ＦＥＭは一様になる。
【００４６】
デカルト座標系に変換（
【数６】
2015111401.tif 000008
を使って、)すれば
【数７】
2015111401.tif 000009
ここで、ｓｉｎθｊ＝０のときは、Ｂθ＝０から、Ｂｘ＝Ｂｙ＝０であることに注意する。磁気モーメントｍｘ、ｍｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場
【数８】
2015111401.tif 000010
も同様に求まる。
【数９】
2015111401.tif 000011
【００４７】
【数１０】
2015111401.tif 000012
における磁場
【数１１】
2015111401.tif 000013
（Ｎ個のスピンからの寄与）は
【数１２】
2015111401.tif 000014
ここで、
【数１３】
2015111401.tif 000015
は全外場を意味する。Ｐ−１４とＰ−１５を連立することで, 磁気モーメントｍｐ が求まる。磁性体外部に生じる磁場は、Ｐ−１５であり、磁性体内部に生じる磁場は、
【数１４】
2015111401.tif 000016
Ｐ−１４は球体表面に生じる磁化による反磁場が考慮され、
【数１５】
2015111401.tif 000017
を満足する。
【００４８】
Ｐ−１３はマトリクスが優対角でない(ｊ＝ｐが除外されている) から収束は遅い。そこで、以下の運動方程式から解を得る。
【数１６】
2015111401.tif 000018
ここで、ｍｖ＝５×１０−２ｄｔ２は仮想質量、γ＝ｍｖ／（１０ｄｔ）は、減衰係数である。
【００４９】
ＦＩＲＥ（非特許文献１参照）を追加すればさらに高速化が可能。以下にＦＩＲＥのコードを示す。なお、粒子数８０２、残差＜１０−８における計算で１０．０５[ｓ］である。ＦＩＲＥ無しでは、残差が１０−５以下にならなかった。
【数１７】
2015111401.tif 000019
【００５０】
以下にLegendre 関数と倍Legendre 関数の特別な場合を列挙する。（ｘ：＝ｃｏｓθｊ）
【数１８】
2015111401.tif 000020
【００５１】
以下に、ｘ軸方向にＢｏｘ＝０．１［Ｔ］の一様な外場においた磁性球（帯磁率が９９９）の磁化
【数１９】
2015111401.tif 000021
を計算した結果を示す。球を構成する粒子（原子）数は４００９個とし、ｆｃｃ構造に配置した。厳密解は
【数２０】
2015111401.tif 000022
である。図４、図５、図６、図７に多重極展開の１次〜４次までの結果を示す。共に中心を通過する断面を表示した。多重極展開の次数を上げると, 厳密解と良く合う事が判る。次数が低いと収束が悪い。また計算時間の大半は球面調和関数の呼び出しである。
【００５２】
図４は、ｎ＝１次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３４６２３１[Ｔ］である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。残差１ｅ−１３以下で、計算時間は約８０分であった。
【００５３】
図５は、ｎ＝２次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３２１８９１[Ｔ]である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。残差１ｅ−８以下で、計算時間は約３ 時間であった。残差は１ｅ−９以下に落ちなかった。
【００５４】
図６は、ｎ＝３次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３１２４４４[Ｔ]である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる.
【００５５】
図７は、ｎ＝４次までの展開に対応する計算結果を示す図である。厳密解μｏＭｘ＝０．２９９１[Ｔ］に対し計算値は平均０．３１２２２２[Ｔ]である。磁化は外場Ｂｘ＝０．１[Ｔ]に平行かつ一様であることが確認できる。約５時間後において、残差は１ｅ−９以下に落ちなかった。
【００５６】
図８は、球面調和関数の次数と計算値との関係を示す図である。次数を上げると厳密解に近づくことが分かる。
【００５７】
解析装置１００で使用される解析手法の原理を、次のように説明することもできる。
【００５８】
【数２１】
2015111401.tif 000023
Ｍｐは原子核の質量、ｍｅは電子の質量である。ｅは原子核と電子の電荷(の絶対値)である。最終的に電荷は電流に置き換わり、電子の質量は原子核の質量と和して原子の質量Ｍとなるから、陽には現れない。δＶは磁性球の体積、ｍｖは仮想質量、λｉはラグランジェの未定常数である。厳密に拘束を満足させるには、ＳＨＡＫＥ法に代表される収束計算が必要であるが、本手法ではλｉ＝１として、正解の周りの微小振動を誤差として扱う。
【００５９】
スピンの向きをｚ軸に取り、ａを原子半径とする。格子点ｊにあるスピンＳｊｚ（磁気モーメントｍｊｚ）を発生させる円電流は
【数２２】
2015111401.tif 000024
であるから、磁場はマクロなベクトルポテンシャル：
【数２３】
2015111401.tif 000025
から以下の様に求まる。
【数２４】
2015111401.tif 000026
ここで、
【数２５】
2015111401.tif 000027
ｎ＝０のみを残すと一様に磁化：
【数２６】
2015111401.tif 000028
した球の作る磁場と同等になる。
【００６０】
デカルト座標系に変換(
【数２７】
2015111401.tif 000029
を使って、）すれば
【数２８】
2015111401.tif 000030
ここで、ｓｉｎθｊ＝０のときは、Ｈθ＝０から、Ｈｘ＝Ｈｙ＝０であることに注意する。磁気モーメントｍｘ、ｍｙがｘ軸、ｙ軸にそれぞれ平行なときの磁場
【数２９】
2015111401.tif 000031
も同様に求まる。
【数３０】
2015111401.tif 000032
【００６１】
【数３１】
2015111401.tif 000033
におけるモーメントが感じる外磁場
【数３２】
2015111401.tif 000034
(Ｎ個のスピンからの寄与) は
【数３３】
2015111401.tif 000035
誘導磁場
【数３４】
2015111401.tif 000036
は、ｐ点が導体内部である場合、
【数３５】
2015111401.tif 000037
ｐ点が導体外部である場合、
【数３６】
2015111401.tif 000038
Ｃ１、Ｃ２は、例えば、非特許文献２に示される式において
【数３７】
2015111401.tif 000039
となる極限と比較したときの補償係数であり、Ｃ１＝１／３、Ｃ２＝３／５の値をとる。
【００６２】
ここで、Ｑ−２１の右辺第１項（またはＱ−２１−２）は、時間変化する外部磁場により各粒子に誘起される誘導磁化Ｍａｎｔを表す項である。また、Ｑ−２１の右辺第２項およびＱ−２２の右辺は、Ｑ−２１−４式で表される各粒子の誘導磁化に基づいた磁気モーメントｍａｎｔ同士の相互作用により得られる磁場を表す項である。誘導磁場をこのように表すことで、時間変化する外部磁場による影響を考慮した動磁場解析を実現することができる。
【００６３】
なお、繰り込みに際し、電気伝導度σは
【数３８】
2015111401.tif 000040
であり、今はδ＝２としている。
【００６４】
Ｑ−２１−２右辺のLagrange 微分
【数３９】
2015111401.tif 000041
は、粒子法であるから
【数４０】
2015111401.tif 000042
で置き換えることが出来る。ただし、磁化の並進移動は粒子の移動により考慮されるが、磁化ベクトルの回転は別に考慮する必要がある(ランダウ=リフシッツ、電磁気学、§５１)。
【００６５】
磁性体内部に生じる磁場は、
【数４１】
2015111401.tif 000043
であり、全外場ＨＯ（ｘｐ）と、磁化Ｍ（Ｈ）との関係は、
【００６６】
【数４２】
2015111401.tif 000044
である。Ｑ−２６とＱ−２７を連立することで、磁気モーメント
【数４３】
2015111401.tif 000045
が求まる。また、磁性体外部の磁場は、
【数４４】
2015111401.tif 000046
となる。
【００６７】
非線形磁化とヒステリシス
Ｑ−２７とＱ−２４から磁場
【数４５】
2015111401.tif 000047
を消去する：
【数４６】
2015111401.tif 000048
Ｑ−２８を解くことで
【数４７】
2015111401.tif 000049
が得られる（図９にはＢ＝ｆ（Ｈ）の場合を示す）。
【００６８】
図９は、ヒステリシスカーブを示す図である。ヒステリシスカーブとＢ＋２μｏＨ＝３Ｂｏの交点からＭを求める。今
【数４８】
2015111401.tif 000050
を永久磁化とするＢ−Ｈカーブ：
【数４９】
2015111401.tif 000051
が与えられたとする。
【数５０】
2015111401.tif 000052
と連立し
【数５１】
2015111401.tif 000053
を得る。一般のＢ＝ｆ（Ｈ）に対しては、Newton-Laphson法で解く。
【００６９】
磁場の粒子化法
Ｑ−１９はマトリクスが優対角でない（ｊ＝ｐが除外されている）から収束は遅い。そこで、以下の運動方程式から解を得る。
【数５２】
2015111401.tif 000054
ここで、ＨＯ（ｘｐ）は、Ｑ−２６により表される。また、ｍｖは仮想質量である。推奨値は２５０×ｄｔ×ｄｔである。ｍｖを十分小さく取れば、正解の周りの微小な減衰振動が起こり、誤差を小さくできる。
【００７０】
模型を考える：
【数５３】
2015111401.tif 000055
減衰項を持つ強制振動の解(Landau and Lifshitz,Mechanics,Sec.26) は
【数５４】
2015111401.tif 000056
【数５５】
2015111401.tif 000057
を考えれば、
【数５６】
2015111401.tif 000058
と成る。即ちｍｖを十分小さく取れば、解くべき本来の方程式
【数５７】
2015111401.tif 000059
の解に一致する。
【００７１】
離散化
蛙飛び法（速度はｎとｎ＋１の中間点における評価）に従い離散化すれば、収束回数をｓと書き、
【数５８】
2015111401.tif 000060
ここで、誘導磁場Ｈｉｎｄ（ｘｉ）の離散化において、
【数５９】
2015111401.tif 000061
ととる。このとき、対角要素のみ未来を使えば、陽的解法に出来、約３０％程度効率改善になる。計算の大半は特殊関数の呼び出しである。
【００７２】
ＦＩＲＥによる高速化：静磁場解析にのみ有効
ＦＩＲＥ（非特許文献１参照）を追加すればさらに高速化が可能。以下にＦＩＲＥのコードを示す。なお、粒子数８０２、残差＜１０−８における計算で１０．０５[ｓ]である。ＦＩＲＥ無しでは、残差が１０−５以下にならなかった。
【数６０】
2015111401.tif 000062
【００７３】
精度検証：球体の一様磁化
図１０は、ビーズにより構成される導体球の解析モデルを示す図である。精度検証では、ｘ軸方向にＨｅｘｔ，ｘ＝ＨＯｓｉｎ（２πｆｔ）となる空間一様な外部磁界を印加した。ここで、ＨＯ＝７．９５８×１０４［Ａ／ｍ］、ｆ＝２００［Ｈｚ］、導体球の半径Ａ＝１０［ｍｍ］、電気電導率σ＝５９×１０６［Ｓ／ｍ］とした。導体球を構成する粒子（原子）数は、６０９９個とし、ｆｃｃ構造に配置した。
【００７４】
図１１は、導体球表面での時間変化する磁界の計算値を示すグラフであり、導体表面近傍に位置する粒子の重心位置（ｒｘ＝９．８９［ｍｍ］、ｒｙ＝ｒｚ＝０）上の磁界Ｈｘの時間変化を示す。なお、本図では、シミュレーションによる計算結果とともに、外部磁界および厳密解による値も示している。誘導磁界の影響により、導体内部の磁界が外部磁界に対して位相が遅れて変化する様子を再現することができており、計算結果と厳密解の傾向は良く一致することが示された。
【００７５】
図１２は、計算値の位相誤差と振幅誤差を示すグラフであり、ｘ軸上に位置する粒子の重心位置における磁界の位相誤差εｐｈａｓｅと振幅誤差εａｍｐｌｉｔｕｄｅを示す。位相誤差εｐｈａｓｅおよび振幅誤差εａｍｐｌｉｔｕｄｅは、
【数６１】
2015111401.tif 000063
と定義した。ここで、θｅｘａｃｔは厳密解の位相、θＭＢＭは解析結果の位相である。また、Ｈｅｘａcｔは厳密解により得られた磁界、ＨＭＢＭは解析結果の磁界である。なお位相誤差は、Ｈｅｘａcｔが最大となる位相θｅｘａｃｔと、ＨＭＢＭが最大となる位相θＭＢＭとを比較した。図１２より、ｘ軸上における解析結果と厳密解の磁界の位相誤差は最大で５．８６％、振幅誤差は６．３８％となり、厳密解と比較して精度の高い解析結果が得られることが示された。
【００７６】
Legendre関数の級数表示
以下にLegendre関数と倍Legendre関数の特別な場合を列挙する。（ｘ：＝ｃｏｓθｊ）
【数６２】
2015111401.tif 000064
【００７７】
粒子に働く力の導出
粒子（電子と核）と、粒子磁場間の相互作用を記述するLagrangeanは
【数６３】
2015111401.tif 000065
である。ベクトルポテンシャル
【数６４】
2015111401.tif 000066
と電磁場の関係は、
【数６５】
2015111401.tif 000067
である。
【数６６】
2015111401.tif 000068
を使い、以下の運動方程式を得る。
【数６７】
2015111401.tif 000069
原子Ｍ＝Ｍｐ＋ｍｅに対する運動方程式は、
【数６８】
2015111401.tif 000070
となる。右辺：Lorentz力を平均電流密度
【数６９】
2015111401.tif 000071
、電気モーメントを
【数７０】
2015111401.tif 000072
で置き換えると
【数７１】
2015111401.tif 000073
Maxwellの方程式とＱ−２４から
【数７２】
2015111401.tif 000074
最終的に原子に働く力は、電気分極を無視し、
【数７３】
2015111401.tif 000075
となる。
【００７８】
粒子に働く力は、直接Maxwellの応力テンソルからも求まる。
【数７４】
2015111401.tif 000076
球に局在した電流が作るモーメント
【数７５】
2015111401.tif 000077
に作用する力は、最低の次数では
【数７６】
2015111401.tif 000078
である。従い、ポテンシャルエネルギーは、電気分極によるポテンシャルエネルギーを合わせて、
【数７７】
2015111401.tif 000079
となるが、渦電流や真電流を包含しない。
【００７９】
【数７８】
2015111401.tif 000080
結局、
【数７９】
2015111401.tif 000081
を計算すれば磁性球に働く磁力が求まる。
【００８０】
磁化による電流
【数８０】
2015111401.tif 000082
と渦電流
【数８１】
2015111401.tif 000083
を完全に切り分けることは恐らく出来ないが、それぞれ、
【数８２】
2015111401.tif 000084
【００８１】
微係数の計算
以下に必要な微分係数を書き出す。
【数８３】
2015111401.tif 000085
等からｒｐｊ＞ａのみ列挙すると（なお、
【数８４】
2015111401.tif 000086
はｓｉｎθに比例するから、θ＝０で発散は無い）、
【数８５】
2015111401.tif 000087
【数８６】
2015111401.tif 000088
【００８２】
磁気モーメントが作る磁化電流の計算
ｚ軸に平行に磁化したとき
今一度、ｚ軸に磁化したときの磁場を書く
【数８７】
2015111401.tif 000089
磁化電流を
【数８８】
2015111401.tif 000090
と書けば、磁化により生じる磁力
【数８９】
2015111401.tif 000091
は
【数９０】
2015111401.tif 000092
である。先ず、
【数９１】
2015111401.tif 000093
を求める。
【数９２】
2015111401.tif 000094
【００８３】
ｘ軸に平行に磁化したとき
ｘ軸に磁化したときの磁場を書く：
【数９３】
2015111401.tif 000095
【数９４】
2015111401.tif 000096
磁化電流
【数９５】
2015111401.tif 000097
を求める。
【数９６】
2015111401.tif 000098
【００８４】
ｙ軸に平行に磁化したとき
ｙ軸に磁化したときの磁場を書く：
【数９７】
2015111401.tif 000099
磁化電流
【数９８】
2015111401.tif 000100
を求める。
【数９９】
2015111401.tif 000101
【００８５】
誘導磁界が作る渦電流の計算
導体内部の誘導磁場は、Ｑ−２１より、
【数１００】
2015111401.tif 000102
となる。渦電流は
【数１０１】
2015111401.tif 000103
なお、上式Ｒ−５２−２は、Ｒ−１２において、球面調和関数の展開次数ｎ＝０とした場合と同じである。
【００８６】
粒子に働く力の計算
ｍｘ、ｍｙ、ｍｚと渦電流による寄与を加えれば、全電流が得られる。
【数１０２】
2015111401.tif 000104
座標
【数１０３】
2015111401.tif 000105
に粒子に働く力は
【数１０４】
2015111401.tif 000106
最後に
【数１０５】
2015111401.tif 000107
の微分から出てくる項を付加する必要がある。これは磁場の誤差
【数１０６】
2015111401.tif 000108
から生じる見かけの力であり、運動と磁場の全エネルギーを保存させる。
【数１０７】
2015111401.tif 000109
【００８７】
（第２の実施の形態）
上述の第１の実施の形態では、Ｑ−２１の右辺第１項に示される誘導磁化と、右辺第２項に示される磁気モーメントの相互作用とにより表される誘導磁場の式（Ｑ−２１）を用いて動磁場解析を行う場合を示した。本実施の形態では、Ｑ−２１の式を用いる代わりに、ベクトルポテンシャルの微分方程式から導出されるベクトルポテンシャルの厳密解を利用して誘導磁場を定式化する。以下、第１の実施の形態との相違点を中心に述べる。
【００８８】
ベクトルポテンシャルＡについて、以下の微分方程式が成立する。
【数１０８】
2015111401.tif 000110
上記微分方程式の解は、
【数１０９】
2015111401.tif 000111
となる。ここで、
【数１１０】
2015111401.tif 000112
として、各粒子に分割すれば、各粒子の位置ｒｉでのベクトルポテンシャルは、
【数１１１】
2015111401.tif 000113
となる。ここで、ｖｊは各粒子の体積である。なお、ｊの和に際し、ベクトル
【数１１２】
2015111401.tif 000114
が粒子表面を横切る場合には和を取らない。また、粒子径ａが満たすべき条件は、
【数１１３】
2015111401.tif 000115
である。
【００８９】
このとき、各対象物を粒子に分割することにより、ベクトルポテンシャルＡの並進対称性を失うおそれがある。そこで、本発明者は、下記Ｓ−４に示されるゲージ変換を導入することにより、ベクトルポテンシャルＡの並進対称性を回復できることを見いだした。なお、このゲージ変換を本明細書において「ビーズゲージ」ともいう。
【数１１４】
2015111401.tif 000116
ここで、ベクトルｇｉは、計算対象とする粒子群（粒子数：Ｎｉ）に関する局所重心ベクトルであり、
【数１１５】
2015111401.tif 000117
で表される。したがって、ベクトルｇｉは、粒子群の重心位置を表している。また、Ｓ−４に示すビーズゲージは、磁束密度Ｂｊと局所重心ベクトルｇｉの外積を用いて記述されると言える。
【００９０】
上記のビーズゲージについて、磁束密度ＢとベクトルポテンシャルＡの関係式
【数１１６】
2015111401.tif 000118
を求めれば、
【数１１７】
2015111401.tif 000119
となる。したがって、粒子数Ｎｉが十分に大きければ、並進対称性が満たされることがわかる。
【００９１】
ここで、Ｓ−４に示すゲージ変換から、Ｓ−３の右辺に含まれるベクトルポテンシャルの時間微分
【数１１８】
2015111401.tif 000120
を磁束密度Ｂの時間微分を用いて以下のように表すことができる。
【数１１９】
2015111401.tif 000121
さらに、ベクトルｄｊｉを、
【数１２０】
2015111401.tif 000122
とすると、ベクトルポテンシャルの時間微分は、
【数１２１】
2015111401.tif 000123
と表される。したがって、Ｓ−３にＳ−９を代入して積分を実行すれば、ベクトルポテンシャルＡの厳密解を磁束密度Ｂの時間微分を用いて記述できる。
【００９２】
ベクトルポテンシャルＡの厳密解を導出するため、まず、磁束密度Ｂのｚ成分の時間微分である、
【数１２２】
2015111401.tif 000124
のみが作用する場合の積分を実行する。これにより、ベクトルポテンシャルＡは、
【数１２３】
2015111401.tif 000125
となる。磁束密度Ｂのｘ成分およびｙ成分の時間微分である、
【数１２４】
2015111401.tif 000126
の寄与も加えれば、
【数１２５】
2015111401.tif 000127
が得られる。
【００９３】
解析対象とする導体内部を記述する誘導磁場Ｈは、
【数１２６】
2015111401.tif 000128
の関係式より、
【数１２７】
2015111401.tif 000129
となる。ここで、
【数１２８】
2015111401.tif 000130
と置換すると、誘導磁場Ｈの式として、
【数１２９】
2015111401.tif 000131
が得られる。
【００９４】
磁場演算部１２１は、粒子系により生成される誘導磁場を、Ｓ−１３に示す誘導磁場の式を用いて演算する。誘導磁場の時間変化値を得るための計算は、上述の第１の実施の形態で示したものと同様の方法を用いればよい。磁場演算部１２１は、得られた誘導磁場から電流密度を計算してもよい。また、力演算部１２２は、得られた磁場および電流の値から、各粒子に作用する電磁力を計算してもよい。
【００９５】
本実施の形態において、Ｓ−５に示す局所重心ベクトルｇｉは、互いに同電位となる粒子群に対して導入され、電気的に互いに絶縁される複数の粒子群に対しては、それぞれ別の局所重心ベクトルが導入される。例えば、解析対象とする粒子系において、互いが電気的に絶縁される複数の部分を有する場合、それぞれの部分を構成する粒子群に対して重心位置の異なる局所重心ベクトルが適用される。また、解析対象とする物体が時間経過とともに変形して複数の物体に分裂し、分裂した部分同士が互いに絶縁した状態となる場合には、それぞれの部分を構成する粒子群に対して重心位置の異なる局所重心ベクトルを適用する。
【００９６】
したがって、解析対象の粒子系に、電気的に互いが絶縁される第１粒子群および第２粒子群を含む複数の粒子群が配置される場合、磁場演算部１２１は、それぞれの粒子群に対応する誘導磁場の式を用いて解析をする。例えば、第１粒子群における誘導磁場は、第１粒子群の重心位置を表す第１局所重心ベクトルを用いたゲージ変換により導出される第１の誘導磁場の式を用いて演算される。また、第２粒子群における誘導磁場は、第２粒子群の重心位置を表す第２局所重心ベクトルを用いたゲージ変換により導出される第２の誘導磁場の式を用いて演算される。より具体的には、Ｓ−１３に示す誘導磁場の式において、右辺第１項に含まれるｄｉｊがそれぞれの粒子群での計算において異なる値をとることとなる。
【００９７】
本実施の形態に係る誘導磁場の式（Ｓ−１３）を用いた解析手法について精度検証を行った。第１の実施の形態と同様に、図１０に示す導体球の解析モデルを用い、ｘ軸方向にＨｅｘｔ，ｘ＝ＨＯｓｉｎ（２πｆｔ）となる空間一様な外部磁界を印加した。計算条件は、ＨＯ＝７．９４８×１０４［Ａ／ｍ］、ｆ＝１００［Ｈｚ］、導体球の半径Ａ＝１０［ｍｍ］、電気電導率σ＝５９×１０６［Ｓ／ｍ］とした。導体球を構成する粒子（原子）数は、６０９９個とし、ｆｃｃ構造に配置した。
【００９８】
図１３は、導体球表面での磁界、電流密度および電磁力の計算値を示すグラフである。本図は、導体表面近傍に位置する粒子の重心位置における、（ａ）磁界Ｈｘ、（ｂ）電流密度ｊｚ、（ｃ）電磁力Ｆｙの時間変化を示す。本図では、厳密解の値を実線で、計算値を破線で示している。図示されるように、本実施の形態に係る解析手法を用いると、計算結果と厳密解がよく一致することがわかった。
【００９９】
以上、実施の形態に係る解析装置１００の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。
【０１００】
実施の形態では、数値演算部１２０において粒子の位置と速度の両方を演算する場合について説明したが、これに限られない。例えば、数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法もあり、本実施の形態に係る技術的思想をそのような手法に適用してもよい。
【符号の説明】
【０１０１】
１００  解析装置、  １０２  入力装置、  １０４  ディスプレイ、  １０８  粒子系取得部、  １１０  磁気モーメント付与部、  １１４  粒子データ保持部、  １１８  表示制御部、  １２０  数値演算部、  １２１  磁場演算部、  １２２  力演算部、  １２４  粒子状態演算部、  １２６  状態更新部、  １２８  終了条件判定部。
【図１】
2015111401.tif 000132
【図２】
2015111401.tif 000133
【図３】
2015111401.tif 000134
【図８】
2015111401.tif 000135
【図９】
2015111401.tif 000136
【図１１】
2015111401.tif 000137
【図１２】
2015111401.tif 000138
【図１３】
2015111401.tif 000139
【図４】
2015111401.tif 000140
【図５】
2015111401.tif 000141
【図６】
2015111401.tif 000142
【図７】
2015111401.tif 000143
【図１０】
2015111401.tif 000144