SENAC Vol. 55, No. 1（2022. 1）― 15 ―［共同研究成果］

流体解析における時間発展計算の AI サロゲートモデルの追求――― 直感による究極的な単純化はどこまで有効か？
松岡 浩
技術士事務所 AI コンピューティングラボ
筆者らは、東北大学サイバーサイエンスセンターの共同研究公募制度により、平成 30 年度から
「リカレントニューラルネットワークによる高解像度流体解析コードの開発」を行ってきた。すでに、SENAC Vol. 54 No.1(2021.1)において、筆者がめざしている「“格子ガス法流体解析モデル”「流体計算をと“ニューラルネットワーク”の融合[[1]」について、その趣旨を述べた。すなわち、行う際に設定する格子点に“多層パーセプトロン”が内蔵されていると考え、格子ガス法に基づ
く流体計算の時間発展計算を微小整数だけを用いるニューロン計算で行い、その結果を疎視化す
ると、ナビエ・ストークス(NS)方程式を解いた結果とほぼ同等の結果を得る」という手法の可能性について述べた。この手法は、NS 方程式による時間発展計算をニューラルネットワークの計算で代替するという意味で、本稿では、“AI サロゲートモデル”と呼ぶことにする。
これからの“ものづくりをめぐる計算機シミュレーション”の究極的な目標は、マルチスケー
ル、マルチフィジクス、リアルワールド(不完全な計測情報でも機能するリアルタイムデータ同化
等)、カーボンニュートラル(温室効果ガス排出削減のための超低消費電力計算等)の同時実現であると思われる。“AI サロゲートモデル”は、この目標に近づく有力な１候補であると考えられる
が、まだまだ力不足である。ここでは、この“AI サロゲートモデル”を極力単純化することによって、計算機の設計や性能に対する要求条件を大きく緩和し、少しでも上記目標の達成を容易にする可能性を追求した。
なお、ここで述べる内容は、学術論文的なものの要点解説ではない。まだ、研究途上であるた
め、これまでの共同研究を通じて、筆者が考えるに至った“AI サロゲートモデル”の方向性を、読み物的に紹介するものとした。また、計算法の詳細は、過去の SENAC に掲載した内容[[1-5]とほ【補足】に記
とんどが共通であり、本稿では、AOBA-A で行った計算の代表的な計算条件を末尾の【[6] ”載するにとどめた。計算法の基礎理論は、Teixeira が考案した“４次元 FCHC54 速度モデル[に基づく“格子ガス法[[7-10]”である。筆者なりの理論解説は、文献[[11]に掲載している。
１．ものづくりをめぐる計算機シミュレーションへの要求の激化
ものづくりをめぐる計算機シミュレーションへの要求は、時代とともにハードルが急速に高まっている。
計算機性能の飛躍的な向上とあいまって、はじめに強調されたのは、「実験や試作を最小限の数にできる“高精度なマルチスケール and/or マルチフィジクスな事象解析”の実現」であった。
次に、Society 5.0 に代表される高度な IoT 時代の到来とともに“リアルワールド問題”への対応を可能にする“不完全な計測情報でも機能するリアルタイムデータ同化(計測データ等のシミュ
レーション計算への実時間反映)”も重要な課題になった。しかし、はじめの時点では、スーパーコンピュータによる大規模計算が必要な“高精度なマルチスケール and/or マルチフィジクスな事
象解析”と、計測情報を入力しながら超高速計算にフィードバックする“リアルタイムデータ同化”の同時実現は、技術的にハードルが高すぎる遠い目標であり、それぞれ別の達成目標とみなされていた。その後、特に東日本大震災を契機にして、津波等に対する防災対策強化などの観点から、大規模計算とリアルタイムデータ同化を両立させる方向の研究開発が盛んになったように思う。また、人工知能(主に、ニューラルネットワーク)の第３期ブームが到来し、深層学習(DeepLearning)等の多方面への応用が注目されるようになり、大規模計算の超高速化にも成果をあげつ― 16 ―

SENAC Vol. 55, No. （1 2022. 1）

つある。また、これらの動きに並行して、スーパーコンピュータの計算性能は、半導体の微細加工技術による性能向上が限界に近付きつつあるとは言え、着実に高速化を続けてきた。さらに、現時点では応用分野が限られているとはいえ、量子コンピュータの実運用も開始された。こうした背景から、すでに、“マルチスケール”、“マルチフィジクス”、“リアルワールド(不完全な計測情報でも機能するリアルタイムデータ同化等)”の同時実現は、単なる妄想ではなく現実味のある目標になってきている。
そして今、世界は、“カーボンニュートラル最優先の時代”を迎えた。ものづくりについては、製品の設計、試作から、その利用・廃棄・リサイクルに至るライフサイクル全体での温室効果ガ
スの排出量削減が要求される。ものづくりをめぐる設計・試作・製造・利用(運転)などのために活用する計算機シミュレーションに対しても、温室効果ガスの排出につながる火力発電の利用を抑制するため、“超低消費電力化”の要求が今まで以上に重視されることになろう。
以上は一般論であったが、ここで具体的な例をひとつとりあげてみる。
たまたま、電気新聞を見ていたところ、「2021 年 11 月 24 日、スペインのシーメンス・ガメサ・
リニューアブル・エナジーが、フィンランド最大の陸上風力発電所向けに、風車 69 基を受注した。」との記事があった。これは、ロータ直径が約 170ｍ、タワーの高さはおよそ 150ｍということである。風車のブレード各部における空力計算に必要な解像度を大きめに見積もって仮に cm オーダーであったとしても、ロータ直径とのスケールの違いは４桁ある。また、ウインドファーム
においては、ある風車の後流は、下流側の風車の存在やその風車の後流と重なり、さらに下流の風車に対しても複雑な後流となって到達する。精度の高い発電量の予測には、地形を考慮した風況予測に加えて、このような後流の状況をウインドファーム全体について解析する必要がある。
このため、少なくとも数 km～数 10km 以上のスケールまでカバーするマルチスケール解析が要求される。また、複雑に変化する後流の中に存在する風車の構造解析では、弾性変形をともなう強度解析や、寿命の長い安全運転を実現するための疲労解析を行う必要がある。すなわち、流体と複雑な構造物の間の連成を伴うマルチフィジクス解析が必要になる。さらに、いろいろな風況の出現頻度の予測ができたとしても、リアルワールドで生じる時々刻々の風況予測は不可能である。
結局、運転中における疲労損傷の蓄積が小さくなるような運転制御を実現するには、現地センタ
ーの計測情報と連携したリアルタイムシミュレーションが重要になってくる。そして、最後に、
カーボンニュートラルの実現のために導入されるウインドファームについて、その最適配置、強度設計、運転制御等に必要なコンピュータシミュレーションで大電力を消費してしまうのでは意
味がない。
すなわち、
“ものづくりをめぐる計算機シミュレーション”の今後の究極的な目標は、マルチス
ケール、マルチフィジクス、リアルワールド、カーボンニュートラルの同時実現である。
２．“AI サロゲートモデル”という方向性への期待“マルチスケール”、“マルチフィジクス”、“リアルワールド(リアルタイムデータ同化)”、“カ
ーボンニュートラル(超低消費電力)”という技術的なハードルは、いずれも極めて高い。スーパ
ーコンピュータの性能が著しく向上してきたとは言え、これらを同時に解決できるほどの実用的な方策は見つかっていない。現時点では、ハードウェアとソフトウェアに関する個々の技術を積み上げて、少しずつでも目標に近づく努力をするほかはない。言い換えれば、目標達成には、まだかなりの時間を要すると思われる。他方、それでも、“ものづくりをめぐる次世代の数値シミュ
レーション技術”として、何か抜本的に大きな効果のある対応策はないものだろうか？と思いをめぐらすのが人情である。
以上のような状況の中、“AI サロゲートモデル”という方向性は、かなり効果のある対応策であると考えられる。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 17 ―

一般に、“AI サロゲートモデル”という言葉からは、どのようなものを想像されるであろうか？
ただし、ここでは、冒頭で述べたように、本稿では、「NS 方程式による時間発展計算をニューラル
ネットワークの計算で代替するモデル」という狭い意味で、この言葉を用いている。
いまどきの AI は、ほとんどが“深層学習(DL: Deep Learning)”のイメージであろう。従って、
通常想像されるイメージは、例えば、以下のようなものではないだろうか。
【通常想像される“AI サロゲートモデル”のイメージ例】
マルチスケール解析を克服するため、まずは、非常に大規模な計算ができるというスーパー
コンピュータの“Capability Computing 能力”を活かし、NS 方程式による高解像度な時間発展計算を行う。これにより、１回の計算でも、「時刻ｔにおける詳細な流速空間分布」と「時刻ｔ＋⊿ｔにおける詳細な流速空間分布」のセットを非常に多数得ることができる。
次に、中規模な計算であれば非常に多数の計算を並列実行できるというスーパーコンピュー
タの“Capacity Computing 能力”を活かして、上記の計算で得られた各時刻の流速空間分布(速度ベクトルの各成分を表す３つの数値の巨大な３次元配列)を深層学習(DL)により次元圧縮する。次元圧縮によって、必要な記憶容量を大きく削減できる。ただし、次元圧縮で得られた数
値配列は、ニューラルネットワークの計算によって機械的に学習されたものであるので、その数値分布をみても物理的な意味は理解できない。
最後に、次元圧縮した時刻ｔと時刻ｔ＋⊿ｔにおける流速空間分布のセットを入出力データとしてニューラルネットワークに学習させる。このときの入出力データセットは非常に多数あるが、次元圧縮されているので、ひとつひとつの学習に要する計算規模は比較的小さい。
以上の結果、次元圧縮された流速空間分布の時間発展計算法が、ニューラルネットワークの計算法として学習されたことになる。これが、NS 方程式による時間発展計算を代替する“AI サロゲートモデル”のひとつのイメージである。
なお、この“AI サロゲートモデル”を用いて時間発展計算を行えば、直接得られるデータは、
やはり意味のわからない次元圧縮されたデータである。しかし、時間発展計算の過程で、我々がそのスナップショットを知りたい時はいつでも、次元圧縮の逆計算によってもとの高解像度レベルで物理的な意味のわかる流速空間分布を復元することができる。
以上のような技術は、深層学習(DL)の登場によってはじめて実用性が認識された AI 技術である。今後とも、その様々な活用によってこれまでにない重要な成果を期待できるであろう。
しかしながら、“マルチスケール”、“マルチフィジクス”、“リアルワールド(リアルタイムデー
タ同化)”、“カーボンニュートラル(超低消費電力)”の同時達成は、極めて高いハードルである。
この克服を目指すのであれば、何かさらに抜本的なアイデア、特に、上記の例で示した通常の“AIサロゲートモデル”のイメージをはるかに超えた“究極的な単純化”が必要だと思われる。
そこで、筆者は、以下、“AI サロゲートモデル”を究極的に単純化した姿について考察した。
３．“AI サロゲートモデル”の究極的単純化
今回の考察では、以下の３つの究極的な単純化について検討した。
（１）異なる物理現象の表現を単純化できないか？⇒“サブマクロ粒子”の挙動で代替
ものづくりをめぐる現実の流体解析の問題では、純粋な流体解析だけが要求される場合はほとんどない。構造物との熱流体力学的な相互作用をはじめ、相変化や、燃焼などの化学反応も同時並列して計算するべき問題も多い。これらの解析では、異なった物理現象が相互に関係しあって時間発展していくものと認識され、それぞれの物理現象に対して成立する支配方程式や経験式を

― 18 ―

SENAC Vol. 55, No. （1 2022. 1）

同時に満足するように連成計算がなされる。
他方、「すべての物質は、原子・分子からできており、異なるように見える物理現象も、結局は、“ミクロ粒子”の挙動で説明できるはず。」という自然科学的な知見がある。これは、マルチフィジクス解析をシングルフィジクス解析で処理するという単純化の方向を連想させる。ただし、「原子・分子スケールの計算から積み上げて、ものづくりをめぐる通常の工学スケールでの解析結果を得るという方法は、計算すべきマルチスケールの幅を極端に高めてしまう。」という重大な欠点がある。このため、ここでは、“ミクロ粒子(＝原子・分子)”ではなく、仮想的な“サブマクロ粒子”の挙動による物理現象の代替を考える。このことが、“単純化されたサロゲートモデル”を追求するうえで、一番重要な最初の選択になると考える。
【注】ものづくりをめぐる流体解析を行う場合、最終的に求めたい解像度のスケールにおける速度分布は、通常、“マクロな流体塊”の重心速度分布である。例えば、格子ガス法の場合、仮想粒子がもつ運動量などの物理量を、疎視化するために設定したある時空間領域で合計(平均)して時空間のスケールを一段上げることにより、“マクロな流体塊”がもつ物理量を求める。ここでは、このように、「一段スケールアップ(疎視化)しただけでマクロな流体挙動が得られるようなスケールで仮想する粒子」という意味を強調するため、“サブマクロ粒子”と呼ぶことにした。
また、“サブマクロ粒子モデル”の採用は、マルチフィジクスをシングルフィジクス化できるというメリット以外に、以下のようなメリットがある。すなわち、もともと流体解析だけというシングルフィジクスの場合を考えても、流体が激しく動きまわり、界面をもって砕け散るような事象については、通常の CFD では計算処理が困難である。このような事象に対しても、“サブマクロ粒子”のモデルは、一般的に対応しやすいとされている。
なお、現在の“サブマクロ粒子”のモデルを概観すると、一番注目されている手法は、「粒子法」であると筆者は思う。実際、従来手法では解析が困難な様々なマルチフィジクス事象について粒子法による解析事例が多数発表されている。また、「格子ボルツマン法」もかなりの成果を積み上げてきており、「格子ガス法」を呑み込んだ感じもある。以上は、筆者の個人的感想であるが、「粒子法」、「格子ボルツマン法」、「格子ガス法」は、「一定規模の原子や分子の粒子集団を仮想的な“サブマクロ粒子”に対応させてその挙動解析を行う」という点は共通であり、多様な
マルチフィジクスへ応用のアイデアは、相互に活用することができるであろう。
（２）サブマクロ粒子の状態記述を単純化できないか？⇒とりうる状態の数を微小整数に抑制“サブマクロ粒子”は、もともとが“多数の原子や分子の集団”である。このため、それがもつ物理量としては、いろいろなものを想定できる。具体的には、サブマクロ粒子内におけるミクロ粒子の粒子数密度、質量、運動量、角運動量、エネルギー、電荷などの分布やその合計値などである。また、それぞれの物理量がとりうる数値は、一般には、巨大な整数値や実数値になるので、とりうる状態の数は、ほとんど無限である。
しかし、ここでは、究極的な単純化を行う場合を考察しているので、まずは、サブマクロ粒子の挙動を特徴づける物理量を、“運動量”と“エネルギー”に限定してみる。
【注 1】仮想粒子に“質量”を担わせないモデルは、マルチフィジクスへの応用可能性を拡げる
可能性がある。例えば、弾性体内部の応力については、質量を伴わない“運動量の流れ(テンソル)”と解釈できる。また、運動量を担う同じ仮想粒子が、流体と弾性体の境界で相互に運動量を交換できるので、流体と構造物の間のより自然なシミュレーションを実行できると思われる。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 19 ―【注 2】仮想粒子に“質量”を担わせないモデルは、仮想粒子が存在しない状態「０」について、新しい解釈を可能にする。すなわち、「０」は、仮想粒子が存在しない状態を直ちに意味するのではなく、単に、仮想粒子が静止している状態も含む。この静止状態の物理的な解釈は、「多数のミクロ粒子の集団が滑らかに移動できずに微小振動しながら渋滞している状態」のイメージになる。
さて、前述した３つ手法(粒子法、格子ボルツマン法、格子ガス法)は、“サブマクロ粒子”の挙動計算方法について、それぞれ異なる特徴をもつ。ここでは、究極的な単純化という意味で、浮動小数点表示による実数演算が不要であるというメリットを重視し、“格子ガス法[[7-10]”に基づく“サブマクロ粒子”のモデルを選択した。
第５章に後述する試計算では、Teixeira が考案した“４次元 FCHC54 速度モデル[[6]”を使用している。このモデルでは、54 種類の“サブマクロ粒子(仮想粒子)”を設定する。その内訳は、エネルギーが０で静止状態の仮想粒子が６種類、エネルギーが１で運動量の向きが異なる仮想粒子が 24 種類、エネルギーが２で運動量の向きが異なる仮想粒子が 24 種類である。また、「ひとつの格子点には、同じ種類のサブマクロ粒子は、高々ひとつまでしか存在できない」という格子ガス法に特有な通常の仮定をする。この仮定によって、時間発展計算時の演算が著しく単純化され、記憶容量も極めて小さく抑えられる。
ここで重要なことは、「格子ガス法」の選択によって、サブマクロ粒子がとりうる状態の数を「54」という微小な整数に抑えるという単純化が可能になることである。
【注】「54」という数値は、通常の格子ガス法の応用としては、数値が大きいように思われるかもしれない。しかし、“４次元 FCHC54 速度モデル”は、格子ガス法が正しく機能するための条件(ガリレイ変換に対する不変性、応力テンソルの等方性、非圧縮性流体に関する NS 方程式への漸近性などの条件)を自然な形で満足できることが、Teixeira の論文[[6]によって示されている。この点を重視して、筆者は、専ら 54 速度モデルを利用している。
（３）近傍粒子の状態考慮を単純化できないか？⇒考慮範囲を最小近傍に限定し方向ごとに処理
オリジナルな格子ガス法では、各格子点にいろいろな向きから到着した仮想粒子は、格子点上で瞬時に衝突散乱を起こし、いろいろな向きに出発していく。このとき、到着粒子と出発粒子の間で、粒子数、運動量およびエネルギーの総和が保存される。すなわち、出発粒子が出ていく向きのパターンは、到着粒子の向きのパターンのみが考慮された衝突散乱規則によって決定される。
ただし、粒子数や運動量の保存則を満たす出発粒子のパターンが複数存在した場合には、その選択は等確率でなされる。
他方、格子ガス法の幅広い応用に向けて、相分離[[3]や相変化などを伴う特別なシミュレーショ
ンを行う場合がある。この応用では、最隣接格子点までか、場合によっては、格子点間距離の数倍遠方までの範囲で、仮想粒子間に働く力を考慮して出発粒子のパターンを決定するモデルが採用されている。なお、粒子法の場合は、“影響半径”の概念がこれに相当し、それぞれの粒子から“影響半径”以内の距離に存在する粒子からの影響を考慮している。このように、近傍粒子の状態考慮は、幅広い応用の実現にあたっては不可欠な選択となっている。
そこで、ここでは、衝突散乱直後の出発粒子は、近傍に存在する仮想粒子の影響を受けて最終的にたどり着く格子点が変化しうるモデルとした。
[6]”では、4 次元面心超立方体(FCHC)格子を
Teixeira が考案した“4 次元 FCHC54 速度モデル[

― 20 ―

SENAC Vol. 55, No. （1 2022. 1）

用い、各格子点上に存在できる仮想粒子の数は、前述のとおり、最大 54 個までである。仮想粒子が持ちうる速度を４次元方向の成分で表現すると、静止粒子(0,0,0,0)が６個までと、運動粒子が、エネルギー１の粒子：(±1,±1,0,0),(±1,0,±1,0),(±1,0,0,±1),(0,±1,±1,0),(0,±1,0,±1),(0,0,±1,±1)で 24 個、エネルギー２の粒子：(±1,±1,±1,±1),(±2,0,0,0),(0,±2,0,0),(0,0,±2,0),(0,0,0,±2)で 24 個となる。従って、ある格子点から１回の時間ステップで到達することができる１番近い格子点は、４次元空間において、エネルギー１の粒子に関して 24
個、エネルギー２の粒子に関しても同じく 24 個である。これら 48 個の最隣接格子点の状態は、それぞれの最隣接格子点に存在しうる 54 個の仮想粒子の存否情報(存在するか否かの情報)で表現できる。従って、ある格子点の近傍状態を求めるには、48×54＝2592 個の仮想粒子の存否情報を扱う必要があり、かなり複雑になることが予想される。
ここでは、究極的な単純化を目指す観点から、考慮する近傍の仮想粒子の範囲を極力せばめてみた。具体的には、最隣接格子点より手前で、格子点間距離の中央位置までを考慮すべき“最小近傍”とした。格子ガス法では、静止粒子以外のすべての仮想粒子は、時刻ｔから時刻ｔ＋⊿ｔの間に格子点間を同期して移動する。ということは、時刻ｔ＋⊿ｔ/2 の時点では、すべての運動
粒子は格子点間の中央付近を移動中と考えることができる。いま、時刻ｔ－⊿ｔ/2 の時点において、上記の“最小近傍”の範囲に存在しうる仮想粒子を考える。それは、①注目している格子点を時刻ｔ－⊿ｔに出発して離れていく仮想粒子か、あるいは、②注目している格子点に時刻ｔに到着する予定で近づきつつある仮想粒子、のいずれかである。従って、存否情報を扱う必要がある仮想粒子の数は、48 個の向きがあっても各向きに２個ずつなので、合計で 48×2＝96 個に激減させることができる。
また、ここでは、さらなる単純化を行う。具体的には、96 個の近傍粒子の存否情報を複合して扱うのではなく、１方向(互いに反対の２つの向き)ごとに、その方向の“最小近傍”を移動中の仮想粒子だけを観察する。この場合、１方向の“最小近傍”には、最大でも４個の仮想粒子しか存在しない。これら４個の最小近傍粒子の存否パターンがある特定の条件を満足した場合に、注目している方向のうち、ある一方の向きに飛び出そうとしている仮想粒子があれば、その速度を反転させて反対向きに出発させるという操作を確率的に行う。ここでは、このモデルを“仮仮想粒子反転連行モデル”を呼ぶ。
４．“仮想粒子反転連行モデル”の概要
液体の密度は、固体の密度とそれほど変わらない。従って、液体中では、固体中と同様に、粒子はいたるところに密集して存在する。ある粒子の近傍には、その粒子を取り囲むように粒子が存在するので、その粒子の自由な行動は許されない。このイメージをヒントにして、以下のように考える。
個々の格子点における到着粒子の衝突散乱によって、たまたまある向きに飛び出そうとする出発粒子があったとしても、その前後には仮想粒子、すなわち、ミクロ粒子の巨大集団が存在する。
前方に存在する仮想粒子の速さが遅ければ、出発粒子は、前方の仮想粒子に追突し、跳ね返されて減速される。後方にそれより速い仮想粒子が存在すれば、今度は、後方の仮想粒子に追突されて加速される。結果的に、出発粒子は、その前後に存在する仮想粒子による挟まれ具合いに応じて速さを変化させつつ、その流れの中に“連行”されてしまうであろう。
【注 1】上記のような連行過程のイメージにおいては、一般的に、「“速度の向き”が変化する揺動」も考えられる。しかし、ここでは、前述のとおり、考慮する仮想粒子の運動方向を１方向に限定している。また、“速度の向き”は、時間平均あるいは集合平均では変化しないことから、個々の“連行”においても、出発粒子の速度が他の方向に変化する揺動はないものと仮定した。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 21 ―

しかし、上述した「連行過程における個別粒子の“減速”や“加速”のイメージ」をそのままモデル化することはできない。なぜなら、１方向に限定した場合、運動粒子が持ちうる速さの選択
は、運動エネルギーが１または２の二者択一であり、個別粒子に微妙な“減速”や“加速”を行わせることはできないからである。そこで、“仮想粒子”が「もともとは粒子集団である」ことを思い出し、集団としての“速さ”の微妙な変化は、「出発粒子がもつ速度をそのまま保持しながら同じ向きに連行される確率」と、「出発粒子の速度が前後の仮想粒子との相互作用によって反転させられ、反対向きに同じ速さで連行される確率」の比率を用いて表現することにする。
【注 2】「仮想粒子の完全なる速度反転が生じる」という本モデルの仮定について、「多数の粒子が同時に反転するという事象が実際には起こりにくい」と思われることから、不適切な仮定ではないか？という疑問がわくかもしれない。しかし、一般に、サブマクロ粒子(仮想粒子)の挙動に関する仮定は、それを疎視化する時空間の範囲で合計(平均)したときに、マクロな流体塊の重心挙動を適切に表現できれば問題ないとされるべきである。よく用いられる“バウンズバック”の境界条件は、この代表的な事例であろう。
ただし、個々の“仮想粒子反転連行操作”は、仮想粒子の運動量を変化させる。特に何かの外力を模擬しようとしている場合でなければ、疎視化時空間の範囲で運動量を合計(平均)したときに、総和の保存が必要である。そこで、今回のモデルでは、ある向きへの連行と逆向きへの反転連行を、必ず交互に行うこととした。これより疎視化スケールでは、運動量保存則が成立している。
５ ． 試 計算に よ る 単純化 の 効 果に関 す る 考察（ １ ） １方向 ご と の最小 近 傍 仮想粒 子 の 挙動の 考 慮 だけで 何 が できる か ？
「サブマクロスケールにおける仮想粒子の挙動に対する単純な操作が、疎視化によって得られるマクロスケールの流体挙動において、どのような変化を発現させるか？」を定性的に観察するため、以下の試計算を行った。
はじめに、比較の基準を作るため、各格子点における衝突散乱の後、ある向きに飛びだした出発粒子が、「最小近傍粒子の影響を受けずにそのまま出発した場合」の円柱後流の様子を計算した。
図１に示す。
その結果を図ここでは、静止流体中に置かれていた無限長円柱に対して、時刻ゼロで、円柱軸に垂直な方向に
急に流体が動き出した場合の過渡変化を計算している。円柱後流の渦領域は、長い尾の形状になって一時的に円柱直径の約７倍の長さまで伸びたが、その後ゆっくり縮んで最終的には、約５倍の長さで安定した。
次に、「１方向ごとに最小近傍粒子の存否状態
を調べ、それに応じて、出発粒子の飛び出す向き
が前方か後方かの確率を意図的に偏らせた場合」に生じる円柱後流の様子を比較した。
図 １ ．近 傍 粒 子の影 響 を 考慮し な い 通常の 格 子 ガス 法 に よる円 柱 後 流の変 化
１方向の最小近傍粒子は、注目している格子点から見ると、“ある向き”と“その逆向き”の両

― 22 ―

SENAC Vol. 55, No. （1 2022. 1）

方に、それぞれ到着直前の粒子と出発直後の粒子が存在しうる。従って、最大４個の仮想粒子が存在しうる。最小近傍粒子の存否パターンは、1 方向あたり 16 パターンということになる。
このうち、４パターンは、格子点から見て、前後に存在する最小近傍粒子の配置が全く同じになる。この場合、最小近傍粒子の配置パターンを考慮しても、前後に飛び出す仮想粒子の出発確
率の間に意図的な偏りを生じさせる合理的な理由はない。そこで、この４パターンを除いた 12 パ
ターンすべてについて、その配置パターンが生じたときに、その直後の格子点における衝突散乱において「前方に飛び出す粒子がなく後方に飛び出そうとする粒子があれば、それを５％の確率で反転させる」という操作を行ってみた。
図２に示す。12 パターンのうち、３パターンでは、すぐにカルマン渦が発生した。
この結果を図また、長い尾が２パターン、短い尾が１パターン、非常に短い尾が１パターン、残りの５パター
ンでは、レイノルズ数が非常に小さい場合に生じる“円柱前後の対称流れ”が直ちに現れた。
これらの結果を分析すると、「１方向ごとに存在しうる４個の最小近傍粒子の合計の運動量が、ある向きに大きく、かつ、それに続けて同じ向きに仮想粒子を反転加速する確率が高いほど、す
なわち、周辺の粒子の流れに“連行される”と想像できる３パターンの場合に、円柱後流は、長い尾からカルマン渦に移行している」ということがわかる。これは、レイノルズ数が増加したときの変化である。試計算では、円柱の寸法は不変であり、入口のマクロ流速と各部の粒子数密度もほぼ一定に保たれているので、流体の動粘性が小さくなったものと解釈できる。従って、“１方向ごとの最小近傍仮想粒子の挙動を考慮する超単純なモデル”であっても、流体粘性を変化させる効果を発現でき、これを流体粘性の制御法として活用できる可能性がある。
【注】ここで、興味深い点がある。通常、格子ガス法で流体粘性を下げるためには、出発粒子を運動量勾配が正になる向きにできるだけ出発させるという方策がとられる[[13]。ところが、上記のモデルの場合、仮想粒子を反転させることによって出発確率を増す向きが、格子点前後における運動量の勾配がゼロ(下図の第４行第３列)あるいは負の向き(下図の第４行第５列)であっても、カルマン渦が発生し、流体粘性の低下が観察されている。
右図の第１行目は、前方の最小近傍に存在しうる到
着直前粒子と出発直後粒子の配置を示し、第１列目は、後方の最小近傍に存在する
到着直前粒子と出発直後粒子の配置を示している。
これら前後の配置パターンが同時に成立したときに、衝突散乱後の仮想粒子が出発する向きを、意図的に、前方を向く確率を増加させた場合の円柱後流の様子を示している。
なお、シミュレーションの計算体系等の詳細を【補足】
に示す。

前方
後方
カルマン
渦へ発展
円柱前後
対称流れ

長い尾

カルマン
渦へ発展

非常に
短い尾

長い尾

円柱前後
対称流れ

短い尾

円柱前後
対称流れ

円柱前後
対称流れ
カルマン
渦へ発展

円柱前後
対称流れ

図２ ． ⊿t / 2 過去 の 最 小近傍 仮 想 粒子の 配 置 に応じ て 衝 突散乱 直 後 の仮想 粒 子
を 前 方へ反 転 加 速する 場 合 の条件 と そ の結果 生 じ る円柱 後 流 の変化

流体解析における時間発展計算の
AI サロゲートモデルの追求― 23 ―

なお、上記の渦発生過程は、サブマクロスケールで以下のように解釈できる。
「各格子点の最小近傍には、合計 96 個の仮想粒子が存在しうる。このうち、48 個の到着粒子は、別の格子点から時々刻々やってくるので、これらを含めた 96 個の仮想粒子の重心速度は、激しく揺らいでいるであろう。従って、この１格子点ごとの重心速度(＝最大 96 個の仮想粒子小集団の重心速度)は、その格子点を含むマクロな疎視化領域(流体塊)における平均値、すなわち、“マクロな流速”を中心にしてそれなりに揺らいでいると考えられる。
１格子点ごとの重心速度が、たまたまある向きに揺らいだ時に、そのタイミングに合わせて続
けて同じ向きに仮想粒子が出発するような仮想粒子の反転連行操作を確率的にしかも前方と後方に向けて交互に行えば、この揺らぎは拡大を続け、やがて、疎視化時空間領域を超えたマクロスケールでの流体振動を発現するであろう。この流体振動は、より大きなマクロスケールの流体領域の重心を中心にいろいろな方向に発展していくので、そこに流体振動の位相差や流体領域の並進運動があれば、“目に見える渦”という形で視覚されるようになる。」
このように、「マクロな流体挙動は、サブマクロスケールで生じる仮想粒子の“連行”の強度、
すなわち、『格子点間隔スケールでの仮想粒子小集団の重心揺動が続けて同じ向きに生じる確率』に支配される。」ということができる。
45
また、渦発生過程とマクロな速度勾配の関係は、以下の通りに解釈できるであろう。
「一般に、マクロなスケールにおいて、速度勾配が一定の値に維持されている状態を考えてみる。マクロな流体塊の重心速度は、みな主流であるマクロな流速と同じ向きを向いている。このとき、マクロな46
流体塊の中をサブマクロなスケールで観察すると、各格子点の“最小近傍”に存在する仮想粒子小集団の重心速度は、主流と同じ向きではなくマクロな主
流ベクトルを中心にしていろいろな向きに揺らいでいると考えられる。この揺動状況を分析すると、速47
度勾配がある場所では、到着粒子の揺動主方向と出発粒子の揺動主方向は互いに乖離していることがわかる。この状態に、前述の“仮想粒子反転連行操作”を施してやると、到着粒子が来る頻度が少ない方向における“出発粒子の揺動”が拡大され、やがて、疎視化領域のスケールに達する。
48
このとき、主流速度の速い側から遅い側へ移動する流体塊は、主流に垂直でその流体塊の移動と同じ向きのマクロな流速速度成分をもつとともに、移動先の主流速度より常に少し速い主流方向のマクロな流速速度成分をもつであろう。逆に、主流速度の遅
49
い側から速い側へ移動する流体塊は、主流に垂直で
上記とは逆向きのマクロな流速速度成分をもつとともに、移動先の主流速度より常に遅い主流方向のマ
クロな流速速度成分をもつであろう。すなわち、主流に垂直なマクロな流体塊の往復振動は、その中心が主流に乗って下流に移動しながら楕円運動する流体の流れ、すなわち、渦の発生を意味している。そして、“仮想粒子反転連行操作”を行う確率を上げれば、速度勾配が生じる場所からの、このような渦の発生頻度を上げることができる。」

50

― 24 ―

試計算では、円柱の境界層における速度勾配からの渦発生が、“仮想粒子反転連行操作”の実行確率を上げることによって、より盛んになる様子が観測された。
図３に円柱境界層から発生して拡大していく渦の例(スナップショットの時刻ステップ 45～56)を示す。
なお、このような円柱から発生する渦が、このあと、ローラ渦や縦渦をうみ、複雑に重なった構造をもつ後流へと変化していくことが知られている[[12]。

SENAC Vol. 55, No. （1 2022. 1）

5152（ ２ ） どこま で 乱 流を模 擬 で きるか ？
乱流のように激しく変化する流れを計算機で安定的に計算していく観点からは、格子ガス法には重要なメリットがある。格子ガス法の時間発展計算は、ビット演算または微小整数演算のみで行うことができ、時間発展計算の個々の過程では、浮動少数点演算に起因する誤差の発生を避けることができる。このため、どんなに激しい流れについても、途中で計算不能になることはなく、必ず、何かの計算結果を安定的に求めることができる。
また、乱流状態は、大きなレイノルズ数の条件下で生じる。そして、多くの工学的なものづくり解析において、非常に大きなレイノルズ数条件下での解析が要求される。ここでは、どんな大きなレイノルズ数の解
析も理論上は実現可能なように、まず、流体粘性を負にできるシミュレーション方法を探求した。“仮想粒子反転連行モデル”の適用確率を調整することによって負の流体粘性を発現できることが、試計算によって定性的に確認できている。なお、格子ガス法による“負の粘性発現”を乱流解析に利用するアイデアについては、Rothman の論文[[13]を参照されたい。

53

54

55

56

図 ３ ． 円柱境 界 層 からの 渦 の 生成
さらに、実際にシミュレーション計算を実行する場合には、“負の流体粘性を発現させる方法”と“正の流体粘性を発現させる方法”を組み合わせて、模擬したい乱流状態に対応した小さい正の流体粘性を安定的に維持する具体的な方法が必要である。
図４に示す。
この実現を目指した試計算の例を図図２
具体的には、「衝突散乱直後の仮想粒子がその前後の流れにのみこまれる“連行参入”」が図の第２行第３列に示す近傍配置の場合に生じ、「前後の流れからはじき出される“連行離脱”」
図２の第３行第４列に示す配置の場合に生じるモデルを仮定した。そして、“連行参入”と“連が図行離脱”の確率をいろいろ変化させた。図４から、流体粘性が大きく変化している様子がわかる。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 25 ―

ケース 1
連行参入確率：20.0%
連行離脱確率：10.0％

ケース 4
連行参入確率：28.3%
連行離脱確率：1.7％

ケース 2
連行参入確率：25.0%
連行離脱確率：5.0％

ケース 5
連行参入確率：28.4%
連行離脱確率：1.6％

ケース 3
連行参入確率：28.0%
連行離脱確率：2.0％

ケース 6
連行参入確率：28.6%
連行離脱確率：1.4%

図 ４ ． 負の粘 性 発 現と正 の 粘 性発現 の つ り合い に よ る流体 粘 性 制御の 試 み（衝突散乱後の仮想粒子の“連行参入”が負の粘性を、“連行離脱”が正の粘性を発現させる。）

― 26 ―

SENAC Vol. 55, No. （1 2022. 1）

６．“仮想粒子反転連行モデル”がもつ潜在的なメリットとまとめ（１）理論的な観点
乱流の模擬を実現するにあたって、ナビエ・ストークス方程式を解く通常の数値流体力学(CFD)
による解法との関係が気になるところである。
まず言えることは、Teixeira の“４次元 FCHC54 速度モデル”による格子法流体解析の結果は、
非圧縮性流体に関するナビエ・ストークス方程式に基づく解析結果と比較して、マッハ数に関す[6] 。
る３次の精度で一致することが、理論的に示されている[
さらに、以下のとおり、CFD で必要になるような“複雑な乱流モデル”を考える必要はないと思われる。
“サブマクロ粒子”は、多数のミクロ粒子の集団を仮想的にひとつの巨大な粒子のイメージにまとめたものである。従って、実際には、あるサブマクロ粒子に属するミクロ粒子は、近傍のサブマクロ粒子に属するミクロ粒子と相互作用を起こすので、一般に、サブマクロ粒子がもつ運動量などの物理量は揺らいでいるはずである。原子や分子を想定した粒子に関する物理法則は、通常、このような内部構造に起因した物理量の揺らぎを仮定しないミクロ粒子を適用対象にして設定されたものと考えられる。従って、サブマクロ粒子の物理量についても同じ法則をそのまま適用していいとは限らない。より具体的には、以下のとおりである。
まず、原子や分子のミクロスケールで、物理学的な理論モデルを仮定する。ここでは、格子ガス法を考察しているので、この理論モデルは“ボルツマン方程式”であろう。計算機シミュレー
ションでは、離散化した格子上で計算を行わざるをえないので、実際には、“格子ボルツマン方
程式”をミクロスケールの超微細格子上で計算することが、物理学的理論モデルを忠実に解く方法だと考えられる。従って、本来は、ミクロスケールの超微細格子で“格子ボルツマン方程式”に基づく粒子挙動を計算し、これを合計(平均)してマクロスケールの流体挙動を得る必要がある。
しかし、このような計算は、常識的な時間では実行できない。そこで、格子ガス法では、サブマ
クロ粒子が“格子ボルツマン方程式”を満足すると仮定して、サブマクロスケールから計算を開始する。この結果が、ミクロスケールから時間をかけて計算したと仮定した場合の結果に、理論的に一致するか否か？が問題である。
Teixeira の“４次元 FCHC54 速度モデル”については、サブマクロスケールで平均された仮想粒子の数密度が、ミクロ粒子の数密度が満足すべきと仮定した“格子ボルツマン方程式”を精度よく満足することが示されている[6]。また、このことは、「サブマクロスケールにおいて個々の格子点の物理量が激しく揺らぎ大きなノイズが入っているように見えても、その数値をそのまま合計(平均)して疎視化すれば、必ず適切な結果を得ることができる」ことを意味している。
（２）AI 的な観点
本稿で掲げた具体的な主題は、「NS 方程式による時間発展計算を超高速に実行できる“AI サロ
ゲートモデル”を構築する」ことであった。このため、格子ガス法における仮想粒子のイメージ
から出発することにした。特に、格子ガス法における仮想粒子の衝突散乱と並進移動のイメージ
は、「“仮想粒子が移動するパスをシナプス結合にもつニューロン”が多数結合したニューラルネ
ットワーク」を自然に連想させる。具体的には、格子ガス法の時間発展過程は、リカレント型多[3]。このため、“誤差逆伝搬法”による学習がリアル
層パーセプトロンの時間展開計算に等しい[
タイムデータ同化を実現する基本的手法として利用できる。また、このニューラルネットワークの計算は、入出力と重みの値を微小整数の範囲に限定できるので、演算の簡素化と記憶容量の節約により、高速で低消費電力な計算の実現に大きく寄与できるであろう。さらに、１格子点に関する時間発展計算を１ビット幅で実行できるので超並列計算が可能であることも、従来から述べている[[2-4]。また、ニューラルネットワークの基本構造を、格子ガス法の物理的イメージから導流体解析における時間発展計算の
AI サロゲートモデルの追求― 27 ―

いているので、計測情報を学習して得た新しい確率頻度の物理的な意味を考察できる。すなわち、“説明可能な AI”になるので、次の改良に我々の思考を反映させやすくなる。
（３）将来展望のまとめ
今後の工学的なものづくりをめぐる計算機シミュレーションの将来を展望したとき、「マルチス
ケール、マルチフィジクス、リアルワールド(不完全な計測情報でも機能するリアルタイムデータ
同化等)、カーボンニュートラル(温室効果ガス排出削減のための超低消費電力等)の同時実現」が究極的な達成目標になってくると考えられる。これを実現する道のりは長いかもしれないが、本稿で述べたような“サブマクロ粒子イメージに基づく単純化された AI サロゲートモデル”は、その前進に大きく貢献できると思われる。特に、確率的なビット演算による仮想粒子の衝突散乱計
算は、“ゲート型量子計算機”の応用に適していると思われる。上記の究極的な目標は、結局、量子計算機によって達成されるのかもしれない。
しかし、まずは、足元の第一目標として、幅広いレイノルズ数の流体シミュレーションを可能にするため、“ゼロに近い正粘性を安定的に発現できる方法”が必要であり、この方法を可能な限り単純なモデルで確立し、しっかり検証しておくことが望まれる。ここでは、まだ不十分であるが、その方向のヒントになる“仮想粒子反転連行モデル”を示した。

謝辞
本稿で述べた研究課題の実行にあたっては、“AI サロゲートモデル”の候補となる“モデル”を多数考案し、それぞれのモデルについて、多数のパラメータスタディを行う必要があった。このとき、異なる計算モデルを試すたびに計算コードのチューニングに多くの時間をさくことは非効率的で悩ましい。この点、東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュ
ータ AOBA-A は、特別なチューニングをしないでも十分な計算速度を確保することができ、大変有益であった。また、利用にあたって同センター関係各位のご親切なご指導とご協力をいただき、心から感謝する次第である。今後とも、このような使い勝手のよいベクトル型スーパーコンピュ
ータのさらなる開発導入と同センターの有意義な活動の継続を期待している。
参考文献[1] 松岡,“格子ガス法流体解析モデルとニューラルネットワークの融合”, SENAC Vol.54
No.1, pp.39-49, 2021
[2] 松岡, 菊池,“多速さ格子ガス法実用化展開への手がかり”, SENAC Vol.49 No.4, pp.1-15,
2016[3] 松岡,菊池,“コンパクトな計算機によるリアルタイム流体解析の実現に向けて”, SENAC
Vol.51 No.2, pp.1-10, 2018[4] 松岡,菊池,“仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性”, SENAC
Vol.52 No.2, pp.18-27, 2019[5] 松岡,菊池,“リカレントニューラルネットワークによる実世界流れ場解析用時間発展計算モデルの探求”, SENAC Vol.53 No.1, pp.25-33, 2020
[6] Chistopher M. Teixeira, “Continuum Limit of Lattice Gas Fluid Dynamics”, MIT,1993
[7] B. Haaslacher, U. Frisch, Y. Pomeau, “Lattice Gas Automata for the Navier-Stokes
Equation”, Physical Review Letters Vol.56, No.14, pp.1505-1508, 1986
[8] Uriel Frisch, Dominique d'Humiè res, Brosl Hasslacher, Pierre Lallemand, Yves
Pomeau,Jean-Pierre Rivet, “Lattice Gas Hydrodynamics in Two and Three Dimensions”,
Complex Systems,１(1987), pp.649-707, 1987
[9] Michel Hénon, “Viscosity of a Lattice Gas”, Complex Systems, 1(1987), pp.763-789,

― 28 ―

SENAC Vol. 55, No. （1 2022. 1）

1987[10] Hudong Chen, Chris Teixeira, Kim Molvig, “Digital physics approach to computational
fluid dynamics: Some basic theoretical features”, International Journal of Modern
Physics C, Vol.8, No.4 (1997), pp.675-684, 1997[11] 松岡,“ビット演算による CFD と等価な高精度流体解析手法”, RIST News No.64, pp.1728, 2018[12]日野,“乱流の科学”, pp597-637, 2020, 朝倉書店（ISBN 978-4-254-20161-1）[13] Daniel H. Rothman, “Negative-Viscosity Lattice Gases”, Journal of Statistical
Physics, Vol.56, Nos. 3/4, 1989
円 柱 後 流 の過渡 変 化 計算に よ る 流体粘 性 の 変化の 観 察 に利用 し た 計算体 系 等【補足】円
粘性制御の効果を確認するため、円柱後流のシミュレーション計算を行った。粘性に変化があれば、後流の様子が変化するはずである。
具体的には、東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュータ AOBA-Aの 16CPU(128 コア)を利用し、MPI による並列計算で約 12 分の小規模計算で求められる円柱後流の３次元過渡変化計算（約 3900 万格子点で 25600 時刻ステップ、円柱軸に垂直なある平面上の流体挙動）を試みた。
計算条件等の詳細を下下図に示す。
３次元空間中で多数の格子点をＸ
ＹＺの各方向に並べ、直方体形状の格Z

子点配列を作る。各格子点は、その内

流体入口

流体出口
X

流れ

部に４次元目の座標として R＝R

Z

0,1,2,3 の位置を識別できる自由度をXもつとする。これが左図であり、４次
元面心超立方体格子を３次元空間へ

Y

Y

Z

投影した姿である。今回の数値シミュ
レーションでは、３次元縮退格子として、Ｘ方向に 1920 個、Ｙ方向に 1280
X個、Ｚ方向に 16 個の格子点を配置した。
格子点間隔の長さを 1 にとると、円柱の直径は 40 に設定した。

Y[ 過 渡変 化シ ミュ レ ーシ ョンの 条 件]
シミュレーション計算を開始する時刻ステップ０の時点で、各格子点には、そこに存在できる仮想粒子の最大数の 20％の数の仮想粒子をランダムな向きの速度で配置する。この結果、疎視化して得られるマクロな流速はゼロであり、流体は、直方体形状の中で静止している。
次に、時刻ステップ１の時点から、＋Ｘ向きの速度をもつ仮想粒子をＸ＝０の位置から注入していく。
すると、時刻ステップが進むにつれて、流体全体が＋Ｘ向きのマクロな速度をもつようになる。このとき、＋Ｘ側の先にある直方体出口においては、出口直前に存在する格子点上の仮想粒子配置を、出口直後に存在する格子点の仮想粒子配置にコピーして、出口におけるマクロな流速の勾配がゼロになるという境界条
件を近似的に実現した。また、±Ｙ方向と±Ｚ方向には、周期的境界条件を適用した。この流れの中の入り口に近い位置に、“Ｚ方向の中心軸をもつ無限大の長さの円柱”を置き、その後流に生じる流体挙動を計算した。

図 ． 円 柱後流 の 過 渡変化 シ ミ ュレー シ ョ ンを行 う 計 算