CN104239661A - 一种大规模油藏数值模拟计算的方法 - Google Patents

Abstract

本发明提供了一种大规模油藏数值模拟计算的方法，属于油气田开发领域。本方法包括：(1)在主处理器上定义变量；(2)在主处理器上进行数据准备和输入；(3)主处理器将数据传递给各处理器；(4)主处理器对模拟区域进行网格剖分，得到模拟区域的网格剖分结果；(5)在主处理器上采用区域分解法将模拟区域分解成多个子区域；(6)主处理器给各处理器分配存储空间；(7)主处理器读取各块子区域数据和连接信息并分配数据，各处理器接收属于自己处理器的数据；(8)主处理器和其他各处理器之间进行数据交换。

Description

一种大规模油藏数值模拟计算的方法

技术领域

本发明属于油气田开发领域，具体涉及一种大规模油藏数值模拟计算的方法。

背景技术

大规模油藏数值模拟越来越制约着油藏数值模拟的发展，在1998年，Argonne国家实验室和Texas大学利用IBM SP成功地试算了400万网格点、3200万个未知数的超大规模油藏数值模拟问题。在中国，大庆油田有限责任公司勘探开发研究院于2000年左右引进了Orogin2000高性能计算机及一套ParallelVIP2000商业版并行软件，该系统上运行一次116万网格点规模、300口井、30年拟合历史数据，需要70小时以上；其他石油公司也陆续引进了几套ParallelVIP、Parallel Eclipse软件，由于国外并行机及并行软件的技术壁垒和高价政策，引进的软硬件系统难以有效地进行百万级网格点的数值模拟；中国科学院软件中心研究院开发了一套油藏数值模拟并行软件PRIS，并有一定的后处理并行显示能力，在计算机能力和效率方面，PRIS有着其他专业软件不可比拟的独特优势。2000年，在国家863项目的支持下，通过与大庆勘探研究的合作，PRISV1.1在Linux微机群(16P-III500+Ethemet100Mb)上完成了116万网格点规模、300口井、30年拟合历史数据的并行模拟，一次模拟历时64h，与大庆引进的国外商业并行机、商业并行软件性能相当；同年，PRISE V1.1在曙光3000高性能并行机上模拟同样的数据，一次模拟仅需要13.7h；2001年，PRISE V2.0在Linux微机群(16P-III500+Ethemet100MMb)上模拟一次116万网格点、291口井、31.5年历史拟合的数据为例，在LSSC-II万亿次量级的Linux微机群(256节点，512Xeon2GCPU+1G Memory/CPU，Myrinet2000)上，使用128CPU进行模拟，需要消耗时间0.77h(单节点单进程情形)或0.87h(单节点双进程情形)，每个处理器使用的内存不超过50Mb，程序运行一次所发生的I/O总量约为40Gb，其中读取文件的总量约占500Mb。就并行计算领域的科学研究而言，百万量级的数据模拟，需要硬盘总量约150Gb；如果百万量级的数据用于实际生产，需要为可视化程度提供大量数据，需要硬盘总量至少为300Gb。

2003年，中国科学院完成了基于Linux微机群的百万量级网格点、20年历史拟合的实际数据的有效模拟计算。目前千万量级网格点油藏数值模拟问题目前的巨大挑战。

近年来，随着储层地质建模技术的高速发展，综合了地震、测井、地质和油藏等众多学科的储层地质模型，其规模已经达到几百万甚至数千万个网格节点，基本到达了足够精细的程度。这种精细油藏模拟对于研究剩余油分布，促进我国相当部分已进入高含水后期老油田的高产、稳产具有重要意义，它的计算规模实际需要达到百万、千万量级，历史拟合时间一般超20年。然而，由于油田应用部门在高性能计算机、高效率并行软件使用上的滞后，由于生产实际对油藏数值模拟计算时效性的追求，目前油藏数值模拟的规模还只在10～30万个网格节点水平，远远达不到的精细的程度。正是由于受计算机能力的制约，高精度的地质模型必须进行粗化处理，才能用于油藏模拟，这导致大量地质信息的损失，达不到精细描述油藏流体分布的目的，使用油田中后期的开发调整方案及实施效果受到很大影响。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种大规模油藏数值模拟的计算方法，实现大规模油藏的快速模拟。

本发明是通过以下技术方案实现的：

一种大规模油藏数值模拟的计算方法，包括以下步骤：

(1)在主处理器上定义变量(因为是并行计算，分主处理器和其他处理器，所有处理器包括这两部分)；

(2)在主处理器上进行数据准备和输入；

(3)主处理器将数据传递给各处理器；

(4)主处理器对模拟区域进行网格剖分，得到模拟区域的网格剖分结果；

(5)在主处理器上采用区域分解法将模拟区域分解成多个子区域；

(6)主处理器给各处理器分配存储空间；

(7)主处理器读取各块子区域数据和连接信息并分配数据，各处理器接收属于自己处理器的数据；

(8)主处理器和其他各处理器之间进行数据交换；

(9)在每个处理器上设置各自的局部方程，并分别采用牛顿迭代法求解所述局部方程，求解得到各个子区域上的物理参数的值；所述局部方程包括压力求解方程和饱和度求解方程；

(10)在所有处理器上利用步骤(9)获得的物理参数的值计算边界处及共有的物理参数(根据边界条件或共有网格附近参数利用压力、饱和度进行计算)，并更新物理参数(即利用新的物理参数的值替换旧的值)，这些物理参数不仅体现在局部方程中，还体现在全局方程中；全局方程是指所有局部方程的总和，全局是指整个模拟计算区域，局部是指每个处理器的模拟计算区域；

(11)判断压力求解方程和饱和度求解方程是否收敛，如果是，则转入步骤(12)，如果否，则返回步骤(8)；

(12)判断时间t是否大于计算时间T_max，如果是，则转入步骤(13)，如果否，则返回步骤(8)；

(13)将从压力求解方程和饱和度求解方程计算得到的结果(主要包括压力场、饱和度场等数据)输出到主处理器上；

(14)结束。

所述步骤(1)中的变量包括油压、油饱和度、气饱和度或饱和压力。

所述步骤(2)中的数据包括地层参数、油气参数、模拟区域的剖分网格和计算参数(包括相对渗透率、毛管力等)。

所述步骤(6)中采用共享存储、分布式存储和混合分布式存储这三种存储结构中的一种来给各处理器分配存储空间。

所述步骤(7)中分配数据是指将子区域与处理器一一对应，即每个处理器计算一个模拟子区域的数据，并且将网格均匀分配到各个处理器上，且子区域边界上的网格连接数最小。

所述步骤(7)中的连接信息是指两个相邻网格之间的正交距离和接触面积。

所述步骤(9)中的压力求解方程和饱和度求解方程是根据子区域中各个网格所对应的质量与能量守恒方程得到的。

所述步骤(10)中的物理参数包括流体的粘度、密度、流度、压缩系数。

所述步骤(11)中判断压力求解方程和饱和度求解方程是否收敛是这样实现的：判断两次迭代后压力求解方程的结果的差值和两次迭代后饱和度求解方程的结果的差值是否小于油藏模拟允许的误差限，如果两个差值中有一个小于误差限，则判断为不收敛，否则判断为收敛。

所述步骤(12)中的计算时间是指模拟区块开发时间。

与现有技术相比，本发明的有益效果是：实现了大规模油藏的快速数值模拟，在大规模油藏数值模拟并行计算中取得了超加速比。本发明应用领域为大规模油藏的开发和该地质结构的地下水的开发等领域。该发明能开展大规模油藏(如大油田)的整体精细数值模拟，提高了模拟效率，必将推动该类油藏数值模拟的进步。

附图说明

图1，大规模油藏数值模拟并行计算流程；用于确定软件编程设计方案，开展油藏数值模拟并行计算研究与应用。

图2，大规模油藏数值模拟并行计算共享存储；与图3、图4的存储结构进行比较，确定大规模油藏数值模拟并行计算采用的存储结构。

图3，大规模油藏数值模拟并行计算分布式存储；与图2、图4的存储结构进行比较，确定大规模油藏数值模拟并行计算采用的存储结构。

图4，大规模油藏数值模拟并行计算混合分布式共享存储；与图2、图3的存储结构进行比较，确定大规模油藏数值模拟并行计算采用的存储结构。

图5-1，模拟区域的分解

图5-2，模拟区域的连接信息的压缩存储；说明区域分解如何在大规模油藏数值模拟并行计算中的应用以及连接信息的压缩存储格式。

图6-1，跨CPU通信方式之一。说明处理器间的数据通信方式。

图6-2，跨CPU通信方式之二。

图7，网格剖分实例。

图8，并行计算的加速比曲线。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明提出了一种新的进行大规模油藏数值模拟的新的方法。大规模油藏数值模拟目前计算速度慢，甚至无法计算，为此发明了本方法。本发明的技术思路是：大规模油藏数值模拟的并行计算流程，优选实现并行计算的存储结构，然后提出并行计算的任务分解方法，继而进行组装雅克比矩阵，并给出了线性方程组的并行解法，也对计算过程中的通信方式进行了说明，最后通过实施例验证了方法正确性，实现了大规模油藏的快速模拟。

大规模油藏数值模拟计算采取并行计算方法，其并行思路如图1所示，首先是在主处理器上定义计算中的主变量(主变量包括油压，油饱和度，气饱和度或饱和压力)、进行数据准备和输入(根据软件(如中石化石油勘探开发研究院的Karstsim软件)的输入数据格式，准备输入文件，包括地层参数，油气参数，模拟区域的剖分网格，计算参数等)，并将数据传递给各处理器，进行网格剖分(通过其他预处理软件(如美国XYZ，Inc.公司的网格剖分软件TrueGrid))对模拟区域进行网格剖分，具体剖分方式根据不同的问题而定(如对于缝洞型油藏，其存在溶洞、裂缝和基质，对溶洞部分应采用六边形网格，可以设置网格的排列方向。对于裂缝一般为细长矩形网格，对于基质可以采用任意网格，如图7所示的网格剖分实例，该模型采用非规则网格进行剖分，基质和裂缝采用四边形网格，溶洞采用六边形网格。网格总数16100个，约4万个链接。)，采用区域分解法分解计算区域(采用METIS软件包，根据各子区域的计算量、通讯量基本平衡，实现最少通讯量的综合考虑的原则，对模拟区域进行自动分解为子区域)，给各处理器分配存储空间，读取各块子区域网格数据和连接信息(指两个相临网格之间的正交距离，接触面积等信息)并分配数据(分配数据是指将单元区域与处理器的对应。)，这期间其他处理器接收属于自己处理器的数据，然后主处理器与其他处理器进行数据交换(每个子区域的计算需要其相邻子区域的信息才能进行，比如位于不同子区域相连接的两网格在各自的子区域中计算时都需要对方的主变量值。)，在每个处理器上设置各自的局部方程(子区域中各个网格所对应的质量守恒方程(如方程(1)-(3))，并在该处理器上采用牛顿迭代法求解方程(对每个处理器的局部方程进行求解)，更新物理参数(次要变量，如流体的黏度，密度等)，判断是否收敛(根据给定的标准判断是否已收敛，标准为主变量(压力、饱和度)两次迭代得到的差值小于油藏模拟允许的误差限)，如果不收敛继续进行处理器间通信，如果收敛，判断是否达到计算时间(这根据模拟区块开发时间，如果模拟区块开发时间为10年，这计算时间就设置10年；如果要预测20年，看看20年后的采收率、剩余油分布情况，这时计算时间就是20年)，未达到则进入下一个时间步计算，达到了则把计算结果输出到主处理器上。

质量守恒方程

气相： $\frac{\∂}{\∂t}\left\{\mathrm{\φ}(S_o{\stackrel{\‾}{\mathrm{\ρ}}}_{\mathrm{dg}}+S_g{\mathrm{\ρ}}_g)\right\}=-\▿\·({\stackrel{\‾}{\mathrm{\ρ}}}_{\mathrm{dg}}{\stackrel{\→}{V}}_o+{\mathrm{\ρ}}_g{\stackrel{\→}{V}}_g)+q_g---\left(1\right)$

水： $\frac{\∂}{\∂t}\left(\mathrm{\φ}{S}_w{\mathrm{\ρ}}_w\right)=-\▿\·\left({\mathrm{\ρ}}_w{\stackrel{\→}{V}}_w\right)+q_w---\left(2\right)$

油： $\frac{\∂}{\∂t}\left(\mathrm{\φ}{S}_o{\stackrel{\‾}{\mathrm{\ρ}}}_o\right)=-\▿\·\left({\stackrel{\‾}{\mathrm{\ρ}}}_o{\stackrel{\→}{V}}_o\right)+q_o---\left(3\right)$

这里，β相(β＝g为气，＝w为水，＝o为油)的达西速度如下定义：

${\stackrel{\→}{V}}_{\mathrm{\β}}=-\frac{{\mathrm{kk}}_{\mathrm{r\β}}}{{\mathrm{\μ}}_{\mathrm{\β}}}(\▿P_{\mathrm{\β}}-{\mathrm{\ρ}}_{\mathrm{\β}}g\▿D)$

ρ_β是β相在油藏条件下的密度；

是在油藏条件下脱去溶解气的油相密度；

是在油藏条件下油相中溶解气的密度；

φ是油层的有效孔隙度；

μ_β是β相的粘度；

S_β是β相的饱和度；

P_β是β相的压力；

q_β是地层β组分每单位体积汇点/源点项；

g是重力加速度；

k是油层的绝对渗透率；

k_rβ是β相的相对渗透率；

D是深度。

2)并行计算存储结构

KARSTSIM-MP(并行版油、气、水多相地下流体模拟器)可应用于共享存储(如图2所示)、分布式存储(如图3所示)和混合分布式存储(如图4所示)的存储结构。共享存储的存储方式允许多台服务器访问(写)同一个分区，实现各台主机可相互高效率的共享数据和传递数据，并实现数据保存和保护的功能。分布式存储通过网络使用不同角落机器上的磁盘空间，并将这些分散的存储资源构成一个虚拟的存储设备，数据分散的存储在各个角落。混合分布式共享存储集共享存储和分布式存储于一体的存储方式。本发明方法可使用于所有的这些存储方式。

3)网格剖分和区域分解

形成一套有效的和高效率的非结构网格剖分方法对于并行计算的成功是关键的一步。并行计算的区域分解法如图5-1和图5-2所示，首先，为了实现更好的数值模拟性能，并行模拟需要将网格均匀分配到不同的处理单元或处理器，也就是说，分配到每个处理单元(处理单元就是指除主处理器外的其它处理器)的网格数大致相同；第二，区域边界上的网格连接数达到最小，实现通讯量的平衡和最少(采用METIS软件包所提供的算法来实现的)，这样做的目的就是减少处理器相互通信所花费的时间。

对于KARSTSIM-MP模拟，模型区域采用一组一维、二维或三维的网格块和网格连接代表的界面信息(相连两网格中心点到其接触界面的距离及两网格之间的接触面积)(即前述的连接信息)进行表征(如图5-1和表1所示)。网格系统视为非结构网格，并通过网格的连接信息构造出网格连接信息的数组，模型网格的连接信息通过压缩进行存储。

表1

KARSTSIM-MP任务并行原则就是计算负载平衡和通信量最小化，即实现最小的平均(或总的)通信量。KARSTSIM-MP采用METIS包进行分区(采用一种特殊有效的优化算法，实现计算量、通讯量平衡并最小化，具体内容见参考文献(美国明尼苏达大学计算机系George Karypis和Vipin Kumar1998年编写的“非结构图形分区、网格分区、方程组计算分组软件包”技术报告)，同时还可以根据处理器数量进行自动分区或者根据具体条件人工分区。整个模拟区域被离散成很多网格块，如果有n个处理器用于模拟，那么这些网格块或模拟区域被分成n个部分(也就是图5-1所示的3个处理器对应3个部分)。

4)组装雅克比矩阵

KARSTSIM-MP采用积分有限差分法处理空间离散，形成一组强耦合的非线性代数方程组，并采用牛顿迭代法进行线性化处理。在每步牛顿迭代步中，雅可比矩阵初次通过数值差分进行组装，然后利用不同的预处理程序对线性方程系统进行线性迭代求解。

每个处理器内部的雅克比矩阵以分布式可变块行(DVBR)格式存储(指的是矩阵存储格式，DVBR是一种科学计算通用的稀疏矩阵存储方式，与图2、3、4无关)。对角块(矩阵是由“块”组成的，每个“块”由三行三列组成的。对角块指位于对角线上的那些“块”)首先存储在每个块行(一行一行逐块存储)上，然后所有的矩阵块逐行存储。每个矩阵块的标量元素(块中的各个元素)以列优先级存储(一列一列逐列存储)。雅克比矩阵的数据结构由一个实型向量和五个整型向量组成。

5)线性方程组的并行解法

利用AZTEC软件包(AZTEC只用于解线性方程组，用牛顿迭代法对非线性方程线性化后，解线性方程中用到)提供的线性迭代方法求解每个牛顿迭代步形成的线性方程系统(即图1中的“所有处理器：采用牛顿迭代法求解方程”)。AZTEC软件包提供几种不同的求解器和预处理器供用户选择。求解全局线性方程系统是由所有的处理器共同分担，即每个处理器负责处理各自区域的方程。为了实现并行的解决方案，处理器之间通信要快捷顺畅，以保证相邻网格之间的数据交换。此外，全局通信对于计算规范化向量以检验计算的收敛情况也是非常必要的。

在并行计算过程中，时间步长会根据迭代的收敛速率进行自动调整。对于KARSTSIM-MP，第一个处理器(主处理器，名为PEO)收集来自其他处理器的必要数据后，就会计算时间步长。不同处理器迭代收敛的速率不同，只有当所有的处理器达到收敛准则以后才进入下一个时间步长的计算。

6)跨CPU通信

跨CPU通信包括(图6-1和图6-2所示)相邻CPU的通信和主CPU和其他CPU之间的通信。线性方程的求解和组装雅可比矩阵结束时需要用到相邻CPU的通信(即图1中的“所有处理器：更新物理参数”)，而模拟前初始条件和参数的分配、计算结果输出以及需要全局变量时就会用到主CPU和其他CPU之间的通信。

并行算法的核心部分就是将模拟区划分为不同的区域，然后将这些区域分配到不同的处理器，不同区域之间又通过网格块之间的链接来实现跨CUP通信。此外，对于全局通信，需要考虑所有处理器的贡献值，并计算向量的规范性，以检查并行计算的收敛情况。除了用以解决线性方程系统的线性求解器例程之间的通信，相邻处理器之间的通信对于及时更新主变量是非常必要的。当处理器调用子程序时，就会执行交换相应网格组的向量元素。在时间步和牛顿迭代过程中，外部变量的交换对于包含初始变量的向量组是必要的。此外，本发明将“无堵塞”通信引入到Aztec软件包和牛顿迭代，以更好地改善通信方案，“无堵塞”通信是通过MPI“无堵塞”通讯函数来实现的。

为了验证该方法的正确性，设计一均质三维模型，模型大小为100m×100m×100m。将其离散成100×100×100网格。无限大底水(储层底部下面为含水层，认为是“无限大的水体)，中心一口生产井，X、Y、Z方向上的渗透率均为30mD，孔隙度为0.2，顶部深度为1m，岩石压缩系数为1.0×10^-51/Pa，原始地层压力20MPa，采油速度为25方/天。利用KARSTSIM-MP分别采用不同处理器进行模拟，计算10年，获得的并行计算的特性如表2所示，并行计算的加速比曲线如图8所示，可见并行计算大大提高了计算速度，最终获得了超加速比，即并行计算的加速比(就是图8中的总的加速比)大于理想加速比。

表2

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (10)

1.一种大规模油藏数值模拟计算的方法，其特征在于：所述方法包括以下步骤：

(1)在主处理器上定义变量；

(2)在主处理器上进行数据准备和输入；

(3)主处理器将数据传递给各处理器；

(4)主处理器对模拟区域进行网格剖分，得到模拟区域的网格剖分结果；

(5)在主处理器上采用区域分解法将模拟区域分解成多个子区域；

(6)主处理器给各处理器分配存储空间；

(7)主处理器读取各块子区域数据和连接信息并分配数据，各处理器接收属于自己处理器的数据；

(8)主处理器和其他各处理器之间进行数据交换；

(9)在每个处理器上设置各自的局部方程，并分别采用牛顿迭代法求解所述局部方程，求解得到各个子区域上的物理参数的值；所述局部方程包括压力求解方程和饱和度求解方程；

(10)在所有处理器上利用步骤(9)获得的物理参数的值计算边界处及共有的物理参数，并更新物理参数；

(11)判断压力求解方程和饱和度求解方程是否收敛，如果是，则转入步骤(12)，如果否，则返回步骤(8)；

(12)判断时间t是否大于计算时间T_max，如果是，则转入步骤(13)，如果否，则返回步骤(8)；

(13)将从压力求解方程和饱和度求解方程计算得到的结果输出到主处理器上；

(14)结束。

2.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(1)中的变量包括油压，油饱和度，气饱和度或饱和压力。

3.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(2)中的数据包括地层参数、油气参数、模拟区域的剖分网格和计算参数。

4.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(6)中采用共享存储、分布式存储和混合分布式存储这三种存储结构中的一种来给各处理器分配存储空间。

5.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(7)中分配数据是指将子区域与处理器一一对应，并且将网格均匀分配到各个处理器上，且子区域边界上的网格连接数最小。

6.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(7)中的连接信息是指两个相邻网格之间的正交距离和接触面积。

7.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(9)中的压力求解方程和饱和度求解方程是根据子区域中各个网格所对应的质量与能量守恒方程得到的。

8.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(10)中的物理参数包括流体的粘度、密度、流度、压缩系数。

9.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(11)中判断压力求解方程和饱和度求解方程是否收敛是这样实现的：判断两次迭代后压力求解方程的结果的差值和两次迭代后饱和度求解方程的结果的差值是否小于油藏模拟允许的误差限，如果两个差值均小于误差限，则判断为收敛，否则判断为不收敛。

10.根据权利要求1所述的大规模油藏数值模拟计算的方法，其特征在于：所述步骤(12)中的计算时间是指模拟区块开发时间。