JP2016534436A - 等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子 - Google Patents

Abstract

格子速度セットにおいて流体容積中の粒子の移動をシミュレートするステップであって、移動が粒子間の衝突を引き起こし、シミュレートされた移動に基づいて、容積内の特定位置における粒子の粒子相対速度を決定するステップであって、粒子相対速度が（ｉ）容積内の特定位置において容積の流れがゼロの下で測定される粒子の絶対速度と（ｉｉ）容積内の特定位置における１又は２以上の粒子の平均速度との間の差であるステップと、粒子相対速度に基づいて、衝突を表す特定の次数の非平衡衝突後分布関数を決定するステップと、を含む方法。【選択図】図１３

Description

（優先権の主張）
本願は、２０１３年７月２４日出願の米国特許仮出願第６１／８５８，０５１号に対して米国特許法第１１９条（ｅ）の下で優先権を主張するものであり、その内容全体は、これにより引用によって組み込まれる。

衝突過程は多粒子系において２つの基本的な動力学的過程の内の１つであり、もう１つは移流過程である。衝突過程は、個々の粒子が相互作用して集団的な挙動を成すために本質的なものである。衝突過程中、質量、運動量及びエネルギは保存則に従って粒子間で交換される。これらの保存則により、関係する粒子間の全体的な質量及び運動量（及び時にはエネルギー）は衝突の前後で不変であることが保証される。

米国特許第５５９４６７１号明細書

一般に、本明細書は、格子速度セット(lattice velocity set)において流体容積中の粒子の移動をシミュレートする技術を説明しており、移動が粒子間の衝突を引き起こし、シミュレートされた移動に基づいて容積内の特定位置における粒子の粒子相対速度(relative particle velocity)を決定するが、その粒子相対速度は（ｉ）容積内の特定位置において容積の流れがゼロの下で測定される粒子の絶対速度と（ｉｉ）容積内の特定位置における１又は２以上（複数）の粒子の平均速度と、の間の差であり、粒子相対速度に基づいて衝突を表す特定の次数(specified order)の非平衡衝突後分布関数(non-equilibrium post-collide distribution function)を決定する。本態様の他の実施形態は、対応するコンピュータシステム、装置、機械可読ハードウェア記憶装置、及び１又は２以上のコンピュータ記憶装置に記録されたコンピュータプログラムが含まれ、各々が本方法の動作及び特徴を実行するように構成されている。１又は２以上のコンピュータのシステムは、作動時、システムにその動作を実行させるソフトウェア、ファームウェア、ハードウェア、又はそれらの組み合わせをシステムにインストールすることによって、特定の処理又は動作を実行するように構成することができる。１又は２以上のコンピュータプログラムは、データ処理装置による実行時にデータ処理装置に特定の処理又は動作を実行させる命令を備えることによってその動作を実行するように構成することができる。

前述の及び別の実施形態は、各々選択的に以下の特徴のうちの１又は２以上を単独で又は組み合わせて含むことができる。特に、１つの実施形態は、以下の特徴の全てを組み合わせて含むことができる。特徴は、１又は２以上のコンピュータシステムによって流体力学運動を粒子速度の或る次数までサポートすることができる格子速度セットを提供するステップを含み、シミュレートすることは、１又は２以上のコンピュータシステムでシミュレートすることを含む。また、特徴は、格子速度セットのサポートされる次数(supported order)が、非平衡衝突後分布関数の特定の次数より小さくかつこれとは異なること、非平衡衝突後分布関数の特定の次数が、粒子速度の次数によって決定されることを含む。

また、特徴は、容積内の特定位置における１又は２以上の粒子の平均速度がその特定位置における特定種類の粒子の平均速度を有することを含む。また、特徴は、格子速度セットが格子ボルツマン法と関連する状態ベクトルのセットであることを含む。また、特徴は、非平衡衝突後分布関数が、（ｉ）予め定義された物理量に対する非平衡モーメントを維持し、（ｉｉ）未定義の物理量に対する非平衡モーメントを特定の次数まで除去することを含む。また、特徴は、特定の次数が流体速度の格子音速に対する比率に関連した指数値(exponential value)であり、格子速度セットが指数値をサポートすることを含む。また、特徴は、格子速度セットが、格子に限定される空間で運動量状態のセットを有することを含む。また、特徴は、粒子相対速度が、容積内の特定位置において容積の流れがゼロの下で測定される粒子の絶対速度から容積内の特定位置における１又は２以上の粒子の平均速度を差し引いたものであることを含む。また、特徴は、非平衡衝突後分布関数がガリレイ不変フィルタ処理演算子(Galilean invariant filtered operator)であることを含む。また、特徴は、非平衡衝突後分布関数に基づいて流体容積中での粒子の衝突過程をモデル化するステップを含む。また。特徴は、非平衡衝突後分布関数が格子速度セットに対するマッハ数に関して１次のガリレイ不変の衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
ｘは容積内の特定位置であり、ｔは時間における特定時点であり、ｉは格子速度セットでの格子速度の指数であり、Ｔ₀は一定格子温度であり、ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける粒子間の平均速度であり、Ｉは２階の単位テンソルであり、τは衝突緩和時間であり、ｗ_iは一定重み係数であり、Π^neqは非平衡運動量流束である。

また、特徴は、非平衡衝突後分布関数が、流体力学モーメントに関して無限次をサポートする、格子速度セットに対する衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
ｘは容積内の特定位置であり、ｔは時間における特定時点であり、ｉは格子速度セットでの格子速度の指数であり、Ｔ₀は一定格子温度であり、Ｉは２階の単位テンソルであり、τは衝突緩和時間であり、
は相対的な粒子速度であり、ρは流体密度であり、
は平衡分布関数であり、
は非平衡運動量流束である。また、特徴は、非平衡衝突後分布関数が、流体力学モーメントに関して２次のサポートを提供する、格子速度セットに対するマッハ数に関して２次のガリレイ不変の衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
ｘは容積内の特定位置であり、ｔは時間における特定時点であり、ｉは格子速度セットでの格子速度の指数であり、Ｔ₀は一定格子温度であり、ｃ_iは衝突前の粒子の速度ベクトルであり、ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける粒子間の平均速度であり、Ｉは２階の単位テンソル(second rank unity tensor)であり、τは衝突緩和時間であり、ｗ_iは一定重み係数であり、
は非平衡運動量流束である。また、特徴は、予め定義された物理量が特定の容積中の流体質量、特定の容積中の流体運動量、及び特定の容積中の流体エネルギを有することを含む。

また、特徴は、非平衡衝突後分布関数が、エネルギ流束に関連する衝突演算子
であることを含み、衝突演算子は次式に従って定義され、
ｘは容積内の特定位置であり、ｔは時間における特定時点であり、ｉは格子速度セットでの格子速度の指数であり、Ｔ₀は一定格子温度であり、Ｉは２階の単位テンソルであり、τは衝突緩和時間であり、
は相対的な粒子速度であり、
は平衡分布関数であり、
は非平衡エネルギ流束である。

上述の技術の実装は、方法又は過程、システム又は装置、又はコンピュータアクセス可能媒体上のコンピュータソフトウエアを含むことができる。

システム及び方法及び技術は、多相流れに対するＳｈａｎ−Ｃｈｅｎ方法と格子ボルツマン方式とのような様々なタイプの数値シミュレーション手法を使用して実施することができる。格子ボルツマン方式に関する追加の情報を本明細書で以下に説明する。しかし、本明細書に説明するシステム及び技術は、格子ボルツマン方式を使用するシミュレーションに限定されず、他の数値シミュレーション手法に適用することができる。

システム及び技術は、格子ボルツマン方式を使用する格子ガスシミュレーションを使用して実施することができる。従来の格子ガスシミュレーションは、各格子部位での限られた数の粒子を仮定し、粒子は、ビットのショートベクトルによって表される。各ビットは、特定の方向に移動する粒子を表している。例えば、ベクトル中の１ビットは、特定の方向に沿って移動する粒子の存在（１に設定時）又は不在（０に設定時）を表すことができる。そのようなベクトルは、６ビットを有することができ、例えば、値１１００００は、２つの粒子がＸ軸に沿って反対方向に移動し、かつＹ軸及びＺ軸に沿って移動する粒子がないことを示している。１組の衝突規則が、各部位での粒子間の衝突の挙動を支配する（例えば、１１００００ベクトルは、００１１００ベクトルになることができ、Ｘ軸に沿って移動する２つの粒子間の衝突は、Ｙ軸に沿って離れる２つの粒子を生成したことを示す）。規則は、ビットに対して置換を実行する（例えば、１１００００を００１１００に変形する）ルックアップテーブルに状態ベクトルを供給することによって実施される。粒子は、その後に、隣接部位に移動される（例えば、Ｙ軸に沿って移動する２つの粒子は、Ｙ軸に沿って左右に隣接部位に移動されると考えられる）。

強化されたシステムにおいて、各格子部位での状態ベクトルは、粒子エネルギ及び移動方向の変動を与えるためにより多くのビットを含み（例えば、亜音速流れに対して５４ビット）、完全な状態ベクトルの部分集合を伴う衝突規則が使用される。更に強化されたシステムにおいて、１つよりも多い粒子が、各格子部位又はボクセル（これらの２つの用語は、本明細書を通して交換可能に使用される）で各運動量状態で存在することが許される。例えば、８ビット実装において、０〜２５５個の粒子は、特定のボクセルで特定の方向に移動している可能性がある。状態ベクトルは、１組のビットである代わりに、１組の整数であり（例えば、１組の８ビットバイトが、０〜２５５の範囲の整数を与える）、その各々は、与えられた状態にある粒子の数を表する。

更なる強化において、格子ボルツマン方法（ＬＢＭ）は、従来の計算流体力学（「ＣＦＤ」）手法で可能であるよりも深いレベルで複合的な幾何学形状における３Ｄ不安定圧縮性乱流過程をシミュレートするために流体のメゾスコピック表現を使用する。ＬＢＭ方法の概要を以下に与える。

ボルツマンレベルのメゾスコピック表現(Boltzmann-Level Mesoscopic Representation)
流体システムは、いわゆる「メゾスコピック」レベルに対する運動方程式によって表すことができることが統計物理学で公知である。このレベルに対しては、個々の粒子の詳細な運動は決定する必要はない。これに代えて、流体の特性は、単一粒子位相空間、
を使用して定義される粒子分布関数によって表され、ｘは、空間座標であり、νは、粒子速度座標である。質量、密度、流量速度、及び温度のような典型的な流体力学的な量は、粒子分布関数の単純なモーメントである。粒子分布関数の力学は、ボルツマン方程式に従う。
方程式（１）
ここで、Ｆ（ｘ，ｔ）は、（ｘ，ｔ）で外部又は自己矛盾なく発生された物体力を表している。衝突項Ｃは、様々な速度及び位置の粒子の相互作用を表している。衝突項Ｃに対して特定の形態を指定することなく、上述のボルツマン方程式は、希薄ガスの公知の状況（ボルツマンによって本来構成されたような）だけでなく全ての流体システムに適用可能であることを強調することが重要である。

一般的に、Ｃは、二点相関関数の複雑な多次元積分を含む。分布関数ｆのみで閉じたシステムを形成するために、並びに効率的な計算のために、最も便利かつ物理的に一貫した形態の１つは、公知のＢＧＫ演算子である。ＢＧＫ演算子は、衝突の詳細がたとえ何であろうとも、分布関数は、衝突：
方程式（２）
を通じて
によって与えられる明確な局所平衡に近づくという物理的議論に従って構成され、ここで、パラメータτは、衝突を通じた平衡までの固有緩和時間を表している。粒子（例えば、原子又は分子）を扱うと、緩和時間は、典型的には定数として取られる。「混成」（水力−動力学）表現において、この緩和時間は、歪み率及び乱流運動エネルギなどのような流体力学的変数の関数である。すなわち、乱流は、局所的に決定された固有特性を有する乱流粒子（「渦」）のガスとして表すことができる。

ボルツマン−ＢＧＫ方程式の数値解は、ナビエ−ストークス方程式の解に優るいくつかの計算上の利点を有する。第１に、複雑な非線形項又は高次空間導関数が方程式にないことを直ちに認識することができ、すなわち、移流不安定性に関する問題がほとんどない。このレベルの記述では、方程式は、圧力を扱う必要がないので局所的であり、これは、アルゴリズム並列化に対してかなりの利点を提供する。２次空間導関数を有する拡散演算子がないという事実と共に線形移流演算子の別の望ましい特徴は、流体偏微分方程式（「ＰＤＥ」）の数学的な条件ではなく、実際に粒子が固体面と真に相互作用する方法を模倣するように非スリップ面又はスリップ面のような物理的境界条件を実現する際の容易さである。直接の利益の１つは、固体面上でインタフェースの移動を処理する問題がないということであり、これは、格子ボルツマンベースのシミュレーションソフトウエアが失敗なく複合乱流空気力学をシミュレートすることを可能にするのを補助する。これに加えて、有限粗度面のような境界からのある一定の物理特性も、力に組み込むことができる。更に、ＢＧＫ衝突演算子は、純粋に局所的であり、一方、自己矛盾のない物体力の計算は、近−隣接情報のみを通して達成することができる。その結果、ボルツマン−ＢＧＫ方程式の計算は、並列処理に実質的に適応させることができる。

格子ボルツマン方式
連続体ボルツマン方程式を解くことは、それが位置及び速度位相空間における積分微分方程式の数値評価を伴うということで有意な課題を表している。位置だけでなく速度位相空間も離散化することができることが観測された時に大規模な簡素化が行われ、それは、ボルツマン方程式の解のための効率的な数値アルゴリズムをもたらした。流体力学的な量は、高々最近隣接情報に依存する単純和の項で書くことができる。従来的には、格子ボルツマン方程式の方式は、速度の離散集合：
に対する粒子の推移を規定する格子ガスモデルに基づいていたが、この方程式は、連続体ボルツマン方程式の離散化としての第１の原理から系統的に導出することができる。その結果、ＬＢＥは、格子ガス手法に関連付けられた公知の問題を被らない。従って、位相空間における連続体分布関数ｆ（ｘ，ν，ｔ）を取り扱う代わりに、離散速度指数をラベル付けする下付き文字を有する離散分布の有限集合
を追跡することが必要なだけである。巨視的記述の代わりにこの動態方程式を取り扱う重要な利点は、システムの位相空間の増大が、問題の局所性によってオフセットされるということである。

対称性の考慮に起因して、速度値の集合は、それらが構成空間内で張られた時にある一定の格子構造を形成するように選択される。そのような離散システムの力学は、形態：
を有するＬＢＥに従い、ここで、衝突演算子は、上述のように通常はＢＧＫ形態を取る。平衡分布形態の適正な選択により、格子ボルツマン方程式は正しい流体力学及び熱−流体力学を生じさせることを理論的に示すことができる。すなわち、
から導出された流体力学モーメントは、巨視的限界においてナビエ−ストークス方程式に従う。これらのモーメントは、
方程式（３）
して定義され、ここで、ρ、ｕ、及びＴは、それぞれ、流体密度、速度、及び温度であり、Ｄは、離散化された速度空間の次元（物理空間次元に全く等しくない）である。

他の特徴及び利点は、図面及び特許請求の範囲を含めて、以下の説明から明らかとなる。

ＬＢＭモデルの速度成分を例示する。 ＬＢＭモデルの速度成分を例示する。 物理過程シミュレーションシステムが従う手順の流れ図である。 マイクロブロックの斜視図である 図３のシステムで使用する格子構造を例示する。 図３のシステムで使用する格子構造を例示する。 可変分解能技術を例示する。 可変分解能技術を例示する。 表面のファセットにより影響を受ける領域を例示する。 ボクセルから表面への粒子の移動を例示する。 表面から表面への粒子の移動を例示する。 表面動力学を実行するための手順の流れ図である。 特定次数の非平衡衝突後分布関数を決定するするプロセスの流れ図である。 特定次数の非平衡衝突後分布関数を決定するするシステムの構成要素のブロック図である。

Ａ．特定の非平衡移動を保持し他の非平衡移動を除去する衝突演算子
格子ボルツマン・シミュレーションなどのシミュレーションシステムでは、シミュレーション空間は多数の離散点に分割され、それらは直線で連結されて離散的な数の点及び方向を与える。シミュレーションはまた、時間ステップの離散集合に制約される。このようなシステムでは、シミュレーションが現実世界の流れを近似するために、多数の異なる量が保存されなければならない。例えば、システムは質量、運動量及びエネルギを保存する。従って、シミュレーションは、適切な質量流束、運動量流束、及びエネルギ流束を有するように構成する必要がある。これらの保存量はその流束と共に、真実の物理世界と関係するシミュレーションシステムにおいて本質的なモーメントである。しかしながら、これらの量を保存する場合に、シミュレーションは離散速度空間（例えば、粒子が所定の時間ステップで進むことのできる方向及び距離の離散集合）に起因して非意図的に付加的なモーメント量を作り出すことがある。非意図的に生成された量（本明細書では、意図しない又は望ましくない不変量、保存される又は非平衡モーメントとして言及する）は、シミュレーション結果に悪影響を与えることがある。例えば、このような望ましくない量は、誤った流体の動的挙動と計算結果の数値的不安定性とをもたらす場合がある。

非意図的に生成される不変量の影響を低減するために、保存される物理量に関してのみ非平衡モーメントを保持するが、残りの非平衡モーメントを所望の次数まで除去する、衝突演算子を本明細書に説明する。本明細書では、衝突演算子の次数は格子マッハ数（流体速度と格子音速の比）の指数を用いて定義される。望ましくない非平衡モーメントをフィルタ処理することを確実にするためのトレードオフは、計算時間と処理手続きの増大である。本明細書に説明する衝突演算子では、理論的な形式は運動量及びエネルギの非平衡流束の両方に対する任意の次数まで体系的に構築される。従って、これらの衝突演算子は、質量、運動量及びエネルギの保存を満たし、正確な質量流束、運動量流束及びエネルギ流束を選択した次数まで保証するが、同時に、望ましくない非平衡モーメントのすべてを選択した次数まで除去する。本シミュレーションシステムでは、高速流又は低粘性を取り扱う際に、フィルタ処理のスキームに関してより高次オーダーを選択することが、より低次オーダーの場合に比べて非物理的効果又は影響をより少なくする。

Ｂ．シミュレーション空間のモデル化
ＬＢＭベースの物理過程シミュレーションシステムでは、離散速度ｃ_iの集合で評価される分布関数値ｆ_iによって、流体の流れを表現することができる。分布関数の力学は、方程式４により支配され、ここにｆ_i（０）は平衡分布関数として知られ、
方程式（４）
として定義され、ここで、
である。
方程式（５）
この方程式は、分布関数ｆ_iの時間発展を説明する、よく知られた格子ボルツマン方程式である。左辺はいわゆる「流動過程」による分布の変化を表す。流動過程は、流体ポケットが或るグリッド位置から出発して、それから速度ベクトルの１つに沿って次のグリッド位置へ移動する過程である。その地点で、「衝突演算子」つまり近接した流体ポケットの開始流体ポケットへの影響を計算する。流体は別のグリッド位置に移動することだけが可能なので、全速度の成分全体は共通した速度の倍数であるように速度ベクトルの適切な選択が必要である。

第１の方程式の右辺は、流体ポケット間の衝突による分布関数の変化を表す前述の「衝突演算子」である。本明細書で使用する衝突演算子の特別な形式は、Ｂｈａｔｎａｇａｒ、Ｇｒｏｓｓ及びＫｒｏｏｋ（ＢＧＫ）によるものである。衝突演算子は、分布関数を、「平衡」形である第２の式により与えられる規定値へ向かうように強いる。

このシミュレーションから、密度ρ及び流速ｕなどの従来の流体変数は、式（３）の単純和として得ることができる。ここで、ｃ_i及びｗ_iの集団的な値はＬＢＭモデルを規定する。ＬＢＭモデルは、大規模に実現可能なコンピュータ・プラットフォームで効率的に実行することができ、時間的に不安定な流れ及び複雑な境界条件に対して高い堅牢性で走らせることができる。

ボルツマン方程式から流体システムに対する巨視的な運動方程式を得るための標準的な技術は、Ｃｈａｐｍａｎ−Ｅｎｓｋｏｇ法であり、そこでは完全なボルツマン方程式の逐次近似がなされる。

流体システムでは、密度の小さな擾乱が音速で伝播する。気体系では、音速は一般に温度で決まる。流れにおける圧縮性の影響の重要性は、特徴的な速度と音速との比によって評価され、それはマッハ数として知られる。

図１を参照すると、第１モデル（２Ｄ−１）１００は２１個の速度を含む２次元モデルである。これら２１個の速度の内、１個（１０５）は動いていない粒子を表し、３組の４速度は、正規化速度（ｒ）（１１０−１１３）、正規化速度の２倍（２ｒ）（１２０−１２３）又は正規化速度の３倍（３ｒ）（１３０−１３３）で、格子のｘ軸又はｙ軸に沿って正方向又は負方向に動いている粒子を表し、２組の４速度は、正規化速度（ｒ）（１４０−１４３）又は正規化速度の２倍（２ｒ）（１５０−１５３）で、２組の４速度は、格子のｘ軸とｙ軸の両方に対して動いている粒子を表す。

さらに図２に図示するように、第２モデル（３Ｄ−１）２００は３９個の速度を含む３次元モデルであり、各速度は図２の矢印の内の１つにより表されている。これら３９個の速度の内、１個は動いていない粒子を表し、３組の６速度は、正規化速度（ｒ）、正規化速度の２倍（２ｒ）又は正規化速度の３倍（３ｒ）で、格子のｘ軸、ｙ軸又はｚ軸に沿って正方向又は負方向に動いている粒子を表し、８個は、正規化速度（ｒ）でｘ、ｙ、ｚの全３格子軸に対して動いている粒子を表し、１２個は、正規化速度の２倍（２ｒ）でｘ、ｙ、ｚの格子軸の内の２つに対して動いている粒子を表す。

例えば、３Ｄ−２は１０１個の速度を含み、２Ｄ−２モデルは３７個の速度を含むなど、さらに複雑なモデルも使用することができる。

３次元モデル３Ｄ−２に関しては、１０１個の速度の内、１個は動いていない粒子を表し（グループ１）、３組の６速度は、正規化速度（ｒ）、正規化速度の２倍（２ｒ）又は正規化速度の３倍（３ｒ）で、格子のｘ軸、ｙ軸又はｚ軸に沿って正方向又は負方向に動いている粒子を表し（グループ２、４及び７）、３組の８速度は、正規化速度（ｒ）、正規化速度の２倍（２ｒ）又は正規化速度の３倍（３ｒ）で、ｘ、ｙ、ｚの全３格子軸に対して動いている粒子を表し（グループ３、８及び１０）、１２個は、正規化速度の２倍（２ｒ）でｘ、ｙ、ｚの格子軸の内の２つに対して動いている粒子を表し（グループ６）、２４個は、正規化速度（ｒ）及び正規化速度の２倍（２ｒ）でｘ、ｙ、ｚの格子軸の内の２つに対して動いており、残る１軸に対しては動いていない粒子を表し（グループ５）、２４個は、正規化速度（ｒ）でｘ、ｙ、ｚの格子軸の内の２つに対して動いており、残る１軸に対しては正規化速度の３倍（３ｒ）で動いている粒子を表す（グループ９）。

２次元モデル２Ｄ−２に関しては、３７個の速度の内、１個は動いていない粒子を表し（グループ１）、３組の４速度は、正規化速度（ｒ）、正規化速度の２倍（２ｒ）又は正規化速度の３倍（３ｒ）で、格子のｘ軸又はｚ軸に沿って正方向又は負方向に動いている粒子を表し（グループ２、４及び７）、２組の４速度は、正規化速度（ｒ）又は正規化速度の２倍（２ｒ）で、格子のｘ軸とｙ軸の両方に対して動いている粒子を表し、８速度は、正規化速度（ｒ）で格子のｘ軸とｙ軸の内の１つに対して、他方の軸に対しては正規化速度の２倍（２ｒ）で動いている粒子を表し、８速度は、正規化速度（ｒ）で格子のｘ軸とｙ軸の内の１つに対して、他方の軸に対しては正規化速度の３倍（３ｒ）で動いている粒子を表す。

上記のＬＢＭモデルは、２次元及び３次元の流れの数値シミュレーション用の効率的で堅牢な、離散速度の動力学的モデルの特定のクラスを提供する。この種のモデルには、離散速度とこれらの速度に関係する重みとの特定の組が含まれる。速度は速度空間においてデカルト座標のグリッド点と一致し、これにより、離散速度モデル、特に格子ボルツマンモデルとして知られる種類の正確で効率的な実行が容易になる。このようなモデルを用いて、流れを高い忠実度でシミュレートすることができる。

図３を参照すると、物理過程シミュレーションシステムは、流体の流れなどの物理過程をシミュレートするために手順３００に従って作動する。シミュレーションに先立って、シミュレーション空間をボクセルの集合としてモデル化する（ステップ３０２）。通常、シミュレーション空間はコンピュータ支援設計（ＣＡＤ）を用いて生成される。例えば、ＣＡＤプログラムは、風洞に置かれた微小デバイスを描くために使用することができる。その後に、適切な分解能を備える格子構造を加えるために、及びシミュレーション空間内の物体及び表面に対処するために、ＣＡＤプログラムにより作成されたデータを処理する。

格子の分解能は、シミュレートするシステムのレイノルズ数に基づいて選択することができる。レイノルズ数は、流れの粘度（ν）、流れ中の物体の特徴的な長さ（Ｌ）、及び流れの特徴的な速度に関係し、
方程式（６）
である。

物体の特徴的な長さは、物体の大規模な特徴を表す。例えば、微小デバイスの周りの流れがシミュレートされる場合には、微小デバイスの高さは、特徴的な長さであると見なすことができる。物体の小領域（例えば、自動車のサイドミラー）周りの流れが関心領域である場合には、シミュレーションの分解能を増大させることができ、又は関心領域の周りには分解能が増大した領域を使用することができる。ボクセルの寸法は、格子の分解能が増加するにつれて減少する。

状態空間は、ｆ_i（ｘ，ｔ）として表され、ｆ_iは、時間ｔにおいて３次元ベクトルｘにより示される格子部位で状態ｉにある単位容積当たりの要素又は粒子の数（つまり、状態ｉの粒子密度）を表す。既知の時間増分に対しては、粒子数を単にｆ_i（ｘ）と呼ぶ。１つの格子部位の全ての状態の組み合わせをｆ（ｘ）として示す。

状態の数は、各エネルギ準位内の可能な速度ベクトルの数により決定される。速度ベクトルは、３次元ｘ、ｙ、及びｚを備える空間において整数値の線速度から構成される。状態の数は、複数の種シミュレーションに対して増大される。

各状態ｉは、特定エネルギ準位での異なる速度ベクトルを表す（つまり、エネルギ準位０、１又は２）。各状態の速度ｃ_iは、以下のように３つの次元の各々におけるその「速度」で示す。
方程式（７）

エネルギ準位０状態は、いずれの次元でも動いていない、つまりｃ_stopped＝（０、０、０）である静止粒子を表す。エネルギ準位１状態は、３つの次元の内の１つで±１の速度、他の２つの次元では速度０を有する粒子を表す。エネルギ準位２状態は、全ての３つの次元で±１の速度、又は３つの次元のうちの１つで±２の速度及び他の２つの次元で速度０を有する粒子を表す。

３つのエネルギ準位の可能な順列の全てを生成することにより、合計３９個の可能な状態（１個のエネルギ０状態、６個のエネルギ１状態、８個のエネルギ３状態、６個のエネルギ４状態、１２個のエネルギ８状態、及び６個のエネルギ９状態）が得られる。

各ボクセル（つまり、各格子部位）は、状態ベクトルｆ（ｘ）により表される。状態ベクトルはボクセルのステータスを完全に規定し、３９個のエントリを含む。３９個のエントリは、１個のエネルギ０状態、６個のエネルギ１状態、８個のエネルギ３状態、６個のエネルギ４状態、１２個のエネルギ８状態、及び６個のエネルギ９状態に対応する。この速度セットを使用することによって、システムは、実現された平衡状態ベクトルに関してマックスウェル−ボルツマン統計を生成することができる。

処理効率のために、ボクセルをマイクロブロックと呼ばれる２ｘ２ｘ２ボリュームにグループ化する。マイクロブロックは、ボクセルの並列処理を可能にし、かつデータ構造に関連づけされたオーバーヘッドを最小化するように構成される。マイクロブロック内のボクセルの略記法は、Ｎ_i（ｎ）として規定され、ｎは、マイクロブロック内の格子部位の相対位置を表し、ｎ∈｛０，１，２，．．．，７｝である。マイクロブロックを図４に示す。

図５Ａ及び５Ｂを参照すると、表面Ｓはシミュレーション空間においてファセットＦ_αの集合として表され、
方程式（８）
ここで、αは、特定のファセットを列挙する指標である。ファセットは、ボクセルの境界に限定されず、通常は、ファセットが比較的小さな数のボクセルに影響を与えるように、ファセットを隣接するボクセルの大きさの程度に又はそれより僅かに小さくサイズ設定される。表面動力学を実行するために、ファセットに特性を割り当てる。特に、各ファセットＦ_αは、単位法線（ｎ_α）、表面積（Ａ_α）、中心位置（ｘ_α）、及びファセットの表面動力学特性を説明するファセット分布関数（ｆ_i（α））を有する。

図６を参照すると、処理効率を向上させるために、異なるレベルの分解能をシミュレーション空間の異なる領域で使用することができる。通常は、物体６５５の周りの領域６５０が最大関心領域なので、最高分解能でシミュレートされる。粘性の影響は物体からの距離と共に減少するので、低レベルの分解能（つまり、拡張されたボクセルボリューム）を使用して、物体６５５から増大する距離で離間する領域６６０、６６５をシミュレートする。同様に、図７に図示するように、より低レベルの分解能を使用して、物体７７５のより有意でない特徴部の周りの領域７７０をシミュレートすることができるが、最高レベルの分解能を使用して、最も有意な特徴部（例えば、先端面及び後端面）の周りの領域７８０をシミュレートする。周辺領域７８５は、最低レベル分解能と最大のボクセルとを用いてシミュレートする。

Ｃ．ファセットに影響を受けるボクセルの特定
図３を再び参照すると、シミュレーション空間をモデル化した状態で（ステップ３０２）、１又は２以上のファセットにより影響を受けるボクセルを特定する（ステップ３０４）。ボクセルは、多くの点でファセットの影響を受ける場合がある。最初に、１又は２以上のファセットと交差するボクセルは、交差しないボクセルと比べて低減された容積を有するという点で影響を受ける。これは、ファセット及びファセットにより表される表面下にある物質がボクセルの一部を占めるので発生する。分画因子Ｐ_f（ｘ）は、ファセットの影響を受けないボクセルの部分（つまり、流れがシミュレートされる流体又は他の物質が占める可能性がある部分）を示す。交差しないボクセルに対しては、Ｐ_f（ｘ）は１に等しい。

粒子をファセットに移動するか又は粒子をファセットから受け取ることによって１又は２以上のファセットと相互作用するボクセルも、ファセットにより影響を受けるボクセルとして特定される。ファセットと交差する全てのボクセルは、ファセットから粒子を受け取る少なくとも１つの状態と、粒子をファセットへ移動する少なくとも１つの状態とを含むことになる。ほとんどの場合に、付加的なボクセルは、このような状態を含むことになる。

図８を参照すると、０以外の速度ベクトルｃ_iを有する各状態ｉに対して、ファセットＦ_αは、平行六面体Ｇ_iαにより画定される領域から粒子を受け取るか又はこの領域へ粒子を移動するが、この平行六面体は速度ベクトルｃ_iとファセットの単位法線とのベクトルドット積（｜ｃ_iｎ_i｜）の大きさにより規定される高さとファセットの表面積Ａ_αにより規定される底面とを有するので、平行六面体Ｇ_iαの容積Ｖ_iαは次のようになる。
方程式（９）

ファセットＦ_αは、状態の速度ベクトルがファセットに向けられた時に（｜ｃ_iｎ_i｜＜０）、容積Ｖ_iαから粒子を受け取り、状態の速度ベクトルがファセットから離れるように向けられた時に（｜ｃ_iｎ_i｜＞０））、その領域へ粒子を移動する。以下に説明するように、この式は、別のファセットが平行六面体Ｇｉαの一部、すなわち、内部のコーナのような非凸面特徴部の近くに発生する可能性がある状態を占有した時には修正する必要がある。

ファセットＦ_αの平行六面体Ｇ_iαは、複数のボクセルの内の一部又は全てと重なる場合がある。ボクセル又はその部分の数は、ボクセルのサイズに対するファセットのサイズ、状態のエネルギ、及び格子構造に対するファセットの向きに依存する。影響を受けるボクセルの数は、ファセットのサイズと共に増加する。従って、ファセットのサイズは、上述のように、通常はファセット近くに位置するボクセルのサイズと同程度に又はそれより小さいように選択される。

平行六面体Ｇ_iαと重なるボクセルＮ（ｘ）の一部分は、Ｖ_iα（ｘ）として定義する。この項を使用して、ボクセルＮ（ｘ）とファセットＦ_αとの間を移動する状態ｉ粒子の流束Γ_iα（ｘ）は、ボクセル（Ｎ_i（ｘ））における状態ｉ粒子の密度にボクセルと重なる領域の容積（Ｖ_iα（ｘ））を乗じたものに等しい。
方程式（１０）

平行六面体Ｇ_iαが１又は２以上のファセットと交差する場合、以下の条件が正しく、
方程式式（１１）
ここで、第１の和は、Ｇ_iαと重なる全ボクセルに相当し、第２の項は、Ｇ_iαと交差する全てのファセットに相当する。平行六面体Ｇ_iαが別のファセットと交差しない場合、この式は以下に帰する。
方程式（１２）

Ｄ．シミュレーションの実行
１又は２以上のファセットにより影響を受けるボクセルが特定された状態で（ステップ３０４）、シミュレーションを開始するためにタイマを初期化する（ステップ３０６）。シミュレーションの各時間増分中に、粒子のボクセルからボクセルへの移動は、粒子の表面ファセットとの相互作用に対応する移流段階（ステップ３０８−３１６）によりシミュレートされる。次いで、衝突段階（ステップ３１８）は、各ボクセル内の粒子の相互作用をシミュレートする。次に、タイマを増分する。増分されたタイマがシミュレーションの完了（ステップ３２２）を示さない場合、移流段階及び衝突段階（ステップ３０８−３２０）が繰り返される。増分されたタイマがシミュレーションの完了（ステップ３２２）を示す場合、シミュレーションの結果を記憶及び／又は表示する（ステップ３２４）。

１．表面に対する境界条件
表面との相互作用を正確にシミュレートするために、各ファセットは、４つの境界条件を満たす必要がある。第１に、ファセットが受け取る粒子の合計質量は、ファセットが移送する粒子の合計質量に等しくなければならない（つまり、ファセットに対する正味の質量流束は０でなければならない）。第２に、ファセットが受け取る粒子の合計エネルギは、ファセットが移送する粒子の合計エネルギに等しくなければならない（つまり、ファセットに対する正味のエネルギ流束は０でなければならない）。これら２つの条件は、各エネルギ準位（つまり、エネルギ準位１及び２）での正味質量流束が０に等しいことを必要とすることにより満たすことができる。

他の２つの境界条件は、ファセットと相互作用する粒子の正味運動量に関連する。本明細書では滑り面として言及する、表面摩擦のない表面に対しては、正味接線方向運動量流束は０でなければならず、正味法線方向運動量流束はファセットでの局所的な圧力に等しくなければならない。従って、受け取った合計運動量と移送した合計運動量の、ファセットの法線ｎ_αに対して垂直な成分（つまり、接線方向成分）は等しくなければならないが、受け取った合計運動量と移送した合計運動量の、ファセットの法線ｎ_αに対して平行な成分（つまり、法線方向成分）の差はファセットでの局所的な圧力に等しくなければならない。非滑り面に対しては、表面の摩擦は、ファセットが受け取る粒子の合計接線方向運動量に対するファセットが移送する粒子の合計接線方向運動量を摩擦量に関連する係数だけ減少させる。

２．ボクセルからファセットへの集まり
粒子と表面間の相互作用をシミュレートする場合の第１のステップとして、粒子は、ボクセルから集めてファセットへ供給する（ステップ３０８）。前述のように、ボクセルＮ（ｘ）とファセットＦ_αとの間の状態ｉ粒子の流束は、以下のようになる。
方程式（１３）

この式から、ファセットＦ_α（ｃ_iｎ_α＜０）に向けられた各状態ｉに対して、ボクセルからファセットＦ_αに供給される粒子数は、以下のようになる。
方程式（１４）

Ｖ_iα（ｘ）が０以外の値を有するボクセルのみを合計する必要がある。前述のように、ファセットのサイズは、Ｖ_iα（ｘ）が少数のボクセルに対してのみ０以外の値を有するように選択される。Ｖ_iα（ｘ）及びＰ_f（ｘ）が非整数値を取る場合があるので、Γ_α（ｘ）は、実数として記憶及び処理される。

３．ファセットからファセットへ移動
次に、粒子をファセット間で移動させる（ステップ３１０）。ファセットＦ_αの流入状態（ｃ_iｎ_α＜０）の平行六面体Ｇ_iαが別のファセットＦ_βと交差する場合、ファセットＦ_αが受け取る状態ｉの粒子の一部は、ファセットＦ_βから到来することになる。特に、ファセットＦ_αは、以前の時間増分中にファセットＦ_βにより生成される状態ｉ粒子の一部を受け取ることになる。この関係を図１０に示すが、ファセットＦ_βと交差する平行六面体Ｇ_iαの部分１０００は、ファセットＦ_αと交差する平行六面体Ｇ_iβの部分１００５に等しい。前述のように、交差する部分は、Ｖ_iα（β）として示す。この項を使用して、ファセットＦ_βとファセットＦ_αとの間の状態ｉ粒子の流束は、以下のように説明することができ、
方程式（１５）
ここで、Γ_i（β，ｔ−１）は、前回の時間増分中にファセットＦ_βによって生成された状態ｉ粒子の尺度である。この式から、ファセットＦ_α（ｃ_iｎ_α＜０）に向けられた各状態ｉに対して、他のファセットによりファセットＦ_αに供給される粒子数は、
方程式（１６）
であり、ファセットへの状態ｉ粒子の全流束は、以下のようになる。
方程式（１７）

ファセット分布関数とも呼ばれるファセットに対する状態ベクトルＮ（α）は、ボクセル状態ベクトルのＭエントリに対応するＭエントリを有する。Ｍは、離散格子速度の数である。ファセット分布関数Ｎ（α）の入力状態は、それらの状態への粒子の流束を容積Ｖ_iαで割ったものに等しく設定され、ｃ_iｎ_α＜０に対して以下のようになる。
方程式（１８）

ファセット分布関数は、ファセットからの出力流束を生成するためのシミュレーションツールであり、必ずしも実際の粒子を表すものではない。正確な出力流束を生成するために、他状態の分布関数に値を割り当てる。外向きの状態は、内向きの状態をポピュレートするために前述の技術を使用してポピュレートされ、ｃ_iｎ_α≧０に対して
方程式（１９）
となり、ここで、Γ_iother（α）は、Γ_iIN（α）を生成するための前述の技術を用いるが、流入状態（ｃ_iｎ_α＜０）以外の状態（ｃ_iｎ_α≧０）にこの技術を適用して決定する。代替技術では、前回の時間ステップからΓ_iOUT（α）の値を用いて、Γ_iother（α）を生成することができるので、
方程式（２０）
となる。

平行状態（ｃ_iｎ_α＝０）に対しては、Ｖ_iα及びＶ_iα（ｘ）は共に０である。Ｎ_i（α）の式において、Ｖｉα（ｘ）は分子に現れ（Γ_iother（α）の式から）、Ｖ_iαは、分母に現れる（Ｎ_i（α）の式から）。従って、平行状態に対するＮ_i（α）は、Ｖ_iα及びＶ_iα（ｘ）が０に近づく時のＮ_i（α）の極限値として決定される。

０速度を有する状態（つまり、静止状態及び状態（０、０、０、２）及び（０、０、０、−２））の値は、温度と圧力に関する初期条件に基づいてシミュレーションの最初に初期化される。これらの値は、その後、時間と共に調整される。

４．ファセットの表面力学の実行
次に、各ファセットが前述の４つの境界条件を満たすように、表面力学を実行する。ファセットに対して表面力学を実行するための手順を図１１に示す。最初に、ファセットＦ_αに垂直な合計運動量は、ファセットにおける粒子の合計運動量Ｐ（α）を確定することによって、以下のように決定され（ステップ１１０５）、全てのｉに対して
方程式（２１）
これから、法線方向運動量Ｐ_n（α）は、以下のように決定される。
方程式（２２）

その後、この法線方向運動量は、Ｎ_n-（α）を生成するために押込／引抜技術（ステップ１１１０）を用いて除去する。この技術に従って、法線方向運動量にのみ影響を及ぼすように状態間で粒子を移動させる。この押込／引抜技術は米国特許第５，５９４，６７１号に記載されており、この特許は、引用により組み込まれている。

その後、Ｎ_n-（α）の粒子を衝突させてボルツマン分布Ｎ_n-β（α）を生成する（ステップ１１１５）。流体力学の実行に関して以下で説明するように、ボルツマン分布は、１組の衝突規則をＮ_n-（α）に適用することによって実現することができる。

その後、ファセットＦ_αに関する流出流束分布を流入流束分布及びボルツマン分布に基づいて決定する（ステップ１１２０）。最初に、流入流束分布Γ_i（α）とボルツマン分布との差を以下のように決定する。
方程式（２３）

この差分を用いると、流出流束分布は、ｎ_αｃ_i＞０に対して
式（２４）
であり、ここで、ｉ^*は、状態ｉと反対の方向を有する状態である。例えば、状態ｉが（１、１、０、０）の場合、状態ｉ^*は、（−１、−１、０、０）である。表面摩擦及び他の要因を考慮するために、流出流束分布を以下のようにさらに精緻化することができる。
ｎ_αｃ_i＞０に対して、
方程式（２５）
ここで、Ｃ_fは、表面摩擦関数であり、ｔ_1αは、ｎ_αに垂直な第１の接線ベクトルであり、ｔ_2αは、ｎ_α及びｔ_1αに垂直な第２の接線ベクトルであり、ΔＮ_j,1及びΔＮ_j,2は、状態ｉ及び表示される接線ベクトルのエネルギ（ｊ）に対応する分布関数である。分布関数は以下の方程式に従って決定する。
方程式（２６）
ここで、ｊは、エネルギ準位１に対して１、エネルギ準位２に対して２に等しい。

Γ_iOUT（α）の方程式の各項の関数は以下の通りである。第１及び第の２項は、衝突がボルツマン分布を生み出すことに有効な範囲で法線方向運動量流束の境界条件を実行するが、接線方向運動量流束偏差を含む。第４及び第５の項は、不十分な衝突による離散効果又は非ボルツマン構造に起因して生ずる場合があるこの偏差を補正する。最後に、第３の項は、表面上の接線方向運動量の所望の変化を強いるために特定の量の表面分画を追加する。摩擦係数Ｃ_fの生成を以下に説明する。ベクトル演算を伴う全ての項は、シミュレーションの開始に先立って計算することのできる幾何形状因子であることに留意されたい。

これにより、接線方向速度を以下のように決定する。
方程式（２７）
ここで、ρは、ファセット分布の密度である。
方程式（２８）

前と同様に、流入流束分布とボルツマン分布との差を以下のように決定する。
方程式（２９）

従って、流出流束分布は、以下のようになる。
方程式（３０）
この方程式は、前回の技術によって求めた流出流束分布の最初の２行に対応するが、特異な接線方向流束に対する補正を必要としない。

あらゆる手法を用いても、結果として生じる流束分布は、運動量流束の条件の全てを満たし、つまり、
方程式（３１）
であり、ここで、ｐ_αは、ファセットＦ_αにおける平衡圧力であり、ファセットに粒子を与えるボクセルの平均密度及び温度値に基づいており、ｕ_αは、ファセットの平均速度である。

質量及びエネルギの境界条件が満たされることを保証するために、入力エネルギと出力エネルギとの差を、以下のように各エネルギ準位ｊに対して測定する。
方程式（３２）
ここで、指標ｊは、状態ｉのエネルギを示す。その後、このエネルギ差を使用して差の項を生成する。
ｃ_jiｎ_α＞０に対して、
方程式（３３）
この差の項は、流束が以下になるように流出流束を修正するために使用する。
ｃ_jiｎ_α＜０に対して、
方程式（３４）
この演算は、質量及びエネルギの流束を補正するが、接線方向の運動量流束を不変のままにする。この調整は、流れがファセットの近傍でほぼ均一で平衡状態に近い場合は小さい。この調整の後に、結果として得られた法線方向運動量流束は、近傍の不均一性又は非平衡特性に起因する補正を加えて近傍の平均特性に基づいて平衡圧力である値に僅かに変更される。

５．ボクセルからボクセルへの移動
再び図３を参照すると、粒子を３次元直線格子に沿ってボクセル間で移動させる（ステップ３１４）。このボクセルからボクセルへの移動は、ファセットとは相互作用しないボクセル（つまり、表面近くに位置しないボクセル）で実行される唯一の移動演算である。典型的なシミュレーションでは、表面と相互作用するほどには近くに位置していないボクセルは、ボクセルの大多数を構成する。

異なる状態の各々は、３次元ｘ、ｙ及びｚの各々において整数速度で格子に沿って移動する粒子を表す。整数速度は、０、±１及び±２を含む。速度の符号は、粒子が対応する軸に沿って移動している方向を示す。

表面と相互作用しないボクセルに関して、移動演算は、計算的にはきわめて簡単である。各時間増分中に、状態の全粒子数を現在のボクセルから目的地ボクセルへ移動させる。同時に、目的地ボクセルの粒子をボクセルから固有の目的地ボクセルへ移動させる。例えば、＋１ｘ方向（１、０、０）に移動しているエネルギ準位１の粒子は、現在のボクセルからｘ方向に＋１、他の方向に対して０のボクセルに移動させる。粒子は、結局、移動（１、０、０）の前に有していたのと同じ状態で目的地ボクセルに行く。ボクセル内の相互作用は、他の粒子及び表面との局所的な相互作用に基づいて、その状態に対して粒子総数を変化させることになる。そうでない場合は、粒子は同じ速度及び同じ方向で格子に沿って移動し続けることになる。

移動演算は、１又は２以上の表面と相互作用するボクセルに対しては少し複雑となる。これにより、１又は２以上の分画粒子は、ファセットに移送される。このような分画粒子のファセットへの移送により分画粒子はボクセル内に残る。これらの分画粒子は、ファセットで占有されるボクセルへ移送される。例えば、図９を参照すると、ボクセル９０５の状態ｉ粒子の部分９００がファセット９１０に移動した場合（ステップ３０８）、残りの部分９１５は、ボクセル９２０へ移動させるが、そこにはファセット９１０が位置しており、状態ｉの粒子はそこからファセット９１０へ向けられる。従って、状態粒子数が２５に等しく、Ｖ_iα（ｘ）が０．２５に等しい場合（つまり、ボクセルの４分の１が平行六面体と交差する）、６．２５個の粒子は、ファセットＦ_αに移動させることになり、１８．７５個の粒子は、ファセットＦ_αにより占有されたボクセルへ移動させることになる。複数のファセットが単一のボクセルと交差することができるので、１又は２以上のファセットにより占有されたボクセルＮ（ｆ）に移送される状態の粒子数は、以下のようになる。
方程式（３５）
ここで、Ｎ（ｘ）は、ソースボクセルである。

６．ファセットからボクセルへの散乱
次に、各ファセットからの流出粒子をボクセルに散乱させる（ステップ３１６）。本質的に、このステップは、粒子をボクセルからファセットへ移動させた収集ステップの逆である。ファセットＦ_αからボクセルＮ（ｘ）へ移動する状態ｉの粒子数は、以下の通りである。
方程式（３６）
ここで、Ｐ_f（ｘ）は、部分的なボクセルの容積減少に対処する。これから、各状態ｉに対して、ファセットからボクセルＮ（ｘ）に向けられる全粒子数は以下のようになる。
方程式（３７）

粒子をファセットからボクセルまで散乱させ、周囲のボクセルから流入した粒子と結合して結果を整数化した後に、特定のボクセルの特定の方向がアンダーフローするか（負になる）、又はオーバーフロー（８ビット実装で２５５を超える）する可能性がある。これは、これらの量が許容範囲の値に適合するように切り捨てられた後に、質量、運動量、及びエネルギの取得又は損失をもたらす。このような事象の発生を防ぐために、範囲を越えた質量、運動量、及びエネルギが、問題のある状態の切り捨てに先立って蓄積される。状態が属するエネルギに関して、得られる（アンダーフローによる）又は失われる（オーバーフローによる）値に等しい量の質量は、同じエネルギを有しかつそれ自体オーバーフロー又はアンダーフローを受けないランダムに（又は連続的に）選択された状態に追加される。この質量及びエネルギの追加から生じる追加の運動量は蓄積され、切り捨てられた運動量に追加される。同じエネルギ状態にだけ質量を追加することにより、質量カウンタが０に到達した場合に質量及びエネルギの両方が補正される。最後に、運動量アキュムレータが０に戻るまで、運動量は、押込／引抜技術を利用して補正される。

７．流体力学の実行
最後に、流体力学を実行する（ステップ３１８）。このステップは、微視的動力学又はボクセル内演算と呼ぶことができる。同様に、移流手順は、ボクセル内演算と呼ぶことができる。また、以下に説明する微視的動力学演算を使用して、ファセットで粒子を衝突させてボルツマン分布を生成することができる。

流体力学は、ＢＧＫ衝突モデルとして知られる特定の衝突演算子によって格子ボルツマン方程式において保証される。この衝突モデルは、実際の流体システムの分布の動力学を模倣している。衝突過程は、方程式１の右辺及び式２によりうまく説明することができる。移流ステップの後で、流体システムの保存量、具体的には、密度、運動量、及びエネルギは、方程式３を用いて分布関数から取得する。これらの量から、方程式（２）のｆ^eqで説明される平衡分布関数は、方程式（４）により完全に特定される。速度ベクトルセットｃ_i及び重みの選択は表１に列挙され、方程式２と一緒に、巨視的挙動が正確な流体力学の方程式に従うことを保証する。

Ｅ．衝突過程
関連する流体物理を再現するために、格子ボルツマンシステムでの衝突過程は、実際の流体システムと同様の基本的な役割を果たし、同様の基本的な保存要件に従う。便宜上、以下の方程式の番号は、方程式１．１）から始まる。
は、「衝突前」粒子分布関数（つまり、時間ｔにおいて位置ｘにある単位容積中の粒子数）とし、この分布関数は衝突後に
に変わる（つまり、「衝突後」分布関数）。質量の保存は以下の方程式を満たす。
ここで、ρ（ｘ、ｔ）は、位置ｘ及び時間ｔでの全速度ベクトル値の粒子数密度に等しい流体密度である。（１．１）（及び以下の方程式）の合計は、格子ボルツマンモデルにおける可能な粒子速度ベクトル値の全てに対するものである。運動量の保存は、
によって与えられ、ここでｕ（ｘ，ｔ）は流体速度であり、単に、位置ｘ及び時間ｔにおける粒子間の平均速度である。

粒子の運動エネルギも保存される特定の流体システム（多粒子で構成される）に関しては、（１．１）及び（１．２）に加えて、次の関係も定義される。
ここで、
であり、定数Ｄは粒子運動の次元である。Ｔ（ｘ，ｔ）は、ｘ及びｔにおける流体システムの温度である。

保存則から明らかであるが、ρ（ｘ、ｔ）及びｕ（ｘ，ｔ）（及びエネルギ保存システムではＴ（ｘ、ｔ））は、衝突過程中には不変である。所定のρ（ｘ，ｔ）及びｕ（ｘ，ｔ）（及びＴ（ｘ，ｔ））に対して、平衡分布関数とも呼ばれる、特別な形式の粒子分布関数
が存在する。平衡分布関数は、式（１．１）及び（１．２）（及び（１．３））に定義されるものと同じ質量及び運動量（及びエネルギー）を有し、ρ（ｘ，ｔ）及びｕ（ｘ，ｔ）（及びＴ（ｘ，ｔ））の関数として完全に決定される。

分布関数の合計（粒子速度ベクトル値に対する）に関する量は、一般に、分布関数のモーメントと呼ばれる。質量、運動量、及びエネルギの巨視的保存量に対応する、（１．１）−（１．３）の３つの基本的なモーメントの他に、それらの流束に対応するモーメントは、同様に重要である。これらは、衝突後の分布関数の合計によって定めることができる。
ここで、Π（ｘ，ｔ）及びＱ（ｘ，ｔ）はそれぞれ、位置ｘ及び時間ｔにおける運動量及びエネルギの流束テンソルである。（１．２）の定義を想起されたい。運動量流束及びエネルギ流束だけが上記の巨視的保存量（（１．１）−（１．３）で定義される）に無関係である。

平衡状態での分布に関して（つまり、（１．５）及び（１．６）で代わりに
を用いる）、平衡状態の運動量流束及びエネルギ流束は、それぞれ基本的な（連続）気体分子運動理論から公知の形態
を有し、ここで、ｐ（ｘ，ｔ）は、圧力であり、理想気体の法則に基づいてｐ（ｘ，ｔ）＝ρ（ｘ，ｔ）Ｔ（ｘ，ｔ）であり、（１．７）の「Ｉ」は２階の単位テンソルを示す。Ｔｒ［ ］はトレース演算である。格子ボルツマン法での中心テーマは、（１．７）及び（１．８）の形態を取り戻すことである。

（１．５）及び（１．６）のモーメントは、平衡部及び非平衡部の項で以下のように表すこともできる。

明らかに、
である。非平衡運動量流束は、流体システムでの輸送挙動を決定する際に重要な役割を果たす。例えば、非平衡運動量流束（１．９）は、ニュートン流体システムにおいて粘性を直接決定するが、（１．１０）の非平衡エネルギ流束は、温度伝導率を決定する。粒子分布関数から無数のモーメントを構成することができるが、（１．１）−（１，１０）によるモーメントだけが物理的な流体システムにおいて巨視的に関連する。

本システムは、平衡状態からの偏差の程度を
として示し、従って、物理的に安定な衝突過程は、常にこの偏差を低減する方向に働く。すなわち、平衡状態からの衝突後の偏差は、衝突前の偏差よりも小さい。

実際の多粒子システムでの衝突は非常に複雑な場合がある。１つの最も簡単な衝突モデル（演算子）は、保存則要件と（１．１１）の平衡状態への収束要件との両方を満たす。「ＢＧＫ」衝突過程は、数学的には以下のように表され、
衝突演算子は、以下の通りである。

３つの分布関数の全てが同じρ（ｘ，ｔ）及びｕ（ｘ，ｔ）（及びエネルギ保存システムに対してはＴ（ｘ，ｔ））を与えるので、質量及び運動量（及びエネルギー）の保存（つまり、（１．１）−（１．３））は、自動的に満たされる。さらに、（１．１２）から、衝突後分布関数に関する平衡状態からの偏差は、衝突前分布関数に関する平衡状態からの偏差に係数（１−１／τ）で比例する。従って、収束条件（１．１１）は、パラメータ値τ（衝突緩和時間と呼ぶ）が１／２よりも大きい限り満たされ、衝突後偏差は、τ＝１の時にゼロになる。

ＢＧＫ関係式（１．１２）及び（１．１３）から、方程式（１．９）及び（１．１０）は、以下のように書き換えることができ、
ここで、Π^neq（ｘ，ｔ）及びＱ^neq（ｘ，ｔ）はそれぞれ、以下に定義される衝突前の非平衡状態の運動量流束及びエネルギ流束である。

τの値が選択された状態で、ＢＧＫ衝突演算子（１．１２）を用いて格子ボルツマン流体における動粘性率値は、以下のように決定される。

ここで、定数Ｔ₀は標準格子温度であり、格子単位系での時間増分が１となるように、
いわゆる格子単位規則を用いる。ＢＧＫ衝突演算子は、格子ボルツマンモデルにおいて最も一般的に用いられる演算子である。ＢＧＫ衝突演算子は、過去２０数年で種々レベルの成功を収めた。一方で、ＢＧＫ演算子にはいくつかの本質的な制限がある。単一プラントル数（＝粘性率と熱拡散率の比）の他に、（１．９）−（１．１０）の基本的な問題に加えて、１つの問題は、τが１でない場合は、全ての衝突後非平衡モーメントが生成される点にある。実際にＢＧＫに関して、衝突前後の非平衡モーメントは、以下の関係を示し、
ここで、Ｍ_n′（ｘ，ｔ）及びＭ_n（ｘ，ｔ）はそれぞれ、衝突前及び衝突後の任意のｎ次モーメントを表す。

Ｆ．フィルタ衝突演算子(Filtered Collision Operator)
任意の格子ボルツマンモデルの最も一般的な特徴は、有限かつ一定の粒子速度ベクトル値セットを有することである。所定の格子ボルツマンモデルに対して、定数ベクトル値セットが特定される。その結果、離散速度セットから構成された粒子分布関数のモーメントの有限セットだけが、実際の流体におけるそれらの対照物を取り戻すことができる。実際の流体のモーメントを任意の所定の次数まで取り戻すための一般的な構成は、厳密に規定される。格子ボルツマンベクトル値の異なるセットは、物理的モーメントを異なる次数までサポートすることができる。例えば、いわゆるＤ３Ｑ１９及びＤ３Ｑ１５は、平衡運動量流束と低マッハ数の制限下での線形偏差の（ニュートンの）非平衡運動量流束までのモーメントをサポートするだけである。一方では、Ｄ３Ｑ３９等のいわゆる高次格子ボルツマンモデルは、有限マッハ数及びクヌーセン数での線形偏差を超える平衡エネルギ流束及び非平衡運動量流束までモーメントをサポートすることができる。

物理的に関連するモーメントは、平衡及び非平衡運動量、並びにエネルギのモーメント及びそれらの流束だけなので、正確な物理的流体の挙動を実現するためには、非平衡運動量流束（場合によってはエネルギー流束）以外の全ての衝突後非平衡モーメントがゼロになるように、衝突演算子を設計することが望ましい。具体的には、非平衡運動量流束及びエネルギ流束は、以下の方程式で表されるが、
他の衝突後非平衡モーメントはゼロである。

（１．１９）のＢＧＫ衝突は、全ての非平衡モーメントを生成することを想起されたい。また、（２．１）−（２．３）を実現する衝突演算子は、非本質的な非平衡モーメントをフィルタ除去するので、フィルタ衝突演算子と呼ばれる。フィルタ衝突演算子は、他の非平衡モーメントを全て除去しながら、非平衡運動量流束を（２．１）のように維持する。このフィルタ衝突演算子の明示的表現は以下の通りである。
ここで、定数ｗ_iは、特定の格子速度セットが選択された状態で決まる重み係数値である。格子速度セットは、格子に制限される空間での微視的速度又は運動量の離散セットである。重み係数値のセットは、異なる格子ボルツマン粒子速度セットに対して異なり、その目的は、期待される次数までモーメントの等方性を実現することである。新しい（フィルタ処理された）衝突演算子（２．４）は、ＢＧＫと同じ衝突後非平衡運動量流束（１．９）をもたらし、自動的にＢＧＫと同じ（１．１８）の粘性作用を引き起こす。

特定の格子システムに関して、１ではないτ_eを有する方程式（２．２）に示すように、非ゼロの非平衡エネルギ流束を得ることが望ましい。このような目的を実現するために示された特定の衝突演算子の形式は以下の通りである。

値τ_eは、τと等しい必要はないので、プラントル数は、ＢＧＫの場合とは対照的に単一の値に限定されないことに留意されたい。運動量流束とエネルギ流束との間のモーメント直交性により、より一般化されたフィルタ衝突演算子形式は（２．４）及び（２．５）の２つの形式の直和により構成されるので、（２．１）及び（２．２）の非平衡流束は、両方とも自動的に満たされるが、残りの非平衡モーメントはゼロになる。

所定の格子速度セットではサポートされない不必要な非平衡モーメントをフィルタ除去することによって、一般化された衝突演算子（（２．４）、（２．５）又はこれらを組み合わせたもの）は、ＢＧＫ衝突演算子よって著しく改善された流体流れ等方性を示し、ＢＧＫ衝突演算子と同様に所望の運動量流束及びエネルギ流束を維持する。一般化された衝突演算子の形式（２．４）及び（２．５）（及びこれらを組み合わせたもの）の態様は、ナビエ−ストークス流体の粘性率及び熱拡散率に適用可能であるだけでなく、有限なクヌーセン数を含むより広範な流体レジームにおいて正確な流体力学を保証する。

Ｇ．Ｎ次のガリレイ不変フィルタ演算子を導出する手順
本開示に適合するシステムは、どれだけの量の速度がシミュレーションによってサポート可能かに基づいて、格子マッハ数の特定の指数を有するフィルタ衝突演算子を生成する。このフィルタ衝突演算子では、絶対値（例えば、速度又はエネルギの）ではなく相対値（例えば速度又はエネルギーの）を使用する。本明細書では、不必要なモーメントを含まない衝突演算子に対する正しい理論形式を適切に構成する手順を説明する。非平衡モーメントの総数は、シミュレーションによってサポート可能な速度の総数に対応する。従って、フィルタ衝突演算子は、所定の格子ボルツマンモデル粒子速度セットによってサポートされる非平衡モーメントを表すので、非平衡モーメントは物理的世界で実際に起こるものに対応する。相対速度を用いることにより、シミュレーションの速度によってサポートされる所望の非平衡状態の寄与部を備えたフィルタ衝突演算子の構成が可能となる。

すなわち、非平衡状態の運動量流束及びエネルギ流束を相対的な粒子速度及びエネルギを用いて決定することにより、シミュレーションにおいて格子ボルツマンモデルの格子速度セットによりサポートされない、より高次の項をフィルタ除去する（例えば、除外する）フィルタ衝突演算子を構成することができる。１つの実施例では、フィルタ衝突演算子は拡張形式で特定されるので、多数の格子速度に関して、ユーザ又はシステムは、どの次数までの項を保存する必要があり、どの次数を超える項を除去する必要があるか（フィルタ衝突演算子から）を知る方法を有する。そうするために、拡張部（例えば、フィルタ衝突演算子）を相対速度で表す必要がある。

１つの実施例では、速度モデルは、粒子速度の有限セットである。従って、粒子移動の本当の物理モーメントだけが、特定の次数までこのモデルで正確に表すことができる。ガリレイ不変衝突演算子を有しかつより速い流速で速度をシミュレートするために、衝突演算子の形式は相対速度に基づいており、ここでは粒子速度は自身の流速に対して測定される。少なくともこの理由から、パラメータｃ_iを以下の方程式３．５で使用するのではなく、パラメータｃ_i’（相対速度又は相対エネルギーを表す）を方程式３．５で使用する。方程式３．５でのパラメータｃ_iの使用は、流速（例えば、ｕ（ｘ，ｔ））の累乗でこの衝突演算子形式を展開することを許さないので、局所的な流速に対する相対速度に関してガリレイ不変性に従う適切な衝突演算の構成をもたらさない。しかしながら、方程式３．５でのパラメータｃ_i’の使用は、流速（例えば、ｕ（ｘ，ｔ））の累乗（例えば、Ｎ次）でこの衝突演算子形式を展開することを可能にするので、Ｎ次のガリレイ不変フィルタ演算子をもたらす。

１つの実施例では、ＢＧＫモデルの下での衝突演算子は、無限次数の非平衡モーメントを含む。特定の速度モデルは無限次数まで正確にモーメントをサポートできないので、特定の次数までの非平衡モーメントは、非物理的なものとなる。従って、このシステムは、高次の非平衡モーメントを除去して物理的人工物を含まないようにする必要がある。例えば、１９速度モデルに関連するモーメントは、一次までである。従って、より高次のモーメントは不適切である。しかしながら、ＢＧＫ演算子を用いて研究シミュレーションを実行する場合、モデルによりサポートされるものを超える非平衡モーメントが含まれる。従って、本明細書に説明するフィルタ衝突演算子は、物理的世界で起こっているものにより密接に対応する。

前述のフィルタ衝突演算子形式は、非常に低い流速の流れ状況に対応するのみである。ガリレイ不変性の基本原理から、物理的流体の統計的特性の全ては、平均流束に対する粒子相対速度だけの関数である必要がある。具体的には、物理的流体システムにおける粒子非平衡分布関数及び関連の非平衡モーメントは、絶対速度ｃ_i（格子ボルツマンシステムのゼロフロー基準座標系で測定される）ではなく、相対速度（ｃ_i−ｕ（ｘ，ｔ））にのみ依存するべきである。従って、（１．１４）及び（１．１５）の代わりに、より物理的に意味のある非平衡の運動量流束及びエネルギ流束は、それぞれ以下の通りである。
ここで、粒子相対速度及びエネルギは、
で与えられる。上記方程式（３．３）に示すように、粒子絶対速度は、相対速度で置換される。興味深いことに、質量及び運動量の保存により、非平衡運動量流束（３．１）は（１．１６）と同じであることが分かる。一方で、非平衡エネルギ流束（３．２）は、（１．１７）に変換できない。

本開示に適合するシステムは、非平衡運動量流束に関するガリレイ不変衝突演算子を生成する。気体分子運動論の基本法則に従って、流速の不均一性に起因する、１次（ｌｅａｄｉｎｇ ｏｒｄｅｒ）の非平衡粒子分布関数は、以下の通りである。
ここで、Ｕは、連続体気体分子運動論での局所的平均流れに対する粒子相対速度である。この概念に端を発して、本明細書に説明したシステムは、ボルツマン法において、ガリレイ不変性に一致する非平衡分布関数に類似した式を特定する。本システムは、完全に対応する非平衡衝突後分布関数として、以下の明示形式を特定する。
ここで、Π^neq（ｘ，ｔ）は、（１．１４）又は（３．１）で与えられる。上記の方程式３．２に示すように、分布関数は、平衡成分及び非平衡成分を含む。さらに、上記の方程式は、粒子速度セットを無限次数まで必要とする。（３．５）の非平衡分布関数は、任意のマッハ数に対する完全形式である。

それぞれ非平衡及び平衡の分布関数に対する（３．５）及び（３．６）の完全形式は、格子速度セットが全次数まで正確な流体力学モーメントをサポートする場合にのみ実現可能である。換言すると、粒子速度値の有限セットを備える任意の所与の速度セットに対して、完全形式は実現できない。しかしながら、以下に詳細に説明するように、これらを所定の有限格子速度セットを用いて対応する次数Ｎまで実現することが可能である。

式（３．５）及び（３．６）は、ｕ（ｘ，ｔ）の累乗で多項式形式に展開することができる。平衡分布関数（３．６）は、エルミート多項式級数を用いて表現できることはとはよく知られている。
ここで、Ｖ（ｘ，ｔ）^[n]は、ベクトルＶ（ｘ，ｔ）ｎ次直積に関する簡略表記である。
ｎ次エルミート多項式Ｈ⁽ⁿ⁾（ξ_i）は、標準（スカラー）ｎ次エルミート関数のｎ階テンソル一般化である。裏付けはないが、非平衡の寄与分（３．５）は、エルミート多項式級数で表現することもできる。

エルミート多項式の直交特性を用いると、ｍ≦Ｎに対してｕ^m（ｘ、ｔ）に比例する項を維持しながら、２つの無限級数（３．７）及び（３．８）を単純に切り捨てることによって、システムは、これらの形式に対するＮ次近似を得る。１つの速度モデルは、各格子点をその第１及び第２近傍と連結する１９速度のＤ３Ｑ１９立方格子である。具体的には、Ｄ３Ｑ１９（又はＤ３Ｑ１５）型の格子速度セットに対して、システムは（３．６）のｆ_i ^eq（ｘ，ｔ）をｕ³（ｘ，ｔ）まで、（３．５）のＣ_i（ｘ，ｔ）をｕ（ｘ，ｔ）の第１（線形）累乗まで展開する。

上記の方程式３．９に示すように、この衝突演算子は、１９速度をサポート可能なシミュレーションシステムに対して正確である。さらに、この衝突演算子は、１９速度をサポート可能なシミュレーションに対して物理的世界で起こるものに対応するように、（３．５）のＣ_i（ｘ，ｔ）を第１（線形）累乗まで展開した形式である。この衝突演算子は、方程式１．１２で使用して分布関数を修正することができる。上記の方程式３．９において、ｘは容積内の特定位置、ｔは時間における特定時点、ｉは格子速度セットでの格子速度の指数、Ｔ₀は一定格子温度、ｃ_iは衝突前の粒子速度ベクトル、ｕ（ｘ，ｔ）は時間ｔ及び特定位置ｘにおける各粒子間の平均速度、Ｉは２階の単位テンソル、τは衝突緩和時間、ｗ_iは一定重み係数、及びΠ^neqは非平衡運動量流束である。

別の速度モデルはＤ３Ｑ３９格子であり、これは３９速度までサポートする。Ｄ３Ｑ３９に関して、展開は、ｕ（ｘ，ｔ）の２乗又は３乗まで行うことができる。ｕ²（ｘ，ｔ）まででの切り捨ては、以下に明示する通りである。

上記の方程式３．１０に示すように、衝突演算子の展開形式は、２乗を超える、従ってシミュレーションでサポートされない高次のモーメントをその展開形式中に含まないことにより、不必要なモーメントを除外する。本実施例では、ｘは容積内の特定位置、ｔは時間における特定時点、ｉは格子速度セットでの格子速度の指数、Ｔ₀は一定格子温度、ｃ_iは衝突前の粒子速度ベクトル、ｕ（ｘ，ｔ）は時間ｔ及び特定位置ｘにおける各粒子間の平均速度、Ｉは２階の単位テンソル、τは衝突緩和時間、ｗ_iは一定重み係数、及びΠ^neqは非平衡運動量流束である。

従って、システムは、（３．５）からｕ（ｘ，ｔ）と無関係の項だけを含むので、（２．４）の衝突演算子形式を０次近似（つまり、（３．５）のＣ_i ⁽⁰⁾（ｘ，ｔ））として再解釈する。Ｄ３Ｑ３９等のより高次の格子速度セットに関しては、ｕ⁵（ｘ，ｔ）までの平衡分布関数（３．６）の項を維持できるが、ｕ³（ｘ，ｔ）までの非平衡分布関数（３．５）の項は維持できない。一般に、格子速度セットが流体力学モーメントに十分に対応する次数のサポートを提供する場合、この体系的な手順は、任意次数（ｕ（ｘ，ｔ）の累乗）まで行うことができる。所与の有限次数においてガリレイ不変性は厳密には満たされないが、より高い格子速度セットを使用する場合、及び（３．５）及び（３．６）においてより高次の展開形式を用いる場合、誤差はより高次に向かう。

上記では、流体速度の累乗での任意の所定の次数に対する一般化されたフィルタ衝突演算子を構成するために、システムは、体系的な手順を生成して実行する。特に、衝突後非平衡流束に関しては、既存のフィルタ演算子（２．４）に対する１次及び２次の補正は、（３．９）及び（３．１０）に明示的に示される。フィルタ衝突演算子は、所望のモーメント以外の非平衡モーメントをフィルタ除去する。フィルタ演算子（２．４）及び（３．９）（又は（３．１０）は、非平衡運動量流束を維持する目的を果たす。一方で、（３．９）（及び（３．１０））は、ＢＧＫ及び（２．４）と同じ非平衡運動量流束を提供するが、この衝突演算子は、改善された数値的安定性及びガリレイ不変性を実現する。

同じ手順をエネルギ流束に関係するフィルタ衝突演算子に対しても策定することができる。一般的な完全ガリレイ不変の形態は、（３．５）に類似しており、相対速度で表される。
ここで、Ｗ^neq（ｘ，ｔ）は、（３．１）及び（３．２）のＱ^neq（ｘ，ｔ）及びｕ（ｘ，ｔ）Π^neq（ｘ，ｔ）の適切な線形結合である。同じ手順に従って、所定の格子速度セットによって十分にサポートされるｕ（ｘ，ｔ）の累乗での任意の有限次数の形式は、体系的に得ることができる。運動量流束とエネルギ流束との間のモーメント直交性により、所望の粘性率及び熱拡散率を独立して実現できる一般的な衝突後演算子形式は、（３．５）及び（３．１１）の追加として（適切な展開形式で）簡単に生成される。上記の方程式（３．１１）では、ｘは容積内の特定位置、ｔは時間における特定時点、ｉは格子速度セットでの格子速度の指数、Ｔ₀は一定格子温度、Ｉは２階の単位テンソル、τは衝突緩和時間、ｃ_i’（ｘ，ｔ）は粒子相対速度、ｆ_i ^eqは平衡分布関数、及びＷ^neqは非平衡エネルギ流束である。

図１２を参照すると、本開示に適合するシステムは、非平衡衝突後分布関数を決定するために処理１２００を実行する。作動時、本システムは、粒子速度の特定の次数まで流体力学運動をサポートする格子速度セットを提供する（そうでない場合には取得する）（１２０２）。例えば、本システムは、１次までの流体力学運動の１９速度をサポートするＤ３Ｑ１９を取得する。本実施例では、格子速度セットに関するサポートされる次数（例えば１９速度）は、非平衡衝突後分布関数の特定の次数より小さくこれとは異なり（例えば、１次又は線形）、非平衡衝突後分布関数の特定の次数は、粒子速度の次数によって決定される。すなわち、１９速度をサポートする速度モデルに対して、非平衡衝突後分布関数は線形であるか又は１次までサポートされるだけである（１９速度は、１次の流速（つまり、ｕ（ｘ，ｔ））と関係するので）。本実施例において、本システムは、特定の速度モデルに対して流速のサポートされる次数にアクセスするように構成される。本実施例において、本システムは、この情報を有する一覧表又は別の関数にアクセスする。

１つの実施例では、格子速度セットは、格子ボルツマン法と関係する状態ベクトルセットである。本実施例において、状態ベクトルは、格子部位での粒子間の衝突挙動を表す一連の２進ビットである。別の実施例において、格子速度セットは、格子に限定される空間で運動量状態のセットから成る。

本システムは、格子速度セットにおいて、移動が粒子間の衝突を引き起こす、流体容積中の粒子の移動をシミュレートする（１２０４）。シミュレーション過程は前述の通りである。シミュレートされた移動に基づいて、容積内の特定位置における粒子の相対速度を決定し（１２０６）、粒子相対速度は、（ｉ）容積内の特定位置において容積の流れが０の下で測定される粒子の絶対速度と（ｉｉ）容積内の特定位置における１又は２以上の粒子の平均速度との間の差である。１つの実施例では、粒子相対速度は、容積内の特定位置において容積の流れが０の下で測定される粒子の絶対速度から容積内の特定位置における１又は２以上の粒子の平均速度を差し引いたものを含む。１つの実施例では、容積内の特定位置における１又は２以上の粒子の平均速度は、その特定位置における特定タイプの粒子の平均速度から成ることが含まれる。例えば、流体容積には多様な異なるタイプの粒子が含まれる場合がある。本実施例において、本システムは、少なくとも特定タイプの粒子から成る部分セットの平均速度を決定するように構成される。前述のように、本システムは、式（３．３）に基づいて粒子の相対速度を決定する。

また、本システムは、粒子の相対速度に基づいて、衝突を表す特定の次数の非平衡衝突後分布関数を決定する（１２０８）。１つの実施例では、非平衡衝突後分布関数は、（ｉ）予め定義された物理量に対する非平衡モーメントを維持し、（ｉｉ）未定義の物理量に対する非平衡モーメントを特定の次数まで除去する。本実施例では、非平衡衝突後分布関数は、非平衡衝突後分布関数の展開形式においてこれらの予め定義された物理量を表す項を含むことにより、予め定義された物理量に対する非平衡モーメントを維持する。非平衡衝突後分布関数は、その展開形式において未定義の物理量を表す項を含まないことにより、例えば、方程式（３．７）、（３．８）に示される無限級数を切り捨てて特定の次数に比例する項を維持することにより、未定義の物理量に対する非平衡モーメントを特定の次数まで除去する。本実施例では、特定の次数は、流体速度の格子音速に対する比に関係した指数値であり、格子速度セットは、指数値をサポートする。

図１３は、ネットワーク環境１３００の構成要素のブロック図である。ネットワーク環境１３００にはさらにシステム１３０２が含まれ、そのシステムにはメモリ１３０４、バスシステム１３０６、及びプロセッサ１３０８が含まれる。メモリ１３０４は、ハードドライブと、ダイナミック・ランダム・アクセス・メモリ、機械可読ハードウェア記憶装置、機械可読媒体又は他タイプの固定機械可読記憶装置などのランダムアクセスメモリとを備えることができる。例えば、データバス及びマザーボードを含むバスシステム１３０６は、システム１３０２の構成要素間のデータ通信を確立し制御するために使用することができる。プロセッサ１３０８は、１又は２以上のマイクロプロセッサ及び／又は処理装置を含むことができる。一般に、プロセッサ１３０８は、データの受信と記憶、及びネットワーク上の通信（図示せず）が可能な任意の適切なプロセッサ及び／又は論理回路を含むことができる。

システム１３０２は、サーバ、分散型計算システム、デスクトップ型コンピュータ、ラップトップ、携帯電話機、ラックマウント・サーバなどのデータ受信が可能な種々の計算装置のいずれでもよい。システム１３０２は、単一のサーバとすること、或いは同じ場所に又は別の場所にある一群のサーバとすることができる。図示するシステム１３０２は、入力／出力（「Ｉ／Ｏ」）インタフェース１３１０を経由してデータを受信することができる。Ｉ／Ｏインタフェース１３１０は、イーサネット（登録商標）、無線ネットワーク・インタフェース、光ファイバ・ネットワーク・インタフェース、モデムなどの、ネットワーク上のデータを受信可能な任意タイプのインタフェースとすることができる。システム１３０２は、速度モデル、シミュレーションデータなどを記憶するように構成することができるデータ格納部１３１２との通信のために設定される。

本明細書に説明する技術を用いて、システムは、非平衡衝突後分布関数、例えばガリレイ不変フィルタ演算子を生成するように説明される。これらの技術を用いて、種々のタイプの非平衡衝突後分布関数は、例えば、方程式（３．５）、（３．９）、（３．１０）及び（３．１１）に示すように生成される。生成された非平衡衝突後分布関数は、流体容積中の粒子衝突過程をモデル化する際に、例えば方程式（１．１２）に示すように使用する

実施形態は、デジタル電子回路、又はコンピュータハードウエア、ファームウェア、ソフトウェア、又はその組合せに実施することができる。本明細書に説明する技術の装置は、プログラマブルプロセッサによる実行のために機械可読媒体（例えば、ストレージデバイス）に有形に具現化又は記憶されたコンピュータプログラム製品に実施することができ、方法アクションは、入力データに対して演算して出力を生成することにより、本明細書に説明する技術の演算を実行する命令のプログラムを実行するプログラマブルプロセッサによって実行することができる。本明細書に説明する技術は、データ及び命令をデータストレージシステム、少なくとも１つの入力デバイス、及び少なくとも１つの出力デバイスから受信して、データ及び命令をそれらに送信するように結合された少なくとも１つのプログラマブルプロセッサを含むプログラマブルシステム上で実行可能である１つ又はそれよりも多くのコンピュータプログラムで実行することができる。各コンピュータプログラムは、高レベル手順又はオブジェクト指向プログラミング言語、又は必要に応じてアセンブリ又はマシン語で実行することができ、いずれにしても、言語は、コンパイル又は解釈された言語とすることができる。

適切なプロセッサには、一例として、汎用及び専用マイクロプロセッサの両方が挙げられる。一般的に、プロセッサは、命令及びデータを読取専用メモリ及び／又はランダムアクセスメモリから受信することになる。一般的に、コンピュータは、データファイルを記憶する１つ又はそれよりも多くの大容量ストレージデバイスを含むことになり、そのようなデバイスには、内蔵ハードディスク及び取外し可能ディスクのような磁気ディスク、光磁気ディスク、光ディスクが挙げられる。コンピュータプログラム命令及びデータを有形に具現化するのに適するストレージデバイスには、一例として、ＥＰＲＯＭ、ＥＥＰＲＯＭ、及びフラッシュメモリデバイスのような半導体メモリデバイス、内部ハードディスク及び着脱式ディスクのような磁気ディスク、光磁気ディスク、及びＣＤ−ＲＯＭディスクを含む全ての形態の不揮発性メモリが挙げられる。上述のいずれも、ＡＳＩＣ（特定用途向け集積回路）によって補足されるか、又はＡＳＩＣに組み込むことができる。

いくつかの実施を説明した。それにも関わらず、本発明の精神及び範囲から逸脱することなく様々な修正を行うことができることは理解されるであろう。従って、他の実施も以下の特許請求の範囲内である。

１３００ ネットワーク環境
１３０２ システム
１３０４ メモリ
１３０６ バスシステム
１３０８ 処理装置
１３１０ インタフェース
１３１２ データ格納部

Claims (48)

1. 格子速度セットにおいて流体容積中の粒子の移動をシミュレートするステップであって、前記移動が前記粒子間の衝突を引き起こすものである、ステップと、
   前記シミュレートされた移動に基づいて、
   前記容積内の特定位置における粒子の粒子相対速度を決定するステップであって、前記粒子相対速度が（ｉ）前記容積内の前記特定位置において前記容積の流れがゼロの下で測定される前記粒子の絶対速度と（ｉｉ）前記容積内の前記特定位置における１又は２以上の前記粒子の平均速度との間の差である、ステップと、
   前記粒子相対速度に基づいて、前記衝突を表す特定の次数の非平衡衝突後分布関数を決定するステップと、
   を含む方法。
2. １又は２以上のコンピュータシステムによって粒子速度の或る次数まで流体力学運動をサポートする格子速度セットを提供するステップをさらに含み、前記シミュレートするステップは、前記１又は２以上のコンピュータシステムによってシミュレートするステップである、請求項１に記載の方法。
3. 前記格子速度セットに関する前記サポートされる次数は、前記非平衡衝突後分布関数の前記特定の次数より小さくかつこれとは異なり、
   前記非平衡衝突後分布関数の前記特定の次数は、前記粒子速度の前記次数によって決定される、請求項２に記載の方法。
4. 前記容積内の前記特定位置における前記１又は２以上の前記粒子の前記平均速度は、前記特定位置における特定タイプ粒子の平均速度を含む、請求項１に記載の方法。
5. 前記格子速度セットは、格子ボルツマン法に関連する状態ベクトルのセットである、請求項１に記載の方法。
6. 前記非平衡衝突後分布関数は、（ｉ）予め定義された物理量に対する非平衡モーメントを維持し、（ｉｉ）未定義の物理量に対する非平衡モーメントを前記特定の次数まで除去する、請求項１に記載の方法。
7. 前記特定の次数は、前記流体速度の格子音速に対する比率に関連した指数値であり、前記格子速度セットは、前記指数値をサポートする、請求項１に記載の方法。
8. 前記格子速度セットは、格子に限定される空間における運動量状態のセットを含む、請求項１に記載の方法。
9. 前記粒子相対速度は、前記容積内の前記特定位置において前記容積の流れがゼロの下で測定される前記粒子の前記絶対速度から前記容積内の前記特定位置における前記１又は２以上の前記粒子の前記平均速度を差し引いたものである、請求項１に記載の方法。
10. 前記非平衡衝突後分布関数は、ガリレイ不変フィルタ処理演算子である、請求項１に記載の方法。
11. 前記非平衡衝突後分布関数に基づいて、前記流体容積中での前記粒子の衝突過程をモデル化するステップをさらに含む、請求項１に記載の方法。
12. 前記非平衡衝突後分布関数は、流体力学モーメントに対して１次のサポートを提供する格子速度セットに対するマッハ数に関する１次のガリレイ不変衝突演算子
    であり、前記衝突演算子は、
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    ｃ_iは衝突前の前記粒子の速度ベクトルであり、
    ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける前記粒子間の平均速度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    ｗ_iは一定重み係数であり、
    は非平衡運動量流束である、請求項１に記載の方法。
13. 前記非平衡衝突後分布関数は、流体力学モーメントに関して無限次数のサポートを与える、格子速度セットに対する衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    は粒子相対速度であり、
    ρは流体密度であり、
    は平衡分布関数であり、
    は非平衡運動量流束である、請求項１に記載の方法。
14. 前記非平衡衝突後分布関数は、流体力学モーメントに関して２次のサポートを提供する、格子速度セットに対するマッハ数に関して２次のガリレイ不変衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    ｃ_iは衝突前の前記粒子の速度ベクトルであり、
    ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける前記粒子間の平均速度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    ｗ_iは一定重み係数であり、
    は非平衡運動量流束である、請求項１に記載の方法。
15. 前記予め定義された物理量は、前記特定の容積中の前記流体の質量、前記特定の容積中の前記流体の運動量、及び前記特定の容積中の前記流体のエネルギを含む、請求項１に記載の方法。
16. 前記非平衡衝突後分布関数は、エネルギ流束に関連する衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    は相対的な粒子速度であり、
    は平衡分布関数であり、
    は非平衡エネルギ流束である、請求項１に記載の方法。
17. オペレーションを実行する１又は２以上の処理装置で実行可能な命令を記憶する１又は２以上の機械可読ハードウェア記憶装置であって、前記オペレーションは、
    格子速度セットにおいて流体容積中の粒子の移動をシミュレートするステップであって、前記移動が前記粒子間の衝突を引き起こすものである、ステップと、
    前記シミュレートされた移動に基づいて、
    前記容積内の特定位置における粒子の粒子相対速度を決定するステップであって、前記粒子相対速度が（ｉ）前記容積内の前記特定位置において前記容積の流れがゼロの下で測定される前記粒子の絶対速度と（ｉｉ）前記容積内の前記特定位置における１又は２以上の前記粒子の平均速度との間の差である、ステップと、
    前記粒子相対速度に基づいて、前記衝突を表す特定の次数の非平衡衝突後分布関数を決定するステップと、
    を含む１又は２以上の機械可読ハードウェア記憶装置。
18. 前記オペレーションは、粒子速度の或る次数まで流体力学運動をサポートする格子速度セットを提供するステップをさらに含む、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
19. 前記格子速度セットに関する前記サポートされる次数は、前記非平衡衝突後分布関数の前記特定の次数より小さくかつこれとは異なり、
    前記非平衡衝突後分布関数の前記特定の次数は、前記粒子速度の前記次数によって決定される、請求項１８に記載の１又は２以上の機械可読ハードウェア記憶装置。
20. 前記容積内の前記特定位置における前記１又は２以上の前記粒子の前記平均速度は、前記特定位置における特定タイプ粒子の平均速度を含む、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
21. 前記格子速度セットは、格子ボルツマン法に関連する状態ベクトルのセットである、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
22. 前記非平衡衝突後分布関数は、（ｉ）予め定義された物理量に対する非平衡モーメントを維持し、（ｉｉ）未定義の物理量に対する非平衡モーメントを前記特定の次数まで除去する、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
23. 前記特定の次数は、前記流体速度の格子音速に対する比率に関連した指数値であり、前記格子速度セットは、前記指数値をサポートする、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
24. 前記格子速度セットは、格子に限定される空間における運動量状態のセットを含む、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
25. 前記粒子相対速度は、前記容積内の前記特定位置において前記容積の流れがゼロの下で測定される前記粒子の前記絶対速度から前記容積内の前記特定位置における前記１又は２以上の前記粒子の前記平均速度を差し引いたものである、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
26. 前記非平衡衝突後分布関数は、ガリレイ不変なフィルタ処理演算子である、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
27. 前記非平衡衝突後分布関数に基づいて、前記流体容積中での前記粒子の衝突過程をモデル化するステップをさらに含む、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
28. 前記非平衡衝突後分布関数は、流体力学モーメントに対して１次のサポートを提供する格子速度セットに対するマッハ数に関する１次のガリレイ不変衝突演算子
    であり、前記衝突演算子は、
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    ｃ_iは衝突前の前記粒子の速度ベクトルであり、
    ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける前記粒子間の平均速度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    ｗ_iは一定重み係数であり、
    は非平衡運動量流束である、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
29. 前記非平衡衝突後分布関数は、流体力学モーメントに関して無限次数のサポートを提供する格子速度セットに対する衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    は粒子相対速度であり、
    ρは流体密度であり、
    は平衡分布関数であり、
    は非平衡運動量流束である、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
30. 前記非平衡衝突後分布関数は、流体力学モーメントに関して２次のサポートを提供する、格子速度セットに対するマッハ数に関して２次のガリレイ不変衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    ｃ_iは衝突前の前記粒子の速度ベクトルであり、
    ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける前記粒子間の平均速度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    ｗ_iは一定重み係数であり、
    は非平衡運動量流束である、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
31. 前記予め定義された物理量は、前記特定の容積中の前記流体の質量、前記特定の容積中の前記流体の運動量、及び前記特定の容積中の前記流体のエネルギを含む、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
32. 前記非平衡衝突後分布関数は、エネルギ流束に関連する衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    は相対的な粒子速度であり、
    は平衡分布関数であり、
    は非平衡エネルギ流束である、請求項１７に記載の１又は２以上の機械可読ハードウェア記憶装置。
33. １又は２以上の処理装置と、
    オペレーションを実行する１又は２以上の処理装置で実行可能な命令を記憶する１又は２以上の機械可読ハードウェア記憶装置と、
    を備えるシステムであって、前記オペレーションは、
    格子速度セットにおいて流体容積中の粒子の移動をシミュレートするステップであって、前記移動が前記粒子間の衝突を引き起こすものである、ステップと、
    前記シミュレートされた移動に基づいて、
    前記容積内の特定位置における粒子の粒子相対速度を決定するステップであって、前記粒子相対速度が（ｉ）前記容積内の前記特定位置において前記容積の流れがゼロの下で測定される前記粒子の絶対速度と（ｉｉ）前記容積内の前記特定位置における１又は２以上の前記粒子の平均速度との間の差である、ステップと、
    前記粒子相対速度に基づいて、前記衝突を表す特定の次数の非平衡衝突後分布関数を決定するステップと、
    を含む、システム。
34. 前記オペレーションは、粒子速度の或る次数まで流体力学運動をサポートする格子速度を提供するステップをさらに含む、請求項３３に記載のシステム。
35. 前記格子速度セットに関する前記サポートされる次数は、前記非平衡衝突後分布関数の前記特定の次数より小さくかつこれとは異なり、
    前記非平衡衝突後分布関数の前記特定の次数は、前記粒子速度の前記次数によって決定される、請求項３４に記載のシステム。
36. 前記容積内の前記特定位置における前記１又は２以上の前記粒子の前記平均速度は、前記特定位置における特定タイプ粒子の平均速度を含む、請求項３３に記載のシステム。
37. 前記格子速度セットは、格子ボルツマン法に関連する状態ベクトルのセットである、請求項３３に記載のシステム。
38. 前記非平衡衝突後分布関数は、（ｉ）予め定義された物理量に対する非平衡モーメントを維持し、（ｉｉ）未定義の物理量に対する非平衡モーメントを前記特定の次数まで除去する、請求項３３に記載のシステム。
39. 前記特定の次数は、前記流体速度の格子音速に対する比率に関連した指数値であり、前記格子速度セットは、前記指数値をサポートする、請求項３３に記載のシステム。
40. 前記格子速度セットは、格子に限定される空間における運動量状態のセットを含む、請求項３３に記載のシステム。
41. 前記粒子相対速度は、前記容積内の前記特定位置において前記容積の流れがゼロの下で測定される前記粒子の前記絶対速度から前記容積内の前記特定位置における前記１又は２以上の前記粒子の前記平均速度を差し引いたものである、請求項３３に記載のシステム。
42. 前記非平衡衝突後分布関数は、ガリレイ不変フィルタ処理演算子である、請求項３３に記載のシステム。
43. 前記非平衡衝突後分布関数に基づいて、前記流体容積中での前記粒子の衝突過程をモデル化するステップをさらに含む、請求項３３に記載のシステム。
44. 前記非平衡衝突後分布関数は、流体力学モーメントに対して１次のサポートを提供する格子速度セットに対するマッハ数に関する１次のガリレイ不変の衝突演算子
    であり、前記衝突演算子は、
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    ｃ_iは衝突前の前記粒子の速度ベクトルであり、
    ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける前記粒子間の平均速度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    ｗ_iは一定重み係数であり、
    は非平衡運動量流束である、請求項３３に記載のシステム。
45. 前記非平衡衝突後分布関数は、流体力学モーメントに関して無限次数のサポートを提供する格子速度セットに対する衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    は粒子相対速度であり、
    ρは流体密度であり、
    は平衡分布関数であり、
    は非平衡運動量流束である、請求項３３に記載のシステム。
46. 前記非平衡衝突後分布関数は、流体力学モーメントに関して２次のサポートを提供する、格子速度セットに対するマッハ数に関して２次のガリレイ不変衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    ｃ_iは衝突前の前記粒子の速度ベクトルであり、
    ｕ（ｘ，ｔ）は時間ｔでの特定位置ｘにおける前記粒子間の平均速度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    ｗ_iは一定重み係数であり、
    は非平衡運動量流束である、請求項３３に記載のシステム。
47. 前記予め定義された物理量は、前記特定の容積中の前記流体の質量、前記特定の容積中の前記流体の運動量、及び前記特定の容積中の前記流体のエネルギを含む、請求項３３に記載のシステム。
48. 前記非平衡衝突後分布関数は、エネルギ流束に関連する衝突演算子
    であり、前記衝突演算子は
    に従って定義され、
    ｘは前記容積内の前記特定位置であり、
    ｔは時間における特定時点であり、
    ｉは前記格子速度セットでの格子速度の指数であり、
    Ｔ₀は一定格子温度であり、
    Ｉは２階の単位テンソルであり、
    τは衝突緩和時間であり、
    は相対的な粒子速度であり、
    は平衡分布関数であり、
    は非平衡エネルギ流束である、請求項３３に記載のシステム。