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

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

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

最適性規準法の続き

　図1の*Bで示した「Dが1に近づくようにρ_iを更新」の方法を説明します。λがマイナス値であることを意識して式6を眺めてみましょう。有限要素法解析の結果、ある要素の応力が大きな値であったとします。式6の分母が大きくなり、Dを1に近づけるにはρ_iを大きくしなければなりません。反対に要素の応力が小さな値のときはρ_iを小さくしなければなりません。参考文献［1］は密度法ではなく「均質化法」なのですが、次の世代のρ_iを求める式があります。しかし、ρ_iの変化量が大き過ぎると反復計算が発散するのでρ_iの変化量にリミッターをかけています。ほとんどの場合、リミッターに引っ掛かりますので、リミッターの数値をρ_iの変化量とします。次式でρ_iを更新します。

式10

式11

　もう1つリミッターをかけます。ρ_iは1以上にはなれず、かつマイナス値にもなれません。次式となります。

式12

式13

　この辺りは「Excel」のマクロプログラムに仕込んでおらず、シートのセルに式として設定しています。任意に数値やアルゴリズムを変えられると思います。

　図1の*Cで示した「λを更新」の方法を説明します。式7が成立するようなλを求めます。式7を眺めると、ρ_iが大きくなれば、gも大きくなります。図2左図のようになります。

　次に式5が1と等しいとして、その式を変形します（式14）。

式14

　λとρ_iの関係は、nが2か3で、マイナスの反比例です。図2中央図のようになります。そうすると、gとλの関係は図2右図になるのでしょうか。

図2　g、ρi、λの関係［クリックで拡大］

　式7、つまりg＝0が目標です。図2右図から、gがマイナスのときはλを増やし、gがプラスのときはλを減らせばよいということになります。次式でλを更新します。

式15

式16

　前述の操作でgの符号が反転すればしめたものです。以降は図3のように、1次式で新しいλを求めて、g／V_oが許容値（1％で十分です）以下になれば、λが求まったことになります。あと1つ、ρ_iは0.05［-］刻みで変化するので、λを変えてもgが変化しなくなることがあります。よって、gが変化しなくなったらλが求まったことにして、図1の内側の反復計算を終了します。

図3　λの求め方［クリックで拡大］

　図1の*Dで示した「各ρ_iが0ないし1になったか」の判定と計算終了の判断は人間が行います。図4に最適化した例を示します。グレーの領域があります。全ての要素のρ_iが完全に0か1になるまで計算を続けると、大切な小骨がなくなってしまいます。最後の判断は人間で行いましょう。

図4　最適化した例［クリックで拡大］