CN102667804A - 使用均匀化混合有限元建模地质性质的方法和系统 - Google Patents

Abstract

本发明提供一种储层的烃管理的方法。该方法包括生成储层的模型，该模型包含在非结构化计算网格中多个均匀化的混合有限元。非结构化计算网格可被粗化从而在模型中形成多个更粗糙计算网格。对流扩散地下过程可在最粗糙计算网格上被评估。结果可从最粗糙计算网格转移到最精细计算网格，并且可从模型预测烃储层的性能参数。预测的性能参数可用于储层的烃管理。

Description

使用均匀化混合有限元建模地质性质的方法和系统

相关申请的交叉参考

本申请要求美国临时申请No.61/263,633的权益，该申请提交于2009年11月23日，标题为Method and System for Modeling GeologicProperties Using Homogenized Mixed Finite Elements，其整个内容为全部目的通过引用并入。

技术领域

本技术的示范实施例涉及估算如通过非结构化栅格表示的非均匀地层内对流扩散地下过程的参数的方法和系统。

背景技术

本节意图介绍可和本技术示范实施例关联的技术的各种方面。据信该讨论有助于提供框架从而促进本技术特别方面的更好理解。因此，应理解本节应如此阅读，并且不必需承认为现有技术。

现代社会非常依赖使用烃用于燃料和化学原料。一般在可称为“储层”的地下岩层中发现烃。从储层去除烃依赖岩层的许多物理性质，例如含烃岩石的渗透率、烃流过岩层的能力，以及烃存在的百分比。经常地，数学模型被用来定位烃并优化烃的生产。数学模型使用地下过程的数值模型预测这样的参数：如生产率、最优钻井位置、烃位置等。

地下过程例如流体流动动态、热流和多孔介质中压力分布的数值建模包括求解对流扩散型的数学方程。在许多这样的应用中，输入数据例如渗透率或导热率通过实验观察获得或通过使用某个理论模型推测。这个输入数据可在可被称为“精细地质网格”的高分辨率网格上呈现。然而，对于大多数应用，关于精细地质的信息量超过实际计算能力，使这样的仿真在计算上是花费高的或困难的。结果，大多数计算可仅在具有较低分辨率的网格上执行。较低分辨率网格可被称为“粗糙计算网格”。

精细地质网格和存储计算网格的分辨率之间的失配意味着规程必须被设计为将精细地质网格上原始输入数据的全部或部分转换到粗糙计算网格的分辨率。该规程被称为尺度粗化(up-scaling)。

存在针对具有变化的复杂度的尺度粗化规程的许多不同方法，范围从较简单(并且经常较不准确)的平均化技术到包括具有不同组边界条件的多个局部问题的更复杂(并且计算上昂贵)技术。例如这些的尺度粗化方法已经被证明相当成功。然而，尺度粗化的方法不为在使用粗糙计算模型研究复杂对流扩散过程时存在的尺度粗化解提供数值准确度的先验估计。

已经提出换算来自地下过程的数据的各种根本不同的多尺度(multi-scale)方法从而直接适应精细尺度描述。与尺度粗化相反，多尺度方法以具有原分辨率的完全问题为目标。尺度粗化方法通常是基于通过最大化局部操作来分解感兴趣的长度尺度和时间尺度。在采用混合有限元法的一些途径中，原问题被分解为两个子问题。首先，使用数值Greens函数按照粗糙尺度求解精细尺度，然后，在将精细尺度信息纳入粗糙尺度基础函数之后求解粗糙尺度问题。见于例如T.Arbogast，S.L Btyant，Numerical Subgrid upscaling for waterfloodsimulations，SPE 66375和M.Peszynska，M.F.Wheeler，I.Yotov，Mortarupscaling for multiphase flow in porous media，ComputationalGeosciences，2002，v.6，No.1，pp.73-100。另一途径采用有限元法构造捕捉小尺度的特定基函数。再次，通过粗糙元的边界条件假设实现定位。参见例如T.Hou，X.H.Wu，A multiscale finite element methodfor elliptic problems in composite materials and porous media，J.Comp.Phys.，1997，v.134，pp.169-189；同样见于例如J.Aarnes，S.Krogstad，k.Lie，A hierarchical multiscale method for two-phase flow based uponmixed finite elements and nonuniform coarse grids，Multiscale Modellingand Simulation，v.5，pp.337-363。

近年来，已经开发基于变化原理的新尺度粗化算法，这些变化原理准确并有效捕捉多尺度介质的效应。见于例如S.P.MacLachlan，J.D.Moulton，Multilevel upscaling through variational coarsening，WaterResources Research，v.42，No.2，W02418。提到的全部多尺度技术向原精细尺度问题提供比具有惯例尺度粗化应用的标准技术更准确的解。然而，这些多尺度方法是为结构化的、主要是矩形的网格开发的。非结构化栅格的使用在数值离散化和尺度粗化方法上放置特定约束。见于例如美国专利No.6,826,520。

在许多情况下，感兴趣的域可表示为堆叠在一起的一组不同厚度的层。地质层可沿垂直或倾斜表面断裂并退化，创造所谓的尖灭(pinch-outs)。尖灭被定义为具有零厚度的地质层的部分。地下环境的地质复杂度和准确度需求在为解决地下问题而考虑的数值方法上强加严格约束。另外，大多数实际问题不仅需要准确确定主要变量(例如压力或温度)，同样需要准确确定其通量(能量、流体、热流的流速)。当前，可应用于大多数地下问题的仅两个离散化方法是有限体积法和混合有限元法。

美国专利No.6,823,297公开多尺度有限体积(MSFV)法，从而用起因于多孔介质中单或多相流动的多个空间尺度解决椭圆问题。在其应用中的主要困难是其依赖分层Voronoi(沃罗诺伊)网格的构建，沃罗诺伊网格的构建对任意三维域或具有内部地质特征(例如断层、尖灭等)的域可能是不可能的。构建这种层级的问题不在本专利中考虑，并且可表示其使用的限制。

使地质数据尺度粗化的有希望的数值离散化方法是混合有限元法，其局部质量守恒，在有不均匀介质的情况下准确，并向主要未知量和通量提供准确近似。然而，混合有限元法不可直接应用于在地下应用中普遍的由非结构化多面体栅格覆盖的域。因此，在不规则或非结构化多面体栅格和任意三维域上使地质数据尺度粗化的技术是有用的。

发明内容

本技术的示范实施例提供使用均匀化混合有限元建模地质性质的方法。该方法包括投影储层特征到水平面上从而形成投影，并创造解析投影中期望特征的二维非结构化计算网格。二维非结构化计算网格被投影到边界面上以定义最精细计算网格。生成至少一个更粗糙计算网格，其中更粗糙计算网格包括多个计算单元。多个计算单元的每个包括多个更精细单元。生成关联多个计算单元的每个的多个计算面，其中计算面的每个包含多个更精细面。使第一未知量关联多个计算单元的每个，并使第二未知量关联多个计算面的每个。在最精细计算网格上推导出宏观复合(macro-hybrid)混合有限元离散化。执行迭代粗化规程从而从最精细计算网格转移已知信息到最粗糙计算网格。求解矩阵方程从而为最粗糙计算网格中多个计算单元的每个的每个第一未知量获得值。同样求解矩阵方程从而为最粗糙计算网格中多个计算面的每个的每个第二未知量获得值。执行迭代恢复(restoration)规程从而恢复主未知量的值到多个更精细单元的每个，并恢复次未知量的值到多个更精细面的每个。

投影储层的特征可包括投影尖灭边界、断层线或井位置到水平面中。投影可以是非正交的和/或倾斜的。

多个二维非结构化分层网格的每个可包括正方形、多边形、四边形或三角形或其任何组合。进一步地，多个计算单元的每个可包括箱体、六角体、棱柱体、四面体或棱锥体。

第一未知量可对应储层的物理性质，例如流体压力或温度。第二未知量可对应通量的法向分量。

最精细计算网格可逼近感兴趣层的边界面。物理性质可在最精细计算网格上定义。物理性质可包括渗透率和/或导热率。该方法可包括执行均匀化混合有限元规程以便在棱柱体网格上求解扩散方程。

本技术的另一示范实施例提供使用均匀化混合有限元建模地质性质的系统。该系统可包括处理器和包括数据库的存储介质，该数据库包括储层数据。该系统也包括存储代码的机器可读介质，该代码经配置以引导处理器投影储层特征到水平面上从而形成投影，并创造解析投影中期望特征的二维非结构化计算网格。代码也可经配置而引导处理器投影二维非结构化计算网格到边界面上，以定义最精细计算网格，并生成至少一个更粗糙计算网格，其中更粗糙计算网格包括多个计算单元，并且多个计算单元的每个包含多个更精细单元。代码也可引导处理器生成关联多个计算单元的每个的多个计算面，其中每个计算面包含多个更精细面。代码也可引导处理器使第一未知量与多个计算单元的每个关联，并引导第二未知量与多个计算面的每个关联，在最精细计算网格上推导出宏观复合混合有限元离散化，并通过粗化规程迭代从而从最精细计算网格转移已知信息到最粗糙计算网格。代码可引导处理器求解矩阵方程从而为最粗糙计算网格中多个计算单元的每个的每个第一未知量获得值，求解矩阵方程从而为最粗糙计算网格中多个计算面每个的第二未知量的每个获得值，并通过恢复规程迭代从而恢复主未知量的值到多个更精细单元的每个，并恢复次未知量的值到多个更精细面的每个。该系统也可包括显示器，其中机器可读介质包括经配置在显示器上生成储层图像的代码。储层数据可包括静毛比(net-to-gross ratio)、孔隙率、渗透率、地震数据、AVA参数、AVO参数或其任何组合。

本技术的另一示范实施例提供储层的烃管理的方法。该方法包括生成在非结构化计算网格中包含多个均匀化混合有限元的储层模型，并粗化非结构化计算网格从而在模型中形成多个更粗糙计算网格。在最粗糙计算网格上评估对流扩散地下过程，并且结果从最粗糙计算网格转移到最精细计算网格。从模型预测烃储层的性能参数，并且预测的性能参数被用于储层的烃管理。

该方法可包括投影储层特征到水平面上从而形成投影，并创造解析投影中期望特征的二维非结构化计算网格。二维非结构化计算网格可投影到边界面上以定义最精细计算网格。可生成至少一个更粗糙计算网格，其中更粗糙计算网格包括多个计算单元，并且多个计算单元的每个包含多个更精细单元。多个计算面关联多个计算单元的每个，其中每个计算面包含多个更精细面。第一未知量可关联多个计算单元的每个，并且第二未知量可关联多个计算面的每个。可在第一计算网格上推导出宏观复合混合有限元离散化，并且可执行迭代粗化规程从而从最精细计算网格转移已知信息到最粗糙计算网格。可求解矩阵方程从而为最粗糙计算网格中多个计算单元的每个的第一未知量的每个获得值。同样可求解矩阵方程从而为最粗糙计算网格中多个计算面的每个的第二未知量的每个获得值。可执行迭代恢复规程从而恢复主未知量的值到多个更精细单元的每个，并恢复次未知量的值到多个更精细面的每个。

储层的烃管理可包括例如烃采掘、烃生产、烃勘探、鉴别潜在烃资源、鉴别井位置、确定井注入率、确定井采掘率、鉴别储层连通度，或其任何组合。性能参数可包括例如生产率、压力、温度、渗透率、传递率、孔隙率、烃成分或其任何组合。

另一示范实施例提供有形的计算机可读介质，其包括经配置以引导计算机执行涉及粗化模型的各种操作的代码。该代码可经配置以投影储层特征到水平面上从而形成投影，并创造解析投影中期望特征的二维非结构化计算网格。代码也可经配置而投影二维非结构化计算网格到边界面上，以定义逼近边界面的最精细计算网格，并生成至少一个更粗糙计算网格，其中更粗糙计算网格包含多个计算单元，并且多个计算单元的每个包含多个更精细单元。代码也可经配置以生成关联多个计算单元的每个的多个计算面，其中每个计算面包含多个更精细面，并使第一未知量与多个计算单元的每个关联，并使第二未知量与多个计算面的每个关联、代码也可经配置而在最精细计算网格上推导出宏观复合混合有限元离散化，并通过粗化规程迭代从而从最精细计算网格转移已知信息到最粗糙计算网格，并求解矩阵方程从而为最粗糙计算网格中多个计算单元每个的第一未知量的每个获得值。代码也可经配置以求解矩阵方程从而为最粗糙计算网格中多个计算面每个的第二未知量的每个获得值，并通过恢复规程迭代从而恢复主未知量的值到多个更精细单元的每个，并恢复次未知量的值到多个更精细面的每个。代码也可经配置以引导处理器显示储层的表示。

附图说明

通过参考下面详细描述和附图可更好地理解本技术的优点，其中：

图1是过程流程图，其根据本技术的示范实施例示出在非结构化计算网格上粗化地质模型的方法；

图2是示范储层的顶视图，其根据本技术的示范实施例示出在储层上最精细计算网格的平面投影；

图3是示范储层的顶视图，其根据本技术的示范实施例图解说明计算网格粗化的第一级的平面投影；

图4是示范储层的顶视图，其根据本技术的示范实施例示出计算网格粗化的另一级的平面投影；

图5是示范储层的顶视图，其根据本技术的示范实施例示出计算网格粗化的另一级的平面投影；

图6是示范储层的顶视图，其根据本技术的示范实施例示出计算网格粗化的最终级的平面投影以创建最粗糙的计算网格；

图7是示范储层的透视图，其根据本技术的示范实施例示出计算网格到层的边界面上的垂直投影；

图8是储层计算域的透视图，其根据本技术的示范实施例图解说明地质层之间的界面；

图9是图示，其根据本技术的示范实施例示出分割为子域(Ω_i，i＝1，...，10)的域(Ω)的二维表示；

图10A和10B是示意图，其根据本技术的示范实施例示出两个粗糙计算网格单元(E∈Ω_H)分割为多个精细计算网格单元(e∈Ω_h)；

图11A和11B是示意图，其根据本技术的示范实施例示出垂直四边形面分割为子面；

图12是示意图，其根据本技术的示范实施例示出垂直三角面F分割为子面F_l，l＝1，...，4；

图13是示意图，其根据本技术的实施例示出粗糙棱柱体划分为四个精细棱柱体1302；以及

图14是计算机系统的框图，在该计算机系统上可实现执行本技术的实施例的处理操作的软件。

具体实施方式

在下面详细描述部分中，连同优选实施例描述本技术的特定实施例。然而，就下面描述针对本技术的特别实施例或特别用途来说，其意图仅用于示范目的并简单提供示范实施例的描述。因此，本技术不限于下面描述的特定实施例，相反，这样的技术包括落入附属权利要求的真实精神和保护范围内的全部替换、修改和等效物。

起初，并为了容易参考，阐述用于本申请中使用的某些术语及其用于上下文中的含义。就本文使用的术语不在下面定义来说，应给予其如在至少一个印刷出版物或已授权专利中反映的本领域技术人员给予的最广泛定义。进一步地，由于服务于相同或相似目的的全部等效物、同义词、新发展和术语或技术被认为在本权利要求的保护范围内，因此本技术不受下面示出的术语的用法限制。

“粗化”指通过使单元更大，例如表示储层中的更大空间，减少仿真模型中单元的数目。可执行粗化的过程被称为“放大”。粗化经常被用来通过在生成或运行仿真模型之前减少地质模型中的单元数目而降低计算成本。

“共同尺度模型”指其中地质模型尺度相似于仿真模型尺度的条件。在此情况下，地质模型的粗化不在仿真之前执行。

“计算机可读介质”或“有形机器可读介质”如在此使用的，指参与提供指令到处理器以便执行的任何有形存储介质。这样的介质可包括但不限于非易失性介质和易失性介质。非易失性介质包括例如非易失随机存取存储器(NVRAM)或磁盘或光盘。易失性介质包括动态存储器例如主存储器。计算机可读介质的普通形式包括例如软盘、软磁盘、硬盘、硬盘阵列、磁带，或任何其它磁介质、磁光介质、只读光盘(CD-ROM)、任何其它光介质、随机存取存储器(RAM)、可编程只读存储器(PROM)、EPROM、FLASH-EPROM、固态介质如存储卡、任何其它存储芯片或卡式盒，或计算机可从其读取数据或指令的任何其它有形介质。在计算机可读介质被配置为数据库时，要理解数据库可以是任何类型的数据库，例如关系数据库、分层数据库、面向对象的数据库，等等。

“示范”在此专门用来意指“用作例子、实例或说明”。在此描述为“示范”的任何实施例不解释为优于其它实施例或比其它实施例有优势。

“烃管理”包括烃采掘、烃生产、烃勘探、鉴别潜在烃资源、鉴别井位置、确定井注入和/或采掘率、鉴别储层连通度，烃资源的获取、处置和/放弃，在烃管理决策之前的回顾，以及任何其它烃相关行为或活动。

“渗透率”是岩石传输流体通过岩石的互连孔隙空间的能力。可使用达两定律测量渗透率：Q＝(kΔPA)/(μL)，其中Q＝流速(cm³/s)，ΔP＝跨长度为L(cm)和横截面积为A(cm²)的圆柱体的压降(atm)，μ＝流体粘度(cp)，并且k＝渗透率(达两)。渗透率测量值的习惯单位是毫达两。术语“相对可渗透”关于地层或其部分被定义为10毫达两或更多(例如10或100毫达两)的平均渗透率。术语“相对低渗透率”关于地层或其部分被定义为低于约10毫达两的平均渗透率。不渗透层一般具有低于约0.1毫达两的渗透率。

“孔隙体积”或“孔隙度”被定义为以百分比表达的孔隙空间体积对材料总体积的比率。孔隙度是储层岩石的流体存储容量的度量。总或绝对孔隙度包括全部孔隙空间，而有效孔隙度仅包括互连孔隙并对应可用于消耗的孔隙体积。

向量函数的“最低阶RT(Raviart-Thomas)”有限元空间基于计算空间的分割，例如分割为四面体。参见例如F.Brezzi和M.Fortin，Mixed and hybrid finite element methods，1991，或在：Handbook ofNumerical Analysis，vol.2，1991，pp.523-639中的J.E.Roberts和J.M.Thomsa，Mixed and hybrid Methods。

“地质模型”是地下土方量，例如石油储层或沉积盆地的基于计算机的表示。地质模型可采取许多不同形式。取决于环境，为石油应用建立的描述性或静态地质模型可以是被分配地质和/或地球物理学性质例如岩石学、孔隙度、声阻抗、渗透率或水饱和(这样的性质在此总称为“储层性质”)的单元的3D阵列形式。许多地质模型受地层面或结构面(例如洪泛面、层序界面、流体接触面、断层)和边界(例如相变)的约束。这些表面和边界在可能具有不同储层性质的模型内定义区域。

“储层仿真模型”或“仿真模型”指代实际烃储层的特定数学表示，其可被认为是特殊类型的地质模型。仿真模型用来进行目标是确定最有利润的操作策略的关于油田的未来性能的数值试验。管理烃储层的工程师可创造可能具有变化复杂度的许多不同仿真模型，以便量化储层的过去性能并预测其未来性能。

“放大”指例如通过使性质在某范围内平均化，或通过使用更少数目的测量或计算性质值的点，使生产数据的计算网格聚合到更粗糙计算网格中的过程。该规程降低制作储层模型的计算成本。

“传递率”指给定压降时的在单位粘度的两个点之间的体积流率。传递率是连通度的有用度量。在储层中任何两个分隔(断块或地质带)之间，或在井和储层(或特别的地质带)之间，或在注入器和生产井之间的传递率都可对理解储层中的连通度有用。

概述

本技术的示范实施例公开估算不均匀地层内对流扩散地下过程的参数的方法，该地层被表示为一组堆叠在一起的，并由具有分层组织结构的非结构化栅格覆盖的不同厚度层。这些技术为任意多面体栅格上的扩散型方程利用混合有限元法，参见Yu.Kuznetsov和S.Repin，New mixed finite element method on polygonal and polyhedral meshes，Russ.J.Numer.Anal.Math.Modelling，2003，v.18，pp.261-278(其为使用混合有限元建模这样的过程提供背景知识)。同样参见O.Boiarkine，V.Gvozdev，Yu.Kuznetsov和S.Maliassov，Homogenizedmixed finite element method for diffusion equations on prismatic meshes，Russ.J.Numer.Anal.Math.Modelling，2008，v.23，N5，pp.423-454，以及K.Lipnikov，J.D.Moulton，D.Svyatskiy，A multilevel multiscalemimetic method for two-phase flows in porous media，Journal ofComputational Physics 2008，v.227，pp.6727-6753(其提供应用称为模拟离散化的不同离散化方法的技术，该模拟离散化在许多方面类似于多尺度环境中的混合有限元法，并因此提供用于地质性质的尺度粗化的技术)。

放大问题可以在分层组织的多边形栅格(在下文中称为“计算网格”)的总体上考虑。使用各种规程，信息可按照层次从最精细计算网格系统转移到最粗糙计算网格。代数方程系统可然后在最粗糙计算网格上被求解，由此减少运算的计算需求。使用粗化规程的相反规程，将与最粗糙计算网格上的解有关的信息传播回到(原始)最精细计算网格。对于单棱柱计算网格的情况，地质应用中热传输方程的准确建模的这样方法的方法学和实施详情在专利申请No.PCT/US2008/080515中描述，该申请提交于2008年10月20日，并且标题为“Modeling Subsurface Processes on Unstructured Grid”。

图1是过程流程图，其根据本技术的示范实施例图解在非结构化计算网格上粗化地质模型的方法。该方法一般由参考号100指代。该方法在框102开始，其中地质和几何特征例如尖灭边界、断层线或井位置投影到水平面中。在示范实施例中投影正交执行。在其它实施例中投影可以是非正交的或倾斜的。

如在框104指出的，可创造二维非结构化计算网格从而在该平面上解析期望特征。在其它实施例中，这可以是二维非结构化计算网格的分层顺序。计算网格可由矩形、多边形、四边形或三角形组成。在示范实施例中，生成精细矩形一致性网格，从而覆盖投影域的全部特征。矩形一致性网格可与在其上提供材料数据的最精细计算网格相同大小。

如在框106指出的，二维计算网格(或计算网格的层次)可被投影回到层的边界表面上，从而构造棱柱形计算网格。计算网格含有单元，该单元可包括例如箱体、六角体、棱柱体、四面体、棱锥体和其它任何三维立方体及其组合。因此，这种方式建立的最精细计算网格逼近层的边界表面，并定义感兴趣的最精细计算网格。因此，在该网格上定义物理性质例如渗透率或导热率。

如在框108指出的，可生成至少一个更粗糙计算网格(或构造的计算网格的层次结构)从而获得自嵌入、逻辑连接的更粗糙计算网格的总体。在粗糙计算网格上的每个单元(或计算体积)可称为宏单元并包括更精细计算网格上单元的总体。在示范实施例中，粗化不均匀执行，从而保持精细三角测量接近一些地质或几何特征，但获得远离这些特征的更粗糙分辨率。

如在框110指出的，可创造计算面的层次并将其和计算单元的层次关联。可称为宏面的在粗糙计算网格上的每个计算面由更精细计算网格上的(微)面的总体组成。

在框112，第一未知量可关联每个计算单元，第一未知量被认为位于该单元的中心。第一未知量一般表示单元的物理性质，例如除了别的之外，压力、温度或烃含量。第二未知量可关联每个单元的每个面，并被认为位于面中心。第二未知量表示单元之间通量的法向分量，例如跨单元面的热或质量流。然后可为最精细网格推导出宏观复合混合有限元离散化，如通过框114指出的。在框116，递归粗化/均匀化规程可在通量有限元向量函数的法向分量上使用，从而从最精细计算网格向最粗糙计算网格转移已知信息和物理性质。下面关于图2-10更详细讨论粗化规程。

空间离散化在最粗糙计算网格上产生稀疏矩阵方程，该方程可被称为“尺度粗化的”方程。在框118，可在最粗糙计算网格上为第一和第二未知量求解稀疏矩阵方程。在框120，在最粗糙计算网格上计算的解可然后用于递归规程，从而将解函数和通量向量函数的值恢复到作为宏单元的组成部分的更精细计算网格。继续迭代直到达到最精细计算网格。这将解从最粗糙计算网格有效转移到最精细计算网格。

图2是示范储层的顶视图，其根据本技术的示范实施例示出在储层之上的最精细计算网格的平面投影。储层网格的投影一般由数字200指代。如在图2中示出的，投影是二维计算网格，该网格可以是在储层上叠加的均匀三角形栅格。在地下地层中的对流扩散型问题中，输入数据可关联最精细计算网格的节点(单元交点)或单元。地质和几何特征例如尖灭边界、断层线202或井位置204可例如使用正交投影来投影到水平面中。

图3是储层的顶视图，其根据本技术的示范实施例图解二维计算网格粗化的第一级。如在该图中示出，粗化在具有显著特征的区域302中，例如井的投影202和断层的投影204中可以是不均匀的。尽管二维计算网格被图解为三角形的栅格，但可使用任何数目的其它形状，包括正方形、矩形和其它类型的多面体。

图4是储层的顶视图，其根据本技术的示范实施例示出计算网格粗化的另一级。如关于图3讨论，最精细计算网格可在显著特征例如井202和断层204附近保留。

图5是储层的顶视图，其根据本技术的示范实施例示出创造更粗糙计算网格的粗化的另一级。图6是储层的顶视图，其根据本技术的示范实施例示出创造最粗糙计算网格的粗化的最终级。尽管图2-6示出粗化的连续级，但该方法可应用于计算网格的任意层次序列，例如应用于图2、4和6中计算网格的序列。

图7是示范储层的透视图，其根据本技术的示范实施例示出计算网格到层的边界表面上的垂直投影。储层一般由参考号700表示。如在图7中示出的，投影构造具有单元的棱柱形计算网格，该网格可以是三棱柱、四面体、棱锥体、六角体、箱体或任何其它三维多面实体。这样建立的非结构化棱柱形计算网格逼近全部层的边界表面。一旦构建棱柱形计算网格，那么其可被递归粗化从而生成一系列更粗糙的棱柱形网格。每个粗糙棱柱形网格表示感兴趣的原始物理域，尽管其含有少于最精细棱柱形网格的信息。

问题公式

如果G是具有规则成形边界

(即分段平滑并且在段之间的角度大于0)的R²中的域，那么计算域Ω可被定义如下：

Ω＝{(x，y，z)∈R³：(x，y)∈G，Z_min(x，y)≤z≤Z_max(x，y)}    方程1

其中Z_min(x，y)和Z_max(x，y)是平滑表面。设N_z为正整数并且z＝Z_i(x，y)，i＝0，...，N₂，是在G上定义的单值连续函数，以使：

在

中Z₀(x，y)≡Z_min(x，y)

在

中Z_i-1(x，y)≤Z_i(x，y)

方程2

在

中Z_Nz(x，y)≡Z_max(x，y)

图8是示范储层的计算域透视图，其根据本技术的示范实施例图解地质层之间的界面。计算域一般由参考号800指代。方程1和2可用来定义地质层之间的界面802。即，计算域(Ω)800可被分成为N_z个子域804(条带或层)，对于全部i＝1，...，N_z，其被定义如下：

Ω_i＝{(x，y，z)∈Ω：(x，y)∈G，Z_i-1(x，y)≤z≤Z_i(x，y)}    方程3

可假设子域Ω_i 804满足锥条件，即子域804的边界没有奇点(零度角，等等)，并且另外全部集合：

$G_{i,i-1}=\{(x,y)\∈\stackrel{\‾}{G}:Z_{i-1}(x,y)=Z_i(x,y),(x,y)\∈G\}$ 方程4

没有多边形或由有限数目多边形构成。

图9是根据本技术的示范实施例的分割为子域902(Ω_i，i＝1，...，10)的示范域(Ω)900的垂直剖面的二维表示。子域Ω_i-1和Ω_i之间的界面(或表面)904可由I_i-1，i表示，并且集合

P_i-1，i＝{(x，y，z)：Z＝Z_i-1(x，y)＝Z_i(x，y)，(x，y)∈G}    方程5

可被描述为“尖灭”906，i＝1，...，N_z。如在此使用的，尖灭906，P_i-1，i可具有与P_i-2，i-1或与P_i，i+1，或与该两者的非零交点。对应集合P_i-1，i的边界可由

表示。为简化讨论，可假设尖灭(P_i-1，i)906是简单连接的集合。

粗糙棱柱形网格的定义

如果G_H是G中一致的粗糙三角形计算网格，即G_H中任两个不同三角形具有共同边缘、共同顶点或不相互接触，那么一组连续分段线性表面904可在Ω900中被定义为：

z=Z_H，k(x，y)，             方程6

其中Z_H，k(x，y)是定义表面904的单值函数，k＝0，...，K，并且K是正整数。可假设顶表面908(Z_H，0(x，y))和底表面910(Z_H，K(x，y))与原始域的顶部边界912(Z_min(x，y))和底部边界914(Z_max(x，y))一致或重合。即：

在G中Z_H，0(x，y))＝Z_min(x，y)，Z_H，K(x，y)＝Z_max(x，y)

方程7

以及

在G中Z_H，k-1(x，y)≤Z_H，k(x，y)，k＝1，...，K    方程9

两个主要假设可然后施加于表面904{Z_H，k}的集合上。第一，假设表面Z_H，k不穿过表面Z_i，即对于任何整数k，1≤k≤K，存在整数i，1≤i≤N_z，以使对于源自G的全部(x，y)

Z_i-1(x，y)≤Z_H，k(x，y)≤Z_i(x，y)              方程10

第二，假设如果表面904 Z_H，k，1≤k≤K，满足第一假设的不等式，那么任何两个邻近表面904 Z_H，k-1和Z_H，k，1≤k≤K，除尖灭906 P_i-1，i，1≤i≤N_z之外不创造共同的尖灭906。即：

Z_H，k-1(x，y)≤Z_H，k(x，y)                方程11

对于源自G\G_i-1，i的全部(x，y)。术语Z_H，k，0≤k≤K可用来表示“侧向”粗糙网格表面904。

Ω900中的粗糙计算网格Ω_H由粗糙网格表面z＝Z_H，k(x，y)，0≤k≤K和无限棱柱体组{E_G×(-∞，+∞)}的交点定义，其中E_G是G_H中的特殊三角形。Ω_H是一致的并由粗糙计算网格单元(宏单元)E构成。表面z＝Z_i(x，y)和z＝Z_H，k(x，y)可假设为对于G_H中每个单元E_G是平坦的。在做出这些假设之后，每个计算网格单元E∈Ω_H是具有两个“侧向”和三个垂直非零面的“垂直”棱柱体，或在存在一个或两个零垂直边缘时是退化的“垂直”棱柱体。退化的计算网格单元是棱锥体(一个垂直边缘是零)或四面体(两个垂直边缘是零)。为简化运算，表面Z_i，0≤i≤N_z和Z_H，k，0≤k≤K可针对G_H中每个计算网格单元E_G假设为“几乎平坦”。因此，它们可通过对于G_H中每个E_G平坦的表面以合理准确度近似。

精细棱柱网格的定义

至于粗糙计算网格，精细计算网格可在连续分段线性表面集合的帮助下在900中定义：

z＝Z_h，j(x，y)，    方程12

其中Z_h，j(x，y)是定义表面的单值函数，j＝0，...，J并且J是正整数。可假设顶表面908(Z_h，0(x，y))和底表面910(Z_h，J(x，y))分别与原始域的顶部边界912(Z_min(x，y))和底部边界914(Z_max(x，y))重合。因此，即：

在G中Z_h，0(x，y)＝Z_min(x，y)，Z_h，J(x，y)＝Z_max(x，y)

                                        方程13

以及

在G中Z_h，j-1(x，y)≤Z_h，j(x，y)，j＝1，...，J    方程14

也可假设全部表面{Z_i}和{Z_H，k}属于表面{Z_h，j}的集合。因此，表面{Z_h，j}，j＝1，...，J不穿过来自集合{Z_i}和{Z_H，k}的任何表面。然后，可假设任何两个邻近表面Z_h，j-1和Z_h，j，1≤j≤J属于相同的层Ω_j(例如层902)，仅在点(x，y)∈G_i-1，i中重合，即

Z_h，j-1(x，y)≤Z_h，j(x，y)               方程15

对于源自G\G_i-1，i的全部(x，y)。

G_h可被考虑为G中的一致三角形网格，以使每个三角形E_G∈G_H是G_h中三角形的联合或并集。换句话说，G_h是粗糙网格G_H的三角形精细化。因此，Ω中的精细网络Ω_h可然后由精细网格表面z＝Z_h，j(x，y)，j＝0，...，J和无限垂直棱柱体的集合{e_G×(-∞，+∞)}的交点定义，其中e_G是G_h中的特别三角形。Ω_h是一致的并由精细计算网格单元(微单元)e构成。可然后假设表面z＝Z_h，j(x，y)，(x，y)∈e_G对于每个三角形e_G∈G_h是平坦的。因此，每个网格单元e∈Ω_h是具有两个“侧向”面和三个垂直面的“垂直”棱柱体，或四边形棱锥体，或四面体。根据精细网格Ω_h的定义，每个粗糙单元E∈Ω_H是精细网格单元e∈Ω_h的并集。

图10A和10B是图解，其根据本技术的示范实施例示出两个粗糙计算网格单元(E∈Ω_H)分割为多个精细计算网格单元(e∈Ω_h)。如在图10A中所示，棱柱体1002(E∈Ω_H)可被分割为八个精细棱柱体1004(e∈Ω_h)。图10B图解棱锥体1006(E∈Ω_H)分割为六个精细棱柱体1008和两个精细棱锥体1010(e∈Ω_h)。粗线1012指示E与Ω_i-1和Ω_i之间界面边界I_i-1，i的交点。

微分方程

本技术的示范实施例可用于表示对流扩散过程的粗糙化计算网格。为简便，该规程可通过使用示范的3D扩散型方程来描述：

在Ω中 $-\▿\·(K\▿p)+c\·p=f$ 方程16

其中p是未知函数(其可以是例如计算单元中的压力)，K＝K(x)是扩散张量，c是非负函数，f是源函数，并且

是有界计算域。可假设K是均匀正定矩阵，而且域Ω的边界

被分割为两个非重叠的集合Γ_D和Γ_N。

方程16用方程17中示出的边界条件补充：

在Γ_D上p＝g_D

在Γ_N上 $(K\▿p)\·n+\mathrm{\σ}\·p=g_N$ 方程17

其中n是到Γ_N的向外单位法向向量，σ是非负函数，并且g_D和g_N是给定函数。可假设在方程16和17中表达的问题具有唯一解。

在方程16-17中的微分问题可由等效一阶系统替换：

在Ω中 $u+K\▿p=0$

在Ω’中 $\▿\·u+c\·p=f$ 方程18

在Γ_D上p＝g_D

在Γ_N上-u·n+σ·p＝g_N    方程19

因此，由方程18和19说明的问题可被描述为由方程16和17说明的问题的混合公式。注意以这种方式主要未知量p及其通量u被同时近似。

变分混合公式

在方程18和19中说明的微分问题的变分混合公式可表述如下：寻找

p∈L₂(Ω)，以及λ∈L₂(Γ_N)，以使

$\underset{\mathrm{\Ω}}{\∫}\left(K^{-1}u\right)\·\mathrm{vdx}-\underset{\mathrm{\Ω}}{\∫}p(\▿\·v)\mathrm{dx}+\underset{{\mathrm{\Γ}}_N}{\∫}\mathrm{\λ}(v\·n)\mathrm{ds}=-\underset{{\mathrm{\Γ}}_D}{\∫}{g}_D(v\·n)\mathrm{ds}$

$-\underset{\mathrm{\Ω}}{\∫}(\▿\·u)\mathrm{qdx}-\underset{\mathrm{\Ω}}{\∫}c\·\mathrm{pqdx}=-\underset{\mathrm{\Ω}}{\∫}\mathrm{fqdx}$

$\underset{{\mathrm{\Γ}}_N}{\∫}(u\·n)\mathrm{\μds}-\underset{{\mathrm{\Γ}}_N}{\∫}\mathrm{\σ\λ\μds}=\underset{{\mathrm{\Γ}}_N}{\∫}{g}_N\mathrm{\μds},$ 方程20

对于全部

q∈L₂(Ω)以及μ∈L₂(Γ_N)，其中

${\hat{H}}_{\mathrm{div}}\left(\mathrm{\&Omega;}\right)=\{v:v\&Element;{\left[L_2\left(\mathrm{\&Omega;}\right)\right]}^3,\&dtri;\&CenterDot;v\&Element;L_2\left(\mathrm{\&Omega;}\right),\underset{\&PartialD;\mathrm{\&Omega;}}{\&Integral;}{|v\&CenterDot;n|}^2\mathrm{ds}<\&infin;\}.$

在该公式中，λ是压力函数p＝p(x)到Γ_N上的约束，并且边界条件是自然的。

在Γ_N上σ＝0的情况下，变分混合公式可用下面形式表述：在Γ_N上寻找

u·n＝-g_N并且p∈L₂(Ω)，以使

$\underset{\mathrm{\&Omega;}}{\&Integral;}\left(K^{-1}u\right)\&CenterDot;\mathrm{vdx}-\underset{\mathrm{\&Omega;}}{\&Integral;}p(\&dtri;\&CenterDot;v)\mathrm{dx}=-\underset{{\mathrm{\&Gamma;}}_D}{\&Integral;}{g}_D(v\&CenterDot;n)\mathrm{ds}$

$\underset{\mathrm{\&Omega;}}{\&Integral;}(\&dtri;\&CenterDot;u)\mathrm{qdx}-\underset{\mathrm{\&Omega;}}{\&Integral;}c\&CenterDot;\mathrm{pqdx}=-\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{fqdx}}^{\&prime;}$ 方程21

对于全部

在Γ_N上v·n＝0以及q∈L₂(Ω)。在下面特定描述中，考虑方程20中的公式，尽管全部结论可应用于方程21中的公式而没有通用性损失。

宏观复合混合公式

如果Ω被分割为m个不重叠多面体子域E_s，具有边界

和在边界

s，t＝1，...，m之间的界面，那么在示范实施例中，子域可被认为是粗糙计算单元E∈Ω_H或精细计算单元e∈Ω_h。可假设全部非零界面Γ_st是简单连接的分段平坦表面的段，s，t＝1，...，m。因此，全部非零界面Γ_st的并集可由Γ表示，即Γ＝∪Γ_st，并且Γ_N与E_s的交点可由Γ_Ns，s＝1，...，m表示。

令项V_s＝H_div(E_s)，Q_s＝L₂(E_s)，Λ_N，s＝L₂(Γ_N，s)以及Λ_st＝L₂(Γ_st)为在E_s中定义的向量函数u和函数p的空间，并分别令函数λ在Γ_N，s上定义和函数λ在Γ_st上定义，s，t＝1，...，m。新空间然后可被定义为：

V＝V₁×V₂×...×V_m

Q＝Q₁×Q₂×...×Q_m

Λ_N＝Λ_N，1×Λ_N，2×...×Λ_N，m。    方程22

${\mathrm{\&Lambda;}}_{\mathrm{\&Gamma;}}=\underset{1\&le;s\&le;t\&le;m}{\mathrm{\&Pi;}}{\mathrm{\&Lambda;}}_{\mathrm{st}}$

Λ＝Λ_Γ×Λ_N

因此，在方程18和19中示出的微分问题的宏观复合混合公式可以理解为：寻找(u，p，λ)∈V×Q×Λ，以使对于任何(v，q，μ)∈V×Q×Λ，满足E_s中的方程：

$\underset{E_s}{\&Integral;}\left(K^{-1}{u}_s\right)\&CenterDot;v_s\mathrm{dx}-\underset{E_s}{\&Integral;}{p}_s(\&dtri;\&CenterDot;v_s)\mathrm{dx}+\underset{{\mathrm{\&Gamma;}}_s}{\&Integral;}\mathrm{\&lambda;}(v_s\&CenterDot;n_s)\mathrm{ds}=-\underset{{\mathrm{\&Gamma;}}_{D,s}}{\&Integral;}{g}_D(v_s\&CenterDot;n_s)\mathrm{ds}$

$-\underset{E_s}{\&Integral;}(\&dtri;\&CenterDot;u_s)q_s\mathrm{dx}-\underset{E_s}{\&Integral;}c\&CenterDot;p_s{q}_s\mathrm{dx}=-\underset{E_s}{\&Integral;}{\mathrm{fq}}_s\mathrm{dx}$ 方程23

其中s＝1，...，m，其中在Γ_st上正常通量的连续度的变分方程：

$\underset{{\mathrm{\&Gamma;}}_{\mathrm{st}}}{\&Integral;}[u_s\&CenterDot;n_s+u_t\&CenterDot;n_t]{\mathrm{\&mu;}}_{\mathrm{st}}\mathrm{ds}=0,s,t=1,...,m,$ 方程24

以及其中诺依曼(Neumann)边界条件的变分方程：

$\underset{{\mathrm{\&Gamma;}}_{N,s}}{\&Integral;}(u_s\&CenterDot;n_s){\mathrm{\&mu;}}_{N,s}\mathrm{ds}=\underset{{\mathrm{\&Gamma;}}_{N,s}}{\&Integral;}{g}_N{\mathrm{\&mu;}}_{N,s}\mathrm{ds},$ 方程25

s＝1，...，m。这里，n_s是到

的单位向外法线，并且

分别是边界

的非狄利克雷和狄利克雷部分，s＝1，...，m。因此，u_s∈V_s和p_s∈Q_s分别是E_s中u∈V和p∈Q的函数分量，s＝1，...，m。

在棱柱网格上的混合有限元法和在微单元上“对流扩散”有限元空间 的定义

如前面讨论，计算网格Ω_h由是垂直棱柱体或棱锥体或四面体的元素{e_k}构成。为将针对在方程23-25中示出的问题的混合有限元(MFE)法公式化，必须定义空间V_h、Q_h和Λ_h的有限元子空间。

为定义通量向量函数的有限元空间，可假设每个棱柱计算网格单元e∈Ω_h被分割为三个四面体Δ₁、Δ₂和Δ₃，并且每个棱锥计算网格单元e∈Ω_h被分割为两个四面体Δ₁和Δ₂。因此，RT₀(e)可基于上面的e分割为四边形的分割表示向量函数的经典最低阶Raviart-Thomas有限元空间。

如果e是具有s个平坦面f_i，i＝1，...，s的Ω_h中的计算网格单元，那么对于“垂直”棱柱体和棱锥体s＝5并且对于四面体s＝4。通量矢量函数的在e上的有限元空间V_h(e)可然后定义为：

V_h(e)＝{v_h：v_h∈RT₀(e)，v_h·n_e＝const_i在f_i上，

这里，n_e是到e的边界

的向外单位法线。

有限元空间Q_h(e)可由下面的表达式针对解函数p定义：

Q_h(e)＝{q_h：q_h＝常数在e中}

如果E是Q_H中的宏单元，那么E可假设为是微单元e∈Q_h的并集。有限元空间V_h(E)可被定义为向量函数v_h的集合，其满足两个条件。第一，对于任何

v_h|_e∈v_h(e)，并且第二，v_h的法向分量在任何两个邻近精细网格单元e，

之间的界面上连续。因此，解函数的有限元空间Q_h(E)可由下面定义：

Q_h(E)＝{q_h：q_h|_e∈Q_h(e)对于全部}。

针对在分割为宏单元E_s，s＝1，...，m的Ω_h上的通量向量函数和解函数的全部有限元空间可用相似于方程22的方式分别定义为V_h＝V_h，1×V_h，2×...×V_h，m和Q_h＝Q_h，1×Q_h，2×...×Q_h，m。这里V_h，s＝V_h(E_s)并且Q_h，s＝Q_h(E_s)，s＝1，...，m。进一步地，拉格朗日乘子的有限元空间Λ_h≡Λ_h(Γ∪Γ_N)由下面定义：

Λ_h＝{λ_h：λ_h|_f≡const_f在Q_h中的任何面f上，以使

}。

在Ω _h 上的宏观复合混合有限元法

在方程23-25中说明的问题的宏观复合混合有限元离散化可被描述为：寻找(u_h，p_h，λ_h)∈V_h×Q_h×Λ_h，以使对于任何(v，q，μ)∈V_h×Q_h×Λ_h，满足E_s中的方程：

$\underset{E_s}{\&Integral;}\left(K^{-1}{u}_{h,s}\right)\&CenterDot;v_s\mathrm{dx}-\underset{E_s}{\&Integral;}{p}_{h,s}(\&dtri;\&CenterDot;v_s)\mathrm{dx}+\underset{{\mathrm{\&Gamma;}}_s}{\&Integral;}{\mathrm{\&lambda;}}_h(v_s\&CenterDot;n_s)\mathrm{ds}=-\underset{{\mathrm{\&Gamma;}}_{D,s}}{\&Integral;}{g}_D(v_s\&CenterDot;n_s)\mathrm{ds}$

$-\underset{E_s}{\&Integral;}(\&dtri;\&CenterDot;u_{h,s})q_s\mathrm{dx}-\underset{E_s}{\&Integral;}c\&CenterDot;p_{h,s}{q}_s\mathrm{dx}=-\underset{E_s}{\&Integral;}{\mathrm{fq}}_s\mathrm{dx}$ 方程26

s＝1，...，m，其中在Γ_st上正常通量的连续度的变分方程：

$\underset{{\mathrm{\&Gamma;}}_{\mathrm{st}}}{\&Integral;}[u_{h,s}\&CenterDot;n_s+u_{h,t}\&CenterDot;n_t]{\mathrm{\&mu;}}_{\mathrm{st}}\mathrm{ds}=0,s,t=1,...,m,$ 方程27

以及其中诺依曼边界条件的变分方程：

$\underset{{\mathrm{\&Gamma;}}_{N,s}}{\&Integral;}(u_{h,s}\&CenterDot;n_s){\mathrm{\&mu;}}_{N,s}\mathrm{ds}=\underset{{\mathrm{\&Gamma;}}_{N,s}}{\&Integral;}{g}_N{\mathrm{\&mu;}}_{N,s}\mathrm{ds},$ 方程28

s＝1，...，m。

由方程26-28说明的有限元问题导致代数方程：

$M_s{\stackrel{\&OverBar;}{u}}_s+B_s^T{\stackrel{\&OverBar;}{p}}_s+C_s^T\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{g}}_{D,s}$

$B_s{\stackrel{\&OverBar;}{u}}_s-{\mathrm{\&Sigma;}}_s{\stackrel{\&OverBar;}{p}}_s={\stackrel{\&OverBar;}{f}}_s,$ 方程29

s＝1，...，m，由以下代数方程补充

$C\&CenterDot;\left[\begin{array}{c}{\stackrel{\&OverBar;}{u}}_1\\ .\\ .\\ .\\ \stackrel{\&OverBar;}{u_m}\end{array}\right]={\stackrel{\&OverBar;}{g}}_N$ 方程30

方程29和30表示在Ω_H中邻近的宏单元之间界面上的法向通量的连续度条件和在Γ_N上的诺依曼边界条件。这里，M_s是正方形n_u，s×n_u，s对称正定矩阵(在通量空间中的质量矩阵)，B_s是矩形n_p，s×n_p，s矩阵，C_s ^T是矩形n_u，s×n_λ矩阵，∑_s是对角n_p，s×n_p，s矩阵，其中n_u，s＝dimV_h，s并且n_p，s＝dim Q_h，s，s＝1，...，m，以及n_λ＝dimΛ_h。

在方程29和30中说明的系统可用更紧凑形式呈现为：

$\left(\begin{array}{ccc}M& B^T& C^T\\ B& -\mathrm{\&Sigma;}& 0\\ C& 0& 0\end{array}\right)\left(\begin{array}{c}\stackrel{\&OverBar;}{u}\\ \stackrel{\&OverBar;}{p}\\ \stackrel{\&OverBar;}{\mathrm{\&lambda;}}\end{array}\right)=\left(\begin{array}{c}\stackrel{\&OverBar;}{g_D}\\ \stackrel{\&OverBar;}{f}\\ \stackrel{\&OverBar;}{g_N}\end{array}\right),$ 方程31

其中 $M=M_1\&CirclePlus;M_2\&CirclePlus;...\&CirclePlus;M_m$ 和 $B=B_1\&CirclePlus;B_2\&CirclePlus;...\&CirclePlus;B_m$ 是m×m分块对角矩阵，

C＝(C₁...C_m)， $\stackrel{\&OverBar;}{u}=\left(\begin{array}{c}{\stackrel{\&OverBar;}{u}}_1\\ .\\ .\\ .\\ {\stackrel{\&OverBar;}{u}}_m\end{array}\right),$ $\stackrel{\&OverBar;}{p}=\left(\begin{array}{c}{\stackrel{\&OverBar;}{p}}_1\\ .\\ .\\ .\\ {\stackrel{\&OverBar;}{p}}_m\end{array}\right),$ 以及 $\stackrel{\&OverBar;}{\mathrm{\&lambda;}}=R^{n_{\mathrm{\&lambda;}}}.$

在具有粗糙有限元空间的Ω _H 上的离散化

令E_s，s＝1，...，m，为具有面F_s，j，j＝1，...，r_s的Ω_H中的宏单元，其中r_s等于5或4(例如，表示棱柱或棱锥单元)。先前定义的有限元空间V_h(E_s)可然后呈现为(r_s+1)个子空间的直和：

$V_h\left(E_s\right)=W_{h,s,1}\&CirclePlus;W_{h,s,2}\&CirclePlus;...\&CirclePlus;W_{h,s,r_s}\&CirclePlus;W_{h,s,\mathrm{int}},$ 方程32

其中空间W_h，s，j关联面F_s，j，j＝1，...，r_s上的法向通量的自由度(DOF)，并且空间W_h，s，int关联E_s内部中法向通量的内部自由度。

如果{W_j，i}是在W_h，s，j，j＝1，...，r_s中的基底，并且{w_int，i}是W_h，s，int中的基底，那么对于向量函数v∈V_h(E_s)，关于这些基底的自由度的向量v可由下面的方程呈现：

$v^T=(v_1^T,v_2^T,...,v_{r_s}^T,v_{\mathrm{int}}^T).$ 方程33

如果

是W_h，s，j，j＝1，...，r_s的子空间，并且

是W_h，s，int的子空间，那么

j＝1，...，r_s中的基底可表示为

并且

中的基底可表示为

因此，对于属于该空间的向量函数

${\hat{V}}_{h,s}\&equiv;{\hat{V}}_h\left(E_s\right)={\hat{W}}_{h,s,1}\&CirclePlus;{\hat{W}}_{h,s,2}\&CirclePlus;...\&CirclePlus;{\hat{W}}_{h,s,r_s}\&CirclePlus;{\hat{W}}_{h,s,\mathrm{int}},$ 方程34

关于这些基底的自由度的向量

可由下面的方程呈现：

$v_c^T=(v_{c,1}^T,v_{c,2}^T,...v_{c,r_s}^T,v_{c,\mathrm{int}}^T).$ 方程35

V_h(E_s)可被称为精细有限元空间，并且

可称为粗糙有限元空间。对于该表示，底指数c用于自由度的粗糙空间。

上面基底的选择唯一定义(r_s+1)变换矩阵R_s，j，j＝1，...，r_s和R_s，int，以使：

v_j＝R_s，j·v_c，j，j＝1，...，r_s和v_int＝R_s，int·v_c，int，方程36

空间V_h(E_s)和

的变换矩阵可由下面的方程定义：

$R_s=R_{s,1}\&CirclePlus;R_{s,2}\&CirclePlus;...\&CirclePlus;R_{s,r_s}\&CirclePlus;R_{s,\mathrm{int}},$ 方程37

进一步地，向量函数的全局粗糙空间可由下面的方程引入：

V_H＝V_H，1×V_H，2×...×V_H，m，             方程38

其中s＝1，...，m。因此，有限元向量函数v∈V_H的自由度向量v和v_c分别关联V_h和V_H中的基底，满足变换：

v＝R·v_c，                                 方程39

其中R是m×m分块对角矩阵R＝R₁×R₂×...×R_m。

为定义方程23-35中拉格朗日乘子的粗糙有限元空间，可假设Ω_H中两个邻近宏单元E和E’之间界面

上的有限元空间和

的正轨迹重合。更特定地，可假设

和

中选择的基底向量函数的F上正轨迹也重合。

如果F表示在Ω_H中两个邻近宏单元E和E’之间的界面，那么不损失通用性，该界面可关联E的面F_E，1。

的粗糙有限元子空间同样关联F以及基底向量函数

i＝1，...，n_F，其中n_F＝dim一组函数：

i＝1，...，n_F，可然后在F上定义。通过构建，这些函数是线性独立的。然后，空间被定义为

即集合

是Λ_H(F)中的基底。由于假设源自

和源自

的向量函数的F上的正轨迹组重合，因此在E和E’中的F上定义的空间Λ_H(F)也重合。

拉格朗日乘子的有限元空间Λ_H可被定义为在Γ∪Γ_N上定义的分段常数函数λ_H的集合，以使在Ω_H中邻近宏单元E和E’之间全部界面F上，以及在属于Γ_N的宏面F上λ_H|_F∈Λ_H(F)。如果F＝F_s，1是宏单元E_s∈Ω_H，s＝1，...，m的面，那么其可示出对于空间Λ_H(F)和Λ_H的后面定义，在空间Λ_h(F)和Λ_H(F)之间的变换矩阵R_λ，F与方程37中的变换矩阵R_s，1一致，并且在全局空间Λ_h和Λ_H之间的变换矩阵可由

定义，其中在Γ∪Γ_N上的全部不同宏面上采取直接求和

因此，空间Λ_h和Λ_H的自由度的向量λ和λ_c满足变换：

λ＝R_λ·λ_c。                        方程40

最终，解函数的粗糙空间可被简单定义为：

Q_H＝Q_h。                              方程41

在方程23-25中表示的宏观复合混合有限元离散化可然后理解为：

寻找(u_H，p_H，λ_H)∈V_H×Q_H×Λ_H，以使对于任何(v，q，μ)∈V_H×Q_H×Λ_H满足E_s中的方程：

$\underset{E_s}{\&Integral;}\left(K^{-1}{u}_{H,s}\right)\&CenterDot;v_s\mathrm{dx}-\underset{E_s}{\&Integral;}{p}_{H,s}(\&dtri;\&CenterDot;v_s)\mathrm{dx}+\underset{{\mathrm{\&Gamma;}}_s}{\&Integral;}{\mathrm{\&lambda;}}_H(v_s\&CenterDot;n_s)\mathrm{ds}=-\underset{{\mathrm{\&Gamma;}}_{D,s}}{\&Integral;}{g}_D(v_s\&CenterDot;n_s)\mathrm{ds}$

$-\underset{E_s}{\&Integral;}(\&dtri;\&CenterDot;u_{H,s})q_s\mathrm{dx}-\underset{E_s}{\&Integral;}c\&CenterDot;p_{H,s}{q}_s\mathrm{dx}=-\underset{E_s}{\&Integral;}{\mathrm{fq}}_s\mathrm{dx}$ 方程42

s＝1，...，m，其中在Γ_st上正常通量的连续度的变分方程：

$\underset{{\mathrm{\&Gamma;}}_{\mathrm{st}}}{\&Integral;}[u_{H,s}\&CenterDot;n_s+u_{H,t}\&CenterDot;n_t]{\mathrm{\&mu;}}_{\mathrm{st}}\mathrm{ds}=0,s,t=1,...,m,$ 方程43

以及其中诺依曼边界条件的变化方程：

$\underset{{\mathrm{\&Gamma;}}_{N,s}}{\&Integral;}(u_{H,s}\&CenterDot;n_s){\mathrm{\&mu;}}_{N,s}\mathrm{ds}=\underset{{\mathrm{\&Gamma;}}_{N,s}}{\&Integral;}{g}_N{\mathrm{\&mu;}}_{N,s}\mathrm{ds},$ 方程44

s＝1，...，m。

在方程42-44中的有限元问题导致代数方程：

${\hat{M}}_s{\stackrel{\&OverBar;}{u}}_s+{\hat{B}}_s^T{\stackrel{\&OverBar;}{p}}_s+{\hat{C}}_s^T\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{g}}_{D,s}$

${\hat{B}}_s{\stackrel{\&OverBar;}{u}}_s-{\mathrm{\&Sigma;}}_s{\stackrel{\&OverBar;}{p}}_s={\stackrel{\&OverBar;}{f}}_s,$ 方程45

s＝1，...，m，其由以下代数方程补充：

$\hat{C}\&CenterDot;\left[\begin{array}{c}{\stackrel{\&OverBar;}{u}}_1\\ .\\ .\\ .\\ \stackrel{\&OverBar;}{u_m}\end{array}\right]={\stackrel{\&OverBar;}{g}}_N.$ 方程46

方程46表示在Ω_H中邻近的宏单元之间界面上的法向通量的连续度条件以及在Γ_N上的诺依曼边界条件。在这些方程中，是正方形

对称正定矩阵(在通量空间中的质量矩阵)，

是矩形

矩阵，

是矩形

矩阵，∑_s是对角

矩阵，其中

s＝1，...，m，以及

由于通过定义Q_H＝Q_h并特别地Q_H，s＝Q_h，s，s＝1，...，m，

那么在方程45中的矩阵∑_s与在方程29中的矩阵∑_s一致，s＝1，...，m。

在方程45和46中的系统可用更紧凑形式呈现为：

$\left(\begin{array}{ccc}\hat{M}& {\hat{B}}^T& {\hat{C}}^T\\ \hat{B}& -\mathrm{\&Sigma;}& 0\\ \hat{C}& 0& 0\end{array}\right)\left(\begin{array}{c}\stackrel{\&OverBar;}{u}\\ \stackrel{\&OverBar;}{p}\\ \stackrel{\&OverBar;}{\mathrm{\&lambda;}}\end{array}\right)=\left(\begin{array}{c}\stackrel{\&OverBar;}{g_D}\\ \stackrel{\&OverBar;}{f}\\ \stackrel{\&OverBar;}{g_N}\end{array}\right),$ 方程47

其中 $\hat{M}={\hat{M}}_1\&CirclePlus;{\hat{M}}_2\&CirclePlus;...\&CirclePlus;{\hat{M}}_m$ 和 $\hat{B}={\hat{B}}_1\&CirclePlus;{\hat{B}}_2\&CirclePlus;...\&CirclePlus;{\hat{B}}_m$ 是m×m分块对角矩阵，

$\hat{C}=({\hat{C}}_1...{\hat{C}}_m),$ $\stackrel{\&OverBar;}{u}=\left(\begin{array}{c}{\stackrel{\&OverBar;}{u}}_1\\ .\\ .\\ .\\ {\stackrel{\&OverBar;}{u}}_m\end{array}\right),$ $\stackrel{\&OverBar;}{p}=\left(\begin{array}{c}{\stackrel{\&OverBar;}{p}}_1\\ .\\ .\\ .\\ {\stackrel{\&OverBar;}{p}}_m\end{array}\right),$ 以及 $\stackrel{\&OverBar;}{\mathrm{\&lambda;}}=R^{{\hat{n}}_{\mathrm{\&lambda;}}}.$

在这里给出的定义下，可示出：

${\hat{M}}_s=R_s^T{M}_s{R}_s$

${\hat{B}}_s=B_s{R}_s$ 方程48

${\hat{C}}_s^T=R_s^T{C}_s^T{R}_{\mathrm{\&lambda;}}$

s＝1，...，m，因此产生的全局矩阵可由下面的方程表示：

$\hat{M}=R^T\mathrm{MR}$

$\hat{B}=\mathrm{BR}$ 方程49

$\hat{C}=R_{\mathrm{\&lambda;}}^T\mathrm{CR}$

在Ω _H 上的均匀化的离散化

可用先前定义的粗糙有限元空间为方程42-44的离散化推导等效代数系。使用该定义，可能重解释和拉格朗日乘子关联的自由度的定义。这可通过考虑Ω_H中的宏单元E_s，s＝1，...，m并选择

中的基底向量函数

执行。然后对应该基底函数的在方程26中示出的有限元方程可由下面方程给出：

$\underset{E_s}{\&Integral;}\left(K^{-1}{\hat{u}}_h\right)\&CenterDot;{\hat{w}}_h\mathrm{dx}-\underset{E_s}{\&Integral;}{\hat{p}}_h(\&dtri;\&CenterDot;{\hat{w}}_h)\mathrm{dx}+\underset{{\mathrm{\&Gamma;}}_s}{\&Integral;}{\mathrm{\&lambda;}}_h({\hat{w}}_h\&CenterDot;n_s)\mathrm{ds}=-\underset{{\mathrm{\&Gamma;}}_{D,s}}{\&Integral;}{g}_D({\hat{w}}_h\&CenterDot;n_s)\mathrm{ds}.$ 方程50

如果基底函数

属于关联法向通量的内部自由度的空间，或属于关联其中施加狄利克雷边界条件的外部边界的空间，那么方程50的第三项放弃。

可选择属于Γ_s的宏单元E_s的面F≡F_s，j，j＝1，...，r_s，并且

可被选为

的基底向量函数的一个，其中假设在F上的其面平均值非负，即

因此，可对应该基底函数的在方程50中提供的有限元方程可用等效形式写为：

$\underset{E_s}{\&Integral;}\left(K^{-1}{\hat{u}}_h\right)\&CenterDot;{\hat{w}}_h\mathrm{dx}-\underset{E_s}{\&Integral;}{\hat{p}}_h(\&dtri;\&CenterDot;{\hat{w}}_h)\mathrm{dx}+d_w\&CenterDot;{\mathrm{\&lambda;}}_w=0,$ 方程51

其中：

${\mathrm{\&lambda;}}_w=\frac{1}{d_w}\underset{F}{\&Integral;}{\mathrm{\&lambda;}}_h({\hat{w}}_h\&CenterDot;n_s)\mathrm{ds}.$

在该形式中，λ_w可被解释为关联特定选择的基底向量函数

的拉格朗日乘子λ_h的新自由度。如果假设在方程34中的基底向量函数 $\left\{{\hat{w}}_{j,i}\right\}\&Element;{\hat{W}}_{h,s,j}$ 满足条件：

$d_{j,i}=\underset{{\&PartialD;E}_s}{\&Integral;}{\hat{w}}_{j,i}{n}_s\mathrm{ds}>0,j=1,...,r_s,$

那么，对于给定的宏单元E_s，方程45中的第一方程组可写作

${\hat{M}}_s{\stackrel{\&OverBar;}{u}}_s+{\hat{B}}_S^T{\stackrel{\&OverBar;}{p}}_s+{\hat{C}}_S^T\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{g}}_{D,s}$ 方程53

由于拉格朗日乘子的自由度的不同解释，其具有相同向量

但具有不同矩阵

和不同向量

通量向量的分量可分割为属于E_s的边界

的其分量(由指数“F”表示)，以及与关于E_s的内部自由度关联的其分量(由指数“I”表示)。然后方程53可以以2×2块形式呈现：

$\left(\begin{array}{cc}{\hat{M}}_{s,F}& {\hat{M}}_{s,\mathrm{FI}}\\ {\hat{M}}_{s,\mathrm{IF}}& {\hat{M}}_{s,I}\end{array}\right)\left(\begin{array}{c}{\stackrel{\&OverBar;}{u}}_{s,F}\\ {\stackrel{\&OverBar;}{u}}_{s,I}\end{array}\right)+\left(\begin{array}{c}{\hat{B}}_{s,F}^T\\ {\hat{B}}_{s,I}^T\end{array}\right){\stackrel{\&OverBar;}{p}}_s+\left(\begin{array}{c}{\hat{C}}_{s,F}^T\\ 0\end{array}\right)\mathrm{\&lambda;}=\left(\begin{array}{c}{\stackrel{\&OverBar;}{g}}_{D,{\mathrm{\&Gamma;}}_s}\\ 0\end{array}\right),$ 方程54

其中

是的子向量，对应Γ_s上的面。

可考虑描述方程54中矩阵

的简单例子，其中域Ω由两个宏单元E₁和E₂构成，具有界面F和边界

在此情况下，

和

的非零块等于具有对角元素的相同对角t×t矩阵D_F：

$d_{F,i}=\underset{F}{\&Integral;}{\hat{w}}_{1,i}{n}_1\mathrm{ds},$

其中

是在关联F的E₁上的基底向量函数，i＝1，...，t，n₁是关于E₁向外的到F的单位法线，并且t是E₁上基底向量函数的总数。因此，在该例子中，在方程54中的矩阵

可表示为：

${\left({\hat{C}}_{1,F}^T\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)}_F={\left({\hat{C}}_{2,F}^F\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)}_F=D_F\stackrel{\&OverBar;}{\mathrm{\&lambda;}}.$

可引入矩阵Q_F，其中Q_F具有项：

$q_{F,\mathrm{ij}}=\underset{F}{\&Integral;}({\hat{w}}_{1,i},n_1)({\hat{w}}_{1,j}\&CenterDot;n_1)\mathrm{ds},i,j=0,...,t.$

使用针对拉格朗日乘子的在方程52中示出的自由度的解释，可能通过以下变换连接方程54中的向量

与方程47中的向量

${\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_{\mathrm{new}}=D_F^{-1}{Q}_F{\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_{\mathrm{old}}.$ 方程55

在方程55中示出的变换可用对角矩阵D_λ和分块对角矩阵Q_λ延伸到网格Ω_H，其中Ω_H中邻近宏单元之间的每个界面或属于Γ_N的Ω_H中宏单元的每面一个对角块。可注意在粗糙计算网格Ω_H上，有限元子空间满足与最精细计算网格Ω_H的约束相似的约束。特别地，元向量函数在两个邻近宏单元之间的界面上，以及宏单元边界与边界的诺依曼部分(宏诺依曼面)的交点上具有恒定法向分量。由于每个宏面通过更精细网格上计算面的组合形成，因此等于界面和诺依曼面总数的有限元子空间的大小随着网格在分层结构中前进(从更精细到更粗糙网格)减小。

使用在方程52中呈现的拉格朗日乘子λ的自由度的定义，以及在方程54中引入的通量函数的新分割，Ω_H中每个宏单元E_s的方程45和46的系统可用下面的代数形式呈现：

${\hat{M}}_{s,F}{\stackrel{\&OverBar;}{u}}_{s,F}+{\hat{M}}_{s,\mathrm{FI}}{\stackrel{\&OverBar;}{u}}_{s,I}+{\hat{B}}_{s,F}^T{\stackrel{\&OverBar;}{p}}_s+{\hat{C}}_{s,F}^T\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{g}}_{D,{\mathrm{\&Gamma;}}_s}$

${\hat{M}}_{s,\mathrm{IF}}{\stackrel{\&OverBar;}{u}}_{s,F}+{\hat{M}}_{s,I}{\stackrel{\&OverBar;}{u}}_{s,I}+{\hat{B}}_{s,I}^T{\stackrel{\&OverBar;}{p}}_s=0$ 方程56

${\hat{B}}_{s,F}{\stackrel{\&OverBar;}{u}}_{s,F}+{\stackrel{\&OverBar;}{B}}_{s,I}{\stackrel{\&OverBar;}{u}}_{s,1}-{\mathrm{\&Sigma;}}_s{\stackrel{\&OverBar;}{p}}_s={\stackrel{\&OverBar;}{f}}_s$

s＝1，...，m，由下面的方程：

$\hat{C}\stackrel{\&OverBar;}{u}={\stackrel{\&OverBar;}{g}}_N$ 方程57

针对在Ω_H上的宏单元之间界面上法向通量的连续度和在Γ_N上的诺依曼边界条件补充。应注意内部和边界面之间的分离是下面呈现的方法的显著方面。

用于本技术的示范实施例中的均匀化算法由两个主要步骤构成。在第一步骤，假设矩阵 $\left(\begin{array}{cc}{\hat{M}}_{s,I}& {\hat{B}}_{s,I}^T\\ {\hat{B}}_{s,I}& -{\mathrm{\&Sigma;}}_s\end{array}\right)$ 非奇异，那么子向量

和可在由方程56描述的系统中消除。例如，如果方程18中的系数c在E_s，s＝1，...，m中不等于0，那么这些矩阵是非奇异的。

在消除步骤后下面的代数系统导致：

${\tilde{M}}_{s,F}{\stackrel{\&OverBar;}{u}}_{s,F}+{\hat{C}}_{s,F}^T\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{D,s},s=1,...,m,$ 方程58

其由方程57补充。

在第二步骤，可引入新自由度。该自由度是限于宏单元E_s的主变量p_H，s的值。因此，在方程57和58中表示的系统可由下面的系统替换：

$M_{H,s}{\stackrel{\&OverBar;}{u}}_{H,s}+B_{H,s}^T{\stackrel{\&OverBar;}{p}}_{H,s}+{\hat{C}}_{H,s}^T{\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_H={\stackrel{\&OverBar;}{g}}_{D,H,s}$

${\hat{B}}_{H,s}{\stackrel{\&OverBar;}{u}}_{H,s}-{\mathrm{\&Sigma;}}_{H,s}{\stackrel{\&OverBar;}{p}}_{H,s}={\stackrel{\&OverBar;}{f}}_{H,s}$ 方程59

其中s＝1，...，m。这由下面的方程补充：

$C_H{\stackrel{\&OverBar;}{u}}_H={\stackrel{\&OverBar;}{g}}_{N,H}$ 方程60

其中

并且C_H＝(C_H，1...C_H，m)。

方程59中的矩阵B_H，s和∑_H，s以及值f_H，s可通过在Ω_H上的宏单元E_s，s＝1，...，m上积分方程18中守恒定律方程

来推导。令Γ_s为属于Ω∪Γ_N的E_s的边界的部分，并且令{w_s，1，...，w_s，qs}为关联s＝1，...，m的V_H，s中全部基底向量函数的集合。然后E_s上的有限元守恒定律以下面代数方程的形式获得：

$\underset{i=1}{\overset{q_s}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{s,i}{u}_{H,s,i}+{\mathrm{\&Sigma;}}_{H,s}{p}_{H,s}=-f_{H,s},$

其中：

${\mathrm{\&gamma;}}_{s,i}=\underset{{\&PartialD;E}_s}{\&Integral;}{w}_{s,i}\&CenterDot;n_s\mathrm{ds},i=1,...,q_s$

${\mathrm{\&Sigma;}}_{H,s}=\underset{E_s}{\&Integral;}\mathrm{cdx}=\left|E_s\right|\&CenterDot;\stackrel{\&OverBar;}{c}$

$u_{H,s,i}={\mathrm{\&gamma;}}_{s,i}^{-1}\underset{\&PartialD;E_s}{\&Integral;}u\&CenterDot;n_s\mathrm{ds},i=1,...,q_s,$

$p_{H,s}={\mathrm{\&Sigma;}}_{H,s}^{-1}\underset{E_s}{\&Integral;}c\&CenterDot;\mathrm{pdx}$

$f_{H,s}=-\underset{E_s}{\&Integral;}\mathrm{fdx}=-\left|E_s\right|\&CenterDot;\stackrel{\&OverBar;}{f}$

以及s＝1，...，m。在这些方程中，γ_s，i是第i个边界面的面积，|E_s|是宏单元的体积，并且在上面的杆表示单元E_s的体积平均值。进一步地，u_H，s，i，i＝1，...，q_s和p_H，s表示通量矢量函数的新自由度和E_s中解函数的均匀化自由度。因此，方程59中的矩阵B_H，s、M_H，s和向量可由下面的方程分别定义：

$B_{H,s}=-({\mathrm{\&gamma;}}_{s,1}...{\mathrm{\&gamma;}}_{s,q_s})\&Element;R^{1\&times;q_s},$ 方程61

$M_{H,s}={\tilde{M}}_{s,F}-\frac{1}{{\mathrm{\&Sigma;}}_{H,s}}{B}_{H,s}^T{B}_{H,s},$ 以及                    方程62

${\stackrel{\&OverBar;}{g}}_{D,H,s}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{D,s}-\frac{1}{{\mathrm{\&Sigma;}}_{H,s}}{B}_{H,s}^T{f}_{H,s},$ 方程63

在方程61-63中呈现的条件下，方程58中的方程等效于在方程59中示出的方程，并且可通过从方程59消除p_H，s，s＝1，...，m来获得。通过方程59和60表示的系统可被称为在具有有限元空间V_H、Q_H和Λ_H的粗糙网格Ω_H上方程18和19的均匀化的离散化。

通量矢量函数的正交分量的粗化算法

如果E是Ω_H中的宏单元，那么不损失通用性，可假设E具有与子域(或如关于图8和9讨论的地质层)Ω₁，Ω₂，...，Ω_t的非零交点，其中t是正整数，1≤t≤N_z。为描述粗化算法，可施加三个主要假设。第一，尖灭的边界属于Ω_H中宏单元E的侧向边缘的联合。第二，与Ω_l，l＝1，...，t的边界的交点是平坦的。第三，扩散张量K在每个Ω_l中恒定，即在E∩Ω_l，l＝1，...，t中K≡K_l。

根据这些假设，E的侧面是三角形并属于Ω₁和Ω_t内部，或属于这些子域的边界。为此，可假设V_H(E)中有限元向量函数的正交分量在E的顶部和底部侧面上恒定。这可以是粗化算法的第一步骤。

如果F是E的垂直面，并且F_l表示F与子域902Ω_l，l＝1，...，t的交点，那么面F是四边形和三角形(若有的话)的并集，如在图11和12中示出的。

图11A和11B是图解，其根据本技术的示范实施例示出垂直四边形面分割为子面。垂直四边形面一般被称为F 1100。在图11A中，F 1100被分割为子面F_l 1102，l＝1，...，4。在图11B中，F 1100被分割为子面F_l 1104，l＝1，...，5。图12是图解，其根据本技术的示范实施例示出垂直三角形面F分割为子面F_l，l＝1，...，4。E与Ω_l，l＝1，...，t的交点可由F_l表示。进一步地，在E_l-1和E_l，l＝2，...，t之间的边界/界面可由Γ_l-1，l表示。

在粗化算法的第二步骤中，有限元通量向量函数的正交分量恒定的条件可施加在每个子面F_l，l＝1，...，t上。对于该步骤矩阵R_j≡R_F，例如在方程37中示出，是t×t分块对角矩阵，其中对角块是全部分量等于1的列向量。导致的有限元通量向量函数空间在每个子面F_l，l＝1，...，t上具有恒定正交分量。这些分量可由ξ_l，l＝1，...，t表示。

在第三粗化步骤中，分段恒定向量场v_l＝(v_l，1，v_l，2，v_l，3)^T可在每个E_l，l＝1，...，t中引入，服从下面的条件：

2.在F_l，l＝1，...，t上    v_l·n_E＝ξ_l，    方程64

3.在Γ_l-1，l，l＝2，...，t上    (v_l-1-v_l)·n_l-1，l＝0，

其中n_l-1，l是从E_l-1引导到E_l，l＝2，...，t的Γ_l-1，l上的单位法线。另一条件可施加在上面的分段恒定向量函数v上。特定地，可假设分段平滑连续函数ψ＝ψ(x)存在，以使：

在E中 $v=-K\&CenterDot;\&dtri;\mathrm{\&psi;}.$

上面的假设暗示针对函数v的以下约束组：

在Γ_l-1，l，l＝2，...，t上[K^-1(v_l-1-v_l)]×n_l-1，l＝0，方程65即，假设向量函数K^-1v的切向分量在Γ_l-1，l，l＝2，...，t上连续。

条件可在代数系统中组合：

$N\stackrel{\&OverBar;}{v}=\stackrel{\&OverBar;}{\mathrm{\&xi;}},$ 方程66

其中N是n_ξ×n_v矩阵，向量具有大小n_v＝3t，具有分量v_l，1，v_l，2，v_l，3，l＝1，...，t，以及向量

具有大小n_ξ＝3(t-1)+t＝4t-3，具有等于ξ_l，l＝1，...，t的t个分量并且没有其它3(t-1)个分量。

可注意到如果向量v_l，l＝1，...，t中的一个已知或指定(例如向量v₁)，那么全部其它向量(即向量v₂，v₃，...，v_t)可使用在方程64中示出的条件2和在方程65中呈现的条件唯一定义。可遵循矩阵N的秩不可大于3。

因此，针对方程66中系统恒定的条件是ξ_l，l＝1，...，t的任何t-3个值应呈现为具有独立于值ξ_l，l＝1，...，t的系数的三个其它值的线性组合。对于t＞0，集合ξ_l，l＝1，...，t的粗化的代数公式可假设方程66中矩阵N的秩等于3来明确推导。

特定情况的粗化算法

在下面讨论两个特定情况的粗化算法。不损失通用性，可假设面F正交于坐标轴x₁。

对于第一情况，可假设界面边界Γ_l-1，l，l＝2，...，t平行于坐标平面(x₁，x₂)，并且扩散矩阵由下面的公式定义

$K_i=\left(\begin{array}{ccc}{K}_{l,1}& 0& 0\\ 0& K_{l,22}& K_{l,23}\\ 0& K_{l,32}& K_{l,33}\end{array}\right),$ l＝1，...，t。

代数分析可用来示出在此情况下，向量函数v的F上的法向分量ξ_l通过下面的关系连接

k_l-1，1ξ_l-1＝k_l，1ξ_l，l＝2，...，t。

如果挑选ξ₁为独立变量，那么其它t-1个分量由递归公式定义

$\frac{{\mathrm{\&xi;}}_l}{{\mathrm{\&xi;}}_{l-1}}=\frac{k_{l-\mathrm{1,1}}}{k_{l,1}},l=2,...,t.$

并且是R^t中的向量列的对应变换矩阵R_F等于 $R_F={(1,k_{\mathrm{1,1}}{k}_{\mathrm{2,1}}^{-1}{k}_{\mathrm{2,1}}{k}_{\mathrm{3,1}}^{-1},...,k_{t-\mathrm{1,1}}{k}_{t,1}^{-1})}^T.$

第二重要情况可基于界面边界Γ_l-1，l，l＝2，...，t互相平行并正交于F的垂直边缘，并且扩散矩阵具有关于界面Γ_l-1，l，l＝2，...，t的特殊结构的假设。例如，对于分段恒定标量扩散张量，方程66中矩阵N的秩等于2。

多级方法

在前面章节中仅针对两个网格：精细网格Ω_h和粗糙网格Ω_H的情况描述均匀化算法。该方法也可延伸到网格的分层序列

Ω_h，0，Ω_h，1，...，Ω_h，L，

其中L≥1是整数，Ω_h＝Ω_h，0是精细网格，Ω_H＝Ω_h，L是粗糙网格。网格Ω_h，l通过应用上面描述的粗化算法从网格Ω_h，l-1，l＝1，...，L获得。如描述的，粗化算法可以是任意的。例如，在一个特别情况下，每个“粗糙”棱柱体E∈Ω_h，l由八个“精细”棱柱体e∈Ω_h，l-1，l＝1，...，L构成，如在图10中所示的，而在粗化的另一例子中，每个“粗糙”棱柱体E∈Ω_h，l可由四个“精细”棱柱体e∈Ω_h，l-1，l＝1，...，L构成，如在图13中示出的。

图13是示意图，其根据本技术的实施例图解说明粗糙棱柱体1300划分为四个精细棱柱体1302。

多级方法可向在此呈现的技术提供一些进一步优点。特定地，如果在先前章节中描述的粗化算法被直接应用于网格Ω_h＝Ω_h，0和Ω_H＝Ω_h，L，那么将大尺寸矩阵求逆。例如，如果没有尖灭，那么每个网格单元E∈Ω_H由4^L个网格单元e∈Ω_h构成。为在E中进行代数粗化，尺寸等于E中网格单元e之间界面数目的矩阵必须被求逆。该数目具有O(4^L)的量级。例如，如果L＝3，那么大于约64(即4³)的矩阵被求逆。

为减少计算负荷，二级方法可用多级方法代替。在多级方法中，精细网格系统可在网格Ω_h＝Ω_h，0上构造，然后表示Ω_H＝Ω_h，1，并应用该算法从而在Ω_H上获得粗糙系统。因此，没有大于尺寸四的矩阵被求逆。粗化规程可然后在网格Ω_h，1上重复。在此情况下，Ω_h＝Ω_h，1可被定义为新精细网格，并且Ω_H＝Ω_h，2可被定义为新粗糙网格。通过应用相同的算法L次，最终系统转为新粗糙网格Ω_H＝Ω_h，L。粗化算法在下面进一步详细描述。

粗化算法的初始化可通过考虑Ω_h＝Ω_h，0为精细网格来执行。可然后应用复合混合公式(即，可在精细网格单元之间的全部界面上以及在边界Γ_N的诺依曼部分上引入拉格朗日乘子)。这导致下面形式的代数系统：

$\left(\begin{array}{ccc}{M}_0& B_0^T& C_0^T\\ B_0& -{\mathrm{\&Sigma;}}_0& 0\\ C_0& 0& 0\end{array}\right)\left(\begin{array}{c}{u}_0\\ p_0\\ {\mathrm{\&lambda;}}_0\end{array}\right)=\left(\begin{array}{c}{F}_{u,0}\\ F_{p.0}\\ F_{\mathrm{\&lambda;},0}\end{array}\right).$

这里，次标0表示全部矩阵在网格Ω_h，0上定义。然后，通过求逆对角矩阵∑₀执行p₀，从而按照(u₀，λ₀)获得系统：

$\left(\begin{array}{cc}{A}_0& C_0^T\\ C_0& 0\end{array}\right)\left(\begin{array}{c}{u}_0\\ {\mathrm{\&lambda;}}_0\end{array}\right)=\left(\begin{array}{c}{\hat{F}}_{u,0}\\ F_{\mathrm{\&lambda;},0}\end{array}\right),$

其中 $A_0=M_0+B_0^T{\mathrm{\&Sigma;}}^{-1}{B}_0,$ 并且 ${\hat{F}}_{u,0}=F_{u,0}+B_0^T{\mathrm{\&Sigma;}}^{-1}{F}_{p,0}.$ 可注意质量矩阵M₀是分块对角矩阵。因此，矩阵A₀同样是分块对角的并可逐单元估算。

在初始化之后，可执行代数粗化。这可通过令l，l＝1，...，L为整数，并假设代数系统可在下面形式的新精细网格Ω_h＝Ω_h，l-1上构造来完成。

$\left(\begin{array}{cc}{A}_{l-1}& C_{l-1}^T\\ C_{l-1}& 0\end{array}\right)\left(\begin{array}{c}{u}_{l-1}\\ {\mathrm{\&lambda;}}_{l-1}\end{array}\right)=\left(\begin{array}{c}{\hat{F}}_{u,l-1}\\ F_{\mathrm{\&lambda;},l-1}\end{array}\right).$

当然，该系统可针对l＝1来构造，针对其它数目l＝2，...，L构造，下面的规程可递归应用。自由度u_l-1和λ_l-1可分为两个集合：

$u_{l-1}=\left(\begin{array}{c}{u}_{l-1,\mathrm{\&Gamma;}}\\ u_{l-1,i}\end{array}\right),$ 和 ${\mathrm{\&lambda;}}_{l-1}=\left(\begin{array}{c}{\mathrm{\&lambda;}}_{l-1,\mathrm{\&Gamma;}}\\ {\mathrm{\&lambda;}}_{l-1,i}\end{array}\right).$

其中第二集合(子向量u_l-1，i和λ_l-1，i)包括对应关于Ω_h，l中的单元在内部的Ω_h，l-1的面的自由度。剩余自由度可被纳入第一集合中。然后，可进行下面两个步骤从而粗化网格。在第一步骤中，内部自由度u_l-1，i和λ_l-1，I被消除。在第二步骤中，自由度u_l＝R^Tu_l-1，Γ和λ_l＝R_λ ^Tλ_l-1，Γ的变换可被执行。结果，系统：

$\left(\begin{array}{cc}{A}_l& C_l^T\\ C_l& 0\end{array}\right)\left(\begin{array}{c}{u}_l\\ {\mathrm{\&lambda;}}_l\end{array}\right)=\left(\begin{array}{c}{\hat{F}}_{u,l}\\ F_{\mathrm{\&lambda;},l}\end{array}\right)$

在新粗糙网格Ω_H＝Ω_h，l上按照(u_l，λ_l)获得。应注意质量矩阵M_l可通过下面的公式逐单元(关于Ω_h，l)计算：

$M_l=A_l-B_l^T{\mathrm{\&Sigma;}}_l^{-1}{B}_l,$

其中矩阵∑_l是对角的。

在粗化之后，可引入项p_L。在重复粗化规程L次之后，代数系统：

$\left(\begin{array}{cc}{A}_L& C_L^T\\ C_L& 0\end{array}\right)\left(\begin{array}{c}{u}_L\\ {\mathrm{\&lambda;}}_T\end{array}\right)=\left(\begin{array}{c}{\hat{F}}_{u,L}\\ F_{\mathrm{\&lambda;},L}\end{array}\right),$

按照关联最粗糙网格Ω_H＝Ω_h，L的(u_L，λ_L)获得。可引入向量p_L从而获得相同的宏观复合混合系统：

$\left(\begin{array}{ccc}{M}_L& B_L^T& C_L^T\\ B_L& -{\mathrm{\&Sigma;}}_L& 0\\ C_L& 0& 0\end{array}\right)\left(\begin{array}{c}{u}_L\\ p_L\\ {\mathrm{\&lambda;}}_L\end{array}\right)=\left(\begin{array}{c}{F}_{u,L}\\ F_{p.L}\\ F_{\mathrm{\&lambda;},L}\end{array}\right).$

为求解上面的系统，分块对角矩阵M_L被求逆(尽管大部分的块具有尺寸5，但尖灭可导致尺寸4的块)。子向量u_L然后被消除。结果，在粗糙网格上以下面形式获得对称正系统

$\left(\begin{array}{cc}{S}_p& S_{\mathrm{p\&lambda;}}\\ S_{\mathrm{\&lambda;p}}& S_{\mathrm{\&lambda;}}\end{array}\right)\left(\begin{array}{c}{p}_L\\ {\mathrm{\&lambda;}}_L\end{array}\right)=\left(\begin{array}{c}{G}_{p,L}\\ G_{\mathrm{\&lambda;},L}\end{array}\right).$

求解该系统从而获得解对(p_L，λ_L)。然后，通量的自由度向量通过下面的公式重构

$M_L{u}_L=F_{u,L}-B_L^T{p}_L-C_L^T{\mathrm{\&lambda;}}_L$

最终，可在精细网格Ω_h＝Ω_h，0上恢复解。在该步骤，可假设三个一组的解(u_l，p_l，λ_l)对于某个l，l＝1，...，L已知。然后，可重回代数粗化规程从而在网格Ω_h，l-1上获得三个一组的解(u_l-1，p_l-1，λ_l-1)。可重复该规程L次从而获得关联最精细网格Ω_h＝Ω_h，0的期望的三个一组的解(u₀，p₀，λ₀)。

系统

在此讨论的技术可在计算装置上实施，例如在图14中图解的计算装置。图14图解示范计算机系统1400，在示范计算机系统1400上可实施执行本技术实施例的处理操作的软件。中央处理单元(CPU)1401耦合到系统总线1402。在实施例中，CPU 1401可以是任何通用CPU。只要CPU 1401(或示范系统1400的其它组件)支持如在此描述的本发明操作，那么本技术不受CPU 1401(或示范系统1400的其它组件)的架构约束。CPU 1401可根据实施例执行各种逻辑指令。例如，CPU1401可执行机器级指令，以便根据连同图1在上面描述的示范操作流程执行处理。作为特定例子，CPU 1401可执行机器级指令以便执行图1的方法。

计算机系统1400也可包括可以是SRAM、DRAM、SDRAM等的随机存取存储器(RAM)1403。计算机系统1400可包括可以是PROM、EPROM、EEPROM等的只读存储器(ROM)1404。RAM 1403和ROM1404保存用户和系统数据和程序，如在本领域中众所周知的。程序可包括存储在RAM 1404上的代码，该代码可用来根据本技术的实施例用均匀化混合有限元建模地质性质。

计算机系统1400也可包括输入/输出(I/O)适配器1405、通信适配器1414、用户接口适配器1408和显示适配器1409。I/O适配器1405、用户接口适配器1408和/或通信适配器1411可在某些实施例中使用户能够与计算机系统1400交互以便输入信息。

I/O适配器1405可连接总线1402到计算机系统1400的存储装置(若干)1406，例如硬盘驱动器、紧凑光盘(CD)驱动器、软盘驱动器、磁带驱动器、闪存驱动器、USB连接存储器等中的一个或更多。存储装置可在RAM 1403不足以满足和存储本技术的实施例的操作数据关联的存储需求时利用。例如，计算机系统1400的存储装置1406可用来存储这样的信息：计算网格、中间结果和组合的数据集，和/或根据本技术的实施例使用或生成的其它数据。

通信适配器1411适于耦合计算机系统1400到网络1412，网络1412可使信息能够经网络1412从系统1400输出和/或输入到系统1400，网络1412例如是互联网或其它广域网、局域网、公共或私有交换电话网、无线网络或上述任何组合。用户接口适配器1408耦合用户输入装置，例如键盘1413、指示装置1407和话筒1414和/或输出装置，例如(若干)扬声器1415到计算机系统1400。显示适配器1409由CPU 1401驱动从而控制显示装置1410上的显示，例如从而显示处于分析的目标区域的信息，例如根据某些实施例显示计算网格、储层或目标区域的生成的表示。

应该意识到本技术不限于在图14中图解的计算机系统1400的架构。例如，任何合适的基于处理器的装置可用于实施本技术的实施例的全部或一部分，包括而不限于个人计算机、膝上计算机、计算机工作站和多处理器服务器。此外，实施例可在专用集成电路(ASIC)或超大规模集成(VLSI)电路上实施。事实上，本领域技术人员可利用能够根据实施例执行逻辑操作的任何数目的合适结构。

尽管本技术可易受各种修改和可替换形式，但在上面讨论的示范实施例仅作为例子示出。然而，应再次理解本技术不意图限于在此公开的特别实施例。当然，本技术包括落入所附权利要求的真实精神和范围内的全部替换、修改和等效。

Claims (20)

1.一种使用处理器用均匀化混合有限元建模地质性质的方法，包含：

投影储层的特征到水平面上从而形成投影；

创造解析所述投影中期望特征的二维非结构化计算网格；

投影所述二维非结构化计算网格到边界面上以定义最精细计算网格；

生成至少一个更粗糙计算网格，其中所述至少一个更粗糙计算网格包含多个计算单元，并且所述多个计算单元的每个包括多个更精细单元；

生成关联所述多个计算单元的每个的多个计算面，其中所述计算面的每个包含多个更精细面；

使第一未知量关联所述多个计算单元的每个，并使第二未知量关联所述多个计算面的每个；

在所述最精细计算网格上推导出宏观复合混合有限元离散化；

通过粗化规程迭代从而将已知信息从所述最精细计算网格转移到最粗糙计算网格；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算单元的每个的所述第一未知量的每个获得值；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算面中每个的所述第二未知量的每个获得值；以及

通过恢复规程迭代从而将所述主未知量的所述值恢复到所述多个更精细单元的每个，并将所述次未知量的所述值恢复到所述多个更精细面的每个。

2.根据权利要求1所述的方法，其中投影所述储层的所述特征包含投影尖灭边界、断层线或井位置到所述水平面中。

3.根据权利要求2所述的方法，其中所述投影是非正交的和/或倾斜的。

4.根据权利要求1所述的方法，其中所述二维非结构化计算网格包含正方形、多边形、四边形或三角形或其任何组合。

5.根据权利要求1所述的方法，其中所述多个计算单元包含箱体、六角体、棱柱体、四面体、棱锥体或其任何组合。

6.根据权利要求1所述的方法，其中所述第一未知量对应所述储层的物理性质。

7.根据权利要求1所述的方法，其中所述第二未知量对应通量的法向分量。

8.根据权利要求1所述的方法，其中所述最精细计算网格逼近感兴趣层的边界面。

9.根据权利要求8所述的方法，其中物理性质在所述最精细计算网格上定义。

10.根据权利要求9所述的方法，其中所述物理性质包含流体压力、温度、渗透率、导热率或其任何组合。

11.根据权利要求1所述的方法，包含执行均匀化混合有限元规程以便在计算网格上求解扩散方程。

12.一种用均匀化混合有限元建模地质性质的系统，包含：

处理器；

包含数据库的存储介质，所述数据库包含储层数据；以及

包含代码的机器可读介质，所述代码经配置以引导处理器：

投影储层的特征到水平面上从而形成投影；

创造解析所述投影中期望特征的二维非结构化计算网格；

投影所述二维非结构化计算网格到边界面上以定义最精细计算网格；

生成至少一个更粗糙计算网格，其中所述更粗糙计算网格包含多个计算单元，并且所述多个计算单元的每个包含多个更精细单元；

生成关联所述多个计算单元的每个的多个计算面，其中所述计算面的每个包含多个更精细面；

使第一未知量与所述多个计算单元的每个关联，并使第二未知量与所述多个计算面的每个关联；

在所述最精细计算网格上推导出宏观复合混合有限元离散化；

通过粗化规程迭代从而将已知信息从所述最精细计算网格转移到最粗糙计算网格；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算单元每个的所述第一未知量的每个获得值；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算面每个的所述第二未知量的每个获得值；以及

通过恢复规程迭代从而恢复所述主未知量的所述值到所述多个更精细单元的每个，并恢复所述次未知量的所述值到所述多个更精细面的每个。

13.根据权利要求12所述的系统，进一步包含显示器，其中所述机器可读介质包含经配置在所述显示器上生成所述储层的图像的代码。

14.根据权利要求12所述的系统，其中所述储层数据包含静毛比、孔隙率、渗透率、压力、温度或其任何组合。

15.一种储层的烃管理的方法，包含：

生成储层的模型，该模型包含在非结构化计算网格中的多个均匀化混合有限元；

粗化所述非结构化计算网格从而在所述模型中形成多个更粗糙计算网格；

在最粗糙计算网格上评估对流扩散地下过程；

将结果从所述最粗糙计算网格转移到最精细计算网格；

从所述模型预测所述烃储层的性能参数；以及

使用所述预测的性能参数进行所述储层的烃管理。

16.根据权利要求15所述的方法，进一步包含：

投影所述储层的特征到水平面上从而形成投影；

创造解析所述投影中期望特征的二维非结构化计算网格；

投影所述二维非结构化计算网格到边界面上以定义所述最精细计算网格；

生成更粗糙计算网格，其中所述更粗糙计算网格包含多个计算单元，并且所述多个计算单元的每个包含多个更精细单元；

生成关联所述多个计算单元的每个的多个计算面，其中所述计算面的每个包含多个更精细面；

使第一未知量与所述多个计算单元的每个关联，并使第二未知量与所述多个计算面的每个关联；

在所述最精细计算网格上推导出宏观复合混合有限元离散化；

通过粗化规程迭代从而将已知信息从所述最精细计算网格转移到所述最粗糙计算网格；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算单元每个的所述第一未知量的每个获得值；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算面每个的所述第二未知量的每个获得值；以及

通过恢复规程迭代从而恢复所述主未知量的所述值到所述多个更精细单元的每个，并恢复所述次未知量的所述值到所述多个更精细面的每个。

17.根据权利要求15所述的系统，其中所述储层的所述烃管理包含烃采掘、烃生产、烃勘探、鉴别潜在烃资源、鉴别井位置、确定井注入率、确定井采掘率、鉴别储层连通度或其任何组合。

18.根据权利要求15所述的系统，其中所述性能参数包含生产率、压力、温度、渗透率、传递率、孔隙率、烃成分或其任何组合。

19.一种有形计算机可读介质，包含经配置而引导处理器执行下面操作的代码：

投影储层的特征到水平面上从而形成投影；

创造解析所述投影中期望特征的二维非结构化计算网格；

投影所述二维非结构化计算网格到边界面上以定义逼近所述边界面的最精细计算网格；

生成至少一个更粗糙计算网格，其中所述更粗糙计算网格包含多个计算单元，并且所述多个计算单元的每个包含多个更精细单元；

生成关联所述多个计算单元的每个的多个计算面，其中所述计算面的每个包含多个更精细面；

使第一未知量与所述多个计算单元的每个关联，并使第二未知量与所述多个计算面的每个关联；

在所述最精细计算网格上推导出宏观复合混合有限元离散化；

通过粗化规程迭代从而将已知信息从所述最精细计算网格转移到最粗糙计算网格；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算单元每个的所述第一未知量的每个获得值；

求解矩阵方程从而为所述最粗糙计算网格中所述多个计算面每个的所述第二未知量的每个获得值；以及

通过恢复规程迭代从而恢复所述主未知量的所述值到所述多个更精细单元的每个，并恢复所述次未知量的所述值到所述多个更精细面的每个。

20.根据权利要求19所述的有形机器可读介质，包含经配置而引导所述处理器显示储层表示的代码。

Images