CFD におけるペタフロップス級計算機への期待

中橋和博 1
1. はじめにマイクル・クライトンは小説“Time Line”の書き出しで、「19 世紀末の人々は飛行機の出現は予想できただろうが、数万機の飛行機が常時空中を飛び交っているという 100 年後の世界はとうてい想像できなかったであろう」と記している。1903 年にライト兄弟が空への旅立ちの扉を開けて以来、様々な飛行機が空へと羽ばき、今日では音の速さ近くで飛行するジェッ ト 旅 客 機 が 人 々 の 重 要 な 交 通 手 段 と な っ て い る 。 そ し て 2004 年 ９ 月 に は 米 国 のSpaceShipOne（図 1）が大気圏を越えて宇宙の入り口まで飛翔した。過去 100 年の飛行機の発展は、人類の長い歴史のなかでも特筆すべきものの一つであろう。

図 1. SpaceShipOne (//www.scaled.com/)

飛行機は、空気のお陰でその重い機体を空中に浮かべる。ところが同時にその空気のせいで発生する抵抗と闘ってもおり、その闘いこそが飛行機の発展の歴史でもある。今日のジェット旅客機も技術的に熟成したかに見えるが、エアライン間の厳しいコスト競争や環境への配慮から性能改善のための努力は未だに衰えをみせない。2005 年春に初飛行したエアバ ス
A380 や、開発中のボーイング B787、エアバス A350 等は既存の旅客機よりも２割の燃費改善をセールスポイントにしていることからも分かろう。
航空機開発においての空力設計にはライト兄弟の時代から風洞試験が主役を演じてきた。
しかし近年の数値流体力学（Computational Fluid Dynamics; CFD）と高速計算機の進歩により，計算科学の重要性が急激に大きくなっている。図 2 は，主要な航空機の初飛行と CFD および計算機の発展とを並べて書いたものである。筆者の研究室では、図 2 中の unstructuredmesh(非構造格子)を用いた CFD コードを過去十年程をかけて開発してきた。飛行機の全体形状周りの流れの詳細計算は、以前は何ヶ月もの事前作業と膨大な計算時間が必要であったが、非構造格子の採用により格子生成に要する時間が大幅に短縮された。並列計算機による時間
1

東北大学大学院工学研究科航空宇宙工学専攻
1

図２．航空機開発と数値流体力学

CaseMemory# of nodes

Executing Time# ofiterations

Machine1PE

Parallel35

15 minutes(128 PE’s)

ONERA M6(Euler)

500M

234,621

1,000

SGI Origin
2000

Regional Jet(NS)4GB

2.2 million

3,000

NEC SX-4

3.3

2.5 hours(32 PE’s)

ONERA M5(NS)10GB

6.7 million

10,000NEC SX-7

1 week16 hours(16 PE’s)

Racing car(NS)

3.8GB

1.5 million

10,000NEC SX-7

3 days

10 hours(8 PE’s)

図３．並列計算機による計算時間短縮の効果2短 縮 [1] の 効 果 も 大 き い （ 図 3 ）。 現 在 、 こ の TAS コ ー ド （ TAS =Tohoku university
Aerodynamics Simulation）[2]は東北大学内の研究グループだけでなく、国内外のメーカーや研究所でも航空機開発等に用いられている。
さて、成熟したかに見える CFD ではあるが、まだ多くの課題を抱えている。２０年ほど前、
CFD が発展すると風洞がいらなくなるという意見を述べた人がいて風洞運転担当者を不 安がらせたことがあるが、風洞試験は未だに航空機設計に不可欠である。無くなるどころか更に大型の風洞設備が作られつつある。
ここ では、航 空に 関連 した現 在の CFD が抱える課題について考え、それを解決するための筆者なりのアプローチを近未来のペタフロップス級のスーパーコンピュータへの期待をも込めて書いてみたい。
2. CF D の信頼性
今日の CFD は、巡航状態にある飛行機の揚力係数を非常に精度良く予測する。揚力係数より一桁小さな抵抗係数でさえも十分な格子点数を用いれば１％以内の精度で予測できるまでに至ったのは CFD を研究している者にとっても驚きである。そのため、CFD を用いた空力解析は既に航空機の性能解析には不可欠なツールとなり、それを用いた最適設計も近年とくに発 達が著しい。図 1 の SpaceShipOne も CFD のみで空力設計が行われたとも聞く。
しかし、高い予測精度を得るには格子をどの程度まで細かくすればよいか、どの乱流モデルが適当か、あるいは境界層の乱流遷移をどのよ うにモデル化するか等のノウハウ的な要素に 未だに頼っているのも現状の CFD であり、常に風洞試験による計算結果の確認が求められる。図 4 に示すように、翼の発生できる最大揚力の予測（失速点予測）も未だに精度良くできない。離陸や着陸では最大揚力近くで飛行するために、失速点予測の信頼性欠如は致命的であるものの、これに関しては２０年前 からの進展は小さい。このこともＣＦＤが風洞試験に 取って代われない理由である。

Exp .

Ｃｌ vs. α

図 4. NACA663-018 翼型の迎え角変化に対する揚力係数3空気の運動は Navier-Stokes 方程式で記述され、この方程式を厳密に解けばʻ流れが見えるʼはずである。しかしながら、航空機の翼表面流れのレイノルズ数は 10 の 7 乗のオーダーであり、極めて細かなスケールの渦が不規則に運動する乱流境界層を伴う。この微細な非定常運動までを現在の計算機で解像することは、計算時間と記憶容量の制限からほぼ不可能である。それゆえ、現在一般的に用いられている手法は、Navier-Stokes 方程式を微小時間で平均化して不規則変動を除外した Reynolds-Averaged Navier-Stokes(RANS)方程式を解くものである。この際、時間平均操作によって乱流応力と呼ばれる未知数が現れるため、これを「乱流モデル」を用いて補っている。このアプローチにより、巡航時の航空機の揚力や抵抗は前述のように高い精度で予測できるまでになった。しかしながら、この乱流応力は流れの条件によって大きく変わり、特に大きく剥がれた流れになると乱流モデルの定量的信頼性は著しく落ちる。乱流モデルの不完全性が計算の信頼性の限界を作っており、それが解決しない限りＣＦＤの適用範囲も限られたままであり、風洞の必要性も無くなくらなない。
現在、時間平均をせず、少なくとも格子で解像できる大きな渦は計算し、格子サイズ以下の渦運動は空間平均するという Large-Eddy Simulation (LES)の研究が盛んになっている。空間平均に伴うモデル化が必要であるものの、RANS よりも物理モデルの影響を減らそうとの努力である。そして、さらには物理モデルを用いない Direct Numerical Simulatio n (DNS)も乱流現象の理解のための研究に適用されつつある。しかし、飛行機全体周りの流れにこの DNSを適用するには、現在の計算機の少なくとも数桁倍もの演算能力を必要とし、まだかなり遠い未来であろうと言われてきた。とはいうものの、そろそろ物理モデルに頼るアプローチから抜けなくては、いつまで経っても CFD は風洞実験に頼らざるを得ないことになる。
3. 近未来の計算機と CFD
スーパーコンピュータの Top500 という話題が、最近では新聞記事にも取り上げられるようになった。連立一次方程式の解法プログラムである LINPACK というベンチマークソフトによってシステムの速度計測がなされ、上位の 500 機をリストにしたものであるが、昨年まで首位だっ た地球シミ ュレータが 今年は４位 まで下がり 、現在のト ップは Blue Gene/L で65,536 個のプロセッサを用いて 136.8TFLOPS を達成したとのことである。５年以内にはペタフロップスを超すスーパーコンピュータも現れるという。
計算機の性能は１８ヶ月で２倍向上すると言われているが、過去２０年の発達を見るとそれ以上の成長を遂げている。今後も今のペースを当面は保ちそうであり、計算機のユーザー側としては頼もしい。かつてスーパーコンピュータの最大のユーザーは航空 CFD であったが 、今日では必ずしもそうでは無くなっている。しかし、流れの数値計算は前述のように未だに未完成であり、新たなブレークスルーが求められている。ペタフロップス級のコンピュータを期待して CFD の次のステップを考えるべき時期であろう。
CFD は過去３０年程で急速な発展をした比較的新しい学問である。しかし、計算機の発達はそれ以上に展開が速く、CFD のアルゴリズム研究もしばしば計算機アーキテクチャーの進歩を追っかけることにもなった。ベクトル演算器を念頭においたアルゴリズム、メモリーの4急速な低価格化、キャッシュメモリーの利用、そして並列計算等であり、計算プログラム
そ のものは過去の遺産を引きずりつつ改良を加えてきている。と こ ろ が 、 最 近 の 計 算 機 は 益 々 多 数 の CPU を 用 い る 方 向 に 加 速 し つ つ あ り 、 更 に Grid
computing は CFD コードに大幅な変更を要求しつつある。これまでの蓄積を保持しながらも計算方法の枠組みを根本的に見直すべき時期に来ていることを示唆しているように思える。
また先に述べたように、CFD の信頼性を高めるには物理モデルへの依存性をできるだけ小さくすることが重要であるが、そのためには計算格子による解 像度を飛躍的に高める必要がある 。計算機の発達は大規模計算を可能にするものの結果も膨大となる。その膨大なデータから有用な物理情報を如何に効率よく取り出すかも今後の大規模計算に向けて考えなくならない課題である。
4. 大規模 CFD のための考察
今後の計算機の更なる発展を見込んで、大規模計算のための アプローチに要求される事項を 挙げてみた。
１． 複雑な形状を容易に取り扱えること。
２． 格子サイズが境界層や後流、衝撃波等の様々な流れの局所特性長に容易に適合できること。
３． 大規模な並列計算機に容易に対応できること。
４． 空間高次精度解法を容易に適用できること。
５． 大規模データの後処理が容易であること。
６． 計算アルゴリズムなり計算コードの構造が単純であること。
上記の１〜３は明らかであろう。４の高次空間精度も、今後の CFD の一つのターゲットである風切り音の発生源解析には強く求められる。５の後処理も、近未来の数十億、あるいは数兆の格子点を用いた計算では最も重要な課題の一つである。また、これまで以上に多彩な計算機アーキテクチャー、計算環境で使うことになることが予想され、計算コードの保守や改良、ポータビリティーのためには計算アルゴリズムなり計算コードの単純さは重要になる。
また、特に大学では学生が常時入れ替わるため、年輩の学生から下級生への計算コードの引き継ぎにおいても ６は重要な項目である。
現 在の CFD のアプローチを計算格子でもって区別すると、図 5 に示すように(a) Structuredgrid （ 構 造 格 子 、 あ る い は 境 界 適 合 格 子 、 BFC ＝ Boundary-Fitted Coordinate Grid ） , (b)Unstructured grid（非構造格子）, (c) Cartesian grid（直交格子）、に分けられる。航空の CFD
では、壁面摩擦ま で精度良く計算する必要性から物体面に格子線が沿った格子を用いる構造
格子 のアプローチが主流であった。しかしながら 、計算機の発達により飛行機全体まわりの３次 元計算が頻繁に行われるようになると構造格子の生成の困難さが問題となり、その解決手段とし て(b)の非 構造格子が最近では主流になりつつある。直交格子は、構造格子が考案される 以前の CFD 初期のアプローチであるが、その格子生成 の容易さから最近になり見直されつつ ある。

5(a) Structured grid(b) Unstructured grid(c) Cartesian grid

図５. 流体計算に用いられる計算格子
では、先に挙げた要求事項を満たすに図 5 のどのアプローチを採るべきであろうか。構造
格子のアプローチは１の複雑形状の取り扱いで難点がある。非構造格子は複雑な形状に対して優位であり、また２の格子サイズを局所的に細かくできる利点も持つ。しかしながら、非構造格子では格子の形状融通性のために空間精度はせいぜい格子幅の２次のオーダーである。現在、より高次な計算方法の研究が進められているが、計算時間やメモリーの大幅な増加を伴う。また、格子の融通性を認めているために、ときには局所的に歪んだ格子も扱わざるを得なくな るが、そのために局所的に空間精度が悪くなることがある。非構造なデータ構造 であるために、５の大規模データの後処理も困難が予想される。
直交格子法は、計算空間を x, y, z 方向に等間隔で分割するため単純ではあるが、局所的に格 子 サ イ ズ を 小 さ く し よ う と す る と 様 々 な 工 夫 が 必 要 で あ る 。 AMR(Adaptive MeshRefinement)と呼ばれる方法は部分的に直交格子を細分化するものであるが、凝った細分化手法の導入は直交格子の最大の利点である単純さを失うこととになる。また、航空 CFD のように翼表面で高い精度が要求される場合、物体上での境界条件の取り扱いも未だに大きな研究課題である。表面近くに別の格子を用いたり、あるいは Navier-Stokes 方程式側で修正するなどの方法が提案されているが、いずれの方法も複雑さをもたらす。これらアルゴリズムの複雑化は、並列計算や空間高次精度、さらには後処理をも難しくすることとなる。
しかし、直交格子は基本的 に単純であり、上記要求事項のうちの ２、そして物体での境界条 件を複雑化を招かない形で解決することができれば、今後の高性能計算機での超大規模計算には最も魅力的なものである。そのため、ここでは直交格子を基本としたアプローチを考えてみた。
5. Building-Cube Method
5.1 基本コンセプト
直交格子の利点は単純さであり、また格子を細かくしたり粗くすることが容易であることにある。ただし、格子細分化を局所的に不規則に用いるとデータ構造やソルバーの複雑化を招くため、決められたルールのもとで格子細分化を行うこととする。 ここでは、計算空間に様々なサイズの立方体ブロック（cube）を積み上げ(図 6)、個々の cube 内には高密度な等間隔直交格子を用いる計算法を考える。さまざまなサイズの cube を用意することで流れの局所特性長さに適合した格子密度を可能にし、また cube 内に等間隔直交格子を用いることで高精6

度化を可能にするとともに、アルゴリズムの単純さを確保する。また、空間を cube で分割することは、大規模並列化や大規模データの後処理でも有利である。Cube を積み木状に積み上げることから Building-Cube method（直交格子積み上げ法）と名付けた[3, 4]。

図６. 4 要素翼型に対する BCM 格子. (a) Cube 枠（赤線）. (b) 物体境界ちかくの直交格子拡大図（4 本に 1 本だけ描画）.

図７. 物体境界の表現（赤線は実の境界）

直交格子のもう一つの課題である物体境界条件も、 アルゴリズムやデータ構造の単純さを保つために図 7 に示すように格子セルのオン・オフで表現する。ただし、曲面境界を精度よく表現するには階段幅を十分に小さくしなければならない。DNS は境界層内の最小渦までをも解像しようとするもので、そのような計算を目指そうとすると自ずと壁近くで細かな格子が必要となり、そのような格子幅では階段状表現も実際の壁表面粗さに近くなり問題ないだ
ろうとの考えである。
Cube 形状の小領域を積み上げるこの方法では以下の利点が挙げられる。

7

１． 直交格子を基本に用いるため複雑な形状に対する格子生成が容易である。
２． Cube サイズを変更することで複雑さを伴わずに格子解像度を部分的に容易に変更できる。
３． 計算領域を Cube 状の小領域に分割することで大規模な並列計算に直接拡張が可能である。
４． 計算コードは、Cube 内の計算を行うサブルーチンと、cube 間の情報交換や計算順序などを扱う部分とに分けられ、単純なプログラム構造にできる。
５． Cube 内で等間隔直交格子を用いることで空間高次精度の計算法を採用できる。
６． 階段状境界条件の採用で格子データそのものが単純となり、更にデータ圧縮を用いることで非常にコンパクトにできる。
７． 直交格子上で計算結果を持つため後処理が容易あり、かつ画像圧縮技術の応用によるデータ圧縮が可能である。
８． 直交格子のために格子点の間引きが可能であり、それを利用した計算時間短縮が可能である。
5.2 NACA0012 翼型まわりの計算
NACA0012 翼型周りの流れの計算に応用した結果を図 8 に示す。計算条件はレイノルズ数
5x10 4 、 一様流マッハ数 0.5、迎え角３度である。Cube数は 522、cube内には 64x64 の直交格子を用いた。結果として約 200 万の格子セルを用いた計算である。また、壁近くの最小格子幅は翼弦長基準で 4.58x10 -4 である。計算は空間４次、時間２次精度で行い、乱流モデルは用いてない。
この比較的低いレイノ ルズ数では翼面上の境界層は層流を保つ。上面では境界層が途中から 層流剥離をするが、これを密度分布（図 8）で見ると剥離に伴い渦が周期 的に発生している様子が観察される。さらに興味深いのは翼後縁での現象である。剥離によって生じた渦が後縁に達した際、渦中心の低圧と翼下面の高圧との大きな圧力差により、翼下面から後縁を回 り込んで上方へと流れが生じ、それが上面の渦とは反対の回転を有する渦を生成する。また、この際に圧力波が発生して音速で遠方へと伝播していく様子が観察されよう。
翼後縁での風切り音は航空機の着陸時や発電用の風車で問題となっているものの、その発生メカニズムは完全には明らかになっていない。今回の計算は２次元であり、一方、翼の後縁 近くや下流では流れは３次元性を強めるために本計算結果だけでの物理現象を議論するのは危険であるが、細かで等方的格子を用いることで詳細な流れの観察が可能であることを示し ている。
5.3 ４ 要素翼型
NASA supercritical four-element airfoilを、レイノルズ数=2.83x10 6 、マッハ数 0.201、迎え角8.16° の条件で計算した。図 6 にBCM格子を示す。 Cubeの数は 481、物体近くの最小cubeのサ
イ ズは翼弦長を１として 2.9x10 -2 、cube内には 256x256 の等間隔直交格子を用いており、8従って壁近くの最小格子幅は 1.4x10 -4 である。総格子セル数は約３０００万になる。図 6 の下の図は前縁スラットの後縁近くの拡大図である。図では翼面が階段状で表現されている様子が見える。しなしながら、この図は格子を４本毎にのみ描いたものであり、実際の計算に用いた格子では物体面の階段幅は図 6 の 4 分の 1 となっている。計算は空間４次精度、時間２次精度で乱流モデル等を用いずに行った。
図 9 に瞬間マッハ数分布の図を示す。前縁スラットの上面では流れが急激に膨張し、一様流マッハ数が 0.201 であるにもかかわらずスラット上面の局所的な流れはマッハ 0.8 近くまで加速されている。また、最も膨張した直後あたりから境界層が乱れ、その渦が下流へと放出されている様子が再現されている。主要素翼の上面でも同様に前縁の少し下流から境界層9が乱れている。この様子を瞬間密度分布の図（図 10）で見ると、細かな渦が生成されて下流へと転がっている様子が見て取れる。従来の計算では、壁に沿う方向には引き延ばされた格子を用いるため、この様な細かな渦は捕らえられなかったが、本計算では縦横に同じ幅の格子を用いているため、細かな渦の生成と下流への流れ下りの様子が再現されていると言える。
翼面を階段状で表現しているために細かな渦が生成されたとの疑いもでるが、その場合は上面だけでなく下面でも乱れが生じるはずである。しかし、下面では境界層の乱れは主要素翼の下流凹部剥離域以外では観察できない。このことから、上面での渦の発生は数値的なものではなく、物理的な要因であると言えよう。
他の興味深い点は、前縁スラットの後縁から下流へと流れだした渦は小さな随伴渦を伴っている。これは、先の NACA0012 の比較的低いレイノルズ数での計算結果と定性的に同じ現象が起きていることを示している。飛行機の着陸時の騒音源として、前縁スラットが最も大きいと観察結果では報告されている。スラットの下面の剥離域が騒音源と指摘している報告もあるが、その騒音源として前縁スラットの後縁が主であるのではと示唆する結果である。
図 12 に、時間平均による圧力係数分布を実験値[5]と比較している。本計算では乱流モデルを用いてないにもかかわらず、実験値と良く一致している。翼の圧力分布は境界層の排除厚さで決まるため、今回のように非常に細かな格子を用いることで乱流モデルを用いずに圧力分布を精度良く予測できることを示していると言えよう。また、同時に高密度でかつ等方的格子を用いることの重要性を示唆していると言える。
6. おわりに現在の CFD の課題と今後の計算機の更なる性能改善、それを期待しての新しいアプローチについて議論した。また、その一つの方法として Building-Cube Method を簡単に説明し、翼型周りの流れの計算結果を示した。この計算法では、全ての方向に等間隔で、かつ物体境界近くでは極めて細かな格子を用いることで物理モデルを用いずに計算することを特徴としている。計算結果では、境界層が途中から乱れる様子、またその乱れは小さな渦の流れ下りとして捕らえられること等の詳細な流れを再現した。もちろん、高レイノルズ数流れは３次元性が重要であり、それを２次元計算で議論は危険である。しかし、定量的にも圧力係数分布が実験値と良く一致することから、高密度で等方的格子を用いることの重要性を示したと言えよう。
今回の４要素翼の計算には SX-7 の 16CPU を用いて数日の計算時間を要している。ベクトル化の改善や最適な時間刻みの選択などで計算時間はまだ短縮の余地がある。また、並列化が容易であることから多並列計算機を用いるこ とで、数時間以内で結果を得ることも可能とな ろう。Building-Cube 法による３次元計算は、有限スパン長の翼等の限られた形状では現行の計算機で可能であるが、航空機全機周りの計算にはペタフロップス級の演算性能を持つ計算機の出現を待たなくてはならない。逆に、そのような計算機ではより興味深い流れが観察できるとの期待も大きい。
B uilding-Cube Method の問題点は、膨大な格子点を扱うことである。確かに物理モデルを

10

図９. 4 要素翼型の瞬間マッハ数分布

図 10. 4 要素翼型の瞬間密度分布

Comp.
(Time-Averaged)

図 11. ４要素翼の時間平均圧力分布.

図 12. 時間平均圧力係数分布の実験値との比較 .
用いないで精度良い計算を行おうとすると詳細な格子は不可欠であるが、しかし航空機周りの流れをこのアプローチで計算するとなると 10 12 以上の格子点が必要となる。１０年後の計算機でも可能かも知れないが、工学的なツールとしてはまだ数十年先まで待たなくてはならない可能性がある。現在は、生真面目に細かな格子を用いた理想的な計算法を追求しているが、それをブレークダウンし て工学的に使えるツールとする方法の検討も今後必要であり、そのことで新たなアイデアも生まれてこよう。
現在、超音速の旅客機や空気抵抗の更に小さな飛行機、あるいはより静かな飛行機のため

11

の研究が進められている。また、GPS と自動操縦技術によって誰もが自家用機で飛び回る世界も２０年後には可能とも言われている。流体機械の最たるものとしての飛行機が、ペタフロップス級の計算機により更に高度に発展することを期待したい。
2. 謝辞
本研究は、東北大学情報シナジーセンターとの共同研究として行われた。SX7 上での計算コードの並列化・最適化および実行に際して情報シナジーセンターから適切かつ有益な助言を頂いた。この場を借りて感謝したい。
3. 参考文献[1] 中橋和博、藤田健、“非構造格子に基づく流体数値計算法の並列化”、東北大学情報シナジーセンター、SENAC 2002 年[2] Nakahashi, K., Ito, Y., Togashi, F., “Some challenges of realistic flow simulations byun structured grid CFD,” Int. J. for Numerical Methods in Fluids, Vol. 43, 2003, pp.769-78 3.
[ 3] Nakahashi, K., “Building-Cube Method for Flow Problems with Broadband Characteristic
Length,” Computational Fluid Dynamics 2002, edited by S. Armfield, R. Morgan, K. Srinivas,Springer, 2003, pp.77-81.
[ 4] Nakahashi, K., “High-Density Mesh Flow Computations with P re-/Post-Data Co mpressions,”
AIAA[5 ]05-4876, 17th AIAA Computational Fluid Dynamics Conference, Toronto, June 2005.

E. Omar, T. Zierten, M. Habn, E. Szpiro, and A. Mabal, NASA CR-2215, 1979.

12