均質化設計法

均質化設計法
藤井大地
（東京大学）
位相最適化
？
従来の考え方
T
境界形状を変化させて最適な形状・位相を
求める。
ΓT
t
Ω
ΓD
b
境界形状を変化させる問題点
T
T
解析が進むにつれて，有限要素メッシュが
異形になり，再メッシュが必要になる。
位相が変化する問題への適応が難しい。
ΓT
ΓT
t
t
Ω
ΓD
b
Ω
ΓD
b
領域の拡張と特性関数の導入
1
χΩ ( x) = ®
¯0
if
if
x∈Ω
x∉Ω
1985 Muratら
χΩ = 1
χΩ = 0
t
ΓT
Ω
ΓD
b
D
材料のON/OFF問題
整数計画問題
微分不可能問題
位相最適化の手法
T
連続緩和法（整数条件の緩和）
T
T
T
均質化設計法（ミクロ的な材料のON/OFF）
密度法（材料定数が密度のべき乗に比例）
離散最適化問題として解く方法
T
T
T
遺伝的アルゴリズムを利用した方法
セル・オートマトンを利用した方法
ESO法（Evolutionary Structural Optimization)
均質化設計法(HDM)
マクロ構造
ミクロ構造
x2
a
b
x1
θ
1
y2
y1
Unit cell
質量制約条件下で外力の仕事量の最小化
均質化設計法の目的関数
T
従来の最適化問題
T
T
応力の制約条件の下での重量の最小化
均質化設計法
T
T
重量制約条件の下でのコンプライアンス（外力
の仕事量）の最小化
さらに一般的には，最小ポテンシャルエネル
ギーの最大化
なぜ応力の制約を用いない ?
T
連続体ではいたるところで応力は発散 !
境界角点
集中荷重点
支持点
なぜコンプライアンスの最小化 ?
T
T
T
T
コンプライアンスを最小化することで構造
物の剛性が最大化される。
応力は局所的に発散してもコンプライアン
スは有限値となる。
コンプライアンスと応力との間には関係が
あり，コンプライアンスが最小化されれば，
最大応力も抑えられる。
感度解析が非常に容易になる。
最小ポテンシャルエネルギーの
最大化の方がより一般的
T
T
T
外力が与えられる問題では，コンプライア
ンスの最小化
変位が与えられる問題では，コンプライア
ンスの最大化
外力が与えられる問題も，変位が与えられ
る問題も，最小ポテンシャルエネルギーの
最大化問題となる。
均質化設計法の最適化問題
ª
1
T
ª
ºº
T
max « min «Π ( v 0 ) = ³ ( ∂ x v 0 ) D H ( α )( ∂ x v 0 ) dD − ³ v 0 t d Γ » »
ΓT
α∈L
2 D
¼¼
¬ v0 ∈K ¬
{
}
K = {v0 v0 = g
L = α 0 ≤ a ≤ 1, 0 ≤ b ≤ 1, ³ ρ dD ≤ mS ,
D
1
D
ΓT
1
D
θ
b
a
y2
x2
ΓD
x1
Global macro coordinates
Unit cell
y1
Local micro coordinates
on Γ D }
均質化された弾性マトリックス
D H ( a, b,θ ) = R (θ ) D H ( a, b ) R (θ )
T
T
D H ( a, b ) = min ª« ³ ( I − ∂ y χ ) DE ( I − ∂ y χ ) dY º»
χ
¬Y
¼
periodic
ª cos 2 θ sin 2 θ cos θ sin θ º
«
»
R (θ ) = « sin 2 θ cos 2 θ − cos θ sin θ »
« − sin 2θ sin 2θ
cos 2θ »¼
¬
1
D
ΓT
1
D
θ
b
a
y2
x2
ΓD
x1
Global macro coordinates
Unit cell
y1
Local micro coordinates
χ ：ミクロ構造の特性変位関数
有限要素法による離散化
min ª¬dT K 0d º¼ ,
α∈L
{
L = α 0 ≤ a ≤ 1, 0 ≤ b ≤ 1,
要素剛性マトリックス
K 0 e = Te
T
³
De
³
D
ρ dD ≤ mS
}
x2
T
B 0 D H B 0 dDe Te
x1
a
b
y2
1
θ
y1
均質化弾性マトリックスの離散化
DH = ³ DdY − XT K1X
Y
K 1e X e = q e
K 1e = ³ B DB1dYe
Ye
T
1
q e = ³ B DdYe
Ye
T
1
Χ ij (4)
Χ ij (3)
Χ ij (1)
Χ ij (2)
i=1,2
j=1,2,3
Finite element
χ = N1X e
∂ y χ = B1X e
1
1
ユニットセル
計算の効率化
T
均質化弾性マトリックスDHのデータベース化
T
T
T
穴の大きさa,bに関してデータベースを作成する。
データベースの値を利用して，その間の値は補
間関数を用いて補間する。
ユニットセルの角度θは，その要素の応力
の主軸方向に一致させる。
均質化設計法の解
グレースケール
チェッカーボード
フィルタリング法
T
境界線長さの制約法（Perimeter control）
T
T
T
チェッカーボードを防ぎ，シンプルな位相形状を
求めることができる。
グレースケールを防ぐ処理を付加する必要が
ある。
重力制御法（gravity control）
T
チェッカーボードもグレースケールも同時に防ぎ，
シンプルな位相形状を求めることができる。
フィルタリングの効果
重力制御関数
N
G = ¦ gi
i =1
mi
¦ m = ¦¦ ( ρ ⋅ ρ
N
i =1
N
i
i =1 j =1
i
j
+ ρi ⋅ ρ j )
N
¦m
i =1
i
ρi = 1 − ai bi , ρi = 1 − ρi , 0 < G ≤ 1
i
mi = 4
gi = 4
gi = 3
gi = 2
gi = 1
gi = 0 gi = 0.25
gi = 1
gi = 1
gi = 1
gi = 1
g i = 1 gi = 0.875 gi = 0.75 gi = 0.625
gi = 0
gi = 1
gi = 2
gi = 3
gi = 4 gi = 3.25
gi = 0.5 gi = 0.75
gi = 2.5
gi = 1.75
要素密度が0, 0.5, 1の場合のgravity control 関数gi値
gi = 1
g i = 0.5
gi = 1
フィルタリングの導入法
min ª¬C ( α ) = UT KU º¼
α = {a1 , a2 , , aN , b1 , b2 , , bN }
α
subject to :
N
W = ¦ (1 − ai bi ) ≤ W ,
i =1
0 ≤ α i ≤ 1,
i = 1, , 2 N
制約条件として加える
G≥G
mi
ただし， G ( α ) = ¦¦ ª¬ ρi ρ j + (1 − ρi ) (1 − ρ j ) º¼
N
i =1 j =1
N
¦m
i =1
i
最適化問題の解法
T
逐次線形計画法(SLP)
T
T
T
最適性規準法(OC)
T
T
T
制約条件の扱いが易しくロバスト性が高い
対称問題を解いても対称な位相が求まらない
制約条件の数は限られるが収束が速い
対称問題では対称な最適位相が得られる
凸線形化法(CONLIN）
T
T
多数の制約条件を扱え，収束も速い
対称問題では対称な最適位相が得られる。
最適性規準法による解法
ラグラジアンの定義
§ N
·
L ( α ) = C ( α ) − Λ ¨ ¦ (1 − ai bi ) − W ¸ − Λg G − G ( α )
© i =1
¹
歪みエネルギーの２倍
重力制御関数
質量制約条件
（平均ｺﾝﾌﾟﾗｲｱﾝｽ）
制約条件
(
)
a
α = {a1 , a2 , , aN , b1 , b2 , , bN }
：各要素のミクロ構造ユニットセルの穴の大きさ
穴の角度θは要素中心の応力の主軸方向
b
1
θ
最適性規準から得られる更新式
ラグラジアン最小化の条件
N
N
§ ∂C ( α )
∂G ( α ) ·
§ ∂C ( α )
∂G ( α ) ·
§ N
·
δ L(α) = ¦¨
+ bi Λ +
Λg ¸ δ ai + ¦ ¨
+ ai Λ +
Λg ¸ δ bi − ¨ ¦ (1 − aibi ) − W ¸ δΛ − ( G − G ( α ) ) δΛg = 0
∂ai
∂ai
∂bi
∂bi
i =1 ©
i =1 ©
© i =1
¹
¹
¹
ª (k )
∂G ( α ( k ) ) ( k ) º
(k )
« bi Λ +
Λg »
(k )
a
∂
«
» (k )
i
= «−
(k )
» ai
∂C ( α )
«
»
(k )
«¬
»¼
∂ai
β
ai
( k +1)
( k +1)
Λ
k +1)
¦ (1 − a
N
i =1
i
(k )
bi
(k )
)
β
º (k )
» Λ
¼
β
(k )
ª (k )
º
G
∂
α
(
)
(k )
(k )
« ai Λ +
Λg »
(k )
∂bi
«
» (k )
bi
= «−
(k )
»
∂C ( α )
«
»
(k )
«¬
»¼
∂bi
β
bi (
ª1
=«
¬W
Λg (
k +1)
ª G º
» Λg ( k )
=«
(k )
«¬ G ( α ) »¼
更新式の導出法
∂C ( α )
∂ai
+ bi Λ +
ª
∂G ( α ) º
« bi Λ +
Λg »
∂ai
«−
» =1
«
»
∂C ( α ( k ) )
«
»
(k )
«¬
»¼
∂ai
∂G ( α )
∂ai
Λg = 0
β
ª
∂G ( α ) º
« bi Λ +
Λg »
∂ai
«−
» = ai
«
»
ai
∂C ( α ( k ) )
«
»
(k )
«¬
∂ai
¼»
(k )
ª (k )
º
G
α
∂
(
)
(k )
(k )
« bi Λ +
Λg »
(k )
∂ai
» a (k )
= «« −
» i
∂C ( α ( k ) )
«
»
(k )
«¬
∂ai
¼»
β
ai
( k +1)
β：更新幅を制御　するべき乗係数
変数制約条件の考慮
( i = 1,, N ) ,
0 ≤ ai ≤ 1, 0 ≤ bi ≤ 1
ai (
bi (
k +1)
k +1)
Λ( k +1)
Λg
{
}
= min max {0, sai( k ) } , 1
{
ただし，
( )
}
= min max {0, sbi( k ) } , 1
°
= min ®0,
°¯
( k +1)

°
= min ®0,
°¯
ª 1
«
¬ mS
¦ (1 − a
N
i
i =1
(k )
Λ ≤ 0, Λg ≤ 0
bi
(k )
)
β
º ( k ) ½°
» Λ ¾
¼
°¿
sai( k )
(k )
ª (k )
º
α
G
∂
(
)
(k )
(k )
« ai Λ +
Λg »
(k )
b
∂
» b (k )
i
= «« −
» i
∂C ( α ( k ) )
«
»
(k )
«¬
»¼
∂bi
β
β
½
ª G º
(k ) °
«
»
Λ
¾
g
(k )
«¬ G ( α ) »¼
°¿
β
(k )
ª
º
α
∂
G
k ) (k )
k)
(
(
« bi Λ +
Λg »
«
∂ai
» (k )
= «−
» ai
(k )
∂C ( α )
«
»
(k )
«
»
∂ai
¬
¼
sbi( k )
設計変数のムーブリミット
{
ai
( k +1)
bi
( k +1)
}
max (1 − ζ ) ai ( k ) ,0
°
° (k )
= ® sai
°
(k )
°̄min (1 + ζ ) ai ,1
{
}
{
}
max (1 − ζ ) bi ( k ) ,0
°
° (k )
= ® sbi
°
(k )
+
min
1
ζ
b
,1
(
)
i
°̄
{
}
if
sai
(k )
{
}
≤ max (1 − ζ ) ai ,0
(k )
if
{
}
min {(1 + ζ ) a ,1} ≤ s
if
sbi
if
max (1 − ζ ) bi ,0 ≤ sbi
if
if
max (1 − ζ ) ai ,0 ≤ sai
(k )
(k )
i
(k )
{
(k )
{
}
≤ min (1 + ζ ) ai ,1
(k )
(k )
ai
}
≤ max (1 − ζ ) bi ,0
(k )
{
}
min {(1 + ζ ) b ,1} ≤ s
(k )
(k )
i
(k )
{
}
≤ min (1 + ζ ) bi ,1
(k )
bi
ζ ：設計変数の変動幅を制約するムーブリミット
(k )
各ステップ内での更新
( )
ª
º
G
∂
α
(
)
« b ( ) Λ( ) +
Λ ( )»
k
k
( k +1)
ai
s
«
= «−
«
«
¬
i
k +1
k +1
∂ai
∂C ( α ( k ) )
(k )
g
ª (k )
∂G ( α
( k +1)
+
« ai Λ
(k )
∂
b
i
= «−
«
∂C ( α ( k ) )
«
(k )
«¬
∂bi
)Λ
∂ai
» (k )
» ai
»
»
¼
(k )
(k )
sbi( k +1)
β
g
β
( k +1)
º
»
» b (k )
» i
»
»¼
制約条件をアクティブにする。ただし，感度係数は更新しない
最適性規準法のパラメータ
T
T
T
更新式のべき乗係数β
設計変数のムーブリミット ζ
設計変数，ラグランジェ乗数の更新回数
T
T
感度係数を更新する外側ループの繰返し数
感度係数を更新しない内側ループの繰返し数
β = 0.25, ζ = 0.1
外側ループの繰返し数　40回
内側ループの繰返し数 Λ に関しては100回
Λ g に関しては5回
解析例（MBBはり）
L/6
L /150
L
フィルタリングの効果
(a) G = 0, G = 0.73
(c) G = 0.85, G = 0.89, C / C0 = 1.08
(b) G = 0.8, G = 0.87, C / C0 = 1.02
(d) G = 0.90, G = 0.90, C / C0 = 1.18
(a) mS = 0.3, G = 0.85, G = 0.85, C / C0 = 1.05
(b) mS = 0.4, G = 0.85, G = 0.87, C / C0 = 1.08
(c) mS = 0.5, G = 0.85, G = 0.89, C / C0 = 1.08
(d) mS = 0.6, G = 0.8, G = 0.89, C / C0 = 1.03
無償公開ソフト（連続緩和法３兄弟）
T
骨組構造の位相最適化ソフト（Otto）
T
T
２次元連続体の位相最適化ソフト（Isler）
T
T
グランドストラクチャー法
均質化設計法
３次元連続体の位相最適化ソフト（Gaudi）
T
密度法
Ottoの位相最適化手法
－グランドストラクチャー法－
節点配置
グランドストラクチャー
最適位相
Islerの位相最適化手法
－均質化設計法－
マクロ構造
ミクロ構造
x2
a
b
x1
1
y2
y1
θ
Unit cell
Gaudiの位相最適化手法
－密度法ー
Element
ρ = 1.0
ρ =0
ρ = 0.5
Islerの解析例（Single load）
Islerの解析例（Single Load）
Islerの解析例（Multi Load）
Islerの解析例（Multi Load）
Islerの解析例（Distributed load）
Islerの解析例（Distributed load）
Islerの解析例（Distributed load）
Islerの解析例（Single load）
周辺単純支持
中央集中荷重
Islerの解析例（Distributed load）
周辺単純支持
鉛直等分布荷重
まとめ
T
公開ソフトは，下記のホームページからダウン
ロード可能。
http://www.nasl.t.u-tokyo.ac.jp/~dfujii/homepage.htm
T
実務での利用を考えておられる方は下記ソフト
をご購入下さい。
T
T
MSC Nastran　Opti-shape
Altea Opti-struct

Download Report