積分方程式法による大規模電磁界数値解析の実用化に関する研究

(1)

積分方程式法による大規模電磁界数値解析の実用化に関する研究

Large-Scale Electromagnetic Field Analysis by Using Integral Equation Method

for Practical Electric Machine Design

2008 年 2 月

早稲田大学大学院理工学研究科電気・情報生命専攻

コンピュータ援用電磁工学研究

高橋康人

(2)

(3)

電子計算機の目覚しい発達と様々な電磁界数値計算技術の研究開発は，電気機器の設計や物理現象の解明に画期的な進歩をもたらした。最近は，相当に複雑な形状や構造の解析まで可能になっている。しかし，設計対象や解析対象がますます複雑化，大規模化し，数値解析におけるさらなる高精度化と計算時間の短縮化が要求されている状況にある⁽¹⁾。

電磁界数値解析分野においては，多媒質問題や不均質・非線形・異方性などの複雑な媒質特性を有する問題に適用可能であるなど，その汎用性の高さから有限要素法⁽²⁾が広く普及しているが，空気領域の要素分割が必要であるため，移動物体や無限領域の取扱いに難点がある。一方，境界要素法⁽³⁾に代表される積分方程式法は，複雑な媒質特性の処理は不得手であるが，空気領域の要素分割が不要なため，移動物体や無限領域の解析が容易である。

これらの特徴は相補う関係にあり，対象に応じて解析手法を使い分けることが，電磁界数値解析に基づく電気機器設計の高度化において重要と考えられる。電磁現象は本質的にその影響が無限遠にまで及ぶため，電気機器本体のみならず，周囲の広範な空気領域まで考慮しなければならない。このような特徴を有する電磁界問題においては，積分方程式法を活用する利点は非常に大きい。しかしながら，積分方程式法は密行列を扱うため膨大な計算時間と記憶容量を要する欠点があり，複雑形状を有する実機を対象とした大規模解析は困難とされていた。

ところが近年，必要な精度を確保しつつ多体間相互作用の近似値を高速に計算する手法である高速多重極法（FMM）⁽⁴⁾⁽⁵⁾や，マトリクス自体の代数的な操作により近似行列を作成するAdaptive Cross Approximation（ACA）⁽⁶⁾⁽⁷⁾などの行列圧縮アルゴリズムが様々な分野で適用され，積分方程式法を用いた大規模解析が報告されるようになってきた^(8)-(10)。しかし，

これらの手法を積分方程式法と組み合わせ電気機器設計に展開・実用化する際には，単なる数値計算の大規模化だけでは不十分であり，渦電流，磁性材料の磁気特性の非線形性・

異方性，積層鉄芯や磁気シールドのようなアスペクト比の非常に高い電気機器部材などといった機器解析特有の困難さを克服することが必須となる。すなわち，形状のみならず包含される物理現象が複雑でありかつ要求される計算の精度・時間の制約が厳しい実規模レベルの複雑な問題への適用に際し，上述の高速大規模数値計算技術を組み合わせた電磁界解析におけるさらなる技術の深度化が必要と考えられる。

このような背景のもと，本研究では積分方程式法を活用し，実用的な電気機器設計に耐えうる実規模電磁界解析手法の開発を目的として，上記行列圧縮アルゴリズムの導入による積分方程式法の大規模高速化を図るとともに，電磁現象特有の困難さを解決する種々の数値計算技術を開発・提案し，その有効性を明らかにした。

(8)

1.2 本論文の構成

本論文は，6章で構成されており，各章の要約は以下のとおりである。

第1章では，本研究の背景と目的について述べる。

第 2 章では，電気機器を設計する上で有用な情報を高精度に抽出する一手法として，渦電流場における導体表面電荷の大規模解析について述べる。導体表面電荷は導体中に発生する渦電流分布に大きく依存するため，例えば機器内における導体形状の変化が機器性能に及ぼす影響の評価をはじめ，様々な機器特性の情報抽出を行う上で有効に活用できる⁽¹¹⁾。この導体表面電荷は，通常の辺要素有限要素法の解析結果に後処理を施し，クーロンゲージを付加した導体表面上の電気スカラポテンシャルから，積分方程式法を用いて算出される。したがって，離散化後に得られるマトリクスが密行列となり，膨大な計算時間と記憶容量が必要となる。これが，回転機のような複雑形状へ導体表面電荷解析を適用する際の大きな障壁となっていた。そこで本章では，実規模問題に適用可能な導体表面電荷解析技術の確立を目的とした種々の検討を行った。まず，代表的行列圧縮アルゴリズムである高速多重極法およびACAを導入し，導体表面電荷解析の大規模化を達成した。また，さらなる高速化を達成するため，高速多重極法において解析モデルの対称性・周期性を活用する方法を新たに提案した。最終的に，複雑形状を有する実規模モデルとして，かご形誘導電動機を対象とした大規模導体表面電荷解析を実行し，2次導体の表面電荷分布を初めて明らかにした。

第 3 章では，近年の電磁環境問題を背景にその重要性がさらに高まりつつある磁気シールドを対象として，磁気特性の非線形性を考慮可能な積分方程式法である磁気モーメント

法^(12)(13)の大規模高速化，および非線形反復解法における収束特性の安定化について述べる。

磁気シールド問題は，磁性体の厚さとモデルサイズとの間に高いスケール比を有し，かつ開領域・非線形問題となる。そのため，通常の有限要素法で高精度な解を得るためには膨大な要素数を必要とする。また，高スケール比の影響で要素が扁平になることは避けられず，反復解法の収束特性に難点が生じる可能性がある。一方，積分方程式法の一種である磁気モーメント法は，磁性体領域のみの要素分割でよく，磁気特性の非線形性を考慮可能であり，かつ無限領域を容易に扱うことができるため，磁気シールドを対象とした磁界解析において非常に有効な手法であると考えられる。しかし，非対称密行列を扱うため，膨大な計算時間と記憶容量を必要とする難点があった。また，磁気特性の非線形性が強い場合，ニュートン・ラフソン法のような非線形反復解法の収束特性が非常に悪化することも考えられる。そこで本章では，まず高速多重極法およびACAを磁気モーメント法に導入し，

非線形磁気シールド解析の大規模化を達成した。次に，多くの反復回数を要する非線形解法の収束特性をより安定化させることを目的に，直線探索⁽¹⁴⁾を積分方程式法に適用した。

その結果，非線形反復毎に適切な修正係数を求めることが可能となり，積分方程式法において非線形収束特性のロバスト性が向上することを明らかにした。行列圧縮アルゴリズム，

(9)

直線探索および磁気モーメント法の特長を活かした開発手法により，安定な収束特性を有する大規模かつ高速な非線形静磁界解析の実用化を達成した。

第 4 章では，解析対象を渦電流場問題にまで拡張し，無限領域や移動物体を容易に考慮可能な境界要素法の特長と有限要素法の高い汎用性を併せ持つ有限要素・境界要素併用法

(15)-(18)の大規模高速化について述べる。本手法は，空気領域の要素分割が不要で，電磁現象

特有の様々な困難さ（開領域，非線形・異方性・多媒質，移動境界）にも対応可能であり，

汎用性および解析モデル作成の観点から非常に有効な解析手法である。しかし，境界要素法が生成する密行列の影響により膨大な計算時間と記憶容量を必要とするため，大規模問題への適用は困難であった。また，有限要素・境界要素併用法の全体マトリクスは異なる性質を有する小行列の集合であり，反復解法の収束特性に難点があった。したがって，さらなる高速化のためには，適切な前処理による収束特性の改善が必須である。そこで本章では，種々の有限要素・境界要素併用法の定式化に対して境界要素法部分に高速多重極法を導入し，渦電流問題・非線形問題を対象とした電磁界解析の大規模化を達成した。また，

さらなる高速化を目指し，高速多重極法の特長を活かした反復前処理⁽¹⁹⁾と不完全LU分解を組み合わせた新たな前処理方法を提案し，その有効性を具体的な数値解析例に基づき確認した。さらに，磁気ベクトルポテンシャルと磁界の強さを変数とする有限要素・境界要素併用法において，磁界の接線成分を物理的に等価な磁気スカラポテンシャルに変換することで高速多重極法の演算量を削減する高速化アルゴリズムの開発も行った。

第 5 章では，上記開発手法を用いて電気機器解析に残された重要課題の一つである積層鉄芯の磁界解析を試みた。電気機器設計においては，その損失を低減し高効率化を図るため，高精度な積層鉄芯の電磁界解析が要望されている。一般に有限要素法のような領域法のみでは，周囲の空気領域を考慮した上で積層構造を直接メッシュ分割した解析を実行することは，ギャップ要素⁽²⁰⁾・2重節点⁽²¹⁾などの工夫を施しても膨大な要素数を必要とするため，記憶容量および計算時間の観点から実用上困難である。その膨大な計算コストを削減するため，均質化法^(22)-(24)などの積層鉄芯近似モデル化手法が同じく有限要素法ベースで提案されているが，実機解析における計算コストや精度などについては，不明瞭な部分も残されている。そこで，代表的な積層鉄芯近似モデル化手法を永久磁石同期電動機の鉄損解析に適用し，均質化法の実用性をまず明らかにする。一方，本研究で開発した有限要素・

境界要素併用法は，有限要素法の単一利用とは異なり空気領域の要素分割が必要ないため，

解析精度を損なうことなく要素数を大幅に削減することができる。そのため，電磁鋼板および絶縁ギャップを多数層に分割し，積層構造を詳細に模擬することが可能である。そこで，第 4 章で開発した高速多重極法を導入した有限要素・境界要素併用法により，実機レベルと同等の積層枚数・アスペクト比を有する積層鉄芯の厳密な 3 次元解析を実施し，各種の積層鉄芯近似モデル化手法を比較するための基準となる解を導出した。このような積層鉄芯の詳細な 3 次元解析は過去に例がない。得られた基準解との比較により，磁束密度の変化が非常に急峻な鉄芯端部においても，有限要素法ベースの均質化法を適用した際に

(10)

見られる精度劣化を生じることなく，開発手法では高精度な解析が実現可能であることを明らかにした。

第6章では，本研究で得られた知見を総括し，各章ごとに成果の要約を示す。

(11)

参考文献

(1) 電気学会：「電磁界解析における高速大規模数値計算技術」，技術報告，no. 1043 (2006).

(2) 高橋則雄：「三次元有限要素法 ― 磁界解析技術の基礎」，電気学会 (2006).

(3) 加川幸雄，榎園正人，武田毅：「電気・電子境界要素法基礎と応用」，森北出版 (2001)．

(4) L. Greengard and V. Rokhlin, “A New Version of the Fast Multipole Method for the Laplace Equation in Three Dimensions,” Acta Numerica, vol. 6, pp. 229-269 (1997).

(5) H. Cheng, L. Greengard, and V. Rokhlin, “A Fast Adaptive Multipole Algorithm in Three Dimensions,” J. Comput. Phys., vol.155, pp.468-498 (1999).

(6) S. Kurtz, O. Rain, and S. Rjasanow, “The Adaptive Cross-Approximation Technique for the 3-D Boundary-Element Method,” IEEE Trans. Magn., vol. 38, no. 2, pp. 421-424 (2002).

(7) M. Bebendorf and S. Rjasanow, “Adaptive Low-Rank Approximation of Collocation Matrices,” Computing, vol. 70, pp. 1-24 (2003).

(8) 小林昭一編：「波動解析と境界要素法」, 京都大学学術出版会 (2000).

(9) 宅間董，濱田昌司：「数値電界計算の基礎と応用」, 東京電機大学出版局 (2006).

(10) A. Buchau, W. M. Rucker, O. Rain, V. Rischmuller, S. Kurz, and S. Rjasanow, “Comparison between different approaches for fast and efficient 3-D BEM computations,” IEEE Trans.

Magn., vol. 39, no. 3, pp.1107-1110 (2003).

(11) Y. Fujishima and S. Wakao, “Surface Charge Analysis in Eddy Current Problems,” IEEE Trans. Magn., vol. 39, no. 3, pp. 1123-1126 (2003).

(12) 野島洋一，大場章人：「磁気モーメント法による三次元磁界解析の高速化」，日本シミ

ュレーション学会第9回計算電気・電子工学シンポジウム論文集，vol. 9, pp. 113-118 (1988).

(13) 薮内利親，朝井敏文，矢野博幸：「磁気モーメント法と表面磁荷法の長所を併せ持つ

多面体要素の開発」，電気学会静止器・回転機合同研究会資料，SA-05-05，RM-05-05 (2005).

(14) K. Fujiwara, Y. Okamoto, A. Kameari, and A. Ahagon, “The Newton-Raphson Method Accelerated by Using a Line Search – Comparison between Energy Functional and Residual Minimization,” IEEE Trans. Magn., vol. 41, no. 5, pp.1724-1717 (2005).

(15) S. Wakao and T. Onuki, “Novel Boundary Element Formulation in Hybrid FE-BE Method for Electromagnetic Field Computation,” IEEE Trans. Magn., vol. 28, no. 2, pp. 1162-1165 (1992).

(16) S. Wakao, T. Onuki, and M. Shimazaki, “3D eddy current analysis by the hybrid FE-BE method using magnetic field intensity H,” IEEE Trans. Magn., vol. 28, no. 5, pp. 2259-2261 (1992).

(17) S. Wakao and T. Onuki, “Electromagnetic field computations by the hybrid FE-BE method using edge elements,” IEEE Trans. Magn., vol. 29, no. 2, pp. 1487-1490 (1993).

(12)

(18) T. Onuki, S. Wakao, and T. Yoshizawa, “Eddy Current Computations in Moving Conductors by the Hybrid FE-BE Method,” IEEE Trans. Magn., vol. 31, no.3, pp.1436-1439 (1995).

(19) S. Hamada and T. Takuma, “Effective Precondition Technique to Solve a Full Liner System for the Fast Multipole Method,” IEEE Trans. Magn., vol. 39, no. 3, pp.1666-1669 (2003).

(20) T. Nakata, N. Takahashi, K. Fujiwara, and Y. Shiraki, “3-D Magnetic Field Analysis Using Special Elements,” IEEE Trans. Mag., vol. 26, no. 5, pp. 2379-2381 (1990).

(21) Y. Kamiyama and T. Onuki, “3-D Eddy Current Analysis by the Finite Element Method Using Double Nodes Technique,” IEEE Trans. Mag., vol. 32, no. 3, pp. 741-744 (1996).

(22) H. Kaimori, A. Kameari, and K. Fujiwara, “FEM Computation of Magnetic Field and Iron Loss in Laminated Iron Core Using Homogenization Method,” IEEE Trans. Magn., vol. 43, no. 4, pp. 1405-1408 (2007).

(23) K. Muramatsu, T. Okitsu, H. Fujitsu, and F. Shimanoe, “Method of Nonlinear Magnetic Field Analysis Taking Into Account Eddy Current in Laminated Core,” IEEE Trans. Magn., vol. 40, no. 2, pp. 896-899 (2004).

(24) J. Gyselinck, R. V. Sabariego, and P. Dular, “A Nonlinear Time-Domain Homogenization Technique for Laminated Iron Cores in Three-Dimensional Finite Element Models,” IEEE Trans. Magn., vol. 42, no. 4, pp. 763-766 (2006).

(13)

2.1 緒言

電気機器を設計する上で有用な情報を高精度に抽出する一手法として，渦電流場における導体表面電荷の大規模解析について述べる。導体表面電荷は導体中に発生する渦電流分布に大きく依存するため，例えば機器内における導体形状の変化が機器性能に及ぼす影響の評価をはじめ，様々な機器特性の情報抽出を行う上で有効に活用できる⁽¹⁾。この導体表面電荷は，準定常場における辺要素有限要素法の解析結果に後処理を施し，クーロンゲージを付加した導体表面上の電気スカラポテンシャルから，積分方程式法を用いて算出される。

そのため，一般的に離散化後に得られるマトリクスが密行列となり，膨大な計算時間と記憶容量が必要となる。これが，回転機のような複雑形状へ導体表面電荷解析を適用する際の大きな障壁となっていた。表面電荷解析の適用範囲を複雑な形状を有する実機設計への活用にまで引き上げるためには，高速化・省メモリ化が必須となる。

そこで本章では，実規模問題に適用可能な導体表面電荷解析技術の確立を目的とした種々の検討を行った。まず，代表的行列圧縮アルゴリズムである高速多重極法^(2)-(5)および

ACA^(6)-(7)について概説し，これらの大規模数値解析手法を組み合わせた表面電荷法を開発し

た。その際，高速多重極法および ACA の比較を行い，両者の得失を明らかにした。また，

誘導電動機などの電気機器はほとんどの場合，対称性もしくは周期性を有しているため，

この対称性・周期性を活用して計算時間の削減を図ることは極めて重要となる。そこで，

さらなる高速化を達成するため，高速多重極法において解析モデルの対称性・周期性を活用する方法を新たに提案し，その有効性を検証した。最終的に，複雑形状を有する実規模モデルとして，かご形誘導電動機を対象とした大規模導体表面電荷解析を実施し，二次導体の表面電荷分布を明らかにした⁽⁸⁾。

2.2 高速多重極法

2.2.1 高速多重極法の概略計算手順

本論文では，ラプラス場を扱う。ラプラス場においては，任意の点のポテンシャルを多重極展開と局所展開を用いて，無限級数で記述することができる。ツリー法⁽⁹⁾や高速多重極法は，点電荷の集合にセル構造を導入し，その構造を利用して分割統治計算を進める。その際，無限級数で表されたポテンシャルの多重極展開係数と局所展開係数を効率的に計算することにより，必要な精度を確保しつつ多体間相互作用の近似値を高速に計算する。以下，高速多重極法の計算手順について，概要を述べる。

(14)

まず，点電荷の集合にセル構造を導入する。図2.1 に，2次元空間に 15個の点電荷が分布している場合の例を示す⁽⁴⁾。境界要素法などに高速多重極法を導入する際には，三角形や四角形など要素の集合に対してセル構造を作成するが，要素重心を点電荷とみなせば同様に考えることができる。しかし，境界要素は点電荷とは異なり，有限の大きさを有する。

そこで本論文では，解析領域中の境界要素に対し，以下のようにセル構造を作成した⁽¹⁰⁾。 1. 全重心点群を内包する立方体を考え，その立方体と中心が同一で1辺の長さをC₁倍に拡

大した立方体をルートセルとする。

2. 注目しているセルを同じ大きさの8個の立方体領域に分割し，8個の分割領域のうち重心点を内包する領域にのみ新規セルを定義する。このとき，分割するセルを親セル，親セルの分割により新たに定義されたセルを子セルとする。

3. セル分割に際して，新規セル中心から新規セル内包境界要素の各頂点までの全距離を調べ，新規セルの多重極展開半径の1/C₂倍以上のものがない場合のみ新規分割を行う。存在する場合は（境界要素が新規セルから大きくはみ出すならば），新規分割は取りやめ，

もとのセルをリーフセルとする。

4. セル分割に際して，新規セルの点群数がすべてC₃個以下ならば，新規分割は取りやめて，

元のセルをリーフセルとする。

ただし，C₁，C₂，C₃は1以上の定数である。

セル間の遠方と近傍は，セル構造を基準に定義される。2つのセル間に1つ（第1近傍）

または2つ（第2近傍）のセルが挟まれている状態を，2つのセルがお互いに遠方にあるという。図2.2に，2つのレベルのセルについて，第1近傍を基準としてその遠方と近傍のセルの関係を示す⁽⁴⁾。注目しているセルから見ると，親セルの遠方は遠方にあり，親セル近傍の中にも子セルと同レベルの遠方セルが存在する。親セルの近傍セルの中にリーフセルがあるときには，そのセルは子セルにとっても近傍となり，その類別は孫にも引き継がれる。

得られたセル構造に基づいて，できるだけ多くの点電荷からの影響を多重極・局所展開表現を用いて一括計算する。以下に，計算手順を簡単に述べる。

1. リーフセル中心に多重極展開を定義する（MP）。

∑

=

−

=

^k + i

i i m An n A

n i i m

An

Y

q R M

1

] 1 [ ] [ ]

[

ρ ( α , β )

(2.1)

2. 子セルで定義された多重極展開を親セルに移動して加え合わせ，これをルートセルまで繰り返す（M2M）。

∑ ∑

= =− − +

+

− +

+

− +

−+

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

= ^j

n n

n

m j n

A k j m k m k

A m k

n B j

m n j

B m k

n j m n B

k j

A R i

M R Y

A A M

0 1

] [ ]

[ 1

] [

) , ( ρ

β ρ α

(2.2)

(15)

1 2 4 3

6 5 7

8

10 11

12 13 14 15

Level 0

Level 1

Level 2

Level 3 9 8 12 13 14 15 6 5 4 3

10 11 7 1 2

leaf empty cell

9

図2.1 2次元空間におけるセル構造

Parent cell of Near cells of Far cells of Far cells of Parent cell of Near cells of Far cells of Far cells of

図2.2 遠方・近傍セルの定義

(16)

3. 遠方セルの多重極展開を自分のセルの局所展開に変換して加え合わせる（M2L）。

∑ ∑

⁺

=

− +

−

= − +

+

− +

+

−

−+

−

⎟⎟⎠

⎜⎜ ⎞

⎝

⎟⎟ ⎛

⎠

⎜⎜ ⎞

⎝

⎛

= ^j ^p −

j n

j n k

m n j

A j

B m n k m m k

k m

j An m

An k m

j n k j j n k

B j

R A R

i

M Y

A L A

) (

)

( 1

] [ ]

[ ]

[

) , ( )

1 (

ρ ρ

β

α

(2.3)

4. 親セルに局所展開がある場合には，その局所展開を自分に移動し，加え合わせる（L2L）。

∑ ∑

⁻

= =− +

+ + + +

−

+ + +

+

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

⎟⎟⎠

⎜⎜ ⎞

⎝

− ⎛

=^p ^j

n n

n

m j

B m k

n j k m m k

A m k

n B j

m n n j

A m n k j n j

B k

j

A R i

L R Y

A A L

0

] [ ]

[ ]

[

) , ( )

1 (

ρ β ρ α

(2.4)

5. リーフセルにおいて，(2.5)式より求める局所展開の寄与（LC）と，直接計算により求める近傍セル内の点電荷の影響を加え合わせる。

∑ ∑

= =−

⎟⎟

⎠

⎜⎜ ⎞

⎝

=

^p

⎛

n

m An n

A n

n m

m

An

Y

R L r

0 [ ] [ ]

4 1

φ π

(2.5)

ここで，M_n^m，L_n^mは多重極・局所展開係数，Y_n^mは文献(2)で定義される球面調和関数，mは次数nに対して-n ≤ m ≤ nを満たす整数，iは虚数単位，pは展開係数の打ち切り上限数，A_n^m はmとnにより与えられる定数⁽²⁾，添字A, Bは変換前後の極中心位置，(ρ, α, β)はBからA を観測したときの極座標，R_A，R_Bは極A およびBに対応する球面の陽に与えた展開半径である⁽¹¹⁾。図2.3にM2M，M2L，L2L変換公式の概念図を，図2.4に高速多重極法の概略計算手順の概念図を示す⁽⁴⁾。また，(2.3)式に対応するM2L変換マトリクスの概念図を図2.8(a) に示す。

通常の高速多重極法では，(p+1)²個の展開係数から(p+1)²個の展開係数への変換を行うた

めO(p⁴)の演算量が必要となり，次数の増加に伴い計算時間が急激に増大する。

2.2.2 擬似粒子法

高速多重極法における各変換公式の複雑さは，定式化に球面調和関数が含まれていることに起因する。そこで，球面調和関数を用いないAndersonの方法⁽¹²⁾や擬似粒子法⁽¹³⁾などの手法も提案されている。擬似粒子法では，多重極・局所展開係数を場の展開半径に対応する球面上に配置された擬似点電荷群により表現する。その際，粒子の空間位置は固定して電荷の大きさだけ変えるものとし，擬似粒子の空間配置としてt-design点⁽¹⁴⁾を採用する。擬似粒子法は，変換公式が単純で種類も少ないため，通常の高速多重極法と比較してコーディングが容易となる。

t-design点は単位球面S上の座標点r_i (1 ≤ i ≤ K)の集合であり，t次以下の任意の多項式g(r) に対して次式が成り立つ。

(17)

A( ρ , α , β ) R

_A

R

_B

B

] [A

M

]

M

[B

(a) M2M

A( ρ , α , β )

B R

_B

R

_A

] [A

M L

[B]

(b) M2L

A( ρ , α , β ) R

_B

R

_A

B

] [A

L

]

L

[B

(c) L2L

図2.3 展開係数変換公式の概念図

(18)

MP M2M

M2L

L2L

LC MP

M2L

Direct Direct

Calculation Calculation MP

M2M

M2L

L2L

LC MP

M2L

Direct Direct

Calculation Calculation

図2.4 高速多重極アルゴリズムの概念図

(19)

∫ ∑

=

= ^K

i

S g i

ds K g

1

) 4 (

)

(r π r _(2.6)

(2.6)式を数値面積分公式とみなすと，t-design点は重みがすべて4π/Kである積分点群を意味する。擬似粒子法における各種変換公式の無限級数は，採用したt-design点に応じてt /2次以下で打ち切るとよい⁽¹²⁾。

図2.5に，擬似粒子法による多重極・局所展開の概念図を示す。(2.2)式に対応するM2M 変換では，子セルでの擬似点電荷qを，(2.7)式を用いて親セルでの擬似点電荷q’_jに変換する（図2.5 (a)）。リーフセルにおいて多重極展開を求める際も，(2.7)式を用いて点電荷qの寄与を擬似点電荷q’_jに変換する。また，(2.3)式，(2.4)式に対応するM2L，L2L変換では，

遠方セルもしくは親セルにおける擬似点電荷qを，(2.8)式を用いて着目しているセルでの擬似点電荷q’_jに変換する（図2.5(b)）。

∑

^∞

=

⎟⎠

⎜ ⎞

⎝ + ⎛

′ =

0

) (cos )

1 2 1 (

n

j n n

j P

q a K n

q

ρ γ

_(2.7)

∑

^∞

=

+

⎟⎟⎠

⎜⎜ ⎞

⎝ + ⎛

=

0

1

) (cos )

1 2 1 ( '

n

j n n

j a P

q K n

q γ

ρ ^(2.8)

ここで，γ _jは擬似点電荷qと擬似点電荷q’_jがなす角度，aは展開半径，ρ は極中心Cと擬似点電荷qの距離，Pnはルジャンドル多項式，nはその次数である。任意のP(r,θ,φ)，r<a でのポテンシャルφ は，次式で求められる（図2.5(c)）。

∑ ∑

=

∞

=

+

⎟⎠

⎜ ⎞

⎝

′ ⎛

= ^K

j n n j

n

j P

a r r q

1 0

1

) 4 (cos

1 γ

φ π

(2.9)

ここで，γ _jは擬似点電荷q’_jと点Pがなす角度である。

また，境界要素法においては，電気双極子の作る電位と同形式（2重層ポテンシャル）の多重極展開が必要となる。球面の中心を基準とした擬似粒子の位置ベクトルを a_k，その強さをq’_k，双極子をp = qd = qdn_d，その位置ベクトルをρ ，擬似粒子と双極子のなす角をγ _k とすると，最終的に双極子の擬似粒子展開式は次式のように求められる（図2.5(d)）⁽¹⁵⁾。

⎭⎬

⎫

⎩⎨

⎧ ′ − ′

⎟ ⋅

⎠

⎜ ⎞

⎝ + ⎛

′ = ^∞ ₋

∑

=

⁽ ² ¹ ⁾ ^(cos ⁾

¹

^(cos ⁾

1

k n

n

k P P

a n a

K

q p

γ

γ ρ ρ

ρ

n_d a_k ρ

(2.10)

擬似点電荷群に変換された後は，電気双極子場も点電荷と同様に扱うことができる。

擬似粒子法は演算量がO(p⁴)のため， pが小さい場合には有用であるものの，次数の増加にともなう計算時間の大幅な増加は避けられない。

(20)

C q a

ρ γ

_j

q'j

C a

ρ

γ

_j

q

q'j

C a

ρ

γ

_j

q

q'j

(a) M2M (b) M2L, L2L

C P a

r γ

_j

q'j

C a

ρ γ

_j

p n

_d

q'j_j

q'

(c) Evaluation of potential (d) Multipole expansion of dipole 図2.5 擬似粒子法の概念図^(10)(15)

(21)

2.2.3 回転変換版高速多重極法

この演算量を削減する手法として，M2M，M2L および L2L 変換に回転変換を適用する

(2)(3)(5)(16)。球面調和関数の定義軸を図 2.6 のように回転させることにより，Ynm(0,0)=0，

Yn0(0,0)=1 となることから，(2.2)式から(2.4)式で表される変換公式を(2.11)式から(2.13)式のように単純化することができる。

] [

0 1

1 0

]

[ A

k n j j

n j n

A k j

j

B k

n j n B

k

j

M

A R A R A

M

₋

= − +

+

∑

−

⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎛

⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎛

= ρ

ρ

(2.11)

] 1 [ 0

] [

) 1 (

A k

j n p

j

n j n

A j

B n

k n j k j n j k B

k

j

M

R A R

A

L

⁺

A

₋

= − + +

+

− +

∑

−

⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎟⎟ ⎛

⎠

⎜⎜ ⎞

⎝

⎛

= −

ρ

^(2.12)

] [ 0

0

] [

) 1 (

A k

n j p

j

n j

B k

n j

A k j n n

B k

j

L

A R A R A

L

⁻ ⁺ ₊

=

+

∑

⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎛

⎟⎟ ⎠

⎜⎜ ⎞

⎝

− ⎛

= ρ

ρ

(2.13)

展開係数の変換前後に行う座標系の回転・逆回転には，それぞれO(p³)の演算量が必要であり，全体としてもO(p³)の演算量となる。回転操作の詳細については，文献(2)(3)(5)(16)を参照されたい。(2.12)式に対応するM2L変換マトリクスの概念図を図2.8(b)に示す。

2.2.4 対角変換版高速多重極法

対角変換版高速多重極法(2)(3)(5)(16)では，回転変換版高速多重極法の中で最も演算量を要するM2L変換において指数関数展開を導入することにより，さらなる高速化を図る。対角変換，回転変換の概念図を図2.7に示す。多重極展開係数Mを指数関数展開係数Wに変換した後，指数関数変換で極移動（W からV）を実行し，最終的に指数関数展開係数Vを局所展開係数Lに逆変換する。これらの変換式を(2.14)式から(2.16)式に示す。

n

e k M e M p

m n

nm p

p m

m im j

k R

R R

R m n m n e M

i

W ^j ⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛ +

− −

=

∑ ∑

=

−

=

− ^α

λ

)!

( )!

) (

(

_(2.14)

⎭⎬

⎫

⎩⎨

⎧ − − −

− −

=

^e ^j

M j L e

M k L e

M

k L R

y y R

x i x R

z z j k j

k

W e e

V

^λ ^λ ^cos^α ^sin^α ^(2.15)

(22)

∑ ∑

= =

−

⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎛ − +

= −

⁽ ⁾

1 1

) )! (

( )!

(

1

^s^ε

λ ω

α

k

M

j

j k m im

k k n

e k L mn

k

j

V e M i

R R m

n m

L n

(2.16)

ここで，Reはセル辺長，(xM, yM, zM), (xL, yL, zL)は多重極・局所展開の極位置である。また，

s(ε ), λ_k, ω_k, Μ_k, α_jは，指数関数変換を行う際の数値積分公式に必要な変数である⁽²⁾⁽³⁾。対角変換におけるM2Lでは，多重極・局所展開の指数関数展開への変換・逆変換にそれぞれO(p³)，

(2.16)式の変換にO(p²)の演算量を必要とする。また，指数関数変換を行うためには，座標系の回転変換が必要であり，回転・逆回転操作にそれぞれO(p³)の演算量が必要となる。図2.8(c) に，(2.16)式に対応する変換マトリクスの概念図を示す。

2.2.5 ^{対称性の活用方法}

対称性を考慮する領域には，実際の解析領域と全く同じセル構造を対称面・軸・点に対して鏡像を取ることにより定義する。対称領域からの寄与は，M2L 変換時に必要となる。

そこで，対称領域中の多重極展開係数Mˆ_n^mを，実解析領域中の対応するセルで定義された多重極展開M_n^mを用いて表すことを考える。Mˆ_n^mは，実領域のセル中心に定義されたM_n^mと，対称面・軸・点に対して反転した座標位置に定義されている。そのため，Mˆ_n^mは球面調和関数の定義式より，M_n^mを用いて次のように表すことができる⁽¹⁷⁾。

1. xz平面に関して対称： Mˆ_n^m=M_n^m 2. z軸に関して対称： Mˆ_n^m=(−1)^mM_n^m 3. yz平面に関して対称： Mˆ_n^m=(−1)^mM_n^m 4. xy平面に関して対称： Mˆ_n^m=(−1)^m⁺ⁿM_n^m 5. x軸に関して対称： Mˆ_n^m=(−1)^m⁺ⁿM_n^m 6.原点に関して対称： Mˆ_n^m =(−1)ⁿM_n^m 7. y軸に関して対称： Mˆ_n^m =(−1)ⁿM_n^m

ここで，_M_n^mは_M_n^mの複素共役を表す。図2.9に，本手法の概念図を示す。

解析モデルにおいて，電荷密度などの物理量分布が対称であれば展開係数は等しくなり，

反対称であれば符号を反転すればよい。よって，対称領域における多重極展開係数Mˆ_n^mは，

特別な演算をすることなく容易に求めることができる。得られたMˆ_n^mを用いて回転変換版 M2Lを実行すれば，対称領域からの寄与分を計算できる。

一方，擬似粒子法においては，対称領域における擬似粒子を実空間の擬似粒子の対称面・

軸・点に関して鏡像をとることにより定義すれば，対称領域における擬似粒子と実粒子の成す角は実空間のものと等しくなる。したがって，対称領域の擬似粒子の電荷の強さの絶対値は，実空間のそれと一致する。対称領域からのM2L計算の際にも実空間と同じ変換行列を用いるため，M2Mと同様の変換を行うことで対称領域の擬似粒子を実空間の相対位置に移動した後，M2Lを実行する⁽⁸⁾。

(23)

Rotation

α ^＝ β ^＝ 0 ( ρ , )

x

y

z z’

x’

y’

0 0 , ) , ( ρ

x

y

z z’

x’

y’

0 0 ,

x

β α

y z

) , , ( ρ α α , β β

x

β α

y z

) , , ( ρ α α , β β

図2.6 回転変換における座標系の回転

O(p

⁴

)

M L

M ˆ L ˆ

W V

O(p

³

)

O(p

²

) Rotation:

O(p

³

) Rotation:

O(p

³

)

Exponential translations Plane wave

expansion: O(p

³

) Plane wave

expansion: O(p

³

)

図2.7 対角変換，回転変換の概念図

(24)

* * * * * * * * * * * * * * * * * * * * * * * * *

L M

(a) M2L translation matrix

* * * * *

* * * *

* * * * *

* * * *

* * *

* * * * *

* * * *

* * *

* *

* * * * *

* * * *

* * *

* *

*

M ˆ

L ˆ

(b) M2L Translation matrix after rotating coordinate system

*

V W

(c) Diagonal translation matrix

図2.8 M2L変換マトリクスの概念図

(25)

擬似粒子法や回転変換版高速多重極法は，遠方セルに定義されたMを自分のセルのLへ変換するのに対し，対角変換版高速多重極法では自分のセルに定義されたMを遠方セルの L に変換する。対角変換版高速多重極法に関しても，上記と同様のアルゴリズムでMˆ_n^mを求め，対称領域からのM2Lを対角変換で考慮することは可能であるが，それほど効率的ではないと考えられる。よって，対角変換版高速多重極法においては，対称領域からのM2L変換は回転変換で行うことにする。

2.2.6 高速多重極法の積分方程式法への適用と反復前処理

積分方程式法における面積分・体積積分を，数値面積分公式に応じた積分点に離散化すると，要素上に配置された積分点個の等価な点電荷に置き換えることができる。したがって，前述の点電荷に対する高速多重極法がそのまま適用可能となる。図 2.10に，積分点に応じた点電荷に置き換えた場合の概念図を示す⁽¹⁸⁾。

積分方程式法によって得られる連立一次方程式Ax = bは，1つの積分方程式（＝1行分の式）にはすべての未知変数が寄与し，係数行列AはN×Nの密行列になる。ここで，「すべての未知変数が寄与する」というのは，N体間相互作用の総当り計算と同じ状況である。よって，積分方程式法に適用された高速多重極法は，係数行列Aと任意ベクトルx との積を高速計算する手法と等価である。そこで，連立方程式の解法として反復法を適用する。

高速多重極法を導入した積分方程式法の前処理方法として，反復前処理⁽¹⁹⁾を採用する。

以下に，その概要を述べる。反復解法を用いた連立方程式の求解には，係数行列 A と任意ベクトル b_arbの積Ab_arbの計算が必要である。その収束特性を高めるため，通常は近似行列 A_apprを用いた前処理A_appr^-1b_arbを行う。反復前処理では，A_apprx_arb= b_arbに反復法を適用し，

任意ベクトルx_arb= Aappr-1b_arbを求めることで前処理を実行する。つまり，Ax = bの解を求めるループ（主反復）の中に前処理部分であるx_arb= Aappr-1b_arbを求めるループ（副反復）が存在し，入れ子状の構造となる。副反復は前処理用であるから求める値は近似値でよい。そのため，高速多重極法の計算精度を故意に低下させることで高速化を図ることができる。

このとき，A_apprは遠方からの寄与も含めた A の大域的な近似行列となり，良好な収束特性が期待できる。

密連立一次方程式の反復解法としては，GMRES法⁽²⁰⁾やBi-CGSTAB法⁽²⁰⁾が一般的に良く用いられている。GMRES 法は 1 反復あたりの行列ベクトル積の回数が 1 回だけであり，

Bi-CG系統の反復法の2回に比べ計算時間の観点から優位だと考えられる。しかし，GMRES

法では反復ごとに生成されるクリロフ部分空間基底ベクトルを破棄せずメモリに格納する必要があり，反復回数が増えるたびに計算量およびメモリ量が急激に増加する。反復前処理を導入した反復法では，主反復はたかだか数回の反復で収束するが，副反復内では多くの反復回数が必要とされるため，GMRES法は副反復内の解法としては適さないと考えられる。以上の理由から，本論文においては高速多重極法の場合，主反復用の反復解法には GMRES法を，副反復用の反復解法にはBi-CGSTAB2法⁽²⁰⁾を採用した。

(26)

Inverting with respect to a symmetric plane

Real space Image space

Direct calculation

M2L

図2.9 対称性の考慮方法

Discretization with numerical integral formula

Boundary element

Discretization with numerical integral formula

Boundary element

図2.10 高速多重極法の積分方程式法への適用例⁽¹⁸⁾

積 分 方 程 式 法 に よ る 大 規 模 電 磁 界 数 値 解 析 の 実 用 化 に 関 す る 研 究