SENAC Vol. 54, No. 1（2021. 1）― 39 ―［共同研究成果］

格子ガス法流体解析モデルとニュー ラルネットワークの融合
松岡 浩
技術士事務所ＡＩコンピューティングラボ
筆者らは、東北大学サイバーサイエンスセンターの共同研究公募制度により、平成 30 年度から３年計画で「リカレントニューラルネットワークによる高解像度流体解析コードの開発」を行っている。本稿では、その２年目に当たる平成 31 年度(令和元年度)以降に行った研究のうち、“格子ガス法流体解析モデル”と“ニューラルネットワーク”の融合により実現を目指す３つの目標：
１．不完全な実世界情報でも利用可能なリアルタームデータ同化
２．どんなに激しい乱流でも安定的に計算可能な超並列整数計算
３．極めて単純な要素モデルからでも自己組織化可能な粘性制御について、その趣旨とこれまで試みた基本的な方法を概説する。
なお、説明の都合上、１年前の SENAC Vol. 53 No.1(2020.1)に掲載した研究成果「リカレント
ニューラルネットワークによる実世界流れ場解析用時間発展計算モデルの探求」と重複する内容の記載がある。これは、説明の飛躍をなくして読みやすくするためであり、ご容赦願いたい。
１．不完全な実世界情報でも利用可能なリアルタームデータ同化をめざして（１）C. M. Teixeira が考案した高精度な格子ガスモデルから出発する。
本研究では、C. M. Teixeira が考案した４次元面心超立方体格子による「３速さ格子ガス法 54
速度モデル[1]」を出発点にして、ニューラルネットワークの計算原理を導入する。ニューラルネ
ットワークは、入力信号の一部が欠落したり、間違っていたりしても、それなりに適切な答えを出力することができる。この特徴は、「計測システムを設置できない場所の物理場情報の欠落や、計測システムの故障等による間違った物理場情報の混在が一部にあったとしても、それなりに適切なシミュレーション結果を出力し続けることができる」という実世界応用への可能性を広げる。
また、Teixeira の 54 速度モデルに基づく格子流体解析の結果は、ナビエ・ストークス方程式を解く標準的な数値流体力学(非圧縮性流体)の結果と比較した場合、マッハ数に関する３次の精度まで一致することが Teixeira の論文[1]において証明されており（cf.関連する重要論文[2]）、
理論的にもそれなりの精度を期待できるものである。（cf.関連する解説[3,4]）
図 1 に示す。
なお、格子ガス法流体解析の一般的な概要を図

連続体流れ場（連続体の流動）

仮想粒子の格子点挙動（格子点上での衝突散乱と格子点間の並進移動）

格子ガス法流体解析では、流体が存在する空間中に規則正しく格子点を配置し、特定の質量と運動量を担う多数の仮想粒子が、①格子点上での「衝突散乱」、②格子点間の「並進移動」
を繰り返しながら移動していく様子を疎視化(平均)して、マクロな流体物理量(密度、運動量、速度等)の挙動を得る。

図１．格子ガス法流体解析の一般的な概要

― 40 ―

SENAC Vol. 54, No. （1 2021. 1）（２）格子ガス法の衝突散乱過程は多層パーセプトロンの計算に等しい。
格子ガス法流体解析では、流体が存在する空間(流体場)に一定間隔で多数の格子点を規則正しく配置する。次に、一定の質量を担う多数の“仮想粒子”が、これらの格子点上で互いに衝突散乱を繰り返しながら格子点間を次々と並進移動していく様子を計算で求める。そして、ある一定の時間が経過するごとに、近傍の格子点に存在する仮想粒子の質量や運動量を合計(疎視化)して、流体場各部におけるマクロな流体物理量(密度、流速等)の時間変化(スナップショット)を導く。
このとき、各格子点では、短い時間間隔⊿τが刻まれるたびに、近傍に存在するいくつかの格子点から仮想粒子が飛来する。これらの“到着粒子”は当該格子点上で他の到着粒子と衝突散乱を起こし、瞬時に飛行の向きを変えていろいろな近傍格子点に向かう“出発粒子”になる。このモデル上の仮定から、“到着粒子”も“出発粒子”も、近傍格子点間の移動をぴったり時間⊿τで行う必要があるため、仮想粒子がもちうる速度は、連続的ではなくいくつかの“離散ベクトル値”になる。Teixeira の 54 速度モデルでは、これらの離散ベクトル値を 54 個考え、ひとつの格子点には、それぞれが最大１個まで存在できるとする。従って、ひとつの格子点に存在できる仮想粒子の最大数は 54 個である。ただし、その内訳は、①６個の“静止粒子”、②速さが c で向きが異なる 24 個の“遅い運動粒子”、③速さが 2c で向きが異なる 24 個の“速い運動粒子”であり、静止粒子だけは例外で、同じ速度ゼロをもつ仮想粒子がひとつの格子点に複数存在してもよい。
ひとつの格子点における 54 種類の仮想粒子の存在状態は、54 ビットのビット列パターンで完全に表現できる。例えば、種類 i（i＝1,2,･･･,54）の仮想粒子が存在する場合、i 番目のビットを「1」にし、存在しない場合は i 番目のビットを｢0｣にすればよい。この表現を用いれば、各格子点で生じる仮想粒子の衝突散乱は、「到着粒子の存否を表す 54 ビット列パターン」を、「出発粒子の存否を表す 54 ビット列パターン」へ変換する関数機能に他ならない。一般に、このような関数機能は、“多層パーセプトロン”で実現できる。従って、格子ガス法による時間発展計算を実行する立場からは、「流体場に存在する各格子点は、仮想粒子の衝突散乱過程を計算する“多層パー
セプトロン”を内蔵している」と見なすことができる。
ここで、念のため、２次元格子ガス法 FHP モデル[5]における仮想粒子の衝突散乱過程を例に
図２に衝突散乱過とり、その過程を計算できる具体的な多層パーセプトロンを構成してみよう。図程の１例とその計算を実現できる多層パーセプトロンの構成例を示す。
31４青丸：格子点（-1+1+1-1+1-1）
1 -1
2
1
3 -1
4
1
5
1
-1
6
到着粒子の存否
入力ビット列2 次元格子ガス法 FHP モデルでは、正三角形格1子を用いる。各格子点の近傍格子点は６個あるが、

6
出発粒子の

上での衝突散乱の結果、1,3,6 の近傍格子点に向か

４

6
5
5
赤丸：仮想粒子
到着粒子のビットパターン2

3

2

ビットパターン（+1-1-1+1-1+1）

入力層

中間層 出力層

左図では 2,4,5 から仮想粒子が３個飛来し、格子点って３個の仮想粒子が出発している。衝突散乱の前後で粒子数と運動量が保存されている。

1
-1
1
-1
-1
1

1
2
3
4
5
6

多層パーセプトロン

出発粒子の存否

内蔵格子点

出力ビット列左上図の衝突散乱過程の計算
は、左図の３層パーセプトロ
ンで実現できる。この例では、仮想粒子の存否を｢1or0｣ではなく｢1or－1｣に対応させているので注意。計算内容は、本文中で説明を参照されたい。

図２．２次元格子ガス法 FHP モデルにおける衝突散乱過程の１例とその計算を実現できる多層パーセプトロンの例

格子ガス法流体解析モデルとニューラルネットワークの融合― 41 ―

図２に示す多層パーセプトロンは、３層からなる。“入力層”は到着粒子の存否を表す入力値(±１)をそのまま中間層に伝える。中間層は入力ビット列のパターンを識別する。例えば、中間層の一番上に位置する綠色ニューロンにおいて、赤線入力に+1、黄線入力に－1 の重みを設定したとすれば、上から黄(-1)→赤(+1)→黄(-1) →赤(+1) →赤(+1) →黄(-1)の重みパターンになる。
このとき、入力層への入力ビットパターンがこの重みパターンと一致するときにのみ積和が最大値(6)になる。この状態を、中間層ニューロンの活性化関数を“しきい値＝6 の階段関数”にする
ことによって識別する。そして、この中間層ニューロンの出力は、出力層のうち活性化すべきニューロン 1,3,6 にだけ接続しておく。出力層のニューロンは、接続されているシナプスの入力の重みをすべて+1 に設定し、活性化関数を“しきい値＝1 の階段関数”にする。このことで、いずれかひとつでも入力が+1 になれば活性化（出力＝+1）する。実際の適用時には、多数の衝突規則を実現する必要があり、それらが干渉しないように、このような階層を多段に重ねて適用する。
なお、衝突規則を確率的に適用する場合は、入力層のニューロンを増やし、そこに±１の入力を所定の確率でランダムに与えればよい。
（ ３ ） 格子ガ ス 法 の時間 発 展 過程は リ カ レント 型 多 層 パ ー セ プ ト ロ ン の 時 間 展 開 計 算 に 等 し い 。
格子ガス法流体解析における時間発展過程の全体もニューラルネットワークの計算で実現できることを示す。はじめに、仮想粒子の挙動をイメージしながら、①“格子点”⇔“ニューロン”、②“仮想粒子の存否”⇔“情報の内容(ここでは、±１とする。)”、③“仮想粒子の格子点間移動経路”⇔“シナプス結合”という対応を考える。これから、ただちに“相互結合型ニューラルネットワーク”の構造を連想できる。このとき、仮想粒子は、格子点間を相互に移動できるので、シナプス結合において、情報はニューロン間を双方向に流れうる。
次に、時間発展計算の手順を明示的に表現するため、上記の“相互結合型ニューラルネットワーク”を時間方向に展開する。このため、「時刻ステップｔにおけるニューロン」と「時刻ステッ
プｔ＋1 におけるニューロン」を別のものとして扱い階層的に並べる。そして、「時刻ステップｔ
におけるニューロン」の出力を「時刻ステップｔ＋1 におけるニューロン」の入力に伝えるシナプス結合を配置する。時間発展計算は、この繰り返しで実現できるため、「時刻ステップｔ＋1 におけるニューロン」の出力は、「時刻ステップｔにおけるニューロン」の入力にフィードバックされる。こうして、“リカレント型多層パーセプトロン”が得られる。このとき、各シナプス結合に
図３に示す。
おける情報の流れは、双方向ではなくひとつの向きになっている。イメージを図並進移動と時間推移を可能に

2 次元の場合する１時刻ステップ先の異なるのイメージ

格子点間をつなぐｼﾅﾌﾟｽ結合

時刻ステップｔ
衝突散乱

衝突散乱
到着粒
子ビッ
ト列

計算用多

出発粒

層ﾊﾟｰｾﾌﾟﾄ

子ビッ

ﾛﾝ

ト列

時刻ステップｔ+1
衝突散乱

並進
移動

到着粒子
スナップショット出力

到着粒子

出発粒子

図３．格子ガス法の時間発展過程を計算するリカレント型多層パーセプトロン

出 発 粒子

― 42 ―

SENAC Vol. 54, No. （1 2021. 1）（４）“誤差逆伝搬法”がリアルタイムデータ同化を実現するヒントになる。
“多層パーセプトロン”には、“誤差逆伝搬法”という強力な“教師あり学習法”が開発されている。実流体の過渡変化実測データまたは高精度なシミュレーションデータがあれば、例えば、ある場所のマクロな流体物理量の時系列データを教師データとして利用することができる。すな
わち、各格子点に内蔵された多層パーセプトロンによる衝突散乱過程の繰り返し計算によって得られたマクロな流体物理量の時間変化が、実測データまたは高精度なシミュレーションデータに合うように、多層パーセプトロンの重みを学習させることができる。学習後の重みが、例えば、
仮想粒子をある向きに加速するなど運動量保存則を一時的に破る操作に相当する場合もある。しかしながら、局所的な時空間平均では保存則を守るようにプログラムすることが可能である。そして、流体場のある場所の格子点に対するこのような操作の効果は、仮想粒子の衝突散乱過程を通じて短時間のうちに周辺の流体挙動に伝わり、物理的に矛盾のない新たな局所平衡状態をつく
りだす。これは、「リアルタイムデータ同化」を実現する有力な手法開発のヒントになる。
（５）全ての機能を含めたリカレント型多層パーセプトロンを構成する。
格子ガス法流体解析における“衝突散乱”と“並進移動”の計算は、すべての場所で全く同じというわけではない。例えば、流体中に不浸透な静止した固体壁が存在する場合、通常、“粘着条件”としてバウンズバック条件が課せられる。この場合、到着粒子の衝突散乱計算によって得られた出発粒子(＝到着粒子に比べて分布の凸凹が緩やかになった“緩和粒子”)の存否計算の結果は捨てられ、到着粒子の速度を反転させた存否分布を出発粒子の分布とする（下図の緑線）。
さらに、実世界の流速計測情報をリアルタイムでフィードバックするような場合、外部環境条
件に合うように個別粒子を巧みに加速するなど達人的な操作(“達人操作”)を行って、いわゆる“データ同化”(計測融合)を実現したい場合もある。また、後述する粘性制御のため、１時刻ス
テップ前の出発粒子の存否情報を用いる場合もある（下図の赤線）。
図４を得る。
これらをまとめて図
計測データ等の外部環境条件入力ス
並進移動と時間推移
１時刻ステップ未来に到着する格子点との間をつなぐｼﾅﾌﾟｽ結合

到着粒
子ビッ
ト列

衝突散乱

緩和粒

個別粒子へ

出発粒

の計算を

子ビッ

の達人操作

子ビッ

行う多層

ト列

の計算を行

ト列

ﾊﾟｰｾﾌﾟﾄﾛﾝ

う多層ﾊﾟｰｾﾌﾟﾄﾛﾝ

１時刻ステップ前の出発粒子の存否情報をリカレントする！
疎視化データ等のスナップショット出力

図４．全ての機能を含めたリカレント型多層パーセプトロン

格子ガス法流体解析モデルとニューラルネットワークの融合― 43 ―

２．どんなに激しい乱流でも安定的に計算可能な超並列整数計算を目指して（１）“整数型多層パーセプトロン”による計算は乱流計算に貢献できる。
これまで考察してきたことを振り返ると、格子ガス法流体解析の時間発展計算過程は、“多層パ
ーセプトロン”による計算の組合せで全てを実行できそうである。
図２に示した３層パーセプトロンの例を考察すると、
特に、図①入力層への入力パターンが重みパターンと一致するときにのみ積和が最大値になるので、中間層ニューロンの活性化関数のしきい値を当該最大値に設定して入力パターンを選別する。
②活性化関数を“しきい値＝1 の階段関数”にすることで、いずれかひとつでも入力が+1 にな
ればニューロンを活性化させて、出力＝+1 を実現する。
という２つの操作を組みあわせていることがわかる。これらの操作は他の“衝突散乱”の計算の場合にも広く適用できるものである。
また、“並進移動”については、シナプス結合のつなぎ方の問題であり、計算機上ではメモリ間の情報移動に過ぎないと考えても良いし、あるいは、格子点間にすべてのシナプス結合が存在していてどこの重みを“ゼロでなく+1”にするかの問題であると考えてもよい。
以上のことから、格子ガス法流体解析の時間発展計算の全過程は、各ニューロンの“入出力値”、“重み”、“階段型活性化関数のしきい値”をすべて微小な整数に設定した“整数型多層パーセプ
トロン”によって、微小整数の加減乗算のみで実行できることがわかる。
このことは、実数計算による時間発展計算を行う場合と比較して、計算の高速化と記憶容量の節約を可能にする。そして、流体計算の観点から最も重要なメリットは、実数計算の場合に生じる打切り誤差等の発生がなく、どんなに激しく変化する流れに対しても安定的に答えを出すことができる点にある。従って、乱流を解明する研究への応用が潜在的に期待される。

３．極めて単純な要素モデルからでも自己組織化可能な粘性制御をめざして（１）格子ガスモデルのスケールで大胆な粘性発現機構を考えることも意義があるかも。
格子ガス法は、それが考案された初期、等方性やガリレイ不変性などの基本的要件を満たさないことが指摘されたが、いずれも改良された手法で克服されている。現在、格子ガス法が数々のメリットを持ちながらも、実用的な流体解析手法として一般的でない唯一の理由は、発現できる
流体粘性を広い範囲で制御することが容易ではないからであろう。なお、格子ガス法の発展形態
としては、実数計算を用いる“格子ボルツマン法”が成果を挙げている。しかし、筆者としては、実数計算が不要であるメリットを活かしたい。そこで、整数演算のメリットを残したままで、極めて単純な要素還元モデルから広範囲の流体粘性を自己組織化的に発現できる手法を探求した。
流体粘性は、結局のところ、ファンデルワールス力等の複雑な分子レベルの力学的相互作用から発現する。従って、流体粘性を高精度に扱うためには、分子レベルのスケールをシミュレーシ
ョンする分子動力学を用いるのが自然なアプローチであろう。従って、分子レベルよりも桁違いに大きい格子ガスモデルのスケールで意図的な相互作用モデルを考案しても、そのモデルが発現する粘性効果は“大胆な近似”にしかならないとも思われる。
他方、ナビエ・ストークス方程式を解く王道的な“数値流体力学”は、マクロな連続体モデルであるにも拘わらず、粘性流体の挙動をかなり高精度にシミュレーションできる。また、格子ガ
スモデルと同様のスケールにおいて粒子間相互作用を設定する“粒子法”が開発されているが、いろいろな流体解析事例で成功を収めている。従って、格子ガスモデルのスケールで“大胆な近似”を導入しても、例えば統計平均的な効果により、マクロな疎視化スケールでは、現実の流体

― 44 ―

SENAC Vol. 54, No. （1 2021. 1）

挙動を高精度に模擬できる粘性発現機構になっている可能性があることも否定はできない。
（２）格子流体の粘性を自由に制御するには“等確率の常識”を超える必要がある。
「格子ガス法がマクロなスケールでどのような流体粘性を発現するか？」という問題は、より高いレイノルズ数領域の流体解析に格子ガス法を適用するため、「格子ガス法では、どこまで小さい値の流体粘性を発現できるか？」という点から重要である。
格子ガス法が提案された当初から、粘性に関する重要な論文がいくつも発表されている。例えば、Hénon は、「Viscosity of a Lattice Gas（1987）」[6]において格子ガス法によって発現される流体粘性の公式を導き、格子流体が発現する粘性は、格子点間隔を一定にした場合、衝突規則の詳細に依存して変化することを示した。また、（全ての衝突散乱が等確率で生じる場合）粘性は必ず正の値になることを明らかにしている。
格子ガス法の粘性制御に関する論文例としては、Chen, Teixeria, Molvig は、「Digital physics
approach to computational fluid dynamics: Some basic theoretical features (1997)」[7]において、通常の衝突散乱規則に“過緩和過程”を追加することで粘性をより小さくできることを示している。他方、Rothman は、「Negative-Viscosity Lattice Gases(1989)」[8]で負の粘性を発現させる方法を具体的に示すとともに、乱流解析への応用可能性を示唆している。いずれにしても、仮想粒子の衝突散乱過程を完全に等確率で行うのではなく、「到着粒子に特定の衝突規則だけを頻繁に適用したり、出発粒子を特定の向きだけに頻繁に放出させたり」というある種の“負のエントロピー発現操作”が提案内容である。すなわち、格子流体の粘性を自由に制御するには、“等確率の常識”を超える必要がある。
（３）幅広い流体粘性を発現させるため個別粒子への“達人操作”を導入してみる。
筆者らは、これまで“等確率の常識”を超えるいくつかの方法[9,10,11]を試してきた。例えば、１年前の SENAC Vol.53 No.1(2020.1)に掲載した研究成果[11]においては、「仮想粒子が並進する
際の経路は直線ではなくジグザグに揺動する経路であると仮定し、同じ速度をもつ個別粒子がときどき接近しすぎてジグザグ経路を乗り換える際に、前回の接近の際に粒子が存在した側のジグザク経路を頻度多く選択する」というモデルを示した。このような“等確率の常識”を超える方法に共通する目標は、「望まれる粘性をぴたりと発現できるように、仮想粒子の衝突散乱過程において個別粒子の運動に影響を与える何等かの巧みな“達人操作”を見つける」という点にある。
筆者らの研究ではまだ本目標を達成していないが、本稿では、上述の１年前の試みとはまた別の試みを示す。ただし、「流体粘性は、流体速度の時間相関に関わるものである」という非平衡系
統計力学の知識を参考にして、１時刻ステップ前の出発粒子の存否情報を“個別粒子に対して行
う達人操作”を行うタイミングの判断基準のひとつにする点は１年前の試みから変えていない。
“個別粒子への達人操作”は、“衝突散乱過程”のあと「ある向き(順向き)
今回試みた具体的な“に出発するはずだった仮想粒子(緩和粒子)を、ある条件が整ったときにだけ確率的に逆向きに出発させる」というものである。このときの確率は、すべての向きに対して同じ値を適用する。
ただし、この操作は、あるひとつの粒子を意図的に逆向きに加速したことを意味する。従って、局所的な時空間平均で運動量保存則が破れないように、今回の研究では「順向きの加速と逆向きの加速は必ず交互に行う」ものとした。
（４）粘性制御のための達人操作も多層パーセプトロンの計算で実現できることを確認する。
今回試みた粘性制御のための“個別粒子への達人操作”も、多層パーセプトロンの計算で実現
図５に示す。
できることを確認しておく。その１例を図

格子ガス法流体解析モデルとニューラルネットワークの融合― 45 ―

図５において、入力層は入力値±１をそのまま中間層に伝えるだけの役割を果たす。中間層は入力ビットパターンを識別する。赤線は+1、青線は－1 の重み設定を意味する。例えば、中間層の一番上に位置する黄色ニューロンにおいては、その７つの入力線は、上から赤(+1)→青(-1)→
赤(+1)→赤(+1) →青(-1) →赤(+1) →青(-1)の重みパターンになっている。この中間層の重みパ
ターンに一致するビットパターンが入力層にきた場合にのみ積和が最大値(７)になる。この状態
を中間層ニューロンの活性化関数を“しきい値＝7 の階段関数”にすることによって識別する。
このニューロンの出力は、出力層のうち活性化すべきニューロンである上のニューロンにだけ接続しておく。出力層のニューロンは、接続されているシナプスの入力重みをすべて+1 に設定し、活性化関数を“しきい値＝0 の階段関数”にする。このことで、いずれかひとつでも入力が+1 になれば活性化（出力＝+1）する。
図５の例では、中間層のニューロンは、上から順に、後述する「順向き局所主流巻込操
なお、図作」,「逆向き局所主流巻込操作」,「順向き近傍粒子引合操作」,「逆向き近傍粒子引合操作」を行うべき条件を識別するものである。

入力層
達人操作の適用確率で+1 の値をラ

中間層

出力層

≧7 で発火≧0 で発火≧7 で発火≧0 で発火ンダムに与える
順向き出発粒子

順向き緩和粒子の存否入力±１
逆向き緩和粒子の存否入力±１の存否出力±１

逆向き出発粒子の存否出力±１

順向き到着粒子の存否入力±１≧7 で発火逆向き到着粒子の存否入力±１≧7 で発火

１時刻ステップ前の順向き出発粒子の存否入力±１
１ 時 刻 ス テ ッ プ前の 逆 向 き 出 発 粒子の存否入力±１

図５．粘性制御を行う多層パーセプトロン（５）円柱後流の過渡変化計算により粘性変化の兆候を観察する。
上述した粘性制御の効果を確認するため、円柱後流のシミュレーション計算を行った。粘性に変化があれば、後流の様子が変化するはずである。
具体的には、東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュータ SX-ACE
を利用し、１CPU(４コア)による 40 分程度の小規模計算で求められる円柱後流の３次元過渡変化
計算（約 530 万格子点で 12800 時刻ステップ、円柱軸に垂直なある平面上の流体挙動）を試みた。
計算条件等の詳細を図図６に示す。

― 46 ―

SENAC Vol. 54, No. （1 2021. 1）Z

３次元空間中で多数の格子点をＸＹ

流体出口

流体入口

X

流れ

R

Z

Ｚの各方向に並べ、直方体形状の格子X

点配列を作る。各格子点は、その内部に４次元目の座標として R＝0,1,2,3

Yの位置を識別できる自由度をもつとYZする。これが左図であり、４次元面心
超立方体格子を３次元空間へ投影した姿である。今回の数値シミュレーシXョンでは、３次元縮退格子として、Ｘ方向に 864 個、Ｙ方向に 384 個、Ｚ方向に 16 個の格子点を配置した。

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

図６．円柱後流の過渡変化シミュレーションを行う計算体系
シミュレーション結果は、ある Z 断面上で疎視化された運動量分布を計算し、ParaView で過渡変化の動画とした。本稿では、この中からスナップショット画像を抜粋して表示する。
（６）流体粘性を変える“局所主流巻込操作(仮称)”の適用効果を観る。
“個別粒子への達人操作”として、「局所主流巻込操作（仮称）」なるものを考え、確率的に適用してみたので、以下、概要を示す。
マクロな速度がゼロで静止している流体が存在する場所にある格子点上では、移動することが可能なすべての向きから等方的に仮想粒子が到着し、衝突散乱後、すべての向きに等方的に仮想粒子が出発していると考えられる。もし、流体がある向き（順向き）にマクロな速度をもっていれば、順向き速度成分が正の速度をもった仮想粒子の方が、逆向き速度成分が正の速度をもった仮想粒子よりも、頻繁に到着し、かつ、頻繁に出発しているであろう。
次に、個々の格子点において、一直線上で相対する前後２つの向きへの仮想粒子の移動を考える。ある瞬間に順向きの到着粒子と出発粒子がともに存在し、逆向きの到着粒子と出発粒子がともに存在しない場合、局所的な流れの主流が瞬間的に順向きであると考える。このとき、格子点から逆向きに出発しようとした仮想粒子（緩和粒子）の一部について、主流に巻き込まれて速度を反転させ順向きの格子点に向かって並進していくように操作する。これを、ここでは、“局所主流巻込操作（仮称）”と呼ぶことにする。この場合のシミュレーション結果を図図７に示す。

格子ガス法流体解析モデルとニューラルネットワークの融合

巻 込 操 作 適 用 率 0%

― 47 ―

巻 込 操 作 適 用 率 3%

(LGA06234-50)

巻 込 操 作 適 用 13%

巻 込 操 作 適 用 率 27 %(LGA06225-50)(LGA06228 50)

図７．“局所主流巻込操作(仮称)”の適用確率を変化させた場合の円柱後流の変化（挙動が定常化した時点におけるある瞬間(12800 時刻ステップ目)の疎視化された運動量）“局所主流巻込操作(仮称)”の適用確率を、0%､3%、13％、27%と変化させた場合、円柱後流の過渡変化挙動がかなり変化する。ここでは、挙動が定常化した時点におけるある瞬間（12800 時
刻ステップ目）の疎視化された運動量を示す。同一格子点幅で同じ時刻ステップ経過後の状態であるが、シミュレーションを行った適用確率の範囲では、適用確率を上げるにしたがって生じる
円柱後流の変化は、レイノルズ数が増加した場合に現れる変化に近く、定性的には、流体粘性が低下していることがわかる。
（７）流体粘性を変える“近傍粒子引合操作(仮称)”の適用効果を観る。
次に、“個別粒子への達人操作”として、「近傍粒子引合操作（仮称）」なるものを考え、確率的に適用してみた場合の結果を示す。
この場合も、まず、個々の格子点において、一直線上で相対する前後２つの向きへの仮想粒子の移動を考える。ある瞬間に順向きの出発粒子と逆向きの到着粒子がともに存在し、逆向きの出発粒子と順向きの到着粒子がともに存在しない場合、局所的には、その瞬間、順向き前方の方が後方よりも多くの仮想粒子が存在していたと考える。このとき、格子点から逆向きに出発しようとした仮想粒子（緩和粒子）の一部について、前方に多く存在する粒子からの粒子間引力が勝った結果、その速度を反転させ順向きの格子点に向かって並進していくように操作する。これを、ここでは、“近傍粒子引合操作（仮称）”と呼ぶことにする。この場合のシミュレーション結果を図８に示す。
“近傍粒子引合操作(仮称)”の適用確率を、0.1%､0.5%、2%、7％と変化させた場合、死水域の面積が変化する。また、このとき、円柱に衝突してバウンズバックする仮想粒子の X 成分運動量の変化から、円柱が流体から受けている抗力を求めた。この相対値は、それぞれ、5709, 4131, 3207,
2511 であり、死水域の減少とともに抗力が大きく減少した。円柱の配置体系、入口流入速度、格子点間隔等はすべて同じであるにも関わらず、抗力が大きく変化したことから、流体粘性が大きく変化したと考えられる。

― 48 ―

SENAC Vol. 54, No. （1 2021. 1）

引 合 操 作 適 用 率 0. 1%

引 合 操 作 適 用 率 0. 5%(LGA060615 50)(LGA060613-50)

引合操作適用率 2%

引合操作適用率 7%(LGA060610 50)(LGA06065 50)

図８．“近傍粒子引合操作(仮称)”の適用確率を変化させた場合の円柱後流の変化（挙動が定常化した時点におけるある瞬間(12800 時刻ステップ目)の疎視化された運動量）

４．まとめ（１）今回試みた格子流体の粘性制御法に関する数値シミュレーションの結果について
今回試みた格子流体の粘性制御法は、『時刻ｔにおける衝突散乱の結果各方向に出発しようとする“緩和粒子”に対して、当該格子点に時刻ｔに到着した同方向の仮想粒子と時刻ｔ－１に出発した同方向の仮想粒子の有無に関する情報だけを観察し、ある条件が整った際にある確率で速度を反転させて出発させてやる』というものである。この方法の適用確率をいろいろ変えて円柱後流の挙動の違いを観察した結果、円柱の配置体系、入口流入速度、格子点間隔等はすべて同じであるにもかかわらず、円柱後流のうずの発生挙動や円柱に働く抗力が大きく変化することを確認した。現時点では定性的な評価にとどまるが、今回試みた方法は、流体粘性をかなり大きく変化させる効果がある。従って、格子ガス法において、高レイノルズ数領域の流体挙動を比較的少な
い格子点数でも短時間に計算できる手法を開発する重要なヒントになるものと考えられる。
なお、流体粘性の制御については、各格子点における各方向の速度をもつ到着粒子と出発粒子の頻度をモニタリングすることでその時点の流体粘性を把握し、そのとき保持したい流体粘性の値に近づける操作をリアルタイムフィードバックで行うことを想定している。
（２）格子ガス法流体解析モデルとニューラルネットワークの融合について
本稿では、格子ガス法流体解析におけるすべての計算が、整数型多層パーセプトロンの組み合
わせ計算として実行できる可能性を示した。
“整数型”にこだわった理由は、実数計算に比べて計算速度の向上やメモリ節約の点で有利なこともあるが、格子ガス法がもつ最も重要なメリットである「誤差が発生しない安定した時間発展計算の実現」を大切にしたかったからである。これは、将来、乱流の解明に貢献できるかもしれないという期待感がある。
また、筆者らが行ってきたこれまでの研究では、ニューラルネットの学習原理を活用した「不完全な実世界情報でも利用可能なリアルタームデータ同化」という第一に掲げた目標を達成できていない。今後、この方向のシミュレーション事例を積み上げ、“格子ガス法流体解析モデル”と“ニューラルネットワーク”の融合が実用的な流体解析に貢献できることを示せたら幸いである。

格子ガス法流体解析モデルとニューラルネットワークの融合― 49 ―

謝辞
本稿で述べた研究課題では、“格子ガス法流体解析モデル”と“ニューラルネットワーク”の融合を目指して、いろいろなアイデアを多数試すことが非常に重要であった。この時期に、個々の計算モデルを試すたびに計算コードのチューニングに多くの時間をさくことは非効率的で悩ましい。この点、東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュータ SX-ACE
は、特別なチューニングをしないでも十分な計算速度を確保することができ、大変助かった。今後、さらに進化したベクトル機である SX-Aurora TSUBASA の利用が楽しみである。
また、利用にあたって同センター関係各位のご親切なご指導とご協力をいただき、心から感謝する次第である。今後とも同センターの有意義な活動を継続的に発展させられることを期待する。
参考文献[1] Chistopher M. Teixeira, “Continuum Limit of Lattice Gas Fluid Dynamics”, MIT, 1993
[2] 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
[3] 松岡, 菊池,“多速さ格子ガス法実用化展開への手がかり”, SENAC Vol.49 No.4, pp.1-15,
2016[4] 松岡,“ビット演算による CFD と等価な高精度流体解析手法”, RIST News No.64, pp.17-28,
2018[5] 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
[6] Michel Hénon, “Viscosity of a Lattice Gas”, Complex Systems, 1(1987), pp.763-789,
1987[7] 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
[8] Daniel H. Rothman, “Negative-Viscosity Lattice Gases”, Journal of Statistical Physics,
Vol.56, Nos. 3/4, 1989
[9] 松岡,菊池,“コンパクトな計算機によるリアルタイム流体解析の実現に向けて”, SENAC
Vol.51 No.2, pp.1-10, 2018[10] 松岡,菊池,“仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性”, SENAC
Vol.52 No.2, pp.18-27, 2019[11] 松岡,菊池,“リカレントニューラルネットワークによる実世界流れ場解析用時間発展計算モデルの探求”, SENAC Vol.53 No.1, pp.25-33, 2020