(57)【要約】【課題】構造物の周波数応答性を迅速に計算する。【解決手段】主構造であるシェル要素１０を非設計要素とし、補強材としての梁要素１２を設計要素とする。シェル要素１０と梁要素１２の動剛性を用い、目的関数を単一または複数の周波数の振動騒音特性として、体積の制約条件下でトポロジー最適化を行う。主構造として、特定の部分構造に縮退した動剛性を用いることもできる。【選択図】図１
2010067022.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを用いた構造物の設計方法であって、
コンピュータに、
主構造に対して隣接する節点間の全てに補強材として梁要素を配置し、前記主構造の特性として剛性、質量、減衰を含む動剛性を与え、前記主構造を非設計要素、前記梁要素を設計要素に設定するステップと、
単一または複数の周波数応答特性を目的関数とし、体積の制約条件下において前記設計要素を更新して前記目的関数が最小となる配置を探索するステップと、
を実行させることを特徴とする構造物の設計方法。
【請求項２】
請求項１記載の方法において、
前記主構造は、特定の部分構造に縮退した動剛性であることを特徴とする構造物の設計方法。
【請求項３】
コンピュータに構造物のモデルを設計させるプログラムであって、
前記プログラムはコンピュータに、
主構造に対して隣接する節点間の全てに補強材として梁要素を配置し、前記主構造の特性として剛性、質量、減衰を含む動剛性を与え、前記主構造を非設計要素、前記梁要素を設計要素に設定するステップと、
単一または複数の周波数応答特性を目的関数とし、体積の制約条件下において前記設計要素を更新して前記目的関数が最小となる配置を探索するステップと、
を実行させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は構造物の設計方法に関し、特に梁要素を用いたトポロジー最適設計方法に関する。
【背景技術】
【０００２】
ＣＡＥは自動車等の製品開発に広く適用されている。そこでは１００万自由度を超える計算モデルが高い精度を保証するために活用され、試作車両の軽減を可能としコストや期間の短縮に貢献している。しかし、計算結果が目標を満足していない場合、その要因を特定し効果的な対策案を創出することに多くの時間と労力を必要とする。特に、構造物の振動問題の場合、最小の質量増加で対策部位とその対策構造を見出すことが重要な課題となっている。一般に、最適化問題は、寸法最適化問題、形状最適化問題、トポロジー最適化問題に分類することができる。振動問題では、形状最適化や位相最適化の適用は部品レベルに留まっているのが現状である。大規模構造の振動問題では、寸法最適化問題が利用されている。
【０００３】
トポロジー最適化方法を振動問題に適用した研究には、平板等の基本構造を対象とした固有値問題や周波数応答問題などがある。
【０００４】
一方、離散化構造要素のはり要素を用いたトポロジー最適設計問題も剛性問題や固有値問題で扱われている。下記の特許文献１には、隣接節点及び非隣接節点同士を６軸の剛性をもつ梁要素で結合させた設計モデルを作成し、この設計モデルについて最適化を行う目的関数である平均コンプライアンスとして印加する力と変位の内積を設定し、所定の体積制約条件または重量制約条件を満足する中で、設計モデルにおける梁要素の剛性を変化させ、平均コンプライアンスについての感度を検出し、検出した感度に応じて平均コンプライアンスを最小化することで、設計モデルの剛性が最大となる形状に最適化する技術が開示されている。
【０００５】
また、大規模問題を簡便に解くために、縮退技術により自由度を小さくする手法も知られている。モーダル座標系による縮退法と、物理座標系を維持したままの縮退法に大別される。下記の特許文献２には、後者の方法が開示されている。
【０００６】
【特許文献１】特許第３５５１９１０号公報
【特許文献２】特開２００１−１２６０８７号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
寸法最適化問題では、計算上の制約から設計変数を限定するため、設計空間上の全ての可能性を考慮した上での最適値とならない問題がある。
【０００８】
また、連続体要素を用いたトポロジー最適設計法では、車体構造等の大規模モデルに適用した場合に、多くの処理時間を必要とするため、実用的なツールとならない問題がある。
【０００９】
一方、梁要素を用いたトポロジー最適設計法は剛性問題や固有値問題に適用されているが、周波数応答問題には全く適用されていない。また、基本構造が主体であり、概念設計段階では有効な方法ではあるが、詳細設計段階には適さない問題がある。
【００１０】
本発明の目的は、梁要素を用いたトポロジー最適設計方法において、任意の周波数における周波数応答特性を解析することで最小の質量増加で対策部位とその対策構造を見出すことができる方法を提供することにある。
【課題を解決するための手段】
【００１１】
コンピュータを用いた構造物の設計方法であって、コンピュータに、主構造に対して隣接する節点間の全てに補強材として梁要素を配置し、前記主構造の特性として剛性、質量、減衰を含む動剛性を与え、前記主構造を非設計要素、前記梁要素を設計要素に設定するステップと、単一または複数の周波数応答特性を目的関数とし、体積の制約条件下において前記設計要素を更新して前記目的関数が最小となる配置を探索するステップとを実行させることを特徴とする。
【００１２】
本発明において、前記主構造は、特定の部分構造に縮退した動剛性としてもよい。
【発明の効果】
【００１３】
本発明によれば、対象構造の動剛性を用いることで、周波数応答特性を得ることができる。また、本発明によれば、規定の質量増加で最大の効果が得られる対策部位とその構造特性が得られる。また、本発明によれば、梁要素を設計要素とすることで、安定かつ高速な最適化演算を行うことができる。また、本発明によれば主構造を非設計要素とすることで、全体剛性を保持したまま、任意の周波数特性だけの改善が可能である。
【発明を実施するための最良の形態】
【００１４】
以下、図面に基づき本発明の実施形態について説明する。
【００１５】
図１に、本実施形態における構造要素で構成された対象領域１を示す。対象領域は、シェル要素１０と梁要素１２から構成される。梁要素１２は、隣接節点間及び非隣接節点間同士を結び、６軸の剛性を有する。ある節点（ノード）１４に加振し、任意の周波数に対する評価節点１６の振幅ｗを目的関数とする。必要に応じ、対象領域１の中に固定領域１８を設定する。この構造物の離散化された運動方程式は次式となる。
【００１６】
【数１】
2010067022.tif 000003
【００１７】
ここで、｛Ｆ｝、｛Ｕ｝はそれぞれ加振ベクトル、変位ベクトルであり、[Ｍ]、[Ｃ]、[Ｋ]はそれぞれ質量、減衰、剛性の行列、[Ｂ]は動剛性行列、ωは角振動数を表す。相互平均コンプライアンスlを用いると、評価節点１６の応答は次式で示される。
【００１８】
【数２】
2010067022.tif 000004
また、目的関数ｗ（ω）は次式となる。
【００１９】
【数３】
2010067022.tif 000005
ここで、｛Ｇ｝、｛Ｖ｝はそれぞれ評価節点１６の単位加振ベクトルとその力が作用したときの変位ベクトルを示す。また、l*はlの共役値を示す。
【００２０】
次に、図１に示すようなシェル要素１０がベース構造となる固定設計領域に梁要素１２をグランドストラクチャ法に基づき配置する。グランドストラクチャ法では、隣接する節点間に全て梁要素１２を配置する。なお、必要に応じ、非隣接の節点間までも拡張して梁要素１２を配置することも可能である。トポロジー最適化においては、固定設計領域の設定とその固定設計領域に配置する構造要素の有無を判定するには、次式の特性関数に基づく。
【００２１】
【数４】
2010067022.tif 000006
【００２２】
ここで、Ωは固定設計領域、Ωｄは固定設計領域から最適化により選択させた最適構造を形成する領域、ｘは固定設計領域Ω内の位置座標である。この特性関数を用いることで、最適設計問題は要素配置問題に置き換えることができる。この特性関数を用いると離散化された変数を取り扱わねばならず、計算上の取り扱いは不可能となる。そこで、密度法的な考え方に基づいた次式により特性関数を連続関数に置き換える。
【００２３】
【数５】
2010067022.tif 000007
【００２４】
ここで、ρｖは正規化された設計変数である。また、ρはρｖにペナルティを与えるパラメータである。上式を用いることで、最適設計問題を固定設計領域Ω内における要素の連続的な分布問題に置き換えることができる。本実施形態では、梁要素１２のみを設計要素として取り扱い、シェル要素１０は非設計要素とする。上式を最大時の断面積Ａｍａｘ、長さＬの梁要素１２に適用すれば、その体積Ｖは次式となり、ρｖの値により梁要素１２の配置が決定できる。この正規化された面積密度ρｖを各要素の設計変数として、最適化問題を設定する。
【００２５】
【数６】
2010067022.tif 000008
【００２６】
図２に、薄肉断面の梁要素１２を示す。中実断面と比較して、薄肉断面は、曲げ剛性及びねじり剛性をより高くすることができる。そのため曲げやねじり剛性の寄与が大きい場合、少ない付加質量で必要な低減効果が得られる可能性がある。直径と肉厚との関係は、薄肉比γにより次式で定義する。
【００２７】
【数７】
2010067022.tif 000009
【００２８】
ここで、γ＝１は中実断面を示す。上式を用いて、断面積Ａ、断面二次モーメントＩｙ、Ｉｚ、断面二次極モーメントＪ、質量ｍは次式で定義される。
【００２９】
【数８】
2010067022.tif 000010
【００３０】
【数９】
2010067022.tif 000011
【００３１】
【数１０】
2010067022.tif 000012
【００３２】
【数１１】
2010067022.tif 000013
【００３３】
ここで、ρｍは要素質量密度を示す。これらの梁の断面特性から動剛性行列を作成する。さらに、主構造となるシェル要素１０の動剛性行列を重ね合わせた全体行列を作成し、周波数応答解析から目的関数の値を導出する。
【００３４】
数式（１）の動剛性行列をベースのシェル要素１０と梁要素１２の動剛性行列[Ｂ０]、[Ｂａ]で示すと次式となる。
【００３５】
【数１２】
2010067022.tif 000014
正規化された設計変数ρｖに対する相互平均コンプライアンスの感度は、次式で示される。
【００３６】
【数１３】
2010067022.tif 000015
この感度は、梁要素１２の動剛性Ｂａに対する設計変数の感度と２種類の応答ベクトル｛Ｖ｝、｛Ｕ｝から導出することができる。さらに、振幅に対する感度は、次式で示される。
【００３７】
【数１４】
2010067022.tif 000016
【００３８】
この感度を最適化アルゴリズムに適用し、設計変数を更新することで最適解を導出する。なお、上記のベースとなる動剛性行列として縮退した行列を用いることも可能である。動剛性行列を縮退する方法として伝達関数による縮退を用いる。その理由として、剛性付加の効果を直接評価できること、及び縮退後も振動応答の出力結果の精度が保証されているからである。縮退する自由度間の伝達関数が保持されていれば、縮退前後の入力と出力の関係は自由度の数に依存せず数値誤差以外の差は存在しない。剛性値の感度は振動応答の出力結果から求められるため、最適化計算を行う上で縮退後の出力結果の精度が重要である。
【００３９】
動剛性行列を縮退する自由度の領域Ｏと消去する自由度の領域ｂに分離し、次式のように表す。
【００４０】
【数１５】
2010067022.tif 000017
上式は領域に外力が作用しないことから、縮退した動剛性行列は次式で表される。
【００４１】
【数１６】
2010067022.tif 000018
図１に示す固定領域を対象に、角周波数ωｊの評価節点１６の振幅値をωｊとする。目的関数Ｗは、考慮するｍラインの周波数の重み付き線形和で表す。ｃjは各周波数の重み係数を表す。シェル要素１０は非設計要素とし、ｎ本の梁要素１２を設計要素とする。以下に、最適化問題の定式化を示す。評価関数Ｗを最小化するものである。
【００４２】
【数１７】
2010067022.tif 000019
但し、
【００４３】
【数１８】
2010067022.tif 000020
【００４４】
【数１９】
2010067022.tif 000021
【００４５】
【数２０】
2010067022.tif 000022
【００４６】
【数２１】
2010067022.tif 000023
である。ここで、ＶΩ、ＶUはそれぞれ梁要素１２の総体積と体積制約の上限値である。｛Ｆ｝jと｛Ｇ｝jが同一である場合には、目的関数は平均コンプライアンスに退化し、加振点応答の振幅値となる。
【００４７】
図３に、本実施形態における最適化の手順をフローチャートとして示す。まず、固定領域Ωに梁要素１２の初期値と非設計要素のシェル要素１０の配置を図１のように与え、任意の周波数ωｊに対する動剛性行列[Ｂ]を作成する（Ｓ１０１）。あるいは、縮退した動剛性行列を作成してもよい。
【００４８】
次に、設計変数ρｖ，ｉに対する動剛性行列の感度を算出する（Ｓ１０２）。ここでは、ペナルティを与えるパラメータｐは初期値の１とする。
【００４９】
次に、任意の周波数において、釣り合い方程式を解き、応答ベクトル｛Ｕ｝j、｛Ｖ｝jを求める（Ｓ１０３）。
【００５０】
次に、相互平均コンプライアンスｌj、目的関数Ｗ、体積制約ＶΩの値を求める（Ｓ１０４）。
【００５１】
目的関数を算出した後、算出された目的関数が収束したか否か、つまり最小値に達したか否かを判定する（Ｓ１０５）。収束したか否かは、前回に算出した値との差分を算出し、差分値が所定の許容値より小さいか否かで判定できる。
【００５２】
目的関数が収束した場合にはその時点で処理を終了する。一方、目的関数が収束していない場合には、相互平均コンプライアンス、目的関数、体積制約の設計変数ρｖ，ｉに対する感度を算出する（Ｓ１０６）。
【００５３】
感度を算出した後、得られた感度をもとに、最適化手法を用いて設計変数ρｖ，ｉを更新して（Ｓ１０７）、再びＳ１０１以降の処理を繰り返し実行する。Ｓ１０１以降の処理は、Ｓ１０５で収束と判定されるまで繰り返し実行される。
【００５４】
トポロジー最適化問題では、最適化手法として最適性基準法、逐次線形計、逐次凸関数近似法等が用いられる。特に、逐次凸関数近似法の１つであるＣＯＮＬＩＮ（Convex Linearization）は剛性問題に適用され、解の安定性と速い収束性が実証されている。ＣＯＮＬＩＮの基本的な考え方は、目的関数または制約条件である関数ｈを次式のパラメータηｉを用いて近似することにある。本実施形態でも、最適化手法にはこのＣＯＬＩＮを用いるのが好適である。
【００５５】
【数２２】
2010067022.tif 000024
【００５６】
剛性問題では目的関数の設計変数ρｖ，ｉに対する感度は負となる。このため上式の２番目の式が適用される。他方、設計変数ρｖ，ｉに関して目的関数は反比例に近い関数となる。つまり、パラメータに対して線形に近い近似となるため高い収束性を実現できる。
【００５７】
以下、本実施形態の処理を具体的に示す。
【００５８】
図４に示すように、平板に補強材を体積制約下で配置する場合を対象とする。４つの頂点を単純支持し（図中三角点で示す）、中心に面外方向の単位周期外力を与える（図中矢印で示す）。材料は鉄を想定し、ヤング率２．１×１０-5Ｎ／ｍｍ2、ポアソン比０．３、質量密度７．８６×１０-6ｋｇ／ｍｍ3とする。ベースのシェル要素の板厚を５ｍｍ、その上に配置した梁要素１２を設計要素とする。設計変数である正規化した梁の断面積の初期値をρｖｉ＝１０-4とする。また、薄肉比γ＝０．２、最大直径ｄｍａｘ＝１０ｍｍとする。シェルの要素数は１６、梁の要素数は７２である。５００Ｈｚ以下の任意の周波数において加振点の振幅を目的関数とする。付加する梁要素１２の体積を制約条件とし、最大体積の２５％以下とする。なお、梁要素１２の剛性、質量、減衰を含む動剛性を用いるが、シェル要素で構成される全体の質量４．７ｋｇに対して梁要素単体の最大質量は２８ｇであり１％に満たない。また、要素単体の固有値も１０ＫＨｚ以上と高い、つまり対象とする５００Ｈｚ以下の周波数特性への影響は小さいと考えられる。そこで、計算上、梁要素１２の剛性のみを考慮し、梁要素１２の質量を実質的に無視することができる。
【００５９】
図５に、初期状態の加振点コンプライアンスを示す。６８Ｈｚ、３７１Ｈｚにピークが生じる。静的変形モード（０Ｈｚ）とピークを形成するモード１（６８Ｈｚ）、モード２（３７１Ｈｚ）の変形モードを歪みエネルギ分布とともに図６（ａ）、（ｂ）、（ｃ）に示す。
【００６０】
図６（ａ）は静的モード、図６（ｂ）はモード１、図６（ｃ）はモード２の歪みエネルギ分布である。静的モードとモード１の変形は同等であるが、歪みエネルギ分布が異なることがわかる。また、２次のピークはモードと歪みエネルギの双方でいずれも異なる。図６に示す３つのモードの周波数について最適化問題の解を探索した。
【００６１】
０Ｈｚ、６８Ｈｚ、３７１Ｈｚを対象とした最適化結果をＯｐｔ１、Ｏｐｔ２、Ｏｐｔ３とし、図７に各周波数の最適化結果を示す。また、図８にＯｐｔ１の遷移図、図９にＯｐｔ２の遷移図、図１０にＯｐｔ３の遷移図をそれぞれ示す。これらの図において、ρｖｉが０．９５以上となる梁要素１２のみを示す。
【００６２】
図１１に、最適化結果に基づくピーク周波数の変化を示す。また、初期値に対する正規化したピーク周波数を図１２に示す。ピークレベルの低減の良否は、ピーク周波数の変化で表すことができる。また、図１３に、最適結果に基づくピーク値を示す。また、図１４に、初期値に対する正規化したピーク値を示す。最適配置の結果は、歪みエネルギ分布の高い部位に対応し、各モードにおいて妥当な補強が与えられている。Ｏｐｔ１とＯｐｔ２において、最適化の結果は一部異なるが、ピークの周波数は同等である。これは、変形モードや歪みエネルギ分布に大きな相違がなく、同等な効果となったものと考えられる。Ｏｐｔ３はモード２のピーク周波数を最も高くしているが、モード１は最も変化していない。これらの結果より各周波数において最適配置は異なり、高い低減効果を得るためには対象とする周波数で最適化を実行する必要がある。また、それぞれの最適化結果は、歪みエネルギ分布に対応しており、力学的にも妥当な結果であると考えられる。
【００６３】
このように、本実施形態では、対象構造の動剛性特性を用いることで、静的応答計算と同等に周波数応答計算が可能となる。また、グランドストラクチャ法に基づき、設計空間の全てのトポロジー（レイアウト）を設計変更の候補とすることができ、これにより従来の限定された設計変数による最適化から全ての可能性を考慮した最適化が可能となる。
【００６４】
また、本実施形態では、主構造であるシェル要素を非設計変数とすることで、主体となるトポロジー（レイアウト）が維持される。これにより、安定した最適化計算が可能である。主構造までも設計変数とした場合、極端なトポロジーの変更が解の安定性を乱し、適切な解が得られない。また、同様の設計空間でもシェル要素やソリッド要素を使用することで自由度が２乗、３乗で増加する。これらと比較すると、１次元要素の梁要素は自由度を小さくすることができる。また、逐次凸関数近似法を用いることで、高速な収束を実現することができる。
【００６５】
また、１００万自由度を超えるような車体構造に適用する場合でも、縮退技術により精度を低下させることなく自由度を縮小することができる。
【００６６】
さらに、本実施形態では既知である主構造は非設計要素とするため、初期の剛性は確保されている。これに対して、任意の周波数における補強材の最適配置問題として解くため、剛性確保と振動低減の両立が可能である。そもそも車両構造では全体構造が衝突や剛性で決まるため、振動問題は最適な補強材配置問題として扱うことで実用的な解が得られる。
【００６７】
本実施形態の設計方法は、具体的にはパーソナルコンピュータやワークステーションなどのコンピュータに処理プログラムをインストールし、コンピュータ上で実行することができる。本実施形態の設計モデルは、シェル要素と梁要素のような構造力学に基づいた要素を使用しているため、大容量のメモリは不要であり、例えばパーソナルコンピュータ上で表計算ソフトウェアフロントエンドとして用いることも可能である。なお、処理プログラムはＣＤ−ＲＯＭやＤＶＤ、ハードディスク、半導体メモリ等、電磁気的、光学的、あるいは化学的に情報を保持できる任意の媒体に記録することができる。処理プログラムのインストールは、例えばＣＤ−ＲＯＭに処理プログラムを記録し、ＣＤ−ＲＯＭからコンピュータのハードディスクに処理プログラムを供給することで行うことができる。もちろん、当初からコンピュータのハードディスクあるいはＲＯＭに処理プログラムを記憶させ、構造物モデル設計専用のコンピュータとして用いてもよい。
【図面の簡単な説明】
【００６８】
【図１】実施形態の構造要素を示す説明図である。
【図２】梁要素の説明図である。
【図３】実施形態の処理フローチャートである。
【図４】平板に補強材を配置する場合の説明図である。
【図５】初期状態の加振点コンプライアンス説明図である。
【図６】静的モード、モード１、モード２の歪みエネルギ分布説明図である。
【図７】静的モード、モード１、モード２の最適結果を示す配置説明図である。
【図８】静的モードの最適化に伴う遷移図である。
【図９】モード１の最適化に伴う遷移図である。
【図１０】モード２の最適化に伴う遷移図である。
【図１１】モード１、モード２の最適化によるピーク周波数変化を示す表図である。
【図１２】モード１、モード２の最適化による正規化ピーク周波数を示すグラフ図である。
【図１３】モード１、モード２の最適化によるピークレベル変化を示す表図である。
【図１４】モード１、モード２の最適化による正規化ピークレベルを示すグラフ図である。
【符号の説明】
【００６９】
１０  シェル要素、１２  梁要素。
【図１】
2010067022.tif 000025
【図２】
2010067022.tif 000026
【図３】
2010067022.tif 000027
【図４】
2010067022.tif 000028
【図５】
2010067022.tif 000029
【図７】
2010067022.tif 000030
【図８】
2010067022.tif 000031
【図９】
2010067022.tif 000032
【図１０】
2010067022.tif 000033
【図１１】
2010067022.tif 000034
【図１２】
2010067022.tif 000035
【図１３】
2010067022.tif 000036
【図１４】
2010067022.tif 000037
【図６】
2010067022.tif 000038