明細書 :構造最適化装置、構造最適化方法及び構造最適化プログラム

発行国	日本国特許庁(JP)
公報種別	特許公報(B2)
特許番号	特許第5377501号 (P5377501)
登録日	平成25年10月4日(2013.10.4)
発行日	平成25年12月25日(2013.12.25)
発明の名称または考案の名称	構造最適化装置、構造最適化方法及び構造最適化プログラム
国際特許分類	G06F  17/50        (2006.01)
ＦＩ	G06F 17/50 604A G06F 17/50 612H
請求項の数または発明の数	5
全頁数	22
出願番号	特願2010-528686 (P2010-528686)
出願日	平成21年7月17日(2009.7.17)
国際出願番号	PCT/JP2009/062974
国際公開番号	WO2010/029810
国際公開日	平成22年3月18日(2010.3.18)
優先権出願番号	2008233176
優先日	平成20年9月11日(2008.9.11)
優先権主張国	日本国(JP)
審査請求日	平成24年7月10日(2012.7.10)
特許権者または実用新案権者	【識別番号】504132272 【氏名又は名称】国立大学法人京都大学
発明者または考案者	【氏名】山田 崇恭 【氏名】西脇 眞二 【氏名】泉井 一浩 【氏名】吉村 允孝
個別代理人の代理人	【識別番号】100121441、【弁理士】、【氏名又は名称】西村 竜平 【識別番号】100113468、【弁理士】、【氏名又は名称】佐藤 明子 【識別番号】100154704、【弁理士】、【氏名又は名称】齊藤 真大
審査官	【審査官】加舎 理紅子
参考文献・文献	特開２００３－１５０６５１（ＪＰ，Ａ） 特開２００５－２５８８１３（ＪＰ，Ａ） Lin He et al.，Incorporating topological derivatives into shape derivatives based level set methods，Journal of Computational Physics，２００７年 ７月 １日，Volume 225, Issue 1，p.891-909 Martin Burger et al.，Incorporating topological derivatives into level set methods，Journal of Computational Physics，２００４年 ２月１０日，Volume 194, Issue 1，p.344-362 山田崇恭 他，レベルセット法に基づくコンプライアントメカニズムの構造最適化，日本計算工学会論文集，２００８年 １月１８日，Vol.2008,No.20080001
調査した分野	Ｇ０６Ｆ １７／５０ ＣｉＮｉｉ ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）
特許請求の範囲	【請求項1】 構造物の設計領域を示す設計領域データを格納する設計領域データ格納部と、 初期構造が設定された設計領域の各部が構造を形成する物体領域、空洞を形成する空洞領域、又はそれら領域の境界であるかを示し、物体領域を表す値及び空洞領域を表す値の間の所定値が前記境界を表すレベルセット関数を示すレベルセット関数データを格納しているレベルセット関数データ格納部と、 所定の制約条件下において、剛性等の構造物の性能を目標値に近づけるように前記レベルセット関数を更新して、前記物体領域内のトポロジー変化（形態変化）及び前記物体領域と前記空洞領域との境界の移動を同時に行わせるレベルセット関数更新部と、を具備する構造最適化装置。 【請求項2】 構造物の設計領域を示す設計領域データを格納する設計領域データ格納部と、 初期構造が設定された設計領域の各部が構造を形成する物体領域、空洞を形成する空洞領域、又はそれら領域の境界であるかを示し、物体領域を表す値及び空洞領域を表す値の間の所定値が前記境界を表すレベルセット関数を示すレベルセット関数データを格納しているレベルセット関数データ格納部と、 所定の制約条件下において、剛性等の構造物の性能を目標値に近づけるように前記レベルセット関数を更新して、前記物体領域と前記空洞領域との境界を移動させるとともに、前記レベルセット関数の更新に伴う前記物体領域内のトポロジー変化（形態変化）を許容して前記物体領域内に新たな空洞領域を形成し、その新たな空洞領域と物体領域との境界を移動させるレベルセット関数更新部とを具備し、 前記レベルセット関数更新部が、レベルセット関数を変数とする関数族、物体領域におけるエネルギー密度、空洞領域におけるエネルギー密度、及び界面エネルギー密度により示されるエネルギー汎関数から、エネルギー汎関数最小化原理に従って、前記レベルセット関数の時間発展を示す反応拡散方程式を算出し、当該反応拡散方程式を用いて前記レベルセット関数を時間発展させることにより、前記レベルセット関数を更新するものである構造最適化装置。 【請求項3】 前記レベルセット関数更新部が、前記レベルセット関数を更新した結果得られる構造の構造的な複雑さを示す複雑度が予め設定された複雑度となるように、前記レベルセット関数を更新するものである請求項１記載の構造最適化装置。 【請求項4】 入力手段により、構造物の設計領域を示す設計領域データをコンピュータに入力する設計領域データ入力ステップと、 入力手段により、初期構造が設定された設計領域の各部が構造を形成する物体領域、空洞を形成する空洞領域、又はそれら領域の境界であるかを示し、物体領域を表す値及び空洞領域を表す値の間の所定値が前記境界を表すレベルセット関数を示すレベルセット関数データを前記コンピュータに入力するレベルセット関数入力ステップと、 前記コンピュータに設けられたレベルセット関数更新部により、所定の制約条件下において、剛性等の構造物の性能を目標値に近づけるように前記レベルセット関数を更新して、前記物体領域内のトポロジー変化（形態変化）及び前記物体領域と前記空洞領域との境界の移動を同時に行わせるレベルセット関数更新ステップと、を具備する構造最適化方法。 【請求項5】 構造物の設計領域を示す設計領域データを格納する設計領域データ格納部と、 初期構造が設定された設計領域の各部が構造を形成する物体領域、空洞を形成する空洞領域、又はそれら領域の境界であるかを示し、物体領域を表す値及び空洞領域を表す値の間の所定値が前記境界を表すレベルセット関数を示すレベルセット関数データを格納しているレベルセット関数データ格納部と、 所定の制約条件下において、剛性等の構造物の性能を目標値に近づけるように前記レベルセット関数を更新して、前記物体領域内のトポロジー変化（形態変化）及び前記物体領域と前記空洞領域との境界の移動を同時に行わせるレベルセット関数更新部と、としての機能をコンピュータに備えさせる構造最適化プログラム。
発明の詳細な説明	【技術分野】 【0001】 本発明は、構造最適化装置、構造最適化方法及び構造最適化プログラムに関するものである。 【背景技術】 【0002】 従来、構造最適化の方法としては、寸法最適化、形状最適化、トポロジー最適化がある。 【0003】 この中の形状最適化は、外形形状を設計変数として、感度情報に基づいて前記外形形状を更新することにより最適構造を得る方法であり、実用的な方法として機械産業において広く利用されている。 【0004】 しかしながら、形状最適化は、最適構造の穴の数などの形態を変更することが難しいという欠点をもち、構造の大幅な性能向上を期待できない。 【0005】 これに対して、トポロジー最適化は、最適設計問題を材料分布問題に置き換えて解くことにより、最適構造の形態変更を可能とし、構造の大幅な性能向上を期待できる。 【0006】 しかしながら、トポロジー最適化はグレースケールなどの数値不安定性問題を生じることがある。 【0007】 そして近年、新しい構造最適化方法として、レベルセット法に基づく構造最適化が提案されている（非特許文献１）。 【0008】 この方法は、外形形状を１次元高位のレベルセット関数で表現し、形状と形態の変更をレベルセット関数値の変化に置き換えて最適構造を得る。これにより従来のトポロジー最適化とは異なり、最適構造の輪郭を常に明確に表現することが可能であり、グレースケールなどの問題を生じることがないという長所を有する。 【0009】 しかしながら、レベルセット法は、移流方程式に基づいてレベルセット関数を更新させるため、外形形状（物体領域）に穴が創出されるようなトポロジー変化（形態変化）を許容しないことを前提としている。 【0010】 これに対して、非特許文献２に示すように、レベルセット法に基づく構造最適化の過程において、適宜、目的関数のトポロジカルデリバティブに基づいて、恣意的に穴を創出する方法が提案されている。また、非特許文献３に示すように、最適化の過程において、トポロジカルデリバティブの値が所定の閾値である物体領域（構造を形成する領域）を空洞領域（空洞を形成する領域）に置き換えることにより、外形形状（物体領域）に穴が創出されるようなトポロジー変化（形態変化）を許容する方法が提案されている。 【0011】 しかしながら、いずれの方法も、穴の数又は閾知の設定などのパラメータの依存性が極めて高く、適切にそれらのパラメータを設定しなければ、物理的に妥当な最適構造が得られないという問題がある（非特許文献４参照） 【先行技術文献】 【0012】 【非特許文献1】Wang, M. Y., Wang, X. and Guo, D., A level Set Method forStructural Topology Optimization, Computer Methods in Applied Mechanics andEngineering, Vol.192, (2003), pp.227-246. 【非特許文献2】Yamasaki, S., Nishiwaki, S., Yamada, T., Izui, K., andYoshimura, M., A Structural Optimization Method Based on the Level Set MethodUsing A New Geometry-based Re-initialization Scheme, International Journal for NumericalMethods in Engineering,Vol. 18, (2008), pp.487-505. 【非特許文献3】Park, K. S., and Youn, S. K., Topology Optimization ofshell structures using adaptive inner-front (AIF) level set method, Structuraland Multidisciplinary Optimization, Vol.36, (2008), pp.43-58. 【非特許文献4】Yamada, T., Nishiwaki, S., Izui, K., and Yoshimura, M., AStudy of Boundary Setting in the Design Domain for Structural OptimizationBased on the Level Set Method, Transaction of the Japan Society for Industrialand Applied Mathematics, (submitted). 【発明の概要】 【発明が解決しようとする課題】 【0013】 そこで本発明は、上記問題点を一挙に解決するためになされたものであり、最適構造の明確な形状表現を可能にし、物体領域でのトポロジー変化を許容する等の自由度の高い構造最適化を可能にすることをその主たる所期課題とするものである。 【課題を解決するための手段】 【0014】 すなわち本発明に構造最適化装置は、以下の構成（１）～（３）を具備することを特徴とする。 【0015】 （１）構造物の設計領域を示す設計領域データを格納する設計領域データ格納部。 （２）初期構造が設定された設計領域の各部が構造を形成する物体領域、空洞を形成する空洞領域、又はそれら領域の境界であるかを示し、物体領域を表す値及び空洞領域を表す値の間の所定値が前記境界を表すレベルセット関数を示すレベルセット関数データを格納しているレベルセット関数データ格納部。 （３）所定の制約条件下において、剛性等の構造物の性能を目標値に近づけるように前記レベルセット関数を更新して、前記物体領域と前記空洞領域との境界の形態及び形状を同時に変化させるレベルセット関数更新部。 【0016】 このようなものであれば、物体領域を表す値及び空洞領域を表す値の間の所定値が、物体領域及び空洞領域の境界を示すレベルセット関数を用いていることにより、最適構造の明確な形状表現を可能にすることができる。また、物体領域及び空洞領域間の境界を移動させるとともに、物体領域内でのトポロジー変化を許容して新たに生じる境界を移動させるようにレベルセット関数を更新することにより、物体領域でのトポロジー変化を許容する等の自由度の高い構造最適化を可能にすることができる。 【0017】 レベルセット関数の更新方法の具体的な実施の態様としては、前記レベルセット関数更新部が、レベルセット関数を変数とする関数族、物体領域におけるエネルギー密度、空洞領域におけるエネルギー密度、及び界面エネルギー密度により示されるエネルギー汎関数を、エネルギー汎関数最小化原理に従って、前記レベルセット関数の時間発展を示す反応拡散方程式を算出し、当該反応拡散方程式を用いて前記レベルセット関数を時間発展させることにより、前記レベルセット関数を更新するものであることが望ましい。 【0018】 また、前記レベルセット関数更新部が、前記レベルセット関数を更新した結果得られる構造の構造的な複雑さを示す複雑度が予め設定された複雑度となるように、前記レベルセット関数を更新するものであることが望ましい。これならば、無数に存在する局所最適解の中で、設計者（ユーザ）の意図した構造の複雑さ（つまり細かさ）を有する最適構造を創成することができる。 【0019】 また本発明に係る構造最適化方法は、構造物の設計領域を定める設計領域設定ステップと、初期構造が設定された設計領域の各部が構造を形成する物体領域、空洞を形成する空洞領域、又はそれら領域の境界であるかを示し、物体領域を表す値及び空洞領域を表す値の間の所定値が前記境界を表すレベルセット関数を定めるレベルセット関数設定ステップと、所定の制約条件下において、剛性等の構造物の性能を目標値に近づけるように前記レベルセット関数を更新して、前記物体領域と前記空洞領域との境界の形態及び形状を同時に変化させるレベルセット関数更新ステップと、を具備する。 【0020】 さらに本発明に係る構造最適化プログラムは、構造物の設計領域を示す設計領域データを格納する設計領域データ格納部と、初期構造が設定された設計領域の各部が構造を形成する物体領域、空洞を形成する空洞領域、又はそれら領域の境界であるかを示し、物体領域を表す値及び空洞領域を表す値の間の所定値が前記境界を表すレベルセット関数を示すレベルセット関数データを格納しているレベルセット関数データ格納部と、所定の制約条件下において、剛性等の構造物の性能を目標値に近づけるように前記レベルセット関数を更新して、前記物体領域と前記空洞領域との境界の形態及び形状を同時に変化させるレベルセット関数更新部と、としての機能をコンピュータに備えさせることを特徴とする。 【発明の効果】 【0021】 このように本発明によれば、最適構造の明確な形状表現を可能にし、物体領域でのトポロジー変化を許容する等の自由度の高い構造最適化を可能にすることができる。 【図面の簡単な説明】 【0022】 【図1】本実施形態に係る構造最適化装置の機器構成図である。 【図2】同実施形態に係る構造最適化装置の機能構成図である。 【図3】同実施形態の最適化アルゴリズムを示すフローチャートである。 【図4】レベルセット法における物体領域、空洞領域及び境界を示す図である。 【図5】設計問題１及び設計問題２の設計領域と境界条件を示す図である。 【図6】Ｋ１（φ）＝１の初期構造、最適化過程の構造及び最適構造を示す図である。 【図7】Ｋ２（φ）＝ｅｘｐ（－φ２）の初期構造、最適化過程の構造及び最適構造を示す図である。 【図8】初期構造に関する比較を行うためのＣａｓｅ１の初期構造、最適化過程の構造及び最適構造を示す図である。 【図9】初期構造に関する比較を行うためのＣａｓｅ２の初期構造、最適化過程の構造及び最適構造を示す図である。 【図10】設計問題１における複雑度係数に関する比較を示す図である。 【図11】設計問題２における複雑度係数に関する比較を示す図である。 【図12】等断面制約を付加した場合の最適化問題の設計領域と境界条件を示す図である。 【図13】等断面制約を付加しない場合及び付加した場合の最適構造を示す図である。 【図14】熱伝達問題の設計領域と境界条件を示す図である。 【図15】熱伝達問題の初期構造及び複雑度係数別の最適構造を示す図である。 【図16】内部発熱問題の設計領域及び境界条件を示す図である。 【図17】内部発熱問の複雑度係数別の最適構造を示す図である。 【発明を実施するための形態】 【0023】 次に、本発明の一実施形態について図面を参照して説明する。なお、図１は本実施形態の構造最適化装置１００の機器構成図であり、図２は構造最適化装置１００の機能構成図であり、図３は構造最適化装置１００の動作を示すフローチャートである。 【0024】 ＜装置構成＞ 本実施形態に係る構造最適化装置１００は、予め設定された設計領域内に、所定の制約条件下において所望の性能を有する構造物を創成するものであり、図１に示すように、ＣＰＵ１０１に加えて揮発メモリやＨＤＤ等の記憶装置１０２を備え、さらにマウスやキーボード等の入力手段１０３、解析モデルや演算結果の出力するためのディスプレイやプリンタ等からなる出力手段１０５を接続するための入出力インターフェイス１０４等を有した汎用乃至は専用のコンピュータである。 【0025】 そして、所定のプログラムを前記記憶装置１０２にインストールし、そのプログラムに基づいてＣＰＵ１０１や周辺機器を協働させることにより、この構造最適化装置１００は、図２の機能構成図に示すように、設計領域データ格納部１、レベルセット関数データ格納部２、境界条件データ格納部３、解析データ格納部４、レベルセット関数更新部５、演算結果出力部６等としての機能を発揮する。 【0026】 以下各部について説明する。 【0027】 設計領域データ格納部１は、構造物の設計領域（当該設計領域を要素分割する構造格子（メッシュ）情報を含む。）を示す設計領域データを格納するものである。なお、設計領域データは、例えばユーザが入力手段１０３を用いることにより入力される。 【0028】 レベルセット関数データ格納部２は、初期構造などの設計領域内にある構造を特定するためのレベルセット関数を示すレベルセット関数データを格納するものである。レベルセット関数とは、初期構造が設定された設計領域の各部が、構造を形成し、物体により占められた物体領域（物体相）、空洞を形成する空洞領域（空洞相）、又はそれら領域の境界であるかを示すものであり、物体領域を表す値及び空洞領域を表す値の間の所定値が、物体領域及び空洞領域の境界を表す。なお、レベルセット関数データは、例えばユーザが入力手段１０３を用いることにより入力される。 【0029】 境界条件データ格納部３は、前記設計領域（以下、固定設計領域ともいう。）の境界条件を示す境界条件データを格納するものである。具体的な境界条件としては、例えば設計領域の拘束条件、初期構造に作用する荷重等の外力（表面力）などである。なお、境界条件データは、例えばユーザが入力手段１０３を用いることにより入力される。 【0030】 解析データ格納部４は、最適構造を求めるときに用いる体積制約に必要な制約の条件値及び変位場の解析に必要なヤング率やポアソン比の値等の材料定数である解析データを格納している。なお、解析データは、例えばユーザが入力手段１０３を用いることにより入力される。 【0031】 レベルセット関数更新部５は、所定の制約条件下において、剛性、固有振動数等の構造物の性能を目標値に近づけるようにレベルセット関数を更新して、物体領域と空洞領域との境界を移動させるとともに、レベルセット関数の更新に伴う物体領域内のトポロジー変化（形態変化）を許容して物体領域内に新たな空洞領域（穴）を形成し、その新たな空洞領域と物体領域との境界を移動させるものである。 【0032】 具体的にレベルセット関数更新部５は、レベルセット関数を変数とする関数族、物体領域におけるエネルギー密度、空洞領域におけるエネルギー密度、及び界面エネルギー密度により示されるエネルギー汎関数を、エネルギー汎関数最小化原理に従って、前記レベルセット関数の時間発展を示す反応拡散方程式を算出し、当該反応拡散方程式を用いて前記レベルセット関数を時間発展させることにより、前記レベルセット関数を更新する。なお、レベルセット関数更新部５の具体的な機能については後述する。 【0033】 演算結果出力部６は、レベルセット関数更新部５により更新されたレベルセット関数の演算結果、つまり最適構造の形状を出力するもので、この実施形態ではディスプレイ１０４に表示する。 【0034】 ＜構造最適化装置１００の動作＞ 次にこのように構成した構造最適化装置１００の動作について図３を参照して説明する。 【0035】 まず、設計領域データ、レベルセット関数データ、境界条件データ及び解析データを入力すべく入力手段１０３をユーザが操作する。ここで、ユーザにより入力されるレベルセット関数データは、初期構造を示す初期のレベルセット関数を示すものである。 【0036】 このようにして入力された各データをデータ受付部（図示しない）が受け付け、設計領域データを設計領域データ格納部１に、レベルセット関数データをレベルセット関数データ格納部２に、境界条件データを境界条件データ格納部３に、解析データを解析データ格納部４に格納する（ステップＳ１）。 【0037】 次に、レベルセット関数更新部５が、レベルセット関数に対応する目的汎関数（剛性、固有振動数等の構造物の性能を目標値に近づけるための汎関数）と制約汎関数（制約条件を示す汎関数）とを有限要素法を用いて計算する。 【0038】 そして、レベルセット関数更新部５は、目的汎関数値が収束したか否かを判断し（ステップＳ３）、収束したならば、最適解が得られたと判断して最適化を終了し、その時のレベルセット関数を演算結果出力部６に出力する。一方、収束しない場合には、有限要素法を用いてレベルセット関数を更新する（ステップＳ４）。このとき、レベルセット関数更新部５は、制約条件（例えば体積制約）を満たすか否かを判断し（ステップＳ５）、更新後に制約条件が満たされていれば、再び、ステップＳ２に戻る。一方、制約条件が満たされていれば、後述する体積修正法を用いて、レベルセット関数を修正し（ステップＳ６）、ステップＳ２に戻る。なお、初期構造の体積が体積制約と大きく異なる場合、２００ステップ程度をかけて体積が減少するように体積修正を行う。 【0039】 以下に本実施形態に係る手法をレベルセット関数更新部５の機能とともに詳述する。 【0040】 ＜レベルセット法に基づく構造最適化問題＞ 物体により占められている領域Ω(以下、物体領域)の存在が許容される固定領域Ｄ（以下、固定設計領域）において、物体領域の構造最適化について考える。レベルセット法は、固定設計領域Ｄにおいて、レベルセット関数と呼ばれるスカラー関数φ（ｘ）を導入する。そして、そのレベルセット関数φのゼロ等位面によって物体境界を陰的に表現する方法である（図４参照）。 【0041】 本実施形態では、後述のように、フェーズフィールド法の考え方に基づいてレベルセット関数φを時間発展させるため、レベルセット関数φに上限値と下限値を設定し、次式で定義する。 【0042】 【数1】 【0043】 なお、物体領域内部において、レベルセット関数値が正の実数値であれば、同一の形状を表現することになるので、上述の制約を設けることは許容される。 【0044】 レベルセット法による形状表現を用いて、構造最適化問題を次式で定義する。 【0045】 【数2】 【0046】 ここで、Ｆは目的汎関数、Ｇは制約汎関数、Ｖｍａｘは許容される体積の上限値である。また、目的汎関数Ｆは、エネルギー汎関数であり、物体領域における自由エネルギー密度、空洞領域における自由エネルギー密度及び境界エネルギー密度により表わされる。さらに、自由エネルギー密度は、レベルセット関数を変数とする関数族Ω（φ）により与えられる。なぜなら、レベルセット関数φの値により各点の相が識別され、それにより自由エネルギー密度が与えられるからである。すなわち、自由エネルギー密度はレベルセット関数φの陽関数ではなく、物体形状Ωの陽関数である。また、レベルセット関数φから物体形状Ωへの写像について考えてみると、それは一意ではないからである。 【0047】 トポロジー最適化は、微小の空隙が至るところに散らばっているような多孔質のものや、リブが極めて短い間隔で並ぶ板構造、混合体などの一般化された構造を最適設計解として許容している。そのため、数値計算における空間的な離散化の大きさより小さい構造を表現することはできないので、構造の構造的な複雑さを示す複雑度を予め設定する。もしくは、最適化問題の緩和をしなければならない。実用上の観点から、無限に微細な構造を排除し、任意に構造の複雑度を設定することが望ましい。 【0048】 この複雑度の設定方法は、フェーズフィールド法の考え方に基づき、境界エネルギーを考慮したエネルギー汎関数を最小化させることにより、構造の複雑度（微細度）を陰的に考慮することにする。すなわち、境界エネルギー項の拡散項が微細な構造を排除する機能を有することによって、構造の複雑度が設定される。本実施形態においては、フェーズフィールド法の境界エネルギー表現に着目し、陰的に構造の複雑度の設定を行う。すなわち、式（２）、式（３）の構造最適化問題を、次式に示す境界エネルギー及び目的汎関数の和を最小化させる最適化問題に置き換えることにより、上述の問題を解決する方法を提案する。 【0049】 【数3】 【0050】 ここで、目的汎関数Ｆの第二項が境界エネルギーを表すことになる。レベルセット関数φが、フェーズフィールド法のフェーズフィールド変数と同様なプロファイルを持つことで、境界のエネルギーを表現することになるため、式（１）においてレベルセット関数φの上限値及び下限値を設定している。また、τは、境界エネルギーと目的汎関数Ｆとの比を与えるパラメータであり、複雑度係数と呼ぶことにする。 【0051】 ＜反応拡散方程式＞ 最適化問題、式（４）及び式（５）の最適解φを直接求めることは困難であるため、最適化問題を、レベルセット関数φの時間発展方程式を解く問題に帰着させる。通常、外形形状感度に基づき、次式に示す移流方程式に従って、レベルセット関数φを時間発展させることにより、形状最適化を行う。 【0052】 【数4】 【0053】 ここで、ＶＮは物体境界の移流速度の法線成分であり、最適化問題の定式化より得られる外形形状感度を、これに置き換える。この方法は、基本的には形状最適化の方法であるため、物体領域に空洞領域が創出されるようなトポロジー変化が許容されない問題を持つ。そのため、恣意的に物体領域中に穴を創出する方法が提案されているが、パラメータの設定に大きく依存し、安定的に最適解を得るのは難しい。さらには、熱・構造の連成問題、静電場・熱・構造の連成問題においては、極めて数値的に不安定であり、解の収束性が悪いことが、多くの検討結果からわかっている。そのため、レベルセット関数φの滑らかさが保障されており、かつ、物体領域に穴が創成されるようなトポロジー変化を許容した時間発展方程式を解く問題に帰着させることが望ましい。 【0054】 そこで本実施形態では、上述の問題を解決する方法として、界面拡散に基づいて界面の移流を表現するフェーズフィールド法の考え方に着目し、フェーズフィールド法と同じ考えにより、エネルギー汎関数最小化原理によって、レベルセット関数φに関する時間発展方程式である反応拡散方程式を導出する。すなわち、次式に示すように、レベルセット関数φを時間発展させる駆動力は、目的汎関数Ｆの勾配に比例するものとする。 【0055】 【数5】 【0056】 ここで、Ｋ（φ）＞０は比例定数、δＦ／δφは目的汎関数Ｆの汎関数微分を表す。上式に式（４）を代入すれば次式が得られる。 【0057】 【数6】 【0058】 ここで、Ｈ（φ）はヘビサイド関数を表す。また、形状Ω（φ）に対するレベルセット関数φの感度ｄΩ／ｄφは、物体領域Ωにおいてレベルセット関数φが一定となるようにすれば、定数Ｃとして考えればよいことがわかる。したがって、上式を次式で置き換えることができる。 【0059】 【数7】 【0060】 比例定数Ｋと複雑度係数τがパラメータであることから、定数Ｃはパラメータとして取り扱うことができ、レベルセット関数φのプロファイルが急峻になるようにＣを設定すればよい。また、境界条件については、予め物体領域境界であることが指定されている境界∂ＤＮ（以下、非設計境界という。）においては、ディリクレ境界条件、その他の境界においてはノイマン境界条件を与えることにより、固定設計領域Ｄ外部からの影響がないことを表現する。このとき、時間発展方程式系は次式となる。 【0061】 【数8】 【0062】 なお、フェーズフィールド法では、フェーズフィールド変数の関数を被積分関数とした領域積分の勾配に比例するものとしていることに対して、本実施形態では、領域変動による勾配に比例して、レベルセット関数φを時間発展させる点で異なる。その結果として、フェーズフィールド法においては、時間発展の過程において、界面が移流するのみであることに対して、本実施形態で提案する方法では、物体領域内に空洞領域が創成されるようなトポロジー変化を許容している。また、従来の方法で用いられていた移流方程式（６）は、レベルセット関数φの滑らかさが保障されていないことに対して、時間発展方程式（１０）は反応拡散方程式と呼ばれる偏微分方程式であり、拡散項が含まれることにより、レベルセット関数φの滑らかさが保障されていることを注記しておく。 【0063】 体積制約については、ラグランジュ未定乗数法により体積制約を課す方法を用いる。このとき、ラグランジュ乗数をλとすれば、時間発展方程式（９）は次式となる。 【0064】 【数9】 【0065】 ラグランジュ乗数λについては、制約条件が活性であるとき、 【0066】 【数10】 【0067】 である。したがって、式（１１）、式（１２）よりラグランジュ乗数λは次式で与えられる。 【0068】 【数11】 【0069】 また、制約条件が非活性であるときは、次式となる。 【0070】 【数12】 【0071】 以上より、ラグランジュ乗数λの値を求め、式（１０）に基づいてレベルセット関数φを初期値から更新させていき、目的汎関数Ｆが収束したときのレベルセット関数φが最適構造を与える。 【0072】 ここで、式（１０）に基づいてレベルセット関数φを時間発展させる場合、目的汎関数Ｆが単調減少することを、次式より確かめることができる。 【0073】 【数13】 【0074】 なお、フェーズフィールド法においては、エネルギー密度は識別関数の関数としているため、境界∂Ｄが移流することにより、エネルギー汎関数が単調減少する。それゆえ、時間発展の過程において、物体領域に空洞領域が創出されるようなトポロジー変化は許容されない。これに対し、本実施形態の方法では物体形状の表現を関数族Ω（φ）を導入することにより、物体形状の変化に対する勾配に比例して、レベルセット関数φが時間発展する。その結果、境界∂Ｄの移動に加えて、恣意的ではない自動的なトポロジー変化を許容していることを注記しておく。 【0075】 ＜平均コンプライアンス最小化問題＞ 線形弾性体で構成される物体領域と空洞領域の混在が許容される固定設計領域Ｄに対し、境界Γｕを完全拘束し、境界Γｔに表面力ｔ、物体領域に物体力ｂを作用させる構造問題を考える。ただし、境界Γｕは、固定設計領域境界∂Ｄに固定されているものとする。このとき、体積制約下で、平均コンプライアンスを最小化させる構造最適化問題は次式のように記述される。 【0076】 【数14】 【0077】 ここで、上記式中の各表記は次式で定義される。 【0078】 【数15】 【0079】 さらに、εはひずみテンソル、Ｆは弾性テンソル、Ｖｍａｘは体積の制約値を表し、Ｕは以下の式にて定義される変位関数空間である。 【0080】 【数16】 【0081】 なお、表面力が作用する境界Γｔは常に物体境界である必要があるため、境界Γｔを非設計境界∂ＤＮとする。 【0082】 次に、上の最適化問題に関するＫＫＴ条件を導くとともに、その結果から、レベルセット関数φの更新に必要な関数ｆ（ｘ）を与える。上述の定式化より、ラグラジュアンＦ’（Ω）は、ラグランジュ乗数λ、随伴変位場ｖを用いて以下のように記述される。 【0083】 【数17】 【0084】 式（２３）を用いてＫＫＴ条件を導けば、次式となる。 【0085】 【数18】 【0086】 このＫＫＴ条件を満たすレベルセット関数φが最適解（最適構造）の候補となる。しかし、これらを満たすレベルセット関数φを直接求めることは困難であるため、適当な初期構造となるレベルセット関数φを与え、式（１１）を用いてレベルセット関数φを更新することにより、ラグラジュアンＦ’（Ω）及び境界エネルギーの和を減少させる。 【0087】 以下、式（１１）を用いてレベルセット関数φを更新する際に必要となるラグラジュアンの目的汎関数を与える関数ｆ’（ｘ）を導出する。ここで、随伴変位場ｖを次式と定義することにより、最適化問題を自己随伴問題とする。 【0088】 【数19】 【0089】 式（２３）に式（２５）を代入し、平衡方程式（２４）を用いれば、ラグラジュアンＦ’（Ω）は次式となる。 【0090】 【数20】 【0091】 したがって、目的汎関数を与える被積分関数ｆ’（ｘ）は次式となる。 【0092】 【数21】 【0093】 ＜反応拡散方程式の数値解法＞ 次式に示すように、式（１０）を時間方向に対して、差分法による離散化を行う。 【0094】 【数22】 【0095】 ここで、Δｔは時間差分である。なお、式（１５）の第一式の拡散項τ∇２φを更新後の値とすることにより、更新後のレベルセット関数φに対して、境界エネルギーによる拡散効果を考えることができる。したがって、全更新ステップに対して拡散効果を考慮していることになる。 【0096】 次に、有限要素法を用いて空間方向に離散化するために、式（２８）の弱形式を導出すると、次式となる。 【0097】 【数23】 【0098】 ここでΦ～は、次式にて定義されるものであり、レベルセット関数φに関する関数空間である。 【0099】 【数24】 【0100】 式（２９）を有限要素法を用いて離散化すれば、次式となる。 【0101】 【数25】 【0102】 ここで、Φ（ｔ）は、時刻ｔにおける各節点のレベルセット関数値により構成されるベクトルであり、マトリックスＴ及びベクトルＹは次式で与えられる。 【0103】 【数26】 【0104】 本実施形態では、陰解法によりレベルセット関数φを更新するため、ＣＦＬ条件による時間増分Δｔの制限がない。その結果、時間増分Δｔを増加させていくことにより、時刻ｔが増加するとともに境界移動速度が徐々に遅くなる問題を回避することが可能となる。なお、最初の更新ステップ２０回程度においては、レベルセット関数φの変動の最大値が１程度となるように時間増分Δｔを設定し、徐々に時間増分を増大させることにする。 【0105】 ＜オイラー座標系に基づく変位場の近似解析法＞ レベルセット関数φを更新させるためには、物体領域において変位場の解析をする必要がある。ラグランジュ座標系に基づいて、物体領域の解析を行う場合、各更新ステップ毎に物体形状が異なるため、逐次メッシュを生成する必要があるうえ、離散的に得られた変位場を固定設計領域Ｄに写像しなければならない問題を持つ。さらに、最適化過程において、変位固定部と連結していない物体領域が生じた場合、その物体領域では剛体モードを持つ問題がある。そこで、本実施形態では、オイラー座標系に基づき、空洞領域では小さな縦弾性係数を持つ物体とみなし、境界近傍では縦弾性係数が空間的に滑らかに変化するものとして近似することにする。すわなち、平衡方程式（２４）を次式に置き換えることにより、物体形状を陽に抽出することなく変位場の解析を行う。 【0106】 【数27】 【0107】 ここでＨｅ（φ）は次式で記述される関数である。 【0108】 【数28】 【0109】 ここで、ｄは物体領域の縦弾性係数に対する空洞領域の縦弾性係数の相対値（比）、ｔは材料定数の遷移幅を示す値であり、どちらも十分に小さい正の値として設定する。また、体積制約関数Ｇ（Ω）に関しても、次式の近似式を用いて求めることにする。 【0110】 【数29】 【0111】 ＜実装法の詳細について＞ パラメータＫ（φ）は、レベルセット関数φの関数であるため、境界∂Ｄの移動と構造の形態変化の優先度の度合を決めるパラメータ（易動度）であることがわかる。すなわち、境界近傍におけるＫ（φ）を大きく設定することは境界∂Ｄの移動を、トポロジー変化（形態変化）に対して優先させることに相当する。多くの計算結果から、パラメータＫ（φ）の設定によらず、同一の最適構造が得られることがわかった。この点に関して、Ｋ１（φ）＝１（一定）と、φ＝０を中心としたガウス分布Ｋ２（φ）＝ｅｘｐ（－φ２）に対して、数値例においてその結果を示す。なお、Ｋ１（φ）の場合は、φによらず一定であるため、境界∂Ｄの移動とトポロジー変化との優先度合が同程度、Ｋ２（φ）の場合はφ＝０近傍の値が大きいため、境界∂Ｄの移動を優先させることになる。 【0112】 また、ラグランジュ未定乗数法により体積制約を考慮する場合、各更新ステップにおける体積制約に対する数値誤差は微小ではあるものの、繰り返し計算を行うことにより、その数値誤差が累積される問題を持つ。そのため、体積制約を満たさない場合、レベルセット関数φに修正を加えることにより、体積修正を行う必要がある。本実施形態では、次式を満たす微小値Δφ（ｘ）の最小値を二分法を用いて求め、新たにφ（ｘ）＋Δφ（ｘ）をレベルセット関数φの値とする。 【0113】 【数30】 【0114】 ＜数値例＞ 次に、本実施形態の構造最適化装置１００を用いて、上述した方法の妥当性を検討する。 【0115】 図５（ａ）に設計問題１の設計領域Ｄと境界条件を示す。固定設計領域Ｄは８×１０－２ｍ×６×１０－２ｍの長方形領域とし、設計領域Ｄを要素長５×１０－４ｍの構造格子で要素分割した。また、許容される体積の最大値Ｖｍａｘは、固定設計領域Ｄの４０％とした。境界条件は、左端を完全変位拘束し、右端中央に下向きの表面力を作用させる。 【0116】 図５（ｂ）に設計問題２の設計領域Ｄと境界条件を示す。固定設計領域Ｄは８×１０－２ｍ×６×１０－２ｍの長方形領域とし、設計領域Ｄを要素長５×１０－４ｍの構造格子で要素分割した。また、許容される体積の最大値Ｖｍａｘは、固定設計領域Ｄの５０％とした。境界条件は、下端左側及び右側を完全変位拘束し、下端中央に下向きに表面力を作用させる。 【0117】 変位場の解析には、設計問題１及び設計問題２のいずれの場合においても、４節点のアイソパラメトリック４角形平面応力要素、レベルセット関数φの更新には、４節点のアイソパラメトリック４角形要素を用いた。解析モデルの材料は、いずれの場合も、等方性材料を想定し、材料定数は縦弾性係数を２１０ＧＰａ、ポアソン比を０．３とした。さらに式（３５）において定義された空洞領域の縦弾性係数の物体領域の縦弾性係数に対する比ｄを１×１０－６、陰的な遷移幅の設定値ｔを０．１と設定した。 【0118】 ＜パラメータＫ（φ）に関する比較＞ 最初に、パラメータＫ（φ）に対する比較を行った。 【0119】 初期構造を固定設計領域Ｄの上部２／３を占めている構造とし、複雑度係数τを０．０７、Ｃ＝１／ｍａｘ（ｆ（ｘ））と設定し、設計問題１に対して最適化を行った。 【0120】 Ｋ１（φ）＝１とした場合の初期構造、最適化過程の構造及び最適構造を図６に示す。同様に、Ｋ２（φ）＝ｅｘｐ（－φ２）とした場合の初期構造、最適化過程の構造及び最適構造を図７に示す。図６及び図７中、黒色部分が物体領域であり灰色部分が空洞領域である。図６及び図７の結果を比較すると、最適化過程の様子が異なり、境界∂Ｄの移動を優先させている図７の結果の方が、速く境界∂Ｄが移動している様子が分かる。また、得られた最適構造は同一の構造であり、明瞭かつ物理的に妥当な結果が得られていることが分かる。 【0121】 ＜初期構造に関する比較＞ 次に、異なる初期構造を設定し、設計問題１に対して最適化を行った。 【0122】 一方（Ｃａｓｅ１）の初期構造は、図８の（ａ）に示すように、固定設計領域Ｄの全領域が物体領域に占められている構造とし、他方（Ｃａｓｅ２）の初期構造は、図９の（ａ）に示すように、Ｃａｓｅ１の初期構造において穴が２つ空いた構造とした。また、Ｃａｓｅ１及びＣａｓｅ２のいずれも、複雑度係数τは０．０７、その他のパラメータはそれぞれ、Ｋ１（φ）＝１、Ｃ＝１／ｍａｘ（ｆ（ｘ））とした。図８及び図９の（ｂ）及び（ｃ）に最適化過程、（ｄ）に最適構造を示す。得られた結果は同一の最適構造であり、物理的にも妥当な構造であることが分かる。 【0123】 ＜複雑度係数に関する比較＞ 次に、複雑度係数τの設定値を変化させて、設計問題１及び設計問題２に対して最適化を行った。 【0124】 ここで、固定設計領域Ｄの全領域が物体領域により占められた構造を初期構造とし、Ｋ１（φ）＝１、Ｃ＝１／ｍａｘ（ｆ（ｘ））とした。また、設計問題１において、複雑度係数τは、Ｃａｓｅ１の設定値を０．５とし、Ｃａｓｅ２の設定値を０．０５とし、Ｃａｓｅ３の設定値を０．０３とした。一方、設計問題２において、Ｃａｓｅ１の設定値を０．０１とし、Ｃａｓｅ２の設定値を０．００５とし、Ｃａｓｅ３の設定値を０．０００１とした。なお、いずれの設計問題においても、Ｃａｓｅ２の方が複雑度係数τの値が小さいので、より複雑な（詳細な）構造を最適解として許容していることになる。 【0125】 図１０に設計問題１におけるＣａｓｅ１、Ｃａｓｅ２及びＣａｓｅ３の最適構造を示し、図１１に設計問題２におけるＣａｓｅ１、Ｃａｓｅ２及びＣａｓｅ３の最適構造を示す。これらの図により、いずれの場合も物理的に妥当な最適構造が得られていることが分かった。さらに、複雑度係数τの大きさによって、最適構造の複雑度が異なることが分かる。したがって、本実施形態で提案する方法は、複雑度係数の設定により、構造の複雑度を陰的に設定可能であり、最適構造の複雑度を定性的に考慮可能な方法であることを確かめることができた。 【0126】 ＜製造上の制約（等断面制約）を付加した場合の最適化＞ 図１２に等断面制約を設けた場合の最適化における設計領域及び境界条件を示す。固定設計領域Ｄは、２．０ｍ×０．８ｍ×０．１５ｍの直方体領域とし、その下部に２．０ｍ×０．８ｍ×０．０５ｍの直方体形状をなす非設計領域を設定した。境界条件は、長手方向両端部を完全変位拘束し、非設計領域の下面全面に下向きに表面力を作用させた。 【0127】 このとき、最適構造がＸ３方向に等断面となるように幾何学的制約を付加した。具体的には、複雑度係数τを異方性とし、等断面としたい方向の複雑度係数をその他の複雑度係数よりも十分に大きく設定する。つまり、各軸方向における複雑度係数τ（Ｘ１）、τ（Ｘ２）、τ（Ｘ３）において、τ（Ｘ３）をτ（Ｘ１）及びτ（Ｘ２）に比べて十分大きく設定する（τ（Ｘ３）≫τ（Ｘ１）、τ（Ｘ２））。 【0128】 図１３に（Ａ）「等断面制約なし」の場合の最適構造、及び（Ｂ）「等断面制約あり」の場合の最適構造を示す。これらの図より、複雑度係数τを異方性（τ（Ｘ３）≫τ（Ｘ１）、τ（Ｘ２））とすることによって、等断面制約を付加した最適構造を得ることが分かる。 【0129】 ＜熱拡散問題＞ 上記の比較例等においては、剛性最大化問題に適用した場合について説明したが、熱拡散問題に適用することもできる。熱拡散問題としては、熱伝導問題、熱伝達問題及び内部発熱問題が考えられる。 【0130】 以下に具体例として、熱伝達問題及び内部発熱問題について説明する。 【0131】 ＜熱拡散最大化問題の定式化＞ 線形熱伝導体で構成される物体領域と空洞領域で構成される固定設計領域Ｄに対し、境界Γｔにおいて温度Ｔ０で温度規定、境界Γｈにおいて熱伝達係数ｈの熱伝達境界、境界Γｑにおいて熱流束ｑの熱流束境界、固定設計領域Ｄに対して内部発熱Ｑが与えられている熱拡散最大化問題について考える。但し、境界Γｔ及び境界Γｑは、固定設計領域Ｄの境界∂Ｄ上に設定しているものとする。また、境界Γｈは固定設計領域内部において、物体領域の境界上で設定され、設計変数であるレベルセット関数の値に依存して決定される設計変数依存性の境界条件となる。このとき、熱拡散最大化問題は、次式に示す全ポテンシャルエネルギー最大化問題として定式化される。 【0132】 【数31】 【0133】 ここで、上式中の各表記は次式で定義される。 【0134】 【数32】 【0135】 さらに、κは熱伝導テンソル、Ｔａｍｐは周囲温度を表し、Ｕｔは以下の式にて定義される温度関数空間である。 【0136】 【数33】 【0137】 なお、式（４３）におうては、全ポテンシャルエネルギーにマイナス符号を付加し、目的関数の最小化問題として定式化している。 【0138】 ＜熱拡散問題の数値例（熱伝達問題）＞ 図１４に設計領域と境界条件を示す。固定設計領域は１×１０－２ｍ×１×１０－２ｍの正方形領域とし、設計領域を要素長５．０×１０－５ｍの構造格子で要素分割した。図に示すように、左端下方と下端左方をＴ０＝５０℃で温度設定し、その他の境界では断熱境界を与える。また、固定設計領域内部において、熱伝達率ｈ＝０．１Ｗ／（ｍ２Ｋ）、周囲温度Ｔａｍｐ＝２５℃の下で、熱伝達境界を与える。許容される体積の最大値Ｖｍａｘは固定設計領域の６０％とした。また、熱伝達境界を定義するためには、初期構造に物体領域と空洞領域との境界が少なくとも一箇所必要となるため、図１５（ａ）に示すように、初期構造を左下の角を中心とした１／４円弧とした。 【0139】 図１５の（ｂ）～（ｅ）に、複雑度係数τを５×１０－４、１×１０－４、５×１０－５、１×１０－５に設定した場合に得られた最適構造を示す。図１５より複雑度係数の設定値により、フィン形状の複雑さが変化することが分かった。また、いずれの場合も物理的に妥当な最適構造が得られていることが分かる。 【0140】 ＜熱拡散問題の数値例（内部発熱問題）＞ 図１６に内部発熱問題の設計領域及び境界条件を示す。固定設計領域は１×１０－２ｍ×１×１０－２ｍの正方形領域とし、設計領域を要素長２．５×１０－５ｍの構造格子で要素分割した。図１６に示すように、上端中央をＴ０＝２５℃で温度設定し、その他の境界では断熱境界を与える。また、固定設計領域内部において、熱量Ｑ＝１．０×１０－７Ｗ／ｍ３を与え、許容される体積の最大値Ｖｍａｘは固定設計領域の４０％とした。初期構造は、固定設計領域の全領域が物体に占められた構造とした。ここでは、複雑度係数と最適構造との関係を調べるために、複雑度係数を変化させて最適化を行った。 【0141】 図１７の（ａ）～（ｄ）に、複雑度係数τを５×１０－５、１×１０－５、５×１０－６、１×１０－６に設定した場合に得られた最適構造を示す。なお、設計領域は左右対称であるため、右半分のみ最適化解析を行った。これらの結果により、この設計問題においても、複雑度係数の設定により、最適構造の複雑さが変化していることが分かる。また、いずれの場合も物理的に妥当な最適構造が得られていることが分かる。 【0142】 ＜本実施形態の効果＞ このように構成した本実施形態に係る構造最適化装置１００によれば、物体領域Ωを表す値及び空洞領域を表す値の間の所定値が、物体領域Ω及び空洞領域の境界∂Ｄを示すレベルセット関数φを用いていることにより、最適構造の明確な形状表現を可能にすることができる。また、物体領域Ω及び空洞領域間の境界∂Ｄを移動させるとともに、物体領域Ω内でのトポロジー変化を許容して新たに生じる境界∂Ｄを移動させるようにレベルセット関数φを更新することにより、物体領域Ωでのトポロジー変化を許容する等の自由度の高い構造最適化を可能にすることができる。 【0143】 なお、本発明は前記実施形態に限られるものではなく、その趣旨を逸脱しない範囲で種々の変形が可能であることは言うまでもない。 【0144】 例えば、前記実施形態の構造最適化装置は、剛性最大化問題又は熱拡散最大化問題に適用したものであったが、その他固有振動数最大化問題など、種々の構造問題に適用することができる。 【産業上の利用可能性】 【0145】 本発明を適用することにより、最適構造の明確な形状表現を可能にし、物体領域Ωでのトポロジー変化を許容する等の自由度の高い構造最適化を可能にすることができる。 【符号の説明】 【0146】 １００・・・構造最適化装置 Ｄ ・・・設計領域 Ω ・・・物体領域 ∂Ω ・・・領域の境界 φ ・・・レベルセット関数 １ ・・・設計領域データ格納部 ２ ・・・レベルセット関数データ格納部 ５ ・・・レベルセット関数更新部
図面	【図1】 【図2】 【図3】 【図4】 【図5】 【図6】 【図7】 【図8】 【図9】 【図10】 【図11】 【図12】 【図13】 【図14】 【図15】 【図16】 【図17】