令和２年度高性能汎用計算機高度利用事業「富岳」成果創出加速プログラム「航空機フライト試験を代替する近未来型設計技術の先導的
実証研究」
成果報告書

令和 3 年 5 月 28 日
国立大学法人東北大学
国立研究開発法人宇宙航空研究開発機構
河合

宗司

目次
1.補助事業の目的 .......................................................................... 1
2.令和２年度（当該年度）の実施内容 ........................................................ 1
2-1．当該年度（令和 2 年度）の事業実施計画................................................ 1
2-2．実施内容（成果） ................................................................... 2(1) 航空機実機複雑形状での空力課題に対する壁面モデルを用いた準直接解析 ................ 2(2) 実飛行レイノルズ数での空力課題に対する主翼基本形状の準直接解析 ................... 10(3) 主翼空力構造連成課題に対する多設計変数・多目的最適化手法の試行 ................... 14
2.3. 活動（研究会の活動等） ............................................................. 172.4. 実施体制 ........................................................................... 18
別添１

学会等発表実績 ................................................................... 19

補助事業の名称「富岳」成果創出加速プログラム
航空機フライト試験を代替する近未来型設計技術の先導的実証研究

１．補助事業の目的
航空機の設計・開発において、開発下流段階の実機フライト試験など実機製造後でないと評価できずにいる設計上の重要課題（空力課題（最大揚力・失速予測など）や主翼空力構造連成課題など）がある。このため開発下流段階で重大な手戻りが発生するリスクが生じ、開発期間の長期化等の開発リスクとなっている。本補助事業では、高度な計算科学により、ものができる前の設計初期段階でこれらの重要課題を評
価可能とすることで、航空機実機フライト試験を代替する新たな設計プロセス技術を先導的に実証することを目的とする。また先の重点課題では、前後処理を含めて高速かつ高精度に流体現象の本質をとらえる準直接解析技術を開発した。本事業では、この準直接解析を航空機実機複雑形状かつ 10 の 7 乗オーダ
ーの実飛行レイノルズ数環境下に適用し、その有効性を示し、国内航空機開発メーカーとの連携のもと民間旅客機設計開発プロセスへの展開を試みる。

２．令和２年度（報告年度）の実施内容
２－１．当該年度（令和２年度）の事業実施計画（１）航空機実機複雑形状での空力課題に対する壁面モデルを用いた準直接解析（１ａ）直交格子基盤ソルバーFFVHC-ACE における壁面モデル LES (Large Eddy Simulation) を開発
し、本手法を主翼基本形状解析に適用し、参照データベースと比較することでその解析精度を検証する。また「富岳」の本格利用が開始されれば、本手法により実機複雑形状解析が可能であることを示す予備解析を実施する。
＜実施機関：国立大学法人東北大学＞（１ｂ）直交格子基盤ソルバーFFVHC-ACE における物体適合レイヤー格子法を開発し、本手法を用いた壁面モデル LES を主翼基本形状解析に適用し、参照データベースと比較することでその解析精度を検証する。また「富岳」の本格利用が開始されれば、本手法により実機複雑形状解析が可能であることを示す予備解析を実施する。
＜実施機関：国立研究開発法人宇宙航空研究開発機構＞（１ｃ）直交格子基盤ソルバーFFVHC-ACE の高速化・機能強化および利便性向上を実施する。FFVHCACE の大規模並列性能評価とノード性能の高速化、およびデータハンドリングの効率化を実施する。
＜実施機関：国立研究開発法人宇宙航空研究開発機構＞（２）実飛行レイノルズ数での空力課題に対する主翼基本形状の準直接解析「富岳」の本格利用が開始されれば、高次精度・構造格子ソルバーを用いて実飛行レイノルズ数・主翼基本形状の大規模 LES 解析が可能であることを示す予備解析を実施する。
＜実施機関：国立大学法人東北大学＞-1-

（３）主翼空力構造連成課題に対する多設計変数・多目的最適化手法の試行「富岳」の本格利用が開始されれば、主翼空力構造連成最適化解析に向けた第一段階である、主翼
内部構造変数（構造解析に用いる部材厚み等）を設計変数とした最適化解析が可能であることを示
す予備解析を実施する。
＜実施機関：国立大学法人東北大学＞
プロジェクトの総合的推進
プロジェクト全体の連携を密としつつ円滑な運営のため、研究協力機関との実施者会議や統括会議などを開催し、研究協力機関や連携機関の連携・調整にあたる。国内外の関連課題との連携や、産業界の実ニーズの把握をタイムリーに取り込むため、海外連携大学研究者および連携企業研究者と定期的に交流する。これらの目的のために、研究会や学会企画セッションなどを実施する。また、プロジェクトで得られた成果は論文発表やホームページでのデータベースの公開など積極的に公開する。若手研究員（ポスドク等）については、有能な人材を確保し、育成する計画を継続する。これに伴い、若手研究員の連携、将来のステップアップまで見据えた登用や人材育成の取り組みを継続していく。

２－２. 実施内容（成果）（１）航空機実機複雑形状での空力課題に対する壁面モデルを用いた準直接解析（１ａ）直交格子基盤ソルバーFFVHC-ACE における壁面モデル LES の開発および検証
直交格子基盤ソルバーFFVHC-ACE に対する壁面モデル LES の実装・検証
本年度は、物体非適合格子用の壁面モデル LES を FFVHC-ACE (FrontFlow/Violet Hierarchical Cartesian Aeronautics based on Compressible-fluid Equations)を開発し、その検証を実施した。壁面モデル LES では、
壁面せん断応力𝜏𝑤 を壁面平行方向の運動量流束として与えることで乱流境界層を再現している。しかし
物体非適合格子では、保存則が厳密に満たされないため、物体境界から過剰な運動量流入・流出が起こり、壁面モデルが正しく働かなくなる。この物体非適合格子で保存則が満たされない問題を改善するため、本壁面モデル LES では、Fig. 1a.1 のように、壁面モデルに入力を与える点 (Image Point、図中 IP)よ
り下側で速度が線形となるように壁面で滑り速度を与える。すなわち、IP より下側の速度分布は𝑈(𝑌) = 𝑈𝐼𝑃 −𝑑𝑈| (𝑌 − 𝑌)
𝑑𝑌 𝐼𝑃 𝐼𝑃

(1a.1)と表される。ここで、𝑈は壁面平行方向速度、𝑌は壁面垂直方向の座標である。なお、本研究では壁面モ
デルには代数モデル[1a.1]を使用し、速度勾配𝑑𝑈/𝑑𝑌は壁面モデルの解析的な微分により与える。また、
上記の滑り速度によって速度勾配が減少し、境界層内のせん断応力の不足が生じる。そこで、
𝑚𝑜𝑑𝑒𝑙
𝜏𝑖𝑗
= 𝜏𝑤 (𝑡𝑖 𝑛𝑗 + 𝑛𝑖 𝑡𝑗 )𝑓𝑠 ,𝑌𝐼𝑃 − 𝑌
𝑓𝑠 = max [, 0.0]
𝑌𝐼𝑃(1a.2)というモデルせん断応力項を陽的に加えることで、乱流境界層内の全せん断応力バランスを維持できるように設計している。ここで、𝑡𝑖 、𝑛𝑖 はそれぞれ壁面平行方向、壁面垂直方向の単位ベクトルである。
上記の壁面モデル LES を実装した FFVHC-ACE の検証計算として、平板乱流境界層の解析を実施した。
解析対象は Fig. 1a.2 に示したような平板上に発達する乱流境界層であり、領域長さは2𝐿、一様流速度、-2-

長さ𝐿基準のレイノルズ数5.0 × 106である。また、𝑅𝑒𝑥 = 1.0 × 105の地点に突起(trip wire)を設けることで乱流遷移を誘起している。計算格子は格子幅がΔ𝑥/𝐿 = 6.3 × 10−3 (grid 1)およびΔ𝑥/𝐿 = 3.2 × 10−3 (grid2)の 2 ケースで解析を実施した。これらの格子幅は、計算結果を比較する𝑅𝑒𝑥 = 9.5 × 106の地点における境界層厚みに対して、それぞれ約 1/40、1/80 点の大きさとなっている。また、IP の長さを 3.5Δ𝑥および4.5Δ𝑥とし、IP 長さに対する感度も調査する。各計算においては、空間離散化には 3 次精度 MUSCL を適用した風上法を用い、時間積分法には 3 次精度 total-variation-diminishing Runge–Kutta 法を用いる。また、
サブグリッドスケールモデルは用いない。
Figure 1a.3 に𝑅𝑒𝑥 = 9.5 × 106の地点における無次元速度分布およびレイノルズせん断応力分布を示す。
各ケースで風上粘性によるとみられる若干の log-layer mismatch が見られるが、格子・IP 長さによらず概ね速度分布は log-law に沿う。また、各ケースで−𝑢′ 𝑣 ′ の境界層内ピーク値がおおよそ𝑢𝜏2 (𝑢𝜏 ≡ 𝜏𝑤 /𝜌𝑤 は摩擦速度、𝜌𝑤 は壁面での密度)に等しくなる。これは、壁面モデルが正しく実装できていること、つまり壁面で与えた𝜏𝑤 が、全せん断応力バランスを介して流体側に正しく作用していることを示している。

Fig. 1a.1 Partial-slip boundary condition.

Fig. 1a.2 Computational settings for flat-plate boundary layer cases.

(a) Mean streamwise velocity(b) Reynolds shear stressFig. 1a.3 Mean streamwise velocity and Reynolds shear stress profiles at 𝑅𝑒𝑥 = 9.5 × 106 .
Red lines, grid 1 (𝑌𝐼𝑃 = 3.5Δ𝑥); blue lines, grid 1 (𝑌𝐼𝑃 = 4.5Δ𝑥);
green lines, grid 2 (𝑌𝐼𝑃 = 4.5Δ𝑥 ); black dashed line in (a), log-law 𝑢+ = log(𝑦 + ) /0.41 + 5.1.
また、主翼基本形状周りの流れ場での壁面モデル LES 解析を実施した。解析対象は A-airfoil、翼弦長・-3-一様流速ベースのレイノルズ数は 2.1×106、
迎角は 13.3°である。
本解析対象については、
Gleyzes & Capbern[1a.2]による風洞試験や Asada & Kawai [1a.3]による LES 解析が行われている。解析結果の詳しい分析はこれからであるが、計算が安定に実施できており、境界層の乱流化や後縁での境界層剥離等が再現できることを確認した(Fig. 1a.4 参照)。さらに、航空機実機複雑形状(JAXA standard model: JSM)の壁面モデル LES
試計算にも成功している。平均空力翼弦長・一様流速ベースのレイノルズ数は 1.7×106、迎角は 18.48°であり、総セル数は約 36 億である。Figure 1a.5 に示すように、このような航空機離着陸時の複雑形状周り（スラット、フラップ、それらの支持装置など）であっても、FFCHC-ACE による壁面モデル LES 計算が安定に実施できることも予備解析で確認した。

Fig. 1a.4 Wall-modeled LES of 2D airfoil flow using Fig. 1a.5 Wall-modeled LES of flow around full-aircraft
FFVHC-ACE.

configuration (JSM) using FFVHC-ACE.

直交格子基盤ソルバーFFVHC-ACE に対する KEEP スキームの実装・実用化
これまでの FFVHC-ACE では、対流流束評価に風上スキームが用いられていた。風上スキームは複雑形状まわりの流れ場解析を安定に行うことができ、圧縮性流体解析ソルバーで広く一般的に用いられている手法である。一方、安定な解析を実現するために大きい数値粘性を加えているのも事実であり、乱流の支配的スケールを格子で直接解像する LES に適しているとは言い難い。それに対し、我々の研究グルー
プでは近年、kinetic energy and entropy preserving (KEEP) scheme を提案している[1a.4]。KEEP スキームは
中心差分系のスキームであるため、数値粘性（拡散誤差）が存在しない。更に中心差分系スキームで問題となる数値不安定性を、運動エネルギーとエントロピーの保存を離散的に満たすことで解決している。
KEEP スキームは階層型直交格子にも展開され、中心差分系スキームによる複雑形状まわりの高精度解析が期待できるまでとなっている[1a.5]。
本プロジェクトでは、FFVHC-ACE の産業界への展開も考え、数値計算手法としては、安定ではあるが数値粘性の影響が無視できない既存の風上スキーム、および数値計算手法として攻めた LES 解析が実施可能な KEEP スキームの 2 通りを実装し、ユーザーが選択できるソルバーを構築している。ここでは、階層型直交格子 KEEP スキームを FFVHC-ACE へ実装し、その性能評価を行った上で、主翼基本形状さ-4-

らには実機複雑形状の LES 解析への展開を試みる。今年度は、階層型直交格子 KEEP スキームの実装および性能評価を実施した。性能評価は、直径ベースのレイノルズ数が𝑅𝑒𝐷 = 1、000（後流が乱流化する条件）における球まわりの流れ場解析で実施した。総セル数は約 1 億の解析である。その結果、階層型直交格子 KEEP スキームを用いることで、数値粘性を付加していないにもかかわらず安定な解析が実現できることを確認した。一方で、格子細分化レベルが切り替わるブロック境界（以下、階層境界と呼ぶ）で、計算が破綻するまでではないが、非物理的な数値振動が生じることも明らかとなった（Fig. 1a.6 (a)参照）。
またこの数値振動は、階層境界の取り扱いが原因であることも明らかとなっている。
この数値振動の問題を解決するために、本研究では数値的 low-pass フィルターで高周波成分の数値振動除去を行うことを考えた。フィルターは数値解析では広く用いられているものの、階層型直交格子に直接用いると階層境界で保存則を満たさなくなるという問題がある。そこで、本研究では流束形式のフィル
ター[1a.6]を階層型直交格子に展開し、階層境界で保存則を満足するフィルターを開発し、FFVHC-ACE へ実装した。2𝑛次精度の流束形式フィルターは以下の式で与えられる。
𝑛−1𝑞̃|𝑚 Δ𝑉 = 𝑞|𝑚 Δ𝑉 − (𝑓|𝑚+1 − 𝑓|𝑚−1 ) Δ𝑆,
2

2(1a.3)

𝑓|𝑚+1 Δ𝑆 = ∑ 𝑏𝑘 (𝑞|𝑚+1+𝑘 − 𝑞|𝑚−𝑘 )Δ𝑉 .
2

𝑘=0

ここで、𝑞は保存変数、𝑞̃はフィルター後の保存変数、𝑓はフィルター流束、𝑏𝑘 は係数、Δ𝑆はセル境界面の面積、Δ𝑉はセル体積であり、下添字𝑚はセル番号を示している。階層境界での保存則は、細かいセルからのフィルター流束の和が粗いセルへのフィルター流束に等しくなるように設定することで満たしている。
一般に、高次精度フィルターは解像できている低波数に大きな悪影響を与えずに高波数のみを減衰させることができるが、ステンシル幅が大きいため多くの階層境界や壁面によるステンシル切断を有する階
層型直交格子や複雑形状に適用するには実装が煩雑となる。そこで本研究では、比較的低次精度な 4 次
精度フィルターを用い、フィルター強さを調節可能にすることを考えた。フィルター強さは、流束形式であることを活用し、フィルター流束を𝛼倍(0 ≤ 𝛼 ≤ 1)することで調節した。なお、理論的なフーリエ解析
でフィルター流束を𝛼倍すると、単純にフィルター強さが𝛼倍されることを確認している。Figure 1a.6 (b)に提案するフィルターを用いたときの球まわりの流れ場解析の結果を示す。フィルターなしで生じていた階層境界での数値振動を除去することに成功した。これらの成果により、階層型直交格子 KEEP スキ
ームとフィルターを実装した FFVHC-ACE で、高精度な LES 解析が実現できると期待できる。

(a) Without low-pass filter(b) With the proposed low-pass filter

Fig. 1a.6 Streamwise velocity contours around a sphere (𝑅𝑒𝐷 = 1000) obtained by FFVHC-ACE with KEEP
scheme.

-5-参考文献[1a.1] S. R. Allmaras and F. T. Johnson, “Modifications and clarifications for the implementation of the Spalart-Allmaras turbulencemodel,” Seventh international conference on computational fluid dynamics (ICCFD7). (2012).
[1a.2] C. Gleyzes and P. Capbern, “Experimental study of two AIRBUS/ONERA airfoils in near stall conditions. Part I: Boundary layers,”
Aerospace science and Technology, Vol. 7, pp. 439-449, (2003).
[1a.3] K. Asada and S. Kawai, “Large-eddy simulation of airfoil flow near stall condition at Reynolds number 2.1 × 106,” Physics of
Fluids, Vol. 30, pp. 085103, (2018).
[1a.4] Y. Kuya, K. Totani and S. Kawai, “Kinetic energy and entropy preserving schemes for compressible flows by split convective
forms,” Journal of Computational Physics, Vol. 375, pp. 828-853, (2018).
[1a.5] Y. Kuya and S. Kawai, “A stable and non-dissipative kinetic energy and entropy preserving (KEEP) scheme for non-conforming
block boundaries on Cartesian grids,” Computers & Fluids, Vol. 200, pp. 104427, (2020).
[1a.6] T. Shiroto, S. Kawai and N. Ohnishi, “Finite-volume-concept-based Pade-type filters,” Journal of Computational Physics, Vol.
349, pp. 215-219, (2017).

（１ｂ）直交格子基盤ソルバーFFVHC-ACE における物体適合レイヤー格子法の開発および検証
直交格子基盤ソルバーFFVHC-ACE における物体適合レイヤー格子法の開発を行った。直交格子と物体適合レイヤー格子を組み合わせた解析を実施するための必須ツールとして、新しく領域分割ツールを開発した。物体適合レイヤー格子の流体解析ソルバー（以後、レイヤーソルバーと呼ぶ）には、本プロジェク
トで開発した平衡壁面モデルを実装している。平衡壁面モデルでは壁から垂直方向に数セル離れた Image
Point(IP)と呼ばれる点での物理量を入力とし、壁における摩擦力（𝜏𝑤 ）を計算し、この値をレイヤーソル
バーで利用することになる。並列計算を行うために物体適合レイヤー格子も領域分割されることになるが、壁と対応する IP が別領域に分割されると、新たにプロセス間通信が必要となる。これを回避するた
めに、壁と対応する IP が同じ領域になるように領域分割を行うツールを開発した。また、IP における物理量の計算手法として、一般的には周囲の物理量から補間を行う手法が用いられるが、補間に用いる物理量定義点が領域外や物体内部に存在する場合は補間に利用できず、その様な場合を考慮する必要がある
が処理が複雑となり高速化の妨げになる。ここでは、あらかじめ IP の最近傍点を探索しておき、その点にける物理量の微分値から IP の値を外挿する方法を適用した。
開発したプログラムを用いて、基礎的な形状として平板および主翼基本形状に対して解析精度の検証を開始した。Fig. 1b.1 に主翼基本形状（2 次元断面形状は A-airfoil）ならびに平板での解析例を示す。

-6-

(b) Flat plane (iso-surfaces of density colored by(a) 2D airfoil (Mach number distributions)streamwise velocity)(c) Flat plane flow (streamwise velocity distributions, upper, BCM grid; lower, layer grid)Fig. 1b.1 Flow Simulations for basic geometries using layer grid.
（１ｃ）直交格子基盤ソルバーFFVHC-ACE の高速化・機能強化および利便性向上
直交格子基盤ソルバーFFVHC-ACE の高速化チューニング
ノード性能の高速化に関しては、RIST の高度化支援で実施した 4 つのカーネルチューニングを参考に、
対流項流束の差分を計算するルーチンの高速化チューニングを FFVHC-ACE に適用した。当該箇所は 3重ループのストリーミング計算であったため、多重ループの 1 重化、ループブロッキング、コンパイラによるループ分割を適用した。Fig. 1c.1 に高速化チューニング適用前後での CPU 性能解析レポートによる分析結果を示す。この図で、横軸はあるプロセス内の各スレッドを、縦軸は色で示されたそれぞれの処理に掛かった経過時間を示している。ループ 1 重化によるメモリアクセス性能の向上、ループブロッキン
グによる L2 キャッシュ性能の向上、ループ分割による浮動小数点演算待ちの減少などのチューニングによる改善効果がわかる。これらの高速化チューニングにより当該箇所の処理速度が 1.5 倍高速化した。
ノード並列性能の評価を行った。結果を Fig. 1c.2 に示す。500 ノードで並列化効率が 50%前後とあまり良くない結果が得られた。本計測では弱スケーリングでの計測を行っているが、本来なら通信処理（図では mpi）のオーバーヘッドにより並列化効率が減少するのが一般的であるが、ここでは通信処理以外の処理で、特にプロセス間のばらつきが原因であることがわかった。多ノード（具体的には数百ノード以上）
でのプロセス間の処理時間のばらつきが発生する原因は調査中である。

-7-

(a) Before tuning (asis)(b) After tuning
Fig. 1c.1 Results of high-speed tuning (analysis by CPU performance analysis report).

-8-

Fig. 1c.2 Evaluation of parallel performance (weak-scaling).
物体非適合・直交格子基盤ソルバーFFVHC-ACE における物体表面積分手法の構築・検証
物体に働く力など物体表面での物理量の積分や表面での物理量分布の予測精度の検証を基礎的な 2 次
元形状に対して行った。特に物体表面における積分手法に関しては、従来用いる階段状近似ではなく物体
表面の法線ベクトルを考慮した手法を開発し[1c.1]、その精度の検証を行った。Fig. 1c.3 は薄翼（NACA4402）に対し、a) 表面積の誤差、および b) 面ベクトルの総和（物体が閉じているかどうか）を示している。表面積に関しては階段状近似（図における STEP）ではフラクタル性により格子を細かくしても正しい値に収束しないが、提案手法（図における CUTD）では、格子サイズを細かくすることで正しく求めることが可能である事がわかる。一方、面ベクトルの総和に関しては、階段状近似では誤差の範囲で 0 となり、提案手法でも格子を細かくすることで 1 次精度で 0 になることがわかる。以上のように本手法では、直交格子法が苦手な薄い物体でも問題なく表面積分が可能であることを確認した。
100

102

10-2
1

10-4

10-6
2 1/2

100

CUTD(circle)
10-8

STEP(circle)
CUTD(naca4402)

2((・Sx) +(・Sy) )Area error [%]

10

10-1

STEP(naca4402)

-10

10

1st-order

10-12

10-2

10-14

CUTD(circle)
STEP(circle)
CUTD(naca4402)

10-16

STEP(naca4402)
1st-order
10-3
10-4

10-3

10-2

10-18
10-4

10-1

10-3

dx (cell size)

10-2

10-1

dx (cell size)(a) Errors of surface area(b) Sum of surface normal vectors(c) 2D airfoil (NACA4402)
Fig. 1c.3 Evaluation of surface integration method-9-Fig. 1c.4 に NACA4402 の表面圧力（𝐶𝑝 ）分布および表面摩擦力（𝐶𝑓 ）分布に関して、従来用いられる物
体適合格子を用いた解析結果（UPACS）と比較した結果を示す。BCM (L08, L10, L12, L14) が本手法による結果で、L は直交格子法で用いる 8 分木の階層を表し、大きくなるほど格子が細かくなることを示している。L10 以上で物体適合格子を用いた解析結果と良く一致していることがわかる。

-3

-0.5

BCM(L14)

BCM(L14)

BCM(L12)

BCM(L12)

BCM(L10)-2

BCM(L10)

BCM(L08)

BCM(L08)

-0.4

UPACS

UPACS

-1

Cp

Cp

-0.3

0
-0.2

1
-0.1

2

-0.01

0

0.01

0.02

0.03

0.04

0
0.95

0.05

0.96

0.97

(a) 𝐶𝑝 distributions near leading edge
0.5

0.99

1

1.01

(b) 𝐶𝑝 distributions near trailing edge
0.3

BCM(L14)
BCM(L12)

BCM(L14)
BCM(L12)

BCM(L10)
0.4

0.98
x/c

x/c

0.25

BCM(L08)

BCM(L10)
BCM(L08)

UPACS
0.2

UPACS

0.3

Cf

Cf

0.15

0.2

0.1

0.05

0.1

0

0

-0.01

0

0.01

0.02

0.03

0.04

-0.05
0.95

0.05

0.96

0.97

0.98

0.99

1

1.01

x/c

x/c

(c) 𝐶𝑓 distributions near leading edge(d) 𝐶𝑓 distributions near trailing edge

Fig. 1c.4 Comparisons of surface physical quantity distributions for NACA4402.
参考文献[1c.1] Ryoji Takaki, Soshi Kawai, Yuma Fukushima, Yoshiharu Tamaki, Seiji Tsutsumi and Hisaichi Shibata, “Development of a highspeed and high-precision turbulent flow solver using hierarchical cartesian grids,” Proceedings of the 51st Fluid Dynamics Conference/ the 37th Aerospace Numerical Simulation Symposium, pp165-171, JAXA-SP-19-007, (2012).

（２）実飛行レイノルズ数での空力課題に対する主翼基本形状の準直接解析
航空機実飛行レイノルズ数における LES 解析
航空機実飛行レイノルズ数(翼弦長𝑐、一様流速𝑢∞ ベースのレイノルズ数𝑅𝑒𝑐 = 1.0 × 107)、失速直前の条件(迎角α = 13。3∘)、主翼基本形状(A-Airfoil)における通常の壁面まで解像する大規模 LES 解析を実施- 10 -

した。本年度は、9 月までに計算コード、結果処理環境の整備や粗い格子における試計算を進め、追加計
算資源の配分を受け 10 月以降に本計算の実施を進め、当初計画以上に研究が進展した。計算格子は Fig.
2.1 に示した C 型構造格子であり、スパン方向領域サイズは過去の計算結果に基づき十分だと考えられる
大きさ(コード長の 4.8%)としている。また、壁面平行方向の格子幅については、壁面近傍の乱流境界層のストリーク構造を十分に解像できるよう、壁面粘性スケールでΔ𝜉 + < 25 (コード方向)、 Δ 𝜂 + < 13 (スパン方向)となるように設定した。この時、各方向の格子点数は 18,937 (コード方向)、3,204 (スパン方向)、
629 (垂直方向)、総格子点数は約 382 億点である。計算には「富岳」8,640 ノードを使用し、MPI (34、560
分割領域)及び OpenMP (12 スレッド)を併用して計算を並列化している。空間離散化には 6 次精度コンパクト差分法 [2.1]、エイリアンシングエラー等による数値不安定性を抑制するために 8 次精度 3 重対角フ
ィルタ [2.2]を使用した。時間積分法は 2 点後退差分を適用した陰解法[2.3, 2.4]であり、内部反復回数は 5回である。なお、本計算はサブグリッドスケール乱流モデルを用いていない。
Figure 2.2 に瞬時流れ場(速度勾配テンソル第二不変量𝑄の等値面)を示す。翼前縁付近では層流から乱流への遷移が見られる他、翼の後半部分では逆圧力勾配によると考えられる境界層の急速な成長が見られる。次に、壁面摩擦係数𝐶𝑓 と壁面圧力係数𝐶𝑝 の翼面上分布について、過去に実施した風洞試験スケールレ
イノルズ数(𝑅𝑒𝑐 = 2.1 × 106 )における LES 解析[2.5]と比較する。なお、時間方向の計算の収束は現在確認中であるが、以下では本年度中に取得できた𝑇𝑎𝑣𝑒 𝑢∞ /𝑐 = 0.6 (𝑇𝑎𝑣𝑒 : 平均時間)分の結果について Fig. 2.3に示す。𝐶𝑓 の分布から分かるように、𝑅𝑒𝑐 = 2.1 × 106の時と比べて乱流遷移点(𝑥/𝑐 = 0.1付近の急激な𝐶𝑓の増加)が前方に移動しており、層流剥離(𝐶𝑓 < 0の領域)が観察されない。同時に、後縁付近でも剥離がほぼ消失する。この結果として翼型の圧力分布は非粘性解(図中黒点線)に近づき、𝑅𝑒𝑐 = 2.1 × 106の時と比べて前縁付近の負圧ピークが強くなり、後縁での圧力回復の促進が見られる。またその結果、レイノルズ数が風洞試験レベルから実機飛行レベルに上がることで揚力が高くなる。これらの観察結果は過去の実
験的研究 [2.6, 2.7]における典型的なレイノルズ数効果の傾向と一致しており、本解析によって得られる
詳細な乱流データベース(統計平均量・瞬時流れ場)がレイノルズ数効果を検証する際の良いデータベースとなることが期待される。本研究では、当初の計画以上に研究が進展している。当初計画では、実飛行レ
イノルズ数・主翼基本形状の大規模 LES 解析が可能であることを示す予備解析を実施する予定であったが、予備解析を完了し、「富岳」の追加リソース配分により計画を前倒しし、本解析を開始している。

Fig. 2.1 Computational grid around A-airfoil (every 50 grid points are shown).

- 11 -

Fig. 2.2 Iso-surfaces of the Q-criterion colored by streamwise velocity in near-stall condition at 𝑅𝑒𝑐 = 1.0 × 107.
-6

Cp

-4

0

0.015

0.4

-2

0.8

0.9

x/c

Cf

0.2

-3

Cp

0.02

-0.2

-5

1

0.01

0.005

-1
0

0

1
0

0.2

0.4

x/c

0.6

0.8

1

0

0.2

0.4

x/c

0.6

0.8

1

Fig. 2.3 Mean surface pressure (left) and skin friction coefficient (right) along the airfoil. Red lines, LES at 𝑅𝑒𝑐 =1.0 × 107; blue lines, LES at 𝑅𝑒𝑐 = 2.1 × 106; Black dashed line, inviscid solution obtained by XFOIL [2.8].
データ科学を用いた超大規模時系列データ解析手法の構築
上記のように、今まで困難であった航空機実飛行レイノルズ数（𝑅𝑒𝑐 = 1.0 × 107）の通常の壁面まで解像する LES 解析が「富岳」を用いて実現可能となってきている。しかし、このような超大規模解析で得られる複雑流体現象の超大規模データから、どのように流体現象を理解するかについては、未だ多くの困難がある。例えば、乱流統計量（平均流速やレイノルズ応力など）を解析し流体現象を理解する方法が広く用いられているが、統計量だけでは流体構造や現象のきっかけなどを理解することは難しい。一方で、瞬間場を解析する方法では、細かい構造から大きな構造まで含む瞬間場から、現象を支配する主要な構造を理解することは難しい（特にデータが超大規模になると困難さを極める）。そこで本研究では、このよ
うな超大規模データを解析する方法論として、データ科学的アプローチの開発・適用を試みている。具体的には、POD [2.9]や DMD [2.10]を用いた主要構造の低次抽出、および Resolvent 解析[2.11]を用いた
Input/Output 解析手法の確立を実施している。今年度は、POD/DMD ソルバーを構築し、遷音速バフェッ
ト壁面モデル LES 解析データへ適用することで、その妥当性と適用可能性を検証した。また、Resolvent
解析ソルバーも構築し、ポアズイユ流れへ適用することでその妥当性検証を行った。
Figure 2.4 に遷音速バフェットの壁面モデル LES 解析で得られた圧力変動と、それに対する DMD 解析
結 果 を 示 す 。 遷 音 速 バ フ ェ ッ ト の 壁 面 モ デ ル LES 解 析 は 「 京 」 を 用 い て 行 わ れ 、 格 子 点 数 が4603×565×169≒4.4 億の大規模解析であり、本 DMD 解析ではそのデータに対してスパン方向に平均をとった 2 次元データ(格子点数 4603×169≒77.8 万)を用いた。また、DMD 解析で抽出する周波数モードはバ- 12 -

フェット周波数（ストローハル数𝑆𝑡 = 0.065）のモードとした。各モードは、遷音速バフェット現象に特徴的な剥離せん断層上部の圧力変動などが細かい構造を取り除かれた状態で抽出されている。また、
POD 解析で支配的な構造を抽出すると、DMD 解析と同様の流れ場構造が得られることも確認している。
これより、大規模解析で得られたデータから重要な流体構造を抽出できる可能性が示唆された。なお、より簡単な円柱まわりの流れ場にも POD/DMD 適用した際には特徴的なカルマン渦の構造が抽出できており、本ソルバーの妥当性は検証済である。
また、Fig. 2.5 にポアズイユ流れに対して Resolvent 解析を行った結果を示す。Fig. 2.5 (a)は、forcing モード(input)の周波数に対する resolvent gain (output への最大増幅率)を示しており、Sun らによる解析結果[2.12]も合わせて示している。本解析結果は Sun らの結果とよく一致しており、本解析の妥当性が示される。また、Fig. 2.5 (b)は resolvent gain が最大となる周波数の forcing/response モードを示しており、forcing
モードのような擾乱を流れに加えると、response モードのような増幅が得られることを意味する。Nishioka
らの実験[2.13]においても同様な擾乱の増幅が見られており、これからも本 Resolvent 解析ソルバーの妥当性が確認される。以上より、Resolvent 解析ソルバーの基盤構築が完了し、今後は遷音速バフェットのLES 解析データなど大規模データへと適用する予定である。

(a) Original flow (Wall-modeled LES)(b) DMD mode

Fig. 2.4 Pressure fluctuations of transonic airfoil buffet obtained by wall-modeled LES and DMD analysis。

(a) Resolvent gain(b) Forcing mode
Fig. 2.5 Resolvent analysis for Poiseuille flow.

- 13 -(c) Response mode

参考文献[2.1] S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, Vol. 103, No. 1,
pp. 16–42, (1992).
[2.2] D. V. Gaitonde and M. R. Visbal, “Pade-type higher-order boundary filters for the Navier–Stokes equations,” AIAA Journal, Vol.
38, No. 11, pp. 2103–2112, (2000).
[2.3] S. Obayashi, K. Fujii and S. Gavali, “Navier–Stokes simulation of wind-tunnel flow using LU-ADI factorization algorithm,”
NASA-TM-100042, NASA, (1988).
[2.4] N. Iizuka, “Study of Mach number effect on the dynamic stability of a blunt re-entry capsule, ” Ph.D. thesis, University of Tokyo,(2006).
[2.5] K. Asada and S. Kawai, “Large-eddy simulation of airfoil flow near stall condition at Reynolds number 2.1 × 106, ” Physics of
Fluids, Vol. 30, No. 8, pp. 085103, (2018).
[2.6] E. N. Jacobs and A. Shermann, “Airfoil section characteristics as affected by variations of the Reynolds Number, ” NACA Report
No. 586, (1936).
[2.7] R. M. Pinkerton, “The variation with Reynolds number of pressure distribution over an airfoil section, ” NACA Report No. 613,
(1938).
[2.8] M. Drela, “XFOIL: An analysis and design system for low Reynolds number airfoils,” Low Reynolds number aerodynamics,
Springer, Berlin, Heidelberg, pp. 1-12, (1989).
[2.9] J. L. Lumley, “The structure of inhomogeneous turbulent flows,” Proceedings of the International Colloquium on the Fine Scale
Structure of the Atmosphere and Its Influence on Radio Wave Propagation, edited by A. M. Yaglam, and V. I. Tatarsky, Doklady
Akaademii Nauk SSSR, Nauka, Moscow, (1967).
[2.10] P. J. Schmid and J. Sesterhenn, “Dynamic mode decomposition of numerical and experimental data,” Journal of Fluid Mechanics,
vol. 656, pp. 5-28, (2010).
[2.11] B. J. McKeon and A. S. Sharma, “A critical layer model for turbulent pipe flow,” Journal of Fluid Mechanics, vol. 658, pp. 336382, (2010).
[2.12] Y. Sun, “Global stability analysis and control of compressible flow over rectangular cavities,” PhD Thesis, Florida State University,(2017).
[2.13] M. Nishioka, S. Iida and Y. Ichioka, “An experi- mental investigation of stability of plane Poiseuille flow”, Journal of Fluid
Mechanics, vol. 72, pp. 731–751, (1975).

（３）主翼空力構造連成課題に対する多設計変数・多目的最適化手法の試行
多設計変数・多目的最適化手法の構築
多設計変数・多目的関数で構成される最適化問題の高速解法を目指して、遺伝的アルゴリズム（GA）と深層ニューラルネットワーク（DNN）からなるデータ駆動型の最適化手法（Fig. 3.1）を開発した。本
開発手法は、実験計画による既存のサンプルデータを学習することで、設計変数・目的関数間の主要な論
理構造を DNN 上に表現する。これにより、最適化問題を単純化（例えば、次元縮約）できる。さらに、
この DNN を用いて、任意の設計変数値に対する目的関数値を瞬時に予測することができるため、最適化計算に要する時間を短縮できる。DNN をベースとして求められた最適解のデータを既存のデータに追加- 14 -

し、これを再度学習することで DNN を更新する。以上のデータ駆動のプロセスを繰り返すことで、大規模な最適化問題の高速解法に繋げる。

Fig. 3.1 Data-driven optimization approach with GA and DNN.
テスト関数問題（2 目的関数，30 設計変数）を解いて得られた Pareto 最適解の結果を Fig. 3.2 に示す。
本問題は 2 つの目的関数（f1 および f2）を最小化する問題であり、プロットの左下方向に位置する解を探索できるほど、最適化手法の性能が優れていることを意味する。また、真の Pareto 最適解（Pareto 面）
は Fig. 3.2(a)に示す位置に存在する。Fig. 3.2(b)で、青のプロットは本開発手法（GA+DNN）、赤のプロッ
トは従来手法（GA 単独）の結果を示す。従来手法で探索された解は真の Pareto 面に未到達の段階である
一方、本開発手法で探索された解は既に真の Pareto 面に到達している。この結果から、本開発手法について、多設計変数問題の解法として優れた性能（探索の高速化、計算コストの削減効果）を確認することができた。

(a) True Pareto front(b) Solusions searched by two different optimizers

Fig. 3.2 Numerical experiments in a test function problem (minimizing f1 and f2).
主翼内部構造変数を設計変数とした最適化解析の定式化
次年度に実施を予定している、空力変数を固定し、主翼内部構造変数（構造解析に用いる部材に関する変数等）を設計変数とした多設計変数・多目的最適化の試行解析に向けて、先行的に Fig. 3.3 に示す
主翼構造を対象として、最適化解析の定式化を検討した。またその結果、下記の定式化を導出すること
ができた。
⚫

設計変数

- 15 -

⚫

➢

リブの配置（30 設計変数）➢

翼形状（6 設計変数、ただし最初の最適化計算では固定）

目的関数➢

構造重量の最小化➢

揚力最大化（設計変数に応じて釣り合い状態の構造変形は変化し、空力荷重も変化する為）

Fig. 3.3 Main wing box structures.
またこれに併せて、空力構造連成解析のための計算ツールを整備した。本計算ツールの概要を Fig. 3.4に示す。構築したツールでは、まず翼平面形状と翼内部構造（桁配置、リブ配置・本数）、また使用する材料（本解析では炭素繊維強化複合材 CFRP に着目する）を定める。なお、次年度に行う翼形を固定した最適化解析では、翼内部構造は設計変数として前節の DNN をベースとした GA により更新されていく。その後，空力解析により得られた空力荷重を構造解析に引き渡し、与えられた荷重に耐えられるよう構造部材の寸法を繰り返し更新（構造サイジング）する。全ての部材の安全余裕が所望の範囲（本解析では 0.45-0.55）に収まった段階でサイジングを終了とし、その時の構造の変形を反映させた空力解析へと移る。以上のように空力・構造解析を繰り返し行い、構造変形が収束（各接点の変形量が 10 の-5 乗以下）
した段階で、空力と構造の釣り合い状態の設計（静空弾設計）が得られたとみなす。

Fig. 3.4 Flowchart of aerostructural sizing for wing box structures
Fig. 3.5 に、本ツールを用いて行った予備解析の結果を示す。予備解析では、CFRP の炭素繊維の変化- 16 -

に応じて、構造重量がどのような違いを示すかを明らかにした。具体的には、T700S、T800S、T1100Gの繊維を用いた CFRP（樹脂には#3900-2B を仮定）
による空力構造連成解析及び構造サイジングを行い、Fig. 3.5(a)に示すように、最も剛性の高い繊維である T1100G を用いた場合に最軽量の構造が得られることが示された。この時、特に下面の板厚が薄くなることによって T1100G を用いた場合に構造重量が軽くなるが、逆に上面では剛性の高い繊維が厚い板厚となり、興味深い結果となった（Fig.3.5(b)）
。以上のように、繊維剛性と構造重量の間には正の相関が見られ、構築したツールにより妥当な結果を得ることが出来たと言える。以上の予備解析結果に基づき、本ツールを用いて部材配置と翼平面形状を設計変数とした最適化問題へ取り組んでいくことが可能であると判断出来たため、次年度以降これに取り組んでいく。

(a) Structural weight of different CFRPs(b) Thickness of skin with different CFRPs

Fig. 3.5 Preliminary result of aerostructural sizing of wing box structures

２－３．活動（研究会の活動等）
本課題を含め「富岳」成果創出加速プログラム「領域 3 産業競争力の強化」の 3 課題合同で「第 4 回
HPC ものづくり統合ワークショップ」を令和 2 年 11 月 27 日（金）10:00〜17:00 にオンラインで開催した（3 課題は、「富岳」を利用した革新的流体性能予測技術の研究開発（課題代表者：加藤 千幸）、スー
パーシミュレーションと AI を連携活用した実機クリーンエネルギーシステムのデジタルツインの構築と活用（課題代表者：吉村 忍）、および本課題）。
「富岳」利用の期待セッションでは、「航空機フライト試験を代替する世界初の実飛行レイノルズ数・航空機全機 LES 解析に向けて」という題目で、課題代表者の河合宗司から本課題の概要と研究成果について講演した。また「富岳」利用の展望セッションでは、
協力機関代表者の高木亮治から「FFVHC-ACE の高速化に向けた取り組み」という題目で、開発している
直交格子基盤ソルバーFFVHC-ACE の高速化チューニングについて講演した。参加機関等は 95 機関（官庁・団体・研究機関 11 機関、大学関係 18 機関、民間企業 62 社、その他 4）で、参加人数は 157 名であり多数の方々にご参加いただけた。
また日本航空宇宙学会主催の第 53 回流体力学講演会／第 39 回航空宇宙数値シミュレーション技術シンポジウムでは、企画セッション「高レイノルズ数流れと空気力学」を企画し、本課題から 3 件の講演を予定していたが、新型コロナウイルス感染症の影響で学会が中止となった。
プロジェクト全体の連携を密としつつ円滑な情報共有や運営のため、6 月、9 月、12 月、3 月に研究協
力機関（JAXA）および連携機関（三菱重工グループ）と全体推進会議を開催した。またデータ科学を用いた超大規模時系列データ解析手法の構築で連携している海外連携機関（UCLA）とも定期的な実施者会- 17 -

議を開催し、連携を密とした。

２－４．実施体制
業務項目

担当機関

担当責任者

（１）実機複雑形状の壁面モデ
ルLES解析
（１ａ）FFVHC-ACEにおける壁

国立大学法人東北大学

河合宗司（１ｂ）レイヤー格子法を用い国立研究開発法人宇宙航空研究高木亮治

たFFVHC-ACEにおける壁面モデ

開発機構

面モデルLESと実機形状予備解析ルLESと実機形状予備解析（１ｃ）FFVHC-ACE の高速

国立研究開発法人宇宙航空研究高木亮治

化・機能強化・利便性向上

開発機構（２）実飛行レイノルズ数で

国立大学法人東北大学

河合宗司

国立大学法人東北大学

下山幸治の主翼基本形状の LES の予備
解析（３）主翼空力構造連成多設
計変数・多目的最適化手法の試行- 18 -

別添１ 学会等発表実績
査読付き学術論文
1.

Y. Tamaki, Y. Fukushima, Y. Kuya and S. Kawai, “Physics and modeling of trailing-edge
stall phenomena for wall modeled large-eddy simulation,” Physical Review Fluids, Vol. 5,
Issue 7, pp. 074602, (2020).

解説論文
1.

河合宗司, “最新の圧縮性 LES 研究と｢富岳｣で航空機開発に第二のイノベーションを, ” 計算工学ナビ, Vol. 20, pp. 6, (2021).

2.

河合宗司, “実飛行レイノルズ数・航空機前期 LES 解析に向けて, ” 計算工学, Vol. 26, No. 1,
pp.4195-4199, (2021).

学会発表
1.

玉置義治, 河合宗司, “Rec~107 の翼型失速流れの LES 解析：「富岳」を用いた高精度乱流デー
タベースの構築,” 第 34 回数値流体力学シンポジウム（オンライン開催）, 2020 年 12 月 21 日.

2.

淺田啓幸, 玉置義治, 久谷雄一, 河合宗司, “階層型直交格子 KEEP スキームにおけるハンギン
グノード面での高次補間法の提案と実用性,” 第 34 回数値流体力学シンポジウム（オンライン開催）, 2020 年 12 月 21 日.

3.

高木亮治, 河合宗司, 久谷雄一, 玉置義治，“埋め込み境界法における二次元薄翼表面での解析
精度の検証”, 流体力学講演会/航空宇宙数値シミュレーション技術シンポジウム 2020 オンライ
ン, 2020 年 9 月 30 日.

4.

高木亮治, 河合宗司, 久谷雄一, 玉置義治，““階層型等間隔直交構造格子方の精度検証 ―直交格子 CFD ワークショップ―”, 流体力学講演会/航空宇宙数値シミュレーション技術シンポジウ
ム 2020 オンライン, 2020 年 9 月 30 日.

5.

Koji Shimoyama, “Data-Driven Analysis and Design in Mechanical and Aerospace
Engineering,” International Joint Conference on Aerospace and Aviation (ICASA) Mechanical and Advance Materials Research Exchange (MATRIX) 2020, online, 26
November 2020 (keynote).

- 19 -