JP6967176B1 - 非対角優位不定係数行列のための圧力ソルバーを用いた貯留層シミュレーション - Google Patents

Abstract

コンピュータの性能は、コンピュータ処理による貯留層シミュレーション中に、貯留層モデルのグリッドセル間の圧力分布の決定中に改善される。収束は、条件が処理中に含まれる係数行列を不定にするときに得ることが困難であることが分かる。不定係数行列は、蒸気液体平衡に関連する貯留層内の物理的条件のために、または不適切な導関数により数値的に生成された非物理的条件のために、発生する可能性がある。貯留層シミュレーション中に収束が困難になると、時間ステップ切断の従来の修正措置が回避される。時間ステップ切断は、数百万、数十億またはそれ以上の数の未知のパラメータ値を含むモデルを有する非常に大きな貯留層に対して、コンピュータの使用および時間の点で非常にコストがかかることが分かっている。

Description

本発明は貯留層シミュレーション、特に、流体の生産中の貯留層内の流体流および圧力分布の両方を考慮に入れたシミュレーション方法論を含む。

石油・ガス産業では、探鉱・生産目的のコンピュータ化されたシミュレーション、モデリング、分析のために大量のデータを処理する必要がある。例えば、地下炭化水素貯留層の開発は、典型的には貯留層のコンピュータシミュレーションモデルの開発および分析を含む。貯留層、およびその流体の存在の現実的なシミュレーションモデルは、炭化水素貯留層からの最適な将来の石油およびガスの回収を予測するのに役立つ。石油・ガス会社は、石油貯蔵を開発する能力を強化するための重要なツールとして、シミュレーションモデルに依存するようになってきた。

地下炭化水素貯留層は、典型的には石油流体混合物および水の両方を含有する複雑な岩石層である。貯留層流体内容物は通常、２つ以上の流体相で存在する。貯留層流体中の石油混合物は、これらの岩石層に掘削され、仕上げられた油井によって生成される。時には、石油流体の回収率を向上するために、水またはガスのいずれかまたは両方のような流体も、これらの岩石層に注入される。

貯留層シミュレーションは、多孔質体シミュレーションにおけるフローの一般領域に属する。しかしながら、貯留層シミュレーションは通常、高圧および高温下にある地下地層中の複数の炭化水素成分および複数の流体相を含む。これらの炭化水素流体および含まれる地下水の化学相挙動は、これらのシミュレータにおいて考慮されなければならない。

シミュレーションモデルは、岩石層および油井の特定の幾何学的形状を記述する体積データ、流体および岩石特性などの貯留層特性データ、ならびに問題の油田またはガス田の特定の貯留層に関する生産および注入履歴を含む。シミュレーションモデルは、データ処理システムＳ上で実行される一連のコンピュータプログラムであるシミュレータ（貯留層シミュレータとして知られる）によって形成される。

これらのモデルを実行する貯留層シミュレータは、コンピュータにより実施される数値的方法論、または基礎となる数学的モデルのコード化されたアルゴリズムおよびデータ構成である。これらの炭化水素貯留層内の流体移動の物理的現象を表す数学的モデルは、流体の生産注入によって誘発されるこれらの貯留層内の過渡的な多相、多成分流体流、および物質バランス挙動、ならびに貯留層流体の圧力−体積−温度（ＰＶＴ）関係を表す非線形偏微分方程式のシステムである。

貯留層シミュレータは、体積をグリッドブロックとしても知られている隣接セルに細分することによって、地下貯留層および含まれる周辺の多孔質岩石層における多相多成分流体流および物質バランスをシミュレートする。したがって、シミュレーションモデルでは、貯留層が多数の個々のセルに編成される。セルまたはグリッドブロックは、基礎となる数学モデルが適用される基本的な有限体積である。セルの数は、シミュレーションに必要な分解能および問題の貯留層のサイズによって変化する。

当業界で巨大な貯留層として知られているタイプのような、数十億バレルの原始埋蔵量（ＯＯＩＰ）を有することができる大きな貯留層の場合、グリッドセルの数は、数億から１０億を超える数とすることができる。このセルの数は、流動力学、地層岩多孔率および透過率の不均質性、ならびに貯留層内の多くの他の地質学的および堆積上の複雑性を表すのに十分な分解能を有するために必要とされる。このサイズの貯留層のシミュレーションは、ギガセル貯留層シミュレーションと呼ぶことができる。図１は、１００において、貯留層シミュレーションのための地下貯留層の例示的な構造化貯留層グリッドモデルを示す。

炭化水素貯留層シミュレーションにおける課題は、費用効果的な方法で回収を最大化するために最新の技術の使用を必要とする。ＧｉｇａＰＯＷＥＲＳの商標で提供されるような貯留層シミュレーションが、文献に記載されている。例えば、Ｄｏｇｒｕ，Ａ．Ｈ．ｅｔ ａｌによる記事“Ａ Ｎｅｘｔ−Ｇｅｎｅｒａｔｉｏｎ Ｐａｒａｌｌｅｌ Ｓｉｍｕｌａｔｏｒ ｆｏｒ Ｇｉａｎｔ Ｒｅｓｅｒｖｏｉｒｓ”，ＳＰＥ １１９２７２，２００９ ＳＰＥ Ｒｅｓｅｒｖｏｉｒ Ｓｙｍｐｏｓｉｕｍ，Ｔｈｅ Ｗｏｏｄｌａｎｄｓ，Ｔｅｘａｓ，ＵＳＡ，Ｆｅｂ ２−４，２００９の議事録、およびＤｏｇｒｕ，Ａ．Ｈ．、Ｆｕｎｇ，Ｌ．Ｓ．、Ｍｉｄｄｙａ，Ｕ．、Ａｌ−Ｓｈａａｌａｎ，Ｔ．Ｍ．、Ｂｙｅｒ，Ｔ．、Ｈｏｙ，Ｈ．Ｈａｈｎ，Ｗ．Ａ．、Ａｌ−Ｚａｍｅｌ，Ｎ．、Ｐｉｔａ，Ｊ．、 Ｈｅｍａｎｔｈｋｕｍａｒ，Ｋ．、Ｍｅｚｇｈａｎｉ，Ｍ．、Ａｌ−Ｍａｎａ，Ａ．、 Ｔａｎ，Ｊ．、Ｄｒｅｉｍａｎ，Ｔ．、Ｆｕｇｌ，Ａ．、Ａｌ−Ｂａｉｚ，Ａ．による記事“Ｎｅｗ Ｆｒｏｎｔｉｅｒｓ ｉｎ Ｌａｒｇｅ Ｓｃａｌｅ Ｒｅｓｅｒｖｏｉｒ Ｓｉｍｕｌａｔｉｏｎ”，ＳＰＥ １４２２９７，２０１１ ＳＰＥ Ｒｅｓｅｒｖｏｉｒ Ｓｉｍｕｌａｔｉｏｎ Ｓｙｍｐｏｓｉｕｍ，Ｔｈｅ Ｗｏｏｄｌａｎｄｓ，Ｔｅｘａｓ，ＵＳＡ，Ｆｅｂ ２１−２３，２０１１の議事録、を参照されたい。ＧｉｇａＰＯＷＥＲＳ貯留層シミュレーションは、シナリオ当たり数百ギガバイト（ＧＢ）のフットプリントを利用しながら、後処理のための１０億セルバリアを超える微細スケールでのグリッドシミュレーションが可能である。

多数の炭化水素貯留層およびかなりの埋蔵量を有する会社がシミュレーションを実行する総数は年間数万を超え、これらのデータを提供するためには１ペタバイト以上の高性能ストレージが必要とされる。

多相多成分系の過渡解は、貯留層の初期条件から時間ステップのシーケンスにおける質量とエネルギー保存則の進展を含む。各時間ステップに対して、各有限体積に対する非線形離散方程式系は、一般化されたニュートン法として知られているものを用いて線形化される。

当業界では、これはニュートン反復法または非線形反復法と呼ばれる。各ニュートン反復法において、ヤコビ行列として知られる行列と残差として知られる右辺ベクトルが、系の主要変数の変化に対して解くために使用される、線形方程式系が構築される。

貯留層に関する他の要因は、ギガ規模の貯留層シミュレーションの複雑さをさらに増加させた。地層中の地下岩は、透過率または多孔率のいずれにおいても均質ではない。非常に不均質な貯留層では、貯留層圧力分布の計算において反復法が収束しない場合がある。これにより、シミュレーション処理にさらなる複雑さがもたらされる。

線形方程式系を解くための現在の業界の実践は、前処理つき反復法を介して行われている。前処理の例は例えば、Ｓａａｄの“Ｉｔｅｒａｔｉｖｅ Ｍｅｔｈｏｄｓ ｆｏｒ Ｓｐａｒｓｅ Ｌｉｎｅａｒ Ｓｙｓｔｅｍｓ”，Ｊａｎｕａｒｙ，２０００，Ｃｈ．１０，Ｐｐ．２６５−３１９に記載されている。計算上コストの高い前処理法は、計算コストを増大させる。さらに、前処理法は、正確な解への収束を保証するわけではない。

前処理行列は一般に、貯留層内の未知の圧力と飽和分布に対して良い推定値を提供することができる。そして、この初期推定値は、時間ステップに対する正しい圧力と飽和分布への反復によって改善される。前処理行列として使用される１つの方法は、その単純さと、反復あたりの比較計算コストが低いことから使用される、Ｇａｕｓｓ Ｓｅｉｄｅｌ（ガウス・ザイデル）法として知られている。ガウス・ザイデル前処理法は、点ガウス・ザイデルあるいは線形ガウス・ザイデル前処理として知られる形式をとることがある。次のシミュレータ時間ステップのための貯留層行列のグリッドセル圧力の直接解が決定され、残差が評価される。この処理は、グリッドについて計算された残差が指定された許容範囲未満になるまで続く。ガウス・ザイデル方法論の例は、Ｏｌｓｅｎ−Ｋｅｔｔｌｅの“Ｎｕｍｅｒｉｃａｌ Ｓｏｌｕｔｉｏｎ ｏｆ Ｐａｒｔｉａｌ Ｄｉｆｆｅｒｅｎｔｉａｌ Ｅｑｕａｔｉｏｎｓ”，Ｍａｒｃｈ，２０１１，Ｃｈａｐｔｅｒ ４．２，Ｐｐ．２８−３３，Ｔｈｅ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｑｕｅｅｎｓｌａｎｄに記載されている。しかしながら、ガウス・ザイデル前処理法は、指定された許容範囲内に残差があるグリッドセル圧力を得る上で、特に非常に不均質な地層岩石特性を有する貯留層に対して、しばしば非常に遅い収束を示すことが分かっている。

圧力方程式は貯留層シミュレーションにおいて重要な役割を果たす。貯留層シミュレーション中に決定された圧力条件は、貯留層のモデル中の貯留層グリッドセルに対して関心のある指定された時間ステップでの圧力条件を表す。貯留層シミュレーション処理の間、物理的または数値的処理の理由のいずれかにより、係数行列は、シミュレーションの間、非対角優位または不定になり得る。これは、反復法を発散させる。これまで知られているように、現在の実践は時間ステップサイズを縮小し、係数行列を再構成し、縮小したサイズの時間ステップで処理を繰り返すことである。こうした状況では、シミュレーション処理が、コンピュータの稼働時間および多数の未知数に対するコストの点で非常に高価になる可能性があり、これは任意の大型の貯留層、特に巨大な貯留層の場合に当てはまる。

簡単に説明すると、本発明は、貯留層からの予測生産の時間ステップのシーケンスについて、貯留層内の圧力分布を決定して、貯留層内のシミュレートされた流体流を決定する収束性を改善した、生産貯留層内の流体流の貯留層シミュレーションのコンピュータ実施方法を提供し、貯留層は貯留層セルの３次元グリッドに組織化され、貯留層シミュレーションは、貯留層内の油井からの実際の流体圧力および流体生産に基づいて、時間ステップのシーケンスについて実行され、コンピュータ実施方法は、貯留層セルの３次元グリッドから油井へのシミュレートされた流体流および油井のシミュレートされた油井生産速度を決定する。

コンピュータ実施方法は、時間ステップのシーケンスにおける、ある時間ステップについて、既知の貯留層属性の初期コンピュータ行列と、実際の貯留層流体生産測定値の初期コンピュータベクトルとを形成する。次に、既知の貯留層属性の初期コンピュータ行列が対角優位であるかどうかの判定が行われる。

対角優位である場合、圧力分布は、貯留層のグリッドセルの個々の中で圧力分布の形成された初期測定値に基づいて、貯留層シミュレーションの時間ステップについて、貯留層のグリッドセルの個々の中で決定される。次いで、貯留層のグリッドセル内の推定流体流量を、貯留層のグリッドセル内の決定された推定圧力分布に基づいて、貯留層シミュレーションの時間ステップについて決定する。次いで、コンピュータは貯留層のグリッドセル内の決定された流体流量が収束したかどうかを判定し、収束した場合、時間ステップについて貯留層内の流体流の貯留層シミュレーションが終了する。

既知の貯留層属性の初期コンピュータ行列が対角優位でない場合、近似解析プレコンディショナ（前処理法）が生成される。次に、生成された近似解析前処理法を既知の貯留層属性の初期コンピュータ行列に適用して、貯留層のグリッドセル内の圧力分布の測定値を形成する。Ｋｒｙｌｏｖベクトルは、貯留層のグリッドセル内の圧力分布の形成された測定値から生成され、貯留層圧力分布の係数ベクトルは圧力分布の形成された測定値から生成されたＫｒｙｌｏｖベクトルから決定される。

貯留層のグリッドセル内の推定流体流量を、貯留層圧力分布の決定された係数ベクトルに基づいて、貯留層シミュレーションの時間ステップに対して決定した。そして、貯留層圧力分布の決定された係数ベクトルに基づいて、貯留層のグリッドセル内の推定流体流量が収束したかどうかを決定する。収束した場合、貯留層内の流体流の貯留層シミュレーションは、その時間ステップについて終了する。

本発明はまた、貯留層からの予測生産の時間ステップのシーケンスについて、貯留層内の圧力分布を決定して、シミュレートされた流体流を決定する、収束性を改善した生産貯留層内の流体流の貯留層シミュレーションを実行する、新しい改善されたデータ処理システムを提供し、貯留層は、貯留層セルの３次元グリッドに組織化され、貯留層シミュレーションは、貯留層内の油井からの実際の流体圧力および流体生産に基づいて、時間ステップのシーケンスについて実行される。

データ処理システムは、貯留層セルの３次元グリッドから油井へのシミュレートされた流体流および油井のシミュレートされた油井生産速度を決定するプロセッサを含む。プロセッサは、時間ステップのシーケンスにおける、ある時間ステップについて、既知の貯留層属性の初期コンピュータ行列と、実際の貯留層流体生成測定値の初期コンピュータベクトルとを形成する。次に、既知の貯留層属性の初期コンピュータ行列が対角優位であるかどうかの判定がプロセッサにおいて行われる。

対角優位である場合、圧力分布は、貯留層のグリッドセルの個々の中の圧力分布の形成された初期測定値に基づいて、貯留層シミュレーションの時間ステップについて、プロセッサによって、貯留層のグリッドセルの個々の中で決定される。次いで、貯留層のグリッドセル内の推定流体流量は、貯留層のグリッドセル内の決定された推定圧力分布に基づいて、貯留層シミュレーションの時間ステップについて、プロセッサによって決定される。次いで、プロセッサは、貯留層のグリッドセル内の決定された流体流量が収束したかどうかを判定し、収束した場合、その時間ステップについてのプロセッサによる貯留層内の流体流の貯留層シミュレーションが終了する。

既知の貯留層属性の初期コンピュータ行列が非対角優位である場合、近似解析プレコンディショナ（前処理法）がプロセッサによって生成される。そして、生成された近似解析前処理法は、プロセッサによって、既知の貯留層属性の初期コンピュータ行列に適用されて、貯留層のグリッドセル内の圧力分布の測定値を形成する。Ｋｒｙｌｏｖベクトルは、貯留層のグリッドセル内の形成された圧力分布の測定値からプロセッサによって生成され、貯留層圧力分布の係数ベクトルは、形成された圧力分布の測定値から生成されたＫｒｙｌｏｖベクトルから決定される。

貯留層のグリッドセル内の推定流体流量は、貯留層圧力分布の決定された係数ベクトルに基づいて、貯留層シミュレーションの時間ステップについてプロセッサによって決定される。次いで、貯留層圧力分布の決定された係数ベクトルに基づいて、貯留層のグリッドセル内の推定流体流量が収束したかどうかを決定する。収束した場合、プロセッサによる貯留層内の流体流の貯留層シミュレーションは、その時間ステップについて終了する。

貯留層シミュレーションのための地下貯留層構造化グリッドのコンピュータ化モデルの等角図である。

簡略化された例示的な二次元貯留層モデルのセグメントの平面図である。

本発明による貯留層シミュレーションのための不均質な岩石層における圧力分布を決定するためのコンピュータネットワークの概略図である。

図３のコンピュータネットワークのアプリケーションサーバまたはコンピュータノードの概略図である。

不定貯留層係数行列の存在下での、本発明による貯留層シミュレーションのための不均質な岩石層における圧力分布を決定するための、ワークフローの機能ブロック図またはフローチャートである。 不定貯留層係数行列の存在下での、本発明による貯留層シミュレーションのための不均質な岩石層における圧力分布を決定するための、ワークフローの機能ブロック図またはフローチャートである。 不定貯留層係数行列の存在下での、本発明による貯留層シミュレーションのための不均質な岩石層における圧力分布を決定するための、ワークフローの機能ブロック図またはフローチャートである。

従来技術の方法と比較した、本発明による圧力分布の比較結果を得るための簡略化された貯留層モデルの概略図である。

図８の貯留層モデルに対する圧力分布圧力分布決定時の収束誤差の比較結果の表示図である。

従来技術の方法と比較した、本発明による圧力分布の比較結果を得るための簡略化された貯留層モデルの概略図である。

従来技術の方法と比較した、本発明による図１０の貯留層モデルに対する圧力分布圧力分布決定時の収束誤差の比較結果の表示図である。

従来技術の方法と比較した、本発明による圧力分布の比較結果を得るための簡略化された貯留層モデルの概略図である。

従来技術による処理結果に対する本発明による圧力分布決定における収束誤差の比較結果の表示図である。

数値貯留層シミュレーションでは、定式化に基づいて、圧力方程式が流れ方程式の本体に明示的に形成されるか、含まれるかのいずれかで形成される。流体速度、従って、油、水およびガスの流れは、圧力分布から計算される。したがって、圧力場が適切に計算されることが重要である。

貯留層シミュレーションでは、未知数を解くため、流れ方程式を線形化し、係数行列を形成する。係数行列は一般にヤコビ行列として知られている。係数行列またはヤコビ行列は、貯留層シミュレーションによって決定される変数に関する圧力などの未知数の一次導関数を含む。ヤコビ行列の一例を、以下に示す式（１）に示す。

地下の流体貯留層は岩石特性において非常に不均質である。これらの不均質性は、貯留層内の流体流と圧力分布を記述する偏微分方程式の係数に反映される。

多くの場合、シミュレータは、シミュレーション処理中に時間ステップで適切な圧力解を見つけるのが困難なため、収束問題に直面する。したがって、正確に計算された圧力解は、シミュレータが時間領域で進行し、妥当な計算コストで貯留層性能を決定するために重要である。

本発明は、貯留層シミュレータの圧力方程式が、グリッドセル内の圧力条件の決定に基づくパラメータに存在する不定係数行列と収束するために開発された。不定係数行列は、蒸気液体平衡に関連する物理的条件をシミュレートすることにより、または計算された不適切な導関数により数値的に作成された非物理的条件のいずれかにより、発生し得る。

以前の反復法は不定係数行列に対して発散する。本発明は、不定係数行列に遭遇した場合に収束する。本発明は、時間ステップサイズを縮小することなく、貯留層シミュレーションプロセスを継続することを可能にする。時間ステップサイズの縮小は、約数百万および数十億以上の多数の未知数が存在するコンピュータ処理時間および貯留層シミュレーション実行に関してコストがかかる。本発明は、数百万、数十億またはそれ以上の未知数を含む非常に大きな方程式系にとって非常に高価であり得る時間ステップ切断を回避する。本発明は、新しい反復法によって不定貯留層係数行列で圧力条件を決定し、時間ステップサイズを切断することなく、貯留層シミュレータがシミュレーションを継続することを可能にする。

圧力方程式は、貯留層シミュレーションにおける流体流をシミュレートする上で重要な役割を果たす。選択された数値解スキームによって、別個の圧力方程式が形成されるか、同じ未知数が他の貯留層パラメータ未知数と同時に解かれるかのいずれかである。本発明によれば、貯留層シミュレータは別個の圧力方程式を形成し、貯留層シミュレーションのための次の時間ステップ圧力は、式（１）に示されるように一般的な行列形式で表される一次方程式系を解くことによって解かれる。

貯留層に対する圧力分布決定のコンテキストにおける式（１）に対して、Ａは既知の係数行列であり、ｂは実際の貯留層生産測度の既知の右辺ベクトルであり、ｐは次の時間ステップに対する未知の圧力ベクトルである。係数行列Ａは、後述するように、流体流を表わす偏微分方程式の離散化から生じるｎ×ｎの疎行列である。式（１）において、ｂはｎ×１行ベクトルであり、これは通常、古い時間ステップ情報および油井生産／注入量を含み、ここで、ｎは未知の圧力ベクトルの総数である。

［用語］
参照および理解を容易にするために、本発明による様々なパラメータ、測定値、および測定単位を表す式で使用される用語のリストを以下に示す。

［微分形式および離散形式の圧力方程式］
一例として、貯留層内の流体の気泡点より上で動作する、流体特性が一定の正方形の平坦な単相石油貯留層Ｖ（図２）を考えてみる。線形化された、わずかに圧縮性の単相流体に対する貯留層圧力ｐを表わす偏微分方程式は、次のように与えられる。

ここで、ｐは、貯留層内で動作する流れ強度ｑを有する供給源またはシンク（油井）による未知の圧力分布である。パラメータｋは、貯留層モデル１００における次元座標のｘ、ｙ、ｚ方向における変化する透過率を表す。パラメータＣ_ｔは、流体と岩石の総圧縮率で、パラメータｔは時間である。この式の形式は、体積バランス関係を定義するものに相当する。

式（２）は、所与の外側境界条件および内側境界条件、ならびに貯留層内の圧力分布に関する初期条件に対して、ｐ（ｘ、ｙ、ｚ、ｔ）に対する一意解が存在することに基づいている。

簡略化された例（ｘおよびｙのみの２次元を考える）について、式（２）は以下のようになる。

貯留層モデルＶ（図２）のような簡略化された２次元貯留層モデルが、グリッドブロック（制御体積）に編成または分割されると仮定する。より一般化された例では、貯留層モデルは、ｘ次元では多数のグリッドブロックとしてのＮ_ｘ、ｙ次元では多数のグリッドブロックがＮ_ｙで構成される。２次元貯留層内のグリッドブロックの総数はｎで、ｎ＝Ｎ_ｘ＊Ｎ_ｙである。

図２のモデルＶでは、ｉおよびｊがそれぞれｘおよびｙ方向の付番シーケンス指標として割り当てられる。ｘ方向において、ｉ＝１，２．．．Ｎ_ｘであり、ｙ方向において、ｊ＝１，２．．．Ｎ_ｙである。以下の式（４）は、モデルＶについて、式（３）の微分方程式に対する離散化された系の有限差分形式を表す。式（３）の未知関数ｐを解くために、モデルＶのｘ−ｙ領域を図２に示す計算要素に分割すると、式（４）は次のように表される。

ここで、Ｔは、グリッドブロック（セル）（ｉ，ｊ）とその隣接ブロックとの間のリンクフローを表わす伝達率である。貯留層モデルにおける伝達率とは、流動流体の粘度を補正した地層グリッドセル要素の導電率を表す尺度である。Ｖ_ｐｉ，ｊは要素（ｉ，ｊ）の細孔容積であり、Ｃ_ｔ はその要素（ｉ，ｊ）の岩石と流体の総圧縮率である。

式（４）において、上付き文字ｎは古い時間ステップを示し、上付き文字（ｎ＋１）は新しい時間ステップを示し、Δtは新しい時間ｔ（ｎ＋１）と古い時間ｔ（ｎ）との間の時間ステップサイズであり、Ｖｐ（ｉ，ｊ）（はグリッドブロックの細孔容積である。

貯留層モデルＶを例示すると、座標（ｉ，ｊ）におけるセル１１０の西側伝達率は、 T_i-1/2,jで表され、ここで、（ｉ−１／２，ｊ）は座標（ｉ，ｊ）におけるグリッドブロック１１０の西側面１１２を意味する。同様に、T_i+1/2,jは同じセル１１０に対する東側面１１４の伝達率である。比較可能に、セル１１０の北側面１１６の伝達率はT_{i,j +1/2}で表され、セル１１０の南側面１１８の伝達率はT_{i,j -1/2・}で表される。

一般に、ｘ次元におけるセルブロックに対するグリッド内の面の伝達率Ｔは、以下のように定義される。

ここで、

はグリッドブロック（ｉ＋１，ｊ）と（ｉ，ｊ）との間の調和的に平均化された透過率であり、Ａはグリッドブロックの界面領域であり、ｘは、グリッドブロックのｘ座標である。式（５）において、Ｃ_１＝１．１２７Ｅ−３であり、これは単位換算係数である。換算係数Ｃ_１は、立方センチメートル／秒の計量単位のメートル系単位での流体流量を、実用油田単位、すなわちバレル／日に変換する。したがって、伝達率Ｔは、貯留層モデルにおけるｘまたはｙ距離がｆｔで表され、面積Ａがｆｔ^２で表され、透過率ｋが計測単位のダルシーではなくミリダルシーで表される場合、バレル／日／ｐｓｉの測定の体積、時刻、および圧力の単位で表される。

貯留層内の各グリッドブロックに対する時間ステップｎでの圧力分布は、前の時間ステップの結果として知られている。新しい時間ステップについての各グリッドブロックに対する未知の圧力は、セルについて式（４）を記述し、式（６）の行列形式に並べ替えることにより、線形方程式系を形成することにより得られる。

グリッドブロック（ｉ，ｊ）について式（４）を並べ替えると、以下の式（６）が得られる。

便宜上、式（６）の東側面指定子（ｉ＋１／２，ｊ）は東側面を示す略語「Ｅ」に置き換えられ、西側面指定子（ｉ−１／２，ｊ）は西側面を示す略語「Ｗ」に置き換えられ、北側面指定子（ｉ，ｊ＋１／２）は北側面を示す略語「Ｎ」に置き換えられ、南側面指定子（ｉ，ｊ−１／２）は南側面を示す略語「Ｓ」に置き換えられている。

そして、時間ステップ（ｎ＋１）における未知の圧力ベクトルに対する行列形式の式（６）は、以下の形式の式（７）として行列形式で表すことができる。

式（６）と同様に、便宜上、式（７）において、東側面指定子（ｉ＋１／２，ｊ）は東側面を示す略語「Ｅ」に置き換えられ、西側面指定子（ｉ−１／２，ｊ）は西側面を示す略語「Ｗ」に置き換えられ、北側面指定子（ｉ，ｊ−１／２）は北側面を示す略語「Ｎ」に置き換えられ、南側面日指定子（ｉ，ｊ＋１／２）は南側面を示す略語「Ｓ」に置き換えられている。油田単位では、式（７）はバレル／日／ｐｓｉ圧力変化、または立方センチメートル／秒／大気圧差の単位である。

式（７）の右辺において、ｑ_ｉは所与のグリッドブロックにおける生産または注入速度であり、項

は、グリッドブロックｋの細孔容積Ｖ_pと、総圧縮率Ｃ_Ｔと、時間ステップサイズΔt当たりのグリッドブロックｋの古い時間ステップ圧力ｐ_ｋとの積である。

式（７）の線形方程式系は、用語の項に従って定義されるパラメータで表される式（１）の行列形式で以下のように編成される。

ここで、式（８）において、［Ｔ］は式（１）の係数行列［Ａ］に対応する流体伝達率に関する係数行列であり、

は式（１）の右辺行列［ｂ］に対応するソース項、細孔容積、総圧縮率、および時間ステップサイズからなる式（８）の右辺ベクトルであり、ｐは、時間ステップ（ｎ＋１）における未知の圧力ベクトルである。

［貯留層モデル］
図１の１００で示されるような実際の例示的な地下貯留層の構造化貯留層グリッドモデルは、単一の未知の圧力ｐを各要素が有するグリッドブロック（計算要素）で構成されている。貯留層モデルは、数百万の要素に分割される。したがって、行列Ａは非常に大きく、すなわち、１００万×１００万の疎行列とすることができる。貯留層モデルに対する係数行列Ａのサイズが大きいために、ガウス消去法として知られるプロセスによるような解法の直接法は計算的に非常に費用がかかり、遅い。したがって、式（１）を解くために、反復解法が通常使用される。

［反復法への挑戦］
反復法は、行列が強く対角優位である場合に良好に機能する。これはまた、行列のすべての固有値が同じ符号（正または負のいずれか）を有することを意味する。行列が対角優位性を失う場合（固有値が符号を変える場合）、式（１）の係数行列Ａは不定となる。これまで知られている限り、反復法は、不定行列に収束しない。

［行列用語および反復法］
行列の固有値
ｎ行ｎ列を有する式（１）による正方行列［Ａ］の固有値は、行列［Ａ］の行列式に対する関係の次式の根として定義される。

その結果、次数ｎの多項式が次のように定義される。

行列［Ａ］の固有値は、式（１０）の多項式表現を満たす根である。係数行列は、貯留層内の圧力分布を解くために貯留層シミュレータで使用される、式（１）の同じ行列である。λによって表される行列［Ａ］の固有値は、式（５）、式（６）で定義される伝達率Ｔ、要素（ｉ，ｊ）の細孔容積Ｖ_ｐｉ,ｊ，および全圧縮率Ｃ_ｔなどの変数を含む。

［定値行列および不定行列］
正定値行列は、そのような行列の固有値のそれぞれがゼロより大きい行列である。負定値行列は、その行列Ａの固有値のそれぞれがゼロより小さい行列である。不定行列は、いくつかの固有値が正で、いくつかの固有値が負の行列である。ある行列の１つの固有値がその行列の固有値の残りと符号が異なる場合、その行列は不定である。

実際には、シミュレーション実行中に、任意のグリッドブロックの全圧縮率が負になると、係数行列Ａは不定になり、貯留層シミュレーションの反復法は発散する。上述したように、グリッドブロックに対する負の総圧縮率は、物理的な理由、またはある位相平衡計算における岩石および流体の微分値の組み合わせによるものである可能性がある。

［貯留層シミュレーションの反復法のための収束問題を有する不定行列の例］
圧力方程式の式（８）の係数行列［Ｔ］は通常、対角優位であり、従って正定値である。したがって、任意の反復法は、異なる収束率で真の解に収束する。収束は、係数行列Ａの要素の大きさによって影響を受ける。係数行列Ａが強く対角優位である場合、行列の行におけるすべての要素の絶対値の合計が、同じ行の対角項（数）の絶対値よりもはるかに小さい。

非対角要素と対角要素の和の差が小さい場合には、行列は弱対角優位と呼ばれる。強対角優位行列は、数回のコンピュータ処理反復内で真の解に収束する。一方、弱対角優位行列は収束に多くの反復を必要とする。しかし、これらのタイプの行列のそれぞれは、場合によっては１つ以上の行列行またはグリッドブロックにおいて対角優位性を失う。このような場合、係数行列［Ｔ］が不定となることがある。典型的な反復法は、不定係数行列が存在するときには収束しないであろう。

行列形式に再配置されたセルの各グリッドブロック線形方程式の未知の圧力を次のように表現する式（６）を再度呼び出す。

式（６）において、要素（ｉ，ｊ）の行列行は、以下のように表される。

すべての伝達率値Ｔ＞０であるので、絶対値の行列［Ｔ］の対角項は、

についてすべての対角外項の合計よりも大きい。または、

この条件が行列［Ｔ］のすべての行（すべてのグリッドブロック（ｉ，ｊ））について真である場合、係数行列［Ｔ］は正定値であり、固有値は同じ符号を有する。

式（４）で表されるように、圧力方程式行列において

が負の場合には、全圧縮率は負、つまりＣ_ｔ＜０である。係数行列は不定になり、少なくとも１つの固有値は他の固有値とは反対の符号を有する。従来の反復法は、これが起こると収束しない。

シミュレーション実行中に、発散、または収束の失敗が起こる場合、シミュレータは通常、時間ステップサイズを切断して、新しい係数行列Ｔを形成する。新しい係数行列Ｔを形成することは、通常、非常にコストがかかり、特に、未知の圧力分布の多数のグリッドブロックを有する貯留層モデルの場合、非常にコストがかかる。したがって、時間ステップサイズを切断しないことが重要である。本発明は、ある状況では不定であっても、同じ係数行列Ｔを使用する収束法を提供する。

式（８）の形の解を生成するために、システムは前処理されなければならない。例えば、先に述べたＯｌｓｅｎ−Ｋｅｔｔｌｅの“Ｎｕｍｅｒｉｃａｌ Ｓｏｌｕｔｉｏｎ ｏｆ Ｐａｒｔｉａｌ Ｄｉｆｆｅｒｅｎｔｉａｌ Ｅｑｕａｔｉｏｎｓ”に記載されているような、ガウス・ザイデル反復前処理法を使用することができる。

前処理も反復法であるため、式（１）を解くのに役立たない。これは、不定係数行列が存在すると、反復計算処理が発散するためである。したがって、本発明によれば、近似解析行列［Ｂ］前処理（近似解析解）は、１つ以上の不定係数行列が存在する場合にグリッドセル圧力分布を決定する。この処理は後述するように、図７による処理のステップ２２８の間に実行される。

本発明は貯留層モデルのグリッドセル間の圧力分布の決定中に、データ処理システムＳ（図３）の性能を改善する。収束は、条件が不定係数行列を存在させる場合、得難いことが分かる。これまで知られているように、圧力分布の決定のための以前の反復法は、不定係数行列が存在するときに発散する。

貯留層シミュレーション中の発散に対する時間ステップを切断する従来の処理行為は回避される。本発明は、時間ステップサイズを切り取ることなく、貯留層シミュレーションプロセスを継続することを可能にする。時間ステップの切断は、数百万、数十億またはそれ以上の数の未知のパラメータ値を含むモデルを有する非常に大きな貯留層に対して、コンピュータの使用および時間の点で非常にコストがかかることが分かっている。

本発明によれば、式（１）を行列形式で解くためのコンピュータ処理による圧力分布の決定は、Ｋｒｙｌｏｖベクトルとして知られるものの線形結合として表すことができ、次のように表すことができる。

式（１３）において、合計として表されるベクトルの数の上限、すなわちベクトルの数ｍは、理想的には決定されるべき貯留層グリッド要素内の圧力分布に関する未知の値の総数に等しい。しかしながら、後のセクションで説明するように、未知の値の総数は、Ｋｒｙｌｏｖベクトルの数よりもはるかに少ないことが分かっている。

以降のセクションでは、Ｋｒｙｌｏｖベクトルとして線形独立ベクトルＶ_ｉを構成し、本発明によるコンピュータ処理によって圧力分布を決定するための係数Ｃ_ｉを選択する方法を説明する。

近似解析解は、直接解析解法を必要とせずに、全く同じ方法で不定行列システムを解くための前処理行列として使用される。これは、以下に説明するように、いくつかの例によって実証可能である。

［データ処理システムＳ］
本発明による処理およびワークフローは、様々な高性能コンピュータ（またはＨＰＣ）クラスタコンピュータ／プロセッサノードコンピュータシステム上での展開に適している。これらは通常、マルチコアアーキテクチャの複数のＣＰＵを含む複数のコンピュートノードを備えたラックマウント型ハードウェアである。ノードは通常、従来の低遅延高帯域幅ネットワーク、スイッチ、およびルータと相互接続される。

本発明は、ニューヨーク州アーモンクのＩｎｔｅｒｎａｔｉｏｎａｌ Ｂｕｓｉｎｅｓｓ Ｍａｃｈｉｎｅｓ（ＩＢＭ）または他のソースから入手可能なものなど、任意の従来のタイプの適切な処理容量のメインフレームコンピュータ上で実行することができることを理解されたい。

このシミュレーションシステムと共に使用するための典型的なＨＰＣ環境は、今日のマルチノード、マルチＣＰＵ、マルチコア計算クラスタである。このようなクラスタの一例は、図３のデータ処理システムＳのＣに図示されている。クラスタＣは、ルータサーバまたはサーバ１５４によって矢印１５２で示すように並列にデータが提供される複数のコンピュータノード１５０（図３および図４）から構成されている。必要に応じて、この目的のためにいくつかのこのようなルータサーバを使用することができる。上述のタイプの元のシミュレーションまたは入力データは、適切な数のデータ記憶／ファイルサーバ１５６に記憶される。データ記憶／ファイルサーバ１５６は、データ処理システムＳによって使用されるための操作命令、制御情報およびデータベース記録を記憶する。

ルータサーバ１５４は、メモリに記憶されたコンピュータコード１５５の制御下で動作し、記憶サーバ１５６からの入力シミュレーションデータ、ならびに矢印１５８で示すシミュレーション処理結果を、クラスタＣのコンピュータノード１５０との間で並列に転送する。

本発明によるプログラムコード１５５は、サーバまたはサーバ１５４に、貯留層シミュレーション結果のインデックス、順序付け、処理、および転送を行わせる、非一時的コンピュータ動作可能命令の形態である。典型的には、データ処理システムＳは、ネットワーク１５９によってシステムに接続された適切な従来型のワークステーション１５７の一式を含む。ワークステーション１５７は、貯留層技術者が貯留層シミュレーション結果を観察するように、出力データおよび表示を提供する。貯留層シミュレーション結果に基づいて、技術者は油井および貯留層の生産および性能を監視および調整するために、必要に応じて評価および措置を講じることができる。

クラスタＣのコンピュータノード１５０は、コンピュータノード１５０のメモリ１６４に記憶されたコンピュータコードまたはプログラム製品１６２の命令の下で並列に動作する、図４に示されたタイプの複数のプロセッサまたはコア１６０を含む。本発明によるプログラムコード１６２は、データプロセッサ１６０に貯留層シミュレーションモジュール１６６内で貯留層シミュレーションを実行させる非一時的コンピュータ動作可能命令の形態である。

プログラムコード１５５および１６２はデータ処理システムＳの機能を制御し、その動作を指示する特定の順序付けられた動作のセットを提供するマイクロコード、プログラム、ルーチン、またはシンボリック・コンピュータ動作可能言語の形態であり得ることに留意されたい。プログラムコード１５５および１６２の命令は、サーバ１５４またはプロセッサノード１５０のメモリ、またはコンピュータディスケット、磁気テープ、従来のハードディスクドライブ、電子読取り専用メモリ、光記憶装置、またはそこに記憶された非一時的なコンピュータ使用可能媒体を有する他の適切なデータ記憶装置に記憶することができる。プログラムコード１６２はまた、図示のように、コンピュータ可読媒体としてサーバ１５６のようなデータ記憶装置上に収容することができる。

［データ処理システムＳによる処理］
本発明は、貯留層１００（図１）のような、貯留層の不均質な岩石層の３次元にわたる圧力分布および圧力値の決定法を改善する。本発明はさらに、係数行列が不定になり、従来の前処理が収束しなくなる場合に、地下の貯留層１００に対する貯留層シミュレーション中に、貯留層モデルを通して圧力決定のための迅速な初期圧力測定を可能にする。

本発明は、貯留層シミュレーションの時間ステップサイズにおける望ましくない切断を必要とせずに、それらの状況における貯留層シミュレーション中の収束問題を回避する。前述したように、時間ステップサイズの切断は、コンピュータ処理操作および複雑性および貯留層シミュレーション処理時間長を増加させる。決定される初期圧力分布は、貯留層シミュレーションのための貯留層内の正しい圧力プロファイルを表す。

初期圧力の推定値は、簡略化された貯留層モデルにおける解析解に基づく実際の解に近いことが例によって証明される。このようにして決定された初期圧力推定値は、貯留層シミュレータが指定された精度限界内で所望の解に収束することを可能にする。また、本発明による初期圧力分布の決定は、点ガウス・ザイデルまたは線形ガウス・ザイデルのいずれかまたは両方の前処理のような他の前処理法よりも少ないコンピュータ処理反復で収束する。

図５は、貯留層の不均質な岩石層の３次元にわたる圧力分布の改善された決定法に基づく貯留層シミュレーションのための、本発明によるシーケンスの好ましい実施形態を概説するワークフローＦの概略図である。

ステップ２００の間、貯留層および生産データを読み取ることによってシミュレーションが開始する。ステップ２００のデータの間に読み込まれた貯留層データは、貯留層の幾何学的情報、すなわち、そのサイズ、ｘ、ｙ、およびｚ方向の範囲（長さ）、透過率分布などの貯留層特性、多孔率分布、層の厚さ、相対透過率データ、毛管圧データ、流体密度表などの流体特性データ、地層体積係数表、粘度表、油井の位置、および貯留層内の油井穿孔の位置を含む。生産データには、前のステップで定義された油井について測定または指定された石油、水、およびガス生産速度が含まれる。生産データには、各油井の最低坑底圧力も含まれる。貯留層シミュレータもステップ２００中に初期化され、シミュレーション日および時間ステップがゼロに設定される。

ステップ２０２の間、貯留層データおよび貯留層特性は式（１）に従って、貯留層モデル１００の貯留層グリッドセルブロックのためのＡ行列として表される形式に、データ処理システムＳ内で編成される。ステップ２０２の間、貯留層内の油井についての生産データは、式（１）に従って、貯留層モデル１００の油井のｂベクトルについても編成される。

ステップ２０４の間、係数行列Ａは、行列Ａ内の係数のうちの任意の１つ以上の係数の値が負であるかどうかを判定するために評価される。負でない場合、処理はステップ２０６に進む。

ステップ２０６の間、貯留層モデル１００のグリッドセルに対する最終的な圧力分布が決定される。これは、例えば、本出願人の、２０１８年４月２６日に出願された“Ｄｅｔｅｒｍｉｎｉｎｇ Ｐｒｅｓｓｕｒｅ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｉｎ Ｈｅｔｅｒｏｇｅｎｅｏｕｓ Ｒｏｃｋ Ｆｏｒｍａｔｉｏｎｓ ｆｏｒ Ｒｅｓｅｒｖｏｉｒ Ｓｉｍｕｌａｔｉｏｎ”という先願の米国特許出願第１５／９６３，２５１号の技術に従って行うことができ、これは参照により本明細書に組み込まれる。

また、貯留層内の圧力分布を決定するために、ステップ２０６の間に他の多くの反復法を使用することができることを理解されたい。貯留層内の圧力分布を決定するために正定値係数行列が存在する場合に使用できるこのような反復法の例は、Ｋｌａｕｓ Ｓｔuｂｅｎの“Ｓｏｌｖｉｎｇ Ｒｅｓｅｒｖｏｉｒ Ｓｉｍｕｌａｔｉｏｎ Ｅｑｕａｔｉｏｎｓ” Ｎｉｎｔｈ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｆｏｒｕｍ ｏｎ Ｒｅｓｅｒｖｏｉｒ Ｓｉｍｕｌａｔｉｏｎ，Ａｂｕ Ｄｈａｂｉ，Ｕｎｉｔｅｄ Ａｒａｂ Ｅｍｉｒａｔｅｓ，ｐｇ．１−５３；Ｄｅｃｅｍｂｅｒ ２００７に記載されている。

ステップ２０６の間の反復ソルバーによるデータ処理システムＳにおける処理は、収束が達成されるまで継続する。次いで、そのように決定された圧力は、貯留層のグリッドセル内の流体流を決定するために、ステップ２０６の間に貯留層シミュレーションの一部として提供される。そして、ステップ２０６の間に決定された複数次元における貯留層セル内の圧力分布は、貯留層シミュレーションステップ２０８の実行のために利用可能である。

ステップ２０８では、貯留層シミュレータによって貯留層シミュレーションが実行され、ステップ２０６から得られた貯留層モデル１００のセルについての決定された圧力分布と、貯留層生産データとの間の連続する時間ステップを使用して、式（７）に示されるように、複数次元で編成された貯留層セル内の多相流体に関連する、貯留層グリッドセル流体組成、圧力、流量、相対的存在および関連する地層パラメータを決定し、次の流体流動貯留層シミュレーション時間ステップ（ｎ＋１）についての圧力分布の推定値を計算する。

貯留層シミュレーションステップ２０８の間、ステップ２０６の間に決定された圧力推定値は、データ処理システムＳの貯留層シミュレーションモジュール１６６内の流体流貯留層シミュレーションにおいて使用される。ステップ２０８の間の貯留層シミュレーションは、ステップ２１０で示されるように、反復時間ステップのシーケンスで圧力として使用されるべき数値解法によって圧力を決定するために、複数の反復で実行される。貯留層シミュレーションは、先に述べたようないくつかの既知のタイプのうちの１つであり得る。

そして、ステップ２０８の間に得られた貯留層シミュレーション結果は、ステップ２１２の間に、１つ以上の油井、したがって貯留層の生産または性能のいずれか、あるいは生産および性能の両方を監視し、調整するための基準となる。ステップ２１２の貯留層エンジニアリングアクティビティは例えば、貯留層内の油井の生産からの流体の生産を調整する、貯留層の注入油井への二次回収流体の注入を調整または追加する、および貯留層内の追加の油井を穿孔する、という貯留層生産管理オペレーションのうちの１つ以上を含む。

ステップ２０４の間に、係数行列Ａ内の１つ以上の係数が負である場合、ステップ２２０（図５および図６）において、本発明に従って処理が開始される。ステップ２０２から生じる貯留層モデル１００のグリッドブロックセルに対する初期グリッドブロック圧力はステップ２２２の間に概略的に示されるように、処理のために利用可能にされる。

［システムのプレコンディショニング］
データ処理システムＳの式（１）、すなわち、［Ａ］［ｐ］＝［ｂ］を解くためには、実質的な反復回数内で収束を達成するために前処理を行わなければならない。前処理が用いられない場合、解は、収束のために過剰な反復回数を要する。これは、コンピュータ処理時間およびコンピュータ処理コストの観点から、高価になる可能性がある。計算的に安価な前処理の１つは、ガウス・ザイデル前処理法として知られているものである。しかし、不定係数行列で議論されてきたようなガウス・ザイデル前処理法が存在する。

図６および図７は、ステップ２２０による処理を概略的に示しており、各時間ステップに対して貯留層圧力分布を得て、貯留層モデルの各グリッドブロック（ｉ，ｊ）に対する圧力ｐを求める。各グリッドブロックに対する圧力の決定は式（１）に従って［Ａ］［ｐ］＝［ｂ］で表される線形方程式の行列システムに対してデータ処理システムＳにおいて実行され、ここで、Ａは既知の係数行列であり、ｂは行列システムの右辺に対する既知の大きさのベクトルであり、ｐは未知のグリッドブロック圧力のベクトルである。

ステップ２２０の間にこのようなシステムを解くために、データ処理システムＳ内に前処理行列が形成される。これは、行列式（１）に行列［Ｂ］の逆数を掛けることによって行われ、行列［Ｂ］は行列［Ａ］の近似である前処理行列を表す。本発明によれば、前処理行列［Ｂ］は、逆行列を有するように選択され、行列［Ａ］を乗算したときに計算するためのコンピュータ処理時間の点で安価である。本発明による前処理システムは、後述するようにＫｒｙｌｏｖベクトルを一緒に加算することによって実行される各グリッドブロックの圧力に対する解を可能にする。

ステップ２２２の間に、貯留層モデル行列のグリッドセルに対する初期近似圧力ｐ［ｒ（ｘ，ｙ）］が決定されると、ステップ２２４で概略的に示されるように、式（１）の行列式関係の圧力が前処理される。

ステップ２２４による処理の間に、最初のステップ２２８として、近似分析前処理法［Ｂ］が得られ、データ処理システムＳが前処理法を決定する。本発明により形成された前処理法は、近似解析解前処理法と呼ばれる。

行列［Ｂ］はＡに対する近似であり、式（１）に対する「前処理行列」と広く呼ばれている。行列［Ｂ］はその逆数が式（１）の右辺によって乗算されるときに、貯留層内の未知の圧力分布に対して近似が提供されるような形式のものである。本発明では、近似解析解から行列［Ｂ］の要素を求める。行列［Ｂ］はそれが容易に可逆であるように、係数行列［Ａ］に対する合理的に良い近似を表す。式（１）を左から逆数［Ｂ^−１］で乗算すると、式（１４）が得られる。

これは、式（１５）のように略すことができる。

ここで、

全圧縮性Ｃ_Ｔの尺度は、次のように示される細孔（岩石）と流体圧縮性の線形結合である。

ここで、変数名と単位は用語の項で説明されている。単一のソースP_max
による最大圧力値（または最小値）は、ステップ２２４の間に決定される。

簡略化のために、古い時間ステップｎにおける初期貯留層圧力はゼロであると仮定する。式（２０）および（２１）に従って圧力の値を等式化すると、未知の最大圧力P_maxは、式（２２）に従って表すことができる。

式（１６）と式（２２）は、貯留層内の新しい時間ステップにおけるおおよその圧力分布を表わしている。本発明の目的のために、式（２）による時間ステップは、単一の生産井についての式（１）の解に対する近似と考えることができる。この決定は、貯留層の透過率行列が不定であるときに収束しない反復前処理法の使用の要件を回避する。

［複数の油井］
貯留層が複数の生産井または注入井を含む場合、新しい時間ステップでの圧力分布は、再び古い時間ステップ圧力がゼロであると仮定する簡略化のために、重ね合わせによって近づけることができる。貯留層内にｎ_ｗ個の油井がｑ_ｉの割合であると仮定し、ｉ＝１，…，ｎ_ｗ，であるとき、単独で作用する各油井ｉｗの最大圧力値ｐ_{ｍaｘ（ｉｗ）}は式（５）と式（６）から求められる。油井ｉｗから測定された任意の半径ｒにおける近似圧力値ｐ［ｒ（ｘ，ｙ）］は、単一の生産井の各々について決定された個々の圧力解を複合解に重ね合わせることによって決定される。

ここで、ｒ（x（ｉｗ）,ｙ（ｉｗ））は、油井（ｉｗ）の座標である。

ステップ２２８の結果として近似解析前処理法が得られた後、ステップ２３０が実行されて、Ｋｒｙｌｏｖベクトルが生成される。ステップ２３０の間にＫｒｙｌｏｖベクトルに対する係数が決定され、その結果、ステップ２３４の間に決定された係数が一緒に加算されて、ステップ２２６の貯留層モデルの各グリッドブロック（ｉ、ｊ）に対する圧力を得ることができる。

ステップ２３０（図７）の間のデータ処理システムＳは、グリッド係数行列からＫｒｙｌｏｖベクトルとして知られるものを形成する。ステップ２３０の間のベクトルの生成は、式（９）から生成されたベクトルのＫｒｙｌｏｖ部分空間で行われる。ステップ２３０では、Ｋは式（８）、したがって式（１）から生成された各Ｋｒｙｌｏｖベクトルの集合を表す。Ｋｒｙｌｏｖベクトル形式のＫ集合は、次のように表される。

式（２３）を用いてv₂について解くと、式（２６）の右辺Av₁は、式（２３）の生産速度ｑと同様に扱うことができる。
第ｉのベクトルについては、

Ｋｒｙｌｏｖベクトルは、行列乗算を実行する必要性を回避する。データ処理システムＳにおいて、本発明に従って得られるＫｒｙｌｏｖベクトルは、貯留層グリッドセル圧力測定のための利点を提供する。係数行列［Ａ］は行列対角線の中の特定の部分だけがゼロではないが、大部分の行列がゼロの要素を持つという点で疎である。行列のこれらの部分は値としてゼロを有するので、コンピュータ記憶装置を行列のゼロ要素に割り当てる必要はない。

しかしながら、式（１）の解に対して、行列［Ａ］は反転されなければならない。しかし、生産有意の貯留層に対する行列［Ａ］のいかなる逆行列も完全行列になり、そこではすべての行列要素はゼロではない。このような行列の反転プロセスは非常に困難であり、計算集約的であり、高価である。計算および記憶は極めて高く、すなわち、１００万×１００万グリッドセル完全行列の形態の貯留層行列は、ほとんどのコンピュータでは不可能なコンピュータメモリ内の１兆の記憶位置を必要とする。Ｋｒｙｌｏｖベクトルは疎行列ベクトル乗算のみの使用を可能にし、その結果の各々はＮ個の未知数のベクトルをもたらし、したがって、はるかに少ないコンピュータメモリ記憶装置しか必要としない。Ｋｒｙｌｏｖベクトルは、計算機コストおよび記憶要求はほとんどなく、行列式（１）に対する解の式を提供する。

ベクトルを互いに垂直にするために、直交化および正規直交化が形成され、正規化は、ベクトル長を１にするだけである（正規化）。直交ベクトルの総和で表現される未知ベクトル（解ベクトル）を表現するために、ベクトルが互いに垂直（直交）になるように、直交化および正規直交化が行われる。

ステップ２３２中のＫｒｙｌｏｖベクトルＫの変換は、グラム・シュミット（Ｇｒａｍ Ｓｃｈｍｉｄｔ）またはアーノルディ（Ａｒｎｏｌｄｉ）ベクトル正規直交化法として知られているものを使用するなど、いくつかの従来の方法で行うことができる。このような変換例は、Ｊｉｍ ＬａｍｂｅｒｓのＭａｔ ４１５／５１５： Ｇｒａｍ−Ｓｃｈｍｉｄｔ Ｏｒｔｈｏｇｏｎａｌｉｚａｔｉｏｎ Ｌｅｃｔｕｒｅ ３ Ｎｏｔｅｓ（２０１４）ｈｔｔｐ：／／ｗｗｗ．ｍａｔｈ．ｕｓｍ．ｅｄｕ／ｌａｍｂｅｒｓ／ｍａｔ４１５／ｌｅｃｔｕｒｅ３．ｐｄｆに記載されている。

［未知係数Ｃ_ｉの計算］
ステップ２３４は、貯留層セルモデルのグリッドブロック間の圧力分布を決定するための未知係数Ｃ_ｉの決定を含む。ステップ２３４の目的のために、式（１）による貯留層グリッドブロック内の圧力分布のデジタル化された形式に対する解は、以下の形式で表される。

なお、式（２８）の直交ベクトルｕ_ｉは、実際には式（１）の近似解である。ベクトルｕの値が完全な解である場合には、式（２８）において、係数の数は１だけであり、ｍ＝１であり、解を表現するのに１つのベクトルだけが必要である。したがって、ｕベクトルの値に対する解の合理的に正確な推定値を得ることが重要である。次に、式（９）の前処理システムの剰余は、以下のように定義される。

式（２８）を式（１５）に代入し、以下の式（３０）において二次関数ｊを定義する。

式（３１）の線形システムは通常、より小さいサイズ（Ｋｒｙｌｏｖ系から生成される線形独立ベクトルの数は通常、ｎよりはるかに少ない）を持つべきであり、それゆえ、直接法−ガウス消去法による処理によって解くことができる。その例がＷｉｋｉｐｅｄｉａのガウス消去法１〜８ページ、ｈｔｔｐｓ：／／ｅｎ．ｗｉｋｉｐｅｄｉａ．ｏｒｇ／ｗｉｋｉ／Ｇａｕｓｓｉａｎ＿ｅｌｉｍｉｎａｔｉｏｎで説明されている。

ステップ２３４の間に決定された圧力推定値は、ステップ２２６の間に、貯留層シミュレーションのために、ステップ２０８の間に油井および貯留層の生産および性能を監視および調整するために、貯留層セル内の多相流体に関連する貯留層グリッドセル流体組成、圧力、流量、相対的存在、および関連する地層パラメータを決定するために説明された方法で使用される。

［結果］
本発明の第１の例は、係数行列Ａの固有値を容易に決定可能である不定行列の場合を例示するために、小さな問題上で実証される。

［例１−２次元不均質貯留層、５×５グリッド］
図８は、従来技術の方法の結果と比較した、本発明による圧力分布の比較結果を得るための、簡略化された貯留層モデル３００の概略図である。例示的な貯留層３００は石油貯留層であり、僅かに圧縮性流体の気泡点圧力以上で動作する。３，０００フィート×３，０００フィートの正方形の面積を有し、流れの境界条件がなく、複数の油井を注入または生産することができる。この例では、中央に位置する３０２で示されるように、単一の注入油井が存在すると仮定される。したがって、面積グリッドサイズはΔｘ＝Δｙ＝６００フィートであり、Δｚ＝３０フィートの厚さであり、セルの数Ｎ_ｘ＝５、Ｎ_ｙ＝５であった。

［貯留層および流体データ］
例示的なシミュレーションでは、初期貯留層圧力は３，０００ｐｓｉであり、単相流体粘度＝１ｃＰである。全圧縮率は、注入油井が配置される中心の位置３０２を除いて、貯留層全体について３Ｅ−６ １／ｐｓｉと仮定される。中央セルにおける全圧縮率C_tは、−３００Ｅ−６ １／ｐｓｉと仮定した。以下に記載されるように、貯留層モデル３００のＸおよびＹ次元にランダムな透過率および空隙率分布が存在した。

このように係数行列Ａのすべての固有値は、正の符号を有する第４の固有値（１.０１）を除いて、負の符号を有する。したがって、係数行列Ａは不定である。

［シミュレーション］
比較目的のための３つの前処理法、点ガウス・ザイデル前処理法、線ガウス・ザイデル前処理法、および本発明による方法論を用いて、貯留層シミュレータを１つの時間ステップについて実行して、この実施例の圧力分布を３０日の時点で決定した。

ガウス消去法によって決定される図８の簡略化されたモデルについて決定される圧力分布、ならびに本発明による選択された反復回数の後の圧力分布は、以下のように表１、表２および表３に記載される。

図９は、貯留層モデル３００についての従来技術による処理結果に対する、本発明による圧力分布決定における収束誤差の比較結果の表示である。点ガウス・ザイデル（ＰＧＳ）の収束挙動は３０４で示され、線ガウス・ザイデル（ＬＧＳ）法の収束挙動は３０６で示される。本発明による収束挙動は、図９の３０８に示されている。

図９に見られるように、縦軸にバレル／日（ｂ／ｄ）の最大残差が示され、横軸は反復回数を示す。示されるように、本発明が収束する間に、点ガウス・ザイデル法と線ガウス・サイデル法の両方が発散する。

例えば、１０回の反復の終わりに、最大残差（誤差）は点ガウス・ザイデル法では１，５７１，８４６ｂ／ｄであり、線ガウス・ザイデル法では６，６５７，１０２，８２９ｂ／ｄであり、本発明ではわずか２．５ｂ／ｄである。本発明は正しい解に対する収束解を提供するが、対照的に、点ガウス・ザイデル法および線ガウス・ザイデル法は提供しない。さらに、本発明は、完全行列ではなく、疎行列を使用する。したがって、本発明は、何億もの未知数を有する行列を解くために収束する有効な代替案を提供する。これは、貯留層行列処理解のためのコンピュータ処理および動作時間の節約であり、収束することなくむしろ発散する反復を繰り返し実行し続ける必要がない。また、さらに複雑なコンピュータ処理を導入する時間ステップ切断の従来のアプローチは必要なく、実行される必要がない。これは、本発明によるコンピュータ操作のさらなる改善である。

［例２−１００×１００グリッド］
図１０は、実施例２についての簡略化された例示的な貯留層モデル４００の概略図である。油井は、モデル４００の中央（５０，５０）付近の４０２に位置する。位置４０２における油井は実施例２の目的のために１，０００ｂ／ｄを注入している。モデル４００内の貯留層の長さは、各ｘ、ｙ方向において６０，０００フィートであり、実施例１と同じグリッドサイズをもたらす。セル（５０，５０）および（５１，５０）の総圧縮率Ｃ_Ｔは−３００Ｅ−６，１／ｐｓｉであり、モデル４００の残りのセルの総圧縮率は３Ｅ−６ １／ｐｓｉ総圧縮率である。実施例２の貯留層モデル４００は、実施例１において同じ垂直厚さおよび統計的特性を有する不均質貯留層として構成される。

図１１は、従来技術の点ガウス・ザイデルまたは線ガウス・ザイデル前処理法に従って得られた処理結果に対する、本発明に従って決定された貯留層モデル４００に対する圧力分布決定における収束誤差の比較結果の表示である。

図１１の４０４と４０６に示すように、点ガウス・ザイデル法および線ガウス・ザイデル法の両方が発散した。４０８で示されるように、本発明は、残留誤差を、６５ｂ／ｄの誤差に対応するグリッドブロック細孔容積の１．３パーセントに低減することができた。本発明の反復法によるさらなる反復は、定常的であるが残差のわずかな減少を示す。この例はさらに、本発明が正しい解を見つけることができる一方で、他の反復法が収束せず、貯留層シミュレーション実行を完了する試みにおいて時間ステップ切断が必要になることを実証する。

［例３−１００万のセルの不均質な問題（１，０００×１，０００）］
図１２は実施例３の簡略化された例示的な貯留層モデル５００の概略図であり、油井５０２はモデル５００の中心に示されているように配置され、４つの追加の油井５０４、５０６、５０８および５１０はモデル５００の各方形枠の中心に配置されている。これらの５つの油井の各々はこの実施例の目的のために１，０００ｂ／ｄで注入している。負の圧縮性グリッドブロックは、現在、油井５０２の西、東、北、南の１ブロックのグリッドセルと、モデル５００の中心の油井ブロックに対して確立されている。残りのパラメータ、垂直厚さおよび統計的特性は、実施例２および３と同じである。

図１３は、モデル５００に対する圧力分布決定における収束誤差の比較結果の表示である。５３０および５３２に見られるように、点ガウス・ザイデル法および線ガウス・ザイデル法の両方が発散する。本発明は、５３４で示されるような収束挙動を示す。初期残差（最大ノルム）１,０００ｂ／ｄは、グリッドブロック細孔容積の１２．８％に相当する１０回の反復で７０７．８ｂ／dに低減された。さらなる反復は残差を減少させることが分かり、本発明による処理は収束する。

この場合も、図１３の５３０および５３２に示すように、点ガウス・ザイデル法と線ガウス・ザイデル法の両方が発散した。この追加の例は、本発明が正しい解を見つけることができる一方で、他の反復法が収束せず、貯留層シミュレーション実行を完了する試みにおいて時間ステップ切断を必要とすることをさらに実証する。

したがって、図９、図１１および図１３のそれぞれにおいて、従来の点ガウス・ザイデルおよび線ガウス・ザイデル反復前処理法の両方が発散したことが分かる。具体的には、圧力計算における誤差が反復の数が増加することにつれて増加した。従って、これらの従来の前処理法は正しい圧力分布をもたらすことができない。一方、図９、図１１および図１３に示すように、実施例１、実施例２および実施例３のそれぞれにおける本発明は、連続反復の結果、圧力測定の誤差が低減されていることが分かる。したがって、本発明は正確な圧力解を提供し、反復の数が増加することにつれて収束する。

したがって、本発明は、貯留層シミュレーション動作中のコンピュータの機能を改善することが分かる。本発明は、不定係数行列が存在すると判定された場合に、貯留層シミュレーションのための貯留層モデル内の圧力分布の決定中に、コンピュータが収束する能力を提供する。本発明は、時間ステップ切断のような状況において、従来技術の修復作用を回避する。

本発明は、貯留層モデリングおよびシミュレーションの分野における平均的な知識を有する人が本明細書の本発明に記載された結果を再現し、得ることができるように、十分に説明されてきた。それにもかかわらず、本明細書の本発明の主題である技術分野の当業者は、本明細書の要求に記載されていない修正を実行して、これらの修正を決定された構造および方法論に適用することができ、またはその使用および実施において、以下の特許請求の範囲に記載された事項を必要とし、そのような構造およびプロセスは、本発明の範囲内に包含されるものとする。

添付の特許請求の範囲に記載される本発明の精神または範囲から逸脱することなく、上記に詳細に記載された本発明の改良および修正がなされ得ることに留意され、理解されたい。

Claims (17)

1. 貯留層からの予測生産の時間ステップのシーケンスについて、前記貯留層内の圧力分布を決定して前記貯留層内のシミュレートされた流体流を決定する、収束性を改善した生産貯留層内の流体流の貯留層シミュレーションのコンピュータ実施方法であって、前記貯留層は、貯留層セルの３次元グリッドに組織化され、前記貯留層シミュレーションは、貯留層内の油井からの実際の流体圧力および流体生産に基づいて、時間ステップのシーケンスについて実行され、前記コンピュータ実施方法は、貯留層セルの３次元グリッドから油井へのシミュレートされた流体流および前記油井のシミュレートされた油井生産速度を決定し、前記コンピュータ実施方法が、
   （ａ）前記時間ステップのシーケンスにおける、ある時間ステップについて、既知の貯留層属性の初期コンピュータ行列と実際の貯留層流体生産測定値の初期コンピュータベクトルとを形成するステップと、
   （ｂ）既知の貯留層属性の前記初期コンピュータ行列が対角優位であるかどうかを判定するステップと、
   （ｃ）対角優位である場合、前記貯留層のグリッドセルの個々の中の圧力分布の形成された初期測定値に基づいて、貯留層シミュレーションの時間ステップについての、貯留層のグリッドセルの個々の中の圧力分布を決定するステップと、
   （ｄ）前記貯留層のグリッドセル内の決定された推定圧力分布に基づいて、前記貯留層シミュレーションの時間ステップについての、前記貯留層のグリッドセル内の推定流体流量を決定するステップと、
   （ｅ）前記貯留層のグリッドセル内の決定された流体流量が収束したかどうかを決定するステップと、
   （ｆ）収束した場合、前記時間ステップについて、前記貯留層内の流体流の貯留層シミュレーションを終了するステップと、
   （ｇ）既知の貯留層属性の前記初期コンピュータ行列が対角優位でない場合、近似解析プレコンディショナ（前処理法）を生成するステップと、
   （ｈ）生成された近似解析前処理法を既知の貯留層属性の前記初期コンピュータ行列に適用して、前記貯留層のグリッドセル内の圧力分布の測定値を形成するステップと、
   （ｉ）前記貯留層のグリッドセル内の圧力分布の形成された測定値からＫｒｙｌｏｖベクトルを生成するステップと、
   （ｊ）前記圧力分布の形成された測定値から生成されたＫｒｙｌｏｖベクトルから、貯留層圧力分布の係数ベクトルを決定するステップと、
   （ｋ）貯留層圧力分布の決定された係数ベクトルに基づいて、前記貯留層シミュレーションの時間ステップに対する、前記貯留層のグリッドセル内の推定流体流量を決定するステップと、
   （ｌ）前記貯留層圧力分布の決定された係数ベクトルに基づいて、前記貯留層のグリッドセル内の推定流体流量が収束したかどうかを決定するステップと、
   （ｍ）収束した場合、前記時間ステップについて、前記貯留層内の流体流の前記貯留層シミュレーションを終了するステップと、
   のコンピュータ処理ステップを含む、コンピュータ実施方法。
2. 前記貯留層シミュレーションに基づく前記油井のうちの少なくとも１つの生産を調整するステップをさらに含む、請求項１に記載の方法。
3. 前記貯留層シミュレーションに基づく前記油井のうちの少なくとも１つの性能を調整するステップをさらに含む、請求項１に記載の方法。
4. 前記決定された流体流量が収束した場合に前記貯留層シミュレーションを終了する後続の時間ステップにインクリメントするコンピュータ処理ステップをさらに含む、請求項１に記載の方法。
5. 前記後続の時間ステップについて、コンピュータ処理ステップ（ａ）〜（ｍ）を繰り返すステップをさらに含む、請求項４に記載の方法。
6. 生成されたＫｒｙｌｏｖベクトルを直交Ｋｒｙｌｏｖベクトルに変換するコンピュータ処理ステップをさらに含む、請求項１に記載の方法。
7. 前記直交Ｋｒｙｌｏｖベクトルの直交正規化を実行するコンピュータ処理ステップをさらに含む、請求項６に記載の方法。
8. 貯留層圧力分布の係数ベクトルを決定する前記コンピュータ処理ステップは、直交Ｋｒｙｌｏｖベクトルの直交正規化を実行するステップの後に、貯留層圧力分布の係数ベクトルを決定するコンピュータ処理ステップを含む、請求項７に記載の方法。
9. 貯留層圧力分布の係数ベクトルを決定する前記コンピュータ処理ステップは、線形システムのガウス消去法を実行することによって係数ベクトルの前記線形システムを処理するコンピュータ処理ステップを含む、請求項１に記載の方法。
10. 貯留層からの予測生産の時間ステップのシーケンスについて、前記貯留層内の圧力分布を決定して、前記貯留層内のシミュレートされた流体流を決定する、収束性を改善した生産貯留層内の流体流の貯留層シミュレーションを実行するデータ処理システムであって、前記貯留層は、貯留層セルの３次元グリッドに組織化され、前記貯留層シミュレーションは、前記貯留層内の油井からの実際の流体圧力および流体生産に基づいて、時間ステップのシーケンスについて実行され、前記コンピュータ実施方法は、前記貯留層セルの３次元グリッドから油井へのシミュレートされた流体流および前記油井のシミュレートされた油井生産速度を決定し、前記データ処理システムが、
    （ａ）プロセッサであって、
    （１）時間ステップのシーケンスにおける、ある時間ステップについて、既知の貯留層属性の初期コンピュータ行列と実際の貯留層流体生産測定値の初期コンピュータベクトルとを形成するステップと、
    （２）既知の貯留層属性の前記初期コンピュータ行列が対角優位であるかどうかを判定するステップと、
    （３）対角優位である場合、前記貯留層のグリッドセルの個々の中の圧力分布の形成された初期測定値に基づいて、貯留層シミュレーションの時間ステップについて、前記貯留層のグリッドセルの個々の中の圧力分布を決定するステップと、
    （４）前記貯留層のグリッドセル内の決定された推定圧力分布に基づいて、前記貯留層シミュレーションの時間ステップについて、前記貯留層のグリッドセル内の推定流体流量を決定するステップと、
    （５）前記貯留層のグリッドセル内の決定された流体流量が収束したかどうかを決定するステップと、
    （６）収束した場合、前記時間ステップについて、前記貯留層内の流体流の貯留層シミュレーションを終了するステップと、
    （７）前記既知の貯留層属性の初期コンピュータ行列が非対角優位である場合、近似解析プレコンディショナ（前処理法）を生成するステップと、
    （８）生成された近似解析前処理法を前記既知の貯留層属性の初期コンピュータ行列に適用して、前記貯留層のグリッドセル内の圧力分布の測定値を形成するステップと、
    （９）前記貯留層のグリッドセル内の形成された圧力分布の測定値からのＫｒｙｌｏｖベクトルを生成するステップと、
    （１０）形成された圧力分布の測定値から生成されたＫｒｙｌｏｖベクトルから貯留層圧力分布の係数ベクトルを決定するステップと、
    （１１）貯留層圧力分布の決定された係数ベクトルに基づいて、前記貯留層シミュレーションの時間ステップについて、前記貯留層のグリッドセル内の推定流体流量を決定するステップと、
    （１２）前記貯留層圧力分布の決定された係数ベクトルに基づいて、前記貯留層のグリッドセル内の推定流体流量が収束したかどうかを決定するステップと、
    （１３）収束した場合、前記時間ステップについて、前記貯留層内の流体流の前記貯留層シミュレーションを終了するステップと、
    のコンピュータ処理ステップを実行する、プロセッサと、
    （ｂ）前記貯留層シミュレーション結果の表示を提供するためのワークステーション
    からなる、データ処理システム。
11. 前記貯留層シミュレーション結果を記憶するメモリをさらに含む、請求項１０に記載のデータ処理システム。
12. 前記プロセッサが、前記決定された流体流量が収束した場合に前記貯留層シミュレーションを終了する後続の時間ステップにインクリメントすることをさらに含む、請求項１０に記載のデータ処理システム。
13. 前記プロセッサが、前記後続の時間ステップについてコンピュータ処理ステップ（ａ）〜（ｍ）を繰り返すことをさらに含む、請求項１２に記載のデータ処理システム。
14. 前記プロセッサが、生成されたＫｒｙｌｏｖベクトルを直交Ｋｒｙｌｏｖベクトルに変換することをさらに含む、請求項１０に記載のデータ処理システム。
15. 前記プロセッサが、直交Ｋｒｙｌｏｖベクトルの直交正規化を実行することをさらに含む、請求項１４に記載のデータ処理システム。
16. 貯留層圧力分布の係数ベクトルを決定する際、前記プロセッサが、前記直交Ｋｒｙｌｏｖベクトルの直交正規化を実行した後に、貯留層圧力分布の係数ベクトルを決定する、請求項１５に記載のデータ処理システム。
17. 貯留層圧力分布の係数ベクトルを決定する際、前記プロセッサが、前記線形システムのガウス消去法を実行することによって、係数ベクトルの線形システムを処理する、請求項１５に記載のデータ処理システム。

Images