長岡工業高等専門学校研究紀要 第 53 巻（2017）

論 文

トポロジー最適化シミュレーション－剛性最大化構造に対する考察
および3Dプリンタによる造形－
倉橋 貴彦1・吉原 健太2・丸岡 宏太郎3・井山 徹郎4
1 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
2 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
3 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
4 機械工学科（Department of Mechanical Engineering, National Institute of Technology, Nagaoka College）（現：3東芝機械株式会社）

Numerical Simulation for Topology Optimization-Considerations for stiffness maximum structureand 3D printing of stiffness maximum structureTakahiko KURAHASHI1, Kenta YOSHIHARA2, Kotaro MARUOKA3,and Tetsuro IYAMA4
Abstract
In this study, numerical simulation for topology optimization and 3D printing for stiffness maximum
structure are carried out. The stiffness maximum structure expressed by optimal density distribution isiteratively computed by using the gradient of the Lagrange function with respect to the density. In thiscomputation, the filtering technique for the gradient is introduced, and the update of the density value iscarried out based on optimality criteria method. In numerical experiments, considerations for the stiffness
maximum structure is also performed.
Key Words : topology optimization, FEM, stiffness maximum structure, Lagrange function, filtering
technique, optimality criteria method.

１．序論
近年の数値計算における研究動向の一つとして，2016 年度に開催された国際会議” Asian Congress ofStructural and Multidisciplinary Optimization 20161)(ACSMO2016)”や，” European Community on
Computational Methods in Applied Sciences 20162)(ECCOMAS 2016)”を見ると，トポロジー最適化理論 3),4)に関する研究報告が多数見受けられることがわかる．ACSMO 2016（3 日間）においてはTopology optimization というキーワードを含むセッションは 40 セッションのうち 20 セッション，ECCOMAS2016（5 日間）において Topology
optimization に関連した発表は約 60 件あり，plenarylecture の発表のトピックとしても題材にされている．トポロジー最適化の計算は，構造内における密- 1 -倉橋 貴彦・吉原 健太・丸岡 宏太郎・井山 徹郎

度分布の最適化シミュレーションを行うものであり，外力による仕事量の最小化等を目的に最適な構造形状（外力による仕事量の最小化を目的とした場合は剛性最大化形状）を算定することができる．本
研究では，トポロジー最適化の計算に基づき 5)，得られた形状に対して 3D プリンタにより構造の作成を試みる．また，剛性最大化構造について，同体積の材料を用いた異なる構造モデルとの比較により考察を行う．

ここに，式(4)において{－λ}T[K]={f}T となり， [K]T{－λ} ={f}と書き表すことができる．剛性行列[K]は対称行列であることから，式(1)に示す有限要素方程式より{－λ}={u}が成り立つことになる．そのため，式(4)は式(5)のように書き表すことができる．

∂J * 1 T ∂[K ]= {λ}{u}∂ρ i 2∂ρ i=

２．トポロジー最適化理論に基づく定式化

1 T T ∂[D]{λ} [B][B]{u}2∂ρ i=−
２．１

ラグランジュ関数の各要素密度に関する勾配ベクトルの誘導(5)

1 T T ∂[D]{u} [B][B]{u}2∂ρ i

また，ヤング率が密度（単位は無次元）の関数で与化すると，式(1)に示す有限要素方程式が得られる． えられるものと考え，ヤング率を式(6)のように定義する．Eminは計算不安定性を回避するパラメータ
ここに，要素タイプは四角形の一次要素とする．
であり，本検討ではE0=1MPa，Eminは=1×10-3MPaと[K]は剛性行列，{u}は変位ベクトル，{f}は外力ベク設定する．pはペナルティーパラメータあり，本検トルを示す．
討ではp=3と設定する．
線形弾性体の支配方程式を有限要素法により離散[K ]{u} = { f }(1)

外力による仕事量の最小化を目的とし，評価関数Jを式(2)のように定義する．

J=1{ f }T {u} = 1 {u}T [K ]{u}.
22(2)E i (ρ i ) = (E 0 − E min )ρ ip + E min

式(6)を踏まえると，式(5)は式(7)のように書き表すことができる．結果として，ラグランジュ関数の各要素密度に対する勾配ベクトルは式(7)により計算することができる．

ここに，式(1)の制約条件のもと式(2)の最小化を行うため，式(2)に示すラグランジュ関数を導入する．
{λ}は随伴変数ベクトルを示す．

J* =()

1{ f }T {u}+ {λ}T ([K ]{u}− { f })2

(7)勾配ベクトルに対するフィルタリング処理の導入
式(7)による勾配ベクトルを使用したトポロジー
最適化シミュレーションでは，チェッカーボード状に最適化構造が得られることが報告されており，本検討においては，式(8)に示す勾配ベクトルの修正式により勾配ベクトルに対してフィルタリング処理を導入する．

⎞1 ⎛ T ∂[K ]⎜⎜ {λ}{u}+ { f }T + {λ}T [K ] ∂{u}⎟⎟∂ρ i∂ρ i ⎠2⎝(⎤⎡0 ⎥⎢1 ν1TT
⎢ν 10 ⎥[B]{u}× {u} [B]1 −ν 2 ⎢1 −ν ⎥⎥⎢0 02 ⎦⎣
２．２⎞⎞
∂J * 1 ⎛⎜ T ∂{u}T ⎛ ∂[K ]{u}+ [K ] ∂{u}⎟⎟ ⎟⎟= ⎜{ f }+ {λ} ⎜⎜∂ρ i 2 ⎝∂ρ i∂ρ i ⎠ ⎠⎝ ∂ρ i=1∂J *= − p(E 0 − E min )ρ ip −12∂ρ i(3)

ラグランジュ関数の要素密度ρi（単位は無次元とする．）に対する勾配ベクトルは式(4)のように書くことができる．iは要素番号を示す．

(6))

(4)
- 2 -

トポロジー最適化シミュレーション∂J∂ρ i*( k )∂J∑ w(x , y )ρ ∂ρ=ρ ∑ w(x , y )jj∈M iijjj∈M i

*j

(8)

i

(x − x ) − (y − y )
2j

2i

ji

( )

= ρ i( k ) Ai( k )

(9)

( )

最適性基準法に基づく各要素における密度の更新式について
本研究では，最適性基準法に基づく密度の更新式を反復計算において用いる．まず，領域内の体積Vtotal に対する制約条件を課した式(10)に示す修正ラグランジュ関数を導入する．式(10)における k は反復回数を示す．

)

(10)

修正ラグランジュ関数の第一変分を計算すると式(11)のように表すことができる．

(i = 1,....mx)( )( )⎧max(0, ρ − m) in ρ ( k ) A( k ) η ≤ max(0, ρ − m)iiii⎪⎪
( k +1)(k ) (k ) ηρi= ⎨ min(1, ρi − m) in ρi Ai≤ max(0, ρi + m)
⎪ (k ) (k ) ηinotherwise⎪⎩ ρi Ai
(i = 1,....mx)

２．３(( )η(12)

本検討では，式(12)に基づく更新の際，設計変数に対して上下限値を設け，式(13)による条件付きの更新を行うことにする．

本論文では，Rは1.5と設定する．

J **( k ) = J *( k ) + Λ( k ) V ρ ( k ) − Vtotalη( )j

ここに，重みパラメータ w(xj)は i 番目の要素から半径 R 内に含まれる要素を対象とし，i 番目の要素からの距離に応じた係数値とする（式(9)）．また，式(8)における Mi は，i 番目の要素を中心に半径 R内に含まれる要素の集合を表す．

w(x j , y j ) = R −⎞⎛∂J *( k )⎟⎜⎟⎜ρ∂iρ i( k +1) = ρ i( k ) ⎜(k ) ⎟⎟⎜ − Λ( k ) ∂V ρ⎟
⎜ρ
∂i
⎠
⎝

(13)
ここに，パラメータ m は密度の更新における更新幅を示す．また，式(12)におけるラグランジュの未定乗数 Λ は式(10)における平均密度の制約条件を満たすように計算が行われる．
２．４ 最適密度分布算定のための計算の流れ
以下にトポロジー最適化における最適密度分布算定のための計算アルゴリズムを示す．

1. トポロジー最適化に必要な計算条件の入力を行う．

δJ**(k )( )∂J *( k ) ( k )∂V ρ ( k )=δρi + Λ( k )δρi( k ) = 0∂ρi∂ρi(11)

2. 有限要素解析により線形弾性解析を行う．
3. 収束判定：|J(l+1)-J(l)|/ J(1)<ε を計算し，収束判定
定数 ε より小さい場合は計算を終了する．そう
でなければ，次のステップへ進む．

ここに，ラグランジュ関数の各要素密度に対する勾配ベクトルは，式(8)に示すフィルタリング処理後の勾配ベクトルとする．この第一変分における各変数をもとに，式(12)に示す各要素における密度の更新式を導入する．本研究においては，ダンピングパラメータ η を 0.5 に設定する．

4. ラグランジュ関数の各要素密度に対する勾配ベクトルの計算を行う．

5. 前のステップにおいて求まった勾配ベクトルに対してフィルタリング処理を行う．

6. 最適性基準法により密度分布の更新を行う．ま- 3 -た，この計算におけるラグランジュの未定乗数
Λ は二分法により求める．反復回数を更新し，ステップ 2 に戻る．

倉橋 貴彦・吉原 健太・丸岡 宏太郎・井山 徹郎

３．剛性最大化構造の算定に対する検討

３．２ 剛性最大化構造の３D プリンティング
3.1節において求まった剛性最大化構造の造形を

３．１ トポロジー最適化シミュレーション
本検討では，図-1に示す計算モデル，境界条件に対して，剛性最大化問題に関するトポロジー最適化シミュレーションを行う．計算条件を表-1に示す．
ここでは平面応力状態を仮定し，二次元領域における変形解析を行う．本解析では，2.4節に示す計算ループを30回行い，最終的に得られた密度分布を図-1示す．結果より，チェッカーボード状の密度分布は得られず，部材の有無を適切に表した剛性最大化構造を得ることができた．

目的に，SolidWorksによる作図を行う．作図に際して， 3.1節の結果をSolidWorks上においてトレースすることにより構造の内部の空洞や外形形状の形取り を 行 う （ 図 -3 ） ． ま た ， 比 較 対 象 用 に ，SolidWorksにより作成した剛性最大化形状（図-3）と体積が同等になるように，異なる形のトラス形状を用意する（図-4，図-5）．

P=1N

Design domain

30mm

図-3

剛性最大化構造モデル図

40mm

図-1
表-1

計算モデル図計算条件

要素総数（x方向×y方向）
メッシュサイズΔx，Δy無次元密度の初期値
ポアソン比
収束判定定数

図-2

1200(40×30)
1mm
0.3
0.3
10-5

トポロジー最適化の計算により得られた剛性最大化構造

図-4

図-3 のモデルと同体積の構造モデル 1図-5図-3 のモデルと同体積の構造モデル 2

次に，本節において用意した構造の変形形状を調べるために，SolidWorksによる構造解析を行う．最適密度分布の結果より，ヤング率の分布を式(6)より算定し，全要素におけるヤング率から算出した平均値を構造解析のヤング率とする．またポアソン比は3.1節と同様に設定し，横弾性率Gはヤング率E、ポアソン比 ν の値よりG=E/2(1+ν )により計算したものを与える．材料物性値を表-2に示す．境界条件
は．図-3～図-5の構造に対して，下面の両端を縦方- 4 -トポロジー最適化シミュレーション

向のみ固定し，上面の中心に集中荷重1Nを与える．

以上のことを検証するために，図-3の構造の荷重点近傍において図-2の結果をもとに荷重点近傍の部
表-2 材料物性値
材にヤング率の高いものを用意し，2種類のヤング0.82
ヤング率，MPa率を用いることで構造解析を行った．部材のヤング0.30
ポアソン比，率を変えた構造のモデル図および材料物性値を図-9，0.32
横弾性率，MPa表-3に示す．以上の条件で行った構造解析の結果を図-10に示す．結果より最も変形の小さい構造モデ
各構造モデルに対する変形図を図-6～図-8に示す． ルとなっていることがわかる．
結果より，剛性最大化構造モデルの変形量が一番小さいものの，どの構造モデルにおいても結果の差異部材 2は小さいものとなった．また，結果より剛性最大化構造モデルの最大変位は，荷重点において得られていることがわかる．これは，3.1節において求められた剛性最大化構造において，荷重点近傍における高密度部分，すなわち，高い剛性となる部分を無視し，ヤング率を一律としたことが原因と考えられる．

部材 1

図-9

剛性最大化構造モデル図（部材ごとに剛性を変更した構造モデル）
表-3

図-6

（部材1）ヤング率，MPa（部材2）ヤング率，MPaポアソン比，（部材1）横弾性率，MPa（部材2）横弾性率，MPa

剛性最大化構造モデルに対する変形図

図-7

材料物性値0.8239.630.300.3215.24

図-10 図-12 の構造モデルに対する変形図

構造モデル 1 に対する変形図

このような形でモデル化した構造は，現在においては3Dプリンタによる造形も可能であり，図-11，図-12のように作成することができる．また，3Dプリンタによる造形は金属によるものも可能であり，本解析モデルのように密度分布や，材料を変えることで剛性を調節しながら最適構造を作製することで，同材料コストにおいても，変形量を最小とする構造の作製が可能となると考えている．
図-8構造モデル 2 に対する変形図- 5 -倉橋 貴彦・吉原 健太・丸岡 宏太郎・井山 徹郎

謝辞：本研究を進めるにあたり，長岡技術科学大学学長戦略的経費（⑤材料科学領域：未来を牽引する高信頼性材料分野）の援助を受けた．また，本論文の解析結果は，九州大学情報基盤研究開発センターの高性能アプリケーションサーバSR16000を使用し計算を行ったものである．ここに謝意を表す．
図-11 剛性最大化構造モデルに対する樹脂構造

参考文献
1)

ACSMO2016,
available from < http://www.acsmo2016.com/>

2)

ECCOMAS Congress 2016,
available from < https://www.eccomas2016.org/>3)

西脇眞二，Frecker, M. I.，閔相載，菊池昇，”柔軟性を考慮した構造の最適化（第 1 報

図-12 構造モデル 1 に対する樹脂構造定式化とコンプライアントメカニズムへの応用）”，日本機械学会論文集(C 編), Vol. 63, No. 612, pp. 2657-2664, 1997.
4)

４．結論

山田崇恭，山崎慎太郎，西脇眞二，泉井一浩，吉村允孝，”レベルセット法に基づくコンプライントメカ
ニ ズ ム の 構 造 最 適 化”， 日 本 計 算 工 学 会 論 文 集 ，No.20080001, 2008.

本研究では，トポロジー最適化の計算に基づき，弾性係数の大きさの違いによる剛性最大化形状の比5)較を行った．定式化においては，外力仕事による評

西脇眞二，泉井一浩，菊池昇，”トポロジー最適化”，丸善出版，2013．

価関数に対して，線形弾性体の変形に対する支配方程式を考慮したラグランジュ関数を導入し，ラグランジュ関数の各要素密度に対する勾配ベクトルの誘導を行った．本論文において，紹介した方法により適切に密度分布を最適にした剛性最大化形状を得ることができ，3D プリンタによる造形モデルに対する構造解析の検討を行った．結果として，密度分布により得られた形状も重要であるが，密度分布の高い箇所（剛性の高い箇所）が最大変位に最も関連していることを確認することができた．
本論文において紹介を行ったトポロジー最適化は，現実的に形成が不可能と考えられる構造形態の結果を算出することもあることから，従来，解析結果が机上の空論となることがあったが，この解析方法は 3D プリンタ技術と最も親和性が高く，剛性最大化等の最適形状も現状においては形成可能となってきている．将来的に，
3D プリンタの解像度の向上し，使用する材料を自由に変えることが可能となった場合には，このトポロジー最適化技術は最も威力を発揮するものと考えられる．

- 6 -(2017. 7. 21 受付)