流体有限体积仿真的方法和系统

摘要

本发明提供一种建模碳氢化合物储层的方法，其包括使用多级混合多尺度有限体积（MMMFV）方法在精细非结构化网格上推导可计算网格。推导可计算网格包括计算用于压力的第一代数多级基础函数，构建互相作用区域，以及生成主网格。用于速率近似的第二代数多级基础函数被计算并且主网格被离散。使用可计算的网格仿真碳氢化合物储层。至少部分基于仿真的结果在非瞬态计算机可读介质内产生物理碳氢化合物储层的数据表示。

说明书

相关申请的交叉引用

本申请要求申请号为No.61/330,012美国临时专利申请的权益，该 申请提交于2010年4月30日，标题为METHOD AND SYSTEM FOR  FINITE VOLUME SIMULATION OF FLOW，其全部内容包括在此作为 参考。

技术领域

本技术的示例性实施例涉及提供多孔介质中流体的有限体积仿真 的方法和系统。

背景技术

本部分意图介绍可以与本技术的示例性实施例关联的本领域各方 面。该讨论据信帮助提供框架从而促进本技术的特别方面的更优理解。 因此，应理解本节应据此阅读并且不必需作为现有技术的承认来阅读。

现代社会非常依赖为燃料和化学原料使用碳氢化合物。碳氢化合 物一般在一般称为储层的地下岩层中发现。从储层去除碳氢化合物取 决于岩层的许多物理性质，尤其例如含碳氢化合物岩石的渗透度、碳 氢化合物流过岩层的能力，以及碳氢化合物存在的百分比。

经常地，称为“仿真模型”的数学模型用来仿真碳氢化合物储层 以便定位碳氢化合物并最优化碳氢化合物的生产。仿真模型是一类计 算流体动力学仿真，其中支配流过多孔介质和被连接设施网络的多相、 多成分流体流动的一组偏微分方程（PDE）可以逼近并求解。这是其 中最优化特别碳氢化合物生产策略的迭代的、时间步进的过程。

储层仿真模型的性质，例如渗透度或孔隙度，经常是高度非均匀 的并且可以变化。变化在从最小到可以与储层大小比较的最大尺度的 全部长度尺度中。使用非常精细的栅格离散化从而捕捉非均匀性的这 样模型的计算机仿真是计算上非常昂贵的。然而，忽视非均匀性可以 导致错误结果。

为实现合理的计算性能，储层性质经常放大尺度。例如，调和化 技术可以用来在较粗糙的仿真栅格上定义储层性质。在不同的变体中 该技术以合理的成功广泛用于学术界和工业界。然而，放大尺度可以 具有众多缺点。例如，放大尺度可以不为带有不可分离尺度的问题良 好运作，如在下面更详细讨论。进一步地，放大尺度可以不完全捕捉 储层流动中的全局特征。在复杂流动建模时，放大尺度也可以具有评 估误差的困难。

非均匀或多尺度现象可以分为两类：可分离现象和不可分离现象。 对于多孔介质，具有介质性质可以在其上变化的多尺度是普遍的。例 如，如果特征跨储层连续变化，那么储层中的尺度可以是不可分离的。 放大尺度方法可以用可分离尺度求解模型，但可以不能正确求解不可 分离模型。

除不可分离现象之外，具有可以延伸穿过储层的巨大部分的特征 如沟渠、长裂缝和断层。这些可以称为全局特征和信息。当前的放大 尺度方法可以丢失可以对流动仿真显著的全局特征的影响，即可以没 有放大尺度算法可以用来生成更复杂模型的信息。因此，结果可以不 是非常现实的。

可以代替放大尺度或除放大尺度之外使用的另一途径是多尺度仿 真。在多尺度仿真中，计算仍在较粗糙栅格上执行，但精细栅格信息 用来构造一组基础函数，该组基础函数可以用于映射精细栅格性质到 粗糙栅格。多尺度仿真可以比在精细栅格上的仿真快若干数量级，提 供具有可比较质量的解。一种多尺度方法是多尺度有限元法（MsFEM）， 见于例如T.Y.Hou和X.H.Wu，A multiscale finite element method for  elliptic problems in composite materials and porous media，J.Comput. Phys.，134：169-189，1997。该技术向其它数值多尺度法共享相似性， 例如变分多尺度有限元法，见于T.Arbogast和K.Boyd，Subgrid  upscaling an mixed multiscale finite elements，SIAM Num.Anal.，44： 1150-1171，2006；也见于T.J.R.Hughes，G.R.Feijoo，L.Mazzei和J.-B. Quincy，The variational multiscale method-A paradigm for computational  mechanics，Comput.Meth.Appl.Mech.Eng.，166：3-24，1998。该技 术也可以相似于由W.E和B.Engquist，The heterogeneous multi-scale  methods，Comm.Math.Sci.，1（1）（2003），pp.87--133提出的非均匀 多尺度法。已发展的另一途径是混合多尺度有限元法MMsFEM，其局 部质量守恒并可以应用到多相仿真。见于Z.Chen和T.Y.Hou，A mixed  multiscale finite element method for elliptic problems with oscillatory  coefficients，Math.Comp.，72：541-576，2002。如这些参考示出，对 多尺度方法的使用的深入研究已执行多于十年，尤其在学术界。尽管 具有将多尺度方法应用于更实际和复杂问题的趋势，但用多尺度方法 求解的仿真问题已简化并且是理论多于实践的。

上面讨论的多尺度方法使用局部信息并为可分离尺度良好执行。 然而，仅使用局部信息的多尺度方法可以受到谐振误差影响，该谐振 误差可以是多尺度仿真中的支配误差。谐振误差作为跨仿真单元的函 数的振荡示出，增加接近仿真网格中单元边缘的量值。谐振误差通常 与特性长度尺度和粗糙网格大小之间的比率成比例。在特性长度尺度 的量值与粗糙网格大小的量值明显不同时，该比率为小。在特性长度 尺度接近粗糙网格大小时，该比率变大。因此，误差为特性长度尺度 的整个范围示出，并且误差对于不同的特性长度尺度不同。见于例如 Y.Efendiev，V.Gunting，T.Y.Hou和R.Ewing，Accurate multiscale finite  element methods for two-phase flow simulations，J.Comp.Phys.，220（1）： 155-174，2006。使用一些受限全局信息可以对基础函数的构造有用， 以便发展减小或去除谐振误差的多尺度方法，见于J.E.Aarnes，Y. Efendiev和L.Jiang，Mixed multiscale finite element methods using  limited global information，SIAM MMS，7（2）：655-676，2008。此外， 使用全局信息的MsFEM可以应用于没有尺度分离的问题。

工业界和学术界研究者已为相似技术，多尺度有限体积法（MSFV） 报告结果。见于例如P.Jenny，S.H.Lee和H.A Tchelepi，Multi-scale  finite-volume method for elliptic problems in subsurface flow simulation， J.Comput.Phys.，187（2003），pp.47-67；也见于P.Jenny，S.H.Lee 和H.A Tchelepi，Adaptive multiscale finite-volume method for multiphase  flow and transport，SIAM MMS，3（2004），pp.50-64。

例如，Jenny等人的美国专利No.7,496,488披露用于地层流动仿真 的多尺度有限体积法。在该方法中，多尺度有限体积（MSFV）法用来 求解起因于多孔介质中的单或多相流动的带有多个空间尺度的椭圆问 题。该方法在粗糙网格上有效捕捉小尺度的效应，是守恒的，并正确 处理张量导磁率。基础思想是构造捕捉微分算子的局部性质的传递率。 这为有限体积解算法导致多点离散化方案。MSFV法的传递率优选仅 作为预处理步骤构造一次并可以局部计算。

Lee等人的美国专利申请公开No.2008/0208539披露使用包括黑油 建模的多尺度有限体积法的用于储层仿真的方法、装置和系统。多尺 度有限体积（MSFV）法在重力和毛细管力存在的情况下将非线性的不 可混溶三相可压缩流动建模。与MSFV构架一致，使用全隐式按序算 法分离地且不同地处理流动和传送。使用算子分裂算法求解压力场。 压力的通解分解成椭圆部分、浮力/毛细管力支配部分，以及带有源/ 汇和蓄积的不均匀部分。MSFV法用来计算椭圆分量的基础函数，捕 捉压力场中的长范围相互作用。在最初粗糙栅格上速率场的直接构造 和传送问题的解在调节物理机制中提供灵活性。MSFV法计算合适压 力场，包括过程尺度压力方程的解；构造精细尺度流量；以及计算相 传送方程。

如在上面参考中描述，MSFV法已用于结构化笛卡尔栅格，并紧 密涉及在石油工业中受欢迎的多点流量近似方案（MPFA）。MSFV依 靠都用于粗糙尺度压力的两组基础函数。

然而，具有使用MSFV的有限体积离散化的现有途径的若干问题。 例如，为基础函数求解问题可以是计算昂贵的。进一步地，MSFV没 有延伸到非构造化栅格，并且全局信息不用于MSFV。可以需要几何 信息从而施加边界条件，这对于非结构化栅格非常难以实施。另外， 多尺度基础仅使用局部信息，并且不可以求解全局特征如沟渠、不可 渗透的页岩障碍、裂缝等。

发明内容

本技术的示例性实施例提供用于将碳氢化合物储层建模的方法。 该方法包括使用多级混合多尺度有限体积（MMMFV）法从精细非结 构化网格导出可计算网格。导出可计算网格包括计算用于压力的第一 代数多级基础函数、构造互相作用区、生成主网格、计算用于速率近 似的第二代数多尺度基础函数，以及使用可计算网格仿真碳氢化合物 储层。该方法也包括至少部分基于仿真结果在非瞬时计算机可读介质 中生成物理碳氢化合物储层的数据表示。

该方法也可以包括通过选择辅助网格、对辅助网格上的解求近似， 以及使用提供的全局信息计算一个或多于一个速率基础函数，在第二 代数多尺度基础函数的构造中使用全局信息。进一步地，在该方法中 第一代数多级基础函数、第二代数多尺度基础函数或该两者可以至少 部分基于离散调和函数。

在该方法中，导出可计算网格可以包括计算主网格和一个或多于 一个互相作用区，并选择近似空间，其中该近似空间包括P和U。导 出可计算网格也可以包括选择测试空间，其中该测试空间包括Q和V， 并且选择离散梯度和散度算子。

在该方法中，计算第一代数多级基础函数可以包括在主网格中粗 糙可计算单元的中心点将多个边界条件设置成一、在主网格中粗糙可 计算单元的边界将该多个边界条件设置成零，以及使用该多个边界条 件求解基础函数。

在该方法中，计算粗糙基础可以包括在可计算网格∑_h中用直线[a_i， b]连接两点a_i和b。令最接近直线[a_i，b]的∑_h中的点为a_i+1，连接生成 从a_i到a_i+1的直线，并且如果a_i+1≠b，那么设置i＝i+1。计算粗糙基 础可以包括选择另一连接并重复，例如直到全部连接完成。

在该方法中，计算粗糙基础可以包括选择互相作用区的中心点， 选择互相作用区的一个或多于一个边缘的中心点，将互相作用区的中 心点连接到边缘的中心点，以及连接一个或多于一个边缘从而形成多 边形的面。

在该方法中，构造互相作用区可以包括通过构造三角测量的每个 边缘来构造粗糙点的三角测量。构造三角测量的每个边缘可以包括在 可计算网格∑_h中用直线[a_i，b]连接两个粗糙点a_i和b，并且令∑_h中的 点为a_i+1，其中该点最接近直线[a_i，b]。进一步地，构造三角测量的每 个边缘可以包括连接生成从a_i到a_i+1的直线；并且如果a_i+1≠b，那么 设置i＝i+1，并重复直到a_i+1=b。构造三角测量的每个边缘也可以包 括选择另一连接并重复。最终，构造三角测量的每个边缘可以包括通 过选择互相作用区的中心点，选择互相作用区的一个或多于一个边缘 的中心点，将互相作用区的中心点连接到边缘的中心点，以及连接多 边形的一个或多于一个面，来构造三角测量的一个或多于一个面。

该方法可以包括至少部分基于数据表示管理源自物理碳氢化合物 储层的碳氢化合物的生产。管理生产可以包括将注入井转变成生产井、 将生产井转变成注入井、改变生产率、钻探新井到储层，或其任何组 合。在该方法中，速率近似的质量可以通过计算特殊向量基础函数， 并为粗糙单元求解小型局部系统来改善。

本技术的另一示例性实施例提供用于从碳氢化合物储层生产碳氢 化合物的方法。该方法包括在非结构化网格上使用多级混合多尺度有 限体积法仿真碳氢化合物储层。仿真储层可以包括为压力计算第一代 数多级基础函数、构造互相作用区、生成主网格、为速率近似计算第 二代数多尺度基础函数，以及至少部分基于仿真结果从碳氢化合物储 层生产碳氢化合物。

在该方法中，生产碳氢化合物可以包括钻探一个或多于一个井到 碳氢化合物储层，其中该井包括生产井、注入井或该两者。进一步地， 生产碳氢化合物可以包括设置源自碳氢化合物储层的生产率。

本技术的另一示例性实施例提供用于仿真多孔介质中的流动的混 合多尺度有限体积法。该混合多尺度有限体积法包括收集关于碳氢化 合物储层的全局信息、基于饱和度方程导出椭圆问题、计算多个非结 构化粗糙网格和对应于该多个非结构化粗糙网格的多个压力基础函 数、求解该多个基础函数从而为多个非结构化粗糙网格的每个重构速 率向量场，以及迭代仿真直到达到最终时间帧。该仿真可以包括为多 个非结构化粗糙网格的每个中的多个计算单元求解压力方程、为多个 粗糙网格的每个计算总速率、为多个非结构化粗糙网格的每个求解椭 圆问题、为多个非结构化粗糙网格的每个计算相位速率，以及为多个 非结构化粗糙网格的每个计算成分运输。

混合多尺度有限体积法可以包括将为多个非结构化粗糙网格的每 个计算的相位速率组合，从而为碳氢化合物储层获得相位速率。混合 多尺度有限体积法也可以包括将为多个粗糙网格的每个计算的成分运 输组合，从而为碳氢化合物储层获得相位速率。

本技术的另一示例性实施例提供用于仿真碳氢化合物储层的系 统。该系统包括处理器和非瞬时存储设备，其中该存储设备包括碳氢 化合物储层的数据表示，其中该数据表示是由混合多尺度有限体积法 确定的非结构化可计算网格。该系统也包括内存设备，其中该内存设 备包括代码从而引导处理器为压力计算第一代数多级基础函数、构造 互相作用区、生成主网格、为速率近似计算第二代数多尺度基础函数、 使用主网格仿真碳氢化合物储层，以及至少部分基于仿真结果更新数 据表示。在该系统中，处理器可以包括多处理器集群。

附图说明

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

图1是根据本技术的示例性实施例的储层的示意图；

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

图3是根据本技术的示例性实施例的用于将储层建模的工作流程 的过程流程图；

图4是根据本技术的示例性实施例的用于为一个牛顿步的顺序公 式计算近似解的方法的过程流程图。

图5是称为Voronoi网格的可以用于示例性实施例的第一类网格的 图解；

图6是称为矩形网格的可以用于示例性实施例的第二类网格的图 解；

图7是称为四边形网格的可以用于示例性实施例的第三类网格的 图解；

图8根据本技术的示例性实施例示出非结构化的二维三角形网格；

图9是根据本技术的示例性实施例示出双线性基础函数的图表；

图10根据本技术的示例性实施例示出可以用来图解建立下级更粗 糙基础函数的进程的非结构化的二维三角形网格；

图11是根据本技术的示例性实施例的计算单元的图解，其图解用 于计算基础函数的技术；

图12是根据本技术的示例性实施例的从上面讨论的进程导致的多 尺度基础函数的例子；

图13是根据本技术的示例性实施例的代数多尺度函数的例子；

图14是根据本技术的示例性实施例的能量最小化基础函数的图 解；

图15是根据本技术的示例性实施例的精细可计算网格的图解，其 图解在过程中选择的粗糙点；

图16是根据本技术的示例性实施例的在精细可计算网格内粗糙点 的三角测量，该精细网格带有通过将粗糙点的每个与线段连接创造的 六个三角形；

图17是根据本技术的示例性实施例的矩形可计算网格的图解，其 示出三角测量的近似结果；

图18是根据本技术的示例性实施例的精细矩形可计算网格的图 解，其示出主网格的生成；

图19是根据本技术的示例性实施例的带有三角形互相作用区的可 计算网格的图解，其示出恒速场；

图20是根据本技术的示例性实施例的带有四边形互相作用区2000 的可计算网格；

图21是根据本技术的示例性实施例的示出不均匀基础函数的图 表；

图22是根据本技术的示例性实施例的互相作用区的图解，其示出 精细总速率的计算；

图23是根据本技术的示例性实施例的过程流程图，其示出用于使 用在此描述的技术执行仿真的方法；以及

图24是可以用于本技术的示例性实施例的示例性集群计算系统的 框图。

具体实施方式

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

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

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

“计算机可读介质”或“非瞬时计算机可读介质”如在此使用， 指代参与提供指令到处理器以便执行的任何非瞬时存储和/或传输介 质。这样的介质可以包括但不限于非易失性介质和易失性介质。非易 失性介质包括例如NVRAM或磁或光盘。易失性介质包括动态存储器 例如主存储器。计算机可读介质的共同形式包括例如软盘、软磁盘、 硬盘、硬盘阵列、磁带，或任何其它磁介质、磁光介质、CD-ROM、 任何其它光介质、RAM、PROM、EPROM、FLASH-EPROM、固态介 质如存储器卡、任何其它存储器芯片或盒式磁带，或计算机可以从其 读取数据或指令的任何其它有形介质。

如在此使用，“显示”或“正在显示”包括导致显示物理物体的 图形表示的直接行动，以及促进显示物理物体的图形表示的任何间接 行动。间接行动包括提供用户通过其能够影响显示器的网站、超链接 到这样的网站，或与执行这样的直接或间接行动的实体协作或合作。 因此，第一方可单独操作或与第三方卖主协作，从而使信息能够在显 示器设备上生成。显示器设备可以包括适合显示参考图像的任何设备， 无限制例如虚拟现实显示器、3d显示器、CRT监视器、LCD监视器、 等离子设备、平板设备或打印机。显示器设备可以包括已通过使用任 何常规软件校准的设备，该常规软件意图用于评估、校正和/或改善显 示效果（例如已使用监视器校准软件调整的彩色监视器）。代替（或 除此之外）在显示器设备上显示参考图像，与本发明相容的方法可以 包括向对象提供参考图像。“提供参考图像”可以包括创造或由实体 的、电话的或电子的输送分配参考图像、提供经由网络访问参考，或 创造或向对象分配软件，该软件经配置在包括参考图像的对象的工作 站或计算机上运行。在一个例子中，提供参考图像可以包括使对象能 够以硬拷贝的形式经打印机获得参考图像。例如，信息、软件和/或指 令可以传输（例如经数据存储设备或硬拷贝电子或物理传输）和/或以 其它方式可用（例如经网络），以便促进对象使用打印机打印参考图 像的硬拷贝形式。在这样的例子中，打印机可以是已通过使用任何常 规软件校准的打印机，该常规软件意图用于评估、校正和/或改善打印 结果（例如已使用色彩校正软件调整的彩色打印机）。

“示例性”在此专门用来意谓“用作例子、实例或说明”。在此 描述为“示例性”的任何实施例不解释为超过其它实施例优选或有利。

“流动仿真”定义为使用仿真模型将传送通过物理系统的质量（通 常，流体例如油、水和气体）或能量仿真的数值方法。物理系统可以 包括三维储层模型、流体性质以及井的数目和位置。流动仿真可以使 用或提供策略（经常称为井管理策略）以便控制注入和生产率。这些 策略可以用来通过用注入流体（例如水和/或气体）代替生产流体来维 持储层压力。在流动仿真正确再创造过去的储层性能时，其称为“历 史匹配的”，并且更高程度的置信度置于其能力中从而预测储层中的 未来流体行为。

“渗透度”是岩石传递流体通过岩石的互连孔隙空间的能力。可 以使用达西定律测量渗透度：Q=（k ΔP A）/（μL），其中Q=流速 （cm³/s），ΔP=跨长度L（cm）和截面积A（cm²）的圆柱体的压降 （atm），μ=流体粘度（cp），并且k=渗透度（达西）。渗透度测 量值的习惯单位是毫达西。术语“相对可渗透”关于地层或其部分定 义为10毫达西或更多（例如100毫达西）的平均渗透度。

“孔隙度”定义为以百分比表达的孔隙空间的体积对材料总毛体 积的比率。孔隙度是储层岩石的流体存储容量的测量。孔隙度优选从 岩心、声测井记录、密度测井记录、中子测井记录或电阻率测井记录 获得。总或绝对孔隙度包括全部孔隙空间，而有效孔隙度仅包括互连 孔隙并对应可用于消耗的孔隙体积。

“储层”或“储层地层”定义为包括砂岩、石灰岩、白垩、煤和 一些类型页岩的产油层（例如，碳氢化合物生产层）。产油层可以在 厚度上从小于一英尺（0.3048m）到数百英尺（数百m）变化。储层地 层的渗透度为生产提供潜力。

“储层性质”和“储层性质值”定义为表示含有储层流体的岩石 的物理属性的量。术语“储层性质”如在本申请中使用，包括可测量 和描述性的属性。可测量储层性质值的例子包括孔隙度、渗透度、含 水饱和度和裂缝密度。描述性储层性质值的例子包括外观、岩石学（例 如砂岩或碳酸盐），以及沉积环境（EOD）。储层性质可以构成储层 构架从而生成储层模型。

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

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

“井”或“井孔”包括套管的、套管并胶结的或裸眼的井孔，并 可以是任何类型的井，包括但不限于生产井、试验井、探井等。井孔 可以是垂直的、水平的、在垂直和水平之间任何角度的、偏离的或不 偏离的，及其组合，例如带有不垂直组成的垂直井。井孔通常钻探并 然后通过在井孔内安置套管柱来完成。常规地，通过使水泥流通进入 在套管柱的外表面和井孔面之间定义的环带，套管柱胶结到井面。套 管柱一旦在井内水泥中嵌入，则穿孔从而允许跨感兴趣间隔的管道内 外之间的流体连通。该穿孔允许处理化学品（或物质）从套管柱里面 流入周围地层，以便刺激流体的生产或注入。后来，穿孔用来从地层 接收碳氢化合物的流动，因此它们可以通过套管柱输送到地面，或为 储层管理或处置目的允许流体连续注入。

概述

有限体积（FV）法是用于储层仿真的高效的离散化方法。FV法可 以为可接受计算成本提供合理的准确度，并可以延伸到非结构化栅格。 在本技术的示例性实施例中，FV法的多尺度一般化应用到非结构化栅 格。由于使用多尺度算法的离散化构造，因此该方法是计算有效的。

该方法使用多级途径建立粗糙栅格基础函数。离散调和函数可以 用于粗糙栅格基础，这可以减小或消除边界条件的影响，并使本技术 可应用于非结构化栅格。可以然后以提供最优近似，尤其例如能量最 小化的方式在每级上建立支集。本技术使用精细的和粗糙的非结构化 栅格，并在过程中包括全局信息。

图1是根据本技术的示例性实施例的储层102的示意图100。储层 102，例如石油或天然气储层，可以是可以通过从地面110钻探井104、 106和108通过覆盖层112进入的地下地层。井104、106和108可以 偏离，例如方向地钻探从而沿着储层102。进一步地，井可以分支从而 增加可以从储层抽取的碳氢化合物的量，如为井104和108示出。井 104、106和108可以具有带有穿孔120（指示为靠近井的点）的众多 区域从而允许碳氢化合物从储层102流入井104、106和108以便移除 到地面。储层102可以具有可以约束或增强碳氢化合物流动的一个或 多于一个断层114分割区，例如区域116和118。

储层102的仿真模型或仿真器可能发现在井104、106和108附近 发生的最大改变，以及其它储层特征，例如断层114。因此，在这些特 征的每个附近保持精细结构是有用的。

图2是根据本技术的示例性实施例的储层顶视图，其示出可计算 网格200在储层上的平面投影。尽管可计算网格200示作计算单元（或 块）202的二维栅格从而简化问题的解释，但应理解实际可计算网格 200可以是包括储层的计算单元202的三维结构。进一步地，计算单元 202可以具有任何大小或形状，导致非结构化栅格。一般地，计算单元 202是表示储层中物理位置的仿真模型内单独二维或三维位置。计算单 元202可以具有关联属性，例如孔隙度，该关联属性在整个计算单元 202上假设为单值，并赋值到计算单元202的中心。例如通过具有赋值 到与邻近计算单元202的共享边沿的流量性质，计算单元202可以与 邻近计算单元202相互作用。流量性质可以包括由在邻近计算单元202 之间的温度或压力差驱动的热或质量转移。

例如通过将不接近井或其它储层特征的计算单元202组合，可计 算网格200可以在具有较少显著改变的区域中粗化。相似地，如在图2 中示出，可计算网格200可以在井或其它储层特征，例如第一井204， 或其它储层特征例如第二井206、第三井208、断层210，或可以示出 大于其它区域的改变的任何其它特征附近保持精细网格结构。可计算 网格200可以用来将储层建模，如关于图3进一步讨论。

用于将储层建模的工作流程

图3是根据本技术的示例性实施例的用于将储层建模的工作流程 300的过程流程图。尽管求解过程的离散化（粗化）和隐性水平（其规 定变量，例如压力或饱和度在公式中隐性或显性处理）变化，但仿真 模型可以用与工作流程300相似的方式执行。仿真模型可以通过解析 用户输入数据在方框302开始。输入数据可以包括问题公式、带有在 每个栅格块定义的物理性质的离散化进入栅格块的地质模型，该物理 性质包括岩石性质（例如渗透度）和流体性质（例如传递率）。在方 框304，例如可以从基本方程计算用于仿真的边界条件。

在方框306，线性解算器可以使用雅可比矩阵为仿真生成近似解。 在方框308，从近似解计算物理性质。在方框310，将计算性质与先前 计算性质或与测量性质比较，从而例如通过检查收敛来确定是否达到 希望的准确度。在示例性实施例中，通过确定计算性质自从最近迭代 （其可以指示收敛）没有显著改变，因此做出确定。例如，收敛可以 指示当前计算解是否在先前计算解的0.01%、0.1%、1%、10%或更多 之内。如果没有达到希望的准确度，那么过程流程返回到方框306从 而执行线性解算器的另一迭代。

如果在方框310已达到希望的准确度，那么过程流程进展到方框 312，在方框312结果生成并且时间以希望的时间步进增加。该结果可 以存储在有形的、机器可读的介质例如数据库上的数据结构中，以便 后来呈现，或该结果可以在生成之后立即显示或打印。时间步进可以 是例如一天、一周、一个月、一年、5年、10年或更多，至少部分取 决于仿真的希望时间长度。在方框314，新时间与仿真希望的长度比较。 如果仿真已达到希望的时间长度，那么仿真在方框316结束。如果时 间没有达到希望的长度，那么流程返回到方框304，从而继续下个增加。 仿真时间帧可以是一个月、一年、五年、十年、二十年、五十年或一 个世纪或更多，取决于希望的仿真结果使用。

计算近似解

在方框306计算近似解可以包括在一系列步骤求解物理模型。例 如，在N_P相中的N_C成分的流动的顺序隐式公式遵循在图1中示出的 压力方程。

$\phi [\frac{1}{\phi }·\frac{\mathrm{d\phi }}{\mathrm{dp}}-\frac{1}{v_T}\left(\frac{\partial v_T}{\partial p}\right)]\frac{\partial p}{\partial t}=\underset{i=1}{\overset{N_c}{\Sigma }}{V}_{\mathrm{Ti}}{\Theta }_i$方程1

在方程1中，□表示孔隙度，并且V_Ti表示成分i的偏摩尔体积。进一 步地，其中U_i表示成分i的摩尔流率，并且R_i表示成 分i的源和汇。成分i的摩尔流率在方程2中定义。

$U_i=\underset{J=1}{\overset{N_P}{\Sigma }}{x}_{\mathrm{ij}}{\xi }_j{v}_j$方程2

在方程2中，x_ij是相j中分量i的摩尔分数，ξ_j是相j的摩尔密度，并 且v_j表示相位速率。广义达西定律给出相位速率和相压力之间的关系， 如在方程3中示出。

$v_j=-\frac{K{k}_{\mathrm{rj}}}{{\mu }_j}(▿P_j-{\rho }_{j}g▿D),$方程3

在方程3中，K表示绝对渗透度，k_rj表示相j的相对渗透度，μ_j表示相 j的粘度，ρ_j表示相j的质量密度，g表示重力加速率常数，并且D表 示深度。对于带有蒸汽相（v）和水相（a）的液相（l），相压力通过 毛细压力P_c，j，j=v，a，涉及底部液体压力：P_j=P+P_c，j，j=v，a， P=P_l。

饱和度定义为比率其中V_j是相j的体积，并且V_P是孔隙体积。 饱和度方程如在方程4中示出。

$\frac{\partial }{\partial t}\left(\phi S_j\right)+\phi S_j{c}_j\frac{\partial p}{\partial t}=\underset{i=1}{\overset{N_C}{\Sigma }}{V}_{\mathrm{ji}}{\Theta }_j,$j=1，2，...，N_P。方程4

在方程4中，V_ji表示关于分量i的相j的偏摩尔体积，并且c_j表示相j 的压缩率。这些方程可以用来使用在图4中示出的方法求解近似。

图4是根据本发明的示例性实施例的用于为一个牛顿步的顺序公 式计算近似解（即图3的方框306）的方法400的过程流程图。方法 400用压力方程的求解在方框402开始（方程1）。在方框404，使用 方程2、3计算在这些压力的流体流动的总速率。饱和度方程（方程4） 在方框406显性或隐性求解。在方框408，计算每个单元202（图2） 中相位速率、成分运输和每个分量的量。方法400可应用于复合与黑 油模型。

混合有限体积法

有限体积离散化可以为线性椭圆方程确定，即在类时间欧拉后向 法中的隐式离散化执行之后。考虑在方程5中示出的方程组。

$-\mathrm{div}(K▿p\left(x\right))=f\left(x\right)$在Ω中，

方程5

p(x)=0在上。

如果那么方程5可以重写为在图6中示出的方程组。

$K^{-1}u+▿p=0$方程6

div(u)=f(x)

在方程6中的第一方程可以乘以向量函数v，并且在方程6中的第二方 程可以与标量函数q相乘。方程的积分导致问题：找出u，u∈U和p， p∈P以满足方程7。

$\left(\begin{array}{cc}\underset{\Omega }{\int }{K}^{-1}u·v+\underset{\Omega }{\int }▿p·v=0& \forall v\in V,\end{array}\right)$方程7

$\underset{\Omega }{\int }\mathrm{div}\left(u\right)q=\underset{\Omega }{\int }f\left(x\right)q,\forall q\in Q.$

如果空间U，V和P，Q适当挑选，那么在方程7中存在的问题等 效于在方程5中存在的问题。例如，本领域技术人员认识到U是函数 的空间，该空间的散度是平方可积分函数。进一步地，该V是平方可 积分向量函数的空间；P是带有平方可积分第一导数的函数的空间；以 及Q是平方可积分函数的空间。可以通过求解近似离散问题找出对该 问题的近似解u_h和p_h。第一步是用有限维子空间U_h，V_h和P_h，Q_h替 代有限维空间U，V和P，Q。接下来，连续算子（梯度）和div（散 度）与离散算子和div_h近似。因此，离散问题是找出u_h∈U和p_h∈ P_h以满足方程8。

$\left(\begin{array}{cc}\underset{\Omega }{\int }{K}^{-1}{u}_h·v_h+\underset{\Omega }{\int }{▿}_h{p}_h·v_h=0& \forall v_h\in U_h\end{array}\right),$方程8

$\underset{\Omega }{\int }{\mathrm{div}}_h\left(u_h\right)q_h=\underset{\Omega }{\int }f\left(x\right)q_h,\forall q_h\in Q_h.$

混合有限体积栅格的例子

在下面提出可计算网格的若干例子。在这些例子的全部中，具有 主计算单元和双计算单元（其可以称为“互相作用区”）。如果描述 主网格和对偶网格、定义近似（测试和试验）空间，并且指定离散算 子，那么完全定义如由方程8描述的特别的混合有限体积法。

图5是称为Voronoi网格500的可以用于示例性实施例的第一类网 格的图解。在Voronoi网格500中，主计算单元是Voronoi体积502， 并且互相作用区是三角形504。通过将每个Voronoi体积502的中心点 506连接到与Voronoi体积502共享面510的邻接Voronoi体积508， 来形成三角形504。认为在中心点506之间的直线512垂直于面510， 并因此平行于法向向量n。d_ij是连接节点I和j的间隔的长度。认为适 当参数例如c、p和s在中心点506存在。

P_h是在Voronoi体积双倍的Delaunay三角形上的线性分段连续函 数的空间。这些是三角形504。Q_h是Voronoi体积502上分段常数的空 间。进一步地，U_h=V_h是带有连续法向分量的三角形上分段向量常数 的空间。可以与Voronoi网格500一起使用的算子包括以及 在一些实施例中，Voronoi网格500可以由 Donald网格替代。在Donald网格中，可以通过将三角形502重心连接 到三角形502边缘的中心来形成控制体积502。

图6是称为矩形网格600的可以用于示例性实施例的第二类网格 的图解。矩形网格600具有主计算单元602和互相作用区604。通过将 计算单元602的中心点608连接的直线606创造互相作用区604。P_h= Q_h是单元中心网格，即单元602上分段常数函数的空间。进一步地， U_h=V_h是带有连续法向分量的双三角形上分段向量常数的空间。

可以与矩形网格600一起使用的算子包括其中[]表明在计算单元602之间边缘610的不同侧面上的值之间的差 （跳跃），例如另一可能MFV法可以经构 造假设P_h是分段双线性函数的空间。在此实例中，对于每个 单元中心608具有一个双线性基础函数。基础函数在该特别点是一， 并在单元外是零。

图7是称为四边形网格700的可以用于示例性实施例的第三类网 格的图解。在四边形网格700中，主计算单元是四边形体积702，并且 互相作用区包括通过连接四边形702的中心点706形成的四个多边形 704。P_h是在每个四边形702内创造的三角形708上分段线性函数的空 间。每个三角形708的一个顶点在单元中心点706中，而其它两个顶 点710在四边形702的边缘712上。进一步地，Q_h是四边形702上分 段常数的空间，U_h是带有连续法向分量的多边形704上分段向量常数 的空间，并且V_h是多边形704上分段常数向量的空间。因此，U_h和 V_h可以由从三角形708的顶点710绘制的向量714表示。

可以与四边形网格700一起使用的算子包括和另一版本具有由四边形702和多边形704 的互相作用区形成的四边形子网格上的分段双线性函数的空间P_h（例 如，互相作用区的四分之一）。可以然后为多边形704的中心添加另 外方程。在此实例中

上面例子的全部可以延伸到三维网格。特别例子是三维Voronoi 栅格，所谓2.5-D PEBI（中垂线）栅格，该栅格在xy平面中是Voronoi 的并且在z方向上是棱柱，并由一致或变形的平行六面体构成。

多级基础函数

在示例性实施例中，使用多级技术建立空间P_h。因此，下面讨论 可应用于有限体积和有限元离散化。在多尺度基础的多级构造中，如 在此讨论，用于仿真的计算花费和精细网格上未知数的数目成比例。 相反，不使用多级技术构造的标准多尺度基础具有与至少精细网格上 未知数数目的平方成比例的计算花费。为简化该概念的解释，可以在 以下附图的解释中假设具有精细非结构化网格，并且P_h的基础在该网 格上定义。

图8根据本技术的示例性实施例示出非结构化的二维三角形网格 800。基础函数可以为有限元离散化在每个顶点802定义。基础函数的 值在顶点802是一，并且在其支集区804之外是零。三个基础函数的 支集区804在图8中示出。函数是分段线性的，即在每个三角形上是 线性的。有限元法的基础函数在每个单元上是常数。

图9是根据本技术的示例性实施例示出双线性基础函数的图表 900。对于矩形网格600（图6），函数是双线性的。在该图表900中， 轴902表示网格中的物理位置。基础函数的峰值904在计算单元的中 心。

在多级进程中建立下级更粗糙基础

图10根据本技术的示例性实施例示出可以用来图解建立下级更粗 糙基础函数的进程的非结构化的二维三角形网格1000。一个或多于一 个粗糙点1002可以在网格上选择。粗糙点是有限元或有限体积的单元 中心的顶点，该顶点是更粗糙基础函数的“中心”，即其中基础函数 是一的点。本领域技术人员认识到具有可以用于选择粗糙点和粗糙网 格的众多算法。通常，粗糙网格和在其上定义的离散化必须能够求解 储层中的流动。因此，可以选择粗糙点从而表示可以具有储层中最高 级改变例如井等等的点。

一旦选择粗糙点，那么与每个粗糙点1002关联的精细点1004可 以在网格1000上鉴别。可以注意在精细网格上的一些函数可以属于若 干粗糙网格点的支集。一旦已鉴别精细点，那么可以在每级计算基础 函数。进程可以递归重复直到最粗糙基础函数具有足够小的尺寸。

粗化的目标是可观减小粗糙网格的大小，因此计算花费也减少。 另一方面，产生的粗糙基础必须足够精细从而准确近似问题的解。为 平衡这两个对立目标，可以使用一些进一步的准则。由于基础函数仅 计算一次或非常不频繁地计算，因此基础的质量可以检查。这可以用 不同方式完成，例如通过将模型椭圆问题的精细解与粗糙解比较。这 样的进程允许适应地粗化网格，同时提供良好近似性质并与原精细网 格的大小比较减小最终粗糙网格的大小。基础的质量可以取决于选择 粗糙网格点和怎样构造基础函数的方式。

图11是根据本技术的示例性实施例的计算单元1100的图解，其 图解用于计算基础函数的技术。具有计算基础函数的若干不同方式。 粗糙基础函数的支集包括数个三角形1102，其中中心点1104构成每个 三角形1102的顶点中的一个。方程5可以在每个三角形1102上但用 指定的不同边界条件求解。基础函数的值可以在中心点1104设置成一 并在每个边界1106设置成零。一维问题可以然后沿每个边缘1108求 解。一维问题的解可以然后用作二维基础函数的边界条件。用于三维 基础函数的进程首先为边缘求解一维问题，并然后使用该解作为面的 二维问题的边界条件。二维问题的解提供基础函数的边界条件。

图12是根据本技术的示例性实施例的从上面讨论的进程导致的多 尺度基础函数的例子。然而，上面概述的计算进程为网格使用几何信 息，并因此可以是麻烦的并且计算密集的。

在示例性实施例中，可以考虑可替换途径。使用该途径计算的基 础函数可以称为代数多尺度。关于先前的进程，边界条件可以在支集 的中心点1104中设置成一，并且在边界1106上设置成零。方程5可 以然后用上面指定的边界条件求解，沿每个边缘消除个别一维问题的 解。

图13是根据本技术的示例性实施例的代数多尺度函数的例子。可 以注意代数多尺度函数的计算比关于图12讨论的解容易。给定全局精 细矩阵A，行和列可以为粗糙点的支集中的点选择，从而形成矩阵A_c， 如在图9中示出。

$b=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right),$方程9

或等效地（c=1）：

$\left(\begin{array}{cc}{a}_{f_1{f}_1}& a_{f_1{f}_2}\\ a_{f_2{f}_1}& a_{f_2{f}_2}\end{array}\right)\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right)=-\left(\begin{array}{c}{a}_{f_{1}c}\\ a_{f_{2}c}\end{array}\right)$

方程9表示在为特别粗糙点将右手侧设置成一并为全部其它点设置成 零之后有待求解的线性问题。为每个粗糙点求解一个线性问题。

可以注意几何信息不用来计算这些基础函数。因此，代数多尺度 基础可以对计算非结构化网格上的基础函数有用。一旦已求解局部问 题，那么它们可以改变尺度从而确保基础函数和为一。这可以使用方 程10中的公式执行。

$f_c^i\left(x\right)=\frac{f_c^i\left(x\right)}{\underset{j}{\Sigma }{f}_c^j\left(x\right)}.$方程10

其它途径可以用来在当前技术的示例性实施例中计算基础函数。 例如，全部粗糙函数可以这样同时计算，以便最小化基础函数的能量。 该进程可以称为能量最小化基础。能量最小化基础函数可以数学表达， 如在方程11中示出。

$\min \frac{1}{2}\underset{i=1}{\overset{m}{\Sigma }}{f}_c^i{\left(x\right)}^T{A}_c^i{f}_c^i\left(x\right),$以使方程11

图14是根据本技术的示例性实施例的能量最小化基础函数的图解。上 面描述的技术可以一起引入从而实施混合多尺度有限体积法。

混合多尺度有限体积（MsMFV）法

对网格段求近似

在实施MsMFV法中的第一步是对在给定网格Σ_h和对偶网格Ψ_h上的线段[a，b]求近似。如下面描述，a，b∈Σ_h或a，b∈Ψ_h。可 以使用两种算法。第一算法是在二维网格中线段的近似。第二方法对 平面多边形，例如在三维网格中的平面多边形求近似。

为简化用于线段近似的第一算法的解释，可以假设a，b∈Σ_h， 并且起始点是a。该算法在此称为“第一算法”。首先，点a和b可以 与直线连接，提供设置I=0和a=a₀。如果还没有挑选，那么选择来 自a_i和最接近[a，b]的Σ_h的点的连接。该点可以设置成a_i+1。将a_i连接 到a_i+1。如果a_i+1≠b那么设置i＝i+1。为点a_i+1重复相同步骤。如果 已达到点b，那么线段的近似完成。尽管第一算法对构造二维网格有用， 但需要另外的工具以便构造三维网格。

第二算法为三维网格提供多边形表面的近似。该算法在此称为“第 二算法”。给定多边形表面及由第一算法构造的其边缘的近似，可以 选择多边形表面的一个顶点。在该顶点之外具有至少两条线段。选择 属于精细网格的面（多边形表面）的两条线段。如果在该多边形表面 中具有之前没有近似并在其之外具有至少两个线段的点，那么重复相 同进程。如果源自全部点的线段都已近似，那么近似完成。

一旦在网格上的多边形已近似，那么MsMFVM可以构造。起始， 可以生成粗糙基础函数，如上面关于图10-13描述。

图15是根据本技术的示例性实施例的精细可计算网格1500的图 解，其图解在过程中选择的粗糙点1502。可以然后构造对偶网格，即 互相作用区。这可以通过建立粗糙点1502的三角测量来执行。该三角 测量包括在二维网格中的三角形和在三维网格中的四面体或棱柱。

图16是根据本技术的示例性实施例的在精细可计算网格1500内 粗糙点1502的三角测量，该精细网格1500带有通过将粗糙点1502的 每个与线段1604连接创造的六个三角形1602。可以然后使用第一算法 对三角测量的每个边缘（线段1604）求近似。该点是精细网格单元的 中心，并且我们在互相作用区（对偶网格）的边缘上移动。即，精细 网格点的一些（单元的中心）可以选择为粗糙网格点。接下来，基础 函数建立。最终，通过首先建立边缘（使用第一算法），并在此之后 建立面（使用第二算法），来建立对偶网格。

图17是根据本技术的示例性实施例的矩形可计算网格1500的图 解，其示出三角测量的近似结果。是每个原三角形1602（图16）的边 缘的线段1604由线段1702近似，创造对应于三角形1602的六个多边 形1704。在精细网格1500非结构化时，由线段1604近似边缘（线段 1604）变得更准确。由近似边缘1702跨越的面可以使用第二算法近似。

生成主网格

图18是根据本技术的示例性实施例的精细矩形可计算网格1500 的图解，其示出主网格（由粗线1802指示）的生成。粗糙网格点1804 用来使用上面讨论的进程构造粗糙互相作用区1806（由虚线1808描 绘）。通过首先在每个粗糙互相作用区1806里面选择点1810来创造 主网格。接下来，主精细网格到选择点1810的最接近顶点可以鉴别。 最接近点鉴别为互相作用区1806的中心点1812。如果在精细网格（例 如1814）中具有处在到选择点1816的相同距离的多于一个点，那么第 一点可以挑选为中心点1812。在每个互相作用区1806中鉴别中心点 1812之后，最接近粗糙互相作用区1806的每个边缘（虚线1808）的 中心的精细网格上的点选为该边缘的“中心”。第一算法可以然后用 来将互相作用区的中心连接到边缘的中心。中心点1812变成精细主网 格（由粗线1802指示）中的顶点。

接下来，考虑主网格（由粗线1802指示）的边缘。如果问题是二 维的，那么这是粗糙主网格。对于三维问题，这些是粗糙主网格的面 的边缘。第二算法可以用来从边缘构造面。

构造速率场

对于为速率场计算基础函数具有两个可能情况。第一情况是在互 相作用区相似于三角形时，如在图19中示出。第二情况是在互相作用 区是四边形时，如在图20中示出。通常，四边形情况可以分为两个三 角形，但经常更方便的是保持四边形。进一步地，具有构造速率场的 多个途径。例如，可以使用速率的常数近似。由于速率可以变化，因 此多尺度速率近似也可以使用。

图19是根据本技术的示例性实施例的带有三角形互相作用区1900 的可计算网格的图解，其示出恒速场。粗糙中心点1901选择并用来生 成粗糙互相作用区1900，如关于图15-17讨论。恒速场是速率的最简 单近似，并且其使用上面提出的技术构造。粗线1902是粗糙单元的边 界，细线1904是精细单元的边界，并且虚线1906是互相作用区的边 界。速率场（由箭头1908指示）在互相作用区1900的每个部分1910 中是常数，并且具有连续法向分量。计算的速率是跨互相作用区1900 的平均速率。在许多情况下提供充足准确度，同时避免源自更复杂计 算的开销。然而，常数近似不考虑系数的变化，并且在一些状况下采 取更全面途径可以是有用的。

图20是根据本技术的示例性实施例的带有四边形互相作用区2000 的可计算网格。使用上面讨论的技术，更粗糙中心点2002选择并用来 生成互相作用区2000。在计算单元之间的边缘2004由粗线指示。流动 的速率2006由箭头指示。

速率的更优近似可以通过使用更优基础函数构造。为每个互作用 子区使用两个基础函数。例如，对于带有边界e₁、e₂、l₁₁、l₁₂的互作用 子区K₁₄，使用混合有限体积或混合有限元离散化解决在方程8中说明 的问题，如在方程12中示出。

在K₁₄中，-div(v₁)=0，  在K₁₄中，-div(v₁)=0，

在e₄上，v₁·n=1，      在e₄上，v₂·n=0，

在l₁₂上，v₁·n=-1，    在l₁₂上，v₂·n=0，      方程12

在e₁上，v₁·n=0，      在e₁上，v₂·n=1，

在l₁₁上，v₁·n=0，     在l₁₁上，v₂·n=-1

可替换地，可以通过在两个邻近子区中求解方程12构造基础函数。注 意如果在每个子区K₁₄中的每个精细单元中系数相同，那么基础函数 v₁和v₂是常数，并与上面讨论的混合有限体积栅格的例子中的基础一 致。因此，如在图20中示出，局部速率不再相同。因此，速率基础中 的每个函数构造成分段常数，并具有连续的法向分量。

在上面提出的本技术的示例性实施例中，已构造主网格和对偶网 格与近似空间P_h和U_h。测试空间Q_h挑选为主网格的体积上分段常数 的空间。空间V_h是带有连续法向分量的分段常数向量函数的空间。可 以在实施例中使用其它空间，例如P_h=Q_h和U_h=V_h。这些空间可以具 有相同的近似性质，但需要稍多的计算。

全局信息的使用

关于储层的全局信息可以用来改善粗糙网格解的质量。例如，全 局信息可以用于速率的基础函数的计算。如果通过求解简化问题或从 地质信息获知流动的行为，那么向量场w可以在比粗糙网格更精细的 网格上确定。因此，可以选择中间辅助网格，即比粗糙网格更精细， 但比精细网格更粗糙的辅助网格。粗糙基础函数可以使用在方程13中 示出的公式在边缘e上计算。

$v·n_{|e_i}=\frac{w·n_{|e_i}}{\underset{e}{\int }w·n},$i＝1，2，...，k，          方程13

其中w已在中间网格上计算。

在方程13中，边缘分为k个部分。这改变基础函数的形状。

图21是根据本技术的示例性实施例的示出不均匀基础函数的图 表。在该例子中，左流量（a）是均匀的。然而，通过比较，右流量（b） 是均匀的，在第二半部中具有多出三倍的流动。一般地，这可以对应 于沿一侧具有更可渗透岩层的地层。在仿真期间基础函数的形状改变 流动的行为。

注意并因此可以使用在方程14中示出的修改问题来计 算方程8中的基础函数。

在K₁₄中，-div(v₁)=0，在K₁₄中，-div(v₂)=0，

在e₄上，$\underset{e_4}{\int }{v}_1·n=1,$在e₄上，v₂·n=0，

在l₁₂上，$\underset{l_{12}}{\int }{v}_1·n=-1,$在l₁₂上，v₂·n=0，方程14

在e₁上，v₁·n=0，在e₁上，$\underset{e_1}{\int }{v}_2·n=1,$

在l₁₁上，v₁·n=0，在l₁₁上，$\underset{l_{11}}{\int }{v}_2·n=-1$

一致速率的计算

图22是根据本技术的示例性实施例的互相作用区的图解，其示出 精细总速率的计算。无关于使用恒速或变速，第一步是使用如在此讨 论由混合多尺度有限体积法产生的近似u计算总速率。这对一些情况 可以是足够的，但在希望精细网格上的更准确总速率时，结果可以后 处理从而得到更高准确度。第一途径是求解方程5，但用Neumann边 界条件，即如在方程15中示出。

$-\mathrm{div}(K▿p\left(x\right))=f\left(x\right)$在K中，

方程15

$-K▿p\left(x\right)·n=u_{H\{\partial K}·n$在上。

方程15在精细网格上使用混合有限体积或混合有限元法。该途径可以 是计算昂贵的。另一方式是使用混合有限体积构架求解粗糙问题。首 先可以使用在方程16中的公式计算基础函数v₁₂、v₁₃、v₂₄、v₃₄。

在K₁∪K₂中，-div(v₁₂)=0，在K₁∪K₃中，-div(v₁₃)=0，

在l₁₂上，v₁₂·n=1，在l₁₃上，v₁₃·n=1，

在l₁₃∪l₂₄上，v₁₂·n=0，在l₁₂∪l₃₄上，v₁₃·n=0，

在e₁₁∪e₁₂上，v₁₂·n=0，在e₁₁∪e₁₂上，v₁₃·n=0，

在e₂₁∪e₂₂上，v₁₂·n=0，在e₃₁∪e₃₂上，v₁₃·n=0，

                                                方程16

在K₂∪K₄中，-div(v₂₄)=0，在K₃∪K₄中，-div(v₃₄)=0，

在l₂₄上，v₂₄·n=1，在l₃₄上，v₃₄·n=1，

在l₁₂∪l₃₄上，v₂₄·n=0，在l₁₃∪l₂₄上，v₃₄·n=0，

在e₂₁∪e₂₂上，v₂₄·n=0，在e₃₁∪e₃₂上，v₃₄·n=0，

在e₄₁∪e₄₂上，v₂₄·n=0，在e₄₁∪e₄₂上，v₃₄·n=0

相似地，计算在外侧边缘上速率的基础函数。在方程16中的计算 仅完成一次。它们唯一定义空间U。空间V可以与空间U相同，或仅 具有分段常数。空间P是分段常数或双线性函数的空间，并且空间Q 由分段常数构成。下个步骤是在方程15的离散化版本中求解u和p。 这通过为每个粗糙单元求解小方程组来执行。计算的精细网格流量在 全部精细网格中一致。进程不限于速率，但可以用于改善任何向量变 量的精细网格近似和一致性。这可以是总速率、相位速率或感兴趣的 一些其它向量变量。

用于多孔介质中流动的仿真的多级混合多尺度有限体积法的应用。

图23是根据本技术的示例性实施例的过程流程图，其示出用于使 用在此描述的技术执行仿真的方法2300。在该示例性实施例中，一般 目标是用由方程1和方程4描述的模型执行储层仿真。这是在图4中 示出的过程流程图的延伸。该方法在方框2302用问题的分析开始。该 分析包括为跨储层的流动的行为收集全局信息。在这点上，可以关于 使用全局信息示例性有益做出决定。全局信息可以从简化仿真搜集， 或可以由专家提供。

在方框2304，导出基于方程4的椭圆问题。该椭圆问题是方程4 的简化，通过跳过时间导数来创造。多级网格在方框2306计算，如在 此讨论。在方框2308，为多级网格计算基础。在方框2310，计算与压 力关联的均匀速率场。如果在方框2312做出希望更准确速率场（例如， 带有变化压力）的决定，那么在方框2314计算更准确的基础函数。在 任一情况下，流程进展到方框2316。在方框2316，评估压力基础函数 的质量。这可以通过确定近似怎样良好来执行，例如通过对源自储层 的流动数据的历史匹配来执行。如果基础函数中的任何不与希望一样 良好，那么可以重计算它们。流程然后进展到方框2318。

方框2318标志迭代仿真过程的开始。在方框2318，使用在此描述 的多级混合多尺度有限体积法将在方程1中示出的压力方程离散化， 并然后求解。系数的一些可以在每个时间步进用源自先前时间步进的 数据重计算。在方框2320，计算总速率。为执行该函数，使用源自方 框2310的解计算粗糙网格总速率流量。然后，如希望，那么使用在方 框2130或在方框2314计算的速率基础函数计算精细网格总速率。

在方框2322，求解在方程4中示出的饱和度方程。饱和度方程可 以在粗糙网格上隐性求解或在精细网格上显性求解。适应性解也是可 能的。例如，在饱和度的隐性粗糙解之后，其中饱和度改变高于给定 容限的粗糙单元可以选择，并且可以使用源自粗糙网格解的边界条件 求解精细网格问题。可以执行若干精细时间步进，直到达到粗糙时间 步进。精细时间步进可以选择，因此精细网格显性饱和度问题是稳定 的，并且多个精细网格时间步进等于粗糙网格时间步进。

在方框2324可以计算相位速率，并且在方框2326计算成分运输。 这些计算可以用与在方框2324的总速率计算相似的方式执行。可以使 用如在方框2322描述的适应性算法。不同的相和分量可以使用不同的 适应性区域。

在方框2328，时间以选择的时间步进增加。该时间步进可以取决 于仿真的总时间。例如，总仿真可以覆盖1年、5年、20年、40年或 更多。该时间步进可以覆盖1个月、2个月、6个月、1年或更多。在 方框2330，将增加的时间与为仿真设置的初始时间比较。如果时间指 示还没有达到仿真的结束，那么过程流程返回到方框2316从而继续另 一迭代。如果仿真已结束，那么过程流程进展到方框2332。在方框2332， 仿真的结果存储在非瞬时计算机可读介质中，并可以向用户呈现，例 如通过用来计算物理储层的视觉表示。在结果呈现之后，在方框2334， 进程结束。

本技术在碳氢化合物采集中的使用

在示例性实施例中，方法2300可以用来在储层中定位碳氢化合物， 或调整源自碳氢化合物油气田的生产。这包括基于多级混合多尺度有 限体积模型生成储层仿真模型。源自油气田的碳氢化合物生产的控制 可以至少部分基于从储层模型获得的结果来调整。

调整源自油气田的碳氢化合物生产的控制可以包括改变注入压 力、将注入井转变成生产井、将生产井转变成注入井、钻探更多新井 到储层，或其任何组合。统计地质学测量值可以包括传递率、孔隙体 积、排放体积、在井之间或在井和表示储层一部分的单元之间的最小 累积反传递率、传递时间，或其任何组合。

示例性集群计算系统

图24是可以用于本技术的示例性实施例的示例性集群计算系统 2400的框图。图解的集群计算系统具有四个计算单元2402，其每个都 可以为仿真模型的部分执行计算。然而，本领域技术人员认识到由于 可以选择任何数目的计算配置，因此本技术不限于该配置。例如，微 小仿真模型可以在单个计算单元2402例如工作站上运行，而巨大仿真 模型可以在具有10个、100个、1000个以至更多的计算单元2402的 集群计算系统2400上运行。在示例性实施例中，计算系统2402的每 个都为单个子域或计算单元的群运行仿真。然而，计算单元2402的分 配可以用任何数目的方式执行。例如，多个子域可以分配到单个计算 单元2402，或多个计算单元2402可以分配到单个子域，取决于在每个 计算单元2402上的负载。

集群计算系统2400可以从一个或多于一个客户端2404经由网络 2406访问，例如通过高速网络接口2408访问。网络2406可以包括局 域网（LAN）、广域网（WAN）、互联网或其任何组合。客户端系统 2404的每个可以具有计算机可读存储器2410以便操作代码和程序的存 储，包括随机访问存储器（RAM）和只读存储器（ROM）。操作代码 和程序可以包括用来实施在此讨论的方法中的全部或任何部分的代 码，例如关于图23讨论。进一步地，非瞬时计算机可读介质可以保存 物理碳氢化合物储层的数据表示，例如由多级混合多尺度有限体积 （MMMFV）在非结构化网格上生成的该数据表示。客户端系统2404 也可以具有其它非瞬时计算机可读介质，例如存储系统2412。存储系 统2412可以包括一个或多于一个硬盘驱动器、一个或多于一个光盘驱 动器、一个或多于一个闪存驱动器、这些单元的任何组合，或任何其 它合适存储设备。存储系统2412可以用于代码、模型、数据和用来实 施在此描述的方法的其它信息的存储。例如，数据存储系统可以保存 至少部分使用多级混合多尺度有限体积（MMMFV）在非结构化网格 上生成的物理碳氢化合物储层的数据表示。

高速网络接口2408可以在集群计算系统2400中耦合到一条或更 多通信总线，例如通信总线2414。通信总线2414可以用来将指令和数 据从高速网络接口2408通信到集群存储系统2416，并通信到集群计算 系统2400中计算单元2402的每个。通信总线2414也可以用来在计算 单元2402和存储阵列2416之间通信。除了通信总线2414之外，高速 总线2418可以存在从而提高计算单元2402和/或集群存储系统2416 之间的通信速率。

集群存储系统2416可以具有一个或多于一个有形的、计算机可读 的介质设备例如存储阵列2420，以便存储数据、视觉表示、结果、代 码或其它信息，例如关于图23的方法的实施和源自该方法的结果的信 息。存储阵列2420可以包括硬盘驱动器、光盘驱动器、闪存驱动器、 全息存储阵列或任何其它合适设备的任何组合。

计算单元2402的每个可以具有处理器2422和关联的本地有形计 算机可读介质，例如存储器2424和存储2426。存储器2424可以包括 ROM和/或RAM，该ROM和/或RAM用来存储代码，例如用来引导 处理器2422实施在图23中图解的方法的代码。存储2426可以包括一 个或多于一个硬盘驱动器、一个或多于一个光盘驱动器、一个或多于 一个闪存驱动器，或其任何组合。存储2426可以用来为中间结果、数 据、图像或包括用来实施图23的代码的与操作关联的代码提供存储。

本技术不限于在图24中图解的架构或单元配置。例如，任何合适 的基于处理器的设备可以用来实施本技术的实施例的全部或一部分， 无限制包括个人计算机、网络个人计算机、膝上计算机、计算机工作 站、GPU、移动设备，以及带有（或没有）共享存储器的多处理器服 务器或工作站。此外，实施例可以在专用集成电路（ASIC）或超大规 模集成（VLSI）电路上实施。事实上，本领域技术人员可以利用能够 根据实施例执行逻辑操作的任何数目的合适结构。

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