JP6825038B2 - 相対浸透率シミュレーションのための質量交換モデル - Google Patents

Description

〔優先権の主張〕
本出願は、２０１３年５月１６日出願の「相対浸透率シミュレーションにおける履歴効果を捕捉するための質量シンク／ソースモデル（Mass Sink/source Model for Capturing History Effects in Relative Permeability Simulations）」という名称の米国特許仮出願出願番号第６１／８２４，１００号に対する「３５ ＵＳＣ §１１９（ｅ）」の下での優先権を主張するものであり、その内容全体は、これにより引用によって組み込まれる。

本記述は、多孔質媒体（porous media）を通る多種流れ（multi-species flow）に対する相対浸透率の決定／推定を含む多孔質媒体を通る多種流れのコンピュータシミュレーションのような物理過程のコンピュータシミュレーションに関する。

地下で見出される炭化水素は、典型的に岩石層に存在する。これらの岩石層は、通常はある程度多孔質であり、多孔質媒体として分類することができる。この多孔質媒体からの炭化水素抽出は、炭化水素と混和しない流体（例えば、水）を使用して行われることが多い。多孔質媒体からの炭化水素抽出を理解するために、多孔質媒体及び多孔質媒体を通る流れを特徴付けるのにシミュレーションが使用される。

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

一般的に、本明細書は、シミュレーションが異なる飽和度（saturations）で実行される多孔質媒体を通る多種流れをシミュレートするための技術を説明する。本明細書に説明するシミュレーションは、選択された位置の部分集合（例えば、ボクセル）において流体質量（fluid mass）が異なる種の間で交換される質量シンク／ソース方法を使用する。このようにして飽和度を変えることができ、一方、以前の飽和条件からの流体種の空間分布の多くは不変のままであり、従って、１つの飽和条件でのシミュレーション結果からの情報は、別の飽和条件でのシミュレーションに影響を与える。

１つの一般的態様において、多孔質媒体を通る多種流れをシミュレートするコンピュータ実施型方法は、容積中に第１の種及び第２の種を含む多相流体（multi-phase fluid）の活動をシミュレートし、容積中の流体の活動が、容積内の要素の移動をモデル化するようにシミュレートされる段階と、第１の飽和値に対して容積中の第１の種の相対浸透率を決定する段階と、識別された交換領域（exchange regions）から第１の種の第１の質量を除去し、第１の種の第１の質量を第２の種の第２の質量で置換して容積に対する飽和値を修正する段階とを含む。本方法はまた、修正された飽和値に基づいて多相流体の活動をシミュレートする段階と、修正された飽和値に基づいて容積中の第１の種の相対浸透率を決定する段階とを含む。

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

本方法はまた、容積（volume）内の第１の種の移動を示す値に基づいて容積内の交換領域を識別する段階を含むことができる。

交換領域を識別する段階は、１組の対流ボクセル（convective voxels）を識別する段階を含むことができる。

容積中の対流ボクセルを識別する段階は、第１の種の高い流量（flow rate）を有するボクセルを識別する段階を含むことができる。

容積中の対流ボクセルを識別する段階は、速度が閾値を超えるボクセルを識別する段階を含む。

容積中の対流ボクセルを識別する段階は、ボクセルをそれらに関連付けられた速度に基づいてランク付けする段階と、最大速度を有するボクセルの一部分を選択する段階とを含むことができる。

第１の種の第１の質量を除去する段階及び第１の種の第１の質量を第２の種の第２の質量で置換する段階は、第１の種の第１の質量を第２の種の等しい質量で置換する段階を含むことができる。

識別された交換領域から第１の種の第１の質量を除去する段階及び第１の種の第１の質量を第２の種の等しい質量で置換する段階は、望ましい飽和値が得られるまで複数の時間ステップにわたって質量交換過程を実行する段階を含むことができる。

多相流体の活動をシミュレートする段階は、有効浸透率が収束するまで多相流体の活動をシミュレートする段階を含むことができる。

第１の種は、油とすることができ、第２の種は、水とすることができる。

容積は、多孔質媒体容積とすることができる。

多相流体の活動をシミュレートする段階は、重力駆動シミュレーションを使用して多相流体の活動を周期的ドメイン内でシミュレートする段階を含むことができる。

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

上述の技術の実装は、方法又は過程、システム又は装置、又はコンピュータアクセス可能媒体上のコンピュータソフトウエアを含むことができる。

一部の追加の態様において、炭化水素リザーバの生産性を推定する方法は、炭化水素リザーバ中の少なくとも油及び水を含む多相流体の活動をシミュレートし、炭化水素リザーバ中の流体の活動が、炭化水素リザーバ内の要素の移動をモデル化するようにシミュレートされる段階と、第１の飽和値に対して炭化水素リザーバ中の油の相対浸透率及び水の相対浸透率を決定する段階と、識別された交換領域から油の第１の質量を除去し、油の第１の質量を水の第２の質量で置換して容積に対する飽和値を修正する段階と、修正された飽和値に基づいて多相流体の活動をシミュレートする段階と、修正された飽和値に基づいて炭化水素リザーバ中の油の相対浸透率及び水の相対浸透率を決定する段階と、決定された相対浸透率に基づいて炭化水素リザーバから生成することができる油の量を推定する段階とを含む。

上述の技術の実装は、方法又は過程、システム又は装置、又はコンピュータアクセス可能媒体上のコンピュータソフトウエアを含むことができる。

システム及び方法及び技術は、多相流れ（multi-phase flow）に対するＳｈａｎ−Ｃｈｅｎ方法と格子ボルツマン方式（Lattice Boltzmann formulation）とのような様々なタイプの数値シミュレーション手法を使用して実施することができる。格子ボルツマン方式に関する追加の情報を本明細書で以下に説明する。しかし、本明細書に説明するシステム及び技術は、格子ボルツマン方式を使用するシミュレーションに限定されず、他の数値シミュレーション手法に適用することができる。

システム及び技術は、格子ボルツマン方式を使用する格子ガスシミュレーションを使用して実施することができる。従来の格子ガスシミュレーションは、各格子部位での限られた数の粒子を仮定し、粒子は、ビットのショートベクトルによって表される。各ビットは、特定の方向に移動する粒子を表している。例えば、ベクトル中の１ビットは、特定の方向に沿って移動する粒子の存在（１に設定時）又は不在（０に設定時）を表すことができる。そのようなベクトルは、６ビットを有することができ、例えば、値１１００００は、２つの粒子がＸ軸に沿って反対方向に移動し、かつＹ軸及びＺ軸に沿って移動する粒子がないことを示している。１組の衝突規則が、各部位での粒子間の衝突の挙動を支配する（例えば、１１００００ベクトルは、００１１００ベクトルになることができ、Ｘ軸に沿って移動する２つの粒子間の衝突は、Ｙ軸に沿って離れる２つの粒子を生成したことを示す）。規則は、ビットに対して置換を実行する（例えば、１１００００を００１１００に変形する）ルックアップテーブルに状態ベクトルを供給することによって実施される。粒子は、その後に、隣接部位に移動される（例えば、Ｙ軸に沿って移動する２つの粒子は、Ｙ軸に沿って左右に隣接部位に移動されると考えられる）。

強化されたシステムにおいて、各格子部位での状態ベクトル（state vector）は、粒子エネルギ及び移動方向の変動を与えるためにより多くのビットを含み（例えば、亜音速流れに対して５４ビット）、完全な状態ベクトルの部分集合を伴う衝突規則が使用される。更に強化されたシステムにおいて、１つよりも多い粒子が、各格子部位又はボクセル（これらの２つの用語は、本明細書を通して交換可能に使用される）で各運動量状態で存在することが許される。例えば、８ビット実装において、０〜２５５個の粒子は、特定のボクセルで特定の方向に移動している可能性がある。状態ベクトルは、１組のビットである代わりに、１組の整数であり（例えば、１組の８ビットバイトが、０〜２５５の範囲の整数を与える）、その各々は、与えられた状態にある粒子の数を表する。

更なる強化において、格子ボルツマン方法（ＬＢＭ）は、従来の計算流体力学（「ＣＦＤ」）手法で可能であるよりも深いレベルで複合的な幾何学形状における３Ｄ不安定圧縮性乱流過程をシミュレートするために流体のメゾスコピック表現を使用する。ＬＢＭ方法の概要を以下に与える。

ボルツマン−レベルメゾスコピック表現
流体システムは、いわゆる「メゾスコピック」レベルに対する動態方程式（kinetic equations）によって表すことができることが統計物理学で公知である。このレベルに対しては、個々の粒子の詳細な運動は決定する必要はない。これに代えて、流体の特性は、単一粒子位相空間、

を使用して定義される粒子分布関数によって表され、ｘは、空間座標であり、νは、粒子速度座標である。質量、密度、流量速度、及び温度のような典型的な流体力学的な量は、粒子分布関数の単純なモーメントである。粒子分布関数の力学は、ボルツマン方程式に従う。

ここで、Ｆ（ｘ，ｔ）は、（ｘ，ｔ）で外部又は自己矛盾なく発生された物体力を表している。衝突項Ｃは、様々な速度及び位置の粒子の相互作用を表している。衝突項Ｃに対して特定の形態を指定することなく、上述のボルツマン方程式は、希薄ガスの公知の状況（ボルツマンによって本来構成されたような）だけでなく全ての流体システムに適用可能であることを強調することが重要である。

一般的に、Ｃは、二点相関関数の複雑な多次元積分を含む。分布関数ｆのみで閉じたシステムを形成するために、並びに効率的な計算のために、最も便利かつ物理的に一貫した形態の１つは、公知のＢＧＫ演算子である。ＢＧＫ演算子は、衝突の詳細がたとえ何であろうとも、分布関数は、衝突：

を通じて

によって与えられる明確な局所平衡に近づくという物理的議論に従って構成され、ここで、パラメータτは、衝突を通じた平衡までの固有緩和時間を表している。粒子（例えば、原子又は分子）を扱うと、緩和時間は、典型的には定数として取られる。「混成」（水力−動力学）表現において、この緩和時間は、歪み率及び乱流運動エネルギなどのような流体力学的変数の関数である。すなわち、乱流は、局所的に決定された固有特性を有する乱流粒子（「渦」）のガスとして表すことができる。

ボルツマン−ＢＧＫ方程式の数値解は、ナビエ−ストークス方程式の解に優るいくつかの計算上の利点を有する。第１に、複雑な非線形項又は高次空間導関数が方程式にないことを直ちに認識することができ、すなわち、移流不安定性に関する問題がほとんどない。このレベルの記述では、方程式は、圧力を扱う必要がないので局所的であり、これは、アルゴリズム並列化に対してかなりの利点を提供する。２次空間導関数を有する拡散演算子がないという事実と共に線形移流演算子の別の望ましい特徴は、流体偏微分方程式（「ＰＤＥ」）の数学的な条件ではなく、実際に粒子が固体面と真に相互作用する方法を模倣するように非スリップ面又はスリップ面のような物理的境界条件を実現する際の容易さである。直接の利益の１つは、固体面上でインタフェースの移動を処理する問題がないということであり、これは、格子ボルツマンベースのシミュレーションソフトウエアが失敗なく複合乱流空気力学をシミュレートすることを可能にするのを補助する。これに加えて、有限粗度面のような境界からのある一定の物理特性も、力に組み込むことができる。更に、ＢＧＫ衝突演算子は、純粋に局所的であり、一方、自己矛盾のない物体力の計算は、近−隣接情報のみを通して達成することができる。その結果、ボルツマン−ＢＧＫ方程式の計算は、並列処理に実質的に適応させることができる。

格子ボルツマン方式
連続体ボルツマン方程式を解くことは、それが位置及び速度位相空間における積分微分方程式の数値評価を伴うということで有意な課題を表している。位置だけでなく速度位相空間も離散化することができることが観測された時に大規模な簡素化が行われ、それは、ボルツマン方程式の解のための効率的な数値アルゴリズムをもたらした。流体力学的な量は、高々最近隣接情報に依存する単純和の項で書くことができる。従来的には、格子ボルツマン方程式の方式は、速度の離散集合：

に対する粒子の推移を規定する格子ガスモデルに基づいていたが、この方程式は、連続体ボルツマン方程式の離散化としての第１の原理から系統的に導出することができる。その結果、ＬＢＥは、格子ガス手法に関連付けられた公知の問題を被らない。従って、位相空間における連続体分布関数ｆ（ｘ，ν，ｔ）を取り扱う代わりに、離散速度指数をラベル付けする下付き文字を有する離散分布の有限集合Ｆ（ｘ，ｔ）を追跡することが必要なだけである。巨視的記述の代わりにこの動態方程式を取り扱う重要な利点は、システムの位相空間の増大が、問題の局所性によってオフセットされるということである。

対称性の考慮に起因して、速度値の集合は、それらが構成空間内で張られた時にある一定の格子構造を形成するように選択される。そのような離散システムの力学は、形態：

を有するＬＢＥに従い、ここで、衝突演算子は、上述のように通常はＢＧＫ形態を取る。平衡分布形態の適正な選択により、格子ボルツマン方程式は正しい流体力学及び熱−流体力学を生じさせることを理論的に示すことができる。すなわち、

から導出された流体力学モーメントは、巨視的限界においてナビエ−ストークス方程式に従う。これらのモーメントは、

して定義され、ここで、ρ、ｕ、及びＴは、それぞれ、流体密度、速度、及び温度であり、Ｄは、離散化された速度空間の次元（物理空間次元に全く等しくない）である。

他の特徴及び利点は、図面及び特許請求の範囲を含む以下の説明から明らかであろう。

行き止まり細孔が取り付けられた垂直細孔の周期的システムを示す図である。 質量交換を使用せず、従って、履歴効果が含まれない図１に示す例示的システム内の細孔における影響を示し、初期流体分布が、全ての単一ボクセルに対して流体ドメイン全体にわたって密度を均一に分配する過程に基づいている図である。 質量交換を使用せず、従って、履歴効果が含まれない図１に示す例示的システム内の細孔における影響を示し、初期流体分布が、全ての単一ボクセルに対して流体ドメイン全体にわたって密度を均一に分配する過程に基づいている図である。 質量交換を使用せず、従って、履歴効果が含まれない図１に示す例示的システム内の細孔における影響を示し、初期流体分布が、全ての単一ボクセルに対して流体ドメイン全体にわたって密度を均一に分配する過程に基づいている図である。 質量交換を使用した図１に示す例示的システム内の細孔における履歴効果を示す図である。 質量交換を使用した図１に示す例示的システム内の細孔における履歴効果を示す図である。 質量交換を使用した図１に示す例示的システム内の細孔における履歴効果を示す図である。 ２つのＬＢＭモデルの速度成分を示す図である。 ２つのＬＢＭモデルの速度成分を示す図である。 物理過程シミューレーションシステムが従う手順の流れ図である。 マイクロブロックの斜視図である。 図６のシステムによって使用される格子構造の図である。 図６のシステムによって使用される格子構造の図である。 可変分解能技術を示す図である。 可変分解能技術を示す図である。 面のファセットによって影響を受ける領域を示す図である。 ボクセルから面への粒子の移動を示す図である。 面から面への粒子の移動を示す図である。 面力学を実行するための手順の流れ図である。 異なるサイズのボクセル間のインタフェースを示す図である。 可変分解能条件下のファセットとの相互作用をシミュレートするための手順の流れ図である。 異なる飽和度レベルでの相対浸透率を決定するための手順の流れ図である。 質量交換のための手順の流れ図である。

Ａ．相対浸透率シミュレーションへの手法
多孔質媒体を通る多種流れの複合流体流れシミュレーションを完了する時に、履歴効果を捕捉することが有益である可能性がある。例えば、岩石層からの炭化水素抽出をシミュレートする時に、流体が地層に注入された時に移動しない炭化水素の捕捉されたポケット又は部分を決定するために履歴効果を含めることが有益である可能性がある。

この点に関して、多孔質媒体を通る例示的な多種流れは、実際のリザーバ岩石試料を通る多種多孔質媒体流れのコンピュータシミュレーションを含むことができる。例えば、地下で見出される炭化水素は、典型的には岩石層に存在する。これらの岩石層は、通常は部分的に多孔質であり、多孔質媒体として分類することができる。この多孔質媒体からの炭化水素抽出は、典型的には水のような炭化水素と混和しない流体を使用して実施される。多孔質媒体からの炭化水素抽出を理解するために、本明細書に説明するシステム及び方法は、多孔質媒体及び媒体を通る流れを特徴付ける。そのようなシミュレーションを実行する際には、システム内の異なる種の相対浸透率が重要である。絶対浸透率が単一スカラー値であり、細孔幾何学的形状に固有であるのに対して、多種多孔質媒体流れ解析における相対浸透率は、どの程度簡単に各種が基準種（多くの場合、水）の容積率（飽和）の関数として移動するかを示す結果の集合である。「相対的」という用語は、複数の種が存在する場合の種の浸透率が絶対浸透率によって正規化されることを示すために使用される。飽和値の重要な範囲にわたる各種の浸透率結果が得られた時に相対浸透率プロットがもたらされる。そのようなシミュレーションにおいて、多孔質媒体試料は、典型的にはＣＴ又はマイクロＣＴ走査過程を通してデジタル的に捕捉される。得られた画像の解析及び処理を使用して、流れ数値シミュレーションにおける使用のために多孔質の岩石幾何学的形状の３次元デジタル表現を生成することができる。その後に、多種流れを取り扱うことができる数値的方法を使用して、様々な飽和条件に対して試料を通る流れをシミュレートして相対浸透率値を含む流れ挙動の予想を取得することができる。「定常状態」手法において、各飽和レベルに関して、シミュレーションは、有効浸透率が収束するまで実行されることになる。収束後に、質量ソースシンクモデルをオンにして、シミュレーションは、望ましい飽和レベルがもたらされるまで実行されることになる。

本明細書に説明するシステム及び方法は、多孔質媒体を通る多種（例えば、少なくとも２つの不混和性成分）流れの相対浸透率の予想を提供する。一例示的数値的手法は、格子ボルツマン方式（ＬＢＭ）であり、この方法に関して、いくつかの異なる技術が、多種流れをシミュレートすることができるように通常の単一流体アルゴリズムを拡張する方法に対して存在する。これらの技術の１つは、Ｓｈａｎ−Ｃｈｅｎ相互作用力方法として公知である。一部の実装において、Ｓｈａｎ−Ｃｈｅｎベースの多種ＬＢＭ手法を使用して、多孔質媒体を通る多種流れをシミュレートして相対浸透率を含む流れ挙動の重要な態様を予想する。

図１は、行き止まり細孔が取り付けられた垂直細孔の図を示している。図１に示すように、デジタル多孔質媒体を通る多種流れの数値シミュレーションを実行する際に、シミュレーションを設定する１つの方法は、周期的ドメインを流れ方向に使用することである（例えば、周期的境界条件を含むことにより）。この方法では、流体が出入りする入口又は出口はなく、むしろ、領域の一端１２を出る流体は、他端１４で直ちに入る。従って、片側でシミュレーションセルを去るあらゆる物体／種は、反対側で入る。より具体的には、物体は、シミュレーション空間の出口端（例えば、１２）でボクセルを通過した時に、同じ速度及び方向で入口（例えば、１４）で対応するボクセルで再び現れる。流体をシミュレーション領域に通して移動するために、物体力が望ましい流れ方向に印加される。この物体力は、重力と類似のものである（かつ多くの場合にそのようにラベル付けされる）。周期的な設定条件の１つの潜在的な利点は、それが、流体が多孔質媒体とバルク流体領域の間のインタフェースを通って移動する時に発生する「端効果」を回避するということであり、これらの端効果は、シミュレーションからの予想浸透率に好ましくない影響を与える可能性がある。

上述のように、本明細書に説明する多相システムをモデル化する時に、飽和値の範囲にわたってシステム内の種の各々の相対浸透率に及ぼす影響を決定することが重要であり、その理由は、相対浸透率が、飽和値が異なると異なるからである。例えば、約０％の飽和、約２０％の飽和、約４０％の飽和、約６０％の飽和、及び約８０％の飽和に対して相対浸透率をシミュレート／決定することが望ましい可能性があるが、飽和レベルの選択は、岩と、同じくユーザが選択した過程、例えば、事前に固定されたＳｗ対適応Ｓｗとに依存する。ユーザは、可能な飽和レベルの一部だけをシミュレートするように選択する場合がある。

飽和度を修正する（例えば、システムにおいて様々な種の比率を修正する）１つの方法は、システム（例えば、非周期性システム）の入口に種の望ましい量の各々を単に追加することである。しかし、この方法は、周期的システムの安定性及び計算効率上の利点をもたらすものではなく、その理由は、物体力が、領域全体を最終状態に向けて同時に引っ張るからである。従って、飽和値を修正すると同時にシステムを周期的システムとしてシミュレートすることを可能にする方法が望ましいと考えられる。

飽和度を修正する別の方法は、複数の別々のシミュレーションを実行することである。各々は、望ましい飽和値で初期化される。図２Ａ〜図２Ｃは、異なる飽和度で実行される１組の重力駆動の周期的構成シミュレーションの図を示している（例えば、２０％を図２Ａに示し、４０％を図２Ｂに示し、５５％を図２Ｃに示す）。図２Ａ〜図２Ｃに示すシミュレーションに示す例示的システムにおいて、システムは、行き止まり細孔１８を含む。流体が入口１４から出口１２まで流れる時に、行き止まりポート１８に捕捉された流体は捕捉されたままである（例えば、移動しない）。この例において、異なる飽和度をシミュレートするために、各飽和値のシミュレーションは、同じ初期条件下でかつ他の飽和値でのシミュレーションの結果とは独立して実行されることになる。これは、１つのシミュレーションにおいて存在する異なる流体種が、どこで他の飽和レベルシミュレーション（飽和値の望ましい範囲に関する）のいずれにも影響を及ぼさないかに関する情報を意味する。その結果、履歴／ヒステレシスの影響は、相対浸透率予想において捕捉されていない。この例で分るように、履歴／ヒステレシスを欠くそのようなシミュレーションは、全体的にシステム内の流体と同じ飽和値を仮定すると、行き止まり細孔１８における流体をもたらす。しかし、そのような結果は、行き止まり細孔１８内の流体は捕捉されるので物理的に取得不可能である。

本明細書に説明するシステム及び方法は、異なる飽和度レベルでのシミュレーションが以前にシミュレートした飽和レベルからの履歴情報を保持する周期的システムを提供する。図３Ａ〜図３Ｃは、本明細書に説明する質量交換方法を使用して実施する異なる飽和度で実行する１組の重力駆動の周期的構成シミュレーションの図を示している（例えば、２０％を図３Ａに示し、４０％を図３Ｂに示し、５５％を図３Ｃに示す）。図３Ａ〜図３Ｃに示すシミュレーションに示す例示的システムにおいて、システムは、行き止まり細孔１８を含む。しかし、質量交換過程は、相対浸透率予想において履歴／ヒステレシス効果を捕捉し、行き止まり細孔１８に捕捉された流体は捕捉されたままであり（例えば、種は行き止まり細孔１８に存在するが動かない）、行き止まり細孔１８内の流体の相対的な飽和は変化しないようになっている。例えば、システムが８０％の油及び２０％の水で初期化された場合に、６０％の油及び４０％の水での第２のシミュレーションが実行された時に、細孔１８内の行き止まり内の油の百分率は８０％の初期化百分率のままであり、一方、システムの他の領域における油の百分率は、６０％よりも小さい値に低減され、全体的な飽和値は、６０％であるようになっている。物理的に意義がある方法で飽和度を修正するために、飽和値は、物理的に意義がある位置で第１の種（例えば、油）の一部を第２の種（例えば、水）で置換することによって修正される。

質量ソース／シンク方法は、ある一定の位置で異なる種間で流体質量を交換する。一部の例において、これは、一定の熱力学的圧力を維持しながら行われる。一部の例において、質量交換が発生する領域は、置換される種の流量に基づいて選択される。例えば、油及び水システムにおいて、システム内の最も対流的なゾーン又は閾値を超える流量を有するゾーンのような比較的高い流量がある領域を質量交換が発生する領域として選択することができる。交換領域が識別された状態で、識別された交換領域において、飽和レベルを修正するために油の一部が領域から除去されて水によって置換される。

油及び水を上述の第１及び第２の種として説明したが、本明細書に説明する方法は、あらゆる組の複数の異なる種に適用することができる。

多種多孔質媒体流れをシミュレートするこの手法は、マサチューセッツ州バーリントン所在のＥｘａ Ｃｏｒｐｏｒａｔｉｏｎから市販されるＰｏｗｅｒＦＬＯＷシステムのような格子ボルツマン方式（ＬＢＭ）に基づいて、時間明示的ＣＦＤ／ＣＡＡ解法に関連して使用することができる。巨視的連続体方程式の離散化に基づく方法と異なり、ＬＢＭは「メゾスコピック」ボルツマン動態方程式から開始して巨視的流体力学を予知する。航空音響学及び純粋音響学上の問題に対して得られる圧縮性非定常解は、様々な複雑流れの物理特性を予知するのに使用することができる。ＬＢＭベースのシミューレーションシステムを以下で全体的に説明し、次に、そのようなモデル化手法をサポートするために流体流れシミュレーションに関連して使用することができるスカラー解法手法に関して説明する。

Ｂ．モデルシミュレーション空間
ＬＢＭシステム物理過程シミューレーションシステムにおいて、流体流れは、１組の離散速度ｃ_iで評価される分布関数値ｆ_iによって表すことができる。分布関数の力学は、方程式４によって支配され、ここでｆ_i（０）は、平衡状態分布関数として公知であり、

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

この方程式は、分布関数ｆ_iの時間推移を説明する公知の格子ボルツマン方程式である。左側は、いわゆる「ストリーミング過程」に起因する分布の変化を表している。ストリーミング過程は、流体のポケットがグリッド位置で始まり、その後に、次のグリッド位置に速度ベクトルの１つに沿って移動する時である。その時点て、「衝突演算子」、すなわち、流体の開始ポケットに及ぼす流体の近くのポケットの影響を計算する。流体は、別のグリッド位置にのみ動くことができるので、全ての速度の全ての成分が共通の速度の倍数であるように速度ベクトルの適正な選択が必要である。

第１の方程式の右側は、流体のポケット間の衝突に起因する分布関数の変化を表す上述の「衝突演算子」である。ここに使用される衝突演算子の特定の形態は、Ｂｈａｔｎａｇａｒ、Ｇｒｏｓｓ、及びＫｒｏｏｋ（ＢＧＫ）によるものである。この形態は、強制的に分布関数を「平衡状態」の形態である第２の方程式によって与えられる規定の値にする。

このシミュレーションから、質量ρ及び流量速度ｕのような従来の流体変数は、方程式（３）における単純な合計として得られる。ここで、ｃ_i及びｗ_iの集合値は、ＬＢＭモデルを定義する。ＬＢＭモデルは、拡張可能なコンピュータプラットフォーム上で効率的に実行され、かつ時間非定常流及び複雑な境界条件に対して高い堅牢性で実行することができる。

ボルツマン方程式から流体システムに対する巨視的運動方程式を取得する標準的な手法は、完全なボルツマン方程式の連続した近似が取られるＣｈａｐｍａｎ−Ｅｎｓｋｏｇ方法である。

流体システムにおいて、密度の小さい外乱は音の速度で進む。ガスシステムにおいて、音の速度は、一般的に温度に依存する。流れにおける圧縮率の影響の重要度は、マッハ数として公知である固有速度及び音速の比率によって測定される。

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

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

１０１個の速度を含む３Ｄ−２モデル及び３７個の速度を含む２Ｄ−２モデルのようなより複合的なモデルを使用することができる。

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

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

上述のＬＢＭモデルは、２次元及び３次元の両方における流れの数値シミュレーションの効率的かつ堅牢な離散速度動態モデルの特定クラスを与える。この種類のモデルは、それらの速度に関連付けられた離散速度及び重みの特定の組を含む。速度は、離散速度モデル、特に格子ボルツマンモデルとして公知である種類の正確かつ効率的実施を容易にする速度空間内の直交座標のグリッド点と一致する。そのようなモデルを使用して、流れを高い忠実度でシミュレートすることができる。

図６を参照すると、物理過程シミューレーションシステムは、手順３００に従って作動し、流体流れのような物理過程をシミュレートする。シミュレーション前に、シミュレーション空間をボクセルの集合としてモデル化する（段階３０２）。典型的には、シミュレーション空間は、コンピュータ支援設計（ＣＡＤ）プログラムを使用して生成される。例えば、ＣＡＤプログラムを使用して風洞に位置決めされたマイクロデバイスを示すことができる。その後に、ＣＡＤプログラムによって生成されたデータを処理して、適切な分解能を有する格子構造を追加し、かつシミュレーション空間内の物体及び面に対処する。

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

物体の固有長は、物体の大規模特徴を表している。例えば、マイクロデバイスの周りの流れがシミュレートされる場合に、マイクロデバイスの高さは、固有長であると考えることができる。物体（例えば、自動車のサイドミラー）の小さい領域の周りの流れが関心領域である時に、シミュレーションの分解能は、増大させることができ、又は分解能が増大した領域は、関心領域の周りに使用することができる。ボクセルの寸法は、格子の分解能が増加する時に減少する。

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

状態の数は、各エネルギレベル内の可能な速度ベクトルの数によって決定される。速度ベクトルは、３次元ｘ、ｙ、及びｚを有する空間内の整数線速度から構成される。状態の数は、複数の種シミュレーションに対して増大される。

各状態ｉは、比エネルギレベル（すなわち、エネルギレベル０、１、又は２）での異なる速度ベクトルを表している。各状態の速度ｃ_iを、以下のように３次元の各々における「速度」に示している。
ｃ_i＝（ｃ_i,x，ｃ_i,y，ｃ_i,z）． 方程式（７）

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

３つのエネルギレベルの可能な置換の全てを生成することにより、合計３９個の可能な状態が得られる（１個のエネルギゼロ状態、６個のエネルギ１状態、８個のエネルギ３状態、６個のエネルギ４状態、１２個のエネルギ８状態、及び６個のエネルギ９状態）。

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

処理効率を得るために、ボクセルは、マイクロブロックと呼ばれる２ｘ２ｘ２容積にグループ分けされる。マイクロブロックは、ボクセルの並列処理を可能にし、かつデータ構造に関連付けられたオーバヘッドを最小にするように構成される。マイクロブロック内のボクセルの略記法は、Ｎｉ（ｎ）として定義され、ｎは、マイクロブロック内の格子部位の相対位置を表し、ｎ∈｛０，１，２．．．，７｝である。マイクロブロックを図７に示す。

図８Ａ及び図８Ｂを参照すると、面Ｓは、ファセットＦ_αの集合：
Ｓ＝｛Ｆ_α｝ 方程式（８）
としてシミュレーション空間（図８Ｂ）内で表現され、ここで、αは、特定のファセットを含む指標である。ファセットは、ボクセル境界に制限されず、典型的には、ファセットが比較的小さい数のボクセルに影響を与えるように、ファセットの近くのボクセルのサイズの程度であり、又はサイズよりも僅かに小さくサイズ設定される。特性が、面力学を実行するためにファセットに割り当てられる。特に、各ファセットＦ_αは、単位法線（ｎ_α）、面積（Ａ_α）、中心位置（ｘ_α）、及びファセットの面の動的特性を説明するファセット分布関数（ｆ_i（α））を有する。

図９を参照すると、分解能の異なるレベルをシミュレーション空間の異なる領域に使用して処理効率を改善することができる。典型的には、物体６５５の周りの領域６５０が、最大関心領域であり、従って、最高分解能でシミュレートされる。粘性の影響は物体からの距離と共に減少するので、分解能（すなわち、拡張されたボクセル容積）の減少するレベルを使用して、物体６５５から増大する距離で離間している領域６６０、６６５をシミュレートする。同様に、図１０に示すように、分解能のより低いレベルを使用して、物体７７５のより有意ではない特徴部の周りの領域７７０をシミュレートすることができ、一方、分解能の最高レベルを使用して、物体７７５の最も有意な特徴部（例えば、先端面及び後端面）の周りの領域７８０をシミュレートする。遠くの領域７８５は、分解能の最低レベル及び最も大きいボクセルを使用してシミュレートする。

Ｃ．ファセットによって影響を受けるボクセルの識別
図６を再び参照すると、シミュレーション空間がモデル化された状態で（段階３０２）、１つ又はそれよりも多くのファセットによって影響を受けるボクセルが識別される（段階３０４）。ボクセルは、いくつかの点で影響を受ける場合がある。最初に、１つ又はそれよりも多くのファセットによって交差されるボクセルは、ボクセルが非交差ボクセルに対して低減された容積を有するということにおいて影響を受ける。これは、ファセット及びファセットによって表される面の下にある材料がボクセルの一部を占有するので発生する。分画係数Ｐ_f（ｘ）は、ファセットによって影響されないボクセルの部分（すなわち、流れがシミュレートされている流体又は他の材料によって占有される可能性がある部分）を示している。非交差ボクセルに関して、Ｐ_f（ｘ）は、１に等しい。

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

図１１を参照すると、ゼロ以外の速度ベクトルｃ_iを有する各状態ｉに関して、ファセットＦ_αは、平行六面体Ｇ_iαの容積Ｖ_iαが、
Ｖ_iα＝｜ｃ_iｎ_α｜Ａ_α 方程式（９）
に等しいように、速度ベクトルｃ_i及びファセットの単位法線ｎ_αのベクトルドット積の大きさ（｜ｃ_iｎ_i｜）によって定義された高さ及びファセットの面積Ａ_αによって定義される基部を有する平行六面体Ｇ_iαによって定義された領域から粒子を受け入れるか、又は該領域に粒子を移動する。

ファセットＦ_αは、状態の速度ベクトルがファセットに向けられた時に（｜ｃ_iｎ_i｜＜０）、粒子を容積Ｖ_iαから受け入れ、状態の速度ベクトルがファセットから離れるように向けられた時に（｜ｃ_iｎ_i｜＞０）、粒子を領域に移動する。以下に説明するように、この式は、別のファセットが平行六面体Ｇ_iαの一部、すなわち、内部のコーナのような非凸面特徴部の近くに発生する可能性がある状態を占有した時には修正しなければならない。

ファセットＦ_αの平行六面体Ｇ_iαは、複数のボクセルの各部分又は全部に重なり合う場合がある。ボクセル又はその部分の数は、ボクセルのサイズに対するファセットのサイズ、状態のエネルギ、及び格子構造に対するファセットの向きに依存する。影響を受けるボクセルの数は、ファセットのサイズと共に増加する。従って、ファセットのサイズは、上述のように、典型的にはファセットの近くに位置するボクセルのサイズの程度か、又はそれよりも小さいように選択される。

平行六面体Ｇ_iαによって重ね合わされるボクセルＮ（ｘ）の部分は、Ｖ_iα（ｘ）として定義される。この項を使用して、ボクセルＮ（ｘ）とファセットＦ_α間で移動する状態ｉ粒子の流束Γ_iα（ｘ）は、ボクセル（Ｖ_iα（ｘ））との重なりの領域の容積ｘボクセル（Ｎ_i（ｘ））において状態ｉの密度粒子に等しい。
Γ_iα（ｘ）＝Ｎ_i（ｘ）Ｖ_iα（ｘ） 方程式（１０）

平行六面体Ｇ_iαが１つ又はそれよりも多くのファセットによって交差される時に、以下の状態が当て嵌まる。
Ｖ_iα＝ΣＶ_α（ｘ）＋ΣＶ_iα（β） 方程式（１１）
ここで、第１の合計は、Ｇ_iαによって重ね合わされた全てのボクセルに対処し、第２の項は、Ｇ_iαと交差する全てのファセットに対処する。平行六面体Ｇ_iαが別のファセットによって交差されていない時に、この式は、以下になる。
Ｖ_iα＝ΣＶ_iα（ｘ）． 方程式（１２）

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

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

他の２つの境界条件は、ファセットと相互作用する粒子の正味運動量に関連する。スリップ面と本明細書で呼ぶ表面摩擦のない面に関して、正味接線運動量流束は、ゼロに等しくなければならず、正味法線運動量流束は、ファセットでの局所圧力に等しくなければならない。従って、ファセットの法線ｎ_αに垂直である結合された受け入れられた運動量及び移動された運動量の成分（すなわち、接線成分）は等しくなければならず、ファセット（すなわち、法線成分）の法線ｎ_αに平行である結合された受け入れられた運動量及び移動送された運動量の成分間の差は、ファセットでの局所圧力に等しくなければならない。非スリップ面に関して、面の摩擦は、摩擦の量に関連する係数によってファセットによって受け入れられた粒子の結合された接線運動量に対するファセットによって移動された粒子の結合された接線運動量を低減する。

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

これから、ファセットＦ_α（ｃ_iｎ_α＜０）に向けられる各状態ｉに関して、ボクセルによってファセットＦ_αに供給される粒子数は、以下の通りである。

方程式（１４）

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

３．ファセットからファセットへの移動
次に、粒子をファセット間で移動する（段階３１０）。ファセットＦ_αの流入状態（ｃ_iｎ_α＜０）の平行六面体Ｇ_iαが別のファセットＦ_βによって交差される場合に、ファセットＦ_αによって受け入れられた状態ｉ粒子の一部は、ファセットＦ_βから到来することになる。特に、ファセットＦ_αは、以前の時間増分中にファセットＦ_βによって生成される状態_i粒子の一部を受け入れることになる。この関係を図１３に示すが、ファセットＦ_βによって交差される平行六面体Ｇ_iαの部分１０００は、ファセットＦ_αによって交差される平行六面体Ｇｉβの部分１００５に等しい。上述のように、交差された部分は、Ｖ_iα（β）として示している。この項を使用して、ファセットＦ_βとファセットＦ_αの間の状態_i粒子の流束は、以下として説明することができる。
Γ_iα（β，ｔ−１）＝Γ_i（β）Ｖ_iα（β）／Ｖ_iα， 方程式（１５）
ここで、 Γ_i （β，ｔ−１）は、前回の時間増分中にファセットＦ_βによって生成された状態ｉ粒子の尺度である。これから、ファセットＦ_α（ｃ_iｎ_α＜０）に向けられる各状態ｉに関して、他のファセットによってファセットＦ_αに供給される粒子数は、

方程式（１６）
であり、ファセットへの状態ｉ粒子の全流束は、以下の通りである。

方程式（１７）

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

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

並行した状態（ｃ_iｎ_α＝０）に関して、Ｖ_iα及びＶ_iα（ｘ）は、ゼロである。Ｎ_i（α）の式において、Ｖ_iα（ｘ）は、分子に現われ（Γ_iOTHER（α）の式から）、Ｖ_iαは、分母に現れる（Ｎ_i（α）の式から）。従って、並行した状態のＮ_i（α）をＶ_iα及びＶ_iα（ｘ）がゼロに近づく時にＮｉ（α）の限界として決定する。

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

４．ファセット面力学の実行
次に、各ファセットが上述の４つの境界条件を満たすために面力学を実行する（段階３１２）。ファセットの面力学を実行するための手順を図１４に示す。最初に、ファセットＦ_αに垂直な結合された運動量を以下としてファセットで粒子の結合された運動量Ｐ（α）を決定することによって決定する（段階１１０５）
全てのｉに関して、

方程式（２１）

これから、通常の運動量Ｐ_n(α)を以下として決定する。
Ｐ_n（α）＝ｎ_α・Ｐ（α）． 方程式（２２）

その後に、Ｎ_n-（α）を生成するために押し／引き技術を使用してこの法線運動量を排除する（段階１１１０）。この技術に従って粒子を法線運動量だけに影響を与える方法で状態間で移動する。押し／引き技術は、米国特許第５，５９４，６７１号明細書に説明されており、この特許は、引用によって本明細書に組み込まれている。

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

その後に、ファセットＦ_αの流出流束分布を流入流束分布及びボルツマン分布に基づいて決定する（段階１１２０）。最初に、流入流束分布Γ_i（α）とボルツマン分布との差を以下として決定する。
ΔΓ_i（α）＝Γ_iIN（α）−Ｎ_n-βi（α）Ｖ_iα． 方程式（２３）

この違いを使用して、流出流束分布は、以下の通りであり、ｎ_αｃ_i＞０に対して
Γ_iOUT（α）＝Ｎ_n-βi（α）Ｖ_iα−．Δ．Γ_i*（α） 方程式（２４）
ここで、ｉ^*は、状態ｉと反対の方向を有する状態である。例えば、状態ｉが（１、１、０、０）の場合に、状態ｉ^*は、（−１、−１、０、０）である。表面摩擦及び他の要素に適合するように、流出流束分布を以下に更に精緻化することができる。
ｎ_αｃ_i＞０に関して、

方程式（２５）

ここで、Ｃ_fは、表面摩擦の関数である、ｔ_iαは、ｎ_αに垂直である第１の接線ベクトルであり、ｔ_2αは、ｎ_α及びｔ_1αに垂直である第２の接線ベクトルであり、ΔＮ_j,1及びΔ_Nj,2は、状態ｉ及び表示される接線ベクトルのエネルギ（ｊ）に対応する分布関数である。分布関数を以下に従って決定する。

方程式（２６）

ここで、ｊは、エネルギレベル１状態に対して１、及びエネルギレベル２状態に対して２に等しい。

Γ_iOUT（α）の方程式の各項の関数は、以下の通りである。第１及び第２の項は、衝突がボルツマン分布を生成する際に有効であった範囲で法線運動量流束境界条件を実施するが、接線運動量流束偏差を含む。第４及び第５の項は、分離性効果又は不十分な衝突に起因する非ボルツマン構造に起因して生じる可能性があるこの偏差を補正する。最後に、第３の項は、面上の接線運動量流束の望ましい変化を実施するために指定の量の表面分画を追加する。摩擦係数Ｃ_fの生成を以下に説明する。ベクトル演算を伴う全ての項がシミュレーションを始める前に計算することができる形状因子であることに注意されたい。

これから、接線速度を以下として決定する。
ｕ_i（α）＝（Ｐ（α）−Ｐ_n（α）ｎ_α）／ρ， 方程式（２７）
ここで、ρは、ファセット分布の密度である。

方程式（２８）

前と同様に、流入流束分布とボルツマン分布の差を以下として決定する。
ΔΓ_i（α）＝Γ_iIN（α）−Ｎ_n-βi（α）Ｖ_iα． 方程式（２９）

その後に、流出流束分布は、以下になる。
Γ_iOUT（α）＝Ｎ_n-βi（α）Ｖ_iα−ΔΓ_i*（α）＋Ｃ_f（ｎ_αｃ_i）［Ｎ_n-βi*（α）−Ｎ_n-βi（α）］Ｖ_iα， 方程式（３０）
この方程式は、前回の技術によって求められた流出流束分布の最初の２つの方向に対応するが、変則的な接線流束の補正は不要である。

あらゆる手法を使用して、得られる流束分布は、運動量流束条件の全てを満たし、すなわち、以下のようになる。

方程式（３１）

ここで、ｐ_αは、ファセットＦ_αでの平衡圧であり、ファセットに粒子をもたらすボクセルの平均化された密度及び温度値に基づいて、ｕ_αは、ファセットでの平均速度である。

質量及びエネルギ境界条件が満たされることを保証するために、入力エネルギと出力エネルギの差を以下として各エネルギレベルｊに対して測定する。

方程式（３２）

ここで、指数ｊは、状態ｉのエネルギを示している。その後に、このエネルギ差を使用して差の項を生成する。
ｃ_jiｎ_α＞０に関して、

方程式（３３）

この差の項は、流束が以下になるように流出流束を修正するのに使用される。
ｃ_jiｎ_α＞０に対して
Γ_αjiOUTf＝Γ_αjiOUT＋δΓ_αji 方程式（３４）
この演算は、接線運動量流束を不変のままにすると同時に質量及びエネルギ流束を補正する。この調節は、流れがファセットの近傍においてほぼ均一かつ平衡状態に近い場合は小さい。得られた法線運動量流束は、調節の後に、近傍の不均一性又は非平衡特性に起因する補正を加えて近傍平均特性に基づいて平衡圧である値に僅かに変更される。

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

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

面と相互作用しないボクセルに関して、移動演算は、計算機的には全く簡単である。状態の母集団全体をあらゆる時間増分中に現在のボクセルからの目的地ボクセルに移動する。同時に、目的地ボクセルの粒子をボクセルから固有の目的地ボクセルに移動する。例えば、＋１ｘ及び＋１ｙ方向（１、０、０）に移動しているエネルギレベル１粒子は、現在のボクセルからｘ方向に＋１だけ上のかつ他の方向に対して０であるボクセルに移動する。粒子は、移動（１、０、０）の前に有していた同じ状態で目的地ボクセルで終わる。ボクセル内の相互作用は、他の粒子及び面との局所相互作用に基づいてその状態に対して粒子数を一部の場合に変えることになる。そうでない場合に、粒子は、同じ速度及び方向で格子に沿って移動し続けることになる。

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

方程式（３５）

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

６．ファセットからボクセルへの散乱
次に、流出粒子を各ファセットからボクセルに散乱させる（段階３１６）。本質的に、この段階は、粒子をボクセルからファセットに移動した集める段階の逆である。ファセットＦ_αからボクセルＮ（ｘ）に移動する状態ｉ粒子の数は、以下である。

方程式（３６）

ここで、Ｐ_f（ｘ）は、部分的なボクセルの容積低減に対処する。これから、各状態ｉに関して、ファセットからボクセルＮ（ｘ）に向けられる粒子の総数は、以下である。

方程式（３７）

ファセットからボクセルまで粒子を散乱させ、周囲のボクセルから流入した粒子と結合して結果を整数化した後に、特定のボクセルの特定の方向はアンダーフローする（負になる）か、又はオーバーフローする（８ビット実装において２５５を超える）可能性がある。それによってこれらの量が値の許された範囲に収まるように打ち切られた後に質量、運動量、及びエネルギの利得又は損失が結果的に発生することになる。そのような発生に対して保護するために、限度よりも大きい質量、運動量、及びエネルギが、問題のある状態の打切り前に蓄積される。状態が属するエネルギに関して、（アンダーフローのために）得られるか又は（オーバーフローのために）失われる値に等しい量の質量は、同じエネルギを有し、かつそれ自体オーバーフロー又はアンダーフローを受けないランダムに（連続的に）選択された状態に追加される。質量及びエネルギのこの追加から生じる追加の運動量は蓄積され、打切りからの運動量に追加される。同じエネルギ状態にのみ質量を追加することにより、質量カウンタがゼロに到達した時に質量及びエネルギが補正される。最後に、運動量アキュムレータがゼロに戻るまで運動量を押し／引き技術を使用して補正する。

７．流体力学の実行
最後に、流体力学を実行する（段階３１８）。この段階は、マイクロ力学又はボクセル内演算と呼ぶことができる。同様に、移流手順は、ボクセル間演算と呼ぶことができる。また、以下に説明するマイクロ力学演算を使用して、ファセットで粒子を衝突させてボルツマン分布を生成することができる。

流体力学は、ＢＧＫ衝突モデルとして公知である特定の衝突演算子によって格子ボルツマン方程式モデルにおいて保証される。この衝突モデルは、実際の流体システムにおける分布の力学を模倣している。衝突過程は、方程式１及び方程式２の右側によって十分に説明することができる。移流段階後に、流体システムの保存量、具体的には、密度、運動量、及びエネルギを方程式３を使用して分布関数から取得する。これらの量から、方程式（２）の

によって説明された平衡状態分布関数は、方程式（４）によって完全に指定される。方程式２と共に表１に両方とも列挙された速度ベクトルセットｃ_i、重みの選択により、巨視的挙動が正しい流体力学方程式に従うことが保証される。

Ｅ．可変分解能
図１５を参照すると、（図９及び１０に示して上述したような）可変分解能は、異なるサイズのボクセルを使用し、粗ボクセル１２０００及び精細ボクセル１２０５と以下で呼ぶ。（以下に説明する内容は、２つの異なるサイズを有するボクセルに言及する。説明する技術は、分解能の追加のレベルをもたらすためにボクセルの３つ又はそれよりも多くの異なるサイズに適用することができることを認められたい。）粗ボクセル及び精細ボクセルの領域間のインタフェースは、可変分解能（ＶＲ）インタフェース１２１０と呼ぶ。

可変分解能が面又は面の近くに使用される時に、ファセットは、ＶＲインタフェースの両側でボクセルと相互作用することができる。これらのファセットは、ＶＲインタフェースファセット１２１５（Ｆ_αIC）又はＶＲ精細ファセット１２２０（Ｆ_αIF）として分類される。ＶＲインタフェースファセット１２１５は、ＶＲインタフェースの粗側上に位置決めされ、かつ精細ボクセル内に延びる粗平行六面体１２２５を有するファセットである。（粗平行六面体は、ｃ_iが粗ボクセルの寸法に従って寸法決めされる平行六面体であり、一方、精細平行六面体は、ｃ_iが精細ボクセルの寸法に従って寸法決めされる平行六面体である。）ＶＲ精細ファセット１２２０は、ＶＲインタフェースの精細側上に位置決めされ、かつ粗ボクセル内に延びる精細平行四面体１２３０を有するファセットである。インタフェースファセットに関連する処理はまた、粗ファセット１２３５（Ｆ_αC）及び精細ファセット１２４０（Ｆ_αF）との相互作用を伴う場合がある。

ＶＲファセットの両方のタイプに関して、面力学は、精細規模で実行され、上述のように作動する。しかし、ＶＲファセットは、粒子がＶＲファセットへ及びＶＲファセットから移流する方法に関して他のファセットと異なっている。

ＶＲファセットとの相互作用は、図１６に示す可変分解能手順１３００を使用して取り扱う。この手順の殆どの段階は、非ＶＲファセットとの相互作用に対して上述の類似の段階を使用して実施される。手順１３００は、精細時間ステップに各々対応する２相を含む粗時間ステップ（すなわち、粗ボクセルに対応する期間）中に実行される。ファセット面力学は、各精細時間ステップ中に実行される。こういう理由から、ＶＲインタフェースファセットＦ_αICは、黒色ファセットＦ_αICb及び赤色ファセットＦ_αICrとそれぞれ呼ぶ２つの同一にサイズ設定されかつ向けられた精細ファセットと見なされる。黒色ファセットＦ_αICbは、粗時間ステップ内の第１の精細時間ステップに関連し、一方、赤色ファセットＦ_αICrは、粗時間ステップ内の第２の精細時間ステップに関連する。

最初に、粒子を第１の面から面への移流段階によってファセット間で移動する（流入させる）（段階１３０２）。粒子は、ファセットＦ_αから延び、かつファセットＦ_βの後部にある精細平行四面体（図１５、１２４５）の阻止されない部分を差し引いたファセットＦ_βの後部にある粗平行四面体（図１５、１２２５）の阻止されない部分の容積に対応するＶ_-αβの重み係数で黒色ファセットＦ_αICbから粗ファセットＦ_βCに移動される。精細ボクセルのｃ_iの大きさは、粗ボクセルのｃ_iの大きさの半分である。上述のように、ファセットＦ_αの平行六面体の容積は、以下として定義される。
Ｖ_iα＝｜ｃ_iｎ_α｜Ａ_α． 方程式（３８）

従って、ファセットの面積Ａ_αは粗平行六面体と精細平行六面体間では変化しないので、かつ単位法線ｎ_αが常に１の大きさを有するので、ファセットに対応する精細平行四面体の容積は、ファセットの対応する粗平行四面体の容積の半分である。

粒子は、ファセットＦ_αから延び、かつファセットＦ_βの後部にある精細平行四面体の阻止されない部分の容積に対応するＶ_-αβの重み係数で、粗ファセットＦ_αCから黒色ファセットＦ_βICbに移動される。

粒子は、Ｖ_αβの重み係数で赤色ファセットＦ_αICrから粗ファセットＦ_βCに、及びＶ_-αβの重み係数で粗ファセットＦ_αCから赤色ファセットＦ_βICrに移動される。

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

最後に、粒子は、同じ重み係数で精細ファセットＦ_αIF又はＦ_αFから他の精細ファセットＦ_βIF又はＦ_βFに、及びファセットＦ_αから延び、かつファセットＦ_βの後部にある粗平行四面体の阻止されない部分の容積に対応するＶ_Cαβの重み係数で粗ファセットＦ_αCから他の粗ファセットＦ_cに移動される。

粒子が面から面に移流された後に、粒子を第１の集める段階においてボクセルから集める（段階１３０４〜１３１０）。粒子を精細ファセットＦ_αＦに関して、精細平行四面体を使用して精細ボクセルから（段階１３０４）、粗ファセットＦ_αCに対して粗平行四面体を使用して粗ボクセルから集める（段階１３０６）。その後に、粒子を精細平行四面体を使用して粗ボクセル及び精細ボクセルの両方から黒色ファセットＦ_αIRbに関して、及びＶＲ精細ファセットＦ_αIFに対して集める（段階１３０８）。最後に、粒子を赤色ファセットＦ_αIRrに対して粗平行四面体と精細平行四面体との差を使用して粗ボクセルから集める（段階１３１０）。

次に、精細ボクセル、又はＶＲファセットと相互作用する粗ボクセルを精細ボクセルの集合に分割させる（段階１３１２）。単一粗時間ステップ内で精細ボクセルに粒子を伝達することになる粗ボクセルの状態が分割される。例えば、ファセットによって交差されていない粗ボクセルの適切な状態は、図７のマイクロブロックのように向けられる８つの精細ボクセルに分割される。１つ又はそれよりも多くのファセットによって交差された粗ボクセルの適切な状態は、ファセットによって交差されていない粗ボクセルの一部分に対応する完全な及び／又は部分的な精細ボクセルの集合に分割される。分割から生じる粗ボクセル及び精細ボクセルの粒子密度Ｎ_i（ｘ）は等しいが、精細ボクセルは、粗ボクセルの分画係数、及び他の精細ボクセルの分画係数と異なる分画係数Ｐ_fを有することができる。

その後に、面力学を精細ファセットＦ_αIF及びＦ_αFに対して（段階１３１４）、及び黒色ファセットＦ_αICbに対して実行する（段階１３１６）。力学を図１４に示して上述した手順を使用して実施する。

次に、粒子は、実際の精細ボクセル及び粗ボクセルの分割から生じた精細ボクセルを含む精細ボクセル間を移動する（段階１３１８）。粒子が移動された状態で、粒子を精細ファセットＦ_αIF及びＦ_αFから精細ボクセルに散乱させる（段階１３２０）。

また、粒子を黒色ファセットＦ_αICbから（粗ボクセルを分割することで生じる精細ボクセルを含む）精細ボクセルに散乱させる（段階１３２２）。粒子は、ボクセルが面の存在がない限りその時間に粒子を受け入れていた場合は精細ボクセルに散乱される。特に、粒子は、ボクセルが（粗ボクセルの分割から生じる精細ボクセルではなくて）実際の精細ボクセルである時に、ボクセルＮ（ｘ）を１速度単位よりも大きいボクセルＮ（ｘ＋ｃ_i）が実際の精細ボクセルである時に又はボクセルＮ（ｘ）を１速度単位を超えるボクセルＮ（ｘ＋ｃ_i）が粗ボクセルの分割から生じる精細ボクセルである時にボクセルＮ（ｘ）に散乱される。

最後に、流体力学を精細ボクセルに実行することによって第１の精細時間ステップを完了する（段階１３２４）。流体力学が実行されるボクセルは、粗ボクセルを分割することで生じる精細ボクセルを含まない（段階１３１２）。

手順１３００では、第２の精細時間ステップ中に類似した段階を実行する。最初に、粒子を第２の面から面移流段階において面から面を移動する（段階１３２６）。粒子は、黒色ファセットから赤色ファセットに、黒色ファセットから精細ファセットに、精細ファセットから赤色ファセットに、及び精細ファセットから精細ファセットに移流される。

粒子が面から面に移流された後に、粒子を第１の集める段階においてボクセルから集める（段階１３２８〜１３３０）。粒子を赤色ファセットＦ_αIRrに対して精細平行四面体を使用して精細ボクセルから集める（段階１３２８）。また、粒子を精細ファセットＦ_αF及びＦ_αFに対して精細平行四面体を使用して精細ボクセルから集める（段階１３３０）。

その後に、上述のように、面力学を粗ファセットＦ_αCに対して（段階１１３４）、精細ファセットＦ_αIF及びＦ_αFに対して（段階１３３２）、及び赤色ファセットＦ_αICrに対して実行する（段階１３３６）。

次に、粒子が、精細ボクセル及び粗ボクセルを表す精細ボクセルへ及びそこから移動されるように、粒子を精細分解能を使用してボクセル間で移動する（段階１３３８）。その後に、粒子が粗ボクセルへ及びそこから移動されるように、粒子を粗分解能（段階１３４０）を使用してボクセル間で移動する。

次に、組み合わせた段階において、粗ボクセルを表す精細ボクセル（すなわち、粗ボクセルを分割させることで生じる精細ボクセル）が粗ボクセルに組み込まれる間に粒子をファセットからボクセルに散乱させる（段階１３４２）。この組み合わせた段階において、粒子は、粗平行四面体を使用して粗ファセットから粗ボクセルに、精細平行四面体を使用して精細ファセットから精細ボクセルに、精細平行四面体を使用して赤色ファセットから精細又は粗ボクセルに、及び粗平行四面体と検索平行六面体間の差を使用して黒色ファセットから粗ボクセルに散乱される。最後に、流体力学を精細ボクセル及び粗ボクセルに実行する（段階１３４４）。

Ｆ．質量交換
上述のように、様々なタイプのＬＢＭは、多孔質媒体を通る多種流れの背景として機能する流体流れの解法に適用することができる。一部のシステムにおいて、相対浸透率予想をもたらすために、１組の重力駆動周期的構成シミュレーションを流体質量が特定の位置で異なる種の間で交換される質量シンク／ソース法を使用して異なる飽和度で実行することができる。この方法では、飽和度を変えることができ、一方、以前の飽和条件からの流体種の空間分布の多くは不変のままであり、従って、１つの飽和条件でのシミュレーション結果からの情報は、別の飽和条件でのシミュレーションに影響を与える。例えば、シミュレーションが一時点にわたって進行し、その後に、油を特定の位置で水と交換される場合に、それによって飽和値が変化し、一方、全ての他の位置での油及び水の分布は保存される。より具体的には、各水飽和率で、シミュレーションは、有効浸透率が収束するまで実行されることになる。収束後に、質量交換が始まることになり、次の望ましい水飽和率レベルがもたらされるまでシミュレーションは実行されることになる。望ましい水飽和率レベルがもたらされた状態で、最終の交換はオフにされることになり、シミュレーションは、有効浸透率が収束するまで実行されることになる。従って、システムは、質量交換シミュレーションと収束シミュレーションの間で交替する。

図１７は、多孔質媒体を通る多種流れのシミュレーション中にシステムに対して飽和値を修正するために質量交換過程を使用する例示的な過程の流れ図である。過程は、様々な飽和値に対して相対浸透率予想を生成する。相対浸透率曲線は、炭化水素リザーバの生産性を推定する重要な入力である。２つの終点、すなわち、非還元性水飽和率（Ｓｗｉ）及び残留物油飽和率（Ｓｏｒ）により、（ＳｏｒとＳｗｉとの飽和値の差を取ることによって）生成することができる油の量、並びに抽出することができない（Ｓｏｒ）油の量の推定が可能である。相対浸透率曲線の形状及び水と油の相対浸透率曲線の関係は、水が存在する場合に油を移動することがいかに容易であるかに関する知識をもたらし、かつ石油及びガス業界を通して使用されるリザーバモデル化システムの極めて重要な入力である。

過程は、多孔質媒体に関する情報を取得して（１４１０）多孔質媒体の３Ｄデジタル表現を生成する（１４１２）ことで始まる。例えば、多孔質媒体試料は、典型的にはＣＴ又はマイクロＣＴ走査過程を通してデジタル的に捕捉することができる。得られた画像の解析及び処理を使用して、流れ数値シミュレーションにおける使用に対して多孔質媒体幾何学的形状（例えば、多孔質岩）の３次元デジタル表現を生成することができる。

その後に、コンピュータシステムは、重力駆動周期的システムを初期化する（１４１４）。重力駆動周期的システムにおいて、物体力を望ましい流れ方向に印加してシミュレーションを通して流体を移動する。この方法では、流体が出入りする入口又は出口はなく、むしろ、領域の一端を出る流体は、他端で直ちに入る。システムはまた、特定の飽和レベルに又は特定のＳｗに対応する特定の流体分布に初期化される。入口及び出口を接続するために、岩幾何学的形状は、流れ方向にミラーリングされる。これは、出口での細孔プロファイルが入口の幾何学的形状と同一であることを意味する。浸潤タイプの相対浸透率計算を開始する２つの方法がある。この説明において、両方の方法では、水及び油が２つの不混和性の流体であり、岩は水で濡れていると仮定する。第１の方法において、細孔全体を油（Ｓｗ＝０％）で初期化し、Ｓｗが区分的に増加されるように水を質量シンク／ソース交換機構を使用して岩内に導入する。その後に、第１の水流量がゼロ以外になった時に非還元性水飽和率を飽和Ｓｗとして決定する。より正確な手法は、排水シミュレーションを実行することである。そのようなシミュレーションにおいて、岩を水（Ｓｗ＝１００％）で最初に完全飽和させる。その後に、終圧がリザーバ条件下の圧力に対応するように油を区分的に増加する圧力で片側で細孔に押し込む。次の圧力増分が印加される前に、入口での圧力は油侵入が停止する（収束する）ように小さい増分で増大させる。そのような排水試験の結果が、毛管圧力と対応する水飽和率との関係である。最大圧力での水飽和率は、非還元性水飽和率Ｓｗｉであると考えて、その後に、相対浸透率試験の開始点として使用することができる。水成分の分布を相対浸透率試験の初期水分分布として使用することができる。

システムの初期化後に、定常状態が得られるまでシミュレーションを実行する（１４１６）。シミュレーションは、典型的には重要な飽和領域の１つの極の飽和値から開始する。この飽和極値で、典型的には、単一種のみが、多孔質媒体を通って流れ、一方、他の種は存在するが動かない（捕捉される）。この状態のシミュレーションが安定した状態に到達した時に相対浸透率を決定する（１４１８）。各成分の相対浸透率を成分の有効浸透率の比率÷絶対浸透率

として計算し、αは、水又は油成分を示している。絶対浸透率は、同じ岩幾何学的形状の単相流れシミュレーションにおいて決定することができ、又は非還元性水飽和率Ｓｗｉでの有効浸透率として定義することができる。
有効浸透率の方程式は、以下のように書くことができる。

動的粘性を動粘性と密度の積

で置換し、流束（flux）を細孔（pore）内の測定された平均速度と多孔率（porosity）の積

で置換すると、以下になる。

重力駆動周期的システムに関して、圧力水頭（pressure head）

を重力と密度の積

で置換することができる。この式を有効浸透率の方程式に代入すると、以下になる。

この方程式をシミュレーション領域内の全てのボクセルに対して計算し、その後に、シミュレーション領域全体にわたって平均化する。

システムは、その後に、シミュレートされる追加の飽和レベルがあるか否かを決定する（１４２０）。シミュレートされる追加の飽和レベルがある場合に、システムは、質量が種間で交換される過程を始める１４２２。この過程において、各質量交換位置で、１つの種を除去して、他の種の等しい質量を注入する。この質量交換過程では、飽和率を開始値から新しい値に変える。そのような交換過程の１つの例を図１８に関してより詳細に説明する。一般的に、質量交換は、次の望ましい飽和条件をもたらすのに十分な容積にわたって適切な高い流量位置で行われる（１４２４）。飽和条件が満たされた時に、質量交換をオフにし（１４２６）、再び定常状態に到達するまでシミュレーションは続ける（１４１６）。

飽和率は、各時間ステップで更新されず、むしろ、収束に到達するのに十分な時間ステップが実行された後にのみ更新される。この過程は、望ましい１組の飽和値にわたって繰り返される。このようにして、本明細書に説明するシステム及び方法は、飽和値間に流体種分布の影響を取り込みながら飽和値の範囲全体を横切ることを可能にする。このようにして、履歴効果が相対浸透率結果に捕捉され、シミュレーション手法は、対応する物理試験状態において起こることをより正確に再現することができる。相対浸透率を予想するシミュレーションが、低コスト化、かつ物理試験よりも時間短縮になる可能性を前提として、本発明の質量シンク／ソース過程を使用する多種多孔質媒体シミュレーション方法は、石油炭鉱及び生産業界に価値があると予想される。

油／水シミュレーションの１つの特定の例において、質量交換は、以下の方程式に基づくことができる。

油成分Ｓ_oilのシンク項は、１／時間ステップという単位を有する交換率αの関数、臨界速度を制御するＨｅａｖｉｓｉｄｅ関数、Ａｔｗｏｏｄ数判断基準を使用して油成分に属するボクセルを識別するＨｅａｖｉｓｉｄｅ関数、及び油成分の密度である。大きい交換率αは、結果的に時間ステップ当たりの交換された質量のより大きい画分になることになる。速度ベクトルに適用されるＨｅａｖｉｓｉｄｅ関数は、特定の臨界速度ｕ₀が対流ゾーンの一部である容積として特定の液量を識別するために超えられるべきであることを保証する。Ａｔｗｏｏｄ数は、成分を識別するのに使用されるパラメータである。この数は、水及び油に対応する−１と＋１の間で変動する。第２のＨｅａｖｉｓｉｄｅ関数は、交換に考慮される液量がほとんど油で満たされることを保証する。Ａｔｗｏｏｄ数Ａ₀の閾値は、典型的には−１の近くにあり、例えば、０．９５よりも小さい。水成分Ｓ_waterの質量ソース項の大きさは、油成分ｘ−１の質量シンク項の大きさに等しく設定される。しかし、一部の例において、ソース項の質量は、これに代えて、例えば、単位容積当たりの定圧を維持する異なる制約によって決定することができる。

シミュレートされた多孔質岩石試料における質量交換の位置の選択が重要である。物理定常相対浸透率実験において、飽和率は、流体の対応する組合せを多孔質媒体に注入することによって設定される（すなわち、望ましい飽和値をもたらすことになる混合物比率を用いて）。例えば、浸潤過程と呼ばれるものにおいて、最初に、水の方が容積率が高い混合物比率を有する水油混合物を油で殆どは満たされる岩石試料に注入する。この注入により、いくつかの油で満たされた領域を水で置換することによって多孔質媒体内の水の飽和が増大する。この物理状況において、最高局所流量を有する多孔質媒体領域は、水に取って代わられる油を有する第１の領域であり、入れ替えはまた、時間と共により低い局所流量を有する領域において起こる。

図１８は、質量交換に向けて位置を決定し、ターゲット飽和レベルに到達するために質量交換を実行する例示的な過程を示している。シミュレーションにおいて上述の物理過程の影響を再現するために、交換の領域（質量シンク／ソースが適用される場所）は、置換される種の比較的高い流量がある領域（例えば、上述の例における高い油流れの領域）であるように選択される。より具体的には、システムは、ボクセルの速度ベクトルに関する情報を取得し（１５１０）、速度ベクトル情報に基づいて質量交換のためにボクセルを識別する（１５１２）。１つの例において、臨界速度値は、システムによって記憶することができ、閾値を超える速度を有するあらゆるボクセルを質量交換に対して選択することができる。別の例において、１組の速度ベクトルをランク付けすることができ、最高速度値を有するボクセルの定められた百分率（例えば、約１０％）を選択することができる。別の例において、交換率は、局所速度の関数として定義することができ、例えば、速く移動する油小塊は、遅く移動する油小塊と比較してより速く交換されることになる。従って、システムは速度ベクトルに基づいて対流ゾーンを識別し、識別された対流ゾーンにおいて質量交換を実行する。質量交換は、システムの非対流ゾーン内で実行されない。動かない（例えば、細孔内に捕捉されたか又はアクセス不能）流体は置換されず、その理由は、そのような領域における速度は、低いか又はゼロの近くになるからである。

交換領域が識別された状態で、低減されるように指定された種は、その領域から除去され、増大されるように指定された種の１つ（又はそれよりも多く）によって入れ替わる（１５１４）。交換領域において、質量交換は、システムにおいて増大されるよう望ましい種を入力するソース、及び低減されるように望ましい種を除去するシンクをもたらすことによって実行される。等しい質量交換の場合に、ソースは、シンクによってボクセルから除去される第２の種の質量に等しい第１の種の質量を入力する。上述のように、流体密度を全ての密度分布ｆ_i

にわたる和のように計算する。全ての伝播後の密度分布を以下の方法でスケーリングする。

Ｓ_oilが負であり、これは、油の密度が浸潤過程において低減されることを意味することに注意されたい。平均水密度が平均油密度と同一であると仮定して、水成分の密度分布の組を同様に修正する。

例えば、油及び水システムにおいて、水は、油に取って代わる。交換は、複数の時間ステップにわたって行われ、ボクセルの質量の小さい部分（例えば、０．１〜０．３％、０．０１％〜０．０３％）のみが、単一時間ステップ中に変更される。他の例においても、より大きい交換率が可能である。従って、１０％だけシステムの飽和率を修正するのに、シミュレーションにおいて１０００回超の時間ステップが掛かる可能性がある。質量の一部を入れ替えた後に、システムは、流体流れシミュレーションの１つ又はそれよりも多くの時間ステップを実行し（１５１６）、ターゲット飽和値に到達したか否かを決定する。到達した場合に、システムは、質量交換を終了する（１５２０）。到達していない場合に、システムは、段階１５１４で種を交換することに戻る。別の例において、システムは、交換過程全体にわたって飽和度をモニタし、いつターゲット飽和値に到達したかを検出することになる。

上述の例において、システムは、識別されたボクセルの各々において同量の質量を入れ替えたが、一部の例において、交換量は、ボクセル当たりで速度ベクトルに基づくことができる。例えば、質量は、より大きい速度ベクトルを有する領域においてより大きい割合で交換することができる。一部の例において、質量交換率は、ボクセルに対して速度ベクトルに正比例する。

上述の例において、速度ベクトルを使用して質量交換の対象である対流ゾーンを決定したが、いくつかの追加の例においてシステム／ボクセル内の成分の移動を示す他の値を使用することができる。速度の代わりに、位置判断基準を交換過程を可能にする代替判断基準として導入することができる。そのような位置情報は、別のシミュレーションから又は孔径分布解析から取得することができる（例えば、浸出経路内の流体のみを交換する）。

実施形態は、デジタル電子回路、又はコンピュータハードウエア、ファームウエア、ソフトウエア、又はその組合せに実施することができる。本明細書に説明する技術の装置は、プログラマブルプロセッサによる実行のために機械可読媒体（例えば、ストレージデバイス）に有形に具現化又は記憶されたコンピュータプログラム製品に実施することができ、方法アクションは、入力データに対して演算して出力を生成することにより、本明細書に説明する技術の演算を実行する命令のプログラムを実行するプログラマブルプロセッサによって実行することができる。本明細書に説明する技術は、データ及び命令をデータストレージシステム、少なくとも１つの入力デバイス、及び少なくとも１つの出力デバイスから受信して、データ及び命令をそれらに送信するように結合された少なくとも１つのプログラマブルプロセッサを含むプログラマブルシステム上で実行可能である１つ又はそれよりも多くのコンピュータプログラムで実行することができる。各コンピュータプログラムは、高レベル手順又はオブジェクト指向プログラミング言語、又は必要に応じてアセンブリ又はマシン語で実行することができ、いずれにしても、言語は、コンパイル又は解釈された言語とすることができる。適切なプロセッサには、一例として、汎用及び専用マイクロプロセッサの両方が挙げられる。一般的に、プロセッサは、命令及びデータを読取専用メモリ及び／又はランダムアクセスメモリから受信することになる。一般的に、コンピュータは、データファイルを記憶する１つ又はそれよりも多くの大容量ストレージデバイスを含むことになり、そのようなデバイスには、内蔵ハードディスク及び取外し可能ディスクのような磁気ディスク、光磁気ディスク、光ディスクが挙げられる。コンピュータプログラム命令及びデータを有形に具現化するのに適するストレージデバイスには、一例として、ＥＰＲＯＭ、ＥＥＰＲＯＭ、及びフラッシュメモリデバイスのような半導体メモリデバイス、内部ハードディスク及び着脱式ディスクのような磁気ディスク、光磁気ディスク、及びＣＤ−ＲＯＭディスクを含む全ての形態の不揮発性メモリが挙げられる。上述のいずれも、ＡＳＩＣ（特定用途向け集積回路）によって補足されるか、又はＡＳＩＣに組み込むことができる。

いくつかの実施を説明した。それにも関わらず、本発明の精神及び範囲から逸脱することなく様々な修正を行うことができることは理解されるであろう。従って、他の実施も以下の特許請求の範囲内である。

１２ 領域の一端
１４ 領域の他端
１８ 行き止まり細孔

Claims (24)

1. 非一時的コンピュータ可読媒体に有形に具現化され、実行された時に物理過程流体流れをシミュレートし、シミュレートされる前記物理過程流体流れのターゲット特性を決定する命令を記録したコンピュータプログラム製品であって、前記コンピュータプログラム製品は、コンピュータをして、
   多孔質媒体のデジタル表現を受信し、
   前記多孔質媒体の前記デジタル表現において多種流体の活動をシミュレートし、前記多種流体は第１の流体種（第１の種）及び第２の流体種（第２の種）を含み、前記多種流体の前記活動は前記多孔質媒体内の前記多種流体の要素の動きをモデル化し、
   前記多孔質媒体内の前記第１の種の流体の閾値を超える流量がある領域を交換領域として識別し、
   前記識別された交換領域において、前記第１の種の第１の質量を除去し、前記第２の種の第２の質量で置換して、前記識別された交換領域における前記多孔質媒体内の前記第２の種に対する前記第１の種の特性値を修正し、
   前記修正された多種流体の活動をシミュレートし、かつ
   前記修正された特性値に基づいて前記シミュレートされた物理過程の前記ターゲット特性を決定する、
   ように構成される、
   ことを特徴とするコンピュータプログラム製品。
2. 前記多孔質媒体内の前記第１の種の移動を示す値に基づいて前記多孔質媒体内の前記交換領域を識別する命令を更に含むことを特徴とする請求項１に記載のコンピュータプログラム製品。
3. 前記多孔質媒体の前記デジタル表現は、１組のボクセルとして表わされた格子構造を含むシミュレーション空間でモデル化され、前記プログラムは、
   前記多孔質媒体内の１組の対流ボクセルを識別するためのコンピュータ命令の１つ以上の構成を含み、前記コンピュータ命令の１つ以上の構成は、
   前記第１の種の高い流量を有するボクセルを識別するための構成、
   速度が閾値を超えるボクセルを識別するための構成、
   ボクセルをそれらに関連付けられた速度に基づいてランク付けし、かつ最大速度を有する前記ボクセルの一部分を選択するための構成、の１つ以上を含むことを特徴とする請求項１に記載のコンピュータプログラム製品。
4. 修正された前記特性は、前記第１の種の圧力であり、前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の飽和値であり、前記製品はさらに、
   前記第１の種の圧力及び前記第２の種の前記飽和値から前記ターゲット飽和値に到達したか否かを決定するための命令をさらに含む請求項１に記載のコンピュータプログラム製品。
5. 修正された前記特性は、前記第１の種の浸透率であり、前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の相対浸透率であり、前記製品はさらに、
   前記ターゲット相対浸透率に到達したか否かを決定するための命令をさらに含む請求項１に記載のコンピュータプログラム製品。
6. 修正された前記特性は、前記第２の種に対する前記第１の種の飽和値であり、
   前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の飽和値であり、前記製品はさらに、
   前記ターゲット飽和値に到達したか否かを決定するための命令をさらに含む請求項１に記載のコンピュータプログラム製品。
7. 修正された前記特性は、重力に関する圧力値及び前記第２の種に対する前記第１の種の密度であり、
   前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第１の種の圧力水頭の値であり、前記製品はさらに、
   前記ターゲット圧力水頭の値に到達したか否かを決定するための命令をさらに含む請求項１に記載のコンピュータプログラム製品。
8. 前記ターゲット特性値に到達していなかったときに前記修正された多種の流体の活動をシミュレートするための命令をさらに含む請求項１に記載のコンピュータプログラム製品。
9. 重力駆動シミュレーションを使用して前記多種の流体の前記活動をシミュレートするための命令をさらに含む請求項１に記載のコンピュータプログラム製品。
10. 前記修正された飽和値に基づいて前記多孔質媒体中の前記第１の種の相対浸透率を決定するための命令をさらに含む請求項１に記載のコンピュータプログラム製品。
11. 物理過程流体流れをシミュレートするためのコンピュータシステムであって、前記システムは、
    １つ以上のプロセッサ装置、
    前記１つ以上のプロセッサ装置に接続されたメモリ、
    コンピュータプログラムを記憶する１つ以上の記憶装置であって、前記コンピュータプログラムは、実行されたときに、前記システムが物理過程流体流れをシミュレートし、シミュレートされる前記物理過程流体流れのターゲット特性を決定する、命令を含むものであり、前記コンピュータプログラムは、コンピュータをして、
    多孔質媒体のデジタル表現を受信し、
    前記多孔質媒体の前記デジタル表現において多種流体の活動をシミュレートし、前記多種流体は第１の流体種（第１の種）及び第２の流体種（第２の種）を含み、前記多種流体の前記活動は前記多孔質媒体内の前記多種流体の要素の動きをモデル化し、
    前記多孔質媒体内の前記第１の種の流体の閾値を超える流量がある領域を交換領域として識別し、
    前記識別された交換領域において、前記第１の種の第１の質量を除去し、前記第２の種の第２の質量で置換して、前記識別された交換領域における前記多孔質媒体内の前記第２の種に対する前記第１の種の特性値を修正し、
    前記修正された多種流体の活動をシミュレートし、かつ
    前記修正された特性値に基づいて前記シミュレートされた物理過程の前記ターゲット特性を決定する、
    ように構成される、
    ことを特徴とするシステム。
12. 前記コンピュータプログラムは、前記多孔質媒体内の前記第１の種の移動を示す値に基づいて前記多孔質媒体内の前記交換領域を識別する命令を更に含むことを特徴とする請求項１１に記載のシステム。
13. 前記多孔質媒体の前記デジタル表現は、１組のボクセルとして表わされた格子構造を含むシミュレーション空間でモデル化され、前記プログラムは、
    前記多孔質媒体内の１組の対流ボクセルを識別するためのコンピュータ命令の１つ以上の構成を含み、前記コンピュータ命令の１つ以上の構成は、
    前記第１の種の高い流量を有するボクセルを識別するための構成、
    速度が閾値を超えるボクセルを識別するための構成、
    ボクセルをそれらに関連付けられた速度に基づいてランク付けし、かつ最大速度を有する前記ボクセルの一部分を選択するための構成、の１つ以上を含むことを特徴とする請求項１１に記載のシステム。
14. 修正された前記特性は、前記第１の種の圧力であり、前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の飽和値であり、前記システムはさらに、
    前記第１の種の圧力及び前記第２の種の前記飽和値から前記ターゲット飽和値に到達したか否かを決定するための命令をさらに含む請求項１１に記載のシステム。
15. 修正された前記特性は、前記第１の種の浸透率であり、前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の相対浸透率であり、前記システムはさらに、
    前記ターゲット相対浸透率に到達したか否かを決定するための命令をさらに含む請求項１１に記載のシステム。
16. 修正された前記特性は、前記第２の種に対する前記第１の種の飽和値であり、
    前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の飽和値であり、前記システムはさらに、
    前記ターゲット飽和値に到達したか否かを決定するための命令をさらに含む請求項１１に記載のシステム。
17. 修正された前記特性は、重力に関する圧力値及び前記第２の種に対する前記第１の種の密度であり、
    前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第１の種の圧力水頭の値であり、前記システムはさらに、
    前記ターゲット圧力水頭の値に到達したか否かを決定するための命令をさらに含む請求項１１に記載のシステム。
18. 物理過程流体流れをシミュレートし、シミュレートされる前記物理過程流体流れのターゲット特性を決定する方法であって、前記方法は、
    多孔質媒体のデジタル表現を受信する段階と、
    前記多孔質媒体の前記デジタル表現において多種流体の活動をシミュレートさせ、前記多種流体は第１の流体種（第１の種）及び第２の流体種（第２の種）を含み、前記多種流体の前記活動は前記多孔質媒体内の前記多種流体の要素の動きをモデル化する段階と、
    前記多孔質媒体内の前記第１の種の流体の閾値を超える流量がある領域を交換領域として識別する段階と、
    前記識別された交換領域において、前記第１の種の第１の質量を除去し、前記第２の種の第２の質量で置換して、前記識別された交換領域における前記多孔質媒体内の前記第２の種に対する前記第１の種の特性値を修正する段階と、
    前記修正された多種流体の活動をシミュレートする段階と、かつ
    前記修正された特性値に基づいて前記シミュレートされた物理過程の前記ターゲット特性を決定する段階と、
    を含む方法。
19. 前記多孔質媒体内の前記第１の種の移動を示す値に基づいて前記多孔質媒体内の前記交換領域を識別する段階を更に含む請求項１８に記載の方法。
20. 前記多孔質媒体の前記デジタル表現は、１組のボクセルとして表わされた格子構造を含むシミュレーション空間でモデル化され、前記方法は、
    前記多孔質媒体内の１組の対流ボクセルを識別する段階を含み、識別する段階は、
    前記第１の種の高い流量を有するボクセル、
    速度が閾値を超えるボクセル、
    ボクセルをそれらに関連付けられた速度に基づいてランク付けすること、の１つ以上を含むことを特徴とする請求項１８に記載の方法。
21. 修正された前記特性は、前記第１の種の圧力であり、前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の飽和値であり、前記方法はさらに、
    前記第１の種の圧力及び前記第２の種の前記飽和値から前記ターゲット飽和値に到達したか否かを決定する段階をさらに含む請求項１８に記載の方法。
22. 修正された前記特性は、前記第１の種の浸透率であり、前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の相対浸透率であり、前記方法はさらに、
    前記ターゲット相対浸透率に到達したか否かを決定する段階をさらに含む請求項１８に記載の方法。
23. 修正された前記特性は、前記第２の種に対する前記第１の種の飽和値であり、
    前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第２の種に対する前記第１の種の飽和値であり、前記方法はさらに、
    前記ターゲット飽和値に到達したか否かを決定する段階をさらに含む請求項１８に記載の方法。
24. 修正された前記特性は、重力に関する圧力値及び前記第２の種に対する前記第１の種の密度であり、
    前記シミュレートされた物理過程流体流れの前記ターゲット特性は、前記第１の種の圧力水頭の値であり、前記方法はさらに、
    前記ターゲット圧力水頭の値に到達したか否かを決定する段階をさらに含む請求項１８に記載の方法。

Images