— 1 —SENAC Vol. 49, No. 4（2016. 10）［共同研究成果］

多速さ格子ガス法実用化展開への手がかり―“テシャラの方法”から“マクスウェルの達人操作”へ―

松岡 浩
東北大学電気通信研究所(客員)
一般財団法人高度情報科学技術研究機構
菊池 範子
株式会社サイエンス・サービス

1.

はじめに

もう 20 年以上昔(1993 年)になるが、テシャラ(Chistopher M. Teixeira)は、彼の MIT における博士論文“Continuum Limit of Lattice Gas Fluid Dynamics[1]”において、CFD より格段にすぐれた多速さ格子ガス法を提案した。同じ格子点規模の流体シミュレーションを行うのに、CFDに比べて、1/10 のメモリサイズで 100 倍近い並列高速計算ができるという。筆者らがこれを知ったのは、10 年前のことであるが、以後、東北大学サイバーサイエンスセンターの SX シリーズを中心にしたベクトル型スーパーコンピュータを利用させて頂き、その実用化研究に取り組んできた。本稿では、テシャラの方法の概念と、筆者らが考える今後の方向性について紹介する。従って、学術論文の概要紹介とは異なり、一種の読み物として掲載させていただくものである。なお、
本稿の本文は松岡が執筆し、掲載したシミュレーション事例に関する計算及び可視化はすべて菊池が分担した。

２.

本稿の概要

格子ガス法は、流体が存在する空間中に格子を張り、多数の仮想粒子が、その格子点上のみで他の仮想粒子と衝突して進行の向きを変えながら格子点間を移動していく様子を平均化(疎視化)
して流体の挙動を模擬する方法である。このとき、仮想粒子がもつ質量・運動量・エネルギーが衝突の前後で保存されるような粒子衝突を想定する限り、その挙動は、自然界におけるある条件
下の流体挙動とかなり似かよったものになる。しかしながら、多少の違いが生じる。この主な原因は、仮想粒子がもちうる速度が、流体分子がもちうる速度とは異なり、制約があるからである。
すなわち、仮想粒子は、格子点が存在する位置の間しか移動できないので、速度の大きさも向きも離散的な値をとり、自然界の流体分子のように連続的な値をとることができない。この制約によって、格子ガス法が導くマクロな挙動は、連続流体を仮定している CFD が導く流体挙動とは多少異なったものになる。
“テシャラの方法”は、多速さ格子ガス法において、異なる速さをもつ仮
想粒子どうしの衝突頻度をうまく調節することによって、この“多少の差異”を解消する方法である。これにより、多速さ格子ガス法は、CFD を代替できる精度をもったシミュレーション手法となりうる。本稿では、多少の差異を解消するテシャラの方法の概念を身近な大きさのシミュレーション空間に対して考察してみる。このとき付記したシミュレーション事例は、概算によって数値の大きさを感覚的につかむことが目的である。
テシャラの方法は、個別の粒子を識別して、それに応じた操作を行い、自然界では生じないマ
クロな挙動を自己組織化により生成する。これは、負のエントロピーを生み出す“マクスウェルの悪魔”の仕事に似ている。ここでは、上記の“多少の差異”をピッタリとキャンセルする“達人技”のイメージを強調して、このような個別粒子に対する意図的な操作を“マクスウェルの達人操作”と名付ける。今後、この種のシミュレーション手法の実用化展開に当たって、自由な発想による“マクスウェルの達人操作”の導入がひとつの重要な手がかりになるであろう。

— 2 —

SENAC Vol. 49, No. 4（2016. 10）３.

格子ガス法による仮想世界と現実空間の対応

3.1

現実空間における４次元直交格子の配置

格子ガス法では、格子が張られた仮想世界の中を、仮想粒子が動きまわる。この仮想世界と現実世界の対応を実感できるように、あえて、具体的な数値と単位を使う。
まず、身近な大きさをもつ空間として、人ひとりがゆったり入れるぐらいの空間：“一辺(L)が2.56m の立方体（L×L×L）”（L＝2.56m）を考える。次に、この空間を“格子間隔(⊿ℓ)が 1cm の直交格子”（⊿ℓ＝1cm）で埋め尽くす。直交格子の x,y,z 軸は、立方体空間の各辺に平行になるよ
うに配置する。従って、３次元格子点は、x,y,z 方向にそれぞれ約 256 個並ぶので、３次元格子点の総数は全部で 256×256×256≒1678 万個となる。我々がシミュレーションで知りたい流体挙動の解像度(⊿L)は、もちろん目的によるが、ここでは、“手の指の長さ程度の解像度 8cm”（⊿L＝8cm）としておく。ただし、ここで注意が必要である。テシャラは、３次元シミュレーション用の直交格子として、対称性に優れた“４次元面心超立方体(FCHC: Face-Centered Hyper Cubic)格子”を用いている。従って、上記のすべての３次元格子点には、第４番目の自由度をもった座標を割りあてる必要がある。イメージ的には、各３次元格子点には、その格子点を中心にして半径 r が異なる多層の球殻が配置されており、仮想粒子は、半径 r が異なる球殻に別々に存在しうるものとする。こうして、４次元格子点の位置座標は（x,y,z,r）で与えられることになる。ここで、４次元目の座標位置 r も長さの単位で測ることができ、４次元格子点は、x,y,z 座標と同様に 1cm間隔で配置される。座標軸 r のとりうる範囲は、現実のシミュレーション対象である３次元物体の制約を受けないので任意性があるが、ここでは、r＝0,1,2,3cm の４点のみから成る“循環座標”としておく。すなわち、仮想粒子が、r＝3cm の位置にいて、さらに r のプラス向きに移動したい
ときは、r 座標は r＝0cm につながる。逆に、r＝0cm の位置にいて、さらに r のマイナス向きに移動したいときは、r 座標は r＝3cm につながると考える。

3.2

仮想粒子をめぐる物理量の設定

仮想粒子の速さについては、「静止粒子」、「遅い粒子」、「速い粒子」の３種類を仮定する。また、
��“仮想粒子の質量(m)はどれも同じで 1μg”(m ＝ 10 kg)とする。格子ガス法では、多数の仮想粒子が、同じ時間間隔(⊿τ秒)で時刻ステップがひとつ刻まれるたびに、格子点から格子点へと飛び移る。ここでは、この移動過程に要する時間⊿τを千分の 1 秒(⊿τ＝0.001sec)に設定する。
仮想粒子がもちうる速度については、注意が必要である。４次元 FCHC 格子を採用しているので、各格子点は４次元空間中に存在するため、格子点の間を飛び回る仮想粒子は、速度成分を４つも��� ＝(���� ,����� ,����� ,����� ) m/sec の形になる。ここで、
つことになる。すなわち、“４次元速度”�����
下付き添え字の「j」は、仮想粒子がもつ速度の大きさ(速さ)の違い、すなわち、エネルギーの違いを表す。j＝0 なら｢静止粒子｣、j=1 なら｢遅い粒子｣、j=2 なら｢速い粒子｣とする。また、「i」
は速度の向きの違いを表す。ただし、「静止粒子」の「i」は、単に複数の静止粒子の区別を表す。
仮想粒子の速さは、「格子点間距離 ℓ[1cm]を時刻ステップの刻み間隔⊿τ[0.001sec]で移動する速さｃ」のオーダーである。ここでは、この c＝ℓ/⊿τ＝1cm/0.001sec＝10 m/sec を基準にし��� を、次のように具体的に与えてしまう。エネルギーも求まる。
て、仮想粒子��がもちうる速度 �����

・「静止粒子（j=0,i=1～8）」がもちうる速度： 粒子がもつエネルギー(運動エネルギー)は、もちろんε� ＝0[J]。

������
��� ＝（0,0,0,0）m/sec、静止粒子は複数想定してもよい。通常、１個または６個の静止粒子を仮定することが多いが、ここでは、８個としておく。

・「遅い粒子（j=1,i=1～24）」がもちうる速度：

�
�
��
�
��
������
������
エネルギーは、ε� ＝ m(��� ･��� )＝mc ＝10 ･10 ＝10 [J]。
�— 3 —

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —������
��� ＝（±c,±c,0,0）m/sec，（±c,0,±c,0）m/sec，（±c,0,0,±c）m/sec，（0,±c,±c,0）m/sec，（0,±c,0,±c）m/sec，（0,0,±c,±c）m/sec

・「速い粒子（j=2,i=1～24）」がもちうる速度：の 24 種類。
����������
������
エネルギーは、ε� ＝ �m(��� ･��� ) ＝ 2mc ＝2×10 [J]。

������
�
（0,±2c,0,0）m/sec，（0,0,±2c,0）m/sec，（0,0,0,±2c）m/sec，
�� ＝（±2c,0,0,0）m/sec，（±c,±c,±c,±c）m/secの 24 種類。

仮想粒子が、４次元格子点上のどこかに存在すれば、上記の速度をもって移動する限り、仮想粒子は設定された４次元 FCHC 格子から逸脱することはない。ただし、体系の周辺境界や固体境界等では、特別な移動規則を与える必要がある。また、仮想粒子は、格子点上でのみ、他の仮想粒子と衝突して移動の向きや速さを変えることができる。このとき、衝突の前後で、衝突に関与する仮想粒子合計の粒子数(質量)、運動量及びエネルギー(運動エネルギー)が保存される衝突規則を適用する。衝突には、｢遅い粒子｣どうし、｢速い粒子｣どうし、｢遅い粒子｣と｢速い粒子｣の衝突などのように、“仮想粒子数のエネルギー分布を変えない衝突”と、｢静止粒子｣と｢速い粒子｣が２個の｢遅い粒子｣と入れ替わる“仮想粒子数のエネルギー分布を変える衝突”の２種類がある。
最後に、重要な仮定として、「各４次元格子点においては、静止粒子を除き、複数の仮想粒子が同じ４次元速度をもつことはできない」という制約を課す。

3.3

マクロな挙動を観察する時空間スケール

仮想粒子の動きを平均化(以下、“疎視化”という。)してマクロな挙動を導出したい。まず、疎視化する範囲を決める。3.1 に前述したとおり、現実の３次元空間の中で求める解像度は、⊿L＝8cm とした。従って、間隔⊿ℓ＝1cm で４次元空間中に配置された格子点については、(x 方向 8cm)
×(y 方向 8cm)×(z 方向 8cm)×(r 方向 4cm)の４次元空間中に存在する 8×8×8×4＝2048 個の４次元格子点ごとに、そこに存在する仮想粒子の動きをまとめればよい。最後が「×4」になっているのは、第 4 番目の次元である r 方向には、４次元格子点が r=0,1,2,3 の位置に４つ存在するだ
けで、しかも、これらはすべて同じ３次元位置(x,y,z)に属しているからである。
また、時間についても、⊿τ＝0.001 秒という短い時間では、各仮想粒子は、⊿ℓ＝1cm のオー
ダーの距離しか移動できず、⊿L＝8cm の解像度が求められるマクロなオーダーでは大きな変化は生じない。特に、L＝2.56m の立方体内部における大きなスケールの挙動を追跡する場合は、⊿τ＝0.001 秒よりもかなり長い時間間隔⊿T 秒でスナップショットをとっていけばよいであろう。このスナップショット間隔⊿T の選択には、かなりの任意性があるが、ここでは、約 0.5 秒（⊿T＝0.512sec）とした。これは、仮想粒子の移動を 512 時刻ステップだけ刻める時間であり、仮想粒子が速さ c＝10m/sec で、立方体空間の一辺方向に他の仮想粒子と全く衝突せずに直進できたとしたら、二辺分の長さだけ進行できる時間である。例えば、初期状態から、⊿T＝0.512 秒ごとに 13
回のスナップショットを観察すれば、一辺が 2.56m の立方体空間における、約 6 秒間(6.144sec)のマクロな挙動を観察できることになる。

3.4

サブマクロな世界の記述

現実世界の流体を構成する流体粒子は、実際は、分子レベルの非常にミクロな世界を運動している。他方、仮想粒子が動き回る直交格子は⊿ℓ＝1cm の世界であり、現実の流体粒子は、そのよ
うな格子に動きを制約されることはない。この意味で、仮想粒子が動き回る世界は、現実にはありえない“仮想世界”である。しかし、この仮想世界における仮想粒子の動きを疎視化すれば、解像度⊿L＝8cm 程度のレベルで、現実のマクロな世界における流体挙動が得られることを期待しているわけである。ここで、仮想粒子が動き回る仮想世界は、スケール的に見ると、マクロな世界のすぐ下の“サブマクロな世界”であるということができる。

— 4 —

SENAC Vol. 49, No. 4（2016. 10）

サブマクロな世界の状態を完全に記述するには、各４次元格子点に、仮想粒子��が実際に存在するか否か？を示してやればよい。前述したように、各４次元格子点(x,y,x,r)には、①複数個の「静止粒子」（ここでは、これを仮に 8 個とした。）、②24 個のいろいろな向きの速度を持ちうる「遅い粒子」、③24 個のいろいろな向きの速度を持ちうる「速い粒子」が存在しうる。従って、ひとつの４次元格子点には、最大 8+24+24＝56 個の仮想粒子が存在しうる。ここでは、仮想粒子��が存在するか否か？を「1」または「0」の２値数(��� )で表現する。ただし、�� は、最大 56 個存在しうる仮想粒子の識別番号である。
��� (x, y, z, t): ４次元 FCHC 格子上の位置(x, y, z, t)にある４次元格子点において、時刻ステップ
ｔに仮想粒子��が存在する場合「1」、存在しない場合「0」の値をとる２値数。これにより、
格子ガス法において仮想粒子が動き回る４次元 FCHC 格子という“サブマクロな世界”の状態を完全に記述することができる。

3.5

マクロな状態量の導出

サブマクロな世界を記述する状態量 ��� から、マクロな世界を記述する状態量を導出する。はじめに、“質量密度ρ”と“運動量密度 P”を次式で ��� から定義する。
[kg/�� ]・“質量密度”： ρ(x, y, z, t)≡m･f ∑������ m･��� (x, y, z, t)
ここで、m：仮想粒子１個の質量[kg/個]

1μg＝10�� kg

f：４次元単位体積(１�� )あたりに存在する４次元格子点の数[個/�� ]。

４次元格子点は、４次元空間中に⊿ℓ＝1cm の間隔で並んでいるので、４次元の単位体積中には、f＝(1/⊿ℓ)･(1/⊿ℓ)･(1/⊿ℓ)･(1/⊿ℓ)＝10� 個/�������・“運動量密度”： ����������������������
P(x, y, z, t)≡m･f ∑���������� ･��� (x, y, z, t)}
�����
ここで、�
�� ：

[kg/㎥ sec]

仮想粒子�� がもつ４次元速度ベクトル[m/sec]

成分で書けば、 �� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)}[kg/㎥ sec]�� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)}[kg/㎥ sec]�� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)} [kg/㎥ sec]�� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)}[kg/㎥ sec]�����
ここで、���� , ���� , ���� , ���� は、各々��� のｘ成分,ｙ成分,ｚ成分,ｒ成分［m/sec］である。

次に、“マクロな速度”は、“質量密度”と“運動量密度”から次式で定義しておく。
����������������������
y, z, t)／ρ(x, y, z, t)
・“マクロな速度”： ����������������������
U(x, y, z, t)≡P(x,

��
�����
＝〔m･f ∑���������� ･��� (x, y, z, t)}〕／{m･f ∑���� ��� (x, y, z, t)}�������＝ 〔∑��
�������� ･��� (x, y, z, t)}〕／{∑���� ��� (x, y, z, t)}[m/sec]・・・mや f に依存しない。

��
成分で書けば、�� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]���� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]���� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]���� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —

3.6

— 5 —

マクロな速度の大きさの事例評価“仮想格子流体”が一様に流れている状況を考察する。このため、ｘ＝0 の面に存在する全ての４次元格子点から＋ｘの向きに、４次元速度が(ｃ,0,0,-ｃ)（ｃ,0,0,ｃ）（2ｃ,0,0,0）（ただし、ｃ＝10m/sec）となる３種類の仮想粒子を適当な頻度(確率 1/4)で送りこんでみよう。このときの初期条件としては、立方体内部のすべての４次元格子点において、1/4 の存在確率でランダ
ムに仮想粒子が存在しているものとする。このとき、質量密度は、次式の値となり、数値的には、常温常圧の空気の密度にほぼ等しい。もちろん、ｍの値の設定の仕方で、どんな物質のどんな状態の質量にも、値を対応させることができる。
��
kg×10� 個/�� ×(56×1/4)＝1.4[kg/�� ]ρ(x, y, z, t) ＝ m･f ∑������ �� (x, y, z, t)＝10

また、ｙ方向とｚ方向には周期的境界条件を課し、ｘ方向出口では、出口直前に位置する格子
点状態 ��� を出口直後に位置する格子点状態 ��� に上書きする方法で境界条件を与えた。

時刻ステップ＝0 の初期状態では、いろいろな向きの速度をもつ仮想粒子が、等方的に同じ確率(0.25)で各格子点に存在する。従って、最初は、マクロな速度はゼロである。しばらくすると、＋ｘ向きの運動量を得て＋ｘ向きの“流れ”が発生する。⊿τ＝0.001 秒の時刻刻みを約 1000 回(1024 回)繰り返した１秒後の時点では、多少のゆらぎはあるが、マクロな速度としては、全体が一様に＋ｘ向きに流れている。シミュレーション結果は省略するが、このときの仮想粒子の存否
データ：��� (x,y,z,r)からｘ向きのマクロな速度を算出してみると、約 0.5m/sec（数値データの意味は省略するが、110×10／2300≒0.48m/sec）であることがわかる。
一般に、４次元格子点において、仮想粒子がもつ速度の向きをそろえてある向きのマクロな速度をできるだけ大きくしようとしても、その最適な向きに存在できる仮想粒子は高々１個である。
また、最適な向きをあきらめて似たような向きの速度をもつ仮想粒子を増やしていくと、最適な向きとはますます異なった周辺向きの速度をもった仮想粒子の存在場所しか空いていなくなる。
従って、個々の仮想粒子の速度をある特定の向きに束ねるには限界がある。この結果、マクロな速度の大きさＵm/sec は、ｃ＝10m/sec よりもかなり小さい値になる。すなわち、Ｕ<<ｃという関係が多くの場合に成立している。

４.波型せん断流の減衰解析
4.1

はじめに

第３章で、格子ガス法による仮想世界を現実世界に対応させて、４次元速度が(ｃ,0,0,-ｃ)（ｃ,0,0,ｃ）（2ｃ,0,0,0）（ただし、ｃ＝10m/sec）となる３種類の仮想粒子を確率 1/4 の頻度で＋ｘ向きに送りこむと、１秒後には、一辺が 2.56m の立方体中における静止した“仮想流体”が、
約 0.5m/sec のマクロな速度を得ることを予想した。格子ガス法によるサブマクロな世界では、質量、運動量、エネルギーが保存されている。従って、それを疎視化して得られた“仮想流体”も、これらが保存されることで、実際の流体に似た挙動を示すであろう。しかし、このときの質量密度が、たまたま常温常圧の空気の密度に等しいからと言って、この“仮想流体”が常にマクロな世界における実際の空気の流れを模擬できていると即決はできない。ここでは、もう少し複雑な状況、すなわち、粘性が関与する状況を考察してみる。

4.2

波型せん断流の事例計算

ここでも、まずは、疎視化したときに「＋ｘ向きにマクロな速度 �� (x, y, z, t)＝約 0.5m/sec で一様に流れているマクロな状態」を再現できる“仮想格子流体”を考える。次に、ある時点で、瞬間的に適当な力を加えて、マクロな速度のｙ方向成分 �� (x, y, z, t) を波型にしてみる。これによってｙ方向に粘性によるせん断力が働くことになる。この瞬間の３次元画像を図１に示す。こ— 6 —

SENAC Vol. 49, No. 4（2016. 10）の状態から時刻ステップ⊿τで時刻刻みを進めていくと、マクロな速度のｙ成分 ܷ௬ の速度勾配が減少しながら、マクロな速度のｘ成分 ܷ௫ に乗って、＋ｘ向きに移動する“波型せん断流の減衰過程”がスナップショットで得られる。この減衰から、粘性の定量値等を求めることができる。

図１

瞬間的に生じさせた波型せん断流（ParaView による可視化）

シミュレーションによって得られた“波型せん断流の減衰過程”を図２に示す。図２には、比較のために、格子点間隔が⊿ℓ＝1cm の場合以外に、⊿ℓ＝2cm と粗くした場合、及び、⊿ℓ＝0.5cm、0.25cm と細かくした場合の結果を一緒に示してある。また、図２において、各ケース名の下に付記したνの数値は、動粘性係数の大ざっぱな概算値である。この算出方法は、4.5 に後述する。
一般に、格子点間隔を短くしていけば、仮想粒子の平均自由行程が短くなり、粘性は小さくなる。実際、図２から、格子点間隔が短い「ケース４」に近づくほど、波型せん断流の減衰の仕方
が緩慢になっており、粘性が小さくなっていることを見てとれる。
ケース１：
格子点間隔⊿ℓ＝2cm の場合
ν≒0.04 [㎡/秒]

格子点数：210 万
計算時間：
SX-ACE 32 ノード
約 14 秒
ケース２：
格子点間隔⊿ℓ＝1cm の場合

200

T=0

100
0

T=1
T=2
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31

-100

200

T=0
T=1

100
0
-100
-200

T=4
T=5

-200

ν≒0.0264 [㎡/秒]

格子点数：1678 万
計算時間：
SX-ACE 32 ノード
約 130 秒

T=3

T=2
1 3 5 7 9 1113151719212325272931

T=3
T=4
T=5

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —

ケース３：

200

格子点間隔⊿ℓ＝0.5cm の場合

T=0
T=1

100

ν≒0.0216 [㎡/秒]0

格子点数：1 億 3422 万

計算時間：
SX-ACE 32 ノード
約 27 分 30 秒

-100

T=2
1 3 5 7 9 1113151719212325272931

200

格子点間隔⊿ℓ＝0.25cm の場合
ν≒0.0056 [㎡/秒]
格子点数：10 億 7374 万

計算時間：
SX-ACE 32 ノード
約 8 時間 43 分
図２

-100

T=4

T=0
T=1

100
0

T=3
T=5

-200

ケース４：

— 7 —

T=2
1 3 5 7 9 1113151719212325272931

-200

T=3
T=4
T=5

波型せん断流の減衰過程（マクロな速度の y 方向成分 �� の変化）

（横軸は x 方向の位置[m]、縦軸はマクロな速度のｙ方向成分[m/sec]で、値は相対値である。）（T=0：静止状態、T=1：一様＋X 向き流れ、T=2：�� に対して波型速度分布を瞬間的に付与、T=2 以降：減衰。）

4.3“仮想格子流体”が導くマクロな挙動が満足する方程式

サブマクロな世界における４次元 FCHC 格子上の仮想粒子の衝突・移動という“仮想格子流体”のモデルは、それを疎視化してマクロなレベルで眺めると、現実の流体挙動に似た状態を再現できる。しかし、実際には、“多少の差異”が生じる。すなわち、実際の流体現象を正確に表現できると言われているナビエ・ストークス方程式の解との間には、ある差異が生じる。
テシャラの論文にも示されているが、ここで想定している“仮想格子流体”を疎視化した場合に、マクロなレベルで得られる方程式を導出してみると、次式の形になる。これは、ナビエ・ス
トークスとは少し異なる。
��＋�Ｕ･∇Ｕ ≒ νΔＵ��
これを、4.2 の事例計算で仮定した初期条件：
�� (�� �� �� �) = �� 、�� (�� �� �� �) = �� ��� ��、�� (�� �� �� �) = 0の下で解くと、次の過渡変化の理論解が得られる。
これが、今回想定した“仮想格子流体”の動きを疎視化した後に得られるマクロな挙動である。
�� (�� �� �� �) = ����� (�� �� �� �) = �� ���� � ･���(�� � ��� ��)

4.4

�� (�� �� �� �) = 0“仮想格子流体”に関するパラメータ�の評価“仮想格子流体”が満足する方程式に含まれるパラメータ�の値については、シミュレーショ
ン結果として得られた“波型せん断流の減衰過程”における位相の解析から、定量的な評価を行うことができる。

— 8 —

SENAC Vol. 49, No. 4（2016. 10）

上述の �� (�� �� �� �) の理論解によれば、スナップショット取得時刻Ｔにおける振幅ピーク位置
は、������ � ��� �Ｔ�＝1 を満たすｘで与えられる。実際、�� � ��� �� ＝��� から一番最初のピーク
位置を知ることができ、�������� ＝�� ��＋���� であることがわかる。同様にして、スナップショッ
ト取得時刻Ｔ＋⊿Ｔにおける最初のピーク位置は、�������(��⊿�) ＝�� �(� � ⊿�)＋���� である。そのピーク位置の差は、�������(��⊿�) － �������� ＝�� �⊿� となり、この値は k やＴには依存しない。ここで、�� はシミュレーション結果から測定(�� ＝0.48m/sec)でき、⊿Ｔは既知(⊿Ｔ＝0.512 秒)であるから、ピーク位置の差の測定値から � の値を求めることができる。
シミュレーション結果より �������(��⊿�) － �������� ≒0.16ｍであるから、次の値を得る。

�＝�������(��⊿�) － ��������
�� ⊿�

＝

�����

(���������)�(��������)≒ 0.65

さて、前述した“仮想格子流体”が満足する方程式：
��＋�Ｕ･∇Ｕ = νΔＵ��は、�=1 の場合、“一定密度・一定温度の場合の非圧縮性流体に関するナビエ・ストークス方程式”に等しくなる。すなわち、今回計算した“仮想格子流体”の挙動は、�≠1 であるので、実在の流体の挙動とは � の差に応じた差異がある。後述するように、テシャラは、静止粒子と速い粒子が衝突して２個の遅い粒子になる衝突の頻度を制御することによって、この差異を解消する方法を考案した。

4.5“仮想格子流体”に関するパラメータ�の評価“仮想格子流体”が満足する方程式に含まれるパラメータλについても、シミュレーション結果として得られた“波型せん断流の減衰過程”における振幅の解析から、定量的な評価を行うこ�とができる。前述の �� (�� �� �� �) の理論解によれば、波型せん断流の振幅ピーク値の減衰は、� ��� �に従う。ここで、νは“仮想格子流体の動粘性係数[㎡/秒]”、k は波数[1/m]、t は経過時間[秒]である。
“波数”とは、波型せん断流の初期条件である正弦波の単位長さあたりの位相進行量である。今回の波型せん断流の初期条件では、立方体を構成するＸ方向の一辺 L=2.56m の長さの中に２つの波を配置したので、k＝4π/L＝π/0.64[1/m]であることがわかる。
�
さて、あるスナップショット取得時刻Ｔにおける振幅ピーク値の理論解は �� � ��� � で、その次�のスナップショット取得時刻Ｔ＋⊿Ｔにおける振幅ピーク値は �� � ��� (��⊿�) である。従って、そ�の減衰比は、�� やＴに依存せず、� ��� ⊿� で与えられることがわかる。ここで、k と⊿Ｔは既知であるから、シミュレーション結果から測定できる振幅ピーク値の減衰比から仮想格子流体の動粘性係数νの値を求めることができる。
いま、格子点間隔⊿ℓ＝1cm の場合を例にとる。波型せん断流の減衰過程における最初のスナッ
プショット取得時刻Ｔでの最初の振幅の相対ピーク値は約 150 で、Ｔ＋⊿Ｔでの相対ピーク値は
約 100 になっている。従って、波型せん断流の振幅ピーク値の減衰比は、100/150≒0.67 である。
�
これより、0.67≒� ��� ⊿� からνの値を概算することができ、ν≒0.0264[㎡/秒]という値を得る。これは、空気の動粘性係数のレベルよりもかなり大きい。従って、テシャラの方法によって、
�=1 の補正がなされ、これが何かの条件下における実在流体の挙動を模擬できていたとしても、粘性値を見ると、それは空気ではなく、うすい“はちみつ”ぐらいに対応するものだとわかる。

4.6“仮想格子流体”を考えるアプローチからの視点

さて、岩波講座「現代工学の基礎」の“計算熱流体力学”に“物理現象のメゾスコピックモデル”に関して以下の記述[2]がある。

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —— 9 —「流体のように、多数の分子とその相互作用でできている系は、分子間の相互作用と相互作用の間の相関がすぐに薄れてしまうという特性(分子カオス)をもつ。このため全体としてのマクロな挙動は、相互作用のミクロな詳細には依存しなくなる。相互作用の詳細は、マクロな係数には影響するが、方程式の代数構造には現れなくなる。このことを逆に考えれば、多数の分子相互作用の集積として現れる物理現象を、連続体近似でもなければ分子レベルでもない異なるレベルでモデル化することが可能であることがわかる。マクロなレベルに疎視化したときに同じマクロな挙動を示すように仮想的なレベルで相互作用をモデル化できれば、それをマクロな現象のモデルとして使えることになる。
」（cf.計算熱流体力学、笠木伸英、松本洋一郎、大橋弘忠、岩波講座、現
代工学の基礎、p108）
上記の記述の下線部を前述の“波型せん断流の解析”の場合にあてはめると、以下のことに対応していることがわかる。すなわち、“仮想格子流体”を疎視化して得られるマクロな挙動も、“実在流体”も、その支配方程式の代数構造は、
��＋�Ｕ･∇Ｕ = νΔＵ��という共通の形をしている。違いは、マクロな係数 � のみに現れたわけである。

５.マクロな係数の差異を解消するテシャラの方法
5.1

はじめに

マクロな係数の差異を解消する“テシャラの方法”の概念を述べる。本来は、サブマクロな世界における仮想粒子の運動から、マクロな挙動が導出される全過程を詳細に追跡する必要がある。
しかし、その導出過程はかなり長く複雑で、この部分を厳密に記載した参考論文([1]、[3]、[4]等)も当然存在するので、テシャラの方法が導入されるところまでだけの“大筋”を示す。具体的には、「局所平衡状態における仮想粒子分布を導出する式(後述の⑧式)」までを記述し、以後の過程は省略した。なお、“大筋”も、筆者なりの独自解釈で簡素化してあり、文献どおりではない。

5.2

サブマクロな世界における仮想粒子分布の記述

時刻 t に４次元空間内のある格子点位置(x,y,z,r)に到着した仮想粒子の“到着粒子分布”を����������������������������
�(���� �� �� �) 、衝突後の“出発粒子分布”を ������������������������������(�� �� �� �� �) と置く。これらの状態ベクトルの ji 成分

をそれぞれ ��� (�� �� �� �� �)、�� �� (�� �� �� �� �)と書き、仮想粒子 ji が存在するときは「1」、存在しないときは「0」と定義する。ここで、仮想粒子 ji とは、エネルギーが j で速度の向きが i の仮想粒子を意味する。ただし、静止粒子の場合(j=0)は速度の向きの区別はないので、i は、単に粒子の区別を表す。今、ji は 56 通りを考えるので、状態ベクトルは 56 ビットの“ビット列”で表現でき����������� (�＝1る。一般に、任意の 56 ビットのビット列について、全ての可能なビットパターンを、�(�)��～ 2 のいずれか)と書き、その ji 成分を ��� (�)と表現すれば、次式が成り立つ。
��� (�� �� �� �� �) ＝����〔��� (�) × ( �������������������������
�(�� �� �� �� �) �� �����������
�(�) )〕

���

・・・①式

����������������������������
������������������������������ �� �� �) �� ������������
�(�) ) は、２つのベクトル �(���� �� �� �) と ������������
�(�) が一致したと
ただし、上式で、( �(��きのみ｢1｣で、一致しないときは｢0｣になる２値数を表すものとする。

5.3“仮想粒子衝突規則”の表現

各格子点位置に到着した個々の仮想粒子は、そこに到着した他の仮想粒子と衝突を起こし、運動量やエネルギーを交換する。この結果、“到着粒子分布”は、別の“ビット列”で表現できる“出発粒子分布”に変換される。このときの具体的な変換内容を、“仮想粒子衝突規則”と呼び、���������� � �(�)����������) によって表す。ξ(�(�)���������� � �(�)����������) は、生じうるすべての入力ビット列のひとつひとつ
ξ(�(�)
������������
�(�) （�＝1～ 2 のいずれか）に対して、衝突後、それをただひとつの出力ビット列 �����������
�(�) （�＝1

— 10 —

SENAC Vol. 49, No. 4（2016. 10）���������� � ����������
～ 2�� のいずれか）に変換する規則を与える。ξ(�(�)�(�)) の値は、変換してよい入力ビット列���������� と出力ビット列 �(�)���������� が与えられたときのみ｢1｣になり、それ以外は｢0｣になる２値数である。
�(�)�(�)（�＝1～ 2�� ）に対して次式が成り立つ。
従って、任意の出力ビット列 �������������́ (�� �� �� �� �) �� ���������
�(�) � ＝ �� �����������������������

������������� � �����������ξ ��(�)�(�)� × � ����������������������������
�(�� �� �� �� �) �� �����������
�(�) ��

・・・②式“仮想粒子衝突規則”は、概念的には、2�� × 2�� の巨大な数値表になる。なお、すべての仮想粒子 ji が衝突を起こさず、素通りする場合も衝突の特別な場合である。この場合、ビット列に変化はないので、数値表の対角位置に｢1｣が並ぶ。また、数値表における｢1｣｢0｣の位置は、確率的に変動してもよい。すなわち、同じ入力ビット列に対して、確率的に異なるビット列を出力する“仮想粒子衝突規則”も可能である。いずれにしても、a,b＝1～ 2�� としたとき、任意の a に対して、
�

���

���

���������� � ����������
ξ(�(�)
�(�)) ＝ 1 ・・・③式 が成立する。

5.4

衝突過程における仮想粒子分布の変化

時刻 � での格子点位置(x,y,z,r)における“出発粒子分布”�����������������������������́ (�� �� �� �� �) の ji 成分を求める。
�́ �� (�� �� �� �� �)＝����〔��� (�) × ( ��������������������������́ (�� �� �� �� �) �� ����������
�(�) )〕

���

＝�

���

〔��� (�) × �

���

���

���

（∵ 前記①式の関係を適用）���������� � ����������
�(�)� × ( ����������������������������
�(�� �� �� �� �) �� �����������
�(�) )�〕�ξ ��(�)・・・④式 （∵ 前記②式を適用）

���������� � ����������
�(�)� × ( ����������������������������
�(�� �� �� �� �) �� �����������
�(�) )�〕�ξ ��(�)－④式から①式を引くことで、衝突過程における仮想粒子 ji の粒子分布の変化を求める。
�́ �� (�� �� �� �� �) － ��� (�� �� �� �� �)

＝�
＝�

���

〔��� (�) × �

���

���
���
���

�

���

���

－ �

＝�

���

���

�

���
���

���

���

�

���

〔��� (�) × ��

���

���

〔��� (�) × ( �������������������������
�(�� �� �� �� �) �� �����������
�(�) )〕

���

���

���

����������������������������
���������� � ������������ �� �� �) �� ������������
�(�) )〕 （∵ 値が 1 の③式を代入しても同じ。
）
ξ(�(�)
�(�))� × ( �(������������������������������
���������� � ����������
���� (�) × ξ ��(�)�(�)� × ( �(���� �� �� �) �� �����������
�(�) )� （∵ 上式の第２項で、�＝�，�＝�, �＝� とおく。
）
〔��� (�) × ��

���

���

�

����������������������������
���������� � ����������
���� (�) × ξ ��(�)�(�)� × ( �(���� �� �� �) �� �����������
�(�) )�

���

���

－ �

＝�

���

���

���

���������� � ����������
�(�) )〕
�(�))� × ( ����������������������������
�(�� �� �� �� �) �� �����������
ξ(�(�)

����������������������������
���������� � ����������
�(��� (�)－��� (�)) × ξ ��(�)�(�)� × ( �(���� �� �� �) �� �����������
�(�) )�

・・・⑤式

������������������������������ �� �� �) �� �����������
�(�) ) の具体的な評価式を考える。通常は、次式を用いる。
ここで、⑤式に含まれる ( �(�����(�� �� �� �� �) �� ���������
�(�) ) ＝ ∏��･�1 � ��� (�� �� �� �� �)�
( ��������������������������〔���� (�� �� �� �� �)�

� (�)����� (�)〕 ・・・⑥式

このとき、右辺に現れる ��� と ��� はいずれも｢1｣または｢0｣の値しかとらない２値数であり、「1� ＝1，1� ≡1，0� ＝0、0� ≡1」であることに注意すれば、右辺は、��� と ��� の全ビットが一致したときのみ｢1｣で、それ以外のときは｢0｣になることがわかる。このような“一致判定”は、２値数を扱う限りは、論理演算式で行う方が簡— 11 —

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —単である。しかし、ここでは、後で集合平均をとるので、��� が確率として｢0 と 1 の間の実数値｣をとった場
合にもそのまま拡張できる評価式を用いたい。このため、⑥式の右辺の式表示を用いる。

5.5“仮想粒子衝突規則”の適用頻度を変えるテシャラの方法“テシャラの方法”を用いる場合、“仮想粒子衝突規則”の入力ビット列を作る前に、“到着粒子分布”に対してある特別な操作を行う。具体的には、“エネルギー ｊ の粒子が存在することを隠すマスク” �� と“エネルギー ｊ の粒子が存在しないことを隠すマスク” �� を導入する。�� と �� は２値数であり、通常は｢1｣で何の効果も発揮しない。マスク効果を発揮させるときのみ、値を「0」にする。シミュレーション計算では、ある一定の確率で「0」の値をとらせる。この場合、⑥式の表示は、次のとおりになる。
���(�� �� �� �� �) �� ���������
�(�) ) ＝∏��･��� (� � ��� (�� �� �� �� �))�( ����������������������
�������� ��� (�� �� �� �� �)�

� (�)����� (�)]・・・⑦式もし �� ＝0 であれば、静止粒子が存在(��� =1)していても、上式右辺で �� ��� =0 なので、入力ビット列ではその存在を検知できない。従って、「静止粒子」＋「速い粒子」⇒「遅い粒子」＋「遅い粒子」 という衝突変化は生じなくなる。また、�� ＝0 であれば、静止粒子が存在していない(��� =0)場合でも、上式右辺で�� (� � ��� )=0 なので、入力ビット列では存在していないことを検知できない。従って、上記と逆の「遅い粒子」＋「遅い粒子」⇒「静止粒子」＋「速い粒子」 という衝突変化は生じなくなる。最終的には、“「�� が 1 に���� 」の比”： �� ≡ ��� /����� が上記の衝突変化のバランスを決定することになる確率��� 」と「�� が 1 になる確率�
なる。すなわち、�� の値を調節して、局所平衡状態における仮想粒子数のエネルギー分布を制御できる。

5.6

局所平衡状態における仮想粒子分布を導出する式

��� (�� �� �� �� �) や �� �� (�� �� �� �� �) は、サブマクロな世界にある４次元 FCHC 格子上の時空間で、目���������� →まぐるしく｢1｣または｢0｣の値で変動している。また、5.3 に述べたとおり、“仮想粒子衝突規則” ξ(�(�)����������
�(�)) も、その数値表における｢1｣や｢0｣の位置が確率的に変化してもよい。ここで、このような確率的に変動する２値数を含む⑤式について、統計力学にいう“集合平均”をとる。詳しくは省略するが、局所平衡状態では「⑤式の集合平均＝０」の式が成り立つとして、これから、マクロな速度がＵであるときの仮想粒子分布を導出する。具体的に⑤式にテシャラの⑦式を代入して、集合平均をとると次式を得る。

0 ＝ �
＝�

���
���

���
���

�

���

���

���

�

�����������������������
����������������������������������
���������� → ����������
�(�)� × ( ����������������������������
�(�) )�
�(��� (�)－��� (�)) × ξ ��(�)�(�� �� �� �� �) �� ������������ (�)����� (�)����������������������
����������) ∏��������������������
���������� → �(�)� ������������������ �� ･���� ･(� � ����� (�) � ��� (�)}ξ(�(�)〕 ・・・⑧式
�� (�� �� �� �� �))��������� ･��� (�� �� �� �� �)�

���

ここで、集合平均をとるとき、これまで２値数として扱っていたいくつかの数は、0 以上 1 以下の実数値をとりうる数としての扱いに変更する必要が生じる。上式では、この数に上付きバーを付して区別した。ま���������� → ����������
����������������������������
�(�) )、�� と ��� 、 �� と (� � ��� ) の確率的変動は互いに独立なので、た、ξ(�(�)�(�)) と ( �(���� �� �� �) �� �����������
この積の集合平均をとると、その結果は、それぞれの集合平均の積になる。
⑧式から長い式変形過程を経て最終的に「マクロな速度がＵm/sec の局所平衡状態に対応する“仮想粒子数のエネルギー分布”」が以下のとおり求まる。ただし、u ≡(Ｕ/ｃ)<< 1 とする。
��������
��� (u)＝�＋����������
�����(�)��� �� ＋��(�)＋�� �(�)＋���� ･��� と定義するが、非熱流体解析用３速さモデルの場合は、�� = �� = 1 として、
ここで、�� ≡ ��� /��� のみを調節パラメータとして用いる。上式から一般に、�� の値を操作して、仮想粒子数のエネ
ルギー分布を変化させられることがわかる。マクロな質量、運動量、運動量フラックス等はすべてこの粒子分布から導出されるので、ナビエ・ストークス方程式等のマクロな支配方程式を正しく満足するように �� の値を調節して、マクロな係数の差異を解消することができる。

— 12 —

６.

SENAC Vol. 49, No. 4（2016. 10）“マクスウェルの達人操作”の展開

6.1 “マクスウェルの達人操作”
分子や原子の動きを個別に自由に操って、負のエントロピーが生じるような過程を生じさせる
想像上の小さな生物？で“マクスウェルの魔物”と言われるものがある。テシャラが ݃＝1 を達成するために用いた方法は、数多く存在する仮想粒子の衝突過程の中から、静止粒子が介在する衝突に注目し、静止粒子の“存在”や“非存在”をマスクによって隠すことによって、意図的にその発生頻度(確率)をある適切な値に制御するものであった。これによって、自然界では生じないようなある特別な“仮想粒子数の局所平衡エネルギー分布”をサブマクロな世界で実現し、それを疎視化したマクロな世界では、݃＝1 となる絶妙なバランスを得ることができる。ここでは、“多少の差異”をピッタリとキャンセルする“達人技”のイメージを強調して、このような個別粒子に対する意図的な操作を“マクスウェルの達人操作”と名付けることにする。
テシャラの方法は、“マクスウェルの達人操作”のひとつであると見なすことができる。本章では、この“マクスウェルの達人操作”のさらなる展開について考える。
6.2 負の粘性の発現
仮想粒子の衝突規則は、通常、すべての効果を合計したとき、“正の粘性”を発現させる[4]。
もし、“負の粘性”を発現させる特別な方法を見つければ、これを“正の粘性”を発現させる通常の衝突規則と確率的に組み合わせて、それを疎視化したときにマクロな世界で発現する“粘性”の値を自由に制御できるであろう。このことは、格子ガス法において、格子点数を増やさなければいけない２つの理由（①解像度を上げたい、②粘性を小さくしたい）のうち後者を解消し、格子点数は、もはや、求める解像度の要求からだけ決めることができるようになる。このことは、結果的に、必要とされる計算機パワーの大幅な節約を可能にする。
Daniel H. Rothman は、1989 年の論文[5]“Negative-Viscosity Lattice Gases”において、既に、“負の粘性”を発現させるひとつの方法を見出し、２次元の流体シミュレーションで粘性が実際に小さくなる効果を確認している。このとき用いられた方法は、仮想粒子に「赤」と「青」の色をつけてその違いに応じた処理を行う“着色法”の応用である。仮想粒子に対して通常の衝突をさせた後に、その格子点近傍の赤色粒子と青色粒子がもつ運動量の差の勾配が強調される向きに出発する仮想粒子に赤色を着色する。ここでは、“マクスウェルの達人”は、個々の粒子の動きをみて色を塗る着色屋さんの仕事をしている。
さて、粘性は、速度勾配があるところで、その効果を現す。例えば、ある短い時間、流線がある場所で局所的にある弧を描いたとする。いま、この流線上に位置する流体塊の重心速度がＵm/sec であり、弧の少し外側の流体塊の重心速度はＵより少し大きい速度Ｕ+εで、弧の少し内側の流体塊の重心速度はＵより少し小さい速度
Ｕ－εであったと仮定する。このとき、速度Ｕの向き(流線の接線
方向プラス向き)に横軸をとり、流線に垂直な方向に縦軸をとって、その近傍にある格子点間で、仮想粒子がどのようにやりとりされているか？を想像すると
図３のようになる。
図３

局所平衡的に速度勾配が維持された場所の仮想粒子分布(相対確率)

— 13 —

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —図３では、議論を簡単にするため、２次元平面上の流れを想定し、仮想粒子の移動可能な向き
は、右上、右下、左上、左下の４つのみを仮定した。注目している流線Ｕ上の真ん中の格子点における、局所平衡状態時の仮想粒子の到着分布及び衝突後の出発分布が、図３中央の黄色点線枠
内に記載されている。すなわち、到着分布では、右上に向かう仮想粒子が左下の格子点から相対確率６の頻度で到着、右下に向かう仮想粒子が左上の格子点から相対確率７の頻度で到着、左上に向かう仮想粒子が右下の格子点から相対確率３の頻度で到着、左下に向かう仮想粒子が右上の格子点から相対確率３の頻度で到着していることを示す。また、これらの仮想粒子が衝突して移動の向きを変えた後が出力分布であり、仮想粒子は右上の格子点に向かって相対確率７の頻度で出発、右下の格子点に向かって相対確率６の頻度で出発、左上の格子点に向かって相対確率２の頻度で出発、左下の格子点に向かって相対確率３の頻度で出発していることを示している。
ここで、a. 速い側から到着する仮想粒子の寄与は、加速の寄与：減速の寄与＝７：２であるから、正身＋５の加速効果。
b. 遅い側から到着する仮想粒子の寄与は、加速の寄与：減速の寄与＝６：３であるから、正身＋３の加速効果。
c. 速い側へ出発する仮想粒子の寄与は、加速の寄与：減速の寄与＝２：７であるから、正身－５の減速効果。
（↑これは、速い側の流体に対しては、正身＋５の加速効果を与えている。
）
d. 遅い側へ出発する仮想粒子の寄与は、加速の寄与：減速の寄与＝３：６であるから、正身－３の減速効果。
（↑これは、遅い側の流体に対しては、正身＋３の加速効果を与えている。
）となり、流線上の運動量は収支がつり合って局所的な平衡状態を保持している。
以上の考察から、速度勾配が維持されている局所平衡状態では、仮想粒子の衝突前後における速度分布は、流線軸に対して反転されていることがわかる。
“負の粘性”が発現している状態とは、速度勾配がどんどん急激化していく変化途上の状態であるから、上記で考察したつり合い状態を破って、例えば、以下のような状態を作り出せば、“負の粘性”が発現するものと考えられる。
とりあえず、a.と b.は同じとして、c.と d.を次の e.と f.に変える。
e. 速い側へ出発する仮想粒子の寄与を、加速の寄与：減速の寄与＝１：８とし、正身－７の減速効果にする。
（↑これは、速い側の流体に対しては、正身＋７の加速効果を与える。）
f. 遅い側へ出発する仮想粒子の寄与を、加速の寄与：減速の寄与＝４：５とし、正身－１の減速効果にする。
（↑これは、遅い側の流体に対しては、正身＋１の加速効果を与える。）
流線を挟んで、速い側の流体に対する正身の加速効果は以前より増大し、遅い側の流体に対する正身の加速効果は以前より減少しているので、より大きな運動量勾配を作りだす方向に平衡状態が変化していくであろう。また、このときの新しい平衡状態に対応して、はじめに仮定した a.
と b.の状態も更新されていくことになる。以上の考察を図示すると、図４のようになる。
負粘性発現のために修正した

出発確率(相対値)

到着確率(相対値)

出発確率(相対値)

７

６
３

２

８
１
流線の向き

２

３
７
図４

６

４

５

負粘性を発現させるための仮想粒子の出発分布（相対確率）

— 14 —

SENAC Vol. 49, No. 4（2016. 10）

すなわち、“負粘性”を発現するには、通常の衝突規則によって得られた各３次元格子点における出発仮想粒子の速度分布(確率)を図４のように凹凸強調の修正をしてやればよいことがわかる。
図５は、筆者らが行った負粘性発現シミュレーションである。前述の“一辺が 2.56m の立方体”の３次元空間(ただし、全方向に周期的境界条件を適用)を想定した。
“負の粘性”の発現によって、静止状態から自発的に渦が発生して拡大していく様子がわかる。筆者らの方法では、「同じ３次元
格子点の異なる４次元格子点間で、それぞれの４次元格子点に存在する仮想粒子群の重心速度を測定し、その差(揺らぎ幅)を拡大する反発力(揺動力)を与える仮想粒子の速度ペアを見つけてそれらを反転させる」という操作を行ったが、これが“マクスウェルの達人”がした仕事である。

(a) 負粘性発現による渦の拡大過程（Z 軸に垂直な断面上）図５(b)ある時点の３次元可視化図
4 次元 FCHC 格子上で負粘性発現衝突規則を適用し疎視化して得たマクロな速度

6.3“マクスウェルの達人”のさらなる進化
Robert P. Bosch, Jr.は、彼の MIT における 1995 年の論文[6]“A Multi-grid Algorithm for
Lattice Gases”において、格子ガス法におけるマルチグリッド法を提案している。これは、CFD
における BCM(Building Cube Method)に似たものであり、物体に近くて解像度を上げたい領域では細かな立方格子を用い、物体から遠く離れた解像度をあまり必要としない領域では、粗い立方多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —— 15 —

格子を用いる方法である。これも、必要とされる計算機パワーを大幅に節約する方法である。両者の格子が接続される領域では、細かい格子の世界のみを動き回る軽い質量の仮想粒子８個と、粗い格子の世界のみを動き回る重い質量の仮想粒子 1 個が相互に変換される。しかし、同じ向きの速度をもつ軽い仮想粒子が８個集まるのを待っている間、質量や運動量が保存しなくなる。この過不足分の保存量を補うため、速度���で正の質量 m をもつ通常の仮想粒子に対して、速度は同�
じ���を有するが、負の質量:－m と負のエネルギー:－ � mv � をもつ“質量反転粒子”を導入している。
ここでの“マクスウェルの達人”の仕事は、２つの異なる格子点間隔の接続領域にやってくる仮
想粒子を観察して、必要に応じた数の“質量反転粒子”を生成または消滅させることである。
このように、格子ガス法における仮想粒子衝突規則には、それを疎視化したときに発現するマクロな挙動を支配する DNA が集約されており、衝突規則に対して“マクスウェルの達人”が行う操作を自由に発想することによって、まだまだ未知の可能性があることを感じさせる。さらに、
今や、“達人”を超える“人工知能（AI）”の時代が到来している。実のところ、仮想粒子の衝突規則に対するいろいろな操作が、それを疎視化したとき、どのようなマクロな挙動になって現れるのか？その因果関係を知ることは非常に難しい。しかし、他方、仮想粒子の衝突規則は、到着粒子の存在分布を表すビット列を、出発粒子の存在分布を表すビット列に対応させるものであるから、簡単にニューラルネットで入出力を扱うことができる。また、このとき、周辺格子点における仮想粒子の存在分布の影響も考慮する必要があれば、そこのビット列を入力に追加するだけでよい。いずれにしても、このとき用いるニューラルネットの基本構造は、流体中にあるすべての格子点において共通のはずである。また、格子ガス法の最大の特徴は、誤差が蓄積しないビッ
ト演算を使うことができるので、どんな激しい変化を示す流体挙動であってもその計算を安定に行い答えを出すことができるという点にある。もし、乱流について、深層学習を行ったらどうなるであろうか？人工知能を手に入れた“マクスウェルの達人”は、乱流状態における仮想粒子衝
突規則の操作方法を事前に知らなくても、乱流の測定データという自然界の教師データから何か
本質的なものを直接学びとるかもしれない。
これからの“AI コンピューティング時代”が楽しみである。

謝辞
本稿で述べた今後の研究展望は、これまで長期にわたり東北大学サイバーサイエンスセン
ターのスーパーコンピュータ(特に、ＳＸシリーズのベクトルコンピュータ)を利用すること
によってはじめて得られたものである。利用にあたっては同センター関係各位のご親切なご指導とご協力をいただき、心から深く感謝する次第である。
参考文献[1] Chistopher M.Teixeira, “Continuum Limit of Lattice Gas Fluid Dynamics”, MIT, 1993
[2] 笠木伸英, 松本洋一郎, 大橋弘忠, 岩波講座,現代工学の基礎, 計算熱流体力学, p108
[3] 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
[4] Michel Hé non, “Viscosity of a Lattice Gas”, Complex Systems 1(1987), pp.763-789, 1987
[5] Daniel H. Rothman, “Negative-Viscosity Lattice Gases”, Statistical Physics, Vol.56,
Nos.3/4(1989), pp.517-524, 1989
[6] Robert P.Bosch,Jr.,“A Multi-grid Algorithm for Lattice Gases”,MIT,1995