(57)【特許請求の範囲】
【請求項１】
構造物が周期性を有する微視構造により構成されると仮定し、前記微視構造からなる所定の有限要素から前記構造物の最適位相形状を生成する位相最適化システムであって、前記構造物の有限要素モデルデータおよび位相最適化制御データを入力する入力手段と、前記微視構造は格子状の該微視構造の外形を構成する外辺と内形を構成する内辺とからなり、異なる長さの前記内辺または同一の長さの前記内辺を有する複数の各前記微視構造が存在し、異なる長さの前記内辺または同一の長さの前記内辺を有する複数の各前記微視構造に均質化法を適用して巨視平均的な応力−歪マトリックスデータを作成するデータ作成手段と、複数の代表的な前記内辺の長さを用いて作成した前記応力−歪マトリックスデータに各前記内辺の長さに応じた補間を成し、各要素の主軸方向に座標変換をした応力−歪マトリックスを算出する第１の演算手段と、第１の演算手段により算出された前記応力−歪マトリックスを用いて有限要素解析を行い、目的関数および該目的関数の感度計算を行う第２の演算手段と、前記有限要素解析により算出された前記目的関数と、あらかじめ定められた閾値との関係を評価し収束判定する収束判定手段と、収束判定の結果に応じて最適化規準法で各前記内辺の長さを更新する更新手段と、前記データ作成手段、前記第１の演算手段、前記第２の演算手段、前記収束判定手段および前記更新手段による処理を前記収束判定手段により収束判定がなされるまで繰り返す手段とを具備し、前記最適位相形状を生成することを特徴とする位相最適化システム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、均質化法および有限要素法を用いた最適構造設計に使用されるシステムおよび方法に関する。
【０００２】
【従来の技術】
位相最適化の方法は、１９８８年デンマーク工科大学のマーティン・フィリップ・ベンソー（Ｍａｒｔｉｎ  Ｐｈｉｌｉｐ  Ｂｅｎｄｓｏｅ）とミシガン大学の菊池昇によって提示され（Ｇｅｎｅｒａｔｉｎｇ  Ｏｐｔｉｍａｌ  Ｔｏｐｏｌｏｇｉｅｓ  ｉｎ  Ｓｔｒｕｃｔｕｒａｌ  Ｄｅｓｉｇｎ  ｕｓｉｎｇ  ａ  Ｈｏｍｏｇｅｎｉｚａｔｉｏｎ  Ｍｅｔｈｏｄ．  Ｍａｒｔｉｎ  ＰｈｉｌｉｐＢｅｎｄｓｏｅ  ａｎｄ  Ｎｏｂｏｒｕ  Ｋｉｋｕｃｈｉ、ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓ  ｉｎ  Ａｐｐｌｉｅｄ  Ｍｅｃｈａｎｉｃｓ  ａｎｄ  Ｅｎｇｉｎｅｅｒｉｎｇ７１（１９８８）１９７−２２４．）、鈴木克幸と菊池昇によって２次元問題の完成を見た（Ａ  ｈｏｍｏｇｅｎｉｚａｔｉｏｎ  ｍｅｔｈｏｄ  ｆｏｒ  ｓｈａｐｅ  ａｎｄ  ｔｏｐｏｌｏｇｙ  ｏｐｔｉｍｉｚａｔｉｏｎ．  Ｋａｔｓｕｙｕｋｉ  Ｓｕｚｕｋｉ  ａｎｄ  Ｎｏｂｏｒｕ  Ｋｉｋｕｃｈｉ、Ｃｏｍｐｕｔｅｒ  Ｍｅｔｈｏｄｓ  ｉｎ  Ａｐｐｌｉｅｄ  Ｍｅｃｈａｎｉｃｓ  ａｎｄ  Ｅｎｇｉｎｅｅｒｉｎｇ９３（１９９１）２９１−３１８．）。
【０００３】
２次元問題に関するこの方法は、現在では３次元問題に拡張されている。すなわち、連続体の位相最適化方法の概要は、まず、構造物の設計対象領域を、有限要素と称する小さな立方体または直方体に分割したモデルを生成する。したがって、構造物は、有限個の要素に分割される。この有限要素はさらに、図２に示すように、立方体または直方体の空洞を有している。したがって、設計対象領域は、微視構造、すなわち有限要素が、周期的に多数並んだ多孔質体のモデルに仮定できる。次に、このモデルにおいて、各有限要素からなる微視構造の空洞のサイズをコントロールすることにより、構造物の体積を制約して、最も剛性の高い位相形状を算出する。
【０００４】
【発明が解決しようとする課題】
しかしながら、従来の３次元の位相最適設計方法は、各有限要素で残された材料を充填された状態との比で密度表示すると、微視構造で算出した位相形状に中間密度が多く見られ、結果の判断が簡単ではない。このため、より鮮明な形状を生成する位相最適設計方法が望まれていた。
【０００５】
本発明は、位相最適化の結果に現れる中間密度をより少なくして、形状の認識がよりしやすい位相最適化システムおよび位相最適化方法を提供することを目的とする。
【０００６】
【課題を解決するための手段】
本発明による位相最適化システムは、構造物が周期性を有する微視構造により構成されると仮定し、前記微視構造からなる所定の有限要素から前記構造物の最適位相形状を生成する位相最適化システムであって、前記構造物の有限要素モデルデータおよび位相最適化制御データを入力する入力手段と、前記微視構造は格子状の該微視構造の外形を構成する外辺と内形を構成する内辺とからなり、異なる長さの前記内辺または同一の長さの前記内辺を有する複数の各前記微視構造が存在し、異なる長さの前記内辺または同一の長さの前記内辺を有する複数の各前記微視構造に均質化法を適用して巨視平均的な応力−歪マトリックスデータを作成するデータ作成手段と、複数の代表的な前記内辺の長さを用いて作成した前記応力−歪マトリックスデータに各前記内辺の長さに応じた補間を成し、各要素の主軸方向に座標変換をした応力−歪マトリックスを算出する第１の演算手段と、第１の演算手段により算出された前記応力−歪マトリックスを用いて有限要素解析を行い、目的関数および該目的関数の感度計算を行う第２の演算手段と、前記有限要素解析により算出された前記目的関数と、あらかじめ定められた閾値との関係を評価し収束判定する収束判定手段と、収束判定の結果に応じて最適化規準法で各前記内辺の長さを更新する更新手段と、前記データ作成手段、前記第１の演算手段、前記第２の演算手段、前記収束判定手段および前記更新手段による処理を前記収束判定手段により収束判定がなされるまで繰り返す手段とを具備し、前記最適位相形状を生成することを特徴とする。
【０００９】
【発明の実施の形態】
本発明による最適位相形状の生成について、図面を参照して詳細に説明する。
【００１０】
図３を参照して、本発明による最適位相形状の生成手順を説明する。
【００１１】
まず、設計対象となる構造物、すなわち橋梁や船舶などの建造物や、自動車部材などの機械部品が形成される領域に適用する有限要素データを作成する（ステップＳ１１）。本説明において、直方体を有限個に分割したものを有限要素という。有限要素データには、節点座標情報、要素結合情報、要素と材料定数の対応情報、要素と物理定数の対応情報、材料定数情報、物理定数情報、拘束条件情報、荷重条件情報が含まれる。この有限要素データを作成するには、米国オハイオ州所在のＳＤＲＣ（Ｓｔｒｕｃｔｕｒａｌ  Ｄｙｎａｍｉｃｓ  Ｒｅｓｅａｒｃｈ  Ｃｏｒｐｏｒａｔｉｏｎ）製のＦＥＭＡＰなど、市販の有限要素モデル作成ソフトウェアが利用できる。具体的には、次の１）から７）の手順で行う。
【００１２】
１）設計領域を設定し、この３次元領域を複数の六面体などの多面体ブロックに分割する。本実施例では、この設計領域内で、最適化された構造物が形成されることになる。たとえば、構造物の生成が予想される領域を覆域にして、接合された数個の六面体ブロックを想定し、これを設定領域に設定する。同様に四面体または五面体のブロックを想定してもよい。比較的単純な形状の構造物であれば、設計領域を分割せず、１つの六面体ブロックを設計領域として設定することも可能である。
【００１３】
２）分割された六面体ブロックごとに、六面体ブロックの各頂点の座標値（ｘ，ｙ，ｚ）を設定する。
【００１４】
３）この六面体ブロックのいずれかの頂点を原点に設定した局所座標系を導入して、全体座標系から局所座標系へ座標変換する。さらに、各座標軸方向の分割数を設定する。
【００１５】
４）設定した分割数に従ってこの六面体ブロックを分割し、局所座標空間で立方体または直方体の要素を生成する。生成する情報は、ｘｙｚの各軸ごとの、分割された頂点（以下、節点と呼ぶ）の番号（１から始まる連続番号）と座標値（ｘ，ｙ，ｚ）、および生成された複数の立方体または直方体を構成する頂点の番号列（要素結合情報）である。
【００１６】
５）各要素の各座標を、局所座標系から全体座標系に座標変換する。
【００１７】
６）各要素に材料特性を割当てる。
【００１８】
分割された他の六面体ブロックについても、同様に２）から６）の工程を繰り返す。すべての六面体ブロックの要素を合わせると、設計領域がすべて要素によって規定でき、設計領域が有限要素モデル化されたことになる。
【００１９】
７）生成された有限要素モデルに、境界条件、すなわち拘束条件および荷重を設定する。具体的には、拘束条件は拘束する節点に対してのみ、決められた各自由度ごとに、拘束するかどうかを入力する。荷重は、集中、分布にかかわらずプログラムの中で載荷する節点の値に置き換えたものを自由度分与える。
【００２０】
次に、最適化制御データ、すなわち、要素毎の設計／非設計フラッグ、制約体積、繰り返し計算を始めるステップ番号、繰り返し計算を終了させるステップ番号、収束判定に使う閾値、設計変数の変動許容値（ｍｏｖｉｎｇ  ｂｏｕｎｄ）の入力をする（ステップＳ１２）。
【００２１】
図１に示す通り、本実施例で設定する要素の微視構造は格子状であって、外郭を形成する外辺と内郭を形成する内辺とを有する。すなわちこの微視構造は、外辺長さに対して内辺長さが最小であるとき、充填された立方体となり、外辺長さに対して内辺長さが最大であるとき、空洞化される。また、外辺長さに対して内辺長さがこれらの中間の値をとるとき、六面体状の要素の各辺が、格子状となって残るので、フレームベースユニットセルと称する。
【００２２】
外辺の長さを１に正規化し、内辺の代表長さを０．０、０．２、０．４、０．６、０．８、１．０の６種類に設定した場合について説明するが、他の値を設定して微視構造を生成してもよい。したがって、以下の演算では、６種類の内辺の代表長さを組み合わせた合計２１６種類の微視構造形状について計算を行い、さらに、それを補間することによって、代表長さの中間値を得ている。
【００２３】
各微視構造の外辺は、図１に示すように、各微視構造ごとのｘ、ｙ、ｚの各軸ごとに、１の長さを有する。また、各微視構造の内辺は、各微視構造ごとのｘ、ｙ、ｚの各軸ごとに、長さａ、ｂ、ｃを有する。各微視構造において、各内辺の長さ（ａ，ｂ，ｃ）は一般に、それぞれ異なる値を有する。この場合、この微視構造は、各面に長方形状の孔を有する形状になる。また、微視構造内の内辺長さ（ａ，ｂ，ｃ）は、各微視構造ごとに、同一の長さに設定して、すなわち、この微視構造が正方形状の孔を有すると設定して演算を行うこともできる。この場合、各微視構造ごとに３種ある長さ（ａ，ｂ，ｃ）に代わって、各軸に係る長さを（ａ，ａ，ａ）に設定できるので、演算において扱うパラメータが少なくなり、入力データ量および演算に要する負荷がより小さくなる。
【００２４】
次に、この微視構造に均質化法を適用して、巨視平均的な応力−歪マトリックスデータを作成する（ステップＳ２１）。具体的には、次の１）から３）の手順で行う。
【００２５】
１）微視構造内で成り立つ次式を満足する特性変位（Ｃｈａｒａｃｔｅｒｉｓｔｉｃ  ｄｅｆｏｒｍａｔｉｏｎ）  χ∈ＶＹを求める。特性変位とは、微視構造の工学歪６成分に対応した基本変形モードのことである。
【００２６】
【数１】
0003861259.tif 000002
【００２７】
ここで、Ｅｉｊｋｌは弾性テンソル、ｘはマクロ座標、ｙはマクロ座標のある場所の近傍を構成する微視構造を形成する座標系を示す。Ｙは微視構造の領域を示し、νは微視構造が持つ仮想変位で、領域Ｙ内の周期関数（Ｙ−ｐｅｒｉｏｄｉｃ）とする。
【００２８】
２）式（１）を解いて得られた特性変位χを次式に代入して、巨視的に均質化された弾性テンソルを得る。
【００２９】
【数２】
0003861259.tif 000003
【００３０】
３）テンソルの対称性を利用し、工学的な応力−歪みマトリックスを得る。
【００３１】
【数３】
0003861259.tif 000004
【００３２】
つぎに、内辺の長さに応じた巨視平均的な応力−歪マトリックスおよび巨視平均的な応力−歪マトリックスの感度をベジエ補間によって算出する（ステップＳ３１）。ベジエ補間のほか、ラグランジェ補間などの補間によって算出することもできる。感度とは、内辺の長さが微小に変化した場合の応力−歪マトリックスの変化率である。内辺の長さが前記２１６種類に合わない中間値になる場合は５次のベジエ関数を補間近似関数として用い補間する。ベジエ関数は以下の通りである。
【００３３】
【数４】
0003861259.tif 000005
【００３４】
ここで、ｍは使用するベジエ関数の次数（本実施例では５を使用）、Ｑｉは制御点の値，ｘは０．０から１．０の間の内辺の長さを示す。また、感度を計算するためのベジエ関数の微分はＰ（ｘ）をｘで微分して作成する。
【００３５】
後述の収束判定において、収束が認められないときは、ステップＳ３１以降の計算が繰り返されるが、繰返し計算２回目以降はさらに、有限要素解析により求まった各有限要素中心で計算した主応力により主軸を決定し、ベジエ補間を成した応力−歪マトリックスに主軸に合わせる座標変換を行い巨視平均的な応力−歪マトリックスを算出する（ステップＳ３２）。
【００３６】
すなわち、前述のベジエ補間により算出した応力−歪マトリックスを［ＤＨ０］とし、有限要素解析で求まる各要素の主応力により決まる座標変換マトリックスを［Ｔ］とすれば、回転した座標系における応力−歪マトリックス［ＤＨ］は次のように計算する。
【００３７】
【数５】
0003861259.tif 000006
【００３８】
巨視平均的な応力−歪マトリックスの感度も同様に求める。
【００３９】
ステップＳ３１（繰り返し２回目以降はステップＳ３１およびＳ３２）で算出された各有限要素の巨視平均的な応力−歪マトリックスを用いて対象モデルを有限要素解析し、目的関数である平均コンプライアンス（総歪みエネルギーの２倍の値）と、その感度を、次のように計算する（ステップＳ４０）。
【００４０】
有限要素モデルは、節点の番号と座標値、要素を構成する節点番号のテーブル、モデルの材料定数、物理定数、拘束条件、節点に付加される荷重値、等によって構成され、解かれるべき節点の変位（Ｕｘ，Ｕｙ，Ｕｚ）と要素内歪の関係を次のように定義する。
【００４１】
【数６】
0003861259.tif 000007
【００４２】
ここで、｛ｕ｝は仮想変位ベクトル、｛ε｝は要素内の歪ベクトル、｜Ｂ｜は歪−変位関係マトリックスである。また要素内における仮想変位によって起こる歪と応力の関係は次式で決定される。
【００４３】
【数７】
0003861259.tif 000008
【００４４】
ここで、｜Ｄ｜は応力−歪関係マトリックスと呼ばれ、線形弾性問題では等方性材料の場合は１種類のヤング率とポアソン比のみで決定できる。本実施例では後述する均質化法を用いて計算した内辺のサイズに応じた巨視平均的な応力−歪マトリックスを用いることとする。｛σ｝は要素内応力成分である。
【００４５】
仮想仕事の原理を適用して、ポテンシャルエネルギーを停留させる条件から次の関係式が導かれる。
【００４６】
【数８】
0003861259.tif 000009
【００４７】
ここでＫは剛性マトリックスと呼ばれる。ｆは外荷重ベクトルとする。この方程式を解いて得られた節点変位ベクトルｕから前述した関係に従って歪および応力を計算し、これらの値から目的関数である次式の平均コンプライアンスを計算する。
【００４８】
【数９】
0003861259.tif 000010
【００４９】
ここでΩは要素内領域を示す。目的関数の感度は同様にして、前記の手続きで計算した巨視平均的な応力−歪マトリックスの感度を使い次式のように計算する。
【００５０】
【数１０】
0003861259.tif 000011
【００５１】
上述した有限要素解析により算出された目的関数とあらかじめ定められた閾値との関係を評価し収束判定する（ステップＳ５０）。具体的には、以下のように行う。
【００５２】
第ｎ回目に計算された目的関数（平均コンプライアンス）をΠｎ、第ｎ−１回目に計算された目的関数をΠｎ−１とし、最適化制御データとして入力した収束判定の閾値をεとすると、収束判定は次のように行う。
【００５３】
【数１１】
0003861259.tif 000012
【００５４】
（１１）式が満足されるか、最適化制御データとして入力した打ち切り回数に早く到達したほうを採用する。
【００５５】
収束せず、かつ、計算打ち切り回数に達しない場合は、最適化規準法により異なる各内辺または同一内辺の長さを更新し（ステップＳ６０）、ステップ３１、３２、４０、５０および６０を、収束判定が収束となるか、計算打ち切り回数に達するまで繰り返す。最適化規準法は次のような手順で行う。
【００５６】
１）ラグランジュ未定乗数Λを導入してラグランジュアンを次のように定義する。
【００５７】
【数１２】
0003861259.tif 000013
【００５８】
ここで、ａは設計変数（内辺の長さ：実際は３種類または１種類）、Πは目的関数（ここでは平均コンプライアンス）、Ｖ（ａ）は内辺のサイズによって変化した設計領域全体の体積、Ｖｃは最適化制御データとして入力した制約体積を示す。
【００５９】
２）このラグランジュアンが停留することが最適化の必要条件になるから、（１２）式の変分を取ると、
【００６０】
【数１３】
0003861259.tif 000014
【００６１】
となる。δａ、δΛの任意性から次式が成り立つ。
【００６２】
【数１４】
0003861259.tif 000015
【００６３】
これを利用して、設計変数ａｉおよびラグランジュ未定乗数Λに関し以下の更新ルールを適用する。
【００６４】
【数１５】
0003861259.tif 000016
【００６５】
【数１６】
0003861259.tif 000017
【００６６】
ここでηは適当な減速パラメータで、数値実験により０．７５としたが、他の適切な値を設定してもよい。ｋは繰り返し回数を示す。
【００６７】
上式Ｅｉ（ｋ）、Ｗｉ（ｋ）は繰り返し計算の中で常に１になるように設計変数（内辺の長さ）、ラグランジュ未定乗数を更新する。但し、算出した設計変数は上限値が１．０、下限値は０．０と決められているので、たとえば１．２など、どちらかに達した場合は、上下限値を採用する。さらに、設計変数の変動幅が大きくなると収束が不安定になることから、変動幅にも最適化制御データで入力をした値で制約（ｍｏｖｉｎｇ  ｂｏｕｎｄ）を与える。
【００６８】
【実施例】
図４ないし１２を参照して、２つのモデルに、従来例による位相最適化方法および本発明による位相最適化方法を適用した計算例について説明する。
【００６９】
図４に、短い片持ち梁モデルを有限要素によって分割した例を示す。片側全面を固定し、反対側の中心部に下向きの荷重を載荷した短い片持ち梁モデルについて、制約体積を設計領域の４０％として計算した位相を図５、図６および図７に示す。
【００７０】
計算モデルは、寸法：０．２×０．１１×０．１５［ｍ］、ヤング率：９８０００［Ｍｐａ］、ポアソン比：０．３、総節点数：４０３０、総要素数：３３００、荷重は９８［Ｎ］とした。
【００７１】
図５を参照すると、従来の方法による計算例では、表面に中間密度を残すものの、濃淡は見てとれる。これはユニットセルに直方体の穴を空け、主応力方向に穴の向きを回転することで非等方性材料が形成されるが、結果的にはこの操作が中間密度の減少に寄与していると思われる。これに対して、図６を参照すると、本発明による計算例では、同じ繰り返し数にもかかわらずほとんどの要素で材料は充填状態と空隙状態に分離している。
【００７２】
両手法共に材料マトリックスを主軸に合わせて回転しているので，結果の差は剪断剛性の強さに起因すると思われる。
【００７３】
図１２に従来例および本実施例による目的関数の履歴を示す。両手法を比較すると、ユニットセル形状の違いによる収束性の差はほとんど認められなかった。
【００７４】
図８に、直方体の設計領域（１／２対称モデル）の下から約１／３の水平帯状部分に分布荷重を載荷した、橋梁をイメージしたモデルについて位相を計算した例を示す。支持条件は橋脚ができるであろう下部の１箇所と橋桁の端部を完全拘束し，対称面に対称条件を付加した。
【００７５】
制約条件は設計領域の１０％の体積とし、分布荷重を載荷した部分を含む要素２層分は始めから材料を充填状態にして、計算中も一切変更をしない非設計領域とした。
【００７６】
計算モデルはミニチュアモデル、すなわち縮小した計算モデルである。このモデルの寸法は０．３３×０．０６×０．２３［ｍ］、ヤング率は２０５８００［Ｍｐａ］、ポアソン比は０．３、総節点数は５７１２、総要素数は５５４、また、荷重は単位面積あたり９．８［Ｍｐａ］である。
【００７７】
図９、１０および１１に、このモデルの計算例を示す。この計算例は、前述した短い片持ち梁の計算例に比べて、さらに大きな効果が現れている。従来の方法は、設計領域から９０％の材料を削り取ったにもかかわらず、鉄橋を構成する主要な部材を判断するに足る情報を算出したが、かなり曖昧な部分が残った。本発明の方法はさらに鮮明な位相を算出した。剪断剛性の差がこのモデルではさらに大きな影響を与えたと推測できる。両手法とも我々が日常目撃する鉄道橋全体のイメージは再現している。
【００７８】
【発明の効果】
実施例でも解るように、本発明によれば、図２に示す従来のユニットセルを使用した場合に比べ、位相最適化の結果に現れる中間密度の部分を大幅に減らし、認識しやすい形状を得ることができる。
【図面の簡単な説明】
【図１】本発明によるフレームベースユニットセルを示した図である。
【図２】従来のユニットセルを示した図である。
【図３】位相最適化方法の概要を示したフローチャートである。
【図４】短い片持ち梁の有限要素モデルを示した図である。
【図５】従来の均質化法による位相を示した図である。
【図６】本実施例による位相を示した図である。
【図７】本実施例による位相を示した平面図である。
【図８】橋をイメージした有限要素モデルを示した図である。
【図９】従来の均質化法による位相を示した図である。
【図１０】本実施例による位相を示した図である。
【図１１】本実施例による位相を示した図である。
【図１２】目的関数の履歴を示したグラフである。
【符号の説明】
ａ、ｂ、ｃ    内辺の長さ
【図１】
0003861259.tif 000018
【図２】
0003861259.tif 000019
【図３】
0003861259.tif 000020
【図４】
0003861259.tif 000021
【図５】
0003861259.tif 000022
【図６】
0003861259.tif 000023
【図７】
0003861259.tif 000024
【図８】
0003861259.tif 000025
【図９】
0003861259.tif 000026
【図１０】
0003861259.tif 000027
【図１１】
0003861259.tif 000028
【図１２】
0003861259.tif 000029