JP7402125B2 - 陽的数値拡散問題を安定化させる不規則空間グリッドによる物理流体のコンピュータシミュレーション - Google Patents

Description

本明細書は、物理的流体流れ（physical fluid flow）などの物理的プロセスのコンピュータシミュレーションに関する。

計算流体動力学は、物理的環境における流体流れを分析およびシミュレートするためのコンピュータ実施数値分析技法を含む流体力学の一部門である。

いわゆる「格子ボルツマン法」（ＬＢＭ）は、計算流体動力学で使用される有利な技法である。格子ボルツマン系の基礎をなす動力学は、ボルツマン方程式に従う多くの粒子の運動を含む運動論の基礎物理学に属する。基本的なボルツマンの運動系には２つの基本的な動的プロセス、すなわち、衝突および移流がある。衝突過程は、保存則に従い、平衡まで緩和する粒子間の相互作用を含む。移流過程は、粒子の微視的速度に従ってある場所から別の場所への粒子の移動をモデル化することを含む。

「格子ボルツマン法」（ＬＢＭ）を含む実際の計算流体力学シミュレーション問題に見られる１つの共通の側面は、複雑な形状、例えば、表面および体積の離散化が原因の不規則グリッドに関わる問題である。

拡散支配的な物理現象の数値シミュレーションは、伝導性熱伝達、質量拡散、電気伝導などを含む用途で普通に使用される。これらの現象の支配方程式は、非定常（不安定）拡散項および体積ソース項を含む偏微分方程式（ＰＤＥ）のセットとして定式化される。数値解法は、対象の空間ドメインを離散化することと、時間内に解法を進めるために時間積分技法を適用することとを含む。空間離散化手順は、通常、高度に自動化されたグリッド生成ツールを使用して達成され、一方、時間離散化（時間内、すなわち、時間ステップサイズ内に解法を進めるために時間積分技法を適用するための）は、受け入れ可能な計算コストで数値解の安定および精度を保証するように選ばれる。

特に、時間前進スキームのＣｏｕｒａｎｔ－Ｆｒｉｅｄｒｉｃｈｓ－Ｌｅｗｙ（ＣＦＬ）制約と一般に呼ばれる安定特性は、解に大幅な不安定性を導入することなく使用することができる最大時間ステップサイズを決定する。２つのタイプの時間進行スキーム、すなわち、陰的（implicit）および陽的（explicit）スキームが、一般に、使用される。

陰的方法は、構築によってＣＦＬ制約を満たし、したがって、大きい時間ステップが、解を不安定にすることなく使用され得る（しかしながら、あまりにも大きい時間ステップは、通常、不正確な結果をもたらす）。陰的方法は、行列係数の大規模システムの解を必要とし、したがって、それらの実際の実施を非自明にし、計算上高価にする。

他方、陽的方法は、実施することが非常に簡単であり、計算上安価（反復当たり）であり、高度に並列化可能である。しかしながら、陽的方法は、ＣＦＬ制約を満たす必要がある。陽的拡散スキームのこのＣＦＬ制約は、（κΔ＿ｔ）／（Δ＿ｘ＾２）によって与えられるＣＦＬ数が特定の限界（Ｏ（１）である）未満であることを要求し、ここで、κは、拡散率であり、Δ＿ｘは、最小空間グリッドのサイズであり、Δ＿ｔは、時間ステップである。言い換えれば、空間グリッドサイズΔ＿ｘが、係数Ｆだけドメイン内のどこかで減少する場合、時間ステップΔ＿ｔは、数値の安定を維持するために係数Ｆ²だけ減少しなければならないことになる。

したがって、陽的方法は、シミュレーション性能に著しく影響を与える小さいサイズの要素を有する空間グリッドでは極端に小さい時間ステップを必要とすることがある。これは、そのような小さいサイズの要素の数がシミュレーションドメイン内で非常に限定されている場合でさえ当てはまり、その理由は、ドメイン全体内の最小要素がＣＦＬ条件を決定し、その結果、時間ステップサイズを決定するからである。

複雑な形状を含む実際の問題では、不規則なグリッドの使用が、表面および体積の離散化では避けられない。これらのグリッド上で、Δ＿ｘは、著しく変化することがあり、陽解法の使用は、ＣＦＬ制約によって必要とされる極端に小さい時間ステップに起因して非常に非効率的になることがある。

それゆえに、陽解法の実践者は、空間グリッドの品質を改善しようとして大量の時間および労力を費やし、問題の緩和を試みる。その場合にも、現実的な空間ドメインの離散化からすべての小さいサイズの要素を除去することはほとんど不可能であり、その結果、小さい時間ステップ（少なくとも局所的に）が陽的解を安定化させる唯一の方法である。

一態様によれば、固体の表面のまわりの流体流れをシミュレートするためのシステム、コンピュータ実施方法、およびコンピュータプログラム製品は、ボクセルの集合として表された格子構造と、物理的オブジェクトの表現とを含むシミュレーション空間のモデルを受け取ることであり、ボクセルが、物理的オブジェクトの表面を説明するための適切な分解能を有する、受け取ることと、流体の体積中の粒子の移動をシミュレートすることであり、粒子の移動が粒子間の衝突を引き起こす、シミュレートすることと、２つのボクセル間の面（ファセット）を識別することであり、面のうちの少なくとも１つが安定条件に違反している、識別することと、２つのボクセルの近傍の空間的に平均化された温度勾配を使用して修正熱流束を計算することであり、面の少なくとも１つが安定条件に違反している、計算することと、後続のボクセルへの粒子の移流操作を実行することとを含む。

開示された技法は、要素のうちの少なくとも１つが、制約、例えば、Ｃｏｕｒａｎｔ－Ｆｒｉｅｄｒｉｃｈｓ－Ｌｅｗｙ（ＣＦＬ）制約）に違反している場合、２つの隣接する要素間の熱流束計算に修正を導入する。熱流束計算へのこれらの修正は、２つの要素の物質および幾何学的特性、ならびに要素のすぐ近くの対象の量の既存の状態に依存し、２つの要素のサイズに関係なく数値解を安定させるのを支援し、空間－時間精度を保証する。２つの隣接する要素が大きい（それゆえに、ＣＦＬ制約を満たす）場合、新しく提案する流束計算は、陽解法実施態様になり、それは、この新規の手法が陽的手法と一致しており、さらに、陽的手法で上記の不備を克服していることを意味する。

態様は、方法、コンピュータプログラム製品、１つまたは複数の機械可読ハードウェアストレージデバイス、装置、およびコンピューティングシステムを含む。

他の特徴および利点は、図面を含む以下の説明および特許請求の範囲から明らかになる。

流体流れのシミュレーションのためのシステムを示す図である。 熱流束計算を用いた格子ボルツマンモデルに基づく操作シミュレーションを示す流れ図である。 要素の少なくとも１つが制約に違反している場合の２つの隣接する要素間の熱流束計算への修正の態様を含む流れ図である。 要素の少なくとも１つが制約に違反している場合の２つの隣接する要素間の熱流束計算への修正の態様を含む流れ図である。 要素の少なくとも１つが制約に違反している場合の２つの隣接する要素間の熱流束計算への修正の態様を含む流れ図である。 要素の少なくとも１つが制約に違反している場合の２つの隣接する要素間の熱流束計算への修正の態様を含む流れ図である。 要素の少なくとも１つが制約に違反している場合の２つの隣接する要素間の熱流束計算への修正の態様を含む流れ図である。 ユークリッド空間で表わされた２つのＬＢＭモデルの速度成分を示す図である（先行技術）。 ユークリッド空間で表わされた２つのＬＢＭモデルの速度成分を示す図である（先行技術）。 ２つの隣接の要素間の修正熱流束計算を使用する物理的プロセスシミュレーションシステムが従う手順の流れ図である。 マイクロブロックの斜視図である（先行技術）。 図１のシステムによって使用される格子構造の図である（先行技術）。 図１のシステムによって使用される格子構造の図である（先行技術）。 可変分解能技法を示す図である（先行技術）。 可変分解能技法を示す図である（先行技術）。 粒子の移動を示す図である（先行技術）。 表面のファセットによって影響される領域を示す図である（先行技術）。 表面動力学を示す図である（先行技術）。 表面動力学を実行するための手順の流れ図である。

モデルシミュレーション空間
ＬＢＭベース物理的プロセスシミュレーションシステムでは、流体流れは離散速度ｃ_iのセットで評価される分布関数値ｆ_iによって表される。分布関数の動力学は式Ｉ．１によって支配される。
ｆ_i（ｘ＋ｃ_i，ｔ＋１）＝ｆ_i（ｘ，ｔ）＋Ｃ_i（ｘ，ｔ） 式（Ｉ．１）
この式は、分布関数ｆ_iの時間発展を記述するよく知られた格子ボルツマン方程式である。左辺は、いわゆる「流動過程」による分布の変化を表す。流動過程は、流体のポケットが、あるメッシュ位置から出発して、次いで、複数の速度ベクトルのうちの１つに沿って次のメッシュ位置まで移動する場合である。その時点で、「衝突係数」、すなわち、流体の近傍ポケットの流体の開始ポケットへの影響が計算される。流体は別のメッシュ位置にしか移動できず、そのため、すべての速度のすべての成分が共通速度の倍数であるように、速度ベクトルの適切な選択が必要である。

第１の式の右側は、流体のポケット間の衝突による分布関数の変化を表す前記の「衝突演算子」である。衝突演算子の特定の形態は、限定はしないが、Ｂｈａｔｎａｇａｒ、Ｇｒｏｓｓ、およびＫｒｏｏｋ（ＢＧＫ）演算子の形態とすることができる。衝突演算子は、分布関数を、「平衡」形態である第２の式によって与えられる規定値に近づける。

ＢＧＫ演算子は、衝突の詳細に関係なく、分布関数が、衝突

を介して、｛ｆ^eq（ｘ，ν，ｔ）｝で与えられる明確な局所平衡に近づくという物理的論拠によって構築され、ここで、パラメータτは、衝突を介した平衡までの固有緩和時間を表す。

このシミュレーションから、質量ρおよび流速ｕなどの従来の流体変数は、式（Ｉ．３）における単純な和として得られる。以下を参照されたい。

対称性の考慮により、速度値のセットは、配置空間にわたるときに特定の格子構造を形成するように選択される。そのような離散系の動力学は、形態
ｆ_i（ｘ＋ｃ_i，ｔ＋１）＝ｆ_i（ｘ，ｔ）＋Ｃ_i（ｘ，ｔ）
を有するＬＢＭ方程式に従い、ここで、衝突演算子は、通常、上述のようなＢＧＫ形態をとる。平衡分布形態の適切な選択によって、格子ボルツマン方程式が、正しい流体力学および熱流体力学結果をもたらすことを理論的に示すことができる。すなわち、ｆ_i（ｘ，ｔ）から導出された流体力学的モーメントは、巨視的極限ではナビエ－ストークス方程式に従う。これらのモーメントは、
ρ（ｘ，ｔ）＝Σ_iｆ_i（ｘ，ｔ）；ρ（ｘ，ｔ）ｕ（ｘ，ｔ）＝Σ_iｃ_iｆ_i（ｘ，ｔ） 式（Ｉ．３）
として定義され、ここで、ρおよびｕは、それぞれ、流体密度および速度である。

ｃ_iおよびω_iの集合値はＬＢＭモデルを定義する。ＬＢＭモデルは、拡張可能なコンピュータプラットフォームで効率的に実施することができ、時間的に非定常な流れおよび複雑な境界条件に対して優れた頑健性で実行することができる。

流体系の巨視的な運動方程式をボルツマン方程式から得る標準的な技法は、完全なボルツマン方程式の逐次近似が取り入れられるＣｈａｐｍａｎ－Ｅｎｓｋｏｇ法である。流体系では、密度の小さい擾乱は音速で移動する。ガス系では、音速は、一般に、温度によって決定される。流れにおける圧縮性の影響の重要性は、マッハ数として知られる、特性速度と音速との比によって測定される。

従来のＬＢＭベース物理的プロセスシミュレーションシステムのさらなる説明に関しては、解釈が、米国特許出願公開第２０１６－０１８８７６８号に言及されており、その全内容が、参照により本明細書に組み込まれる。

図１を参照すると、例えば物理的オブジェクトの表現に関して、流体流れをシミュレートするためのシステム１０が示される。この実施態様のシステム１０は、クライアント－サーバアーキテクチャに基づき、大規模並列コンピューティングシステム１２として実装されたサーバシステム１２と、クライアントシステム１４とを含む。サーバシステム１２は、メモリ１８と、バスシステム２２と、インタフェース２０（例えば、ユーザインタフェース／ネットワークインタフェース／ディスプレイまたはモニタインタフェースなど）と、メッシュおよびシミュレーションエンジン３４と一緒にシミュレーションプロセス３０を行う処理デバイス２４とを含む。

メモリ１８には、メッシュ準備エンジン３２およびシミュレーションエンジン３４がある。図１はメッシュ準備エンジン３２をメモリ１８内に示しているが、メッシュ準備エンジンは、サーバ１２とは異なるシステム（例えば、システム１４または別のシステム）で実行されるサードパーティアプリケーションとすることができる。メッシュ準備エンジン３２がメモリ１８で実行されるかまたはサーバ１２とは異なるシステムで実行されるかにかかわらず、メッシュ準備エンジン３２はユーザ指定のメッシュ定義２８を受け取り、メッシュ準備エンジン３２は、メッシュを準備し、準備したメッシュをシミュレーションエンジン３４に送る。

シミュレーションエンジン３４は、粒子衝突相互作用モジュール３４ａと、粒子境界モデルモジュール３４ｂと、移流操作を実行する移流モジュール３４ｃとを含む。システム１０は、２Ｄおよび／または３Ｄメッシュおよびライブラリを格納するデータレポジトリ３８にアクセスする。移流モジュール３４ｃは、以下で論じるように、２つの隣接する要素間の流束計算の修正に従って移流操作を実行するサブモジュール３６を含む。

シミュレーションエンジンでシミュレーションを実行する前に、シミュレーション空間が、ボクセルの集合としてモデル化される。一般に、シミュレーション空間は、コンピュータ支援設計（ＣＡＤ）プログラムを使用して生成される。例えば、ＣＡＤプログラムを使用して、風洞に位置づけられたマイクロデバイスを描くことができる。その後、ＣＡＤプログラムによって作り出されたデータを処理して、適切な分解能を有する格子構造を追加し、シミュレーション空間内のオブジェクトおよび表面を明らかにする。

次に、図２を参照すると、物理的オブジェクトの表現に関する流体流れをシミュレートするプロセスが示される。本明細書で論じる例では、物理的オブジェクトは翼である。しかしながら、物理的オブジェクトは任意の形状とすることができ、特に、平面および／または曲面を有することができるので、翼の使用は単に例示である。プロセスは、例えば、クライアントシステム１４から、またはデータレポジトリからの検索によって、シミュレートされる物理的オブジェクトのためのメッシュを受け取る３５ａ。他の実施形態では、外部システムまたはサーバ１２のいずれかが、ユーザ入力に基づいて、シミュレートされる物理的オブジェクトのメッシュを生成する。プロセスは、検索したメッシュから幾何学量を事前計算し３５ｂ、検索したメッシュに対応する事前計算された幾何学量を使用して、動的格子ボルツマンモデルシミュレーションを実行する３５ｃ。格子ボルツマンモデルシミュレーションは、粒子分布３４ａの発展３５ｄのシミュレーションと、エンジン３４ｂ（図１）から作り出された境界決定（図２には図示せず）による、ＬＢＭメッシュにおける次のセルｑ（上にバー）への粒子の移流３５ｅとを含む。移流３５ｃプロセスは、ＣＦＬ制約違反３７ａをテストし、見つかった場合、エンジン３６（図１）を使用して流束計算３７ｂに修正を適用する。

不規則な空間グリッドに関する拡散問題の陽の数値スキームの安定化
この説明のために、陽的オイラースキームおよび有限体積定式化が仮定される。以下の説明において、対象の量は温度であり、支配方程式は熱伝導方程式である。数値スキームは、要素のすべての面の熱流束を計算する必要がある。続いて、これらの流束は、合計されて、検討中の要素の温度を更新するために使用される。２つの面共有隣接要素αおよびβを考える。熱伝導のフーリエの法則によれば、熱流束は、

であり、ここで、

は、共通面の熱伝導率であり、

は、共通面に垂直な温度勾配であり、「ｍ」は、量が時間ステップ「ｍ」で評価されることを明示するために使用される。熱流入α（熱流出αの代わりに）が考慮されるので、フーリエの法則の一般に使用される形式の負号は省かれる。使用される温度勾配は、特に、異なるサイズの要素が存在する状態で、粒子移流の滑らかさを確実するように計算される。２つの隣接する要素αおよびβがＣＦＬの制約を満たす場合、時間ステップｍからｍ＋１への間に共通面を横切るエネルギー移送は、熱流束に共通面の面積Ａ^αβおよび時間ステップサイズΔ_tを乗算することによって得られ、すなわち、

である。

従来の手法では、時間ステップの終わりにおける要素αの最終温度は、αへの正味エネルギー移送（すべての面からのエネルギー移送の合計）から計算される。

上述の式（３）は、温度変化が正味熱流束に比例し、要素のサイズ∀^αに反比例する、すなわち、小さい要素では、同じ正味エネルギー移送がより大きい温度変化をもたらすことを述べていることに留意されたい。

図３を参照すると、流束修正計算プロセス３６は、時間ステップの決定を含む。流束計算プロセス３６において、時間ステップがテストされる３６ａ。時間ステップΔ_tが、２つの要素の少なくとも一方でＣＦＬ制約に違反する３６ｂほど大きくない場合、上述の形態（式３）は、解が数値的不安定になる可能性が低いことになり、したがって、上述の形態の流束計算を使用することができる。時間ステップΔ_tが、２つの要素の少なくとも一方でＣＦＬ制約に違反する３６ｃほど大きい場合、上述の形態（式３）は、解が数値的不安定になることがある。この場合、修正流束計算手法３６ｄが適用される。

一般性を失うことなしに、修正手法は以下のように説明することができる。

要素αは隣接する要素βよりも小さく、それゆえに、ＣＦＬ制約が少なくとも要素αで違反されると仮定する。この数値不安定は、要素α（そのサイズが小さいと仮定されている）では、

を計算するために使用される温度勾配が、時間ステップΔ_tの期間の全体を通して一定値のままであると仮定することが正しくないという理由で生じる。上記のように、同じ正味エネルギー移送では、小さい要素の温度変化は大きく、したがって、標準の陽的時間積分は、一定温度勾配の仮定が有効であることを保証するために時間ステップを減少させることを必要とする。所与のΔ_tでは、この問題は、問題に非定常性が存在する限り存続し、すべての要素に関して入って来るおよび出て行く流束がすべて互いに正確に釣り合うときの定常状態でのみ解決する。

図４を参照すると、修正手法の一部として、流束修正された流束計算は、式（１）におけるように定義された項

を、以下の２つの部分に細分する４９ａ。２つの部分は、（１ａ）αの温度発展（上述の合計において）に向けて使用されることになる付加流束（ａｐｐｌｉｅｄ ｆｌｕｘ）４９ｂ

と、（１ｂ）要素αの温度変化なしに、インタフェースαβの反対側に送出されることになる平衡流束４９ｃ

とであり、それらは、以下のように、

として表され、ここで、項

およびΔＧは、

によって与えられ、そして、

として表される。

上述の式において、幾何学的フィーチャは、温度勾配の計算で使用される距離ｄ^αβと、要素体積∀^αおよび∀^βと、共通面（隣接する要素によって共有される）の面積Ａ^αβとによって表される。材料特性は、

および

（ここで、ρは密度を示し、Ｃ_pは比熱を示す）、ならびに熱伝導率

によって説明される。加算における流束項

および

は、それぞれ、要素αおよびβの異なる面に存在する可能性がある流束の推定量を提供する。

２つの流束項

および

の物理的解釈は、以下の通りである。付加流束

は、数値不安定を導入することなく要素αの温度発展に向けて使用することができる全流束

（式（１）により与えられる）の一部分を表す。この形態は、要素αおよびβからなる孤立系のための第１原理から導出され、この系の近くで連続的に発展する温度場の影響の推定量を含むことができる。この理由で、この項は、要素の幾何学的／熱的特性

ならびに環境とのこれらの要素の相互作用（ΔＧ）の両方への依存を示す。

式（４）において、

は、依然として、以前の時間ステップに他の面で観察された流束に依存しており、それによって、現在の時間ステップの間にそれらの面で継続している相互作用の推定量を提供することに留意すべきである。強硬な過渡的問題では、これは、

と

との間にずれをもたらす。

平衡流束と呼ばれる第２の項

は、このずれを説明しており、それはαを通ってインタフェースαβの反対側の隣接する要素に伝えられる。この平衡流束は、ＣＦＬ制約を満たすのに十分な大きさの要素に置かれるときのみ温度発展に付加され、それまで、流束方向に沿って連続的に送出される。

上述のスキームは、２つの要素αとβとの間のインタフェースにおいて、小さい要素αの温度発展に有効な流束の量を正確に制御しながら、全流束の正しい量

が組み入れられることを厳格に保証する。全体として、このスキームは、数値安定、ならびに良好な空間および時間精度を維持することができる。最後に、定常状態では、付加流束は全流束に等しくなり、すなわち、

であり、結果として、平衡流束は０に等しくなり、すなわち、

であることに留意すべきである。

修正流束計算プロセスの特徴
上述の説明で言及したように、そして図５～図７に示すように、修正流束計算手法３６は、いくつかのアルゴリズムのプロセスを必要とする。修正流束計算プロセスは、熱流束の修正された定義が適用されるべきである面を識別する５２ａことを含む。ＣＦＬ条件に違反していないすべての面では、熱流束の標準定義が使用される５２ｂ（式３）。ＣＦＬ条件に違反しているすべての面において、例えば、少なくとも１つがＣＦＬ条件基準に違反している２つの要素間のいかなる面においても、修正熱流束計算が使用される５２ｃ。全熱流束量は、解の滑らかさを保証するために、検討中の要素の近傍の空間的に平均化された温度勾配を使用して計算される。（対照的に、標準熱流束は、従来の差分形式に基づいて計算された温度勾配を利用する。）修正熱流束計算プロセス５２は、修正熱流束を２つの部分、すなわち、付加流束４９ｂの項および平衡流束４９ｃの項に分割する５２ｄことを含む。次いで、流束計算は、分割に従って実行される５２ｅ。

次に、図６を参照すると、流束の計算は、検討中の要素の温度発展式で常に使用される付加流束計算５２ｆ（式３）を使用する。

次に、図７を参照すると、平衡流束は、含まれる要素のサイズに応じて温度発展に使用されることもされないこともある。要素がＣＦＬ制約に違反するほど小さい５２ｇ場合、平衡流束が、計算され５２ｈ、流束の方向に隣接する要素に送出され５２ｉ、隣接する要素がＣＦＬ制約に違反しないほど大きい場合、平衡流束が、温度発展に向けて使用される５２ｊ。平衡流束は、平衡化流束が、最終的に、平衡化流束を温度発展に向けて適用する十分大きい要素（ＣＦＬ制約に違反しないほどの）に移送されるまで、流束の方向に沿って連続的に送出される。

アルゴリズムの独自性および利点
要素サイズが変化する不規則なグリッドの拡散問題における数値不安定の問題を克服するためのいくつかの手法が知られている。最も一般に使用される手法は、そのようなシナリオを減じるために、グリッド生成ツールに追加の制約を強制することである。その場合にも、この問題を完全には避けることができないので、安定を保証するほど十分に小さいグローバル時間ステップを使用するか、または小さいグリッド要素に出会ったとき局所サブサイクリングを使用することが一般的な方法である。第１の手法（小さいグローバル時間ステップ）は、空間グリッドのどこかに小さいサイズの要素が単一で生じた場合でさえ計算コストを事実上増加させ、一方、第２の手法（局所サブサイクリング）は、アルゴリズムおよびその実施の複雑さを増加させる。代替の手法は、陽解法の代わりに陰解法を使用するか、または少なくとも陽解法の使用を小さい要素の近くの局所領域に制限することである。この陰解法手法は、実施の複雑さ、ならびに大規模並列処理コンピュータ実施に都合のよくない連立方程式をもたらす解の非局所的性質を欠点として有する。

対照的に、修正流束計算手法は、いくつかの明確な利点を提供する。修正流束計算手法は、グリッドの最小要素のサイズではなく、時間精度の考察に基づいて選択された単一時間ステップサイズの使用を可能にする。これは、考えられるすべてのシナリオに対して、計算コストおよび実施の容易さの点から大きな利点である。修正流束計算手法は、２つの隣接する要素の幾何学的特性への依存性を有し、それにより、修正流束計算手法は要素のサイズに関係なく機能することを保証する。したがって、グリッド生成プロセスへの通常の制約（グリッド品質、サイズなど）を大幅に緩和することができる。項の数学的形式が簡単であり、反復を含まないので、様々な項を計算する計算コストは妥当である。これは、高い計算コストを有することがある既存の手法（時間ステップサイズを縮小するか、またはハイブリッド陰解法－陽解法を使用する）とは際立って対照的である。加えて、定式化の体積特性のために、このスキームは、多くの適用での要件である正確な保存性を維持する。修正流束計算は、依然として、本質的に陽的であり、検討中の要素から小さい距離内の要素からの情報を必要とする。陽とは、オリジナルのシステム実施態様による既存の計算システムへの最低限の変更が必要とされることを意味する。オリジナルの陽的方法の並列化特性は保持され、したがって、修正流束計算手法は、大規模並列処理コンピュータ実施態様を利用することができる。

したがって、修正流束計算手法は、既存の方法と比較して利点がある。流束計算修正手法の上述の説明は熱伝導の有限体積定式化に基づいていたが、この手法は、実際は、多くの拡散支配問題に適用可能である。

モデルシミュレーション空間
ＬＢＭベース物理的プロセスシミュレーションシステムにおいて、流体流れは、離散速度ｃ_iのセットで評価された分布関数値ｆ_iによって表される。分布関数の動力学は、式Ｉ１によって支配され、ここで、ｆ_i（０）は、平衡分布関数として知られており、

として定義され、ここで、

である。

この式は、分布関数ｆ_iの時間発展を記載するよく知られた格子ボルツマン方程式である。左辺は、いわゆる「流動過程」による分布の変化を表す。流動過程は、流体のポケットが、あるメッシュ位置から出発して、次いで、複数の速度ベクトルのうちの１つに沿って次のメッシュ位置まで移動する。その時点で、「衝突係数」、すなわち、流体の近傍ポケットの流体の開始ポケットへの影響が、計算される。流体は、別のメッシュ位置にしか移動することができず、そのため、すべての速度のすべての成分が共通速度の倍数となるように、速度ベクトルの適切な選択が必要である。

第１の式の右側は、流体のポケット間の衝突による分布関数の変化を表す前記の「衝突演算子」である。衝突演算子の特定の形態は、Ｂｈａｔｎａｇａｒ、Ｇｒｏｓｓ、およびＫｒｏｏｋ（ＢＧＫ）演算子の形態である。衝突演算子は、分布関数を、「平衡」形態である第２の式によって与えられる規定値に近づける。

ＢＧＫ演算子は、衝突の詳細に関係なく、分布関数が、衝突

を介して、｛ｆ^eq（ｘ，ν，ｔ）｝によって与えられる明確な局所平衡に近づくという物理的論拠によって構築され、ここで、パラメータτは、衝突を介した平衡までの固有緩和時間を表す。粒子（例えば、原子または分子）を扱う場合、緩和時間は、一般に、定数と見なされる。

このシミュレーションから、質量ρおよび流速ｕなどの従来の流体変数は、式（Ｉ３）における簡単な加算として得られる。

ここで、ρ、ｕ、およびＴは、それぞれ、流体の密度、速度、および温度であり、Ｄは、離散化された速度空間次元である（物理的空間次元と必ずしも等しくない）。

対称性の考慮により、速度値のセットは、配置空間にわたるときに特定の格子構造を形成するように選択される。そのような離散系の動力学は、形態
ｆ_i（ｘ＋ｃ_i，ｔ＋１）－ｆ_i（ｘ，ｔ）＝Ｃ_i（ｘ，ｔ）
を有するＬＢＭ式に従い、ここで、衝突演算子は、通常、上述のようなＢＧＫ形態をとる。平衡分布形態の適切な選択によって、格子ボルツマン方程式が、正しい流体力学および熱流体力学をもたらすことを理論的に示すことができる。すなわち、ｆ_i（ｘ，ｔ）から導出された流体力学的モーメントは、巨視的極限ではナビエ－ストークス方程式に従う。これらのモーメントは、上述の式（Ｉ３）によって定義される。

ｃ_iおよびｗ_iの集合値はＬＢＭモデルを定義する。ＬＢＭモデルは、拡張可能なコンピュータプラットフォームで効率的に実施することができ、時間的に非定常な流れおよび複雑な境界条件に対して優れた頑健性で実行することができる。

流体系の巨視的な運動方程式をボルツマン方程式から得る標準的な技法は、完全なボルツマン方程式の逐次近似が取り入れられるＣｈａｐｍａｎ－Ｅｎｓｋｏｇ法である。流体系では、密度の小さい擾乱は音速で移動する。ガス系では、音速は、一般に、温度によって決定される。流れにおける圧縮性の影響の重要性は、マッハ数として知られる、特性速度と音速との比によって測定される。

流体流れのシミュレーションを行うためにＣＡＤプロセスとともに使用することができるＬＢＭベースシミュレーションシステムの一般的な議論が、以下で提供される。

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

図９を参照すると、第２のモデル（３Ｄ－１）２６０、すなわち、３９個の速度を含む３次元モデルが示されており、各速度は、図９の矢印のうちの１つによって表される。これらの３９個の速度のうち、１つは、移動していない粒子を表し、６つの速度の３組は、格子のｘ軸、ｙ軸、またはｚ軸に沿って正方向または負方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、または正規化速度の３倍（３ｒ）のいずれかで移動している粒子を表し、８つは、ｘ、ｙ、ｚ格子軸の３つのすべてに対して正規化速度（ｒ）で移動している粒子を表し、１２個は、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度の２倍（２ｒ）で移動している粒子を表す。

１０１個の速度を含む３Ｄ－２モデル、および３７個の速度を含む２Ｄ－２モデルなどのさらに複雑なモデルを使用することもできる。３次元モデル３Ｄ－２では、１０１個の速度のうち、１つは、移動していない粒子を表し（グループ１）、６つの速度の３組は、格子のｘ軸、ｙ軸、またはｚ軸に沿って正方向または負方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、または正規化速度の３倍（３ｒ）のいずれかで移動している粒子を表し（グループ２，４、および７）、８つの速度の３組は、ｘ、ｙ、ｚ格子軸の３つのすべてに対して正規化速度（ｒ）、正規化速度の２倍（２ｒ）、または正規化速度の３倍（３ｒ）で移動している粒子を表し（グループ３、８、および１０）、１２個は、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度の２倍（２ｒ）で移動している粒子を表し（グループ６）、２４個は、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度（ｒ）および正規化速度の２倍（２ｒ）で移動しており、残りの軸に対して移動していない粒子を表し（グループ５）、２４個は、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度（ｒ）で、および残りの軸に対して正規化速度の３倍（３ｒ）で移動している粒子を表す（グループ９）。

２次元モデル２Ｄ－２では、３７個の速度のうち、１つは、移動していない粒子を表し（グループ１）、４つの速度の３組は、格子のｘ軸またはｙ軸のいずれかに沿って正方向または負方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、または正規化速度の３倍（３ｒ）のいずれかで移動している粒子を表し（グループ２，４、および７）、４つの速度の２組は、ｘ格子軸とｙ格子軸の両方に対して正規化速度（ｒ）または正規化速度の２倍（２ｒ）で移動している粒子を表し、８つの速度は、ｘ格子軸およびｙ格子軸の一方に対して正規化速度（ｒ）で、および他の軸に対して正規化速度の２倍（２ｒ）で移動している粒子を表し、８つの速度は、ｘ格子軸およびｙ格子軸の一方に対して正規化速度（ｒ）で、および他の軸に対して正規化速度の３倍（３ｒ）で移動している粒子を表す。

上述のＬＢＭモデルは、２次元と３次元の両方における流れの数値シミュレーションのための効率的で頑健な離散速度動力学モデルの特定のクラスを提供する。この種類のモデルは、離散速度とその速度に関連する重みとの特定のセットを含む。速度は、速度空間におけるデカルト座標のグリッド点と一致し、それは、離散速度モデル、特に、格子ボルツマンモデルとして知られている種類の正確で効率的な実施を容易にする。そのようなモデルを使用して、高い忠実度で流れをシミュレートすることができる。

図１０を参照すると、物理的プロセスシミュレーションシステムは、上述で論じたようにＣＡＤプロセスを使用して流体流れなどの物理的プロセスをシミュレートするための手順２７０に従って動作する。シミュレーションの前に、シミュレーション空間がボクセルの集合としてモデル化される（ステップ２７２）。一般に、シミュレーション空間は、コンピュータ支援設計（ＣＡＤ）プログラムを使用して生成される。例えば、ＣＡＤプログラムを使用して、風洞に位置づけられたマイクロデバイスを描くことができる。その後、ＣＡＤプログラムによって作り出されたデータを処理して、適切な分解能を有する格子構造を追加し、シミュレーション空間内のオブジェクトおよび表面を明らかにする。

物理的プロセスシミュレーションシステムは、上述で論じた修正流束計算プロセスを使用して、手順２７０に従って動作する。格子の分解能は、シミュレートされているシステムのレイノルズ数に基づいて選択することができる。レイノルズ数は、流れの粘性（ν）、流れ中のオブジェクトの特性長（Ｌ）、および流れの特性速度（ｕ）に関連し、
Ｒｅ＝ｕＬ／ν 式（Ｉ４）
である。

オブジェクトの特性長は、オブジェクトの大きいスケールのフィーチャを表す。例えば、マイクロデバイスのまわりの流れがシミュレートされる場合、マイクロデバイスの高さを特性長であると見なすことができる。オブジェクト（例えば、自動車のサイドミラー）の小さい領域のまわりの流れが対象である場合、シミュレーションの分解能を向上させることができ、または分解能の向上した区域を対象の領域のまわりで使用することができる。ボクセルの寸法は、格子の分解能の向上につれて減少する。

状態空間は、ｆ_i（ｘ，ｔ）として表され、ここで、ｆ_iは、時間ｔでの３次元ベクトルｘによって表わされる格子サイトにおける状態ｉの単位体積当たりの要素または粒子の数（すなわち、状態ｉの粒子の密度）を表す。既知の時間増分に対して、粒子の数は単にｆ_i（ｘ）と呼ばれる。格子サイトのすべての状態の組合せは、ｆ（ｘ）として表わされる。

状態の数は、各エネルギーレベル内の可能な速度ベクトルの数によって決定される。速度ベクトルは、３つの次元、ｘ、ｙ、およびｚを有する空間における整数の直線速度からなる。状態の数は、多数の種のシミュレーションでは増加する。

各状態ｉは、特定のエネルギーレベル（すなわち、エネルギーレベル０、１、または２）における異なる速度ベクトルを表す。各状態の速度ｃ_iは、以下のように３次元の各々の「速度」で示される。
ｃ_i＝（ｃ_ix，ｃ_iν，ｃ_iz） 式（Ｉ５）

エネルギーレベル０の状態は、いずれの次元にも移動していない停止した粒子を表し、すなわちｃ_stopped＝（０，０，０）である。エネルギーレベル１の状態は、３つの次元のうちの１つで±１の速度および他の２つの次元で０の速度を有する粒子を表す。エネルギーレベル２の状態は、３つの次元すべてで±１の速度、または３つの次元のうちの１つで±２の速度および他の２つの次元で０の速度を有する粒子を表す。

３つのエネルギーレベルの可能な並べ換えのすべてを生成すると、合計で３９個の可能な状態が与えられる（１つのエネルギー０の状態、６つのエネルギー１の状態、８つのエネルギー３の状態、６つのエネルギー４の状態、１２個のエネルギー８の状態、および６つのエネルギー９の状態）。

各ボクセル（すなわち、各格子サイト）は、状態ベクトルｆ（ｘ）によって表される。状態ベクトルは、ボクセルの状態を完全に定義し、３９個のエントリを含む。３９個のエントリは、１つのエネルギー０の状態、６つのエネルギー１の状態、８つのエネルギー３の状態、６つのエネルギー４の状態、１２個のエネルギー８の状態、および６つのエネルギー９の状態に対応する。この速度セットを使用することによって、システムは、達成される平衡状態ベクトルに対してマクスウェル－ボルツマン統計を作り出すことができる。

処理効率のために、ボクセルは、マイクロブロックと呼ばれる２×２×２体積にグループ化される。マイクロブロックは、ボクセルの並列処理を可能にし、データ構造に関連するオーバーヘッドを最小にするように編成される。マイクロブロックのボクセルの簡易表記は、Ｎ_i（ｎ）として定義され、ここで、ｎは、マイクロブロック内の格子サイトの相対位置を表し、ｎ∈｛０，１，２，・・・，７｝である。

マイクロブロックが図１１に示される。

図１２Ａおよび図１２Ｂを参照すると、表面Ｓ（図１２Ａ）は、シミュレーション空間（図１２Ｂ）にファセットＦ_αの集合として表されており、
Ｓ＝｛Ｆ_α｝ 式（Ｉ６）
であり、ここで、αは、特定のファセットを挙げるインデクスである。ファセットは、ボクセル境界に制限されるのではなくて、ファセットが比較的少数のボクセルに影響を与えるように、一般に、ファセットに隣接するボクセルのサイズと同等かまたはそれよりもわずかに小さいサイズにされる。表面動力学を実施するために、特性がファセットに割り当てられる。特に、各ファセットＦ_αは、単位法線（ｎ_α）、表面積（Ａ_α）、中心位置（ｘ_α）、およびファセットの表面動力学特性を記述するファセット分布関数（ｆ_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αが１つまたは複数のファセットと交差するとき、以下条件は真である。
Ｖ_iα＝ΣＶ_α（ｘ）＋ΣＶ_iα（β） 式（Ｉ９）
ここで、第１の合計は、Ｇ_iαが重なったすべてのボクセルに相当し、第２の項は、Ｇ_iαと交差するすべてのファセットに相当する。平行六面体Ｇ_iαが別のファセットと交差しない場合、この式は、
Ｖ_iα＝ΣＶ_iα（ｘ） 式（Ｉ１０）
である。

Ｄ．シミュレーションの実行
１つまたは複数のファセットによって影響されるボクセルが識別された（ステップ２７４）後、シミュレーションを始めるために、タイマーが初期化される（ステップ２７６）。シミュレーションの各時間増分の間に、ボクセルからボクセルへの粒子の移動は、粒子と表面ファセットとの相互作用を引き起こす移流段階（ステップ２７８～２８６）によってシミュレートされる。次に、衝突段階（ステップ２８８）は、各ボクセル内の粒子の相互作用をシミュレートする。その後、タイマーが増分される（ステップ２９０）。増分されたタイマーが、シミュレーションの完了（ステップ２９２）を示さない場合、移流および衝突段階（ステップ２７８～２９０）が繰り返される。増分されたタイマーが、シミュレーションの完了（ステップ２９２）を示す場合、シミュレーションの結果が、格納および／または表示される（ステップ２９４）。

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

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

２．ボクセルからファセットへの収集
粒子と表面との間の相互作用をシミュレートする際の第１のステップとして、粒子が、ボクセルから収集され、ファセットに供給される（ステップ２７８）。上記のように、ボクセルＮ（ｘ）とファセットＦ_αとの間の状態ｉの粒子の流束は、
Γ_iα（ｘ）＝Ｎ_i（ｘ）Ｖ_iα（ｘ） 式（Ｉ１１）
である。

この式から、ファセットＦ_αに向けられた状態ｉ（ｃ_iｎ_α＜０）ごとに、ボクセルがファセットＦ_αに供給する粒子の数は、
Γ_iαV→F＝Σ_XΓ_iα（ｘ）＝Σ_XＮ_i（ｘ）Ｖ_iα（ｘ） 式（Ｉ１２）
である。

Ｖ_iα（ｘ）が０でない値を有するボクセルのみが合計されなければならない。上記のように、ファセットのサイズは、Ｖ_iα（ｘ）が少数のボクセルでのみ０でない値を有するように選択される。Ｖ_iα（ｘ）およびＰ_f（ｘ）は非整数値を有することがあるので、Γ_α（ｘ）は実数として格納および処理される。

３．ファセットからファセットへの移動
次に、粒子がファセット間で移動される（ステップ２８０）。ファセットＦ_αの流入状態（ｃ_iｎ_α＜０）の平行六面体Ｇ_iαが、別のファセットＦ_βと交差する場合、ファセットＦ_αが受け取った状態ｉの粒子の一部は、ファセットＦ_βから来ることになる。特に、ファセットＦ_αは、以前の時間増分の間にファセットＦ_βによって作り出された状態ｉの粒子の一部を受け取ることになる。

次に、図１７を参照すると、以前の時間増分の間にファセットＦ_βによって作り出された状態ｉの粒子の関係が示される。図１７において、ファセットＦ_βと交差する平行六面体Ｇ_iαの一部３８０が、ファセットＦ_αと交差する平行六面体Ｇ_iβの一部３８２と等しいことが示される。上記のように、交差部分は、Ｖ_iα（β）として表される。この項を使用して、ファセットＦ_βとファセットＦ_αとの間の状態ｉの粒子の流束は、
Γ_iα（β，ｔ－１）＝Γ_i（β）Ｖ_iα（β）／Ｖ_iα 式（Ｉ．１３）
として記述することができ、ここで、Γ_i（β，ｔ－１）は、以前の時間増分の間にファセットＦ_βによって作り出された状態ｉの粒子の分量である。この式から、ファセットＦ_α（ｃ_iｎ_α＜０）に向けられた状態ｉごとに、他のファセットによってファセットＦ_αに供給される粒子の数は、
Γ_iαF→F＝Σ_βΓ_iα（β）＝Σ_βΓ_i（β，ｔ－１）Ｖ_iα（β）／Ｖ_iα 式（Ｉ．１４）
であり、ファセットへの状態ｉの粒子の全流束は、
Γ_iIN（α）＝Γ_iαF→F＋Γ_iαF→F＝Σ_xＮ_i（ｘ）Ｖ_iα＋Σ_βΓ_i（β，ｔ－１）Ｖ_iα（β）／Ｖ_iα 式（Ｉ．１５）
である。

ファセット分布関数とも呼ばれるファセットの状態ベクトルＮ（α）は、ボクセル状態ベクトルのＭ個のエントリに対応するＭ個のエントリを有する。Ｍは、離散格子速度の数である。ファセット分布関数Ｎ（α）の入力状態は、それらの状態への粒子の流束を体積Ｖ_iαで除算したものに等しく設定され、ｃ_iｎ_α＜０に対して、
Ｎ_i（α）＝Γ_iIN（α）／Ｖ_iα 式（Ｉ．１６）
である。

ファセット分布関数は、ファセットからの出力流束を生成するためのシミュレーションツールであり、必ずしも実際の粒子を表していない。正確な出力流束を生成するために、分布関数の他の状態に値を割り当てる。内向き状態を投入するための上述の技法を使用して、外向き状態が投入され、ｃ_iｎ_α≧０では、
Ｎ_i（α）＝Γ_iOTHER（α）／Ｖ_iα 式（Ｉ．１７）
であり、ここで、Γ_iOTHER（α）は、Γ_iIN（α）を生成するための上述の技法を使用するが、流入状態（ｃ_iｎ_α＜０）以外の状態（ｃ_iｎ_α≧０）にこの技法を適用して決定される。代替の手法では、Γ_iOTHER（α）は、以前の時間ステップからのΓ_iOUT（α）の値を使用して生成することができ、その結果、
Γ_iOTHER（α，ｔ）＝Γ_iOUT（α，ｔ－１） 式（Ｉ．１８）
となる。

平行状態（ｃ_iｎ_α＝０）では、Ｖ_iαとＶ_iα（ｘ）は共に０である。Ｎ_i（α）の式において、Ｖ_iα（ｘ）は、分子に現われ（Γ_iOTHER（α）の式から）、Ｖ_iαは、分母に現われる（Ｎ_i（α）の式から）。その結果、平行状態のＮ_i（α）は、Ｖ_iαおよびＶ_iα（ｘ）が０に近づくときのＮ_i（α）の極限として決定される。０速度を有する状態（すなわち、静止状態、ならびに状態（０，０，０，２）および（０，０，０，－２））の値は、温度と圧力の初期条件に基づいてシミュレーションの初めに初期化される。次いで、これらの値は、時間とともに調節される。

４．ファセット表面動力学の実行
次に、表面動力学が、ファセットごとに、上述で論じた４つの境界条件を満たすように実行される（ステップ２８２）。ファセットに表面動力学を実行するための手順が、図１８に示される。最初に、ファセットＦ_αに垂直な合計運動量が、ファセットにおける粒子の合計運動量Ｐ（α）を決定することによって、すべてのｉに対して、

として決定される（ステップ３９２）。この式から、法線運動量Ｐ_n（α）は、
Ｐ_n（α）＝ｎ_α・Ｐ（α） 式（Ｉ．２０）
として決定される。

次いで、この法線運動量がプッシュ／プル技法を使用して除去されて（ステップ３９４）、Ｎ_n-（α）が作り出される。この技法によって、粒子は、法線運動量にのみ影響を与えるように状態間で移動される。プッシュ／プル技法は、米国特許第５，５９４，６７１号に説明されており、それは参照によって組み込まれる。

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

ファセットＦ_αの流出流束分布が、流入流束分布、ＣＦＬ制約違反を考慮に入れるための修正流束計算、およびボルツマン分布に基づいて、決定される（ステップ３９８）。

最初に、流入流束分布Γ_i（α）とボルツマン分布との間の差が、
ΔΓ_i（α）＝Γ_iIN（α）－Ｎ_n-βi（α）Ｖ_iα 式（Ｉ．２１）
として決定される。

この差を使用して、流出流束分布は、
ｎ_αｃ_i＞０に対して、
Γ_iOUT（α）＝Ｎ_n-βi（α）Ｖ_iα－．Δ．Γ_i*（α） 式（Ｉ．２２）
であり、ここで、ｉ^*は、状態ｉと反対の方向を有する状態である。例えば、状態ｉが（１，１，０，０）である場合、状態ｉ^*は（－１，－１，０，０）である。表面摩擦および他の要因を考慮するために、流出流束分布は、以下のようにさらに精密にすることができる。ｎ_αｃ_i＞０に対して、
Γ_iOUT（α）＝Ｎ_n-Bi（α）Ｖ_iα－ΔΓ_i*（α）＋
Ｃ_f（ｎ_α・ｃ_i）－［Ｎ_n-βi*（α）－Ｎ_n-βi（α）］Ｖ_iα＋
（ｎ_α・ｃ_i）（ｔ_1α・ｃ_i）ΔＮ_j,1Ｖ_iα＋
（ｎ_α・ｃ_i）（ｔ_2α・ｃ_i）ΔＮ_j,2Ｖ_iα 式（Ｉ．２３）
であり、ここで、Ｃ_fは、表面摩擦の関数であり、ｔ_iαは、ｎ_αに垂直な第１の接線ベクトルであり、ｔ_2αは、ｎ_αとｔ_1αの両方に垂直な第２の接線ベクトルであり、ΔＮ_j,1およびΔＮ_j,2は、状態ｉおよび表示された接線ベクトルのエネルギー（ｊ）に対応する分布関数である。分布関数は、

に従って決定され、ここで、ｊは、エネルギーレベル１の状態では１に等しく、エネルギーレベル２の状態では２に等しい。

Γ_iOUT（α）の式の各項の役割は、以下の通りである。第１および第２の項は、衝突がボルツマン分布を作り出すのに有効な範囲で法線運動量流束の境界条件を強制するが、接線運動量の流束異常を含む。第４および第５の項は、不十分な衝突に起因する離散の影響または非ボルツマン構造のために生じることがあるこの異常を補正する。最後に、第３の項は、表面の接線運動量流束の所望の変化を強制するために特定の量の表面分画を追加する。摩擦係数Ｃ_fの生成について以下で説明する。ベクトル操作を含むすべての項は、シミュレーションを始める前に計算することができる幾何学的因子であることに留意されたい。

この式から、接線速度は、
ｕ_i（α）＝（Ｐ（α）－Ｐ_n（α）ｎ_α）／ρ 式（Ｉ．２５）
として決定され、ここで、ρは、ファセット分布の密度

である。

前と同様に、流入流束分布とボルツマン分布との間の差は、
ΔΓ_i（α）＝Γ_iIN（α）－Ｎ_n-βi（α）Ｖ_iα 式（Ｉ．２７）
として決定される。

次いで、流出流束分布は、
Γ_iOUT（α）＝Ｎ_n-βi（α）Ｖ_iα－ΔΓ_i*（α）＋Ｃ_f（ｎ_αｃ_i）［Ｎ_n-βi*（α）－Ｎ_n-βi（α）］Ｖ_iα 式（Ｉ．２８）
になり、これは、以前の技法によって決定された流出流束分布の最初の２つのラインに対応するが、異常な接線流束への補正を必要としない。

いずれの手法を使用しても、結果として生じる流束分布は、運動量流束条件のすべてを満たし、すなわち、

であり、ここで、ｐ_αは、ファセットＦ_αの平衡圧力であり、ファセットに粒子を供給するボクセルの平均密度および温度値に基づき、ｕ_αは、ファセットにおける平均速度である。

質量およびエネルギー境界条件が確実に満たされるために、入力エネルギーと出力エネルギーとの間の差が、エネルギーレベルｊごとに、

として測定され、ここで、インデクスｊは、状態ｉのエネルギーを表す。次いで、このエネルギー差を使用して、差の項を生成し、ｃ_jiｎ_α＞０に対して、

である。この差の項を使用して、ｃ_jiｎ_α＞０に対して、
Γ_αjiOUTf＝Γ_αjiOUT＋δΓ_αji 式（Ｉ．３２）
になるように流出流束を修正する。この操作は、接線運動量流束を不変のままにしながら、質量およびエネルギー流束を補正する。この調節は、流れがファセットの近傍でほぼ一様で平衡に近い場合、わずかである。調節の後に結果として生じる法線運動量流束は、近傍の平均特性に、近傍の非均一または非平衡特性による補正を加えたものに基づいて、平衡圧力である値までわずかに変更される。ＣＦＬ制約が違反されている場合、プロセスは、修正流束計算手法を、図１０のプロセスに含まれるいずれかの流束計算に適用する２８５。

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

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

表面と相互作用しないボクセルでは、移動操作は、計算上極めて簡単である。状態の集団全体が、時間増分ごとの間に、現在のボクセルから目的地のボクセルに移動される。同時に、目的地のボクセルの粒子は、そのボクセルからそれ自体の目的地のボクセルに移動される。例えば、＋１ｘおよび＋１ｙ方向（１，０，０）に移動しているエネルギーレベル１の粒子は、現在のボクセルからｘ方向に＋１および他の方向に０にあるボクセルに移動される。粒子は、最後には、移動（１，０，０）の前に有していたものと同じ状態で目的地のボクセルに行き着く。ボクセル内の相互作用は、他の粒子および表面との局所相互作用に基づいて、その状態に対して粒子数を変更する可能性があることになる。そうでない場合、粒子は、同じ速度および方向で格子に沿って移動し続けることになる。

移動操作は、１つまたは複数の表面と相互作用するボクセルでは、わずかに複雑になる。これにより、１つまたは複数の分画粒子がファセットに移送されることになる。そのような分画粒子をファセットに移送することにより、分画粒子がボクセルに残ることになる。これらの分画粒子は、ファセットによって占められるボクセルに移送される。

図１６を参照すると、ボクセル３６２の状態ｉの粒子の一部３６０が、ファセット３６４に移動される（ステップ２７８）場合、残りの一部３６６は、ファセット３６４が配置されているボクセル３６８に移動され、そこから、状態ｉの粒子は、ファセット３６４に向けられる。したがって、状態集団が２５に等しく、Ｖ_iα（ｘ）が０．２５に等しい（すなわち、ボクセルの４分の１が平行六面体Ｇ_iαと交差する）場合、６．２５個の粒子がファセットＦ_αに移動されることになり、１８．７５個の粒子がファセットＦ_αによって占められたボクセルに移動されることになる。多数のファセットが単一のボクセルと交差するので、１つまたは複数のファセットによって占められたボクセルＮ（ｆ）に移送される状態ｉの粒子の数は、

であり、ここで、Ｎ（ｘ）はソースボクセルである。

６．ファセットからボクセルへの散乱
次に、各ファセットからの流出粒子がボクセルに散乱される（ステップ２８６）。本質的に、このステップは、粒子がボクセルからファセットに移動される収集ステップの逆である。ファセットＦ_αからボクセルＮ（ｘ）に移動する状態ｉの粒子の数は、

であり、ここで、Ｐ_f（ｘ）は、部分ボクセルの体積減少に相当する。この式から、状態ｉごとに、ファセットからボクセルＮ（ｘ）に向けられる粒子の総数は、

である。

ファセットからボクセルに粒子を散乱させ、それを、周囲のボクセルから移流した粒子と組み合わせ、その結果を整数化した後、特定のボクセルの特定の方向が、アンダーフローする（負になる）か、またはオーバーフローする（８ビット実装では２５５を超える）可能性がある。これは、これらの量が許容範囲の値に適合するように切り捨てられた後、質量、運動量、およびエネルギーの増加または減少のいずれかをもたらすことになる。そのような現象の発生を防ぐために、違反している状態を切り捨てる前に、境界を外れた質量、運動量、およびエネルギーは蓄積される。状態が属するエネルギーに対して、増加した値（アンダーフローのために）または減少した値（オーバーフローのために）に等しい量の質量は、同じエネルギーを有しかつそれ自体オーバーフローまたはアンダーフローの影響を受けないランダムに（または逐次に）選択された状態に追加される。この質量およびエネルギーの追加から生じる追加の運動量は、蓄積され、切捨てによる運動量に追加される。同じエネルギー状態にのみ質量を追加することによって、質量カウンタが０に達すると、質量とエネルギーの両方が補正される。最後に、運動量は、運動量アキュムレーターが０に戻るまでプッシュ／プル技法を使用して補正される。

７．流体動力学の実行
流体動力学が実行される（ステップ２８８）、図１０。このステップは、マイクロ動力学またはイントラボクセル操作と呼ぶこともできる。同様に、移流手順は、インターボクセル操作と呼ぶこともできる。以下で説明するマイクロ動力学操作を、さらに、ファセットで粒子を衝突させるのに使用して、ボルツマン分布を作り出すことができる。

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

可変分解能
可変分解能（米国特許出願公開第２０１３／０１５１２２１号に論じられているような）が、さらに、使用されてもよく、異なるサイズのボクセル、例えば、粗いボクセルおよび細かいボクセルを使用することになる。

固有過渡的格子ボルツマンベース物理学を活用することによって、システムは、現実の条件を正確に予測するシミュレーションを実行することができる。例えば、エンジニアは、変更の影響が設計および経費にとって最も重要である場合、プロトタイプが構築される前の設計プロセスにおいて製品性能を早期に評価する。システムは、ＣＡＤ幾何学を使用して、空気力学、空気－音響、および熱管理シミュレーションを正確に効率的に実行する。システムは、次のような用途に対処するためにシミュレーションを実行することができる。空気力学（空気力学効率；車両操作；汚れおよび水管理；パネル変形；運転動力学）、空力音響（温室風騒音；底部風騒音；ギャップ／シール騒音；ミラー、ホイッスル、およびトーナル騒音；サンルーフおよび窓バフェッティング；通過／ 都市騒音；冷却ファン騒音）、熱管理（冷却気流；熱保護；ブレーキ冷却；運転サイクルシミュレーション；キーオフおよびソーク；電子機器およびバッテリ冷却；ＲＯＡ ／吸気ポート）、空調（室内快適性；ＨＶＡＣユニット＆配電系統性能；ＨＶＡＣシステムおよびファン騒音；霜取りおよび解凍）、パワートレイン（ドライブトレイン冷却；排気システム；冷却ジャケット；エンジンブロック）、汚れおよび水管理（ピラーオーバーフロー、ちりおよびほこりの蓄積、タイヤスプレー）など。

本明細書に記載された主題および機能動作の実施形態は、デジタル電子回路、有形に具現されたコンピュータソフトウェアもしくはファームウェア、コンピュータハードウェア（本明細書に開示されている構造および構造等価物を含む）で、またはこれらのうちの１つまたは複数の組合せで実現することができる。本明細書に記載されている主題の実施形態は、１つまたは複数のコンピュータプログラム（すなわち、データ処理装置による実行のため、またはデータ処理装置の動作の制御のために有形の非一時的プログラムキャリアに符号化されたコンピュータプログラム命令の１つまたは複数のモジュール）として実現することができる。コンピュータ記憶媒体は、機械可読ストレージデバイス、機械可読ストレージサブストレート、ランダムもしくはシリアルアクセスメモリデバイス、またはそれらの１つまたは複数の組合せとすることができる。

「データ処理装置」という用語は、データ処理ハードウェアを参照し、例として、プログラマブルプロセッサ、コンピュータ、または多数のプロセッサもしくはコンピュータを含むデータを処理するためのすべての種類の装置、デバイス、および機械を包含する。装置は、専用論理回路（例えば、ＦＰＧＡ（フィールドプログラマブルゲートアレイ）またはＡＳＩＣ（特定用途向け集積回路））とすることもでき、またはそれをさらに含むことができる。ハードウェアに加えて、装置は、コンピュータプログラムのための実行環境を作り出すコード（例えば、プロセッサファームウェア、プロトコルスタック、データベース管理システム、オペレーティングシステム、またはそれらの１つまたは複数の組合せを構成するコード）をオプションとして含むことができる。

プログラム、ソフトウェア、ソフトウェアアプリケーション、モジュール、ソフトウェアモジュール、スクリプト、またはコードと呼ぶまたは記載することもできるコンピュータプログラムは、コンパイル型もしくはインタープリタ型言語、または宣言型もしくは手続き型言語を含む任意の形態のプログラミング言語で書くことができ、スタンドアロンプログラムとして、またはコンピューティング環境での使用に適したモジュール、コンポーネント、サブルーチン、もしくは別のユニットとしてを含む任意の形態で配備することができる。コンピュータプログラムは、ファイルシステムにおけるファイルに対応することができるが、対応する必要はない。プログラムは、他のプログラムまたはデータ（例えば、マークアップ言語ドキュメントに、当該のプログラムに専用の単一ファイルに、または多数の協調ファイル（例えば、１つまたは複数のモジュール、サブプログラム、もしくはコードの一部を格納するファイル）に格納された１つまたは複数のスクリプト）を保持するファイルの一部に格納することができる。コンピュータプログラムは、プログラムが、１つのコンピュータで実行されるか、または１つのサイトに配置されるかもしくは多数のサイトにわたって分散され、データ通信ネットワークによって相互接続された多数のコンピュータで実行されるように配備することができる。

コンピュータプログラムの実行に適したコンピュータは、汎用もしくは専用マイクロプロセッサまたは両方、あるいは他の種類の中央処理装置に基づくことができる。コンピュータプログラム命令およびデータを格納するのに適したコンピュータ可読媒体は、例として、半導体メモリデバイス（例えば、ＥＰＲＯＭ、ＥＥＰＲＯＭ、およびフラッシュメモリデバイス）、磁気ディスク（例えば、内部ハードディスクまたはリムーバブルディスク）、光磁気ディスク、およびＣＤ－ＲＯＭとＤＶＤ－ＲＯＭディスクを含む、媒体およびメモリデバイス上のすべての形態の不揮発性メモリを含む。プロセッサおよびメモリは、専用論理回路によって補足されるかまたは専用論理回路に組み込まれてもよい。

本明細書に記載されている主題の実施形態は、バックエンド構成要素（例えば、データサーバとしての）を含むか、またはミドルウェア構成要素（例えば、アプリケーションサーバ）を含むか、またはフロントエンド構成要素（例えば、ユーザが本明細書に記載されている主題の実施態様と対話することができるグラフィカルユーザインタフェースまたはウェブブラウザを有するクライアントコンピュータ）を含むコンピューティングシステムで、あるいは１つまたは複数のそのようなバックエンド、ミドルウェア、またはフロントエンド構成要素の任意の組合せで実施することができる。システムの構成要素は、任意の形態または媒体のデジタルデータ通信（例えば、通信ネットワーク）によって相互接続することができる。通信ネットワークの例は、ローカエリアネットワーク（ＬＡＮ）およびワイドエリアネットワーク（ＷＡＮ）（例えば、インターネット）を含む。

コンピューティングシステムは、クライアントおよびサーバを含むことができる。クライアントおよびサーバは、通常、互いに離れており、一般に、通信ネットワークを通して対話する。クライアントとサーバの関係は、それぞれのコンピュータ上で実行するとともに互いにクライアント－サーバ関係を有するコンピュータプログラムによって生じる。いくつかの実施形態では、サーバは、クライアントとして働くユーザデバイス（例えば、ユーザデバイスと対話するユーザにデータを表示し、ユーザからユーザ入力を受信するための）にデータ（例えば、ＨＴＭＬページ）を送信する。ユーザデバイスで生成されたデータ（例えば、ユーザ対話の結果）は、サーバにおいてはユーザデバイスから受信され得る。

主題の特定の実施形態が説明された。他の実施形態は、以下の特許請求の範囲内にある。例えば、特許請求の範囲で詳述される動作は、異なる順序で実行され、それでもなお、望ましい結果を達成することができる。１つの例として、添付図に示されているプロセスは、望ましい結果を達成するのに、図示の特定の順序または逐次的順序を必ずしも必要としない。場合によっては、マルチタスキングおよび並列処理が有利であることがある。

１０ システム
１２ サーバシステム
１４ クライアントシステム
１８ メモリ
２０ インタフェース
２２ 移流操作
２４ 処理デバイス
２８ メッシュ定義
３０ シミュレーションプロセス
３２ メッシュ準備エンジン
３４ シミュレーションエンジン
３４ａ 粒子衝突相互作用モジュール
３４ｂ 粒子境界モデルモジュール
３４ｃ 移流モジュール
３６ サブモジュール、エンジン
３８ データレポジトリ

Claims (20)

1. 物理的オブジェクトのまわりの流体流れをシミュレートするためのコンピュータ実施方法であって、前記方法は、
   １つまたは複数のコンピューティングシステムによって、ボクセルの集合として表された格子構造と、前記物理的オブジェクトの表現とを含むシミュレーション空間のモデルを受け取るステップであって、前記ボクセルが、前記物理的オブジェクトの表面を説明するための分解能を有する、受け取るステップと、
   前記１つまたは複数のコンピューティングシステムによって、流体の体積中の粒子の移動をシミュレートするステップであって、前記粒子の前記移動が前記粒子間の衝突を引き起こす、シミュレートするステップと、
   前記１つまたは複数のコンピューティングシステムによって、ボクセルの前記集合中の後続のボクセルへの前記粒子の移流のための時間ステップ値をテストするステップと、
   前記１つまたは複数のコンピューティングシステムによって、テストされた前記時間ステップ値に基づいて、２つのボクセル間の面を識別するステップであって、前記面のうちの少なくとも１つが安定条件に違反している、識別するステップと、前記安定条件に違反している前記面のうちの少なくとも１つに対して、
   前記コンピューティングシステムによって、前記２つのボクセルを含む領域の空間的に平均化された温度勾配を使用して、前記安定条件を満たす修正流束を計算するステップであって、前記面の前記少なくとも１つが前記安定条件に違反している、計算するステップと、
   前記コンピューティングシステムによって、前記修正流束を使用して後続のボクセルへの前記粒子の移流操作を実行するステップと、
   前記１つまたは複数のコンピューティングシステムによって、前記シミュレーションの結果を記憶または表示するステップと、
   を含む、コンピュータ実施方法。
2. 計算された前記修正流束は計算された修正熱流束である、請求項１に記載のコンピュータ方法。
3. 前記修正熱流束を計算するステップが、
   前記コンピューティングシステムによって、付加流束を計算するステップと、
   前記コンピューティングシステムによって、平衡流束を計算するステップと
   をさらに含む、請求項２に記載のコンピュータ方法。
4. 前記ボクセルの所定の１つに対して、前記計算された付加流束が、前記ボクセルの前記所定の１つのための温度発展を計算するために使用される、請求項３に記載のコンピュータ方法。
5. 前記平衡流束が、前記ボクセルのサイズに応じて前記温度発展の前記計算に使用される、請求項４に記載のコンピュータ方法。
6. 前記平衡流束は、前記ボクセルの前記サイズが制約を満たすほどに大きい場合、前記温度発展に使用される、請求項５に記載のコンピュータ方法。
7. 前記コンピューティングシステムによって、前記平衡流束を、流束流の方向の１つまたは複数の隣接するボクセルに送出するステップ
   をさらに含む、請求項４に記載のコンピュータ方法。
8. 物理的オブジェクトのまわりの流体流れをシミュレートするための装置であって、前記装置は、メモリ、１つまたは複数のプロセッサデバイスを含み、前記１つまたは複数のプロセッサデバイスは、
   ボクセルの集合として表された格子構造と、前記物理的オブジェクトの表現とを含むシミュレーション空間のモデルを受け取ることであって、前記ボクセルが、前記物理的オブジェクトの表面を説明するための分解能を有する、受け取ることと、
   流体の体積中の粒子の前記シミュレーション空間を通じた移動をシミュレートすることであって、前記粒子の前記移動が前記粒子間の衝突を引き起こす、シミュレートすることと、
   ボクセルの前記集合中の後続のボクセルへの前記粒子の移流のための時間ステップ値をテストすることと、
   テストされた前記時間ステップ値に基づいて、２つのボクセル間の面を識別することであって、前記面のうちの少なくとも１つが安定条件に違反している、識別することと、前記安定条件に違反している前記面のうちの少なくとも１つに対して、
   前記２つのボクセルを含む領域の空間的に平均化された温度勾配を使用して、前記安定条件を満たす修正流束を計算することであって、前記面の前記少なくとも１つが前記安定条件に違反している、計算することと、
   前記修正流束を使用して後続のボクセルへの前記粒子の移流操作を実行することと、
   前記シミュレーションの結果を記憶または表示することと、
   を実行するように構成されている、装置。
9. 計算された前記修正流束は計算された修正熱流束であり、前記条件は安定特性を含む、請求項８に記載の装置。
10. 前記装置は、前記修正熱流束を付加流束成分および平衡流束成分を計算することによって計算する、ようにさらに構成されている、請求項９に記載の装置。
11. 前記ボクセルの所定の１つに対して、前記計算された付加流束が、前記ボクセルの前記所定の１つのための温度発展を計算するために付加される、請求項１０に記載の装置。
12. 前記平衡流束を、流束流の方向の１つまたは複数の隣接するボクセルに送出する、ようにさらに構成されている、請求項１０に記載の装置。
13. 前記平衡流束を、前記平衡流束が十分に大きいボクセルに移送されるまで流束流の前記方向に沿って連続的に送出し、前記平衡流束は、前記十分に大きいボクセルでの温度発展に向けて適用される、請求項１２に記載の装置。
14. １つまたは複数の処理デバイスによって実行可能である、物理的オブジェクトのまわりの流体流れをシミュレートするための１つまたは複数の機械可読ハードウェアストレージデバイスに有形的に記憶されたコンピュータプログラムであって、前記コンピュータプログラムはコンピュータに、
    ボクセルの集合として表された格子構造と、前記物理的オブジェクトの表現とを含むシミュレーション空間のモデルを受け取ることであって、前記ボクセルが、前記物理的オブジェクトの表面を説明するための分解能を有する、受け取ることと、
    流体の体積中の粒子の前記シミュレーション空間を通じた移動をシミュレートすることであって、前記粒子の前記移動が前記粒子間の衝突を引き起こす、シミュレートすることと、
    ボクセルの前記集合中の後続のボクセルへの前記粒子の移流のための時間ステップ値をテストすることと、
    テストされた前記時間ステップ値に基づいて、２つのボクセル間の面を識別することであって、前記面のうちの少なくとも１つが安定条件に違反している、識別することと、前記安定条件に違反している前記面のうちの少なくとも１つに対して、
    前記２つのボクセルを含む領域の空間的に平均化された温度勾配を使用して、前記安定条件を満たす修正流束を計算することであって、前記面の前記少なくとも１つが前記安定条件に違反している、計算することと、
    前記修正流束を使用して後続のボクセルへの前記粒子の移流操作を実行することと、
    前記シミュレーションの結果を記憶または表示することと、
    を実行させるように構成されている、コンピュータプログラム。
15. 計算された前記修正流束は計算された修正熱流束であり、前記条件は安定特性を含む、請求項１４に記載のコンピュータプログラム。
16. 前記コンピュータに、前記修正熱流束を付加流束成分および平衡流束成分を計算することによって計算すること、をさらに実行させる、請求項１５に記載のコンピュータプログラム。
17. 前記ボクセルの所定の１つに対して、前記計算された付加流束が、前記ボクセルの前記所定の１つのための温度発展を計算するために付加される、請求項１６に記載のコンピュータプログラム。
18. 前記コンピュータに、前記平衡流束を、流束流の方向の１つまたは複数の隣接するボクセルに送出すること、をさらに実行させる、請求項１７に記載のコンピュータプログラム。
19. 前記コンピュータに、前記平衡流束を、前記平衡流束が十分に大きいボクセルに移送されるまで前記流束の前記方向に沿って連続的に送出し、前記平衡流束は、前記十分に大きいボクセルでの温度発展に向けて適用されること、をさらに実行させる、請求項１８に記載のコンピュータプログラム。
20. 前記条件が、安定分布を維持するために使用され得る最大時間ステップサイズを決定する時間進行スキームのＣｏｕｒａｎｔ－Ｆｒｉｅｄｒｉｃｈｓ－Ｌｅｗｙ（ＣＦＬ）制約である、請求項１７に記載のコンピュータプログラム。

Images