(57)【要約】【課題】  コンピュータを用いて構造最適化を実行する場合に於いて、支持構造のない構造体の最適化が可能な手法を提供すること。【解決手段】  本発明の構造最適化演算装置は、設計変数により構造が決定される構造体の剛性行列の擬似逆行列を用いて、荷重作用領域に於ける局所コンプライアンス行列の固有値及び固有ベクトルと、固有ベクトルを荷重ベクトルとしたときの構造体に於ける変位ベクトルとを算出し、設計変数に対する固有値の感度を用いて、設計領域に対する構造体の体積制約割合を満たした状態にて固有値が最小化されるよう設計変数を更新する各処理を所定の終了条件が成立するまで反復し、設計変数により決定される構造体の最新の構造を出力又は表示する。固有値の感度は、固有値が大きいほど重みが大きくなる固有値の重み付き和の感度であり、設計変数は、固有値の重み付き和が最小化されるよう更新されてよい。【選択図】        図４
2020016936.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを用いた数値演算処理により設計領域内にて設定された条件の下で設計変数にて構造が決定される構造体にして支持構造のない構造体の最適な構造を見出す構造最適化演算装置であって、
前記設計領域の寸法、前記構造体に於いて荷重の作用する荷重作用領域の範囲、前記構造体を構成する材料の弾性特性を表す物性定数及び前記設計領域の体積に対する前記構造体の体積の割合である体積制約割合を含む条件を設定する条件設定手段と、
前記設計変数を用いて前記構造体の剛性行列を算出する剛性行列算出手段と、
前記剛性行列の擬似逆行列を算出する擬似逆行列算出手段と、
前記擬似逆行列を用いて前記構造体の前記荷重作用領域に於ける局所コンプライアンス行列を算出する局所コンプライアンス行列算出手段と、
前記局所コンプライアンス行列の少なくとも一つの固有値と対応する固有ベクトルとを算出する固有値解析手段と、
前記固有ベクトルと前記擬似逆行列とを用いて前記固有ベクトルを荷重ベクトルとしたときの前記構造体に於ける変位ベクトルを算出する変位ベクトル算出手段と、
前記変位ベクトルを用いて前記設計変数に対する前記少なくとも一つの固有値の感度を算出する固有値感度算出手段と、
前記体積制約割合を満たした状態で前記感度を用いて前記少なくとも一つの固有値が最小化されるよう前記設計変数を更新する設計変数更新手段と、
所定の終了条件が成立するまで、前記更新された設計変数を用いて、前記剛性行列算出手段、前記擬似逆行列算出手段、前記局所コンプライアンス行列算出手段、前記固有値解析手段、前記変位ベクトル算出手段、前記固有値感度算出手段及び前記設計変数更新手段の各処理を反復させる反復処理手段と、
前記設計変数により決定される構造体の最新の構造を出力又は表示する構造体出力手段と
を含む装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピュータを用いた構造体の構造最適化演算を実行する装置に係り、より詳細には、コンピュータによる数値解析を通じて、任意の設計条件・境界条件下に於いて最適な構造（所与の力学的条件下にて剛性が最大化される構造）を見出す装置に係る。
【背景技術】
【０００２】
コンピュータを用いた有限要素法などの数値計算技術の発展により、任意の構造体の、任意の設計条件・境界条件下に於ける最適な構造又は形状をコンピュータ上で求める「構造最適化」演算が実施されるようになっている。かかるコンピュータを用いた構造最適化の手法によれば、車両のボディーを初め、種々の機械器具、建築物などの種々の構造体を実際に製作する前に、構造体が要求される特性を備えているか否かの把握、或いは、要求される若しくは好ましい性能を備えた構造体の設計案の取得が可能となる点で有利である。
【０００３】
そのような「構造最適化」の一つとしては、「トポロジー最適化（位相最適化）」と称される手法が知られている（特許文献１、２、非特許文献１、２）。かかる「トポロジー最適化」では、端的に述べれば、有限要素法等による数値解析とその結果から得られる平均コンプライアンス等の目的関数又は評価関数を最適化する演算処理とを通じて、或る与えられた領域（設計領域）内に於いて、任意の設計条件・境界条件の下での最適な構造、例えば、任意の力学的条件下での剛性が最大化される構造、を与える材料分布が生成される。より具体的には、例えば、特許文献１では、初めに、多数の、それぞれ微小な孔を有する有限要素から構成される任意の寸法の設計領域を設定した上で、任意の境界条件と体積制約条件の下、各有限要素の孔の寸法を変化させながら設計領域内にて有限要素により形成される構造体の平均コンプライアンスとその（各要素の孔の大きさに対する）感度を反復して演算し、最小の平均コンプライアンスを与える有限要素内の孔の寸法の分布を見出し、剛性の最大化された構造体の形状を生成する手法が開示されている。特許文献２は、材料の密度がそれぞれ割り当てられた多数の有限要素から構成される任意の寸法の設計領域に於いて、任意の境界条件と体積制約条件の下、各有限要素の材料密度を変化させながら有限要素により形成される構造体の平均コンプライアンスとその（各要素の密度に対する）感度を反復して演算し、平均コンプライアンスを最小化する材料の密度分布を生成し、剛性の最大化された構造体の形状を生成する手法を開示している。特許文献３は、任意の構造体に於いて外力を受ける部分の構造最適化を行うべく、構造最適化の対象となる部分と構造体のその他の部分との間で荷重が伝達されるように構造体の数理モデルを設定することを提案している。
【０００４】
また、任意の構造体のコンプライアンスを最小化するように構造体の形状を変化させる形状最適化演算も種々知られている。例えば、特許文献４では、形状最適化演算に於ける計算の処理負担を低減させるべく、フェーズフィールド法を適用して任意の構造体の形状の最適化演算を実行する構成が提案されている。
【０００５】
更に、構造体に対して任意の荷重が作用した場合、即ち、構造体に対して作用する荷重が不定である場合に剛性が最大化される最適構造を求める「ロバスト構造最適化」も提案されている（非特許文献３、４）。かかるロバスト構造最適化に於いては、構造体の或る領域に作用する荷重の値と分布が不定な場合の最適な構造として、想定される最も危険な荷重（最悪荷重）に対して剛性を最大化する構造体の構造が演算により求められる。より具体的には、最悪荷重が付与されたときのコンプライアンスと荷重の分布は、構造体に於いて荷重の作用する領域（荷重作用領域）内の荷重分布を表すベクトルＦｌと荷重作用領域の変位分布を表す変位ベクトルＵｌとを関係づける方程式：
Ｕｌ＝Ｃ・Ｆｌ      …（１）
に於ける局所的なコンプライアンス行列Ｃの（荷重ベクトルＦｌのノルム||Fl||＝１のときの）最大固有値とこれに対応する固有ベクトルにより与えられることが理論的に見出されているので、ロバスト構造最適化では、実際には、式（１）の局所コンプライアンス行列Ｃの最大固有値を最小化するように、構造体に於ける材料分布や形状の最適化演算が実行される。このロバスト構造最適化演算によれば、任意の構造体に於いて、任意の領域に最も変形を大きくする最悪荷重が作用した場合に最大の剛性を呈し変形が最も小さくなる構造が算定できることとなるので、如何なる分布にて荷重が作用したしても高い剛性を呈する最適な構造が提供できることが期待される。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００２−１８９７６０
【特許文献２】特開２００４−３４８６９１
【特許文献３】特開２０１４−１４９７３３
【特許文献４】特開２０１０−１０８４５１
【非特許文献】
【０００７】
【非特許文献１】計算力学レクチャーコース  トポロジー最適化  一般社団法人日本計算工学会編  西脇眞二・泉井一浩・菊池昇著  丸善出版  ２０１３年
【非特許文献２】カイ・リュウ、アンドレス・トバル（Kai  Liu  Andr´es  Tovar）  ストラクト・マルチディスク・オプティム（Struct  Multidisc  Optim）  2014年  50:1175−1196
【非特許文献３】新居他３名、日本機械学会論文集（Ａ編）77巻775号  2011年  ４７２−４８２頁
【非特許文献４】下田他３名、日本機械学会論文集81巻832号(DOI:10.1299/transjsme.15-00353)  2015年
【発明の概要】
【発明が解決しようとする課題】
【０００８】
ところで、現実の任意の構造体に荷重が作用する状態を考えたとき、構造体を支持する構造或いは構造体を固定する支点の位置などが不定であることが多い。例えば、ユーザーの手からの荷重を受ける携帯電話の外殻など、実際の構造体は、どのように保持或いは支持されるかは不定であり、そのような場合、構造体に於ける支点の位置は不定となる。
【０００９】
しかしながら、上記の如き従前のロバスト構造最適化又はその他の構造最適化の方法に於いては、構造体に於いて支持構造或いは支点を明確に規定する拘束条件を与えなければ、演算を実行できないようになっている。上記の如き構造最適化演算では、端的に述べれば、最適構造を見出すために必要な最小化されるべきコンプライアンス若しくはその他の評価関数又はその感度（材料分布又は形状の変化に対する評価関数の変化、設計変数についての評価関数の導関数）の決定のために、構造体に於ける変位ベクトルＵと荷重ベクトルＦとを関連付ける剛性方程式：
ＫＵ＝Ｆ      …（２）
に於ける剛性行列Ｋの逆行列Ｋ-1を用いて、上記の局所コンプライアンス行列Ｃの算出、並びに、変位ベクトルＵの算出が必要となる。ところが、構造体を支持する構造又は支点を与える拘束条件がない場合には、剛性行列Ｋの逆行列Ｋ-1が存在せず、局所コンプライアンス行列Ｃと変位ベクトルＵの演算が実行できないこととなっている。これは、物理的には、構造体を支持する構造がない状態では、構造体が剛体として自由に平行又は回転運動し得ることとなり（以下、「剛体運動」と称する。）、構造体内の各点の変位が式（２）の解として一意に決定されないことに相当する。
【００１０】
この点に関し、本発明の発明者が研究したところ、支持構造又は支点のない構造体モデルの場合でも、局所コンプライアンス行列Ｃの算出並びに変位ベクトルＵの算出に於いて、剛性行列Ｋの逆行列Ｋ-1に代えて、剛性行列Ｋの「擬似逆行列（ムーア・ペンローズの逆行列または一般逆行列とも称される。）」Ｋ＋を用いると、構造体に於ける剛体運動を除いた変位ベクトル、かかる変位ベクトルにより生ずるコンプライアンス及びその感度を算出することができ、最小化されたコンプライアンスを与える最適な性能を有する構造体の構造又は形状が得られること、即ち、構造の最適化が可能であることが見出された。かかる知見は、本発明に於いて利用される。
【００１１】
かくして、本発明の課題は、コンピュータを用いた数値解析技術により、任意の構造体に任意の荷重が付与される場合に構造体の剛性が最大化される最適な構造を見出す構造最適化演算を実行する場合に於いて、支持構造又は支点のない構造体の構造最適化が可能な手法を提供することである。
【課題を解決するための手段】
【００１２】
本発明によれば、上記の課題は、コンピュータを用いた数値演算処理により設計領域内にて設定された条件の下で設計変数にて構造が決定される構造体にして支持構造のない構造体の最適な構造を見出す構造最適化演算装置であって、
前記設計領域の寸法、前記構造体に於いて荷重の作用する荷重作用領域の範囲、前記構造体を構成する材料の弾性特性を表す物性定数及び前記設計領域の体積に対する前記構造体の体積の割合である体積制約割合を含む条件を設定する条件設定手段と、
前記設計変数を用いて前記構造体の剛性行列を算出する剛性行列算出手段と、
前記剛性行列の擬似逆行列を算出する擬似逆行列算出手段と、
前記擬似逆行列を用いて前記構造体の前記荷重作用領域に於ける局所コンプライアンス行列を算出する局所コンプライアンス行列算出手段と、
前記局所コンプライアンス行列の少なくとも一つの固有値と対応する固有ベクトルとを算出する固有値解析手段と、
前記固有ベクトルと前記擬似逆行列とを用いて前記固有ベクトルを荷重ベクトルとしたときの前記構造体に於ける変位ベクトルを算出する変位ベクトル算出手段と、
前記変位ベクトルを用いて前記設計変数に対する前記少なくとも一つの固有値の感度を算出する固有値感度算出手段と、
前記体積制約割合を満たした状態で前記感度を用いて前記少なくとも一つの固有値が最小化されるよう前記設計変数を更新する設計変数更新手段と、
所定の終了条件が成立するまで、前記更新された設計変数を用いて、前記剛性行列算出手段、前記擬似逆行列算出手段、前記局所コンプライアンス行列算出手段、前記固有値解析手段、前記変位ベクトル算出手段、前記固有値感度算出手段及び前記設計変数更新手段の各処理を反復させる反復処理手段と、
前記設計変数により決定される構造体の最新の構造を出力又は表示する構造体出力手段と
を含む装置によって達成される。
【００１３】
上記の構成に於いて、「構造体」とは、車両のボディーを初め、種々の機械器具、建築物又はそれらを構成する任意の部材などとして使用される構造体であってよい。「構造体の最適な構造を見出す」とは、通常の構造最適化演算、例えば、トポロジー最適化演算、形状最適化演算などの場合と同様に、構造体に荷重が付与される場合の構造体の変形の程度が最小となる構造、より具体的には、構造体のコンプライアンスが最小となる構造を見出すことを意味している（なお、本明細書に於いて、「構造体の構造」というときには、構造体の形状、トポロジー或いは材料密度の分布をいうものとする。）。従って、上記の「設計領域」は、通常の構造最適化演算と同様に設定されてよく、構造体の構造を決定する「設計変数」は、構造最適化演算の態様によって適宜選択されてよい。構造最適化演算の手法としては、具体的には、例えば、所謂「密度法」を用いたトポロジー最適化演算（例えば、特許文献２、非特許文献１、２）が採用されてよく、その場合には、「設計変数」は、「設計領域」を任意に格子状に分割して得られる複数の節点を結んで構成される多面体（又は多角形）の要素に於いて割り当てられる０〜１の値を取る「材料密度」であり、演算過程に於いて更新されて、好ましくは、０又は１に収束していくことが期待される変数であってよい。そして、「設計変数」が「材料密度」である場合には、材料密度が実質的に非０である要素（例えば、小数点第一位を四捨五入して１になる要素）の集合が成す形状が、「設計領域」内にて生成される目的の「構造体」の形状となる。また、設計変数は、構造体の輪郭を画定する線又は面を決定する変数であってもよい。
【００１４】
上記の本発明の装置に於ける「条件設定手段」では、設計領域の寸法、構造体に於いて荷重の作用する領域（荷重作用領域）の範囲、構造体を構成する材料の弾性特性を表す物性定数、設計領域の体積に対する構造体の体積の割合である体積制約割合を含む構造最適化演算に必要な初期条件が設定される。しかしながら、特に本発明の装置に於いては、最適化の対象となる構造体が、支持構造のない構造体、即ち、外部からの拘束を受けず、自由な剛体運動が許される状態にある構造体であるので、従前の構造最適化演算に際して通常設定される構造体に於ける変位に対する拘束条件又は境界条件は一切設定されなくてよいことは理解されるべきである。また、荷重条件について、構造体に於ける荷重作用領域の具体的な範囲は指定されるが、荷重作用領域内に於ける荷重の値或いは荷重の分布は、従前のロバスト構造最適化演算の場合と同様に不定のままであってよく、これらの条件は、後述の固有値解析演算の結果によって設定される。構造体を構成する材料の弾性特性を表す物性定数は、典型的には、ヤング率とポアソン比であってよいが、これに限定されない。「体積制約割合」は、最適化演算に於ける拘束条件となる（構造体の運動を拘束するものではない。）。
【００１５】
上記の装置に於いて、
「剛性行列算出手段」にて算出される構造体の「剛性行列」は、上記の式（２）に於ける剛性行列Ｋと同様のものであってよく、具体的には、構造体に格子状に設定される複数の要素に於ける点（節点）の各々に作用する荷重を成分とする荷重ベクトルと各節点に於ける変位を成分とする変位ベクトルとを関係付ける剛性行列であってよい。各要素の弾性は、材料の弾性特性を表す物性定数と各要素の設計変数（例えば、材料密度）によって決定されるので、剛性行列は、各要素の設計変数の関数となり、本発明による最適化演算の過程で設計変数が更新されると共に剛性行列の各成分の値も変化することとなる。
「擬似逆行列算出手段」にて算出される構造体の剛性行列の「擬似逆行列」は、既に述べた如き、ムーア・ペンローズの逆行列または一般逆行列とも称される行列であり、剛性行列の成分を用いて所定のアルゴリズムにより一義的に算出される。具体的には、擬似逆行列は、プログラム言語に用意された関数コマンドを利用して算出されてよく、或いは、任意の近似方法を用いて算出することも可能である。
「局所コンプライアンス行列算出手段」にて算出される「局所コンプライアンス行列」とは、上記の式（１）に於ける局所コンプライアンス行列Ｃと同様のものであり、具体的には、荷重作用領域内に於ける各節点の荷重を成分とする局所的な荷重ベクトルと各節点の変位を成分とする局所的な変位ベクトルとを関係づける行列である。かかる局所コンプライアンス行列Ｃは、後の実施形態の欄にて詳細に説明される如く、本来は、構造体の全節点の荷重を表す荷重ベクトルと荷重作用領域内の局所的な荷重ベクトルとを関連付ける変換行列と構造体の剛性行列の逆行列とを用いて算出されるところ、本発明の装置に於いては、既に述べた如く、支持構造のない構造体が最適化の対象となっているので、構造体の剛性行列の逆行列が存在していない。そこで、「局所コンプライアンス行列算出手段」では、構造体の剛性行列の逆行列に代えて、構造体の剛性行列の擬似逆行列を用いて、局所コンプライアンス行列が算出される。
【００１６】
次いで、「固有値解析手段」に於いては、上記の局所コンプライアンス行列の固有値と固有ベクトルとが算出されるところ、かかる算出処理は、任意の態様にて実行されてよく、典型的には、プログラム言語に用意された関数コマンドが利用されてよい。なお、好適には、固有値解析手段に於いて得られる固有ベクトルはそのノルム＝１のものが算出される（以下、特に断らない限り、固有ベクトルは、ノルム＝１のベクトルであるものとする。）。局所コンプライアンス行列の固有値と固有ベクトルは、典型的には、複数組算出され、本発明の装置に於いては、最大固有値とその固有ベクトルだけではなく、値の大きいものから複数の組の固有値と固有ベクトルが、後の処理に用いられてよい。なお、ここで得られた固有値と固有ベクトルは、それぞれ、荷重作用領域内の変位ベクトルが荷重ベクトルに平行な場合のコンプライアンスと荷重ベクトルを表しており、固有値が大きいほど、対応する固有ベクトルが荷重ベクトルとして作用したときのコンプライアンスが大きく、変形の程度が大きくなることを示している。また、後の実施形態の欄に於いて説明される如く、支持構造のない構造体の剛性行列の擬似逆行列を用いて算出される局所コンプライアンス行列の固有ベクトルは、構造体に於いて荷重の合計とモーメントの合計とが０となり、剛体運動を発生しない荷重ベクトルを与え、局所コンプライアンス行列の固有値は、そのときのコンプライアンスとなる。
【００１７】
更に、上記の本発明の装置の「変位ベクトル算出手段」では、上記の固有ベクトルと擬似逆行列とを用いて構造体全体に於ける各節点の変位を表す「変位ベクトル」算出される。かかる「変位ベクトル」は、本来であれば、式（２）に於ける荷重ベクトルＦに固有ベクトルを代入し、式（２）の解Ｕを演算することにより、即ち、
Ｕ＝Ｋ−１Ｆ      …（３）
を演算することにより算出されるが、既に述べた如く、本発明に於いては、式（２）の剛性行列Ｋの逆行列Ｋ−１が存在しない。そこで、本発明では、剛性行列Ｋの逆行列Ｋ−１に代えて、剛性行列Ｋの疑似逆行列Ｋ＋を用いて、
Ｕ＝Ｋ＋Ｆ      …（４）
を演算することにより、Ｕとして「変位ベクトル」が算出される。なお、「変位ベクトル」は、局所コンプライアンス行列の複数得られた固有ベクトルのそれぞれに対して、算出されてよい。荷重ベクトルから変位ベクトルを算出する具体的な演算手法は、剛性行列の逆行列に代えて擬似逆行列を使用することの他は、任意の手法にて、例えば、この分野に於いて通常利用されている有限要素法の手法にて実行されてよい。上記の如く局所コンプライアンス行列Ｃの固有ベクトルを荷重ベクトルとして疑似逆行列Ｋ＋を用いて算出された変位ベクトルＵは、構造体に於ける剛体運動を除いた構造体内に於いて対称的に生ずる変位を表すこととなる。
【００１８】
そして、上記の装置に於いて、変位ベクトルを用いて「固有値感度算出手段」にて設計変数に対する固有値、即ち、構造体のコンプライアンスの感度が算出され、その感度を用いて「設計変数更新手段」にて、構造体のコンプライアンスである固有値が最小化されるよう設計変数が更新される。即ち、本発明の装置に於いて、「固有値」（又はそれらの線形和）が構造最適化演算に於ける「評価関数」となる。また、ここにおいて、「感度」とは、設計変数の変化に対する固有値の変化量、即ち、設計変数についての固有値の微分に相当する量であり、変位ベクトルと、設計変数についての剛性行列の成分の微分値とから算出可能である。設計変数の更新は、従前の構造最適化演算に於いて用いられている手法（例えば、最適性基準法、逐次線形計画法、逐次凸関数近似法など）と同様であってよく、設計変数は、評価関数である固有値の感度を用いて、構造体が、初めに設定された体積制約割合（設計領域内にて最終的に形成される構造体の体積の、設計領域の体積に対する割合）を満たすことを数学的な拘束条件として、構造体の剛性が最大化されるよう更新されることとなる。
【００１９】
上記の本発明の装置に於いて、反復処理手段による更新された設計変数を用いた、剛性行列算出手段、擬似逆行列算出手段、局所コンプライアンス行列算出手段、固有値解析手段、変位ベクトル算出手段、固有値感度算出手段及び設計変数更新手段の各処理の反復は、所定の終了条件が成立するまで実行される。ここで、「所定の終了条件」は、目的の構造体の形状が得られるように任意に設定されてよいところ、具体的には、例えば、（ｉ）反復処理回数が任意に設定されてよい回数に到達した場合、（ii）構造体にて算出される評価関数（固有値又はそれらの線形和）の値が収束した場合（処理サイクル毎の評価関数の変化量が所定の小数以下となった場合）、（iii）更新される設計変数の値が収束した場合（処理サイクル毎の設計変数の変化量の最大値が所定の小数以下となった場合）などが選択されてよい。更に、本発明の装置に於いては、上記の如く「構造体出力手段」が設けられ、演算処理過程に於いて得られた設計変数により表される構造体の構造又は形状（例えば、材料密度の分布、有意な材料密度を有する有限要素の分布）がコンピュータのモニターの画面上や紙面にてグラフィックス表示されるようになっていてよい。
【００２０】
かくして、上記の本発明の装置では、基本的な構成としては、従前のロバスト構造最適化演算の場合と同様に、設計領域内に於ける構造体の構造を決定する設計変数を逐次的に更新しながら、局所コンプライアンス行列の固有値解析により構造体のコンプライアンスに相当する局所コンプライアンス行列の固有値とその固有値を与える固有ベクトルを決定し、その固有ベクトルを荷重ベクトルとして有限要素法等によって各要素の節点の変位ベクトルを算出し、かかる変位ベクトルを用いて評価関数である局所コンプライアンス行列の固有値が最小化されるように設計変数の値を探索し、そこで得られた設計変数にて表される形状が目的の構造体の形状として得られることとなる。かかる一連の演算処理を実行する構成に於いて、最適化の対象となる構造体が支持構造のない構造体である場合、構造体の剛性行列の逆行列がなく、局所コンプライアンス行列の算出及び変位ベクトルの算出ができないところ、本発明の装置では、後により詳細に説明される如く、構造体に支持構造のない場合でも算出することのできる構造体の剛性行列の擬似逆行列を用いると、支持構造のない構造体に於いて剛体運動を除いた変位ベクトル、コンプライアンス及びその感度を算出することができるという知見に基づき、局所コンプライアンス行列の算出及び変位ベクトルの算出に於いて、擬似逆行列を用い、これにより、支持構造のない構造体に於いて剛体運動を除いた構造体内で生ずる変位によるコンプライアンスを最小化する構造最適化演算の実行を可能にしている。
【００２１】
ところで、上記の如く、局所コンプライアンス行列の固有値と固有ベクトルは複数組得られるところ、最適化の対象とされた構造体に支持構造があり、構造体の剛体運動が制限されている場合には、構造体の構造が最適化されていく過程で局所コンプライアンス行列の最大固有値は、常に或る一つの固有ベクトルにより与えられ、その一つの固有ベクトルが連続的に変化することとなる。即ち、支持構造のある構造体の場合には、構造最適化の過程に於いて、複数の固有値と固有ベクトルの組のうちで最大固有値を与える固有ベクトルが交代することはなく、最初に固有値を最大にした固有ベクトルが、その後の反復演算の間も、対応する固有値が最大となったまま連続的に変化することとなる。従って、変位ベクトルの算出が反復して実行される間に於いて、常に最大固有値の固有ベクトルを荷重ベクトルとして選択しておけば、算出される変位ベクトル及びそれを用いて算出される評価関数の感度も連続的に変化し、その感度を用いて設計変数の更新を実行することにより、構造体の構造は連続的に変化して最適な構造に収束することが可能となっていた。
【００２２】
一方、支持構造のない構造体の場合には、構造最適化の過程に於いて、局所コンプライアンス行列の複数の固有値と固有ベクトルの組のうちで最大固有値を与える固有ベクトルの交代が生じ得ることが明らかとなった。そして、かかる固有ベクトルの交代が生ずると、最大固有値を与える固有ベクトルが不連続に変化し、構造体の構造が最適な構造へ収束しづらくなる状況が発生した。そこで、本発明の装置に於いては、局所コンプライアンス行列の複数の固有値と固有ベクトルの組のうち、固有値の大きい方から複数個の組を選択し、それらの選択された固有値が全体的に最小化されるように設計変数の更新が実行されるようになっていてよい。そのために、本発明の装置に於いては、構造最適化の評価関数として、固有値の大きい方から選択された複数の固有値の線形和が採用されてよい。具体的な装置の構成に於いては、固有値の大きい方から複数個の組の固有ベクトルを選択した後、選択された組の固有ベクトルの各々に対して変位ベクトルの算出及び感度の算出が実行され、構造体全体の固有値の感度として、各固有値の感度の線形和が算出され、その構造体全体の固有値の感度を用いて、選択された固有値の線形和が評価関数として最小化されるように設計変数の更新が実行されてよい。なお、複数の固有値及びそれらの感度の線形和は、重み付き和であってよく、重み付き和に於ける各項の重みは、上記の如く、個々の固有値が大きいほど重み付き和に於ける対応する固有値に対する重みが大きくなるよう設定されてよい。このように重みを設定すると、より大きな固有値、即ち、より大きなコンプライアンスを与える荷重ベクトルに対する剛性を高くして、かかる荷重ベクトルに対する変形を最小化することが可能となる。
【００２３】
かくして、上記の本発明の装置の固有値感度算出手段に於いて、固有値を二つ以上用いるときには、固有値の感度が、局所コンプライアンス行列の複数の固有値のうちの値の大きい方から選択された複数の固有値の設計変数に対する感度の線形和であってよく、設計変数更新手段に於いて、固有値を二つ以上用いるときには、局所コンプライアンス行列の複数の固有値のうちの値の大きい方から選択された複数の固有値の線形和が最小化されるよう設計変数が更新されるようになっていてよい。また、上記の固有値及びその感度の線形和は、重み付き和であってよく、その場合の各項の重みは、個々の固有値が大きいほど重み付き和に於ける対応する固有値に対する重みが大きくなるよう設定されていてよい。上記の如く値の大きい複数の固有値が全体的に最小化されるように設計変数の更新を実行する場合には、構造体の構造がより安定的に最適な構造へ収束されることとなる。なお、設定された条件によって、最大固有値を与える固有値の交代がない場合或いはそのように想定される場合には、評価関数を最大固有値とし、設計変数に対する感度は、最大固有値の感度であってもよく、そのような場合も本発明の範囲に属することは理解されるべきである。（その場合には、局所コンプライアンス行列の固有値解析に於いて、最大固有値とそれに対応する固有ベクトルが選択される。）
【発明の効果】
【００２４】
上記の本発明によれば、支持構造又は支点のない構造体の場合、即ち、構造体の剛性行列に逆行列が存在しない場合でも構造最適化が可能な手法を提供される。かかる本発明の装置によれば、構造体に作用する荷重が不定である場合だけでなく、構造体を支持する構造が不定である場合でも、任意の荷重が付与される場合に構造体の剛性が最大化される最適な構造を見出すことが可能となると期待される。また、その際、支持構造のある構造体の場合と異なり、局所コンプライアンス行列の最大固有値を与える固有ベクトルが設計変数の更新に伴う構造体の構造の変化に伴って交代することがあり、かかる固有ベクトルの交代があると、構造体の構造の収束が困難となることがあるところ、本発明の装置は、局所コンプライアンス行列の複数の固有値が全体的に最小化されるように評価関数の感度及び設計変数の更新を実行するようになっていてよく、これにより、最大固有値を与える固有ベクトルの交代があっても、構造体の構造を収束させて、最適な構造が得られるようになっていてよい。
【００２５】
本発明のその他の目的及び利点は、以下の本発明の好ましい実施形態の説明により明らかになるであろう。
【図面の簡単な説明】
【００２６】
【図１】図１は、本発明の装置による構造最適化演算が実行されるコンピュータを模式的に表した図である。
【図２】図２（Ａ）は、構造体における弾性システム（二次元）を模式的に表した図であり、図２（Ｂ）は、図２（Ａ）のシステムに於いて構造体に支持構造がない場合に構造体に於いて作用する荷重の態様のいくつかの例を示したものである。
【図３】図３は、支持構造がない直方体の構造体に於いて底面（ｙ＝０の面）にのみ不定荷重が作用すると仮定したモデルの模式図である。
【図４】図４（Ａ）は、図３のモデルに於いて算出される局所コンプライアンス行列の複数の固有ベクトルのそれぞれから得られるベクトル成分の分布の例を示している。（ａ）〜（ｆ）は、固有値が大きい固有ベクトルの順に並べられている。図４（Ｂ）は、図３のモデルから構造最適化演算を実施した後の図４（Ａ）のそれぞれに対応するベクトル成分の分布の例を示している。
【図５】図５は、構造最適化演算による設計変数の更新前後の構造体の局所コンプライアンス行列の複数の固有ベクトルφｉとそれらに対応する固有値λｉの大きさの関係の変化を説明する図である。図５（Ａ）は、支持構造のある構造体の場合であり、図５（Ｂ）は、支持構造のない構造体の場合である。
【図６】図６は、本発明の装置による構造最適化演算の処理をフローチャートの形式にて表した図である。
【図７】図７は、本発明の装置による支持構造のない構造体の構造最適化演算によって得られた構造体の形状の例である。図示の例は、設計変数の更新を１００回実行した状態のものである。
【符号の説明】
【００２７】
１…コンピュータ本体
２…コンピュータ端末
３…モニター
４…キーボード、マウス（入力装置）
【発明を実施するための最良の形態】
【００２８】
コンピュータ装置の構成
本発明による構造最適化演算を実行する装置は、この分野で通常使われている形式の、図１に例示されている如き、コンピュータ装置１に於けるコンピュータ・プログラムの作動により実現されてよい。コンピュータ装置１には、通常の態様にて、双方向コモン・バスにより相互に連結されたＣＰＵ、記憶装置、入出力装置（Ｉ／Ｏ）が装備され、記憶装置は、本発明の演算で使用する演算処理を実行する各プログラムを記憶したメモリと、演算中に使用されるワークメモリ及びデータメモリを含んでいる。また、実施者によるコンピュータ装置１への指示及び計算結果その他の情報の表示及び出力は、コンピュータ装置１に接続されたコンピュータ端末装置２を通じて為される。コンピュータ端末装置２には、通常の態様にて、モニター３とキーボード並びにマウスといった入力装置４が設けられ、プログラムが起動されると、実施者は、プログラムの手順に従ってモニター３上の表示に従って、入力装置４を用いてコンピュータ装置１に各種の指示及び入力を行うとともに、モニター３上にてコンピュータ装置１からの演算状態及び演算結果等を視覚的に確認することが可能となる。なお、図示していないその他の周辺機器（結果を出力するプリンタ、計算条件及び演算結果情報等を入出力するための記憶装置など）がコンピュータ装置１及びコンピュータ端末装置２に装備されていてよいことは理解されるべきである。コンピュータ装置１を用いて、以下に述べる各種の処理又は演算を実行する際には、通常の態様にて、各種の処理又は演算に必要なプログラムが起動され、実施者は、コンピュータ端末装置２に於いて、プログラムに於いて準備された入力手順に従って、演算に必要なデータ、演算時の計算条件その他の各種設定を入力し、演算が開始される。そして、演算の実行中又は終了後に、演算結果が、適宜、コンピュータ端末装置２を通じて出力可能となる。
【００２９】
上記のコンピュータ装置１に於いて、本発明による構造最適化演算を行う場合、端的に述べれば、まず、入力装置４からの実施者の入力操作によって、最適化された構造体が生成される設計領域の寸法、構造体を形成する材料の特性値（ヤング率、ポアソン比）及び体積制約割合を含む設計条件や構造体に作用する荷重条件を含む境界条件などの構造最適化演算を実行するために必要な初期条件の設定が実行される。そうすると、コンピュータ装置１は、構造最適化演算処理の準備として、プログラムメモリに記憶されたプログラムを用いて、設計領域を、典型的には格子状に、複数の要素に分割する処理を実行すると共に、得られた要素の各々への設計変数の初期値の割り当てを実行する。かくして、かかる演算処理の準備が完了すると、コンピュータ装置１は、プログラムメモリに記憶されたプログラムに従い、（ａ）構造体の剛性行列Ｋ、擬似逆行列Ｋ＋、局所コンプライアンス行列Ｃの算出、（ｂ）局所コンプライアンス行列Ｃの固有値解析、（ｃ）各要素の各節点に生ずる変位の算出、（ｄ）設計領域内にて形成される構造体の評価関数の感度（各要素の設計変数の変化に対する評価関数の変化量）の算出、（ｅ）評価関数が最適化されるように各要素の設計変数を更新する処理を含む一連の演算処理を、後述の終了条件が成立するまで反復して実行し、これにより、設計変数により形成される構造体の構造又は形状が、逐次的に最適化された状態に近づいていくことが期待される。また、一連の演算処理の結果として得られた構造体の構造又は形状、設計変数の分布、局所コンプライアンス行列Ｃの固有値の推移、固有ベクトルの成分の分布などのグラフィック表示がモニター３上に表示されてよい。
【００３０】
ロバスト構造最適化演算の概要
本発明の装置の実施形態に於ける構造最適化演算では、基本的には、ロバスト構造最適化演算の手法（非特許文献３、４）に従って、不定の荷重の存在下に於いて、剛性を最大にする、即ち、変形を最小にする構造体の構造を最適な構造として求めるために、構造体の変形程度を表す構造体のコンプライアンスを最大にする荷重（最悪荷重）の存在下に於いて、構造体のコンプライアンスが最小化されるように、設計領域内の構造体の体積が所与の体積制約割合の条件を満たすことを条件として、構造体の構造を決定する設計変数の更新する演算処理が実行され、更新された設計変数により表される構造体の構造が目的の構造体の構造として見出されるようになっていてよい。
【００３１】
より具体的には、まず、ロバスト構造最適化演算に於ける最小化の対象となる評価関数となるコンプライアンスｌは、
ｌ＝ＦＴ・Ｕ      …（５）
により与えられる。ここに於いて、ＦとＵは、構造体に於いて格子状の要素を設定した場合の要素の節点に於ける荷重を成分とする荷重ベクトル（上付きの添え字Ｔは、転置を表す。）と変位を成分とする変位ベクトルであり、既に述べた如く、構造体に於ける剛性行列Ｋを用いた式（２）［ＫＵ＝Ｆ］にて互いに関係付けられる。上記の式（５）のコンプライアンスｌに於いて、荷重ベクトルＦが０以外の値を持つのは、構造体に於いて荷重の作用する「荷重作用領域」内の節点のみとなるので、構造体のコンプライアンスｌは、荷重作用領域内の荷重分布を表す荷重ベクトルＦｌと荷重作用領域の変位分布を表す変位ベクトルＵｌを用いて
ｌ＝ＦｌＴ・Ｕｌ      …（６）
により表される。そして、荷重ベクトルＦｌと変位ベクトルＵｌは、式（１）により局所コンプライアンス行列Ｃを用いて互いに関係付けられるので、コンプライアンスｌは、
ｌ＝ＦｌＴ・Ｃ・Ｆｌ      …（７）
により与えられる。このコンプライアンスｌに関し、その値が最大となるとき、即ち、構造体に任意の荷重が作用する場合にその変形が最大となるときは、荷重ベクトルＦｌが局所コンプライアンス行列Ｃの最大固有値を与える固有ベクトルとなっているときであり、荷重ベクトルＦｌのノルム||Fl||＝１としたとき、コンプライアンスｌは、局所コンプライアンス行列Ｃの最大固有値となることが、Ｒａｙｌｅｉｇｈ−Ｒｉｔｚの定理により理論的に示されている（非特許文献３、４参照）。従って、不定の荷重に対して構造体の剛性を最大化する構造最適化演算に於いては、荷重ベクトルが局所コンプライアンス行列Ｃの最大固有値を与える固有ベクトルとなっているときのコンプライアンスｌを最小化するように構造体の構造を決定する設計変数が更新されることとなる。
【００３２】
ところで、式（１）に於ける局所コンプライアンス行列Ｃに関して、荷重作用領域内の荷重ベクトルＦｌと構造体全体の荷重ベクトルＦとを関連づける変換行列Ｈを
Ｆ＝Ｈ・Ｆｌ      …（８ａ）
Ｆ＝ＨＴ・Ｆ      …（８ｂ）
と置くと、荷重作用領域内の変位ベクトルＵｌは、変換行列Ｈと構造体の剛性行列Ｋの逆行列Ｋ−１を用いて、
Ｕｌ＝ＨＴ・Ｕ＝ＨＴ・Ｋ−１・Ｆ＝ＨＴ・Ｋ−１・Ｈ・Ｆｌ      …（９）
と表されるので、式（１）から、局所コンプライアンス行列Ｃは、
Ｃ＝ＨＴ・Ｋ−１・Ｈ      …（１０）
により与えられる。なお、変換行列Ｈは、設計領域に対して指定される荷重作用領域の範囲に基づいて決定される。
【００３３】
従って、ロバスト構造最適化演算では、設定された初期条件を用いて式（２）の剛性行列Ｋを決定した後、その逆行列Ｋ−１を用いて式（１０）により局所コンプライアンス行列Ｃを算出し、局所コンプライアンス行列Ｃの固有値解析を実行し、そこで得られた最大固有値を与える固有ベクトルを荷重ベクトルＦとして用いて
Ｕ＝Ｋ−１・Ｆ      …（１１）
により、式（２）を解いて構造体の変位ベクトルＵを算出し、しかる後に、算出された変位ベクトルＵを用いて、コンプライアンスｌを最小化するように設計変数の更新が実行される。
【００３４】
支持構造のない構造体の構造最適化演算
（ａ）支持構造のない構造体に於ける変位ベクトルの算出-擬似逆行列Ｋ＋を用いた演算
「発明の概要」の欄に於いて既に述べた如く、現実の任意の構造体に於いては、ユーザーの手からの荷重を受ける携帯電話の外殻など、構造体を支持し或いは固定する構造や支点の位置が不定な場合がある。従って、構造最適化演算に於いては、荷重だけでなく、支持構造も不定である場合、即ち、支持構造のない状態の構造体についても、構造最適化演算が実行できることが好ましい。
【００３５】
しかしながら、構造体に於いて支持構造がない場合、構造体は自由な平行移動及び回転（剛体運動又は剛体変位）が許されるので、或る荷重ベクトルＦが与えられた場合の式（２）［ＫＵ＝Ｆ］の解である変位ベクトルＵは、一意に定まらなくないこととなる。即ち、支持構造がない構造体に於いて、数学的には、剛性行列Ｋの逆行列Ｋ−１が存在しないため、上記のロバスト構造最適化演算の実行に際して、逆行列Ｋ−１を用いる式（１０）による局所コンプライアンス行列Ｃの算出及び式（１１）による変位ベクトルＵの算出を実行することができないこととなっている。
【００３６】
ところで、式（２）［ＫＵ＝Ｆ］の如き連立一次方程式に於いて、方程式の係数を成分とする行列Ｋに逆行列Ｋ−１が存在しない場合に擬似逆行列Ｋ＋を用いることにより方程式の解を得ることが可能である。擬似逆行列Ｋ＋の性質から、式（２）［ＫＵ＝Ｆ］が解を有するとき、解Ｕは、任意のベクトルＵｗを用いて、
Ｕ＝Ｋ＋・Ｆ＋（Ｉ−Ｋ＋・Ｋ）Ｕｗ      …（１２）
により表される（Ｉは、単位行列である。）。ここで、構造体の剛性方程式である式（２）［ＫＵ＝Ｆ］が静的な状態を表し、荷重ベクトルＦによって平行及び回転について加速運動が生じていないとすると、式（１２）の右辺の第２項（Ｉ−Ｋ＋・Ｋ）・Ｕｗは、剛性行列Ｋを乗算すると、擬似逆行列Ｋ＋の定義から
Ｋ・（Ｉ−Ｋ＋・Ｋ）・Ｕｗ＝０      …（１３）
となり、構造体に対して力が作用していない状態でも生じ得る構造体の剛体運動を表す変位ベクトルであると考えられる。一方、式（１２）の右辺の第１項Ｋ＋・Ｆは、この場合、剛体運動以外の構造体内に生ずる変位であると考えられ、従って、荷重ベクトルＦは、剛体運動を生じさせない構造体内に於いて対称的に変位を生ずる荷重ベクトルとなっていると考えられる（もし対照的でない変位があったとすると、支持構造がないので、剛体運動となっているはずである。）。即ち、荷重ベクトルＦが、図２に模式的に例示されている如く、構造体に対して荷重の合計とモーメントの合計とがそれぞれ０となるように作用している荷重ベクトルであるときには、式（１２）の右辺の第１項Ｋ＋・Ｆが構造体内部の対称的な変位を表していることとなる。従って、式（９）に於いて、逆行列Ｋ−１に代えて擬似逆行列Ｋ＋を用いた場合、即ち、変位ベクトルＵｌが
Ｕｌ＝ＨＴ・Ｋ＋・Ｆ＝ＨＴ・Ｋ＋・Ｈ・Ｆｌ      …（１４）
により表された場合、荷重ベクトルＦが剛体運動を生じない荷重ベクトルであるとき、Ｋ＋・Ｆが剛体運動による変位を含まないので、変位ベクトルＵｌは、構造体内部の対称的な変位となっていると考えられる。
【００３７】
更に、上記の式（１４）を参照して、式（１０）に代えて、局所コンプライアンス行列Ｃを、擬似逆行列Ｋ＋を用いて
Ｃ＝ＨＴ・Ｋ＋・Ｈ      …（１５）
により決定した場合に、或る荷重ベクトルＦｌｉ（ただし、||Flｉ||＝１）が局所コンプライアンス行列Ｃの一つの固有値λｉを与える固有ベクトルφｉであるとき、即ち、
Ｕｌｉ＝Ｃ・Ｆｌｉ＝λｉＦｌｉ（Ｃ・φｉ＝λｉφｉ）      …（１６）
であるとき、荷重ベクトルＦｌｉは、構造体内部の対称的な変位ベクトルＵｌｉに平行となるので、固有ベクトルφｉは、構造体に対して荷重の合計とモーメントの合計とがそれぞれ０となるように作用している荷重ベクトルを与えることとなる。実際、図３に模式的に描かれている如き、直方体に於いて、その底面（ｙ＝０）を荷重作用領域として設定して（ｙ＝０の面に不定荷重が作用すると設定したことに相当する。）、式（１５）による局所コンプライアンス行列Ｃの固有ベクトルを算出すると、図４（Ａ）の（ａ）〜（ｆ）に描かれている如く、固有ベクトルの成分により表されるｙ＝０の面内に於ける荷重分布は、荷重の合計とモーメントの合計とがそれぞれ０となる滑らかな分布となることが見出された。即ち、局所コンプライアンス行列Ｃの固有ベクトルを荷重ベクトルとした場合、荷重ベクトルと変位ベクトルとは平行となるので、構造体には対称的な変位が生ずることとなる。
【００３８】
かくして、本実施形態に於ける支持構造のない構造体の構造最適化演算では、剛性行列Ｋの逆行列Ｋ−１ではなく、擬似逆行列Ｋ＋を用いることにより、局所コンプライアンス行列Ｃは、式（１５）により算出され、局所コンプライアンス行列Ｃの固有ベクトルは、構造体内に於いて対称的な分布の荷重ベクトルＦを与え、変位ベクトルＵが、
Ｕ＝Ｋ＋・Ｆ      …（４）
により算出される（Ｕは、構造体内に対称的な変位分布を与えるベクトルとなる。）。
【００３９】
（ｂ）支持構造のない構造体に於けるコンプライアンスの変化
局所コンプライアンス行列Ｃに対しては複数組の固有値と固有ベクトルが得られるところ、支持構造のある構造体の場合、図５（Ａ）に模式的に描かれている如く、構造最適化演算の過程に於いて、最初の段階から、一つの固有ベクトル（φ０）が、徐々に変化しながら、最大のコンプライアンス、即ち、最大固有値（λ０）を与えるようになっている。構造体が支持構造を持つ場合、固有ベクトルにより与えられる荷重分布が対称的となる条件がなく、固有ベクトルに一方向に向かう傾向のあるベクトル（φ０）が存在し、そのような固有ベクトル（φ０）が最大固有値（λ０）を与える荷重ベクトルとなると、支点などの支持構造とその付近に大きな変形を生み、コンプライアンスを著しく増大させ、かくして、構造最適化演算の反復が進んでも、常に最大固有値を与える固有ベクトルであり続けることとなる。即ち、構造最適化演算の過程に於いて、複数存在する固有ベクトルのうち、最大固有値を与える固有ベクトルは通常交代しないので、最大固有値を与える荷重分布及び変位分布は徐々に変化することとなる。従って、従前のロバスト構造最適化演算では、構造体のコンプライアンスを最小化する設計変数の更新に於いて使用する変位ベクトルＵの算出のために、常に、局所コンプライアンス行列Ｃの固有値解析により得られた固有値のうちの最大のものを与える固有ベクトルを荷重ベクトルとして選択していれば、構造体の構造は、安定的に最適な構造へ収束することとなっていた。
【００４０】
一方、支持構造のない構造体の場合、（擬似逆行列Ｋ＋を用いて得られた）局所コンプライアンス行列Ｃの固有ベクトルに於いては、図２或いは図４の例示からも理解される如く、複数得られる固有ベクトルは、それぞれの態様にて、構造体に於いて荷重の合計及びモーメントの合計が共に０となる対称的な荷重分布を保存しながら、構造体の構造最適化演算の過程に於いて徐々に変化することとなり、いずれも、一方向に向かうベクトルになる傾向がなく、変形が局所的に一方向に生ずることがない（図５（Ｂ）に示されているように、固有値λ０を与える一方向に向かう固有ベクトルφ０が存在しない。）。こういった支持構造のない構造体に於ける性質により、構造最適化演算の過程では、図５（Ｂ）に描かれているように、最大固有値を与える固有ベクトルが分布の態様の異なるベクトルへ代わってしまう現象、即ち、最大固有値を与える固有ベクトルの交代（最大固有値がλ１からλ２へ代わるため、最大固有値を与える固有ベクトルがφ１からφ２に交代する。）が生じ得ることが見出された。かかる最大固有値を与える固有ベクトルの交代が生ずると、最大固有値を与える固有ベクトルを用いて算出される変位ベクトルの分布が不連続的に変化し、これにより、構造体のコンプライアンスを最小化するべく更新される設計変数も不連続的に変化することとなり、構造体の構造が安定的に最適な構造へ収束し難くなる場合が生じた。
【００４１】
そこで、本実施形態では、局所コンプライアンス行列Ｃから得られる複数の固有ベクトルにより与えられる複数の固有値が全体的に最小化されるように設計変数の更新が実行されてよい。具体的には、局所コンプライアンス行列Ｃの固有値解析を実行した後、まず、固有値の大きい方から所定の数Ｎの固有値λｉと固有ベクトルφｉの組が選択され、
λｔ＝ｗ１λ１＋ｗ２λ２＋…＋ｗＮλＮ      …（１７）
により与えられる固有値の線形和（重み付き和）λｔが評価関数として最小化されるように設計変数の更新が実行されてよい。ここにおいて、ｗｉは、線形和λｔに於ける各項λｉについての重みである。かかる重みｗｉは、具体的には、例えば、下記のように設定されてよい。
ｗｉ＝（λｉ−λＮ＋１）μ／Σ（λｉ−λＮ＋１）μ      …（１８）
（Σは、ｉ＝１〜Ｎまでの総和であり、μはべき指数である。）
式（１８）の重みｗｉによれば、固有値λｉが大きいほど、線形和λｔに於ける寄与が大きくなり、最小化に於ける寄与が増大することとなるので、構造体の構造がより速やかに最適な構造へ収束することが期待される。かくして、設計変数に対する最小化される評価関数の感度は、各固有値の感度∂λｉ／∂ρｅを用いて（ρｅは要素ｅに割り当てられた設計変数）
ｓｅ＝−Σｗｉ（∂λｉ／∂ρｅ）      …（１９）
（Σは、ｉ＝１〜Ｎまでの総和）
により与えられる。
【００４２】
構造最適化演算手法の概要（密度法を用いたトポロジー最適化演算の概要）
本発明の装置の実施形態に於ける具体的な構造最適化演算は、上記の如く、擬似逆行列を利用し、評価関数として複数の固有値の線形和を採用しながら、任意のトポロジー構造最適化や形状最適化などの任意のアルゴリズムに従って実行されてよいことは理解されるべきである。本実施形態に於いては、例えば、以下に概説される密度法（Density-based  Approach）を用いたトポロジー最適化演算が採用されてよい。
【００４３】
密度法を用いたトポロジー最適化演算に於いては、具体的には、まず、三次元空間に於いて任意の寸法の設計領域（例えば、図３の如き直方体の領域など）が設定され、設計領域を格子状に複数の要素に分割し、各要素に対する設計変数として各要素内材料の密度（材料密度）ρｅ（ｅは、要素番号）が
０≦ρｅ≦１      …（２０）
の条件を満たすよう設定され、材料密度ρｎのときの要素の材料のヤング率Ｅｎが、
Ｅｅ（ρｅ）＝ρｅＰＥｏ      …（２１）
と設定される。ここで、Ｐは、ペナルティパラメータと称する値（＞１）であり、Ｅｏは、材料密度ρｅ＝１のときの材料のヤング率である。なお、後に説明される変位ベクトルＵの演算に於ける特異性を回避するために、ヤング率Ｅｎは、
Ｅｅ（ρｅ）＝Ｅmin＋ρｅＰ(Ｅｏ−Ｅmin)      …（２２）
と修正されて用いられてよい（非特許文献１、２）。ここで、Ｅminは、材料密度ρｅ＝０のときの仮想的な材料のヤング率である。そして、評価関数として構造体のコンプライアンスが採用され、有意な材料密度ρｅを有する要素から構成される構造体のコンプライアンスが最小化されるように要素の材料密度ρｅを反復して更新する演算処理が実行され、設計領域内に於ける有意な材料密度を有する要素の分布の形状が、目的の構造体の形状、即ち、剛性ができるだけ高くなる構造体の形状、として見出される。本実施形態に於いては、最適化の評価関数である構造体のコンプライアンスは、構造体の局所コンプライアンス行列Ｃの固有値又は固有値の線形和（式（１７））となる。
【００４４】
なお、材料密度ρｎの更新は、設計領域の体積に対する構造体の体積の割合Ｖflacが所与の体積制約割合の条件を満たすように実行される。設計領域内の構造体の体積は、要素ｅの体積ｖｅを用いて、Σρｅ・ｖｅ（ここで、Σは、全要素についての総和である。）により与えられるので、材料密度ρｅは、体積制約割合Ｖflacの条件として、
Σρｅ・ｖｅ≦Ｖflac・Ｖｔ      …（２３）
（ここで、Ｖｔは、設計領域の体積）
が満たされるように更新されることとなる。また、演算結果に於いて、各要素の材料密度の値ができるだけ０又は１に収束するように材料密度ρｅには、任意のフィルター処理を施した値が用いられてよい。
【００４５】
トポロジー最適化演算の具体的な処理手順に於いては、まず、初期状態の設定として、設計領域内にて生成される構造体に対する拘束条件、構造体に対して作用する荷重条件、体積制約割合（設計領域の体積に対する生成されるべき構造体の体積の割合）、材料密度ρｅの初期値、構造体を構成する材料の特性値（ヤング率、ポアソン比）を含む構造最適化演算のための初期条件が実施者により設定される。そして、設定された設計領域に対して、任意のプログラムを用いて要素の設定が為される。材料密度ρｅの初期値は、設計領域内全域に於いて均一な値であってよい。また、本実施形態に於いては、最適化の対象は、支持構造のない構造体であるので、構造体に於ける自由な剛体運動を制限する拘束条件は設定されない。更に、本実施形態は、ロバスト構造最適化の手法に従って荷重条件を設定するので、構造体に作用する荷重条件については、荷重作用領域は設定されるが、かかる領域内での荷重分布は局所コンプライアンス行列Ｃの固有値解析により得られた固有ベクトルによって決定されるので、荷重の具体的な値及びその分布は初期条件としては設定されない。
【００４６】
次いで、下記の（ａ）、（ｂ）、（ｃ）の演算処理が所定の終了条件が成立するまで反復して実行される。
【００４７】
（ａ）荷重条件の設定−構造体に最悪荷重が作用している条件を見出すために、まず、設計領域内の各要素の材料密度ρｅから各要素のヤング率Ｅｅ（ρｅ）を計算し、かかるヤング率を用いて設計領域内に形成されている構造体の剛性行列Ｋが決定される。剛性行列Ｋに於いて使用される各要素のヤング率Ｅｅには、式（２１）又は（２２）が用いられる。そして、剛性行列Ｋを用いて、擬似逆行列Ｋ＋、局所コンプライアンス行列Ｃが算出され、局所コンプライアンス行列Ｃの固有値解析を実行し、複数組の固有値λｉ及び固有ベクトルφｉのうち、固有値の大きい方から所定の数Ｎの組が選択され、選択された固有ベクトルφｉがそれぞれ荷重ベクトルＦｉに設定される。なお、擬似逆行列Ｋ＋及び局所コンプライアンス行列Ｃは、任意のプラグラムに従って計算されてよい（プログラム言語に用意されている関数を利用してもよく、任意のアルゴリズムによって算出されてもよい。）。また、局所コンプライアンス行列Ｃの固有値解析も任意のプラグラムに従って計算されてよい（プログラム言語に用意されている関数を利用してもよく、任意のアルゴリズムによって算出されてもよい。）
【００４８】
（ｂ）変位ベクトルの算出−上記の荷重条件の下、設計領域内に於ける各要素の各節点の変位が算出される。具体的な算出方法は、有限要素法等に従った任意のプログラムを用いたものであってよい。既に述べた如く、本実施形態に於いては、最適化の対象となる構造体に支持構造がなく、構造体の全要素の剛性行列Ｋの逆行列Ｋ−１が存在しないので、各要素の各節点の変位は、剛性行列Ｋの擬似逆行列Ｋ＋を用いて、式（４）によって算出される。また、最適化の評価関数として、局所コンプライアンス行列Ｃの複数の選択された固有値λｉの線形和が用いられる場合には、各節点の変位は、複数の選択された固有ベクトルの各々を荷重ベクトルとして用いた場合のそれぞれについて算出される（即ち、評価関数がＮ個の固有値の線形和（式（１７））であるときには、Ｎ個の対応する固有ベクトルを荷重ベクトルとした場合に対応して、Ｎ個の変位ベクトルが算出されることとなる。）。
【００４９】
（ｃ）材料密度の更新−上記の演算により得られた設計領域内に全要素の節点変位ベクトルＵと全要素の剛性行列Ｋとを用いて、式（１９）により構造体の材料密度ρｅに対する評価関数の感度（材料密度ρｅについての評価関数の微分）が算出され、更に、材料密度ρｅに対する構造体の評価関数の感度を用いて、体積制約割合条件（式（２３））が満たされることを拘束条件として構造体の評価関数が最小化されるように新たな材料密度ρｅを算出する演算処理、即ち、材料密度ρｎの更新処理が実行される。材料密度ρｅの更新処理は、例えば、最適性基準法、逐次線形計画法、逐次凸関数近似法などの手法に従って実行されてよい。具体的には、材料密度ρｅについての各固有値の感度∂λｉ／∂ρｅは、式（２２）を参照して、
∂λｉ／∂ρｅ＝−ＵｉＴ[pρｅp-1(Ｅｏ−Ｅmin)Ｋ０ｅ]Ｕｉ      …（２４）
により与えられる。ここで、Ｕｉは、固有値λｉを与える固有ベクトルφｉを荷重ベクトルとした場合の変位ベクトルであり、Ｋ０ｅは、要素ｅの単位剛性行列（ヤング率＝１と設定した剛性行列）である。従って、固有値λ１〜λＮを選択したときには、式（２４）、式（１８）及び式（１９）を用いて、各要素ｅの設計変数に対する評価関数の感度ｓｅが算出される。そして、例えば、新しい材料密度ρｅ(k+1)は、それまでの材料密度ρｅ(k)を用いて（kは、反復処理の回数）、
ρｅ(k+1)＝Ｂｅ(k)η・ρｅ(k)      …（２５ａ）
或いは、
ρｅ(k+1)＝max(ρｅ(k)-m,  min(ρｅ(k)+m,  Ｂｅ(k) η・ρｅ(k)))…（２５ｂ）
の形式にて与えられる。ここで、Ｂｎ(k)は、感度ｓｅを用いて、
Ｂｅ(k)＝Ｈｂ(ｓｅ)[Λ（∂ｖ（Ρ(k)）／∂ρｅ）]−１      …（２６）
により与えられる。ここで、Ｈｂ(ｓｅ)は、ρｅを滑らかにするフィルター関数であり、∂ｖ（Ρ(k)）／∂ρｅは、式（２３）の体積制約条件
ｖ（Ρ(k)）＝Σρｅ(k)・ｖｅ−Ｖflac・Ｖｔ≦０
のρｅについての微分であり（∂ｖ（Ρ(k)）／∂ρｅ＝ｖｅ）、Λは、ラグランジュの未定係数であり、二分法などによって体積制約条件（式（２３））を満たす値が決定されてよい。ｍは、1回の更新に於ける材料密度の変化限度であり、ηは、Ｂｅ(k)のべき指数であり、＜１である。なお、Ｂｅ(k)ηの演算負荷を低減するために、
Ｂｅ(k)η〜１＋η（Ｂｅ(k)−１）      …（２７）
による近似演算が用いられてよい。具体的な演算処理のアルゴリズムは、例えば、非特許文献１、２などに記載されているものが使用されてよい。
【００５０】
上記の（ａ）〜（ｃ）の演算処理は、（ｃ）の処理にて新たな材料密度ρｎ(k+1)が得られる度に、それを材料密度ρｎ(k)として用いて、繰り返し実行される。かかる反復処理を実行されるにつれて、材料密度ρｎの分布は、体積制約割合条件（２３）を満たしつつ、構造体の評価関数を最小化する方向、即ち、構造体の剛性が高くなる方向に変化していくことが期待される。そして、かかる処理が反復されていく過程で、好ましくは、各要素の材料密度ρｎの値は、０又は１に収束され、これにより、設計領域内に於いて有意な材料密度（ρｎ≒１）を有する要素の分布の形状が体積制約割合を満たす条件下で剛性が高くなった構造体の形状として見出されることとなる。なお、（ａ）〜（ｃ）の反復処理を終了するための所定の終了条件としては、（ｉ）反復処理回数が任意に設定されてよい回数に到達した場合、（ii）構造体の評価関数の値が収束した場合（処理サイクル毎の評価関数の変化量が所定の小数以下となった場合）、（iii）更新される材料密度ρｎの値が収束した場合（処理サイクル毎の材料密度の変化量の最大値が所定の小数以下となった場合）などが選択されてよい。（ii）又は（iii）を終了条件として採用する場合には、結果として得られる構造体の形状は、体積制約割合条件の下で、最適化されていることが期待される。
【００５１】
装置の作動
図６を参照して、本実施形態の装置に於ける構造最適化演算の処理手順は、次のように実行されてよい。ロバスト最適化演算を実行するために必要な初期条件の設定が実行される（ステップ１）。既に触れた如く、これらの初期条件の設定は、演算の実施者により入力装置４を通じて実行される。しかる後、コンピュータの作動により、まず、剛性行列Ｋの演算（ステップ２）、擬似逆行列Ｋ＋の演算（ステップ３）、擬似逆行列Ｋ＋を用いた局所コンプライアンス行列Ｃの演算（ステップ４：式（１５））が順に実行され、局所コンプライアンス行列Ｃの固有値解析が実行されて、局所コンプライアンス行列Ｃの所定数の組の固有値λｉと固有ベクトルφｉが選択される（ステップ５）。局所コンプライアンス行列Ｃの固有値と固有ベクトルが選択されると、選択された固有ベクトルφｉをそれぞれ荷重ベクトルとして、疑似逆行列Ｋ＋を用いて、有限要素法による設計領域内の各要素の各節点の変位（変位ベクトルＵｉ）が演算される一方（ステップ６：式（４））、選択された固有値λｉを用いて、評価関数感度に於ける各固有値の重みｗｉが演算される（ステップ７：式（１８））。そして、固有ベクトル毎の変位ベクトルＵｉと剛性行列Ｋの成分を用いて式（２４）により、各固有値λｉの材料変数に対する感度（∂λｉ／∂ρｅ）が算出され、それらの感度と重みｗｉを用いて、各要素ｅの評価関数（固有値の線形和）の感度ｓｅが算出される（ステップ８：式（１９）、固有値感度算出手段）。そうすると、構造体の評価関数の感度ｓｅを用いて、上記の態様にて、トポロジー最適化の手法に従った材料密度（密度関数）の更新演算が実行され、設計領域内の新たな材料密度の分布が得られることとなる（ステップ９：式（２５ａ）又は式（２５ｂ））。かくして、一連の演算処理が完了すると、設計領域内の材料密度の分布、有意な材料密度を有する要素の分布、又は、構造体の形状がコンピュータのモニター３上に表示されてよい（ステップ１０）。
【００５２】
上記の一連の処理が実行されると、既に述べた如き所定の終了条件が成立しているか否かが判定され（ステップ１１）、所定の終了条件が成立していなければ、ステップ２〜ステップ１０の実行が繰り返され、所定の終了条件が成立すると、演算処理が終了する。
【００５３】
計算実験例
上記に説明した本発明の装置の有効性を検証するために、図３に例示する直方体を設計領域１０とし、荷重作用領域として、直方体の底面（ｙ＝０の面）に上下方向に不定荷重が作用する領域とし、上記の処理手順に従って、計算実験を実行し、構造体が最適化されることを確認した。なお、設計領域の寸法は、ｘ×ｙ×ｚ＝１４×７×１８と設定し、一つの要素の体積は、ｖｅ＝１と設定し、体積制約割合Ｖflacは、０．３と設定した。また、構造体を構成する材料のヤング率は、Ｅｏ＝１、Ｅmin＝１０−９とし、ポアソン比は、０．３とし、ペナルティパラメータは、Ｐ＝３とした。また、式（１７）の評価関数λｔに於いて選択する固有値と固有ベクトルの組数は、Ｎ＝４とし、式（１８）に於いて、μ＝１とした。更に、式（２５ａ）、（２５ｂ）に於いて、η＝０．１５、ｍ＝０．０５と設定した。
【００５４】
まず、図４（Ａ）は、既に説明されている如く、構造最適化演算開始前の図３に示す直方体に於ける局所コンプライアンス行列Ｃの固有値解析によって得られた固有ベクトルの成分のｙ＝０の面内に於ける分布を示している（初期条件に於いて、ｙ＝０の面にのみ荷重が生ずるように設定しているので、ｙ＝０以外の領域の荷重は全て０である。）。図４（Ａ）の（ａ）〜（ｆ）の分布は、値が大きい方から順に六つの固有値を与える固有ベクトルの分布である。上記の如く、構造最適化演算開始前の段階で、支持構造のない構造体に於いて算出される局所コンプライアンス行列Ｃの固有ベクトルは、構造体に於いて一方向の運動を発生させない対称的な荷重分布、即ち、荷重の合計とモーメントの合計との両方が０となる荷重分布を与えることが示された。なお、構造体の荷重作用領域に於ける変位ベクトルは、局所コンプライアンス行列Ｃの固有ベクトルと平行であるので、構造体内に生ずる変位も対称的になることが理解される。
【００５５】
一方、図４（Ｂ）は、図３に示す直方体に於いて本実施形態による構造最適化演算に於いて設計変数の更新を１００回実行した後の構造体に於ける局所コンプライアンス行列Ｃの固有値解析によって得られた、値が大きい方から六つの固有値を与える固有ベクトルを示している。図４（Ｂ）に於ける（ａ）〜（ｆ）は、図４（Ａ）の（ａ）〜（ｆ）の固有ベクトルに対応する固有ベクトルの分布である。図４（Ｂ）に示した固有ベクトルの分布に於いても局所コンプライアンス行列Ｃの固有ベクトルは荷重の合計とモーメントの合計との両方が０となる対称的な分布となった。また、図４（Ｂ）と図４（Ａ）とを比較して、構造最適化後の固有ベクトルは、構造最適化前の固有ベクトルに比して振幅又は起伏が緩やかであることが把握される。これは、構造最適化により、構造体の剛性が高まり、コンプライアンスが低減された、即ち、変形の程度が小さくなったことを示している。
【００５６】
図７は、図３に示す直方体から上記の手順に従って設計変数の更新を１００回実行した後に得られた構造体の形状を示している。図に於いては、材料密度が０．５以上になった要素が表示されている。色の濃さは、材料密度の高さに対応する。図示の如く、不定荷重の条件下にて支持構造のない構造体に於いて剛性がより高く且つ体積制約条件を満たす形状が得られ、本発明の装置により支持構造のない構造体の構造最適化が可能であることが示された。
【００５７】
以上の説明は、本発明の実施の形態に関連してなされているが、当業者にとつて多くの修正及び変更が容易に可能であり、本発明は、上記に例示された実施形態のみに限定されるものではなく、本発明の概念から逸脱することなく種々の装置に適用されることは明らかであろう。
【図１】
2020016936.tif 000003
【図２】
2020016936.tif 000004
【図３】
2020016936.tif 000005
【図４】
2020016936.tif 000006
【図５】
2020016936.tif 000007
【図６】
2020016936.tif 000008
【図７】
2020016936.tif 000009