令和 3 年度高性能汎用計算機高度利用事業「富岳」成果創出加速プログラム「『富岳』を利用した革新的流体性能予測技術の研究開発」
成果報告書

令和 4 年 5 月 28 日
国立大学法人東京大学 生産技術研究所
加藤

千幸

目次
補助事業の名称 ........................................................................................................................................ - 1 １．補助事業の目的 ................................................................................................................................. - 1 ２．令和３年度（報告年度）の実施内容 ................................................................................................ - 1 ２－１．当該年度（令和３年度）の事業実施計画.................................................................................. - 1 ２－２．実施内容（成果） ...................................................................................................................... - 3 (1) 数値曳航水槽の実現と省エネデバイスによる推進効率の向上 ...................................................... - 3 (2) 細隙部を含めた多段遠心ポンプの内部流れの Wall-Resolved LES............................................. - 18 (3) 圧縮機サージの直接解析............................................................................................................... - 27 (4) リアルワールド自動車空力性能の予測 ......................................................................................... - 40 (5) リアルワールド自動車空力音予測 ................................................................................................ - 50 (6) プロジェクトの総合的推進 ........................................................................................................... - 75 ２－３．活動（研究会の活動等） ......................................................................................................... - 76 ２－４．実施体制................................................................................................................................... - 77 別添１ 学会等発表実績 ........................................................................................................................ - 79 --2-

補助事業の名称「富岳」成果創出加速プログラム「富岳」を利用した革新的流体性能予測技術の研究開発

１．補助事業の目的「富岳」の有する高い計算性能を十二分に引き出すことができるアプリケーション・ソフトウェア（以下、単にアプリケーション）を駆使することにより、壁面近傍の微細な乱流渦の直接数値計算（Wall-Resolved LES（Large-Eddy Simulation））、および、これらの乱流渦の効果を従来よりもはる
かに高精度にモデル化する大規模な流体計算（Wall-Modeled LES）をエンジンルーム・タイヤハウス
内の流れも含めた自動車まわりの流れの解析、および、軸封部や戻り流路などの細隙部内の流れも含めた流体機械の内部流れの解析に適用し、このような大規模数値計算による流れの予測技術は、自動車の風洞試験（空力性能・騒音）や流体機械の性能試験（一般性能・吸込み性能）を完全に代替えし得る計算精度を有することを証明する。また、この予測技術を用いて、自動車の空力・騒音開発で問題となる、走行・操縦安定性に対する空力現象の寄与や実走行状態における空力音の発生機構、遠心圧縮機の運転範囲を拡大する上で重要となるサージの発生機構など、製品開発上、重要な現象でありながら従来は経験的に扱われていた複雑な流体現象を解明する。

２．令和３年度（報告年度）の実施内容
２－１．当該年度（令和３年度）の事業実施計画
エネルギー産業の心臓部となる「ターボ機械」、および、輸送産業の中核となる「自動車」を対象として、「富岳」あるいは「富岳」の時代における HPC を利用することにより、ものづくりの在り方を抜本的に変革できることを証明するために、5 つの実証研究テーマを設定し、大学等の研究者と民間企業の技術者・研究者が強力に連携して、本課題を実施する。以下に具体的な事業内容について記す。
(1) 数値曳航水槽の実現と省エネデバイスによる推進効率の向上
本実証研究テーマでは、「富岳」とポスト「京」重点課題⑧で研究開発した LES 解析ソフトウェア
である FrontFlow/blue（FFB）を用いて、300 億要素から最大 1,200 億要素の Wall-Resolved LES
計算により、船のまわりの乱流境界層を完全に解像する予測を実現し、曳航水槽試験を完全に数値
シミュレーションに代替えし得ることを実証するとともに、船の推進効率の向上のための省エネデ
バイスの動作メカニズムを解明し、さらなる高効率な省エネデバイスの開発に貢献する。
令和 3 年度は、「富岳」において、船を対象とした Wall-Modeled LES による数値シミュレーショ
ンを実施する。
(2) 細隙部を含めた多段遠心ポンプの内部流れの Wall-Resolved LES
本実証研究テーマでは、「富岳」を利用して、30～100 μm の最小渦スケールまで計算格子により直接解析する、Wall-Resolved LES を実施し、内部流れやその結果として決まる水力性能（全揚程・水力トルク・水力効率）の完全な予測を実現するとともに、細隙部内部流れの挙動やそれが性
能や信頼性に与える影響を明らかにし、ポンプ設計の高度化に貢献する。さらに、上記の 1/100 程-1-度の計算格子を用いた Wall-Modeled LES も実施し、Wall-Resolved LES の結果と計算精度、計算コストなどを比較することにより、Wall-Modeled LES の実用化を図る。
令和 3 年度は、「富岳」において、細隙部を含めた多段遠心ポンプの内部流れの Wall-Resolved
LES を実施する。
（協力機関：国立大学法人九州大学）(3) 圧縮機サージの直接解析
本実証研究テーマでは、プラント、およびガスタービンなどの実機に用いられる遷音速圧縮機を対象として、圧縮機本体だけでなく、それが設置される配管系まで含めたシステム全体を計算領域とし、圧縮機羽根車の失速現象を再現できるほど短い時間刻みを設定するとともに、システム全体にわたる長周期の流体振動現象を捉え得るほど膨大な時間ステップ数にわたって DES（Detached
Eddy Simulation）による非定常三次元流動解析を「富岳」上で実施することにより、圧縮機サー
ジの初生（マイルドサージ）から、逆流を伴うディープサージに至る非定常流動メカニズムを解明し、圧縮機サージの予測技術を確立する。
令和 3 年度は、「富岳」において、圧縮機サージの直接解析を実施する。
（協力機関：国立大学法人九州大学、国立大学法人岩手大学）(4) リアルワールド自動車空力性能の予測
本実証研究テーマでは、電気自動車や自動走行車等の、次世代自動車の設計・開発に貢献するための HPC シミュレーション技術の構築と実証を行うことを目的として、実走行状態の自動車の空
力性能（空気抵抗、操安性、横風安全性等）の、車体形状再現性を数ミリまで高めた高解像度・長
時間シミュレーションによる評価を実現する。
令和 3 年度は、「富岳」において、リアルワールド自動車空力性能予測のための解析を実施する。
（協力機関：国立大学法人神戸大学、国立大学法人山梨大学、国立研究開発
法人理化学研究所）(5) リアルワールド自動車空力音予測
本実証研究テーマでは、電気自動車や自動走行車等の、次世代自動車の設計・開発に貢献するための HPC シミュレーション技術の構築と実証を行うことを目的として、ボンネット隙間やフロン
トグリルから発生する狭帯域音、ピラーやドアミラーから発生する広帯域音、空力・構造振動・音響連成による車内騒音を対象に、実走行時の空力音予測を行う。
令和 3 年度は、「富岳」において、リアルワールド自動車空力音予測のための解析を実施する。
（協力機関：国立大学法人豊橋技術科学大学、国立大学法人神戸大学、国立研究開発法人理化学研究所）-2-

(6) プロジェクトの総合的推進
プロジェクト全体の連携を密としつつ円滑に運営していくため、プロジェクトの推進や実証研究
テーマ間の連携のための会議等を適宜開催し、参画の協力機関・連携機関との連携・調整にあたる。
特に、プロジェクト全体の進捗状況を確認し、計画の合理化の検討等を行うなど、プロジェクトの効果的・効率的推進に資する取組みを実施する。また、プロジェクトで得られた成果については、
可能な限り積極的に公開して今後の展開に資するとともに、ものづくり産業での早期戦力化を支援する。
令和 3 年度は、令和 2 年度に引き続き、上記を順次実施しながら、本課題の有効な研究活動のための支援を実施し、研究開発全体を円滑に実施するための調整を行う。
２－２．実施内容（成果）(1) 数値曳航水槽の実現と省エネデバイスによる推進効率の向上(a) 研究成果の要旨
本実証研究テーマでは、「富岳」とポスト「京」重点課題⑧で研究開発した LES 解析ソフトウェア
である FrontFlow/blue（FFB）を用いて、300 億要素から最大 1,200 億要素の Wall-Resolved LES
計算により、船のまわりの乱流境界層を完全に解像する予測を実現し、曳航水槽試験を完全に数値
シミュレーションに代替えし得ることを実証するとともに、船の推進効率の向上のための省エネデ
バイスの動作メカニズムを解明し、さらなる高効率な省エネデバイスの開発に貢献する。
令和 3 年度は、「富岳」において、船を対象とした Wall-Modeled LES による予測精度が高いことを確認するとともに、数値曳航試験のための予備計算を実施し、今後の課題を明確にした。また、
省エネデバイスのメカニズム解明のための実証研究において、プロペラ後方のデバイスの影響でプロペラ後方の流場のみならず、プロペラ上方から前方の流場にまで影響を及ぼしていることが分かり、今後の分析方針を明確にすることが出来た。
In this research, we will show that towing tank tests can be completely replaced by numericalsimulations by resolving turbulence boundary layer developed on the hull surface, and we willcontribute to the development of even more efficient energy-saving devices for a hull by clarifying
their mechanism. To this end, we will use FrontFlow/blue (FFB), which has been developed forsupercomputer “Fugaku”, for performing Wall-Resolved LES with 30~120 billion computational
elements.
In this fiscal year, we confirmed that accurate predictions of hull resistance are possible by wallmodeled LES, and by performing test computations, we clarified issues for realizing numericaltowing tests. We also analyzed unsteady flow fields around the propeller in computations for selfpropulsion tests in detail, for clarifying mechanisms that energy-saving devices enhance the selfpropulsion performance. By this analysis, we found that not only the wake of the propeller but
also flow fields ahead and upstream of the propeller could be altered by the energy-saving devices,which is important findings for clarifying these mechanisms.

-3-

(b) 本研究の内容
1) FrontFlow/blue（FFB）の開発
1-1) FFB の並列効率向上のための調査の概要
前年度までの FrontFlow/blue（FFB）の開発では、メモリへロードするデータの局所化、ループ構造の改良等により、スーパーコンピュータ「富岳」におけるコアカーネルのメモリスループット
を 333~803 GB/sec まで上げ、ノード単体性能 179 GFLOPS（ピーク性能比 5.3 %）を達成した[1]。
並列性能については、立方体領域を直列に連結し、両隣の領域とのみ通信する、通信コストが低い
weak scale ベンチマークテストでは、「富岳」のほぼ全系を用いて 85 %の並列化効率を達成している。Figure 2.2.1-1 における「富岳」におけるメモリスループット、Fig. 2.2.1-2 に「富岳」における weak scale テストの結果を示す。
上記の簡易形状を対象とした weak scale ベンチマークテストでは、高い並列化効率を確認することができたが、実用計算では並列化効率が低下する課題が確認された。また、FFB は、独立した領域関連でデータを受け渡すオーバーセット計算機能をサポートする[2]。オーバーセット計算では、
通常の領域境界に関する通信（以外、隣接通信）の他にオーバーセット処理のための通信（以外、
オーバーセット通信）が必要になり、これに起因して並列性能がさらに低下する課題が確認された。
今年度は「富岳」における FFB の並列性能を向上させるため、Tofu インターコネクト D における実際の通信経路を考慮した通信コストと通信性能の関係を調査した。

Fig. 2.2.1-1 Single-node memory throughputs of core-kernels of FFB flow solver on Fugaku.

Fig. 2.2.1-2 Sustained performance of FFB flow solver in weak-scale benchmark tests with 2
million grids per a computational node of Fugaku.

-4-1-2) 通信経路を考慮した通信コスト
FFB の隣接通信の通信情報（どのプロセスがどのプロセスにどれだけの情報を送受信するか）は、領域分割するための中間ファイルである Domain Decomposition Description（DDD）ファイルに格納されているため、このファイルを参照することにより、実行前に通信情報を把握できる。
「富岳」の計算ノードは 4 個の Core Memory Group（CMG）で構成されているので、並列計算の際にはノードに 4 プロセスが割り当てられ、ノードは論理上の空間（以下、Node 空間と称する）に割り当てられる。ユーザーは Node 空間の種類（離散、メッシュ、トーラス）、次元および形状（各方向の大きさ）を指定できる。物理的には、ユーザーが指定した Node 空間を実現するよう、
12 ノード（48 CMG）で構成される Tofu の 1 つのグループ（Tofu 装置）に、ノードが割り当てられる。Tofu 装置は格子状に配置されており、遠くの Tofu 装置との通信は、通信経路上にある Tofu
装置を介して行われる。Tofu 装置は実行時に Tofu-D システムにより自動的に割り当てられる。割り当てられた Tofu 装置の集合をここでは Tofu 空間と称する。Figure 2.2.1-3 に Tofu 装置、（Tofu
装置内に 12 個ある）ノードおよび（ノード内に 4 個ある）CMG（Core Memory Group）の関係を示す。
「富岳」における通信コストを正確に把握するためには、上記した FFB の通信情報を用いて、Tofu 空間における通信経路を考慮して、割り当てられた Tofu 装置が実際に実行する通信の回数やデータ転送量を積算する必要がある。Node 空間から、Tofu 空間への割付については、ユーザーは制御できないが、実行時に問い合わせることにより、各ノードの Tofu 座標 (X, Y, Z, A, B, C) の情報を得ることができる。ここで (X, Y, Z) は Tofu 装置の位置に対応した座標値であり、(A, B, C) はTofu 装置あたり 12 個あるノードの Tofu 装置内の位置に対応した座標値である。
FFB の通信情報と Tofu 座標を照合することにより、全てのプロセス間通信に関して、Tofu 空間における経路を考慮した通信コストがわかり、それを積算することができる。具体的には、Tofu 装置
をつなぐ各方向の 12 本のリンクに対して、通信の回数と転送データ量（バイト数）を積算する。なお、Tofu 空間における異なる 2 点をつなぐ経路は複数あるが、Tofu-D システムでは、XYZ の優先度で経路が決定される。2 点間の座標差を、Nx、Ny、Nz とすると、通信データははじめに X 方向に Nx 移動し、次に Y 方向に Ny 移動し、最後に Z 方向に Nz 移動し、通信先の Tofu 装置に到達する。通信先の Tofu 装置に到達したデータは必要があれば当該 Tofu 装置内で bca 軸上を移動する。
ただし、通信時間という観点では、Tofu 装置内のデータの移動は Tofu 装置間のデータの移動に対して無視できるものと推定されるため、以下では Tofu 装置間のデータの移動のみに関して議論する。この場合、Nx+Ny+Nz 回の通信が発生するので、Nx+Ny+Nz を Tofu 距離と称する。Figure
2.2.1-4 に Tofu 距離の定義を示す。Figure 2.2.1-4 の Tofu 装置を Node 空間におけるノードと置き換えれば、同様に Node 距離を定義することができる。

-5-

Fig. 2.2.1-3 Relationship of Tofu-unit, compute node and core memory groups (CMG) in tofu
system.

Fig. 2.2.1-4 Definition of node distance and tofu distance.
1-3) 隣接通信(1-3-1) 1 段階分割と 2 段階分割
FFB による並列計算では、プロセス数が決まれば計算格子をそのプロセス数に領域分割する。
「富岳」による計算では、ノードに 4 個ある各 CMG に対しプロセスを割りあてるため、プロセス数はノード数の 4 倍になる。通常、全計算領域を一度にプロセス数に分割するが（1 段階分割）、一度、ノード数に分割し、分割された各領域を 4 分割することにより、領域番号が連続した 4 領域が同じノードに配置されるため（2 段階分割）、通常の 1 段階分割と比較し並列化効率が向上することが確認されている。Figure 2.2.1-5 に、Cavity flow を対象とし、ノードあたりの要素数を 200 万とした
場合の、「富岳」における weak scale ベンチマークテストの結果を示す。1 段落分割に対し、2 段落分割の性能は最大 8 %程度向上している。

-6-

Fig. 2.2.1-5 Weak scale benchmark test with 2 million elements per compute node of Fugaku.
(1-3-2) 通信コストと通信性能
立方体領域を計算領域とする Cavity flow を対象に、領域分割方法を変更し、通信コストと並列性能の関係を調査した。総要素数は 20 億であり、「富岳」の 1,000 ノードを用いて、ノード空間として、10×10×10 の 3D トーラス配置を指定した。領域分割は上記のとおり、はじめにノード数（=1,000）に分割し、その後に各領域を 4 分割した。2 段階分割の初めの分割として、20 億要素からなる計算メッシュを 1,000 分割する際に、以下の 3 とおりの方法をテストした。
◼ Case A
Metis ライブラリを用いて分割する。領域境界面積が最小化されるが、隣接する領域が、Node 空間において近くに配置される保証はない。
◼ Case B
計算領域を 10×10×10 のブロック状に分割する。
（以下、この分割をブロック分割と称す）。これにより、Node 空間における通信は隣接するノードに限定される。
◼ Case C
Metis ライブラリを用いて分割した後に、各領域の領域番号がブロック分割した場合に近くなるよう領域番号を並び替える。領域境界面積は Case A と同じであるが、ブロック分割に近づける並び替えをすることにより、隣接する領域が Node 空間において近くなることにより、通信コストの低減が期待できる。
Node 空間における通信コストを比較するため、各ノードの通信回数を、Node 距離毎に全ノードにわたって積算しプロットした結果を Fig. 2.2.1-6 に示す。Case B はノード空間における隣接ノー
ドとのみ通信するので最大 Node 距離は 3 であり（ノード座標が 3 方向とも異なる隣接するノードとの距離は 3 となる）
、Fig. 2.2.1-6 でも Node 距離が 4 以上の通信は存在しない。Case A、Case C
は全体の通信回数は Case B よりも小さいものの、Node 距離 4 以上の通信が多く存在する。Case B（ブロック分割）に近づくよう領域番号を並び変えた Case C については、最大の Node 距離が 6 であり、Node 距離 7 以上の通信が多数ある Case A よりも通信コストが低いことが期待される。
Table 2.2.1-1 に、3 ケースの時間ステップあたりの計算時間、通信時間および演算性能の比較を示す。各ケースの通信時間は、ノード間の通信処理がない、1 ノード計算（計算規模 200 万要素、以下、Case 0 と称す）の計算時間との差分により算出した。Case 0 の計算時間は、使用した流れソル
バーFFB ver. 67 の演算性能推定値 144 GFLOPS/node から算出した。一方、Case A、B および C-7-の演算性能は、測定した計算時間から算出した。Case A の 123.1 GFLOPS に対し、Case B は 137.9
GFLOPS であり、ブロック分割により Node 距離が短くなり、12 %の並列性能の向上を確認することができた。一方、Case C に関しては、Node 距離から推定する限り、Case A よりも性能向上することが期待されたが、若干悪化することがわかった。

Fig. 2.2.1-6 Comparisons of maximum node distance.
Table 2.2.1-1 Comparison of time for computation and communication, and sustained
performance.
ケース

ノード数

Case 0Case A
Case B
Case C

1
1,000
1,000
1,000

要素数
200 万
20 億
20 億
20 億

計算時間
[sec/step]
0.278
0.325
0.290
0.331

通信時間
[sec/step]
0.000
0.047
0.012
0.053

実効性能
[GFLOPS/node]
144.0
123.1
137.9
120.8

Case A と Case C については、Node 距離によっては並列性能の違いが説明できないことがわかったので、より直接的に通信コストを反映する Tofu 距離を 3 ケースで比較した。10×10×10 のノード
空間に割り当てられたプロセスは 5×4×5 の Tofu 空間に割り当てられる。Table 2.2.1-2 に割り当てられた Tofu 空間全体で最も通信量が多いリンク（以下、最混雑リンク）の通信量を示す。前述したとおり、隣接する Tofu 装置を結ぶリンクは 12 本あるので、Table 2.2.1-2 にはリンク毎に最混雑リ
ンクの通信量およびそれらの最大値がまとめられている。各方向の通信量を比較すると、X 方向の通信量が多くなっているが、これは、前述したとおり、Tofu-D システムにおける通信が、XYZ の順に通信が行われるためである。ケース間で通信量を比較すると、Case B の通信量が少なく、Case A、C の最大通信量には大きな差がないことがわかる。これは、Case B の実効効率が最良であり、Case A、C の実効性能に大きな差がないことと整合する。各ケースの最混雑リンクの通信量
はそれぞれ、5.44 MB （Case A）、0.48 MB（Case B）、4.67MB（Case C）であり、Case A、C の
最混雑リンクの通信量は Case B と比較して、およそ 10 倍程度多い。Table 2.2.1-2 に示す各ケースあたりの通信時間（sec/step）を比較すると、Case B に対し、Case A、C の通信時間はおよそ 4 倍になっている。最混雑リンクの通信量が 10 倍程度になると、通信時間は 4 倍程度になることより、これらの量は比例関係にあるわけではないが、最混雑リンクの通信量の増加が通信時間の要因とな-8-

っていることは明らかである。このことから、Tofu 空間において、各通信の Tofu 距離を減らし、最混雑リンクの通信量を減らすことが、高い並列性能を達成するために重要であることがわかった。
Table 2.2.1-2 Maximum size of communication data [MB] in the X, Y and Z directions in Tofu
space.

リンク 1
リンク 2
リンク 3
リンク 4
リンク 5
リンク 6
リンク 7
リンク 8
リンク 9
リンク 10
リンク 11
リンク 12
最大値

X 方向
5.10
4.61
4.78
5.44
4.28
5.15
4.37
5.31
4.71
5.25
4.57
4.86
5.44

Case AY 方向 Z 方向
4.10
1.89
4.08
1.71
4.92
1.84
4.40
1.79
2.94
1.34
2.93
0.95
4.93
1.65
4.69
1.56
4.01
1.86
4.26
1.86
2.31
1.38
2.48
1.27
4.93
1.89

X 方向
0.48
0.48
0.48
0.48
0.48
0.48
0.48
0.48
0.48
0.48
0.48
0.48
0.48

Case BY 方向 Z 方向
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46
0.45
0.46

X 方向
2.65
3.07
2.97
3.17
3.30
4.43
2.34
3.91
2.77
3.16
3.50
4.67
4.67

Case CY 方向 Z 方向
1.41
1.78
1.69
1.86
2.13
1.40
1.62
1.93
1.55
1.87
1.70
1.70
1.47
2.00
1.39
1.59
1.52
1.61
2.49
1.88
1.37
1.65
1.49
1.52
2.49
2.00

1-4) オーバーセット通信
オーバーセットを含む計算として、風車まわりの流れ計算を対象に通信コストの分析をした（Case D）。Table 2.2.1-3 に Case D のテスト条件および計算性能をまとめる。1-3)に示したオーバーセットを含まないケースの実効性能は「富岳」の 1,000 ノードを用いて、およそ 120~140
GFLOPS/node であったが、Case D では 15.8 GFLOPS/node まで実効性能が低下している。
Table 2.2.1-3 Computational conditions and sustained performance of Case D.
計算条件

計算性能

Case 名
解析対象
要素数
ノード数
領域数
Node 空間形状
計算時間 sec/step
実効性能
GFLOPS/node-9-Case D
風車まわりの流れ
6.3 億
1,536
6,144
3D, 8×12×6
0.520
15.8

Table 2.2.1-4 Maximum size of communication data [MB] in the X, Y and Z directions in Tofu
space in Case D.

リンク 1
リンク 2
リンク 3
リンク 4
リンク 5
リンク 6
リンク 7
リンク 8
リンク 9
リンク 10
リンク 11
リンク 12最大値

隣接通信
X 方向 Y 方向
2.33
2.01
1.99
2.80
3.52
2.26
1.80
1.98
3.57
2.84
1.57
1.34
1.89
2.38
2.15
2.45
3.93
2.87
2.29
1.78
4.49
2.76
2.05
1.43
4.49
2.87

Z 方向
1.64
2.34
1.66
2.39
2.01
2.15
2.11
2.89
2.86
2.39
3.02
2.23
3.02

節点オーバーセット通信
X 方向 Y 方向 Z 方向
0.84
0.35
0.01
0.76
0.29
0.01
0.68
0.37
0.01
0.45
0.32
0.01
1.11
0.49
0.02
0.99
0.57
0.01
0.58
0.50
0.02
0.58
0.36
0.02
0.83
0.52
0.01
1.38
0.96
0.02
0.66
0.37
0.02
0.71
0.60
0.01
1.38
0.96
0.02

1-3)と同様、最混雑リンクにおける通信量を隣接通信、オーバーセット通信それぞれについて算出した結果を Table 2.2.1-4 に示す。最混雑リンクの通信量はオーバーセット通信よりも隣接通信のほうが多いため、オーバーセット通信による実効性能の極端な低下の理由は最混雑リンクの通信量増
大では説明できないことがわかった。
隣接通信とオーバーセット通信の通信パターンの違いを確認するため、各プロセスの通信の Tofu
距離の最大値と積算値を Fig. 2.2.1-7 にプロットした。Tofu 距離の最大値は、どちらの通信においても 0~10 の間に分布している。積算値については、隣接通信では概ね数十程度であるのに対し、オーバーセット通信では大半のプロセスが 30 以下であるが特定のプロセスで極端に大きな値をとっていることがわかる。これは、オーバーセット通信では、一部の特定のプロセスが遠くの Tofu 装置と通信していることを示しており、このような状況が通信時間の増大につながっていると推察される。特に、遠くの Tofu 装置との通信においてはホップ数が多くなるが、「富岳」の各リンクのルー
タは自身の CPU のデータ転送に 50%のパケット（1 回の通信で転送されるデータであり、「富岳」のディフォルトのパケット長は 1.9 kB）を割り当て、残りの 50%のパケットを他の CPU のデータ転送に割り当てるため、ホップ数が多い通信はリンクにおいて割り当てられるパケットが少なくなり、極端に長い待ち時間が発生している可能性がある。
一部のプロセスにのみ遠くの Tofu 装置との通信が集中しないよう、オーバーセット通信にかかわるプロセスを Tofu 空間において近くに配置することにより、並列性能向上の向上が期待できる。具体的には、準備計算により、オーバーセット通信にかかわる領域を特定し、これらのノードをNode 空間において近く局在するよう領域番号を並び替える方法が考えられる。ここで、Tofu 空間
ではなく、Node 空間での局在とする理由は、ユーザーは Tofu 空間への割り当てを制御できないためである。

- 10 -

(a) Maximum value.

(b) Integrated value.
Fig. 2.2.1-7 Maximum and integrated Tofu-distance of each process.
2) LES 用壁面モデルの開発
本課題では実証研究テーマ 1「数値曳航水槽の開発と省エネデバイスのメカニズムの解明」および実証研究テーマ 2「細隙部を含めた多段遠心ポンプの内部流れの Wall-Resolved-LES」において、主として、乱流境界層中の生成スケールの乱流渦まで計算格子によって直接その挙動を予測する
LES 解析手法である、Wall-Resolved LES（以下、WR-LES）の実証研究を実施している。これまでの調査研究によって、WR-LES を適用することによって船体の推進抵抗や遠心ポンプの水力性能の定量的な予測が可能である見通しが得られている。一方、WR-LES は粘性長さスケールに比例する乱流境界層中の内層まで直接計算する手法であるため、多くの計算資源が必要となる。これに対して、Wall-Modeled LES（以下、WM-LES）は乱流境界層中の内層は直接計算せずに、壁面摩擦応力をモデル化して与える手法であり、WM-LES によって、LES 計算に必要となる計算資源を飛躍的に少なくすることが可能であり、その実用化にも大きな期待が集まっている。
令和 2 年度に実施した調査研究では、壁面摩擦応力を計算するための参照点の壁面からの距離を適切に設定することによって、船体表面に発達する乱流境界層や船体の摩擦抵抗・圧力抵抗などを定量的に予測できる見通しを得た。令和 3 年度は WM-LES の適用に関する信頼性を向上させるために、WR-LES により計算される渦構造と WM-LES により計算される渦構造を比較するとともに、これらの計算によって予測される平均流速分布や変動流速分布を定量的に比較した。Figure 2.2.1-8に WM-LES によって予測された、船体表面近傍の瞬時の流れ場における渦構造を示す。対象として船型は KVLCC2 であり、また、船長と船速によって定義されるレイノルズ数は 2.3×106 である。

- 11 -

Fig. 2.2.1-8 Snapshot of vortical structures in WM-LES of flow around KVLCC2. Reynolds
number is 2.3×106.
令和 2 年度の研究成果を参考にして、壁面から粘性長さスケールの 100 倍の位置に壁面摩擦応力を計算するための参照点を設置した。渦構造は速度勾配テンソルの第 2 不変量の等値面によって可視化した。参照とした WR-LES は約 40 億を用いた計算であるのに対して、WM-LES は約 5,000 万要素の計算結果であるにも係わらず、WR-LES の結果（図示せず）と定性的には同様な渦構造が計算されていることが確認された。この WM-LES では WR-LES の 4 倍の時間刻み幅を設定できることも考慮すると、WM-LES は WR-LES の 256 分の 1 の計算資源で計算できることを確認した。境
界層内の平均速度分布および変動速度分布を比較した結果を Fig. 2.2.1-9 に示す。

- 12 -

(a) Average velocity distribution.

(b) Fluctuation velocity distribution.
Fig. 2.2.1-9 Comparison between velocity distribution in WM-LES (left) and that in WR-LES(right) of flow around KVLCC2. Reynolds number is 2.3×106.
Fig. 2.2.1-9(a)は平均速度分布、Fig. 2.2.1-9(b)は変動速度分布を示し、また、それぞれの左図は
WM-LES の計算結果、右図は WR-LES の計算結果をそれぞれ示す。いずれの速度も船速によって無次元化されている。各図の横軸は船長で無次元化した壁面からの距離を表し、凡例は船尾からの距離を船長で無次元化した値である。境界層内層を計算しない計算手法である WM-LES によっても壁面摩擦応力を計算するための参照点を適切に設定すれば、摩擦抵抗や圧力抵抗だけでなく、
WR-LES とほぼ同じ境界層速度分布が得られることが確認された（Fig. 2.2.1-9(a)）。一方、WMLES による変動速度分布は WR-LES による変動速度よりも大きく、特に、壁面近傍の変動速度の最大値は 1.5 倍程度に増加している（Fig. 2.2.1-9(b)）。変動速度分布を示した図では縦軸のスケー
ルが異なっていることに注意されたい。これらのことから、WM-LES による実用計算に関して次のことが明らかとなった。摩擦抵抗やその船長方向の分布に関して、WM-LES は WR-LES とほぼ同

- 13 -

じ結果を与えているため（令和 2 年度の成果報告書において報告済み）、壁面摩擦応力に関しては正
しく与えられていることがわかる。壁面摩擦応力が正しく与えられていれば、粘性長さスケールの100 倍程度のスケールの渦しか計算していない WM-LES によっても境界層内の平均速度分布は正
しく計算される。ただし、壁面近傍の最大渦スケールが大きくなるため、変動速度は WR-LES よりも大きくなる。今後、乱流境界層中の乱流エネルギーの収支のバランスなどを計算することによって、この点に関してさらに調査することが望まれる。
3) 数値曳航水槽
本研究課題では、曳航水槽試験の代替技術として数値曳航水槽を確立するため、船体に発達する
乱流境界層を解像する LES 解析を実施することにより、曳航水槽試験の計測精度に匹敵する 1 %以内の予測精度で船体抵抗を予測できることを実証してきた[3]。数値曳航水槽の利点は、試験に匹敵する精度を有する船体抵抗の予測に加え、試験で取得することにできない詳細な流れ場データを取得できること、これらを詳細に分析することにより船体まわりの流れ場に関する新たな発見ができることにある。本年度は、船長 6.7747 m、主流速度および船長を基準とするレイノルズ数 5×106 の船体モデルの抵抗試験に対応した抵抗計算を実施し、船体表面の乱流構造を分析した。LES 解析に使用した計算メッシュの要素数は約 720 億であり、「富岳」の 48,144 ノードを用いて計算し、計算
時間は 1.547 sec/step であり、船長および主流速度で規格化した単位無次元時間を計算するのに34.4 時間を要した。
Figure 2.2.1-10 に船首近傍および船尾近傍の第 2 不変量の等値面を示す。等値面はヘリシティで色付けされている。特に、この船型では底面の比較的大きなスケールの縦渦が安定的に存在することがわかった（Fig. 2.2.1-10(a)）
。これらの縦渦はレイノルズ数が大きい場合は弱くなるが、レイノ
ルズ数が小さい場合は顕著であり船体抵抗にも影響を与えることが確認された。また、船尾には、既知のビルジ渦以外の喫水付近に大規模な縦渦が存在することが確認された。

- 14 -

(a) Near the bow.

(b) Near the stern.
Fig. 2.2.1-10 Iso-surface of second invariant colored by helicity near hull surface.
FFB は計算実行時にメモリ上で計算メッシュを細分化する Refine 機能[4]をサポートしており、計算メッシュの要素数が数億を超える大規模計算を実行する場合、通常、全てのデータは可視化せず
Refine 前のベースメッシュに定義される、間引かれたデータを可視化している。本年度は、720 億
要素のデータを間引かずに可視化することを試行した。このために、「富岳」の PRE-POST ノード（6 TB=128 GB×48 メモリ）を使用した。Figure 2.2.1-11 に 720 億要素データを間引かずに使用した、渦構造の可視化結果を示す。渦構造は、Fig. 2.2.1-10 と同様、ヘリシティで色付けされた第 2不変量の等値面を用いている。データを間引かずに可視化することによって、乱流境界層中の微細な渦構造の存在を観察することができる。

- 15 -

(a) Near the bow.

(b) Enlarged view of (a).

(c) Enlarged view of (b).
Fig. 2.2.1-11 Iso-surface of second invariant colored by helicity near hull surface visualized with72 billion grids data.
4) 省エネデバイス
省エネデバイスが船尾近傍の流れ場および自航要素にあたえる影響を調査するため、Fig. 2.2.1-12に示すとおり、付加物をつけない裸穀（Naked）、舵フィンのみを付けた場合、舵バルブのみを付けた場合、および、舵フィン・舵バルブを付けた場合の 4 種類の条件の自航計算を実施した。自航性能は、裸穀と比較し、舵フィンのみを付けた場合は悪化するが、舵バルブをつけた場合は改善した。舵バルブ、舵フィンの両方をつけた場合はさらに改善することがわかった。Figure 2.2.1-13 にプロペラ後流における主流方向速度の瞬時分布、Fig. 2.2.1-14 にプロペラ近傍の中央段目における主流方向速度の瞬時分布を示す。

- 16 -

Fig. 2-2-12 Configuration of energy-saving devices to a hull and their population performance.

Fig. 2-2-13 Distributions of instantaneous streamwise velocity at cross section in wake ofpropeller.

Fig. 2-2-14 Distributions of instantaneous streamwise velocity at center cross section.

- 17 -

Figure 2.2.1-13 と Fig. 2.2.1-14 の左下に示すとおり、舵フィンのみのケースではハブ渦の減衰が遅れることがわかり、このハブ渦が自航要素の悪化の原因である可能性がある。舵フィンがあるとハブ渦の発生する理由についてはわかっておらず今後、流れ場を詳細に分析する。一方、舵バルブをつけた場合はハブ渦の発生がなく、この場合は、舵フィンにより、自航要素が改善しており、ハブ渦の有無により、自航要素に対する舵フィンの影響が変わることがわかった。
舵バルブが与える船尾近傍の流れ場に与える影響として、プロペラ直上（Fig. 2.2.1-14 の黄色の四角で囲った領域）における逆流が強くなることが確認できる。これは舵バルブによってプロペラ周りの流れが斜流化した結果と推察され、このこととプロペラ性能向上、自航要素の改善は整合している。
上記のように、省エネデバイスのメカニズムの解明に関して、今年度は、舵フィンや舵バルブが流れ場に与える影響として、ハブ渦の発生やプロペラまわり流れの斜流化といった現象を抽出した。
[1] Kato, C., Yamade, Y., Nagano, K., Kumahata, K., Minami, K. and Nishikawa, T.: Towardrealization of numerical towing-tank tests by wall-resolved large eddy simulation based on
32 billion grid finite-element computation, Proceedings of the International Conference.
[2] Kato, C., Kaiho, M. and Manabe, A.: An overset finite-element large-eddy-simulation methodwith application to turbomachinery and aeroacoustics, Transactions of ASME, Journal ofApplied Mechanics, Vol. 70 (2003), pp. 32-43.
[3] T. Nishikawa, Y. Yamade, M. Sakuma, and C. Kato: Application of Fully-resolved Large Eddy
Simulation to KVLCC2 －Bare Hull Double Model at Model Ship Reynolds Number, Journalof the Japan Society of Naval Architects and Ocean Engineers, vol. 16, pp. 1-9. 2012.
[4] 吉村忍、徳永健一、杉本振一郎、奥田洋司、末光啓二、加藤千幸、山出吉伸、吉村英人：並列モデル細分化による大規模有限要素解析の効率化、日本計算工学会論文集（2013）、Paper No.
20130012。
(2) 細隙部を含めた多段遠心ポンプの内部流れの Wall-Resolved LES(a) 研究成果の要旨
本実証研究テーマでは、「富岳」を利用して、30～100 μm の最小渦スケールまで計算格子により直接解析する、Wall-Resolved LES を実施し、遠心ポンプの内部流れやその結果として決まる水力
性能（全揚程・水力トルク・水力効率）の完全な予測を実現するとともに、細隙部内部流れの挙動
やそれが性能や信頼性に与える影響を明らかにし、ポンプ設計の高度化に貢献する。さらに、上記の 1/100 程度の計算格子を用いた Wall-Modeled LES も実施し、Wall-Resolved LES の結果と計算
精度、計算コストなどを比較することにより、Wall-Modeled LES の実用化を図る。
令和 3 年度は、「富岳」において、細隙部を含めた多段遠心ポンプの内部流れの Wall-Resolved
LES を実施するため、解析格子を改善することによって、幅 0.17mm 程度の細隙部近傍で発生する計算の不安定性を解消し、最大 400 億要素規模の Wall-Resolved LES の準備が完了した。約 6.6 億
要素の 1 回リファイン計算と 53 億要素の 2 回リファインを行い（1 部の計算は途中経過）、試験結果と比較し予測精度を検証した。その結果、締切流量（流量零）近くの右上がり揚程特性を予測でき、LES はこの右上がり揚程特性の原因となる内部流れを正確に予測できていることを確認した。

- 18 -

しかしながら設計流量から大流量側ではリファインにより実験結果との差が拡大していることがわかった。これは、オーバーセット境界での圧力誤差が発生することが原因であることを確認し、オーバーセット境界での補間計算方法を修正することとした。
In this research, we will realize a highly accurate prediction of the internal flow of a centrifugal
pump and its performance, clarify flow mechanism in the leakage passages and their influenceon the performance and reliability of the pump, and finally contribute to advancing pump design.
To this end, we perform Wall-Resolved LES of internal flows of pump, where dynamics of thesmall vortices with the scale of 30~100 micro-meters are directly computed by using “Fugaku”.
We will also perform wall-modeled LES, which requires 1/100 times smaller number of thecomputational elements than those for Wall-Resolved LES. We will realize practical use of wallmodeled LES by comparing the accuracy and the cost of computations by Wall-Model LES and
Wall-Resolved LES.
The preparation of the Wall-Resolved LES with up to 40 billion elements has almost been
completed by the research activity in this fiscal year. An instability of the calculated flow fieldnear the small clearance of 0.17 mm in the balance drum domain was suppressed by investigating
the effect of the mesh quality on the computed flow fields and thereby by improving the mesh
quality. The once-refined calculations with about 660 million elements and the twice-refined
calculations with about 5.3 billion elements were performed, the latter is still underway, and theprediction accuracy was verified by the comparisons with the test data. As a result, it was
confirmed that LES could predict the unstable head characteristic near the shut-off flow condition.
However, at the designed and larger flow rates, the total head predicted by the LES was higher
than the measured data, and the over-prediction of the total pump head was larger for the twicerefined computations than for the once-refined computations. One of the main reason of thisproblem was the unphysical pressure increase at the overset interface. To solve this problem it
was decided to modify the interpolation calculation method.
(b) 本研究の内容
1) 遠心ポンプ内部流れ解析の計算モデル
スーパーコンピュータ「富岳」用にチューニングした FrontFlow/blue（FFB）流れ解析ソルバーを用いて、遠心ポンプ内部流れの LES 解析を実施した。FFB の性能向上に関する研究成果に関しては実証研究テーマ 1 に記述したとおりである。Figure 2.2.2-1 に遠心ポンプ内部流れ解析の計算モデルを示す。Figure 2.2.2-1 に示すとおり、この LES 解析では軸スラストを評価するため、入口配管、羽根車、ディフューザ、出口配管に加え、front/back shroud、balance drum、return pipe
といった細隙部を考慮した計算モデルを用いている。前年度まで用いていた計算モデルは上記した
流路ごとに用意された 10 個の独立の計算領域から構成されており、これらをオーバーセット法により連結させていたが、今年度の研究では、計算の安定性を向上させるため独立計算領域を 4 つにまとめた計算モデルを用いた。ベースとした計算メッシュの要素数はおよそ 8 千万である。FFB は計- 19 -

算実行時に計算格子を細分化する機能（refine 機能）をサポートしており、refine 機能を用いることによって、1 回 refine で約 6.6 億要素、2 回 refine で約 53 億要素、3 回 refine で約 421 億要素の計算メッシュを用いた計算が実行できる。本研究では 421 億要素の計算メッシュによって、細隙部も含め乱流境界層を解像する WR-LES を実施することを計画としている。今年度は 8 千万要素、
6.6 億要素および 53 億要素の計算を実行したが、53 億要素の計算では低流量側（20 %設計流量）での解析では計算が発散する問題が生じており、長時間の計算は実施できておらず原因を調査中である。ここでは、8 千万要素、6.6 億要素の計算結果、ならびに 53 億要素の計算結果の一部を報告する。以降、8 千万要素の計算を R0、6.6 億要素の計算を R1、53 億要素の計算を R2 と称する。
R0 および R1 は「富岳」の 2,052 ノードを用い、羽根車 1 回転あたり各々30 分、2 時間で計算した。R2 の計算時間は、「富岳」の 4,104 ノードを用いて、羽根車 1 回転あたりおよそ 15 時間を要した。

Fig. 2.2.2-1 Computational model of LES of pump internal-flow (top: overview, bottom: impeller).

- 20 -

Balance drum 部や front shroud 部にはポンプの漏れ流れを少なく維持するため、半径方向に 0.17mm 程度の小さい隙間部分（クリアランス部とよぶ）が設置されている。独立した計算領域を 4 つにまとめた後の計算において、計算要素の軸方向幅Δz と半径方向幅Δr のアスペクト比Δz/Δr が過大である場合には Balance drum のクリアランス部からその下流の広い領域に流出した後、非現実的な流れが発生することが分かった。Figure 2.2.2-2 に、Balance drum 部の形状のみを取り出した計算によって、アスペクト比の影響を調査した結果を示す。アスペクト比が 200 と大きい場合にはクリアランス部から流出した流れは大部分で逆流（コンターの赤い部分）していることがわかる。
これに対しアスペクト比Δz/Δr＝5.0 の場合にはクリアランス部から広い領域に流出する紫色の噴流部が観察され、正常な流れ状態になっていることがわかる。以下に示すポンプ全体の計算結果はBalance drum 部のアスペクト比を適切に修正した計算メッシュを使用して求めたものである。

Fig. 2.2.2-2 Effect of aspect ratio of element in the balance piston region.
2) ポンプ性能の計算結果
ポンプ性能として、全揚程（上）とポンプパワー（下）の流量特性を比較した結果を Figure
2.2.2-3 に示す。本ポンプでは、締切流量（流量零）近くの Q/QD<0.5 で、全揚程の右上がり特性が発生するが、LES はこの右上がり特性を捉えることができており、右上がり特性の原因となる内部流れを正確に予測していると考えられる。一方、LES は Q/QD>0.4 で全揚程（上）を過大評価しており、流量が大きいほど LES と実験結果の差が大きい。また、R0 より R1、R2 の方が実験との差- 21 -

が大きくなっている。ポンプパワー（下）の流量特性については、LES は全体的にやや過大評価しているものの、概ね再現できている。

Fig. 2.2.2-3 Comparisons of pump performance (top: total pressure, bottom: pump power).
Figure 2.2.2-4 に軸スラストの比較を示す。R0 では全体的に軸スラストを過小評価しているが、計算メッシュをリファインした R1、R2 では軸スラストの予測結果は実験値に近づく傾向にあり、
軸スラストの予測精度は向上していることが確認できる。この理由は、後述する細隙部の流れ場の予測精度向上に起因していると考えられる。

- 22 -

Fig. 2.2.2-4 Comparisons of axial thrust.
以下、LES が全揚程を過大評価する原因について調査した結果を述べる。Figure 2.2.2-5(a) は
Q/QD=1.0 でのポンプ入口から出口までの静圧の変化を静圧係数で示したものである。ここで PMP
INL、IMP INL、IMP OUT、DIF INL、DIF OUT、PMP OUT はそれぞれポンプ入口、羽根車入口、羽根車出口、ディフューザ入口、ディフューザ出口、ポンプ出口を表す。これを見ると R1 ではディフューザ出口よりもポンプ出口の静圧が大きくなっていることが分かる。Figure 2.2.2-5(b)はディフューザ領域と吐出ケーシング領域間に設置されたオーバーセット境界での静圧差を示している。計算メッシュを 1 段階リファインした R1 では境界部で静圧が非物理的に増加しており、オーバーセット境界での流れ場の補間計算に問題があることがわかった。本問題については補間方法
を圧力による補間方法に変更することで解決できる見込みである。

- 23 -

Fig. 2.2.2-5 Comparison of pump head and static pressure change at the overset interface
between the vaned diffuser and discharge casing domains.
上記のオーバーセット境界における静圧増加量（overset rise）とそれを差し引いた場合の全揚程係数（corrected）を試験結果とともに Fig. 2.2.2-5(c)に示す。オーバーセット境界での静圧増加量を差し引いても全揚程の解析値は Q/QD>0.6 では実験値よりも大きい。
この要因をさらに分析するため、Table 2.2.2-1 に羽根車のオイラーヘッドと各部の全圧損失を示す。全損失に対する各部の損失の割合は 20 %流量以外ではディフューザ部の損失割合が大きい。羽根車の損失は 100%流量で最小となっているが、ディフューザ損失は 60 %流量付近で最小値を示し、ディフューザベーンへの負の入射角が過大となることにより大流量側での損失増加が大きい。
R0 と R1 の比較により、計算メッシュのリファインにより羽根車の損失は増加しているが一方でディフューザ部の損失は減少していることがわかる。本研究対象である低比速度の遠心ポンプでは、羽根車やディフューザ部に発生する境界層や低エネルギー流体の蓄積による流路のブロッケージが性能に及ぼす影響が大きく、計算メッシュの解像度が十分でない場合には羽根車のオイラーヘッドを過大評価し、静止流路部の損失を過小評価すると考えられる。今後オーバーセット境界での圧力補間方法の改善と、さらなるリファイン計算によりその改善効果を確認する。

- 24 -

Table 2.2.2-1 Computed Euler’s head and decomposed losses.

ERH: Euler’s head, IMP: impeller, DIF: diffuser, DIS: discharge, OHR: other, TTL: total

3) 計算メッシュのリファインによるポンプ内部流れの変化
Figure 2.2.2-6 に 100 %流量における羽根車とディフューザ幅方向中央部での渦度分布を示す。計
算モデルは R0（左）と R2（右）である。R0，R2 の両ケースで羽根車内では負圧面に沿って渦度が発達していることがわかる。計算メッシュのリファインを施した R2 では羽根出口に向かって流路ピッチ中央付近の渦度が増加するが、圧力面側では渦度が小さいままであり、負圧面と比べ乱流境界層の明確な発達が見られない。これにより羽根車内のブロッケージを過小評価していることが考えられる。この問題の対策として、羽根車に流入する流れに乱れを与えることによって、羽根車
動翼の圧力面の境界層遷移を促進させることを検討している。

Fig. 2.2.2-6 Vorticity distributions of the pump mid-span plane.

- 25 -

4) 羽根車前面間隙部とライナーリング部の流れ羽根車に作用する軸スラスト力を適切に予測するためには、羽根車とケーシング間の間隙部（Front chamber、Rear chamber）と半径方向隙間が約 0.17 mm のライナーリング部の流れを適切に予測する必要がある。本研究では、1)に示したとおり、ベースメッシュを 3 回リファインした
421 億要素の計算メッシュでライナーリング隙間などの細隙部の乱流境界層も解像できるようにベースメッシュがデザインされている。
Figure 2.2.2-7 に R1 の場合の羽根車前面間隙部の周方向速度分布（上）と、半径方向静圧分布（下）を設計流量（Q/QD=1.0、左）と 20 %流量（Q/QD=0.2、右）の場合について示す。Fig. 2.2.27(a)に示した間隙部内の周方向速度分布はおおよそその半径位置における羽根車回転速度の 50 %となることが知られている。設計流量での周方向速度分布は外周側で最も大きい分布となっているが、20 %流量ではディフューザ側から間隙部に流入する逆流の影響により外周側で周方向速度が大きい部分（赤部）小さい部分（青部）が交互に発生し、間隙部の周方向速度を減速する作用を及ぼす。
Fig. 2.2.2-7(b)に示した半径方向静圧分布の比較では、LES の計算結果は実験値よりも小さい値を示している。半径方向静圧分布は間隙部に持ち込まれる周方向速度（角運動量）と逆流などの効果、ライナーリング部での漏れ流量の予測精度に関係するが、計算メッシュのリファインを進めること
によって予測精度が向上できると考えている。

Fig. 2.2.2-7 Circumferential velocity contours (top) and pressure distributions (bottom) in thefront gap of the impeller.

- 26 -

Figure 2.2.2-8 にライナーリング部における隙間幅方向の軸方向速度、周方向速度分布を R1 および R2 の場合について示す。R1 では隙間方向に要素数が 10 と少ないため隙間部内の乱流は計算できず、z=−0.805、−0.835 では放物線的な軸方向速度分布となっている。R2 の場合にはこれがやや改善し、同位置での軸方向速度分布を見ると幅中央部での平坦部分が増えていることが分かる。今後、3 回リファインした計算（R3）を実施する予定であるが、ライナーリング部の流れのレイノルズ数は小さいため、R3 においてもライナーリング部の境界層は完全な乱流には遷移せずに、遷移境界層に留まる可能性もある。

Fig. 2.2.2-8 Velocity distributions in the liner ring clearance.
(3) 圧縮機サージの直接解析(a) 研究成果の要旨
本実証研究テーマでは、プラント、およびガスタービンなどの実機に用いられる遷音速圧縮機を対象として、圧縮機本体だけでなく、それが設置される配管系まで含めたシステム全体を計算領域とし、圧縮機羽根車の失速現象を再現できるほど短い時間刻みを設定するとともに、システム全体にわたる長周期の流体振動現象を捉え得るほど膨大な時間ステップ数にわたって DES（Detached
Eddy Simulation）による非定常三次元流動解析を「富岳」上で実施することにより、圧縮機サー
ジの初生（マイルドサージ）から、逆流を伴うディープサージに至る非定常流動メカニズムを解明し、圧縮機サージの予測技術を確立する。
令和 3 年度は、前年度に引き続き、ターボ機械用 DES 解析ソフトのチューニングを実施し、最終的にオリジナルソフトと比較して約 90 倍の高速化を実現した。本アプリケーションを遷音速遠心圧
縮機サージの直接解析に適用した結果、逆流を伴わずに流量および圧力が長周期に変動する現象が観測され、シミュレーションによってマイルドサージ現象を再現することに世界で初めて成功した。
これにより、本シミュレーション技術によってマイルドサージの発生を予測できることを実証した。

- 27 -

さらに、解析された流れ場を詳細に分析することによって、マイルドサージ発生時の圧縮機内部の詳細な非定常流動現象を明らかにした。
In this study, we will develop the unsteady three-dimensional flow analysis technology of thesurge phenomenon for transonic compressors used in plants and gas turbines to elucidate theunsteady flow mechanism that leads the inception of the surge (the mild surge), and eventually
results in the deep surge, which is accompanied by reverse flows. This analysis is based on the
DES (detached eddy simulation) and considers not only the compressor but also the entire system,
where the compressor is installed. The surge simulations require a huge number of time steps tocapture both of short-cycle stall phenomena in the compressor and a long-period surge
phenomenon throughout the system, and we aim to achieve such simulations by using “Fugaku”.
In this fiscal year, we tuned the DES analysis software for turbomachinery and achieved a
speed-up of 90 times, compared with the original version. Subsequently, the software was applied
to the surge simulation of a transonic centrifugal compressor. Consequently, long-period
oscillations of the flow rate and pressure were observed when the flow rate was still positive. Thisphenomenon had been known as the mild surge by the experiments, but was never be reproduced
in direct simulations. We demonstrated that the present simulation method can predict theoccurrence of the mild surge in the first time in the world. Furthermore, we elucidated detailed
unsteady flow phenomena in the compressor when the mild surge occurs.
(b) 本研究の内容
1) 計算モデル
車両用ターボ過給機に用いられる遷音速遠心圧縮機を解析対象とした。Figure 2.2.3-1 に供試遠心圧縮機の概観を示す。本遠心圧縮機は、オープンタイプのインペラ、羽根なしディフューザ、およびスクロールから構成されている。Table 2.2.3-1 に主要諸元を示す。供試圧縮機インペラの翼枚数
はフルブレードが 6 枚、スプリッタブレードが 6 枚の計 12 枚である。インペラとケーシングの間には 0.28 mm の翼端隙間を有する。
供試遠心圧縮機のサージ試験装置を Fig. 2.2.3-2 に示す。供試遠心圧縮機は、ターボ過給機のタービン側にスクリュ圧縮機から圧縮空気を供給することにより駆動される。タービン側に供給する圧縮空気の流量を調整することによって、圧縮機の回転数を最大 160,000 rpm までの範囲で任意に設定可能である。本研究では、回転数 140,000 rpm の条件で計算を行った。圧縮機流量は圧縮機の吐出側に配置したプレナムの下流に位置する流量調整弁により調整される。本試験装置では、供試圧縮機の吐出側配管系のプレナム容積を 3 とおりに変化させることができるが、本研究では最小容積の条件を対象とした。
上記のサージ試験装置を模擬した解析モデルを Fig. 2.2.3-3 に示す。サージの直接解析では圧縮機本体だけでなく、圧縮機が組み込まれたシステム全体を考慮する必要がある。したがって、本解析モデルでは、供試圧縮機に加え、その上流側および下流側の配管が考慮されている。ディープサー
ジ発生時には順流と逆流が周期的に繰り返されることから、逆流時に上流側配管の入口から吐き出- 28 -

された高温流体が順流時に再度吸い込まれる現象が生じる。この現象を再現するために、 Fig.
2.2.3-3 に示すように、入口配管の上流には試験装置の外部空間を模擬した領域を設置している。外部領域は半球状の形状とし、その半径は圧縮機外径の 20 倍とした。また、サージ初生の流動現象を正確に再現するため、実際の試験と同様に、計算中に流量を漸減させる。流量制御は、実際の流量調節弁を模擬して行われる。流量調節弁の流路形状を計算格子によって直接再現し、計算途中に開度を変化させる方法は計算コストが非常に大きいため、流量調節弁については抵抗係数としてモデル化される。
本解析モデルについて計算格子を作成した。格子生成には市販ソフト Pointwise を使用し、格子は
すべて 6 面体の構造格子とした。また、粘性長さスケールで無次元化した壁面上の最小格子幅 y+はy+<1 を満足するように設定した。Table 2.2.3-2 に作成した計算格子について各領域のセル数の内訳を示す。サージ初生における圧縮機内部の小スケールの流動現象を捉えるため、解析領域全体の総セル数の約 4 割を圧縮機部に割り当てている。Figure 2.2.3-4 にインペラ部の格子図を示す。フルブレードおよびスプリッタブレードのコード方向にそれぞれ 320 セルおよび 200 セルを配置し、フルブレードとスプリッタブレードの各翼間の周方向には 120 セルを設けた。また、ハブからケーシングまでのスパン方向には 180 セルを設定し、そのうち翼端すき間には 40 セルを割り当てた。イ
ンペラ部が約 1.0 億セル、ディフューザ部が約 1.0 億セル、スクロール部が約 1.8 億セルであり、圧縮機本体の合計は約 3.8 億セルである。また、入口配管は約 2.3 億セル、出口配管は約 1,100 万セル、外部領域は約 2.9 億セルとなった。インペラ前縁部の逆流を捕捉するため、下流側よりも上流側配管の格子解像度が比較的高く設定されている。解析領域全体の総セル数は約 9 億である。本研究で実施される DES では、乱れの長さスケールに応じて、RANS または LES で計算される領域が決定される。前年度までに、RANS/LES 領域を可視化することにより、壁面近傍の境界層に対応した領域のみが RANS 領域であり、その他の主流領域は LES として計算されていることを確認している。以上のように、本計算格子は十分な格子解像度を有しており、圧縮機内に生じる複雑な大規模はく離流れを高精度に解析可能である。

Fig. 2.2.3-1 Test centrifugal compressor.

- 29 -

Table 2.2.3-1 Specifications of the test compressor.
Number of full blades

6

Number of splitter blades

6

Impeller inlet diameter [mm]39.9

Impeller outlet diameter [mm]52.5

Tip clearance [mm]0.28

Diffuser width [mm]3.42

Drive shaft speed [rpm]140,000

Fig. 2.2.3-2 Schematic diagram of the experimental equipment for compressor surge.

- 30 -

外部領域

20D

圧縮機
入口配管

出口配管

Fig. 2.2.3-3 Computational domain.
Table 2.2.3-2 Number of cells in computational grid.
Impeller
Diffuser
ScrollInlet duct
Outlet duct
External area
Total

100,608,000
100,224,000
178,791,000
228,960,000
11,491,200
288,000,000
908,074,200

(a)

180

320
200
Rotation

116

116

Fig. 2.2.3-4 Computational grid in the impeller.

- 31 -

2) 並列計算の方法と効果（性能）
本研究で実施する圧縮機サージの直接計算では、サージ発生時における圧縮機内部の非定常流動まで詳細に再現する。そのため、圧縮機羽根車の失速現象を再現できる程度に、短い時間刻みを設定して計算を実施する。その結果、サージについてシステム全体にわたる長周期の流体振動現象を捉えるためには、膨大な時間ステップ数にわたって非定常三次元流動解析を実施することになる。
サージ 10 周期を計算する場合、7,500 step/rev×200 rev/cycle×10 cycle より、15,000,000 ステップの計算が必要であると見積もられる。上式において、step、rev、cycle はそれぞれ時間ステップ数、羽根車回転数、および、サージのサイクル数をそれぞれ表す。仮に、アプリケーションの時間
ステップあたりの計算時間が 0.1 sec/step であれば、総計算時間は 17 日となり、現実的な時間内で計算を完了できることになる。そこで、0.1 sec/step を目標に、アプリケーションであるターボ機械
用 DES 圧縮性流れ解析のチューニングを前年度より継続して実施した。なお、アプリケーションのチューニングにあたっては一般財団法人高度情報科学研究機構（RIST）の高度化支援を利用した。
前年度までに、ファイル出力のプロセス毎に分散化、圧縮性流体解析の主要部である高コスト箇所についての OpenMP によるスレッド並列化、sliding mesh 処理に使用している Akima 補間サブ
ルーチンについての OpenMP スレッド並列化および call 回数の削減などの対策を実施しており、オリジナルのソルバーと比較して、約 60 倍の高速化が達成されている。今年度は、さらに sliding
mesh ルーチン内の集団通信性能の分析および最適化、特にランクマッピング等の設定による性能改善について検討を行った。サブコミュニケータによる通信処理に mpi_allreduce を使用すること
で、通信時間が 1/2.7 に短縮された。一方、ノード形状の変更やランクマップの指定では高速化の効果は得られないことがわかった。通信データサイズ 40 MB 弱に要する通信時間について、IntelMPI Benchmarks（IBM）の計測結果と比較した結果、概ね IBM 結果と相応の通信性能が得られ、最適化はほぼ限界に達していることが確認された。
チューニングを施したコードについて、本計算と同様の条件で実行時間を計測し、性能を評価した。評価結果を Table 2.2.3-3 に示す。Asis に比べて、最新バージョンの Tune 3 は約 90 倍の高速化が実現されたことが確認できる。実行時間は 0.34 sec/step であり、目標とした 0.1 sec/step には至っていないが、Asis から大幅に実行時間は改善されており概ね目標は達成されている。
Table 2.2.3-3 Performance evaluation of tuning codes with the surge simulation.
実行時間[sec/step]

Asis との速度比

Asis

31.53

Tune 1 (2020/8)

2.91

10.8

Tune 2 (2021/3)

0.53

59.5

Tune 2 (2022/1)
Tune 3 (2022/1)

0.44
0.34

71.7
92.7

- 32 -

3) マイルドサージの解析
チューニングされたコードを用いて、初めに、逆流を伴わないサージ形態であるマイルドサージの計算を試みた。マイルドサージを再現するために、Fig. 2.2.3-5 の性能線図に示すように、出口境界条件に導入したバルブの抵抗係数 Kt を調整し、高流量側から流量を徐々に低下させていった。
Figure 2.2.3-6 に Kt=1.8 および Kt=2.0 の条件の流量の計算履歴を示す。Kt=1.8 の条件では流量変動の振幅は±10 %以下であったが、Kt=2.0 の条件では振幅が±50 %程度の長周期の流量変動が観測された。Figure 2.2.3-7 には、出口配管の圧力および入口配管の流量について Kt=2.0 における計算結果とマイルドサージ発生時の実験結果との比較を示す。計算結果は実験結果に比べて、流量および圧力ともに振幅が大きく、周波数は若干小さいが、定性的にはマイルドサージ現象を再現できていると考えられる。定量的な予測精度については、バルブモデルにおいて計算を安定化されるために導入した緩和係数などが影響している可能性があり、今後の課題であることが確認された。
2.4

Exp.
Sim.

Total pressure ratio [-]

0.2
2.2

Kt=1.8Kt=2.0Kt=1.26

2

1.8
1.6
1.4
1.2

0

0.05

0.1

0.15

0.05

Mass flow rate [kg/s]

0.2

Fig. 2.2.3-5 Performance curves of the compressor.
200

配管入口
圧縮機入口
配管出口

150

Mass flow rate [%]Mass flow rate [%]

200

100

50

0

配管入口
圧縮機入口
配管出口

150

100

50

0
0

1000

2000

3000

0

1000

t* [-]

2000

t* [-](a) Kt = 1.8.

(b) Kt = 2.0.

Fig. 2.2.3-6 Time histories of the flow rate.

- 33 -

3000

240

1.6
1.4

230
Pressure [kPa]

Mass flow [-]

1.2
1

0.8
0.6
0.4

DES

0.2
0

210
200
DES

190

Exp.

0

220

Exp.

180

0.05

0.1
Time [s]

0.15

0.2

0

(a) Mass flow rate in the inlet pipe.

0.05

0.1
Time [s]

0.15

0.2

(b) Pressure in the outlet pipe.

Fig. 2.2.3-7 Comparison of simulation and experimental results.
Kt=2.0 の計算結果を用いて、マイルドサージ発生時の圧縮機内の流れ場を調べた。ただし、サージ 1 周期はインペラが 160 回転する時間に相当し、これらの全てを可視化するのはデータ量が膨大となり困難であるため、Fig. 2.2.3-8 に示すように、サージ 1 周期中の 8 つの代表的な流量フェーズ
について可視化した。

Mass flow rate [%]

200

B C

150

D
E

A

F

100

H

50

0
3500

G
4000

4500

5000

5500

t* [-]

Fig. 2.2.3-8 Representative points for flow visualization.
Figure 2.2.3-9 に各流量フェーズについて圧縮機 3 回転分の時間平均および周方向平均した子午面速度分布を示す。いずれの流量フェーズにおいてもシュラウド面付近にインペラとディフューザ
にわたる逆流領域が確認される。この逆流領域はサージによる流量変化に伴って、拡大・縮小を繰り返すことが明らかになった。すなわち、高流量フェーズ（B-D）には逆流領域は縮小し、低流量フェーズ（F-H）には逆流領域がスパン方向に拡大する。Figure 2.2.3-10 には 90 %スパン断面における相対マッハ数分布を示す。高流量フェーズでは低速領域は翼間後半部のみに存在するが、低流量フェーズでは低速領域がインペラ前縁の上流まで拡大することがわかる。Figure 2.2.3-11 に各流量点についてフルブレードの 90 %スパン位置における翼負荷のコード方向分布を示す。流量変動に- 34 -

伴って、フルブレード翼端付近の翼負荷が変動することがわかる。低流量フェーズでは翼先端前半部で失速していることが明らかになった。

-0.4

Shroud-0.2

Full blade
Splitter blade

Hub

Diffuser(a) point A (100 %).

(b) point B (120 %).

(c) point C (135 %)(d) point D (120 %).

(e) point E (100 %).

(f) point F (75 %).

(g) point G (65 %).

(h) point H (75 %).

Fig. 2.2.3-9 Meridional velocity contours.

- 35 -

1.2-0.0(a) point A (100 %).

(b) point C (135 %).

(c) point E (100 %).

(d) point G (65 %).

Fig. 2.2.3-10 Relative Mach number contours in the cross-section located at the 90 % span.

Fig. 2.2.3-11 Chordwise load distributions of the full blade at 90 % span.
A 点について 4 つの周方向位置におけるディフューザ内の半径方向速度分布とスクロール内の周
方向速度分布を Fig. 2.2.3-12 に示す。ディフューザ内には、シュラウド付近に逆流領域が認められる。一方、スクロール内には，舌部の上流に位置する θ=0°および 270°の断面において部分的に逆流している領域が観察される。Figure 2.2.3-13 に A 点（中間流量），C 点（最小流量）および，G 点（最大流量）について θ=0°の断面におけるディフューザ内の半径方向速度分布とスクロール内の周- 36 -

方向速度分布を比較している。上記のスクロース内の逆流領域は高流量フェーズでは縮小し、低流量フェーズでは拡大することが明らかになった。

Cθ

（+）

θ = 180 °

周方向速度

Cθ[-]
Cr[-]
-0.2
-0.0

（-）
θ = 90°θ = 270 °半径方向速度

-0.5

-0.0

逆流領域
θ=0°

θ = 90°θ = 180 °θ = 270 °θ=0°Fig. 2.2.3-12 Radial velocity distributions in the diffuser and circumferential velocity
distributions in the scroll for phase A.

(a) point A (100 %).

(b) point C (135 %).

(c) point G (65 %).

Fig. 2.2.3-13 Radial velocity distributions in the diffuser and circumferential velocity
distributions in the cross-section located at θ=0°.

- 37 -

4) ディープサージの解析
次に、逆流を伴うディープサージの解析を実施した。マイルドサージの計算ではバルブの抵抗係
数 Kt を 2.0 に設定していたが、Fig. 2.2.3-14 に示すように、これを 3.0、4.0、5.0 と徐々に大きくすることによって流量を低下させていった。その結果、Kt=5.0 の条件における流量履歴（Fig.
2.2.3-15）において、配管入口および圧縮機入口断面の流量がマイナスになる瞬間が観測され、ディープサージが発生したと考えられる。

2.42.42.4

Kt=3.02.22.22.2

Kt=2.0Kt=4.0pressure ratio
Total pressure[-]
ratio [-]
Total
Total pressure ratio [-]22 2

1.81.81.8

Kt=5.01.6
1.61.6
1.4
1.41.41.2
1.21.2

Exp.Exp.
Exp.
DESDES
DES
RANS
RANS
RANS
DES_lissajous
DES_lissajous
DES_lissajous

11 1
-0.2
-0.2-0.2

0
-0.1
00 0
0.1
-0.1-0.1
0.10.1
Mass
flow
rate
[kg/s]
Mass
flow[kg/s]Massflowraterate[kg/s]

0.2
0.20.2

Fig. 2.2.3-14 Performance curves of the compressor.
Figure 2.2.3-15 に示す A～D の 4 つの時刻における入口配管の中心断面における軸方向速度分布を Fig. 2.2.3-16 に示す。時刻 A は圧縮機入口の流量が負の値になる直前であり、時刻 B は圧縮機入口の逆流流量がピークの瞬間、時刻 C は配管入口の逆流流量がピークの瞬間、時刻 D は圧縮機入口
および配管入口の流量が正の値に戻った直後である。Figure 2.2.3-16 に示すように、入口配管において逆流は壁面付近のみに存在することが明らかになった。この逆流はインペラの下流から発生しており、旋回速度成分を有するため、遠心力によってこのような偏りが発生していると考えられる。
Figure 2.2.3-17 に 4 つの時刻について θ=0°の断面におけるディフューザ内の子午面速度分布とスクロール内の周方向速度分布を比較している。同図より、配管入口および圧縮機入口の逆流流量がピークに達する瞬間には、スクロールの断面全体にわたって逆流が発生していることがわかった。

- 38 -

A BC

500

D

400

Mass Flow [%]

300
200
100
0
-100
200

qin

配管入口
qout
圧縮機入口
配管出口

-200

Mass flow rate [%]

-300

150

-400

q_hw

-500

100 0

500

1000
t* [-]

1500

2000

50

Fig. 2.2.3-15 Time histories of the flow rate.

0
0

Cz[-]

1000

2000

3000

t* [-]

0.0

-0.2

配管入口

A
圧縮機

B
C
①

② ③

④

D
①
④ inlet pipe.
③
Fig. 2.2.3-16 Axial velocity distributions
in②the central cross-sectionof the
配管入口(a) time A.

(b) time B.

(c) time C.

(d) time D.

Fig. 2.2.3-17 Meridional velocity distributions in the diffuser and circumferential velocity
distributions in the cross-section located at θ=0°.

- 39 -

(4) リアルワールド自動車空力性能の予測(a) 研究成果の要旨
本実証研究テーマでは、電気自動車や自動走行車等の、次世代自動車の設計・開発に貢献するための HPC シミュレーション技術の構築と実証を行うことを目的として、実走行状態の自動車の空
力性能（空気抵抗、操安性、横風安全性等）を車体形状再現性を数ミリまで高めた高解像度・長時間シミュレーションによって評価する。
令和 3 年度は、「富岳」において、リアルワールド自動車空力性能予測のための解析を実施し、自動車マルチボディダイナミクスと車両空力の双方向連成解析システムおよび大気変動風発生モジュールを開発し、「富岳」を駆使した、高解像度かつ長時間のシミュレーションを必要とする、大気変動風を考慮した実走行時の燃費および操安性評価のためのテスト計算を実施し、定性的な妥当性を確認した。その結果、一様流中における車体固定状態とドライバー操舵あり、さらには変動風中
での車体固定ケースとドライバー操舵ありの 4 ケースに関して、空気抵抗（燃費）が、操縦安定性に無視できない最大 9%の変化が生ずることがわかった。さらに、レーンチェンジ時において変動
風がドライバー操舵に与える影響が明らかになった。
The purpose of this research is to build and validate the framework applying HPC technology
to contribute to the design and development of next-generation vehicles such as electric vehicles
or autonomous vehicles. To reach this goal, the evaluation of the aerodynamic characteristics,
such as the drag force, the vehicle maneuverability, and the crosswind safety, of the vehicle inpractical driving conditions, must be done with a high-resolution computational mesh for
precisely representing the vehicle shape in millimeter-scale and long-term simulations arenecessary, all of which is carried out in this research.
In this fiscal year, continued from the last year, a coupling simulation framework of automotive
multi-body dynamics and vehicle aerodynamics, and a module for atmospheric turbulence
generator have been developed to realize the real-world evaluation of vehicle aerodynamics
performance. The developed simulation methods were applied to evaluate real-world fuelconsumption rate and maneuverability of a hatchback vehicle in on-road conditions. In these test
simulations, various evaluations were made by changing the incoming turbulence characteristics
and driver’s steering action and corresponding vehicle’s motion. Namely, we compared four caseswith and/or without the incoming turbulence and with and/or without the vehicle’s posturechange on-road condition by the driver’s steering action, and found that the aerodynamic dragchanges by as much as 9 percent among these cases. It was also investigated that how incoming
turbulence affects the driver’s steering action.

- 40 -

(b) 本研究の内容
1)「富岳」での基盤ソフトウェア CUBE の性能評価
本年度は、よりプロダクトランに近いセッティングにおいて、基盤ソフトウェア CUBE の性能測定を実施した。まず、CUBE の非圧縮性流体ソルバーを使用し、自動車モデルを使用したスケーリ
ングの測定を行った。測定時の CUBE の実行パラメータは Table 2.2.4-1 のとおりである。Table
2.2.4-2 と Fig. 2.2.4-1 にプログラム全体（時間積分ループ）のウィークスケーリングの測定結果を示す。51,200 ノードのスケーリングは 100 ノードの場合との比で 74.11 %を示し、階層格子を使用した実問題のセッティングにおける性能としては良好な値となった。
次に、CUBE の圧縮性流体ソルバーを使用し、自動車モデルを使用した性能評価を行った。測定
時の CUBE の実行パラメータは Table 2.2.4-3 のとおりである。時間積分ループ内の各処理のコスト分布を Fig. 2.2.4-2 に示す。以前性能ベンチマーク用に測定した Cavity flow の問題の時と異なり、粘性項や対流項といった主要計算部のコストよりも、LUSGS（Lower-Upper Symmetric
Gauss Seidel）法とよばれる時間発展スキームのコストが大きなもの（44.2 %）となることが確認できた。さらに、この LUSGS 法のサブルーチンを対象に、性能チューニングを実施した。LUSGS
法のサブルーチン内のコストを分析し、最も大きいコスト（サブルーチン全体の 26.7 %）を持つ
Sweeping のループを切り出し、この部分をチューニングのターゲットとした。Table 2.2.4-4 と Fig.
2.2.4-3 にカーネルループの性能値を示す。このループは非常に長いループボディ（約 250 行）の中にループの反復に関する依存性を持つ処理を含んでおり、コンパイラによる SIMD 化やソフトウェ
アパイプライニング化が実行されなかったが、配列の次元の変更、ループの一重化、ループ分割、条件分岐の削除といった最適化により、最終的には 1.45 倍の性能向上が得られた。
Table 2.2.4-1 Simulation parameter setups.
Solver type
Domain size
Number of cubes
Number of cells per
cube
Total number of cells
Minimum cell size
Time step size
Integration time
Reynolds numberTime integration
schemePoisson solver

Incompressible
18.75m×8.75m×3.75m
32,341
16×16 ×16 = 4,096
32,341×4,096 = 132,468,736
1.2207 mm
2×10-6 sec
2×10-4 sec (100 steps)
2.572×106
Crank-Nicolson
Red-Black SOR

Table 2.2.4-2 Scaling of time-stepping loop.
Number of
nodes
100
800
6,400
51,200

Number of
cubes
32,341
260,3312,082,648
16,661,184

Number of cells
132,468,736
1,066,315,776
8,530,526,208
68,244,209,664- 41 -

Elapsedtime
22.92 s32.63 s33.46 s30.92 s

Scaling100.00 %70.23 %
68.49 %
74.11 %

Fig. 2.2.4-1 Elapsed time against number of compute nodes.
100%

Subroutine name90%tstep_lusgs_slts_mix

44.2 %

flux_roe_lm

17.6 %

field_intf_unpack_e

12.6 %

cns_absorbing_bc

5.7 %field_intf_unpack_edge_e

5.1 %vscf_ce2nd

4.8 %field_intf_q_e

2.0 %

field_intf_unpack_corner_e

1.7 %field_intf_pack_e

1.4 %vscf_ce2nd
field_intf_unpack_edge_e
80%
cns_absorbing_bc
70%

field_intf_unpack_e

60%
50%

flux_roe_lm

Cost ratio

40%
30%
20%

tstep_lusgs_slts_mix

10%0%Fig. 2.2.4-2 Cost distribution in time-stepping loop.
Table 2.2.4-3 Simulation parameter setups.
Solver type
Domain size
Number of cubes
Number of cells per
cube
Total number of cells
Minimum cell size
Time step size
Integration time
Reynolds numberTime integration
schemeViscous term

Compressible
20m × 10m × 5m
32,341
16 × 16 × 16 = 4,096
32,341 × 4,096 = 132,468,736
1.2207 mm
1.0 × 10-4 sec
1.3 × 10-3 sec (13 steps)
1.1 × 106
Lower-Upper Symmetric-GaussSeidel (LUSGS) implicit method
2nd order central difference- 42 -

Table 2.2.4-4 Execution performance of kernel loop.
Original
Exchange array
dimension
Loop collapse &manual fission
Auto loop fission
Extract if
condition

Exec. time
5.08 sCalc. performance
19.7 GFLOPS (2.57%)

Mem. bandwidth
5.92 GB/s (2.31%)

4.50 s22.23 GFLOPS (2.90%)

6.24 GB/s (2.44%)

4.05 s

11.63 GFLOPS (1.52%)

6.56 GB/s (2.56%)

3.84 s12.29 GFLOPS (1.60%)

7.19 GB/s (2.81%)

3.49 s

14.97 GFLOPS (1.95%)

8.01 GB/s (3.13%)

Fig. 2.2.4-3 Tuning results of kernel loop.
2) 大気変動風発生モジュールの開発
昨年度に引き続き，変動風作成モジュールの開発を行った。変動風は Fig. 2.2.4-4 に示すとおり、複数の翼型を車両上流にスパン方向に並べ、ピッチングおよびヒービング運動を組み合わせて与えて作成する。翼のサイズや配置については、既存の変動風生成装置を参考にしている。それぞれの運動の周波数と振幅、さらに位相差を変化させて、生成する乱れの積分長さ、スペクトル、乱流強度のコントロールを行う。昨年度は特に作成した変動の積分長さについて、実走行時のデータと比較して十分な結果を得ることができなかった。本年度は特に翼まわりの格子解像度の観点からこの問題を精査し、2 mm 程度の解像度を確保することで翼表面の不自然な剥離が解消され、これによ- 43 -

り十分な積分長さを有する変動風の作成に成功した。Figure 2.2.4-5 に作成した変動風の一例として瞬時速度分布として示す（図には、ピッチ角 6 度で周波数が 0.1 Hz、0.5 Hz、1 Hz の 3 ケースについて示している。
）。ここでは翼に与えるピッチ角の正弦波振動を 2 度から 6 度まで変化させ、振動周波数を 0.1 Hz から 1 Hz まで変化させて変動風を作成している。昨年度見られた翼からの非物理的剥離に伴う変動の減衰はここでは見られず、与えた翼ピッチ運動に応じた変動が再現されているのが確認される。このようにして得られた変動風の性状を翼から 3 メートル下流（車体の前縁を設置する場所）で計測した変動風の振れ幅（ヨー角）、乱流強度、積分長さ、主流方向積分長さと積分長さとの比として Fig. 2.2.4-6 に示す。ピッチ角周波数を 0.1 Hz 程度に設定することで、実走行試験で得られているのと同程度の数メートルの積分長さの変動が作成されている。なお、Fig.
2.2.4-5 で確認できるとおり、翼からの剥離が抑制された影響で乱流の性状にスパン方向の「むら」が発生している。この問題については、次年度の課題として引き続き検討を進める必要がある。

Fig. 2.2.4-4 Setups of the turbulence generator.

- 44 -

Fig. 2.2.4-5 Flow fields (velocity magnitude) reproduced by turbulence generator (top: =6 deg.,
f=0.1 Hz; middle: =6 deg., f=0.5 Hz; bottom: =6 deg., f=1 Hz).

- 45 -

Fig. 2.2.4-6 Reproduced “yaw variance”, “turbulence intensity”, “the overall integral lengthscale”, and “ratio of streamwise integral length scale and overall integral length scale”,
as a function of imposed sinusoidal pitching oscillation of the airfoils (amplitude andfrequency).
3) 自動車運動連成空力解析システムの概要
まず、昨年度開発した車両運動連成空力解析システムの概要を示す。車両運動解析には、
University of Parma、University of Wisconsin-Madison を中心にオープンソースで開発が進められている Project Chrono（Chrono）を用いた。Chrono は車両運動解析用のライブラリ群の開発が積極的に進められており、商用のソフトウェアとほぼ同等の精度を有することが確認されている。
Chrono におけるマルチボディダイナミクス解析の基礎方程式には、オイラーの運動方程式を用いた。ここで𝑴は剛体の質量、𝒓̈ は加速度、𝑱は慣性モーメント、𝝎̈は角加速度、𝑭は外力、𝝉はトルク、
̅ 𝑇 (𝜱)は高速条件による反力、および、𝝀はラグランジュの未定乗数をそれぞれ示している。
𝜱𝑇𝑟 , 𝜫𝑴𝒓̈ + 𝜱𝑻𝒓 𝝀 = 𝑭
̅ 𝑻 (𝜱)𝝀 = 𝝉𝑱𝝎̈ + 𝜫タイヤと地面との接地面に対しては、前後力（Longitudinal force）、横力（Lateral force）、セル
フアライニングトルク（Self Aligning Torque: SAT）が発生しており、これらの力により、車両は進行方向を変化させている。しかし、これらの力はタイヤの物性値や空気圧などさまざまな影響を受けるため、非線形特性を有しており、車両運動解析などでは，タイヤ力をモデル化したものを用いた解析が行われている。本研究ではタイヤモデルとして Pacejka89 を用いた。Pacejka89 モデル
はタイヤ試験機を用いて実測された値𝑎𝑖 (𝑖 = 0~13)、𝑏𝑗 (𝑗 = 0~10)、𝑐𝑘 (𝑘 = 0~17)を用いてタイヤ力のモデル化を行っており、スリップ角𝛼とスリップ比𝜅を元に作用するタイヤ力を決定する。
自動車はエンジンで生み出された回転エネルギーをトランスミッションによって適切なトルクと回転速度に変速し、タイヤに伝達することで駆動力を得ている。本研究では、モデル化したパワートレインを用いることで、ドライバーからの入力の影響を含めたパワートレインの再現を行っている。エンジンはクランクシャフトを通じて，トランスミッションとつながっており，トランスミッ
ションはプロペラシャフトを通じてディファレンシャルギアにつながっている。ここにエンジンの- 46 -

トルクマップを考慮し、プロペラシャフトの回転数とギアの変速比からエンジンの回転数が求まり、トルクカーブより、回転数に応じた最小・最大トルクを得ることができる。
ドライバーモデルはステアリングを操作し、車両の進行方向を制御する操舵モデルと、一定速度を維持するためにアクセル、ブレーキ操作をする加減速モデルから構成されている。操舵モデルには入力として走行コースを 3 次ベジエ曲線によって与える。時刻 t におけるドライバー視点の位置を𝑫𝒑𝒕 、ドライバーの視線が向けられている点をSentinel point （𝑺𝒑𝒕 ）、𝑺𝒑𝒕 から最も近いベジエ曲
線上の点をTarget point（𝑻𝒑𝒕）とする。時刻 t における視線ベクトル𝑺𝒕と，ターゲットベクトル𝑻𝒕は次式で表すことができる。
𝑺𝒕 = 𝑺𝒑𝒕 − 𝑫𝒑𝒕
𝑻𝒕 = 𝑻𝒑𝒕 − 𝑫𝒑𝒕
よって，ドライバーの視線と，目標コースとの差errは次式で表すことができる。
𝒆𝒓𝒓 =(𝑺𝒕 × 𝑻𝒕 ) ∙ [𝟎, 𝟎, 𝟏]|𝑻𝒑𝒕 − 𝑺𝒑𝒕 ||(𝑺𝒕 × 𝑻𝒕 ) ∙ [𝟎, 𝟎, 𝟏]|

errを偏差として PID 制御を行うと、ドライバーのステアリング入力𝜏𝑠𝑟𝑡 （ − 1 ≤ 𝜏𝑠𝑟𝑡 ≤ 1）は次式
で表すことができる。
𝒕̇𝝉𝒔𝒓𝒕 = 𝑲𝒑 𝐞𝐫𝐫 + 𝑲𝒊 ∫ 𝒆𝒓𝒓 ∙ 𝒅𝒕 + 𝑲𝒅 𝒆𝒓𝒓𝒐

マルチボディダイナミクス解析を行うステアリング機構に𝜏𝑠𝑟𝑡 強制変位として与えることにより、ドライバーのステアリング操作を車両運動に反映させる。ステアリング操作のモデルの概要を Fig.
2.2.4-7 に示す。

Fig. 2.2.4-7 Model of Steering control.
自動車運動連成空力解析フレームワークの有用性を確認するために用いた解析モデルを Fig. 2.2.48 に示す。対象は、実在のステーションワゴン車とし、本田技術研究所より提供された実開発で用いられる形状再現性を有する CAD を用いて解析を行った。

- 47 -

Fig. 2.2.4-8 Target model.
4) リアルワールド燃費評価
実走行時の燃費評価を目的として、変動風作成モジュールと自動車運動連成空力解析フレームワークを用い、変動風下で 100 km/h で直進走行するシミュレーションを行った。Figure 2.2.4-9 に走行中の車体周りの瞬時流速を示す。Figure 2.2.4-10 に一様流下・変動風下での、CUBE のみでの解析（すなわち車両運動を考慮しない）と Chrono との連成解析（直進走行するためにドライバーの修正操舵とそれに伴う微小な車両運動を考慮）の、それぞれの場合における標準化された空気抵抗係数の値（𝐶𝐷 値）を示す。なお、翼を置いた場合には、そのウェイクによって車体周囲に圧力勾配が発生する為、車体を置かずに翼だけでの解析を別に行い、圧力勾配によって発生する力や車体に実際に当たる流速を求め、𝐶𝐷 値を算出した。結果を比較すると、変動風下では𝐶𝐷 値が 6 %程度増加している。特に、Chrono との連成解析によって車両運動を再現した場合では、一様流下に対して変動風下では𝐶𝐷 値が約 9 %増加した。この要因として、車体のヨー角変動とドライバーの操舵入力の変化が挙げられる。変動風下での方がヨー角変動、操舵入力が共に大きくなっていることが分かる。このことから、変動風によって車体運動が大きくなることに伴い、𝐶𝐷 値が増加するという現象を捉えられたと考えられる。
今後、変動風モジュールのパラメータを変更し、積分長さや乱流強度を変化させた場合に、空気抵抗にどのような影響を与えるのかを評価する必要がある。

Fig. 2.2.4-9 Snapshots of the velocity profile w/o (left) and w/ (right) incoming turbulence.

- 48 -

Fig.2.2.4-10 𝐶𝐷 in each case.

Fig.2.2.4-11 Yaw angle time history.
5) リアルワールド操安評価
ここでは 4)と同じフレームワークを用い、100 km/h で走行する車両が解析開始 1 秒後から 3 秒間
でレーンチェンジ運動をする際の解析を行った。Figure 2.2.4-12 にドライバーの操舵入力・ヨー角の時間変化と変化が比較的大きい部分を拡大した図を示す。レーンチェンジ運動の間で操舵入力がほぼ一定となる時間（1.8 s~2.0 s）で変動風を与えたことによる変化が大きく現れている。その要因として、Fig. 2.2.4-13 で示す、1.8 s から 2.0 s でのヨーイングモーメントとヨー角速度の変化が考えられる。変動風によって生じたヨーイングモーメントの変化にドライバーが反応し、修正操作を行った結果、ステアリング入力量に差が生じたと考えられる。
4)の結果からも、変動風が燃費・操縦安定性評価に影響を与えることが見てとれた。特に変動風の精度向上が可能になると、このフレームワークを活用して、実走行環境において燃費や操安性に寄与する空力パーツの開発、その物理メカニズムの解明が可能になると期待される。

- 49 -

Fig.2.2.4-12 Steering input (top) and Yaw angle (bottom).

Fig.2.2.4-13 Yawing moment (left) and Yaw rate velocity (right) between 1.8 s and 2.0 s.
(5) リアルワールド自動車空力音予測(a) 研究成果の要旨
本実証研究テーマでは、電気自動車や自動走行車等の、次世代自動車の設計・開発に貢献するための HPC シミュレーション技術の構築と実証を行うことを目的として、ボンネット隙間やフロン
トグリルから発生する狭帯域音、ピラーやドアミラーから発生する広帯域音、空力・構造振動・音響連成による車内騒音を対象に、実走行時の空力音予測を行う。

- 50 -

令和 3 年度は、「富岳」において、リアルワールド自動車空力音予測のための解析を実施た。自動車の狭帯域音を予測するために、圧縮性流体ソルバーCUBE を使用し、フロントグリルから発生するフィードバック音の直接解析を実施した。その結果、音源構造に関する新たな知見として、剥離したせん断層がエンジンフードと干渉して音が発生している様子が確認された。また、自動車の広帯域音予測のため、熱流体解析システム FFB と空力音響解析システム FFB-ACOUSTICS を組み合わせた分離解法と比較して低計算コストを実現する LBM 法（格子ボルツマン法）に基づいた流体
解析システム FFX を用いて車体まわりの空力音響計算を実施した。FFX によって得られる空気抵
抗係数のレイノルズ数依存性は、分離解法を用いた結果と比較して、同程度の精度で得られ、FFXを用いた解析の妥当性を検証した。そして、FFX を用いた直接解析によって、フロントグリル、ドアミラー、およびホイールハウスに強い音源があることがわかった。
The purpose of this research theme is to establish HPC simulation technology and demonstratethat it can contribute to the design and development of next-generation automobiles including
electric vehicles. For this purpose, the narrow-band noise generated from the hood gap and thefront grille, the broadband noise generated from the pillars and the door mirrors are simulated.
Moreover, the interior noise are predicted by coupling the external aerodynamics-aeroacousticalsimulation, the structural-vibration simulation, and the interior acoustics simulation.
In this fiscal year, continued from the last year, the narrow-band feedback noise emitted between
the vehicle’s body and the engine hood and the broad-band noise generated from the wholevehicle’s body were simulated, as a part of the real-world vehicle noise simulations. For the
narrow-band noise, by increasing the grid resolution around the gap of the engine hood, thepredicted noise spectra showed a better agreement with the measurements than the last year’s
results. For the broad-band noise prediction, direct simulation of the noise was performed by FFX
based on the Lattice-Boltzmann Method (LBM) and the results were compared with thoseobtained by the conventional splitting method of the flow field simulated by FFB and theacoustics field simulated by FFB-ACOUSTICS. It was confirmed that the results of the direct
lattice-Boltzmann method show a good agreement with the conventional splitting method basedon the Navier-Stokes equations for the prediction of the Reynolds-number dependence of the
aerodynamic drag. The direct simulation by the FFX also revealed that the front grill, the door
mirrors, and the wheel-houses compose the dominant source of the vehicle’s external noise.
(b) 本研究の内容
1) 狭帯域音予測
自動車の狭帯域音を予測するため、昨年度に引き続き、Fig. 2.2.5-1 に示す車両まわりの音の直接解析を行った。流れの中に微小な段差・隙間が置かれると、発生した空力音が渦自体に影響を及ぼすフィードバック音が生じることが知られているが、一般にこのような空力的な自励音は上流の流れが層流である場合に起きることや、発生した渦が後方の障害物のエッジに衝突するような場合に発生することがわかっている。しかし、自動車のように形状が複雑な場合、上流の流れ（境界層）- 51 -

の状態を確認することが難しく、また、本モデルのようにボンネット内のエンジン等の圧力損失に
より主流流れが変化し、はく離した流れの衝突箇所が予測できない場合がある。自動車から放出される狭帯域音を予測する場合、音が発生する部位だけでなく、車両全体の流れ場を考える必要があ
るのはこのためである。本研究では、この問題を明らかにするため、狭帯域音が発生する箇所だけ
でなく、エンジンルームを含む車両全体の解析を行った。エンジンルームを含む複雑形状を解析するために直交格子による解析アプリである CUBE を使用した。また、空力音を解析するために圧縮性ナビエ・ストークス方程式を解き、音を直接求めることを試みた。音の発生個所であるボンネッ
ト付近は小さなキューブを用いて計算格子の解像度を上げ、車両まわりの解像度は車体まわりの比較的大規模な流れが再現できるように調整した。

Fig. 2.2.5-1 Computational grid for narrow-band sound analysis.
計算条件としては、実車においてボンネットから狭帯域音が発生した 60 km/h に車両速度を設定した。昨年度、最小格子サイズ 0.2 mm（格子数 20 億）としてテスト計算を行った結果、狭帯域音の発生の再現には成功したものの、音圧レベルは実験とは一致しなかった。本年度はこの原因を精査した結果、車体表面から剥離したせん断層とボンネットへの衝突、それに伴う音波の発生メカニ
ズムが実験とは定性的に異なっていたことがわかった。そこでこの問題を解決するために、ボンネット近傍の一部の解像度をさらに 0.05 mm に高めた解析を行い、空間解像度がせん断層挙動に与える影響について検討を行った（Fig. 2.2.5-2 参照）。
Figure 2.2.5-3 にその結果を示す。空間解像度が 0.1 mm の場合（下）は、車体表面から剥離したせん断層は、下流のエンジンフードと強く干渉することなくエンジンルーム内に侵入し、エンジンルーム内の構造物に衝突することで空力音が発生しているが、0.05 mm の場合（上）は、剥離したせん断層がエンジンフードと干渉して音が発生している様子が確認される。一方、エンジンフードと車体のギャップは車体横方向全般にわたっており、この領域全体に 0.1 mm の格子解像度を確保した場合、メッシュ数は数兆に及び、「富岳」全系を用いてもこのような解析は難しい。そこで現- 52 -

在、0.2 mm 程度の格子解像度でもこの現象を再現できるよう、より数値粘性の少ない（すなわち剥離せん断層の厚さが過大評価されない）数値スキームの実装を行っている。

Fig. 2.2.5-2 Mesh generation to increase the spatial resolution around the engine hood to capture
the narrow-band noise and reproduced flow field.

Fig. 2.2.5-3 Reproduced narrow-band noise and its dependence on the mesh resolution (top:0.05mm, bottom: 0.1mm). Velocity (middle) and pressure (right) fields.

- 53 -

2) 広帯域音予測2-1) FFX を用いた空力音響解析
空力音を解析する手法として、流体解析と音響解析を別々に行う分離解法、および、圧縮性ナビ
エ・ストークス方程式を解くことによって空力音を求める直接解法が利用されてきたが、分離解法
では、音源のモデル化に課題があり、直接解法は計算コストの課題があった。格子ボルツマン法は、
圧縮性ナビエ・ストークス方程式を直接解くよりも計算コストが小さく、また並列化しやすいことから HPC に適した解析手法であることが知られている。格子ボルツマン法は空力音の直接計算手法であることからも自動車空力音の解析に適した手法である。また、直交階層格子を用いることによって、CAD データから自動的に解析格子を生成することが容易なこともあり、産業利用にも適した解析手法である。
本年度は、HPC に適した大規模な格子ボルツマン法に基づいた汎用流体解析システム FFX を開発
し、自動車周りの流れ場の解析を試みた。
格子ボルツマン法汎用流体解析システム FFX の開発
2020 年度までに開発した格子ボルツマン法に基づく汎用流体解析システム FFX の計算精度、ロバ
スト性、利便性を向上させるため、衝突項、粒子モデル、サブグリットスケールモデル、流体力計
算機能、インターセクトデータ作成機能に関して、拡張・改良し、基礎検証を実施した。ここでは、これらの機能について報告するとともに、開発したコードを用いた検証結果を報告する。なお、
2020 年度までの開発において、ノードあたり 4300 万格子を用いた Weak Scale ベンチマークテストを実施し、最大で「富岳」の 32,768 ノードを用いて、4.4 PFLOPS を達成し、1.4 兆グリッドを用いた計算を 1 時間ステップあたり 0.08 秒で計算できることを確認している。
a) 衝突モデル
LBM は直交格子上に定義された粒子の分布関数の時間発展を計算する。分布関数𝑓𝑖 (𝑡, 𝑥𝛼 )の時間発
展は粒子の移動（以降、並進と称する）および衝突で表される。衝突は格子における分布関数の変化を介して粒子間の密度，運動量およびエネルギーの交換を表す。分布関数の時間発展を記述する
方程式を格子ボルツマン方程式という。以下に格子ボルツマン方程式を示す。
𝑓𝑖 (𝑡 + Δ𝑡, 𝑥𝛼 + 𝑐𝑖,𝛼 Δ𝑡) = 𝑓𝑖 (𝑡, 𝑥𝛼 ) + Ω𝑖 [𝑓𝑖 (𝑡, 𝑥𝛼 )](2.2.5-1)

ここでΩ𝑖 は衝突項であり、FFX のプロトタイプでは、以下に示す BGK モデル[1]が実装されていた。
BGK モデルは、各粒子の衝突項が、その粒子の分布関数の非平衡成分のみに依存しており、Single
Relaxation Time（SRT）とも称される。本年度は、衝突モデルとして、SRT に加え、Multi
Relaxation Time（MRT）[2]を実装した。式の記述は割愛するが、MRT では、速度空間における分布関数をモーメント空間に変換し、モーメント空間における非平衡成分を速度空間に変換し直す
ことにより衝突項を計算する手法であり、演算数が増えるものの計算の安定性が向上することが報告されている。
衝突モデルと計算精度、安定性の関係を調査するため、初期速度の変動値と領域の長さを基準とするレイノルズ数 121 の一様等方性乱流の計算を実施した。Figure 2.2.5-4 にエネルギースペクト
ルの比較、Fig. 2.2.5-5 に乱流構造の可視化結果を示す。FFX の結果は、SRT、MRT ともに、DNSによる結果[3]と一致しており、MRT は SRT と同等の計算精度を有することが確認できた。

- 54 -

Fig. 2.2.5-4 Energy spectra of homogeneous isotropic flow computed by DNS [3] and FFX with
SRT and MRT.

(a) SRT.

(b) MRT.

Fig. 2.2.5-5 Iso-surface of second invariant of homogeneous isotropic flow（left：SRT, right:
MRT).
一般に、LBM による計算はレイノルズ数が高くなると、緩和時間係数が 0.5 に近い値となり、計算が不安定になることが多い。計算の安定性は、壁面の有無、速度勾配等、さまざまな条件にも影響を受ける。壁面や平均的な速度勾配の影響のない条件下での安定性を把握するため、FFX を用いてレイノルズを上げた条件に対して一様等方性乱流のベンチマーク計算を実行したところ、SRT、MRT の安定限界（安定に計算できる緩和時間係数の最低値）はそれぞれ、0.507、0.504 であり、MRT のほうが若干安定であった。

- 55 -

b) 粒子モデル
LBM は直交格子上に定義された粒子の分布関数の時間発展を計算する。密度、速度等のマクロ変数は分布関数を積分することにより計算できる。各粒子の速度は、速度と時間刻みの積が格子幅になるように定義される。したがって，ある時刻ある格子に存在する粒子は、次の時間ステップでは隣の格子に移動する。このため、対流や伝播現象を正確に予測することができる。粒子モデルとして、3 次元の場合、自分の点を含め 15 点、19 点あるいは 27 点の隣接点に対応した速度を定義する
選択肢があり、これらをそれぞれ 3D15Q、3D19Q、3D27Q と称する。FFX のプロトタイプには、
D3Q15 が実装されていたが、本年度はこれに加え D3Q27 による計算機能を追加した。粒子モデルと計算精度、安定性の関係を調査するため、a)に示した一様等方性乱流の計算を実施した。粒子モデルには SRT を用いた。Fig. 2.2.5-6 にエネルギースペクトルの比較、Fig. 2.2.5-7 に乱流構造の可視化結果を示す。どちらの粒子モデルのエネルギースペクトルも DNS による結果に近く、同等の計算精度があることがわかったが、D3Q27 のほうが高波数領域において若干乱流エネルギーが小さく、これ対応して、乱流構造の可視化結果において、SRT の結果にある細かな構造が消えている。
安定限界を調べるため、レイノルズ数を徐々にあげた計算を実行したところ、D3Q15 は、緩和時間
係数 0.5059 の計算では計算が発散したのに対し、D3Q27 では緩和時間係数を 0.5015 まで下げても計算は発散しなかった。このことより、D3Q27 は、D3Q15 と比較して計算の安定性が向上することが確認された。

Fig. 2.2.5-6 Energy spectra of homogeneous isotropic flow computed by DNS[3] and FFX with
D3Q15 and D3Q27.

- 56 -

(a) SRT, D3Q15.

(b) SRT, D3Q27.

Fig. 2.2.5-7 Iso-surface of second invariant of homogeneous isotropic flow（left：D3Q15, right:
D3Q27).
c) サブグリットスケールモデル
FFX には SubGridScale（SGS）モデルとして標準スマゴリンスキーモデルが実装されている。ス
マゴリンスキーモデルでは、渦粘性を計算する際にひずみ速度が必要となる。FFX のプロトタイプ
では、分布関数から速度を計算し、これの差分を取ることによりひずみ速度を求めていたが、分布関数の非平衡成分からひずみ速度を計算する方法[4]に変更した。以下にひずみ速度テンソルと分布関数の非平衡成分の関係を示す。
𝑆𝛼𝛽 = −

3𝑒𝑞
∑ 𝑐𝑖𝛼 𝑐𝑖𝛽 (𝑓𝑖 − 𝑓𝑖 )2𝜏

(2.2.5-2)𝑖

速度場の差分による方法では隣接格子の情報が必要となるため、並列計算においては適切な通信処理をしなければひずみ速度の連続性を保証できない。通信処理コストの増大を避けるため、この通信処理を省略していたため、ひずみ速度およびこれから計算される渦粘性に CUBE 境界において不連続が生じ、これに起因した振動が生じることが確認されていた。分布関数の非平衡成分を用いる手法では自格子の情報のみを参照するので不連続の問題は生じず、振動の問題が解消された。
ひずみ速度の計算手法の妥当性を確認するため、a)および b)に示した一様等方性乱流の分布関数を用いて、二通りの手法でひずみ速度を計算した。Fig. 2.2.5-8 は、横軸に速度場の差分、縦軸に分布関数の非平衡成分から計算したひずみ速度をプロットした。分布関数による計算は、分布関数のChapman-Enskog 展開に基づく近似式であるため、速度場の差分により計算したひずみ速度と完全には一致しないが、高い相関があることが確認できた。

- 57 -

Fig. 2.2.5-8 Strain of homogeneous isotropic flow computed by two different method (horizontal
axis: by velocity field, longitude axis: by distribution function).
d) 流体力計算機能
ナビエ・ストークス方程式に基づく流体解析手法において壁面を取り扱う場合、壁上に定義される速度にゼロをセットする。LBM における壁面は、隣接する計算格子を遮断することによって表現される。壁がない隣り合う格子は並進処理により分布関数の情報が引き渡されるが、壁がある場
合は、その代わりに、壁面で粒子が跳ね返る（Bounce Back）。格子間にある壁面の位置により、粒子の跳ね返りを計算する手法である Interrelated Bounce Back（以下、IBB）[5]が、FFX にも実装されている。
α方向にインターセクトがある格子 A に対し、IBB に基づきインターセクト処理に伴い、格子 A に生じる運動量変化より、インターセクト面に作用する流体力𝐹𝑖 を求める。α方向の分布関数𝑓𝛼 は、インターセクトが無ければ並進処理により隣の格子にコピーされるが、α方向にインターセクトがある場合はどこからも参照されずに、隣の格子の分布関数で置き換えられる。これは、計算領域全体の運動量バランスの観点からみると、格子 A が運動量Δ3 𝜌0 𝑐𝑐𝑖,𝛼 𝑓𝛼 を失ったと解釈することができる。
ここで、𝑐𝑖,𝛼 は、無次元粒子速度である。α方向にインターセクトがある場合、α方向の逆方向（これをβ方向とする）の分布関数𝑓𝛽 は IBB に基づき計算され、𝑓𝛽𝐼𝐵𝐵 に書き換えられる。これによる格子 A における運動量変化は、Δ3 𝜌0 𝑐𝑐𝑖,𝛽 (𝑓𝛽𝐼𝐵𝐵 − 𝑓𝛽 )となる。
LBM の時間刻みΔ/𝑐の間に、格子 A に作用する流体力 −𝐹𝑖 （インターセクト面に作用する流体力×−1）が作用することにより、上記の運動量変化が生じたと考えれば以下の関係が成立する。
Δ−𝐹𝑖 × ( ) = −Δ3 𝜌0 𝑐𝑐𝑖,𝛼 𝑓𝛼 − Δ3 𝜌0 𝑐𝑐𝑖,𝛽 𝑓𝛽𝑐(2.2.5-3)

整理すると、流体力として以下を得る。
𝐹𝑖 = 𝜌0 𝑐 2 Δ2 [𝑐𝑖,𝛼 𝑓𝛼 − 𝑐𝑖,𝛽 𝑓𝛽 ](2.2.5-4)

これを全壁面について足し合わせれば物体に作用する流体力を計算することができる。

- 58 -

FFX に実装された流体力計算機能の妥当性は角柱まわり流れの計算、車体まわり流れの計算において検証しており、これについては後述する。
e) インターセクトデータ作成機能
前述したとおり、LBM における壁面は計算格子間の遮断によって表現される。本資料では遮断に係る情報（計算格子のインデックスおよび格子と壁面の距離）をインターセクトデータと称する。
2020 年までに開発した FFX のプロトタイプでは、インターセクトデータを前処理プログラムで作成し、流れ計算の実行時にこれを入力していたが、本年度、解析システムの利便性向上のため、インターセクトデータ作成機能を流れソルバーに実装した。ナビエ・ストークス方程式に基づく流体解析で、複雑形状を含む解析対象を扱う場合、形状ファイルの修正や計算格子データの作成に多大な時間と手間を必要とする場合があるが、FFX による解析では、STL データ（三角形データ）さえ用意すれば、計算実行時に三角形データを入力するだけで、容易に計算を実行できるようになった。
三角形データからインターセクトデータを作成する際に、FFX プロトタイプでは三角形の面を点群でモデル化し、点群と隣接格子を結ぶ線（以下、格子線）の距離が一定値以下であればインター
セクトしているとみなす方法をとっていた。点群を用いるメリットは、点群が位置する計算格子を特定できるため検索処理が不要であり、点群あたりの演算回数が少ないことにあるが、格子解像度とともに、三角形から作成する点群の解像度をあげる必要があるため、大規模解析では点群の数が膨大になるデメリットがあった。これを解消するため、インターセクトデータを作成する方法を三角形と格子線の接触判定をするアルゴリズムに変更した。三角形と格子線の接触判定は、三角形を含む平面と格子線を含む直線の交点を求めるための 3×3 の連立方程式の求解に帰着する。この計算
では三角形あたり複数の計算格子に対し接触判定を必要とする検索問題になるが、隣接格子のパターンが決まっているため（D3Q15 であれば 7 パターン、D3Q27 であれば 13 パターン）、3×3 の連立方程式の逆行列をあらかじめ計算しておくことにより処理を高速化することができる。
開発したインターセクトデータを用いて、複雑形状を含む車体まわり流れの計算を実行した結果、複雑形状を含む解析対象に対し 1 兆格子規模の計算において、30 秒程度でメッシュデータを自動作成できることを確認している。
角柱まわりの流れ解析（解析精度の検証）
まず、解析精度を検証するため、角柱まわりの流れ場を解析した。
一辺の長さが 2 cm の矩形角柱を一様流れ（30 m/s）中に設置し、非圧縮性ナビエ・ストークス方
程式を用いた解析（FFB）と格子ボルツマン法（FFX）による解析を行った。代表寸法を角柱の辺の長さとしたレイノルズ数は 4.0×104 である。角柱表面の格子サイズは FFB では 0.04 mm、格子ボ
ルツマン法では 0.08 mm とした。解析空間は 0.8 m×0.8 m×0.04 m であり、角柱のスパン方向の境
界条件は周期境界条件とした。FFB の解析要素数は約 5,500 万、FFX の格子点数は 536 億である。
FFX における格子レイノルズ数は 3.1×103、格子ボルツマン法における計算安定性に影響を及ぼす
緩和時間係数は 0.5001、乱流モデルに用いたスマゴリンスキー定数は 0.1 である。また音速は 340
m/s とした。
Figure 2.2.5-9 から Fig. 2.2.5-12 に FFX により解析した角柱まわりの流れ場を示す。Figure 2.2.59 は主流速度場、Fig. 2.2.5-10 は渦渡場、Fig. 2.2.5-11 は速度場の発散を示す。Figure 2.2.5-9、Fig.

- 59 -

2.2.5-10 より角柱下流にカルマン渦が形成されていることがわかる。また、Fig. 2.2.5-11 の速度場の発散は音響場に相当することから、角柱の上下に位相が反転した音場が形成されていることがわかる。角柱下流の渦領域では、渦度変動による強い変動が観察される。自動車の車内騒音予測を行う場合、ドアラミー近傍にはこのような音場と圧力場が混在した領域が見られるが、車内音を解析する場合は、両者を分離する必要がある。Figure 2.2.5-12 に示す密度場でも同様に流体運動に伴う強い密度（変動）と音場として遠方まで伝達する密度（変動）が観察される。

Fig. 2.2.5-9 Flow pattern around a square cylinder.

Fig. 2.2.5-10 Vorticity field around a square cylinder.

- 60 -

Fig. 2.2.5-11 Distribution of divergence field around a square cylinder.

Fig. 2.2.5-12 Density field around a square cylinder.
Figure 2.2.5-13 に FFB との比較結果を示す。渦渡場の全体的な傾向が良く一致していることがわかる。

- 61 -

FFX (LBM)

FFB (NS-Eq)

Fig. 2.2.5-13 Comparison of vorticity fields simulated with FFX and FFB.
Figure 2.2.5-14 に流体力の比較結果示す。抗力、揚力ともに非圧縮性ナビエ・ストークス方程式
による計算とよく一致していることがわかる。Table 2.2.5-1 に実験結果の比較示す。Figure 2.2.59 から Fig. 2.2.5-12 の結果は FFX Ver.52.02 による結果である。開発過程で計算手法にいくつかの改良を加えているが、現状の FFX Ver.52.02 では平均抗力 Cd、変動揚力 Cl’ともに実験と非常に良く一致することが確認された。音響場の定量的な比較については現在、解析データ長が不足していることから求めることができなかったが、解析を続けることにより定量的なデータを得ることができると考えている。

Fig. 2.2.5-14 Comparison of time histories of aerodynamic force on a square cylinder.

- 62 -

Table 2.2.5-1 Comparison of aerodynamic forces on square cylinder.

St
Cd
Cl’

Exp.

FFB(NS)

0.134
2.22
1.45

0.125
2.06
1.37

FFX(LBM)
Ver.48.02
0.119
2.15
1.26

FFX(LBM)
Ver.50.03
0.123
2.10
1.56

FFX(LBM)
Ver.52.02
0.123
2.14
1.43

「富岳」を用いて本解析（約 500 億格子）に必要な解析時間は、ノード数 6,400 の場合、37 時間
20 分であった。ノード当たりの解析速度は 3 GFLOPS、ノード全体で 20 TFLOPS であった。解析
速度が低い理由して本解析ではノード当たり 1 領域で計算を実施しているためである。本来は 1
CMG あたり 1 領域となるようにノード当たり 4 領域に領域分割することが望ましいが、今回はメモリの問題から領域分割とノードを対応させて計算を行っている。予備計算により 1 領域 1 CMG
で割り当てた場合、1 領域 1 ノードの場合の 3 倍程度の速度になることを確認していることから、領域分割を適切に行えば、解析速度は 10 GFLOS 程度となり、解析時間は 12 時間程度となると予想される。なお、プリ処理に要する時間は角柱の CAD データ作成に必要な時間を除けば 6 分 14 秒である。ポスト処理に必要な時間はデータダウンロードに必要な時間を除けば約 1 分である。
角柱解析を行うことにより FFX の精度検証を実施し、ナビエ・ストークス解析及び実験結果と一致する結果が得られることを確認した。
車両解析
自動車まわりの流れ場解析として、本田技術研究所から提供された車両データ（空力解析用に実車のCAD データを修正したもの）を使用して、FFX および FFB による解析を実施した。
解析は車両速度 100 km/h で実施し、タイヤおよび地面は固定条件で実施した。FFB の解析は、解析要
素数 54 億であり、これまでに風洞実験結果との比較等を行い、空気抵抗、圧力変動等の結果が実験と良く一致することを確認済みのデータである。LBM で自動車のような複雑形状物体を解析する場合、格子サイ
ズを流れ場の渦スケール等を考慮して変更する階層型直交格子を用いることが多いが、今回は、FFX の基
本性能を確認するため、格子の大きさを解析領域ですべて同一としたシングルグリッドで解析した。なお、マルチスケール解析に評価については別途実施している。シングルグリッドの解析では格子点数をあまり増やすことはできないが、LBM 解析自体は解析速度が速く、メモリ使用量も小さいので、格子サイズ 6.25 mmの解析では格子数 47 億、格子サイズ 2 mm の場合で 1,274 億であった。1,274 億の解析を行うためのプリ処理にかかる時間は 48 分 37 秒であった。
Figure 2.2.5-15 に FFB で解析した車両まわりの流れを示す。時速 100 km で走行する車両周りの中央断面の流れは、リアウインドシールドの中央部分でわずかにはく離が見られるがほとんどはく離はなく、リアエ
ンドの止水領域も床下流れが下流側で合流することにより車両の半分程度の長さの領域にしか見られない。
ドアミラーや A ピラーなどの流れがはく離するような部位を除くと、車体表面の渦スケールは数十ミクロン程度と予想される。仮に車両代表寸法を 1 m とし、車体を平板と見立てれば、時速 100 km で走行する場合の乱流境界層の厚さは約 18 mm であり、無次元境界層距離 y+は 50 m と見積もることができる。今回実施した FFX での解析では格子解像度は 6.25 mm および 2 mm であるため、乱流モデルを用いているとはいえ、この格子解像度では解析することはできない、格子解像度 6.25 mm で解析可能なレイノルズ数は実車のレイノルズ数 Re0 のおおよそ 1/300 程度、2 mm の場合で 1/40 程度である。ただし、この場合、車両表- 63 -

面の境界層は層流であり、ルーフ上で遷移すると予想されるため、実際の車両まわりの流れとはやや異なることとなる。

Fig. 2.2.5-15 Flow pattern around an automobile (FFB, Re = Re0).
これらのことを勘案し、実レイノルズ数における解析と合わせてレイノルズ数を 1/300 および 1/50 とした解析を行った。ここでスマゴリンスキー定数は 0.15、音速は 340 m/s である。
Figure 2.2.5-16 に格子解像度 6.25 mm の結果を示す。Figure 2.2.5-15 とは異なり、ルーフ後端で流れがはく離していることがわかる。レイノルズ数を実レイノルズ数に設定しても格子解像度が不足していため、
低レイノルズ数の流れになっていることがわかる。Figure 2.2.5-17 に FFB で解析したレイノルズ数 1/200 の解析結果を示す。FFX の解析結果はこの流れ場に近く、設定した格子解像度（6.25 mm）では実レイノルズ
数の 1/200 から 1/300 程度の流れが再現されているものと考えられる。

- 64 -

(a) Re=Re0/300.

(b) Re = Re0.
Fig. 2.2.5-16 Flow fields around an automobile simulated by FFX.

Fig. 2.2.5-17 Flow field around an automobile simulated by FFB, Re=Re0/200.
Figure 2.2.5-18 に音場の解析結果を示す。フロントグリル部とフロントシールドの付け根部分から強い音が発生していることがわかる。車両下流には角柱解析と同様に流れ場による強い音場が見られるが、これは音として遠方に伝播する成分ではないことから車内騒音を予測する場合は、遠方まで伝播する音の寄与と流体成分の寄与を分離する必要がある。

- 65 -

Fig. 2.2.5-18 Sound field around an automobile simulated by FFX, Re=Re0/300.
Figure 2.2.5-19 に FFB で解析した実レイノルズ数の 1/50 の流れ場を示す。Figure 2.2.5-20 は FFX による解析結果（Re=Re0/40）である。両者が非常に近い流れ場になっていることがわかる。Figure 2.2.5-21
および 2.2.5-22 は実レイノルズ数の解析結果である。FFX による解析結果は実レイノルズ数に近いものであるがリアシールド後半で流れがはく離してしまい、FFB の結果とは車両下流側で異なるものとなることがわかる。マルチグリッドを適用し、解像度を上げればこの問題は改善されると予想される。

Fig. 2.2.5-19 Flow field around an automobile simulated by FFB, Re=Re0/50.

Fig. 2.2.5-20 Flow field around an automobile simulated by FFX, Re=Re0/40.

- 66 -

Fig. 2.2.5-21 Flow field around an automobile simulated by FFB, Re=Re0.

(a) Centerline.

(b) Around rear-view mirror.
Fig. 2.2.5-22 Flow field around an automobile simulated by FFX, Re=Re0.
Figure 2.2.5-23 に速度勾配テンソルの第 2 不変量の等値面を示す。微細な渦構造が解析されていることがわかる。格子ボルツマン法では直交格子を用いていること、本解析では格子解像度を一定としていることから車両下流領域まで微細な渦構造を解析することができていることがわかる。

- 67 -

Fig. 2.2.5-23 Iso surface of the second invariant of velocity gradient tensor around an automobile
simulated by FFX, Re=Re0/40.
Figure 2.2.5-24 に車両まわりの空力音場を示す。解像度 6.25 mm の場合と同様にセンターライン
上では、ボンネット、フロントシールドの付け根部分に強い音源が見られる。また、車両サイドの可視化結果からドアミラー、フロントホイールハウスに音源があることがわかる。Figure 2.2.5-25
はドアミラー周りの音場をより詳しく見るため、速度場の発散を示したものである。車体表面近傍には流れに起因する強い音場（圧力変動）が観察される。ドアミラーを取り巻くように円弧上に薄い輪が見られるが、これは遠方まで到達するいわゆる音波である。LBM 解析により音場と流れ場を計算することができることを確認することができた。

- 68 -

Fig. 2.2.5-24 Sound fields around an automobile simulated by FFX, Re=Re0/40.

Fig. 2.2.5-25 Distribution of divergence of velocity vector around an automobile simulated by
FFX, Re=Re0/40.
Figure 2.2.5-26 に車両空気抵抗の結果を示す。格子解像度 6.25 mm の場合はレイノルズ数 1/300 の結果、格子解像度 2 mm の場合は、レイノルズ数 1/40 の結果とほぼ等しい流体力が得られている。レイノル- 69 -

ズ数の最も大きなケース（実レイノルズ数）の解析では、設定したレイノルズ数が実レイノルズ数であっても得られる流れは格子解像度に依存した流れとなっているため、より低いレイノルズ数の空気抵抗に相当する結果となっていることがわかる。

Fig. 2.2.5-26 Reynolds number dependence of aerodynamic drag coefficient of vehicle andevaluation of accuracy of LBM.
解析時間は解像度 2 mm の場合で約 30 時間であった。解析速度は 15 GFLOPS/node であり、解析領
域を 1 CMG に割り振った場合は 45 GFLOPS/node であった。
より詳細な解析を実施するため、エンジンルームなどを含む車両 CAD データ（Figure 2.2.5-27）を使用した場合、プリ処理に必要な時間は 13 時間 32 分 15 秒であった。

Fig. 2.2.5-27 Preparation status of detailed LBM analysis including engine compartment,
wheels, suspension etc.
2-2) 分離解法を用いた解析
空力音を予測という観点からは、LBM のような空力音の直接解法は非常に有益であるが、自動車車内騒音のように車室内に音を伝達する部材であるウィンドシールドや車体壁が変動流れ場の内部にある場合、流- 70 -

れ場による圧力変動と音場による圧力変動が車内騒音（板の曲げ振動）にどのように寄与しているかを理解する必要があり、そのような観点からは分離解法が優れている。予測するという点だけで考えると音場と流れ場の圧力変動を分ける必要はないが、振動伝達機構を考える場合、両者の寄与を分離する必要がある。そこで、本研究では車内騒音に対する音場と流体力の影響を調べるため、以下の解析を実施した。
Figure 2.2.5-28 に示すように風洞ノズル下面と高さが一致するように厚さ 20 mm の鉄板を厚さ 15 cm のコンクリート製の升の上に設置して、自動車の車室内を模擬した。箱の内部には吸音材が詰められており、内部での音の反射を抑えた。鉄板の一部分は自動車の窓を模擬したアクリル板からなり、この窓部以外は、質量が大きく、剛性も高いことから振動が極力抑制されるようになっている。アクリル窓の上流側に矩形ブリックを置き、流れを変動させることにより、アクリル板の表面に流体圧力変動、空力音響変動を与え、アクリル板を介して車室内に伝播する音を計測する。この装置を模擬した解析を実施し、車内音の伝播解析を行った。
流体解析には FFB を使用し、FFB の解析結果より音源となる流体圧力変動の時間 2 回微分項を算出し前年度に開発した Acoustic Perturbation Equation（APE）法を用いて外部空力音を算出した。非圧縮流れ解析で得られた流体圧力変動と分離解法で得られた空力音の圧力変動を加振力として、アクリル板の振動を FEM 解析により求めた。

Fig. 2.2.5-28 Picture of experimental set up of wind tunnel for vehicle interior noise.
Figure 2.2.5-29 に流れ場と窓に作用する流体圧力の瞬時分布を示す。Figure 2.2.5-30 に角柱から放出された空力音を、APE を用いて求めた結果を示す。

- 71 -

Fig. 2.2.5-29 Distribution of Instantaneous flow around a square bump and the surface pressure
on the windshield.

Fig. 2.2.5-30 Sound field radiated form a square bump.
Figure 2.2.5-31 に振動解析結果を示す。黒い実線は実験により得られた振動スペクトルであり、流体圧力
変動及び空力音の両方の影響が合わさったもののである。赤い実線は流体加振力に起因した振動、青い
実線は音場に起因した振動である。この実験では比較的低い周波数の加振力が大きく、アクリル板に対する音場のコインシデンス周波数（2 kHz）より低い周波数領域での実験であるため、流体加振力の影響が大き
いことがわかる。

- 72 -

Fig. 2.2.5-31 Comparison of hydrodynamic and acoustic excitation forces on the vibration
spectrum of a flat plate.
このように同じ圧力変動であっても流体の圧力変動と空力音の圧力変動が板の振動に及ぼす影響は異なる。この違いは流体圧力変動と音響圧力波の伝播速度の違いにより両者の波数領域が異なること、パネルの曲げ振動の波数が周波数の平方根に比例するため、流体加振力と音響加振力の寄与が大きくなる波数が異なることに起因する。すなわち波数区間で見ると、流体加振力は低い周波数帯域でパネル振動の波数に接近し、音響加振力はより高い周波数帯域でパネル振動の波数に近づく。二つの物理量の波数が近く、波数区間での分布関数が似ているほど振動しやすくなることから、波数がどの周波数で一致するかが重要となる。また波数空間での分布関数は、振動と音響は（波であることから）複素波数空間の原点を中心とした円形状に広がるのに対して、流体の波数スペクトルは主流があることから一方向に偏ったものとなり、分布関数の形が異なることから加振力としての効率が悪い。これらを模式的に表したものが Figure 2.2.5-32 である。流れ場と音場の波数・周波数スペクトルを算出もしくはモデル化し、これとパネルの曲げ振動の波数・周
波数スペクトルの相関を求め、相関値が大きい場合は、加振力が効率よく振動となり、相関が小さければ振動が生じない。低マッハ数流れでは、流体加振力は空力音の加振力の 1,000 倍近いパワーを持つが、車内音への寄与は小さいと言われるのはこのためである。

- 73 -

Fig. 2.2.5-32 Simulation model for vehicle interior noise using modal joint acceptance.
このことを定式化したのが式(2.2.5-5)である。

𝐽𝑟 (𝝎)2 =̂ 𝑟 (𝒌,𝜔)|2 𝑑𝑘𝑥 𝑑𝑘𝑦∬ 𝛷𝑝𝑝 (𝒌,𝜔) |𝛷
4𝜋 2 𝑆0 𝐴2）(2.2.5-5)

板の曲げ振動の加速度はモーダルジョイントアクセプタンスを用いて以下のように記述することができる。

|𝛇𝒓̇ | =𝝎𝟐 𝑨𝟐 𝑱𝒓 (𝝎)𝟐 𝑺𝟎
𝒎𝒓 |𝒀𝒓 |𝟐(2.2.5-6)

この式から板から音への放射音を計算すると

𝟏
𝑰𝒓 = 𝝆𝟎 𝑪𝟎 𝟐
𝟒𝝅

̂𝒓 |𝟐
|∅
∬
𝒅𝒌𝒙 𝒅𝒌𝒚 |𝜻𝒓̇ |𝟐𝟐(|𝒌|/𝒌)𝟏−
𝒂𝒄𝒐(2.2.5-7)|𝒌|<𝒌𝒂𝒄𝒐となる。したがって、車両周りの流れ場や空力音が車内音として伝わる場合、流体圧力変動、空力音の波
数・周波数スペクトルが重要であり、それとパネル自体の曲げ振動の波数・周波数スペクトルの関係が重要となる。式(2.2.5-5)から式(2.2.5-7)の個々の変数の数値を数値計算で求め、それらをデータベース化し、機械学習により推定することで車内騒音を予測することが可能となると考えられる。あるいはこれらの変数をモデル化することによりサロゲートモデルを作成することが可能であると考えらる。
Figure 2.2.5-33 は一般的な車両の窓ガラスの物性値をもとにパネルの波数・周波数スペクトルを求め、乱流境界層と拡散音場をもとに外部加振力をモデル化して車内騒音を推定した結果である。周波数 1 kHz 付
近以上から外部空力音が卓越し、低周波数では流体加振力の寄与が大きいことがわかる。外部入力の波数・周波数スペクトルの分布関数のモデル化、パネルの振動モードの推定を行う必要があるが、これらをモデル化することで車内音の特性を評価することが可能であることがわかる。大規模解析と機械学習を組み合
わせ、サロゲートモデルを構築することが今後の課題である。

- 74 -

Fig. 2.2.5-33 Estimation results using surrogate models of vehicle interior noise withwavenumber frequency spectra.
[1] Y. H. Qian, D. D'Humires and P. Lallemand, Europhysics Letters, Volume 17 (1992), Number
6.
[2] D. d’Humires et al., “Multiple relaxation time lattice Boltzmann models in three dimensions”,
Phil. Trans. R. Soc. A 360, (2002), 437 451.
[3] 店橋護、岩瀬識、宮内敏雄、日本機械学会論文集 B 編 / 65 巻 638 号、pp. 3237-3243 (1999)。
[4] B. Stahl, B. Chopard, and J. Latt, "Measurements of wall shear stress with the lattice
Boltzmann method and staircase approximation of boundaries", Computers & Fluids 39(2010) 1625-1633.
[5] M’Hammed Bouzidi, Mouaouia Firdaouss Firdaouss and Pierre lallemard, "Momentum
transfer of Blotzman-lattice fluid with boundaries", Physics of Fluid 13 (2001), 3452.
(6) プロジェクトの総合的推進
プロジェクト全体の連携を密としつつ円滑に運営していくため、プロジェクトの推進や実証研究
テーマ間の連携のための会議等を適宜開催し、参画の協力機関・連携機関との連携・調整にあたる。
特に、プロジェクト全体の進捗状況を確認し、計画の合理化の検討等を行うなど、プロジェクトの効果的・効率的推進に資する取組みを実施する。また、プロジェクトで得られた成果については、
可能な限り積極的に公開して今後の展開に資するとともに、ものづくり産業での早期戦力化を支援する。
令和 3 年度は、令和 2 年度に引き続き、上記を順次実施しながら、本課題の有効な研究活動のための支援を実施し、研究開発全体を円滑に実施するための調整を行った。そして、得られた成果については、HPCI ものづくり統合ワークショップや「富岳」流体予測革新プロジェクトシンポジウムを開催し、そこで公開した。

- 75 -

以下に代表的な活動を示す。
1) 「『富岳』成果創出加速プログラム」第 5 回 HPC ものづくり統合ワークショップ
日時：令和 3 年 10 月 22 日
場所：Webex Events
参加者：192 名
2) 第 2 回「富岳」流体予測革新プロジェクトシンポジウム
日時：令和 4 年 3 月 9 日
場所：東京大学生産技術研究所コンベンションホール、および、Webex Events
参加者：304 名（含、現地参加 10 名）
２－３．活動（研究会の活動等）(1) アウトリーチ活動
年月日

行

事

R3.6.11~12

東京大学生産技術研究所キャンパス公開

R3.9.2~3

国際フロンティア産業メッセ 2021

R3.10.22「『富岳』成果創出加速プログラム」第 5 回 HPC ものづくり統合ワークショップ

R3.10.29

RIST 成果報告会

R3.11.15~18SC21

R4.3.9

第 2 回「富岳」流体予測革新プロジェクトシンポジウム

R4.3.14~15「富岳」成果創出加速プログラム 研究交流会

R4.3.29「富岳」成果創出加速プログラムシンポジウム(2) 「富岳」流体予測革新プロジェクト推進会議
年月日

行

事

R3.5.26

第 5 回「富岳」流体予測革新プロジェクト推進会議

R3.9.29

第 6 回「富岳」流体予測革新プロジェクト推進会議

R3.12.3

第 7 回「富岳」流体予測革新プロジェクト推進会議

R4.3.30

第 8 回「富岳」流体予測革新プロジェクト推進会議(3) 「富岳」高性能シミュレーション技術交流会
年月日
R4.1.21

行事

第 2 回「富岳」高性能シミュレーション技術交流会- 76 -

(4) アプリケーション性能向上ワーキンググループ
年月日

行

事

R3.4.20

第 13 回アプリケーション性能向上ワーキンググループ

R3.5.20

第 14 回アプリケーション性能向上ワーキンググループ

R3.6.15

第 15 回アプリケーション性能向上ワーキンググループ

R3.7.14

第 16 回アプリケーション性能向上ワーキンググループ

R3.8.24

第 17 回アプリケーション性能向上ワーキンググループ

R3.9.22

第 18 回アプリケーション性能向上ワーキンググループ

R3.10.27

第 19 回アプリケーション性能向上ワーキンググループ

R3.12.3

第 20 回アプリケーション性能向上ワーキンググループ

R4.1.14

第 21 回アプリケーション性能向上ワーキンググループ

R4.2.16

第 22 回アプリケーション性能向上ワーキンググループ

R4.3.25

第 23 回アプリケーション性能向上ワーキンググループ

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

担当機関

担当責任者

(1) 数値曳航水槽の実現と省エネデバイスによる推進効率の向上(2) 細隙部を含めた多段遠心ポンプの内部流
れ の Wall-Resolved
LES(3) 圧縮機サージの直接
解析

国立大学法人東京大学 生産技術研究所

革新的シミュレーション研究センターセンター長／教授 加藤千幸

国立大学法人九州大学 大学院
工学研究院
国立大学法人東京大学 生産技術研究所
国立大学法人九州大学 大学院
工学研究院
国立大学法人岩手大学 理工学部(4) リアルワールド自動 国立大学法人神戸大学 大学院
車空力性能の予測システム情報学研究科
国立研究開発法人理化学研究
所 計算科学研究センター
国立大学法人山梨大学 大学院
総合研究部- 77 -

教授 渡邉聡革新的シミュレーション研究センターセンター長／教授 加藤千幸
教授 古川雅人
准教授 山田和豊
教授 坪倉誠チームリーダー 坪倉誠
教授 岡澤重信

(5) リアルワールド自動 国立大学法人豊橋技術科学大
車空力音予測学 大学院工学研究科
国立大学法人神戸大学 大学院
システム情報学研究科
国立研究開発法人理化学研究
所計算科学研究センター
国立大学法人東京大学 生産技術研究所(6) プロジェクトの総合 国立大学法人東京大学 生産技
的推進
術研究所

- 78 -

教授 飯田明由
教授 坪倉誠チームリーダー 坪倉誠
革新的シミュレーション研究センターセンター長／教授 加藤千幸
革新的シミュレーション研究センターセンター長／教授 加藤千幸

別添１ 学会等発表実績
１．学会等における口頭・ポスター発表
発表した成果（発表題目、
口頭・ポスター発表の別）「富岳」の時代におけるものづくりとシミュレーション〜期待と課題〜
富岳」の時代の HPC とものづくり
The role of a professionalsocietyofmechanicalengineers after and/or with
Covid-19 pandemic
High-Performance
ComputingandIts
Industrial Applications inFugaku Era
多重プロセス型モデルを用いた Clark Y-11.7%翼まわりのキャビテーション流れの解析
多重プロセス型モデルを用いた多段遠心ポンプ最終段モデルのキャビテーション解析
多重プロセス型モデルを適用した単独翼周りのキャビテー
ションの Wall-Resolved LES

発表した場所（学会等名）

発表し
た時期

国内・
外の別

加藤千幸

燃焼システム用次世代 CAE
コンソーシアム総会

2021 年
5月

国内

加藤千幸

2021 年
10 月
2021 年
10 月

国内

Chisachi Kato

VINAS
OnlineUsers
Conference 2021
69th
PSME
National
Conference （ フ ィ リ ピ ン 機
械学会年次大会）

Chisachi Kato

IAHR-Asia 2021

2021 年
11 月

国外

高峯大輝，津田伸
一，渡邉聡

日本機械学会第 99 期流体工学部門講演会2021 年
11 月

国内

山田晃久，高峯大輝，津田伸一，渡邉聡
津田伸一，高峯大輝，渡邉聡

日本機械学会第 99 期流体工学部門講演会2021 年
11 月

国内

キャビテーションに関するシンポジウム（第 20 回）

2021 年
12 月

国内発表した場所（学会誌・雑誌等名）

発表し
た時期

国内・
外の別

2022 年
3月
2021 年
9月

国内

発表者氏名

国外

２．学会誌・雑誌等における論文掲載
掲載した論文（発表題目）
圧縮機サージングの非定常三次元数値解析
A one-sided direct forcing
immersed boundary method
using moving least squares

発表者氏名
古川雅人

ターボ機械

Rahul
Bale,
Amneet Pal Singh
Bhalla, Boyce E.
Griffith,
MakotoTsubokura

Journal of Computational
Physics- 79 -

国外