― 18 ―SENAC Vol. 52, No. （2 2019. 4）［共同研究成果］

仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性―

整数型格子ボルツマン法における連続的な粘性制御機能の実現に向けて―松岡 浩
一般財団法人高度情報科学技術研究機構
東北大学電気通信研究所（客員）菊池 範子
技術士事務所ＡＩコンピューティングラボ
筆者らは、東北大学サイバーサイエンスセンターとの共同研究公募制度により、平成 27 年度から 29 年度まで「連続感度解析の実現を目指した整数型格子ボルツマン法流体解析手法の開発」を行ってきた。ここでは、本研究を通じて得られた知見のうち、“仮想粒子の並進過程に干渉効果を加味することによって、数値シミュレーションで模擬できる流体粘性の範囲を幅広く調整できる可能性”について述べる。円柱後流の解析事例から定性的ではあるが、大きく変化する不規則な流動が発生し、その変化は流体粘性が小さくなり乱流化したように見えた。

1．はじめに筆者らが本研究で最終的に目指すものは、昨年４月の SENAC 誌[1]に記載したとおり、「流体中を運動する物体が係わる工学システムのものづくり設計において、斬新な設計アイデアの探索に役立つ“リアルタイム応答流体シミュレーション設計ツール”を実現すること」である。
より具体的には、「ものづくり設計のために必要とされる流体シミュレーションの解像度を維持しながら、①流体中または流体周辺に存在する物体や構造物の位置や形状（“設計上の境界条件”）、②工学的に制御できるある特定の場所の流体運動の速さや向き（“運転上の環境条件”）、などを連続的かつ気ままに変化させながら、その効果を実際の物理現象の実時間に近い応答としてビジュアルに観察できる、コンパクトな規模の計算機による超高速流体シミュレーションシステムの構築」を目指している。
この最終目標への道のりは長いが、筆者らは“整数型格子ボルツマン法（多速さ格子ガス法[2]）”を基礎にして、上記のような“リアルタイム応答流体解析”を実現するために必要な新しい機能をひとつひとつ付加しつつある。本稿で述べる新しい機能は、“連続的な粘性制御機能”であるが、現時点では定性的な評価にとどまっている。今後さらなる改良も必要であるため、本稿の表題では、「可能性」という表現に留めている。
なお、本稿の本文は松岡が執筆し、掲載した数値シミュレーション事例に関する計算及び可視化は菊池が分担した。

２．“整数型格子ボルツマン法（多速さ格子ガス法）”の特徴“リアルタイム応答流体解析”を実現していく場合、筆者は、ナビエ･ストークス方程式を差分化して解く従来の王道的な CFD 手法よりも、よりミクロな視点でボルツマン方程式を離散化した格子ボルツマン方程式に従い、かつ、整数型の確率変数を用いて時間発展を追跡していく“整数型格子ボルツマン法”とも呼ぶべき「多速さ格子ガス法[2]」を用いた方が有利だと考えている。
以下、本手法の概要や特徴は、昨年４月の SENAC 誌等[1,3]にも記載したところであるが、本稿における説明の流れ上不可欠な部分だけを繰り返し記載させていただく。
本手法は、流体が存在する空間中に格子を張り、多数の仮想粒子が、その格子点(サイト)上の仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性― 19 ―

みで他の仮想粒子と衝突して進行の向きを変えながら格子点(サイト)間を走行していく様子を平均化(疎視化)して流体の挙動を模擬する方法である。このとき、仮想粒子がもつ質量・運動量・エネルギーが衝突の前後で保存されるような粒子衝突を想定する限り、その挙動は、自然界におけるある条件下の流体挙動とかなり似かよったものになる。しかしながら、仮想粒子は、格子点(サイト)が存在する位置の間しか移動できないので、速度の大きさも向きも離散的な値をとり、自然界の流体分子のように連続的な値をとることができない。この制約によって、格子ガス法が導くマクロな挙動は、連続流体を仮定している CFD が導く流体挙動とは多少異なったものになる。しかしながら、テシャラが考案した方法［4］に従い、異なる速さをもつ仮想粒子どうしの衝突頻度をうまく調節すれば、この“多少の差異”を解消することができる。これにより、“多速さ格子ガス法”は、CFD を代替できる精度をもったシミュレーション手法となりえる。[5]
筆者が、この“整数型格子ボルツマン法（多速さ格子ガス法）”を“リアルタイム応答流体解析”を実現するツールの最有力候補と考える理由は、本手法が次の特徴をもつからである。
① 状態変化を各速度の仮想粒子の存否情報に分割して 1 ビット幅で追跡し超並列計算ができる。
② 時間発展計算を誤差の蓄積がないビット演算で行いどんな激しい流れでも安定に計算できる。
③ 複雑な形状で時間変化する境界条件でも比較的容易にシミュレーション計算に取込みできる。

３．円柱後流の数値シミュレーションによる考察
3.1 仮想粒子が飛び移る“サイト”の位置座標の設定
３次元空間が静止流体で満たされ、その中に無限の長さをもつ円柱が存在していたとする。時刻ステップゼロにおいて、円柱軸に垂直な方向に流体が動きはじめた場合、円柱の後流がどのような挙動を示すかを調べる。
これは典型的な流体解析の問題で幾何学的な形状は２次元であるが、レイノルズ数が少し大きくなると、すぐに３次元構造の流れが発生することが知られている。このため、本質的に３次元計算が必要である。しかし、ここで、格子ガス法によって３次元の流体シミュレーションを行う
場合の特殊事情がある。すなわち、仮想粒子の速度ベクトルから構成されるテンソルの等方性を向上させる必要があり、そのひとつの方法として、筆者らは“４次元面心超立方体(FCHC)格子”を用いることにしている。このため、３次元問題であっても４次元空間において仮想粒子の動きを追跡する必要がある。
そこで、まずは３次元直交座標系（X,Y,Z）を設定する。円柱軸と平行にＺ軸をとり、流体が動き出す向きに＋Ｘ軸をとる。Ｙ軸は、Ｘ軸にもＺ軸にも垂直であるから、結果的に円柱断面のある直径と平行な方向になる。また、上述のとおり、仮想粒子は４次元空間中を飛び移るので、その位置を特定するには、通常の３次元位置座標（X,Y,Z）に加えて第４次元目の位置座標（R）が必要になる。３次元位置座標（X,Y,Z）が変化する範囲は、現実に流体が存在する３次元空間中の位置に対応して定まるが、R 座標については複数の位置座標値をもてばよく、現実に流体が存在する位置の広がりからは何の制約も受けない。このため、ここでは最も簡単に、R の位置座標値は２個の値(±ℓ/2)のみをとるものとする。ここで、ℓ（＞0）は、シミュレーションの時間発展計
算において時刻ステップがひとつ刻まれる時間⊿t の間に仮想粒子が１方向に進む距離である。
仮想粒子のその方向の速度成分が c であれば、ℓ = c ⊿t となる。R = -ℓ/2 の座標位置から出発した仮想粒子が R 方向成分速度+c をもっていれば、次の時刻ステップで仮想粒子は R = +ℓ/2 の座標位置に到着する。これが R = +ℓ/2 の座標位置から出発した仮想粒子であれば、次の時刻ステップで R = +3ℓ/2 の座標位置に到着すべきであるが、“周期座標”を仮定して、R = +3ℓ/2 ではな
く R = -ℓ/2 の座標位置に到着するものとする。同様に、R = -2/ℓの座標位置から出発した仮想粒子が R 方向成分速度-c をもっていれば、次の時刻ステップでは、R = -3ℓ/2 の座標位置に到着すべきであるが、周期座標の仮定により R = +ℓ/2 の座標位置に到着するものとする。
以上の考察も含め、４次元空間中において仮想粒子が飛び移れる場所については、以下のとお― 20 ―SENAC Vol. 52, No. （2 2019. 4）

りまとめることができる。
①４次元空間中に直交座標（X,Y,Z,R）を設定し、一辺が ℓ の“４次元超立方体セル（2×2×2×2）”を、その中心が X,Y,Z の各方向に間隔 2ℓ で並ぶように配置する。R 方向には、R=0 の位置にセルの中心がくるようにして１列の奥行き分だけを並べる。
②すべての仮想粒子は、それがもつ速度に従って４次元超立方体セルの頂点の位置（以下、この位置を“サイト”と呼ぶ。具体的な座標位置は（X±ℓ/2,Y±ℓ/2,Z±ℓ/2,±ℓ/2）である。）を飛び移って並進移動することができる。また、“サイト”においてのみ、他の仮想粒子と衝突することができ、質量、運動量及びエネルギーの保存則を満たしながら、その速度を変化させて散乱することができる。
③各サイトの位置を識別するためには、そのサイトが属する４次元超立方体セルの中心位置（X,Y,Z,0）と、そのセルの各頂点に付された 1～16 の位置識別番号を指定すればよい。
図１に“サイト”の位置を「◎」で示す。「◎」の内側の丸は４次元目の座標位置が R = -ℓ/2、「◎」の外側の丸は４次元目の座標位置が R = +ℓ/2 のサイトに対応している。
内側から 5,13Z軸

X軸

内側から6,14内側から 1,9内側から 2,10
ひとつの4 次元超立方体
セルとサイト番号

内側から 8,16Y軸

内側から 7,15内側から 3,11内側から 4,12Z軸
X軸Y軸
図１．仮想粒子が飛び移ることができる４次元サイトの配置

仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性― 21 ―

3.2 シミュレーション計算の体系規模
今回おこなったシミュレーション計算では、図１におけるセルの配列は、Ｘ方向に 512 個、Ｙ方向に 256 個、Ｚ方向に 128 個である。ひとつのセルには３次元空間で考えると 8 個の頂点があるので、頂点数としては、約１億 3400 万個の規模に対応している。
疎視化については、スナップショットをとるべき時刻ステップがやってくるたびに、ＸＹＺ方向に 4×4×4 個並ぶ 64 個のセルをひと塊にしてその中に存在する仮想粒子の運動量の合計を求め、128×64×32 個の流速ベクトルを導出した。ただし、３次元ベクトルは紙面上で表現しにくいため、図３に後述するシミュレーション結果は、Ｚ軸に垂直なある平断面上における流速ベクトル分布の時間変化を示すこととした。
3.3 仮想粒子の存在状態の記述と並列ビット演算
格子ガス法のシミュレーションでは、通常「ひとつのサイトには、複数個の静止粒子と、ある特定の指定された速度をもつ運動粒子がそれぞれの速度ごとに高々１個まで存在できる。」と仮定する。ここで「特定の指定された速度」というのは、解析の目的に応じて適切なモデルを選定する。今回のシミュレーション計算の目的は、流体粘性を幅広い範囲で連続的に制御できる方法を探ることにあり、定量的な厳密性は追求しない。このため、仮想粒子がもちうる速度としては、もっとも単純な“単速さ面心超立方体(FCHC)格子 24 速度モデル”を採用した。このモデルでは、仮想粒子は次の 24 種類の４次元速度をもつことができる。
(±c,±c,0,0)、(±c,0,±c,0)、(±c,0,0,±c)、(0,±c,±c,0)、(0,±c,0,±c)、(0,0,±c,±c）
仮想粒子は４次元空間内を動くのでその速度成分は４つある。また、仮想粒子の速さは、すべて同じで大きさで、�(±c)� ＋(±c)� ＋0� ＋0� ＝ √2 c である。（c:基準となる仮想粒子の速さ）

各サイトにおける仮想粒子の存在状態は、仮想粒子の速度 D（＝1,･･･,24）ごとに、その粒子が存在する場合は｢1｣、存在しない場合は｢0｣とし、１ビット幅で表現できる。ここでは、並列計算を効率よく行うため、同じ速度 D をもつ仮想粒子について、セルごとに 16 個の頂点サイトにおける仮想粒子の存否情報を 16 ビット幅のビット列情報にまとめる。そして、時間発展計算に必要な演算は、この 16 ビット列を複数個まとめて 32 ビット幅や 64 ビット幅等で一括計算を行う。

3.4 仮想粒子の並進移動過程における干渉効果の加味
時間発展計算の過程において、仮想粒子は、図１に示した４次元空間中のサイトを飛び移りながら移動していく。あるサイトに到着した仮想粒子は、“衝突散乱過程”に入り、同じサイトに到着した別の速度をもつ仮想粒子と衝突してその向きを変える。次に、この仮想粒子がそのサイトからの“出発粒子”になり、次の時刻ステップがくるまでの間、その粒子がもつ速度に従って次のサイト位置に到着するまで“並進移動過程”が続く。
今回の新しい試みは、この“並進移動過程”に、仮想粒子の干渉効果を導入することである。
しかし、できるだけモデルを複雑にしたくはない。できるだけミクロなレベルで要素還元的なしくみを考えて、そこから自己組織化によりマクロな挙動を発現させたい。
さて、仮想粒子の動きをアニメーション的に想像すると、“衝突散乱過程”を経てたまたま同じ向きに出発することになった複数の仮想粒子は、同じセルの各サイトから同じ向きに飛び出してくる。これらの仮想粒子は、相互にある一定の近距離を保ちながら時刻ステップ間隔⊿t のあいだ同じ向きに同じ速さで並進移動していく。このとき、「もし粒子間力が働けば、仮想粒子の並進軌道が曲げられ、到着先のサイトの位置が確率的にずれてしまうことになるだろう」と想像してみる。この想像が現実の物理現象のある側面の特徴を反映できているか否かは別にして、ここでは「同じセルの異なるサイトから同じ速度で出発する仮想粒子の個数分布を観察して、出発粒子の個数分布の凸凹を強調した分布に置き換える操作を行う」こととした。これは、仮想粒子の並進移動過程に仮想粒子間の相互作用による干渉効果を加味するひとつの試みである。
また、できるだけ要素還元的なしくみにするため、図２に示すように、各セルをその中に存在

― 22 ―

SENAC Vol. 52, No. （2 2019. 4）

する６組の平行４枚組超平面に分けて考えることにする。すると、４枚の平行な超平面上にある４つ正方形の各頂点からは合計で最大４個の仮想粒子が同じ速度をもって飛び出すことになる。
さらに、この速度に垂直な方向をもつ正方形対角線に属する２つの頂点に注目し、その４組の２頂点間で仮想粒子をやりとりして出発粒子の個数分布の凸凹を強調することとした。具体的には、相対する４頂点における仮想粒子の存否分布の AND をとって出発粒子の少ない側の４頂点における仮想粒子の存否分布を置き換え、４頂点における仮想粒子の存否分布の OR をとって出発粒子の多い側の４頂点における仮想粒子の存否分布を置き換えることにした。さらに、この置き換えを確率的に行い、その確率の大小が円柱後流の様相にどのような違いをもたらすかを観察した。
４次元セル

６組の平行４枚組超平面とその垂直軸上記の平行４枚組超平面のそれぞれにおいて、左図のように仮想粒子がもつ速度に垂直な方向をもつ正方形対角線に属する２つの頂点に注目し、その４組の２頂点間で仮想粒子をやりとりして出発粒子個数分布の凸凹を強調する。粒子数の多い側の４頂点を OR 操作で得た粒子配置で、少ない側の４頂点を AND 操作で得た粒子配置で置き換える。左図の例では、この操作で粒子数分布が３：２から４：１に変化し、分布の凸凹が強調されている。

図２．４次元セルに含まれる６組の平行４枚組超平面上における出発粒子分布の凸凹強調操作

3.5 時間発展計算における仮定
時間発展計算を開始する時刻ステップ 0 の時点で、各サイトには、そこに存在できる仮想粒子の最大数の 20％の数の仮想粒子がランダムな向きに配置されていると仮定した。この結果、疎視化して得られるマクロな流速はゼロであり、流体は、直方体形状の中で静止している。次に、時刻ステップ１の時点から、＋Ｘ向きの速度をもつ仮想粒子をＸ＝0 の位置から注入していく。すると、時刻ステップが進むにつれて、流体全体が＋Ｘ向きのマクロな速度をもつようになる。このとき、＋Ｘ側の先にある直方体出口においては、出口直前に存在するサイト上の仮想粒子配置を、出口直後に存在するサイトの仮想粒子配置にコピーして、出口におけるマクロな流速の勾配がゼロになるという境界条件を近似的に実現した。また、±Ｙ方向と±Ｚ方向には、周期的境界条件を適用した。そして、この流れの中の入り口に近い位置に、“Ｚ方向の中心軸をもつ無限大の長さの円柱”を置き、その後流に生じる流体挙動を計算した。
3.6 シミュレーション結果の説明
図３に、３次元的に生じている円柱後流の流速変化のうち、Ｚ軸に垂直なある平断面上における流速ベクトル分布の時間変化を示した。はじめ静止していた流体中におかれた円柱の後流には、流体の速度が増すにつれていろいろな変化が生じる。図３では、紙面の都合から、こうした変化がほぼ似た様相の繰り返しになった時点で、その変化の特徴をよく表す連続した 2 枚のスナップショットを抜粋した。（注：スナップショット間には何回の揺らぎがある。）流速ベクトルの色が赤色なら高速、黄色なら中速、緑色なら低速、青色なら静止～超低速であることを表している。
“出発粒子数分布の強調操作”を行う確率をゼロから上げていくと、円柱後流の様相が大きく変化し不規則な渦発生に至る。その変化は、円柱からの層流剥離→再付着→乱流剥離と見られる様相を含み、流体粘性が小さくなって乱流化したように見えた。
これらの計算は、東北大学サイバーサイエンスセンターの SX-ACE32 ノード(128 ベクトルコア)を利用して行い、どのケースも約 39000 時刻ステップの計算を約 55 分で実行している。

仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性― 23 ―

― 24 ―

SENAC Vol. 52, No. （2 2019. 4）

図３．円柱後流の時系列スナップショットの典型的な連続２画面（ある Z 断面の流速分布）
図３(a)出発粒子分布の強調操作なし（連続 2 画面）

図３(b)出発粒子分布の強調操作 1％（連続 2 画面）

仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性

図３(c)出発粒子分布の強調操作 2％（連続 2 画面）

図３(d)出発粒子分布の強調操作 3％（連続 2 画面）― 25 ―

― 26 ―

図３(e)出発粒子の強調操作 4％（連続 2 画面）

図３(f)出発粒子の強調操作 5％（連続 2 画面）

SENAC Vol. 52, No. （2 2019. 4）

仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性

図３(g)出発粒子分布の強調操作 6％（連続 2 画面）

図３(h)出発粒子分布の強調操作 7％（連続 2 画面）― 27 ―