JP2004348691A

JP2004348691A - 連続的材料分布を用いた位相最適化手法

Info

Publication number: JP2004348691A
Application number: JP2003175543A
Authority: JP
Inventors: Keizo Ishii; 惠三石井; Kazumi Matsui; 和己松井; Kenjiro Terada; 賢二郎寺田; Shinji Nishiwaki; 眞二西脇
Original assignee: QUINT CORP
Current assignee: QUINT CORP
Priority date: 2003-05-19
Filing date: 2003-05-19
Publication date: 2004-12-09
Anticipated expiration: 2023-05-19
Also published as: JP3924701B2

Abstract

【課題】構造物の鮮明な最適位相形態を生成するシステムを提供する。
【解決手段】対象構造の有限要素モデルデータ（１．１）および位相最適化制御データ（１．２）を入力する。連続的材料分布を実現するために各有限要素の頂点である節点において設計変数として定義した材料密度を用いて応力−歪マトリックスデータを作成する（２）。前記算出された応力−歪マトリックスを用いて有限要素解析を行う（３）。目的関数、制約関数およびそれぞれの感度（４．１）、（４．２）を計算する。前記有限要素解析および感度解析により算出された目的関数と予め定められた閾値との関係を評価し収束判定する（５）。収束判定の結果に応じて設計変数である節点における材料密度を更新する（６）。前記材料密度を更新した後、応力−歪マトリックスの算出、有限要素解析、感度解析、収束判定、および設計変数の更新を繰り返す（７）。
【選択図】図５

Description

【０００１】
【発明の属する技術分野】本発明は、有限要素法を用いた最適構造設計に使用される手法および装置である。
【０００２】
【従来の技術】連続体の位相最適化手法は１９８８年デンマーク工科大学のＭａｒｔｉｎＰｈｉｌｉｐＢｅｎｄｓｏｅとミシガン大学のＮｏｂｏｒｕＫｉｋｕｃｈｉによって提示された（ＧｅｎｅｒａｔｉｎｇＯｐｔｉｍａｌＴｏｐｏｌｏｇｉｅｓｉｎＳｔｒｕｃｔｕｒａｌＤｅｓｉｇｎｕｓｉｎｇａＨｏｍｏｇｅｎｉｚａｔｉｏｎＭｅｔｈｏｄ．，ＭａｒｔｉｎＰｈｉｌｉｐＢｅｎｄｓｏｅａｎｄＮｏｂｏｒｕＫｉｋｕｃｈｉ，ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓｉｎＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓａｎｄＥｎｇｉｎｅｅｒｉｎｇ７１（１９８８）１９７−２２４．：文献１）。この方法は対象構造物を含む大きな設計領域を無数の小さな穴が無限に空いた多孔質体と見做し、穴の空き方に周期性を仮定する。これに均質化法という応用数学の理論を適用して、ユニットセルと呼ばれる周期の最小単位の形状に周期境界条件を課して解析することにより、各有限要素でそれぞれ穴の状態に応じて等価になる均質な応力−歪マトリックスを算出する。そして、総量を制約された材料で最大の剛性を持つ形状を求める最適設計問題を、多孔質体の穴の大きさを制御することで、領域内に材料を最適に配置し、最大の剛性を持つ形状を求めた。彼らが用いた穴の形状は２次元問題で正方形であった。この手法は１９９１年ＫａｔｓｕｙｕｋｉＳｕｚｕｋｉとＮｏｂｏｒｕＫｉｋｕｃｈｉによって更に改良され、ユニットセルの穴の形状を、正方形を含む長方形にすることと、長方形の穴の長い方の辺を力のかかる方向に回転することにより、より効率の良い異方性材料の配置が可能となり、少ない繰り返し計算でより質の高い位相形状を算出することに成功した。（Ａｈｏｍｏｇｅｎｉｚａｔｉｏｎｍｅｔｈｏｄｆｏｒｓｈａｐｅａｎｄｔｏｐｏｌｏｇｙｏｐｔｉｍｉｚａｔｉｏｎ．，ＫａｔｓｕｙｕｋｉＳｕｚｕｋｉａｎｄＮｏｂｏｒｕＫｉｋｕｃｈｉ，ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓｉｎＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓａｎｄＥｎｇｉｎｅｅｒｉｎｇ９３（１９９１）２９１−３１８．：文献２）。しかしながら、複雑な問題に適用した場合に、必ずしも穴の大きさは空隙、または充填状態に分かれるとは限らず、中間の状態が多く残る場合や、チェッカーボードと呼ばれる空隙、充填の状態を交互に繰り返す位相形状を得ることもしばしばであった。図１に位相最適化によく用いられる梁の問題を示す。図１上の丸のマーク（１３）は対称拘束条件を表し、このモデルが上下に１／２対称であることを示す。図１左下の三角形に丸が２つついているマーク（１４）は、このモデルの左下端を垂直方向にはすべりを許容し、水平方向には拘束することを表している。図１の右上の矢印（１２）は右から水平方向に荷重が負荷されることを示す。この条件で、体積を領域全体の３０％に制約して最適化した結果が図２である。図２の下部に示すスケールは、ユニットセル内の材料の充填状態を空隙状態が０、充填状態を１としたときに、両者間を１０等分し、それを色の濃さで示したものである。薄い色の状態や、チェッカーボード状の模様が残っているのが見て取れる。この状態を改善するために、手法を踏襲しながらユニットセルの穴の形状を図４に示すような３次元の格子型にして、より鮮明な位相形状を得る手法も発明者の一人により提示されている（フレームベース・ユニットセルを用いた位相最適化の研究、石井惠三、青村茂、菊池昇、日本機械学会論文集（Ｃ編）６７巻６５４号、４９９−５０６：文献３）が、多くの場合、算出された結果に画像処理フィルターを適用しなければ、結果をそのまま利用することは困難である。
【０００３】
【発明が解決しようとする課題および目的】前記技術で算出した位相形状を各有限要素で残った材料と、充填された状態との比で密度表示すると、図２に示すように空隙状態と充填状態との中間的な密度や、チェッカーボードと呼ばれる充填状態と空隙状態が交互に配置される分布が頻繁に現れ、設計者が算出した位相形状を基に構造形状を決定するためには、複雑な画像処理フィルターを最適化の繰り返し処理の中に組み込むなどの高度の技量を必要とした。本発明は一切のフィルター処理を用いることなく従来の位相最適化の計算結果に現れる中間密度およびチェッカーボードの分布を極力なくし、この手法を利用する技術者に形状の認識がよりしやすいシステムを提供することを目的とする。
【０００４】
【課題を解決するための手段】上記目的を達成するために、本発明による位相最適化システムは、従来の技術では有限要素分割されたモデルの各有限要素で、要素内では一定となる均質化された材料定数を定義するのに対し、本発明によるシステムは各有限要素の頂点である節点において材料密度を定義し、各有限要素内では節点の密度値と、要素で定義された補間近似関数を用いて設計変数（密度値）が要素内で連続に分布することを仮定する。これは、要素内で一定となる設計変数を定義していた従来の手法は、要素間での設計変数の連続性を阻害し、これがチェッカーボードを生む原因になると予測したことによる。
【０００５】
本発明によるシステムは、構造物を含む設計領域内に最適な連続的材料分布を形成して構造物の最適位相形状をコンピュータによって生成するシステムであって、対象構造の有限要素モデルデータおよび位相最適化制御データの入力手段と、連続的材料分布を実現するために各有限要素の頂点である節点において設計変数として定義した材料密度を用いて応力−歪マトリックスデータを作成する手段と、前記算出された応力−歪マトリックスを用いて有限要素解析を行う手段と、目的関数、制約関数およびそれぞれの感度を計算する感度解析の手段と、前記有限要素解析および感度解析により算出された目的関数と予め定められた閾値との関係を評価し収束判定する手段と、収束判定の結果に応じて設計変数（節点における材料密度）を更新する手段と、前記材料密度を更新した後、前記応力−歪マトリックスを算出する手段と、前記有限要素解析を行う手段と、前記感度解析を行う手段と、前記収束判定する手段と、前記設計変数を更新する手段とを繰り返す手段を具備することを特徴とする最適位相形状を算出するシステムである。
【０００６】
更に、構造物を含む設計領域内に最適な連続的材料分布を形成して構造物の最適位相形状をコンピュータによって生成するシステムに、対象構造の有限要素モデルデータおよび位相最適化制御データを入力する手順と、連続的材料分布を実現するために各有限要素の頂点である節点において設計変数として定義した材料密度を用いて応力−歪マトリックスデータを作成する手順と、前記算出された応力−歪マトリックスを用いて有限要素解析を行う手順と、目的関数、制約関数およびそれぞれの感度を計算する感度解析の手順と、前記有限要素解析および感度解析により算出された目的関数と予め定められた閾値との関係を評価し収束判定する手順と、収束判定の結果に応じて設計変数である節点における材料密度を更新する手順と、前記材料密度を更新した後、前記応力−歪マトリックスを算出する手順と、前記有限要素解析を行う手順と、前記感度解析を行う手順と、前記収束判定する手順と、前記設計変数を更新する手順とを繰り返す手順と、を実行させるためのコンピュータプログラムである。
【０００７】
更に、構造物を含む設計領域内に最適な連続的材料分布を形成して構造物の最適位相形状をコンピュータによって生成するコンピュータプログラムを記録したコンピュータ読み取り可能な記憶媒体であって、システムに、対象構造の有限要素モデルデータおよび位相最適化制御データを入力する手順と、連続的材料分布を実現するために各有限要素の頂点である節点において設計変数として定義した材料密度を用いて応力−歪マトリックスデータを作成する手順と、前記算出された応力−歪マトリックスを用いて有限要素解析を行う手順と、目的関数、制約関数およびそれぞれの感度を計算する感度解析の手順と、前記有限要素解析および感度解析により算出された目的関数と予め定められた閾値との関係を評価し収束判定する手順と、収束判定の結果に応じて設計変数である節点における材料密度を更新する手順と、前記材料密度を更新した後、前記応力−歪マトリックスを算出する手順と、前記有限要素解析を行う手順と、前記感度解析を行う手順と、前記収束判定する手順と、前記設計変数を更新する手順とを繰り返す手順とを実行させるためのコンピュータプログラムを記憶したことを特徴とするコンピュータ読み取り可能な記憶媒体である。
【０００８】
【発明の実施形態】本発明の実施形態を図面により詳細に説明する。図５に本発明の最適位相形状を生成するプログラムのフローチャートを示す。
【０００９】
図５の（１．１）において、マイクロソフト社のオペレーティングシステム、Ｗｉｎｄｏｗｓを採用したパーソナルコンピュータで稼動する、ＥｌｅｃｔｒｏｎｉｃＤａｔａＳｙｓｔｅｍｓ社（アメリカ合衆国）の有限要素法プリポストプロセッサＦＥＭＡＰで作成した有限要素データ、即ち節点座標情報、要素結合情報、要素と材料定数の対応情報、要素と物理定数の対応情報、材料定数情報、物理定数情報、拘束条件情報、荷重条件情報を図５に示す本発明の最適位相形状を生成するプログラムに入力する。
【００１０】
図５の（１．２）において、要素毎の設計／非設計情報（設計する部分と設計しない部分を異なる番号で区別する）、制約体積（設計領域内に残す材料の体積が設計領域全体積に占める割合）、繰り返し計算を始めるステップ番号、繰り返し計算を終了させるステップ番号、収束判定に用いる閾値、設計変数の変動幅を図５に示す本発明の最適位相形状を生成するプログラムに入力する。尚、これらの値は設計の条件や設計者が解の精度を考慮して任意に決定できる。
【００１１】
図５の（２）において、設計変数として定義した要素内の任意の節点ｉにおける材料密度ρ_ｉを用いて、要素内の任意の場所における材料密度ρ（ｘ）を以下のように計算する。

ここでｍは要素を構成する節点の数、Ｎ_ｉは要素で定義された補間近似関数とする。
次に、密度に依存する応力−歪マトリックスは、材料が充填された状態の応力−歪マトリックスＤ_０と材料密度ρ（ｘ）のｐ乗の積と仮定して算出する。

ここで、Ｄ^Ｈは後の有限要素解析に用いる応力−歪マトリックス、ｘは領域内の任意の場所を示す座標、ρは設計変数として仮定した仮想的な材料密度で値は０（ゼロ）より大きく１より小さい。ｐは任意の正の実数で本発明では数値実験により２．５を用いた。
【００１２】
図５の（３）において、前項において算出された各節点における応力−歪マトリックスを用いて対象モデルを有限要素解析し、平均コンプライアンス（領域内の総歪みエネルギーの２倍の値）を計算する。
有限要素モデルは、節点の番号と座標値、要素を構成する節点番号のデーブル、モデルの材料定数、物理定数、拘束条件、節点および要素に負荷される荷重によって構成され、未知数である節点の変位ｕ（Ｕｘ，Ｕｙ，Ｕｚ）と要素内の歪εの関係を次のように定義する。

ここで、ｕは要素を構成する節点における仮想変位ベクトル、εは要素内の歪ベクトル、Ｂは要素で仮定した補間近似関数の微分で構成する歪−変位関係マトリックスである。また要素内における仮想変位ｕによって起こる歪εと応力σの関係は次のように決定される。

ここで、Ｄ_０は応力−歪関係マトリックスと呼ばれ、線形弾性問題では等方性材料の場合はヤング率とポアソン比のみで決定される。便宜上式（４）の右辺の歪εに掛かる係数をＤ^Ｈ＝ρ^ｐＤ_０と定義する。σは要素内応力ベクトルである。本発明では従来の手法が要素内で一定なＤ^Ｈを仮定していたのに対し、式（２）で表される各節点で計算される応力−歪マトリックスに、要素で定義した補間近似関数を適用して、要素内で連続したＤ^Ｈを仮定する。
仮想仕事の原理を適用して、ポテンシャルエネルギーを停留させる条件から次の関係式が導かれる。

ここでＫは各要素の剛性マトリックスを重ね合わせた全体剛性マトリックス、ｆは同様に要素の節点荷重成分を重ね合わせた全体荷重ベクトルとする。またΩは設計領域を示す。この平衡方程式を解いて得られた節点変位ベクトルｕから前述した関係に従って歪および応力を計算し、これらの値から目的関数である次式の平均コンプライアンスを計算する。

【００１３】
図５の（４．１）において、前項で算出された目的関数の感度は、各要素で次のように計算される。
目的関数である平均コンライアンスの節点ｉにおける設計変数ρ_ｉに関する感度は次式のようになる。

【００１４】
図５の（４．２）において、同様に制約関数ｇ（体積）の節点ｉにおける設計変数ρ_ｉに関する感度は次式のように計算される。

【００１５】
図５の（５）において、上記有限要素解析を行う手段により算出された目的関数とあらかじめ定められた閾値との関係を評価し収束判定する。
第ｎ回目に計算された目的関数値（平均コンプライアンス）をΠ_ｎ、第ｎ−１回目に計算された目的関数値をΠ_ｎ−１とし、最適化制御データとして入力した収束判定の閾値をｒとすると、収束判定は次のように行う。

式（９）が満足されるか、最適化制御データとして入力した繰り返し計算を終了させるステップ番号に早く到達したほうの結果を採用する。
【００１６】
図５の（６）は、図５の（５）において収束しなかった場合の処理を示す。即ち、収束しなかった場合は、最適化規準法、凸線形化最適化法、逐次線形化法などの数値最適化手法により節点の密度値を更新する。
本発明では最適化規準法を使用した。最適化規準法は次のような手順で行う。１）ラグランジュ未定乗数Λを導入してラグランジュ関数を次のように定義する。

ここで、ρは設計変数（仮想的な材料密度）、Πは目的関数（ここでは平均コンプライアンス）、Ｖ（ρ）は仮想的な材料密度による設計領域の体積、Ｖｃは最適化制御データとして入力した制約体積を示す。
２）このラグランジュ関数が停留することが最適化の必要条件になるから、式（１０）の変分を取ると、

となる。δρ，δΛの任意性から次式が成り立つ。

これを利用して、設計変数ａ_ｉおよびラグランジュ未定乗数Λに関し以下の更新ルールを適用する。

ここでηは適当な減速パラメータで、数値実験により０．７５とした。ｋは計算繰り返し回数を示す。
式（１３）および式（１４）中の右辺の括弧内の値が繰り返し計算の中で常に１になるように設計変数、ラグランジュ未定乗数を更新する。但し、算出した設計変数は上限値が１．０、下限値は０．０と決められているので、どちらかに達した場合は、上下限値を採用する。さらに、設計変数の変動幅が大きくなると収束が不安定になることから、変動幅にも最適化制御データで入力をした値（ｍｏｖｉｎｇｂｏｕｎｄ）で制約を与える。本発明では０．１５を用いた。
【００１７】
図５の（６）の処理、即ち設計変数を更新した後に、図５の（７）に示すように、（２）、（３）、（４）、（５）および（６）を収束判定が収束となるか、繰り返し計算を終了させるステップ番号に達するまで繰り返し、最適な位相形状を得る。
【００１８】
本発明の有効性を示すために２つの具体例をあげる。
（実施例１：）
図１に示す垂直方向１５０．０ｍｍ、水平方向２５．０ｍｍ、厚さ１．０ｍｍの梁（１１）の対称性を考慮した１／２モデルの対称境界上（１３）に荷重（１２）ｐ＝１．０Ｋｇｆが負荷された場合の位相最適化を行なう。拘束条件は、対称境界上に垂直方向に拘束を与え、左端下（１４）に水平方向の拘束を与える。材料定数はヤング率Ｅ＝２１０００．０Ｋｇｆ／ｍｍ^２、ポアソン比ｖ＝０．３を使用した。モデルの有限要素分割は横×縦×奥行き：７５×２５×１、即ち要素辺長１ｍｍの立方体要素を用いてモデルを構築した。
最適化の条件は、目的関数として平均コンプライアンスの最小化、制約条件は領域に残す材料を全領域の３０％とする。また、設計／非設計情報は全領域を設計の対象とした。繰り返し計算を終了させるステップ番号は５０とし収束判定誤差は１０^−６とした。
図２に従来の要素毎に設計変数を定義し、応力−歪マトリックスの作成にＳｕｚｕｋｉ、Ｋｉｋｕｃｈｉ（文献２）の長方形の穴を３次元の直方体に拡張したユニットセルを用い、均質化法を適用して計算した結果を示す。図２の下部に示すスケールは穴を除いた材料の占める体積を、元の要素の体積で除した体積比、即ち要素の材料充填率として、空隙から充填までを０から１の間に正規化した値を１０個のランクに色分けしたものである。
そして図３に本発明によるシステムで計算した結果を示す。図３の下部に示すスケールは要素を構成する各節点において計算された設計変数（密度値）の要素内平均値を１０個のランクに色分けしたものである。
両者を見比べると本発明によるシステムで求めた位相形状は明らかに従来のものより簡素な形状で、しかも中間値やチェッカーボード状の位相形状も残っていない。
【００１９】
（実施例２：）
図６に示すミッチェルトラスと呼ばれる、横４０．０ｍｍ、縦５５．０ｍｍ、厚さ１．０ｍｍの直方体の板（６１）の上端から垂直方向に１５．０ｍｍ、左端から横方向に２０．０ｍｍの位置に中心を持つ半径１０．０ｍｍの穴をあけ、穴の周上（６２）（図６中の太い輪で示す部分）を拘束する。そして下部左端から２０．０ｍｍの板中心に１．０Ｋｇｆの荷重（６３）を左向きに（図６中矢印の方向に）載荷する。材料定数はヤング率Ｅ＝２１０００．０Ｋｇｆ／ｍｍ^２、ポアソン比ｖ＝０．３を使用した。
最適化の条件は、目的関数は平均コンプライアンスの最小化、制約条件は領域に残す材料を全領域の１５％とする。また、設計／非設計情報は全領域を設計の対象とした。繰り返し計算を終了させるステップ番号は５０とし収束判定誤差は１０^−６とした。
図７に従来の手法により、要素で設計変数を定義し、応力−歪マトリックスの作成にＳｕｚｕｋｉ、Ｋｉｋｕｃｈｉ（文献２）の長方形の穴を３次元の直方体に拡張したユニットセルを用い、均質化法を適用して計算した結果を示す。図７の下部に示すスケールは穴を除いた材料の占める体積を、元の要素の体積で除した要素の材料充填率である。そして図８に本発明の手法で計算した結果を示す。図８の下部に示すスケールは要素を構成する各節点において計算された設計変数（密度値）の要素内平均値を１０個のランクに色分けしたものである。
【発明の効果】前記２種類の実施例において、図２および図７を見ても明らかなように、チェッカーボード状の位相形状が算出され、これらの結果を設計に利用するためには複雑なフィルター処理が避けられない。しかし，本発明によればフィルター処理を全く用いなくとも、従来の手法を使用した場合に比べ、位相最適化の結果に現れる中間密度やチェッカーボードを作らず、認識しやすい形状を得ることができる。
【図面の簡単な説明】
【図１】図１は梁の１／２対称モデルの有限要素分割と境界条件である。
【図２】図２は図１に示すモデルを従来の手法で計算して得られた位相形状である。
【図３】図３は図１に示すモデルを本発明による手法で計算して得られた位相形状である。
【図４】図４は文献３で使用されたユニットセル形状である。
【図５】図５は本発明によるシステムのフローチャートである。
【図６】図６はミッチェルトラスと呼ばれるモデルの有限要素分割と境界条件である。
【図７】図７は図６に示すモデルを従来の手法で計算して得られた位相形状である。
【図８】図８は図６に示すモデルを本発明による手法で計算して得られた位相形状である。
【符号の説明】
１１梁モデル
１２梁モデルに載荷される荷重
１３梁モデルの１／２対称境界の拘束条件
１４梁端部の拘束条件
６１ミッチェルトラスモデル
６２ミッチェルトラス円周部の拘束条件
６３ミッチェルトラスに載荷される荷重

Claims

構造物を含む設計領域内に最適な連続的材料分布を形成して構造物の最適位相形状をコンピュータによって生成するシステムであって、対象構造の有限要素モデルデータおよび位相最適化制御データの入力手段と、連続的材料分布を実現するために各有限要素の頂点である節点において設計変数として定義した材料密度を用いて応力−歪マトリックスデータを作成する手段と、前記算出された応力−歪マトリックスを用いて有限要素解析を行う手段と、目的関数、制約関数およびそれぞれの感度を計算する感度解析の手段と、前記有限要素解析および感度解析により算出された目的関数と予め定められた閾値との関係を評価し収束判定する手段と、収束判定の結果に応じて設計変数である節点における材料密度を更新する手段と、前記材料密度を更新した後、前記応力−歪マトリックスを算出する手段と、前記有限要素解析を行う手段と、前記感度解析を行う手段と、前記収束判定する手段と、前記設計変数を更新する手段とを繰り返す手段を具備することを特徴とする最適位相形状を生成するシステム。
構造物を含む設計領域内に最適な連続的材料分布を形成して構造物の最適位相形状をコンピュータによって生成するシステムに、対象構造の有限要素モデルデータおよび位相最適化制御データを入力する手順と、連続的材料分布を実現するために各有限要素の頂点である節点において設計変数として定義した材料密度を用いて応力−歪マトリックスデータを作成する手順と、前記算出された応力−歪マトリックスを用いて有限要素解析を行う手順と、目的関数、制約関数およびそれぞれの感度を計算する感度解析の手順と、前記有限要素解析および感度解析により算出された目的関数と予め定められた閾値との関係を評価し収束判定する手順と、収束判定の結果に応じて設計変数である節点における材料密度を更新する手順と、前記材料密度を更新した後、前記応力−歪マトリックスを算出する手順と、前記有限要素解析を行う手順と、前記感度解析を行う手順と、前記収束判定する手順と、前記設計変数を更新する手順とを繰り返す手順と、を実行させるためのコンピュータプログラム
構造物を含む設計領域内に最適な連続的材料分布を形成して構造物の最適位相形状をコンピュータによって生成するコンピュータプログラムを記録したコンピュータ読み取り可能な記憶媒体であって、システムに、対象構造の有限要素モデルデータおよび位相最適化制御データを入力する手順と、連続的材料分布を実現するために各有限要素の頂点である節点において設計変数として定義した材料密度を用いて応力−歪マトリックスデータを作成する手順と、前記算出された応力−歪マトリックスを用いて有限要素解析を行う手順と、目的関数、制約関数およびそれぞれの感度を計算する感度解析の手順と、前記有限要素解析および感度解析により算出された目的関数と予め定められた閾値との関係を評価し収束判定する手順と、収束判定の結果に応じて設計変数である節点における材料密度を更新する手順と、前記材料密度を更新した後、前記応力−歪マトリックスを算出する手順と、前記有限要素解析を行う手順と、前記感度解析を行う手順と、前記収束判定する手順と、前記設計変数を更新する手順とを繰り返す手順とを実行させるためのコンピュータプログラムを記憶たことを特徴とするコンピュータ読み取り可能な記憶媒体