JPH10501356A

JPH10501356A - 物理的過程のシミュレーションにおける粘性の減少

Info

Publication number: JPH10501356A
Application number: JP8501339A
Authority: JP
Inventors: ケニスアール．トラウブ; トーマスフレデリックジュニアナイト; キムモルビグ; クリストファーエム．テイクセイラ
Original assignee: エグザコーポレーション
Priority date: 1994-06-08
Filing date: 1995-06-07
Publication date: 1998-02-03
Also published as: AU2703195A; EP0764306A1; WO1995034041A1; EP0764306A4; US5606517A

Abstract

(57)【要約】物理的過程をシミュレーションする方法（１００）を実現するコンピュータである。この方法は各ボクセルの数に応じた状態ベクトルをメモリ内に記憶する方法からなる。各状態ベクトルは複数の整数からなり、この整数は、実現可能な幾つかの運動量状態のうちの特定の運動量状態に対応し、さらにこの特定の運動量状態を有する要素の数を表す。各整数は２つ以上の実現可能な値を有する。ここで提供する方法は、異なる運動量状態の要素間の相互作用をモデル化するものであり、相互作用演算(１２０〜１３４)、シミュレーションされる物理的過程の粘性を変更させる粘性修正演算（１３６〜１４８）、そして新しいボクセルに要素の運動を反映させる移動演算（１１４）を状態ベクトルに施す。

Description

【発明の詳細な説明】物理的過程のシミュレーションにおける粘性の減少発明の背景本発明は、物理的過程、例えば流体の流れをシミュレーションするためのコンピュータシステムに関する。高いレイノルズ数の流れをシミュレーションする従来の方法は、ナビエ-ストークス（Navier-Stokes）の微分方程式の離散解を求める方法であって、高い精度の浮動小数点演算は巨視的物理量（例えば密度、温度、流速）を表す変数について離散的な空間の各位置で実行するものである。最速かつ演算力のある利用可能なコンピュータを用いられてきたが、未だ非常に制限され厳密でない結果しか得られていない。実行処理を扱いやすくするために、分解能が非常に粗い格子を用いても、粗い分解能では連続した浮動小数点の算術演算を実行する際の固有の丸め誤差が累算されて容認出来ない誤差が解に生じる。一般に格子気体（すなわちセルラ(cellular)）オートマトンとして知られているものについて微分方程式による解析に置換するための努力が長い間なされている。この格子気体オートマトンにおいてはナビエ-ストークス(Navire-Stokes)の方程式を解くことにより得られる巨視的レベルのシミュレーションが、格子上のサイト間を移動する粒子の運動を実行する微視的レベルのモデルによって置換されている。ここにおける目的は、正確な巨視的結果（すなわち、ナビエ-ストークスの方程式によって定められる密度、温度、その他）を導き得る粒子の相互作用及び運動の微視的レベルのモデルを見いだすことであった。これまでの格子気体のシミュレーションは、ビットの短いベクトルによって表される粒子を有する各格子サイトでの限られた数の粒子を仮定している。各ビットは特定の方向に移動する粒子を表す。例えば、ベクトルの１つのビットは、特定の方向に沿って運動する１つの粒子が存在するか（１と設定）存在しないか（０と設定）を表す。このような１つのベクトルは６ビットを有し、例えば、1100 00とは２つの粒子がＸ軸に沿って反対方向に移動し、Ｙ軸及びＺ軸に沿って移動する粒子は存在しないことを表す。衝突の規則の組は、各サイトの粒子間の衝突の挙動を支配する（例えば、Ｘ軸に沿って移動する２つの粒子間の衝突は、Ｙ軸に沿って離れていく２つの粒子を生成することは110000ベクトルは001100ベクトルになることを表す。）。この規則は、ビットの順列を実行(例えば110000を001 100に変換する)する参照用テーブルに状態ベクトルを与えることで実行される。粒子は隣接したサイトに移動する（例えば、Ｙ軸に沿って移動する２つの粒子はＹ軸に沿った右及び左の隣接したサイトに移動する）。モルビグ(Molvig)氏他は改善した格子気体の手法を教える。この手法において、特に、粒子のエネルギー及び運動方向の変化を知るために、より多くのビット（例えば、超音波流に対して５４ビット）が各格子サイトの状態ベクトルに加えられ、さらに、完全な状態ベクトルの部分を表す衝突規則が用いられる。これはモルビグ(Molvig)氏他によるPCT/US91/04930、離散的運動理論、格子気体の動力学及び流体力学の基礎の研究会の予稿集(Proceedings of the Workshop on Discrete Kinetic Theory，Lattice Gas Dynamics，and Found ations of Hydrodynamics World Scientific Publishing Co.，Pte.，Ltd.，Sin gapore(1989))のモルビグ(Molvig)氏他による「３次元格子気体流体における離散的な人工構造の除去(Removing the Discreteness Artifacts in 3D Lattice-G as Fluids)」、及びSpringer Proceedings in Physics，Vol．46，複合物理系のセルオートマトンとモデル化（Cellular Automata and Modeling of Complex Ph ysical Systems)，Springer-Verlag Berlin，Heidelberg（1990）の中のモルビグ(Molvig)氏他による「実在流体の動力学に対する多種格子気体のオートマトン (Multi-species Lattice-Gas Automata for Realistic Fluid Dynamics)」に記載されている（これらの開示内容が本願明細書に組み入れられる)。モルビグ(Mo lvig)氏他による教示されるこれらの改善等は、最初の実用的な格子気体のコンピュータシステムを提示される。流体の流れのモデル化で不正確な初期の格子気体モデルを作った離散的な人工物(artifact)は除去された。離散的な人工物を回避するもう１つの方法は、格子-ボルツマンモデルと呼ばれるもので、格子気体の手法のブール変数を実変数に置き換えるものである。これはチェン(Chen)氏他による「磁性流体力学のシミュレーションに対する格子ボルツマンモデル(Lattice Boltzmann Model for Simulation of Magnetohydrodyn amics)」Physical Review letters，Vol．67，n．27，30 Dec．1991及びキアン( Qian)氏他の「ナビエ-ストークス方程式の格子ＢＧＫモデル(Lattice BGK Model s for Navier-Stokes Equation)」Europhysics Letters，Vol．17，pp．479-484 ，1 February 1992に記載されている。格子における各サイトの各状態に対する個々の粒子の存在を監視するよりむしろ、格子-ボルツマンモデルは各々の上述の状態の粒子の分布関数を監視する。キアン(Qian)氏他の説明によれば、格子−ボルツマンの方法で用いる実数は緩和の方法の適用を許すものであり、これは単純に以下のように記述できる。ここで、Ｎ(ｔ)は時刻ｔでの量、Ｎ_eは平衡状態(ボルツマン)での量であり、ω は０から２の間の値を取る緩和パラメータである。この方法は０＜ω＜１のとき副緩和(subrelaxation)といい、１＜ω＜２のとき過緩和(over relaxatio n)という(ω＝１のとき、Ｎ(ｔ＋１)は単にＮ_eとなり、緩和は起らない。)。キアン(Qian)氏他はωが増加するとシミュレーションした流体の剪断粘性は減少することも記述している。発明の概要本発明は、各状態を表す整数値からなる状態ベクトルによって、各格子がサイト（ボクセル(voxel)）が表される格子気体系の粘性を修正することを特徴とする(例えば、各ボクセル、０〜２５５の要素は特定のエネルギーを有する特定の方向に移動できる。)。系の粘性が減少すると格子の有効密度が増加し、システムがある物理的過程をシミュレーションできる分解能は劇的に増加する。系の粘性を増加させると、高粘性流体のシミュレーションが可能となる。粘性は、剪断力(すなわち、流体の流れの方向と平行な方向に作用する力)に対するの流体の抵抗力の大きさである。実在する流体では、粒子の速度が平均値に近づこうとさせる原因となる流体内の隣接する粒子間の相互作用の結果が粘性である。格子系においては、あるボクセル内の粒子の正味(net)の速度が隣接したボクセル内の粒子の正味(net)の速度に近づこうとさせる傾向が原因となる特定のボクセル内の粒子間の相互作用の結果が粘性である。実在する粒子によって占められる物理空間より十分に大きいシミュレーションされる空間の領域を各ボクセルは表すので、ボクセル間の相互作用の結果である粘性は実在する流体の分子粒子間相互作用の結果よりも大きくなる(すなわち、各ボクセル間の相互作用を「平均」した結果は、各分子粒子間の相互作用の結果よりも大きい空間領域に作用する。)。典型的には、系の粘性は緩和の手法を用いて修正される。キアン(Qian)氏他によって述べられている緩和の手法は、実数を用いる系に注目しているので、整数をベースとしたシステムに適応することはすぐに明確とはならない。実際に、キアン(Qian)氏他によって述べられている整数値からなる状態ベクトルを用いた格子気体系へ緩和手法を直接適用しても解くことはできない。これは実数値の緩和パラメータを整数値に乗算すると１つの整数値ではなく実数値となるからである。乗算によって生成された実数値は整数値に切り捨てられるので、このような切り捨てはこの系の質量、運動量及びエネルギーをもはや保存しない。本発明は、実数値の緩和パラメータによって修正される整数値の状態ベクトルを可能にし、切り捨て誤差を発生させずに質量、運動量及びエネルギーを保存する。格子系の粘性は格子密度を増加させると減少する(すなわち、各ボクセルで表されるシミュレーションされる空間の量(quantity)の減少させることによって) 、さらに過緩和を使用して減少させることもできる。粘性νは緩和パラメータω によって表される。ここで、Ｔは流体の温度である。このようにして、例えば、１つの緩和パラメータ（ν＝Ｔ／２）に関して、１．８(ν＝Ｔ／１８)の緩和パラメータは、格子内の粘性を９分の１に減少させる。本発明は、物理的過程をシミュレーションする格子系の能力を改善することを約束する。格子密度が増加すると同じ効果が生じるために(すなわち、格子系の粘性が減少する)、過緩和を使用することは格子密度を有効的に増加させる。従って、過緩和を使用することは、特定の分解能(若しくは、特定のプロセッサは物理系をシミュレーションできる分解能)を有する物理系をシミュレーションするために必要な処理にめざましい効果を発揮する。例えば、３次元の格子の有効密度を１０倍に増加させると、特定のレベルの分解能で格子を有する物理系をシミュレーションするのに必要な処理をほぼ１００００の１に減少する（すなわち、過緩和を実行するのに必要な追加処理を１０の３乗ほど少なくし、与えられた速度の流体をシミュレーションするのに必要な時間は１０分の１に減らせる。）。本発明は物理的過程をシミュレーションする新しいコンピュータシステムを実現する。１つの要素(例えば、特定のエネルギーを有する特定の方向に移動する多くても１つの要素)が各格子サイト（ボクセル）にある格子気体モデルの代わりに、本発明は、各ボクセルで幾つかの状態の各々の多数の要素が存在する多数要素の手法を用いている（例えば、０〜２５５個の要素が特定の方向に移動できる）。状態ベクトルはビットの組であったがその代わりにある状態での状態数を表す各々の整数の組である（例えば、０から２５５までの範囲の整数を与える８ビットの組である１バイト）。ボクセル内の異なる運動量状態の要素間の相互作用をモデル化するために、コンピュータシステムは状態ベクトルの相互作用演算を実行する。典型的には、これらの相互作用演算はボクセルの状態ベクトルに対する規則の組を用いて各ボクセルに対して実行するもので、この各規則は状態ベクトルから特定の整数の組を修正する。シミュレーションされた物理的過程の粘性を変更させるために、コンピュータシステムは状態ベクトルの粘性修正演算を実行する。これらの演算は、概して相互作用演算の後に実行し、相互作用演算の間に使用する規則と似たものか若しくは同じ規則の組を使用する。同じ規則が使用される場合、相互作用を演算する間の第１の量及び粘性修正演算をする間の第２の量によって、この法則は状態ベクトルを修正し、緩和パラメータによって、第１の量と第２の量が関係づけられている。相互作用演算に用いられている規則は質量、運動量及びエネルギーを保存するため、この方法は、これらの特性が粘性修正演の間で保存されること保証する。第２の量が緩和パラメータから得られる実数を第１の量に乗算することによって決定されるとき、第２の量が整数値であることを保証するために、乗算のこの結果は切り捨てられる。切り捨て演算によって、系に統計的バイアスを導入しないように、０から１までの間の乱数値は切り捨て演算の前に乗算の結果に加算される。１より大きく２未満の値の緩和パラメータを用いると、格子系の粘性は減少する。緩和パラメータが２に近づくにつれて、シミュレーションされた系の粘性が０に近づき、系は不安定になる。動摩擦の形である粘性は系のゆらぎを減衰させる傾向がある。このようにして、これらのゆらぎが抑制されずに系に広がっていくので、粘性がなくなったときには不安定性が発生する。１．８以下の緩和パラメータを用いると、一般的に不安定性を避けることができることが見いだされた。相互作用及び粘性の修正演算を実行した後に、コンピュータシステムは状態ベクトルについて、新しいボクセルに対する要素の移動を反映する移動演算を施す。図面の簡単な説明図１は物理的過程のシミュレーションシステムによる手順のフローチャートである。図２は微小ブロックの斜視図である。図３は過緩和を伴う衝突を実行する手順のフローチャートである。図４はすべり表面の動力学を実行するシステムのブロックダイヤグラムである。図５は正反射の図である。図６は物理的過程のシミュレーションシステムの機能ユニットのブロックダイヤグラムである。図７は図６のシステムの微視的動力学ユニットのブロックダイヤグラムである。図８は図７の微視的動力学ユニットの１つのボクセルデータパスのブロックダイヤグラムである。好ましい実施例の詳細な説明図１を参照すると、手順１００に従って、物理的過程シミュレーションシステムは演算していることが解る。最初に、タイマーが初期化される（ステップ１０２）。次に、格子内の特定のボクセル（若しくは位置）を示す１つのボクセルのカウンターは、格子内の第１のボクセルを示すように初期化される（ステップ１０４）。初期化後、システムはボクセルの数（総数）によって指示されるボクセルに対応する状態ベクトルをロードする（ステップ１０６）。状態ベクトルはボクセルの状態を完全に定義し、４９（ビット）若しくはより多くの多数ビットの入力からなり、その各々は整数値に対応する。これらの４９（ビット）の入力は、静止状態、第１のエネルギー準位における２４方向のベクトル及び第２のエネルギー準位における２４方向のベクトルに対応する。４９（ビット）の入力のみが必要とされるが、実施例は６つの静止状態を与えので、５４（ビット）の入力を用いる。６つの静止状態は十分な数の静止した「スロット」があることを保証するために用いる。もちろん、この同じ効果は、４９（ビット）の入力の実施例で静止状態に対応する単一の入力のビット数を増やすことによっても成し遂げることができる。多数ビットの入力を用いることで、このシステムはボクセル状態を定義する単一ビット入力を用いるシステムの改善したパフォーマンスを提供する。特に、多くの応用に適さないフェルミーディラック統計を単に作り出す異なる単一ビットのシステムと異なり、このシステムはマックスウエル−ボルツマン統計を作り出す。状態ベクトルをロードした後、このシステムは状態ベクトルの全ての内部ボクセルの演算を実行する(ステップ１０８)。内部ボクセルの演算は他のボクセルに関する情報を必要としない演算である。例えば、流体シミュレーションにおいて、内部ボクセル演算がボクセル内の粒子間の衝突を説明することができる。内部ボクセル演算を完全にするために、このシステムはボクセル数を増加させる(ステップ１１０)。新しいボクセル数（総数）が格子内のボクセル数を越えないならば(ステップ１１２)、システムは次のボクセルの状態ベクトルをロードし (ステップ１０６)、処理を続行する。新しいボクセル数（総数）が格子内のボクセル数を越えるならば(ステップ１１２)、システムは内部ボクセル演算を実行する(ステップ１１４)。内部ボクセル演算は１つ以上のボクセルからの情報を必要とする演算である。例えば、流体シミュレーションシステムにおいて、内部ボクセル演算はボクセル間の粒子の運動を説明する。内部ボクセル演算の実行後、システムは時間を増加し(ステップ１１６)、最初のボクセルを示すボクセル数（総数）を再初期化し(ステップ１０４)、処理を続行する。システムの実施例の演算を以下に詳しく記す。明らかにするために、上述のシステムは、逐次演算するように記述してきた。しかし、以下に述べるように、このシステムは、他の格子システムのように、理想的には平行に演算するのに適している。例えば、内部ボクセル演算は、同時に多数のボクセルについて実行することができる。同様に、内部ボクセル演算を必要とする全てのボクセルについての内部ボクセル演算が完了する間、この内部ボクセル演算は他の内部ボクセル演算を同時に実行することができる。 1993年３月１２日出願の米国出願第08/030,573号、1991年７月１２日出願のPC T出願第PCT/US91/04930号、1991年１２月２０日出願の米国出願第07/812,881号、1990年７月１２日出願の米国出願第07/555,754号及び1993年１２月１０日出願の米国出願第08/165,293号の開示内容が本願明細書に組み込まれる。幾つかの計算の演算を記述する前に、各ボクセルに対する基準状態ベクトルからなる簡単なデータ構造を簡単に記述する必要がある。これは計算に必要となる大部分が演算する基本要素である。各格子サイト若しくはボクセル(これらの２つの用語はこの明細書中においては完全に互換性のあるものとして用いられる) は音速以下の単一種類のシミュレーションに対しては５４個の状態からなる。この状態の数は音速に近い流れ若しくは多数種類のシミュレーションに対して拡張される。この明細書において、状態空間は以下の表記で表される。ここで、Ｎ_iは時間ステップｔでの３次元ベクトルｘで示される格子サイトにおける状態ｉの粒子の数を表す。この状態数は各エネルギー準位で可能な速度ベクトルの数によって決定される。この速度ベクトルは４次元空間（ｘ、ｙ、ｚ及びｗ）内の整数の線速度からなる。この４番目の次元ｗは３次元空間に射影されたもので、従って３次元格子内の実際の速度を示すものではない。音速以下の流れに対してはｉの範囲は０から５３の範囲である。各状態は特定のエネルギー準位での異なる速度ベクトルを表す。各状態の速度は、以下のように各々の４次元内の「速さ」を示す。エネルギー準位０の状態は停止した粒子として知られており、これら（粒子）はある次元内で移動するものではなく、例えば、Ｃ _stopped＝（０,０,０,０）である。エネルギー順位１の状態は４次元のうちの２つは±１を、他の２つは０の速度を有する。さらに、エネルギー準位が２の状態は４次元全てが０か、若しくは４次元のうちの１つは±２で他の３つは０の速度を有する。３つのエネルギー準位の可能な順列の全てを作り出すことは４９の可能な状態の全てを与える（１個のエネルギー（準位）０の状態、２４個のエネルギー（準位）１の状態、２４個のエネルギー（準位）２の状態）。加えて、音速以下の流れの状態空間は、４９個の代わりに５４個の全状態数を与え、１に対立するものとして全部で６個の停止した状態を維持する。ボクセルはマイクロブロックと呼ばれる小さい２×２×２の体積にグループ化される。このマイクロブロックはデータ構造と関係するオーバーヘッドを最小化するのと同様に格子サイトの平行処理を最適化するようにまとめられている。マイクロブロック内の格子サイトに対する簡略化した表記は以下のように定義され、この明細書では至る所で用いられている。ここで、ｘはマイクロブロック内の格子サイトの相対位置を表す。マイクロブロックを図２に示す。微視的動力学（内部ボクセル演算）微視的動力学の演算はボクセル内の純粋に発生する物理的相互作用の組である。この演算の種類（クラス）は流体の粒子及び種々のタイプの物体表面間の物理的相互作用を説明するための流体の状態空間の順列を可能にする。通常の衝突通常の衝突は速度及び方向が変化するように互いに衝突する粒子を許す演算である。粒子は速度ベクトルを決定する状態に置かれるので、粒子の速度及び方向の変化は異なる状態にその粒子を移すことでなされる。典型的な衝突は演算は２つの対の入力（入射）状態ベクトル（全部で４つ）及び同様に２つの対の出力（反射）状態ベクトルからなる。基本的な衝突演算は２つの「流入」粒子が衝突し、２つの「流出」粒子に状態が変化する。流入及び流出の対は常に質量、運動量及びエネルギーが保存しなければならない。従って、５４の状態での可能な４組は「適法」な衝突の組であるとは限らない。この基本的衝突演算は、本来、双方向的なものであり、従って、「流入」及び「流出」状態は、局所的な衝突状態の数に依存して起る衝突のときに決定される。２つの対は選択され、局所的な密度に依存しており、１対（流入）は他の（流出）対の元となる。基本的な衝突演算は以下のように記述され、ここで、SignOfは括弧内の演算の符号（±１）のみを返す関数である。このSi gnOf演算子は値が０のときは、＋１を返す。粒子の数はｉ及びｊの状態からｋ及びｌの状態へ散乱し、Ｎ_scattは小さい正の定数δ（デルタ）を、衝突演算Ｃの符号に乗算することで決定される。このδは衝突則のリスト内に状態指数と一緒に記されている。ｉとｊ若しくはｋとｌの状態の対はより大きい積を有し、より小さい積の原状態の粒子対となる。Ｎ_scattが負ならば、状態ｋ及びｌの状態から状態ｉ及びｊの状態に遷移する。全てが４である状態は同じボクセルでの粒子を表す。衝突の全ては特定のサイトに局所的な状態情報にのみに依存する。状態指数ｉ，ｊ，ｋ及びｌは、ｉとｊの各状態の粒子は、ｋとｌの各状態の粒子と同じ全運動量及び全エネルギーをもつように決定される。全て４の指数は４つの異なる状態を表さなければならなず、全て４の状態は同じエネルギー準位でなければならない。通常の衝突の例として、次の初期状態を与える。ｉ，ｊ，ｋ及びｌの状態の選択から見られるように、iとｊの対はｘ次元（方向）は＋２でｙ、ｚ及びｗ（次元（方向））は０である正味(net)の運動量を有する。加えて、２つの粒子はエネルギー準位１である。同様のことがｋとｌの対にも当てはまる。第１のステップの後に、上で述べた衝突処理においてはＣは− １に等しくなる(SignOf[(２５×４０)−(３０×５０)]＝−１)。負号はｋとｌの対がｉとｊの対の元の粒子であることを示す。衝突則に記されているようにδが４であることから、Ｎ_scattは−４と計算される。Ｎ_scattはｉとｊの状態から減じられ、ｋとｌに加えられる。ここで、衝突演算は終了し、４つの新しい出力状態の集合が生成される。ｋとｌの状態の中からｉとｊの状態に同一の粒子数を遷移することにより、質量も保存される。上述した衝突演算においては、ある状態での粒子数のオーバーフロー及びアンダーフローの可能性がある。質量、運動量及びエネルギーの保存はこのシミュレーション環境において最優先され、チェックしていないならば、ある状態のオーバーフローは運動量及びエネルギーと同様に質量の損失にもなることを注意しなければならない。同様に、アンダーフローが発生することも有り得、この場合には質量が生成されたり、消滅しなかったりする。従って、衝突演算は演算する４組の状態の質量を保存することが必要である。衝突を伴うどのような状態においても、演算がオーバーフロー若しくはアンダーフローのどちらかを引き起こすならば、粒子のあらゆる交換を防止することでこれは成し遂げられる。エネルギー交換衝突エネルギー交換衝突は２つの流出粒子が２つの流入粒子よりも異なるエネルギー準位である点を除いて、前の節で述べた通常の衝突と同様に実行される。音速以下の流れに対しては、３つのエネルギー順位である０（停止）、１、及び２があるのみである。エネルギーを保存させるために、可能なエネルギー交換衝突は２つのエネルギー１の粒子からなる１つの対にのみ起り、他の対はエネルギー２の粒子及び停止した粒子からなる。これらの衝突を特別にする点は通常の衝突の最大限の割合(rate)では起らないことである。従って、エネルギー衝突のクラスは特定の割合(rate)でのみ起ることが可能である。この割合は２つの成分で記述され、これは進行（forward）及び逆行（inverse）である。これらの割合(rate)は整数として記述され、以下のような衝突処理で実行される。ここで、Ｒ₁→₂及びＲ₂→₁は異なるエネルギー準位の状態間の粒子の交換を制御する割合(rate)である。エネルギー交換衝突は、２つのエネルギー１の粒子（状態ｉ及びｊ）が１つのエネルギー２の粒子及び１つのエネルギー０の粒子（状態ｋ及びｌ）若しくはこの反対も起るように衝突するとして構成される。通常の衝突と同様に、質量、運動量及びエネルギーは常に保存する。これらの割合(rate) は温度を基にして決定される。エネルギー交換衝突で用いられるこの２つの割合Ｒ₁→₂及びＲ₂→₁は流体の温度を基にして計算される。しかし、流体の温度はシミュレーション、特に熱伝導を含むシミュレーションの全体にわたって一定である必要はない。これらの割合は局所的な温度を反映するためにシミュレーションの間で動的に最新のものにしなければならない。上のこれらの２つの割合は独立して計算されない。２つの割合によって与えられる適切な情報はこれらの比のみである。全体的なRateは比であり、ここで、Rateは次に記すように温度Ｔから計算される。積の項を展開して、次式を得る。音速以下の流れに対して満足する温度の範囲は１／３と２／３との間である。温度Ｔが０．５より小さいときはRateは１より小さく、温度が０．５より大きいときは１より大きい。このRateから２つのエネルギーの割合が決定される。しかし、エネルギー交換の割合は許される精密さの範囲で変倍化されなければならない。オーバーフロー及びアンダーフローの条件に関する同じ関係はこれらの衝突演算にも適用される。これらの条件はエネルギー交換衝突を可能にし、その上、禁止もここでしなければならない。過緩和(over-relaxation)（粘性減衰）図３を参照すると、物理的過程のシミュレーションシステムは手順１２０に従って各流体ボクセルに対する粘性減衰を伴う衝突を実行していることが解る。システムがボクセルを平衡状態（ボルツマン）にするために、手順１２０はボクセルに対する通常衝突の演算を実行する衝突の段階(ステップ１２２〜１３４)及びシステムが粘性減衰を実行する粘性減衰の段階(ステップ１３６〜１４８)からなる。衝突の段階の初期で、システムは規則カウンターを１つの値に初期化する(ステップ１２２)。次に、規則カウンターに対応する規則によって状態ベクトルの状態に影響するシステムが決定し、規則は４つの状態に影響することを仮定して、システムは変数ｉ，ｊ，ｋ及びｌにこれらの同一状態をロードする(ステップ１２４)。(しかし、以下の議論は各規則が４つの状態に正確に影響することを仮定する。システムは状態の異なる数に影響する規則に対応する。)ｉ，ｊ，ｋ及びｌをロードした後、上で議論した通常の衝突演算若しくはエネルギー交換衝突演算のどちらかを用いて、システムはＮ_i、Ｎ_j、Ｎ_k及びＮ_lの関数として散乱した粒子の数Ｎ_scattを決定する（ステップ１２６）。このシステムはＮ_i及びＮ_jからＮ_scattを引き、Ｎ_k及びＮ_lにＮ_scattを加える(ステップ１２８)(若しくはアンダーフロー又はオーバーフローを避けるため元の変化が省かれるならば０となる)。上で議論したように、ある影響された状態でアンダーフロー若しくはオーバーフローの条件が粒子の交換の原因となるならば、システムは粒子の交換を実行しない。次に、システムは規則カウンター(ｒ)に対応する配列入力(Ｎ_scatt[ｒ])にＮ_scattを記憶させ (ステップ１３０)、規則カウンターを増やす(ステップ１３２)。増加した規則カウンターが規則の最大数を越えないならば(ステップ１３４)、この最大数は好ましくは２７６であり、システムは増加した規則カウンターに対応する規則に基づいたｉ，ｊ，ｋ及びｌに値をロードし(ステップ１２４)、上述で述べた衝突の段階を続行する。増加した規則カウンターが規則の最大数を越えたならば、システムは粘性減衰の段階に入る。最初に、システムは１つの値に規則カウンタを再初期化する(ステップ１３６)。次に、システムは粘性減衰によって散乱する粒子の数であるＮ_V _R を決定する。ここで、ωは緩和パラメータで、Ｎ_scatt[ｒ]は規則カウンタに対応する配列入力で、「noise」は０と１の間の乱数である（ステップ１３８）。「noise」の使用することで力Ｎ_VRを整数にする値切り捨ての演算が特定の方向にＮ_VRを統計学的に偏らせないことを保証する。Ｎ_VRを決定した後、システムは変数ｉ，ｊ，ｋ及びｌに規則カウンターに対応する規則に作用する状態の識別情報（identity）をロードする(ステップ１４０)。システムはＮ_i及びＮ_jからＮ_VRを引きＮ_scatt をＮ_k及びＮ_lに加える(ステップ１４２)。衝突演算に関して述べたように、この交換は影響を及ぼしている状態においてアンダーフローもしくはオーバーフローの条件となるならば、システムは粒子の交換を実行しない。次に、システムは規則カウンタを増やし(ステップ１４４)、この増加した規則カウンタが規則の最大数を越えないならば(ステップ１４６)、増加した規則カウンタに対応する規則に対するＮ_VRを決定し(ステップ１３８)、上述のように粘性減衰のステージと共に続行する。別な方法で、システムは衝突処理を終了する(ステップ１４８)。非すべりの表面動力学上で述べた全ての衝突演算は境界条件によってかなり複雑になる。物体の表面がボクセルを占めるとき、全ての通常衝突及びエネルギー交換衝突の演算は特定のサイトにおける衝突を防止しなければならない。そこでは衝突演算の特別な形が起こらねばならない。これらの演算は粒子と他の粒子が衝突するのではなく、その代わり物体の表面での流入粒子と衝突する。運動量の一部は物体の表面に分け与えるので、これらの衝突は流体の運動量を保存させない。境界衝突がどのように起こるかで決まる表面のタイプの異なる幾つかのクラスがある。非すべり表面での衝突である第１のクラスはこの節で詳細に説明し、他の表面クラスは以下で述べる。非すべり表面は流入粒子が入射してきた方向とは反対方向に同じ速さで表面で跳ね返る反射(bounce back)という特質をもつ。非すべり表面衝突のこのタイプは「反射」としても知られている。反射衝突の効果は物体の表面で静止する流体を持ち込むことである。ボクセルが非すべり物体の表面上にあるならば、ここで、ｉは５４個の状態の１つであり、―ｉは状態ｉと反対方向に同じ速さを有する状態を表す。粒子速度の状態空間は対称であるので、あらゆる状態ｉに対して反対方向で同じ速さの状態である−ｉが存在する。例えば、状態（１,０, −１,０）がｉであるなら、−ｉ状態は（−１,０,１,０）である。 −ｉ状態のオーバーフローは起こりえない、反射が起こった後、そのボクセルでの粒子は新しい外側の方向に応じて移動するので−ｉ状態のオーバーフローは起こりえず、このことは次の時間ステップでこれらのボクセルは空になることを保証する。音速以下の流れに対する５４個の状態のうち３次元空間内で実際に移動するのは４６状態のみである。これはｗ次元での速度成分のみを有する２つのエネルギー準位２の粒子(０,０,０,２)及び(０,０,０,−２)に６つの静止した粒子を加えることによるものである。対称性により、移動状態の半分のみ（２３（の状態））がある表面についての内部状態（内側への状態）とみなされる。反対状態である他の２３（の状態）は外側に向いた状態である。音速以下の流れに対して、状態−ｉは２３の状態のばらばらになった(disjoint)組の１部分であるので、状態ｉは５４の全状態のうちの１部分の２３状態である。すべり表面の動力学高いレイノルズ数で、物体の表面での流れは物体表面での複雑な流れの薄い層の領域である乱流の境界層を形成する。高レイノルズ数の流れをシミュレーションするために、この層を完全に分解する必要がある。しかし，本当にこれをするための分解能の必要条件は天文学的なものとなる。すべり表面は境界層の外側表面を表すのに用いられる手法であり、これは乱流境界層を分解することを避けるシミュレーションを可能にする。図４を参照すると、すべり表面に適用する微視的動力学は手順１５０によって用いられていることが解る。手順１５０のステップは以下で議論する。折り返し防止(アンチエイリアシング)(ステップ１５２) 折り返し防止とは格子と結び付いていない表面に対するより大きな分解能を与えるためのすべり表面を用いる手法である。すべり表面に対して、ボクセルを増加して表面を分離化する(discretizing)ことは境界層をモデル化するのに十分な、正確さを与えない。非すべり表面では全ての流入粒子は反対方向に反射されるので、折り返し防止は非すべり表面を必要としない。しかし、すべり表面にぶつかる粒子は表面で反射されずに、表面上(表面線(surface line))ですべる(conti nue down)。表面の勾配はぎざぎざのある多くの段差からなる階段状である(slop in the surface)ので、接線方向に移動する粒子はこれらのぎざぎざに突き当たる。折り返し防止の手法とは２つの部分つまり流体部分及び表面部分にボクセルを分けることを可能にすることを用いる。２つの部分の折り返し防止されたボクセルの分布は幾何学によって決定され、物体に費やされるボクセルの割合に依存する。折り返し防止の処理は一時的に２つの分離したボクセルつまり流体ボクセル及び表面ボクセルに折り返し防止されたボクセルに分ける。これらの一時的なボクセルの各々に置かれた粒子の数は流体であるボクセルの割合を基にしたものである。計算を以下に示す。流体の割合であるＰ_fは因子Ｐ_scaleによって整数に変倍されて０と１の間の実数である。Randは同じＰ_factorによって変倍されて０と１の間の乱数である。これは完全な表現が再規格化されFloor関数(operation)によって負の無限大に端数を丸められるときに失われる精度を戻すために用いられる。端数の変倍は整数領域での表現に変換する必要がある。折り返し防止をされたボクセルのあらゆる粒子の状態ｉに対して、２つの分離した状態が上述の処理によって生成される。これらの２つの独立した状態は流体及びすべり-表面のボクセルとして分離して扱われる。流体部分は前に述べたように衝突する流体ボクセルとして処理される。表面部分は残りのサブセクションで述べるように表面ボクセルとして厳密に処理される。全ての微視的動力学ルーチンが用いられた後、２つの独立した状態は１つに再び組み合わせられる。反射(Bounce-Back)の割合(ステップ１５４) すべり表面は鏡面のような表面における反射する流入粒子の割合を設定することでもできる。前に述べたように表面が非すべりであるならば、すべらない流入粒子は反射する。反射の割合を設定することは１００％のすべりを与える非常に低い皮膜動摩擦の代わりに、皮膜動摩擦の範囲で特長づけられるすべり表面を可能にする。流入粒子の反射の割合のメカニズムは以下の式で記述される。反射の割合であるＰ_bbは、因子Ｐ_scaleで整数に変倍された０と１の間の実数である。Randは同じＰ_scale因子によって変倍された０と１の間の乱数である。これはこの表現がFloor関数(演算)によって負の無限大に再規格化され端数を丸められるとき失う精度を戻すために用いられる。Ｎ_bbは状態ｉから状態−ｉに反射される粒子数である。非すべり表面の微視的動力学の節で述べたように、Ｎ_-iは状態ｉの反射状態の数である。Ｎ_iは状態ｉの状態数であり、ｉは２３個の流入状態の組の一部である。実際の２３個の状態の組は表面の垂線の方向によって決定される。Ｎ_iは鏡面上での反射のようにすべり表面で反射する粒子のみが含まれる。散乱(ステップ156) の出射角度で表面で跳ね返る特質をすべり表面は有することが解る。すべり表面の目的は流入粒子の運動量の接線成分を保つことである。流入粒子の垂直な運動量は保たれる、すなわち、元の垂直方向の運動量の大きさの２倍分変化する。この過程は鏡面反射としても知られている。粒子の状態空間は離散的であるので、ある表面の垂線は連続的な尺度で変動する必要がある。従って、正確な垂線を基準にして正しい角度で流出する粒子を規定するある１つの状態はないといってもよい。この理由のため、正確な角度となる本当の運動量を得るために、多数の流出状態の重み付けした平均を用いる必要がある。流入粒子を正確に反射させるために、現在実行している散乱は３つから６つの外部状態(out-state)を用いる必要がある。音速以下の流れに対して、２３個の流出状態の組に散乱される必要のある２３個の流入状態がある。各々の流入状態ｉに対して、以下の演算の組は流出状態ｊの各々に対して実行されなければならない。２３個の流入状態ｉのどれか１つに対して３から６のｊの状態がある。ここで、ｉは流入粒子の状態であり、ｊは流入粒子の割合Ｐｊ(＊)は変倍化されている流出状態の１つである。散乱の割合は因子Ｐ_scaleによって変倍化された０と１（を含む）の間の実数である。Randは同じＰ_scale因子によって変倍化された０と１の間の乱数である。流入粒子数（総数）Ｎ_iは(上式のステップ２で示したように)流出状態の１つに各散乱で減ぜられるので、次の出力状態に適用する割合はこれを説明しなければならない。例えば、１つの流入状態が３つの流出状態に等しく散乱されるならば、Ｎ_iに用られる割合は１／３未満である。その代わりに１／３、１／２そして１になる。このようにして、すべての流入粒子は完全に散乱されることを保証する。上で述べた散乱処理の最後の２つのステップはあらゆる状態ｊに対して実行される必要はない。これらの演算は実行される特定の散乱則に依存する。３番目の演算は超過エネルギーの量Exsが散乱処理にわたって累算されることによるものである。多くの散乱法則はあいにくエネルギーを保存しない。後で修正されるので、散乱処理の間で生成されるこの超過エネルギーは累算される。全てのボクセルに対して累算されるので、Exsは下付き文字のｊを含まないことを覚えておかねばならない。散乱処理の最後の演算は通常の外部状態(out-state)に加えて共同状態(co-sta te)の使用を規定する。幾つかの場合、外部状態と共同状態の結合は望むべき運動量を生成する必要がある。この式に記載されているように、共同状態の使用は質量の超過を生成する。２倍もの多くの粒子が内部状態(in-state)から減じられるように外部状態及び共同状態に置かれる。共同状態の数の増加に加えて、超過の質量カウンタMxsは生成される超過の粒子の数によって増加する。後ですべり処理において修正されるので、散乱処理の結果として生成される超過の質量は累算される。超過のエネルギーカウントを有するように、この値は全てのボクセルに対して累算される。外部状態及び共同状態の組は２３個の内部状態に重複する。各２３の内部状態を散乱するために用いられる最小の３つの外部状態があるので、この重複は起る。これにより散乱処理の間に幾つかの外部状態若しくは共同状態で起るオーバーフローの可能性が現れる。オーバーフローが起るならば、用いることができる２つの可能な戦略がある。最初に、オーバーフローの条件をもつ表面のボクセルに対する微視的動力学処理は非すべり表面に逆戻りさせることができる。この場合、その時間-ステップの間に散乱処理によってなされる仕事の全てを消す必要がある。全ての流入粒子は非すべり表面の微視的動力学によって反射される場合に、ボクセルは非すべり表面として仮われる。反射に逆戻りすることは質量及びエネルギーを保存する系を可能にするが、あいにくすべり表面に予想される正確な運動量は保存されない。第２の方法は、オーバーフローが状態カウントをリセットさせる代わりに、状態粒子カウントを最大値で保留する(停止する、止める(clamp))ことである。この場合、質量、エネルギー及び運動量は保存されないが、反射の可能性に逆戻りするよりも、結果はすべりの表面の特性にかなり近づく。押込み/引出し(pushing/pulling)(ステップ１５７) すべり表面に対する結果を処理する次のステップは「押込み及び引出し」と呼ばれるものである。流入粒子が表面で反射された後、流入粒子の接線方向の運動量は保存するので、この散乱処理は慎重に構成しなければならない。しかし、散乱処理は正確な垂直方向の運動量を生成することを保証しない。垂直方向に流体を押込み若しくは引出しすることによって垂直方向の運動量を正確にするために用いられる。垂直方向の運動量が小さければ流体を押込み、大きければ引き出す。この手順は２つの部分に分解する。最初の段階で、垂直方向の運動量の誤差を決定せねばならない。第２の段階は、第１の部分で計算された量によって流体が実際に押込まれたり引出されたりする微視的動力学である。垂直方向の運動量の過剰分(若しくは不足分)は、散乱する前の垂直方向の運動量を比較することによって決定される。散乱後の垂直方向の運動量は散乱前の垂直方向の運動量と等しいが反対方向である。垂直方向の過剰NSの計算は以下に記す。びafterの識別子は散乱処理が起る前後の垂直方向の運動量であることを示す。o ut及びinの識別子は内部へ(表面の内部へ)及び外側へ(表面から離れて)垂直方向に向いた運動量であることを示す。上に述べたように、垂直の過剰分は運動量の前後の差である。しかし、散乱後では、全て外側に散乱された状態になるので、全ての内側に向いた状態は空であることが知られている。従って、散乱後では内側及び外側に向いた垂直の運動量は０になることが知られており、上式から省略される。内側及び外側の垂直の運動量は反対方向であり、従って、反対の符号を有する。このため、差を得るために、内側に向いた垂直方向の運動量は外側に向いた垂直方向の最終的な運動量に加えられる。３つの成分は以下のようにして計算される。れている。状態速度ベクトルはｃ _iで表されている。これは４つの成分をもち、その各々は状態空間の４次元の１つである。前に述べたように、格子気体のアルゴリズムは各々の格子方向の各々での整数の速さで単に移動する粒子に制限する。押込み/引出し処理は現象をかなり複雑する。この処理は上で計算された垂直方向の過剰分をより多く取り除くために設定する。この目的は垂直方向の運動量のみ影響されるように、確かな外部状態(out-state)から他の外部状態(out-stat e)に粒子を遷移することによって成し遂げることである。これらの外部状態(out -state)の対は前もって決定されたリストに含まれる。外部状態(out-state)の対のリストの３つの組がある。粒子は正確な割合の３つのリスト間の法則に従って遷移しなければならない。この割合は表面のこの処理は以下の擬コード(pseudocode)によって記述される。ここで、押込み/引出し則(これは２つの外部状態(out-state)の対である)は１つ以上の粒子を原状態である外部状態(out-state)の１つから到達状態である他の外部状態(out-state)に遷移させる原因となる。空の原状態若しくは満杯な到達状態の数の故に、ある規則が妥当でないならば、規則のリストから他の規則が選択される。リスト内の規則が妥当でないならば、押込み/引出し処理は垂直方向の運動量を完全に修正できない。３つのリストの各々の規則が妥当である総数はＰ₀、Ｐ₁及びＰ₂のあらかじめ決められた重みによって制御される。押込み/引出し規則は状態の対を記述し、原状態及び到達状態を特定する。２つの状態が同じエネルギー準位でない場合、これは頻繁に起るが、規則が妥当であることは運動量を変化させるだけだはなくエネルギーをも変化させる。押込み /引出し規則を妥当なものにするのは以下の処理によって成される。ここで、δは付号のある整数で、処理が押込みであるならば正で、処理が引出しならば負である。概して、δは±１であるが、しかし、垂直の過剰分Ｎ_sが大きい場合には、垂直方向の過剰分の運動量より早く修正するために増加する。Ｎ_sは原状態の総数であり、Ｎ_dは到達状態の総数である。Ｅ_n は規則の数ｎに対して原状態から到達状態に１つの粒子の遷移するのに起因するエネルギーの内部変化である。全エネルギーの変化は散乱処理で用いられた同じ過剰エネルギーカウンターExsを累算する。押込み/引出し規則はエネルギー差 (delta)Ｅ_nと同様に状態ｓ及びｄを規定する。Ｎ_s及びＮ_dで規定される２つの状態の総数がオーバーフロー若しくはアンダーフローするならば、規則が妥当であることから完全に保存しなければならない。冷却(cooling)(ステップ158) 冷却処理は散乱及び押込み/引出し処理の両者が生成した過剰エネルギーを取り除く処理である。冷却処理はボクセルでのエネルギー準位を減ずる幾つかの冷却規則を何度も選択することによって、全ての過剰エネルギーExsを取り除こうとするものである。冷却規則は２つの原状態及び到達状態の４つの状態を規定する。２つの原状態は到達状態と同じ運動量を有するが、２つの原状態のエネルギーの合計は２つの到達状態の合計よりも高いエネルギーを有する。両方の原状態から両方の到達状態に遷移する粒子は、この処理での質量及び運動量は保存される一方、エネルギーは減少する。冷却規則が可逆ならば加熱(heating)もなされることをここでのべなければならない。冷却規則を妥当なものにすることは以下に示すように為される。ここで、Ｎ_s1及びＮ_s2は２つの原状態の総数であり、Ｎ_d1及びＮ_d2は２つの到達状態の総数である。この状態ｓ１、ｓ２、ｄ１及びｄ２は冷却規則によって規定される。冷却規則によって取り除かれるエネルギーの総計は、２つのグループ間のエネルギー差Enを掛けた原状態から到達状態に移動される粒子の数δに等しい。このエネルギー差はあらゆる規則ｎを明確にする。冷却規則は過剰エネルギーExsが取り除かれるまで冷却規則のリストから処理される。リストに与えられている冷却規則のどれも妥当でなければ、冷却処理は目標に到達することができずに終了する。冷却規則に明記された４つの状態のどれもが状態カウントのアンダーフロー若しくはオーバーフローを発生させない。冷却規則が選択された結果として４つの状態のどれもがオーバーフロー若しくはアンダーフローを起こすのであれば、規則は選択されない。ダイエッティング(Dieting)(ステップ１６０) ダイエッティング処理は上で述べた冷却処理と非常に似ている。質量を減少させるダイエッティング規則が多く選択されることによって、散乱処理の間に累積された過剰質量mxsを取り除く。ダイエッティング規則は２つの原状態及び１つの到達状態を明示する。２つの状態のエネルギー及び運動量の和は、１つの到達状態と同じエネルギー及び運動量となる。エネルギー及び運動量が保存される一方、粒子を各々の原状態から取り除き、さらに粒子を１つの到達状態にすることはボクセルの質量を減少させる原因となる。以下にダイエッティング規則の選択の手法を記す。ここで、Ｎ_s1及びＮ_s2は２つの原状態の総数であり、Ｎ_dは到達状態の総数である。状態ｓ１、ｓ２及びｄはダイエッティング規則によって明示される。ダイエッティング規則によって取り除かれる質量の量は原状態から到達状態に移動した粒子数δに等しい。過剰質量Mxsが取り除かれるまで、ダイエッティング規則はダイエッティング規則のリストから処理される。どのダイエッティング規則も妥当でないならば、目標に到達することができずにダイエッティング処理は終了する。ダイエッティング規則に明示される３つの状態のどれも状態カウントをアンダーフロー若しくはオーバーフローを発生させない。ダイエッティング規則が選択された結果として、これらのどの状態もオーバーフロー若しくはアンダーフローを起こすのであれば、規則は選択されない。表面衝突(ステップ162) 表面ボクセルを処理する最後のステップは制限された衝突の規則の組を実行するものである。この表面衝突の処理は流体ボクセルに対して記述したものと同じである。１つの例外は、衝突を許されている状態の組が出て行く状態のみに制限されていることである。移流(中間・相互(Inter)ボクセル演算) 粒子の微視的動力学に加えて、直線で囲まれた３次元格子に沿って、粒子も移動(移流)する。各々の分離した状態は、３次元のｘ,ｙ及びｚの各々(の方向)での整数の速さを有する格子に沿って移動する粒子を表す。整数の速さは０、±１及び±２をからなる。速さの符号は軸に沿って移動する粒子の方向を示す。これらの線速は０から４までの範囲のエネルギー準位を有する粒子に対応する。０から２までのエネルギー準位のみが音速以下の流れに必要であり、全てが５のときは音速に近い流れのシミュレーションに必要となる。移動演算は計算上非常に単純である。各時間ステップで、状態の全体の総数は現在のボクセルから到達するボクセルに移動しなければならない。同時間に、到達するサイトでの粒子は現在の位置からそれ自身の到達するサイトに移動している。例えば、＋１ｘ及び＋１ｙ方向(１,１,０,０)に移動している１つのエネルギー準位が１である粒子は、現在の格子サイトからＸ方向に＋１、ｙ方向に＋１だけ移動しなければならない。この粒子は(１,１,０,０)の移動の前に同じ状態の到達格子サイトで終わる。微視的動力学の次の段階が他の粒子及び表面の局所的な相互作用を元にした状態に対する粒子カウントを変化させ得る。そうでなければ、全ての粒子は、同じ速さ及び方向で格子をたどって移動し続ける。与えられた格子サイトに粒子を配置する全ての移動演算は、ボクセルが再び衝突する前に起らなければならない。全ての微視的動力学演算は局所的な１つのボクセルに実行され、微視的動力学演算の処理が開始される前に全て完了する必要はない。静止した状態は移動状態にならないことをここで記すべきである。加えて、粒子は３次元内で移動状態になるのみである。ｗ次元の０でない状態の値は移動されるべき格子サイトに作用しない。例えば、−１ｚで＋１ｗの状態ベクトル(０, ０,−１,１)である１つのエネルギー準位が１である粒子は、同じ状態ベクトル( ０,０,−１,１)でサイト−１の粒子として移動しなければならない。全く移動しない(０,０,０,２)及び(０, ０,０,−２)の２つのエネルギー準位が２である状態もある。図６を参照すると、上述したシステムは、メモリバンク１４及び規則記憶１６と結合したアプリケーション特有の集積回路(「ＡＳＩＣ」)１２からなる機能ユニット１０を用いて実行できることが解る。システムは単一の機能ユニット１０で実行されるが、ＡＳＩＣ１２は、数百若しくは数千もの機能ユニット１０が共に接続され、パフォーマンスを改善するために平行に演算するように設計されている。ＡＳＩＣ１２は、命令セットを少なくしたコンピュータ(「ＲＩＳＣ」)プロセッサ２０、命令キャッシュ２２及びＲＩＳＣインターフェース論理をそれ自身からなるプロセッサ制御ユニット（「ＰＣＵ」）１８からなる。プロセッサ制御ユニット１８は中心となるコントローラとしての役割を果し、ベクトル的な命令を幾つかの機能ユニットに発する。この機能ユニットはコプロセッサとしての役割を果し、プロセッサ制御ユニット１８からの命令を受け取り、タスクを完了するために命令の定義済みのシーケンスの組を実行する。プロセッサ制御ユニットの命令は、メモリバンク１４にダウンロードされ、そこからプロセッサ制御ユニット１８によって実行される。これらの命令はプロセッサ制御ユニット１８を制御し、プロセッサ制御ユニット１８はＡＳＩＣ１２内の他の機能ユニットに命令を発する。ＰＣＵ１８に対して命令を記憶することに加えて、メモリバンク１４は関係づけられたデータ構造とオーバーヘッドの情報とを一緒にＡＳＩＣ１２に割り当て各ボクセルの粒子の状態を記憶する。メモリバンク１４によって与えられる記憶量はＡＳＩＣ１２が全シミュレーション体積の１部分として割り当てることできる分解能の量を必要とする。パイプラインアクセスモードで高いデータ幅の容量のため、標準的な非同期のＤＲＡＭが選ばれてきたが、メモリバンク１４は標準的な商品であるシンクロナスＤＲＡＭの組から構成されている。微視的動力学の制御論理２６及び微視的動力学のデータの経路２８を含んでいる微視的動力学ユニット(「ＭＤＵ」)は、粒子の動力学を処理する機能ユニットである。これは粒子間の互いの相互作用も表面を有する粒子からもなる。微視的データ経路２８はボクセルデータの計算を実行するために必要なハードウエアからなる。微視的動力学制御論理(回路)２６は、データ経路を設定し、規則記憶装置１６から制限される規則のセットに基づいたデータを選択する。規則記憶装置１６は、速いアクセス時間及び変更される微視的動力学法則をメモリにロードする能力を提供するためにＳＲＡＭから構成されている。マルチーポートＲＡＭ３０は微視的動力学データ経路２８が現在処理している粒子の状態及び関連するデータのレジスタファイルとしての役割を果す。６ポートのスタティックＲＡＭ３０は、微視的動力学データ経路２８に対して２５６つまり６４ビットの入力及び２つの読み込みポートと２つの書き込みポートを有する。加えて、ＲＡＭ３０は、マイクロブロックユニット３２に対する読み込みポート及び書き込みポートを提供する。ポートのこの３つのセットは、統計（処理）と同様に、ＲＡＭ３０及びメモリアクセスユニット３４の間で同時に発生するデータの移動に用いられる。メモリアクセスユニット３４は、ＡＳＩＣ１２内で粒子の元と到達先のデータの移動の中心であり、メモリバンク１４を制御する。図７及び８を参照すると、１つのボクセルを処理するために必要な８ビットデータ経路２９と同一の８個のコピーを提供することによって、平行に微視的動力学データ経路２８は８個のボクセルを処理していることが理解できる。微視的動力学制御論理回路２６はSIMD(単一命令、多データ)コントローラとして構成され、全ての８個のサイトは独立したボクセルデータに同じ演算を実行するために設定する。バスインターフェースユニット（「ＢＩＵ」）３６は、例えば、システム制御に用いられる「ホスト」プロセッサの一般的な目的のように、外部プロセッサにＡＳＩＣ１２に接続される。ＢＩＵ３６の主な機能は、ＡＳＩＣ１２に初期化情報を与え、さらにＡＳＩＣ１２からホストプロセッサに統計的な情報を送るためにシステムのホストプロセッサのために経路を提供することである。加えて、ＢＩＵ３６は、メモリバンク１４がホストプロセッサにアクセス可能であるので、メモリアクセスユニット３４に経路を提供する。移流ユニット(「ＡＵ」)３８は、データ経路を提供し、立方格子をたどって、粒子が移流する（移動する）ことを可能にするように制御し、もう１つのＡＳＩＣ１２が持っている到達ボクセルに移動するこれらの粒子の保持もする。もう１つのＡＳＩＣ１２が持っているボクセルに粒子を送り渡すために、ＡＵ３８は通信ポートをもっている。マイクロブロックユニット３２はＭＤＵで処理されたデータの統計計算するための責任を果す。ボクセルデータがＭＤＵ内で処理された後、マイクロブロックユニット３２は統計（データ）を生成するためにボクセルデータを用いる。このデータはメモリアクセスユニット３４を介してメモリバンク１４に送り戻される。ここでの統計とは表示及び解析する目的でシミュレーションを実行することによりデータを取り出す方法である。以下で詳細にある瞬間的な値の計算に関しての全ての統計計算を記す。しかし、全てのこれらの値は多数の時間-ステップについて累算もされる。累算できる時間-ステップの正確な数は、値が保持できるまでの正確さに依存する。システムが直接サポートする４つの統計量のタイプは、状態質量、質量、エネルギー及び運動量である。以下に示すように、これらの統計量はマイクロブロックを基にして計算される。多数のマイクロブロックが一緒になった統計量を加えて、より大きい体積に対しての統計量を累算することが可能でもある。加えて、各々の統計演算は、マスクを与えることにより、より小さいマイクロブロックの統計量を集計する能力をも提供する。状態質量はマイクロブロックのセルの２×２×２配列の中に特定の状態を占める粒子の数である。以下のように決定される。ここで、Ｎ_i(ｘ)は、ｘで表されるマイクロブロックの位置での８個のボクセルのうちの１つの状態の数である。マイクロブロック全体にわたって状態質量Ｄ_i は状態ごとに累算される。このＭ(ｘ)は、論理和演算によって論理的「１」に設定されたときボクセルｘの累算を可能にするマスクである。Ｍ(ｘ)がボクセルｘに「０」にセットするならば、０はそのボクセルに対して累算される。質量は、マイクロブロック内のすべての粒子を用いるのであって、与えられた状態の粒子を用いないので状態質量とは異なる。以下に述べるように質量は状態ベクトルから計算される。ここで、Ｎ_i(ｘ)は、ｘで表されるマイクロブロックの位置での８個のボクセルのうちの１つの状態の数である。上で見ることができるように、質量Ｎは与えられたマイクロブロックに対する全ての状態での粒子の和である。流体のエネルギーは質量に速度の２乗を掛けた２分の１であるとして決定される。格子気体モデルにおいて、エネルギー準位は離散的な速度状態の組に制限されるために量子化される。各分離化した状態は離散的なエネルギー準位の１つとなる。マイクロブロックのエネルギーは以下のように単純に計算される。ここで、Ｎ_i(ｘ)は、ｘで表されるマイクロブロックの位置での８個のボクセルのうちの１つの状態の総数である。Ｕ_Totは特定のマイクロブロックでの全エネルギーであり、Ｅ_iは状態ｉの離散的なエネルギー準位である。流体の運動量は４成分からなり、各々４次元の１つである。運動量は質量×速度である。この運動量は以下のように計算される。あるいは、この式はベクトル表記で明記することができる。ここで、Ｎ_i(ｘ)は、ｘで表されるマイクロブロックの位置での８個のボクセルのうちの１つの状態の総数である。ｘ次元での状態ｉに対する格子速度はＣ_xi によって示される。状態ｉに対する４次元の速度ベクトルはｃ _iによって示される。他の統計量はこれまでに計算してき統計量から直接得られ、これらは速度、温度及び力を含む。流体速度は４成分からなり、各々４次元の１つである。各成分は速さによって与えられる状態での粒子の数を掛けることで決定される。その次元での状態が示す速さによって与えらる状態の粒子の数を掛け、粒子の全数で和を割って、各成分は決定される。これは運軌量を質量で割っても表すことができる。この４成分は次のように計算される。ここで、ｕ_xはx次元での平均の格子速度である。同様のことはｕ_y、ｕ_z及びｕ_w にも当てはまる。Ｎは上で記したように密度である。速度は次のようにベクトル表示でも表すことができる。温度の計算は計算するのにより複雑な巨視的データの１つである。この値は以下のように質量、エネルギー及び速度から計算される。ここで、Ｕ_Totはエネルギー、Ｎは質量、ｕは流体速度である（全て上で記したもの）。力の計算は流体が物体の表面に及ぼしている(assert)力を決定するのに用いられる。これは表面で粒子の衝突が発生する前の運動量から衝突後の表面境界での運動量を減じることによって決定される。この力も４次元のベクトル量である。この計算をベクトル表示すると以下のようになる。（衝突の）前後の運動量ベクトルＰ ^Before及びＰ ^Afterは、上の運動量の節で述べたように計算される。表面での流体の衝突のために失われる運動量を計算することによって、力は物体の表面に分け与える有効運動量を測定することと同じである。この運動量の計算は全てのマイクロブロック内のボクセルについて、上で記したように計算される。しかし、力は物体の表面に存在するこれらのボクセルに関係するのみであり、全てのマイクロブロック内のボクセルは表面にいる必要はない。しかし、厳密に流体であるこれらのボクセルに対して、衝突の微視的動力学は運動量を保存することを保証する。（衝突の）前後の運動量に流体セルの運動量に含むことは正確な力の計算を可能にすることを無効にする。代わりに、マスクＭ(ｘ)は運動量の計算で累算されることから、いくつかの流体ボクセルを不可能（計算対象外）に設定することができる。図６を再び参照すると、ＡＳＩＣ１２のプログラムモデルは主なものであり、ＡＳＩＣ１２からなる特定目的のコプロセッサの１組にベクトル的な命令を発行するＲＩＳＣプロセッサ２０であることが理解できる。ＲＩＳＣプロセッサ２０は、流体の処理の重要な内部制御ループ（critical inner control loops）を記憶できるＰＣＵ１８内に命令キャッシュ２２を有している。ＲＩＳＣプロセッサ２０はメモリバンク１４に記憶される全てのデータに直接アクセスする。ＲＩＳＣプロセッサ２０は３２ビット整数の標準的な組のアルゴリズム及び論理演算を実行することができる。要するに、ＲＩＳＣプロセッサは標準プロセッサが符号化するのと同様にメモリ内のプログラムを実行するために符号化することができる。メモリマップはプロセッサの４ＧＢアドレス空間内で３つのセグメントを生成するＡＳＩＣ１２に確立されている。 1.)ＭＡＵによって制御される外部メモリバンクに直接アクセスする２ＧＢ(マイナス８ＭＢ)。 2.)機能ユニットに所有されデコードされるレジスタ及び命令をマップされるＡＳＩＣ上のメモリの８ＭＢ。 3.)ＢＩＵを介して接続されるチップ外のホストメモリアクセスの２ＧＢ。ＰＣＵ内の論理（回路）はＲＩＳＣプロセッサによって発せられるアドレスをデコードし、送られるべき場所を決定する。局所的なメモリの要求はＭＡＵに送られる。機能ユニット命令からなるレジスタをマップされたメモリの読み書きは、適当な機能ユニットに送られる。ホストメモリの参照はホストシステムへの要求を送るＢＩＵに送られる。ＭＡＵからメモリ書き換えを要求できる４つの機能ユニットがあり、（それぞれ）プロセッサ制御ユニット、バスインターフェースユニット、マイクロブロックユニット及び移流ユニットである。ＭＡＵはこれらの４つのユニットからの要求間で仲裁し、単一及び多数のワード転送を可能にする。ＭＡＵはバースト転送をサポートするＤＭＡエンジンを保持する。ＭＡＵからの許可を受けて、転送が完了するまで、ユニットはＭＡＵに接続し、データを移動したり受け取ったりする。ＲＩＳＣプロセッサは幾つかの機能ユニット及びメモリアクセスユニット間の大部分のデータパケットの転送を始める。このモデルの機能ユニット若しくはコプロセッサはシステム若しくはＡＳＩＣレベルの処理の理解を必要としない。ユニットはＲＩＳＣエンジンによって処理されるデータ及びタスクが与えられ、そのデータを局所的に処理する。例えば、メモリアクセスユニットによって制御されるメモリからデータの特定のブロックを取り出すために、プロセッサはマイクロブロックユニットに命令しマルチポートＲＡＭに出す（書き込む）。処理が完了した後、ＲＩＳＣエンジンは、ＭＰＲからのデータブロックを取り戻すためにマイクロブロックユニットに命令し、メモリ内の特定の位置に戻す。メモリ要求を開始するＰＣＵに加えて、ＭＡＵメモリ転送の２つのチップ外の元があり、これはＢＩＵを経てホストシステムによって開始されるメモリ要求及びＡＵの通信ポートを経て他のプロセッサから受けられるデータである。ＲＩＳＣエンジン及び局所的に記憶保持されているデータブロックからの単一の命令に基づく数十から数百、数千のサイクルのタスクのルーチンを処理することがコプロセッサは可能である。データ構造を最新のものにしたり次に処理するセグメントを決定するための幾つかの計算をすることに加えて、命令を発し他のユニットの状態をチェックするのに十分自由な時間をＲＩＳＣエンジンは与えられる。書き込みをマップされたメモリを介してＲＩＳＣプロセッサは機能ユニットに命令を発する。機能ユニットに対して命令の演算コードのほかに、書き込みアドレスは機能ユニットを規定する。書き込みによって与えらる３２ビットのデータは演算コードのパラメータである。各演算コードはデータバスをどのように解釈するかを規定する。各機能ユニットは一度に１つの命令を処理できるのみである。１つのユニットは前に発した命令が処理中であるならば、幾つかの新しい命令はＰＣＵで保留される。機械の現在の状態を決定するために各々の機能ユニットの状態とレジスタを読み込むのにもＲＩＳＣプロセッサは自由である。命令をマップされたメモリに発するための必要な複雑な制御ループからＲＩＳＣプロセッサの負担を軽くするために、幾つかの命令をＰＣＵ内で待機するように命令待ち行列は規定されている。ユニットが新しい命令を受けつけるのにできるだけ自由であるように、待ち行列の先頭の命令は適切な機能ユニットで処理される。ユニットが処理中であるならば、待ち行列は命令を発しないでその後の命令は実行されない。この待機のメカニズムは、機能ユニット間の依存性が待機させられている命令を整理して符号化することを可能にする。いつでも、プロセッサは命令待ち行列を回避して自由であり、命令を直接発し、機能ユニットからの情報を読み込む。幾つかの機能ユニットがボクセルデータの処理を利用される方法を示す第１段階の概観をこの節で示す。内部のループは流体マイクロブロック(８個のボクセル)を更新する１つの時間-ステップを実行するように設定する。１つの時間-ステップの更新は衝突、移動及び１つのマイクロブロック内の８個のボクセルに対する統計量の集計に必要なすべての計算からなる。音速以下の乾燥空気のシミュレーションは、２７６個の衝突規則のリストが全てのボクセルに用られ、さらに適切な到達サイトに移動する全ての粒子にも用いられることを要する。マイクロブロックの全５４個の状態が規則記憶(Rules S tore)に与えられている衝突法則で処理される前にがマルチポートＲＡＭにロードされることを衝突処理は必要とする。移流処理は全て状態が６方向(＋ｘ,−ｘ ,＋ｙ,−ｙ,＋ｚ,−ｚ)のそれぞれに移動することを必要とする。しかし、６方向全てにマイクロブロックを移し、次のマイクロブロックを処理することは全く役に立たない。何故なら、データが行き先のマイクロブロックにロードされるまで、一時的にマイクロブロックの外に移動させる粒子を記憶する必要があるからである。移流を処理する効率のより良い方法は、マイクロブロックのグループ上の他方向に粒子を移動させる前に、多数のマイクロブロック上で一度に６方向のうちの１方向に移動する粒子のみを処理することである。他の実施例は次の請求項に含まれる。

───────────────────────────────────────────────────── フロントページの続き (81)指定国ＥＰ(ＡＴ，ＢＥ，ＣＨ，ＤＥ，ＤＫ，ＥＳ，ＦＲ，ＧＢ，ＧＲ，ＩＥ，ＩＴ，ＬＵ，ＭＣ，ＮＬ，ＰＴ，ＳＥ)，ＯＡ(ＢＦ，ＢＪ，ＣＦ，ＣＧ，ＣＩ，ＣＭ，ＧＡ，ＧＮ，ＭＬ，ＭＲ，ＮＥ，ＳＮ，ＴＤ，ＴＧ)，ＡＰ(ＫＥ，ＭＷ，ＳＤ，ＳＺ，ＵＧ)，ＡＭ，ＡＴ，ＡＵ，ＢＢ，ＢＧ，ＢＲ，ＢＹ，ＣＡ，ＣＨ，ＣＮ，ＣＺ，ＤＥ，ＤＫ，ＥＥ，ＥＳ，ＦＩ，ＧＢ，ＧＥ，ＨＵ，ＪＰ，ＫＥ，ＫＧ，ＫＰ，ＫＲ，ＫＺ，ＬＫ，ＬＲ，ＬＴ，ＬＵ，ＬＶ，ＭＤ，ＭＧ，ＭＮ，ＭＷ，ＭＸ，ＮＯ，ＮＺ，ＰＬ，ＰＴ，ＲＯ，ＲＵ，ＳＤ，ＳＥ，ＳＧ，ＳＩ，ＳＫ，ＴＪ，ＴＴ，ＵＡ，ＵＺ，ＶＮ (72)発明者モルビグキムアメリカ合衆国マサチューセッツ州 01742 コンコードモニュメントストリート 1200 (72)発明者テイクセイラクリストファーエム. アメリカ合衆国マサチューセッツ州 02139 ケンブリッジハンコックストリート 118

Claims

【特許請求の範囲】１．物理的過程をシミュレーションする方法を実現するコンピュータであって、（１）複数のボクセルの各々の状態ベクトルを記憶し、各状態ベクトルは複数の整数からなり、ベクトルの各整数はボクセルにおける実現可能な複数の運動量状態のうちの１つの特定の運動量状態に対応し、各整数は特定の運動量状態を有する要素の数を表し、各整数は２つの以上の実現可能な値を有し、（２）異なる運動量状態の要素間の相互作用をモデル化した状態ベクトルの相互作用演算を実行し、（３）シミュレーションされた物理的過程の粘性を変化させる状態ベクトルの粘性修正演算を実行し、（４）新しいボクセルに要素の運動を反映させる状態ベクトルの移動演算を実行するステップからなる。２．物理的過程をシミュレーションするコンピュータシステムであって、（１）複数のボクセルの各々に対する１つの状態ベクトルをメモリに記憶する手段と、複数の整数からなる状態ベクトルの各々と、ボクセルで実現可能な複数の運動量状態のうちの１つの特定の運動量状態に対応するベクトルの整数の各々と、特定の運動量状態を有する要素の数を表す整数の各々と、実現可能な２つ以上の値を有する整数の各々と（２）異なる運動量状態の要素間の相互作用をモデル化する状態ベクトルの相互作用演算を実行する手段と、（３）シミュレーションされた物理的過程の粘性を変化させる状態ベクトルの粘性修正演算を実行する手段と、（４）新しいボクセルの要素に運動を反映させる状態ベクトルの移動演算を実行する手段とからなる。３．請求項１若しくは２に記載されている発明であって、粘性修正演算は緩和手法に従ってシミュレーションされる物理的過程の粘性を変化させる発明。４．請求項３に記載されている発明であって、相互作用演算は粘性修正演算が実行される前にボクセルの状態ベクトルに実行する発明。５．請求項４に記載されている発明であって、相互作用演算及び粘性修正演算の２つは規則の組を状態ベクトルに適用し、各規則は状態ベクトルからの特定の整数の組を修正する発明。６．請求項５に記載されている発明であって、規則は第１の量によって状態ベクトルの整数を修正するための相互作用演算をする間及び第２の量によって状態ベクトルの整数を修正するための相互作用演算をする間に適用され、第１の量及び第２の量は緩和パラメータによって関連づけられる発明。７．請求項６に記載されている発明であって、第２の量は緩和パラメータから得られる数による第１の量を掛けることで決定される発明。８．請求項７に記載されている発明であって、第１の量は相互作用演算をする間に各規則を決定及び記憶され、第２の量は緩和パラメータから得られる数を記憶されている第１の量に乗算することで決定する発明。９．請求項７に記載されている発明であって、第２の量は乗算のステップの結果を切り捨てることによって決定され、そのため第２の量は整数となる発明。 10．請求項９に記載されている発明であって、乗算のステップの結果が切り捨てられる前に、乱数値を乗算のステップの結果に加えられる発明。 11．請求項１０に記載されている発明であって、乱数値は０と１の間の値である発明。 12．請求項６に記載されている発明であって、緩和パラメータは１より大きくかつ２未満である発明。 13．請求項１２に記載されている発明であって、緩和パラメータは１．８以下である発明。 14．請求項３に記載されている発明であって、粘性修正演算はシミュレーションされた物理的過程の粘性を減少させる発明。 15．請求項１若しくは請求項２に記載されている発明であって、要素はボクセル間で運動できる粒子を表し、相互作用演算によってモデル化される相互作用は粒子間若しくは粒子及び境界表面間での衝突であり、移動演算が反射する運動は隣接したボクセル間の粒子の運動である発明。 16．請求項１５に記載されている発明であって、整数値はボクセルでの１つの運動量状態の粒子の数を表す発明。 17．請求項１６に記載されている発明であって、相互作用演算は１つの状態ベクトルの整数の部分集合に作用する衝突規則からなる発明。 18．請求項１７に記載されている発明であって、衝突規則の１つは、状態ベクトルが同じボクセルにおける運動量状態の第１の対を表す整数を有するとき、運動量状態の第２の対に複数の粒子をある条件のもとで、質量、運動量若しくはエネルギーのような１つ以上の物理的不変量の結合した値が第２の対と同じである値を第１の対が有するように移動させ、散乱規則が運動量状態の第２の対に粒子を移動させるとき、質量、運動量若しくはエネルギーのような１つ以上の物理的不変量の結合した値が第２の対と同じである値を保持しつづけるように、粘性修正演算は運動量状態の第１の対から運動量状態の第２の対に粒子の追加の量を移動させる散乱規則である発明。 19．請求項１８に記載されている発明であって、粘性修正演算によって移動される粒子の追加の量は、緩和パラメータから１を減算し、散乱規則によって移動される粒子の数を減算のステップの結果に乗算することによって決定される発明。 20．請求項１９に記載されている発明であって、散乱規則によって移動される粒子の数は相互作用演算の間に確定され記憶され、追加の量は記憶された粒子の数を減算ステップの結果に掛けることで決定される発明。 21．請求項１９に記載されている発明であって、追加の量は乗算のステップの結果を切り捨てることによって決定され、そのため追加の量は１つの整数となる発明。 22．請求項２１に記載されている発明であって、乗算のステップの結果が切り捨てられる前に、乱数値は乗算ステップの結果に加えられる発明。 23．請求項２２に記載されている発明であって、乱数値は０と１の間の値である発明。 24．請求項３に記載されている発明であって、規則の組は第１の量によって状態ベクトルの整数を修正する相互作用演算をする間及び第２の量によって状態ベクトルの整数を修正する相互作用演算をする間に適用され、第１及び第２の量は緩和パラメータによって関係づけられている発明。 25．請求項２４に記載されている発明であって、第２の量は緩和パラメータから得られる数に第１の量を乗算され、乱数値を乗算の結果に加算され、さらに加算の結果を切り捨てられることで決定され、そのため第２の量は１つの整数である発明。