2017, Vol. 38
油藏数值模拟技术在油气田开采开发中发挥着越来越重要的作用.如何建立可靠的油藏模型, 从而使模拟结果准确可信是油藏数值模拟前处理工作中的重要研究内容.油藏模型的可靠性取决于它采用的网格形式及其对地下地质情况(断层、不整合、尖灭等)的逼近程度.油藏网格的建立目标是在遵循地质特点的前提下尽量保持网格稀疏且正交, 以确保模拟过程的正确和高效[1].
油藏网格模型的研究起步较早, 出现过多种网格形式.Zakrevsky[2]对网格进行了较为全面的划分, 并且概述了网格剖分的一般流程.油藏模型网格主要分为结构化网格和非结构化网格.其中结构化网格一般是六面体形式, 以i, j, k方式有序组织, 如笛卡尔网格、角点网格等; 非结构化网格有四面体网格、PEBI网格.在诸多网格研究中, 以Pillar Grid为代表的角点网格曾是研究的热点, 并且应用也最为成熟[3].Pillar Grid的特点是网格走向可以沿着断层线、边界线或尖灭线, 也就是说网格可以扭曲.这样就克服了笛卡尔正交网格的不灵活性, 可以方便地模拟断层、边界和尖灭.随着技术的发展和油藏模拟精度的不断提高, Pillar Grid的缺陷逐渐暴露出来.首先, 由于Pillar几何位置的限制, 使用它建立复杂断层网时不得不简化断层; 其次, Pillar Grid的非正交性一方面给数值模拟过程中传导率的计算带来难度, 另一方面对结果精度也有影响.文献[4]提出一种网格生成方法.该方法的内容是在整体上采用笛卡尔网格形式, 通过切割六面体来模拟断层, 最终形成XY方向上正交、局部被切割的网格模型.这种网格能在表达断层形态的同时避免因扭曲而损失正交性, 但是在断层较多的模型中会产生大量多面体网格, 同样影响流体模拟运算.文献[5]分析了Pillar Grid的局限后, 提出一种更适用流体模拟的垂向竖直网格, 并且用实验证明相比于Pillar Grid, 它在网格粗化和流体模拟上有更好的效果, 但是文中并未给出这种网格的生成方法.
此外, 还有许多研究致力于非结构化网格.四面体网格生成方法成熟, 但是由于其结构化差, 计算太复杂, 基本未见有用它进行油藏数值模拟方面的报道.PEBI(perpendicular bisectors)网格灵活且局部正交, 能够降低网格走向对结果的影响.但是使用PEBI网格进行数值模拟求解还不够成熟, 应用还不广泛.
虽然针对油藏数值模拟网格形式的研究多样, 但在当前应用领域中, 最适用于模拟计算的仍然是相互正交的六面体网格.其原因有三:1) 相比于多面体或PEBI网格, 六面体网格连通面的数量更少, 因而需要求解较少的流体模拟方程, 这意味着计算效率的提高; 2) 六面体的正交性保证了数值模拟计算结果的稳定性; 3) 大多数商业流体模拟软件仍然只支持六面体网格.目前, 在生成方法上, 除了Pillar Grid和PEBI网格, 针对其他网格的研究较少.为此, 本文研究一种六面体网格模型生成方法, 其特点是在XY方向上规则, 采用阶梯形式逼近断层、边界等不连续面, 避免网格扭曲或因切割而产生多面体.同时基于该研究提出了一套从构造模型建立到油藏网格模型生成的工作流程.
1 工作流程网格生成的一般思路是离散化的方法:首先建立构造地质模型, 并将地质界面按沉积顺序划分序列, 分别离散化为多个四边形, 然后将相邻序列面的所有对应四边形连接成为六面体.然而在复杂地质构造情况下, 由于断层(正、逆)以及不整合面的存在, 在特定位置的地层序列会发生变化, 往往不能确定此处的网格属于哪一层.(在图 1a中, 地层S2, S3分别尖灭于S1和S4.在仅有地层面的情况下, 若以顶面为网格分层依据, 则G1, G2被认定属于同一层; 反之若以底面为依据, 则G2, G3被认定属于同一层.在图 1b中, F1, F2组成“Y”型断层, 在仅有断层面和地层面的情况下, 无论以网格顶面或是底面为分层依据, 顶底面都是断层的网格G1无法确定其地层属性).因此在生成网格之前, 有必要建立宏观拓扑模型, 表达地质对象之间的相互关系.
|
图 1 网格分层不确定性示意图 Fig.1 Uncertainty of grid stratigraphic division (a)－尖灭；(b)－断层. |
基于上述分析, 本文提出一套可行的油藏网格生成方法, 包括构造模型建立、宏观拓扑模型建立以及网格剖分(见图 2).该方法能处理包括逆断层、Y型断层和不整合面在内的复杂情况.在网格剖分的过程中, 这些不连续面与地层面等同处理, 生成的网格在不连续处正好是阶梯形式, 满足流体模拟的要求.后文将以不规则三角网为例, 对具体过程进行介绍.
|
图 2 油藏网格模型生成方法流程图 Fig.2 Flow chart of the grid generation method for the reservoir |
三维构造建模主要表达地下空间的几何形态和空间关系, 是建立三维地质模型的基础和核心部分之一.在本文的油藏网格生成方法中, 它是模型的基本骨架, 为网格的生成提供边界约束.因此, 构造建模直接决定油藏模型的准确程度.在现有的构造建模方法中, 复杂交互建模方法由于提供了很高程度的人工干预[6], 并且可将建模用户的地质专业经验融入到建模过程中, 能够处理任意复杂地质情况下三维构造模型的建立, 往往在建模生产中作为最终解决方案.本文的构造建模采取这种思路, 其主要步骤包括断层建模[7-9]和地层建模[10].
2.2 宏观拓扑模型如前文所述, 为了在网格生成过程中确定其地层属性, 需要在构造模型的基础上建立宏观拓扑模型.本文提出使用G-maps(generalized maps)建立宏观拓扑.
宏观拓扑[11](macroscopic-topology)表达地质界面的拓扑关系, 如断层之间的主次关系、断层地层之间的错断关系等, 其实质是使用数学模型来表达地质关系.G-maps是最早由Lienhardt引入的具有代数学基础的拓扑框架[12-13], 是一种表达拓扑结构的有效方式, 它能表达任意维度的流形对象.有许多学者研究使用G-maps表达拓扑结构, 但大都是针对微观数据结构[14], 几乎没有将它应用于三维地质模型宏观拓扑结构的研究.
在G-maps中, 点、边、面和体分别被称为0维单元、1维单元、2维单元和3维单元.它将所有几何单元划分成一系列抽象元素Dart, G-maps就是包含这些基本元素Dart的集合, 以及对这些Dart元素的对合运算αi.Dart可以理解为半边, 在一个二维的G-maps中, 可用三元组d=(F, E, V)来表示, 其中F, E, V分别代表半边所在的面片、边和节点.为了方便地实现拓扑关联信息的查找, 实现这样的对合运算:α0(d)代表映射到与d的面片和边相同但节点不同(即F, E相同, V不同)的相邻Dart; α1(d)代表映射到与d的面片和节点相同但边不同(F, V相同, E不同)的相邻Dart; α2(d)代表映射到与d的节点和边相同但面片不同(即V, E相同, F不同)的相邻Dart(见图 3a).基于这个思想, 有以下定义.
|
图 3 G-maps示意图 Fig.3 Sketch map of G-maps (a)－简单二维图形的2-G-maps；(b)－以d1起始获得的3种Dart集合. |
n-G-maps:其中G=(D, α0, α1, α2, …, αn), 由若干Dart组成的非空、有限集合D和n+1个αi组成，并且满足
● ∀i∈{0, …, n-2}, ∀j∈{i+2, …, n}, αi, αj是一个对合;
● 对于D上的每一个Dart d有
|
Orbit:orbit[α0, α1, α2, α3](d)描述了从Dart d开始, 通过连续按任意顺序使用α0, α1, α2, α3可以“获得”的所有Dart的集合.一个i-orbit=[{αj|j≠i, 0≤j≤n}](d), 在一个3-G-maps中, 1-orbit[α0, α2, α3](d)包含一个边的所有Dart; 2-orbit[α0, α1, α3](d)包含一个面的所有Dart; 3-orbit[α0, α1, α2](d)包含一个体的所有Dart(见图 3b).
Embedding:每个k-embedding与k-orbit相对应, 代表1个k维几何体, 例如0-embedding代表1个点, 1-embedding代表1条边.利用它可以将相关联的几何信息嵌入到拓扑模型中.
本文提出使用G-maps建立宏观拓扑模型的过程在上一步每层地层模型建立之后, 除了首层(按沉积序列最早的一层)需要创建一个3-G-maps之外, 每新生成一层地层需更新这个3-G-maps.在最终的拓扑模型中, 带地层属性的地质体作为3-embedding, 与三维几何体对应; 地质界面作为2-embedding, 与二维几何面对应.这样与地质界面求交得出的网格, 可以通过其对应二维几何面的拓扑转换操作找到与地质体对应的三维几何体, 就能确定该网格的地层属性.该方法具体步骤介绍如下:
步骤1 建立拓扑框架.搜索本层以及上一层地层面的边界E和这些边界的交点V, 这些边界在几何上由若干封闭环线组成, 交点将这些封闭环线分割为若干弧段.然后可以对这些点线进行抽象, 形成拓扑分界线集合E*和拓扑分界点集合V*, 完成该地层的拓扑框架.
步骤2 提取几何框架.根据构造建模所遵循的拓扑原则, 两个界面有接触关系, 则必有其中一个界面边界在另一个界面上.因此, 根据搜索出的边界线, 可以将地质界面在不连续处拓扑分离为一系列的几何面片集合S.
步骤3 创建本层3-G-maps.在步骤1建立的拓扑框架基础上, 以E*, V*为数据原型, 构建Dart和运算α, 形成本层的3-G-maps模型.
步骤4 嵌入几何信息.对步骤3的3-G-maps模型追踪k-orbit, 分别将边界线E、面片S以及该层地层属性信息作为k-embedding(k=1, 2, 3) 嵌入模型.若本层为初始层, 则步骤结束, 本层3-G-maps模型即目前全局3-G-maps; 否则, 进入步骤5.
步骤5 更新全局3-G-maps.为了把新一层地层的拓扑信息加入, 需要更新目前全局3-G-maps模型, 构造建模的拓扑一致原则保证这一步能顺利实施.更新内容主要是新加入Dart和几何嵌入对象, 为此实现针对3-G-maps的“缝合”操作.另外, 地层尖灭导致上层地层面被拓扑分离为新的面片, 上层拓扑模型的二维几何嵌入也要更新.这一步骤结束后, 若本层不是最后一层, 则在下一层地层建模后重新进入步骤1;否则, 方法结束, 输出最终的G-maps拓扑模型.
2.3 网格剖分如前文所述, 本方法中的网格剖分最终目的是生成六面体网格, 避免网格变形和切割对数值模拟带来的影响.其特点是在XY平面上规则, 对于不连续面的表达则是阶梯形式.基于这个目标, 本文提出基于“无限面”的网格剖分方法, 核心思想是将层面离散化为一系列切向“无限”平面, 使用这些切平面和纵向的侧面组合为六面体网格.该方法将地层面和断层面等同对待, 还能自动解决包含逆断层、Y形断层在内的复杂问题.
首先说明一般性思路:在确定了合适的坐标系和水平分辨率后, 可在XY平面上建立规则四边形网格, 以它为骨架在竖直方向延伸形成网格列(以下简称Column), 然后从宏观拓扑模型中取出地质界面按地层序列与Column求交, 得到横向顶底面，将Column划分为一个个网格.然而, 考虑到地质界面部分落入某个Column的情况, 网格顶底面可能“缺失”或者“二义”.本文分别以“切分”和“粘合”两步解决这个问题.
在上述一般思路中, 地质界面是否落入Column由它与Column的4个网格柱求交决定, 如图 4所示, 该网格柱称为初始柱(以下简称InitPL).为防止二义性出现, 取4个InitPL的几何中心, 定义为中心柱(以下简称CenPL), 地质界面是否落入Column由CenPL与地质界面求交唯一决定.若未与当前层面求得交点, 则认为该层不在Column内; 反之, 求出层面在此交点的切面.因为切面为平面，在三维空间可以无限延伸，本文称为无限面.无限面因为“无限”性质, 它与InitPL必定存在交点, 也就可用无限面在交于Column内的部分逼近本层网格横向面.
|
图 4 网格划分算法示意图 Fig.4 Schematic diagram of grid generation algorithm |
层面的离散化, 必然使得网格失去了横向连续性, 因此需要把同一层网格横向面“粘合”.粘合思路较多, 最简单的是记录下每个Column上的交点信息, 取每个InitPL上的同层交点的几何平均值作为这个InitPL与本层的唯一交点.
3 实例论证为验证本文方案正确可行, 笔者选取某区域的局部模块, 在自主开发的软件系统中进行了实验.该模块的数据经过解译后有两层地层, 都被1条逆断层错断.按上述流程, 分别进行构造建模、拓扑建模和网格剖分.图 5所示结果表明该方法能生成符合要求的网格模型.
|
图 5 应用效果图 Fig.5 Application effect (a)－构造模型；(b)－提取地质界面边界获得的宏观拓扑结构; (c)－宏观拓扑存储的几何要素(地质界面, 地质体)；(d)－网格模型. |
油藏数值模拟运算都是基于网格系统进行的, 一套好的网格系统可以保证模拟结果的正确性, 能满足复杂油藏数值模拟的要求.本文从这一点出发, 分析了多种网格的应用现状及各自优缺点, 总结出油藏数值模拟对网格的要求, 并针对性地对网格生成方法进行了研究, 用实例验证方法的可行性.本文的不足之处在于没有研究对该网格的粗化方法, 然后对粗化后的网格模型进行油藏数值模拟实验.这将是下一步研究的重点.
| [1] | Abdou M K, Pham H D, Al-Aqeell A S. Impact of grid selection on reservoir simulation[J]. Journal of Petroleum Technology, 1993, 45(7): 664–669. DOI:10.2118/21391-PA |
| [2] | Zakrevsky K E. Geological 3D modelling[M]. Houten: EAGE Publisher, 2011: 1-100. |
| [3] | Fremming N P.3D geological model construction using a 3D grid[C]//ECMOR VIII-8th European Conference on the Mathematics of Oil Recovery.Catania, 2002. |
| [4] | Haussling H J, Coleman R M. A method for generation of orthogonal and nearly orthogonal boundary-fitted coordinate systems[J]. Journal of Computational Physics, 1981, 43(2): 373–381. DOI:10.1016/0021-9991(81)90129-7 |
| [5] | Gringarten E J, Haouesse M A, Arpat G B, et al.Advantages of using vertical stair step faults in reservoir grids for flow simulation[C]//SPE Reservoir Simulation Symposium.Woodlands, Texas:Society of Petroleum Engineers, 2009. |
| [6] | Pan M, Li Z L, Gao Z B, et al. 3-D geological modeling—concept, methods and key techniques[J]. Acta Geologica Sinica (English Edition), 2012, 86(4): 1031–1036. DOI:10.1111/j.1755-6724.2012.00727.x |
| [7] |
李兆亮, 潘懋, 杨洋, 等.
三维复杂断层网建模方法及应用[J]. 北京大学学报(自然科学版), 2015, 51(1): 79–85.
( Li Zhao-liang, Pan Mao, Yang Yang, et al. Research and application of the three-dimensional complex fault network modeling[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2015, 51(1): 79–85. ) |
| [8] | Tertois A L, Mallet J L. Editing faults within tetrahedral volume models in real time[J]. Geological Society, 2007, 292(1): 89–101. DOI:10.1144/SP292.5 |
| [9] | Cherpeau N, Caumon G, Lévy B. Stochastic simulations of fault networks in 3D structural modeling[J]. Comptes Rendus Geoscience, 2010, 342(9): 687–694. DOI:10.1016/j.crte.2010.04.008 |
| [10] | Euler N, Sword J C H, Dulac J C.A new tool to seal a 3D earth model:a cut with constraints[M]//SEG Technical Program Expanded Abstracts 1998.Tulsa: Society of Exploration Geophysicists, 1998:710-713. |
| [11] | Caumon G, Collon-Drouaillet P, De Veslud C L C, et al. Surface-based 3D modeling of geological structures[J]. Mathematical Geosciences, 2009, 41(8): 927–945. DOI:10.1007/s11004-009-9244-2 |
| [12] | Lienhardt P. Topological models for boundary representation:a comparison with n-dimensional generalized maps[J]. Computer-Aided Design, 1991, 23(1): 59–82. DOI:10.1016/0010-4485(91)90082-8 |
| [13] | Lienhardt P. N-dimensional generalized combinatorial maps and cellular quasi-manifolds[J]. International Journal of Computational Geometry & Applications, 1994, 4(3): 275–324. |
| [14] | Levy B, Mallet J L.Cellular modelling in arbitrary dimension using generalized maps[R].Nancy:ISA-GOCAD (Inria-Lorraine/CNRS), 1999. |