スーパーコンピュータを用いた音速域航空機の空力的研究
松島 紀佐， 山崎 渉， 中橋 和博
東北大学大学院工学研究科航空宇宙工学専攻

概要
本研究では音速域という複雑な現象下での空力最適設計を、形状選択手法として遺伝的アルゴリズムを、また、流れ解析手法として非構造格子流れ解析コードを用いておこなった。遺伝的アルゴリズムは、シミュレーションによるパラメトリックスタディを多数回（通常数百回から数千回程度）行うことが本質である。東北大学情報シナジーセンターのスーパーコンピ
ュータ SX-7 を利用し、計算を効率的に行うことで、衝撃波による抵抗を25％低減する形状を得ることが出来た。この形状最適化過程を解析することで衝撃波抵抗低減のメカニズムについて、いくつかの知見が得られた。

1. はじめに「より速く」というのは人々の航空機への要望の一つである。速い航空機としてまず頭に浮かぶのはコンコルドのような超音速機である。しかしながら、超音速機はソニックブームやオゾン層破壊といった環境への悪影響やエネルギー消費の低効率といった難しい問題を抱えている。このような環境やエネルギー問題の解決が必要な超音速機の代わりに、近年、未来型航空機として乗客 200 人クラスの音速機の研究が米国ボーイング社などで行われた[1]。遡って、50 年代や 70 年代にも音速付近で航行するための翼型や航空機の研究が行われたこともあった[2,3]。音速機とは巡航速度がマッ
ハ 0.95〜0.98 の旅客機であり、現在就航している遷音速旅客機より 15％程度高速となる。この高速化の飛行時間短縮により、機体の運用効率が増すため、航空会社にとっては航空機保有台数が少なくてすむというメリットがある。
音速付近の流れは衝撃波の挙動が複雑で，優れた空力性能を得るためには衝撃波の制御が重要となってくる。そのためには、実際的な形状で高精度の大規模流れ解析計算を行った上で、空力形状の最適化を行う必要がある。本研究ではこのような複雑な現象下での空力最適設計を、形状選択手法として遺伝的アルゴリズムを、また、流れ解析手法として非構造格子 Computational Fluid Dynamics(CFD)として国内外で広く知られている Tokoku University Aerodynamic Simulation（TAS）コード[4,5]を用いておこなう。両者を計算機上で有機的に連携させて空力最適化を試みる。遺伝的アルゴリズ
ムは、シミュレーションによるパラメトリックスタディを多数回（通常数百回から数千回程度）行うことが本質である。本研究のように、高精度の大規模流れ解析を必要とする最適化にはスーパーコンピュータ、並列計算が必要不可欠である。今回、東北大学情報シナジーセンターのスーパーコンピュータ SX-7 を利用させて頂き、効率的に短期間で、衝撃波による抵抗を低減する形状を得ることが出来た。以下、この形状最適化の手法や過程について、パラメトリックスタディのための効率の良い並列 JOB構成法を含めて述べる。また、基本(初期)形状と最適化後のそれぞれの形状周りの流れ場を可視化し比較することで、抵抗低減の物理的メカニズムについての解析を行う事が出来た。その具体的手法についても述べる。

2. 最適化2.1 基本形状
基本形状は航空宇宙技術研究所（現宇宙航空研究開発機構）で開発された超音速実験機 NEXST-1[6,7]を用いる。以前の研究により、音速域での Euler 計算の妥当性やNEXST-1 モデルの音速機初期形状としての有効性等を確認しているので[8,9]、それを踏まえた上で主翼形状の空力最適化を行う。Fig. 1 に基本形状(Baseline)及び、計算のための座標系とある平面内での空間格子を示す。加えて Fig. 1 には基本形状航空機がマッハ 0.98 で航行した際の流れ場の圧力分布も示されている。翼の後縁に衝撃波が停留しているのが分る。

YX

Z

Fig.1 NEXST-1 Geometry and Pressure Contour at M ∞ = 0 . 98 .
2.2 最適化問題定義
最適化は NEXST-1 主翼形状に対して行う。形状最適化問題の定義はまとめると、・設計変数：翼断面形状と捻り角の修正（計 66 変数）
・目的関数：M∞=0.98 での圧力抵抗係数CDPの最小化
・制約条件：揚力係数CL=0.26、翼平面形固定、主翼の内部体積一定ということになる。形状をコントロールする設計変数を Fig. 2 に示した。主翼は左右対称であるので、片側の半スパン形状で考える。15％位置までは胴体で覆われているので、それより外側の 6 つのスパン断面を選択する。翼面上の黒い線分で示されているのがその位置である。最も外側は 100％半スパン位置である Tip となっている。選択したスパン位置の翼断面形状 (Fig.2 の一点鎖線で抽出した形状)を 10 個の自由設計

変数で変化させる。このとき翼型前縁と後縁は固定である。更に捻り(Twist)角を変化させるための設計変数を 1 個定義している。つまり、各翼断面に 11 個の設計変数があり、断面数が 6 つであるので、設計変数は 66 となる。選択したスパン断面以外の場所の断面形状はスプライン関数補間で決定している。
10 Airfoil ShapeVariables

Twist Variables

100% (Tip)
6 Design stations
15% (Fixed)
Fig.2 Design Parameters.

2.3 最適化手法
音速域における空力最適化では、強い非線形性が想定されるため、勾配計算を必要としない遺伝的アルゴリズムを最適化手法として採用した[10]。遺伝的アルゴリズム
では設計変数を任意に変化させた候補をいくつか(通常数十)つくり、その個体の性能(目的関数の値)をシミュレーションで決定する。目的関数の値のより良い個体を選び、それらをブレンドして次の世代の候補を作成し、再び各個体に対しシミュレーションを行う。このような処理を反復することで性能の悪いものが淘汰され最終的に最適解が残るわけである。この手法の特徴は、シミュレーション回数が個体候補数(数十)×世代数(数十)となり、計算負荷が大きい事である。特に、本研究における流体計算の場合は 1 個体のシミュレーションに、SX-7 スーパーコンピュータを利用して、約 1 時間かかる。これは、非粘性の Euler 方程式を解く計算の場合である。粘性の考慮を入れた Navier-Stokes(NS)計算では境界層の解像度を保証する細かい格子を用いなければならないため約 20 時間かかる。実際の流体挙動を記述する NS 計算で各個体の性能評価を行うのが、理想的ではあるが、現状の計算機資源では、現実的ではない。本研究では、境界層の効果が最適化に与える影響は少ないと考えられるので[8,9]、最適化過程における性能評価シミュレーションには Euler 計算を用いる。ただし、最終的に最適化された形状に対し NS 計算を行い、最適化の信頼性を確認する。

2.4 遺伝的アルゴリズムによる最適化計算(計算機利用法)
遺伝的アルゴリズムにおける進化計算では候補集団の個体数を 30 とし 70 世代の計算（計 2100 回の流れ場のシミュレーション）を行った。その結果、後述するように25％に達する大幅な空力性能の向上が図られた。この進化計算を経済的に進めるため、各世代における各個体の性能評価シミュレーションを並列に行い計算機利用の効率化を図った。簡単な考え方ではあるが、その効率的利用法について、本節で述べる。
本最適化の過程で行った並列計算は、実行されるプログラムを並列化するのでなく、単体の CPU 仕様のプログラムを並列に複数の CPU で稼働させるものである。計算条
件を変化させてデータを収集するパラメトリックスタディや本研究のような進化計算における各世代の個体評価を効率的に且つ自動的に行うための方法である。Fig.3(次頁)に概念図を示す。
単体仕様のプログラム名を”SIMULATION”とする。並列実行を実現するにはまず、”SIMULATION”のプログラム実行を駆動するマスター(制御)部分と、そのマスターから駆動されるスレーブ的な”SIMULATION”実行部分からなる JOB 投入のマクロ的な枠組みを作成する。各 CPU で稼働する”SIMULATION”にはそれぞれ異なった識別子を持った、入力・出力ファイル、さらには計算条件などのパラメターファイルを割り当てなければならない。”SIMULATION”プログラムに、稼働させるシミュレーショ
ン条件のあった適切な各種ファイル名を設定するための処理を書き加える必要がある
が、これらの処理は簡単に行える。
例えば 10 並列で行った場合、進化計算の 1 世代分の個体性能評価 Euler 計算が 3 時間で終了する。これを単体 CPU で処理した場合は 30 時間かかる。実際問題として待ち時間なども考慮すると、並列計算によるスループットの効率向上は多大である。

2.5 最適化の過程
Fig.4 に最適化の開始から終了まで、最適化の目的関数である抗力係数値の変化の履歴を示す。各世代 30 個体のもつ抗力係数値の平均、最大、最小値をプロットしている。
初期 Baseline 形状の抗力係数値を太い実線で示している。世代を重ねるごとに目的関数の最小化が進行している様子がわかる。70 世代付近で最小値の変化がなくなり、平均値と最小値が同等になってきたので、70 世代で最適解に収束したと見なした。
本最適化においては可視化ツール等も積極的に活用し、最適化プロセスの状況の観測を行なった。例えば、次節で解析する Shock Function という関数を基にした衝撃波可視化図（Fig. 7 参照）などである。最適化の進行に伴い、外翼の衝撃波が消失していく様子が確認できた。

単一 CPU
入力
ファイル

計算条件
パラメーターPROGRAMSIMULATION

出力
ファイル

複数 CPU

並列実行
MASTER Controller

Drive

入力
ファイル１

計算条件１
パラメーターSLAVESIMULATIONCPU1

MASTER

出力
ファイル１

Drive
SIMULATIONPROGRAM

Drive

入力
ファイル n

計算条件 nパラメーター

SLAVE

SIMULATION・・・・・・・・

CPU n

Fig.3 Outline of Parallel Simulation Job Processing

出力ファイル n

0 .0 2 6A verage

M in im u m0 .0 2 4M a x im u m

I n it ia l

CDP0 .0 2 20 .0 2 00 .0 1 80 .0 1 60 .0 1 4010Fig.4

203040G e n e ra tio n506070

History of Optimization Process.

3. 最適化結果
3.1 空力性能
Table. 1 に最適化による空力性能の向上についてまとめた。初期値である Baseline 形状における空力性能値と比較した。また、その比較は、Euler 計算と NS 計算の両方で行った。揚力係数（CL）一定の下で抗力が 37 Counts(0.0037)減少した。
Table. 1 The Aerodynamic Performance of Initial and Optimal Airplanes
Wing volume
Beseline

1.000

Optimal1.003

Evaluation
Euler
NS

A.O.A.
3.48°
3.46°

CL
0.26
0.26

Euler

4.02°

0.26

NS

4.24°

0.26

CDP
0.01855
0.019650.01481-37 counts
0.01596-37 counts

CDf0.01078

L/D
14.02
8.56

-

17.56

0.01093

9.68+13.1%

3.2 抵抗低減についての考察
今回の最適化を通して、低抵抗の翼設計の指針を得るため、抵抗が低減したメカニ
ズムのついての考察を行った。まず、形状の変化について、Baseline と最適化後の形状を比較する。Fig.5 は、主翼上面の等高線を描いたものである。最適化により翼面上の凹凸が顕著になっている。このように凹凸をつけることで圧力波の圧縮と膨張が効果的におこるようにし、3 次元的に衝撃波を緩和することで造波抵抗を減少させている。外翼部の 2 つのディップにより衝撃波が等エントロピー的な圧縮波に弱められていることが、Fig .7 の衝撃波関数の可視化図を合わせて見ることにより分る。

Fig. 6 より、最適形状では翼端付近で捻り角が負になり､より捻り下げの効いた形状に進化している事がわかった。また、この捻り下げの効果で、外翼部分では、上面のクレスト位置が後方に移動することでの推力増加が顕著であった。
Fig. 5 からは、正のキャンバーが効いた形状に進化している事がわかる。つまり、キャンバーを増大させることで揚力を得、そのトレードオフで、実質的な翼の迎角が減少するように捻り角分布が変化しているといえる(Fig. 6 参照)。この迎角の減少も抗力低減に寄与している。

Dip
OptimalBaseline

Fig.5 Contour Maps of the Height of the Wing Upper Surface
3Twist Angle[deg]210-1-2-3Initial

Optimal-4-50204060
Semi-Span[%]80100

Fig.6 Twist Angle Distribution in the Span-wise Direction.

次に現象の観点から衝撃波の挙動を調べ、抵抗低減のメカニズムを考えてみる。
Baseline 形状と最適化形状で発生する衝撃波を(r関数 Fshock = V ⋅ ∇P) (a ⋅ ∇P ) を用いて数値化し[11]可視化する。

ここで V は速度ベクトル、∇P は圧力の勾配、a は音速である。
F shock = 1 なる等値面を可視化したものが Fig.7 である。最適化形状においては、設計により外翼の衝撃波が消失している様子が確認できる。
Fig.8 は 40％セミスパン断面内での関数 F shock の等高線分布である。圧縮波となっている部分のみを描いている。最適化された翼面からは弱い衝撃波が広範囲に渡って生成されている様子がわかる。つまり、衝撃波の上流領域で等エントロピー的な圧縮を起こすことで、発生する衝撃波を弱め、抵抗低減に至るという事がいえる。

Baseline

OptimalFig. 7 Visualization of Shock Locations Using a Shock Function.

Fig.8 Visualization of Compression Region at 40% Semi-span Station.
3.3 最適化結果の有効性検討
今回の最適化は、マッハ数＝0.98 の単一設計点で行った。最適化の有効性を確認する目的で、マッハ数に関するオフデザイン特性と粘性効果の影響に関する調査を行った。
0.028
InitialOptimal
0.026

CDP

0.024
0.022
0.020
0.018
0.016
0.014
0.7

0.9

1.1

1.3

1.5

Mach number

Fig.9 Drag Bucket at C L = 0.26 .
-0 .80 .0 0
In itia lO p tim a l-0 .4-0 .1 00 .0

Cp(y/c)-0 .0 5-0 .1 50 .4-0 .2 0-0 .2 5-0 .10 .30 .7

1 .10 .8
(x /c )

Fig.10 Cp Distribution at 40% Semi-Span Station by Euler Simulation
0 .0 0-0 .8
In itia l

O p tim a l-0 .4-0 .1 00 .0Cp

(y/c)

-0 .0 5-0 .1 50 .4-0 .2 00 .8-0 .2 5-0 .1

0 .30 .7

1 .1

(x /c)

Fig.11 Cp Distribution at 40% Semi-span Station by NS Simulation.

まず、オフデザイン特性について Fig. 9 に示す。揚力係数 0.26 としてマッハ数を変化させて抗力係数をプロットしたものである。本研究の最適化により遷音速域(マッハ数0.7〜1.2)全般にわたって抗力が低減されたことが確認された。
次に、粘性の影響が局所的で線形的である事を示す。Fig. 10 と Fig. 11 に、初期形状と最適化形状が実現する表面 Cp 分布を比較したグラフを示す。濃い線が最適化された形状に関するもので、薄い線が初期形状である。これらのグラフは、40％半スパン位置でのものであるが全てのスパン断面で傾向を代表するものである。Fig.10 が Euler計算による比較、Fig.11 が NS 計算による比較である。最適化による圧力分布変化の傾向は両計算結果でよく一致している。破線で囲んだ垂直(に近い)衝撃波に NS 計算とEuler 計算の定量的な差が見られるのみである。NS 計算結果は境界層の影響で衝撃波が弱められている。しかし、最適化による相対的な変化量は Table.1 の結果からもわかるように、Euler と NS の両者で差がないと考えられる。以上より、Euler 計算を用いた本手法による最適設計の総合的な有効性が確認されたといえる。

4． 結論
本研究では NEXST-1 モデルを初期(Baseline)形状として、66 個の設計変数を用いた大きな設計空間内で、主翼形状の空力最適化を音速域で行った。最適化手法は遺伝的アルゴリズムである。また、空力性能の評価は非構造格子 CFD コード“TAS コード”を用いて行った。特に主翼の断面形状に注目して最適化を行った。マッハ数 0.98 とし、揚力一定の拘束条件の下、L/D の向上を目的とした最適化を施した。70 世代後の最適解を得ることが出来、L/D の 25％程度の改善が可能となった。この最適化を通して、以下の成果が合わせて得られた。
計算機ツールに関して、
１．東北大学情報シナジーセンターのスーパーコンピュータ SX-7 を利用して、パラメトリックスタディシミュレーションを効率的・経済的に実行する枠組を作成した。
２．空力形状最適化による形状変化が流れ場の物理現象とどう連携しているのかを明らかにするための解析方法を提案した。その一例として、波動(衝撃波)抵抗低減
過程の監視のため衝撃波関数を用いた可視化による解析法を提案した。
高遷音速域における空力性能の優れた翼形状の一般的条件の手がかりとして、
１．翼断面形状のキャンバーを大きくし捻り角を減少させる事で下面(高圧面)クレスト位置を前方に寄せて下面の推力を増加させる。
２．同様に、捻り角を減少させる事で上面クレスト位置を後方に移動させ、上(低圧面での推力を増加させる。衝撃波が後縁付近にある場合に有効。
３．衝撃波上流で等エントロピー的な圧縮がなされる様に形状を設計する事により衝撃波の強さを弱めて造波抵抗の減少を図るべきである。
といった知見を得た。
また、このような形状変化が非設計点での性能を向上させることも確認した。

謝辞
本研究は東北大学情報シナジーセンターとの共同研究として行われました。スーパーコンピュータ SX-7 上でのプログラム最適化に関し、センターより適切な助言を頂いた事を感謝致します。

References
1)
2)
3)
4)
5)
6)
7)
8)
9)

10)
11)

Boeing Company web page http://www.boeing.com
Becker, V. J.: The High Speed Frontier, NASA SP-445, pp. 36-59, pp. 84-117, 1980.
Haines, B.: The Aerodynamic Design of the Wing–Fuselage for a Near-Sonic Transport– Early Thoughts and Results Revisited, AIAA Paper 2002-0517, 2002.
Ito, Y. and Nakahashi, K.: Surface Triangulation for Polygonal Models Based on CAD
Data, International Journal for Numerical Methods in Fluids, Vol.39, Issue 1,
pp.75-96, 2002.
Sharov, D. and Nakahashi, K.: Reordering of Hybrid Unstructured Grids forLower-Upper Symmetric Gauss-Seidel Computations, AIAA Journal, Vol.36, No.3,
pp.484-486, 1998.
Sakata, K.: Supersonic Experimental Airplane (NEXST) for Next Generation SST
Technology, AIAA Paper 2002-0527, 2002.
Matsushima, K., Iwamiya, T. and Ishikawa, H.: Supersonic Inverse Design of Wings
for the Full Configuration of Japanese SST, Proceedings of the 22nd ICAS Congress,
ICAS 2000 CONGRESS, pp.213.1-213.8, 2000.
Matsushima, K., Yamazaki, W. and Nakahashi, K.: Transonic Design of SST,
Proceedings of IUTAM Symposium TRANSSONICUM IV, Germany, 2002., Kluwer
Academic Publishers, pp. 317-324, 2003.
Wataru Yamazaki, Kisa Matsushima and Kazuhiro Nakahashi‘Aerodynamic Optimization of NEXST-1 SST Model at Near-Sonic Regime’
42nd AIAA Aerospace Sciences Meeting and Exhibit, AIAA-2004-34, Reno, Nevada,January, 2004.
ESTECO Company web page http://www.esteco.it/
Paparone, L. and Tognaccini, R.: Computational Fluid Dynamics-Based Drag
Prediction and Decomposition, AIAA Journal, Vol.41, No.9, pp.1647-1657, 2003.