在非结构化栅格上对地下过程进行建模

摘要

本发明的实施例包含形成棱柱体栅格以及使用该棱柱体栅格和混合有限元分析来求解对流-扩散问题。可以通过在模型的平面上提供三角形网格来形成所述棱柱体栅格。接着对所述网格进行粗化以使较不满意的网格变大。接着将所述粗化的网格进行投影以形成所述棱柱体栅格。接着所述栅格的每一单元被分配多个自由度。所述栅格的混合有限元分析产生矩阵，接着求解该矩阵以得到所述对流-扩散问题的解。

说明书

相关申请的交叉引用

该申请要求2007年12月14日提交的题为“MODELING-SUBSURFACEPROCESSES ON UNSTRUCTURED GRID”的美国临时专利申请61/007,761的利益，通过引用将其整体合并到本文。

技术领域

该申请一般地涉及计算机建模，具体地说，涉及用非结构化栅格对地下过程进行建模。

背景技术

在地质勘探中，期望获得关于地球表面下存在的各种岩层和结构的信息。这类信息可以包括测定的地质层、密度、孔隙度、组成成分等。接着，这个信息用来对地下盆地进行建模，即使用获得的数据预测碳氢化合物储藏的位置并且帮助碳氢化合物的提取。

非结构化栅格对于对复杂的地质结构(例如地下盆地)中的物理过程进行建模具有许多有吸引力的特性。这种栅格也可以用于其他工业中，例如航空工业或汽车工业。所关心的盆地或域可以被建模或表示为堆叠在一起的不同厚度的层的集合。地质层可沿着垂直表面或倾斜表面断裂并且退化，生成所谓的尖灭(pinch-out)。尖灭被定义为具有接近零厚度的地质层的部分。栅格应该考虑这种复杂性以产生地质层的良好模型。非结构化栅格提供比结构化栅格更好的模型。非结构化栅格可以包括由它们的顶点定义的多面体元素或单元的集合并且具有完全任意的拓扑图。例如，栅格的顶点可以属于许多个单元并且每一单元可以具有任何数目的边或面。

可以通过对流-扩散类型的数学方程来描述许多物理地下过程。这类过程的示例可以是多孔介质中的流体流、温度分布和/或压力分布。用于石油勘探的重要过程是温度分布或热建模。热建模包含从地壳下的岩浆并且穿过沉积层和源岩运动的热量。源岩是包含在石油和其他碳氢化合物材料的地层中的岩石。石油和/或其他碳氢化合物材料从源岩中被排出并且迁移到其他地方。碳氢化合物的质量由源岩和它们周围区域上承受的温度和压力条件确定。质量也受源岩和其当前位置之间的迁移路径的温度和压力条件影响。因此，盆地在其整个历史中的压力和温度条件是重要的。

为了更精确地对过程进行建模，重要的是不仅对主变量(例如压力或温度)进行建模，而且对在任何给定的表面之上的主变量的通量或能量、流体的流速等进行建模。有多种已知的方法用于对这些过程进行建模，例如有限差分、有限体积或有限元方法。在这些方法中，在考虑物理过程的地方，域被栅格覆盖。接着，通过以下方法来在栅格上逼近域：在栅格单元的指定位置处引入称为自由度的一组未知量并且为每一位置导出代数方程，所述方程将在那个位置的自由度与其他自由度连接起来。对于上面所提到的不同方法，导出这类方程的方式和自由度的位置是不同的，但是所有这些方法具有共同的特征，即，它们仅包含主变量，例如温度或压力。

为了计算通量，感兴趣的研究人员将首先使用上面所描述方法中的一个来计算所期望的主变量，接着使用数值微分来计算主变量的通量。在规则栅格(例如矩形或平行六面体栅格)上精确的所有现有数值微分方法在非结构化栅格上是不精确的并且在计算上非常耗时，尤其是如果考虑物理过程的域是高度不均匀的或异质的。此外，使用有限差分方法求解对流-扩散问题的方法需要笛卡尔(Cartesian)栅格，因此不可应用于不得不使用非结构化栅格的许多地下应用。能够对复杂几何形状进行建模的有限元方法不具有局部守恒属性并且不能够被应用在许多地下过程中。相反地，有限体积方法是局部守恒的并且能够被应用在局部正交的非结构化栅格的子集上。然而，当非结构化栅格不具有局部正交属性时，有限体积方法提供不精确的解。因此，从上面所提到的所有三种方法中，没有一个可应用于描述以非结构化栅格建模的盆地中的地下对流-扩散过程。

有另一数学方法来同时逼近主未知量和它们的通量，叫做混合有限元方法，其在以下文献中描述：F.Brezzi和M.Fortin的“Mixed and hybrid finite elementmethods”，Springer Verlag，柏林，1991年。这种方法被证明在存在不均匀介质的情况下是局部质量守恒的、精确的，并且提供对主未知量和通量两者的精确逼近。直到最近，混合有限元方法才能够被直接应用到由非结构化多面体栅格覆盖的域，所述非结构化多面体栅格对于地下应用是非常常见的。在以下文献中提出了针对在任意多面体栅格上的扩散类型方程的混合有限元方法的新变型：Yu.Kuznetsov和S.Repin的“New mixed finite element method on polygonal andpolyhedral meshes”，Russian Journal of Numerical Analysis and MathematicalModeling，第18卷，261-278页，2003年。

发明内容

本发明涉及对在经过地质时期形成的盆地中的能量传递和/或压力分布提供精确模型的系统和方法。在任何给定的时间，盆地被表示为堆叠在一起的不同厚度的层的集合。在盆地中的某些位置，层的厚度退化到零，形成尖灭。本发明的实施例使用棱柱体网格和混合有限元分析来对盆地中的各种过程进行建模，这些过程包括能量传递(例如热能量)和压力。因此，本发明的实施例求解主未知量(例如温度或压力)和从未知量(例如温度通量或压力通量)两者。下述方面中的一个或更多个可以用来提供经过地质时期形成的盆地中的能量传递和/或压力分布(例如物理过程)的精确模型。该模型可以用来解释现代储油层，并且继而依赖于该模型，从而基于模型的模拟结果来控制碳氢化合物开采活动。基于从由下述方面中的一个或多个产生的模拟盆地模型中解释的结果，可以控制碳氢化合物的开采，例如，可以控制来自地面设施的开采速率，可以有策略地打井，和/或一般性地表征储油层。

在一个一般方面，一种用于在计算机上对物理区域进行建模的方法，其中所述物理区域包括多个地层，所述方法包括：接收定义所述物理区域的至少一个物理特性的数据；在所述物理区域的模型的平面上提供三角形网格，其中所述网格包括多个单元；以非均匀方式对所述三角形网格进行粗化，使得较不满意的单元较大；以及将粗化的三角形单元在与所述物理区域中的所述平面正交的方向上进行投影以形成棱柱体栅格，其中根据所述地层，粗化的三角形网格的每一单元被分割成子单元。

这个方面的实现可以包括下述特征中的一个或更多个。例如，粗化可以包括将所述数据投影在平面上以及使用所述数据来确定哪些单元是较不满意的。所述模型可以包括对所述物理区域中的物理特征建模的建模特征，以及其中使用可以包括将优先值分配给每一单元，其中基于每一单元是否接近建模特征和该建模特征的类型来确定该值。提供三角形网格可以包括在所述平面上提供矩形网格；以及沿着至少一个对角线分割所述矩形网格的每一单元。粗化可以包括通过消除两个相邻三角形的公共边来合并所述两个相邻三角形。所述棱柱体栅格可以包括多个棱柱体单元、多个棱锥体单元和多个四面体单元。所述方法可以用来对所述物理区域中的物理过程的至少一个通量进行建模，所述方法还包括为每一子单元中的通量分配多个自由度；将混合有限元分析应用到每一子单元以产生矩阵；以及求解该矩阵以确定所述区域中的物理过程的通量。

分配可以包括对于每一单元，为所述物理过程分配一个自由度；以及为所述单元的每一面分配另一自由度。应用可以包括使用恒定划分法来形成有限元空间。所述物理过程可以是对流-扩散过程。所述物理过程可以是温度和压力中的一个并且所述物理区域是地下地质盆地。所述物理过程可以包含碳氢化合物材料的形成。所述物理过程可以包含碳氢化合物材料的运动。可以从来自传感器的信息导出所述数据，所述传感器测量所述物理区域的至少一个物理特性。

在另一一般方面，一种用于在计算机上对物理过程和该物理过程的通量进行建模的方法，其包括形成对物理区域建模的非结构化棱柱体栅格，其中所述物理过程在所述物理区域内发生并且所述棱柱体栅格包括多个单元；对于每一单元，为所述物理过程和所述通量分配多个自由度；将混合有限元分析应用到每一单元以产生矩阵；以及求解该矩阵以确定所述区域中的所述物理过程和所述通量。

这个方面的实现可以包括下述特征中的一个或更多个。例如，形成可包括在所述物理区域的模型的平面上提供三角形网格，其中所述网格包括多个单元；以非均匀方式对所述三角形网格进行粗化，使得较不满意的单元较大；以及将粗化的三角形单元在与所述物理区域中的所述平面正交的方向上进行投影以形成所述棱柱体栅格。所述棱柱体栅格可以包括多个棱柱体单元、多个棱锥体单元以及多个四面体单元。分配可以包括对于每一单元，为所述物理过程分配一个自由度；以及对于每一单元，为所述单元的每一面分配另一自由度。应用可以包括使用恒定划分法来形成有限元空间。所确定的物理过程和通量可以用来影响所述物理区域中的改变。所述物理过程可以是温度和压力中的一个并且所述物理区域是地下地质盆地。

在另一一般方面，一种具有计算机可读介质的计算机程序产品，所述计算机可读介质具有在其上记载的计算机程序逻辑，所述计算机程序逻辑用于对物理区域中的物理过程和该物理过程的通量进行建模，所述计算机程序产品包括用于形成对所述物理区域进行建模的非结构化棱柱体栅格的代码；用于将混合有限元分析应用到所述棱柱体栅格以产生矩阵的代码；以及用于求解该矩阵从而确定所述区域中的物理过程和通量的代码。

这个方面的实现可以包括下述特征中的一个或更多个。例如，用于形成的代码可以包括这样的代码：用于在所述物理区域的模型的平面上提供三角形网格的代码，其中所述网格包括多个单元；以非均匀方式对所述三角形网格进行粗化，使得较不满意的单元是较大的；以及将粗化的三角形网格在与所述物理区域中的所述平面正交的方向上进行投影以形成所述棱柱体栅格。所述棱柱体栅格可以包括多个单元，以及所述用于应用的代码可以包括对于每一单元，为所述物理过程分配一个自由度；对于每一单元，为所述单元的每一面分配另一自由度；以及使用恒定划分法来形成有限元空间。所述用于求解的代码可以包括使用预调的共轭梯度分析来求解所述矩阵。

本发明的实施例通过将某些或大多数地质和几何特征，例如尖灭边界投影到水平平面中来工作。注意，投影可以是非正交的或倾斜的。接着，本发明的实施例在那个平面上生成分解所有期望的特征的非结构化栅格。注意到，栅格可以由多边形、四边形、三角形或其组合组成。接着，本发明的实施例将获得的栅格投影回到所有层的所有边界表面上，从而构造棱柱体栅格。所述棱柱体栅格可以包括多个单元，这些单元可以是棱柱体、四面体形状、棱锥体或其组合。注意，非结构化棱柱体栅格逼近所有层的边界表面。

接着本发明的实施例可以通过以下步骤来工作：对于主未知量，在单元中心处为每一单元关联一个自由度，并且对于通量的法向分量，在面中心处为所述单元的每一面关联一个自由度。接着本发明的实施例使用混合有限元方法，例如Yu.Kuznetsov和S.Repin的方法来使问题离散化。空间离散化产生稀疏矩阵方程。接着本发明的实施例可以求解该矩阵方程，以得到主未知量和在所述单元的各面处的通量的法向分量两者。因此，本发明的实施例在没有极大地扩展需要求解的未知量的数目的情况下，提供更精确的建模。

本发明的实施例可以提供覆盖所述物理区域的水平平面的三角形网格来形成所述棱柱体栅格。接着可以以非均匀方式对所述网格进行粗化，使得较不满意的单元是较大的，而使满意的单元保留更细的形式。接着将粗化的网格在所述物理区域中的垂直方向上进行投影以形成所述棱柱体网格。

为了可以更好地理解后面的本发明的详细描述，前面已经相当广泛地概述了本发明的特征和技术优势。在下文中将描述本发明的另外特征和优势，其构成本发明的权利要求的主题。本领域技术人员应意识到，公开的概念和具体的实施例可以容易地用作修改或设计其他结构的基础，其他结构用于实现本发明的相同目的。本领域技术人员也应该认识到，这类等效的构造不偏离随附的权利要求中陈述的本发明的精神和范围。当结合附图考虑时，从下述描述中将更好地理解被认为是本发明的特性的新颖特征(关于其组织架构和操作方法两者)以及其他目标和优势。然而，应清楚地明白，提供每一附图是仅为了解释和描述的目的，无意作为本发明的限制的定义。

附图说明

为了更彻底地理解本发明，现在结合附图参考下述描述，在附图中：

图1描绘根据本发明的实施例的被分割成层的域。

图2A和2B描绘根据本发明的实施例的域和用矩形网格覆盖的域。

图3描绘根据本发明的实施例的非均匀粗化的三角形栅格。

图4描绘根据本发明的实施例形成的不同类型的单元。

图5描绘根据本发明的实施例形成的3D棱柱体栅格。

图6描绘本发明的实施例使用的四面体单元。

图7描绘本发明的实施例使用的棱锥体单元。

图8A-8D根据本发明的实施例描绘本发明的实施例使用的、被分成三个四面体的棱柱体单元。

图9描绘根据本发明的实施例的分割相邻单元的独立面；

图10描绘根据本发明的实施例的形成棱柱体栅格的方法；

图11描绘根据本发明的实施例的求解矩阵的方法；以及

图12描绘适合使用本发明的计算机系统的框图。

具体实施方式

本发明的实施例对于对地下油田进行建模是有用的。在此所描述的实施例的示例可以参考这类油田。然而，实施例可以用来对包含其他材料和/或过程的的其他域进行建模。例如，可以包含其他碳氢化合物材料，例如煤炭。本发明的实施例对于采矿或掘坑道是有用的。本发明的实施例可以用于其他域类型，例如大气，并且对于对天气、温度和/或污染进行建模是有用的。另一域可以是海洋，并且实施例可以用来测量声音、温度、含盐度和/或污染度。可以使用本发明的实施例来对任何类型的分层域进行建模。可以使用本发明的实施例来对通过对流-扩散过程而运动的任何类型的材料进行建模。可以使用本发明的实施例来对域或材料中存在的任何类型的通量进行建模。

如前面所述的，可以将本发明的实施例应用到任何对流-扩散过程。下面是3D对流-扩散类型方程的示例：

$-▿·(K▿p)+c·p=f$在Ω中           (1.1)

这里p是未知函数(称为压力)，K＝K(x)是扩散张量，c是非负函数，f是源函数，以及是有界计算域。假设K是均匀正定矩阵，并且域Ω的边界被分割成两个不重叠的集合Γ_D和Γ_N。

用边界条件补充方程(1.1)

p＝g_D    在Γ_D上

$(K▿p)·n+\sigma ·p=g_N$在Γ_N上      (1.2)

这里n是Γ_N的外向单位法向向量，σ是非负函数，并且g_D和g_N是给定函数。假设方程(1.1)-(1.2)具有唯一解。

可以由等价一阶系统代替偏微分方程(1.1)-(1.2)：

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

$▿·u+c·p=f$在Ω中         (1.3)

p＝g_D  在Γ_D上

-u·n+σ·p＝g_N     在Γ_N上       (1.4)

方程(1.3)-(1.4)是方程(1.1)-(1.2)的混合公式表示。注意，以此方式，可以同时逼近主未知量p及其通量u。

如前面所述的，本发明的实施例可以在不同的域工作。因此，设G是R²中具有规则成形的边界的域，例如分段光滑并且段之间的角度大于0。设计算域Ω被定义为如下：

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

这里Z_min(x，y)和Z_max(x，y)是光滑表面。

设m是正整数并且z＝Z_i(x，y)，i＝0，...，m，是定义在G上的单值连续函数，使得

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

Z_i-1(x，y)≤Z_i(x，y)       在中$i=\overline{1,m}$

Z_m(x，y)≡Z_max(x，y)在中

这些函数定义地质层之间的界面。换句话说，计算域Ω可以被分成m个子域(带或层)，其被定义如下：对于所有i＝1，...，m，

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

图1描绘了将计算域Ω100分割成多个子域或层101-107的示例。注意，图1以分解图描绘了不同的层，然而对于计算，不需要将层隔开。还注意到，假设子域Ω_i满足锥条件，这里子域的边界不具有奇点(零角度等)，并且此外，所有集合

$P_i=\{(x,y)\in \bar{G}:Z_{i-1}(x,y)=Z_i(x,y)\}$

由有限数目的多边形组成。用表示对应的集合P_i的边界。

图1描绘了盆地的不同地层。可以通过各种技术例如地层分析和/或地震反演，使用传感器来测量盆地的各种特性来确定用来形成图1的不同层的数据。

可以以下述方式将集合P_i的边界投影到平坦平面。对于来自的任何给定的点(x，y，z)，投影的点具有坐标(x，y，0)。所有这类点组成像图2A那样的闭合线的集合，其用来生成平面三角测量，在图2B和图3中示出了这种示例。

微分方程(1.3)-(1.4)的变分混合公式可以被写成如下：找到p∈L₂(Ω)和λ∈L₂(Γ_N)，使得对于所有q∈L₂(Ω)和μ∈L₂(Γ_N)

$\underset{\Omega }{\int }\left(K^{-1}u\right)·\mathrm{vdx}-\underset{\Omega }{\int }p(▿·v)\mathrm{dx}+\underset{{\Gamma }_N}{\int }\lambda (v·n)\mathrm{ds}=-\underset{{\Gamma }_D}{\int }{g}_D(v·n)\mathrm{ds}$

$\left(\begin{array}{cc}-\underset{\Omega }{\int }(▿·u)\mathrm{qdx}-\underset{\Omega }{\int }c·\mathrm{pqdx}& =-\underset{\Omega }{\int }\mathrm{fqdx}\end{array}\right)---\left(1.5\right)$

$\left(\begin{array}{cc}\underset{{\Gamma }_N}{\int }(u·n)\mathrm{\mu ds}& -\underset{{\Gamma }_N}{\int }\mathrm{\sigma \lambda \mu ds}=-\underset{{\Gamma }_N}{\int }{g}_N\mathrm{\mu ds}\end{array}\right)$

这里

${\hat{H}}_{\mathrm{div}}\left(\Omega \right)=\{v:v\in {\left[L_2\left(\Omega \right)\right]}^3,▿·v\in L_2\left(\Omega \right),\underset{\partial \Omega }{\int }{|v·n|}^2\mathrm{ds}<\infty \}.$

注意，λ是压力函数p＝p(x)在Γ_N上的约束。在这个公式中，所有边界条件是自然边界条件。

在Γ_N上σ＝0的情况下，可以以不同的形式将变分混合公式写为如下：找到u·n＝-g_N(在Γ_N上)和p∈L₂(Ω)，使得对于所有v·n＝0(在Γ_N上)和q∈L₂(Ω)

$\underset{\Omega }{\int }\left(K^{-1}u\right)·\mathrm{vdx}-\underset{\Omega }{\int }p(▿·v)\mathrm{dx}=-\underset{{\Gamma }_D}{\int }{g}_D(v·n)\mathrm{ds}$

$\underset{\Omega }{\int }(▿·u)\mathrm{qdx}+\underset{\Omega }{\int }c·\mathrm{pqdx}=\underset{\Omega }{\int }\mathrm{fqdx}---\left(1.6\right)$

在下述分析中，考虑方程(1.5)，然而也可以将结论应用到方程(1.6)而不失一般性。

本发明的实施例使用棱柱体栅格，其在对对流-扩散地下过程进行建模方面提供许多有吸引力的特性。在许多情况下，域可以被表示为堆叠在一起的不同厚度的层的集合。棱柱体栅格满意地表示了堆叠层中的分片地质结构和非结构化几何特征。通过地质统计信息的后处理来提供2D几何数据。通常，数据是与2D细矩形栅格的节点或单元关联的材料属性。注意，可能有数百万个节点。节点中存在材料数据意味着其是计算节点，而不存在材料数据则意味着其是局外节点(node-outlier)。计算节点的集合定义了计算域。

地质层之间的界面是地质统计数据，并且可能互相相交，导致拓扑上不正确的情形。底部界面是最低的地质层的底部边界，并且也由地质统计数据表示。地质层可以沿着垂直表面断裂并且退化。尖灭被定义为厚度模不大于用户定义的阈值δ≥0的地质层的部分。断层折线被定义为底部地质界面和断层的相交处。

为了简化算法的描述，在本发明的实施例中要使用的栅格应该满足某些非常自然的要求。目标棱柱体栅格应该是1D栅格和2D三角测量的逻辑积，所述1D栅格和2D三角测量的节点形成不具有局外节点的地质统计矩形栅格的子集，其中节点密度等于所关心的用户指定区域中原始栅格中的节点密度。除此之外，在用户定义的断层和井以及自动检测到的尖灭附近必须精炼(改进)三角测量。此外，2D栅格中三角形的数目不应该大于用户定义的数目。关于棱柱体栅格，棱柱体的侧面不得不逼近地质界面并且形成形状规则的2D三角测量。

下面所描述的过程是可以用来构造棱柱体栅格的过程的示例。注意，可以使用其他过程。此外，那种类型的棱柱体栅格不限制混合有限元方法。首先是为了将尖灭投影在表示井的底部地质界面、断层折线和点上，对所述2D规则三角测量的产生进行改进。接下来，将2D三角测量投影在由函数Z_i(x，y)，i＝1，...，m所定义的表面上，以形成得到的3D棱柱体栅格。在下面的段落中进一步描述这个过程。

为了形成2D规则三角形栅格，示例性过程可以开始于矩形栅格。给定节点在x-和y-方向上的节点坐标，产生覆盖域G的矩形一致网格G^h，其由具有带材料数据的四个节点中的至少一个的单元组成。例如，图2A描绘域200，图2B描绘覆盖该域200的矩形栅格201。不失一般性，假设网格G^h由正方形组成，即在x-和y-方向上的网格尺寸相等，h_x＝h_y＝h。

接下来，根据本发明的实施例，每一矩形单元被其对角线分割成两个三角形。下面描述可以用来形成三角形的一个过程，注意，可以使用其他过程。可以根据下述规则来进行两个可能的对角线之间的选择。设每一矩形单元被分配一整数，该整数等于其节点的最小x-和y-索引的和。对于具有偶数的单元，分割对角线具有带最小x-和y-索引的单元的节点。对于具有奇数的单元，选择另一对角线。对于具有给定的材料数据分布的给定的节点集合，上述过程唯一地指定三角测量。改变对角线的方向减少了栅格朝向的问题。作为形成三角形的替换过程，可以通过使用两个对角线来将每一矩形单元分割成四个三角形。

将产生的三角测量投影在底部地质界面上，如公式(1.2)那一段中所描述的。注意，在域G中可以定义感兴趣的区域ω_i，这里栅格的修改是不必要的或不期望的。

设P_i^h表示G^h的矩形元素的最大子集，其属于P_i。如果没有属于P_i的G^h的元素，但是有属于P_i的G^h的顶点，则这个顶点被认为属于P_i^h。那么集合被定义为

${\partial P}^h=\underset{i=1}{\overset{m}{\cup }}{{\partial P}_i}^h$

这被称为“尖灭”投影。根据该定义，“尖灭”投影是属于G^h的边和顶点的子集。

接下来，根据本发明的实施例，将不同的优先级分配给三角形。将优先级或整数标记分配给细栅格的每一三角形。优先级的值控制粗化过程。

在开始时，将零优先级分配到所有三角形。对于其闭合线与断层相交的三角形，它们的优先级被改变为1。为了找到与断层相交的三角形，可以使用下述方法。首先，从断层三角测量中提取三角形，所述断层三角测量与最底部地质界面相交。其次，检查每一提取的三角形是否与细栅格三角形相交。对于其闭合线与尖灭相交的三角形，它们的优先级被改变为2。这些三角形被定义为违反下述条件的三角形：在所有三角形节点中，地质层的厚度或大于δ或小于-δ。对于其闭合线包含井点的三角形，它们的优先级被改变为3。对于属于用户定义的感兴趣的区域ω_i的三角形，它们的优先级被改变为4。对于满足若干条件的三角形，它们的优先级被改变为最大优先级值。

在已经分配优先级之后，可以对三角形栅格进行非均匀粗化。栅格的细部分可以具有大量的同样小的三角形。这些区域是更期望的，因为它们包含更多的信息，包括感兴趣的地质特征，例如井、断层、尖灭和/或用户指示为期望的特征。栅格的粗化部分与栅格的较细部分不是同样满意的。栅格可以具有粗化范围，这里最粗化的区域指示具有较不满意或不满意质量的部分，并且不具有粗化的区域指示最期望的区域。在最粗化的区域和不粗化的区域之间的粗化区域指示具有某些期望方面的区域。

粗化是一系列合并三角形的过程。例如，可以通过消除它们的公共边来将两个三角形合成一个。这个过程包括两个阶段。首先，标记特定三角形用于进行粗化。其次，对它们进行粗化。应该注意到，栅格一致性可能引起对未标记三角形的粗化。每一粗三角形继承了两个合并的三角形的最大优先级。除了优先级外，每一三角形还被分配表示级别的另一整数。初始细栅格的任何三角形具有级别1。可以将粗化应用到相同级别j的一对三角形，并且得到具有级别j+1的更粗的三角形。

下面是根据本发明的实施例的粗化过程的一个示例。注意，可以使用其他过程。

粗化过程可以被描述为下述循环：

1)设置k＝1。

2)由具有零优先级的三角形形成集合M，零优先级三角形的粗化不会引起对具有非零优先级的其他三角形进行粗化。

3)如果M为空，则到6。

4)对来自M的三角形进行粗化。

5)到2。

6)如果新栅格中三角形的数目不大于用户定义的阈值Nt_usr，则停止。

7)由非零优先级不大于k的三角形形成集合M，其粗化不会引起对优先级大于k的其他三角形进行粗化。

8)对来自M的三角形进行粗化。

9)如果k≤3，则设置k＝k+1，否则设置k＝1。

10)到1。

输出三角形栅格具有与投影的输入矩形网格的特定节点重合的节点和感兴趣区域中的细三角形，以及朝着井点、断层折线和尖灭进行改进的三角形。图3描绘了上述栅格粗化过程的示例。得到的栅格300描绘了非均匀区域，其具有最粗三角形301、粗化三角形302和细三角形303。注意，细三角形可能具有某些粗化或没有粗化。注意，在图3中，通过使用上面所描述的两对角线方法来形成三角形。

在粗化之后，可以形成3D棱柱体栅格。设e_h是G的三角测量中的三角形，并且a^(k)＝(a_x^(k)，a_y^(k))，k＝1，2，3是e_h的顶点。考虑在R³中的三条垂直线(x，y)＝a^(k)，k＝1，2，3，，并且用a^(k，i)，k＝1，2，3表示它们与表面z＝Z_i(x，y)，i＝0，...，m的交界处。于是，具有位于相邻表面上的顶点的任何多面体或是垂直棱柱体(所有6个顶点是不同的)、或是棱锥体(两个顶点重合，例如对应的顶点a^(k)，k＝1，2，3，属于)、或四面体(两对顶点重合，即集合的两个顶点a^(k)，k＝1，2，3，属于P^h)。图4描绘了具有三种类型的棱柱体的3D棱柱形栅格的部分，所述三种类型的棱柱体即垂直棱柱体401、棱锥体402和四面体403。换句话说，棱锥体是一个边消失的棱柱体，而四面体是两个边消失的棱柱体。

对中的所有三角形执行上述操作将提供域Ω的分区Ω_h。在特定的单元中，通过分段三角形表面z＝Z_i^(h)(x，y)来逼近每一表面z＝Z_i(x，y)，所述分段三角形表面由棱柱体的顶部(底部)三角形以及棱锥体和四面体的特定面组成。图5描绘了3D棱柱体栅格500的示例。

在完成3D棱柱体栅格之后，可以对栅格进行混合有限元分析。在前面的部分中指出，栅格Ω_h包括元素{E_k}，其或是垂直棱柱体，或是棱锥体，或是四面体。为了针对方程(1.5)用公式表示混合有限元(MFE)方法，应该定义空间L₂(Ω)和L₂(Γ_N)的有限元子空间。

有限元空间包括函数p_h，其在每一栅格单元上是恒定的。有限元空间包括函数λ_h，其在Ω_h中的栅格单元E_k与边界部分Γ_N的每一交界处上是恒定的。这些交界处可以是四边形或是三角形。

混合有限元方法中的一个问题是设计空间的有限元子空间V_h。为了计算效率，应该仅考虑那些有限元向量函数，它们在相邻单元E_k和E_l(k＞l)之间的界面Γ_kl上以及在单元E_k与边界Γ_N的交界处上具有恒定的法向分量。空间的有限元子空间的维数等于不同界面{Γ_kl}和的总数目。可以基于在以下文献中描述的“恒定划分(div-constant)”法来构造这个有限元空间：Yu.Kuznetsov和S.Repin的“New mixed finite element method on polygonal andpolyhedral meshes”，Russian Journal of Numerical Analysis and MathematicalModeling，第18卷，261-278页，2003年。

对于四面体单元T，有限元空间V_h|_T与经典的最低阶Raviart-Thomas有限元空间RT₀(T)一致(参见F.Brezzi和M.Fortin的“Mixed and hybrid finite elementmethods”，Springer Verlag，柏林，1991年)。有限元向量值函数w_h∈RT₀(T)具有四个自由度(DOF)，即

$w_h\left(x\right)=\underset{i=1}{\overset{4}{\Sigma }}{w}_i{\phi }_i\left(x\right)$

这里φ_i(x)是与四面体T的面γ_i(i＝1，2，3，4)关联的基向量函数。

用γ_j表示的四面体T的面与顶点A_j相对，即，γ₁是面A₂A₃A₄，γ₂是面A₃A₄A₁，γ₃是面A₄A₁A₂，γ₄是面A₁A₂A₃。设n_j是在面γ_j上的外向单位向量，并且h_j是从顶点A_j到面γ_j上的垂直线长度，j＝1，2，3，4。在图6中示出了这种四面体单元600。

在四面体T上最低阶Raviart-Thomas元素的空间可以被定义为

RT₀(T)＝span{φ₁，φ₂，φ₃，φ₄}

这里基向量函数φ_i满足条件并且δ_ij是kronecker符号，如果i＝j，则δ_ij等于1，否则其等于0。

简单明了的计算显示了基函数可以被明确地定义为

${\phi }_i\left(x\right)=\frac{1}{h_i}(x-x^{\left(i\right)})\equiv \frac{\left|{\gamma }_i\right|}{3\left|T\right|}(x-x^{\left(i\right)})$

这里x⁽ⁱ⁾是顶点A_i的坐标，i＝1，2，3，4。

对于棱锥体单元P，当单元P∈Ω_h是四边形棱锥体时，使用“恒定划分”法来构造有限元空间V_h|_P。

用γ_j，j＝1，2，3，4，5来表示棱锥体P的面，即，γ₁是面A₂A₃A₅，γ₂是面A₁A₃A₄，γ₃是面A₁A₂A₄A₅，γ₄是面A₁A₂A₃，以及γ₅是面A₃A₄A₅。在图7中示出了这种棱锥体单元700。

为了描述“恒定划分”法，将棱锥体P分成两个四面体T₁和T₂。可以以两种不同的方式来完成，并且两种方式都提供下面所描述的可工作的算法。因此，不失一般性，假设棱锥体P被划分成两个四面体T₁＝A₁A₂A₃A₄和T₂＝A₂A₃A₄A₅。设n_i表示到面γ_i的单位外向法向向量，i＝1，2，3，4，5，并且n₆表示到四面体T₁和T₂之间的界面γ₆＝A₂A₃A₄的、从T₁指向T₂的单位法向向量。

在每一四面体上，构造向量函数的经典的最低阶Raviart-Thomas空间，确定基函数(k＝1，2)的集合，并且P中的向量场u_h被定义为如下：

这个表达式显示出向量函数u_h在四面体T₁和T₂中的每一个上是线性的，其属于空间H_div(P)，并且满足所要求的条件

$u_h{|}_{{\gamma }_j}·n_j=u_j,j=\mathrm{1,2,3,4,5}.$

为了确定在界面γ₆上通量的法向分量的未知值u₆，下述条件应该是可操作的

$▿·u_h\equiv \mathrm{const}$(常数)在P中             (2.2)

从u_h的定义(2.1)中获得T₁和T₂中u_h的散度的表达式。将散度算子应用于(2.1)的两边，得到

$▿·u_h{|}_{T_1}=u_6▿·{\phi }_1^{\left(1\right)}+u_2▿·{\phi }_2^{\left(1\right)}+u_3▿·{\phi }_3^{\left(1\right)}+u_4▿·{\phi }_4^{\left(1\right)}---\left(2.3\right)$

以及

$▿·u_h{|}_{T_2}=-u_6▿·{\phi }_4^{\left(2\right)}+u_5▿·{\phi }_1^{\left(2\right)}+u_3▿·{\phi }_2^{\left(2\right)}+u_1▿·{\phi }_3^{\left(2\right)}---(2.4)$

由于因此得到

$u_6=\frac{u_5▿·{\phi }_1^{\left(2\right)}+u_3▿·{\phi }_2^{\left(2\right)}+u_1▿·{\phi }_3^{\left(2\right)}-u_2▿·{\phi }_2^{\left(1\right)}-u_3▿·{\phi }_3^{\left(1\right)}-u_4▿·{\phi }_4^{\left(1\right)}}{▿·{\phi }_1^{\left(1\right)}+▿·{\phi }_4^{\left(2\right)}}---\left(2.5\right)$

可以以不同的方式得到在假设(2.2)下u₆的值。第一，应用Stokes公式

$\underset{P}{\int }▿·u_h\mathrm{dx}=\underset{\partial P}{\int }{u}_h·\mathrm{nds}$

由下式确定的值

$▿·u_h=\frac{u_1\left|{\gamma }_1\right|+u_2\left|{\gamma }_2\right|+u_3\left|{\gamma }_3\right|+u_4\left|{\gamma }_4\right|+u_5\left|{\gamma }_5\right|}{\left|P\right|}$

这里|γ_i|是对应的面γ_i的面积，|P|是棱锥体的体积。第二，从(2.4)确定u₆的值：

$u_6=\frac{u_5▿·{\phi }_1^{\left(2\right)}+u_3▿·{\phi }_2^{\left(2\right)}+u_1▿·{\phi }_3^{\left(2\right)}-▿·u_h}{▿·{\phi }_4^{\left(2\right)}}---\left(2.6\right)$

在这里，重新构造在棱锥体P上的基函数其满足条件i，j＝1，2，3，4，5。

由于在棱锥体P的面上的基函数的法向分量的值是已知的，因此其法向分量u₆⁽ⁱ⁾，i＝1，2，3，4，5在内部面γ₆上的值可以从公式(2.5)和(2.6)中得到。因此，由以下给出基函数的明确的表达式：

对于具有三角形底的棱柱体单元，设∏∈Ω_h。使用“恒定划分”法来构造有限元空间V_h|_∏。

用γ_j来表示棱柱体∏的面，即，γ₁是面A₂A₃A₅A₆，γ₂是面A₁A₃A₄A₆，γ₃是面A₁A₂A₄A₅，γ₄是底面A₁A₂A₃，以及γ₅是顶面A₄A₅A₆。在图8A中示出了这种棱柱体单元800。

为了根据实施例应用“恒定划分”法，将棱柱体∏分成三个四面体T₁、T₂和T₃。在图8B-8D中示出了分割的四面体801、802、803。不失一般性，假设棱柱体∏被分成四面体T₁＝A₁A₂A₃A₆、T₂＝A₁A₄A₅A₆和T₃＝A₁A₂A₅A₆。在图8B-8D中示出了分割的四面体801、802、803。

设n_i，i＝1，2，3，4，5表示到面γ_i的外向单位法向向量，n₆表示到四面体T₁和T₃之间的界面γ₆＝A₁A₂A₆的、从T₁指向T₃的单位法向向量，并且最后，n₇表示到四面体T₂和T₃之间的界面γ₇＝A₁A₅A₆的、从T₂指向T₃的单位法向向量。在每一四面体上构造经典的最低阶Raviart-Thomas有限元空间，确定基函数k＝1，2，3的集合，并且∏中的向量场u_h被定义为如下：

根据这个表达式，向量函数u_h在四面体T₁、T₂和T₃中的每一个上是线性的，属于空间H_div(∏)，并且满足所要求的条件

$u_h{|}_{{\gamma }_j}·n_j=u_j,j=\mathrm{1,2,3,4,5}.$

对于分别在附属界面γ₆和γ₇上未知法向向量u₆和u₇的自然选择来自于以下条件：在∏上u_h的散度是恒定的，即

$▿·u_h=\mathrm{const}$在∏上         (2.8)

以类似于针对棱锥体单元所讨论的方式得到未知值u₆和u₇。首先，使用Stokes公式来确定在∏上的值：

$\underset{\Pi }{\int }{\mathrm{div}u}_h\mathrm{dx}=\underset{\partial \Pi }{\int }{u}_h·\mathrm{nds}$

因此

$▿·u_h=\frac{u_1\left|{\gamma }_1\right|+u_2\left|{\gamma }_2\right|+u_3\left|{\gamma }_3\right|+u_4\left|{\gamma }_4\right|+u_5\left|{\gamma }_5\right|}{\left|\Pi \right|}---2.9$

这里|γ_i|是对应的面γ_i的面积并且|∏|是棱柱体的体积。其次，u₆和u₇的值计算如下：

$u_6=\frac{▿·u_h-u_1▿·{\phi }_1^{\left(1\right)}-u_2▿·{\phi }_2^{\left(1\right)}-u_4▿·{\phi }_4^{\left(1\right)}}{▿·{\phi }_3^{\left(1\right)}}---\left(2.10\right)$

以及

$u_7=\frac{▿·u_h-u_5▿·{\phi }_1^{\left(2\right)}-u_2▿·{\phi }_3^{\left(2\right)}-u_3▿·{\phi }_4^{\left(2\right)}}{▿·{\phi }_2^{\left(2\right)}}---\left(2.11\right)$

在这里，重新构造在棱柱体∏上的基函数其满足下述条件：i，j＝1，2，3，4，5。

由于在棱柱体∏的面上的基函数的法向分量的值是已知的，因此可以通过公式(2.10)和(2.11)来分别确定法向分量u₆⁽ⁱ⁾和u₇⁽ⁱ⁾，i＝1，2，3，4，5在内部面γ₆和γ₇上的值。因此，由下面给出基函数的明确的表达式：

强调“恒定划分”法的下述区别是合情合理的。将栅格单元分割成四面体不依赖于其相邻单元的分割，例如，两个相邻单元的公共的四边形面可以以不同的方式从不同的边进行划分。例如，图9描绘了两个相邻单元901、902，其具有以不同的方式分割的邻接的面903、904。与四面体栅格相比，这个特征简化了棱柱体栅格的产生和其上的边界值离散化问题。

对于一般的棱柱体单元，上面所描述的针对具有三角形底的棱柱体单元的方法可以被推广到覆盖一般的棱柱体∏。每一棱柱体可以被独立地分成四面体并且使用在前面的部分中描述的“恒定划分”法来构造有限元空间V_h|_∏。

上面的段落描述了对于通量的法向分量，单元的每一个面在面中心处具有一个自由度。对于每一主未知量，每一单元在单元中心也被分配关联的一个自由度。例如，(多个)主未知量可以包括周围的温度和/或流体压力。

下述段落描述了将混合有限元(MFE)分析的杂化变型用在3D棱柱体栅格的单元上，其中引入了与棱柱体的面关联的额外的自由度。这些自由度允许将一个单元的通量与另一单元的通量隔开。这些自由度在数学文献中被称为拉格朗日(LaGrange)乘子。由于原始的变量(例如一个棱柱体单元的单元中心的温度(或压力)和在边界上的通量)脱离任何其他单元的温度和通量，引入额外的自由度允许简化数值问题的结构。因此，可以消除那些未知量，这允许减少未知量的数目。

根据上面提供的用于自由度的定义，可以引入MFE方法如下：找到u_h∈V_h，p_h∈L_h以及使得对于所有v∈V_h，q∈L_h以及

$\underset{\Omega }{\int }\left(K^{-1}{u}_h\right)·\mathrm{vdx}-\underset{\Omega }{\int }{p}_h(▿·v)\mathrm{dx}+\underset{{\Gamma }_N}{\int }{\lambda }_h(v·n)\mathrm{ds}=-\underset{{\Gamma }_D}{\int }{g}_D(v·n)\mathrm{ds}$

$\left(\begin{array}{cc}-\underset{\Omega }{\int }(▿·u_h)\mathrm{qdx}-\underset{\Omega }{\int }c·p_h\mathrm{qdx}& =-\underset{\Omega }{\int }\mathrm{fqdx}\end{array}\right)---\left(3.1\right)$

$\left(\begin{array}{cc}\underset{{\Gamma }_N}{\int }(u_h·n)\mathrm{\mu ds}& -\underset{{\Gamma }_N}{\int }\sigma {\lambda }_h\mathrm{\mu ds}=-\underset{{\Gamma }_N}{\int }{g}_N\mathrm{\mu ds}\end{array}\right)$

有限元问题(3.1)得到线性代数方程系

$A\left(\begin{array}{c}\bar{u}\\ \bar{p}\\ \bar{\lambda }\end{array}\right)=\left(\begin{array}{c}{\bar{g}}_D\\ \bar{f}\\ {\bar{g}}_N\end{array}\right)---\left(3.2\right)$

其中鞍点矩阵

$A=\left(\begin{array}{ccc}M& B^T& C^T\\ B& -D& 0\\ C& 0& -\Sigma \end{array}\right)---\left(3.3\right)$

这里M＝M^T是正定矩阵，并且D＝D^T和∑＝∑^T或是正定矩阵或是半正定矩阵。可以显示出，系(3.2)具有唯一解。

很好地开发出了用于具有对称鞍点矩阵的代数系统的迭代方法。然而，起因于多面体栅格上MFE离散化的、用于鞍点矩阵的高效预调技术仍然是关注点。对于高效预调，对称正定矩阵是更好的对象。可以通过使用混合有限元问题(3.1)的杂化来将系(3.2)变换为具有对称正定矩阵的等价系。在接下来段落中，这个方法被描述为求解问题(3.1)的优选方式。

对于混合有限元分析的杂化，设E_k是Ω_h中的栅格单元，并且V_h^(k)和L_h^(k)分别是有限元空间V_h和L_h在E_k上的约束。此外，生成新的有限元空间Λ_h，其是在栅格单元之间的界面Γ_kl上以及在栅格单元与Γ_N的边界部分交界处上定义的函数λ_h＝λ_h(x)的空间。在每一界面上，函数λ_h∈Λ_h等于常量。

为了引入混合杂化有限元(MHFE)问题，不得不定义两个额外的有限元空间以及许多双线性形式和线性泛函。两个新的有限元空间是

${\hat{V}}_h=\underset{k=1}{\overset{n}{\Pi }}{V}_h^{\left(k\right)}$和${\hat{L}}_h=\underset{k=1}{\overset{n}{\Pi }}{L}_h^{\left(k\right)}$

这里n是Ω_h中单元E_k的数目。注意到空间V_h^(k)中的任何一个的维数最多是5，并且每一L_h^(k)的维数等于1。

对于元素u，v∈V_h，p，q∈L_h以及λ，μ∈Λ_h，引入下述双线性形式：

$a(u,v)=\underset{k=1}{\overset{n}{\Sigma }}\underset{E_k}{\int }\left(K^{-1}{u}_k\right)·v_k\mathrm{dx}$

$b(v,p)=-\underset{k=1}{\overset{n}{\Sigma }}\underset{E_k}{\int }{p}_k▿·v_k\mathrm{dx}\equiv -\underset{k=1}{\overset{n}{\Sigma }}{\hat{p}}_k\underset{E_k}{\int }▿·v_k\mathrm{dx}$

$c(v,\lambda )=\underset{l>k}{\overset{n}{\underset{k=1}{\Sigma }}}\underset{{\Gamma }_{\mathrm{kl}}}{\int }\lambda (v_k·n_{\mathrm{kl}})\mathrm{ds}-\underset{l<k}{\overset{n}{\underset{k=1}{\Sigma }}}\underset{{\Gamma }_{\mathrm{lk}}}{\int }\lambda (v_k·n_{\mathrm{lk}})\mathrm{ds}+\underset{k=1}{\overset{n}{\Sigma }}\underset{{\partial E}_k\cap {\Gamma }_N}{\int }\lambda (v_k·n)\mathrm{ds}$

$d(p,q)=\underset{k=1}{\overset{n}{\Sigma }}\underset{E_k}{\int }c·p_k{q}_k\mathrm{dx}\equiv \underset{k=1}{\overset{n}{\Sigma }}{\hat{c}}_k{\hat{p}}_k{\hat{q}}_k$

$\sigma (\lambda ,\mu )=\underset{k=1}{\overset{n}{\Sigma }}\underset{{\partial E}_k\cap {\Gamma }_N}{\int }\mathrm{\sigma \lambda \mu ds}\equiv \underset{k=1}{\overset{n}{\Sigma }}{\sigma }_{\mathrm{kN}}{\hat{\lambda }}_N{\hat{\mu }}_N$

这里u_k，和是p_k，的值，是λ，μ∈Λ_h的值，n_kl是到Γ_kl的、从E_k指向E_l(k＜l)的单位法向向量，并且

${\hat{c}}_k=\underset{E_k}{\int }\mathrm{cdx},$${\sigma }_{\mathrm{kN}}=\underset{{\partial E}_k\cap {\Gamma }_N}{\int }\mathrm{\sigma ds}$

此外，定义下述线性泛函

$l_D\left(v\right)=-\underset{k=1}{\overset{n}{\Sigma }}\underset{{\partial E}_k\cap {\Gamma }_D}{\int }{g}_D(v_k·n)\mathrm{ds}$

$l_f\left(q\right)=-\underset{k=1}{\overset{n}{\Sigma }}\underset{E_k}{\int }\mathrm{fqdx}$

$l_N\left(\mu \right)=-\underset{k=1}{\overset{n}{\Sigma }}\underset{{\partial E}_k\cap {\Gamma }_N}{\int }{g}_N\mathrm{\mu ds}$

这里和μ∈Λ_h。

借助上述定义，有限元问题(3.1)的等价混合杂化公式变为如下：找到和λ_h∈Λ_h，使得对于所有的和μ∈Λ_h：

a(u_h，v)+b(v，p_h)+c(v，λ_h)＝l_D(v)

b(u_h，q)-d(p_h，q)＝l_f(q)     (3.4)

c(u_h，μ) -σ(λ_h，μ)＝l_N(μ)

有限元问题(3.4)得到线性代数方程系

$A·\left(\begin{array}{c}\bar{u}\\ \bar{p}\\ \bar{\lambda }\end{array}\right)=\left(\begin{array}{c}{\bar{g}}_D\\ \bar{f}\\ {\bar{g}}_N\end{array}\right)---(3.5)$

其中鞍点矩阵

$A=\left(\begin{array}{ccc}M& B^T& C^T\\ B& -D& 0\\ C& 0& -\Sigma \end{array}\right),$这里

是具有对称正定子矩阵M_k，k＝1，...，n的块对角矩阵，D是对角正定或半正定矩阵，以及∑是对角半正定矩阵。通过(3.4)中的线性泛函来定义右手边子向量的分量。

矩阵A具有非常有用的表达式这里

$A_i=\left(\begin{array}{ccc}{M}_i& B_i^T& C_i^T\\ B_i& -D_i& 0\\ C_i& 0& -{\Sigma }_i\end{array}\right)---\left(3.6\right)$

是用于单元E_i的局部鞍点矩阵以及N_i是对应的组合矩阵(assemble matrix)。

相应地，系(3.5)的右手边可以被写为如下：

$\left(\begin{array}{c}{\bar{g}}_D\\ \bar{f}\\ {\bar{g}}_N\end{array}\right)=\underset{i=1}{\overset{n}{\Sigma }}{N}_i\left(\begin{array}{c}{\bar{g}}_{D,i}\\ {\bar{f}}_i\\ {\bar{g}}_{N,i}\end{array}\right)$

重要的是，发现可以通过将局部混合有限元离散化应用于下述问题来获得矩阵A_i以及子向量和

$u_i+K▿p_i=0$在E_i中

$▿·u_i+c·p_i=f_i$在E_i中

p_i＝g_D    在上

-u_i·n+σp_i＝g_N    在上

$u_i·n_i^k=0$在上，如果

这里Γ_i^k，k＝1，...s_i是单元E_i的面。

(2.6)中的矩阵M_i、B_i、C_i、D_i和∑_i的重要属性可以被展现在内部单元E_i上，即设单元E_i具有s_i个面，那么是对称正定矩阵，则

$B_i=-\left(\begin{array}{cccc}\left|{\Gamma }_i^1\right|& \left|{\Gamma }_i^2\right|& ...& \left|{\Gamma }_i^{s_i}\right|\end{array}\right)\in R^{1\times s_i}$

$C_i=\mathrm{diag}\left(\begin{array}{cccc}\left|{\Gamma }_i^1\right|& \left|{\Gamma }_i^2\right|& ...& \left|{\Gamma }_i^{s_i}\right|\end{array}\right)\in R^{s_i\times s_i}$

$D_i=c_i{V}_{E_i}\in R^{1\times 1}$

这里V_Ei是单元E_i的体积并且对于内部单元矩阵∑_i＝0。

设于是这个属性适用于来自(3.6)的任何A_i。

相关地，注意主变量和i＝1，...，n仅在单个单元内被连接。所以，可以容易排除这些未知量：

$\bar{u}=M^{-1}({\bar{g}}_D-C\bar{\lambda }-B\bar{p})---\left(3.7\right)$

以及

$\bar{p}={(B^T{M}^{-1}B+D)}^{-1}(B^T{M}^{-1}{\bar{g}}_D-\bar{f}-B^T{M}^{-1}C\bar{\lambda })---\left(3.8\right)$

由于矩阵M、B和D的结构，矩阵B^TM^-1B+D是对角矩阵，因此它是可逆的。使用关系(3.7)和(3.8)来将系(3.5)变换为系：

$S\bar{\lambda }=\bar{\xi }---\left(3.9\right)$

这里

S＝C^TM^-1C-C^TM^-1CB(B^TM^-1B+D)^-1B^TM^-1C+∑

以及

$\bar{\xi }=C^T{M}^{-1}{\bar{g}}_D-C^T{M}^{-1}B{(B^T{M}^{-1}B+D)}^{-1}(\bar{f}+B^T{M}^{-1}{\bar{g}}_D)-{\bar{g}}_N$

矩阵S被称为“浓缩矩阵(condensed matrix)”。这个矩阵是对称的和正定的，除了在Neumann边界条件的情况下，S是半正定的，但是其具有简单的核恒定(kernel constant)向量。这个矩阵在本质上是全局的，其将所有的节点或单元连接在一起。可以同时求解大的线性方程系。

可以将任何迭代方法应用到求解具有该矩阵的线性方程系(3.9)。一个方法是预处理的共轭梯度方法(PCG)，然而可以使用其他的方法。注意，在半正定性的情况下，应该在与核正交的子空间中执行PCG。在求解系(3.9)之后，可以分别使用方程(3.8)和(3.7)来逐个元素地局部恢复主未知量和

矩阵S也可以被表示为这里

$S_i=C_i^T{M}_i^{-1}{C}_i-C_i^T{M}_i^{-1}{C}_i{B}_i{(B_i^T{M}_i^{-1}{B}_i+D_i)}^{-1}{B}_i^T{M}_i^T{C}_i+{\Sigma }_i$

以及是对应的组合矩阵。(3.9)的右手边具有类似的表达式。

注意，本发明的实施例可以以单个主未知量操作，例如温度或压力，及其相关的通量。其他实施例可以以多于一个的主未知量操作。

根据本发明的各种实施例，上面概述的各种过程和方法可以被组合在一个或更多个不同的方法中，用在一个或更多个不同的系中，用在一个或更多个不同的计算机程序产品中。

例如，一个示例性方法1000可以用于形成棱柱体栅格，如图10中所示。使用正交投影将感兴趣的地质特征和几何特征(例如尖灭边界、断层线或井的位置)投影到水平面中(1001)。产生细的符合矩形的网格，该网格覆盖与在其上提供材料数据的细栅格相同大小的投影域的所有特征(1002)。矩形栅格被分割成三角形(1003)。在栅格上的各种线和点可以表示例如断层线和井位置。以非均匀方式对三角形进行粗化(1004)。期望保持某些地质或几何特征附近的细三角测量，但是远离这些特征具有较粗分辨率将允许更容易的分析。这种栅格仅由三角形组成。将粗化的栅格垂直投影在所有层的所有边界表面上以形成棱柱体栅格(1005)。这种栅格将包含可以是三角棱柱体、四面体或棱锥体的单元。以这种方式建立的非结构化棱柱体栅格逼近所有层的边界表面。注意，在对流-扩散地下问题中，输入数据与数百万个节点关联。

另一方法1100可以用于求解针对地质盆地的对流-扩散问题，如图11中所示。形成对盆地建模的栅格(1101)。注意，可以通过图10中所示的方法1000来形成栅格，或可以使用另一方法。该方法对于主未知量在单元中心处为每一栅格单元关联一个自由度，以及对于通量的法向分量在面中心处为单元的每一面关联一个自由度(1102)。使用混合有限元法来分析具有关联的自由度的栅格(1103)。这个分析产生稀疏矩阵方程。接着该方法可以求解矩阵方程，得到(多个)主未知量和(多个)未知量的通量在单元的面处的法向分量两者(1104)。

注意到，可以以硬件、软件和/或固件和/或其任何组合来实施在此描述的功能中的任何一个。当以软件实施时，本发明的元素本质上是执行必要任务的代码段。程序或代码段可以被存储在计算机可读介质中或通过计算机数据信号发送。“计算机可读介质”可以包括能够存储或传递信息的任何介质。计算机可读介质的示例包括电子电路、半导体存储器器件、ROM、闪存、可擦除ROM(EROM)、软盘、压缩盘CD-ROM、光盘、硬盘、光纤介质等。计算机数据信号可以包括能够在传输介质(例如电子网络信道、光纤、空气、电磁、RF链接等)上传播的任何信号。可以经由计算机网络(例如因特网、内部网等)下载代码段。

图12示出了适于使用本发明的计算机系统1200。中央处理器(CPU)1201被耦合到系统总线1202。CPU 1201可以是任何通用CPU，例如英特尔奔腾(IntelPentium)处理器。然而，本发明不被CPU 1201的架构所限制，只要CPU 1201支持在此描述的发明操作。总线1202被耦合到随机存取存储器(RAM)1203，其可以是SRAM、DRAM或SDRAM。ROM 1204也被耦合到总线1202，其可以是PROM、EPROM或EEPROM。正如本领域中熟知的，RAM 1203和ROM 1204保存用户和系统的数据和程序。

总线1202也被耦合到输入/输出(I/O)控制卡1205、通信适配器卡1211、用户接口卡1208和显示卡1209。I/O适配器卡1205将存储设备1206(例如硬盘驱动器、CD驱动器、软盘驱动器、磁带驱动器中的一个或多个)连接到计算机系统。I/O适配器1205可以被连接到打印设备，其允许系统打印信息(文件、照片、文档等)的纸质副本。注意，打印设备可以是打印机(例如喷墨打印机、激光打印机等)、传真机或复印机。通信卡1211适于将计算机系统1200耦合到网络1212，其可以是电话网络、局域网(LAN)和/或广域网(WAN)、以太网和/或因特网中的一个或多个。用户接口卡1208将用户输入设备(例如键盘1213和指针设备1207)耦合到计算机系统1200。用户接口卡1208也可以将声音经由扬声器输出到用户。由CPU 1201驱动显示卡1209以控制在显示设备1210上的显示。

虽然已经详细地描述了本发明及其优势，但是应明白，在不偏离由随附的权利要求所限定的本发明的精神和范围的情况下，在此可以进行各种改变、替代和变换。此外，无意将本申请的范围限制在说明书中所描述的过程、机器、制造、物质的组成成分、装置、方法和步骤的特定实施例。正如本领域任何一个技术人员从本发明的公开中容易理解的，根据本发明可以使用当前现有的或后来开发的、执行与在此所描述的对应的实施例实质上相同的功能或实现实质上相同的结果的过程、机器、制造、物质的组成成分、装置、方法或步骤。因此，随附的权利要求意在将这种过程、机器、制造、物质的组成成分、装置、方法或步骤包括在它们的范围内。