令和３年度高性能汎用計算機高度利用事業「富岳」成果創出加速プログラム「「富岳」が拓く Society 5.0 時代のスマートデザイン」
成果報告書

令和４年５月３０日
理化学研究所計算科学研究センター
坪倉 誠

目次
１．

補助事業の目的 .......................................................................................................................... - 1 -２．令和３年度（報告年度）の実施内容 ................................................................................................ - 1 ２－１．当該年度（令和３年度）の事業実施計画 .............................................................................. - 1 ２－２．実施内容（成果） .................................................................................................................. - 3 ２－３．活動（研究会の活動等） ..................................................................................................... - 41 ２－４．実施体制 ............................................................................................................................... - 45 --2-

補助事業の名称「富岳」成果創出加速プログラム「富岳」が拓く Society 5.0 時代のスマートデザイン

１．補助事業の目的
少子高齢化に伴う労働生産性の低下、確実に来襲する新興ウイルスによるパンデミックや予期せぬ自然災害への対応、さらにはカーボンニュートラル社会の実現等、激変する社会に対応して今後も我が国がものづくり大国としてその国際競争力を維持していくためには、単なる性能試験の代替としてのシミュ
レーションを脱却し、新たなデジタルエンジニアリング技術を創生することが急務である。ここでは、「富岳」の性能を十二分に活用できるアプリケーションソフトウェアによる高性能大規模シミュレーシ
ョンと、多目的最適化、機械学習、データ同化等のデータ科学を融合させることで、サイバー空間とフィジカル空間を高度に結合させた次世代デジタル設計システムを構築する。即ち、深層学習を用いたビッグデータ解析とは異なり、限定されたデータ群に物理シミュレーション技術を効率的に統合することで、演繹的なシミュレーション技術と帰納的な機械学習技術を高度に融合した最適設計システムを開発する。
開発したシステムは、我が国の機関産業分野を構成する自動車・機械・土木・建築分野に対して以下の四つの課題を選定し、その有用性を実証することで成果創出を狙う。構築したシステムは、我が国の基幹産
業分野を構成する自動車･機械･土木･建築分野の応用課題への展開を目指し、自動車空力設計、都市・建築性能設計、室内環境設計、カーボンフリーガスタービン設計に適用することで、その有用性を実証する。研究開発はそれぞれのターゲットに対応した産学コンソーシアム（既に設置済み）を活用し、産業界での実用化を想定したシステムの設計から実証、さらには産業界メンバーへの技術教育を通した HPC 人
材育成を経て、社会実装までを連携して行う。これにより、我が国を代表する基幹産業分野において、Society 5.0 時代のスマートデザインを実現する。

２．令和３年度（報告年度）の実施内容
２－１．当該年度（令和３年度）の事業実施計画
我が国の基盤産業を構成する機械･材料･土木･建築分野の実用的な課題より、特に経済波及効果の大き
い四つの実証テーマを設定し、大学等の研究者と産業界のエンジニアが強力に連携して、課題を推進する。以下、それぞれのテーマの本年度の実施内容の詳細を示す。
(1) サブ課題 A「意匠空間を考慮した AI 支援多目的最適化による自動車空力デザイン」
性能設計空間と意匠デザイン空間を融合させた新たな空力多目的最適化システムを構築し、実際の自動車空力設計に適用することで、その有用性を実証する。多次元設計空間での最適化を実現するために、遺伝的アルゴリズムに基づく多目的最適化を自動車設計で実用化するフレームワークを構築する。
キャパシティコンピューティングによる多目的最適化で問題となる計算負荷低減技術として、機械学
習によりサロゲート（代理）モデルやリダクション（縮約）モデルを構築し、多目的設計の大幅な時間短縮を目指す。
本年度は、AI 支援多目的最適化フレームワークの設計、基盤アプリの「富岳」チューニング、基礎-1-

流れを対象とした数千万セル規模の縮約モデル構築技術の開発、多目的最適化ベンチマークモデルの構築と解析の実施、及び意匠空間の数理モデル化とターゲット車体の選定を行う。
（協力機関：神戸大学、広島大学）(2) サブ課題 B「変容する都市・建築の自然擾乱対応の性能設計」
データサイエンスに基づきメタボリズム（新陳代謝）的に変容する都市・建築モデルをサイバー空間
上に再現し、都市・建築の離散的空間での換気性能、グリーンインフラに基づく気温調整機能を定量化する。気象モデルと乱流シミュレーションを融合した広範囲の周波数帯の乱流変動を再現し、極端気象（スーパータイフーン、竜巻、豪雨）下での都市・建築の挙動を把握することで防災性能を予測する。
Society 5.0 の理念を踏襲しながら、自然起源の擾乱の都市・建築へのインパクト推定技術の社会貢献をめざし、コンソーシアムを基本として、産業界を巻き込んで実務的展開と不確定性に基づく性能設計の理念と方法を構築する。
本年度は、気象学的なスケールでの自然擾乱の作成、広域都市解析に基づく流入条件データベースの基本データ作成、風・熱・湿度の連成解法のモデル開発、及びメタボリズムに基づく都市・建築の変容に関する調査を行う。
（協力機関：東京工業大学、横浜国立大学、近畿大学、お茶の水女子大学）(3) サブ課題 C「新興ウイルス感染症にロバストで健康・快適・サステイナブルなポストコロナ時代の室
内環境設計」
各種の空気感染性汚染物質の生成から室内空気中での輸送動態、さらに体内動態と免疫システム応
答･生理反応までを包括的に解析し、可視化する統合的システムを開発し、「富岳」に実装する。室内の快適性と感染症に対するロバストネスを両立させる多目的最適化による設計手法を提案する．
本年度は、統合シミュレーション技術の開発（人体幾何形状と呼吸器系幾何形状を精緻に再現した数
値人体モデルの開発）、ベンチマーク問題の設定と解析モデルの作成、統合シミュレーションの実施、
多目的最適化フレームワークの設計、及び基盤アプリの「富岳」チューニングを行う。
（協力機関：九州大学、神戸大学）(4) サブ課題 D「マルチコンポーネント統合シミュレーションによるカーボンフリーガスタービン設計」
ガスタービンを構成する圧縮機、燃焼器、タービンを完全に統合して LES で一気に解く、世界初、真のマルチコンポーネント統合シミュレーション技術の開発を行う。本統合シミュレーションが、フラッシュバックや燃焼振動等の異常燃焼、さらには、窒素酸化物（NOx）、一酸化炭素（CO）、およびすす等の様々なエミッションの正確な予測に有用であることを示す。本統合シミュレーション（サイバー空間）結果から、異常燃焼発生時の物理メカニズムを解明し、さらに異常燃焼の兆候を自動で探知するた
めの A I センシング技術を開発することによって、これらの知見、技術が実機（フィジカル空間）の設計や最適操作条件選定に貢献しうることを実証する。
本年度は、対象とするガスタービンの選定・詳細確認、圧縮機・燃焼器・タービンの各コンポーネン
トの解析実施（試計算）、乱流モデルの調査、選定、及び水素、アンモニア燃焼詳細機構の調査および
1 次元解析による精度検証を行う。

-2-

（協力機関：京都大学、九州大学）(5) プロジェクトの総合的推進
以上のサブ課題 A～D の進捗を管理し、課題間の連携を促進するために、各課題の代表者を集めて統括する推進委員会を設置した。推進委員会は、３か月に一度程度の頻度で理研が主催して委員会を開催する。また、年に一度の成果報告会や、他の「産業競争力の強化」課題と連携した統合ワークシ
ョップを企画・開催する。
２－２．実施内容（成果）(1) サブ課題 A「意匠空間を考慮した AI 支援多目的最適化による自動車空力デザイン」(a) 研究成果の要旨
本年度は、AI 支援多目的最適化フレームワークの設計、基盤アプリの「富岳」チューニング、基礎流れを対象とした数千万セル規模の縮約モデル構築技術の開発、多目的最適化ベンチマークモデ
ルの構築と解析の実施、及び意匠空間の数理モデル化とターゲット車体の選定を行った。この結果、基盤アプリ「CUBE」の「富岳」へのチューニングでは、単純流れ解析における単体性能で理論性能比 6.4%、実用解析における並列性能では 51,200 ノード並列の Weak Scaling 性能で並列化効率 74%を達成した。縮約モデル構築技術では、流れ場モード分解型の畳み込みニューラルネットワーク（MD-CNN-AC）でより多数のモードに基づく CNN を高速に構築するため、富岳上での並列化、分散学習化を実現した。また、100 台超の市販 SUV 車両群の 3D モデルを対象に形状特性の統計デー
タ解析を行い、その特徴量に基づく設計変数を定義することができた。さらにデザイナーの意図と相関のある形状の特徴量を抽出することで、形状特性に基づく意匠性空間の数理モデル化の目途が立った。
In this fiscal year, we designed an AI-assisted multi-objective optimization framework,
tuned the fundamental application on the supercomputer "Fugaku", developed a technology
for constructing a reduced model with tens of millions of cells for fundamental flow,
constructed and analyzed a multi-objective optimization benchmark model, and mathematicallymodeled the design space and selected target vehicle bodies. As a result, the tuning of
the basic application "CUBE" on "Fugaku" achieved 6.4% of theoretical performance insingle-core performance in simple flow analysis and 74% of parallelization efficiency inpractical analysis with 51,200 nodes parallel weak-scaling performance. Regarding theconstruction technique of reduced order model, parallelization and distributed learning
were achieved on “Fugaku” for faster construction of CNNs in a mode decomposing
convolutional neural network autoencoder (MD-CNN-AC) based on a larger number of modes offlow field. Moreover, statistical data analysis of shape characteristics was performed on
3D models of over 100 commercial SUVs (Sport Utility Vehicles), and design variables basedon their features were defined. Furthermore, by extracting the shape features that
correlate with the designer's intention, the prospect of mathematical modeling of the-3-

design space based on shape characteristics was established.
(b) 本研究の内容
本年度(2021 年度)の研究成果を項目別に示す．
(b-1) AI 支援多目的最適化フレームワークの設計
既往の成果をベースに、意匠空間を考慮した自動車空力形状の多目的最適化問題への適用を想定
して、Fig. 2.2.1-1 に示す AI 支援多目的最適化フレームワークを設計した。車体形状を表現する複数のパラメータを設計変数とし（パラメトリック最適化）、最適化手法にはチェビシェフ距離と選好領域と呼ばれる考え方を用いた多目的進化アルゴリズム CHEETAH/R (CHEbyshev-Epsilon optimizer
AlgoritHm)を使用する。設計変数に応じた車体形状は、基準車体形状を市販 CAE ソフトウェア（DEP
MeshWorks）の形状モーフィング機能によって変形させることで生成し、形状変更に伴う解析格子の再生成作業を必要としない複雑現象統一解析コード CUBE を用いた数値解析によって目的関数を評価する。そして、設計変数とする車体形状パラメータを、デザイン的な意図と相関のある形状特性とすることで、意匠特性を考慮した形状最適化の実現を目指す。また、機械学習技術による最適解探索の効率化として、進化計算の過程で実施した数値解析結果を基にガウス過程回帰モデルによるサロゲートモデルを構築し、目的関数の評価に要する計算コストの低減を図る。
これらの要素技術に基づくフレームワークの原型は、既に Feasibility Study を終えており、本
年度は、サロゲートモデル構築時のカーネル関数とハイパーパラメータのフィッティングに使用する最適化アルゴリズム、および目的関数のハイブリッド評価において数値解析の適用を判定する指標について検討した。また、設計変数となる形状パラメータ定義の検討として、項目(b-4)で後述する市販車両の形状特性の統計分析も行った。

START

Multi-objective evolutionary algorithm“CHEETAH/R”

Parametric geometry morphing(commercial software “DEP MeshWorks”)

Generate parameter sets
for new generation
Parameters(text)

Parameters(text)

Shift to the hybrid evaluation of
objective function after some generations.

Morphing geometry
based on the parameters
Geometry (STL file)

Aerodynamics simulation byComplex Unified Simulation Framework“CUBE”

Evaluate objectivefunction by asurrogate modelYes

Is CFD simulationnecessary?

No

Evaluate objective function
by CFD simulation
Objective functions(text)

No

Reached tothe specified generation?

Yes
END

Fig. 2.2.1-1 Framework of AI supported multi objective optimization.
(b-2) 基盤アプリの「富岳」チューニング
統一解法フレームワーク CUBE について、プロダクトランに近いセッティングを使用し、富岳上における性能測定を実施した（基盤アプリが共通であるため、本課題は、成果創出加速プログラム：「「富岳」を利用した革新的流体性能予測技術の研究開発」と連携して実施）。

-4-

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

Incompressible

Domain size

18.75m×8.75m×3.75m

Number of cubes

32,341

Number of cells per cube

16×16 ×16 = 4,096

Total number of cells

32,341×4,096 = 132,468,736

Minimum cell size

1.2207 mm

Time step size

2×10-6 sec

Integration time

2×10-4 sec (100 steps)

Reynolds number

2.572×106

Time integration scheme

Crank-Nicolson

Poisson solver

Red-Black SOR

Table 2.2.1-2 Scaling of time-stepping loop.
Number of

Number of

nodes

cubes

Number of cells

Elapsed time

Scaling

100

32,341

132,468,736

22.92 s100.00 %800

260,331

1,066,315,776

32.63 s70.23 %

6,400

2,082,648

8,530,526,208

33.46 s68.49 %

51,20016,661,184

68,244,209,664

30.92 s74.11 %

-5-

Fig. 2.2.1-2 Elapsed time against number of nodes.

Subroutine name

100%
90%
80%
70%

vscf_ce2nd
field_intf_unpack_edge_e
cns_absorbing_bc

tstep_lusgs_slts_mix

44.2 %

flux_roe_lm

17.6 %

field_intf_unpack_efield_intf_unpack_e

12.6 %

cns_absorbing_bc

5.7 %field_intf_unpack_edge_e

5.1 %vscf_ce2nd

4.8 %field_intf_q_e

2.0 %

field_intf_unpack_corner_e

1.7 %field_intf_pack_e

1.4 %

60%
50%

flux_roe_lm

40%
30%
20%
10%
0%

Cost ratio

tstep_lusgs_slts_mix

Fig. 2.2.1-3 Cost distribution in time-stepping loop.

-6-Table 2.2.1-3 Simulation parameter setups.
Solver type

Compressible

Domain size

20m × 10m × 5m

Number of cubes

32,341

Number of cells per cube

16 × 16 × 16 = 4,096

Total number of cells

32,341 × 4,096 = 132,468,736

Minimum cell size

1.2207 mm

Time step size

1.0 × 10-4 sec

Integration time

1.3 × 10-3 sec (13 steps)

Reynolds number

1.1 × 106Time integration scheme

Lower-Upper Symmetric-Gauss-Seidel(LUSGS) implicit method

Viscous term

2nd order central difference

Table 2.2.1-4 Execution performance of kernel loop.
Exec. time

Calc. performanceMem. bandwidth

Original

5.08 s19.7 GFLOPS (2.57%)

5.92 GB/s (2.31%)

Exchange array

4.50 s22.23 GFLOPS (2.90%)

6.24 GB/s (2.44%)

4.05 s

11.63 GFLOPS (1.52%)

6.56 GB/s (2.56%)

Auto loop fission

3.84 s12.29 GFLOPS (1.60%)

7.19 GB/s (2.81%)

Extract if condition

3.49 s

14.97 GFLOPS (1.95%)

8.01 GB/s (3.13%)

dimension
Loop collapse &manual fission

Fig. 2.2.1-4 Tuning results of kernel loop.

-7-

(b-3) 基礎流れを対象とした数千万セル規模の縮約モデル構築技術の開発
非定常三次元流れ場に対して機械学習によるモード分解を行い，主要なエネルギーモードを抽出することで流れ場の特徴的な流れ構造を抽出し，その主要モードから流れ場を再現するための縮約モデルの構築に取り組んでいる。今年度は、富岳で実装済みの大規模分散学習の流れ場再現精度の向上を目的とした機械学習手法の高度化を実施した。我々が使用した、3 次元流れ場を対象としたモー
ド分割用のネットワークの構造を Fig. 2.2.1-5 に示す。今回、このネットワークをレイノルズ数 1,000の 3 次元円柱周り流れ場に適用し、縮約モデルによる流れ場の再現精度を評価する。
高精度シミュレーションの流れ場及び機械学習によって再構成された流れ場のスナップショット
を Fig. 2.2.1-6 に示す．機械学習によって再構成された流れ場は，高精度シミュレーションによって表現された円柱後流の複雑な流れ構造が再現できていないが，学習の反復回数(エポック数)を増やすことにより，より複雑な構造が見えてくることが分かる．

3D flow fieldsnapshots

3D flow fieldsnapshots

Latent vector

a1

Input

Decodera2

Output+Encoder
DecoderFig. 2.2.1-5 Network structure for 3D flow

Fig. 2.2.1-6 Comparison of flow field

field mode decomposition.

snapshot.

Fig. 2.2.1-7 は，高精度モデルの流れ場及び機械学習によって再構成された流れ場を対象に POD による分割を行い，各モードについてエネルギーの大きい順にソートしたものである．図の縦軸は各モードの持つ(規格化された)エネルギーを表している．機械学習による再構成は，全て 7 つのモードを使用しているが，学習の反復回数を 5000, 8000, 12000 としたものをプロットしている．高精度モデルに比べると，機械学習によって再構成された流れ場は十数モード程度でエネルギーが低下し始めるが，反復回数を増やすほど，多数のモードに至るまで相対的に大きなエネルギーを保持できていることがわかる．
最後に，本研究で使用した分散学習スクリプトの並列計算性能を示す．Fig. 2.2.1-8 は学習全体及びそれに含まれる各処理のスケーリングを示したものである．また，Table 2.2.1-5 は学習全体の演算性能及びスケーリングの数値を示している．最大規模となる 215 のモード数(20,736 ノード)を使用した学習全体の演算性能は 6.43PFLOPS であり，また 7 つのモードを使用した学習を基準としたそのスケーリングは 88.7%と良好な値を示している．

-8-

Fig. 2.2.1-7 Energy distribution

Fig. 2.2.1-8 Energy distribution of each POD mode.

of each POD mode.
Table 2.2.1-5 Weak scaling of entire training loop.
Modes

Nodes

Cores

Elapsed time

PetaFLOPS

Scaling

7768

36,864

196.6 s0.26

100.0%

26

2,592

124,416200.2 s

0.85

93.9%

63

6,144

294,912204.7 s

1.95

90.7%

215

20,736995,328208.2 s

6.43

88.7%

(b-4) 意匠空間の数理モデル化に向けた市販車形状の分析
デザイン的意図の特徴軸により構成される意匠空間を考えるため、その特徴軸に対応する車体形状特性の定義、定量化を検討した。ここでは、市場で販売されている車体形状が様々なデザイン的意
図を含むと仮定し、意匠性に影響を与えうるとされる形状特性の主成分分析を行うことで、代表的なデザイン的意図の方向性に対応した車体形状特性の抽出を試みた。主成分分析による車体形状特性の統計データ解析の対象とした、124 台の市販 SUV 車両群の 3 次元 CG(Computer Graphics)モデル（Hum3D.com 社製）の側面視画像を Fig. 2.2.1-9 に示す。また、意匠性に影響を与えうる車体外観形
状特性としては、車両空力への影響も大きいと考えられるマクロな形状特性を考慮することとし、各部位の寸法比率に基づくプロポーション要素に着目した。先行研究（原田・森，1993；南本・原田，
2018）を参考に、SUV 車両の側面視形状におけるプロポーション要素を Fig. 2.2.1-10 に示す 22 項目
設定し、各車両の側面視画像からその値（比率）を算出して、その主成分を求めた。

Fig. 2.2.1-9 Sideview of the investigated SUVs.

-9-

①ルーフ長さ/全長⑫全高/全長②キャビン長さ/全長⑬ノーズ高/全高③フード長さ/全長⑭ノーズスラント量/フード長さ⑮基調ライン基点高さ/全高⑯ベルトライン高/全高④L2/全長⑤ホイールベース/全長⑥基調ライン基点長さ⑰FR バンパー下端/全高/ホイールベース⑦前輪-L の長さ
⑱RR バンパー下端/全高/ホイールベース⑧FR オーバーハング/全長 ⑲サイドシル下端/全高⑨RR オーバーハング/全長 ⑳基調ライン角⑩RR オーバーハング/FR オーバーハング⑪H1/L1㉑RR ウィンドウ傾斜角㉒FR ウィンドウ傾斜角

Fig. 2.2.1-10 Proportional elements of an SUV in side view.
その結果得られた、
64 台の SUV 車両の車体形状から抽出された側面視プロポーションの主成分と、その主成分が表す形状差の傾向を Fig. 2.2.1-11 に示す。第 1 主成分は約 4 割の寄与度を持ち、Fig.
2.2.1-11(a)中に示す多くのプロポーション要素を含むものであった。第 1 主成分の主成分得点順に代
表的な車体を並べた Fig. 2.2.1-11(b)では、主成分得点が小さいほど車体上面がより滑らかで流線形のような外形を持つ傾向が示されている。また、Fig. 2.2.1-11(a)に示すように、第 2 主成分として側面視形状における重心位置の後方への偏り、第 3 主成分として床下空間の大きさに関連するプロポーション要素の組み合わせが、それぞれ抽出された。

(a) principal components of the sideview proportion(b) geometrical trend explained by the first principal component score
Fig. 2.2.1-11 Principal components of the sideview proportion of commercial SUVs.

- 10 -

さらに主成分として抽出された形状の特徴量とデザイン的な意図との関係を調査するため、SD 法によるアンケート調査を設計した。その実現可能性調査では、人が車体画像から受ける印象（デザイン的な意図）が、形状の特徴量（主成分得点）と有意な相関を持ちうることが確認された。以上の成果により、形状特性に基づく意匠性空間の数理モデル化の目途が立った。
(b-5) 多目的最適化ベンチマークモデルの構築とターゲット車体の選定
形状特性分析を行った市販 SUV 車両群の 3D モデルを対象として、車両周囲流れの時系列数値シミ
ュレーションを行い、同モデル群が構成する空力特性空間の構築を行った。流体解析には CUBE を使用し、実車スケールにおいて最小格子幅 7mm、総格子数 2 億 7 千万セルの解析格子を用いて、実時間
2 秒間に相当する陰的 LES 解析を実施した。Fig. 2.2.1-12 に解析結果の平均抗力係数で規格化した抗力・揚力係数の散布図と、代表的なプロットに対応する車両中央断面の流速分布を示す。図より、特に高抵抗な車体ではキャビン形状や車高に特徴がみられる。なお、使用した 3D モデルデータは CG 用途のため、車体下部を中心に車両空力性能に影響する付加部品やタイヤ周りの構造が省略されており、実際の車両が持つ空力性能とは異なっている可能性が高い。しかし、それらの形状はいずれの車両データについても省略や簡略化をされていることから、デザイン性と関係性の深い車体上部形状と空力性能との関係を調査するうえでは有効であると考える。

Fig. 2.2.1-12 Scatter plots of normalized drag and lift coefficients,and streamwise velocity distribution in the center of the vehicle.
さらに、前項と本項で示した車体形状のプロポーション特性と空力係数特性の両情報から、平均的な車体形状と抗力・揚力係数の双方で好適値を示す車体を、今後の研究において多目的最適化の対象とするターゲット車体として選定した。Fig. 2.2.1-13 に車体の側面視プロポーションの第 1・第 2 主成分の主成分得点と抗力・揚力係数の散布図を示す。図中の赤色のプロットがターゲット車体に相当する。

- 11 -

(a) Principal component scores(b) Aerodynamic coefficients

Fig. 2.2.1-13 Scatter plots of aerodynamic coefficients and principal component scores ofsideview proportion of SUVs (Red plot means the target vehicle).

(2) サブ課題 B「変容する都市・建築の自然擾乱対応の性能設計」(a) 研究成果の要旨
本年度は、気象学的なスケールでの自然擾乱の作成、広域都市解析に基づく流入条件データベースの基本データ作成、風・熱・湿度の連成解法のモデル開発、及びメタボリズムに基づく都市・建築の変容に関する調査を行った。具体的には、極端気象を対象として、気象モデルに基づく流入変動風の作成し、台風時・突風時の都市・建築の機能性・安全性を検証した。各種統一解法を建築物に適用し、洪水時の耐水性の評価、過度な風荷重時の弾塑性挙動を分析した。さらに非構造格子 LES を用いて猛暑下での複雑な市街地の暑熱リスクの評価を進め、複雑な形態での熱環境を明らかにした。また、未来の変容した都市を対象とする強風防災の可能性を確認した。
In this fiscal year, we created data of meteorological-scale disturbances and prepared
fundamental database of inflow conditions based on broad-area urban simulation, developed
coupling models for wind, heat and humidity, and researched transformation of urban areaand building based on metabolism theory. Specifically, inflow turbulence based on
meteorological models were created for extreme weather events, and the functionality andsafety of urban areas and buildings were examined in the events of typhoons and strong
wind gusts. Various unified solution methods were applied to buildings to evaluate theirflood resistance and to analyze their elastoplastic behavior under excessive wind loads.
Using unstructured grid LES, the thermal risk under extreme heat were evaluated in thethermal environment of the urban area with complicated shapes. Furthermore, the study also
investigated the potential for high-wind disaster prevention for future transformed urban
area.
(b) 本研究の内容(b)-1

気象擾乱を再現した流入風の作成・建築物の変動風圧力推定- 12 -

(b)-1-1 気象庁非静力学モデルによる 2019 年台風 19 号の再現解析
本研究では東京都心で 41.5m/s の最大瞬間風速を記録した 2019 年台風 19 号を対象に、実都市における強風の実態分析、および建築物の精緻な風圧力による強風リスク評価のための LES を「気象モデル・工学 LES ハイブリッド解析」[2-1]の枠組みにより実施した。なお本項および次項の気象モ
デルによる解析は、「富岳」成果創出加速プログラム「防災・減災に資する新時代の大アンサンブル
気象・大気環境予測」と協力し、研究を実施している。気象場の台風構造の再現は、再現気象庁非静力学モデル(NHM)により 1km、 250m、 50m の 3 段階の one-way ネスティングを用いて計算を行った。
最外解析領域 Domain 1 は最接近の 6 時間前まで近づけることで上陸位置の再現性や上陸後の移動速度の低下に改善が見られた。ただし最接近時の経路は実際の経路より 40-50km 程度北西側となっている。最内領域 Domain 3 では解析におけるピーク時刻前後 60 分間の計算を行った。上流地形の影
響や乱流発達を再現できるよう南北の格子数を増やし、房総半島の山地を含む領域に範囲変更した
結果、解析時間のほとんどの時間帯で都心上空に乱流が維持され，幅数 km 程度の大きな主流方向に伸びた筋状の風速変動の中に、更に細かな数百 m スケールの風向や風速の揺動が生じている様子が確認された。東京管区気象台の位置での時刻歴から台風構造に基づく 10 秒から 20 秒程度の風速変動の再現が確認された。解析での平均風速は、決定論的解析での経路・中心気圧の再現性の限界から、実観測に比べ 6m/s 程度高くなった。
(b)-1-2

NHM アンサンブルによる強風の極値の評価と極大ケースの解析

気象モデルによる実現象の解析は、数値解析や初期・境界条件の精度、同化に利用した観測結果の誤差等に起因する不確実性が含まれている。本研究では事象の確率分布を表現しうる 1000 メンバー
規模の気象モデルの大アンサンブルの結果を用いて、気象擾乱の不確実性を考慮した市街地の BCMLES 解析を実施することを目標としている。本年度は準備段階として、台風の気象モデル大アンサン
ブル結果の風速に基づく整理、アンサンブルから抽出した 5 ケースの 1km 解像度までのダウンスケ
ーリングの試計算を実施した。
気象モデル大アンサンブルのデータは、2019 年台風 19 号について、アンサンブル変分法（EnVar）
により得られた 5km 解像の NHM 解析 1000 メンバーであり、気象研究所より提供を受けた。アンサン
ブル結果は、東京管区気象台の位置における最大風速に基づいて整理された。地上 10m での風速の最頻値は 10.5m/s から 11m/s であり、観測地点の高さや地表面粗度に基づく換算から実測値(地上35.3m、16.8m/s)と妥当に対応していることが判断された。また台風経路と対応する最大風速の軌跡
は御前崎付近から東京湾の間の広範囲で摂動しており、風速値も比較的分布を持つことが確認された。
次に、抽出した少数のメンバーの特性をダウンスケーリングにより確認した。1000 メンバーのうち東京接近時に最も高風速となった 3 ケースを抽出し計算した結果、いずれのケースも実況の西側を通過し、経路降水や気圧の分布がよく類似する傾向を示した。最接近時の中心気圧が 955hPa に達し関東を通過する台風としては最大クラスのものが再現されたことが確認された。一方、最も低風速の 2 ケースを抽出し計算した結果、経路はより東側を通り、進行速度が著しく衰退した状態で東京へ最接近する傾向を示した。以上の結果から、様々な風速・経路での気象場の特性を抽出し分析する
ことの可能性が確認された。また、今回の解析では各メンバーの経路や気圧・降水分布は、初期場の- 13 -

不確実性の影響も認められたものの、境界値に用いた全球アンサンブル解析（特異ベクトル法(SV)・
50 メンバー）の影響を大きく受けている可能性が確認された。したがって、50 メンバーの境界条件からガウス分布に従う 1000 メンバー分の境界条件に再構築した上でアンサンブル解析を再度実施
し、その結果に基づく気象場構造の解析や BCM-LES による地表近傍の強風性状に関する統計的分析を進めることで結果を改善しうることを確認した。
(b)-1-3

気象モデル解析に対する高周波変動成分の再生成計算

気象モデル解析された流れ場を市街地 LES の流入条件として利用するためには気象モデルで減衰しやすい高周波変動成分を再生成させることが重要である。乱流エネルギーカスケードによる高周波成分の再生成手法[2-1]を節(b)-1-1 の気象モデルの解に適用した。東京管区気象台のある大手町と高層密集街区の渋谷を含む水平 12 km×7 km、鉛直約 2 km の領域において格子解像度 25m で計算を行い、流入風を 12 分間分取得した。変動付加前の気象場に基づいて、さらに高波数の 3 次元的な乱れの構造が再生成されたことが確認された。
(b)-1-4

都市広域解析による BCM-LES

市街地で発生するピーク風速や風圧を解析するための気象擾乱を流入条件とする BCM-LES は、理
研 R-CCS により開発された CUBE[2-3]により広域解析(3 km × 6 km × 1.5km、解像度約 1.5 m)と詳細解析（950 m × 950 m × 950、約 93 cm）の 2 段階で気象擾乱を東・南・西の 3 面から与えて、実施した。CUBE の格子系は階層型の直交マルチブロック格子であり、個々の建物・地形を再現する
地表の複雑境界には immersed boundary(IB)法を用いている。
大手町での解析から、今回の台風（2019 年台風 19 号）では、平均速度の垂直方向分布において、通常の境界層と同様のべき指数型に近い形であったが、主流直交鉛直断面の可視化から直径 1 km 程度までの大きさの異なるロール状の縦渦が存在し局所的に運動量を地表付近に輸送するプロセスが確認された。スペクトルからは、-5/3 乗則に従った変動が高度 300 m で 0.5 Hz 程度、高度 40 m では 1Hz 以上の周波数まで再現される一方、低周波数側には気象擾乱によると推察されるエネルギーがあり、広帯域の乱流流入風が得られた。LES 解析と実測値の東京管区気象台での 10 秒間隔の風向・風速の時刻歴の比較により、時間的な変動パターンについて定性的な類似性が確認された。市街地内の異なる高さの平均風速と最大瞬間風速の分布特性は、高さ 2018 年台風 21 号時の大阪での結果[24]と符合した。地上 85m 付近での最大瞬間風速の分布からは細い筋状のピーク領域が確認され、上空から局所的に降下した高速流体塊の移流により発生している可能性が推察された。都市キャノピー内部の地上 35m 付近では、建物により作られた乱れの影響がより顕著となった。地上 10m 付近では最大瞬間風速 40m/s を超える特に強い風の発生は高層建築物の周辺に集中し、剝離せん断や吹きおろし、縮流による効果が顕著であることを示した。また高層建物周辺で風が通りに沿って流下する場所では、離れた場所でも瞬発的に強風が発生する様子が確認された。一方、本台風においては結果の定量性については課題が見られた。台風経路・勢力の再現性により平均風速は気象モデルのみの結果と同様 5m/s 程度、実況より高くなった。また、変動は風向、風速ともに観測に比べて小さく、最大瞬間風速は 33m/s であり実測でのピーク値 40.5 m/s は十分に再現できていない。実測のピーク風速の発生状況を確認すると 1 分間程度、最低・最高風速が上昇し強風が継続したことが確認される- 14 -

ことから、気象場のマイクロスケールの構造の影響（台風に伴うロール構造やダウンバースト的構
造、壁雲内の鉛直渦[2-5]）あるいは流れの複雑な市街地での低頻度な強風パターンの発生などが要因として考えられる。今後アンサンブル解析を活用し気象モデルや市街地解析の不確実性に関する検討を進める。
第二段階の解析では対象建物に作用する風圧力を解析した。上流側の隣接建物の影響を受ける当該建物における風圧分布が適切に再現された。また最大負圧を経験した高層棟上層階前縁の点 P の圧力係数の時刻歴からは、解像度 93cm の解析では-2.5 を超えるピークが発生した。解像度 1.5 m においても同様の変動傾向を持つ圧力波形が得られたが振幅は大幅に小さかった。このことから、建物
幅の 1/50 から 1/100 の範囲の格子解像度の変化は、直交格子・IB 法による風圧解析に重大な影響を示すことを確認した。
また負圧のピーク前後の非定常場からは、前後の数分間、上空の気象場の強風の高度が建物高さま
で下がり比較的強い負圧が継続して作用したこと、前方の隣接建物を越える風によって、対象建物の高層棟上層階において逆円錐渦が発生し、鉛直方向の風速シアにより後方に引伸ばされて強化されて負圧のピークが発生したことが推察された。
(b)-2 広域市街地解析に基づく流入条件データベース作成と実建物解析への適用(b)-2-1

広域解析に基づく流入変動風データベースの取得とのその乱流特性

建築物の風圧力の予測においては、滑面もしくは粗面を仮定した乱流境界層が使用され、粗面においては粗度ブロック状で生成された乱流を流入条件として用いる。しかしながら、東京の市街地では
高層建築物群などを含んだ様々な粗度条件を有しており、流入条件も都市境界層内の様々な粗度条
件の影響を受けると考えられる。さらに、台風のような強風時には大気境界層内には大きなスケールで変動する気象擾乱が存在し、建物の風圧力特性に対して影響を及ぼすと考えられる。
本サブ課題では、東京都心部に暴風をもたらした台風を対象として、WRF-LES によって得られた気象場に基づく流入条件を用いた CUBE での広域解析（乱流境界層の流入条件、2017 年台風 21 号の気象場に基づく流入条件）を実施し、将来変容する都市のレジリエンス強化のための広域の強風のデー
タベースを構築した。
広域解析では GIS データと国土地理院の標高データにより、東京中心部 9km×9km の解析領域内にある建物形状と地形の起伏を再現している。最小空間解像度は 2.25m、要素数は約 9.6 億である。階層型直交格子に基づく解析を実施しており、使用コードは理化学研究所計算科学研究センター(RCCS)が開発した階層型直交格子のコード CUBE である。
流入境界条件には乱流境界層に基づく流入条件と WRF-LES に基づいて作成された気象擾乱を含む流入条件を用い、2 ケースの広域解析を実施した。広域解析では、風圧力を予測する実建物の風上側の位置に PlaneA、B にサンプリング面を設け、風速の時空間データを保存した。
本研究では保存したデータに基づき、WRF-LES によって得られた実都市の流入変動風のデータベースにおける風速の鉛直プロファイルを確認した結果、本広域解析ではコリオリ力は考慮されておらず v 方向の流速（𝑣/𝑢 ）で 1～1.5 程度のずれが生じているものの、流入面の鉛直プロファイルは流入面から約 2.5 ㎞、5km 下流にある PlaneA、B でも維持されていることを確認した。さらに、PlaneBにおける u、v、w 成分におけるスペクトル密度分布を確認した結果、乱流境界層に基づく流入変動風- 15 -

のケース、気象擾乱を含む流入条件のケースともに、0.1Hz 程度まで 5/3 乗則の変動が維持されており、実スケールで数秒程度の変動が得られていることが確認できる。
(b)-2-2

流入変動風を用いた CUBE 解析に基づく実建物の耐風性能の予測(b)-2-1 で作成した流入変動風のデータベースは詳細解析領域での流入条件として接続され、実都市の粗度条件と気象擾乱に基づく変動の影響を考慮した風圧力の予測を行う。PlaneB の流入条件を接続した CUBE での実建物詳細解析（空間解像度 0.25m、格子数約 5 億）であり、複雑な対象建物ファサードの周りの流れ場を再現している。対象建物では壁面から離れた梁や柱上の部材を有する複雑なファサード形状をしており、空間解像度 0.25m の解析では柱の間を抜ける微細な流れ場が再現されており、建物の隅角部からはファサード形状の細かいスケールに対応した微細な乱流構造が確認された。以上のような、流入変動風のデータベースを活用した実建物の解析により、実気象場に近い条件での解析と実際の都市形状を再現した耐風性能の予測が実現されることが期待される。
(b)-3 極端気象被災推定のための各種統一解法の適用(b)-3-1 建物における過度な風荷重時の弾塑性挙動の評価
台風時・突風時の機能性・安全性を確保するためには、建築物に作用する風圧力を予想するだけではなく、外装材に作用する応力や挙動を把握する必要がある。本サブ課題では実建物における過度な風荷重時の弾塑性挙動の評価を実施した。対象とする外装カーテンウォール（CW）は渋谷地区の超高層建物に取り付けられているものと想定し、まずは渋谷地区を対象とした非構造格子に基づく実市
街地の CFD 解析を実施した。解析領域は、2.8 ㎞×1.8 ㎞×1.2 ㎞で風向は南南西、総格子数は約 1億である。本解析により超高層建物壁面で負圧のピークが確認された高層部と低層部の出隅部分に位置する外装 CW を対象に FEM 解析を実施して損傷度を評価した。解析に用いた CW は建築工事標準
仕様書・同解説 JASS カーテンウォール工事に提示されているものを対象に、動的応力解析用の FEMモデルを作成した。FEM 解析の境界条件は、上下を周期境界条件、左右を対称境界条件とし、CW に作用する風圧力は、再現期間 500 年相当の風とした。FEM 解析で得られた CW の反力の時刻歴波形と風圧力の時刻歴波形とを比べると、形状が変化していることから動的な影響を考慮する必要性が確認された。次に、CW の接合部品の疲労損傷の検討を行うために、接続部品を六面体でメッシュ分割し、FEM 解析で得られた CW の反力を接合部品の入力値として与えて、応力解析を実施して接合部品の応力状態を確認した。相当応力最大の所で塑性化、主応力最大部分では塑性化しないとの結果が得られた。
(b)-3-2 CUBE による気液二相流解析に基づく建物浸水評価
既往の洪水解析では、1 次元または 2 次元の洪水解析モデルが使用されることが一般的である。ただし、交差点・下水道・公園・建築物の周囲やその内部等の水の 3 次元的な流動および構造に作用す
る非定常流体力を評価するには、ナビエ・ストークス方程式に基づく 3 次元解析が必要であるが、実際の構造物の複雑な形状を再現するためには膨大な計算コストを要する。また、複雑形状に対する非構造格子の生成には形状データの補正等が必要となることからモデルの作成にも膨大な時間を要する。そのため、詳細な構造物形状を再現した 3 次元洪水解析の研究例は極めて少ない。そこで、界面- 16 -

捕捉法として VOF 法を用いたビルディング・キューブ法に基づくセル中心有限体積法を開発し、3 次
元洪水解析により建築物内部への流入挙動を推定した。
本年度は解析領域境界から流速 3m/s で流入する高さ 2m のダムブレイク流を仮定する。建築物として、階段と地下空間を有する構造を仮定し、窓からの流入を想定する。なお、本研究では建築物を埋め込み境界法によりモデル化する。解析領域は直交メッシュにより階層的に分割し、最小セルサイ
ズは 78.125mm、総セル数は 14,876,672 である。数値計算では、908 プロセス×4 スレッドを用いた。
結果として、ダムブレイク流が建築物の周囲に流れ込み、構造と水流が衝突することで、複雑な水の流動が再現された。特に、建築物の壁面に水流が衝突し跳ね上がり、窓から水が建物内に侵入する挙動が確認できる。また、非定常な水の流動によって、窓枠上部と下部にそれぞれ正圧と負圧が生じていることが確認できた。
以上のような水の 3 次元的な流動および建物壁に生じる非定常な流体力は、
1 次元や 2 次元の洪水解析モデルでは評価が困難であり、3 次元解析によってのみ評価することができるものである。
今後は、本手法を渋谷駅地下空間の浸水評価に適用する予定である．氾濫水の流入を仮定した試行
解析として、総セル数は 651,587,584、最小セルサイズは 48.828125 mm、氾濫水を 3m/s で流入させ、富岳の 3977 ノードを使用して計算した。今後、さらに計算を進め、渋谷駅地下空間における氾濫水の流動形態を解析する予定である(b)-4 WRF-LES と熱放射・LES の連成解析に基づく暑熱リスクの評価方法の構築(b)-4-1 WRF-LES による 2018 年猛暑の気象解析
2018 年 7 月 23 日に関東・東海を中心に記録的猛暑が観測され、埼玉県熊谷市で日本観測史上最高となる 41.1℃を観測した。また、東京都青梅市でも都内観測史上初の 40℃超えとなる 40.8℃を観測し、東京管区気象台では観測史上第 3 位となる 39.0℃を観測した。熊谷の高温発生に関する既往の研究では、Nishi, Kusaka [2-6]は北陸地方からの北西風によるフェーン現象がその一因として考えられると述べているが、東京都市域での高温発生要因やその発生過程の詳細については明らかにされていない。
同日の東京都市域では、熊谷での高温発生要因の一つであるフェーン現象による移流特性によって温度が上昇し、その後、海風が南東側から侵入し温度が低下したと考えられている。ここでは、この記録的猛暑日における東京都の都市域での高温発生状況や温熱環境の変化を WRF および WRF-LESで再現し、東京、青梅、熊谷の気象台観測値との比較を行う[2-7]。
領域気象モデル WRF Version 3.7.1 を用いて、WRF および WRF-LES の解析を実施した。解析対象は関東地方で高温が観測された 2018 年 7 月 23 日 13 時 (JST) 頃とし、解析期間は 2018 年 7 月 22 日
21 時から 24 日 0 時までの 27 時間とした。本年度は熊谷、青梅、東京の各観測データ（10 分値）と、WRF（Domain1）の気温、風速、風向を比較し、さらに東京管区気象台における WRF-LES（Domain4）の結果と観測データ（10 分値）の風向風速、気温の時系列変化を比較した。11 時から 13 時 20 分頃までの気温変化は、WRF-LES と観測値で同様の結果を示した。風向変化を併せて確認すると、観測では13 時過ぎから南風に切り替わり、冷涼な海風の侵入により気温が低下したものと推察される。WRFLES では 13 時 20 分頃まで北西～北北西の風向きとなっており、その後 13 時 30 分～14 時にかけて風向変動を伴いながらも 14 時過ぎには南～南東風に収束している。WRF-LES の結果は、北側の風が- 17 -

吹いている時間帯での気温上昇、および南～南東風（海風）の侵入による気温低下の過程について、観測値を概ね再現できているものと考えられる。
東京における風向の鉛直分布に関する時間変化を確認した。11 時、12 時には上空 2000m まで北北西となっているが、13 時には高さ 100～800m において北北西～北西の範囲で高さごとに風向変化がみられ、14 時には地上から高さ約 400m までが北北東かつ上空では北北西～北を示し、15 時には高さ 1200m 付近まで南寄りの風（南東～南風）となっている。このことから、14 時から 15 時の間に低い位置から海風の侵入が始まり、15 時頃には海風が東京に到達して南寄りの風となっていることがわかる。
高さ 100m 地点での風速、温位分布を確認すると、南風（冷涼な海風）の侵入に伴う気温変化について、風速・温位の空間的な分布とともに、東京管区気象台における変化の過程に対応することが確認できた。
(b)-4-2

熱放射・LES の連成解析に基づく暑熱リスクの評価

気象解析結果を用いて，猛暑日における東京エリアの非定常温熱場を LES で再現する。気象擾乱を考慮した流入変動風は，河合・田村のﾌｨﾙﾀﾘﾝｸﾞ/ﾘｽｹｰﾘﾝｸﾞ手法を介して WRF-LES から生成・接続することで作成した。作成した流入変動風のデータについては，風速，温位の時空間 DB を境界条件として接続され，非構造格子に基づく熱放射・LES の連成解析により，キャノピー内の風速分布や表面温度分布が明らかにされている。このような連成解析方法により，猛暑時の気象場の擾乱を再現しな
がら，市街地における風速・温度のピーク特性を明らかにし，暑熱リスクを評価することが可能となった。
(b)-5 渋谷を対象とした都市変容に関する換気性能と耐風性能の評価
気候変動に伴って激甚化・頻発化する自然災害や新型ウイルスの世界的パンデミック等、現在直面している様々な社会的リスクに対し、これからの都市・建築の在り方としてより高度な防災性能と環境性能が求められている。これまでのアーバンデザインは、経済的発展を前提とした意匠性・快適性の高い空間/形態の創出が主であったが、今後は、極端気象インパクトに対するレジリエンスやwith/after コロナ時代への順応性にも配慮した計画でなければならない。特に、建築物の形状・配置・密度から特徴づけられる都市形態は、空気の流れと相互に影響を及ぼし合うことから、台風等による暴風・強風リスク、猛暑時の暑熱リスク、ウイルス拡散時の感染リスク等に密接に関わっている。
したがって、風工学的視点から新たな都市形態を模索し、防災・衛生工学的知見を蓄積することは、
サステイナブルな都市デザイン手法の社会実装に繋がるものと期待できる。
ここでは，開発が急速に進行している「渋谷地区」をベースエリアとして選定する。多角的な視点
から将来あるべき都市形態のパターンを考案し、LES で得られた非定常乱流場に基づいて、都市の強風防災や換気性能について考察する。
(b)-5-1

特定街区を対象とした変容モデル（その１）

渋谷駅周囲の超高層建物が林立するエリアに着目し、その近傍における特定の中低層密集街区を対象として、当該街区に含まれる複数棟の建物群を単一の建築物に置換する。ここでは、既存街区に- 18 -

対して系統的に形態を変化させる「変容モデル」と位置付ける。対象エリアは、周囲にストリーム、
スクランブルスクエア、ヒカリエ等が位置する三角形の街区（渋谷駅東側）であり、建物群の体積を保持するとともに既存街区形態を反映させて２種の変容パターンを作成した。モデル(1)-1 は既存建物群を高層建物に集約したものであり、元の街区ラインを反映して三角柱形状とした。高層建物に集約したことによりその周囲には大規模なオープンスペースが設けられている。モデル(1)-2 は街区の外形ラインに沿って外壁面を有する低層建物であり、建築計画的配慮から中庭型のオープンスペー
スを設けている。
CUBE による計算結果に基づき、高さ 80m および 20m における各ケースの最大瞬間風速分布を検討した。現状モデルの複雑な街区形態に比べ、高層建物に集約されたモデル(1)-1 では建物側方において強い剥離流が形成されるため、地上付近における建物周囲のオープンスペースで強風域が発生す
るが、同時に建物背後で減速されることから、高さ 80m では後流域に位置する建物周囲で現状モデ
ルに比べ増速域が縮小される事が確認された。中庭型低層建物に置換されたモデル(1)-2 では、高さ80m では現状モデルと概ね同様の分布であるが、高さ 20m では後流域の道路における増速が顕在化する結果となった。要因として、当該エリアの建物高さの均一化により、その上部の流れがスムーズになったためと考えられる。以上より、街区形態の変容パターンに応じて、その周囲において強風防災
上注意すべき空間が顕著に異なることが確認された。特に形態が変化する街区の側方や後方においては、既存建物の耐風性能だけでなく、空間用途に応じて居住環境障害を新たに発生させない配慮が重要となる。
(b)-5-2

特定街区を対象とした変容モデル（その２）

渋谷駅周囲における地域特性を考慮した変容エリアとして、同駅南側に位置する、旧来からの桜並
木の名所「さくら坂」と、広域再開発が進む「桜丘エリア」に挟まれた街区を選定した。当該街区は、渋谷駅から約 100m の距離にありながら 1981 年（新耐震基準制定）以前の建物が過半数を占めているため現実的な開発ニーズが想定されるとともに、周囲を超高層建築物（含 建設中）が囲んでいる
ため、風工学的視点での防災・環境計画が重要となることが推察される。ここでは、当該エリアに対し、中低層建物による構成を前提として以下の各変容パターンが考案された。
【A シリーズ：コの字型プラン】
中規模の中庭空間をコの字平面の建物で囲む基本プラン A-1 に対し建物の高低差を設ける。高い
建物について、主風向（南北）方向に直交（壁状に）配置させる A-2、平行に配置する A-3 を計画した。
【B シリーズ：主風向平行プラン】
主風向に沿う方向で、細長い平面を持つ複数の建物を配置する基本プラン B-1 に対し、建物に高低差を設ける B-2 を計画した。
【C シリーズ：中央大型空地プラン】
当該街区の中央に大規模なオープンスペースを設けるものとし、大まかなボリュームスタディを行った C-1、空地形状を楕円形として建物を分散させた C-2、周囲の高層建物から吹き降ろされた気流を積極的に導入するため C-2 壁面の一部を斜面とした C-3 を計画した。

- 19 -

上記変容パターンより、特に都市の換気性能に着目して、現状／A-2／A-3／C-3 の各モデルについて CUBE による計算を実施した。風向は夏期の主風向 S とした。計算結果に基づく平均流れ場に対し、地上約 35m で発生させた各ケースの流線図を検討した。当該エリアにおいては、上流側の２棟（インフォスタワー、桜丘 B 棟）によって縮流された流れが作用するとともに、下流側の大規模建物（桜丘 A 棟）の左右に分かれて流下してゆく点は共通しているが、各ケースによって変容街区内における流れのパターンに大きく変化が生じていることが確認された。現状モデルでは比較的スムー
ズに街区内を流れが通過する一方、流れに直交して高い建物が配置されている A-2 では水平・鉛直
方向ともに明瞭な渦構造が、高い建物が平行配置されている A-3 では建物間の縮流が、それぞれ顕在化する結果となった。中央に楕円形のオープンスペースを有する C-3 では風上建物からの剥離流が空地内に入り、曲面上を沿う流れが支配的となった。
今後の課題として、これら多様な流れ場を対象として都市の換気性能を適切に定量評価できる手法を構築するとともに、効果的に換気を促進できる都市形態を追及する。
参考文献[2-1] Masaharu Kawaguchi,

Tetsuro Tamura, Hidenori Kawai: Analysis of tornado and near-ground turbulence using a hybrid meteorological model/engineering LES method. Int. J. Heat
Fluid Fl., Vol. 80, No. 108464, 2019.
[2-2] 「富岳」成果創出加速プログラム「防災・減災に資する新時代の大アンサンブル気象・大気環
境予測」実績報告書, 2022.
[2-3] Niclas Jansson, Rahul Bale, Keiji Onishi, Makoto Tsubokura: CUBE: A scalable
framework for large-scale industrial simulations, Int. J. High Perform. Comput. Appl., Vol.
33 (4), 678-698, 2019.
[2-4] 川口真晴,田村哲郎,益子渉: 気象擾乱構造を再現した LES による顕著な台風時の実都市の強風評価, 第 34 回数値流体力学シンポジウム, E06-5, 2021.
[2-5] 河合，田村：メソ気象モデルにより再現された風速場への高周波変動成分付加方法の提案 空間フィルタリング・リスケーリング手法の適用-，日本建築学会構造系論文集，第 85 巻，第 767 号，pp.19-27，2020[2-6] Nishi, A. and Kusaka, H.：Effect of Foehn Wind on Record-Breaking High Temperature
Event (41.1°C) at Kumagaya on 23 July 2018, SOLA, Vol.15, (2019), pp. 17-21.
[2-7] 新井舞子，河合英徳，田村哲郎：2018 年猛暑日における東京都市域での弱風・温熱場予測に関する研究，第 35 回数値流体シンポジウム，2021(3) サブ課題 C「新興ウイルス感染症にロバストで健康・快適・サステイナブルなポストコロナ時代の室
内環境設計」(a) 研究成果の要旨
本年度は、統合シミュレーション技術の開発（人体幾何形状と呼吸器系幾何形状を精緻に再現した- 20 -

数値人体モデルの開発）、ベンチマーク問題の設定と解析モデルの作成、統合シミュレーションの実
施、多目的最適化フレームワークの設計、及び基盤アプリの「富岳」チューニングを行った。具体的には、鼻腔･口腔から気管支までの上気道モデル、気管から気管支第 16 分岐までの下気道モデル、
ガス交換を担う肺胞サブモデルを個別に開発した上で、各モデルを統合した数値気道プロトタイプ
モデルの開発を進めると共に、室内環境における感染者と気道モデルを取り込んだ被感染者間での感染リスク統合シミュレーションのベンチマークテストの実施を行った。尚、多目的最適化フレーム
ワークについては、サブ課題 A と連携して進め、フレームワーク自体は完成した（成果の詳細はサブ課題 A を参照）。
In this fiscal year, we developed (i) integrated simulation methods (in silico humanmodel that precisely reproduces human body morphology and respiratory tract geometry) and(ii) a multi-objective optimization framework, (iii) benchmark tests and analytical models,and conducted (iv) integrated numerical analyses and, and (v) these were implemented andoptimized on the supercomputer “Fugaku”. Specifically, an upper airway model from thenasal /oral cavities to the bronchi, a lower airway model from the trachea to the 16th
bronchial tubes, and an alveolar sub-model responsible for gas exchange were developed
separately, and a prototype in silico airway model integrating each model was developed.
Benchmark testing of an integrated numerical analysis to predict the airborne transmission
and infection risk was conducted for simple indoor environmental condition. In addition,a benchmark test of infection risk assessment in a music hall was also conducted. The
framework of multi-objective optimization framework was developed in collaboration with
sub-project A (For the detail of the achievements, refer to sub-project A).
(b) 本研究の内容
本年度(2021 年度)の研究成果を項目別に示す。
(b-1) 上気道から下気道までを再現した数値気道モデル(プロトタイプモデル)の作成
九州大学医系地区部局観察研究倫理審査委員会にて承認済みの胸部 CT データの 2 次利用研究「胸部 CT データを用いた in silico 人体肺モデルの精度向上をめざす研究」(許可番号 21101-11)にて匿名化された胸部の高精細 CT データの提供を受け、鼻腔・口腔から気管支第 16 分岐までを含む上気道･下気道統合型の in silico 肺モデルを作成した。
胸部 CT データは成人男性の上気道から全肺を含む高分解能 CT(0.5mm 厚) であ り、 Mimics(Materialize) に て DICOM デ ー タ よ り 気 道 幾 何 形 状 の 抽 出 と 表 面 ス ム ー ジ ン グ 、 引 き 続 い て
ANSYS/Meshing にて流体解析用のメッシュ生成を行った。上気道から下気道気管支第 16 分岐まで再現した気道モデルの形状を Fig. 2.2.3-1 に示す。また、定常呼吸 7.5L/min を鼻孔面に与えた条件で気道内流れ場解析を行った結果を Fig. 2.2.3-2 に示す。PIV による測定結果(Phuong and Ito, 2015)と比較することで、SST k-ωモデルもしくは低 Re 型 k-εモデル(Abe-Nagano-Kondo モデル)の 2 種
類の乱流モデルを用いた解析結果は十分な精度で流れ解析が可能であることを確認した。また、定常
呼吸 7.5L/min の条件で、鼻孔面位置で dp=10μm の粒子を 10,000 個発生させた場合の気道内粒子沈- 21 -

着分布の解析結果を Fig. 2.2.3-3 に示す。上気道では鼻腔ならびに咽頭部で粒子沈着量が多くなる傾向があり、下気道では分岐位置にて沈着量が多くなる傾向が確認できる。

Fig. 2.2.3-1 in silico Airway geometry.

Fig. 2.2.3-2 Airflow profile at bronchial pathway.

Fig. 2.2.3-3 Particle deposition distributions in airway (dp=10μm).

- 22 -

Fig. 2.2.3-4 Deposition fraction (%) as a function of Inertia Parameter.
呼吸流量を変化させた場合の吸入粒子の上気道に対する沈着割合(Deposition fraction)を算出した上で、Kelly らの既往研究と比較した結果を Fig. 2.2.3-4 に示す。本研究で作成した in silico 気
道モデルを用いた気道内粒子沈着量予測結果は in vivo ならびに in vitro 実験結果と良く対応した結果を示している。
上気道から下気道(気管支第 16 分岐)までを再現して作成した新たな in silico 気道モデルのプロ
トタイプモデルは、気道内流れ場ならびに粒子拡散場解析の点で、既往研究による実験結果を非常に良く対応する結果となり、今後の経気道曝露評価への展開に十分使用可能であることを確認した。
(b-2) ガス交換を再現する肺胞サブモデルの作成
本年度(2021 年度)は(b-1)で開発した in silico 気道モデルの第 16 分岐以降を再現するサブモデ
ルとして、酸素―二酸化炭素ガス交換を再現する肺胞サブモデルを作成した。肺胞は第 23 分岐位置にあると想定し、16 分岐から 23 分岐位置までの死腔域は単純に移流方程式で空気輸送の時間遅れを再現した。肺胞サブモデルでは、CO₂分圧の時間変化を再現するため、肺胞毛細血管からの拡散を考慮すると共に、呼吸に伴う肺胞体積の変化をモデル化して作成した。肺胞サブモデルと in silico 気
道モデルとの統合の概要を Fig. 2.2.3-5 に示す。また、呼吸量変化に伴う肺胞体積変化のモデル化結果を Fig. 2.2.3-6 に示す。
肺胞 CO2 分圧のマスバランス式、肺胞 CO2 分圧の拡散体積、肺毛細管内の CO2 分圧分布、の基本
モデル式を以下に示す。

VA (t )
・dPA (t ) Qres (t ) PA (t )  ( Patm  PH 2O )(Vdif (t )  Qres (t )Cres (t ))
dt

(3-1)

DL n
 ( Pc ,i (t )  PA (t ))n i 1(3-2)

V dif - 23 -

Vc ,1
Vc ,idPc ,1 (t )dtdPc ,i (t )
dt



DL
Pc,1 (t )  PA (t )   Qblo Pc,v (t )  Pc,1 (t )
n

(3-3)



DL
Pc,i (t )  PA (t )   Qblo Pc,i 1 (t )  Pc,i (t ) i = 2, 3, …, n
n

(3-4)

Fig. 2.2.3-5 Outline of gas-exchange model at alveolus.

Fig. 2.2.3-6 Time series of alveolar volume according to the transient breathing profile.
肺胞サブモデルを in silico 気道モデルに実装した上で、人工気候室内で人体由来の CO2 発生量を測定した先行実験(Kuga et al, 2021)を対象として数値解析を実施した。呼吸１サイクル分の非定常解析結果を Fig. 2.2.3-7 に示す。CO2 濃度約 400ppm の室内空気を鼻孔から吸入し、気道内輸送の後に肺胞サブモデルでガス交換し、呼出時には 50,000ppm を超える CO2 濃度となり鼻孔面から室内側に吐出される。呼吸２サイクル分に対して、鼻孔面での吸入 CO2 濃度の変化、肺胞モデルでのCO2 分圧変化、肺毛細管の CO2 分圧変化、CO2 拡散体積の変化に関する解析結果を Fig. 2.2.3-8 に示す。

- 24 -

Fig. 2.2.3-7 Time series of distribution of CO2 concentration exhaled by the CSP model
connected with the airway.

Fig. 2.2.3-8 Results of integrated analysis with CO2 gas exchange model with and without- 25 -

the airway model. (a) inhaled and exhaled CO2 concentration in nasal openings, (b) CO2
partial pressure in alveoli, (c) CO2 partial pressure in blood flowing to the artery, and (d)CO2 diffusing volume.
本年度(2021 年度)は CO2 ガス交換に着目した肺胞サブモデルを作成した。
(b-3) 室内環境中での空気感染リスク解析のベンチマークテスト
感染者からの咳飛沫・飛沫核の拡散が対面する非感染者へ曝露リスク解析の予備的な解析として、
3m 立方の単純室内に人体表面幾何形状と気道形状を再現した 2 体の数値人体モデル(in silico 人体モデル)が Physical distance(1m と 2m)の間隔で向かい合って立っている状況を想定した(Fig. 2.2.39)。1 体は感染者であり、もう 1 体は非感染者である。室内は、下部前方から Uin=0.05 m/s で外気が流入し、人体熱負荷により加熱された空気が上部後方から自由流出する置換換気方式を想定している。人体発熱は Fanger の熱的中立モデルを適用した。また、感染者が吐出した飛沫・飛沫核による非感染者の吸入曝露量および気道内沈着分布までを連続して解析するため、鼻腔開口面を介して非感染者を対象とした人体モデルに(b-1)で作成した in silico 気道モデルの上気道モデル部分をシームレスに統合した。
咳プロファイルは、Gupta et al の被験者実験による咳の呼出流量の結果を基に cough profile を作成し、口腔開口面に設定した。本予備解析では、飛沫･飛沫核の空気力学径は 1μm~10μm の微粒子を想定し、統計的有意性を鑑みて、対象粒子は口腔開口面で 45,000 個発生させる条件とし、蒸発などによる粒径変化は無視した条件とした。
結果の表示は割愛するが、今回の解析では置換換気を想定したため、人体周辺に熱上昇流および温度成層が形成された。Fig. 2.2.3-10 に感染者から吐出する飛沫･飛沫核の室内分布の 30 秒間の飛跡(流れ場 URANS と粒子 Lagrange 解析の連成結果)を示す。咳によって呼出された粒子は慣性力にて対面する非感染者顔面に向けて輸送される。Fig. 2.2.3-11 には、対面する非感染者の気道壁面への飛沫・飛沫核の沈着分布を示す。感染者が咳をしてから 75 秒後の時点で、感染者と非感染者間の距離
が 1m 離れて定常吸入する場合は、非感染者の人体皮膚表面への沈着割合が 2.89％、気道壁面への沈着割合は 2.64％であった。対して 2m 離れて定常吸入する場合は、非感染者の人体皮膚表面への沈着
割合は 0.016％、気道壁面への沈着割合は 0.13％となり、経気道感染の抑制には物理的距離の確保が重要となることを確認した。

Fig. 2.2.3-9 Room model Analyzed.

- 26 -

(a)(b)

Fig. 2.2.3-10 Particle Path lines for 30sec. (a) 1m physical distance, (b) 2m physicaldistance.

Fig. 2.2.3-11 Particle deposition distributions in upper airway.

(4) サブ課題 D「マルチコンポーネント統合シミュレーションによるカーボンフリーガスタービン設計」(a) 研究成果の要旨
本年度は、対象とするガスタービンの選定・詳細確認、圧縮機・燃焼器・タービンの各コンポーネントの解析実施（試計算）、乱流モデルの調査、選定、及び水素、アンモニア燃焼詳細機構の調査および 1 次元解析による精度検証を行った。その結果、川崎重工業製ガスタービンを対象に、圧縮機、燃焼器、タービン、各コンポーネントの計算メッシュを作成し、試計算を実施した。また、水素およびアンモニアの燃焼詳細機構の違いが、それらの燃焼特性に及ぼす影響を、1 次元解析により明らかにすることが出来た。
In this fiscal year, the selection of the target gas turbine, the confirmation of its
details, the trial numerical simulations of each component, namely compressor, combustor
and turbine, the investigation and selection of turbulence model, and the examinations of- 27 -

validity and accuracy of the detailed reaction mechanisms for hydrogen and ammonia
combustion using one-dimensional numerical simulations were conducted. The trial numericalsimulations based on large-eddy simulations (LES) were successfully carried out, and theone-dimensional numerical simulation showed how the difference in the reaction mechanism
affected the combustion characteristics of hydrogen and ammonia.
(b) 本研究の内容(b-1) ガスタービン解析
本解析では、非構造格子に適用可能な FrontFlow /red large-eddy simulation (LES)を、京都大学と（株）数値フローデザインにより改良したものを用い[4-1,4-2]、本年度は、簡易メッシュを用いた試計算を実施した。
1. 計算対象
Fig. 2.2.4-1 に、対象とした川崎重工業株式会社の産業用ガスタービンを示す。本ガスタービン
は、圧縮機、燃焼器、タービンの 3 つのコンポーネントから構成される。以下に、各コンポーネントの概要を示す。

Fig. 2.2.4-1 Target gas turbine.
1-1 圧縮機
圧縮機の概形およびメッシュを Fig. 2.2.4-2 に示す。また、1 段目のインペラの表面メッシュをFig. 2.2.4-3 に、2 段目のインペラの表面メッシュを Fig. 2.2.4-4 に示す。圧縮機の解析格子の節
点数、要素数はそれぞれ 0.12 億、0.71 億である。

- 28 -

Fig. 2.2.4-2 Geometry and mesh of compressor.

Fig. 2.2.4-3 Mesh on 1st impellor of compressor.

- 29 -

Fig. 2.2.4-4 Mesh on 2nd impellor of compressor.
1-2 燃焼器
燃焼器の概形およびメッシュを Fig. 2.2.4-5 に示す。
また、メインバーナへの燃料供給系統を Fig.
2.2.4-6 に、追焚きバーナの形状を Fig. 2.2.4-7 に、追焚きバーナの中心を通る断面のメッシュをFig. 2.2.4-8 に示す。燃焼器の解析格子の節点数、要素数はそれぞれ 0,23 億、1.3 億である。

Fig. 2.2.4-5 Geometry and mesh of combustor.

- 30 -

Fig. 2.2.4-6 Fuel supply system of combustor.

Fig. 2.2.4-7 Supplemental burner

Fig. 2.2.4-8 Mesh around auxiliary burner.

- 31 -

1-3 タービン
タービンの翼列および流路内のメッシュを Fig. 2.2.4-9 に、後流側の翼面のメッシュの様子をFig. 2.2.4-10 に示す。タービンの解析格子の節点数、要素数はそれぞれ 0.34 億、1.1 億である。

Fig. 2.2.4-9 Geometry and mesh of turbine.

Fig. 2.2.4-10 Mesh on the exit plane of turbine.
2. 計算方法、モデル、および条件
本解析では、予混合燃焼をベースとした Flamelet Generated Manifold (FGM)法[4-3]を用いた LESを実施する。FGM 法の支配方程式は、以下の連続の式、運動量の保存式、混合分率の保存式、反応進
行変数の保存式、エンタルピーの保存式、および NOx の生成予測のための輸送方程式で構成される。

質量保存式：- 32 -

+

𝜌̅ 𝑢

(1)

=0

運動量保存方程式：
+

̅

+

+

(2)

𝜌̅ 𝐷

+

(3)

=−

混合分率輸送式：+=

反応進行度変数輸送式：+=+ 𝜌̅ 𝑊

(4)

+

+

(5)

𝜌̅ 𝐷

+

𝜌̅ 𝐷

+

エンタルピー輸送式：
+

=

̅

+𝑢

̅

NO 輸送式：+=+ 𝜌̅ 𝑊(6)

ここで、上付きの―は LES フィルタを表し、～はファーブル平均を意味する。上記輸送方程式より求めた流れ場の混合分率及び反応進行度変数を元に、1 次元予混合燃焼解析から得られた flamelet データベースを参照し、温度や、反応速度、化学種組成を求めることで乱流燃焼場の予測を行う。
圧縮機およびタービンの動翼を含む領域は非慣性系として扱い、運動量保存式にはコリオリ力および遠心力を含む。なお、圧縮機およびタービンでは反応が起こらないと仮定し、反応進行度変数輸
送式の生成項を 0 とする。また、タービンは既燃ガスとなる値を初期値と流入値として設定する。
圧縮機、燃焼器、タービンの解析条件を以下に示す。
2-1 圧縮機
圧縮機の解析条件を表 2.2.4-1 に示す。
Table 2.2.4-1 Numerical conditions for compressor.
乱流モデル

LES(標準スマゴリンスキモデル)

移流項の補間法

1 次精度風上

時間積分法

1 次オイラ陰解法

時間刻み

2.0×10－7 秒

壁面の条件

速度：対数則 温度：断熱

流入条件

圧力固定(全圧 0.1MPa)

流出条件

圧力(全圧 1.0MPa)、流量固定(8.0kg/秒)- 33 -

22,000rpm

回転数

2-2 燃焼器
燃焼器の解析条件を表 2.2.4-2 に示す。
Table 2.2.4-2 Numerical conditions for combustor.
乱流モデル

LES(標準スマゴリンスキモデル)

移流項の補間法

2 次精度風上

時間積分法

1 次オイラ陰解法

時間刻み

2.0×10－6 秒

壁面の条件

速度：対数則 温度：断熱

燃料

都市ガス 13A

酸化剤

空気

流入条件

酸化剤：348℃、1MPa、7.4kg/s
燃料：25℃、1MPa
パイロット：1.1×10－3kg/s
メイン：0.1kg/s
追焚き：2.8×10－2kg/s

2-3 タービンタービンの解析条件を表 2.2.4-3 に示す。
Table 2.2.4-3 Numerical conditions for turbine.
乱流モデル

LES(標準スマゴリンスキモデル)

移流項の補間法

1 次精度風上

時間積分法

1 次オイラ陰解法

時間刻み

2.0×10－7 秒

壁面の条件

速度：対数則 温度：断熱

流入条件

圧力固定(全圧 1.0MPa)

流出条件

圧力(全圧 0.1MPa)、流量固定(8.0kg/秒)

回転数

22,000rpm

2-4 並列計算の方法と効果（性能）
並列計算はフラット MPI による。圧縮機は 1920 並列、160 ノード(12CPUs/ノード)で解析を行い、1,000 ステップあたりの計算時間は 42 分であった。燃焼器は 2,400 並列、100 ノード(24CPUs/ノード)、で解析を行い、1,000 ステップあたりの計算時間は 30 分であった。タービンは 9,600 並列、800
ノード(12CPUs/ノード)で解析を行い、1,000 ステップあたりの計算時間は 2 時間であった。

- 34 -

3. 研究成果
3-1 圧縮機
Fig. 2.2.4-11～2.2.4-13 に、圧縮機の回転軸を通る断面、1 段目のインペラの中心を通る断面、および 2 段目のインペラの中心を通る断面の流れ場の様子をそれぞれ示す。詳細な評価、精度検証は実施中であるが、安定して計算を行えることが確認できている。

Fig. 2.2.4-11 Flow field in compressor (cross section on the rotation axis).

Fig. 2.2.4-12 Flow field in compressor (cross section on the center of 1st impellor).

- 35 -

Fig. 2.2.4-13 Flow field in compressor (cross section on the center of 2nd impellor).
3-2 燃焼器
Fig. 2.2.4-14 に、燃焼室内の瞬時燃焼場の様子を示す。メイン燃料は、燃焼器上流部のダクト内で吸入空気と混合し比較的均一な燃料濃度で燃焼室内へ流入している。二次空気流入位置より上流
で希薄条件となっているが、保炎されていることが確認できる。二次空気及び追焚き用の燃料の流入により燃焼室出口では混合分率がさらに小さな値となっているが、これは実機の想定とも一致する。

- 36 -

Fig. 2.2.4-14 Instantaneous reacting flow field in combustor.
3-3 タービン
Fig. 2.2.4-15 に、タービンの回転軸を通る断面の流れ場の様子を示す。詳細な評価、精度検証は実施中であるが、安定して計算が行えることが確認できている。

- 37 -

Fig. 2.2.4-15 Flow field in turbine.
4. まとめと今後の課題
本研究では、川崎重工業株式会社の産業用ガスタービンを対象に、それを構成する圧縮機、燃焼器、タービン、3 つのコンポーネントについて、それぞれ個別に LES 解析を実施した。燃焼器に関しては、
希薄燃焼条件で保炎すること、実機で想定された当量比を再現することを確認した。また、圧縮機、タービンに関しては、安定に計算を実行可能であることを確認した。
(b-2) 水素およびアンモニアの燃焼評価
1. 計算方法、モデル、および計算条件
1-1 計算方法
本解析では、水素燃焼およびアンモニア燃焼用の詳細反応機構の精度評価を行う。着火遅れ時間、
層流燃焼速度、消炎伸長率の評価には、それぞれ定積断熱自着火の 0 次元数値解析、予混合自由伝播
火炎の 1 次元数値解析、予混合対向流火炎の 1 次元数値解析を行う。また、いずれの数値解析においても、計算コードとして FlameMaster [4-4]を用いる。
1-2 計算条件
水素燃焼においては、Konnov 反応機構[4-5]（14 化学種 75 段反応）
、Conaire 反応機構[4-6]（9 化
学種 19 段反応）
、Miller 反応機構[4-7]（9 化学種 20 段反応）
、および Zhou et al. [4-8]と Glarborg
et al.[4-9]の詳細反応機構の組み合わせから、C 原子を含む化学種と反応を除いた詳細反応機構[410]（33 化学種 208 段反応、以降、ZG 反応機構と呼ぶ）の 4 種類の詳細反応機構について、着火遅れ
時間、燃焼速度を実験値と比較し精度評価を行う。定積断熱自着火の 0 次元数値解析においては、当量比 0.1（H2: 0.81 vol%, O2: 4.03vol%, Ar: 95.16vol%）と当量比 0.5（H2: 3.47vol%, O2: 3.47vol%,N2: 93.06vol%）の 2 種類の組成について、初期圧力を 1atm、初期温度を 950K から 2000K まで変化

- 38 -

させて解析を行う。予混合自由伝播火炎の 1 次元数値解析においては、圧力を 1atm、初期温度を 300Kとし、当量比を 0.3 から 2.2 まで変化させて解析を行う。
アンモニア燃焼においては、Otomo 反応機構[4-11]（31 化学種 204 段反応）
、Okafor 反応機構[412]（59 化学種 356 段反応）、GRI-Mech3.0[4-13]（53 化学種 325 段反応）
、UCSD 反応機構[4-14]（70
化学種 318 段反応）の 4 種類の詳細反応機構について、着火遅れ時間、層流燃焼速度、消炎伸長率を実験値と比較し精度評価を行う。定積断熱自着火の 0 次元数値解析においては、初期圧力 1.4、 11,30atm の 3 条件について、当量比を 1.0 として、初期温度を 1600K から 2500K まで変化させて解析を行う。予混合自由伝播火炎の 1 次元数値解析においては、圧力を 0.1MPa、初期温度を 300K とし、当量比を 0.8 から 1.2 まで変化させて解析を行う。予混合対向流火炎の 1 次元数値解析においては、
圧力 0.1、 0.5MPa の 2 条件について、当量比 0.8、 1.0、 1.2 の 3 条件を考慮し、初期温度を 300Kとし、火炎伸長率をパラメータとして増加させ、消炎するまで解析を行う。
2. 研究成果
2-1 水素燃焼
Fig. 2.2.4-16 に、4 種類の詳細反応機構を用いて得られた水素の着火遅れ時間を実験値[4-15]と比較して示す。図より、当量比 0.1 においては、いずれの詳細反応機構でも、実験値を精度良く予測することが分かる。一方、当量比 0.5 においては、10000/T が小さい、すなわち高温条件下において、Miller 反応機構と ZG 反応機構は実験値を過大評価することが分かる。また、Conaire 反応機構とKonnov 反応機構を比較すると、Konnov 反応機構の方が実験値を精度良く予測することが分かる。
Fig.
2.2.4-17 に、数値解析により得られた水素の層流燃焼速度を実験値[4-16～4-19]と比較して示す。
図より、Miller 反応機構は実験値を過大評価し、その他の詳細反応機構は実験値を精度良く予測することが分かる。以上の結果より、水素/空気予混合燃焼を対象とする数値解析においては、Konnov
反応機構、または Conaire 反応機構が適すると考えられる。
2-2 アンモニア燃焼
Fig. 2.2.4-18 に、4 種類の詳細反応機構を用いて得られたアンモニアの着火遅れ時間を実験値[420]と比較して示す。図より、いずれの圧力条件下においても、GRI-Mech3.0 は実験値を大きく過大評価することが分かる。その他の詳細反応機構は僅かに実験値を過大評価するものの、比較的精度は高く、本解析で比較した 4 種類の詳細反応機構の中では、Otomo 反応機構の精度が最も高いことが分かる。Fig. 2.2.4-19 に、数値解析により得られたアンモニアの層流燃焼速度を実験値[4-21]と比較して示す。図より、UCSD 反応機構は実験値を過大評価し、その他の詳細反応機構は同程度の精度で層流燃焼速度を予測することが分かる。以上の結果より、Otomo 反応機構と Okafor 反応機構は着火
遅れ時間と層流燃焼速度の両方を精度良く予測することが分かる。そのため、Fig. 2.2.4-20 では、
これら二つの詳細反応機構を用いて得られた消炎伸長率を実験値[4-22]と比較する。図より、いずれの圧力条件においても Otomo 反応機構は実験値を精度良く予測し、Okafor 反応機構は実験値を過大評価することが分かる。以上の結果より、アンモニア/空気予混合燃焼を対象とする数値解析においては、着火遅れ時間、層流燃焼速度、消炎伸長率いずれについても実験値を精度良く再現する Otomo
反応機構が適すると考えられる。

- 39 -

Fig. 2.2.4-16 Ignition delay time (H2).

Fig. 2.2.4-17 Laminar burning velocity (H2).

Fig. 2.2.4-18 Ignition delay time (NH3).

Fig. 2.2.4-19 Laminar burning velocity (NH3).

Fig. 2.2.4-20 Extinction flame stretch rate (NH3).

3. まとめと今後の課題
本研究では、水素燃焼およびアンモニア燃焼の 0 次元および 1 次元の数値解析により、最適な詳- 40 -

細反応機構を選定した。
参考文献[4-1] H. Moriai et al., J. Eng. Gas Turb. Pow., 135, 091503, (2013).
[4-2] S. Tachibana et al., Combust. Flame, 162, 2621-2637, (2015).
[4-3] J.A. van Oij, and L.P.H. de Goey, Combust. Sci. Tech., 161, 113-137, (2000).
[4-4] H. Pitsch, FlameMaster, a C++ computer program for 0D combustion and 1D laminarflame calculations.
[4-5] A. A. Konnov, Combust. Flame, 203, 14-22, (2019).
[4-6] M. Ó Conaire et al., Int. J. Chem. Kinet., 36, 603-622, (2004).
[4-7] J. A. Miller, C. T. Bowman, Prog. Energy and Combust. Sci., 15, 287-338, (1989).
[4-8] C. W. Zhou et al., Combust. Flame, 197, 423-438, (2018).
[4-9] P. Glarborg et al., Prog. Energy and Combust. Sci., 67, 31-68, (2018).
[4-10] A. Durocher et al., Proc. Combust. Inst., 38, 2093-2100, (2021).
[4-11] A. Otomo et al., Int. J. Hydrog. Energy, 43, 3004-3014, (2018).
[4-12] E. C. Okafor et al., Combust. Flame, 187, 185-198, (2018).
[4-13] Gregory P. Smith et al., http://www.me.berkeley.edu/gri_mech/[4-14] San Diego Mechanism web page, http://combustion.ucsd.edu.
[4-15] A. Kéromnès et al., Combust. Flame, 160, 995-1011, (2018).
[4-16] D. R. Dowdy et al., Proc. Combust. Inst., 23, 325-332, (1991).
[4-17] O. C. Kwon, G. M. Faeth, Combust. Flame, 124, 590-610, (2001).
[4-18] A.E. Dahoe, J. Loss Prev. Process Ind., 18, 152-166, (2005).
[4-19] M. Kuznetsov et al., Int. J. Hydrog. Energy, 37, 17580-17588, (2011).
[4-20] O. Mathieu, E. L. Petersen, Combust. Flame, 162, 554-570, (2015).
[4-21] A. Hayakawa et al., Fuel, 159, 98-106, (2015).
[4-22] S. Colson et al., J. Therm. Sci. Technol., 11, JTST0048, (2016).
(5) プロジェクトの総合的推進
以上のサブ課題 A～D の進捗を管理し、課題間の連携を促進するために、各課題の代表者を集めて統括する推進委員会を設置した。推進委員会は、３か月に一度程度の頻度で理研が主催して委員会を開催した。当初、年に一度の成果報告会や、他の「産業競争力の強化」課題と連携した統合ワークシ
ョップの企画・開催を計画したが、本年度は成果報告会については単独で実施せず、「富岳」成果創
出加速プログラムシンポジウム並びに技術交流会で報告を行った。
２－３．活動（研究会の活動等）(1) アウトリーチ活動
年 月 日

行

- 41 -

事

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

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

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

行

事

R3.9.7

第 1 回「富岳」スマートデザインプロジェクト推進会議

R3.12.6

第 2 回「富岳」スマートデザインプロジェクト推進会議

R3.3.8

第 3 回「富岳」スマートデザインプロジェクト推進会議(3) コンソーシアム活動
年 月 日

行

事

R3.10.8

都市・建築 CFD コンソーシアム 総会

R3.11.5

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

R3.12.24

HPC を活用した自動車用次世代 CAE コンソーシアム(4) ワーキンググループ
年 月 日

行

事

R3.10.14,
R3.11.10,
R4.1.11,

サブ課題Ａ サブ課題全体研究打合せ（全 5 回）R4.2.21,R4.3.31R3.9.3,R3.9.10,R3.9.17,R3.9.24,
R3.10.1,
R3.10.8,

サブ課題Ａ 最適化打ち合わせ（全 18 回）

R3.10.29,
R3.11.5,
R3.11.26,
R3.12.3,
R3.12.10,- 42 -

R4.1.7,
R4.1.17,R4.1.28,
R4.2.4,R4.2.18,
R4.3.11,R4.3.25
R3.9.17,R3.9.24,
R3.10.8,
R3.10.15,
R3.10.22,
R3.10.29,
R3.11.12,
R3.11.19
R3.11.26,

サブ課題Ａ デザイン・空力モデル化打ち合わせ（全 17 回）

R3.12.3,
R3.12.10,
R4.1.7,R4.1.17,R4.1.28,
R4.2.4,R4.2.18,
R4.3.28
R3.11.5,
R3.12.20,
R4.2.7,

サブ課題 B サブ課題全体研究打合せ（全 4 回）R4.3.7
R3.11.4,
R3.12.3,
R4.1.14,

サブ課題 B 極端気象流入変動風整備 WG（全 5 回）R4.2.22,R4.3.30
R3.10.4,
R3.11.4,
R3.12.8,

サブ課題 B 台風時の都市・建築における機能性・安全性 WG（全 6 回）

R4.1.21,R4.2.28,- 43 -

R4.3.31
R3.10.13,
R3.11.24,
R3.12.17,

サブ課題 B 猛暑時の熱環境推定・安全性指標 WG（全 5 回）R4.2.24,R4.3.23
R3.10.21,
R3.11.19,
R3.12.24,

サブ課題 B 早期発展 WG（全 5 回）R4.1.27,
R4.3.4R3.9.9,R3.9.16,R3.9.24,
R3.9.30,
R3.10.15,
R3.10.21,
R3.11.24,
R3.12.23,

サブ課題 C サブ課題全体研究打合せ（全 14 回）

R4.1.6,R4.1.13,R4.2.10,R4.2.24,R4.3.7,R4.3.25
R3.9.2,R3.10.6,
R3.11.12,

サブ課題 C CREST 連携会議（全 5 回）

R3.12.8,
R4.1.27
R3.9.8,R3.10.21,
R3.11.18,
R3.12.20,

サブ課題 D サブ課題全体研究打合せ（全 7 回）R4.1.21,R4.2.7,R4. 3.14- 44 -

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

担当機関

担当責任者

（１）サブ課題 A「意匠空間理化学研究所計算科学研究センを考慮した AI 支援多目的最ター適化による自動車空力デザイ

神戸大学大学院

教授 坪倉誠広島大学大学院

助教 中島卓司

東京工業大学大学院教授 坂田弘安

チームリーダー 坪倉誠ン」（２）サブ課題 B「変容する都市・建築の自然擾乱対応の性能設計」

研究員 田村哲郎
横浜国立大学大学院

教授 山田貴博

近畿大学

教授 道岡武信

お茶の水女子大学

准教授 河合英徳（３）サブ課題 C「新興ウイ九州大学総合理工学研究院

教授 伊藤一秀

ルス感染症にロバストで健

理化学研究所計算科学研究セン

チームリーダー 坪倉誠

康・快適・サステイナブルなター

ポストコロナ時代の室内環境

神戸大学大学院

教授 坪倉誠（４）サブ課題 D「マルチコ京都大学大学院教授 黒瀬良一

ンポーネント統合シミュレー

九州大学総合理工学研究院

教授 渡邊裕章

設計」

ションによるカーボンフリー
ガスタービン設計」- 45 -

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

発表者氏名発表した場所

発表し

国内・（学会等名）た時期

外の別第49回可視化情報シンポジ

2021 年「富岳」を用いた飛沫・エ
アロゾルの可視化と新型コ
ロナ感染リスク評価による

坪倉 誠

ウム

9月

国内

社会啓発（特別講演）
スパコンによるウイルス飛沫エアロゾルの飛散予測とそのリスク低減対策（招待第23回日本咳嗽学会学術大

坪倉 誠

会

2021 年
9月

国内

講演）「富岳」で目指すSociety
5.0時代のスマートデザイン

日本機械学会第34回計算力

坪倉 誠

学講演会

（招待講演）

2021 年
9月

国内

スパコン「富岳」による飛沫・エアロゾル感染リスク
評価のデジタルトランスフ

第38回「センサ・マイクロ坪倉 誠

マシンと応用システム」シォーメーション（招待講

ンポジウム

2021 年
11月

国内

演）
The International
HPC Response to the
COVID-19 Pandemic(Panelist)

Conference for High

MakotoTsubokura

Performance Computing,
Networking, Storage, and2021 年
11月

国外

Analysis (SC21)
Digital Transformation ofDroplet/Aerosol Infection
Risk Assessment Realized
on “Fugaku” for the Fightagainst COVID-19(Finalist of the ACM
Gordon Bell Special Prize
for High Performance
Computing-Based COVID-

Kazuto Ando,
Rahul Bale,The International

ChungGang Li,

Conference for High

Satoshi

Performance Computing,Matsuoka, Keiji

Networking, Storage, andOnishi, Makoto

Analysis (SC21)

2021 年
11月

国外

Tsubokura

19 Research
スパコン富岳による飛沫エ
アロゾル感染リスク評価の情報処理学会第84回全国大

坪倉 誠

会

デジタルトランスフォーメ- 46 -

2022 年
3月

国内

ーション（招待講演）
Digital transformation of
the manufacturing process
and smart design in theSociety 5.0 era realized on
the supercomputer

MakotoTsubokura

Supercomputing Asia 2022

2022 年
3月

国外“Fukuga” (Opening
Plenary talk)
スパコン「富岳」によるウ
イルス飛沫エアロゾル拡散
シミュレーションと感染リ

第92回日本衛生学会学術総

坪倉 誠

会

2022 年3月

国内

スク評価（招待講演）
スパコン富岳による飛沫・
エアロゾル感染リスク評価のデジタルトランスフォー

坪倉 誠

日本薬学会第142年会

2022 年
3月

国内

メーション（招待講演）「富岳」を用いた飛沫・エ
アロゾル飛沫解析と新型コ
ロナ感染リスク評価による

産総研構造設計コンソーシ

坪倉 誠

アム

2021 年
9月

国内

社会啓発
富岳が拓くSociety 5.0時代のスマートデザイン

VINAS Online Users

坪倉 誠

Conference 2021

スパコン「富岳」による飛沫・エアロゾル感染リスク
評価のデジタルトランスフ

第38回「センタ・マイクロ坪倉 誠

マシンと応用システム」シンポジウムォーメーション
スパコン「富岳」による新型コロナウイルス飛沫感染

坪倉 誠

第20回全国福祉村サミットの予測と対策坪倉 誠

革新的ものづくり

Assessment System for
the Viral Droplet Infection

2021 年
11月

2021 年
11月

国内

国内

国内

記念講演, 国際電気通信電

2021 年

機連合（ITU）「世界情報社

11月

国内会・電気通信日のつどい」

Development of theIntegrated Risk

10月

日本ITU協会賞特別賞受賞「富岳」を用いた新型コロ
ナウイルス対策への貢献と2021 年

JST/STI for Global
Challenges: International

MakotoTsubokuraon a Supercomputer andResearch Collaboration
against the COVID-19
Crisis- 47 -

2021 年
12月

国外its Social Implementation
スーパーコンピュータ「富全日本科学技術協会 新型岳」によるウイルス飛沫・コロナウイルス感染防止にエアロゾル感染リスクの評

坪倉 誠

向けた地域プラットフォー

価と対策提案 ～飲食店でム形成シンポジウム第４回の感染対策を例に～

ワークショップ

2021 年
12月

国内

計算科学振興財団・神戸商
富岳が拓くSociety 5.0時代のスマートデザイン

坪倉 誠

工会議所主催, スーパーコ

2022 年

ンピュータ・ソリューショ

2月

国内

ンセミナー
日本自動車技術会流体技術

Society 5.0時代のものづく
りと新型コロナ感染リスク
評価のデジタルトランスフ

部門委員会シンポジウム
坪倉 誠「自動車開発を支える最新の空力技術および空力音響ォーメーション

のデジタルトランスフォー

3月

国内

技術」

富岳を用いたウイルス飛沫
エアロゾル感染リスク評価

2022 年（一社）センサイト協議会
坪倉 誠

主催 センサイトWEBセミ
ナー

メーション

2022 年
3月

国内

スーパーコンピュータ「富岳」によるウイルス飛沫・
エアロゾルのシミュレーシ

坪倉 誠

ョンと新型コロナ感染防止

日本機械学会関西支部シニ

2022 年

ア会第16期総会特別講演

3月

国内

対策提案
スーパーコンピュータ「富サロゲートモデルを併用した多目的最適化フレームワ

坪倉 誠

ーク

岳」成果創出加速プログラ

2022 年

ムシンポジウム 富岳百景

3月

国内

技術交流会

「富岳」が拓くSociety 5.0
時代のものづくり

坪倉 誠「富岳」成果創出加速プロ

2022 年

グラム公開シンポジウム

3月

サロゲートモデルを併用し

濱田 玄太郎, 犬井

第 35 回数値流体力学シンポ

2021 年た自動車空力性能の 多目

太 一 , 清水 圭 吾 ,ジウム、オンライン。

12 月

的最適化シミュレーション

平岡 武宜, 大山

Takuji

Conference

2022 年NAKASHIMA,

Aerodynamik”,

TaichiGermany (online).

国内
国内

聖, 中島 卓司, 中
村 優佑
Multi-Objective
Aerodynamic

Shape

Optimization of Simplified

INUI,- 48 -“FahrzeugMunich,

1月

国外Automobile Using HPC

Yusuke
NAKAMURA,
Keigo

Shimizu,

Takenori Hiraoka,
HiroakiFukumoto, Akira
Oyama, Takahide
Nouzawa, Makoto
Tsubokura
第 35 回数値流体力学シン

2021 年

ポジウム

12 月

第 35 回数値流体力学シン

2021 年

ポジウム

12 月

第 35 回数値流体力学シン

2021 年

ポジウム

12 月

第 35 回数値流体力学シン

2021 年

ポジウム

12 月

第 35 回数値流体力学シン

2021 年

ポジウム

12 月

Keynote Speech, ISHVAC

2021 年

2021, Seoul, South Korea

11 月

Hanyu Li, Kazuki

IEICES

2021 年

andKuga,

NguyenKyushu University

the

Dang

Khoa,

台風下の実市街地における

川口真晴

強風・風圧のピーク 特性に国内関する LES 解析（口頭発表）
強風時における広域市街地

河合英徳の LES に基づく流入 変動

国内

風データベースの作成（口頭発表）
WRF-LES を用いた 2019 年

中島慶梧

台風 19 号時の都市の 地表

国内

近傍流れのシミュレーショ
ン（口頭発表）2018 年猛暑日における東京

新井舞子

都市域での弱風・温熱場予国内測に関する研究（口頭発表）
洪水時の建物室内への浸水

西口浩司

評価に関する BCM の適用

国内（口頭発表）
Integrated CFD and In
SilicoModeling

Kazuhide Itofor

Airborne

Transmission

Analysisin

国外

IndoorEnvironment (口頭)
Effectsof

InitialConditions
Parameterson

Prediction of SARS-CoV-2

Kazuhide Ito

Viral Load in the Upper
Respiratory Tract Basedon Host-Cell Dynamics (口- 49 -

2021,

IGSES,

12 月

国内

頭)

In Silico Human Model Kazuhide Ito

2021

International 2021 年 国外

based Indoor Design forConference

Healthy Environment ( 口Building and Low Carbon

頭)Technology, China, 2021,
Xi'an

of

Green 11 月

UniversityofArchitectureandTechnology
Computersimulated

SCONNA

Kazuhide Ito2022,

2022 年

Oklahoma, USA, 2022

1月

第 2 回「富岳」流体予測革新

2022 年

プロジェクトシンポジウム

3月

第 5 回 CAE ワークショップ

2022 年の実現に貢献する燃焼数値～CAE シミュレーションの3月

シミュレーション（口頭）

最新動向～

personforpredicting

国外

inhalation exposure andairborne infection (口頭)
カーボンニュートラル社会

黒瀬良一の実現に向けた燃焼数値シ

国内

ミュレーションの役割（口頭）
カーボンニュートラル社会

黒瀬良一

国内

２．学会誌・雑誌等における論文掲載
掲載した論文（発表題目）

発表者氏名

Topology-free immersedboundary method forincompressible turbulence
flows: An aerodynamic
simulation for “dirty” CAD

immersed boundary
method using moving
least squares

発表し

国内・

（学会誌・雑誌等名）

た時期

外の別

2021年

国外

2021年

国外

2021年

国外

Computer methods in
Keiji Onishi,
MakotoTsubokurageometry
A one-sided direct forcing

発表した場所applied mechanics andengineering, 378, 113734.
(DOI:
10.1016/j.cma.2021.113734)

Rahul Bale,
Amneet Pal Singh
Bhalla, Boyce E.
Griffith, Makoto
Tsubokura

Journal of Computational
Physics, 440, 110359.
(DOI:
10.1016/j.jcp.2021.110359)

Eulerian finite volume

Tokimasa

International Journal forformulation using

Shimada, Koji

Numerical Methods in

Lagrangian marker

Nishiguchi, Rahul

Engineering

particles forBale, Shigenobu(DOI: 10.1002/nme.6896)- 50 -

incompressible fluid-

Okazawa, Makoto

structure interaction

Tubokura

problems
Building Cube Methodを利用した車体フレームのボ
クセルトポロジー最適化

和田有司, 嶋田宗将, 西口浩司, 大高雅史, 岡澤重信,
坪倉誠

日本計算工学会論文集,2021巻, 20210019.
(DOI:

2021年

国内

2021年

国内

2021年

国内

2022 年

国外

国外

10.11421/jsces.2021.20210
019)「富岳」による飛沫・エア

シミュレーション（日本シ

ロゾル飛散シミュレーショ

坪倉誠

ンと小規模飲食店における

ミュレーション学会誌）,
40(3).

感染リスク評価について
スパコン「富岳」が実現す
るSociety 5.0時代の自動車

宙舞（自動車技術会中部支坪倉誠部報）, 89号.

ものづくり
and

Kazuki

Kuga,

Indoorexhaled air re-inhalation

Pawel

Wargocki,

e13003rate

Kazuhide Ito(DOI: 10.1111/ina.13003)

Exhaled Aerosol and Jet

Hana

Aerosol and Air Quality 2022 年

Flow

MehrdadResearch,

Khamooshi,

210338Breathing
under

zone

transient

Air,

2022,

32,

conditions assessed with acomputer-simulated
personCharacterization

During Nasal Sneezing

Salati,

Jingliang

Dong,

Kazuhide

Ito,

David

Fletcher,Sara Vahaji, Kiao
Inthavong- 51 -

22(4),2022,(DOI:
10.4209/aaqr.210338)