トポロジー最適化 〜密度法の定式化〜：フリーFEMソフトとExcelマクロで形状最適化（5）（4/5 ページ）

原理原則を押さえていれば、高額なソフトウェアを用意せずとも「パラメトリック最適化」「トポロジー最適化」「領域最適化」といった“形状最適化”手法を試すことができる！　本連載ではフリーのFEM（有限要素法）ソフトウェア「LISA」と「Excel」のマクロプログラムを用いた形状最適化にチャレンジする。連載第5回では、トポロジー最適化における密度法の定式化に取り組む。

[高橋良一／ＲＴデザインラボ 代表，MONOist]

密度ρ_iの導入と有限要素法ソフト「LISA」との相性

　密度法ではヤング率を次式で定義します。添字iはi番目の要素という意味です。密度が小さいとヤング率も小さくなり、密度が1のときヤング率は元の材料のヤング率E_oに等しくなります。これをラグランジュの未定乗数法による式に導入しましょう。

式33

　上式のヤング率を有限要素法ソフトに入力したときの、ひずみと応力がどうなるかを考察します。ここでは話を簡単にするため、図7に示す平面応力三角形要素を使います。

図7　平面応力三角形要素［クリックで拡大］

　荷重ベクトル｛f_i｝、変位ベクトル｛u_i｝とそれらを関連付ける剛性マトリクス［k_i］を以下に示します。ここでのiは図7のi節点ではなく要素番号です。

式34

式35

　連載第4回の剛性マトリクスの式を以下に示します。

式36

　式33から、素材のヤング率はE_o、それぞれの要素の最適化後のヤング率はE_iなので、ヤング率E_oを使った剛性マトリクス［k_i（E_o）］は式37、ヤング率E_iを使った剛性マトリクス［k_i（E_i）］は式38となります。

式37

式38

　ヤング率としてE_oを使った場合と、E_iを使った場合のひずみを比較してみましょう。まず、E_oを使った場合です。式35を左から［k_i（E_o）］⁻¹を掛けます。肩の「−1」は逆マトリクスを意味します。

式39

　変位ベクトルからひずみを求めます。連載第4回で出てきた、変位−ひずみマトリクス［B］を使います。

式40

　今度は、E_iを使った場合です。式35を左から［k_i（E）］⁻¹を掛けます。

式41

　そして、式38を代入します。

式42

　変位−ひずみマトリクス［B］は、要素の形状で決まるのでヤング率が変わっても変化はありません。E_iを使った場合のひずみは次式となります。

式43

　ヤング率E_i＝ρ_iⁿE_oを使った場合、ひずみはヤング率E_oを使った場合の「1／ρ_iⁿ」倍になります。ρ_iは0〜1の実数なのでひずみは大きくなります。そりゃそうですね。ヤング率が小さくなったら変形は大きくなって、ひずみも大きくなります。

　応力はどうでしょうか。調べてみましょう。連載第4回で出てきた応力−ひずみマトリクス［D］を使います。ヤング率がE_oの場合の応力は式44で表されます。

式44

　ヤング率がE_iの場合は以下となります。

式45

　式43を代入しましょう。

式46

　式44と式46は等しくなりました。応力はヤング率が変わっても変化はありません。静定はりの曲げモーメントは形状と荷重だけで決まるので、ヤング率が変わっても曲げ応力は変わらないのと同じですね。

　ヤング率を変えても応力が変わらない性質を利用しましょう。有限要素法ソフトLISAが出力する応力値を使って密度法で最適化問題を解くとき、式33で示したようにヤング率をいろいろと変化させるのですが、ヤング率を変えても目的関数の元となる応力値は変わりません。よって、LISAが出力する応力値をそのまま使うことができます。