特許 7296216 高速流のための格子ボルツマンベースのソルバ

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2023-06-14

(45)【発行日】2023-06-22

(54)【発明の名称】高速流のための格子ボルツマンベースのソルバ

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20230615BHJP

【ＦＩ】

G16Z99/00

【請求項の数】 21

【外国語出願】

(21)【出願番号】P 2019027119

(22)【出願日】2019-02-19

(65)【公開番号】P2019215848

(43)【公開日】2019-12-19

【審査請求日】2022-02-18

(31)【優先権主張番号】62/632,584

(32)【優先日】2018-02-20

(33)【優先権主張国・地域又は機関】US

(31)【優先権主張番号】16/274,403

(32)【優先日】2019-02-13

(33)【優先権主張国・地域又は機関】US

(73)【特許権者】

【識別番号】512293770

【氏名又は名称】ダッソー システムズ シムリア コーポレイション

(74)【代理人】

【識別番号】100094569

【弁理士】

【氏名又は名称】田中 伸一郎

(74)【代理人】

【識別番号】100103610

【弁理士】

【氏名又は名称】▲吉▼田 和彦

(74)【代理人】

【識別番号】100109070

【弁理士】

【氏名又は名称】須田 洋之

(74)【代理人】

【識別番号】100067013

【弁理士】

【氏名又は名称】大塚 文昭

(74)【代理人】

【識別番号】100086771

【弁理士】

【氏名又は名称】西島 孝喜

(74)【代理人】

【識別番号】100109335

【弁理士】

【氏名又は名称】上杉 浩

(74)【代理人】

【識別番号】100120525

【弁理士】

【氏名又は名称】近藤 直樹

(74)【代理人】

【識別番号】100139712

【弁理士】

【氏名又は名称】那須 威夫

(72)【発明者】

【氏名】プラディープ ゴパラクリシュナン

(72)【発明者】

【氏名】ラオヤン チャン

(72)【発明者】

【氏名】フードン チェン

【審査官】宮地 匡人

(56)【参考文献】

【文献】特表２０１６－５２８６２８（ＪＰ，Ａ）

【文献】特表２０１５－５０７７７８（ＪＰ，Ａ）

【文献】中国特許出願公開第１０６５２９０６３（ＣＮ，Ａ）

【文献】蔦原 道久，格子ボルツマン法の基礎と応用，日本機械学会論文集（B編），2011年12月，Vol.77 No.784，pp.2367-2378

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ１６Ｚ ９９／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

コンピュータ上で流体流をシミュレートする方法であって、
メッシュにわたる粒子の動きをモデル化するように、前記メッシュにわたる流体の活動をシミュレートするステップと、
前記メッシュ内の各メッシュ位置の状態ベクトルの組をコンピュータアクセス可能メモリに記憶するステップであって、前記状態ベクトル各々が、対応するメッシュ位置における可能な運動量状態のうちの特定の運動量状態に対応する複数のエントリを含む、ステップと、
近隣のメッシュ位置から流入する分布の組を衝突演算のために収集し、
前記コンピュータによって各位置のスカラー値を計算し、
前記衝突演算と熱源の追加との積として流出する分布を求め、
前記コンピュータが前記粒子の次のメッシュ位置への移流を一定期間にわたって実行することによって前記流体流を修正する、
ことによって、前記流体流のエントロピーの時間発展をシミュレートするステップと、
を含むことを特徴とする方法。

【請求項2】

前記流体流の活動をシミュレートするステップは、第１の離散格子速度の組に部分的に基づいて前記流体流をシミュレートするステップを含み、前記方法は、第２の離散格子速度の組に部分的に基づいて前記エントロピー量の時間発展をシミュレートするステップをさらに含む、
請求項１に記載の方法。

【請求項3】

前記第２の離散格子速度の組は、前記第１の離散格子速度の組と同じ格子速度である、
請求項２に記載の方法。

【請求項4】

前記コンピュータによって、流入する格子セットベクトルからの高次誤差項を計算するステップと、
前記高次誤差項の平均値を求めるステップと、
前記高次誤差項の前記平均値を衝突演算子から減算するステップと、
をさらに含む、請求項１に記載の方法。

【請求項5】

さらなる熱源項が計算されて第２の状態に追加され、前記方法は、
前記コンピュータによって、流体粘性による加熱及び流体伝導による加熱の影響を計算するステップと、
前記コンピュータによってエントロピー拡散を計算し、該エントロピー拡散を前記さらなる熱源項から除去するステップと、
をさらに含む、請求項１に記載の方法。

【請求項6】

前記方法は、前記コンピュータによって前記メッシュ内の前記メッシュ位置の物理量の組を計算するステップをさらに含む、
請求項１に記載の方法。

【請求項7】

格子ボルツマンエントロピーソルバが二次速度項を回避する、
請求項１に記載の方法。

【請求項8】

前記衝突演算子は、速度に二次項を含まない非平衡計算を伴う、
請求項４に記載の方法。

【請求項9】

特定のエントロピーｓを表すためのさらなる格子ベクトルの組ｑ_iを提供することによってエントロピーを解決するステップをさらに含み、これらの状態の時間発展が、

【数1】

によって与えられ、式中のｑ_iは格子ベクトルであり、ｘは方向であり、ｃ_iは状態の速度であり、Δｔは時間ｔの変化であり、

は平衡における格子ベクトルであり、τ_qは緩和時間であり、

は非平衡寄与であり、ｐは圧力であり、Ｑ_sは前記さらなる熱源である、
請求項５に記載の方法。

【請求項10】

前記衝突演算子は、

【数2】

に関連する、
請求項８に記載の方法。

【請求項11】

コンピュータ上で流体流をシミュレートするように構成されたコンピュータシステムであって、
１又は２以上のプロセッサ装置と、
前記１又は２以上のプロセッサ装置に動作可能に結合されて計算命令を記憶するメモリと、
を備え、前記計算命令は、前記１又は２以上のプロセッサ装置に、
メッシュにわたる粒子の動きをモデル化するように、前記メッシュにわたる流体の活動をシミュレートすることと、
前記メッシュ内の各メッシュ位置の状態ベクトルの組を前記メモリに記憶することであって、前記状態ベクトル各々が、対応するメッシュ位置における可能な運動量状態のうちの特定の運動量状態に対応する複数のエントリを含む、ことと、
近隣のメッシュ位置から流入する分布の組を衝突演算のために収集し、
前記コンピュータによって各位置のスカラー値を計算し、
前記衝突演算と熱源の追加との積として流出する分布を求め、
前記コンピュータが前記粒子の次のメッシュ位置への移流を一定期間にわたって実行することによって前記流体流を修正する、
ための命令によって、前記流体流のエントロピーの時間発展をシミュレートすることと、
を行わせる、
ことを特徴とするシステム。

【請求項12】

前記流体流の活動をシミュレートするための前記命令は、
第１の離散格子速度の組に部分的に基づいて前記流体流をシミュレートするための命令と、
第２の離散格子速度の組に部分的に基づいて前記エントロピー量の時間発展をシミュレートするための命令と、
を含む、請求項１１に記載のシステム。

【請求項13】

前記第２の離散格子速度の組は、前記第１の離散格子速度の組と同じ格子速度である、
請求項１２に記載のシステム。

【請求項14】

流入する格子セットベクトルからの高次誤差項を計算するための命令と、
前記高次誤差項の平均値を求めるための命令と、
前記高次誤差項の前記平均値を衝突演算子から減算するための命令と、
をさらに含む、請求項１１に記載のシステム。

【請求項15】

さらなる熱源項が計算されて第２の状態に追加され、前記システムは、
流体粘性による加熱及び流体伝導による加熱の影響を計算するための命令と、
エントロピー拡散を計算し、該エントロピー拡散を前記さらなる熱源項から除去するための命令と、
をさらに含む、請求項１１に記載のシステム。

【請求項16】

前記メッシュ内の前記メッシュ位置の物理量の組を計算する命令をさらに含む、請求項１１に記載のシステム。

【請求項17】

格子ボルツマンエントロピーソルバが二次速度項を回避する、
請求項１１に記載のシステム。

【請求項18】

前記衝突演算子は、速度に二次項を含まない非平衡計算を伴う、
請求項１４に記載のシステム。

【請求項19】

特定のエントロピーｓを表すためのさらなる格子ベクトルの組ｑ_iを提供することによってエントロピーを解決するための命令をさらに含み、これらの状態の時間発展が、

【数3】

によって与えられ、式中のｑ_iは格子ベクトルであり、ｘは方向であり、ｃ_iは状態の速度であり、Δｔは時間ｔの変化であり、

は平衡における格子ベクトルであり、τ_qは緩和時間であり、

は非平衡寄与であり、ｐは圧力であり、Ｑ_sは前記さらなる熱源である、
請求項１５に記載のシステム。

【請求項20】

前記衝突演算子は、

【数4】

に関連する、
請求項１８に記載のシステム。

【請求項21】

非一時的ハードウェア記憶装置に有形的に記憶された、コンピュータ上で流体流をシミュレートするための実行可能命令を含むコンピュータプログラムであって、前記実行可能命令は、
メッシュにわたる粒子の動きをモデル化するように、前記メッシュにわたる流体の活動をシミュレートすることと、
前記メッシュ内の各メッシュ位置の状態ベクトルの組をメモリに記憶することであって、前記状態ベクトル各々が、対応するメッシュ位置における可能な運動量状態のうちの特定の運動量状態に対応する複数のエントリを含む、ことと、
近隣のメッシュ位置から流入する分布の組を衝突演算のために収集し、
前記コンピュータによって各位置のスカラー値を計算し、
前記衝突演算と熱源の追加との積として流出する分布を求め、
前記コンピュータが前記粒子の次のメッシュ位置への移流を一定期間にわたって実行することによって前記流体流を修正する、
ための命令によって、前記流体流のエントロピーの時間発展をシミュレートすることと、
を行うようにコンピュータシステムを構成する、
ことを特徴とするコンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

〔合衆国法典第３５編第１１９条（ｅ）に基づく優先権の主張〕
本出願は、２０１８年２月２０日に出願された「高速流のための格子ボルツマンベースのソルバ（Ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｂａｓｅｄ ｓｏｌｖｅｒ ｆｏｒ ｈｉｇｈ ｓｐｅｅｄ ｆｌｏｗｓ）」という名称の米国仮特許出願第６２／６３２，５８４号の合衆国法典第３５編第１１９条（ｅ）に基づく優先権を主張するものであり、この文献の内容はその全体が引用により本明細書に組み入れられる。

【背景技術】

【0002】

格子ボルツマン法（ＬＢＭ）は、流体シミュレーションのための計算流体力学（ＣＦＤ）法の一種である。ナビエ・ストークス方程式を解く代わりに、Ｂｈａｔｎａｇａｒ－Ｇｒｏｓｓ－Ｋｒｏｏｋ（ＢＧＫ）などの衝突モデルを用いて離散ボルツマン方程式を解いてニュートン流体の流れをシミュレートする。限定数の粒子にわたる流動過程及び衝突過程をシミュレートすることにより、固有の粒子反応が、さらに大きな質量にわたって適用可能な粘性流挙動の縮図を明らかにする。

【先行技術文献】

【特許文献】

【0003】

【文献】米国特許出願公開第２０１３／０１５１２２１号明細書

【文献】米国特許出願公開第２０１６／０１８８７６８号明細書

【文献】米国特許第９，５７６，０８７号明細書

【文献】米国特許第５，５９４，６７１号明細書

【発明の概要】

【課題を解決するための手段】

【0004】

ある態様によれば、コンピュータ上で流体流をシミュレートする方法が、メッシュにわたる粒子の動きをモデル化するように、メッシュにわたる流体の活動をシミュレートするステップと、メッシュ内の各メッシュ位置の状態ベクトルの組をコンピュータアクセス可能メモリに記憶するステップであって、状態ベクトル各々が、対応するメッシュ位置における可能な運動量状態のうちの特定の運動量状態に対応する複数のエントリを含む、ステップと、近隣のメッシュ位置から流入する分布の組を衝突演算のために収集し、コンピュータによって各位置のスカラー値を計算し、衝突演算と熱源の追加との積として流出する分布を求め、コンピュータが粒子の次のメッシュ位置への移流（advection）を一定期間にわたって実行することによって流れを修正することによって、流れのエントロピーの時間発展（time evolution）をシミュレートするステップと、を含む。

【0005】

上記の態様には、以下の特徴のうちの１つ又は２つ以上を含めることができる。

【0006】

流体流の活動をシミュレートするステップは、第１の離散格子速度の組に部分的に基づいて流体流をシミュレートするステップを含み、方法は、第２の離散格子速度の組に部分的に基づいてエントロピー量の時間発展をシミュレートするステップをさらに含む。第２の離散格子速度の組は、第１の離散格子速度の組と同じ格子速度である。方法は、コンピュータによって、流入する格子セットベクトルからの高次誤差項を計算するステップと、高次誤差項の平均値を求めるステップと、高次誤差項の平均値を衝突演算子から減算するステップと、をさらに含む。さらなる熱源項が計算されて第２の状態に追加され、方法は、コンピュータによって、流体粘性による加熱及び流体伝導による加熱の影響を計算するステップと、コンピュータによってエントロピー拡散（entropy diffusion）を計算し、エントロピー拡散をさらなる熱源項から除去するステップと、をさらに含む。方法は、コンピュータによってメッシュ内のメッシュ位置の物理量の組を計算するステップをさらに含む。方法は、二次速度項を回避する格子ボルツマンエントロピーソルバを提供する。衝突演算子は、速度に二次項を含まない非平衡計算を伴う。方法は、特定のエントロピー（specific entropy）ｓを表すためのさらなる格子ベクトルの組ｑ_iを提供することによってエントロピーを解決するステップをさらに含み、これらの状態の時間発展が、

【数1】

によって与えられ、式中のｑ_iは格子ベクトルであり、ｘは方向であり、ｃ_iは状態の速度であり、Δｔは時間ｔの変化であり、

は平衡における格子ベクトルであり、τ_qは緩和時間（relaxation time）であり、

は非平衡寄与であり、ｐは圧力であり、Ｑ_sはさらなる熱源である。衝突演算子は、

【数2】

に関連する。

【0007】

他の態様は、コンピュータプログラム製品と、コンピュータ及びデータ処理システムなどの装置とを含む。

【0008】

上記の態様のうちの１つ又は２つ以上は、以下の利点のうちの１つ又は２つ以上を提供することができる。エントロピーソルバは、Ｍａ²（速度の二乗）に比例するという理由で特に高速流の場合に不安定性などの誤差をもたらすはずの項である（ρｖｖ）を回避する。従って、このエントロピーソルバは、高速用途で安定する。格子ボルツマン法の文脈では、高度に正確な過渡的結果と圧縮性の影響とを必要とする高速用途を解決するためにＬＢベースのエントロピーソルバが有効である。高温度比（high temperature ratio）用途も、このＬＢエントロピーソルバの独自の衝突モデルによって安定性が高まるので、その恩恵を受ける。

【図面の簡単な説明】

【0009】

【図1】圧縮性流れ（compressible flows）のためのエントロピーソルバを含む流体流（fluid flows）シミュレーションシステムを示す図である。

【図2】エントロピーソルバを含む格子ボルツマンモデルシミュレーションを構築する動作を示すフローチャートである。

【図3】エントロピーソルバを含む格子ボルツマンモデルシミュレーションを用いたシミュレーション動作を示すフローチャートである。

【図4A】流体力学ベースの方法のシミュレートした流れを示す図である。

【図4B】格子ボルツマンベースの方法のシミュレートした流れを示す図である。

【図5】音圧レベル（ＳＰＬ）対周波数のグラフである。

【図6A】シミュレートした指向角対周波数及び音圧レベル（simulated directivity angle vs. frequency and sound pressure level）のプロットである。

【図6B】シミュレートした指向角対周波数及び音圧レベルのプロットである。

【図7】ユークリッド空間で表す１つのＬＢＭモデルの速度成分を示す図である（先行技術）。

【図8】ユークリッド空間で表す別のＬＢＭモデルの速度成分を示す図である（先行技術）。

【図9】物理プロセスシミュレーションシステムが従う手順のフローチャートである。

【図10】マイクロブロックの透視図である（先行技術）。

【図11A】図１のシステムが使用する格子構造の図である（先行技術）。

【図11B】図１のシステムが使用する格子構造の図である（先行技術）。

【図12】可変分解能技術を示す図である（先行技術）。

【図13】可変分解能技術を示す図である（先行技術）。

【図14】粒子の動きを示す図である（先行技術）。

【図15】表面のファセットの影響を受ける領域を示す図である（先行技術）。

【図16】表面動力学のフローチャートを示す図である（先行技術）。

【図17】表面動力学の実行手順のフローチャートである（先行技術）。

【図18】ボクセルから表面への粒子の動きを示す図である。

【図19】流体シミュレーションを表すスクリーンショットである（先行技術）。

【発明を実施するための形態】

【0010】

添付図面及び以下の説明には、本発明の１又は２以上の実施形態の詳細を示す。説明及び図面、並びに特許請求の範囲からは、他の特徴、目的及び利点が明らかである。

【0011】

流体流をシミュレートする１つの方法は、いわゆる格子ボルツマンモデル（ＬＢＭ）である。ＬＢＭベースの物理プロセスシミュレーションシステムでは、分布関数の時間発展を表す周知の格子ボルツマン方程式（以下の方程式１を参照）を用いて一連の離散速度で評価した分布関数値によって流体流が表現される。この分布関数は、以下の流動過程（streaming process）及び衝突過程（collision process）という２つの過程を伴う。

【数3】

（方程式１）

【数4】

（方程式２）
式中、ξ及びｕは、それぞれ無次元格子ベクトル及び速度であり、ｆ_i（ｘ，ｔ）は、時点ｔにおける位置ｘの格子方向ｃｉに沿った粒子分布関数であり、

は平衡分布であり、τは特徴的緩和時間である。

【0012】

流動過程は、あるメッシュ位置において流体ポケットが開始した後に、速度ベクトルの１つに沿って次のメッシュ位置に移動する時点である。この時点で、「衝突係数」、すなわち開始流体ポケットに対する近隣の流体ポケットの影響を計算する。流体は別のメッシュ位置に移動することしかできず、従って全ての速度の全ての成分が共通速度の倍数になるように正しい速度ベクトルの選択が必要である。

【0013】

第１の方程式の右辺は、流体ポケット同士の衝突に起因する分布関数の変化を表す上述した「衝突演算子」である。特定の形の衝突演算子は、Ｂｈａｔｎａｇａｒ－Ｇｒｏｓｓ－Ｋｒｏｏｋ（ＢＧＫ）演算子である。衝突演算子は、「平衡」形態である第２の方程式によって与えられる規定値に分布関数を至らしめる。

【0014】

ＢＧＫ演算子は、衝突の詳細にかかわらず、分布関数が以下の式に従って明確に定義された局所平衡に近づくように構築され、

【数5】

【数6】

ここでのパラメータτは、衝突を通じた平衡までの特徴的な緩和時間を表す。粒子（例えば、原子又は分子）を扱う場合、この緩和時間は、通常は定数として解釈される。

【0015】

このシミュレーションから、質量ρ及び流体速度ｕなどの従来の流体変数が方程式２の単純和として得られる。

【数7】

（方程式３）
ここでのρ及びｕは、それぞれ上述した流体密度及び流体速度である。

【0016】

対称性を考慮した結果、速度値の組は、構成空間内に広がった時に特定の格子構造を形成するように選択される。このような離散系の動特性は、以下の形を有する格子ボルツマン方程式（ＬＢＥ）に従い、
ｆ_i（ｘ＋ｃ_iΔｔ，ｔ＋Δｔ）－ｆ_i（ｘ，ｔ）＝Ｃ_i（ｘ，ｔ）
通常、ここでの衝突演算子は上述したＢＧＫ形態を取る。

【0017】

正しい平衡分布形態の選択により、格子ボルツマン方程式が正しい流体力学をもたらすことを理論的に示すことができる。すなわち、ｆ_i（ｘ，ｔ）から導出される流体力学的モーメントは、巨視的極限においてナビエ・ストークス方程式に従う。これらのモーメントは、上記の方程式（３）によって定められる。

【0018】

ｃ_i及びｗ_iの集合的な値は、ＬＢＭモデルを定める。ＬＢＭモデルは、スケーラブルなコンピュータプラットフォーム上に効率的に実装することができ、時間的に不安定な流れ及び複雑な境界条件の場合に高いロバスト性を持って動作することができる。

【0019】

ボルツマン方程式から流体系の巨視的運動方程式を得るための標準的な方法は、完全なボルツマン方程式の逐次近似を行うチャップマン・エンスコッグ（Ｃｈａｐｍａｎ－Ｅｎｓｋｏｇ）法である。流体系では、わずかな密度の変動が音速で伝わる。気体系では、一般に音速が温度によって決まる。流れにおける圧縮率の影響の重要性は、マッハ数として知られている特性速度と音速との比率によって測定される。従来のＬＢＭベースの物理プロセスシミュレーションシステムのさらなる説明については、上記で引用により組み入れた米国特許出願公開第２０１６／０１８８７６８号に測定値が示されている。

【0020】

熱格子ボルツマンモデル
方程式１及び方程式２によって表されるＬＢＭは、均一温度の低速流を解く場合にのみ適用可能である。熱流では、単独のソルバを導入してエネルギー方程式を解くことが慣習であり、計算された温度は格子ボルツマン法に逆結合される。温度方程式は、有限差分法、又はさらなる格子ベクトルの組の導入のいずれかによって解くことができる。

【0021】

固有のＬＢＭベースの温度ソルバを使用する。このＬＢＭベースの温度ソルバでは、さらなる格子ベクトルの組を用いて特定のスカラー、すなわち流れベクトルに同伴して１つのメッシュ点から別のメッシュ点に移動する温度を表す。このソルバには複数の利点があり、米国特許出願公開第２０１３／０１５１２２１号にさらに詳細な説明が記載されており、この文献の内容はその全体が引用により本明細書に組み入れられる。

【0022】

本明細書では、後述するエントロピーソルバに同様の輸送機構を使用する。

【0023】

温度が分かると、この温度を流体ＬＢＭソルバに結合する。従来の方法は、ＬＢＭ（方程式１）の時間発展において、熱流に以下の平衡分布（方程式４）を使用する。

【数8】

（方程式４）
ここでのξ、ｕ及びθは、それぞれ無次元格子ベクトル、無次元速度及び無次元圧力である。しかしながら、高温度比用途では、θ－１の項が負になることによってソルバが不安定になる。

【0024】

比較的広い範囲の流れ条件で衝突中の平衡分布を正にする結合法の具体的な詳細は、Ｐｒａｄｅｅｐ Ｇｏｐａｌａｋｒｉｓｈｎａｎ他による「ハイブリッド熱格子ボルツマン法のための温度結合アルゴリズム（Ｔｅｍｐｅｒａｔｕｒｅ Ｃｏｕｐｌｉｎｇ Ａｌｇｏｒｉｔｈｍ ｆｏｒ Ｈｙｂｒｉｄ Ｔｈｅｒｍａｌ Ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ Ｍｅｔｈｏｄ）」という名称の同時係属中の米国特許出願公開第２０１６／０１８８７６８号に詳述されており、この文献の内容はその全体が引用により本明細書に組み入れられる。

【0025】

方程式１に示すＢＧＫ衝突演算子は、比較的低い安定範囲を有する。方程式１の演算子は、低マッハ（マッハ＜０．６）及び低温度比用途では有用であるが、高マッハ（マッハ＞０．６）及び高温度比用途では不安定になる傾向にある。航空宇宙分野における高亜音速及び遷音速用途に適用可能な温度ソルバは、比較的高い安定範囲を有するガリレイ不変演算子に基づくさらに安定した衝突演算子を必要とする。このような衝突演算子の具体的な詳細は米国特許第９，５７６，０８７号に詳述されており、この文献の内容はその全体が引用により本明細書に組み入れられる。

【0026】

高速圧縮性流れのエネルギー保存
高速流、特にマッハ数が０．３を上回る流れでは、流体の圧縮性及び熱変化粘性に起因する流体エネルギーに対する影響を無視することができない。この圧縮性の影響は、以下の方程式５によって与えられる圧縮性流れのエネルギー方程式によって明らかになる。

【数9】

（方程式５）
ここでのｅは内部エネルギーであり、Ｖは速度であり、ｋは熱伝導率であり、Φは、粘性に起因して機械的エネルギーが熱エネルギーに変化する速度を表す散逸関数である。非圧縮性流れでは、方程式５の左辺の第３項である圧力作用項と粘性散逸項とが脱落し、従って方程式５が、比較的解きやすい典型的な対流拡散スカラー方程式に変化する。

【0027】

しかしながら、有限差分法又はその他の手法を用いて方程式５全体を解くと、発散速度項の存在に起因して不安定性が生じる。この項を避けるために、実行可能性の高いオプション、すなわち以下のエントロピー方程式（方程式６）を解くことによってエネルギー保存を強制する。

【数10】

（方程式６）
ここでのｓは、この系のエントロピーである。

【0028】

圧縮性流れのための格子ボルツマンベースのエントロピーソルバ
エントロピー方程式である方程式６は、有限体積法又は有限差分法のような従来の方法を用いて解くことができる。複雑な問題のために、メッシュが複数の分解能レベル及び複雑な境界条件を有する。従来の方法では、これらの状況の勾配を計算することが困難であることにより、強いメッシュ依存性及び数値的アーチファクトが導入される。

【0029】

本明細書では、これらの複雑な問題に対処することができる格子ボルツマン（ＬＢ）ベースのエントロピーソルバ３４ａ（図１）について説明する。エントロピーの解決は、特定のエントロピー「ｓ」を表すために導入されるさらなる第２の格子ベクトルの組ｑｉを使用する。これらの状態の時間発展は以下によって与えられ、

【数11】

（方程式７）
ここでのｐは圧力であり、Ｑ_sは、例えばさらなる熱源などのさらなるソースである。方程式７及び方程式１は、１つのメッシュ点から別のメッシュ点にエントロピー項を伝え、これらの２つの方程式の積として以下のエントロピー保存方程式（方程式８）が得られる。

【数12】

（方程式８）
ここでのＳは総エントロピーである。米国特許出願公開第２０１６／０１８８７６８号に記載されるように、この方法は、ＬＢＭベースの温度ソルバに使用されている機構と同様の、移流する流れに対するスカラー項として総エントロピーを輸送する。このエントロピー輸送方法には、異なるメッシュレベル付近及び複雑な境界付近の比較的滑らかな遷移を含む複数の利点がある。

【0030】

方程式７の右辺の第２項は、いわゆる正規化衝突演算子である。（方程式１によって与えられる）ＢＧＫ衝突演算子では、「非平衡部分」として知られている実際の状態と平衡状態との間の差分が計算される。この非平衡部分は、実際の状態が平衡状態まで緩和するように緩和される。

【0031】

しかしながら、この演算子は、保存にとって不要なあらゆる次数の非平衡効果を伴う。また、あらゆる次数の非平衡効果を伴えば、不安定性が生じることもある。一方、方程式７で使用されている演算子は、以下のように非平衡寄与の一次効果（格子ベクトルを乗算することによって計算される一次モーメント）しか考慮していない。

【数13】

（方程式９）
その後、このフィルタリングされた非平衡部分にガリレイ不変演算子ｃ_i－Ｖを乗算する。

【0032】

ここで図１を参照しながら、高速の圧縮性流れのためのＬＢベースのエントロピーソルバ３４ｂを含むシステム１０について説明する。この実装のシステム１０は、クライアントサーバアーキテクチャ又はクラウドベースアーキテクチャに基づき、超並列コンピューティングシステム１２（スタンドアロン又はクラウドベース）として実装されるサーバシステム１２と、クライアントシステム１４とを含む。サーバシステム１２は、メモリ１８と、バスシステム１１と、インターフェイス２０（例えば、ユーザインターフェイス／ネットワークインターフェイス／ディスプレイ又はモニタインターフェイスなど）と、処理装置２４とを含む。メモリ１８内には、メッシュ準備エンジン３２及びシミュレーションエンジン３４が存在する。

【0033】

図１では、メッシュ準備エンジン３２をメモリ１８内に示しているが、メッシュ準備エンジンは、サーバ１２とは異なるシステム上で実行されるサードパーティアプリケーションとすることもできる。メッシュ準備エンジン３２は、メモリ１８内で実行されるか、それともサーバ１２とは異なるシステム上で実行されるかにかかわらず、ユーザ指定のメッシュ定義３０を受け取り、メッシュを準備してシミュレーションエンジン３４に送信（及び／又は記憶）する。シミュレーションエンジン３４は、衝突相互作用モジュール３４ａと、高速流のためのＬＢベースのソルバモジュール３４ｂと、境界モジュール３４ｃと、移流粒子衝突相互作用モジュール３４ｄとを含む。システム１０は、２Ｄ及び／又は３Ｄメッシュ（デカルト座標及び／又は曲線）、座標系及びライブラリを記憶するデータリポジトリ３８にアクセスする。

【0034】

次に、図２に、物理的対象表現の周囲の流体流をシミュレートするプロセス４０を示す。本明細書で説明する例では、物理的対象が翼形部である。しかしながら、物理的対象はあらゆる形状とすることができ、具体的には平面及び／又は（単複の）曲面を有することもできるので、翼形部の使用は例示にすぎない。プロセス４０は、例えばクライアントシステム１４から、シミュレートする物理的対象のメッシュ（又はグリッド）を受け取り、又はデータレポジトリ３８から検索する（４２）。他の実施形態では、外部システム又はサーバ１２のいずれかが、シミュレートする物理的対象のメッシュをユーザ入力に基づいて生成する。このプロセスは、検索したメッシュから幾何学量を事前計算し（４４）、検索したメッシュに対応する事前計算された幾何学量を用いて動的格子ボルツマンモデルシミュレーションを実行する（４６）。格子ボルツマンモデルシミュレーションは、粒子分布の発展（evolution of particle distribution）のシミュレーション４６ａと、エントロピーソルバの評価４６ｂと、ＬＢＭメッシュにおける次のセルへの粒子及びエントロピーの移流（advection）４６ｃとを含む。

【0035】

図３を参照して分かるように、シミュレーションプロセス４６は、例えば高速流３４ｂのためのＬＢベースのソルバモジュール（図１）に適する修正された格子ボルツマン方程式（ＬＢＥ）に従って粒子分布の発展をシミュレートする。プロセス４６（図２を参照）は、衝突演算４６ａ、エントロピーの時間発展のシミュレーション４６ｂの実行後に、ＬＢＭ空間における次のセルへの粒子の移流４６ｃを実行する。

【0036】

エントロピーの時間発展のシミュレーション４６ｂは、近隣のメッシュ位置から流入する分布の組を衝突演算から収集するステップ（５２）と、エントロピーなどの特定のスカラーを表す、流れベクトル上に存在して１つのグリッド点から別のグリッド点に移動するさらなる格子ベクトルの組を収集するステップ（５４）とを含む。エントロピーの時間発展のシミュレーションは、メッシュ内の各位置におけるスカラー値の組をコンピュータによって計算するステップ（５６）も含む。その後、コンピュータは、衝突演算方程式７とさらなる熱源Ｑ_sとの積として流出する分布を求める（５８）。

【0037】

衝突演算子は、流入する格子セットベクトル（ｌａｔｔｉｃｅ ｓｅｔ ｖｅｃｔｏｒｓ）からの高次誤差項と、法線方向衝突演算子から減算（６２）された対応する平均値とをコンピュータが計算（６０）することによって安定する。流体粘性によって流体を加熱して流体伝導によって流体を加熱するさらなる熱源の影響をコンピュータが計算した結果であるさらなる熱源項を計算して（６４）第２の状態に追加する（６６）。コンピュータによってエントロピー拡散を計算し、計算されたエントロピー拡散の結果をさらなる熱源項から除去する（６８）。

【0038】

拡散による熱源の処理
さらなる熱源が存在しなければ、ＬＢＭベースのスカラーソルバは以下のような典型的なスカラー対流拡散方程式になり、

【数14】

（方程式１０）
ここでの値ｋ_sは、τ_qで使用される緩和時間に依存する。（以下に示す方程式６と方程式１０との間の差分に対応する）さらなる熱源項は、エントロピー方程式の正しい挙動を回復させる。

【数15】

（方程式１１）

【0039】

高速流のための強化されたエントロピー衝突演算子
上述した米国特許第９，５７６，０８７号に示されるように、流れ状態のための特殊な衝突演算子について説明する。この演算子は安定の幅が広く、高マッハ用途に適用することができる。

【0040】

しかしながら、高マッハ数の流れでは、エントロピーのための衝突演算子が不安定になる。

【0041】

流体及びエントロピー状態の発展にチャップマン・エンスコッグの展開理論を適用（方程式１及び方程式７）して、ナビエ・ストークス方程式及びエントロピースカラー方程式を回復させる。拡散項は、エントロピーソルバを不安定にし得る衝突演算及び誤差に起因して生じる。

【0042】

この例では、非平衡成分の一次計算（方程式９）を考慮し、テイラー級数展開及びＬＢＭの基本格子要件を使用することによって一次非平衡成分を以下のように近似させることができる。

【数16】

（方程式１２）

【0043】

上記の式では、エントロピーの勾配（ｐδ）に圧力を乗算したものが唯一必要な成分であり、さらなる項（ρｖｖ）は誤差を導入する。項（ρｖｖ）は、Ｍａ²（速度の二乗）に比例するため、特に高速流では不安定性の主因である。

【0044】

この項を削除するために、上記の非平衡計算を以下のように修正する。

【数17】

（方程式１３）

【0045】

上記の非平衡成分は、高次の誤差項を排除して、高速用途にとって安定したエントロピーソルバを提供する。方程式（１１）及び方程式（１３）を方程式（７）に代入すると、以下によって完全な形のエントロピー展開が与えられる。

【数18】

（方程式１４）

【0046】

ＬＢベースのエントロピーソルバ機能
上述したＬＢエントロピーソルバは、高度に正確な過渡的結果と圧縮性の影響とを必要とする高速用途を解決するために有効である。高温度比用途も、このＬＢエントロピーソルバの独自の衝突モデルによって安定性が高まるので、その恩恵を受ける。

【0047】

この節では、特に高速、可変メッシュ分解能、複雑な幾何的形状及び高温度比の流れに関する圧縮性流れのためのＬＢベースのエントロピーソルバの潜在的利点を示す基準研究を提示する。シミュレーション結果を、実験結果及び典型的な有限差分ＦＤベースのエントロピーソルバと比較する。

【0048】

ここで図４Ａ及び図４Ｂに、大気温度よりも２．７０２倍高いジェットノズル温度のマッハ０．９０１のジェットノズル流について行ったシミュレーションの結果を示す。ジェットノズルによって生じる正確なノイズを予測するには、より正確な過渡的結果を得ることが必要である。

【0049】

図４Ａには、ＦＤ法の圧力振動５０ａを示し、図４Ｂには、ＬＢベースの方法の圧力振動５０ｂを示す。

【0050】

ＦＤ法では、ノズル５２ａの出口付近に、可変分解能に起因してジェットノズル５２ａの下流の流れ５６ａに影響する複雑な境界条件を伴う高レベルのノイズ５４ａが生じる。ＬＢベースの方法では、ジェットノズル５２ｂの出口のノイズ５４ｂレベルがきれいに、すなわち低くなっており、ノズルの下流の微細構造５６ｂも維持されている。

【0051】

次に、図５に（デシベルで示す音圧レベル（ＳＰＬ）対Ｈｚで示す周波数の）グラフを示す。このグラフは、ＦＤ法及び（エントロピーソルバを用いた）ＬＢ法のシミュレートした音圧レベルを実験結果と比較したものである。ＦＤ法では、数値的アーチファクトに起因して、特に１ＫＨｚを越える周波数において実験結果よりもはるかに高いノイズレベルが予測されるのに対し、ＬＢベースの方法は、比較的良好に実験結果のものに匹敵する。

【0052】

次に、図６Ａ及び図６Ｂに、音圧レベルマップの予測における、実験結果と比べたＬＢベースのソルバの精度を示す。

【0053】

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

【0054】

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

【0055】

１０１個の速度を含む３Ｄ－２モデル及び３７個の速度を含む２Ｄ－２モデルなどのさらに複雑なモデルを使用することもできる。

【0056】

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

【0057】

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

【0058】

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

【0059】

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

【0060】

格子の分解能は、シミュレートされるシステムのレイノルズ数に基づいて選択することができる。レイノルズ数は、流れの粘度（ｖ）、流れにおける物体の特性長（Ｌ）及び流れの特性速度（ｕ）に関連する。
Ｒｅ＝ｕＬ／ｖ （方程式Ｉ．５）

【0061】

物体の特性長は、物体の大規模な特徴を表す。例えば、マイクロデバイスの周囲の流れをシミュレートする場合には、マイクロデバイスの高さを特性長と見なすことができる。小さな物体領域（例えば、自動車のサイドミラー）の周囲の流れに関心がある時には、シミュレーションの分解能を高めるか、或いは関心領域の周囲に高分解能の領域を使用することができる。格子の分解能が高まるにつれてボクセルの次元は減少する。

【0062】

状態空間はｆ_i（ｘ，ｔ)として表され、ここでのｆ_iは、時点ｔにおいて３次元ベクトルｘによって示される、ある格子サイトにおける状態ｉの単位容積当たりの要素又は粒子の数（すなわち、状態ｉの粒子の密度）を表す。既知の時間増分では、粒子の数が単純にｆ_i（ｘ）として示される。格子サイトの全ての状態の組み合わせは、ｆ（ｘ）として示される。

【0063】

状態の数は、各エネルギーレベル内の考えられる速度ベクトル数によって決まる。速度ベクトルは、ｘ、ｙ及びｚの３次元を有する空間内の整数の線形速度から成る。多重種シミュレーションでは状態の数が増加する。

【0064】

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

【0065】

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

【0066】

３つのエネルギーレベルの考えられる順列を全て生成すると、全部で３９個の可能な状態（１つのエネルギー０状態、６つのエネルギー１状態、８つのエネルギー３状態、６つのエネルギー４状態、１２のエネルギー８状態及び６つのエネルギー９状態）が得られる。

【0067】

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

【0068】

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

【0069】

図１０にマイクロブロックを示す。

【0070】

図１１Ａ及び図１１Ｂも参照すると、表面Ｓ（図１１Ａ）をシミュレーション空間（図１１Ｂ）内にファセットＦαの集合として表している。
Ｓ＝｛Ｆ_α｝ 方程式（Ｉ．７）
ここでのαは、特定のファセットを列挙する指数である。ファセットはボクセル境界に制限されず、比較的少数のボクセルにファセットが影響を与えるように、典型的にはファセットに隣接するボクセルのサイズと同程度又はそれよりもわずかに小さなサイズを有する。ファセットには、表面動力学を実装する目的で特性が割り当てられる。具体的に言えば、各ファセットＦ_αは、単位法線（ｎ_α）と、表面積（Ａ_α）と、中心位置（ｘ_α）と、ファセットの表面動特性を表すファセット分布関数（ｆ_i（α））とを有する。

【0071】

図１２を参照すると、シミュレーション空間の異なる領域内で異なるレベルの分解能を用いて処理効率を高めることができる。通常は、物体６５５の周囲の領域６５０に最も関心があり、従ってこの領域を最高分解能でシミュレートする。粘度の効果は物体からの距離と共に減少するので、減少レベルの分解能（すなわち、拡大されたボクセル容積）を用いて、物体６５５から増加する距離に配置された領域６６０、６６５をシミュレートする。

【0072】

同様に、図１３に示すように、低レベルの分解能を用いて物体７７５のそれほど重要でない特徴の周囲の領域７７０をシミュレートする一方で、最高レベルの分解能を用いて物体７７５の最も重要な特徴（例えば、先端面及び後端面）の周囲の領域７８０をシミュレートすることもできる。中心から遠い領域７８５は、最低レベルの分解能及び最大ボクセルを用いてシミュレートされる。

【0073】

ファセットの影響を受けるボクセルの識別
再び図９を参照すると、シミュレーション空間をモデル化したら（３０２）、１又は２以上のファセットの影響を受けるボクセルを識別する（３０４）。ボクセルは、複数の形でファセットの影響を受けることができる。まず、１又は２以上のファセットが交わるボクセルは、交わっていないボクセルに比べて容量が減少するという点で影響を受ける。この理由は、ファセットと、ファセットが表す表面の下にある材料とがボクセルの一部を占有するからである。分数係数（ｆｒａｃｔｉｏｎａｌ ｆａｃｔｏｒ）Ｐ_f（ｘ）は、ファセットの影響を受けないボクセルの部分（すなわち、流れをシミュレートする対象である流体又は他の材料が占有できる部分）を示す。交わっていないボクセルでは、Ｐ_f（ｘ）が１に等しい。

【0074】

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

【0075】

図１４を参照すると、ファセットＦ_αは、非ゼロの速度ベクトルｃ_iを有する各状態ｉについて、速度ベクトルＣ_iとファセットの単位法線ｎ_αとのベクトルドット積（｜ｃ_iｎ_i｜）の大きさによって定められる高さと、ファセットの表面積Ａ_αによって定められる底面とを有する平行六面体Ｇ_iαによってその容積Ｖ_iαが次式に等しくなるように定められた領域から粒子を受け取り、又はこの領域に粒子を移動させる。
Ｖ_iα＝｜ｃ_iｎ_α｜Ａ_α 方程式（Ｉ．８）

【0076】

ファセットＦ_αは、状態の速度ベクトルがファセットの方に向いている（｜ｃ_iｎ_i｜＜０の）時には容積Ｖ_iαから粒子を受け取り、状態の速度ベクトルがファセットから離れる方に向いている（｜ｃ_iｎ_i｜＞０の）時にはこの領域に粒子を移動させる。後述するように、内角などの非凸面特徴の近傍で起こり得る状態である、別のファセットが平行六面体Ｇ_iαの一部を占有する時には、この式を修正しなければならない。

【0077】

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

【0078】

平行六面体Ｇｉαが重なり合ったボクセルの部分Ｎ（ｘ）は、Ｖ_iα（ｘ）として定められる。この項を使用すると、ボクセルＮ（ｘ）とファセットＦ_αとの間を移動する状態ｉの粒子の流束Γ_iα（ｘ）は、ボクセルの状態ｉの粒子の密度（Ｎ_i（ｘ））と、ボクセルが重なり合った領域の容積（Ｖ_iα（ｘ））とを乗算したものに等しい。
Γ_iα（ｘ）＝Ｎ_i（ｘ）Ｖ_iα（ｘ） 方程式（Ｉ．９）

【0079】

平行六面体Ｇ_iαに１又は２以上のファセットが交わる時には、以下の条件が当てはまり、
Ｖ_iα＝ΣＶ_α（ｘ）＋ΣＶ_iα（β） 方程式（Ｉ．１０）
ここでの第１の加算は、Ｇ_iαが重なり合った全てのボクセルに相当し、第２項は、Ｇ_iαに交わる全てのファセットに相当する。平行六面体Ｇ_iαに別のファセットが交わらない時には、この式が以下のように変形する。
Ｖ_iα＝ΣＶ_iα（ｘ） 方程式（Ｉ．１１）

【0080】

シミュレーションの実行
再び図９を参照すると、１又は２以上のファセットの影響を受けるボクセルが識別されたら（３０４）、タイマを初期化してシミュレーションを開始する（３０６）。シミュレーションの各時間増分中には、粒子と表面ファセットとの相互作用を考慮した移流段階（３０８～３１６）によってボクセルからボクセルへの粒子の移動をシミュレートする。次に、上述したエントロピーソルバの使用に部分的に基づいて、衝突段階（３１８）によって各ボクセル内の粒子の相互作用をシミュレートする。その後、タイマを増分する（３２０）。増分したタイマによってシミュレーションが完了した旨が示されない（３２２）場合には、移流段階と衝突段階と（３０８～３２０）を繰り返す。増分したタイマによってシミュレーションが完了した旨が示された（３２２）場合には、シミュレーションの結果を記憶及び／又は表示する（３２４）。

【0081】

図９に示すように、エントロピー状態の移流には「並列経路」を使用する。この並列経路は、一般に並列動作３０８’～３１８’の後にエントロピーソルバ３４ｂを用いた流れの評価を含む。エントロピー状態は、上述したボクセルからボクセルへの移動を表す３１４と同様に、３次元直線格子３１４’に沿ってボクセル間で移動する。増分されたタイマがシミュレーションの完了（３２２）を示さない限り、このプロセスを繰り返す。

【0082】

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

【0083】

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

【0084】

２．ボクセルからファセットへの収集
粒子と表面との間の相互作用をシミュレートする最初のステップとして、ボクセルから粒子を集めてファセットに提供する（３０８）。上述したように、ボクセルＮ（ｘ）とファセットＦ_αとの間の状態ｉの粒子の流束は以下の通りである。
Γ_iα（ｘ）＝Ｎ_i（ｘ）Ｖ_iα（ｘ） 方程式（Ｉ．１２）

【0085】

この式から、各状態ｉがファセットの方を向いている（ｃ_iｎ_α＜０の）場合、ボクセルによってファセットＦ_αに提供される粒子の数は以下の通りである。
Γ_iαＶ→Ｆ＝ΣｘΓ_iα（ｘ）＝Σ_xＮ_i（ｘ）Ｖ_iα（ｘ） 方程式（Ｉ．１３）

【0086】

Ｖ_iα（ｘ）が非ゼロの値を有するボクセルのみを加算しなければならない。上述したように、ファセットのサイズは、Ｖ_iα（ｘ）が少数のボクセルのみについて非ゼロの値を有するように選択される。Ｖ_iα（ｘ）及びＰ_f（ｘ）は非整数値を有することもできるので、Γ_α（ｘ）は実数として記憶され処理される。

【0087】

３．ファセットからファセットへの移動
次に、粒子がファセット間を移動する（３１０）。ファセットＦαが流入状態（ｃ_iｎ_α＜０）にある平行六面体Ｇ_iαに別のファセットＦ_βが交わる場合、ファセットＦ_αが受け取る状態ｉの粒子の一部はファセットＦ_βから流入するようになる。具体的に言えば、ファセットＦ_αは、前回の時間増分中にファセットＦ_βによって生成された状態ｉの粒子の一部を受け取るようになる。

【0088】

この関係を図１５に示しており、ここでは、ファセットＦ_βが交わる平行六面体Ｇ_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α
方程式（Ｉ．１６）

【0089】

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

【0090】

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

【0091】

平行状態（ｃ_iｎ_α＝０）では、Ｖ_iα及びＶ_iα（ｘ）がいずれもゼロである。Ｎ_i（α）を求める式では、Ｖ_iα（ｘ）が（Γ_iOTHER（α）の式からの）分子に現れ、Ｖ_iαが（Ｎ_i（α）の式からの）分母に現れる。従って、Ｖ_iα及びＶ_iα（ｘ）がゼロに近付くと、平行状態のＮ_i（α）は、Ｎ_i（α）の極限として決定される。ゼロ速度を有する状態の値（すなわち休止状態、並びに状態（０，０，０，２）及び（０，０，０，－２））は、シミュレーションの最初に温度及び圧力の初期条件に基づいて初期化される。その後、これらの値は時間をかけて調整される。

【0092】

４．ファセット表面動力学の実行
次に、各ファセットが上述した４つの境界条件を満たすように表面動力学を実行する（図９、３１２）。

【0093】

図１６に、ファセットの表面動力学の実行手順を示す。最初に、ファセットにおける粒子の合計運動量Ｐ（α）を以下のように求めることによって、ファセットＦ_αに垂直な合計運動量を求める（１１０５）。
全てのｉについて、

【数19】

方程式（Ｉ．２０）
この式から、以下のように法線運動量Ｐ_n（α）を求める。
Ｐ_n（α）＝ｎ_α・Ｐ（α）方程式（Ｉ．２１）

【0094】

次に、プッシング／プリング法を用いてこの法線運動量を排除（１１１０）してＮ_α-（α）を取り出す。この手法によれば、粒子は、法線運動量にのみ影響を与える形で状態間を移行する。このプッシング／プリング法は、引用により組み入れられる米国特許第５，５９４，６７１号に記載されている。

【0095】

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

【0096】

次に、流入流束分布とボルツマン分布とに基づいてファセットＦ_αの流出流束分布を求める（１１２０）。まず、流入流束分布Γ_i（α）とボルツマン分布との差分を以下のように求める。
ΔΓ_i（α）＝Γ_iIN（α）－Ｎ_n-βi（α）Ｖ_iα 方程式（Ｉ．２２）

【0097】

この差分を使用すると、流出流束分布は、
ｎ_αｃ_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αは、ｎ_αとｔ_iαの両方に垂直な第２の接線ベクトルであり、ΔＮ_j,1及びΔＮ_j,2は、状態ｉのエネルギー（ｊ）及び指示される接線ベクトルに対応する分布関数である。分布関数は、次式に従って決定され、

【数20】

方程式（Ｉ．２５）
ここでのｊは、エネルギーレベル１状態では１に等しく、エネルギーレベル２状態では２に等しい。

【0098】

Γ_iOUT（α）の方程式の各項の関数は以下の通りである。第１項及び第２項は、ボルツマン分布を生じる上で衝突が効果的である限りは法線運動量流束境界条件を強制するが、接線運動量流束の例外を含む。第４項及び第５項は、この不十分な衝突に起因する離散性効果又は非ボルツマン構造を原因として生じ得る例外を補正する。最後に、第３項は、表面上の接線運動量流束の所望の変化を強制するように一定量の表面摩擦を追加する。摩擦係数Ｃ_fの生成については以下で説明する。なお、ベクトル操作に関与する全ての項は、シミュレーションの開始前に計算できる幾何学的因子である。

【0099】

この式から、接線速度が以下のように求められ、
ｕ_i（α）＝（Ｐ（α）－Ｐ_n（α）ｎ_α）／ρ 方程式（Ｉ．２６）
ここでのｐは、以下のようなファセット分布の密度である。

【数21】

方程式（Ｉ．２７）

【0100】

上記と同様に、流入流束分布とボルツマン分布との差分を以下のように求める。
ΔΓ_i（α）＝Γ_iIN（α）－Ｎ_n-βi（α）Ｖ_iα 方程式（Ｉ．２８）

【0101】

すると、流出流束分布は以下のようになり、
Γ_iOUT（α）＝Ｎ_n-βi（α）Ｖ_iα－ΔΓ_i*（α）＋Ｃ_f（ｎ_αｃ_i）［Ｎ_n-βi*（α）－Ｎ_n-βi（α）］Ｖ_iα
方程式（Ｉ．２９）
この式は、先の手法によって求めた流出流束分布の最初の２行に対応するが、異常な接線流束のための補正を必要としない。

【0102】

いずれかの方法を使用すると、結果として得られる流束分布は運動量流束条件を全て満たし、すなわち、

【数22】

方程式（Ｉ．３０）
となり、ここでのｐ_αは、ファセットＦαにおける平衡圧であって、ファセットに粒子を提供するボクセルの平均密度及び温度値に基づき、ｕ_αは、ファセットにおける平均速度である。

【0103】

質量及びエネルギー境界条件が確実に満たされるように、各エネルギーレベルｊについて入力エネルギーと出力エネルギーとの差分を以下のように測定する。

【数23】

方程式（Ｉ．３１）
ここでの添字ｊは、状態ｉのエネルギーを示す。次に、このエネルギー差を用いて差分項（ｄｉｆｆｅｒｅｎｃｅ ｔｅｒｍ）を生成すると、
ｃ_jiｎ_α＞０の場合、以下のようになる。

【数24】

方程式（Ｉ．３２）
この差分項を用いて流出流束を修正すると、流速は、
ｃ_jiｎ_α＞０の場合、以下のようになる。
Γ_αjiOUTf＝Γ_αjiOUT＋δΓ_αji 方程式（Ｉ．３３）
この演算は、接線運動量流束をそのままにした状態で質量及びエネルギー流束を訂正する。ファセットの近傍において流れがおおよそ均一であって平衡に近い場合には、この調整はわずかなものである。結果として得られる調整後の法線運動量流束は、近傍の平均特性に近傍の非均一性又は非平衡特性に起因する補正をプラスしたものに基づく平衡圧である値へとわずかに変化する。

【0104】

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

【0105】

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

【0106】

表面と相互作用しないボクセルの場合、この移動操作は計算的に極めて単純である。毎時間増分中に、状態の集団全体をその現在のボクセルから宛先ボクセルに移動させる。同時に、宛先ボクセルの粒子を、そのボクセルから独自の宛先ボクセルに移動させる。例えば、＋１ｘ及び＋１ｙ方向に移動しているエネルギーレベル１の粒子（１，０，０）は、その現在のボクセルから、ｘ方向に＋１上回り他の方向が０であるボクセルに移動させる。この粒子は、最後には移動前の状態と同じ状態の宛先ボクセル（１，０，０）に行き着く。ボクセル内の相互作用は、他の粒子及び表面との局所的相互作用に基づいてその状態の粒子総数を変化させることがある。変化しなければ、粒子は格子に沿って同じ速度で同じ方向に移動し続ける。

【0107】

この移動操作は、１又は２以上の表面と相互作用するボクセルでは若干複雑になる。この結果、１又は２以上の端数粒子（ｆｒａｃｔｉｏｎａｌ ｐａｒｔｉｃｌｅｓ）がファセットに移動することがある。このような端数粒子のファセットへの移動により、ボクセルに端数粒子が留まるようになる。これらの端数粒子を、ファセットに占有されたボクセルに移動させる。

【0108】

図１７を参照すると、ボクセル９０５の状態ｉの粒子の一部９００がファセット９１０に移動した時には（３０８）、残り部分９１５をファセット９１０が位置するボクセル９２０に移動させ、ここから状態ｉの粒子をファセット９１０に向ける。従って、状態集団が２５に等しく、Ｖ_iα（ｘ）が０．２５に等しい（すなわち、ボクセルの４分の１が平行六面体Ｇ_iαに交わる）場合には、６．２５の粒子がファセットＦ_αに移動し、１８．７５の粒子を、ファセットＦ_αに占有されたボクセルに移動させる。１つのボクセルに複数のファセットが交わることもあるので、１又は２以上のファセットに占有されたボクセルに移動する状態ｉの粒子の数Ｎ（ｆ）は以下のようになり、

【数25】

方程式（Ｉ．３４）
ここでのＮ（ｘ）はソースボクセルである。

【0109】

６．ファセットからボクセルへの散乱
次に、各ファセットからの流出粒子をボクセルに散乱させる（３１６）。基本的に、このステップは、ボクセルからファセットに粒子が移動する収集ステップの逆である。ファセットＦαからボクセルＮ（ｘ）に移動する状態ｉの粒子の数は以下のようになり、

【数26】

方程式（Ｉ．３５）
ここでのＰ_f（ｘ）は、部分的なボクセルの容積縮小に相当する。この式から、各状態ｉについて、ファセットからボクセルに向かう粒子の総数Ｎ（ｘ）は以下のようになる。

【数27】

方程式（Ｉ．３６）

【0110】

粒子をファセットからボクセルに散乱させ、これらの粒子を周囲のボクセルから移流した粒子と組み合わせて結果を整数化した後には、一定のボクセルの一定の方向がアンダーフロー（負になる）又はオーバーフロー（８ビット実装において２５５を超える）のいずれかになることがある。この結果、これらの量を許容できる値の範囲に収まるように切り捨てた後に、質量、運動量及びエネルギーが増加又は減少するようになる。このような発生を防ぐために、有害な状態を切り捨てる前に、限界を超えた質量、運動量及びエネルギーを蓄積する。状態を有するエネルギーでは、（アンダーフローに起因して）得られた又は（オーバーフローに起因して）失われた値に等しい質量が、同じエネルギーを有するランダムに（又は順次に）選択された状態に再び追加され、これ自体がオーバーフロー又はアンダーフローに影響されることはない。この質量及びエネルギーの追加によって生じるさらなる運動量を蓄積して、切り捨てからの運動量に追加する。同じエネルギー状態のみに質量を追加することにより、質量カウンタがゼロに達した時に質量及びエネルギーがいずれも補正される。最後に、運動量積算器がゼロに戻るまでプッシング／プリング法を用いて運動量を補正する。

【0111】

７．流体力学の実行
流体力学を実行する（図９、３１８）。このステップは、マイクロ動力学操作又はボクセル内操作と呼ぶことができる。同様に、移流手順は、ボクセル間操作と呼ぶことができる。後述するマイクロ動力学操作を用いて、ファセットにおける粒子を衝突させてボルツマン分布を生成することもできる。

【0112】

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

【0113】

図１８を参照すると、システムは、エントロピーの決定が衝突演算から分離された粒子輸送の分布を生成する（１２００）。この例では、分布関数が移流を含み、エントロピー決定を用いて粒子衝突を増強するのではなく、エントロピー決定を用いて移流を増強することにより、エントロピー決定が分布部分に含まれる。一般に、移流は、粒子の（例えば、１つの領域から別の領域への水平方向の）輸送を含む。

【0114】

動作中、システムは、格子速度の組において、一定量の流体中の粒子間の衝突を引き起こす粒子の輸送をシミュレートする（１２０２）。システムは、粒子の輸送の分布関数に基づいて粒子分布も生成し（１２０４）、この分布関数は、熱力学ステップ及び粒子衝突ステップを含み、熱力学ステップは、粒子衝突ステップから実質的に独立して分離する。

【0115】

図１２の例では、システムが、特定の時点ｔにおける流体容積中の特定の位置ｘの衝突についての衝突後分布関数ｆ_i’（ｘ，ｔ）も求め（１２０４）、ｆ_i’（ｘ，ｔ）＝ｆ_i（ｘ，ｔ）＋Ｃ_i（ｘ，ｔ）であり、ここでのＣｉは衝突演算子であり、ｆ_iは衝突前の粒子の分布関数である。

【0116】

システムは、上記の方程式７に従って、さらなる格子ベクトルの組ｑ_iを用いて特定のエントロピーｓを表すことによってエントロピー及び輸送を解決する（１２０６）。方程式７及び方程式１は、温度などのスカラーの移流について米国特許出願公開第２０１６／０１８８７６８号に記載されているように、１つのメッシュ点から別のメッシュ点にエントロピーを伝えて、流体流の移流と共に移流するスカラーとして効果的にエントロピーを提供する。この方法は、異なるメッシュレベル付近及び複雑な境界付近の滑らかな遷移をもたらす。

【0117】

可変分解能
（米国特許出願公開第２０１３／０１５１２２１号に記載されるような）可変分解能を使用することもでき、この可変分解能は、例えば粗ボクセル及び微小ボクセルなどの異なるサイズのボクセルを使用する。

【0118】

図１９に、流体シミュレーションのスクリーンショットを示す。上述したエントロピーソルバ法を用いた流体シミュレーションは、同様の流体シミュレーション表現と、いずれかの慣習的な対応する計算データ出力とをもたらす。しかしながら、上述したエントロピーソルバ法を用いたこのような流体シミュレーションは、例えば実際の物理的対象などの曲面を有する物体をモデル化する際に、他の方法よりも高速で、及び／又は少ない計算リソースを用いて流体シミュレーションを実行できるとともに、高度に正確な過渡的結果と圧縮性の影響とを必要とする用途の場合に比較的安定したエントロピーソルバを維持しながら高マッハ流体流における物体をモデル化するために使用することもできる。

【0119】

ハードウェア及びソフトウェア実装
システム１０は、一般に様々な異なる電気的又は電子的計算又は処理装置のうちのいずれか１つとして実装することができ、上述した様々なステップのあらゆる組み合わせを実行して、開示する熱管理システムの様々なコンポーネントを制御することができる。

【0120】

システム１０は、一般に及び任意にプロセッサ（又は複数のプロセッサ）、メモリ、記憶装置及び入力／出力装置のうちのいずれか１つ又は２つ以上を含むことができる。これらのコンポーネントの一部又は全部は、図１に示すようなシステムバスを用いて相互接続することができる。プロセッサは、実行命令を処理することができる。いくつかの実施形態では、プロセッサがシングルスレッドプロセッサである。特定の実施形態では、プロセッサがマルチスレッドプロセッサである。通常、プロセッサは、メモリ又は記憶装置に記憶された命令を処理して、ユーザインターフェイスのためのグラフィック情報を入力／出力装置上に表示し、上述した様々なモニタリング及び制御機能を実行することができる。本明細書に開示するシステムに適したプロセッサは、汎用及び専用マイクロプロセッサ、並びに唯一のプロセッサ、或いはあらゆる種類のコンピュータ又はコンピュータ装置の複数のプロセッサのうちの１つを含む。

【0121】

メモリは、システム内の情報を記憶し、揮発性又は不揮発性メモリなどのコンピュータ可読媒体とすることができる。記憶装置は、システム１０のための大容量ストレージを提供することができる。一般に、記憶装置は、コンピュータ可読命令を記憶するように構成されたいずれかの非一時的有形媒体を含むことができる。例えば、記憶装置は、内蔵ハードディスク及びリムーバブルディスクなどの磁気ディスクと、光磁気ディスクと、光ディスクとを含むコンピュータ可読媒体及び関連するコンポーネントを含むことができる。コンピュータプログラム命令及びデータを有形的に具体化するのに適した記憶装置は、一例としてＥＰＲＯＭ、ＥＥＰＲＯＭ及びフラッシュメモリデバイスなどの半導体メモリデバイスと、内蔵ハードディスク及びリムーバブルディスクなどの磁気ディスクと、光磁気ディスクと、ＣＤ－ＲＯＭ及びＤＶＤ－ＲＯＭディスクとを含む全ての形態の不揮発性メモリを含む。本明細書に開示するシステムのプロセッサ及びメモリユニットは、ＡＳＩＣ（特定用途向け集積回路）によって補完し、又はＡＳＩＣに組み込むことができる。

【0122】

入力／出力装置は、システム１０のための入力／出力動作を行い、キーボード及び／又はポインティングデバイスを含むことができる。いくつかの実施形態では、入力／出力装置が、グラフィカルユーザインターフェース及びシステム関連情報を表示するためのディスプレイユニットを含む。

【0123】

様々な測定機能、モニタリング機能、制御機能及び通信機能を実行するためのコンポーネントを含む本明細書で説明した特徴は、デジタル電子回路に、或いはコンピュータハードウェア、ファームウェア又はこれらの組み合わせに実装することができる。方法ステップは、例えば機械可読記憶装置などの情報担体内で有形的に具体化された、プログラマブルプロセッサ（例えば、システム１０の）によって実行されるコンピュータプログラム製品に実装することができ、プログラマブルプロセッサがこのようなプログラム命令を実行して上述したステップ及び機能のいずれかを実行することによって機能を実行することができる。１又は２以上のシステムプロセッサによる実行に適したコンピュータプログラムは、命令を実行するプロセッサ又は他のコンピュータ装置に上述した様々なステップを含む特定の動作を実行させる、直接的又は間接的に使用できる命令セットを含む。

【0124】

本明細書に開示したシステム及び方法と共に使用するのに適したコンピュータプログラムは、コンパイラ型言語又はインタープリタ型言語を含むあらゆる形態のプログラミング言語で書くことができ、スタンドアロンププログラムとしての形態、或いはモジュール、コンポーネント、サブルーチン又はコンピュータ環境での使用に適した他のユニットとしての形態を含むあらゆる形態で展開することができる。

【0125】

複数の実装について説明した。それでもなお、本開示の趣旨及び範囲から逸脱することなく様々な修正を行うことができると理解されるであろう。従って、以下の特許請求の範囲には他の実装も含まれる。

【図1】

【図2】

【図3】

【図4A】

【図4B】

【図5】

【図6A】

【図6B】

【図7】

【図8】

【図9】

【図10】

【図11A】

【図11B】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】