CN111339671B - 页岩储层双向流-固耦合数值计算方法 - Google Patents

Abstract

本发明公开页岩储层双向流‑固耦合数值计算方法，该方法是针对应力敏感性页岩储层的流‑固耦合模拟仿真数值计算方法进行改进，提出一种新的迭代形式，在原有的弱耦合数值计算方法基础上，附加一个循环，通过该附加的循环将室内岩芯力学实验数据及每个时间步计算出的体应变进行相互校正，直至满足收敛条件后跳出附加循环，进入下一个时间步的计算，基于本发明的数值算法，可以达到岩芯室内力学实验数据与理论计算进行相互校正的目的，更准确地表征地质力学响应；本发明的方法在保持原有数值计算方法优点的同时，强化了相关物理过程，在有效捕抓物理现象的同时，保证计算效率及稳定性，为流‑固耦合数值计算提供一些新方法和新思路。

Description

页岩储层双向流-固耦合数值计算方法

技术领域

本发明所属技术领域为非常规油气开发中页岩储层模拟仿真数值计算领域，具体涉及页岩储层双向流-固耦合数值计算方法。

背景技术

作为一种新型的非常规天然气资源，近些年页岩气受到高度的关注。页岩储层中有机孔、无机孔及应力敏感的狭缝形裂隙共同发育，页岩储层复杂的孔隙结构及流体赋存状态给页岩储层模拟仿真带来了极大的挑战，针对应力敏感性储层的模拟，目前主要的研究方法主要分为如下两类：1、基于相关的解析或经验统计模型结合实时的储层压力响应数据校正孔隙介质物性参数(孔隙度、渗透率)从而刻化应力敏感性储层；2、基于经典pore-thermal-elastic理论及应力平衡方程，构建地质力学模型，并耦合流体方程求解。方法1只需求解流体方程，计算效率高，但过于简化了相关物理现象，该方法无法反映应力等地质力学物理量的动态演化及相关影响。方法2基于pore-thermal-elastic理论，需要求解流体方程及应力平衡方程，严格耦合了应力场，真正意义上实现了多物理耦合，但方法2较方法1计算成本会相对较高。双向分布迭代耦合方法被广泛应用于严格耦合应力场，流-固信息交互的媒介为孔隙度和渗透率，仅通过渗透率来耦合流体方程和应力平衡方程的方式被称为弱耦合模式，由于弱耦合模式只通过渗透率耦合，该方式具备较好的数值稳定性及计算效率，但弱化了相关物理过程。针对应力敏感性储层的模拟，现有技术存在的问题主要集中在如何有效平衡物理现象客观描述及数值计算效率、稳定性及准确性之间的关系。现有技术或弱化了流-固耦合过程中涉及到的物理现象或需要付出高昂的计算成本以真实还原相关过程中涉及到的物理现象。

发明内容

为解决现有技术中存在的问题，本发明的目的在于提供页岩储层双向流-固耦合数值计算方法，本发明能够有效抓捕、还原双向流-固耦合过程中涉及到的物理现象，同时保证计算效率及稳定性。本发明旨在改进弱耦合模式在页岩储层数值仿真中的应用，在保持该方式数值计算稳定性和效率优势的同时进一步强化相关物理过程。

本发明采用的技术方案如下：

页岩储层双向流-固耦合数值计算方法，包括如下过程：

在弱耦合数值计算方法基础上，附加一个循环，通过附加的循环将室内岩芯力学实验数据及每个时间步计算出的体应变进行相互校正，直至满足收敛条件后跳出附加循环，进入下一个时间步的计算。

优选的，附加循环的过程包括如下步骤：

S1，通过油藏模型(即流体物质平衡方程)计算出孔隙压力；

S2，通过应力平衡方程和孔隙压力计算平均应力，利用平均应力计算有效应力；

S3，基于有效应力与孔隙半径的实验数据得到当前孔隙半径r_p1；

S4，通过有效应力与体形变的关系及体变形与孔隙半径关系得到当前孔隙半径r_p2；

S5，计算r_p1与r_p2之差的绝对值，如果绝对值满足收敛条件，则以r_p1或r_p2更新表观渗透率，然后将更新后的表观渗透率传递给油藏模型，进入下一个时间步的计算；

如果如果r_p1与r_p2之差的绝对值不满足收敛条件，则通过附加循环迭代求解，直至r_p1与r_p2之差的绝对值满足收敛条件为止，然后将r_p1与r_p2的均值传递给油藏模型，基于表观渗透率从而更新渗透率，进入下一个时间步的计算。

优选的，S1中，流体物质平衡方程为：

其中：ρ_g为自由气体密度；ρ_a为吸附气体密度；u_g为达西流速；φ为油藏孔隙度；t为时间；

达西流速u_g为：

其中：φ为油藏孔隙度；μ为气体粘度；

为压力梯度；K_a为页岩储层的表观渗透率；

页岩储层的表观渗透率K_a为：

其中：K_a：页岩储层的表观渗透率；μ为气体粘度；ρ_avg为气体平均密度；P_avg为毛管中的平均压力；R为气体常数；α为切向动量调节系数；T为温度；M为气体分子量；r_p为实时孔隙等效半径；

以毛管中的平均压力P_avg作为孔隙压力。

优选的，S2中，应力平衡方程为：

其中：σ′为有效应力；σ_mean为平均应力；α为Biot常数；P为孔隙压力；v为泊松比；F_b为体力；

有效应力σ′为：

σ′＝σ_mean-αP

优选的，S3中，当前孔隙半径r_p1为：

r_p1＝r₀*γ

其中：r₀为初始孔隙半径；γ为孔径折减系数。

优选的，S4中，体应变与有效应力关系为：

σ′＝Kε_v

其中：K为模量；σ′为有效应力；ε_v为体应变；

体应变与孔隙半径具有以下关系：

上式中：V_p为孔隙体积；V_B ⁰为初始视体积；n_p为毛管数量；L₀为毛管长度；φ为真实孔隙度；φ₀为初始真实孔隙度；C_p为孔隙压缩系数；r_p2为前孔隙半径。

本发明具有如下有益效果：

本发明的页岩储层双向流-固耦合数值计算方法是针对应力敏感性页岩储层的流-固耦合模拟仿真数值计算方法进行改进，提出一种新的迭代形式，在原有的弱耦合数值计算方法基础上，附加一个循环，通过该附加的循环将室内岩芯力学实验数据及每个时间步计算出的体应变进行相互校正，直至满足收敛条件后跳出附加循环，进入下一个时间步的计算，基于本发明的数值算法，可以达到岩芯室内力学实验数据与理论计算进行相互校正的目的，更准确地表征地质力学响应；本发明的方法在保持原有数值计算方法优点的同时，强化了相关物理过程，在有效捕抓物理现象的同时，保证计算效率及稳定性，为流-固耦合数值计算提供一些新方法和新思路。

附图说明

图1为本发明页岩储层双向流-固耦合数值计算方法的流程图。

图2为有效应力与孔径折减系数的关系图。

具体实施方式

下面结合实施例来对本发明做进一步的说明。

本发明主要针对基于经典pore-thermal-elastic理论及应力平衡方程，构建地质力学模型，并耦合流体方程求解的弱耦合模式进行改进，提出一种新的迭代形式，在原有的数值计算方法基础上，附加一个循环，通过该附加的循环将室内岩芯力学实验数据及每个时间步计算出的体应变进行相互校正，直至满足收敛条件后跳出附加循环，进入下一个时间步的计算，改进的新方法在保持原有数值计算方法优点的同时，强化了相关物理过程，进一步完善了地质力学在双向耦合过程中的响应

本发明页岩储层双向流-固耦合数值计算方法的过程如下：

首先本发明通过流体物质平衡方程和孔隙压力计算出孔隙压力，然后通过应力平衡方程计算平均应力，从而得到有效应力，基于有效应力与孔隙半径的实验数据得到当前孔隙半径r_p1，同时通过有效应力与体形变关系及有效应力与孔隙半径关系得到当前孔隙半径r_p2，如果r_p1与r_p2之差的绝对值足够小(即满足收敛条件)，则将新的孔隙半径(取r_p1或r_p2)更新表观渗透率，然后将更新后的表观渗透率传递给油藏模型，进入下一个时间步的计算；如果如果r_p1与r_p2之差的绝对值不满足收敛条件，则通过附加循环迭代求解直至满足收敛条件为止，基于该数值算法，可以达到岩芯室内力学实验数据与理论计算进行相互校正的目的，更准确地表征地质力学响应。本实施例页岩储层双向流-固耦合数值计算方法具体计算流程图如图1所示，具体包括如下过程：(1)流体物质平衡：

页岩储层的表观渗透率：

上式中：K_a：页岩储层的表观渗透率，单位为m²；

μ：气体粘度，单位为pa.s；

ρ_avg：气体平均密度，单位为kg/m；

R：气体常数；

P_avg为毛管中的平均压力，单位为Pa；

α：切向动量调节系数；

T：温度，单位为K；

M：气体分子量，单位为kg/mol；

r_p：实时孔隙等效半径，m；

达西定律：

上式中：u_g：达西流速，单位为m/s；

φ：油藏孔隙度；

μ：气体粘度，单位为Pa.s；

压力梯度，单位为Pa/m；

物质平衡方程：

上式中：ρ_g：自由气体密度，单位为kg/m³；

ρ_a：吸附气体密度，单位为kg/m³；

压力梯度，单位为Pa/m；

φ：油藏孔隙度；

(2)地质力学模型：

应力平衡方程：

上式中：σ′：有效应力，单位为Pa；

σ_mean：平均应力，单位为Pa；

α：Biot常数；

P：孔隙压力，单位为Pa；

v：泊松比；

F_b：体力(重力)，单位为Pa；

有效应力：

σ′＝σ_mean-αP

上式中：σ′：有效应力，单位为Pa；

体应变与有效应力关系：

σ′＝Kε_v

上式中：K：体模量，单位为Pa；

ε_v：体应变；

体应变与孔隙半径关系：

上式中：V_p：孔隙体积，单位为m³；

V_B ⁰：初始视体积，单位为m³；

n_p：毛管数量；

L₀：毛管长度，单位为m；

φ：真实孔隙度；

φ₀：初始真实孔隙度；

C_p：孔隙压缩系数，单位为Pa^-1；

r_p2：孔隙半径，单位为m；

有效应力与孔隙半径关系-室内岩芯力学实验数据：

r_p1＝r₀*pore radius reduction coefficient＝r₀*γ

上式中：r₀：初始孔隙半径，单位为m；

r_p1：孔隙半径，单位为m；

γ为孔径折减系数，孔径折减系数可根据图2所示的关系进行选取，有效应力越大，孔径折减系数越小，表明实际的孔隙半径也越小。

Claims (1)

1.页岩储层双向流-固耦合数值计算方法，其特征在于，包括如下过程：

在弱耦合数值计算方法基础上，附加一个循环，通过附加的循环将室内岩芯力学实验数据及每个时间步计算出的体应变进行相互校正，直至满足收敛条件后跳出附加循环，进入下一个时间步的计算；

附加循环的过程包括如下步骤：

S1，通过油藏模型计算出孔隙压力；

S2，通过应力平衡方程和孔隙压力计算平均应力，利用平均应力计算有效应力；

S3，基于有效应力与孔隙半径的实验数据得到当前孔隙半径r_p1；

S4，通过有效应力与体形变的关系及体变形与孔隙半径关系得到当前孔隙半径r_p2；

S5，计算r_p1与r_p2之差的绝对值，如果绝对值满足收敛条件，则以r_p1或r_p2更新表观渗透率，然后将更新后的表观渗透率传递给油藏模型，进入下一个时间步的计算；

如果r_p1与r_p2之差的绝对值不满足收敛条件，则通过附加循环迭代求解，直至r_p1与r_p2之差的绝对值满足收敛条件为止，然后将r_p1与r_p2的均值传递给油藏模型，基于表观渗透率从而更新渗透率，进入下一个时间步的计算；

S1中，油藏模型为：

其中：ρ_g为自由气体密度；ρ_a为吸附气体密度；u_g为达西流速；φ为油藏孔隙度；t为时间；

达西流速u_g为：

其中：φ为油藏孔隙度；μ为气体粘度；▽P为压力梯度；K_a为页岩储层的表观渗透率；

页岩储层的表观渗透率K_a为：

其中：K_a：页岩储层的表观渗透率；μ为气体粘度；ρ_avg为气体平均密度；P_avg为毛管中的平均压力；R为气体常数；α为切向动量调节系数；T为温度；M为气体分子量；r_p为实时孔隙等效半径；

以毛管中的平均压力P_avg作为孔隙压力；

S2中，应力平衡方程为：

其中：σ′为有效应力；σ_mean为平均应力；α为Biot常数；P为孔隙压力；v为泊松比；F_b为体力；

有效应力σ′为：

σ′＝σ_mean-αP；

S3中，当前孔隙半径r_p1为：

r_p1＝r₀*γ

其中：r₀为初始孔隙半径；γ为孔径折减系数；

S4中，体应变与有效应力关系为：

σ′＝Kε_v

其中：K为模量；σ′为有效应力；ε_v为体应变；

体应变与孔隙半径具有以下关系：

上式中：V_p为孔隙体积；V_B ⁰为初始视体积；n_p为毛管数量；L₀为毛管长度；φ为真实孔隙度；φ₀为初始真实孔隙度；C_p为孔隙压缩系数；r_p2为前孔隙半径。

Images