通过重复使用地质模型的流动解来尺度粗化储藏层模型

摘要

提供一种模拟物理过程的方法，该过程诸如在多孔介质中的流体流动，该方法通过执行介质中的密网格计算并且在随后的粗网格计算中重复使用密网格的解来进行。对于地下地层中的流体流动，可以使用该方法优化根据地质模型形成的尺度粗化的计算网格。该方法降低了优化网格以模拟物理过程的成本，该物理过程由扩散方程来数学描述。

说明书

【0001】本申请要求2006年7月7日提交序列号为60/819,186的美国 临时申请的优先权。

背景技术

【0002】本发明一般涉及为扩散过程(通过扩散方程来精确描述的过 程)的计算机模拟产生数值模型的领域，该扩散过程诸如在多孔介质 中的流体流动。更具体地，提供一种将地下地层的精密标度的地质模 型尺度粗化(upscale)到适于计算机模拟的一组网格的方法。

技术领域

【0003】在多孔介质中的流体流动的计算机模拟被广泛用于石油工业、 水文地理学和用于修复受污染地下水的环境研究。模拟预测经常对资 产经济评估、油气资产的消耗计划和政府政策有显著影响。

【0004】在自然界被称为“储藏层”的复杂地下结构中发现油气储藏， 诸如石油和天然气。储藏层由各种类型的具有不同物理属性(诸如多 孔性和渗透性)的多孔介质(岩石)构成。这些属性可以在短距离上 有大幅变化。储藏层中的流体流动由物理属性确定。

【0005】地质学家提供的地质属性随机变化的建模技术，允许建模者 使用大量数据创建地表下模型，用覆盖地表下体积的三维“网格”表示 该地表下模型。因为地质模型网格中的大量单元，所以在地质模型比 例下对感兴趣的各种情况执行储藏层模拟是不实际的。同样，复杂的 属性分布使得简单的平均渗透技术变得无用。因此，“尺度粗化”(用 于流动计算的粗网格构成)已经成为储藏层模拟的主要部分。

【0006】汇编描述岩石属性和地质结构的数据是朝向准确模拟储藏层 中流体流动的关键性步骤。汇编数据的微孔模型包括在每个单元中定 义的岩石属性(如多孔率和渗透率)。地质单元形成了储藏层的非重叠 分区。

【0007】地质微孔模型可以包括数以百万计的地质单元来描述储藏层， 所以对于许多感兴趣情形下储藏层流体运动的直接模拟是受成本限制 的。因此，从经济角度考虑，需要将具体的地质模型转换成具有较少 自由度的粗糙模拟模型，以便以可接受的成本执行储藏层模拟。这个 转换被称为“按比例尺度粗化(scaleup)”和“尺度粗化(upscaling)”。 D.Stern(“Practical Aspects of Scaleup of Simulation Models，”J.Pet.Tech， Sept.2005，pp.74-82)和L.J Durlofsky(“Upscaling and Gridding of Fine Scale Geologic Models for Flow Simulation，”文章发表于8^th Int’l Forum on Reservoir Simulation，Stressa，Italy，June，2005)已经公开了“尺度粗 化”的最新观点。(参考 http://ekofisk.stanford.edu/faculty/durlofskypub12.html)。

【0008】尺度粗化包括建立比地质网格粗糙的模拟网格并且将在地质 网格中限定的属性转换到模拟网格。一旦模拟网格被定义，转换地质 属性一般需要计算地质属性的特定平均值来填充模拟网格。对于一些 属性，诸如多孔率，用适当权重来简单平均就已经足够。为了尺度粗 化渗透率，已经证明基于流动的平均方法是最好的方法。Durlofsky (2005)评述了此方法并且最新的基于流动的渗透率尺度粗化的数学 分析是Wu等人给出。(“Analysis of Upscaling Absolute Permeability，” Discrete and Continuous Dynamical System-Series B，Vol.2，No.2，2002)。

【0009】基于流动的尺度粗化需要求解精密标度网格上的单相达西 (Darcy)流动方程。大多数现有方法需要把密网格与粗糙模拟网格对 齐。最近，He描述了一种使用在密网格上限定的流动解来尺度粗化模 拟网格传输率的方法，并且其中密网格不需要与模拟网格对齐。(C.He， “Structured Flow-based Gridding and Upscaling for Reservoir Simulation，” PhD Thesis，Stanford University，Stanford CA，December，2004)。White 和Home提出了一种算法，当密网格标度(fine-grid scale)下的渗透率 不均匀且各向异性时，该算法可以计算传输率的尺度粗化值。 (“Computing Absolute Transmissibility in the Presence of Fine-Scale Heterogeneity，”paper SPE 16011，Ninth SPS Symposium on Reservoir Simulation，Society of Petroleum Engineers，209-220(1987))。

【0010】如Stern(2005)和Durlofsky(2005)所述，成功的尺度粗化 通常需要能够直接捕捉相关的不均匀性的模拟网格。通常需要迭代过 程，其包括建立多个模拟网格来确定“最佳”网格。这个过程称为网 格最优化。建立多个模拟网格需要对地质模型的反复尺度粗化。为了 尺度粗化渗透率，在密网格标度上产生流动解是最耗时和费钱的步骤。 由于其高成本，网格的自动优化不可行；实际上，在实践中也很少手 动改变模拟网格。因此，模拟模型通常不具有最佳准确度，其可能产 生与地质模型不一致的推算。需要一种方法，其允许更快和更低成本 的网格优化。

发明内容

【0011】本发明提供在尺度粗化期间的更快和更低成本的网格优化。 本发明的主要特征是在地质模型上重复使用直接计算的流动解，因而 避免重复尺度粗化过程中的集约计算部分。反复使用这些流动解来尺 度粗化不同模拟网格的渗透率。通过重复使用流动解，可以更有效地 执行对不同模拟网格的尺度粗化。因此，手动改变模拟网格不再是受 耗时因素限制的，并且自动优化网格成为现实。当通过线性偏微分方 程描述物理模型或当物理模型需要基于非线性方程的数学模型时，可 以使用该方法，如在多孔介质中的多相流体流动的情形下。一种计算 机执行的方法，其用以将感兴趣区域的物理属性从精密标度网格尺度 粗化到多个粗网格，其中所述属性的值是已知的，所述属性与所述感 兴趣区域内的扩散过程相关联，所述方法包括：

(a)选择感兴趣区域的体积，所述体积是是所述感兴趣区域的至少一 部分；

(b)将所述体积细分成多个密网格单元，从而为该体积形成精密标度 网格，并且为多个密网格单元中的每一个获得物理属性值；

(c)使用物理属性的精密标度值，求解表示所选体积的精密标度网格 上的扩散过程的扩散方程，从而产生全局解；

(d)保存所述全局解；

(e)将所选体积细分成具有至少一个粗网格单元的第一粗网格，其中 多个密网格单元比该至少一个粗网格单元大；

(f)选择一个粗网格单元，并且确定所述多个密网格单元中的哪些被 整体或部分地包含在所选的粗网格单元中，其中对于被部分包含 的网格单元使用预先选定的标准；

(g)通过为包括在所选粗网格单元内的密网格单元检索并使用全局 解，为所选的粗网格单元计算物理属性的尺度粗化值；

(h)重复步骤(f)-(g)，为至少一个其他粗网格单元计算物理属性 的尺度粗化值，其中所述至少一个其他粗网格单元从感兴趣区域 的选定体积中的所述至少一个粗糙单元中选择；及

(i)为至少又一个粗网格重复步骤(e)-(h)，其中为每个粗网格使 用全局解。

【0012】在本发明的一些实施例中，通过将所选体积细分成可以重叠 两个或多个部分而产生全局解，并且分别求解在每个子体积内的扩散 方程，其中这些解在子体积之间可兼容。此外，本发明方法不必用于 尺度粗化到粗网格，但取而代之的是可以用来尺度粗化到任何描述的 两个或多个粗化体积。

【0013】在本发明的另一个实施例中，提供一种计算机执行的方法， 其用于将地表下区域的物理属性从密网格中的单元处已知的值按比例 缩放到多个不同单元，其中所述属性与该区域中的扩散过程相关联。 该按比例缩放方法包括：(a)选择地表下区域的一个体积，所述体积 是地表下区域的至少一部分；(b)将所述体积细分成精密标度单元， 从而为该体积形成精密标度网格，并且为每个精密标度单元获得物理 属性值；(c)求解所选体积的精密标度网格上的扩散方程从而产生全 局解，其中通过将该体积细分成两个或多个子体积、分别求解每个子 体积内的扩散方程并且在子体积边界匹配所述解而产生所述全局解； (d)将全局解保存在计算机存储器或数据存储器中；(e)在所述两个 或多个子体积中的一个中定义一个不同单元，所述不同单元与精密标 度单元的大小或形状不同；(f)确定哪些密网格单元可以被整体或部分 地包含在该不同单元中，其中对于被部分包含的网格单元使用预先选 定的标准；(g)通过为包括在不同单元内的密网格单元检索并使用全 局解，从而为不同单元计算物理属性的按比例缩放值(scaled value)； (h)重复步骤(e)-(g)，将从计算机存储器或数据存储器中调出的 全局解应用到每个不同单元，从而为地表下区域的所述体积中的至少 又一个不同单元计算物理属性的按比例缩放值。

【0014】在本发明的另一个实施例中，提供一种方法，该方法用于从 地表下地层中提取油气。该方法包括：获得地表下区域的地质模型， 所述模型为精密标度网格提供介质物理属性的离散值，其中该精密标 度网格覆盖构成地表下区域的至少一部分的所选体积；以及获得所述 物理属性的尺度粗化模型以适用于在储藏层模拟程序中使用。尺度粗 化模型通过下面方法而制成：(i)使用物理属性的精密标度值，来求解 表示所选体积上的精密标度网格上的扩散过程的扩散方程，从而产生 一个全局解；(ii)将所述全局解保存在计算机存储器或数据存储器中； (iii)将该体积细分成第一粗网格，所述粗网格具有的单元比精密标度 网格具有的单元少；(iv)选择一个粗网格单元，并且确定哪些密网格 单元被整体或部分地包含在所选粗网格单元内，其中对于被部分包含 的网格单元使用预先选定的标准；(v)通过为包括在所选粗网格单元 内的密网格单元检索并使用所述全局解，从而为所选粗网格单元计算 物理属性的尺度粗化值；(vi)重复步骤(iv)-(v)，为在地表下区域 的所选体积中的所选的其他粗网格单元计算物理属性的尺度粗化值； (vii)为至少又一个粗网格重复步骤(iii)-(vi)，其中将从计算机存 储器或数据存储器中检索出的所述全局解应用到每个粗网格；及(viii) 基于预先确定的网格优化标准来选择优选的粗网格。该提取方法进一 步包括：至少部分基于储藏层模拟从地表下区域提取油气，其中该储 藏层模拟通过使用优选的粗网格上的物理属性的尺度粗化模型而产 生。

【0015】在本发明的又一个实施例中，提供一种物理属性的尺度粗化 方法。该方法包括：(a)计算至少一个方程的精密标度解，该方程描 述所关注介质内的物理扩散过程(the physics of a diffusive process)，其 中为精密标度网格中每个精密标度单元确定精密标度解，并且该精密 标度解被存储在存储器中；(b)为至少部分所关注介质构建粗网格， 其中粗网格包括多个单元；(c)在所关注介质中形成至少一个尺度粗 化的体积，其中尺度粗化的体积用来计算粗网格中所关注介质内的物 理扩散过程；(d)在精密标度解组和至少一个粗网格之间构建映射， 其中映射包括使至少一个精密标度单元与多个粗网格单元中的一个相 关；(e)为与粗网格单元相关的每个精密标度单元从存储器中检索精 密标度解；(f)使用精密标度解为粗网格计算尺度粗化的物理属性；(g) 为所关注的介质的另外部分构建至少一个另外的粗网格，其中另外的 粗网格包括多个单元；(h)使用从步骤(a)得到的精密标度解，为至 少一个另外的粗网格迭代地重复步骤(c)-(f)。

附图说明

【0016】通过参考下面具体描述和附图可以更好地理解本发明和其优 点：

图1显示尺度粗化物理属性的示例性方法的图表；

图2图示说明执行流体计算的精密标度地质模型；

图3A图示说明具有两层的粗糙结构化网格；

图3B显示结构化的、正交的密网格单元和与密网格对齐的粗网格单 元；

图4A-4C图示说明非结构化的尺度粗化的体积；

图4A显示传统的基于单元的尺度粗化的体积，其中所述尺度粗化的体 积是简化的粗网格单元；

图4B显示菱形的尺度粗化的体积，其用于水平连接，该水平连接被设 计成两个单元间的连接模型；

图4C显示垂直连接的尺度粗化的体积，其包围邻近单元之间的连接；

图5A图示说明不规则粗糙的尺度粗化的体积，其被添加在密网格上；

及

图5B图示说明密网格单元的中心可以用来联系密网格单元与具体的 尺度粗化的体积。

具体实施方式

【0017】将联系其优选实施例对本发明进行描述。然而，在某种程度 上，下面详细描述是针对本发明的某一具体实施例或具体使用，这仅 仅是为了直观说明，并不能被解释为对本发明范围的限制。相反，这 是为了覆盖所有包括在本发明思想和范围内的替换方案、修改方案和 等效方案，这些方案由附属权利要求限定。

【0018】在此处公开的本发明是一种新的尺度粗化程序，其允许为多 个粗糙标度网格重复使用精密标度(全局)解。参考图1，该程序包括 下面步骤：

【0019】步骤1：计算一组求解合适方程的精密标度解，该方程描述所 关注区域内的物理扩散过程。线性压力边界条件可以用来产生流动解。 可以如Wu等人(2002)所述，使用其他边界条件。

【0020】步骤2：构建适于每个所关注区域的粗网格。

【0021】步骤3：根据此粗网格，形成尺度粗化体积以计算在粗网格上 具体的所关注物理属性。

【0022】步骤4：在密网格模型和尺度粗化体积之间构建映射。

【0023】步骤5：对于每个尺度粗化的体积，为每个精密标度单元检索 精密标度解，其中所述每个精密标度单元通过映射与尺度粗化体积相 联系。

【0024】步骤6：使用这些精密标度解，为每个尺度粗化体积计算尺度 粗化的属性。

【0025】步骤7：为新的粗网格重复步骤2-6。为执行在新的粗网格上 的尺度粗化，不重新计算精密标度解。而是从存储设备中检索，诸如 从计算机存储器或磁盘。

【0026】在此公开的方法和先前方法之间的一个差异是在不同粗网格 上重复使用精密标度(全局)解的能力。先前方法需要为每个新的粗 网格计算精密标度解。

【0027】在本发明的一个示例性实施例中，描述了穿过多孔介质的达 西流模型内渗透率的尺度粗化。应当理解尽管该方法用于渗透率和流 体流动，但是其还可以应用于扩散方程描述的其他物理过程，其中该 扩散方程是指：

$a\frac{\partial u}{\partial t}+▿·\left(\overrightarrow{b}u\right)+\mathrm{cu}=▿·(\underset{=}{D}▿u)$                (方程1)

其中a，c和是公知的空间和时间函数。该系数的物理意义取决 于方程使用的环境。对于达西流，a与岩石的压缩率和多孔率有关，是渗透率张量。方程1中的变量u是需要从方程中求解的未知量；其对 应于多孔介质流中的压力、饱和度或浓度。

【0028】图2演示步骤1—为在三维内流动穿过多孔介质的流体模型， 计算精密标度解。密网格或地质微孔模型20包括多孔性和渗透性的岩 石属性，该岩石属性应用于每个单元。三维的流动计算结果在21、22、 23中图示说明。对于达西流的三维(“3-D”)模型，需要三种解计算尺 度粗化的渗透率。尽管通常使用数值计算法，诸如使用储藏层模拟领 域中公知的有限差分或有限元法，但是也可以使用任何所需方法来计 算这些流动解。在这种情况下，所关注的区域是整个模型，但是也可 以根据需要将模型分成几个区域，以便于计算。体积的几个部分的解 可以通过使用合适的边界条件被合并从而形成覆盖整个体积的全局 解，即覆盖模型或所关注区域的全局解。然后存储精密标度解以便以 后使用。

【0029】如果密网格是结构化和正交的，则此处公开的方法将会以其 最有效的模式工作。与非结构化网格相比，结构化的网格允许使用更 简单和更有效的方法来处理信息。特别是，在步骤4中使用的映射算 法可以被简化和使其更有效。如果密网格是既结构化又正交的，则可 以利用更简单和更有效的方法来获得实际问题的解。这些思路的重要 性随着密网格的单元(即解)数量而迅速增加，尤其是在3-D应用中。

【0030】如先前所述，如果表示物理过程的密网格的尺寸非常大以至 于其变得不利于计算整个模型的解，则可以将模型分成几个区域，并 分别获得每个区域的解。优选地，选择比粗网格单元的尺寸大许多的 重叠区域和区域尺寸。此选择可以帮助减少边界条件对局部(区域) 解的影响，并且还可以使区域包含较大标度的特性(Wu等人2002)。 粗网格也可以是非结构化网格。

【0031】图3A-3B图示说明步骤2的实施例—为图2中显示的精密标 度地质模型构建粗结构化网格。在图3A中，粗网格30由粗网格单元 31组成，其仅可被称为结构化的和正交化的密网格单元的特殊组的联 合。

【0032】图3B显示结构化、矩形化的密网格35和两个粗网格单元37， 其中两个粗网格单元37与密网格35对齐。对于这种简单情况，存在 有效离散，即双点通量有限容积近似法。此外，由于网格结构，求解 程序会更有效。因为粗网格单元与密网格对齐，粗网格和密网格之间 的映射是微不足道的，不会产生取样误差。

【0033】在步骤3中，为粗网格计算尺度粗化的体积。尺度粗化的体 积是为解决问题所关注的具体体积。对于储藏层模拟而言，这些体积 一般与粗网格单元或连接有关。

【0034】本发明的方法对于结构化网格和非结构化网格同样适用。图 4A-4C图示说明步骤3的示例性实施例—在粗非结构化网格中形成尺 度粗化的体积。图4A显示传统的基于单元的尺度粗化的体积40和其 内的密网格单元41。粗体积40是粗网格单元的近似，其中粗体积40 内的尺度粗化属性是被计算出的。对于除了在垂直维度上具有层状结 构外，具有非结构化或Vornoi区域网格的粗非结构化网格而言，如果 使用有限差分法来获得粗网格的流动解，则优选图4B中的尺度粗化的 体积44和图4C中尺度粗化的体积45。然而，对于通常非结构化网格 的有限差分法，可以使用基于单元的尺度粗化的体积或基于两个邻近 单元的联合的尺度粗化的体积。尺度粗化的体积允许直接计算传递率， 传递率为有限差分法的主要参数。可以使用序列号为No.6,826,520的 美国专利所公开的方法来计算传递率。本领域的技术人员可能知道其 他方法。对于其他的数值离散方法，可能需要不同的尺度粗化的体积。

【0035】图5A-5B图示说明步骤4的示例性实施例—确定与每个尺度 粗化的体积相关的密网格单元的映射。在图5A中，示出的粗尺度粗化 的体积51被附加到密标度网格50上。在图5B中，一种优选的方法被 描述以确定例如密网格单元52是否与粗尺度粗化的体积51相关(即 认为密网格单元52是否被包括在粗尺度粗化的体积51内)。用这种方 法，如果密网格单元中心53位于粗尺度粗化的体积51内，则密网格 单元52与尺度粗化的体积51相关。Durlofsky(2005)和序列号为 No.6,826,520的美国专利讨论了这种部分包含的方法或标准。也可以使 用本领域公知的其他方法。可以使用许多不同的方法构建密网格和粗 网格之间的映射。例如，可以使用本领域公知的计算几何和网格划分 的几何算法。

【0036】序列号为6,106,561的美国专利展示了一种产生网格的合适方 法。也可以使用本领域公知的其他产生网格的方法。关于这个主题有 许多参考文献，诸如Handbook of Grid Generation(J.F.Thompson et al.， CRC Press，1999)。

【0037】如图1中的步骤6的示例，可以考虑多孔介质的渗透率情形， 该渗透率在模拟储油层以便于从储油层中提取油气的过程中非常重 要。在这种情形下，速度和压力梯度是多孔介质中达西流的精密标度 (即全局)解的分量。因此，为步骤1中计算的三维中的每一维从数 据存储器中检索压力梯度和速度(步骤5)。对于达西流动方程的渗透 属性，Wen和Gomez-Hernandez(“Upscaling Hydraulic Conductivity in Heterogeneous Media，”J.Hydrology 183，9-32(1996))已经证实粗网格 渗透率属性可以由下式表示：

$⟨\overrightarrow{v}⟩=-{\underset{=}{K}}^{*}⟨▿P⟩$        (方程2)

其中是尺度粗化的体积中精密标度速度的体积加权平均值，是 精密标度压力梯度的体积加权平均值，是粗标度渗透率。为每个流 动解计算这些平均值。应当注意速度和压力梯度是矢量，渗透率被表 示为张量。这是优选三个不同解的原因；三个解并且每个解三个方程 (每一个对应一个矢量分量)允许计算具有九个分量的粗标度渗透率 张量。

【0038】一旦密网格解可用，则有多种计算粗网格有效属性的方法。 Durlofsky论述了这些方法。优选的方法是使用体积平均法和方程2。

【0039】如果需要一个新的粗网格来提高性能，则在本发明方法中不 需要重新计算精密标度解。如图1所示，必须重新构建密网格的尺度 粗化的体积及其到密网格的映射，并且精密标度解仅仅是在基于新的 粗网格限定的新尺度粗化的体积上重新取样。然后可以将使用新的粗 网格计算的结果与使用第一粗网格计算的结果相比较。根据不同的粗 网格计算的线性结果或单相流动结果可以与基于地质微孔模型的全局 流动解进行比较，以选择优选的粗网格。可以重复这个过程直到找到 最优化的粗网格。之后，根据这些比较优选的粗网格可以用于基于非 线性方程的数学模型，如多孔介质中流动的多相流体的情形中。

示例1

【0040】使用具有1400万个单元的地质模型来计算密网格，其中58 万个单元是激活单元。获得模型中单相达西流的速度和压力的全局解。 使用本发明方法的一个实施例，对粗网格的初始尺度粗化需要60分钟 的计算时间，并且通过检索和重复使用精密标度解的结果，尺度粗化 一个重新网格化的模型仅需要7分钟。相比之下，以前使用的典型方 法均需要125分钟来尺度粗化初始模型和重新网格化的模型。两个粗 网格都具有40,500个激活单元。

示例2

【0041】使用具有750万个单元的地质模型来计算密网格，其中几乎 所有的单元都被激活。获得模型中单相达西流的速度和压力的全局解。 使用本发明的方法，对粗网格的初始尺度粗化需要390分钟的计算时 间，且通过检索和重复使用精密标度解的结果，尺度粗化一个重新网 格化的模型仅需要20分钟。相比之下，使用先前方法均需要150分钟 来尺度粗化初始模型和重新网格化的模型。两个粗网格都具有87,000 个激活单元。

【0042】示例1的模型尺寸在现行实践中比较常遇到。对于任一尺寸 的模型，优化此处公开的方法与目前使用的方法相比，能进一步提高 减少时间和成本的优势。使用该公开的方法，清楚的是，为重新网格 化的解所减少的大量时间可以使一系列手动重新网格化的解或应用自 动重新网格化的解变得可行。

【0043】尽管按照尺度粗化模拟网格来描述本发明，但是应当理解此 处描述的方法同样适用于未形成网格的样品体积组，即这些体积未形 成地表下区域的非重叠分区。可以随机选择样本体积或根据规则式样 选择样本体积。本发明允许更快和以更低成本统计不同的样品体积组。 还应当注意本发明方法不需要样品体积比精密标度网格单元大。本发 明对比精密标度网格单元小的粗网格单元(尺度粗化的体积)同样有 效。

【0044】尽管按照多孔介质中的流体流动来描述本发明，但是应当理 解扩散方程所描述的其他物理现象的模拟也可以使用此处描述的方法 来实践。例如，可以使用本发明方法模拟固体中的热扩散和液体中的 分子扩散。在这些情形下，可以使用上述步骤将类似于渗透率的物理 属性从密网格计算尺度粗化到粗网格计算。

【0045】上述应用是针对本发明具体实施例，用以解释本发明。然而， 对本领域的技术人员来说，显然可以对此处描述的实施例进行许多修 改和变化。所有此类修改和变化都在本发明范围内，如所附权利要求 限定。

参考文献：

1.D.Sterm，“Practical Aspects of Scaleup of Simulation Models，”JPT(9 月，2005)74.

2.L.J.Durlofsky，“Upscaling and Gridding of Fine Scale Geological Models for Flow Simulation，”Proceedings of the 8th International Forum on Reservoir Simulation(6月20-25，2005)，Stresa，Italy.

3.Ph.Renard and G.de Marily，“Calculating Effective Permeability：A Review，”Advances in Water Resources，20(1997)，253-278.

4.X.H.Wen and J.J.Gomez-Hernandez，“Upscaling Hydraulic Conductivity in Heterogeneous Media，”Journal of Hydrology，183(1996)， 9-32.

5.X.H.Wu，Y.Efendiev，and T.Y.Hou，“Analysis of Upscaling Absolute Permeability，”Discrete and Continuous Dyanmical Systems-Series B，2 (2002)，185-204.

6.S.A.Khan and A.G.Dawson，“Method of Upscaling Permeability for Unstructured Grids，”U.S.Patent No.6,826,520，B1(11月30日，2004).

7.C.He，“Structured Flow-based Gridding and Upscaling for Reservoir Simulation，”PH.D.Thesis(2004)，Stanford University，Stanford，CA.

8.M.de Berg，M.van Kreveld，M.Overmars and O.Schwarzkopf， “Computational Geometry：Algorithms and Applications，”Springer，1997.

9.P.G.Ciarlet，“The Finite Element Method for Elliptic Problems，” NorthHolland，1978.

10.K.Aziz and A.Settari，“Petroleum Reservoir Simulation，”Elsevier， 1979.

11.C.L.Farmer.“Simulation Gridding Method and Apparatus Including a Structured Areal Gridder Adapted for Use by a Reservoir Simulator.”U.S. Patent No.6,106,561(8月22日，2000).

12.J.F.Thompson，B.K.Soni，and N.P.Weatherill.“Handbook of Grid Generation.”CRC Press，1999.

13.C.D.White and R.N.Horne，“Computing Absolute Transmissibility in the Presence of Fine-Scale Heterogeneity，”paper SPE 16011，Ninth SPE Symposium on Reservoir Simulation，Society of Petroleum Engineers， 209-220(1987).