JP3861259B2

JP3861259B2 - Phase optimization system and phase optimization method

Info

Publication number: JP3861259B2
Application number: JP2000404449A
Authority: JP
Inventors: 惠三石井; 昇菊池; 茂青村
Original assignee: 株式会社くいんと
Priority date: 2000-12-19
Filing date: 2000-12-19
Publication date: 2006-12-20
Anticipated expiration: 2020-12-19
Also published as: JP2002189760A

Description

【０００１】
【発明の属する技術分野】
本発明は、均質化法および有限要素法を用いた最適構造設計に使用されるシステムおよび方法に関する。
【０００２】
【従来の技術】
位相最適化の方法は、１９８８年デンマーク工科大学のマーティン・フィリップ・ベンソー（ＭａｒｔｉｎＰｈｉｌｉｐＢｅｎｄｓｏｅ）とミシガン大学の菊池昇によって提示され（ＧｅｎｅｒａｔｉｎｇＯｐｔｉｍａｌＴｏｐｏｌｏｇｉｅｓｉｎＳｔｒｕｃｔｕｒａｌＤｅｓｉｇｎｕｓｉｎｇａＨｏｍｏｇｅｎｉｚａｔｉｏｎＭｅｔｈｏｄ．ＭａｒｔｉｎＰｈｉｌｉｐＢｅｎｄｓｏｅａｎｄＮｏｂｏｒｕＫｉｋｕｃｈｉ、ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓｉｎＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓａｎｄＥｎｇｉｎｅｅｒｉｎｇ７１（１９８８）１９７−２２４．）、鈴木克幸と菊池昇によって２次元問題の完成を見た（Ａｈｏｍｏｇｅｎｉｚａｔｉｏｎｍｅｔｈｏｄｆｏｒｓｈａｐｅａｎｄｔｏｐｏｌｏｇｙｏｐｔｉｍｉｚａｔｉｏｎ．ＫａｔｓｕｙｕｋｉＳｕｚｕｋｉａｎｄＮｏｂｏｒｕＫｉｋｕｃｈｉ、ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓｉｎＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓａｎｄＥｎｇｉｎｅｅｒｉｎｇ９３（１９９１）２９１−３１８．）。
【０００３】
２次元問題に関するこの方法は、現在では３次元問題に拡張されている。すなわち、連続体の位相最適化方法の概要は、まず、構造物の設計対象領域を、有限要素と称する小さな立方体または直方体に分割したモデルを生成する。したがって、構造物は、有限個の要素に分割される。この有限要素はさらに、図２に示すように、立方体または直方体の空洞を有している。したがって、設計対象領域は、微視構造、すなわち有限要素が、周期的に多数並んだ多孔質体のモデルに仮定できる。次に、このモデルにおいて、各有限要素からなる微視構造の空洞のサイズをコントロールすることにより、構造物の体積を制約して、最も剛性の高い位相形状を算出する。
【０００４】
【発明が解決しようとする課題】
しかしながら、従来の３次元の位相最適設計方法は、各有限要素で残された材料を充填された状態との比で密度表示すると、微視構造で算出した位相形状に中間密度が多く見られ、結果の判断が簡単ではない。このため、より鮮明な形状を生成する位相最適設計方法が望まれていた。
【０００５】
本発明は、位相最適化の結果に現れる中間密度をより少なくして、形状の認識がよりしやすい位相最適化システムおよび位相最適化方法を提供することを目的とする。
【０００６】
【課題を解決するための手段】
本発明による位相最適化システムは、構造物が周期性を有する微視構造により構成されると仮定し、前記微視構造からなる所定の有限要素から前記構造物の最適位相形状を生成する位相最適化システムであって、前記構造物の有限要素モデルデータおよび位相最適化制御データを入力する入力手段と、前記微視構造は格子状の該微視構造の外形を構成する外辺と内形を構成する内辺とからなり、異なる長さの前記内辺または同一の長さの前記内辺を有する複数の各前記微視構造が存在し、異なる長さの前記内辺または同一の長さの前記内辺を有する複数の各前記微視構造に均質化法を適用して巨視平均的な応力−歪マトリックスデータを作成するデータ作成手段と、複数の代表的な前記内辺の長さを用いて作成した前記応力−歪マトリックスデータに各前記内辺の長さに応じた補間を成し、各要素の主軸方向に座標変換をした応力−歪マトリックスを算出する第１の演算手段と、第１の演算手段により算出された前記応力−歪マトリックスを用いて有限要素解析を行い、目的関数および該目的関数の感度計算を行う第２の演算手段と、前記有限要素解析により算出された前記目的関数と、あらかじめ定められた閾値との関係を評価し収束判定する収束判定手段と、収束判定の結果に応じて最適化規準法で各前記内辺の長さを更新する更新手段と、前記データ作成手段、前記第１の演算手段、前記第２の演算手段、前記収束判定手段および前記更新手段による処理を前記収束判定手段により収束判定がなされるまで繰り返す手段とを具備し、前記最適位相形状を生成することを特徴とする。
【０００９】
【発明の実施の形態】
本発明による最適位相形状の生成について、図面を参照して詳細に説明する。
【００１０】
図３を参照して、本発明による最適位相形状の生成手順を説明する。
【００１１】
まず、設計対象となる構造物、すなわち橋梁や船舶などの建造物や、自動車部材などの機械部品が形成される領域に適用する有限要素データを作成する（ステップＳ１１）。本説明において、直方体を有限個に分割したものを有限要素という。有限要素データには、節点座標情報、要素結合情報、要素と材料定数の対応情報、要素と物理定数の対応情報、材料定数情報、物理定数情報、拘束条件情報、荷重条件情報が含まれる。この有限要素データを作成するには、米国オハイオ州所在のＳＤＲＣ（ＳｔｒｕｃｔｕｒａｌＤｙｎａｍｉｃｓＲｅｓｅａｒｃｈＣｏｒｐｏｒａｔｉｏｎ）製のＦＥＭＡＰなど、市販の有限要素モデル作成ソフトウェアが利用できる。具体的には、次の１）から７）の手順で行う。
【００１２】
１）設計領域を設定し、この３次元領域を複数の六面体などの多面体ブロックに分割する。本実施例では、この設計領域内で、最適化された構造物が形成されることになる。たとえば、構造物の生成が予想される領域を覆域にして、接合された数個の六面体ブロックを想定し、これを設定領域に設定する。同様に四面体または五面体のブロックを想定してもよい。比較的単純な形状の構造物であれば、設計領域を分割せず、１つの六面体ブロックを設計領域として設定することも可能である。
【００１３】
２）分割された六面体ブロックごとに、六面体ブロックの各頂点の座標値（ｘ，ｙ，ｚ）を設定する。
【００１４】
３）この六面体ブロックのいずれかの頂点を原点に設定した局所座標系を導入して、全体座標系から局所座標系へ座標変換する。さらに、各座標軸方向の分割数を設定する。
【００１５】
４）設定した分割数に従ってこの六面体ブロックを分割し、局所座標空間で立方体または直方体の要素を生成する。生成する情報は、ｘｙｚの各軸ごとの、分割された頂点（以下、節点と呼ぶ）の番号（１から始まる連続番号）と座標値（ｘ，ｙ，ｚ）、および生成された複数の立方体または直方体を構成する頂点の番号列（要素結合情報）である。
【００１６】
５）各要素の各座標を、局所座標系から全体座標系に座標変換する。
【００１７】
６）各要素に材料特性を割当てる。
【００１８】
分割された他の六面体ブロックについても、同様に２）から６）の工程を繰り返す。すべての六面体ブロックの要素を合わせると、設計領域がすべて要素によって規定でき、設計領域が有限要素モデル化されたことになる。
【００１９】
７）生成された有限要素モデルに、境界条件、すなわち拘束条件および荷重を設定する。具体的には、拘束条件は拘束する節点に対してのみ、決められた各自由度ごとに、拘束するかどうかを入力する。荷重は、集中、分布にかかわらずプログラムの中で載荷する節点の値に置き換えたものを自由度分与える。
【００２０】
次に、最適化制御データ、すなわち、要素毎の設計／非設計フラッグ、制約体積、繰り返し計算を始めるステップ番号、繰り返し計算を終了させるステップ番号、収束判定に使う閾値、設計変数の変動許容値（ｍｏｖｉｎｇｂｏｕｎｄ）の入力をする（ステップＳ１２）。
【００２１】
図１に示す通り、本実施例で設定する要素の微視構造は格子状であって、外郭を形成する外辺と内郭を形成する内辺とを有する。すなわちこの微視構造は、外辺長さに対して内辺長さが最小であるとき、充填された立方体となり、外辺長さに対して内辺長さが最大であるとき、空洞化される。また、外辺長さに対して内辺長さがこれらの中間の値をとるとき、六面体状の要素の各辺が、格子状となって残るので、フレームベースユニットセルと称する。
【００２２】
外辺の長さを１に正規化し、内辺の代表長さを０．０、０．２、０．４、０．６、０．８、１．０の６種類に設定した場合について説明するが、他の値を設定して微視構造を生成してもよい。したがって、以下の演算では、６種類の内辺の代表長さを組み合わせた合計２１６種類の微視構造形状について計算を行い、さらに、それを補間することによって、代表長さの中間値を得ている。
【００２３】
各微視構造の外辺は、図１に示すように、各微視構造ごとのｘ、ｙ、ｚの各軸ごとに、１の長さを有する。また、各微視構造の内辺は、各微視構造ごとのｘ、ｙ、ｚの各軸ごとに、長さａ、ｂ、ｃを有する。各微視構造において、各内辺の長さ（ａ，ｂ，ｃ）は一般に、それぞれ異なる値を有する。この場合、この微視構造は、各面に長方形状の孔を有する形状になる。また、微視構造内の内辺長さ（ａ，ｂ，ｃ）は、各微視構造ごとに、同一の長さに設定して、すなわち、この微視構造が正方形状の孔を有すると設定して演算を行うこともできる。この場合、各微視構造ごとに３種ある長さ（ａ，ｂ，ｃ）に代わって、各軸に係る長さを（ａ，ａ，ａ）に設定できるので、演算において扱うパラメータが少なくなり、入力データ量および演算に要する負荷がより小さくなる。
【００２４】
次に、この微視構造に均質化法を適用して、巨視平均的な応力−歪マトリックスデータを作成する（ステップＳ２１）。具体的には、次の１）から３）の手順で行う。
【００２５】
１）微視構造内で成り立つ次式を満足する特性変位（Ｃｈａｒａｃｔｅｒｉｓｔｉｃｄｅｆｏｒｍａｔｉｏｎ） χ∈Ｖ_Ｙを求める。特性変位とは、微視構造の工学歪６成分に対応した基本変形モードのことである。
【００２６】
【数１】

【００２７】
ここで、Ｅ_ｉｊｋｌは弾性テンソル、ｘはマクロ座標、ｙはマクロ座標のある場所の近傍を構成する微視構造を形成する座標系を示す。Ｙは微視構造の領域を示し、νは微視構造が持つ仮想変位で、領域Ｙ内の周期関数（Ｙ−ｐｅｒｉｏｄｉｃ）とする。
【００２８】
２）式（１）を解いて得られた特性変位χを次式に代入して、巨視的に均質化された弾性テンソルを得る。
【００２９】
【数２】

【００３０】
３）テンソルの対称性を利用し、工学的な応力−歪みマトリックスを得る。
【００３１】
【数３】

【００３２】
つぎに、内辺の長さに応じた巨視平均的な応力−歪マトリックスおよび巨視平均的な応力−歪マトリックスの感度をベジエ補間によって算出する（ステップＳ３１）。ベジエ補間のほか、ラグランジェ補間などの補間によって算出することもできる。感度とは、内辺の長さが微小に変化した場合の応力−歪マトリックスの変化率である。内辺の長さが前記２１６種類に合わない中間値になる場合は５次のベジエ関数を補間近似関数として用い補間する。ベジエ関数は以下の通りである。
【００３３】
【数４】

【００３４】
ここで、ｍは使用するベジエ関数の次数（本実施例では５を使用）、Ｑｉは制御点の値，ｘは０．０から１．０の間の内辺の長さを示す。また、感度を計算するためのベジエ関数の微分はＰ（ｘ）をｘで微分して作成する。
【００３５】
後述の収束判定において、収束が認められないときは、ステップＳ３１以降の計算が繰り返されるが、繰返し計算２回目以降はさらに、有限要素解析により求まった各有限要素中心で計算した主応力により主軸を決定し、ベジエ補間を成した応力−歪マトリックスに主軸に合わせる座標変換を行い巨視平均的な応力−歪マトリックスを算出する（ステップＳ３２）。
【００３６】
すなわち、前述のベジエ補間により算出した応力−歪マトリックスを［ＤＨ０］とし、有限要素解析で求まる各要素の主応力により決まる座標変換マトリックスを［Ｔ］とすれば、回転した座標系における応力−歪マトリックス［ＤＨ］は次のように計算する。
【００３７】
【数５】

【００３８】
巨視平均的な応力−歪マトリックスの感度も同様に求める。
【００３９】
ステップＳ３１（繰り返し２回目以降はステップＳ３１およびＳ３２）で算出された各有限要素の巨視平均的な応力−歪マトリックスを用いて対象モデルを有限要素解析し、目的関数である平均コンプライアンス（総歪みエネルギーの２倍の値）と、その感度を、次のように計算する（ステップＳ４０）。
【００４０】
有限要素モデルは、節点の番号と座標値、要素を構成する節点番号のテーブル、モデルの材料定数、物理定数、拘束条件、節点に付加される荷重値、等によって構成され、解かれるべき節点の変位（Ｕｘ，Ｕｙ，Ｕｚ）と要素内歪の関係を次のように定義する。
【００４１】
【数６】

【００４２】
ここで、｛ｕ｝は仮想変位ベクトル、｛ε｝は要素内の歪ベクトル、｜Ｂ｜は歪−変位関係マトリックスである。また要素内における仮想変位によって起こる歪と応力の関係は次式で決定される。
【００４３】
【数７】

【００４４】
ここで、｜Ｄ｜は応力−歪関係マトリックスと呼ばれ、線形弾性問題では等方性材料の場合は１種類のヤング率とポアソン比のみで決定できる。本実施例では後述する均質化法を用いて計算した内辺のサイズに応じた巨視平均的な応力−歪マトリックスを用いることとする。｛σ｝は要素内応力成分である。
【００４５】
仮想仕事の原理を適用して、ポテンシャルエネルギーを停留させる条件から次の関係式が導かれる。
【００４６】
【数８】

【００４７】
ここでＫは剛性マトリックスと呼ばれる。ｆは外荷重ベクトルとする。この方程式を解いて得られた節点変位ベクトルｕから前述した関係に従って歪および応力を計算し、これらの値から目的関数である次式の平均コンプライアンスを計算する。
【００４８】
【数９】

【００４９】
ここでΩは要素内領域を示す。目的関数の感度は同様にして、前記の手続きで計算した巨視平均的な応力−歪マトリックスの感度を使い次式のように計算する。
【００５０】
【数１０】

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

【００５４】
（１１）式が満足されるか、最適化制御データとして入力した打ち切り回数に早く到達したほうを採用する。
【００５５】
収束せず、かつ、計算打ち切り回数に達しない場合は、最適化規準法により異なる各内辺または同一内辺の長さを更新し（ステップＳ６０）、ステップ３１、３２、４０、５０および６０を、収束判定が収束となるか、計算打ち切り回数に達するまで繰り返す。最適化規準法は次のような手順で行う。
【００５６】
１）ラグランジュ未定乗数Λを導入してラグランジュアンを次のように定義する。
【００５７】
【数１２】

【００５８】
ここで、ａは設計変数（内辺の長さ：実際は３種類または１種類）、Πは目的関数（ここでは平均コンプライアンス）、Ｖ（ａ）は内辺のサイズによって変化した設計領域全体の体積、Ｖｃは最適化制御データとして入力した制約体積を示す。
【００５９】
２）このラグランジュアンが停留することが最適化の必要条件になるから、（１２）式の変分を取ると、
【００６０】
【数１３】

【００６１】
となる。δａ、δΛの任意性から次式が成り立つ。
【００６２】
【数１４】

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

【００６５】
【数１６】

【００６６】
ここでηは適当な減速パラメータで、数値実験により０．７５としたが、他の適切な値を設定してもよい。ｋは繰り返し回数を示す。
【００６７】
上式Ｅ_ｉ ^（ｋ）、Ｗ_ｉ ^（ｋ）は繰り返し計算の中で常に１になるように設計変数（内辺の長さ）、ラグランジュ未定乗数を更新する。但し、算出した設計変数は上限値が１．０、下限値は０．０と決められているので、たとえば１．２など、どちらかに達した場合は、上下限値を採用する。さらに、設計変数の変動幅が大きくなると収束が不安定になることから、変動幅にも最適化制御データで入力をした値で制約（ｍｏｖｉｎｇｂｏｕｎｄ）を与える。
【００６８】
【実施例】
図４ないし１２を参照して、２つのモデルに、従来例による位相最適化方法および本発明による位相最適化方法を適用した計算例について説明する。
【００６９】
図４に、短い片持ち梁モデルを有限要素によって分割した例を示す。片側全面を固定し、反対側の中心部に下向きの荷重を載荷した短い片持ち梁モデルについて、制約体積を設計領域の４０％として計算した位相を図５、図６および図７に示す。
【００７０】
計算モデルは、寸法：０．２×０．１１×０．１５［ｍ］、ヤング率：９８０００［Ｍｐａ］、ポアソン比：０．３、総節点数：４０３０、総要素数：３３００、荷重は９８［Ｎ］とした。
【００７１】
図５を参照すると、従来の方法による計算例では、表面に中間密度を残すものの、濃淡は見てとれる。これはユニットセルに直方体の穴を空け、主応力方向に穴の向きを回転することで非等方性材料が形成されるが、結果的にはこの操作が中間密度の減少に寄与していると思われる。これに対して、図６を参照すると、本発明による計算例では、同じ繰り返し数にもかかわらずほとんどの要素で材料は充填状態と空隙状態に分離している。
【００７２】
両手法共に材料マトリックスを主軸に合わせて回転しているので，結果の差は剪断剛性の強さに起因すると思われる。
【００７３】
図１２に従来例および本実施例による目的関数の履歴を示す。両手法を比較すると、ユニットセル形状の違いによる収束性の差はほとんど認められなかった。
【００７４】
図８に、直方体の設計領域（１／２対称モデル）の下から約１／３の水平帯状部分に分布荷重を載荷した、橋梁をイメージしたモデルについて位相を計算した例を示す。支持条件は橋脚ができるであろう下部の１箇所と橋桁の端部を完全拘束し，対称面に対称条件を付加した。
【００７５】
制約条件は設計領域の１０％の体積とし、分布荷重を載荷した部分を含む要素２層分は始めから材料を充填状態にして、計算中も一切変更をしない非設計領域とした。
【００７６】
計算モデルはミニチュアモデル、すなわち縮小した計算モデルである。このモデルの寸法は０．３３×０．０６×０．２３［ｍ］、ヤング率は２０５８００［Ｍｐａ］、ポアソン比は０．３、総節点数は５７１２、総要素数は５５４、また、荷重は単位面積あたり９．８［Ｍｐａ］である。
【００７７】
図９、１０および１１に、このモデルの計算例を示す。この計算例は、前述した短い片持ち梁の計算例に比べて、さらに大きな効果が現れている。従来の方法は、設計領域から９０％の材料を削り取ったにもかかわらず、鉄橋を構成する主要な部材を判断するに足る情報を算出したが、かなり曖昧な部分が残った。本発明の方法はさらに鮮明な位相を算出した。剪断剛性の差がこのモデルではさらに大きな影響を与えたと推測できる。両手法とも我々が日常目撃する鉄道橋全体のイメージは再現している。
【００７８】
【発明の効果】
実施例でも解るように、本発明によれば、図２に示す従来のユニットセルを使用した場合に比べ、位相最適化の結果に現れる中間密度の部分を大幅に減らし、認識しやすい形状を得ることができる。
【図面の簡単な説明】
【図１】本発明によるフレームベースユニットセルを示した図である。
【図２】従来のユニットセルを示した図である。
【図３】位相最適化方法の概要を示したフローチャートである。
【図４】短い片持ち梁の有限要素モデルを示した図である。
【図５】従来の均質化法による位相を示した図である。
【図６】本実施例による位相を示した図である。
【図７】本実施例による位相を示した平面図である。
【図８】橋をイメージした有限要素モデルを示した図である。
【図９】従来の均質化法による位相を示した図である。
【図１０】本実施例による位相を示した図である。
【図１１】本実施例による位相を示した図である。
【図１２】目的関数の履歴を示したグラフである。
【符号の説明】
ａ、ｂ、ｃ内辺の長さ[0001]
BACKGROUND OF THE INVENTION
The present invention relates to systems and methods used for optimal structural design using homogenization and finite element methods.
[0002]
[Prior art]
The method of phase optimization was presented in 1988 by Martin Philip Bendsoe of the Danish Institute of Technology and Noboru Kikuchi of the University of Michigan. Computer Methods in Applied Mechanics and Engineering 71 (1988) 197-224., Katsuyuki Suzuki and Noboru Kikuchi saw the completion of a two-dimensional problem (A homogenization method for topologies). Katsyuyuki Suzuki and Noboru Kikuchi, Computer Methods in Applied Mechanics and Engineering 93 (1991) 291-318.).
[0003]
This method for two-dimensional problems is now extended to three-dimensional problems. That is, in the outline of the phase optimization method for a continuum, first, a model in which the design target area of the structure is divided into small cubes or cuboids called finite elements is generated. Therefore, the structure is divided into a finite number of elements. The finite element further has a cubic or cuboid cavity, as shown in FIG. Therefore, the design target region can be assumed to be a microscopic structure, that is, a model of a porous body in which a large number of finite elements are periodically arranged. Next, in this model, by controlling the size of the cavity of the microscopic structure composed of each finite element, the volume of the structure is constrained and the most rigid phase shape is calculated.
[0004]
[Problems to be solved by the invention]
However, in the conventional three-dimensional phase optimum design method, when the density is displayed in a ratio with the state in which the material remaining in each finite element is filled, there are many intermediate densities in the phase shape calculated by the microscopic structure, Judging results is not easy. For this reason, a phase optimum design method for generating a clearer shape has been desired.
[0005]
It is an object of the present invention to provide a phase optimization system and a phase optimization method in which the intermediate density appearing in the result of phase optimization is reduced and the shape can be easily recognized.
[0006]
[Means for Solving the Problems]
The phase optimization system according to the present invention assumes that the structure is constituted by a microscopic structure having periodicity, and generates a phase optimum for the optimal phase shape of the structure from a predetermined finite element composed of the microscopic structure. An input means for inputting finite element model data and phase optimization control data of the structure, and the microscopic structure has an outer shape and an inner shape constituting the outer shape of the microscopic structure in a lattice shape. A plurality of the microscopic structures having different inner lengths or the same inner side, and different inner lengths or the same length. Data creation means for creating macroscopic average stress-strain matrix data by applying a homogenization method to each of the plurality of microscopic structures having the inner side, and using a plurality of representative lengths of the inner side The stress-strain matrix created by And a first computing means for calculating a stress-strain matrix obtained by performing interpolation according to the length of each inner side in the data and performing coordinate transformation in the principal axis direction of each element. The finite element analysis is performed using the stress-strain matrix, the second calculation means for calculating the objective function and the sensitivity of the objective function, the objective function calculated by the finite element analysis, and a predetermined function Convergence determining means that evaluates the relationship with the threshold and determines convergence; updating means that updates the length of each inner edge by the optimization criterion method according to the result of convergence determination; the data creating means; calculation means, said second arithmetic means, and wherein the processing by the convergence determining means and the updating means and means for repeated until convergence determination is made by the convergence determining unit generates the optimum phase shape That.
[0009]
DETAILED DESCRIPTION OF THE INVENTION
The generation of the optimum phase shape according to the present invention will be described in detail with reference to the drawings.
[0010]
With reference to FIG. 3, the procedure for generating the optimum phase shape according to the present invention will be described.
[0011]
First, finite element data to be applied to a structure to be designed, that is, a region where a mechanical part such as a building such as a bridge or a ship or an automobile member is formed (step S11). In this description, a cuboid divided into a finite number is called a finite element. The finite element data includes node coordinate information, element coupling information, element-material constant correspondence information, element-physical constant correspondence information, material constant information, physical constant information, constraint condition information, and load condition information. In order to create the finite element data, commercially available finite element model creation software such as FEMAP manufactured by SDRC (Structural Dynamics Research Corporation) in Ohio, USA can be used. Specifically, the following steps 1) to 7) are performed.
[0012]
1) A design area is set, and this three-dimensional area is divided into a plurality of polyhedral blocks such as a hexahedron. In this embodiment, an optimized structure is formed within this design area. For example, assuming a region where the structure is expected to be generated and covering several hexahedral blocks that are joined, this is set as the setting region. Similarly, tetrahedral or pentahedral blocks may be assumed. If the structure has a relatively simple shape, it is possible to set one hexahedral block as the design region without dividing the design region.
[0013]
2) For each divided hexahedron block, the coordinate value (x, y, z) of each vertex of the hexahedron block is set.
[0014]
3) A local coordinate system in which one of the vertices of this hexahedral block is set as the origin is introduced, and coordinate conversion is performed from the global coordinate system to the local coordinate system. Furthermore, the number of divisions in each coordinate axis direction is set.
[0015]
4) Divide the hexahedron block according to the set number of divisions, and generate cube or cuboid elements in the local coordinate space. The information to be generated includes the number (vertical number starting from 1) and coordinate values (x, y, z) of the divided vertices (hereinafter referred to as nodes) for each axis of xyz, and a plurality of generated cubes Or it is the number sequence (element connection information) of the vertices constituting the rectangular parallelepiped.
[0016]
5) The coordinates of each element are transformed from the local coordinate system to the global coordinate system.
[0017]
6) Assign material properties to each element.
[0018]
Similarly, the steps 2) to 6) are repeated for the other divided hexahedron blocks. When the elements of all hexahedral blocks are combined, the design area can be defined by the elements, and the design area is modeled as a finite element.
[0019]
7) Set boundary conditions, that is, constraint conditions and loads, in the generated finite element model. Specifically, whether or not to restrain is input for each of the determined degrees of freedom only for the nodes to be restrained. Regardless of concentration and distribution, the load is replaced with the value of the node loaded in the program, giving the degree of freedom.
[0020]
Next, the optimization control data, that is, the design / non-design flag for each element, the constraint volume, the step number for starting the iterative calculation, the step number for ending the iterative calculation, the threshold used for convergence determination, and the design variable fluctuation tolerance ( (moving bound) is input (step S12).
[0021]
As shown in FIG. 1, the microscopic structure of the elements set in the present embodiment is a lattice shape, and has an outer side that forms an outer shell and an inner side that forms an inner shell. In other words, this microscopic structure becomes a filled cube when the inner side length is minimum with respect to the outer side length, and is hollowed out when the inner side length is maximum with respect to the outer side length. The Further, when the inner side length takes an intermediate value with respect to the outer side length, each side of the hexahedral element remains in a lattice shape, and hence is referred to as a frame base unit cell.
[0022]
The case where the length of the outer side is normalized to 1 and the representative length of the inner side is set to six types of 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0 will be described. However, the microscopic structure may be generated by setting other values. Therefore, in the following calculation, calculation is performed for a total of 216 types of microscopic structure shapes that combine the representative lengths of the six types of inner sides, and further, an intermediate value of the representative length is obtained by interpolation. Yes.
[0023]
As shown in FIG. 1, the outer side of each microscopic structure has a length of 1 for each axis of x, y, z for each microscopic structure. In addition, the inner side of each microscopic structure has lengths a, b, and c for each axis of x, y, z for each microscopic structure. In each microscopic structure, the length (a, b, c) of each inner side generally has a different value. In this case, the microscopic structure has a shape having a rectangular hole on each surface. Further, the inner side length (a, b, c) in the microscopic structure is set to the same length for each microscopic structure, that is, when the microscopic structure has a square hole. It is also possible to set and perform calculations. In this case, instead of the three types of lengths (a, b, c) for each microscopic structure, the lengths related to the respective axes can be set to (a, a, a). Thus, the amount of input data and the load required for calculation are further reduced.
[0024]
Next, a homogenization method is applied to the microscopic structure to create macroscopic average stress-strain matrix data (step S21). Specifically, the following steps 1) to 3) are performed.
[0025]
1) Characteristic deformation (Characteristic deformation) χ∈V _Y that satisfies the following equation established in the microscopic structure is obtained. Characteristic displacement is a basic deformation mode corresponding to six engineering strain components of a microscopic structure.
[0026]
[Expression 1]

[0027]
Here, E _ijkl is an elastic tensor, x is a macro coordinate, and y is a coordinate system that forms a microscopic structure that forms the vicinity of a place where the macro coordinate exists. Y indicates a region of the microscopic structure, and ν is a virtual displacement of the microscopic structure, which is a periodic function (Y-periodic) in the region Y.
[0028]
2) The characteristic displacement χ obtained by solving equation (1) is substituted into the following equation to obtain a macroscopically homogenized elastic tensor.
[0029]
[Expression 2]

[0030]
3) An engineering stress-strain matrix is obtained using the symmetry of the tensor.
[0031]
[Equation 3]

[0032]
Next, the sensitivity of the macroscopic average stress-strain matrix and the macroscopic average stress-strain matrix corresponding to the length of the inner side is calculated by Bezier interpolation (step S31). In addition to Bezier interpolation, it can also be calculated by interpolation such as Lagrange interpolation. Sensitivity is the rate of change of the stress-strain matrix when the length of the inner side changes slightly. When the inner side length is an intermediate value that does not match the 216 types, interpolation is performed using a fifth-order Bezier function as an interpolation approximation function. The Bezier function is as follows.
[0033]
[Expression 4]

[0034]
Here, m is the order of the Bezier function to be used (5 is used in this embodiment), Qi is the value of the control point, and x is the length of the inner side between 0.0 and 1.0. The differential of the Bezier function for calculating the sensitivity is created by differentiating P (x) by x.
[0035]
In the convergence determination described later, when convergence is not recognized, the calculation after step S31 is repeated, but after the second iteration, the principal axis is further calculated by the principal stress calculated at the center of each finite element obtained by finite element analysis. A macroscopic average stress-strain matrix is calculated by performing coordinate transformation in accordance with the principal axis of the stress-strain matrix determined and subjected to Bezier interpolation (step S32).
[0036]
That is, if the stress-strain matrix calculated by Bezier interpolation is [DH0] and the coordinate transformation matrix determined by the principal stress of each element obtained by finite element analysis is [T], the stress-strain in the rotated coordinate system is The matrix [DH] is calculated as follows.
[0037]
[Equation 5]

[0038]
The sensitivity of the macroscopic average stress-strain matrix is determined in the same manner.
[0039]
The target model is subjected to finite element analysis using the macroscopic average stress-strain matrix of each finite element calculated in step S31 (steps S31 and S32 after the second iteration), and the average compliance (total strain energy) that is an objective function And the sensitivity thereof are calculated as follows (step S40).
[0040]
A finite element model is composed of node numbers and coordinate values, a table of node numbers that make up the elements, model material constants, physical constants, constraint conditions, load values added to the nodes, etc. The relationship between displacement (Ux, Uy, Uz) and intra-element distortion is defined as follows.
[0041]
[Formula 6]

[0042]
Here, {u} is a virtual displacement vector, {ε} is a strain vector in the element, and | B | is a strain-displacement relationship matrix. The relationship between strain and stress caused by virtual displacement in the element is determined by the following equation.
[0043]
[Expression 7]

[0044]
Here, | D | is called a stress-strain relationship matrix. In the case of an isotropic material, | D | can be determined by only one kind of Young's modulus and Poisson's ratio. In this embodiment, a macroscopic average stress-strain matrix corresponding to the size of the inner side calculated using the homogenization method described later is used. {Σ} is an in-element stress component.
[0045]
Applying the principle of virtual work, the following relational expression is derived from the condition for stopping potential energy.
[0046]
[Equation 8]

[0047]
Here, K is called a stiffness matrix. f is an external load vector. Strain and stress are calculated from the nodal displacement vector u obtained by solving this equation according to the relationship described above, and the average compliance of the following equation, which is an objective function, is calculated from these values.
[0048]
[Equation 9]

[0049]
Here, Ω indicates the element internal region. Similarly, the sensitivity of the objective function is calculated as follows using the sensitivity of the macroscopic average stress-strain matrix calculated in the above procedure.
[0050]
[Expression 10]

[0051]
Convergence is determined by evaluating the relationship between the objective function calculated by the finite element analysis described above and a predetermined threshold value (step S50). Specifically, this is performed as follows.
[0052]
Suppose that the objective function (average compliance) calculated at the nth time is Π _n , the objective function calculated at the (n−1) th time is Π _n−1, and the convergence judgment threshold value input as optimization control data is ε. The convergence determination is performed as follows.
[0053]
[Expression 11]

[0054]
The one that satisfies the expression (11) or that has reached the number of times of censoring input as optimization control data earlier is adopted.
[0055]
If it does not converge and does not reach the number of calculation censoring, the different internal sides or the length of the same internal side is updated by the optimization criterion method (step S60), and steps 31, 32, 40, 50 and 60 are performed. This is repeated until the convergence judgment is converged or the number of calculation aborts is reached. The optimization criteria method is as follows.
[0056]
1) Lagrangian multiplier Λ is introduced to define Lagrangian as follows.
[0057]
[Expression 12]

[0058]
Here, a is a design variable (inner side length: actually three or one type), Π is an objective function (here, average compliance), and V (a) is a volume of the entire design region that changes depending on the size of the inner side. , Vc represents the constraint volume input as optimization control data.
[0059]
2) Since stopping this Lagrangian is a necessary condition for optimization, taking the variation of equation (12),
[0060]
[Formula 13]

[0061]
It becomes. The following expression holds from the arbitraryness of δa and δΛ.
[0062]
[Expression 14]

[0063]
Using this, the following update rule is applied to the design variable a _i and Lagrange undetermined multiplier Λ.
[0064]
[Expression 15]

[0065]
[Expression 16]

[0066]
Here, η is an appropriate deceleration parameter and is set to 0.75 by a numerical experiment, but other appropriate values may be set. k indicates the number of repetitions.
[0067]
The design variables (inner side length) and Lagrange undetermined multipliers are updated so that the above formulas E _i ^(k) and W _i ^(k) are always 1 in the iterative calculation. However, since the upper limit value of the calculated design variable is determined to be 1.0 and the lower limit value is determined to be 0.0, the upper and lower limit values are adopted when the value reaches either 1.2, for example. Furthermore, since the convergence becomes unstable when the fluctuation range of the design variable becomes large, a constraint (moving bound) is given to the fluctuation range by the value input by the optimization control data.
[0068]
【Example】
A calculation example in which the phase optimization method according to the conventional example and the phase optimization method according to the present invention are applied to the two models will be described with reference to FIGS.
[0069]
FIG. 4 shows an example in which a short cantilever model is divided by finite elements. FIG. 5, FIG. 6 and FIG. 7 show the phases calculated for a short cantilever model in which the entire surface of one side is fixed and a downward load is loaded on the opposite center, assuming that the constrained volume is 40% of the design area.
[0070]
The calculation model is: dimensions: 0.2 × 0.11 × 0.15 [m], Young's modulus: 98000 [Mpa], Poisson's ratio: 0.3, total number of nodes: 4030, total number of elements: 3300, load is 98 [N].
[0071]
Referring to FIG. 5, in the calculation example according to the conventional method, although the intermediate density is left on the surface, the shading can be seen. This is because the unit cell is made with a rectangular parallelepiped hole and the direction of the hole is rotated in the principal stress direction to form an anisotropic material. As a result, this operation contributes to the reduction of the intermediate density. I think that the. On the other hand, referring to FIG. 6, in the calculation example according to the present invention, the material is separated into a filled state and a void state in most elements despite the same number of repetitions.
[0072]
Since both methods rotate the material matrix along the main axis, the difference in the results seems to be due to the strength of the shear stiffness.
[0073]
FIG. 12 shows the history of objective functions according to the conventional example and this embodiment. When both methods were compared, there was almost no difference in convergence due to the difference in unit cell shape.
[0074]
FIG. 8 shows an example in which the phase is calculated for a model simulating a bridge in which a distributed load is loaded on a horizontal strip portion of about ３ from the bottom of a rectangular parallelepiped design region (1/2 symmetric model). As for the support condition, one part of the lower part where the pier would be able to be formed and the end of the bridge girder were completely constrained, and the symmetry condition was added to the symmetry plane.
[0075]
The constraint condition was a volume of 10% of the design area, and the two layers of elements including the portion where the distributed load was loaded were filled with the material from the beginning, and a non-design area that was not changed during the calculation.
[0076]
The calculation model is a miniature model, that is, a reduced calculation model. The size of this model is 0.33 × 0.06 × 0.23 [m], Young's modulus is 205800 [Mpa], Poisson's ratio is 0.3, the total number of nodes is 5712, the total number of elements is 554, and the load Is 9.8 [Mpa] per unit area.
[0077]
9, 10 and 11 show calculation examples of this model. This calculation example has a greater effect than the calculation example of the short cantilever described above. Although the conventional method calculated 90% of material from the design area, it calculated enough information to determine the main members that make up the steel bridge. The method of the present invention calculated a clearer phase. It can be inferred that the difference in shear stiffness had a greater effect in this model. Both methods reproduce the image of the entire railway bridge we witness everyday.
[0078]
【The invention's effect】
As can be seen from the embodiments, according to the present invention, compared with the case where the conventional unit cell shown in FIG. 2 is used, the intermediate density portion appearing in the phase optimization result is significantly reduced, and an easily recognizable shape is obtained. be able to.
[Brief description of the drawings]
FIG. 1 shows a frame base unit cell according to the present invention.
FIG. 2 is a diagram illustrating a conventional unit cell.
FIG. 3 is a flowchart showing an outline of a phase optimization method.
FIG. 4 shows a finite element model of a short cantilever beam.
FIG. 5 is a diagram showing a phase obtained by a conventional homogenization method.
FIG. 6 is a diagram showing a phase according to the present embodiment.
FIG. 7 is a plan view showing a phase according to the present embodiment.
FIG. 8 is a diagram showing a finite element model in which a bridge is imaged.
FIG. 9 is a diagram showing a phase obtained by a conventional homogenization method.
FIG. 10 is a diagram showing a phase according to the present embodiment.
FIG. 11 is a diagram showing a phase according to the present embodiment.
FIG. 12 is a graph showing a history of objective functions.
[Explanation of symbols]
a, b, c inner side length

Claims

A phase optimization system that generates an optimal phase shape of the structure from a predetermined finite element composed of the microscopic structure on the assumption that the structure is composed of a microscopic structure having periodicity. The finite element model data and the phase optimization control data are input means, and the microscopic structure is composed of an outer side constituting the outer shape of the lattice-like microscopic structure and an inner side constituting the inner shape. There are a plurality of the microscopic structures having the inner side of the length or the inner side of the same length, and the plurality of the microscopic structures having the inner side of different lengths or the inner side of the same length. Data creation means for creating macroscopic average stress-strain matrix data by applying a homogenization method to a microscopic structure, and the stress-strain matrix data created using a plurality of representative inner side lengths In accordance with the length of each inner side. And performing a finite element analysis using the first calculation means for calculating a stress-strain matrix coordinate-transformed in the principal axis direction of each element and the stress-strain matrix calculated by the first calculation means Second function means for calculating sensitivity of the objective function and the objective function, convergence determination means for evaluating convergence by evaluating a relationship between the objective function calculated by the finite element analysis and a predetermined threshold value; Updating means for updating the length of each inner side according to an optimization criterion method according to a result of convergence determination, the data creation means, the first calculation means, the second calculation means, and the convergence determination means And a means for repeating the processing by the updating means until the convergence judgment is made by the convergence judgment means, and generating the optimum phase shape.