　上式が満足するような領域Ωsを求めればよいのですが、領域Ωsがスプライン曲線などでパラメトリックに定義されていればできそうですが、今回の領域Ωsはそうではないので求められそうにあません。トポロジー最適化でやったように、領域の形状を少しずつ変化させていって、ラグランジュ関数Lが最小になる形を求めていきます。繰り返し作業になります。

　ラグランジュ関数の導関数が、形状勾配関数Gを係数関数とする速度場の一次形式で与えられたので、参考文献［1］［2］［3］で提案されている「力法（traction method）」と名付けられた方法が使えます。力法は、k回目の領域変動を表す速度場V^（k）を次の支配方程式に基づいて解析する方法です。

式28

　この支配方程式を満足する領域変動Vは、凸性が保証されている問題において、ラグランジュ関数Lを必ず減少させるとのことです。ということは今回の問題は、領域の形Ωsを直接求めるのではなく、目的関数l（v）が小さくなるような形状変更の方向Vを求める問題に変わりました。毎回微小移動量Vを求めて、節点座標を移動させていくことになります。

　式28は、仮想変位の原理の式と同じ形です。つまり、弾性変形問題を解けばその変位量がVとなります。式21の形は体積力fがないとした式7とそっくりです。つまり、弾性変形問題の外力はGとなります。ある条件を満たせばGは次式となるとのことですが、一番大切なこの式の導出過程を参考文献［1］［2］［3］から見つけることができず、文献を読んで試みたのですが筆者の力量不足で導出できませんでした……。よって、いきなり掲載することにします。うまく説明できなくて少し残念です。

式29

　この式29の第1項はひずみと応力の積です。トポロジー最適化でも使いましたね。第2項はラグランジュ乗数なのですが、体積制約式3を満足させるための調整定数として扱います。nは領域の輪郭上の外向き法線ベクトルです。

　式29は平均コンプライアンス最小化から出発した剛性最大化の解を導く式です。しかし、次回以降で手を動かすと気が付くのですが、この式は「応力の大きい所を少し膨らせればよい。しかし、応力値に従って膨らませると領域表面が波打つので、弾性変形解を次世代の領域として使って、波打ちを抑える作戦」と読み取ることができます。剛性最大化と応力均一化は“お隣さん”なのだと思っています。

　以上が力法の説明ですが、歯切れが悪くて申し訳ありません。しかし、結果としてこの方法が有効だったので使うことにします。

　SF漫画・アニメ『銀河鉄道999』の中で、メーテルが999号のことを「遠い外宇宙の滅亡した科学惑星の遺跡や異邦人から手に入れた、まだ人類には理解できないけれど、とにかく使うことだけはできる資料を基に作られた機関車……」と説明していますが（参考文献［4］）、それを聞いた気分でExcelシートを作りました。