第 3 章 解析モデルの形状変更に適した メッシュモーフィング手法

メッシュモーフィング手法

3.1   緒言

本章では，形状変更に適したメッシュモーフィング手法について述べる．

電気機器設計では，数値解析と最適化手法とを組み合わせて形状最適化設計が行われる ことが多く，電磁界数値解析では可動部の移動を考慮した過渡応答解析を行うことも多い ため，メッシュモデルの変形作業の効率化に対する要求は多い．しかし，電磁界数値解析 では機器周辺の空間部分もメッシュモデル化する必要であるため，一般に解析対象のメッ シュとその周辺の空間メッシュの整合性を保ちつつ要素の扁平を極力抑えてメッシュモデ ルを変形することは困難である．その対応策として，メッシュモデルを直接変形するメッ シュモーフィング手法がある[3.1]-[3.11]．メッシュモーフィング手法とは，メッシュモデル をベースに形状変更を行うことであり，一度ベースとなるメッシュモデルを作成すれば，

そのベースメッシュを元に形状の異なるメッシュモデルを自由に作成することができる．

メッシュモーフィング手法では，物性境界部分の節点に変形量を与え，ある支配方程式に 従ってその物性境界の内外の節点の移動量を算出し，節点分布が全体的に滑らかになるよ うに節点配置を変更する．その支配方程式の種類にはラプラス方程式（ラプラス法）

[3.1]-[3.4]やバネの釣り合い式（バネ法）[3.5]-[3.8]，構造弾性方程式（弾性法）[3.9]-[3.11]

などがある．標準的なメッシュモーフィング手法では，形状変更に伴って節点の座標値が 変わるのみで節点の接続関係や節点数，要素数は変わらないため，有限要素法[3.12]で必要 な境界条件や各種条件の再設定も不要である．この利便性のため，流体解析や構造解析の 分野でのメッシュモーフィング手法の適用事例は多く，その手法も様々なものが提案され ているが，電磁界解析分野での適用事例はまだ少ない．

そこで本章では，空間を含む電磁界数値解析特有のメッシュモデルで既存手法を検証し，

数値実験から既存手法では要素の扁平化に対する高い耐久性と少ない処理時間を同時に満 たすことが困難であることを明らかにした．さらに，その解決策としてラプラス法とバネ 法を適切に組み合わせたハイブリッド手法を開発した．ハイブリッド手法は，計算量が少 ないラプラス法をベースとした手法であり，メッシュ変形の過程で局所的に要素の著しい 扁平が発生したエリアを，耐久性の高いバネ法で修復する方法である．提案手法は処理時 間の短縮と耐久性を同時に期待できるが，バネ法で修復するエリアの範囲を制御するパラ メータがその性能を大きく左右する．数値実験では，変形量の大きいメッシュモデルを例

題に取り上げ，提案手法の制御パラメータの特性について明らかにした上で，提案手法の 有効性を検証した．

3.2   メッシュモーフィング手法

メッシュモーフィング手法には様々な手法が提案されており，物性境界の変形に合わせ てその内外の節点移動を決定する支配方程式としてラプラス方程式やバネの釣り合い式，

構造弾性方程式がある．以降では，これらの支配方程式に基づくメッシュモーフィング手 法をそれぞれラプラス法，バネ法，弾性法と呼ぶこととする．本節ではこの 3 つの従来手 法について概説する．

3.2.1 ラプラス方程式を基にしたメッシュモーフィング手法

ラプラス法は，全節点の移動量ベクトルを次式のラプラス方程式に基づいて決定する方 法である．

2

r = 0 r

∇ q

^(3.1)

ここで，

q r

は次式で定義される各節点の移動量ベクトルである．

ξ

0

ξ

^{r r} r = _m−

q ^(3.2)

ここで，

ξ

^r_mは変形後の各節点の位置ベクトル，

ξ

^r₀は変形前の各節点の位置ベクトルである．

全節点の

q r

は，変形境界面上の節点の移動量ベクトルを境界条件として，式(3.1)を節点要素 有限要素法[3.13]で離散化して得られる次式の連立一次方程式を反復解法で解くことにより 得られる．

[ ]

M

{ } { }

q = 0 (3.3)

ここで，[M]は正定値の対称マトリクス，{q}は解ベクトル，{0}は右辺ベクトルである．な お，移動量ゼロの物性境界に囲まれた領域など全く変形の影響が及ばない部分の節点の移 動量ベクトルについては，予め式(3.3)の解ベクトルから除外される．本章では，連立一次方 程式の反復解法に ICCG 法[3.14]を用いる．ラプラス法による処理手順を図 3.1に示す．メ ッシュの変形は，数回の移動量に分割して行う．図3.1の処理は，移動節点が最終位置に到 達するまで繰り返される．ある移動量でメッシュを変形させて一部の要素が反転もしくは 著しく扁平し，計算不能なメッシュとなってしまった場合には，移動した全節点の位置を 一つ前の位置まで戻し，変形境界面上の節点の移動量を半分にして全節点の移動量ベクト ルを再計算する．移動量を半分にした場合，計算回数（図3.1における反復回数）は倍とな る．計算回数を少なくし，全処理時間を短縮するためにも，メッシュの変形による要素の 反転，扁平が発生しにくい耐久性が要求される．

ラプラス法は後述する手法と比べて演算量が少なく処理時間は小さいものの，移動節点 の移動量が大きい場合に要素が反転してしまうなどメッシュの変形に対する耐久性に問題 がある．ラプラス法における耐久性向上の工夫は幾つか提案されている[3.1][3.4]ものの，あ らゆるメッシュ，変形量に対応できるものではないと考えられる．

Determine deformation using Laplace method

Update positions of grids Collapsed mesh

exists

YES

NO

Move boundary to final position

YES

End Update matrix using

updated meshes

NO

Reset displacement vectors as boundary condition Start

Determine deformation using Laplace method

Update positions of grids Collapsed mesh

exists

YES

NO

Move boundary to final position

YES

End Update matrix using

updated meshes

NO

Reset displacement vectors as boundary condition Start

図3.1 ラプラス法によるメッシュモーフィング処理手順

3.2.2 バネの釣り合い式を基にしたメッシュモーフィング手法

バネ法は要素辺を線形ばねとみなし，変形境界面以外の節点の位置をばね力の釣り合い によって決定する方法である．線形ばねのばね力は次式で定義される．

{ } [ ]

linear^ij

{ }

^ij

ij

linear

k q

F =

(3.4)

{ } {

z^j

}

^T

j y j x i z i y i x ij

linear F F F F F F

F = , , , , , ^(3.5)

{ } {

z^j

}

^T

j y j x i z i y i x

ij q q q q q q

q = , , , , , ^(3.6)

ここで，{F^ijlinear}は要素辺eijの節点i，jに作用するばね力，{q^ij}は節点i，jの移動量ベクト ルである．[k^ijlinear]は要素辺 eijに関する線形ばねの弾性係数マトリクスであり，図 3.2 に示 す要素辺eijの場合，次式のように定義される．

[ ] ^{[ ] [ ]}

[ ]

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎣

⎡

−

−

−

=

=

φ φ θ

φ θ

φ φ θ

φ θ

cos sin sin

sin cos

cos sin sin

sin cos

1

2

R

R l R

k ^T

ij ij

linear

(3.7)

[k^ijlinear]の値は，e_ij長さ lijの二乗の逆数に比例するため，要素が扁平して節点が接近するほ

ど線形ばねが硬くなる仕組みとなる．

X

Y Z

i

j

θ Φ

X

Y Z

i

j

θ

Φ

この仕組みにより，要素の変形による節点同士の衝突は回避できるものの，要素辺や要 素面同士の交差を回避することはできない．そこで，三角形要素の角にねじりばねを追加 する．図3.3に，2次元三角形要素において線形ばねとねじりばねを組み合わせた例を示す．

ねじりばねは各節点に追加され，そのばね係数は次式のように要素辺のなす角の関数で定 義される．

jik i jik

C

i

θ sin

2

= 1

(3.8)

ここで，C^jikiは三角形tijkの節点iにおけるねじりばねのばね係数，θ^ijkiは節点iの角 jikの なす角である．θ^ijkiが 0°もしくは180°に近づくに従い，C^jikiの値は無限大へと近づき，

ねじりばねが硬くなる結果，要素辺や要素面同士の交差を回避する仕組みとなる．図3.3の 三角形tijkのねじりばね弾性係数マトリクス[k^ijktorsion]は次式の通りとなる．

[ ] [ ] [ ][ ] [ ]

[ ] ⎥ ⎥ ⎥

⎦

⎤

⎢ ⎢

⎢

⎣

⎡

=

⎥ ⎥

⎥

⎦

⎤

⎢ ⎢

⎢

⎣

⎡

−

−

−

−

−

−

−

−

−

−

−

−

=

=

ijk ki ijk j ijk i ijk

kj ki ki kj kj

kj ki

ki

jk jk

ji jk jk ji ji

ji

ik ik

ij ij

ij ik ij ik ijk

ijk T ijk

ijk ijk

torsion

C C C C

a a b b a

b a

b

a b

a a b b a

b

a b

a b

a a b b M

M C M k

0 0

0 0

0 0

(3.9)

i j

k

l

ij

1

l

jk

1 l

ki

1

ijk

θ

i

sin

2

1

ijk

θ

j

sin

2

1

ijk

θ

k

sin

2

1

i j

k

l

ij

1

l

jk

1 l

ki

1

ijk

θ

i

sin

2

1

ijk

θ

j

sin

2

1

ijk

θ

k

sin

2

1

図3.3 三角形要素における線形ばねとねじりばねの弾性係数

ここで，aij，bijは次のように定義し，xiは節点iのx座標，yiは節点iのy座標，xjは節点j のx座標，yjは節点jのy座標である．

ij i j ij

ij i j ij

l y b y

l x a x

= −

= −

(3.10)

線形ばねとねじりばねの組み合わせを 3 次元要素に拡張する場合，ねじりばねに相当す る壁面を要素内に挿入することを考える．図 3.4に，テトラ要素 Tpqrs内に三角形の壁面tsjq

を挿入した例を示す．tsjqは Tpqrsの要素面 Spqrに垂直となるように置かれている．先に説明 した2次元三角形要素のねじりばねの仕組みでこのtsjqの扁平を抑制することで，Tpqrs自体 の扁平を抑制することが可能となる．三角形の壁面は，テトラ要素の場合，図3.5に示すよ うに12枚挿入することが可能である．図 3.5において，色づけされた三角形が壁面，太線 で囲まれた三角形がその壁面と垂直となるテトラ要素面を表す．

q p

r s

j

q p

r s

j

図3.4 テトラ要素におけるねじりばね相当の三角形壁面

図3.5 テトラ要素における三角形壁面の種類

以上より，epqに関わる線形ばねおよびねじりばねのばね力の合計は，次式で表される．

{ } [ ] { } [ ] [

torsion^sjq

] { }

^pqrs

e pqrs pq pq

pq linear pq

total

k q B S k q

F

pq pqrs

∑

⊃

+

=

τ

]

[

_(3.11)

ここで，[B^prrspq]はテトラ要素Tpqrsの各節点の移動量ベクトルのうち要素辺epqに関わる成分 のみを抽出する 0，１で構成されるマトリクス，[k^sjqtorsion]は三角形の壁面 tsjqのねじりばね 係数マトリクス，{q^pqrs}は節点 p，q，r，s の移動量ベクトルである．[S]は移動量ベクトル {q^sjq}を{q^pqrs}に変換するマトリクスであり，次式のように表される．

[ ]

[ ]

( )

⎥ ⎥

⎥ ⎥

⎦

⎤

⎢ ⎢

⎢ ⎢

⎣

⎡ −

=

⎪ ⎪

⎭

⎪ ⎪

⎬

⎫

⎪ ⎪

⎩

⎪ ⎪

⎨

⎧

⎪ =

⎭

⎪ ⎬

⎫

⎪ ⎩

⎪ ⎨

⎧

0 0

0 0

0 0

0 1

0

I

I I I S

q q q q S q

q q

s r q p

q j s

λ

λ

^(3.12)

ここで，λは図3.4におけるpjの長さlpjとprの長さlprの比lpj/lprであり，Iは3×3の単位 行列である．[S]により，辺epr上に位置する節点jの移動量ベクトルqjは，節点p，rの移動 量ベクトル qp，qrから線形補間される．最終的には，全節点におけるばね力の釣り合い式 を考え，次式を節点要素有限要素法で離散化して得られる連立一次方程式をICCG法で解く．

{ } { }

∑

∈

=

epq

p pq

total

for all nodes q

F 0

_(3.13)

なお，バネ法による処理手順はラプラス法と同じであり，節点の移動量を決定する支配方 程式が違うのみである．

バネ法は大変形に対して耐久性が高い傾向にあるが，演算量がラプラス法より多く，処 理時間の面で課題が残る．

3.2.3 構造弾性方程式を基にしたメッシュモーフィング手法

弾性法はメッシュモデル全体を線形弾性体とみなし，変形境界面以外の節点の位置を構 造弾性方程式に基づいて決定する方法であり，その支配方程式は次式の通りとなる．

1. 変位とひずみの関係式

x w z u z

w

z v y w y

v

y u x v x

u

zx z

yz y

xy x

∂ + ∂

∂

= ∂

∂

= ∂

∂ + ∂

∂

= ∂

∂

= ∂

∂ + ∂

∂

= ∂

∂

= ∂

γ ε

γ ε

γ ε

, , ,

(3.14)

2. 応力のつりあいの式

0 0 0

∂ = + ∂

∂ + ∂

∂

∂

∂ = + ∂

∂ + ∂

∂

∂

∂ = + ∂

∂ + ∂

∂

∂

z y x

z y

x

z y

x

zy z zx

yz y

yx xy xz x

τ σ τ

τ σ τ

τ τ σ

(3.15)

3. 応力とひずみの関係式

( )

( )( )

( )

( )

( )

⎪⎪

⎪⎪

⎭

⎪⎪⎪

⎪

⎬

⎫

⎪⎪

⎪⎪

⎩

⎪⎪⎪

⎪

⎨

⎧

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎣

⎡

−

− −

− −

− −

−

−

−

−

−

− +

= −

⎪⎪

⎪⎪

⎭

⎪⎪⎪

⎪

⎬

⎫

⎪⎪

⎪⎪

⎩

⎪⎪⎪

⎪

⎨

⎧

zx yz xy z y x

zx yz xy z y x

E

γ γ γ ε ε ε

ν ν ν

ν ν ν ν

ν ν

ν ν

ν ν

ν ν

ν ν ν

ν ν

ν

τ τ τ σ σ σ

1 2

2 0 1

0 0

0 0

1 0 2

2 0 1

0 0

0

0 1 0

2 2 0 1

0 0

0 0

0 1 1

1

0 0

1 0 1 1

0 0

1 0 1 1

2 1 1

1

(3.16)

ここで，εは垂直ひずみ，γはせん断ひずみ，u，v，wはx，y，z方向の変位，σは垂直応 力，τはせん断応力，Eはヤング率，νはポアソン比である．これらの式から仮想仕事の原 理に基づいて導かれる次式の積分方程式を節点要素有限要素法で離散化し，変位を未知数 とした連立一次方程式をICCG法で解く．

( )

∫∫∫ ^ε

^*x

^σ

^x

^{+ ε}

^*y

^σ

^y

^{+ ε}

^*z

^σ

^z

^{+ γ}

^*xy

^τ

^xy

^{+ γ}

^*yz

^τ

^yz

^{+ γ}

^zx*

^τ

^zx

^{dxdydz = 0}

^(3.17)

ここで，ε^*は仮想垂直ひずみ，γ^*は仮想せん断ひずみである．なお，境界条件の扱いにつ いては，ラプラス法，ばね法と同じとなる．

弾性法では，ヤング率，ポアソン比の材料特性を各要素別にコントロールすることで要 素の極端な扁平を抑える．例えば，要素体積の逆数をヤング率として与えた場合，要素が 扁平して体積が小さくなるとその要素が硬くなり，それ以上の扁平を抑えることができる．

弾性法はバネ法に比べて演算量が少なく，ラプラス法よりメッシュの扁平に対する耐久性 も高く安定した手法ではあるが，制御パラメータとなる各要素の材料特性の決定方法に十 分なノウハウが求められる．例えば，単にヤング率を体積の逆数とした場合，体積が変わ らない要素のゆがみに対しては要素の硬度が変わらず変形を抑えることができないため，

更なる工夫が必要となる．なお，弾性法の処理手順もラプラス法，バネ法と同様となる．

3.2.4 ハイブリッド型メッシュモーフィング手法

前節で概説したラプラス法，バネ法，弾性法は，メッシュの変形に対して要素の極端な 扁平や反転を回避する耐久性，および全体の処理時間の短縮の点においてそれぞれ一長一 短があり，耐久性と処理時間の短縮を一つの手法で満たすことは難しいと考えられる．そ こで，これらの手法を適切に組み合わせて，耐久性と処理時間の短縮を同時に満たすこと とを考える．ここでは，アルゴリズムの明瞭性からラプラス法とバネ法を結合したハイブ リッド法を提案する．

図3.6にハイブリッド法の処理手順を示す．ハイブリッド法のメインフローはラプラス法 であり，図3.6の太線の矢印で示された流れとなる．ハイブリッド法は，その処理時間のほ とんどをラプラス法が占めるため，全体の処理時間を小さくすることができる．しかし，

ラプラス法だけでは局所的な要素の扁平や反転を避けることは難しい．そこで，要素の扁 平や反転が生じた局所的な部分にのみバネ法を適用しメッシュの修復を行う．バネ法が適 用される要素群は限られているため，メッシュ修復による処理時間の増加は小さい．ラプ ラス法による全体処理時間の短縮とバネ法によるメッシュ修復を組み合わせることにより，

ハイブリッド法では処理速度の短縮と，メッシュの変形に対する耐久性の両方が期待でき る．

Determine deformation using gradient field manipulation method

Update positions of grids Collapsed mesh

exists

YES Pick up

failed areas

NO Determine deformation

for local area using linear-torsion spring

analogy method

Move boundary to final position

YES

End Update matrix using

updated meshes

NO

NO

Reset displacement vectors as boundary condition

YES Start

Collapsed mesh exists Determine deformation

using gradient field manipulation method

Update positions of grids Collapsed mesh

exists

YES Pick up

failed areas

NO Determine deformation

for local area using linear-torsion spring

analogy method

Move boundary to final position

YES

End Update matrix using

updated meshes

NO

NO

Reset displacement vectors as boundary condition

YES Start

Collapsed mesh exists

図3.6 ハイブリッド法によるメッシュモーフィング処理手順

ハイブリッド法のメッシュ修復において，極端に扁平，反転した要素のみを対象とする だけでは修復はできないため，その周辺の正常に変形できた要素も対象に含める必要があ る．図3.7にメッシュ修復のための要素選択について概念図を示す．極端に扁平，反転した 要素は確実にメッシュ修復の対象として選出されるが，同時にその周辺の正常に変形でき た要素はランク付けされ，ランクの高い順に対象要素として選出されていく．ハイブリッ ド法では，どの程度周辺要素を含めるかを制御するパラメータが必要となる．この制御パ ラメータは，全体の処理時間およびメッシュ修復成功の確率を左右する重要なパラメータ となる．本章では，メッシュ修復の対象となる要素数に上限を設け，その上限はメッシュ モーフィングの対象となる全要素数に対する比率pで決定する．

×100

=

t r

N

p N (3.18)

ここで，Nrはメッシュ修復の対象となる要素数，Ntはメッシュモーフィングの対象となる 全要素数である．比率 p が大きい場合，メッシュ修復のエリアが広くなるためメッシュ修 復の成功率が上がり，全体の反復回数が減る可能性が高くなる．しかし，1回のメッシュ修 復のための処理時間が大きくなるため，全体の処理時間の増加が懸念される．一方，比率p が小さい場合，1回のメッシュ修復のための処理時間は小さくできるものの，メッシュ修復 の成功率が下がり，全体の反復回数が増加することが懸念される．つまり，メッシュモデ ルに応じて適切な制御パラメータ値を選択する必要がある．

Collapsed element

Intact element with high rank

Non-target element Intact element with low rank (2)

(3)

(1) (3)

(2) (2)

(2) (2) (2) (2)

(2) (2) (2) (2) (2)

(2) (1)

(1)

(1) (1) (1)

(1) (1) (1)

(1) (1) (1)

(1) (1) (1) (1) (1) (1) (1)

(1) (1) (1) (1) (1) (1)

(1)

Collapsed element

Intact element with high rank

Non-target element Intact element with low rank (2)

(3)

(1) (3)

(2) (2)

(2) (2) (2) (2)

(2) (2) (2) (2) (2)

(2) (1)

(1)

(1) (1) (1)

(1) (1) (1)

(1) (1) (1)

(1) (1) (1) (1) (1) (1) (1)

(1) (1) (1) (1) (1) (1)

(1)

3.3   数値解析による検証

3.3.1 メッシュモーフィング手法の特性評価

各種メッシュモーフィング手法の電磁界解析特有のメッシュモデルへの適用可能性およ び基本的な特性を評価するために，単純形状のメッシュモデルを用いて，ラプラス法，バ ネ法，およびハイブリッド法の処理速度，変形後のメッシュの全体品質を比較する．ここ では，ボール磁石の単純な形状変形を扱う．メッシュモデルには磁石とその周辺空気領域 が混在したものを用い，全てテトラ要素で分割されているものを用いる．ボール磁石はメ ッシュモーフィングにより，縦方向に0.5倍，横方向に2.0倍に変形する．なお，各種メッ シュモーフィング手法の反復解法であるICCG法の収束判定については，相対残差ノルムが 1.0e-10以下に達した時点で収束とみなした．使用計算機には，Quad-core Intel Xeon processor 5130(2.0GHz)，2.0GByte RAMを用いている．

図3.8 に，要素数10,592 のメッシュモデルをハイブリッド法で変形した場合の変形前後 のメッシュモデルを示す．ボール磁石周辺に極端に扁平した要素もなく，周辺の節点が滑 らかに移動していることがわかる．

a) Original mesh

b) Deformed mesh a) Original mesh

b) Deformed mesh

図3.8 ボール磁石モデルのメッシュ図

表3.1に，各手法の数値解析結果の比較を示す．ここでのハイブリッド法の制御パラメー タは20%とした．反復回数は図3.1および図3.6のループの繰り返し回数を表し，その回数 が少ない程処理時間は短く，安定していると言える．また，メッシュの全体品質は次式に よって計算される．

∑

=

=

_N

i i

all

Q Q N

1

1

^(3.19)

ここで，N は要素数，Q_iはテトラ要素の内接球，外接球の半径比に基づいて計算される各 要素の要素品質値[3.15]である．Qiは0.0〜1.0の値を取り，正四面体に近い程1.0に近づき，

要素品質が良いと評価する．表 3.1より，ラプラス法は他の手法に比べて反復回数が多く，

要素数が多くなる程その回数が増えており，このメッシュモデルにおいては不安定と言え る．バネ法は要素数の増加に対する反復回数の増加が少なく安定性は良いものの，1反復あ たりの処理時間が大きいため，比較した手法の中で最も処理時間を要している．一方，ハ イブリッド法は要素数によらず安定性，処理時間共に最も良い結果を示している．また，

全体のメッシュ品質についてもバネ法とほとんど変わらず，メッシュの扁平に対する耐久 性は高いと言える．

表3.1  ボール磁石モデルにおける数値解析結果 Number of

elements Laplace method Spring method Hybrid method Number of steps 5 1 1 Total CPU time (s) 1.20 2.80 2.34 10,592

Joint quality 0.50 0.53 0.53 Number of steps 9 3 1 Total CPU time (s) 13.91 37.84 12.61 84,736

Joint quality 0.20 0.45 0.43 Number of steps 21 3 1 Total CPU time (s) 328.53 377.31 193.98 677,888

Joint quality 4.12e-15 0.39 0.38

次に，ハイブリッド法の制御パラメータの特性を明らかにするために，同メッシュモデ ルを用いて制御パラメータを変更する検証を行った．図3.9に制御パラメータに対する全体 の処理時間および要素の最低品質値の関係を示す．ここでは，要素数677,888のメッシュモ デルを用いている．制御パラメータ値を小さくすると，メッシュ修復エリアが狭くなるた めに全体処理時間が少なくなるものの，修復されずに品質の悪い要素が残っていることが わかる．一方，制御パラメータ値を大きくすると，メッシュ修復エリアが広くなりメッシ ュ品質は改善されるものの，全体の処理時間が増加する傾向が見られる．全体処理時間と 要素の最低品質値の間には，概ねトレードオフの関係があり，制御パラメータにはモデル ごとに適切な値を選択する必要があることが明らかとなった．

0 50 100 150 200 250 300

5 10 15 20 25 30

Parameter value

CPU Time, sec

0.000 0.002 0.004 0.006 0.008 0.010 0.012

Minimum quality

CPU time Minimum quality

×

図3.9 ハイブリッド法の制御パラメータとCPU時間およびメッシュ品質

3.3.2 電磁界解析用メッシュモデルへの適用

実用的な電磁界解析用メッシュモデルに対する検証として，図 3.10，図3.11に示す非接 触給電装置モデル，磁気ヘッドモデル[3.16]を扱い，それぞれコイル位置の変更やヘッド先 端形状の変形を各種メッシュモーフィング手法で行う．前節の検証と同様に，各種手法で 処理速度およびメッシュの全体品質を比較する．ICCG法の収束判定値や使用計算機は同じ ものを使用する．

Z

X Y

Core

Power cable

Secondary coil Z

X Y

Core

Power cable

Secondary coil

図3.10 非接触給電装置モデル

Z

X Y Yoke Upper coil

Lower coil

Z

X Y Yoke Upper coil

Lower coil

まず，図3.10 に示す非接触給電装置モデルにおいて，コアおよび2次コイルの位置をY 軸方向に5mm，10mm，15mmと移動させた場合について検証を行う．変形量を大きくする に従い，メッシュモーフィングには高い耐久性が求められる．メッシュモデルには，要素 数860,493のテトラ要素メッシュを用いる．

図3.12，図3.13にハイブリッド法でコアおよび2次コイルの位置を15mm移動させた場 合の変形前後のメッシュモデル，表3.2に各手法の数値解析結果を示す．ハイブリッド法の 制御パラメータには 5%を用いており，このメッシュモデルでは 5%が最適値となった．移 動量の増加に伴い，処理時間はどの手法も増加しているが，ハイブリッド法はその増加量 が小さく，メッシュの品質劣化も小さいことがわかる．ハイブリッド法はメッシュ修復が 不要な場合，ラプラス法のみで動作する．そのため，移動量が 5mm，10mm の時のハイブ リッド法の処理時間はラプラス法に対して高々10 秒程度の違いとなっている．一方，移動 量が 15mm の場合はメッシュ修復が必要となるため，内部でラプラス法とバネ法が混合で 動作することになるが，反復回数が小さく抑えられ，結果として他の手法より全体処理時 間が短くなる．ハイブリッド法は状況に応じて処理手順を適切に切り替え，処理時間が最 速となるように動作していると言える．

表3.2  非接触給電装置モデルにおける数値解析結果 Displace-

ment(mm) Laplace method Spring method Hybrid method Number of steps 1 1 1 Total CPU time (s) 104.08 241.00 112.36 5

Joint quality 0.78 0.78 0.78 Number of steps 1 3 1 Total CPU time (s) 104.14 395.78 111.89 10

Joint quality 0.78 0.78 0.78 Number of steps 18 5 5 Total CPU time (s) 346.09 554.34 185.33 15

Joint quality 7.36e-6 0.77 0.77

a) Original mesh

b) Deformed mesh a) Original mesh

b) Deformed mesh

図3.12 非接触給電装置モデルのメッシュ図（モデル周辺図）

a) Original mesh

b) Deformed mesh a) Original mesh

b) Deformed mesh

図3.13 非接触給電装置モデルのメッシュ図（モデル拡大図）

次に，図3.11に示す磁気ヘッドモデルにおいて，先端形状をY軸方向に引き伸ばした場 合について検証を行う．メッシュモデルには，要素数603,114のテトラ要素メッシュを用い る．

図3.14，図3.15にハイブリッド法で変形させた場合の変形前後のメッシュモデル，表3.3 に各手法の数値解析結果を示す．ハイブリッド法の制御パラメータには 5%を用いており，

このメッシュモデルでは5%が最適値となった．このメッシュモデルでは，ラプラス法のみ が変形に失敗するという結果となった．ハイブリッド法は，バネ法より反復回数が多いも のの，全体の処理時間はバネ法より少なく，メッシュ品質もバネ法と同じとなった．

表3.3  磁気ヘッドモデルにおける数値解析結果

Laplace method Spring method Hybrid method Number of steps Failure 1 3 Total CPU time (s) - 131.60 90.98

Joint quality - 0.51 0.51

a) Original mesh

b) Deformed mesh a) Original mesh

b) Deformed mesh

図3.14 磁気ヘッドモデルのメッシュ図（モデル周辺図）

a) Original mesh

b) Deformed mesh a) Original mesh

b) Deformed mesh

図3.15 磁気ヘッドモデルのメッシュ図（モデル拡大図）

3.3.3 並列化計算の検討

各手法におけるICCG法の収束性はいずれも良好であり，数十〜数百回程度で収束する．

プログラムの各サブルーチンにおける処理時間を計測した結果，そのほとんどはICCG法の ためのマトリクスのセットアップ処理が占めているがわかった．

そこで本章では，更なる高速化を目的として，そのセットアップ処理を並列処理によっ て加速させること検討した．ここでは Intel 社が提供する並列化ライブラリ Intel Threading Building Blocks2.2(TBB)[3.17]を用いて，セットアップ処理を並列化させることを試みた．セ ットアップ処理とは，各要素単位で有限要素法に基づきローカルマトリクスを計算し，そ れをグローバルマトリクスへ足し込む処理である．ローカルマトリクスの計算は各要素で 独立して行えるため並列化による処理速度の向上が期待できる．一方，グローバルマトリ クスへの足し込みでは同じメモリ領域にアクセスする場合があるため，各スレッドが適切 に同期を取らなければいけない．TBB を用いると，その同期処理を容易に実装することが でき，またセットアップ処理中のforループの並列化についてもプログラムソースコードを ほとんど変更することなく最適な並列化を実装することができる．

本章における並列化の検討では，Quad-core Intel Xeon processor 5130(2.0GHz)，2.0GByte RAMの計算機を用いて並列処理を行った．並列計算では4コア全てを用い，並列化時のFor ループの分解にはライブラリの自動分解機能（auto_partitioner）を用いている．なお，メッ シュモデルやICCG法の収束判定，ハイブリッド法の制御パラメータについては，非並列時 と全く同じものを使用する．

表 3.4 に並列処理を行った場合のバネ法とハイブリッド法の全体処理時間と速度向上率 を示す．なお，速度向上率は，非並列時の処理時間に対する並列時の処理時間の比率で計 算される．ここでは，並列化性能が期待できるバネ法とハイブリッド法のみを対象とした．

バネ法はハイブリッド法よりセットアップ処理が全体の処理時間に占める割合が大きいた め，ハイブリッド法より高い並列性を示している．しかし，それでもほとんどの場合にお いてハイブリッド法の方が全体処理時間は小さく，ハイブリッド法の処理速度の高さが示 されている．

表3.4  並列計算による処理時間と速度向上率

Spring method Hybrid method Models

CPU Time(s) Speed up ratio CPU Time(s) Speed up ratio Ball with 10,592

elements 1.48 1.89 1.50 1.56

Ball with 84,736

elements 19.83 1.91 10.58 1.19

Ball with 677,888

elements 230.14 1.64 145.34 1.33 IPS device with

5mm disp. 164.03 1.47 118.22 0.95 IPS device with

10mm disp. 245.89 1.61 120.59 0.93 IPS device with

15mm disp. 317.98 1.74 185.11 1.00 Magnetic head 93.61 1.41 83.67 1.09

3.4 まとめ

本章では，形状最適化設計や電磁界数値解析における可動部の移動などメッシュモデル の変形作業の効率化に対する高い要求のもと，メッシュモデルの変形に対する耐久性が高 い，高速なメッシュモーフィング手法を開発し，電磁界解析特有の空気領域を含む複雑な メッシュモデルへの適用を可能とした．本章で得られた成果をまとめると，以下のように なる．

(1) 数値解析において代表的なメッシュモーフィング手法である，ラプラス方程式を支配方 程式に持つ手法（ラプラス法）とバネの釣り合い式を支配方程式に持つ手法（バネ法）

を電磁界解析特有の空気領域を含むメッシュモデルに適用し，既存手法の適用性を評価 した．数値解析による検証の結果，既存手法ではメッシュモデルの変形に伴う扁平要素 の発生を回避する耐久性と全体の処理時間の短縮を同時に満たすことが困難であるこ とを明らかにした．

(2) 耐久性の向上と処理時間短縮を目的として，既存手法であるラプラス法とバネ法を組み 合わせるハイブリッド法を開発した．ハイブリッド手法は演算量の少ないラプラス法を ベースとした手法であり，処理時間を少なくすることが期待できる．しかし，ラプラス 法だけではメッシュモデル変形時の局所的な要素の扁平や反転を回避することができ ないため，ハイブリッド法では耐久性の高いバネ法を用いて修復を行う．ハイブリッド 法の耐久性および処理時間はバネ法が対象とする修復エリアの大きさに左右されるた め，修復するエリアの要素数を全体要素数の比率から決定する制御パラメータを提案し，

その制御パラメータの特性を数値解析により明らかにした．また，ハイブリッド法を実 用的な電磁界解析用メッシュモデルに適用し，提案手法は従来手法より耐久性および全 体処理時間において高いパフォーマンスを示した．

(3) メッシュモーフィング手法全てに共通して，各種手法の支配方程式を節点有限要素法で 離散化して得られる連立一次方程式の係数行列のセットアップに処理時間の大半を要 していることを数値解析の結果から見出し，その処理時間短縮を目的として並列化を検 討した．並列化ではIntel社が提供するThreading Building Blockライブラリを使用し，4 コアによる並列計算において，バネ法で最大で 1.91 倍，ハイブリッド法で最大で 1.56 倍程度まで処理時間を短縮することを可能とした．