CN111339702B - 油藏数值模拟等效井筒半径计算方法 - Google Patents

Abstract

本发明涉及油气藏数值模拟，特别是一种油藏数值模拟等效井筒半径计算方法，包括如下步骤：S1.根据数值离散方法对地质模型进行网格离散化，标记井筒所在的网格结点；S2.根据数值离散方法对非稳态单相渗流的压力方程进行数值离散，在网格上开展数值计算，获得当前时刻井筒所在网格结点的压力；S3.将井筒网格结点的压力代入到无限大地层非稳态渗流的解析解中，计算相应的等效井筒半径；S4.增加时间，使t_o＝t_o+△t,重复步骤S2和步骤S3,得到不同时刻的等效井筒半径。本发明基于无限大地层非稳态渗流的解析解，提出了一种计算等效井筒半径的通用方法，该方法可以针对任何的网格和任何的数值离散方法计算等效井筒半径。

Description

油藏数值模拟等效井筒半径计算方法

技术领域

本发明涉及油气藏数值模拟，特别是一种油藏数值模拟等效井筒半径计算方法。

背景技术

油藏数值模拟是利用数学方法描述油气藏中流体的渗流特征，是研究油气藏及天然气水合物藏开采宏观渗流特征、预测开采动态行为的技术手段和基本工具。

在油气藏开采数值模拟中，井筒的处理是模拟最关键的部分。由于井筒半径通常在0.1米量级，而油气藏的尺度可以达到几公里甚至几十公里，如果把真实的井筒半径考虑到模型中，则会造成网格量的急剧增加，消耗大量的计算资源。同时，对于多相流动问题，过大的网格差异会造成饱和度计算出错。因此，在实际油藏数值模拟过程中，通常要对井筒进行处理。

目前，常用的处理方法采用以下两种。第一种，将井筒半径扩大，然后将扩大后的井筒处的压力利用稳态渗流公式转换到实际的井底压力，井筒扩大后，模型的网格差异减小，网格量也相应减少。第二种，忽略井筒半径，井筒处理为源汇项，将井所在网格的压力利用等效井筒半径与实际的井底压力关联起来，该方法在模型中没有井筒，二维问题中将井筒简化为一个点，三维问题中将井筒简化为一条线，地质建模和数值模拟都十分方便，是现有油藏数值模拟中最为常用的方法。

忽略井筒半径、采用源汇项的方式处理井筒时，等效井筒半径的正确与否直接关系到整个模拟的正确性。等效井筒半径的计算与采用的数值离散方法和网格系统有关。WELL FLOW MODELS FOR VARIOUS NUMERICAL METHODS(2009,INTERNATIONAL JOURNAL OFNUMERICAL)一文中详细给出了二维问题有限差分法在正方形网格系统、有限元法在Denaulay三角网格系统和矩形网格系统以及控制体积有限元法(Control Volume BasedFinite Element Method,简写为CVFEM)的等效井筒半径的计算公式。然而，上述计算公式是二维模型的结果，对于三维模型不适用，尤其是三维非结构网格的有限元方法中，等效井筒半径的确定较为困难。

为满足油气藏数值模拟中等效井筒半径的计算，本发明公开一种基于数值模拟的计算等效井筒半径的通用方法。

发明内容

本发明的目的在于克服现有技术存在的上述缺陷，针对油气藏数值模拟中等效井筒半径在复杂网格中计算的难题，提出了一种油藏数值模拟等效井筒半径计算方法,该方法可以针对任何的网格和任何的数值离散方法计算等效井筒半径。

本发明的技术方案是：一种油藏数值模拟等效井筒半径计算方法，其中，包括如下步骤。

S1.根据数值离散方法对地质模型进行网格离散化，标记井筒所在的网格结点；

S2.根据数值离散方法对非稳态单相渗流的压力方程进行数值离散，在网格上开展数值计算，获得当前时刻井筒所在网格结点的压力；

S3.将井筒网格结点的压力代入到无限大地层非稳态渗流的解析解中，计算相应的等效井筒半径:

无限大地层非稳态渗流的解析解为：

其中，r表示距离井筒的距离,单位m；h为储层厚度，单位m；Ei为指数积分函数；

将井筒网格结点的压力p_o及对应的时刻t_o代入方程(8)中可得：

其中，r_o表示压力p_o对应的等效井筒半径；

在特定的时刻t_o，方程(9)的位置量只有r_o，采用牛顿迭代法计算r_o；

设函数

则求解r_o的牛顿迭代公式为：

其中，d为迭代步数,

表示函数

对r_o求导；使用方程(11)迭代，当|r_o ^{d +1}-r_o ^d|<10^-6时，计算满足精度；

S4.增加时间，使t_o＝t_o+Δt,重复步骤S2和步骤S3,得到不同时刻的等效井筒半径，即等效井筒半径与时间的关系,等效井筒半径将趋于一个稳定值，该稳定值即为最终的等效井筒半径。

步骤S1中，其具体处理过程包括以下步骤：

S1.1根据储层厚度和储层范围建立油藏地质模型，该模型中的井筒通过线段来表示；

S1.2根据井的坐标绘制表示井筒的直线，井筒位于储层中心，设井筒纵向上贯穿整个储层，油藏厚度为Z₀米，则井筒的线段为点(0,0,0)-(0,0,Z₀)；

S1.3用表示井筒的直线限定网格结点，网格结点均设置在表示井筒的线段上，将油藏地质模型离散，划分网格单元，网格单元为四面体；

S1.4记录表示井筒的网格结点的编号，表示井筒的网格结点的编号为按照从底部到顶部的顺序依次为0、2、3、4、5、…、a、1。

步骤S2中，其具体处理过程包括以下步骤：

S2.1根据网格系统，将井筒的流量在与井筒相关的单元上转换为源汇强度；

井筒流量Q的单位是m³/s，源汇强度q的单位是s^-1，将井筒流量转换为源汇强度的转换方法是将井筒流量在井筒单元上进行平均，井筒单元定义为有边在井筒上的单元；

相邻两网格结点之间的线段上的井筒单元数量为b个，将所有表示井筒的线段上的单元相加，得到井筒单元的数量，并设为m,第i个井筒单元的体积为V_i,则源汇强度计算公式为：

S2.2将非稳态单相渗流的压力方程在网格上离散，形成线性方程组；

非稳态单相渗流的压力方程为：

储层的外边界为封闭条件：

其中φ为储层孔隙度；c_t为储层综合压缩系数,单位Pa^-1；p为储层压力，单位Pa；k为储层渗透率，单位m²；μ为储层流体粘度,单位Pa·s；t为时间，单位s；q为源汇强度，单位s^-1，只在井所在的位置存在；n为外边界Γ的外法线方向；

采用有限元的方法对方程进行离散：计算单元的刚度方程，取网格系统中的一个单元，将压力方程(2)两端同乘以压力的变分δp并在单元上积分、降阶，可得：

其中V为单元体积,单位m³；s表示在储层外边界Γ上的单元的边；

压力和压力的变分δp有如下的插值形式：

其中，N_i为单元结点i的插值形函数，p_i为单元结点i的压力值；

代入压力的插值形式，得到单元刚度方程为：

其中，i,j,k,l分别表示四面体单元的四个顶点，

x_i,y_i,z_i为四面体顶点i的坐标；x_j,y_j,z_j为四面体顶点j的坐标；x_k,y_k,z_k为四面体顶点k的坐标；x_l,y_l,z_l为四面体顶点l的坐标；Δt为时间步长，单位s，n+1表示当前时刻，n表示前一时刻，且有Δt＝tⁿ⁺¹-tⁿ；

将所有单元的刚度方程叠加形成线性方程组

[A]{p}＝{f} (7)

其中[A]表示线性方程组的系数矩阵，由单元刚度方程(6)中的左端系数叠加而成；{p}表示所有结点的列向量；{f}表示荷载列向量，由单元刚度方程(6)右端的值叠加而成；

S2.3求解线性方程组(7)，获得所有网格结点上的压力，根据标记的井筒网格结点为0、2、3、4、5、…、a、1的压力值，该压力值表示井筒网格结点的压力p。

本发明的有益效果是：

本发明基于无限大地层非稳态渗流的解析解，提出了一种计算等效井筒半径的通用方法，该方法可以针对任何的网格和任何的数值离散方法计算等效井筒半径。

附图说明

图1为油藏地质模型示意图；

图2为油藏地质模型网格图；

图3为井筒线段和井筒单元示意图；

图4为等效井筒半径与时间关系曲线；

图5为不同渗透率条件下的等效井筒半径与时间关系曲线。

图中：1井筒；2模型外边界；3网格结点Ⅳ；4网格结点Ⅴ；5网格结点Ⅵ。

具体实施方式

为了使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。

在以下描述中阐述了具体细节以便于充分理解本发明。但是本发明能够以多种不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广。因此本发明不受下面公开的具体实施方式的限制。

本发明公开了一种油藏数值模拟等效井筒半径计算方法，该方法包括如下步骤。

第一步，根据选用的数值离散方法对地质模型进行网格离散化，标记井筒所在的网格结点。其具体处理过程包括以下步骤。

(一)根据储层厚度和储层范围建立油藏地质模型，如图1所示，该模型中的井筒1通过线段来表示；

(二)根据井的坐标绘制表示井筒的直线，井筒位于储层中心，设井筒纵向上贯穿整个储层，油藏厚度为Z₀米，则井筒的线段为点(0,0,0)-(0,0,Z₀)。本实施例中，假设油藏厚度10米，储层范围100米，则此时井筒的线段为点(0,0,0)-(0,0,10)。

(三)用表示井筒的直线限定网格结点，即网格结点必须在表示井筒1的线段上，将油藏地质模型离散，划分网格单元，网格单元为四面体，划分网格单元后的油藏地质模型网格图如图2所示。

(四)记录表示井筒的网格结点的编号，图2中表示井筒的网格结点的编号为按照从底部到顶部的顺序依次为0、2、3、4、5、…、a、1。本实施例中，表示井筒网格结点的变化为按照从底部到顶部的顺序依次为0、2、3、4、5、6、7、8、9、10、1。

第二步，根据选用的数值离散方法对非稳态单相渗流的压力方程进行数值离散，在网格上开展数值计算，获得当前时刻井筒所在网格结点的压力。其具体处理过程包括以下步骤。

(一)根据网格系统，将井筒的流量在与井筒相关的单元上转换为源汇强度。

井筒流量Q的单位是m³/s，源汇强度q的单位是s^-1，将井筒流量转换为源汇强度的转换方法是将井筒流量在井筒单元上进行平均，从而得到源汇强度，井筒单元定义为有边在井筒上的单元。

图3为选取了连接网格结点Ⅳ3和网格结点Ⅴ4的井筒线段、连接网格结点Ⅴ4和网格结点Ⅵ5的井筒线段、以及与这两条线段有关的井筒单元。网格结点Ⅳ3和网格结点Ⅴ4线段上的井筒单元数量为b个，网格结点Ⅴ4和网格结点Ⅵ5线段上的井筒单元数量为b个，本实施例中所采用的井筒单元数量为四个，因此井筒单元组成四边体。将所有表示井筒的线段上的单元相加，得到井筒单元的数量，并设为m,则第i个井筒单元的体积为V_i,则源汇强度计算公式为：

(二)将非稳态单相渗流的压力方程在网格上离散，形成线性方程组。

非稳态单相渗流的压力方程为：

储层的外边界为封闭条件：

其中φ为储层孔隙度；c_t为储层综合压缩系数,单位Pa^-1；p为储层压力，单位Pa；k为储层渗透率，单位m²；μ为储层流体粘度,单位Pa·s；t为时间，单位s；q为源汇强度，单位s^-1，只在井所在的位置存在；n为外边界Γ的外法线方向；

为了说明本方法，采用有限元的方法对方程进行离散。计算单元的刚度方程，取网格系统中的一个单元，将压力方程(2)两端同乘以压力的变分δp并在单元上积分、降阶，可得：

其中V为单元体积,单位m³；s表示在储层外边界Γ上的单元的边。

压力和压力的变分δp有如下的插值形式：

其中，N_i为单元结点i的插值形函数，p_i为单元结点i的压力值。

代入压力的插值形式，可得单元刚度方程为：

其中源汇项

只在井筒线段的点上存在，由于一个单元只可能有两个点在井筒上，故源汇项在两个点上平均分配。

其中，i,j,k,l分别表示四面体单元的四个顶点，

其中x_i,y_i,z_i为四面体顶点i的坐标；x_j,y_j,z_j为四面体顶点j的坐标；x_k,y_k,z_k为四面体顶点k的坐标；x_l,y_l,z_l为四面体顶点l的坐标；Δt为时间步长，单位s，n+1表示当前时刻，n表示前一时刻，且有Δt＝tⁿ⁺¹-tⁿ。

将所有单元的刚度方程叠加形成线性方程组

[A]{p}＝{f} (7)

其中[A]表示线性方程组的系数矩阵，由单元刚度方程(6)中的左端系数叠加而成；{p}表示所有网格结点的列向量；{f}表示荷载列向量，由单元刚度方程(6)右端的值叠加而成。

(三)求解线性方程组(7)，获得所有网格结点上的压力，根据标记的井筒网格结点获得井筒网格结点处的压力。其具体处理过程如下所述。

求解线性方程组(7)，获取所有网格结点上的压力值，得到井筒网格结点即图2中结点号为0、2、3、4、5、…、a、1的压力值，该压力值表示井筒网格结点的压力p。

第三步，将井筒网格结点的压力代入到无限大地层非稳态渗流的解析解中，计算相应的等效井筒半径。

无限大地层非稳态渗流的解析解为：

其中，r表示距离井筒的距离,单位m；h为储层厚度，单位m；Ei为指数积分函数。

将井筒网格结点的压力p_o及对应的时刻t_o代入上述方程中可得：

其中，r_o表示压力p_o对应的等效井筒半径。

在特定的时刻t_o，方程(9)的位置量只有r_o，采用牛顿迭代法计算r_o。

设函数:

则求解r_o的牛顿迭代公式为：

其中，d为迭代步数,

表示函数

对r_o求导；使用方程(11)迭代，当|r_o ^d+1-r_o ^d|<10^-6时，认为计算满足精度。

第四步，增加时间，重复第二步和第三步，得到等效井筒半径与时间的关系。等效井筒半径将趋于一个稳定值，该稳定值即为最终的等效井筒半径。其具体处理过程如下所述。

增加时间，使t_o＝t_o+Δt,重复步骤S2和步骤S3,可以得到不同时刻的等效井筒半径，如图4所示。本实施例中，从图中可以看出，随着时间推移，等效井筒半径趋于稳定值0.173m，该稳定值即为最终的等效井筒半径。

由于等效井筒半径只与离散格式和网格有关，为了验证上述方法，改变储层的渗透率，重复计算，得到如图5所示的等效井筒半径与时间的关系，可以看出不同渗透率值下，最终稳定的等效井筒半径的值都是一致的。

需要说明的是，本实施例使用了有限元方法在三维四面体网格上进行计算，但本发明的方法不限定在有限元方法和四面体网格上，对于有限差分法、有限体积法等任意离散方法及二维和三维的任意网格形式都可适用；同时，上述方法计算的等效井筒半径也可适用于多相渗流的情况。

以上对本发明所提供的油藏数值模拟等效井筒半径计算方法进行了详细介绍。本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以对本发明进行若干改进和修饰，这些改进和修饰也落入本发明权利要求的保护范围内。对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (2)

1.一种油藏数值模拟等效井筒半径计算方法，其特征在于包括如下步骤：

S1.根据数值离散方法对地质模型进行网格离散化，标记井筒所在的网格结点；

S2.根据数值离散方法对非稳态单相渗流的压力方程进行数值离散，在网格上开展数值计算，获得当前时刻井筒所在网格结点的压力；

S3.将井筒网格结点的压力代入到无限大地层非稳态渗流的解析解中，计算相应的等效井筒半径:

无限大地层非稳态渗流的解析解为：

其中，r表示距离井筒的距离,单位m；h为储层厚度，单位m；Ei为指数积分函数；

将井筒网格结点的压力p_o及对应的时刻t_o代入方程(8)中得到：

其中，r_o表示压力p_o对应的等效井筒半径；

在特定的时刻t_o，方程(9)的位置量只有r_o，采用牛顿迭代法计算r_o；

设函数

则求解r_o的牛顿迭代公式为：

其中，d为迭代步数,

表示函数

对r_o求导；使用方程(11)迭代，当|r_o ^d+1-r_o ^d|<10^-6时，计算满足精度；

S4.增加时间，使t_o＝t_o+△t,重复步骤S2和步骤S3,得到不同时刻的等效井筒半径，即等效井筒半径与时间的关系,等效井筒半径将趋于一个稳定值，该稳定值为最终的等效井筒半径；

步骤S2中，其具体处理过程包括以下步骤：

S2.1根据网格系统，将井筒的流量在与井筒相关的单元上转换为源汇强度；

井筒流量Q的单位是m³/s，源汇强度q的单位是s^-1，将井筒流量转换为源汇强度的转换方法是将井筒流量在井筒单元上进行平均，井筒单元定义为有边在井筒上的单元；

相邻两网格结点之间的线段上的井筒单元数量为b个，井筒单元组成四边体，将所有表示井筒的线段上的单元相加，得到井筒单元的数量m,第i个井筒单元的体积为V_i,则源汇强度计算公式为：

S2.2将非稳态单相渗流的压力方程在网格上离散，形成线性方程组；

非稳态单相渗流的压力方程为：

储层的外边界为封闭条件：

其中φ为储层孔隙度；c_t为储层综合压缩系数,单位Pa^-1；p为储层压力，单位Pa；k为储层渗透率，单位m²；μ为储层流体粘度,单位Pa·s；t为时间，单位s；q为源汇强度，单位s^-1，只在井所在的位置存在；n为外边界Γ的外法线方向；

采用有限元的方法对方程进行离散：计算单元的刚度方程，取网格系统中的一个单元，将压力方程(2)两端同乘以压力的变分δp并在单元上积分、降阶，得到：

其中V为单元体积,单位m³；s表示在储层外边界Γ上的单元的边；

压力和压力的变分δp有如下的插值形式：

其中，N_i为单元结点i的插值形函数，p_i为单元结点i的压力值；

代入压力的插值形式，得到单元刚度方程为：

其中，i,j,k,l分别表示四面体单元的四个顶点，

x_i,y_i,z_i为四面体顶点i的坐标；x_j,y_j,z_j为四面体顶点j的坐标；x_k,y_k,z_k为四面体顶点k的坐标；x_l,y_l,z_l为四面体顶点l的坐标；△t为时间步长，单位s，n+1表示当前时刻，n表示前一时刻，且有△t＝tⁿ⁺¹-tⁿ；

将所有单元的刚度方程叠加形成线性方程组

[A]{p}＝{f} (7)

其中[A]表示线性方程组的系数矩阵，由单元刚度方程(6)中的左端系数叠加而成；{p}表示所有结点的列向量；{f}表示荷载列向量，由单元刚度方程(6)右端的值叠加而成；

S2.3求解线性方程组(7)，获得所有网格结点上的压力，根据标记的井筒网格结点为0、2、3、4、5、…、a、1的压力值，该压力值表示井筒网格结点的压力p。

2.根据权利要求1所述的油藏数值模拟等效井筒半径计算方法，其特征在于步骤S1中，其具体处理过程包括以下步骤：

S1.1根据储层厚度和储层范围建立油藏地质模型，该模型中的井筒通过线段来表示；

S1.2根据井的坐标绘制表示井筒的直线，井筒位于储层中心，设井筒纵向上打开全部储层，油藏厚度为Z₀米，则井筒的线段为点(0,0,0)-(0,0,Z₀)；

S1.3用表示井筒的直线限定网格结点，网格结点均设置在表示井筒的线段上，将油藏地质模型离散，划分网格单元，网格单元为四面体；

S1.4记录表示井筒的网格结点的编号，表示井筒的网格结点的编号为按照从底部到顶部的顺序依次为0、2、3、4、5、…、a、1。

Images