第 4 章 扁平要素混在時の求解過程での収束性改善手法

扁平要素混在時の求解過程での収束性改善手法

4.1 緒言

本章では，代数マルチグリッド法を前処理とした共役勾配法[4.1][4.2]による扁平要素混在 時の収束性改善手法について述べる．

電気機器の高効率化を目的とした電磁界数値解析では要求される計算精度が非常に高い ため，機器の形状を忠実にメッシュモデルで再現する必要があり，扁平な部材も簡略化せ ずに扱わざるを得ない．特に電気機器の鉄損軽減のために一般に用いられる積層鉄心にお いて，その内部の磁束密度分布を正確に得るには積層鋼板を 1 枚ずつメッシュモデル化す る必要があり，そのメッシュは極端に高アスペクト比の扁平要素を含みやすい．扁平要素 は数値解析で解くべき連立一次方程式の性質を悪化させるため，方程式の反復解法の収束 性が悪化し，結果として計算時間の大幅な増加を招く場合がある．電磁界数値解析では一 般に不完全コレスキー分解を前処理とした共役勾配法（ICCG法）[4.3]が反復解法として用 いられるが，扁平要素を多く含むメッシュモデルでは，著しく収束性が悪化する傾向にあ る．許容される計算時間の制約を満たす要素数の中で，機器周辺の空間も考慮した上で全 体的に扁平な要素の度合いを極力少なくしたメッシュを生成することは，一般技術者には 容易ではない．

一方，近年，大規模メッシュモデルに対する高速解法としてマルチグリッド法が注目さ れている．マルチグリッド法の計算時間は理想的には計算自由度数に比例するため，計算 時間が計算自由度数の1.5乗程度に比例すると言われるICCG法に比べて高速化が期待でき る[4.4]．

そこで本章では，積層構造を有する電気機器の高速解法を目指し，反復解法の前処理と して代数マルチグリッド法を用いることにより ICCG 法に比べて計算時間の悪化を緩和で きることを見出し，またその収束性に大きく影響する制御パラメータを明らかにした．

4.2   電磁界解析における代数マルチグリッド法

4.2.1 代数マルチグリッド法

マルチグリッド法は，ある基礎方程式を有限要素法で離散化して得られる次式の連立一 次方程式の反復解法として用いられる．

[ ]

K

{ } { }

A = b ^(4.1)

ここで，[K]は有限要素法における連立一次方程式の係数行列，{A}は解ベクトル，{b}は右 辺ベクトルである．マルチグリッド法は，階層化された複数の係数行列を用いて効率よく 誤差を除去する手法であり，その係数行列および階層間で情報を射影する際に用いる行列

（補間行列）を代数的に作成する代数マルチグリッド法（以下，AMG法と呼ぶ）[4.5][4.21]

とメッシュモデルの情報を用いて作成する幾何マルチグリッド法（以下，GMG法と呼ぶ）

[4.6][4.7]に大別される．GMG法で用いる階層化された係数行列は粗密の異なるメッシュモ デルから生成するため，ユーザーは予め複数のメッシュモデルを作成しておかなければい けない．また，GMG法における反復解法の収束性は各階層用に用意したメッシュモデルに 依存するため，その生成ノウハウが求められる．ただし，その階層化された係数行列は物 理的にも明確であり，また形状関数に基づいて作成される補間行列は各階層間で受け渡さ れる残差および摂動量をメッシュ上での連続性を損なうことなく補間できるため，経験的 にAMG法に比べて安定的に収束解が得られることが多い．一方，AMG法はGMG法と比 べて経験的に反復解法の収束性が劣る場合が多い．また，電磁界解析特有の辺要素有限要 素法に適用する場合，AMG法の適用例が多い構造解析や流体解析における節点要素有限要 素法とは異なる特殊な工夫が必要となる．ただし，AMG法はGMG法のように予め複数の 粗密メッシュモデルをユーザーが用意する必要がなく，また既存のソルバープログラムの ICCG法のサブルーチンをAMG法のサブルーチンに置き換えるだけで開発することができ る．その利便性は，ユーザーのメッシュ生成負荷の軽減という点で非常に重要となる．そ こで本章では，AMG法について検討を行う．

AMG法で代表的なVサイクル処理手順を図4.1に示す．ここで，{A}は求めるべき解ベ クトル，{b}は右辺ベクトル，[K]は各階層での係数行列，{r}は各階層での残差ベクトル，

{e}は各階層での摂動量ベクトルであり，右下の数字は階層の深さを示すインデックスであ る．インデックス値が小さいほど下位階層であることを示し，下位階層になるほど係数行 列の規模は小さくなる．まず，最上層の連立一次方程式についてヤコビ法やガウス・ザイ デル法などの定常反復法を数回だけ反復させて暫定的な解と残差を得る．この処理をスム ージングと呼び，スムージングで用いる反復解法をスムーザーと呼ぶ．上位階層で得られ た残差は下位階層へと射影され，その射影残差を右辺ベクトルとして下位階層で同様にス ムージングを行い，その下位階層における解と残差を得る．この上位階層から下位階層へ

の残差の射影を Restriction と呼び，その射影で用いる Restriction 行列を[R]と表記する．

Restrictionとスムージングは最下層に到達するまで繰り返され，最下層では十分に解くべき 連立一次方程式のサイズが小さくなるため，LU分解[4.8]などの直接法で求解する．下位階 層で得られた解は上位階層に射影することで上位階層での摂動量となる．この下位階層か ら上位階層への摂動量の射影をProlongationと呼び，その射影で用いるProlongation行列を [P]と表記する．最下層からProlongationとスムージングを最上層に達するまで繰り返すこと により，各階層における解，つまり摂動量が順次更新され，最上層に到達した時点で求め るべき解の誤差が除去されたこととなる．これは，スムーザーに用いられる定常反復法が 誤差の空間的な高調波成分を速やかに除去する性質を利用したものである．上位階層にお いて除去されにくい誤差の低周波成分は下位階層に射影することで高周波成分とみなすこ とができ，スムーザーで効率よく誤差を除去することが可能となる．なお，[R]および[P]は 一般に [R]=[P]^Tとすることが多い．

{ } [ ] ^{{ }}

n n n

n R r

r₋₁ = ₋₁

[ ]

K2

{ } { }

A2 = b2

[ ]

K1

{ } { }

e1 = r1

[ ]

K₀

{ } { }

e₀ = r₀

Fine LevelCoarse Level

: Smoothing

: Solve

: Restriction

: Prolongation { } { }

e_n = e_n +

[ ]

P_nⁿ−1

^{{ }}

e_n−1

{ }

e₀ =

[ ]

K₀ ⁻¹

{ }

r₀

{ } [ ] ^{{ }}

n n n

n R r

r₋₁ = ₋₁

[ ]

K2

{ } { }

A2 = b2

[ ]

K1

{ } { }

e1 = r1

[ ]

K₀

{ } { }

e₀ = r₀

Fine LevelCoarse Level

: Smoothing

: Solve

: Restriction

: Prolongation { } { }

e_n = e_n +

[ ]

P_nⁿ−1

^{{ }}

e_n−1

{ }

e₀ =

[ ]

K₀ ⁻¹

{ }

r₀

図4.1  Vサイクルマルチグリッド法

Agg1

Agg3

Agg2

A C

B

Agg1

Agg3

Agg2

1

2

3

4 5 6

7

Fine Level Coarse Level

Agg1

Agg3

Agg2

A C

B

Agg1

Agg3

Agg2

A C

B

Agg1

Agg3

Agg2

1

2

3

4 5 6

7

Agg1

Agg3

Agg2

1

2

3

4 5 6

7

Fine Level Coarse Level

図4.2  節点要素における領域分割

節点有限要素法における[P]の作成方法の 1 つに，領域分割による手法[4.9]がある．上位 階層の節点を図4.2に示すようなAggregateと呼ばれる複数の節点集合の領域に分割し，そ のAggregate1 つ1つを下位階層の節点とみなして，下位階層の節点値を上位階層の節点値 に補間するPを作成する．Aggregateを作成するにあたり，上位階層の節点をCとFの節点 グループに分類する．Cは下位階層の節点として残る節点グループ，Fは上位階層にのみ存 在する節点グループと定義する．上位階層の節点グループ全体をΩとすると，F≡Ω-Cが成 り立つ．Cの中でも節点i に値を補間する節点グループをCiと定義すると，[P]は次式を満 足するように作成される．

( ) ⎪⎩

⎪ ⎨

⎧

∈

= ∑ ∈

∈

F i if u w

C i if u

Pu

Ci

j

c j ij c i i

c (4.2)

ここで，ui

cは下位階層の節点iでの値，wijは下位階層の節点iの値に対する上位階層の節点 jの値を補間する際の重みである．

以下に，Ciの選出方法およびwijを決定方法について述べる．

AMGのスムーザーとして用いられる定常反復法における解ベクトルの更新は，次式のよ うに表すことができる．

{ } ^A

ⁿ⁺¹

= [ ] [ ] ^M

⁻¹

^N { } ^A

ⁿ

+ [ ] ^M

⁻¹

{ } ^f

^(4.3)

ここで，{Aⁿ}は定常反復法の反復 n 回目における解ベクトル，[M]および[N]は係数行列[K]

を分解することにより得られる係数行列であり，ガウス・ザイデル法の場合，[M]=[D]+[L]，

[N]=-[U]となる．なお，[U]は[K]の上三角行列，[D]は対角行列，[L]は下三角行列であり，

[K]=[U]+[D]+[L]が成り立つ．

{r}={b}-[K]{A}の{r}を残差ベクトル，[K]{e}={r}の{e}を摂動量ベクトルとすると，定常 反復法による摂動量の更新は，次式のように表すことができる．

{ } ^e

ⁿ⁺¹

⁼ ( [ ] [ ] [ ] ^{I − M}

⁻¹

^K ) { } ^e

^{n (4.4)}

ここで，{eⁿ}は定常反復法の反復n回目における摂動量ベクトル，[I]は単位行列である．定 常反復法の収束性が悪化した場合，摂動量はほとんど更新されず

{ } { } ^e

ⁿ⁺¹

^{≈ e}

^{n の状態と考え}

られる．この関係式を式(4.4)に代入すると，

[ ] K { } e

ⁿ

≈ { } 0

が成り立つことがわかる．この関 係式より，節点iにおける誤差eiは次式のように表される．

ii j i

j ij

i

k

e k e

∑

≠

−

≈

^(4.5)

ここで，k_ijは係数行列のi行j列の値である．つまり，定常反復法において収束性が悪化し ている状態では，eiはkij/kiiの値に応じてejの影響を受けると考えられ，kij/kiiの値が大きい ほど節点iと節点jは誤差の伝播という点で影響力の強い接続関係にあると言える．そこで，

強い接続関係にある節点iおよびjを次式Siのように定義し，Ci=Si∩Cとすることで接続の 強い節点から値を補間するような[P]を得ることを考える[4.10]．

( )

{

^ij _{k i} ^ik

}

i

j i k k

S ≡ ≠ − ≥ −

max

≠

: θ

_(4.6)

ここで，θはユーザーパラメータであり，Aggregateの大きさや[P]の補間精度を決定する重 要なパラメータとなる．θには，一般的に0.2程度の値が用いられる．

式(4.5)を次式のようにまとめることで，wijの値が明らかになる．

ii C j

j ij

ii C k

k ik

i

k

e k k

e k

e

^{i i}

∑

∑

∈

−

∉

−

≈

^(4.7)

⎪⎪

⎩

⎪⎪

⎨

⎧ ∈

=

∈

∀

∑

∑

∈

∈ otherwise

k e k

S j if e

e C j

i i

C n

jn C n

n jn

i i

j

i, (4.8)

例えば，図4.3に示すような節点グループがあり，節点の添え字が節点番号，黒丸がCの 節点，白丸がFの節点，太線が節点1から見てS1の関係にある節点とつながる要素辺とし た場合，式(4.5)は図4.3に示す通りとなり，w12，w16を読み解くことができる．

3

4

5 6

7 2

1

Node in C

Node in F

Strong connection Weak connection

( ) ( )

11

16 12

6 16 2 12 14 1

13 6

16 2 12

1

k

k k

e k e k k

e k e

k e k

e ⎭ ⎬ ⎫

⎩ ⎨

⎧ ⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎛

+

− + +

− +

−

−

≈

3

4

5 6

7 2

1

Node in C

Node in F

Strong connection Weak connection

( ) ( )

11

16 12

6 16 2 12 14 1

13 6

16 2 12

1

k

k k

e k e k k

e k e

k e k

e ⎭ ⎬ ⎫

⎩ ⎨

⎧ ⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎛

+

− + +

− +

−

−

≈

図4.3  節点要素におけるProlongationの重みの作成例

なお，係数行列の各行番号が節点番号に対応するため，各行の非ゼロ要素の列番号から節 点同士の接続情報を読み取ることができ，特にメッシュモデルの情報は必要ない．[P]を決 定することにより，L-1階層目の係数行列[KL-1]はL階層目の係数行列[KL]から次式のように 生成することができる．

[ ] ^K [ ] ^P [ ] ^K

L

[ ] ^P

T

L₋₁

=

^(4.9)

ここで生成される KL-1においても節点の接続情報を読み取ることができるため，同様に[P]

を作成し，[KL-2]，[KL-3]…の生成を再帰的に繰り返すことで，最終的には階層化された複数 の係数行列を得ることができる．この操作をコースニングと呼ぶ．なお，[K0]が最下層の係 数行列であり，一方最上層の係数行列は本来解くべき連立一次方程式の係数行列と一致す る．

なお，CおよびFは次式を満たすように選ばれることが望ましい．

0

, ∈ ∩ ≠

∈ F each j S

_i

is either in C or S

_j

C

_i

i

each

For

(4.10)

式(4.10)に従うCおよびFを選出するアルゴリズムは，図4.4の通りとなる．

( )

endfor

do U S j all for endfor

endfor

do U S k all for

j U U j F F Set

do S j all for

i U U i C C Set

with U i i Find

S Set

U F

C Set

U until repeat

j j

i k k

j T i

i T

i i

1

: 1

: ,

: ,

max ,

, 0 ,

0 0

−

=

∩

∈ +

=

∩

∈

−

= +

=

∈

−

= +

=

∈

=

Ω

=

=

=

=

λ λ

λ λ

λ λ

図4.4  CおよびFの選定アルゴリズム

ここで，Uは一時的に節点を保持しておく作業用スタックであり，スタック内の全節点がC もしくは F に振り分けられて，スタックが空になるまで上記のアルゴリズムは繰り返され る．まず，全節点について，節点iに強く依存する節点の数λiをセットする．なお，節点m から見て強い接続にある節点nは，必ずしも節点nから見ても節点mが強い接続であると は限らないため，節点mから見て節点nが強い接続にある場合，「節点nは節点mに強く依 存する」と呼び，また節点nから見て節点mが強い接続にある場合，「節点nは節点mに強 く影響する」と呼ぶ．次に，λiが最大となる節点iを見つけ，それをCに加え，スタックU から節点iを削除する．同様に節点iに強く依存する節点jはFに加えられ，スタックUか ら節点jは削除される．さらに，まだスタックUに残っていて節点jに強く影響する節点に ついては，その節点のλの値を+1 だけ更新する．ここまでの処理で i に影響する節点がま だスタックUに残っている場合には，その節点のλの値を-1だけ更新する．更新されたλ の値に従ってλが最大となる節点i を見つけ，また同様にCおよびFへの振り分けを繰り 返す．ただし，このアルゴリズムだけでは完全に式(4.10)を満たすことはできないため，随 時Fの節点をCの節点に変換する処理を行う[4.10]．Fの節点iが別のFの節点jに強く依 存する場合，jに強く影響する節点全てがCiに含まれるかどうかをチェックして，含まれな い場合には，節点jを一旦Cに変換し，Ciに加える．ここで接続関係が更新されるため，節 点iに強く依存する全てのFの節点を再チェックし，それが全てCiの節点に強く依存して

いるならば，正式に節点 jをCとして採用し，次の節点 Fに移る，一方，Ci以外の節点に 強く接続する節点が一つでもある場合には，節点jはFに戻される．

4.2.2   静磁界解析における代数マルチグリッド法

ここでは，領域分割に基づいたAMG法を静磁界解析に適用することを考える．静磁界解 析における基礎方程式は，マクスウェル方程式から磁気ベクトルポテンシャルAを用いて，

A法で次式のように表される．

M J

A r r r

×

∇ +

=

×

∇

×

∇

₀

0

1

µ

^(4.11)

ここで，μ0は真空の透磁率，J0は強制電流密度，Mは磁化である．式(4.11)を辺要素有限要 素法で離散化して得られた連立一次方程式の反復解法として代数マルチグリッド法を適用 する．

[ ] ^K

e

{ } { } ^A = ^b

^(4.12)

ここで，[K_e]は式(4.11)を辺要素有限要素法で離散化して得られる連立一次方程式の係数行 列，{A}は解ベクトル，{b}は右辺ベクトルである．

節点要素有限要素法では係数行列の構成要素から節点同士の接続関係を読み解くことが できるため，係数行列を元にAggregateを作成することができる．しかし，電磁界解析で用 いられる辺要素有限要素法では上位階層と下位階層の辺同士を関係づけるProlongation行列 が必要であり，かつ係数行列[Ke]には節点の接続情報が含まれないため，係数行列から Aggregateを作成することができず，節点有限要素法における AMGの手法をそのまま電磁 界解析に適用することはできない．そこで，メッシュモデルの情報を反映した擬似的な係 数行列[Kn

TOP]を用意する．[Kn

TOP] は節点の接続情報と要素辺の長さから次式に基づいて擬 似的に作成される係数行列（擬似行列）である[4.11]．

( ) _⎪ ^⎪ ( )

⎩

⎪ ⎪

⎨

⎧

=

−

≠

−

=

∑

≠

j i K

j edge i

K

j i

ij top n ij ij

top n

1

1

(4.13)

ここで，||edgeij||は節点 i，j の距離を表す．[Kn

TOP]の作成方法は任意であるため，任意の支 配方程式を節点有限要素法で離散化して得られる係数行列を用いることもできる．ただし，

順次下位階層の擬似行列を作成する過程で偶然にゼロの要素が発生すると，その行と列に 対応する節点の接続関係が失われてしまう．そのため，[Kn

TOP]には優対角正定値行列を用い ることが望ましい．この擬似行列に対しては節点有限要素法におけるAMG法を適用するこ とができるため，上位階層から順次下位階層の擬似行列を作成することができる．下位下

層の擬似行列[Kn

L-1]は，上位階層の擬似行列[Kn

L]と節点要素有限要素法におけるProlongation 行列[Pn]を用いて次式のように生成される．

[ ] [

⁻¹

⁼

^{→ −1}

] [ ][

n^L→L−1

]

L n L T L n L

n

P K P

K

^(4.14)

擬似行列は擬似的にメッシュとみなすことができるため，[Kn L]と[Kn

L-1]から上位階層と下 位階層の辺同士の関係を読み解き，辺要素有限要素法における Prolongation 行列[Pe]を次の ように作成する．

( ) ( ( ) ( ) )

( ) ( )

( )

⎪⎩

⎪⎨

⎧

=

−

= +

=

otherwise i ind i ind j

i ind i ind j P_{e ij}

0

, 1

, 1

1 2

2

1 (4.15)

ここで，上位階層の辺i は節点 i1，i2を両端に持ち，辺の向きがi1→i2である場合，i=(i1,i2) と表す．下位階層の辺jについても同様に表すものとする．辺同士の対応関係では向きが重 要となるため，Prolongation行列 は-1，0，+1で構成される．この辺要素有限要素法におけ る[Pe]の作成方法はReitzingerによって提案された手法[4.11]であり，次式を満足することが 証明されている．

H n h H

H

e

grad u grad P u

P =

^(4.16)

ここで，uHは下位階層における節点値，gradH，gradhはそれぞれ下位階層，上位階層にお ける勾配である．また，この関係式を図式化したものを図4.5 に示す．QH，Qhは下位階層 または上位階層の擬似的なメッシュに対して節点有限要素法のスカラー形状関数で離散化 した場合に定義される関数空間，VH，Vhは下位階層または上位階層の擬似的なメッシュに 対して辺要素有限要素法のベクトル形状関数で定義される関数空間である．下位階層にお いて節点値の勾配で表現されるベクトル量は，[Pn]を経由するパスと[Pe]を経由するパスの2 通りで上位階層の関数空間に補間できることがわかる．つまり，ベクトル量の非回転成分 は，[Pe]により上位階層に補間されることが保証され，その補間精度は非回転成分を節点値 に変換して[Pn]で補間したものと等価であることを示している．電磁界解析において，ベク トル量の非回転成分は係数行列の基底ベクトルにより定義される関数空間のカーネルであ り，そのカーネルの誤差は反復解法の収束性に大きく影響することが知られている[4.12]．

Reitzingerの提案手法は，そのカーネルの誤差が補間により増幅しないことを保証した手法 である．

Q

_H

V

_H

Q

_h

V

_h

grad

_H

grad

_h

P

_n

P

_e

Q

_H

V

_H

Q

_h

V

_h

grad

_H

grad

_h

P

_n

P

_e

図4.5  [Pn]と[Pe]の関係

12 11 10 9 8 7 6 5 4 3 2 1

0 0 0

1 0 0

1 0 0

0 0 1

0 0 1

0 0 0

1 0 0

0 0 0

0 1 0

0 0 0

0 0 1

0 0 0

C B A

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢

⎣

⎡

+ + +

−

+

−

−

Fine Level Coarse Level

Agg1

Agg3

A Agg2 C

B Agg1

Agg3

Agg2 1

2

3

4

5 7 6

8

9 10 11 12

[P

_e

]:=

12 11 10 9 8 7 6 5 4 3 2 1

0 0 0

1 0 0

1 0 0

0 0 1

0 0 1

0 0 0

1 0 0

0 0 0

0 1 0

0 0 0

0 0 1

0 0 0

C B A

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢

⎣

⎡

+ + +

−

+

−

−

Fine Level Coarse Level

Agg1

Agg3

A Agg2 C

B Agg1

Agg3

A Agg2 C

B Agg1

Agg3

Agg2 1

2

3

4

5 7 6

8

9 10 11 12 Agg1

Agg3

Agg2 1

2

3

4

5 7 6

8

9 10 11 12

[P

_e

]:=

図4.6  辺要素有限要素法におけるProlongation行列

図4.6に辺要素有限要素法におけるProlongation行列の作成例を示す．なお，擬似行列に おける要素辺の向きは，節点番号の小さい方から大きい方を正とする．

あ と は 下 位階 層 の 擬 似行 列[Kn

L]を 作 成 す る 場合 と 同 じ よう に 辺 要 素有 限 要 素 法の Prolongation行列[Pe]を元に階層化された辺要素有限要素法の係数行列[Ke

L]を作成していく．

[ ] [

⁻¹

⁼

^{→ −1}

] [ ][

e^L→L−1

]

L e L T L e L

e

P K P

K

^(4.17)

なお，最上層の[Ke

TOP]には最初に与えられた係数行列を使用する．

4.2.3   過渡応答磁界解析における代数マルチグリッド法

ここでは，領域分割に基づいたAMG法を過渡応答磁界解析に適用することを考える．過 渡応答磁界解析における基礎方程式は，その表現方法に磁気ベクトルポテンシャル A のみ を用いる A 法と，磁気ベクトルポテンシャル A と電気スカラポテンシャルφを用いる A- φ法の2通りがある．一般に過渡応答磁場解析においてはA-φ法の方がA法より計算収束 性がよいことが多く[4.13]，本節においてもA-φ法によるAMG法の定式化を行う．A-φ法 による基礎方程式は次式で表される．

M t J

A A ⎟= +∇×

⎠

⎜ ⎞

⎝

⎛ +∇ + ∂

×

∇

×

∇ ₀

0

1

σ φ

µ

^(4.18)

= 0

⎭ ⎬

⎫

⎩ ⎨

⎧ ⎟

⎠

⎜ ⎞

⎝

⎛ + ∇

− ∂

⋅

∇ σ φ

t

A

(4.19)

ここで，μ0は真空の透磁率，J0は強制電流密度，Mは磁化，σは電気伝導率である．A- φ法では式(4.18)と式(4.19)の連立方程式を考え，辺要素有限要素法と節点要素有限要素法で 離散化して得られる次式の連立一次方程式にAMG法を適用する．

⎭⎬

⎫

⎩⎨

=⎧

⎭⎬

⎫

⎩⎨

⎥⎧

⎦

⎢ ⎤

⎣

⎡

0 f A K K

K K

nn ne

en ee

φ

^(4.20)

ここで，[Kee]，[Ken]は式(4.18)の辺要素有限要素法における連立一次方程式の係数行列，

[Kne]，[Knn]は式(4.19)の節点要素有限要素法における連立一次方程式の係数行列，{A}はベク トルポテンシャルAの解ベクトル，{φ}は電気スカラポテンシャルφの解ベクトル，{f}は 荷重ベクトルである．A-φ法に AMG 法を適用する場合，要素辺に未知数を置く A の Prolongation 行列[Pe]と節点に未知数を置くφの Prolongation 行列を[Pn]同時に考慮する Prolongation行列が必要となる．電流ベクトルポテンシャルTと磁気スカラポテンシャルω を用いるT-ω法による電磁界有限要素法解析において，辺要素と節点要素を用いたAMG法 は既に提案されており[4.14]，本章ではそれと同等の手法をA-φ法に適用する．本章で用い るA-φ法におけるProlongation行列[Pe+n]を次式に示す．

[ ]

_⎥

⎦

⎢ ⎤

⎣

=⎡

+

n e n

e P

P P

0

0 _(4.21)

[Pe]は前節で述べた辺要素有限要素法特有のProlongation行列であり，[Pn]は[Pe]を作成する 過程で用いた擬似行列を元とする節点有限要素法における Prolongation 行列である．なお，

Restriction行列[Re+n]は [Re+n]=[Pe+n]^Tとする．静磁界解析におけるA法と過渡応答磁界解析 におけるA-φ法では最上層の係数行列やProlongation行列，Restriction行列は異なるものの，

AMG法の処理手順自体は全く同じとなる．

4.2.4   代数マルチグリッド法を前処理とした共役勾配法

前処理付き共役勾配法（PCG 法）[4.15]は，前処理行列[M]を用いて解くべき連立一次方 程式を同値で性質の良い方程式に変換してから解く方法であり，[M]^-1が変換前の係数行列 をよく近似している場合に収束性が向上する．辺要素有限要素法で離散化して得られた連 立一次方程式に ICCG 法を適用する場合，係数行列[Ke]に不完全コレスキー分解を用いて [M]^-1を求め，図4.7に示すようにPCG法の処理の中で{r^~}=[M]^-1{r}の計算を行う．一方，

AMG法をPCG法の前処理として用いる場合，[M]にはそのまま[Ke]を用いて図4.8に示すよ うにPCG法の処理の中で[M]{r^~}={r}の連立一次方程式を解くが，この反復解法にAMG法 のVサイクル処理を適用する．以降，PCG法の前処理にAMG法を用いる方法をAMGCG 法と呼ぶ．辺要素有限要素法による電磁界解析では解くべき連立一次方程式の係数行列が 特異となるため，VサイクルAMG法の全ての階層において一意に解が決定されるとは限ら ない．そのため，実際には次式のように[M]に[Ke]の対角部分をわずかに補正して代数的に 正則化させたものを用いる．

[ ] [ ] [ ] ^M = ^K

e

+ α ^D

^(4.22)

ここで，[D]は[K_e]の対角成分のみを抽出した対角行列である．αには 1.0×10^-6〜1.0×10^-8 程度の小さい値が用いられ，[M]と[Ke]が極度に変わらないようにする．なお，AMG法の最 上層以下の係数行列は，この正則化された最上層の係数行列を基に作成される．

この正則化により，静磁界解析においてもベクトル量の非回転成分の誤差は無視できな くなる．そのため，Reitzingerの提案手法が有効となる．

電磁界解析におけるマルチグリッド法のスムーザーには，ハイブリッドスムーザー[4.16]

や AFW スムーザー[4.17]など電磁界解析特有の収束性悪化を改善するものが提案されてい る．しかし，ハイブリッドスムーザーを用いるには[Ke]に対してクーロンゲージによる正則 化処理が必要となる．また，ハイブリッドスムーザー，AFWスムーザー共に効率よく計算 を行うには[Ke]の他に補助的な行列を用意する必要があり，係数行列と同等サイズのメモリ を重複して確保しておかなければいけない．一方，スムーザーに標準的な定常反復法であ るヤコビ法やガウス・ザイデル法を用いる場合は，最もメモリを消費する最上層の係数行 列については上述の代数的な正則化処理に従って[Ke]の対角部分を補正した値を保存する ベクトル列のみを確保しておくだけでよいためメモリの使用効率がよい．また，標準的な 定常反復法はハイブリッドスムーザーやAFWスムーザーなどの特殊なスムーザーに比べ収 束性では劣るものの，安定的に収束することが経験的にわかっている．よって，全体の計 算時間および処理速度に大きく影響するメモリ使用効率に着目して，本章では係数行列の 対称性も加味してスムーザーにSSOR法を用いている．

( )

( )

( )

( )

for end

p r

p

r r

r r

r M r

p K r

r

p x

x

p K p

r r

do b r

until n

for p

r M r

x K b r x

n n n n

n n

n n n

n n

n e n n n

n n n n

n e n

n n n

n e

β β

α α α

ε

+

=

=

=

−

=

+

=

=

≤

=

=

=

−

=

=

+

+ +

+

− +

+ +

−

~

~ ,

~ ,

~

,

~ ,

: ,

1 , 0 0

~ 0

1

1 1

1 1 1

1 1 0

0 1 0

0 0

0

L

IC:

IC:

( )

( )

( )

( )

for end

p r

p

r r

r r

r M r

p K r

r

p x

x

p K p

r r

do b r

until n

for p

r M r

x K b r x

n n n n

n n

n n n

n n

n e n n n

n n n n

n e n

n n n

n e

β β

α α α

ε

+

=

=

=

−

=

+

=

=

≤

=

=

=

−

=

=

+

+ +

+

− +

+ +

−

~

~ ,

~ ,

~

,

~ ,

: ,

1 , 0 0

~ 0

1

1 1

1 1 1

1 1 0

0 1 0

0 0

0

L

IC:

IC:

図4.7  ICCG法の処理手順

( )

( )

( )

( )

for end

p r

p

r r

r r

r r

M

p K r

r

p x

x

p K p

r r

do b r

until n

for p

r r M

x K b r x

n n n

n

n n

n n n

n n

n e n n

n

n n n

n

n e n

n n n

n e

β β

α α α

ε

+

=

=

=

−

=

+

=

=

≤

=

=

=

−

=

=

+

+ +

+ +

+ +

~

~ ,

~ ,

~

,

~ ,

: ,

1 , 0 0

~ 0

1

1 1

1 1

1 1 0

0 0

0 0

0

L

AMG:

AMG:

( )

( )

( )

( )

for end

p r

p

r r

r r

r r

M

p K r

r

p x

x

p K p

r r

do b r

until n

for p

r r M

x K b r x

n n n

n

n n

n n n

n n

n e n n

n

n n n

n

n e n

n n n

n e

β β

α α α

ε

+

=

=

=

−

=

+

=

=

≤

=

=

=

−

=

=

+

+ +

+ +

+ +

~

~ ,

~ ,

~

,

~ ,

: ,

1 , 0 0

~ 0

1

1 1

1 1

1 1 0

0 0

0 0

0

L

AMG:

AMG:

図4.8  AMGCG法の処理手順

4.3   数値解析による検証

4.3.1   扁平要素に対する代数マルチグリッド法を前処理とした

共役勾配法の収束性の評価

積層構造を有する電気機器の新たな高速解法として AMGCG法を評価するため，まず簡 易的な形状のメッシュモデルについて，収束までに要した反復回数およびCPU時間をICCG 法と比較する．メッシュモデルには，積層構造を考慮しないメッシュと積層構造を考慮し たメッシュの 2 種類を用いる．積層構造を考慮しない場合，コア材は単純なブロック材と なるため，極端に扁平した要素を含むことなく数十万程度の標準的なメッシュで計算する ことができる．一方，積層構造を考慮した場合，全長が数mm〜数十mmの解析モデルに対 して0.1mm 未満の鋼板間の空隙をメッシュで表現しなければいけないため，その空隙近傍 のメッシュにひきずられて全体の要素数が百万を超える莫大なものとなる．その上，積層 鋼板およびその空隙は数層のメッシュで分割する必要があるため，極端にアスペクト比の 大きい要素が多数生成される．アスペクト比の大きい要素は，反復解法の収束性を悪化さ せることが多いことが知られている．そのため，積層構造を考慮したメッシュモデルでの 計算は，反復解法の収束性の悪さとメッシュモデル規模の大きさの両方から，積層構造を 考慮しない場合に比べて極端に計算時間がかかる傾向にある．

数値解析による検証では，反復解法の収束判定については，ICCG法，AMGCG法共に相 対残差ノルムを CG 法の反復ごとにチェックし，その値が 1.0×10^-8以下となった時点で収 束とみなした．また，電磁鋼板の非線形性を考慮するため，ニュートン・ラプソン法[4.18]

を用いているが，その収束判定については，解ベクトルの2ノルムの変動率が 1.0×10^-3以 下になった時点で収束とみなした．AMGCG法については，スムーザーとしてSSOR法を用 い，その加速係数は1.2，反復回数は1回とした．Vサイクルの繰り返し回数は1回とし，

階層数は最下層の未知数が 100 未満になるように調整する．最下層では十分に未知数が少 ないため，解法にはLU分解を用いる．

図 4.9 に積層鉄心ベンチマークモデルのメッシュモデル[4.22]を示す．電磁鋼板は図 4.9 のz方向に積層されており，鋼板の厚みはギャップを含めて0.5mmである．

表4.1に，積層構造を考慮しない場合と考慮した場合のメッシュモデルにおけるICCG法 と AMGCG法の収束性の比較を示す．ここで，速度比はAMGCG法の CPU時間に対する ICCG 法の CPU時間の比率とした．なお，積層構造を考慮した場合の積層鋼板の占積率は 96%とした．表4.1から，AMGCG法の反復回数がICCG法に比べて圧倒的に少ないことが わかる．前処理の計算コストは，IC法に比べて AMG 法の方が数倍大きいため，速度比は 反復回数の比率に比べて小さくなる．しかし，積層構造を考慮しない場合と考慮した場合 のどちらにおいても AMGCG法はICCG法より高い計算速度を示している．積層鋼板を考

慮したメッシュモデルは，考慮しないメッシュモデルに比べて要素数は 1.76倍程度の増加 であるが，CPU時間はICCG法，AMGCG法共に10倍以上増加しており，扁平要素が混在 することにより収束性が著しく悪化することがわかる．ただし，ICCG法とAMGCG法の速 度比を比較すると，扁平要素を含む前に比べて，含んだ方が速度比の差が大きくなってい ることがわかる．

Lamination Factor 90%

Lamination Factor 50%

0.5mm a) Solid iron core mesh model

b) Laminated iron core mesh model

0.5mm Iron core

Coil

Lamination Factor 90%

Lamination Factor 50%

0.5mm a) Solid iron core mesh model

b) Laminated iron core mesh model

0.5mm Iron core

Coil

図4.9  積層鉄心ベンチマークモデルのメッシュ図（1/8部分モデル）

表4.1  積層鉄心ベンチマークモデルにおける収束特性 Mesh type Number of

elements Solver type Total Number of CG iterations

Total CPU

time(s) Speed ratio ICCG 2,113 437.0

Standard 829,960

AMGCG 522 393.3 1.11 ICCG 165,703 57,103.0 Lamination 1,464,285

AMGCG 15,292 20,558.5 2.78 Computer used : Intel Xeon 3.4GHz / 8GByte RAM

0 100 200 300 400 500 600

0 500 1000 1500 2000 2500 3000

Aspect ratio of an element in the gap

CPU Time, sec

ICCG AMGCG

図4.10  要素のアスペクト比に対するAMGCG法の計算時間

そこで，メッシュの扁平に対する ICCG法およびAMGCG法の収束性の変化を検証する ために，積層鋼板を考慮したメッシュで占積率を変更した時の ICCG法およびAMGCG法 の計算時間を比較する．ここでは，積層鋼板および積層鋼板間の空隙の分割層数は 3 で一 定とし，節点移動により占積率を変更する．つまり，全体要素数は占積率に関わらず常に 一定であり，占積率が小さい場合，扁平要素数が減少し，占積率が大きい場合，扁平要素 数が増加する．図4.9に，占積率が50%と90%の場合の積層鋼板のメッシュの拡大図を示す．

図4.10に積層鋼板間の空隙に位置する要素のアスペクト比とICCG法およびAMGCG法 の計算時間を示す．図 4.10は，各占積率における積層鋼板間の空隙に位置する要素のアス ペクト比を確認し，そのアスペクト比と CPU 時間の相関を取ったものとなる．ICCG 法は 要素の扁平に対して急激に計算時間が増加しているが，AMGCG法はICCG法に比べてその 増加が緩やかであることがわかる．この図から，AMGCG法はICCGに比べて扁平要素混在 時の収束性劣化に対する耐久性が高いことがわかる．

次に，収束性に大きく影響すると考えられるAggregateのパラメータθについて，AMGCG 法の処理時間との関係を検証する．図4.11にθの値を0.1〜0.7まで変化させた時の収束ま でに要した計算時間を示す．検証モデルには，同様に積層鉄心ベンチマークモデルを用い，

非線形反復 1 回目の収束に要した時間を計測した．θの値と計算時間には概ねトレードオ フの関係があり，θ=0.56が計算時間最小となる．θの値が小さい場合，Aggregateのサイズ は小さくなりPeは密な行列となる．そのため，Peによる補間精度が向上し，全体の反復回 数は減るものの，1反復あたりの計算時間が長くなる．一方，θの値が大きい場合，Aggregate のサイズは小さくなりPeは粗い行列となる．そのため，Peによる補間精度が劣化し，反復 回数が増加するものの，1反復あたりの計算時間が短くなる．今回の検証では，θ=0.66 を 境に急激に収束性が悪化し，反復回数が数倍以上となった．

扁平要素をほとんど含まない通常のメッシュモデルではθの値に0.25 程度を用いること が多いが，扁平要素を多く含む特殊なメッシュモデルにおいては，その最適値が異なるこ とがわかる．θの値の選び方によっては計算時間が 30%以上異なることから，θが収束性 を大きく左右する制御パラメータであることが明らかになった．

2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000

0.10 0.20 0.30 0.40 0.50 0.60 0.70

Theta

CPU time, sec

図4.11  制御パラメータ値に対するAMGCG法の計算時間

4.3.2   積層構造を有する電気機器モデルへの適用

本節では，実用的な形状のメッシュモデルに対する AMGCG法の性能を静磁界解析およ び過渡応答磁界解析で評価する．

図4.12にリアクトルモデルと積層鋼板部分のメッシュ拡大図を示す．電磁鋼板は図4.12 の z方向に積層されており，鋼板の厚みはギャップを含めて 0.5mm，占積率は96%として いる．電磁鋼板には 50H1300 相当の材料特性を用い，スペーサーおよびコイルの材料特性 は空気として扱う．ここではコイルに 750AT の一定電流を流した場合の静磁界解析を行っ た．

表4.2に，ICCG法およびAMGCG法の収束性の比較を示す．AMGCG法の合計反復回数 はICCG法の1/16であり，ICCG法より収束性が良いことがわかる．AMGCG法はICCG法 に比べて1反復あたりの演算量は多いが，合計CPU時間は1/4以下となっており良好な計 算速度を示している．

図4.13および図4.14にコアにおける磁束密度の積層方向成分の分布を示す．積層構造を 考慮する効果を確認するため，積層構造を考慮しないメッシュモデルでの結果を併載する．

積層構造を考慮しないメッシュモデルの結果では，磁束がコア内部まで浸透しているのに 対し，積層構造を考慮したメッシュモデルでは積層構造の効果で，ほとんど磁束が内部に 浸透していないことがわかる．

表4.2  リアクトルモデルにおける収束特性 Number of

elements Solver type Total Number of CG iterations

Total CPU

time(s) Speed ratio ICCG 46,535 19,582.9

1,755,810

AMGCG 2,804 4,542.1 4.31 Computer used : Intel Xeon 3.4GHz / 8GByte RAM

Insulation Gap Coil

Iron core Spacer

Z X Y

Insulation Gap Coil

Iron core Spacer

Z X Y

図4.12  リアクトルモデル

a) Solid iron core

b) Laminated iron core

(T)

a) Solid iron core

b) Laminated iron core

(T)

図4.13  鉄心コア表面における磁束密度の積層方向成分の分布

0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18

0 5 10 15

Distance, mm

Magnetic Flux Density Bz, T

NoLamination

Lamination

図4.14  鉄心コア端部の磁束密度の積層方向成分

図4.15 に電気学会の回転機のバーチャルエンジニアリングのための電磁界解析技術調査 専門委員会で検証された埋込み磁石構造電動機モデル[4.19]を元に作成したIPMモータモデ ルとそのローターおよびステーターのメッシュ拡大図を示す．このモデルではローターの 積厚がステーターより5mmだけオーバーハングしているため，ローター端部において積層 方向への漏れ磁束の影響が現れる．積層鋼板の厚さは 0.5mm，占積率は 96%としている．

鋼板および磁石の材料特性は電気学会で検証された埋込み磁石構造電動機モデルと同じと し，ローター初期位置における無負荷静磁界解析を行った．

表4.3にICCG法とAMGCG法の収束性の比較を示す．AMGCG法の合計反復回数はICCG 法の1/13，合計CPU時間はICCG法の1/3であり，良好な収束性，計算速度を示している．

図4.16，図4.17にステーター表面の磁束密度の積層方向成分の分布と，ギャップに面す る面の中央部における磁束密度の積層方向成分の比較を示す．ローター端部における積層 構造の重要性を確認するため，積層を考慮したモデルと合わせて積層を考慮しないソリッ ドローターモデルの結果と比較する．ソリッドローターモデルではオーバーハング部で積 層方向に流れ込む磁束が大きいのに対し，積層モデルでは積層構造の効果により積層方向 の磁束が制約されていることがわかる．

X Y

Z X

Y

Z

Stator Rotor Magnet

X Y

Z X

Y

Z X

Y

Z X

Y

Z

Stator Rotor Magnet

図4.15  埋込み型磁石構造電動機モデル

表4.3  埋込み型磁石構造電動機モデルにおける収束特性 Number of

elements Solver type Total Number of CG iterations

Total CPU

time(s) Speed ratio ICCG 94,780 90,264.0

4,655,600

AMGCG 6,982 31,078.3 2.90 Computer used : Intel Xeon 3.4GHz / 8GByte RAM

b) Laminated rotor a) Solid rotor

(T)

b) Laminated rotor a) Solid rotor

(T)

図4.16  ローター表面中央部における磁束密度の積層方向成分の分布