CN108241777A - 基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法 - Google Patents

Abstract

本发明涉及海洋天然气水合物资源开采数值模拟技术领域，具体涉及基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，包括S1.根据地质资料，建立水合物降压开采储层的地质模型，包括井筒内边界和储层外边界；S2.在地质模型上建立非结构网格系统；S3.建立水合物降压开采的多场耦合数学模型，包括气相压力方程、水相饱和度方程、水合物分解方程、初始条件、内边界条件和外边界条件；S4.在建立的非结构网格系统上，利用Galerkin有限元方法建立单元刚度方程，对多场耦合数学模型进行离散求解，获得压力场和饱和度场的空间分布；S5.根据求解得到的网格结点上的压力场和饱和度场，利用单元刚度方程计算气相和水相的速度场。该方法提高了计算速度和精度。

Description

基于非结构网格有限元法计算水合物沉积物中渗流速度场的 方法

技术领域

本发明涉及海洋天然气水合物资源开采数值模拟技术领域，具体涉及一种基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法。

背景技术

天然气水合物分布范围广，资源量巨大，能量密度高，是有望成为满足未来人类能源需求的高效清洁能源。降压开采时，水合物储层中的气水渗流规律是制定合理降压方案和工作制度的重要依据。水合物开采过程是一个复杂的传热传质过程，涉及到水合物的分解、气水两相渗流和热传递过程，各物理过程间相互耦合，因此，数值模拟技术是计算水合物储层内气水渗流过程，制定降压方案，优化开发技术指标的关键技术。

数值模拟技术是将计算区域进行网格剖分，在离散的网格系统上进行数值计算，获得各网格结点上的物理量的值，形成物理量的空间演化特征。如何在离散的网格系统上利用离散点的值计算速度矢量场是数值模拟技术面临的一项技术难题。在数值模拟技术中，通常使用两种方式来计算渗流速度场。一种是在求出压力场后，直接利用达西定律，基于不同方向的离散点的压力梯度来计算速度场；另外一种是将达西定律的方程作为求解的基本方程，与质量连续方程形成方程组，速度作为基本的求解变量直接计算速度场。前一种方法是最为直接的方法，计算量较小，但该方法主要适用于规则结构网格系统的计算，而对于复杂的井型和储层，结构网格系统不适用；后一种方法计算精度高，但由于将速度作为未知量，大大增加了求解的未知量的数量，计算过程需要求解大型的线性方程组，资源消耗大大提升，计算速度急剧降低。

因此，为满足水合物开采数值模拟快速高效的需求，计算降压开采过程中储层的渗流速度场，为认识水合物开采的传热传质规律提供技术手段，本发明公开一种基于非结构网格有限元法计算水合物降压开采渗流速度场的方法。

发明内容

本发明针对目前的水合物开采渗流速度场计算方法无法适用于非结构网格以及计算速度慢的问题，提出一种基于非结构网格有限元法的水合物沉积物中渗流速度的计算方法，该方法可以在非结构网格上直接利用结点压力计算渗流速度场，而无需求解线性方程组，可快速高效的获得水合物沉积物中的渗流速度场，解决数值模拟技术中速度场求解的难题。

为了实现上述目的，本发明采用如下技术方案：基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，包括如下步骤，

S1.根据地质资料，建立水合物降压开采储层的地质模型，包括井筒内边界和储层外边界；

S2.在地质模型上建立非结构网格系统；

S3.建立水合物降压开采的多场耦合数学模型，包括气相压力方程、水相饱和度方程、水合物分解方程、初始条件、内边界条件和外边界条件；

S4.在建立的非结构网格系统上，利用Galerkin有限元方法建立单元刚度方程，对多场耦合数学模型进行离散求解，获得压力场和饱和度场的空间分布；

S5.根据求解得到的网格结点上的压力场和饱和度场，利用单元刚度方程计算气相和水相的速度场。

进一步地，所述步骤2的具体包括：

S21.将地质模型的井筒内边界和储层外边界的几何信息输入；

S22.在井筒内边界和储层外边界上设置取样点，形成边界的离散点；

S23.以边界离散点为基础，对整个模型进行网格离散化处理得到非结构网格系统。

进一步地，所述步骤S2中，离散化之前，先将井筒附近网格进行加密处理。

根据权利要求3所述的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，其特征在于，在井筒附近网格进行加密处理具体为：设置地质模型的井筒内边界上的取样点的间距小于储层外边界上的取样点的间距。

进一步地，所述步骤S4具体包括：

S41.取非结构网格系统上的某一单元，计算单元刚度方程；

S42.计算总体刚度方程；

S43.根据边界条件求解总体刚度方程。

进一步地，所述步骤S5具体包括：

S51.对单元刚度方程进行变形；

S52.计算单元边界上的速度；

S53.单元边界上速度的分解；

S54.单元边界速度的合成；

S55.绘制速度场。

本发明的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，克服了传统有限元计算压力场和饱和度场中的惯性思维，基于单元刚度方程中包含的速度信息，在单元上利用结点压力值和饱和度值直接计算速度场，而不需要求解大型线性方程组，大大提高了计算速度，且计算精度高。

附图说明

图1为水合物降压开采储层地质模型图；

其中，1-地质模型井筒内边界；2-储层外边界；

图2为水合物降压开采储层模型网格图；

图3为三角形单元边界上的速度示意图；

图4为三角形单元边界速度分解示意图；

图5为气相速度场图；

图6为水相速度场图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本发明适用于基于有限元的二维和三维非结构网格，如二维的三角形网格、三维的四面体网格，下面以二维三角形网格为例对本发明的方法进行详细的说明：

本发明的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，包括如下步骤：

S1.根据地质资料，建立水合物降压开采储层的地质模型，包括井筒内边界和储层外边界；如图1所示，井筒半径为0.1米，储层圆形半径为10米。

S2.在地质模型上建立非结构网格系统；

现有技术中，非结构网格系统的建立以及离散化的实现有多种方法，在本实施例中，非结构网格系统的建立以及离散化的实现是借助Netgen开源的网格划分程序实现的。

具体处理的步骤为：

S21.将地质模型的井筒内边界和储层外边界的几何信息输入；

S22.在井筒内边界和储层外边界上设置取样点，形成边界的离散点；

S23.以边界离散点为基础，对整个模型进行网格离散化处理得到非结构网格系统。

为了确保计算精度，离散化之前，先将井筒附近网格进行加密处理，在Netgen开源的网格划分程序中实现的具体步骤为：设置地质模型的井筒内边界上的取样点的间距小于储层外边界上的取样点的间距。

如图2所示，由Netgen开源的网格划分程序处理后，获取的三角形网格单元数量为1202，结点数为619。

S3.建立水合物降压开采的多场耦合数学模型，包括气相压力方程、水相饱和度方程、水合物分解方程、初始条件、内边界条件和外边界条件；具体为：

(1)气相压力方程：

(2)水相饱和度方程：

(3)水合物分解方程：

(4)初始条件：

p_g(t＝0)＝p_i，S_w(t＝0)＝S_wi，T(t＝0)＝T_i (5)

(5)内边界条件：

(6)外边界条件：

式中符号含义如下：

ρ_g和ρ_w分别为气相和水相密度，kg/m³；S_g和S_w分别为气相和水相饱和度；k_rg和k_rw分别为气相水两相渗流的气相相对渗透率和水相相对渗透率；k为沉积物多孔介质的绝对渗透率，m²；μ_g和μ_w分别为甲烷气体和水的黏度，Pa·s；Z为气体偏差因子；R_g为气体常数，8.314J/(mol·K)；T为温度，K；φ为含水合物沉积物孔隙度；ρ_h为水合物密度，kg/m³；S_h为水合物饱和度；m_h为单位体积水合物分解速率，kg/(m³·s)；m_g为单位体积水合物分解的产气速率，kg/(m³·s)；m_w为单位体积水合物分解产水速率，kg/(m³·s)；为甲烷的摩尔质量，

0.016kg/mol；p_e为水合物平衡压力，Pa；p_g为气体压力，Pa；p_c为毛管力，Pa；k_reac为水合物分解速率常数，mol/(m²·Pa·s)；A_rs为单位体积储层水合物分解的表面积，m^-1；p_i为储层初始压力，Pa；S_wi为储层初始含水饱和度；T_i为储层初始温度，K；p_wf为井底压力，Pa；Γ_i表示井筒边界；Γ_o表示外边界；n_o表示外边界的法线方向；

S4.在建立的非结构网格系统上，利用Galerkin有限元方法建立单元刚度方程，对多场耦合数学模型进行离散求解，获得压力场和饱和度场的空间分布；

获得压力场和饱和度场的空间分布，即得到网格结点上的气相压力值和水相饱和度值。所述步骤S4具体包括：

S41.取非结构网格系统上的某一单元，计算单元刚度方程；

S42.计算总体刚度方程；

S43.根据边界条件求解总体刚度方程。

其中，气相压力值的具体求解过程如下：

(a)取非结构网格系统上的某一单元，计算气相压力的单元刚度方程；

取图2所示的网格系统中的某一个单元(图3)，将气相压力方程(1)两端同乘以压力变分并在三角形单元上积分、降阶：

其中，δp_g为气相压力的变分。

代入压力的插值函数计算刚度方程，压力在单元上的插值为：

式中，N为单元形函数，s为图3所示三角形单元的边界；为单元边界法线方向的压力梯度。

最终单元刚度方程为：

其中，

在图3所示的三角形单元上，代入插值形函数，可得气相压力的单元刚度方程为：

其中，L_i为图3所示三角形单元的jk边的边长。

(b)计算气相压力的总体刚度方程；

将图3所示的所有网格的单元刚度方程按照单元结点编号进行叠加。方程(11)中的最后两项为单元边界上的流量。需要说的是对于结点i来说，和这两项分别表示表示图3所示单元边界ij和ik上的流量；对于结点j，则表示单元边界ji和jk上的流量；对于结点k，则表示单元边界ki、kj上的流量。叠加过程中，方程(11)的最后两项在相邻单元的公共边上相互抵消。

(c)根据边界条件求解气相压力的总体刚度方程；

叠加后形成总体刚度方程的线性方程组，结合相应的边界条件求解线性方程组获得网格结点的气相压力值。

其中，水相饱和度值的具体求解过程如下：

(a)取非结构网格系统上的某一单元，计算水相饱和度的单元刚度方程；

取图2所示的网格系统中的某一个单元(图3)，将水相饱和度方程(2)两端同乘以饱和度变分并在三角形单元上积分、降阶：

代入饱和度的插值函数计算单元刚度方程，饱和度在单元上的插值为：

S_w＝N{S_w}^e，δS_w＝N{δS_w}^e (13)

最终饱和度的单元刚度方程为：

上式中，

在图3所示的三角形单元上，代入插值形函数，可得水相饱和度的单元刚度方程为：

(b)计算水相饱和度的总体刚度方程；

将图2所示的所有网格的单元刚度方程按照单元结点编号进行叠加，形成最终的总体刚度方程。方程(15)中的最后两项为单元边界上的流量，叠加过程中，方程(15)的最后两项在相邻单元的公共边上相互抵消。

(c)根据边界条件求解水相饱和度的总体刚度方程；

结合相应的边界条件求解总体刚度方程的线性方程组获得网格结点的水相饱和度值。

S5.根据求解得到的网格结点上的压力场和饱和度场，利用单元刚度方程计算气相和水相的速度场。具体步骤如下：

S51.对单元刚度方程进行变形；

S52.计算单元边界上的速度；

S53.单元边界上速度的分解；

S54.单元边界速度的合成；

S55.绘制速度场。

其中，气相速度场的计算，主要包括以下步骤：

(a)对气相压力单元刚度方程进行变形；

通过第一步至第四步，图2所示的网格结点上的气相压力值为已知量。取图3所示的单元分析，其单元刚度方程为方程(11)，其中方程(11)的最后一项为图3所示单元边界上的气相流量(气相速度与边界长度的乘积)，该项在单元刚度方程叠加过程中被抵消，而没有求解。在图3的单元上所示将方程(11)中的最后一项单独提出，做如下变形：

(b)计算单元边界上的气相速度

变形后的方程(16)右端项全部为已知量，将左端的气相速度作为未知量，将方程(16)右端项用M_i表示，则在图3所示的三角形单元上有如下方程：

求解方程(17)，可得图3所示的三角形单元三条边上的速度值：

对图2所示的所有单元采用公式(18)进行计算，即可得到所有三角形单元三条边上的速度值，速度为正的方向定义为三角形边的外法线方向。

(c)单元边界上气相速度的分解

将三角形边上的速度值分解为x和y两个方向。以ij边为例，说明其分解方法。如图4所示，做出ij边对应的三角形结点k关于ij边的对称点s，以及ks连线在ij上的交点p，则ij边上的速度v_gk分解为x和y方向为：

式中，v_gkx为ij边x方向的气相速度，v_gky为ij边y方向气相速度，x_s，y_s和x_p，y_p分别为s点和p点的坐标。

(d)单元边界气相速度的合成

对于具有公共边的两个三角形单元，利用两个三角形的单元刚度方程可在公共边上计算得到两个速度值，将两个速度值的平均值作为该公共边上的速度值。

(e)绘制气相速度场

最终利用tecplot软件做出气相速度场如图5所示。

水相速度场的计算与气相速度的计算步骤和方法完全一致，不同在于(a)步使用的是水相饱和度的单元刚度方程(15)，将方程(15)变形如下：

将方程(20)的右端项设为M_i，则水相速度的计算步骤与气相速度计算步骤的(b)、(c)、(d)、(e)完全一致，最终用tecplot做出水相速度场如图6所示。

本发明的核心是基于三角形非结构网格的有限元方法来计算渗流速度场。主要的步骤是在有限元求解得到三角形结点上的压力值和饱和度值以后，将单元刚度方程变形，利用单元刚度方程中包含的边界速度信息计算得到速度场。

而本领域中，在常规有限元计算时，积分降阶的方程为：

为了更加方便的处理边界条件，在计算单元刚度方程之前，方程中的右端边界项是作为整个求解区域的整体考虑，直接利用模型的边界条件进行处理，即对于本质边界条件，这一项中的δp_g＝0，而对于自然边界条件，直接将代入到方程中。因此，在常规有限元计算时，形成总体刚度方程后，右端的边界项不直接出现在方程中，而是用求解区域的边界条件代替。

与常规做法不同的是，本发明是在一个三角形网格单元上进行计算，将边界项作为单元的边界项，而不是整个求解区域的边界，充分利用单元刚度方程中的右端项的中包含的速度信息，将这一项作为未知量计算，获得了每个单元上的速度，从而计算得到整个区域的速度场，克服了本领域技术人员的惯性思维，实现了不用求解线性方程组而计算速度场的目的，大大提高了计算速度，且计算精度高。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (6)

1.基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，其特征在于，包括如下步骤，

S1.根据地质资料，建立水合物降压开采储层的地质模型，包括井筒内边界和储层外边界；

S2.在地质模型上建立非结构网格系统；

S3.建立水合物降压开采的多场耦合数学模型，包括气相压力方程、水相饱和度方程、水合物分解方程、初始条件、内边界条件和外边界条件；

S4.在建立的非结构网格系统上，利用Galerkin有限元方法建立单元刚度方程，对多场耦合数学模型进行离散求解，获得压力场和饱和度场的空间分布；

S5.根据求解得到的网格结点上的压力场和饱和度场，利用单元刚度方程计算气相和水相的速度场。

2.根据权利要求1所述的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，其特征在于，所述步骤2的具体包括：

S21.将地质模型的井筒内边界和储层外边界的几何信息输入；

S22.在井筒内边界和储层外边界上设置取样点，形成边界的离散点；

S23.以边界离散点为基础，对整个模型进行网格离散化处理得到非结构网格系统。

3.根据权利要求1或2所述的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，其特征在于，所述步骤S2中，离散化之前，先将井筒附近网格进行加密处理。

4.根据权利要求3所述的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，其特征在于，在井筒附近网格进行加密处理具体为：设置地质模型的井筒内边界上的取样点的间距小于储层外边界上的取样点的间距。

5.根据权利要求4所述的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，其特征在于，所述步骤S4具体包括：

S41.取非结构网格系统上的某一单元，计算单元刚度方程；

S42.计算总体刚度方程；

S43.根据边界条件求解总体刚度方程。

6.根据权利要求5所述的基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法，其特征在于，所述步骤S5具体包括：

S51.对单元刚度方程进行变形；

S52.计算单元边界上的速度；

S53.单元边界上速度的分解；

S54.单元边界速度的合成；

S55.绘制速度场。