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

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

-1-

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

-2-

目次
目次 ............................................................................................................................................................ i
１．補助事業の目的 ................................................................................................................................. 1
２．令和 2 年度（報告年度）の実施内容 ................................................................................................ 1
２－１．当該年度（令和 2 年度）の事業実施計画 .............................................................................. 1
２－２．実施内容（成果） ................................................................................................................... 2(1) 各アプリケーションの「富岳」向けチューニング .................................................................... 2(2) ラボスケール石炭ガス化炉（電中研炉）スーパーシミュレーション ....................................... 8(3) ベンチスケール石炭ガス化炉（三菱実験炉）スーパーシミュレーション .............................. 22(4) パイロットスケール超臨界圧 CO2 ガスタービン（東芝エネルギーシステムズ製）スーパーシミュレーション ................................................... 25(5) 単体及びタンデム配置大型風車のスーパーシミュレーション ................................................ 30(6) 洋上ウィンドファーム全体のスーパーシミュレーション ....................................................... 41(7) 個別アプリケーションと AI の連携活用 .................................................................................. 53(8) 産業界コンソーシアムや外部機関等との連携.......................................................................... 68(9) 成果発表、広報・アウトリーチ活動 ........................................................................................ 78
２－３．活動（研究会の活動等） ...................................................................................................... 80
２－４．実施体制 ............................................................................................................................... 80
別添１ 学会等発表実績 ........................................................................................................................ 81i補助事業の名称「富岳」成果創出加速プログラム
スーパーシミュレーションと AI を連携活用した実機クリーンエネルギーシステムのデジタルツインの構築と活用１．補助事業の目的
ポスト「京」重点課題⑥プロジェクトにおいて「富岳」での活用を念頭に開発準備を進めてきたアプリ
ケーション群（ADVENTURE、FFB、FFR-Comb、REVOCAP_Coupler、RC HPC 版）に基づくマル
チフィジクス・マルチスケール統合シミュレーション（本事業においてスーパーシミュレーションと称する） と AI（機械学習）を連携活用し、Society5.0 を支える電力システムの主要クリーンエネルギーシ
ステムとなる CCS（Carbon Capture and Storage）技術の適用に適した次世代火力発電システム、及び洋上ウィンドファームのデジタルツインを構築し、実機エネルギーシステムの開発、設計、運用の最適化に活用できるようにする。本事業において構築・活用されるアプリケーション群は、いずれも我が国産業界における実活用の実績も豊富な高度に汎用化された力学系 CAE アプリケーションであり、上記のクリーンエネルギーシステムのみならず、燃焼機器、伝熱機器、流体機械等の幅広い実機に活用可能である。
また、本事業は、様々な企業群、産業界コンソーシアムとも強く連携しながら遂行される。以上のことから、本事業は、「富岳」を活用した我が国の産業競争力の強化にも幅広く大きく貢献する。

２．令和 2 年度（報告年度）の実施内容
２－１．当該年度（令和 2 年度）の事業実施計画(1) 次世代火力発電システム及び洋上ウィンドファームのスーパーシミュレーション
・ポスト「京」重点課題⑥から解析を進めている、ラボスケール石炭ガス化炉（電中研炉）の解析について、運転条件を変更した追加解析を行い、結果をとりまとめ、国際ジャーナル論文を執筆し投稿する。
（東大・九大・京大）
・ポスト「京」重点課題⑥から構築を進めている、ベンチスケール石炭ガス化炉（三菱実験炉）の解析
モデルと粗計算格子メッシュ（燃焼部は完成済、冷却管付き炉容器部、連成界面部）を完成し、モデルの検証を行う。
（東大・九大）・ポスト「京」重点課題⑥と⑧において、それぞれのコデザインのターゲットアプリであった(a)ADVENTURE 及び(b)FFB について、令和元年に引き続き、「富岳」実機上で性能検証を継続し、同時に最終チューニングを行い、「富岳」上での実機解析の準備を整える。（東大・豊橋技科大）・ポスト「京」重点課題⑥において「富岳」での適用を念頭に、「京」及び名大 FX100、九大 ITO、京大 Cray XC40 等において解析機能拡張及び高速化チューニングを行ってきた (a)FFR-Comb、(b)REVOCAP_Coupler に関して、「富岳」実機上での性能検証を行い、「富岳」向けチューニングを完成させ、「富岳」上での実機解析の準備を整える。
（東大・京大・九大）
・HPCI 第二階層マシンを用いて、粗計算格子のメッシュを用いて、(a)ベンチスケール石炭ガス化炉（三菱実験炉）、(b)パイロットスケール超臨界圧 CO2 ガスタービン（東芝エネルギーシステムズ製）、(c)
洋上ウィンドファーム（Alpha Ventus）の全体と(d)大型風車の後流を介した相互干渉のマルチフィジク
ス・マルチスケール統合シミュレーション（スーパーシミュレーション）の予備解析を行う。
（東大・京

1

大・九大・豊橋技科大）(2) 個別シミュレーションと AI（機械学習）の連携活用
・各個別アプリケーション（(a)ADVENTURE、(b)FFB、(c)FFR-Comb、(d)RC HPC 版）の解析結果と AI（機械学習）を連携活用し、それぞれのサロゲートモデル構築に関する準備研究を行う。
（東大・京大・九大・豊橋技科大）(3) 総合的事業推進（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) 各アプリケーションの「富岳」向けチューニング
(1-1) ADVENTURE
令和 2 年度は並列非線形解析コード ADVENTURE_Solid2（以下 ADV_Solid2）の「富岳」向けチューニングの実施に向けた事前検討として ADV_Solid2 プロトタイプ実装を用いて各コアアルゴリズム（BDD、Local Solver 等）の計測及び最適化を実施した。その結果 BDD に関しては、線形代数ソルバーの OpenMP 対応化を進め、また、コース行列のスカイライン記憶のブロック化と BLAS、 LAPACK ラ

2

イブラリの効率的な利用等により、億レベル自由度の大規模板曲げ問題においてピーク性能比 20%程度を実現した。また、Local Solver に関しては、BDD と同様にこれまでの 1 領域 1 コアの割当てを 1 領域
複数コアに変更し、それに伴い並列スカインソルバーを導入することにより、億レベル自由度の大規模板曲げ問題においてピーク性能比 15%を実現した。これらの成果を ADV_Solid2 に移植し、引き続き「富岳」上での最適化を進めた。
一方、並列熱伝導解析コード ADVENTURE_Thermal（以下 ADV_Thermal）の「富岳」向けチュー
ニングに向け、各コアアルゴリズム（BDD、Local Solver 等）の計測を行った。まず、ADV_Thermal の公開版を「富岳」上でコンパイルし、実行できることを、いくつかの数千万要素規模の問題を用いて確認した。その結果、ADV_Thermal に実装されている BDD は、ADV_Solid2 に実装されている BDD と比較して効きがよくなくかつ計算負荷が軽いことがわかった。そこで、ADV_Thermal に関しては、最も計算量が多い、「部分領域ソルバーの前進・後退代入及び行列ベクトル積」の高速化に主眼をおくこととした。
(1-2) FFB
FrontFlow/blue(FFB)は Large Eddy Simulation（LES）をベースとする並列有限要素法流れ解析コードである。FFB の最大の特長は、乱流中の微小な渦のダイナミクスを直接計算することにより乱流を高精度に計算できることにある。大型発電用風車周りの流れにおける微小な渦のダイナミクスを計算する
ために、高解像度の計算格子が必要となり、「富岳」を用いて大規模な計算を短時間に実行することが求められる。このために「富岳」の性能を最大限引き出すコードの高速化技術が求められる。FFB を「富岳」で高速に動作させるため、「富岳」の開発と並行し、「コデザイン」のもと、東京大学生産技術研究所
および理化学研究所が連携し、FFB を「富岳」で高速動作させるための高速化技術が開発された。
「「富岳」成果創出加速プログラム『「富岳」を利用した革新的流体性能予測技術の研究開発』」では、上記の高速化技術を FFB に実装し、「富岳」上で高い計算性能を確認している[1-3]。本研究テーマでは、「富岳」
を用いた風車まわり流れの実証計算の準備として高速化された FFB をベースに風車まわり流れ計算に必要となる機能を取り入れカスタマイズした。FFB の高速化の基本部分は、「「富岳」を利用した革新的流体性能予測技術の研究開発」において実施された。ここでは、高速化された FFB の「富岳」における計算性能と、風車まわり流れ計算のために拡張した機能について報告する。
（1-2-1）高速化した FFB の「富岳」における計算性能
まず、「「富岳」を利用した革新的流体性能予測技術の研究開発」（代表者：加藤千幸）において行われた FFB の高速化の基本部分を記す。
「富岳」の CMG（Core Memory Group）に 50 万要素の直方体領域を割り当て、これを直列につなげた計算領域を対象とする Weak Scale ベンチマークテストを実施した。
使用したノード数は 158,976 ノード（7,630,848 コア）である。高い並列性能を維持するため、通信の隠蔽処理も実装している。すなわち、通信節点部分の演算を先行処理し、通信命令を発行するとともに、通信処理と同時に内部の演算を実行する計算方法を採用している。Fig.1.2-1 に Weak Scale ベンチマーク
テストの結果を示す。並列性能は 48 ノード計算時を基準としているが、300 万コアあたりまで並列化効
率 85%を維持している。全体の実行性能は、22.6 PFLOPS を達成した。この性能は流体解析や構造解析
といった連続体系のアプリケーションの性能を評価するための HPCG ランキングにおいて「富岳」が達3成した計算性能 16.00 PFLOPS（2 位の SUMMIT は 2.93 PFLOPS）を上回る性能（ただし HPCG は倍精度計算、FFB は単精度計算）である。

Fig.1.2-1. Parallel-computing performance of FFB run on Fugaku in weak-scale benchmark test.
（1-2-2）風車まわり計算のための FFB のカスタマイズ
風車の出力は、大気境界層や流入風の乱流強度によって変化するため、風車まわり計算においても大気境界層や流入風の乱流強度をユーザ入力により制御する機能が必要であった。このため、これらの機能を「富岳」用にカスタマイズした FFB に実装し、その動作テストを実施した。流入に乱れを入れる機能の動作確認として、Fig.1.2-3 に示すような簡易なパイプ内部流れを対象に動作確認を行い、パイプ内部流れにおける主流速度の変動スペクトルを確認した。緑は入口、青は入口から 0.1D（D はパイプ直径）の内部点における主流速度の変動スペクトルである。図に示す通り、入口近く（緑のデータ）では、パワースペクトルにおいてピークが生じており非物理的なスペクトルになっているが、パイプ内部の変動スペ
クトルに関しては-5/3 則（コルモゴロフのマイナス 3 分の 5 乗則）に従う変動スペクトルが確認されることにより、実装した機能を用いると計算領域内に適切な乱れが生成されていることを確認した。また、
風車まわり流れにおける動作確認として、NREL（国立再生可能エネルギー研究所 （The National
Renewable Energy Laboratory））5MW 風車単体を対象としたテスト計算を実施した。Fig.1.2-4 に風車
上流 1D における主流方向速度の鉛直分布を示す。乱れを流入した速度場（緑）は、境界条件で与えた大気境界層分布（赤）にそって、複数の波数を含む乱れが付加されていることより、大気境界層、乱流乱れのどちらに関しても、意図した分布や乱れが形成されており、機能が正常に実装されていることを確認し
た。

4

wall

+

inlet ×free*

Fig.1.2-3. Test for function to set intensity of turbulence at inlet boundary (above: computationalmodel, below: power spectra of streamwise velocity).

5

(a)

(b)

(c)

(d)

Fig.1.2-4. Vertical profiles of streamwise velocity in LES of flow around wind turbine (NREL-5MW)
with considerations of atmospheric boundary layer and inlet turbulence at Z/D= (a) -1.0, (b) -0.5, (c)0.5, (d) 1.0. (Note that wind turbine is located at Z/D=0).
文献[1] K. Kumahata, S. Inoue, and K. Minami, “Kernel performance improvement for the FEM-based
fluid analysis code on the K computer,” Procedia Computer Science, vol. 18, pp. 2496-2499, 2013.K.
[2] K. Kumahata, K. Minami, and N. Maruyama, “High-performance conjugate gradientperformance improvement on the K computer,” International Journal of High Performance
Computing Applications, vol. 30, no. 1, pp. 55-70, 2015.

6

[3] C. Kato, Y. Yamade, K. Nagano, K. Kumahata, K. Minami and T. Nishikawa, “Toward
Realization of Numerical Towing-Tank Tests by Wall-Resolved Large Eddy Simulation based on 32
Billion Grid Finite-Element Computation”, Proceeding SC20, 2020.
(1-3) FFR-Comb
FFR-Comb に対して、「富岳」向けのチューニングを実施した。主に Bi-CGSTAB ルーチンの precondition 処理について、LU 分解をローカル化し、対角項の探査ループについて上三角ループと下三角
ループを一体化するよう修正を行った。
ベンチマークとして石炭ガス化反応計算の事例を使用し、チューニング実施後の FFR-Comb の性能と「京」における FFR-Comb の性能を比較した。石炭ガス化反応計算の計算格子は 1 億 1800 万要素である。計算条件として、乱流モデルに Dynamic LES を、流体モデルにゼロマッハ数近似を使用し、Euler
二相流を計算した。
「富岳」と「京」それぞれについて 100 ステップの試計算を行い、プロファイル情報を取得した。
「京」では 9216 並列（1152 ノード、ノード当たり 8 プロセス）、「富岳」では 9216 並列（192ノード、ノード当たり 48 プロセス）で計算した。
ノード当たりの GFLOPS は「京」
で 1.47 GFLOPS/node、「富岳」では 7.95 GFLOPS/node となり、「京」と比較すると 5.4 倍の性能となった（Table1.3-1）。
また「京」と「富岳」それぞれで仮にシステム全体を使用したと仮定すると、「京」は 82944 ノード、「富岳」は 158976 ノードであるため、「富岳」では 5.4×(158976/82944)=10.4 倍の性能が得られる試算となった。
Table 1.3-1. Performance comparison of FUGAKU and K computer
GFLOPS/node

倍率

並列数

「京」

1.47

1.0

9216 並列(1152 ノード×8 プロセス)

「富岳」

7.95

5.4

9216 並列(192 ノード×48 プロセス)「京」との性能比較のために Euler 二相流の解析を行ったが、「富岳」上では Lagrange 粒子の計算性能がポイントとなるため粒子計算のテストを行った。ベンチマークとして微粉炭燃焼計算の事例を使用
し、「富岳」上での Lagrange 粒子計算の有無による性能を比較した。微粉炭燃焼計算の計算格子は 105
万要素である。計算条件として、乱流モデルに Dynamic LES を、流体モデルにゼロマッハ数近似を使用
し、6 か所の粒子投入口から順次粒子を投入し計算を行った。
「富岳」上で粒子の計算を行わない場合と粒子の計算を行う場合について、それぞれ 1000 ステップの試計算を行い、プロファイル情報を取得した。
いずれも 192 並列で計算した。ノード当たりのプロセス数はメモリの制約から最大 24 プロセスとなった(粒子あり 1)。ノード内プロセス数の影響を確認するためノード当たり 4 プロセスでも計算を行った(粒子あり 2)。粒子あり 2 の方が粒子あり 1 に比べてノード内データのメモリ転送が速く行われたものと考えられる。また、投入粒子数の影響を確認するため投入粒子数をステップ当たり 1 粒子だけではなく、100 粒子の計算も行った(粒子あり 3)。それぞれの結果を Table 1.3-2 にまとめた。
「富岳」上では粒子あり計算の負荷により粒子なし計算に比べて GFLOPS 値が 3 割から 5 割減少することが確認できた。以上
から実スケール計算を「富岳」で行う上での粒子計算の取り扱いに関する知見を得ることができた。

7

Table1.3-2. Particle calculation performance on FUGAKU
GFLOelapse time [s]
倍率
並列数
PS
粒子なし

479.2

84.7

1.0

粒子あり 1

768.1

57.8

0.68

粒子あり 2

649.8

68.3

0.81

粒子あり 3

1626.141.3

0.49

192 並列(8 ノード×24 プロセス)
192 並列(8 ノード×24 プロセス)
1 particle/step injection
192 並列(48 ノード×4 プロセス)
1 particle/step injection
192 並列(8 ノード×24 プロセス)
100 particle/step injection(1-4) REVOCAP_Coupler「京」で使用していた FFR-Comb、ADV_Thermal と REVOCAP_Coupler を連携させた双方向連成
解析モデルをベースとして、「富岳」上でビルドを行うとともに、動作確認を行った。この解析モデルでは「京」対応上、FFR-Comb を起動した後に、FFR-Comb から REVOCAP_Coupler と ADV_Thermalを起動する仕組みを採用していたが、「富岳」の MPI 仕様では、MPI プロセス内で MPI_Comm_spawn
等を呼び出す際に、必要以上にメモリを獲得しメモリ不足になるという課題が判明し、途中で止まってしまった。この課題の「富岳」システム側では 2021 年 9 月頃にならないと解決しないということなので、
FFR-Comb から REVOCAP_Coupler と ADV_Thermal を呼び出す処理を変更するという対応を進めることとし、「富岳」上で動作することを確認した。
(2) ラボスケール石炭ガス化炉（電中研炉）スーパーシミュレーション(2-1) FFR-Comb による単体解析
令和 2 年度は、令和元年度ラボスケール石炭ガス化炉の 100 万要素の簡易格子を用いて燃焼モデルについてチューニングされた FFR-Comb を用い、
約 1 億要素の詳細格子における Verifivation & Validation（V&V）を実施した。
(2-1-1) 解析方法
Eulerian-Lagrangian 法による固気二相反応流の LES[1,2]解析を行う。気相の支配方程式は、連続の式、運動量保存式、エネルギー保存式、及び化学種保存式であり、これらを解いた。SGS 応力モデルには Dynamic Smagorinsky モデル[3]を用いた。固相については、個々の粒子の BBO 方程式をパーセル近似により解き、流体と粒子の干渉効果については、PSI-CELL 法[4]により考慮した。ガス化炉内の化学反応としては、石炭の熱分解による揮発分放出、石炭チャーの固気（ガス化）反応、及び気相反応を考慮した。揮発分放出過程を Eq. (1)により表し、揮発分組成については、Umemoto らの実験データ[5]を用いた。揮発分放出速度については、Eq. (2)に示す単一一次反応モデルにより考慮し、その反応速度パラメ
ータは FLASHCHAIN モデル[6]により見積もった。
𝐶𝑜𝑎𝑙 → 𝛼1 𝐶𝐻4 + 𝛼2 𝐻2 + 𝛼3 𝐶𝑂 + 𝛼4 𝐶𝑂2 + 𝛼5 𝐻2 𝑂 + 𝛼6 𝐶ℎ𝑎𝑟

(1)

𝑑𝑉𝑖

(2)

𝑑𝑡

𝐸

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

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

n th-order モデル[7]を用いた。それぞれの反応経路の反応速度パラメータは文献[8]の値を用いた。
1

𝐶ℎ𝑎𝑟 + 2 𝑂2 → 𝐶𝑂

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

(5)

𝑑𝑥
𝑑𝑡

𝑘11 𝑃𝐶𝑂2= 𝜂(1 − 𝑥)√1 − Ψln(1 − 𝑥) (1+𝑘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)𝐶𝐻4 + 𝑂2 → 𝐶𝑂 + 2𝐻2
1

2

𝐻2 + 𝑂2 → 𝐻2 𝑂

(8)

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

𝐶𝐻4 + 𝐻2 𝑂 → 𝐶𝑂 + 3𝐻2(10)

𝐶𝑂 + 𝐻2 𝑂 ↔ 𝐶𝑂2 + 𝐻2

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

(12)

2
1
2

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

9

(2-1-3) 解析結果
Fig. 2.1-2 に、燃焼モデルチューニング前後の単体計算、及び連成計算による中心軸上ガス温度分布の計算結果を示す。同図より、チューニングにより単体計算の解析精度が格段に向上していることを確認できる。
Fig. 2.1-3 に、チューニング前後の単体計算、及び連成計算による出口ガス組成の計算結果を示す。出
口ガス組成についても、中心軸上ガス温度と同様に、チューニングにより単体計算の解析精度が格段に向上していることを確認できる。
なお、同図において連成計算が実験値から大きく外れているが、この連成計算はチューニング前の
FFR-Comb コードを用いたものであるためである。燃焼モデルチューニング後の FFR-Comb コードを用いた連成解析の精度は大きく改善するものと思われる。この連成解析の再計算は(2-2)にも詳述するように「京」停止後は九大 ITO を用いて計算していたが、計算に時間がかかりなかなか終了しないため、途中で「富岳」に切り替えたが、(1-4)に述べた理由に解析モデルが動かなかったため、連成解析における複数コードの起動法を変更した。。
(2-1-4) 運転条件変更時の V&V
次に、ラボスケール石炭ガス化炉の運転条件を変更した際の解析精度について検討を行った。解析条件
は、（一財）電力中央研究所から提供された試験データに基づいて設定した。Table 2.1-1 に運転条件についてまとめる。ここで、表中 Case C は、FFR-Comb の燃焼モデルチューニングに用いた条件であり、Case A と B は、ガス化剤中の CO2 濃度の影響を調べるために、CO2 濃度をそれぞれ 0 %と 15 %に設定したものである。Case A はいわゆる空気吹き条件に一致している。また、各条件は炉内温度が等しくなるように、酸素比が調整されている。
Fig. 2.1-4 に、Case A、B、及び C の中心軸上ガス温度分布の実験値と単体計算結果との比較を示す。
炉内温度が等しくなるように酸素比が調整されているため、運転条件間でのガス温度の差異は大きくはないものの、いずれの運転条件についても中心軸上ガス温度分布の計算結果は実験値と良い一致を見ていることがわかる。
Fig. 2.1-5 に、Case A、B、及び C の出口ガス組成の実験値と単体解析結果の比較を示す。出口ガス組成については、ガス化剤中の CO2 濃度が異なることから、石炭のガス化反応性に大きく影響を及ぼし、運転条件間での差異が大きいが、計算結果はこの傾向を良く捉えられていることがわかる。
Fig. 2.1-6～8 に、Case A、B、及び C のガス温度、及び主要ガス成分の炉内分布を示す。上述のとおり、運転条件の変更によってもガス温度分布に大きな差異は見られないものの、Case A→B→C と、コン
バスタ内の CO2 濃度が上昇するに従い、各主要ガス成分も変化する様子が捉えられている。
以上により、ラボスケール石炭ガス化炉単体計算における V&V は一定の成果を収めることができたものと考える。

10

(a) Computational(b) Computational grid.

domain.

Fig. 2.1-1. Computational domain and grid.

Fig. 2.1-2. Comparison of gas temperature distribution in the axial direction.

11

Fig. 2.1-3. Comparison of product gas composition at the exit.
Table 2.1-1. Cases performed.
Cases

Unit

Case A

Case B

Case C

Oxygen ratio

-

0.52

0.55

0.54

Coal input

kg/h

100

100

100

O2

25

25

25

CO2

Vol%
Vol%

0

15

25

N2

Vol%

75

60

50

Fig. 2.1-4. Comparison of gas temperature distribution in the axial direction for Cases A, B, andC.

12

Fig. 2.1-5. Comparison of product gas composition at the exit for Cases A, B and C.

(a) Temperature(b) H2(c) CO(d) CO2(e) H2O
Fig. 2.1-6. Gas temperature and concentrations of major gases for Case A.

13

(a) Temperature(b) H2(c) CO(d) CO2(e) H2O
Fig. 2.1-7 Gas temperature and concentrations of major gases for Case B.

(a) Temperature(b) H2(c) CO(d) CO2(e) H2O
Fig. 2.1-8 Gas temperature and concentrations of major gases for Case C.
文献[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 simulationof 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 turbulenceand 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 gas-

14

droplet 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 andH2O 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 coal gasifier,Fuel, 85, pp.1935-1943, 2006.
[8] H. Watanabe, K. Tanno, H. Umetsu, S. Umemoto, Modeling and simulation of coal gasification
on an entrained flow coal gasifier with a recycled CO2 injection, Fuel, 142, pp.250-259, 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 discrete-ordinates
method, J. Thermophys., 2, pp.309-316, 1988.

(2-2) FFR-Comb – REVOCAP_Coupler – ADV_Thermal の連成解析と V&V
令和元年度に実施した FFR-Comb の燃焼モデルチューニング結果を受け、(2-1)節と同様に燃焼モデルのチューニングを行った FFR-Comb を使用し、FFR-Comb ⇔ REVOCAP_Coupler ⇔ ADV_Thermal
による電中研炉双方向連成解析を実施した。連成にあたり、FFR-Comb（流体側）からは熱流束を送信し、ADV_Thermal（構造側）からは壁面温度を送信するものとし、熱流束と温度の通信は各 Time step ごとに実施するものとした。
計算は九大 ITO 上において実施し、並列数は Table 2.2-1 のように設定した。
Table 2.2-1. number of process
並列数
FFR-comb

2304 並列(64 ノード×36 プロセス)

REVOCAP_Coupler

1697 並列(48 ノード×36 プロセス)

ADV_Thermal

540 並列(15 ノード×36 プロセス)

合計

4541 並列(127 ノード×36 プロセス)

令和 2 年度末の時点において本連成計算は途中過程の 39000step までとなっており、結果の評価にはさらに計算を続行する必要があるが、計算に問題が生じていないかを確認するため同一の Time step における単体解析と連成解析の結果を比較した。比較は Fig. 2.2-1～3 に示すように、断面内温度分布、壁面近傍（第一格子点）温度分布、中心軸上温度分布に関して行った。いずれの比較においても、連成解析の結果において若干温度が低い結果となったものの定性的な分布は近似しており、計算は問題なく進行していると判断される。引き続き計算を続行し、単体解析の結果と合わせ総合的に V&V を行う。

15

Fig. 2.2-1. Temperature distribution (cross section).

16

Fig. 2.2-2. Temperature distribution (neighborhood of the wall).

Fig. 2.2-3. Temperature distribution (center line).
令和 2 年度はさらに、「京」で実施していた FFR-Comb ⇔ REVOCAP_Coupler ⇔ ADV_Thermal に17

よるラボスケール炉（電中研炉）の燃焼流-伝熱連成解析を九大 ITO 及び「富岳」において実施した。九
大 ITO では先に述べたように「京」で使用した双方向連成解析システムをそのまま使用して、正常に動作することを確認した。一方、「富岳」においては、(1-4)節でも述べたように「富岳」の MPI 処理系に問題があり、ソフトウェアの起動処理に関して改修が必要となった。従来の連成解析システムでは、「京」のジョブ実行時のファイルステージング仕様に対応するため、FFB-Comb 内部から MPI_Comm_spawn関数を使用して REVOCAP_Coupler 及び ADV_Thermal を起動する処理を行っていた。しかし、「富岳」の MPI 処理系では MPI プロセス内で MPI_Comm_spawn を使用するとメモリの異常使用を引き起こすバグがあることが判明し、メモリ不足でジョブが異常終了する問題が発生した。「京」と違って「富岳」
ではファイルステージングは行われないため、FFB-Comb 内部から MPI_Comm_spawn 関数を使用してREVOCAP_Coupler 及び ADV_Thermal を起動する必要がなくなった。そこで、「富岳」上での双方向連成解析実行にあたっては、MPI_Comm_spwan による起動処理を廃止し、ジョブスクリプトから FFRComb、RECOVAP_Coupler、及び ADV_Thermal を連続的にバックグラウンド実行するように改修を行い、「富岳」上で動作することを確認した。
(2-3) ADV_Thermal の伝熱・冷却管解析
ADV_Thermal による伝熱・冷却管連成解析コード中の、連成境界面での物理量交換処理にバグがあり、一部の冷却管について不適切な冷却水の温度分布が出力される問題が発見されたため、改修を行った。
使用したモデルは、テスト作業で使用したラボスケール石炭ガス化炉（電中研炉）の炉容器の 3 次元
ソリッドメッシュ(四面体 1 次要素、要素数 155,999,061、節点数 25,510,852)をそのまま使用した。境界条件として、燃焼室内壁表面部分に燃焼に伴う熱移動を模した熱流束境界条件を与えた。また、４本の冷却管モデルには流速一定の条件を課したほか、全ての冷却管の冷却水流入温度として 24.77℃の温度規定
境界条件を与えた。炉容器側の初期条件として(2-2) 節に述べた FFR-ADV 連成解析で得られた 17,999
ステップ目の節点温度及び節点熱流束ベクトルを与えた。炉容器側の時間進展パラメータとして、Δ
t=0.001[s]を与えた。冷却管ソルバー（Discontinuous Galerkin[DG]法で離散化）内のサブサイクルとして 1 時間ステップ当たり 1,000 サブサイクルの計算を全ての冷却管に対して行うように設定した。連成
解析のため、冷却管表面にソルバー内部で熱伝達(参照温度=流体温度、熱伝達係数 0.062)を模擬した熱流束境界条件を動的に与えた。
Fig.2.3-1 ～ Fig.2.3-4 に改修前のコードで得られた冷却水の温度分布をグラフ化したものを示す。冷
却水の温度分布のグラフは、いずれも横軸が 1D モデルの節点座標(単位:mm)、縦軸が冷却材温度(単位:℃)である。ADV_Thermal に組み込まれている 1D 冷却管ソルバーは相変化には対応していないため、流体温度が 100℃以上になった状態での気化及び 0℃以下になった状態での凍結に関しては無視した計算結果である点に注意が必要である。また、18,000 時間ステップ目の温度及び熱流束を初期条件として計算している関係上、ステップ番号も 18,000 からの連番で示した。

18

Fig. 2.3-1. Temperature distribution of coolant in cooling pipe A at time step 20,000.

Fig. 2.3-2. Temperature distribution of coolant in cooling pipe B at time step 20,000.

19

Fig. 2.3-3. Temperature distribution of coolant in cooling pipe C at time step 20,000.

Fig. 2.3-4. Temperature distribution of coolant in cooling pipe D at time step 20,000.
冷却水の温度分布は、炉容器内の燃焼部からの熱によって、定性的には冷却管入口（Fig.2.3-1~4 の左側）から出口（Fig.2.3-1~4 の右側）に向かって連続的に上昇するのが自然である。しかし、特に Fig.2.3-

20

4 においてはのこぎりの歯のような不連続な温度分布や最大で-140℃程度の明らかに異常な低温度の箇所がみられるなど、コード改修により一部の問題は解消したものの、いまだに不適切な温度分布が出力されているため、対応を継続中である。
(2-4) ADV_Solid2 の構造(非弾性・熱応力)解析(2-4-1) はじめに高温での金属材料においては、ひずみ硬化とともに熱回復によって、一定ひずみ速度あるいは一定応力下でクリープ現象が生じる。有限要素法においてこの非線形材料挙動を解析するには、適切な弾塑性・ク
リープ構成則の導入が必要となる。これまでの開発において Ohno-Wang（OW）モデルおよび Norton 則に基づく弾塑性・クリープ則が ADVENTURE_Solid2 に導入されている。
ベンチスケール石炭ガス化炉（三菱実験炉）の解析においては高温運転状況において、時間変動する熱荷重を考慮に入れる必要が生じる。令和 2 年度の開発では、ヤング率のような材料物性値が温度に依存する際の初期試行応力の計算方法を改良し、熱弾塑性計算の応力計算の精度を向上するための検討および開発を行った。
(2-4-2) OW モデルの試行応力計算の改良設計(2-4-2-1) 従来アルゴリズムの課題
ADV_Solid2 の非線形解析における応力計算ルーチンは以下のような式により応力を更新していた。

1(𝑡+ Δ𝑡)2: Δ𝜺𝑒𝝈(𝑡+Δ𝑡) = 𝝈𝒕 + 𝑫(1)

1(𝑡+ Δ𝑡)2: (Δ𝜺 − Δ𝜺𝑝 −Δ𝜺𝑡ℎ )= 𝝈𝒕 + 𝑫

Δ𝜺𝑡ℎ = 𝛼(2)

1(𝑡+ Δ𝑡)2
𝛥𝑇𝑰

ここで、𝑫は弾性テンソル、Δ𝜺は全ひずみ増分、Δ𝜺𝑝 は塑性ひずみ増分、Δ𝜺𝑡ℎ は熱ひずみ増分、𝛼は線膨
1張係数、𝛥𝑇は温度増分、𝑰は単位テンソルである。時刻𝑡 + Δ𝑡の値で𝑫, αをサンプリングしている理由は
2

応力積分の時間積分を中点則で行っているためである。
このような定式化は簡易に実装ができる反面、１荷重増分の間にひずみが変化せず増分熱ひずみがゼロとなるような状況において、温度が変化してヤング率のみが変わった場合には、応力がヤング率に応じて変化しないという課題が存在していた。これは(1)式では常に試行応力の初期値は𝝈𝒕 を元にしているためである。このことにより、ひずみの変化がない場合にヤング率が小さくなった時、𝝈(𝑡+Δ𝑡)は𝝈𝒕 より小さくなることが不可能になっていた。

(2-4-2-2) 改良アルゴリズム
前項の課題を解決するには、応力更新を全量形にして実装する方法が考えられる。

21

𝝈(𝑡+Δ𝑡) = 𝑫(𝑡+Δ𝑡) : 𝜺𝑒(𝑡+Δ𝑡)(𝑡+Δ𝑡)=𝑫(𝑡+Δ𝑡): (𝜺(3)
𝑝 (𝑡+Δ𝑡)−𝜺𝑡ℎ (𝑡+Δ𝑡)−𝜺)𝜺𝑡ℎ (𝑡+Δ𝑡) = 𝛼 (𝑡+Δ𝑡) (𝑇 (𝑡+Δ𝑡) − 𝑇 (0) )𝑰(4)

この場合、𝜺は全ひずみ、𝜺𝑝 は塑性ひずみ、𝜺𝑡ℎ は熱ひずみのそれぞれ全量形となっている。また、本定式化では中点則を用いていないため、弾性テンソル𝑫および線膨張係数αは時刻𝑡 + Δ𝑡の値をサンプリン
グするものとする。またこの場合、線膨張係数αは熱ひずみ－温度曲線の割線の傾きを想定したものとなっている。この式ではひずみ一定で保持し、ヤング率を小さくすれば𝝈(𝑡+Δ𝑡) は𝝈(𝑡)に比べて小さくなることがわかる。
(3)式のままでは現状の ADV_Solid2 の応力更新の定式化との差異が大きくなってしまうため、現状のコードに沿う形で(3)式の効果を取り入れる方法を考える。(1)式において𝝈𝒕 を𝑫(𝑡+Δ𝑡) : 𝑫(𝑡)−1 𝒕

𝝈 に置き換えると、−1

𝝈(𝑡+Δ𝑡) = 𝑫(𝑡+Δ𝑡) : 𝑫(𝑡) 𝝈𝒕 + 𝑫(𝑡+Δ𝑡) : Δ𝜺𝑒= 𝑫(𝑡+Δ𝑡) : (𝜺𝑡 − 𝜺𝑝 𝑡 −𝜺𝑡ℎ 𝑡 ) + 𝑫(𝑡+Δ𝑡) : (Δ𝜺 − Δ𝜺𝑝 −Δ𝜺𝑡ℎ )(5)= 𝑫(𝑡+Δ𝑡) : (𝜺(𝑡+Δ𝑡) − 𝜺𝑝 (𝑡+Δ𝑡) −𝜺𝑡ℎ (𝑡+Δ𝑡) )を得る。これは(3)式と同様のものとなっていることがわかる。これより、既存のコードを活かしながら(3)式を実装するには(1)式の𝝈𝒕 を𝑫(𝑡+Δ𝑡) : 𝑫(𝑡)−1 𝒕

𝝈 に置き換えればよいことがわかる。これを行うコード

部分は試行応力の初期値を計算する部分となっている。令和 2 年度の改良では(6)式のような試行応力の計算スキームを追加実装した。
𝑡+𝛥𝑡 (𝑇 𝑡+𝛥𝑡 ):𝝈𝑡+𝛥𝑡(𝑫𝑡 (𝑇 𝑡 ))−1 : 𝝈𝑡𝑡𝑟𝑦0 = 𝑫

(6)

𝑡+Δ𝑡
𝑡+𝛥𝑡
𝝈𝑡+Δt
(𝑇 𝑡+Δ𝑡 ): 𝚫𝜺𝑡𝑟𝑦 = 𝝈𝑡𝑟𝑦0 + 𝑫(2-4-3)まとめ令和 2 年度の開発では ADV_Solid2 に対して応力積分の試行応力計算アルゴリズムの改良を行った。
これまでの試行応力計算の方法では、ひずみを保持したまま温度変化に伴いヤング率が変化した際の応力の変化について対応できないという課題が存在していた。そこで令和 2 年度の開発では応力計算の試
行応力の初期値計算を改良することでヤング率が変化した際の応力計算の精度を上げることが可能になった。
(3) ベンチスケール石炭ガス化炉（三菱実験炉）スーパーシミュレーション(3-1) FFR-Comb による単体解析
令和元年度に、ラボスケール石炭ガス化炉（電中研炉）の簡易格子（100 万要素格子）を用いてチューニングされた FFR-Comb を用い、ベンチスケール石炭ガス化炉（三菱実験炉）の 1 億要素簡易格子を用いて V&V を実施したところ、出口ガス組成の主要成分であるメタンの予測精度に改善の余地があることを見出した。そこで令和 2 年度は、メタンの予測精度向上をねらった反応モデルチューニングを約 600

22万要素の超簡易格子により実施するとともに、1 億要素簡易格子における V&V を再度実行し、その有効性を確認した。
(3-1-1) 解析方法
解析方法は、(2-1-1)節で述べたものとほぼ同様である。ただし、唯一、下記のメタン水蒸気改質反応について、順反応に加えて逆反応も考慮するようにコードを修正した。
（旧）CH4 + H2O ⇒ CO + 3H2(1)（新）CH4 + H2O ⇔ CO + 3H2(2)(3-1-2) 解析結果
Fig. 3.1-1 に、解析の結果得られたガス化炉（三菱実験炉）出口ガス組成の比較を示す。ただし、値は実験値との相対値である。旧反応モデルではメタンの生成を考慮できなかったため、ガス化炉出口においてメタン濃度はゼロであったが、Eq. (2)のように、今回新反応モデルとしてメタン水蒸気改質反応に逆反応を考慮することにより、ガス化炉出口のメタン濃度は実験値に近い値を予測することができるようになった。この反応モデルの改良により、生成ガスの発生量を適切に捉えることができるものと考えられる。
次に Fig. 3.1-2 に、ガス化炉内に 2 点設置された熱電対によるガス温度の実験値との比較を示す。ただし、熱電対位置、及びガス温度は代表値により正規化してある。ガス温度については、旧反応モデルによっても良い精度で予測が可能であったが、同図よりさらに精度が向上していることがわかる。
以上により、FFR-Comb の解析モデル（反応モデル）の改良を実施し、1 億要素の簡易格子による予備解析は完了した。

Fig. 3.1-1. Comparison of product gas composition at the exit.

23

Fig. 3.1-2. Comparison of gas temperature at two axial points.
(3-2) FFR-Comb – REVOCAP_Coupler – ADV_Thermal による連成解析
FFR-Comb ⇔ REVOCAP_Coupler ⇔ ADV_Thermal によるベンチスケール石炭ガス化炉（三菱実験炉）の燃焼流－伝熱双方向連成解析を「富岳」において実施した。
ベンチスケールモデルのうちの炉容器サブモデルは、要素数 406,317,033、節点数 67,112,837 のメッシュ（４面体一次要素）である。伝熱解析では 1 節点 1 自由度なので、全体の自由度は 67,112,837 となる。経験的に ADV_Solid2 の階層型領域分割法による解析では、1 部分領域当たり 1000～5000 自由度程
度で最も収束性がよくなることがわかっている。今回の解析ではその中間の 3000 自由度を目標値として領域分割パラメータを設定した。その結果、パート数は 35 ノード×48 コア=1,680 パート、1 パート当たりの部分領域数は 14 となった。これらの領域分割パラメータから、改めて部分領域当たりの自由度数の概数を再計算したところ、67,112,837 / 1,680 / 14 = 2854 自由度となった。また、FFR-Comb が解析を担うベンチスケールモデルのうちの炉内部の燃焼流サブモデルについても、別途分割数の検討を行い、106,642,914 要素 / 20,680,539 節点のメッシュを 2340 プロセスに分割することとした。
燃焼流－伝熱の連成界面は燃焼室内壁部とした。FFR-Comb による燃焼流解析で得られた連成界面における熱流束値を FFR-Comb から REVOCAP_Coupler に送信し、REVOCAP_Coupler で補間して
ADV_Thermal に渡された値を熱流束境界条件として ADV_Thermal において伝熱計算が行われるようにモデル作成を行った。また、燃焼流サブモデル、伝熱サブモデルともに、それぞれ FFR-Comb 及びADV_Thermal による単体解析の結果をもとに解析条件を設定した。
これらの準備を行い設定した解析モデル及び解析システムを「富岳」上で構築し、計算を行ったところ、(1-4)節で述べたように、「富岳」の MPI 処理系に問題があり、動かなかった。
「富岳」の MPI 処理系の修正は令和 3 年 9 月頃になるということだったので、解析システム側でソフトウェアの起動処理に関して改修を行うこととした。すなわち、従来のシステムでは、「京」のジョブ実行時のファイルステージン
グ仕様に対応するため、FFB-Comb 内部から MPI_Comm_spawn 関数を使用して REVOCAP_Coupler
及び ADV_Thermal を起動する処理を行っていた。
しかし、
「富岳」
の MPI 処理系では MPI_Comm_spawn

24

を使用するとメモリの異常使用を引き起こすバグがあり、メモリ不足でジョブが異常終了する問題が発生した。「富岳」ではファイルステージング機構がないため、FFB-Comb 内部から MPI_Comm_spawn関数を使用して REVOCAP_Coupler 及び ADV_Thermal を起動する必要がなくなった。そこで、
MPI_Comm_spwan による起動処理を廃止し、ジョブスクリプトから FFR-Comb、
RECOVAP_Coupler、及び ADV_Thermal を連続的にバックグラウンド実行するように改修を行った。
(4) パイロットスケール超臨界圧 CO2 ガスタービン（東芝エネルギーシステムズ製）スーパーシミュレーション(4-1) 目的
令和 2 年度は、
乱流燃焼モデルの一種である Flamelet 法に AI（機械学習）を連携させた Flamelet/ANN（Artificial Neural Network）法を確立し、この手法を用いた LES をパイロットスケール（商用機の 10分の 1 程度）の東芝エネルギーシステムズ株式会社製 50MWth 級超臨界 CO2 タービン用燃焼器内の乱流燃焼場[1-4]に適用することによって、同手法の有用性を確認することを目的とした。なお、令和 2 年
度は HPCI 第二階層マシン（京大 Cray）を用いたため、燃焼器の上流部のみを計算対象とし、要素数を最小限に抑えた計算を実施した。
(4-2) 計算方法
(4-2-1) 燃焼モデルと支配方程式
本解析では、Flamelet 法として火炎の浮き上がりや消炎を再現可能な Flamelet/progress-variable（FPV）法[5,6]をベースとし、壁面からの熱損失及び実在実気体を考慮可能に拡張したものを用いた[7]。
支配方程式は、連続の式、運動量保存式（Navier-Stokes 方程式）
、混合分率 Z の保存式、反応進行変数

C の保存式、及びエンタルピーh の保存式で構成される。
(4-2-2) 3 種吹込みガス燃焼の Flamelet 法
従来の Flamelet 法では、燃料と酸化剤の二つの独立したガス系統からの流入混合による燃焼のみに対応可能であり、燃料、酸化剤に加えて冷却用 CO2 の吹込み等、三つ以上の独立したガス系統からの吹込みがある条件に適用することができなかった。そこで、今回、3 種吹込みガス条件に適用可能な Flamelet

法の開発を行った。その概念を Fig. 4-1 に示す。

Fig. 4-1. Concept of three stream Flamelet model [8].
従来の Flamelet 法では、燃料由来成分の質量分率を示す混合分率 Z を輸送することで、燃焼場における燃料と酸化剤の混合具合を求め、その値を基に Flamelet データベースを参照し、燃焼場を再現する。

25

これに対して、今回開発を行った 3 種吹込みガス燃焼の Flamelet 法では、燃料由来の質量分率（図中のFuel）を Z1、Cooling CO2 の質量分率 (図中の Coflow)を Z2 と定義し、これら二つの輸送方程式を解く
ことで、3 系統の流入ガスの混合具合を判断する。今、酸化剤の質量分率を Z0 として、Cooling CO2(Coflow)を酸化剤の一部と定義すると、酸化剤に含まれる Cooling CO2 (Coflow)の割合 W は、
𝑍

(4-1)

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

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

𝐶 = 𝑌𝐻2𝑂 + 𝑌𝐻2で与える。
(4-2-3) ニューラルネットワーク（ANN）による Flamelet モデルの構築

今後益々増加しうる Flamelet データベースの情報量の増加に対応するため、多次元の Flamelet 空間"2 , 𝑊
̃ , 𝐶̃ , ∆𝐻̃ ）の関数で表現すにおける多様体であるデータベースを基底ベクトル（本解析においては𝑍̃, 𝑍̃
ることを考える。こうすることで流体計算時の必要なメモリがデータベースそのものからパラメータベ
クトルだけとなり、メモリの大幅な節約となる。すなわちこれは Flamelet データベースに対する回帰問題である。
機械学習という枠組みの中でも回帰問題を解くための手法は数多く存在するが、ここではディープラーニングの手法、すなわちニューラルネットワークを用いた関数の表現を採用する。実際に本作業で用いるニューラルネットワークでは、入力層のパーセプトロンの数はデータベース空間の基底ベクトル数である 5、出力層のパーセプトロンの数は 1 である。隠れ層をどのように設計するかはハイパーパラメータ
であってそれ自体が最適化の対象となり得る。一般的に層が深いほど、層あたりのパーセプトロンが多いほど再現可能な関数形状は多くなり、データベースの再現性が高くなることが考えられる。しかし一方で、ネットワークを大きくすれば燃焼計算においてデータベースの基底ベクトルから燃焼状態を推定する際の計算負荷が高くなる。そこで、ネットワークの構築においてはデータベースの一定の再現性を担保しつつもネットワークができるだけ小さくなるように留意した。本解析では、学習の収束条件として、データベースの値とネットワークによって作成した関数の値の RMSE が、データベースの物理量の標準偏差の 4.5%未満となるようにした。この条件の下で試行錯誤によって、ネットワーク構成は隠れ層が 7 層、１層あたりのパーセプトロンは 9 とした。各パーセプトロンの出力は活性化関数 Swish[9]で変換され、オプティマイザーRadam[10]によって重みが更新される。
(4-2-4) 計算対象及び計算格子
Fig. 4-2 に本解析で対象とした超臨界燃焼器の形状（青枠）を示す。対象としたのは、東芝エネルギーシステムズ製 50MWth 級検証試験用超臨界 CO2 タービン及び燃焼器[1-4]である。本ガスタービン燃焼
器は、酸化剤／燃料(CH4)／再循環 CO2 により燃焼が行われる。実機の運転条件はタービン入口圧力約26

30MPa の超高圧条件を想定しているが[4]、今回対象とした燃焼試験は設備容量の制限から最大 10MPa
で行っている。また、本解析では計算負荷の低減のため、赤枠で囲った燃焼器上流部のみを対象に解析を行う。

Fig. 4-2. Computational domain.
Fig. 4-3 に解析格子形状の詳細の様子を示す。点数及びセル数は、それぞれ 930 万点及び 3,350 万個である。解析格子下流部には計算安定性のための遷移領域を設け、逆流を抑えるために出口を絞っている。

Fig. 4-3. Computational grids.

27

(4-2-5) 計算条件
本解析では、NuFD/FrontFlow/Red（FFR-Comb）[5-7,11,12]を用いた熱流体―構造連成の LES を実施する乱流モデルとしては、標準 Smagorinsky モデルを用いた。圧縮性流体スキームを用い、実在気体モデルとして Soave-Redlich-Kwong(SRK)モデル[13]を用いた。計算は、Table 4-1 に示す四つの条件について行った。ケース１では、2 段総括反応モデル（以下に示す 2 段の Westbrook CH4 反応モデル[14]）を用いた。
𝐶𝐻4 + 3⁄2 𝑂2 ⇒ 𝐶𝑂 + 2𝐻2 𝑂
𝐶𝑂 + 1⁄2 𝑂2 ⇔ 𝐶𝑂2
Table 4-2 に Case 2-4 で用いたデータベースの概要を示す。Case 2 では、データベース容量の問題か
ら W の分割を 2 と最低限の分割数としているものの、データベース容量は 1.1GB と非常に大きくなる
ことがわかる。Case 4 のように W の分割数を６とさらに上げた場合、通常のデータベースを用いるとデータ容量はおおよそ 3.3GB となり通常の大規模計算機の 1 コア当たりのメモリ容量を超え、解析が不可能となることが予想される。一方、メモリ容量の問題を解決するため ANN 法を用いた機械学習を用いた"2 , 𝐶̃ , ∆ℎ̃ ) の 4 次元データベースを W の分割分用意し、それをCase 3, 4 では、各 W に対して𝜙̃ = 𝜙̃(𝑍̃, 𝑍̃

W 方向に補間するようにしている。その結果、データベースに必要なデータは各層間の重み関数パラメータのみであり、データ容量がそれぞれ 0.48×10-3GB(480kB)、1.44×10-3GB(1440kB)と大幅にデータ
ベース容量を削減することができていることが確認できる。

Case

Table 4-1 Cases performed in this study.
Combustion model

Reaction Mechanism

1

2 step global reaction

2step global

2

3 injection non-adiabatic FPV(W =2)

GRI-3.0

3

3 injection non-adiabatic FPV/ANN(W =2)

GRI-3.0

4

3 injection non-adiabatic FPV/ANN(W=6)

GRI-3.0

Case

Table 4-2 Flamelet database.
Dimension

2
3
4

5 dimensions:"2
̃1 × 𝑍̃
̃
̃
𝑍
×
𝑊
×
𝐶
× ∆ℎ = 100 × 25 × 2 × 100 × 5
1
4 dimensions×2:"2
̃1 × 𝑍̃̃̃(𝑍
1 × 𝐶 × ∆ℎ) × 𝑊 = 100 × 25 × 100 × 5 × 2
4 dimensions × 6:
"2
̃1 × 𝑍̃̃̃(𝑍
1 × 𝐶 × ∆ℎ) × 𝑊 = 100 × 25 × 100 × 5 × 6

Size
1.1Gb
0.48E-3Gb(480kb)
1.44E-3Gb(1440kb)(4-3) 結果及び考察
Fig. 4-4 に、燃焼器内ガス温度の瞬間分布を比較して示す。同図から、Case 3（3 種吹込み非断熱 FPV／ANN 法）の結果は、Case 2（3 種吹込み非断熱 FPV 法）の結果と概ね一致しており、ANN の採用が良好に機能していることがわかる。Case2,3 の W を 2 分割(W = 0.0,1.0)から、Case 4 では、W =0.0,0.1,0.2,0.3,0.4,1.0 と分割数を 6 に増やしたことで、解析結果の精度改善が確認できた。しかし、冷28

却用の再循環 CO2 が多く含まれ燃焼が緩慢となる W = 0.4～1.0 の間のデータが現在ないため、その Wの範囲での結果の精度が不十分であり反応が十分に完了せず未燃燃料が出口で見られた。

(a) Case 1(b) Case 2(c) Case 3(d) Case 4
Fig. 4-4. Comparison of instantaneous distribution of gas temperature.
(4-4) まとめ令和 2 年度は、超臨界燃焼場に適用可能な LES/Flamelet/ANN 法の確立に成功し、その手法が有効であることを、パイロットスケールの東芝エネルギーシステムズ株式会社製 50MWth 級超臨界 CO2 タービ
ン用燃焼器内の乱流燃焼場に適用することによって実証した。
文献[1] Sasaki, T., Itoh, M., Maeda, H., Tominaga, J., Saito, D., Niizeki, Y., “Development of Turbine andCombustor for a Semi-closed Re-cuperated Brayton Cycle of Supercritical Carbon Dioxide”, In Proc.
of the ASME 2017 Power Conference Joint With ICOPE-17, POWER-ICOPE2017-3419, 2017.

29

[2] Allam, R., Fetvedy, J., Forrest, B., Jones, C., Nomoto, H., Itoh, M., “A novel, high-efficiency, oxyfuel power plant with low-cost electricity production and 100% capture of carbon dioxide”, POWERGEN International 2013.
[3] 岩井保憲, 伊東正雄, “超臨界 CO2 サイクル発電用ガスタービン燃焼器”, 東芝レビュー, 70 (5), 2015.
[4] Iwai, Y., Itoh, M., Morisawa, Y., Suzuki, S., Cusano, D., Harris, M., “Development approach tothe 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 simulationof turbulent spray combustion in a subscale aircraft jet engine combustor - Predictions of NO and sootconcentrations -", 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 andnumerical investigation of thermo-acoustic instability in a liquid-fuel aero-engine combustor atelevated pressure: validity of large-eddy simulation of spray com-bustion", Combustion and Flame,
162, 2621-2637, 2015.
[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, 13091317. 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 of asupercritical 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.
[14] Westbrook, C.K., Dryer, F.L., “Simplified Reaction Mechanisms for the Oxidation of
Hydrocarbon Fuels in Flames", Combustion Science and Technology, 27, 31-43, 1981.
(5) 単体及びタンデム配置大型風車のスーパーシミュレーション(5-1) 解析条件の整理
FFB 及び RC HPC 版によるシミュレーションの解析条件を以下に示す。

30

問題 A：NREL5MW 風車１台を配置したモデルを用いて、FFB から ADVENTURE 用流体力時刻歴を出力する（累積疲労損傷解析の感度解析）
Table5.1-1～4 に各問題における解析条件を示す。解析順に示す番号は、解析の優先順位を表す。ここで、NREL5MW 風車の運転カーブ（回転数(RotSpeed,rpm)、ブレードピッチ(BIPitch1)、発電機トルク(GenTq,kN-m)、周速比(TSR)）を Fig.5.1-1 に示す。解析条件は、Fig. 5.1-1 の運転カーブに基づいて決めた。

風速
ロータ回転数ブレードピッチ
周速比(TSR)ヨー角

Table5.1-1. Simulation conditions for case A.
[m/s]
6.0
10.0
11.4（Reference）
[rpm]
8.0
12.0
12.1
[deg.]
0.0
0.0
0.0
[-]
8.8
7.92
7.0
0.0
[deg.]
15.0

ウィンドシア
乱流強度
解析順

0.1
0.19
(Cat.B)③

0.1
0.16(Cat.B)
②

0

0.1

0.0

-

①-1

①-2

0.15
(Cat.B)
①-3

12.0
12.1
4.0
6.65

0.1
0.15
(Cat.B)④

Fig. 5.1-1. NREL 5MW operational curve.
問題 B：NREL5MW１台風車モデルの Wake 計算のための FFB、RC-HPC 解析条件
前方配置の１台とするが、後方に伝わる Wake の様子を計算できるように後方に解析領域を長くとる。

風速
ウィンドシア
乱流強度
解析順Table5.1-2. Simulation conditions for case B.
[m/s]
11.4
0.1
0.15
①

31

問題 C：NREL5MW２台風車のタンデム配置モデルの FFB、RC-HPC 解析条件
事前に問題ＡやＢから得られた後方風車直前位置の中心風速を使って、後方風車の回転速度を変更した解析とする。ただし、計算中に流入風速に応じた回転速度の変更はしない。

風速
ウィンドシア
乱流強度
離隔距離
解析順

Table5.1-3. Simulation conditions for case C.
[m/s]
11.4
0.1
0.15
3D
①

7D
③

問題 D：Alpha Ventus における RC HPC 版の解析条件

風速
ウィンドシア
乱流強度
風車レイアウト
解析順

Table5.1-4. Simulation conditions for case D.
[m/s]
11.4
0.1
0.15
Alpha Ventus (AV10)①(5-2) 発電量評価フロー
発電量評価フローを、以下に示す。(A)及び(B)は従来の工学評価に用いられているものであり、(C)がFFB を用いたものである。Table 5.2-1 に発電量評価で用いたパラメータを示す。
(A) 工学モデル①ウェイクのない流れ場を求める。
②風速階級・風向別頻度分布に基づき、スラストカーブを用いて、速度欠損を算出する。
③風速階級・風向別頻度分布に基づき、パワーカーブを用いて、年間発電量を求める。
(B) AD/AL モデル①ウェイクを含めた流れ場を求める。
②風速階級・風向別頻度分布に基づき、パワーカーブを用いて、年間発電量を求める。
(C) FFB①ウェイクを含めた流れ場を求める。
②パワー係数を求める③風速階級・風向別頻度分布に基づき、パワー係数に発電機効率を乗じて出力を算出し、年間発電量を
求める。

32

Table 5.2-1. Parameters for AEP(Annual Energy Production) calculation.
年平均風速
[m/s]
8.5
ワイブル係数（形状係数 k）[-]
2.0
ワイブル係数（尺度係数 c）
[m/s]
9.6
定格出力
[kW]
5,000
空気密度
[kg/m3]
1.225
ロータ直径
[m]
126
2
ロータ面積
[m ]
12445
発電機効率
[%]
94.4
(5-3) FFB による単体及びタンデム解析
本研究テーマでは令和元年度までは「京」を用いて NREL5MW 風車まわりの流れ解析を実施してきた。結果を Fig.5.3-1 に示す。実線は翼素運動量理論に基づく予測結果（FAST というソフトウェアを使用）、水色は Fine グリッド（細かいメッシュ）による結果、緑は Coarse グリッド（粗いメッシュ）による結果である。図が示す通り、Fine Grid は FAST による予測結果とほほ一致しており、これまでの研究成果として FFB を用いて風車の性能を精度良く計算できることが確認されているが、羽根車 1 回転あたり、「京」800 ノードを用いて、60 時間程度の時間がかかっており、長時間の計算や条件を変更した計算を多数実行できないことが課題であった。令和 2 年度は、(1-2)節に示したチューニングされた FFB を「富岳」で動作させることにより、NREL5MW 風車の単体およびタンデムの流れ解析を実行した。Table
5.3-1 に「京」および「富岳」の計算時間の比較を示す。Table 5.3-1 に示す通り、チューニングしたコー
ドを「富岳」で動作させることにより、計算時間は 13%程度に短縮できることを確認することができた。
また、(6)章「洋上ウィンドファーム全体のスーパーシミュレーション」で述べる RC HPC 版における
風車 Wake モデルの高度化のためのリファレンスデータとして、NREL5MW 風車を風車間距離 3D（Dはブレード直径）離して直列に配置したタンデム計算を羽根車 40 回転分計算し、最後の 10 回転分の流れ場を平均したデータを作成した。Fig.5.3-2 にタンデム計算における風車 wake の流れ場の可視化結果を示す。Fig.5.3-3 および Fig.5.3-4 に風車後流のサンプリング位置および流速分布を示す。
Table 5.3-1. Comparison of calculation time of NREL wind turbine flow calculation
タンデム

「京」

「富岳」

ノード数

800

1,152

羽根車 1 回転あたりの計算時間

60 時間/rev.

7.7 時間/rev.

33

Fig.5.3-1. Comparison of thrust coefficients of NREL wind turbines

Fig.5.3-2. Mainstream velocity distribution of wind turbine wake in tandem calculation

34

Fig.5.3-3. Sampling position of wind turbine wake in tandem calculation

Fig.5.3-4. Vertical distribution of mainstream velocity time average (30-40 rev.) In wind turbine
wake in tandem calculation35

(5-4) FFB ➡ REVOCAP_Coupler ➡ ADV_Solid2 による一方向連成解析
NREL5MW 風車の単機モデルについて、流入風条件が異なる以下の 4 ケースの一方向連成解析を実施した。
(A) TSR=7 / 境界層なし / 乱れ成分なし(B) TSR=7 / 境界層あり / 乱れ成分なし(C) TSR=7 / 境界層なし / 乱れ成分あり(D) TSR=7 / 境界層あり / 乱れ成分あり(5-3)節に示した FFB による解析結果のうち、風車 1 回転分のブレード表面に相当する解析領域境界における流体力履歴をブレード表面の表面力履歴に変換し、これを境界条件として風車ブレードの応力解析を行った。流体力ベクトル→表面力ベクトルの変換には、Fig.5.4-1 に示すように REVOCAP_Couplerとファイルフォーマット変換のためのユーティリティプログラムを使用した。REVOCAP_Coupler の仕様上、FFB のプロセス数と ADV のプロセス数の合計と同じ数の MPI プロセスによる並列実行が必要であるため、REVOCAP_Coupler の実行は九大 ITO で行った。それ以外のユーティリティプログラムは、
ファイル I/O を主体とした逐次実行プログラムであるため、吉村研究室のワークステーションにて実行した。また、ADV_Solid2 による応力解析も同ワークステーションで実施した。Fig. 5.4-2、Fig.5.4-3、

FFB

及び Fig.5.4-4 にブレード先端における X、Y、Z 方向変位履歴を 4 ケースまとめたものを示した。

MESH_RCAP.P*

SURF.P*

ffb_partilst.txt

カプラ+ユーティリティ

ffbsurftrac2rcap_all
trac_*.dat

cpr.conf

cpr
trac_*.dat

coupling.dat

rcap2advtrac_all

ADV

adv_trac_*.txt
model/advhddm_in_*.adv

advsolid-p

Fig. 5.4-1. Operation Flow of conversion from fluid force to surface traction.

36━ (A) no boundary layer / no turbulence━ (B) with boundary layer / no turbulence━ (C) no boundary layer / with turbulence━ (D) with boundary layer & turbulenceFig. 5.4-2. X displacement at blade tip.

━ (A) no boundary layer / no turbulence━ (B) with boundary layer / no turbulence━ (C) no boundary layer / with turbulence━ (D) with boundary layer & turbulenceFig. 5.4-3. Y displacement at blade tip.

37

━ (A) no boundary layer / no turbulence━ (B) with boundary layer / no turbulence━ (C) no boundary layer / with turbulence━ (D) with boundary layer & turbulenceFig. 5.4-4. Z displacement at blade tip.
また、(5-3)節で NREL5MW 風車を 2 台直列に並べたタンデム解析モデルについても、単機風車モデルと同様に流体力→表面力変換を行った。応力解析にあたっては前方配置と後方配置の 2 台の風車ブレードをそれぞれ個別のモデルとして解析を行った。前方側と後方側のそれぞれの風車の位置に合わせて、
REVOCAP_Coupler のメッシュの座標変換パラメータを調整して、ブレード表面の表面力履歴を個別に変換する予定であったが、コードの不備によるデータ不整合が発生したため、対応を開始した。現在対応中である。
(5-5) ADV_Solid2 及び ADV_Fatigue による累積損傷解析
ADV_Fatigue に実装された工学モデルに基づく疲労寿命モデルを用いることで NREL5MW 風車ブレードの累積損傷解析を行った。疲労寿命モデルは主に、①４点レインフローカウントアルゴリズム、②線
形 Goodman ダイアグラム、③Miner 則の 3 つにより構成される。Fig.5.5-1 に本疲労解析のフローチャートを示す。(5-4)節で説明したように FFB による解析から得られた流体力履歴を境界条件として風車ブレードの動的構造解析を ADV_Solid2 により実施する。この構造解析から得られる、応力の時刻歴応答
が ADV_Fatigue による累積疲労解析の入力データとなる。このデータをもとに①レインフローカウント
アルゴリズムにより疲労サイクル数を求める。次に②Goodman ダイアグラムと③Miner 則によりブレー
ド全体の疲労損傷の分布を推定する。Fig.5.5-2 に②線形 Goodman ダイアグラムを示す。最後に、Paraview により累積疲労分布の可視化を行う。

38

Fig. 5.5-1. Flowchart of the fatigue life model.

Fig. 5.5-2. Linear Goodman diagram.
(5-4)節で示された、4 ケース(境界層の有無 2 ケースと乱れ成分の有無 2 ケースを考え、2×2=4 ケース)の解析により得られた応力の時刻歴応答をもとに、疲労損傷解析を行った。仮にこの荷重履歴が標準運用期間である 20 年間続いたと仮定した場合のブレードの累積疲労損傷を Fig.5.5-3 から 5.5-6 に示した。
各ケースにおいて累積疲労損傷の最大値とその発生箇所が一連の疲労解析により明らかになった。

39

Fig. 5.5-3. Fatigue life distribution in case (a) (no boundary layer/ no turbulence).

Fig. 5.5-4. Fatigue life distribution in case (b) (with boundary layer/ no turbulence).

Fig. 5.5-5. Fatigue life distribution in case (c) (with boundary layer and turbulence).

40

Fig. 5.5-6. Fatigue life distribution in case (d) (no boundary layer/ with turbulence).
また、はじめに ADV_Fatigue は Windows でのみ動作する Fortran プログラムとして開発されたが、今後 ADV_Solid2 による解析と合わせてスパコン上で実行することを考え、Linux で動作する C プログラムとして再実装を行った。
(6) 洋上ウィンドファーム全体のスーパーシミュレーション(6-1) タンデム風車の RC HPC 版と FFB の比較
一般的に、風車ブレードの回転に伴い、その後流域には Fig.6.1-1 に示されるように風速の欠損や時間・空間変動が形成される。これらの流動現象は風車ウェイクと呼ばれる。特に、複数の風車群から構成される洋上の大規模ウィンドファームにおいては、上流側に位置する風車ウェイクに起因して、下流側に位置する風車の発電量低下や風車内外の故障などが生じる可能性がある。そのため、上流側に位置する風車が形成する風車ウェイクの挙動を正しく評価することが不可欠である。風車ウェイクモデルの概要とその比較を Fig. 6.1-2 に示す。RC HPC 版では、風車ブレードの回転の影響をアクチュエータラインモデル
により再現している。アクチュエータラインモデルによるアプローチでは、風車固有のブレードに関する
いくつかのパラメータ、具体的には、風車ブレードの任意断面位置におけるコード長、ねじれ角度、抗力
係数、揚力係数が必要である。

41

(a) Visualization by wind tunnel experiment(b) Conceptual diagram
Fig. 6.1-1. Wind turbine wake.

42

Fig. 6.1-2. Overview of wind turbine wake model and its comparison.
RC HPC 版では、数値計算手法にデカルト座標系不等間隔スタガード格子に基づく(有限)差分法を用いる。LES の SGS モデルには、計算安定性に優れ、かつ壁面減衰関数を必要としない混合時間スケール
モデルを採用する。陽的フィルタ操作にはシンプソン則を適用する。速度場と圧力場のカップリングアル
ゴリズムには、陽解法を基礎とした部分段階法を採用する。圧力に関するポアソン方程式は SOR 法により緩和計算する。空間項の離散化に関して、対流項は補間法による 4 次精度中心差分を基礎とし、4 階微分の数値粘性項を付加する修正 3 次精度風上差分を用いる。ここで、数値粘性項の重みは 0.5 とし、その影響は十分に小さくした。一般的な 3 次精度風上差分の河村-桑原スキームでは 3.0 である。残りの空間
項は 2 次精度中心差分で近似した。
本計算は九州大学応用力学研究所のベクトル型スーパーコンピュータ「SX-Aurora TSUBASA」、九州
大学情報基盤研究開発センターのスーパーコンピュータシステム「ITO」、理化学研究所のスーパーコン
ピュータ「富岳」を用いた。
最初に、FFB による大規模計算と比較するために実施した RC HPC 版の結果を Fig. 6.1-3 及び Fig.
6.1-4 に示す。本計算における格子分解能は各方向においてΔx=0.5m, Δy=Δz=2.0m の等間隔とした。
各方向の格子数は Nx=1261, Ny=251, Nz=231 であり、総格子点数は約 7300 万点である。無次元時間刻
みは 5.0×10-4 であり、無次元時間 25～50 において時間平均を施した。
Fig. 6.1-3 には、瞬間場(無次元時間 50)に関してスパン中央断面における流れの可視化(side view)を示す。Fig. 6.1-4 には、Fig. 6.1-3 に対応した風車下流側から可視化した主流方向(x)速度成分の空間分
布(rear view)を示す。最初に、Fig. 6.1-3 (side view)に注目する。下流側に設置された風車は、上流側に設置された風車のブレード、ナセル、タワーが形成するウエイクの影響を強く受けていることが見て取れる。Fig. 6.1-3 (b)に示す渦度の可視化結果では、下流側に設置された風車から形成されている翼先端渦(tip vortex)に、上流側に設置された風車からの翼先端渦が干渉し、それらが合体しながら流下している様子も観察される。

43

Fig. 6.1-4 (rear view)では、Fig. 6.1-4 (b)に示す Measurement station B を含む断面内での結果が、Fig. 6.1-4 (a)に示す Measurement station A を含む断面内の結果と比較してより複雑な挙動を示していることが分かる。

(a) Spatial distribution of streamwise (x) velocity component.

(b) Spatial distribution of vorticity.
Fig. 6.1-3. Flow visualization in the central cross section of the span, Instantaneous field, Side view.

44

(a) Cross section including Measurement station A (B) Cross section including Measurement
station B
Fig. 6.1-4. Spatial distribution of streamwise (x) velocity component visualized from the downstreamside of the wind turbine, Instantaneous field, Rear view.
次に Fig. 6.1-5 には、(5)章「単体及びタンデム配置大型風車のスーパーシミュレーション」で述べた
FFB における主流方向(x)速度成分の空間分布の可視化(瞬間場)を示す。ここで、FFB の計算における無次元時間刻みは 1.57×10-4 であり、Fig. 6.1-5 に示す無次元時間は約 45 である。風車ブレードは 14.4 回
転後の結果である。FFB では「富岳」を用いた大規模計算を実施しているため、可視化結果にも非常に細かな変動が再現されている様子が伺える。Side view に注目すると、
RC HPC 版の計算結果と同様、下流側に設置された風車は、上流側に設置された風車のブレード、ナセル、タワーが形成するウェイクの影響を強く受けており、風車ウェイクが強く干渉している様子が確認される。その結果として、Rear view
においてもウェイク内の気流構造に有意な差異が確認される。
Fig. 6.1-6 には、主流方向(x)速度成分の瞬間場の比較(赤色の実線：FFB、青色の実線:RC)を示す。先に述べたように、赤色の実線で示す FFB では、非常に細かな変動が再現されている。一方、青色の実線で示す RC では、FFB のような細かな変動は再現されていないものの、FFB の傾向を良好に再現していることが分かる。
Fig. 6.1-7 には、主流方向(x)速度成分の瞬間場の比較を示す。瞬間場の比較ではあるものの、RC HPC
版の計算結果と FFB の計算結果の定性的な挙動は類似した結果を得た。

45

Fig. 6.1-5. Visualization of spatial distribution of streamwise (x) velocity component in FFB.

Fig. 6.1-6. Comparison of streamwise (x) velocity component.

46

(a)Measurement station A(b)Measurement station B

Fig. 6.1-7. Comparison of streamwise (x) velocity component.
(6-2) RC HPC 版による Alpha Ventus 解析
次に、アップウィンド型の 5MW の実機スケール風車のウエイクシミュレーションについて述べる。
ここでは、Fig.6.2-1 に示すドイツで初めてのオフショア・ウィンドファーム（大規模洋上風力発電所）
である「Alpha Ventus Offshore Wind Farm」を対象とした。このサイトには、5MW の大型風車が 12基設置されている。Fig. 6.2-1 に示すように、風車のローター直径や風車ハブ高さに若干の違いはあるものの、Fig. 6.2-1 の右図に示す風車スケールを対象にした。
最初に Fig. 6.2-2 に示すように RC HPC 版の性能等を評価するために、風車の離隔距離を実際よりも短くしたモデル（Mini Alpha Ventus モデルと称する）を作成し、これを対象に計算を行った。その結果も Fig. 6.2-2 に示す。同図に示す計算結果は渦度の空間分布を示したものであり、12 基の風車群からそれぞれ風車ウェイクが形成され、下流側ではそれらのウェイクが複雑に相互干渉している様子が明確に再現された。

47

Fig. 6.2-1. Alpha Ventus Offshore Wind Farm.

48

Fig. 6.2-2. Mini Alpha Ventus Offshore Wind Farm model construction and numerical results.
(6-3) RC HPC 版のチューニング
RC HPC 版は、将来の応用展開を見越し、「富岳」をはじめとした様々な計算機において高速に動作す
るプログラム開発を進めてきた。高速計算の核となるのは、非圧縮流解析において計算に時間がかかる圧力のポアソン方程式の反復解法部分を効率的に解く手法である Parallel Cyclic Reduction（PCR）法である[1]。PCR 法については、これまでマルチコアアーキテクチャである Intel Xeon Skylake-SP、「富岳」と同様なメニーコアアーキテクチャである SGI UV300、キャッシュ付きのベクトル型プロセッサである
NEC SX-Aurora TSUBASA を用いたベンチマークを実施し、いずれのアーキテクチャにおいても高い実行性能を達成できるアルゴリズムであることを確認した。令和 2 年度は、PCR 法のさらなる高速化、「富岳」でのチューニングと実行性能の確認を行った。
PCR 法は、三重対角行列の解法で最も効率的である LU 分解を元にした Thomas Algorithm を並列処理する解法である。Thomas Algorithm は逐次アルゴリズムであり、高速計算のために必須な SIMD 演
算も利かない。その原因はループの依存性であり、PCR 法ではその依存性を除去するために、Fig. 6.3-1に示すような「縮約」操作を繰り返し行う。１回の「縮約」により、与えられた１組の N 元連立一次方程式は依存関係のない２組の N/2 元連立一次方程式に分類される。この操作を繰り返すと、１組→２組→４組→・・・→N 組の１個の方程式となる。最終段階では方程式は自明に解けることになる。

49

Fig. 6.3-1. Reduction operation in Parallel Cyclic Reduction method.
Fig. 6.3-2 に PCR 法のカーネルコードを示す。関数への引数ｓは縮約の回数を示しており、縮約が進むと参照するメモリの距離が離れることがわかる。考えている連立一次方程式の数は Z 軸の格子点数な
ので、せいぜい 512 程度である。単精度の場合、d, a, c, dn, an, cn の６配列を利用する実装となっているので、計算に必要な記憶領域は 4Byte × 6 × 512 = 12kB 程度である。テストした計算機の L1 キャ
ッシュサイズは 32kB 程度であるため、計算に必要なデータはほぼ L1 キャッシュに載る。したがって、
参照のストライドは大きくなるが、全て L1 キャッシュ内でアクセス可能となる。実際に Skylake-SP での実行時のプロファイル結果からは、Fig. 6.3-3 のように PCR 法は他の古典的な方法に比べて圧倒的に
L1 キャッシュのヒット率が高いことが確認でき、アルゴリズムの有効性の根拠を示している。
次に、PCR 法をさらに高速化するため、PCR の縮約回数を最終段から２段手前で止める実装を評価した。これは、縮約を繰り返すことによる精度低下を抑制することと、２段手前で得られる 4×4 の方程式
を直接反転することにより B/F 値をさらに低下させるねらいがある。Fig. 6.3-4 にその効果を示す。ここではアウタースイープを昇順に走査する pcr_esa と Red-Black オーダリングで走査する pcr_rb_esa を比較した。Original と 4×4 を Cramer の方法で直接反転する実装を比較している。Cramer の方法はいずれもオリジナルよりも性能が改善されていることがわかる。

50

Fig. 6.3-2. Pseudo code of PCR.

Fig. 6.3-3. Cache hit rate for 36 cores on Skylake-SP that has 32kB L1 cache.

51

Fig. 6.3-4. Performance improvement by direct inversion of 4x4 matrices.
次に、「富岳」でのプロファイリングに基づくチューニングを実施した。実施した内容は、不連続なメ
モリアクセスが生じる箇所でプリフェッチを指示と、ソフトウェアパイプライニング最適化の２点である。測定結果を Fig. 6.3-5 に示す。演算性能の測定は、「富岳」のプロファイラと PMlib による手動演算
数カウントによる方法を載せている。プロファイラは投機実行などによる余計な演算器の動作も含んでいるため、プログラムループ中の演算数をカウントした測定よりも３割ほど多くなる点に注意を要する。
１CMG（12 スレッド）の場合も４CMG（48 スレッド）の場合もオリジナルよりも 1.23 倍以上の演算性能であることがわかる。メモリスループットは 50GB/s 程度でバンド幅を使い切っていないが、本アルゴリズムは L1 律速であるため、余裕がある。チューニング後の 48 コア実行時のハードウェアカウンタの値は推定１TFLOIPS となり、ピーク性能の 29％となる。
以上のように、提案する PCR 法は L1D キャッシュ負荷が大きく、低 B/F 値、動作周波数の高い CPU
向く、キャッシュを有効利用できる、メニーコアの並列性を引き出せるアルゴリズムであり、「富岳」においても高い性能を発揮することが確認できた。

52

Fig. 6.3-5. Performance results on an A64FX.
本提案手法は、圧力のポアソン方程式だけでなく、時間積分の陰解法に現れる連立一次方程式にも適用することができ、今後 RC HPC 版に組み込んでいく予定である。
文献[1] Kenji Ono, Toshihiro Kato, Satoshi Ohshima, Takeshi Nanri, “Scalable Direct-Iterative Hybrid
Solver for Sparse Matrices on Multi-Core and Vector Architectures,” HPCAsia2020: Proceedings ofthe International Conference on High Performance Computing in Asia-Pacific Region, pp.11–21,
https://doi.org/10.1145/3368474.3368484, 2020.
(7) 個別アプリケーションと AI の連携活用(7-1) ADVENTURE
ADV_Thermal が生成する大規模数値計算結果を用いて AI（以下、深層学習）により予測器を生成[1,2]することを目的としている。予測器は、以下の目的のために利用する。
（目的１）大規模数値シミュレーションを行わないでも傾向を予測する（目的２）超高速予測により解の存在パラメータ範囲を推定する
つまり、大規模数値シミュレーションを実施する前の予備検討を行う簡易シミュレーターを構築することが目標となる。数値シミュレーションの結果を用いて予測器を生成するためにはこれまでのデータサイエンス分野で前提[3]としている、データ数が一様に膨大である、予測対象が事前にモデル化可能である、といった仮説（前提）が一般論として成立ない。物理現象を予測するためのデータは次の事項を前提とするのが工学的な観点で合理性があり本研究の調査対象でもある。
（仮説１）計算結果は数十から数千程度しか存在しない（極めて限られた数の学習データ）（仮説２）予測範囲は事前に決定可能である

53

以下、ADV_Thermal との連携を考慮した発熱体の位置条件（境界条件の変化）により測定点（固定）位置の予測が可能かどうかを調べた。この調査にあたっては次の点に留意した。
（留意点１）仮説１から過学習を起こさないための正則化技術[4]をどこまで適用するかの調査（留意点２）パラメータ空間の設定と精度の悪い箇所の関係をある程度理解[5]（留意点３）出現頻度が低いデータは予測精度が悪い（AI 及び機械学習の特徴）[6]（留意点４）１水準のデータでは学習ごとに予測結果が激しく変化する（線形では２水準、非線形では最低３水準必要）[7]
これらの考えの有効性を確認するために、ケース A は空間に配置された発熱体の位置のみを変化させて学習した。ケース B では最も影響度の高い発熱体 HE1 の発熱量を３水準として学習データを増加させた。また、正則化は様々な方法があるが、シミュレーションのデータの場合には、画像認識や分類における正則化技術がそのままでは適用できない。そのため本検討ではデータ拡張を実施した。基本的にデー
タ拡張はオリジナルのデータの入力側に数％の誤差を混入させる。誤差評価関数は二乗和誤差を用いる
ため、誤差を追加した結果その誤差の和が０となるように誤差を付加する。この点を守らないと異なる値が正解であるかのように学習する。
(7-1-1) 解析対象ケース A の解析条件
発熱体を Fig. 7.1-1 のように配置した。HE1 と 3 は固定とし、HE2 と 4 は破線内を動かし温度分布が変化する因子とした。この理由は、汎用プログラムと同等の任意な対象を予測するのではなく特定の現象の特定の対象に対してコストを低くするための必要最低限の条件を設定することが重要であるためである。ここでは、発熱体の位置のみに着目することにした。必要に応じて、発熱体の発熱量も因子にすることは可能であるが、せいぜい 3 水準程度でよいと考える。この理由は、機械学習が速やかに学習を終了するために必要なためである。

×(12)

×(7)

×(5)
各発熱体発熱量
HE1: 10.0 W
HE2: 1.5 W
HE3: 1.5 W
HE4: 1.5 W

HE3

×(1)

×(13)

×(11)

×(10)

×(3)
HE1

×(9)HE2
×(4)

×(2)
HE4

×(8)

×(6)

Fig. 7.1-1. Positions of heating elements：HE1 and HE3 are fixed. HE2 and HE4 are moved withindashed squares for making variations of thermal distribution. Cross mark, which is ×, of #1 to 8
are for training and Cross marks of 9 to 13 are for validation.

54

(7-1-2) 解析対象ケース A の畳み込みニューラルネットワークの入力データ設計
Fig. 7.1-2 に入力データマトリックスと畳み込みニューラルネットワーク（CNN）の構成を示す。正則化のための様々な技術は CNN において実証されてきた。
そのため、この研究においても回帰問題を CNNで予測する。画像認識と異なる点は、入力となるマトリックスが画像に比べて小さく、また、色情報ではなく物理量や隣接関係を表すトポロジーなどになる点である。データを準備する段階で行うべきでない
ことは、二値しかとらない値（たとえば男性と女性、あるとない、などの二分類されるデータ）は入力として扱わないことである。これらの分類は意味のある連続量として置き換える必要がある。理由は、回帰問題であるのに二値入力があった場合に学習させるためのコストが膨大になるためである。これらを考慮して決めた入力マトリックスを Fig. 7.1-2(a)に示す。活性化関数の特性を考慮して空欄には 0 を入れる。解析対象ケース A では、次のようにデータを生成した。
・発熱体配置：7P2=42 通り（HE1,HE3 は固定）
・発熱量：1 通り
・学習点：8（×の 1～8）・予測点：5（×の 9～13）・出力：温度 T（熱伝導解析）7P2×1×8＝336 ケースのオリジナルデータ
・データ拡張：10 倍を実行 合計 3696 ケース

座標位置データ（発熱量[W]）
HE3
IC3
HE1
IC1

HE4
IC4

HE2
IC2

熱源座標値
x1

y1

x2

y2

x3

y3

x4

y4

X

dx1

dy1

dx3

dy3

Y

dx2

dy2

dx4

dy4

学習 or 予測点

X,Yとxi,yiとの距離(a) Input data matrix: 0 is assigned at all of blanks. 0 reproduce 0 in activation functions of ReLU
and tanh.

axbx1

axbxc1

axbxc2

axbxc3

Zero padding / (m,n)filter
ResNetd=axbxc3

e 1: tt: temperature at apoint(b) Configuration of convolutional neural network: a=b=6, c1=c2=c3=500, e=240
Fig. 7.1-2. Convolutional neural network configuration and hyper parameters and input parameters
design for training case A.

55

(7-1-3) 解析対象ケース A の学習結果
Fig.7.1-3 と 7-1-4 に学習結果を示す。Fig.7.1-3 の見方は、検証データを×の位置の順番に並べた。最初の予測温度が高い部分は×9、その次の高い部分が×10 となり、最も遠方である×13 が約 170 番目以降の分布になる。基本的には入力に対してグラフ化することが困難なため、正解を青実線としてその参照
解と一致すれば精度の高い予測ができたと判断する。この結果からわかることは、×10 を除く全てにおいて予測結果が悪く、特に、低温部分の予測は高めの予測結果となっている。Fig. 7.1-4 により、10％以上の誤差がある予測値が 45％以上あることがわかる。このような結果になった理由として、×9 と×10
は空間パラメータとして学習した結果のため補間領域となっていることが考えられる。一方で×11 以降
は完全に補外領域で今回想定した学習結果からは微分量まで学習できずに補外領域に入った後は活性化関数が飽和する値をほぼそのまま出力していると理解できる。微分量を学習させるためには、発熱のパターンを増やすことが１つの方法として有効である可能性がある。したがって、解析対象ケース B では最も影響度の高い HE１の発熱量を変化させることが有効である可能性があると考察する。

✕９の

✕10 の

✕11 の

✕12 の

✕13 の

予測結果

予測結果

予測結果

予測結果

予測結果

と正解

と正解

と正解

と正解

と正解

Fig. 7.1-3. Evaluation result of case A: Temperature is normalized between 0.0 to 1.0. Blue solid line
is reference temperature. Red solid line and dots are predicted temperature. Under 0.7, whichmeans evaluation positions are far from heating elements, predicted t doesn’t agree with reference t.

56

Fig. 7.1-4. Frequency and percentage of error ranges: over 10% error range accounts for 45%.
(7-1-4) 解析対象ケース B の解析条件
Fig. 7.1-5 に解析対象ケース B の発熱体配置及び測定点、発熱体の発熱量を示す。解析対象ケース A との違いは HE1 の発熱量を 5W、10W、15W の 3 ケースとした点である。336 ケースの 3 倍のケース数となり、1008 ケースがオリジナルのデータ数となる。

×(12)

×(7)

×(5)
各発熱体発熱量
HE1: 5,10,15 W
HE2: 1.5 W
HE3: 1.5 W
HE4: 1.5 W

HE3

×(1)

×(13)

×(11)

×(10)

×(3)
HE1

×(9)HE2
×(4)

×(2)
HE4

×(8)

×(6)

Fig. 7.1-5. Positions of heating elements：HE1 has 3 levels of heating. HE2 and HE4 are moved
within dashed squares for making variations of thermal distribution. Cross mark, which is ×, of #1
to 8 are for training and Cross marks of 9 to 13 are for validation.

57

(7-1-5) 解析対象ケース B の畳み込みニューラルネットワークの入力データ設計
Fig. 7.1-6 に入力データマトリックスと CNN の構成を示す。解析対象ケース A とは入力データマトリックスが同じで、CNN の構成を Fig. 7.1-6(b)のように変更した。変更は学習が最も進むように、ハイパーパラメータである Fig. 7.1-6(b)の変数を、c1=560, c2=280, c3=140, e=160 とした。
(7-1-6) 解析対象ケース B の学習結果
Fig. 7.1-7 と 7.1-8 に学習結果を示す。HE1 の発熱量を３ケースとして学習データを増やすことにより、大幅な精度の向上が見られた。基本的には補外領域で×13 に相当する部位の予測結果（Fig. 7.1-7 の横軸 500 以降の予測値）は発熱量が低いケースでは良好な一致を示すが、それ以外では徐々に予測値の誤差が増えることがわかる。なぜ、発熱量のケース数を増やすことが温度場の予測精度の向上につながるのかについては、より多くのパラメータにより補間能力の向上により微分量の精度も向上し、微分量がネットワーク内部に陰に構成されたと理解する。このことより、ある程度の補外領域における予測が可能になるとの仮説[4]を間接的に説明できると考えている。

座標位置データ（発熱量[W]）
HE3
IC3
HE1
IC1

HE4
IC4

HE2
IC2

熱源座標値
x1

y1

x2

y2

x3

y3

x4

y4

X

dx1

dy1

dx3

dy3

Y

dx2

dy2

dx4

dy4

学習 or 予測点

X,Yとxi,yiとの距離(a) Input data matrix: 0 is assigned at all of blanks. 0 reproduce 0 in activation functions of ReLU
and tanh.

axbx1

axbxc1

axbxc2

axbxc3

Zero padding / (m,n)filter
ResNetd=axbxc3

e 1: tt: temperature at apoint(b) Configuration of convolutional neural network: a=b=6, c1=560, c2=280, c3=140, e=160
Fig. 7.1-6. Convolutional neural network configuration and hyper parameters and input parameters
design for training case B.

58

Fig. 7.1-7. Evaluation result of case B: Temperature is normalized between 0.0 to 1.0. Blue solid line
is reference temperature. Red solid line and dots are predicted temperature. 3 levels reinforce
ability of interpolation in the far area from heating elements.

Fig. 7.1-8. Frequency and percentage of error ranges:
over 10% error range accounts for less than 2%.

59

(7-1-7) まとめ
本検討により以下の結論を得た。
・オリジナルの発熱量のみではデータ量不足で過学習が生じる・パターンを与えると精度が大幅に上昇する
・以前より使用している入力設計で学習・予測が可能
・今回の対象とする温度場の予測では、1000 ケース程度は必要であり、ある程度温度のバリエーションも必要である・データ的に補外に該当する部分が原因で学習不良が生じる
・データ量の単純な増加でも傾向を表現することは可能であることを示した
・より精度を出すには精度の悪い予測に対して、出現頻度などの統計量を詳細な確認が必要である
今後は、これらの知見を用いて、ある程度のモデルを前提[8]とした手法を検討し、物理法則をもつ意
味のあるデータ量[9,10]を推定し、実機への適用を目指す。
文献[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] 機械工学年鑑 2020, https://www.jsme.or.jp/kikainenkan2020/chap02/#a08（参照日 2021 年 4 月
1 日）[9] A. Choromanska, et al., The loss surfaces of multilayer networks, arXiv:1412.0233 [cs.LG], 13
pages, 2015[10] G. Montúfar, et al., On the Number of Linear Regions of Deep Neural Networks, Twenty-eighth
Conference on Neural Information Processing Systems, 9 pages, 2014(7-2) FFR-Comb(4)章「パイロットスケール超臨界圧 CO2 ガスタービン（東芝エネルギーシステムズ製）スーパーシ
ミュレーション」を参照のこと。

60

(7-3) FFB
洋上風車のような大規模な装置の流体解析は「富岳」を用いて可能になったとはいえ、計算負荷が大き
いこと、風車を運用する事業者にとって負担であることには変わりがない。そこで大規模解析で得られたデータをもとに風車周りの流れを推定する手法の開発が必要となる。令和 2 年度は、これまでの解析で得られた風車後流のデータベースを用いて、風車下流の流れ場を推定することを試みる。大規模解析により得られた多数のデータをもとにサロゲーションモデルを作製するには、解析データを機械学習によって推定する方法と簡便化された数学モデルを用いる方法が考えられるが、今回はまずは物理的な背景に基づく数学モデルを作製し、モデルにおけるパラメータを大規模解析データにより調整することとした。
風車周りの流れ場は Fig.7.3-1 に示すような流管モデルとしてモデル化される。

Actuator disc

z

A0

A

y
x

T

Aw

p0 ,U 0pw , U w
Fig.7.3-1. Control volume of horizontal axis wind turbine rotor modelled by an actuator disc.
このとき、風車下流の速度欠損は

U= C ( ( x)) f (r ,  ( x))U0と表すことができる。ここで関数 f を

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

２つのガウス関数の和として表すこととする。
このとき風車に作用する推進力は
r =

T = U 02C ( ( x))  (e DG +e DG −
r =0

r =

+

−

+
+
−
C ( ( x)) 2 DG+(e+ e2 DG + 2e DG + DG ))rdr2

M =  (e DG + e DG )rdr
+

−

r =0

r =

(

)

+
−
+
−
1

N =   e2 DG + e2 DG + 2e DG + DG  rdr2

r =0 

61

T = U02C( ( x))(M − C( ( x)) N ) .

M 2 − 1/ 2 NCT D 2
2N.

C  ( ( x)) = M と表すことができる。
− r02

M = 2 ( x)e2

N =  ( x )e2

2 2 ( x )− r02
 2 ( x)+
r0
+ 2 r0 ( x) erf  2 ( x)  , r r0 ( x)erf  0 2  ( x)  .

2


2

  
2
z
−
z
(
)
r
1
y


h
0

  
 exp  − 2   2 ( x, y ) +  2 ( x, z ) +  ( x, y, z )   yz  
1
f ani (r ,  ( x, y, z )) = 
2
2



2


( z − zh ) − r0   y2
 + exp  − 1 + 2   2 ( x, y )  2 ( x, z )  ( x, y, z )   

y
z



  


さらに風車下流の乱流の非等方度を考慮すると

  2 − r02
D2 2
 r0  
1 − 1 − CT  =  2 e 2 + 2 r0 erf 


8 

 2  
 
2

1/2

2
2


−r2
−r2
 2 − r02

 r0     2 202
 r0   1  2  20
 r 
2
 2 e + 2 r0 erf  −  2 e + 2 r0 erf  −  e +r0 erf  0   CT D 2 




2
    2    
 2   2 




2
 2 − r20r 2 e  +r0 erf  0  2   
.

が得られる。これらの式に現れる未知項を数値解析データを用いて学習させ、サロゲートモデルとして利用する。解析結果等を用いてパラメータを推定した一例を下記に示す。

kr = 0.5405; sout = 3.6077; k y = 0.0133; k z = 0.0123;  y = 0.2335;  z = 0.2119
開発したモデルにより風車後流を推定した結果を Fig.7.3-2 に示す。
提案した手法により風車近傍(x/D=1.4)から遠方後流（x/D=9）まで精度良く近似できることがわかる。
従来もこの様なモデルは提案されていたが、今回はガウス関数を二つ組み合わせることにより風車近傍と遠方後流を同時に解析することが可能となった。また、パラメータの同定には上流側については数値解
析データ、下流側については実験を参照することにより、広い範囲にわたってフィッティング可能なパラメータを設定することが可能となった。
次に風車上流の境界層を考慮した計算事例を示す。Fig.7.3-3 に風車上流の速度分布、Fig.7.3-4 に風車

62

下流の速度場を示す。境界層を考慮した鉛直方向分布についても推定することが可能となった。
Fig.7.3-5、6 にサロゲートモデルによって推定した風車下流の速度分布を示す。現在のモデルは上流流
れの乱流の影響、後流の蛇行に関するモデルが適用されていないため、これらの現象をモデル化するための検討を行う必要がある。また、風向変化の影響についてもモデル化を検討する。

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

63

Fig.7.3-3. Vertical profile of incoming normalized ABL wind velocity.

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

65

Fig.7.3-5. The contour of cross-sectional wake velocity at several downstream distances from theturbine.

66

Fig.7.3-6. The contour of vertical wake velocity at several spanwise positions.
(7-4) RC HPC 版
洋上風車の解析では、信頼性の高い予測結果を得るために、様々な風向の境界条件や諸条件を考慮して総合的な観点から設置検討を実施する必要がある。このため、多数のシミュレーションが必要になるが、高精度化のために計算規模を大きくすると限られた計算資源や時間の中で結果を得ることができない。
また、将来産業界への展開を考えた場合、より小さな規模で精度のよい計算ができる方が望ましい。洋上67

風車の風車配置の問題については、後流（ウェイク）干渉を精度良く再現することが鍵である。令和 2 年
度は、高精度の計算と小規模計算は成立しにくいトレードオフの関係にあるが、少数の大規模計算から粗い格子でも細かい格子の高精度な後流干渉の状況を再現するサロゲートモデルを開発することにより、
洋上風車シミュレーションの活用シナリオを実現するアプローチについて検討した。
サロゲートモデルの開発には、風車の作用により後流を再現する後流モデルを基本として、アンサンブ
ル計算を想定する粗い格子の格子幅パラメータを導入したモデル探索を行う。探索方法として、少数の高精度なシミュレーション結果を参照解として、その解を粗い格子で再現するモデル方程式を遺伝的アルゴリズムを用いた方程式探索を行う。探索する項は、元のナビエストークス方程式には現れない付加的な項であり、その項が細かい格子と粗い格子の差を表現･再現する。
方程式発見の手法については、小野らが開発した微分項を考慮した方程式発見の方法[1]を用いる。
文献[1] K.Ono, and I. Koga, “Rediscovery of governing equations from simulation data using Genetic
Programming,”Transactions of JSCES, Paper No.20201004, 2020.
(8) 産業界コンソーシアムや外部機関等との連携(8-1) 燃焼系産業界コンソについて
燃焼システム用次世代 CAE コンソーシアムの第 5 回総会（令和 2 年 9 月 24 日）においてプロジェク
ト代表者の吉村が基調講演を行った。本プロジェクトで開発中のスーパーシミュレーション技術の紹介を行うとともに、ディスカッションを通して燃焼系産業界のニーズの聴収を行い、その結果を精査し、本シミュレーション技術の開発に反映させた。
(8-2) 風力産業界コンソについて(8-2-1) 国内風力産業コンソーシアムの検討
本プロジェクトでめざす洋上風力コンソーシアム構築にあたり、協力を要請する団体やコンソーシア
ムで考慮すべき問題等について検討するため、既存の団体について調査し、それぞれの特徴や参考とすべ
き点をまとめた。
(8-2-1-1) 日本風力発電協会（JWPA）
日本風力発電協会（以下、JWPA と呼ぶ）は、日本の風力発電産業団体で、設立から 10 年経ち、その会員企業は 400 社を超える。JWPA では、中・長期導入目標実現に向けた提言、最新情報の収集・発信、業界の課題を抽出・検討する専門部会の設置、国内外関係団体との交流・協調、普及・啓発活動等を行っている。
専門部会は、政策部会、環境部会、技術部会、系統部会、国際部会及び広報部会の 6 つの部会から成る。技術部会が技術的な調査研究を担っており、その下に洋上技術調査 WG、風況・検討 WG が設置され、会員企業から、洋上風力や風況の専門家、担当者が参加をしている。
JWPA は、風力産業コンソーシアム組成にあたり、相談をする先として、最も有力な産業団体の一つと考えられる。不定期に技術セミナーや技術紹介等の開催を実施しており、会員外からの技術紹介の場を68

設けることが可能である。
(8-2-1-2)日本風力エネルギー学会（JWEA）
日本風力エネルギー学会（以下、JWEA と呼ぶ）は、日本の風力エネルギー学術団体で、設立から 44
年経ち、団体会員 14 社、個人会員等 400 名を超える。JWEA では、風力エネルギー利用に関する基礎と応用について、科学技術の振興と普及を目的とし、会員相互の連絡・親睦、国内外の研究者、研究団体との交流を行っている。
年 4 回の機関誌「風力エネルギー」の発行、年 1 回のシンポジウムの開催を行っている。
また、JWEA では研究会の開催をしており、これまでに「空中風力発電研究会」、「ブレード技術研究会」が実施された。2021 年からは「風車ウェイク技術研究会」、「風車運用寿命延長技術研究会」、「垂直
軸風車空力弾性解析技術研究会」が新たに発足する。
JWEA の会員は、産学から成り、企業の研究開発部門に所属する会員等は、スーパーシミュレーショ
ンの直接の利用者として、コンソーシアムメンバーの候補となりうると考えられる。また、スーパーシミュレーションに期待する現場のニーズを聞く場を設けることも可能である。
(8-2-1-3) IoS-OP プラットフォーム
IoS-OP（Internet of Ships Open Platform）は船舶の運行データを、データ提供者の利益を損なわずに、ステークホルダー間での共有や、造船所やメーカー等への利用権販売、各種サービスへの提供を可能とすべく、海事業界内で合意されたルールと、データセンターで構成された共通基盤であり、その運用を
日本海事協会の子会社である株式会社シップデータセンターが担っている。
データの流通に関わる部分を協調領域とし、データを活用したイノベーション、新規サービスの開発な
どの競争領域に注力できることを目指し、業務改善を目指して自らが強くなり、新しいビジネスモデルを探求することで、業界全体としてデジタル時代における新たな海事クラスターの形をつくり、次世代につなぐことを目的として活動している。Fig.8.2-1 に IoS-OP コンソーシアムテストベッドサービス活用シ
ナリオを示す。
船舶 IoT データの基盤を、洋上風力向けに拡張させることで、効率よく、短期間で洋上風力業界におけるビッグデータの活用の可能性があり得ると考えられる。

Fig.8.2-1. IoS-OP Consortium testbed service utilization scenario.

69

(8-2-1-4) JIP（Joint Industry Project）
欧州では洋上風力発電の大幅な低コスト化を実現してきたが、そこで実施された技術開発手法の例として、英国の CarbonTrust による市場プル型アプローチを採用した Offshore Wind Accelerator（OWA）
や Floating Wind Joint Industry Project （FLOJIP）が挙げられる。JIP（Joint Industry Program、
産業界共同技術開発プログラム）とは、技術を必要とする顧客（洋上風力発電事業者）が集まり、政府機関と共に技術開発資金を分担し、技術開発者（大学、研究所、メーカー、ゼネコン、マリコン）と密接に協力関係を構築することで、市場のニーズを反映した、レバレッジ効果の高い技術開発を実現するプログラムである。
技術開発の対象となるのは、技術開発水準／技術成熟度（TRL：Technology Readiness Level ）が 4（実証試験の実施）以上の技術とされており、商用利用を対象とした迅速な技術開発の実現を目指している。技術プル型アプローチによる開発に比べ、JIP を採用した市場プル型アプローチによる技術開発は、
「市場性」
、
「促進化」
、
「サプライチェーン」、「習熟」の四つの観点から、Fig.8.2-2 に示す特徴をもつ。

• イノベーションの迅速な商用化• 最もコスト効率のよい解
決方法を選択• 産業界は市場のニーズを踏まえた開発が可能• サプライチェーン確立及び強化促進• サプライチェーンへの投資増加• 産業界における事業開発
の加速

市場プル

促進

サプライ
チェーン

習熟

• リスクとファイナンスコストの低減• 風力産業の競争力強化

Fig.8.2-2. Features of JIP.
この JIP が成功するカギは、JIP 参加者間の関係構築（ステークホルダー・エンゲージメント）にあり、参加者が多いほどレバレッジ効果の高い技術開発となる。Fig. 8.2-3 は、JIP 参加のメリットを示す。

70

Fig. 8.2-3. Benefits of participating in JIP.
市場プル型アプローチ

利点

・発電事業者主体の組織

・コスト目標が明確

・汎用性・普遍性のある技術開発要素の抽出

・開発時から認証機関が関与・共通のモデルウインドファーム及びシナリオに基づいた研究開発することで技術の普及が速まる。

・マネージメントチームが準備する技術情報をベースに議論・一般化した情報のため評価しやすい。

・情報は参加者のみで共有
・知財は技術開発者（メーカー・マリコンなど）が取得
高いレバレッジ効果

利点

・発電事業者及び国による共同出資（国＋参加者の・発電事業者・開発者双方で案分）

開発リスクを軽減できる。

・発電事業者のニーズ（技術・コスト）とサプライチェーンのマッチングが高まる。
高い習熟効果

利点・発電事業者のエンジニアリングスキルが向上する。

・オーナーズエンジニアリン
グのスキルが向上する。

・発電事業者がリスクや利点をこれまで以上に判断できるようになる。
これまで、OWA の中では、ウェイクモデルの開発、ウィンドファーム制御、ブロッケージ効果などの研究テーマに取り組んできている。発電事業者がプロジェクトに直接参加することにより、検証サイトやデータの提供等が円滑に行われている。Fig. 8.2-4 にウェイクモデルチャレンジのベンチマーク結果の例を示す。ここでは、隣接するウィンドファーム間のウェイク予測に関して、従来のウェイクモデルでは過小評価をするという課題が報告されている。

71

Fig. 8.2-4. Wake model challenging benchmark results.
一方、国内では、現在、NEDO「次世代浮体式洋上風力発電システム実証研究（浮体式洋上風力発電低コスト化技術開発）」において、浮体式洋上風力発電の FS 調査の実施にあたり、JIP 方式の技術開発が検討されており、市場プル型の技術開発の国内展開に関する上での課題の抽出などが行われている。
スーパーシミュレーションシステムの実用レベルの普及にあたり、JIP 方式を参考にすることによって、実用化の促進、迅速化に繋がる可能性があると考える。
(8-2-2) 洋上風力コンソーシアムの立ち上げに向けた取り組み
次世代 CAE を活用する風力産業コンソーシアム（仮称）を立ち上げるため、令和 3 年 3 月 31 日にオンラインにて洋上風力ワークショップを開催し、コンソーシアム立ち上げの説明と意見交換を行った。その際に、風力業界関係者を含むワークショップ参加者に本プロジェクトで開発中のスーパーシミュレー
ション技術の紹介を行うとともに、風力産業界のニーズの聴収を行った。その結果を精査し、本シミュレーション技術の開発に反映させ、令和 3 年度上半期中にコンソーシアムの準備会合を開催するめどを立てた。国内外の企業・大学、新聞社からも参加があり、洋上風力発電の技術やコンソーシアムへの関心の高さが伺われた。
(8-3) IEA Wind Task31 へ参画準備
IEA Wind Task31 への参画準備として、Task31 の現在の活動状況を調査した。
(8-3-1) IEA Wind の概要
IEA は、1974 年に OECD の枠内に設立され、2020 年末時点で日本をはじめとする 37 カ国が加盟し、「四つの E」（エネルギー安全保障の確保（Energy Security）、経済成長（Economic Development）、環境保護（Environmental Awareness）、世界的なエンゲージメント（Engagement Worldwide））を目標に掲げ活動をしている。

72

IEA の下に 38 の技術協力プログラムがあり、IEA Wind TCP は、そのうちの一つになる。各国の最新の研究開発、政策動向等についての情報収集が可能となるとともに、IEA Wind の成果（技術推奨基準）
が IEC の国際規格に発展するケースが多くある。二つの国際活動は、各々独立した体制・組織により行われているが、技術的には相互補完的関係にある。
IEA Wind の活動方針は、年 2 回開催される執行委員会（Executive Committee：ExCo）において審議される。国際協力活動は、研究開発テーマごとに設置される Task を通じて行われ、2021 年 3 月時点
で 18 の Task が活動している（Table 8.3-1 及び Fig.8.3-1）。日本は、そのうちの 9 の Task に参加している。近年、分散型エネルギー、寿命延長評価及びデジタリゼーションが 2019 年から、風車翼の腐食が2021 年から新規タスクとして立ち上げられた。
国内では、国内の IEA Wind に関する意思決定機関として国内委員会が設置され、各 Task へ派遣する
有識者、参加する Task、分科会の設置等を審議（年間 3 回程度実施）している。
Table 8.3-1. List of IEA Wind Task.
Task 11

基礎技術情報交換

Task 19

寒冷気候における風力発電

Task 25

産総研

風力大量導入時の電力システムの設計

京大、東理大、電中研、東大、横国大と運用

Task 26

風力発電のコスト

東大

Task 28

風力発電事業の社会受容性

名大、名経大、法政大、ISEP

Task 29
Task 30
Task 31

風洞実験計測データ解析と風車空力モ
デルの改良
洋上風車動的解析コード比較・検証

東大、日本海事協会、RCCM

ウィンドファーム流れモデルのベンチ

東大、WINC、豊技科大、みずほ情報総研

マーク

Task 32

風計測ライダーの風力発電への導入

Task 34

風力発電における環境影響の低減策

Task 36

風力発電予測Task 37

風力発電システムの全体設計

Task 39

低騒音風車

Task 40

ダウンウィンド風車技術

Task 41

分散型エネルギー

Task 42

寿命延長評価Task 43

デジタリゼーション

Task 46

風車翼の腐食

三菱電機、CTC、産総研、GPI

九大、WEIT、日立、産総研、東大

産総研（2021 年 3 月現在、太字が日本参加 Task。右列が参加機関。
）

73

Fig. 8.3-1. Strategic priority areas and tasks 2015-2020.
(8-3-2) 専門家会議
専門家会議（Topical Expert Meeting：TEM。以下、TEM と呼ぶ）は、専門家間の技術情報交換を通じて風車技術を振興することを目的とした会議で、Task11 の活動として行われる。Table 8-3-2 に最近の
TEM のトピックを示す。TEM として議論された研究開発テーマが新たな Task として発展するケースが多くある。近年の TEM の開催国を見ると、デンマークとアメリカがほとんどで、この 2 国が風力の技
術開発のイニシアティブを取っている状況といえる。

74

Table 8.3-2. Recent TEM topics (bold letters are new tasks).
年

No.
102
101

2020

100

99

98
97
96
2019

95

94
93
92
2018

91

90

トピック

開催日

Airborne Wind Energy（空中風車）
HybridPower

PlantsChallenges

andOpportunities（ハイブリッド発電）
Aviation System Cohabitation（航空システムとの両立）
Floating Offshore Wind Arrays（浮体式洋上風
車配置）
Erosion of Wind Turbine blades（風車翼の腐食）
Wind Farm Controls（ウィンドファーム制御）
Wind Plant Decommissioning, Repowering andRecycling（風力発電所の撤去、リプレース、リサイクル）
Improving the Reliability and Availabilityof Electrical Infrastructure Components（電気設備部品の信頼性及び稼働率向上）
Large Component Testing for Ultra-long Wind
Turbine blades（超長尺風車翼の大規模要素試験）2020.9.23,24

オンライン

2020.8.24-26

オンライン

2020.12.8,9

オンライン

2020.7.15,17,20

オンライン

2020.2.6,7

デンマーク

2019.9.25

オランダ

2019.11.14,15

イタリア

2019.4.8,9

英国

2019.2.25,26

米国

Wind Turbine Life Extension（風車寿命延長） 2018.12.13
Wind Energy Digitalization（風力エネルギーのデジタリゼーション）
Durability and Damage Tolerant Design ofWind Turbine Blades（風車翼の耐久性及び耐損
傷性設計）
Community and Distributed Wind（分散風力コ
ミュニティ）

開催国

デンマーク

2018.10.4,5

米国

2018.6.12-14

米国

2018.3.26-28

デンマーク※本プロジェクトに関連するトピックについて、青字・水色塗潰しで示す。
Table 8.3-2 において、本プロジェクトと関連するトピックを青字・水色塗潰しで示した。近年、風車の大型化及びライフタイムを通した発電コストの低減に関する技術開発が重要となっており、これらはスーパーシミュレーションが活用される可能性の高い技術分野であると言えるが、現在、日本から
Task42 寿命延長及び Task43 デジタリゼーションには参加していない。なお、Wind Europe（欧州風力エネルギー協会）の定義では、風車基礎を変える場合をリパワリング、変えない場合を寿命延長と呼ぶ。
欧州と比べると、風速が低く、乱流強度が高い日本において、高効率、低発電コストの実現をするためには、寿命予測、寿命延長、風車最適設計、最適風車配置及び最適制御を行うための高精度な疲労評価モ
デル、発電量評価モデル及びそれらのモデルを開発するためのスーパーシミュレーションが必要と考えられる。また、風力発電の主力電源化を実現するためには、浮体式風力発電の導入も欠かせない。今後、浮体式洋上風車配置、風車・基礎一体解析の技術開発も重要な分野と考えられる。浮体式洋上風車配置に75

おいては、風車動揺を考慮する必要があるため、着床式風車に比べて、風車後流予測、風車制御は、より複雑になり、モデルの確立には、スーパーシミュレーションと AI の連携活用が期待される。
(8-3-3) Task 31 流れモデル及びウェイクモデルのベンチマーク
Task31 は、流れモデル及びウェイクモデルのベンチマークプログラムで、モデルの精度向上、モデル評価手順の合意形成をし、国際的な V&V フレームワークを構築することを目的とした国際研究活動である。研究範囲は、ウィンドファームスケール、シングルウェイクから始まり、マルチウェイク、メソスケール（気象スケール）、風車翼スケール、ロータスケールと、その対象を拡張してきた。
以下に各 Phase における実施内容を示す。本タスクは、Phase3（2021 年 5 月まで）で終了予定である。
Phase1

2011 年～2014 年：マイクロスケールモデル、V&V フレームワークの定義

Phase22015 年～2018 年：スケールを拡張したフレームワークの構築

Phase32018 年～2021 年：国際的 V&V 戦略、国際規格を通した産学連携、オープンサイエンスベンチマーク

ベンチマークのデータは、風洞からフィールドに至り、各国の研究機関から多数のデータが提供された。Fig. 8.3-2 にベンチマークの一覧、Fig. 8.3-3 に検証用データの種別、長所・短所を示す。
ベンチマークプログラムの実施にあたっては、ベンチマークコミュニティの設置・運営、データ取り扱
いに関するルール（秘密情報レベルとアクセシビリティ）の策定、ベンチマークシステムの構築が行われ、これらをガイドラインとしてとりまとめた。また、各ベンチマークのテストケースのベストプラクテ
ィス（優良事例）のガイドラインも作成した。Fig. 8.3-4 及び Fig. 8.3-5 に、ベンチマーク検証用データのアクセシビリティレベルの定義及びベンチマークプロセスを示す。

76

Fig. 8.3-2. Task31 Benchmark list.

Fig. 8.3-3. Types of verification data.

77

Fig. 8.3-4. Accessibility of validation data.

Fig. 8.3-5. Data accessibility and benchmarking process.
今後、「富岳」による洋上風力産業界向けのベンチマークプラットフォームを検討する上で、IEA Wind
Task31 で実施した検討事項及びモデル評価手順のガイドラインは、参考になると考えられる。
(9) 成果発表、広報・アウトリーチ活動
本プロジェクト全体として令和 2 年 10 月 27 日に国内シンポジウム「クリーンエネルギー「富岳」プ
ロジェクトシンポジウム」（第 1 回、ポスト「京」重点課題⑥から数えると第 5 回目）をオンラインにて開催し、本プロジェクトのねらい、途中段階の成果を発表し、社会への広報・普及を行った。洋上風力発
電のみならず、開発中の燃焼器の技術や AI 連携に対しても期待が寄せられた。また令和 3 年 3 月 31 日にオンラインにて洋上風力ワークショップを開催し、コンソーシアム立ち上げの説明と意見交換を行っ78

た。国内外の企業・大学、新聞社からも参加があり、洋上風力発電の技術やコンソーシアムへの関心の高さが伺われた。
そのほか令和 2 年 11 月 27 日に東京大学生産技術研究所革新的シミュレーション研究センター主催、
本プロジェクト及び「航空機フライト試験を代替する近未来型設計技術の先導的実証研究」（課題代表機関：東北大学大学院工学研究科）の共催で、「第 4 回 HPC ものづくり統合ワークショップ」をオンラインにて開催した。
「「富岳」を利用した革新的流体性能予測技術の研究開発」（課題代表機関：東京大学生
産技術研究所）との協働であり、ポスト「京」重点課題⑥と同⑧として活発に行っていた研究の協働を引き継いでいる。

79

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

行事

2020 年 5 月 25 日

第 1 回全体会

2020 年 10 月 27 日

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

2020 年 11 月 27 日

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

2020 年 12 月 18 日

第 2 回全体会

2021 年 3 月 31 日

クリーンエネルギー「富岳」プロジェクト 洋上風力ワークショップ

２－４．実施体制
業務項目（1）次世代火力発電システム及び洋上ウィンドファームのス

担当機関

担当責任者

東京都文京区本郷 7-3-1 国立大学法教授 吉村 忍人東京大学大学院工学系研究科

ーパーシミュレーション（2）個別シミュレーションとAI（機械学習）の連携活用

京都府京都市西京区京都大学桂 C3
棟教授 黒瀬 良一

国立大学法人京都大学（3）総合的事業推進
福岡県福岡市西区元岡 744

国立大教授 渡邊 裕章

学法人九州大学大学院総合理工学研究院
愛知県豊橋市天伯町字雲雀ケ丘１番
１

教授 飯田 明由

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

大学院工学研究科
東京都目黒区駒場 4-6-1 国立大学法教授 加藤 千幸人東京大学生産技術研究所
福岡県春日市春日公園 6-1

国立大

准教授 内田 孝紀

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

80

教授 小野 謙二

別添１

学会等発表実績

1． 学会誌・雑誌等における論文掲載
No. 掲載した論文（発表 発表者氏名
題目）

1

2
3

4

5

6

7

Performancecomparisonsofgeometric
multigrid solvers
andbalancingdomaindecomposition
solvers
流体構造制御連成
解析の開発
Modelingandsimulation of coalgasification on an
entrainedflow
coal gasifier
A cohesive zone
based
DE/FE
coupling approachforinterfacialdebondinganalysisoflaminated glass
FiniteElement
Based Analysis of
Bio-Heat Transfer
in Human Skin,
During Burn andAfterwards
DirectForcing/Fictitious
Domain 法 に よ る
流体剛体連成解析における補間分散関数の改良
AComparison
between Artificial
Intelligence
MethodandStandard
Diagnosis
Methods for Power
Transformer
Dissolved
Gas
AnalysisUsing
Two
PublicDatabases

発表した場所（学会誌・ 発 表 し た
雑誌名等）
時期

国
内・
国際
の別

Hongjie
Zheng,
Ryuji
Shioya,Hiroshi
Kawai,
SebastianReiter,Gabriel Wittum

Computingand
Visualizationin
Science volume 23,
Articlenumber:5(2020), pp.18, DOI:10.1007/s00791020-00323-4, 27

2020/4/1

国際

金子栄樹、吉村忍

第 25 回計算工学講演会
講演論文
Advanced
PowderTechnology

2020/6/1

国内

2020/6/10

国際○

Wei Gao, Xin Liu,
ShunhuaChen,
Tinh Quoc Bui,
Shinobu Yoshimura

TheoreticalandAppliedFracture
Mechanics(108 巻)

2020/7/1

国際○

A.M.M.Mukaddes,
Ryuji Shioya, Masao
Ogino, Dipon Roy,
Rezwan Jaher

International JournalofComputational
Methods, DOI:10.1142/
S0219876220410108,
2041010, pp.1-14, 11

2020/8/1

国際○

川上幸亮、三目直登、 日本計算工学会論文集
金子栄樹、吉村忍

2020/8/7

国内○

Hongjie
Zheng,
Ryuji Shioya

2020/9/1

国際○

Hiroaki Watanabe,
Kurose Ryoichi

IEEJ Transactions on
ElectricalandElectronic
Engineering, Vol. 15,
Issue 9, pp.1305-1311,
DOI:10.1002/tee.23197,

81

査読
（有
りの
場合
○を
記
入）
○

8

9

10

11

12

13

14

15

Numerical
investigationofpremixed flamewall interaction:
Effectivenessof
insulation wall onheat loss reduction
Investigation
of
lean
spray
combustion
instability in a
back
step
combustor using
LES
Directionalmigration of an
impinging droplet
on a surface withwettability
difference
CNN に基づく鶏肝臓画像における脂肪肝分類に対する
転移学習の効果検
証
Robustfluidstructureinteraction
analysisforparametric studyof flapping motion
Unsteady flamelet
modellingfor
N2H4/N2O4 flame
accompaniedbyhypergolic ignition
andthermaldecomposition
Directional
rebounding of adroplet impinging
hydrophobic
surfaceswithroughness
gradientsFlameselfinteractions in an
open turbulent jet
spray flameR.
Kai,
Takahashi,
KuroseA.
R.

Journal of Thermal
ScienceandTechnology, 15, Paper
No.20-00390 (11 pages)(2020)

2020/10/1

国際○

J. Nagao, A. L.
Pillai, R. Kurose

Journal of Thermal
ScienceandTechnology, 15, Paper
No.20-00330 (12 pages)(2020)

2020/11/1

国際

○

Z.
Yuan,
Matsumoto,KuroseM.
R.

PhysicalReview
Fluids, 5, 113605 (19pages) (2020)

2020/11/1

国際○

鄭宏杰, 中林靖, 増田正人, 西宏起, 山
中大介, 高橋伸一郎,
伯野史彦, 宮内陽介,
岡崎敬, 吉田一寛,
塩谷隆二
Giwon
Hong,
Shigeki
Kaneko,Naoto
Mitsume,Tomonori Yamada,Shinobu
Yoshimura
K. Konishi, R. Kai,
R. Kurose

日本計算工学会論文集,Vol.
2020-1,
p.20201003,
ISSN:1347-8826,
DOI:10.11421/jsces.
2020.20201003
Finite Elements inAnalysis and Design

2020/11/1

国内○

2020/12/9

国際○

Applications in Energy
andCombustion
Science, 5, 100022(16pages) (2021)

2021/1/1

国際

○

Z.
Yuan,
Mtaumoto,
KuroseM.
R.

International Journalof Multiphase Flow,
138, 103611 (10 pages)(2021)

2021/2/1

国際○

S.P. Malkeson, U.
Ahmed, A.L. Pillai,
N. Chakraborty, R.
KurosePhysics of Fluids, 33,
035114(11
pages)(2021)

2021/3/1

国際

○

82

2．学会等における口頭・ポスター発表
No. 発表した成果（発表題 発表者氏名
目、口頭・ポスター発表
の別）
1

2

3

4

5

6

7

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

発表した時
期

国内・ 招 待
国際
講演
の別
（○
を記
入）

日本流体力学会年会
2020, オンライン

2020/9/19

国内「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)
FR-Comb による CO2 渡邊裕章
「スーパーシミュレ
回収型 IGCC 用ガス化
ー シ ョ ン と AI を
炉解析、口頭連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)
ADV_Solid の富岳チュ 山田知典
「スーパーシミュレーニング＆非弾性解析、ー シ ョ ン と AI を
口頭連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)洋上風力発電のロード 今井博「スーパーシミュレ
マップと「富岳」への期
ー シ ョ ン と AI を待、口頭連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)
RIAM-COMPACT 富 内田孝紀、小野謙 「スーパーシミュレ
岳チューニングと二ー シ ョ ン と AI を
風車ウエイク解析、口頭連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)FronFlow/blue （FFB ) 加藤千幸
「スーパーシミュレ
Oprimizationforー シ ョ ン と AI を
FugakuandIts
連携活用した実機クResutls、口頭
リーンエネルギーシ

2020/10/1

国内

2020/10/1

国内

2020/10/1

国内

2020/10/1

国内

2020/10/1

国内

2020/10/1

国内

Numerical simulationof atomization andevaporationofkerosenejetincrossflow
underoscillating flowrate conditions, 口頭
燃焼・ガス化炉のスーパ
ーシミュレーション
全体会用、口頭J.
Wen,
Konishi,
Kurose

K.
R.

吉村忍

83

8

洋上ウィンドファーム
スーパーシミュレーシ
ョン
流体構造連成＆疲労損
傷解析、口頭陳順華、吉村忍

9洋上ウィンドファームのマルチフィジクスシ
ミュレーション、口頭

飯田明由

10

洋上風車後流の大規模
乱流解析、口頭

飯田明由

11

富岳における
ADVENTURE の性能
最適化、口頭LES 解析を用いた二次
元翼の遠方後流の予測、
口頭
ス パ コ ン 版 RIAMCOMPACT の開発と風車ウエイクシミュレー
ションへの適用－風洞
模型スケールから実機
スケールまで－、口頭・ポスター
Finite Element Based
Thermal-Solid
Coupling Analysis ofConcrete Slab, oral

河合浩志

15

Designoptimization
with finite element
method and artificial
intelligence, oral

Hongjie Zheng,
Kawai Hiroshi,
Ryuji Shioya16

輻射が壁面衝突噴霧火

村 松 宏 起 , A. L.

12
13

14

水本達也、鈴木康
方、加藤千幸内田孝紀、小野謙二、飯田明由、吉村忍、加藤千幸、
山出吉伸、今村博、
植田祐子
Abul Mukid Md.
Mukaddes, Ryuji
Shioya,Masao
Ogino, Kaworu
Yodo

ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
第 1 回シンポジウム(online)
第 4 回 HPC ものづ
くり統合ワークショップ(online)
第 4 回 HPC ものづ
くり統合ワークショップ(online)
日本機械学会第 98期流体工学部門講演会
第 42 回風力エネル
ギー利用シンポジウ
ム、オンラインInternational
ConferenceonComputational
EngineeringandScience for Safety
and Environmental
Problems(COMPSAFE2020),
Kobe,
Japan(Online)
International
ConferenceonComputational
EngineeringandScience for Safety
and Environmental
Problems(COMPSAFE2020),
Kobe,
Japan(Online)
第 59 回燃焼シンポ

84

2020/10/1

国内

2020/10/1

国内

2020/11/1

国内

2020/11/1

国内

2020/11/11

国内

2020/11/26

国内

2020/12/1

国際

2020/12/1

国際

2020/12/2

国内

炎の熱伝達特性に及ぼす影響、口頭Large Eddy Simulation
を用いたエタノール噴
霧浮き上がり火炎における着火挙動の解析、口頭
Enhancingthestability of Cassie
wetting state using
secondary structures、
口頭
Numerical simulationof atomization andcombustion of liquid
fuel in Sydney piloted
need burner、口頭Structural and fatigue
analysis of a 5 MW
windturbinecomposite
blade
consideringwind
effects, oral

Pillai, 黒瀬良一

ジウム, オンライン

村上太一, Y. Hu,
黒瀬良一

第 58 回燃焼シンポ
ジウム, オンライン

2020/12/2

国内

Z.
Yuan,
Matsumoto,KuroseM.
R.

日本流体力学会 第34 回数値流体力学
シンポジウム, オンライン

2020/12/22

国内

J.Wen, Y. Hu, J.
Iino, T. Nishiie,
R. Kurose

日本流体力学会 第34 回数値流体力学
シンポジウム, オンライン

2020/12/22

国内

Shunhua Chen,
Shinobu
Yoshimura,
KaworuYodo,
Naoto Mitsume,
Yasunori Yusa,
Tomonori
Yamada2021/1/1

国際

21

Artificial Intelligence
Simulation to Predictof Liver Lipid Levels,oral

R. Shioya,
Zheng

2021/1/1

国際

22

Simulationofthe
Human Skin Burns
and Evaluation ofDifferent Treatments
Using
Adventure
System, oral

A. A. Mukaddes,
S. Mridul, R.
Shioya14th
World
CongressonComputational
Mechanics and 8th
European Congress
on Computational
Methods in Applied
SciencesandEngineering
STS28 - Japan
Session2
Computational
mechanics
forintegrityofindustrial
infrastructures,
14th
World
CongressonComputational
Mechanics (WCCM
XIV)
8th
European Congress
on Computational
Methods in Applied
ScienceandEngineering(ECCOMAS 2020),
Paris(Online)
MS368 - OpenSource
Scientific
Software and theApplication to RealWorld
Problems,
14th
World
CongressonComputational
Mechanics (WCCM
XIV)
8th
European Congress
on Computational

2021/1/1

国際

17

18

19

20

H.

85

23

Super-simulations ofClean Energy Systems,
口頭
マスクの流体力学

Yoshimura
Shinobu25

高温反応場を対象とした数値シミュレーショ
ン技術の現状と今後の課題 －石炭燃焼・ガス化を中心にして－、口頭

黒瀬良一

26

クリーンエネルギー「富岳」プロジェクト概要紹
介、口頭

吉村忍

27

洋上風力発電の現状と今後の展望、口頭

今井博

28

Task31：ウィンドファ
ーム流れモデルのベン
チマーク、口頭植田祐子

29

洋上ウィンドファームのスーパーシミュレー
ションの全体像、口頭

飯田明由

30

RIAM-COMPACT の 小野謙二
「富岳」チューニング、
口頭

24

飯田明由

Methods in Applied
ScienceandEngineering(ECCOMAS 2020),
Paris(Online)
The 3rd R-CCS
International
Symposium, online
第 4 回 CAE ワーク
ショップ、オンライン
一般社団法人 石
炭 ・ 炭素資源有効利用研究会 第 1回研究会 ～脱炭素化に向けた今後の炭素資源利用技術～,
エッサム神田ホール、東京「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
洋上風力ワークショップ(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
洋上風力ワークショップ(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
洋上風力ワークショップ(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
洋上風力ワークショップ(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」

86

2021/2/15

国際

2021/3/1

国内

2021/3/9

国内

2021/3/31

国内

2021/3/31

国内

2021/3/31

国内

2021/3/31

国内

2021/3/31

国内

○

○

31

RIAM-COMPACT による風車ウエイク解析、口
頭

内田孝紀

32

洋上大型風車ブレードの流体構造連成・疲労損
傷解析、口頭陳順華、吉村忍33「富岳」の時代のシミュ
レーションとものづく
り、口頭

加藤千幸

34

FFR-CombとADVENTURE 双 方 向連成による石炭ガス化
炉のスーパーシミュレーション、口頭吉村忍

3．受賞等
No. 名称1

第 42 回風力エネルギー
利用シンポジウムポス
ター賞
ス パ コ ン 版 RIAMCOMPACT の開発と風車ウエイクシミュレー
ションへの適用－風洞
模型スケールから実機
スケールまで－洋上風力ワークショップ(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
洋上風力ワークショップ(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
洋上風力ワークショップ(online)「スーパーシミュレー シ ョ ン と AI を連携活用した実機クリーンエネルギーシ
ステムのデジタルツ
インの構築と活用」
洋上風力ワークショップ(online)
燃焼システム用次世
代 CAE コンソーシ
アム 第 5 回総会、
オンライン

2021/3/31

国内

2021/3/31

国内

2021/3/31

国内

2020/9/24

国内授賞機関(学会名
等)

受賞した時
期

内田孝紀、小野謙二、 第 42 回風力エネ
飯田明由、吉村忍、 ル ギ ー 利 用 シ ン
加藤千幸、山出吉伸、 ポジウム（日本風今村博、植田祐子
力エネルギー学会）

2020/11/26

国内・
国際の
別
国内

受賞者氏名

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

名称
日付
クリーンエネルギー「富岳」 2020/7/16
プロジェクト ウェブサイト公開

説明重点課題⑥のサイトをリニューア
ル。プロジェクトのイベントや成果を公開する。

87

備考

備考

5．広報活動等(ワークショップ・研究会等の開催)
No.
1
2
3

名称
クリーンエネルギー「富岳」プロジェ
クト 第 1 回シンポジウム
クリーンエネルギー「富岳」プロジェ
クト 洋上風力ワークショップ
第 4 回 HPC ものづくり統合ワークショップ（東京大学革新的シミュレーション研究センターとの共催

開催日時
2020/10/27

開催場所
オンライン

参加者（人数）692021/3/31

オンライン

249

2020/11/27

オンライン

157

88