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

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

千幸

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

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

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

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

-2-

（協力機関：国立大学法人豊橋技術科学大学、国立大学法人神戸大学、国立研究開発法人理化学研究所）(6) プロジェクトの総合的推進
プロジェクト全体の連携を密としつつ円滑に運営していくため、プロジェクトの推進や実証研究
テーマ間の連携のための会議等を適宜開催し、参画の協力機関・連携機関との連携・調整にあたる。
特に、プロジェクト全体の進捗状況を確認し、計画の合理化の検討等を行うなど、プロジェクトの効果的・効率的推進に資する取組みを実施する。また、プロジェクトで得られた成果については、
可能な限り積極的に公開して今後の展開に資するとともに、ものづくり産業での早期戦力化を支援する。
令和 4 年度は、令和 3 年度に引き続き、上記を順次実施しながら、本課題の有効な研究活動のための支援を実施し、研究開発全体を円滑に実施するための調整を行う。得られた成果については、「富岳」流体予測革新プロジェクトシンポジウム等を開催し、そこで公開する。
２－２．実施内容（成果）(1) 数値曳航水槽の実現と省エネデバイスによる推進効率の向上(a) 研究成果の要旨
本実証研究テーマでは、「富岳」とポスト「京」重点課題⑧で研究開発した LES 解析ソフトウェア
である FrontFlow/blue（FFB）を用いて、300 億要素から最大 1,200 億要素の Wall-Resolved LES
計算により、船のまわりの乱流境界層を完全に解像する予測を実現し、曳航水槽試験を完全に数値
シミュレーションに代替えし得ることを実証するとともに、船の推進効率の向上のための省エネデ
バイスの動作メカニズムを解明し、さらなる高効率な省エネデバイスの開発に貢献する。
令和 4 年度は、数値シミュレーションが、曳航水槽試験を代替えし得ることを実証するとともに、省エネデバイスの動作メカニズムの解明に貢献できることを検証した。具体的には、数値曳航水槽の実現のサブテーマに関しては、数値シミュレーションによる予測結果と国内複数社の曳航水槽試
験結果との比較検証を行い、両者が良好に一致することを確認した。省エネデバイスによる推進効率の向上のサブテーマに関しては、舵フィンは単独で用いると性能が悪くなるが、舵バルブと一緒に用いると性能が良くなるといったケースについて、その流体力学的メカニズムを調査した。その結果、プロペラシャフト中心から左舷のやや下側に放出されるハブ渦の影響により、その効率が低下していることがわかった。具体的には、舵フィンは単独で用いると、ハブ渦の負圧によって効果を発揮しないことがわかった。また、ハブ渦が大幅に拡散される舵バルブと一緒に用いると舵フィ
ンはその効果を発揮し、性能が良くなることもわかった。
In this research, we will show that towing tank tests can be completely replaced by numericalsimulations by resolving turbulent boundary layer on the hull surface, and we will contribute to
the development of even more efficient energy-saving devices for a hull by clarifying theirmechanism. To this end, we will use FrontFlow/blue (FFB), which has been developed for-3-

supercomputer “Fugaku”, for performing Wall-Resolved LES with 30~120 billion computational
grids.
In this fiscal year, we have demonstrated that the numerical simulation has the capability of
replacing towing tank tests and contribute to clarifying the hydrodynamical mechanism ofenergy-saving devices. Specifically, regarding the sub-theme of realizing a numerical towing tank,
we compared and verified the results of numerical simulations with the results of towing tanktests conducted by several domestic shipbuilding companies and confirmed that thecomputational results were in good agreement with the experimental results. Regarding the subtheme of improving propulsion efficiency through energy-saving devices, we investigated thehydrodynamical mechanism of the cases where the rudder fins by themselves did not perform as
intended but they performed as expected when used with a rudder bulb. As a result, it was foundthat the rudder fins are ineffective when used alone due to the negative pressure of the hub vortex.
Furthermore, it was observed that when the rudder fins are used with a rudder bulb, whichsignificantly scatters the hub vortex, the rudder fins are more effective and perform better.
(b) 本研究の内容
数値曳航水槽の実現のサブテーマに関して、曳航水槽試験と計算の全抵抗値（CT）を比較したものを Fig. 2.2.1-1 に示す。実験と計算はそれぞれ、船長 4.4 m と 8.1 m の相似模型船を対象とした。
図中実線と点線で示しているものは摩擦抵抗係数を使った全抵抗値の、4.4 m と 8.1 m の相似模型
船を対象としたそれぞれの実験値からの外挿ライン（3 次元外挿法）である。使っている摩擦抵抗
係数は ITTC（International Towing Tank Conference、国際試験水槽会議）で推奨されている式と、Schoenherr の式の 2 種類をプロットしている。8.1 m の計算はほぼ実験と一致しているが、4.4m の計算は 4.4 m の実験より若干大きい値が得られ、8.1 m の実験値を使った Schoenherr の外挿
ラインと一致し、ITTC の外挿ラインより小さくなった。実験で 4.4 m と 8.1 m の外挿ラインが一致しない理由は 8.1 m では船体表面全面で乱流になっているが、4.4 m では層流領域が残っていると考えられる。これまでは全面乱流となるような高レイノルズ数のケースを想定して計算手順を構築していたため、乱流促進スタッドのまわりの格子は乱流になりやすいような簡易的な計算格子を使っており、今回の計算でも 4.4 m の計算結果の可視化からは、船体全面で乱流であることが確認されており、8.1 m の実験値を使った Schoenherr の外挿ラインと一致しているのはこのためである。そこで、本研究では改めて、乱流促進スタッドのまわりの格子を密にした計算と従来の格子の比較を行った。スタッド付近の渦の様子を Fig. 2.2.1-2 に示す。この計算からは局所細密化した計算ではスタッド直後に層流領域が残っていることがわかった。しかしながら、スタッド周辺以外への影響はほとんどないことも同時に確認され、船体全体の抵抗値は、格子の細密化によってほとん
ど影響しない（0.1 %程度）ことがわかった。今回の調査により、層流領域が残ると予想される対象については従来考えられていた格子密度、すなわち(x+,y+,z+)=(30,3,10)より密な格子が必要であるという知見が得られた。

-4-

Fig. 2.2.1-1 Comparison between experiment and calculation of total resistance (CT) for
similarity model ship of 4.4 m and 8.1 m.

Fig. 2.2.1-2 Computational grid (upper panels) which are used in previous studies (left) andpresent studies (right) and calculated vortices (lower panels).
省エネデバイスによる推進効率向上に関するサブテーマについては、Fig. 2.2.1-3 に示すように 4種類の省エネ付加物を対象とした。左から付加物が何もついていないもの、舵にバルブがついたもの、舵にフィンがついたもの、舵にバルブとフィンがついたものである。舵にフィンがついたもの
は、何もついていないものに比べて推進効率が悪いことが実験的に確認されていた一方で、バルブとフィンを併用することによりバルブ単体より効果が出ることも実験的には確認されていた。本研-5-

究では効率的に省エネデバイスの設計が行えるように、計算で流場を確認し、フィン単体で効果が出なかった理由を調べた。その結果、主な原因はハブ渦であることがわかった。Fig. 2.2.1-4 に 4 状
態のハブ渦の様子を示す。この図ではハブ渦に注目するため、プロペラは可視化していない。Fig.
2.2.1-4(a)の舵に何もついていない付加物なしの状態ではハブ渦はプロペラシャフト中心から左舷のやや下側に放出されていることがわかる。この理由については船体形状の影響から船尾では上昇流があり、ハブ渦はその回転方向（後ろから見て時計回り）から左向きの力を受け、その後、舵付近を通過するが、左舷舵付近ではプロペラ回転流によって下側に負圧があり、下側に引き寄せられたと考えられる。Fig. 2.2.1-4(c)の舵フィンはプロペラ回転流を回収するのが狙いであり、左舷側下面
では正圧を期待しているが、ハブ渦によって負圧になり、左舷側につけたフィンが有効に働いていないことがわかった。一方で、Fig. 2.2.1-4(b)の舵バルブがついている状態ではハブ渦が拡散していることがわかる。よって、Fig. 2.2.1-4(d)の舵バルブと舵フィンがついた状態では舵バルブによるハブ渦の拡散によって舵フィンが有効に働いていることがわかった。以上が、舵フィン単体と舵バル
ブと舵フィンを同時に装着したケースの違いであるが、舵フィン単体でも効果をより有効に働かせるためには、ハブ渦を弱めることが有効であることがわかったため、プロペラハブにハブ渦を弱めるための既存の省エネデバイスを追加するか、ハブ後端と舵の距離を調整することにより、舵フィンの改善が期待できるという知見が得られた。

Fig. 2.2.1-3 Configuration of energy-saving devices (ESD).

-6-

Fig. 2.2.1-4 Hub vortex structure and pressure distribution at stern. (a) without energy-saving
devices, (b) with rudder bulb, (c) with rudder fins, and, (d) with rudder bulb and rudder
fins.
(2) 細隙部を含めた多段遠心ポンプの内部流れの Wall-Resolved LES(a) 研究成果の要旨
本実証研究テーマでは、「富岳」を利用して、30～100 μm の最小渦スケールまで計算格子により直接解析する、Wall-Resolved LES を実施し、内部流れやその結果として決まる水力性能（全揚程・水力トルク・水力効率）の完全な予測を実現するとともに、細隙部内部流れの挙動やそれが性
能や信頼性に与える影響を明らかにし、ポンプ設計の高度化に貢献する。さらに、上記の 1/100 程
度の計算格子を用いた Wall-Modeled LES も実施し、Wall-Resolved LES の結果と計算精度、計算コストなどを比較することにより、Wall-Modeled LES の実用化を図る。
令和 4 年度は、細隙部内部流れの挙動やそれが性能や信頼性に与える影響を明らかにし、ポンプ設計の高度化に貢献できることを検証した。その結果、ポンプの主流部に加えて隙間幅 0.18 mm 程度の細隙部における縦渦も捉える LES 解析により、遠心ポンプの性能および流体力を定量的に予測できる見込みを得た。また、ポンプの設計流量ならびに部分流量（20 %設計流量）において、壁面近傍の詳細な渦構造も含めて流れ場を明らかにした。部分流量ではポンプ主流と羽根車とケーシング間の間隙部の非定常な干渉による軸スラスト力の変動を捉えることができた。このことによって、ポンプ軸系の設計の高度を図ることが可能になり、信頼性の向上に寄与できる。本研究対象である
遠心ポンプはプラントなどで使用される実機ポンプ形状相当であり、本研究結果は直接産業界で参照され、ポンプの設計・開発に活用することが可能である。

-7-

In this research, we will realize highly-accurate prediction of the internal flow of a pump andits performance, clarify the flow mechanism in the leakage passages and its influence on thepump performance, and finally contribute to the improvements of the pump design. To this end,
we will perform the wall-resolved LES of the internal flows of the pump where dynamics of thesmall vortices with the scale of 30~100 μm are directly computed by using “Fugaku”. We will also
perform wall-modeled LES, which requires 1/100 less computational grids than those for the wallresolved LES. We will realize practical use of the wall-modeled LES by comparing the accuracy
and cost of the computations between the wall-model LES and the wall-resolved LES.
In this fiscal year, we clarified the behavior of the flow inside the small gaps in the pump suchas the radial clearances at the liner ring and the balance piston, by using LES analysis which
captured longitudinal vortices in the narrow gaps with a gap width of about 0.18 mm. These
results verified that it can contribute to the advancement of pump design and expected tocontribute to a quantitative prediction of the performance and hydrodynamic force of centrifugal
pumps. In addition, we clarified the flow field including the detailed vortical structures near thewall at the design flow rate and a partial flow rate (20 % design flow rate) of the pump. At thepartial flow rate, the axial thrust force fluctuation due to unsteady interference between the main
flow of the pump and the gap flow between the impeller and the casing was detected. The
centrifugal pump, which is the subject of this research, has a shape equivalent to the actual pumpused in plants, etc., and the results of this research can be directly referenced in the industrials,
and can be utilized in the design and development of pumps.
(b) 本研究の内容
1) 遠心ポンプ内部流れ解析の計算モデル
Fig. 2.2.2-1 に遠心ポンプ内部流れ解析の計算モデルを示す。Fig. 2.2.2-1 に示すとおり、本 LES（Large Eddy Simulation）解析では軸スラストを評価するため、入口配管、羽根車、ディフュー
ザ、出口配管に加え、front/back chamber、balance piston、return pipe といった細隙部を考慮した計算モデルになっている。計算領域は Domain 1 から Domain 4 の 4 つの独立の計算領域から構成され、これらをオーバーセット法により連結させた。
LES 解 析 ソ ル バ ー は 「 FrontFlow/blue 」 v67.01 (FFB ） を 用 い 、 SGS モ デ ル は Dynamic
Smagorinsky モデルを適用した。
本遠心ポンプのベースの計算格子の要素数はおよそ 8,600 万である。FFB は計算実行時に計算格子を細分化する機能（refine 機能）をサポートしており、これを用いると、1 回 refine で 6.9 億要
素、2 回 refine で 55 億要素、3 回 refine で約 440 億要素の計算を実行することができる。本研究では 3 回 refine 時に 440 億要素で細隙部も含め乱流境界層を解像するものとして格子規模を計画したが、計算の安定性などの課題により本研究期間での実施には至らなかった。以降、8,600 万要素の計算を R0、6.9 億要素の計算を R1、55 億要素の計算を R2 と称する。ポンプ性能特性は 120 %設
計流量（設計流量（QD）に対する流量（Q）の比 Q/QD=1.2）から 20 %設計流量（Q/QD=0.2）まで
20 %流量刻みで解析した。計算は 100 %流量（Q/QD=1.0）を除き、R0 では羽根車 10 回転分、R1-8-

では羽根車 10 回転分、R2 では羽根車 5 回転分を行った。Q/QD=1.0 では R0 では羽根車 50 回転分、R1 では羽根車 20 回転分、R2 では羽根車 10 回転分を行った。これらの計算は「富岳」4,096 ノードを用いて実施した。

Fig. 2.2.2-1 Computational model of LES of pump internal-flow.
2) FrontFlow/blue（FFB）の開発
細隙部を含めたポンプ内部流れの LES 解析を「富岳」で効率的に実施するため、流体解析コード
FrontFlow/blue（FFB）の並列効率を向上させるための領域分割ツールの開発およびその効果の検証を実施した。
「富岳」における FFB の単体性能は現状 160 GFLOPS/ノード（ブーストモード）
であるが、1,000 ノード使用時の並列性能は 50 %程度である（Fig. 2.2.2-2）。この並列性能は、オーバーセット（以下、OS）計算を含まない場合であり、OS 計算を含む場合は、さらに数分の一程度まで性能が低下している。
「富岳」における FFB の実証計算の効率をあげるため、OS 計算の並列性能を向上させる必要があり、これを実現するため領域分割方法の改良を実施した。FFB では領域分割の際に Recursive Graph Bisection（以下、RGB 法）を用いている。RGB 法では要素の接続関係をグラフ化しこれを分割するが、グラフには OS の接続関係が考慮されていないため、隣接通信と比較して、OS 通信のコストが上がってしまうという課題があった。これを解消するため、グ
ラフ作成時に、OS データの親子を結びつける要素（ブリッジ要素）を加える改良をした（Fig.
2.2.2-3）。

-9-

Fig. 2.2.2-2 Sustained performance of FFB in weak-scale benchmark tests on Fugaku.

Fig. 2.2.2-3 Dummy elements (Bridge-elements) to represent overset connection when making
graph.
OS を含むテストとして、船舶用プロペラまわりの流れ計算を行った。計算格子の要素数は全体で1,720 万要素であり、プロペラまわりの領域は 990 万要素、外部領域は 730 万要素である。
「富岳」の 96 ノードを用いて、隣接通信、OS 通信の通信時間および全体の計算時間を比較した。Fig.
2.2.2-4 に通信パターンの比較を示す。従来方法での OS 通信パターンには特定の領域に OS 通信が集中しており、これが OS 通信の時間が増大し並列性能が低下する要因であった。グラフを一体化した新たな方法では、OS 通信の偏りが解消されていることがわかる。一方、隣接間通信について
は、新たな分割方法によるパターンのほうがばらついていることより、通信コストがあがっている
ことが予想される。Table 2.2.2-1 に測定した通信時間（隣接通信・OS 通信）および全体の計算時間を示す。グラフの一体化により、隣接通信は悪化した（0.088 s→0.151 s）ものの、OS 通信は期
待どおり高速化され（0.781 s→0.272 s）、全体の計算時間は 40 %程度短縮された（1.230 s→0.784s）。グラフの一体化により、オーバーセット計算が高速化できることが確認された。

- 10 -

(a) Previous method.

(b) Present method.
Fig. 2.2.2-4 Comparison of communication patterns of adjacent communication (left) andoverset communication (right).
Table 2.2.2-1 Comparison of communication time and total calculation time.
分割方法
従来方法
グラフの一体化

通信時間 [sec/step]
隣接通信
OS 通信
0.088
0.781
0.151
0.272

計算時間
[sec/step]
1.230
0.784

上記の技術を実装した FFB を用いて、細隙部を含めたポンプ内部流れの LES 解析を実施した。
3) ポンプ性能特性計算結果
Fig. 2.2.2-5 にポンプ性能の比較結果を全揚程（上）と軸動力（下）の流量特性について示す。
LES は Q/QD>0.4 で全揚程（上）を過大評価しており、流量が大きいほど LES と実験結果の差が大きい。また、R0 より R1、R2 の方が実験との差が大きくなっている。設計流量における 2 回 refine
時の解析結果は実験値よりもポンプ揚程を 5 %程度過大に評価し、また軸動力も 4 %程度過大に評価している。

- 11 -

Fig. 2.2.2-5 Comparisons of pump performance (upper panel: total pressure, lower panel: pump
power).
本ポンプでは、Q/QD<0.5 の全揚程の右上がり特性が発生するが、LES はこの右上がり特性を捉えることができており、右上がり特性の原因となる内部流れを適切に予測していると考えられる。
Fig. 2.2.2-6 には refine なし（R0）、1 回 refine 時（R1）ならびに 2 回 refine 時（R2）の設計流量における羽根車内の壁面近傍の渦構造を示す。ここで HUB は羽根車のハブ面、PS と SS はそれぞれ羽根の圧力面と負圧面を表す。羽根車が作用する仕事を精度良く予測するためには各羽根面上の乱流境界層を適切に解像する必要があるが、図が示す壁面近傍の渦構造を見ると、R2 では負圧面および圧力面近傍を除くハブ面の境界層には縦渦が発達していることが捉えられ、乱流に遷移していると判断されるが、圧力面にはそのような縦渦は確認できず、乱流に遷移していない様相を示している。

- 12 -

圧力面には縦渦
は確認できない
圧力面の境界層が薄く、ブロッケージが小さくなるため羽根のオイラーヘッドを過大評価している、あるいは損失を過小評価
している可能性
ポンプ揚程の過
大評価の要因Fig. 2.2.2-6 Vortex structures on the impeller walls.
上記の分析結果より、2 回 refine の格子解像度で軸動力・揚程を過大評価する理由として、ブロッ
ケージの過小評価による軸動力の過大評価（流路断面積が大きく出口速度が小さくなるため）、損
失の過小評価による揚程の過大評価が可能性として推察される。この点については今後、格子解像度の影響度に加え、羽根車に流入する乱れおよび偏流の影響を調査していく計画である。
4) ポンプ軸スラスト力予測結果
Fig. 2.2.2-7 に軸スラストの比較を示す。実験結果では大流量側から Q/QD=0.8 まではスラスト力は増加し、さらに流量を減じていくと Q/QD=0.4 まではスラスト力が減少したのち締切までは
Q/QD=0.2 でやや増加し再度減少するといった複雑な特性を示す。一方で解析結果は Q/QD=0.8 までのスラスト力増加傾向は捉えておらず、流量に対してフラットな結果となっている。この要因としては、間隙部の乱流境界層の解像度不足によって間隙部の速度分布および圧力分布の予測精度が不足していることが考えられる。

- 13 -

Fig. 2.2.2-7 Comparisons of axial thrust.
5) ポンプ細隙部の流れ
羽根車に作用する軸スラスト力の予測には Front chamber や Rear chamber ならびにライナーリ
ングやバランスピストン部の半径方向隙間などの細隙部の流れを精度良く予測する必要がある。
Fig. 2.2.2-8 は Q/QD=1.0、Q/QD=0.2 での Front chamber と Rear chamber 内の半径方向静圧分布の予測値を実験結果と比較し示したものである。Q/QD=1.0 では LES 解析の結果は半径方向の圧力勾配が大きく、内径側では実験値との差が大きくなっている。Q/QD=0.2 では Front chamber、Rear chamber ともに LES 解析の結果が小さいことがわかる。本研究では 1)で述べたように、格子
リファインを 3 回施すことで細隙部を含めた乱流境界層を解像するものとして格子規模を計画したものであり、解像度の改善によりこれらの予測精度が向上するものと考えている。

Fig. 2.2.2-8 Comparison of pressure distributions in the front and rear chamber.

- 14 -

Fig. 2.2.2-9 には前面ライナーリング隙間部の軸方向速度と周方向速度の分布をライナーリング入口側（z=-0.835）、ライナーリング中央部（z=-0.805）、ライナーリング出口部（z=-0.775）について示した。1 回リファインの場合の(a)では、ライナー入口部の軸方向速度はケーシング側の(rRi)/h=0.8～1.0 ではライナー入口の面取り部を流れが急激に曲がる影響で軸方向速度が小さく、(rRi)/h=0.1～0.7 では軸方向速度が一定である。流れがライナーリング中央部と出口部に進むと軸方向速度分布は層流的な放物線状の軸方向速度分布を示す。一方で 2 回リファインの(b)ではライナー
リング中央部では(r-Ri)/h=0.4～0.8 の領域で速度分布が 1 回リファイン時よりも平坦になり、ケーシング側の速度勾配が大きくなっていることがわかる。また、周方向速度分布についてもケーシン
グ側の速度勾配が大きくなっていることがわかる。

Fig. 2.2.2-9 Velocity distribution in the liner ring clearance.
Fig. 2.2.2-10 と Fig. 2.2.2-11 にはライナーリング部のケーシング側（静止壁側）と羽根車側（回転壁側）ならびにライナーリング隙間（幅 0.18 mm 程度）入口側、中央部、出口側の渦度分布を 1
回リファイン時、2 回リファイン時について示す。1 回リファイン時ではライナーリング入口側で縦渦の発生を示す、渦度の大きい領域が周方向に縞状に分布しているが、流れがライナーリング中央部、出口部に向かうにつれ渦が消滅していることがわかる。2 回リファイン時にはケーシング側では縦渦構造が隙間出口まで維持されているが、羽根車側では渦度の変動は残るものの、縦渦構造は消滅していることがわかる。
ライナーリング隙間入口の速度が約 15 m/s であることから、隙間部壁面に発生する縦渦の渦直径は 0.0492 mm、渦間隔は 0.246 mm と予測される。2 回リファイン時（R2）での結果はそれぞれ約
0.38 mm、0.84 mm であり、渦間隔で考えると少なくとも約 3.5 倍の格子数が必要と考えられる。

- 15 -

Fig. 2.2.2-10 Vorticity distribution near the outer wall, inner wall and on the cross sections ofliner ring clearance - R1 case.

Fig. 2.2.2-11 Vorticity distribution near the outer wall, inner wall and on the cross sections ofliner ring clearance - R2 case.

- 16 -

6) ポンプ内部流れの変化
Fig. 2.2.2-12 と Fig. 2.2.2-13 には Q/QD=1.0、Q/QD=0.2 での半径方向速度成分のコンターを R0、
R1、R2 に対して示した。Q/QD=1.0 では流路内に大きな逆流領域は存在せず、非定常性の強い領域
としてはディフューザスロート部からディフューザ負圧面に沿って発達する半径方向速度が正の領域とディフューザ後縁ウェーク部の半径方向速度が負の領域が交互に放出されているディフューザ
ベーン後縁部下流の流れがあげられる。羽根車内では負圧面側の速度の大きい領域と圧力面側の速度の小さい領域が明確に分かれ、その境界では速度の乱れが大きい。また、負圧面の羽根車出口部
では速度分布の乱れが増大している。ディフューザスロート部では流れは加速され、羽根車がスロ
ート部を通過するタイミングで加速が大きくなるなどの非定常干渉による流れの変動も確認できる。
ベーン前縁部ではブロッケージによる局所的な強い半径方向内向き流れが発生していることがわかる。格子リファインによりこれらの流れの特徴に変化はないが、速度の乱れの空間スケールは格子規模に応じて小さくなっていることがわかる。

Fig. 2.2.2-12 Radial-velocity distribution on the mid-span plane of impeller and diffuser at
Q/QD=1.0.
Fig. 2.2.2-13 に示した Q/QD=0.2 では Q/QD=1.0 と比較して流れの様相は大きく変化し、羽根車とディフューザ間の領域で多くの渦が発生していることがわかる。また、紫色で示されるディフュー
ザから羽根車に逆流する領域が大きく発達するが、その領域は羽根車とディフューザの相対位置に応じて周方向に偏在し、羽根車の回転とともに発生位置は変化する。

- 17 -

Fig. 2.2.2-13 Radial-velocity distribution on the mid-span plane of impeller and diffuser at
Q/QD=0.2.
Fig. 2.2.2-14 には Q/QD=1.0 と Q/QD=0.2 での羽根車とディフューザ間の隙間部（Gap A）を通じ
Front chamber 内に漏れ出る流れと Front chamber 内流れの干渉の様子を示した。Q/QD=1.0 では
ディフューザシュラウド壁面近くの速度ベクトルを見るとディフューザ部では逆流は発生しておらず、ディフューザ前縁付近では Gap A を通じ Front chamber から漏れ出る流れが発生するものの、漏れ流れと Front chamber 内流れの干渉は小さく、Front chamber 内の流れはスムースである。一方、Q/QD=0.2 ではディフューザ部で発達する逆流により Gap A から Front chamber への漏れ流れが強くなり Front chamber 内に渦を発生させる。

Fig. 2.2.2-14 Flow field showing leakage flow from the main flow to front chamber through the“gap A” and the effects on the flow in the front chamber.

- 18 -

Fig. 2.2.2-15 には実験で確認された Q/QD=1.0 と Q/QD=0.2 での軸スラスト力変動の波形と振幅ス、ペクトルを示す。軸スラスト力変動の横軸は羽根車回転周期 TN [s]（1,600 min-1 で TN=0.0375 s）
振幅スペクトルの横軸は変動周波数を羽根車回転周波数 fn（26.7 Hz）との比で表したものである。
これらを見ると Q/QD=0.2 では Q/QD=1.0 と比べ、軸スラスト変動には羽根車回転周波数以下の周波数の変動が増加し、振幅も二乗平均値で 1.9 倍ほど大きくなっていることがわかる。このような
Q/QD=0.2 における羽根車に作用する流体力変動の増加は LES 解析結果に示される羽根車とディフ
ューザ間での渦の増加やディフューザから羽根車への逆流領域の発達、Gap A からの漏れ流れとFront chamber 内流れの干渉の増加に起因すると考えられる。

Fig. 2.2.2-15 Measured axial thrust fluctuations and FFT analysis results.
7) まとめ本実証研究テーマでは、「富岳」により高品質の解析格子を用いた計算を実施し、細隙部を含めた
遠心ポンプの性能および流体力を定量的に予測できる見通しを得ることができた。具体的には、格子のリファインにより、ポンプの主流部に加えてライナーリングにおける隙間幅 0.18 mm 程度の細隙部における縦渦も捉える LES 解析により、遠心ポンプの性能および流体力を定量的に予測できる見込みを得た。
ポンプ性能の予測に加え、ポンプの設計流量ならびに部分流量（20 %設計流量）において、壁面
近傍の詳細な渦構造も含めて流れ場を明らかにした。部分流量ではポンプ主流と羽根車とケーシン
グ間の間隙部の非定常な干渉による軸スラスト力の変動を捉えることができた。

- 19 -

本研究対象である遠心ポンプはプラントなどで使用される実機ポンプ形状相当であり、本研究結
果は直接産業界で参照され、ポンプの設計・開発に活用することが可能である。
(3) 圧縮機サージの直接解析(a) 研究成果の要旨
本実証研究テーマでは、プラント、およびガスタービンなどの実機に用いられる遷音速圧縮機を対象として、圧縮機本体だけでなく、それが設置される配管系まで含めたシステム全体を計算領域とし、圧縮機羽根車の失速現象を再現できるほど短い時間刻みを設定するとともに、システム全体にわたる長周期の流体振動現象を捉え得るほど膨大な時間ステップ数にわたって DES（Detached
Eddy Simulation）による非定常三次元流動解析を「富岳」上で実施することにより、圧縮機サー
ジの初生（マイルドサージ）から、逆流を伴うディープサージに至る非定常流動メカニズムを解明し、圧縮機サージの予測技術を確立する。
令和 4 年度は、圧縮機サージの初生から逆流を伴うディープサージに至る非定常流動メカニズムを解明し、圧縮機サージの予測に利用可能であることを検証した。前年度のマイルドサージの計算に引き続き、逆流を伴うディープサージの計算を実施した。ディープサージの発生作動点やサージ発
生時の周波数および圧力変動レベルは実験結果と良く一致し、世界で初めてディープサージの直接シミュレーションに成功した。これにより、本シミュレーション技術によってディープサージの発生を予測できることを実証した。さらに、ディフューザ内で発生した旋回失速が突然、全失速に至ることなど、計算結果からディープサージ発生時の、圧縮機内部の非定常流動現象を詳細に明らかにした。
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 reaching from the surge inception (the mild surge) to the deep surge
including the reverse flows. This analysis is based on the DES (detached eddy simulation) andconsiders not only the compressor but also the entire system where the compressor is installed.
The surge simulations require a huge number of time steps to capture both of the short-cycle stall
phenomena in a compressor and the long-period surge phenomenon throughout the system, andsuch simulations will be achieved by using “Fugaku”.
In this fiscal year, we demonstrated that the simulation with “Fugaku” could predict the deep
surge, and clarified the unsteady flow mechanism from the surge inception (the mild surge) tothe deep surge accompanied by the reverse flow. Following the simulation of the mild surge donein the previous year, the simulation of the deep surge was conducted. The operating point of the
deep surge inception, and the frequency and amplitude of the pressure fluctuation at the deep
surge were in a good agreement with the experimental results. This is the first time in the world
the unsteady flow phenomena of the deep surge have been successfully simulated. We
demonstrated that the present simulation method could predict the occurrence of the deep surge.

- 20 -

Furthermore, we elucidated detailed unsteady flow phenomena in the compressor during thedeep surge(b) 本研究の内容
1) 計算モデル
車両用ターボ過給機に用いられる遷音速遠心圧縮機を解析対象とした。Fig. 2.2.3-1 に供試遠心圧縮機の概観を示す。本遠心圧縮機は、オープンタイプのインペラ、羽根なしディフューザ、およびスクロールから構成されている。Table 2.2.3-1 に主要諸元を示す。供試圧縮機インペラの翼枚数はフルブレードが 6 枚、スプリッタブレードが 6 枚の計 12 枚である。インペラとケーシングの間には0.28 mm の翼端隙間を有する。

Fig. 2.2.3-1 Test centrifugal compressor.
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]Tip clearance [mm]

52.5
0.28

Diffuser width [mm]3.42

Drive shaft speed [rpm]140,000

供試遠心圧縮機のサージ試験装置を Fig. 2.2.3-2 に示す。供試遠心圧縮機は、ターボ過給機のタービン側にスクリュ圧縮機から圧縮空気を供給することにより駆動される。タービン側に供給する圧縮空気の流量を調整することによって、圧縮機の回転数を最大 160,000 rpm までの範囲で任意に設定可能である。本研究では、回転数 140,000 rpm の条件で計算を行った。圧縮機流量は圧縮機の吐出側に配置したプレナムの下流に位置する流量調整弁により調整される。本試験装置では、供試圧縮機の吐出側配管系のプレナム容積を 3 とおりに変化させることができるが、本研究では最小容積の条件を対象とした。

- 21 -

Fig. 2.2.3-2 Schematic diagram of the experimental equipment for compressor surge.
上記のサージ試験装置を模擬した解析モデルを Fig. 2.2.3-3 に示す。サージの直接解析では圧縮機本体だけでなく、圧縮機が組み込まれたシステム全体を考慮する必要がある。したがって、本解析モデルでは、供試圧縮機に加え、その上流側および下流側の配管が考慮されている。ディープサー
ジ発生時には順流と逆流が周期的に繰り返されることから、逆流時に上流側配管の入口から吐き出された高温流体が順流時に再度吸い込まれる現象が生じる。この現象を再現するために、Fig.
2.2.3-3 に示すように、入口配管の上流には試験装置の外部空間を模擬した領域を設置している。外部領域は半球状の形状とし、その半径は圧縮機外径の 20 倍とした。また、サージ初生の流動現象を正確に再現するため、実際の試験と同様に、計算中に流量を漸減させる。流量制御は、実際の流量調節弁を模擬して行われる。流量調節弁の流路形状を計算格子によって直接再現し、計算途中に開度を変化させる方法は計算コストが非常に大きいため、流量調節弁については抵抗係数としてモデル化される。

- 22 -

20D

外部領域

圧縮機
入口配管

出口配管

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

- 23 -

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

180

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

Full blade
Splitter blade
320200

Rotation

116

116

Fig. 2.2.3-4 Computational grid in the impeller.
2) 並列計算の方法と効果（性能）
本研究で実施する圧縮機サージの直接計算では、サージ発生時における圧縮機内部の非定常流動まで詳細に再現する。そのため、圧縮機羽根車の失速現象を再現できる程度に、短い時間刻みを設定して計算を実施する。その結果、サージについてシステム全体にわたる長周期の流体振動現象を捉えるためには、膨大な時間ステップ数にわたって非定常三次元流動解析を実施することになる。
サージ 10 周期を計算する場合、7,500 step/rev×200 rev/cycle×10 cycle より、15,000,000 ステップの計算が必要であると見積もられる。上式において、step、rev、cycle はそれぞれ時間ステップ数、羽根車回転数、および、サージのサイクル数を表す。仮に、アプリケーションの時間ステップ
あたりの計算時間が 0.1 sec/step であれば、総計算時間は 17 日となり、現実的な時間内で計算を完了できることになる。そこで、アプリケーションであるターボ機械用 DES 圧縮性流れ解析についてチューニングをこれまで継続して実施してきている。
最適化は前年度まででほぼ限界に達していることが確認された。コードの性能は Table 2.2.3-3 に示す通りである。チューニング後の実行時間は 0.34 sec/step であり、目標とした 0.1 sec/step には至っていないが、概ね目標は達成されと判断した。したがって、以上でアプリケーションのチュー
ニングを完了とした。なお、アプリケーションのチューニングにあたっては一般財団法人高度情報科学研究機構（RIST）の高度化支援を利用した。

- 24 -

Table 2.2.3-3 Performance evaluation of tuning codes with the surge simulation.
実行時間[sec/step]

Asis との速度比

Asis

31.53
2.91

10.8

Tune 2 (2021/3)

0.53

59.5

Tune 3 (2022/1)

0.34

92.7

Tune 1 (2020/8)

3) マイルドサージの解析
初めに、逆流を伴わないサージ形態であるマイルドサージの計算を実施した。Fig. 2.2.3-5 は数値解析結果および実験結果から得られた供試圧縮機の性能特性である。横軸は質量流量、縦軸は全圧比を示しており、それぞれマイルドサージ発生時の質量流量と全圧比で無次元化されている。DES
解析結果については抵抗係数 Kt をそれぞれ 0.87、1.26、1.4、1.6、1.8 とした場合の結果である。
なお、質量流量が 0.6-1.4 の範囲における拡大図も図中に示した。拡大図に示した赤、黄、緑、青の実線はそれぞれの抵抗係数における計算過程で得られた全圧比および質量流量の変動を示しており、
橙色四角シンボルは実線で示した結果を算術平均したものである。また、ひし形シンボル付きの破線は実験結果であり、参考のため三角シンボル付きの実線で RANS 解析結果も併せて示している。
DES 解析から得られた性能特性は実験結果と比較して最大で 4 %程度高く評価されているが、性能特性の変化は RANS 解析と比較して定性的によく捉えられている。DES 解析結果において、抵抗係数を Kt=0.87 から徐々に増加させることにより流量が低下し、Kt=1.4 で圧力上昇が最大となっている。さらに抵抗係数を上げて流量を絞ると、性能特性は右肩上がりの分布を示しており、全圧比および質量流量の変動が大きくなっている。また、抵抗係数 Kt=1.8 において、逆流を伴わずに流量や圧力の長周期変動が発生するマイルドサージが確認される。質量流量は平均値から最大約 25 %だけ変動している。全体図中に示した灰色の領域は実験でマイルドサージの初生が計測された範囲である。繰り返し実験を実施する中で初生位置に変動が見られるものの、DES 解析でマイルドサー
ジが確認された流量はその範囲内であり、妥当な結果であるといえる。

- 25 -

Fig. 2.2.3-5 Experimental performance curve with numerical results obtained at some
operating points.
マイルドサージが発生した Kt=1.8 の DES 解析結果について実験結果との比較を行う。Fig. 2.2.36 に熱線流速計位置における質量流量の時間変化を示す。同図に灰色で示した領域は後述する流れ
場の可視化範囲である。実験結果および DES 解析結果は時刻 0.0 から 0.51 s における時間平均値よりそれぞれ無次元化している。DES 解析結果では過渡的な状態も含めると約 8 周期の流量変動が認められる。実験結果と比較すると、DES 解析で得られた質量流量の変動振幅は概ね一致しているが、サージング周期は DES 解析結果の方がわずかに長い。Fig. 2.2.3-7(a)に入口計測管での壁面圧力の時間変化を示す。実験結果は高応答圧力センサで計測されたものであり、DES 解析結果はセンサ設
置位置に相当する壁面の十数セルの圧力値を平均したものである。なお、図中縦軸に示した圧力は
時刻 0.0 から 0.51 s における平均圧力からの偏差である。同図より、入口配管における壁面で圧力が周期的に変動しており、Fig. 2.2.3-6 と同様、DES 解析よってマイルドサージが再現されていることが確認できる。圧力振幅を比較すると、実験結果と DES 解析結果は概ね一致している。また、
Fig. 2.2.3-7(b)は、Fig. 2.2.3-7(a)に示した壁面圧力の時間変化を周波数解析した結果である。実験
結果および DES 解析結果のいずれでもマイルドサージに相当する周波数にピークが存在する。実験結果のピーク周波数は 19.0 Hz であるのに対し、DES 解析結果のピーク周波数は 17.3 Hz とわず
かに低い値を示しているものの、圧縮機上流においてマイルドサージを定量的に捉えられている。
また、Fig. 2.2.3-8(a)および(b)に出口計測管での壁面圧力の時間変化および周波数解析結果を示す。
振幅を比較すると、入口計測管と同様、実験結果と解析結果は概ね一致している。周波数解析結果
において、実験結果では 18.0 Hz のピークを示しているのに対し、DES 解析結果は 17.3 Hz と入口
計測管と同様の結果が確認される。以上より、圧縮機下流においてもマイルドサージを定量的に捉えられている。

- 26 -

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

Fig. 2.2.3-7 Comparison of signals obtained atinlet-side in the experiment and in thesimulation: (a) signals (b) frequency
spectra.

Fig. 2.2.3-8 Comparison of signals obtained atoutlet-side in the experiment and in thesimulation: (a) signals (b) frequency
spectra.

次に、マイルドサージ発生時の非定常流動現象について示す。Fig. 2.2.3-9 は熱線流速計位置における DES 解析結果の質量流量の時間変化のうち、ある 1 周期について拡大したもの（Fig. 2.2.3-6中に図示）である。Fig. 2.2.3-9 中に示した代表的な 5 流量点（Point A-E）における瞬時の流れ場を比較する。なお、図中に示した各流量点の色は、後述する Fig. 2.2.3-12、Fig. 2.2.3-13 および Fig.
2.2.3-14 で示す各流量点の色と対応している。

- 27 -

Fig. 2.2.3-9 Closeup of time history of mass flow rate and positions of analysis points.
Point A-E の流量点について、瞬時の流れ場を周方向に平均した子午面流れ場を Fig. 2.2.3-10 に示す。いずれも流れ方向の格子線を基準に符号を定義した子午面速度 Cm で色付けされており、流線も併せて表示している。子午面速度は軸方向速度成分と半径方向速度成分を合成した速度、すなわ
ち子午面内の速度成分であり、その大きさは大気状態の音速で無次元化されている。いずれの流量点においてもナット着脱面に相当する位置で剥離が生じており、ディフューザ内ではシュラウド壁
面近傍、ハブ側壁面近傍で互い違いに剥離していることが確認される。一方、Fig. 2.2.3-9 中で示した流量変動に応じて、シュラウド壁面付近の逆流域の大きさに変動が見られる。Point A ではフル
ブレード上流のインペラ入口付近まで逆流領域が広がっているのに対し、流量の増加に伴い Point
B では翼端付近の逆流領域が縮小し、シュラウド壁面付近の流体が加速している。一方、流量が最小値を迎える Point D では逆流領域がインペラ入口上流にまで及び、大規模な再循環領域を形成している（Fig. 2.2.3-11）。また、下流側ではディフューザ内のシュラウド壁面近傍で逆流域が拡大している。その後、流量の回復に伴い、インペラ入口上流まで発達していた再循環領域が Point E では縮小しており、ハブ側の流体が加速している。以上より、圧縮機の流量変動に伴い、翼端付近で再循環領域が形成、発達、縮小、消失を繰り返していることが確認された。

Fig. 2.2.3-10 Meridional view of the distribution of streamlines at several time steps during
the mild surge cycle.

- 28 -

Fig. 2.2.3-11 Streamlines around impeller at Point D (minimum mass flow rate) colored bymeridional velocity.
Fig. 2.2.3-12(a)にインペラ入口（Position I）、Fig. 2.2.3-12(b)にインペラ出口（Position II）における子午面速度分布を示す。Position I および II については Fig. 2.2.3-10 における Point C の子午面に示してある。子午面速度は計算領域入口における音速で無次元化したものである。インペラ入口のスパン高さ 0 %（ハブ面）から 60 %の位置において、子午面速度は流量の変動に関わらず概ね一致している。しかし、スパン高さが 60 %から 100 %（シュラウド面）の位置においては、子午面速度が大幅に変動しており、流量が最も低い Point D および流量が回復途中の Point E では、スパン高さ 90 %から 100 %の位置において大規模な逆流を伴っている。一方、インペラ出口における子午面速度分布に注目すると、Point D および Point E のみシュラウド壁面近傍で逆流が見られるものの、その領域はインペラ入口の逆流領域と比べると小さい。また、流量の変動に関わらずシュラ
ウド壁面近傍を除くスパン高さ全体にわたって子午面速度分布は定性的に一致しており、流量に応じて速度分布が一様に変動している。また、Fig. 2.2.3-13 に各流量点の流れ場についてピッチ平均を施したフルブレードにおける 90 %スパン高さ位置の翼負荷分布を示す。縦軸に示したΔp は正圧面と負圧面の静圧差で定義し、計算領域入口における密度およびフルブレード先端周速に基づく動圧で無次元化している。コード方向 30 %から 100 %（フルブレード後縁）においては、圧縮機内の流量変動に関わらず翼負荷分布が概ね一致しているが、コード方向 0 %（フルブレード前縁）から30 %では著しい変動が認められる。特に、シュラウド壁面近傍で逆流を伴っていた Point D およびE では前縁付近では失速を伴っている。一方、Fig. 2.2.3-14 にスプリッタブレードにおける 90 %ス
パン高さ位置の翼負荷分布を示す。フルブレードの翼負荷と異なり、スプリッタブレードではコー
ド長全体にわたって翼負荷の変動が認められる。フルブレードと同様に、前縁付近の変動が最も顕著であり、Point D および E で失速を伴っている。

- 29 -

(a) Position I (near impeller inlet).

(b) Position II (near impeller outlet).

Fig. 2.2.3-12 Spanwise meridional velocity distribution.

Fig. 2.2.3-13
Chordwise blade loading
distributions of full blade at 90 % spanheight.

Fig. 2.2.3-14
Chordwise blade loading
distributions of splitter blade at 90 %span height.

翼端付近の流れ場を比較するため、Fig. 2.2.3-15(a)にスパン高さ 95 %における子午面速度分布を示す。なお、図中には逆流領域のみを表示しており、Point A と C、Point D と E は流れ場が類似しているため、ここでは Point B、C および D についてのみ示す。また、Fig. 2.2.3-15(b)には Fig.
2.2.3-15(a)で示した流量点における渦コア構造を示している。渦構造は無次元ヘリシティ Hn（渦度ベクトルと速度ベクトルの内積を両ベクトルの大きさで無次元化した値）で色付けされており、ス
パン高さ 95 %位置における曲面に沿ったエントロピー分布も併せて表示している。着目する渦を表示するため、Point D のエントロピー分布のみスパン高さ 90 %の位置を示しているが、スパン高さ
95 %と定性的な分布は同様である。まず、流量が最大となる Point B では、逆流領域が比較的小さ
く、全流路において翼端漏れ渦はスパイラル型の渦崩壊を伴っている様子が見られる。エントロピ
ー分布より、それぞれの渦崩壊は各流路内でのみ振れ回っており、隣接流路には干渉していない。
一方、Point C では逆流領域が Point B に比べて拡大し、FB3 における漏れ流れの一部はシュラウ
ド壁面近傍の逆流に巻き込まれて上流側へと移行しており、この渦を起点として一部の流路（図中
FB3-FB5）では隣接翼を乗り越えて spillage が生じている。流量が最小となる Point D では、逆流領域が著しく拡大し、全流路において隣接翼を乗り越えた逆流領域を形成しており、FB4 において

- 30 -

のみ翼端漏れ渦の渦崩壊が認められる。その他のフルブレード負圧面側には明確な翼端漏れ渦が確認されない一方で、竜巻状の剥離渦が発生している。Point D について、FB2 周辺の渦構造およびスパン高さ 90 %位置における静圧分布を Fig. 2.2.3-16 に示す。同図より、フルブレード先端の負圧面に足を持つ竜巻上の剥離渦が確認される。なお、Fig. 2.2.3-16 中に見られる翼（FB2）負圧面から上流に伸びる低圧域は、この剥離渦に対応したものである。このような竜巻状の剥離渦の発生により流路が閉塞され、Fig. 2.2.3-13 中の Point D における前縁付近の翼負荷はその他の流量点に比べて大幅に低下していると考えられる。以上のように、マイルドサージ発生中は翼端付近の流動変動が顕著である。特に圧縮機内の流量変動に伴い渦構造が常に変化しており、翼前縁付近で失速と回復を繰り返している。

Fig. 2.2.3-15 Instantaneous reversal-velocity distributions at 95 % span height (left) and vortex
structures and entropy distributions at 95 % span height (right).

- 31 -

Fig. 2.2.3-16 Vortex structure with streamline (top) and static pressure distribution (bottom)
at 90 % span height around FB2 (shown in Fig. 2.2.3-15) in Point D.
4) ディープサージの解析
次に、逆流を伴うディープサージの解析を実施した。Fig. 2.2.3-17 に DES 解析結果、RANS 解析
結果および実機試験結果から得られた性能特性を示す。横軸は質量流量、縦軸は全圧比を示しており、いずれも DES 解析においてディープサージが発生した Kt=7.5 の質量流量および全圧比で無次元化をしている。同図では、DES 解析結果は白抜きの赤丸、RANS 解析は白抜きの緑丸、実験結果
は白抜きの黒丸で示されている。なお、全圧比は入口計測管と出口計測管で計測された静圧と、熱線流速計位置における動圧を足し合わせることで簡易的に算出した。実験での質量流量に関しては、入口配管に設置された熱線流速計で取得した流速、大気状態の密度および熱線位置における配管断面積の積で定義した。DES 解析は抵抗係数 Kt を 0.87 から開始し、少しずつ Kt を増加させること
で圧縮機内に流入する空気の流量を減少させた。Fig. 2.2.3-17 にはディープサージが発生した
Kt=7.5 までの結果を記載した。ここで、それぞれの流量点における質量流量および全圧比は計算過程で得られるそれぞれの値を算術平均したものである。図中に水色で示す質量流量域は実験でディープサージの初生点が確認された領域を示している。DES 解析におけるディープサージの初生点はこの水色で示された範囲内に存在していることから、妥当な結果であると考えられる。

- 32 -

Fig. 2.2.3-17 Experimental and numerical performance curve.
Fig. 2.2.3-18 に抵抗係数 Kt=7.5 における熱線流速計位置での質量流量履歴を示す。時刻 0.0 から
0.18 s までは逆流を伴わずに流量が変動している。時刻 0.18 s 以降では急激な逆流が発生し、質量流量が周期的に変動していることから、ディープサージが発生していると考えられる。

Fig. 2.2.3-18 Time histories of the flow rate.
Fig. 2.2.3-19 には、ディープサージ発生時の DES 解析結果、および実験結果から得られたリサー
ジュ曲線を示している。いずれも DES 解析におけるディープサージ 2 周期分に相当する 0.25 秒分（Fig. 2.2.3-18 中の時刻 0.37 から 0.62 s）の結果を示している。また、実験結果における質量流量は熱線流速計による流量計測結果をもとにしているため、逆流領域を捉えることができないことを留意されたい。Fig. 2.2.3-19 より、実験結果および DES 解析結果のリサージュ曲線は定性的によく一致していることが確認される。特に、ディープサージの周期のうち高流量側から低流量側に向かう過程（図中①）では定量的にも良い一致を示しており、両者ともにそれぞれの性能特性の曲線に沿うような変動をすることも共通している。一方で、逆流から高流量側に流量が回復する過程（図中②）では、DES 解析から得られた全圧比が実験結果に比べて過少評価されていることが分かる。

- 33 -

Fig. 2.2.3-19 Experimental and numerical Lissajous curve of deep surge with performance
curve.
Fig. 2.2.3-20 および Fig. 2.2.3-21 は、入口計測管および出口計測管に設置した高応答圧力センサ
位置での壁面圧力変動履歴及びその周波数解析結果である。実験結果は高応答圧力センサで計測されたものであり、DES 解析結果はセンサ位置に相当する壁面の十数セルの圧力値を平均化したものである。いずれも 100 kHz でサンプリングを行った。ただし、抵抗係数の変更に伴う過渡的な影響を考慮し、DES 解析の周波数解析には、時刻 0.37 s 以降のディープサージ約 7 周期の結果のみを参照している。また、Fig. 2.2.3-20 及び Fig. 2.2.3-21 の縦軸は時間平均値からの圧力偏差である。

Fig. 2.2.3-20 Time histories and spectra ofpressure fluctuation at inlet duct.

- 34 -

Fig. 2.2.3-21 Time histories and spectra ofpressure fluctuation at outlet duct.

Fig. 2.2.3-20 より、入口計測管の壁面において、実験結果では約 9 周期の圧力変動が認められる一方で、DES 解析では約 7 周期の圧力変動のみ確認される。DES 解析では実験結果に比べて圧力を
最大約 4 kPa 過大に評価しており、周期も長く捉えられている。一方、圧力変動の波形は実験結果と定性的に一致しており、DES 解析結果からもディープサージに相当する周波数のピークが存在することから、圧縮機上流においてディープサージが再現されていることが示された。
Fig. 2.2.3-21 より、出口配管の壁面においても入口配管の壁面と同様に約 7 周期の圧力変動が認められる。DES 解析では実験結果に比べて圧力を最大 8.2 kPa 程度過大評価しており、周期も長く捉えられている。一方、圧縮機上流と同様に圧力変動の波形は定性的には一致しており、DES 解析結
果からもディープサージに相当する周波数のピークが存在することから、圧縮機下流においてディープサージが再現されていることが示された。また、周波数解析から、DES 解析で得られたディー
プサージの周波数は実験結果と比較して 2.2 Hz だけ低くなっていることがわかった。
Fig. 2.2.3-22 は Fig. 2.2.3-18 に示した熱線流速計位置における質量流量変動のうち、ある 1 周期を取り出したものである。次に、図中に示した代表的な 6 流量点（A~F）における瞬時の流れ場を比較する。解析点 A～F における 6 つの流れ場に対して、周方向平均処理を施した。Fig. 2.2.3-23に入口配管およびスクロール部分を含むインペラ周辺の子午面流れ場を示す。また、解析点 D および E についてはインペラ上流の入口配管にまで特徴的な流れ場の変動を伴っていることから、インペラ領域を含む入口配管の子午面流れ場を Fig. 2.2.3-24 に示している。参考のため、解析点 D と Eの中間点である流れ場（以下、解析点 DE と呼称）についても併せて示した。ここで、子午面速度は流れ方向の格子線を基準に符号を定義した。Fig. 2.2.3-23 より、圧縮機の流量変動に応じて流れ場が大きく変動していることが確認される。まず、解析点 A において、シュラウド壁面近傍では逆流が発生しており、その領域は入口配管にまで及んでいる。流量の増加とともに逆流領域は徐々に縮小し、流量が最大となる解析点 B においてはほとんど確認されない。その後、流量は低下し、解析点 C ではシュラウド壁面付近において、フルブレード先端側を起点とした逆流領域が発達し、流れがハブ側へと偏っている。流量が解析点 A を下回った解析点 D に至るころには、C で発生していたシュラウド側での逆流領域がスパン方向に拡大するだけでなく、大規模な再循環領域が入口配管にまで発達している（Fig. 2.2.3-24 の D）
。また、ディフューザ内においてハブ側に発生していた逆流領域がインペラ出口近くまで拡がっている。逆流発生直前の解析点 E では解析点 D で見られた再循環領域がさらに上流へと伸長している（Fig. 2.2.3-24 の E）。インペラ内部およびインペラ上流で逆流領域が発達した結果、インペラ出口からディフューザまでの子午面速度が著しく低下している。
逆流の発生する解析点 F では、完全に流路が閉ざされており、シュラウド壁面付近のみならず、インパラ領域やディフューザ領域の流れが完全に閉塞されている。Fig. 2.2.3-25 にインペラ入口断面を始点とする流線を示す。インペラ側へはほとんど流入することなく、すべての流体が入口配管側へと逆流する様子が確認される。

- 35 -

Fig. 2.2.3-22 Closeup of time history of mass flow rate (shown in Fig. 2.2.3-18) and positionsof analysis points.

Fig. 2.2.3-23 Meridional view of the distribution of streamlines around impeller region atseveral time steps during the deep surge cycle.

- 36 -

Fig. 2.2.3-24 Meridional view of the distribution of streamlines around suction pipe at several
time steps from analysis point D to E.

Fig. 2.2.3-25 Streamlines around impeller colored by meridional velocity at Ponit F.
ここで、ディープサージ発生時の流量が減少する過程における流れ場を調査するため、解析点 B、
C、D、E の 4 点に対して渦構造の抽出を行った。Fig. 2.2.3-26 には同解析点のフルブレードスパン
高さ 90 %位置における翼負荷分布を示す。なお、翼負荷は翼の正圧面と負圧面の静圧の差で定義しており、計算領域入口における密度およびフルブレード先端周速に基づく動圧で無次元化している。
加えて、同図中には、Fig. 2.2.3-24 で示した解析点 D および解析点 E の中間点 DE の結果も示している。また、Fig. 2.2.3-27 は無次元ヘリシティで色付けした渦構造である。インペラおよびディフ
ューザ領域のスパン高さ 90 %面の静圧分布も併せて示している。Fig. 2.2.3-27 より、解析点 B では
フルブレード前縁の翼端部から下流に向かって翼端漏れ渦が生じており、フルブレード負圧面に沿って低圧領域が広がっている。その結果、コード長全体にわたって翼負荷が高くなっていることが確認される。解析点 C ではフルブレード前縁部の翼端から巻き上がった漏れ渦がスパイラル型の渦崩壊を生じている。これは Fig. 2.2.3-23 で見られたシュラウド壁面近傍の逆流領域と対応しており、
漏れ渦の崩壊に伴うブロッケージ効果で流れがハブ側へ傾いた結果、コード長 20 %より下流側では解析点 B に比べて翼負荷が低くなっていることが考えられる。一方、解析点 D においては、2 ピッチにおいて Fig. 2.2.3-27 に示されるような翼負圧面に足を持つ竜巻上の剥離渦がフルブレード前縁- 37 -

付近に存在することが確認された。竜巻状の剥離渦は、翼端漏れ渦のスパイラル型渦崩壊によって放出された渦構造が翼負圧面の境界層とともに巻き上がったもので、大規模なブロッケージを形成する。その結果、流路が閉塞されるとともに急激に翼負荷が低下し、失速に至ったと考えられる。
同様に、解析点 E についても 3 つのフルブレードにて Fig. 2.2.3-27 に示すような負圧面に足を持つ
竜巻状の剥離渦が確認され、シュラウド面近傍での失速を伴っている。これらはシュラウド面近傍で生じる大規模な逆流によって巻き上げられた端面境界層から形成された剝離渦がシュラウド面に足を持ったものであると考える。

Fig. 2.2.3-26 Chordwise blade loading distributions of full blade.at 90 % span height.

- 38 -

Fig. 2.2.3-27 Meridional view of the distribution of streamlines around suction pipe at several
time steps from analysis point B to E.
また、Fig. 2.2.3-26 より、解析点 D は流量の減少割合が緩やかになる変曲点である。解析点 D から解析点 E までの区間では、フルブレード前縁付近で常に失速を生じており（Fig. 2.2.3-26）、翼先端の失速領域を起点とする再循環領域が概ね一定の速度で入口配管側へと発達している（Fig.
2.2.3-24）。逆流直前までこのような傾向が確認されていることから、翼先端における失速は再循
環領域の拡大に寄与しており、再循環領域の発達はシステム全体に及ぶ逆流の発生に密接な関わりがあると考えられる。

- 39 -

Fig. 2.2.3-27 より解析点 B から解析点 E へと流量が減少していく中で、ディフューザ領域の流体は徐々に圧力が上昇していることも確認されており、ディープサージの流量減少過程における特徴的な現象であると考える。
5) まとめと今後の課題
遷音速遠心圧縮機のマイルドおよびディープサージのシミュレーションに世界で初めて成功するとともに、サージ発生限界流量、および、サージ発生時の周波数および圧力変動レベルの予測が可能であることを実証した。さらに、サージ発生時の非定常流動現象を明らかにし、圧縮機の性能向上のための知見が得られた。今後は、サージ初生の流動現象を解析することにより、サージの初生
プロセスの解明に取り組む。
(4) リアルワールド自動車空力性能の予測(a) 研究成果の要旨
本実証研究テーマでは、電気自動車や自動走行車等の、次世代自動車の設計・開発に貢献するための HPC シミュレーション技術の構築と実証を行うことを目的として、実走行状態の自動車の空
力性能（空気抵抗、操安性、横風安全性等）の、車体形状再現性を数ミリまで高めた高解像度・長
時間シミュレーションによる評価を実現する。
令和 4 年度は、実走行状態の自動車の空力性能（空気抵抗、操安性、横風安全性等）をシミュレー
ションにより評価が可能であることを検証し、車体を固定した風洞実験条件に対して、接近流の変動や車体の揺動により空気抵抗が最大約 10 %変化することを明らかにした。また、レーンチェンジ
時のドライバー操舵を考慮したシミュレーションを実施し、特に接近流の変動が車体操安性に与える影響について検討を行った。これにより、設計開発現場におけるリアルワールド評価の重要性を実証すると共に、環境条件に応じた変動風の性状の再現といった、その利用指針を得た。
The purpose of this research is to build and validate the framework applying HPC technology
to contribute to the design and development of the next-generation vehicles such as electricvehicles or autonomous vehicles. To reach this goal, the evaluation based on the analysis of
aerodynamic characteristics of the vehicle in the actual running conditions (drag force, vehicle
maneuverability, crosswind safety ... et al.), high resolution for precisely representing the vehicle
shape in millimeter-scale and a long-term simulation is necessary.
In this fiscal year, we verified whether it is possible to evaluate the aerodynamic performance(aerodynamic drag, steering stability, crosswind safety, etc.) of automobiles in the real driving
conditions using simulations developed by us. As a result, it was clarified that the aerodynamic
drag changes by up to about 10 % due to the fluctuations in the approaching flow and the car
body’s 6 DoF (Degree of Freedom) motion, compared to the wind-tunnel-test conditions where theapproaching flow is uniform and the car body is fixed. In addition, we conducted a simulationthat takes into account the driver steering when the automobile changes the lanes, andinvestigated the effects of the fluctuations in the approaching flow on the drivability. As a result,- 40 -

the importance of the real-world evaluation at the design and development site was verified, anda guideline for its use was obtained.
(b) 本研究の内容
1) 解析対象
本研究の解析対象には、本田技術研究所より提供を受けた市販自動車「SHUTTLE」を再現したモデルを用いた。Fig. 2.2.4-1 のように、全長 4.440 m、全幅 1.695 m、全高 1.545 m となっている。また、Fig. 2.2.4-2 に示すとおり、床下形状やエンジンルーム内部も詳細に再現されており、実開発で用いられる形状再現性を有する CAD となっている。
3)リアルワールド燃費評価では、100 km/h で直進走行する運動を対象とした。その際に変動風の特徴の差異による影響を解析するため、Table 2.2.4-1 に示すように車幅方向（以下、y 方向）の積
分長さの異なる２つの変動風（Case 1、Case 2）を用意した。Fig. 2.2.4-3 に示す各 Case の瞬時ヨ
ー角場のように、Case 2 では、車体前方にあたる中央の翼列の格子解像度を上げることで、比較的長周期の乱流を生成している。
4-1)レーンチェンジ解析では、100 km/h で 3 m 右側にレーンチェンジする運動を解析した。その際に与えた変動風は、Table 2.2.4-2 に示すように、実際の走行試験で得られたデータと y 方向の積
分長さが近い値が得られた条件を用いている。

Fig. 2.2.4-1 Target model.

Fig. 2.2.4-2 Details of under-floor and engine bay of target model.

- 41 -

Table 2.2.4-1 Natural wind parameters of Case 1 vs. Case 2.

Case 1Case 2x

Turbulence Intensity [%]
y
z
7.24
8.17
6.52
2.90
4.54
1.31

Turbulence Length Scale [m]
x
y
z
0.18
0.10
0.07
0.18
3.36
0.04

Fig. 2.2.4-3 Snapshots of yaw angle Case 1(right) and Case 2 (left).
Table 2.2.4-2 Natural wind parameters of experimental vs. simulation.

Sim.
Exp.

Turbulence Intensity [%]
x
y
z
2.04
2.70
0.79
6.07
4.39
2.26

Turbulence Length Scale [m]
x
y
z
3.74
10.09
0.10
22.45
11.34
0.60

2) 変動風モジュールの開発
昨年度に引き続き、変動風作成モジュールの開発を行った。変動風の作成方法は Fig. 2.2.4-4 に示
すとおり、複数の翼型を車両上流にスパン方向に並べ、ピッチング運動を組み合わせて与えて作成する。翼のサイズや配置については、既存の変動風生成装置を参考にしている。それぞれの運動の周波数と振幅を変化させて、生成する乱れの積分長さ、スペクトル、乱流強度のコントロールを行う。

Fig. 2.2.4-4 Setups of the turbulence generator.

- 42 -

昨年度、与えた単一周波数の翼ピッチ運動に応じた変動が再現されたモデルに作成した変動風の周期条件を適用した結果を Fig. 2.2.4-5 に示す（図には、ピッチ角 6 度で周波数が 2 Hz と 8 Hz の
2 ケースについて示している）。周期条件が無い翼からの後流では、壁面効果により振れ幅（ヨー角）の減衰が見られ、周期条件を適用した変動が再現されているのが確認される（Fig. 2.2.4-5 左、横軸は翼端からの距離、縦軸はヨー角 β）。このようにして得られた変動風の性状を翼から 3 m 下流（車体の前縁を設置する場所）で計測した変動風のヨー角分布図（Fig. 2.2.4-5 右、横軸はヨー角 β、縦軸は当該ヨー角 β の発生頻度）に示すように低周波数領域の変動風がよく再現されている。

Fig. 2.2.4-5 Comparison of yaw angle time history and histogram at 3 m behind the airfoil w/
and w/o periodic boundary conditions. 2 Hz (Top) and 8 Hz (Bottom).
本年度は変動風のヨー角、乱流強度および積分長さを制御する観点から、Fig. 2.2.4-6 に示された提案のように、車両実験で計測された各方向の風速から得られた、車体に対するヨー角 β を増幅させたり、周波数帯を変更して最大周波数を制限した。この際、得られたヨー角 β のプロファイルは翼のピッチ角に適用され、この解法を通じて特定の成分だけを制御することが可能となり、実験の乱流特性を再現することに成功した。

- 43 -

Fig. 2.2.4-6 Schematic of the experimental equipment and reproducing yaw angles profiles(Left). Methodologies for controlling oscillatory flow: Amplification yaw angle (Right top)and Low-pass filter (Right bottom).
Fig. 2.2.4-7 に既存変動風データの、車体に対するヨー角を増幅させることによる変動風の性状分布を示す（図には、翼のピッチ角の増幅比 f が 1、1.25、1.5、2.0 の 4 ケースに翼から 1、3、5 m
下流について示している）。ヨー角 β と乱流強度は増幅比によって増加することが確認できるが、積分長は増幅比 1.25、1.5 に実走行試験で得られているのと同程度の数 m の積分長さの変動が作成され、それ以上の増幅比で減少する傾向が確認できる。

Fig. 2.2.4-7 Reproduced “yaw variance”, “turbulence intensity”, and “the overall integral lengthscale” as amplification ratio.
Fig. 2.2.4-8 に既存変動風データの、車体に対するヨー角に特定周波数以下の成分でフィルターリ
ングさせる（ローパスフィルター）ことによる揺れ幅を適用して変動風の性状分布を示す（図には、最大周波数が 2、10、20、30 の 4 ケースに翼から 1、3、5 m 下流について示している）。ここでは- 44 -

増幅比と同じようにヨー角 β と乱流強度もフィルター周波数によって増加することが確認できる。
一方、積分長はフィルター周波数 10、20 に実走行試験で得られているのと同程度の数 m の積分長さの変動が作成されるが、他の周波数でそれより非常に増加することが確認できる。

Fig. 2.2.4-8 Reproduced “yaw variance”, “turbulence intensity”, and “the overall integral lengthscale” as cut-off frequency of low-pass filter at a function of natural wind oscillation of the
airfoils.
二つのパラメータを適切に制御することにより、特定環境条件の変動風の生成が可能になった。
今後の課題として、環境条件に応じた変動風の性状を再現することに検討が必要である。
3) リアルワールド燃費評価
実走行時の燃費評価を目的に、上記の変動風モジュールと、車体運動やドライバー操舵を再現する自動車運動空力連成フレームワークを用い、直進走行のシミュレーションを行った。Fig. 2.2.4-9に各場合における車体まわりの瞬時流速を、Fig. 2.2.4-10 に各場合での空気抵抗係数（CD 値）の値を示す。なお、変動風生成のため翼を計算領域に配置した場合、翼のウェイクによって車体周囲に圧力勾配が発生するため、車体を配置しない翼のみの解析を別に行い、圧力勾配によって発生する力や車体への接近流の流速の変化を補正し、CD 値を算出した。Fig. 2.2.4-10 に示す CD 値の比較か
ら、変動風下の方が増加している。さらに、積分長さが長い場合では増加幅がより大きいことがわかる。これらの要因として、Fig. 2.2.4-11 に示す車体のヨー角変動の変化が挙げられる。ヨー角変動においても CD 値と同様な増加傾向が見られることから、変動風によって車体運動が大きくなり、CD 値が増加するという現象を捉えられたと考えられる。また、積分長さの大きい長周期の乱れによってより車体運動が大きくなり、それに伴って本解析では最大で約 10 %CD 値が大きくなることが確認できた。

- 45 -

Fig. 2.2.4-9 Snapshots of velocity profile w/o incoming turbulence (left) and Case 1(center) andCase 2 (right).

Fig. 2.2.4-10 CD in each case.

Fig. 2.2.4-11 Yaw angle time history.
4) リアルワールド操安性評価4-1) レーンチェンジ解析
ここでは、実走行時の操縦安定性を評価するために、レーンチェンジ解析を行った。Fig. 2.2.4-12（縦軸はステアリングの最大舵角に対する割合）にドライバー操舵の時間変化を示す。図のように、変動風によってドライバー反応に差が生じたことが見て取れ、特に変化が大きい 2 s から 2.5 s の間では、最大の差が見られる時刻では、ステアリング入力量が 0.1 %から 0.3 %となり 3 倍の入力量の差が捉えられた。その要因として、ヨーイングメーメントが挙げられる。Fig. 2.2.4-13 に 2 s から
2.5 s のヨーイングモーメントを示す。図のように、変動風によって変化が生じており、それによりFig. 2.2.4-14 に示すようにヨー角に影響を与えている。以上のように、変動風の影響によって現れ- 46 -

たヨーイングモーメントの変化から生じたヨー角変動を修正するために、最大で 3 倍のドライバー
操作の差が現れることが確認できた。

Fig. 2.2.4-12 Steering input time history.

Fig. 2.2.4-13 Yawing moment time history.

Fig. 2.2.4-14 Yaw angle time history.
4-2) 追い越し解析
これまでのリアルワールド操安性評価では、車両 1 台での解析を行ってきた。しかし、実際の走行では他車からの空力的影響も無視できない。特に追い越し時には、追い越す車両と追い越される車- 47 -

両の 2 台の相互作用によって操縦安定性が低下することが知られている。そこで、4-1)で用いられるフレームワークに、追い越し解析が可能となるようアップデートを行った。ここでは、追い越し
解析に先立ち行った追い抜き解析での結果を報告する。Fig. 2.2.4-15 のように、高速道路上での追い抜きを想定し、80 km/h で走行中の車両を、100 km/h で走行中の車両が右側車線から 4 秒間で追い抜くという条件で解析を行った。なお、今回の解析では追い抜かれる車両に着目し、この車両でのみ、車両の多体運動解析が可能なオープンソースソフトウェア Chrono との連成解析を行った。
Fig. 2.2.4-16 に、追い抜かれた車両がとった経路と、ドライバー応答についての図を示す。約1.5~2 s の間で 2 台の車両は最も接近し、この間に追い抜かれる車両は不安定な走行経路をたどった。また、ドライバーは経路修正のため、経路とは逆向きのステアリング入力を行ったことが確認できた。Fig. 2.2.4-17 には、追い抜き前（0.6 s）と追い抜き中（2.16 s）の瞬時流れ場を示す。この図より、追い抜き前にはそれぞれの車両で独立していた後流が、追い抜き中には複雑に混合していることが確認できた。これらより、追い抜きによる影響が再現されており、追い抜き解析が可能となったことが確認できた。また、条件の変更によって追い越し解析も可能である。
さまざまな条件下で、この追い越し解析を活用し、他車の影響を考慮した操安性評価を行うこと
で、これまでは確認できなかった現象の発見や、そのメカニズム解明が可能になると期待される。

Fig. 2.2.4-15 Schematic of overtaking analysis.

Fig. 2.2.4-16 Path of overtaken vehicle (left) and steering input (right).

- 48 -

Fig. 2.2.4-17 Snapshots of the velocity field at 0.6 sec (left) and 2.16 sec (right).
(5) リアルワールド自動車空力音予測(a) 研究成果の要旨
本実証研究テーマでは、電気自動車や自動走行車等の、次世代自動車の設計・開発に貢献するための HPC シミュレーション技術の構築と実証を行うことを目的として、ボンネット隙間やフロン
トグリルから発生する狭帯域音、ピラーやドアミラーから発生する広帯域音、空力・構造振動・音響連成による車内騒音を対象に、実走行時の空力音予測を行う。
令和 4 年度は、ボンネット隙間やフロントグリルから発生する狭帯域音、ピラーやドアミラーから発生する広帯域音、空力・構造振動・音響連成による車内騒音を対象に、実走行時の空力音予測が可能であることを検証し、以下の知見を得た。狭帯域音の解析について、格子解像度を 1 mm 以下にすることによりボンネット隙間からの狭帯域音の予測が可能である見通しが得られ、狭帯域音の発生機構が明らかなとなった。また、広帯域音の解析について、格子ボルツマン法（LBM）を用いて車両まわりの広帯域空力音を従来の分離解法と同程度の精度で予測できることを確認した。空力・構造振動・連成解析により車内騒音を定量的に評価し、流体加振、音響加振の寄与を明らかにした。
The purpose of this research theme is to establish and demonstrate HPC simulation technology
to contribute to the design and development of the next-generation automobiles including electricvehicles. For this purpose, the narrow-band noise generated from the hood gap and the front
grille, the broadband noise generated from the pillars and the door mirrors will be simulated.
Moreover, the interior noise caused by the coupling of aerodynamics, structural vibration, andacoustics will be predicted.
In this fiscal year, we verified that it is possible to predict aerodynamic noise in the actual
driving conditions: the narrow-band noise generated from the hood gap and the front grille, thebroad-band noise generated from the pillars and the door mirrors, and the interior noise caused
by aerodynamic, structural vibration, and acoustic coupling. The following findings wereobtained: The analysis of the narrow-band noise showed that the narrow-band noise from thehood gap can be predicted by reducing the grid resolution to 1 mm or less, and the mechanism of
the narrow-band noise generation was clarified. In addition, it was confirmed that the lattice
Boltzmann method (LBM) can predict the broadband aerodynamic noise around the vehicle withthe same level of accuracy as the conventional decoupled solution method for analysis of the- 49 -

broadband sound. Aerodynamic, structural vibration, and coupled analysis were used toquantitatively evaluate the interior noise of the vehicle, and the contributions of fluid excitation
and acoustic excitation were clarified.
(b) 本研究の内容
1) 狭帯域音予測
自動車の狭帯域音を予測するため、昨年度に引き続き Fig. 2.2.5-1 に示す車両まわりの音の直接解析を行った。流れの中に微小な段差・隙間が置かれると、発生した空力音が渦自体に影響を及ぼす
フィードバック音が生じることが知られているが、一般にこのような空力的な自励音は上流の流れが層流である場合に起きることや、発生した渦が後方の障害物のエッジに衝突するような場合に発生することがわかっている。しかし、自動車のように形状が複雑な場合、上流の流れ（境界層）の状態を確認することが難しく、また、本モデルのようにボンネット内のエンジン等の圧力損失により主流流れが変化し、はく離した流れの衝突箇所が予測できない場合がある。自動車から放出される狭帯域音を予測する場合、音が発生する部位だけでなく、車両全体の流れ場を考える必要があるのはこのためである。本研究では、この問題を明らかにするため、狭帯域音が発生する箇所だけでなく、エンジンルームを含む車両全体の解析を行った。エンジンルームを含む複雑形状を解析するために直交格子による解析アプリケーションである CUBE を使用した。また、空力音を解析するた
めに圧縮性ナビエ・ストークス（NS）方程式を解き、音を直接求めることを試みた。音の発生個所であるボンネット付近は小さなキューブを用いて計算格子の解像度を上げ、車両まわりの解像度
は車体まわりの比較的大規模な流れが再現できるように調整した。

- 50 -

Fig. 2.2.5-1 Computational grid for narrow-band sound analysis. Each box has 163 cells.
計算条件としては、実車においてボンネットから狭帯域音が発生した 54 km/h に車両速度を設定した。音の発生源と考えられるボンネット付近の解像度を高くし、最小格子サイズは 0.2 mm（格子数 20 億）である。昨年度までの解析においては、CUBE の標準圧縮性スキームである ROE スキームを用いて解析を行っていた。しかし、0.2 mm オーダの格子サイズでは、車体表面から剥離したせん断層がエンジンフードと強く干渉することなくエンジンルーム内に侵入し、それに伴う音波の発生メカニズムが実験とは異なっていた。これは ROE スキームの強い散逸性に起因すると考えられる。ボンネット付近の最小格子サイズを局所的に 0.05 mm とした解析では実験と対応した流れが得られたが、これを車体横方向全域に適用すると数兆オーダの格子数となり、「富岳」を用いてもこのような解析は難しい。そこで、0.2 mm オーダの格子サイズで実験と対応した流れを得るた
めに、本年度は散逸の小さい KEEP（Kinetic Energy and Entropy Preserving）スキーム [1] を実装し、狭帯域音解析への適用を行った。KEEP スキームは、各種フラックスが一次保存量だけでなく、非圧縮、非粘性の極限で保存量である運動エネルギーとエントロピーの保存も満たすように離散化されている。そのため、数値粘性なしで長時間安定計算することが可能である。ただし、KEEP スキームは二次精度の中心差分で離散化されているため、位相誤差による数値振動が問題となる。これを抑えるため、淺田らによって提案された数値散逸フィルター [2] を同時に実装した。
このフィルターは、淺田らの論文内でも KEEP スキームと併用されているものである。フィルター
は、NS 方程式で解く全ての一次保存量𝑞𝑞に対して以下のようなフィルタフラックス𝑓𝑓を用いて適用される。ここで、𝑏𝑏0 = −3/16、𝑏𝑏1 = 1/16であり、フィルター係数𝛼𝛼(= 0~1)の大きさによって散逸強度を調整することができる。

- 51 -

𝑞𝑞�𝑖𝑖 = 𝑞𝑞𝑖𝑖 − 𝛼𝛼�𝑓𝑓𝑖𝑖+1/2 − 𝑓𝑓𝑖𝑖−1/2 �
𝑓𝑓𝑖𝑖+1/2 = [𝑏𝑏0 (𝑞𝑞𝑖𝑖+1 − 𝑞𝑞𝑖𝑖 ) + 𝑏𝑏1 (𝑞𝑞𝑖𝑖+2 − 𝑞𝑞𝑖𝑖−1 )]非粘性 Taylor-Green 渦問題に対して、CUBE へ実装した KEEP スキームと、昨年度までの解析で用いていた低マッハ数用 ROE スキーム（ROE_LM）の比較を行った。その結果を Fig. 2.2.5-2 に示す。純粋な KEEP スキーム（𝛼𝛼 = 0）では完全に運動エネルギーとエントロピーが保存されていることが確認できる。一方、ROE_LM スキームでは、運動エネルギーが減衰し、エントロピーが増大する。KEEP スキームにフィルターを適用し、𝛼𝛼を大きくしていくと運動エネルギーの減衰率は大きくなっていくが、ROE スキームより散逸が強い（𝛼𝛼 ≳ 0.01）条件においても、エントロピーの保存性はよいことが分かった。

Fig. 2.2.5-2 Kinetic energy and entropy evolution in Taylor-Green-Vortex problem. M0 means
mach number (0.05), ρk, ρs means kinetic energy and entropy. The index 0 means initialvalue.
狭帯域音解析においては、クーラン数を同等とした非粘性 Taylor-Green 渦問題において、
ROE_LM スキームより散逸強度が低いフィルター係数（𝛼𝛼 = 0.2）で解析を行った。Fig. 2.2.5-1 の0.2 mm 規模格子における解析結果を Fig. 2.2.5-3 に示す。ROE_LM スキームでは、車体表面から剥離した流れはエンジンフードと強く干渉することなくエンジンルーム内に侵入し、エンジンルー
ム内の構造物に衝突して音が発生しているが、KEEP スキームではせん断層の厚みが増し、エンジ
ンフードの先端から音が発生している様子が確認できる。このように、0.2 mm 規模の格子でも実験と対応した流れ場で音の解析が可能となった。
その一方で、KEEP スキームではエンジンルーム内に高速の流れが発生している。これは圧縮性の IBM（Immersed Boundary Method）によって発生した漏れが、数値散逸が弱くなったことにより顕在化したものだと考えられる。このような漏れは車体まわりでも見られ、音圧レベルにも影響を与えてしまう。これを抑えるための IBM の改良および壁モデルの導入が必要である。

- 52 -

Fig. 2.2.5-3 The velocity magnitude (top) and sound pressure (bottom) near the bonnet. The
ROE_LM (left) and KEEP + filter (right).
2) 広帯域音予測2-1) CUBE を用いた直接空力音解析
自動車の広帯域音を予測するために、圧縮性流体ソルバーCUBE を使用し、A ピラーやドアミラーから発生する風切り音の直接解析を実施した。渦と音を直接解析で捉えるための格子解像度指標を確立し、自動車風切り音の予測に重要な A ピラー縦渦、サイドミラー後流を捉えることができた。
現状の格子（Fig. 2.2.5-4）では、Fig. 2.2.5-5 に示したように高周波数側の音圧スペクトルは 2.5
kHz まで捉えられることを確認した。

Fig. 2.2.5-4 The cube voxel distributions (left) and the volume mesh (right) around the A-pillar
and side mirror for the aeroacoustics noise simulations.

- 53 -

Fig. 2.2.5-5 Vortical structures of an instantaneous flow field visualized by iso-surfaces of the
Q-criterion colored with velocity magnitude around the A-pillar and side mirror (left) andsound pressure level (SPL) at a point location on driver's side-window surface (right).
多項式近似で速度変動を与える変動風生成モジュールを CUBE に組み込み、自動車実走行時の風速・風向変動と乱流特性（乱流強度や渦スケール）の再現に取り組んだ。さらに、変動風を模擬した自動車の風切り音解析を実施し、Fig. 2.2.5-6 に示すサイドガラス表面で測定した圧力の時間周波数特性（音圧スペクトログラム）から、バサバサ音の音響特性を調査した。その結果、バサバサ音
は進行方向の風速変動と相関があり、A ピラー縦渦の強弱に依存することがわかった。
しかし、シミュレーションでは自動車メーカーが問題視する 1 kHz 以上の高周波側の変調を捉えていないので、今後は音圧レベルだけなく、変調による音の違いを見るための解析を実施する予定である。また、Fig. 2.2.5-7 に、このような変調を伴う風切り音の発生メカニズムを解明するため、
A ピラー／サイドミラー周辺の剥離流れの計算結果を示す。

Fig. 2.2.5-6 The time history of the sound pressure (top) and the spectrogram (bottom).

- 54 -

Fig. 2.2.5-7 Flow streamlines around the A-pillar and side mirror.
2-2) FFX を用いた直接空力音解析
格子ボルツマン法を用いて車両まわりの空力音を解析するシステムを構築し、車両モデルを用いた検証を行った。
2-2-1) FFX の開発
広帯域音解析のため、本課題では、Lattice Boltzmann Method（以下、LBM）に基づく流体解析
コード FFX を開発した。令和 3 年度までに主要計算部分を含むプロトタイプの開発は完了しており、令和 4 年度は入力データ作成機能および計算結果出力機能を拡張し、FFX を流体解析システム
として完成させた。
a) Refine 機能の実装
マルチ要素作成プログラム mltcub は物体近傍キューブを Refine することでマルチレベルキュー
ブを作成する機能をサポートしていたが、これに加え、ユーザが指定したバウンディングボックス
内のキューブを Refine する機能を追加した。マルチレベルキューブの作成方法として以下 3 とおりの方法でマルチレベルキューブを作成できるよう機能拡張した。
 Bottom-up：物体形状に近い領域を Refine する。
（従来方法）
 Bounding Box（BB）：BB により、Refine 領域を指定する。
 BB + Bottom-up：はじめに BB で Refine し、その後 Bottom-up で Refine する。接続情報の作成アルゴリズムの都合で順番を逆にすることはできない。
機能拡張した周辺プログラム mltcub の動作確認として、キューブ数が 30×20×4 であり、XY 断面の中央に角柱を設置したベースのシングルレベルキューブに対し、以下の 3 とおりのマルチレベ
ルデータを作成した。拡張した機能が正常動作することを確認した。
 マルチキューブ 1：Bottom-up、3 レベル- 55 -

 マルチキューブ 2：BB、3 レベル
 マルチキューブ 3：マルチレベル CUBE を初期値に Bottom-up、2 レベル（総レベル 4）作成されたキューブのキューブ数を Table 2.2.5-1 に、キューブ配置を Fig. 2.2.5-8 に示す。
Table 2.2.5-1 Number of multi-level cubes for computing flow around square cylinder.

Lv. 1
Lv. 2
Lv. 3
Lv. 4
合計

マルチレベル
CUBE1
2,048
512
2,304
4,864

マルチレベル
CUBE2
6,144
2,432
2,000
10,576(a) マルチレベル CUBE1

マルチレベル
CUBE3
7,680
5,184
2,432
2,000
17,296(b) マルチレベル CUBE2(c) マルチレベル CUBE3
Fig. 2.2.5-8 Multi-level cubes for computing flow around square cylinder.

- 56 -

b) 入力データの作成機能の拡張
LBM は直交格子に基づいているため格子の幾何情報をメモリに保持する必要がなく、グリッドあ
たりのデータ数が少なくメモリ容量が小さく大規模の計算が容易である利点がある。粒子モデルとして D3Q27 を用いた場合でもグリッドあたりの word 数は 50 程度（倍精度であれば、0.4 kB/grid)である。「富岳」の計算ノードあたりの有効メモリ容量を 28 GB（32 GB より、MPI 等のシステムが使用する分として 4 GB を差し引いた値）とすると、ノードあたり 7,000 万グリッドを乗せることができ、1 万ノードを用いれば 7,000 億グリッドの計算が可能である。大規模解析が容易に実行できることより、後処理のユーザ負担が大きくなる。たとえば、1,000 億グリッド程度の計算は「富岳」数千ノードを用いて比較的容易に計算することが可能であるが、FFX が出力する流れ場ファイルは 1.6 TB（100 G-grid×4 変数×4、可視化用流れ場は単精度で出力）になり、このデータの転送・変換・可視化をする負担が大きい。FFX による計算結果の可視化の作業負担を軽減するため、
流れソルバーに流れ場の可視化用データの出力機能として、ブロック出力機能および表面データ出力機能を実装した。Fig. 2.2.5-9 に FFX システムのデータフローを示す。上記の 2 機能に実装による、物体表面の圧力データおよびユーザが指定する任意のブロックデータにおける流れ場を出力することが可能となり、ユーザの可視化における作業負担が大幅に軽減された。Fig. 2.2.5-10、Fig.
2.2.5-11 にそれぞれ表面データ出力機能、ブロックデータ出力機能を用いて可視化した車体表面まわり流れの可視化例を示す。

Fig. 2.2.5-9 Data flow in FFX-system.

- 57 -

Fig. 2.2.5-10 Distributions of density on vehicle surface computed by FFX.

Fig. 2.2.5-11 Distribution of streamwise momentum on center plain computed by FFX.
c) 基礎検証：角柱まわりの流れ解析
前年度に検証解析を行った角柱まわりの流れ解析を引き続き実施し、格子ボルツマン法による空力音解析における課題について検討を行った。格子ボルツマン法の標準モデルある BGK 方程式から導かれ単一緩和時間モデル（Single Relaxation Time、SRT）は、速度場の発散項にかかる物理定数である体積粘性率が NS 方程式の係数よりも 3 倍大きいにも関わらず、発散項の不安定性から発散しやすいことがわかった。この問題を解決するため、格子ボルツマン法では多緩和時間モデル（Multiple Relaxation Time、MRT）が使用されるようになってきている。MRT は発散項（もしくは音響場）の物理定数である体積粘性率を実際よりも大きくすることにより計算の安定化を図っている。この効果により非圧縮性の流れ場は安定した解が得られるが、音場については減衰が大きくなるという弊害がある。Fig. 2.2.5-12 に角柱まわの流れ解析結果（レイノルズ数 Re=4×104）を示す。流れ場については SRT、MRT ともにほぼ同じ結果が得られるが、Fig. 2.2.5-13(a)、(b)に示すように音響場については結果が異なる。SRT 解析では音場に非物理的な振動が発生するが、MRT では非物理的な振動が生じない。ただし、音場全体に減衰の効果が作用するため、高周波の乱流成分による音は過小評価となる。

- 58 -

Fig. 2.2.5-12 Velocity field around a square cylinder by LBM analysis.

(a) SRT(b) MRT

Fig. 2.2.5-13 Sound field around a square cylinder by LBM analysis.
Fig. 2.2.5-14 に格子ボルツマン法による音波の減衰計算を行った結果を示す。赤い実線は格子ボル
ツマン方程式から求めた減衰を考慮した音波の方程式の理論解と解析結果を比較した結果である。
縦軸は音波（密度振幅）の減衰割合、横軸は音が放射されてからの時間を周期で無次元化した値である。音源からの距離と見做すこともできることから、音源からの距離と考えることができる。格子ボルツマン法により理論解と一致する結果を得ることができるが、解析解自体が一般的な NS 方
程式の解析解とは体積粘性率の違いの分だけ異なる。式=p exp ( −α x )γ −1  2
4π 2 
=
α
µ
+
χ
+
κ f
(
)

ρ C p 2 ρ c 3 

(2.2.5-1)は音波の減衰式である。ここでｘは音源からの距離、α は減衰係数、c は音速、ρ は平均密度、γ は比熱比、κ は熱拡散率、Cp は定圧比熱、f は周波数である。粘性係数 μ と体積粘性率 χ の和により減衰の大きさが決まるが、NS 方程式では粘性係数と体積粘性率の和は粘性係数の 4/3 であるのに対- 59 -

して、SRT は粘性係数の 2 倍、MRT では粘性係数の 10 倍程度の値となり、音波の減衰が非常に大きくなる。LBM の計算安定性を考慮すると体積粘性率を大きくし、音場（エネルギー）に関する緩和時間を長くする必要があるが、前述したように音場自体は高周波数の変動が抑えられ、なだらかになる点に注意が必要である。

Fig. 2.2.5-14 Attenuation characteristics of sound waves by LBM analysis.
LBM で高レイノルズ数流れを安定に計算するには、多緩和時間モデルを用いる必要があることから、体積粘性率を大きくした場合の音の減衰が解析結果に及ぼす影響の評価を行なった。
Fig. 2.2.5-14 は自動車の解析において一般的な計測条件に近い、車体から 5 m 離れた点における減衰の影響を評価した結果である。NS 方程式により得られる音場に対して、MRT の解析結果の誤差を周波数別に求めた（Fig. 2.2.5-15）。凡例は音の周波数[Hz]を表す。その結果、5 m 離れた位置の音の減衰は 2 kHz では 0.15 %以下、10 kHz で 3.6 %、20 kHz で 13.8 %であり、音圧に換算するとそれぞれ 0.01 dB、0.3 dB、1.3 dB である。この程度の音圧レベルの誤差であれば、実用上はあまり問題がないが減衰が生じるということは音場の位相も変化するため、フィードバック音のように音の位相差が問題となる場合は注意が必要である。実用上はほとんど問題ないが、格子ボルツマ
ン法では、体積粘性率の大きさを大きくすることにより解析の数値安定性を得ている点についてユーザが認識する必要がある。NS 方程式においても数値粘性の影響があることを考えると、LBM では体積粘性率の影響として数値的に影響が既知である点はむしろ良いことなのかもしれない。

- 60 -

Fig. 2.2.5-15 Evaluation of excess attenuation in LBM analysis for NS equation.
LBM 解析を用いて音場を解析する場合、境界における音の反射を抑える必要がある。本研究では平衡分布関数を境界付近で指定値（無反射条件）に漸近させる形で音の反射を抑制した。式2+ ∆t )fi (x + ci ∆t , t =x 1fi (x, t ) −  fi (x, t ) − fi eq (x, t )  − α  d   f i eq (x, t ) − f i t arg et (x, t ) τ
 Ld 

(2.2.5-2)の右辺第 2 項に反射を抑制した場合の平衡分布関数（ftargert）を定義することで境界付近での反射を抑制する。ここで feq は局所平衡分布関数、τ は緩和時間係数，α が減衰率，xd は吸収境界からの距
離，Ld は吸収境界の幅である。Fig. 2.2.5-16 にこの無反射境界領域を適用した場合と適用しない場
合の角柱まわりの音圧スペクトルを示す。無反射境界条件を設定しない場合は、解析初期に発生した強い圧力変動が解析領域内で反射を繰りかえすため、無次元周波数 10 以上の高い周波数の変動が支配的になる。無反射境界領域を適用することにより音場を解析することが可能となった。Fig.
2.2.5-17 に角柱まわりの音場の解析結果を示す。

- 61 -

Fig. 2.2.5-16 Comparison of sound pressure spectra with and without non-reflective boundary
region.

Fig. 2.2.5-17 Sound field around a square cylinder.
d) 基礎検証：フロントグリル模型の音響共鳴音
Fig. 2.2.5-18 に自動車のフロントグリルを模擬した風洞実験の模式図と発生音の予測結果を示す。
平行平板間に生じる気柱共鳴音と流れ場から発生する空力音が連成することにより強い空力音が発生することが知られている。解析結果には強いピーク性の音が生じていることがわかる。空力音
は音源の空間的な構造に依存するため、スパン方向の空間相関を解析領域内で考慮した場合と考慮
しない場合を仮定し、空間相関の補正を行った。Fig. 2.2.5-19 にフロントグリルまわりの音場を示す。実験結果と解析結果はピーク周波数が一致すること、相関長補正を行った結果は実験値と一致- 62 -

することを確認した（Fig. 2.2.5-20）。このことから、LBM 解析により音場の共鳴と流れ場の連成問題が精度良く解析できることを確認することができた。

Fig. 2.2.5-18 Schematic diagram of wind-tunnel experiment simulating front grille (left) andpredicted sound pressure spectrum (right).

Fig. 2.2.5-19 Aerodynamic resonance sound field around the front grille model.

- 63 -

Fig. 2.2.5-20 Sound pressure spectrum of the front grille model.
e) 基礎検証：マルチキューブ解析
直交格子系の数値解析では格子点数が膨大になる問題があるため、格子サイズを階層的に変えることが一般的である。前述のように、FFX では物体近傍の格子を階層的に変える Refine と領域内の格子サイズを指定する Bounding Box（BB）の二つのマルチキューブ機能が実装された。Fig.
2.2.5-21 は BB によるマルチキューブ解析の事例である。マルチキューブ解析を適用することにより計算格子数を 1/10 から 1/100 程度にできることを確認した。

Fig. 2.2.5-21 An example of multicube analysis using bounding boxes.

- 64 -

2-2-2) 車両解析
マルチキューブ解析を行うことにより自動車の解析規模をより小さくすることが可能となる。本
研究では、マルチレベルを 3～5 に設定することにより車両まわりの流れと音を約 1,500 億格子で解析できることを確認した。プリ処理にかかる時間を Table 2.2.5-2 に示す。短時間で STL データから解析データを生成できることを確認した。実車のエンジン等を含む詳細データでも 5 時間以内でプリ処理を終えることができる。CFD 用であれば 2 時間以内で処理が可能である。
Table 2.2.5-2 Pre-processing time for LBM analysis.
Pre-processing timeMulti-Cube
Generation

Grid
points[billion]

Grid Length/Characteristic
Length[mm]

LEVEL

3.3

0.31/201

4 sec

17 sec21 sec21.3

0.04/ 20

4

1 min 31 sec

1 hour 46
min 25 sec

1 hour 47
min 56 sec

0.3

0.47/30

42 min 40 sec

8 min 26 sec

11 min 6 sec

4.2
143.2
3.5

1.56/1600
0.12/1600
2.5/50001
3
5

4 min 59 sec
9 min 5 sec
3 min 52 sec

Car (CFD)

155.3

0.78/5000

3

6 min 50 sec

Car (fulldetails)

151.3

0.78/5000

3

9 min 44 sec

9 min 20 sec
16 min 59 sec
23 min 40 sec
1 hour 18
min 13 sec
4 hour 38
min 18 sec

14 min 19 sec
26 min 4 sec
27 min 32 sec
1 hour 25
min 3 sec
4 hour48
min.2 sec

CASE
Cylinder(Re 150)
Cylinder(Re 4×104)
Front
grilleBicycleBicyclesCar (CFD)

STL

Total

Fig. 2.2.5-22 にマルチキューブ解析の実施例を示す。車両まわりに複数のサイズのキューブを配置
し、LBM 解析を実施する。音場の解析結果を Fig. 2.2.5-23 に示す。ホイールハウス、ドアミラーの近くに強い音場が存在することがわかる。Fig. 2.2.5-24 に音圧スペクトルの解析結果を示す。風
洞実験結果と 1 kHz 付近までスペクトルが良く一致することを確認した。今回の解析は最小格子サ
イズが 0.8 mm 程度であるため、高周波数まで解析するにはメッシュサイズが不十分であったが、LBM を用いて車両まわりの空力音を定量的に評価できることがわかった。Fig. 2.2.5-25 に音の減衰特性を示す。解析空間は地面を除き、無反射境界となっているため、音の減衰はほぼ自由音場に近いものとなっていることを確認した。

- 65 -

Fig. 2.2.5-22 Cube arrangement for LBM analysis (Multicube: Level 5).

- 66 -

(a) TOP VIEW(b) SIDE VIEW
Fig. 2.2.5-23 Sound field distribution around the vehicle.

- 67 -

Fig. 2.2.5-24 Comparison of experimental and LBM sound pressure spectra.

Fig. 2.2.5-25 Attenuation characteristics of the sound field around the vehicle.
車内騒音への影響を調べるため、車体表面（屋根部）とドアミラーわきのサイドシールドの圧力場に対して波数・周波数スペクトル解析を行った。Fig. 2.2.5-26 の灰色の部分のデータを使用した。
Fig. 2.2.5-27 に波数・周波数スペクトル解析結果を示す。屋根部分では流れ場（境界層）と音場の圧力が分離されているが、ドアラミー下流では剥離に伴う渦が見かけの音場を作るため、流れ場の圧力変動の一部が音場の波数に接近している。この現象は非圧縮解析を用いた分離解法でも確認されていたが、空力音の直接解析である LBM でも確認することができた。このことからドアラミー
後流のはく離はサイドシールド面にみかけの音場を作るため、車内音への寄与が大きくなると考えられる。

- 68 -

Fig. 2.2.5-26 Extraction of pressure data for wavenumber frequency analysis.

(a) Roof(b) Door Mirror

Fig. 2.2.5-27 Wavenumber and frequency spectrum.
2-2-6) 空力・振動・音響解析
流れに起因する車内騒音は、車体周りの流れ場、音場を加振力として車体パネルやガラスが振動し、その振動により車内に騒音が発生すると考えられる。そこで風洞に平板を設置し、平板下部に
車室内を模擬した空間を設け、平板に作用する流体力、音場が室内にどのような音場を形成するかを調べた。平板である窓のまわりの流れ場解析から窓に作用する音場と流体圧力場を算出し、窓の壁面に作用する圧力変動から振動及び車内騒音の解析を行った。Fig. 2.2.5-28 に流れ場と車内音の解析結果を示す。時間領域の流体解析、振動解析、音響解析を実施することより車内騒音を予測することが可能となった。Fig. 2.2.5-29 は実験結果との比較である。この条件では流体力学的な圧力変動（凡例の HPF、Hydrodynamical-Pressure Forcing）が支配的なため、振動、車内音ともに実験 結果 は流 体力 学的な圧 力に 起因 する ものと一 致し 、外 部の 空力音の 寄与 （凡 例の APF 、
Acoustical-Pressure Forcing）が小さいことがわかる。しかし、400 Hz 付近では空力音の寄与が大きくなることがわかる。このように周波数や条件によって寄与が異なる点が車内騒音の予測の難しい点であるが、数値解析よりその寄与を評価することができることを確認した。

- 69 -

Fig. 2.2.5-28 Flow field and internal sound around a flat plate.

Fig. 2.2.5-29 Flat plate vibration and internal sound due to flow and sound.
[1] 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,
200, 104427, (2020).
[2] H. Asada and S. Kawai, "Conservative low-pass filter with compact stencils for hierarchical
Cartesian mesh," Computers & Fluids, 252, 105769, (2023).
(6) プロジェクトの総合的推進
プロジェクト全体の連携を密としつつ円滑に運営していくため、プロジェクトの推進や実証研究
テーマ間の連携のための会議等を適宜開催し、参画の協力機関・連携機関との連携・調整にあたる。
特に、プロジェクト全体の進捗状況を確認し、計画の合理化の検討等を行うなど、プロジェクトの効果的・効率的推進に資する取組みを実施する。また、プロジェクトで得られた成果については、
可能な限り積極的に公開して今後の展開に資するとともに、ものづくり産業での早期戦力化を支援する。
令和 4 年度は、令和 3 年度に引き続き、上記を順次実施しながら、本課題の有効な研究活動のための支援を実施し、研究開発全体を円滑に実施するための調整を行った。そして、得られた成果については、HPCI ものづくり統合ワークショップや「富岳」流体予測革新プロジェクトシンポジウムを開催し、そこで公開した。
以下に代表的な活動を示す。

- 70 -

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

行

事

R4.6.10~11

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

R4.9.1~2

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

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

R4.11.14~17SC22

R4.12.21「富岳」成果創出加速プログラムシンポジウム

R5.3.1

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

R5.3.7~8「富岳」成果創出加速プログラム 研究交流会(2) 「富岳」流体予測革新プロジェクト推進会議
年月日

行

事

R4.6.15

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

R4.10.7

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

R4.12.8

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

R5.3.7

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

行事

第 3 回「富岳」高性能シミュレーション技術交流会(4) アプリケーション性能向上ワーキンググループ
年月日

行

事

R4.5.27

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

R4.7.28

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

R4.9.21

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

R4.11.22

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

R5.1.18

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

R5.3.31

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

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

担当機関

担当責任者

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

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

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

国立大学法人九州大学 大学院
工学研究院
国立大学法人東京大学 生産技術研究所(3) 圧縮機サージの直接
国立大学法人九州大学 大学院
解析
工学研究院
国立大学法人岩手大学 理工学部(4) リアルワールド自動 国立大学法人神戸大学 大学院
車空力性能の予測システム情報学研究科
国立研究開発法人理化学研究
所 計算科学研究センター
国立大学法人山梨大学 大学院
総合研究部(5) リアルワールド自動 国立大学法人豊橋技術科学大
車空力音予測学 大学院工学研究科
国立大学法人神戸大学 大学院
システム情報学研究科
国立研究開発法人理化学研究
所計算科学研究センター
国立大学法人東京大学 生産技術研究所(6) プロジェクトの総合 国立大学法人東京大学 生産技
的推進
術研究所

- 72 -

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

別添１ 学会等発表実績
１．学会等における口頭・ポスター発表
発表した成果（発表題目、
口頭・ポスター発表の別）

Digital transformation of
the manufacturing process
and smart design in theSociety 5.0 era realized onthesupercomputer“Fugaku” (Opening Plenary
talk)Engineering
Tubulence
Simulation on Massively
ParallelEnvironments
towardNext-Generation
Computer-Aided
Engineering(Invited
lecture)Digital transformation andsmart
designinmanufacturing
process
realizedonthesupercomputer“Fugaku”(Invited lecture)
CFD
Analysis
andVisualization of Cavitating
Bubbly Flow Structure in aCentrifugal Pump Using
Multi-Process
Cavitation
Model
Real-World Aerodynamics
Assessment for Road Vehicle
Development Realized ontheSupercomputer"Fugaku"「富岳」の時代の大規模流体解析の現状と将来展望
Engineering Applications ofLarge-scale
Fluid-flow
Computations「富岳」を用いた遷音速遠心
圧縮機のマイルドサージの大規模 DES 解析「富岳」の時代のシミュレー
ションとものづくり
多重プロセス型キャビテーシ
ョンモデルにおける輸送方程式に関する考察：遠心ポンプ
モデルを対象としたケースス
タディ

発表者氏名

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

発表し
た時期

国内・
外の別

Makoto Tsubokura

Supercomputing Asia 2022

2022 年
5月

国外

Makoto Tsubokura

The twelfth international
symposium on Turbulenceand Shear Flow Phenomena(TSFP12)2022 年
7月

国外

Makoto Tsubokura

The
11th
International
ConferenceofComputational
Fluid
Dynamics (ICCFD11)2022 年
7月

国外

Akihisa Yamada,
Taiki
Takamine,
Shin-ichi
Tsuda,
Satoshi Watanabe

13th Pacific Symposium on
Flow Visualization
and
Image
Processing(PSFVIP13)

2022 年
8月

国内

Ryo
Miyazaki,
Akinori
Yoshikawa, Rahul
Bale,
Kunihiko
Yoshitake,
Fortunato Nucera,
Takashi Yoshino,
Makoto Tsubokura
加藤千幸

World
CongressonComputational Mechanics,
WCCM-APCOM 2022

2022 年
8月

国外

第 19 回海洋エネルギーシン
ポジウム 2022（OES2022）
34th
Workbenchon
Sustained
Simulation
Performance
第 50 回日本ガスタービン学会定期講演会2022 年
9月
2022 年
10 月

国内

2022 年
10 月

国内

第 31 回日本国際工作機械見本市
日本機械学会第 100 期流体工学部門講演会2022 年
11 月
2022 年
11 月

国内

Chisachi Kato

真部魁人、古川雅人、草野和也、山田和豊、香西和人、兼山直輝加藤千幸
山田晃久、高峯大輝、津田伸一、渡邉聡- 73 -

国内

国内

実走行状態における空気抵抗の予測シミュレーション

自然風を考慮した自動車レーンチェンジ運動空力双方向連
成解析「富岳」を用いた大規模
DES による遷音速遠心圧縮機のマイルドサージ解析「富岳」を用いた大規模
DES による遷音速遠心圧縮機のディープサージ解析「富岳」の時代の HPC の産業応用と今後の展開

吉 川 暁 典 、 Rahul
Bale
、
Kim
Sangwon 、 吉 武 邦
彦 、 Fortunato
Nucera、吉野崇、坪倉誠
相 本 翔 太 、 Rahul
Bale
、
Kim
Sangwon 、 吉 武 邦
彦 、 Fortunato
Nucera、吉野崇、
坪倉誠
真部魁人、兼山直輝、古川雅人、草野和也、山田和豊
兼山直輝、真部魁人、古川雅人、草野和也、山田和豊
加藤千幸- 74 -

第 36 回数値流体力学シンポ
ジウム、C02-3

2022 年
12 月

国内

第 36 回数値流体力学シンポ
ジウム、C02-2

2022 年
12 月

国内

第 36 回数値流体力学シンポ
ジウム、C07-4

2022 年
12 月

国内

第 36 回数値流体力学シンポ
ジウム、C07-5

2022 年12 月

国内「富岳」EXPANDS～可能性を拡張する～

2023 年1月

国内