CN100489558C

CN100489558C - 用于求解隐式储层仿真矩阵的方法

Info

Publication number: CN100489558C
Application number: CNB2005800186194A
Authority: CN
Inventors: J·W·瓦特三世
Original assignee: Exxon Production Research Co
Current assignee: ExxonMobil Upstream Research Co
Priority date: 2004-06-07
Filing date: 2005-04-13
Publication date: 2009-05-20
Anticipated expiration: 2025-04-13
Also published as: WO2005121840A1; US7672818B2; NO20070090L; EP1759226A1; CA2569102A1; US20070255779A1; CN1965248A

Abstract

一种用于求解矩阵方程AX＝B的方法，其中A代表块稀疏矩阵(205)，B代表右手侧块向量(205)，X代表解块向量(298)。在一个实施方式中，所述方法包括接收块稀疏矩阵(205)和右手侧块向量(205)，从块稀疏矩阵(205)和右手侧块向量(205)构造简化的变换块向量(205)，并且使用简化的变换块稀疏矩阵(265)和简化的变换剩余块向量(275)，来求解所述解块向量(298)。

Description

用于求解隐式储层仿真矩阵的方法

技术领域

【0001】本发明的实施方式一般涉及地下储层中碳氢化合物的开采和开发，更优选地是涉及一种用于预测地下具有碳氢化合物的地层的状态的改进过程。

背景技术

【0002】储层仿真是从储层模型的性能推测实际储层状态的过程。因为石油储层中的质量传递和流体流动过程十分复杂，所以储层仿真是用计算机来进行的。执行计算以仿真储层的计算机程序称为储层模拟器或模拟程序。储层仿真的目的是充分理解发生在石油储层中的复杂的化学、物理和流体流动过程，以能够预测储层的未来状态并最大化碳氢化合物的回收。储层模拟器能解决通常用任何其它方式不可解决的储层问题。例如，储层模拟器能预测储层管理决策的结果。储层仿真通常指储层中流动的流体力学，但从更广泛的意义讲，它也指包括储层，地表设施和任何相关的重要工作的全部石油系统。

【0003】图1示意地说明了在石油储层的储层仿真的一个实例中的四个基本步骤。第一步(步骤1)是基于储层中发生的化学，物理和流体流动过程构造实际储层的数学模型。该数学模型包括一组非线性偏微分方程。第二步(步骤2)包括在时间和空间上对所述储层进行离散化。通过将储层划分成合适的栅格单元来对空间进行离散化，每一栅格单元具有一组非线性有限差分方程。第三步(步骤3)是线性化出现在非线性有限差分方程中的非线性项，并基于此线性化，构造组合为矩阵方程的线性代数方程。第四步(步骤4)是求解组合为矩阵方程的线性代数方程。仿真在一系列时间步长内进行，且在每一时间步长中都执行步骤3和步骤4。仿真提供了储层状态的预测，其能使石油工程师预测包括储层产生的速率的储层性能。在该模型已用于仿真的回收过程后，模型的精确性可通过对比储层的历史进行验证。

【0004】另一方面，已经提出了许多种仿真方法。选择的方法会影响解的稳定性和精确性。在每一时间步长基础上，一些方法所需的计算工作比另一些方法要多。各方法的不同主要在于其处理储层变量(如压力和饱和度)随时间变化的方式不同。大多数方法包括对下述两个程序的变化：

1)显式过程所使用的迁移率(mobility)和毛细压力被计算成每一时间步长开始时的饱和度的函数。通过在先时间步长的计算，可得知该饱和度。假设在一个时间步长内迁移率和毛细压力保持同一个值，其与该时间步长开始时的值相同。

2)隐式过程所使用的迁移率和毛细压力被计算成每一时间步长结束时的饱和度的函数。直到完成该时间步长的计算，才能得知这些值。因此，必须使用迭代过程才能确定他们的值。

【0005】全隐式方法通常使用在隐式过程中。该方法是无条件稳定的，这是因为它隐式地(implicitly)处理压力和饱和度。流动速率是通过每一时间步长结束时的阶段压力和饱和度计算得到的。此方法中的饱和度不能降到小于零，这是因为只有在每一时间步长结束时流体为具有流动性时，它才能流动。只有在饱和度大于零时，流体才是可移动的。流动速率，压力和饱和度的解的计算包括使用合适的迭代方法求解非线性方程。一旦求解了压力和饱和度，这些项将继续使用压力和饱和度的新值进行更新。当满足收敛判据时，该迭代过程终止。

【0006】完全隐式方法的主要缺点是所需的计算时间的量。根据计算成本，该方法通常在单井模型或部分储层模型中是令人满意的，但在全部储层模型中使用时却非常昂贵。为了减少所需的计算，人们尝试了各种方法，但基本都是以接受这样一种方法为代价的，该方法不允许完全隐式方法的时间步长的长度。为减少计算时间，已经提出了顺序隐式方法，自适应隐式方法，和级联方法。然而，这些方法也有其自身的缺陷。在完全隐式方法中，消耗最长计算时间的是方程的求解步骤(图1的步骤4)。此步骤一般消耗约四分之三的总计算时间。

【0007】因此，仍需要在计算上更为高效的方法来求解完全隐式储层仿真中出现的线性代数方程。

发明内容

【0008】本发明的各种实施方式涉及用于求解矩阵方程AX＝B的方法，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量。在一个实施方式中，所述方法包括接收所述块稀疏矩阵和所述右手侧块向量，根据所述块稀疏矩阵来构造简化的变换块稀疏矩阵，根据所述块稀疏矩阵和所述右手侧块向量来构造简化的变换剩余块向量，并使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量求解所述解块向量。

【0009】在另一实施方式中，该方法包括根据所述块稀疏矩阵来构造简化的变换块稀疏矩阵，根据所述块稀疏矩阵和所述右手侧块向量来构造简化的变换剩余块向量，使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量求解简化的变换解变化块向量，将所述简化的变换解变化块向量转换成具有一个或更多个质量未知数变化和一个或更多个压力未知数变化的解变化块向量，以及将所述解变化块向量加上所述解块向量的当前估计值，来更新所述解块向量。

附图的简要描述

【0010】图1示出了在示例储层仿真过程中的基本步骤的示意图。

【0011】图2图示说明根据本发明的一个或更多个实施方式的一种在矩阵中求解一个或更多个线性代数方程的方法。

【0012】图3图示说明本发明的实施方式可以实施到的一个或更多个计算机网络。

具体实施方式

绪论和定义

【0013】下面将对本发明进行详细描述。每一附加权利要求定义一个单独的发明，为防止侵权，认为每个发明包括与权利要求中说明的各部分或各限制等同的部分。根据上下文，在某些情况下，下文所使用的“发明”只指某些特定的实施方式。在其它情况下，所使用的“发明”指在一个或更多个权利要求中，但并不一定是所有的权利要求中叙述的主旨。下文将对每个发明进行详细描述，包括对特定的实施方式，形式和例子进行详细描述，但是本发明并不局限于这些实施方式，形式或例子，本领域一般技术人员通过这些实施方式及本专利提供的信息并结合现有的信息和技术能理解并使用专利的发明。这里使用的术语见以下定义。对在权利要求中使用但下文没有定义的术语，应该给该术语赋予最广泛的定义，即相关领域人员已赋予该术语如出现在一个或更多个印刷出版物或发行专利的定义。

【0014】这里使用的术语“栅格单元”的定义是：组成三维储层模型一部分的一个单元或块。这样，三维储层模型可包括许多栅格单元，范围从几十几百到成千上万乃至几百万个栅格单元。在某些情况下，每一栅格单元可代表三维储层模型的一个明确指定的部分。栅格单元的全部组合可构成一个地质模型，该模型代表了感兴趣的地表下地球体积。每一栅格单元优选为代表地表下的唯一部分。这种栅格单元优选为不互相重叠。栅格单元的尺寸或维数被优选为在一栅格单元内储层性质相对为同质的，但并不产生太多的栅格单元。这些栅格单元具有范围从小于一米到几百米的边。优选地，每一栅格单元在平面图上为正方形或矩形，且具有恒定或可变的厚度。然而，可以认为可以选择使用其它的形状。可以把栅格单元看作是具有渗透性边的搅拌均匀的油箱(well-stirred tanks)。因此，栅格单元的内容可认为一致地分散在栅格单元内，流体流入或流出的速度可由栅格单元的边的渗透性和相邻栅格单元间的压力差决定。这样，数学问题就简化为相邻栅格单元间的流量计算。

【0015】这里使用的术语“奇异值分解”的定义是：用于将矩形阵分解为三个因式的数学方法。给定一个M×N阶矩阵A，这里M>N，通过使用至少某些类型的奇异值分解，矩阵A被改写成：

[\begin{matrix} A \end{matrix}] = [\begin{matrix} U \end{matrix}] [\begin{matrix} W \end{matrix}] [\begin{matrix} V^{T} \end{matrix}]

这里W为具有“奇异值”元的对角阵。矩阵U的列和矩阵V^T的行(或，相当于矩阵V的列)是正交的。举例来说，也就是如果U_i是矩阵U的一列，U_j是矩阵U的另一列，那么

U_i ^TU_i＝1

U_i ^TU_j＝0，j≠i

【0016】对方阵V而言，该等式也成立，即矩阵V的逆为矩阵V^T。矩阵U的列，矩阵W的对角元，和矩阵V^T的行通常被排列为使w₁>w₂>…>w_N，其中w_i为矩阵W的i行的对角元。矩阵W的对角元素为矩阵A的奇异值。为说明本发明的各种实施方式，将它们放置成相反的顺序会更方便，即w₁<w₂<…<w_N，并相应记录U和V的列。矩阵V叫作右矩阵，包含在其中的向量叫作右奇异向量。

【0017】基于上述定义，矩阵A可由包含矩阵U的最后几列的子集，矩阵W的元，和矩阵V^T的行进行近似。这样，矩阵A可被改写成：

[\begin{matrix} A \end{matrix}] \approx [\begin{matrix} U_{s} \end{matrix}] [W_{s}] [\begin{matrix} V_{s}^{T} \end{matrix}] .

【0018】这里使用术语“时间步长”定义为储层的寿命被离散成的时间增量。对至少某些类型的时间步长，储层模拟器计算每一栅格单元在许多时间步长的每一时间步长中的变化(流量，压力等等)。一般来说，各条件只在时间步长的开始或末端定义，在一个时间步长中的任何中间时刻不作任何定义。因此，每一栅格单元内的条件会从一个时间步长到另一个时间步长发生剧烈变化。通常，将时间步长选择为足够小以将这些剧烈变化的范围限制在可接受的限度内。时间步长的大小取决于对精确度的考虑和稳定性约束。通常，时间步长越小，解也就越精确，然而，时间步长越小，需要的计算工作也就越多。

【0019】这里使用的术语“单位矩阵”的定义是：任意维数的方阵，其从西北到东南对角上的元素全为1，其它的元素全为0。任何方阵被其维数的单位矩阵乘等于其本身。

【0020】这里使用的术语“体积约束”为一等式，该等式基于一定时间内一个栅格单元必须包含填充栅格单元所需的一定量的流体的原则。例如，如果储层包含液体碳氢化合物和水相，那么：

液体碳氢化合物体积+水体积＝栅格单元体积

【0021】这里使用的术语“未知数”定义为未知变量，其由组合在矩阵方程的线性代数方程进行求解。对包含油和水的储层，求解线性代数方程的各种未知变量，包括对压力和质量的求解。其它量可从这些变量导出。

【0022】这里使用的术语“块稀疏矩阵”的定义是：元素基本全为空或基本全为0的矩阵和其子矩阵为余项的矩阵。然而，沿其对角的所有子矩阵存在，即不为空。例如，块稀疏矩阵A可表示成：

A = [\begin{matrix} A_{11} & A_{12} & \cdot \cdot \cdot & A_{1 n} \\ A_{21} & A_{22} & \cdot \cdot \cdot & A_{2 n} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ A_{n 1} & A_{n 2} & \cdot \cdot \cdot & A_{nn} \end{matrix}]

【0023】每一子矩阵，如，A_ij可表示成：

A_{ij} = [\begin{matrix} a_{11}^{(ij)} & a_{12}^{(ij)} & \cdot \cdot \cdot & a_{1, m_{j} + 1}^{(ij)} \\ a_{21}^{(ij)} & a_{22}^{(ij)} & \cdot \cdot \cdot & a_{2, m_{j} + 1}^{(ij)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ a_{m_{i} + 1,1}^{(ij)} & a_{m_{i} + 1,2}^{(ij)} & \cdot \cdot \cdot & a_{m_{i} + 1, m_{j} + 1}^{(ij)} \end{matrix}]

这里m_i为i行中质量平衡方程的数目，m_j为j列中质量变化未知数的数目。每一子矩阵具有m_i+1行m_j+1列。第一m_i行包含质量平衡方程的系数。最底部的一行(m_i+1行)包含与体积约束方程相关的系数。第一m_j列的系数被配置为被质量变量乘，最右边的一列(m_j+1列)的系数被配置为被压力变量乘。同样地，每一子矩阵具有特殊的结构。

【0024】质量变量有几种形式。通常，质量变量是一个牛顿迭代中的变化。但也可以是一个时间步长的变化或时间步长末端的质量，和不变化。它们也可以以其它的质量的度量形式表示，如饱和度或防波堤/质量分率。类似地，压力变量通常是一个牛顿迭代中的变化，但也可以是时间步长末端的压力。从一组变量可以很容易地转换成另一组变量的这个意义而言，这些变量的选择是可互换的。为了简便起见，接下来的描述使用的是时间步长末端的质量和压力。

【0025】对于非对角子矩阵，这里i≠j，A_ij可以表示成：

这里f和

系数涉及栅格块之间的流量。f系数占据A_ij的第一m_j列，被配置为被质量变化乘。因此，它们也被称作质量变化项系数。

系数占据A_ij的m_i+1列，被配置为被压力变化乘。因此，

系数也被称作压力变化项系数。f系数用于奇异值分解，而

系数并不用于奇异值分解。

A_ij也可表示成：

这里

F_{ij} = [\begin{matrix} f_{11}^{(ij)} & f_{12}^{(ij)} & \cdot \cdot \cdot & f_{1, m_{j}}^{(ij)} \\ f_{21}^{(ij)} & f_{22}^{(ij)} & \cdot \cdot \cdot & f_{2, m_{j}}^{(ij)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ f_{m_{i}, 1}^{(ij)} & f_{m_{i}, 2}^{(ij)} & \cdot \cdot \cdot & f_{m_{i}, m_{j}}^{(ij)} \end{matrix}]

且

Φ_{ij} = [\begin{matrix} φ_{1}^{(ij)} \\ φ_{2}^{(ij)} \\ \cdot \cdot \cdot \\ φ_{m_{i}}^{(ij)} \end{matrix}]

【0026】对角子矩阵A_ii可以表示成：

A_{ii} = [\begin{matrix} 1 + f_{11}^{(ii)} & f_{12}^{(ii)} & \cdot \cdot \cdot & f_{1, m_{i}}^{(ii)} & φ_{1}^{(ii)} \\ f_{21}^{(ii)} & 1 + f_{22}^{(ii)} & \cdot \cdot \cdot & f_{2, m_{i}}^{(ii)} & φ_{2}^{(ii)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ f_{m_{i}, 1}^{(ii)} & f_{m_{i}, 2}^{(ii)} & \cdot \cdot \cdot & 1 + f_{m_{i}, m_{i}}^{(ii)} & φ_{m_{i}}^{(ii)} \\ {\tilde{V}}_{1}^{(i)} & {\tilde{V}}_{2}^{(i)} & \cdot \cdot \cdot & {\tilde{V}}_{m_{i}}^{(i)} & c_{i} \end{matrix}]

这里f和

系数涉及栅格块间的流量，

为流体分体积，且c涉及压缩率。

A_ii也可以表示成：

A_{ii} = [\begin{matrix} I + F_{ii} & Φ_{ii} \\ {\tilde{V}}_{i} & c_{i} \end{matrix}]

这里

F_{ii} = [\begin{matrix} f_{11}^{(ii)} & f_{12}^{(ii)} & \cdot \cdot \cdot & f_{1, mi}^{(ii)} \\ f_{21}^{(ii)} & f_{22}^{(ii)} & \cdot \cdot \cdot & f_{2, m_{i}}^{(ii)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ f_{m_{i}, 1}^{(ii)} & f_{m_{i}, 2}^{(ii)} & \cdot \cdot \cdot & f_{m_{i}, m_{i}}^{(ii)} \end{matrix}]

Φ_{ii} = [\begin{matrix} φ_{1}^{(ii)} \\ φ_{2}^{(ii)} \\ \cdot \cdot \cdot \\ φ_{m_{i}}^{(ii)} \end{matrix}],

且I为单位矩阵。

的元为多相分体积。它们将栅格单元中流体的体积和其中的组分的量相联系。c_i项涉及流体压缩率。它将栅格单元的流体体积和其压力相联系。

【0027】这里使用的术语“列矩阵”的定义是：在块稀疏矩阵的列中包含子矩阵的矩阵。例如，A的第j列矩阵为

A_{j} = [\begin{matrix} A_{1 j} \\ A_{2 j} \\ \cdot \cdot \cdot \\ A_{n 1, j} \end{matrix}]

【0028】这里使用的术语“解块向量”的定义是：由子向量X₁...X_n组成的块向量。例如，解块向量X可以表示成：

X = [\begin{matrix} X_{1} \\ X_{2} \\ \cdot \cdot \cdot \\ X_{n} \end{matrix}]

每一子向量可表示成：

X_{i} = [\begin{matrix} x_{1}^{(i)} \\ x_{2}^{(i)} \\ x_{m_{i}}^{(i)} \\ x_{m_{i} + 1}^{(i)} \end{matrix}]

或

X_{i} = [\begin{matrix} m_{1}^{(i)} \\ m_{2}^{(i)} \\ \cdot \cdot \cdot \\ m_{m_{i}}^{(i)} \\ p^{(i)} \end{matrix}]

或

X_{i} = [\begin{matrix} M_{i} \\ P_{i} \end{matrix}]

这里m或M代表质量未知数，p或P代表压力未知数。

【0029】在某些例子中，与解块向量X相比，可能更关心解变化向量δX。解变化向量δX可表示成：

{δX}_{i} = [\begin{matrix} δ x_{1}^{(i)} \\ {δx}_{2}^{(i)} \\ δ x_{m_{i}}^{(i)} \\ δ x_{m_{i} + 1}^{(i)} \end{matrix}]

或

{δX}_{i} = [\begin{matrix} δ m_{1}^{(i)} \\ {δm}_{2}^{(i)} \\ \cdot \cdot \cdot \\ {δm}_{m_{i}}^{(i)} \\ δ p^{(i)} \end{matrix}]

或

{δX}_{i} = [\begin{matrix} {δM}_{i} \\ δ P_{i} \end{matrix}]

【0030】这里使用的术语“右手侧块向量”的定义是：方程AX＝B的右边。例如，右手侧块向量B的子向量可表示成：

B_{i} = [\begin{matrix} b_{1}^{(i)} \\ b_{2}^{(i)} \\ \cdot \cdot \cdot \\ b_{m_{i}}^{(i)} \\ b_{m_{i} + 1}^{(i)} \end{matrix}]

或

B_{i} = [\begin{matrix} b_{M, 1}^{(i)} \\ b_{M, 2}^{(i)} \\ \cdot \cdot \cdot \\ b_{{M, m}_{i}}^{(i)} \\ b_{V}^{(i)} \end{matrix}]

或

B_{i} = [\begin{matrix} B_{M, i} \\ B_{V, i} \end{matrix}]

这里M下标表示质量平衡右边，V下标表示体积约束右边。

【0031】这里使用的术语“剩余块向量”的定义是：R＝B-Ax。例如，剩余块向量R的子向量可表示成：

R_{i} = [\begin{matrix} r_{1}^{(i)} \\ r_{2}^{(i)} \\ \cdot \cdot \cdot \\ r_{m_{i}}^{(i)} \\ r_{m_{i} + 1}^{(i)} \end{matrix}]

或

R_{i} = [\begin{matrix} r_{M, 1}^{(i)} \\ r_{M, 2}^{(i)} \\ \cdot \cdot \cdot \\ r_{{M, m}_{i}}^{(i)} \\ r_{V}^{(i)} \end{matrix}]

或

R_{i} = [\begin{matrix} R_{M, i} \\ R_{V, i} \end{matrix}]

这里M下标表示质量平衡剩余，V下标表示体积约束剩余。

【0032】这里使用的术语“质量平衡方程”的定义是：栅格单元的内容和流进栅格单元及流出栅格单元之间的数学关系。它基于物质不从该系统产生也不从该系统消失的假设。储层流体中的每一化学组分必须在每一栅格单元中满足质量平衡。例如，对于甲烷组分，在一预定的时间步长末端的一特定栅格单元中的甲烷的质量必须满足下述方程：

新时刻甲烷质量＝旧时刻甲烷质量+甲烷流进的质量-甲烷流出的质量

类似的关系施加于该储层流体中出现的任何其它化学组分。

【0033】整个储层模型的质量平衡方程可改写成：

IM+F_MM+F_PP＝B_M (1)

【0034】F_M和F_P矩阵涉及栅格单元之间的流动。M包含质量子向量，P是压力的向量，B_M包含质量平衡右手侧子向量。如果假设一维模型具有五个栅格单元，矩阵F_M可被改写成如下的块三对角线矩阵：

F_{m} = [\begin{matrix} F_{m 11} & F_{m 12} \\ F_{m 21} & F_{m 22} & F_{m 23} \\ F_{m 32} & F_{m 33} & F_{m 34} \\ F_{m 43} & F_{m 44} & F_{m 45} \\ F_{m 54} & F_{m 55} \end{matrix}]

其中，一般

F_{Mij} = - \underset{i}{Σ} F_{Mij}

F_P和F_m具有相同的结构形式。我们可以将方程(1)写成：

[\begin{matrix} M_{1} \\ M_{2} \\ M_{3} \\ M_{4} \\ M_{5} \end{matrix}] + [\begin{matrix} F_{M 11} & F_{M 12} \\ F_{M 21} & F_{M 22} & F_{M 23} \\ F_{M 32} & F_{M 33} & F_{M 34} \\ F_{M 3} & F_{M 44} & F_{M 45} \\ F_{M 54} & F_{M 55} \end{matrix}] [\begin{matrix} M_{1} \\ M_{2} \\ M_{3} \\ M_{4} \\ M_{5} \end{matrix}] + F_{P} P = B_{M}

【0035】这里使用的术语“弗罗比尼斯范数(Frobenius norm)”的定义是：矩阵中系数的平方和的平方根。例如，矩阵为：

A = [\begin{matrix} a_{11} & a_{12} & \cdot \cdot \cdot & a_{1 n} \\ a_{21} & a_{22} & \cdot \cdot \cdot & a_{2 n} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ a_{m 1} & a_{m 2} & \cdot \cdot \cdot & a_{mn} \end{matrix}]

其弗罗比尼斯范数可表示成：

{| | A | |}_{F} = \sqrt{Σ_{i = 1}^{m} Σ_{j = 1}^{n} a_{ij}^{2}}

弗罗比尼斯范数提供了矩阵A的奇异值分解中的最大奇异值的有用估计值。

特定实施方式

【0036】下面描述各种特定实施方式，权利要求中也对其中的至少一些实施作了叙述。

【0037】在至少一个特定的实施方式中，一种用于求解矩阵方程AX＝B的方法，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量，所述方法包括：接收块稀疏矩阵和右手侧块向量；根据块稀疏矩阵来构造简化的变换块稀疏矩阵；根据块稀疏矩阵和右手侧块向量来构造简化的变换剩余块向量；以及使用简化的变换块稀疏矩阵和简化的变换剩余块向量求解解块向量。

【0038】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵。

【0039】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；和对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵。

【0040】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；和对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵，其中对角阵包括一个或更多个以增序排列的奇异值。

【0041】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵，其中对角阵包含一个或更多个以增序排列的奇异值；和舍弃小于预定阈值或门限值的一个或更多个奇异值以生成简化的对角阵。预定的阈值通常设定为0.01和0.1之间的数值，但也可使用其它的值。

【0042】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵；舍弃小于预定阈值的一个或更多个奇异值以生成简化的对角阵；和用简化的对角阵乘左矩阵以生成临时的列矩阵。

【0043】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵；舍弃小于预定阈值的一个或更多个奇异值，以生成简化的对角阵；用简化对角阵乘左矩阵生成具有多个质量变化项系数的临时列矩阵；和将临时列矩阵的质量变化项系数组合到一个或更多个对应于列矩阵的质量变化项系数的临时质量变化项系数子矩阵。

【0044】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵；舍弃小于预定阈值的一个或更多个奇异值以生成简化的对角阵；用简化对角阵乘左矩阵以生成具有多个质量变化项系数的临时列矩阵；将临时列矩阵的质量变化项系数组合到一个或更多个临时质量变化项系数子矩阵，其对应于列矩阵的质量变化项系数；和用右矩阵的转置左乘每一临时质量变化项系数子矩阵，以产生变换的临时质量变化项系数子矩阵。

【0045】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵；舍弃小于预定阈值的一个或更多个奇异值以生成简化的对角阵；用简化对角阵乘左矩阵以生成具有更多个质量变化项系数的临时列矩阵；将临时列矩阵的质量变化项系数组合到一个或更多个临时质量变化项系数子矩阵，其对应于列矩阵的质量变化项系数；用右矩阵的转置左乘每一临时质量变化项系数子矩阵，以生成变换的临时质量变化项系数子矩阵；和用右矩阵的转置左乘块稀疏矩阵中包含压力变化项系数的每一子向量，产生变换的压力变化项系数子向量。

【0046】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，构造简化的变换块稀疏矩阵进一步包括构造具有与块稀疏矩阵相同的块结构和子矩阵的变换的块稀疏矩阵。

【0047】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；和对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵；舍弃小于预定阈值的一个或更多个奇异值以生成简化的对角阵；用简化对角阵乘左矩阵以生成具有更多个质量变化项系数的临时列矩阵；将临时列矩阵的质量变化项系数组合到一个或更多个临时质量变化项系数子矩阵，其对应于列矩阵的质量变化项系数；用右矩阵的转置左乘每一临时质量变化项系数子矩阵，以生成变换的临时质量变化项系数子矩阵；用右矩阵的转置左乘块稀疏矩阵中包含压力变化项系数的每一子向量，以产生变换的压力变化项系数子向量；和用右矩阵右乘块稀疏矩阵中的每一流体分体积子向量，以生成变换的流体分体积向量。

【0048】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；和对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵；舍弃小于预定阈值的一个或更多个奇异值以生成简化的对角阵；用简化对角阵乘左矩阵以生成具有更多个质量变化项系数的临时列矩阵；将临时列矩阵的质量变化项系数组合到一个或更多个临时质量变化项系数子矩阵，其对应于列矩阵的质量变化项系数；用右矩阵的转置左乘每一临时质量变化项系数子矩阵，以生成变换的临时质量变化项系数子矩阵；用右矩阵的转置左乘块稀疏矩阵中包含压力变化项系数的每一子向量，以产生变换的压力变化项系数子向量；用右矩阵右乘块稀疏矩阵中的每一流体分体积子向量，以生成变换的流体分体积向量；和从一个或更多个变换的临时质量变化项系数子矩阵，变换的压力变化项系数子向量和变换的流体分体积向量构造变换的块稀疏矩阵。

【0049】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，矩阵方程代表一维或更多维的流体流动方程的系统，其具有一个或更多个压力变化项和一个或更多个质量变化项，其中块稀疏矩阵包含一个或更多个压力变化项系数和一个或更多个质量变化项系数，和其中构造简化的变换块稀疏矩阵包括：将来自块稀疏矩阵的一个列中的质量变化项系数组合到列矩阵；和对列矩阵执行奇异值分解以生成左矩阵，对角阵和右矩阵；舍弃小于预定阈值的一个或更多个奇异值以生成简化的对角阵；用简化对角阵乘左矩阵以生成具有更多个质量变化项系数的临时列矩阵；将临时列矩阵的质量变化项系数组合到一个或更多个临时质量变化项系数子矩阵，其对应于列矩阵的质量变化项系数；用右矩阵的转置左乘每一临时质量变化项系数子矩阵，以生成变换的临时质量变化项系数子矩阵；用右矩阵的转置左乘块稀疏矩阵中包含压力变化项系数的每一子向量，以产生变换的压力变化项系数子向量；用右矩阵右乘块稀疏矩阵中的每一流体分体积子向量，以生成变换的流体分体积向量；从一个或更多个变换的临时质量变化项系数子矩阵，变换的压力变化项系数子向量和变换的流体分体积向量构造变换的块稀疏矩阵；和消去对应于变换的块稀疏矩阵中舍弃的奇异值的一个或更多个变换流体分体积，以生成简化的变换块稀疏矩阵。

【0050】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，计算列矩阵的弗罗比尼斯范数，其被称作列矩阵范数或列矩阵值。如果列矩阵值等于或小于预定的阈值，对列矩阵中不执行奇异值分解，奇异值被设定为0，右矩阵被设定为等于单位矩阵。结果，因为右矩阵等于单位矩阵，使用右矩阵的某些计算被跳过或被简化。

【0051】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，简化的变换块稀疏矩阵包括一个或更多个简化的变换对角子矩阵和一个或更多个简化的变换非对角子矩阵，其中每一简化的变换对角子矩阵仅在每一变换的对角子矩阵的底部r_i+1行和最右边r_i+1列内包括质量变化项系数和压力变化项系数。

【0052】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，简化的变换块稀疏矩阵包括一个或更多个简化的变换对角子矩阵和一个或更多个简化的变换非对角子矩阵，其中每一简化的变换对角子矩阵仅在每一变换的对角子矩阵的底部r_i+1行和最右边r_i+1列内包括质量变化项系数和压力变化项系数，和每一简化的变换非对角子矩阵仅在每一变换的非对角子矩阵的底部r_j+1行和最右边r_j+1列内包括质量变化项系数和压力变化项系数。

【0053】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，构造简化的变换剩余块向量包括构造变换的剩余块向量。

【0054】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，构造变换的剩余块向量包括构造具有变换的质量平衡剩余子向量和变换的体积约束剩余子向量的变换的剩余块向量。

【0055】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，构造简化的变换剩余块向量包括：用右矩阵的转置左乘质量平衡剩余子向量，以生成变换的质量平衡剩余子向量。

【0056】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，求解解块向量包括使用简化的变换块稀疏矩阵和简化的变换剩余块向量来求解简化的变换解变化块向量。

【0057】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，求解解块向量包括：使用简化的变换块稀疏矩阵和简化的变换剩余块向量来求解简化的变换解变化块向量；和将简化的变换解变化块向量转换成解变化块向量。

【0058】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，求解解块向量包括：使用简化的变换块稀疏矩阵和简化的变换剩余块向量来求解简化的变换解变化块向量；将简化的变换解变化块向量转换成解变化块向量；和将解变化块向量加到解块向量的当前估计值，以更新解块向量。

【0059】在等同为上述的方法或本文其它地方描述的方法的特定的实施方式中，求解解块向量包括：使用简化的变换块稀疏矩阵和简化的变换剩余块向量来求解简化的变换解变化块向量；和将简化的变换解变化块向量转换成具有一个或更多个质量未知数变化和一个或更多个压力未知数变化的解变化块向量。

【0060】在至少一个特定的实施方式中，一种用于求解矩阵方程AX＝B的方法，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量，所述方法包括：根据块稀疏矩阵来构造简化的变换块稀疏矩阵；根据块稀疏矩阵和右手侧块向量来构造简化的变换剩余块向量；使用简化的变换块稀疏矩阵和简化的变换剩余块向量来求解简化的变换解变化块向量；将简化的变换解变化块向量转换成具有一个或更多个质量未知数变化和一个或更多个压力未知数变化的解变化块向量；和将解变化块向量加到解块向量的当前估计值，以更新解块向量。

附图中的具体实施方式

【0061】现在将描述示于附图中的具体实施方式。

【0062】图2图示说明根据本发明的一个实施方式的用于求解矩阵方程形式的一个或更多个线性代数方程的方法200的流程图。在步骤205中，接收块稀疏矩阵A和块向量B。根据隐式传输矩阵方程AX＝B，使用块稀疏矩阵A和块向量B来求解解块向量X。

【0063】在步骤210，块稀疏矩阵A的j列中的f系数被组合到块矩阵M_j。组合到块矩阵M_j的f系数可表示成：

M_{j} = [\begin{matrix} f_{11}^{(jj)} & f_{12}^{(jj)} & \cdot \cdot \cdot & f_{1, m_{j}}^{(jj)} \\ f_{21}^{(jj)} & f_{22}^{(jj)} & \cdot \cdot \cdot & f_{2, m_{j}}^{(jj)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ f_{m_{j}, 1}^{(jj)} & f_{m_{j}, 2}^{(jj)} & \cdot \cdot \cdot & f_{m_{j}, m_{j}}^{(jj)} \\ f_{11}^{(i_{1} j)} & f_{12}^{(i_{1} j)} & \cdot \cdot \cdot & f_{1, m_{j}}^{(i_{1} j)} \\ f_{21}^{(i_{1} j)} & f_{22}^{(i_{1} j)} & \cdot \cdot \cdot & f_{2, m_{j}}^{(i_{1} j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ f_{m_{i_{1}}, 1}^{(i_{1} j)} & f_{m_{i_{1}}, 2}^{(i_{1} j)} & \cdot \cdot \cdot & f_{m_{i_{1}}, m_{j}}^{(i_{1} j)} \\ f_{11}^{(i_{2} j)} & f_{12}^{(i_{2} j)} & \cdot \cdot \cdot & f_{1, m_{j}}^{(i_{2} j)} \\ f_{21}^{(i_{2} j)} & f_{22}^{(i_{2} j)} & \cdot \cdot \cdot & f_{2, m_{j}}^{(i_{2} j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ f_{m_{i_{2}}, 1}^{(i_{2}, j)} & f_{m_{i_{2}}, 2}^{(i_{2}, j)} & \cdot \cdot \cdot & f_{m_{i_{2}}, m_{j}}^{(i_{2}, j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \end{matrix}]

或

M_{j} = [\begin{matrix} F_{jj} \\ F_{i_{1} j} \\ F_{i_{2} j} \\ \cdot \cdot \cdot \end{matrix}],

其中F_jj，F_ij，F_ij为在稀疏矩阵A的j列中的包含f系数的子矩阵。j列代表待处理的第一感兴趣的列。

【0064】在步骤215，对矩阵M_j执行奇异值分解以生成矩阵U_j，W_j，和V_j，其中U_j代表包含左奇异向量的第一矩阵或左矩阵，其中W_j代表包含以增序排列(即，从小到大)的奇异值的第二矩阵或对角阵，其与通常次序相反，并且其中V_j代表包含右奇异向量的第三矩阵或右矩阵。奇异值分解确保M_j＝U_jW_jV_j ^T，其中V_j ^T是右矩阵V_j的转置。左矩阵U_j和右矩阵V_j的列被排列成使对角阵W_j的每个元被配置为乘左矩阵的适当的列和转置的右矩阵V_j ^T的适当的行。左矩阵U_j可表示成：

U_{j} = [\begin{matrix} u_{11}^{(jj)} & u_{12}^{(jj)} & \cdot \cdot \cdot & u_{1, m_{j}}^{(jj)} \\ u_{21}^{(jj)} & u_{22}^{(jj)} & \cdot \cdot \cdot & u_{2, m_{j}}^{(jj)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ u_{m_{j}, 1}^{(jj)} & u_{m_{j}, 2}^{(jj)} & \cdot \cdot \cdot & u_{m_{j}, m_{j}}^{(jj)} \\ u_{11}^{(i_{1} j)} & u_{12}^{(i_{1} j)} & \cdot \cdot \cdot & u_{1, m_{j}}^{(i_{1} j)} \\ u_{21}^{(i_{1} j)} & u_{22}^{(i_{1} j)} & \cdot \cdot \cdot & u_{2, m_{j}}^{(i_{1} j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ u_{m_{i_{1}}, 1}^{(i_{1} j)} & u_{m_{i_{1}}, 2}^{(i_{1} j)} & \cdot \cdot \cdot & u_{m_{i_{1}}, m_{j}}^{(i_{1} j)} \\ u_{11}^{(i_{2} j)} & u_{12}^{(i_{2} j)} & \cdot \cdot \cdot & u_{1, m_{j}}^{(i_{2} j)} \\ u_{21}^{(i_{2} j)} & u_{22}^{(i_{2} j)} & \cdot \cdot \cdot & u_{2, m_{j}}^{(i_{2} j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ u_{m_{i_{2}}, 1}^{(i_{2}, j)} & u_{m_{i_{2}}, 2}^{(i_{2}, j)} & \cdot \cdot \cdot & u_{m_{i_{2}}, m_{j}}^{(i_{2}, j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \end{matrix}]

对角阵W_j可表示成：

W_{j} = [\begin{matrix} w_{1}^{(j)} \\ w_{2}^{(j)} \\ \cdot \cdot \cdot \\ w_{m_{j}}^{(j)} \end{matrix}]

右矩阵V_j可表示成：

V_{j} = [\begin{matrix} v_{11}^{(j)} & v_{12}^{(j)} & \cdot \cdot \cdot & v_{1, m_{j}}^{(j)} \\ v_{21}^{(j)} & v_{22}^{(j)} & \cdot \cdot \cdot & v_{2, m_{j}}^{(j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ v_{m_{j}, 1}^{(j)} & v_{m_{j}, 2}^{(j)} & \cdot \cdot \cdot & v_{m_{j}, m_{j}}^{(j)} \end{matrix}]

【0065】奇异值分解的更详细的描述见上面的定义部分。应该指出为了增强此方法，可在块215中对矩阵M_j的弗罗比尼斯范数进行计算。计算出的弗罗比尼斯范数称为列矩阵范数。如果弗罗比尼斯范数大于预定的阈值，对矩阵M_j执行奇异值分解以生成如上所述的矩阵U_j，W_j和V_j。然而，如果弗罗比尼斯范数不大于预定的阈值，将V_j设定为等于单位矩阵，将W_j的元设定为0，并且不计算U_j的元。即，由于U_j的元不影响结果，因此该方法减少了计算量。因此，由于某些计算被跳过或被简化，使该方法得以增强。

【0066】在步骤220，舍弃对角阵W_j中小于预定阈值的奇异值，以生成简化的对角阵

在一个实施方式中，通过将奇异值设定为0而将其舍弃。在对角阵W_j中剩余(或保留的)的没有被舍弃的奇异值的数量可称作r_j。这样，简化的对角阵

可表示成：

W_{j}^{(r)} = [\begin{matrix} 0 \\ \cdot \cdot \cdot \\ 0 \\ w_{m_{j} + 1 - r_{j}}^{(j)} \\ \cdot \cdot \cdot \\ w_{m_{j}}^{(j)} \end{matrix}]

【0067】在步骤225，左矩阵U_j乘以简化的对角阵

以生成临时的列矩阵

临时列矩阵的右边r_j列包含非0元，临时列矩阵的剩余元为0。临时列矩阵

可表示成：

M_{j}^{(t)} = [\begin{matrix} 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{1, m_{j} + 1 - r_{j}}^{(jj)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{1, m_{j}}^{(jj)} \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{2, m_{j} + 1 - r_{j}}^{(jj)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{2, m_{j}}^{(jj)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{m_{j}, m_{j} + 1 - r_{j}}^{(jj)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{m_{j}, m_{j}}^{(jj)} \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{1, m_{j} + 1 - r_{j}}^{(i_{1} j)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{1, m_{j}}^{(i_{1} j)} \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{2, m_{j} + 1 - r_{j}}^{(i_{1} j)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{2, m_{j}}^{(i_{1} j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{m_{i_{1}}, m_{j} + 1 - r_{j}}^{(i_{1} j)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{m_{i_{1}}, m_{j}}^{(i_{1} j)} \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{1, m_{j} + 1 - r_{j}}^{(i_{2} j)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{1, m_{j}}^{(i_{2} j)} \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{2, m_{j} + 1 - r_{j}}^{(i_{2} j)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{2, m_{j}}^{(i_{2} j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{m_{i_{2}}, m_{j} + 1 - r_{j}}^{(i_{2} j)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{m_{i_{2}}, m_{j}}^{(i_{2} j)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \end{matrix}] .

【0068】在步骤230，将临时列矩阵

中的系数组合到对应于子矩阵(如，F_ij)的原始集的块中，该子矩阵包含组合到矩阵M_j的块稀疏矩阵A的j列的f系数。即，临时列矩阵

可表示成：

M_{j}^{(t)} = [\begin{matrix} F_{jj}^{(t)} \\ F_{i_{1} j}^{(t)} \\ F_{i_{2} j}^{(t)} \\ \cdot \cdot \cdot \end{matrix}],

其中，临时子矩阵

可表示成：

F_{ij}^{(t)} = [\begin{matrix} 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{1, m_{j} + 1 - r_{j}}^{(ij)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{1, m_{j}}^{(ij)} \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{2, m_{j} + 1 - r_{j}}^{(ij)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{2, m_{j}}^{(ij)} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ 0 & \cdot \cdot \cdot & 0 & w_{m_{j} + 1 - r_{j}}^{(j)} u_{m_{i}, m_{j} + 1 - r_{j}}^{(ij)} & \cdot \cdot \cdot & w_{m_{j}}^{(j)} u_{m_{i}, m_{j}}^{(ij)} \end{matrix}] .

【0069】在步骤235，对块稀疏矩阵A的剩余列重复210到230的步骤，以得到块稀疏矩阵A的剩余的临时列矩阵。

【0070】在步骤240，对每个i，如果V_i不为单位矩阵，临时列矩阵

的临时子矩阵

被右矩阵的转置矩阵

左乘，以生成子矩阵

的变换形式。该乘积可表示成：

{\hat{F}}_{ij}^{(t)} = V_{i}^{T} F_{ij}^{(t)} .

对每个i，如果V_i为单位矩阵，该乘积为

{\hat{F}}_{ij}^{(t)} = F_{ij}^{(t)} .

【0071】在步骤245，对每个i，如果V_i不为单位矩阵，来自子矩阵A_ij的包含

系数的子向量Φ_ij被转置的右矩阵

左乘，以生成子向量Φ_ij的变换形式。该乘积可表示成：

{\hat{Φ}}_{ij} = V_{j}^{T} Φ_{ij} .

对每个i，如果V_i为单位矩阵，该乘积为

{\hat{Φ}}_{ij} = Φ_{ij} .

【0072】在步骤250，对每个j，如果V_j不为单位矩阵，流体分体积向量

被V_j右乘以生成流体分体积向量

的变换形式。该乘积可表示成：

{\overset{\hat{~}}{V}}_{j} = {\tilde{V}}_{j} V_{j} .

对每个j，如果V_j为单位矩阵，该乘积为

{\overset{\hat{~}}{V}}_{j} = {\tilde{V}}_{j} .

【0073】在步骤255，构造具有和块稀疏矩阵A相同块结构和子矩阵形式的变换的块稀疏矩阵

这样，变换的非对角子矩阵

其中i≠j，可表示成：

{\hat{A}}_{ij} = [\begin{matrix} {\hat{F}}_{ij}^{(t)} & {\hat{Φ}}_{ij} \\ 0 & 0 \end{matrix}]

变换的对角子矩阵

可表示成：

{\hat{A}}_{ii} = [\begin{matrix} I + {\hat{F}}_{ii}^{(t)} & {\hat{Φ}}_{ii} \\ {\overset{\hat{~}}{V}}_{i} & c_{i} \end{matrix}]

或

【0074】在步骤260，消去对应于变换的块稀疏矩阵

中舍弃的奇异值的变换的分体积。在一个实施方式中，通过从

减去

乘

的结果，消去对应于变换的对角子矩阵中舍弃的奇异值的变换的分体积

其中k＝1，...，m_j-r_j。对变换的非对角子矩阵

进行类似的操作，即从变换的非对角子矩阵

中底部行减去变换的分体积

和

相乘的结果。对变换的块稀疏矩阵

的所有行重复该消去过程。一旦消去对应于变换的块稀疏矩阵

的舍弃的奇异值的变换的分体积，变换的对角子矩阵可表示成：

变换的非对角子矩阵

其中i≠j，可表示成：

【0075】在步骤265，构造简化的变换块稀疏矩阵在简化的变换块稀疏矩阵

的简化的变换对角子矩阵

中，底部r_i+1行仅在最右边r_i+1列包含非零系数。这些系数被放到更小维数的矩阵中，生成简化的变换对角子矩阵其可以表示成：

【0076】简化的变换块稀疏矩阵

的简化的变换非对角子矩阵仅在底部r_i+1行和最右边r_j+1列内包含系数。这样，简化的变换非对角子矩阵

可以表示成：

【0077】在步骤270，通过将剩余块向量R设定为等于右手侧向量B，而对其进行初始化。即，将初始的剩余块向量R⁽⁰⁾设定为等于右手侧向量B。

【0078】在步骤275，构造变换的剩余块向量变换的剩余块向量

可表示成：

{\hat{R}}_{i} = [\begin{matrix} {\hat{R}}_{M, i} \\ {\hat{R}}_{V, i} \end{matrix}] .

变换的质量平衡剩余子向量

可表示成：

{\hat{R}}_{M, i} = V_{i}^{T} R_{M, i}^{(0)} .

变换的体积约束剩余子向量

可通过计算得到，即从体积约束剩余子向量R_V，i中减去对应于变换的块稀疏矩阵

的舍弃的奇异值的变换的分体积乘

的结果，其中k＝1，...，m_i-r_i-1。

【0079】在步骤280，构造简化的变换剩余块向量块向量

，其表示成：

{\hat{R}}_{i}^{(r)} = [\begin{matrix} {\hat{r}}_{m_{i} + 1 - r_{i}}^{(i)} \\ \cdot \cdot \cdot \\ {\hat{r}}_{m_{i}}^{(i)} \\ {\hat{r}}_{m_{i} + 1}^{(i)} \end{matrix}]

其中，简化的变换剩余块向量块向量

的值由变换的质量平衡剩余子向量

得到。

【0080】在步骤282，求解矩阵方程

{\hat{A}}^{(r)} δ {\hat{X}}^{(r)} = {\hat{R}}^{(r)}

的简化的变换解变化块向量

其包含简化的质量未知数和压力未知数的组合。

【0081】在步骤284，通过将每一子向量设定为如下，创建部分变换的解变化向量

δ {\hat{X}}_{i}^{(0)} = [\begin{matrix} 0 \\ \cdot \cdot \cdot \\ 0 \\ δ {\hat{x}}_{m_{i} + 1 - r_{i}}^{(r, i)} \\ \cdot \cdot \cdot \\ δ {\hat{x}}_{m_{i}}^{(r, i)} \\ δ {\hat{p}}^{(i)} \end{matrix}]

零对应于舍弃的奇异值。

【0082】在步骤286，使用部分换的解变化向量和下述方程，来计算变换的剩余块向量

{\hat{R}}^{(0)} = \hat{R} - \hat{A} δ {\hat{X}}^{(0)}

【0083】在步骤288，变换的解变化向量

的值构造成如下所示：

δ {\hat{X}}_{i} = [\begin{matrix} {\hat{r}}_{1}^{(0, i)} \\ \cdot \cdot \cdot \\ {\hat{r}}_{m_{i} - r_{i}}^{(0, i)} \\ δ {\hat{x}}_{m_{i} + 1 - r_{i}}^{(r, i)} \\ \cdot \cdot \cdot \\ δ {\hat{x}}_{m_{i}}^{(r, i)} \\ δ {\hat{p}}^{(i)} \end{matrix}],

这里变换的解变化向量

的下半区的值从部分变换的解变化向量

的下半区得到，变换的解变化向量

的上半区的值从步骤286中确定的变换的剩余块向量得到。

【0084】在步骤290，根据δM_i＝V_iX_m，i计算质量变化δM_i，其中参考步骤215前面的所述，V_i代表包含右奇异向量的矩阵，并且δX_m，i代表

中仅包含质量未知数元的向量。然后，在步骤292用质量变化δM_i组合如下的解变化块子向量δX_i：

{δX}_{i} = [\begin{matrix} {δM}_{i} \\ δ P_{i} \end{matrix}]

【0085】然后，在步骤294使用解变化子向量δX_i，如下更新解子向量：X_i为X_i＝X_i+δX_i。在步骤296，对X的每一子向量重复从290到294的步骤。以这种方式，计算当前解块向量X的估计值。

【0086】在步骤297，根据R＝B-AX计算更新的剩余块向量R。在步骤298，确定更新的剩余块向量R和解块向量X是否满足预定的停止判据。如果答案为否定，返回步骤275的处理。如果答案为肯定，那么该处理结束。

【0087】图3图解说明计算机网络300，在其中本发明的实施方式可被实施。计算机网300包括系统计算机330，其可通过常规的个人计算机或工作站，如基于UNIX的工作站实现。系统计算机330与可为外部硬盘存储设备的磁盘存储器设备329、331和333相通信。可以认为磁盘存储器设备329、331和333为常规的硬盘驱动器，且这样将可由局域网或远程访问实现。当然，尽管图示的磁盘存储器设备329、331和333为分开的设备，也可按照需要使用单个磁盘存储器设备来存储任何或全部的程序指令、测量数据和结果。

【0088】在一个实施方式中，输入数据存储在磁盘存储器设备331中。系统计算机330可从磁盘存储器设备331重新得到合适的数据，以根据对应于这里描述的方法的程序指令来求解隐式储层仿真矩阵方程。程序指令可用计算机编程语言，如C++，Java等编写。程序指令可以存储在如程序磁盘存储器设备333的计算机可读的存储器中。当然，存储程序指令的存储器介质可为用于存储计算机程序的任何常规类型，包括硬盘驱动器，软盘，CD-ROM或其它光介质，磁带等。

【0089】根据优选实施方式，系统计算机330将输出主要显示在图形显示器327上，或可选为通过打印机328显示。系统计算机230可将上述方法的结果存储在磁盘存储器329上，以便以后使用或进一步分析。键盘326和定点设备(如，鼠标，跟踪球，或类似设备)225可随系统计算机330一起提供以能进行交互操作。

【0090】系统计算机330可位于远离储层的数据中心。尽管图3图解说明的磁盘存储器331直接连接到系统计算机330，但也可认为磁盘存储器设备331通过局域网或远程访问是可存取的。而且，尽管磁盘存储器设备329、331图解说明为用于存储输入数据和分析结果的分开的设备，但磁盘存储器设备329、331可在单个磁盘驱动器(与程序磁盘存储器设备333在一起，或者与其相分离)内实现，或以任何其它的常规方式实现，本领域技术人员在参考本说明书后能完全理解它。

Claims

1.一种用于仿真储层中流体流动的计算机实现方法，所述方法包括：

构造流体流动方程系统；

将所述流体流动方程系统写为矩阵方程AX＝B，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量；

接收所述块稀疏矩阵和所述右手侧块向量；

从所述块稀疏矩阵构造简化的变换块稀疏矩阵；

从所述块稀疏矩阵和所述右手侧块向量构造简化的变换剩余块向量；其中构造所述简化的变换块稀疏矩阵包括：

将一个或更多个质量变化项的一个或更多个系数从所述块稀疏矩阵的列组合到列矩阵中；

对所述列矩阵执行奇异值分解以生成一个左矩阵，一个对角阵和一个右矩阵；

舍弃一个或更多个小于预定门限值的奇异值，以生成简化的对角阵；

使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量求解所述解块向量；

基于所述解块向量显示所述储层的流动速率，其中所述矩阵方程代表一维或更多维流体流动方程的系统，所述流体流动方程具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数以及所述质量变化项的所述一个或更多个系数；以及

基于仿真的流动速率管理来自储层的碳氢化合物的产生。

2.根据权利要求1的所述方法，其中构造所述简化的变换块稀疏矩阵进一步包括构造变换的块稀疏矩阵，其具有与所述块稀疏矩阵相同的块结构和子矩阵形式。

3.根据权利要求1的所述方法，其中所述对角阵包括增序排列的一个或更多个奇异值。

4.根据权利要求1的所述方法，其中构造所述简化的变换块稀疏矩阵包括用所述简化的对角阵乘所述左矩阵，以生成临时列矩阵。

5.根据权利要求1的所述方法，其中构造所述简化的变换块稀疏矩阵包括：

用所述简化对角阵乘所述左矩阵，以生成具有多个质量变化项系数的临时列矩阵，和

将所述临时列矩阵的所述质量变化项系数组合到对应于所述列矩阵的所述质量变化项系数子矩阵的一个或更多个临时质量变化项系数子矩阵。

6.根据权利要求5所述的方法，其中构造所述简化的变换块稀疏矩阵包括：所述右矩阵的所述转置左乘每一临时质量变化项系数子矩阵。

7.根据权利要求6所述的方法，其中构造所述简化的变换块稀疏矩阵包括：所述右矩阵的所述转置左乘在所述块稀疏矩阵中包含所述压力变化项系数的每一子向量，以生成变换的压力变化项系数子向量。

8.根据权利要求7所述的方法，其中构造所述简化的变换块稀疏矩阵包括：所述右矩阵右乘所述块稀疏矩中每一流体分体积子向量，以生成变换的流体分体积向量。

9.根据权利要求8所述的方法，其中构造所述简化的变换块稀疏矩阵包括从所述变换的临时质量变化项系数子矩阵，所述变换的压力变化项系数子向量和所述变换的流体分体积向量中的一个或更多个构造变换的块稀疏矩阵。

10.根据权利要求9所述的方法，其中构造所述简化的变换块稀疏矩阵包括消去与所述变换的块稀疏矩阵中的所述舍弃的奇异值对应的一个或更多个变换的流体分体积，以生成所述简化的变换块稀疏矩阵。

11.根据权利要求1所述的方法，其中构造所述简化的变换块稀疏矩阵包括：

为所述列矩阵确定列矩阵范数；

如果所述列矩阵范数等于或小于所述预定的门限值，跳过所述奇异值分解；

如果所述列矩阵范数等于或小于所述预定的门限值，将所述奇异值设定为0；和

如果所述列矩阵范数等于或小于所述预定的门限值，将所述右矩阵设定为所述单位矩阵。

12.根据权利要求11所述的方法，其中确定所述列矩阵范数包括计算所述列矩阵的弗罗比尼斯范数。

13.根据权利要求1所述的方法，其中所述简化的变换块稀疏矩阵包括一个或更多个简化的变换对角子矩阵和一个或更多个简化的变换非对角子矩阵，其中每一简化的变换对角子矩阵仅在每一变换的对角子矩阵的底部r_i+1行和最右边r_i+1列内包括质量变化项系数和压力变化项系数。

14.根据权利要求13所述的方法，其中每一简化的变换非对角子矩阵仅在每一变换非对角阵的底部r_i+1行和最右边r_j+1列内包括质量变化项系数和压力变化项系数。

15.根据权利要求1所述的方法，其中构造所述简化的变换剩余块向量包括构造变换的剩余块向量。

16.根据权利要求1所述的方法，其中构造所述变换的剩余块向量包括构造具有变换的质量平衡剩余子向量和变换的体积约束剩余子向量的变换的剩余块向量。

17.根据权利要求1所述的方法，其中构造所述简化的变换剩余块向量包括：所述右矩阵的所述转置左乘质量平衡剩余子向量，以生成变换的质量平衡剩余子向量。

18.根据权利要求1所述的方法，其中求解所述解块向量包括：使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量求解简化的变换解变化块向量。

19.根据权利要求1所述的方法，其中求解所述解块向量包括：

使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量，来求解简化的变换解变化块向量；和

将所述简化的变换解变化块向量转换成解变化块向量。

20.根据权利要求1所述的方法，其中求解所述解块向量包括：

将所述解变化块向量加到所述解块向量的当前估计值，以更新所述解块向量。

21.根据权利要求1所述的方法，其中求解所述解块向量包括：

将所述简化的变换解变化块向量转换成具有一个或更多个质量未知值变化和一个或更多个压力未知值变化的解变化块向量。

22.一种用于仿真储层中流体流动的计算机实现方法，所述方法包括：

构造流体流动方程系统；

从所述块稀疏矩阵构造简化的变换块稀疏矩阵；

从所述块稀疏矩阵和所述右手侧块向量构造简化的变换剩余块向量；

使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量，来求解简化的变换解变化块向量；

将所述简化的变换解变化块向量变换成具有一个或更多个质量未知数变化和一个或更多个压力未知数变化的解变化块向量；

将所述解变化块向量加到所述解块向量的当前估计值，以更新所述解块向量；

基于所述解块向量显示结果，其中所述矩阵方程代表一维或更多维流体流动方程的系统，所述流体流动方程具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数以及所述质量变化项的所述一个或更多个系数；以及

基于仿真的流动速率管理来自储层的碳氢化合物的产生。

23.一种用于仿真储层中流体流动的计算机实现方法，所述方法包括：

构造流体流动方程系统；

接收所述块稀疏矩阵和所述第一块向量；

从所述块稀疏矩阵构造简化的变换块稀疏矩阵；

从所述块稀疏矩阵和所述第一块向量构造简化的变换剩余块向量；其中构造所述简化的变换块稀疏矩阵包括：

将与所述至少一个质量变化项相关的至少一个系数，从所述块稀疏矩阵的一列组合到一个列矩阵；

对所述列矩阵执行奇异值分解以生成第一矩阵，第二矩阵和第三矩阵；和

舍弃小于预定门限值的每一奇异值以生成简化的对角阵；

使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量，来求解所述解块向量；

基于仿真的流动速率管理来自储层的碳氢化合物的产生。

24.根据权利要求23所述的方法，其中构造所述简化的变换块稀疏矩阵包括：

从所述列矩阵确定列矩阵范数；

如果所述列矩阵范数等于或小于所述预定门限值，跳过所述奇异值分解；

如果所述列矩阵范数等于或小于所述预定的所述门限值，将所述奇异值设定为0；和

如果所述列矩阵范数等于或小于所述预定的所述门限值，将所述第三矩阵设定为所述单位矩阵。

25.根据权利要求24所述的方法，其中确定所述列矩阵范数包括计算所述列矩阵的所述弗罗比尼斯范数。

26.一种系统，其包括：

计算机；

构造流体流动方程的系统的装置；

将所述流体流动方程系统写为矩阵方程AX＝B的装置，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量；

接收矩阵方程AX＝B的块稀疏矩阵和右手侧块向量的装置，其中A代表所述块稀疏矩阵，B代表所述右手侧块向量，X代表解块向量；

从所述块稀疏矩阵构造简化的变换块稀疏矩阵的装置；

从所述块稀疏矩阵和所述右手侧块向量构造简化的变换剩余块向量的装置；其中为构造所述简化的变换块稀疏矩阵，构造流体流动方程的系统的所述装置包括：

将一个或更多个质量变化项的一个或更多个系数从所述块稀疏矩阵的列组合到列矩阵的装置；

在所述列矩阵执行奇异值分解以生成一个左矩阵，一个对角阵和一个右矩阵的装置；

舍弃一个或更多个小于预定门限值的奇异值以生成简化的对角矩阵的装置；

用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量求解所述解块向量的装置；

耦连到所述计算机的图形显示器，并被配置为将基于所述解块向量的结果显示在所述图形显示器上；以及

基于仿真的流动速率管理来自储层的碳氢化合物的产生的装置。

27.根据权利要求26所述的系统，其中所述矩阵方程代表一维或更多维的流体流动方程系统，其具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数，和所述质量变化项的所述一个或更多个系数。

28.根据权利要求26所述的系统，其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为构造具有与所述块稀疏矩阵相同的块结构和子矩阵的变换的块稀疏矩阵。

29.根据权利要求26所述的系统，其中所述矩阵方程代表一维或更多维的流体流动方程系统，其具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数，和所述质量变化项的所述一个或更多个系数，并且其中所述对角矩阵包括以递增次序排列的一个或更多个奇异值。

30.根据权利要求26所述的系统，其中所述矩阵方程代表一维或更多维的流体流动方程系统，其具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数和所述质量变化项的所述一个或更多个系数，并且其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为用所述简化的对角矩阵乘所述左矩阵，以生成临时列矩阵。

31.根据权利要求26所述的系统，其中所述矩阵方程代表一维或更多维的流体流动方程系统，其具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数和所述质量变化项的所述一个或更多个系数，并且其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为：

用所述简化的对角矩阵乘所述左矩阵，以生成具有多个质量变化项系数的临时列矩阵；和

将所述临时列矩阵的所述质量变化项系数组合到与所述列矩阵的所述质量变化项系数对应的一个或更多个临时质量变化项系数子矩阵。

32.根据权利要求26所述的系统，其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为用所述右矩阵的所述转置左乘每一临时质量变化项系数子矩阵，以生成变换的临时质量变化项系数子矩阵。

33.根据权利要求32所述的系统，其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为用所述右矩阵的所述转置左乘所述块稀疏矩阵中包含所述压力变化项系数的每一子向量，以生成变换的压力变化项系数子向量。

34.根据权利要求33所述的系统，其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为用所述右矩阵右乘所述块稀疏矩阵中的每一流体分体积子向量，以生成变换的流体分体积向量。

35.根据权利要求34所述的系统，其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为从所述变换的临时质量变化项系数子矩阵、所述变换的压力变化项系数子向量和所述变换的流体分体积向量中的一个或多个来构造变换的块稀疏矩阵。

36.根据权利要求35所述的系统，其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为消去对应于所述变换的块稀疏矩阵中所述丢弃的奇异值的一个或更多个变换的流体分体积，以生成所述简化的变换块稀疏矩阵。

37.根据权利要求26所述的系统，其中为了构造所述简化的变换块稀疏矩阵，所述程序指令被配置为：

为所述列矩阵确定列矩阵范数；

38.根据权利要求37所述的系统，其中为了确定所述列矩阵的范数，所述程序指令被配置为计算所述列矩阵的弗罗比尼斯范数。

39.根据权利要求26所述的系统，其中所述简化的变换块稀疏矩阵包括一个或更多个简化的变换对角子矩阵和一个或更多个简化的变换非对角子矩阵，其中每一简化的变换对角子矩阵仅在每一变换的对角子矩阵的所述底部r_i+1行和最右边r_i+1列包含质量变化项系数和压力变化项系数。

40.根据权利要求39所述的系统，其中每一简化的变换非对角子矩阵仅在每一变换的非对角子矩阵的所述底部r_i+1行和最右边r_j+1列包含质量变化项系数和压力变化项系数。

41.根据权利要求26所述的系统，其中为了构造所述简化的变换剩余块向量，所述程序指令被配置为构造变换的剩余块向量。

42.根据权利要求26所述的系统，其中为了构造所述简化的变换剩余块向量，所述程序指令被配置为构造具有变换的质量平衡剩余子向量和变换的体积约束的剩余子向量的变换的剩余块向量。

43.根据权利要求26所述的系统，其中为了构造所述简化的变换剩余块向量，所述程序指令被配置为用所述右矩阵的所述转置左乘质量平衡剩余子向量，以生成变换的质量平衡剩余子向量。

44.根据权利要求26所述的系统，其中为了求解所述解块向量，所述程序指令被配置为使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量来求解简化的变换解变化块向量。

45.根据权利要求26所述的系统，其中为了求解所述解块向量，所述程序指令被配置为：

使用所述简化的变换块稀疏矩阵和所述简化的变换剩余块向量来求解简化的变换解变化块向量；和

将所述简化的变换解变化块向量转换成解变化块向量。

46.根据权利要求26所述的系统，其中为了求解所述解块向量，所述程序指令被配置为将所述解变化块向量加到所述解块向量的当前估计值，以更新所述解块向量。

47.根据权利要求26所述的系统，其中为了求解所述解块向量，所述程序指令被配置为：

将所述简化的变换解变化块向量转换成其具有一个或更多个质量未知数变化和一个或更多个压力未知数变化的解变化块向量。

48.根据权利要求26所述的系统，其包括耦连到所述计算机的打印机，其被配置为向用户提供所述解块向量。

49.一种用于求解储层仿真的矩阵方程AX＝B的计算机实现方法，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量，所述方法包括：

接收所述块稀疏矩阵和所述右手侧块向量；

从所述块稀疏矩阵构造简化的变换块稀疏矩阵；

从所述块稀疏矩阵和所述右手侧块向量构造简化的剩余块向量，其中构造所述简化的变换块稀疏矩阵包括：

将一个或更多个质量变化项的一个或更多个系数从所述块稀疏矩阵的列组合到列矩阵；

对所述列矩阵执行奇异值分解，以生成一个左矩阵，一个对角矩阵和一个右矩阵；和

舍弃一个或更多个小于预定的门限值的奇异值，以生成简化的对角矩阵；

基于所述解块向量显示所述储层仿真的流动速率，其中所述矩阵方程代表一维或更多维流体流动方程的系统，所述流体流动方程具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数以及所述质量变化项的所述一个或更多个系数；和

基于仿真的流动速率管理来自储层的碳氢化合物的产生。

50.一种用于求解代表储层的矩阵方程AX＝B的计算机实现方法，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量，所述方法包括：

从所述块稀疏矩阵构造简化的变换块稀疏矩阵；

将所述简化的变换解变化块向量转换为具有一个或更多个质量未知值变化和一个或更多个压力未知值变化的解变化块向量；

将所述解变化块向量加上所述解块向量的当前估计值，以更新所述解块向量；

基于所述解块向量显示结果，其中所述矩阵方程代表一维或更多维流体流动方程的系统，所述流体流动方程具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数以及所述质量变化项的所述一个或更多个系数；和

基于仿真的流动速率管理来自储层的碳氢化合物的产生。

51.一种用于求解代表储层的矩阵方程AX＝B的计算机实现方法，其中A代表块稀疏矩阵，B代表右手侧块向量，X代表解块向量，所述方法包括：

接收所述块稀疏矩阵和所述第一块向量；

从所述块稀疏矩阵构造简化的变换块稀疏矩阵；

将与所述至少一个质量变化项相关的至少一个系数从所述块稀疏矩阵的列组合到列矩阵；

对所述列矩阵实行奇异值分解，以生成第一矩阵，第二矩阵和第三矩阵；和

丢弃小于预定门限值的每一奇异值，以生成简化的对角矩阵；

基于所述解块向量显示流动速率，其中所述矩阵方程代表一维或更多维流体流动方程的系统，所述流体流动方程具有一个或更多个压力变化项和所述一个或更多个质量变化项，其中所述块稀疏矩阵包含所述压力变化项的一个或更多个系数以及所述质量变化项的所述一个或更多个系数；和

基于仿真的流动速率管理来自储层的碳氢化合物的产生。