CN112752894A - 具有用于非对角占优的不定系数矩阵的压力求解器的储层模拟 - Google Patents

具有用于非对角占优的不定系数矩阵的压力求解器的储层模拟 Download PDF

Info

Publication number
CN112752894A
CN112752894A CN201980062351.6A CN201980062351A CN112752894A CN 112752894 A CN112752894 A CN 112752894A CN 201980062351 A CN201980062351 A CN 201980062351A CN 112752894 A CN112752894 A CN 112752894A
Authority
CN
China
Prior art keywords
reservoir
computer
simulation
pressure distribution
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201980062351.6A
Other languages
English (en)
Other versions
CN112752894B (zh
Inventor
阿里·海达尔·多鲁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Saudi Arabian Oil Co
Original Assignee
Saudi Arabian Oil Co
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Saudi Arabian Oil Co filed Critical Saudi Arabian Oil Co
Publication of CN112752894A publication Critical patent/CN112752894A/zh
Application granted granted Critical
Publication of CN112752894B publication Critical patent/CN112752894B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V20/00—Geomodelling in general
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B2200/00—Special features related to earth drilling for obtaining oil, gas or water
    • E21B2200/20—Computer models or simulations, e.g. for reservoirs under production, drill bits
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B2200/00—Special features related to earth drilling for obtaining oil, gas or water
    • E21B2200/22—Fuzzy logic, artificial intelligence, neural networks or the like
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Mining & Mineral Resources (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Geochemistry & Mineralogy (AREA)
  • Fluid Mechanics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

在通过计算机处理进行储层模拟期间确定储层模型的网格单元之间的压力分布期间，改进了计算机的性能。当条件导致处理中涉及的系数矩阵变得不定时，可以证明难以获得收敛。不定系数矩阵的出现可能是由于与汽液平衡有关的储层中的物理条件导致的，或者是由于不适当的导数而在数值上产生的非物理条件导致的。当在储层模拟期间收敛变得困难时，避免了传统的先前采取的时间步长切割的校正动作。对于具有涉及数百万、数十亿或更多数量的未知参数值的模型的非常大的储层，在计算机使用和时间方面，时间步长切割已经被证明是非常昂贵的。

Description

具有用于非对角占优的不定系数矩阵的压力求解器的储层 模拟
技术领域
本发明涉及储层模拟，并且尤其涉及在流体生产期间考虑储层中的流体流动和压力分布的模拟方法。
背景技术
在油气工业中，为了勘探和生产的目的，需要处理大量的数据进行计算机化的模拟、建模和分析。例如，地下烃储层的开发通常包括储层的计算机模拟模型的开发和分析。储层的实际模拟模型及其流体的存在有助于预测来自烃储层的最佳未来油气采收。油气公司已经开始依赖于模拟模型，作为增强开发石油储备的能力的重要工具。
地下烃储层通常是含有石油流体混合物和水的复杂岩层。储层流体含量通常存在于两个或更多个流体相中。储层流体中的石油混合物通过钻入这些岩层并在其中完井的井来生产。有时，诸如水或气体或两者的流体也被注入这些岩层中以提高石油流体的采收。
储层模拟属于多孔介质模拟中流动(flow)的一般领域。然而，储层模拟通常涉及处于高压和高温下的地下地质构造中的多种烃组分和多个流体相。这些烃流体和所包含的地下水的化学相行为必须在这些模拟器中考虑。
模拟模型包含描述岩层和井的特定几何形状的体积数据，以及诸如流体和岩石性质的储层性质数据，以及与所讨论的油田或气田的特定储层有关的生产和注入历史。模拟模型由模拟器(称为储层模拟器)形成，该模拟器是在数据处理系统S上运行的一套计算机程序。
运行这些模型的储层模拟器是计算机实现的数值方法，或编码算法和基础数学模型的数据结构。表示这些烃储层中流体运动的物理性质的数学模型是非线性偏微分方程组，其描述了由流体的生产注入引起的这些储层中的瞬时多相、多成分流体流动和物质平衡行为，以及储层流体的压力体积温度(PVT)关系。
通过将体积细分为连续单元(也称为网格块)，储层模拟器模拟地下储层和所包括的周围多孔岩层中的多相多成分流体流动和物质平衡。在模拟模型中，储层因此被组织成多个单独的单元。单元或网格块是应用基础数学模型的基本有限体积。单元的数量根据模拟所需的分辨率和所讨论的储层的尺寸而变化。
对于大的储层，例如在工业中已知为巨型储层的类型，其可以具有数十亿桶的原始石油地质储量(OOIP)，网格单元的数量可以是数亿到十亿以上。需要这一数量的单元以便具有足够的分辨率来表示流动动力学、地层岩石孔隙率和渗透率异质性以及储层内的许多其它地质和沉积复杂性。这种尺寸的储层的模拟可以被称为千兆单元储层模拟。图1在100处示出了用于储层模拟的地下储层的示例结构化储层网格模型。
烃储层模拟中的挑战要求使用最新技术以成本有效的方式使采收最大化。诸如以商标GigaPOWERS提供的储层模拟已经在文献中被描述。例如，参见以下文章：Dogru,A.H等人的“A Next-Generation Parallel Reservoir Simulator for Giant Reservoirs,”SPE119272,proceedings of the 2009SPE Reservoir Simulation Symposium,TheWoodlands,Texas,USA,Feb 2-4,2009；以及Dogru,A.H；Fung,L.S；Middya,U.；Al-Shaalan,T.M.；Byer,T.；Hoy,H.Hahn,W.A.；Al-Zamel,N.；Pita,J.；Hemanthkumar,K.；Mezghani,M.；Al-Mana,A.；Tan,J.；Dreiman,T.；Fugl,A.；Al-Baiz,A.的“New Frontiers in LargeScale Reservoir Simulation,”SPE 142297,Proceedings of the 2011 SPE ReservoirSimulation Symposium,The Woodlands,Texas,USA,Feb 21-23,2011。GigaPOWERS储层模拟能够进行超过十亿单元障碍的精细尺度网格模拟以进行后处理，同时每种情况利用数百千兆字节(GB)覆盖区(footprint)。
具有多个烃储层和可观储量的公司的模拟运行的总数每年超过数万，并且需要一个或多个千万亿字节(petabyte)的高性能存储器来容纳这些数据。
多相多成分系统的瞬态解涉及从储层的初始条件开始的一系列时间步长中的质量和能量守恒的演变。对于每个时间步长，使用所谓的广义牛顿法来线性化每个有限体积的非线性离散方程组。
在工业中，这被称为牛顿迭代或非线性迭代。在每个牛顿迭代中，构建线性方程组，其中使用被称为雅可比矩阵的矩阵和被称为残差的右端向量来求解方程组的主要变量的变化。
关于储层的其它因素进一步增加了千兆级储层模拟的复杂性。地层中的地下岩石在渗透率或孔隙率方面不是同质的。对于高度异质的储层，迭代方法在计算储层压力分布时可能不收敛。这就给模拟处理引入了额外的复杂性。
用于求解线性方程组的当前工业实践通过预调节迭代法。在例如Saad的“Iterative Methods for Sparse Linear Systems,”January,2000,Ch.10,Pp.265–319中描述了预调节的示例。在计算上昂贵的预调节器方法增加了计算成本。此外，预调节法不保证收敛到精确的解。
预调节器通常可以提供对储层内未知压力和饱和度分布的良好估计。然后通过迭代将该初始估计改进为针对该时间步长的正确压力和饱和度分布。一种用作预调节器的方法被称为高斯塞德尔方法(Gauss Seidel method)，由于其简单性和每次迭代的较低的比较计算成本，使用该方法。高斯塞德尔预调节可以采用被称为点高斯塞德尔或线性高斯塞德尔预调节器的形式。确定用于下一个模拟器时间步长的储层矩阵的网格单元压力的直接解，并且评估残差。这个处理继续，直到所计算的网格的残差小于指定的公差。在Olsen-Kettle的“Numerical Solution of Partial Differential Equations”,March,2011,Chapter 4.2,Pp.28–33,The University of Queensland中描述了示例高斯塞德尔方法。然而，已经发现，在获得网格单元压力时，高斯塞德尔预调节经常显示出非常慢的收敛，其中在指定的公差内有残差，尤其是对于具有高度异质地层岩石特性的储层。
压力方程在储层模拟中起重要作用。在储层模拟期间确定的压力条件表示在储层的模型中的储层网格单元的所关注的指定时间步的压力条件。在储层模拟处理期间，由于物理或数值处理原因，系数矩阵在模拟期间可能变得非对角占优或不定。这导致迭代方法发散。就目前所知，当前的实践是减小时间步长大小，并重建系数矩阵，并以大小减小的时间步长重复处理。在这些情况下，模拟处理在针对大量未知量的计算机操作时间和成本方面可能非常昂贵，这是对于任何可调大小的储层的情况，并且特别是对于巨大的储层的情况。
发明内容
简言之，本发明提供一种对生产储层中的流体流动进行储层模拟的计算机实现方法，所述计算机实现方法具有确定所述储层内的压力分布的改进的收敛，以针对从所述储层预计生产的时间步长序列，确定所述储层中的模拟流体流动，所述储层被组织成储层单元的三维网格，针对所述时间步长序列，基于实际流体压力和来自所述储层中的井的流体生产来执行所述储层模拟，所述计算机实现方法确定从储层单元的所述三维网格到所述井的模拟流体流动和所述井的模拟井生产速率。
该计算机实现方法对于所述时间步长序列中的时间步长，构成已知储层属性的初始计算机矩阵和实际储层流体生产测量的初始计算机向量。接着，确定已知储层属性的所述初始计算机矩阵是否对角占优。
如果对角占优，则基于所述储层的所述网格单元中的各个网格单元中的压力分布的所构成的初始测量，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的所述各个网格单元中的压力分布。然后，基于所述储层的所述网格单元中的所确定的估计压力分布，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的估计流体流速。然后，计算机确定所述储层的所述网格单元中的所确定的流体流速是否已收敛；以及如果已收敛，则结束针对所述时间步长的所述储层中的流体流动的所述储层模拟。
如果已知储层属性的所述初始计算机矩阵不是对角占优的，则生成近似分析预调节器。然后，将所生成的近似分析预调节器应用到已知储层属性的所述初始计算机矩阵，以构成所述储层的所述网格单元中的压力分布的测量。根据所述储层的所述网格单元中的压力分布构成的测量，生成Krylov向量，以及基于根据压力分布构成的测量所生成的Krylov向量，确定储层压力分布的系数向量。
基于所确定的储层压力分布的系数向量，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的估计流体流速。然后，确定基于所确定的储层压力分布的系数向量的所述储层的所述网格单元中的所述估计流体流速是否已收敛。如果已收敛，则结束针对所述时间步长的所述储层中的流体流动的储层模拟。
本发明还提供了一种新的改进的数据处理系统，其执行生产储层中的流体流动的储层模拟，所述数据处理系统具有确定所述储层内的压力分布的改进的收敛，以针对从所述储层预计生产的时间步长序列，确定所述储层中的模拟流体流动，所述储层被组织成储层单元的三维网格，针对所述时间步长序列，基于实际流体压力和来自所述储层中的井的流体生产来执行所述储层模拟。
该数据处理系统包括处理器，该处理器确定从储层单元的三维网格到井的模拟流体流动和井的模拟井生产速率。处理器对于所述时间步长序列中的时间步长，构成已知储层属性的初始计算机矩阵和实际储层流体生产测量的初始计算机向量。接着，处理器确定已知储层属性的所述初始计算机矩阵是否对角占优。
如果对角占优，则基于所述储层的所述网格单元中的各个网格单元中的压力分布的所构成的初始测量，针对所述储层模拟的所述时间步长，由处理器确定所述储层的所述网格单元中的所述各个网格单元中的压力分布。然后，基于所述储层的所述网格单元中的所确定的估计压力分布，针对所述储层模拟的所述时间步长，由处理器确定所述储层的所述网格单元中的估计流体流速。然后，处理器确定所述储层的所述网格单元中的所确定的流体流速是否已收敛；以及如果已收敛，则通过处理器结束针对所述时间步长的所述储层中的流体流动的所述储层模拟。
如果已知储层属性的所述初始计算机矩阵不是对角占优的，则由处理器生成近似分析预调节器。然后，由处理器将所生成的近似分析预调节器应用到已知储层属性的所述初始计算机矩阵，以构成所述储层的所述网格单元中的压力分布的测量。由处理器根据所述储层的所述网格单元中的压力分布构成的测量生成Krylov向量，以及基于根据压力分布构成的测量所生成的Krylov向量，确定储层压力分布的系数向量。
基于所确定的储层压力分布的系数向量，针对所述储层模拟的所述时间步长，由处理器确定所述储层的所述网格单元中的估计流体流速。然后，确定基于所确定的储层压力分布的系数向量的所述储层的所述网格单元中的所述估计流体流速是否已收敛。如果已收敛，则通过处理器结束针对所述时间步长的所述储层中的流体流动的储层模拟。
附图说明
图1是用于储层模拟的地下储层结构化网格的计算机化模型的等距视图。
图2是简化的示例二维储层模型的一部分的平面图。
图3是根据本发明的针对储层模拟的用于确定异质岩层中的压力分布的计算机网络的示意图。
图4是图3的计算机网络的应用服务器或计算机节点的示意图。
图5、图6和图7是根据本发明的在存在不定储层系数矩阵的情况下针对储层模拟的用于确定异质岩层中的压力分布的工作流程的功能框图或流程图。
图8是用于获得根据本发明与现有技术方法比较的压力分布的比较结果的简化的储层模型的示意图。
图9是在图8的储层模型的压力分布确定期间的收敛误差的比较结果的显示。
图10是用于获得根据本发明与现有技术方法比较的压力分布的比较结果的简化的储层模型的示意图。
图11是根据本发明的图10的储层模型的与现有技术方法相比的在压力分布确定期间的收敛误差的比较结果的显示。
图12是用于获得根据本发明与现有技术方法比较的压力分布的比较结果的简化的储层模型的示意图。
图13是根据本发明的压力分布确定中的收敛误差与根据现有技术的处理结果的比较结果的显示。
具体实施方式
在数值储层模拟中，基于公式(formulation)，明确地形成压力方程，或者将压力方程包括在流动方程的主体中。根据压力分布计算流体速度以及因此计算油、水和气体的流量。因此，正确地计算压力场是重要的。
在储层模拟中，流动方程被线性化，并且系数矩阵被形成以求解未知量。该系数矩阵通常被称为雅可比矩阵。系数或雅可比矩阵包含诸如压力的未知量相对于通过储层模拟确定的变量的一阶导数。下面阐述的方程(1)中示出了示例雅可比矩阵。
地下流体储层在岩石性质上是高度异质的。这些异质性反映在描述储层中的流体流动和压力分布的偏微分方程的系数中。
由于在模拟处理期间在一个时间步长上难以找到适当的压力解，模拟器经常面临收敛问题。因此，正确计算的压力解对于模拟器在时域中进行并以合理的计算成本确定储层性能是关键的。
本发明已经被开发用于储层模拟器的压力方程，以针对在网格单元中的压力条件的确定所基于的参数中存在的不定系数矩阵收敛。不定系数矩阵可以由于模拟与汽液平衡相关的物理条件而出现，或者由于计算的不适当导数而在数值上产生的非物理条件而出现。
以前的迭代方法对于不定系数矩阵有所差异。本发明在遇到不定系数矩阵的情况下收敛。本发明允许储层模拟过程继续而不减少时间步长。在计算机处理时间和储层模拟运行方面，时间步长的减小是昂贵的，其中存在大量的未知量，大约数百万和数十亿以及更多。本发明避免了时间步长切割，对于包含百万、十亿或更多数量的未知量的非常大的方程组，时间步长切割可能是非常昂贵的。本发明通过新的迭代方法确定具有不定储层系数矩阵的压力条件，并允许储层模拟器继续模拟而不切割时间步长大小。
压力方程在储层模拟中模拟流体流动方面起到重要作用。根据所选择的数值求解方案，可以形成单独的压力方程，或者与其它储层参数未知量同时求解相同的未知量。根据本发明，储层模拟器形成单独的压力方程，并且通过求解以如方程(1)所示的一般矩阵形式表示的线性方程组，来求解储层模拟的下一个时间步长压力：
[A][p]＝[b] (1)
对于方程(1)，在储层的压力分布确定的情况下，A是已知的系数矩阵，b是实际储层生产测量的已知右端向量，p是下一个时间步长的未知压力向量。系数矩阵A是由描述流体流动的偏微分方程的离散化产生的稀疏n×n矩阵，如将描述的。在方程(1)中，b是n×1行向量，其通常包括旧时间步长信息和井生产/注入速率，其中n是未知压力向量的总数。
术语
为了便于参考和理解，下面列出了方程中使用的术语的列表，其表示根据本发明的各种参数、测量和测量单位。
P＝压力，磅/平方英寸(psi)
Figure BDA0002988214720000081
pmax＝最大压力(psi)
t＝时间(天)
Δt＝时间变化(天)
L＝距离，英尺(ft.)
x＝x方向上的距离(ft.)
y＝y方向上的距离(ft.)
Δx＝x方向上的网格块尺寸
Δy＝y方向上的网格块尺寸
Δz＝z方向上的网格块尺寸
q＝生产率，桶每天(b/d)
r＝离点x*的径向距离(ft.)
Vp＝孔隙体积，桶(1桶＝5.615立方英尺，1ft.＝30.48厘米)
CT＝总压缩系数(1/psi)
Co＝油压缩系数(1/psi)
Cw＝水压缩系数(1/psi)
Cg＝气压缩系数(1/psi)
Cp＝孔隙压缩系数(1/psi)
Figure BDA0002988214720000091
K＝渗透率(mD)
μ＝粘度，厘泊(cP)
<x，y>＝xTy＝向量x和向量y的内积
微分形式和离散化形式的压力方程
作为一个示例，考虑具有恒定流体性质的方形扁平单相油储层V(图2)，在油储层中的流体的泡点以上操作。描述线性化的、可略微压缩的单相流体的储层压力p的偏微分方程如下给出：
Figure BDA0002988214720000092
其中p是由于在储层中操作的具有流动强度q的源或水槽(井)而导致的未知压力分布。参数k表示储层模型100中在x、y、z维度坐标方向的变化的渗透率。参数Ct是总的流体和岩石压缩系数，参数t是时间。该表达式的形式与定义体积平衡关系的表达式的形式相当。
方程(2)基于对于储层中的压力分布的给定的外部和内部边界条件以及初始条件，存在p(x,y,z,t)的唯一解。
对于简化的示例(仅考虑x和y中的两个维度)，方程(2)变为
Figure BDA0002988214720000093
假设这种简化的二维储层模型，例如储层模型V(图2)，被组织或划分为网格块(控制体积)。对于更一般化的示例，储层模型由x维中的数量为Nx的网格块和y维中数量为Ny的网格块组成。二维储层中的网格块的总数为n，其中n＝Nx×Ny。
在图2的模型V中，i和j分别被指定为x和y方向上的编号序列索引。在x方向上，i＝1，2,......,Nx；以及在y方向上，j＝1，2,......,Ny。下面的方程(4)对于模型V表示用于方程(3)的微分方程的离散化方程组的有限差分形式。为了求解方程(3)中的未知函数p，模型V的x-y域被分成图2所示的计算元素。方程(4)可以表示如下：
Figure BDA0002988214720000101
其中T是描述网格块(单元)(i,j)与其相邻网格块之间的链接流动的传输率。储层模型中的传输率是表示针对流动流体的粘度校正的地层网格单元元素的电导率的度量。
Figure BDA0002988214720000104
是元素(i,j)的孔隙体积。ct是该元素(i,j)的总岩石和流体压缩系数。
在方程(4)中，上标n表示旧的时间步长，上标(n+1)表示新的时间步长，△t是新时间t(n+1)与旧时间t(n)之间的时间步长大小，且Vp(i，j)是网格块的孔隙体积。
在储层模型V作为示例的情况下，坐标(i,j)处的单元110的西部传输率由Ti-1/2，j表示，其中，(i-1/2,j)是指坐标(i,j)处的网格块110的西面112。类似的，Ti+1/2，j是同一单元110的东面114的传输率。以可比较的方式，单元110的北面116的传输率由Ti，j+1/2表示，且单元110的南面118的传输率由Ti，j-1/2表示。
通常，网格中的面在x维度上到单元块的传输率T被定义为：
Figure BDA0002988214720000102
其中
Figure BDA0002988214720000103
是在网格块(i+1,j)和(i,j)之间的谐波平均渗透率，A是网格块的界面面积，并且x是网格块的x坐标。在方程(5)中，C1＝1.127E-3，其是单位转换因子。转换因子C1将公制度量单位为立方厘米每秒的流体流速转换成实际的油田单位(桶每天)。因此，当储层模型中的x或y距离以ft表示时，传输率T以体积、时间和桶/天/psi的压力度量单位表示；面积A以ft2表示；而渗透率k以度量单位毫达西而不是达西表示。
储层中每个网格块在时间步长n处的压力分布作为前一个时间步长的结果，是已知的。通过针对单元写方程(4)并重新排列成方程(6)的矩阵形式，构成线性方程组，可以获得每个网格块在新时间步长处的未知压力。
重新排列用于网格块(i,j)的方程(4)得到如下方程(6)：
Figure BDA0002988214720000111
为了方便起见，方程(6)中的东面指示符(i+1/2,j)已经被缩写“E”代替以指示东面；西面指示符(i-1/2,j)被缩写“W”代替以指示西面；北面指示符(i,j+1/2)被缩写“N”代替以指示北面；以及南面指示符(i,j-1/2)被缩写“S”代替以指示南面。
然后，对于在时间步长(n+1)的未知压力向量的矩阵形式的方程(6)可以以矩阵形式表示为以下形式的方程(7)：
Figure BDA0002988214720000112
与方程(6)一样，为了方便起见，方程(7)中的东面指示符(i+1/2,j)已经被缩写“E”代替以指示东面；西面指示符(i-1/2,j)被缩写“W”代替以指示西面；北面指示符(i,j+1/2)被缩写“N”代替以指示北面；以及南面指示符(i,j-1/2)被缩写“S”代替以指示南面。在油田单位中，方程(7)是以桶每天每psi压力变化，或立方厘米每秒每大气压差，为单位。
在方程(7)的右侧，qi是给定网格块处的生产或注入速率，并且项
Figure BDA0002988214720000113
是网格块k的孔隙体积vp、总压缩系数CT与每时间步长大小△t的网格块k的旧时间步长压力pk的乘积。
方程(7)中的线性方程组以方程(1)的矩阵形式组织，根据术语部分定义的参数表示为：
Figure BDA0002988214720000121
其中，在方程(8)中，[T]是与方程(1)中的系数矩阵[A]对应的流体传输率方面的系数矩阵；
Figure BDA0002988214720000122
是方程(8)的右端向量，其包括与方程(1)中的右端矩阵[b]对应的源项、孔隙体积、总压缩系数和时间步长大小；p是时间步长(n+1)处的未知压力向量。
储层模型
实际示例地下储层的结构化储层网格模型(例如图1中的100所示)由网格块(计算元素)组成，其中，每个元素具有单个未知压力p。储层模型被分成数百万个元素。因此，矩阵A可以非常大，即，百万乘百万个稀疏矩阵。由于储层模型的系数矩阵A的尺寸较大，所以例如通过称为高斯消元法的过程的直接求解方法在计算上非常昂贵且缓慢。因此，通常使用迭代求解方法来求解方程(1)。
迭代方法的挑战
如果矩阵是强对角占优的，则迭代方法工作良好。这也意味着矩阵的所有特征值具有相同的符号(正或负)。当矩阵失去对角占优时(当特征值改变符号时)，方程(1)中的系数矩阵A变得不定。就目前所知，迭代方法针对不定矩阵不收敛。
矩阵术语和迭代方法
矩阵的特征值
根据方程(1)的具有n行和n列的方阵[A]的特征值被定义为以下矩阵[A]的行列式的关系表达式的根，如下：
det(A-λI)＝0 (9)
其产生如下定义的n次多项式：
λn+c1λn-1+…cn-1＝0 (10)
矩阵[A]的特征值是满足方程(10)的多项表达式的根。系数矩阵是在方程(1)中用在储层模拟器中以求解储层中的压力分布的同一矩阵。λ表示的矩阵[A]的特征值包括变量，如方程(5)、方程(6)限定的传输率T；元素(i,j)的孔隙体积
Figure BDA0002988214720000131
和总压缩系数Ct。
正定矩阵和不定矩阵
正定矩阵是其中这种矩阵的特征值中的每一个大于零的矩阵。负定矩阵是其中矩阵A的每个特征值小于零的矩阵。不定矩阵是其一些特征值为正而一些特征值为负的矩阵。如果矩阵的一个特征值在符号上与该矩阵的其余特征值不同，则该矩阵是不定的。
实际上，在模拟运行期间，如果任何网格块的总压缩系数变为负，则系数矩阵A变得不定，并且储层模拟的迭代方法发散。如上所述，网格块的负总压缩系数可能是由于物理原因或某些相平衡计算中岩石和流体导数值的组合导致的。
用于具有储层模拟的迭代方法的收敛问题的不定矩阵的示例
压力方程的方程(8)的系数矩阵[T]通常是对角占优的，因此是正定的。因此，任何迭代方法将以不同的收敛速率收敛到实解。收敛性受到系数矩阵A的元素的大小影响。当矩阵的一行中的所有元素的绝对值的总和远小于同一行的对角项的绝对值(数)时，系数矩阵A是强对角占优的。
在非对角元素的和与对角元素的和之间的差较小的情况下，则该矩阵被称为弱对角占优的。强对角占优矩阵在几个计算机处理迭代内收敛到实解。另一方面，弱对角占优矩阵需要许多迭代来收敛。然而，这些类型的矩阵中的每一个在一些情况下失去一个或多个矩阵行或网格块中的对角优势。在这种情况下，系数矩阵[T]可能变得不定。当存在不定系数矩阵时，典型的迭代方法将不收敛。
回到方程(6)，其将用于被重新布置成矩阵形式的单元的每个网格块线性方程组的未知压力表达为：
Figure BDA0002988214720000141
在方程(6)中，元素(i,j)的矩阵行被表示为：
Figure BDA0002988214720000142
由于所有的传输率值T＞0，因此对于
Figure BDA0002988214720000143
矩阵[T]中的对角项的绝对值大于所有非对角线项的和；或
Figure BDA0002988214720000144
如果这个条件对于矩阵[T]的所有行(所有网格块(i，j))都为真，则系数矩阵[T]是正定的，并且特征值具有相同的符号。
在如方程(4)中所表示的
Figure BDA0002988214720000145
在压力方程矩阵中为负的情况下，总压缩系数为负，即Ct＜0。系数矩阵变得不定，其中至少一个特征值具有与其它特征值相反的符号。当这种情况发生时，传统的迭代方法不收敛。
如果在模拟运行期间出现发散或不收敛，则模拟器通常切割时间步长大小并形成新的系数矩阵T。形成新的系数矩阵T通常是非常昂贵的，尤其是对于包括具有未知压力分布的大量网格块的储层模型。因此，重要的是不切割时间步长大小。本发明提供了一种使用相同系数矩阵T的收敛方法，即使在某些情况下它可能是不定的。
为了生成方程(8)形式的解，必须对系统进行预调节。可以使用高斯塞德尔迭代预调节，例如在前面提到的Olsen-Kettle的“Numerical Solution of PartialDifferential Equations”中描述的。
由于预调节器也是迭代方法，因此它无助于求解方程(1)。这是因为迭代计算过程在存在不定系数矩阵的情况下将发散。因此，根据本发明，当存在一个或多个不定系数矩阵时，近似分析矩阵[B]预调节器(近似分析解)确定网格单元压力分布。如将要描述的，在根据图7的处理的步骤228期间执行该处理。
本发明在确定储层模型的网格单元之间的压力分布期间改进了数据处理系统S(图3)的性能。当条件导致存在不定系数矩阵时，可能证明难以获得收敛。迄今为止，已知当存在不定系数矩阵时，用于确定压力分布的先前迭代方法已经发散。
避免了传统的储层模拟期间的时间步长切割对发散的处理动作。本发明允许储层模拟过程继续而不减少时间步长大小。对于具有涉及数百万、数十亿或更多数量的未知参数值的模型的非常大的储层，在计算机使用和时间方面，时间步长切割已经被证明是非常昂贵的。
利用本发明，通过计算机处理求解矩阵形式的方程(1)来确定压力分布可以表示为被称为Krylov向量的线性组合，表示如下：
Figure BDA0002988214720000151
在方程(13)中，表示为总和的向量数量的上限，即向量数量m，理想地等于待确定的储层网格元素中的压力分布的未知值的总数。然而，如将在稍后部分中解释的，已经发现未知值的总数远小于Krylov向量的数目。
随后的部分解释了如何构造线性独立向量vi作为Krylov向量，并选择系数ci，以通过根据本发明的计算机处理确定压力分布。
近似分析解用作预调节器以完全相同的方式求解不定矩阵系统，而不需要直接分析解方法。如将解释的，这通过几个实例来证明。
数据处理系统S
根据本发明的处理和工作流程适于部署在各种高性能计算机(或HPC)集群计算机/处理器节点计算机系统上。这些通常是具有若干计算节点的机架安装的硬件，其包含具有多核架构的多个CPU。节点通常与传统的低延迟高带宽网络、交换机和路由器互连。
应当理解，本发明可以在具有任何常规类型的适当处理能力的大型计算机(例如可以在从纽约阿蒙克的国际商业机器公司(IBM)或其它来源购得的大型计算机)上执行。
用于与该模拟系统一起使用的典型HPC环境是当今的多节点、多CPU、多核计算集群。在图3的数据处理系统S中，在C处示出了这样的集群的示例。集群C由多个计算机节点150(图3和4)形成，(一个或多个)路由器服务器154如箭头152所示并行地向这些计算机节点提供数据。如果需要，可以使用几个这样的路由器服务器来实现这个目的。上述类型的原始模拟或输入数据被存储在适当数量的数据存储/文件服务器156中。数据存储/文件服务器156存储操作指令、控制信息和数据库记录，以供数据处理系统S使用。
路由器服务器154在存储在存储器中的计算机代码155的控制下操作，并且向集群C的计算机节点150并行地传送来自存储服务器156的输入模拟数据以及如箭头158所示的模拟处理结果，以及传送来自集群C的计算机节点150的所述输入模拟数据以及模拟处理结果。
根据本发明的程序代码155是非暂态计算机可操作指令的形式，其使得(一个或多个)服务器154索引、排序、处理和传送储层模拟结果。典型地，数据处理系统S包括一组适当的、传统类型的工作站157，其通过网络159连接到系统。工作站157提供输出数据和显示，使得储层工程师观察储层模拟结果。基于储层模拟结果，工程师能够根据需要评估和采取行动以监测和调整井和储层的生产和性能。
集群C的计算机节点150包括图4所示类型的多个处理器或核160，其在存储于计算机节点150的存储器164中的计算机代码或程序产品162的指令下并行操作。根据本发明的程序代码162是非暂态计算机可操作指令的形式，其使得数据处理器160在储层模拟模块166中执行储层模拟。
应当注意，程序代码155和162可以是微码、程序、例程或符号计算机可操作语言的形式，其提供控制数据处理系统S的功能和指导其操作的特定的一组有序操作。程序代码155和162的指令可以存储在服务器154或处理器节点150的存储器中，或者存储在计算机磁盘、磁带、常规硬盘驱动器、电子只读存储器、光学存储装置或其上存储有非暂时性计算机可用介质的其他适当的数据存储装置上。程序代码162也可以包含在诸如服务器156的数据存储装置上，作为计算机可读介质，如图所示。
由数据处理系统S进行的处理
本发明提供了对储层(例如储层100(图1))的异质岩层的三维上的压力分布和压力值的改进的确定。当系数矩阵变得不定并导致传统的预调节不收敛时，本发明还允许在地下储层100的储层模拟期间通过储层模型进行用于压力确定的即时初始压力测量。
本发明避免了在那些情况下储层模拟期间的收敛问题，而不需要储层模拟的时间步长大小的不期望的减小。如已经讨论的，时间步长大小的减小增加了计算机处理操作和复杂性以及储层模拟处理时间长度。所确定的初始压力分布表示储层模拟的储层内的正确压力分布曲线。
通过示例证明初始压力的估计接近基于简化储层模型中的分析解的实际解。所确定的初始压力估计因此使得储层模拟器能够在指定的精度限制内收敛到期望的解。根据本发明的初始压力分布的确定还以比其它预调节方法(如点高斯塞德尔或线性高斯塞德尔预调节器中的任一个或两者)更少的计算机处理迭代收敛。
图5是工作流程F的示意图，其概述了用于基于储层的异质岩层的三维上的压力分布的改进确定的储层模拟的根据本发明的序列的优选实施例。
在步骤200期间，通过读取储层和生产数据开始模拟。在步骤200期间读取的储层数据包括储层几何信息(其尺寸、x、y和z方向上的范围(长度))；诸如渗透率分布、孔隙率分布、层厚度的储层属性；相对渗透率数据；毛细管压力数据；诸如流体密度表、地层体积因子表、粘度表的流体属性数据；井的位置；和储层内的井穿孔的位置。生产数据包括针对在前一步骤中限定的井测量的或指定的油、水和气生产速率。生产数据还包括每个井的最小井底压力。在步骤200期间，储层模拟器也被初始化，将模拟日和时间步长设置为零。
在步骤202期间，根据方程(1)，在数据处理系统S中将储层数据和储层属性组织为由储层模型100的储层网格单元块的A矩阵表示的形式。在步骤202期间，同样，根据方程(1)，针对储层模型100的井的b向量，组织储层中的井的生产数据。
在步骤204期间，评估系数矩阵A以确定矩阵A中的任何一个或多个系数的系数值是否为负。如果不是这种情况，则处理进行到步骤206。
在步骤206期间，确定储层模型100的网格单元的最终压力分布。这可以例如根据申请人的在2018年4月26日提交的题为“Determining Pressure Distribution inHeterogeneous Rock Formations for Reservoir Simulation”在先共同未决美国专利申请No.15/963,251的技术来完成，该申请通过引用并入本文。
还应当理解，在步骤206期间可以使用许多其他迭代方法来确定储层中的压力分布。在Stüben,Klaus的“Solving Reservoir Simulation Equations,NinthInternational Forum on Reservoir Simulation,Abu Dhabi,United Arab Emirates,pg.1-53；December 2007中描述了当存在正定系数矩阵时可以用于确定储层中的压力分布的这种迭代方法的示例。
在步骤206期间，由迭代求解器进行的数据处理系统S中的处理继续，直到实现收敛。然后在步骤206期间，提供如此确定的压力作为储层模拟的一部分，以确定储层的网格单元中的流体流动。然后，在步骤206期间确定的多个维度上的储层单元中的压力分布可用于执行储层模拟步骤208。
在步骤208中，通过储层模拟器对由步骤206产生的储层模型100的单元的确定的压力分布和储层生产数据执行储层模拟，一系列时间步长被用于确定储层网格单元流体成分、压力、流速、相对存在和与在多个维度上组织的储层单元中的多相流体有关的相关地层参数，如方程(7)所示，以计算下一个流体流动储层模拟时间步长(n+1)的压力分布的估计。
在储层模拟步骤208期间，在步骤206期间确定的压力估计被用于数据处理系统S的储层模拟模块166中的流体流动储层模拟中。在步骤208期间的储层模拟被多次迭代地执行，以通过数值解确定压力，该压力将被用作一系列迭代时间步长的压力，如步骤210所示。储层模拟可以是若干已知类型中的一种，例如先前提到的那些。
在步骤208期间获得的储层模拟结果然后是在步骤212期间监测和调整一个或多个井的生产或性能、或生产和性能两者的基础，并因此是监测和调整储层的生产或性能、或生产和性能两者的基础。步骤212的储层工程活动包括例如以下储层生产管理操作中的一个或多个：调节来自储层中的生产井的流体的生产；调节或增加二次采收流体到所述储层的注入井中的注入；以及在储层中钻额外的井。
在步骤204期间，如果系数矩阵A中的一个或多个系数是负的，则在步骤220中根据本发明开始处理(图5和图6)。如步骤222中示意性地示出的，使得由步骤202产生的储层模型100的网格块单元的初始网格块压力可用于处理。
预调节系统
为了在数据处理系统S中求解方程(1)，即[A][p]＝[b]，必须进行预调节以在实际迭代次数内实现收敛。在不使用预调节器的情况下，为了收敛，解需要过多次数的迭代。这在计算机处理时间和计算机处理成本方面可能是昂贵的。计算上较便宜的预调节器之一是所谓的高斯塞德尔预调节器。但是存在已经讨论过的具有不定系数矩阵的高斯塞德尔预调节。
图6和图7示意性地示出了根据步骤220的处理，以获得每个时间步长的储层压力分布，从而确定储层模型的每个网格块(i,j)的压力p。在数据处理系统S中，对于根据方程(1)表示为[A][p]＝[b]的线性方程的矩阵组，执行每个网格块的压力的确定，其中A是已知系数矩阵，b是矩阵组的右端的已知大小的向量，并且p是未知网格块压力的向量。
为了在步骤220期间求解这种方程组，在数据处理系统S中构成预调节器矩阵。这通过将矩阵方程(1)乘以矩阵[B]的逆来完成，其中矩阵[B]表示作为矩阵[A]的近似的预调节矩阵。根据本发明，选择预调节矩阵[B]，使其具有逆矩阵，并且就计算机处理时间而言，当乘以矩阵[A]时进行计算是廉价的。根据本发明的预调节器系统允许对每个网格块的压力求解，如将描述的，这通过将Krylov向量加在一起来执行。
一旦在步骤222期间确定了储层模型矩阵的网格单元的初始近似压力p[r(x，y)]，则如在步骤224示意性地指示的，预调节方程(1)的矩阵方程关系的压力。
在根据步骤224的处理期间，作为初始步骤228，获得近似分析预调节器[B]，其中，数据处理系统S确定预调节器。根据本发明构成的预调节器被称为近似分析解预调节器。
矩阵[B]是A的近似，并且被广泛地称为方程(1)的“预调节矩阵”。矩阵[B]是这样的形式，当其逆乘以方程(1)的右端时，为储层中的未知压力分布提供了紧密近似。利用本发明，从近似分析解中求出矩阵[B]的元素。矩阵[B]表示系数矩阵[A]的相当好的近似，使得它容易可逆。将方程(1)从左边乘以逆[B-1]，得到方程(14)：
B-1Ap＝B-1b (14)
其可以缩写为方程(15)：
Figure BDA0002988214720000201
在根据步骤228的处理期间，由于在二维平面(x,y)中新的时间步长处储层模型中的网格单元的120(图2)所示位置处的源r*而产生的压力可以由下式表示
Figure BDA0002988214720000202
假设r*较小，则方程(2)可以写为
Figure BDA0002988214720000203
其中：
Figure BDA0002988214720000204
总压缩系数CT的测量是如下所示的孔隙(岩石)和流体压缩系数的线性组合：
CT＝Cp+SoCo+SwCw+SgCg (19)
其中在术语部分中解释了变量名称和单位。在步骤224期间确定由于单个源而产生的最大压力值(或最小值)pmax。
通过在数据处理系统S中在二维储层域(通常为x、y)上积分，根据方程(3)确定当前时间步长的平均储层压力
Figure BDA0002988214720000205
平均储层压力
Figure BDA0002988214720000206
还可以通过方程(3)的简单求和在网格块上近似，如下所示
Figure BDA0002988214720000211
其中距离点x*，y*的径向距离定义为
Figure BDA0002988214720000212
并且Vp是网格块的孔隙体积，由
Figure BDA0002988214720000213
定义。
使用体积平衡，平均储层压力
Figure BDA0002988214720000214
由下式给出：
Figure BDA0002988214720000215
为了简单起见，假设旧时间步长n处的初始储层压力为零。使根据方程(20)和(21)的压力值相等，未知的最大压力pmax可以根据方程(22)表示：
Figure BDA0002988214720000216
方程(16)结合方程(22)描述了储层中新时间步长处的近似压力分布。为了本发明的目的，根据方程(2)的时间步长可以被认为是对用于单个生产井的方程(1)的解的近似。该确定避免了需要使用迭代预调节器，当储层的传输率矩阵不定时，迭代预调节器不收敛。
多个井
如果储层包含多个生产井或注入井，则新时间步长处的压力分布可以通过叠加来近似，再次为了简化的目的，假设旧时间步长压力为零。假设储层中存在具有速率qi，i＝1，...，nw的nw个井，单独作用的每个井iw的最大压力值pmax(iw)由方程(5)和方程(6)确定。通过将为每个单个生产井确定的各个压力解叠加到复合解中，确定在从井iw测量的任何半径r处的近似压力值p[r(x，y)]。
Figure BDA0002988214720000217
其中r(x(iw)，y(iw))是井(iw)的坐标。
在作为步骤228的结果获得近似分析预调节器之后，执行步骤230以生成Krylov向量。在步骤230期间确定Krylov向量的系数，使得在步骤234期间确定的系数然后可以加在一起，以在步骤226中获得储层模型中的每个网格块(i,j)的压力。
在步骤230(图7)期间，数据处理系统S根据网格系数矩阵构成所谓的Krylov向量。在步骤230期间的向量的生成是在根据方程(9)生成的向量的Krylov子空间中。在步骤230期间，K表示根据方程(8)以及因此根据方程(1)生成的Krylov向量中的每一个的集合。Krylov向量格式中的集合K表示如下：
Figure BDA0002988214720000221
Figure BDA0002988214720000222
向量1，v1
Figure BDA0002988214720000223
这意味着
Bv1＝b
使用方程(23)，向量b包括井生产率q，以求解p。
第二Krylov向量，v2
Figure BDA0002988214720000224
这意昧着
Figure BDA0002988214720000225
使用方程(23)来求解v2，方程(26)的右端Av1可以作为方程(23)中的井生产率q被处理。
对于第i个向量，
Figure BDA0002988214720000226
使用方程(23)求解vi，Avi-1可以作为每个网格块的生产率q被处理。在方程(23)中，m是向量空间的维数。向量空间中的维数m远小于方程(1)中的未知量n的数目。还可以认为，m是用于构造方程(1)的解的向量的数量，超过该数量，v个向量变为线性相关的。
Krylov向量避免了执行矩阵乘法的需要。根据本发明在数据处理系统S中获得的Krylov向量为储层网格单元压力确定提供了优势。系数矩阵[A]是稀疏的，因为只有某些矩阵对角线是非零的，但是大多数矩阵具有零元素。由于矩阵的这些部分具有零值，因此不需要为矩阵的零元素分配计算机存储。
然而，为了求解方程(1)，矩阵[A]必须被求逆。然而，对于具有生产意义的储层的矩阵[A]的任何逆变成全矩阵，其中所有矩阵元素都是非零的。这种矩阵的求逆过程非常困难，计算强度大，并且成本高。计算和存储将是极高的，即，1百万×1百万网格单元全矩阵形式的储层矩阵需要计算机存储器中的1万亿个存储位置，这对于大多数计算机是不可能的。Krylov向量使得能够仅使用稀疏矩阵向量乘法，其结果各自产生N个未知量的向量，且因此需要少得多的计算机存储器存储。Krylov向量提供了矩阵方程(1)的解的表达式，其具有很少的计算机成本和存储要求。
在步骤232期间，由在步骤230期间构成的单独的Krylov向量vi组成的Krylov向量K被转换为标准正交向量矩阵，该标准正交向量矩阵由从单独的Krylov向量{v1，v2，...vm}的原始向量{v1，v2，...vm}中的对应的标准正交向量集合{u1，u2，...um}组成。
形成正交化和标准正交化以使向量彼此垂直，而标准正交化仅使向量长度一致(归一化)。为了表示将以正交向量的和来表示的未知向量(其为解向量)，执行正交化和标准正交化，以使向量彼此垂直(正交)。
在步骤232期间Krylov向量K的转换可以以几种常规方式完成，例如通过使用被称为Gram Schmidt或Arnoldi向量正交化。这种转换的一个示例在Lambers,Jim“Mat 415/515:Gram-Schmidt Orthogonalization Lecture 3 Notes"(2014)http://www.math.usm.edu/lambers/mat415/lecture3.pdf中描述。
未知系数ci的计算
步骤234包括确定未知系数ci，其用于确定储层单元模型的网格块之间的压力分布。为了步骤234的目的，根据方程(1)的储层网格块中的数字化形式的压力分布的解以下面的形式表示：
Figure BDA0002988214720000231
应当注意，方程(28)中的标准正交向量ui实际上是方程(1)的近似解。在向量u的值是理想解的情况下，则在方程(28)中，系数的数量仅为一，m＝1，并且仅需要一个向量来表示解。因此，重要的是获得u个向量的值的解的相当准确的估计。接下来，方程(9)的预调节方程组的余数被定义为：
r＝b-Ap (29)
将方程(28)代入方程(15)并定义下列方程(30)中的二次函数j
Figure BDA0002988214720000241
关于未知系数c1，c2，....cm的二次函数j的最小值出现在j关于每个系数的一阶导数为零处。执行这些操作得到未知系数c1，c2，....cm的线性方程组，如下：
Figure BDA0002988214720000242
方程(31)中的线性方程组通常应该具有更小的尺寸(从Krylov方程组生成的线性独立向量的数量通常远小于n)，因此可以通过利用直接方法高斯消元法处理来求解，其示例在Wikipedia,“Gaussian Elimination”pgs.1-8；https://en.wikipedia.org/wiki/Gaussian_elimination描述
在步骤234期间确定的压力估计在步骤226期间用于储层模拟，其方式已经描述为确定用于在步骤208期间监测和调整井和储层的生产和性能的储层网格单元流体成分、压力、流速、与储层单元中的多相流体相关的相对存在和相关地层参数。
示例
结果
本发明的第一示例是针对小问题来说明不定矩阵的情况，其中可以容易地确定系数矩阵A的特征值。
示例1——二维异质储层，5×5网格
图8是简化的储层模型300的示意图，用于获得根据本发明的压力分布与现有技术方法的结果相比较的比较结果。示例性储层300是油储层，其在可略微压缩的流体的泡点压力之上操作。它具有3000ft×3000ft的正方形面积，没有流动边界条件，并且可以有多个井注入或生产。对于该示例，假设存在位于中心的单个注入井，如302所示。因此，面积网格大小为Δx＝Δy＝600ft，其中Δz＝30ft厚度，并且其中单元数目Nx＝5，Ny＝5。
储层和流体数据
针对示例模拟，初始储层压力为3000psi，单相流体粘度＝1cP。假设总压缩系数对于整个储层(除了注入井所处的中心处的位置302)为3E-61/psi。假设中心井处的总压缩系数CT为-300E-61/psi。随机渗透率和孔隙率分布在储层模型300的X和Y维度上被呈现为如下：
随机渗透率和孔隙率分布
Min Kx＝0.2677mD；max Kx＝49.8mD；average Kx＝24.6mD
Min Ky＝5.01mD，max Ky＝49.95mD，average Ky＝30.75mD
孔隙率分布
Figure BDA0002988214720000251
要生成的Krylov向量的最大数目＝5
系数矩阵A的特征值：
-0.388E+01 -0.269E+01 -0.148E+01 0.101E+01 -0.562E+01
0.485E+01 -0.164E+01 -0.851E+00 -0.635E+00 -0.184E+01
-0.198E+01 -0.304E+00 -0.335E+01 -0.320E+01 -0.404E+01
0.295E+01 -0.953E+00 -0.157E+00 -0.287E+01 -0.142E+01
0.461E+01 -0.633E+01 -0.422E+01 -0.219E+01 -0.367E+00
如所看到的，除了具有正号的第四特征值(1.01)之外，系数矩阵A的所有特征值具有负号。因此系数矩阵A是不定的。
模拟
使用三种预调节方法用于比较目的：执行点高斯塞德尔预调节；线高斯塞德尔预调节；和根据本发明的方法，将储层模拟器运行一个时间步长，以确定该示例在30天的时间的压力分布。
在下面的表1、表2和表3中阐述了对于如通过高斯消元法确定的图8的简化模型确定的压力分布，以及在根据本发明的选定次数的迭代之后的压力分布：
表1——直接法(高斯消元法)的精确解
精确井压力＝2702.05psi
时间＝30天的压力分布。
0.27768E+04 0.27675E+04 0.27540E+04 0.27599E+04 0.27657E+04
0.27812E+04 0.27695E+04 0.27224E+04 0.27526E+04 0.27597E+04
0.27805E+04 0.27651E+04 0.27021E+04 0.27376E+04 0.27540E+04
0.27796E+04 0.27693E+04 0.27292E+04 0.27438E+04 0.27542E+04
0.27784E+04 0.27695E+04 0.27413E+04 0.27515E+04 0.27552E+04
表2——本发明在迭代10次的解
井压力＝2708.78psi
0.27643E+04 0.27578E+04 0.2748EE+04 0.27551E+04 0.27612E+04
0.27671E+04 0.27595E+04 0.27245E+04 0.27516E+04 0.27579E+04
0.27656E+04 0.27562E+04 0.27088E+04 0.27400E+04 0.27549E+04
0.27642E+04 0.27593E+04 0.27320E+04 0.27456E+04 0.27553E+04
0.27645E+04 0.27592E+04 0.27411E+04 0.27523E+04 0.27561E+04
表3——本发明在迭代＝120次时的解
井压力＝2702.84psi
0.27753E+04 0.27663E+04 0.27534E+04 0.27593E+04 0.27653E+04
0.27795E+04 0.27683E+04 0.27226E+04 0.27525E+04 0.27596E+04
0.27787E+04 0.27639E+04 0.27028E+04 0.27379E+04 0.27542E+04
0.27776E+04 0.27681E+04 0.27296E+04 0.27441E+04 0.27544E+04
0.27767E+04 0.27682E+04 0.27413E+04 0.27516E+04 0.27553E+04
图9是针对储层模型300根据本发明的压力分布确定中的收敛误差与根据现有技术的处理结果的比较结果的显示。点高斯塞德尔(PGS)的收敛行为在304处指示，线高斯塞德尔(LGS)方法的收敛行为在306处指示。根据本发明的收敛行为在图9中的308处示出。
如图9所示，纵轴上示出了以桶每天(b/d)为单位的最大残差，横轴示出了迭代次数。如所示，点高斯塞德尔方法和线高斯塞德尔方法两者发散，而本发明收敛。
例如，在10次迭代结束时，对于点高斯塞德尔方法，最大残差(误差)是1571846b/d，对于线高斯塞德尔是6657102829b/d，对于本发明仅是2.5b/d。本发明提供了一种正确解的收敛解，而相反，点高斯塞德尔方法和线高斯塞德尔方法不提供。此外，本发明使用稀疏矩阵而不是全矩阵。因此，本发明提供了一种有效的备选方案，其针对求解具有数亿个未知量的矩阵是收敛的。这节省了储层矩阵处理解的计算机处理和操作时间，并且不需要继续执行不收敛而是发散的重复迭代。此外，引入进一步复杂的计算机处理的时间步长切割的常规方法不是必需的，并且不需要被执行。这是根据本发明的计算机操作的进一步改进。
示例2——100×100网格
图10是示例2的简化的示例储层模型400的示意图。井位于402，在模型400的中心(50，50)附近。位置402处的井用于示例2的目的，注入1000b/d。模型400中的储层长度在x、y每个方向上为60000ft，产生与示例1相同的网格大小。单元(50，50)和(51，50)的总压缩系数CT为-300E-6，1/psi，并且模型400的剩余单元的总压缩系数为3E-61/psi总压缩系数。示例2的储层模型400被配置为具有与示例1中相同的垂直厚度和统计性质的异质储层。
图11是根据本发明确定的储层模型400的压力分布确定中的收敛误差与根据现有技术的点高斯塞德尔或线性高斯塞德尔预调节器获得的处理结果的比较结果的显示。
如图11中的404和406所示，点高斯塞德尔方法和线高斯塞德尔方法都有所发散。如408所示，本发明能够将残差减少到网格块孔隙体积的1.3％，其对应于65b/d的误差。根据本发明的迭代的进一步的迭代表现出稳定，但是残差减小的少。该示例进一步验证了本发明可以找到正确的解，而其他迭代方法不收敛，从而需要时间步长切割以尝试完成储层模拟运行。
示例3——百万个单元异质问题(1000×1000)
图12是示例3的简化的示例储层模型500的示意图。井502如所示位于模型500的中心，并且四个额外的井504、506、508和510位于模型500的每个象限的中心。用于该示例的目的，这五个井中的每一个都是以1000b/d注入。如所示针对网格单元建立负压缩系数网格块，井502的西、东、北、南的一个块以及模型500的中心处的井块。其余的参数、垂直厚度和统计性质与示例2和示例3中的相同。
图13是模型500的压力分布确定中的收敛误差的比较结果的显示。如在530和532可见，点高斯塞德尔方法和线高斯塞德尔方法两者都有所发散。本发明示出了如534所示的收敛行为。初始残差(最大范数)1000b/d在10次迭代中已经减小到707.8b/d，其对应于12.8％的网格块孔隙体积。可以看到进一步的迭代减少残差，并且根据本发明的处理收敛。
再次，如图13中的530和532所示，点高斯塞德尔方法和线高斯塞德尔方法都有所发散。该附加示例进一步验证了本发明可以找到正确的解，而其他迭代方法不收敛，需要时间步长切割以尝试完成储层模拟运行。
在图9、图11和图13的每一个中，可以看出传统的点高斯塞德尔和线高斯塞德尔迭代预调节器方法都有所发散。具体地，压力计算中的误差随着迭代次数的增加而增加。因此，这些传统的预调节器方法不能产生正确的压力分布。另一方面，可以看出，在本发明如图9、图11和图13所示的示例1、示例2和示例3的每一个中，压力确定中的误差由于连续迭代而减小。因此，本发明提供了正确的压力解，并且随着迭代次数的增加而收敛。
因此可以看出，本发明改进了在储层模拟操作期间计算机的功能。本发明提供了在针对储层模拟确定储层模型中的压力分布期间在确定存在不定系数矩阵的情况下收敛的计算机能力。本发明避免了在这种时间步长切割情况下的传统现有技术补救措施。
本发明已经被充分地描述，使得在储层建模和模拟领域具有一般知识的人员可以再现和获得本发明在此提到的结果。尽管如此，作为本发明主题的技术领域的任何技术人员可以进行本文请求中未描述的修改，以将这些修改应用于确定的结构和方法；或者在其使用和实践中，需要所附权利要求中要求保护的主题；这些结构和方法应被涵盖在本发明的范围内。
应当注意和理解，在不背离所附权利要求中阐述的本发明的精神或范围的情况下，可以对上面详细描述的本发明进行改进和修改。

Claims (17)

1.一种对生产储层中的流体流动进行储层模拟的计算机实现方法，所述计算机实现方法具有确定所述储层内的压力分布的改进的收敛，以针对从所述储层预计生产的时间步长序列，确定所述储层中的模拟流体流动，所述储层被组织成储层单元的三维网格，针对所述时间步长序列，基于实际流体压力和来自所述储层中的井的流体生产来执行所述储层模拟，所述计算机实现方法确定从储层单元的所述三维网格到所述井的模拟流体流动和所述井的模拟井生产速率，所述计算机实现方法包括以下计算机处理步骤：
(a)对于所述时间步长序列中的时间步长，构成已知储层属性的初始计算机矩阵和实际储层流体生产测量的初始计算机向量；
(b)确定已知储层属性的所述初始计算机矩阵是否对角占优；
(c)如果是对角占优的，则基于所述储层的所述网格单元中的各个网格单元中的压力分布的所构成的初始测量，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的所述各个网格单元中的压力分布；
(d)基于所述储层的所述网格单元中的所确定的估计压力分布，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的估计流体流速；
(e)确定所述储层的所述网格单元中的所确定的流体流速是否已收敛；以及
(f)如果已收敛，则结束针对所述时间步长的所述储层中的流体流动的所述储层模拟；以及
(g)如果已知储层属性的所述初始计算机矩阵不是对角占优的，则生成近似分析预调节器；
(h)将所生成的近似分析预调节器应用到已知储层属性的所述初始计算机矩阵，以构成所述储层的所述网格单元中的压力分布的测量；
(i)根据所述储层的所述网格单元中的压力分布构成的测量生成Krylov向量；
(j)基于根据压力分布构成的测量所生成的Krylov向量，确定储层压力分布的系数向量；
(k)基于所确定的储层压力分布的系数向量，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的估计流体流速；
(l)确定基于所确定的储层压力分布的系数向量的所述储层的所述网格单元中的所述估计流体流速是否已收敛；以及
(m)如果已收敛，则结束针对所述时间步长的所述储层中的流体流动的所述储层模拟。
2.根据权利要求1所述的方法，还包括基于所述储层模拟来调整所述井中的至少一个的生产的步骤。
3.根据权利要求1所述的方法，还包括基于所述储层模拟来调整所述井中的至少一个的性能的步骤。
4.根据权利要求1所述的方法，还包括当所述确定的流体流速已经收敛时，递增到结束所述储层模拟的随后的时间步长的计算机处理步骤。
5.根据权利要求4所述的方法，还包括对所述随后的时间步长重复计算机处理步骤(a)到(m)的步骤。
6.根据权利要求1所述的方法，还包括将所生成的Krylov向量转换为正交Krylov向量的计算机处理步骤。
7.根据权利要求6所述的方法，还包括对所述正交Krylov向量执行标准正交化的计算机处理步骤。
8.根据权利要求7所述的方法，其中，确定储层压力分布的系数向量的计算机处理步骤包括：在对所述正交Krylov向量执行标准正交化的步骤之后，确定储层压力分布的系数向量的计算机处理步骤。
9.根据权利要求1所述的方法，其中，确定储层压力分布的系数向量的计算机处理步骤包括：通过执行系数向量的线性方程组的高斯消元法来处理所述线性方程组的计算机处理步骤。
10.一种数据处理系统，其执行生产储层中的流体流动的储层模拟，所述数据处理系统具有确定所述储层内的压力分布的改进的收敛，以针对从所述储层预计生产的时间步长序列，确定所述储层中的模拟流体流动，所述储层被组织成储层单元的三维网格，针对所述时间步长序列，基于实际流体压力和来自所述储层中的井的流体生产来执行所述储层模拟，计算机实现方法确定从储层单元的所述三维网格到所述井的模拟流体流动和所述井的模拟井生产速率，所述数据处理系统包括：
(a)处理器，其执行以下计算机处理步骤：
(1)对于所述时间步长序列中的时间步长，构成已知储层属性的初始计算机矩阵和实际储层流体生产测量的初始计算机向量；
(2)确定已知储层属性的所述初始计算机矩阵是否对角占优；
(3)如果是对角占优的，则基于所述储层的所述网格单元中的各个网格单元中的压力分布的所构成的初始测量，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的所述各个网格单元中的压力分布；
(4)基于所述储层的所述网格单元中的所确定的估计压力分布，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的估计流体流速；
(5)确定所述储层的所述网格单元中的所确定的流体流速是否已收敛；以及
(6)如果已收敛，则结束针对所述时间步长的所述储层中的流体流动的所述储层模拟；以及
(7)如果已知储层属性的所述初始计算机矩阵不是对角占优的，则生成近似分析预调节器；
(8)将所生成的近似分析预调节器应用到已知储层属性的所述初始计算机矩阵，以构成所述储层的所述网格单元中的压力分布的测量；
(9)根据所述储层的所述网格单元中的压力分布构成的测量生成Krylov向量；
(10)基于根据压力分布构成的测量所生成的Krylov向量，确定储层压力分布的系数向量；
(11)基于所确定的储层压力分布的系数向量，针对所述储层模拟的所述时间步长，确定所述储层的所述网格单元中的估计流体流速；
(12)确定基于所确定的储层压力分布的系数向量的所述储层的所述网格单元中的所述估计流体流速是否已收敛；以及
(13)如果已收敛，则结束针对所述时间步的所述储层中的流体流动的所述储层模拟；以及
(b)工作站，用于提供所述储层模拟结果的显示。
11.根据权利要求10所述的数据处理系统，还包括：
存储器，其存储所述储层模拟结果。
12.根据权利要求10所述的数据处理系统，还包括当所确定的流体流速已经收敛时，所述处理器递增到结束所述储层模拟的随后的时间步长，
13.根据权利要求12所述的数据处理系统，还包括所述处理器针对所述随后的时间步长重复计算机处理步骤(a)至(m)。
14.根据权利要求10所述的数据处理系统，还包括所述处理器将所生成的Krylov向量转换为正交Krylov向量。
15.根据权利要求14所述的数据处理系统，还包括所述处理器对所述正交Krylov向量执行标准正交化。
16.根据权利要求15所述的数据处理系统，其中，所述处理器在确定储层压力分布的系数向量时，在对所述正交Krylov向量执行标准正交化之后，确定储层压力分布的系数向量。
17.根据权利要求15所述的数据处理系统，其中，所述处理器在确定储层压力分布的系数向量时，通过系数向量的线性方程组的高斯消元法来处理所述线性方程组。
CN201980062351.6A 2018-09-24 2019-09-12 具有用于非对角占优的不定系数矩阵的压力求解器的储层模拟 Active CN112752894B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US16/139,326 2018-09-24
US16/139,326 US11461514B2 (en) 2018-09-24 2018-09-24 Reservoir simulation with pressure solver for non-diagonally dominant indefinite coefficient matrices
PCT/US2019/050875 WO2020068442A1 (en) 2018-09-24 2019-09-12 Reservoir simulation with pressure solver for non-diagonally dominant indefinite coefficient matrices

Publications (2)

Publication Number Publication Date
CN112752894A true CN112752894A (zh) 2021-05-04
CN112752894B CN112752894B (zh) 2023-02-28

Family

ID=68069882

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201980062351.6A Active CN112752894B (zh) 2018-09-24 2019-09-12 具有用于非对角占优的不定系数矩阵的压力求解器的储层模拟

Country Status (6)

Country Link
US (1) US11461514B2 (zh)
EP (1) EP3857022B1 (zh)
JP (1) JP6967176B1 (zh)
CN (1) CN112752894B (zh)
CA (1) CA3112183C (zh)
WO (1) WO2020068442A1 (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11461514B2 (en) * 2018-09-24 2022-10-04 Saudi Arabian Oil Company Reservoir simulation with pressure solver for non-diagonally dominant indefinite coefficient matrices
CN113012276B (zh) * 2021-01-27 2021-09-24 中国科学院空天信息创新研究院 基于辐射度的地表高分辨率光谱信息遥感反演方法
CN113049471B (zh) * 2021-03-23 2021-10-08 中国石油大学(北京) 一种碳酸盐岩层序地层的孔隙度演化过程的恢复方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102165413A (zh) * 2008-09-30 2011-08-24 埃克森美孚上游研究公司 自适应迭代求解器
US20120221303A1 (en) * 2011-02-24 2012-08-30 Schlumberger Technology Corporation System and Method For Performing Reservoir Simulation Using Preconditioning
US9063882B1 (en) * 2010-09-09 2015-06-23 Sas Ip, Inc. Matrix preconditioners for simulations of physical fields
EP2960430A2 (en) * 2014-06-27 2015-12-30 Services Petroliers Schlumberger Multilevel monotone constrained pressure residual multiscale techniques
CN106326591A (zh) * 2016-08-31 2017-01-11 西南石油大学 水力压裂过程中裂缝内压裂液的压力场获取方法及装置

Family Cites Families (31)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1759226A1 (en) * 2004-06-07 2007-03-07 ExxonMobil Upstream Research Company Method for solving implicit reservoir simulation matrix equation
US20080167849A1 (en) * 2004-06-07 2008-07-10 Brigham Young University Reservoir Simulation
US7774725B1 (en) 2005-11-04 2010-08-10 Purdue Research Foundation Computationally efficient modeling and simulation of large scale systems
BR112012010094A2 (pt) 2009-10-28 2016-05-31 Chevron Usa Inc método de volume finito em multiescala para uso na simulação de um modelo geológico de escala fina de um reservatório de subsuperfície, sistema para uso em simulação de um modelo geológico de escala fina de um reservatório de subsuperfície, e, software
CA2785569A1 (en) * 2010-02-02 2011-08-11 Hector Klie Multilevel percolation aggregation solver for petroleum reservoir simulations
US8682628B2 (en) * 2010-06-24 2014-03-25 Schlumberger Technology Corporation Multiphase flow in a wellbore and connected hydraulic fracture
US9540911B2 (en) * 2010-06-24 2017-01-10 Schlumberger Technology Corporation Control of multiple tubing string well systems
BR112012031776A2 (pt) 2010-07-26 2016-11-01 Exxonmobil Upstream Res Co método e sistema para simulação multinível paralela
US8924186B1 (en) * 2010-09-09 2014-12-30 Sas Ip, Inc. Simulations of physical systems for multiple excitations
US8583411B2 (en) * 2011-01-10 2013-11-12 Saudi Arabian Oil Company Scalable simulation of multiphase flow in a fractured subterranean reservoir as multiple interacting continua
US8437999B2 (en) * 2011-02-08 2013-05-07 Saudi Arabian Oil Company Seismic-scale reservoir simulation of giant subsurface reservoirs using GPU-accelerated linear equation systems
US9164191B2 (en) * 2011-02-09 2015-10-20 Saudi Arabian Oil Company Sequential fully implicit well model for reservoir simulation
US20130085730A1 (en) * 2011-10-04 2013-04-04 Gareth J. Shaw Preconditioner for reservoir simulation
EP2783069B1 (en) * 2011-11-22 2015-09-09 Saudi Arabian Oil Company Coupled pipe network - reservoir modeling for multi-branch oil wells
WO2013119906A1 (en) * 2012-02-09 2013-08-15 Saudi Arabian Oil Company Multi-level solution of large-scale linear systems in simulation of porous media in giant reservoirs
US9208268B2 (en) 2012-02-14 2015-12-08 Saudi Arabian Oil Company Giga-cell linear solver method and apparatus for massive parallel reservoir simulation
EP2845143A4 (en) * 2012-05-30 2016-09-28 Landmark Graphics Corp OIL OR GAS PRODUCTION USING COMPUTER SIMULATION OF OIL OR GAS FIELDS AND PRODUCTION PLANTS
US20140025720A1 (en) 2012-07-17 2014-01-23 Nvidia Corporation Solving linear equation systems with multiple right hand sides by krylov subspace expansion
US9798698B2 (en) * 2012-08-13 2017-10-24 Nvidia Corporation System and method for multi-color dilu preconditioner
US20150338550A1 (en) * 2012-11-20 2015-11-26 Stochastic Simulation Limited Method and system for characterising subsurface reservoirs
US20140207426A1 (en) * 2013-01-18 2014-07-24 The Board Of Trustees Of The Leland Stanford Junior University Simulation of phenomena characterized by partial differential equations
US20160177674A1 (en) * 2013-08-27 2016-06-23 Halliburton Energy Services, Inc. Simulating Fluid Leak-Off and Flow-Back in a Fractured Subterranean Region
US20150186562A1 (en) * 2013-12-30 2015-07-02 Halliburton Energy Services, Inc Preconditioning a Global Model of a Subterranean Region
US20160202389A1 (en) * 2015-01-12 2016-07-14 Schlumberger Technology Corporation H-matrix preconditioner
US10242136B2 (en) * 2015-05-20 2019-03-26 Saudi Arabian Oil Company Parallel solution for fully-coupled fully-implicit wellbore modeling in reservoir simulation
US10838093B2 (en) 2015-07-02 2020-11-17 Exxonmobil Upstream Research Company Krylov-space-based quasi-newton preconditioner for full-wavefield inversion
US10209685B2 (en) * 2015-10-29 2019-02-19 Mitsubishi Electric Research Laboratories, Inc. Method and apparatus for preconditioned model predictive control
WO2018015972A1 (en) * 2016-07-18 2018-01-25 Indian Institute Of Science Eigen augmentation methods for electromagnetic modelling and simulation
WO2018089059A1 (en) * 2016-11-08 2018-05-17 Landmark Graphics Corporation Selective diffusion inclusion for a reservoir simulation for hydrocarbon recovery
US11499412B2 (en) * 2017-11-24 2022-11-15 Total Se Method and device for determining hydrocarbon production for a reservoir
US11461514B2 (en) * 2018-09-24 2022-10-04 Saudi Arabian Oil Company Reservoir simulation with pressure solver for non-diagonally dominant indefinite coefficient matrices

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102165413A (zh) * 2008-09-30 2011-08-24 埃克森美孚上游研究公司 自适应迭代求解器
US9063882B1 (en) * 2010-09-09 2015-06-23 Sas Ip, Inc. Matrix preconditioners for simulations of physical fields
US20120221303A1 (en) * 2011-02-24 2012-08-30 Schlumberger Technology Corporation System and Method For Performing Reservoir Simulation Using Preconditioning
EP2960430A2 (en) * 2014-06-27 2015-12-30 Services Petroliers Schlumberger Multilevel monotone constrained pressure residual multiscale techniques
CN106326591A (zh) * 2016-08-31 2017-01-11 西南石油大学 水力压裂过程中裂缝内压裂液的压力场获取方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
康忠良等: "适用于混合网格的并行GMRES+LU-SGS方法", 《空气动力学学报》 *
张理涛: "线性方程组与鞍点问题的迭代法与预处理技术研究", 《中国优秀博硕士学位论文全文数据库（博士） 基础科学辑》 *

Also Published As

Publication number Publication date
JP2021534517A (ja) 2021-12-09
JP6967176B1 (ja) 2021-11-17
CA3112183C (en) 2022-03-15
EP3857022B1 (en) 2023-07-12
US20200097622A1 (en) 2020-03-26
US11461514B2 (en) 2022-10-04
CN112752894B (zh) 2023-02-28
CA3112183A1 (en) 2020-04-02
WO2020068442A1 (en) 2020-04-02
EP3857022A1 (en) 2021-08-04

Similar Documents

Publication Publication Date Title
US11066907B2 (en) Sequential fully implicit well model with tridiagonal matrix structure for reservoir simulation
US10762258B2 (en) Parallel solution for fully-coupled fully-implicit wellbore modeling in reservoir simulation
EP2783069B1 (en) Coupled pipe network - reservoir modeling for multi-branch oil wells
CN112752894B (zh) 具有用于非对角占优的不定系数矩阵的压力求解器的储层模拟
US10175386B2 (en) Sequential fully implicit well model with tridiagonal matrix structure for reservoir simulation
US10822925B2 (en) Determining pressure distribution in heterogeneous rock formations for reservoir simulation
Yoon et al. Hyper-reduced-order models for subsurface flow simulation
US10113400B2 (en) Sequential fully implicit well model with tridiagonal matrix structure for reservoir simulation
EP3423672A1 (en) Sequential fully implicit well model with tridiagonal matrix structure for reservoir simulation
Batycky et al. Reservoir simulation
de Menezes Development of a multipurpose reservoir simulator based on a plugin architecture
Temizel et al. An Analysis of Multiplicative Schwarz Procedure in Coupling of Reservoir and Networks in Next-Generation Reservoir Simulators
Angelim et al. NUMERICAL FLOW SIMULATION IN NATURALLY FRACTURED PETROLEUM RESERVOIRS USING A MPFA-STREAMLINE FORMULATION
EP3423676A1 (en) Sequential fully implicit well model with tridiagonal matrix structure for reservoir simulation
Tambue EXPONENTIAL EULER TIME INTEGRATOR FOR ISOTHERMAL INCOMPRESSIBLE TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant