基于组分方程的油气运移聚集数值模拟方法及系统与流程

1.本发明属于数值模拟技术领域，更具体地，涉及一种基于组分方程的油气运移聚集数值模拟方法及系统。


背景技术：

2.油气运移聚集分析是石油地质研究中的重要内容，而油气运移聚集的数值模拟是盆地模拟研究中的重要组成部分，可用于预测油气的聚集场所和聚集量。在现有的盆地模拟中，所采用的油气运移聚集数值模拟方法主要包括流线法、侵入渝渗法和达西渗流方程求解法三种。
3.其中，流线法是最简单的且计算速度最快的方法，通常被用于快速的筛选和评价。流线法典型地应用于高渗透储层或者疏导层。在应用时，通常采用设置渗透率门限值的方式来区分疏导层和低渗透层或者直接人工定义给出疏导层和低渗透层。当油气进入高渗透层时，浮力被认为是驱动力，油气垂向运移到下一个上覆的岩石单元直至遇到毛管进入压力较高的低渗透层，阻止其进一步的运动。流线法不考虑时间或流动速率。
4.侵入渝渗法通过绘制每一个模型网格体的突破毛管压力分布图来模拟地下油气的运移。虽然侵入渝渗法不是一个物理上的精确的数值模拟方法，但其产生的运移类型与实验吻合地很好，对多种地质场景能给出满意的结果，并且数值方法简单，容易编程。另外，侵入渝渗法对于储层和非储层采取相同的处理方式，因此也能够模拟油气在低渗透岩性中的运移。
5.由于流线法和侵入渝渗法采用的都是基于静力的模型，因此，流线法和侵入渝渗法中均没有时间或流动速率的概念，与粘度和渗透率等参数无关。同时，流线法和侵入渝渗法都是将油气的运移过程和聚集过程分开处理，油气在网格体中的饱和度通常都是给一个固定的值，因此这两种方法都不能精确地模拟储层含油饱和度的分布。另外，这两种方法都无法处理油气在源岩中的运聚或排出过程。
6.达西渗流方程求解法在盆地模拟中能处理自然界中范围广泛的岩石渗透率。达西渗流方程求解法的计算时间随模拟网格数的增加而迅速增加，因此，通常使用较低精度的网格。达西渗流方程求解法是三种方法中唯一能够模拟真实油气运移和聚集速率的方法，能够模拟异常压力，同时也是模拟深盆气等动态气藏形成和发展的唯一方法，并且支持油气的向下运移。因此，达西渗流方程求解法具有独特性和不可替代性。
7.在盆地模拟中，基于达西渗流的油气运移聚集数值模拟的实现方法很多，例如petromod盆地模拟软件采用基于有限元的求解方法，而更多的是基于有限体积法的模拟技术(wendebourg&harbaugh,1997；袁益让和韩玉笈，2008；石广仁等，2010；郭秋麟等，2015)，在有限体积解法中，前两个文献是求解油、水两相问题的，后两个文献求解的是黑油模型。
8.由于盆地中烃类存在的相态的多样性，而且由于热演化的程度不同导致油气组分变化很大，因此，基于组分守恒的达西渗流求解的数值模拟变得很有必要。多组分达西渗流求解更多地见于油藏模拟中，通常采用基于组分摩尔守恒的偏微分方程，并基于组分的摩
尔比作为基本变量进行求解(郎兆新，程林松，1990；张茂林等,1991；叶继根，吴向红，2000；张茂林等，2002；chen，2007)。
9.在盆地模拟中会遇到各种类型的盆地，并且模拟地层由浅到深的埋藏过程，因此生成的烃类组分也随着源岩深度的逐渐增加由重烃组分为主逐渐向轻烃组分为主演化，同时已经生成的烃类也会随着地层温度的增高进行烃类的热裂解，因此在整个盆地模拟过程中烃类的组分含量变化通常是非常大的，通常用于描述单一油气相态的用于油藏模拟的油水或气水两相模型或黑油模型等都无法满足现有盆地模拟的对油气两种相态同时存在或在演化过程中相态变化的多种情况进行模拟的需求，因此，亟需一种能够处理烃类组分和相态变化的数值模拟模型。


技术实现要素：

10.本发明的目的在于解决现有基于达西渗流的油气运移聚集数值模拟方法无法准确地对盆地演化中组分含量多变的烃类进行模拟的问题。
11.为了实现上述目的，本发明提供一种基于组分方程的油气运移聚集数值模拟方法及系统。
12.根据本发明的第一方面，提供了一种基于组分方程的油气运移聚集数值模拟方法，所述基于组分方程的油气运移聚集数值模拟方法包括以下步骤：
13.建立基本方程和辅助方程，所述基本方程为基于组分质量平衡的偏微分方程组，所述偏微分方程组包括累积项、达西渗流项、压实项和源汇项，所述辅助方程包括相平衡方程、饱和度约束方程、毛细管压力方程和质量分数约束方程；
14.采用有限体积法将所述偏微分方程组离散化，得到有限体积方程；
15.对所述有限体积方程进行全隐式求解，以获知所述偏微分方程组的主要未知量；
16.根据获知的所述主要未知量，并基于所述辅助方程获知所述偏微分方程组的其余未知量；
17.根据已求解的所述偏微分方程组进行油气运移聚集的数值模拟。
18.作为优选的是，所述基于组分质量平衡的偏微分方程组的表达式为：
[0019][0020]
上式中，设l＝1,
…
,n-1为烃组分，l＝n为水组分；
[0021]
为梯度算子，uo、uw和ug分别为油超压、水超压和气超压，φ为岩石的有效孔隙度，so、sw和sg分别为油饱和度、水饱和度、气饱和度，u
l
为静岩势，k
ro
、k
rw
和k
rg
分别为油相流体的相对渗透率、水相流体的相对渗透率和气相流体的相对渗透率，k为岩石的绝对渗透率，ρo、ρw和ρg分别为油相流体的密度、水相流体的密度和气相流体的密度，μo、μw和μg分别为油相流体的粘度、水相流体的粘度和气相流体的粘度，qw为水的源汇强度，q
l
表示烃组分的源汇强度，σe为有效应力，c
lw
为水相中组分l的质量分数，c
lo
为油相中组分l的质量分数，c
lg
为气相中组分l的质量分数。
[0022]
作为优选的是，所述相平衡方程的表达式为：
[0023][0024][0025]
上式中，k
lgo
为组分l在气相与油相间的平衡常数，k
lgw
为组分l在气相与水相间的平衡常数，f(x)为压力、温度和组分含量的函数。
[0026]
作为优选的是，所述饱和度约束方程的表达式为：
[0027]
so+sw+sg＝1
ꢀꢀꢀ
(4)。
[0028]
作为优选的是，所述毛细管压力方程的表达式为：
[0029]uo-uw＝p
co
(sw)+p
wstatic-p
ostatic
ꢀꢀꢀ
(5)
[0030]ug-uo＝p
cg
(so)+p
ostatic-p
gstatic
ꢀꢀꢀ
(6)
[0031]
上式中，p
wstatic
、p
ostatic
和p
gstatic
分别为静水压力、静油压力和静气压力，p
co
(sw)为油相所受的毛管压力，p
cg
(so)为气相所受的毛管压力。
[0032]
作为优选的是，所述质量分数约束方程的表达式为：
[0033][0034]
上式中，o代表油相，g代表气相，w代表水相。
[0035]
作为优选的是，所述有限体积方程的表达式为：
[0036][0037]
上式中，d
i，j，k
代表目标计算区域的中心点坐标为(i，j，k)的网格单元的控制体积。
[0038]
作为优选的是，所述对所述有限体积方程进行全隐式求解，具体为：
[0039]
所述有限体积方程的隐式迭代格式为：
[0040][0041]
上式中，v
i,j,k
为中心点坐标为(i,j,k)的网格单元的体积，t
n+1
和tn均为时间步，l和l+1均为迭代次数。
[0042]
作为优选的是，所述偏微分方程组的全部未知量包括流体的超压、流体的饱和度和相中的组分质量分数；
[0043]
所述主要未知量包括水的超压和烃组分质量。
[0044]
根据本发明的第二方面，提供了一种基于组分方程的油气运移聚集数值模拟系统，所述基于组分方程的油气运移聚集数值模拟系统包括处理器和存储器，所述处理器执行所述存储器中保存的计算机程序时实现上述任一种油气运移聚集数值模拟方法。
[0045]
本发明的基于组分方程的油气运移聚集数值模拟方法，首先建立用于模拟油气运移聚集的基于组分质量平衡的偏微分方程组和辅助方程；其次采用有限体积法将所述偏微分方程组离散化，得到有限体积方程；再次对所述有限体积方程进行全隐式求解，以获知所述偏微分方程组的主要未知量；从次根据获知的所述主要未知量，并基于所述辅助方程获知所述偏微分方程组的其余未知量；最后根据已求解的所述偏微分方程组进行油气运移聚集的数值模拟。
[0046]
本发明的有益效果在于：
[0047]
一、现有基于达西渗流的油气运移聚集数值模拟方法采用的都是两相流或基于黑油模型方程的解法，因此，无法完全描述实际地质条件下油气相态多变、组分组成差别可能巨大的地质实际。而本发明的基于组分方程的油气运移聚集数值模拟方法采用多组分的解法，因此能对盆地演化中烃类组成变化多样的特点进行准确的模拟。
[0048]
二、现在对于多组分数学模型的建立以及数值求解都是建立在组分摩尔守恒且以组分摩尔比为基本变量进行求解的基础上，而本发明的基于组分方程的油气运移聚集数值模拟方法适应盆地模拟中的实际地质情形，采用基于网格体组分质量平衡的数学模型，并以网格体组分质量为基本变量进行求解，因此求解结果直接反映到网格单元的油气质量分布上。
[0049]
本发明的基于组分方程的油气运移聚集数值模拟系统与上述基于组分方程的油气运移聚集数值模拟方法属于一个总的发明构思，故与上述基于组分方程的油气运移聚集数值模拟方法具有相同的有益效果，在此不再赘述。
[0050]
本发明的其它特征和优点将在随后具体实施方式部分予以详细说明。
附图说明
[0051]
通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。
[0052]
图1示出了根据本发明的实施例1的基于组分方程的油气运移聚集数值模拟方法的实现流程图；
[0053]
图2示出了根据本发明的实施例1的典型六面体单元的示意图；
[0054]
图3示出了根据本发明的实施例1的塔里木盆地塔中斜坡区的地质图，其中，小方框所框选的区域为具体示例1所模拟的区域；
[0055]
图4示出了根据本发明的实施例1的中下寒武地层在两个不同时期的含油饱和度分布图，其中，左图为中下寒武地层在224百万年的含油饱和度分布图，右图为中下寒武地层在206百万年的含油饱和度分布图；
[0056]
图5示出了根据本发明的实施例1的下寒武统地层和中下奥陶统地层在433百万年的含油饱和度分布图，其中，左图为下寒武统地层在433百万年的含油饱和度分布图，右图为中下奥陶统地层在433百万年的含油饱和度分布图。
具体实施方式
[0057]
下面将更详细地描述本发明的优选实施方式。虽然以下描述了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
[0058]
实施例1：本发明实施例的基于组分方程的油气运移聚集数值模拟方法，是一种能够处理复杂地质场景下油气运移和聚集的基于达西渗流模型的三维三相多组分的偏微分方程的有限体积法求解方法。这里的地质场景可以是一个概念化的地质体、一个油气系统或者一个盆地。
[0059]
本发明实施例的基于组分方程的油气运移聚集数值模拟方法中，用于模拟油气运移聚集的偏微分方程组是基于组分的质量守恒而得到的，包含了累积项、达西渗流项、压实项和源汇项。同时考虑的辅助方程有相平衡方程、饱和度约束方程、毛细管压力方程和质量分数约束方程。偏微分方程组待求的变量包括流体的超压、流体的饱和度以及相中的组分质量分数。本发明实施例的基于组分方程的油气运移聚集数值模拟方法通过有限体积法将所述偏微分方程组离散化。由于积分守恒对任意一组控制体积都得到满足，因此整个计算区域也满足守恒。本发明实施例给出了方程各项的有限体积法离散过程和结果。
[0060]
由于盆地模拟涉及的通常是大时间尺度的数值模拟过程，对模拟过程的稳定性要求高，因此，本发明实施例给出的是有限体积方程的全隐式解法。在具体的求解过程中，主要未知变量为水的超压和烃组分的质量，其它的变量通过与压力、温度和烃组分质量间的函数关系而求解。这些函数能通过质量和质量分数关系、辅助方程以及其它如相态的密度、粘度等参数求解函数而形成函数关系。
[0061]
图1示出了本发明实施例的基于组分方程的油气运移聚集数值模拟方法的实现流程图。参照图1，本发明实施例的基于组分方程的油气运移聚集数值模拟方法包括以下步
骤：
[0062]
步骤s100、建立基本方程和辅助方程，所述基本方程为基于组分质量平衡的偏微分方程组，所述偏微分方程组包括累积项、达西渗流项、压实项和源汇项，所述辅助方程包括相平衡方程、饱和度约束方程、毛细管压力方程和质量分数约束方程；
[0063]
步骤s200、采用有限体积法将所述偏微分方程组离散化，得到有限体积方程；
[0064]
步骤s300、对所述有限体积方程进行全隐式求解，以获知所述偏微分方程组的主要未知量；
[0065]
步骤s400、根据获知的所述主要未知量，并基于所述辅助方程获知所述偏微分方程组的其余未知量；
[0066]
步骤s500、根据已求解的所述偏微分方程组进行油气运移聚集的数值模拟。
[0067]
进一步地，本发明实施例的步骤s100中，所述基于组分质量平衡的偏微分方程组的表达式为：
[0068][0069][0070]
上式中，设l＝1,
…
,n-1为烃组分，l＝n为水组分；
[0071]
为梯度算子，uo、uw和ug分别为油超压、水超压和气超压，[pa]；φ为岩石的有效孔隙度，无量纲；so、sw和sg分别为油饱和度、水饱和度、气饱和度，无量纲；u
l
为静岩势，无量纲，等于静岩压力-静水压力；k
ro
、k
rw
和k
rg
分别为油相流体的相对渗透率、水相流体的相对渗透率和气相流体的相对渗透率，无量纲；k为岩石的绝对渗透率，[m2]；ρo、ρw和ρg分别为油相流体的密度、水相流体的密度和气相流体的密度，[kg/m3]；μo、μw和μg分别为油相流体的粘度、水相流体的粘度和气相流体的粘度，[帕斯卡.秒，pa.s]；qw为水的源汇强度，[kg/(s.m3)]；q
l
表示烃组分的源汇强度，[kg/(s.m3)]；σe为有效应力，等于u
l-uw；c
lw
为水相中组分l的质量分数，c
lo
为油相中组分l的质量分数，c
lg
为气相中组分l的质量分数。
[0072]
再进一步地，本发明实施例的步骤s100中，所述相平衡方程的表达式为：
[0073][0074]
[0075]
上式中，k
lgo
为组分l在气相与油相间的平衡常数，k
lgw
为组分l在气相与水相间的平衡常数，f(x)为压力、温度和组分含量的函数。
[0076]
在每两相之间，对每一组分有一个平衡常数，这个平衡常数是压力、温度和组分含量的函数。
[0077]
再进一步地，本发明实施例的步骤s100中，所述饱和度约束方程的表达式为：
[0078]
so+sw+sg＝1
ꢀꢀꢀ
(4)。
[0079]
再进一步地，本发明实施例的步骤s100中，所述毛细管压力方程的表达式为：
[0080]uo-uw＝p
co
(sw)+p
wstatic-p
ostatic
ꢀꢀꢀ
(5)
[0081]ug-uo＝p
cg
(so)+p
ostatic-p
gstatic
ꢀꢀꢀꢀ
(6)
[0082]
上式中，p
wstatic
、p
ostatic
和p
gstatic
分别为静水压力、静油压力和静气压力，p
co
(sw)为油相所受的毛管压力，p
cg
(so)为气相所受的毛管压力。
[0083]
再进一步地，本发明实施例的步骤s100中，所述质量分数约束方程的表达式为：
[0084][0085]
上式中，o代表油相，g代表气相，w代表水相。
[0086]
再进一步地，本发明实施例的步骤s200为建立有限体积方程。有限体积法的基本思想是将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待求解的基于组分质量平衡的偏微分方程组对每一个控制体积积分，并得出一组离散方程，其中的未知数是网格点上的因变量的值。为了求出控制体积的积分，必须获知假定值在网格点之间的变化规律，即假设值的分段分布特征。
[0087]
考虑控制体为六面体的情形，将研究区域剖分成如图2所示的六面体单元，对于单元(i，j，k)，在(i，j，k)控制体积(d
i，j，k
)上，对偏微分方程组进行积分，展开整理得到：
[0088][0089]
下面分别给出针对内部节点的控制体单元的各项：
[0090]
累积项(即储存变化量)：
[0091][0092][0093][0094]
单元边界流入项(即达西渗流项)：
[0095][0096][0097][0098][0099]
其中，
[0100][0101][0102][0103]
压实项：
[0104][0105]
源汇项：
[0106][0107]
记：
[0108][0109]
得到：
[0110]
[0111]
记：
[0112][0112]
δ
x
u＝u
i+1-ui，δ
x
txδ
x
ui＝t
i+1/2
(u
i+1
—ui)+t
i-1/2
(u
i-1-ui)，δyu＝u
j+1-uj，δytyδyuj＝t
j+1/2
(u
j+1-uj)+t
j-1/2
(u
j-1-uj)δzu＝u
k+1-uk，δztzδzuk＝t
k+1/2
(u
k+1-uk)+t
k-1/2
(u
k-1-uk)δtδu＝δ
x
txδ
x
u+δytyδyu+δztzδzu以及
[0113]
得到的方程(18)简写形式：
[0114][0115]
再进一步地，本发明实施例的步骤s300为有限体积方程的全隐式求解。对于任意一个从tn到t
n+1
的时间步长，用全隐式迭代求解时，设对任意变量u，定义两次迭代l和l+1之间u的差值为：
[0116][0117]
即
[0118][0119]
而且，当l＝0时，u0＝un，经过多次迭代，当|u
l+1-u
l
|＜ε时，u
n+1
＝u
l+1
。其中ε为所给的精度要求。则迭代格式可写成：
[0120][0121]
将对流项和累积项展开的有限体积方程的隐式迭代格式如下：
[0122][0123]
利用隐含函数关系，获得关于利用隐含函数关系，获得关于与和的函数关系。这里的mr为每种组分在网格体中的质量。通过一系列的简化处理，可获得精简后的方程组：
[0124][0125]
方程组矩阵元素的具体计算：
[0126]
上面的方程与差分有关的一些项需展开，包括左端项未知项的系数和右端项的r4
li,j,k
、r5
li,j,k
、r7
i,j,k
、r8
i,j,k
，涉及到及(i,j,k)周围6个节点的值，需要展开整理。需要特别注意的是有三种情形：
[0127]
①
对于节点(i,j,k)，与储存量和压实项有关的部分，是直接计算该节点处的未知数的系数和相应的部分右端项；
[0128]
②
与对流项有关的部分，如果表达式涉及到节点(i,j,k)及其周围6个节点，这时就同时计算出7个点未知数的系数或7个点的右端项；
[0129]
③
与对流项有关的部分，如果表达式涉及到节点(i,j,k)和其周围6个节点之间1/2处的值，这时也同时计算出7个点或6个点的右端项；具体情况依照算法的不同有所差异。如果采用上游权方法，则可能涉及到7个点，如果简单用其他方法代替，可能只涉及6个点。
[0130]
总体矩阵的形成：
[0131]
在上述方程中，有未知量定义在1/2处，对位于1/2处的p和sw的取值采用上游权的原则。现在的方程中，待求量分别是(r＝1,2,3,
…
,n-1),因此共有n个待求对于组分质量守恒方程，总体矩阵的形成方式如下：
[0132]
(1)节点按照k＝1,nz；j＝1,ny；i＝1,nx的循环顺序排列；
[0133]
(2)对于给定的(i,j,k)，按照组分编号1,2，
…
,n-1(假设组分n是水组分)分别列出各组分；于是离散方程的个数是(n-1)
×
nx
×
ny
×
nz。
[0134]
按照上述方法构成的总体矩阵方程。矩阵求解采用lu分解方法。
[0135]
以下基于两个具体示例来说明本发明实施例的基于组分方程的油气运移聚集数值模拟方法的效果：
[0136]
(1)塔里木盆地小范围内的油气运移聚集模拟：
[0137]
选择塔里木盆地塔中斜坡区的一小范围进行模拟(图3中小方框内部分)。该地区有两套烃源岩，即中下寒武统和中下奥陶统。生烃模型采用behar et al.(1997)的四组分生烃模型，四组分分别为c1、c2-c5、c6-c14和c15+。油气的相态采用闪蒸计算进行模拟。通过埋藏史、热史和生排烃史的模拟，源岩中排出的四个组分的量作为源项进入到数值模拟模型中进行油气运聚模拟，模拟中的超压来自盆地模拟中埋藏史的超压模拟。图4是模拟后得到的中下寒武地层在两个不同时刻含油饱和度的分布。
[0138]
(2)整个塔里木盆地范围内的油气运移聚集模拟：
[0139]
选择整个塔里木盆地进行模拟。全盆地范围也采用中下寒武统和中下奥陶统两套烃源岩，生烃模型采用behar et al.(1997)的四组分生烃模型，四组分分别为c1、c2-c5、c6-c14和c15+。油气相态采用闪蒸计算进行模拟。通过埋藏史、热史和生排烃史的模拟，源岩中排出的四个组分的量作为源项进入到数值模拟模型中进行油气运聚模拟，模拟中的超压来自盆地模拟中埋藏史的超压模拟。图5是模拟后得到的中下寒武统地层和中下奥陶统地层在433百万年的含油饱和度的分布。
[0140]
本发明实施例提出了一种三相多组分的达西渗流的数值模拟方法。鉴于目前盆地模拟中油气组分在网格单元中通常以质量来表述的特点，本发明实施例的方法以组分质量而不是摩尔比作为基本的变量，进而进行基于有限体积法的全隐式离散求解。本发明实施例提出了一种包括了累积项、达西渗流项、压实项和源汇项的适用于盆地模拟的偏微分方程数值求解方法。本发明实施例采用有限体积法进行离散，并进行全隐式求解。求解中主要的未知量采用水的超压和烃组分的质量，其它变量和参数通过与温度、压力以及组分质量间的关系进行求解，进而达到对整个微分方程组求解的目的。
[0141]
实施例2：在实施例1提出的基于组分方程的油气运移聚集数值模拟方法的基础上，本发明实施例相应地提出了一种基于组分方程的油气运移聚集数值模拟系统，该系统包括处理器和存储器，所述处理器执行所述存储器中保存的计算机程序时实现如实施例1所述的基于组分方程的油气运移聚集数值模拟方法。
[0142]
本发明实施例的基于组分方程的油气运移聚集数值模拟系统与实施例1所述的基于组分方程的油气运移聚集数值模拟方法具有相同的有益效果，为避免重复，这里不再赘述。
[0143]
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。