構造物のトポロジー最適化 について 東京大学藤井大地

構造物のトポロジー最適化

について

東京大学 藤井大地

講演内容

T 力学教育，設計教育に用いることを目的と して開発した骨組，連続体の位相最適化 ツールについて T 位相最適化における問題点とその解決法 について T 位相最適化に関するその他の研究につい て

本日の発表資料は下記に掲載

http://www.nasl.t.u-tokyo.ac.jp/~dfujii/homepage.htm 書き物の中にあります。

力学系科目の衰退

T 構造力学系の科目の削減 T 構造力学→構造解析学 T 骨組，トラスの解法→有限要素法 T なぜ構造力学が衰退したか？ T 解析ソフトの発達で手計算で問題を解く必要 がなくなった。 T 学生にとって魅力ある講義が行われていない。

建築学科の現状

T 建築士の試験があるので，他の学科に比 較すると構造力学教育に力を入れている。 T しかし，全体から見ると，構造系は，計画 系，環境系に押されて縮小されてきている。 T 学生の人気は計画系，環境系に傾き，構造系 への志望学生が年々減少している。 T 力学の講義は一般に難しいと思われて敬遠さ れる。

魅力ある力学教育の必要性

T コンピュータによる解析は万能ではない！ T 結果が正しいかどうかの判断は人間が下す 必要がある。 T そのためには，力学的センスを身につけること が重要。 学生に自然に興味を抱かせるような 魅力ある力学教育が必要

力学教育のツール

T 簡単な実験 T 理解しやすい。 T 実験機材にコストがかかる。 T 教える側の準備がかなり大変。 T コンピュータソフト T 教える側の負担は少ない。 T 学生も興味をもつ。 T 力学を教えることのできる簡単なソフトは少ない。 T 意外にコストがかかる。

本研究の目的

T 力学とデザインを結びつけるようなソフト T 学生の興味を引きつけるようなソフト T 大学の講義・演習等で使えるソフト T やってみたいと思い立ったらすぐに使える ソフト そういうソフトの開発！

本研究の開発ソフト

T 骨組構造の位相最適化ソフト（Otto）

T ２次元連続体の位相最適化ソフト（Isler） T ３次元連続体の位相最適化ソフト（Gaudi）

Ottoの位相最適化手法

－グランドストラクチャー法－

グランドストラクチャー 最適位相

Islerの位相最適化手法

－均質化設計法－

a b θ 1 Unit cell 2 x 1 x 2 y 1 y ミクロ構造 マクロ構造

Gaudiの位相最適化手法

－密度法ー

0 ρ = 0.5 ρ = 1.0 ρ = Element

位相を求める問題から入る

？

以下の設計領域，荷重，境界条件が与えられたとき， 重量○○以下で，より変形の小さい構造物を作れ

最初は思考錯誤

解析ソフトで設計の良否を検討

変位を比較して互いに競争させる

応力の概念も武器に使わせる

軸力図 曲げモーメント図 Ｏｔｔｏによる結果表示

さらに最適化ツールを使う

本当にそうなのか確かめる

応力の観点からも考察する

連続体への拡張でデザインへ

デザインのセンスも身に付く？

キャドを使ってデザイン

３次元だともっと面白いデザイン？

建築デザインにも応用可能

Design domain c b 50cm 50 cm w 0.2w fixed fixed 50cm 門形ラーメンの最適なカタチは？

こんな形なんだ！

こんなラーメンもありえる

Gaudiのような建築家が育つ？

グエル公園内の立体道路（ガウディー）

位相最適化とは？

位相最適化の手法

T 連続緩和法（整数条件の緩和） T 均質化設計法（ミクロ的な材料のON/OFF） T 密度法（材料定数が密度のべき乗に比例） T グランドストラクチャー法 T 離散最適化問題として解く方法 T 遺伝的アルゴリズムを利用した方法 T セル・オートマトンを利用した方法

連続緩和法の利点と欠点

T 過去の優れた最適化問題の解法を利用す ることができる。 T 局所最適解の可能性はあるが必ず最適解 が求まる。 T 大規模問題も容易に解くことができる。 T 数学的知識が必要で，学生がとっつきにく い。 T 多峰性の問題への適用が難しい。

離散最適化法の利点と欠点

T ヒューリスティックな方法は，ルールが単純で，数 学がわからなくてもよい。 T 自然にあるものが最適形態をもっているとすれ ば，生命科学にヒントを得たアルゴリズムを適用 することは自然である。 T 感度計算の必要がなく，多峰性の問題に適用可 能である。 T 最適解かどうかの理論的保証がない。 T 大規模問題への適用が難しい。

なぜ連続緩和法か？

T ヒューリスティックな方向に流れすぎている のではないか？ T ヒューリスティックな方法は結局膨大な時 間がかかる？ T ヒューリスティックな方法は，理論的な保証 がないため論文になりにくい？

連続緩和法の問題点

T 連続緩和したために，材料のON/OFFだけ でなく，ONとOFFの中間的な材料が最適 解に現れる。 T 低自由度の要素を用いた場合，チェッカー ボード状の密度分布が求まる。 フィルタリング法の導入

均質化設計法

(HDM)

a b θ 1 Unit cell 2 x 1 x 2 y 1 y ミクロ構造 マクロ構造 質量制約条件下で歪みエネルギーの最小化

均質化設計法の解

グレースケール

密度法

0 ρ = 0.5 ρ = 1.0 ρ = Element 質量制約条件下で歪みエネルギーの最小化

( )

2 e ρ = ρ e K K

密度法の解

重力制御関数

( ) 1 1 1 1 1 1 , 1 , 0 1 i m N N N N i i i j i j i i i i j i i i i i i G g m m a b G ρ ρ ρ ρ ρ ρ ρ = = = = = = = ⋅ + ⋅ = − = − < ≤ ¦ ¦ ¦¦ ¦ 要素密度が0, 0.5, 1の場合のgravity control 関数g_i値 m_i = 4 i 4 i g = 1 i g = 0 i g = 3 i g = 1 i g = 1 i g = 2 i g = 1 i g = 2 i g = 1 i g = 1 i g = 3 i g = 0 i g = 1 i g = 4 i g = 0.25 i g = 0.875 i g = 3.25 i g = 0.75 i g = 0.625 i g = 1.75 i g = 0.5 i g = 0.75 i g = 2.5 i g = 1 i g = 0.5 i g = 1 i g =

フィルタリングの導入法（

HDM）

( ) { } ( ) 1 2 1 2 1 min , , , , , , , subject to : 1 , 0 1, 1, , 2 T N N N i i i i C a a a b b b W a b W α i N = ª = º ¬ ¼ =

¦

− ≤ ≤ ≤ = U KU    αααα αααα α =α =α =α = ( ) ( )

(

)

1 1 1 1 1 i m N N i j i j i i j i G ρ ρ ρ ρ m = = = ª º = ¦¦_¬ + − − _¼ ¦ αααα ただし， G ≥ G 制約条件として加える 制約条件として加える 制約条件として加える 制約条件として加える

最適化問題の解法

T 最適性規準法(OC) T 設計変数の数に依らず収束が速い。 T 目的関数によっては収束解が得られない場合がある。 T 制約条件の扱いが煩雑。 T 逐次線形計画法(SLP) T 局所解を見つけやすい。対称解析でも対称な解が得られない。 T 制約条件の扱いが簡単で，収束に関するロバスト性が高い。 T 汎用サブルーチンを利用できる。 T CONLIN T 設計変数と設計変数の逆数でテーラー展開した問題を双対法を用いて 解く方法であり，収束が速く，汎用性がある。 T 現在導入を検討中

最適性規準法による解法

( ) ( ) ( )

(

( )

)

1 1 N i i g i L C Λ a b W Λ G G = § · = − _¨ − − _¸ − − ©

¦

¹ α α α αα αα αα α α α ラグラジアンの定義 質量制約条件 重力制御関数_制約条件 歪みエネルギーの２倍 （平均ｺﾝﾌﾟﾗｲｱﾝｽ） {a a1, 2,,aN , ,b b1 2,,bN} α = α = α = α = a b θ 1 ：各要素のミクロ構造ユニットセルの穴の大きさ 　穴の角度θは要素中心の応力の主軸方向

( ) ( ) ( ) ( ) ( ) ( ) ( ( )) 1 1 1 1 0 N N N i g i i g i i i g i i i i i i i C G C G L b a a b a b W G G a a b b δ Λ Λ δ Λ Λ δ δΛ δΛ = = = §∂ ∂ · §∂ ∂ · § · = _¨ + + _¸ + _¨ + + _¸ −_¨ − − _¸ − − = ∂ ∂ ∂ ∂ _{© ¹} © ¹ © ¹ ¦ αααα αααα ¦ αααα αααα ¦ α α α α α α α α ラグラジアン最小化の条件 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 1) ( ) ( ) ( ) k k k k i k g k _i k i _k i k i G b a a a C a β Λ Λ + ª ∂ º « + » ∂ « » = −_{« »} ∂ « » « ∂ » ¬ ¼ αααα αααα ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) k k k k i k g k _i k i _k i k i G a b b b C b β Λ Λ + ª ∂ º « + » ∂ « » = −_{« »} ∂ « » « _∂ » ¬ ¼ αααα αααα ( )1 ( ( ) ( )) ( ) 1 1 1 N k k k k i i i a b W β Λ + Λ = ª º = _« − _» ¬ ¦ ¼ ( ) ( ) ( ) 1 ( ) k k g _k g G G β Λ + ₌ ª_« º_» Λ « » ¬ αααα ¼

最適性規準から得られる更新式

更新式の導出法

( ) ( ) 0 i g i i C G b a Λ a Λ ∂ ∂ + + = ∂ ∂ α α αα αα α α ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 1) ( ) ( ) ( ) k k k k i k g k _i k i _k i k i G b a a a C a β Λ Λ + ª ∂ º « + » ∂ « » = −_{« »} ∂ « » « ∂ » ¬ ¼ αααα αααα ( ) ( )( ) ( ) 1 i g i k k i G b a C a Λ Λ ª _∂ º « + » ∂ «₋ » ₌ « _∂ » « » « ∂ » ¬ ¼ αααα αααα ( ) ( )( ) ( ) i g i i k i k i G b a a a C a β Λ Λ ª _∂ º « + » ∂ «₋ » ₌ « _∂ » « » « ∂ » ¬ ¼ αααα αααα β：更新幅を制御　 　　するべき乗係数

変数制約条件の考慮

( ) 0 ≤ ≤a_i 1, 0 ≤ ≤b_i 1 i =1,,N , Λ ≤ 0, Λ_g ≤ 0 ( )1 { { ( )} } min max 0, , 1 k k i ai a + = s ( )1 { { ( )} } min max 0, , 1 k k i bi b + = s ( )1 ( ( ) ( )) ( ) 1 1 min 0, 1 N k k k k i i i S a b m β Λ + Λ = ­ _{ª º} ½ ° ° = _{® «} − _{» ¾} ¬ ¼ ° ° ¯ ¦ ¿ ( ) ( ) ( ) 1 ( ) min 0, k k g _k g G G β Λ + ­° ª º Λ ½° « » = _{® ¾} « » ° _{¬ ¼} ° ¯ αααα ¿ ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) k k k k i g k k i ai _k i k i G b a s a C a β Λ Λ ª _∂ º « ₊ » « ∂ » = −_{« »} ∂ « » « _∂ » ¬ ¼ αααα αααα ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) k k k k i k g k k i bi _k i k i G a b s b C b β Λ Λ ª _∂ º « + » ∂ « » = −_{« »} ∂ « » « ∂ » ¬ ¼ αααα αααα ただし，

設計変数のムーブリミット

( ) { } ( ) { } ( ) ( 1) ( ) ( ) max 1 ,0 min 1 ,1 k i k k i ai k i a a s a ζ ζ + ­ ₋ ° ° = _® ° ₊ °¯ ( ) { } ( ) { } {( ) } ( ) { } ( ) ( ) ( ) ( ) ( ) ( ) ( ) if max 1 ,0 if max 1 ,0 min 1 ,1 if min 1 ,1 k k ai i k k k i ai i k k i ai s a a s a a s ζ ζ ζ ζ ≤ − − ≤ ≤ + + ≤ ( ) { } ( ) { } ( ) ( 1) ( ) ( ) max 1 ,0 min 1 ,1 k i k k i bi k i b b s b ζ ζ + ­ ₋ ° ° = _® ° ₊ °¯ ( ) { } ( ) { } {( ) } ( ) { } ( ) ( ) ( ) ( ) ( ) ( ) ( ) if max 1 ,0 if max 1 ,0 min 1 ,1 if min 1 ,1 k k bi i k k k i bi i k k i bi s b b s b b s ζ ζ ζ ζ ≤ − − ≤ ≤ + + ≤ ζ ：設計変数の変動幅を制約するムーブリミット

各ステップ内での更新

( ) ( )

( )

( ) ( )

( )

1 1 ( ) ( ) ( 1) ( ) ( ) k k k k i k g k k i ai _k i k i G b a s a C a β Λ + Λ + + ª _∂ º « ₊ » « ∂ » = −_{« »} ∂ « » « _∂ » ¬ ¼ αααα αααα

( )

( )

( ) ( ) _{( 1)} ( 1) ( ) ( ) ( 1) ( ) ( ) k k _k k i k g k k _i bi _k i k i G a b s b C b β Λ ₊ Λ + + ª _∂ º + « » ∂ « » = −_{« »} ∂ « » « ∂ » ¬ ¼ αααα αααα 制約条件をアクティブにする。ただし，感度係数は更新しない

最適性規準法のパラメータ

T 更新式のべき乗係数β T 設計変数のムーブリミット T 設計変数，ラグランジェ乗数の更新回数 T 感度係数を更新する外側ループの繰返し数 T 感度係数を更新しない内側ループの繰返し数

ζ

0.25, 0.1 β = ζ = 外側ループの繰返し数　40回 内側ループの繰返し数 Λ に関しては100回 g Λ に関しては5回

解析例（

MBBはり）

L / 6 L /150 L

フィルタリングの効果（

HDM）

(a) G =0, G =0.73 0 (b) G =0.8, G =0.87, C C/ =1.02 0 (c) G =0.85, G =0.89, C C/ =1.08 0 (d) G=0.90, G=0.90, C C/ =1.18

0 (a) m_S =0.3, G =0.85, G =0.85, C C/ =1.05 0 (b) m_S =0.4, G =0.85, G =0.87, C C/ =1.08 0 (c) m_S =0.5, G =0.85, G =0.89, C C/ =1.08 0 (d) m_S =0.6, G =0.8, G =0.89, C C/ =1.03

フィルタリングの効果（密度法）

(a) G =0, G =0.61, m =0.50 (b) G =0.75, G =0.84, m =0.50, C C/ ₀ =1.04

0

0 (a) m_S =0.3, m=0.30, G=0.9, G=0.89, C C/ =1.49 0 (b) m_S =0.4, m=0.40, G=0.80, G =0.86, C C/ =1.10 0 (c) m_S =0.5, m=0.50, G =0.75, G=0.84, C C/ =1.04 0 (d) m_S =0.6, m=0.60, G=0.75, G =0.84, C C/ =1.02

弾性変形機構の位相最適化

ツール

弾性変形機構とは？

T 弾性変形機構は，従来の剛性を最大化す る構造とは逆に，変形を生み出す構造形 態であり，柔軟性が要求される構造を設計 する上で重要となる。 T 現在，マイクロ構造のメカニズムの設計等 にニーズがある。

絶対変位の最大化と

相対変位の最大化

Extend design domain Ω

A B (a) 拡張された設計領域   (b) 剛性最大化問題の位相 (c) 弾性変形機構の位相 （絶対変位の最大化） (d) 弾性変形機構の位相  （相対変位の最大化）

位相最適化のための

境界条件と荷重条件

Extend design domainΩ

A

B

F

Extend design domainΩ

A

B

F

Extend design domainΩ

A

B

1

Extend design domainΩ

A

B 1

CASE 1 CASE 2

目的関数

( ) ( ) ( ) ( ) ( ) ( ) ( )

( )

{ } 1 1 1 1 1 2( ) 2 1 1 2 1 2( ) 3 3 3 3 3 2( ) 4 4 4 4 4 2( ) ( ) ( ) 1 ( ) ( ) 1 T G A T G B T G A T _G B L F u d L u d L F u d L u d Ω Ω Ω Ω = ⋅ = Ω = ⋅ = Ω = ⋅ = Ω = ⋅ = Ω ³ ³ ³ ³ u u D u u u D u u u D u u u D u ε ε ε ε ε ε ε ε ε ε εε εε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε

( ) ( )

( ) ( )

2 1 4 4 1 1 3 3 Maxmize L L L L § ₋ · ¨ ¸ ¨ ₊ ¸ © ¹ u u u u （A点から入力する仕事量とB点から出力される仕事量の比）

最適化問題の定式化

( )

₍

(

_{( )}( )

₎

)

₍

(

_{( )}( )

₎

)

{ } ( ) ( ) 2 1 4 4 1 1 3 3 1 1 1 1 4 4 4 4 Maxmize , , , , , , , , subject to 1 10 1 10 , 10 1 10 ( 1, , ) N N N N i i i i i L L C L L a a b b W a b W a b i N θ θ = − − − − § ₋ · ¨ ¸ = ¨ ₊ ¸ © ¹ = = − ≤ ≤ ≤ − ≤ ≤ − =

¦

X u X u X X u X u X X X    

均質化設計法の適用

Extend design domainΩ A B 1 x 2 x a b θ 1 U n it c ell 1 y 2 y

解析例１

Extend design domain Ω

A B 60 30 30 1 F =

解析例２

Extend design domainΩ

A B 60 30 30 1 F = CASE 1 CASE 2 CASE 3 Case1 Case2 Case3

解析例３

Extend design domainΩ

A B 60 30 1 F= CASE 1 CASE 2 CASE 3 CASE 4 Case1 Case2 Case3 Case4

材料の内部構造の位相最適化

ツール

研究の背景

T 光造形法等により製造技術の革命的進歩 が起こっている。 T 将来このような製造技術を用いて様々な 新しい複合材料を開発できる可能性があ る。

複合材板のミクロ構造の設計

Design domain (Unit cell) 1 1 y₂ y₁ Periodic microstructures Phase A material Phase B material Phase C material A macrostructure x₂ x₁

複合材板の材料最適設計

1m 1kN/m 1m t=0.01m t=0.01m 1m E (GPa) ν

Gray Material Cast Epoxy resin 3.0 0.25

Black Material E-Glass Fiber 72.4 0.15

応用の可能性

T 環境に優しい材料の開発

T 剛性を落とさずにリサイクル不可能材を吸収

まとめ

T 力学，設計教育に利用することを目的とした位相 最適化ツールの開発について紹介した。 T 連続体の位相最適化における問題点とフィルタ リング法を用いる解決法について示した。 T 位相最適化に関する他の研究として弾性変形機 構，および材料の内部構造の位相を求める方法 を示した。

Otto, Isler, GaudiについてはWebで公開中