(57)【要約】      （修正有）【課題】固体の表面周囲の流量をシミュレートするための方法を提供する。【解決手段】方法は、固体の表面に適合する曲線メッシュを表す座標系を受け取るステップと、格子速度の組を用いて、一定量の流体中の粒子間の衝突を引き起こす粒子の輸送をシミュレートするステップと、粒子衝突の決定と、曲線メッシュに関連する粒子分布の変化とを含む、粒子の輸送の分布関数を実行するステップと、粒子運動量値に適用される制約下で座標系の移流動作を実行するステップと、座標系において決定された粒子運動量値及び空間座標を、曲線空間における運動量値及び空間値に変換することによって、シミュレーションから得られた値を曲線メッシュ上にマッピングするステップと、を含む。【選択図】図２
2020115343.tif 000002
【特許請求の範囲】
【請求項１】
固体の表面周囲の流体流をシミュレートするためのコンピュータ実装方法であって、
前記コンピュータシステムが、前記固体の前記表面に適合する曲線メッシュを表す座標系を受け取るステップと、
格子速度の組を用いて、一定量の流体中の粒子間の衝突を引き起こす前記粒子の輸送をシミュレートするステップと、
粒子衝突の決定と、前記曲線メッシュに関連する粒子分布の変化とを含む、前記粒子の輸送の分布関数を実行するステップと、
前記コンピュータシステムが、粒子運動量値に適用される制約下で前記座標系の移流動作を実行するステップと、
前記コンピュータシステムが、前記座標系において決定された前記粒子運動量値及び空間座標を前記曲線空間における運動量値及び空間値に変換することによって、シミュレーションから得られた値を前記曲線メッシュ上にマッピングするステップと、
を含むことを特徴とするコンピュータ方法。
【請求項２】
前記固体表面は、湾曲した固体表面である、
請求項１に記載のコンピュータ方法。
【請求項３】
前記メッシュは、ユークリッド空間における曲線メッシュである、
請求項１に記載のコンピュータ方法。
【請求項４】
前記座標系は、リーマン幾何学で定められた非ユークリッド空間である、
請求項１に記載のコンピュータ方法。
【請求項５】
前記一般的曲線座標系は、非ユークリッド空間におけるデカルト座標系である選択された座標系への１対１のマッピングを有する、
請求項１に記載のコンピュータ方法。
【請求項６】
粒子移流は、各時間ステップにおいて粒子分布関数が１つのボクセル（格子サイト）から別のボクセルに一意的に移流するように実行される、
請求項３に記載のコンピュータ方法。
【請求項７】
前記粒子分布の変化は、前記曲線メッシュにおける各時間増分において各ボクセルに加えられる力の印加をシミュレートして前記座標系における正確な運動量保存を実行することを前記コンピュータシステムに行わせる、
請求項１に記載のコンピュータ方法。
【請求項８】
一般的曲線メッシュにおける各ボクセルが規定の体積を有し、質量保存が実行されて質量及び運動量密度が規定されるように、体積表現が使用される、
請求項３に記載のコンピュータ方法。
【請求項９】
マッピングするステップは、前記メッシュ及び前記固体表面を非ユークリッド空間における前記デカルト座標系で表す結果として制限される運動量状態を決定するステップをさらに含む、
請求項１に記載のコンピュータ方法。
【請求項１０】
固体の表面周囲の流体流をシミュレートするための装置であって、
メモリと、
１又は２以上のプロセッサ装置と、
を備え、前記１又は２以上のプロセッサ装置は、
前記固体の前記表面に適合する曲線メッシュを表す座標系を受け取り、
格子速度の組を用いて、一定量の流体中の粒子間の衝突を引き起こす前記粒子の輸送をシミュレートし、
粒子衝突の決定と、前記曲線メッシュに関連する粒子分布の変化とを含む、前記粒子の輸送の分布関数を実行し、
粒子運動量値に適用される制約下で前記座標系の移流動作を実行し、
前記座標系において決定された前記粒子運動量値及び空間座標を前記曲線空間における運動量値及び空間値に変換することによって、シミュレーションから得られた値を前記曲線メッシュ上にマッピングする、
ように構成される、
ことを特徴とする装置。
【請求項１１】
前記固体表面は、湾曲した固体表面である、
請求項１０に記載の装置。
【請求項１２】
前記メッシュは、ユークリッド空間における曲線メッシュである、
請求項１０に記載の装置。
【請求項１３】
前記座標系は、リーマン幾何学で定められた非ユークリッド空間である、
請求項１０に記載の装置。
【請求項１４】
前記一般的曲線座標系は、非ユークリッド空間におけるデカルト座標系である選択された座標系への１対１のマッピングを有する、
請求項１０に記載の装置。
【請求項１５】
粒子移流は、各時間ステップにおいて粒子分布関数が１つのボクセル（格子サイト）から別のボクセルに一意的に移流するように実行される、
請求項１２に記載の装置。
【請求項１６】
前記粒子分布の変化は、前記曲線メッシュにおける各時間増分において各ボクセルに加えられる力の印加をシミュレートして前記座標系における正確な運動量保存を実行することを前記コンピュータシステムに行わせる、
請求項１０に記載の装置。
【請求項１７】
一般的曲線メッシュにおける各ボクセルが規定の体積を有し、質量保存が実行されて質量及び運動量密度が規定されるように、体積表現が使用される、
請求項１２に記載の装置。
【請求項１８】
マッピングするステップは、前記メッシュ及び前記固体表面を非ユークリッド空間における前記デカルト座標系で表す結果として制限される運動量状態を決定するステップをさらに含む、
請求項１０に記載の装置。
【請求項１９】
１又は２位以上の処理装置によって実行可能な命令を記憶した１又は２以上の機械可読ハードウェア記憶装置であって、前記命令は、コンピュータに、
前記固体の前記表面に適合する曲線メッシュを表す座標系を受け取ることと、
格子速度の組を用いて、一定量の流体中の粒子間の衝突を引き起こす前記粒子の輸送をシミュレートすることと、
粒子衝突の決定と、前記曲線メッシュに関連する粒子分布の変化とを含む、前記粒子の輸送の分布関数を実行することと、
粒子運動量値に適用される制約下で前記座標系の移流動作を実行することと、
前記座標系において決定された前記粒子運動量値及び空間座標を前記曲線空間における運動量値及び空間値に変換することによって、シミュレーションから得られた値を前記曲線メッシュ上にマッピングすることと、
を実行させる、
ことを特徴とする１又は２以上の機械可読ハードウェア記憶装置。
【請求項２０】
前記メッシュは、ユークリッド空間における曲線メッシュである、
請求項１９に記載の１又は２以上の機械可読ハードウェア記憶装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本明細書は、物理的流体流などの物理的プロセスのコンピュータシミュレーションに関する。
【背景技術】
【０００２】
いわゆる「格子ボルツマン法」（ＬＢＭ）は、計算流体力学において使用される有利な技術である。格子ボルツマンシステムの基礎を成す力学は、ボルツマン方程式に従う多くの粒子の動きを伴う気体分子運動論の基本物理学に存在する。ボルツマン基本運動システムには、衝突及び移流という２つの基本的な動力学プロセスがある。衝突プロセスは、保存則に従う粒子間の相互作用、及び平衡状態への緩和を伴う。移流プロセスは、微視的粒子速度に従う１つの位置から別の位置への粒子のモデリング動作を伴う。
【０００３】
標準的なＬＢＭモデルでは、粒子速度が一定値の離散集合を呈し、後者は、１つの格子サイトから、３次元（３Ｄ）均一立方体デカルトメッシュ（ｕｎｉｆｏｒｍ  ｃｕｂｉｃａｌ  Ｃａｒｔｅｓｉａｎ  ｍｅｓｈ）に対応する単純なブラベー格子上の近隣サイトへの正確なリンクを形成する。
【０００４】
任意座標系（任意メッシュ）に基づいてＬＢＭをメッシュに拡張する様々な試みが行われてきた。主な代表的方法の１つは、１対の格子サイト間の１対１の移流マッピングを弛緩させることである。このような任意メッシュ上では、一般にその元々のメッシュサイトからの移流後の粒子が単一の近隣サイトに付着しない。これらの解決策では、これらのメッシュサイトが補間によって粒子の位置を表す。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】米国特許出願公開第２０１６／０１８８７６８号明細書
【特許文献２】米国特許第５，５９４，６７１号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００６】
均一立方体デカルトメッシュは、多くの応用において重要ではあるが、他の応用に対して基本的限界を示す。例えば、実際の流体シミュレーションでは、シミュレーションが固体幾何学的曲面のものであることが多い。デカルトメッシュは、固体幾何学的曲面への滑らかな適合を示さない。また、通常、実際の物理的流れは、いくつかの空間領域及び方向に集中した小構造を有する。例えば、いわゆる乱流境界層では、壁に垂直なフロースケールが、接線方向又は流体領域の大部分よりもはるかに小さい。この結果、境界層内部の壁に対する垂直方向では、空間分解能の要件が著しく高い。立方体デカルトメッシュは、異なる方向の異なる空間分解能との柔軟性を提供しない。
【課題を解決するための手段】
【０００７】
１つの態様によれば、固体の表面周囲の流体流をシミュレートするためのコンピュータ実装方法が、コンピュータシステムが、固体の表面に適合する曲線メッシュを表す座標系を受け取るステップと、格子速度の組を用いて、一定量の流体中の粒子間の衝突を引き起こす粒子の輸送をシミュレートするステップと、粒子衝突の決定と、曲線メッシュに関連する粒子分布の変化とを含む、粒子の輸送の分布関数を実行するステップと、コンピュータシステムが、粒子運動量値に適用される制約下で座標系の移流動作を実行するステップと、コンピュータシステムが、座標系において決定された粒子運動量値及び空間座標を、曲線空間における運動量値及び空間値に変換することによって、シミュレーションから得られた値を曲線メッシュ上にマッピングするステップと、を含む。
【０００８】
他の態様は、コンピュータプログラム製品、１又は２以上の機械可読ハードウェア記憶装置、装置及びコンピュータシステムを含む。
【０００９】
本明細書に開示する方法は、デカルトメッシュ上で予測される現在のＬＢＭベースのシミュレーションを、曲線空間における非デカルトメッシュフレームワークに拡張する。本明細書に開示する方法は、体積定式化に基づくことによって、質量及び運動量の保存が満たされるようにする。結果として得られる質量の連続方程式は、曲線座標における正しい形を有し、従ってこの手法は、結果として得られる流体力学のアーチファクトを補正するために人工的な質量ソース項（ｍａｓｓ  ｓｏｕｒｃｅ  ｔｅｒｍｓ）を導入する必要がない。また、マニホルド上の連続気体分子運動論（ｃｏｎｔｉｎｕｕｍ  ｋｉｎｅｔｉｃ  ｔｈｅｏｒｙ）と同様に、開示する拡張ＬＢＭの外部ソース項（ｅｘｔｅｒｎａｌ  ｓｏｕｒｃｅ  ｔｅｒｍ）のみが、曲線空間に起因する慣性力の存在に起因する。この慣性力項（ｉｎｅｒｔｉａｌ  ｆｏｒｃｅ  ｔｅｒｍ）は質量に寄与せず、連続気体分子運動論における分析形式に依拠せずに構成される。
【００１０】
この離散空間及び時間における慣性力項は、連続気体分子運動論の力を流体力学的極限において漸近的に回復させる。さらに、この慣性力は、離散空間−時間ＬＢＭモデルにおける根本的なユークリッド空間の正確な運動量保存を実行する。この力項（ｆｏｒｃｅ  ｔｅｒｍ）は、結果として得られる正しいナビエ・ストーク流体力学を粘性順に（ａｔ  ｔｈｅ  ｖｉｓｃｏｕｓ  ｏｒｄｅｒ）生成するために、正しい離散時間モーメントにおいてシステムに運動量を加えるように構成される。
【００１１】
図面を含む以下の説明及び特許請求の範囲から、他の特徴及び利点が明らかになるであろう。
【図面の簡単な説明】
【００１２】
【図１】流体流シミュレーションシステムを示す図である。
【図２】曲線座標において格子ボルツマンモデルを形成する動作を示すフローチャートである。
【図３】シミュレートされる物理的対象に適合する記憶された曲線メッシュを使用した動作を示すフローチャートである。
【図４】曲線座標で表される格子ボルツマンモデルを使用したシミュレーション動作を示すフローチャートである。
【図５】分布関数に適用される制約動作を示すフローチャートである。
【図６Ａ】曲線空間及びデカルト空間で表される形状を有する物理的対象の表現を示す図である。
【図６Ｂ】曲線空間及びデカルト空間で表される形状を有する物理的対象の表現を示す図である。
【図７Ａ】曲線空間及びデカルト空間における粒子の動きの表現を示す図である。
【図７Ｂ】曲線空間及びデカルト空間における粒子の動きの表現を示す図である。
【図８Ａ】曲線空間及びデカルト空間における粒子の動きの別の表現を示す図である。
【図８Ｂ】曲線空間及びデカルト空間における粒子の動きの別の表現を示す図である。
【図９】非ユークリッド空間で表される２つのＬＢＭモデルの速度成分を示す図である。
【図１０】非ユークリッド空間で表される２つのＬＢＭモデルの速度成分を示す図である。
【図１１】物理的プロセスシミュレーションシステムが従う手順のフローチャートである。
【図１２】マイクロブロックの透視図である。
【図１３Ａ】図１のシステムが使用する非ユークリッド空間の格子構造の図である。
【図１３Ｂ】図１のシステムが使用する非ユークリッド空間の格子構造の図である。
【図１４】可変分解能技術を示す図である。
【図１５】可変分解能技術を示す図である。
【図１６】粒子の動きを示す図である。
【図１７】表面のファセットの影響を受ける領域を示す図である。
【図１８】表面動力学のフローチャートである。
【図１９】ボクセルから表面への粒子の動きを示す図である。
【図２０】表面動力学の実行手順のフローチャートである。
【図２１】流体シミュレーションを表すスクリーンショットである。
【発明を実施するための形態】
【００１３】
モデルシミュレーション空間
ＬＢＭベースの物理プロセスシミュレーションシステムでは、離散速度の集合ｃｉにおいて評価される分布関数値ｆｉによって流体流が表される。分布関数の動力学は、以下の方程式Ｉ．１によって制御される。
【数１】
2020115343.tif 000003
方程式（Ｉ．１）
【００１４】
この方程式は、分布関数ｆiの時間発展を表す周知の格子ボルツマン方程式である。左辺は、いわゆる「流動過程」による分布の変化を表す。流動過程は、あるメッシュ位置において流体ポケットが開始した後に、複数の速度ベクトルのうちの１つに沿って次のメッシュ位置に移動する時点である。この時点で、「衝突係数」、すなわち開始流体ポケットに対する近隣の流体ポケットの影響を計算する。流体は別のメッシュ位置にしか移動することできず、従って全ての速度の全ての成分が共通速度の倍数になるように正しい速度ベクトルの選択が必要である。
【００１５】
第１の方程式の右辺は、流体ポケット同士の衝突に起因する分布関数の変化を表す上述した「衝突演算子」である。特定の形の衝突演算子は、限定するわけではないが、Ｂｈａｔｎａｇａｒ、Ｇｒｏｓｓ及びＫｒｏｃｋ（ＢＧＫ）演算子とすることができる。衝突演算子は、「平衡」形態である以下の第２の方程式によって与えられる規定値に分布関数を至らしめる。
【数２】
2020115343.tif 000004
方程式（Ｉ．２）
【００１６】
ＢＧＫ演算子は、たとえ衝突の詳細がどうであろうと、分布関数は、衝突を通じて｛ｆeq（ｘ，ｖ，ｔ）｝によって与えられる明確に定義された局所平衡に近付くという物理的議論に従って構築され、
【数３】
2020115343.tif 000005
ここでのパラメータτは、衝突を介した平衡までの特性緩和時間を表す。
【００１７】
このシミュレーションから、質量ρ及び流体速度ｕなどの従来の流体変数が以下の方程式（Ｉ．３）の単純和として得られる。
【００１８】
対称性を考慮することにより、速度値の集合は、配置空間内でスパニングした時に特定の格子構造を形成するように選択される。このような離散系の動特性は、以下の形を有するＬＢＥ方程式に従い、
【数４】
2020115343.tif 000006
ここでの衝突演算子は、通常は上述したようなＢＧＫの形を取る。平衡分布形態を正しく選択することにより、格子ボルツマン方程式が正しい流体力学及び熱流体力学の結果をもたらすことを理論的に示すことができる。すなわち、ｆi（ｘ，ｔ）から導き出される流体力学的モーメントは、巨視的極限ではナヴィエ・ストークス式に従う。これらのモーメントは、以下のように定義され、
【数５】
2020115343.tif 000007
方程式（Ｉ．３）
ここでのρ及びｕは、それぞれ流体密度及び速度である。
【００１９】
ｃi及びｗiの集合的な値はＬＢＭモデルを規定する。ＬＢＭモデルは、スケーラブルコンピュータプラットフォーム上に効率的に実装して、時間非定常流（ｔｉｍｅ  ｕｎｓｔｅａｄｙ  ｆｌｏｗｓ）及び複雑な境界条件のために高いロバスト性を伴って実行することができる。
【００２０】
ボルツマン方程式から流体系の動きの巨視的方程式を取得する標準技術は、完全なボルツマン式の連続近似を取るチャップマン−エンスコッグ法（Ｃｈａｐｍａｎ−Ｅｎｓｋｏｇ  ｍｅｔｈｏｄ）である。流体系では、密度のわずかな乱れが音速で伝わる。気体系では、一般に音速が温度によって決まる。流れにおける圧縮率の影響の重要性は、特性速度と音速との比率によって測定され、これはマッハ数として知られている。
【００２１】
従来のＬＢＭベースの物理プロセスシミュレーションシステムのさらなる説明については、内容全体が引用により本明細書に組み入れられる米国特許出願公開第２０１６／０１８８７６８号を参照されたい。
【００２２】
図１に、物理的対象表現の周囲の流体流をシミュレートするシステム１０を示す。この実装のシステム１０は、クライアントサーバアーキテクチャ又はクラウドベースアーキテクチャに基づき、超並列コンピューティングシステム１２（スタンドアロン又はクラウドベース）として実装されるサーバシステム１２と、クライアントシステム１４とを含む。サーバシステム１２は、メモリ１８と、バスシステム１１と、インターフェイス２０（例えば、ユーザインターフェイス／ネットワークインターフェイス／ディスプレイ又はモニタインターフェイスなど）と、処理装置２４とを含む。メモリ１８内には、メッシュ準備エンジン３２及びシミュレーションエンジン３４が存在する。
【００２３】
図１では、メッシュ準備エンジン３２をメモリ１８内に示しているが、メッシュ準備エンジンは、サーバ１２とは異なるシステム上で実行されるサードパーティアプリケーションとすることもできる。メッシュ準備エンジン３２は、メモリ１８内で実行されるか、それともサーバ１２とは異なるシステム上で実行されるかにかかわらず、ユーザ指定のメッシュ定義３０を受け取り、メッシュを準備してシミュレーションエンジン３４に送信（及び／又は記憶）する。シミュレーションエンジン３４は、粒子衝突相互作用モジュールと、粒子境界モデルモジュールと、移流動作とを含む。システム１０は、２Ｄ及び／又は３Ｄメッシュ（デカルト座標及び／又は曲線）、座標系及びライブラリを記憶するデータリポジトリ３８にアクセスする。
【００２４】
次に、図２に、物理的対象表現の周囲の流体流をシミュレートするプロセス４０を示す。この例では、物理的対象が翼形部である。物理的対象はあらゆる形状とすることができ、具体的には（単複の）曲面を有するので、気流の使用は例示にすぎない。プロセスは、例えばクライアントシステム１４から、又はデータリポジトリ３８から検索を行うことにより、例えばシミュレートする物理的対象に適合する曲線メッシュなどのメッシュを受け取る（４２）。他の実施形態では、外部システム又はサーバ１２のいずれかが、シミュレートする物理的対象に適合するメッシュをユーザ入力に基づいて生成する。このプロセスは、検索した曲線メッシュから幾何学量を事前計算し（４４）、検索した曲線メッシュに対応する事前計算された幾何学量を用いて動的格子ボルツマンモデルシミュレーションを実行する（４６）。
【００２５】
次に、図３を参照して、曲線格子の構成５０及び事前計算される幾何学量の決定の態様を以下のように説明する。ｑによってｘ＝ｘ（ｑ）が一意的に定まるように（シミュレートする物理的対象に適合する）曲線座標系を選択し、ここでのｑ≡（ｑ1，ｑ2，ｑ3）は、パラメータ化された基本曲線の３つの非共線形合同に沿った座標値である。
【００２６】
システム１０は、空間内の全ての頂点位置が識別され指定されるように与えられた空間内のレイアウトを有する３Ｄ曲線メッシュから開始する（５２）。３次元（３Ｄ）ユークリッド空間内のいずれかの空間点をｘとする。空間点ｘは、曲線メッシュのサイト（すなわち、頂点）上にのみ定められる５４。曲線メッシュ（すなわち、「元々の」曲線メッシュ｛ｘ｝）が与えられると、メッシュ｛ｘ｝上の全ての空間位置ｘｉが識別される。
【００２７】
システム１０は、空間点ｘと曲線メッシュ上のサイト（ｑ）との間の１対１のマッピングである座標系｛ｑ｝を選択する（５６）。曲線メッシュ｛ｘ｝上のあらゆるサイトｘでは、ｘが一意的に定まるようにｘに関連する、ｘ＝ｘ（ｑ）である一意の値ｑが存在する。システム１０は、以下のようにしてメッシュ上にｑの座標値を構築する。
【００２８】
曲線メッシュ上のあらゆるサイトｘ（ｑ）では、正又は負の方向のα番目の（α＝１，２，３）座標曲線に沿ったサイトの最も近い近隣サイトが空間点ｘ±iである。この近隣サイトの一意の座標値をｑ±αとする一意のマッピングｘ±α＝ｘ（ｑ±α）により、ｑ±α及びｑが、α番目の座標成分値のみが一定距離（ｄ０）だけ異なり、ｑ±α＝（ｑ1±α，ｑ2±α，ｑ3±α）かつｑβ±α−ｑβ＝±ｄ0δαβ（式中、α，β＝１，２，３）であるように、ｑの空間変動を選択することが完全に可能である。この例では、定数ｄ０が、一般性の喪失を伴わない単位元（１）となるように選択される。
【００２９】
この構成の下では、座標値｛ｑ｝が、単位元（ｄ０の値）の格子面間隔の単純な均一３Ｄ立方体デカルト格子構造を提供した。この変形した「デカルト」格子｛ｑ｝は、ユークリッド空間における元々の曲線メッシュ｛ｘ｝の変形（屈曲、捻れ及び伸び／圧縮）によってもたらされる。従って、「デカルト」格子｛ｑ｝の位相構造は元々の曲線メッシュ｛ｘ｝と同じであるが、結果として得られるデカルト格子は非ユークリッド空間上に存在する。
【００３０】
曲線メッシュが提供されると、メッシュ上の全ての頂点｛ｘ｝の空間位置が指定され、曲線メッシュ上のいずれか１つの頂点から別の頂点までの距離も完全に決定される。サイトｘ（ｑ）からサイトの近隣のうちの１つであるｘ（ｑ±α）までの距離ベクトルＤ±α（ｑ）（α＝１；２；３）は、以下のように定義される。
【数６】
2020115343.tif 000008
（方程式１）
例えば、
【数７】
2020115343.tif 000009
である。
【００３１】
一般に、１つのメッシュサイトからその最も近い近隣サイトまでの空間距離は、一般的な曲線メッシュの空間的非均一性によって位置毎に変化する。換言すれば、Ｄ±α（ｑ）はｑの関数である。さらに、一般にα番目の座標曲線に沿った正の方向の距離値は、負方向の距離値に等しくない。
【００３２】
距離ベクトルの項では、Ｄα（ｑ）≠−Ｄ-α（ｑ）であることが明らかである。例えば、（方程式１）によって与えられる定義によれば、
【数８】
2020115343.tif 000010
方程式２
となり、曲線メッシュが均一なデカルト格子である時にはどこでも（空間座標値ｑにかかわらず｜Ｄα｜＝Ｃｏｎｓｔとなるように）不等式が等式に変化するだけである。
【００３３】
システム１０は、選択された曲線メッシュの構築された座標系を記憶することができる（５８）。
【００３４】
システム１０は、基本微分幾何学の概念に従って各座標方向に沿って接線基底ベクトルを構築する（６０）。
【数９】
2020115343.tif 000011
（方程式３）
【００３５】
このような構成では、ｇα（ｑ）＝−ｇ-α（ｑ）であるように空間反転対称性が達成される。ユークリッド空間のデカルト座標系とは異なり、一般に曲線座標系の基底接線ベクトルｇα（ｑ）（α＝１，２，３）は正規直交ではない。すなわち、
【数１０】
2020115343.tif 000012
である。
【００３６】
従って、システム１０は、上記の基底接線ベクトルに基づいて定義される以下のような対応する計量テンソルと、
【数１１】
2020115343.tif 000013
（方程式４）
以下のようなｘ（ｑ）を中心とするセルの体積Ｊとを構築し（６２）、
【数１２】
2020115343.tif 000014
（方程式５）
均一なブラベー格子の空間内の定数であるＪ（ｑ）がＪ（ｑ）＞０となるように正しい「掌性（ｈａｎｄｅｄｎｅｓｓ）」を選択する。（本明細書で使用する掌性とは、通常はベクトル分析において参照される、結果として得られるベクトルの方向慣習（典型的には右手の法則）を意味する。）以下の式を検証することができ、
【数１３】
2020115343.tif 000015
（方程式６）
ｄｅｔ［ｇαβ（ｑ）］は、計量［ｇαβ（ｑ）］テンソルの決定因子である。共接線基底ベクトルｇα（ｑ）（α＝１，２，３）を以下のように構築する（６４）。
【数１４】
2020115343.tif 000016
（方程式７）
【００３７】
基底接線ベクトル及び共接線ベクトルは互いに正規直交であり、ここでのδαβはクロネッカーのデルタ関数である。
【００３８】
同様に、逆計量テンソルは以下のように定義され、
【数１５】
2020115343.tif 000017
（方程式８）
逆計量テンソルは、計量テンソルの逆数であり、［ｇαβ（ｑ）］＝［ｇαβ（ｑ）］-1、又は、
【数１６】
2020115343.tif 000018
かつ、
【数１７】
2020115343.tif 000019
である。
【００３９】
上記で基本的幾何学量を定義したら、標準的なデカルト格子上の速度ベクトルと同様に、一般的な曲線メッシュ上に格子ボルツマン速度ベクトルを導入する（６６）。
【数１８】
2020115343.tif 000020
（方程式９）
【００４０】
定数ｃiαは、正の整数、負の整数又はゼロのいずれかであり、３次元座標配列のα番目の成分値ｃi≡（ｃi1，ｃi2，ｃi3）である。例えば、いわゆるＤ３Ｑ１９では、ギリシャ文字の添え字が０から１８にまで及ぶ。
【数１９】
2020115343.tif 000021
【００４１】
正しい完全なナビエ・ストーク流体力学を回復するために、モーメント等方性（ｍｏｍｅｎｔ  ｉｓｏｔｒｏｐｙ）及び標準化条件の組を満たす。これらは、正しい一定加重の組｛ωiｉ＝１，…ｂ｝が存在する時に、格子成分ベクトルの組が最大第６桁のモーメント等方性を認め、すなわち、
【数２０】
2020115343.tif 000022
（方程式１０）
であり、ここでのＴ0は一定温度であり、格子ベクトルの組の選択に依存し、δαβはクロネッカーデルタ関数である。なお、３次元配列はｃiであるのに対し、ｅi（ｑ）は方程式９の格子ベクトルである。
【００４２】
ＬＢＭモデルと共に使用する特定の幾何学量の組を以下のように定義し、
【数２１】
2020115343.tif 000023
（方程式１１）
【数２２】
2020115343.tif 000024
メッシュが均一なデカルト格子である時には、Θβα（ｑ＋ｃi，ｑ）の項がゼロになる。
【００４３】
従って、曲線メッシュが指定されると、上記の全ての幾何学量を完全に決定し、シミュレートする物理的対象に適合する曲線メッシュを記憶することができ（６８）、幾何学量が完全に決定されるので、動的ＬＢＭシミュレーションの開始前にこれらを事前計算することができる（図２を参照）。
【００４４】
図４を参照すると、シミュレーションプロセス４６（図２）は、例えば曲線空間のために適合された修正格子ボルツマン方程式（ＬＢＥ）に従って粒子分布の発展をシミュレートする（８０）。プロセスは、ＬＢＭ非ユークリッド空間内の次のセルｑへの粒子の移流を実行し（８２）、運動量及び空間座標を非ユークリッド空間内の座標系からユークリッド空間内の曲線系に変換することによって、結果として得られた値をシミュレーションから元々の曲線メッシュ上にマッピングし（８４）、結果として得られたマッピングされたユークリッド空間内の値をディスプレイ装置などにおいてレンダリングする（８６）。
【００４５】
図５に、シミュレーションプロセス４６のさらなる詳細を示す。シミュレーションプロセス４６（図２）は、例えば曲線空間のために適合された修正格子ボルツマン方程式（ＬＢＥ）（方程式１２を参照）に従って粒子分布の発展をシミュレートする（８０）（図４）。プロセスは、ＬＢＭ非ユークリッド空間内の次のセルｑへの粒子の移流を実行する（８２）（図４）。プロセスは、粒子移流に制約を課す（８３）。制約プロセス８３は、一般的な曲線メッシュの曲率及び非均一性に関連する有効外力に起因する粒子分布の変化を制限し（８３ａ）、運動量フラックスに制約を課し（８３ｂ）、曲線座標における正しいオイラー方程式及びナビエ・ストーク方程式を流体力学的極限において回復させる規定の（８３ｃ）平衡分布関数を使用する。従って、プロセスは、運動量及び空間座標を非ユークリッド空間内の座標系からユークリッド空間内の曲線系に変換することによって、結果として得られた値をシミュレーションから元々の曲線メッシュ上にマッピングし（８４）（図４）、ユークリッド空間におけるマッピングされた結果としての値をディスプレイ装置などにレンダリングする８６（図４）。
【００４６】
次に、図６Ａに、ユークリッド空間（現実空間）における大まかな３Ｄ曲線メッシュ９３内に翼形部９２を示しており、（明確にするために）２つの次元のみを示す。図示のように、メッシュ９３は、一般に翼形部９２の形状、具体的には外面に適合する。適合の度合いは、図６Ａに示すものよりも高く又は低くすることができる。
【００４７】
図６Ｂには、非ユークリッド空間におけるデカルトメッシュ９５内に翼形部９２（図６Ａ）の表現９２’を示す（２つの次元のみを示す）。この例では、（説明のみを目的として）翼形部９２の表現９２’を長方形固体として示している。図示のように、デカルトメッシュ９５は、一般に翼形部の表現９２’の形状、具体的には外面に適合する。適合の度合いは、図示のものよりも高く又は低くすることができる。
【００４８】
次に、図７Ａ、図７Ｂに、ユークリッド空間（現実空間ではあるが、明確にするために２つの次元しか示していない）における大まかな３Ｄ曲線メッシュ９３内に、流体の粒子の運動量を表すベクトル９６を示す。図示のように、ベクトル９６は、曲線空間９６を直線的に進む流体の粒子を表す。しかしながら、変形デカルト非ユークリッド空間（ｍｏｒｐｈｅｄ  Ｃａｒｔｅｓｉａｎ  ｎｏｎ−Ｅｕｃｌｉｄｅａｎ  ｓｐａｃｅ）では、粒子を表すベクトル９６’が曲線的な動きを行っているように見える。
【００４９】
次に、図８Ａ、図８Ｂに、図７Ａ、図７Ｂのものとは逆に、ユークリッド空間（現実空間ではあるが、明確にするために２つの次元しか示していない）における大まかな３Ｄ曲線メッシュ９３内に、流体の粒子の運動量を表すベクトル９７を示す。図示のように、ベクトル９７は、曲線空間９６を曲線的な動きで進む流体の粒子を表す。しかしながら、変形デカルト非ユークリッド空間では、粒子を表すベクトル９７’が直線的な動きを行っているように見える。次に、これらのプロセスの特定の詳細について説明する。
【００５０】
曲線格子上の体積格子ボルツマンモデル
説明する容量格子ボルツマンモデル手法は、一般に様々な定式化に応用可能であるが、説明例のためにいわゆる「等温ＬＢＭ」の定式化について説明する。
【００５１】
修正分布は、粒子分布の発展の標準分布と同様である。曲線格子上の体積格子ボルツマンモデルは、以下のように提供され、
【数２３】
2020115343.tif 000025
（方程式１２）
ここでのＮi（ｑ，ｔ）は、一般性の喪失を伴わない単位元増分の時間ｔにおけるセルｑの離散方向ｃiに属する粒子の数である。（方程式１２）における項Ωα（ｑ，ｔ）は、局所質量及び運動量保存方程式を満たす衝突項である。
【数２４】
2020115343.tif 000026
（方程式１３）
【００５２】
式１２における余分項（ｅｘｔｒａ  ｔｅｒｍ）δΝi（ｑ，ｔ）は、一般的な曲線メッシュの曲率及び非均一性に関連する有効外力に起因する粒子分布の変化を表す。この余分項は、デカルト格子上の標準的ＬＢＭではゼロになる。
【００５３】
粒子密度分布関数は、以下の式によってΝi（ｑ，ｔ）に関連し、
【数２５】
2020115343.tif 000027
方程式１４
ここでのＪ（ｑ）は、以前に定義されたｑを中心とするセルの体積である。基本流体量は、標準的な流体力学的モーメントによって与えられ、
【数２６】
2020115343.tif 000028
（方程式１５）
ここでのρ（ｑ，ｔ）及びｕ（ｑ，ｔ）は、位置ｑ及び時間ｔにおける流体密度及び速度である。（方程式９）の関係を使用すると、上記の速度モーメントは以下のように書き直され、
【数２７】
2020115343.tif 000029
（方程式１６）
曲線座標系における速度成分値は以下の式によって与えられる。
【数２８】
2020115343.tif 000030
（方程式１７）
【００５４】
表記を単純にするために、３次元流体速度配列をＵ（ｑ，ｔ）≡（Ｕ1（ｑ，ｔ），Ｕ2（ｑ，ｔ）Ｕ3（ｑ，ｔ））として定義し、式１７を同等に以下のように表す。
【数２９】
2020115343.tif 000031
（方程式１８）
【数３０】
2020115343.tif 000032
（方程式１９）
【００５５】
方程式１８は、標準的なデカルト格子ベースのＬＢＭと同じ流体速度の形である。同様に、方程式１３における衝突項の運動量保存は、以下のように書くことができる。
【数３１】
2020115343.tif 000033
（方程式２０）
【００５６】
ＬＢＭでは、衝突項が線形形態を取ることが多く、Ｍij及びｆi（ｑ，ｔ）は、それぞれ衝突行列及び平衡分布関数を表す。特に、いわゆるＢｈａｔｎａｇａｒ‐Ｇｒｏｓｓ‐Ｋｒｏｏｋ（ＢＧＫ）形態は以下に対応し、
【数３２】
2020115343.tif 000034
τは衝突緩和時間である。正しいナビエ・ストーク流体力学を回復させるには、方程式１３及び１９の他に、衝突行列がさらなる条件を満たす必要がある。
【数３３】
2020115343.tif 000035
（方程式２１）
【００５７】
ＢＧＫ形態は、このようなさらなる特性を満たす。上述したように、方程式１２における余分項δΝi（ｑ，ｔ）は、一般的な曲線メッシュの曲率及び非均一性に関連する有効外力に起因する粒子分布の変化を表す。この余分項は、デカルト格子上の標準的ＬＢＭではゼロになる。移流プロセスは、標準的なＬＢＭにおける曲線メッシュ内の１つのサイトから別のサイトへの正確な１対１のホップ、すなわち、
【数３４】
2020115343.tif 000036
（方程式２２）
であり、ここでのＮi’（ｑ，ｔ）は、（ｑ，ｔ）における衝突後の分布である。曲線メッシュに起因して、セルｑから移流する粒子の量は、セルｑ＋ｃiに到達するものに全く等しい（方程式２２を参照）が、運動量は移流中に変化する。
【００５８】
一般に、
【数３５】
2020115343.tif 000037
である。
【００５９】
上記の式では、不等号の左辺が移流プロセスの最後の運動量値であり、右辺がプロセス開始時の値である。不等式が存在する理由は、曲線メッシュに起因して粒子の経路が曲線である（しかも伸長又は圧縮されている）ことによって、移流の最後の速度がその元々の値から変化するようになるためである。このことは、移流プロセス全体を通じて粒子が一定速度を有する均一なデカルト格子とは基本的に異なる。
【００６０】
この結果、全体的な運動量値に以下の不等式が存在し、
【数３６】
2020115343.tif 000038
（方程式２３）
この場合、（２３）における不等号の右辺は、全ての近隣セルから移流した総運動量を表し、左辺は、曲線経路に沿って移流後にセルｑに到達した総運動量である。従って方程式２２及び方程式２３から、全ての近隣セルからセルｑ内への移流を通じた正味の運動量変化が以下の式によって与えられる。
【数３７】
2020115343.tif 000039
（方程式２４）
【００６１】
同様に、セルｑからその全ての近隣セルへの移流を通じた正味の総運動量変化が以下の式によって与えられる。
【数３８】
2020115343.tif 000040
（方程式２５）
【００６２】
その後、以下のδΝi（ｑ，ｔ）に制約が課せられた場合、
【数３９】
2020115343.tif 000041
（方程式２６）
根本的なユークリッド空間について質量保存が維持され、正しい運動量保存が回復する。ここでは、
【数４０】
2020115343.tif 000042
【数４１】
2020115343.tif 000043
である。
【００６３】
具体的に言えば、方程式２６における第１の制約は、 δΝi（ｑ，ｔ）によって導入される質量ソースを有していない。方程式２６における第２の制約は、根本的なユークリッド空間のあらゆる格子サイトｑ及び時間ｔにおいて運動量保存を満たす量に完全に等しい「慣性力」を導入する。このメカニズムは、曲がった空間における連続体運動の理論的説明（ｃｏｎｔｉｎｕｕｍ  ｋｉｎｅｔｉｃ  ｔｈｅｏｒｅｔｉｃ  ｄｅｓｃｒｉｐｔｉｏｎ）に類似する。前の小区分からの記号を直接配置したＦα（ｑ，ｔ）＝χ（ｑ，ｔ）・ｇα（ｑ）としての座標成分形態での書き込みは、以下をもたらし、
【数４２】
2020115343.tif 000044
この場合、（方程式１１）において幾何学関数θβα（ｑ＋ｃi，ｑ）Ｎi’（ｑ，ｔ）が定義される。（２７）及び（１１）から、従来のＬＢＭのように曲線メッシュが規則的な均一デカルト格子である場合には、Ｆα（ｑ，ｔ）がゼロになることが分かる。（２６）における第２の制約は、以下のような座標成分の形で表すことができる。
【数４３】
2020115343.tif 000045
（方程式２８）
【００６４】
完全粘性ナビエ・ストーク方程式（ｆｕｌｌ  ｖｉｓｃｏｕｓ  Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ  ｅｑｕａｔｉｏｎ）を回復させるために、運動量フラックスに対するさらなる制約が以下のように課せられる。
【数４４】
2020115343.tif 000046
（方程式２９）
【数４５】
2020115343.tif 000047
（方程式３０）
【００６５】
【数４６】
2020115343.tif 000048
の特定の形は以下のように選択される。
【数４７】
2020115343.tif 000049
（方程式３１）
【００６６】
この方程式は、方程式２６、２８及び２９のモーメント制約を満たす。曲線座標の正しいオイラー方程式及びナビエ・ストーク方程式を流体力学的極限において回復させるために、平衡分布関数の形を定義する。特に、以下の流体力学モーメントの基礎的条件は以下の通りであり、
【数４８】
2020115343.tif 000050
（方程式３２）
上記では、
【数４９】
2020115343.tif 000051
である。
【００６７】
これらの基礎的条件は、以下の平衡分布形によって満たされる。
【数５０】
2020115343.tif 000052
（方程式３３）
【００６８】
上記の平衡分布形は、マクスウェル・ボルツマン分布の低マッハ数拡張（ｌｏｗ  Ｍａｃｈ  ｎｕｍｂｅｒ  ｅｘｐａｎｓｉｏｎ）に類似するが、曲線座標で表される。実際に、曲線メッシュが均一デカルト格子であってｇαβ= δαβである場合、上記の平衡分布形は、標準的なＬＢＭ平衡分布まで低減する。上述した全ての量及び動力学では、（非ユークリッド）均一デカルト格子｛ｑ｝についてシミュレートした格子ボルツマン方程式１２が、曲線座標を有するナビエ・ストーク流体力学に従う。従って、結果として得られる値は、以下の単純な変換によって元々の曲線メッシュにマッピングすることができる。
【数５１】
2020115343.tif 000053
（方程式３４）
【００６９】
非ユークリッド速度空間におけるデカルト座標
ユークリッド空間における元々の曲線メッシュ｛ｘ｝の変形（屈曲、捻れ及び伸長／圧縮）によって生じる変形「デカルト」格子｛ｑ｝は、上述したように一般的曲線メッシュをユークリッド空間に戻すために、非ユークリッド空間の位相構造における移流時に「デカルト」格子｛ｑ｝の制約が粒子運動量に加えられる限り、物理的物体（ｐｈｙｓｉｃａｌ  ｂｏｄｉｅｓ）の周囲のシミュレーション流体流において従来のＬＢＭのためのデカルト格子（ｘ）と同様に使用することができる。
【００７０】
図９を参照すると、第１のモデル（２Ｄ−１）１００は、２１個の速度を含む２次元モデルである。これらの２１個の速度のうち、１つ（１０５）は動いていない粒子を表し、３組の４つの速度は、格子のｘ軸又はｙ軸のいずれかに沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）（１１０〜１１３）、正規化速度の２倍（２ｒ）（１２０〜１２３）、又は正規化速度の３倍（３ｒ）（１３０〜１３３）のいずれかで動いている粒子を表し、２組の４つの速度は、ｘ格子軸及びｙ格子軸の両方に関して正規化速度（ｒ）（１４０〜１４３）又は正規化速度の２倍（２ｒ）（１５０〜１５３）で動いている粒子を表す。
【００７１】
また、図１０に示すように、第２のモデル（３Ｄ−１）２００は、各速度が図９の矢印のうちの１つによって表される３９個の速度を含む３次元モデルである。これらの３９個の速度のうち、１つは動いていない粒子を表し、３組の６つの速度は、格子のｘ軸、ｙ軸又はｚ軸に沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）のいずれかで動いている粒子を表し、８つの速度は、ｘ、ｙ、ｚ格子軸の３つ全てに関して正規化速度（ｒ）で動いている粒子を表し、１２個の速度は、ｘ、ｙ、ｚ格子軸の２つに関して正規化速度の２倍（２ｒ）で動いている粒子を表す。
【００７２】
１０１個の速度を含む３Ｄ−２モデル及び３７個の速度を含む２Ｄ−２モデルなどのさらに複雑なモデルを使用することもできる。
【００７３】
１０１個の速度から成る３次元モデル３Ｄ−２では、１つが動いていない粒子（グループ１）を表し、３組の６つの速度が、格子のｘ、ｙ又はｚ軸に沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）のいずれかで動いている粒子（グループ２、４及び７）を表し、３組の８つの速度が、ｘ、ｙ、ｚ格子軸の３つ全てに関して正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）で動いている粒子（グループ３、８及び１０）を表し、１２個の速度が、ｘ、ｙ、ｚ格子軸の２つに関して正規化速度の２倍（２ｒ）で動いている粒子（グループ６）を表し、２４個の速度が、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度（ｒ）及び正規化速度の２倍（２ｒ）で動いていて残りの軸に関しては動いていない粒子（グループ５）を表し、２４個の速度が、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度（ｒ）で動いていて残りの軸に関して正規化速度の３倍（３ｒ）で動いている粒子（グループ９）を表す。
【００７４】
３７個の速度から成る２次元モデル２Ｄ−２では、１つが動いていない粒子（グループ１）を表し、３組の４つの速度が、格子のｘ軸又はｙ軸のいずれかに沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）のいずれかで動いている粒子（グループ２、４及び７）を表し、２組の４つの速度が、ｘ及びｙ格子軸の両方に関して正規化速度（ｒ）又は正規化速度の２倍（２ｒ）で動いている粒子を表し、８つの速度が、ｘ及びｙ格子軸の一方に関して正規化速度（ｒ）で動いていて他方の軸に関して正規化速度の２倍（２ｒ）で動いている粒子を表し、８つの速度が、ｘ及びｙ格子軸の一方に関して正規化速度（ｒ）で動いていて他方の軸に対して正規化速度の３倍（３ｒ）で動いている粒子を表す。
【００７５】
上述したＬＢＭモデルは、２次元及び３次元の両方における流れの数値シミュレーションのための特定のクラスの効率的でロバストな離散速度動力学モデルを提供する。この種のモデルは、特定の離散速度の組と、これらの速度に関連する重みとを含む。これらの速度は、特に格子ボルツマンモデルとして知られている種類の離散速度モデルの正確で効率的な実装を容易にする速度空間における（非ユークリッド空間内）デカルト座標のグリッド点と一致する。このようなモデルを用いて、高い忠実度で流れをシミュレートすることができる。
【００７６】
図１１を参照すると、物理プロセスシミュレーションシステムが、手順３００に従って動作して流体流などの物理プロセスをシミュレートする。シミュレーションの前に、シミュレーション空間をボクセルの集合としてモデル化する（ステップ３０２）。通常、シミュレーション空間は、コンピュータ支援設計（ＣＡＤ）プログラムを用いて生成される。例えば、ＣＡＤプログラムを用いて、風洞内に位置するマイクロデバイスを描くことができる。その後、ＣＡＤプログラムによって生成されたデータを処理して適切な分解能の格子構造を追加し、シミュレーション空間内の物体及び表面を考慮する。
【００７７】
格子の分解能は、シミュレートされるシステムのレイノルズ数に基づいて選択することができる。レイノルズ数は、流れの粘度（ν）、流れにおける物体の特性長（Ｌ）及び流れの特性速度（ｕ）に関連する。
Ｒｅ＝ｕＬ／ｖ                                          方程式（Ｉ．５）
【００７８】
物体の特性長は、物体の大規模な特徴を表す。例えば、マイクロデバイスの周囲の流れをシミュレートする場合には、マイクロデバイスの高さを特性長と見なすことができる。小さな物体領域（例えば、自動車のサイドミラー）の周囲の流れに関心がある時には、シミュレーションの分解能を高めるか、或いは関心領域の周囲に高分解能の領域を使用することができる。格子の分解能が高まるにつれてボクセルの次元は減少する。
【００７９】
状態空間はｆi（ｘ，ｔ)として表され、ここでのｆｉは、時点ｔにおいて３次元ベクトルｘによって示される、ある格子サイトにおける状態ｉの単位体積当たりの要素又は粒子の数（すなわち、状態ｉの粒子の密度）を表す。既知の時間増分では、粒子の数が単純にｆi（ｘ）として示される。格子サイトの全ての状態の組み合わせは、ｆ（ｘ）として示される。
【００８０】
状態の数は、各エネルギーレベル内の考えられる速度ベクトル数によって決まる。速度ベクトルは、ｘ、ｙ及びｚの３次元を有する空間内の整数の線形速度から成る。多重種シミュレーションでは状態の数が増加する。
【００８１】
各状態ｉは、特定のエネルギーレベル（すなわち、エネルギーレベル０、１又は２）における異なる速度ベクトルを表す。各状態の速度ｃｉは、以下のように３次元の各次元における「速度」と共に示される。
ｃi＝（ｃix，ｃiy，ｃiz）                                  方程式（Ｉ．６）
【００８２】
エネルギーレベルゼロの状態は、どの次元においても動いていない停止した粒子を表し、すなわちｃ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α（ｘ）                                      方程式（Ｉ．１１）
【００９５】
Ｄ．シミュレーションの実行
１又は２以上のファセットの影響を受けるボクセルが識別されると（ステップ３０４）、タイマを初期化してシミュレーションを開始する（ステップ３０６）。シミュレーションの各時間増分中には、粒子と表面ファセットとの相互作用を考慮した移流段階（ステップ３０８〜３１６）によってボクセルからボクセルへの粒子の移動をシミュレートする。次に、上述したエントロピーソルバの使用に部分的に基づいて、衝突段階（ステップ３１８）によって各ボクセル内の粒子の相互作用をシミュレートする。その後、タイマを増分する（ステップ３２０）。増分したタイマによってシミュレーションが完了した旨が示されない（ステップ３２２）場合には、移流段階と衝突段階と（ステップ３０８〜３２０）を繰り返す。増分したタイマによってシミュレーションが完了した旨が示された（ステップ３２２）場合には、シミュレーションの結果を記憶及び／又は表示する（ステップ３２４）。
【００９６】
１．表面の境界条件
表面との相互作用を正しくシミュレートするために、各ファセットは４つの境界条件を満たさなければならない。まず、ファセットが受け取った粒子の合計質量が、ファセットが移動させた粒子の合計質量に等しくなければならない（すなわち、ファセットへの正味質量流束がゼロに等しくなければならない）。次に、ファセットが受け取った粒子の合計エネルギーが、ファセットが移動させた粒子の合計エネルギーに等しくなければならない（すなわち、ファセットへの正味エネルギー流束がゼロに等しくなければならない）。これらの２つの条件は、各エネルギーレベル（すなわち、エネルギーレベル１及び２）での正味質量流束がゼロに等しくなるよう求めることによって満たすことができる。
【００９７】
他の２つの境界条件は、ファセットと相互作用する粒子の正味運動量に関連する。本明細書では滑り面と呼ぶ表面摩擦がない表面では、正味接線運動量流束がゼロに等しくなければならず、正味法線運動量流束がファセットの局所圧力に等しくなければならない。従って、ファセットの法線ｎαに垂直な受け取った合計運動量の成分と、ファセットの法線ｎαに垂直な移動させた合計運動量の成分と（すなわち、接線成分）が等しくなければならず、ファセットの法線ｎαに平行な受け取った合計運動量の成分と、ファセットの法線ｎαに平行な移動させた合計運動量の成分と（すなわち、法線成分）の差分がファセットの局所圧力に等しくなければならない。非滑り面では、ファセットが受け取った粒子の合計接線運動量に対するファセットが移動させた粒子の合計接線運動量が、表面の摩擦によって摩擦量に関連する倍数だけ減少する。
【００９８】
２．ボクセルからファセットへの収集
粒子と表面との間の相互作用をシミュレートする最初のステップとして、ボクセルから粒子を集めてファセットに提供する（ステップ３０８）。上述したように、ボクセルＮ（ｘ）とファセットＦαとの間の状態ｉの粒子の流束は以下の通りである。
Γ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ｎα≧０）に適用して決定される。別の方法では、以前の時間ステップからのΓiOUT（α）の値を用いて以下のようにΓiOTHER（α）を生成することもできる。
ΓiOTHER（α，ｔ）＝ΓiOUT（α，ｔ−１）                方程式（Ｉ．１９）
【０１０４】
平行状態（ｃiｎα＝０）では、Ｖiα及びＶiα（ｘ）がいずれもゼロである。Ｎi（α）を求める式では、Ｖiα（ｘ）が（ΓiOTHER（α）の式からの）分子に現れ、Ｖiαが（Ｎi（α）の式からの）分母に現れる。従って、Ｖiα及びＶiα（ｘ）がゼロに近付くと、平行状態のＮi（α）は、Ｎi（α）の極限として決定される。ゼロ速度を有する状態の値（すなわち休止状態、並びに状態（０，０，０，２）及び（０，０，０，−２））は、シミュレーションの最初に温度及び圧力の初期条件に基づいて初期化される。その後、これらの値は時間をかけて調整される。
【０１０５】
４．ファセット表面動力学の実行
次に、各ファセットが上述した４つの境界条件を満たすように表面動力学を実行する（ステップ３１２）。図１８に、ファセットの表面動力学の実行手順を示す。最初に、ファセットにおける粒子の合計運動量Ｐ（α）を以下のように求めることによって、ファセットＦαに垂直な合計運動量を求める（ステップ１１０５）。
全てのｉについて、
【数５２】
2020115343.tif 000054
方程式（Ｉ．２０）
この式から、以下のように法線運動量Ｐn（α）を求める。
Ｐn（α）＝ｎα・Ｐ（α）                                  方程式（Ｉ．２１）
【０１０６】
次に、プッシング／プリング法を用いてこの法線運動量を排除（ステップ１１１０）してＮα（α）を取り出す。この手法によれば、粒子は、法線運動量にのみ影響を与える形で状態間を移行する。このプッシング／プリング法は、引用により組み入れられる米国特許第５，５９４，６７１号に記載されている。
【０１０７】
その後、Ｎn-（α）の粒子を衝突させてボルツマン分布Ｎn-β（α）を生じさせる（ステップ１１１５）。以下で流体力学の実行に関して説明するように、ボルツマン分布は、Ｎn-（α）に一連の衝突則を適用することによって達成することができる。
【０１０８】
次に、流入流束分布とボルツマン分布とに基づいてファセットＦαの流出流束分布を求める（ステップ１１２０）。まず、流入流束分布Γi（α）とボルツマン分布との差分を以下のように求める。
ΔΓi（α）＝ΓiIN（α）−Ｎn-βi（α）Ｖiα                方程式（Ｉ．２２）
【０１０９】
この差分を使用すると、流出流束分布は、
ｎαｃi＞０の場合、以下のようになり、
ΓiOUT（α）＝Ｎn-β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αは、ｎαとｔiαの両方に垂直な第２の接線ベクトルであり、ΔＮj,1及びΔＮj,2は、状態ｉのエネルギー（ｊ）及び指示される接線ベクトルに対応する分布関数である。分布関数は、次式に従って決定され、
【数５３】
2020115343.tif 000055
方程式（Ｉ．２５）
ここでのｊは、エネルギーレベル１状態では１に等しく、エネルギーレベル２状態では２に等しい。
【０１１０】
ΓiOUT（α）の方程式の各項の関数は以下の通りである。第１項及び第２項は、ボルツマン分布を生じる上で衝突が効果的である限りは法線運動量流束境界条件を強制するが、接線運動量流束の例外を含む。第４項及び第５項は、この不十分な衝突に起因する離散性効果又は非ボルツマン構造を原因として生じ得る例外を補正する。最後に、第３項は、表面上の接線運動量流束の所望の変化を強制するように一定量の表面摩擦を追加する。摩擦係数Ｃfの生成については以下で説明する。なお、ベクトル操作に関与する全ての項は、シミュレーションの開始前に計算できる幾何学的因子である。
【０１１１】
この式から、接線速度が以下のように求められ、
ｕi（α）＝（Ｐ（α）−Ｐn（α）ｎα）／ρ                    方程式（Ｉ．２６）
ここでのρは、以下のようなファセット分布の密度である。
【数５４】
2020115343.tif 000056
方程式（Ｉ．２７）
【０１１２】
上記と同様に、流入流束分布とボルツマン分布との差分を以下のように求める。
ΔΓi（α）＝ΓiIN（α）−Ｎn-βi（α）Ｖiα              方程式（Ｉ．２８）
【０１１３】
すると、流出流束分布は以下のようになり、
ΓiOUT（α）＝Ｎn-βi（α）Ｖiα−ΔΓi*（α）＋Ｃf（ｎαｃi）［Ｎn-βi*（α）−Ｎn-βi（α）］Ｖiα                                      方程式（Ｉ．２９）
この式は、先の手法によって求めた流出流束分布の最初の２行に対応するが、異常な接線流束のための補正を必要としない。
【０１１４】
いずれかの方法を使用すると、結果として得られる流束分布は運動量流束条件を全て満たし、すなわち、
【数５５】
2020115343.tif 000057
方程式（Ｉ．３０）
となり、ここでのｐαは、ファセットＦαにおける平衡圧であって、ファセットに粒子を提供するボクセルの平均密度及び温度値に基づき、ｕαは、ファセットにおける平均速度である。
【０１１５】
質量及びエネルギー境界条件が確実に満たされるように、各エネルギーレベルｊについて入力エネルギーと出力エネルギーとの差分を以下のように測定する。
【数５６】
2020115343.tif 000058
方程式（Ｉ．３１）
ここでの添字ｊは、状態ｉのエネルギーを示す。次に、このエネルギー差を用いて差分項（ｄｉｆｆｅｒｅｎｃｅ  ｔｅｒｍ）を生成すると、
ｃjiｎα＞０の場合、以下のようになる。
【数５７】
2020115343.tif 000059
方程式（Ｉ．３２）
この差分項を用いて流出流束を修正すると、流速は、
ｃjiｎα＞０の場合、以下のようになる。
ΓαjiOUTf＝ΓαjiOUT＋δΓαji                          方程式（Ｉ．３３）
この演算は、接線運動量流束をそのままにした状態で質量及びエネルギー流束を訂正する。ファセットの近傍において流れがおおよそ均一であって平衡に近い場合には、この調整はわずかなものである。結果として得られる調整後の法線運動量流束は、近傍の平均特性に近傍の非均一性又は非平衡特性に起因する補正をプラスしたものに基づく平衡圧である値へとわずかに変化する。
【０１１６】
５．ボクセルからボクセルへの移動
再び図３を参照すると、３次元直線格子に沿ってボクセル間で粒子を移動させる（ステップ３１４）。このボクセルからボクセルへの移動は、ファセットと相互作用しないボクセル（すなわち、表面付近に存在しないボクセル）に対して行われる唯一の移動操作である。典型的なシミュレーションでは、表面と相互作用するほど十分に表面近くに存在しないボクセルがボクセルの大多数を構成する。
【０１１７】
この分離状態の各々は、３次元のｘ、ｙ及びｚの各々において格子に沿って整数速度で移動する粒子を表す。整数速度は、０、±１及び±２を含む。速度の符号は、対応する軸に沿った粒子の移動方向を示す。
【０１１８】
表面と相互作用しないボクセルの場合、この移動操作は計算的に極めて単純である。毎時間増分中に、状態の集団全体をその現在のボクセルから宛先ボクセルに移動させる。同時に、宛先ボクセルの粒子を、そのボクセルから独自の宛先ボクセルに移動させる。例えば、＋１ｘ及び＋１ｙ方向に移動しているエネルギーレベル１の粒子（１，０，０）は、その現在のボクセルから、ｘ方向に＋１上回り他の方向が０であるボクセルに移動させる。この粒子は、最後には移動前の状態と同じ状態の宛先ボクセル（１，０，０）に行き着く。ボクセル内の相互作用は、他の粒子及び表面との局所的相互作用に基づいてその状態の粒子総数を変化させることがある。変化しなければ、粒子は格子に沿って同じ速度で同じ方向に移動し続ける。
【０１１９】
この移動操作は、１又は２以上の表面と相互作用するボクセルでは若干複雑になる。この結果、１又は２以上の端数粒子（ｆｒａｃｔｉｏｎａｌ  ｐａｒｔｉｃｌｅｓ）がファセットに移動することがある。このような端数粒子のファセットへの移動により、ボクセルに端数粒子が留まるようになる。これらの端数粒子を、ファセットに占有されたボクセルに移動させる。
【０１２０】
図１９を参照すると、ボクセル９０５の状態ｉの粒子の一部９００がファセット９１０に移動した時には（ステップ３０８）、残り部分９１５をファセット９１０が位置するボクセル９２０に移動させ、ここから状態ｉの粒子をファセット９１０に向ける。従って、状態集団が２５に等しく、Ｖiα（ｘ）が０．２５に等しい（すなわち、ボクセルの４分の１が平行六面体Ｇiαに交わる）場合には、６．２５の粒子がファセットＦαに移動し、１８．７５の粒子を、ファセットＦαに占有されたボクセルに移動させる。１つのボクセルに複数のファセットが交わることもあるので、１又は２以上のファセットに占有されたボクセルに移動する状態ｉの粒子の数Ｎ（ｆ）は以下のようになり、
【数５８】
2020115343.tif 000060
方程式３４
ここでのＮ（ｘ）はソースボクセルである。
【０１２１】
６．ファセットからボクセルへの散乱
次に、各ファセットからの流出粒子をボクセルに散乱させる（ステップ３１６）。基本的に、このステップは、ボクセルからファセットに粒子が移動する収集ステップの逆である。ファセットＦαからボクセルＮ（ｘ）に移動する状態ｉの粒子の数は以下のようになり、
【数５９】
2020115343.tif 000061
方程式３５
ここでのＰf（ｘ）は、部分的なボクセルの体積縮小に相当する。この式から、各状態ｉについて、ファセットからボクセルに向かう粒子の総数Ｎ(x)は以下のようになる。
【数６０】
2020115343.tif 000062
方程式３６
【０１２２】
粒子をファセットからボクセルに散乱させ、これらの粒子を周囲のボクセルから移流した粒子と組み合わせて結果を整数化した後には、一定のボクセルの一定の方向がアンダーフロー（負になる）又はオーバーフロー（８ビット実装において２５５を超える）のいずれかになることがある。この結果、これらの量を許容できる値の範囲に収まるように切り捨てた後に、質量、運動量及びエネルギーが増加又は減少するようになる。このような発生を防ぐために、有害な状態を切り捨てる前に、限界を超えた質量、運動量及びエネルギーを蓄積する。状態を有するエネルギーでは、（アンダーフローに起因して）得られた又は（オーバーフローに起因して）失われた値に等しい質量が、同じエネルギーを有するランダムに（又は順次に）選択された状態に再び追加され、これ自体がオーバーフロー又はアンダーフローに影響されることはない。この質量及びエネルギーの追加によって生じるさらなる運動量を蓄積して、切り捨てからの運動量に追加する。同じエネルギー状態のみに質量を追加することにより、質量カウンタがゼロに達した時に質量及びエネルギーがいずれも補正される。最後に、運動量積算器がゼロに戻るまでプッシング／プリング法を用いて運動量を補正する。
【０１２３】
７．流体力学の実行
流体力学を実行する（図１１、３１８）。このステップは、マイクロ動力学操作又はボクセル内操作と呼ぶことができる。同様に、移流手順は、ボクセル間操作と呼ぶことができる。後述するマイクロ動力学操作を用いて、ファセットにおける粒子を衝突させてボルツマン分布を生成することもできる。
【０１２４】
流体力学は、格子ボルツマン方程式モデルにおいて、ＢＧＫ衝突モデルとして知られている特定の衝突演算子によって保証される。この衝突モデルは、実際の流体系における分布の動力学を模倣する。衝突プロセスは、方程式Ｉ．１及び方程式Ｉ．２の右辺によって十分に説明することができる。移流ステップ後には、方程式Ｉ．３を用いた分布関数から流体系の保存量、具体的には密度、運動量及びエネルギーが取得される。これらの量から、方程式（Ｉ．２）にｆeqで示す平衡分布関数が方程式（Ｉ．４）によって十分に規定される。いずれも表１に列挙する速度ベクトルセットｃｉ、重みの選択は、方程式Ｉ．２と相まって、巨視的挙動が正しい流体動力学方程式に従うことを保証する。
【０１２５】
Ｅ．可変分解能
図２０を参照すると、（上述した）可変分解能は、以下では粗ボクセル１２０００及び微小ボクセル１２０５と呼ぶ異なるサイズのボクセルを採用する。（以下の説明では、２つの異なるサイズを有するボクセルを参照するが、説明する技術は、３又は４以上の異なるサイズのボクセルに適用してさらなる分解能レベルを提供することもできると理解されたい。）粗ボクセル及び微小ボクセルの領域間の界面は、可変分解能（ＶＲ）界面１２１０と呼ばれる。
【０１２６】
可変分解能を表面又はその付近で採用すると、ファセットがＶＲ界面の両側のボクセルに相互作用することができる。これらのファセットは、ＶＲ界面ファセット１２１５（ＦαIC）又はＶＲ微小ファセット１２２０（ＦαIF）として分類される。ＶＲ界面ファセット１２１５は、ＶＲ界面の粗側に位置するファセットであり、微小ボクセル内に延びる粗平行六面体１２２５を有する。（粗平行六面体は、粗ボクセルの寸法に従ってｃｉが寸法決めされるものであり、微小平行六面体は、微小ボクセルの寸法に従ってｃｉが寸法決めされるものである。）ＶＲ微小ファセット１２２０は、ＶＲ界面の微小側に位置するファセットであり、粗ボクセル内に延びる微小平行六面体１２３０を有する。界面ファセットに関連する処理は、粗ファセット１２３５（ＦαＣ）及び微小ファセット１２４０（ＦαF）との相互作用を伴うこともできる。
【０１２７】
両タイプのＶＲファセットでは、表面動力学が微小スケールで実行され、上述したように動作する。しかしながら、ＶＲファセットは、ＶＲファセットとの間を粒子が移流する方法に関して他のファセットと異なる。
【０１２８】
図２１に、流体シミュレーションのスクリーンショットを示す。（従来のデカルト手法、又は場合によっては曲線メッシュのための他の手法ではなく）上述した曲線法を用いた流体シミュレーションは、同様の流体シミュレーション表現と、いずれかの慣習的な対応する計算データ出力とをもたらす。しかしながら、上述した曲線法を用いたこのような流体シミュレーションは、例えば実際の物理的対象などの曲面を有する物体をモデル化する際に、他の方法よりも高速に、及び／又は少ない計算リソースを用いて流体シミュレーションを実行することができる。
【０１２９】
複数の実装について説明した。それでもなお、本開示の趣旨及び範囲から逸脱することなく様々な修正を行うことができると理解されるであろう。従って、以下の特許請求の範囲には他の実装も含まれる。
【図１】
2020115343.tif 000063
【図２】
2020115343.tif 000064
【図３】
2020115343.tif 000065
【図４】
2020115343.tif 000066
【図５】
2020115343.tif 000067
【図６Ａ】
2020115343.tif 000068
【図６Ｂ】
2020115343.tif 000069
【図７Ａ】
2020115343.tif 000070
【図７Ｂ】
2020115343.tif 000071
【図８Ａ】
2020115343.tif 000072
【図８Ｂ】
2020115343.tif 000073
【図９】
2020115343.tif 000074
【図１０】
2020115343.tif 000075
【図１１】
2020115343.tif 000076
【図１２】
2020115343.tif 000077
【図１３Ａ】
2020115343.tif 000078
【図１３Ｂ】
2020115343.tif 000079
【図１４】
2020115343.tif 000080
【図１５】
2020115343.tif 000081
【図１６】
2020115343.tif 000082
【図１７】
2020115343.tif 000083
【図１８】
2020115343.tif 000084
【図１９】
2020115343.tif 000085
【図２０】
2020115343.tif 000086
【図２１】
2020115343.tif 000087
【外国語明細書】
2020115343000001.pdf