CN116882218B - 一种油藏数值模拟方法、装置、计算机设备及存储介质 - Google Patents

Abstract

本发明涉及油藏数值模拟技术领域，公开了一种油藏数值模拟方法、装置、计算机设备及存储介质，方法包括：基于井压力创建用于耦合求解油藏的各相压力和各相饱和度的质量守恒方程；创建关于各相压力、各相饱和度和井压力的井方程；初始化包括当前时间步的井压力的初始解向量；利用初始解向量隐式求解质量守恒方程得到下一时间步的各相压力和各相饱和度；利用下一时间步的各相压力和饱和度显式求解井方程得到下一时间步的井压力；将下一时间步的井压力作为当前时间步的井压力，返回隐式求解质量守恒方程的步骤迭代到预设数量时间步得到目标各相压力和目标各相饱和度。本发明同时提高了油藏数值模拟的数值稳定性、时间步长宽松性和负载均衡性。

Description

一种油藏数值模拟方法、装置、计算机设备及存储介质

技术领域

本发明涉及油藏数值模拟技术领域，具体涉及一种油藏数值模拟方法、装置、计算机设备及存储介质。

背景技术

油藏数值模拟是指利用计算机求解油藏数学模型，模拟地下油水流动，给出某时刻油水分布，预测油藏动态的技术。通俗来讲，油藏数值模拟用于求解石油在油藏体中各个位置的含量以及油藏打井处的产量，油产量是石油生产中最关键部分，这决定了石油产出的经济效益。多相流模拟问题是油藏模拟中普遍的流体问题，在油田开发中，地下的油藏通常包含多个相，比如液相(油，oil)、气相(天然气，gas)和水相(水，water)，这些相之间存在相互作用，对于预测和优化油田的开采至关重要。多相流模拟的目的就是通过数学模型和计算方法，模拟这些相之间的运动、相互作用和质量传递过程，以便更好地了解和预测油藏中的流体行为。而油藏中各个位置的石油含量可以通过相的饱和度（饱和度的定义是石油体积占油藏孔隙体积的百分比）直观表示出来，油藏打井处的产量通过相的压力与饱和度两个量代入石油产出量的公式计算得到。

目前，常用的油藏数值模拟方法主要包括全耦合方法和解耦方法两大类，其中全耦合方法将各相的压力方程、饱和度方程以及井方程耦合在一起求解，从而计算得到相的压力和饱和度，虽然全耦合计算稳定性较好，但是井方程计算出的井压力变量个数与油藏网格块对应的变量个数不一致，导致网格内的参数不一致，在通过多个处理器进行并行计算时，网格不能实现均匀分配计算负载给每个处理器，导致耦合方程系统在求解过程中存在负载不均衡的问题。常用的解耦方法是将各相的压力和井方程耦合在一起作为压力方程，剩余其他守恒公式作为饱和度方程，饱和度方程与压力方程采用隐压显饱法或者隐压隐饱法交替求解压力和饱和度。一方面，井方程计算出的变量个数与油藏网格块对应的变量个数依然存在不一致的问题，导致并行计算时负载不均衡；另一方面，解耦方法对求解的时间步长有一定影响，尤其是隐压显饱法受到CFL收敛条件的影响，令时间步长很短，严重影响计算效率；此外，解耦方法相比全耦合方法计算的数值稳定性也较差。所以，亟需一种能够兼顾数值稳定性、时间步长宽松性和负载均衡性的油藏数值模拟新方法。

发明内容

有鉴于此，本发明提供了一种油藏数值模拟方法、装置、计算机设备及存储介质，以解决目前的油藏数值模拟方法不能兼顾数值稳定性、时间步长宽松性和负载均衡性的问题。

第一方面，本发明提供了一种油藏数值模拟方法，所述方法包括：基于井压力创建用于耦合求解油藏的各相压力和各相饱和度的质量守恒方程，所述井压力是注入井的井底压力；创建关于各相压力、各相饱和度和井压力的井方程；初始化包括当前时间步的井压力的初始解向量；利用所述初始解向量隐式求解所述质量守恒方程，得到下一时间步的各相压力和各相饱和度；利用所述下一时间步的各相压力和各相饱和度显式求解所述井方程，得到下一时间步的井压力；将所述下一时间步的井压力作为所述当前时间步的井压力，返回所述隐式求解所述质量守恒方程的步骤，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度。

在一种可选地实施方式中，所述质量守恒方程如下

其中，为孔隙度，、l和v分别表示水相、液相和气相，表示三相之一，、和分别为水相、液相和气相条件下的密度，、和分别为各相的流体速度，、和分别为各相饱和度，、和分别为各相压力，、和分别为各相的注入速率，由井压力计算得到，和分别表示液相和气相中第i个组分的质量分数，和分别表示液相和气相中第i个组分的摩尔分数，和表示液相和气相中第i个组分的逸度，和表示液相和气相逸度的系数值，表示的是计算逸度系数的给定参数，m为各相所含组分的总数量，其中为组分j的摩尔质量，为组分i的摩尔质量，和分别表示液相和气相中第j个组分的摩尔分数，为油藏的绝对渗透率，为重力加速度，、和为相流体的相对渗透率、粘性系数和移动率，为液相压力与相之间的毛细管压力，为压缩系数，和分别表示分子间的吸引力无量纲量和排斥力无量纲量，是液相第k个组分的摩尔分数，和分别为组分i的简化温度和压力值，和是中间变量，、和为计算中间变量和的经验参数，为组分i和k间的二元相互作用系数。

在一种可选地实施方式中，所述井方程如下

其中，是相注入液的注入量，是相流体的注入速率，是第口井的射孔总数，是第口井的第个射孔的段长、和分别是井位置处相流体的绝对和相对渗透率值，和分别为井位置处相流体的粘性系数和密度，表示第口井的井筒半径，为第口井第个射孔所在网格块的排水半径，和分别表示井处的井压力深度和当前射孔位置深度，为井位置处的井压力,为井处相流体的压力，是重力加速度。

在一种可选地实施方式中，所述初始化包括当前时间步的井压力的初始解向量，包括：初始化当前时间步的井压力、当前时间步的各相压力、当前时间步的各相饱和度、当前时间步的液相中各组分的摩尔分数、当前时间步的气相中各组分的摩尔分数，得到所述初始解向量 (,,,,)。

在一种可选地实施方式中，所述利用所述初始解向量隐式求解所述质量守恒方程，得到下一时间步的各相压力和各相饱和度，包括：固定各相注入液的注入量，将所述初始解向量代入所述井方程，显式计算当前时间步的各相注入速率；将所述当前时间步的各相注入速率代入所述质量守恒方程，进行隐式求解，得到下一时间步的各相压力、各相饱和度、液相中各组分的摩尔分数和气相中各组分的摩尔分数。

在一种可选地实施方式中，所述利用所述下一时间步的各相压力和各相饱和度显式求解所述井方程，包括：将所述下一时间步的各相压力、各相饱和度、液相中各组分的摩尔分数和气相中各组分的摩尔分数代入所述井方程，进行显式求解，得到下一时间步的井压力。

在一种可选地实施方式中，通过用牛顿-克雷洛夫-施瓦兹算法对所述质量守恒方程进行隐式求解。

第二方面，本发明提供了一种油藏数值模拟装置，所述装置包括：质量守恒方程创建模块，用于基于井压力创建用于耦合求解油藏的各相压力和各相饱和度的质量守恒方程，所述井压力是注入井的井底压力；井方程创建模块，用于创建关于各相压力、各相饱和度和井压力的井方程；初始解模块，用于初始化包括当前时间步的井压力的初始解向量；隐式求解模块，用于利用所述初始解向量隐式求解所述质量守恒方程，得到下一时间步的各相压力和各相饱和度；显式求解模块，用于利用所述下一时间步的各相压力和各相饱和度显式求解所述井方程，得到下一时间步的井压力；反馈迭代求解模块，用于将所述下一时间步的井压力作为所述当前时间步的井压力，返回所述隐式求解所述质量守恒方程的步骤，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度。

第三方面，本发明提供了一种计算机设备，包括：存储器和处理器，存储器和处理器之间互相通信连接，存储器中存储有计算机指令，处理器通过执行计算机指令，从而执行上述第一方面或其对应的任一实施方式的方法。

第四方面，本发明提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机指令，计算机指令用于使计算机执行上述第一方面或其对应的任一实施方式的方法。

本发明提供的技术方案，具有如下优点：

本发明提供的技术方案将各相的压力参数与饱和度参数耦合到同一组质量守恒方程中，并单独创建井方程。两个方程交叉迭代计算，通过井方程在上一时间步计算得到的井压力，更新质量守恒方程下一时间步输出的各相压力和各相饱和度，然后通过更新的各相压力和各相饱和度返回代入井方程计算下一时间步的井压力，反复迭代直至求解最终时间点的目标各相压力和目标各相饱和度。本实施例通过隐式求解算法求解全耦合的质量守恒方程，令相的压力以及饱和度计算变为全隐式计算，从而继承隐式求解时间步更宽松的特点，得以使用更长的时间步，加快求解效率。同时，压力与饱和度的计算采用全隐式求解使数值求解的结果稳定性更高，减少振荡的问题。之后，当质量守恒方程计算出下一时间步的各相压力和各相饱和度时，单独对井方程进行显式求解，虽然注入井的数量与油藏模型网格的数量并不一一对应，但是单独求解井方程得到的下一时间步的井压力参数并不分配到油藏网格中，被剔除在外。本实施例提供的技术方案相比井方程与相压力方程耦合在一起的计算方法，当使用多个处理器执行并行计算时，不会出现网格计算负载分配不均匀的问题，提高了油藏数值模拟的负载均衡性。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是相关技术中的隐压显饱求解法的流程示意图；

图2是相关技术中的隐压隐饱求解法的流程示意图；

图3是相关技术中多处理器分配网格进行并行计算的结构示意图；

图4是根据本发明实施例的一种油藏数值模拟方法的流程示意图；

图5是相关技术中的三相状态转换的结构示意图；

图6是根据本发明实施例的另一种油藏数值模拟方法的流程示意图；

图7是根据本发明实施例的又一种油藏数值模拟方法的流程示意图；

图8是根据本发明实施例的一种油藏数值模拟装置的结构示意图；

图9是本发明实施例的计算机设备的硬件结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

目前，在油藏数值模拟的相关技术中，为了求解各相的压力和各相的饱和度，常用解耦求解方法和全耦合求解法。其中解耦求解方法有隐压显饱法和隐压隐饱法，其中隐压显饱法将多相流模拟问题分解为两个子问题，隐式求解压力方程和显式求解饱和度方程。对于油藏网格模型，先创建关于该模型的水相（water）、气相（gas）和液相（oil）的数学模型，然后将数学模型分别转换为压力方程和饱和度方程，之后在每个时间步（在数值模拟中，时间步是指将模拟的时间域划分为离散的时间间隔，模拟过程中每个时间间隔称为一个时间步，在每一个时间步内，求解该时刻处的方程或方程组，得到收敛结果后，再继续进行下一时间步的循环），如图1所示，首先进行隐式求解压力方程，以确定网格单元中的压力场，随后利用网格单元已知的压力场进行显式求解饱和方程，以确定各个相的饱和度，并通过交替循环迭代直到收敛，即先通过先隐式求解时间步n到n+1的各项压力，然后利用求解的压力再显示求解时间步n到n+1的饱和度S。压力与饱和度方程的构造，主要是为了获得一个稳定的迭代格式，因为压力量在计算过程中是剧烈变化的，若采用显式求解，迭代很容易发散，因此这部分必须隐式求解，相反饱和度是一个0到1的分数，变化剧烈度要小，迭代的时候可以采用显式或隐式求解。由于在多相流系统中，最关键的是要得到油藏体中的压力和饱和度量，隐压显饱求解法中的压力方程和饱和度方程是按顺序求解的，即首先通过压力方程隐式求解得到压力分布，然后根据得到的压力分布显式求解饱和度方程得到饱和度分布，这两个方程的求解是交替进行的，这种方法的优点是计算机内存占用较小、计算工作量较少、方法简单易行。这种方法的缺点是在处理非均质油藏时（孔隙等油藏参数随空间位置变化发生变化的油藏），对时间步长有严格的限制，时间步长取的稍大仍会出现振荡现象，数值稳定性较差，即压力和饱和度出现负值等情况，以致模拟无法正常进行，这是由于对饱和度的显式处理导致的，无法确保其在复杂非均质油藏中模拟的稳定性和可靠性。

如图2所示，隐压隐饱求解法是在每个时间步内，首先对各个相的压力方程进行隐式求解，待压力方程求解完成后，利用已知的压力场来隐式求解饱和方程，以确定各个相的饱和度。即先通过先隐式求解时间步n到n+1的各项压力，然后利用求解的压力再隐式求解时间步n到n+1的饱和度S。该方法与隐压显饱法的区别在于先隐式求解压力方程，后隐式求解饱和度方程。这种方法的计算复杂度稍高，但是在求解过程中具有更高的灵活性和准确性，在数值稳定性方面优于隐压显饱法方法，并且对时间步的要求更宽松。例如：假设用户想模拟油藏体内10年后的石油分布和产量，很难将时间步设定为10年，通常要将10年划分为不同小的时刻，比如0.5 年，算完一个0.5年后，再用这个时刻的变量更新下一个0.5年的饱和度和压力值，直至到10年时刻，这里的0.5年（每个时间间隔）就称为时间步长。在数学上已有成熟的理论说明，显式求解方法在时间步长上的限制比隐式更为严格，因为显式方法的稳定性受CFL条件数影响，要满足CFL收敛条件，就不得不压缩时间步的长度，否则就会令计算结果振荡，而隐式求解不受限。比如，隐式方法在计算时，时间步长设为0.5年，显式方法在计算时，要想保持整体迭代稳定，它的时间长就只能设置为更小，比如0.01年，如此在迭代求解过程中才能稳定，否则计算会不收敛到求解结果。另外，隐压隐饱法的缺点在于，在处理非均质油藏时可以放宽稳定性限制，但相对来说仍然有一定的稳定性限制，该方法对压力和饱和度相关方程分别隐式求解，虽然在数值稳定性方面稍好，但由于压力和饱和度的分离，仍可能对时间步长施加严格的限制。因为关于井压力的计算会耦合到压力方程中，这两类方法还存在并行处理时负载不均衡的缺点。

关于负载不均衡的缺点，在全耦合求解法中尤其明显，在实际石油工程中，为了确保开采过程的稳定性、生产优化和油藏保护，常常采用控制固定的井底压强或井体流量的方式，对于注入井而言，一般会设定固定流量的注入液注入，之后，井方程作为一个非线性函数，通常与压力方程饱和度方程耦合到同一个方程组中一同求解，以获取每一时刻的各井的井压力、各相压力、各相饱和度，即全耦合求解法。但是，注入井的位置是用户人为设定的，其数量是人为设定的，并不是在每一个油藏网格模型的网格中都设有注入井，在隐式求解过程中，牛顿法对应的雅可比矩阵仅包含稀疏的几行，全耦合的求解算法会导致注入井方程的变量个数与油藏网格块对应的变量个数不一致，如图3所示，即部分网格单元中有井压力参数，部分网格单元中因没有设定注入井则不计算井压力参数但是每个网格单元中都具有各相压力和各相饱和度等参数。因此，采用多个处理器并行求解时，每个处理器虽然平均分配了计算的网格单元，但是网格单元中需要计算的参数不完全相同，导致各个处理器的计算负载不均衡，需要通过填充网格来解决这一问题，即便如此，在多核计算时，仍然很难实现均匀分配计算负载给每个处理器，导致耦合方程系统在求解过程中存在负载不均衡的问题。

除此之外，在解决非线性耦合偏微分方程系统的挑战中，主要困难在于方程之间的相互关联和非线性耦合性质，系统的非线性特性导致全耦合方程的解析解难以获得，因此需要依赖数值方法进行求解，但是在数值求解过程中主要存在复杂性、计算效率低下的问题。

综合上述相关技术存在的问题，本发明旨在提供一种适用于油藏模拟中复杂多相多组分问题的大规模并行求解方法，从而兼顾油藏数值模拟的数值稳定性、时间步长宽松性和负载均衡性。

根据本发明实施例，提供了一种油藏数值模拟方法实施例，需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

在本实施例中提供了一种油藏数值模拟方法，可用于上述的计算机设备，图4是根据本发明实施例的一种油藏数值模拟方法的流程图，该流程包括如下步骤：

步骤S101，基于井压力创建用于耦合求解油藏的各相压力和各相饱和度的质量守恒方程，井压力是注入井的井底压力。

步骤S102，创建关于各相压力、各相饱和度和井压力的井方程；

步骤S103，初始化包括当前时间步的井压力的初始解向量；

步骤S104，利用初始解向量隐式求解质量守恒方程，得到下一时间步的各相压力和各相饱和度；

步骤S105，利用下一时间步的各相压力和各相饱和度显式求解井方程，得到下一时间步的井压力；

步骤S106，将下一时间步的井压力作为当前时间步的井压力，返回隐式求解质量守恒方程的步骤，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度。

具体地，本发明实施例针对油藏网格模型创建对应的数学模型后，对数学模型的处理分为两个部分，第一部分是基于井压力创建的质量守恒方程，一般情况下都是针对水相、气相和液相三相创建，有时也会只针对气相和液相两相创建进行分析。在该方程中，将待求解的各相的压力参数与饱和度参数耦合到一起，用于对各相压力和各相饱和度同时求解。然后，单独创建注入井的井方程，用于单独计算井压力，注入井是在油田开发过程中，为保持或恢复油层压力，在油田边缘或内部钻凿的往油层中注水或注气的井，注入井为油藏流体的产出提供了能量来源，对于注入井而言，一般在数值模拟时会在井方程中设定固定流量的注入液注入，从而方便计算井压力。而油藏模拟的井方程是在石油工程领域用于描述地下油藏中井筒内部流体流动的数学方程，通过油藏模拟的井方程，工程师们能够精确分析和优化井筒操作，包括注水、采油等过程，帮助决策者了解井筒产能、产液率、产气率以及压力变化等关键参数，从而制定更有效的生产策略，最大限度地提高油田的产量和采收率。

本实施例预先初始化至少包括当前时间步的井压力的初始解向量（根据实际创建的方程形式，初始解向量可中灵活添加其他参数），将初始解向量代入质量守恒方程中，从而对于各相压力和各相饱和度参数全耦合的质量守恒方程执行隐式求解，得到下一时间步的各相压力和各相饱和度。井方程相对简单，可利用前述步骤求解得到的下一时间步的各相压力和各相饱和度显式求解井方程，从而得到下一时间步的井压力。上述两个方程交叉迭代进行显隐计算，直至第预设数量个时间步停止迭代，得到最终的目标各相压力和目标各相饱和度。通过本发明实施例提供的求解方案，第一方面可通过隐式求解算法求解全耦合的质量守恒方程，令各相的压力以及饱和度计算变为全隐式计算，从而继承了隐式求解时间步更宽松的特点，得以使用更长的时间步，加快求解效率。第二方面，压力与饱和度的计算采用全隐式求解使数值求解的结果稳定性更高，减少振荡的问题。第三方面，本发明实施例提供的方案将井压力的计算过程从全耦合的方程中剔除出来，当质量守恒方程计算出下一时间步的各相压力和各相饱和度时，还单独对井方程进行显式求解，即使注入井的数量与油藏模型网格的数量并不一一对应，计算得到的下一时间步的井压力参数也不分配到油藏网格中，相比井方程与相压力方程耦合在一起的计算方法，当使用多个处理器执行并行计算时，不会出现网格计算负载分配不均匀的问题，提高了油藏数值模拟的负载均衡性。第四方面，考虑到全耦合计算稳定性高但是复杂性也高，从而计算效率低的特性，虽然本实施例将各相压力和各相饱和度的计算方程耦合在一起，但是拆分井压力的计算方程，实现了另一种迭代交替的计算方法，继承全耦合法数值稳定性的同时也降低了全耦合求解的复杂度，交替计算更简单，保证了油藏数值模拟的效率。

在一些可选地实施方式中，本发明实施例提供的同时耦合各相压力和各相饱和度的质量守恒方程如下：

其中，组分i在液相和气相的质量分数计算如下：

其中，各相流体流动均符合达西定律：

液相压力与相之间的毛细管压力满足如下守恒条件：

液、气两相逸度系数的计算：

分子间的吸引力和排斥力，计算方式如下：

在上述式子中，为孔隙度，、l和v分别表示水相、液相和气相，表示三相之一，、和分别为水相、液相和气相条件下的密度，、和分别为各相的流体速度，、和分别为各相饱和度，、和分别为各相压力，、和分别为各相的注入速率，由井压力计算得到，和分别表示液相和气相中第i个组分的质量分数，和分别表示液相和气相中第i个组分的摩尔分数，和表示液相和气相中第i个组分的逸度，和表示对应相逸度的系数值，表示的是计算逸度系数的给定参数，m为各相所含的组分数，其中为组分j的摩尔质量，为组分i的摩尔质量，和分别表示液相和气相中第j个组分的摩尔分数，为油藏的绝对渗透率，为重力加速度，、和为相流体的相对渗透率、粘性系数和移动率，为液相压力与相之间的毛细管压力，为压缩系数，和分别表示分子间的吸引力无量纲量和排斥力无量纲量，是液相第k个组分的摩尔分数，和分别为简化温度和压力值，和是中间变量，、和为计算中间变量和的经验参数，当i=1时，为前述，当i=2时，为前述，为组分i和k间的二元相互作用系数。

相应的，在本发明实施例中，独立创建的井方程如下：

其中，是相注入液的注入量，是相流体的注入速率，是第口井的射孔总数，是第口井的第个射孔的段长、和分别是井位置处相流体的绝对和相对渗透率值，和分别为井位置处相流体的粘性系数和密度，表示第口井的井筒半径，为第口井第个射孔所在网格块的排水半径，和分别表示井处的井压力深度和当前射孔位置深度，为井位置处的井压力,为井处相流体的压力，是重力加速度。

具体地，本发明实施例假设油藏流体流动过程是在恒定温度条件进行，有限个烃类组分共同组成储层流体，其中组分至多形成如图5所示的三相，烃相(即气相和液相)与水相之间无质量交换，每个组分都可能存在于这些相中的任何一个或全部，从而创建质量守恒方程与井方程，提高有藏数值模拟的准确度。基于上述场景条件，待求解参数除井压力，各相饱和度、和，各相压力、和外，还包括液相和气相的组分摩尔分数和。

相应地，针对上述质量守恒方程和井方程，如图6所示，本发明实施例提供的油藏数值模拟方法具体求解步骤如下：

步骤S603，初始化包括当前时间步的井压力的初始解向量。具体地，详细请参见步骤S103，在此不再赘述。

步骤S603具体包括：

步骤a1，初始化当前时间步的井压力、当前时间步的各相压力、当前时间步的各相饱和度、当前时间步的液相中各组分的摩尔分数、当前时间步的气相中各组分的摩尔分数，得到初始解向量(,,,,)。

步骤S604，利用初始解向量隐式求解质量守恒方程，得到下一时间步的各相压力和各相饱和度。具体地，详细请参见步骤S104，在此不再赘述。

步骤S604具体包括：

步骤a2，固定各相注入液的注入量，将初始解向量代入井方程，显式计算当前时间步的各相注入速率；

步骤a3，将当前时间步的各相注入速率代入质量守恒方程，进行隐式求解，得到下一时间步的各相压力、各相饱和度、液相中各组分的摩尔分数和气相中各组分的摩尔分数。

步骤S605，利用下一时间步的各相压力和各相饱和度显式求解井方程，得到下一时间步的井压力。具体地，详细请参见步骤S105，在此不再赘述。

步骤S605具体包括：

步骤a4，将下一时间步的各相压力、各相饱和度、液相中各组分的摩尔分数和气相中各组分的摩尔分数代入井方程，进行显式求解，得到下一时间步的井压力。

步骤S606，将下一时间步的井压力作为当前时间步的井压力，返回隐式求解质量守恒方程的步骤，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度。具体地，详细请参见步骤S106，在此不再赘述。

步骤S606具体包括：

步骤a5，将下一时间步的井压力作为当前时间步的井压力，重复步骤a2到a4，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度。

具体地，在本发明实施例中，针对上述质量守恒方程和井方程，如图7所示，计算过程为：用户首先定义初始解向量，该初始解向量中包括当前时间步的井压力、当前时间步的各相压力、当前时间步的各相饱和度、当前时间步的液相中各组分的摩尔分数、当前时间步的气相中各组分的摩尔分数。之后，将初始解向量(,,,,)代入井方程的右侧，，其中注入量根据专家经验获取经验参数，从而固定为定值，然后通过上述等式关系显式计算当前时间步各相的注入速率。得到后，将代入第一部分质量守恒方程，通过隐式计算软件同时计算下一个时间步的各相压力，各相饱和度，液相中各组分的摩尔分数，气相中各组分的摩尔分数。之后，继续将(,,,)代入井方程，通过显式求解更新下一个时间步的井压力，继续用更新后的井压力代入质量守恒方程，隐式求解再下一个时间步耦合的各相压力和各相饱和度，通过多次迭代，直至经过预设数量时间步（例如100），最终得到准确的目标各相压力和目标各相饱和度。

通过本实施例提供的方案，隐式求解全耦合的质量守恒方程，令各相的压力以及饱和度计算变为全隐式计算，从而继承了隐式求解时间步更宽松的特点，得以使用更长的时间步，加快求解效率。第二方面，压力与饱和度的计算采用全隐式求解使数值求解的结果稳定性更高，减少振荡的问题。第三方面，本发明实施例提供的方案将井压力的计算过程从全耦合的方程中剔除出来，单独对井方程进行显式求解，即使注入井的数量与油藏模型网格的数量并不一一对应，计算得到的下一时间步的井压力参数也不分配到油藏网格中，相比井方程与相压力方程耦合在一起的计算方法，当使用多个处理器执行并行计算时，不会出现网格计算负载分配不均匀的问题，提高了油藏数值模拟的负载均衡性。第四方面，虽然本实施例将各相压力和各相饱和度的计算方程耦合在一起，但是拆分井压力的计算方程，实现了另一种迭代交替的计算方法，继承全耦合法数值稳定性的同时也降低了全耦合求解的复杂度，交替计算更简单，保证了油藏数值模拟的效率。

在一些可选地实施方式中，本发明实施例通过用牛顿-克雷洛夫-施瓦兹算法对质量守恒方程进行隐式求解。

具体地，在应对如本发明实施例提供的质量守恒方程这类大规模非线性方程组的挑战时，牛顿Newton-克雷洛夫Krylov-施瓦兹Schwarz算法作为一种高效的并行优化求解方法适用性更高。该算法融合了牛顿法、Krylov子空间方法和Schwarz区域分解技术。该算法的具体步骤包括：

1.初始化，给定初始解，设置循环次数k=0；

2.通过下式确定搜索方向，计算非精确的牛顿寻优方向，令

3.通过下式确定搜索步长，利用线搜索方法确定最佳搜索步长

4.通过下式更新解

5.通过下式的停止搜索条件判断是否停止迭代

6.当不满足条件时，设置循环次数k=k+1，返回步骤2。

上式中，，表示在解处的雅可比矩阵,为限制性施瓦兹预条件子,和分别为线性系统求解的相对误差和绝对误差参数，为线搜索控制参数，和分别表示非线性系统求解的相对误差和绝对误差参数，表示偏微分方程离散后的方程系统。

当求解线性方程组时，Krylov子空间方法通过在每次迭代中建立Krylov子空间，并逐步逼近精确解以找到最小残差向量的线性组合。鉴于Krylov子空间方法具备广泛适用性、良好的迭代收敛性和高效的内存利用率等优点，本实施例选择该方法来求解线性系统。同时，为了进一步提高线性系统的求解效率，本实施例还引入施瓦兹预条件子，施瓦兹预条件子是一种基于区域分解的技术，将求解域分割成多个子域，并在每个子域上构建独立的预条件子，从而充分利用线性方程组的结构特点和局部信息，能够有效地加快Krylov子空间方法的收敛速度，减轻整体迭代的计算负担。

在本实施例中，一种限制性施瓦兹预条件子表达式如下：

，

其中和分别为重叠和非重叠子区域的限制算子,表示的转置，矩阵为全局矩阵限制后的子矩阵,全局矩阵A表示为离散后方程系统的雅可比矩阵。表示计算该问题所设定的处理器总个数；表示中的第个处理器序号。综合采用Newton-Krylov-Schwarz求解质量守恒方程，将多个求解方法相结合，充分发挥各自的优点，能够高效地求解复杂的质量守恒方程，并具备良好的可扩展性和并行计算能力，提高并行计算的速度。

在一个实际应用场景实施例中，通过本发明实施例提供的方案进行油藏数值模拟，基于PETSc工具包在计算机上实现了本发明实施例所提出的显隐交替求解方法，经过两相六组分测试实验验证后，结果表明石油产出率数值结果与MRST提供的参照解完全匹配，且数值结果稳定。该方法支持多核计算，适用于油藏模拟中多相多组分问题的大规模并行求解方法。

在本实施例中还提供了一种油藏数值模拟装置，该装置用于实现上述实施例及优选实施方式，已经进行过说明的不再赘述。如以下所使用的，术语“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

本实施例提供一种油藏数值模拟装置，如图8所示，包括：

质量守恒方程创建模块801，用于基于井压力创建用于耦合求解油藏的各相压力和各相饱和度的质量守恒方程，井压力是注入井的井底压力。详细内容参见上述方法实施例中步骤S101的相关描述，在此不再进行赘述。

井方程创建模块802，用于创建关于各相压力、各相饱和度和井压力的井方程。详细内容参见上述方法实施例中步骤S102的相关描述，在此不再进行赘述。

初始解模块803，用于初始化包括当前时间步的井压力的初始解向量。详细内容参见上述方法实施例中步骤S103的相关描述，在此不再进行赘述。

隐式求解模块804，用于利用初始解向量隐式求解质量守恒方程，得到下一时间步的各相压力和各相饱和度。详细内容参见上述方法实施例中步骤S104的相关描述，在此不再进行赘述。

显式求解模块805，用于利用下一时间步的各相压力和各相饱和度显式求解井方程，得到下一时间步的井压力。详细内容参见上述方法实施例中步骤S105的相关描述，在此不再进行赘述。

反馈迭代求解模块806，用于将下一时间步的井压力作为当前时间步的井压力，返回隐式求解质量守恒方程的步骤，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度。详细内容参见上述方法实施例中步骤S106的相关描述，在此不再进行赘述。

上述各个模块和单元的更进一步的功能描述与上述对应实施例相同，在此不再赘述。

本实施例中的一种油藏数值模拟装置是以功能单元的形式来呈现，这里的单元是指ASIC（Application Specific Integrated Circuit，专用集成电路）电路，执行一个或多个软件或固定程序的处理器和存储器，和/或其他可以提供上述功能的器件。

本发明实施例还提供一种计算机设备，具有上述图8所示的一种油藏数值模拟装置。

请参阅图9，图9是本发明可选实施例提供的一种计算机设备的结构示意图，如图9所示，该计算机设备包括：一个或多个处理器10、存储器20，以及用于连接各部件的接口，包括高速接口和低速接口。各个部件利用不同的总线互相通信连接，并且可以被安装在公共主板上或者根据需要以其它方式安装。处理器可以对在计算机设备内执行的指令进行处理，包括存储在存储器中或者存储器上以在外部输入/输出装置(诸如，耦合至接口的显示设备)上显示GUI的图形信息的指令。在一些可选的实施方式中，若需要，可以将多个处理器和/或多条总线与多个存储器和多个存储器一起使用。同样，可以连接多个计算机设备，各个设备提供部分必要的操作(例如，作为服务器阵列、一组刀片式服务器、或者多处理器系统)。图9中以一个处理器10为例。

处理器10可以是中央处理器，网络处理器或其组合。其中，处理器10还可以进一步包括硬件芯片。上述硬件芯片可以是专用集成电路，可编程逻辑器件或其组合。上述可编程逻辑器件可以是复杂可编程逻辑器件，现场可编程逻辑门阵列，通用阵列逻辑或其任意组合。

其中，存储器20存储有可由至少一个处理器10执行的指令，以使至少一个处理器10执行实现上述实施例示出的方法。

存储器20可以包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需要的应用程序；存储数据区可存储根据计算机设备的使用所创建的数据等。此外，存储器20可以包括高速随机存取存储器，还可以包括非瞬时存储器，例如至少一个磁盘存储器件、闪存器件、或其他非瞬时固态存储器件。在一些可选的实施方式中，存储器20可选包括相对于处理器10远程设置的存储器，这些远程存储器可以通过网络连接至该计算机设备。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。

存储器20可以包括易失性存储器，例如，随机存取存储器；存储器也可以包括非易失性存储器，例如，快闪存储器，硬盘或固态硬盘；存储器20还可以包括上述种类的存储器的组合。

该计算机设备还包括通信接口30，用于该计算机设备与其他设备或通信网络通信。

本发明实施例还提供了一种计算机可读存储介质，上述根据本发明实施例的方法可在硬件、固件中实现，或者被实现为可记录在存储介质，或者被实现通过网络下载的原始存储在远程存储介质或非暂时机器可读存储介质中并将被存储在本地存储介质中的计算机代码，从而在此描述的方法可被存储在使用通用计算机、专用处理器或者可编程或专用硬件的存储介质上的这样的软件处理。其中，存储介质可为磁碟、光盘、只读存储记忆体、随机存储记忆体、快闪存储器、硬盘或固态硬盘等；进一步地，存储介质还可以包括上述种类的存储器的组合。可以理解，计算机、处理器、微处理器控制器或可编程硬件包括可存储或接收软件或计算机代码的存储组件，当软件或计算机代码被计算机、处理器或硬件访问且执行时，实现上述实施例示出的方法。

虽然结合附图描述了本发明的实施例，但是本领域技术人员可以在不脱离本发明的精神和范围的情况下做出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。

Claims (8)

1.一种油藏数值模拟方法，其特征在于，所述方法包括：

基于井压力创建用于耦合求解油藏的各相压力和各相饱和度的质量守恒方程，所述井压力是注入井的井底压力；

创建关于各相压力、各相饱和度和井压力之间数值关系的井方程；

初始化包括当前时间步的井压力的初始解向量；

利用所述初始解向量隐式求解所述质量守恒方程，得到下一时间步的各相压力和各相饱和度；

利用所述下一时间步的各相压力和各相饱和度显式求解所述井方程，得到下一时间步的井压力；

将所述下一时间步的井压力作为所述当前时间步的井压力，返回所述隐式求解所述质量守恒方程的步骤，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度；

所述质量守恒方程如下

其中，为孔隙度，、l和v分别表示水相、液相和气相，表示三相之一，、和分别为水相、液相和气相条件下的密度，、和分别为各相的流体速度，、和分别为各相饱和度，、和分别为各相压力，、和分别为各相的注入速率，由井压力计算得到，和分别表示液相和气相中第i个组分的质量分数，和分别表示液相和气相中第i个组分的摩尔分数，和表示液相和气相中第i个组分的逸度，和表示液相和气相逸度的系数值，表示的是计算逸度系数的给定参数，m为各相所含组分的总数量，其中为组分j的摩尔质量，为组分i的摩尔质量，和分别表示液相和气相中第j个组分的摩尔分数，为油藏的绝对渗透率，为重力加速度，、和为相流体的相对渗透率、粘性系数和移动率，为液相压力与相之间的毛细管压力，为压缩系数，和分别表示分子间的吸引力无量纲量和排斥力无量纲量，是液相第k个组分的摩尔分数，和分别为组分i的简化温度和压力值，和是中间变量，、和为计算中间变量和的经验参数，为组分i和k间的二元相互作用系数；

所述井方程如下

其中，是相注入液的注入量，是相流体的注入速率，是第口井的射孔总数，是第口井的第个射孔的段长、 和分别是井位置处相流体的绝对渗透率值和相对渗透率值，和分别为井位置处相流体的粘性系数和密度，表示第口井的井筒半径，为第口井第个射孔所在网格块的排水半径，和分别表示井处的井压力深度和当前射孔位置深度，为井位置处的井压力，为井处相流体的压力，是重力加速度。

2.根据权利要求1所述的方法，其特征在于，所述初始化包括当前时间步的井压力的初始解向量，包括：

初始化当前时间步的井压力、当前时间步的各相压力、当前时间步的各相饱和度、当前时间步的液相中各组分的摩尔分数、当前时间步的气相中各组分的摩尔分数，得到所述初始解向量。

3.根据权利要求2所述的方法，其特征在于，所述利用所述初始解向量隐式求解所述质量守恒方程，得到下一时间步的各相压力和各相饱和度，包括：

固定各相注入液的注入量，将所述初始解向量代入所述井方程，显式计算当前时间步的各相注入速率；

将所述当前时间步的各相注入速率代入所述质量守恒方程，进行隐式求解，得到下一时间步的各相压力、各相饱和度、液相中各组分的摩尔分数和气相中各组分的摩尔分数。

4.根据权利要求3所述的方法，其特征在于，所述利用所述下一时间步的各相压力和各相饱和度显式求解所述井方程，包括：

将所述下一时间步的各相压力、各相饱和度、液相中各组分的摩尔分数和气相中各组分的摩尔分数代入所述井方程，进行显式求解，得到下一时间步的井压力。

5.根据权利要求3所述的方法，其特征在于，通过牛顿-克雷洛夫-施瓦兹算法对所述质量守恒方程进行隐式求解。

6.一种油藏数值模拟装置，其特征在于，所述装置包括：

质量守恒方程创建模块，用于基于井压力创建用于耦合求解油藏的各相压力和各相饱和度的质量守恒方程，所述井压力是注入井的井压力；

井方程创建模块，用于创建关于各相压力、各相饱和度和井压力的井方程；

初始解模块，用于初始化包括当前时间步的井压力的初始解向量；

隐式求解模块，用于利用所述初始解向量隐式求解所述质量守恒方程，得到下一时间步的各相压力和各相饱和度；

显式求解模块，用于利用所述下一时间步的各相压力和各相饱和度显式求解所述井方程，得到下一时间步的井压力；

反馈迭代求解模块，用于将所述下一时间步的井压力作为所述当前时间步的井压力，返回所述隐式求解所述质量守恒方程的步骤，直至经过预设数量时间步的迭代，得到目标各相压力和目标各相饱和度；

所述质量守恒方程如下

其中，为孔隙度，、l和v分别表示水相、液相和气相，表示三相之一，、和分别为水相、液相和气相条件下的密度，、和分别为各相的流体速度，、和分别为各相饱和度，、和分别为各相压力，、和分别为各相的注入速率，由井压力计算得到，和分别表示液相和气相中第i个组分的质量分数，和分别表示液相和气相中第i个组分的摩尔分数，和表示液相和气相中第i个组分的逸度，和表示液相和气相逸度的系数值，表示的是计算逸度系数的给定参数，m为各相所含组分的总数量，其中为组分j的摩尔质量，为组分i的摩尔质量，和分别表示液相和气相中第j个组分的摩尔分数，为油藏的绝对渗透率，为重力加速度，、和为相流体的相对渗透率、粘性系数和移动率，为液相压力与相之间的毛细管压力，为压缩系数，和分别表示分子间的吸引力无量纲量和排斥力无量纲量，是液相第k个组分的摩尔分数，和分别为组分i的简化温度和压力值，和是中间变量，、和为计算中间变量和的经验参数，为组分i和k间的二元相互作用系数；

所述井方程如下

其中，是相注入液的注入量，是相流体的注入速率，是第口井的射孔总数，是第口井的第个射孔的段长、 和分别是井位置处相流体的绝对渗透率值和相对渗透率值，和分别为井位置处相流体的粘性系数和密度，表示第口井的井筒半径，为第口井第个射孔所在网格块的排水半径，和分别表示井处的井压力深度和当前射孔位置深度，为井位置处的井压力，为井处相流体的压力，是重力加速度。

7.一种计算机设备，其特征在于，包括：

存储器和处理器，所述存储器和所述处理器之间互相通信连接，所述存储器中存储有计算机指令，所述处理器通过执行所述计算机指令，从而执行权利要求1至5中任一项所述的方法。

8.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机指令，所述计算机指令用于使计算机执行权利要求1至5中任一项所述的方法。