⽇本応⽤数理学会
三部会連携「応⽤数理セミナー」

電子機器装置の大規模変形応力解析に対する内点法の応用
富士通アドバンストテクノロジ株式会社
稲垣和久
2016年12月27日

発表概要
1. はじめに
2. 内点法による接触解析
1. 定式化と解法

2. 主双対内点法による接触解析
3. Mehrotraの予測子修正子法による接触解析
4. 増分解析におけるwarm-start法の適⽤

3. 大規模並列接触解析

1. マスタースレーブ要素による大規模並列化

4. まとめと今後の課題

1

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

はじめに

2

応用数理セミナー・Copyright 2016 FUJITSU LIMITED富士通アドバンストテクノロジ株式会社
富⼠通グループにおける先⾏共通テクノロジ・開発技術の開発を担う（私はここで）電子機器のCAEに関する技術開発に従事3

応用数理セミナー・Copyright 2016 FUJITSU LIMITEDCAE（Computer Aided Engineering）
計算機を⽤いたシミュレーションで、製品の設計、製造や工程設計の事前検討を支援することにより、下記の効果を得る
設計品質の向上

製品設計期間の短縮
コスト削減

解析
導入

4

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

電子機器のCAE
構造設計上の課題（大元の課題）
堅牢性の確保•製造、運搬、使⽤時に受ける負荷から機能部を保護する要求仕様の実現

・運搬機会の増加
・構成部品⾃体の軽薄短小化→ 構造設計のコストは年々増大•軽薄短小、低コスト•防水などの付加機能

短い開発期間︓企画〜量産約6か月（cf. ⾃動⾞︓約15か月）
携帯電話
落下衝撃や接触による
画面割れ/コネクタ破損

ノートＰＣ
背面カバーからの圧迫による
液晶ディスプレイの割れ
開閉時の負荷によるヒンジ破損
底面側からの圧迫によるチップ部品剥離

電子機器の装置障害事例
5

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

電子機器のCAE
解析における課題（技術編）

近年解析項目・対象装置の拡大により、
解析手番の問題が顕在化している
特にモバイル機器等の接触解析が課題•複雑・大規模な接触問題となるため解析時間が⻑い•部品点数が多いためモデル作成工数がかかる•圧迫評価箇所多数のためケース数が多い
圧迫解析モデル
構成の例（ノートPC）

圧迫箇所の例6

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

電子機器のCAE
解析における課題（運⽤編）
商⽤コードによる解析環境では•ライセンス費⽤が高額になる•大規模計算実⾏時の並列性能が得られない•機能不⾜への対処が難しい

課題解決のカギ = FrontISTR

並列有限要素法ライブラリhecMWを基盤とする
大規模構造解析ソフトウェア
http://www.multi.k.u-tokyo.ac.jp/FrontISTR•領域分割による高並列計算•接触解析含む充実した非線形解析機能•オープンソフトウェア…コードを⾃由にカスタマイズ可・ライセンスコスト0•開発元…東大奥田研究室・FrontISTR研究会

本日の講演では、FrontISTRを活用による上記課題の解決に向けた研究開発について紹介する
7

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

本講演のコンセプト「数理」「コンピューティング」「実⽤」の3「ある」講演
数理

コンピューティング•連続体⼒学•有限要素法（FEM）•内点法•大規模並列FEM•接触解析の並列化

実⽤①元々私は数理の⼈間（修⼠までは数理情報学専
攻で、数値積分法の収束性など研究をしていました）•電子機器向けCAE•解析ツールとして提供中②富⼠通（株）入社後、電子機器向けの大規模解析技術開発に従事③解析における課題解決に、数理専攻時代に知った内点法が使えそう…と思い、解析技術の研究開発を⾏ってきました
8

応用数理セミナー・Copyright 2016 FUJITSU LIMITED発表概要
1. はじめに
2. 内点法による接触解析

背景・基礎

1. 定式化と解法

2. 主双対内点法による接触解析
3. Mehrotraの予測子修正子法による接触解析
4. 増分解析におけるwarm-start法の適⽤

3. 大規模並列接触解析

1. マスタースレーブ要素による大規模並列化

数理
実⽤
コンピュー
ティング

4. まとめと今後の課題

9

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

内点法による接触解析

10

応用数理セミナー・Copyright 2016 FUJITSU LIMITED発表概要
1. はじめに
2. 内点法による接触解析

背景・基礎

1. 定式化と解法

2. 主双対内点法による接触解析
3. Mehrotraの予測子修正子法による接触解析
4. 増分解析におけるwarm-start法の適⽤

3. 大規模並列接触解析

1. マスタースレーブ要素による大規模並列化

数理
実⽤
コンピュー
ティング

4. まとめと今後の課題

11

応用数理セミナー・Copyright 2016 FUJITSU LIMITED基礎方程式（静解析・摩擦なし）
連続体の釣合方程式を…
div

on
on
on

⋅

,
,
,

# 釣合方程式# 変位境界条件# 荷重境界条件

接触・結合制約のもとで解く
接触境界条件⋅≡on,

結合拘束（ねじ止め等の模擬）
#

"

#

on

)

*

非貫入制約（KKT条件）
0,# 接触ギャップ
0,# 法線方向接触圧
0,!

!

︓接触ギャップ（貫入時は負値）- ギャップが0で接触圧が発生(>0)- ギャップが存在(>0)して接触圧が0のいずれかであることを意味している!
$%)&'( への投影点と変位が一致
12

応用数理セミナー・Copyright 2016 FUJITSU LIMITED基礎方程式の離散化
有限要素法による離散化

釣合方程式から仮想仕事式（弱形式）を導出
)

+,

1
3! 0

内⼒

: ./ d2

)

+,

1
3! 0

体積⼒

: . d2

表面⼒

)

+,

3! 45

1

⋅ . d6

接触⼒

)

+,

3! 47

1

⋅ . d6

各仮想仕事を、有限要素離散化したものに置き換えた式を解く
詳細な定義、導出は本講演では割愛するが、最終的に8!!
節点変位ベクトルによって以下のように近似される
8!)
8!9

!

•

内⼒ベクトル
9••

外⼒ベクトル

接触⼒ベクトル
13

)

8!!
8!)
:
8!9
⋮
8<!
8<)
節点変位ベクトル 8<9
応用数理セミナー・Copyright 2016 FUJITSU LIMITED接触の離散化
Node to Surfaceの離散化

接触面対をマスター面、スレーブ面に区分
スレーブ面の節点が、各マスター要素面に貫入しない
接触⼒と接触制約…スレーブ節点ごとに１つ
0, ,= 0, ,= ,= 0 6,=

1,2, … , AB

接触スレーブ節点の数
結合拘束も同様にスレーブ節点ごとに１つ与えられ、
スレーブ節点の変位がマスター要素面の投影点の変位と一致するよう拘束される#"
0 6 1,2, … , AC=
結合スレーブ節点の数

14

応用数理セミナー・Copyright 2016 FUJITSU LIMITEDNode to Surface接触
離散化された接触条件
接触ギャップ• D = ︓スレーブ節点座標,= ≡D = "D #

⋅

xD
=
=

• D # ︓マスター投影点座標•

D#• HI︓投影点の要素面内座標

xD #マスター要素面= ︓法線方向単位ベクトル

マスター投影点座標= ︓スレーブ節点

投影点<L+ FG HI D JK

D JM

G3!

• FG ︓要素面の形状関数（各節点の重み）• D JK ︓要素面構成節点座標

15

Q)

D JN

x HI

QR! , QR)

D JP
Q!

D JO

応用数理セミナー・Copyright 2016 FUJITSU LIMITEDNode to Surface接触
離散化された接触条件
VW

<X

,

+
=3!

S=

,= S=

S= ︓接触方向の射影ベクトル

このとき

S*= D

,=

結合拘束条件
="

⋮
"F! =
"F) =

<L

S*=

Y= ︓結合制約の係数ベクトル

16

マスタースレーブの構成節点以外の成分は0U=
初期ギャップ

JK

G3!

マスター節点⋮#I= " + FG H

Y*=

S*= T

スレーブ節点=

接触⼒ベクトル

Y=

"F!
0
0
"F)
0
0

1
0
0

0 0
1 0
0 1

0
"F!
0
0
"F)
0

スレーブ節点
0
0
"F!
0
0
"F)

マスター節点応用数理セミナー・Copyright 2016 FUJITSU LIMITED離散化された接触問題
最終的に解くべき式は次の通り
Z

<X

"V

S*=

,=

"+
=3!

U=
0,

*
,= S=
Y*=

,= S=

,

# 接触⼒を含む釣合方程式# 変位境界条件

0,6

U=0, 6 1,2, … , AB
1,2, … , AC .

# 非貫入制約# 結合制約

特に微小変形を仮定し、変位境界条件処理を⾏った場合<X\ "V"+
S*=

,=

U=
0,

*
,= S=
Y*=

=3!

,= S=

,

# 釣合方程式
K:剛性⾏列（正定値対称・疎⾏列）

0,

6

U=

0, 6 1,2, … , AB
1,2, … , AC .
17

# 非貫入制約# 結合制約
応用数理セミナー・Copyright 2016 FUJITSU LIMITED制約付き最適化問題としての接触問題
離散化された接触問題は、次の制約付き最適化問題のKKT条件になっている
1 *minimize `,`≡\ " V*
2*
S=
U= 0 6 1,2, … , AB ,Y*=
0 6 1,2, … , AC .

# ポテンシャルエネルギー# 非貫入条件# 結合拘束条件

有限要素接触解析では、上の問題に最適化問題の解法を適⽤して解を得る

18

応用数理セミナー・Copyright 2016 FUJITSU LIMITED非貫入条件の処理
非貫入制約の処理

FrontISTRでは有効制約法によって等式制約問題に帰着する•貫入ギャップが0となるべきスレーブ節点集合
ab ⊂ bdを、試⾏錯誤的に決定する•各試⾏の中では等式制約問題。
FrontISTRではLagrange乗数法
または拡張Lagrange乗数法で解く

19

\

\W e

6 ∈ abが6 ∉ abがV"ZV,= g 0→ abから外す,= i 0→ abに加える

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

電子機器接触解析の課題
積層する薄板部品のため接触範囲が広く、試⾏錯誤的な接触状態の決定が困難
3000
接
触
点
数

接触点の変化数が収束しにくい↓

1000
500

1 2

状態推定反復回数

電子機器モデルでの14 15 接触点数の推移例（⾃由度と⽐べて）多数の接触制約を含む接触問題について
効率的な解法を開発する必要がある
20

応用数理セミナー・Copyright 2016 FUJITSU LIMITED既存の接触問題の解法※実装○＝自由度数の変化無主な接触解析手法
有効接触点
の探索方針

変数

制約条件の
処理

制約

行列

実装

実用上の
課題

有効制約法
と
等式制約問
題の解法を
組み合わせ

主変数
のみ
※主変数
⇔変位，
双対変数
⇔接触力

Penalty法

×近似的

×悪条件

○容易

×制約精度
の確認要

拡張Lagrange法 △近似的・
（乗数法）
反復改良可

△Penalty次
第

○容易

×反復が増
える

変数直接消去○共有節点相
当

×煩雑

×線形のみ

主変数と Lagrange乗数法 ○厳密
双対変数

△不定値

×煩雑

主変数
のみ

不等式制約
を直接扱う

○厳密

バリア関数法

×近似的

×悪条件

○容易

×制約精度
の確認要

接触力関数／
ユーザ定義関数

×近似的

×悪条件

○容易

×汎用性無
×制約精度
の確認要

主変数と (主双対)内点法
双対変数

△近似的・
調整可

×悪条件

○容易

21

応用数理セミナー・Copyright 2016 FUJITSU LIMITED既存の接触問題の解法※実装○＝自由度数の変化無主な接触解析手法
有効接触点
の探索方針
有効制約法
と
等式制約問
題の解法を
組み合わせ

変数

制約条件の
処理

制約

主変数
Penalty法
×近似的のみ※主変数 拡張Lagrange法△近似的・⇔変位， （乗数法）
反復改良可
双対変数
○厳密
⇔接触力 変数直接消去
電子機器の接触解析には

行列

実装

実用上の
課題

×悪条件

○容易

×制約精度
の確認要

△Penalty次
第

○容易

×反復が増
える

○共有節点相
当

×煩雑

×線形のみ

接触状態を組合せ的に扱わない方法が主変数と Lagrange乗数法 ○厳密△不定値
適していると考えらえる

×煩雑

双対変数

不等式制約
を直接扱う

主変数
のみ

バリア関数法

×近似的

×悪条件

○容易

×制約精度
の確認要
・大規模問題向けの解法として注目されている

接触力関数／
×近似的
×悪条件
・手法の改良が進められている
ユーザ定義関数

×汎用性無
×制約精度
→電子機器向けの接触解析手法として検討するの確認要

主変数と (主双対)内点法
双対変数

△近似的・
調整可

22

×悪条件

○容易

○容易

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

内点法接触の研究状況
内点法による接触解析

P.W. Christensen, A. Klarbring, J.S. Pang, and N. Str omberg. Formulation and comparison ofalgorithms for frictional contact problems. International Journal for Numerical Methods inEngineering, Vol. 42, No. 1, pp. 145-173, 1998.

主双対内点法接触（弾性体と剛体）の定式化の初出。semismooth Newton法と比較し、収束性安定性ともに内点法は劣ると結論。ただし数値実験に用いたモデルは⾃由度、接触制約共に⼩さい
P. G Tanoh, Y Renard, and D Noll. Computational experience with an interior point algorithm forlarge scale contact problems. Optimization Online, 2004.
大規模（約2万DOF）モデルで主双対内点法接触の高い収束性を数値実験により確認

宮村倫司, 寒野善博, 大崎純. 内点法による接触問題の求解に関する検討. 計算⼒学講演会講演論⽂集, Vol.
2006, No. 19, pp. 69-70, nov 2006.

解法は主内点法であるが弾性体同士（node to node）の接触問題へ拡張。平板同士の接触問題で
内点法接触の有効性を確認
T. Miyamura, Y. Kanno, and M. Ohsaki. Combined interior-point method and semismooth newton
method for frictionless contact problems. International Journal for Numerical Methods inEngineering, Vol. 81, No. 6, pp. 701-727, 2010.
途中から有効制約法にシフトして収束性を改善する提案

23

応用数理セミナー・Copyright 2016 FUJITSU LIMITED内点法の研究状況
内点法（収束改善に関する議論）
主双対内点法の提案

M. Kojima, S. Mizuno, and A. Yoshise. Progress in mathematical programming interior-point andrelated methods. chapter A Primal-dual Interior Point Algorithm for Linear Programming, pp. 2947. Springer-Verlag New York, Inc., New York, NY, USA, 1988.
線形計画問題に対する主双対内点法の提案
収束の高速化S. Mehrotra. On the implementation of a primal-dual interior point method. SIAM Journal on
optimization, Vol. 2, No. 4, pp. 575-601, 1992.
予測子・修正子法の提案。主双対内点法より実際の問題を早く解けるとして、多くの商用コードで採用される。

良い初期解の利⽤
H. Y. Benson and D. F. Shanno. An exact primal-dual penalty method approach to warmstartinginterior-point methods for linear programming. Comput. Optim. Appl., Vol. 38, No. 3, pp. 371399, 2007.
H.Y. Benson and D.F. Shanno. Interior-point methods for nonconvex nonlinear programming:
regularization and warmstarts. Computational Optimization and Applications, Vol. 40, No. 2, pp.
143-189, 2008
良い初期解を利用して反復回数を削減する

24

応用数理セミナー・Copyright 2016 FUJITSU LIMITED本研究の位置づけ（接触解法）
内点法接触の電子機器接触への拡張と収束性改善の提案
接触

収束性

線形解析
Node to Rigid接触
Node to Node接触

低高Christensen,1998
Tanoh,2006（主双対内点法）

Miyamura,2010
有効制約シフト線形解析
本講演2.2
本講演2.3
Node to Surface接触
予測子修正子法(Mehrotra, 1992)による高速化
結合接触の考慮＋電子機器で検証
大変形・増分解析
Node to Surface接触
結合接触の考慮＋電子機器で検証

本講演2.4warm-start法(Benson 2008)による高速化

25

応用数理セミナー・Copyright 2016 FUJITSU LIMITED発表概要
1. はじめに
2. 内点法による接触解析

背景・基礎

1. 定式化と解法

2. 主双対内点法による接触解析
3. Mehrotraの予測子修正子法による接触解析
4. 増分解析におけるwarm-start法の適⽤

3. 大規模並列接触解析

1. マスタースレーブ要素による大規模並列化

数理
実⽤
コンピュー
ティング

4. まとめと今後の課題

26

応用数理セミナー・Copyright 2016 FUJITSU LIMITED内点法(1/2)
以下の不等式制約付き最適化問題を解く

基本的な考え方

対数バリア関数の導入•最適点を制約境界から遠ざける• j → 0とすると を最小化する問題の解は元問題の解に収束する

27

応用数理セミナー・Copyright 2016 FUJITSU LIMITED内点法(2/2)
対数バリア関数のイメージ

n

基本戦略

次の手続きを反復し、元問題の最適値を得る
1.lm のNewton方向に進む 2. μを小さくする
1.で線形方程式を解く部分の負荷が最も高い
28

応用数理セミナー・Copyright 2016 FUJITSU LIMITED主双対内点法による接触解析
接触問題の支配方程式（線形弾性体・摩擦無微小滑り）

結合拘束をPenalty法で無制約化し、主双対内点法を適⽤する

対応︓

29

応用数理セミナー・Copyright 2016 FUJITSU LIMITED主双対内点法による接触解析
バリア問題の１次の最適条件…(1)

ここでo$ は 1)式を満たす変数として定義される•接触⼒,$ に対応する量•主双対内点法では独⽴変数として扱う

上の系に対して「Newton方向 Δ , Δq の計算」と「jの縮小」処理を反復する

30

応用数理セミナー・Copyright 2016 FUJITSU LIMITED主双対内点法による接触解析
Newton方向の計算

通常のNewton法の手続きと同じで、 , qの各成分で微分する

31

応用数理セミナー・Copyright 2016 FUJITSU LIMITED主双対内点法による接触解析
Newton方向の計算（プログラミング的な話）

FrontISTRで実装する場合、2個目の式のΔo$ について解くことにより、縮約した形にすると簡単にできる• FrontISTRは、ペナルティ法由来の接触非ゼロ項 S$ S*$ を剛性⾏列（ブロック疎⾏列）に⾜しこむ機能を持っている• 一方、内点法の接線⾏列を縮約すると⾏列がペナルティ法と同じ非ゼロ配置になる
ここを とすると、すべての接触制約がactiveの場合のペナルティ法のマトリクスと一致する• 実装の際は、接触非ゼロ項の係数を内点法のものに差し替えるだけでよい（⾏列構造を再構築しなくてよいので、とても楽）32

応用数理セミナー・Copyright 2016 FUJITSU LIMITED主双対内点法による接触解析
解析手続き

初期値設定
１．縮約した⾏列を解いてNewton方向を求める（最も負荷が高い）

２．内点に留まるため減速︓ステップ幅を{倍３．μを|倍縮小する（但し減速率に応じて調整）︓j ← {|1"{ j直接法ソルバを⽤いる場合、メインループの反復回数でほぼ処理時間が決まる

33

応用数理セミナー・Copyright 2016 FUJITSU LIMITED数値実験
主双対内点法(IPM)と有効制約法＆Penalty法(ASM)による解析の手続きをFrontISTRに実装
数値解（変形量）及びステップ数を⽐較
ステップ数は各手法におけるメインループの反復数

⽐較条件
初期

1:σ

μの更新式︓

Newton方向• 減少率の設定︓結合ペナルティ︓1E+6
接触ペナルティ（有効制約法）︓1E+7•段階的に値を大きく取って必要な値を決定した

34

σ：内点制約を満たすた
めの減速の割合
η：μの減少率（入力パラメータ）

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

解析モデル(1/3)
検証⽤小規模モデル２点

モデル１︓1,215⾃由度，接触点63，結合点9モデル２︓11,283⾃由度，接触点909，結合点205

35

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

解析モデル(2/3)
スマートフォンケースモデル

142,596⾃由度，接触点2,563，結合点535
曲げ試験時のパネル部でのひずみ評価が目的
勘合部での接触状態が結果に影響する

36

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

解析モデル(3/3)
ノートパソコン液晶ガラス背面圧迫解析

523,426⾃由度，接触点50,658，結合点2,492
内部の液晶ディスプレイに規定値以上の応⼒が発生しないか⾒る

37

応用数理セミナー・Copyright 2016 FUJITSU LIMITED解析結果(1/4)
変形（鉛直方向変位）

内点法(IPM)と有効制約法(ASM)で一致を確認

0.14

-1.5

0.43

-4.24

モデル3,4の鉛直方向変位コンター：左・有効制約法、右・内点法
38

応用数理セミナー・Copyright 2016 FUJITSU LIMITED解析結果(2/4)
反復回数

内点法(IPM)の反復回数は|に依存•ηが小さいほど収束は速いが、モデル4はη=0.5を境に収束性低下

単純なモデル1を除いては、最良のパラメータを⽤いた場合で約20％から50％程度反復回数減
iter

iter

60
40

48
38

22

21

14
7
η

※ASMはηに依らないη
39

応用数理セミナー・Copyright 2016 FUJITSU LIMITED解析結果(3/4)
バリアパラメータμの収束履歴(1/2)
簡易モデルの場合μは等⽐級数的に減少→ほぼσ=1であり、より大きなηを取れる

モデル1

モデル2

40

応用数理セミナー・Copyright 2016 FUJITSU LIMITED解析結果(4/4)
バリアパラメータμの収束履歴(2/2)
複雑なモデルでは、ηを過剰に小さくすると（減速によって）μの減少が緩和する→適切なηの設定が課題

モデル3

モデル4

41

応用数理セミナー・Copyright 2016 FUJITSU LIMITED内点法はどのようなモデルで有利か
荷重を徐々に変えて接触範囲を変化させ、内点法（η=0.3）と有効制約法の反復回数を⽐較
接触対の個数※に対し、内点法はほぼ一定の反復回数
接触範囲が広い（本モデルでは約200以上で）場合に内点法が有利となる
反復回数

個数
42

※（解析事後にしか計算できない量であるが）
接触の広さの指標としてASMで変化した接触対数を採用した
応用数理セミナー・Copyright 2016 FUJITSU LIMITEDまとめ
Tanohらによる内点法接触を拡張、FrontISTRに実装した結合拘束を考慮

弾性体同⼠のnode to surface接触

電子機器モデルによる数値検証を実施

実機モデルにて有効制約法より約20％〜50％高い収束性を確認．
収束性はパラメータ|に大きく依存
適切なパラメータ設定が必要になる点が課題

接触範囲の広さを変えた場合の収束性について
モデル2で内点法と有効制約法の収束性能の分岐点（約200程度の接触点数）を確認．

43

応用数理セミナー・Copyright 2016 FUJITSU LIMITED発表概要
1. はじめに
2. 内点法による接触解析

背景・基礎

1. 定式化と解法

2. 主双対内点法による接触解析
3. Mehrotraの予測子修正子法による接触解析
4. 増分解析におけるwarm-start法の適⽤

3. 大規模並列接触解析

1. マスタースレーブ要素による大規模並列化

数理
実⽤
コンピュー
ティング

4. まとめと今後の課題

44

応用数理セミナー・Copyright 2016 FUJITSU LIMITED予測子修正子法(Mehtotra, 1989)
予測・修正の2ステップを反復し元問題の最適解を探索【予測ステップ】
・ j 0 としてNewton方向 ∆•, ∆q€ を計算する【jの更新】
・予測ステップでの最大ステップ幅•€ , •€‚ を求める

・予測前後の双対ギャップの⽐を⽤いてjを更新する
j…(†‡•€ ∆•, q,q•€‚ ∆q€ˆj‰Š‹ ここで【修正ステップ】
・主双対内点法における１次の最適条件に,q∆• ∆ƒ, q∆q€ ∆q„ を代入、修正ステップ ∆ƒ, ∆q„ について解く

,q

<7

+ S*$
$3!

U$ o$

(2)修正：中心パスに戻る(1)予測：
μを減らす※修正ステップ量が微小と仮定して得られる線形方程式を解く

45

応用数理セミナー・Copyright 2016 FUJITSU LIMITED予測子修正子法の特徴
1ステップあたりで、⾏列ソルバの実⾏が実質1回で済む

予測ステップと修正ステップの線形方程式について
係数⾏列が一致する（ように手法が設計されている）ため【予測ステップの線形方程式】\<Œ+ Y$ Y*$

S! S• ⋯ S<7

o! S!‘
⋮o<7 S‘<7

•!

$3!

•

•)

Δ•

•

∆q€

•<7

"/

,q

"o! •!
⋮"o<7 •<7【修正ステップの線形方程式】\<Œ+ Y$ Y*$$3!

o! S!‘
⋮o<7 S‘<7

S! S• ⋯ S<7
•!
•

•)

Δƒ

•

∆q„

•<7
46

j " ∆o•! S!‘ ∆•⋮j " ∆o̅<X S‘<7 ∆•
応用数理セミナー・Copyright 2016 FUJITSU LIMITED予測子修正子法の流れ
パラメータ設定 j, “! , “)
, q を与える” 0
Loop ”
・予測ステップ方向 ∆•, ∆o̅$ を計算する
・jの更新
・修正ステップ方向 ∆ƒ, ∆õ$ を計算する・ステップ方向を次のようにとる
Δ , Δq∆• ∆ƒ, ∆q€ ∆q„・内点を保つようステップの減速率{! , {) ∈ 0,1– を決める
・変位を更新する{! ∆ , q —! q
{) ∆q
—!
” ” 1
If 残差i “! and j i ϵ) then stop
Continue loop47

線形方程式を解く
LU分解結果を保存しておく
LU分解結果を流⽤して高速に解く

１反復あたりの計算負荷は増加するが反復回数が減少してトータルの計算時間が減少すればよいと考える

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

数値実験
予測子修正子法(PCM)、主双対内点法(PDIPM)、有効制約法＆Penalty法(ASM)による解析の手続きをFrontISTRに実装
数値解（変形量）及び計算時間、ステップ数を⽐較
ステップ数は各手法におけるメインループの反復数

⽐較条件
初期

μの減少率
結合ペナルティ︓1E+6
接触ペナルティ（有効制約法）︓1E+7

計算環境

CPU︓Intel® Xeon ® CPU X5675 @ 3.07GHz,12core
メインメモリ︓72GB

48

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

解析モデル(追加)ノートパソコンキーボード圧迫解析（相当）キー押下（ここでは面荷重で代⽤）時の沈み量が設定した許容値内に収まるかを⾒る

2,214,459⾃由度、接触候補対数41,039、
結合拘束数8,923
49

応用数理セミナー・Copyright 2016 FUJITSU LIMITED解析結果(1/3)
変形量は各手法とも一致することを確認

モデル5荷重方向変位コンター（左：有効制約法、中：主双対内点法、右：予測子修正子法）
50

応用数理セミナー・Copyright 2016 FUJITSU LIMITED解析結果(2/3)
反復回数

⾃由度の小さいモデル1を除けば予測子修正子法＜主双対内点法＜有効制約法の順
予測子修正子法は主双対内点法の約50％程度かそれ以下
3手法の反復回数比較

PDIPM
η
0.01
0.1
0.3
7
10
1
11
12
40
15
2
18
22
36
19
3
20
20
60
4
50
49
47
135
5
62
60
57
赤…1位、青…2位、黒…3位
ASM

PCM
0.5
20
34
23
44
56

0.7
38
64
39
49
67

5
3
13
25
25

主双対内点法(PDIPM)は最も反復回数の少ないηのケースを採用
51

応用数理セミナー・Copyright 2016 FUJITSU LIMITED解析結果(3/3)
実計算時間

予測子修正子法が反復あたりの計算時間が最も大きいが、他の2 つに⽐べてその
増加は軽微

6000

0

52

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

まとめ
内点法による接触解析を高速化する方法として
Mehrotraの予測子修正子の検討を⾏った
アルゴリズムをFrontISTRに実装
電子機器モデルによる数値検証にて標準的な主双対内点法による接触解析より40％以上実計算時間が短縮されることを確認した

53

応用数理セミナー・Copyright 2016 FUJITSU LIMITED発表概要
1. はじめに
2. 内点法による接触解析

背景・基礎

1. 定式化と解法

2. 主双対内点法による接触解析
3. Mehrotraの予測子修正子法による接触解析
4. 増分解析におけるwarm-start法の適⽤

3. 大規模並列接触解析

1. マスタースレーブ要素による大規模並列化

数理
実⽤
コンピュー
ティング

4. まとめと今後の課題

54

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

内点法接触の増分解析への拡張
増分解析の必要性
モデル3のスマートフォン3点曲げなどは、荷重を段階的に加えたときの変位履歴を取るため、外⼒の大きさを徐々に増大させる一連の解析を実⾏する

内点法接触の場合
標準的な主双対内点法は、近似解を利⽤して反復回数を減らすことができない（２ステップ目以降は実質解き直し）
ウォームスタート法を導入し、増分解析において反復回数を減らす⽅法について検討した
55

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

増分解析の場合の接触問題
全部で = %™ ステップからなる解析の、= %™ ステップ目に解く釣合方程式（微小変形の場合）
\

"V

S*=

U=

,=

0,

*
S
=
,=

Y*=

<X

"+
=3!

,= S=

,

0,
6

U=

0, 6 1,2, … , AB
1,2, … , AC .

（参考）一般的なV の与え方”
⋅VV
F= %™ š‰…›&目標の荷重値Vš‰…›& まで⽐例で増える

最適化問題

1 *
minimize `,`≡\ " V*
2
S*=
U= 0 6 1,2, … , AB ,Y*=
6 1,2, … , AC .
56

# ポテンシャルエネルギー# 非貫入条件# 結合拘束条件

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

内点法による増分解析
内点法のログバリア関数

A•

min l
žœ

0 œ

j

1,2, … , A• ,ログバリア関数

増分解析での内点法接触
Loop ”

1, F= %™
Ÿ!
,j

œ 1

œ

2ステップ目、jがリセットされるとactive（žJ ¡ 0 ）な接触ペアで過大な反発が生じる

jLoop
If 残差とjが0に収束 then stopのNewton方向を計算
1. lm
2. 接触制約を破らないように減速した上で変数を更新する
3. jを縮小する

57

応用数理セミナー・Copyright 2016 FUJITSU LIMITED内点法のウォームスタート(H. Y. Benson et. al., 2005)
ペナルティ項を付加した緩和問題を考える
min l
žœ

0 œ

min l
žœ

œ

1,2, … , A• ,+ ¢œ Qœ

œ 1

œ

0 œ

œ 1

œ

¢J ︓緩和のためのペナルティ項¤J ︓ギャップ関数の上限
QJ ︓緩和のための制御変数F¢"Qœ ≤ žœ ≤ ¤œ œ
Qœ

j

A•

1,2, … , A•

制約境界がžJ 0からžJ QJ 0にずれるため、前ステップの収束解が制約境界上に乗らない

1,2, … , A•

1,2, … , A•

QJ 0のとき、元問題と（上限¤J を除いて）一致
ペナルティ項¢J QJ がQJ を0に近づけるように（緩和問題の解が元問題の解に近づくように）働く
58

応用数理セミナー・Copyright 2016 FUJITSU LIMITED内点法のウォームスタート(H. Y. Benson et. al., 2005)
ペナルティ項を付加した緩和問題を考える
min l
žœ

0 œ

min l
žœ

œ

F¢

0 œ

lRj

+ ¢œ Qœ

œ 1"Qœ ≤ žœ ≤ ¤œ œQœj1,2, … , A• ,

œ

A•
œ 1

≡l
A•

1,2, … , A•

"j + log žœœ 1

1,2, … , A•

œ

F¢

A•

œ 1

œ 1

+ ¢œ Qœ " j + log QœA•Qœ " j + log ¤œ " žœ
œ 1↑のバリア関数を⽤いて内点法の手続きを⾏う

1,2, … , A•

59

応用数理セミナー・Copyright 2016 FUJITSU LIMITED内点法のウォームスタート(H. Y. Benson et. al., 2005)
の停留条件

m

7

J

J

J

J3!
J

J

§J は§J©J は©J

J

J

J

J

J
J

J

J
B

B

j/ ¤J " žJ をみたす変数として、mm"を満たす変数としてª«—¬«-« Ÿª«

新たに定義し、独⽴変数として扱う
Newton方向の導出は上式のヤコビアンを計算する

60

B

H

独⽴変数 ®

¯（頑張って計算すると） Δ に関する縮約された係数⾏列が得られる
応用数理セミナー・Copyright 2016 FUJITSU LIMITEDペナルティパラメータの導入
ペナルティ J

J の設定

十分な大きさの値を設定しないと解に緩和の影響が現れる懸念がある
以下の規則に従って動的にペナルティ数を調整する

初期ペナルティは本研究では以下の通りに取る• °$ は初期ペナルティを調整するパラメータ61

応用数理セミナー・Copyright 2016 FUJITSU LIMITED計算コストについて
ウォームスタートの有/無で、反復あたりの計算コストはほとんど変わらない
内点法接触のホットスポット︓
Newton方向を得るための線形ソルバ実⾏

線形ソルバは基本的に直接法→計算コストは非零成分の配置で決まる

Loop ”

1, A8œ_6 ²³
Ÿ! , jj

Loop
If 残差とjが0に収束 then stop
1. `´のNewton方向を計算
2. 変数を更新する
3. jを縮小する

Newton法の１次近似式において、主変数以外の変数は縮約でき、どちらの手法もソルバに渡す非零配置は同じになる（ペナルティ法などと同じ非零配置）

62

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

数値実験
標準的な主双対内点法（IPM）と内点法warm-start（WS）について結果の精度とステップ数を⽐較
ステップ数は各手法におけるメインループの反復数

⽐較条件
モデル3
初期
ペナルティパラメータ︓
収束閾値︓結合ペナルティ︓1E+663

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

数値実験
解析結果

warm-startあり(WS)、なし(IPM)で変位の一致を確認した

0.344

-3.02
ステップ5変位コンター：左…warm-startあり、右…なし

64

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

数値実験
反復回数

ウォームスタート導入によって約30％程度の反復回数減

総反復回数︓
warm-startなし…95回
warm-startあり…65回

65

応用数理セミナー・Copyright 2016 FUJITSU LIMITEDペナルティパラメータの検証
ペナルティパラメータ $ および緩和項 の初期値を変えた場合の収束性の違いについて検証$ =1,10,100と初期 =0.001,0.01,0.1

初期緩和小

初期緩和大

緩和項を大きく取ると収束に反復回数が費やされる
ペナルティは10が今回最良であったが、依存性が強い
適切な設定手法が今後の課題

66

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

まとめ
接触問題の増分解析に対して，warm-startによって内点法接触の収束性を向上させる方法を検討した．
数値検証ではスマートフォンモデルで標準的な主双対内点法との⽐較を実施，同程度の精度を保ちつつ反復回数を30%程度削減できることを確認した．
一方warm-start法の反復回数はペナルティパラメータ緩和項に依存する
適切な設定（可能であれば汎⽤的な設定）の検討が課題

67

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

大規模並列接触解析

68

応用数理セミナー・Copyright 2016 FUJITSU LIMITED発表概要
1. はじめに
2. 内点法による接触解析

背景・基礎

1. 定式化と解法

2. 主双対内点法による接触解析
3. Mehrotraの予測子修正子法による接触解析
4. 増分解析におけるwarm-start法の適⽤

3. 大規模並列接触解析

1. マスタースレーブ要素による大規模並列化

数理
実⽤
コンピュー
ティング

4. まとめと今後の課題

69

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

電子機器接触解析の課題
問題が大規模（数十万〜数千万⾃由度）
考慮すべき部品点数が多い

解析の高精細化、形状簡易化の工数削減のため

サーバ自重反り解析モデル（28.7M DOF）

ノートＰＣモデル（2.2M DOF）

並列計算機の性能を有効に活⽤する必要がある

システム

node

core/node

メモリ/node

東大FX10

1440

16

32GB

Intel クラスタ

18

12〜18

128GB~

70

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

有限要素構造解析の並列化
領域分割による並列計算（FrontISTRを例に）
全体領域メッシュを小領域に分割し、各プロセスが境界節点の変位のみを通信して小領域の計算をする

始めに領域分割ツール（パーティショナ）を実⾏し、分散領域メッシュを作成→各プロセスが個別に読み込む

領域分割イメージ

71

応用数理セミナー・Copyright 2016 FUJITSU LIMITED既存の並列化手法
一般的な並列接触におけるデータ管理

接触点を探索し、剛性や接触⼒を領域間で共有する

1. slave節点︓⾃領域だけ，master要素面︓全体
2. master要素面︓⾃領域だけ，slave節点︓全体（FrontISTRの並列接触はこの仕様）
3. slave 節点・master 要素面ともに⾃領域だけ

1.2…電子機器モデルのように広い接触面を持つ場合に不適
3…master面の有無を多領域と通信して取得する必要がある

72

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

電子機器モデルでの接触
接触点数は多いが、ほとんど「滑らない」
探索を事前に済ませ、近いペアで固定の通信テーブルを作れば効率的固定要素（マスタースレーブ要素）
での並列計算を提案

ノートＰＣの接触対可視化例（グレーの要素が接触対）73

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

マスタースレーブ要素による大規模並列化
電子機器モデルの幾何的特徴
薄板同⼠の広い接触面を持つ滑りをはじめ、面内方向の相対変位は少ない→ 接触可能性のある対は初期形状でほぼ決まっている

接触面対を「節点集合と面集合の組」で定義せず、事前接触探索で発⾒した「節点-要素面対」からなる要素として定義、連続体ごと領域分割する

74

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

マスタースレーブ要素による大規模並列化
マスタースレーブ要素

slave節点-master要素面の接続関係を与える要素
FrontISTRでは連続体要素と区別なく領域分割できる

個々の接触対によるコネクティビティ（＝非ゼロ成分の配置）はマスタースレーブ要素のそれと一緒
違いは固定か、そうでないか75

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

マスタースレーブ要素の利点・欠点
利点

接触関係が分割時に既知なので、一体の構造体と同様に領域分割に基づく高スケールな並列計算ができる
解析中に接触探索の時間がかからない
プログラムの実装が容易•接触剛性・反⼒の計算は連続体要素の剛性・内⼒計算ルーチンを接触計算におきかえ
ればよい

欠点

滑りの多い問題には対応できない•指定した個数内で隣接する要素も接続できるので、ある程度の滑りは許容できる（右図）

76

応用数理セミナー・Copyright 2016 FUJITSU LIMITED大規模サーバユニットモデル
サーバユニットの⾃重反り解析

ラックのガイドレールにユニットを設置した際の⾃重による撓み量を評価する
⾃由度28,680,309
接触制約646,603
結合拘束208,70277

応用数理セミナー・Copyright 2016 FUJITSU LIMITED大規模サーバユニットモデルでの検討
検討1

反復法線形方程式ソルバと拡張Lagrange乗数法の組み合わせによる
高並列計算
高並列計算時のスケール性能を確認

検討2

内点法vs有効制約法の収束性⽐較
内点法では線形方程式ソルバに直接法を⽤いるため、スケール性能ではなく収束性の⽐較
•実装は一応マスタースレーブ要素で⾏っている

78

応用数理セミナー・Copyright 2016 FUJITSU LIMITED検討1：拡張Lagrange法と反復法
拡張Lagrange法

ペナルティ数1.0，反復法収束閾値1.0E-3，
Newton 法収束閾値1.0E-6．
判定条件•接触状態の変化が接触点数の3%以内•貫入量が1.0E-6 以下 かつ 引張⼒が1.0E-6 以下反復法の設定•AMG前処理付きCG法•収束閾値1.0E-3、反復上限5000回↑未収束であっても計算継続し、外側反復を回す

計算環境

FX10（SPARC64TM IXfx:1.650GHz,16core/node使⽤）…並列数:512,1024,2048
Intel クラスタ( Xeon E5640-2.26GHz, 12core/node使⽤）…並列数:96ノード ※内点法接触と同じ
79

応用数理セミナー・Copyright 2016 FUJITSU LIMITED検討1：拡張Lagrange法と反復法
結果Intelマシンでの内点法接触(4.3H)と⽐べ約5.5倍の計算時間
スケール性能は⽐較的良好（並列化効率70%以上）であり、並列数を増やして高速化・大規模化することは可能
CG法
総計算
並列化効
並列化効
クラスタ
並列数
増速率
1反復時 増速率
時間
率
率
間
512
27.5 473.0 0.924 0.175
512
114.8 878.8 0.858 0.0875
FX10 1,024
1024
1
2,048
9.05 1437.2 0.702 0.0581
1542 0.753
intel
96
23.8
- 0.182
FX10の増速率はアムダール則から算出した
1コア時間からの⽐

サーバユニットモデル
変位Magnitudeコンター
80

応用数理セミナー・Copyright 2016 FUJITSU LIMITED検討2：内点法接触
主双対内点法と有効制約法＆Penalty法による解析の手続きを並列接触FrontISTRに実装
数値解（変形量）及びステップ数を⽐較
⽐較条件
初期

μの減少率
•

結合ペナルティ︓1E+6
接触ペナルティ（有効制約法）︓1E+7•段階的に値を大きく取って必要な値を決定した

計算環境

Xeon E5640-2.26GHz 6コア×2CPU,128GBRAMマシン×8ノード
ノード内スレッド並列計算

81

応用数理セミナー・Copyright 2016 FUJITSU LIMITED検討2：内点法接触
有効制約法はチャタリング発生し解を得られず
内点法…拡張Lagrange法と⽐べ若⼲の誤差

収束閾値が甘い可能性があるが、これ以上厳しいと収束解を得られない

100,000

接触制
約数
10,000
0
反復回数

148

82

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

まとめ
電子機器装置解析向けにマスタースレーブ要素による並列接触解析を検討．

マスタースレーブ要素生成と本手法に基づく並列接触を実⾏するモジュールを作成
実際に大規模な電子機器モデルとして約2,800 万⾃由度のサーバユニットモデルを⽤いて性能を確認した．
•反復法＋拡張Lagrange法接触の性能検証•内点法接触の性能検証

反復法＋拡張Lagrange法接触

並列化効率は70%以上得られた一方で、サーバモデルでは内点法接触と⽐べ約5.5倍の計算時間を要した

内点法と有効制約法の性能⽐較

有効制約法でチャタリングの発生する大規模モデルにおいても内点法は収束解を得られることが分かった
ただし正確な解を得るには収束閾値を厳しくしても収束するよう対応する必要がある．
83

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

まとめ

84

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

まとめ
電子機器構造解析の課題に対し，
接触解析手法
並列化手法の両面から計算効率化について検討した結果を紹介した．

産業上の実問題に対する、数理・コンピューティングの適⽤例として、
ご参考になれば幸いです85

応用数理セミナー・Copyright 2016 FUJITSU LIMITEDまとめのおまけ
少しだけ宣伝弊社では、FrontISTRに電子機器向け
各種機能追加を⾏って活⽤しています

静解析
︓機械的荷重静解析
︓熱荷重•特に接触解析を重点的に改良

プリント基板反り（2013計算工学講演会）

弊社カスタマイズ版FrontISTR「FS-Solver」＋独⾃統合解析環境「FS-Hub」を提供中ですノートPC圧迫・他（2014 OpenCAE学会）•どうぞよろしくお願い致します。

はんだ接合部疲労寿命（2016計算⼒学講演会）

FrontISTR本体開発も⾏っています•回転・トルク荷重機能の開発•⾃動増分・カットバック機能開発（予定）

FS-Hub(Pre)

FS-Solver(Solver)

FS-Hub(Post)＋カスタム部

86

Copyright 2016 FUJITSU LIMITED今後の課題接触解法についてwarm-start法のペナルティ依存性解消に向けた検討
warm-start法と予測子修正子法の併⽤による高速化検討
摩擦接触問題への拡張

87

応用数理セミナー・Copyright 2016 FUJITSU LIMITED

88

応用数理セミナー・Copyright 2016 FUJITSU LIMITED