专利名称：：用来模拟物理系统的特性的方法
技术领域：
：本发明一般涉及模拟物理系统的至少一个特性。在一个方面，本发明涉及一种用来模拟诸如含油气储层之类的物理系统以预测在储层中流体性质和特性的方法。数值模拟作为一种通过使用计算机模拟物理系统的方法广泛用在工业领域。在大多数情况下，希望把发生在物理系统中的运输过程模型化。运输的东西典型的是质量、能量、动量、或其某种组合。通过使用数值模拟，有可能复制和观察物理现象及确定设计参数，而不用使用模型和设备的实际试验。因此能期望，大大地减小设计时间和成本。有较大兴趣的一种模拟类型是一个从真实含油气储层的模型的性能推出该储层的特性的过程。储层模拟的目的在于足够好地理解发生在储层中的复杂化学、物理及流体流动过程，以预测储层的未来特性使油气采收率最大。储层模拟常常指的是在储层内流动的流体动力学，但在较大意义上，储层模拟也能指包括储层、注入井、采油井、表面出油管线、及表面处理设施的整个石油系统。数值模拟的原理在于通过计算机数值解描述一种物理现象的方程。这样的方程一般是普通的微分方程和偏微分方程。作为一种用来数值解这样的方程的方法，已知有有限元法、有限差分法、有限体积法等。在这些方法的每一种中，把要模型化的物理系统划分成较小单元(其一组称作栅格或网格)，并且把在每个单元中连续变化的状态变量由用于每个单元的值的组表示。一个原始微分方程由一组表示在每个较小单位或单元内的质量、能量和/或动量守恒及在单元之间的质量、能量、和/或动量的运动的基本原理的方程替换。这些方程能共有几百万个。对于每个单元由有限数量个值对连续变化值的这种替换称作“离散化”。为了分析按时间变化的现象，必须计算在称作时间步长的分立时间间隔处的物理量，而不顾作为时间函数的连续变化条件。运输过程的时间依赖模型化按时间步长顺序行进。对于大多数运输过程，基本方程是非线性的，因为在一个单元中质量、能量或动量的量和在单元之间质量、能量及动量的运动一般与定义单元物理状态的变量具有非线性关系。在模拟油气储层时，例如模型化储层的方程是非线性偏微分方程，这些方程描述整个储层所有流体的不稳定状态流动，并且在整个储层中与流体随时间的压力和饱和变化有关。为了模拟多种物理系统，希望使用隐函数法，其中在单元之间运输实体的运动取决于在时间步长末端处的解。隐函数法要求在时间步长末端处的未知量全部一起确定。结果，如果方程是非线性的，则一般使用迭代计算未知量。迭代涉及从对未知量的某种初始推测开始并且应用某种重复计算以改进推测，直到在足够数量的迭代之后，使方程满足在某个可接收容许值内。由于每次迭代需要计算时间，所以鼓励使用尽可能减少计算时间的迭代法。已经提出了用来解非线性方程组的多种迭代法。一个例子是熟知的牛顿-莱甫森法。在牛顿-莱甫森迭代中使用的近似导致使在每个单元处的未知量与在其附近处的未知量相关的线性方程组。这些方程组组成一个全局矩阵方程，然后解该全局矩阵方程以得到解的下次估计。如果物理系统的表示是线性的，则得到一个类似矩阵公式。在任一种情况下，矩阵方程一般十分大，并且最好迭代地解。用来解这种矩阵方程的迭代法是一个称作点高斯-塞得尔法的过程。在点高斯-塞得尔法中，一个单元一个单元地计算一个新的解估计。在每个单元处，通过对于该单元解质量、能量、及动量平衡方程得到新的估计，同时把在邻近单元处的未知对应值保持固定在其最新估计处。在该过程中，相邻单元是一个当前单元处于与其连通的单元。一个单元的质量、能量、或动量平衡方程将包含乘以在其邻近处的未知量的项。当对于系统中的所有方程重复这种计算时，创建一个新答案阵列。然后检查该阵列，以确定诸值是否满足单元方程。为了这样作，便利地是对于每个方程定义一个余数(r)。如果新值满足公式，那么所有余数将是零或非常小。如果不是，借助于基于以前迭代的未知量的更新值重复该过程。重复该过程，直到所有余数可接收地接近零。这种类型的迭代法称作“点”迭代法，因为一次一个点或一个单元地进行该方法。如果点高斯-塞得尔法由点逐次超松驰法或PSOR替换，则能得到较快收敛。在PSOR中，在每次迭代处估计解的变化乘以一个超松驰参数ω，该超松驰参数ω必须具有一个在一与二之间的值。PSOR在模拟中的成功应用一般限于较简单的模型。因为PSOR法是其中一次仅计算一个单元的未知值的“显式”方法，所以PSOR法易于减慢收敛。这种缺陷已经导致在模拟方法包括更隐含的努力。用来这样作的一种方法称作线逐次超松驰法(LSOR)。LSOR通过在一个方向上保留隐含对PSOR改进。其中，在把相邻列或行的影响保持在其最新估计的同时，同时解对于列或行的质量、能量、或动量平衡方程。LSOR应用的例子能在(1)Mattax，C.C.和Dalton，R.L.的ReservoirSimulation(储层模拟)，MonographVolume13，SocietyofPetroleumEngineers(石油工程师协会)(1990年)和(2)Aziz，K.和Settari，A.的PetroleumReservoirSimulation(石油储层模拟)，AppliedSciencePublishersLtd，London，1979年中找到。在过去使用的LSOR法主要应用于其中以具有良好定义的行或列正规、构造栅格组织单元的模型。借助于以缺乏这种正规结构的栅格布置的至少一些单元已经提出了多种模型。相信本发明的实施表示LSOR原理对未构造栅格的最初应用。对于未构造栅格与结构栅格相比由解技术的高成本已经减慢了未构造栅的工业使用。对于使用所有类型的单元配置能用来分析物理系统的表示的模拟方法的需要是存在的。本发明的方法用来模拟一种物理系统的至少一个特征，而不顾物理系统是否已经离散成出现在构造或未构造栅格或两者组合中的单元。该方法的第一步骤是，把物理系统离散成彼此相邻布置的多个体积单元，以便在每对相邻单元之间有一条边界。对于每个单元，构造线性方程，这些线性方程把单元的状态变量与其相邻单元的状态变量相联系。以后的步骤是把一个可运输性值与每条边界相联系，并且然后对应于可运输性值的下降顺序来排序边界。边界排序然后用来拓扑地构造一维单元串。对于每个串通过组合与每串的单元有关的线性方程，产生一个矩阵方程。然后通过一次一串地解用于每串的矩阵方程得到单元状态变量的改进估计，直到已经解完用于所有串的矩阵方程。迭代地重复该过程，直到满足收敛条件。这种解产生对于所有单元的状态变量，这些状态变量同时满足对于所有单元的线性方程。由迭代产生的状态变量能用来模拟物理系统的至少一个特性。在一个最佳实施例中，串的建造使用一种提高在串单元中的单元边界处具有高可运输性的串形成的规则。通过参照如下详细描述和附图将更好地理解本发明和其优点，其中类似元件已经给出类似号码并且其中图1是具有五行和十列的两维笛卡儿栅格系统的简化例子，其中单元的几何形状指示在单元之间的耦合强度，在单元之间的耦合对于垂直方向(在列内)的流量是最强的。图2是类似于图1例子的两维笛卡儿栅格系统的简化例子，不同之处在于，对于在水平方向的流动在单元之间的耦合在从左至右的运动中减小，而对于在垂直方向的流动在从左至右的运动中增大。图3是一个两维未构造栅格系统的简化例子，其中单元不具有相同形状，并且在单元之间的耦合不跟随固定图案。图4描绘一个简单两维3单元乘5单元的15个单元栅格，表示在单元之间的可运输性排序。图5描绘图4的栅格，表示把15单元栅格分解成一串单元的初始步骤。图6描绘在图5中所示的串已经切开以形成两个串之后图4的栅格。图7描绘一个简单两维3单元乘6单元的18个单元栅格，表示在单元之间的可运输性排序。图8描绘图7的栅格，表示把18单元栅格分解成一串单元的初始步骤。图9描绘在图8中所示的串已经切开以形成四个串之后图7的栅格。附图不打算从本发明的范围排除是这些具体实施例的通常和期望改进的结果的其他实施例。本发明提供一种用来模拟用偏微分方程数值表示的物理系统的新方法。该方法能在模拟离散成构造栅格、未构造栅格、或两者组合的两维和三维域时使用。它也能用在这样的情形中，其中计算方法产生具有大于三维的拓扑，如在模拟分裂多孔介质时出现的那样。本发明在模拟其中运输现象正在发生的物理系统的一个特性时特别有用。在该专利中所使用的术语“运输现象”在广义上使用，以包括动量运输(粘性流动)、能量运输(热传导、对流、及辐射)、及质量运输(扩散)。本发明能应用于大不相同的领域，如物理、岩石特性化、晶体学、电气工程、生物学、数学、流体力学、及石油工程。在模拟操作中的普通实践是，表示在由方程Mx＝y模拟的物理域上由基本偏微分方程的离散化产生的线性方程组(其中M是一个尺寸n×n的系数矩阵，即n行乘n列，x是表示未知值的尺寸n的一个列向量，y是表示一组已知值的尺寸n的一个列向量)。在模拟操作中的基本操作是解这种线性方程系统。这种操作例如在用于非线性方程解的牛顿-莱甫森法中、以及在普通差分方程的隐式积分期间产生。用来解偏微分方程的传统方法取决于系数矩阵M的块分区。这些解法包括诸如基于直线的松驰之类的迭代技术、诸如附加校正之类的收敛加速方案、及诸如嵌套因式分解之类的预调节。在本发明之前，在由未构造栅格构造块结构产生一个显著问题。本发明的方法通过排序和收集在未构造栅格中的节点以产生在系数矩阵M内的块矩阵结构而克服这种困难，系数矩阵M允许使用基于块的数值解算法，并且同时提高良好的收敛。本发明者已经发现，通过使用根据对于每对相邻单元之间的边界确定的可运输性值的排序建造的单元拓扑串，能解用于物理系统的基本矩阵方程。当形成系数矩阵M时，每串与M中的块有关。在该专利中使用的术语可运输性是指诸如物质、能量、或电荷之类的某一实体在给定时间间隔期间跨过单元边界(或单元连接)的运动的容易度或能力的一种度量。运输的实体例如能是流体的质量或体积、颗粒数量、热能量、辐射或电力。如果模拟的物理系统是油气储层，则在本发明的这种描述中使用的可运输性与透过性同名，透过性是对于熟悉本专业的技术人员熟悉的术语，作为流体在表示多孔介质内的体积的两个相邻单元之间的流动的能力的度量。透过性表示为其中k是多孔介质的有效渗透性，A是在相邻单元之间的边界的面积，及△x是流体在两个单元之间运动时必须走过的平均或特征距离。在实施本发明的方法时，第一步是把物理系统离散成彼此相邻布置的多个体积单元，以便在每对相邻单元之间具有一条边界。使用基于划分物理系统以模型化成较小单位的有限差分、有限体积、有限元、或类似方法，进行离散。随后本发明的描述主要涉及有限差分法。熟悉本专业的技术人员将认识到，联系有限元法或有限体积法也能应用本发明。当借助于有限元法应用它时，单元变成有限元素，并且当借助于有限体积法应用它时，单元变成有限体积。不管使用这些方法的哪一种，他们都把偏微分方程减小到一种有限维系统的代数方程。在储层模拟中，对于每个栅格单元建造表示岩石和用于每种流体的液体性质的有限差分方程。这些方程事实上把要分析的物理系统看作包括多个较小连续单元的体积系统。当使用有限差分和有限体积法时，较小单位一般称作单元或栅格块，而当使用有限元法时，单元一般称作元素。这些单元或元素能总计从小于一百至几百万个。在该专利中，为了表示简单，使用术语单元，但应该理解，如果模拟使用有限元法，则术语元素代替在这种描述中使用的术语单元。在本发明的实施中，单元能具有任何几何形状，如平行六面体(或立方体)或六面体(具有四个长度可以变化的垂直角边缘)、或四面体、平行四边形、梯形、或三角形。栅格能包括以正规、构造图案组织的矩形单元，或者它能包括具有以不规则、示构造图案布局的各种形状的单元，或者它能包括多个构造和未构造图案。能组装几乎假定任何形状的完全未构造栅格。所有单元最好是边界对齐的，由此避免使单元任何边长接触两个其他单元的边长。在该专利中，术语边界有时与术语连接可交换地使用。如果能有物质、能量、或电荷从一个单元至另一个的运动，则两个单元具有一种连接。在一个构造栅格中，每个单元具有它与其连接的固定数量个相邻单元。在一个未构造栅格中，连接的数量能从单元至单元变化。在该方法中的下一步是对于每个单元选择状态变量。状态变量是必须和足够规定系统状态的那些变量。给定状态变量，必须有可能计算单元的所有其他性质。对于储层模拟，状态变量之一几乎总是压力。其他能包括诸如饱和、物种浓度、及物种量。为了简单，随后的讨论把除压力之外的状态变量简单地称作饱和，要理解他们能包括可以不包括饱和的各种物理性质。这些性质能从实际储层数据整体或部分地得到，或者他们能依据进行的储层模拟的类型和实际储层数据的适用性经验地确定或估计。熟悉本专业的技术人员能容易地确定适当状态变量的确定和估计其初始值。本发明的描述假定，解决一个时间依赖问题。然而，有时希望解决稳态问题。在该描述中公开的原理也能应用于稳状问题。象时间依赖问题，稳态问题涉及一次或多次解一个矩阵方程。对于每个单元，构造把一个单元的状态变量与其相邻单元的状态变量相联系的线性方程。构造这些方程以表示在每个单元内质量、能量、或动量的守恒和单元之间的质量、能量、或动量运动的基本原理。在储层模拟中，出现在非线性有限差分方程中的非线性项被线性化，并且基于这种线性化，构造代数方程的一个线性组。这些方程能依据对于模拟操作选择的方法显著变化。已经提出的用来模拟储层的方法主要不同之处在于，他们如何处理储层状态变量(如压力和饱和)按时间变化的方式。在这些方法的多种中，状态变量的值是未知的，直到已经完成用于时间步长的计算。结果，必须使用迭代过程确定他们。一种用来模拟储层的普通使用过程叫做隐式压力显式饱和度法(IMPES法)。在IMPES法中，根据压力在其时间步长末端处的值和饱和在其时间步长处的值计算在相邻单元之间的流动。在这种方法中，在时间步长末端处的压力是相互依赖的，并且必须同时确定。该方法叫做“隐式的”，因为每个压力取决于仅隐含知道的其他量(例如，在时间步长末端处的其他压力)。基本过程是通过转换方程的组合得到一个单压力方程。在压力按时间已经上升之后，明显地更新饱和。在计算饱和之后，能计算新的相对渗透性和毛细管压力；这些明显地用在下个时间步长处。用在储层模拟中的另一个过程叫做全隐式法，该方法隐式地处理压力和饱和。使用在每个时间步长末端处的相位压力和饱和计算流量。流量及压力和饱和解的计算涉及使用适当迭代技术解非线性方程。在解压力和饱和时，这些项的更新继续使用压力和饱和的新值。当满足收敛标准时，迭代过程终止。用在储层模拟中的又一种过程叫做顺序隐式法(SEQ法)。这种方法包括饱和的隐式处理，但不用同时求出压力和饱和。它包括两个步骤。第一步骤以与在IMPES法中进行的准确相同的方式解一组压力方程。该组包括每个单元一个单方程，并且解它在时间步长末端处产生一个完全、新的压力分布。在一个第二步骤中，压力分布用来计算在单元之间每条边界处的所有相位的速度和。这些所有速度在建造一组饱和方程时使用。这组包括在三相情况下的每个单元两个方程和在两相情况下的每个单元一个方程，并且同时解以产生在新时间的饱和。第二步骤是使用线性化隐式速度对于饱和的一个隐式解。通过使用相对渗透性的隐式(时间步长的末端)线性化值和在单元间流体流动项中的毛细管压力，确定在每个单元中的饱和。该方法需要所有饱和方程的同时解。非线性方程的线性化和在解方程时使用的步骤彼此依赖。在线性化的过程中，代数方程将具有依赖于选择的解技术的不同形式。例如，IMPES仅线性化压力依赖项，如比容。比容因此表示为压力的线性函数。SEQ法线性化关于压力的相同压力依赖项，并且它也线性化关于饱和的相位粘滞流动项。全隐式法线性化关于压力的压力依赖项和关于饱和的饱和依赖项(包括相对渗透性和毛细管压力)。有可能以非迭代方式使用这些方法的任一种，其中解线性化方程给出在每个时间步长末端处的解。然而，就全隐式法而论，很少这样做。而是用于时间步长的全隐式解通常使用牛顿-莱甫森迭代得到，其中解线性化方程产生一个近似解。重复牛顿-莱甫森迭代，直到根据预规定收敛标准认为生成的解估计足够精确。熟悉本专业的技术人员能进行适当模拟方法的选择和用来模拟物理系统的适当线性方程的建造。本发明不限于IMPES、全隐式、或SEQ模拟法。在本发明的实施中能使用其他已知的模拟方法、和没有发现的模拟方法。用来建造储层的数学模型的方法的例子在Peaceman，D.W.的数值储层模拟的基础(FundamentalsofNumericalReservoirSimulation)，ElsevierScientificPublishingCompany，Amsterdam(1977年)；和Mattax，C.C.和Dalton，R.L.的储层模拟(ReservoirSimulation)，MonographVolume13，SocietyofPetroleumEngineers(1990年)中描述。本发明的下一步是把一个可运输性值与在单元相邻对之间的每条边界(或连接)相联系。可运输性值与在单元之间的每个连接的耦合强度相对应，这是一种连接如何彼此强烈地耦合两个连接单元的度量。如果两个单元强烈地耦合，则他们彼此具有强烈的连通，在一个单元处状态变量的变化将对另一个单元中的状态变量有显著影响。如果两个单元耦合较弱(弱连接)，则在一个处的变化对另一个没有什么影响。对于使用有限差分的载有流体多孔介质的模拟，耦合强度能认为是连接的可运输性。对于其他物理系统的模拟操作，耦合强度可以与其他已知或容易确定的物理量相对应。对于某种模型化，由矩阵方程的系数能直接确定耦合强度。熟悉本专业的技术人员对于分析的物理系统能够确定在单元之间的耦合强度的适当度量。一旦确定可运输性值，则从具有最大强度的一个至具有最小强度的一个排序单元连接(耦合强度)。在这样做时，耦合强度的联系能以任何适当方式打破。最好，使用一种适当排序过程进行连接强度的排序。一种最佳排序过程使用QUIKSORT算法，该算法在WilliamH.Press，SaulA.Teukolsky，WilliamT.Vetterling和BrianP.Flannery的书“数值方法(NumericalRecipes)”，SecondEdition，CambridgeUniversityPress(1994年)中描述。根据在单元之间可运输性值的排序，然后构造单元的拓扑一维单元串。构造诸串以包含尽可能多的最强连接(最高可运输性值)。在最高排序可运输性值(即最强连接)处开始，在它连接的两个单元之间创建一种串中连接。然后选择次最高排序可运输性值，并且把一个第二串连接放置在它连接的两个单元之间。循环重复该过程，直到对于在单元串中的可能包括已经考虑到所有单元连接。在这种方法中，允许每个单元没有多于两个串连接。如果连接的单元之一已经具有两个串中连接，则该连接不能添加到一个串。每个单元能连接到在相同串中的最多两个其他单元上。因此，没有多于单元相邻的两个能在相同串中。相邻的之一将位于串中上方或在它之后，而一个位于下方或在它之前。几乎所有单元都具有对于两个相邻单元的串中连接。位于串末端处的单元将具有仅对于一个相邻单元的串中连接。少量几个单元可能没有对于任何相邻单元的串中连接。这些单元将形成单独单元串。包括多个单元的串形成拓扑一维直线，但该线在物理上不必是直的。在已经创建诸串之后，串的一些能并且可能接触他们本身。如果串包含连接到串中多于两个其他单元上的一个单元，则串接触它本身。另外，可能创建一种是圆形的串。在用于串构造中的最佳规则中，既不允许圆形串也不允许接触本身的串。如果这些条件之一出现在串中，则切断该串。尽管能应用各种切断过程，但如下描述一种最佳过程。如果串是圆形的，则能在任何处切断串，但最好在串的最低排序串中连接处进行切断。为了切断接触本身的非圆形串，在串的顶部处开始，沿着串进行，及确定在每个单元处一个单元是否接触(连接到)相同串中的另一个单元，该单元在前面但在串内不是紧在它前面。如果一个单元接触在相同串中在前面但在串内不是紧在它前面的另一个单元，则在当前单元与它接触的单元之间的某一地方切断该串。这种辨别过程沿串的单元继续，并且直到辨别到最后单元，该最后单元接触相同串中在前面但在串内不紧在它前面的另一个单元，且辨别第一单元，该第一单元接触相同串中在后面但在串内不紧在它后面的另一个单元。最好在这两个单元之间的最低排序连接处切断该串。继续串的这种分析和切断过程，如需要的那样，直到串没有部分接触它本身。希望的最终结果是一组满足如下规则的串(1)每串除串中的那些之外没有对其本身的连接和(2)没有串是圆形的。如果串不满足这些规则，则截断该串，从而它的确满足规则。在计算机中表示诸单元时，把每个单元分配一个标识它的索引号码。每串将由这些索引的一个排序清单定义，第一索引指示在串开始处的单元，下一个索引指示在串中的下个单元，依此类推，直到指示在串结束处的单元的最后索引。事实上，单元在串中的位置由其索引在该索引清单中的位置指示。最佳串切断过程的更详细描述如下。第一步是切断任何圆形串。首先，必须找到圆形串。这通过一个使用如下过程的消除过程进行。如以上注意到的那样，把单元加索引。从具有最小索引的单元开始，1.检查每个单元，以确定是否已经把它标记为属于一个非圆形串。如果是，则进行到具有次较大索引的单元。2.如果单元还没标记为属于一个非圆形串，则确定它是否具有对于两个其他单元的串连接。如果是，则进行到具有次较大索引的单元。3.如果单元没有串连接，则它属于一个单单元串。把它标记为属于非圆形串，其串添加到串清单，及把单元的索引添加到新串的单元清单。4.如果单元具有一个串连接，则它形成下个串的开始。把单元标记为属于一个非圆形串，把串添加到串清单上，初始化串的单元清单，及然后把单元的索引添加到该单元清单上。从一个单元到下一个通过跟随其串连接而追踪该串，把每个单元标记为属于一个非圆形串，并且把每个单元的索引添加到串的单元清单。当达到不具有对于另一个单元的串连接的单元时，这是串的末端。重复这些步骤，直到已经检查所有单元。在这时，还没标记为属于非圆形串的任何单元属于一个圆形串。一旦已经辨别到一个圆形串，下个步骤就是在其最弱连接处截断每个圆形串。从具有最小索引的单元开始，1.检查每个单元，以确定是否已经把它标记为属于一个非圆形串。如果是，则进行到具有次较大索引的单元。2.如果单元不属于非圆形串，则它属于圆形串。从一个单元到下一个通过跟随其串连接而追踪该串，保持跟踪遇到的最小可运输性和它连接哪两个单元。当达到初始单元时，已经完全横过该圆形。从连接它的两个单元除去具有最小可运输性的串连接。把这两个具有较小索引的单元看作是下个串的开始。把该串标记为属于一个非圆形串，把该串添加到串清单上，初始化串的单元清单，及然后把单元的索引添加到该单元清单。从一个单元到下一个通过跟随其串连接而追踪该串，把每个单元标记为属于一个非圆形串，并且把每个单元的索引添加到串的单元清单。当达到不具有对于另一个单元的串连接的单元时，这是串的末端。3.重复以上过程，直到已经检查所有单元。在这时，所有单元属于非圆形串。下个步骤是截断“接触”本身的；即经非串连接连接到他们本身上的任何串。这一次一串地进行。从串中具有最小索引号码的第一单元处开始，1.确定除其串连接之外的单元连接的任一个是否把它连接到串中的另一个单元上。如果不是，则前进到串的单元清单中的下个单元。2.如果单元具有对于串中其他单元的非串连接，则确定离串开始最近的连接单元。在串中初始化对于当前单元位置的一个位置P1和对于连接单元位置的一个第二位置P2。在这两个位置之间的任何处截断该串。3.运动到串中的下个单元。确定该新当前单元是否具有对于串中其他单元的非串连接。如果是，则把P1设置到当前单元位置。确定最靠近串开始的连接单元。如果它比P2靠近开始，则把P2设置成等于连接单元的位置。4.如果在串中的下个单元的位置是P2，则跳到下面的步骤5。否则，重复步骤3。5.找到在串中在P1处的单元与在P2处的单元之间具有最小可运输性值排序的连接。通过在靠近串开始的连接单元处终止串象征性地截断这种连接。其他连接单元将是在新串中的第一单元。把该串添加到串清单上，初始化串的单元清单，及然后把连接单元的索引添加到该单元清单上。从一个单元到下一个通过跟随其串连接而追踪该串，把每个单元的索引添加到串的单元清单上。当达到不具有对于另一个单元的串连接的单元时，这是串的末端。新串将在串清单的末端处。当一串一串地进行该过程时，最终将达到新的一个。在这点处，可以再截断串。如果是，则创建另一个新串。最终所有串获得处理，在这时所有串将满足一个预定组串建造规则。一旦建造诸串，就通过组成与每串的单元有关的线性方程对于每串产生一个矩阵方程。该矩阵方程的形式与在构造栅格问题中用于LSOR直线的相同。系数矩阵方程包含与在单元与其串中相邻单元之间的流动有关的项对矩阵方程的右手侧产生影响。然后通过一次一串地解用于每串的矩阵方程得到单元的状态变量的改进估计，直到已经解完所有串的矩阵方程。迭代地重复该过程，直到满足收敛条件。进行的迭代基本上与LSOR相同，不同之处在于直线是单元的串而不是常规LSOR的行或列。该方法因此称作串逐次超松驰法。诸串能以任何顺序处理，并且他们能经给定顺序向前运动而被处理，然后经相同顺序向后。这产生一种对称逐次超松驰法。随后的讨论假定常规的而不是对称的逐次超松驰。熟悉本专业的技术人员能够建造该方法的对称逐次超松驰形式。一旦已经创建一组串，就按如下得到在牛顿迭代或时间步长上的解变化。首先，组成用于每串的方程组。然后计算初始余数，如果他们还不知道的话。这些必须包括把串的单元连接到其他串中的单元上的项的影响。然后进行迭代，每次迭代包括如下步骤。1.解串的矩阵方程，把串的当前余数用作右手侧。2.把在步骤1中得到的解变化乘以一个位于一与二之间的超松驰参数ω。如果使用Orthomin加速，则收敛速率通常仅稍微取决于选择的值，最佳值通常在1与1.5之间。如果不使用Orthomin，则最佳值通常稍小于2，并且收敛速率对选择的值更敏感。3.通过把他们乘以量1-ω更新串的余数。4.对于当前串处的解变化更新在连接到当前串上的所有串处的余数。在这样做之后，所有串的余数将与当前解估计相一致。熟悉本专业的技术人员熟悉这样的计算。5.对于每串进行步骤1-4。能以任何顺序处理诸串，但对于每次迭代应该使用相同顺序。6.使用辅助校正选择性地加速收敛，如下面描述的那样。7.使用Orthomin或另一种Krylov子空间法选择性加速收敛，如下面描述的那样。8.检查由收敛测量证实的收敛是否小于预定标准。重复以上迭代步骤1-8，直到得到满意的收敛。在一个最佳实施例中，本发明迭代法的收敛通通借助于常规LSOR使用的相类似的辅助校正提高。一种最佳添加校正的描述在J.W.Watts的、标题为“一种适用于各向异性问题的迭代矩阵解法(AnIterativeMatrixSolutionMethodSuitableforAnisotropicProblems)”的论文中描述，该论文出现在石油工程师协会期刊(SocietyofPetroleumEngineersJournal)卷11、1971年3月第47-51页中。为了应用辅助校正，首先必须通过对每串求和方程建造一个校正矩阵方程。借助于LSOR使用辅助校正的熟悉本专业的技术人员能够建造该方程。使用如下步骤然后能应用辅助校正。(6a).在每串中的单元上求和余数。如果正在使用多于一个物类的守恒方程，则对于这些方程的每一个在串中所有单元上求和余数。(6b).把来自步骤(6a)的求和余数用作右手侧解校正矩阵方程。解对于解的每个未知量将包括对于每串的一种辅助校正。(6c).对于每串把在步骤(6b)中确定的辅助校正添加到在串内每个单元处的未知量上。(6d).对于对于所有串中的所有单元计算新的余数。在另一个实施例中，Orthomin法也能加速收敛或基于正交化和最小化的某种其他方法。Orthomin法属于Krylov子空间法的种类，其中把解映射到一个Krylov子空间上。根据得到的整体解变化应用Orthomin加速过程。这通过把在以上迭代步骤1-5中确定的变化添加到在以上辅助校正步骤(a)至(d)中确定的变化上执行。Orthomin由P.K.W.Vinsome在标题为“Orthomin，一种用来解联立线性方程的稀疏带组的迭代方法(Orthomin，anIterativeMethodforSolvingSparseBandedSetsofSimultaneousLinearEquations)”的论文中描述，论文号SPE5729，呈现在关于储层性能模拟的第四SPE论文集中，LosAngeles1976年2月19-20日。也见Saad，Y.，1989年的“关于超级计算机的Krylov子空间法(Krylovsubspacemethodsonsupercomputers)”，SIAMJ.Sci.Stat.Comput.，10，p.1200-1232。最佳实施例使用Orthomin，但能使用其他加速方法，如GMRES，在由Saad，Y.和Schultz，M.H.，公开的报告“一种用来解非对称线性系统的一般化最小余数算法(AGeneralizedMinimumResidualAlgorithmforSolvingNonsymmetricalLinearSystems)”，TechnicalReport254，YaleUniversity，1993年中描述。Orthomin计算包括如下步骤(7a)计算由Orthomin使用的参数。(7b)使用这些参数，更新解估计。(7c)对于所有串中的所有单元计算新余数。用在步骤6中的校正矩阵方程具有与原始矩阵方程相同的形式。结果，使用以上迭代能解它。这样做涉及诸串的建造串。迭代解产生用于所有单元的状态变量，这些状态变量在与使用的预定收敛标准相对应的精度内同时满足用于所有单元的线性方程。改进解然后能用来模拟物理系统的至少一个特征。如果物理系统是一个储层，则模拟的特征能包括例如油压力、水压力、油饱和、及水饱和。从这些变量能导出其他特征，如油生产率和水生产率。对于多个时间步长能重复迭代计算，并且结果能用来预计物理系统的性质和在其中作为时间的函数出现的运输现象。现在参照附图将描述本发明的方法。作为有助于读者理解本发明的背景信息，就图1和2而论呈现直线逐次超松驰(LSOR)的原理的简短讨论。图3表明未构造栅格系统的一个简化例子，在本发明的方法之前的该例子在模拟操作中没有使用LSOR。图4-9提供栅格系统映象的例子，在描述用来构造适于在模拟中应用LSOR原理的串或单元的直线的最佳过程时参考这些例子。图1表明已经划分成组织成5行(a、b、c、d、及e)和10列(1至10)的50个单元的物理系统的简化两维笛卡儿模型。对于基于图1的单元的模拟，LSOR能应用于形成行或列的单元的线。如果LSOR线是列，并且如果假定模拟计算从左至右进行，则第一步是计算在第一列中的改进解，把第二列中的解保持固定在其当前估计处。LSOR法的第二步计算在第二列中的解，把在第一列中的解保持固定在其在第一步计算的当前估计处，并且也把在第三列中的解保持固定在其当前估计处。以后步骤计算在第三列中、在第四列中的改进解，依此类推，直到计算在所有列中的改进解。该过程构成一次LSOR迭代。重复它，直到得到希望精度的解。在LSOR中，线的方位是重要的。当通过列或通过行进行时LSOR是否收敛得最快主要取决于在行内和在列内单元之间的耦合强度。如果单元之一的状态变化强烈地影响第二个的状态，则在两个单元之间的耦合强烈，而如果在第一单元中的这种变化对第二个没有什么影响，则它较弱。在模拟一个储层时，具有跨过单元之间的边界的较大可运输性的两个单元认为强烈地耦合。如果由位于最强耦合方向上的线进行，则LSOR通常收敛得最快。图1中的单元具有大于高度的宽度的事实指示耦合在列内比在行内强，因为在两个单元之间耦合的强度一般与适用于在他们之间的运输的横截面积直接成正比，而与在他们中心之间的距离成反比。当耦合在列内比在行内强时，即对于1描绘的单元的情形，LSOR在其线是列时比他们是行时收敛得快。在储层模拟时，如果线方位在高运输性的方向上，则迭代收敛速率较快，这对于正规、构造栅格系统常常是沿单元的列取向的单元。已知LSOR收敛能通过应用辅助校正加速。辅助校正当耦合在一个方向上比在其他一个或多个方向上强烈得多时和正在确定诸如在热传导问题中的每个单元的温度之类的单个未知量时最有效。如果LSOR通过列进行，则辅助校正是添加在一列单元中的每个温度上的量。通过求和在一列单元内的方程得到计算辅助校正需要的每个方程，这事实上确定如果把单元列看作单个单元则应用的方程。最强耦合的方向有时能在空间中变化。这种方向变化能由图2中所示的栅格单元表明，图2表明已经划分成组织成5行(a、b、c、d、及e)和15列(1至15)的75个单元的物理系统的两维模型。如在图1中那样，图2单元的几何形状指示耦合的强度。在单元之间的边界越大，在单元之间的耦合越大。在左端(列1)上，耦合在行内最强，而靠近右端(列15)，耦合在列内最强。LSOR方位的任何可能选择表示一种调和。LSOR在这样一种模型中可能缓慢地收敛。图3表明未构造单元栅格的一个简化例子。它叫作未构造的，因为其单元不都具有相同形状，并且其连接性对于所有单元不按照固定图案。该栅格不包含单元的线，即不包含列也不包含行，在这些线上自然应用LSOR。如果LSOR用来解用于这种未构造单元的方程，则必须首先改进LSOR过程。本发明者已经发现一种根据LSOR的原理能用作一种解法的用来产生线(或串)的新颖方法。参照图4、5和6现在将描述用来构造串的一个最佳过程，这些图表明编号201至215的15个单元的拓扑一维映象图。每对相邻单元在其之间具有一条边界，对于15个单元总共有22条边界。按照促进在串中包括尽可能多的具有较大可运输性值的单元的规则建造诸串。首先通过任何适当的装置确定可运输性值，并且排序可运输性值，使最大可运输性值具有一的排序至最小可运输性值具有最低排序。因而从最强单元耦合至最弱排序边界(或连接)。在图4-6中，因为有22条边界，所以可运输性排序在从1至22的范围内。在图4-7中，分配给每条边界的号码表示在每条边界处的可运输性排序。例如，在块207与208之间的边界具有最大可运输性值，并因此分配1的排序。在单元208与209之间的边界具有第二大可运输性值，并且分配2的排序。对于所有22条边界重复该过程。用来构造单元串的最佳规则是拓扑地形成包含尽可能多的最高排序可运输性值的相邻单元的一维体。该规则按如下执行接合在最高排序连接任一侧的两个单元，然后接合在次最高排序连接任一侧的两个单元，及以这种方式循环地进行，总是把接合在最高排序剩余连接任一侧的两个单元，除非在连接任一侧的一个或两个单元已经接合到两个另外的单元上，直到已经取尽连接清单。如果在边界任一侧的任何单元以前接合到两个单元上，则该边界不会形成为建造一个串目的的连接。把该串构造规则应用于图4的栅格，因为在单元207与208之间的连接具有最高可运输性值排序，所以首先接合这两个单元。其次，接合单元208和209，因为在他们之间的连接具有排序号2。排序3位于214与215之间，从而以后接合这两个单元。借助于认为是可能串连接的所有单元连接，循环地继续该过程。尽管在图4中没有表示，但该构造过程可以生成没有对于相邻单元连接的几个单元，在这种情况下，这种独立的单元形成一个单元的串。对于图4的单元使用该过程，结果是图5中所描绘的单串。图5表示由位于一条拓扑一维线上的单元组成的一个串40。串40的分析表示串接触其本身。如在本发明的描述中使用的那样，如果一个串的给定单元与在该串中的一个第二单元具有一条边界，并且第二单元不是紧在给定单元之前或之后的单元，则串接触其本身。因此，按照用于串构造的最佳规则，串40需要切断。使用上述的切断过程，切断过程通过分析在一端开始的串40的单元而开始。从单元207开始，它接触单元212、206、及202；在这些中，单元212最靠近串的开始。位置P1指向单元207，而位置P2指向单元212。其次，考虑单元208。它接触单元203和213。P1现在指向单元208，而P2指向单元213，因为213比212更靠近串的开始。其次，处理单元209。它接触单元204和214；P1现在指向209，而P2指向214。最后，单元210接触单元205。P1现在指向单元210。P2不变，因为单元205不比单元214更靠近串的开始。单元215不接触另一个单元。在单元210与214之间在具有最低可运输值的连接处进行切断。在单元210与214之间的连接具有3和9的可运输性排序。由于最低排序是9，在单元210与单元215之间的连接处进行切断。图6表示在切断图5的串40之后形成两个串50和51的最终结果。参照图6，串50包括单元205、204、203、202、201、206、211、212、213、214、及215的拓扑一维线，而串51包括单元207、208、209、及210的拓扑一维线。在串50和51内，把不在串中的最高排序连接排序为第九(在单元210与215之间)。十一个最高排序连接的九个位于串50和51内。这两个串满足在串中包括尽可能多的具有最高排序可运输性值的连接的目标。上述构造方法学的应用总结在下面关于图4-6的15个单元的表1中。表1的第四列指示串连接是否成为一条串的部分(在执行任何串切断之前)。例如，连接排序号码1与在单元207与208之间的边界相对应，并且因为这是第一串连接，所以在该边界任一侧的单元以前都没有接合到多于一个单元上。该过程应用于连接排序1至22。连接的一些不会成为在串内的连接。例如，参照表1，连接排序号12(在单元206与207之间的边界)不能放置在串内，因为单元206以前接合到单元211和201上。因此认为单元206是满的。换句话说，单元206和207不能是单元串中的相邻单元。类似地，在单元203与208之间的连接排序20不能是在单元串中的连接，因为单元203和208是满的；单元203以前接合到单元202和204上，而单元208以前接合到单元207和209上。图7、8、和9表明编号301至318的18个单元的映象图。每对相邻单元在其之间具有一条边界，对于18个单元总共有27条边界。以与关于图4-6在以上描述的排序过程相类似的方式，把编号1至27的可运输性值的排序分配给边界。图7、8、和9表示与每条边界有关的可运输性值排序。用来建造图5的串40的相同串建造规则用来由图7中描绘的18单元栅格建造串。串建造的结果表示在图8中。形成两个串60和61。串60包括单元313、307、301、302、308、及314的拓扑线，而串61包括单元310、311、312、306、305、304、303、309、315、316、317、及318的拓扑线。使用在产生表1的数据时使用的相同串建造规则，把为建造串60和61的单元排序的逐个单元分析总结在下面的表2中。一旦建造串60和61，下一步就考虑是否需要切断诸串。首先处理串61，因为其第一单元是号码310，而串60从单元313开始。再参照图8，分析串61以确定它是否需要切断。分析从单元310开始。位置P1和P2依次指向单元310和304、311和305、及312和305。在单元305与312之间的最弱耦合处切断串61。对于这种切断的潜在耦合具有16和15的可运输性排序，16是最弱耦合。为了促进包括最高排序可运输性值，在单元305与305之间的连接-最低排序连接-处进行切断，以形成表明在图9中的两个串72和73。通过对于串72中的单元进行类似分析，确定不必进一步切断，因为串72不接触它本身。其次，分析串60确定它是否也必须切断。分析在单元313处开始。位置P1和P2指向单元313和314，并且然后指向307和308。在单元307与308之间，合格切断的连接具有可运输性排序5、22、和9。由于22是这三个的最低排序，所以在单元301与302之间切断串60，以按图9中所示形成两个串70和71。本发明不过度限于为说明目的已经叙述的上文。相反，各种各样的修改和可选择实施例对于熟悉本专业的技术人员是显然的，而不脱离本发明的、在下面叙述的权利要求书中定义的真正范围。表1表权利要求1.一种模拟一个物理系统的至少一个特征的方法，包括步骤(a)把物理系统离散成彼此相邻布置的多个体积单元，以便在每对相邻单元之间具有一条边界；(b)对于每个单元分配状态变量的一个初始估计；(c)对于每个单元建造使其状态变量与同其相邻的单元的状态变量相关的线性方程；(d)把每条边界与一个可运输性值相联系，并且对应于下降可运输性值来排序边界；(e)使用边界排序以拓扑地建造一维单元串；(f)通过组合与每串的单元有关的线性方程产生一个用于每串的矩阵方程；(g)通过一次一串地解用于每串的矩阵方程得到单元状态变量的改进估计，直到已经解完所有串的矩阵方程；(h)重复步骤(g)，直到满足收敛条件，由此得到对于所有单元的状态变量，这些状态变量同时满足对于所有单元的线性方程；及(i)使用在步骤(h)中确定的状态变量，模拟物理系统的至少一个特征。2.根据权利要求1所述的方法，其中串的建造使用一种规则，该规则促进在相同串中的单元之间的单元边界处具有高可运输性值的串的建造。3.根据权利要求1所述的方法，进一步包括步骤(i)对应于可运输性值的相对大小来分级排序边界，以从具有最高可运输性值的边界至具有最低可运输性值的边界的下降顺序排序边界；和(ii)根据第一规则和第二规则建造在步骤(e)中的串，第一规则促进包括具有尽可能多的与单元之间的边界有关的最高排序可运输性值的单元，而第二规则要求串中的任何给定单元与同一单元拓扑串中的不多于两个其他单元具有边界。4.根据权利要求1所述的方法，其中对于在步骤(f)中使用的任何给定串，除紧在给定串中的给定单元之前或之后的单元之外，该给定串中没有单元连接到该串中的单元上。5.根据权利要求1所述的方法，其中在步骤(e)中建造的串是圆形的，该方法进一步包括在单元的圆形串中的最低排序边界处切断串的步骤。6.根据权利要求1所述的方法，其中在步骤(e)中建造的给定串包含一个单元，该单元接触在给定串中既不紧在该单元之前又不紧在该单元之后的另一个单元，该方法进一步包括以下步骤在给定串的一端处开始并且逐个单元地确定在给定串中的一个给定单元是否接触在给定串中既不紧在给定单元之后又不紧在其之前的另一个单元，辨别在串中接触串中另一个单元的最后单元，辨别在所述最后单元与在接触串中另一个单元的最后单元之前的最靠近单元之间的最低排序边界，及在所述最后单元与在接触串中另一个单元的最后单元之前的最靠近单元之间的最低排序边界处切断串，由此由给定串形成两个新串。7.根据权利要求1所述的方法，其中对于线逐次超松驰法建造对于步骤(g)中的对每串组合的矩阵方程。8.根据权利要求7所述的方法，进一步包括把一种辅助校正应用于在每串内步骤(f)的改进估计解的步骤，确定辅助校正，从而其应用使在每串内的余数之和成为零。9.根据权利要求1所述的方法，其中迭代解由一种Krylov加速法加速。10.根据权利要求9所述的方法，其中加速法使用Orthomin法。11.根据权利要求9所述的方法，其中加速法使用GMRES法。12.根据权利要求1所述的方法，其中在步骤(e)中的串建造形成圆形串，一个附加步骤在具有最低排序可运输性值的边界处打断每个圆形串。13.根据权利要求1所述的方法，其中物理系统包括一个含油气储层、从地面向储层延伸的井、在地面处的载油气管线、及油气处理设施。14.根据权利要求1所述的方法，其中物理系统包括一个地下储水层。15.根据权利要求1所述的方法，模拟的特征是在物理系统中的热传递。16.根据权利要求1所述的方法，其中模拟的物理系统是一个含油气层。17.根据权利要求1所述的方法，其中解的矩阵方程由有限差分近似的使用生成。18.根据权利要求1所述的方法，其中解的矩阵方程由有限元近似的使用生成。19.根据权利要求1所述的方法，其中解的矩阵方程由有限体积近似的使用生成。20.根据权利要求1所述的方法，其中体积单元包括多个未构造单元。21.根据权利要求1所述的方法，其中体积单元包括构造和未构造单元。22.根据权利要求1所述的方法，进一步包括步骤(i)辨别具有一个与串的多于两个单元具有一条边界的单元的串；和(ii)切断该串，由此形成两个串。23.根据权利要求1所述的方法，其中串的切断是在具有最低可运输性值位于彼此接触的两个单元之间的边界处，并且单元和在这两个单元之间的边界在一维单元串中。24.根据权利要求1所述的方法，其中步骤(c)的线性方程使单元在一个时间间隔末端处的状态变量与同它相邻的单元也在该时间间隔末端处的状态变量相联系。25.一种预测一个包含多种流体的物理系统的一个特征的方法，包括步骤(a)把物理系统离散成彼此相邻布置的多个体积单元，以便在每对相邻单元之间具有一条边界；(b)对于每个单元分配状态变量的一个初始估计；(c)对于每个单元，建造表示单元中的流体在一个时间间隔上的特性的基本方程，所述方程使用在时间间隔末端计算的流体和运输性质；(d)通过把基本方程线性化建造线性方程；(d)把每条边界与一个可运输性值相联系，并且对应于下降可运输性值来排序边界；(e)使用边界排序以拓扑地建造一维单元串；(f)通过组合与每串的单元有关的线性方程产生一个用于每串的矩阵方程；(g)通过一次一串地解用于每串的矩阵方程得到单元状态变量的改进估计，直到已经解完所有串的矩阵方程；(h)重复步骤(g)，直到满足收敛条件，由此得到对于所有单元的状态变量，这些状态变量同时满足对于所有单元的线性方程；及(h)使用步骤(g)的结果以预测物理系统和它包含的流体在时间间隔末端处的一个特征；及(i)对于多个时间间隔进行步骤(b)至(h)，并且使用结果预测物理系统和它包含的流体的性质作为时间的函数。26.根据权利要求25所述的方法，其中物理系统是一个地下地层。27.根据权利要求25所述的方法，其中地下地层包含油气流体。28.根据权利要求25所述的方法，其中物理系统包括与来自地下含油气地层的油气生产有关的含流体设施。29.根据权利要求25所述的方法，其中含流体设施是表面流动管线和井孔管道。30.根据权利要求25所述的方法，其中步骤(h)的结果用来预测物理系统中流体的压力和饱和。31.根据权利要求25所述的方法，其中单元是有限差分栅格单元，而基本方程是有限差分方程。32.根据权利要求25所述的方法，其中未构造单元。33.根据权利要求25所述的方法，其中构造单元。34.根据权利要求25所述的方法，其中单元是有限元素，并且基本方程是有限元方程。35.根据权利要求25所述的方法，其中单元是有限体积，并且基本方程是有限体积方程。全文摘要提供一种用来模拟一个物理系统的至少一个特征的方法。把物理系统离散成在其之间具有边界的多个体积单元(201－215)。对于每个单元分配状态变量的一种初始估计。对于每个单元,构造把其状态变量与其相邻的单元的状态变量相联系的线性方程。把一个可运输性值与每条边界相联系,并且然后对应于下降可运输性值来排序边界。边界排序用来拓扑地构造一维单元串(50、51)。对于每个串通过组合与单元有关的线性方程,建造用于每串的一个矩阵方程。通过解该矩阵方程确定单元状态变量的改进估计。这些状态变量然后用来模拟物理系统的至少一个特性。文档编号G01V1/28GK1345429SQ00805750公开日2002年4月17日申请日期2000年3月15日优先权日1999年3月31日发明者詹姆斯·W·沃茨三世申请人:埃克森美孚上游研究公司