令和３年度高性能汎用計算機高度利用事業「富岳」成果創出加速プログラム「スーパーシミュレーションと AI を連携活用した実機クリーンエネルギーシステムのデジタルツインの構築と活用」
成果報告書

令和 4 年 5 月 30 日
国立大学法人東京大学
吉村 忍

-1-

本報告書は、文部科学省 令和 3 年度高性能汎用計算機高度利用事業費補助金 「富岳」成果創出加速プログ
ラム「スーパーシミュレーションと AI を連携活用した実機クリーンエネルギーシステムのデジタルツインの構築と活用」の成果をとりまとめたものです。

-2-

目次
１．補助事業の目的 ................................................................................................................................. 1
２．令和３年度（報告年度）の実施内容 ................................................................................................ 1
２－１．当該年度（令和３年度）の事業実施計画 .............................................................................. 1
２－２．実施内容（成果） ................................................................................................................... 2(1) 次世代火力発電システムのスーパーシミュレーションと V&V、デジタルツイン ................... 2(2) 洋上ウィンドファームのスーパーシミュレーションと V&V、デジタルツイン ..................... 23(3) アプリケーションの「富岳」における性能評価とチューニング ............................................ 46(4) シミュレーションと AI（機械学習）の連携活用とデジタルツイン ....................................... 56(5) 産業界等との連携及びコンソーシアムの構築.......................................................................... 70(6) 成果発表、広報・アウトリーチ活動 ........................................................................................ 73
２－３．活動（研究会の活動等） ...................................................................................................... 75
２－４．実施体制 ............................................................................................................................... 75
別添１ 学会等発表実績 ........................................................................................................................ 76i

補助事業の名称「富岳」成果創出加速プログラム
スーパーシミュレーションと AI を連携活用した実機クリーンエネルギーシステムのデジタルツインの構築と活用

１．補助事業の目的
ポスト「京」重点課題⑥プロジェクトにおいて「富岳」での活用を念頭に開発準備を進めてきたア
プリケーション群（ADVENTURE、FFB、FFR-Comb、REVOCAP_Coupler、RC HPC 版）に基づ
くマルチフィジクス・マルチスケール統合シミュレーション（本事業においてスーパーシミュレーシ
ョンと称する）と AI（機械学習）を連携活用し、Society5.0 を支える電力システムの主要クリーンエネルギーシステムとなる CCS（Carbon Capture and Storage）技術の適用に適した次世代火力発電システム、及び洋上ウィンドファームのデジタルツインを構築し、実機エネルギーシステムの開発、設計、運用の最適化に活用できるようにする。本事業において構築・活用されるアプリケーション群
は、いずれも我が国産業界における実活用の実績も豊富な高度に汎用化された力学系 CAE アプリケーションであり、上記のクリーンエネルギーシステムのみならず、燃焼機器、伝熱機器、流体機械等の幅広い実機に活用可能である。また、本事業は、様々な企業群、産業界コンソーシアムとも強く連携しながら遂行される。以上のことから、本事業は、「富岳」を活用した我が国の産業競争力の強化にも幅広く大きく貢献する。

２．令和３年度（報告年度）の実施内容
２－１．当該年度（令和３年度）の事業実施計画(1) 次世代火力発電システム及び洋上ウィンドファームのスーパーシミュレーションと V&V、デジタルツイン
•ポスト「京」重点課題⑥から構築を進めている、ベンチスケール石炭ガス化炉（三菱実験炉）の本格解析を実施し、V&V を実施する。
（東大・九大）•ポスト「京」重点課題⑥と⑧において、それぞれのコデザインのターゲットアプリであった(a)ADVENTURE 及び(b)FFB について、令和 2 年度に引き続き、「富岳」実機上で性能検証を継続し、スケーリングテスト実測を行う。
（東大・豊橋技科大）•ポスト「京」重点課題⑥において「富岳」での適用を念頭に、「京」及び名大 FX100、九大 ITO、
京大 Cray XC40 等において解析機能拡張及び高速化チューニングを行ってきた(a)FFR-Comb、(b) REVOCAP_Coupler に関して、「富岳」実機上での性能検証を行い、スケーリングテスト実測を行う。
（東大・京大・九大）•「富岳」を用いて、 (a)ベンチスケール石炭ガス化炉（三菱実験炉）、(b)パイロットスケー
ル超臨界圧 CO2 ガスタービン（東芝エネルギーシステムズ製）、(c)洋上ウィンドファーム（Alpha Ventus）の全体と(d)大型風車の後流を介した相互干渉のスーパーシミュレーションの解析を行い、V&V を実施する。（東大・京大・九大・豊橋技科大）(2) 個別シミュレーションと AI（機械学習）の連携活用•

各個別アプリケーション（(a)ADVENTURE、(b)FFB、(c)FFR-Comb、(d)RC HPC 版）の解析結果と AI（機械学習）を連携活用し、それぞれのサロゲートモデル構築およびデジタルツインの構築を行う。
（東大・京大・九大・豊橋技科大）(3) 総合的事業推進1

（3-1）産業界等との連携及びコンソーシアムの構築•

IEA Wind 国内委員会と連携し、IEA Wind Task 31 Phase III（ウィンドファーム流れのモデリングと評価の国際共同研究フレームワークのフェーズ III、2018～2021 年）に参画する。
（東大・九大・豊橋技科大）•燃焼システム用次世代 CAE コンソーシアムと本プロジェクト担当者らとの会合を開き、本プ
ロジェクトで開発中のスーパーシミュレーション技術の紹介を行うとともに、燃焼系産業界のニーズの聴収を行い、その結果を精査し、本シミュレーション技術の開発に反映する。
（東大・京大・九大）•

次世代 CAE を活用する風力産業コンソーシアムの立ち上げ作業を継続し、関係者らと本プロ
ジェクト担当者らとの会合を開き、本プロジェクトで開発中のスーパーシミュレーション技術の紹介を行うとともに、風力産業界のニーズの聴収を行い、その結果を精査し、本シミュ
レーション技術の開発に反映する。
（東大・九大・豊橋技科大）•次世代 CAE を活用する風力産業コンソーシアム、
日本風力発電協会、日本海事協会、IEA Wind、
IEA Wind 国内委員会との連携を通して、国内外への開発技術の普及体制を構築する。
（東大・九大・豊橋技科大）（3-2）開発アプリケーションの維持・普及及び成果広報（東大・京大・九大・豊橋技科大）•ADVENTURE、FFB、FFR-Comb、REVOCAP_Coupler の講習会等を開催し、産業界への技術の普及活
動を行う。

•

各分野トップレベルの国際学会での発表を行う。

•

各分野トップレベルの国際ジャーナルへ論文を投稿する。

•

本プロジェクト全体として国内シンポジウム（第２回、ポスト「京」重点課題⑥から数えると第６回目）を開催し、本プロジェクトのねらい、途中段階の成果を発表し、社会への広報・普及を行う。

２－２．実施内容（成果）(1) 次世代火力発電システムのスーパーシミュレーションと V&V、デジタルツイン(1-1) ラボスケール石炭ガス化炉(1-1-1) ADVENTURE_Thermal による伝熱・冷却管解析
ADVENTURE_Thermal（以下 ADV_Thermal）による伝熱・冷却管連成解析コード中の、連成境界面での物理量交換処理にバグがあり、一部の冷却管について不適切な冷却水の温度分布が出力される問題が発見されたため、令和 2 年度に引き続き改修を行った。
これまで使用してきた解析コードの分析を行い、冷却管から炉容器への熱移動処理で、移動方向が反対となるバグを発見し、改修した後、テストモデルによる試解析を行った。令和 2 年度では九州大学スーパーコンピュータ ITO を使って試解析を行っていたが、本番モデルへの適用を見据えて、令和 3 年度は「富岳」上で試解析を実施した。
使用したモデルは、令和 2 年度から引き続きテスト作業で使用したラボスケール石炭ガス化炉（電中研炉）の炉容器の 3 次元ソリッドメッシュ（四面体 1 次要素、要素数 155,999,061、
節点数 25,510,852）である。境界条件として、燃焼室内壁表面部分に燃焼に伴う熱移動を模した熱流束境界条件を与えた。
また、４本の冷却管モデルには流速一定の条件を課したほか、全ての冷却管の冷却水流入温度として24.77℃の温度規定境界条件を与えた。炉容器側の初期条件として FFR-ADV 連成解析で得られた
17,999 ステップ目の節点温度及び節点熱流束ベクトルを与えた。炉容器側の時間進展パラメータとして、Δt=0.001[s]を与えた。冷却管ソルバー（Discontinuous Galerkin[DG]法で離散化）内のサブ2

サイクルとして 1 時間ステップ当たり 1,000 もしくは 10,000 サブサイクルの計算を全ての冷却管に対して行うように設定した。連成解析のため、冷却管表面にソルバー内部で熱伝達（参照温度=流体
温度、熱伝達係数 0.062）を模擬した熱流束境界条件を動的に与えた。
試解析の結果、Fig. 1.1.1-1～Fig. 1.1.1-4 に示すような温度分布が得られた。いずれの図も青が改修前、オレンジが改修後のコードの計算結果である。改修によって、4 本の冷却管全てで冷却管流入口（グラフ左側）から流出口（グラフ右側）にかけて滑らかに温度が単調増加する、定性的に自然な温度分布が得られるようになった。また、これまでは計算途中で異常終了していたが最後まで正常に動作するようになった。

Fig. 1.1.1-1. Temperature distributions in cooling pipe (a) at t=0.616s.

Fig. 1.1.1-2. Temperature distributions in cooling pipe (b) at t=0.616s.

3

Fig. 1.1.1-3. Temperature distributions in cooling pipe (c) at t=0.616s.

Fig. 1.1.1-4. Temperature distributions in cooling pipe (d) at t=0.616s.
加えて、これまでの解析システムは安定性に問題があり細かなΔt の使用が必要であった。炉容器の熱伝導解析の際、冷却水によって奪われる熱量はパイプ表面の温度と冷却水の温度を引数としニュートンの冷却法則を用いて計算されるが、パイプ表面の温度は陽的に処理されていた。この処理が安定性を低下させていることを数学的に明らかにし、陰的に扱うように実装を変更した。コンクリー
トの冷却問題を解くことで、改良後の解析システムの妥当性が確認された。Fig. 1.1.1-5 はジオメトリを表しており、中心に冷却水を流すことで周囲のコンクリートを冷却するという問題設定となっ4ている。改良後の解析システムによる数値解は、先行研究 [1]で示された解析解とよく一致している
ことが、Fig. 1.1.1-6 より確認できる。

Fig. 1.1.1-5. Geometry of water cooling problem for verification.

Fig. 1.1.1-6. Time history of averaged temperature in concrete.
文献[1] B. Zhu, Effect of cooling by water flowing in nonmetal pipes embedded in mass concrete, J.
Constr. Eng. Manage. 125 (1) (1999) 61-68.
(1-1-2) FFR-Comb による約 1 億要素の詳細格子における流体－構造連成解析(1-1-2-1) 解析方法
Eulerian-Lagrangian 法による固気二相反応流の Large-eddy simulation[1,2]を行う。気相の支配方程式は、連続の式、運動量保存式、エネルギー保存式、および化学種保存式を解いた。SGS 応力モ
デルには Dynamic Smagorinsky モデル[3]を用いた。固相については、個々の粒子のいわゆる BBO
方程式をパーセル近似により解き、流体と粒子の干渉効果については、PSI-CELL 法[4]により考慮した。ガス化炉内の化学反応としては、石炭の熱分解による揮発分放出、石炭チャーの固気（ガス化）反応、および気相反応を考慮した。揮発分放出過程を Eq. (1)により表し、揮発分組成については、
Umemoto らの実験データ[5]を用いた。揮発分放出速度については、Eq. (2)に示す単一一次反応モデ
ルにより考慮し、その反応速度パラメータは FLASHCHAIN モデル[6]により見積もった。

5

𝐶𝑜𝑎𝑙 → 𝛼1 𝐶𝐻4 + 𝛼2 𝐻2 + 𝛼3 𝐶𝑂 + 𝛼4 𝐶𝑂2 + 𝛼5 𝐻2 𝑂 + 𝛼6 𝐶ℎ𝑎𝑟

(1)

𝑑𝑉𝑖

(2)

𝑑𝑡

𝐸𝑖

= 𝐴𝑖 exp (− 𝑅𝑇 ) (𝑉𝑖∗ − 𝑉𝑖 )𝑝

ここで、𝛼𝑖 は総和が 1 となる各揮発分の係数である。𝑉𝑖 は揮発した各揮発分量、𝐴𝑖 および𝐸𝑖 は頻度
因子と活性化エネルギー、𝑉𝑖∗は粒子が保有する総揮発分量である。
チャーガス化反応については、Eqs. (3)〜(5)に示す O2 ガス化、CO2 ガス化、および H2O ガス化の3 つの反応経路を考慮した。チャーガス化反応モデルには、CO2 ガス化および H2O ガス化について
は Eq. (6)に示す Langmuir-Hinshelwood モデルに基づく活性点部分共有モデル[5]を用いた。O2 ガ
ス化については nth-order モデル[7]を用いた。それぞれの反応経路の反応速度パラメータは文献[8]の値を用いた。
1

𝐶ℎ𝑎𝑟 + 𝑂2 → 𝐶𝑂

(3)𝐶ℎ𝑎𝑟 + 𝐶𝑂2 → 2𝐶𝑂𝐶ℎ𝑎𝑟 + 𝐻2 𝑂 → 𝐶𝑂 + 𝐻2

(4)
(5)

2

𝑑𝑥
𝑑𝑡

= 𝜂(1 − 𝑥)√1 − Ψln(1 − 𝑥) (1+𝑘𝑘11 𝑃𝐶𝑂2

12 𝑃𝐶𝑂2 +𝑘13 𝑃𝐶𝑂 +𝑎 ⁄𝑐 𝑘22 𝑃𝐻2𝑂 +𝑎 ⁄𝑐 𝑘23 𝑃𝐻2
𝑘21 𝑃𝐻2𝑂

1+𝑏𝑐𝑘12 𝑃𝐶𝑂2 +𝑏𝑐𝑘13 𝑃𝐶𝑂 +𝑘22 𝑃𝐻2𝑂 +𝑘23 𝑃𝐻2

+

)

(6)

ここで、𝑥は粒子内炭素分の反応率、𝜂は粒子内細孔に対する有効拡散係数、Ψは細孔の形態係数、
𝑃𝑖 は各化学種分圧、𝑎、𝑏および𝑐はそれぞれ活性点部分共有効果を表すモデルパラメータである。
気相反応については、Eqs. (7)〜(11)の総括反応を考慮した。気相反応速度の算出にはアレニウス
式の直接クロージャモデルである SSFRRM[9]（Eq. (12)）を用いた。
1(7)

𝐻2 + 𝑂2 → 𝐻2 𝑂

(8)

𝐶𝑂 + 𝑂2 → 𝐶𝑂2(9)

𝐶𝐻4 + 𝐻2 𝑂 → 𝐶𝑂 + 3𝐻2
𝐶𝑂 + 𝐻2 𝑂 ↔ 𝐶𝑂2 + 𝐻2

(10)
(11)

̿̿̿̿̿̿̿̿̿̿̿̿̿̿̿
̿̿̿̿̿̿̿̿̿̿̿̿̿
̅̅̅𝜔̇ 𝑖 = ̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜔̇ (𝜌̅ , 𝑌̃𝑖 , 𝑇̃) + 𝐾1 (𝜔̇(𝜌̅ , 𝑌̃𝑖 , 𝑇̃) − 𝜔̇ (𝜌̿ , 𝑌̃𝑘 , 𝑇̃))

(12)

𝐶𝐻4 + 2 𝑂2 → 𝐶𝑂 + 2𝐻2
1

2
1
2

̅̅̅𝑖 は各反応物の消費速度である。
ここで、𝜔̇
輻射モデルには Discrete Ordinate 法[10]を用いた。
(1-1-2-2) 解析条件
Fig. 1.1.2-1 に、対象とするラボスケール石炭ガス化炉の計算領域ならびに計算格子（バーナ部拡大）を示す。格子数は詳細格子である約 1 億 3 千万要素のものを用いた。ガス化炉の試験条件（解析条 件 ） は 、 石 炭 供 給 量 100 kg/h 、 酸 素 比 0.54 の 条 件 に 対 し て 、 ガ ス 化 剤 組 成 を
O2/CO2/N2=25 %/25 %/50 %に設定した。NS 方程式の対流項の離散化には 2 次中心差分に 5 %の 1
次風上成分をブレンドしたものを用いた。時間進行にはオイラー陰解法を用い、時間刻みは 5×10-6 sとした。
(1-1-2-3) 解析結果
Fig. 1.1.2-2 に、連成計算の温度・化学種分布、Fig. 1.1.2-3 に単体計算の温度・化学種分布を示す。
連成計算では単体計算に比べて温度がやや低下し、化学種分布にも差異が現れていることがわかる。

6

また、Fig.1.1.2-4 に示す中心軸上ガス温度分布から、z=0.8 m より下流のリダクタ―部において連成計算が単体計算よりも実験値に近づいていることがわかる。
Fig. 1.1.2-5 に、連成計算および単体計算による出口ガス組成の計算結果を示す。出口ガス組成は単体計算の方が実験値と近い結果となったが、原因については引き続き調査と詳細な現象理解を進める予定である。

(a) Computational
domain.

(b) Computational grid.
Fig. 1.1.2-1 Computational domain and grid.

(a)(b) H2(c) CO(d) CO2(e) H2O
Temperature
Fig. 1.1.2-2 Gas temperature and concentrations of major gases(Coupling Analysis) .

7

(a)(b) H2(c) CO(d) CO2(e) H2O
Temperature
Fig. 1.1.2-3 Gas temperature and concentrations of major gases(Fluid Analysis) .

Fig. 1.1.2-4 Comparison of gas temperature distribution in the axial direction.

8

Fig. 1.1.2-5 Comparison of product gas composition at the exit.

文献[1] H. Watanabe, K. Tanno, Y. Baba, R. Kurose, S. Komori, Large-eddy simulation of coal
combustion in a pulverized coal combustion furnace with a complex burner, Turbulence, , 6,
pp.1027-1030, 2009.
[2] M. Muto, H. Watanabe, R. Kurose, S. Komori, S. Balusamy, S. Hochgreb, Large-eddy
simulation of pulverized coal jet flame -Effect of oxygen concentration on NOx formation, Fuel,142, pp.152-163, 2015.
[3] P. Moin, K. Squires, W. Cabot, S. Lee, A dynamic subgrid-scale model for compressible
turbulence and scalar transport, Phys. Fluids A, 3, pp.2746-2757, 1991.
[4] C.T. Crowe, M.P. Sharma, D.E. Stock, The Particle-Source-In Cell (PSI-CELL) model for gasdroplet flows, J. Fluids Eng., 99, pp.325-332, 1977.
[5] S. Umemoto, S. Kajitani, S. Hara, Modeling of coal char gasification in coexistence of CO2
and H2O considering sharing of active sites, Fuel, 103, pp.14-21, 2013.
[6] Niksa, S., A.R. Kerstein, FLASHCHAIN theory for rapid coal devolatilization kinetics, 1.
Formulation, Energy Fuels, 5, pp.647-665, 1991.
[7] H. Watanabe, M. Otaka, Numerical simulation of coal gasification in entrained flow coalgasifier, Fuel, 85, pp.1935-1943, 2006.
[8] H. Watanabe, K. Tanno, H. Umetsu, S. Umemoto, Modeling and simulation of coalgasification on an entrained flow coal gasifier with a recycled CO2 injection, Fuel, 142, pp.250259, 2015.
[9] P.E. DesJardin, S.H. Frankel, Large-eddy simulation of a nonpremixed reaction jet:application and assessment of subgrid-scale combustion models, Phys. Fluids, 10, pp.2298-2314,
1998.
[10] W.A. Fiveland, Three-dimensional radiative heat-transfer solutions by the discreteordinates method, J. Thermophys., 2, pp.309-316, 1988.
(1-2) パイロットスケール超臨界圧 CO2 ガスタービン(1-2-1) 目的
令和 3 年度は、令和 2 年度までに確立した乱流燃焼モデルの一種である Flamelet 法に AI（機械9学習）を連携させた Flamelet/ANN（Artificial Neural Network）法の高精度化を行うとともに、本手法を用いた LES（以下、LES/Flamelet/ANN と呼ぶ）をパイロットスケール（商用機の 10 分の 1程度）の東芝エネルギーシステムズ株式会社製 50MWth 級超臨界 CO2 タービン用燃焼器内の乱流燃焼場[1-4]に適用し、試験結果との比較することによって、同手法の有用性を確認することを目的とした。令和 2 年度は HPCI 第二階層マシンを用いた燃焼器上流部のみを対象とした計算を行ったが、令和 3 年度はスーパーコンピュータ「富岳」を用いた燃焼器全体の計算を実施した。
(1-2-2) 計算方法
① 燃焼モデルと支配方程式
本解析では、
Flamelet 法として火炎の浮き上がりや消炎を再現可能な Flamelet/progress-variable（FPV）法[5,6]をベースとし、壁面からの熱損失および実在実気体を考慮可能に拡張したものを用いた[7]。支配方程式は、連続の式、運動量の保存式（Navier-Stokes 方程式）
、混合分率 Z の保存式、
反応進行変数 C の保存式、およびエンタルピーh の保存式で構成される。
② 3 種吹込みガス燃焼の Flamelet 法
従来の Flamelet 法では、燃料と酸化剤の 2 つの独立したガス系統からの流入混合による燃焼のみに対応可能で、燃料、酸化剤に加えて冷却用 CO2 の吹込み等、3 つ以上の独立したガス系統からの吹込みがある条件に適用することができなかった。そこで、令和 2 年度に、3 種吹込みガス条件に適用可能な Flamelet 法の導入を行った[8]。その概念を Fig.1.2-1 に示す。

Fig. 1.2-1. Concept of three stream Flamelet mode.
従来の Flamelet 法では、燃料由来成分の質量分率を示す混合分率 Z を輸送することで、燃焼場における燃料と酸化剤の混合具合を求め、その値を基に Flamelet データベースを参照し、燃焼場を再現する。これに対して、本手法では、燃料由来の質量分率（図中の Fuel）を Z1、Cooling CO2 の質
量分率 （図中の Coflow）を Z2 と定義し、これら 2 つの輸送方程式を解くことで、3 系統の流入ガ
スの混合具合を判断する。今、酸化剤の質量分率を Z0 として、Cooling CO2 (Coflow)を酸化剤の一部と定義すると、酸化剤に含まれる Cooling CO2 (Coflow)の割合 W は、
𝑊=

𝑍2

(a)

𝑍0 +𝑍2とおけるので、3 種吹込みガス燃焼の FPV 法の場合、Flamelet データベースの物理量φは下記で与えられる。
"2 , 𝑊
̃ , 𝐶̃ , 𝛥𝐻̃)𝜙̃ = 𝜙̃(𝑍̃, 𝑍̃(b)

ここで、上付き～はファーブル平均、”は変動値であり、ΔH は壁面からの熱損失を考慮するために用いるエンタルピー差である。また、反応進行変数 C は、(c)

𝐶 = 𝑌𝐻2𝑂 + 𝑌𝐻2
で与える。

10

③ ニューラルネットワーク（ANN）による Flamelet モデルの構築
Flamelet データベースを引用するパラメータの増加に対応するため、多次元の Flamelet 空間で"2 , 𝑊
̃ , 𝐶̃ , ∆𝐻̃ ）の関数で表現するの多様体であるデータベースを基底ベクトル（本解析においては𝑍̃, 𝑍̃
ことを考える。こうすることで流体計算時の必要なメモリがデータベースそのものからパラメータベクトルだけとなり、メモリの大幅な節約となる。即ちこれは Flamelet データベースに対する回帰問題となる。
機械学習という枠組みの中でも回帰問題を解くための手法は多く存在するが、ここではディープ
ラーニングの手法、即ちニューラルネットワークを用いた関数の表現を採用する。本 5 次元 Flamelet
データベースにニューラルネットワークを適用した場合、従来は、入力層のパーセプトロンの数はデータベース空間の基底ベクトル数である 5、出力層のパーセプトロンの数は 1 となり、物理量ごとにANN のデータセットを用意することになる。隠れ層をどの様に設計するかはハイパーパラメータであってそれ自体が最適化の対象となり得る。一般的に層が深いほど、層あたりのパーセプトロンが多いほど再現可能な関数形状は多くなり、データベースの再現性が高くなることが考えられる。また、入力層のパーセプトロンの数が多くなると、各対称条件下において十分な再現精度を保つためには、より多くの層数、各層当たりのパーセプトロンが必要となる。しかし一方で、ネットワークを大きくすれば燃焼計算においてデータベースの基底ベクトルから燃焼状態を推定する際の計算負荷が高くなることから、ネットワークの構築においてはデータベースの一定の再現性を担保しつつもネット
ワークができるだけ小さくなるように留意した。本解析では、本来の Flamelet データベースの 5 次"2 , 𝐶̃ , ∆𝐻
̃ の 4 次元のデータベースに落とし込み、各𝑊
̃ の ANN データベースを
元データベースを、𝑍̃, 𝑍̃̃ 軸に対しては 2 つの ANN データベースを補間することで、再現精度を維持しつつ、ネッ
用意し、𝑊
トワークの複雑さを抑えた。各学習の収束条件として、データベースの値とネットワークによって作成した関数の値の RMSE が、データベースの物理量の標準偏差の 4.5%未満となるようにした。この条件の下で試行錯誤によって、ネットワーク構成は隠れ層が 7 層、１層当りのパーセプトロンは 9 とした。各パーセプトロンの出力は活性化関数 Swish[9]で変換され、オプティマイザーRadam[10]によって重みが更新される。
④ 計算対象および計算格子

Fig.1.2-2 に、本解析で対象とした超臨界燃焼器の概略を示す。対象としたのは、東芝エネルギーシステムズ製 50MWth 級検証試験用超臨界 CO2 タービンおよび燃焼器[1-4]である。本ガスタービン
燃焼器は、酸化剤/燃料(CH4)/冷却用 CO2 により燃焼が行われる。運転条件はタービン入口圧力約30MPa の超高圧条件を想定しているが[4]、今回対象とした燃焼試験は設備容量の制限から最大10MPa で行っている。Fig.1.2-3 に、解析格子形状の詳細の様子を示す。点数およびセル数は、それ
ぞれ 6,160 万点および 1.92 億点である。

11

Fuel
NozzleMeasurement
Point

Liner

Transition
piece
Flow Direction

Fig.1.2-2. Schematic of real scale combustion test rig.

(a) Computational domain(b)Grids
Fig.1.2-3. Targeted combustor.
⑤ 計算条件12

本解析では、NuFD/FrontFlow/Red（FFR-Comb）[5-7,11,12]を用い、乱流モデルとしては、標準
Smagorinsky モデルを用いた。解法には圧縮性流体スキームを用い、実在気体モデルとして SoaveRedlich-Kwong (SRK) モデル[13]を用いた。
Table1.2-1 に 、 Flamelet デ ー タ ベ ー ス の 概 要 を 示 す 。 各 W の デ ー タ ベ ー ス サ イ ズ は 、"2 × 𝐶̃ × ∆𝐻
̃ = 100 × 25 × 100 × 5の分割の場合、従来の Flamelet データベースで 0.58GB であ𝑍̃ × 𝑍̃
り、W の分割を 9 (W=0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,1.0)とした場合、データベース容量は 5.22GB となり通常の大規模計算機の 1 コア当たりのメモリ容量を超え、解析が不可能となることが予想される。一方、メモリ容量の問題を解決するため ANN 法を用いた機械学習を用いた場合、データベースに必要なデータは各層間の重み付け関数パラメータのみで、データ容量がそれぞれ 0.24×10-3GB、2.16×10-3GB の約 1/2400 と大幅にデータベース容量を削減することができることが確認できる。
Table1.2-1. Comparison of flamelet database.

Memory (GB)NW

Conventional database

ANN database

1

0.58

0.24×10-3

9

5.22

2.16×10-3(1-2-3) 結果および考察
Fig.1.2-4 に、本 LES/Flamelet/ANN で得られた、燃焼器内瞬時温度の等値面分布および垂直断面分布を示す。また、Fig.1.2-5 には、燃焼器内瞬時温度、流速、混合分率、冷却用 CO2、および化学
種質量分率（CH4, O2, CO2, CO, H2O）の中心断面分布を示す。図から、以下の挙動が確認できる。
・Fig.1.2-5(c)は奥行き速度を示す。上半分と下半分で色（正負）が異なるので、流れ場が旋回していることがわかる。燃焼器上流に設置された燃料ノズルより流入する酸化剤の一部は、外周側のスワラを通ることにより旋回速度が与えられ、燃焼器内の旋回流を形成する。
・Fig.1.2-5 (g)に示すように、燃料ノズルから流入したメタン燃料は、燃料ノズル内周から流入した少量の酸化剤と反応し即座に消費される。
・Fig.1.2-5 (i)に示すように、燃料ノズル出口中心部分に存在する燃焼領域では、燃料と酸化剤とが燃料過濃条件で燃焼し、大量の CO が生成される。その後、燃料ノズル外周より流入した酸化剤と混合反応することで、CO は速やかに消費される。
・Fig.1.2-5 (b)において、最上流の緑色の領域は負の速度を示しており、燃料ノズル出口下流で逆流領域が形成され、この領域で保炎される。
・Fig.1.2-5 (a)の下流部に示すように、壁面から冷却用 CO2 が注入されるため、壁付近から温度が下がっている。高温の燃焼ガスは、ライナやトランジションピースから流入する冷却用 CO2 により希釈され、燃焼器出口に流出される。また、燃焼器の側面から供給される冷却用 CO2 により壁面近傍では温度が低く保たれる。
Fig.1.2-6 に、出口温度の実験値との比較(ΔT=TLES-TExp)および温度分布を示す。ここで TLES は本
LES/Flamelet/ANN で得られた温度解析結果、TExp は実験で得られた温度である。温度分布の比較から、本 LES/Flamelet/ANN で得られた温度は試験結果を全計測位置において 50℃以内と比較的良好な精度で予測できていることが確認できる。また、本解析は、スーパーコンピュータ「富岳」を用いて 20,000 並列の解析を行い、統計計算を含めた計算を約 3 日で完了した。
以上の結果から、本 LES/Flamelet/ANN が超臨界燃焼場を高精度に予測するうえで有効な手法であると言える。

13

(a) Iso-surface(b) Cross-sectional planes
Fig.1.2-4. Instantaneous distributions of temperature.

Fig.1.2-5. Instantaneous distributions of temperature, velocities, mixture fraction, recycled
CO2 mixture fraction, and mass fractions of CH4, O2, CO2, CO, H2O on combustor central plane.

14

(a)Temperature differences (ΔT=TLES-TExp.) between simulated (TLES) and measured (TExp.)
combustor outlet temperatures(b)Simulated combustor outlet temperature distribution
Fig.1.2-6. Comparison of combustor outlet mean temperature distribution.
(1-2-4) まとめ
令和 3 年度は、令和 2 年度までに有効性を示した LES/Flamelet/ANN を、パイロットスケールの東芝エネルギーシステムズ株式会社製 50MWth 級超臨界 CO2 タービン用燃焼器内の乱流燃焼場全領域に適用し、試験結果と比較することによって、本手法の精度検証を行った 。その結果、
LES/Flamelet/ANN は、スーパーコンピュータ「富岳」を用いて 20,000 並列の計算を行うことにより、燃焼器全領域の大規模乱流燃焼解析を 3 日程度の短期間で完了でき、燃焼器の設計や最適操作条
件の選定において、極めて有効な手段になりうることを確認した。
文献[1] Sasaki, T., Itoh, M., Maeda, H., Tominaga, J., Saito, D., Niizeki, Y., “Development of Turbine
and Combustor for a Semi-closed Recuperated Brayton Cycle of Supercritical Carbon Dioxide”, In
Proc. of the ASME 2017 Power Conference Joint With ICOPE-17, POWER-ICOPE2017-3419,
2017.
[2] Allam, R., Fetvedt, J., Forrest, B., Jones, C., Nomoto, H., Itoh, M., “A novel, high-efficiency,
oxy-fuel power plant with low-cost electricity production and 100% capture of carbon dioxide”,
POWER-GEN International 2013.
[3] 岩井保憲, 伊東正雄, “超臨界 CO2 サイクル発電用ガスタービン燃焼器”, 東芝レビュー, 70 (5),2015.
[4] Iwai, Y., Itoh, M., Morisawa, Y., Suzuki, S., Cusano, D., Harris, M., “Development approachto the combustor of gas turbine for oxy-fuel, supercritical CO2 Cycle”, In Proc. of the ASME Turbo
Expo 2015, Palais des Congres, Montreal, Canada, GT2015-43160, 2015.
[5] Moriai, H., Kurose, R., Watanabe, H., Yano, Y., Akamatsu, F., Komori, S., “Large-eddy
simulation of turbulent spray combustion in a subscale aircraft jet engine combustor - Predictionsof NO and soot concentrations -", Journal of Engineering for Gas Turbines and Power, 135, 091503,
2013.
[6] Tachibana, S., Saito, K., Yamamoto, T., Makida, M., Kitano, T., Kurose, R., “Experimental
and numerical investigation of thermo-acoustic instability in a liquid-fuel aero-engine combustor
at elevated pressure: validity of large-eddy simulation of spray combustion", Combustion and
Flame, 162, 2621-2637, 2015.

15

[7] Kishimoto, A., Moriai, H., Takenaka, K., Nishiie, T., Adachi, M., Ogawara, A., Kurose, R.,“Application of a non-adiabatic flamelet/progress-variable approach to Large Eddy Simulation of
H2/O2 combustion under a pressurized condition”, Journal of Heat Transfer, 139, 124501, 2017.
[8] Ihme, M., and See, Y. C., “LES flamelet modeling of a three-stream MILD combustor:
Analysis of flame sensitivity to scalar inflow conditions”, Proceedings of the Combustion Institute,33, 1309-1317. 2011.
[9] Prajit, R., Barret, Z., Quoc V. L., “Swish: a Self-Gated Activation Function”,
arXiv:1710.05941, 2017.
[10] Liyuan, L., Haoming, J., Pengcheng, H., Weizhu, C., Xiaodong L., Jianfeng, G., Jiawei, H.,“On the Variance of the Adaptive Learning Rate and Beyond”, arXiv:1908.03265, 2019.
[11] Kobayashi, Y., Iwai, Y., Itoh, M., Nishiie, T., Kurose, R., “Large-eddy simulation of fluidstructure interaction for Supercritical CO2 combustor”, The Asian Congress on Gas Turbines 2018
(ACGT2018), Aiina Center, Morioka, Iwate, Japan, 2018.
[12] Jain, P., Iwai, Y., Kobayashi, Y., Itoh, M., Nishiie, T., Kurose, R., “Large-eddy simulation ofa supercritical CO2 combustion field in a realistic combustor”, The 1st R-CCS International
Symposium, Kobe International Conference Center, Kobe, Japan, 2019.
[13] Soave, G., “Equilibrium Constants from a Modified Redlich-Kwong Equation of State,”
Chemical Engineering Science., 27, 1197–1203, 1972.
(1-3) ベンチスケール石炭ガス化炉の燃焼解析
ベンチスケール石炭ガス化炉はラボスケール石炭ガス化炉に比べて約 30 倍の容積比であることから、1 億要素ラボスケール石炭ガス化炉と同等の格子解像度を得るために約 30 億要素のベンチスケ
ール炉解析格子を作成し試計算を実施した。
(1-3-1) 解析格子の作成
ベンチスケール炉 30 億要素格子の作成には、FFR-comb の refiner 機能を使用した。FFR-combの refiner 機能は、対象となる解析格子の格子幅を 1/2 とするため、ベースとなる解析格子の要素数
を 8 倍とすることができる。今回の対象となるベンチスケール炉は 30 億要素規模であることから、
5,000 万要素の格子を作成の上、refiner 機能を２回実施することで 64 倍の 32 億要素格子を作成した。
(1-3-2) 解析方法
解析方法は、(1-1-2-1) 節「解析方法」と同様である。ただし、下記のメタン水蒸気改質反応について、順反応に加えて逆反応も考慮するようにコードを修正した。
（旧）CH4 + H2O ⇒ CO + 3H2(1)（新）CH4 + H2O ⇔ CO + 3H2(2)(1-3-3) 試解析の実施
作成したベンチスケール炉解析格子を使用した単体解析の試解析を実施し、FFR-comb による 30
億要素ラボスケール炉の単体解析が可能であることを確認した。解析時の並列数は、１並列あたりの要素数が１億要素ラボスケール炉と概ね同等となるよう 65,536 並列とした。また、微粉炭粒子を含む固気二相流解析に必要なメモリ量を確保するために 1 ノードあたり 24 コアを使用することとした。したがって、計算に必要なノード数は 2,731 ノードである。試解析の結果、1,000step の解析に必要な計算時間が約 2 時間であり、1step あたり約 7 秒という結果となった。

16

本試解析の結果から、30 億要素ベンチスケール炉の解析に必要な解析時間およびノード時間積を把握することができ、ベンチスケール炉の流体－構造連成解析に向けた知見を得ることができた。
(1-4) ベンチスケール石炭ガス化炉の伝熱・構造解析(1-4-1) 静弾性解析
ベンチスケール炉モデルの構造解析における収束特性の概要を知るため、Table 1.4.1-1 に示した諸条件で静弾性解析を行った。
Table 1.4.1-1 Simulation conditions for static elastic analysis
メッシュ

四面体 1 次要素 (要素数 406,317,033、節点数 67,112,837)

境界条件

変位 : 底面（z= 1500.25）の節点全てを x/y/z3 方向固定
荷重 : 節点温度分布（Fig. 1.4.1-1）による熱応力及び自重

物性値

全要素に高張力鋼 HT80 を想定した物性を設定。
・ヤング率 : 203 [GPa]・ポアソン比 : 0.3・質量密度 : 7.86e-6 [kg/mm3]
・重力加速度 : (0.0, 0.0, -9800.0) [mm/s2]・熱膨張係数 : 11.7e-6 [1/K]
・参照温度 : 26.85 [℃]

領域分割ツール ADVENTURE_Metis Ver.1.1 及び構造解析ソルバー ADVENTURE_Solid Ver.
2.11 （公開版）を使用して、「富岳」30 ノードで計算を行った。領域分割は、1 部分領域の自由度数
を 300-500 程度になるように調整した結果、1,440 部分×300 部分領域（部分領域総数 432,000）に分割した。ソルバーはフラット MPI 並列版（advsolid-p）を 1,440 プロセス並列で実行し、約 5 時間で計算を完了した。
計算結果は LexADV_WOVis を使用して、内部の応力分布がわかるように y<0 の範囲の要素のみ可視化した。計算結果として得られた節点相当応力の最大値は約 45GPa であった。そこで当初は色コンターの範囲を 0-50,000MPa として可視化したところ、炉容器全体の応力分布がほぼ一様になってしまった（Fig. 1.4.1-2）。これは、一部に特異点的な応力集中があり、それ以外の箇所については少し低いオーダーの値となっているためと推測された。そこで、特異点の値を排除した色コンターの範囲を手動で設定するため、相当応力の分布をヒストグラム化した。ヒストグラムは、計算結果のう
ちの節点相当応力を ADVENTURE_Solid 付属の変換ツール hddmmrg を使ってテキスト化した後、ヒストグラム分析ツールによって 100MPa 刻みでカウントした CSV ファイルを作成し、Microsoft
Excel でグラフ化した（Fig. 1.4.1-3）。その結果、99.9%の節点が 4,000MPa 以下の値をとっている
ことが分かったため、色コンターの範囲を 0-4,000MPa として、再度可視化を行った（Fig.1.4.1-4）。
上部リダクタ部と下部ディフューザ・コンバスタ・スラグホール部の応力が比較的高かったため、これらの箇所について炉壁との接続部分がわかるように視点を変えて再度可視化を行った。Fig. 1.4.15 に下部ディフューザ・コンバスタ・スラグホール部、Fig. 1.4.1-6 にリダクタ部の可視化結果をそれぞれ示す。なお、リダクタ部の左側の内壁側および外壁側に高応力部があるが、これは Fig. 1.4.17 に示した通り冷却管流出口が表示されていない y>0 の領域にあり、外壁部と内壁部の温度差から熱応力が発生しているためである。

17

Fig. 1.4.1-1. Distribution of temperature
condition for thermal-stress analysis.

Fig. 1.4.1-2. Distribution of nodal
equivalent stress (range: 0-50,000MPa).

18

9,000,000
8,000,000

6,000,000
5,000,000
4,000,000
3,000,000
2,000,000

1,000,000
0

0
200
400
600
800
1,000
1,200
1,400
1,600
1,800
2,000
2,200
2,400
2,600
2,800
3,000
3,200
3,400
3,600
3,800

Number of nodes

7,000,000

NodalEquivalentStressFig. 1.4.1-3. Histogram of nodal equivalent stress(x axis=nodal equivalent stress(Unit:MPa), y axis=number of nodes).

19

Fig. 1.4.1-4. Distribution of nodal
equivalent stress (range: 0-4,000MPa).

Fig. 1.4.1-5. Distribution of nodal
equivalent stress around diffuser, combustor,and slugholes (range: 0-4,000MPa).

20

Fig. 1-4-1-6. Distribution of nodal equivalent stress around reductor (range: 0-4,000MPa).
冷却水流出口

冷却水流出口
冷却水流出口
応力分布

応力分布
応力分布

温度分布
温度分布
リダクタ部のCAD形状
リダクタ部の(アライドエンジニア
CAD 形状リング様資料より引用)

冷却水流入口
冷却水流入口

Fig. 1-4-1-7. Relationship of locations of inlet and outlet of reductor cooling pipe.
(1-4-2) 連成非線形構造解析
本項目ではベンチスケール炉の連成非線形構造解析を行うための予備的な構造解析を行った。ベ
ンチスケール炉を構成する部材としてメッシュデータに鋼材・TOCAST-13L・HE0152C の 3 種類のボリュームが設定されている。それぞれについて弾性物性値を Table 1.4.2-1 のように設定した。
Table 1.4.2-1. Elastic material properties of bench scale reactor.
TOCAST-13L
HE0152C
鋼材
ヤング率[MPa]

1.93E+05

1.96E+04

2.00E+05

ポアソン比

0.29

0.18

0.294

線膨張係数[1/K]

1.87E-05

7.00E-06

1.70E-05

21

ここで、鋼材は SUS304 を想定した値となっている。TOCAST-13L はキャスタブル耐火物となっており、粉末の耐火物を混練し成形することでガス化炉の耐熱部に使用されているものと考えられる。TOCAST-13L は組成についての情報が見つからなかったため、本解析では文献[1]の耐火レンガの物性値を使用することとした。HE0152C に関しては Ni 耐熱合金と仮定して Inconel718 の物性値を使用した。
鋼材については SUS304 の Ohno-Wang モデルとクリープ物性を材料非線形データとして設定した。初期降伏応力を Table 1.4.2-2、Ohno-Wang の移動硬化パラメータを文献[2]を参考にしてパラメータ設定を Table 1.4.2-3 のように行った。
文献[3]を参考にしたクリープパラメータを Table 1.4.24 に示す。
Table 1.4.2-2. Initial yield stress of SUS304.
2.550E+02
初期降伏応力[MPa]Table 1.4.2-3. Ohno-Wang parameters of SUS304.
r [MPa]
モード
ζ
1
2
3
4

8.000E+03
4.000E+03
3.333E+03
2.000E+03

7.351E+00
2.953E+00
1.531E+00
1.492E+00

Table 1.4.2-4. Creep parameters of SUS304.
A
n
m
温度[℃]
20～7450.000E+00

1.000E+00

0.000E+00

750
900
1000
1200

2.750E-19
4.020E-14
2.030E-09
6.100E-07

6.250E+00
5.200E+00
3.000E+00
2.500E+00

0.000E+00
0.000E+00
0.000E+00
0.000E+00

Δt=1.0s として、各ステップの温度増分を 100℃として均一になるような昇温過程の温度データ
を 4 Step 作成して、「富岳」上でベンチスケール炉の非線形熱応力解析を行った。別途、SUS304 に
弾性物性値を設定し線形の熱応力解析を行い、非線形解析と応力分布を比較した。Step1 と Step2 の相当応力の分布を Fig.1.4.2-1 に、Step3 および 4 の相当応力の分布を Fig.1.4.2-2 にそれぞれ示す。
現状の解析では底面を拘束しているため底面に応力が集中する結果となっている。Fig. 1.4.2-1 および Fig.1.4.2-2 から弾塑性解析は弾性解析に比べて降伏応力に近い値で応力が制限されており、弾塑性解析として妥当な結果になっている。

22

Step1

弾塑性解析

弾性解析

Step2

節点相当
応力分布

弾塑性解析

弾性解析

節点相当
応力分布

Fig.1.4.2-1. Distribution of nodal equivalent stress at step1 and step2.
Step4
Step3

弾塑性解析

節点相当
応力分布

弾性解析

弾塑性解析

弾性解析

節点相当
応力分布

Fig.1.4.2-2. Distribution of nodal equivalent stress at step3 and step4.
文献[1] 宮本 学・尾上 俊雄・成田 貴一 「高炉用耐火物の高温における変形および破壊挙動」鉄と鋼第 67 年第 2 号,pp.313-322 ,1981[2] 谷本 浩一・高橋 由紀夫「熱ラチェット変形解析精度向上のための非弾性構成式の研究」電中
研報告：T93022, 1993
[3] 松井、巨、井上「超高温における SUS304 鋼の非弾性特性と統一型構成式」材料(J. Soc. Mat.
Sci. ,Japan),Vol.41,No.466,pp1153-1159,July 1992
(2) 洋上ウィンドファームのスーパーシミュレーションと V&V、デジタルツイン

23(2-1) Alpha Ventus の解析
スパコン版 RIAM-COMPACT（以下 RC-HPC とする）では、デカルト座標系不等間隔スタガード
格子に基づき、数値計算手法には（有限）差分法を用いる。LES の SGS モデルには、計算安定性に優れ、かつ壁面減衰関数を必要としない混合時間スケールモデルを採用する。陽的フィルタ操作にはシンプソン則を適用する。速度場と圧力場のカップリングアルゴリズムには、陽解法を基礎とした部分段階法を採用する。圧力に関するポアッソン方程式は SOR 法により緩和計算する。空間項の離散化に関して、対流項は補間法による 4 次精度中心差分を基礎とし、4 階微分の数値粘性項を付加する(修正 3 次精度風上差分)。ここで、数値粘性項の重みは 0.5 とし、その影響は十分に小さくした。一
般的な 3 次精度風上差分の河村-桑原スキームでは 3.0 である。残りの空間項は 2 次精度中心差分で近似した。
令和 3 年度は風車ウエイクの新しい可視化技術の検討を行った。本研究で提案する新しい風車ウ
エイク可視化手法のヒントを Fig.2.1-1 に示す。火災現象の移流・拡散現象と風車ウエイク現象の類
似性に着目し、スカラー濃度場を用いて風車ウエイク現象を可視化する手法を提案する。Fig.2.1-2 には、本提案手法を大規模洋上ウィンドファームへ適用した例を示す。ここで、カラーの部分が仮想的に配置したスカラー濃度の発生源である。一方、モノクロ分布は主流方向の風速分布である。カラーのスカラー濃度の空間的な拡がりから、風車ウエイクの影響範囲を視覚的に捉えられている。本提案
手法は、風車ウエイクが相互に干渉した複雑な状況において、注目する風車が形成するウエイクの影
響範囲を正確に把握するのに非常に有効であることが示された。

Fig.2.1-1. Hint for the new wind turbine wake visualization method proposed in this research.

Fig.2.1-2. An example of applying this proposed method to a large-scale offshore wind farm
(RC-HPC).
24

(2-2) NREL5MW 風車の後流解析と V&V、後流分布モデル構築(2-2-1) FFB による NREL5MW 風車の後流解析
NREL5MW（風車直径 D＝126m）下流の流れ場を LES 解析により求め、風車下流のデータを求めた。
「京」コンピュータを用いたプロジェクト（ポスト「京」重点課題 6、課題代表者 吉村忍、平成 26 年～令和 1 年）では、風車下流 1D までの距離で解析を行い、風車の出力特性を求めるのに十分な解像度の解析条件を求めたが、本解析では風車下流距離を 12D まで拡張した。一般的な洋上風車では風車間距離が 7D から 10D とすることが多く、本解析によりこれらの領域の風車後流のデータを詳細に求めた。風車直径が 126m と大きいため、解析領域の出口は 1,512m 下流となる。この領域までの流れ構造を正確にとらえるため、風車後流領域の格子サイズは 30cm 以下とした。すなわ
ち、下流領域を流れ方向に 5,040 分割して、詳細に解析を行った。格子点数は 8.3 億である。風車ブレード表面の境界層メッシュは無次元高さ y+=100,主流方向 z+=500,幅方向 x+=200 である。このメッシュは翼面上の境界層を解像するには十分ではないが、風車の出力特性を求めるには十分な格子解像度であることを確認しており、後流域のメッシュとスムーズに接続することにより、翼面で生じた数 mm 程度の渦が後流域で成長し、下流に進んで行く様子を再現することが可能である。ただし、風車のブレード周速度は、風車に流入する速度の 5 倍から 8 倍程度であり、周速度が大きいことから翼周りのレイノルズ数は 107 オーダーとなる。このように高いレイノルズ数の流れでは、翼前縁の流れを正確に解析しないと翼の揚抗比が正しく得られない。予備計算により風車 1 回転を 6 万ステ
ップ程度で解析すると出力係数が運動量理論及び実験で得られる値に近いことが確認されている。
このため、解析領域入口から出口まで主流が到達するまでに風車は 28 回転（168 万ステップ）の計算を行う必要がある。本解析では助走計算として、風車ブレードを回転させない静的オーバーセット
解析で 180 万ステップ、風車を回転させる動的オーバーセット解析を 35 回転行い、合計で 390 万ス
テップ（実解析時間で 5 分 24 秒間に相当）の解析を行った。
Fig.2.2-1 に風車ブレード中心高さの主流方向速度分布を示す。z/D=0 の位置に風車が設置されている。縦軸は主流速度で無次元化した主流方向速度である。風車上流から流入した流れは、風車の位置で急激に速度を落とし、風車を通過したのち増速していく。解析では風車下流の計算領域を 3D、6D とした結果も合わせて示す。風車下流 6D 付近から主流速度が徐々に回復し、8D から 9D 付近でほぼ主流速度と一致する。このことから風車間距離を 10D 程度とすることが望ましいという従来からの知見と一致する結果が得られた。洋上風車では設備を小さくし、単位面積当たりの発電量を向上させたいことから、2 台目の風車をどこに置くかが重要となる。この解析結果では 5D 付近では主流
速度の 70%程度であることから風車出力は本来得られる出力の 34%程度になると考えられる。7D から 8D とした場合は、50 から 70 程度と予想される。このことから風車設置位置としては少なくとも直径の 7 倍程度とすることが望ましいと考えられる。
Fig.2.2-2 および 3 に風車周りの流れ場のコンター図を示す。風車周りの流れを上面から見た速度分布では、風車後流が蛇行している様子が見える。この解析では水平方向の速度は一様であることから、風車自体の作る流れにより蛇行が生じることがわかる。また、上下方向の分布をみると後流はやや下向きに流れていることがわかる。

25

Fig.2.2-1. Velocity profiles of main flow velocity at wind turbine hub height, TSR = 7.

Horizontal section

Vertical section
Fig.2.2-2. Velocity pattern around a wind turbine, TSR = 7.
Fig.2.2-3 は、風車入口に大気乱流境界層を模擬した流れを導入した結果である。解像度の粗い解析では風車の後流速度欠損が急速に回復するが、本解析の結果では、風車後流の流れは非常に長い距離その構造を保っていることがわかる。Fig.2.2-4 に後流の時系列変化を示す。流れ場が時間ともに蛇行、上下動していることがわかる。
Fig.2.2-5 に主流乱れの影響を示す。グラフは赤線が一様乱流、青線が一様流を示す。速度の測定
位置は風車下流 1D の位置である。また出力特性のヒストリを示す。入口速度乱れは大気乱流に近い5%を導入しているが、乱れの影響は小さく、風車下流の速度欠損はほぼ同程度であった。これにともない風車の出力特性は乱れ強度、速度勾配にはあまり影響せずにほぼ一定であった。風洞実験など
で主流乱れを入れた実験を行うと風車の揚抗比が大きく変化するが、これは風洞実験で対象としている翼のレイノルズ数が小さいためと考えられる。レイノルズ数 105 程度の翼周りの流れでは主流乱れにより前縁にはく離泡ができるなどの大きな変化が観察されるが、本解析ではそのような影響は小さかった。逆に格子解像度が不十分で翼の転向を再現できないような解析では、実際の風車の揚抗比を得ることが難しかった。乱れの影響については渦スケールの影響も考えられるので、今後、渦ス
ケール、乱れ強度をパラメータして解析を行う必要がある。
Fig.2.2-6 に後流の速度分布を示す。後流分布はダブルガウシアン型と呼ばれる二つの正規分布関数を組み合わせた分布となっていることがわかる。二つの山は下流に行くにつれてなだらかになり、最終的には一つの山（ガウス分布）となる。Fig.2.2-7 に後流幅の下流方向変化を示す。一般に地面（水面）により後流幅に制限が加わる鉛直方向よりも水平方向分布の後流幅が大きくなることが知られており、解析結果も同様の傾向が得られた。
大規模乱流解析により後流の非等方性、蛇行などの特徴を捉えることができた。

26

Fig.2.2-3. Wake flow around a wind turbine, TSR = 7.

Fig.2.2-4. Temporal variation of wake distribution behind the wind turbine, TSR = 7.

－一様乱流
－一様流

27

Fig.2.2-5. Effects of inlet turbulence on wind turbine wake flow and power coefficient, TSR =7.

Fig.2.2-6. Horizontal distribution of mean flow velocity in the wake of wind turbine, TSR = 7.

28

Fig.2.2-7. Comparison of horizontal and vertical wake widths, TSR = 7.
（2-2-2）FFB による風車単機解析と RC-HPC との比較
ここでは、FFB による風車単機との比較を目的に行った計算結果について説明する（Fig.2.2-8 を参照）。令和 2 年度の FFB との定量的な比較（風車ハブ高さにおける平均風速分布の比較）と議論を踏まえて、今期はアクチュエータラインモデルにおける抵抗係数の与え方に改良を施した。本計算における格子分解能は各方向においてΔx=0.5m, Δy=Δz=2.0m の等間隔とした。各方向の格子数はNx=1261, Ny=251, Nz=231 であり, 総格子点数は約 7,300 万点である。無次元時間刻みは 5.0×10-4であり、無次元時間 50～100 において時間平均を施した。特に瞬間場からは、風車下流に形成されているウエイク領域がブレードからの翼先端渦、ナセル、タワーが形成する剥離流の影響を強く受けていることが見て取れる。

Fig.2.2-8. Numerical results for a single wind turbine (RC-HPC), upper is the instantaneous
field, lower is the time-averaged field.

29

(2-3) NREL5MW 風車のタンデム解析と V&V(2-3-1) FFB による NREL5MW 風車のタンデム解析
洋上風車では、上流に置かれた風車の後流の影響により下流の風車の性能が十分に得られないことが問題となる。そこで風車を 2 台配置した解析を行った。はじめに風車間距離を風車直径の 3 倍とした場合について解析を行った。通常、風車間距離を直径の 3 倍とすることは少ないが、まずは後流干渉の影響を確認するため、2 台の風車を一直線上に並べて解析を行った。
基本的な解析条件として NREL5MW 機の設計運転条件である周速比 TSR=7 を標準とし、主流乱れがない場合と主流に一様な乱れ（5%）を入れた場合について比較した。格子点数は約 6 億 5 千万である。計算に使用するノード数は「富岳」1,536 ノード、羽根車 1 回転あたりの計算時間は、約 8時間である。これより、羽根車 1 回転あたりに必要な計算リソースは約 1.2 万 NH/rev.であり、280
rev. （7 ケース×40 rev.）の計算に必要な計算リソースは 336 万 NH である。TSR は 5,6,7,8 の 4
ケースについて調べた。
Fig.2.3-1 に入口乱れなしの場合のスラスト係数、出力係数の時系列変化、Fig.2.3-2 に入口乱れありの場合を示す。主流変動が加わることによりスラスト係数、出力係数ともに時間的に変動するが、時間平均値は入口乱れにより大きく変化しないことがわかる。特に 1 台目の風車の平均値は主流の乱れにほとんど影響を受けなかった。この結果は 1 台の風車を解析した場合と同様である。
Fig.2.3-3 に周速比に対するスラスト係数、出力係数の依存性を示す。先に示したように 1 台目の風車は主流乱れの影響をほとんど受けない。一方、2 台目の風車は出力係数がゼロもしくはマイナスとなり流れから出力を取り出せない状態である。これは風車間の距離が 3D と近いためであり、この距離では 2 台目の風車から出力が得られない。主流に乱れがある場合、スラスト係数、出力係数ともにさらに小さくなる傾向がみられる。2 台目の風車は速度欠損の大きな流れが流入しているので、主流の影響は小さいと予想されたが、1 台目の風車よりも主流乱れの影響が大きくなっている。
Fig.2.3-4 にハブ高さにおける主流速度の流れ方向分布を示す。風車に流入する速度を風車の上流
0.5D の位置で定義した場合、1 台目の風車の流入速度はほぼ主流に等しいのに対して、2 台目の風車の流入速度は主流の 40%程度である。単純に風車出力が速度の 3 乗に比例するとすると、風車の出力係数は 1 台目の風車の 6%程度となる。1 台目の風車の出力係数は 42%程度であることから、2 台
目の風車の出力係数は最大でも 2.7%となる。実際に得られた出力係数はさらに低いことから流入速度の減速だけでなく、風車が作る乱流の影響も受けているものと考えられる。

30

Fig.2.3-1. Time variation of thrust and output coefficients (without inlet turbulence) .

31

Fig.2.3-2. Time variation of thrust and output power coefficients (with inlet turbulence) .

32

Fig.2.3-3. Comparison of thrust and output power coefficients of the first and second windturbines.

Fig.2.3-4. Streamwise velocity at the hub height along the axial direction.
Fig.2.3-5 と Fig.2.3-6 に翼面の圧力分布を示す。これまでの議論と同様に 1 台目の風車は主流乱れの影響をあまり受けておらず、翼周りの圧力分布より揚力が発生し、主流から動力を得ているが、2 台目の風車の翼面圧力分布は流体仕事をするような分布になっていないことがわかる。
Fig.2.3-7 と Fig.2.3-8 に風車周りの流れ場を示す。ハブ高さでは主流乱れの影響は小さく見えたが、風車の高さ方向の分布は主流乱れがある場合とない場合で異なることがわかる。この違いが出力係数、スラスト係数の違いを生じさせていると考えられる。

33

1st Trubine
2nd Turbine
Fig.2.3-5. Static pressure distribution on the blade surface (without inlet turbulence) .

1st Trubine
2nd Turbine
Fig.2.3-6. Static pressure distribution on the blade surface (with inlet turbulence).

34

without inlet turbulenc

with inlet turbulenceFig.2.3-7. Time averaged mainstream velocity distribution in the central cross-section, TSR=7.

with inlet turbulence

without inlet turbulenceFig.2.3-8. Instantaneous mainstream velocity distribution in the central cross-section, TSR=7.
風車間距離が近い場合は、下流側の風車から出力を取り出すことができないことがわかった。そこ
で風車間距離を実際のウィンドファームに近い直径の 7.5 倍とし、3 台の風車を並べた解析を行った。
解析規模はベースメッシュが 2.2 億であり、このメッシュを 1 回リファインすると前述の単体解析、タンデム解析と同じ解像度となる。3 台の風車を風車間距離 7.5D で並べた場合、解析空間は 18D となる。このため、この領域を主流が通過するには、風車を最低 50 回させる必要がある。そこで助走35

解析をベースメッシュで行い、流れを十分発達させてから、メッシュのリファインを行い解像度を上げて解析を行った。リファイン後の要素数は約 18 億である。
Fig.2-3-9 に解析空間の概要を示す。風車間距離は異なるがメッシュをリファインした場合の格子解像度がこれまでのものと同様である。
Fig.2.3-10 から 13 に平均速度場、瞬時速度場を示す。第 1 風車、第 2 風車の後流は 2 台目の風車、
3 台目の風車に到達する前にかなり回復していることがわかる。Fiig.2-3-14 から各風車の後流が相似であることがわかる。Fig.2-3-15 にハブ高さにおける主流速度を示す。各風車への流入速度はそれ
ぞれ主流速度の 1.0,0.7,1.0 倍である。このことから出力係数は 0.42, 0.14, 0.42 程度と予想される。
Fig.2-3-16、17 に出力係数を求めた結果を示す。第 1 風車の出力係数は 0.42 であるのに対して、
第 2 風車は 0.35、第 3 風車は 0.39 程度であった。主流の平均速度からすると第 2 風車は予想よりも出力が大きく、第 3 風車は低くなっている。この原因はハブ高さの速度分布だけでは出力係数を算出するには不十分であること、第 3 風車には非常に強い乱れが流入していることも考えられる。また、
この解析は解析時間がまだ不十分であり、出力係数の平均値が時間とともに変化していることからさらに計算を続け、出力特性を評価する必要がある。風車間距離が大きくなると解析規模、解析時間が膨大になり、従来のコンピュータでは LES 解析を行うことが難しかったが、「富岳」をもちいることによりこれらの解析が可能となった。この結果を踏まえ、次年度は計画通り、より大規模なウィン
ドファームを模擬した解析、周速比を変えた解析を行い、風車後流の特性を数値解析より理解すると共に、ウィンドファーム流れの数値データベースを構築し、サロゲートモデルの開発を進めていく。

Fig.2.3-9. Analytical model of three wind turbines in series.

Fig.2.3-10. Average velocity field in wind farms, TSR = 7 (Horizontal section).

36

Fig.2.3-11. Average velocity field in wind farms, TSR = 7 (Vertical section).

Fig.2.3-12. Instantaneous velocity field in wind farms, TSR = 7 (Horizontal section).

Fig.2.3-13. Instantaneous velocity field in wind farms, TSR = 7 (Vertical section).

37

Fig.2.3-14. Bird's-eye view of the flow around the wind farm and the pressure on the water
surface, TSR = 7.

Fig.2.3-15. Velocity profiles of main flow velocity at wind turbine hub height, TSR = 7.

Fig.2.3-16. Comparison of time history of instantaneous output power coefficient of threeturbines, TSR = 7.

38

Fig.2.3-17. Comparison of time averaged output power coefficient of three turbines, TSR = 7.
（2-3-2）RC-HPC による計算結果と FFB による計算結果との比較
Fig.2.3-18 には、瞬間場（無次元時間 50）に関してスパン中央断面における流れの可視化（sideview）を示す。下流側に設置された風車は、上流側に設置された風車のブレード、ナセル、タワーが形成するウエイクの影響を強く受けていることが見て取れる。Fig.2.3-19 には、風車ハブ高さにおける平均風速の比較結果を示す。図中には、RC-HPC の令和 2 年度の結果（Before improvement）と令

和 3 年度の結果（After improvement）の両者を示している。令和 3 年度はアクチュエータライン
モデルにおける抵抗係数の与え方に改良を施した結果、FFB の結果と良好な一致を得ることに成功した。

39

Fig.2.3-18. Numerical results for two wind turbines (RC-HPC), upper is the instantaneous
field, lower is the time-averaged field.

Fig.2.3-19. Comparison of time-averaged wind speeds at wind turbine hub height(measurement line), as shown in Fig. 2.3-18.

(2-4) NREL5MW 風車の累積疲労損傷解析「富岳」上での累積疲労損傷解析環境を整備し、タンデム風車モデルを使って動作試験を行った。
解析環境整備として、これまでワークステーションで実施していた累積疲労損傷解析ソルバー
ADVENTURE_Fatigue を「富岳」に移植した。これにより、流体-構造-疲労損傷の全ての解析を「富岳」で連続して実行する環境を実現した。
動作試験として、2 基の発電用風車（設置間隔=3D=378m）のモデルを作成し、TSR=7、境界層なし、乱れ成分なしの流入条件で流体解析を行い、その結果から後方風車のブレード表面における流体
力履歴を抽出し、構造解析及びその結果を用いた累積疲労損傷解析を実施した。流体力履歴の抽出とデータ形式の変換には REVOCAP_Coupler とファイルフォーマット変換のためのユーティリティプ
ログラムを使用した。Fig. 2.4-1 に構造解析で得られた解析開始から 5.0 秒後の変形（変位を 10 倍に強調した）を可視化したものを、また Fig.2.4-2 に 20 年連続運転した状態における最大疲労の発生個所を可視化したものを示す。これらの成果を踏まえて、流体側の流入条件を変えた 8 ケースについて流体解析の結果を用いた構造解析と累積疲労損傷解析を実施中である。

40

ま た 、 ADVENTURE_Fatigue に よ る ポ ス ト 処 理 と し て の 疲 労 損 傷 評 価 と 平 行 し 、ADVENTURE_Solid によるシミュレーションの中で疲労損傷による剛性低下を考慮可能な損傷モデル（progressive damage model）の検討を進めている。両者の比較によって、多面的な疲労損傷評
価が可能になる。progressive damage model の ADVENTURE_Solid への実装は完了した。引き続き、検証と妥当性確認を実施中である。

Fig.2.4-1 Deformation of wind blade at t=5.0s.

41

Fig. 2.4-2 Maximum fatigue location after 20 years of continuous operation.
(2-5) NREL15MW 風車の解析モデル構築（2-5-1）概要
今まで 5MW 風車で計算しているが、今後国内でも建設が予想される 15MW 風車について、現時点では実物がなく、シミュレーションによって評価するためには新たなメッシュが必要である。そこ
で、米国再生可能エネルギー研究所（NREL）が公開している情報をもとに、風車単機を含む風況解
析のための並列 LES 解析コード FFB 用メッシュを作成した。風況解析メッシュの中に、回転する風車ブレードを考慮したオーバーセットメッシュを組み込むこととし、メッシュを可視化したうえで品質の評価を行った。
（2-5-2）メッシュ生成作業の特長
メッシュを作成した NREL15MW 風車のモデルを Fig.2.5.1-1 に示す。メッシュ作成作業にあたり、ブレードの角度を地面に対し垂直にし、ブレードのピッチ角も一定とした。

42

Fig.2.5-1 NREL 15MW.
NREL 5MW の風車では境界層摩擦速度をもとに無次元化した場合、壁座標 (1000, 200, 2000)としている。境界層厚さ方向を y+=200 に設定した場合、出力係数を求めることができる解像度である。新しい 15MW でもこの解像度を維持できる翼弦（中央付近）の境界層厚さを推定すると、次のとおりになる。
・レイノルズ数 Re

2*107・境界層厚さδ

0.038 m

・摩擦係数

Cf

0.002

・壁面せん断応力 τ
・摩擦速度

U*=

12 Pa (動圧の 0.2%)
3.16 m/s(主流速度の 3%)

ここで、0.95mm が y+=200 に相当する。このことから、翼面の解像度は、翼スパン方向について
平均は 113.92mm、境界層厚さ方向 1mm、流れ方向 10mm のメッシュがあれば出力を求められると想定した。
この解像度でブレードのメッシュを作成すると 30 億程度になるため、FFB のリファイン機能を用いる場合は、2 回のリファインが必要であると考えられる。
ブレード解像度としては、境界層厚さ方向最小厚さ 1mm、境界層厚さ 40mm 程度を解くことを目43

指したメッシュが必要となる。そこで、スパン方向（風車直径）のメッシュ幅平均は 63.21mm 程度、
流れ方向 10mm 程度の解像度をリファインで求める必要があり、今回は 2 回リファインに対応するようなメッシュを作成した。
また、解析領域は次のように設定した。
・x:-1D～3D (全長 4D)・y:-3D～3D（全幅 6D)
・Z:-0.75D～2.25D (全高 3D)
形状表面のメッシュは CAD サーフェスへの投影処理を行っているため、CAD 形状に忠実なメッシュが生成される。このことにより、高精度な計算が可能になる。
（2-5-3）メッシュの可視化および品質チェック
作成したメッシュ図を Fig.2.5-2 から Fig.2.5-4 に示す。

44

Fig.2.5-2. nrel15mw_blade_triple.

Fig.2.5-3. nrel15mw_inner.

45

Fig.2.5-4. nrel15mw_outer.
これらのメッシュについて、Jacobian 値、セルの最大角、体積比により評価を行ったところ、おおむね良好な結果となり、品質の高いメッシュであることが確認できた。このメッシュが完成したことにより、NREL15MW 風車の流体解析を行う準備が整った。
(3) アプリケーションの「富岳」における性能評価とチューニング
(3-1) FFR-Comb
ベンチスケール石炭ガス化炉をラボスケール炉と同程度の格子密度で作成すると、計算格子は 32
億要素規模、領域分割数は 6 万 5 千程度になる。このベンチスケール炉に対して FFR-Comb で燃焼解析を実施するため、FFR-Comb の前処理を実施する必要がある。FFR-Comb の前処理プログラム
は 2 段階あり、1 段階目を pre_metis、2 段階目を prefflow と呼ぶ。
1 段階目の pre_metis ではメインソルバーの MPI 並列計算のための領域分割と通信テーブルの作成などを行う。pre_metis について、処理別にまとめて分割して予備計算をすると、32 億要素格子全ての処理に必要な時間に 37 日以上、実行に必要なメモリ量は 3.5TB と見積もられた。この
pre_metis について、プログラムのアルゴリズムの見直し、配列 allocation の最適化、openMP によるスレッド並列化などコードのチューニングを実施し、プログラムの高速化と使用メモリ量の削減を図った。チューニング後は pre_metis の 32 億要素格子に対する処理時間は計 46 時間まで短縮することができた。ただし「富岳」のプリポストノードの大容量キューは最大 24 時間のため、九州大学の大容量フロントエンドで実施した。

46

2 段階目の prefflow では領域分割部分ごとにコントロールボリュームの体積・面積・面ベクトルなどのジオメトリ情報や壁面距離などを計算する。このプログラムは 1CPU 計算で逐次処理で実装されている。
こちらも予備計算を行い、32 億要素格子の処理が 1 時間当たり 100 領域程度であった。
そこから 6 万 5 千領域全てを処理するのに 650 時間以上必要であると見積もられた。prefflow は領
域ごとのループ構造で処理が独立している部分がほとんどであるため、ループを分割し、分割したループごとに別プログラムとして同時並行して実行できるように高速化改修を実施した。高速化改修
後は同時に 10 個まで並行計算が可能となり、必要な処理時間を 180 時間まで短縮することができた。
FFR-Comb はメインソルバーで計算した後、後処理プログラムを使用して計算結果を可視化ソフ
ト用データに変換する。この後処理プログラムについても 32 億要素規模のデータに対応するため、任意断面を抽出しデータ量削減を行う機能を実装した。
最後にチューニング前後で FFR-Comb の処理に問題がないことを確認した。Fig.2.2-1 にテスト格子を使用したチューニング実施前後の FFR-Comb の解析結果を示す。上がチューニング前、下がチューニング後の図、左が圧力、右が流速 u 成分のコンター図である。チューニング前後で同じ結果が得られることが確認できた。

Fig. 2.2-1. Comparison of test results before and after FFR-Comb tuning.
(3-2) FFB
FrontFlow/blue（FFB）は Large Eddy Simulation（LES）をベースとする並列有限要素法流れ解析コードである。FFB の最大の特長は、乱流中の微小な渦のダイナミクスを直接計算することにより乱流を高精度に計算できることにある。
大型発電用風車周りの流れにおける微小な渦のダイナミクスを計算するために、高解像度の計算
格子が必要となり、「富岳」を用いて大規模な計算を短時間に実行することが求められる。この課題
では、「「富岳」成果創出加速プログラム「「富岳」を利用した革新的流体性能予測技術の研究開発」」
で開発されている、Large Eddy Simulation（LES）をベースとする並列有限要素法流れ解析コード
FrontFlow/blue(FFB)を用いて計算を行う。FFB の最大の特長は、乱流中の微小な渦のダイナミク
スを直接計算することにより乱流を高精度に計算できることにある。
「「富岳」成果創出加速プログラム「「富岳」を利用した革新的流体性能予測技術の研究開発」」
（課
題代表者

東京大学

加藤千幸）では、令和 2 年度までの FFB の開発において、単体性能

179GFLOPS（ピーク性能比 5.3%）を達成するとともに、直方体領域を直列に連結し、両隣の領域とのみ通信する通信コストが低い weak scale ベンチマークテストでは、「富岳」のほぼ全系を用いて85%の並列効率を達成している[1]。しかしながら、実用計算では並列効率が低下する課題が確認されている（Fig.3.2-1）。FFB は、独立した領域関連でデータを受け渡すオーバーセット計算機能をサポートする[2]。オーバーセット計算では、通常の領域境界に関する通信（以外、隣接通信）の他にオ
ーバーセット処理のための通信（以外、オーバーセット通信）が必要になり、これに起因して並列性能がさらに低下する課題が確認されている。そこで、令和 3 年度、上記の隣接通信およびオーバーセット通信の通信コストについて分析し、並列性能向上のための指針を得ている。本節では、その成果の概要を示す。

47

Fig.3.2-1. Weak scale benchmark test with 2 million grids per a computational node of
Fugaku.
「富岳」の計算ノードは 4 個の Core Memory Group (CMG)から構成され、Tofu は 12 個の計算ノードを１グループとして構成される。本報告書では、１グループを Tofu 装置と呼ぶ。Tofu 装置は格子状に配置されており、遠くの Tofu 装置と通信するには、間にある Tofu 装置を介して通信が行われる。したがって、N 個隣の Tofu 装置間の通信コストは、隣接する Tofu 装置間の通信コストのN 倍となる。Fig.3.2-2 に Tofu 装置、（Tofu 装置内に 12 個ある）ノードおよび（ノード内に 4 個ある）CMG（Core Memory Group）の関係を示す。令和 3 年度は、Fig.3.2-2 に示すように、Tofu シ
ステムにおける通信経路を考慮した通信コストと通信性能との関係を調査した。

Fig.3.2-2. Relationship of Tofu-unit, computational node and core memory groups (CMG) in
tofu system.
（1）隣接通信20 億グリッドの計算格子、「富岳」1000 ノードを用いた、Cavity flow 計算において、Tofu システ
ムのリンク（Tofu 装置を結ぶ線）毎に 1 回の通信における通信量を算出した。通信量が最も多くなるリンク（以下、最混雑リンク）の通信量と通信時間を、領域分割方法を変更して調査したところ、最混雑リンクの通信量が多いほど、通信時間が長くなり、通信性能が低下することが確認された。通常の METIS による分割をした場合の最混雑リンクの通信量は、計算領域をブロック状に分割した場合と比較し、10 倍程度多く、この結果、通信に要する時間は 4 倍程度になり、コード全体の実効性能が 1 割程度低いことが確認された。このことから、Tofu 空間において、各通信の Tofu 距離を減らし、最混雑リンクの通信量を減らすことが、高い並列性能を達成するために重要であることがわかっ

48

た。
（2）オーバーセット通信
風車まわり流れを対象に、6.3 億グリッドの計算格子、「富岳」1,536 ノードを用いた計算において、隣接通信およびオーバーセット通信における最混雑リンクの通信量を評価した。隣接通信のみの分
析では、最混雑リンクの通信量の増大とともに通信時間が増大したが、オーバーセット通信は、通信量自体は隣接通信よりも少ないにもかかわらず、通信時間が長くなっていることがわかった。両通信の通信パターンを比較したところ、オーバーセット通信は、大半の計算ノードでは近い Tofu 装置との通信しかしないが、特定の計算ノード間で遠い Tofu 間の通信をしていることがわかり、これが並列性能の低下の要因になっていることが推察された。この分析結果により、一部のプロセスにのみ遠くの Tofu 装置との通信が集中しないよう、オーバーセット通信にかかわるプロセスを Tofu 空間において近くに配置することにより、並列性能向上の向上が期待できる。
文献[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
method with application to turbomachinery and aeroacoustics, Transactions of ASME, Journal ofApplied Mechanics, Vol. 70 (2003), pp. 32-43.
(3-3) ADVENTURE_Solid2（3-3-1）事前検討
令和 3 年度は ADVENTURE_Solid （以下 ADV_Solid）の「富岳」向けチューニングの実施に向け、まず事前検討として ADV_Solid プロトタイプ実装を用いて各アルゴリズム（BDD、Local Solver
等）の計測および最適化を実施した。その結果、BDD に関しては線形代数ソルバーの OpenMP 対応化を進め、また、コース行列のスカイライン記憶のブロック化と BLAS、 LAPACK ライブラリの効率的な利用等により、億レベル自由度の大規模板曲げ問題においてピーク性能比 20%程度を実現した。また、Local Solver に関しては、BDD と同様にこれまでの 1 領域 1 コアの割当てから、1 領域
複数コアに変更し、それに伴い並列スカインソルバーを導入することにより、同じく億レベル自由度の大規模板曲げ問題においてピーク性能比 15%を実現した。これらの成果を ADVENTURE_Solid2に移植し、「富岳」上での最適化を進めることとした。
（3-3-2）計算性能および精度検証
大規模薄肉構造物モデル（一次要素を用いた 1F1 モデル(C1)）について「富岳」上の最適な並列パラメータを調査するため、ノード数、MUMPS プロセス数、パート当たりのサブドメイン数を振って各ケースの時間を計測した。MUMPS プロセス数=60 として ADVENTURE_Solid2 の 1CG 反復における各処理の平均処理時間とサブドメインに含まれる平均節点数の関係を Fig. 3.3-1 に示す。
これより、サブドメインあたりの節点数が小さい時にはコースグリッドの自由度数が増大するため
BDD-DIAG 前処理の求解処理に時間がかかる反面パートの内部求解処理は高速に処理できていることがわかる。サブドメイン節点数が大きい時にはコースグリッドの自由度数が小さくなるため BDDDIAG 前処理の求解処理は高速化されるが、パート内部の求解処理の処理時間が増大する結果となった。

49

各CG反復中の処理時間とサブドメイン節点数の関係(MUMPS=60)
0.4
0.35
平均BDD求解処理時間[s]

0.3
0.25

平均
(Precondition a
resudial vector)

0.2
0.15

平均Partsolve

0.1
0.05
0
0

200

400
600
サブドメインあたりの節点数

800

1000

Fig. 3.3-1. Comparison of each CG process elapse time and average number of the nodes per
subdomain (MUMPS process=60).
また ADVENTURE_Solid2 の C1 モデルの静解析全体の処理時間とサブドメインに含まれる平均
節点数の関係を Fig. 3.3-2 に示す。静解析全体の Elapse Time は変動するものの、サブドメインあたりの節点数=400 付近までおおよそ一定になっていることがわかる。

静解析のElapseTimeとサブドメイン節点数の関係(MUMPS=60)
6000

ElapseTime[s]

5000
4000

3000
2000

1000
0
0

200

400
600
800
サブドメインあたりの節点数

1000

1200

Fig. 3.3-2. Comparison of static analysis elapse time and average number of the nodes per
subdomain (MUMPS process=60).
収束判定基準=1.0×10-4 とした C1 モデルの収束までに行った CG 反復数とサブドメインあたりの平均節点数の関係を Fig. 3.3-3 に示す。Fig. 3.3-3 からサブドメインあたりの平均節点数が大きくな
ると CG 反復回数も増大していることがわかり、これにより Fig. 3.3-2 における節点数=400 付近のElapse Time がほぼ一定になっている原因になっていると考えられる。これは 1 サブドメインが担当する領域が大きくなると BDD-DIAG 前処理の効果が小さくなるためであると考えられる。

50

静解析の総CG反復回数とサブドメイン節点数の関係(MUMPS=60)
14000

CG反復回数

12000
10000
8000
6000

4000
2000
0
0

200

400
600
800
サブドメインあたりの節点数

1000

1200

Fig. 3.3-3. Comparison of CG iteration number and average number of the nodes per
subdomain (MUMPS process=60).
1F1 モデルについて一部の損傷を想定し剛性が変化した場合を考慮した解析を「富岳」上で行った。さらに塑性変形に早い段階で到達するように加速度振幅を変化させている。ケース表を Table
3.3-1 に示す。ここで、データを取得した節点 ID とコンポーネントの各部位を Table 3.3-2 に示す。
Table 3.3-1. Cases performed (C1 / C2 model).
ケース
C1m0a0C1m1a0C2m0a0C2m1a0C2m1a1C2m0a1モデル
四面体 1 次
(C1)

物性値

加速度

「京」と同一オリジナル

四面体 1 次コンクリートヤング率×0.7 倍(C1)（損傷を想定）四面体 2 次
(C2)

「京」と同一オリジナル

四面体 2 次コンクリートヤング率×0.7 倍(C2)（損傷を想定）四面体 2 次コンクリートヤング率×0.7 倍(C2)（損傷を想定）四面体 2 次
(C2)

オリジナル

オリジナル

加速度 1.2 倍

加速度 1.2 倍「京」と同一Table 3.3-2. Node ID of data acquisition points.
節点 ID
コンポーネント部位
3899

Mass 37 炉容器の上部

53887

Mass 37 炉容器の上部

13946255

Mass 5 建屋のスタビライザの上付近
51

C1 モデルについて Fig. 3.3-4 に節点 ID=53,887(Mass 37)の変位差を計算したものを示す。これ
より、変位差は時間的に振動しており一部のヤング率が変化したことによる固有振動数の差が見られることがわかる。
節点ID=53887
差異(京と同一-原子炉建屋
×0.7倍)
節点 ID=53887
差異（「京」と同一－建3.00E+00

屋×0.7 倍）2.00E+00

差異[mm]

1.00E+00
変位x

0.00E+00

変位y
変位z-1.00E+00-2.00E+00-3.00E+00
0

50

100

150

200

250

300

350

400

450

ステップ

Fig. 3.3-4. Displacement difference of C1m0a0-C1m1a0 case at node ID=53887 (C1 model).
四面体 2 次要素(C2)モデルについて（オリジナルケース: C2m0a0）－（コンクリートヤング率×0.7 倍: C2m1a0）として差異を計算したものを Fig. 3.3-5 に示す。C1 モデル同様に固有振動数が変わったことによる時間的な振動が見られている。さらに、節点 ID=13,946,255 についてはヤング率が低下した際に、-Z 方向に沈み込むような変位となった。
C2モデル 節点差異(オリジナル物性 - コンクリ×0.7)
3.00E+00
2.00E+00

差異[mm]

1.00E+00

3899差異x
3899差異y0.00E+00

3899差異z13946255差異x-1.00E+00

13946255差異y13946255差異z-2.00E+00-3.00E+00
0

100

200

300

400

500

600

700

800

900

ステップ

Fig. 3.3-5. Displacement difference of C2m0a0-C2m1a0 case at node ID=3,899 and 13,946,255(C2 model).
C2 モデルについて加速度振幅を 1.2 倍としてコンクリートのヤング率を（オリジナル: C2m0a1）－（0.7 倍: C2m1a1）で変位差を計算した結果を Fig. 3.3-6 に示す。Fig. 3.3-6 より、加速度設定がオリジナルのケースと同様に固有振動数の差により各節点の変位差が時間的に振動していることが52

確認できた。Fig. 3.3-6 における変位の差異が 200~300 ステップ付近で Fig. 3.3-5 から比較すると入力加速度同様に約 1.2 倍になっていることがわかる。
以上のような ADVENTURE_Solid2 の計算性能および精度検証を目的とした大規模薄肉構造物モ
デルに対する解析結果から、現状のコードでは、コースグリッド修正に対する計算コストが高く、それが全体の計算性能を低くしていることがわかった。そこで、部分領域数を減らしコースグリッド問
題の規模を小さくすることを考えた。その結果、部分領域問題の規模が大きくなるため、これについての線形問題の求解にも MUMPS を使えるように解析コードを改良した。プロトタイプコードは完成し、C1 および C2 モデルの解析が可能であることを確認した。
C2モデル加速度 1.2倍 節点差異(オリジナル - コンクリ×0.7倍)
4.00E+00
3.00E+00
2.00E+00

差異[mm]

1.00E+00

3899差異x
3899差異y0.00E+00

3899差異z13946255差異x-1.00E+00

13946255差異y-2.00E+00

13946255差異z-3.00E+00-4.00E+00
0

100

200

300

400

500

600

700

800

900

ステップ

Fig. 3.3-6. Displacement difference of C2m0a1-C2m1a1 case at node ID=3899 and 13946255(C2 model).
(3-4) ADVENTURE_Thermal
ADVENTURE_Thermal（以下 ADV_Thermal）の「富岳」向けチューニングの実施として、各アルゴリズム（BDD、 Local Solver 等）の計測を行った。ADV_Thermal の公開版に対して OpenMPを導入し、ハイブリッド並列化し、いくつかの簡単な数千万要素規模のテスト問題を用いて動作検証を実施した。その結果、ADV_Thermal に OpenMP が正しく動作し OpenMP 導入の導入が完了した。ただし、BDD を実行した際に通信エラーが発生するバグを取り除くことができなかったため、引続きこの問題に対処することとした。また検討の結果、DDM アルゴリズムのボトルネックとなる箇所の性能向上が必要であり、ADVENTURE_Solid での成果を有効活用しつつプログラムをさらに改善することを確認した。
(3-5) REVOCAP_Coupler
令和 2 年度に「富岳」に移植したコードをもとに、「富岳」の MPI 仕様及びファイルシステムの制約に対するチューニングを行い、ラボスケール石炭ガス化炉の燃焼流－伝熱連成解析におけるオン
ライン連成機能、及び洋上ウィンドファームのスーパーシミュレーションにおけるオフライン連成機能の両者ともに良好に動作することを確認した。
(3-6) RC HPC 版
RIAM-Compact（RC） HPC 版は、将来の応用展開を見越し、「富岳」を始めとした様々な計算機53において高速に動作するプログラム開発を進めている。令和 2 年度は、低 B/F のアーキテクチャの特徴をもつ計算機において、非圧縮流解析において高速計算の核となる圧力のポアソン方程式を効率的に解く手法として、Parallel Cyclic Reduction（PCR）法を提案・実装し、その効果を確認してきた。令和 3 年度は、MPI 並列化（FLAT MPI）の作業を実施し、基本的な性能の確認を強スケー
リング評価により実施した。Table 3.6-1 に性能測定の条件を、Fig.3.6-1 に解析したモデルを示す。
Table 3.6-1. Conditions for performance measurement.
Calculation model
Number of grid: nx=605 ny=626 nz=231
Number of windmill: 12
Computer
FUGAKU (A64FX)Compile option-Cpp -Kfast -Kpreex -Kocl -DVALIDATION -DTIMER

Fig.3.6-1. Calculation model (12 windmills).
Table 3.6-2 に測定した並列パターンを示す。Node あたり 48 コアであるので、48 並列以降は
48×Node 数での並列数としている。並列時の領域分割方法として、主流方向 X は分割せず、主流と垂直な YZ 面で領域分割を行う。これは、SIMD 演算器を有効活用するため、メモリ連続な X 方向のメモリアクセスを効率化するためである。分割数 Y と分割数 Z は、通信量を削減するため、ガイドセルの体積が小さくなる（各ランクの担当領域のセル数 Y とセル数 Z が近い値になる）ように決定している。ny > nz であるため、分割数 Y≧分割数 Z となる。

Case
A
B
C
D
E

Table 3.6-2. Measurement cases.
Division# of div.
Y
Z
1
1
1
2
2
1
4
4
1
8
4
1
12
6
1

54

# of node
1
1
1
1
1

F
G
H
I
J
K
L
M

24
48
96
192
384
768
1536
3072

8
12
16
24
32
48
64
96

1
1
2
4
8
16
32
64

1
1
2
4
8
16
32
64

計算時間の測定結果をノード内とノード間に分けて Fig.3.6-2 と Fig.3.6-3 に示す。768 並列以降、台数効果が小さくなっている（並列化効率が下がっている）が、これはガイドセル巾２に対して 768～3,072 並列では、各ランクの担当領域サイズが小さすぎるためである。台数効果を 1~3,072 並列で示すと低並列時の台数効果が分かりにくいため、1~384 並列と 1~3,072 並列の２種類で Fig.3.6-4 とFig.3.6-5 に示す。なお、メモリ使用量については、1～48 並列は 1 ノードにて実行しているため、ほぼ同程度の使用量であったが、96 並列以降は複数ノードでの実行となるため、ノード数に応じてメモリ使用量は減少していることを付記する。

計算時間（ノード間:48-3072並列）
計算時間（ノード間：48-3072 並列）

時間（Sec）
時間 (sec)

時間（Sec）
時間 (sec)

計算時間（ノード内:1-48並列）
計算時間（ノード内：1-48
並列）
10000
9000
8000
7000
6000
5000
4000
3000
2000
1000
0

300

250
200
150
100

50
0

0

10

20

30

40

50

0

60

500

1000

1500

Fig. 3.6-2. Computation time of inside a
node.

定値）

台数効果(理想値）

1000
1000
500500
0

0

3

0

500

500

1000 1500 15002000 2000
2500
5
0

1000

並列数

2500
3000

3000
3500
0

450
450
400
400
350
350
300
300
250
250
200
200
150
150
100
100
5050
0 0 0

4

並列数
Fig. 3.6-4. Scalability (1~384 parallel).

台数効果(測定値）

台数効果（測定値）

台数効果(理想値）

台数効果

台数効果

台数効果（測台数効果

台数効果

台数効果(測定値）

1500
1500

0

3500

並列台数効果(1-384並列）
並列台数効果（1-384 並列）

並列台数効果(1-3072並列）
並列台数効果（1-3072
並列）

2000
2000

3000

Fig. 3.6-3. Computation time between
nodes.

3500
3500

2500
2500

2500

並列数

並列数
並列数

3000
3000

2000

並列数

0

5050

100
100

150 150200

200250

300
250

350
300

400
350

並列数
5

並列数

450
400
0

Fig. 3.6-5. Scalability (1~3072 parallel).

スケーラビリティが低い理由として、今回のモデルは Fig.3.6-1 に示すように 12 台の風車周りの流れを計算するベンチマークとなっているが、12 台の風車の境界条件の実装が効率的でないためで

55

あることがわかっており、この点は今後改修する予定である。
令和 3 年度のまとめとして、・YZ の２次元分割化、メモリ削減および制御部の C++化を行ったフラット MPI のコードに対して、「富岳」において 1～3,072 並列での測定を行った。
・性能測定結果によると、384 並列までは十分な並列性能を得られたと判断できるが、それを上回る並列数に対しては通信のコストが大きく、性能の低下が見られた。今回の評価は強スケーリング評
価であるため、並列数が増加すると相対的に部分領域が小さくなり、十分な性能が得られない傾向となる。具体的には、高並列では、各ランクが担当するグリッド数が 通信の袖巾に対して大きくないため（例：1,536 並列で担当グリッド数 Y 11 個 Z 11 個 通信巾 2)、通信コストの割合が増大したのが原因である。
スレッド並列化を行えば、ランクあたりの担当領域サイズを大きくできるため、通信コストの増大をある程度迎えることが期待できることがわかった。
(4) シミュレーションと AI（機械学習）の連携活用とデジタルツイン(4-1) ADVENTURE と AI 連携活用
ADVENTURE_Thermal（以下、ADV Thermal）が生成する大規模数値計算結果を用いて AI（以下、深層学習）により予測器を生成[1,2]する。予測器は、以下の目的のために利用する。
（目的１）大規模数値シミュレーションを行わないでも傾向を予測（目的２）超高速予測により解の存在パラメータ範囲の推定
つまり、大規模数値シミュレーションを実施する前の予備検討を行う簡易シミュレーター構築することが目標となる。数値シミュレーションの結果を用いて予測器を生成するためにはこれまでのデータサイエンス分野で前提[3]としている、データ数が一様に膨大にある、予測対象が事前にモデル化可能である、といった仮説（前提）は一般論として成立ない。物理現象を予測するためのデータは次の事項を前提とするのが工学的な観点で合理性があり本研究の調査対象でもある。
（仮説１）計算結果は数十から数千程度しか存在しない（極めて限られた学習データ）（仮説２）予測範囲は事前に決定可能
以下、ADV Thermal との連携を考慮した発熱体の位置条件（境界条件の変化）により測定点（固定）位置の予測が可能かどうかを調べた。この調査には次のような点に留意した。
（留意点１）仮説１から過学習を起こさないための正則化技術[4]をどこまで適用するか調査（留意点２）パラメータ空間の設定と精度の悪い箇所の関係をある程度理解[5]（留意点３）出現頻度が低いデータは予測精度が悪い（AI および機械学習の特徴）[6]（留意点４）１水準のデータでは学習ごとに予測結果が激しく変化（線形では２水準、非線形では最低３水準必要）[7]
これらに基づき、令和 2 年度成果報告書に記載の通り温度場の機械学習を実施した。令和 3 年度
は誤差原因を調査し高精度化に取り組んだ。
(4-1-1) 誤差原因の分析
令和 2 年度まで実施した手法を用いて精度向上の評価を行うため、左右対称位置に評価点を設定
し、中心から最も遠い評価点は補外となるように設定した。評価点の位置を Fig. 4.1-1 に示す。評価点は valid●で示される箇所に配置した。全体の分布を表現するための固定された位置の学習用デー
タは、train▲で示される。この配置におけるシミュレーションによる温度分布の１例を Fig. 4.1-2 に示す。実際の学習データの生成は発熱体を移動させて生成する。
⑦
56

①
②

⑥

③

⑤
④

Fig. 4.1-1. Locations of training and validation data.

Fig. 4.1-2. Training example of temperature distribution generated by numerical simulation.
予測結果を Fig. 4.1-3 に示す。また、誤差の評価結果を Table 4.1-1 および 4.1-2 に示す。Table
4.1-1 より全体の約 91%が相対誤差±4%未満となり平均誤差も約 1.6%、絶対誤差では全体の 92.7%が±2℃以内になる。また、決定係数も 0.968 と良好な結果となることから適切な学習が行われていると考える。一方で、相対誤差 10%を超える予測結果が 455 個発生している。特に最大相対誤差は15.2%となっており、現状の精度では代替モデルとしての精度を満足しない。また、Fig. 4.1-10 より高温度域で大きなプラス誤差が確認できる。これらの、誤差の原因を解明すべく t-SNE[1-8]を用いて分析を行う。
Table 4.1-1. Summary of numbers of data sets corresponding with relative error ranges.

57

誤差範囲 [±%]
± 0%～± 2%

個数 [個]
33578

総数 [個]33578

割合 [％]
75.068最大誤差 [%]± 2%～± 4%

7560

41138

16.90115.291± 4%～± 6%

2327

43465

5.202

最小誤差 [%]± 6%～± 8%

690

44155

1.543

0.0

± 8%～± 10%

120

44275

0.268

平均誤差 [%]± 10%～

455

44730

1.017

1.573

Table 4.1-2. Summary of numbers of data sets corresponding with absolute error ranges.

誤差範囲 [±℃]

個数 [個] 総数 [個] 割合 [％]± 0.0℃～± 1.0℃

37000

37000

82.719

最大絶対誤差± 1.0℃～± 2.0℃

4804

41804

10.74

9.550 [℃ ] / 14.692 [％ ]± 2.0℃～± 3.0℃

1512

43316

3.38

最小絶対誤差± 3.0℃～± 4.0℃

722

44038

1.614

0.000 [℃ ] / 0.000 [％ ]± 4.0℃～± 5.0℃

254

44292

0.568

平均絶対誤差± 5.0℃ ～

438

44730

0.979

0.662 [℃ ] / 1.019 [％ ]Fig. 4.1-3. Comparison with reference and prediction temperature for all of validation data
sets: when dots places on red line, prediction would be correct.
誤差の大きさをラベルとして色分け指定をした t-SNE[8]による可視化結果を Fig. 4.1-4 に示す。
t-SNE は多次元量の相対距離を保持しながら任意の次元に写像する手法である。特に誤差が大きい
ラベルのみ描いた結果を Fig. 4.1-5 に示す。誤差 10%以上の測定点は点②、⑤（位置番号は Fig. 4.11 示されている）のみに生じていることがわかる。また、測定点⑤がより密集したクラスタリングになっていることから、誤差が生じる原因がハッキリしていると考える。それぞれについて結果の一例
を Fig.4.1-6 に示す。点②では不一致なだけだが点⑤では測定点が基板上ではなく発熱体から温度を読み取っていることがわかる。

58

Fig. 4.1-4. Visualized result of training data sets with error range labels.

Fig. 4.1-5. Visualized result of limited training data sets with high error range labels.

59

(a) Layout of heat objects in case of over 10% error at the validation position of ②.

(b) Layout of heat objects in case of over 10% error at the validation position of ⑤.
Fig. 4.1-6. Layouts for over 10% error cases: unexpected 360 cases exist at the validation pointof ⑤.
(4-1-2) 学習データのダウンサンプリングによる高精度化
発熱体との干渉が検証データで計 360 データ存在したため、修正データとして該当箇所を除外した入力で再度学習を行った。
ハイパーパラメータは同じであるが、検証データ数のみ 44,370 となる。
Table 4.1-3 および 4.1-4 より、各誤差評価で最大誤差約 5.9%／3.4℃、平均誤差も約 1.03%／0.406℃と前節の問題の解決ができた。また全体の約 99%が誤差±4%および±2℃未満となり、また決定係
数も 0.992 と代替モデルとして十分な精度が得られた。
重要な知見は、一見データに不整合がないように見えても、全体の挙動を表現するルールや法則から逸脱しているデータは学習を阻害する。さらに膨大なデータがあれば異なるルールも包含して学習する可能性はあるが、限られたシミュレーションの結果を学習させることを前提とすると、一般的に期待されている結果は得られないと考えるべきである。

60

Table 4.1-3. Summary of numbers of data sets corresponding with relative error ranges for
reduced training data sets.

誤差範囲 [±%]± 0%～± 2%

個数 [個]
38950

総数 [個]
38950

割合 [％]
87.785最大誤差 [%]± 2%～± 4%

4990

43940

11.2465.870± 4%～± 6%

430

44370

0.969

最小誤差 [%]± 6%～± 8%

0

44370

0

0.000

± 8%～± 10%

0

44370

0

平均誤差 [%]± 10%～

0

44370

0

1.027

Table 4.1-4. Summary of numbers of data sets corresponding with absolute error ranges for
reduced training data sets.

誤差範囲 [±℃]

個数 [個] 総数 [個] 割合 [％]± 0.0℃～± 1.0℃

41191

41191

92.835

最大絶対誤差± 1.0℃～± 2.0℃

2815

44006

6.344

3.366 [℃ ] / 6.119 [％ ]± 2.0℃～± 3.0℃

360

44366

0.811

最小絶対誤差± 3.0℃～± 4.0℃

4

44370

0.009

0.000 [℃ ] / 0.000 [％ ]± 4.0℃～± 5.0℃

0

44370

0

平均絶対誤差± 5.0℃ ～

0

44370

0

0.406 [℃ ] / 0.739 [％ ]Fig. 4.1-7. Comparison with reference and prediction temperature in all of validation data
sets for reduced training data sets: predictions of high accuracy are achieved.
(4-1-3) まとめ
本検討により以下の結論を得た。
•

オリジナルの発熱量のみではデータ量不足で過学習

61

•

パターンを与えると精度が大幅に上昇•

以前より使用している入力設計で学習・予測が可能•

今回の対象とする温度場の予測では、10,000 ケース程度は必要、つまり、ある程度温度のバリエーションも必要•

データ的に補外に該当する部分が原因で学習不良•

データ量の単純な増加でも傾向を表現することは可能•

より精度を出すには精度の悪い予測に対して、出現頻度などの統計量を詳細な確認が必要•

学習を阻害させる要因を可視化し特定することによりさらなる精度向上が達成可能

文献[1] C.M.ビショップ, パターン認識と機械学習 上/下 ベイズ推論による統計的予測, 2012, 丸善[2] 機械工学年鑑 2018, https://www.jsme.or.jp/kikainenkan2018/chap03.html#3-9（参照日 2021
年 4 月 1 日）
[3] E. Kaiser, et al, Sparse identification of nonlinear dynamics for model predictive control in
the low-data limit, Vol.474, Issue 2219, Proc. of the Royal Society A, 2018.
[4] C. Zhang, et al., Understanding deep learning requires rethinking generalization, 5thInternational Conference on Learning Representations, 15 pages, 2017[5] 和田義孝,深層学習によるき裂進展評価～計算力学サロゲートモデルの構築～, 保全学,Vol.182,2019,11-15.
[6] 和田義孝, 深層学習によるサロゲートモデル構築, プラスチック成形加工学会誌, Vol.32, No.3,
pp.83-86, 2020[7] 和田義孝, 設計問題適用を考慮した回帰モデル生成のためのデータ拡張, 第 25 回計算工学会
講演予稿集, 4 pages, 2020
[8] L. van der Maaten, G. Hinton, Journal of Machine Learning Research 9, 2579-2605, 2008
[9] 機械工学年鑑 2020, https://www.jsme.or.jp/kikainenkan2020/chap02/#a08（参照日 2021 年
4 月 1 日）[10] A. Choromanska, et al., The loss surfaces of multilayer networks, arXiv:1412.0233 [cs.LG],
13 pages, 2015[11] G. Montúfar, et al., On the Number of Linear Regions of Deep Neural Networks, Twentyeighth Conference on Neural Information Processing Systems, 9 pages, 2014(4-2) FFB と AI 連携活用
（4-2-1）風車後流のサロゲートモデル「富岳」を用いることにより直径 120m の洋上風車の解析が可能となり、複数台の風車の後流干渉の解析も可能となりつつあるが、計算負荷が大きいという問題は依然として残されている。特に洋上風車のように風車間の距離が長い場合は、風車の相互作用を調べるために長時間の解析が必要となる。風車の設計解析及び運用時の出力予測を行う場合は、この問題が顕在化してくる。そこで、これまでに得られた風車解析データを元に風車周りの流れをモデル化し、風車下流の流れを推定することを試みる。

62

Actuator disc

z

A0

y
x

A

Aw

p0 , U 0
pw , U w
Fig. 4.2-1. Stream-tube control volume of horizontal axis wind turbine rotor modelled by anactuator disc.
Fig. 4.2-1 に示すように風車前後の流れを流線で取り囲んだ流管でモデル化し、風車による速度欠損を式 4-2.1 で表す。

U
= C ( ( x, y, z )) f (r ( y, z ),  ( x, y, z ))U0

(4-2.1)

ここでＣは欠損の大きさを表し、f は速度欠損の形を表すとする。一般に風車の下流の流れは正規分布的な流れとなるが、風車の直下では二つの正規分布を重ね合わせたような形となることから、本研究では、二つの正規分布関数を組み合わせ、それが下流に行くにしたがって一つの正規分布に変化するようなモデルを考え、f を式 4-2.2 の形式で表す。

f (r ( y, z ),  ( x, y, z )) =−( r ( y, z)  r0 )1 DG+e+ e DG , DG  = −2
2 2 ( x, y, z )

(

2

)

(4-2.2)

このとき、風車に作用する推進力 T との間に以下のような関係が成り立つ。

T = U 02C ( ( x, y, z ))  (e DG +e DG −
+

Aw

−

+
+
−
C ( ( x, y, z )) 2 DG+(e+ e 2 DG + 2e DG + DG ))dAw
2(4-2.3)

式中に現れるσ、DG などの変数は未知変数であるが、大規模数値解析で得られたデータを元にこれらの未知変数を推定し、風車の後流を予測する。特に風車直下の二つの正規分布関数が一つになるまでの距離や周流の乱れを考慮し、f を以下のようにモデル化して推定を行った。
2
2





2
2


r0  r0   11  y 2 ( z − zh )1  y 2 ( z − zh )
f (r ,  ( x, y, z )) =  exp  −
+
+
+
exp
−
+
−
2
 2 22
 z2
  
 z2
   y 2   y



 





(4-2.4)（4-2-2）検証結果
Fig. 4.2-2 に検証結果を示す。実験結果[1]、他の研究者のダブルガウシアンモデル[1]、シングルガウシアンモデル[2][3]と本研究で提案した手法の比較結果を示す。風車直下では、シングルガウシア
ンモデルは実験との整合性が悪く、ダブルガウシアンでないと実験と一致しないことがわかる。風車から十分離れた位置ではどのモデルも後流分布を推定することができる。

63

Fig. 4.2-2. Lateral wake velocity profiles at hub height from several downstream distances.

64

主流の速度勾配と乱れ強度を考慮するため、風車上流の流れを Fig. 4.2-3 に示すように速度分布がある場合について検討を行った。Fig. 4.2-4 に示すように提案したモデルにより文献[4]の結果とほぼ同等な流れ場を得ることができた。Fig. 4.2-5 はコンター図で表した結果である。風車下流の任意の位置の流れ場を取得できることがわかる。Fig. 4.2-6、Fig. 4.2-7 に誤差評価の結果を示す。今回提案したモデルは風車直下の流れを精度良く求めることができることが示されている。洋上風車の場合、風車の直下に次の風車を設置しないことから、風車直下の流れ場をモデル化することがあまり行われてきていなかったが、後流の性質を詳細に調べ、蛇行などの問題を検討する際に風車直下の流れが重要であると推察される。大規模流体解析によって得られた風車下流の速度場データを用いることにより従来より高い精度で風車後流をモデル化することが可能となった。

Fig. 4.2-3. Vertical profile of incoming normalized atmospheric boundary layer (ABL) windvelocity [4].

65

66

Fig. 4.2-4. Vertical wake velocity profile at several downstream distances.

(a) DG prediction [4].

(b) Proposed DG prediction
Fig 4.2-5. Proposed DG prediction of contour of streamwise velocity on vertical plane.

67Fig 4.2-6. RMSE of lateral profile distribution of streamwise velocity

Fig. 4.2-7. RMSE of vertical profile distribution of streamwise velocity
文献[1] J. Schreiber, A. Balbaa, and C. L. Bottasso, “Brief communication: A double-Gaussian wake68

model,” Wind Energ. Sci., vol. 5, no. 1, pp. 237–244, Feb. 2020, doi: 10.5194/wes-5-237-2020.
[2] M. Bastankhah and F. Porté-Agel, “A new analytical model for wind-turbine wakes,” Renew.
Energy, vol. 70, pp. 116–123, 2014, doi: 10.1016/j.renene.2014.01.002.
[3] S. Xie and C. Archer, “Self-similarity and turbulence characteristics of wind turbine wakes
via large-eddy simulation,” Wind Energy, vol. 18, no. 10, pp. 1815–1838, 2015, doi:
10.1002/we.1792.
[4] W. Yu-Ting, L. Chang-Yu, H. Chien-En, and L. Shao-Dong, “Investigation of Multiblade
Wind-Turbine Wakes in Turbulent Boundary Layer,” J. Energy Eng., vol. 145, no. 6, p. 4019023,
Dec. 2019, doi: 10.1061/(ASCE)EY.1943-7897.0000625.
(4-3) FFR-Comb と AI 連携活用
(1-2) 節「パイロットスケール超臨界圧 CO2 ガスタービン」を参照のこと。
(4-4) RC HPC 版と AI 連携活用
洋上風車の解析では、中規模の多数のシミュレーションを行い、様々な風向の境界条件や諸条件を考慮して総合的な観点から設置検討を実施する。信頼性の高い予測結果を得るためには、小さな規模の計算でも精度のよい結果が得られることが望ましい。洋上風車の風車配置の問題については、後流干渉を精度良く再現することが鍵である。令和 2 年度は、少数の大規模計算から粗い格子でも細かい格子の高精度な後流干渉の状況を再現するサロゲートモデルの開発方法について検討した。サロゲートモデルの構築にあたり、現象を近似する方程式の構造モデルを探索するが、これにはデータから
目的変数、独立変数と従属変数の同定を行う必要がある。令和 3 年度は、この方程式の構造モデルを探索方法の検討を実施した。
多数のデータ変数間の関係性を調べる方法として、重回帰分析などをもちいた相関係数を用いる方法があるが、うまくいかない場合も多い。そこで、非線形の相関係数としての意味をもつ最大情報係数（MIC）を用いる。最大情報係数は、一般的な相関係数である「ピアソンの積率相関係数」と異なり、線形・非線形を問わずに相関関係を測定できる特徴をもつ。また、最大情報係数は「独立変数とその従属変数」、「同じ独立変数の従属変数２つ」の間で 1 に近づくという特徴がある。さらに、最大情報係数の値は回帰分析における決定係数 R2 に準ずる値をとる。
実験として、独立変数が 1 つ、従属変数が 3 つのローレンツ方程式（次式）のデータを作成し、そのデータから逆にローレンツ方程式の方程式構造が得られるかを調べた。
𝑑𝑥= −𝑝𝑥 + 𝑝𝑦
𝑑𝑡𝑑𝑦
= −𝑥𝑧 + 𝑟𝑥 − 𝑦𝑑𝑡
𝑑𝑧[ 𝑑𝑡 = 𝑥𝑦 − 𝑏𝑧
Fig. 4.4-1 に実験結果を示す。比較のため MIC の結果（左）と重回帰分析の結果（右）を示す。
MIC はその特徴どおり、独立変数とその従属変数、同じ独立変数による従属変数の関係を推定することができた。一方、重回帰分析では関係性を見いだすことはできていない。

69

Fig. 4.4-1. Comparison results of MIC(Left) and Multiple regression analysis(Right).
MIC を用いた方程式構造同定は効果のあることがわかった(5) 産業界等との連携及びコンソーシアムの構築(5-1) 燃焼システム用次世代 CAE コンソーシアムとの連携
令和 2 年度に引き続き、FFR-Comb への機械学習の導入、および「富岳」に対するチューニング
について情報交換を行った。令和 3 年 5 月 21 日、11 月 5 日にコンソーシアムの会合を開催した。
(5-2) 洋上風力関係機関との連携とコンソーシアムの構築（5-2-1）オープンプラットフォーム（5-2-1-1）風力発電運用高度化
近年のデジタル化は加速的に進化し、様々な産業で急速に変革が始まっている。実世界とサイバー世界の融合するインダストリ 4.0 や AI の普及による自動運転、状態診断などが代表的である。
風力の分野でも IoT の活用による運転保守コストの削減のための開発が行われている。Fig.5.2-1
は NEDO 風車運用高度化技術研究開発で検討されたデータプラットフォームである。本事業では、
風車稼働率を 97％以上に向上させる技術を確立させ、国内風車のダウンタイム及び運転維持コストの低減化、及び発電量向上を目指した技術開発を行うことで発電コスト低減を目的としている。

70

Fig.5.2-1. Example of data platform development in wind power [1].
本事業の社会的背景と成果の意義は次のようにまとめられており、以下引用する。
「国内風車メーカー不在の中、国内風車産業展開の一つとして、風車の O&M は事業コストの 35%以上の割合があり、産業経済効果が高く風車低コスト化に貢献できる重要な分野である。他方、O&Mの取り扱い範囲は風力発電設備全般に渡るため、また、技術開発進展も早いことから、関連国内企業の技術的な支援を可能とする社会システムつくり、O&M 企業への新たな知見、技術支援、そして有用な O&M 支援サービス、要素技術を有する企業の参入を後押しすることは、国内の風力発電事業の運用高度化において重要である。
欧州では風力発電メーカー以外でも、人員配置、拡張技術を発電事業に活用できており、低コスト化、信頼性向上へ寄与している。我が国の風力発電事業でも、将来を見据え最新技術の早期導入による産業競争力、発電コストの低減、信頼性の向上を実現できる開発が求められている。」（5-2-1-2）IoS-OP
船舶の分野において、海運データを活用するためのオープンプラットフォームが設立されており、風力産業コンソーシアム設立の参考事例とした。
IoS-OP（Internet of Ships Open Platform）は、海事業界におけるデータ流通を実現し、デジタル時代における新たな海事クラスターの形をつくり、次世代につなぐべく、データの創出・送受信・蓄積・活用など上流から下流までの作業を役割分担し、各社が得意分野に自由に参画できるデータ流
通基盤として設立された。海事産業がデータ利活用に注力し、新たな規制への対応やイノベーションの創出につなげられるよう、船舶 IoT データの流通・共有ルールを整備し、データ流通、共有、活用の拡大を目指している。
Fig.5.2-2 は IoS-OP の船舶データの利活用方法を示している。コンソーシアムの会員から提供されたデータを ShipDC（Ship Data Center）[2]において一元管理し、データの利用規約に従い、協調する領域と競争する領域に分けて、幅広いニーズにこたえるための仕組みとなっている。

71

Fig.5.2-2. Utilization of ship data of IoS-OP.
IoS-OP が、恒久的に秩序ある規則・規約のもと、会員に対するサービス提供を行うためには、公平性、信頼性および独立性のある活動を推進していく必要があるため、ShipDC では IoS-OP コンソーシアムと呼ぶ会員企業により構成される組織を設立し、各種ワーキング・グループ（WG）等の活動を運営している。
ShipDC は、経産省の「公的データ提供要請制度を活用できるデータ共有事業」に認定された初めての機関である[3]。
（5-2-2）IoW（Internet of Wind energy）オープンプラットフォームフォーラムについて
風力発電産業の発展のために、風車運転データをオープン/クローズを明確にした上で、ステークホルダー間で情報共有できるオープンプラットフォームを構築し、適切に運用していくことが必要となる。
風力産業コンソーシアム（仮称）の立ち上げに向けて、東京大学吉村教授、日本海事協会（ClassNK）の関係者と複数回の議論を経て、オープンデータプラットフォームを中心にした産官学のステーク
ホルダーからなる IoW（Internet of Wind energy）オープンプラットフォームフォーラム（IoW フォーラム）を立ち上げるための準備を始めた。以下に、IoW フォーラムの活動目的及び内容を示す。
（5-2-2-1）背景「洋上風力産業ビジョン（第１次）」[4]において、政府は導入目標を、また産業界は国内調達・コ
スト削減目標をそれぞれ掲げ、再エネの導入拡大と産業競争力強化の好循環を目指すこととしている。特に、サプライチェーンの形成等を通じて競争力を高めつつ、今後のアジア展開を見据えて、浮体式の商用化を含め、技術開発を加速化し、世界で戦える競争力を培っていく必要があり、「洋上風力産業の競争力強化に向けて必要となる要素技術を特定・整理し、「技術開発ロードマップ」[5]を令
和 3 年度に策定した。
（5-2-2-2）目的
円滑な国内調達及びコスト削減の推進には、デジタルトランスフォーメーション（DX）による、生産性向上に向けた弛まない取組みが必須であり、シミュレーション技術や 3D-CAD 導入によるデジタル設計手法を駆使し、初期計画、基本設計から詳細設計、生産設計ひいては保守管理、リサイク
ルまでのバリューチェーン全体を通して IoT 技術やデジタルデータの有効活用が極めて需要な要素72

となっている。このためには、データ解析、設計、施工、運転保守等の人材充足に困難な状況が懸念される中、これら関連作業の国内拠点連携やグローバル展開を可能とする、デジタル情報を扱う人材
を業界内外に広く求めてこれを活用することが必須である。今後、DX による生産性の飛躍的な向上が期待される中、デジタルツインなどに必須の設計及び状態監視情報のデジタル化をベースとした共有化及び各種作業のネットワーク化が鍵になり、また熟練技術者の経験及び暗黙知の作業プロセ
スへの取り込みが急務である。 この状況下、データ解析、設計人材、監理技術者等は育成、作業のデジタル化、ネットワーク化による生産性向上に向けた取組みにより、風力業界発展の一助とするところを志すものである。
（5-2-2-3）「IoW オープンプラットフォームフォーラム（仮称）」想定参加者
メンバー：NK、東京大学、日本気象協会、海技研の他、大学・研究所、発電事業者、製造者、
EPC・O&M 事業者
賛助メンバー：JWPA、JWEA、JWES、OEAJ、賛同する産官学から事務局：NK（IDD＋RED）＋WINC
会費：実費のみ（5-2-2-4）活動内容案
平等互恵の精神に基づき、関連作業のデジタル化、ネットワーク化及び人材育成を推進することを目的に、次の事項について連携・協力を推進する。
(1)情報交換及び発信(2)解析・設計・試験・研修サービス及びアプリの提供(3)システム・データベース及びデータプラットフォームの開発(4)海洋データ、設計データ、製造・建設・運転保守データ、撤去・解体データの収集・利活用(5)人材交流、人材育成(6)その他、本目的を達成するために必要な事項（5-2-2-5）準備状況
関係者らと本プロジェクト担当者らとの会合を開き、本プロジェクトで開発中のスーパーシミュ
レーション技術の紹介を行うとともに、風力産業界のニーズの聴収を行い、その結果を精査し、本シ
ミュレーション技術の開発に反映した。洋上風力発電コンソーシアム「Internet of Wind energy フォーラム（仮称）」の令和 4 年度設立に向け、日本海事協会をはじめ関係団体と協議を進め、アカデ
ミアの専門家を交えた第 1 回の準備会合を、令和 4 年 3 月 23 日に開催した。
参照[1] 東大・産総研・中部大、風車運用高度化技術研究開発、2021 年度成果報告会資料、2021[2] https://www.shipdatacenter.com/[3] 経済産業省、公的データ提供要請制度を活用できるデータ共有事業を初めて認定しました、
2019 年 6 月、https://www.meti.go.jp/press/2019/06/20190624001/20190624001.html[4] https://www.meti.go.jp/shingikai/energy_environment/yojo_furyoku/pdf/002_02_02.pdf[5]https://www.enecho.meti.go.jp/category/saving_and_new/saiene/yojo_furyoku/dl/roadmap/ro
admap20210401.pdf(6) 成果発表、広報・アウトリーチ活動
本プロジェクト全体として令和 3 年 11 月 26 日に国内シンポジウム「クリーンエネルギー「富岳」

73プロジェクトシンポジウム」（第 2 回、ポスト「京」重点課題⑥から数えると第 6 回目）をオンラインにて開催し、本プロジェクトのねらい、途中段階の成果を発表し、社会への広報・普及を行った。
特にカーボンニュートラルに向けて「富岳」と本プロジェクトがどのように貢献できるかという観点から特別講演を行い、大学や研究機関のみならず産業界からも多数の参加者が集まった。
そのほか令和 3 年 10 月 22 日に東京大学生産技術研究所革新的シミュレーション研究センター主催、本プロジェクト及び「航空機フライト試験を代替する近未来型設計技術の先導的実証研究」（課
題代表機関：東北大学大学院工学研究科）の共催で、「第 5 回 HPC ものづくり統合ワークショップ」
をオンラインにて開催した。
「「富岳」を利用した革新的流体性能予測技術の研究開発」（課題代表機関：東京大学生産技術研究所）との協働であり、ポスト「京」重点課題⑥と同⑧として活発に行っていた研究の協働を継続させている。

74

２－３．活動（研究会の活動等）
年月日

行事

令和 3 年 4 月 14、19 日

令和 3 年度第 1 回全体会

令和 3 年 10 月 22 日

第 5 回 HPC ものづくり統合ワークショップ

令和 3 年 11 月 26 日

クリーンエネルギー「富岳」プロジェクト第 2 回シンポジウム

令和 3 年 12 月 8、13 日

令和 3 年度第 2 回全体会

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

担当機関

担当責任者

（1）次世代火力発電システ

東京都文京区本郷 7-3-1 国立大学 教授 吉村 忍

ム及び洋上ウィンドファー

法人東京大学大学院工学系研究科

ムのスーパーシミュレーシ
ョン

京都府京都市西京区京都大学桂 C3（2）個別シミュレーション

棟 国立大学法人京都大学教授 黒瀬 良一と AI（機械学習）の連携活
用福岡県福岡市西区元岡 744（3）総合的事業推進

大学法人九州大学大学院総合理工国立教授 渡邊 裕章

学研究院
愛知県豊橋市天伯町字雲雀ケ丘１
番１

教授 飯田 明由

国立大学法人豊橋技術科学

大学大学院工学研究科
東京都目黒区駒場 4-6-1 国立大学

教授 加藤 千幸

法人東京大学生産技術研究所
福岡県春日市春日公園 6-1

国立大

准教授 内田 孝紀

学法人九州大学応用力学研究所
福岡県福岡市西区大字元岡 744 国立大学法人九州大学情報基盤研究
開発センター

75

教授 小野 謙二

別添１

学会等発表実績

別添１ 学会等発表実績
1. 学会誌・雑誌等における論文掲載

No. 掲載した論文（発表題目） 発表者氏名発表した場所

発表した時 国内・国 査読（有（学会誌・雑誌 期

際の別

名等）

りの場合
○ を記
入）

1

Numerical simulation of
kerosene jet in crossflow
atomization andevaporation under theelevated pressure andoscillating air-flow
condition

J. Wen, Y.
Atomization
2021年4月 国際
Hu, R. Kurose and Sprays, 31,
73-87 (2021)〇

2

A hyper-reduction
S. Kaneko, H. the Mechanics 2021年6月 国際
computational method for Wei, Q. He, and Physics ofaccelerated modeling of J.-S. Chen, S. Solids, Volume
thermal cycling- induced Yoshimura
151, June 2021,
plastic deformation
104385〇

3

Numerical analysis of
A. L. Pillai, T. Combustion
heat transfer
Murata, R.
and Flame,characteristics of spray
Kai, R.
Volume 239,
flames impinging on a
Masuda, T.
111615wall under CI engine-like Ikedo, R.
(2021)
conditions
Kurose2021年8月 国際〇

4Experimental andT. Honzawa, Energy and AI, 2021年9月 国際
numerical study of water R. Kai, K.
5, 100076
sprayed turbulent
Hori, M.
(2021)
combustion: Proposal of a Seino, T.
neural network modeling Nishiie,
for five-dimensional
R. Kurose
flamelet approach〇5

High-frequency
electromagnetic field
analysis using pseudoquadruple precision in
subdomain local solver

A. Takei, H.
Kawai, R.
Shioya, T.
YamadaJournal ofAdvanced
Simulation in
Science andEngineering,Vol. 8, No. 2,
pp.194-210(2021)

国際〇

Stability of the Cassie
state in a droplet
impinging on surfaces

Z. Yuan, M.
Matsumoto,R. Kurose

Physics of
2021年12
Fluids, 33,
月
123315 (2021)

国際

〇

6

76

2021年10月

with multiple holes

7

Numerical simulation of H.
Fuel, 318,
2022年1月 国際
binary-component fuelMuramatsu, 123331 (2022)
droplet evaporation using A. L. Pillai, K.
Level Set method
Kitada, R.
Kurose〇8

Direct numerical
Z. Yuan, M.
simulation of droplet
Matsumoto,impingement on surfaces R. Kurose
with hierarchicalstructures

International 2022年2月 国際
Journal ofMultiphase
Flow, 147,
103903 (2022)〇

9

A flamelet LES of
turbulent dense spray
flame using a detailed
high-resolution VOF
simulation of liquid fuelatomization

Combustion
and Flame,
237, 1117412(2022)

2022年3月 国際〇

10

スパコン版リアムコンパク 内田 孝紀、

日本風力エネル 2022年3月 国内〇トによる風車ウエイクの相 小野 謙二、

ギー学会論文集

J. Wen, Y.
Hu, T.
Nishiie, J.
Iino, A.
Masri, R.
KuroseVol.45 No.4
シミュレーション ーウエ 吉村忍、 加藤 pp.71-82, 2022
互干渉に関する大規模数値 飯田 明由、
イク領域内の平均風速欠損 千幸、 山出
量の予測精度検証ー

吉伸、 今村博、 植田祐子

11

Anisotropic doubleB. soesanto, Energy science 2022年4月 国際
Gaussian analytical wake T. Yoshinaga, andmodel for an isolated
A. lidaEngineerng
hodzontal axis windturbine〇2. 学会等における口頭・ポスター発表

No. 発表した成果（発表題目、 発表者氏名
口頭・ポスター発表の別）

発表した場所

発表した時 国内・国 招待講演

（学会名等）

期

際の別

（○を記
入）

1

データに基づく探索的な方 小野謙二

第38回関東CAE 2021年4月 国内

程式モデリング、口頭懇話会（オンラ
イン）2

流体励起振動を利用した圧 金子栄樹、吉 計算工学講演
電エナジーハーベストに関 村忍会、オンラインする連成解析、口頭77

2021年5月 国内○

3

ADVENTURE_Thermalコ 後藤聡太,、武 計算工学講演

2021年5月 国内

ードと非統計学的なアプロ 居周、 金子栄 会、オンラインーチによる不確かさ定量

樹、吉村忍

化、口頭
4

スパコン版リアムコンパク 内田 孝紀、

2021年度日本風 2021年5月 国内トによる風車ウエイクの相 小野 謙二、

工学会年次研究

互干渉シミュレーション、 飯田 明由、発表会(オンラ口頭吉村忍、 加藤 イン)
千幸、 山出吉伸、 今村博、 植田祐子

5

Numerical Study on
S. Kaneko, S. USNCCM, No. 2021年7月 国際
Active Control of Limit
Yoshimura16, Online
Cycle Oscillation with
Energy Harvesting、口頭

6

Optimization of structure J. Nagao, S. Asian Congress 2021年8月 国際of neural network
Terauchi, Y. on Gas
substituting flamelet
Takahara, R. Turbines 2020
approach using Genetic Kurose(ACGT2020),
OnlineAlgorithm, 口頭(Qingdau,
China)

7

Numerical investigation K.Kitada, J.
of evaporation phenomena Wen, R.
of liquid fuel jet
Kuroseatomization in crossflow

8

Numerical investigationof droplet evaporation
with CLSVOF method、ポスター

9

15th Triennial 2021年8月 国際International
Conference on
LiquidAtomizationand Spray
Systems(ICLASS 2021),
Online(Edinburgh,UK)

K. Kitada, J. Asian
2021年9月 国際
Wen, H.
Symposium on
Muramatsu, Computational
A. L. Pillai, R. Heat Transfer
Kuroseand Fluid Flow,
Online(Qingdao,
China)

Applicability of LES for L. Pillai, R.
investigating spray flame Kai, J, Li, R.
impinging on a wall under Kurose
compression- ignition

15th Triennial 2021年9月 国際International
Conference on
Liquid78

engine-likeenvironment、口頭

1011

Atomizationand Spray
Systems(ICLASS 2021),
Online(Edinburgh,
UK)

流体励起振動の能動制御と 金子栄樹、吉 日本機械学会第 2021年9月 国内
エナジーハーベスティング 村忍

34回計算力学講の同時問題、口頭

演会オンライン

非統計学的な不確かさ定量 後藤聡太,、武 日本機械学会第 2021年9月 国内化法とADVENTURE＿

居周、 金子栄 34回計算力学講

Thermalコードとの連携解 樹、吉村忍

演会オンライン

析、口頭
12

複数の並列ソルバーを連携 吉村 忍、山田 日本機械学会第 2021年9月 国内活用した石炭ガス化炉スー 知典、渡邊 裕 34回計算力学講
パーシミュレーション、口 章、黒瀬 良
頭演会オンライン一、吉田 隼也、淀 薫

13

Coupled Analysis for
S. Kaneko, S. JSST, No. 40,
Active Control andYoshimura
Kyoto, Japan
Energy Harvesting from
(Online)Flow-Induced Vibration、

2021年9月 国際

口頭
14Heat Transfer Analysis
with Uncertainty Using
Non-Statistical
Uncertainty
Quantification Method
and Parallelized Heat
Transfer Analysis Tool、

S. Goto, A.
Takei, S.
Kaneko, S.
Yoshimura

JSST, No. 40,
Kyoto, Japan(Online)H. Kawai

Proceedings of 2021年9月 国際
The 40th JSST
AnnualConference:
International
Conference on
Simulation
Technology,
Kyoto, Japan,
Sep. 1-3,
2021

2021年9月 国際

口頭
15High Performance
Implementation of
SkylineSolver for Many Core
Environment、口頭

79

16

17

Numerical simulation ofdroplet impingementon surfaces with microstructures、口頭

Z. Yuan,M.
Matsumoto,R. Kurose

日本流体力学会 2021年9月 国内
年会2021, オン
ライン

水素-酸素超臨界燃焼場に 徳岡泰輝、 甲 日本流体力学会 2021年9月 国内
おける発熱と渦の関連性に 斐玲央、 川村 年会2021, オン
関する数値解析、口頭

隆平、 黒瀬良 ライン
一18

Numerical simulation of H.
Asian
2021年9月 国際
binary fuel droplet
Muramatsu, Symposium on
evaporation with level set A. L. Pillai, K. Computational
Kitada, R.
Heat Transfer
method、ポスター
Kuroseand Fluid Flow,
Online(Qingdao,
China)

19

洋上風車後流の大規模乱流 吉村忍

第5回HPCもの 2021年10

解析、口頭づくり統合ワー 月

国内

クショップ
(online)20「富岳」における

宮村倫司

第5回HPCもの 2021年10

ADVENTUREの性能最適づくり統合ワー 月

化、口頭

クショップ

国内

(online)21

Investigation of the near B. soesanto, The 9th
2021年11
wake behavior of an
T. Yoshinaga, lntemational 月
isolated horizontal axis
A. lidaconference onsustainablewind turbine、口頭
EnergyEngineedngandApplication(ICSEEA 2021)

国際

22

Ghost Fluid 法を⽤いた⼆ 村松宏起、 P. 流体工学部門日 2021年11Abhishek L.、 本機械学会第99 月
成分燃料液滴の蒸発の数値
北田絢也、 黒 期流体工学部門
シミュレーション、口頭
講演会, オンラ
瀬良⼀
イン

国内

23

スカラー濃度場を用いた風 内田 孝紀、高 第43回風力エネ 2021年11

国内

車ウエイクの可視化とウエ 桑 晋、渡邊

ルギー利用シン 月

イクマネージメントへの適 慶一郎、長谷 ポジウム(オン用、口頭

川 聖矢、馬場 ライン)
好孝、村上 礼雄、山崎将80

英、肥高 邦彦

24

チャネル乱流内における予 甲斐玲央、 A. 第59回日本燃焼 2021年11混合火炎-壁相互作用のL. Pillai、 U. シンポジウム,

DNS、口頭Ahmed、 N. オンライン

国内

月

Chakraborty、 黒瀬良一25

LSOR-PCR updates and
PCR-like TPR、口頭

小野謙二、光 サイエンティフ 2021年11田達哉

国内ィック・システ 月
ム研究会
A64FXシステム
アプリ性能検討
WG第6回会合（オンライン）26

データ駆動型のモデリング 小野謙二

理研シンポジウ 2021年11

タ研究、口頭ム: 計算で物事 月

国内を理解する予測する～産業界の実問題に立ち向かうサイエンス～（オンライン）
27

燃焼・ガス化炉のスーパー 吉村忍「スーパーシミ 2021年11シミュレーション 全体会ュレーションと 月

用、口頭AI を連携活用

国内した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)28

FFR-CombによるCO2回収 渡邊裕章

「スーパーシミ 2021年11型IGCC用ガス化炉解析、ュレーションと 月

口頭

AI を連携活用した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)

81

国内

29

ADV_Solidの「富岳」チュ 山田知典

「スーパーシミ 2021年11ーニング＆非弾性解析、口ュレーションと 月

頭AI を連携活用

国内した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)30

洋上風力発電のロードマッ 今井博「スーパーシミ 2021年11プと「富岳」への期待、口ュレーションと 月

頭AI を連携活用

国内した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)31

RIAM-COMPACT 「富内田孝紀、小 「スーパーシミ 2021年11岳」チューニングと風車ウ 野謙二

ュレーションと 月

エイク解析、口頭

AI を連携活用

国内した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)32

FronFlow/blue (FFB )加藤千幸
Oprimization for Fugaku
and Its Resutls、口頭「スーパーシミ 2021年11国内ュレーションと 月
AI を連携活用した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)33

洋上ウィンドファームスー 陳順華、吉村 「スーパーシミ 2021年11
パーシミュレーション流体 忍ュレーションと 月

構造連成＆疲労損傷解析、

AI を連携活用

82

国内

口頭

した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)34

洋上ウィンドファームのマ 飯田明由

「スーパーシミ 2021年11

ルチフィジクスシミュレー

ュレーションと 月

ション、口頭AI を連携活用

国内した実機クリーンエネルギーシ
ステムのデジタ
ルツインの構築と活用」第2回
シンポジウム(online)35

連成現象のスーパーシミュ 吉村忍

第3回スパコン 2021年12レーション -基本的な考え

コロキウム、オ 月

方と事例-、口頭

ンライン

国内

36

Combustion noise
A. L. Pillai, S. 13th Asia2021年12generated by an open
Inoue, T.
Pacific
月lean-premixed low-swirl Shoji, S.
Conference onhydrogen flame: A hybrid Tachibana, T. Combustion
Yokomori, R. 2021
CFD/CAA study、口頭
Kurose(ASPACC2021), Online (AbuDhabi, UAE)

国際

37

超臨界燃焼場におけるLOX 和田 祥、 甲 第35回数値流体 2021年12国内

coreの挙動の数値解析、口 斐玲央、 黒瀬 シンポジウム,
頭
38

良一

○

月

オンライン

燃焼数値シミュレーション 黒瀬良一

日本ガスタービ 2022年1月 国内の基礎と最新動向 －如何

ン学会 第50回に使いこなすか －、口頭

ガスタ－ビンセ○

ミナー, オンラ
イン
39Super-simulation of
S. Yoshimura The 4th R-CCS 2022年2月 国際
Coal Gasification Facility
International
Symposium,
on Fugaku, 口頭
online○40

カーボンニュートラル社会 黒瀬良一

第5回CAEワー 2022年3月 国内○の実現に貢献する燃焼数値

クショップ、オシミュレーション、口頭

ンライン

83

41

カーボンニュートラル社会 黒瀬良一

「富岳」を利用 2022年3月 国内の実現に向けた燃焼数値シした革新的流体

ミュレーションの役割、口性能予測技術の頭研究開発 第２○回「富岳」流体
予測革新プロジェクトシンポジ
ウム、東京大学
生産技術研究所、オンライン42

先見情報を用いない創発的 小野謙二

令和3年度 大規 2022年3月 国内な方程式推定、口頭

模データ解析と○人工知能技術によるがんの起源と多様性の解明
シンポジウム（オンライン）43

並列有限要素法のための領 荻野正雄

第5回 大規模電 2022年3月 国内

域分割法とBDD前処理、

磁界数値解析手

口頭法に関する研究
シンポジウム(LSCEM2022)、八戸工業大学
メディアセンター

44

メニーコア環境におけるス 河合浩志

第5回 大規模電 2022年3月 国内

カイラインソルバーの前進

磁界数値解析手

消去後退代入の高速化、口法に関する研究頭

シンポジウム(LSCEM2022)、八戸工業大学
メディアセンター45

NUMERICAL
F. Ansai, Y.
PREDICTION OF WAKE Suzuki, Y.
OF ANAIRFOIL BY
Miki and C.
USING THE LARGEKatoEDDY SIMULATION、口the 32nd
2022年3月 国際International
Symposium on
Transport
Phenomena○

実機クリーンエネルギープ 渡邊裕章

日本冷凍空調学 2021年9月 国内○

ラントの大規模スーパーシ会 「地球温暖ミュレーション、口頭化に対応するた頭
46

めの先進熱交換

84

技術に関する調査研究」第6回
委員会
3. 受賞等

No

名称

受賞者氏名

授賞機関
(学会名等)

.

受賞した
時期

国内・国

備考

際の別

4. メディアへの情報発信、ウェブサイト等での情報公開

No.

名称

日付

説明

開催日時

開催場所

備考

5. 広報活動等(ワークショップ・研究会等の開催

No. 名称
１

クリーンエネルギー「富岳」プロジェ
クト 第2回シンポジウム

参加者（人数）
参加者271名、参2021年11月26 オンライン 加機関等160機関
日２

第5回HPCものづくり統合ワークショ
ップ（東京大学生産技術研究所 革新的シ
ミュレーション研究センター 主催、「『富岳』を利用した革新的流体性能
予測技術の研究開発」（課題代表機

2021年10月22 オンライン 参加者192名、参関：東京大学生産技術研究所）、「航 日空機フライト試験を代替する近未来型設計技術の先導的実証研究」（課題代
表機関：東北大学大学院工学研究科）、および「『富岳』が拓く
Society 5.0時代のスマートデザイン」（課題代表機関：理化学研究所計算科
学研究センター） と本プロジェクト
の共催）

85

加機関等104機関