JP2004310376A

JP2004310376A - Optimum design of structure

Info

Publication number: JP2004310376A
Application number: JP2003102139A
Authority: JP
Inventors: Teruyoshi Washisawa; 輝芳鷲澤
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 2003-04-04
Filing date: 2003-04-04
Publication date: 2004-11-04
Anticipated expiration: 2023-04-04
Also published as: JP4154271B2

Abstract

<P>PROBLEM TO BE SOLVED: To provide an optimum design method of a structure capable of executing calculation without conducting special processing, even when a total rigidity matrix gets singular, which simplifies a program thereby, and reduces the calculation volume. <P>SOLUTION: In this optimum design method of a structure for finding a solution to a structure optimum design problem formulated as a double optimization problem having the first and second solution finding processes, a state variable vector serves as a displacement in each node, a design variable vector serves as an existence rate of a structural member in each element, the first solution finding process includes the second solution finding process as one step thereof and reads out the stored design variable vector and the state variable vector to update the design variable vector, and the second solution finding process is constituted of a conjugate gradient method, preprocesses a contact point force vector based on the total rigidity matrix (S302) and reads out the design variable vector and state variable vector stored in the second storage means to update the state variable vector (S305). The state variable vector is not initialized at the time point where the second solution finding process is started. <P>COPYRIGHT: (C)2005,JPO&NCIPI

Description

【０００１】
【発明の属する技術分野】
本発明は設計パラメータの最適化問題の解法である構造物の最適設計、特に構造部材のトポロジーを含む形状最適化のための自動設計技術に関するものである。
【０００２】
【従来の技術】
構造トポロジーの最適設計とは、与えられた条件の下で、構造部材の最適なトポロジー及び形状寸法を決定する問題である。以下、構造部材のトポロジー及び形状寸法を設計変関数といい、上記決定問題を設計変関数の最適化問題という。変関数という理由は、トポロジー及び形状寸法が３次元空間の関数になっているからである。設計変関数の最適化問題では、各設計変関数の値に対して、状態変関数の最適化問題を解かなければならない。この意味から、構造トポロジーの最適設計は、内側に状態変関数の最適化問題を持ち、外側に設計変関数の最適化問題を持つ２重構造最適化問題と捉えることが出来る。
【０００３】
内側の状態変関数最適化問題では、従来技術の蓄積により、空間を有限個の要素に分割するという考え方が採用される。特に構造部材の歪エネルギーを評価汎関数としている問題では、その解析手法として有限要素法が一般的である。有限要素法の解法としては一次方程式に対する直接法が採用されている。
【０００４】
一方、設計変関数の最適化問題に対しては、大別すると、以下に示す３種類の方法が提供されている（非特許文献１、或いは、非特許文献２参照）。
１．Ｅｖｏｌｕｔｉｏｎａｒｙ法（以下、Ｅ法）
２．Ｈｏｍｏｇｅｎｉｚａｔｉｏｎ法（以下、Ｈ法）
３．Ｍａｔｅｒｉａｌｄｉｓｔｒｉｂｕｔｉｏｎ法（以下、ＭＤ法）、或いは密度（Ｄｅｎｓｉｔｙ）法（以下、Ｄ法）
Ｅ法では、空間を分割することによって得られる部分空間のそれぞれをセルと称し、その生成と消滅を適当な規則によって繰り返す。構造部材は、最終的に存在しているセルの集合として与えられる。セルが存在するか否かという２つの状態のみを許すことにより、明確な構造部材が得られる。また、評価汎関数の微分情報を用いないので、局所最適解にトラップされないことから、評価汎関数が多峰性の場合に有効である。特許文献１では、Ｅ法の一種である遺伝的探索法を用いた、骨組構造部材の最適化設計装置が提供されている。この最適化設計装置では、蓄積されたノウハウに基づく試行計算を必要とし、従って膨大な設計変数が存在する実設計問題に対応することが出来なかった従来技術の課題を、以下の方法によって解決している。即ち、骨組部材断面寸法などの離散設計変数データの近似式を使用する近似最適化計算装置と、該設計変数データを使用する詳細最適化計算装置を設け、これら２つの計算装置を結合して骨組構造の最適設計装置を構成している。
【０００５】
Ｈ法は、分割された各部分領域に位置する構造要素に、更に微細な構造を仮定し、連続値を取る設計変関数を新たに導入することによって、感度解析の採用を可能にしている。ここで感度解析とは、設計変関数に関する評価汎関数の微分情報を利用した解析手法のことであり、感度解析が可能になれば、勾配法のような反復解法を用いることが出来、Ｅ法のような総当り的手法に比べて、少なくとも局所最適解の探索に係る計算時間は大幅に短縮される（非特許文献３参照）。
【０００６】
ＭＤ法或いはＤ法は、構造部材の存在率を示す０から１の範囲の実数を各要素に割当てることによって、構造部材のトポロジーや形状寸法変化を表現する方法である。構造部材が存在するか否かという離散的な情報を存在率という連続値に置き換えることによって感度解析を可能にしたという意味でＨ法と同様のものであるが、Ｈ法よりパラメータ数が少ない分、モデル化が容易であり適用範囲も広い。
【０００７】
非特許文献４には、Ｄ法による構造物の位相最適化手法が開示されている。本手法の特徴は以下のとおりである。ボクセル有限要素法（空間を等間隔に分割）を用いているため、あらゆる要素に対する要素剛性マトリクスが同じである。従って、要素剛性マトリクスを予め１度だけ計算しておけば、以後の計算に利用できる。更に、要素が規則正しく配置されているため、各要素の節点番号情報を記憶する必要がない。大規模連立一次方程式を解くために、前処理付共役勾配法とＥｌｅｍｅｎｔ−ｂｙ−Ｅｌｅｍｅｎｔ法を組み合わせて用いたことにより、全体剛性マトリクスを組み立てることなく解が求められるので、必要とするメモリ容量が少なくて済む。
【０００８】
均質化法では、１要素に対して６個の設計変数（３次元の場合）が必要になる。更に設計変数が変化するたびに要素剛性マトリクスを再計算しなければならない。一方、構造部材の存在率を密度比で表現する密度法を採用することによって、１要素に対して１つの設計変数でよい。また設計変数の変化が要素剛性マトリクスに影響を与えない。
【０００９】
【特許文献１】
特開平１１−３１４６３１号
【非特許文献１】
Ｓ．Ｂｕｌｍａｎ，Ｊ．Ｓｉｅｎｚ，Ｅ．Ｈｉｎｔｏｎ： ”Ｃｏｍｐａｒｉｓｏｎｓｂｅｔｗｅｅｎａｌｇｏｒｉｔｈｍｓｆｏｒｓｔｒｕｃｔｕｒａｌｔｏｐｏｌｏｇｙｏｐｔｉｍｉｚａｔｉｏｎｕｓｉｎｇａｓｅｒｉｅｓｏｆｂｅｎｃｈｍａｒｋｓｔｕｄｉｅｓ，” ＣｏｍｐｕｔｅｒｓａｎｄＳｔｒｕｃｔｕｒｅｓ，７９，ｐｐ．１２０３−１２１８（２００１）．
【非特許文献２】
Ｙ−Ｌ．Ｈｓｕ，Ｍ−Ｓ．Ｈｓｕ，Ｃ−Ｔ．Ｃｈｅｎ： ”Ｉｎｔｅｒｐｒｅｔｉｎｇｒｅｓｕｌｔｓｆｒｏｍｔｏｐｏｌｏｇｙｏｐｔｉｍｉｚａｔｉｏｎｕｓｉｎｇｄｅｎｓｉｔｙｃｏｎｔｏｕｒｓ，” ＣｏｍｐｕｔｅｒｓａｎｄＳｔｒｕｃｔｕｒｅｓ，７９，ｐｐ．１０４９−１０５８（２００１）．
【非特許文献３】
山川宏：“最適化デザイン，”計算力学とＣＡＥシリーズ９，培風館（１９９６）
【非特許文献４】
藤井，鈴木，大坪：“ボクセル有限要素法を用いた構造物の位相最適化，”ＴｒａｎｓａｃｔｉｏｎｓｏｆＪＳＣＥＳ，ＰａｐｅｒＮｏ．２０００００１０（２０００）．
【００１０】
【発明が解決しようとする課題】
しかしながら、上記従来技術には以下に述べるような問題点があった。
【００１１】
一般に構造最適化問題は、設計変数ベクトルの最適化問題の各反復ステップ毎に、状態変数ベクトルの最適化問題を含む２重の最適化問題として定式化される。設計変数ベクトルの最適化問題を外側の最適化、状態変数ベクトルの最適化問題を内側の最適化と呼ぶと、内側の最適化は、設計変数ベクトルをパラメータとし、即ち設計変数ベクトルを固定し、状態変数ベクトルを求める問題である。これは通常、構造解析と呼ばれ、有限要素法によって線形方程式の解法を用いて解くことが出来る。
【００１２】
ところが、構造が変化し、ある領域の構造部材が存在しなくなる、例えば部材に穴が空いてしまうと、対応する要素の設計変数が０となり、結果的に要素のヤング率が０となる。すると、上記線形方程式の係数行列がフルランクでなくなり、逆行列が計算できないという理由で、直接法では解けなくなる。
【００１３】
そのため従来技術では、設計変数ベクトルが０となる場合に以下のように対処してきた。
【００１４】
係数行列がフルランクになるように、即ちランク数と方程式の数が等しくなるように、方程式を構成し直す。
【００１５】
設計変数ベクトルの値を正確に０とはせずに０に近い微小な値で置き換える。
【００１６】
しかるに、上記前者の方法は、設計変数ベクトルが更新される毎に方程式の更新を行う可能性が生じ、計算時間がかかる。
【００１７】
また、上記後者の方法、即ち要素の設計変数を微小な値に設定するということは、物理的には薄い膜、或いは弱い部材が存在することに相当する。つまり従来法では、物質がない部分を正確に表現していないので、得られた計算結果の精度に疑問が残る。
【００１８】
本発明は、上記従来の問題点に鑑み、全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量を軽減する構造物の最適設計方法を提供する。
【００１９】
【課題を解決するための手段】
上記課題を解決するために、本発明の構造最適設計方法は、状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解工程と、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解工程とを有する２重最適化問題として定式化される構造最適設計問題の解を求める構造最適設計方法において、前記状態変数ベクトルを各ノードにおける変位、前記設計変数ベクトルを各要素における構造部材の存在率とし、前記第１の求解工程は、その１ステップとして前記第２の求解工程を含み、第１の記憶手段に記憶された設計変数ベクトル及び状態変数ベクトルを読み出して設計変更ベクトルを更新し、更新された設計変数ベクトルを該第１の記憶手段に記憶する設計変数更新工程を含み、前記第２の求解工程は、共役勾配法により構成されて、全体剛性行列に基づいて接点力ベクトルに対して前処理を行う前処理工程と、第２の記憶手段に記憶された設計変更ベクトル及び状態変数ベクトルを読み出して状態変数ベクトルを更新し、更新された状態変数ベクトルを該第２の記憶手段に記憶する状態変数更新工程とを含み、前記状態変数ベクトルが前記第２の求解工程が開始される時点で初期化されないことを特徴とする。
【００２０】
【発明の実施の形態】
以下、本発明の構造最適設計法を実現する構造最適設計装置の構成例及び動作例を説明する。
【００２１】
＜本実施形態の構造最適設計装置の構成例＞
図１に、本実施形態の構造最適設計装置の構成例を示す。
【００２２】
図中、２０１は演算・制御用のＣＰＵ、２０２は表示装置、２０３は入力装置、２０４は１次記憶装置、２０５は２次記憶装置、２０６は通信装置、２０７はバスラインである。本実施形態における構成要素は、プログラムとして予め２次記憶装置２０５に格納されており、入力装置２０３或いは通信装置２０６等からのコマンド入力によって１次記憶装置にロードされ、ＣＰＵ２０１により実行されるものである。
【００２３】
図２に、本実施形態の構造最適設計法を実現するために、１次記憶装置２０４及び２次記憶装置２０５に記憶されるデータ及びプログラムの記憶例を示す。
【００２４】
１次記憶装置２０４は、ＲＡＭやＲＯＭからなり、データ記憶領域２１４とプログラム記憶領域２２４を有している。図２には、変更されずに常駐するＯＳやＢＩＯＳのプログラム、変更されないＯＳやＢＩＯＳ用のパラメータなどは図示されていない。データ記憶領域２１４には、本実施形態の構造最適設計を実現するための各種のデータが一時記憶される。この中で再度利用されるデータ／反復して使用されるデータは、一時記憶されることが必須である。保存・保持の必要が無いデータは、一時記憶は必須ではない。
【００２５】
例えば、データ記憶領域２１４には、ｘ（状態変数ベクトル）２１４ａ、ｆ（設計変数ベクトル）２１４ｂ、Ｌ_１（状態変数ベクトルの評価関数）２１４ｃ、Ｌ_２（設計変数ベクトルの評価関数）２１４ｄ、Ｂ１（状態変数のための境界条件）２１４ｅ、Ｗ_０（設計変数のための重量条件）２１４ｆ、Ｘ^（０）やｆ^（０）（初期値）２１４ｇ、λ（ラグランジェ未定定数）２１８ｈ、Ａ（要素剛性行列）２１４ｉ、ｋ（設計変数変更の反復回数）２１４ｊ、Ｃ（感度係数）２１４ｋ、ｒ（残差ベクトル）１２４ｍ、ｐ（探索方向ベクトル）２１５ｎ、ｒ^（０）やｐ^（０）（初期値）２１４ｐ、ｔ（状態変数変更の反復回数）２１４ｑ、ｂ（節点力ベクトル）２１４ｒ、その他のデータ２１４ｓなどが含まれる。一方、プログラム記憶領域２２４には、２次記憶装置２０５或は通信装置２０６からロードされたプログラムが記憶されて、ＣＰＵ２０１により実行される。
【００２６】
２次記憶装置２０５は、例えば、フロッピーディスクやＣＤなどの外部大容量メモリであり、着脱可能であればなお好ましい。２次記憶装置２０５も、データ記憶領域２１５とプログラム記憶領域２２５を有している。
【００２７】
データ記憶領域２１５には、詳細には示さないが、各種の初期値２１５ａや各種の条件２１５ｂ、あるいは２重最適値問題のデータベース２１５ｃが記憶されている。データベース２１５ｃには、標準化されたマトリクスなどが記憶されていて使用されてもよい。
【００２８】
プログラム記憶領域２２５には、本実施形態で設計変数ベクトルにおける求解法として使用される、最適性規準法モジュール２２５ａ、逐次凸関数近似法モジュール２２５ｂ、逐次線形計画法モジュール２２５ｃ、その他の方法のモジュール２２５ｄが記憶され、状態変数ベクトルにおける求解法として使用される、共役勾配法モジュール２２５ｅ、ＧＣＲ法モジュール２２５ｆ、ＧＣＲ（ｋ）法２２５ｇ、Ｏｒｔｈｏｍｉｎ（ｋ）法２２５ｈ、ＧＭＲＥＧ（ｋ）法２２５ｉ、その他の方法のモジュール２２５ｊが記憶されている。このプログラム記憶領域２２５に記憶されたプログラムから選択された方法に対応するモジュールが１次記憶装置２０４のプログラム記憶領域２２４にロードされて、ＣＰＵ２０１により実行され本実施形態の構造最適設計法が実現する。
【００２９】
すなわち、本構造最適設計装置は、状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解手段と、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解手段とを有し、２重最適化問題として定式化される構造最適設計問題を解く構造最適設計装置であって、前記状態変数ベクトルが各ノードにおける変位、前記設計変数ベクトルが各要素における構造部材の存在率であり、前記第２の求解手段は共役勾配法により構成されて、全体剛性行列に基づいて節点力ベクトルに対して前処理を行う前処理手段を含み、前記状態変数ベクトルが前記第２の求解手段の処理が開始される時点で初期化されない。前記前処理手段は、全体剛性行列の対角成分が０となる行あるいは列に対応する接点力ベクトルの成分を０にする。前記第１の求解手段は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかを実行する。
【００３０】
＜本実施形態の構造最適化問題の定式化＞
以下の説明のために、構造最適化問題の定式化を行う。
【００３１】
有限要素法による定式化により変関数は有限次元ベクトルで表現されるものとすると、変関数の評価汎関数は、変数ベクトルの評価関数となる。以下、有限次元ベクトルで表現されているとして記述する。
【００３２】
状態変数ベクトルｘ、設計変数ベクトルｆをそれぞれ列ベクトルとして以下のように書く。
【００３３】
ｘ＝（ｘ_１，ｘ_２，…，ｘ_ｍ）^Ｔ …（式１）
ｆ＝（ｆ_１，ｆ_２，…，ｆ_ｎ）^Ｔ …（式２）
ここで、Ｔは転置を表す。ｘはｍ次元ベクトル、ｆはｎ次元ベクトルである。構造最適化問題において、状態変数ベクトルは各ノードの変位、設計変数ベクトルは各要素内における構造物の存在率より構成される。
【００３４】
例えば、２次元面内での構造最適化問題において、設計領域を縦ｎ_ｙ、横ｎ_ｘに分割すると、
要素数ｎは（ｎ_ｘ×ｎ_ｙ）、ノード数は（ｎ_ｘ＋１）×（ｎ_ｙ＋１）である。
【００３５】
状態変数ベクトルは各ノードの縦方向及び横方向変位の組なので、ｘの次元ｍは２×（ｎ_ｘ＋１）×（ｎ_ｙ＋１）である。
【００３６】
ｘの境界条件は通常、変位拘束条件として与えられる。それをＢ_１と記述する。
Ｂ_１（ｘ）＝０ …（式３）
状態変数ベクトルｘ及び設計変数ベクトルｆに関する評価関数をそれぞれＬ_１、Ｌ_２とすると、次式のように定義される。
Ｌ_１：＝Ｌ_１（ｘ，ｆ）＝（ｂ − Ａｘ）^Ｔ（ｂ − Ａｘ） …（式４）
Ｌ_２：＝Ｌ_２（ｆ，ｘ）＝（１／２）ｘ^ＴＡｘ …（式５）
ここで、Ａは全体剛性行列、ｂは節点力ベクトルである。ただし、ｂは予め与えられるものであり、Ａは設計変数ベクトルｆの関数である。
【００３７】
全体剛性行列は要素ｊに対する要素剛性行列Ａ_ｊの重ね合わせで構成することが出来る。平面歪問題の場合、設計変数ベクトル、即ち構造要素の存在率ｆ_ｊを考慮した要素剛性行列を以下に示す。ただし、２次元有限要素として４ノード要素を採用し、要素変位ベクトルｘ_ｊを次のように構成する。
ｘ_ｊ＝（ｕ_１，ｖ_１，ｕ_２，ｖ_２，… ，ｕ_４，ｖ_４） …（式６）
ここで、ｕ_ｋ、ｖ_ｋはノードｋの水平方向、及び垂直方向の変位である。
【００３８】
（式６）の要素変位ベクトルに対応する要素剛性行列は、次式のように与えられる。
【００３９】
【数１】

【００４０】
上式中
ａ_ｉ，ｊ＝（λ＋２μ）＜∂_ｘｅ_ｉ，∂_ｘｅ_ｊ＞＋ μ＜∂_ｙｅ_ｉ，∂_ｙｅ_ｊ＞ …（式８）
ｂ_ｉ，ｊ＝ λ＜∂_ｘｅ_ｉ，∂_ｙｅ_ｊ＞＋ μ＜∂_ｘｅ_ｊ，∂_ｙｅ_ｉ＞ …（式９）
ｃ_ｉ，ｊ＝ λ＜∂_ｘｅ_ｊ，∂_ｘｅ_ｉ＞＋ μ＜∂_ｘｅ_ｉ，∂_ｙｅ_ｊ＞ …（式１０）
ｄ_ｉ，ｊ＝（λ＋２μ）＜∂_ｙｅ_ｉ，∂_ｙｅ_ｊ＞＋ μ＜∂_ｘｅ_ｉ，∂_ｘｅ_ｊ＞ …（式１１）
ここで、ｅ_ｊはノードｊに対する基底ベクトル、＜，＞はベクトルの内積を表す。∂_ｘはｘに関する偏微分、∂_ｙはｙに関する偏微分である。ξは前述の（非特許文献４）に従い、２とする。また、λ、μはラメ定数と呼ばれる材料定数であり、ヤング率Ｅ及びポアソン比νより次式で計算される。
【００４１】
λ ＝Ｅν／（１＋ν）（１−２ν） …（式１２）
μ ＝Ｅ／｛２（１＋ν）｝ …（式１３）
（式７）より容易に分かるように、ある要素の存在率が０となると要素剛性行列Ａ_ｊのすべての成分が０となる。全体剛性行列は要素剛性行列の重ね合わせで構成される。従って、ノードｊを含む要素の設計変数ｆ_ｊが０になると、全体剛性行列のノードｊに対応する行及び列の全ての成分が０になる。このように、構造最適化問題において全体剛性行列は一般に特異な行列となる。
【００４２】
設計変数に関する拘束条件は、構造最適設計の場合、総重量の上限値Ｗ_０とする。
【００４３】
【数２】

【００４４】
また、設計変数には値域に対する拘束条件が以下のように定められている。
【００４５】
０ ≦ ｆ_ｊ ≦ １ｊ＝１， … ，ｎ …（式１５）
以上の表記より、最適設計問題は、拘束条件（式１４）及び（式１５）の下で、（式５）を極小化する問題として定式化される。
【００４６】
また、状態変数ｘについては、拘束条件（式３）の下で（式４）を極小化する解として得られる。
【００４７】
＜本実施形態の構造最適設計装置の動作例＞
上記定式化に基づき、本実施形態の処理手順例を説明する。
【００４８】
図３には、本実施形態の流れ図が示されている。
【００４９】
図中、ステップＳ１０１はシミュレーション対象となる系の諸元を読み込み、変数の初期化を行う処理である。諸元の読み込みは、入力装置２０３或いは通信装置２０６からの入力データで行っても良いし、予め２次記憶装置２０５にファイルとして格納しておいたデータを読み出して利用することもできる。系の諸元にはｘ、ｆの初期値としてｘ^（０）、ｆ^（０）、境界条件Ｂ、評価関数Ｌ_１、Ｌ_２が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。感度係数の初期値Ｃ_ｊ ^（０）は次式で計算する。
【００５０】
Ｃ_ｊ ^（０）＝ − （ｘ_ｊ ^（０））^Ｔ（ξ（ｆ_ｊ ^（０））^ξ ^−１Ａ_ｊ）ｘ_ｊ ^（０） …（式１６）
ステップＳ１０２では、反復回数ｋを１に初期化する。反復回数ｋは、上述の設計変数の求解工程に対する反復処理に対するものであり、予め定められた終了条件（評価関数Ｌ_１，Ｌ_２の条件）を満たしたら終了する。
【００５１】
ステップＳ１０３では、ラグランジェ未定定数λの更新を行う。ステップＳ１０４では、設計変数ベクトルｆの更新を行う。ステップＳ１０５では、状態変数ベクトルｘの更新を行う。ステップＳ１０６では、感度係数の計算を行う。
【００５２】
上記ステップＳ１０３及びＳ１０４は、拘束条件付最適化問題の解法として採用する方法によって異なる。例えば、この設計変数ベクトルｆ^（ｋ）を算出する求解工程は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかにより実現される。以下、解法として最適性規準法を用いる場合と、逐次凸関数近似法を用いる場合について詳述する。
【００５３】
（最適性規準法）
まず、最適性規準法を用いる場合について説明する。構造最適化問題に対する最適性規準法は（非特許文献４）及び（非特許文献５）等によって公知である。
【００５４】
最適性規準法では、（式４）及び（式５）をラグランジェ未定定数λを用いて、１つの式にまとめる。
【００５５】
Ｌ（ｆ，λ）＝Ｌ_２（ｆ，ｘ） − λ（Ｗ（ｆ） − Ｗ_０） …（式１７）
ステップＳ１０３で実行されるλの更新は、次式を用いて行われる。
【００５６】
λ^{（ｋ＋１）} ＝ｍｉｎ［０，［（１／Ｗ_０）Ｗ（ｆ^（ｋ））］ ^αλ^（ｋ）］ …（式１８）
ただし、上付きのｋ，ｋ＋１は反復回数を表す。また、αはべき乗係数で、通常０．８５程度の値を設定する。
【００５７】
ステップＳ１０４で実行される設計変数の更新は次式による。
【００５８】

ただし、ζは設計変数更新の変動幅の制約値で、０．３程度の値を取る。また、ｓ_ｊ ^（ｋ）は次式で与えられる。
【００５９】
ｓ_ｊ ^（ｋ）＝［（１／λ^{（ｋ−１）}）Ｃ_ｊ ^{（ｋ−１）}］αｆ_ｊ ^{（ｋ−１）} …（式２０）
ここで、Ｃ_ｊ ^（ｋ）は第２の評価関数Ｌ_２の、設計変数ｆ_ｊ ^（ｋ）に関する感度係数と呼ばれる量であり、ステップＳ１０６の処理で計算されるものである。
【００６０】
更に（式１９）で計算された設計変数ｆ_ｊ ^（ｋ）は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、例えば１０^−３程度の値である。
【００６１】
ステップＳ１０５では、状態変数ベクトルｘ^（ｋ）を更新する。
【００６２】
状態変数ベクトルｘ^（ｋ）は、（式４）の評価関数の極小化の過程で得られる。この解法は共役残差法と呼ばれており、後述する。
【００６３】
ステップＳ１０６では、感度係数Ｃ_ｊ ^（ｋ）を次式により計算する。
【００６４】
Ｃ_ｊ ^（ｋ）＝ − （ｘ_ｊ ^（ｋ））^Ｔ（ξ（ｆ_ｊ ^（ｋ））^ξ ^−１Ａ_ｊ）ｘ_ｊ ^（ｋ） …（式２１）
（逐次凸関数近似法）
次に、逐次凸関数近似法を用いた場合の、ステップＳ１０３及びＳ１０４の処理について詳述する。逐次凸関数近似法では設計変数ベクトルの各成分に対してスケーリングを行うことによって、数値計算誤差を避けることができるが、以下では説明を簡単にするため、スケーリングについては記述していない。逐次凸関数近似法の詳細は、（参考文献１）及び（参考文献２）等によって公知である。
【００６５】
（参考文献１）藤井：“パソコンで解く構造デザイン”，丸善（２００２年４月）
（参考文献２）Ｃ．Ｆｌｅｕｒｙ： ”ＣＯＮＬＩＮ，ａｎｅｆｆｉｃｉｅｎｔｄｕａｌｏｐｔｉｍｉｚｅｒｂａｓｅｄｏｎｃｏｎｖｅｘａｐｐｒｏｘｉｍａｔｉｏｎｃｏｎｃｅｐｔｓ，” ＳｔｒｕｃｔｕｒａｌＯｐｔｉｍｉｚａｔｉｏｎ，１，ｐｐ．８１−８９（１９８９）
ステップＳ１０３では、ラグランジェ未定定数を次式により更新する。
【００６６】
λ^（ｋ）＝ λ^{（ｋ−１）} ＋ Δλ^（ｋ） …（式２２）
ただし、初期値λ^（０）は例えば０．０１程度の値を設定する。
【００６７】
Δλ^（ｋ）は次の方程式の解として与えられる。
【００６８】
Δλ^（ｋ）＝２×（１ − Ｗ_０／Ｗ（ｆ^{（ｋ−１）}））λ^{（ｋ−１）} …（式２３）
ステップＳ１０４では、設計変数ベクトルの更新を行う。
【００６９】

ただし、ｇ^（ｋ）＝ −（ｆ_ｊ ^{（ｋ−１）}）^２Ｃ_ｊ ^{（ｋ−１）} …（式２５）
更に、（式２４）で計算された設計変数ｆ_ｊ ^（ｋ）は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、１０^−３程度の値である。
【００７０】
ステップＳ１０５及びＳ１０６の処理は前述のとおりである。
【００７１】
（共役勾配法）
ステップＳ１０５において状態変数ベクトルｘ^（ｋ）を算出する求解工程として、共役勾配法を用いた処理を図４を用いて説明する。
【００７２】
ステップＳ３０１では、シミュレーション対象となる系の諸元を読み込む。系の諸元にはｘ^（ｋ）の初期値ｘ^（０）、設計変数ベクトルｆ^（ｋ）の初期値ｆ^（０）、変位拘束条件Ｂ_１、評価関数Ｌ_１が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。また、残差ベクトルｒ、探索方向ベクトルｐを初期化する。それぞれの初期値をｒ^（０）、ｐ^（０）と表す。
【００７３】
まず、ｒ^（０）は、ｘ^（０）と各ノード（節点）に加わる力をベクトル形式で記述した節点力ベクトルｂとを用いて次式で計算される。
【００７４】
ｒ^（０）＝ｂ − Ａｘ^（０） …（式２６）
探索方向ベクトルｐの初期値ｐ^（０）はｒ^（０）と等しくする。
【００７５】
ステップＳ３０２では前処理を行う。ます、全体剛性行列Ａの対角成分のうち、０となる行あるいは列番号を調べる。０となる成分がなければ何もせずにステップＳ３０３に進む。０となる成分があれば、ｂ（節点力ベクトル）のうち一致する番号の成分を０にする。
【００７６】
ステップＳ３０３では、反復回数ｔを１に初期化する。反復回数ｔは、上述の状態変数の求解工程に対する反復処理に対するものであり、以下のステップＳ３０４からＳ３０８の処理の反復回数が設定値を超えるか、残差ベクトルｒ^（ｔ）のノルムの２乗（ｒ^（ｔ），ｒ^（ｔ））が予め設定された所定値より小さくなったら、反復を終了する。以下、反復回数ｔ回目の値をｘ^（ｔ）のように書く。
【００７７】
ステップＳ３０４では状態ベクトルの更新係数α^（ｔ）を計算する。
【００７８】
α ＝（ｒ^{（ｔ−１）}，ｐ^{（ｔ−１）}）／（ｐ^{（ｔ−１）}，Ａｐ^{（ｔ−１）}） …（式２７）
ステップＳ３０５では、状態変数ベクトルｘを更新する。
【００７９】
ｘ^（ｔ）＝ｘ^{（ｔ−１）} ＋ αｐ^{（ｔ−１）} …（式２８）
ステップＳ３０６では、残差ベクトルｒ^（ｔ）を次式で計算する。
【００８０】
ｒ^（ｔ）＝ｂ − Ａｘ^{（ｔ−１）} …（式２９）
ステップＳ３０７では、探索方向ベクトルｐの更新係数βを計算する。
【００８１】
β ＝（ｒ^（ｔ），ｒ^（ｔ））／（ｒ^{（ｔ−１）}，ｒ^{（ｔ−１）}） …（式３０）
ただし、ｔ＝１のときはβ＝１とする。
【００８２】
ステップＳ３０８では、探索方向ベクトルｐを更新する。
【００８３】
ｐ^（ｔ）＝ｒ^（ｔ）＋ βｐ^{（ｔ−１）} …（式３１）
共役勾配法及びその様々な派生法は、（参考文献３）等によって公知である。
【００８４】
（参考文献３）森，杉原，室田：“線形計算，”岩波講座応用数学，岩波書店（１９９４年）
（本実施形態の構造最適設計の具体例）
本例では、任意の位置に加重を受ける片持ち梁の最適形状自動設計に、本実施形態を適用するものである。説明を簡単にするために平面歪問題に限定する。
【００８５】
図５に示すように、構造部材の存在を可能とする設計領域は長方形４０２であり、有限要素法に従って、該領域を縦ｎ_ｙ、横ｎ_ｘに等間隔に分割する。分割された部分領域をセルと呼び、左下及び右上のセルがそれぞれ（１，１）及び（ｎ_ｙ，ｎ_ｘ）となるように番号付けを行う。
【００８６】
同様に格子点をノードと呼び、左下及び右上のノードが（１，１）及び（ｎ_ｙ＋１，ｎ_ｘ＋１）となるように番号付けを行う。
【００８７】
図中、４０１は支持部材、４０３は荷重ベクトルである。
【００８８】
セル（ｊ，ｋ）には特性関数値ｆ（ｊ，ｋ）が対応する。ここで、特性関数値とはセル（ｊ，ｋ）における構造部材の存在確率を示す０から１の正の実数値をとる変数であり、本実施形態における設計変数ベクトルｆの要素である。
ｆ＝（ｆ（１，１），ｆ（１，２）， …，ｆ（ｎ_ｙ，ｎ_ｘ））^Ｔ …（式３２）
同様にノード（ｊ，ｋ）には横方向変位ｕ（ｊ，ｋ）と縦方向変位ｖ（ｊ，ｋ）が対応する。これらは任意の値を取る実数であり、本実施形態における状態変数ベクトルｘの要素である。
ｘ＝（ｕ（１，１），ｖ（１，１），ｕ（１，２），ｖ（１，２），…，ｕ（ｎ_ｙ＋１，ｎ_ｘ＋１），ｖ（ｎ_ｙ＋１，ｎ_ｘ＋１））^Ｔ…（式３３）
図６に、本具体例に置ける計算結果を示す。図中、黒い領域が構造部材が存在する領域である。本具体施例において、構造部材の設計領域のアスペクト比は縦／横が２対１であり、その解析解は、水平方向に対して±４５度の梁が組み合わせられたものであることが知られている。図６に示した計算結果は、このような解析解とよく一致していることがわかる。
【００８９】
尚、本実施形態では、図１のように、バスライン２０７に各装置が接続した１つのシステムとして説明しているが、これら各要素をネットワークなどを介して接続して複数の装置のシステムで実現することも、あるいは汎用のコンピュータを複数接続して並列に多重処理あるいは分担処理により処理するようにしても良い。
【００９０】
又、本発明は、上記構造最適設計方法を実現するプログラム、及び該プログラムをコンピュータ読み出し可能に記憶する記憶媒体をも含むものである。
【００９１】
【発明の効果】
以上説明したように、本発明によれば、全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量の軽減が実現できる。
【図面の簡単な説明】
【図１】本実施形態の構造最適設計装置の構成例を示す図である。
【図２】図１の１次記憶装置及び２次記憶装置の記憶構成例を示す図である。
【図３】本実施形態の処理手順例を示すフローチャートである。
【図４】図３のステップＳ１０５を更に詳細に示すの処理手順のフローチャートである。
【図５】本実施形態に適応する具体例の問題設定の説明図である。
【図６】図５に係る本実施形態による計算結果を示す図である。[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to an optimal design of a structure, which is a solution to a problem of optimizing design parameters, and more particularly to an automatic design technique for shape optimization including a topology of a structural member.
[0002]
[Prior art]
Optimal design of a structural topology is the problem of determining the optimal topology and geometry of a structural member under given conditions. Hereinafter, the topology and shape dimensions of the structural members are referred to as design variable functions, and the above 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. In the design variable function optimization problem, the state variable function optimization problem must be solved for each design variable function value. In this sense, the optimal design of the structural topology can be regarded 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.
[0003]
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. 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.
[0004]
On the other hand, the following three types of methods have been provided for the optimization problem of the design variable function (see Non-Patent Document 1 or Non-Patent Document 2).
1. Evolutionary method (hereinafter referred to as E method)
2. Homogenization method (hereinafter, H method)
3. Material distribution method (hereinafter, MD method) or density method (hereinafter, D 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. Patent Literature 1 discloses 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 requires a trial calculation based on the accumulated know-how, and therefore solves the problems of the prior art that could not cope with an actual design problem having a large number of design variables by the following method. ing. That is, an approximation optimization calculating device using an approximate expression of discrete design variable data such as a frame member cross-sectional dimension, and a detailed optimization calculating device using the design variable data are provided. It constitutes an optimal design device for the structure.
[0005]
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 partial regions 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 for searching for the local optimum solution is greatly reduced (see Non-Patent Document 3).
[0006]
The MD method or the D method is a method of expressing a topology or a change in shape and dimension of a structural member by assigning a real number in the range of 0 to 1 indicating the existence ratio of the structural member to each element. It is similar to the H method in the sense that sensitivity analysis is enabled by replacing discrete information on whether or not a structural member exists with a continuous value of the existence ratio, but the number of parameters is smaller than that of the H method. It is easy to model and has a wide application range.
[0007]
Non-Patent Document 4 discloses a method for optimizing the phase of a structure using the D method. The features of this method are as follows. Since the voxel finite element method (space is divided at equal intervals) is used, the element stiffness matrices for all elements are 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. By using a combination of the preconditioned conjugate gradient method and the Element-by-Element method to solve a large-scale system of linear equations, the solution can be obtained without assembling the entire rigidity matrix. Less is needed.
[0008]
In the homogenization method, six design variables (in the case of three dimensions) are required 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.
[0009]
[Patent Document 1]
JP-A-11-314631 [Non-Patent Document 1]
S. Bulman, J .; Sienz, E .; Hinton: "Comparisons between algorithms for structural topology optimization optimization using a series of benchmark studies," Computers and Stories, 79. 1203-1218 (2001).
[Non-patent document 2]
YL. Hsu, MS. Hsu, CT. Chen: "Interpreting results from topology optimization using density concentrations," Computers and Structures, 79, pp. 1049-1058 (2001).
[Non-Patent Document 3]
Hiroshi Yamakawa: "Optimization Design," Computational Mechanics and CAE Series 9, Baifukan (1996)
[Non-patent document 4]
Fujii, Suzuki, Ohtsubo: "Topology optimization of structures using voxel finite element method," Transactions of JSCES, Paper No. 20000010 (2000).
[0010]
[Problems to be solved by the invention]
However, the above prior art has the following problems.
[0011]
In general, the structure optimization problem is formulated as a double optimization problem including a state variable vector optimization problem at each iteration step of the design variable vector optimization problem. When the optimization problem of the design variable vector is called the outer optimization and the optimization problem of the state variable vector is called the inner optimization, the inner optimization uses the design variable vector as a parameter, that is, fixes the design variable vector, This is the problem of finding the state variable vector. This is usually called structural analysis and can be solved using a solution of a linear equation by a finite element method.
[0012]
However, when the structure changes and the structural member in a certain region does not exist, for example, when a hole is formed in the member, the design variable of the corresponding element becomes 0, and as a result, the Young's modulus of the element becomes 0. Then, the coefficient matrix of the linear equation does not have a full rank and cannot be solved by the direct method because the inverse matrix cannot be calculated.
[0013]
Therefore, the prior art has dealt with the following when the design variable vector becomes zero.
[0014]
The equations are rearranged so that the coefficient matrix has a full rank, that is, the number of ranks is equal to the number of equations.
[0015]
The value of the design variable vector is replaced with a minute value close to 0 instead of exactly 0.
[0016]
However, in the former method, there is a possibility that the equation is updated every time the design variable vector is updated, and it takes a long calculation time.
[0017]
Further, the latter method, that is, setting the design variable of an element to a minute value corresponds to the presence of a physically thin film or a weak member. That is, since the conventional method does not accurately represent a portion where there is no substance, the accuracy of the obtained calculation result remains questionable.
[0018]
In view of the above-mentioned conventional problems, the present invention can execute calculations without performing special processing even when the overall stiffness matrix becomes singular, thereby simplifying a program and further reducing the amount of calculation. To provide an optimal design method for a structure that reduces the risk of erosion.
[0019]
[Means for Solving the Problems]
In order to solve the above problems, a structure optimization design method according to the present invention includes a first solution step for solving an optimization problem of a first evaluation functional with respect to a state variable vector and a design variable vector; A second solution step for solving a second evaluation functional optimization problem with respect to a design variable vector, the structure optimization design method for solving a structure optimization design problem formulated as a dual optimization problem, A variable vector is a displacement at each node, and the design variable vector is an abundance ratio of a structural member in each element. The first solving step includes the second solving step as one step, and is stored in the first storage means. A design for reading out the stored design variable vector and state variable vector to update the design change vector, and storing the updated design variable vector in the first storage means; A number updating step, wherein the second solving step is configured by a conjugate gradient method, and performs a preprocessing step of performing a preprocessing on the contact force vector based on the overall stiffness matrix; Reading the updated design change vector and state variable vector to update the state variable vector, and storing the updated state variable vector in the second storage means. It is characterized in that it is not initialized at the time when the solution step 2 is started.
[0020]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, a configuration example and an operation example of a structure optimum design apparatus that realizes the structure optimum design method of the present invention will be described.
[0021]
<Example of the configuration of the optimal structure design apparatus of the present embodiment>
FIG. 1 shows an example of the configuration of a structural optimization design apparatus according to the present embodiment.
[0022]
In the figure, 201 is a CPU for calculation and control, 202 is a display device, 203 is an input device, 204 is a primary storage device, 205 is a secondary storage device, 206 is a communication device, and 207 is a bus line. The components in the present embodiment are stored in advance in the secondary storage device 205 as a program, are loaded into the primary storage device by a command input from the input device 203 or the communication device 206, and are executed by the CPU 201. is there.
[0023]
FIG. 2 shows a storage example of data and programs stored in the primary storage device 204 and the secondary storage device 205 in order to realize the structural optimization design method of the present embodiment.
[0024]
The primary storage device 204 is composed of a RAM or a ROM, and has a data storage area 214 and a program storage area 224. FIG. 2 does not show an OS or BIOS program that is resident without being changed, or an OS or BIOS parameter that is not changed. In the data storage area 214, various data for realizing the structure optimum design of the present embodiment is temporarily stored. It is essential that the data used again / the data used repeatedly is temporarily stored. Temporary storage is not essential for data that does not need to be stored and held.
[0025]
For example, the data storage area 214, x (state variable vector) 214a, f (a design variable vector) 214b, (the evaluation function of the state variable vector) _{L 1} 214c, (evaluation function design variable vector) _{L 2} 214d, B1 (Boundary conditions for state variables) 214e, W ₀ (weight conditions for design variables) 214f, X ⁽⁰⁾ and f ⁽⁰⁾ (initial values) 214g, λ (Lagrange undetermined constant) 218h, A ( Element stiffness matrix) 214i, k (number of iterations of changing design variables) 214j, C (sensitivity coefficient) 214k, r (residual vector) 124m, p (search direction vector) 215n, r ⁽⁰⁾ , p ⁽⁰⁾ ( (Initial value) 214p, t (the number of repetitions of state variable change) 214q, b (nodal force vector) 214r, and other data 214s. On the other hand, the program loaded from the secondary storage device 205 or the communication device 206 is stored in the program storage area 224 and executed by the CPU 201.
[0026]
The secondary storage device 205 is, for example, an external large-capacity memory such as a floppy disk or a CD. The secondary storage device 205 also has a data storage area 215 and a program storage area 225.
[0027]
Although not shown in detail, the data storage area 215 stores various initial values 215a, various conditions 215b, or a database 215c for a dual optimal value problem. The database 215c may store and use a standardized matrix or the like.
[0028]
In the program storage area 225, an optimality criterion module 225a, a sequential convex function approximation module 225b, a sequential linear programming module 225c, and a module 225d for other methods are used in the present embodiment as a solution method for design variable vectors. Is stored and used as a solution method in the state variable vector, a conjugate gradient method module 225e, a GCR method module 225f, a GCR (k) method 225g, an Orthomin (k) method 225h, a GMREG (k) method 225i, and other methods. Module 225j is stored. Modules corresponding to the method selected from the programs stored in the program storage area 225 are loaded into the program storage area 224 of the primary storage device 204 and executed by the CPU 201 to realize the structure optimal design method of the present embodiment. .
[0029]
That is, the present structure optimization design apparatus comprises: first solving means for solving an optimization problem of a first evaluation functional for a state variable vector and a design variable vector; and second evaluation for the state variable vector and the design variable vector. A second solving means for solving a functional optimization problem, and a structure optimization design device for solving a structure optimization design problem formulated as a double optimization problem, wherein the state variable vector is The displacement and the design variable vector are the abundance ratio of the structural member in each element, and the second solving means is configured by a conjugate gradient method, and before performing preprocessing on the nodal force vector based on the overall stiffness matrix. A processing means, wherein the state variable vector is not initialized when the processing of the second solving means is started. The preprocessing means sets the contact force vector component corresponding to the row or column where the diagonal component of the overall rigidity matrix is zero to zero. The first solving means executes any one of a sequential linear programming method, an optimality criterion method, and a sequential convex function approximation method.
[0030]
<Formulation of the structural optimization problem of the present embodiment>
For the following description, a structural optimization problem will be formulated.
[0031]
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.
[0032]
The state variable vector x and the design variable vector f are respectively written as column vectors as follows.
[0033]
x = (x ₁ , x ₂ ,..., x _m ) ^T (Equation 1)
f = (f ₁ , f ₂ ,..., f _n ) ^T (Equation 2)
Here, T represents transposition. x is an m-dimensional vector and f is an n-dimensional vector. In the structural optimization problem, the state variable vector is composed of displacement of each node, and the design variable vector is composed of the existence ratio of the structure in each element.
[0034]
For example, in the structure optimization problem in a two-dimensional plane, splitting the design area vertical n _y, the horizontal n _x,
The number n of elements is ( _nx × _ny ), and the number of nodes is ( _nx + 1) × ( _ny + 1).
[0035]
Since the state variable vector is a set of vertical and horizontal displacements of each node, the dimension m of x is 2 × ( _nx + 1) × ( _ny + 1).
[0036]
The boundary condition of x is usually given as a displacement constraint condition. Describe it _{B 1} and.
B ₁ (x) = 0 (Equation 3)
Assuming that the evaluation functions for the state variable vector x and the design variable vector f are L ₁ and L ₂ , respectively, they are defined as follows.
L ₁ : = L ₁ (x, f) = (b−Ax) ^T (b−Ax) (Equation 4)
L ₂ : = L ₂ (f, x) = (１／) × ^T Ax (Equation 5)
Here, A is the overall stiffness matrix, and b is the nodal force vector. Here, b is given in advance, and A is a function of the design variable vector f.
[0037]
The overall stiffness matrix can be configured by superimposing the element stiffness matrix A _{j on} the element j. In the case of the plane distortion problem, a design variable vector, that is, an element rigidity matrix in consideration of the existence ratio f _j of the structural element is shown below. However, a four-node element is adopted as a two-dimensional finite element, and the element displacement vector _xj is configured as follows.
x _j = (u ₁ , v ₁ , u ₂ , v ₂ ,..., u ₄ , v ₄ ) (formula 6)
Here, u _k and v _k are the horizontal and vertical displacements of the node k.
[0038]
The element rigidity matrix corresponding to the element displacement vector of (Equation 6) is given by the following equation.
[0039]
(Equation 1)

[0040]
In the above formula _{a i, j = (λ +} 2μ) <∂ x e i, ∂ x e j> + μ <∂ y e i, ∂ y e j> ... ( Equation 8)
_{b i, j = λ <∂} x e i, ∂ y e j> + μ <∂ x e j, ∂ y e i> ... ( Eq. 9)
_{c i, j = λ <∂} x e j, ∂ x e i> + μ <∂ x e i, ∂ y e j> ... ( Equation 10)
_{d i, j = (λ +} 2μ) <∂ y e i, ∂ y e j> + μ <∂ x e i, ∂ x e j> ... ( Equation 11)
Here, _{e j} is a basis vector for node j, <,> represents an inner product of vectors. _{Ｘ x} is a partial derivative with respect to _x , and ∂ _y is a partial derivative with respect to y. ξ is set to 2 in accordance with the above-mentioned (Non-Patent Document 4). Λ and μ are material constants called lame constants, which are calculated from the Young's modulus E and Poisson's ratio ν by the following equation.
[0041]
λ = Eν / (1 + ν) (1-2ν) (Equation 12)
μ = E / {2 (1 + ν)} (Expression 13)
As is easily understood from (Equation 7), when the existence ratio of a certain element becomes 0, all the components of the element rigidity matrix _Aj become 0. The overall stiffness matrix is formed by superposing element stiffness matrices. Therefore, when the design variable f _{j of the} element including the node j becomes 0, all the components of the row and the column corresponding to the node j of the overall rigidity matrix become 0. As described above, in the structural optimization problem, the overall stiffness matrix is generally a singular matrix.
[0042]
Constraint on the design variables, if the structure optimum design, the upper limit value W ₀ of the total weight.
[0043]
(Equation 2)

[0044]
In addition, the design variables are defined with the constraint conditions for the value range as follows.
[0045]
0 ≦ f _j ≦ 1 j = 1,..., N (Equation 15)
From the above notation, the optimal design problem is formulated as a problem that minimizes (Expression 5) under the constraint conditions (Expression 14) and (Expression 15).
[0046]
The state variable x is obtained as a solution that minimizes (Equation 4) under the constraint condition (Equation 3).
[0047]
<Operation Example of Structural Optimal Design Apparatus of Present Embodiment>
An example of the processing procedure of the present embodiment will be described based on the above formulation.
[0048]
FIG. 3 shows a flowchart of this embodiment.
[0049]
In the figure, step S101 is processing for reading specifications of a system to be simulated and initializing variables. 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 x, ^x as the initial value of ^{f (0), f (0} ), the boundary condition B, and includes an evaluation function _{_L} 1, _L _2. Based on this information, the program secures a necessary variable area in the primary storage device 204 and sets a value. The initial value C _j ⁽⁰⁾ of the sensitivity coefficient is calculated by the following equation.
[0050]
_{^{_{^{C j (0) = - (}}}} x j (0)) T (ξ (f j (0)) ξ -1 A j) x j (0) ... ( Equation 16)
In step S102, the number of repetitions k is initialized to one. The number of iterations k is for the iterative process for the above-described design variable solving process, and ends when a predetermined end condition (conditions of the evaluation functions L ₁ and L ₂ ) is satisfied.
[0051]
In step S103, the Lagrange undetermined constant λ is updated. In step S104, the design variable vector f is updated. In step S105, the state variable vector x is updated. In step S106, a sensitivity coefficient is calculated.
[0052]
The above steps S103 and S104 are different depending on the method adopted as a solution to the optimization problem with constraints. For example, the solution calculation step of calculating the design variable vector f ^(k) is realized by any of a sequential linear programming method, an optimality criterion method, and a sequential convex function approximation method. Hereinafter, the case where the optimality criterion method is used and the case where the successive convex function approximation method is used will be described in detail.
[0053]
(Optimality criteria method)
First, a case where the optimality criterion method is used will be described. The optimality criterion method for the structure optimization problem is known from (Non-Patent Document 4) and (Non-Patent Document 5).
[0054]
In the optimality criterion method, (Equation 4) and (Equation 5) are combined into one equation using Lagrange's undetermined constant λ.
[0055]
L (f, λ) = L ₂ (f, x) −λ (W (f) −W ₀ ) (Equation 17)
The update of λ executed in step S103 is performed using the following equation.
[0056]
^{λ (k + 1) = min} [0, [(1 / W 0) W (f (k))] α λ (k)] ... ( Equation 18)
Here, superscripts k and k + 1 indicate the number of repetitions. Α is a power coefficient, and is usually set to a value of about 0.85.
[0057]
The update of the design variables executed in step S104 is based on the following equation.
[0058]

Here, ζ is a constraint value of the fluctuation range of the design variable update, and takes a value of about 0.3. Further, s _j ^(k) is given by the following equation.
[0059]
s _j ^(k) = [(1 / λ ^(k−1) ) C _j ^(k−1) ] αf _j ^(k−1) (Equation 20)
Here, C _j ^(k) is an amount called a sensitivity coefficient of the second evaluation function L ₂ with respect to the design variable f _j ^(k) , which is calculated in the process of step S106.
[0060]
Further, the design variable f _j ^(k) calculated by (Equation 19) may be forcibly set to 0 when it becomes smaller than a preset value. The preset value is, for example, a value of about 10 ⁻³ .
[0061]
In step S105, the state variable vector x ^(k) is updated.
[0062]
The state variable vector x ^(k) is obtained in the process of minimizing the evaluation function of (Equation 4). This solution is called a conjugate residual method and will be described later.
[0063]
In step S106, the sensitivity coefficient C _j ^(k) is calculated by the following equation.
[0064]
_{^{_{^{C j (k) = - (}}}} x j (k)) T (ξ (f j (k)) ξ -1 A j) x j (k) ... ( Equation 21)
(Iterative convex function approximation method)
Next, the processing of steps S103 and S104 when the successive convex function approximation method is used will be described in detail. In the successive convex function approximation method, a numerical calculation error can be avoided by scaling each component of the design variable vector, but the scaling is not described below for the sake of simplicity. Details of the successive convex function approximation method are known in (Reference Document 1) and (Reference Document 2).
[0065]
(Reference 1) Fujii: "Structural design solved with a personal computer", Maruzen (April 2002)
(Reference 2) C.I. Fleury: "CONLIN, an effective dual optimizer based on convex application concepts," Structural Optimization, 1, pp. 81-89 (1989)
In step S103, the Lagrange undetermined constant is updated by the following equation.
[0066]
λ ^(k) = λ ^(k−1) + Δλ ^(k) (Expression 22)
However, the initial value λ ⁽⁰⁾ is set to, for example, a value of about 0.01.
[0067]
Δλ ^(k) is given as a solution of the following equation:
[0068]
^{Δλ (k) = 2 × (} 1 - W 0 / W (f (k-1))) λ (k-1) ... ( Equation 23)
In step S104, the design variable vector is updated.
[0069]

^{_{^{However, g (k) = - (}}} f j (k-1)) 2 C j (k-1) ... ( Equation 25)
Further, the design variable f _j ^(k) calculated by (Equation 24) may be forcibly set to 0 when it becomes smaller than a preset value. This preset value is a value of about 10 ⁻³ .
[0070]
The processing in steps S105 and S106 is as described above.
[0071]
(Conjugate gradient method)
A process using the conjugate gradient method as a solution step for calculating the state variable vector x ^(k) in step S105 will be described with reference to FIG.
[0072]
In step S301, the specifications of the system to be simulated are read. Initial value ^x of ^{x (k)} is the specifications of the system ^(0), the initial value ^f of the design variable vector ^{f (k)} ^(0), the displacement constraints _{B 1,} which includes an evaluation function _{L 1.} Based on this information, the program secures a necessary variable area in the primary storage device 204 and sets a value. Further, the residual vector r and the search direction vector p are initialized. The respective initial values are represented by r ⁽⁰⁾ and p ⁽⁰⁾ .
[0073]
First, r ⁽⁰⁾ is calculated by the following equation using x ⁽⁰⁾ and a nodal force vector b that describes a force applied to each node (node) in a vector format.
[0074]
r ⁽⁰⁾ = b-Ax ⁽⁰⁾ (Expression 26)
The initial value p ⁽⁰⁾ of the search direction vector p is equal to r ⁽⁰⁾ .
[0075]
In step S302, pre-processing is performed. First, among the diagonal components of the overall stiffness matrix A, a row or column number that becomes 0 is checked. If there is no component that becomes 0, the process proceeds to step S303 without doing anything. If there is a component which becomes 0, the component of the matching number in b (nodal force vector) is set to 0.
[0076]
In step S303, the number of repetitions t is initialized to one. The number of repetitions t is for the repetition processing for the above-described state variable solving process, and the number of repetitions of the processing of the following steps S304 to S308 exceeds a set value or the square of the norm of the residual vector r ^(t). When (r ^(t) , r ^(t) ) becomes smaller than a predetermined value, the repetition is terminated. Hereinafter, the value of the t-th iteration is written as x ^(t) .
[0077]
In step S304, the update coefficient α ^(t) of the state vector is calculated.
[0078]
α = (r ^(t−1) , p ^(t−1) ) / (p ^(t−1) , Ap ^(t−1) ) (Expression 27)
In step S305, the state variable vector x is updated.
[0079]
x ^(t) = x ^(t−1) + αp ^(t−1) (Expression 28)
In step S306, the residual vector r ^(t) is calculated by the following equation.
[0080]
r ^(t) = b-Ax ^(t-1) (Equation 29)
In step S307, the update coefficient β of the search direction vector p is calculated.
[0081]
β = (r ^(t) , r ^(t) ) / (r ^(t-1) , r ^(t-1) ) (Equation 30)
However, when t = 1, β = 1.
[0082]
In step S308, the search direction vector p is updated.
[0083]
p ^(t) = r ^(t) + βp ^(t−1) (formula 31)
The conjugate gradient method and its various derivatives are known from (Reference Document 3) and the like.
[0084]
(Reference 3) Mori, Sugihara, Murota: "Linear Computation," Iwanami Course Applied Mathematics, Iwanami Shoten (1994)
(Specific Example of Structural Optimization Design of the Present Embodiment)
In the present embodiment, the present embodiment is applied to automatic design of the optimum shape of a cantilever that receives a load at an arbitrary position. In order to simplify the explanation, the problem is limited to the plane distortion problem.
[0085]
As shown in FIG. 5, the design area to allow the presence of the structural member is rectangular 402, according to 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, performing numbered as lower left and upper right of the cell is respectively (1,1) and (n _{y, n} _x).
[0086]
Similarly called grid points and node performs numbered as lower left and upper right of the node becomes (1,1) and _{_{(n y + 1, n x}} +1).
[0087]
In the figure, reference numeral 401 denotes a support member, and 403 denotes a load vector.
[0088]
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 probability of the structural member in the cell (j, k), and is an element of the design variable vector f in the present embodiment.
f = (f (1,1), f (1,2),..., f ( _ny , _nx )) ^T (Formula 32)
Similarly, a horizontal displacement u (j, k) and a vertical displacement v (j, k) correspond to the node (j, k). These are real numbers having arbitrary values, and are elements of the state variable vector x in the present embodiment.
x = (u (1,1), v (1,1), u (1,2), v (1,2), ..., u (n y + 1, n x +1), v (n y +1, _nx +1)) ^T (Equation 33)
FIG. 6 shows a calculation result in this specific example. In the figure, a black area is an area where a structural member exists. In this specific example, the aspect ratio of the design area of the structural member is 2: 1 in the vertical / horizontal direction, and the analytical solution thereof is known to be a combination of ± 45 ° beams in the horizontal direction. Have been. It can be seen that the calculation results shown in FIG. 6 are in good agreement with such an analytical solution.
[0089]
In this embodiment, as shown in FIG. 1, a single system in which each device is connected to the bus line 207 is described. However, these components are connected via a network or the like, and a system of a plurality of devices is used. Alternatively, a plurality of general-purpose computers may be connected and processed in parallel by multiplex processing or shared processing.
[0090]
The present invention also includes a program for realizing the above-described structure optimization design method, and a storage medium for storing the program in a computer-readable manner.
[0091]
【The invention's effect】
As described above, according to the present invention, even when the overall stiffness matrix becomes singular, calculations can be performed without performing special processing, thereby simplifying the program and further reducing the amount of calculation. Can be reduced.
[Brief description of the drawings]
FIG. 1 is a diagram illustrating a configuration example of a structural optimization design apparatus according to an embodiment.
FIG. 2 is a diagram illustrating an example of a storage configuration of a primary storage device and a secondary storage device of FIG. 1;
FIG. 3 is a flowchart illustrating an example of a processing procedure according to the embodiment;
FIG. 4 is a flowchart of a processing procedure showing step S105 of FIG. 3 in further detail;
FIG. 5 is an explanatory diagram of a problem setting of a specific example applied to the present embodiment.
FIG. 6 is a diagram showing a calculation result according to the embodiment shown in FIG. 5;

Claims

A first solving step for solving an optimization problem of a first evaluation functional for the state variable vector and the design variable vector; and a second solving step for solving an optimization problem of a second evaluation functional for the state variable vector and the design variable vector. A structural optimization design method for solving a structural optimization design problem formulated as a dual optimization problem having two solution steps.
The state variable vector is the displacement at each node, the design variable vector is the existence rate of the structural member in each element,
The first solving step includes the second solving step as one step, reads out the design variable vector and the state variable vector stored in the first storage means, updates the design change vector, and updates the updated design change vector. A design variable updating step of storing a design variable vector in the first storage means,
The second solving step is configured by a conjugate gradient method, and performs a preprocessing step of performing a preprocessing on the contact force vector based on the overall stiffness matrix, and a design change vector stored in the second storage unit. Reading a state variable vector, updating the state variable vector, and storing the updated state variable vector in the second storage means, wherein the state variable vector starts the second solution step. A structural optimization design method characterized in that it is not initialized at the time of execution.