JP6562307B2 - 層流乱流遷移のモデル化を含む物理的プロセスのコンピュータシミュレーション - Google Patents

Description

本明細書は、流体流及び音響などの物理的プロセスのコンピュータシミュレーションに関する。本明細書は、境界層における層流乱流遷移現象を予測する方法にも関する。

高レイノルズ数の流れは、巨視的物理量（例えば、密度、温度、流速）を表す変数に対して多くの離散的空間位置の各々において高精度浮動小数点算術演算を行うことによってナビエ・ストークス（Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ）微分方程式の離散解を生成することでシミュレートされてきた。別の手法では、微分方程式の代わりに、一般に格子ガス（又はセル）オートマンとして知られているものを導入し、この場合、ナビエ・ストークスの方程式を解くことによってもたらされる巨視的レベルのシミュレーションが、格子サイト間を移動する粒子に関する演算を行う微視的レベルのモデルに置き換えられる。

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

Ｌａｎｇｔｒｙ他著、ＡＳＭＥターボＥｘｐｏ２００４会報、題名「局所的可変を用いた相関に基づく遷移モデル（Ａ Ｃｏｒｒｅｌａｔｉｏｎ−Ｂａｓｅｄ Ｔｒａｎｓｉｔｉｏｎ Ｍｏｄｅｌ Ｕｓｉｎｇ Ｌｏｃａｌ Ｖａｒｉａｂｌｅｓ）」

本明細書は、境界層における層流乱流遷移現象を予測する方法にも関する。本明細書は、予測された層流乱流遷移に基づいて表面上の領域（例えば、ファセット、サーフェル）に適した壁剪断応力値を選択し、この選択された壁剪断応力値を流体流のシミュレーションに使用する方法にも関する。

一般に、本文書では、格子ボルツマン（ＬＢ）法を用いて流体流をシミュレートする技術、及びスカラー輸送方程式を解く技術について説明する。本明細書で説明する手法では、層流乱流境界層遷移を含む流体流をコンピュータ上でシミュレートする方法が、層流境界層の流れのための第１の計算を行うステップと、乱流境界層の流れのための第２の計算を行うステップとを含む。この方法は、第１及び第２の境界層計算の少なくとも一方の結果を基準と比較するステップと、この比較の結果に基づいて、表面及びこの表面近傍の流体の少なくとも一方を表す複数の要素の少なくともいくつかについて、層流境界層の流れのための第１の計算の結果又は乱流境界層の流れのための第２の計算の結果を選択するステップと、ある体積の流体の活動を、この体積内の要素の動きをモデル化するようにシミュレーションを行うステップとを含み、このシミュレーションは、複数の要素のための選択された結果に部分的に基づく。

実施形態は、以下のうちの１つ又はそれ以上を含むことができる。

層流境界層の流れのための第１の計算を行うステップは、層流の壁運動量フラックステンソル特性を計算するステップを含むことができ、乱流境界層の流れのための第２の計算を行うステップは、乱流の壁運動量フラックステンソル特性を計算するステップを含むことができ、複数の要素の少なくともいくつかについて、第１の境界層計算の結果又は乱流境界層の流れのための第２の計算の結果を選択するステップは、層流壁運動量フラックステンソル特性又は乱流壁運動量フラックステンソル特性を選択するステップを含むことができる。

境界層の層流乱流遷移を求めるステップは、表面上の複数のファセットの各々について、第１の境界層計算に基づく第１の測定値と、第２の境界層計算に基づく第２の測定値とを求めるステップと、第１及び第２の測定値の少なくとも一方を基準と比較することにより、複数のファセットの少なくともいくつかについて、流れを層流又は乱流として分類するステップとを含むことができる。

表面上の複数のファセットの少なくともいくつかについて、層流境界層の流れのための第１の計算の結果又は乱流境界層の流れのための第２の計算の結果を選択するステップは、層流として分類されたファセットについて、層流の壁運動量フラックステンソル特性を選択するステップと、乱流として分類されたファセットについて、乱流の壁運動量フラックステンソル特性を選択するステップとを含むことができる。

第１の境界層計算の結果は、層流壁運動量フラックステンソルの測定値を含むことができ、第２の境界層計算の結果は、乱流壁運動量フラックステンソルの測定値を含むことができ、比較は、乱流強度の測定値を含むことができる。

第１の境界層計算を行うステップは、表面上の複数のファセットの各々について、層流壁運動量フラックステンソルの測定値を計算するステップを含むことができ、第２の境界層計算を行うステップは、表面上の複数のファセットの各々について、乱流壁運動量フラックステンソルの測定値を計算するステップを含むことができる。第１及び第２の境界層計算の少なくとも一方の結果を基準と比較するステップは、表面上の複数のファセットの各々について、計算された乱流強度の測定値と、乱流壁運動量フラックステンソルの測定値とを比較するステップを含むことができ、第１の境界層計算の結果又は第２の境界層の結果を選択するステップは、表面上の複数のファセット少なくともいくつかについて、乱流強度の測定値と乱流壁運動量フラックステンソルの測定値との比較に基づいて、計算された乱流壁運動量フラックステンソル特性及び層流壁運動量フラックステンソル特性の一方を選択するステップを含むことができる。

表面上の複数のファセットの各々について、乱流強度の測定値と乱流壁運動量フラックステンソルの測定値とを比較するステップは、乱流強度の測定値が壁運動量フラックステンソルの測定値よりも大きいかどうかを決定するステップを含むことができ、表面上の複数のファセット少なくともいくつかについて、計算された乱流壁運動量フラックステンソル及び層流壁運動量フラックステンソルの一方を選択するステップは、特定のファセットについて、乱流強度の測定値が乱流壁運動量フラックステンソルの測定値よりも大きい場合には、乱流壁運動量フラックステンソルを選択し、乱流強度の測定値が乱流壁運動量フラックステンソルの測定値よりも小さい場合には、層流壁運動量フラックステンソルの測定値を選択するステップを含むことができる。

局所的乱流強度の測定値を計算するステップは、局所的乱流運動エネルギーの値を計算するステップを含むことができる。

所与の壁近傍の流体速度について、乱流壁運動量フラックステンソルの測定値は、層流壁運動量フラックステンソルの測定値よりも大きいものとすることができる。

ある体積の流体の活動のシミュレーションを行うステップは、状態ベクトルに対して、異なる運動量状態の要素間の相互作用をモデルに従ってモデル化する相互作用演算を行うステップと、状態ベクトルの組の第１の移動演算を行って、体積内の新たなボクセルへの要素の移動をモデルに従って反映するステップとを含むことができる。

第２の境界層計算は、速度プロファイル及び壁からの距離に基づいて乱流壁運動量フラックステンソルの測定値を求めるための計算を含むことができる。

この方法は、表面上の複数のファセットの少なくともいくつかについて、層流境界層の流れのための第１の計算の結果、及び乱流境界層の流れのための第２の計算の結果の組み合わせに基づく値を選択するステップを含むこともできる。

この方法は、表面上の複数のファセットの少なくともいくつかについて、乱流壁運動量フラックステンソル特性と層流壁運動量フラックステンソル特性の組み合わせに基づく壁運動量フラックステンソル特性を選択するステップを含むこともできる。

第２の境界層計算は、局所的乱流運動エネルギー及び局所的流体速度に基づいて乱流壁運動量フラックステンソルの測定値を求めるための計算を含むことができる。

表面に隣接する領域のボクセルサイズは、表面から間隔を空けた領域のボクセルサイズに類似することができる。

表面に隣接する領域のボクセルサイズは、表面から間隔を空けた領域のボクセルサイズと同一とすることができる。

いくつかの態様では、コンピュータ可読媒体内で有形的に具体化されるコンピュータプログラム製品が、実行時に層流乱流境界層遷移を含む物理的プロセスの流体流のシミュレーションを行う命令を含むことができる。このコンピュータプログラム製品は、コンピュータに、層流境界層の流れのための第１の計算を行わせ、乱流境界層の流れのための第２の計算を行わせ、第１及び第２の境界層計算の少なくとも一方の結果を基準と比較させ、この比較の結果に基づいて、表面及びこの表面近傍の流体の少なくとも一方を表す複数の要素の少なくともいくつかについて、層流境界層の流れのための第１の計算の結果又は乱流境界層の流れのための第２の計算の結果を選択させ、ある体積の流体の活動を、この体積内の要素の動きをモデル化するようにシミュレーションを行わせるように構成することができ、このシミュレーションは、複数の要素のための選択された結果に部分的に基づく。

実施形態は、以下のうちの１つ又はそれ以上を含むことができる。

層流境界層の流れのための第１の計算を行うための命令は、層流の壁運動量フラックステンソル特性を計算するための命令を含むことができ、乱流境界層の流れのための第２の計算を行うための命令は、乱流の壁運動量フラックステンソル特性を計算するための命令を含むことができ、第１の境界層計算の結果又は乱流境界層の流れのための第２の計算の結果を選択するための命令は、層流壁運動量フラックステンソル特性又は乱流壁運動量フラックステンソル特性を選択するための命令を含むことができる。

境界層の層流乱流遷移を求めるための命令は、表面上の複数のファセットの各々について、第１の境界層計算に基づく第１の測定値と、第２の境界層計算に基づく第２の測定値とを求めるための命令と、第１及び第２の測定値の少なくとも一方を基準と比較することにより、複数のファセットの少なくともいくつかについて、流れを層流又は乱流として分類するための命令とを含むことができる。

表面上の複数のファセットの少なくともいくつかについて、層流境界層の流れのための第１の計算の結果又は乱流境界層の流れのための第２の計算の結果を選択するための命令は、層流として分類されたファセットについて、層流の壁運動量フラックステンソル特性を選択するための命令と、乱流として分類されたファセットについて、乱流の壁運動量フラックステンソル特性を選択するための命令とを含むことができる。

第１の境界層計算の結果は、層流壁運動量フラックステンソル特性の測定値を含むことができ、第２の境界層計算の結果は、乱流壁運動量フラックステンソル特性の測定値を含むことができ、基準は、乱流強度の測定値を含むことができる。

第１の境界層計算を行うための命令は、表面上の複数のファセットの各々について、層流壁運動量フラックステンソルの測定値を計算するための命令を含むことができ、第２の境界層計算を行うための命令は、表面上の複数のファセットの各々について、乱流壁運動量フラックステンソルの測定値を計算するための命令を含むことができ、第１及び第２の境界層計算の少なくとも一方の結果を基準と比較するための命令は、表面上の複数のファセットの各々について、計算された乱流強度の測定値と、乱流壁運動量フラックステンソルの測定値とを比較するための命令を含むことができ、第１の境界層計算の結果又は第２の境界層の結果を選択するための命令は、表面上の複数のファセット少なくともいくつかについて、乱流強度の測定値と乱流壁運動量フラックステンソルの測定値との比較に基づいて、計算された乱流壁運動量フラックステンソル特性及び層流壁運動量フラックステンソル特性の一方を選択するための命令を含むことができる。

いくつかのさらなる態様によれば、物理的プロセスの流体流のシミュレーションを行うためのコンピュータシステムを、層流境界層の流れのための第１の計算を行い、乱流境界層の流れのための第２の計算を行い、第１及び第２の境界層計算の少なくとも一方の結果を基準と比較し、この比較の結果に基づいて、表面及びこの表面近傍の流体の少なくとも一方を表す複数の要素の少なくともいくつかについて、層流境界層の流れのための第１の計算の結果又は乱流境界層の流れのための第２の計算の結果を選択し、ある体積の流体の活動を、該体積内の要素の動きをモデル化するようにシミュレーションを行うように構成することができ、シミュレーションは、複数の要素のための選択された結果に部分的に基づく。

実施形態は、以下のうちの１つ又はそれ以上を含むことができる。

層流境界層の流れのための第１の計算を行うための構成は、層流の壁運動量フラックステンソル特性を計算するための構成を含むことができる。乱流境界層の流れのための第２の計算を行うための構成は、乱流の壁運動量フラックステンソル特性を計算するための構成を含むことができる。第１の境界層計算の結果又は乱流境界層の流れのための第２の計算の結果を選択するための構成は、層流壁運動量フラックステンソル特性又は乱流壁運動量フラックステンソル特性を選択するための構成を含むことができる。

境界層の層流乱流遷移を求めるための構成は、表面上の複数のファセットの各々について、第１の境界層計算に基づく第１の測定値と、第２の境界層計算に基づく第２の測定値とを求めるための構成と、第１及び第２の測定値の少なくとも一方を基準と比較することにより、複数のファセットの少なくともいくつかについて、流れを層流又は乱流として分類するための構成とを含むことができる。

表面上の複数のファセットの少なくともいくつかについて、層流境界層の流れのための第１の計算の結果又は乱流境界層の流れのための第２の計算の結果を選択するための構成は、層流として分類されたファセットについて、層流の壁運動量フラックステンソル特性値を選択するための命令と、乱流として分類されたファセットについて、乱流の壁運動量フラックステンソル特性を選択するための命令とを含むことができる。

第１の境界層計算の結果は、層流壁運動量フラックステンソルの測定値を含むことができ、第２の境界層計算の結果は、壁運動量フラックステンソルの測定値を含むことができ、比較は、乱流強度の測定値を含む。

第１の境界層計算を行うための構成は、表面上の複数のファセットの各々について、層流壁運動量フラックステンソルの測定値を計算するための構成を含み、第２の境界層計算を行うための構成は、表面上の複数のファセットの各々について、乱流壁運動量フラックステンソルの測定値を計算するための構成を含むことができる。第１及び第２の境界層計算の少なくとも一方の結果を基準と比較するための構成は、表面上の複数のファセットの各々について、計算された乱流強度の測定値と、乱流壁運動量フラックステンソルの測定値とを比較するための構成を含むことができる。第１の境界層計算の結果又は第２の境界層の結果を選択するための構成は、表面上の複数のファセット少なくともいくつかについて、乱流強度の測定値と乱流壁運動量フラックステンソルの測定値との比較に基づいて、計算された乱流壁運動量フラックステンソル特性及び層流壁運動量フラックステンソル特性の一方を選択するための構成を含むことができる。

ボルツマンレベルのメゾスコピック表現
統計物理学の分野では、流体系をいわゆる「メゾスコピック」レベルの運動方程式で表せることが周知である。このレベルでは、個々の粒子の詳細な動きを特定する必要はない。その代わりに、流体の特性は、単一の粒子位相空間を用いて定義される粒子分布関数

によって表され、式中、ｘは空間座標、νは粒子速度座標である。質量、密度、流体速度及び温度などの典型的な流体力学的物理量は、粒子分布関数の単純なモーメントである。粒子分布関数の動特性は、以下のボルツマン方程式に従う。

式中、Ｆ（ｘ、ｔ）は、（ｘ、ｔ）における外部無撞着的又は自己無撞着的に生成される体積力を表す。衝突項Ｃは、様々な速度及び位置にある粒子の相互作用を表す。上記のボルツマン方程式は、衝突項Ｃについての特定形態を規定することなく、（ボルツマンによって当初構築されたような）希薄ガスの周知の状況のみならず全ての流体システムに適用可能であることを強調しておくことが重要である。

一般的に、Ｃは、二点相関関数の複雑な多次元積分を含む。分布関数ｆの閉鎖系を形成するためだけでなく、効率的計算のためにも最も好都合で物理的に一貫性のある形態の１つは、周知のＢＧＫ演算子である。ＢＧＫ演算子は、衝突のどのような小さな要素であっても、分布関数は、衝突：

を通じて

によって与えられる明確に定義された局所平衡に近づくという物理的根拠に従って構築される。
式中、パラメータτは、衝突による平衡までの特性緩和時間を表す。粒子（例えば、原子又は分子）の処理では、この緩和時間は通常一定と見なされる。「ハイブリッド」（流体−運動）表現では、この緩和時間は、歪み速度、乱流運動エネルギー及びその他のような流体力学変数の関数である。従って、乱流は、局所的に決定された特性を有する乱流粒子（「渦」）のガスとして表すことができる。

ボルツマン−ＢＧＫ方程式の数値解法は、ナビエ・ストークスの方程式の解法よりも優れた計算上の利点をいくつか有する。第一に、方程式内に複雑な非線形項又は高次空間導関数が存在せず、従って、移流不安定性に関する問題がほとんどないと容易に認識することができる。この記述レベルでは、圧力を扱う必要がないので方程式が局所的であり、このことがアルゴリズムの並列化において大きな利点をもたらす。線形の移流演算子であるという別の望ましい特徴は、二次空間導関数を有する拡散演算子が存在しないことと共に、流体偏微分方程式（「ＰＤＥ」）の数学的条件ではなく、実際の固体面と粒子が実際にどのように相互作用するかを模倣する形で非滑り面又は滑り面などの物理的境界条件を容易に実現できることである。直接的な利点の１つは、固体面上における界面の動きを扱う問題がなく、これにより、格子ボルツマンベースのシミュレーションソフトウェアが複雑乱流空気力学を良好にシミュレートできるようになる点である。また、有限粗面などの境界による特定の物理的特性を有効に組み込むこともできる。さらに、ＢＧＫ衝突演算子は純粋に局所的であるのに対し、自己無撞着性の体積力の算術演算は近接情報のみを介して達成することができる。従って、ボルツマン−ＢＧＫ方程式の演算は、並列処理に効果的に適合することができる。

格子ボルツマン定式化
連続体ボルツマン方程式を解くことは、位置及び速度の位相空間において積分微分方程式の数値的評価を伴う点で重要な課題となる。位置だけでなく速度の位相空間も離散化できることが観測された時に行われる大幅な簡素化の結果、ボルツマン方程式の解法の効率的な数値アルゴリズムがもたらされた。流体力学の物理量は、高々最近傍情報に依存する単純和の項で記述することができる。従来、格子ボルツマン方程式の定式化は、速度

の離散集合上での粒子の発展を規定する格子ガスモデルに基づくものであったが、この方程式は、連続体ボルツマン方程式の離散化としての第一原理から体系的に導くことができる。この結果、ＬＢＥは、格子ガス手法に関連する周知の問題の影響を受けない。従って、位相空間における連続分布関数ｆ（ｘ，ｖ，ｔ）を扱うのではなく、離散速度指数を表記した添字を有する離散分布の有限集合

を辿ることのみが必要とされる。巨視的記述の代わりにこの運動方程式を扱うことの主な利点は、システムの位相空間の増加が問題の局所性によって相殺される点である。

対称性を考慮することにより、速度値の集合は、位置座標空間にわたる時に特定の格子構造を形成するように選択される。このような離散値系の動特性は、

の形を有するＬＢＥに従い、通常、衝突演算子は上述したＢＧＫ形式をとる。平衡分布形式を適切に選択することにより、格子ボルツマン方程式が正しい流体力学特性及び熱流体力学特性をもたらすことを理論的に説明することができる。すなわち、

から導かれる流体力学的モーメントは、巨視的極限においてナビエ・ストークスの方程式に従う。これらのモーメントは、以下のように定義される。

式中、ρ、ｕ、及びＴは、それぞれ流体密度、速度及び温度であり、Ｄは、（物理空間次元とは全く等しくない）離散速度空間の次元である。

その他の特徴及び利点は、図面及び特許請求の範囲を含めて、以下の詳細な説明から明らかになるであろう。

ＬＢＭモデルの速度成分を示す図である。 ＬＢＭモデルの速度成分を示す図である。 物理的プロセスシミュレーションシステムが辿る手順のフローチャートである。 マイクロブロックの斜視図である。 図３のシステムが用いる格子構造を示す図である。 図３のシステムが用いる格子構造を示す図である。 可変分解能技術を示す図である。 可変分解能技術を示す図である。 表面のファセットの影響を受ける領域を示す図である。 ボクセルから表面への粒子の移動を示す図である。 表面から表面への粒子の移動を示す図である。 表面動力学を行う手順のフローチャートである。 異なるサイズのボクセル間の界面を示す図である。 可変分解能条件下でファセットとの相互作用をシミュレートする手順のフローチャートである。 所与の位置における流れが層流であるか、それとも乱流であるかに基づいて、システム内で各ファセット又はサーフェルに壁剪断応力値を割り当てる手順のフローチャートである。

Ａ．壁剪断応力をモデル化する方法
複雑な流体流シミュレーションを完成させる際には、層流境界層の流れと乱流境界層の流れの間の壁剪断応力の違いを考慮することが有益である。

境界層が完全に乱流であり、境界層がゼロ圧力勾配下で発達中の場合には、普遍的な壁法則の適用が有効であり信頼性が高い。しかしながら、高レイノルズ数の壁面流では、必ずしもこの条件が満たされるわけではない。実際に、流体力学装置の前縁部近くでは、発達中の境界層の流れが完全に乱流でなくむしろ層流であることが多く、結果として得られる壁剪断応力の値に影響が及ぶ可能性がある。揚力及び効力のような全体的な流れ特性を予測することに対する影響は、特に流線形のボディでは有意なものとなり得る。従って、本明細書では、固体表面にわたって流れがどこで（及びいつ）層流又は乱流であるかを識別するための方法及びシステムについて説明する。流れがどこで層流であるかを識別することにより、壁剪断応力のためのモデルが、（例えば、特定の位置における表面摩擦／壁剪断応力の値を修正することなどによって）層流状況を適切かつ自動的に考慮することができる。具体的には、表面上の全ての位置について、その位置における流れが層流であるか、それとも乱流であるかを判定するための計算を行うことができ、その位置における流れが層流である場合には、流動力学シミュレーションにおいて壁剪断応力を表す第１の値を使用することができ、その位置における流れが乱流である場合には、流動力学シミュレーションにおいて壁剪断応力を表す第２の異なる値を使用することができる。本明細書で説明するシステム及び方法では、壁剪断応力の値が規定されている限り、格子ボルツマンの境界条件が、任意の幾何学形状の壁の運動量フラックス（すなわち、壁剪断応力）を保障する。本明細書で説明するシステム及び方法は、乱流運動エネルギーレベルと表面位置の壁剪断応力値との比較に基づいて、流れが乱流であるか否かをファセット／サーフェル毎に識別することを含み、乱流運動エネルギーレベルが壁剪断応力値以上である場合には流れが乱流であり、そうでなければ層流である。この流れが層流であるか、それとも乱流であるかについての判定に基づき、適切な壁剪断応力値を特定して表面位置に適用する（例えば、層流の領域には層流の壁剪断応力値を割り当て、乱流の領域には乱流の壁剪断応力値を割り当てる）。

この壁剪断応力をモデル化する手法は、マサチューセッツ州ＢｕｒｌｉｎｇｔｏｎのＥｘａ Ｃｏｒｐｏｒａｔｉｏｎ社から入手可能なＰｏｗｅｒＦＬＯＷシステムなどの、格子ボルツマン法に基づく時間陽的ＣＦＤ／ＣＡＡ解法と併せて使用することができる。ＬＢＭは、巨視的連続体方程式の離散化に基づく方法とは異なり、「メゾスコピック」ボルツマン運動方程式から開始して巨視的流体動力学を予測する。結果として得られる圧縮性非定常解法は、空力音響学及び純音響問題などの様々な複雑な流れの物理的特性を予測するために使用することができる。以下では、ＬＢＭベースのシミュレーションシステムの一般的な考察を行い、その後、このようなモデル化をサポートするための流体流シミュレーションと共に使用できるスカラー解法について考察する。

Ｂ．シミュレーション空間モデル
ＬＢＭベースの物理的プロセスシミュレーションシステムでは、離散速度ｃ_iの集合において評価される分布関数値

によって流体流を表すことができる。分布関数の動特性は、ｆ_i（０）が平衡分布関数として知られている以下のように定められる方程式（４）によって規定され、

式中、

である。

この式は、分布関数

の時間的発展を記述する周知のボルツマン方程式である。左辺は、いわゆる「流動プロセス」に起因する分布の変化を表す。流動プロセスは、流体ポケットがグリッド位置で始まり、次いで速度ベクトルの１つに沿って次のグリッド位置に移動する場合のものである。この時点で、「衝突要因」、すなわち流体の始動ポケットに対する流体の近傍ポケットの作用が算出される。流体は、別のグリッド位置にのみ移動することができ、従って全ての速度の全ての要素が共通速度の倍数になるように速度ベクトルを正しく選択する必要がある。

一次方程式の右辺は、上述の「衝突演算子」であり、流体のポケット間の衝突に起因する分布関数の変化を表す。ここで使用する衝突演算子の特定の形態は、Ｂｈａｔｎａｇａｒ， Ｇｒｏｓｓ及びＫｒｏｏｋ（ＢＧＫ）に起因する。この衝突演算子は、分布関数を二次方程式で与えられる規定値に近づけ、「平衡」形態となる。

このシミュレーションからは、質量ρ及び流体速度ｕなどの従来の流体変数が、式（３）における単純な加算として得られる。ここでは、ｃ_i及びｗ_iの集合値がＬＢＭモデルを定義する。ＬＢＭモデルは、拡張可能なコンピュータプラットフォーム上に効率的に実装され、時間的に不安定な流れ及び複雑な境界条件に関して優れたロバスト性で実行することができる。

ボルツマン方程式から流体系の運動の巨視的方程式を得る標準的技術は、完全なボルツマン方程式の逐次近似を考慮するチャップマン・エンスコッグ（Ｃｈａｐｍａｎ−Ｅｎｓｋｏｇ）法である。

流体系では、密度の微小変動が音速で移動する。ガス系では、一般に音速が温度によって決まる。流れにおける圧縮率の作用の重要性は、マッハ数として知られている、特性速度と音速との比率によって測定される。

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

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

１０１の速度を含む３Ｄ−２モデル、３７の速度を含む２Ｄ−２モデルなどの、より複雑なモデルを用いることもできる。

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

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

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

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

格子の分解能は、シミュレーション中のシステムのレイノルズ数に基づいて選択することができる。レイノルズ数は、流れの粘性（ν）、流れにおけるオブジェクトの特徴的長さ（Ｌ）、及び流れの特徴的速度（ｕ）に関連し、次式となる。

オブジェクトの特徴的長さは、オブジェクトのラージスケールの特徴要素を表す。例えば、マイクロデバイスの周囲の流れをシミュレートしている場合には、マイクロデバイスの高さを特徴的長さと見なすことができる。小領域のオブジェクト（例えば、自動車のサイドミラー）の周囲の流れが対象である場合、シミュレーションの分解能を高めることができ、或いは分解能が高まった範囲を対象領域の周囲で利用することができる。ボクセルの次元は、格子の分解能が高まるにつれて減少する。

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

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

各状態ｉは、特定のエネルギーレベル（すなわち、０、１又は２のエネルギーレベル）における異なる速度ベクトルを表す。各状態の速度ｃiは、以下のように３次元の各々の「速度」で表される。

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

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

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

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

図５Ａ及び図５Ｂを参照すると、表面Ｓ（図３Ａ）が、シミュレーション空間（図５Ｂ）内にファセットＦαの集合（サーフェルとも呼ばれる）として以下のように示されている。

式中、αは特定のファセットを示すインデックスである。ファセットは、ボクセル境界に限定されないが、通常は、ファセットが比較的少数のボクセルに影響を及ぼすように、ファセットに隣接するボクセルのサイズとほぼ同じか、又はそれよりもわずかに小さいサイズにされる。ファセットには、表面動力学を実施する目的で特性が割り当てられる。具体的には、各ファセットＦ_αは、単位法線（ｎ_α）、表面積（Ａ_α）、中心位置（ｘ_α）、及びファセットの表面動力学的特性を記述するファセット分布関数（ｆ_i（α））を有する。

図６を参照して分かるように、処理効率を改善させるために、シミュレーション空間の異なる領域では異なるレベルの分解能を用いることができる。通常は、オブジェクト６５５の周囲の領域６５０が最も重要であり、従って最も高い分解能でシミュレートされる。粘性の作用はオブジェクトからの距離と共に減少するので、オブジェクト６５５から間隔を空けた領域６６０、６６５のシミュレーションには低レベルの分解能（すなわち、拡大したボクセル体積）が用いられる。同様に、図７に示すように、オブジェクト７７５のそれほど重要ではない特徴の周囲の領域７７０のシミュレーションには低レベルの分解能を用いることができ、一方でオブジェクト７７５の最も重要な特徴の周囲の領域７８０のシミュレーションには最高レベルの分解能が用いられる。中心から離れた領域７８５は、最も低レベルの分解能及び最も大きなボクセルを用いてシミュレートされる。

Ｃ．ファセットの影響を受けるボクセルの識別
図３を参照して、シミュレーション空間をモデル化したら（ステップ３０２）、１又はそれ以上のファセットの影響を受けるボクセルを識別する（ステップ３０４）。ボクセルは、数多くの形でファセットの影響を受けることがある。まず、１又はそれ以上のファセットが交差するボクセルは、交差していないボクセルに比べて体積が小さいという点で影響を受ける。この影響が生じるのは、ファセット及びそのファセットによって表される表面下の物質がボクセルの一部を占有するからである。部分因子Ｐ_f（ｘ）は、ファセットの影響を受けないボクセルの部分（すなわち、流体又は流れをシミュレートするその他の物質が占有し得る部分）を示す。交差していないボクセルでは、Ｐ_f（ｘ）は１に等しい。

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

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

ファセットＦ_αは、状態の速度ベクトルがファセットに向けられた時（｜ｃ_iｎ_i｜＜０）に体積Ｖ_iαから粒子を受け取り、状態の速度ベクトルがファセットから離れる方向に向けられた時（｜ｃ_iｎ_i｜＞０）に粒子を領域に移動させる。以下で説明するように、この表現は、平行六面体Ｇ_iαの一部を別のファセットが占有し、内側角部などの非凸形の特徴の近傍で生じ得る条件となる場合には修正しなければならない。

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

平行六面体Ｇ_iαが重なるボクセルＮ（ｘ）の部分は、Ｖ_iα（ｘ）として定義される。この項を用いると、ボクセルＮ（ｘ）とファセットＦ_αの間を移動する状態ｉの粒子の流束Γ_iα（ｘ）は、以下のようにボクセル内の状態ｉの粒子の密度（Ｎ_i（ｘ））にボクセルと重なる領域の体積（Ｖ_iα（ｘ））を乗じたものに等しくなる。

平行六面体Ｇ_iαが１又はそれ以上のファセットと交差している場合には、以下の条件が真である。

式中、第１の総和は、Ｇ_iαが重なる全てのボクセルに相当し、第２項は、Ｇ_iαと交差する全てのファセットに相当する。平行六面体Ｇ_iαが他のファセットと交差しない場合には、この式が次のように短縮される。

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

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

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

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

このことから、ファセットＦ_αに向けられる各状態ｉ（ｃ_iｎ_α＜０）では、ボクセルがファセットＦ_αに提供する粒子の数が以下のようになる。

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

３．ファセットからファセットへの移動
次に、ファセット間で粒子を移動させる（ステップ３１０）。ファセットＦ_αの流入状態（ｃ_iｎ_α＜０）にある平行六面体Ｇ_iαが別のファセットＦ_βと交差する場合、ファセットＦ_αが受け取る状態ｉの粒子の一部は、ファセットＦ_βから得られるようになる。具体的には、ファセットＦ_αは、前回の時間増分中には、ファセットＦ_βによって生成された状態ｉの粒子の一部を受け取る。この関係を図１０に示しており、ここでは、ファセットＦ_βが交差する平行六面体Ｇ_iαの一部１０００が、ファセットＦ_αが交差する平行六面体Ｇ_iβの一部１００５に等しい。上述したように、交差部分はＶ_iα（β）として示している。この項を用いて、ファセットＦ_βとファセットＦ_αの間の状態ｉの粒子の流束を以下のように記述することができる。

式中、Γｉ（β，ｔ−１）は、前回の時間増分中にファセットＦ_βによって生成された状態ｉの粒子の測定値である。このことから、ファセットＦ_αに向けられる各状態ｉ（ｃ_iｎ_α＜０）では、他のファセットがファセットＦ_αに提供する粒子の数が以下のようになる。

また、ファセット内への状態ｉの粒子の全流束は以下のようになる。

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

ファセット分布関数は、ファセットから出力流束を生成するためのシミュレーションツールであり、必ずしも実際の粒子を表すとは限らない。正確な出力流束を生成するには、値を分布関数の他の状態に割り当てる。以下のように、内部の状態を取り込むための上述の技術を用いて外部の状態を取り込む。
ｃｉｎα≧０の場合、

式中、Γ_iOTHER（α）は、Γ_iIN（α）を生成するため上述の技術を用いて求められるが、この技術を流入状態（ｃ_iｎ_α＜０）以外の状態（ｃ_iｎ_α≧０）に適用する。代替の手法では、前回の時間ステップからΓ_iOUT（α）の値を用いて、Γ_iOTHER（α）を以下のように生成することができる。

平行状態（ｃ_iｎ_α＝０）では、Ｖ_iα及びＶ_iα（ｘ）が共に０である。Ｎ_i（α）についての式では、Ｖ_iα（ｘ）が（Γ_iOTHER（α）ついての式の）分子に表れ、Ｖ_iαが（Ｎ_i（α）についての式の）分母に表れる。従って、平行状態のＮ_i（α）は、Ｖ_iα及びＶ_iα（ｘ）がゼロに近づくにつれてＮ_i（α）の極限として求められる。

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

４．ファセット表面動力学の実行
次に、上記で説明した４つの境界条件を満たすように各ファセットについて表面動力学を実行する（ステップ３１２）。図１１に、ファセットの表面動力学を実行する手順を示す。最初に、ファセットＦαに垂直な全ての運動量を、そのファセットにおける粒子の全ての運動量Ｐ（α）を求めることによって以下のように求める（ステップ１１０５）。
全てのｉに関し、

これにより、法線運動量Ｐｎ（α）は以下のように求められる。

その後、プッシュ／プル技術を用いて、この法線運動量を排除して（ステップ１１１０）Ｎ_n-（α）を生成する。この技術に従って、粒子は、法線運動量のみに影響を与えるように状態間を移動する。このプッシュ／プル技術は、米国特許第５，５９４，６７１号に記載されており、この特許は引用により本明細書に組み入れられる。

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

これらの流入流束分布及びボルツマン分布に基づいて、ファセットＦ_αの流出流束分布を求める（ステップ１１２０）。最初に、流入流束分布Γ_i（α）とボルツマン分布の差分を以下のように求める。

この差分を用いて、ｎ_αｃ_i＞０場合、流出流束分布は以下のようになる。

式中、ｉ＊は、状態ｉとは逆の方向を有する状態である。例えば、状態ｉが（１，１，０，０）である場合、状態ｉ＊は、（−１，−１，０，０）である。表面摩擦及びその他の要因を考慮するために、ｎ_αｃ_i＞０である場合、流出流束分布を以下のようにさらに精密化することができる。

式中、Ｃ_fは、（壁剪断応力とも呼ばれる）表面摩擦の関数であり、ｔ_iαは、ｎ_αに垂直な第１の接線ベクトルであり、ｔ_2αは、ｎ_α及びＴ_1αの両方に垂直な第２の接線ベクトルであり、ΔＮ_j,1及びΔＮ_j,2は、状態ｉのエネルギー（ｊ）及び示された接線ベクトルに対応する分布関数である。この分布関数は、以下のように求められる。

式中、ｊは、エネルギーレベル１の状態では１に等しく、エネルギーレベル２の状態では２に等しい。

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

このことから、接線速度は以下のように求められる。

式中、ρはファセット分布の密度である。

上述したように、流入流束分布とボルツマン分布の差分は以下のように求められる。

この時、流出流束分布は以下のようになる。

この式は、従来の技術によって求められる流出流束分布の第１の２つの線に相当するが、異常な接線流束の補正を必要としない。

いずれの手法を用いても、結果として得られる流束分布は運動量流束の条件を全て満たし、すなわち以下のようになる。

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

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

式中、インデックスｊは、状態ｉのエネルギーを意味する。次に、このエネルギーの差分を用いて以下のように差分項を生成する。
ｃ_jiｎ_α＞０である場合、

この差分項を用いて流出流束を修正すると、流束は以下のようになる。
ｃ_jiｎ_α＞０である場合、

この演算は、接線運動量流束を変化させずに質量及びエネルギー流束を補正するものである。流れがファセットの近傍でほぼ均一であり、平衡状態に近い場合には、この調整はわずかなものである。調整後に結果として得られる法線の運動量流束は、近傍平均特性に近傍の非均一又は非平衡特性に起因する補正を加えたものに基づいて、平衡圧力である値までわずかに変化する。

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

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

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

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

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

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

式中、Ｐｆ（ｘ）は、部分ボクセルの体積減少を表す。これにより、各状態ｉについて、ファセットからボクセルＮ（ｘ）に向けられる粒子の総数は以下のようになる。

ファセットからボクセルに粒子を散乱させた後、これらの粒子と周囲のボクセルから運ばれて来た粒子とを結合し、その結果を整数化すると、特定のボクセルにおける特定の方向は、アンダーフロー（負となる）又はオーバーフロー（８ビット実装の場合は２５５を超える）となる可能性がある。これにより、質量、運動量及びエネルギーを許可された値の範囲に収まるように切り捨てた後に、これらの物理量が増加又は減少するようになる。このような発生を防ぐために、境界外の質量、運動量及びエネルギーを、有害な状態を切り捨てる前に集計する。状態が属しているエネルギーについては、（アンダーフローに起因して）増加又は（オーバーフローに起因して）減少した値に等しい質量を、オーバーフロー又はアンダーフローの影響を受けない同じエネルギーを有する無作為に（又は順次的に）選択した状態に加え戻す。この質量及びエネルギーの追加に起因する追加の運動量を集計し、切り捨てによる運動量に加える。同じエネルギー状態にのみ質量を追加することにより、質量カウンターがゼロに達した時に質量及びエネルギーの両方が補正される。最後に、プッシュ／プル技術を用いて、運動量アキュームレータがゼロに戻るまで運動量を補正する。

７．流体力学の実行
最後に、流体動力学を実行する（ステップ３１８）。このステップは、マイクロ動力学又はイントラボクセルの演算と呼ぶことができる。同様に、移流手法は、インターボクセル演算と呼ぶことができる。以下で説明するマイクロ動力学演算は、ファセットにおいて粒子を衝突させてボルツマン分布を生じるために使用することもできる。

格子ボルツマン方程式モデルでは、ＢＧＫ衝突モデルとして知られている特定の衝突演算子によって流体動力学が保証される。この衝突モデルは、実際の流体系における分布の動力学を模倣するものである。この衝突プロセスは、式１及び式２の右辺によって明確に記述することができる。移流ステップ後、式３を用いて、流体系の保存された物理量、特に密度、運動量及びエネルギーを分布関数から取得する。これらの物理量から、式（２）の

によって示される平衡分布関数を式（４）によって十分に規定する。表１に示す速度ベクトル集合ｃ_iと重量の選択は、いずれも方程式２と共に巨視的挙動が正確な流体力学方程式に従うことを確実にする。

Ｅ．可変分解能
図１２を参照して分かるように、（図６及び図７に示す上述したような）可変分解能は、（以下、粗ボクセル１２００及び微細ボクセル１２０５呼ぶ）異なるサイズのボクセルを使用する。（以下の説明では、２つの異なるサイズのボクセルを参照し、説明するる技術は、追加の分解能レベルを提供するように３又はそれ以上の異なるサイズのボクセルに適用することができると理解されたい。）粗ボクセルと微細ボクセルの界面は、可変分解能（ＶＲ）界面１２１０と呼ばれる。

表面又はその近傍で可変分解能を使用する場合、ファセットは、ＶＲ界面の両側でボクセルと相互作用することができる。これらのファセットは、ＶＲ界面ファセット１２１５（Ｆ_αIC）又はＶＲ微細ファセット１２２０（Ｆ_αIF）として分類される。ＶＲ界面１２１５は、ＶＲ界面の粗い側に位置し、微細ボクセル内に延びる粗い平行六面体１２２５を有するファセットである。（粗い平行六面体は、ｃ_iが粗ボクセルの寸法に従って寸法決めされたものであり、微細平行六面体は、ｃ_iが微細ボクセルの寸法に従って寸法決めされたものである。）ＶＲ微細ファセット１２２０は、ＶＲ界面の微細側に位置し、粗ボクセル内に延びる微細平行六面体１２３０を有するファセットである。界面ファセットに関連する処理は、粗いファセット１２３５（Ｆ_αC）及び微細ファセット１２４０（Ｆ_αF）との相互作用を含むこともできる。

両タイプのＶＲファセットに関し、表面動力学は、微細スケールで実行され、上述したような演算を行う。しかしながら、粒子がＶＲファセットとの間で移流する方法に関しては、ＶＲファセットは他のファセットとは異なる。

ＶＲファセットとの相互作用は、図１３に示す可変分解能手順１３００を用いて処理される。この手順のほとんどのステップは、非ＶＲファセットとの相互作用について上述した同等のステップを用いて行われる。手順１３００は、各々が微細な時間ステップに対応する２つの位相を含む粗い時間ステップ（すなわち、粗ボクセルに対応する時間間隔）中に実行される。ファセットの表面動力学は、各微細時間ステップ中に実行される。このため、ＶＲ界面ファセットＦ_αICは、同一サイズ及び同一配向の２つの微細ファセットと見なされ、それぞれが黒ファセットＦ_αICb及び赤ファセットＦ_αICrと呼ばれる。黒ファセットＦ_αICbは、粗い時間ステップ内の第１の微細時間ステップに関連し、赤ファセットＦ_αICrは、粗い時間ステップ内の第２の微細時間ステップに関連する。

最初に、第１の表面間移流段階によって粒子がファセット間を移動（移流）する（ステップ１３０２）。粒子は、ファセットＦ_αから延びてファセットＦ_βの後方に存在する粗い平行六面体（図１２の１２２５）未ブロック化部分の体積から、ファセットＦ_αから延びてファセットＦ_βの後方に存在する微細平行六面体（図１２の１２４５）の未ブロック化部分を差し引いたものに対応する重み付け係数Ｖ_-αβにより、黒ファセットＦ_αICbから粗いファセットＦ_βCに移動する。微細ボクセルのｃ_iの大きさは、粗ボクセルのｃ_iの大きさの２分の１である。上述したように、ファセットＦ_αの平行六面体の体積は以下のように定義される。

従って、ファセットの表面積Ａ_αは、粗い平行六面体と微細平行六面体の間で変化せず、単位法線ｎ_αは常に１の大きさを有するので、ファセットに対応する微細平行六面体の体積は、そのファセットの対応する粗い平行六面体の体積の２分の１である。

粒子は、ファセットＦ_αから延びてファセットＦ_βの後方に存在する微細平行六面体の未ブロック化部分の体積に対応する重み付け係数Ｖ_αβによって粗いファセットＦ_αCから黒ファセットＦ_βICbに移動する。

粒子は、重み付け係数Ｖ_αβによって赤ファセットＦ_αICrから粗いファセットＦ_βCに移動し、重み付け係数Ｖ_-αβによって粗いファセットＦ_αＣから赤ファセットＦ_βICrに移動する。

粒子は、Ｖ_αβの重み付け係数によって赤ファセットＦ_αICrから黒ファセットＦ_βICbに移動する。この段階では、黒から赤への移流は発生しない。また、黒ファセット及び赤ファセットは連続時間ステップを表すので、黒から黒への移流（又は赤から赤への移流）は決して発生しない。同様の理由から、この段階における粒子は、重み付け係数Ｖ_αβによって赤ファセットＦ_αICrから微細ファセットＦ_βIF又はＦ_βFに移動し、同じ重み付け係数によって微細ファセットＦ_αIF又はＦ_αFから黒ファセットＦ_αICbに移動する。

最後に、粒子は、同じ重み付け係数によって微細ファセットＦ_αIF又はＦ_αFから他の微細ファセットＦ_βIF又はＦ_βFに移動し、ファセットＦ_αから延びてファセットＦ_βの後方に存在する粗い平行六面体の未ブロック化部分の体積に対応する重み付け係数ＶＣ_αβによって粗いファセットＦ_αCから他の粗いファセットＦ_Cに移動する。

粒子が表面間を移流した後、第１の収集段階（ステップ１３０４〜１３１０）においてボクセルから粒子を収集する。粒子は、微細ファセットＦ_αFでは微細平行六面体を用いて微細ボクセルから収集され（ステップ１３０４）、粗いファセットＦ_αCでは粗い平行六面体を用いて粗ボクセルから収集される（ステップ１３０６）。次に、黒ファセットＦ_αIRb及び微細ＶＲファセットＦ_αIFでは、微細平行六面体を用いて粗ボクセル及び微細ボクセルの両方から粒子が収集される（ステップ１３０８）。最後に、赤ファセットＦ_αIRrでは、粗い平行六面体と微細平行六面体の差分を用いて粗ボクセルから粒子が収集される（ステップ１３１０）。

次に、微細ボクセル又はＶＲファセットと相互作用する粗ボクセルを微細ボクセルの集合に展開する（ステップ１３１２）。単一の粗い時間ステップ内で粒子を微細ボクセルに送る予定である状態の粗ボクセルを展開する。例えば、ファセットが交差しない適切な状態の粗ボクセルは、図４のマイクロブロックのように配向された８つの微細ボクセルに展開される。１又はそれ以上のファセットが交差する適切な状態の粗ボクセルは、完全微細ボクセル及び／又はいずれのファセットも交差していない粗ボクセルの一部に対応する部分微細ボクセルの集合に展開される。これらの展開から生じた粗ボクセル及び微細ボクセルの粒子密度Ｎ_i（ｘ）は等しいが、微細ボクセルは、粗ボクセルの部分因子及び他の微細ボクセルの部分因子とは異なる部分因子Ｐ_fを有することができる。

その後、微細ファセットＦ_αIF及びＦ_αFについて表面動力学を実行し（ステップ１３１４）た後に、黒ファセットＦ_αICbについて実行する（ステップ１３１６）。動力学は、図１１に示す上述した手順を用いて実行される。

次に、実際の微細ボクセルと粗ボクセルの展開によって生じた微細ボクセルとを含む微細ボクセル間を粒子が移動する（ステップ１３１８）。粒子の移動が完了すると、微細ファセットＦ_αIF及びＦ_αFから微細ボクセルに粒子を散乱させる（ステップ１３２０）。

粒子は、黒ファセットＦ_αICbから微細ボクセル（粗ボクセルを展開することにより生じる微細ボクセルを含む）へも散乱される（ステップ１３２２）。表面が存在しない時点でボクセルが粒子を受け取られる場合、粒子は微細ボクセルに散乱される。具体的には、ボクセルが（粗ボクセルの展開によって生じた微細ボクセルではなく）実際の微細ボクセルである場合、又はボクセルＮ（ｘ）を超える１つの速度単位であるボクセルＮ（ｘ＋ｃｉ）が実際の微細ボクセルである場合、或いはボクセルＮ（ｘ）を超える１つの速度単位であるボクセルＮ（ｘ＋ｃ_i）が粗ボクセルの展開によって生じた微細ボクセルである場合、粒子はボクセルＮ（ｘ）に散乱される。

最後に、微細ボクセル上で流体動力学を実行することによって第１の微細時間ステップが完了する（ステップ１３２４）。流体力学を実行するボクセルは、粗ボクセルを展開したことによって生じた微細ボクセルを含まない（ステップ１３１２）。

手順１３００は、第２の微細時間ステップ中に同様のステップを実行するものである。最初に、第２の表面間移流段階において粒子が表面間を移動する（ステップ１３２６）。粒子は、黒ファセットから赤ファセットへ、黒ファセットから微細ファセットへ、微細ファセットから赤ファセットへ、及び微細ファセットから微細ファセットへ移流される。

粒子が表面間を移流した後、第２の収集段階（ステップ１３２８〜１３３０）においてボクセルから粒子を収集する。粒子は、赤ファセットＦ_αIRrでは微細平行六面体を用いて微細ボクセルから収集される（ステップ１３２８）。微細ファセットＦ_αF及びＦ_αIFでも、微細平行六面体を用いて微細ボクセルから粒子が収集される（ステップ１３３０）。

その後、上述したように、微細ファセットＦ_αIF及びＦ_αFについて（ステップ１３３２）、粗いファセットＦ_αCについて（ステップ１３３４）、及び赤ファセットＦ_αICrについて表面動力学を実行する（ステップ１３３６）。

次に、微細ボクセルと、粗ボクセルを表す微細ボクセルとの間を粒子移動するように、微細分解能を用いてボクセル間で粒子を移動させる（ステップ１３３８）。次に、粗ボクセルとの間を粒子が移動するように、粗い分解能を用いてボクセル間で粒子を移動させる（ステップ１３４０）。

次に、結合ステップにおいて、粗ボクセルを表す微細ボクセル（すなわち、粗ボクセルの展開によって生じた微細ボクセル）が粗ボクセルに合体されている間に、粒子をファセットからボクセルに散乱させる（ステップ１３４２）。この結合ステップでは、粒子が、粗い平行六面体を用いて粗いファセットから粗ボクセルへ、微細平行六面体を用いて微細ファセットから微細ボクセルへ、微細平行六面体を用いて赤ファセットから微細ボクセル又は粗ボクセルへ、及び粗い平行六面体と微細平行六面体との間の差分を用いて黒ファセットから粗ボクセルへ散乱される。最後に、微細ボクセル及び粗ボクセルについて流体動力学を実行する（ステップ１３４４）。

Ｇ．壁剪断応力の決定
上述したように、流体流を解くために様々なタイプのＬＢＭを適用することができる。固体表面にわたる高レイノルズ数の境界層の流れの正確な予測は、これによって流体内で移動するボディに作用する空気／流体力学的力が求められるので、計算空気力学において非常に重要である。境界層の流れの影響を受ける１つの物理量に壁剪断応力（すなわち、表面摩擦）があり、この壁剪断応力は、以下に限定されるわけではないが、圧力の分布及び変動、流れの分離、並びに航空機の翼又は地上車などの流体力学装置の失速特性を含む基本的な流れ特性にさらに影響を与える。従って、流体流シミュレーションは、表面上の接線運動量流束の所望の変化を強制(enforce)するように、（壁剪断応力とも呼ばれる）表面摩擦に部分的に依拠する。表面上の特定の位置における壁剪断応力は、その表面における流体流に影響される。例えば、層流を示す表面領域の壁剪断応力は、乱流を示す表面領域の壁剪断応力よりも小さい。この壁剪断応力の違いにより、境界層における層流乱流遷移現象を予測し、この予測された層流乱流遷移に基づいて、表面上の領域（例えば、ファセット、サーフェル）に適した壁剪断応力値を選択することが有益となり得る。

本明細書で説明するシステム及び方法では、壁剪断応力の値をファセット毎（例えば、サーフェル毎）に割り当て、流体流シミュレーション中の各時間ステップにおいて再決定する。壁剪断応力を計算し、壁近傍のボクセルではなく各サーフェルの動力学において使用する。具体的には、ファセット毎及び時間ステップ毎に、表面位置における流体流が層流であるか、それとも乱流であるかに関する判定を行う。この判定に基づいて、層流の壁剪断応力値又は乱流の壁剪断応力値をファセットに割り当て、流体流シミュレーション中に使用する。サーフェルの剪断応力を計算すると、このサーフェル上の接線方向に反射された粒子の量を求めることができる。なお、壁近傍（例えば、サーフェルの組の近傍）のボクセルにおける粒子分布は、近隣のサーフェルから反射された粒子の影響を受ける。

境界層には、高濃度の非常に小さな流れ構造が存在するので、境界層の流れを数値的に解くための１つの方法は、壁の近傍に十分に高密度の格子点を適用してこれらの構造を求めることである。この手法は、非常に高いレイノルズ数の流れ及び複雑な幾何形状に関する現実の工学的問題のほとんどに関して計算的に困難であり実行不可能である。一方、壁付近の流れが、十分に発達した乱流境界層であると考えられる場合には、その全体的な速度プロファイルを、周知の「不偏的」な壁法則を満たすように十分に検討する。流体力学では、この壁法則が、ある地点における乱流の平均速度は、その地点から「壁」又は流体領域の境界までの距離の対数に比例するとされている。結果として、境界層内部の実際の速度分布を直接計算することなく、このような速度プロファイルから（実際には、壁からの何らかの所与の距離における単一の速度値のみから）局所的壁剪断応力の値を単純に推測することができる。これにより、壁付近の計算コストが大幅に削減される。この間接的手法は、壁付近の流れが伴い、乱流境界層として十分に発達している限り、ロバストかつ正確である。例えば、１９９７年３月２８日に出願された、物理的プロセスのコンピュータシミュレーションという名称の米国特許第５，９１０，９０２号に記載されているように、圧力勾配効果を組み入れることによって曲線的な幾何形状における流れの分離を処理するための拡張も行われており、この特許の内用は引用により本明細書に組み入れられる。高レイノルズ数の境界層流の物理的特性をモデル化することは、一般に乱流壁モデリングと呼ばれる。壁モデルを用いて乱流表面摩擦分布を予測することは、現実の工学的応用で見られる高レイノルズ数の流れの計算流体力学（ＣＦＤ）における常套手段である。

境界層におけるあらゆる場所の流れが完全に乱流であり、境界層がゼロ圧力勾配下で発達中である場合には、普遍的壁法則の適用が有効であり、信頼性が高い。実際には、高レイノルズ数の壁面流では、この条件が必ずしも満たされるとは限らない。実際に、流体力学装置の前縁部近くでは、発達中の境界層の流れが完全に乱流ではなくむしろ層流であることが多い。このことは、到来する流れの特性及び有利な圧力勾配効果などの多くの理由に起因することができる。層流境界層の流れと乱流境界層の流れの主な違いは、結果として生じるこれらの壁剪断応力の値である。同じ壁近傍の流速値では、前者の方が後者に比べて大幅に小さい。結果として、境界層が完全に乱流でない場合には、完全に発達した乱流境界層に基づく壁モデルから予測される壁剪断応力値が、物理的に正しい値よりも大幅に高くなる可能性がある。揚力及び効力のような全体的な流れ特性を予測することに対する影響は、特に流線形のボディに関して有意なものとなり得る。従って、本明細書で説明する方法及びシステムは、固体表面にわたってどこで（及びいつ）流れが層流又は乱流であるかを識別する。

壁剪断応力を計算するための１つの方法は、壁剪断応力と壁近傍の流体速度との関係から、例えば次式に従って壁剪断応力を推測するステップを含むことができる。

式中、Ｕ（ｙ）は、壁から距離ｙの地点で測定した局所的流体速度値であり、

は、その２乗値

が壁剪断応力となる摩擦速度であり、ｆ（ｙ＋）は、ｙ＋に関して定められる無次元流体速度プロファイルであり、

は、壁からの無次元距離であり、ｎｕ（ｖ）は、流体の動粘性率である。

壁剪断応力について解くように書き換えると、以下のようになる。

ｙ＋＞〜１５の場合、無次元流体速度プロファイルｆ（ｙ＋）は、乱流境界層の流れに関する方が層流境界層の流れに関するよりも小さいことが周知である。この結果、同じＵ（ｙ）に関しては、乱流の場合の方が層流の場合よりも剪断応力値が大きい。

ｆ（ｙ＋）の推定値は、以下に従って求めることができる。

層流境界層の場合、
全てのｙ＋値に関して、ｆ（ｙ＋）＝ｙ＋

乱流境界層の場合、
ｙ＋＜１０の場合には、ｆ（ｙ＋）＝ｙ＋
ｙ＋＞１０の場合には、ｆ（ｙ＋）＝（１／ｋ）＊Ｉｎ（ｙ＋）＋β
式（３９）
式中、ｋは約０．４であり、βは約５である。従って、ｙ＋＞１０の場合には、ｆ（ｙ＋）の推定値が周知の「壁法則」に従う。

このように、同じＵ（ｙ）値の場合の乱流剪断応力と層流剪断応力の比率を、所与のｙ＋値について求めることができる。典型的な高レイノルズ数の流れのシミュレーションのいくつかの例では、通常は、壁近傍のボクセルサイズが、約１００のｙ＋値を与えるように選択される。一例としてｙ＋＝１００を使用すると、結果として生じる乱流対層流の比率ｆは約１／６になり、従って剪断応力の比率は約３６になる。

乱流対層流を求める作業は、乱流物理学の全体的困難性に起因して困難になり得る。流れが層流であるか、それとも乱流であるかを判定するための１つの典型的な試みでは、例えば、Ｌａｎｇｔｒｙ他著、ＡＳＭＥターボＥｘｐｏ２００４会報、題名「局所的可変を用いた相関に基づく遷移モデル（Ａ Ｃｏｒｒｅｌａｔｉｏｎ−Ｂａｓｅｄ Ｔｒａｎｓｉｔｉｏｎ Ｍｏｄｅｌ Ｕｓｉｎｇ Ｌｏｃａｌ Ｖａｒｉａｂｌｅｓ）」に記載されるように、一方が間欠性のためのものであり、もう一方が運動量厚さのレイノルズ数に関する遷移開始基準のためのものである２つの輸送方程式に基づくモデルを考案することができる。これらの輸送方程式は、汎用ＣＦＤ法に従う自由な流れの乱流強度及び圧力勾配に基づいて、遷移に関する相関関係を実装するためのフレームワークを形成する。いくつかのベンチマークの流れの事例では、境界層における層流から乱流への遷移を捕捉できることが示された。残念ながら、このモデルでは、壁までの流れ場（すなわち、ｙ＋〜１）に２つのさらなる偏微分方程式（ＰＤＥ）を組み込む必要があることに加え、計算のグリッドの全ての地点において最も近い壁までの距離を計算する必要がある。このことは、計算的にことのほか費用が掛かるだけでなく、複雑な幾何形状においてアルゴリズム的に極めて不便でもある。

本明細書で説明するシステム及び方法は、（Ｌａｎｇａｔｒｙの著書に記載されている拡張のように）壁近傍の流体の乱流モデルへのさらなる拡張を行うのではなく、層流シナリオを適切かつ自動的に考慮する、壁剪断応力のための一般化された壁モデルを提供する。上述したように、乱流ＶＬＥＳモデルは、乱流の動的挙動を予測する上で適切なことが多いと考えられている。本明細書で説明するように、一般に、ＶＬＥＳモデルは、流体内でのモデル化及び壁におけるモデル化の両方を含む。壁モデルは、上述した壁法則に基づくものであり、圧力勾配が存在しない場合の完全に発達した境界層の流れにわたって有効である。しかしながら、実際の流れには圧力勾配が存在し、従って壁モデルは、局所的圧力勾配を考慮するように壁法則を修正した高度境界層法スキームを含むことができる。この流体モデルは、渦粘性を通じた乱流の物理的影響を考慮し、一方が乱流運動エネルギーｋの進化を規定し、もう一方が乱流散逸の進化イプシロンを規定する２つのＰＤＥによって部分的に求められる。この種のｋ−イプシロンモデルは、特にレイノルズの平均ナビエ・ストークス（ＲＡＮＳ）解法では一般的なものである。本発明者らの渦粘性の公式化は、局所的乱流、平均流れ特性及びレイノルズ数に基づく独自の成分も含み、不安的な乱流の平均特性及び時変特性の両方の正確な予測を可能にする。本発明者らのＶＬＥＳ乱流モデルでは、本明細書で説明する新たな遷移モデルスキームにとって重要な乱流運動エネルギーｋの時間的及び空間的分布を正確に予測することもできる。

このモデルが非乱流領域の動的挙動を正確に予測するには、壁モデルにおける壁剪断応力の定義が、層流領域と乱流領域の壁剪断応力の違いを考慮する必要がある。後者は、完全に発達した乱流境界層の壁法則仮説から構成される。本明細書で説明するモデル及びシステムは、ある領域が層流であるか、それとも乱流であるかを判定するために、さらなる流れ特性のために専用のＰＤＥを解いたり、或いは極めて薄い壁近傍領域において格子ボルツマン解法を用いたりする必要がない。また、上述したように、ファセット又はサーフェルベースの格子ボルツマン境界条件が、（例えば、式（３８〜３９）に関連して上述したように）任意の幾何形状の壁における（壁剪断応力を含む）運動量フラックステンソルを、壁剪断応力値が規定されている限り正確かつロバストに定義する。従って、壁剪断応力を適切に割り当てることにより、壁の境界層におけるメッシュ解像度に対してさらに条件を強要する必要がなくなる。従って、壁の境界層におけるメッシュ解像度を、シミュレートするシステムの他の領域のメッシュ解像度と同一又は同様（例えば、１０％以内）にすることができる。

図１４に、所与の位置における流れが層流であるか、それとも乱流であるかに基づいてシステム内の各ファセット又はサーフェルに壁剪断応力値を割り当てる方法を示す。この方法は、壁剪断応力が求められて、流体流シミュレーションの各時間ステップにおいて各ファセット／サーフェルに適用されるように繰り返される。従って、特定のファセット／サーフェルでは、壁剪断応力が、特定の表面位置における流体流がその時点に層流であるか、それとも乱流であるかに基づいて時間関数として異なることができる。同様に、流れが層流から乱流に、又は乱流から層流に遷移する領域では、隣接するファセット／サーフェルが異なる壁剪断応力値を有することができる。

この方法は、固体表面のありとあらゆる局所的領域（例えば、システム内の各ファセット／サーフェル）に関して壁剪断応力値又はその他の壁運動量フラックステンソル特性値などの２組の特性を計算するステップを含む。具体的には、最も単純な例では、システム（１４１０）内の乱流流体流の仮定に基づいて第１の壁剪断応力値を計算する。この計算は、完全に発達した乱流境界層プロファイルのみに基づく。また、壁近傍の層流の仮定に基づいて第２の独立した壁剪断応力値を計算する。後者（すなわち、層流の仮定に基づいて計算した壁剪断応力）は、所与の壁付近の流速値について前者（すなわち、乱流の仮定に基づいて計算した壁剪断応力）よりも実質的に小さな壁剪断応力値を有する。

流動力学シミュレーションで使用すべき実際の壁剪断応力は、壁付近の局所的流れが層流であるか、それとも乱流であるかに応じて、２つの異なる計算した壁剪断応力値のいずれかから選択される。流れが層流であるか、それとも乱流であるかは、第１及び第２の境界層計算の少なくとも一方の結果を基準と比較することによって判定される。１つの例では、この基準に入るキーパラメータが、（乱流運動エネルギーの値に関連する）局所的乱流強度レベルである。この乱流運動エネルギー値は、乱流についてのものの方が層流についてのものよりも高い。従って、基準は、局所的乱流運動エネルギーレベルが、運動量フラックステンソル、例えば乱流の壁剪断応力に基づく特性によって示される壁における加力に関連する別の基本的流体量の値を超えるかどうかに基づく。従って、この方法は、乱流運動エネルギーを、壁付近の乱流状態の仮定に基づいて計算された乱流壁運動量フラックステンソルの測定値と比較するステップを含む。さらに、壁付近の乱流は完全に発達していると仮定し、壁付近の乱流運動エネルギーの測定値を、この計算した乱流壁運動量フラックステンソルの測定値と比較するための基準を適用する（１４１６）。基準が満たされた場合には、所与の表面位置における流れを乱流であると仮定し、方法は、計算した乱流壁剪断応力値をファセット／サーフェルに割り当てる／適用する（１４１８）。一方、基準が満たされない場合には、この表面位置における流れを層流であると仮定し、方法は、計算した層流壁剪断応力値をファセット／サーフェルに割り当てる／適用する（１４２０）。その後、方法は、適用された壁剪断応力を用いて流動力学シミュレーションを行う（１４２２）。この乱流運動エネルギーレベルの測定値は、完全に発達した乱流境界層では、境界層の内部における乱流運動エネルギーはほぼ一定であり、（オーダー１の）一定の定数によって乱流壁剪断応力に比例することが知られているという論拠に基づいて合理的であり望ましい。従って、このロバストな測定値を、流れが乱流であるか否かを識別するために適用することができ、すなわち乱流運動エネルギーレベルが、壁剪断応力値に関連する運動エネルギーレベル以上である場合には乱流であり、そうでなければ層流である。

従って、一般に、本明細書で説明したシステム及び方法は、層流壁剪断応力値を求めるだけでなく、固体表面上のどこにこのような層流壁剪断応力値が存在し付加されているかも自動的に識別する。

上述の例では、運動エネルギーを閾値と比較して、流体流が乱流であるか、それとも層流であるかを判定したが、いくつかの例では、局所的乱流強度を使用することもできる。局所的乱流強度は、乱流運動エネルギー値ｋを局所的流体速度値と組み合わせたものから推測することができる。本明細書で説明した方法及びシステムでは、運動エネルギー又は乱流強度のいずれを使用してもよい。

上述の例では、局所的流れが層流であるか、それとも乱流であるかに基づいて、適用する壁剪断応力値を、乱流壁剪断応力値又は層流壁剪断応力値のいずれかであるように選択したが、いくつかの例では、壁付近の局所的流れが層流であるか、それとも完全に乱流であるか、それともこれらの何らかの中間であるかに応じて、流動力学シミュレーションで使用する壁剪断応力が、２つの異なる計算した壁剪断応力値のいずれか又はこれらの組み合わせから選択される。

このような組み合わせの一例として、局所的乱流運動エネルギーが閾値に近い時には、一方又は他方を選択する代わりに、層流剪断応力と乱流剪断応力を組み合わせた値が使用される。これにより、閾値の近傍における壁モデルスキームの変化を円滑にする「遷移」境界層が形成される。この組み合わせた剪断応力値は、例えば、２つの壁剪断応力値に重み付け係数を適用することにより、純粋な層流値及び乱流値からあらゆる所望の方法で求めることができる。１つの特定の例では、以下の線形的に重み付けした平均化手順を使用することができる。

式中、

は、遷移的壁剪断応力値であり、

は、純粋な層流値であり、

は、純粋な乱流値であり、ｋは、局所的乱流運動エネルギーであり、ｋ₁は、遷移範囲の開始を示し、ｋ₂は、遷移範囲の終了を示す。遷移的（「中間の」）値は、局所的乱流運動エネルギー値ｋが閾値の何らかの割合の範囲内、すなわち±１０パーセント内にある時に使用される。上記で定めた変数を使用すると、ｋ₁＝０．９ｋ₀及びｋ₂＝１．１ｋ₀であることが暗示され、式中のｋ₀は（上述したように、本スキームの局所的乱流剪断応力値の一定倍である）「公称」閾値である。

以上、多くの実装について説明した。それにもかかわらず、特許請求の範囲及びその思想から逸脱することなく、様々な一般化及び拡張を行うことができると理解されるであろう。例えば、本明細書で説明した方法は、壁運動量フラックステンソル以外の係数を含むように拡張することができる。別の例として、本明細書で説明した方法は、層流に基づくものもあれば乱流に基づくものもある複数の境界層計算を行って、これらの計算の少なくとも１つの結果を閾値と比較すると同時に、この比較の結果を用いて、流体流シミュレーションで使用すべき層流境界層計算又は乱流境界層計算の結果を選択するあらゆる状況をカバーするように拡張することもできる。上記の例の少なくともいくつかにおいて説明した壁剪断応力の代わりに、又はこれに加えて使用できる流れ特性の例としては、以下に限定されるわけではないが、空間的及び時間的導関数を含む、壁運動量フラックステンソルの様々なテンソル不変特性、並びに関連する平均的な流れ及び乱流特性の局所的及び非局所的に定められた整数値が挙げられる。また、本発明者らのＶＬＥＳモデルは、局所的乱流運動エネルギーなどの不安定な乱流特性を正確に求めるための効率的な手段のうちの１つの可能性にすぎない。以下の特許請求の範囲は、流れ特性の時間空間分布を適切かつ正確に記載している限り、いずれかの特定の乱流モデルを使用することに依存するものではない。従って、他の様々な一般的な拡張された実装も、以下の特許請求の範囲に含まれる。

１４１０：乱流境界に基づいて壁剪断応力を計算
１４１２：層流境界に基づいて壁剪断応力を計算
１４１４：局所的乱流運動エネルギーを計算
１４１６：基準を満たすか？
１４１８：表面位置に乱流壁剪断応力値を適用
１４２０：表面位置に層流壁剪断応力値を適用
１４２２：適用した壁剪断応力を用いて流動力学シミュレーションを実行

Claims (26)

1. 層流乱流境界層遷移を含む流体流をコンピュータ上でシミュレートする方法であって、 層流境界層の流れのための固体表面上の接線運動量流束の変化についての第１の計算を行うステップであって、該第１の計算は第１の結果を与える、ステップと、
   乱流境界層の流れのための前記固体表面上の接線運動量流束の変化についての第２の計算を行うステップであって、該第２の計算は第２の結果を与える、ステップと、
   前記層流境界層の流れのための前記第１の計算の前記第１の結果及び前記乱流境界層の流れのための前記第２の計算の前記第２の結果の少なくとも１つを基準と比較するステップと、
   前記基準との前記比較に基づいて、前記固体表面を表す複数の要素の少なくともいくつかについて、前記層流境界層の流れのための前記第１の計算の結果又は前記乱流境界層の流れのための前記第２の計算の結果を選択するステップと、
   ある体積の流体の活動のシミュレーションを行うステップと、
   を含み、前記シミュレーションは、前記複数の要素のための前記選択された結果に部分的に基づく、ことを特徴とする方法。
2. 前記層流境界層の流れのための前記第１の計算を行うステップは、前記層流の壁運動量フラックステンソル特性を計算するステップを含み、
   前記乱流境界層の流れのための前記第２の計算を行うステップは、前記乱流の壁運動量フラックステンソル特性を計算するステップを含み、
   前記複数の要素の少なくともいくつかについて、前記第１の計算の結果又は前記第２の計算の結果を選択するステップは、前記層流壁運動量フラックステンソル特性又は前記乱流壁運動量フラックステンソル特性を選択するステップを含む、
   ことを特徴とする請求項１に記載の方法。
3. 前記層流乱流境界層遷移を求めるステップは、
   前記固体表面上の複数の要素の各々について、前記第１の計算に基づく第１の値と、前記第２の計算に基づく第２の値とを求めるステップと、
   前記第１及び第２の値の少なくとも一方を基準と比較することにより、前記複数の要素の少なくともいくつかについて、前記流れを層流又は乱流として分類するステップと、
   を含むことを特徴とする請求項１に記載の方法。
4. 前記第１の計算の結果又は前記第２の計算の結果を選択するステップは、更に、
   層流として分類された要素について、前記層流の壁運動量フラックステンソル特性を選択するステップと、
   乱流として分類された要素について、前記乱流の壁運動量フラックステンソル特性を選択するステップと、
   を含むことを特徴とする請求項３に記載の方法。
5. 前記第１の計算の前記結果は、層流壁運動量フラックステンソルの値を含み、
   前記第２の計算の前記結果は、乱流壁運動量フラックステンソルの値を含み、
   前記比較は、乱流強度の値を含む、
   ことを特徴とする請求項３に記載の方法。
6. 前記第１の計算を行うステップは、前記表面上の複数の要素の各々について、層流壁運動量フラックステンソルの値を計算するステップを含み、前記第２の計算を行うステップは、前記固体表面上の複数の要素の各々について、乱流壁運動量フラックステンソルの値を計算するステップを含み、
   前記比較するステップは、前記表面上の前記複数の要素の各々について、計算された乱流強度の値と、前記乱流壁運動量フラックステンソルの値とを比較するステップを含み、 前記第１の計算の結果又は前記第２の計算の結果を選択するステップは、前記表面上の前記複数の要素の少なくともいくつかについて、前記乱流強度の値と前記乱流壁運動量フラックステンソルの値との前記比較に基づいて、前記計算された乱流壁運動量フラックステンソル特性及び層流壁運動量フラックステンソル特性の一方を選択するステップを含む、
   ことを特徴とする請求項１に記載の方法。
7. 前記表面上の前記複数の要素の各々について、前記乱流強度の値と前記乱流壁運動量フラックステンソルの値とを比較するステップは、前記乱流強度の値が前記層流壁運動量フラックステンソルの値よりも大きいかどうかを決定するステップを含み、
   前記表面上の前記複数の要素の少なくともいくつかについて、前記計算された乱流壁運動量フラックステンソル及び層流壁運動量フラックステンソルの一方を選択するステップは、特定の要素について、前記乱流強度の値が前記乱流壁運動量フラックステンソルの値よりも大きい場合には、前記乱流壁運動量フラックステンソルを選択し、前記乱流強度の値が前記乱流壁運動量フラックステンソルの値よりも小さい場合には、前記層流壁運動量フラックステンソルの値を選択するステップを含む、
   ことを特徴とする請求項６に記載の方法。
8. 局所的乱流強度の値を計算するステップは、局所的乱流運動エネルギーの値を計算するステップを含む、
   ことを特徴とする請求項６に記載の方法。
9. 所与の壁近傍の流体速度について、前記乱流壁運動量フラックステンソルの値は、前記層流壁運動量フラックステンソルの値よりも大きい、
   ことを特徴とする請求項６に記載の方法。
10. 前記体積の流体の活動のシミュレーションを行うステップは、
    状態ベクトルに対して、異なる運動量状態の要素間の相互作用をモデルに従ってモデル化する相互作用演算を行うステップと、
    前記状態ベクトルの組の第１の移動演算を行って、前記体積内の新たなボクセルへの要素の移動を前記モデルに従って反映するステップと、
    を含むことを特徴とする請求項１に記載の方法。
11. 前記第２の計算は、速度プロファイル及び前記固体表面からの距離に基づいて乱流壁運動量フラックステンソルの値を求めるための計算を含む、
    ことを特徴とする請求項１に記載の方法。
12. 前記表面上の前記複数の要素の少なくともいくつかについて、前記層流境界層の流れのための前記第１の計算の前記結果、及び前記乱流境界層の流れのための前記第２の計算の前記結果の加重平均に基づく値を求めるステップをさらに含む、
    ことを特徴とする請求項１に記載の方法。
13. 前記表面上の前記複数の要素の少なくともいくつかについて、前記第１の値及び前記第２の値の組み合わせに基づく壁運動量フラックステンソル特性を求めるステップをさらに含む、
    ことを特徴とする請求項３に記載の方法。
14. 前記第２の計算は、局所的乱流運動エネルギー及び局所的流体速度に基づいて乱流壁運動量フラックステンソルの値を求めるための計算を含む、
    ことを特徴とする請求項１に記載の方法。
15. コンピュータ可読媒体内で有形的に具体化されるコンピュータプログラム製品であって、実行時に層流乱流境界層遷移を含む物理的プロセスの流体流のシミュレーションを行う命令を含み、前記コンピュータプログラム製品は、コンピュータに、
    層流境界層の流れのための固体表面上の接線運動量流束の変化についての第１の計算であって第１の結果を与える第１の計算を行わせ、
    乱流境界層の流れのための前記固体表面上の接線運動量流束の変化についての第２の計算であって第２の結果を与える第２の計算を行わせ、
    前記層流境界層の流れのための前記第１の計算の前記第１の結果及び前記乱流境界層の流れのための前記第２の計算の前記第２の結果の少なくとも１つを基準と比較させ、
    前記基準との前記比較に基づいて、前記固体表面を表す複数の要素の少なくともいくつかについて、前記層流境界層の流れのための前記第１の計算の結果又は前記乱流境界層の流れのための前記第２の計算の結果を選択させ、
    ある体積の流体の活動のシミュレーションを行わせる、
    ように構成され、前記シミュレーションは、前記複数の要素のための前記選択された結果に部分的に基づく、ことを特徴とするコンピュータプログラム製品。
16. 前記層流境界層の流れのための前記第１の計算を行うための命令は、前記層流の壁運動量フラックステンソル特性を計算するための命令を含み、
    前記乱流境界層の流れのための前記第２の計算を行うための命令は、前記乱流の壁運動量フラックステンソル特性を計算するための命令を含み、
    前記第１の計算の結果又は前記第２の計算の結果を選択するための命令は、前記層流壁運動量フラックステンソル特性又は前記乱流壁運動量フラックステンソル特性を選択するための命令を含む、
    ことを特徴とする請求項１５に記載のコンピュータプログラム製品。
17. 前記境界層の前記層流乱流遷移を求めるための命令を含み、
    前記表面上の複数の要素の各々について、前記第１の計算に基づく第１の値と、前記第２の計算に基づく第２の値とを求めるための命令と、
    前記第１及び第２の値の少なくとも一方を基準と比較することにより、前記複数の要素の少なくともいくつかについて、前記流れを層流又は乱流として分類するための命令と、 を更に含むことを特徴とする請求項１５に記載のコンピュータプログラム製品。
18. 前記第１の計算の結果又は前記第２の計算の結果を選択するための命令は、
    層流として分類された要素について、前記層流の壁運動量フラックステンソル特性を選択するための命令と、
    乱流として分類された要素について、前記乱流の壁運動量フラックステンソル特性を選択するための命令と、
    を含むことを特徴とする請求項１７に記載のコンピュータプログラム製品。
19. 前記第１の計算の前記結果は、層流壁運動量フラックステンソル特性の値を含み、
    前記第２の計算の前記結果は、乱流壁運動量フラックステンソル特性の値を含み、
    前記基準は、乱流強度の値を含む、
    ことを特徴とする請求項１７に記載のコンピュータプログラム製品。
20. 前記第１の計算を行うための命令は、前記表面上の複数の要素の各々について、層流壁運動量フラックステンソルの値を計算するための命令を含み、前記第２の計算を行うための命令は、前記固体表面上の複数のファセットの各々について、乱流壁運動量フラックステンソルの測定値を計算するための命令を含み、
    前記第１の計算の結果を前記第２の計算の結果と比較するための命令は、前記表面上の前記複数の要素の各々について、計算された乱流強度の値と、前記乱流壁運動量フラックステンソルの値とを比較するための命令を含み、
    前記第１の計算の結果又は前記第２の計算の結果を選択するための命令は、前記表面上の前記複数の要素の少なくともいくつかについて、前記乱流強度の値と前記乱流壁運動量フラックステンソルの値との前記比較に基づいて、前記計算された乱流壁運動量フラックステンソル特性及び層流壁運動量フラックステンソル特性の一方を選択するための命令を含む、
    ことを特徴とする請求項１５に記載のコンピュータプログラム製品。
21. 物理的プロセスの流体流のシミュレーションを行うためのコンピュータシステムであって、
    層流境界層の流れのための固体表面上の接線運動量流束の変化についての第１の計算であって第１の結果を与える第１の計算を行い、
    乱流境界層の流れのための前記固体表面上の接線運動量流束の変化についての第２の計算であって第２の結果を与える第２の計算を行い、
    前記層流境界層の流れのための前記第１の計算の前記第１の結果及び前記乱流境界層の流れのための前記第２の計算の前記第２の結果の少なくとも１つを基準と比較し、
    前記基準との前記比較に基づいて、前記固体表面を表す複数の要素の少なくともいくつかについて、前記層流境界層の流れのための前記第１の計算の結果又は前記乱流境界層の流れのための前記第２の計算の結果を選択し、
    ある体積の流体の活動のシミュレーションを行う、
    ように構成され、前記シミュレーションは、前記複数の要素のための前記選択された結果に部分的に基づく、ことを特徴とするシステム。
22. 前記システムは、更に、
    前記層流境界層のための前記第１の計算を行って、前記層流の壁運動量フラックステンソル特性を計算し、
    前記乱流境界層の流れのための前記第２の計算を行って、前記乱流の壁運動量フラックステンソル特性を計算するように構成され、
    前記システムは、更に、
    前記第１の計算の結果又は前記第２の計算の結果を選択するにあたって、前記層流壁運動量フラックステンソル特性又は前記乱流壁運動量フラックステンソル特性を選択するように構成される、
    ことを特徴とする請求項２１に記載のシステム。
23. 前記システムは、
    前記境界層の層流乱流遷移を求め、
    前記表面上の複数の要素の各々について、前記第１の計算に基づく第１の値と、前記第２の計算に基づく第２の値とを求め、及び、
    前記第１及び第２の値の少なくとも一方を基準と比較することにより、前記複数の要素の少なくともいくつかについて、前記流れを層流又は乱流として分類するように構成される、ことを特徴とする請求項２１に記載のシステム。
24. 前記システムは、
    層流として分類された要素について、前記層流の壁運動量フラックステンソル特性値を選択し、
    乱流として分類された要素について、前記乱流の壁運動量フラックステンソル特性を選択するする、
    ように更に構成される、特徴とする請求項２３に記載のシステム。
25. 前記第１の計算の前記結果は、層流壁運動量フラックステンソルの値を含み、
    前記第２の計算の前記結果は、壁運動量フラックステンソルの値を含み、
    前記比較は、乱流強度の値を含む、
    ことを特徴とする請求項２３に記載のシステム。
26. 前記システムは、更に、
    前記第１の計算を行って、前記表面上の複数の要素の各々について、層流壁運動量フラックステンソルの値を計算し、前記第２の計算を行って、前記固体表面上の複数の要素の各々について、乱流壁運動量フラックステンソルの値を計算するように構成され、
    前記システムは、更に、
    前記表面上の前記複数の要素の各々について、計算された乱流強度の値と、前記乱流壁運動量フラックステンソルの値とを比較するように更に構成されることにより、前記第１の計算の結果を前記第２の計算の結果と比較するように更に構成され、
    前記表面上の前記複数の要素の少なくともいくつかについて、前記乱流強度の値と前記乱流壁運動量フラックステンソルの値との前記比較に基づいて、前記計算された乱流壁運動量フラックステンソル特性及び層流壁運動量フラックステンソル特性の一方を選択するように更に構成されることにより、前記第１の計算の結果又は前記第２の計算の結果を選択するように更に構成される、
    ことを特徴とする請求項２１に記載のシステム。

Images