Comments
Description
Transcript
橡 Seminar 97-2
7. 2 翼列翼型の設計と多目的遺伝的アルゴリズム 今年 1997 年は、英国の Sir F. Whittle のジェットエンジンが初めて自立運転に成功してから 60 年 目にあたる年である[26]。 その間にたゆまない性能 向上への努力は継続的に 続けられてきた。第一 世 代ターボジェット、第二 世代低バイパス比ターボ ファン、第三世代高バイ パス比ターボファンと発 達 してきたジェットエンジ ンは、今や次世代高バイ パス比ターボファンから 超高バイパス比ターボフ ァ ンの時代に入ろうとして おり[27]、また騒音や排気 を低減させる技術も進歩 し、環境に優しいエン ジ ンとしても発展を遂げて いる。このようなガスタ ービンエンジンの発達の 重要な鍵の一つとなって い るのが圧縮機である。と いうのは、圧縮機がエン ジンの圧縮比や流量を支 配するため、エンジンの 特 性に絶対的な影響を及ぼ し、さらに制御するのに 最も難しいエンジン要素 の一つだからである[11] 。 これからのエンジン開発 には、そのコンパクト化 や軽量化と同様に低燃費 化のため、エンジン圧縮 比 を上げ、圧縮機の効率を 改善しようとする強い動 機づけがある。従って、 このような要求に備える た め、圧縮機の研究や開発 において多大な努力が注 がれており、航空機エン ジン産業界では長い間、 幅 広く深い研究開発プログラムを遂行してきている。 このようなガスタービンエンジンには 2 つの型の圧縮機が広く使われてきた。軸流型と遠心型であ る。軸流型では、主とし てエンジンの軸流方向に 空気が流れる。動静翼列 で段を形成し、流量が大 き く、高効率の可能性をも ち、きわめて高度に発達 している。軸流型は大型 の航空機エンジンで支配 的 となっている。一方、遠 心型では、空気は圧縮機 のロータの中で軸方向か ら半径方向に曲げられる 。 比較的高い段圧縮比が可 能であるため全長が短く 、構造的にも簡単である 。遠心型は小型の航空機 エ ンジンや多くの産業用で使われている。 ガスタービン設計においては、産業用および航空 用を問わず低燃費化が一つの要求項目になって い る[28]。そこで、低燃費化 を図るには、熱サイクル の効率を増すことよりも 全体圧縮比および要素 効 率を高めることが必要と される。さらに航空用で は軽量化も重要な要求項 目である[28]。こうした 高 圧縮比化をより少ない段 数で達成する、すなわち 負荷を高くすることが望 まれるのは当然のことで あ る。 このようなことから、 エンジンの軽量化・低燃 費化という観点において 圧縮機に対する代表的な 要 求項目として、「高圧縮 比化、高効率化、高負荷 化」の３つが挙げられる 。しかしながら、圧縮機 翼 列内の流れ（図 7.13）[29]を考えると、これらを同時に満足させるのはきわめて困難な状況にある。 翼列内では逆圧力勾配の ために境界層が発達しや すく、流入速度が増すと 衝撃波も発生してしまう 。 さらに二次流れやクリア ランス流れなどの流れ場 の三次元性が強くなるた めにかなり複雑な流れ場 と なり、最終的に性能低下 を招く可能性が高くなる 。従って、このような境 界層の発達、衝撃波、そ し て三次元性の流れ等の影響を考慮した設計手法が重要となる。 1970年代には、新しく圧縮機を設計・開発するためにかなりの試行錯誤が成されたが、今や状況は 進展傾向にある。エンジ ンメーカーが新しい圧縮 機を設計するとき、過去 に経験したものと多少の 差 はあってもあまり大きく 違っていなければ、設計 どおり上手くいくことに 疑いをもつことはなくな っ ている。その背景として、近年のCFD技術の急激な進展とコンピュータ技術の著しい進歩が挙げられ る。CFD を設計ツールとし、広範囲な実験データベースとあわせ、圧縮機の流体力学的挙動やその構 造的完全さに対する予測 精度がかなり向上してい る。実際、圧縮機の翼間 内部流動、すなわち翼列 に 関する数値解析を行った研究も多数報告されている。例えば、2 次元翼列の解析[30-32]、3 次元翼列 の解析[33,34]、あるいは動 翼やフラッターを扱っ た非定常解析[35,36]などが 挙げられる。こうし た 多数の翼列解析からみても分かる通り、圧縮機の性能は翼列の性能に左右され、高性能な圧縮機を実 7-22 Inflow Clearance flow Shock wave Clearance vortax Secondary flow Separated flow Secondary flow Corner vortex Rotation Boundary layer Outflow 図7.13 圧縮機翼列内での流れ場の概念図 現するには、高性能な翼 列を構築する必要がある 。特に、二次元翼列に関 する研究開発は最も重要 な 空力的課題の一つであり 、二次元設計法が非常に 重要な位置を占めている 。なぜならば、翼列の空 気 力学的設計では全体の性 能面での設計仕様を考え て、まず翼間流路を決定 しなければならない。こ の ときに性能の良い、つま り圧縮機の代表的要求項 目をできるだけ満足する ような二次元形状を求め 、 最終的にこれらの翼型を スパン方向に積み重ねて 三次元翼列形状を得てい るからである。従って、 二 次元翼列の性能が三次元 形状、もしくは圧縮機全 体の性能に大きく影響し てくる。すなわち、二次 元 翼列設計は圧縮機を設計する上で必要不可欠なものといえる。 翼列翼型の空力設計のほとんどは逆設計法に基 づいている[37]。逆設計法とは予め与えられた圧 力 分布を満たすような形状 を求める方法である。圧 力分布から積分方程式を 解いて直接形状を求める 方 法と、ある形状の圧力分布をCFDで求めてから予め与えられた圧力分布になるように形状を変化させ ていく方法がある。適切 に圧力分布を設定できる 知識を持つ設計者にとっ てこの方法は非常に効果 的 ではあるが、圧力分布を 設定することは一般に難 しい問題である。設定さ れた圧力分布によっては 物 理的に不可能な形状になる場合もありうる。さらに、形状の制約条件も複雑になる。 そこで、本研究では圧 縮機の二次元翼列翼型に 注目し、圧縮機の要求項 目について順設計法に基 づ いた空力形状の最適化を 試みる。順設計法（ある いは順最適化法）とは、 通常の逆設計に対し、形 状 を直接変化させることによ って与えられらた空力的 な目的関数を最大もしく は最小にする方法である 。 ここで、前節を振り返ってみると、3 つの最適化手法（勾配法・焼きなまし法（SA）・遺伝的アル ゴ リズム（GA））を適用した低速の高揚力翼型の空力最適化を行い、この最適化問題に対する適応性に ついて3手法の比較を示した。その結果、3手法のなかでGAを用いた場合が最も高い揚力を発揮する 翼型形状を求めることができ、この空力最適化問題を解くのにGAが最も有効な方法であることが分か 7-23 った。従って、ここで行う翼列翼型の空力最適化問題を解くための手法としてGAを採用する。 翼列にはその用途と流 入マッハ数の範囲に応じ て様々なタイプの形状が 使われている。従来の翼 列 設計では、翼列条件と翼 型特性との相関や翼列配 置において優れた特性を もつ翼型の設計資料の作 成 の目的で、広範囲にわた る系統的な実験的研究が 圧縮機翼列について翼列 風洞により遂行されてい る [38]。亜音速圧縮機翼型に適したNACA65系列翼型（図7.14）や遷音速圧縮機翼型に適した多重円弧 翼（Multiple circular arc airfoil）（図7.15）が代表的である。多重円弧翼とは翼列内部流路における スロート部に相当する位 置で二分し、前縁側と後 縁側とで滑らかに接続す る２つの円弧によって反 り 線が構成されている翼型である。 一方、高速域でも高い性能を発揮させるために、より自由度が高く最適化可能な翼型がCFDを利用 して研究・ 開発され ている。 その代表 的なもの が超臨界 翼列翼型 [39,40]や拡散 制御翼（ Controlled Diffusion Airfoil, CDA）[41-43]である（図 7.16）。CDA は、臨界マッハ数を越える流入条件に対 し、超音速から亜音速へ衝撃波を伴わずに減速し、境界層が剥離しないようなマッハ数分布（図7.17） を与え、これを満足する ように順解法または逆解 法により得られた翼型で ある。従って、損失が小 さ く抑えられていることになる。CDA の形状は前縁付近でキャンバーの曲率が大きく、前縁・後縁の翼 厚が大きくなっている。 本研究では、遺伝的アルゴリズム（GA）を用いて亜音速圧縮機翼型の空力最適化を試みる。実際に 最適化された翼型がどの程度良いものであるかを確かめたい。そこで、CDA を比較の対称翼型として 取り上げ、得られた最適翼型とCDAとの性能比較も行う。 7. 2. 1 数値計算法 翼列翼型の性能値はCFDにより得られる。圧縮機翼列内では逆圧力勾配になるため、翼面上での境 界層の発達は避けられず 、剥離も生じる可能性も ある。また、翼列の性能 、特に損失に関しては境 界 層の発達と深く係わって いるので、計算で境界層 を上手く再現させないと 正しい性能値を評価する こ とができない。そこで、翼列内の流れ場は二次元圧縮性 Navier-Stokes 方程式を基礎方程式として解 析する[44-47]。 対流項は、Wadaらの修正HLLE法[44]を適応し、粘性項は２次の中心差分で近似する。時間方向に は、LU-SGS 法[46]を用いる。この方法は、実質的に行列反転のない、無条件安定な陰解法であり実 用計算によく用いられている 。また定常問題を扱ってい るので、局所時間法（Locally varying time step）を使って収束を早めている[44]。 高レイノルズ数流れにおいて、Thin-layer 近似を用い、乱流現象を数値計算で捕らえるために乱流 モデルを利用する。本研究では乱流モデルとして Baldwin-Lomax の代数モデル[47]を適用し、翼型 において境界層は前縁から乱流としている。 7-24 図7.14 NACA65系翼 図7.15 多重円弧翼 図7.16 拡散制御翼（CDA） 7-25 図7.17 CDAの設計指針 7-26 7. 2. 2 翼列翼型の空力最適化法 圧縮機には作動マッハ 数に応じて様々なものが あり、翼型の形状もそれ ぞれ変わってくる。また 、 圧縮機は動翼（Rotor）と静翼（Stator）により段が形成され、動翼内および静翼内での流れ場は全く 異なるので、両者の形状 も変わる。本研究では亜 音速の静翼に関して最適 化を行い、圧縮機の代表 的 な要求項目「高圧縮比・高効率・高負荷」を目的とした。初めに、3つの項目を別々に評価し、それぞ れに関して空力 的に最適な翼型 形状を求める単 目的遺伝的アル ゴリズム（Single-Objective Genetic Algorithm, SOGA）を試みる。次に、3 つの項目を全て同時に評価し、どの要求性能もよくなるよう な最適形状を求める多目的遺伝的アルゴリズム（Multi-Objective Genetic Algorithm, MOGA）を行 う。MOGA の場合は、他のどの解よりも劣っていない解の集合（パレート最適解）[48]を求めること を目標にしている。そして、パレート最適解の中から最終的に解を 1 つ選ばねばならないが、これは 設計者の意図に反映するものが最適解として決定される。 1. 翼列翼型の形状定義 求める翼列翼型の形状 を定義する。翼型は一般 的に翼厚分布（対称翼型 ）と反り分布（キャンバ ー 線）で構成されている。というのは、翼型の性能がこの 2 つの要素でほぼ与えられるからである。従 って、本研究では対称翼 型とキャンバー線を別々 に形成し、これらを組み 合わせることにより翼型 形 状を決定する（図7.18）。まず、対称翼型の上下面上に14点、キャンバー線上に7点の制御点をそれ ぞれ分布させる。このと き、対称翼型では最大翼 厚となる制御点、キャン バー線の場合は最大キャ ン バーとなる制御点をはじ めに決定し、形成された 対称翼やキャンバー線が 交差しないような制約の も とで残りの点は順次その間に位置するように決め、これらの制御点の座標(x,y)を設計変数とする。 つ まり、1つの制御点はxとyの2つの設計変数から成っている。ただし、前縁・後縁の座標はそれぞれ (0,0), (1,0)に固定する。また対称翼型の場合、上面の制御点が決まれば、これを x 軸について対称移 動させると下面の制御点も自動的に決まる。さらに、後縁点の真上の制御点（x座標1に固定）はCDA のように後縁に厚みのある翼型も表現できるようにするために設けた。結局、全21個の設計変数とな る。 対称翼型およびキャンバー線を表現す る曲線は、制御点を三次 B-spline（Basis-spline）関数に よ って決定する[49]。さらに、三次 B-spline 曲線は制御点を通らないので、指定した翼弦長と最大翼厚 となるように修正を施している。 こうして形成された対称翼型とキャンバー線を組 み合わせることによって求める翼型形状を定義 す る。 2. 単一目的遺伝的アルゴリズム 圧縮機に対する代表的 な要求項目「高圧縮比・ 高効率・高負荷」につい て、各項目ごとに最適化 を 行う（SOGA）。このとき、各世代中の個体数および最大世代数を次のように設定する。 個体数＝64 最大世代数＝80 また、各操作は次のように行う。 7-27 図7.18 翼型の定義 7-28 ・コード化 翼型を定義するときの制御点の座標を設計変数とする。このとき設計変数を 2 進数には変換せず、 実数型データのまま遺伝子として扱う。従って遺伝子の数は合計21個となる。そして、この遺伝子で 表現できる翼型の中から無作為に64個体選び出し、初期個体集団を形成する。 ・評価 性能を評価するとき、翼列内の流れ場は二次元圧縮性 Navier-Stokes 方程式に従う。用いる計算格 子や流れ場の流入条件等は以下のように設定する。ただし、流れ場の条件はCDAに対してOff-design 条件になっている[50]。 ・計算格子：翼間流路のH型格子(201×65) ・流入マッハ数：M1 ＝0.25 ・流入角：β 1 ＝40.0 [deg.] ・入口全圧：Pt1＝1.043×105 [Pa] ・入口全温：T t1＝294 [K] ・レイノルズ数：Re＝0.70×10 6 ・翼列ピッチ：s/c＝0.5988 ・食い違い角：γ＝14.4 [deg.] また境界条件は、 ・流入境界：流入角β 1 を固定し、全温を与え、他の物理量は内点から外挿する。 ・流出境界：流出マッハ数（M2 ＝0.195）を与え、自由流出とする。 ・翼面上：滑り無しの条件とする。 ・格子の接合境界：周期境界を与える。 以上の流れ場の条件は、 文献[50]を参考にしたもの である。また、機械学会 においてもこれと同様 な 条件を用いた翼列解析のワークショップが行われている[51]。 GAを用いた最適化において、最大の問題は計算時間が非常に多いということである。全計算時間の うちほとんどの部分を占 めるのがこの評価である 。従って、より効率的に 、より短時間で評価する こ とが望ましい。そこで、 各個体を同時に全てを評 価してしまうことを考え る。つまり、個体評価の 並 列処理である。このため、航空宇宙技術研究所の所有する高性能並列計算機 NWT（Numerical Wind Tunnel）を利用する。 並列計算機NWTの概要 NWT は 166 台の要素計算機（Processor Element, PE）から成り、各 PE がベクトル演算機を持っ ている。1PE の性能はピークで約 1.7GFLOPS（主記憶 256MB）と非常に高い。また PE 間のデータ 転送に関しても柔軟でありかつ速いネットワークを採用している。 並列 計算を 行う際 には並 列プ ログラ ミング をしな ければ なら ない。 このた め NWT 独 自の NWT Fortran を 用 い る 。こ れ は 通 常 の FORTRAN に 、 DO ルー プ の 分 割 、同 期 に 関 する 翻 訳 指 示 行 7-29 （Directive）を挿入した並列化を指示する言語である。単一のプログラムで記述でき、並列化は簡 潔 かつ容易である。また、 並列化したコードをその まま普通のベクトル計算 機にかけることもできる 。 ここで、スーパーコンピュータについての知識をさらに深めるために文献[52]を挙げておく。 このNWTを用い、各PEに各個体の評価を割り当てる。従って、64台のPEを使用することにより、 並列化してない場合と比べるとGAの評価部分の計算時間が 1／（個体数）、つまり1／64となり 大幅に時間を減らすことが可能である。NWTを利用した計算例は他にも多くあり、例えば文献[53]で は、動静翼列内の非定常流れ解析に関する研究が報告されている。 ・選択 GA の選択法にはいくつかの手法があるが、ここではランキング選択法を採用する。この方法では、 個体集団中における各個体の適応度の大きな個体から順に1, 2, 3,…とランク付けをし、このランクに 応じて各個体の選択確率が決まる。つまり、ランク i の個体の選択確率は Pi = c ⋅ (1 − c )i −1 (7.24) で与えられる。これは、Michalewiczが提案した非線形ランキング選択であり[54]、c はランク1の個 体が選択される確率を表している。ここでは、c = 0.075 と設定した。この値は大きすぎると、ランク 1の評価関数値の高い個体が選択されやすくなるので、その影響で全体の個体が同じものばかりになっ てしまい、逆に小さすぎると、ランク 1 の個体の選択確率も小さくなるので、いつまでたってもよい 個体が増えない。従って、適切な値を設定する必要がある。 このように、ランキン グによる選択確率は直接 適応度によって決まるわ けではなく、設計者の意 図 に基づいて考慮された c によって(7.24)式から決定されるので、選択確率を自由に設定することができ る。さらに、(7.24)式の代わりに別の方法を用いる場合も考えられ、例えば線形の式を用いた線形ラン キング選択法がBakerにより提案されている[55]。 ランクによる選択確率の決定後、この確率の値を各個体の適応度として再定義し、前節のGAで扱っ たルーレット式選択法を用いて親を選択する。 ・交叉 全ての遺伝子を対象にして、2 つの親から互いに対応する遺伝子を 50%の確率で交換する方法を採 る（図 7.19）。つまり、制御点の座標を遺伝子としてるので x 座標どうし、y 座標どうしをそれぞれ 独立に交換する。 ・突然変異 全ての遺伝子に対して突然変異が起こる可能性を 調べ、発生した場合にはその遺伝子に擾乱を加 え るという操作を行う。加える擾乱としては、x座標では-0.2∼+0.2の範囲、y座標では-0.02∼+0.02 の範囲の量をランダムに 決める。また、突然変異 が起こる確率（突然変異 率）は一般的に小さい。 あ まり大きすぎるとランダムサーチになってしまうからである。そこで、突然変異率は1世代目で20%、 100世代目では1%になるように線形的に小さくなるように設定する。これは、世代を重ねるごとに集 団中に良い個体が増えてくるので、その個体をできるだけ破壊しないようにし、解の収束をよくする 7-30 (x 1 1, Parent 2 (x 2 1, y 12 , Offspring 1 (x 2 1, y 11 , (x 1 1, y 12 , Parent 1 Offspring 2 ) (x y 11 , 1 2, ) (x 2 2, ) ( ) ) ( ) y 12 , LL x i1 , yi1 , LL y 22 , LL x i2 , y i2 , LL ) (x , y ), LL(x 1 2 ) (x Probability 50 ％ 2 2, 2 2 2 i , ) ( ) yi2 , LL ) y21 , LL xi1 , yi1 , LL 図7.19 交叉方法 ためである。 ・エリート戦略 ここでは、集団中で最も評価の良い個体と2番目に良い個体を優先的に次世代に保存する。 3. 多目的遺伝的アルゴリズム SOGA では「高圧縮比化・高効率化・高負荷化」の各項目を別々に扱い、1 つの評価関数のみを考 慮した。最も望ましいの はこれらの要求を全て満 足させることであるが、 非常に難しい問題である こ とは既に述べた。しかしながら、MOGA を適用することによりある程度は実現可能である。一般的に 静翼の場合、転向角が増加すると圧縮比も増加するので実質的に圧縮比と効率を高める 2 目的最適化 と考えられる。ところが 、高転向角は境界層剥離 による大きな損失を招く ので圧縮比低下にもつな が る。従って、転向角を大きくするには限界がある。そこでこのパラメータも重要であると考え、SOGA の 3 目的を同時に最適化する MOGA を試みる。MOGA は SOGA と異なり、ベクトルの評価関数を最 適化し、得られた解は単 一ではなくパレート最適 解として知られる解の集 合である。元来、パレー ト は経済学者であり、パレ ート最適は彼の「財が最 も効率的に分配されてい る状態」についての概念 か ら始まった。このとき、 「他の誰かの状態を悪化 させない限り、誰の状態 も改善できない。」[56] 。 これを工学の多目的最適化に対応させたパレート解は、1つの項目を良くするのに残りの2つの項目の うち少なくとも1つが悪くならないと改良できないという意味で最適状態にある。 現在ではGAを用いて多目的最適化問題を解いている研究も報告されている[57-59]。この中で共通 してパレート最適集合を 求める重要な鍵を握って いるのがパレートランキ ング法とシェアリング法 で ある。パレートランキング法は Fonseca と Fleming によって提案されており[60]、他にも Goldberg により考案されたランキング法[8]もある。また、シェアリング法は解が一カ所に集まってしまうこと を防ぎ、解のばらつき（多様性）を保つために用いられる[55]。 7-31 MOGAの設計条件を設定する。コード化、流れ場条件、交叉、突然変異はSOGAの場合と同じであ るが、ここでは選択とエ リート戦略が異なる。こ れらのオペレータにはパ レートランキングとシェ ア リング[55]、さらにベストN選択法[61]を採用した。 各世代中の個体数および最大世代数を以下のように設定する。 個体数＝64 最大世代数＝80 ・パレートランキング GAでは個体集団を用いた探索が進められるので、ある程度評価の良い個体を同時に保ちつつ探索が 可能で、パレート最適集 合を同時に求めることが できる。このために世代 内のパレート最適な個体 を 選択するように拡張されたランキング選択法がある。これには次のような定義を使う。 ［定義］： 集団中に xi ( ) = (x i , y i ), x j = x j , y j が存在し、また最大化すべき評価関数の組を F = (F1 , F2 ) と したとき、 1. 以下の式を満たすとき、 xi は x j に支配されている（ x j の劣解）という。 ( ) ( ) ( ) F1 (xi ) ≤ F1 x j かつ F2 (xi ) ≤ F2 x j 、ただし F(x i ) ≠ F x j (7.25) 2. (7.25)式を満たす x j が存在しない、すなわち ( ) ( ) F1 (xi ) > F1 x j または F2 (xi ) > F2 x j (7.26) のとき、 xi は x j に支配されない（ x j の非劣解）という。つまり、評価関数空間の実行可能領域内 に ある非劣解の集合がパレート最適となる。 個体 xi 、世代 t において、 pit 個の個体で xi が支配されているとき、個体 xi のランクは rank(xi , t ) = 1 + pit (7.27) で決定する。非劣解の個体は全てランク1である。例えば図7.20で示されているような各個体のパレ ートランキングが決まる。 ・シェアリング MOGAを解くためには解がパレート最適であるのと同時に、パレート最適解がパレート集合中に一 様に分布している必要が ある。つまり、解の多様 性を保つことが重要でな る。しかし、確率的な選 択 を行うことによりある特 定の評価の良い個体のみ が選ばれるような偏りを 生ずることがある（遺伝 的 浮動）[55]。この遺伝的浮 動を避けるために、多く の点が集中している部分 での適応度を意図的に 下 げ、孤立している点の適 応度を大きくすることに より適応度を修正する操 作が必要になる。これを シ ェアリングという。シェアリング後修正された F( xi ) を F ′( xi ) としたとき 7-32 1 f2 1 2 1 3 5 2 1 f1 図7.20 パレートランキング F ′( xi ) = ∑ F (x i ) s(d ( xi , x j )) (7.28) j となる。ここで、s は解がどの程度同じ場所に集中しているかという度合いを決めるシェアリング関数 という。もちろん様々なシェアリング関数を考えることが可能であるが、例えば次式で表される[55]。 d α s(d ) = 1 − σ share 0 σ < σ share (7.29) σ ≥ σ share d = d (x i , x j ) は個体 xi , x j 間の距離を表す。一般的に距離 d の測定法により、遺伝子のハミング距 離を用いる遺伝子型シェ アリングとパラメータ空 間でのユークリッド距離 を用いる表現型シェアリ ン グとがある。MOGA では評価関数空間でパレート最適集合を求めるので、後者がよく用いられる。こ こでは、設計変数のユークリッド距離を測定する。また、この方法には新しいGAのシェアリングパラ メータが導入される。ニッチサイズ σ share とべき 指数 α である。ニッチサイズ σ share は、個体集団 の 類似個体の中で互いにど の程度近くにいるときに 評価を下げるかを見積も る定数で、各評価関数の 個 別の最大値と最小値を利用する方法[60]で決定する。べき指数 α はニッチカウント mi = ∑ s(d (x , x )) i j (7.30) j を強めるために0.25とした。 本研究では、各ランク内の個体評価をスケーリングするため、ニッチカウント mi をランクに反映さ 7-33 せ、(7.24)式の確率分布により決まった選択確率に基づいて評価関数を再定義する。このようなシェア リングは 1) 初期収束の防止、2) 多様性の維持、3) 局所解の分布の検出などの目的で SOGA でも利 用されている[55]。 ・ベストＮ選択法 この方法は次世代での個体集団を形成するときに、現世代における N 個の親と N 個の子孫を合わせ て 2N 個の中で、親と子孫との区別なく評価の最も良い個体から順に N 個選ぶという方法である。従 って、常に良い N 個体が次世代に継がれていくので、いわばエリート戦略も取り入れた選択法といえ る。ここでは、ランクの シェアリングを適用して いるため評価関数を考え るのではなく、各個体の ラ ンクを考慮し、ランク1である個体から順に選択する方法をとった。 7. 3 翼列翼型の空力最適化 GAを用いて空力最適化問題を解いた場合、得られた解が本当に最適解に近いものであるかを調べる 方法がないことが多い。従って、GA による最適化が成功したかどうか判断し難い。そこで、GA の最 適解の探索能力および翼型の表現定義の妥当性を確認するためにCDAの圧力分布を利用して、まず逆 問題を解く。ただし、この問題に限りGAの選択法として基本的なルーレット式選択を用いている。 7. 3. 1 CDAの逆設計 はじめに、CDA を用いた翼列内の流れ場を解き、CDA の Cp 分布を求め、これを目標 Cp 分布とす る。もし逆問題が解けたとすれば、得られた翼型形状はCDAになるはずである。従って、逆解法によ りCDAに近い形状が求まればGAによる最適化は成功したことを意味する。 GAの最適化により逆問題を解くために、目標Cp分布 る Cp 分布 (CpCDA ) とGAにより形成された翼型に対す (CpGA ) の作る面積 (SCp )（図7.21）を最小にする問題として定式化する。 逆問題の定式化 面積 S Cp = ∫ Trailing Edge Leading Edge Cp CDA − CpGA dx →最小化 評価関数： F ( = exp − SCp ) (7.31a) (7.31b) (4.1b)式において、F が大きくなるなることは面積 S Cp が小さくなることであり、逆問題を解くことが できたことを意味する。この方法は文献[62]を参考にした。 図7.22に最適化の収束過程を示す。世代とともに F が増加しており、逆問題の解を探索している様 子が分かる。80 世代で最終的に得られた翼型とその Cp 分布を図 7.23 に示す。ここでデカルト座標 (x, y ) を食い違い角 γ だけ回転移動させ、前縁 (ξ c , η c ) = (0 , 0 ) 、後縁(ξ c , η c) = (1, 0 ) とな る座標に変換した。前縁付近でわずかな Cp 分布の差が見られるが、逆設計で得られた翼型形状はCDA と非常によく一致しているといえる。従って、GAにより逆問題を解くことは可能であり、最適解の探 索能力が確認され、翼型の表現方法も適切であることが分かった。 7-34 図7.21 目標の Cp とGAにより作られた翼型の Cp とが作る面積 図7.22 逆問題の収束過程 7-35 図7.23 逆設計で得られた翼型と Cp 分布 7-36 7. 3. 2 翼列翼型最適設計の定式化 CDA の逆設計を通して、GA の最適解の探索能力と翼型表現の妥当性が確かめられた。そこで、は じめに圧縮機の要求項目の1つである「高圧縮比化」について定式化を行う。 高圧縮比化の定式化 圧縮比 P2 P1 →最大化 制約条件：最大翼厚比 t 評価関数： F1 c = (t c )CDA (7.32b) ≥ SCDA (7.32c) V1 ≥ 072 . (7.32d) ：翼型面積 S ：速度比 V2 (7.32a) ∆S − ∆S 2 ∆VR − ∆VR = ( P2 P1 ) × exp × 10 + 2 2 (7.32e) (t c)CDA = 007 . , S CDA = 0047 . c2 ∆ S = S − S , ∆ VR =V V − 072 . CDA 2 1 ここで、 P1 , P2 , (t c )CDA , S CDA はそれぞれ入口静圧、出口静圧、CDA の最大翼厚比、CDA の面積 を表している。(7.32b)式はCDAとの性能比較のため、(7.32c)式は滑らかな翼厚分布のために取り入 れる。(7.32d)式は de Haller の方法[63]という。圧力上昇は翼列内の減速により得られるが、減速の 程度が大きくなると境界 層が発達し、流れが壁面 から剥離して損失が急増 し、性能が劣ってしまう こ とになる。そこで、入口速度 V1 と出口速度V 2 の比V 2 V1 の制約を設ける。0.72 は翼面上の境界層剥 離を抑える経験的な限界値である。 次に、「高効率化」について定式化を行う。効率 が良いということは損失が少ないということと 同 等なので、次のように定式化される。 高効率化の定式化 全圧損失係数 ω = Pt 1 − Pt 2 →最小化 Pt 1 − P1 制約条件：最大翼厚比 t ：翼型面積 S 評価関数： F2 ここで、 c = (t c )CDA ≥ SCDA ∆S − ∆S = −Log(ω ) × exp × 10 2 (7.33a) (7.33b) (7.33c) (7.33d) Pt1 , Pt2 はそれぞれ入口全圧、出口全圧の流量平均値を表す。(7.33d)式では、この問題が ω の最小化問題のため、 − Log(ω ), (ω < 1) を用いることに より最大化問題に変換している。つまり 、 ω が小さくなることは − Log(ω ) が大きくなることを意味する。 7-37 最後に、「高負荷化」について最適化を行う。翼面負荷は翼列翼型の能力度を示す重要な一尺度で、 作動の最適状態や失速限 界を表す[38]。そこで、負 荷パラメータの一つであ る転向角の最大化問題 を 次のように定式化する。 高負荷化の定式化 転向角 ∆β = β 1 − β 2 →最大化 制約条件：最大翼厚比 t ≥ SCDA ：翼型面積 S 評価関数： F3 ここで、 β 1 , c = (t c )CDA ∆S − ∆S = ∆β × exp × 10 2 (7.34a) (7.34b) (7.34c) (7.34d) β 2 はそれぞれ流入角、流出角を表す。 多目的最適化問題の定式化 圧縮機の三大要求項目 「高圧縮比化・高効率化 ・高負荷化」に基づいて 多目的最適化問題を定式 化 する。 1. 圧縮比 P2 P1 →最大化 （制約条件：速度比 V2 V1 ≥ 072 . ） (7.35a) (7.35b) 2. 全圧損失係数 ω →最小化 (7.35c) 3. 転向角 ∆β →最大化 (7.35d) また、1∼3全ての項目に適用される制約条件は 1. 最大翼厚比 t c = (t c )CDA 2. 翼型面積 (7.35e) S ≥ SCDA (7.35f) である。各評価関数には(7.32e, 33d, 34d)をそれぞれ用いる。 7. 3. 3 MOGAによる最適化結果 SOGA の場合と同様、個体評価を並列化して処理したので、計算時間は最大世代数 80 までで約 27 時間要した。図7.24に最適化の収束過程を示す。これは、各世代の集団中におけるランク1の個体数 を表している。世代をおうごとにランク 1 の個体数は増えており、70 世代を越えてからは全個体がラ ンク 1 に収束した。つまり、個体全部がパレート最適解になったわけである。そこで、80 世代目のパ レート最適解の分布（パレート面）を図7.25に示す。各軸は(7.32e, 33d, 34d)式の評価関数の最大値 と最小値を用いて線形的にスケーリングされており、値が1に近いほど評価が良いことを表す。 7-38 Fi scale = Fi − Fi Min , i = 1~ 3 Fi Max − Fi Min (7.36) 図に示されているように、パレート集団は 2 つのグループに分けることができる。圧縮比がより高い 集団と転向角がより高い 集団である。一般的に転 向角が増加すると圧縮比 も増加するが、転向角が あ る地点を越えると圧縮比 は減少してしまう。これ は圧縮比−転向角の図か ら分かる。さらに全圧損 失 係数−転向角の図より、 低損失を実現するには転 向角を小さく抑える必要 があり、両者の性能を上 げ るのは非常に難しい問題である。 図7.26に各性能によるパレート最適解の翼型形状の変化を示す。圧縮比が増すにつれて前縁・後縁 は鋭くなる傾向にある。 また、転向角が増加する と後縁に向かって大きい キャンバーを持つように な り、一方、損失が減少す るにつれて後縁付近のキ ャンバーはほとんど零と なっている。このことか ら も損失−転向角の両者の性能を良くするのは容易でないといえる。 得られたパレート最適解の中から最終的に最適解を 1 つ決定しなければならないが、これは意思決 定主体である人間であり、多目的意志決定システム全体として考えると、その方法論が必要となる[64]。 しかしながら、ここでは 特にそのような方法論は 取り入れず、パレート最 適解の中で比較的圧縮性 能 が良く、効率的な翼型形状を選んだ。この形状と Cp 分布をCDAと比較して図7.27に示す。ここでは この翼型をMOGA翼型と呼ぶことにする。MOGA翼型は圧縮比を増すように比較的鋭い前縁と、CDA のように損失を減らすように鈍い後縁を持つ。 最後に、得られたパレート最適解の中で、各項目において最も良い解と MOGA 翼型、さらに CDA との性能値を示す（表 7.1）。圧縮比が大きい翼型では、全ての項目において CDA より良い性能を示 した。MOGA 翼型は圧縮比・損失では CDA より勝っているが、転向角では若干低い値となった。こ れより、圧縮比と転向角とが完全なる比例的な関係ではないといえる。全体の性能を考慮しても、3つ の性能項目は相反する関係になっていることは一目瞭然である。 以上、新しい多目的遺 伝的アルゴリズムを圧縮 機の翼列翼型の空力最適 化に適用し、圧縮機で代 表 的な「高圧縮比・高効率・高負荷」の 3 大目標を掲げて形状最適化を行った。Fonseca-Fleming のパ レートランキングとシェ アリングに基づき、他の 解より劣らない解の集合 、つまりパレート最適解 を 得ることができた。その結果、3 つの性能において全て CDA よりも勝る翼型がパレート最適解の中に 存在していた。これはパ レート最適解の中で圧縮 比が最も大きい場合の翼 型である。また、高転向 角 と低損失を同時に実現するのは非常に難しい問題であることがパレート集合から明らかとなった。 また 、個々 の性 能項目 につ いて個 別に 単一目 的の最 適化 を行っ た場 合でも 、 SOGA の最 適解 と MOGAの対応する解は、ほぼ同じような性能値を持ち形状もよく似ている。これより、SOGAはMOGA のパレート解の集合の一要素を求めていることがあらためて確認できた。 今回の計算は、NWT 上で 64PE を用いて、１世代約 20 分、80 世代約 27 時間かかった。共同研究 を通じ、この研究の計算 時間を提供していただい た航空宇宙技術研究所に この場を借りて謝意を表 し ておきたい。なお、課金を払ってこの計算をするとすれば、1PE・１時間当たりの使用料が 2,000 円 であるから、約350万円の計算費用がかかる。しかし、すべての解が意味を持っているのであるから、 １パレート解あたりでは ５万４千円にすぎないと もいえる。さて、この値 段は、あなたにとって高 い だろうか、安いだろうか。 7-39 図7.24 MOGAの最適化過程 7-40 図7.25 スケーリングされた評価関数空間におけるパレート面 （３次元立体図と各射影図） 7-41 図7.26 パレート最適な翼型形状の変化 7-42 図7.27 CDAとMOGA翼型の形状及び圧力分布の比較 表7.1 CDAとパレート最適な翼型との性能比較 7-43 参考文献 [1] 芝木 秀 俊, 福 島 雅夫 , FORTRAN77 最 適化プロ グラミ ング, 岩 波コン ピュータ ーサイ エンス , 岩波書店, (1991). [2] F. Liu and A. Jameson, “Multigrid Navier-Stokes Calculations for Three-Dimensional Cascades”, AIAA Journal, Vol. 31, No. 10, pp. 1785-1791, (1993). [3] A. Dadone and P. De Palma, “An Adaptive-Multigrid Upwind Solver for Two- and ThreeDimensional Viscous Flows”, Proceedings of the Third ECCOMAS Computational Fluid Dynamics Conference, Paris, pp. 986-992, (1996). [4] G. N. Vanderplaats, “Approximation Concepts for Numerical Airfoil Optimization”, NASA TP-1370, (1979). [5] G. N. Vanderplaats, Numerical Optimization Techniques for Engineering Design: with Applications, McGraw-Hill, Inc., New York, (1984). [6] W. H. Press, et al., Numerical Recipes in FORTRAN : the art of scientific computing, 2nd ed., Cambridge University Press, Cambridge, (1992). [7] S. Kirkpatric, “Optimization by Simulated-Annealing: Quantitative Studies”, Journal of Statistical Physics, Vol. 34, No. 5/6, pp. 975-986, (1984). [8] D. E. Goldberg, Genetic Algorithms in Search, Optimization & Machine Learning, AddisonWesley Publishing Company, Inc., Reading, (1989). [9] O. Baysal and M. E. Eleshaly, “Aerodynamic Design Optimization Using Sensitivity Analysis and Computational Fluid Dynamics”, AIAA Paper 91-3322-CP, (1991). [10] O. Baysal , M. E. Eleshaly and G. W. Burgreen, “Aerodynamic Shape Optimization Using Sensitivity Analysis on Third-Order Euler Equations”, AIAA Paper 91-1577-CP, (1991). [11] J. L. ケルブロック著， 梶昭次郎訳 , ジェットエンジン概論, 東京大学出版会, 東京, (1991). [12] G. W. Burgreen and O. Baysel, “Aerodynamic Shape Optimization Using Preconditioned Conjugate Gradient Methods”, AIAA Paper 93-3322-CP, (1993). [13] M. Ogot, S. Aly, R. B. Pelz, F. Marconi and M. Silari, “Stochastic Versus Gradient-Based Optimizers for CFD Design”, AIAA Paper 96-0332, (1996). [14] D. Quagliarella and A. Della Cioppa, “Genetic Algorithm Applied to the Aerodynamic Design of Transonic Airfoils”, AIAA Paper 94-1896-CP, pp. 686-693, (1994). [15] K. Yamamoto and O. Inoue, “Applications of Genetic Algorithm to Aerodynamic Shape Optimization”, AIAA Paper 95-165 (1995). [16] I. De Falco, R. Del Balio, A. Della Cioppa and E. Trantino, “Breeder Genetic Algorithm for Airfoil Design Optimization”, Proceedings of the Forth IEEE International Conference on Evolutionary Computation, Nagoya, pp. 71-75, (1996). [17] 伊与田 真志, “遺伝的アルゴリズムによる翼形状の最適化”, 東北大学大学院修士論文, (1996). [18] J. Katz and A. Plotkin, LOW-SPEED AERODYNAMICS from Wing Theory to Panel Methods, McGraw-Hill, Inc., NewYork, (1991). [19] A. M. Kuethe and C-Y Chow, Foundations of Aerodynamics, John Wiley & Sons, New York, (1991). [20] C. A. J. Flecher, Computational Techniques for Fluid Dynamics 2, Springer-Verlag, New 7-44 York, (1991). [21] 三宅 裕, “流体工学におけるパネル法”, 日本機械学会誌, 87-785, pp. 322-327, (1987). [22] R. M. Hicks and P. A. Henne, “Wing Design by Numerical Optimization”, Journal of Aircraft, Vol. 15, pp. 407-412, (1978). [23] 日本機械学会編, 工学問題を解決する適応化・知能化・最適化, 技報堂出版, 東京, (1996). [24] 伊庭 斉志, 遺伝的アルゴリズムの基礎, オーム社, 東京, (1994). [25] G. N. Vanderplats, “ADS-A FORTRAN Program for Automated Design Synthesis, Version 3.00”, Engineering Design Optimization, Inc., Calfornia, (1988). [26] 石澤 和彦,“航空用ガスタービンエンジンの新技術”, 第34回飛行機シンポジウム講演集, 鳥取, pp. 269-272, (1996). [27] 梶 昭次郎, “圧縮機技術の現状と将来の展望”, 日本ガスタービン学会誌, Vol. 22, No. 88, pp. 4-11 (1995). [28] 千葉 薫, 小林 健児, “航空用大型エンジンの圧縮機”, 日本ガスタービン学会誌, Vol. 22, No. 88, pp. 12-19 (1995). [29] B. Lakshminarayana, Fluid Dynamics and Heat Transfer of Turbomachinery, John Wiley & Sons, New York, (1996). [30] N. Satofuka, et al., “Computation of Two-Dimensional Transonic Cascade Flow Using a New Navier-Stokes Solver ”, AIAA Paper 86-1381, (1986). [31] M. Giles, M. Drela and W. T. Thompkins, Jr., “Newton Solution of Direct and Inverse Transonic Euler Equations ”, AIAA Paper 85-1530, (1985). [32] J. D. Denton, “An Improved Time Marching Method for Turbomachinary Flow Calculation”, ASME Paper 82-GT-239, (1982). [33] C. Hah, “Navier-Stokes Calculation of Three-Dimensional Compressible Flow across a Cascade of Airfoils with an Implicit Relaxation Method”, AIAA Paper 86-0555, (1986). [34] W. N. Dawes, “A Numerical Analysis of the Three-Dimensional Viscous Flow in Transonic Compressor Rotor and Comparison with Experiment”, ASME Paper 86-GT-16, (1986). [35] J. N. Scott and W. L. Hankey, Jr., “Navier-Stokes Solutions of Unsteady Flow in a Compressor Rotor”, ASME Paper 86-GT-226, (1986). [36] T. H. Fransson and M. Pandolfi, “Numerical Investigation of Unsteady Subsonic Compressible Flows through an Oscillating Cascade”, ASME Paper 86-GT-304, (1986). [37] U. K. Singh, “A Design Time Marching Method for the Generation of Blade in a Casecade”, ASME Paper 86-GT-1767, (1986). [38] 西山 哲男,“翼型学”,日刊工業新聞社, 東京, (1992). [39] J. M. Sanz, “Design of Supercritical Cascades with High Solidity”, AIAA Journal, Vol. 21, No. 9, pp. 1289-1293, (1983). [40] J. M. Sanz, “Improved Design of Subcritical Cascades Using Complex Characteristics and Boundary-Layer Correction”, AIAA Journal, Vol. 22, No. 7, pp. 950-956, (1984). [41] J. M. Sanz, “Automated Design of Controlled Diffusion Blades”, ASME Paper 88-GT-139, (1988). [42] D. E. Hobbs and H. D. Weingold, “Development of Controlled Diffusion Airfoils for 7-45 Multistage Compressor Application”, Journal of Engineering for Gas Turbines and Power, Vol. 106, pp. 271-278, (1984). [43] H. Rechter, W. Steinert and K. Lehmann, “Comparison of Controlled Diffusion Airfoils with Conventional NACA65 Airfoils Developed for Stator Blade Application in Multistage Axial Compressor”, Journal of Engineering for Gas Turbines and Power, Vol. 107, pp. 494-498, (1985). [44] S. Obayashi and G. P. Guruswamy, “Convergence Acceleration of an Aeroelastic NavierStokes Solver”, AIAA Journal, Vol. 33, No. 6, pp. 1134-1141, (1995). [45] V. Venkatakrishnan, “Preconditioned Convergence Gradient Method for the Compressible Navier-Stokes Equations”, AIAA Journal, Vol. 29, No. 7, pp. 1092-1099, (1991). [46] S. Yoon and A. Jameson, “Lower-Upper Symmetric-Gauss-Seidel Method for the Euler and Navier-Stokes Equations”, AIAA Journal, Vol. 26, No. 9, pp. 1025-1026, (1988). [47] B. S. Baldwin and H. Lomax, “Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows”, AIAA Paper 78-257, (1978). [48] H. Tamaki, H. Kita and S. Kobayashi, “Muti-Objective Optimization by Genetic Algorithms: A Review”, Proceeding of the Forth IEEE International Conference on Evolutionary Computation, Nagoya, pp. 517-522, (1996). [49] 佐藤 義雄, 実習グラフィックス, アスキー出版局, 東京, (1988). [50] Y. Elzar and R. P. Shreeve, “Viscous Flows in a Controlled Diffusion Compressor Cascade with Increasing Incidence”, Journal of Turbomachinery, Vol. 112, pp. 256-266, (1990). [51] 中橋 和博, 荒川 忠一, “翼列まわりの流れの数値解析とその検証”, RC104 流れ解析プログラム 検証研究分科会 研究成果報告書, 日本機械学会編, pp. 129-197, (1994). [52] 日本物理学会編, スーパーコンピュータ, 培風館, 東京, (1985). [53] 松尾 裕一, “動静翼列を過ぎる非同期性非定常流れの並列数値解析”, 第9回数値流体シンポジ ウム講演論文集, 東京, pp. 423-424, (1995). [54] Z. Michalewicz, Genetic Algorithm + Data Structure = Evolution Programs, Second extended edition, Springer-Verlag, Berlin, (1994). [55] 坂和 正敏, 田中 雅博, 遺伝的アルゴリズム, 朝倉書店, 東京, (1995). [56] 奥野 正寛, ミクロ経済学, 日本経済新聞社, 東京, (1996). [57] H. Ishibuchi and T. Murata, “Multi-Objective Genetic Local Search Algorithm”, Proceedings of the Forth IEEE International Conference on Evolutionary Computation, Nagoya, pp. 119-124, (1996). [58] C. Poloni, G. Mosetti and S. Contessi, “Multi Objective Optimization by GAs: Application to System and Component Design”, Proceedings of the Third ECCOMAS Computational Fluid Dynamics Comference, Paris, pp. 258-264, (1996). [59] J. Horn and N. Nafpliotis, “Multiobjective Optimization Using the Niched Pareto Genetic Algorithm”, IlliGAL Report, No. 93005, (1993). [60] C. Fonseca and P. J. Fleming, “Genetic Algolithm for Multiobjective Optimization: Formulation, Discussion and Generation”, Proceedings of the 5th International Conference on Genetic Algorithms, Morgan Kaufmann Publishers, Inc., San Mateo, pp. 416-423, (1993). 7-46 [61] S. Tutsui and Y. Fujimoto, “Forking Genetic Algorithm with Blocking and Shrinking Modes (fGA)”, Proceeding of the 5t h International Conference on Genetic Algorithm, Morgan Kaufmann Publishers, Inc., San Mateo, pp. 206-213, (1993). [62] 山本 研二, “遺伝的アルゴリズムによる空力設計”, 東北大学修士論文, (1995). [63] W. W. Bathie, Fundamentals of Gas Turbines, 2nd Ed., John Wiley & Sons, New York, p. 258, (1996). [64] 北野 宏明, 遺伝的アルゴリズム2, 産業図書, 東京, (1995). GA に関する手頃な入門書としては文献[24]を薦める. また, GA の最近の動向を踏まえ基礎理論から 応用まで特に興味深いトピックを紹介している入門書として, 以下の本をあげておく. メラニ ー・ ミッ チェル 著, 伊庭 斉志監 訳, 遺伝 的アル ゴリ ズム の方法, 東 京電 気大学 出版 局, 東京 , (1997). 7-47