CN109670220A - 一种基于非结构网格的水平井气水两相数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于非结构网格的水平井气水两相数值模拟方法，包括以下步骤：S1：低渗透裂缝性储层水平井的三维地质体生成及网格离散；S2：裂缝性致密储层多相流渗流模型建立；S3：非结构网格单元特征分析和整体数值计算格式建立；S4：考虑水平井筒内流动特征的水平井井点网格拟稳态井模型；S5：对建立的地层渗流与井筒内流动耦合模型进行求解，并对含水饱和度、应力敏感等参数敏感性进行分析。本发明考虑水平井井筒流动情况和低渗透储层微小裂缝发育特征，采用有限控制体积方法，基于非结构网格，获取低渗透裂缝性储层气水两相流动数值解，解决了传统结构网格模拟器出现的网格取向性强、忽略应力敏感及水平井筒内流动特征等问题。

Description

一种基于非结构网格的水平井气水两相数值模拟方法

技术领域

本发明涉及非常规低渗透储层的勘探与开发领域，尤其涉及一种基于非结构网格的水平井气水两相数值模拟方法。

背景技术

随着常规高渗储层开采殆尽，为满足日益增长的能源需求，对于非常规低渗透储层的勘探与开发成为热点和难点。矿场实践表明：对于低渗透致密储层，由于渗透率低，储层致密，有效厚度低，采用水平井技术，扩大泄油面积，成为低渗透储层高效生产的重要手段。同时，处于开发中后期的低渗透气藏，不同程度地见水使井筒流动和地层渗流呈现多相流情况，对数值模拟结果的影响日益突出。油气藏数值模拟技术作为目前定量描述非均质储层多场多相流动规律的重要手段，在开发方案设计，井网部署优化调整，采收率预测等方面有着不可替代的作用。目前，随着计算机硬件的不断发展，以及非常规油气藏的开发，采用非结构网格对复杂流动边界进行精细描述，并结合有限元、有限体积等数值离散方法进行数值求解，已成为油气藏数值模拟领域的主要研究和发展方向，然而传统结构网格模拟器具有网格取向性强、忽略应力敏感及水平井筒内流动特征等问题。

发明内容

本发明考虑水平井井筒流动情况和低渗透储层微小裂缝发育特征，采用有限控制体积方法，基于非结构网格，获取低渗透裂缝性储层气水两相流动数值解，力图解决传统结构网格模拟器出现的网格取向性强、忽略应力敏感及水平井筒内流动特征等问题，提出了一种基于非结构网格的水平井气水两相数值模拟方法。

本发明包括以下步骤：

S1：低渗透裂缝性储层水平井的三维地质体生成及网格离散；

S2：裂缝性致密储层多相流渗流模型建立；

S3：非结构网格单元特征分析和整体数值计算格式建立；

S4：考虑水平井筒内流动特征的水平井井点网格拟稳态井模型；

S5：对建立的地层渗流与井筒内流动耦合模型进行求解，并对含水饱和度、应力敏感等参数敏感性进行分析。

进一步，所述步骤S1还包括以下步骤：

S1.1：根据研究工区的实际地质情况，确定工区边界生成地质体；

S1.2：按照开源非结构网格剖分软件Distmesh的数据结构要求，对地质体数据进行编辑导入，生成网格离散模型。

进一步，所述步骤S2还包括基质***：

考虑应力敏感效应的裂缝***：

k_fe＝k_fie^-γ(p_fi-p_f)

补充方程：

S_w+S_g＝1

P_w＝P_g-P_c(S_w)

P_c＝-B_mln(S_w) (3)

边界条件：

P_c|Γ_out＝P_i (4)

式中符号定义：k_m为基质渗透率，mD；φ_m为基质孔隙度；k_fi为裂缝***原始渗透率，mD；k_fe为应力敏感下裂缝修正渗透率，mD；φ_f为裂缝孔隙度；α为形状因子,1/m²；γ为渗透率模量，MPa^-1；q_sct为标况产量，万方/天；p_w为水相压力，MPa；p_g为气相压力， MPa；p_c为毛细管力，MPa；B_m为毛管力修正系数；p_e为外边界压力，MPa；p_i为原始地层压力，MPa；k_r(·)、μ_(·)、B_(·)、S_(·)分别为相对渗透率、粘度、体积系数和饱和度，包括了基质气相和水相，裂缝气相和水相。

进一步，所述步骤S3还包括以下步骤：

S3.1：选取一个四面体单元，建立单元特征矩阵，定义单元内压力和饱和度的试探解为：

式中符号定义：单元网格中压力试探解；为单元网格中饱和度试探解；N_v为插值函数；v＝1，2，3，4代表四面体的四个顶点。

式中符号定义：V为四面体体积；a，b，c，d为系数。

下一时步的饱和度值：

Pⁿ⁺¹≈P^k+1＝P^k+δP^k (7)

式中符号定义：δ为算子，表示第k次迭代到k+1次迭代后的变化；n为上一时步；n+1为下一时步。

基质单元矩阵为：

裂缝单元矩阵为：

式中符号定义：P、T、W和M分别为压力矩阵、传导率矩阵、窜流项矩阵和时间项矩阵；δP、δT、δW和δM则为相应的算子矩阵；△t为时间步长。将压力和饱和度的试探解带入渗流方程中，可推导出P、T、W和M的具体计算格式为：

(基质)

(裂缝)

(基质)，(裂缝)

S3.2：对离散区域内的每一个四面体网格，建立类似的单元矩阵。假设存在N个网格节点，通过整体叠加，形成关于基质和裂缝***的4N×4N大矩阵：

[K]_4N*4N[δX]_4N*1＝[R]_4N*1 (15)

式中符号定义：K为整体系数矩阵；δX为未知变量变化量矩阵，包括δP_m，δP_f，δS_wm和δS_wf；R为余量项。

进一步，所述步骤S4还包括：

由Peaceman公式，水平井总气水产量可表示为：

式中符号定义：N_w为水平井穿过的离散网格节点；PI为生产指数；p_well为井点节点处水平井筒压力；p_ave为井点节点所在网格的平均压力；

考虑水平井筒内流动特征，采用水平井筒多相流计算公式获得考虑井筒摩阻和流体动量损失的井筒压力沿程梯度：

式中符号含义：ρ_l为液体密度，kg/m³；ρ_g为气体密度，kg/m³；G为气液混合物质量流量，kg/s；v_m为混合物流动速度，m/s；v_sg为气体表观流速，m/s；A为井筒油管截面积， m²；D为油管内径，m。

持液率H_L和摩擦阻力系数f采用Beggs-Brill方法计算。水平井跟端网格点井筒压力 p_well1等于井底压力值p_bh，则结合式(17)可以得到各井点网格处井筒压力值：

将得到的井筒压力值p_welli带入到式(16)中，结合井点网格压力p_avei，得到每时步下的产量值。

进一步，所述步骤S5还包括以下步骤:

S5.1：由有限元基本原理，对于定压外边界条件，通过置大数法，将定值(p_e)向量添加到式(15)载荷矩阵K中；对于封闭外边界条件，添加0向量到式(15)载荷矩阵K中；

S5.2：将式(16)中下一时步总产量转化为并带入到式(15)中，得到考虑内外边界条件的数值计算稀疏矩阵；

S5.3：采用Orthomin共轭梯度法对稀疏矩阵进行求解，通过牛顿-拉夫逊迭代方法，获得基质和裂缝***在一个时步下的压力和饱和度变化值：δp_gf1,…δp_gfN；δp_gm1,…δp_gmN；δS_wf1,…δS_wfN；δS_wm1,…δS_wmN。进而由式(7)和式(8)获得下一时步的值；

S5.4：输出各时步下的离散网格节点压力、饱和度值，以及水平井井底流压和产量值。

本发明的有益效果在于：考虑水平井井筒流动情况和低渗透储层微小裂缝发育特征，采用有限控制体积方法，基于非结构网格，获取低渗透裂缝性储层气水两相流动数值解，解决了传统结构网格模拟器出现的网格取向性强、忽略应力敏感及水平井筒内流动特征等问题。

附图说明

图1是耦合模型建立及数值计算流程图；

图2是矩形封闭气藏水平井三维地质体及网格离散；

图3是四面体网格单元示意图；

图4是不同含水饱和度下模拟结果与商业软件对比图；

图5是应力敏感效应对产量的影响。

具体实施方式

为了使本领域的技术人员更好地理解本发明的技术方案，下面结合附图和具体实施例对本发明作进一步的详细说明。

如图1所示，本发明的主要步骤包括：

S1：低渗透裂缝性储层水平井的三维地质体生成及网格离散；

S2：裂缝性致密储层多相流渗流模型建立；

S3：非结构网格单元特征分析和整体数值计算格式建立；

S4：考虑水平井筒内流动特征的水平井井点网格拟稳态井模型；

S5：对建立的地层渗流与井筒内流动耦合模型进行求解，并对含水饱和度、应力敏感等参数敏感性进行分析。

进一步，所述步骤S1还包括以下步骤：

S1.1：根据研究工区的实际地质情况，确定工区边界生成地质体；

S1.2：按照开源非结构网格剖分软件Distmesh的数据结构要求，对地质体数据进行编辑导入，生成网格离散模型，如图2所示。

进一步，所述步骤S2还包括基质***：

考虑应力敏感效应的裂缝***：

k_fe＝k_fie^-γ(P_fi-P_f)

补充方程：

S_w+S_g＝1

P_w＝P_g-P_c(S_w)

P_c＝-B_ml_n(S_w) (3)

边界条件：

P_c|Γ_out＝P_i (4)

式中符号定义：k_m为基质渗透率，mD；φ_m为基质孔隙度；k_fi为裂缝***原始渗透率，mD；k_fe为应力敏感下裂缝修正渗透率，mD；φ_f为裂缝孔隙度；α为形状因子,1/m²；γ为渗透率模量，MPa^-1；q_sct为标况产量，万方/天；p_w为水相压力，MPa；p_g为气相压力， MPa；p_c为毛细管力，MPa；B_m为毛管力修正系数；p_e为外边界压力，MPa；p_i为原始地层压力，MPa；k_r(·)、μ_(·)、B_(·)、S_(·)分别为相对渗透率、粘度、体积系数和饱和度，包括了基质气相和水相，裂缝气相和水相。

进一步，所述步骤S3对于步骤2建立的非线性渗流方程组，采用有限控制体积方法结合非结构四面体网格建立数值计算格式。具体步骤为：

S3.1：选取一个四面体单元，如图3所示，建立单元特征矩阵，定义单元内压力和饱和度的试探解为：

式中符号定义：单元网格中压力试探解；为单元网格中饱和度试探解；N_v为插值函数；v＝1，2，3，4代表四面体的四个顶点。

式中符号定义：V为四面体体积；a，b，c，d为系数。

采用全隐式计算格式，则下一时步的饱和度值可表示为：

Pⁿ⁺¹≈P^k+1＝P^k+δP^k (7)

式中符号定义：δ为算子，表示第k次迭代到k+1次迭代后的变化；n为上一时步；n+1为下一时步。

基质单元矩阵为：

裂缝单元矩阵为：

式中符号定义：P、T、W和M分别为压力矩阵、传导率矩阵、窜流项矩阵和时间项矩阵；δP、δT、δW和δM则为相应的算子矩阵；△t为时间步长。将压力和饱和度的试探解带入渗流方程中，可推导出P、T、W和M的具体计算格式为：

(基质)

(裂缝)

(基质)，(裂缝)

S3.2：对离散区域内的每一个四面体网格，建立类似的单元矩阵。假设存在N个网格节点，通过整体叠加，形成关于基质和裂缝***的4N×4N大矩阵：

[K]_4N*4N[δX]_4N*1＝[R]_4N*1 (15)

式中符号定义：K为整体系数矩阵；δX为未知变量变化量矩阵，包括δP_m，δP_f，δS_wm和δS_wf；R为余量项。

进一步，所述步骤S4还包括：

由Peaceman公式，水平井总气水产量可表示为：

式中符号定义：N_w为水平井穿过的离散网格节点；PI为生产指数；p_well为井点节点处水平井筒压力；p_ave为井点节点所在网格的平均压力；

考虑水平井筒内流动特征，采用水平井筒多相流计算公式获得考虑井筒摩阻和流体动量损失的井筒压力沿程梯度：

式中符号含义：ρ_l为液体密度，kg/m³；ρ_g为气体密度，kg/m³；G为气液混合物质量流量，kg/s；v_m为混合物流动速度，m/s；v_sg为气体表观流速，m/s；A为井筒油管截面积， m²；D为油管内径，m。

持液率H_L和摩擦阻力系数f采用Beggs-Brill方法计算。水平井跟端网格点井筒压力 p_well1等于井底压力值p_bh，则结合式(17)可以得到各井点网格处井筒压力值：

将得到的井筒压力值p_welli带入到式(16)中，结合井点网格压力p_avei，得到每时步下的产量值。

进一步，所述步骤S5还包括以下步骤:

S5.1：由有限元基本原理，对于定压外边界条件，通过置大数法，将定值(p_e)向量添加到式(15)载荷矩阵K中；对于封闭外边界条件，添加0向量到式(15)载荷矩阵K中；

S5.2：将式(16)中下一时步总产量转化为并带入到式(15)中，得到考虑内外边界条件的数值计算稀疏矩阵；

S5.3：采用Orthomin共轭梯度法对稀疏矩阵进行求解，通过牛顿-拉夫逊迭代方法，获得基质和裂缝***在一个时步下的压力和饱和度变化值：δp_gf1,…δp_gfN；δp_gm1,…δp_gmN；δS_wf1,…δS_wfN；δS_wm1,…δS_wmN。进而由式(7)和式(8)获得下一时步的值；

S5.4：输出各时步下的离散网格节点压力、饱和度值，以及水平井井底流压和产量值。

考虑多相流动情况下的数值计算结果同商业软件模拟结果吻合良好(图4)，说明本方法正确有效。应力敏感系数对产气量的影响如图5所示，在低渗透裂缝性气藏中，裂缝应力敏感效应对水平井产量影响显著，不能忽略。

以上所揭露的仅为本发明较佳实施例而已，当然不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。

Claims (6)

1.一种基于非结构网格的水平井气水两相数值模拟方法，其特征在于，包括以下步骤：

S1：低渗透裂缝性储层水平井的三维地质体生成及网格离散；

S2：裂缝性致密储层多相流渗流模型建立；

S3：非结构网格单元特征分析和整体数值计算格式建立；

S4：考虑水平井筒内流动特征的水平井井点网格拟稳态井模型；

S5：对建立的地层渗流与井筒内流动耦合模型进行求解，并对含水饱和度、应力敏感等参数敏感性进行分析。

2.根据权利要求1所述的一种基于非结构网格的水平井气水两相数值模拟方法，其特征在于：所述步骤S1还包括以下步骤：

S1.1：根据研究工区的实际地质情况，确定工区边界生成地质体；

S1.2：按照开源非结构网格剖分软件Distmesh的数据结构要求，对地质体数据进行编辑导入，生成网格离散模型。

3.根据权利要求1所述的一种基于非结构网格的水平井气水两相数值模拟方法，其特征在于，所述步骤S2还包括基质***：

考虑应力敏感效应的裂缝***：

补充方程：

S_w+S_g＝1

P_w＝P_g-P_c(S_w)

P_c＝-B_mln(S_w) (3)

边界条件：

P_c|Γ_out＝P_i (4)

式中符号定义：k_m为基质渗透率，mD；φ_m为基质孔隙度；k_fi为裂缝***原始渗透率，mD；k_fe为应力敏感下裂缝修正渗透率，mD；φ_f为裂缝孔隙度；α为形状因子,1/m²；γ为渗透率模量，MPa^-1；q_sct为标况产量，万方/天；p_w为水相压力，MPa；p_g为气相压力，MPa；p_c为毛细管力，MPa；B_m为毛管力修正系数；p_e为外边界压力，MPa；p_i为原始地层压力，MPa；k_r(·)、μ_(·)、B_(·)、S_(·)分别为相对渗透率、粘度、体积系数和饱和度，包括了基质气相和水相，裂缝气相和水相。

4.根据权利要求1所述的一种基于非结构网格的水平井气水两相数值模拟方法，其特征在于：所述步骤S3还包括以下步骤：

S3.1：选取一个四面体单元，建立单元特征矩阵，定义单元内压力和饱和度的试探解为：

式中符号定义：单元网格中压力试探解；为单元网格中饱和度试探解；N_v为插值函数；v＝1，2，3，4代表四面体的四个顶点。

式中符号定义：V为四面体体积；a，b，c，d为系数。

下一时步的饱和度值：

Pⁿ⁺¹≈p^k+1＝p^k+δP^k (7)

式中符号定义：δ为算子，表示第k次迭代到k+1次迭代后的变化；n为上一时步；n+1为下一时步。

基质单元矩阵为：

裂缝单元矩阵为：

式中符号定义：P、T、W和M分别为压力矩阵、传导率矩阵、窜流项矩阵和时间项矩阵；δP、δT、δW和δM则为相应的算子矩阵；△t为时间步长。将压力和饱和度的试探解带入渗流方程中，可推导出P、T、W和M的具体计算格式为：

S3.2：对离散区域内的每一个四面体网格，建立类似的单元矩阵。假设存在N个网格节点，通过整体叠加，形成关于基质和裂缝***的4N×4N大矩阵：

[K]_4N*4N[δX]_4N*1＝[R]_4N*1 (15)

式中符号定义：K为整体系数矩阵；δX为未知变量变化量矩阵，包括δP_m，δP_f，δS_wm和δS_wf；R为余量项。

5.根据权利要求1所述的一种基于非结构网格的水平井气水两相数值模拟方法，其特征在于：所述步骤S4还包括：

由Peaceman公式，水平井总气水产量可表示为：

式中符号定义：N_w为水平井穿过的离散网格节点；PI为生产指数；p_well为井点节点处水平井筒压力；p_ave为井点节点所在网格的平均压力；

考虑水平井筒内流动特征，采用水平井筒多相流计算公式获得考虑井筒摩阻和流体动量损失的井筒压力沿程梯度：

式中符号含义：ρ_l为液体密度，kg/m³；ρ_g为气体密度，kg/m³；G为气液混合物质量流量，kg/s；v_m为混合物流动速度，m/s；v_sg为气体表观流速，m/s；A为井筒油管截面积，m²；D为油管内径，m。

持液率H_L和摩擦阻力系数f采用Beggs-Brill方法计算。水平井跟端网格点井筒压力p_well1等于井底压力值p_bh，则结合式(17)可以得到各井点网格处井筒压力值：

将得到的井筒压力值p_welli带入到式(16)中，结合井点网格压力p_avei，得到每时步下的产量值。

6.根据权利要求1所述的一种基于非结构网格的水平井气水两相数值模拟方法，其特征在于：所述步骤S5还包括以下步骤:

S5.1：由有限元基本原理，对于定压外边界条件，通过置大数法，将定值(p_e)向量添加到式(15)载荷矩阵K中；对于封闭外边界条件，添加0向量到式(15)载荷矩阵K中；

S5.2：将式(16)中下一时步总产量转化为并带入到式(15)中，得到考虑内外边界条件的数值计算稀疏矩阵；

S5.3：采用Orthomin共轭梯度法对稀疏矩阵进行求解，通过牛顿-拉夫逊迭代方法，获得基质和裂缝***在一个时步下的压力和饱和度变化值：δp_gf1,…δp_gfN；δp_gm1,…δp_gmN；

δS_wf1,…δS_wfN；δS_wm1,…δS_wmN。进而由式(7)和式(8)获得下一时步的值；

S5.4：输出各时步下的离散网格节点压力、饱和度值，以及水平井井底流压和产量值。