密度法によるトポロジー最適化の準備を行う：フリーFEMソフトとExcelマクロで形状最適化（6）（1/6 ページ）

原理原則を押さえていれば、高額なソフトウェアを用意せずとも「パラメトリック最適化」「トポロジー最適化」「領域最適化」といった“形状最適化”手法を試すことができる！　本連載ではフリーのFEM（有限要素法）ソフトウェア「LISA」と「Excel」のマクロプログラムを用いた形状最適化にチャレンジする。連載第6回では「最適性規準法」を用いて、剛性を最大化＝平均コンプライアンスを最小化する設計パラメータを求める。

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

　これまで式の導出が続きましたが、いよいよ最適化アルゴリズムの本丸です。「最適性規準法（Optimality Criteria法）」（参考文献［1］［2］）」を使って、剛性を最大化＝（イコール）平均コンプライアンスを最小化する設計パラメータρ_iを求めましょう。

⇒ 連載バックナンバーはこちら

参考文献：

• ［1］Katsuyuki Suzuki、Noboru Kikuchi｜A homogenization method for shape and topology optimization、Computer Methods in Applied Mechanics and Engineering 93、P291-318（1991）
• ［2］弓削、菊池｜塑性変形を受ける二次元構造物の最適設計（第1報）／日本機械学会論文集（A編） 62巻 600号 P1910（1996）

最適性規準法

　ラグランジュ関数Lは以下の式でした（式1）。

式1

　これをρ_iで微分したものが0です。

式2

　では、最適性規準法を用いて式2を満足させるρ_iを求めます。式2を式3のように変形します。

式3

　両辺を∂f／∂ρ_iで割ります（式4）。

式4

　そして、連載第5回の結果を代入します（式5）。

式5

　スカスカ状態の見掛けの体積はV_oで、真の体積はaV_oでした。設計パラメータである密度ρ_iの平均値は、ρ_iの定義からaと等しくなります。これから反復計算をしますが、ρ_iの初期値はaとなります。ρ_iは最終的な解ではないので式5は成立しません。次式のように変数Dを導入します（式6）。

式6

　式1をλで微分したものも0でした。次式となります（式7）。

式7

　計算していくと分かることなのですが、λも未知数で反復計算して求めます。反復計算の最中はλの変化に従って、ρ_iも変化します。よって、式7においてはρ_iをλの関数として表現しました。

　計算手順を図1に示します。2重ループになっています。

図1　計算アルゴリズム［クリックで拡大］

　図1の*Aで示した「λの初期値を設定」について説明します。式6において分子にヤング率E_oがあり、鋼の場合E_o＝200×10⁹［Pa］です。分母は応力の2乗です。式6のDを1に近づけるためのλは、10⁹のオーダーか、1桁台か、10⁻⁹のオーダーか見当がつきません。λを反復計算で求めるのですが、やみくもな値を初期値にすると答えが得られません。そこで、以下の発想でλの初期値を決めました。

　式3を変形し、連載第5回の結果を代入します（式8）。

式8

　i＝1、2、3……neと、ne通りのλがあります。これらの平均値をλの初期値としましょう。式9です。λがマイナス値になることは「密度法」ではお約束です。

式9

　式9の出番は図1の大きなループの1回目だけです。2回目以降は前回求まったλを初期値にします。