JP2004199160A

JP2004199160A - Optimum design method

Info

Publication number: JP2004199160A
Application number: JP2002363916A
Authority: JP
Inventors: Teruyoshi Washisawa; 輝芳鷲澤; Akira Asai; 朗浅井; Katsuhiko Shinjo; 克彦新庄; Masayoshi Tachihara; 昌義立原; Nobuhiro Yoshikawa; 暢宏吉川
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 2002-12-16
Filing date: 2002-12-16
Publication date: 2004-07-15
Anticipated expiration: 2022-12-16
Also published as: JP4101046B2

Abstract

<P>PROBLEM TO BE SOLVED: To provide a method for removing a structure element not contributing to an evaluation function in an automatic design technology for shape optimization including especially topology of a structure member, concerning a solution of an optimization problem of a design parameter. <P>SOLUTION: This optimum design method is equipped with a first solution process for solving the optimization problem of a first evaluation function to a state variable vector x by using a design variable vector f as a parameter, and a second solution process (S102) for solving the optimization problem of a second evaluation function L<SB>2</SB>to the state variable vector determined by the first solution process and the design variable vector. The method is provided with an erasure process (S103-108) for erasing a component f(s) corresponding to the structure element s not contributing to the second evaluation function L<SB>2</SB>from the design variable vector f. <P>COPYRIGHT: (C)2004,JPO&NCIPI

Description

【０００１】
【発明の属する技術分野】
本発明は、設計パラメタの最適化問題の解法に関し、特に構造部材のトポロジーを含む形状最適化のための自動設計技術に関するものである。
【０００２】
【従来の技術】
構造トポロジー最適設計とは、与えられた条件の下で最適な、構造部材のトポロジーおよび形状寸法を決定する問題である。以下、構造部材のトポロジーおよび形状寸法を設計変関数といい、上記決定問題を設計変関数最適化問題という。
変関数という理由は、トポロジーおよび形状寸法が３次元空間の関数になっているからである。
【０００３】
設計変関数最適化問題では、各設計変関数の値に対して、状態変関数の最適化問題を解かなければならない。この意味から、構造トポロジー最適設計は、内側に状態変関数の最適化問題、外側に設計変関数の最適化問題を持つ、２重構造最適化問題と捉える事が出来る。内側の状態変関数最適化問題では、従来技術の蓄積により、空間を有限個の要素に分割するという考え方が採用される。
【０００４】
特に構造部材の歪エネルギーを評価汎関数としている問題では、その解析手法として有限要素法が一般的である。有限要素法の解法としては一次方程式に対する直接法が採用されている。一方、設計変関数の最適化問題に対しては、大別すると、以下に示す３種類の方法が提供されている（〔文献１〕或いは〔文献２〕）：
１．Ｅｖｏｌｕｔｉｏｎａｒｙ法（以下、Ｅ法）
２．Ｈｏｍｏｇｅｎｉｚａｔｉｏｎ法（以下、Ｈ法）
３．Ｍａｔｅｒｉａｌｄｉｓｔｒｉｂｕｔｉｏｎ法（以下、ＭＤ法）
Ｅ法では、空間を分割することによって得られる部分空間のそれぞれをセルと称し、その生成と消滅を適当な規則によって繰り返す。構造部材は、最終的に存在しているセルの集合として与えられる。セルが存在するか否かという２つの状態のみを許すことにより、明確な構造部材が得られる。また、評価汎関数の微分情報を用いないので、局所最適解にトラップされないことから、評価汎関数が多峰性の場合に有効である。
【０００５】
〔文献３〕では、Ｅ法の一種である遺伝的探索法を用いた骨組構造部材の最適化設計装置が提供されている。この最適化設計装置では、骨組部材断面寸法などの離散設計変数データの近似式を使用する近似最適化計算装置と、該設計変数データを使用する詳細最適化計算装置を設け、これら２つの計算装置を結合して、膨大な設計変数が存在する実設計問題に対応することが出来るようにしている。
【０００６】
Ｈ法は、分割された各部分空間に位置する構造要素に、更に微細な構造を仮定し、連続値を取る設計変関数を新たに導入することによって、感度解析の採用を可能にしている。ここで感度解析とは、設計変関数に関する評価汎関数の微分情報を利用した解析手法のことであり、感度解析が可能になれば、勾配法のような反復解法を用いることが出来、Ｅ法のような総当り的手法に比べて、少なくとも局所最適解の探索に掛かる計算時間は大幅に短縮される。
【０００７】
ＭＤ法は、構造部材の存在確率を示す０から１の範囲の実数を各要素に割当てることによって、構造部材のトポロジーや形状寸法変化を表現する方法である。
構造部材が存在するか否かという離散的な情報を存在確率という連続値に置き換えることによって感度解析を可能にしたという意味でＨ法と同様のものであるが、Ｈ法よりパラメタ数が少ない分、モデル化が容易であり適用範囲も広い。
【０００８】
〔文献５〕にはＭＤ法による構造物の位相最適化手法が開示されている。本手法の特徴は以下のとおりである：
（１）ボクセル有限要素法（空間を等間隔に分割）を用いているため、あらゆる要素に対する要素剛性マトリクスが同じである。従って、要素剛性マトリクスを予め１度だけ計算しておけば、以後の計算に利用できる。更に、要素が規則正しく配置されているため、各要素の節点番号情報を記憶する必要がない。
【０００９】
（２）大規模連立一次方程式を解くために、前処理付共役勾配法とＥｌｅｍｅｎｔ−ｂｙ−Ｅｌｅｍｅｎｔ法を組み合わせて用いたことにより、全体剛性マトリクスを組み立てることなく解が求められるので、必要とするメモリ容量が少なくて済む。
【００１０】
（３）均質化法では、１要素に対して６個の設計変数（３次元の場合）が必要になる。更に設計変数が変化するたびに要素剛性マトリクスを再計算しなければならない。一方、構造部材の存在率を密度比で表現する密度法を採用することによって、１要素に対して１つの設計変数でよい。また設計変数の変化が要素剛性マトリクスに影響を与えない。
〔文献１〕Ｓ．Ｂｕｌｍａｎ，Ｊ．Ｓｉｅｎｚ，Ｅ．Ｈｉｎｔｏｎ：“Ｃｏｍｐａｒｉｓｏｎｓｂｅｔｗｅｅｎａｌｇｏｒｉｔｈｍｓｆｏｒｓｔｒｕｃｔｕｒａｌｔｏｐｏｌｏｇｙｏｐｔｉｍｉｚａｔｉｏｎｕｓｉｎｇａｓｅｒｉｅｓｏｆｂｅｎｃｈｍａｒｋｓｔｕｄｉｅｓ”，ＣｏｍｐｕｔｅｒｓａｎｄＳｔｒｕｃｔｕｒｅｓ，７９，ｐｐ．１２０３−１２１８（２００１）．
〔文献２〕Ｙ−Ｌ．Ｈｓｕ，Ｍ−Ｓ．Ｈｓｕ，Ｃ−Ｔ．Ｃｈｅｎ：“Ｉｎｔｅｒｐｒｅｔｉｎｇｒｅｓｕｌｔｓｆｒｏｍｔｏｐｏｌｏｇｙｏｐｔｉｍｉｚａｔｉｏｎｕｓｉｎｇｄｅｎｓｉｔｙｃｏｎｔｏｕｒｓ”，ＣｏｍｐｕｔｅｒｓａｎｄＳｔｒｕｃｔｕｒｅｓ，７９，ｐｐ．１０４９−１０５８（２００１）．
〔文献３〕特開２００１−１３４６２８号公報
〔文献４〕山川宏：“最適化デザイン”，計算力学とＣＡＥシリーズ９，培風館（１９９６）
〔文献５〕藤井，鈴木，大坪：“ボクセル有限要素法を用いた構造物の位相最適化”，ＴｒａｎｓａｃｔｉｏｎｓｏｆＪＳＣＥＳ，ＰａｐｅｒＮｏ．２０００００１０（２０００）．
【００１１】
【発明が解決しようとする課題】
しかしながら、従来技術には以下に述べるような課題があった。
【００１２】
〔文献５〕の結果の図に見られるように、総歪エネルギーに寄与しない構造要素領域、即ち浮島領域や突起領域が残っている。最終的にはこれらの構造要素領域は除去されるべきものであるが、その方法については開示されていない。
【００１３】
例えば図４において、要素５１２、５１３、５１４には荷重が加えられていなく、かつ支持もされていない場合には、これら要素は強度上、何の寄与もしていないことになる。即ち総歪エネルギーに寄与しない要素ということになり、構造部材の総重量をなるべく少なくしたいという設計要求から言えば除去すべきである。しかし、単純な幾何学的情報、即ち、特性関数の連結性のみに基づく手法を用いる限り、これらの要素を見出すことは困難である。
【００１４】
【課題を解決するための手段】
上記課題を解決するために、本発明によれば、最適設計方法において、設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程と、前記設計変数ベクトルより、前記第２の評価関数に寄与しない構造要素に対応する成分を消去する消去工程とを備える。
【００１５】
【発明の実施の形態】
以下の説明のために対象とする問題の定式化を行う。
【００１６】
有限要素法による定式化により変関数は有限次元ベクトルで表現されるものとすると、変関数の評価汎関数は、変数ベクトルの評価関数となる。以下、有限次元ベクトルで表現されているとして記述する。
【００１７】
状態変数ベクトルｘ、設計変数ベクトルｆをそれぞれ列ベクトルとして以下のように書く：
（式１）ｘ＝（ｘ（１），ｘ（１），…，ｘ（ｍ））^Ｔ
（式２）ｆ＝（ｆ（１），ｆ（２），…，ｆ（ｎ））^Ｔ
ここでＴは転置を表す。ｘはｍ次元ベクトル、ｆはｎ次元ベクトルである。
【００１８】
ｘおよびｆの境界条件をそれぞれＢ_１、Ｂ_２とする：
（式３）Ｂ_１（ｘ）＝０
（式４）Ｂ_２（ｆ）＝０
状態変数ベクトルｘおよび設計変数ベクトルｆに関する評価関数を、それぞれＬ_１、Ｌ_２とする：
（式５）Ｌ_１＝Ｌ_１（ｘ；ｆ）
（式６）Ｌ_２＝Ｌ_２（ｆ，ｘ）
Ｌ_１はｘを変数ベクトル、ｆをパラメタとする関数、Ｌ_２はｘおよびｆを変数ベクトルとする関数である。
【００１９】
一般の最適化問題では、Ｋ_１個の等式拘束条件とＫ_２個の不等式拘束条件とが課せられることが多い。そこで、設計変数に関するｊ番目の等式拘束条件およびｋ番目の不等式拘束条件を、それぞれＱ_ｊ、Ｒ_ｋとする：
（式７）Ｑ_ｊ（ｆ，ｘ）＝０
（式８）Ｒ_ｋ（ｆ，ｘ）≧０
以上の表記より最適設計問題は、拘束条件（式４）、（式５）、（式６）、（式７）、および（式８）の下で、次式を満足する解として与えられる：
（式９）ｍｉｎ［Ｌ_２］
ただし、状態変数ｘは拘束条件（式３）の下で次式を満足する解として得られる：
（式１０）ｍｉｎ［Ｌ_１］
上記定式化に基づき、本実施形態の処理を説明する。
【００２０】
図１には、実行時の流れ図が示されている。図中、ステップＳ１０１は、シミュレーション対象となる系の諸元を読み込む処理である。諸元の読み込みは、入力装置２０３或いは通信装置２０６からの入力データで行っても良いし、予め２次記憶装置２０５にファイルとして格納しておいたデータを読み出して利用することもできる。系の諸元にはｘ、ｆの初期値、境界条件Ｂ_１、Ｂ_２、評価関数Ｌ_１、Ｌ_２、拘束条件Ｑ_ｊ、Ｒ_ｋが含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保し、値を設定する。
【００２１】
ステップＳ１０２では、（式４）乃至（式９）で定式化された不等式拘束条件付最適化問題の求解工程を実施する。該求解工程には、以下に示すように、幾つかの方法が提供されている：
（１）逐次線形計画法
（２）実行可能方向法
（３）傾斜投影法
（４）一般縮約傾斜法
（５）最適基準法
以下、ステップＳ１０３から１０８が不要要素除去工程である。
【００２２】
まずステップＳ１０３では、要素ｓを１に初期化する。
【００２３】
ステップＳ１０４ではｆ（ｓ）が０かどうかを検査し、そうであればステップＳ１０８へ、そうでなければステップＳ１０５へ進む。
【００２４】
ステップＳ１０５では、ｆ（ｓ）の値を変化させる：
（式１１）ｆ（ｓ）←ｆ（ｓ）＋δ
但しδは、更新したｆ（ｓ）が０以上１以下の値を取るような、０でない任意の実数である。
【００２５】
ステップＳ１０６では、ｆ（ｓ）に関する第２の評価関数の感度∂Ｌ_２／∂ｆ（ｓ）を計算する。感度は次のように導出される：
∂Ｌ_２／∂ｆ＝∂／∂ｆ［（１／２）Ｕ^ＴＡＵ］
剛性方程式ＡＵ＝ｂより
（式１２）∂Ｌ_２／∂ｆ＝∂／∂ｆ［（１／２）ｂ^ＴＵ］＝（１／２）ｂ^Ｔ∂Ｕ／∂ｆ
一方、剛性方程式ＡＵ＝ｂの両辺をｆで変微分すると、荷重ｂはｆに依らず一定であるとすれば、
（∂Ａ／∂ｆ）Ｕ＋Ａ（∂Ｕ／∂ｆ）＝０
より
（式１３）∂Ｕ／∂ｆ＝−Ａ^−１∂Ａ／∂ｆＵ
（式１３）を（式１２）に代入し、Ｕ＝Ａ^−１ｂの関係式を用いると、
（式１４）∂Ｌ_２／∂ｆ＝−（１／２）Ｕ^Ｔ（∂Ａ／∂ｆ）Ｕ
要素ｓの特性関数について書けば
（式１５）∂Ｌ_２／∂ｆ（ｓ）＝−（１／２）Ｕ_ｓ ^Ｔ（∂Ａ_ｓ／∂ｆ（ｓ））Ｕ_ｓ
ただし、Ｕ_ｓは要素ｓに属するノード上の変位より構成されるベクトル、Ａ_ｓはＵ_ｓに対応する要素剛性マトリクスである。
【００２６】
∂Ａ_ｅ／∂ｆ（ｅ）は解析的に解ける場合とそうでない場合がある。解析解がある場合は、その計算工程を利用することができる。例えばＡ_ｓ＝ｆ（ｓ）Ａ_ｓ’のときには∂Ｌ_２／∂ｆ（ｓ）は次式で与えられる：
（式１６）∂Ｌ_２／∂ｆ（ｓ）＝−（１／２）Ｕ_ｓ ^ＴＡ_ｓ’Ｕ_ｓ
∂Ａ／∂ｆ（ｓ）が解析的に解けない場合は、自動微分という技術を利用して計算することが出来る。自動微分技術は〔文献６〕等によって公知である。
〔文献６〕久保田光一，伊理正夫：“アルゴリズムの自動微分と応用”，現代非線形科学シリーズ，コロナ社（１９９８年）
ここで、（式１６）で与えられる感度は、ステップＳ１０２で実行された求解工程の結果、全ての要素に対して０或いはそれに近い微小な値になっている。
【００２７】
要素ｓの特性関数値ｆ（ｓ）が（式１１）に従って変化したとき、（式１６）の値が変化するかどうかについて説明する。
【００２８】
先ず、要素ｓが歪エネルギーの値に寄与している場合、即ち（１／２）Ｕ_ｓ ^ＴＡ_ｓＵ_ｓが０でない場合には、ｆ（ｓ）の変化によってＡ_ｓが変化し、従って要素ｓの歪エネルギー（１／２）Ｕ_ｓ ^ＴＡ_ｓＵ_ｓが変化する。これより、構造解析問題の解として得られるＵ_ｓがｆ（ｓ）の変化の前後で異なる値Ｕ_ｓ’を取る。これより（式１６）は以下のように、負の値を取ることになる：
（式１７）∂Ｌ_２／∂ｆ’（ｓ）＝−（１／２）Ｕ_ｓ’^ＴＡ_ｓ’Ｕ_ｓ’≠ −（１／２）Ｕ_ｓ ^ＴＡ_ｓ’Ｕ_ｓ＝０
一方、要素ｓが歪エネルギーの値に寄与しない、即ち（１／２）Ｕ_ｓ ^ＴＡ_ｓＵ_ｓが０である場合には、Ａ_ｓが正定行列であることからＵ_ｓが０となり、ｆ（ｓ）が変化しても（１／２）Ｕ_ｓ ^ＴＡ_ｓＵ_ｓは変化しない。これより、構造解析問題の解として得られるＵｓがｆ（ｓ）の変化の前後で同じ値を取り、（式１６）は０になる。
【００２９】
このように、ｆ（ｓ）の値を変化させたときの感度値の絶対値が０或いは０に近い微小な値かどうかを検査することによって、該構造要素が不要かどうかを判定することができる。
【００３０】
ステップＳ１０７では、該感度∂Ｌ_２／∂ｆ（ｓ）が０かどうかを検査し、そうであればｆ（ｓ）を０に更新する。そうでなければｆ（ｓ）を更新しない。
【００３１】
ステップＳ１０８ではｓをｓ＋１に更新し、更新したｓがｎを超えたら終了し、そうでなければステップＳ１０４へ進む。
【００３２】
上記方法は、不要要素除去工程を後処理として用いた場合であるが、該第２の求解工程が反復解法の場合には、その反復ループの中で処理しても構わない。この処理について図３を用いて説明する。反復処理として共役勾配法を例に取る。
【００３３】
ステップＳ３０１では初期化を行う。具体的には、該設計変数ベクトルｆを予め与えられた値に設定し、それをｆ_０と書く。
【００３４】
ステップＳ３０２ではｔを１に設定する。
【００３５】
ステップＳ３０３では該設計変数ベクトルｆに関する該第２の評価関数のグラディエントベクトルｇ_ｔのｆ＝ｆ_０における値を計算する：
（式１８）ｇ_ｔ≡∂Ｌ_２／∂ｆ＝（∂Ｌ_２／∂ｆ（１），∂Ｌ_２／∂ｆ（２），…，∂Ｌ_２／∂ｆ（ｎ））^Ｔ
ステップＳ３０４では、ｇ_ｔが等式拘束条件及び不等式拘束条件を満たすように修正する。
【００３６】
ステップＳ３０５では、次式で計算されるのノルムが予め設定された値を超えるかどうか検査し、超えれば処理を終了し、そうでなければステップＳ３０６へ進む：
（式１９）‖ｇ_ｔ‖＝（ｇ_ｔ ^Ｔｇ_ｔ）^−１／２
ステップＳ３０６では次式で定義されるβを計算する：
（式２０）β＝‖ｇ_ｔ‖／‖ｇ_ｔ−１‖
ただしｔ＝１のときはβ＝０とする。
【００３７】
ステップＳ３０７では探索ベクトルｐ_ｔを次式で計算する：
（式２１）ｐ_ｔ＝βｐ_ｔ−ｇ_ｔ
ステップＳ３０８では、ｐ_ｔが等式拘束条件及び不等式拘束条件を満たすように修正する。
【００３８】
ステップＳ３０９では、ｐ_ｔに沿ってライン探索を行い、該第２の評価関数を極小にするｆを見つけ、それをｆ_ｔとする。
【００３９】
ステップＳ３１０では、ｆ_ｔが等式拘束条件及び不等式拘束条件を満たすように修正する。
【００４０】
ステップＳ３１１では、上記不要要素除去工程を実行する。
【００４１】
ステップＳ３１２ではｔをｔ＋１に更新し、予め設定された値を超えたら終了、そうでなければステップＳ３０３へ進む。
【００４２】
尚、ステップＳ３０４、ステップＳ３０８およびステップＳ３１０で実行される処理は、傾斜投影法として〔文献４〕等に開示されている処理を用いることができる。
【００４３】
（実施例）
本実施例は、任意の位置に加重を受ける片持ち梁の最適形状自動設計に上記実施形態を適用するものである。説明を簡単にするために平面歪問題に限定する。
【００４４】
図５に示すように、構造部材の存在を可能とする設計領域は長方形であり、有限要素法に従って、該領域を縦ｎ_ｙ、横ｎ_ｘに等間隔に分割する。分割された部分領域をセルと呼び、左下および右上のセルが（１，１）、（ｎ_ｙ，ｎ_ｘ）となるように番号付けを行う。
【００４５】
同様に格子点をノードと呼び、左下および右上のノードが（１，１）、（ｎ_ｙ＋１，ｎ_ｘ＋１）となるように番号付けを行う。
【００４６】
セル（ｊ，ｋ）には特性関数値ｆ（ｊ，ｋ）が対応する。ここで特性関数値とはセル（ｊ，ｋ）における構造部材の存在率を示す０から１の正の実数値をとる変数であり、本実施例における設計変数ベクトルｆの要素である：
（式２２）ｆ＝（ｆ（１，１），ｆ（１，２），…，ｆ（ｎ_ｙ，ｎ_ｘ））^Ｔ
同様にノード（ｊ，ｋ）には横方向変位ｕ（ｊ，ｋ）と縦方向変位ｖ（ｊ，ｋ）が対応する。これらは任意の値を取る実数であり、本発明における状態変数ベクトルＵの要素である：
（式２３）Ｕ＝（ｕ（１，１），ｖ（１，１），ｕ（１，２），ｖ（１，２），…，ｕ（ｎ_ｙ＋１，ｎ_ｘ＋１），ｖ（ｎ_ｙ＋１，ｎ_ｘ＋１））^Ｔ
同様に剛性マトリクスをＡ、加重ベクトルをｂと書けば、有限要素法の良く知られた結果として、状態変数ベクトルＵは、次式で与えられる評価関数の最適化問題の解として与えられる：
（式２４）Ｌ_１＝（１／２）Ｕ^ＴＡＵ−ｂ^ＴＵ
より具体的には以下の線型方程式の解として与えられることが知られている：
（式２５）ＡＵ＝ｂ
（式２５）の求解法としては、直接解法、反復解法等が提供されている。ただし直接解法を用いる場合は、行列Ａがフルランクでなくなる場合があるので、有限要素分割を再構築し直す必要がある。
【００４７】
設計変数ベクトルｆに対する評価関数Ｌ２は総歪エネルギーで定義される：
（式２６）Ｌ_２＝（１／２）Ｕ^ＴＡＵ
上記（式２６）を最小化するｆを求める問題を構造最適化問題と称する。通常この種の問題には、等式拘束条件として、総重量一定、即ち
【外１】

【００４８】
及び、不等式拘束条件として、特性関数値の取り得る値の範囲に対する拘束、即ち
（式２８）０≦ｆ（ｊ，ｋ）≦１
が伴う。
【００４９】
このような不等式拘束条件付最適化問題は、前述のように、公知の求解法によって解くことが出来る。
【００５０】
次に、本実施例における不要要素除去工程について説明する。
【００５１】
特に図１のステップＳ１０６で実行される感度の計算は次式を用いて行うことが出来る：
（式２９）∂Ｌ_２／∂ｆ（ｊ，ｋ）＝−（１／２）Ｕ_ｊ，ｋ ^ＴＡ_ｊ，ｋＵ_ｊ，ｋ
ただしＵ_ｊ，ｋ、Ａ_ｊ，ｋは、それぞれ、要素（ｊ，ｋ）に属するノード上の変位を成分として持つ要素変位ベクトル、該要素変位ベクトルに対応する要素剛性マトリクスである。
【００５２】
以下、図１に示した方法によって浮島要素の除去を行うことが出来る。
【００５３】
図６には不要要素除去工程を施す前の構造部材形状、図７には不要要素除去工程を施した後の構造部材形状を示す。図６の浮島領域や突起が、図７では消滅しているのが確認できる。これらの例では、感度値の絶対値が正確に０に等しいときのみ、不要要素と判定して除去している。不要要素除去工程の計算時間は、Ｐｅｎｔｉｕｍ（R)III（９３３ＭＨｚ）で４０秒程度である。
【００５４】
尚、本発明は、単一の機器からなる装置に適用しても、複数の機器から構成されるシステムに適用してもよい。また、上述した実施形態の機能を実現するソフトウェアのプログラムコードを記憶した記憶媒体を、装置あるいはシステムに供給し、装置あるいはシステム内のコンピュータが記憶媒体に格納されたプログラムコードを読み出して実行することによって達成してもよい。
【００５５】
更に、装置あるいはシステム内のコンピュータが記憶媒体に格納されたプログラムコードを読み出して実行することによって、上述した実施形態の機能を直接実現するばかりでなく、そのプログラムコードの指示に基づいて、コンピュータ上で稼動しているＯＳなどの処理により、上述の機能を実現される場合も含まれる。
【００５６】
これらの場合、そのプログラムコードを記憶した記憶媒体は本発明を構成することになる。
【００５７】
以下、上記実施形態に係わる本発明の特徴を整理する。
【００５８】
特徴１．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、
前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程と、
前記設計変数ベクトルより、前記第２の評価関数に寄与しない構造要素に対応する成分を消去する消去工程とを備えることを特徴とする最適設計方法。
【００５９】
特徴２．
前記設計変数ベクトルは、各要素における構造部材の存在率であることを特徴とする特徴１に記載の最適設計方法。
【００６０】
特徴３．
前記消去工程では、前記設計変数ベクトルの各成分の値を増減させたとき、当該設計変数ベクトルに関する前記第２の評価関数の感度の絶対値が予め設定された値より小さくなる成分を消去することを特徴とする特徴１に記載の最適設計方法。
【００６１】
特徴４．
前記消去工程は、前記第２の求解工程で計算される感度ベクトルが０である成分に対応する要素であって、かつ構造要素の存在率が０でない要素に対して行われることを特徴とする特徴１に記載の最適設計方法。
【００６２】
特徴５．
前記消去工程は、前記第２の求解工程の反復処理の所定回に１回行われることを特徴とする特徴１に記載の最適設計方法。
【００６３】
特徴６．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解手段と、
前記第１の求解手段で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解手段と、
前記設計変数ベクトルより、前記第２の評価関数に寄与しない構造要素に対応する成分を消去する消去手段とを備えることを特徴とする最適設計装置。
【００６４】
特徴７．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、
前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程と、
前記設計変数ベクトルより、前記第２の評価関数に寄与しない構造要素に対応する成分を消去する消去工程とを備えることを特徴とする最適設計プログラム。
【００６５】
【発明の効果】
以上説明したように、本発明によれば、評価関数値を最適化させる変数を求める際に、評価関数に寄与しない構造要素に対応する成分を消去することができる。
【図面の簡単な説明】
【図１】実施形態における全体処理手順を示すフローチャートである。
【図２】本実施形態に係る装置のブロック構成図である。
【図３】他の実施形態における全体処理手順を示すフローチャートである。
【図４】不要要素の説明図である。
【図５】実施例の問題設定の説明図である。
【図６】不要要素除去前の計算結果を示す図である。
【図７】不要要素除去後の計算結果を示す図である。[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a solution of a design parameter optimization problem, and more particularly to an automatic design technique for shape optimization including a topology of a structural member.
[0002]
[Prior art]
Structural topology optimal design is the problem of determining the optimal structural member topology and geometry under given conditions. Hereinafter, the topology and shape dimensions of the structural members are referred to as design variable functions, and the above-described determination problem is referred to as a design variable function optimization problem.
The reason for the variable function is that the topology and the geometric dimensions are functions of the three-dimensional space.
[0003]
In the design variable function optimization problem, the optimization problem of the state variable function must be solved for each value of the design variable function. In this sense, the structural topology optimal design can be considered as a dual structure optimization problem having an optimization problem of a state variable function inside and an optimization problem of a design variable function outside. The inner state transformation function optimization problem adopts the concept of dividing a space into a finite number of elements by accumulation of the conventional technology.
[0004]
In particular, in a problem in which the strain energy of a structural member is used as an evaluation functional, a finite element method is generally used as an analysis method. As a solution of the finite element method, a direct method for a linear equation is adopted. On the other hand, for the optimization problem of the design variable function, there are roughly three types of methods shown below ([Reference 1] or [Reference 2]):
1. Evolutionary method (hereinafter referred to as E method)
2. Homogenization method (hereinafter, H method)
3. Material distribution method (hereinafter, MD method)
In the E method, each of the subspaces obtained by dividing the space is called a cell, and its generation and extinction are repeated according to appropriate rules. The structural members are given as a set of cells that ultimately exist. By allowing only two states, the presence or absence of a cell, a distinct structural member is obtained. Further, since differential information of the evaluation functional is not used, it is not trapped in the local optimum solution, and thus is effective when the evaluation functional has multimodality.
[0005]
[Reference 3] provides an apparatus for optimizing the design of a skeleton structural member using a genetic search method, which is a type of the E method. This optimization design device includes an approximate optimization calculation device using an approximate expression of discrete design variable data such as a frame member cross-sectional dimension, and a detailed optimization calculation device using the design variable data. Are combined so that an actual design problem in which a large number of design variables exist can be dealt with.
[0006]
The H method makes it possible to adopt sensitivity analysis by assuming a finer structure in the structural elements located in each of the divided subspaces and introducing a new design variable function that takes continuous values. Here, the sensitivity analysis is an analysis method using differential information of an evaluation functional regarding a design variable function. If sensitivity analysis becomes possible, an iterative solution such as a gradient method can be used. Compared with the brute force method as described above, at least the calculation time required for searching for a local optimum solution is greatly reduced.
[0007]
The MD method is a method of expressing a topology or a change in shape and size of a structural member by assigning a real number in the range of 0 to 1 indicating the existence probability of the structural member to each element.
It is similar to the H method in that sensitivity analysis is enabled by replacing discrete information on whether or not a structural member exists with a continuous value of the existence probability, but the number of parameters is smaller than that of the H method. It is easy to model and has a wide application range.
[0008]
[Reference 5] discloses a method of optimizing the phase of a structure by the MD method. The features of this method are as follows:
(1) Since the voxel finite element method (space is divided at equal intervals) is used, the element rigidity matrix for all elements is the same. Therefore, if the element rigidity matrix is calculated only once in advance, it can be used for subsequent calculations. Further, since the elements are regularly arranged, there is no need to store the node number information of each element.
[0009]
(2) A solution is required without assembling the entire stiffness matrix by using a combination of the preconditioned conjugate gradient method and the Element-by-Element method to solve a large-scale simultaneous linear equation. Requires less memory capacity.
[0010]
(3) The homogenization method requires six design variables (in the case of three dimensions) for one element. Further, the element stiffness matrix must be recalculated every time the design variable changes. On the other hand, by adopting the density method of expressing the abundance ratio of the structural members by the density ratio, one design variable may be used for one element. Also, changes in design variables do not affect the element stiffness matrix.
[Reference 1] S.M. Bulman, J .; Sienz, E .; Hinton: “Comparisons between algorithms for structural topology optimization optimization using a series of benchmark studies”, Computers and Structures, 79. 1203-1218 (2001).
[Reference 2] YL. Hsu, MS. Hsu, CT. Chen: “Interpreting results from topology optimization using density contours”, Computers and Structures, 79, pp. 1049-1058 (2001).
[Reference 3] Japanese Patent Application Laid-Open No. 2001-134628 [Reference 4] Hiroshi Yamakawa: "Optimization Design", Computational Mechanics and CAE Series 9, Baifukan (1996)
[Reference 5] Fujii, Suzuki, Otsubo: "Topological optimization of structures using voxel finite element method", Transactions of JSCES, Paper No. 20000010 (2000).
[0011]
[Problems to be solved by the invention]
However, the prior art has the following problems.
[0012]
As can be seen from the result of [Reference 5], a structural element region that does not contribute to the total strain energy, that is, a floating island region or a protrusion region remains. Eventually, these structural element regions should be removed, but the method is not disclosed.
[0013]
For example, in FIG. 4, if the

elements

512, 513, 514 are not loaded and are not supported, they will not contribute any strength. That is, it is an element that does not contribute to the total strain energy, and should be eliminated in view of the design requirement that the total weight of the structural member be reduced as much as possible. However, it is difficult to find these elements as long as a technique based only on simple geometric information, that is, connectivity of characteristic functions is used.
[0014]
[Means for Solving the Problems]
According to the present invention, in order to solve the above problems, in an optimal design method, a first solution step of solving an optimization problem of a first evaluation function for a state variable vector using a design variable vector as a parameter; Contributing to the second evaluation function from the design variable vector, a second solution step for solving an optimization problem of a second evaluation function for the state variable vector obtained in the solution step and the design variable vector. And an erasing step of erasing components corresponding to the structural elements not to be performed.
[0015]
BEST MODE FOR CARRYING OUT THE INVENTION
For the following description, the problem to be addressed is formulated.
[0016]
Assuming that the variable function is represented by a finite-dimensional vector by formulation using the finite element method, the evaluation functional of the variable function becomes the evaluation function of the variable vector. Hereinafter, description will be made assuming that it is represented by a finite-dimensional vector.
[0017]
Write the state variable vector x and the design variable vector f as column vectors, respectively, as follows:
(Equation 1) x = (x (1), x (1),..., X (m)) ^T
(Equation 2) f = (f (1), f (2),..., F (n)) ^T
Here, T represents transposition. x is an m-dimensional vector and f is an n-dimensional vector.
[0018]
Let the boundary conditions of x and f be B ₁ and B ₂ respectively:
(Equation 3) B ₁ (x) = 0
(Equation 4) B ₂ (f) = 0
Assume that the evaluation functions for the state variable vector x and the design variable vector f are L ₁ and L ₂ respectively:
(Equation 5) L ₁ = L ₁ (x; f)
(Equation 6) L ₂ = L ₂ (f, x)
L ₁ is variable vector x, the function to parameters of f, L ₂ is a function whose variable vector x and f.
[0019]
In general optimization problems, K ₁ equality constraints and K ₂ inequality constraints are often imposed. Thus, the j-th equality constraint and the k-th inequality constraint on the design variables are Q _j and R _k , respectively:
(Equation 7) Q _j (f, x) = 0
(Equation 8) R _k (f, x) ≧ 0
From the above notation, the optimal design problem is given as a solution satisfying the following equation under the constraints (Equation 4), (Equation 5), (Equation 6), (Equation 7), and (Equation 8):
(Equation 9) min [L ₂ ]
Where the state variable x is obtained as a solution satisfying the following equation under the constraint condition (Equation 3):
(Equation 10) min [L ₁ ]
The processing of this embodiment will be described based on the above formulation.
[0020]
FIG. 1 shows a flowchart at the time of execution. In the figure, step S101 is processing for reading specifications of a system to be simulated. The reading of the specifications may be performed using input data from the input device 203 or the communication device 206, or data previously stored as a file in the secondary storage device 205 may be read and used. The specifications of the system include initial values of x and f, boundary conditions B ₁ and B ₂ , evaluation functions L ₁ and L ₂ , constraint conditions Q _j and R _k . Based on this information, the program secures a necessary variable area in the primary storage device 204 and sets a value.
[0021]
In step S102, a step of solving the inequality constraint condition optimization problem formulated by (Equation 4) to (Equation 9) is performed. Several methods are provided for the solving step, as shown below:
(1) Sequential linear programming (2) feasible direction method (3) oblique projection method (4) general contracted oblique method (5) optimal reference method Hereinafter, steps S103 to S108 are unnecessary element removing steps.
[0022]
First, in step S103, the element s is initialized to 1.
[0023]
In step S104, it is checked whether f (s) is 0. If so, the process proceeds to step S108; otherwise, the process proceeds to step S105.
[0024]
In step S105, the value of f (s) is changed:
(Equation 11) f (s) ← f (s) + δ
Here, δ is an arbitrary non-zero real number such that the updated f (s) takes a value of 0 or more and 1 or less.
[0025]
In step S106, the sensitivity ∂L ₂ / ∂f (s) of the second evaluation function for f (s) is calculated. The sensitivity is derived as follows:
_{∂L 2 / ∂f = ∂ / ∂f} [(1/2) U T AU]
From stiffness equation AU = b (Equation _{12) ∂L 2 / ∂f = ∂} / ∂f [(1/2) b T U] = (1/2) b T ∂U / ∂f
On the other hand, if both sides of the rigidity equation AU = b are differentiated by f, assuming that the load b is constant independently of f,
(∂A / ∂f) U + A (∂U / ∂f) = 0
From (Equation 13), ∂U / ∂f = −A ⁻¹ ∂A / ∂fU
By substituting (Equation 13) into (Equation 12) and using the relational expression of U = A ⁻¹ b,
(Equation _{14) ∂L 2 / ∂f = -} (1/2) U T (∂A / ∂f) U
If the characteristic function of the element s is described (Equation 15), ΔL ₂ / Δf (s) = − (1/2) U _s ^T (ΔA _s / Δf (s)) U _s
However, U _s is the vector composed of the displacement of the nodes belonging to the element s, A _s is the element stiffness matrix corresponding to the U _s.
[0026]
∂A _e / ∂f (e) can be solved analytically or not. If there is an analytical solution, the calculation process can be used. For example ∂L ₂ / ∂f when the _{_{A s = f (s) A}} s' (s) is given by the following equation:
(Equation _{16) ∂L 2 / ∂f (s} ) = - (1/2) U s T A s' U s
If ∂A / ∂f (s) cannot be solved analytically, it can be calculated using a technique called automatic differentiation. The automatic differentiation technique is known from [Reference 6] and the like.
[Reference 6] Koichi Kubota, Masao Iri: "Automatic Differentiation and Application of Algorithm", Modern Nonlinear Science Series, Corona (1998)
Here, the sensitivity given by (Equation 16) is 0 or a very small value close to all the elements as a result of the solution calculation process executed in step S102.
[0027]
Whether the value of (Equation 16) changes when the characteristic function value f (s) of the element s changes according to (Equation 11) will be described.
[0028]
First, if the element s contributes to the value of the strain energy, i.e. _^(1/2) If U _{s T} A s _{U s} is not 0, _{A s} is changed by a change in f (s), thus strain energy of the element _{^{_{s (1/2) U s T a}}} s U s is changed. From this, U _s obtained as a solution of structural analysis problems take different values U _s' before and after the change of f (s). Thus, (Equation 16) takes a negative value as follows:
(Equation _{17) ∂L 2 / ∂f '(} s) = - (1/2) U s' T A s' U s' ≠ - (1/2) U s T A s' U s = 0
On the other hand, the element s can not contribute to the value of the strain energy, i.e. _{^{_{(1/2) U s T A s}}} U when _s is 0, becomes _{U s} is 0 since _{A s} is positive definite matrix, f (s) is also changed _{^{_{(1/2) U s T A s}}} U s is not changed. Thus, Us obtained as a solution to the structural analysis problem takes the same value before and after the change in f (s), and (Equation 16) becomes zero.
[0029]
As described above, by checking whether the absolute value of the sensitivity value when the value of f (s) is changed is 0 or a minute value close to 0, it is possible to determine whether the structural element is unnecessary. it can.
[0030]
In step S107, it is checked whether or not the sensitivity ΔL ₂ / Δf (s) is 0, and if so, f (s) is updated to 0. Otherwise, f (s) is not updated.
[0031]
In step S108, s is updated to s + 1, and if the updated s exceeds n, the process ends; otherwise, the process proceeds to step S104.
[0032]
The above-described method is a case where the unnecessary element removing step is used as a post-processing. However, when the second solving step is an iterative solution, the processing may be performed in an iterative loop. This processing will be described with reference to FIG. Take the conjugate gradient method as an example of the iterative process.
[0033]
In step S301, initialization is performed. Specifically, it sets in advance given value the design variable vector f, write it as f _0.
[0034]
In step S302, t is set to 1.
[0035]
In step S303 to calculate a value of f = _{f 0} of the gradient vector _{g t} of the evaluation function of the second relates to the design variable vector f:
(Equation _{_{18) g t ≡∂L 2 / ∂f}} = (∂L 2 / ∂f (1), ∂L 2 / ∂f (2), ..., ∂L 2 / ∂f (n)) T
In step S304, g _t is corrected so as to satisfy the equality constraints and inequality constraints.
[0036]
In step S305, it is checked whether the norm calculated by the following equation exceeds a preset value. If it does, the process ends, and if not, the process proceeds to step S306:
(Equation 19) ‖g _t || _{^{_{^{= (g t T g t)}}}} -1/2
In step S306, β defined by the following equation is calculated:
(Equation 20) β = ‖g _t ‖ / ‖g _t-1 ||
However, when t = 1, β = 0.
[0037]
In step S307 to calculate the search vector _{p t} by the following equation:
(Equation _{_{21) p t = βp t -g}} t
In step S308, p _t is corrected so as to satisfy the equality constraints and inequality constraints.
[0038]
At step S309, the perform line search along the p _t, find the f to the minimum evaluation function of the second, to make it as f _t.
[0039]
In step S310, f _t is corrected so as to satisfy the equality constraints and inequality constraints.
[0040]
In step S311, the unnecessary element removing step is performed.
[0041]
In step S312, t is updated to t + 1. If the value exceeds a preset value, the process ends. Otherwise, the process proceeds to step S303.
[0042]
The processing executed in steps S304, S308, and S310 can use the processing disclosed in [Reference 4] or the like as the oblique projection method.
[0043]
(Example)
In the present embodiment, the above-described embodiment is applied to automatic design of an optimum shape of a cantilever which receives a load at an arbitrary position. In order to simplify the explanation, the problem is limited to the plane distortion problem.
[0044]
As shown in FIG. 5, the design area to allow the presence of the structural member is rectangular, in accordance with the finite element method, is divided into equal intervals region vertical n _y, the horizontal n _x. The divided partial region is called a cell, the lower left and upper right of the cell (1,1), and numbered so that the (n y, _n _x).
[0045]
Similarly called grid points and nodes, the lower-left and upper-right nodes _(1,1), and numbered so that the _{(n y + 1, n x} +1).
[0046]
The characteristic function value f (j, k) corresponds to the cell (j, k). Here, the characteristic function value is a variable having a positive real value of 0 to 1 indicating the existence ratio of the structural member in the cell (j, k), and is an element of the design variable vector f in the present embodiment:
(Equation 22) f = (f (1, 1), f (1, 2),..., F ( _ny , _nx )) ^T
Similarly, a horizontal displacement u (j, k) and a vertical displacement v (j, k) correspond to the node (j, k). These are real numbers with arbitrary values and are elements of the state variable vector U in the present invention:
(Equation 23) U = (u (1,1 ), v (1,1), u (1,2), v (1,2), ..., u (n y + 1, n x +1), v ( _{_{n y + 1, n x +1}} )) T
Similarly, if the stiffness matrix is written A and the weight vector is written b, as a well-known result of the finite element method, the state variable vector U is given as a solution to the optimization problem of the evaluation function given by:
(Formula _{^{24) L 1 = (1/2)}} U T AU-b T U
More specifically, it is known to be given as a solution to the following linear equation:
(Equation 25) AU = b
As the solution method of (Equation 25), a direct solution method, an iterative solution method, and the like are provided. However, when the direct solution method is used, the matrix A may not be full rank, and thus it is necessary to reconstruct the finite element division.
[0047]
The evaluation function L2 for the design variable vector f is defined by the total strain energy:
(Equation 26) L ₂ = (１／) U ^T AU
The problem of finding f that minimizes the above (Equation 26) is referred to as a structure optimization problem. Usually, for this type of problem, the equation constraint condition is that the total weight is constant, ie

[0048]
And, as the inequality constraint condition, a constraint on a range of possible values of the characteristic function value, that is, (Equation 28) 0 ≦ f (j, k) ≦ 1
Is accompanied.
[0049]
Such an optimization problem with inequality constraints can be solved by a known solving method as described above.
[0050]
Next, an unnecessary element removing step in this embodiment will be described.
[0051]
In particular, the calculation of sensitivity performed in step S106 of FIG. 1 can be performed using the following equation:
(Equation 29) ∂L ₂ / ∂f (j, k) = − (）) U _{j, k} ^T _{Aj, k} U _{j, k}
However, U _{j, k} and A _{j, k} are an element displacement vector having a displacement on a node belonging to the element (j, k) as a component, and an element rigidity matrix corresponding to the element displacement vector, respectively.
[0052]
Hereinafter, the floating island element can be removed by the method shown in FIG.
[0053]
FIG. 6 shows the shape of the structural member before the unnecessary element removing step is performed, and FIG. 7 shows the shape of the structural member after the unnecessary element removing step. It can be seen that the floating island region and the protrusion in FIG. 6 have disappeared in FIG. In these examples, only when the absolute value of the sensitivity value is exactly equal to 0, it is determined as an unnecessary element and removed. The calculation time of the unnecessary element removing step is about 40 seconds in Pentium (III) III (933 MHz).
[0054]
Note that the present invention may be applied to an apparatus composed of a single device or a system composed of a plurality of devices. Further, a storage medium storing software program codes for realizing the functions of the above-described embodiments is supplied to an apparatus or a system, and a computer in the apparatus or the system reads and executes the program codes stored in the storage medium. May be achieved by:
[0055]
Furthermore, not only the functions of the above-described embodiments are directly realized by reading and executing the program code stored in the storage medium by the computer in the device or the system, but also by executing the program code on the computer based on the instruction of the program code. The above-described functions may be realized by the processing of the OS or the like running on Windows.
[0056]
In these cases, the storage medium storing the program code constitutes the present invention.
[0057]
Hereinafter, features of the present invention according to the above embodiment will be summarized.
[0058]
Features 1.
A first solution step of solving an optimization problem of a first evaluation function for the state variable vector using the design variable vector as a parameter;
A second solution step of solving an optimization problem of a second evaluation function with respect to the state variable vector obtained in the first solution step and the design variable vector;
Erasing a component corresponding to a structural element that does not contribute to the second evaluation function from the design variable vector.
[0059]
Features2.
2. The optimal design method according to Feature 1, wherein the design variable vector is an existence rate of a structural member in each element.
[0060]
Features 3.
In the erasing step, when the value of each component of the design variable vector is increased or decreased, a component in which the absolute value of the sensitivity of the second evaluation function for the design variable vector is smaller than a preset value is eliminated. The optimal design method according to Feature 1, characterized in that:
[0061]
Features 4.
The erasing step is performed on an element corresponding to a component whose sensitivity vector calculated in the second solving step is 0, and an element whose structural element existence rate is not 0. The optimal design method according to Feature 1.
[0062]
Features5.
The optimal design method according to claim 1, wherein the erasing step is performed once in a predetermined number of repetitions of the second solving step.
[0063]
Features 6.
First solving means for solving an optimization problem of a first evaluation function for the state variable vector using the design variable vector as a parameter;
Second solving means for solving an optimization problem of a second evaluation function with respect to the state variable vector obtained by the first solving means and the design variable vector;
An optimizing device comprising: an erasing unit that erases a component corresponding to a structural element that does not contribute to the second evaluation function from the design variable vector.
[0064]
Features 7.
A first solution step of solving an optimization problem of a first evaluation function for the state variable vector using the design variable vector as a parameter;
A second solution step of solving an optimization problem of a second evaluation function with respect to the state variable vector obtained in the first solution step and the design variable vector;
An erasing step of erasing components corresponding to structural elements that do not contribute to the second evaluation function from the design variable vector.
[0065]
【The invention's effect】
As described above, according to the present invention, when determining a variable for optimizing an evaluation function value, it is possible to eliminate components corresponding to structural elements that do not contribute to the evaluation function.
[Brief description of the drawings]
FIG. 1 is a flowchart illustrating an overall processing procedure according to an embodiment.
FIG. 2 is a block diagram of an apparatus according to the embodiment.
FIG. 3 is a flowchart illustrating an overall processing procedure according to another embodiment.
FIG. 4 is an explanatory diagram of an unnecessary element.
FIG. 5 is an explanatory diagram of a problem setting according to the embodiment.
FIG. 6 is a diagram illustrating a calculation result before unnecessary elements are removed.
FIG. 7 is a diagram illustrating a calculation result after removing unnecessary elements.

Claims

A first solution step of solving an optimization problem of a first evaluation function for the state variable vector using the design variable vector as a parameter;
A second solution step of solving an optimization problem of a second evaluation function with respect to the state variable vector obtained in the first solution step and the design variable vector;
Erasing a component corresponding to a structural element that does not contribute to the second evaluation function from the design variable vector.