知財求人 - 知財ポータルサイト「IP Force」
▶ ダッソー システムズ アメリカス コーポレイションの特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2024-05-29
(45)【発行日】2024-06-06
(54)【発明の名称】全エネルギー保存を実施する格子ボルツマンソルバ
(51)【国際特許分類】
G16Z 99/00 20190101AFI20240530BHJP
【ＦＩ】
G16Z99/00
【請求項の数】
18
【外国語出願】
(21)【出願番号】P 2020002949
(22)【出願日】2020-01-10
(65)【公開番号】P2020123325
(43)【公開日】2020-08-13
【審査請求日】2022-12-28
(31)【優先権主張番号】62/790,528
(32)【優先日】2019-01-10
(33)【優先権主張国・地域又は機関】US
(73)【特許権者】
【識別番号】514180812
【氏名又は名称】ダッソー システムズ アメリカス コーポレイション
(74)【代理人】
【識別番号】100094569
【弁理士】
【氏名又は名称】田中 伸一郎
(74)【代理人】
【識別番号】100103610
【弁理士】
【氏名又は名称】▲吉▼田 和彦
(74)【代理人】
【識別番号】100109070
【弁理士】
【氏名又は名称】須田 洋之
(74)【代理人】
【識別番号】100067013
【弁理士】
【氏名又は名称】大塚 文昭
(74)【代理人】
【識別番号】100120525
【弁理士】
【氏名又は名称】近藤 直樹
(74)【代理人】
【識別番号】100139712
【弁理士】
【氏名又は名称】那須 威夫
(74)【代理人】
【識別番号】100141553
【弁理士】
【氏名又は名称】鈴木 信彦
(74)【代理人】
【識別番号】100086771
【弁理士】
【氏名又は名称】西島 孝喜
(74)【代理人】
【氏名又は名称】上杉 浩
(72)【発明者】
【氏名】プラディープ ゴパラクリシュナン
(72)【発明者】
【氏名】フードン チェン
(72)【発明者】
【氏名】ラオヤン チャン
【審査官】宮地 匡人
(56)【参考文献】
【文献】特表２０１５－５０７７７８（ＪＰ，Ａ）
【文献】特表２０１６－５２８６２８（ＪＰ，Ａ）
【文献】特表２０１６－５３４４３６（ＪＰ，Ａ）
【文献】米国特許第０６９１５２４５（ＵＳ，Ｂ１）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ１６Ｚ ９９／００
ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）
ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）
(57)【特許請求の範囲】
【請求項1】
全エネルギー保存を実施するシミュレーションを使用してコンピュータ上で流体流をシミュレートする方法であって、メッシュにわたる粒子の動きをモデル化するように前記メッシュにわたる流体の活動をシミュレートするステップと、
前記メッシュ内の各メッシュ位置の、それぞれが対応するメッシュ位置における可能な運動量状態のうちの特定の運動量状態に対応する複数のエントリを含む一連の状態ベクトルをコンピュータアクセス可能メモリに記憶するステップと、
前記コンピュータが、移流する粒子の状態に、Ｅを比内部エネルギー、νを速度として、比全エネルギーＥ t ＝Ｅ＋ν 2 ／２を計算することによって得られる比全エネルギーの値を加算し、時間間隔にわたって移流しない粒子の状態から前記比全エネルギーの値を減算することによって、前記粒子の前記一連の状態ベクトルを修正するステップと、
を含むことを特徴とする方法。
【請求項2】
前記修正された一連の状態に従って前記粒子が後続のメッシュ位置に移流した後に、前記方法は、前記粒子状態に局所圧力項を加え直してからモーメントを計算して正しい圧力勾配（∇ｐ）をもたらすステップをさらに含む、請求項１に記載の方法。
【請求項3】
前記局所圧力項は、前記コンピュータによって計算されるθ－１項を含む、請求項２に記載の方法。
【請求項4】
前記流体流の活動をシミュレートするステップは、第１の一連の離散格子速度に部分的に基づいて前記流体流をシミュレートするステップを含み、前記方法は、前記第１の一連の離散格子速度と同じ格子速度である、或いは前記第１の一連の離散格子速度とは格子速度の数が異なる第２の一連の離散格子速度に部分的に基づいて、スカラー量の時間発展をシミュレートするステップをさらに含む、請求項２に記載の方法。
【請求項5】
前記コンピュータが、後続のメッシュ位置への前記粒子の移流を一定の時間間隔にわたって実行するステップ、をさらに含み、前記全エネルギーの時間発展をシミュレートするステップは、
近隣メッシュ位置からの流入分布を衝突演算子及びエネルギー演算子について収集するステップと、
前記流入分布に重み付けするステップと、
前記衝突演算子とエネルギー演算子との積として流出分布を決定するステップと、
前記決定された流出分布を伝播するステップと、
を含む、請求項１に記載の方法。
【請求項6】
任意のプラントル数のための正確な剪断応力を回復させるために、エネルギー衝突項は以下によって与えられ、【数1】
【数2】
請求項１に記載の方法。
【請求項7】
前記流入分布が、前記決定された流出分布に等しくなるように、ゼロの正味表面流束境界条件を適用するステップをさらに含む、請求項５に記載の方法。
【請求項8】
活動をシミュレートするステップは、前記コンピュータが、平衡分布及び比全エネルギーＥtの第２の分布関数を適用するステップを含み、前記第２の分布関数は、流れ分布ｆiと共に移流する特定のスカラーとして定義される、請求項１に記載の方法。
【請求項9】
前記第２の分布関数は、エネルギー方程式に対するｆi
の非平衡寄与を考慮して、境界付近の異なる格子分解能にわたる正確な流れ挙動を取得し、分布関数の衝突演算子は以下によって与えられ、【数3】
【数4】
【数5】
【数6】
【請求項10】
実行時に物理的プロセス流体流をシミュレートするとともに、コンピュータシステムに、メッシュにわたる粒子の動きをモデル化するように前記メッシュにわたる流体の活動をシミュレートし、
前記メッシュ内の各メッシュ位置の、それぞれが対応するメッシュ位置における可能な運動量状態のうちの特定の運動量状態に対応する複数のエントリを含む一連の状態ベクトルをコンピュータメモリに記憶し、
移流する粒子の状態に、Ｅを比内部エネルギー、νを速度として、比全エネルギーＥ t ＝Ｅ＋ν 2 ／２を計算することによって得られる比全エネルギーの値を加算し、時間間隔にわたって移流しない粒子の状態から前記比全エネルギーの値を減算することによって、前記粒子の前記一連の状態ベクトルを修正する、
ことを行わせる、ことを特徴とするコンピュータプログラム。
【請求項11】
前記比全エネルギーの値は、圧力が対流するように移流前に追加され、追加された前記比全エネルギーの値を補正するために、停止状態から前記追加された比全エネルギーの値を削除する、請求項１０に記載のコンピュータプログラム。
【請求項12】
前記追加された比全エネルギーの値を削除することによって全エネルギーが保存され、正しい圧力速度項がもたらされる、請求項１１に記載のコンピュータプログラム。
【請求項13】
前記追加された比全エネルギーの値を削除することは、前記停止状態【数7】
【数8】
請求項１２に記載のコンピュータプログラム。
【請求項14】
物理的プロセス流体流をシミュレートするためのコンピュータシステムであって、１又は２以上のプロセッサ装置と、
前記１又は２以上のプロセッサ装置に結合されたコンピュータメモリと、
命令を記憶したコンピュータ可読媒体と、
を備え、前記命令は、実行時に物理的プロセス流体流をシミュレートするとともに、前記コンピュータシステムに、
メッシュにわたる粒子の動きをモデル化するように前記メッシュにわたる流体の活動をシミュレートし、
前記メッシュ内の各メッシュ位置の、それぞれが対応するメッシュ位置における可能な運動量状態のうちの固有の運動量状態に対応する複数のエントリを含む一連の状態ベクトルをコンピュータメモリに記憶し、
移流する粒子の状態に、Ｅを比内部エネルギー、νを速度として、比全エネルギーＥ t ＝Ｅ＋ν 2 ／２を計算することによって得られる比全エネルギーの値を加算し、時間間隔にわたって移流しない粒子の状態から前記比全エネルギーの値を減算することによって、前記粒子の前記状態ベクトルを修正する、
ことを行わせる、
ことを特徴とするコンピュータシステム。
【請求項15】
前記システムは、前記修正された一連の状態に従って前記粒子が後続のメッシュ位置に移流した後に、前記粒子状態に局所圧力項を加え直してからモーメントを計算して正しい圧力勾配（∇ｐ）をもたらすための命令をさらに含む、請求項１４に記載のシステム。
【請求項16】
前記局所圧力項は、前記システムによって計算されるθ－１項を含む、請求項１５に記載のシステム。
【請求項17】
前記流体流の活動をシミュレートするための命令は、後続のメッシュ位置への前記粒子の移流を実行し、
第１の一連の離散格子速度に部分的に基づいて前記流体流をシミュレートし、及び
前記第１の一連の離散格子速度と同じ格子速度である、或いは前記第１の一連の離散格子速度とは格子速度の数が異なる第２の一連の離散格子速度に部分的に基づいて、スカラー量の時間発展をシミュレートする、ための命令を含む、
請求項１５に記載のシステム。
【請求項18】
前記システムに、平衡分布及び比全エネルギーＥtの第２の分布関数を適用させるための命令をさらに含み、前記第２の分布関数は、流れ分布ｆiと共に移流する固有のスカラーとして定義される、請求項１４に記載のシステム。
【発明の詳細な説明】
【技術分野】
【0001】
〔優先権の主張〕
本出願は、２０１９年１月１０日に出願された「全エネルギー保存を実施する格子ボルツマンソルバ（ＬＡＴＴＩＣＥ ＢＯＬＴＺＭＡＮＮ ＳＯＬＶＥＲ ＥＮＦＯＲＣＩＮＧ ＴＯＴＡＬ ＥＮＥＲＧＹ ＣＯＮＳＥＲＶＡＴＩＯＮ）」という名称の米国仮特許出願第６２／７９０，５２８号に対する合衆国法典第３５編第１１９条（ｅ）に基づく優先権を主張するものであり、この文献の内容は全体が引用により本明細書に組み入れられる。
本出願は、２０１９年１月１０日に出願された「全エネルギー保存を実施する格子ボルツマンソルバ（ＬＡＴＴＩＣＥ ＢＯＬＴＺＭＡＮＮ ＳＯＬＶＥＲ ＥＮＦＯＲＣＩＮＧ ＴＯＴＡＬ ＥＮＥＲＧＹ ＣＯＮＳＥＲＶＡＴＩＯＮ）」という名称の米国仮特許出願第６２／７９０，５２８号に対する合衆国法典第３５編第１１９条（ｅ）に基づく優先権を主張するものであり、この文献の内容は全体が引用により本明細書に組み入れられる。
【0002】
本明細書は、物理的流体流（physical fluid flows）などの物理的プロセスのコンピュータシミュレーションに関する。
【背景技術】
【0003】
いわゆる「格子ボルツマン法」（ＬＢＭ）は、計算流体力学において使用される有利な技術である。格子ボルツマンシステムの基礎を成す力学は、ボルツマン方程式に従う多くの粒子の動きを伴う気体分子運動論の基本物理学に存在する。ボルツマン基本運動システムには、衝突及び移流という２つの基本的な動力学プロセスがある。衝突プロセスは、保存則に従う粒子間の相互作用、及び平衡状態への緩和を伴う。移流プロセスは、微視的粒子速度に従う１つの位置から別の位置への粒子の動きをモデル化することを含む。
【0004】
圧縮性流では、質量、運動量及び全エネルギーの保存を満たすために、通常は単一の分布関数を使用する。単一の分布関数を使用してこれらの３つの量を全て解く場合、ＬＢＭに使用される格子は、最大８次のモーメントを満たすべきである。この状況では、モーメント要件が増加し、結果として得られるステンシルサイズも増大するので、数多くの格子速度が必要である。単一の粒子分布関数の使用に伴う他の留意事項としては、分布関数の安定性、運動量及びエネルギーの境界条件、体積力（ｂｏｄｙ ｆｏｒｃｅ）及び熱源などに関する留意事項が挙げられる。
【先行技術文献】
【特許文献】
【0005】
【文献】米国特許出願公開第２０１６－０１８８７６８号明細書
【文献】米国特許第９，５７６，０８７号明細書
【発明の概要】
【課題を解決するための手段】
【0006】
１つの態様によれば、全エネルギー保存を実施するシミュレーションを使用してコンピュータ上で流体流をシミュレートする方法が、メッシュにわたる粒子の動きをモデル化するようにメッシュにわたる流体の活動をシミュレートするステップと、メッシュ内の各メッシュ位置の、それぞれが対応するメッシュ位置における可能な運動量状態のうちの特定の運動量状態に対応する複数のエントリを含む一連の状態ベクトルをコンピュータアクセス可能メモリに記憶するステップと、コンピュータが、メッシュ内のメッシュ位置の一連のエネルギー値を計算するステップと、コンピュータが、粒子の後続のメッシュ位置への移流を一定の時間間隔にわたって実行するステップと、コンピュータが、移流した粒子の状態に比全エネルギーの値（ｓｐｅｃｉｆｉｃ ｔｏｔａｌ ｅｎｅｒｇｙ ｖａｌｕｅｓ）を加算し、時間間隔にわたって移流しなかった粒子の状態から比全エネルギーの値を減算することによって、粒子の状態ベクトルを修正するステップと、を含む。
【0007】
態様は、１又は２以上の機械可読ハードウェア記憶装置上のコンピュータプログラム製品、装置及びコンピュータシステムも含む。
【0008】
上記態様のうちの１つ又は２つ以上は、本明細書で説明する特徴の中でもとりわけ以下の特徴のうちの１つ又は２つ以上を含むことができる。
【0009】
修正された一連の状態に従って粒子が後続のメッシュ位置に移流した後に、態様は、粒子状態に局所圧力項を加え直してからモーメントを計算して正しい圧力勾配（∇ｐ）をもたらすステップを含む。局所圧力項は、コンピュータによって計算されるθ－１項を含む。流体流の活動をシミュレートするステップは、第１の一連の離散格子速度に部分的に基づいて流体流をシミュレートするステップを含み、方法は、第２の一連の離散格子速度に部分的に基づいてスカラー量の時間発展をシミュレートするステップをさらに含む。第２の一連の離散格子速度は、第１の一連の離散格子速度と同じ格子速度である。第２の一連の離散格子速度は、格子速度の数が第１の一連の離散格子速度の数と異なる。全エネルギーの時間発展をシミュレートするステップは、近隣メッシュ位置からの流入分布を衝突演算子及びエネルギー演算子について収集するステップと、流入分布に重み付けするステップと、衝突演算子とエネルギー演算子との積として流出分布を決定するステップと、決定された流出分布を伝播するステップと、を含む。
【0010】
スカラー量の時間発展をシミュレートするステップは、１次項としての衝突演算子とエネルギー演算子との積に由来する運動量流束を効果的に表す衝突演算子に部分的に基づいてスカラー量の時間発展をシミュレートするステップを含む。任意のプラントル数のための正確な剪断応力を回復させるために、エネルギー衝突項は以下によって与えられ、
ここでの
は、非平衡成分のフィルタ処理された２次モーメントである。
【数1】
ここでの
【数2】
は、非平衡成分のフィルタ処理された２次モーメントである。
【0011】
態様は、流入分布が、決定された流出分布に等しくなるように、ゼロの正味表面流束境界条件を適用するステップをさらに含む。流出分布を決定するステップは、ゼロの表面スカラー流束をもたらすように流出分布を決定するステップを含む。スカラー量は、温度、濃度及び密度から成る群から選択されたスカラー量を含む。メッシュ内のメッシュ位置の一連のエネルギー値を計算するステップは、Ｅt＝Ｅ＋ν2／２によって与えられる全エネルギーＥtを計算するステップをさらに含み、ここでのＥはエネルギーであり、νは速度である。
【0012】
活動をシミュレートするステップは、コンピュータが、平衡分布及び比全エネルギーＥtの第２の分布関数を適用するステップを含み、第２の分布は、流れ分布ｆiと共に移流する特定のスカラーとして定義される。第２の分布関数は、エネルギー方程式に対するｆiの非平衡寄与を考慮して、境界付近の異なる格子分解能にわたる正確な流れ挙動を取得し、分布関数の衝突演算子は以下によって与えられ、
ここでの項Ωf、Ωqは、それぞれの衝突演算子を表し、上記方程式において使用される平衡分布は以下の通りである、
【数3】
【数4】
ここでの項Ωf、Ωqは、それぞれの衝突演算子を表し、上記方程式において使用される平衡分布は以下の通りである、
【数5】
【数6】
【0013】
全エネルギーは、圧力が対流するように移流前に項を追加し、追加エネルギーを補正するために停止状態からエネルギーを削除する。追加エネルギーを削除することによって全エネルギーが保存され、正しい圧力速度項がもたらされる。追加エネルギーを削除することは、停止状態
を
となるように修正することによって行われる。
【数7】
を
【数8】
となるように修正することによって行われる。
【0014】
本明細書に開示する方法では、全エネルギーを解くことが、流れ分布と共に移流する特定のスカラーとして定義される第２の分布を使用する。本明細書で説明する流体流をシミュレートする方法は、格子ボルツマン（ＬＢ）法及びスカラーソルバ輸送方程式を使用する。全エネルギー保存を有して高安定域を伴う圧縮性流れのためのＬＢＭソルバを提供し、これを実際の応用に使用することができる。衝突中には、ソルバにおいて使用される平衡分布が正であり、移流中には、温度結合が行われる。
【0015】
図面を含む以下の説明及び特許請求の範囲から、他の特徴及び利点が明らかになるであろう。
【図面の簡単な説明】
【0016】
【図1】流体流シミュレーションシステムを示す図である。
【図2】複雑な流体流シミュレーション及び全エネルギーの計算のための動作を示すフローチャートである。
【図3】複雑な流体流シミュレーション及び全エネルギーの計算のための動作を示すフローチャートである。
【図4】期待結果、予測結果の例示的な計算ベースの比較を示す図である。
【図5】期待結果、予測結果の例示的な計算ベースの比較を示す図である。
【図6】期待結果、予測結果の例示的な計算ベースの比較を示す図である。
【図7】２つのＬＢＭモデルの速度成分を示す図である（先行技術）。
【図8】２つのＬＢＭモデルの速度成分を示す図である（先行技術）。
【図9】物理的プロセスシミュレーションシステムが従う手順のフローチャートである（先行技術）。
【図10】マイクロブロックの透視図である。
【図12】可変分解能技術を示す図である（先行技術）。
【図13】可変分解能技術を示す図である（先行技術）。
【図14】粒子とファセットとの相互作用を理解する上で役立つ平行六面体を示す図である。
【図15】ボクセルから表面への粒子の動きを示す図である（先行技術）。
【図16】表面から表面への粒子の動きを示す図である（先行技術）。
【図17】表面動力学の実行手順のフローチャートである（先行技術）。
【図18】粒子衝突ステップから独立した熱力学ステップを使用して分布関数を生成するプロセスのフローチャートである。
【発明を実施するための形態】
【0017】
以下の図９で説明する手順では、全エネルギー保存を実施する格子ボルツマンソルバを用いたフローシミュレーションプロセスについて説明する。
【0018】
図７、図８及び図１０～図１８では、図にそれぞれ「先行技術」との表記を行っている。これらの図に先行技術との表記を行う理由は、一般にこれらの図が上述した特許で見られるためである。しかしながら、これらの図は、上述した特許出願で見られる通りである。しかしながら、上述した特許出願には、後述する全エネルギー保存を実施する格子ボルツマンソルバについての記載がないため、これらの図は、このＬＢソルバを使用してフローシミュレーションに対して行う修正を一切考慮していない。
【0019】
Ａ．スカラー量を解くための方法
複雑な流体流シミュレーションを達成する際には、流体流を解くのと共に、温度分布、濃度分布及び／又は密度などのスカラー量を同時に解くことが有益となり得る。本明細書で説明するシステム及び方法では、ＬＢＭベースの物理プロセスシミュレーションシステムに基づく流体流のモデル化に（ベクトル量とは対照的に）スカラー量のモデル化を結びつける。モデル化できる例示的なスカラー量としては、温度、濃度及び密度が挙げられる。
複雑な流体流シミュレーションを達成する際には、流体流を解くのと共に、温度分布、濃度分布及び／又は密度などのスカラー量を同時に解くことが有益となり得る。本明細書で説明するシステム及び方法では、ＬＢＭベースの物理プロセスシミュレーションシステムに基づく流体流のモデル化に（ベクトル量とは対照的に）スカラー量のモデル化を結びつける。モデル化できる例示的なスカラー量としては、温度、濃度及び密度が挙げられる。
【0020】
スカラー量のモデル化をシミュレーション流体流に結びつける技術の具体的詳細は、Ｐｒａｄｅｅｐ Ｇｏｐａｌａｋｒｉｓｈｍａｎ他による「熱的混成格子ボルツマン法のための温度結合アルゴリズム（Ｔｅｍｐｅｒａｔｕｒｅ Ｃｏｕｐｌｉｎｇ Ａｌｇｏｒｉｔｈｍ ｆｏｒ Ｈｙｂｒｉｄ Ｔｈｅｒｍａｌ Ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ Ｍｅｔｈｏｄ）」という名称の同時係属中の米国特許出願公開第２０１６－０１８８７６８号に詳述されており、この文献の内容は全体が引用により本明細書に組み入れられる。比較的高い安定域を有するガリレイ不変の衝突演算子（Ｇａｌｉｌｅａｎ－ｉｎｖａｒｉａｎｔ ｃｏｌｌｉｓｉｏｎ ｏｐｅｒａｔｏｒ）に基づく衝突演算子を使用する。この衝突演算子の具体的詳細は米国特許第９，５７６，０８７号に記載されており、この文献の内容は全体が引用により本明細書に組み入れられる。
【0021】
例えば、このシステムを使用して、システム内の対流温度分布を求めることができる。例えば、（複数のボクセルによって表される体積で形成された）システムが熱源を含み、システム内に空気流が存在する場合には、空気流及び熱源との近接性に基づいて、システムの一部の領域が他の領域よりも高温になる。このような状況をモデル化するために、システム内の温度分布を、各ボクセルが関連する温度を有するスカラー量として表すことができる。
【0022】
別の例では、このシステムを使用して、システム内の化学物質分布を求めることができる。例えば、（複数のボクセルによって表される体積で形成された）システムが、汚染爆弾、或いは空気中又は液体中に懸濁した化学物質又は他の粒子などの汚染源を含み、システム内に空気流又は液体流が存在する場合には、流れ及び汚染源との近接性に基づいて、システムの一部の領域が他の領域よりも高濃度になる。このような状況をモデル化するために、システム内の化学物質分布を、各ボクセルが関連する濃度を有するスカラー量として表すことができる。
【0023】
いくつかの応用では、複数の異なるスカラー量を同時にシミュレートすることができる。例えば、システムは、システム内の温度分布及び濃度分布の両方をシミュレートすることができる。
【0024】
スカラー量は、様々な方法でモデル化することができる。例えば、スカラー輸送方程式を解くための格子ボルツマン（ＬＢ）法を使用して間接的にスカラー輸送を解くことができる。例えば、本明細書で説明する方法は、以下の２次巨視的スカラー輸送方程式の間接解をもたらすことができる。
【数9】
【0025】
流体流のための格子ボルツマン関数に加え、輸送スカラーに第２の分布関数の組を導入する。この方法は、体積内の各ボクセルにベクトルを割り当てて流体流を表し、体積内の各ボクセルにスカラー量を割り当てて所望のスカラー変数（例えば、温度、密度、濃度など）を表す。この方法は、正確な保存則を満たす巨視的スカラー輸送方程式を完全に回復させる。この方法は、他の非ＬＢＭ法と比べて、求められるスカラー量の精度を高めると考えられる。また、この方法は、複雑な境界形状を考慮する能力を高めるとも考えられる。
【0026】
このスカラー量をモデル化する方法は、マサチューセッツ州バーリントンのＥｘａ Ｃｏｒｐｏｒａｔｉｏｎ社から市販されているＰｏｗｅｒＦＬＯＷシステムなどの、格子ボルツマン法（ＬＢＭ）に基づく時間陽的ＣＦＤ／ＣＡＡ解法と共に使用することができる。ＬＢＭは、巨視的連続体方程式を離散化することに基づく方法とは異なり、「メゾスコピック」ボルツマン運動方程式から開始して巨視的流体力学を予測する。結果として得られる圧縮非定常解法（ｃｏｍｐｒｅｓｓｉｂｌｅ ａｎｄ ｕｎｓｔｅａｄｙ ｓｏｌｕｔｉｏｎ ｍｅｔｈｏｄ）を使用して、空力音響問題及び純音響問題などの様々な複雑な流れの物理特性を予測することができる。以下では、ＬＢＭベースのシミュレーションシステムの一般的考察を示した後に、このようなモデル化法をサポートするために流体流シミュレーションと共に使用できるスカラー解法について説明する。
【0027】
シミュレーション空間のモデル化
ＬＢＭベースの物理的プロセスシミュレーションシステムでは、流体流が、一連の離散速度ｃiで評価される分布関数値ｆiによって表される。分布関数の動力学は方程式Ｉ１に準拠し、ここでのｆi（０）は、以下のように定義される平衡分布関数として知られており、
方程式（Ｉ１）
式中、
である。
ＬＢＭベースの物理的プロセスシミュレーションシステムでは、流体流が、一連の離散速度ｃiで評価される分布関数値ｆiによって表される。分布関数の動力学は方程式Ｉ１に準拠し、ここでのｆi（０）は、以下のように定義される平衡分布関数として知られており、
【数10】
方程式（Ｉ１）
式中、
【数11】
である。
【0028】
この方程式は、分布関数ｆiの時間発展を表す周知の格子ボルツマン方程式である。左辺は、いわゆる「流動過程」による分布の変化を表す。流動過程は、あるメッシュ位置において流体ポケットが開始した後に、複数の速度ベクトルのうちの１つに沿って次のメッシュ位置に移動する時点である。この時点で、「衝突係数」、すなわち開始流体ポケットに対する近隣の流体ポケットの影響を計算する。流体は別のメッシュ位置にしか移動することできず、従って全ての速度の全ての成分が共通速度の倍数になるように正しい速度ベクトルの選択が必要である。
【0029】
第１の方程式の右辺は、流体ポケット同士の衝突に起因する分布関数の変化を表す上述した「衝突演算子」である。特定の形の衝突演算子は、Ｂｈａｔｎａｇａｒ、Ｇｒｏｓｓ ａｎｄ Ｋｒｏｃｋ（ＢＧＫ）演算子とすることができる。衝突演算子は、「平衡」形態である第２の方程式によって与えられる規定値に分布関数を至らしめる。
【0030】
ＢＧＫ演算子は、たとえ衝突の詳細がどうであろうと、分布関数は、衝突を通じて｛ｆi
eq（ｘ，ｖ，ｔ）｝｛ｆeq（ｘ，ｖ，ｔ）｝によって与えられる明確に定義された局所平衡に近付くという物理的説明に従って構築され、
方程式(Ｉ２)
ここでのパラメータτは、衝突を介した平衡までの特性緩和時間を表す。粒子（例えば、原子又は分子）の取り扱いでは、通常は緩和時間が定数とみなされる。
【数12】
方程式(Ｉ２)
ここでのパラメータτは、衝突を介した平衡までの特性緩和時間を表す。粒子（例えば、原子又は分子）の取り扱いでは、通常は緩和時間が定数とみなされる。
【0031】
このシミュレーションから、質量ρ及び流体速度ｕなどの従来の流体変数が以下の方程式（Ｉ３）の単純和として得られる。
方程式(Ｉ３)
ここでのρ、ｕ及びＴは、それぞれ流体の密度、速度及び温度であり、Ｄは、（必ずしも物理的空間次元とは等しくない）離散化速度空間の次元である。
【数13】
方程式(Ｉ３)
ここでのρ、ｕ及びＴは、それぞれ流体の密度、速度及び温度であり、Ｄは、（必ずしも物理的空間次元とは等しくない）離散化速度空間の次元である。
【0032】
対称性を考慮することにより、速度値の集合は、配置空間内でスパニングした時に特定の格子構造を形成するように選択される。このような離散系の動特性は、以下の形を有するＬＢＥに従い、
ｆi（ｘ＋ｃi，ｔ＋１）－ｆi（ｘ，ｔ）＝ｃi（ｘ，ｔ）
ここでの衝突演算子は、通常は上述したようなＢＧＫの形を取る。平衡分布形態を正しく選択することにより、格子ボルツマン方程式が正しい流体力学及び熱流体力学の結果をもたらすことを理論的に示すことができる。すなわち、ｆi（ｘ，ｔ）から導き出される流体力学的モーメントは、巨視的極限ではナビエ・ストークス方程式に従う。これらのモーメントは、上記の方程式（Ｉ３）によって定義される。
ｆi（ｘ＋ｃi，ｔ＋１）－ｆi（ｘ，ｔ）＝ｃi（ｘ，ｔ）
ここでの衝突演算子は、通常は上述したようなＢＧＫの形を取る。平衡分布形態を正しく選択することにより、格子ボルツマン方程式が正しい流体力学及び熱流体力学の結果をもたらすことを理論的に示すことができる。すなわち、ｆi（ｘ，ｔ）から導き出される流体力学的モーメントは、巨視的極限ではナビエ・ストークス方程式に従う。これらのモーメントは、上記の方程式（Ｉ３）によって定義される。
【0033】
ｃi及びｗiの集合的な値はＬＢＭモデルを規定する。ＬＢＭモデルは、スケーラブルコンピュータプラットフォーム上に効率的に実装して、時間非定常流（ｔｉｍｅ ｕｎｓｔｅａｄｙ ｆｌｏｗｓ）及び複雑な境界条件のために高いロバスト性を伴って実行することができる。
【0034】
ボルツマン方程式から流体系の動きの巨視的方程式を取得する標準技術は、完全なボルツマン式の連続近似を取るチャップマン－エンスコッグ法（Ｃｈａｐｍａｎ－Ｅｎｓｋｏｇ ｍｅｔｈｏｄ）である。流体系では、密度のわずかな乱れが音速で伝わる。気体系では、一般に音速が温度によって決まる。流れの圧縮率の影響の重要性は、特性速度と音速との比率によって測定され、これはマッハ数として知られている。従来のＬＢＭベースの物理プロセスシミュレーションシステムのさらなる説明については、上記で引用により組み入れた米国特許出願公開第２０１６／０１８８７６８号を参照されたい。
【0035】
図１に、物理的対象表現の周囲などの流体流をシミュレートするシステム１０を示す。この実装のシステム１０は、クライアントサーバアーキテクチャに基づき、超並列コンピュータシステム１２として実装されるサーバシステム１２と、クライアントシステム１４とを含む。サーバシステム１２は、メモリ１８と、バスシステム１１と、インターフェイス２０（例えば、ユーザインターフェイス／ネットワークインターフェイス／ディスプレイ又はモニタインターフェイスなど）と、処理装置２４とを含む。メモリ１８内には、メッシュ準備エンジン３２及びシミュレーションエンジン３４が存在する。システム１０は、２Ｄ及び／又は３Ｄメッシュ、座標系及びライブラリを記憶するデータリポジトリ３８にアクセスする。
【0036】
図１では、メッシュ準備エンジン３２をメモリ１８内に示しているが、メッシュ準備エンジンは、サーバ１２とは異なるシステム上で実行されるサードパーティアプリケーションとすることもできる。メッシュ準備エンジン３２は、メモリ１８内で実行されるか、それともサーバ１２とは異なるシステム上で実行されるかにかかわらず、ユーザ指定のメッシュ定義３０を受け取り、メッシュを準備してシミュレーションエンジン３４に送信する。以下で開示するように、シミュレーションエンジン３４は、粒子衝突相互作用モジュールと、粒子境界モデルモジュールと、移流動作を実行する移流モジュールとを含む。シミュレーションエンジン３４は、図３で後述するような全エネルギーソルバも含む。
【0037】
次に、図２に、物理的対象表現の周囲の流体流をシミュレートするプロセス４０を示す。本明細書で説明する例では、物理的対象が翼形部である。しかしながら、物理的対象はあらゆる形状とすることができ、具体的には平面及び／又は（単複の）曲面を有するので、翼形部の使用は例示にすぎない。プロセスは、例えばクライアントシステム１４から、又はデータリポジトリ４２から検索を行うことにより、例えばシミュレートする物理的対象のメッシュを受け取る（４２）。他の実施形態では、外部システム又はサーバ１２のいずれかが、シミュレートする物理的対象のメッシュをユーザ入力に基づいて生成する。このプロセスは、検索されたメッシュから幾何学量を事前計算し（４４）、検索されたメッシュに対応する事前計算された幾何学量を用いて動的格子ボルツマンモデルシミュレーションを実行する（４６）。格子ボルツマンモデルシミュレーションは、移流した粒子の状態に比全エネルギーの値を加算し、一定の時間間隔にわたって移流しなかった粒子の状態から比全エネルギーの値を減算することによって修正された状態ベクトルの組に従う、粒子分布の発展と、ＬＢＭメッシュ内の次のセル
への粒子の移流とのシミュレーションを含む。
【数14】
への粒子の移流とのシミュレーションを含む。
【0038】
全エネルギー保存圧縮性流れソルバ
複雑な流体流シミュレーションを達成する際には、流体流を解くのと共に、温度分布、濃度分布及び／又は密度などのスカラー量を同時に解くことが有益となり得る。特定のスカラーとスカラー分布とを区別することは、ｑ（方程式１０）とｈ（方程式２）との間の差分であり、Ｅ＿ｔである特定のスカラーｑを使用するのに対し、ｈは密度＊Ｅ＿ｔである。説明する方法は、上述した米国特許出願公開第２０１６－０１８８７６８号において説明されているように、移流する流れの上に存在するスカラーのようなものである。
複雑な流体流シミュレーションを達成する際には、流体流を解くのと共に、温度分布、濃度分布及び／又は密度などのスカラー量を同時に解くことが有益となり得る。特定のスカラーとスカラー分布とを区別することは、ｑ（方程式１０）とｈ（方程式２）との間の差分であり、Ｅ＿ｔである特定のスカラーｑを使用するのに対し、ｈは密度＊Ｅ＿ｔである。説明する方法は、上述した米国特許出願公開第２０１６－０１８８７６８号において説明されているように、移流する流れの上に存在するスカラーのようなものである。
【0039】
全エネルギー保存における圧力項の処理
次に、図３を参照して、従来の方法よりも良好な安定域を有して実際の用途に応用できる全エネルギー保存プロセス５０を用いた圧縮性ＬＢＭソルバについて説明する。解く必要がある全エネルギー方程式は、以下によって与えられ、
方程式１
ここでの全エネルギーＥtは、Ｅt＝Ｅ＋ν2／２によって与えられ、ρは密度であり、νは速度であり、ｐは圧力であり、ｋは熱伝導率であり、Ｔは温度であり、τは剪断応力である。この方程式は、高速流のための全計算流体力学ツールによって解く必要がある一般方程式である。
次に、図３を参照して、従来の方法よりも良好な安定域を有して実際の用途に応用できる全エネルギー保存プロセス５０を用いた圧縮性ＬＢＭソルバについて説明する。解く必要がある全エネルギー方程式は、以下によって与えられ、
【数15】
方程式１
ここでの全エネルギーＥtは、Ｅt＝Ｅ＋ν2／２によって与えられ、ρは密度であり、νは速度であり、ｐは圧力であり、ｋは熱伝導率であり、Ｔは温度であり、τは剪断応力である。この方程式は、高速流のための全計算流体力学ツールによって解く必要がある一般方程式である。
【0040】
全エネルギーソルバに伴う複雑さの１つは、方程式１における圧力対流項「ｐ」の加算である。格子ボルツマン法（ＬＢＭ）は、高速流のための全計算流体力学ツールと同様に、この方程式を全エネルギー保存について解く必要がある。ＬＢＭ法には、全ての事例について、とりわけ圧力対流項ｐについてこの方程式を解く上でいくつかの欠点がある。開示する方法は、この方程式を複雑な問題について解くことができ、他の利点も有する。
【0041】
ＬＢＭ法について存在するような圧力対流問題は、以下のように説明することができる。正常なスカラー輸送は以下の方程式の形を取る。
上記方程式の左辺は対流項を含み、右辺は拡散項を有する。全エネルギー方程式とは異なり、上記方程式の全ての項は１つのスカラー変数Ｓしか含んでおらず、このことが移流を単純化する。プロセス５０は、選択された平衡分布形態に従って粒子分布の発展をシミュレートする（５２）。
【数16】
上記方程式の左辺は対流項を含み、右辺は拡散項を有する。全エネルギー方程式とは異なり、上記方程式の全ての項は１つのスカラー変数Ｓしか含んでおらず、このことが移流を単純化する。プロセス５０は、選択された平衡分布形態に従って粒子分布の発展をシミュレートする（５２）。
【0042】
平衡分布（方程式１の選択形態）５２に加え、比全エネルギーＥtを計算するための第２の分布関数を使用する（５４）。プロセス５０は、ＬＢＭにおける次のセルｑへの粒子の移流（５６）を含む。エネルギー保存の実現は、２つの分布の積によって行われる。これらの分布関数の積は、流れ分布ｆiと共に移流する特定のスカラーとして定義される。この第２の分布関数は、エネルギーを解くために単独で第２の分布関数を使用する既存の方法とは異なり、エネルギー方程式に対するｆiの非平衡寄与を考慮して、境界付近の異なる格子分解能にわたる正確な流れ挙動を取得する。
【0043】
分布関数の衝突演算子は以下によって与えられる。
方程式２
方程式３
【数17】
方程式２
【数18】
方程式３
【0044】
方程式２及び方程式３の項Ωf及びΩqは、それぞれの衝突演算子を表すために使用される。方程式５の衝突演算子Ωfについては、米国特許第９，５７６，０８７号において広く説明されており、この文献の内容は全体が引用により本明細書に組み入れられる。上記の方程式で使用される平衡分布は以下の通りである。
方程式４
【数19】
【数20】
方程式４
【0045】
全エネルギーソルバのための分布関数方程式７の上記公式化において明白な観察結果は、平衡分布が、方程式１などにおける圧力項（ｐ）、及びθ－１項、熱流に不可欠なｐ＝ρＲＴを含まない点である。第２の分布は、内部エネルギーＥ＝ＣνＴ及び運動エネルギーν2／２の総和である比全エネルギーである。
【0046】
従って、この分布関数は、衝突演算中に常に正の値をもたらし、従ってこの分布関数は、対応する高い安定性をソルバにもたらす。方程式（２）における衝突後の状態がそのまま移流した場合、方程式４の分布関数は、運動量の等温圧力勾配のみをもたらし、全エネルギー方程式に圧力対流項をもたらさない。平衡分布を使用して移流中及び衝突中に以下のような圧力項を含む他の方法とは対照的に、移流に起因して圧力勾配項及び圧力対流項が生じる一方で、衝突はこれらの項に影響を与えない。
【数21】
【数22】
【0047】
項θ＝ＲＴ／Ｔ0は、圧力に対応して上記の平衡を負の状態にするものであり、従って既存の技術は非常に不安定である。また、第２の分布関数ｈi
eqを全エネルギーについて単独で使用していることにより、異なる格子分解能にわたってノイズも生じる。
【0048】
処理５０は、移流した状態に比全エネルギーを追加する一方で停止状態（非移動状態）から比全エネルギーを削除することによって移流状態を修正する（５８）。
【0049】
具体的には、この正確な運動量及びエネルギー方程式を回復させる修正（５８）では、圧力の導入が移流ステップ中にしか行われない。移流プロセスの前には、メッシュサイトにおける移動状態（∀ｉ≠０）値が、以下の方程式８のように修正される。
方程式５
【数23】
【数24】
方程式５
【0050】
しかしながら、この項（方程式５の第２の方程式）をｑiに加えると、移動状態に追加エネルギーも加わる。この項によって表される追加エネルギーを補正するために、停止状態（粒子が対流しない状態）から同じ量のエネルギーを削除する。プロセスは、これを行うことによってエネルギーを保存し、正確な圧力速度項をもたらす。具体的にいえば、全エネルギーの保存は、停止状態
を
となるように修正することによって行われる。
【数25】
を
【数26】
となるように修正することによって行われる。
【0051】
移流後には、モーメントを計算する前に粒子状態に局所圧力項を加え直し、これによって正しい圧力勾配（∇ｐ）が得られる。
方程式６
【数27】
方程式６
【0052】
移流前（方程式８）及び移流後（方程式６）の修正は、運動量の状態の正しい方程式（又は正しい圧力勾配（∇ｐ））を回復させるとともに、全エネルギーの正しい圧力対流項∇（ｕｐ）も回復させる。圧力項が削除されるので（方程式９）、衝突中に状態が正になって全エネルギーソルバの安定域を改善する。
【0053】
以下で図９において考察するように、移流した粒子の状態に比全エネルギーの値を加算し、その時間間隔にわたって移流しなかった粒子の状態から比全エネルギーの値を減算することによってコンピュータが粒子の状態ベクトルを修正するという移流に対する全エネルギー修正を適用して移流を修正する。
【0054】
エネルギーのための正則化された衝突演算子
流動状態のために使用される衝突演算子Ωfは、上述した米国特許第９，５７６，０８７号において検討されている。以下では、エネルギーのために使用される衝突演算子Ωqについて検討する。１次モーメントを伴う正則化された衝突は以下のように与えられ、
方程式７
ここでの方程式７は、全エネルギーに対するフィルタ処理された非平衡寄与である。エネルギー保存は、方程式（７）の右辺における第２項
である衝突項を以下の状態にする２つの状態の乗算によって満たされる。
方程式８
流動状態のために使用される衝突演算子Ωfは、上述した米国特許第９，５７６，０８７号において検討されている。以下では、エネルギーのために使用される衝突演算子Ωqについて検討する。１次モーメントを伴う正則化された衝突は以下のように与えられ、
【数28】
【数29】
方程式７
ここでの方程式７は、全エネルギーに対するフィルタ処理された非平衡寄与である。エネルギー保存は、方程式（７）の右辺における第２項
【数30】
である衝突項を以下の状態にする２つの状態の乗算によって満たされる。
【数31】
方程式８
【0055】
上記の項のゼロ番目のモーメントは０であり、エネルギーを保存する。チャップマン・エンスコッグの拡張中には、衝突項の１次モーメントが熱拡散の導出に関与し、剪断応力によって機能する。衝突項は、１次近似として以下のように表すことができる。
方程式９
【数32】
方程式９
【0056】
上記の方程式では、両状態
を乗算した結果、運動量流束項ρυυが生じる。この項ρυυは、高速流では数値的不安定性をもたらすとともに、マッハ数依存の拡散ももたらす。これらの影響を軽減するために、衝突演算子を以下のように定義することによって項ρυυを削除する。
方程式１０
【数33】
を乗算した結果、運動量流束項ρυυが生じる。この項ρυυは、高速流では数値的不安定性をもたらすとともに、マッハ数依存の拡散ももたらす。これらの影響を軽減するために、衝突演算子を以下のように定義することによって項ρυυを削除する。
【数34】
方程式１０
【0057】
上記の式は、運動量流束項ρυυを削除して、高マッハ流での数値的安定性を改善する。
【0058】
衝突演算子には多くの方法を利用することができる。これらの方法の１つは「正則化」と呼ばれる。Ｈｕｄｏｎｇ及びＲａｏｙａｎｇによる（係属中の）格子ボルツマンに基づくスカラーについての特許では、「正則化」がさらなるガリレイ不変性特徴（Ｇａｌｉｌｅａｎ ｉｎｖａｒｉａｎｃｅ ｆｅａｔｕｒｅ）と併せて使用されている。しかしながら、この正則化衝突演算子は、高速流との使用にとって理想的なものではない。導出中には、結局のところ演算子がρυυなどの高次項を有するようになる（方程式１２）。ＬＢＭでは、一般に速度値が常に１未満であり、従って低速流ではこれらの誤差項がはるかに小さくなり、無視することもできる。しかしながら、高速流では、速度値がほぼ１程度に高くなり、従って誤差項ρυυが比較的大きくなり、方程式１０によって削除される。
【0059】
任意のプラントル数のための粘性効果
エネルギーソルバで使用される緩和時間は、熱拡散率に等しい粘度を有する剪断応力（τq－０．５ρＲＴによって誤った結果をもたらすτqである。任意のプラントル数のための正確な剪断応力を回復させるために、エネルギー衝突項に以下のさらなる項を追加し、
方程式１１
ここでの
は、非平衡成分のフィルタ処理された２次モーメントである。
エネルギーソルバで使用される緩和時間は、熱拡散率に等しい粘度を有する剪断応力（τq－０．５ρＲＴによって誤った結果をもたらすτqである。任意のプラントル数のための正確な剪断応力を回復させるために、エネルギー衝突項に以下のさらなる項を追加し、
【数35】
方程式１１
ここでの
【数36】
は、非平衡成分のフィルタ処理された２次モーメントである。
【0060】
全エネルギーソルバの能力
全エネルギーソルバは、ＰｏｗｅｒＦＬＯＷ（登録商標）などの既存の流れソルバに含めることができ、亜音速流及び遷音速流に関連する応用を含む幅広い産業用途に使用することができる。
全エネルギーソルバは、ＰｏｗｅｒＦＬＯＷ（登録商標）などの既存の流れソルバに含めることができ、亜音速流及び遷音速流に関連する応用を含む幅広い産業用途に使用することができる。
【0061】
本節では、開示する全エネルギーソルバの、とりわけエネルギー保存面に関連する圧縮性流についての潜在的利点の一部を示すために、いくつかのベンチマークについて説明する。シミュレーション結果を、解析解及び典型的な有限差分ＰＤＥベースのハイブリッドソルバ解と比較する。
【0062】
図４に、前縁からｘ＝１．５ｍの平板にわたる、Ｍａ＝０．７のマッハ速度（Ｍａ）における流れの全温度のグラフを示す。図４に示すように、平板にわたる流れの全エネルギーソルバ（実線６２）は全温度を正確に回復させるのに対し、ハイブリッドソルバ（ＬＢＭ流れソルバ及び有限差分ベースのエントロピーソルバ）（点線６４）は全温度を維持していない。
【0063】
図５には、前縁からｘ＝１．５ｍの平板のＭａ＝０．７のマッハ速度（Ｍａ）における静温度のグラフを示す。このグラフには、ｕ／ｕ∞（ｙ軸）としてのｕに対して温度ＴをＴ／Ｔ∞（ｘ軸）として示す。実線６６として表した開示する全エネルギーソルバをハイブリッドソルバ（ＬＢＭ流れソルバ及び有限差分ベースのエントロピーソルバ）（点線６８）と比較する。この図１には、前縁からｘ＝１．５ｍの平板のＭａ＝０．７のマッハ速度（Ｍａ）における静温度がウォルツの方程式（アスタリスク線７０）に一致することが示されている。
【0064】
図６には、逆圧Ｐb＝０：７５Ｐtの縮小拡大検証（Ｃｏｎｖｅｒｇｉｎｇ－Ｄｉｖｅｒｇｉｎｇ Ｖｅｒｉｆｉｃａｔｉｏｎ）（ＣＤＶ）ノズルの温度シミュレーション結果のグラフを示す。このグラフには、ｘ（ｙ軸）に対して温度ＴをＴ／Ｔ∞（ｘ軸）として示す。実線７２として表した開示する全エネルギーソルバを、ハイブリッドソルバ（ＬＢＭ流れソルバ及び有限差分ベースのエントロピーソルバ）（点線７４）と比較する。開示する全エネルギーソルバは、分析決定プロットの破線７６と良好に一致する。静温度の比較では、有限差分法に基づく圧縮性ソルバが衝突後に温度を回復できていないのに対し、開示する全エネルギーソルバは期待値を達成することが示されている。
【0065】
開示する全エネルギーソルバを用いたシステムを使用して、システム内の対流温度分布を求めることができる。例えば、（複数のボクセルによって表される体積で形成された）システムが熱源を含み、システム内に空気流が存在する場合には、空気流及び熱源との近接性に基づいて、システムの一部の領域が他の領域よりも高温になる。このようなシステムをモデル化するために、システム内の温度分布を、各ボクセルが関連する温度を有するスカラー量として表すことができる。
【0066】
別の例では、このシステムを使用して、システム内の化学物質分布を求めることができる。例えば、（複数のボクセルによって表される体積で形成された）システムが、汚染爆弾、或いは空気中又は液体中に懸濁した化学物質又は他の粒子などの汚染源を含み、システム内に空気流又は液体流が存在する場合には、流れ及び汚染源との近接性に基づいて、システムの一部の領域が他の領域よりも高濃度になる。このような状況をモデル化するために、システム内の化学物質分布を、各ボクセルが関連する濃度を有するスカラー量として表すことができる。
【0067】
いくつかの応用では、複数の異なるスカラー量を同時にシミュレートすることができる。例えば、システムは、システム内の温度分布及び濃度分布の両方をシミュレートすることができる。
【0068】
スカラー量は、様々な方法でモデル化することができる。例えば、スカラー輸送方程式を解くための格子ボルツマン（ＬＢ）法を使用して間接的にスカラー輸送を解くことができる。例えば、本明細書で説明する方法は、上述した以下の２次巨視的スカラー輸送方程式の間接解をもたらすことができる。
【数37】
【0069】
このような配置シミュレーションでは、流体流のための格子ボルツマン関数に加え、輸送スカラーに第２の分布関数の組が導入される。この方法は、体積内の各ボクセルにベクトルを割り当てて流体流を表し、体積内の各ボクセルにスカラー量を割り当てて所望のスカラー変数（例えば、温度、密度、濃度など）を表す。この方法は、正確な保存則を満たす巨視的スカラー輸送方程式を完全に回復させる。この方法は、他の非ＬＢＭ法と比べて、求められるスカラー量の精度を高めると考えられる。また、この方法は、複雑な境界形状を考慮する能力を高めるとも考えられる。
【0070】
このスカラー量をモデル化する方法は、マサチューセッツ州バーリントンのＥｘａ Ｃｏｒｐｏｒａｔｉｏｎ社から市販されているＰｏｗｅｒＦＬＯＷシステムなどの、格子ボルツマン法（ＬＢＭ）に基づく時間陽的ＣＦＤ／ＣＡＡ解法と共に使用することができる。ＬＢＭは、巨視的連続体方程式を離散化することに基づく方法とは異なり、「メゾスコピック」ボルツマン運動方程式から開始して巨視的流体力学を予測する。結果として得られる圧縮非定常解法を使用して、空力音響問題及び純音響問題などの様々な複雑な流れの物理特性を予測することができる。以下では、ＬＢＭベースのシミュレーションシステムの一般的考察を示した後に、このようなモデル化法をサポートするために流体流シミュレーションと共に使用できるスカラー解法について説明する。
【0071】
図７を参照すると、第１のモデル（２Ｄ－１）１００は、２１個の速度を含む２次元モデルである。これらの２１個の速度のうち、１つ（１０５）は動いていない粒子を表し、３組の４つの速度は、格子のｘ軸又はｙ軸のいずれかに沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）（１１０～１１３）、正規化速度の２倍（２ｒ）（１２０～１２３）、又は正規化速度の３倍（３ｒ）（１３０～１３３）のいずれかで動いている粒子を表し、２組の４つの速度は、ｘ格子軸及びｙ格子軸の両方に関して正規化速度（ｒ）（１４０～１４３）又は正規化速度の２倍（２ｒ）（１５０～１５３）で動いている粒子を表す。
【0072】
また、図８に示すように、第２のモデル（３Ｄ－１）２００は、各速度が図８の矢印のうちの１つによって表される３９個の速度を含む３次元モデルである。これらの３９個の速度のうち、１つは動いていない粒子を表し、３組の６つの速度は、格子のｘ軸、ｙ軸又はｚ軸に沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）のいずれかで動いている粒子を表し、８つの速度は、ｘ、ｙ、ｚ格子軸の３つ全てに関して正規化速度（ｒ）で動いている粒子を表し、１２個の速度は、ｘ、ｙ、ｚ格子軸の２つに関して正規化速度の２倍（２ｒ）で動いている粒子を表す。
【0073】
１０１個の速度を含む３Ｄ－２モデル及び３７個の速度を含む２Ｄ－２モデルなどのさらに複雑なモデルを使用することもできる。
【0074】
１０１個の速度から成る３次元モデル３Ｄ－２では、１つが動いていない粒子（グループ１）を表し、３組の６つの速度が、格子のｘ、ｙ又はｚ軸に沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）のいずれかで動いている粒子（グループ２、４及び７）を表し、３組の８つの速度が、ｘ、ｙ、ｚ格子軸の３つ全てに関して正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）で動いている粒子（グループ３、８及び１０）を表し、１２個の速度が、ｘ、ｙ、ｚ格子軸の２つに関して正規化速度の２倍（２ｒ）で動いている粒子（グループ６）を表し、２４個の速度が、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度（ｒ）及び正規化速度の２倍（２ｒ）で動いていて残りの軸に関しては動いていない粒子（グループ５）を表し、２４個の速度が、ｘ、ｙ、ｚ格子軸のうちの２つに対して正規化速度（ｒ）で動いていて残りの軸に関して正規化速度の３倍（３ｒ）で動いている粒子（グループ９）を表す。
【0075】
３７個の速度から成る２次元モデル２Ｄ－２では、１つが動いていない粒子（グループ１）を表し、３組の４つの速度が、格子のｘ軸又はｙ軸のいずれかに沿って正の方向又は負の方向のいずれかに正規化速度（ｒ）、正規化速度の２倍（２ｒ）、又は正規化速度の３倍（３ｒ）のいずれかで動いている粒子（グループ２、４及び７）を表し、２組の４つの速度が、ｘ及びｙ格子軸の両方に関して正規化速度（ｒ）又は正規化速度の２倍（２ｒ）で動いている粒子を表し、８つの速度が、ｘ及びｙ格子軸の一方に関して正規化速度（ｒ）で動いていて他方の軸に関して正規化速度の２倍（２ｒ）で動いている粒子を表し、８つの速度が、ｘ及びｙ格子軸の一方に関して正規化速度（ｒ）で動いていて他方の軸に対して正規化速度の３倍（３ｒ）で動いている粒子を表す。
【0076】
上述したＬＢＭモデルは、２次元及び３次元の両方における流れの数値シミュレーションのための特定のクラスの効率的でロバストな離散速度動力学モデルを提供する。この種のモデルは、特定の離散速度の組と、これらの速度に関連する重みとを含む。これらの速度は、特に格子ボルツマンモデルとして知られている種類の離散速度モデルの正確で効率的な実装を容易にする速度空間における（非ユークリッド空間内）デカルト座標のグリッド点と一致する。このようなモデルを用いて、高い忠実度で流れをシミュレートすることができる。
【0077】
図９を参照すると、物理プロセスシミュレーションシステムが、手順３００に従って動作して流体流などの物理プロセスをシミュレートする。シミュレーションの前に、シミュレーション空間をボクセルの集合としてモデル化する（ステップ３０２）。通常、シミュレーション空間は、コンピュータ支援設計（ＣＡＤ）プログラムを用いて生成される。例えば、ＣＡＤプログラムを用いて、風洞内に位置するマイクロデバイスを描くことができる。その後、ＣＡＤプログラムによって生成されたデータを処理して適切な分解能の格子構造を追加し、シミュレーション空間内の物体及び表面を考慮する。
【0078】
格子の分解能は、シミュレートされるシステムのレイノルズ数に基づいて選択することができる。レイノルズ数は、流れの粘度（ν）、流れにおける物体の特性長（Ｌ）及び流れの特性速度（ｕ）に関連する。
Ｒｅ＝ｕＬ／ν方程式（Ｉ４）
Ｒｅ＝ｕＬ／ν方程式（Ｉ４）
【0079】
物体の特性長は、物体の大規模な特徴を表す。例えば、マイクロデバイスの周囲の流れをシミュレートする場合には、マイクロデバイスの高さを特性長と見なすことができる。小さな物体領域（例えば、自動車のサイドミラー）の周囲の流れに関心がある時には、シミュレーションの分解能を高めるか、或いは関心領域の周囲に高分解能の領域を使用することができる。格子の分解能が高まるにつれてボクセルの次元は減少する。
【0080】
状態空間はｆi（ｘ，ｔ)として表され、ここでのｆiは、時点ｔにおいて３次元ベクトルｘによって示される、ある格子サイトにおける状態ｉの単位体積当たりの要素又は粒子の数（すなわち、状態ｉの粒子の密度）を表す。既知の時間増分では、粒子の数が単純にｆi（ｘ）として示される。格子サイトの全ての状態の組み合わせは、ｆ（ｘ）として示される。
【0081】
状態の数は、各エネルギーレベル内の考えられる速度ベクトル数によって決まる。速度ベクトルは、ｘ、ｙ及びｚの３次元を有する空間内の整数の線形速度から成る。多重種シミュレーションでは状態の数が増加する。
【0082】
各状態ｉは、特定のエネルギーレベル（すなわち、エネルギーレベル０、１又は２）における異なる速度ベクトルを表す。各状態の速度ｃｉは、以下のように３次元の各次元における「速度」と共に示される。
ｃi＝（ｃix，ｃiν，ｃiz） 方程式（Ｉ５）
ｃi＝（ｃix，ｃiν，ｃiz） 方程式（Ｉ５）
【0083】
エネルギーレベルゼロの状態は、どの次元においても動いていない停止した粒子を表し、すなわちｃstopped＝（０，０，０）である。エネルギーレベル１の状態は、３次元のうちの１つの次元において±１の速度を有し、他の２つの次元においてゼロ速度を有する粒子を表す。エネルギーレベル２の状態は、３次元全てにおいて±１の速度を有する粒子、或いは３次元のうちの１つの次元において±２の速度を有し、他の２つの次元においてゼロ速度を有する粒子を表す。
【0084】
３つのエネルギーレベルの考えられる順列を全て生成すると、全部で３９個の可能な状態（１つのエネルギー０状態、６つのエネルギー１状態、８つのエネルギー３状態、６つのエネルギー４状態、１２のエネルギー８状態及び６つのエネルギー９状態）が得られる。
【0085】
各ボクセル（すなわち、各格子サイト）は、状態ベクトルｆ（ｘ）によって表される。状態ベクトルは、ボクセルの状態を完全に定め、３９個のエントリを含む。３９個のエントリは、１つのエネルギー０状態、６つのエネルギー１状態、８つのエネルギー３状態、６つのエネルギー４状態、１２のエネルギー８状態及び６つのエネルギー９状態に対応する。システムは、この速度集合を使用することによって、達成された平衡状態ベクトルのためのマクスウェル－ボルツマン統計を生じることができる。全エネルギーでは、第２の分散ｑi（ｘ，ｔ）に同じ数のエントリ（３９の状態）を使用して比エネルギーを表す。これを状態ベクトルｆ（ｘ）と共に使用して全エネルギー保存を解く。
【0086】
ボクセルは、処理効率のためにマイクロブロックと呼ばれる２×２×２の容量にグループ化される。これらのマイクロブロックは、ボクセルの並行処理を可能にして、データ構造に関連するオーバーヘッドを最小化するように編成される。マイクロブロック内のボクセルの略語はＮi（ｎ）として定義され、ここでのｎは、マイクロブロック内の格子サイトの相対的位置を表し、ｎ∈｛０，１，２，．．．，７｝である。図１０にマイクロブロックを示す。
【0087】
図１１Ａ及び図１１Ｂを参照すると、表面Ｓ（図１１Ａ）をシミュレーション空間（図１１Ｂ）内にファセットＦαの集合として表している。
Ｓ＝｛Ｆα｝ 方程式（Ｉ６）
ここでのαは、特定のファセットを列挙する指数である。ファセットはボクセル境界に制限されず、比較的少数のボクセルにファセットが影響を与えるように、典型的にはファセットに隣接するボクセルのサイズと同程度又はそれよりもわずかに小さなサイズを有する。ファセットには、表面動力学を実装する目的で特性が割り当てられる。具体的に言えば、各ファセットＦαは、単位法線（ｎα）と、表面積（Ａα）と、中心位置（ｘα）と、ファセットの表面動特性を表すファセット分布関数（ｆi（α））とを有する。全エネルギー分散関数ｑi（α）は、ファセットとボクセルの相互作用のための流れ分散と同じ方法で取り扱われる。
Ｓ＝｛Ｆα｝ 方程式（Ｉ６）
ここでのαは、特定のファセットを列挙する指数である。ファセットはボクセル境界に制限されず、比較的少数のボクセルにファセットが影響を与えるように、典型的にはファセットに隣接するボクセルのサイズと同程度又はそれよりもわずかに小さなサイズを有する。ファセットには、表面動力学を実装する目的で特性が割り当てられる。具体的に言えば、各ファセットＦαは、単位法線（ｎα）と、表面積（Ａα）と、中心位置（ｘα）と、ファセットの表面動特性を表すファセット分布関数（ｆi（α））とを有する。全エネルギー分散関数ｑi（α）は、ファセットとボクセルの相互作用のための流れ分散と同じ方法で取り扱われる。
【0088】
図１２を参照すると、シミュレーション空間の異なる領域内で異なるレベルの分解能を用いて処理効率を高めることができる。通常は、物体６５５の周囲の領域６５０に最も関心があり、従ってこの領域を最高分解能でシミュレートする。粘度の効果は物体からの距離と共に減少するので、減少レベルの分解能（すなわち、拡大されたボクセル体積）を用いて、物体６５５から増加する距離に配置された領域６６０、６６５をシミュレートする。同様に、図１３に示すように、低レベルの分解能を用いて物体７７５のそれほど重要でない特徴の周囲の領域７７０をシミュレートする一方で、最高レベルの分解能を用いて物体７７５の最も重要な特徴（例えば、先端面及び後端面）の周囲の領域７８０をシミュレートすることもできる。中心から遠い領域７８５は、最低レベルの分解能及び最大ボクセルを用いてシミュレートされる。
【0089】
Ｃ．ファセットの影響を受けるボクセルの識別
再び図９を参照すると、シミュレーション空間をモデル化したら（ステップ３０２）、１又は２以上のファセットの影響を受けるボクセルを識別する（ステップ３０４）。ボクセルは、複数の形でファセットの影響を受けることができる。まず、１又は２以上のファセットが交わるボクセルは、交わっていないボクセルに比べて容量が減少するという点で影響を受ける。この理由は、ファセットと、ファセットが表す表面の下にある材料とがボクセルの一部を占有するからである。分数係数（ｆｒａｃｔｉｏｎａｌ ｆａｃｔｏｒ）Ｐf（ｘ）は、ファセットの影響を受けないボクセルの部分（すなわち、流れをシミュレートする対象である流体又は他の材料が占有できる部分）を示す。交わっていないボクセルでは、Ｐf（ｘ）が１に等しい。
再び図９を参照すると、シミュレーション空間をモデル化したら（ステップ３０２）、１又は２以上のファセットの影響を受けるボクセルを識別する（ステップ３０４）。ボクセルは、複数の形でファセットの影響を受けることができる。まず、１又は２以上のファセットが交わるボクセルは、交わっていないボクセルに比べて容量が減少するという点で影響を受ける。この理由は、ファセットと、ファセットが表す表面の下にある材料とがボクセルの一部を占有するからである。分数係数（ｆｒａｃｔｉｏｎａｌ ｆａｃｔｏｒ）Ｐf（ｘ）は、ファセットの影響を受けないボクセルの部分（すなわち、流れをシミュレートする対象である流体又は他の材料が占有できる部分）を示す。交わっていないボクセルでは、Ｐf（ｘ）が１に等しい。
【0090】
ファセットへの粒子の移動又はファセットからの粒子の受け取りを行うことによって１又は２以上のファセットと相互作用するボクセルも、ファセットの影響を受けるボクセルとして識別される。ファセットが交わる全てのボクセルは、ファセットから粒子を受け取る少なくとも１つの状態と、ファセットに粒子を移動させる少なくとも１つの状態とを含むようになる。ほとんどの場合、さらなるボクセルもこのような状態を含む。
【0091】
図１４を参照すると、ファセットＦαは、非ゼロの速度ベクトルｃｉを有する各状態ｉについて、速度ベクトルＣiとファセットの単位法線ｎαとのベクトルドット積（｜ｃiｎi｜）の大きさによって定められる高さと、ファセットの表面積Ａαによって定められる底面とを有する平行六面体Ｇiαによってその体積Ｖiαが次式に等しくなるように定められた領域から粒子を受け取り、又はこの領域に粒子を移動させる。
Ｖiα＝｜ｃiｎα｜Ａα 方程式（Ｉ７）
Ｖiα＝｜ｃiｎα｜Ａα 方程式（Ｉ７）
【0092】
ファセットＦαは、状態の速度ベクトルがファセットの方に向いている（｜ｃiｎi｜＜０の）時には体積Ｖiαから粒子を受け取り、状態の速度ベクトルがファセットから離れる方に向いている（｜ｃiｎi｜＞０の）時にはこの領域に粒子を移動させる。後述するように、内角などの非凸面特徴の近傍で起こり得る状態である、別のファセットが平行六面体Ｇiαの一部を占有する時には、この式が修正される。
【0093】
ファセットＦαの平行六面体Ｇiαは、複数のボクセルの一部又は全部に重なり合うことができる。ボクセルの数又はその一部は、ボクセルのサイズに対するファセットのサイズと、状態のエネルギーと、格子構造に対するファセットの配向とに依存する。影響を受けるボクセルの数は、ファセットのサイズと共に増加する。従って、上述したように、通常、ファセットのサイズは、ファセットの近くに位置するボクセルのサイズと同程度又はそれよりも小さくなるように選択される。
【0094】
平行六面体Ｇiαが重なり合ったボクセルの部分Ｎ（ｘ）は、Ｖiα（ｘ）として定められる。この項を使用すると、ボクセルＮ（ｘ）とファセットＦαとの間を移動する状態ｉの粒子の流束Γiα（ｘ）は、ボクセルの状態ｉの粒子の密度（Ｎi（ｘ））と、ボクセルが重なり合った領域の体積（Ｖiα（ｘ））とを乗算したものに等しい。
Γiα（ｘ）＝Ｎi（ｘ）Ｖiα（ｘ） 方程式（Ｉ８）
Γiα（ｘ）＝Ｎi（ｘ）Ｖiα（ｘ） 方程式（Ｉ８）
【0095】
平行六面体Ｇiαに１又は２以上のファセットが交わる時には、以下の条件が当てはまり、
Ｖiα＝ΣＶα（ｘ）＋ΣＶiα（β） 方程式（Ｉ９）
ここでの第１の加算は、Ｇiαが重なり合った全てのボクセルに相当し、第２項は、Ｇiαに交わる全てのファセットに相当する。平行六面体Ｇiαに別のファセットが交わらない時には、この式が以下のように変形する。
Ｖiα＝ΣＶiα（ｘ） 方程式（Ｉ１０）
Ｖiα＝ΣＶα（ｘ）＋ΣＶiα（β） 方程式（Ｉ９）
ここでの第１の加算は、Ｇiαが重なり合った全てのボクセルに相当し、第２項は、Ｇiαに交わる全てのファセットに相当する。平行六面体Ｇiαに別のファセットが交わらない時には、この式が以下のように変形する。
Ｖiα＝ΣＶiα（ｘ） 方程式（Ｉ１０）
【0096】
Ｄ．シミュレーションの実行
１又は２以上のファセットの影響を受けるボクセルが識別されると（ステップ３０４）、タイマを初期化してシミュレーションを開始する（ステップ３０６）。シミュレーションの各時間増分中には、粒子と表面ファセットとの相互作用を考慮した移流段階（ステップ３０８～３１６）によってボクセルからボクセルへの粒子の移動をシミュレートする。次に、上述したエントロピーソルバの使用に部分的に基づいて、衝突段階（ステップ３１８）によって各ボクセル内の粒子の相互作用をシミュレートする。その後、タイマを増分する（ステップ３２０）。増分したタイマによってシミュレーションが完了した旨が示されない（ステップ３２２）場合には、移流段階と衝突段階と（ステップ３０８～３２０）を繰り返す。増分したタイマによってシミュレーションが完了した旨が示された（ステップ３２２）場合には、シミュレーションの結果を記憶及び／又は表示する（ステップ３２４）。
１又は２以上のファセットの影響を受けるボクセルが識別されると（ステップ３０４）、タイマを初期化してシミュレーションを開始する（ステップ３０６）。シミュレーションの各時間増分中には、粒子と表面ファセットとの相互作用を考慮した移流段階（ステップ３０８～３１６）によってボクセルからボクセルへの粒子の移動をシミュレートする。次に、上述したエントロピーソルバの使用に部分的に基づいて、衝突段階（ステップ３１８）によって各ボクセル内の粒子の相互作用をシミュレートする。その後、タイマを増分する（ステップ３２０）。増分したタイマによってシミュレーションが完了した旨が示されない（ステップ３２２）場合には、移流段階と衝突段階と（ステップ３０８～３２０）を繰り返す。増分したタイマによってシミュレーションが完了した旨が示された（ステップ３２２）場合には、シミュレーションの結果を記憶及び／又は表示する（ステップ３２４）。
【0097】
１．表面の境界条件
表面との相互作用を正しくシミュレートするために、各ファセットは４つの境界条件を満たす。まず、ファセットが受け取った粒子の合計質量が、ファセットが移動させた粒子の合計質量に等しい（すなわち、ファセットへの正味質量流束がゼロに等しい）。次に、ファセットが受け取った粒子の合計エネルギーが、ファセットが移動させた粒子の合計エネルギーに等しい（すなわち、ファセットへの正味エネルギー流束がゼロに等しい）。これらの２つの条件は、各エネルギーレベル（すなわち、エネルギーレベル１及び２）での正味質量流束がゼロに等しくなるよう求めることによって満たすことができる。
表面との相互作用を正しくシミュレートするために、各ファセットは４つの境界条件を満たす。まず、ファセットが受け取った粒子の合計質量が、ファセットが移動させた粒子の合計質量に等しい（すなわち、ファセットへの正味質量流束がゼロに等しい）。次に、ファセットが受け取った粒子の合計エネルギーが、ファセットが移動させた粒子の合計エネルギーに等しい（すなわち、ファセットへの正味エネルギー流束がゼロに等しい）。これらの２つの条件は、各エネルギーレベル（すなわち、エネルギーレベル１及び２）での正味質量流束がゼロに等しくなるよう求めることによって満たすことができる。
【0098】
他の２つの境界条件は、ファセットと相互作用する粒子の正味運動量に関連する。本明細書では滑り面と呼ぶ表面摩擦がない表面では、正味接線運動量流束がゼロに等しく、正味法線運動量流束がファセットの局所圧力に等しい。従って、ファセットの法線ｎαに垂直な受け取った合計運動量の成分と、ファセットの法線ｎαに垂直な移動させた合計運動量の成分と（すなわち、接線成分）が等しく、ファセットの法線ｎαに平行な受け取った合計運動量の成分と、ファセットの法線ｎαに平行な移動させた合計運動量の成分と（すなわち、法線成分）の差分がファセットの局所圧力に等しい。非滑り面では、ファセットが受け取った粒子の合計接線運動量に対するファセットが移動させた粒子の合計接線運動量が、表面の摩擦によって摩擦量に関連する倍数だけ減少する。
【0099】
２．ボクセルからファセットへの収集
粒子と表面との間の相互作用をシミュレートする最初のステップとして、ボクセルから粒子を集めてファセットに提供する（ステップ３０８）。上述したように、ボクセルＮ（ｘ）とファセットＦαとの間の状態ｉの粒子の流束は以下の通りである。
Γiα（ｘ）＝Ｎi（ｘ）Ｖiα（ｘ） 方程式（Ｉ１１）
粒子と表面との間の相互作用をシミュレートする最初のステップとして、ボクセルから粒子を集めてファセットに提供する（ステップ３０８）。上述したように、ボクセルＮ（ｘ）とファセットＦαとの間の状態ｉの粒子の流束は以下の通りである。
Γiα（ｘ）＝Ｎi（ｘ）Ｖiα（ｘ） 方程式（Ｉ１１）
【0100】
この式から、各状態ｉがファセットの方を向いている（ｃiｎα＜０の）場合、ボクセルによってファセットＦαに提供される粒子の数は以下の通りである。
ΓiαV→F＝ΣxΓiα（ｘ）＝ΣxＮi（ｘ）Ｖiα（ｘ） 方程式（Ｉ１２）
ΓiαV→F＝ΣxΓiα（ｘ）＝ΣxＮi（ｘ）Ｖiα（ｘ） 方程式（Ｉ１２）
【0101】
Ｖiα（ｘ）が非ゼロの値を有するボクセルのみを加算する必要がある。上述したように、ファセットのサイズは、Ｖiα（ｘ）が少数のボクセルのみについて非ゼロの値を有するように選択される。Ｖiα（ｘ）及びＰf（ｘ）は非整数値を有することもできるので、Γα（ｘ）は実数として記憶され処理される。
【0102】
３．ファセットからファセットへの移動
次に、粒子がファセット間を移動する（ステップ３１０）。ファセットＦαが流入状態（ｃ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αの一部１０００が、ファセットＦαが交わる平行六面体Ｇ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α 方程式（Ｉ．１５）
【0103】
ファセット分布関数とも呼ばれるファセットの状態ベクトルＮ（α）は、ボクセル状態ベクトルのＭ個のエントリに対応するＭ個のエントリを有する。Ｍは、離散格子速度の数である。ファセット分布関数の入力状態Ｎ（α）は、これらの状態になる粒子の流束を体積Ｖiαで除算したものに等しく設定され、
ｃiｎα＜０の場合、以下のようになる。
Ｎi（α）＝ΓiIN（α）／Ｖiα 方程式（Ｉ．１６）
ｃiｎα＜０の場合、以下のようになる。
Ｎi（α）＝ΓiIN（α）／Ｖiα 方程式（Ｉ．１６）
【0104】
ファセット分布関数は、ファセットからの出力流束を生成するためのシミュレーションツールであり、必ずしも実際の粒子を表すものではない。正確な出力流束を生成するには、分布関数の他の状態に値を割り当てる。上述した内向き状態を投入する方法を用いて外向き状態を投入すると、
ｃiｎα≧０の場合、
Ｎi（α）＝ΓiOTHER（α）／Ｖiα 方程式（Ｉ．１７）
となり、ここでのΓiOTHER（α）は、上述したΓiIN（α）を生成するための方法を用いて、ただしこの方法を流入状態（ｃiｎα＜０））以外の状態（ｃiｎα≧０）に適用して決定される。別の方法では、以前の時間ステップからのΓiOUT（α）の値を用いて以下のようにΓiOTHER（α）を生成することもできる。
ΓiOTHER（α，ｔ）＝ΓiOUT（α，ｔ－１） 方程式（Ｉ．１８）
ｃiｎα≧０の場合、
Ｎi（α）＝ΓiOTHER（α）／Ｖiα 方程式（Ｉ．１７）
となり、ここでのΓiOTHER（α）は、上述したΓiIN（α）を生成するための方法を用いて、ただしこの方法を流入状態（ｃiｎα＜０））以外の状態（ｃiｎα≧０）に適用して決定される。別の方法では、以前の時間ステップからのΓiOUT（α）の値を用いて以下のようにΓiOTHER（α）を生成することもできる。
ΓiOTHER（α，ｔ）＝ΓiOUT（α，ｔ－１） 方程式（Ｉ．１８）
【0105】
平行状態（ｃiｎα＝０）では、Ｖiα及びＶiα（ｘ）がいずれもゼロである。Ｎi（α）を求める式では、Ｖiα（ｘ）が（ΓiOTHER（α）の式からの）分子に現れ、Ｖiαが（Ｎi（α）の式からの）分母に現れる。従って、Ｖiα及びＶiα（ｘ）がゼロに近付くと、平行状態のＮi（α）は、Ｎi（α）の極限として決定される。ゼロ速度を有する状態の値（すなわち休止状態、並びに状態（０，０，０，２）及び（０，０，０，－２））は、シミュレーションの最初に温度及び圧力の初期条件に基づいて初期化される。その後、これらの値は時間をかけて調整される。
【0106】
４．ファセット表面動力学の実行
次に、各ファセットが上述した４つの境界条件を満たすように表面動力学を実行する（ステップ３１２）。図１８に、ファセットの表面動力学の実行手順を示す。最初に、ファセットにおける粒子の合計運動量Ｐ（α）を以下のように求めることによって、ファセットＦαに垂直な合計運動量を求める（ステップ１１０５）。
全てのｉについて、
方程式（Ｉ．１９）
この式から、以下のように法線運動量Ｐn（α）を求める。
Ｐn（α）＝ｎα・Ｐ（α） 方程式（Ｉ．２０）
次に、各ファセットが上述した４つの境界条件を満たすように表面動力学を実行する（ステップ３１２）。図１８に、ファセットの表面動力学の実行手順を示す。最初に、ファセットにおける粒子の合計運動量Ｐ（α）を以下のように求めることによって、ファセットＦαに垂直な合計運動量を求める（ステップ１１０５）。
全てのｉについて、
【数38】
方程式（Ｉ．１９）
この式から、以下のように法線運動量Ｐn（α）を求める。
Ｐn（α）＝ｎα・Ｐ（α） 方程式（Ｉ．２０）
【0107】
次に、プッシング／プリング法を用いてこの法線運動量を排除（ステップ１１１０）してＮn-（α）を生成する。この方法によれば、粒子は、法線運動量にのみ影響を与える形で状態間を移行する。このプッシング／プリング法は、引用により組み入れられる米国特許第５，５９４，６７１号に記載されている。
【0108】
その後、Ｎn-（α）の粒子を衝突させてボルツマン分布Ｎn-β（α）を生じさせる（ステップ１１１５）。以下で流体力学の実行に関して説明するように、ボルツマン分布は、Ｎn-（α）に一連の衝突則を適用することによって達成することができる。
【0109】
次に、流入流束分布とボルツマン分布とに基づいてファセットＦαの流出流束分布を求める（ステップ１１２０）。まず、流入流束分布Γi（α）とボルツマン分布との差分を以下のように求める。
ΔΓi（α）＝ΓiIN（α）－Ｎn-βi（α）Ｖiα 方程式（Ｉ．２１）
ΔΓi（α）＝ΓiIN（α）－Ｎn-βi（α）Ｖiα 方程式（Ｉ．２１）
【0110】
この差分を使用すると、流出流束分布は、
ｎαｃ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は、状態ｉのエネルギー（ｊ）及び指示される接線ベクトルに対応する分布関数である。分布関数は、次式に従って決定され、
方程式（Ｉ．２４）
ここでのｊは、エネルギーレベル１状態では１に等しく、エネルギーレベル２状態では２に等しい。
ｎαｃ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は、状態ｉのエネルギー（ｊ）及び指示される接線ベクトルに対応する分布関数である。分布関数は、次式に従って決定され、
【数39】
方程式（Ｉ．２４）
ここでのｊは、エネルギーレベル１状態では１に等しく、エネルギーレベル２状態では２に等しい。
【0111】
ΓiOUT（α）の方程式の各項の関数は以下の通りである。第１項及び第２項は、ボルツマン分布を生じる上で衝突が効果的である限りは法線運動量流束境界条件を強制するが、接線運動量流束の例外を含む。第４項及び第５項は、この不十分な衝突に起因する離散性効果又は非ボルツマン構造を原因として生じ得る例外を補正する。最後に、第３項は、表面上の接線運動量流束の所望の変化を強制するように一定量の表面摩擦を追加する。摩擦係数Ｃfの生成については以下で説明する。なお、ベクトル操作に関与する全ての項は、シミュレーションの開始前に計算できる幾何学的因子である。
【0112】
この式から、接線速度が以下のように求められ、
ｕi（α）＝（Ｐ（α）－Ｐn（α）ｎα）／ρ 方程式（Ｉ．２５）
ここでのρは、以下などのファセット分布の密度である。
方程式（Ｉ．２６）
ｕi（α）＝（Ｐ（α）－Ｐn（α）ｎα）／ρ 方程式（Ｉ．２５）
ここでのρは、以下などのファセット分布の密度である。
【数40】
方程式（Ｉ．２６）
【0113】
上記と同様に、流入流束分布とボルツマン分布との差分を以下のように求める。
ΔΓi（α）＝ΓiIN（α）－Ｎn-βi（α）Ｖiα 方程式（Ｉ．２７）
ΔΓi（α）＝ΓiIN（α）－Ｎn-βi（α）Ｖiα 方程式（Ｉ．２７）
【0114】
すると、流出流束分布は以下のようになり、
ΓiOUT（α）＝Ｎn-βi（α）Ｖiα－ΔΓi＊（α）＋Ｃf（ｎαｃｉ）［Ｎn-βi＊（α）－Ｎn-βi（α）］Ｖiα 方程式（Ｉ．２８）
この式は、先の方法によって求めた流出流束分布の最初の２行に対応するが、異常な接線流束のための補正を必要としない。
ΓiOUT（α）＝Ｎn-βi（α）Ｖiα－ΔΓi＊（α）＋Ｃf（ｎαｃｉ）［Ｎn-βi＊（α）－Ｎn-βi（α）］Ｖiα 方程式（Ｉ．２８）
この式は、先の方法によって求めた流出流束分布の最初の２行に対応するが、異常な接線流束のための補正を必要としない。
【0115】
いずれかの方法を使用すると、結果として得られる流束分布は運動量流束条件を全て満たし、すなわち、
方程式（Ｉ．２９）
となり、ここでのｐαは、ファセットＦαにおける平衡圧であって、ファセットに粒子を提供するボクセルの平均密度及び温度値に基づき、ｕαは、ファセットにおける平均速度である。
【数41】
方程式（Ｉ．２９）
となり、ここでのｐαは、ファセットＦαにおける平衡圧であって、ファセットに粒子を提供するボクセルの平均密度及び温度値に基づき、ｕαは、ファセットにおける平均速度である。
【0116】
質量及びエネルギー境界条件が確実に満たされるように、各エネルギーレベルｊについて入力エネルギーと出力エネルギーとの差分を以下のように測定する。
方程式（Ｉ．３０）
ここでの添字ｊは、状態ｉのエネルギーを示す。次に、このエネルギー差を用いて差分項（ｄｉｆｆｅｒｅｎｃｅ ｔｅｒｍ）を生成すると、
ｃｊｉｎα＞０の場合、以下のようになる。
方程式（Ｉ．３１）
この差分項を用いて流出流束を修正すると、流速は、
ｃjiｎα＞０の場合、以下のようになる。
ΓαjiOUTf＝ΓαjiOUT＋δΓαji 方程式（Ｉ．３２）
この演算は、接線運動量流束をそのままにした状態で質量及びエネルギー流束を訂正する。ファセットの近傍において流れがおおよそ均一であって平衡に近い場合には、この調整はわずかなものである。結果として得られる調整後の法線運動量流束は、近傍の平均特性に近傍の非均一性又は非平衡特性に起因する補正をプラスしたものに基づく平衡圧である値へとわずかに変化する。
【数42】
方程式（Ｉ．３０）
ここでの添字ｊは、状態ｉのエネルギーを示す。次に、このエネルギー差を用いて差分項（ｄｉｆｆｅｒｅｎｃｅ ｔｅｒｍ）を生成すると、
ｃｊｉｎα＞０の場合、以下のようになる。
【数43】
方程式（Ｉ．３１）
この差分項を用いて流出流束を修正すると、流速は、
ｃjiｎα＞０の場合、以下のようになる。
ΓαjiOUTf＝ΓαjiOUT＋δΓαji 方程式（Ｉ．３２）
この演算は、接線運動量流束をそのままにした状態で質量及びエネルギー流束を訂正する。ファセットの近傍において流れがおおよそ均一であって平衡に近い場合には、この調整はわずかなものである。結果として得られる調整後の法線運動量流束は、近傍の平均特性に近傍の非均一性又は非平衡特性に起因する補正をプラスしたものに基づく平衡圧である値へとわずかに変化する。
【0117】
５．ボクセルからボクセルへの移動
再び図９を参照すると、３次元直線格子に沿ってボクセル間で粒子を移動させる（ステップ３１４）。このボクセルからボクセルへの移動は、ファセットと相互作用しないボクセル（すなわち、表面付近に存在しないボクセル）に対して行われる唯一の移動操作である。典型的なシミュレーションでは、表面と相互作用するほど十分に表面近くに存在しないボクセルがボクセルの大多数を構成する。
再び図９を参照すると、３次元直線格子に沿ってボクセル間で粒子を移動させる（ステップ３１４）。このボクセルからボクセルへの移動は、ファセットと相互作用しないボクセル（すなわち、表面付近に存在しないボクセル）に対して行われる唯一の移動操作である。典型的なシミュレーションでは、表面と相互作用するほど十分に表面近くに存在しないボクセルがボクセルの大多数を構成する。
【0118】
この分離状態の各々は、３次元のｘ、ｙ及びｚの各々において格子に沿って整数速度で移動する粒子を表す。整数速度は、０、±１及び±２を含む。速度の符号は、対応する軸に沿った粒子の移動方向を示す。
【0119】
表面と相互作用しないボクセルの場合、この移動操作は計算的に極めて単純である。毎時間増分中に、状態の集団全体をその現在のボクセルから宛先ボクセルに移動させる。同時に、宛先ボクセルの粒子を、そのボクセルから独自の宛先ボクセルに移動させる。例えば、＋１ｘ及び＋１ｙ方向に移動しているエネルギーレベル１の粒子（１，０，０）は、その現在のボクセルから、ｘ方向に＋１上回り他の方向が０であるボクセルに移動させる。この粒子は、最後には移動前の状態と同じ状態の宛先ボクセル（１，０，０）に行き着く。ボクセル内の相互作用は、他の粒子及び表面との局所的相互作用に基づいてその状態の粒子総数を変化させることがある。変化しなければ、粒子は格子に沿って同じ速度で同じ方向に移動し続ける。
【0120】
この移動操作は、１又は２以上の表面と相互作用するボクセルでは若干複雑になる。この結果、１又は２以上の端数粒子（ｆｒａｃｔｉｏｎａｌ ｐａｒｔｉｃｌｅｓ）がファセットに移動することがある。このような端数粒子のファセットへの移動により、ボクセルに端数粒子が留まるようになる。これらの端数粒子を、ファセットに占有されたボクセルに移動させる。
【0121】
移流した粒子の状態に比全エネルギーの値を加算して、及びその時間間隔にわたって移流しなかった粒子の状態から比全エネルギーの値を減算することによって、コンピュータが粒子の状態ベクトルを修正する、移流への全エネルギー修正３１５も示す（図３を参照）。
【0122】
図１５を参照すると、ボクセル９０５の状態ｉの粒子の一部９００がファセット９１０に移動した時には（ステップ３０８）、残り部分９１５をファセット９１０が位置するボクセル９２０に移動させ、ここから状態ｉの粒子をファセット９１０に向ける。従って、状態集団が２５に等しく、Ｖiα（ｘ）が０．２５に等しい（すなわち、ボクセルの４分の１が平行六面体Ｇiαに交わる）場合には、６．２５の粒子がファセットＦαに移動し、１８．７５の粒子を、ファセットＦαに占有されたボクセルに移動させる。１つのボクセルに複数のファセットが交わることもあるので、１又は２以上のファセットに占有されたボクセルに移動する状態ｉの粒子の数Ｎ（ｆ）は以下のようになり、
方程式（Ｉ．３３）
ここでのＮ（ｘ）はソースボクセルである。
【数44】
方程式（Ｉ．３３）
ここでのＮ（ｘ）はソースボクセルである。
【0123】
６．ファセットからボクセルへの散乱
次に、各ファセットからの流出粒子をボクセルに散乱させる（ステップ３１６）。基本的に、このステップは、ボクセルからファセットに粒子が移動する収集ステップの逆である。ファセットＦαからボクセルＮ（ｘ）に移動する状態ｉの粒子の数は以下のようになり、
方程式（Ｉ．３４）
ここでのＰf（ｘ）は、部分的なボクセルの体積縮小に相当する。この式から、各状態ｉについて、ファセットからボクセルに向かう粒子の総数Ｎ(x)は以下のようになる。
方程式（Ｉ．３５）
次に、各ファセットからの流出粒子をボクセルに散乱させる（ステップ３１６）。基本的に、このステップは、ボクセルからファセットに粒子が移動する収集ステップの逆である。ファセットＦαからボクセルＮ（ｘ）に移動する状態ｉの粒子の数は以下のようになり、
【数45】
方程式（Ｉ．３４）
ここでのＰf（ｘ）は、部分的なボクセルの体積縮小に相当する。この式から、各状態ｉについて、ファセットからボクセルに向かう粒子の総数Ｎ(x)は以下のようになる。
【数46】
方程式（Ｉ．３５）
【0124】
粒子をファセットからボクセルに散乱させ、これらの粒子を周囲のボクセルから移流した粒子と組み合わせて結果を整数化した後には、一定のボクセルの一定の方向がアンダーフロー（負になる）又はオーバーフロー（８ビット実装において２５５を超える）のいずれかになることがある。この結果、これらの量を許容できる値の範囲に収まるように切り捨てた後に、質量、運動量及びエネルギーが増加又は減少するようになる。
【0125】
このような発生を防ぐために、有害な状態を切り捨てる前に、限界を超えた質量、運動量及びエネルギーを蓄積する。状態を有するエネルギーでは、（アンダーフローに起因して）得られた又は（オーバーフローに起因して）失われた値に等しい質量が、同じエネルギーを有するランダムに（又は順次に）選択された状態に加え直され、これ自体がオーバーフロー又はアンダーフローに影響されることはない。この質量及びエネルギーの追加によって生じるさらなる運動量を蓄積して、切り捨てからの運動量に追加する。同じエネルギー状態のみに質量を追加することにより、質量カウンタがゼロに達した時に質量及びエネルギーがいずれも補正される。最後に、運動量積算器がゼロに戻るまでプッシング／プリング法を用いて運動量を補正する。
【0126】
７．流体力学の実行
流体力学を実行する（図９、３１８）。このステップは、マイクロ動力学操作又はボクセル内操作と呼ぶことができる。同様に、移流手順は、ボクセル間操作と呼ぶことができる。後述するマイクロ動力学操作を用いて、ファセットにおける粒子を衝突させてボルツマン分布を生成することもできる。
流体力学を実行する（図９、３１８）。このステップは、マイクロ動力学操作又はボクセル内操作と呼ぶことができる。同様に、移流手順は、ボクセル間操作と呼ぶことができる。後述するマイクロ動力学操作を用いて、ファセットにおける粒子を衝突させてボルツマン分布を生成することもできる。
【0127】
流体力学は、格子ボルツマン方程式モデルにおいて、ＢＧＫ衝突モデルとして知られている特定の衝突演算子によって保証される。この衝突モデルは、実際の流体系における分布の動力学を模倣する。衝突プロセスは、方程式１及び方程式２の右辺によって十分に説明することができる。移流ステップ後には、方程式３を用いた分布関数から流体系の保存量、具体的には密度、運動量及びエネルギーが取得される。これらの量から、方程式（２）にｆeqで示す平衡分布関数が方程式（４）によって十分に規定される。いずれも表１に列挙する速度ベクトルセットｃi、重みの選択は、方程式２と相まって、巨視的挙動が正しい流体動力学方程式に従うことを保証する。
【0128】
Ｅ．可変分解能
図１２を参照すると、（上述した）可変分解能は、以下では粗ボクセル１２０００及び微小ボクセル１２０５と呼ぶ異なるサイズのボクセルを採用する。（以下の説明では、２つの異なるサイズを有するボクセルを参照するが、説明する技術は、３又は４以上の異なるサイズのボクセルに適用してさらなる分解能レベルを提供することもできると理解されたい。）粗ボクセル及び微小ボクセルの領域間の界面は、可変分解能（ＶＲ）界面１２１０と呼ばれる。
図１２を参照すると、（上述した）可変分解能は、以下では粗ボクセル１２０００及び微小ボクセル１２０５と呼ぶ異なるサイズのボクセルを採用する。（以下の説明では、２つの異なるサイズを有するボクセルを参照するが、説明する技術は、３又は４以上の異なるサイズのボクセルに適用してさらなる分解能レベルを提供することもできると理解されたい。）粗ボクセル及び微小ボクセルの領域間の界面は、可変分解能（ＶＲ）界面１２１０と呼ばれる。
【0129】
可変分解能を表面又はその付近で採用すると、ファセットがＶＲ界面の両側のボクセルに相互作用することができる。これらのファセットは、ＶＲ界面ファセット１２１５（ＦαIC）又はＶＲ微小ファセット１２２０（ＦαIF）として分類される。ＶＲ界面ファセット１２１５は、ＶＲ界面の粗側に位置するファセットであり、微小ボクセル内に延びる粗平行六面体１２２５を有する。（粗平行六面体は、粗ボクセルの寸法に従ってｃｉが寸法決めされるものであり、微小平行六面体は、微小ボクセルの寸法に従ってｃｉが寸法決めされるものである。）ＶＲ微小ファセット１２２０は、ＶＲ界面の微小側に位置するファセットであり、粗ボクセル内に延びる微小平行六面体１２３０を有する。界面ファセットに関連する処理は、粗ファセット１２３５（ＦαC）及び微小ファセット１２４０（ＦαF）との相互作用を伴うこともできる。
【0130】
両タイプのＶＲファセットでは、表面動力学が微小スケールで実行され、上述したように動作する。しかしながら、ＶＲファセットは、ＶＲファセットとの間を粒子が移流する方法に関して他のファセットと異なる。
【0131】
ＶＲファセットとの相互作用は、図１３に示す可変分解能手順１３００を用いて行われる。この手順のほとんどのステップは、非ＶＲファセットとの相互作用について上述した同等のステップを用いて行われる。手順１３００は、それぞれが微小時間ステップに対応する２つの位相を含む粗時間ステップ（すなわち、粗ボクセルに対応する期間）中に実行される。ファセット表面動力学は、各微小時間ステップ中に実行される。このため、ＶＲ境界面ファセットＦαIcは、それぞれ黒色ファセットＦαIcb及び赤色ファセットＦαICrと呼ばれる２つの同一サイズ及び同一配向の微小ファセットとみなされる。黒色ファセットＦαICbは、粗時間ステップ内の第１の微小時間ステップに関連し、赤色ファセットＦαICrは、粗時間ステップ内の第２の微小時間ステップに関連する。
【0132】
最初に、第１の表面－表面移流段階によって粒子がファセット間を移動（移流）する（ステップ１３０２）。粒子は、黒色ファセットＦαICbから、Ｖ-αβの重み付け係数を有する粗ファセットＦβCに移動し、この重み付け係数Ｖ-αβは、ファセットＦαから広がってファセットＦβの裏側に存在する粗平行六面体（図１２、１２２５）の非ブロック部分から、ファセットＦαから広がってファセットＦβの裏側に存在する微小平行六面体（図１２、１２４５）の非ブロック部分を差し引いた容積に対応する。微小ボクセルのｃiの大きさは、粗ボクセルのｃiの大きさの１／２である。上述したように、ファセットＦαの平行六面体の容積は以下のように定義される。
Ｖiα＝｜ｃiｎα｜Ａα 方程式（Ｉ．３６）
Ｖiα＝｜ｃiｎα｜Ａα 方程式（Ｉ．３６）
【0133】
従って、ファセットの表面積Ａαは粗平行六面体と微小平行六面体との間で変化せず、単位法線ｎαは常に１の大きさを有するので、あるファセットに対応する微小平行六面体の容積は、このファセットの対応する粗平行六面体の容積の１／２である。
【0134】
粒子は、粗ファセットＦαＣから、Ｖαβの重み付け係数を有する黒色ファセットＦβICbに移動し、この重み付け係数Ｖαβは、ファセットＦαから広がってファセットＦβの裏側に存在する微小平行六面体の非ブロック部分の容積に対応する。
【0135】
粒子は、赤色ファセットＦαICrからＶαβの重み付け係数を有する粗ファセットＦβCに、そして粗ファセットＦαCからＶ-αβの重み付け係数を有する赤色ファセットＦβICbに移動する。
【0136】
粒子は、赤色ファセットＦαICrからＶαβの重み付け係数を有する黒色ファセットＦβICbに移動する。この段階では、黒色から赤色への移流は行われない。また、黒色ファセット及び赤色ファセットは連続する時間ステップを表すので、黒色から黒色への移流（又は赤色から赤色への移流）は決して行われない。同様の理由で、この段階の粒子は、赤色ファセットＦαICrからＶαβの重み付け係数を有する微小ファセットＦβIF又はＦβFに、そして微小ファセットＦαIF又はＦαFから同じ重み付け係数を有する黒色ファセットＦαICbに移動する。
【0137】
最後に、粒子は、微小ファセットＦαIF又はＦαFから同じ重み付け係数を有する他の微小ファセットＦβIF又はＦβFに、そして粗ファセットＦαCからＶCαβの重み付け係数を有する他の粗ファセットＦCに移動し、この重み付け係数ＶCαβは、ファセットＦαから広がってファセットＦβの裏側に存在する粗平行六面体の非ブロック部分の容積に対応する。
【0138】
粒子が表面間を移流した後に、第１の収集段階においてボクセルから粒子を収集する（ステップ１３０４～１３１０）。微小ファセットＦαFの粒子は、微小平行六面体を用いて微小ボクセルから収集され（ステップ１３０４）、粗ファセットＦαCの粒子は、粗平行六面体を用いて粗ボクセルから収集される（ステップ１３０６）。次に、微小平行六面体を用いて、粗ボクセル及び微小ボクセルの両方から黒色ファセットＦαIRb及びＶＲ微小ファセットＦαIFの粒子を収集する（ステップ１３０８）。最後に、粗平行六面体と微小平行六面体との差分を用いて、粗ボクセルから赤色ファセットＦαIRrの粒子を収集する（ステップ１３１０）。
【0139】
次に、微小ボクセル又はＶＲファセットと相互作用をする粗ボクセルを一群の微小ボクセルに展開する（ステップ１３１２）。単一の粗時間ステップ内で粒子を微小ボクセルに移動させる状態の粗ボクセルを展開する。例えば、ファセットが交わっていない適切な状態の粗ボクセルは、図４のマイクロブロックのように配向された８つの微小ボクセルに展開される。１又は２以上のファセットが交わっている適切な状態の粗ボクセルは、いずれのファセットも交わっていない粗ボクセルの部分に対応する一群の完全な及び／又は部分的な微小ボクセルに展開される。粗ボクセルと、粗ボクセルの展開によって生じる微小ボクセルとでは粒子密度Ｎi（ｘ）は等しいが、微小ボクセルは、粗ボクセルの分数係数及び他の微小ボクセルの分数係数とは異なる分数係数Ｐfを有することができる。
【0140】
その後、微小ファセットＦαIF及びＦαFについて表面動力学を実行し（ステップ１３１４）、黒色ファセットＦαICbについて表面動力学を実行する（ステップ１３１６）。動力学は、図１１に示して上述した手順を用いて実行される。
【0141】
次に、実際の微小ボクセルと、粗ボクセルの展開によって生じた微小ボクセルとを含む微小ボクセル間で粒子を移動させる（ステップ１３１８）。粒子を移動させたら、微小ファセットＦαIF及びＦαFから微小ボクセルに粒子を散乱させる（ステップ１３２０）。
【0142】
粒子は、黒色ファセットＦαICbから（粗ボクセルの展開によって生じた微小ボクセルを含む）微小ボクセルにも散乱させる（ステップ１３２２）。粒子は、微小ボクセルがその時に表面の存在が無ければ粒子を受け取っていたと思われる場合には微小ボクセルに散乱させる。具体的に言えば、（粗ボクセルの展開によって生じた微小ボクセルではなく）ボクセルが実際の微小ボクセルである時、ボクセルＮ（ｘ）を超える１つの速度単位Ｎ（ｘ）であるボクセルＮ（ｘ＋ｃi）が実際の微小ボクセルである時、或いはボクセルを超える１つの速度単位Ｎ（ｘ）であるボクセルＮ（ｘ＋ｃi）が粗ボクセルの展開によって生じた微小ボクセルである時には、粒子をボクセルＮ（ｘ）に散乱させる。
【0143】
最後に、微小ボクセル上で流体力学を実行することによって第１の微小時間ステップが完了する（ステップ１３２４）。流体力学を実行するボクセルは、粗ボクセルの展開によって生じた微小ボクセルを含まない（ステップ１３１２）。
【0144】
手順１３００は、第２の微小時間ステップ中にも同様のステップを実施する。最初に、第２の表面－表面移流段階において、粒子が表面間を移動する（ステップ１３２６）。粒子は、黒色ファセットから赤色ファセットに、黒色ファセットから微小ファセットに、微小ファセットから赤色ファセットに、そして微小ファセットから微小ファセットに移流する。
【0145】
粒子が表面間を移流した後に、第２の収集段階においてボクセルから粒子を収集する（ステップ１３２８～１３３０）。赤色ファセットＦαIRrの粒子は、微小平行六面体を用いて微小ボクセルから収集される（ステップ１３２８）。微小ファセットＦαF及びＦαIFの粒子も、微小平行六面体を用いて微小ボクセルから収集される（ステップ１３３０）。
【0146】
従って、上述したように、微小ファセットＦαIF及びＦαFについて表面動力学を実行し（ステップ１３３２）、粗ファセットＦαCについて表面動力学を実行し（ステップ１１３４）、赤色ファセットＦαICRについて表面動力学を実行する（ステップ１３３６）。
【0147】
次に、微小ボクセルと粗ボクセルを表す微小ボクセルとの間を粒子が移動するように、微小分解能を用いてボクセル間で粒子を移動させる（ステップ１３３８）。その後、粒子が粗ボクセル間を移動するように、粗分解能を用いてボクセル間で粒子を移動させる（ステップ１３４０）。
【0148】
次に、組み合わせステップにおいて、粒子をファセットからボクセルに散乱させながら、粗ボクセルを表す微小ボクセル（すなわち、粗ボクセルの展開によって生じた微小ボクセル）を粗ボクセルに合体させる（ステップ１３４２）。この組み合わせステップでは、粗平行六面体を用いて粒子を粗ファセットから粗ボクセルに散乱させ、微小平行六面体を用いて粒子を微小ファセットから微小ボクセルに散乱させ、微小平行六面体を用いて粒子を赤色ファセットから微小又は粗ボクセルに散乱させ、粗平行六面体と微小平行六面体との差分を用いて粒子を黒色ファセットから粗ボクセルに散乱させる。最後に、微小ボクセル及び粗ボクセルに流体力学を実行する（ステップ１３４４）。
【0149】
Ｆ．スカラー輸送ソルバ
上述したように、流体流を解くには、スカラー輸送の背景キャリアとしての役割を果たす様々なタイプのＬＢＭを適用することができる。シミュレーション中には、流体流及びスカラー輸送の両方がシミュレートされる。例えば、流体流は、格子ボルツマン（ＬＢ）法を使用してシミュレートされる流れであり、スカラー輸送は、本明細書ではスカラー輸送方程式と呼ばれる一連の分布関数を使用してシミュレートされる。
上述したように、流体流を解くには、スカラー輸送の背景キャリアとしての役割を果たす様々なタイプのＬＢＭを適用することができる。シミュレーション中には、流体流及びスカラー輸送の両方がシミュレートされる。例えば、流体流は、格子ボルツマン（ＬＢ）法を使用してシミュレートされる流れであり、スカラー輸送は、本明細書ではスカラー輸送方程式と呼ばれる一連の分布関数を使用してシミュレートされる。
【0150】
本明細書では、流体流をシミュレートするためのＬＢＭ法の詳細な説明を行っているが、以下は、スカラーシミュレーションと共に使用できる１つの流体流シミュレーション方法の例である。
方程式（Ｉ．３７）
【数47】
方程式（Ｉ．３７）
【0151】
ここで、ｆi（ｘ、ｔ）（ｉ＝１，．．．，１９）は粒子分布関数であり、τは単一の緩和時間であり、ｆi
eq（ｘ，ｔ）は流体速度の３次発展を含む平衡分布関数である。
方程式（Ｉ．３８）
ここではＴ0＝１／３である。離散格子速度ｃiは以下の通りであり、
方程式（Ｉ．３９）
残りの粒子ではｗ0＝１／３であり、デカルト方向の状態ではｗi＝１／１８であり、二重対角方向の状態ではｗi＝１／３６である。ｇi（ｘ，ｔ）は、外部体積力項である。流体力学的量￡及びｕは、粒子分布関数のモーメントである。
方程式（Ｉ．４０）
【数48】
方程式（Ｉ．３８）
ここではＴ0＝１／３である。離散格子速度ｃiは以下の通りであり、
【数49】
方程式（Ｉ．３９）
残りの粒子ではｗ0＝１／３であり、デカルト方向の状態ではｗi＝１／１８であり、二重対角方向の状態ではｗi＝１／３６である。ｇi（ｘ，ｔ）は、外部体積力項である。流体力学的量￡及びｕは、粒子分布関数のモーメントである。
【数50】
方程式（Ｉ．４０）
【0152】
上述したように、流体ソルバは、スカラー輸送情報を生成したスカラー輸送ソルバと共に使用される。従って、スカラー輸送では、流体ソルバに加えて別個の分布関数の組Ｔiが導入される。従って、システムは、システムのボクセル毎に流体流及びスカラー輸送の両方をシミュレートして、流体流を表す状態ベクトルと、スカラー変数を表すスカラー量とを生成する。これらのシミュレーション結果は、コンピュータアクセス可能媒体にエントリとして記憶される。
【0153】
スカラー輸送関数の組は、２次巨視的スカラー輸送方程式の間接的解法をもたらす。Ｔiは、スカラー量の動的発展をモデル化した方程式を提供する。
方程式（Ｉ．４１）
方程式（Ｉ．４２）
方程式（Ｉ．４３）
Ｔiはスカラー分布関数であり、Ｔは解かれるスカラーである。τTは、スカラー拡散率に対応する緩和時間である。緩和時間は、システムが緩和して平衡になるまでにどれほど掛かるかについての尺度を示す。項ｆi、ρ、Ｔ0及びｕは、それぞれ方程式（３８）、（３９）及び（４１）で定義される。
【数51】
方程式（Ｉ．４１）
【数52】
方程式（Ｉ．４２）
【数53】
方程式（Ｉ．４３）
Ｔiはスカラー分布関数であり、Ｔは解かれるスカラーである。τTは、スカラー拡散率に対応する緩和時間である。緩和時間は、システムが緩和して平衡になるまでにどれほど掛かるかについての尺度を示す。項ｆi、ρ、Ｔ0及びｕは、それぞれ方程式（３８）、（３９）及び（４１）で定義される。
【0154】
ｃiによって表される格子速度集合は、スカラーシミュレーションで使用される格子速度の離散集合とすることができる。一般に、スカラーソルバは、基本流体ソルバに付属する追加システムであるため、スカラー分布のための格子速度集合は、流体分布のための格子速度集合と同じである必要はない。例えば、スカラー発展のシミュレーションでは、流体流のシミュレーションよりも少ない数の格子速度を使用することができる。スカラー格子速度集合が流体格子速度集合の部分集合である限り、スカラーのための異なる格子速度集合を適用することができる。例えば、流体シミュレーションに１９速度ＬＢモデルを使用する場合、スカラーシミュレーションには６速度ＬＢモデルを使用することができる。１９速度ＬＢモデルは、６速度ＬＢモデルよりも高次の格子対称性を有するので、後述する例ではスカラーのための同じ１９速度格子モデルを使用する。
【0155】
（例えば、上述したような）標準的な周知のＢＧＫは、全ての次数の非平衡モーメントを含む。全ての非平衡モーメントを含めることは、必ずしも等方的、流体力学的又は物理的に有意義ではないと考えられる。従って、ＢＧＫ正則化／フィルタ処理した衝突演算子の形態を使用する。衝突演算子Φiは、将来的な衝突係数を表す。この衝突演算子は、関連するサポートされた次数のみ（例えば、１次のみ）において非平衡スカラー特性を抽出する。この演算子は、対象モードの維持及び緩和を行う一方で、サポートされていない／望ましくない高次モードに関連する非平衡特性は排除される。この計画は、スカラー輸送の物理的特性（例えば、移流及び拡散）の回復にとって十分なものである。この将来的な衝突演算子の使用は、ノイズを大幅に低減し、より良好な移流挙動を提供し（例えば、ガリレイ不変を高めることができ）、周知のＢＧＫ演算子の他の解と比べて安定すると考えられる。（例えば、方程式４３に示すような）このような形態は、流体力学的範囲におけるスカラー拡散に１次非平衡モーメントのみが寄与することを確実にする。さらに高次の非平衡モーメントは、全てこの衝突プロセスによってフィルタ除去される。上述した衝突演算子Φiの使用は、ＢＧＫにおいて示される数値ノイズの排除及びロバスト性の改善を含む利点をもたらすと考えられる。スカラーＴは独自の平衡としての役割を果たし、スカラー平衡分布関数の複雑な式は不要である。衝突演算子Φiの全体的な計算はかなり効率的である。さらに高次の非平衡モーメントをフィルタ処理すると、高次の均衡解に存在し得るエイリアシングを低減するという利点をさらにもたらすことができると考えられる。
【0156】
方程式（４２）における衝突は、スカラー保存則に従うことが分かる。方程式（４２）の両辺にｆ’i（ｘ，ｔ）＝ｆi（ｘ＋ｃi，ｔ＋１）を乗算すると、
方程式（Ｉ．４４）
となり、この結果、
方程式（Ｉ．４５）
となることが分かり、ここでのＴ’i（ｘ，ｔ）は方程式（４２）の右辺を示す。従って、スカラーを温度とみなした場合、スカラー衝突演算子は局所的なρＴを保存し、このことは局所的エネルギー保存の実現を意味する。Ｔiはｆiと共に伝播するので、移流中にエネルギー分布Ｅ＝ｆiＴiが完全に維持される。従って、ρＴの大域的保存が達成される。さらに、最も注目すべき点として、このスキームは、スカラーＴの均一性について厳密な不変性を維持する。これについては、あらゆる場所で
である場合、背景の流れ場にかかわらずその後の全ての時点においてあらゆる場所でφ（ｘ，ｔ）＝０かつ
であることが容易に分かる。この基本特性は、過去のどの格子ボルツマンスカラーモデルでも実証されていない。
【数54】
方程式（Ｉ．４４）
となり、この結果、
【数55】
方程式（Ｉ．４５）
となることが分かり、ここでのＴ’i（ｘ，ｔ）は方程式（４２）の右辺を示す。従って、スカラーを温度とみなした場合、スカラー衝突演算子は局所的なρＴを保存し、このことは局所的エネルギー保存の実現を意味する。Ｔiはｆiと共に伝播するので、移流中にエネルギー分布Ｅ＝ｆiＴiが完全に維持される。従って、ρＴの大域的保存が達成される。さらに、最も注目すべき点として、このスキームは、スカラーＴの均一性について厳密な不変性を維持する。これについては、あらゆる場所で
【数56】
である場合、背景の流れ場にかかわらずその後の全ての時点においてあらゆる場所でφ（ｘ，ｔ）＝０かつ
【数57】
であることが容易に分かる。この基本特性は、過去のどの格子ボルツマンスカラーモデルでも実証されていない。
【0157】
チャップマン・エンスコッグの拡張を使用すると、方程式（４２）は、以下の２次巨視的スカラー輸送方程式を回復させ、
方程式（Ｉ．４６）
κ＝（τT－１／２）Ｔ0であることが分かる。この均一性不変条件は、ρが∇Ｔの範囲外にあることを確実にする。
【数58】
方程式（Ｉ．４６）
κ＝（τT－１／２）Ｔ0であることが分かる。この均一性不変条件は、ρが∇Ｔの範囲外にあることを確実にする。
【0158】
境界条件
ＬＢＭの１つの実質的利点は、複雑な幾何形状に対応できる点である。任意の幾何形状の無摩擦／摩擦境界条件（ＢＣ）を達成するための一般化された体積ＬＢ表面アルゴリズムでは、質量が保存されて、境界上の接線運動量流束及び法線運動量流束が正確に実現される。局所的詳細釣り合い（ｌｏｃａｌ ｄｅｔａｉｌｅｄ ｂａｌａｎｃｅ）が完全に満たされる。この方法の直接的拡張として、スカラーの任意の幾何形状の断熱（ゼロスカラー流束）ＢＣを導出することができる。断熱ＢＣが実現されると、規定の有限流束ＢＣを達成することができる。
ＬＢＭの１つの実質的利点は、複雑な幾何形状に対応できる点である。任意の幾何形状の無摩擦／摩擦境界条件（ＢＣ）を達成するための一般化された体積ＬＢ表面アルゴリズムでは、質量が保存されて、境界上の接線運動量流束及び法線運動量流束が正確に実現される。局所的詳細釣り合い（ｌｏｃａｌ ｄｅｔａｉｌｅｄ ｂａｌａｎｃｅ）が完全に満たされる。この方法の直接的拡張として、スカラーの任意の幾何形状の断熱（ゼロスカラー流束）ＢＣを導出することができる。断熱ＢＣが実現されると、規定の有限流束ＢＣを達成することができる。
【0159】
他の点別のＬＢとは異なり、表面要素の離散化集合に対して境界条件が実施される。これらの区分的平面要素（ｐｉｅｃｅ－ｗｉｓｅ ｆｌａｔ ｓｕｒｆａｃｅ ｅｌｅｍｅｎｔｓ）は、共に湾曲形状を表す。図１４に示すように、粒子移流中には、各表面要素が隣接する流体セルからの流入粒子を収集する（ステップ１４１０）。流入分布ｆi
in，Ｔi
inに、下にある表面要素からの平行六面体と粒子移動方向におけるセルとの体積重複によって重み付けする（ステップ１４１２）。流入量を受け取った後に、以下の表面スカラーアルゴリズムを適用する（ステップ１４１４）。表面からの流出分布を求めるために、
方程式（Ｉ.４７）
であり、ここでは、
方程式（Ｉ.４８）
Ｐi（≡｜ｎ・ｃi｜Ａ）
ｃi・nｎ＜０，ｃiｎ＝－ｃi・Ｐi＝Ｐi・nである。
【数59】
方程式（Ｉ.４７）
であり、ここでは、
【数60】
方程式（Ｉ.４８）
Ｐi（≡｜ｎ・ｃi｜Ａ）
ｃi・nｎ＜０，ｃiｎ＝－ｃi・Ｐi＝Ｐi・nである。
【0160】
ここでのｎは、流体領域の方を示す表面法線であり、ｃi・nｎ＜０，ｃiｎ＝－ｃi・Ｐi（≡｜ｎ・ｃi｜Ａ）は、表面法線ｎ及び所与の表面要素の面積Ａに関連する粒子方向ｃｉにおける平行六面体の体積であり、明らかにＰi＝Ｐi*である。流出分布は、同じ表面移流プロセスに従って表面要素から流体セルに逆伝播される（ステップ１４１６）。上記の表面スカラー衝突が厳密なゼロの表面スカラー流束を達成することを示すことは困難ではない。流出方向にわたる総和をとると、流出スカラー流束は以下のようになる。
方程式（Ｉ．４９）
【数61】
方程式（Ｉ．４９）
【0161】
Ｐi＝Ｐi*及び方程式（４９）におけるＴinの定義に留意すると、右辺の第２の総和項はゼロである。また、質量流束保存に起因して
であるため、総流出スカラー流束は総流入スカラー流束と同じである。
方程式（Ｉ．５０）
【数62】
であるため、総流出スカラー流束は総流入スカラー流束と同じである。
【数63】
方程式（Ｉ．５０）
【0162】
従って、任意の幾何形状についてゼロの正味表面流束（断熱）ＢＣが完全に満たされる。表面上で外部スカラーソースＱ（ｔ）が指定される場合、方程式（４８）に直接ソース項を加えることができる。
方程式（Ｉ．５１）
【数64】
方程式（Ｉ．５１）
【0163】
境界条件が規定のスカラー量Ｔw（例えば、表面温度）を有する場合、以下に従って表面熱流束を計算することができる。
方程式（Ｉ．５２）
【数65】
方程式（Ｉ．５２）
【0164】
数値検証
図１５～図１８に、ＬＢスカラーソルバの能力をその数値精度、安定性、ガリレイ不変性、及びグリッド方位独立などに関して立証する４つのシミュレーション結果の組を示す。比較として、２つの異なる２次ＦＤスキームを用いた結果、ｖａｎ Ｌｅｅｒ型の流束制限スキーム（ｆｌｕｘ ｌｉｍｉｔｅｒ ｓｃｈｅｍｅ）、及び直接混合スキーム（中心スキームと１次風上スキーム（ｆｉｒｓｔ ｏｒｄｅｒ ｕｐｗｉｎｄ ｓｃｈｅｍｅｓ）との混合）も示す。
図１５～図１８に、ＬＢスカラーソルバの能力をその数値精度、安定性、ガリレイ不変性、及びグリッド方位独立などに関して立証する４つのシミュレーション結果の組を示す。比較として、２つの異なる２次ＦＤスキームを用いた結果、ｖａｎ Ｌｅｅｒ型の流束制限スキーム（ｆｌｕｘ ｌｉｍｉｔｅｒ ｓｃｈｅｍｅ）、及び直接混合スキーム（中心スキームと１次風上スキーム（ｆｉｒｓｔ ｏｒｄｅｒ ｕｐｗｉｎｄ ｓｃｈｅｍｅｓ）との混合）も示す。
【0165】
Ａ．剪断波減衰
第１のテストケースは、一定の均一流体流によって運ばれる温度剪断波減衰である。初期温度分布は、均一分布に空間正弦波変動が加わったものであり、格子波長Ｌ＝１６であり、大きさは
である。ＴＡは定数である。背景平均流の速度は０．２であり、熱拡散率κは０．００２である。このような低分解能及びκでは、数値安定性及び精度が十分に有効となり得る。背景流を伴わない温度減衰では、ＬＢスカラーソルバ及び有限差分法の両方が理論との優れた一致を示す。非ゼロの背景平均流では、ＬＢスカラーソルバが依然として正確に理論に匹敵することができる。しかしながら、ＦＤ結果は顕著な数値誤差を示す。図１５には、格子時間ステップ８１における温度プロファイルを示す。流束制限ＦＤスキームでは数値拡散がはっきりと見られるが、混合ＦＤスキームでは正確な温度プロファイルもその位置も維持することができない。
第１のテストケースは、一定の均一流体流によって運ばれる温度剪断波減衰である。初期温度分布は、均一分布に空間正弦波変動が加わったものであり、格子波長Ｌ＝１６であり、大きさは
【数66】
である。ＴＡは定数である。背景平均流の速度は０．２であり、熱拡散率κは０．００２である。このような低分解能及びκでは、数値安定性及び精度が十分に有効となり得る。背景流を伴わない温度減衰では、ＬＢスカラーソルバ及び有限差分法の両方が理論との優れた一致を示す。非ゼロの背景平均流では、ＬＢスカラーソルバが依然として正確に理論に匹敵することができる。しかしながら、ＦＤ結果は顕著な数値誤差を示す。図１５には、格子時間ステップ８１における温度プロファイルを示す。流束制限ＦＤスキームでは数値拡散がはっきりと見られるが、混合ＦＤスキームでは正確な温度プロファイルもその位置も維持することができない。
【0166】
Ｂ．体積熱源を有する傾斜チャネル
第２のテストケースは、異なる格子配向でのチャネル流の温度分布シミュレーションである。チャネル壁は自由滑り（ｆｒｅｅ－ｓｌｉｐ）であり、流体流は均一を維持し、結果としてＵ0＝０．２である。チャネル幅は５０（格子面間隔）であり、流量Ｒｅは２０００である。熱拡散率κは０．００５である。壁の温度はＴw＝１／３に固定される。バルク流体領域に一定体積熱源ｑ＝５×１０-6を適用する。流れは、流動別方向に周期的境界条件を有し、格子整列状況ではこの実現が容易である。図１６に示すようにチャネル（明色）が傾斜している時には、内外のチャネル境界が座標方向に完全に一致することによって再び流動別方向に周期的境界条件が実現されるようになる。格子非依存を立証するために、傾斜角を２６．５６度になるように選択する。第１のテストケースと同様に、チャネルが格子整列している時には、ＬＢスカラーソルバ及び２つのＦＤスキームを用いた温度分布が解析解と非常に良好に一致する。しかしながら、チャネルが傾斜している時には、ＦＤスキームからの結果が理論から大きく外れる。図１６に、傾斜チャネルにわたる温度分布のシミュレーション結果を示す。ＬＢ結果は、格子の配向に依存しないことがはっきりと示されている。傾斜した境界配向に関する勾配計算の処理の基本的困難性からＦＤ法からの誤差も生じ、従ってさらなる数値アーチファクトが導入されるようになる。ここに示すＢＣではＬＢスカラー粒子移流が正確になるので、格子配向に依存しないスカラー発展を達成することができる。
第２のテストケースは、異なる格子配向でのチャネル流の温度分布シミュレーションである。チャネル壁は自由滑り（ｆｒｅｅ－ｓｌｉｐ）であり、流体流は均一を維持し、結果としてＵ0＝０．２である。チャネル幅は５０（格子面間隔）であり、流量Ｒｅは２０００である。熱拡散率κは０．００５である。壁の温度はＴw＝１／３に固定される。バルク流体領域に一定体積熱源ｑ＝５×１０-6を適用する。流れは、流動別方向に周期的境界条件を有し、格子整列状況ではこの実現が容易である。図１６に示すようにチャネル（明色）が傾斜している時には、内外のチャネル境界が座標方向に完全に一致することによって再び流動別方向に周期的境界条件が実現されるようになる。格子非依存を立証するために、傾斜角を２６．５６度になるように選択する。第１のテストケースと同様に、チャネルが格子整列している時には、ＬＢスカラーソルバ及び２つのＦＤスキームを用いた温度分布が解析解と非常に良好に一致する。しかしながら、チャネルが傾斜している時には、ＦＤスキームからの結果が理論から大きく外れる。図１６に、傾斜チャネルにわたる温度分布のシミュレーション結果を示す。ＬＢ結果は、格子の配向に依存しないことがはっきりと示されている。傾斜した境界配向に関する勾配計算の処理の基本的困難性からＦＤ法からの誤差も生じ、従ってさらなる数値アーチファクトが導入されるようになる。ここに示すＢＣではＬＢスカラー粒子移流が正確になるので、格子配向に依存しないスカラー発展を達成することができる。
【0167】
Ｃ．傾斜チャネルにおける温度伝播
デカルト格子上の非格子整列壁近傍領域には近隣情報が存在しないため、ＦＤベースの方法にとって不可欠な局所勾配を正確に推定することは極めて困難である。さらに、ＦＤ法では風上情報に強く依存するため、スカラー移流の計算がさらに困難になり得る。対照的に、ＬＢスカラーソルバの境界処理は、上述したように正確な局所スカラー流束保存を達成する。このような壁近傍領域のスカラー移流を正確に計算することができる。証明として、３０度傾斜したチャネルにおける高温対流を実施する。固体壁において自由滑り（ｆｒｅｅ－ｓｌｉｐ）かつ断熱のＢＣを実施すると、流体速度は、チャネルに沿って一定のＵ0＝０．０９０９である。熱拡散率κは０．００２である。最初は、入口のＴ＝４／９を除き、あらゆる場所で温度が１／３である。その後に、この温度前線は、均一な背景流体流によって歪みを伴わずに下流位置まで対流するはずである。図１７に、格子時間ステップ２０００におけるチャネルにわたる計算された温度前線分布を示す。ＬＢスカラーソルバの温度前線は、ほぼ真っ直ぐな形状を維持する。一方で、２つのＦＤスキームからの温度前線は、壁近傍領域においてかなりの歪みを示す。ＬＢスカラー結果は、ＬＢスカラーソルバの数値拡散が小さくなったことを意味する最も薄い温度前線を示したとも言える。
デカルト格子上の非格子整列壁近傍領域には近隣情報が存在しないため、ＦＤベースの方法にとって不可欠な局所勾配を正確に推定することは極めて困難である。さらに、ＦＤ法では風上情報に強く依存するため、スカラー移流の計算がさらに困難になり得る。対照的に、ＬＢスカラーソルバの境界処理は、上述したように正確な局所スカラー流束保存を達成する。このような壁近傍領域のスカラー移流を正確に計算することができる。証明として、３０度傾斜したチャネルにおける高温対流を実施する。固体壁において自由滑り（ｆｒｅｅ－ｓｌｉｐ）かつ断熱のＢＣを実施すると、流体速度は、チャネルに沿って一定のＵ0＝０．０９０９である。熱拡散率κは０．００２である。最初は、入口のＴ＝４／９を除き、あらゆる場所で温度が１／３である。その後に、この温度前線は、均一な背景流体流によって歪みを伴わずに下流位置まで対流するはずである。図１７に、格子時間ステップ２０００におけるチャネルにわたる計算された温度前線分布を示す。ＬＢスカラーソルバの温度前線は、ほぼ真っ直ぐな形状を維持する。一方で、２つのＦＤスキームからの温度前線は、壁近傍領域においてかなりの歪みを示す。ＬＢスカラー結果は、ＬＢスカラーソルバの数値拡散が小さくなったことを意味する最も薄い温度前線を示したとも言える。
【0168】
Ｄ．レイリー－ベナール（Ｒａｙｌｅｉｇｈ－Ｂｅｒｎａｒｄ）自然対流
レイリー－ベナール（ＲＢ）自然対流は、数値ソルバの精度検証のための古典的なベンチマークである。ＲＢ自然対流は、ケース設定は単純であるが複雑な物理現象である。レイリー数Ｒａが特定の臨界値を超えると、システムが非流から対流に移行する。図１８に示すような現在の研究は、流体に作用する浮力がαρｇ（Ｔ－Ｔm）として定義されるブシネスク（Ｂｏｕｓｓｉｎｅｓｑ）近似下で行われ、ここでのαは熱膨張率であり、ｇは重力であり、Ｔmは上下境界の平均温度値である。最も不安定な波数は、Ｒaが臨界値Ｒａ＝１７０７．７６２を超えた時にκc＝３．１１７であり、この研究では分解能１０１×５０を使用する。ここで使用されるＰrは０．７１である。ＲBＲＢ対流が確立されると、２つのプレート間の熱伝達が大幅に強化される。この熱伝達の強化は、Ｎｕ＝１＋＜ｕyＴ＞Ｈ／κΔＴによって表される。
レイリー－ベナール（ＲＢ）自然対流は、数値ソルバの精度検証のための古典的なベンチマークである。ＲＢ自然対流は、ケース設定は単純であるが複雑な物理現象である。レイリー数Ｒａが特定の臨界値を超えると、システムが非流から対流に移行する。図１８に示すような現在の研究は、流体に作用する浮力がαρｇ（Ｔ－Ｔm）として定義されるブシネスク（Ｂｏｕｓｓｉｎｅｓｑ）近似下で行われ、ここでのαは熱膨張率であり、ｇは重力であり、Ｔmは上下境界の平均温度値である。最も不安定な波数は、Ｒaが臨界値Ｒａ＝１７０７．７６２を超えた時にκc＝３．１１７であり、この研究では分解能１０１×５０を使用する。ここで使用されるＰrは０．７１である。ＲBＲＢ対流が確立されると、２つのプレート間の熱伝達が大幅に強化される。この熱伝達の強化は、Ｎｕ＝１＋＜ｕyＴ＞Ｈ／κΔＴによって表される。
【0169】
図１９に、流体シミュレーションのスクリーンショットを示す。（従来の方法ではなく）上述した全エネルギー保存を使用した流体シミュレーションは、同様の流体シミュレーションの表現と、いずれかの慣習的な対応する計算データ出力とをもたらす。しかしながら、上述した方法を用いたこのような流体シミュレーションは、例えば実際の物理的対象などの対象をモデル化する際に他の方法よりも高速に及び／又はより少ない計算リソースを使用して流体シミュレーションを実行することができる。
【0170】
本明細書で説明した平衡分布は圧力項を含んでおらず、従ってその値が衝突演算中に常に正であり、比較的強力な安定性をもたらす。衝突状態がそのまま移流した後には、運動量に等温圧力勾配が生じ、全エネルギー方程式に圧力対流項が生じない。全Ｐエネルギーの保存は、停止状態を修正することによって行われ、移流後には、粒子状態に局所圧力項を加え直してからモーメントを計算して正しい圧力勾配をもたらす。
【0171】
複数の実装について説明した。それでもなお、本開示の趣旨及び範囲から逸脱することなく様々な修正を行うことができると理解されるであろう。従って、以下の特許請求の範囲には他の実装も含まれる。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11A】
【図11B】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】