CN101661514B - 一种油藏黑油模型数值模拟系统 - Google Patents

Abstract

本发明涉及一种油藏黑油模型数值模拟系统，属于石油开采技术领域，其结构包括前处理子系统、串并行模拟器和后处理子系统，首先是根据油藏的地质数据、生产动态数据、流体性质数据、岩石渗流特性数据，利用系统的前处理子系统，建立地质模型、动态模型、流体模型和岩石渗流模型，进而形成数值模拟模型，然后利用系统的高效并行模拟器，进行模拟运算，再现油田的开发过程，得到不同时间的油藏模型，再利用系统的后处理模块，以可视化的形式，分析模拟运算的结果及剩余油分布，优化设计油藏进一步开发的调整方案，最后提交当前条件下的最佳方案进行现场实施。

Description

一种油藏黑油模型数值模拟系统

(一)技术领域

本发明涉及一种应用于石油开采领域的油藏黑油模型数值模拟系统，能够提高油藏数值模拟的运算速度和并行效率。

(二)背景技术

油田开发是一项隐蔽的技术工程，具有高投入、高产出、高风险的特点和不可逆性，科学决策尤为重要，油藏数值模拟是数字化油藏、寻找剩余油、优化设计开发方案的主要技术手段，在油田开发中发挥着极为重要的作用。近几年来，随着油田开发的深入，大多数老油田都进入了高含水期和特高含水期，剩余油分布日趋零散、复杂，油藏数值模拟进行剩余油分布研究的主要目标集中在井间、层间和层内，这就要求油藏数值模拟研究向大规模、精细化方向发展，剩余油分布研究向精确化、定量化、动态化方向发展。

目前国际流行的三维三相黑油模型商业化软件，尽管其功能强大，稳定性好，但在实际应用中仍然存在许多不足。一方面胜利油区以多期构造运动形成的断陷湖盆沉积为主、不同时期、不同区域断裂活动差异明显，储层特征呈复杂多样的特点。开发生产过程中，研究人员针对不同油藏对软件提出了不同的需求，但由于引进商业软件只有可执行码，应用过程中的一些现场需求无法通过修改软件得到满足，因而也限制了油藏数值模拟的研究水平。研制自主版权的油藏模拟并行软件可以按需灵活修改软件，扩大应用领域。另一方面，引进商业软件不可避免的受到版权的限制。由于价格昂贵，购买的软件许可较少，难以展开全面的应用。许多软件在引进的同时，还被限定了应用领域，往往被要求不得服务于国外竞争对手的招标项目。这些限制障碍了油藏数值模拟技术的发展和大规模应用，更制约了中石化参与国际竞争的总体战略部署。另外，我国的石油地质工作者，在长期的工作实践中，总结和摸索出了一整套行之有效的适应国内油藏特点的技术方法和工作经验。通过自主版权数值模拟并行软件的研制，可以将其标准化、定量化，通过不断积累形成具有国内技术特色的工具软件，发展成为油藏研究人员的得力助手。

(三)发明内容

本发明的技术任务是针对现有技术的不足，提供一种油藏黑油模型数值模拟系统，该系统提高了油藏黑油模型数值模拟整体并行效率，解决了超大规模的油藏模拟问题，而并行效率的提高又涉及到多方面的技术问题，如：线性方程组解法的并行化问题，软件的并行度问题；整体并行中的程序结构调整问题，整体并行中的数据通讯问题，整体并行中的输入输出问题，油藏的高度非均质问题，软件对微机集群的适应性问题。

本发明解决其技术问题所采用的技术方案是：

一种油藏黑油模型数值模拟系统，包括前处理子系统、串并行模拟器和后处理子系统，其特征在于：根据前期油藏描述的成果，由前处理子系统建立油藏黑油模型数值模拟模型，然后利用串并行模拟器，进行模拟运算，再利用后处理子系统分析模拟的结果，设计优化油藏开发调整方案，最后形成文字报告和汇报多媒体，提交最佳设计方案。

上述油藏黑油模型数值模拟系统，其前处理子系统用于建立数值模拟模型，包括项目管理、方案管理、数据管理、模型定义、网格建立、饱和度函数、流体性质、分区设定、水体描述、动态模型建立、作业控制和输出控制12个模块，项目管理实现项目相关信息的统一组织与管理，具体包括项目建立、项目打开、保存、另存、项目删除、退出6项功能；方案管理实现项目不同研究方案的管理，具体包括新建方案、输入方案描述信息，拷贝方案、保存方案、删除方案5项功能；数据管理实现项目涉及数据的统一管理，具体包括数据输入、数据检查、数据编辑、数据显示、数据分析、数据保存、数据输出多项功能；模型定义，定义所建模型的基本信息，网格建立用于生成表征油藏地质信息的网格模型，具体包括数据加载、编辑、场图插值、网格设计、网格编辑、网格优化、网格赋值、三维图形显示多项功能；饱和度函数主要是建立岩石流体的饱和度函数及毛管压力函数；具体包括相对渗透率曲线和毛管压力曲线导入、曲线显示、曲线编辑以及曲线处理和数据输出；流体性质主要是处理分析油气水的实验分析数据，建立描述流体性质的流体模型；分区设定主要是设定不同饱和度函数、流体性质以及储量统计对应的分区；水体描述用于表征油藏的水体分布及大小；动态模型建立主要是通过油藏开发历史数据的分析与处理，建立反映油藏开发过程的动态模型；作业控制用于设置模型运算的迭代误差控制、时间步长控制；输出控制用于设置输出的内容和频率。

上述油藏黑油模型数值模拟系统，其串并行模拟器以数值模拟模型作为输入文件，进行模拟运算，串并行模拟器的实现包含了以下关键技术：可重叠的空间子区域划分方案、负载平衡优化设计、动态分配内存、网格重排方法、区域分解中分割数据的快速索引、加法Schwarz方法、并行解的邻域校正方法、曲面预估方法、线性方程组求解参数及方法优化、关键技术包括数据压缩通讯、压缩存储技术；这些关键技术的应用提高了模拟运算的速度，解决大规模精细数值模拟研究问题。

上述油藏黑油模型数值模拟系统，其后处理子系统用于分析模拟结果，设计优化油藏开发方案，包括自定义参数场计算、平面视图绘制、曲线图绘制、三维立体视图绘制、油藏纵向切片图绘制、注采流线图绘制、四位一体的油藏多指标综合图绘制、油藏动静态参数定量动态描述8个模块，自定义参数场计算利用已有的常规动静态参数场进行组合计算，生成新的有地质意义的数据场，用户可自定义最多10个计算参数场；平面视图绘制具有两大特色；一是图件绘制齐全且有创新：不仅可以绘制彩色分布图、填充等值线图、等值线图、含水柱状图、累油水饼图，还可以绘制叠加等值线图；二是图件绘制参数多样：包括油藏静态数据、动态数据、原始数据及自定义场数据四类；曲线图绘制可以绘制区块及油水井的累积、瞬时曲线图，分析不同时期区块及油水井的生产状况，设置了灵活的绘图工具：一是绘图的小层、井号及参数可以任意选择；二是绘图参数的X、Y坐标可以任意选择，曲线的条数、绘图的时间段可以任意选择；三是曲线图中的字体大小、字型、曲线的线型、颜色等均可自行选择，方便了用户；三维立体视图绘制包括三维立体彩色分布图、三维立体等值线图两种图件的绘制；通过三维立体图的绘制，可以直观地反映油藏在平面和纵向上的变化特征；油藏纵向切片图绘制可以自动绘制静态数据、动态数据、原始数据及自定义场数据四类参数的油藏切片图；通过油藏切片图的绘制，可以直观地反映油藏在纵向上的变化特征，搞清每口油井钻遇的油层数、有效厚度、含油饱和度、剩余储量等基本数据，同时可以直观地描述剖面线上油藏的非均质变化特征；注采流线图绘制可以根据当前时间压力场绘制注采流线分布图，显示流体流动的方向，分析注采对应状况；根据注水流线的分布图，可以进行注采井网完善状况的判断与分析，确定那些区域是优势流场，那些区域是非优势流场，确定剩余油潜力范围和挖潜方向；四位一体的油藏多指标综合图绘制可以在同一微机屏幕上同时进行4种油藏地质图件的显示和对比分析，从而多角度、多方面的对油藏进行综合分析；油藏动静态参数定量动态描述把油藏储层参数和剩余油描述指标有机地结合起来，在微机上形成细化的网格数据体，全方位、多指标的对油藏动静态参数进行定量描述，并通过时间坐标的设置，实现油藏二维、三维图件的动态显示，通过不断补充数据实现对油藏的跟踪分析和数据化管理。

本发明针对上述技术问题，制定了相应的措施，通过优选线性方程组的求解方法，提高线性方程组算法的可并行性；采用整体并行，提高软件的并行度；调整程序结构，提高软件结构的可并行性；采用压缩通讯、压缩存储技术，减少数据的交换量；采用专门输入输出结点，减小输入输出的影响；建立灵活分区模式，减小负载的不平衡；建立线性方程组并行、整体并行两种模式，提高软件对机器的适应性。

整体并行方案主要是基于区域分解思想，采用多种技术与方法实现软件的整体并行。主要步骤如下：

第一步：区域分解。通过灵活的二维平面上的区域划分，将网格数据和井数据划分到不同的子区域，即并行处理前数据的准备。关键技术包括可重叠的空间子区域划分方案、负载平衡优化设计。

第二步：分区数据的预处理。在区域分解的基础之上，对各个分区的网格、井数据进行并行处理，设置本区域相应的参数。关键技术包括动态分配内存、网格重排方法、区域分解中分割数据的快速索引。

第三步：构建线性方程组。各分区按照本分区网格、井数据的特征，通过相应的计算，并行形成各自的渗流方程组系数矩阵和子区域边界条件更新。

第四步：线性方程组的并行求解。首先，通过初始压力预估技术，根据每个网格在此之前的压力变化情况，预估一个压力初值，从而使初始解曲面更接近于最终精确解；其次，通过线性方程组求解技术，实现分区线性方程组的求解；第三，对分区得到的解分别进行收敛性判断；第四，对于收敛后的解，进行边界网格压力后加速，提高牛顿迭代收敛速度。最后，对整个区域进行牛顿迭代收敛性判断。关键技术包括：加法Schwarz方法、并行解的邻域校正方法、曲面预估方法、线性方程组求解参数及方法优化。

第五步：结果收集整理和输出。利用一个0号进程为根进程负责数据的收集和输出，关键技术包括数据压缩通讯、压缩存储技术。

本发明的油藏黑油模型数值模拟系统与现有技术相比，所产生的有益效果是：黑油模型数值模拟整体并行技术的应用，能够大幅度地提高模拟器的运算速度和模拟能力，在2个CPU情况下，并行效率超过80％，在16个CPU情况下，并行效率也可以达到50％以上，对于200万节点的油藏模型，模拟时间控制在8小时以内。

(四)附图说明

附图1为油藏数值模拟整体并行程序框架示意图；

附图2为油藏模型32×32×20网格系统被分成四个无重叠子区域示意图；

附图3为油藏模型32×32×20网格系统被分成四个有重叠层的子区域示意图；

附图4为并行油藏数值模拟边界流量交换和方程组系数矩阵形成过程示意图；

附图5为并行油藏数值模拟子区域内部网格间流量交换过程示意图；

附图6为并行油藏数值模拟子区域边界层网格间流量交换过程示意图；

附图7为并行油藏数值模拟子区域的上下左右边界上流量交换后过程示意图；

附图8为非对角线元素的排列示意图；

附图9为子区域内部网格与前面的网格流量交流示意图；

附图10为子区域内部网格和左侧边界网格流量交换对方程组系数矩阵的影响示意图；

附图11为子区域内部网格与后面的网格流量交流示意图；

附图12为子区域内部网格和右侧边界网格流量交换对方程组系数矩阵的影响示意图；

附图13为邻域校正加速策略示意图；

附图14为数据收集及输出示意图；

附图15为油藏黑油模型数值模拟系统的整体结构图。

(五)具体实施方式

1.区域分解

在区域分解中，需要将整个油田模拟区域的数据按照一定的分配方法进行分割，每个子进程对分到的每一小块区域进行计算，如果各个子进程之间计算的区域有重叠，就称作重叠型的区域分解方法。主要采用了2项关键技术：

(1)可重叠的空间子区域划分方案

考虑到油气田开发中注水井和生产井的存在，为了降低三维问题渗流方程离散化的复杂度，区域划分方案使用X-Y平面上的二维分割，深度方向不作离散分割。考虑到区域分解后的负载平衡，同时为后续可重叠网格策略提供足够的可扩展性，采用可重叠的空间子区域划分方案。该方法采用了由参数文件提供区域划分方案的形式，增加了程序的灵活性，降低了程序运行时的额外开销。用户能够在数据文件中灵活指定区域划分方法，不仅可以控制在不同方向上分区数目，还可以非常简单的实现子区域间边界网格的任意层数重叠。例如，空间分解参数定义行(1)代表如图2的空间分解方案，子区域间无重叠网格；参数命令行(2)代表图3的空间分解方案，x方向上子区域间有x＝17，18两层重叠网格。

参数命令行(1)：

NXP 2

XP_BOUND 1161732

NYP 2

YP_BOUND 1161732

参数命令行(2)：

NXP 2

XP_BOUND 1181732

NYP 2

YP_BOUND 1161732

(2)负载平衡优化设计

负载平衡问题是并行化处理的首要问题，由于层间地质上的非均质性、动用程度的差异造成在不同时间步、迭代子时间步，甚至是同一时间步，不同的层计算量存在很大差别。同时按多级并行策略的思想，负载平衡问题是整个并行问题的关键。动态负载平衡就是按任务进行状况动态地分配任务。

油藏的非均质性是造成负载不平衡的主要根源，而非均质性不仅表现在层间也表现在层内，不仅来自油藏自身也有后期开发的原因。油藏的大小、形态、渗流能力、动用状况等因素都是造成负载不平衡的原因，解决好负载平衡问题是提高并行效率的又一关键。

从油藏模拟问题看，负载平衡的实现可以有两种方式：一是静态负载平衡，任务开始各PE就分配固定的任务量；二是动态负载平衡，任务的分配按各PE任务的执行情况而定。动态负载平衡的实现需要完全的主控节点，即主控节点只完成任务调度、监控工作。

区域分解是串行软件并行化中常用的方法，它按照一定规则将模拟区域划分成各个子区域，每个子区域分配给一个进程进行求解。

在实际油藏模型中，模拟区域中可能有些网格的渗透率为0(比如岩石网格)，这些网格点上没有流量，也不需要参与计算，这类网格称为无效网格(Invalid Grids)。无效网格在整个模拟区域中的分布是不规则的，因为无效网格不会参与计算，它会影响并行程序的负载分配。井数组在线性方程组的求解中占有一定比重，井也是影响负载的因素。

为了提高程序的并行性能，区域分解技术要遵循以下几条规则：

①分配给各个进程的计算规模要大体相当，即要达到一个好的负载平衡状态；

②区域分解带来的通信开销要尽可能小；

③区域分解造成的系数矩阵缺失要尽量少，避免区域分解后，各个子区域形成的线性方程组性态发生严重改变。

按照上述原则，区域分解方法仅在X，Y方向对模拟区域进行了规则分割，且分割不在Z方向上进行。如果在Z方向进行分割，会将同一口井划分给不同的进程进行计算，迭代过程中不同进程需要频繁交换井的计算结果，大大增加了通信开销，同时也使线性方程组的性态发生改变，另外，采用规则的分割方法也有利于有效地组织和实现并行程序的进程通信，避免由于不规则分割带来的复杂的进程间通信管理。

由于整个计算过程中各个分区的计算负载相对稳定，区域分解算法采用了静态的分解模式，仅在刚读入数据场时进行区域分解，数据场一旦分割之后每个进程的负载就不再改变，这样减少了由于动态网格划分带来的迁移开销。在任务划分算法中，每个进程均有一个负载评价函数E(gridNum)，gridNum为该进程分配到的子区域的有效网格数，任务划分负载平衡的目标就是使各个子区域评价函数E(gridNum)的方差最小。没有把每个进程计算的井数目加入负载评价函数中，这是因为在实际的油藏模拟中，井矩阵计算形成的负载同有效网格计算形成的负载相比是次要因素，负载评价函数中对井的省略不会对任务划分的负载均衡性产生较大影响。

带负载平衡的区域分解算法形式化描述如下：

①确定目标不均衡因子F₀

②按照几何上的均匀分割对模拟区域进行初始划分

③计算各个进程负载的评价函数E(gridNum)，得到各进程评价函数的均值

和方差D，计算不均衡因子 $F=\frac{D}{\stackrel{\‾}{E}}$

④如果F＜＝F₀，转⑥

⑤找到评价函数E(gridNum)同

相差最大的那个进程，在其最大不平衡方向进行一维搜索寻找优化划分，转(3)

⑥找到了满足负载平衡的分区方式，按照此方式建立索引，将全局数据场分配给各进程计算。

2.数据预处理

(1)网格及井数据划分并行实现

①网格数据划分并行实现

网格相关的数据结构的分割有三部分内容：所有网格的数据分割；子区域新加边的数据构造；水区网格的数据分割。涉及的数据结构主要是在ARRCOM.TYT、AQFCOM.TYT、MSTCOM.TYT中定义，ARRCOM.TYT和AQFCOM.TYT中主要是一些数组，后者主要是程序运行中的一些参数。

数组分为两类，一类是记录每个网格点数据的，比如厚度th(lbmax)，压力p(lbmax)等，这些保存每个网格数据的数组可以直接把数组元素从全局的编号映射到局部的编号，另一类是存和区域相关信息的数组，比如每个网格点到活动网格点的映射map(lbmax)，每个网格点的排序后的编号mno(lbmax)等，这些数组中的元素在区域分解之后要重新计算。通过debug跟踪，第一次进入iter34时，需要分割和改变的数据结构见表1：

表1分割和改变的数据结构

在程序代码中，网格数据的分割都是setdomain中完成的，setdomain的接口如下。

输入：

网格的映射信息，包括loc2glob_map，left_edge，right_edge，top_edge，bottom_edge，is_left，is_right，is_top，is_bottom

子区域边界的映射信息left_border_map，right_border_map，top_border_map，bottom_border_map

水区网格的映射信息loc2glob_aqu_map

输出：

数据分割后增加的边界上的数据结构left_border_**，right_border_**，top_border_**，bottom_border_**

输入输出：

pvf，sf.....等数组，要分割的数组，输入时为全局数组，输出时为分割后本子区域的数组；

map，mno，输入时为全局的网格映射关系表，输出是本子区域的网格映射关系表；

memst，MSTCOM中的输入是全局参数，输出是本子区域修改了相关参数后的值；

mmab，水区网格编号到全局网格(包括死节点)的映射；

maqubl，当前所有水区的网格数的累积。比如第一个水区有a个网格，第二个水区有b个网格，则maqubl(1)＝a，maqubl(2)＝a+b。

②井数据划分并行实现

在每个时间段，程序都会从数据文件sim2.dat中读取井的数据，通过一系列处理和转换存储在相关数据结构中，因此，在每次进入一个时间段时，都要把整个区域的井数据进行分割，每次结束一个时间段之前，都要收集每个进程的井数据，构成全局的数据。

与井相关的数据结构都是在WELCOM.TYT中定义的，其中与区域分解相关的数据结构分为三类：井的位置信息，如射孔位置等；与井相关的数据，如井的产率等；井的指示信息，如是否独立井等。其中，井的位置信息在区域分解之后要重新计算，而井的信息和数据数组中的元素，要从全局位置映射到局部位置。数据跟踪的结果见表2：

                 表2数据跟踪的结果

井数据的分割都在welldiv中完成，welldiv的接口如下。

输入：

left_edge，right_edge，top_edge，bottom_edge，每个子区域的左右上下四个边界位置；

local_map，本子区域的map数组，也就是在调用gridmap之后形成的map数组。

输出：

well_proc，每口井对应的子区域(进程)；

glob2loc_well_map，全局井的编号到子区域内局部井的编号的映射(数组的值从0开始)；

loc2glob_well_map，井的局部编号到全局编号的映射(数组的值从0开始)；

local_nwell，本进程的子区域的井的个数；

Nwell，jnwt，MSTCOM.TYT和WELCOM.TYT文件中井的个数的信息。

输入输出：

Mwiloc，每口井的i(x)方向坐标；

Mwjloc，每口井的j(y)方向坐标；

Mipf，每口井，每个射孔层的i(x)方向坐标；

Mjpf，每口井，每个射孔层的j(y)方向坐标；

Mkpf，每口井，每个射孔层的k(z)方向坐标；

Mperf，每口井，每个射孔层的有效网格编号(这里是去掉了死网格点的编号)

以上这些数据，输入为全局的井和网格，输出为局部的井和网格。bhp...等记录每个井的数据的数组，maback....等记录每个井的信息的数组。

(2)动态内存的分配

整个的模拟系统要运行在不同的油田数据模型中，整个区域的网格大小是不确定的，所以在不同模型上，应用schwarz方法时发送/接受的数据个数是不同的，而数组的声明又必须使用常量，这就给缓冲区数组的声明带来了困难，只使用常量声明无法确定缓冲区数组的真实大小。在软件系统的外部文件PARAMS.TYT中，定义了各种同区域大小和分割方式相关的常量，其中每个区域的最大网格数目定义为lbmax＝250,000，从全局区域分割出子区域的数目最多为max_proc＝128。由于要发送/接收数据的进程数目不会超过最大进程数max_proc，同时每个进程要发送的数据数目不会超过最大网格数目的四倍(因为每个网格有四个解要发送)，这样在声明时，可以将根进程缓冲区数组的大小定义成

(lbmax×4)×max_proc

将非根进程的接收/发送数据缓冲区的大小定义成：lbmax*4但这样声明对内存空间的浪费是相当严重的，因为很少有区域划分能将整个区域划分成max_proc＝128个子区域，每个子区域的网格数目也很少能达到

lbmax＝250,000的规模，以目前测试过的最大数据模型为例，该模型有百万规模的网格数，大小为210×210×24，在x方向上分成210个网格，Y方向上分成210个网格，Z方向上分成24个网格，对于给定模型，一般将其等分成49个子区域进行并行计算(子区域如果划分得过多，会造成通信开销增大，降低并行效率)，这样，每个子区域的内部网格(非重叠部分)数目为 $\frac{210\×210\×24}{49}=21600,$ 设重叠部分的网格数目是非重叠部分的网格数目的30％(一般在区域划分时都不会将重叠区域的层数选得过多)，则每个子区域的重叠区网格数目为21600×30％＝6480，那么每个子区域发送/接受数据缓冲区的实际有效数据个数为：6480×4＝25920，而实际上每个子区域的数据缓冲区数组在声明时分配的大小为：lb max×4＝250000×4＝1000000，数据空间的利用率为： $\frac{25920}{1000000}=2.592\%,$ 就是说每个子区域的数据缓冲区中有97.408％的空间被浪费了。根进程的数据缓冲区也类似，根进程数据缓冲区的实际有效数据个数为：(6480×4)×49＝1270080，而在声明时分配给根进程数据缓冲区数组的大小为：

(lbmax×4)×max_proc＝(250000×4)×128＝128000000

数据空间的利用率为： $\frac{1270080}{128000000}=0.99225\%,$ 有接近99％的根进程数据缓冲区空间被浪费。

在高性能计算中，对内存空间的浪费会降低程序执行的效率，所以将根进程数据缓冲区和非根进程数据缓冲区都定义成了动态分配内存，需要在声明的时候将缓冲区BUF_DATA_ALL，RSFT，RS都声明成指针的形式，然后再在使用之前调用FORTRAN90中的ALLOCATE函数为其分配空间。由于在主进程中，每个进程要发送的数据个数已经保存在了数组DATA_SEND_OVL_NUM(NP)中，将数组DATA_SEND_OVL_NUM中的各个元素累加起来，就是主进程要接收/发送的数据个数，也就是要为BUF_DATA_ALL分配的内存空间大小。同理，每个子区域位于重叠部分的网格数目MY_OVERLAP_NUM乘以4就是为该子区域对应进程的发送/接收缓冲区分配的内存空间大小。

同样的动态内存分配技术也可以应用在其它存储空间中，比如前面提到的重叠区网格编号索引OVL2GLOB_INDX，根进程中保存的每个从进程要发送的数据个数、位置偏移等，它们的实际有效大小也是同全局区域的大小和进程数的多少密切相关的，动态内存分配也可以节省这些存储区域的存储空间。

3.构建线性方程组

(1)渗流方程组系数矩阵的形成过程

对三相三组分黑油模型，程序调用填充方程左侧系数矩阵的子程序LHS3项和右端残余项子程序RHS。这是形成方程的重要部分。经过详细分析，形成系数矩阵的子程序主要涉及到四种情况下方程矩阵的形成：

①油饱和度大于零，并且气饱和度大于零的情况；

②油饱和度大于零，并且气饱和度等于零的情况，称为“死油”；

③油饱和度等于零，并且气饱和度大于零的情况，称为“死气”；

④油饱和度等于零，并且气饱和度等于零的情况，这种情况是“全水”。

对于上述四种情况，方程中的部分项要采用不同的计算公式，因而要分别对待。图4显示了形成系数矩阵和右端残余项的子程序结构。在整个方程的形成过程中，只有此处涉及到了网格间的物理量计算，因此，并行程序的子区域边界层交换数据也发生在此处。

(2)子区域边界条件更新并行化实现

空间区域分解之后，每个子区域的边界条件要重新构造，利用边界条件形成方程组的过程要修改，每次迭代之后要更新区域边界上的数据。这些问题是实现牛顿迭代并行化的关键，需要分析利用边界条件生成方程组的过程，设计适合算法的边界数据的存储结构，确定边界条件更新的时间和方法，因此这是对牛顿迭代并行化的难点所在。

在油藏数值模拟串行程序中，模拟区域只有一个整个的边界，串行程序使用这个边界信息构造全局的方程组。在区域分解之后，每个子区域要有自己新的边界数据，并且要使用这个新的边界数据构造本子区域的方程组，因此引入了边界部分的数据结构，每一个和边界部分相关的数组都被分割到对应数据结构中。对位于整个区域中间的子区域来说，这些边界数据对是它周围四个相邻区域中的网格点上的信息，对于位于整个区域边上的子区域来说，这些边界数据包括它相邻区域中的网格点上的信息和原来整个区域的边界上的信息。

①子区域边界层流量交换过程

具体流量交换过程见图5。以区域内部x方向流量交换为例，x方向上从第二个网格开始循环，直到该行的最后一个网格(这里是9)。x＝2时，当前网格标志为P，当前网格前面邻近网格标志为M，该循环处理网格M和P之间的流量交换。如果M或者P是死网格，网格间没有流量交换，则结束该次循环。否则，存在流量交换。网格M和P的活动网格序号分别为MM1，JMM1；排序后的网格序号是IM1，IM。按照前面推导的相关公式，形成网格M和P的流量交换，进而生成方程系数矩阵C和V；x依次递增，直到x＝9时，结束该行的循环；y方向上递增，重复上述过程，直到y等于Y方向上的最大值(这里是9)。

完成内部X和Y方向的流量交换后，开始边界层的流量交换，以x方向上的边界流量交换为例，见图6。在y＝1行左边界上，以左侧邻近子区域的边界层网格为M，以本子区域边界层网格为P，按照内部网格的公式计算流量交换。需要注意的是，网格M的所有物理量，都需要从邻近网格传递过来，是存在于为并行处理而构造的数据结构中的；在y＝1行右边界上，以本子区域边界层网格为M，以左侧邻近子区域的边界层网格为P，按照内部网格的公式计算流量交换。需要注意的是，与左边界相反，从邻近网格传递来的数据在这里是两层网格中的右侧网格。因此，在公式中需要调换物理量的位置，要十分小心；在Y方向上依次递增，直至y等于Y方向上最大值，处理每行左、右边界上的流量交换。如果边界网格是死网格，则不需要处理。

上、下边界层的流量交换过程与此相同，只是要注意下标的处理。经过上下左右边界上的流量交换后(如图7)，子区域方程系数矩阵的某些位置上出现由于边界交换而产生的数值。这些位置不会出现在主对角线C上，因此，C的数据不会变化。新增的系数出现在非对角线元素V的某些位置上。

②子区域边界流量交换对方程系数矩阵的影响

在形成方程系数矩阵过程中，分别在x，y，z三个方向上计算网格间流量交换。下面以x方向上的流量交换为例，讨论对系数矩阵的影响。

a.左边界流量交换对方程组系数矩阵影响

设l₁₁＝1，(x方向)，LD＝1(第一条非对角线元素V₁)。方程组系数矩阵非零对角线元素的排列见图8。令三维空间的当前网格取kk＝1，jj＝1，ii＝2(记做网格2)，该网格与前面的网格1之间作流量交流，网格0代表子区域外部相邻区域的边界层网格，如图9所示。

在计算相邻两个网格的流量交换时，设定左侧网格1除去死网格后的活动网格序号是MM1，在方程组中的排序编号是IM1；右侧网格2除去死网格后的活动网格序号是JMM1，在方程组中的排序编号是M1。因此，在计算网格0和网格1流量交换时，它们的活动网格序号表示为MM1′和JMM1′，方程系数矩阵中的编号分别表示为IM1′和IM′，如图2-9所示。

两个网格1、2形成的系数矩阵F是如下形式，其中每个矩阵元素都是一个3×3的矩阵。

$F=\left[\begin{array}{cc}C(\mathrm{IM}\mathrm{1,1}~9)& V(\mathrm{IM}\mathrm{1,1}~9,\mathrm{LD}+1)\\ V(\mathrm{IM},1~9,\mathrm{LD})& C(\mathrm{IM},1~9)\end{array}\right]---\left(1\right)$

由公式：

$V(\mathrm{IM},2,\mathrm{LD})=J\×(p_o^{\mathrm{MM}1}-p_o^{\mathrm{JMM}1}-\mathrm{UZ}\×\frac{({\mathrm{\γ}}_o^{\mathrm{MM}1}+{\mathrm{\γ}}_0^{\mathrm{JMM}1})}{2})---\left(2\right)$

$V(\mathrm{IM}1,2,\mathrm{LD}+1)=J\×(p_o^{\mathrm{JMM}1}-p_o^{\mathrm{MM}1}-\mathrm{UZ}\×\frac{({\mathrm{\γ}}_o^{\mathrm{MM}1}+{\mathrm{\γ}}_0^{\mathrm{JMM}1})}{2})---\left(3\right)$

可知，第二行第一列元素V(IM，1～9，LD)是由于网格1的存在而对网格2产生的影响。第一行第二列元素V(IM1，1～9，LD+1)是由于网格2的存在而对网格1的存在而对网格1产生的影响。这里的IM和IM1都是经过映射的排序序号，UZ是网格的相对地面的深度。

因此，得到结论如下：(左边)网格IM1对(右边)网格IM的影响，反映在系数矩阵与(右边)网格IM同行、与(左边)网格IM1同列的非对角线元素V(IM，1～9，LD)上；(右边)网格IM对(左边)网格IM1的影响，反映在系数矩阵与(左边)网格IM1同行、与(右边)网格IM同列的非对角线元素V(IM1，1～9，LD+1)上。

在计算子区域内部相邻两个网格间流量交换时，IM是从网格2开始，到x方向最大网格N_x结束；IM1从网格1开始，到N_x-1结束。如果考虑空间分解算法，在子区域人工边界左侧引入相邻区域边界网格0后，由上面的分析得知，会在系数矩阵F外面出现因为边界网格对1号网格影响而产生的元素V(IM′，1～9，LD)，在与网格1同行且与网格0同列的非对角线位置上，如图10所示。将V(IM′，1～9，LD)乘上边界网格IM1′对应的未知变量，并从右端项中减去，就可以恢复原来形状的矩阵，并且该矩阵是考虑了边界影响的。

b.右边界流量交换对方程组系数矩阵影响

同理，在右边界处，我们要分析的是边界网格0对网格2的影响(如图2-11)，即V(IM1′，1，LD+1)。

注意：在边界处理函数中，子区域内部最右边的网格为IM1′(＝2)，边界层网格为IM′(＝0)，因此，右侧边界流量引起的V的变化，反映在数组元素V(IM1′，1～9，LD+1)中，而不是V(IM′，1～9，LD)，如图2-12所示。

4.线性方程组并行求解

(1)加法Schwarz方法

由于经典Schwarz需要先求出一个区域上的解，再以该解作为边界条件求另一个区域上的解，这种串行方法对于解空间较大的情况下效率很低，为此人们提出了并行求解的加法Schwarz方法。

对于椭圆型Dirichlet问题：

$\left\{\begin{array}{cc}\mathrm{Lu}=f& \left(\mathrm{\Ω}\right)\\ u=g& (\∂\mathrm{\Ω})\end{array}\right.---\left(4\right)$

将区域Ω分解为m个子区域，即 $\underset{i=1}{\overset{m}{\∪}}{\mathrm{\Ω}}_i=\mathrm{\Ω},$ 且每个子区域至少和另一个子区域有交集(同该子区域有交集的子区域的多少可以根据问题的复杂程度任意选取)

加法Schward的执行过程描述如下：

(1)取初值，选择初始近似 $u_0\∈H_g^1\left(\mathrm{\Ω}\right),$ n：＝0

(2)求解，并行计算子区域上的边值问题：

$\left\{\begin{array}{c}{\mathrm{Lu}}_i^{n+1}=f,\left({\mathrm{\Ω}}_i\right)\\ u_i^{n+1}=u^n,\left({\∂\mathrm{\Ω}}_i\right),i=1,...m\end{array}\right.---\left(5\right)$

(3)解更新，将u_i ⁿ⁺¹的定义延拓到整个Ω空间，定义：

$\left\{\begin{array}{cc}{\stackrel{2}{u}}_i^{n+1}=u_i^{n+1}& x\∈{\mathrm{\Ω}}_i\\ u^n& x\∈\mathrm{\Ω}\backslash {\mathrm{\Ω}}_i\end{array}\right.---\left(6\right)$

最后取算术平均，本次迭代(第n+1次)得到的解为：

$u^{n+1}=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\stackrel{2}{u}}_i^{n+1}---\left(7\right)$

(2)加速收敛策略

采用区域分块并行方法后，在人为边界处网格之间的影响会被割裂。通过引入并行边界条件后，能够正确形成人为边界附近网格的系数项。但是，在各个进程求解过程中，右端项减去的是上一时刻的边界系数项与未知数的乘积，引入误差，导致用ESPIDO方法作为线性求解器的核心算法而解出的局部区域解的误差较大。虽然严格的误差容限能够提高精度，但是牛顿迭代次数会大幅增加。这是目前面临的最大困难。

从上述分析中可以看出，要使并行算法的迭代次数降低到与串行相当的水平，可以采用两条途径：首先，研究并行程序求解的方程组的校正方法，使之在最大程度上接近真实方程曲面；其次，引入误差方程快速更新技术，使收敛速度加快。

①动态解曲面预估

方程的系数矩阵和右端项是决定方程形态的重要部分。并行程序系数矩阵与真实方程系数矩阵越接近，非线性迭代次数就会越接近串行方法。在初始求解阶段，并行和串行程序形成的方程系数矩阵比较接近，并行系数矩阵缺失的元素占整个矩阵的比例很小，因此影响也比较弱。因此，并行和串行程序的迭代次数相对比较接近。随着计算时间的增加，并行程序系数矩阵缺失部分的比重逐渐增加，对方程形态的影响也增大，误差方程曲面距离真实曲面越来越远，因而导致非线性迭代次数大幅增加。

为了减轻缺失元素的影响，采用了误差方程校正的方法，使并行程序的方程曲面在函数空间中尽量接近真实方程组，以期降低迭代次数。由于系数矩阵和右端项在每次牛顿迭代中都要重新形成，因此，这个形成过程是校正误差方程的关键。一般认为，整个空间压力场是连续的(除了个别特殊网格附近)，而且随时间变化是连续的，甚至是C1，C2连续的。借鉴有限单元法中形函数的概念，可以构造一个压力场时间形函数H(p，t)。有限单元法在空间网格内部用形函数近似逼近真实解，在网格边界上满足一定的连续条件，如C0，C1连续。该方法在时间网格(可能包含几个时间步长)内用时间形函数近似逼近真实解，并在时间网格边界上满足连续条件。基于非线性外推方法，通过对已知压力点预测未来时间网格的压力点数值，能够将误差系数矩阵非常快速的校正到真实函数曲面，大大降低并行程序牛顿迭代的次数，具有明显效果。

最简单的初始解预估方法是线性外推算法。根据现在时刻压力值p^t ₀和过去时刻的压力值p^t _-1，计算下一时刻的压力值p^t ₁，计算公式为：

$p^{t_1}=p^{t_0}+\frac{t_1-t_0}{t_0-t_{-1}}\·(p^{t_0}-p^{t_{-1}})---\left(8\right)$

基于抛物线方程的二次外推预估算法，根据当前时刻压力值p^t ₀和过去两个时刻的压力值p^t _-2，p^t _-1构造二次公式，计算下一时刻的压力值p^t ₁计算公式为：

p^t ₁＝a·(p^t ₁)²+b·p^t ₁+c               (9)

其中系数a，b，c由下式确定：

$\left\{\begin{array}{c}a\\ b\\ c\end{array}\right\}={\left[\begin{array}{ccc}{t}_{-2}^2& t_{-2}& 1\\ t_{-1}^2& t_{-1}& 1\\ t_0^2& t_0& 1\end{array}\right]}^{-1}\·\left\{\begin{array}{c}{p}_{-2}\\ p_{-1}\\ p_0\end{array}\right\}---\left(10\right)$

基于有理式的非线性外推预估算法，根据当前时刻压力值p^t ₀和过去n-1个时刻的压力值构造一个有理式函数，计算下一时刻的压力值p^t ₁。

②并行解的邻域校正方法

由于并行方程系数矩阵缺失元素，迭代中修正方程得到的解曲面不能快速接近真解曲面，导致算法收敛性大幅降低。

为了有效地将区域边界外网格的影响有效计入各区域的迭代过程，我们提出了压力邻域校正加速策略。该方法的核心思想如图13所示。

图中G_b为区域内的边界网格，G_a是与G_b相邻的其它子区域内网格；p_Gb ^k和p_Ga ^k是两个网格的第k次迭代的压力解，而p_Gb和p_Ga是两个网格的压力真解。由于p_Gb ^k和p_Ga ^k基于流量Q＝J(p_Ga-p_Gb)建立的线性相关性，可以用α×p_Ga ^k作为G_a对G_b的影响在G_b所在压力子空间中的量化估计，修正p_Gb ^k，使其更有效地接近真解p_Gb，提高算法收敛性。

$p_{G_b}^k=p_{G_b}^k+\mathrm{\α}\×p_{G_0}^k---\left(11\right)$

该方法的难点在于修正系数α的选择，因为修正系数的作用随计算问题不同而往往有所变化。我们针对不同的加速系数α，在多个模型上根据实际测试发现，a取值从0到0.6，迭代次数呈下降趋势，但是当a≥0.7时，会发生过度修正，使得迭代收敛速度更慢。

5.结果收集整理输出

油藏数值模拟问题并行计算输出特点：

(1)并行化是在每个时间步之内进行的；

(2)在每次进入下一时间步之前，各进程数据相同；

(3)大部分输出发生在每个时间步结束之后。

利用一个0号进程为根进程负责输出，输出之前增加了收集输出数据的函数。数据收集工作由gatherout函数完成。其定义如下：

SUBROUTINE GATHEROUT(PF，SF，RSF，RF，SRSF，TAUF，SIGF，

*                    BWF，BOF，BGF，G，SWCAR，

*                    SGCAR，XCO2F，SWA，

*                    CVEC，RCUMP，RCUMI，FWP，FTP，THP，

*                    MTHPT，MWSTAT)

函数参数中列出的数据结构是可能需要收集的，不是一定收集的，对非根进程来说，这些数据结构是输入的；对根进程来说，这些数据结构在gatherout之后，可能会变成全局的数据结构。具体需要收集哪些数据，由gatherout调用的setgatherlist函数设置。

Gatherout分为以下几部分实现(图14)：

首先调用setgatherlist设置需要收集的数组标志。

非根进程调用packdata将需要发送的数据压缩入消息缓冲区。

调用gatherdata进行进程间通信。

根进程调用unpackdata从消息缓冲区中还原数据到全局数组的对应位置。

Claims (1)

1.一种油藏黑油模型数值模拟系统，包括前处理子系统、串并行模拟器和后处理子系统，其特征在于：根据前期油藏描述的成果，由前处理子系统建立油藏黑油模型数值模拟模型，然后利用串并行模拟器，进行模拟运算，再利用后处理子系统分析模拟的结果，设计优化油藏开发调整方案，最后形成文字报告和汇报多媒体，提交最佳设计方案；

所述前处理子系统用于建立数值模拟模型，包括项目管理、方案管理、数据管理、模型定义、网格建立、饱和度函数、流体性质、分区设定、水体描述、动态模型建立、作业控制和输出控制12个模块，项目管理实现项目相关信息的统一组织与管理，具体包括项目建立、项目打开、保存、另存、项目删除、退出6项功能；方案管理实现项目不同研究方案的管理，具体包括新建方案、输入方案描述信息，拷贝方案、保存方案、删除方案5项功能；数据管理实现项目涉及数据的统一管理，具体包括数据输入、数据检查、数据编辑、数据显示、数据分析、数据保存、数据输出多项功能；模型定义，定义所建模型的基本信息，网格建立用于生成表征油藏地质信息的网格模型，具体包括数据加载、编辑、场图插值、网格设计、网格编辑、网格优化、网格赋值、三维图形显示多项功能；饱和度函数主要是建立岩石流体的饱和度函数及毛管压力函数；具体包括相对渗透率曲线和毛管压力曲线导入、曲线显示、曲线编辑以及曲线处理和数据输出；流体性质主要是处理分析油气水的实验分析数据，建立描述流体性质的流体模型；分区设定主要是设定不同饱和度函数、流体性质以及储量统计对应的分区；水体描述用于表征油藏的水体分布及大小；动态模型建立主要是通过油藏开发历史数据的分析与处理，建立反映油藏开发过程的动态模型；作业控制用于设置模型运算的迭代误差控制、时间步长控制；输出控制用于设置输出的内容和频率；

其后处理子系统用于分析模拟结果，设计优化油藏开发方案，包括自定义参数场计算、平面视图绘制、曲线图绘制、三维立体视图绘制、油藏纵向切片图绘制、注采流线图绘制、四位一体的油藏多指标综合图绘制、油藏动静态参数定量动态描述8个模块，自定义参数场计算利用已有的常规动静态参数场进行组合计算，生成新的有地质意义的数据场，用户可自定义最多10个计算参数场；平面视图绘制具有两大特色；一是图件绘制齐全且有创新：不仅可以绘制彩色分布图、填充等值线图、等值线图、含水柱状图、累油水饼图，还可以绘制叠加等值线图；二是图件绘制参数多样：包括油藏静态数据、动态数据、原始数据及自定义场数据四类；曲线图绘制可以绘制区块及油水井的累积、瞬时曲线图，分析不同时期区块及油水井的生产状况，设置了灵活的绘图工具：一是绘图的小层、井号及参数可以任意选择；二是绘图参数的X、Y坐标可以任意选择，曲线的条数、绘图的时间段可以任意选择；三是曲线图中的字体大小、字型、曲线的线型、颜色等均可自行选择，方便了用户；三维立体视图绘制包括三维立体彩色分布图、三维立体等值线图两种图件的绘制；通过三维立体图的绘制，可以直观地反映油藏在平面和纵向上的变化特征；油藏纵向切片图绘制可以自动绘制静态数据、动态数据、原始数据及自定义场数据四类参数的油藏切片图；通过油藏切片图的绘制，可以直观地反映油藏在纵向上的变化特征，搞清每口油井钻遇的油层数、有效厚度、含油饱和度、剩余储量等基本数据，同时可以直观地描述剖面线上油藏的非均质变化特征；注采流线图绘制可以根据当前时间压力场绘制注采流线分布图，显示流体流动的方向，分析注采对应状况；根据注水流线的分布图，可以进行注采井网完善状况的判断与分析，确定哪些区域是优势流场，哪些区域是非优势流场，确定剩余油潜力范围和挖潜方向；四位一体的油藏多指标综合图绘制可以在同一微机屏幕上同时进行4种油藏地质图件的显示和对比分析，从而多角度、多方面的对油藏进行综合分析；油藏动静态参数定量动态描述把油藏储层参数和剩余油描述指标有机地结合起来，在微机上形成细化的网格数据体，全方位、多指标的对油藏动静态参数进行定量描述，并通过时间坐标的设置，实现油藏二维、三维图件的动态显示，通过不断补充数据实现对油藏的跟踪分析和数据化管理；

所述串并行模拟器以数值模拟模型作为输入文件，进行模拟运算，串并行模拟器的实现包含了以下关键技术：可重叠的空间子区域划分方案、负载平衡优化设计、动态分配内存、网格重排方法、区域分解中分割数据的快速索引、加法Schwarz方法、并行解的邻域校正方法、曲面预估方法、线性方程组求解参数及方法优化、关键技术包括数据压缩通讯、压缩存储技术；这些关键技术的应用提高了模拟运算的速度，解决大规模精细数值模拟研究问题。

Images