JP2010277249A

JP2010277249A - 飛しょう体の形状決定方法及び形状決定装置

Info

Publication number: JP2010277249A
Application number: JP2009127847A
Authority: JP
Inventors: Shogo Harada; 昌悟原田
Original assignee: Toshiba Corp
Current assignee: Toshiba Corp
Priority date: 2009-05-27
Filing date: 2009-05-27
Publication date: 2010-12-09

Abstract

【課題】飛しょう体頭部の形状最適化に有効で効率的な形状決定方法および形状決定装置を提供する。
【解決手段】パラメトリック関数により飛しょう体の頭部の初期形状を設定するステップと、設定された前記初期形状の節点のおのおのについて衝撃波・膨張波法により、圧力を算出するステップと、形状表面の各節点の圧力を表面積により積分して抗力を算出し評価関数とするステップと、前記形状表面の節点を変えて各節点の圧力算出及び評価関数算出を繰り返して、所定の個体数の評価関数を得るステップと、この評価関数の対応する形状について、交叉、複写及び突然変異の処理を行って、次世代の形状を形成するステップと、形成された次世代について、圧力算出、評価関数の算出を行ってこの世代の所定の個体数の評価関数を取得し、これを所定の世代数について繰り返して最終世代の所定個体数の評価関数を得、適切な形状を決定するステップと、を有する。
【選択図】図１

Description

本発明は、流体中を移動する物体の形状を決定する方法に係り、特に飛しょう体の形状決定方法及び形状決定装置に関する。

流体中を物体が移動するとき流体抵抗を受けるので、その物体の形状を最適化して流体抵抗を減らすことは重要である。

例えば、特許文献１には、空気騒音が小さくなるように空気中を移動する物体形状を決定することが開示されている。その方法は、第１過程で初期形状を決定し第２過程で物体表面に数値流体力学用メッシュを作成し第３過程でナビエ・ストークス方程式を用いて物体表面圧力を求め第４過程で発生騒音を求め第５過程で第１過程から繰り返し、第６過程で最小騒音の形状を決定する。

また、特許文献２には、目的パラメータを突然変異や組換えによりパラメータセットの子を作るステップと、その子の質を評価するステップを有する進化戦略に基づく最適化方法が開示されている。更に特許文献３には、複数の設計変数を組み合わせた解候補である個体を遺伝子で表現して集団発生させ、個体の集団に対する遺伝的操作の繰り返しにより解を探索し遺伝的アルゴリズムを用いて流体の流れ場を構成する、ターボ流体機械の設計方法が開示されている。この特許文献３には、適応度の高い個体ほど親として選択される可能性が高くなるように、具体的には模擬２進交叉を用いて次世代候補を生成し、突然変異により次世代解候補の探索を行うことが記載されている。

特開平１１−３１１６３号公報特開２００２−３２４２２５号公報特開２００３−１０６２９８号公報

しかし、上記従来の形状決定方法は抗力を低減するための形状を決定するには過去の経験や既存の形状を元に模型を作成し風洞実験を繰り返し行って抗力を計測し要求の性能を満たす形状を採用する必要があった。特許文献１では、ＣＦＤ（Computaitional Fluid dynamics;計算流体力学）を用いることが記載されているが、この場合においても解析の対象となる形状について同様に試行錯誤を行う必要がある。

本発明はこれら従来の物体形状の決定方法に鑑みてなされたもので、飛しょう体頭部の形状最適化に有効で効率的な形状決定方法および形状決定装置を提供することを目的とする。

上記の課題を解決するために、本発明の請求項１は、パラメトリック関数により飛しょう体の頭部の初期形状を設定する初期形状決定ステップと、この初期形状設定ステップにより設定された前記初期形状の節点のおのおのについて衝撃波・膨張波法により、圧力を算出する圧力算出ステップと、形状表面の各節点の圧力を表面積により積分して抗力を算出し評価関数とする評価関数算出ステップと、前記形状表面の節点を変えて各節点の圧力算出及び評価関数算出を繰り返して、所定の個体数の評価関数を得る評価関数取得ステップと、この評価関数取得ステップにより得られた、評価関数の対応する形状について、交叉、複写及び突然変異の処理を行って、次世代の形状を形成する次世代形状形成ステップと、この次世代形状形成ステップにより形成された次世代について、前記圧力算出ステップによる圧力算出、前記評価関数算出ステップによる評価関数の算出を行ってこの世代の所定の個体数の評価関数を取得し、これを所定の世代数について繰り返して最終世代の所定個体数の評価関数を得、適切な形状を決定する形状決定ステップと、を有することを特徴とする、飛しょう体の形状決定方法を提供する。

本発明の特徴の１つは、形状を表示する曲線としてパラメトリック曲線を用いることであり、他の１つは、節点の圧力を計算するために衝撃波・膨張波法を用いることである。

本発明によれば、飛しょう体頭部の形状最適化に有効で効率的な形状決定方法および形状決定装置が得られる。

本発明の一実施形態の構成を示す図である。本発明の一実施形態の方法を示すフローチャートである。図２のフローチャートにおいて、衝撃波・膨張波法により各節点に圧力を算出する方法の一例を示すフローチャートである。本発明一実施形態の制御点と節点の関係の例を示す図である。本発明一実施形態において衝撃波・膨張波法により節点２における圧力を求める例を説明するための図である。本発明一実施形態において２つの速度の評価関数から最適な形状を求める点を説明するための図である。

以下、本発明の実施の形態について、詳細に説明する。図１に本発明一実施形態の構成を示す。

この形状決定装置は、最初に所定の形状を設定する初期形状設定部１０と、設定された形状に基づいて圧力算定プログラムの入力ファイルを作成する入力ファイル作成部１１と、この入力ファイルに基づいて各節点の圧力を算出する圧力算出部１２と、算出された圧力に基づいて圧力分布を算出する圧力分布算出部１３と、圧力分布から抗力を算出する抗力算出部１４と、この算出された抗力を評価関数とする評価関数決定部１５と、個体が所定数得られたか否かを判定する個体数判定部１６と、所定数の世代のデータが得られたか否かの判定を行う世代数判定部１７と、上記評価関数に基づいて次の形状データを作成する形状データ作成部１８と、前の世代の所定数の個体が得られたとき、次世代のデータを作成する次世代データ作成部１９と、最終世代の個体数のデータが得られたときそれらの評価関数により最適な形状を決定する形状決定部２０とを有する。

この形状決定装置による形状決定の方法を図２及び図３に示すフローチャートに基づいて説明する。この例では飛しょう体の頭部の最適な形状を決定するものとする。

まず、ステップＳ２０１で図１に示す初期形状設定部１０において、飛しょう体の頭部の初期の形状を仮に設定する。この形状として本願の発明では、ｏｇｉｖｅ形状やｐａｒａｂｏｌｉｃ形状などの既存の形状に限定されない、パラメトリック曲線を用いた形状を採用する。パラメトリック曲線の一例としてここではＢｅｚｉｅｒ曲線を用いる。Ｂｅｚｉｅｒ曲線は、一般的にはｎ個の制御点から成る曲線で、その内始点（ｉ＝０）と終点（i＝ｎ−１）以外の制御点を通らない曲線であり、次式により定義される。いまの場合この制御点の数ｎを例えば、４個とする。
但し、ｔはパラメータ値であり、Ｐｘ[ｉ]，Ｐｙ［ｉ］はｉ番目の制御点のｘｙ座標である。_ｎＣ_ｍ＝ｎ！／｛ｍ！（ｎ−ｍ）！｝である。

ｎ＝４としたときの、Bezier曲線の例を図４に実線で示す。横軸にｘ座標をとり、縦軸にｙ座標をとる。図４において、点４１は制御点を示している。

このように、形状を示す曲線は一部であるが、この曲線を中心軸Ｘの周りに回転させることにより飛しょう体頭部の立体形状を表すこととなる。Ｂｅｚｉｅｒ曲線は、始点と終点を基点として基点の制御点を通る必要がないので、基点を通る滑らかな曲線を描くことができる。それゆえ、最適化の過程で制御点の位置が大きく移動しても、滑らかな形状を維持することができる利点がある。

本発明に用いられる形状はＢｅｚｉｅｒ曲線に限られないが、パラメトリック曲線によって描かれた曲線上の任意の数の節点によって形状を表現する。

パラメトリック曲線例えばＢｅｚｉｅｒ曲線により初期形状を決定した後、ステップＳ２０２で世代数Ｇを１とし、ステップＳ２０３で個体数Ｐを１とする。また、この実施形態では、マッハ数が１の場合と２の場合について評価関数を求めるので、まずステップＳ２０３Ｍでは、マッハ数を１とする。

次のステップＳ２０４では、上記のように設定した形状のときの圧力算出プログラムの入力ファイルを作成する。

ステップＳ２０５では、圧力算出部１２において、衝撃波・膨張波法によって各節点の圧力を算出し、ステップＳ２０６では、圧力分布算出部１３において形状表面の各節点の圧力を表面積で積分し抗力算出部１４でその形状のときの抗力を算出し、ステップＳ２０６で評価関数決定部１５において、算出された抗力をその形状のときの評価関数とする。

この抗力の算出の過程について、ステップＳ２０６の図３に示す詳細フローチャートに基づいて更に詳しく説明する。まず、初期形状の制御点（ｎ）を設定して、上記式（１）にＰｘ［ｉ］，Ｐｙ[ｉ]に、ｉ＝０からｉ=（ｎ−１）を代入する。
形状を表現する曲線上のｍ＋１個の節点についてｔ＝０，１／ｍ，２／ｍ，・・・（ｍ−１）／ｍ，１と代入することにより節点の座標ｘ（ｔ），ｙ（ｔ）を求める。節点１から節点ｍ＋１までの圧力を順次求めていく。

例えば制御点の数ｎ＝４の場合、Ｐｘ［０］＝０、Ｐｙ[０]＝０、Ｐｘ[１]＝０、Ｐｙ[１]＝１、Ｐｘ[２]＝２、Ｐｙ[２]＝３、Ｐｘ[３]＝３、Ｐｙ[３]＝３、と任意に設定する。
このとき上記制御点によって決定する曲線上の点である節点が２１個であるとする。上記式（１）にｔ＝０、０．０５、０．１、・・・・０．９、０．９５、１．０を入れて各ｘ（ｔ）、ｙ（ｔ）を求める。

いま、第３節点４２を例に説明する。例えばｘ（ｔ）は、次に示す（２）式となる。またこのとき例えばｔ＝０．１のときのｘ（０．１）は次に示す（３）式となる。

ｙ（０．１）も同様にして求める。結局、第３節点（ｘ（０．１），ｙ（０．１）＝（０．０５７，０．３２７）となる。

このような演算をｔ＝０、０．０５、０．１、・・・・、０．９５、１．０について行い、節点２１個の位置が求められる。

図５に、例えば節点１と節点２の拡大図を模式的に示す。実際には節点３以降もあり節点ｍまで存在する。各節点において、膨張波が発生した前後で２点を設ける。例えば、節点１には膨張波発生前後で２点（１１，１２）あり、それら各々の圧力をｐ１１，ｐ１２とする。また、節点２には膨張波発生前後で２点（２１，２２）あり、それら圧力をｐ２１，ｐ２２とする。

まず、ステップＳ３０１において、節点１の圧力ｐ１を求める。形状先端の節点１の点１１における圧力ｐ１１は、円周周りの圧力とみなして節点１の圧力ｐ１とする。

ところで、仰角や横滑り角のない円錐周りの流れ場は軸対称となる。特に衝撃波が付着している状態においては、流体の特性は円錐頂点から延ばした直線に沿っては変化しない。このような流れは錐状流れ（ｃｏｎｉｃａｌｆｌｏｗ）と呼ばれる。錐状流れの解を求める式が次の（４）式であるＴａｙｌｏｒ−Ｍａｃｃｏｌｌの式である。
上記節点１の圧力ｐ１は、上記Ｔａｙｌｏｒ−Ｍａｃｃｏｌｌの式により求められる。

次に、ステップＳ３０２において、節点ｊ＝１とする。ステップＳ３０３において、節点１における膨張波発生後の圧力ｐ１２を、次に示すＰｒａｎｄｔｌ−Ｍｅｙｅｒの式（（５）式）により求める。

ここで、Ｍｎは点ｎにおけるマッハ数であり、δは対象となる点における流れの傾き角である。またγは比熱比であり、ここでは空気の比熱比１．４を用いる。

次にステップＳ３０４では、節点１の点１２における圧力勾配を求める。この圧力勾配は次式により求められる。

次に、ステップＳ３０５で節点２の圧力を求める。節点２の圧力ｐ２は、点２１における圧力ｐ２１を円錐台周りの圧力とみなして節点２の圧力とすることにより求められる。
ｐ３＝Ｐｃ−（ｐｃ−ｐ２）ｅ^−η
ｐｃは対象となる直線要素の傾きを半頂角とする円錐の圧力である。ηは次の（７）式で表せる。

ステップＳ３０６において節点すべてについて圧力を求められたか確認し、求められてなければ、ステップＳ３０７でｊ＋１をｊとして再びステップＳ３０３に戻り、節点ｊの膨張波後の点ｊ２における圧力ｐｊ２を求める。

このようにして、すべての節点における圧力を求め、図２に示すステップＳ２０６に移る。

ステップＳ２０６では、圧力分布算出部１３及び抗力算出部１４において、形状表面の各節点の圧力を表面積で積分して抗力を算出する。次のステップＳ２０７では、このようにして求めた抗力を評価関数とする。

なお、抗力は小さいほど好ましいので、評価関数としても値が小さいことが好ましい。しかし、常に小さいものだけを短絡的に採用すると、アルゴリズムが硬直化してしまい形状も限定されるので、全体として必ずしも常に好ましい形状に収束するとは限らない。

ステップＳ２０７Ｍで、マッハ数が２かどうか判定する。マッハ数が１の場合には、ステップＳ２０８Ｍでマッハ数を＋１すなわち、２としステップＳ２０４からの処理を繰り返す。

次のステップＳ２０８では、個体数判定部１６において個体数Ｐが５０を超えたか否かを判定する。個体数Ｐが５０を超えない限り、ステップＳ２０９に移り形状表面の節点を変え、ステップＳ２１０で個体数Ｐをプラス１（今の場合Ｐ＝２、すなわち第２個体）として、ステップＳ２０４に移る。

ステップＳでは、入力ファイル作成部１１において、この第２個体の場合の形状について圧力算出プログラムの入力算出ファイルを作成する。ここでも、形状を表現する式として、パラメトリック曲線の１つであるＢｅｚｉｅｒ曲線を用いる。

ステップＳ２０５では圧力算出部１２において、先に図３を用いたような計算を行い、第２個体の形状表面の各節点における圧力を算出する。ステップＳ２０６では、各節点の圧力を表面積で積分して抗力を算出する。次のステップで、評価関数決定部１５においてこの算出された抗力を第３個体の評価関数とする。

このように評価関数を取得し（Ｓ２０７）、形状表面の節点を変え（Ｓ２０９）個体数を１つ増やして（Ｓ２１０）、次の個体の評価関数を得る。形状表面の節点を変え、抗力を求め、評価関数を得る処理を、個体数Ｐが５０になるまで繰り返す。したがって、５０個体に対する評価関数（抗力）が得られる。

ステップＳ２０８で、個体数判定部１６で個体数Ｐが５０を超えたことが判定されると、ステップＳ２１２に移る。ステップＳ２１２では、世代数判定部１７において、世代数Ｇが５０を超えたか否かが判定される。世代数が５０になるまではステップＳ２１３に移り、次世代データ作成部１９において、現世代（Ｇ＝１）の次の世代データが作成される。

次世代データの作成は、遺伝的アルゴリズムを用いる。すなわち、５０個の個体の評価関数を調べ、形状に関して交叉、複写、突然変異の処理を行い、それらの形状について重み付けを行って次の世代の形状を作成する。

複写は、現世代の形状のコピーであり、交叉は父親の形状と母親の形状を織り交ぜて子孫を作る（新しい形状を生成する）処理である。突然変異は、形状の要素の一部をランダムに変えて新しい形状を生成する処理である。重み付けでは、それらの形状に係数を掛けて加算する。評価関数（抗力）が大きい形状は、係数をゼロとして、採用されないものもある。

ステップＳ２１４では、世代数を＋１としてこの次世代の第１個体とする。ステップＳ２０３に戻り、マッハ数１、２についての評価関数を求める。同様にして５０個の個体の評価関数を求める。遺伝的アルゴリズムを行って次世代の形状を求め、これらの処理を５０世代まで繰り返すことによって最終世代の評価関数を得る。ステップＳ２１５では、形状決定部２０において、最終世代の評価関数から最適な形状を決定する。

ところで、通常、飛しょう体頭部の先端は鈍頭であり、後端は円柱形状となっている。そして、形状の表現にＢｅｚｉｅｒ曲線を用いるとすると、次に示す（８）式を条件とする。

図６（ａ）（ｂ）（ｃ）に、評価関数とこの評価関数から最終的な形状を決定する様子を示す。図６（ａ）は、マッハ数が１のときのある世代の評価関数ｆ１に対応する形状を示し、図６（ｂ）はマッハ数が２のときのその世代の評価関数ｆ２を示す。横軸はｘ軸であり縦軸はそのｘ位置における半径ｒを示す。

図６（ｃ）は、ｘ軸にマッハ数１のときの評価関数の値をとりｙ軸にマッハ数２のときの評価関数ｆ２の値を取った評価関数の図である。この図には、各々第１０世代（１０ｔｈ）、第２０世代（２０ｔｈ）、第３０世代（３０ｔｈ）、第４世代（４０ｔｈ）、第５世代（５０ｔｈ）の値を示している。

評価関数は抗力を意味するので小さい値であることが望ましい。したがって、飛しょう体の速度がマッハ１やマッハ２のときがあるとすると、マッハ数が１のときの評価関数ｆ１もマッハ数２のときの評価関数ｆ２も小さくなることが望ましい。図６（ａ）と図６（ｂ）に示す形状のときの５０個の評価関数が各世代について得られる。

図６（ｃ）に着目すると、最初の世代（Ｇ＝１）の評価関数ｆ１，ｆ２を示す５０個の点のデータは、世代が進むに連れて評価関数ｆ１及び評価関数ｆ２が小さくなる方向に移動することになる。

そして、最終世代(第５０世代)においてｆ１及びｆ２が最も小さい値、すなわち図６（ｃ）で最も原点に近い点に対応する形状が、最適の形状を意味することになる。このようにして、図２のステップＳ２１５における処理である、飛しょう体頭部の最適化形状が決定される。

上記実施形態では、飛しょう体の速度が異なる、マッハ数Ｍが１と２の場合にどちらの評価関数も小さい値となる形状を最適な形状と決定する例について説明した。しかし本発明はこの場合に限られず、１つの速度において評価関数を求めこれらから最適な形状を得ることも可能である。また３以上の速度の場合の評価関数を求めそれらの評価関数が小さくなる飛しょう体頭部の形状を最適な形状と決定することも可能である。

また上記実施形態では個体数Ｐが５０個の場合について、また世代数が５０個の場合について説明した。しかし本発明において、個体数、世代数はこれらの値に限られない。

本発明の上記実施形態においては、形状を示す曲線としてＢｅｚｉｅｒ曲線を用いた場合について説明した。しかし、本発明で形状を示す曲線はＢｅｚｉｅｒ曲線に限られず、一般的にはパラメトリック曲線を用いることが可能である。

本発明は実施形態の限られず種々変形して実施可能であり、これらの変形例も本発明の技術思想が同一性を有する範囲内で本発明に含まれる。

１０・・・初期形状決定部、
１１・・・入力ファイル作成部、
１２・・・圧力算出部、
１３・・・圧力分布算出部、
１４・・・抗力算出部、
１５・・・評価関数決定部、
１６・・・個体数判定部、
１７・・・世代数判定部、
１８・・・形状データ作成部、
１９・・・次世代データ作成部、
２０・・・形状決定部。

Claims

パラメトリック関数により飛しょう体の頭部の初期形状を設定する初期形状決定ステップと、
この初期形状設定ステップにより設定された前記形状の節点のおのおのについて衝撃波・膨張波法により、圧力を算出する圧力算出ステップと、
形状表面の各節点の圧力を表面積により積分して抗力を算出し評価関数とする評価関数算出ステップと、
前記形状表面の節点を変えて各節点の圧力算出及び評価関数算出を繰り返して、所定の個体数の評価関数を得る評価関数取得ステップと、
この評価関数取得ステップにより得られた、評価関数の対応する形状について、交叉、複写及び突然変異の処理を行って、次世代の形状を形成する次世代形状形成ステップと、
この次世代形状形成ステップにより形成された次世代について、前記圧力算出ステップによる圧力算出、前記評価関数算出ステップによる評価関数の算出を行ってこの世代の所定の個体数の評価関数を取得し、これを所定の世代数について繰り返して最終世代の所定個体数の評価関数を得、適切な形状を決定する形状決定ステップと、
を有することを特徴とする、飛しょう体の形状決定方法。
前記評価関数取得ステップでは複数の速度における評価関数を求め、前記形状決定ステップでは上記複数の評価関数の値がいずれも小さいものを最適な評価関数とし、この評価関数に対応する形状を最適な形状と決定することを特徴とする請求項１記載の飛しょう体の形状決定方法。
前記パラメトリック曲線として、Ｂｅｚｉｅｒ曲線を用いることを特徴とする請求項１又は２記載の飛しょう体の形状決定方法。
前記圧力算出ステップでは、個体の最初節点の圧力を、Ｔａｙｌｏｒ−Ｍａｃｃｏｌｌの式により求めることを特徴とする請求項１ないし３のいずれか一項記載の飛しょう体の形状決定方法。
前記圧力算出ステップでは、第２以降の節点の圧力を、Ｐｒａｎｄｔｌ−Ｍｅｙｅｒの式により求めることを特徴とする請求項１ないし４のいずれか一項記載の飛しょう体の形状決定方法。
パラメトリック関数により初期形状を設定する初期形状決定手段と、
この初期形状設定手段により設定された前記初期形状の節点のおのおのについて衝撃波・膨張波法により、圧力を算出する圧力算出手段と、
形状表面の各節点の圧力を表面積により積分して抗力を算出し評価関数とするする評価関数算出手段と、
前記形状表面の節点を変えて各節点の圧力算出及び評価関数算出を繰り返して、所定の個体数の評価関数を得る評価関数取得手段と、
この評価関数取得手段により得られた、評価関数の対応する形状について、交叉、複写及び突然変異の処理を行って、次世代の形状を形成する次世代形状形成手段と、
この次世代形状形成手段により形成された次世代について、前記圧力算出手段による圧力算出、前記評価関数算出手段による評価関数の算出を行ってこの世代の所定の個体数の評価関数を取得し、これを所定の世代数について繰り返して最終世代の所定個体数の評価関数を得、適切な形状を決定する形状決定手段と、
を有することを特徴とする、飛しょう体の形状決定装置。
前記評価関数取得手段では、複数の速度における評価関数を求め、前記形状決定手段では上記複数の評価関数の値がいずれも小さいものを最適な評価関数とし、この評価関数に対応する形状を最適な形状と決定することを特徴とする請求項６記載の飛しょう体の形状決定装置。
前記パラメトリック曲線として、Ｂｅｚｉｅｒ曲線を用いることを特徴とする請求項６又は７記載の飛しょう体の形状決定装置。
前記圧力算出装置では、個体の最初の節点の圧力を、Ｔａｙｌｏｒ−Ｍａｃｃｏｌｌの式により求めることを特徴とする請求項６ないし８のいずれか一項記載の飛しょう体の形状決定装置。
前記圧力算出装置では、第２以降の節点の圧力を、Ｐｒａｎｄｔｌ−Ｍｅｙｅｒの式により求めることを特徴とする請求項６ないし９のいずれか一項記載の飛しょう体の形状決定装置。