デジタルツインを実現するCAEの真価

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

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

» 2022年03月24日 11時00分 公開

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

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

最適性規準法

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

式1 式1

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

式2 式2

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

式3 式3

 両辺を∂f/∂ρiで割ります(式4)。

式4 式4

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

式5 式5

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

式6 式6

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

式7 式7

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

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

計算アルゴリズム 図1 計算アルゴリズム[クリックで拡大]

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

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

式8 式8

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

式9 式9

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

       1|2|3|4|5|6 次のページへ

Copyright © ITmedia, Inc. All Rights Reserved.