形状最適化方法、形状最適化装置、及び、プログラム
【課題】形状最適化における計算の処理負担を低減させる。
【解決手段】形状最適化装置１０は、まず、初期形状Ω０に対応するφの初期値を決定する（ステップＳ２０）。次に、有限要素法により変位ｕ、固有モードφｋ、随伴状態ｐを導出する（ステップＳ２１）。続いて、目標関数を導出する（ステップＳ２２）。そして、目標関数の感度を導出し（ステップＳ２３）、また二重井戸型ポテンシャルｆ（φ）を導出して、反応拡散方程式を有限体積法により解く（ステップＳ２４）。そして、以上の処理を目標関数が収束するまで繰り返す。
【解決手段】形状最適化装置１０は、まず、初期形状Ω０に対応するφの初期値を決定する（ステップＳ２０）。次に、有限要素法により変位ｕ、固有モードφｋ、随伴状態ｐを導出する（ステップＳ２１）。続いて、目標関数を導出する（ステップＳ２２）。そして、目標関数の感度を導出し（ステップＳ２３）、また二重井戸型ポテンシャルｆ（φ）を導出して、反応拡散方程式を有限体積法により解く（ステップＳ２４）。そして、以上の処理を目標関数が収束するまで繰り返す。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、形状最適化方法、形状最適化装置、及び、プログラムに関する。
【背景技術】
【０００２】
偏微分方程式で特定の物理問題が記述される領域において、特定の物理量の最適化を目的とし、その領域形状を変化させる構造最適化問題を考える。このような問題は、機械製品の性能向上を目指す上で非常に重要な技術であり、様々な解法が提案されている。これらの解法の一般的な区分として、対象となる領域の境界を移動することにより最適化を行う形状最適化がある。また、対象の形状のみならず形態も最適化可能なトポロジー最適化が挙げられる。
【０００３】
しかし、一般的な構造最適化の方法においては二つの問題点がある。一つは領域を支配する偏微分方程式の数値計算を行う際に、有限要素法を用いるのが一般的であるが、領域の形状が更新される度に有限要素メッシュを再生成する必要があり、計算の処理負担が非常に大きくなることである。
【０００４】
もう一つは、スプライン関数等で領域を陽的に表現するため、境界同士の接触による穴の消失等の境界のトポロジー変化を表現できないことである。
【０００５】
トポロジー最適化は、多くの数学者によって基礎理論が構築されている。この方法は領域の形状のみならず形態も最適化可能であり、様々な工学的問題に適用され、成功を収めている。しかし、基本的には最適構造が材料密度分布で表されるため、そこから明確な構造を得ることが難しく、圧力や熱放射等の境界形状に依存するような境界条件の設定が困難という問題点がある。
【０００６】
近年、これらの問題を解消する方法として、レベルセット法を用いた形状最適化が提案されている（例えば非特許文献１）。レベルセット法は平均曲率流の数値計算法として提案され、現在では数値流体力学、画像処理等の様々な分野における移動境界問題に適用されている。レベルセット法では、領域の境界が陰関数を用いて表現されるため、オイラー型のメッシュを用いて数値計算が可能であり最適化の過程で計算量が一定である。また、境界同士の接触、穴の消失といった形態の変化が考慮可能である等の利点がある。
【０００７】
ここで、レベルセット法に基づく形状最適化における計算処理について簡単に説明する。図１３は、レベルセット法に基づく形状最適化処理を説明するためのフローチャートである。図１３に示すように、先ず、初期形状、各種初期値を決定する（ステップＳ１０）。そして、決定された初期形状について、有限要素法を用いて場の問題を解く（ステップＳ１１）。これにより、初期形状の変位や応力などを導出する。そして、その結果に基づいて目標関数を導出する（ステップＳ１２）。さらに、感度νを導出する（ステップＳ１３）。導出されたνを有限要素法を用いて滑らかにする（ステップＳ１４）。次に、レベルセット関数の初期化を行う（ステップＳ１５）。続いて、ハミルトン・ヤコビ方程式を解く（ステップＳ１６）。そして、計算値が収束したか否かを判定し（ステップＳ１７）、このようなステップＳ１１〜Ｓ１６を計算値が収束するまで繰り返す。
【０００８】
【非特許文献１】G．Allaire, F．Jouve， A．Toader， Structural optimization using sensitivity analysis and a level-set method、 Journal of Computational Physics 194 (2004) pp．363-393
【発明の開示】
【発明が解決しようとする課題】
【０００９】
しかし、レベルセット法における発展方程式であるハミルトン・ヤコビ方程式は移流方程式であり、数値計算上不安定になり易い。ハミルトン・ヤコビ方程式において安定した解を得るためには、最適化の繰り返し計算が必要となり、さらに、必ずしも安定した解が得られるわけではなかった。そのため、形状最適化にレベルセット法を適用しても依然として計算の処理負担は大きかった。
【００１０】
本発明は、上記実情に鑑みてなされたものであって、形状最適化における計算の処理負担を低減させることを目的とする。
【課題を解決するための手段】
【００１１】
上記目的を達成するため、本発明の第１の観点に係る形状最適化方法は、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化方法であって、
コンピュータが、
有限要素法により前記形状の物理量を導出し、
前記導出した物理量に基づいて目標関数を導出し、
前記導出した目標関数を最小化するために用いる当該目標関数の感度を導出し、
前記導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義し、
前記反応拡散方程式を有限体積法で解くことにより、前記導出した目標関数を最小化する
最小化処理を実行することを特徴とする。
【００１２】
また、前記形状が収束するまで、前記最小化処理を繰り返してもよい。また、前記物理量が、線形弾性場あるいは材料非線形性及び幾何学的非線形性を有する連続体の弾性場における変位あるいは固有振動数、熱伝導場における定常温度、ストークス流れ場の定常流速及び圧力、ポテンシャル流れ場定常流速及び圧力、音場における定常音圧、電磁場問題における電界のいずれかであってもよい。
【００１３】
本発明の第２の観点に係る形状最適化装置は、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置であって、
有限要素法により前記形状の物理量を導出する物理量導出手段と、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段と、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段と、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段と、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段と、
を備えることを特徴とする。
【００１４】
本発明の第３の観点に係るプログラムは、コンピュータを、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置として機能させるためのプログラムであって、
コンピュータを、
有限要素法により前記形状の物理量を導出する物理量導出手段、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段、
として機能させることを特徴とする。
【発明の効果】
【００１５】
本発明によれば、レベルセット法を用いた形状最適化計算に比べて、有限要素法の計算回数を減らすことができるので、形状最適化計算の処理負担を低減できる。
【発明を実施するための最良の形態】
【００１６】
本発明の実施の形態にかかる形状最適化方法を実行する形状最適化装置１０について図面を参照しながら説明する。本発明の形状最適化方法はフェーズフィールド法を用いた形状最適化方法である。より具体的には、フェーズフィールド法を用い、オイラー型メッシュでの計算が可能で、形態変化にも対応可能であり、かつ数値計算上安定した形状最適化方法である。
【００１７】
図１は、形状最適化装置１０の構成を示すブロック図である。形状最適化装置１０は、ＣＰＵ（Central Processing Unit）１１、入力装置１２、記憶装置１３、表示制御部１４から構成される。形状最適化装置１０は、例えば、汎用または専用のコンピュータである。
【００１８】
ＣＰＵ１１は、形状最適化装置１０全体の動作を制御する。また、ＣＰＵ１１は、記憶装置１３の所定領域に格納されるプログラムを読み込むことにより、各種の処理を実行する。特に本発明においては、ＣＰＵ１１は、フェーズフィールド法を用いた形状最適化処理を実行する。
【００１９】
入力装置１２は、キーボードやマウスなどから構成され、ユーザの操作に応じた操作信号をＣＰＵ１１に入力する。
【００２０】
記憶装置１３は、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、ハードディスク装置、その他の記憶媒体などから構成される。記憶装置１３は、フェーズフィールド法を用いた形状最適化処理のプログラムや必要なデータを記憶する。また、記憶装置１３は、ＣＰＵ１１のワークエリアとして機能する。
【００２１】
表示制御部１４は、ＣＰＵ１１から伝送された表示信号をディスプレイ２０で表示可能な形式に復号し、復号した信号をディスプレイ２０に入力する。
【００２２】
ディスプレイ２０は、ＬＣＤ（Liquid Crystal Display）などから構成され、形状最適化装置１０から伝送された信号に基づいて各種画面を表示する。例えば、本発明の形状最適化装置１０により最適化される形状のイメージなどを表示する。
【００２３】
以下フェーズフィールド法、及び、形状最適化装置１０が実行するフェーズフィールド法を用いた形状最適化処理について説明する。
【００２４】
フェーズフィールド法は、相転移における移動境界問題をシミュレートするための方法として開発された。フェーズフィールド法においては、領域全体に相の状態を表すフェーズフィールド関数を設定し、異なる相にそれぞれ対応する値を設定する。これにより、領域内のある箇所における相の状態が関数の値により決定される。異なる相間の境界では、フェーズフィールド関数は滑らかに補間され、その領域はインターフェイスと呼ばれる。領域内に二つの相が存在するとき、フェーズフィールド関数をφ（ｘ）とし、二つの相にそれぞれφ＝０、φ＝１とした場合の例を図２に示す。物理的な意味としては、フェーズフィールド関数は領域内の各箇所における相の平均化された値を表す。このフェーズフィールド関数を用いて、二つの相に関するVan der Waalsの自由エネルギーは以下の式（１）ように表される。
【００２５】
式（１）
【数１】
【００２６】
式（１）の第一項は平均場理論における干渉エネルギーを表し、第二項は二重井戸型ポテンシャルを表す。図２のように二つの相にそれぞれφ＝０、φ＝１が対応する場合は、井戸型ポテンシャルはｆ’（０）＝ｆ’（１）＝０という性質を持つ。この二重井戸型ポテンシャルは、それぞれの相において、自由エネルギーが極小値をとることを意味する。しかし、異なる相間で、フェーズフィールド関数が急激に変化し勾配が大きくなることは、第一項の値の増加に繋がるため、その間は０＜φ＜１の値をとる微小なインターフェイス領域が存在し、関数は滑らかに補間されることになる。
【００２７】
次に、このフェーズフィールド関数の時間変化を決定する発展方程式について述べる。最初に、フェーズフィールド関数が領域内でその総和が保存される量である場合を考える。この保存系のフェーズフィールド関数を改めてＣと表すと、領域内で以下の連続の式が成り立つ。
【００２８】
式（２）
【数２】
【００２９】
ここで、ｊは流束ベクトルである。さらに、閉じた系での自由エネルギーが増加しないと仮定すると、以下の式が成り立つ。
【００３０】
式（３）
【数３】
【００３１】
式（２）を式（３）に代入して、流入のない境界条件を適用すると、次式が得られる。
【００３２】
式（４）
【数４】
【００３３】
これが非正であるためには、ｊがＦに関して線形であると仮定すると、式（５）が得られる。
【００３４】
式（５）
【数５】
【００３５】
ただし、Ｍ（Ｃ）＞０とする。これを式（４）に代入すると、式（６）が得られる。
【００３６】
式（６）
【数６】
【００３７】
説明の簡単のためＭを定数とし、Ｆに式（１）を用いると、式（７）が得られる。
【００３８】
式（７）
【数７】
【００３９】
式（７）は、Ｃａｈｎ−Ｈｉｌｌｉａｒｄ方程式と呼ばれ、保存場に関するフェーズフィールド関数の発展方程式として用いられる。
【００４０】
続いて、フェーズフィールド関数が領域内でその総和が保存されない量である場合について考える。フェーズフィールド関数をφとすると、式（３）は非保存の場合についても成り立つので、以下の式が成り立つ。
【００４１】
式（８）
【数８】
【００４２】
これが任意のφで満たされるためには、下記の式が満たされなければならない。
【００４３】
式（９）
【数９】
【００４４】
式（９）はＡｌｌｅｎ−Ｃａｈｎ方程式と呼ばれ、非保存場のフェーズフィールド関数の発展方程式として用いられる。
【００４５】
続いて、フェーズフィールド法とペリメータ最小化問題について述べる。ここでは、Ｃａｈｎ−Ｈｉｌｌｉａｒｄ方程式で記述される系のように、領域内で秩序変数が保存される場合の解の存在について述べる。先の例で挙げたように、フェーズフィールド関数φ（Ω；［α，β］）のαとβとがそれぞれ異なる相に対応する場合に、以下の最小化問題を考える。
【００４６】
式（１０）
【数１０】
【００４７】
ここで、Ｗ（φ）は図３（ａ）のような井戸型ポテンシャルに、アフィン変換ｃ１φ＋ｃ２を加えた、極小値が等しい二重井戸型ポテンシャルである。この変換は保存系においては定数になるため、式（１０）の最小化問題に影響を与えない。
【００４８】
式（１１）
【数１１】
ただし、εは０に近い微小な値とする。
【００４９】
この問題の解は以下のように表される。
【００５０】
式（１２）
【数１２】
【００５１】
ここで、φ（Ω；｛α，β｝）は領域Ωで最小の境界Ｓを有する関数であり、φ
１はφ１（∞）＝０となる関数で、図４に示すような異なる相間でのインターフェイスを与える。式（１０）において、第一項はφをαとβに近づける効果があり、第二項は不要なインターフェイス領域に制約を与える効果がある。式（１２）の解を式（１０）に代入して考えると、両項がε−１のオーダーになるため、ε→０のとき、式（１０）は以下に示すインターフェイス最小化問題と解釈可能である。
【００５２】
式（１３）
【数１３】
【００５３】
ここで、Ｐｅｒ（Ａ，Ω）はＡのΩにおけるペリメータを表す。これにより、φは［α，β］の値をとる滑らかな関数で定義されているのに対し、ε→０のとき式（１０）の解は結果として｛α，β｝の不連続な値となる。
【００５４】
次に、フェーズフィールド法における界面移動について説明する。式（９）のＡｌｌｅｎ−Ｃａｈｎ方程式に従い、発展する系においては、二重井戸型ポテンシャルの効果により、十分時間がたてばインターフェイス領域が形成され、それが法線方向に次式で現される速度で移動することが知られている。
【００５５】
式（１４）
【数１４】
【００５６】
式（１４）において、αは二重井戸型ポテンシャルによって定まる定数であり、κはインターフェイスの曲率を表す。二重井戸型ポテンシャルの二つの極小値が等しい場合には、αは０となり、式（１４）は完全な平均曲率運動を表す。極小値間に差がある場合には、インターフェイス領域におけるポテンシャルの値が小さくなるような方向に速度αが定まる。
【００５７】
本発明では、式（１４）のインターフェイスの放線方向に関する性質を用いて、フェーズフィールド法を用いた形状最適化方法の定式化を提供する。
【００５８】
（フェーズフィールド関数を用いた領域表現）
ここでは、Ａｌｌｅｎ−Ｃａｈｎ方程式に支配される系において、インターフェイスの移動が式（１４）で表されることを利用し、フェーズフィールド法を用いたオイラー型形状最適化法を提案する。最初に、ｄ次元実数空間に属する、体積が一定値Ｖとなる領域Ωの最適化問題を考える。このとき、図５（ａ）に示すように、Ωを包括するような設計領域Ｄを考え、その中の許容される解の集合Ｕａｄを以下のように定義する。
【００５９】
式（１５）
【数１５】
【００６０】
このとき、領域Ωの形状最適化問題は以下のように定義される。
【００６１】
式（１６）
【数１６】
【００６２】
ここで、Ｊ（Ω）は領域Ωで定義される任意の目標関数である。ここでは、この問題の数値解法として、設計領域Ｄ全体にフェーズフィールド関数φを定義し、領域Ωを陰的に表すことを考える。即ち、図５（ｂ）に示すように、ΩＡ⊂ΩとなるΩＡと、ΩＢ⊂Ｄ＼ΩとなるΩＢを定義し、さらに、領域ＤのうちΩＡとΩＢの間に、微小な厚みを持つ領域ξが存在すると仮定し、ξ⊂Ｄ＼（ΩＡ∪ΩＢ）を定義する。即ち、以下の関係が成り立つと仮定する。
【００６３】
式（１７）
【数１７】
【００６４】
このとき、これらの領域ΩＡとΩＢとはフェーズフィールド関数を用いて以下のように表される。
【００６５】
式（１８）
【数１８】
【００６６】
以上の定義において、境界∂Ωはξ内に存在することがわかるが、その位置は不定である。即ち、式（１８）によって完全に領域Ωが厳密に表されるわけではない。しかし、領域ξは微小な厚みを持つ領域であり、その厚みが十分小さい場合は、と見なすことができる。このとき、∂Ω≒ξとなり、領域Ωが式（１８）により、φで近似的に表されることになる。
【００６７】
式（１８）のように、［０，１］範囲で定義される関数を用いた領域表現方法は、均質化法に基づくトポロジー最適化でも用いられているが、その物理的意味は全く異なる。即ち、均質化法に基づくトポロジー最適化においては、材料密度分布関数で表された領域の最適化問題は、特性関数で表される領域の最適化問題の弱収束解と厳密にその物理的意味が解釈可能である。しかし、通常、式（１）で表されるように関数の勾配を目標関数に含むような場合には均質化法を適用することができない。即ち、関数φの物理的な意味は不明瞭であるが、ここでは領域表現のための陰関数としてのみ定義する。
【００６８】
（目標関数について）
領域Ｄにおける、線形弾性体問題を考える。領域ΩＡが弾性テンソルＡを有する材料で構成され、領域ΩＢを非常に弱い材料が占めていると仮定し、領域Ωにおける弾性テンソルＡ＊を以下のように定義する。
【００６９】
式（１９）
【数１９】
【００７０】
ここで、領域ξにおいては、材料の存在が明確に定義できないため、弾性テンソルの厳密な導出が不可能である。さらに、領域Ｄの境界∂Ｄのうち、ディリクレ境界条件が設定される境界をΓＤ、ノイマン境界条件が設定される境界をΓＮとする。今、領域Ｄが境界∂Ｄで完全固定され、領域Ｄに作用する体積力が０であるとし、境界ΓＮに分布荷重ｇが作用するものとすると、領域Ｄにおける線形弾性問題は弱形式で以下のように表される。
【００７１】
式（２０）
【数２０】
式（２１）
【数２１】
式（２２）
【数２２】
【００７２】
ここで、ｕは変位ベクトル、ｖはテスト関数、ｅはひずみテンソルを表す。
【００７３】
本発明では、線形弾性問題として、剛性最大化問題と、特定箇所の変位を目標値に近づけるメカニズム創生問題を考える。コンプライアンス最小化問題においては、以下で表される平均コンプライアンスを目標関数とする。
【００７４】
式（２３）
【数２３】
【００７５】
コンプライアントメカニズム創生問題では、以下で表される、変位と変位の目標値との最小自乗誤差を目標関数とする。
【００７６】
式（２４）
【数２４】
【００７７】
さらに、領域ΩＡを構成する材料の質量密度をρとし、領域ΩＢを占める材料の質量密度を非常に小さな値に設定すると、領域Ｄにおける質量密度ρ−は以下のように表される。
【００７８】
式（２５）
【数２５】
【００７９】
このとき、領域Ｄが境界∂Ｄで完全固定された場合の固有振動数問題は弱形式で以下のように表される。
【００８０】
式（２６）
【数２６】
【００８１】
ここで、ωｋはｋ次の固有振動数を表し、φｋはそれに対応する固有モードベクトルを表す。本発明では、固有振動数最適化問題に対し、以下の、固有振動数の平方の重み付き和を目標関数とする。
【００８２】
式（２７）
【数２７】
【００８３】
ここで、固有振動数ｗｋは以下の式により表される。
式（２８）
【数２８】
【００８４】
以上の目標関数に関する最小化問題に、式（１５）に含まれる体積制約を、ラグランジュ未定乗数μを乗じて加え、最終的に最適化問題は以下のように定式化される。
【００８５】
式（２９）
【数２９】
【００８６】
（フェーズフィールド関数の発展方程式）
式（１８）におけるフェーズフィールド関数は仮想的な時間ｔ＞０を設定し、以下の方程式に従い、更新されるものとする。
【００８７】
式（３０）
【数３０】
【００８８】
この式は上述したように、ｆ（φ）が二重井戸型ポテンシャルであり、φの総体積が一定となる場合には、以下のように、表面エネルギーとポテンシャルｆ（φ）の領域Ｄでの総和と領域ΩＡのペリメータの最小化問題と解釈できる。
【００８９】
式（３１）
【数３１】
【００９０】
ここで、ｆ（φ）を以下のように定義する。
【００９１】
式（３２）
【数３２】
【００９２】
ここで、ｗ（φ）はｗ（０）＝ｗ（１）＝ｗ’（０）＝ｗ’（１）を満たす関数、即ち等しい極小値を持つ二重井戸型ポテンシャルである。Ｗ（ｘ）はポテンシャルの高さを決定する係数である。また、ｇ（φ）はｇ（０）＝０、ｇ（１）＝１、ｇ’（０）＝ｇ’（１）＝０を満たす関数であり、Ｈ（ｘ）はその係数である。この項により、ポテンシャルｆ（φ）の極小値に差が与えられる。上述したように、式（３２）においては、領域ξにおいて、関数φはポテンシャルｆ（φ）の極小値が低い方に変化する。ここで、図６のように、Ｈ（ｘ）を以下のようにする。
【００９３】
式（３３）
【数３３】
【００９４】
式（３３）のＪ−’’の項を無視すると、位置ｘにおけるφの変化は以下のように表される。
【００９５】
式（３４）
【数３４】
【００９６】
ここでη（ｔ）＞０はＨ（ｘ）の値に基づき決定される、φの時間ｔでの変化量を表す値である。これは、領域ξにおいて、φが通常の勾配法と同様に、目標関数のφに関する感度Ｊ−’（ｘ，φ）に従い更新されることを示している。即ち、この定式化により、目標関数の減少が通常の勾配法と同様に議論可能なことを示す。
【００９７】
（感度解析）
本発明では、前節で述べたように、式（３０）のポテンシャルの決定に目標関数ｆ（φ）の関数φに関する感度を用いる。ここでは、目標関数Ｊ’の関数φに関する感度Ｊ−’（ｘ，φ）を導出する。最初に、式（２３）の平均コンプライアンスの感度を導出する。随伴状態ｐを用いて、ラグランジアンＬ（φ）を定義する。
【００９８】
式（３５）
【数３５】
【００９９】
このＬ（φ）の微小変化量を以下のように導出する。
【０１００】
式（３６）
【数３６】
【０１０１】
式（３６）においける右辺の第二項、第三項は等式制約となる式（２０）でｖ＝∂ｐ／∂φとした際に等しいので、値が0となる。さらに、第四項、第五項は同式において、ｕ＝ｐ、ｖ＝∂ｕ／∂φとした場合に等しいため、同様に値が０になる。よってｕとｐの支配方程式が等しく、ｐ＝ｕとなることから、この問題は自己随伴問題であり、目標関数Ｊ（φ）のφに関する感度は以下のように表される。
【０１０２】
式（３７）
【数３７】
【０１０３】
続いて、メカニズム創生問題で使用される、式（２４）の目標変位と変位との最小自乗誤差の感度を導出する。最初に、ラグランジアンＬ（φ）を以下のように定義する。
【０１０４】
式（３８）
【数３８】
【０１０５】
このＬ（φ）の微小変化量を以下のように導出する。
【０１０６】
式（３９）
【数３９】
【０１０７】
ただし、Ｃ０は以下で表される定数である。
【０１０８】
式（４０）
【数４０】
【０１０９】
右辺第二項、第三項は式（３６）と同様に値が０となる。第四項、第五項は随伴状態ｐが以下の式を満たすものとすると、値が０になる。
【０１１０】
式（４１）
【数４１】
【０１１１】
よって、目標関数Ｊ（φ）のφに関する感度は以下のように表される。
【０１１２】
式（４２）
【数４２】
【０１１３】
ただし、随伴状態ｐは式（４２）を満たす。
【０１１４】
続いて、式（２７）に含まれるｋ次の固有値の感度を導出する。随伴状態ｐを用いて、ラグランジアンＬ（φ）を以下のように定義する。
【０１１５】
式（４３）
【数４３】
【０１１６】
このＬ（φ）の微小変化量を以下のように導出する。
【０１１７】
式（４４）
【数４４】
【０１１８】
式（４４）の右辺の第五項、第六項は等式制約となる式（２６）でｖ＝∂ｐ／∂φとした際に等しいので、値が０となる。さらに、第七項、第八項は同式において、ｕ＝ｐ、ｖ＝∂φｋ／∂φとした場合に等しいため、同様に値が０になる。よってｕとｐの支配方程式が等しく、ｐ＝ｕとなることから、この問題は自己随伴問題である。さらに、φｋが正規化されており、∫Ｄρ｜φｋ｜２ｄｘ＝０とすると、第三項、第四項は相殺される。よって、ｋ次の固有値λｋのφに関する感度は以下のように表される。
【０１１９】
式（４５）
【数４５】
【０１２０】
（アルゴリズム）
以上をふまえ、形状最適化装置１０は、図７に示すような形状最適化処理を実行する。まず、ＣＰＵ１１は、初期形状Ω０に対応するφの初期値を決定する（ステップＳ２０）。次に、有限要素法により変位ｕ、固有モードφｋ、随伴状態ｐを導出する（ステップＳ２１）。続いて、目標関数を導出する（ステップＳ２２）。そして、式（４２）を用いて目標関数の感度を導出する（ステップＳ２３）。次に、導出した目標関数及びその感度に基づいて二重井戸型ポテンシャルｆ（φ）を式（３２）から導出して、式（３０）の反応拡散方程式（フェーズフィールド方程式）を有限体積法により解く（ステップＳ２４）。そして、以上の処理を目標関数が収束するまで繰り返す。なお、ステップＳ２１で導出する物理量は、線形弾性場あるいは材料非線形性及び幾何学的非線形性を有する連続体の弾性場における変位あるいは固有振動数、熱伝導場における定常温度、ストークス流れ場の定常流速及び圧力、ポテンシャル流れ場定常流速及び圧力、音場における定常音圧、電磁場問題における電界のいずれであってもよい。
【０１２１】
（実施例）
本実施例では、式（３０）に用いる関数を以下のように設定する。
【０１２２】
式（４６）
【数４６】
【０１２３】
さらに、関数ｗ（φ）の高さを決定する係数Ｗ（ｘ）はＷ（ｘ）＝（０．１／０．５４）×Ｈ（ｘ）とする。
【０１２４】
さらに、式（１９）で表される、領域ξにおける弾性定数には、トポロジー最適化法におけるＳＩＭＰ法で使用されている、ペナルティ係数ｐを定義し、関数φをｐ乗した値を材料の元の弾性テンソルに乗じる方法により以下のように導出する。
【０１２５】
式（４７）
【数４７】
【０１２６】
また、式（２５）で表される、領域ξにおけるρ−の値は、ρを係数とするφの線形関数として以下のように表す。
【０１２７】
式（４８）
【数４８】
ここで、有限体積法に関して若干補足する。離散化された各要素の体積ＣＶについて以下の保存則が成り立つとする。
【０１２８】
式（４９）
【数４９】
本実施例では以下のように、拡散項は陰解法による離散化を行い、反応項は湧き出し項として扱う。
【０１２９】
式（５０）
【数５０】
【０１３０】
ここで、φｎＰは時間ｎにおける要素Ｐのφの値、ＶＰは要素Ｐの体積、Σｉは要素Ｐの境界に関する総和、Ａｉは要素Ｐの境界の長さもしくは面積、ｄＡＰは隣り合う要素同士の中心間距離を示す。この式をｔについて解くことにより、各ｔにおけるφの値が導出される。なお、時間ｔの更新幅Δｔは、勾配法における最急降下幅に相当するため、以下のように更新されるものとする。ｋ回目の繰り返し計算における関数φの値をφｋ、目標関数の値をＪ−ｋ（φｋ）、時間ｔの更新幅をΔｔｋとすると、
１．｛Ｊ−ｋ（φｋ）−Ｊ−ｋ−１（φｋー１）｝／Ｊ−ｋ−１（φｋー１）≦Ｊ−ｔｈｒｅのとき、Δｔｋ−１＝ｍａｘ（１．５Δｔｋ，Δｔｍａｘ）として、ｋ＋１回目の繰り返し計算を行う。
２．｛Ｊ−ｋ（φｋ）−Ｊ−ｋ−１（φｋー１）｝／Ｊ−ｋ−１（φｋー１）＞Ｊ−ｔｈｒｅのとき、Δｔｋ＝０．５Δｔｋとして、φｋを再導出し、ｋ回目の繰り返し計算をやり直す。そして、上記１．の条件を満たすまで、これを繰り返す。
【０１３１】
簡単な数値例を用いて本手法の妥当性を検証する。説明の簡単のため、各物理量は無次元化された単位量として扱う。いずれの数値例においても、材料のヤング率Ｅは１．０、ポアソン比νは０．３、質量密度ρは１．０とする。また、式（４７）のパラメータｐの値は3とする。さらに、式（３０）の計算においては、有限要素法によるｆ（φ）の計算と、有限体積法によるφの更新を１ステップずつ交互に行うものとする。なお、関数φは各有限要素ごとに設定され、有限要素法と有限体積法においては共通のメッシュを使用するものとする。
【０１３２】
（剛性最大化問題）
最初に、式（２３）の平均コンプライアンスの最小化を目的とした、二次元片持ち梁の最適化を行う。設計領域（２００×４００）を図８に示す。設計領域の左端を完全固定し、右端中央に垂直荷重１．０を作用させる。
【０１３３】
（初期案が最適形状に与える影響）
初期案が最適形状に与える影響について検討を行う。図８の設計領域を１００×２００個の正方形４節点アイソパラメトリック要素で分割する。図９にそれぞれ、穴の個数が（ａ）９個と（ｂ）１７個と（ｃ）３９個の場合の初期解を示す。穴の半径はそれぞれ（ａ）４６と（ｂ）３０と（ｃ）２２である。式（３０）における係数κは０．５、式（３３）における係数ｈは２．０とし、時間ｔの最大更新幅Δｔｋは１．０とした。体積制約は総体積の約４０％と設定し、式（２９）のラグランジュ未定乗数を随時更新した。
【０１３４】
図１０にそれぞれ、初期案の穴の個数が（ａ）（ｂ）9個と（ｃ）（ｄ）17個と（ｅ）（ｆ）39個の場合の最適解及び最適化途中の形状を示す。いずれの初期解を使用した場合においても、最適解に初期解依存性が確認できる。これは、本手法では穴の生成が不可能なためである。しかし、初期案、最適化途中、最適解の図を比較してもわかるように、穴の消失に限り形態変化が可能である。図１１に平均コンプライアンスの収束履歴を、図１２に目標関数の収束履歴を示す。いずれにおいても、滑らかな収束が得られていることがわかる。
【０１３５】
なお、この発明は上記実施形態に限定されず、種々の変形及び応用が可能である。例えば、上記実施形態では、ＣＰＵ１１が実行するプログラムは、予め記憶装置１３に格納されているものとして説明したが、外部記憶媒体や通信網を介して形状最適化装置１０の所定領域に格納されるものであってもよい。
【産業上の利用可能性】
【０１３６】
以上説明したように、本発明によれば、形状最適化においてフェーズフィールド法を適用することにより、計算の処理負担の大きい有限要素法の計算を減らしたので、形状最適化の計算の処理負担を低減することができる。
【図面の簡単な説明】
【０１３７】
【図１】形状最適化装置の構成を示すブロック図である。
【図２】フェーズフィールド法を説明するための図である。
【図３】井戸型ポテンシャルを説明するための図である。
【図４】異なる相間でのインターフェイスを説明するための図である。
【図５】領域の一例を示す図である。
【図６】井戸型ポテンシャルを説明するための図である。
【図７】フェーズフィールド法に基づく形状最適化処理を説明するためのフローチャートである。
【図８】実施例における設計領域を示す図である。
【図９】実施例の設計領域の初期解を示す図である。
【図１０】実施例の設計領域の最適解及び最適化途中の形状を示す。
【図１１】実施例の平均コンプライアンスの収束履歴を示す図である。
【図１２】実施例の目標関数の収束履歴を示す図である。
【図１３】レベルセット法に基づく形状最適化処理を説明するためのフローチャートである。
【符号の説明】
【０１３８】
１０ 形状最適化装置
１１ ＣＰＵ
１２ 入力装置
１３ 記憶装置
１４ 表示制御部
２０ ディスプレイ
【技術分野】
【０００１】
本発明は、形状最適化方法、形状最適化装置、及び、プログラムに関する。
【背景技術】
【０００２】
偏微分方程式で特定の物理問題が記述される領域において、特定の物理量の最適化を目的とし、その領域形状を変化させる構造最適化問題を考える。このような問題は、機械製品の性能向上を目指す上で非常に重要な技術であり、様々な解法が提案されている。これらの解法の一般的な区分として、対象となる領域の境界を移動することにより最適化を行う形状最適化がある。また、対象の形状のみならず形態も最適化可能なトポロジー最適化が挙げられる。
【０００３】
しかし、一般的な構造最適化の方法においては二つの問題点がある。一つは領域を支配する偏微分方程式の数値計算を行う際に、有限要素法を用いるのが一般的であるが、領域の形状が更新される度に有限要素メッシュを再生成する必要があり、計算の処理負担が非常に大きくなることである。
【０００４】
もう一つは、スプライン関数等で領域を陽的に表現するため、境界同士の接触による穴の消失等の境界のトポロジー変化を表現できないことである。
【０００５】
トポロジー最適化は、多くの数学者によって基礎理論が構築されている。この方法は領域の形状のみならず形態も最適化可能であり、様々な工学的問題に適用され、成功を収めている。しかし、基本的には最適構造が材料密度分布で表されるため、そこから明確な構造を得ることが難しく、圧力や熱放射等の境界形状に依存するような境界条件の設定が困難という問題点がある。
【０００６】
近年、これらの問題を解消する方法として、レベルセット法を用いた形状最適化が提案されている（例えば非特許文献１）。レベルセット法は平均曲率流の数値計算法として提案され、現在では数値流体力学、画像処理等の様々な分野における移動境界問題に適用されている。レベルセット法では、領域の境界が陰関数を用いて表現されるため、オイラー型のメッシュを用いて数値計算が可能であり最適化の過程で計算量が一定である。また、境界同士の接触、穴の消失といった形態の変化が考慮可能である等の利点がある。
【０００７】
ここで、レベルセット法に基づく形状最適化における計算処理について簡単に説明する。図１３は、レベルセット法に基づく形状最適化処理を説明するためのフローチャートである。図１３に示すように、先ず、初期形状、各種初期値を決定する（ステップＳ１０）。そして、決定された初期形状について、有限要素法を用いて場の問題を解く（ステップＳ１１）。これにより、初期形状の変位や応力などを導出する。そして、その結果に基づいて目標関数を導出する（ステップＳ１２）。さらに、感度νを導出する（ステップＳ１３）。導出されたνを有限要素法を用いて滑らかにする（ステップＳ１４）。次に、レベルセット関数の初期化を行う（ステップＳ１５）。続いて、ハミルトン・ヤコビ方程式を解く（ステップＳ１６）。そして、計算値が収束したか否かを判定し（ステップＳ１７）、このようなステップＳ１１〜Ｓ１６を計算値が収束するまで繰り返す。
【０００８】
【非特許文献１】G．Allaire, F．Jouve， A．Toader， Structural optimization using sensitivity analysis and a level-set method、 Journal of Computational Physics 194 (2004) pp．363-393
【発明の開示】
【発明が解決しようとする課題】
【０００９】
しかし、レベルセット法における発展方程式であるハミルトン・ヤコビ方程式は移流方程式であり、数値計算上不安定になり易い。ハミルトン・ヤコビ方程式において安定した解を得るためには、最適化の繰り返し計算が必要となり、さらに、必ずしも安定した解が得られるわけではなかった。そのため、形状最適化にレベルセット法を適用しても依然として計算の処理負担は大きかった。
【００１０】
本発明は、上記実情に鑑みてなされたものであって、形状最適化における計算の処理負担を低減させることを目的とする。
【課題を解決するための手段】
【００１１】
上記目的を達成するため、本発明の第１の観点に係る形状最適化方法は、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化方法であって、
コンピュータが、
有限要素法により前記形状の物理量を導出し、
前記導出した物理量に基づいて目標関数を導出し、
前記導出した目標関数を最小化するために用いる当該目標関数の感度を導出し、
前記導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義し、
前記反応拡散方程式を有限体積法で解くことにより、前記導出した目標関数を最小化する
最小化処理を実行することを特徴とする。
【００１２】
また、前記形状が収束するまで、前記最小化処理を繰り返してもよい。また、前記物理量が、線形弾性場あるいは材料非線形性及び幾何学的非線形性を有する連続体の弾性場における変位あるいは固有振動数、熱伝導場における定常温度、ストークス流れ場の定常流速及び圧力、ポテンシャル流れ場定常流速及び圧力、音場における定常音圧、電磁場問題における電界のいずれかであってもよい。
【００１３】
本発明の第２の観点に係る形状最適化装置は、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置であって、
有限要素法により前記形状の物理量を導出する物理量導出手段と、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段と、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段と、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段と、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段と、
を備えることを特徴とする。
【００１４】
本発明の第３の観点に係るプログラムは、コンピュータを、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置として機能させるためのプログラムであって、
コンピュータを、
有限要素法により前記形状の物理量を導出する物理量導出手段、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段、
として機能させることを特徴とする。
【発明の効果】
【００１５】
本発明によれば、レベルセット法を用いた形状最適化計算に比べて、有限要素法の計算回数を減らすことができるので、形状最適化計算の処理負担を低減できる。
【発明を実施するための最良の形態】
【００１６】
本発明の実施の形態にかかる形状最適化方法を実行する形状最適化装置１０について図面を参照しながら説明する。本発明の形状最適化方法はフェーズフィールド法を用いた形状最適化方法である。より具体的には、フェーズフィールド法を用い、オイラー型メッシュでの計算が可能で、形態変化にも対応可能であり、かつ数値計算上安定した形状最適化方法である。
【００１７】
図１は、形状最適化装置１０の構成を示すブロック図である。形状最適化装置１０は、ＣＰＵ（Central Processing Unit）１１、入力装置１２、記憶装置１３、表示制御部１４から構成される。形状最適化装置１０は、例えば、汎用または専用のコンピュータである。
【００１８】
ＣＰＵ１１は、形状最適化装置１０全体の動作を制御する。また、ＣＰＵ１１は、記憶装置１３の所定領域に格納されるプログラムを読み込むことにより、各種の処理を実行する。特に本発明においては、ＣＰＵ１１は、フェーズフィールド法を用いた形状最適化処理を実行する。
【００１９】
入力装置１２は、キーボードやマウスなどから構成され、ユーザの操作に応じた操作信号をＣＰＵ１１に入力する。
【００２０】
記憶装置１３は、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、ハードディスク装置、その他の記憶媒体などから構成される。記憶装置１３は、フェーズフィールド法を用いた形状最適化処理のプログラムや必要なデータを記憶する。また、記憶装置１３は、ＣＰＵ１１のワークエリアとして機能する。
【００２１】
表示制御部１４は、ＣＰＵ１１から伝送された表示信号をディスプレイ２０で表示可能な形式に復号し、復号した信号をディスプレイ２０に入力する。
【００２２】
ディスプレイ２０は、ＬＣＤ（Liquid Crystal Display）などから構成され、形状最適化装置１０から伝送された信号に基づいて各種画面を表示する。例えば、本発明の形状最適化装置１０により最適化される形状のイメージなどを表示する。
【００２３】
以下フェーズフィールド法、及び、形状最適化装置１０が実行するフェーズフィールド法を用いた形状最適化処理について説明する。
【００２４】
フェーズフィールド法は、相転移における移動境界問題をシミュレートするための方法として開発された。フェーズフィールド法においては、領域全体に相の状態を表すフェーズフィールド関数を設定し、異なる相にそれぞれ対応する値を設定する。これにより、領域内のある箇所における相の状態が関数の値により決定される。異なる相間の境界では、フェーズフィールド関数は滑らかに補間され、その領域はインターフェイスと呼ばれる。領域内に二つの相が存在するとき、フェーズフィールド関数をφ（ｘ）とし、二つの相にそれぞれφ＝０、φ＝１とした場合の例を図２に示す。物理的な意味としては、フェーズフィールド関数は領域内の各箇所における相の平均化された値を表す。このフェーズフィールド関数を用いて、二つの相に関するVan der Waalsの自由エネルギーは以下の式（１）ように表される。
【００２５】
式（１）
【数１】
【００２６】
式（１）の第一項は平均場理論における干渉エネルギーを表し、第二項は二重井戸型ポテンシャルを表す。図２のように二つの相にそれぞれφ＝０、φ＝１が対応する場合は、井戸型ポテンシャルはｆ’（０）＝ｆ’（１）＝０という性質を持つ。この二重井戸型ポテンシャルは、それぞれの相において、自由エネルギーが極小値をとることを意味する。しかし、異なる相間で、フェーズフィールド関数が急激に変化し勾配が大きくなることは、第一項の値の増加に繋がるため、その間は０＜φ＜１の値をとる微小なインターフェイス領域が存在し、関数は滑らかに補間されることになる。
【００２７】
次に、このフェーズフィールド関数の時間変化を決定する発展方程式について述べる。最初に、フェーズフィールド関数が領域内でその総和が保存される量である場合を考える。この保存系のフェーズフィールド関数を改めてＣと表すと、領域内で以下の連続の式が成り立つ。
【００２８】
式（２）
【数２】
【００２９】
ここで、ｊは流束ベクトルである。さらに、閉じた系での自由エネルギーが増加しないと仮定すると、以下の式が成り立つ。
【００３０】
式（３）
【数３】
【００３１】
式（２）を式（３）に代入して、流入のない境界条件を適用すると、次式が得られる。
【００３２】
式（４）
【数４】
【００３３】
これが非正であるためには、ｊがＦに関して線形であると仮定すると、式（５）が得られる。
【００３４】
式（５）
【数５】
【００３５】
ただし、Ｍ（Ｃ）＞０とする。これを式（４）に代入すると、式（６）が得られる。
【００３６】
式（６）
【数６】
【００３７】
説明の簡単のためＭを定数とし、Ｆに式（１）を用いると、式（７）が得られる。
【００３８】
式（７）
【数７】
【００３９】
式（７）は、Ｃａｈｎ−Ｈｉｌｌｉａｒｄ方程式と呼ばれ、保存場に関するフェーズフィールド関数の発展方程式として用いられる。
【００４０】
続いて、フェーズフィールド関数が領域内でその総和が保存されない量である場合について考える。フェーズフィールド関数をφとすると、式（３）は非保存の場合についても成り立つので、以下の式が成り立つ。
【００４１】
式（８）
【数８】
【００４２】
これが任意のφで満たされるためには、下記の式が満たされなければならない。
【００４３】
式（９）
【数９】
【００４４】
式（９）はＡｌｌｅｎ−Ｃａｈｎ方程式と呼ばれ、非保存場のフェーズフィールド関数の発展方程式として用いられる。
【００４５】
続いて、フェーズフィールド法とペリメータ最小化問題について述べる。ここでは、Ｃａｈｎ−Ｈｉｌｌｉａｒｄ方程式で記述される系のように、領域内で秩序変数が保存される場合の解の存在について述べる。先の例で挙げたように、フェーズフィールド関数φ（Ω；［α，β］）のαとβとがそれぞれ異なる相に対応する場合に、以下の最小化問題を考える。
【００４６】
式（１０）
【数１０】
【００４７】
ここで、Ｗ（φ）は図３（ａ）のような井戸型ポテンシャルに、アフィン変換ｃ１φ＋ｃ２を加えた、極小値が等しい二重井戸型ポテンシャルである。この変換は保存系においては定数になるため、式（１０）の最小化問題に影響を与えない。
【００４８】
式（１１）
【数１１】
ただし、εは０に近い微小な値とする。
【００４９】
この問題の解は以下のように表される。
【００５０】
式（１２）
【数１２】
【００５１】
ここで、φ（Ω；｛α，β｝）は領域Ωで最小の境界Ｓを有する関数であり、φ
１はφ１（∞）＝０となる関数で、図４に示すような異なる相間でのインターフェイスを与える。式（１０）において、第一項はφをαとβに近づける効果があり、第二項は不要なインターフェイス領域に制約を与える効果がある。式（１２）の解を式（１０）に代入して考えると、両項がε−１のオーダーになるため、ε→０のとき、式（１０）は以下に示すインターフェイス最小化問題と解釈可能である。
【００５２】
式（１３）
【数１３】
【００５３】
ここで、Ｐｅｒ（Ａ，Ω）はＡのΩにおけるペリメータを表す。これにより、φは［α，β］の値をとる滑らかな関数で定義されているのに対し、ε→０のとき式（１０）の解は結果として｛α，β｝の不連続な値となる。
【００５４】
次に、フェーズフィールド法における界面移動について説明する。式（９）のＡｌｌｅｎ−Ｃａｈｎ方程式に従い、発展する系においては、二重井戸型ポテンシャルの効果により、十分時間がたてばインターフェイス領域が形成され、それが法線方向に次式で現される速度で移動することが知られている。
【００５５】
式（１４）
【数１４】
【００５６】
式（１４）において、αは二重井戸型ポテンシャルによって定まる定数であり、κはインターフェイスの曲率を表す。二重井戸型ポテンシャルの二つの極小値が等しい場合には、αは０となり、式（１４）は完全な平均曲率運動を表す。極小値間に差がある場合には、インターフェイス領域におけるポテンシャルの値が小さくなるような方向に速度αが定まる。
【００５７】
本発明では、式（１４）のインターフェイスの放線方向に関する性質を用いて、フェーズフィールド法を用いた形状最適化方法の定式化を提供する。
【００５８】
（フェーズフィールド関数を用いた領域表現）
ここでは、Ａｌｌｅｎ−Ｃａｈｎ方程式に支配される系において、インターフェイスの移動が式（１４）で表されることを利用し、フェーズフィールド法を用いたオイラー型形状最適化法を提案する。最初に、ｄ次元実数空間に属する、体積が一定値Ｖとなる領域Ωの最適化問題を考える。このとき、図５（ａ）に示すように、Ωを包括するような設計領域Ｄを考え、その中の許容される解の集合Ｕａｄを以下のように定義する。
【００５９】
式（１５）
【数１５】
【００６０】
このとき、領域Ωの形状最適化問題は以下のように定義される。
【００６１】
式（１６）
【数１６】
【００６２】
ここで、Ｊ（Ω）は領域Ωで定義される任意の目標関数である。ここでは、この問題の数値解法として、設計領域Ｄ全体にフェーズフィールド関数φを定義し、領域Ωを陰的に表すことを考える。即ち、図５（ｂ）に示すように、ΩＡ⊂ΩとなるΩＡと、ΩＢ⊂Ｄ＼ΩとなるΩＢを定義し、さらに、領域ＤのうちΩＡとΩＢの間に、微小な厚みを持つ領域ξが存在すると仮定し、ξ⊂Ｄ＼（ΩＡ∪ΩＢ）を定義する。即ち、以下の関係が成り立つと仮定する。
【００６３】
式（１７）
【数１７】
【００６４】
このとき、これらの領域ΩＡとΩＢとはフェーズフィールド関数を用いて以下のように表される。
【００６５】
式（１８）
【数１８】
【００６６】
以上の定義において、境界∂Ωはξ内に存在することがわかるが、その位置は不定である。即ち、式（１８）によって完全に領域Ωが厳密に表されるわけではない。しかし、領域ξは微小な厚みを持つ領域であり、その厚みが十分小さい場合は、と見なすことができる。このとき、∂Ω≒ξとなり、領域Ωが式（１８）により、φで近似的に表されることになる。
【００６７】
式（１８）のように、［０，１］範囲で定義される関数を用いた領域表現方法は、均質化法に基づくトポロジー最適化でも用いられているが、その物理的意味は全く異なる。即ち、均質化法に基づくトポロジー最適化においては、材料密度分布関数で表された領域の最適化問題は、特性関数で表される領域の最適化問題の弱収束解と厳密にその物理的意味が解釈可能である。しかし、通常、式（１）で表されるように関数の勾配を目標関数に含むような場合には均質化法を適用することができない。即ち、関数φの物理的な意味は不明瞭であるが、ここでは領域表現のための陰関数としてのみ定義する。
【００６８】
（目標関数について）
領域Ｄにおける、線形弾性体問題を考える。領域ΩＡが弾性テンソルＡを有する材料で構成され、領域ΩＢを非常に弱い材料が占めていると仮定し、領域Ωにおける弾性テンソルＡ＊を以下のように定義する。
【００６９】
式（１９）
【数１９】
【００７０】
ここで、領域ξにおいては、材料の存在が明確に定義できないため、弾性テンソルの厳密な導出が不可能である。さらに、領域Ｄの境界∂Ｄのうち、ディリクレ境界条件が設定される境界をΓＤ、ノイマン境界条件が設定される境界をΓＮとする。今、領域Ｄが境界∂Ｄで完全固定され、領域Ｄに作用する体積力が０であるとし、境界ΓＮに分布荷重ｇが作用するものとすると、領域Ｄにおける線形弾性問題は弱形式で以下のように表される。
【００７１】
式（２０）
【数２０】
式（２１）
【数２１】
式（２２）
【数２２】
【００７２】
ここで、ｕは変位ベクトル、ｖはテスト関数、ｅはひずみテンソルを表す。
【００７３】
本発明では、線形弾性問題として、剛性最大化問題と、特定箇所の変位を目標値に近づけるメカニズム創生問題を考える。コンプライアンス最小化問題においては、以下で表される平均コンプライアンスを目標関数とする。
【００７４】
式（２３）
【数２３】
【００７５】
コンプライアントメカニズム創生問題では、以下で表される、変位と変位の目標値との最小自乗誤差を目標関数とする。
【００７６】
式（２４）
【数２４】
【００７７】
さらに、領域ΩＡを構成する材料の質量密度をρとし、領域ΩＢを占める材料の質量密度を非常に小さな値に設定すると、領域Ｄにおける質量密度ρ−は以下のように表される。
【００７８】
式（２５）
【数２５】
【００７９】
このとき、領域Ｄが境界∂Ｄで完全固定された場合の固有振動数問題は弱形式で以下のように表される。
【００８０】
式（２６）
【数２６】
【００８１】
ここで、ωｋはｋ次の固有振動数を表し、φｋはそれに対応する固有モードベクトルを表す。本発明では、固有振動数最適化問題に対し、以下の、固有振動数の平方の重み付き和を目標関数とする。
【００８２】
式（２７）
【数２７】
【００８３】
ここで、固有振動数ｗｋは以下の式により表される。
式（２８）
【数２８】
【００８４】
以上の目標関数に関する最小化問題に、式（１５）に含まれる体積制約を、ラグランジュ未定乗数μを乗じて加え、最終的に最適化問題は以下のように定式化される。
【００８５】
式（２９）
【数２９】
【００８６】
（フェーズフィールド関数の発展方程式）
式（１８）におけるフェーズフィールド関数は仮想的な時間ｔ＞０を設定し、以下の方程式に従い、更新されるものとする。
【００８７】
式（３０）
【数３０】
【００８８】
この式は上述したように、ｆ（φ）が二重井戸型ポテンシャルであり、φの総体積が一定となる場合には、以下のように、表面エネルギーとポテンシャルｆ（φ）の領域Ｄでの総和と領域ΩＡのペリメータの最小化問題と解釈できる。
【００８９】
式（３１）
【数３１】
【００９０】
ここで、ｆ（φ）を以下のように定義する。
【００９１】
式（３２）
【数３２】
【００９２】
ここで、ｗ（φ）はｗ（０）＝ｗ（１）＝ｗ’（０）＝ｗ’（１）を満たす関数、即ち等しい極小値を持つ二重井戸型ポテンシャルである。Ｗ（ｘ）はポテンシャルの高さを決定する係数である。また、ｇ（φ）はｇ（０）＝０、ｇ（１）＝１、ｇ’（０）＝ｇ’（１）＝０を満たす関数であり、Ｈ（ｘ）はその係数である。この項により、ポテンシャルｆ（φ）の極小値に差が与えられる。上述したように、式（３２）においては、領域ξにおいて、関数φはポテンシャルｆ（φ）の極小値が低い方に変化する。ここで、図６のように、Ｈ（ｘ）を以下のようにする。
【００９３】
式（３３）
【数３３】
【００９４】
式（３３）のＪ−’’の項を無視すると、位置ｘにおけるφの変化は以下のように表される。
【００９５】
式（３４）
【数３４】
【００９６】
ここでη（ｔ）＞０はＨ（ｘ）の値に基づき決定される、φの時間ｔでの変化量を表す値である。これは、領域ξにおいて、φが通常の勾配法と同様に、目標関数のφに関する感度Ｊ−’（ｘ，φ）に従い更新されることを示している。即ち、この定式化により、目標関数の減少が通常の勾配法と同様に議論可能なことを示す。
【００９７】
（感度解析）
本発明では、前節で述べたように、式（３０）のポテンシャルの決定に目標関数ｆ（φ）の関数φに関する感度を用いる。ここでは、目標関数Ｊ’の関数φに関する感度Ｊ−’（ｘ，φ）を導出する。最初に、式（２３）の平均コンプライアンスの感度を導出する。随伴状態ｐを用いて、ラグランジアンＬ（φ）を定義する。
【００９８】
式（３５）
【数３５】
【００９９】
このＬ（φ）の微小変化量を以下のように導出する。
【０１００】
式（３６）
【数３６】
【０１０１】
式（３６）においける右辺の第二項、第三項は等式制約となる式（２０）でｖ＝∂ｐ／∂φとした際に等しいので、値が0となる。さらに、第四項、第五項は同式において、ｕ＝ｐ、ｖ＝∂ｕ／∂φとした場合に等しいため、同様に値が０になる。よってｕとｐの支配方程式が等しく、ｐ＝ｕとなることから、この問題は自己随伴問題であり、目標関数Ｊ（φ）のφに関する感度は以下のように表される。
【０１０２】
式（３７）
【数３７】
【０１０３】
続いて、メカニズム創生問題で使用される、式（２４）の目標変位と変位との最小自乗誤差の感度を導出する。最初に、ラグランジアンＬ（φ）を以下のように定義する。
【０１０４】
式（３８）
【数３８】
【０１０５】
このＬ（φ）の微小変化量を以下のように導出する。
【０１０６】
式（３９）
【数３９】
【０１０７】
ただし、Ｃ０は以下で表される定数である。
【０１０８】
式（４０）
【数４０】
【０１０９】
右辺第二項、第三項は式（３６）と同様に値が０となる。第四項、第五項は随伴状態ｐが以下の式を満たすものとすると、値が０になる。
【０１１０】
式（４１）
【数４１】
【０１１１】
よって、目標関数Ｊ（φ）のφに関する感度は以下のように表される。
【０１１２】
式（４２）
【数４２】
【０１１３】
ただし、随伴状態ｐは式（４２）を満たす。
【０１１４】
続いて、式（２７）に含まれるｋ次の固有値の感度を導出する。随伴状態ｐを用いて、ラグランジアンＬ（φ）を以下のように定義する。
【０１１５】
式（４３）
【数４３】
【０１１６】
このＬ（φ）の微小変化量を以下のように導出する。
【０１１７】
式（４４）
【数４４】
【０１１８】
式（４４）の右辺の第五項、第六項は等式制約となる式（２６）でｖ＝∂ｐ／∂φとした際に等しいので、値が０となる。さらに、第七項、第八項は同式において、ｕ＝ｐ、ｖ＝∂φｋ／∂φとした場合に等しいため、同様に値が０になる。よってｕとｐの支配方程式が等しく、ｐ＝ｕとなることから、この問題は自己随伴問題である。さらに、φｋが正規化されており、∫Ｄρ｜φｋ｜２ｄｘ＝０とすると、第三項、第四項は相殺される。よって、ｋ次の固有値λｋのφに関する感度は以下のように表される。
【０１１９】
式（４５）
【数４５】
【０１２０】
（アルゴリズム）
以上をふまえ、形状最適化装置１０は、図７に示すような形状最適化処理を実行する。まず、ＣＰＵ１１は、初期形状Ω０に対応するφの初期値を決定する（ステップＳ２０）。次に、有限要素法により変位ｕ、固有モードφｋ、随伴状態ｐを導出する（ステップＳ２１）。続いて、目標関数を導出する（ステップＳ２２）。そして、式（４２）を用いて目標関数の感度を導出する（ステップＳ２３）。次に、導出した目標関数及びその感度に基づいて二重井戸型ポテンシャルｆ（φ）を式（３２）から導出して、式（３０）の反応拡散方程式（フェーズフィールド方程式）を有限体積法により解く（ステップＳ２４）。そして、以上の処理を目標関数が収束するまで繰り返す。なお、ステップＳ２１で導出する物理量は、線形弾性場あるいは材料非線形性及び幾何学的非線形性を有する連続体の弾性場における変位あるいは固有振動数、熱伝導場における定常温度、ストークス流れ場の定常流速及び圧力、ポテンシャル流れ場定常流速及び圧力、音場における定常音圧、電磁場問題における電界のいずれであってもよい。
【０１２１】
（実施例）
本実施例では、式（３０）に用いる関数を以下のように設定する。
【０１２２】
式（４６）
【数４６】
【０１２３】
さらに、関数ｗ（φ）の高さを決定する係数Ｗ（ｘ）はＷ（ｘ）＝（０．１／０．５４）×Ｈ（ｘ）とする。
【０１２４】
さらに、式（１９）で表される、領域ξにおける弾性定数には、トポロジー最適化法におけるＳＩＭＰ法で使用されている、ペナルティ係数ｐを定義し、関数φをｐ乗した値を材料の元の弾性テンソルに乗じる方法により以下のように導出する。
【０１２５】
式（４７）
【数４７】
【０１２６】
また、式（２５）で表される、領域ξにおけるρ−の値は、ρを係数とするφの線形関数として以下のように表す。
【０１２７】
式（４８）
【数４８】
ここで、有限体積法に関して若干補足する。離散化された各要素の体積ＣＶについて以下の保存則が成り立つとする。
【０１２８】
式（４９）
【数４９】
本実施例では以下のように、拡散項は陰解法による離散化を行い、反応項は湧き出し項として扱う。
【０１２９】
式（５０）
【数５０】
【０１３０】
ここで、φｎＰは時間ｎにおける要素Ｐのφの値、ＶＰは要素Ｐの体積、Σｉは要素Ｐの境界に関する総和、Ａｉは要素Ｐの境界の長さもしくは面積、ｄＡＰは隣り合う要素同士の中心間距離を示す。この式をｔについて解くことにより、各ｔにおけるφの値が導出される。なお、時間ｔの更新幅Δｔは、勾配法における最急降下幅に相当するため、以下のように更新されるものとする。ｋ回目の繰り返し計算における関数φの値をφｋ、目標関数の値をＪ−ｋ（φｋ）、時間ｔの更新幅をΔｔｋとすると、
１．｛Ｊ−ｋ（φｋ）−Ｊ−ｋ−１（φｋー１）｝／Ｊ−ｋ−１（φｋー１）≦Ｊ−ｔｈｒｅのとき、Δｔｋ−１＝ｍａｘ（１．５Δｔｋ，Δｔｍａｘ）として、ｋ＋１回目の繰り返し計算を行う。
２．｛Ｊ−ｋ（φｋ）−Ｊ−ｋ−１（φｋー１）｝／Ｊ−ｋ−１（φｋー１）＞Ｊ−ｔｈｒｅのとき、Δｔｋ＝０．５Δｔｋとして、φｋを再導出し、ｋ回目の繰り返し計算をやり直す。そして、上記１．の条件を満たすまで、これを繰り返す。
【０１３１】
簡単な数値例を用いて本手法の妥当性を検証する。説明の簡単のため、各物理量は無次元化された単位量として扱う。いずれの数値例においても、材料のヤング率Ｅは１．０、ポアソン比νは０．３、質量密度ρは１．０とする。また、式（４７）のパラメータｐの値は3とする。さらに、式（３０）の計算においては、有限要素法によるｆ（φ）の計算と、有限体積法によるφの更新を１ステップずつ交互に行うものとする。なお、関数φは各有限要素ごとに設定され、有限要素法と有限体積法においては共通のメッシュを使用するものとする。
【０１３２】
（剛性最大化問題）
最初に、式（２３）の平均コンプライアンスの最小化を目的とした、二次元片持ち梁の最適化を行う。設計領域（２００×４００）を図８に示す。設計領域の左端を完全固定し、右端中央に垂直荷重１．０を作用させる。
【０１３３】
（初期案が最適形状に与える影響）
初期案が最適形状に与える影響について検討を行う。図８の設計領域を１００×２００個の正方形４節点アイソパラメトリック要素で分割する。図９にそれぞれ、穴の個数が（ａ）９個と（ｂ）１７個と（ｃ）３９個の場合の初期解を示す。穴の半径はそれぞれ（ａ）４６と（ｂ）３０と（ｃ）２２である。式（３０）における係数κは０．５、式（３３）における係数ｈは２．０とし、時間ｔの最大更新幅Δｔｋは１．０とした。体積制約は総体積の約４０％と設定し、式（２９）のラグランジュ未定乗数を随時更新した。
【０１３４】
図１０にそれぞれ、初期案の穴の個数が（ａ）（ｂ）9個と（ｃ）（ｄ）17個と（ｅ）（ｆ）39個の場合の最適解及び最適化途中の形状を示す。いずれの初期解を使用した場合においても、最適解に初期解依存性が確認できる。これは、本手法では穴の生成が不可能なためである。しかし、初期案、最適化途中、最適解の図を比較してもわかるように、穴の消失に限り形態変化が可能である。図１１に平均コンプライアンスの収束履歴を、図１２に目標関数の収束履歴を示す。いずれにおいても、滑らかな収束が得られていることがわかる。
【０１３５】
なお、この発明は上記実施形態に限定されず、種々の変形及び応用が可能である。例えば、上記実施形態では、ＣＰＵ１１が実行するプログラムは、予め記憶装置１３に格納されているものとして説明したが、外部記憶媒体や通信網を介して形状最適化装置１０の所定領域に格納されるものであってもよい。
【産業上の利用可能性】
【０１３６】
以上説明したように、本発明によれば、形状最適化においてフェーズフィールド法を適用することにより、計算の処理負担の大きい有限要素法の計算を減らしたので、形状最適化の計算の処理負担を低減することができる。
【図面の簡単な説明】
【０１３７】
【図１】形状最適化装置の構成を示すブロック図である。
【図２】フェーズフィールド法を説明するための図である。
【図３】井戸型ポテンシャルを説明するための図である。
【図４】異なる相間でのインターフェイスを説明するための図である。
【図５】領域の一例を示す図である。
【図６】井戸型ポテンシャルを説明するための図である。
【図７】フェーズフィールド法に基づく形状最適化処理を説明するためのフローチャートである。
【図８】実施例における設計領域を示す図である。
【図９】実施例の設計領域の初期解を示す図である。
【図１０】実施例の設計領域の最適解及び最適化途中の形状を示す。
【図１１】実施例の平均コンプライアンスの収束履歴を示す図である。
【図１２】実施例の目標関数の収束履歴を示す図である。
【図１３】レベルセット法に基づく形状最適化処理を説明するためのフローチャートである。
【符号の説明】
【０１３８】
１０ 形状最適化装置
１１ ＣＰＵ
１２ 入力装置
１３ 記憶装置
１４ 表示制御部
２０ ディスプレイ
【特許請求の範囲】
【請求項１】
予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化方法であって、
コンピュータが、
有限要素法により前記形状の物理量を導出し、
前記導出した物理量に基づいて目標関数を導出し、
前記導出した目標関数を最小化するために用いる当該目標関数の感度を導出し、
前記導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義し、
前記反応拡散方程式を有限体積法で解くことにより、前記導出した目標関数を最小化する
最小化処理を実行することを特徴とする形状最適化方法。
【請求項２】
前記形状が収束するまで、前記最小化処理を繰り返す
ことを特徴とする請求項１に記載の形状最適化方法。
【請求項３】
前記物理量が、線形弾性場あるいは材料非線形性及び幾何学的非線形性を有する連続体の弾性場における変位あるいは固有振動数、熱伝導場における定常温度、ストークス流れ場の定常流速及び圧力、ポテンシャル流れ場定常流速及び圧力、音場における定常音圧、電磁場問題における電界のいずれかである
ことを特徴とする請求項１または２に記載の形状最適化方法。
【請求項４】
予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置であって、
有限要素法により前記形状の物理量を導出する物理量導出手段と、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段と、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段と、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段と、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段と、
を備えることを特徴とする形状最適化装置。
【請求項５】
コンピュータを、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置として機能させるためのプログラムであって、
コンピュータを、
有限要素法により前記形状の物理量を導出する物理量導出手段、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段、
として機能させることを特徴とするプログラム。
【請求項１】
予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化方法であって、
コンピュータが、
有限要素法により前記形状の物理量を導出し、
前記導出した物理量に基づいて目標関数を導出し、
前記導出した目標関数を最小化するために用いる当該目標関数の感度を導出し、
前記導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義し、
前記反応拡散方程式を有限体積法で解くことにより、前記導出した目標関数を最小化する
最小化処理を実行することを特徴とする形状最適化方法。
【請求項２】
前記形状が収束するまで、前記最小化処理を繰り返す
ことを特徴とする請求項１に記載の形状最適化方法。
【請求項３】
前記物理量が、線形弾性場あるいは材料非線形性及び幾何学的非線形性を有する連続体の弾性場における変位あるいは固有振動数、熱伝導場における定常温度、ストークス流れ場の定常流速及び圧力、ポテンシャル流れ場定常流速及び圧力、音場における定常音圧、電磁場問題における電界のいずれかである
ことを特徴とする請求項１または２に記載の形状最適化方法。
【請求項４】
予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置であって、
有限要素法により前記形状の物理量を導出する物理量導出手段と、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段と、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段と、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段と、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段と、
を備えることを特徴とする形状最適化装置。
【請求項５】
コンピュータを、予め与えられた物体または場の形状に、所定の負荷をかけた際の最適化形状を求める形状最適化装置として機能させるためのプログラムであって、
コンピュータを、
有限要素法により前記形状の物理量を導出する物理量導出手段、
前記物理量導出手段が導出した物理量に基づいて目標関数を導出する目標関数導出手段、
前記目標関数導出手段が導出した目標関数を最小化するために用いる当該目標関数の感度を導出する感度導出手段、
前記目標関数導出手段が導出した目標関数を最小化する方向に、形状の境界を移動させるために用いる二重井戸型ポテンシャルと、当該二重井戸型ポテンシャルを設定した反応拡散方程式と、を定義する定義手段、
前記定義手段が定義した前記反応拡散方程式を有限体積法で解くことにより、前記目標関数導出手段が導出した目標関数を最小化する手段、
として機能させることを特徴とするプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【公開番号】特開２０１０−１０８４５１（Ｐ２０１０−１０８４５１Ａ）
【公開日】平成２２年５月１３日（２０１０．５．１３）
【国際特許分類】
【出願番号】特願２００８−２８２５２５（Ｐ２００８−２８２５２５）
【出願日】平成２０年１０月３１日（２００８．１０．３１）
【出願人】（５０４１３６５６８）国立大学法人広島大学 (924)
【Ｆターム（参考）】
【公開日】平成２２年５月１３日（２０１０．５．１３）
【国際特許分類】
【出願日】平成２０年１０月３１日（２００８．１０．３１）
【出願人】（５０４１３６５６８）国立大学法人広島大学 (924)
【Ｆターム（参考）】
[ Back to top ]