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

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

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

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

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

・ADVENTURE の解析結果と AI（機械学習）を連携活用し、石炭ガス化炉の炉容器のデジタルツインを構
築する。
（東大）
・Flamelet 法と AI（機械学習）を連携活用した高精度かつ高効率のサロゲート燃焼モデルを構築し、それを FFR-Comb に実装し、超臨界圧 CO2 ガスタービンのデジタルツインを構築する。
（京大）
・FFB の解析結果と AI（機械学習）を連携活用し、高精度後流分布モデル（サロゲートモデル）を構築する。
（東大、豊橋技科大）・FFB と RC HPC 版の解析結果と AI（機械学習）を連携活用し、RC HPC 版に実装される工学風車モデル
である Line Actuator Model のチューニングを行い、RC HPC 版の高精度化を図り、洋上ウィンドファー
ムのデジタルツインを構築する。（東大・九大・豊橋技科大）（３）総合的事業推進（３－１）産業界等との連携及びコンソーシアムの構築
・燃焼システム用次世代 CAE コンソーシアムと本プロジェクト担当者らとの会合を開き、本プロジェク
トで開発中のスーパーシミュレーション技術の紹介を行うとともに、燃焼系産業界のニーズの聴収を行い、その結果を精査し、本シミュレーション技術の開発に反映する。
（東大・京大・九大）
・2021 年度に開始した次世代 CAE を活用する風力産業コンソーシアム（Internet of Wind energy[IoW]
フォーラム（仮称））の立ち上げ作業を継続し、関係者らと本プロジェクト担当者らとの会合を開き、本プロジェクトで開発中のスーパーシミュレーション技術の紹介を行うとともに、風力産業界のニーズの聴収を行い、その結果を精査し、本シミュレーション技術の開発に反映する。
（東大・九大・豊橋技科大）・IoW フォーラム（仮称）、日本海事協会、日本気象協会、日本風力発電協会との連携を通して、国内外への開発技術の普及体制を構築する。（東大・九大・豊橋技科大）（３－２）開発アプリケーションの維持・普及及び成果広報（東大・京大・九大・豊橋技科大）・ADVENTURE、FFB、FFR-Comb、REVOCAP_Coupler の講習会等を開催し、産業界への技術の普及活動を行う。
・各分野トップレベルの国際学会での発表を行う。
・各分野トップレベルの国際ジャーナルへ論文を投稿する。
・本プロジェクト全体として国内シンポジウム（第３回、ポスト「京」重点課題⑥から数えると第７回目）を開催し、本プロジェクトのねらい、途中段階の成果を発表し、社会への広報・普及を行う。
２－２．実施内容（成果）(1) 次世代火力発電システムのスーパーシミュレーション(1-1) ラボスケール石炭ガス化炉と V&V（1-1-1）連成解析のチューニング
令和 4 年度は、ラボスケール石炭ガス化炉の約 1 億要素の詳細格子における連成解析のチューニング
を実施していたが、「京」や九大 ITO 上では発生していなかったのだが、「富岳」に連成解析システムを移行し計算したところ、計算の途中からバーナーコーナー部に非物理的な高温部が発生したことから計算安定化のためにこの高温発生への対処を行った。

-2-

（1-1-1-1）解析方法
Eulerian-Lagrangian 法による固気二相反応流の Large-eddy simulation (LES)[1,2]を行う。気相の支配方程式は、連続の式、運動量保存式、エネルギー保存式、及び化学種保存式を解いた。SGS 応力モデ
ルには Dynamic Smagorinsky モデル[3]を用いた。固相については、個々の粒子の BBO 方程式をパーセ
ル近似により解き、流体と粒子の干渉効果については、PSI-CELL 法[4]により考慮した。ガス化炉内の化学反応としては、石炭の熱分解による揮発分放出、石炭チャーの固気（ガス化）反応、及び気相反応を考慮した。揮発分放出過程を Eq. (1)により表し、揮発分組成については、Umemoto らの実験データ[5]を用いた。揮発分放出速度については、Eq. (2)に示す単一一次反応モデルにより考慮し、その反応速度パラメータは FLASHCHAIN モデル[6]により見積もった。
(1)

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

(3)𝐶ℎ𝑎𝑟 + 𝐶𝑂 → 2𝐶𝑂(4)𝐶ℎ𝑎𝑟 + 𝐻 𝑂 → 𝐶𝑂 + 𝐻(5)= 𝜂(1 − 𝑥) 1 − Ψln(1 − 𝑥)

⁄

⁄

+
(6)

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

(7)

𝐻 + 𝑂 →𝐻 𝑂

(8)𝐶𝑂 + 𝑂 → 𝐶𝑂

(9)

𝐶𝐻 + 𝐻 𝑂 → 𝐶𝑂 + 3𝐻

(10)

-3-

(11)

𝐶𝑂 + 𝐻 𝑂 ↔ 𝐶𝑂 + 𝐻
𝜔̇ = 𝜔̇ 𝜌̅ , 𝑌 , 𝑇 + 𝐾(12)𝜔̇ 𝜌̅ , 𝑌 , 𝑇 − 𝜔̇ 𝜌̿ , 𝑌 , 𝑇

ここで、𝜔̇ は各反応物の消費速度である。
輻射モデルには Discrete Ordinate 法[10]を用いた。
（1-1-1-2）解析条件
Fig.1.1-1 に、対象とするラボスケール石炭ガス化炉の計算領域ならびに計算格子（バーナ部拡大）を示す。格子数は詳細格子である約 1 億 3 千万要素のものを用いた。ガス化炉の試験条件（解析条件）は、
石炭供給量 100 kg/h、酸素比 0.54 の条件に対して、ガス化剤組成を O2/CO2/N2=25 %/25 %/50 %に設定した。
NS 方程式の対流項の離散化には 2 次中心差分に 5 %の 1 次風上成分をブレンドしたものを用いた。
時間進行にはオイラー陰解法を用い、時間刻みは 1×10-4 s とした。
（1-1-1-3）ヒートスポット発生への対処
Fig.1.1-2 に、連成解析時に途中で発生したヒートスポットを示す。流体側格子から突出したバーナーのエッジ部においてこのような非物理的な高温部が発生し、計算が不安定化することから、FFR-Combの連成解析用通信部にリミッターを付加し、通信を行う温度・熱流束の上限値を設定することで発散を抑止した。
（1-1-1-4）解析結果
Fig.1.1-3 に、
150,000step における連成解析のガス温度及び主要化学種の炉内分布を示す。
150,000step
では初期条件の炉内ガスが出口周辺部に残っており、まだ準定常状態に達していないので引き続き計算が必要であるが、ノード時間の限界からこの時点で本年度の計算は終了となった。
Fig.1.1-4 に、連成計算と単体解析における中心軸上ガス温度分布の計算結果及び実験値を示す。連成
解析は 150,000step の結果を、単体解析は連成解析と同じ 150,000step と最終 step である 320,000stepの双方の結果を記載した。図より、連成解析の 150,000step では出口付近での中心軸上温度が実験値と離れているが、単体解析の 150,000step とは近い値となっている。単体解析の 320,000step では実験値と良好な一致を示していることから、本連成解析でも計算を進めれば同様に良好な結果を得ることが出来ると推察している。

-4-

(a) Computational domain.

(b) Computational grid.

Fig.1.1-1 Computational domain and grid.

Fig.1.1-2 Heat Spot at the edge of burner.

-5-

(a) Temperature(b) H2(c) CO(d) CO2(e) H2O

Fig.1.1-3 Gas temperature and concentrations of major gas components

Fig.1.1-4 Comparison of gas temperature distribution in the axial direction.
参考文献[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-6-

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
gas-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 coalgasifier, 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.
（1-1-2）伝熱・冷却管連成解析コードのチューニングと精度検証
ADVENTURE_Thermal による伝熱・冷却管連成解析コード中の連成境界面での物理量交換処理に関するバグを令和 3 年度に修正し、コンクリートの冷却問題を解き精度検証を行った。令和 4 年度は、この修正済みのコードを用いて、ラボスケール石炭ガス化炉の冷却問題の解析を行った。その作業の中で連成境界面が正しく抽出されていないというバグを発見し、その修正に取り組んだ。また、ラボスケー
ル石炭ガス化炉モデルを用いて解析すると、数値不安定性が現れてしまい、それが解消できずにいたため、モデルの再作成を行った。その結果、妥当な解を得られるようになった。以下に、連成境界面の抽出に関する作業、モデル再作成及び冷却解析について述べる。
（1-1-2-1）連成境界面の抽出に関する作業
本研究のモデルは 4 本の冷却パイプ（パイプ A~D）を有する。それゆえ、連成境界面はパイプ A～D
由来の 4 種に分かれる。従来のコードではその区別を考慮していなかった。具体的には、4 種の連成面上の節点を抽出し、一つのグループ S にまとめ、次にモデルが持つすべての面(全要素数×4[四面体要素であるため])について、面を構成する 3 節点がすべてグループ S に属するかどうかを判定する。もし属する場合には、その面を連成面と判定するというルゴリズムであった。このため、余分な連成面がカウントされていた。修正後のコードでは、パイプ A(B, C, or D)の連成面上の節点をグループ Sa (Sb, Sc or Sd)にまとめるようにした。そして、前述の連成面判定作業において 3 節点がそれぞれのグループ Sa(Sb, Sc
or Sd)に属するかどうかを調べるように変更した。その結果として、Fig.1.1-5 に示すように正しく連成面が抽出できるようになった。

-7-

Fig.1.1-5. Extraction of interfaces between gasifiation reactor and four cooling pipes（1-1-2-2）モデル再作成及び冷却解析の様子「富岳」に領域分割ツール ADVENTURE_Metis をインストールした。領域分割前の入力ファイルは
ローカルの PC によって作成済みであり、それを「富岳」上にアップロードした。この結果、任意のパート分割数、部分領域分割数の領域分割済み入力ファイルを「富岳」上で用意できるようになった。なお、
この領域分割計算は「富岳」のジョブとして実行される。次に示す計算例は、パート分割数 20、部分領
域分割数 400 のケースの結果である。
４本の冷却管モデルには実験に合わせて、流速一定の条件を課し、全ての冷却管の冷却水流入温度として 298.15 K（室温）の温度規定境界条件を与えた。炉容器側の初期の温度条件として一様温度 1273.15
K を課した。炉容器側の時間ステップとして、Δt=0.1[s]を与えた。Discontinuous Galerkin(DG)法で離散化した冷却管ソルバー内のサブサイクルとして、炉容器側の 1 時間ステップ当たりに 1,000 サブサイ
クルの計算を全ての冷却管に対して行うように設定した。簡単のため、炉容器における連成面以外の境界にはすべて断熱条件を課した。Fig.1.1-6 に炉容器の温度分布を示す。4 本の冷却パイプは燃焼部付近に配置されているため、その部分で顕著な温度減少（冷却）がみられる。また、時間が経過するにつれて、冷却が進行している様子も観察される。

-8-

Fig.1.1-6. Temperature distributions in gasification reactor at t = 25, 50, 75 and 100 s
また Fig.1.1-7 に Pipe A~D の冷却水の温度分布を示す。横軸は各パイプの入り口からの距離を示して
いる。

-9-

Fig.1.1-7. Time variations of temperature of cooling water along pipe A, B, C and D
どのパイプでも出口に近づくにつれ冷却水の温度が上昇している。また、冷却が進行するにつれ、交換する熱量が減少していくため、冷却水の温度が低下する。また、冷却開始１分経過後の入り口と出口の温度差はいずれのパイプでもほぼ 100°強であり、定性的に妥当な結果を得ることができた。
(1-2) パイロットスケール超臨界圧 CO2 ガスタービンと V&V(1-2-1) 目的
令和 4 年度は、昨年度までに確立した乱流燃焼モデルの一種である Flamelet 法に AI（機械学習）を連携させた Flamelet/ANN（Artificial Neural Network）法の高精度化を行うともに、本手法を用いたLES（以下、LES/Flamelet/ANN と呼ぶ）をパイロットスケール（商用機の 10 分の 1 程度）の東芝エ
ネルギーシステムズ株式会社（東芝 ESS）製 50MWth 級超臨界 CO2 タービン用燃焼器内の乱流燃焼場[1-4]に適用し、試験結果と比較することによって、同手法の有用性を確認することを目的とした。具体的には、昨年度は、出口圧力 8.5MPa の燃焼試験を対象に、「富岳」を用いた燃焼器全体の計算を実施
し試験で得られた出口温度分布との比較を行うことで、本手法の精度検証を行った[5]。令和 4 年度は、
実際のガスタービンの運転で想定されるタービン入口圧力約 28.5MPa 及び、出口温度 1,150℃の運転条
件[4]を対象として、燃焼器内乱流燃焼場の予測解析を行った。
(1-2-2) 計算方法①

燃焼モデルと支配方程式

本解析では、Flamelet 法として火炎の浮き上がりや消炎を再現可能な Flamelet/progress-variable（FPV）法[6,7]をベースとし、壁面からの熱損失及び実在実気体を考慮可能に拡張したものを用いた[8]。
支配方程式は、連続の式、運動量の保存式（Navier-Stokes 方程式）
、混合分率 Z の保存式、反応進行変数 C の保存式、及びエンタルピーh の保存式で構成される。
② 3 種吹込みガス燃焼の Flamelet 法
従来の Flamelet 法では、燃料と酸化剤の 2 つの独立したガス系統からの流入混合による燃焼のみに対応可能で、燃料、酸化剤に加えて冷却用 CO2 の吹込み等、3 つ以上の独立したガス系統からの吹込みがある条件に適用することができなかった。そこで令和 3 年度に、3 種吹込みガス条件に適用可能な
Flamelet 法の導入を行った[9]。その概念を Fig.1.2-1 に示す。

Fig.1.2-1. Concept of three stream Flamelet mode.

- 10 -

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

(1)

𝑊=

と置けるので、3 種吹込みガス燃焼の FPV 法の場合、Flamelet データベースの物理量φは下記で与えられる。
(2)

𝜙 = 𝜙 𝑍, 𝑍 " , 𝑊 , 𝐶 , 𝛥𝐻

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

(3)

+𝑌

で与える。
③ ニューラルネットワーク（ANN）による Flamelet モデルの構築
Flamelet データベースを引用するパラメータの増加に対応するため、多次元の Flamelet 空間での多様体であるデータベースを基底ベクトル（本解析に於いては𝑍, 𝑍 " , 𝑊 , 𝐶 , ∆𝐻 ）の関数で表現することを考える。こうすることで流体計算時の必要なメモリがデータベースそのものからパラメータベクトルだけとなり、メモリの大幅な節約となる。即ちこれは Flamelet データベースに対する回帰問題となる。
機械学習という枠組みの中でも回帰問題を解くための手法は多く存在するが、ここではディープラー
ニングの手法、即ちニューラルネットワークを用いた関数の表現を採用する。本 5 次元 Flamelet データ
ベースにニューラルネットワークを適用した場合、従来は、入力層のパーセプトロンの数はデータベー
ス空間の基底ベクトル数である 5、出力層のパーセプトロンの数は 1 となり、物理量ごとに ANN のデー
タセットを用意することになる。隠れ層をどの様に設計するかはハイパーパラメータであってそれ自体が最適化の対象となり得る。一般的に層が深いほど、層あたりのパーセプトロンが多いほど再現可能な
関数形状は多くなり、データベースの再現性が高くなることが考えられる。また、入力層のパーセプト
ロンの数が多くなると、各対象条件下において十分な再現精度を保つためには、より多くの層数、各層当たりのパーセプトロンが必要となる。しかし一方で、ネットワークを大きくすれば燃焼計算においてデータベースの基底ベクトルから燃焼状態を推定する際の計算負荷が高くなることから、ネットワーク- 11 -

の構築においてはデータベースの一定の再現性を担保しつつもネットワークができるだけ小さくなるよ
うに留意した。本解析では、本来の Flamelet データベースの 5 次元データベースを、𝑍, 𝑍 " , 𝐶 , ∆𝐻の 4 次
元のデータベースに落とし込み、各𝑊 の ANN データベースを用意し、𝑊 軸に対しては 2 つの ANN デー
タベースを補間することで、再現精度を維持しつつ、ネットワークの複雑さを抑えた。各学習の収束条件として、データベースの値とネットワークによって作成した関数の値の RMSE が、データベースの物
理量の標準偏差の 4.5%未満となるようにした。この条件の下で試行錯誤によって、ネットワーク構成は隠れ層が 7 層、１層当りのパーセプトロンは 9 とした。各パーセプトロンの出力は活性化関数 Swish[10]で変換され、オプティマイザーRadam[11]によって重みが更新される。
④ 計算対象および計算格子
Fig.1.2-2 に、本解析で対象とした超臨界燃焼器の概略を示す。対象としたのは、東芝エネルギーシス
テムズ（東芝 ESS）製 50MWth 級検証試験用超臨界 CO2 タービン及び燃焼器[1-5]である。本ガスター
ビン燃焼器は、酸化剤/燃料(CH4)/冷却用 CO2 により燃焼が行われる。運転条件はタービン入口圧力約
28.5MPa 及び、出口温度 1,150℃ である[4]。Fig.1.2-3 に、解析格子形状の詳細の様子を示す。点数及
びセル数は、それぞれ 5,880 万点及び 1.93 億点である。

Fig.1.2-2. Supercritical CO2 turbine and combustor for demonstration plant(a) Computational domain- 12 -

(b)Grids
Fig.1.2-3. Targeted Combustor.
⑤ 計算条件本解析では、NuFD/FrontFlow/Red（FFR-Comb）[5-8,12,13]を用い、乱流モデルとしては、標準
Smagorinsky モ デ ル を 用 い た 。 解 法 に は 圧 縮 性 流 体 ス キ ー ム を 用 い 、 実 在 気 体 モ デ ル と し てSoave-Redlich-Kwong（SRK）モデル[14]を用いた。
Table1.2-1 に、Flamelet データベースの概要を示す。
各 W のデータベースサイズは、
𝑍 × 𝑍 " × 𝐶 × ∆𝐻 =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と大幅にデータベース容量を削減できることが確認できる。
Table 1.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 から得られた、燃焼器内瞬時温度の中心断面及び垂直断面瞬時分布を示す（28.5 MPa）
。また、Fig.1.2-5 には、出口圧力 8.5MPa 及び 28.5MPa 条件での燃焼器内瞬時温度、流速、及び化学種質量分率（O2, CH4, H2O, O2, CO2, CO）の中心断面瞬時分布を示す。Fig.1.2-6 には、
出口圧力 8.5MPa 及び 28.5MPa 条件でのバーナ近傍の時間平均温度、流速、及び化学種質量分率（CH4,
H2O, CO）の中心断面分布を示す。これらの図から、以下の挙動が確認できる。
・燃焼器上流に設置された燃料ノズルより流入する酸化剤の一部は、外周側のスワラを通ることにより旋回速度が与えられ、燃焼器内の旋回流を形成する。
・燃料ノズルから流入したメタン燃料は、燃料ノズル内周から流入した酸化剤と反応し即座に消費され- 13 -

る。
・燃料ノズル出口中心部分に存在する燃焼領域では、燃料と酸化剤とが燃料過濃条件で燃焼し、大量のCO が生成される。その後、燃料ノズル外周より流入した酸化剤と混合反応することで、CO は速やかに消費される。
・燃料ノズル出口下流で逆流領域が形成され、この領域で保炎される。
・高温の燃焼ガスは、ライナやトランジションピースから流入する冷却用 CO2 により希釈され、燃焼器
出口に流出される。また、この冷却用 CO2 により壁面近傍では温度が低く保たれる。
・運転条件の違いにより、28.5MPa 条件は、8.5MPa 条件と比べて、燃焼温度が高くなる。
・令和 3 年度の出口圧力 8.5MP 条件での出口温度分布の試験結果との比較から、全計測位置において50℃以内と比較的良好な精度で予測できていることが確認できた。また、出口圧力 28.5MPa 条件解析で得られた出口断面平均温度は 1,140℃と、想定される 1,150℃と良好な一致が見られた。
本解析は、「富岳」を用いて 20,000 並列の解析を行い、統計計算を含めた計算を約 3 日で完了した。
以上の結果から、本 LES/Flamelet/ANN が超臨界燃焼場を高精度に予測するうえで有効な手法であると言える。

Fig.1.2-4. Instantaneous distributions of temperature in 28.5 MPa condition.

- 14 -

Fig.1.2-5. Comparison of instantaneous distributions of temperature, velocities, and mass fractionsof O2, CH4, H2O, O2, CO2, CO on combustor central plane.

- 15 -

Fig.1.2-6. Comparison of time-averaged distributions of temperature, streamwise velocities, andmass fractions of CH4, H2O, CO in the vicinity of the burner; the black lines show streamlines.

- 16 -

(1-2-4) まとめ
LES/Flamelet/ANN をパイロットスケールの東芝 ESS 製 50MWth 級超臨界 CO2 タービン用燃焼器内の乱流燃焼場に適用し、試験結果と比較することによって、本手法の精度検証を行った。本研究での成果を下記に示す。
・「富岳」を用いて、パイロットスケール超臨界 CO2 タービン（東芝 ESS 製）の、実圧(28.5MPa)に関するスーパーシミュレーションを実施し、実機の高度化に向けた知見を抽出可能である旨の見通しを得た。
・昨年度の結果も含め燃焼器出口における温度分布を実測値及び想定値と比較することにより、本シミュレーションが十分に高い精度で行われていることを確認した。
・Flamelet 法と AI（機械学習）を連携活用した高精度かつ高効率のサロゲート燃焼モデルを構築し、それを FFR-Comb に実装することによって、超臨界 CO2 タービンのデジタルツイン技術を構築した。
本手法により、シミュレーションの予測精度を保ったままで、Flamelet 法で用いるデータベースの容量、及びシミュレーションの使用メモリを大幅に削減することに成功した。
参考文献[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 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] Kasuya,H., Iwai, Y., Itoh, M., Morisawa,Y., Nishiie, T., Kai, R., Kurose, R., “LES/flamelet/ANNof oxy-fuel combustion for a supercritical CO2 power cycle”, Applications in Energy and Combustion
Science, 12, 100083, 2022.
[6] 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 - Predictions ofNO and soot concentrations –“, Journal of Engineering for Gas Turbines and Power, 135, 091503,
2013.
[7] 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 combustion”, Combustion and Flame,
162, 2621-2637, 2015.
[8] Kishimoto, A., Moriai, H., Takenaka, K., Nishiie, T., Adachi, M., Ogawara, A., Kurose, R.,- 17 -

“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.
[9] 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.
[10] Prajit, R., Barret, Z., Quoc V. L., “Swish: a Self-Gated Activation Function”, arXiv:1710.05941,
2017.
[11] 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.
[12] Kobayashi, Y., Iwai, Y., Itoh, M., Nishiie, T., Kurose, R., “Large-eddy simulation of
fluid-structure interaction for Supercritical CO2 combustor”, The Asian Congress on Gas Turbines
2018 (ACGT2018), Aiina Center, Morioka, Iwate, Japan, 2018.
[13] 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.
[14] Soave, G., “Equilibrium Constants from a Modified Redlich-Kwong Equation of State,”
Chemical Engineering Science., 27, 1197–1203, 1972.
(1-3) ベンチスケール石炭ガス化炉の燃焼流解析
令和 3 年度に実施したベンチスケール石炭ガス化炉の反応モデルチューニング及び 1 億要素簡易格子における V&V の結果を踏まえ、令和 4 年度は 30 億要素のベンチスケール炉石炭ガス化炉詳細格子を使用した燃焼流解析を実施した。
(1-3-1)解析方法
ベンチスケール炉は容積比でラボスケール炉の約 30 倍であることから、ラボスケール炉１億要素詳細
格子と同等の空間解像度での解析を実施するためには約 30 億要素の解析格子が必要となる。解析格子は昨年度に FFR-Comb の refiner 機能により 5,000 万要素格子から作成した 32 億要素格子を使用した。
解析条件は、1-1-1 節のラボスケールガス化炉解析と同様である。ただし、下記のメタン水蒸気改質反応について、順反応に加えて逆反応も考慮する。
（旧）CH4 + H2O ⇒ CO + 3H2(1)（新）CH4 + H2O ⇔ CO + 3H2(2)(1-3-2)解析結果
Fig.1.3-2 に、ガス化炉内に 2 点設置された熱電対によるガス温度の実験値との比較を示す。ただし、熱電対位置、及びガス温度は NDA の関係から代表値により正規化してある。Z=2.43 においては計算値
が実験値よりもやや高くなっているが、本年度の解析においてはノード時間の上限到達により炉内が定常に到達したとみなせる状態まで計算を進めることができなかったが、計算を続行すればより実験値に- 18 -

近づく可能性がある。また、炉内の初期状態のガスが完全に抜けきっていないことから出口ガス組成の比較は行えなかった。

Fig. 1.3-2 Comparison of gas temperature at two axial points.

(1-4) ベンチスケール石炭ガス化炉の燃焼流・伝熱連成解析
ベンチスケール石炭ガス化炉の燃焼流・伝熱連成解析の解析条件設定を行った。
(1-4-1) 炉容器側の解析条件設定
炉容器側の解析条件は、(A)物性値、(B)燃焼流解析との連成界面における連成条件、(C)冷却水熱流体
解析との連成界面における連成条件、及び(D)非連成面（外気との接触面）における通常の境界条件から構成される。(A)物性値に関しては、メッシュ作成時の一部形状簡略化の影響を考慮しながら実機の設計図面を基に各部材に対する物性値の割当を決定した。(B)については、これまで実施してきたラボスケー
ル石炭ガス化炉の連成条件を踏襲し、燃焼流ソルバーFFR-Comb で得られた連成界面上の熱流束を炉容器伝熱ソルバーADVENTURE_Thermal に渡し、逆に ADVENTURE_Thermal で得られた連成界面上の温度を FFR-Comb に渡すこととした。連成界面は燃焼流解析と炉容器伝熱解析で用いる格子形状の違いを考慮しながら、燃焼室内壁部に設定することとした。なお、燃焼室内壁部には本来微粉炭噴射口などの突出部が存在しているが、今回作成したメッシュでは簡略化のためこれらの突出部を省略したため、突出部の連成界面については考慮しなかった。(C)に関しては、スラグホール冷却管以外の冷却管は部材
間の 2 面で熱移動を行うように簡略化したため、実機実験における流量及び入口温度と出口温度の差か
ら炉容器から冷却水への移動熱量を計算し、それを 2 面の表面積で割ることで得られた定常熱流束を通常の境界条件として与えた。一方、スラグホール冷却管については冷却水熱流体との連成解析を行うた
め 、 ラ ボス ケー ル 石 炭 ガ ス 化炉 で 行 っ た 炉 容 器 伝 熱 ・ 冷 却 水 熱 流 体 連 成 解 析 の 手 法 を 導 入 し 、ADVENTURE_Thermal で得られた連成界面（冷却管表面）の温度を１次元熱流体ソルバーに渡し、１- 19 -

次元熱流体ソルバーで得られた熱流束を ADVENTURE_Thermal に渡すこととした。また、参考情報として、他の冷却管同様に実機実験における計測値から、定常状態における炉容器から冷却水への移動熱量を求めた。(D)については、炉容器が外気に露出する表面形状を特定し、実機実験の計測値から得られた参照温度と熱伝達係数を与える熱伝達境界条件を設定した。
(1-4-2) 燃焼流側の解析条件設定
燃焼流側の解析条件は、(E)燃焼流の流体解析条件、(F)炉容器解析との連成界面における連成条件から構成される。(E)については、(1-3-1)節における燃焼流単体解析の解析条件を踏襲した。(F)については、(1-4-1)節にて既述の通り、FFR-Comb から ADVENTURE_Thermal に熱流束を渡すこととしており、熱流束計算の実機試験の値から熱伝達率は 1,000W/m2・K と見積もった。連成界面のうち、熱流体モデル
から炉容器側に突出する灯油燃焼部の管路に関しては、側面と端部を非連成面とし境界面の熱伝達条件
は断熱条件とした。炉容器構造体モデルから熱流体モデル側に突出するコンバスタバーナー、チャーバーナー、リダクタバーナーの管路に関しては、側面と端部を連成界面とした。
(1-4-3) 連成境界面情報ファイルの出力
連成解析の実施のため、燃焼流解析側の前処理において連成境界面情報ファイルの出力を実施した。
本年度作業においては、連成境界面情報ファイルの出力作業を実施中に「富岳」のノード時間が上限に到達してしまったことから本連成解析への移行は行えなかった。
(1-5) ベンチスケール石炭ガス化炉の非弾性熱応力解析
本項目ではベンチスケール石炭ガス化炉について 5 サイクルの繰り返し温度サイクルを設定した熱弾塑性解析を行った。ベンチスケール炉は鋼材と、２種類の耐火材 TOCAST-13L、HE0152C の 3 つの材質から構成されており、このうち非弾性モデルとして鋼材部に Ohno-Wang モデルを適用した。なお、
ベンチスケール炉の実機では各材質の間には緩衝材が配置されており、力学的に非連続となっているた
め熱応力が発生しないようになっている。しかし、ベンチスケール炉の炉容器モデルでは、上述した３種類の各部材が結合されたモデル化を行っていたため、物性値の違いによる熱応力が発生しやすい設定となっている。炉容器部に関して実機と解析モデルに相違点があることを前提として、ベンチスケール
炉レベルの大規模複雑形状モデルにおいて、並列構造解析コード ADVENTURE_Solid3 の非弾性解析機
能が正常に機能するかどうかの検証を行った。なお、ベンチスケール炉よりもさらに巨大となる商用炉
では、炉容器部分には耐火材は採用されず、冷却管が組み込まれた鋼構造のみによって構成されることを、追記しておく。
(1-5-1)

ベンチスケール炉の熱弾塑性モデル設定

各部材の弾性設定値は Table 1.5-1 のように設定した。ここで、鋼材は SUS304 を想定した値としている。耐火材である TOCAST-13L 及び HE0152C はキャスタブル耐火物となっており、一般的な耐火レン
ガの物性値を文献[1]の物性データを参考に設定した。
Table 1.5-1. Elastic material properties of bench scale reactor- 20 -

ヤング率[MPa]鋼材

TOCAST-13L

HE0152C

1.93E+05

1.96E+04

1.96E+04

0.29

0.18

0.18

1.87E-05

7.00E-06

7.00E-06

ポアソン比
線膨張係数[1/K]

鋼材の非弾性物性値は文献[2]を参考にして、初期降伏応力を Table 1.5-2、Ohno-Wang モデルの移動
硬化パラメータを表 Table 1.5-3 のように設定した。ここで、本解析の温度範囲ではクリープの寄与は小さいものと想定されるため Ohno-Wang モデルのみを設定した。
Table 1.5-2. Initial yield stress and work hardening parameter of SUS304
初期降伏応力[MPa]2.550E+02

加工硬化係数[MPa]3.860E+03

Table 1.5-3. Ohno-Wang model’s parameters of SUS304
モード

ζ

r

[MPa]

1

8.000E+03

7.351E+00

24.000E+03

2.953E+00

3

3.333E+03

1.531E+00

42.000E+03

1.492E+00

実物のベンチスケール炉は拘束による熱応力が発生しないように、中央部付近のフランジ部分で吊り下げた状態で設置されており、水平周方向には自由に変形できるようになっている。そこで、本非弾性
熱応力解析についても、吊り下げた状態を再現するために周方向フリーの境界条件を Fig.1.5-1 のように課した。

- 21 -

Fig. 1.5-1. Boundary condition of bench scale reactor(1-5-2)

ベンチスケール炉の温度サイクル設定

ベンチスケール炉の連成解析が実施できれば、そこで得られた温度分布やその履歴を入力として非弾性熱応力解析を実施する予定であった。しかし、(1-4-3)で述べたようにペンチスケール炉の連成解析を実施できなかった。そこで、代わりに、ベンチスケール炉の炉容器に、低温側=300K、高温側=1000K の温度変化を５サイクル繰り返す条件で伝熱解析を行い、そこで得られた温度分布の時間変化を与えて非弾性熱応力解析をすることとした。Fig.1.5-2 に伝熱解析から得られた温度分布を示す。この最低温度(A)と最高温度(B)の時刻の間は 16 分割して線形補間して温度分布を計算し、最低温度から最高温度を経て最低温度に戻る温度サイクルを 5 サイクル繰り返すことにより温度データを作成した。Fig.1.5-3 に温度サイクルの設定を示す。

初期温度・最低温度時刻(A)最高温度時刻(B)

節点温度分布
[K]

300K 一様
Fig. 1.5-2. Temperature distribution of bench scale reactor- 22 -

温度サイクル設定
1200

B

1000
温度[K]

800
600
400

A

200
0
0

3600

7200

10800 14400 18000 21600 25200 28800 32400 36000
時刻[s]Fig. 1.5-3. Thermal cycle of bench scale reactor(1-5-3)

ベンチスケール炉の熱弾塑性解析結果

Fig.1.5-3 の温度サイクルで熱弾塑性解析を行い、t=3,600s から t=32,400s までの最高温度点における節点相当応力及び節点相当塑性ひずみの分布を Fig.1.5-4 から Fig.1.5-8 に示す。これらの図より、最高温度点における熱応力が高くなる部分が次第に広がっていることがわかる。さらに、相当塑性ひずみを見るとスラグホール部の相当塑性ひずみが温度サイクルごとに増加する結果となった。ただし、この結果については、先に述べたように、拘束条件、温度場、物性値、各部材の接合状況などについて実機を十分に反映していないことを十分に考慮する必要がある。
時刻[s]

3,600

時刻[s]

節点相当応力節点相当塑性ひ分布[MPa]ずみ分布[-]

3,600

Fig. 1.5-4. Distributions of nodal equivalent stress and equivalent plastic strain at t=3,600s- 23 -

時刻[s]

10,800

時刻[s]節点相当応力節点相当塑性ひ分布[MPa]ずみ

10,800分布[-]

Fig. 1.5-5. Distributions of nodal equivalent stress and equivalent plastic strain at t=10,800s
時刻[s]18,000

時刻[s]節点相当応力節点相当塑性ひ分布[MPa]ずみ分布[-]

18,000Fig. 1.5-6. Distributions of nodal equivalent stress and equivalent plastic strain at t=18,000s
時刻[s]25,200

時刻[s]- 24 -

25,200節点相当応力節点相当塑性ひ分布[MPa]ずみ分布[-]

Fig. 1.5-7. Distributions of nodal equivalent stress and equivalent plastic strain at t=25200s
時刻[s]32,400

時刻[s]節点相当応力節点相当塑性ひ分布[MPa]ずみ分布[-]

32,400Fig. 1.5-8. Distributions of nodal equivalent stress and equivalent plastic strain at t=32,400s(1-5-4)

スラグタップ冷却管の応力-ひずみ線図評価

冷却管に発生する非弾性ひずみの評価を行うため、スラグタップ直下の冷却管について応力－ひずみ線図をプロットした。冷却管付近のサンプリング点を Fig.1.5-9 に示す。サンプリング点について応力－
ひずみ線図をサンプリングした結果を Fig.1.5-10 から Fig.1.5-12 に示す。これを見るとひずみが比較的大きくなっている成分は xx、yy、zz、yz となった。このうち、xx、zz、yz では応力－ひずみ線図のヒス- 25 -

テリシスが大きくループ状になっていることがわかる。これより、冷却管で塑性変形が増加している過程がとらえられていると考えられる。この結果に関しても、実機では熱応力が発生しないように各部材の間に緩衝材が配置されているのに対して、本解析ではその効果が反映できていないため、本解析では拘束条件、温度場、物性値、各部材の接合状況などは実機を十分に反映できていないことに留意する必要がある。

Fig. 1.5-9. Sampling point for stress-strain curve

節点ID=9545636

応力-ひずみ線図(xx)

節点ID=9545636

600
500

300

400

200

200

応力yy[MPa]

応力xx[MPa]

300
100
0
-100

100
0
-100

-200
-300

-200

-400
-500
-0.0005

応力-ひずみ線図(yy)40000.00050.001

0.0015-300-0.003

0.002

ひずみxx

-0.002-0.001

0

0.001
ひずみyy

Fig. 1.5-10. Stress-strain curve (xx and yy components)- 26 -

0.002

0.003

0.004

節点ID=9545636

応力-ひずみ線図(zz)

節点ID=9545636

120

3

100

2.5

80

2

応力xy[MPa]60

応力zz[MPa]

応力-ひずみ線図(xy)

40
20
0

1.5
1
0.5
0

-20

-0.5

-40
-60
-0.016-0.014-0.012-0.01-0.008-0.006-0.004-0.002-1-5.00E-05

00.00E+00

5.00E-05

ひずみzz

1.00E-04

1.50E-04

2.00E-04

ひずみxy

Fig. 1.5-11. Stress-strain curve (zz and xy components)

節点ID=9545636

節点ID=9545636

応力-ひずみ線図(yz)

応力-ひずみ線図(zx)

0.8

150

0.6

100

応力zx[MPa]

応力yz[MPa]

0.4
50

0

0.2
0
-0.2

-50

-0.4

-100
0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

ひずみyz

-0.6
-1.00E-05 0.00E+00 1.00E-05 2.00E-05 3.00E-05 4.00E-05
ひずみzx

Fig. 1.5-12. Stress-strain curve (yz and zx components)- 27 -

5.00E-05 6.00E-05

参考文献[1] 宮本 学・尾上 俊雄・成田 貴一 「高炉用耐火物の高温における変形および破壊挙動」鉄と鋼 第
67 年第 2 号,pp.313-322 ,1981[2] 谷本 浩一・高橋 由紀夫「熱ラチェット変形解析精度向上のための非弾性構成式の研究」電中研報
告：T93022, 1993(2) 洋上ウィンドファームのスーパーシミュレーション(2-1) NREL5MW 風車の後流解析と V&V、後流分布モデル構築
令和 4 年度は、令和 3 年度に実施した NREL5MW（風車直径 D＝126 m, 出力 5MW）の解析（風車
下流距離 12D）をさらに長時間実施（実時間 10 分）し、後流モデルを作成するためのデータベースを構築した。Fig.2.1-1 は LES 解析による風車のハブ高さ位置における速度分布（主流方向）である。このデータを元に後述（4-3 節「後流分布モデル化に向けた FFB と AI 連携活用」）する Wake モデルを開発し、実際の洋上風車における風車下流速度の実験結果との比較を行った。Fig.2.1-2 に実験結果との比較を示す。〇印が Lidar を用いた実測結果（Areva M5000 Multibrid、D＝116 m、5MW）、黒の実線が本解析をもとに開発した Wake モデルによる計算結果、
赤及び青線が従来の後流モデルによる予測結果である。
本研究で開発した後流モデルは実際の洋上ウィンドファームの速度場を広い範囲で推定することが可能であることがわかる。従来の Wake モデルは洋上ウィンドファームの標準的な風車設置間隔に相当する
遠方領域（X/D >6）の速度分布を予測可能であったが、風車近傍(X/D < 5 以下)の予測ができなかった。
本研究では LES による詳細解析により風車近傍の流れ場データをもとに Wake モデルを構築したことにより風車近傍の予測を行うことが可能となった。これらの結果より、精度の高い解析が実施できたことを確認した。

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

Fig.2.1-2. Wake centerline velocity behind AREVA M5000 turbine within a full-wake region(2-2) NREL5MW のタンデム配置風車の後流解析
NREL5MW を 2 台及び 3 台並べたタンデム配置風車の解析を実施し、風車の相互干渉の影響を調べた
結果を Table2.2-1 に示す。なお、この解析では、２ないし３台の風車の回転速度はすべて同一とした。
ただし、現実には、後方配置の風車は Wake の影響により、回転速度が前方配置の風車よりも遅くなる
ことが予想される。今回の解析では、そうした効果が考慮されていないことに注意を要する。
今回の解析では、風車間距離が狭い（3D）場合は、下流側の風車は十分な出力が得られなくなることが確認された。風車 3 台を並べた解析を令和 3 年度に引き続き実施した。前年度の解析は解析時間が実時間に換算して 2 分程度であったが、今回の解析では実時間で 10 分程度行い、3 台目の風車下流まで流れを発達させることができた。解析結果より 2 台目の風車の性能低下と比較して、3 台目の風車の性能低下のほうが小さいことが確認された。実際の風車でも 3 台目以降の風車性能の低下が 2 台目と比較して小さいことが知られており、解析結果はこのことを裏付けているものと考えられる。
Fig.2.2-1 及び Fig.2.2-2 に風車下流の平均速度分布、乱れ強度分布を示す。2 台目の風車に流入する主流速度は 1 台目の後流の影響により減速していることがわかる。2 台目の風車出力が低いため、2 台目の風車下流の速度欠損が小さくなり、3 台目の風車の流入速度は 2 台目よりも大きくなっている。風車に流入する乱れ強度については 2 台目、3 台目はほぼ同じ程度であることから流入乱れの影響は平均速度の影
響よりも小さいと考えられる。1 台目の風車の流入乱れ（大気乱流）は 8%程度であるのに対して、下流側の風車への流入乱れは 10%以上であり、大気乱流よりも大きな乱れが流入している。特に 3 台目の風車に対してはブレード先端（回転半径）付近で非常に強い乱れが流入しているにも関わらず、風車性能
が 2 台目よりも低下していないことか確認された。
風車の流入乱れは流入速度に対するものであるため、10%以上と大きく見積もられるが、ブレードに対する相対速度は風車の周速度が速いことから流入乱れの影響は相対的に小さくなる。このことからも平均速度の減少のほうが風車性能に及ぼす影響が大きいと考えられる。
Table2.2-1. Analysis conditions for tandem wind turbines- 29 -

周速比

一様流

4

0.36/0.1

一様流

5

0.41/0.0

一様流

6

一様流

7

0.48/0.0

大気乱流

70.48/0.0

大気乱流

7

大気乱流

5

風車配置

風車 2 台（風車間隔 3D)

解析規模

要素数 6.5 億

風車 3 台要素数 18 億（風車間隔 7.5D)（ベース 2.2 憶）

回転回数

Cp

流入条件

40 回転

0.45/0.0

0.46/0.39/0.42
100 回転Fig.2.2-1. Geometric averaging of the wind velocity potential behind the triple turbines- 30 -

Fig.2.2-2. Geometric averaging of the Tix behind the triple turbines.
複数台の風車を配置した解析により風車の Wake が風車性能に及ぼす影響を確認することが可能とな
り、Wake による流入乱れの影響よりも平均速度分布の影響が大きいことがわかった。風車間距離が 7.5倍程度の場合、下流側風車の流入速度は 1 台目の 70％程度であった。単体風車の解析より 2 台目の風車の流入速度が主流速度に回復するのは X/D>10 であり、風車間距離は Wake のことだけを考えるのであれば風車直径の 10 倍程度必要となる。実際の設計ではケーブル長やメンテナンス性、経済性を考慮する必要があり、単純に X/D>10 とすることはできないが、流体力学的には、風車 Wake の速度欠損の影響を考慮すると X/D>10 が望ましいことがわかった。
(2-3) RC HPC 版と FFB による NREL5MW 風車の解析結果の比較
令和 3 年度までの RC HPC 版と FFB との定量的な比較（風車ハブ高さにおける平均風速分布の比較）と議論に基づいて、RC HPC 版に実装されている風車のブレード回転を模擬するアクチュエータライン
モデルにおける抵抗係数の与え方に改良を施した。計算の諸条件に関して、本計算における格子分解能
は各方向においてΔx=0.5m, Δy=Δz=2.0m の等間隔とした。各方向の格子数は Nx=1261、Ny=251、Nz=231 であり、総格子点数は約 7,300 万点である。無次元時間刻みは 5.0×10-4 であり、無次元時間 50～100 において時間平均を施した。瞬間場の観察からは、風車下流に形成されているウエイク領域がブレードからの翼先端渦、ナセル、タワーが形成する剥離流の影響を強く受けていることが確認された。特に令和４年度は、令和３年度に見い出した風車ブレード根本付近の抵抗係数の数値は固定（一定）とした上で、その設定範囲の影響を検討した。その結果、Fig.2.3-1 に示すように、令和３年度に設定した影
響範囲の結果が FFB とは良好な一致を示すことが明らかになった。

- 31 -

Fig.2.3-1. Comparison of average wind speed at wind turbine hub height(2-4) RC HPC 版の Actuator Line Model のパラメータチューニングと Alpha Ventus の解析と V&V
ここでは、RC HPC 版を「Alpha Ventus Offshore Wind Farm (Fig.2.4-1)」に適用した結果を議論する。本研究では、7 号機風車の Wake 領域内の風速分布（風車ハブ高さ）に関して、実測データ[1]と比較した。
Fig.2.4-2 に示すように、7 号機の風車 Wake 計測は赤色の点線が示す領域内で Galion Lidar を用いて実施された結果である。
本計算では風向を固定した場合と水平断面内において風向を周期的に変化させた場合を検討した（α =0、2.5、5.0、6.25、7.5、10.0 deg の計 6 ケース）。
Fig.2.4-3 には、風向変化 α = 6.25 deg における計算結果（風車ハブ高さにおけるスカラー水平風速の可視化）を示す。Fig.2.4-3(a)に示す瞬間場に注目すると、風車下流に形成された Wake は大きく蛇行し、揺らいでいる様子が見てとれる。こうした非常な Wake 挙動を反映して、Fig.2.4-3(b)に示す時間平均場の結果では、風車下流に向かって Wake 領域が水平方向に徐々に拡大する結果が得られている。
Fig.2.4-4 には、風車ハブ高さ（91.5m）におけるスカラー平均風速の比較を示す。この図を吟味すると、周期的に変動させる風向角度を徐々に増加させることで、実測値の傾向に近づくことが示された。
以上、実スケールの「Alpha Ventus Offshore Wind Farm」の 7 号機風車に着目した Wake シミュレーションでは、流入風の周期的変動を考慮する手法が実測値の挙動を再現するのに非常に有効であることが明らかになった。

- 32 -

Fig.2.4-1. Alpha Ventus Offshore Wind Farm

Fig.2.4-2. Full scale Alpha Ventus Offshore Wind Farm layout- 33 -

(a) Momentary field(b) time average field
Fig.2.4-3. Calculation results for the full-scale Alpha Ventus Offshore Wind Farm No. 7 wind turbine,
wind direction change α = 6.25 deg

Fig.2.4-4. Comparison of scalar average wind speed at hub height (91.5m) of No. 7 turbine atfull-scale Alpha Ventus Offshore Wind Farm
参考文献
[1] Gallacher, Daniel and More, G., Lidar measurements and visualisation of turbulence and wakedecay length, European Wind Energy Association Conference and Exhibition 2014, EWEA 2014.
(2-5) NREL5MW 風車の累積疲労損傷解析
令和 3 年度に、それまでワークステーションで実施していた累積疲労損傷解析コード
ADVENTURE_Fatigue を「富岳」に移植した。これにより、流体-構造-疲労損傷の全ての解析プロセス
を「富岳」上で連続して実行する環境を実現した。動作試験として、NREL5MW の風車 2 基を離間距離- 34 -

3D（=378m）で配置したタンデム配置風車の FFB による LES 解析から得られたブレードに作用する流体力分布を与えたブレードの動的応答解析、その応力履歴から累積疲労損傷分布を求める解析を行った。流体力履歴の抽出とデータ形式の変換には REVOCAP_Coupler とファイルフォーマット変換のためのユーティリティプログラムを使用した。解析ケースは、一様分布の流入風を仮定し、TSR=5、6、7、8 に対して、入力風に乱れ無しと乱れ有りの解析の合計８ケースである。また、今回実施した２基のタンデム解析では、本来各風況ケースに対して、さらに前方配置風車のみならず、Wake の影響を受ける後方配置風車に作用する流体力を抽出した解析も可能である。しかし、(2-2)節で述べたように、今回行った離間距離 3D のタンデム解析では、後方配置風車の回転速度が前方配置風車の回転速度と同一（同位相）に設定していたため、後方配置風車へ流入風速度の低下を考慮できていない。このため、今回得られた
後方配置風車に作用する流体力は実体を反映しているとは言えない。そこで、今回は、前方配置風車に生じる累積疲労損傷解析のみを実施した。
次に具体的な解析手順と解析結果について順に述べる。まず、各解析ケースにおいて、FFB による LES
解析では、数 10 回転分の解析を行い、流体場がほぼ定常になった時の風車２回転分の解析結果から、ブレードに作用する流体力を抽出し、その値を REVOCAP_Coupler で変換し、ブレードの構造解析モデルに流体力として負荷した。なお、ADVENTURE_Solid によるブレードの動的構造解析を行うにあたり、
荷重ゼロの状態から数十ステップかけて、所定の入力荷重レベルまで徐々に増加させた。この結果、流体力を負荷することによる人工的な衝撃力の発生を回避した。
風車２回転分の動的応答解析結果から得られた各積分点での応力の時刻歴データに対して、Miner 則と Goodman 線図を用いて疲労損傷の評価を行った。
Fig.2.5-1～2.5-4 に入力風に乱れを導入した４ケース(TSR=5, 6, 7, 8)の年間累積疲労損傷値の分布を示す。ここで、損傷度が 10-5 以上の要素のみを色コンターで表示し、それ以外の要素は灰色で半透明表示している。色コンターは青が 10-5、赤が 10-2 以上である。また、各図において、上の図は全体形状がわかるように少し斜め上から見た視点の表示となっており、下の図は、ブレード両面への高疲労箇所の分布がわかるように xz 平面鉛直方向視点の図となっている。これらの図から、NREL5MW ブレードにおいては、翼根から計った翼長方向の座標で、y=7.539m と y=40.12m 付近で相対的に大きな疲労損傷が生じていることがわかる。この部分は、どちらもブレードの構造不連続部に相当している。また、TSR=8
から 7 になる段階で、累積疲労損傷値が少し減少するが、TSR=7 から 6、5 となるにしたがって、累積
疲労損傷値が増加することが示された。

- 35 -

Fig. 2.5-1. Distribution of annual accumulated fatigue damage for turbulent inflow of TSR=5.
(Maximum fatigue damage = 0.858)

Fig. 2.5-2. Distribution of annual accumulated fatigue damage for turbulent inflow of TSR=6.
(Maximum fatigue damage = 0.158)

Fig. 2.5-3. Distribution of annual accumulated fatigue damage for turbulent inflow of TSR=7.
(Maximum fatigue damage = 0.0019)- 36 -

Fig. 2.5-4. Distribution of annual accumulated fatigue damage for turbulent inflow of TSR=8.
(Maximum fatigue damage = 0.0088)(2-6) NREL15MW 風車の解析モデル構築
これまで NREL5MW 風車で計算してきたが、今後国内でも建設が予想される 15MW 風車について、現時点では実物がなく、シミュレーションによって評価するには新たなメッシュが必要である。そこで、
令和３年度には、米国再生可能エネルギー研究所（NREL）が公開されている情報をもとに、風車単機を含む風況解析のための並列 LES 解析コード FFB 用メッシュを作成した。具体的には、風況解析メッシ
ュの中に、回転する風車ブレードを考慮したオーバーセットメッシュを組み込むこととし、メッシュを可視化した上で品質の評価を行った。
令和 4 年度は、NREL15MW 風車のブレードメッシュ（６面体要素）を構築した。構築したメッシュ
モデルの規模は、ボリューム数 128、総要素数 137,310、総節点数 273,272、総自由度数 819,816、要素は６面体要素である。Fig.2.6-1 は CAD モデルを、Fig.2.6-2 にメッシュを示す。

Fig. 2.6-1. CAD model of a blade of NREL15MW wind turbine.

- 37 -

Fig. 2.6-2. Mesh of a blade of NREL15MW wind turbine.
(3) アプリケーションの「富岳」における性能評価とチューニング
(3-1) FFR-Comb(3-1-1) 可視化ファイル生成モジュール ffr2viz の高速化
FFR-Comb はメインソルバーでの計算後、可視化ファイル生成モジュール ffr2viz を使用して計算結果を可視化ソフト用データに変換する。この ffr2viz について、数十億要素規模の結果処理を効率化するた
め、昨年度実装した任意断面のみの可視化ファイル出力機能に加えて任意境界面の出力機能及び並列実
行機能を実装した。改修後の ffr2viz では、メインソルバー実行時に領域分割を行った各領域の計算結果の中から、設定した断面または境界面の存在する領域のみを選択して結果ファイルの生成を並列実行し、可視化ソフト上において結果ファイルを統合する。上記の改修によって、従来は 24 時間以上必要であった可視化ファイルの生成が 30 分前後の処理時間へと短縮された。これにより、「富岳」を使用した数十億要素規模解析の結果処理が大幅に効率化された。
(3-1-2) メインソルバーfflowHPC のメモリ軽量化「富岳」での Flamelet 計算におけるメモリ量削減を目標とし、Flamelet/ANN モデルの実装と消費メ
モリ量の計測を行った[1]。
従来の Flamelet テーブルを用いた解析では、テーブルの値を線形補間して出力するため、計算時にはテーブル内の全ての点のデータを保持する必要がある。そのため、テーブルのサイズに応じたメモリ容量が必要となり、「富岳」の単位コアあたりのメモリが 1.3GB であることから、高次元のテーブルではしばしば計算が不可能となる。そこで、テーブルを ANN に学習させる方法を用いて計算時のメモリ軽量化を図った。今回、ANN には典型的な全結合型ニューラルネットワークを用いる。ANN は、元のテーブ- 38 -

ルの入力パラメータである(Z, Z^', C, Δh)を入力層とし、温度、密度、化学種濃度などのテーブルの物理量を出力パラメータとする。計算の前にテーブルを ANN に学習させることにより、テーブルが非線形関数として表現されるため、計算時にはテーブルのデータを保持することなく計算できメモリ容量の大幅な削減が期待できる。
従 来 型 Flamelet テ ー ブ ル と ANN/Flamelet テ ー ブ ル の デ ー タ 容 量 を Table3.2-1 に 示 す 。
ANN/Flamelet テーブルは従来型 Flamelet テーブルの約 1/2,000 のデータ容量であり、Flamelet 計算時の使用メモリを大幅に削減することができる。10,000step 計算した際の計算時間は従来型 Flamelet テー
ブルが 33 時間、ANN/Flamelet テーブルが 31 時間で大きな差は無かったが、ANN/Flamelet テーブルは使用メモリ量が少ないことから、1 ノードあたりの使用 CPU 数を増やして同一計算資源での計算速度を向上させることも可能である。
Table 3.2-1. Comparison of memory requirements between conventional and ANN libraries
Memory(MB)Conventional library

ANN library

5220

2.160

(3-1-3) 前処理モジュール prefflow の高速化
前処理モジュール prefflow に関して、回転機械等に使用する periodic 境界面や sliding 境界面といった２つの面のペアリングが必要な境界面のペアリング面探索アルゴリズム及びメモリ消費を改修した。
上記の改修により、ペアリング面を持つ格子に対しての prefflow 実行速度が約 2 倍となった。
参考文献[1]H. Kasuya, Y. Iwai, M. Itoh, Y. Morisawa, T. Nishiie, R. Kai,R. Kurose, LES/flamelet/ANN ofoxy-fuel combustion for a supercritical CO2 power cycle, Applications in Energy and Combustion
Science, 12,10083. 2022(3-2) FFB
風車まわり流れ解析を高速化するために改良された FFB の領域分割ツールを風車まわり流れ解析に適用した。
領域分割ツールは、「『富岳』を利用した革新的流体性能予測技術の研究開発（JPMXP1020200311）」
で開発されたものであり、本報告ではその概要を示す。FFB で領域分割をする場合、要素の接続情報をグラフ化しこれを分割する Recursive Graph Bisection (RGB)を用いる。グラフデータにはオーバーセッ
トデータの接続情報が反映されないため、RGB で領域分割すると、オーバーセット通信が特定の領域に局在してしまいこれが並列性能の低下の要因となる課題があった。これを解決するため、グラフを作成する際に、オーバーセットデータを参照し、グラフデータにオーバーセットの接続情報を反映できるように改良した。これを実現するため、グラフデータにオーバーセットにより接続される要素をつなぐブリッジ要素を設け、これをグラフに追加した（Fig.3.2-1）。
改良した領域分割ツールの効果を検証するため、NREL5MW 風車を直列に 3 台つなげた計算モデルに- 39 -

対しテスト計算を実施した。計算格子の総要素数は約 2.2 億要素であり、「富岳」
400 ノード（1,600 CMG）
を用いて計算し、50 ステップの計算時間の比較をした。領域分割は従来の RGB によるケース（以下、ケース A）とグラフ作成時にブリッジ要素を用いる改良方法によるケース（以下、ケース B）の 2 通りを実施した。Fig.3.2-2 に隣接通信及びオーバーセット通信の通信パターンを示す。隣接通信に関しては、
ケース A と比較し、ケース B のほうが通信相手は広範囲に広がっていることから、通信コストが増大している。一方、オーバーセット通信では、ケース A で特定の領域に多数の領域と通信する領域が局在している傾向が、ケース B で改善されていることから、通信性能の改善が期待できる。
上記の 2 ケースの計算時間はケース A が 1.230 sec/step、ケース B が 0.957 sec/step であった。なお、
速度発散場の最大値を 2 ケースで比較し、ほぼ一致することより、分割方法によって計算結果が変わらないことを確認済みである
本ベンチマークテストでは、グラフの一体化により、計算が 20%高速化されていることが確認された。
さらに、分割数を増やすと通信負荷が増すので、グラフ一体化の効果が大きくなることが期待される。

Fig.3.2-1. Bridge elements that connect oversets when graphing(a)隣接通信- 40 -

(b)オーバーセット通信
Fig.3.2-2. Comparison of communication patterns (left: case A, right: case B)(3-3) ADVENTURE_Solid3(3-3-1) ADVENTURE_Solid 3.0FS の概要
ADVENTURE_Solid は階層型領域分割法（HDDM）をベースとした陰的な有限要素法に基づく並列
構造解析コードであり、オープンソースの CAE システムである ADVENTURE システムに含まれる解析
モジュールのひとつである[1]。非公開の ADVENTURE_Solid FS 版はオープンソース版の機能を強化したものである。
「富岳」成果創出加速プログラムではこの FS 版の機能を強化しており、最終的に完成した版を ADVENTURE_Solid3 FS と呼ぶ。機能強化の内容は、構成則等の解析機能の追加と、48 コアのCPU を持つ約 16 万個の計算ノードを 6 次元メッシュ／トーラスネット―ワークである TofuD インターコネクトで結合した構成の「富岳」向けの線形ソルバー部分の改良である。その概要を以下に示す。
(3-3-2)解析機能(3-3-2-1) 風車ブレードの構造解析向けの解析機能
複合材料積層板のモデル化のために積層要素を開発した[2]。積層要素の各層は直交異方性体であり、弱形式の積分については、
層内方向は Gauss 求積、積層方向は合成 Simpson 公式による求積としている。
Fig.3.3-1 に積層要素の概念図、Fig.3.3-2 に風車ブレードの解析モデルを示す。
風車のモデル化には、流体力、回転する重力、遠心力を入力できるようにした。風車ブレードモデル
を実際に回転させるのではなく、Fig.3.3-3 のように力学的に等価な回転しない風車ブレードモデルを考えた。流体力、回転する重力、遠心力を境界条件・物体力として入力する機能を開発した。
大型風車ブレードの累積疲労損傷解析のために疲労損傷による剛性低下を考慮した解析が必要となる。
そこで、複合材料構造物への適用の実績がある Progressive Damage Model を実装した。Fig.3.3-4 に解析例を示す。

- 41 -

Fig.3.3-1. Layered element

Fig.3.3-2. Analysis model of wind turbine blade

Fig.3.3-3. Modeling of centrifugal force- 42 -

Fig.3.3-4. Progressive damage model(3-3-2-2) 石炭ガス化炉の構造解析向けの解析機能
石炭ガス化炉のような高温環境下で使用される構造物で使われる材料では、弾塑性挙動とクリープ挙動の両方が発生する可能性がある。そこで、移動硬化熱弾塑性（Ohno-Wang モデル）とクリープ（Norton-Bailey モデル）を組み合わせた複合構成則を開発した。さらに、これを ADVENTURE_Solid
コードに実装するための分離型応力積分手法に基づくフレームワークを開発した[3]。物性値の温度依存にも対応している。また、熱伝導解析コード ADVENTURE_Thermal の非定常熱伝導解析温度データを引き継ぐことで構造・熱連成解析を行うことができる。すなわち、温度場に基づき熱ひずみと温度依存性のある物性値を計算する。高温環境下で使われることを想定しているため、非弾性ひずみによる発熱は無視できると考え、考慮していない。そのため、熱伝導解析の後に構造解析を行う片連成解析となっている。Fig.3.3-5 に入力データの構造を示す。Fig.3.3-6 に設定できるオプションを示す。
Fig.3.3-7 にスラスタモデルの四面体メッシュを示し、Fig.3.3-8 に境界条件を示す。ここでは 3 種類の材料パラメータ（material #1, #2, #3）を設定している。Fig.3.3-9 に時刻 50 s における温度分布を示し、Fig.3.3-10 に Fig.3.3-8 の点 A と B における温度の時刻歴を示す。Fig.3.3-11 に一定の相当応力の値を593.5 MPa としたときの試験結果を再現する実験式、及び、3 種類のパラメータに対して Norton-Bailey
モデルにより求めた相当クリープひずみの時刻歴を示す。Fig.3.3-12 に各パラメータに対して得られた相当応力の Fig.3.3-8 の点 C から A の線上における分布を示す。Ohno-Wang モデルのみでクリープを考慮しないケース、material #1～#3 に対してクリープ（Norton-Bailey モデル）のみとしたケースとOhno-Wang モデルと組み合わせたモデル（複合モデル）によるケースの合計 7 ケースについてプロットしている。material #1 については、複合モデルと Ohno-Wang モデルの結果がほぼ一致しており、Norton-Bailey モデルの結果と異なっている。そのため、material #1 の複合モデルによる結果は塑性優
位な結果であるといえる。一方、material #3 については、複合モデルと Norton-Bailey モデルの結果がほぼ一致しており、Ohno-Wang モデルの結果とは異なっている。そのため、material #3 の複合モデル
による結果はクリープ優位な結果であるといえる。material #2 については、Ohno-Wang モデル、複合
モデル、Norton-Bailey モデルの結果は全て異なっており、複合モデルでは塑性変形とクリープ変形が同時に起こっていることがわかる。

- 43 -

移動硬化パラメータ設定

OWモデル設定:###############################################
OhnoWangHardeningTable 8
OWモデル移動硬化パラメータテーブル################################################
fega_type=Voidformat=f8f8f8
温度テーブルを使用する(1)table_id=0
かしない(0)かを指定する
depend_temperature=1
numOWmode=4
numTempTbl=2
300.0 1110.0 115.0#温度1 ζ1 r1
300.0 370.0 57.4
#温度1 ζ2 r2
300.0 200.0 24.3
#温度1 ζ3 r3
300.0 125.0 13.
#温度1 ζ4 r4
温度ごとに(ζ, r)の500.0 1110.0 115.0 #温度2 ζ1 r1
組を記述する500.0 370.0 57.4
#温度2 ζ2 r2
500.0 200.0 24.3
#温度2 ζ3 r3
500.0 125.0 13.3
#温度2 ζ4 r4################################################
OhnoWangCreep 1
OWモデル################################################
fega_type=AllElementConstant
テーブル参照format=f8
material_id=3
elastopl_table=0
owhard_table=0
テーブル参照creep_table=02.000E+01 #ref_temp

テーブル参照################################################
OhnoWangElastoplTable 3
OW弾塑性パラメータテーブル################################################
fega_type=Voidformat=f8f8f8f8f8f8table_id=0
depend_temperature=1
温度テーブルを使用する(1)かしない(0)かを指定する
300.0 1.25e5 0.30 1.0e3 1.0e9 0.0 #temp yng poi hp ys alpha400.0 1.25e5 0.30 1.0e3 1.0e9 1.0e-5 #temp yng poi hp ys alpha450.0 1.25e5 0.30 1.0e3 1.0e9 1.0e-5 #temp yng poi hp ys alpha

温度

ヤング率

ポアソン比

加工硬化

降伏応力

クリープ設定線膨張係数################################################
OhnoWangCreepTable 3
OWモデルクリープパラメータテーブル################################################
fega_type=Voidformat=f8f8f8f8
温度テーブルを使用する(1)table_id=0
かしない(0)かを指定する
depend_temperature=1
300.0 0.0 1.000E+00 0.0 #温度1 クリープA1 クリープN1 クリープM1
400.0 0.0 1.000E+00 0.0 #温度2 クリープA2 クリープN2 クリープM2
500.0 0.0 1.000E+00 0.0 #温度3 クリープA3 クリープN3 クリープM3

温度ごとにクリープパラメータ
A,n,mを指定する

Fig.3.3-5. Input data of combined elastic-plastic and creep model

熱弾塑性クリープ解析実行オプション

熱弾塑性クリープ解析出力オプション

メインオプション-creep-time-step numstep-time-result

サブオプション

--time-delta dt

メインオプション

内容
熱弾塑性クリープ解析を有効にする。
熱弾塑性クリープ解析を行う
タイムステップ数を指定する。
熱弾塑性クリープ解析時の時間刻みを
指定する。

出力物理量
サブオプション
(-time-result,-result共通)

サブオプション

内容-thermal

伝熱解析の結果引継ぎを有効にする。

-thermal-file filename

伝熱解析結果のADVファイル名を指定する。
伝熱解析結果のディレクトリを指定する。

-thermal-dir dirname

内容
各時間ステップの出力物理量を指定する。

-result

熱弾塑性クリープ解析 伝熱解析引継ぎオプション :ADVENTURE_Thermal の結果を引き継ぎ可能
メインオプション

サブオプション

最終ステップの出力物理量を指定する。
--owcrp-backstress

背応力(合計)を要素で出力する。

--owcrp-backstress-n

背応力(合計)を節点で出力する。

--owcrp-backstress-i

背応力(合計)を積分点ごとに出力する。

--eqplstra相当塑性ひずみを要素で出力する。

--eqplstra-n相当塑性ひずみを節点で出力する。

--eqplstra-i相当塑性ひずみを積分点で出力する。

--eq-crp-stra相当クリープひずみを要素で出力する。

--eq-crp-stra-n相当クリープひずみを節点で出力する。

--eq-crp-stra-i

相当クリープひずみを積分点ごとに出力する。

Fig.3.3-6. Options of the analysis code for combined elastic-plastic and creep model

Fig.3.3-7. Sectional view of the mesh of the thruster model- 44 -

Fig.3.3-8. Boundary conditions applied to the thruster model

Fig.3.3-9. Distributions of the temperature at

Fig.3.3-10. Transition of the temperature at

the end of the heating process (50 s)

points A and B

Fig.3.3-11. Equivalent creep strain evolution

Fig.3.3-12. Equivalent stress at the final time

under a constant equivalent stress of 593.5

step on the line from point C to point A

MPa(3-3-2-3) 鉄筋コンクリート構造の解析機能
風車や石炭ガス化炉を支える構造物には鉄筋コンクリート（RC）構造が使われる可能性があるため、防災科学技術研究所で開発された RC 向けの構成則[4]を組み込んだ。この構成則は拡張 Drucker-Prager
モデルをベースとしており、圧縮側の圧壊を損傷モデルで表し、引張クラックを主方向による座標系（クラック座標系）の応力とひずみの関係においてクラックひずみを導入することでモデル化している。引
張クラックモデルでは側圧の影響（ポアソン効果）も考慮している。Fig.3.3-13 に示す簡単なテストモデ
ルにより Verification を行った。その結果を Table 3.3-1 に示す。

- 45 -

1mm

1mm

5.0e-3mm
押し込みor
引っ張り

1mm

Z

Y
X

Z拘束
XYZ
拘束

Z拘束

YZ
拘束

Fig.3.3-13. One element for the verification of the concrete model
Table 3.3-1. Results of the analyses using the one element model for various cases
押し込み or 引っ張り

ケース

圧壊

クラック

NR 反復収束 要素応力値[MPa]

1

押し込み

なし

なし

〇

4.552E+01

2

押し込み

あり

なし

〇

8.604E+00

3

引っ張り

なし

あり

〇

4.540E+00

4

引っ張り

あり

あり

〇

3.317E+00

5

押し込み

あり

あり

〇

8.604E+00

6

引っ張り

なし

なし

〇

4.557E+00

7

押し込み

なし

あり

〇

4.552E+01

8

引っ張り

あり

なし

〇

3.329E+00(3-3-3) 線形ソルバーの「富岳」[5]向けチューニング[6](3-3-3-1) BDD 法
Mandel[7]により提案されたバランシング領域分割（Balancing Domain Decomposition; BDD）法は、
部分構造型領域分割法に基づく、線形問題に対する前処理付きの反復解法である。ここではその定式化の概要を示す。BDD 法の計算過程では元の線形問題よりも小規模ないくつかの線形問題を解く必要がある。実装方法を検討するためにこれらを整理する。
(1) BDD 法の定式化
有限要素法に基づき離散化された静的線形構造解析における釣り合い式は、次のような連立一次方程式となる。
Ku  p(1)

ここに、 K は剛性行列、 u は節点変位ベクトル、 p は節点荷重ベクトルである。
ここでは、領域分割法の中でも部分領域内部の自由度を消去する部分構造型領域分割法を用いる。ま- 46 -

ず、解析領域を N 個のオーバーラップのない部分領域に分割する。部分領域 k に対する剛性行列を K  k  、
k 

節点変位を u  k  、節点荷重を p としたとき、これらの行列、ベクトルを次式のように部分領域内部の自由度と部分領域間境界の自由度に分割する。
 K IIk kK     k 
 K BIuIk  pIk  K IBk  
k 
k 
u

p


、
、


 k  
k
k 
K BB
u B p B (2)

ここに、添字 I は部分領域内部の自由度を、添字 B は部分領域間境界上の自由度を表す。領域間境界(k )全体の自由度数を nB とする。また、部分領域 k に関係する領域間境界自由度の自由度数を nB とする。
部分領域の内部の自由度に関する釣り合いは次式のように各部分領域において独立に成立する。
K II  uI   K IB uB   pI 
k

k

k

k

k

(3)

一方、部分領域間境界自由度の釣り合いは次式で表される。
N

N K u
k 1

(k )
B

(k )
BI

(k )
I

N

N

k 1

k 1

k ) (k )  N (Bk ) K (BBu B   N (Bk ) p (Bk )

(4)

(k )

ここに、NB は部分領域 k に関する領域間境界自由度を全ての部分領域に対する領域間境界自由度に変換するブーリアン行列である。式(3)は次のように内部節点変位について解くことができる。

uIk   K IIk   p   K   u  
1

k

I

k
IB

k
B

(5)

式(5)を式(4)に代入すると次式のようにまとめられる。

SuB  pB(6)

ここに、S は Schur complement、uB は部分領域間境界上の全ての節点に対する節点変位ベクトル、p B
は静的縮約された部分領域間境界上の節点荷重ベクトルである。行列 S は nB  nB 型行列となる。部分領(k )域 k に対するローカル Schur complement を S とすると、全体の Schur complement である S は次式となる。

- 47 -NS   N (Bk ) S ( k ) N (Bk )T(7)

k)
S(k )  K(BB
 K(BIk ) K(IIk ) 1K(IBk )

(8)

k 1

ここに、

である。
式(6)は前処理付き共役勾配法で解くものとする。Mandel[7]が提案した BDD 法はこの問題に対する効
果的な前処理手法であり、ADVENTURE_Solid にも実装されている。その前処理行列は次式となる。
1
1
M BDD  I  RK C-1R T S  M NN I  SRK C-1R T (9) RK -1C R T1

ここに、 MNN は Neumann-Neumann 前処理のための前処理行列である。 R と R

Tはそれぞれマルチグリッド法の prolongation 行列と restriction 行列に相当する。固体、構造解析では行列 R は nB  6N 型行列となる。 K C は次式によって計算される。

KC  RTSR(10)

行列 K C は Coarse grid（コース）行列と呼ばれ、剛性行列 K を部分領域の剛体運動を利用して低次元で近似した行列となる。 K C を係数行列とした線形問題をコース（グリッド）問題と呼ぶ。コース問題を解いた結果を用いて CG 法の収束性を高める処理は、前処理の一種でありコースグリッド修正と呼ばれる。初期残差に対してコースグリッド修正をするとそれ以降はコース空間の成分が取り除かれるため、式(9)の第 1 項の最初のコースグリッド修正と最後の項を省略できる。このとき、前処理行列は、
1
1M BDD  I  RK -1C R T S  M NN

(11)

となる。
1

Neumann-Neumann 前処理を表す MNN は、
N

1
M NN
  N (Bk ) D ( k )  S ( k )  D ( k )T N (Bk )T


k 1

(12)

と表される。ここに、D( k ) は領域間境界自由度を共有する部分領域の数の逆数を対角成分とした対角行列であり、重み行列と呼ばれる。  S( k )  はローカル Schur complement の一般逆行列を表す。荻野[8]は一般逆行列を使う代わりに次式のようにペナルティ項を加えて正則化し、逆行列を使う方法を提案し- 48 -

た。

S   S(k ) 

(k )





 max diag  K ( k )  I

1

(13)

(2) 解くべき線形問題の整理
BDD 法の計算過程で解く必要がある線形問題を整理する。以下では、一般的なベクトルを a 、 b あるいはそれらに添字を加えた記号で表す。
CG 法 の 計 算 過 程 に お け る 残 差 ベ ク ト ル の 計 算 、 式 (10) の コ ー ス 行 列 の 作 成 、 式 (11) のNeumann-Neumann 前処理の後に S との積の計算がある。そのために、次式のようなローカル Schur(k )complement S とベクトルの積の計算が必要である。

b(Bk )  S(k )a(Bk )(14)

S(k) を陽に作る代わりに次式のように本来のディリクレ境界に加えて部分領域境界自由度もディリクレ境界として与えた問題（ディリクレ問題）を解く。ディリクレ境界の処理は、ディリクレ境界自由度を除く方法ではなく、次式のようにその自由度に関する係数行列の対角項に 1 を入れる方法で処理する。
K IIk  OO  u (I k )   K IBk  b (Bk ) I BB  u (Bk )   a (Bk ) (15)

ここに、 I BB は単位行列を表す。左辺の係数行列の逆行列を陽に求めるのではなく、線形ソルバーを使(k )う。次に、 bB を次のように求める。
u ( k ) kK BB   (I k ) u B 

b(Bk )  K BIk

(16)

一方、Neumann-Neumann 前処理では、次式のように式(13)の一般逆行列を近似した行列とベクトルの積の計算が必要である。



 a



b (Bk )  S ( k )  max diag  K ( k )  I

(k )

1

(k )
B

(17)

ここでも S を陽に作る代わりに次の方程式（ノイマン問題）を線形ソルバーにより解く。

- 49 -

K IIk  kK BIK IBk  a ( k )   0 
  I    (k ) kK BB  max diag  K ( k )  I  a (Bk )  b B 






(18)

-1

最後に、コースグリッド修正に含まれる K C とベクトルの積についても、逆行列を陽に作る代わりに線形ソルバーを用いて計算する。すなわち、b  K-1Ca(19)

KCb  a(20)を求めるために、を線形ソルバーにより解く。
以上のように、BDD 法の実装では式(15)、(18)、(20)の求解の実装方法がポイントとなる。式(15)、(18)は部分領域単位の局所的な計算であるため、ローカル問題と呼ぶ。式(20)の計算はコース問題と呼ぶ。
(3-3-3-2) ADVENTURE_Solid におけるバランシング領域分割法の従来の実装方法(1) 階層型領域分割法（HDDM）
ADVENTURE_Solid は Yagawa and Shioya[9]が提案した階層型領域分割法をベースとして並列実装をしている。Fig.3.3-14 のようにメッシュはまず Part と呼ばれる領域に分割され、各 Part がさらに部
分領域に分割される。階層型領域分割は領域分割ツール ADVENTURE_Metis で行う。その中で、Partへの分割には ParMetis[10]が、Part の部分領域への分割には Metis[11]が使われている。

Fig.3.3-14. Hierarchical domain decomposition(2) 従来の実装方法
ADVENTURE_Solid FS 版では式(15)、(18)のローカル問題は直接法の一種であるスカイライン法のシリアル処理版で解いている。一方、式(20)のコース問題はスパース直接法ソルバーMUMPS[12]で解いて- 50 -いる。ADVENTURE_Solid は MPI により並列化されており、1 個の Part を 1 個の MPI プロセスに割り当てている。マルチコア CPU を使う場合には、1 個の MPI プロセスの中で OpenMP によるスレッド
並列化を行っており、1 個の Part の中にある複数の部分領域に対する式(15)、(18)のローカル問題をコア数と同数のスレッドに割り当てることで並列に計算する。この並列化では理想的なスケーラビリティが得られている。一方、コース問題に使う MUMPS も MPI 並列化されているが、超並列計算機では Part数と同数の MPI プロセス数は非常に大きいため、Part の情報を MUMPS に適した数の MPI プロセスに集約し、少ない MPI プロセス数で並列計算を行っている。また、MUMPS で使われる BLAS をスレッド並列版とすることで複数のコアを使ったマルチスレッド並列化ができる。
(3) 従来版による原子力発電所モデルの解析「富岳」に実装した ADVENTURE_Solid FS 版の従来版により、四面体二次要素による福島第 1 原子力発電所 1 号機モデルの静的解析を「富岳」により行う。このモデルの自由度数は 1,497,322,665（約15 億）である。なお、「富岳」の 1 個の計算ノードは 4 個の CMG と呼ばれる単位により構成され、48個のコアとメインメモリの関係は均質ではない。ここでは 1 個の CMG に 1 個の MPI プロセスを割り当てて計算する。したがって、1 個の MPI プロセスには 12 個のコアが割り当てられる。
Table 3.3-2 に 5 種類の領域分割のケースを示す。ここでは Case A と B に対して従来版により解析を行った結果を示す。Table 3.3-3 にそれぞれのケースにおけるコース問題の規模とローカル問題の規模を示す。コース問題の規模は全部分領域数×6（剛体モードの数）となる。ローカル問題の規模は、全自由度数/全部分領域数で概算したものである。実際には部分領域間境界上の自由度が重複することや、式(15)のディリクレ問題よりも式(18)のノイマン問題の方が係数行列の非零成分が少ないという差異はあるものの、計算のオーダーの見積りにはなっている。
「富岳」上の従来版により計算性能に対するパラメータ
スタディをした結果、従来版では Case A の領域分割によって最も速い計算速度が得られている。
Fig.3.3-15 に CG 法の 1 反復ステップ（1CG）の計算におけるコース問題の計算時間とその他の計算
時間の割合を示す。Case A のローカル問題の数は Case B の 10 倍以上で、ローカル問題の規模は 1/10
以下となっている。ローカル問題の求解に用いているスカイライン法は自由度に対してスケーラブルな
解法ではないため、全てのローカル問題に対する計算時間は Case A の方が圧倒的に短い。しかし、コース問題の規模が大きいため、MUMPS によるコース問題の求解の時間が全体の半分近くを占めている。
MUMPS に割り当てる MPI プロセス数は 96～192 程度であり、全体の 6,144 個の MPI プロセスの中のほんの一部に過ぎないため、MUMPS の計算の間はほとんどの MPI プロセスは何もしていないことになる。一方、Case B では 1 個の Part の部分領域数を 1 個に減らすことでコース問題の規模を小さくしている。コース問題とローカル問題の規模は同程度である。ただし、コース問題は MUMPS で、ローカル
問題はスカイラインソルバーで解いているため、計算時間の大半はローカル問題の求解時間になっている。全体の計算速度は Case A よりも大幅に遅く、収束まで計算していない。
プロセスとスレッドの割り当て、及び、1CG および全体の計算時間と CG 法の反復回数は後の Table
3.3-4、5 にそれぞれ改良版の情報と共に示す。CG 法の収束判定では、後の改良版による解析も含めて相対残差の閾値を 1.0  10  3 としている。
(3-3-3-3) ローカル問題にも MUMPS を使った実装- 51 -

(1) 実装
3-3-3-2 (3)節の数値実験の結果から、部分領域数を調整することで、コース問題とローカル問題の規模を同程度にすることができ、多くの計算ノードを有効活用できる可能があることがわかった。しかし、その結果、ローカル問題の規模が大きくなり、スカイライン法を使うと計算時間が長くなる。そこで、ローカル問題にも MUMPS を使えるように解析コードを改良する。
Table 3.3-2. Cases of domain decomposition
Number of subdomains per

Total number ofa Partsubdomains

6,144

24

147,456

Case B

12,2881

12,288

Case C

6,144

4

24,576

Case D12,2882

24,576

Case E

24,0001

24,000

Cases

Number of Parts

Case A

Table 3.3-3. Size of the coarse grid problem and average size of the local problem for each caseApproximate size of the coarse gridApproximate and average size of the

problem (DOFs)

local problem (DOFs)

Case A

880,000

10,000Case B

73,000122,000

Case C

150,000

61,000Case D

150,000

61,000Case E

140,000

63,000Cases

6144 parts, 24 subdomains/partTime for MUMPS
Time for local solver and other(a) Case A- 52 -

12288 parts, 1 subdomains/partTime for MUMPS
Time for local solver and other(b) Case B
Fig.3.3-15. The percentages of computation time for the coarse grid problem and of the total
computation time for local problems and miscellaneous procedures using the previous version
1 個の部分領域に関するローカル問題を MPI により分散並列化されている MUMPS で解く場合、1 個の 部 分 領 域 に 複 数 の MPI プ ロ セ ス （ 以 下 、 プ ロ セ ス ） を 割 り 当 て る 必 要 が あ る 。 従 来 の
ADVENTURE_Solid では、Fig.3.3-16 のように部分領域を管理する Part に 1 個のプロセスを割り当てていたため、プロセスの割り当て方法を見直す必要がある。
MPI では集団的操作の対象となるプロセスの集まりをコミュニケータとして定義する。Fig.3.3-17 に改良した ADVENTURE_Solid（改良版）における MPI のコミュニケータの配置を示す。従来の実装では Part 数と同数のプロセスで構成される Master Group のみであったのに対し、今回の実装では Worker
Group を加える。MPI_COMM_WORLD は全プロセスを表すコミュニケータである。PLU コースグリ
ッド MUMPS コミュニケータは、コースグリッド問題の求解に使う MUMPS 用である。MPC MUMPS
コミュニケータは、MPC を含む問題において拘束条件に関する射影を行う MUMPS 用のコミュニケータである。このコミュニケータは複数個配置することもある。以上のコミュニケータは Master Groupの中に置く。MUMPS ローカルコミュニケータは、ローカル問題の求解に使う MUMPS 用であり、Partの数と同じ数になるように生成する。コミュニケータに 2 個以上のプロセスが含まれる場合には、Master
Group のプロセスだけでは足らないので、Worker Group のプロセスを使う。1 個の Part に複数の部分領域がある場合には、MUMPS のオブジェクトは部分領域毎に生成するが、計算は順番に行うものとしてコミュニケータは同じもの（各 Part に一つ配置）を使う。
(2) 計算性能評価
Table 3.3-2 の Case C, D, E について、改良版で解析を行う。Table 3.3-4 に各ケースのプロセスとス
レッドの割り当てを、Table 3.3-5 に計算時間と CG 法の反復回数を示す。改良版ではローカル問題に対して多くの MPI プロセスが必要なこと、
また、
スレッド並列が部分領域毎の処理に対するものではなく、スレッド並列版の BLAS に対するものであるため、並列性能が低いことから、1 個の計算ノードに 8 個の MPI プロセスを割り当てている。Case C, D, E の全部分領域数はほぼ同じであり、違いは 1 個の Partあたりの部分領域数である。
Case C と D について、Case C の CG 法の反復数が Case D と同じであると仮定した上で並列化効率- 53 -

を求めると、66.4%となり良好な結果となっている。1CG ループの計算時間で比較すると、Case E では
さらに高速になっている。しかし、Case E では収束までの CG 法の反復回数が多く、全計算時間は従来
版の最速値である Case A の計算時間よりも遅くなっている。BDD 法の性質により、部分領域数が少な
くなるとコースグリッド修正の性能が低下して反復回数が増大する。それだけでなく、1 個の Part あたりの部分領域数が 1 個の場合、収束性がより低下する現象がみられる。部分領域数を少し変えたケース
でも同様の傾向となった。Part あたりの部分領域数が 1 の場合、領域分割は ParMetis のみで行うことになるので、ParMetis による領域分割の品質が、Metis によるものと比べてやや劣っている可能性が考えられる。
(3-3-3-4) まとめ
ここではまず、階層型領域分割法に基づく並列有限要素構造解析コード ADVENTURE_Solid（ver. 2）
を「富岳」で実行した場合の問題点について分析した。次に、数万個の計算ノードを用いても十分な並列化効率が得られるように、ローカル問題をスパース直接法ソルバーMUMP で解くような実装を行い（ADVENTURE_Solid（ver. 3））、「富岳」におけるその計算性能を評価した。

MPI_COMM_WORLDParent Group(0 .. npart-1)
PLU コースグリッド
MUMPSコミュニケータMPC
MUMPSコミュニケータ

Fig.3.3-16. Allocation of communicator in the previous version
MPI_COMM_WORLD
Master Group(0 .. npart-1)

Worker Group(npart .. nproc-1)

Parent Group(0 .. npart-1)
PLU コースグリッド
MUMPSコミュニケータMPC
MUMPSコミュニケータMUMPSローカル

MUMPSローカル

コミュニケータ1

コミュニケータ1

MUMPSローカル

MUMPSローカル

コミュニケータ2

コミュニケータ2

Fig.3.3-17. Allocation of communicator in the improved version
Table 3.3-4. Assignment of processes and threads (Cases A and B are solved by the previous version、and Cases C, D, and E are solved by the improved version）
Cases

Number of nodes

Number of

Number of MPI

Number of MPI

Number of

MPI processes

processes perprocesses for the

threads per MPI- 54 -

subdomaincoarse gridprocess

problemCase A

1,536

6,144

1

192

12

Case B

3,072

12,2881

16

12Case C

3,072

24,576424

6Case D

6,144

49,152424

6Case E

12,00096,000

4

24

6

Table 3.3-5. Computation time and the number of iterative steps of the CG method
Cases

Computation time for 1CG

Number of iterative

Total computation time（s）（s）

steps of the CG method

Case A

0.55

678

461

Case B

3.72－－

Case C

0.86

1,260

1,193

Case D

0.53

1,024

594

Case E

0.38

1,321

524

(3-3-4) ローカル問題ソルバー向け並列スカイラン法
ADVENTURE_Solid の「富岳」向けチューニングを令和 3 年度より継続して実施した。特に BDD 法のローカル問題ソルバーにおいてスカイライン法の並列化技法を検討し、疎行列データのブロック化と短冊(strip)化を導入した（Fig.3.3-18）。短冊化の導入により、行列演算のスレッド並列化において並列の粒度を大きくとることが可能となり、SIMD 幅拡大など将来的な計算機への対応も可能となった。また、BDD 法のプロセス内スレッド並列化アルゴリズムとして部分領域単位並列化と行列内のブロック／
短冊単位並列化のハイブリッドが可能となり、階層性を持つコア／メモリなど将来的な計算機への対応も可能となった。
また、BDD 法におけるコース問題ソルバーにおいて並列スパース直接法ソルバーMUMPS を長らく利用してきたが、BDD 法では前進後退代入の効率向上が重要となるものの、MUMPS はチューニングが困難という問題があった。そこで、コース問題ソルバーとして前述した並列スカイライン法も選択できる実装を行った。簡易形状モデルを用いた性能評価（Table 3.3-6）では、数万自由度規模のコース問題における前進後退代入が数十ミリ秒、数十万自由度規模では数百ミリ秒となり、十分な実用性を持つことが示された。
さらに、これまでの「富岳」向けチューニングに関する研究成果を並列スカイライン法ライブラリとして整備した。ADVENTURE_Solid FS 版を同ライブラリと結合し、ローカル問題ソルバーを並列スカ
イライン法ソルバに置き換えて約 15 億自由度規模モデル解析の性能評価を行った（Table 3.3-7）。なお、
開発ライブラリはメモリ使用量に関する最適化が十分でないこともあり、コース問題には並列スカイラ
イン法を適用せず MUMPS を用いた。表より、従来実装でも解ける問題と計算機規模の場合には高速化の実現には至らなかったが、わずかに遅い程度の性能までは達成できた。

- 55 -

Fig.3.3-18. Strip/block mining for the skyline matrix
Table 3.3-6. Performance evaluation of parallel skyline solver on A64FX
Number of

DOFs per

BDD

subdomains subdomain

Number ofComputational time [sec]iterations computing
nodesTime per BDD
iteration [msec]

Factorization

Solution

Local

Coarseproblemproblem

6,144

160×103

1,040

1,536

48

204

185

11

8,192

120×103

897

2,048

63

112

109

15

16,384

3

61×10

610

4,096

33

59

59

38

32,76831×103

433

4,096

338

60

48

90

65,53615×103

313

2,048

148

100

89

232

Table 3.3-7. Performance evaluation of parallel skyline solver in solving 1.5 billion DOFs model on
1,536 computing nodes
Local problem

ThreadCoarseComputational time [sec]solverparallelization

problem

Factorization

SolutionsolverSkylineSubdomain-wise

MUMPS

56

1,122

Block skyline

Subdomain-wise

MUMPS

154

1,489

Block skyline

Block-wise

MUMPS

537

6,122

参考文献
[1] Home page of ADVENTURE Project, https://adventure.sys.t.u-tokyo.ac.jp[2] 遊佐泰紀、宮村倫司、山田知典、吉村忍：領域分割法を用いた風車ブレード構造の並列有限要素解
析、計算工学講演会論文集、Vol. 23、D-09-02、2018.
[3] Yasunori Yusa, Tomoshi Miyamura, Jun Yin, Kuniaki Koike, Takashi Ikeda, and Tomonori Yamada,“Partitioned coupling framework to combine a kinematic hardening plasticity model and a creep model for
structures in a high-temperature environment,” International Journal for Numerical Methods in Engineering, Vol.

- 56 -

123, Issue24 (30 December 2022), pp. 6170-6196 (DOI: 10.1002/nme.7108) (Published online: 26 August 2022)[4] 山下拓三、藤原淳、宮村倫司、「ポアソン効果を考慮した引張クラックモデルの開発と RC 構造物の大規模有限要素解析への適用」、日本計算工学会論文集、Vol. 2022（2022 年 7 月 28 日）、
Paper No.20220006(DOI: 10.11421/jsces.2022.20220006)[5] Web site of Fugaku, https://www.r-ccs.riken.jp/en/fugaku/[6] 宮村倫司、小池邦昭、吉村忍、「並列有限要素構造解析コード ADVENTURE_Solid の富岳向けチューニング」、計算工学講演会論文集、Vol. 28（2023 年 5 月 31 日－6 月 2 日（発表 2 日））
、つくば（つくば国際会議場）[7] Mandel, J: Balancing domain decomposition, Communications on Num. Methods in Eng. Vol. 9, pp. 233-241,
1993.
[8] 荻野正雄、塩谷隆二、金山寛、田上大助、吉村忍：バランシング領域分割法による並列弾性有限要
素解析、日本機械学会論文集 A 編、Vol.69, No.685, pp.1360-1367, 2003.
[9] Yagawa, G., Shioya, R.: Parallel finite elements on a massively parallel computer with domain
decomposition, Comput. Syst. Engng., Vol. 4, pp. 495–503, 1993.
[10] Karypis, G, Kumar, V: Parallel Multilevel K-Way Partitioning Scheme for Irregular Graphs,

Tech. Rep., Dept. of Comp. Sci. Univ. of Minnesota, TR 96-036, 1996.
[11]Karypis. G, Kumar, V: Multilevel K-Way Partitioning Scheme for Irregular Graphs, Tech. Rep.,

Dept. of Comp. Sci. Univ. of Minnesota, TR 95-064, 1995.
[12] Home page of MUMPS: a parallel sparse direct solver, https://mumps-solver.org/(3-4) ADVENTURE_Thermal
令和 4 年度は ADVENTURE_Thermal の「富岳」向けチューニングの実施に向け、まず、事前検討として ADVENTURE_Solid プロトタイプ実装を用いて各アルゴリズム（BDD、Local Solver 等）の計測及び最適化を昨年度より継続して実施した。その結果 BDD に関しては、線形代数ソルバーの OpenMP
対応化を進め、また、コース行列のスカイライン記憶のブロック化と BLAS、LAPACK ライブラリの効率的な利用等により、億レベル自由度の大規模板曲げ問題においてピーク性能比 20%程度を実現した。
また、Local Solver に関しては、BDD と同様にこれまでの 1 領域 1 コアの割当て 1 領域複数コアに変更
し、それに伴い並列スカインソルバーを導入することにより、同じく、億レベル自由度の大規模板曲げ問題においてピーク性能比 15%を実現した。これらの成果を ADVENTURE_Solid2 に移植し、引続き「富岳」上での最適化を進めた。
ADVENTURE_Thermal の「富岳」向けチューニングの実施として、各アルゴリズム（BDD, Local
Solver 等）の計測をはじめとする検討を行った。また、ADVENTURE_Thermal の公開版に対してOpenMP を導入し、ハイブリッド並列化し、いくつかの簡単な数億要素規模のテスト問題を用いて動作検証を実施した。その結果、ADVENTURE_Thermal に OpenMP が正しく動作し OpenMP 導入の導入が完了した。令和 3 年度から問題となっていた、BDD を実行した際に通信エラーが発生するバグを取り除くデバッグを今年度も引続き実施したが、これについては問題解決には至らなかった。しかしながら、
ADVENTURE_Thermal の場合は、特に数億要素以上の大規模非定常解析は BDD 前処理の効果は限定的となるため、本研究においては目的を達成できたと考える。このほかに、DDM アルゴリズムのボトル- 57 -

ネックとなる箇所の性能向上に向け、ADVENTURE_Solid の成果も有効活用しつつプログラムをさらに改善し、ピーク性能比の向上を実現した。本研究の ADVENTURE_Thermal のチューニング開発の成果
は、広い分野への適用が期待できるものとなった。
(3-5) REVOCAP_Coupler
令和 4 年度は、令和 2 年度に「富岳」に移植したコードをもとに、「富岳」の MPI 仕様及びファイル
システムの制約に対するチューニングを行い、令和 3 年度から引き続きラボスケール石炭ガス化炉の燃
焼流－伝熱連成解析におけるオンライン連成機能、及び洋上ウィンドファームのスーパーシミュレーシ
ョンにおけるオフライン連成機能の両者ともに良好に動作することを確認した。
(3-6) RC-HPC 版
RC-HPC 版は産業分野での利用・展開も考慮し、「富岳」だけでなく、第２階層群の計算資源を始め多様な環境での実行を視野に入れて開発を進めてきた。本年度は、今後の導入が進むと予想される電力性
能あたりの性能が高い GPU 並列計算機での実行を想定した、マルチ GPU 並列版を開発し、精度検証を行った。
基礎方程式や解法などは従来と同様に LES によりモデル化された非圧縮 Navier-Stokes 方程式をFractional Step 法で解く。開発及び性能評価を実施した計算機環境は九州大学のスーパーコンピュータ
ITO サブシステムＢで、Fig.3.6-１に示すようにノードあたり４GPU を搭載する。並列計算時には計算領域を複数のサブ領域を分割し、１サブ領域を１プロセスに割り当てる。領域内の計算は GPU で行い、CPU は MPI 通信と GPU の制御を行う。GPU での並列化方法は OpenACC のディレクティブによる。
計算に最も時間を要する圧力の Poisson 方程式には Red-Black SOR 法を用いた。

Fig.3.6-1 GPU configuration on ITO subsystem B at Kyushu University.
並列計算においては、プログラムの並列化率が低いと逐次処理の時間が増え、ノード数が増加すると性能向上が大幅に阻害される。逐次処理の主要因は通信時間であるため、通信隠蔽を用いることにより、並列化効率を高めた。Fig.3.6-2 に示す通信アルゴリズムにより、計算処理を実行している間に通信を同時に行い、全体の実行時間を短縮する。

- 58 -

Fig.3.6-2 Overlap communication to avoid performance degradation.
開発したコードの性能を弱スケーリングと強スケーリングで評価した。まず、弱スケーリングの測定では、１プロセスあたり 500x350x150 セル（26M）を割り当て、分割は主流方向に１次元分割とした。
Fig.3.6-3 に弱スケーリング時のスピードアップ性能を示す。テストでは、4 ノード（16GPU）で約４億セルの規模の計算を想定した。測定結果からは良好なスケーラビリティが確認できる。このときの効率
は Fig.3.6-4 に示すように、95%以上を維持している。

Fig.3.6-3 Weak scaling performance.

Fig.3.6-4 Efficiency for weak scaling test.

- 59 -

Fig.3.6-5 および 3.6-6 には、それぞれ強スケーリング時のスピードアップと効率を示す。強スケーリ
ングテストでは、1GPU のメモリに収まる 800x 350 x 200(56M セル)を割り当てた。強スケーリングではプロセスあたりの計算量が減少し、相対的に通信時間の影響が大きくなる。このため、弱スケーリン
グに比べるとスピードアップも効率も低下しているが、文献[1]などと比較しても良い結果を得ている。

Fig.3.6-5 Strong scaling performance.

Fig.3.6-6 Efficiency for strong scaling test.
次に、開発した GPU 版をベースに、風車計算に必要なアクチュエータラインモデルおよびアクチュエ
ータディスクモデルを組み込み、予測精度を実験値および CPU 版の RIAM-COMPACT（RC）の結果[2]と比較・評価した。計算条件としては風車１基の場合と２基の場合を評価したが、ここでは２基の結果について示す。840x470x470（1.85 億）格子を用いて、無次元時刻 40 まで計算し、無次元時刻 25～40（風車は４回転）の区間で時間平均操作を行った。TSR(Tip Speed Ratio)は４とした。計算した瞬時値の可視化結果を Fig.3.6-7 に示す。食い違い配置された２基の風車の周りの流れ場を渦度の等値面で可視化している。Fig.3.6-8 には速度成分Ｕ（x 方向成分）の時間平均場と比較する測定位置を示している。
Fig.3.6-9 において、x/D=6 の線上でデータをサンプリングし、風洞実験の値と比較している。Simulation
No1+ No.2 の結果と風洞実験の結果を比較すると、良い一致を示していることがわかる。

- 60 -

Fig.3.6-7 Flow field around two windmills visualized by vorticity contour. Color means magnitude
of velocity component of x-direction.

Fig.3.6-8 Time averaged U-velocity field around two windmills.

Fig.3.6-9 Comparison of U-velocity component at measurement line x/D=6 in Fig.3.6-8.
以上示したように、本年度はマルチ GPU コードを開発し、性能評価を行った。その結果、弱スケーリ
ング・強スケーリングとも良好なスケーラビリティを確認した。また、GPU コードにアクチュエータラ
インモデル及びアクチュエータディスクモデルを実装し、実験結果と比較してよい一致を得た。

- 61 -

参考文献[1] Christoph Riesinger, Arash Bakhtiari, Martin Schreiber, Philipp Neumann, Hans-Joachim
Bungartz. Holistic Scalable Implementation Approach of the Lattice Boltzmann Method for
CPU/GPU Heterogeneous Clusters. 2017. Computation 5(4).
[2] Uchida, T. Effects of Inflow Shear on Wake Characteristics of Wind-Turbines over Flat Terrain.
Energies 2020, 13, 3745. https://doi.org/10.3390/en13143745.
(4) 個別シミュレーションと AI（機械学習）の連携活用(4-1) 石炭ガス化炉適用に向けた ADVENTURE_Thermal と AI 連携活用
令和 4 年度は、令和 3 年度に達成した成果[1]を「富岳」に実装し、動作確認を行った。
「富岳」に導入されている Python と PyTorch を利用し、実装及び学習時間の評価を行った。Pythonと、関連パッケージのバージョンは以下の通りである。
Python ver.3.8.2
mpi4py ver.3.0.3
pandas ver.1.2.2
numpy ver.1.19.0
scipy ver.1.5.2
h5py ver.2.8.0
fapp ver.1.0.0
PyTorch と、関連パッケージのバージョンは以下の通りである。
PyTorch ver.1.7.0
Horovod ver.0.20.3
oneDNN ver.2.1.0L01_aarch64
計算には、リソースグループ small を使用した[2]。使用ノード数を 8、1 ノード当たりのプロセス数を8 として計算を実行した。学習に用いたデータセットの概要を以下に示す。
学習のサンプルデータ数 : 1,098
検証のサンプルデータ数 : 80
バッチサイズ: 128

学習エポック数: 150

コード実行時のおおよその時間と利用メモリはそれぞれ 15m55s であり、ノードごとの最大使用メモ
リは 3586.1 MiB であった。使用した 8 ノードの、ノードごとのメモリ使用量は次の通りである。
Node #1 3586.1 MiB, Node #2 3409.5 MiB, Node #3 3481.2 MiB, Node #4 3331.3 MiB
Node #5 3360.9 MiB, Node #6 3419.4 MiB, Node #7 3472.2 MiB, Node #8 3343.9 MiB
同じデータセットおよび同じ学習エポックを NVIDIA A6000 で実行した場合に比べて約 5 倍の高速化が達成されている。また、現状で最も早い NVIDIA H100 と比べると予想される計算速度は 1/20 程度で- 62 -

ある。並列化による高速化で対応する方法が現実的な学習時間削減の方法である。
１つの解析データから複数のデータセットを生成し高精度に予測する手法を提示した。多くのデータセットを追加することで設計に資する温度場の予測は可能であるとの結論を得た。
参考文献[1] https://www.rist.or.jp/sc/sc_report.htm, https://www.rist.or.jp/sc/report/r03/f303_r03.pdf[2] https://www.hpci-office.jp/materials/r03b_fugaku_seminar_system.pdf(4-2) 洋上風力適用に向けた ADVENTURE_Solid と AI 連携活用
洋上風力適用のため、ADVENTURE_Solid により動的応力解析が実施されている。この計算は、複数の条件の異なる流体解析の結果から得られた翼表面に与えられた流体力に基づきブレードの変形を時間
ステップごとに行う必要がある。この計算コストは設計評価においては現実的ではなく、設計における評価は迅速に物理的な傾向を損なわない範囲の精度で予測できることが重要である。このような背景のもと、流体力からブレードの代表的な位置の変位及びねじりを予測するサロゲートモデルの構築を実施した。
(4-2-1) 予測器構築のための機械学習の特徴と対象の分析
Table 4.2-1 に示すように、物理現象を学習するための技術開発が進められている。機械学習においては、正則化技術が汎化性能を高めるための技術として分類されている。対象によって有効に働く場合や複数の効果的な組み合わせが存在したりと明確な指針があるわけではない。これらの正則化は、物理問
題を学習するように制約を与えているとも理解できる。そこで、データドリブンな考えで汎化性能を高める方法を弱い制約と呼ぶことができる。その理由は、現象を特定することなく多くの場合で適用可能な方法であるためである。一方で、現象を特定する、また、ある現象を表現するようなデータセットを意図的に生成することなどは、強い制約と呼ぶことができる。これらの観点は、マージナルゲインを追及する機械学習においては、技術の組み合わせや相乗効果の評価など正しく評価するための重要な指針になる。
Table 4.2-1 Weak and strong constraints for generalized prediction- 63 -

Algorithm
・Weight Decay
・Weight Sharing
Weakconstraints= ・Dropout(ensemble method)
regularization ・Batch Normalizationetc.
・Graph Neural Network

Strongconstraints(adding several constraints to represent
physical phenomenon, i.e. Boltzmann
Machine=Ising model, etc.)

・Physical Informed Neural
Networketc.

Data
・Data Augmentation(oversampling)
・Downsampling・Addition of interaction
terms (The additional factors
linearize between input and output
values by feature engineering)(4-2-2) 学習に用いた入力因子、予測量及びネットワーク構成
ブレードの表面を長手方向（y 軸方向）にスライスした 10 点×20 スライス＝200 点の物理量を与える
または評価する位置として指定した。Fig. 4.2-1 には全体の評価節点の分布を示す。また、予測対象は、
Gr#0, #5, #10 の位置の平均変位及びねじり角を予測対象とする。また、Fig. 4.2-2 は、学習に用いる入力因子としての流体力、変位の代表点を示す。また、入力因子には翼の回転位相も含む。これらは時系列として予測が進むため、予測する現時刻を t としたときに流体力に関しては、t, t-1, t-2 を因子とする。
また、変位に関しては変位を予測するため、t-1, t-2, t-3 を因子とする。回転角は流体力と同様に t, t-1, t-2の回転位相を因子とする。このような入力因子を採用する理由は、１つの大規模計算結果から可能な限
りデータセットを生成する必要があるためである。また、このような差分形式の入力であれば、現象に対するルールそのものの学習を強く促すような構成になっている。これらの入力因子の総数は、963 になる。入力規模が大きい場合は、CNN を用いて入力をベクトル列ではなくマトリックス上に適切に配置したほうが学習の効率が上がり、結果として予測精度の向上が期待できる。このことから、Fig. 4.2-3 に示
すように、入力因子をマトリックス状に配置した。出力は連続量で 4×3 要素のマトリックスとした。
Fig. 4.2-4 に、畳み込みニューラルネットワークの構成とそのハイパーパラメータを示す。画像認識では Max Pooling は有効であるが、回帰問題では Max Pooling が有効でないこともある。今回のデータでは、Pooling は行っていない。また、フィルタサイズも(2,2)と小さいものとし、隣接する値のみでフィル
タ処理を行うように設定した。

- 64 -

Fig. 4.2-1. Nodal groups for prediction: #0, #5 and #10 are candidates for relative displacement of ablade.

Fig. 4.2-2. Nodes for input factors which are composed as time sequence: 16 groups are chosen as
input factors.

- 65 -

Range of time step n : 20 – 1219 (1200 steps)
min=18 or 21
Time stepN,t=n-1

0

Rotationangle ω,t=n-1

Time step Rotation
angle ω,N,t=n-2
t=n-2

Gr# 0

Gr# 0

Gr# 0

Gr# 0

Gr# 0

Gr# 0

Gr# 0

averaged

averaged

averaged

averaged

averaged

averaged

averagedf, t=nf, t=n-1

f, t=n-2θ, t=n-1,n-2, n-3u, t=n-1u, t=n-2u, t=n-3

…

…

…

…

…

…

…

Gr# 15

Gr# 15

Gr# 15

Gr# 15

Gr# 15

Gr# 15

Gr# 15averaged

averaged

averaged

averaged

averaged

averaged

averagedf, t=nf, t=n-1

f, t=n-2θ, t=n-1,n-2, n-3u, t=n-1u, t=n-2u, t=n-3

mout=3
nout=3 or 4

nin=17

Time step Rotation
angle ω,N,t=n
t=n

Gr#0 averaged u
Gr#5 averaged u
Gr#10 averaged u
Rot. Ang. θi: Gr#0, 5, 10

Fig. 4.2-3. Input and output parameters for CNN: this illustration is for one datasets and 1,200
datasets are generated from 20th to 1219th time step.

Fig. 4.2-4. Configuration of a network and hyper parameters for convolutional neural network.
(4-2-3) 予測対象及び結果
ケース１として比較的短い区間（20 タイムステップ）の予測を行った。ケース２として、長い区間（40
タイムステップ）の予測も実施した。時系列データはブレードの１つのケースしかないため、Fig. 4.2-5
及び Fig. 4.2-6 に示す予測対象は、学習データから取り除き、未学習データとして学習を行った。ケース
１及びケース２ともに補間の予測は可能であった。Fig. 4.2-7 と 4-2-8 に示すとおり、それぞれの平均誤差は、0.519%と 0.641%であった。
ケース３として補外領域の補間能力の評価を実施した。最終 80 ステップを検証データとして学習デー
タには含めない。x 方向変位と z 方向変位をそれぞれ Fig. 4.2-9 と 4-2-10 に示す。太い実線で示された部分が検証データとなる。Fig. 4.2-11 に予測結果を示す。補外領域の予測においても平均 1.16%の誤差- 66 -

で予測した。この領域においては、取りうる最大値そのものも学習データに含まれておらずルールに相当する関係を内部に構築していることを期待できる。その証明には、他の評価技術が必要であることとより多くの学習データによる予測結果の評価が必要である。

Fig. 4.2-5. Interpolation validation for short time step ranges.

Fig. 4.2-6. Interpolation validation for long time step range.

- 67 -

Fig. 4.2-7. Prediction result of case 1: short time step range.

Fig. 4.2-8. Prediction result of case 2: long time step range.

- 68 -

Fig. 4.2-9. Extrapolation validation of displacement in x direction.

Fig. 4.2-10. Extrapolation validation of displacement in z direction.

- 69 -

Fig. 4.2-11. Prediction result of displacement of extrapolation during last 80 time steps.
ブレードはたわみと同時にねじりも生じる。平均ねじり角を Fig. 4.2-12 に示すように求め、その値を表が Gr#0, 5, 10 において予測を行った。本予測は、補外領域で評価した。ねじりは時間とともに増加す
る傾向は変位と同じであるがピーク近傍で振動を示している。この様子を Fig. 4.2-13 に示す。変位同様に平均 1.16%の誤差で予測した。予測結果を Fig. 4.2-14 に示す。Group#10 の最終ステップ近傍では若干の差異は認められるが、おおむね許容される精度で予測ができたと判断する。

Fig. 4.2-12. Definition of average torsional angle by several measurement points.

- 70 -

Fig. 4.2-13. Torsional angles of group #0, #5 and #10: gradual transition according to increment oftime step.

Fig. 4.2-14. Prediction result of extrapolation of torsional angle during last 80 time steps.

- 71 -

(4-2-4) まとめ
本検討により以下の結論を得た。
・時系列データにおいて、後退差分的な対応で補間及び一部未学習データの予測に成功
・形状に関わるデータも入力に入れることにより汎用性の獲得は可能だが、本データでは異なる翼形状データがないため、変位に幾何的な意味を持たせた・形状は断面二次モーメント等の幾何量を積極区的に因子・1 ケースのみのデータにおいても、未学習変位、ねじれ角の予測成功＝本物理現象の予測の可能性高・学習データのコントロールによりサロゲート構築
・粗解析（パラメータ空間の拡大）と詳細解析（精度向上）
時系列データを用いて学習データ数を増加させ補外領域までの予測を可能とする機械学習の方法を提示しその可能性について示した。多くのデータを用いることで疲労設計に資するサロゲートモデルの構築が可能であることを示した。
(4-3) 後流分布モデル化に向けた FFB と AI 連携活用
洋上風車用の大型風車の解析は「富岳」を用いても長時間の解析時間が必要であるため、解析結果をもとにモデル化を行い、設計等に活用しやすくする必要がある。従来より Wake モデルとして後流の性質をもとに風車下流の流れを予測することが行われてきているが、本研究では、流体力学モデルの未知項、実験定数に相当する部分を解析結果をもとに推定し、汎用性の高いモデルの構築を試みた。Fig.4.3-1に風車下流の流れ場の模式図を示す。風車下流には速度欠損が生じるが、その形状は風車近傍と風車から十分離れた位置では異なる。従来は風車から十分離れた位置での速度欠損をモデル化することで予測が行われていたが、本研究では、風車近傍の流れ場を基本データとし、下流への発達過程、水平方向、
鉛直方向の速度分布の違い（鉛直方向の場合、海面があることにより上下非対称な流れとなる）を考慮したモデルを構築し、モデルのなかに出てくるいくつかの係数（従来は実験定数や経験定数として処理されていた）を長時間の LES 解析により得られたデータより推定した、Fig.4.3-2 に後流速度分布の変化の模式図を示す。本研究では、風車近傍では後流の速度欠損が二つの山を持つ 2 重ガウシアン分布となる（ダブルガウシアンモデル）を使用した。風車から十分離れた位置では１つ山のガウシアン分布となることから、2 つのガウシアン分布が一つになるまでの過程をモデル化する必要がある。この仮定を後流が線形に広がっていくと仮定してモデル化した。
Fig.4.3-3 に水平方向と鉛直方向の主流方向速度の変化の様子を示す。これらのデータを元に後流のモデルを構築した。
Fig.4.3-4 に開発した高精度 Wake 分布モデルを実験結果と比較し、誤差評価を行なった結果を示す。
提案モデルは従来モデルよりも高い精度で後流を予測することができることがわかる。特に風車近傍については従来のモデルよりも誤差が 1/3 程度となっている。ここで Full Wake は洋上風車の風車間隔として一般的な風車直径の 6 倍から 10 倍の範囲の後流を示している。この範囲における後流予測の精度が高いことが風車設置位置を決めるうえで重要であることから、本研究によって得られた高精度 Wake 分
布モデルは従来のモデルと比較して高い精度を有することがわかる。
Fig.4.3-5 は 3 台の風車解析に本解析モデルを適用した結果である。▲印が高精度 Wake 分布モデルに- 72 -

よる推定結果、実線が LES 解析による結果である。3 台目の風車に流入する速度分布を過大評価しているが、これは速度欠損を受けた 2 台目の風車の出力を過大評価したことによるものと考えられる。高精度 Wake 分布モデルは風車係数をもとに計算されるため、2 台目の風車の出力係数が速度欠損や乱れの影
響を十分に反映できていないためと考えられる。この問題を解決するには下流側の風車の出力係数をモデル化する必要がある。LES 解析結果もとに下流側の風車の出力低下をモデル化し、高精度 Wake 分布
モデルに組み込むことが今後の課題である。

Fig. 4.3-1. Wake development behind a HAWT and its transition mechanism.

Fig. 4.3-2. Illustration of the wake development behind a HAWT using the DG analytical wakemodel.

- 73 -

Fig. 4.3-3. The mean wake velocity in the vertical and lateral directions.

Fig. 4.3-4. Average RMSE against lidar measurement within the near-wake, far-wake, and full-wakeregions at the incoming reference velocity of 11 m/s.

Fig. 4.3-5. A Geometric averaging of wake velocity downstream of the triple turbines.

- 74 -

洋上ウィンドファームの解析手法として、これまでの LES 解析結果をもとに風車下流の運動量欠損をモデル化し、解析負荷の小さな乱流解析手法である格子ボルツマン法に組み込むモデルを開発することを試みた。Fig.4.3-6 に解析結果の一例を示す。NREL5MW の風車を 9 台配置し、解析を行った結果である。風車下流の格子解像度は従来の LES と同程度である。解析時間は LES 解析に 1/30 以下であり、複数風車の解析を非常に高速に行うことが可能であることを確認した。風車自体は解析していないため、
風車下流の運動量欠損を LES 解析等を用いてモデル化する必要があるが、機械学習等を用いて風車近傍
流れのサロゲートモデルを構築することができれば、格子ボルツマン法と組み合わせてウィンドファー
ムの解析が高速に行えることを確認した。

Fig.4.3-6. A case study of an offshore wind farm analysis combining the lattice Boltzmann method
and velocity deficit data near a wind turbine.
(4-4) 燃焼モデル構築のための FFR-Comb と AI 連携活用
燃焼モデル構築のための FFR-Comb と AI 連携活用については、 パイロットスケール超臨界圧 CO2
ガスタービンの解析と V&V において取り組んでいるため、詳細は(1-2)節「パイロットスケール超臨界圧
CO2 ガスタービンと V&V」を参照されたい。
(4-5) 洋上風力適用に向けた RC HPC 版と AI 連携活用
洋上のウインドファームの特性を評価するためには、風速・風向などのパラメータを様々に変化させ、総合的な状況を勘案した評価が必要である。そのためには、多数の計算結果を利用することになり、１ケースあたりの計算時間を可能な限り短くする必要がある。LES に基づく計算は現時点では最も信頼性- 75 -

の高い結果が得られるが、そのためには大規模なシミュレーションが要求され、計算時間が長くなり、現実的ではない。そこで、高い信頼性をもつ少数の結果から、比較的小規模な計算結果を用いて軽い計
算負荷で妥当な結果を与えるサロゲートモデルの開発が指向されている。機械学習や回帰モデルなど
様々な方法が提案されているが、ここではモデルの根拠や妥当性が検討できる方程式の形式のモデルを探索する方法に取り組んでいる。Equation Discovery のアプローチの中でも、探索するパラメータ空間の範囲をできるだけ制限せずに多様な表現を試行できる、遺伝的プログラミング手法を採用し、データからデータを表現する方程式を見いだす方法論を開発している。これまで１次元の非線形バーガーズ方
程式については、そのデータから元の方程式を見いだすことができることがわかっている[1]。また、ある程度のノイズ耐性も期待できる。開発したプロトタイプシステムは、サロゲートモデルを探索できる可能性を示しているが、実際の三次元の大規模データからモデル方程式を発見するためには計算時間がかかることが問題であることがわかっている。これは、プロトタイプシステムの開発にあたって、効率を優先して Python で開発したためである。Python は実行速度が遅く、効率的な並列処理が難しい点が問題である。そこで、Python と同等な開発効率であり、かつ実行速度もＣ言語並みである Julia 言語を用いたリファクタリングに取り組んだ。Julia による方程式探査システムは、以前の Python によるシス
テムと同等な結果を示すことを確認した。また、並列化の実装にも取り組み、動作を確認している。
参考文献[1] K.Ono, and I. Koga, “Rediscovery of governing equations from simulation data using Genetic
Programming,”Transactions of JSCES, Paper No.20201004, 2020.
(5) 産業界等との連携及びコンソーシアムの構築(5-1) 燃焼システム用次世代 CAE コンソーシアムとの連携
燃焼システム用次世代 CAE コンソーシアムと本プロジェクト担当者らとの会合を 2022 年 4 月 22 日
及び 2022 年 11 月 25 日に、それぞれ理研（神戸）とオンラインのハイブリッドにて開催し、本プロジェ
クトで開発中のスーパーシミュレーション技術の紹介を行うとともに、燃焼系産業界のニーズの聴収を行い、その結果を精査し、本シミュレーション技術の開発に反映した。
(5-2) 洋上風力関係機関との連携とコンソーシアムの構築
令和 3 年度に開始した次世代 CAE を活用する風力産業コンソーシアム（Internet of Wind energy[IoW]
フォーラム（仮称））の立ち上げ作業を継続し、関係者らと本プロジェクト担当者らとの会合を開き、本プロジェクトで開発中のスーパーシミュレーション技術の紹介を行うとともに、風力産業界のニーズの聴収を行い、その結果を精査し、本シミュレーション技術の開発に反映した。また、IoW フォーラム（仮称）、日本海事協会、日本気象協会、日本風力発電協会との連携を通して、国内外への開発技術の普及体制を構築した。なお、IoW 構築の背景、想定される参加者、活動内容の案については、本プロジェクトの令和 3 年度の成果報告書※1 を参照されたい。
IoW 構築の取組の周知と意見収集を目的として、令和 4 年 7 月 27 日に第 1 回のセミナーを開催した※2。本セミナーは約 500 名の視聴申し込みを受け、好評のうちに終了した。セミナーの動画は 2022 年 9

月 1 日から 10 月 31 日までの期間にてオンデマンド配信された。その成果を受けて、さらに適切な周知- 76 -

期間と関係各所への調整を経てコンソーシアム設立を行うべく、日本海事協会との取り組みを続けることとした。

※1 https://www.rist.or.jp/sc/sc_report.htm
https://www.rist.or.jp/sc/report/r03/f303_r03.pdf※2 https://www.classnk.or.jp/hp/ja/hp_news.aspx?id=8102&type=press_release&layout=1(6) 成果発表、広報・アウトリーチ活動
本プロジェクト全体として「クリーンエネルギー「富岳」プロジェクトシンポジウム」（第 3 回、ポス
ト「京」重点課題⑥から数えると第 7 回）を東京大学山上会館大会議室とオンラインのハイブリッド形式にて開催し、プロジェクトのねらいと成果を発表し、社会への広報・普及を行った※1。カーボンニュ
ートラルに向けて「富岳」と本プロジェクトがどのように貢献できるか、また産業競争力強化という観点からのセクションを設けるなど、プロジェクト終了後の発展を見据えた企画としたことや、洋上風力発電に対する期待の高まりを受け、大学や研究機関のみならず産業界からも多数の参加者が集まった。
そのほか東京大学生産技術研究所革新的シミュレーション研究センター主催、本プロジェクト及び「航空機フライト試験を代替する近未来型設計技術の先導的実証研究」（課題代表機関：東北大学大学院工学研究科）の共催で、「第 6 回 HPC ものづくり統合ワークショップ」を生産技術研究所コンベンション
ホールとオンラインのハイブリッド形式にて開催した※2。
「「富岳」を利用した革新的流体性能予測技術の研究開発」（課題代表機関：東京大学生産技術研究所）との協働であり、「「富岳」成果創出加速プログ
ラム」の前身であるポスト「京」重点課題⑥と同⑧として活発に行っていた協働を継続し、事業を推進した。
※1 https://postk6.t.u-tokyo.ac.jp/event/symposium202301.html※2 https://www.fugaku-pj.iis.u-tokyo.ac.jp/event/220927_WSNo6_Top-1.html- 77 -

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

行事

2022 年 7 月 4、5、28 日

2022 年度 全体会

2022 年 9 月 27 日

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

2023 年 1 月 27 日

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

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

担当機関

（1）次世代火力発電システム及 東京都文京区本郷7-3-1 国立び洋上ウィンドファームのス

大学法人東京大学大学院工学系

ーパーシミュレーション

研究科

担当責任者
教授 吉村 忍（2）個別シミュレーションとAI（機械学習）の連携活用（3）総合的事業推進

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

C3棟 国立大学法人京都大学
福岡県福岡市西区元岡744 国教授 渡邊 裕章

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

教授 飯田 明由１番１ 国立大学法人豊橋技術科学大学大学院工学研究科
東京都目黒区駒場4-6-1 国立教授 加藤 千幸

大学法人東京大学生産技術研究所
福岡県春日市春日公園6-1 国准教授 内田 孝紀

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

教授 小野 謙二

別添１

学会発表等実績

1．学会誌・雑誌等における論文掲載（すべて査読あり）
No.

掲載した論文（発表題目）

発表者氏名

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

発表時期

国内
／国
際

1

A Scalable Parallel Partition
Tridiagonal
Many-Core

SolverandLow

T. Mitsuda, K. Ono

3

IEEE

for

Parallel

B/F

Processing

Processors
22022

Internationaland2022/6

国際

2022/6

国内

2022/8

国際

2022/8

国際

2022/8

国際

2022/8

国際

2022/8

国内

Distributed
Symposium

Workshops (IPDPSW)

移動音の移動方向知覚に関する

瀧澤哲, 工藤彰洋, 武日本シミュレーション学会論文誌,

研究 ー適応法による検討ー

居周Vol. 14, No. 2, pp.71-77

High-Frequency

KentoElectromagnetic

Field

Onaka

and

Amane Takei

Proceedings of The 41th JSST
AnnualConference:

Analysis in Microwave Oven

International

Conferenceonby Paralle Finite Element

Simulation

MethodFukuoka, Japan, Aug. 31- Sep.

Technology,

2, 2022
4

FE Analysis of Numerical

Shin-ichiro

Human Body Model with 100

Sugimoto,

Amane

AnnualMillion

in

Takei

Masao

International

–

Ogino

DOFs

Electromagnetic
Heat

Field

Conduction

andProceedings of The 41th JSSTSound

Conference

Simulation

CoupledonTechnology,

Fukuoka, Japan, Aug. 31- Sep.

Problem5

Conference:

2, 2022
field

experimental

evaluationofenvironment

AkihiroKudo

andAmane Takei

Proceedings of The 41th JSST
Annualwith ADVENTURE_Sound

International

Conference:
Conference

SimulationonTechnology,

Fukuoka, Japan, Aug. 31- Sep.
2, 2022
6

Developmentof

ParalleAmane Takei

Proceedings of The 41th JSST

Microwave analysis software:

AnnualADVENTURE_Fullwave

International

Conference:
Conference

SimulationonTechnology,

Fukuoka, Japan, Aug. 31- Sep.
2, 2024
7

AI と AR を用いた平面図の 3 次

岡谷夏実, 塩谷隆二,

日 本 計 算 工 学 会 論 文 集 , Vol.

元モデル化システムの開発 25

中林靖, 多田光利

2022, No. 20220009

Aug. 2022.

- 79 -8

LES study on the breakup

S. Wada, R. Kai, R.

Proceedings of the Combustion

mechanism of LOX core in

KuroseInstitute, in press

LES study of stabilization

Y. Hu, R. Kai, J.

Proceedings of the Combustion

mechanism in lifted ethanol

Wen, T. Murakami,

Institute, in press

spray flames

Y. Jiang, R. Kurose

LOX/GH2

2022/9国際

2022/9国際

2022/9国内

2022/11

国際

2022/12

国際

2022/12

国際

2023/1

国際

2023/2/19

国際

supercritical

combustion
9

10

Non-Intrusive

Polynomial

Chaos 法を用いた効率的な不確後藤聡太, 金子栄樹,

計 算工 学 会論 文 集 , Paper No.

武居周, 吉村忍

20220013

M. Nomura, A. Takei

Proceedings

かさ定量化法11

Accurate

and

Fast

of

IEEE

Electrostatic Analysis Using

CEFC2022, Denver, Colorado,

Mesh

Nov. 24-26, 2022

SmoothingandGeometric

Multi-Gridfor

Numerical

Human

Body

Model
12

LES/flamelet/ANN of oxy-fuel

H. Kasuya, Y. Iwai,

Applications in Energy andcombustion for a supercritical

M.

Combustion

CO2 power cycle

Morisawa, T. Nishiie,

Itoh,

Y.

Science,

12,International

Journalfor

Tomoshi Miyamura,

Numerical

Methods

inJun

KuniakiEngineering, Vol. 123, No. 24,
pp. 6170–6196

100083 (2022).

R. Kai, R. Kurose
13

Partitioned

coupling

frameworktocombinekinematicahardening

Yasunori

Yusa,

Yin,

plasticity model and a creep

Koike,

Takashimodel for structures in a

Ikeda,

Tomonori

high-temperature

Yamadaenvironment
14

Balancingdomaindecomposition

methodTomoshiMiyamura

for

and

an

Yoshimuralarge-scale

analysisofassembly

structurehavingmillionsofmultipoint

linearwake

expansionfunctionfortheShinobu

Computer Methods in Applied
Mechanics

andEngineering,Vol. 405 (15 February 2023)
115846 (29 pages)

constraints
15

A

double-Gaussian
wake model

analytical

Qidun Maulana Binu

EnergySoesanto,

Engineering

Tsukasa

Yoshinaga, Akiyoshi
Iida- 80 -

Scienceand16

Investigation of convergence

Kento

improvement to speed-up of

Amane Takei

theOhnaka,Proceedings of Compumag2023,

2023/5

国際

2023/5

国際

発表した時期国内

Kyoto, Japan, May. 22-26, 2023,

high-frequency

electromagnetic field analysis
solver17

An

Efficient

Amane

Proceedings of Compumag2023,

Takei,

Shigeki

Kyoto, Japan, May. 22-26, 2023

Large-Scale High-Frequency

Kaneko,ShinobuElectromagnetic

Yoshimura

Quantification

Uncertainty

Sota

Method

for

Field

Goto,

Problems

2．学会等における口頭・ポスター発表No.

発表した成果（発表題目、口頭・発表者氏名発表した場所（学会名等）

ポスター発表の別）

／国
際

1

2

並列幾何マルチグリッド法に基第 34 回「電磁力関連のダイナ

武居周、野村政宗

づ く高 速・ 高 精度 数値 人体 解

ミ ク ス 」 シ ン ポ ジ ウ ム

析、口頭（SEAD34）, 2022.

高周波電磁界解析コード：

第 34 回「電磁力関連のダイナ

大中健登、武居周ADVENTURE_FullWave に

ミ ク ス 」 シ ン ポ ジ ウ ム

よる電子レンジ内環境解析、口（SEAD34）, 2022.

2022/5/11

国内

2022/5/11

国内

2022/5/31

国内

2022/5/31

国内

2022/5/31

国内

2022/6/2

国内

2022/6/2

国内

頭
3

電磁界-熱伝導連成問題におけ杉本振一郎、 武居周、

第 27 回計算工学講演会論文

る数値人体モデルの有限要素

荻野正雄

集, 2022.

大規模高周波電磁波問題に対後藤聡太、 武居周、

第 27 回計算工学講演会論文する効率的な不確かさ定量化

金子栄樹、 吉村忍

集, 2022.

武居周、 野村政宗

第 27 回計算工学講演会論文

解析、口頭
4

法、口頭
5

高精度と高速性を両立 する幾集, 2022.

何マルチグリッド法に基づく人体内静電界解析、口頭
6

スカイラインソルバーにおける前

日本計算工学会 第 27 回計

河合浩志

進消去後退代入の OpenMP 並

算工学講演会列化、口頭
7

粒子充填向けマルチエージェン

日本計算工学会 第 27 回計荻野正雄

ト強化学習における学習手法の算工学講演会

検討、口頭- 81 -

89

燃 焼の 数値 シ ミ ュ レー シ ョ ン ：企業に役立つ応用数理手法

ザインを目指して、口頭（招待講の研究会 第 44 回技術セミ演）

ナー

ESRGAN-Based

Chenyue

Jiao,

The China Visualization

for

Chongke Bi, Lu Yang,

and

Large-Scale Volume Data,

Zhen Wang, Zijun Xia,

Conference 2022

oral

Kenji Ono

Development of Huge-Scale

Amane Takei

Microwave

11

日本応用数理学会 ものづくり

Society 5.0 時代のスマートデ

Visualization

10

黒瀬良一

VisualMethod(WCCM),

ADVENTURE_Fullwave,

Yokohama, Japan, Jul. 31-

oral

Aug. 5, 2022
Electrostatic

Field

Analysis with Unstructured
Numerical
Model

Human

Using

Geometric

2022/7/1

国際

2022/7/31

国際

2022/7/31

国際

2022/7/31

国際

2022/7/31

国際

2022/7/31

国際

Computational

software:

Fast

国内

Analytics

16th World Congress on

analysis2022/6/21

Masamune

Nomura,Amane Takei

17th World Congress on
Computational

Body

Method(WCCM),

Parallel

Yokohama, Japan, Jul. 31-

MultigridAug. 5, 2022

Method, oral
12

Effect of

Metabolic Heat

A M M Mukaddes,

18th World Congress on

Generation, Blood Perfusion

Kanta

Computationaland

Ryuji Shioya, Amane

Method(WCCM),

Takei

Yokohama, Japan, Jul. 31-

Environment

Temperature on the Body

Purkaysta,

Temperature-an Approach ofAug. 5, 2022

Finite Element Simulation,
oral
13

Performance evaluation of

Akihiro Kudo, Amane

19th World Congress onparallelwave-soundTakei

Computational

analysis

software:

Method(WCCM),

ADVENTURE_Sound, oral

Yokohama, Japan, Jul. 31Aug. 5, 2022

14

Fast

Electrostatic

Field

Analysis with Unstructured
Numerical
Model

Human

Using

Geometric

Masamune

Nomura,Amane Takei

20th World Congress on
Computational

Body

Method(WCCM),

Parallel

Yokohama, Japan, Jul. 31-

MultigridAug. 5, 2022

Method, oral

- 82 -

15

Uncertainty

Quantification

Sota

Using

Non-Intrusive

Amane

21st World Congress on

Takei,

ShigekiComputational

Polynomial Chaos Method

Kaneko,Shinobu

Method(WCCM),

forYoshimura

Large-Scale

Electromagnetic

Goto,

2022/7/31

国際

2022/7/31

国際

2022/7/31

国際

Yokohama, Japan, Jul. 31-

Wave

Aug. 5, 2022

Analysis, oral
16

FE Analysis of Numerical

Shin-ichiro Sugimoto,22nd World Congress on

Human Body Model with

Amane Takei, Masao

Computational

100

Ogino

Method(WCCM),

MillionDOFs

inHigh-Frequency

Yokohama, Japan, Jul. 31-

Electromagnetic Field - Heat

Aug. 5, 2022

Conduction

CoupledProblem, oral
17

Coupledsimulation

vibration

and

of

Masao

sound

Amane

Yokoyama,
Takei,

radiation of violin in large

Yoshidome,

space, oral

Yagawa23rd World Congress on

Ryo

Computational

Genki

Method(WCCM),
Yokohama, Japan, Jul. 31Aug. 5, 2022

18

BDD preconditioner for a
diagonal-scaled

Masao Ogino

WCCM-APCOM 2022

2022/7/31

国際

Hiroshi Kawai

WCCM-APCOM 2022

2022/7/31

国際

Zheng,

18th World Congress on

2022/7/31

国際

Kaori

Computational

2022/7/31

国際

2022/9/27

国内

Schur

complement system, oral
19

Efficient implementation of
skyline solver for many coreand GPU environment, oral

20

Motion Feature Analysis of

Hongjie

Noh Dance using ANN based

ChiekoMotion Capture Technology,

Harada, Ryuji Shioya

Kato,

oral

Method(WCCM),
Yokohama, Japan, Jul. 31Aug. 5, 2022

21

Developmentof

3D

NatsumiOkatani,

18th World Congress onvisualization system from

Ryuji Shioya, Yasuhi

Computational

2D Plane Figures, oral

Nakabayashi,

Method(WCCM),

Terutoshi Tada

Yokohama, Japan, Jul. 31Aug. 5, 202222

日本の洋上ウィンドファームの第 6 回 HPC ものづくり統合ワ

吉村忍

開発・運用における HPC の役

ークショップ、東京大学生産技割、口頭術研究所 コンベンションホール/オンライン- 83 -

23

高精度と高速性を両立した人体

野村政宗、 武居周

内静電界解析、口頭

電子情報通信学会エレクトロ2022/10/13

国内

2022/10/13

国内

2022/11/16

国内

2022/11/16

国内

2022/11/16

国内

2022/11/16

国内

2022/11/16

国内

2022/11/16

国内

2022/11/16

国内

2022/11/16

国内

2022/11/16

国内

ニクスソサエティ，エレクトロニ
クスシミュレーション研究会（信学技報），EST 研 R4 第 3 回,
EST, 2022.

24高周波電磁解析への PUFEM

武居周

電子情報通信学会エレクトロ導入についての基礎検討、口ニクスソサエティ，エレクトロニ

頭

クスシミュレーション研究会（信学技報），EST 研 R4 第 3 回,
EST, 2022.

25

弾塑性構成則とクリープ構成則

宮 村倫 司 、 遊 佐 泰 紀、

日本機械学会第 35 回計算を組み合わせるための分離連

殷 峻、 小池 邦 昭 、 池 田

力学講演会、オンライン

成フレームワークによる高温下貴、山田知典のスラスタノズルの解析、口頭
26

弾塑性構成則とクリープ構成則

宮 村倫 司 、 遊 佐 泰 紀、

日本機械学会計算力学部門

を組み合わせるための分離連

殷 峻、 小池 邦 昭 、 池 田

第 35 回計算力学講演会講演

成フレームワークによる高温下貴、山田知典論文集（2022 年 11 月 16－18
日（発表 16 日），オンラインのスラスタノズルの解析、口頭
27

28

領域分割法に基づくバイオリン

谷 翔太 、 工 藤 彰 洋 、 武

日本機械学会 第 35 回計算の大規模音響解析の基礎検

居 周、 横山 真 男 、 矢 川

力 学 講 演 会 （ CMD2022 ） ,

討、口頭

元基

2022.

室内環境の並列流体解析の基諸 井翔 、 杉 山 博 信 、 藤

日本機械学会 第 35 回計算

本検討、口頭元昭一、武居周

力 学 講 演 会 （ CMD2022 ） ,
2022.

29

高周波電磁界-熱伝導連成問題

杉本振一郎、 武居周、

日本機械学会 第 35 回計算における 2.7 億自由度数値人体

荻野正雄

力 学 講 演 会 （ CMD2022 ） ,2022.

モデルの有限要素解析、口頭
30

大規模並列マイクロ波解析、口日本機械学会 第 35 回計算

武居周

力 学 講 演 会 （ CMD2022 ） ,

頭

2022.
31

スカイラインソルバーにおける前

日本機械学会 第 35 回計算

河合浩志

進消去後退代入の性能ベンチ

力学講演会

マーク、口頭32

異種混合プロセッサにおける並日本機械学会 第 35 回計算

荻野正雄

列 FEM コードの性能評価、口

力学講演会

頭、口頭
33

Improve
distribution

theuneven

代晨、 増田 正人 、 鄭日本機械学会 第 35 回計算

of宏杰、 塩谷 隆二 、 中

力 学 講 演 会 （ CMD2022 ） ,- 84 -

microexpression recognition

林靖

2022.

劉思涵、 塩谷隆二、

日本機械学会 第 35 回計算

中林靖

力 学 講 演 会 （ CMD2022 ） ,data, oral
34

Human

body

structure

of

anime

extraction

characters based on deep

2022/11/17

国内

2022/11/18

国内

2022/12/1

国内

2022/12/3

国内

2022/12/3

国内

2022/12/21

国内

2023/1/26

国内

2023/1/27

国内

2022.

learning, oral
35

CulturalAnalyzationofMyanmar Ethnic Group by
the

Artificial

ココアウン、 増田正人、

日本機械学会 第 35 回計算

塩谷隆二、 中林靖

力 学 講 演 会 （ CMD2022 ） ,intelligence2022.

technology, oral
36

スパコン版 RIAM-COMPACT

内田 孝紀、小野 謙

第 44 回風力エネルギー利用の開発と風車ウエイクシミュレー

二、飯田 明由、吉村

シンポジウムションへの適用－ドイツ Alpha

忍、加藤 千幸、山出Ventus 洋上風力発電所の実測

吉伸、今村 博、植田

データの再現性について－、口祐子

頭・ポスター
37

38

AI と AR を用いた学習支援シス

岡谷夏実、 塩谷隆二、

日本機械学会 2022 年度技術

テムの開発、口頭

中林靖、 多田光利と社会部門講演会

油入変圧器異常自動診断シス

鄭宏杰、 塩谷隆二

日本機械学会 2022 年度技術テムのための GUI 開発、口頭
39

カーボンニュートラル社会実現と社会部門講演会
吉村忍

第２回 スーパーコンピュータに向けた富岳の役割 －洋上「富岳」成果創出加速プログラ

風力開発を例としてー、口頭

ム シンポジウム「富岳百景」（招待講演）―「富岳」3 年の軌跡―、オンライン

40

領域分割型並列有限要素法に

西嶋幹太、 後藤聡太、

電子情報通信学会エレクトロ

基づく電子回路基板の電磁界

村山敏夫、 武居周ニクスソサエティ，エレクトロニ

解析に関する基礎検討、口頭

クスシミュレーション研究会（信学技報），EST 研 R4 第 4 回,
EST, 2023

41

タンデム配置大型風車群の飯田明由「スーパーシミュレーションとLES 乱 流 解 析 と 大 型 風 車

AI を連携活用した実機クリー

Wake モデルの開発、口頭

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン- 85 -

42

洋上ウィンドファーム全体の内田孝紀「スーパーシミュレーションとLES 乱流解析、口頭

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

AI を連携活用した実機クリー
ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン

43

大型風車ブレードの流体構造

吉村忍「スーパーシミュレーションと連成振動と累積疲労損傷解AI を連携活用した実機クリー

析、口頭

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン

44

AI 連携による大型風車ブレード

和田義孝「スーパーシミュレーションとAI を連携活用した実機クリーのサロゲートモデル開発、口頭

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン
45

パイロットスケール超臨界 CO2

黒瀬良一「スーパーシミュレーションとタービン燃焼器の燃焼乱流解

AI を連携活用した実機クリー

析と AI 連携、口頭

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン

46

ベンチスケール石炭ガス化炉の渡邊裕章「スーパーシミュレーションとAI を連携活用した実機クリー

燃焼乱流解析、口頭

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン
47

ラボスケール石炭ガス化炉の燃吉村忍「スーパーシミュレーションと焼 ・ 乱 流・ 伝 熱・ 冷却 ・ 構造 解AI を連携活用した実機クリー

析、口頭

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン

48

並列乱流解析コード FFB、口頭

山出吉伸「スーパーシミュレーションとAI を連携活用した実機クリー- 86 -

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン
49

並 列 燃 焼 乱 流 コ ー ド

吉田隼也「スーパーシミュレーションとFFR-Comb、口頭

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

2023/1/27

国内

2023/2/10

国内

AI を連携活用した実機クリー
ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン

50

並 列 構 造 解 析 コ ー ド

山田知典「スーパーシミュレーションとADVENTURE_Solid、口頭AI を連携活用した実機クリー
ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン

51

ウィンドファーム乱流場解析コー
ド

ス

パ

コ

ン

小野健二

「スーパーシミュレーションとAI を連携活用した実機クリー

版

RIAM-COMPACT、口頭

ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン

52

並 列 連 成 解 析 カ プ ラ ー

吉村忍「スーパーシミュレーションと

REVOCAP_Coupler、口頭AI を連携活用した実機クリー
ンエネルギーシステムのデジ
タルツインの構築と活用」第 3
回シンポジウム、東京大学山上会館大会議室/オンライン53

領域分割型並列有限要素法に

溝口奈七子、武居周電子情報通信学会エレクトロ

基づく人体 fullwave 電磁界解

ニクスソサエティ，エレクトロニ

析、口頭

クスシミュレーション研究会（信学技報），EST 研 R4 第 4 回,
EST, 2024

54

富岳を用いた乱流解析の産業

2022 年度 産応協対話交流

飯田明由

利用事例（洋上風車と自動車空会セミナー

力音の解析事例）、口頭（共催：日本応用数理学会ものづくり研究会）
社会課題解決に向けた数理- 87 -

科学の進展、オンライン

55

Parallelwave

analysis

sound

based

hierarchicalondomainAmane Takei, Akihiro

ICAROB2023, Oita, Japan,

Kudo,

Feb. 9-12, 2023

Makoto2023/2/12

国際

2023/2/12

国際

2023/3/1

国際

Sakamoto

decomposition method, oral
56

Parallelfull-wave

Amane Takei, Nanako

ICAROB2023, Oita, Japan,

field

Mizoguchi,

Kento

Feb. 9-12, 2024

on

Ohnaka,MakotodomainSakamoto

electromagnetic
analysisbasedhierarchicaldecomposition method, oral
57

3 Performance Evaluation of
AoS

and

SoA

Incompressible

QI JI, Kenji ONO

for

2023

Conferenceon

Advanced Topics and Auto

Fluid

TuningSimulation on GPUs, oralinHigh-Performance
Scientific Computing

3．受賞等
No.

名称

受賞者氏名

授賞機関(学会名等)

受賞時期

国内／
国際

1

Honorable

MentionAward

2

日 本シミュレーシ ョ ン

Chenyue Jiao, Chongke Bi,

The China Visualization and2022 年 7

Lu Yang, Zhen Wang, Zijun

Visual Analytics Conference

月

Xia, Kenji Ono

2022

武居周、工藤彰洋

日本シミュレーション学会

学会論文賞
3

2022 年 8

国際

国内

月

第 44 回 風力エネルギ

内田 孝紀、小野 謙二、飯田

第 44 回風力エネルギー利用シ

2022

ー利用シンポジウム

明由、吉村 忍、加藤 千幸、

ンポジウム

12 月

ポスター賞

山出 吉伸、今村 博、植田 祐
子- 88 -

年

国内