Solution of Nonlinear Seepage Model for Fracture Well Groupin Low Permeability Reservoirs
-
摘要:
精细描述低渗透油藏中流体流速与与压力梯度的非线性关系，是准确计算低渗透油藏压裂井组产量的基础。为此，在描述低渗透油藏非线性渗流特征的基础上，建立了低渗透油藏和压裂裂缝耦合的非线性数学模型，该模型根据渗流特征将渗流过程分为非线性渗流阶段和拟线性渗流阶段进行计算。利用Taylor展开对非线性数学模型进行线性化处理，建立了有限差分方程组，并编制了计算机求解程序。算例分析表明：采用非线性数学模型计算出的地层中压力和饱和度的分布符合地层实际情况；五点法井网压裂井组注水井的裂缝导流能力会随着裂缝闭合而降低，注水效果变差，导致油井产量降低。研究结果表明，低渗透油藏和压裂裂缝耦合的非线性数学模型可以较准确地描述低渗透油藏中流体流速与压力梯度的非线性关系，为准确计算低渗透油藏压裂井组产量奠定基础，为低渗透油藏注水开发提供指导。
Abstract:Having a closely detailed description of the nonlinear relationship between flow velocity and pressure gradient in low permeability reservoir is necessary for accurately developing the frac design, and calculating the production of a group (or unit) of wells that have been hydraulically fractured. Therefore, based on the description of the nonlinear seepage characteristics of low permeability reservoir, a nonlinear mathematical model of coupling low permeability reservoir and hydraulic fractures was established, which divided the seepage process into the nonlinear seepage stage and quasi-linearity stage according to the seepage characteristics. The Taylor expansion was used to linearize the nonlinear mathematical model, and established the finite difference equations, and then formed the computer solving model. The results of example analysis showed that the distributions of formation pressure and saturation calculated by the nonlinear mathematical model were in line with the actual situations of the stratum; the fracture flow conductivity of injection well in the fractured five-spot well pattern decreased with the formation closure, which led to poor water injection effect and low oil well production. Thus, the fracture design should be modified in accordance with the study’s results. The study results indicated that the nonlinear mathematical model and hydraulic fracture coupling could accurately describe the nonlinear relationship between flow velocity and pressure gradient in low-permeability reservoir. This breakthrough establishes a foundation to calculate the production of fractured well group in low-permeability reservoir accurately, and provides a guidance for water flooding development of low permeability reservoir.
-
随着我国中浅层油气藏开发进入中后期，油气勘探目的层逐渐转向深层（大于 4500 m）和超深层（大于6000 m）[1–4]，深层超深层油气藏具有高温、高压和高应力的特点[5–6]，显著影响储层岩石孔隙结构[7]，从而影响油气渗流[8–9]。数字岩心是数值模拟研究油气微观渗流规律的载体，重构数字岩心是开展油气微观渗流数值模拟的重要步骤。现有重构方法分为物理实验法和数值模拟法，其中物理实验法是指用CT等成像技术结合图像处理算法构建三维数字岩心[10–11]；数值模拟法包括随机模拟算法[12]、过程模拟法[13–14]和深度学习算法[15–17]，其中随机模拟算法、深度学习算法均依赖CT扫描图像，但是现有CT扫描时承压能力仅为30～40 MPa，不能满足深层高应力的需求[18]。过程模拟法则通过模拟岩石颗粒的沉积、压实和成岩过程，构建能反映应力影响的数字岩心。离散元法（discrete element method, DEM）是实现过程模拟法的手段[19–21]，DEM建模需要解决以下问题：一是分割CT图像中相互接触的颗粒，通常采用分水岭算法[22]；二是颗粒轮廓的提取与分析，通常用一个连续函数表示轮廓曲面，包括超二次椭球函数[23]、非均匀有理B样条函数[24]、傅里叶级数[25]和球面谐波函数[26]，其中在三维情况下应用最为成熟的是球面谐波函数分析[27]；三是将颗粒轮廓转化为用于离散元法计算的颗粒，一种是重叠离散元颗粒簇（overlapping discrete element cluster, ODEC）方法填充轮廓形成的clump（团簇）作为颗粒计算[28-30]，另一种是直接使用提取的轮廓作为颗粒计算，其中前者使用最为广泛[31]。
根据上述分析，笔者提出了基于离散元法考虑深层超深层油气藏高应力影响的数字岩心构建方法。首先，基于常温常压下的扫描结果，采用球面谐波函数分析得到形状和粒径分布，并转化为Clump构建模板库，根据粒径分布从模板库中选取Clump进行建模，实现常温常压数字岩心从图像到离散元模型的转化；然后，施加应力边界模拟深层超深层的高地应力，并使用PFC3D进行模拟计算；最后，根据颗粒的位置和半径利用体素化算法转化为数字岩心，分析了数字岩心孔隙度、孔隙结构和渗透率在不同应力下的变化规律，并利用Bentheim算例验证了考虑高应力影响的数字岩心重构方法的可行性。研究结果为深层超深层油气藏孔隙尺度模拟提供了技术途径。
1. 数字岩心构建方法及原理
1.1 离散元法理论基础
离散元法的基本原理是使用相互接触的颗粒建模，每个颗粒的运动均符合牛顿第二定律，采用显式中心差分格式求解每个颗粒的运动，由各颗粒的运动得到整体力学性质。颗粒i的平动和转动方程分别为：
{m}_{i}\frac{{\rm{d}}{\boldsymbol{v}}_{i}}{{\rm{d}}t}+{\beta }_{\mathrm{g}}{\boldsymbol{v}}_{i}\left(t\right)={\boldsymbol{F}}_{\mathrm{t}i} (1) {I}_{i}\frac{{\rm{d}}{\boldsymbol{\omega }}_{i}}{{\rm{d}}t}+{\beta }_{\mathrm{g}}{{\boldsymbol{\omega}} }_{i}\left(t\right)={\boldsymbol{M}}_{i} (2) 式中：mi为第i个颗粒的质量，kg；Ii为第i个颗粒分转动惯量，kg·m²；vi为第i个颗粒速度，m/s； ωi为第i个颗粒角速度，rad/s；βg为整体阻尼，kg/s；Fti为第i个颗粒合外力，N；Mi为第i个颗粒合力矩，N·m。
使用显式有限差分算法更新每个颗粒的速度和转动角度，并更新颗粒的位置。显式方法中，时间步的选择是保证数值稳定性的关键，PFC3D中定义了临界时间步长：
\Delta {t}_{\text{crit}}=\sqrt{\frac{m}{k}} (3) 式中：Δtcrit为临界时间步长，s；m为颗粒的质量，kg；k为接触弹簧的刚度，N/m。
颗粒位置更新后，需要重新确定颗粒的接触关系，常用的接触搜索算法为基于网格的接触搜索算法。确定新的接触关系后，根据颗粒之间的本构关系计算颗粒间的接触力。模拟岩石材料的力学行为常用的接触模型为线性平行黏结模型（见图1）。
线性平行黏结模型中，颗粒间接触力和力矩的计算公式为：
{\boldsymbol{F}}_{\mathrm{c}}={\boldsymbol{F}}^{1}+{\boldsymbol{F}}^{\mathrm{d}}+\bar{\boldsymbol{F}} (4) {\boldsymbol{M}}_{\mathrm{c}}=\bar{\boldsymbol{M}} (5) \,其中\qquad\qquad\qquad\;\; \;\;{\boldsymbol{F}}^{1}={\boldsymbol{F}}_{\mathrm{n}}^{1}+{\boldsymbol{F}}_{\mathrm{s}}^{1}\qquad\qquad (6) {\boldsymbol{F}}^{\mathrm{d}}={\boldsymbol{F}}_{\mathrm{n}}^{\mathrm{d}}+{\boldsymbol{F}}_{\mathrm{s}}^{\mathrm{d}} (7) 式中：Fc为接触力，N；Mc为接触力矩，N·m；Fl为线性力，N；Fd为阻尼力，N；
{\bar {\boldsymbol{F}}} 为平行黏结力，N；{\bar {\boldsymbol{M}} } 为平行黏结力矩，N·m；{\boldsymbol{F}}_{\mathrm{n}}^{\mathrm{l}} 为法向线性力，N；{\boldsymbol{F}}_{\mathrm{s}}^{\mathrm{l }} 为切向线性力，N；{\boldsymbol{F}}_{\mathrm{n}}^{\mathrm{d}} 为法向阻尼力，N；{\boldsymbol{F}}_{\mathrm{s}}^{\mathrm{d}} 为切向阻尼力，N。1.2 球面谐波函数分析及颗粒轮廓重构
球面谐波函数（spherical harmonic function）分析是傅里叶级数展开分析的一种，是用球面谐波函数线性组合表达式确定的闭合曲面来表征颗粒轮廓形状：
\mathit{r}(\theta,\varphi)=(x\left(\theta,\varphi\right),y\left(\theta,\varphi\right),z\left(\theta,\varphi\right))^T=\sum_{n=0}^{n_{\mathrm{m}\mathrm{a}\mathrm{x}}}\sum_{m=-n}^n\mathit{c}_n^m\mathrm{Y}_n^m\left(\theta,\varphi\right) (8) \,其中 \qquad \qquad\;\;\;\;{\boldsymbol{c}}_{n}^{m}={\left({c}_{x,n}^{m},{c}_{y,n}^{m},{c}_{z,n}^{m}\right)}^{{\rm{T}}} \qquad\qquad (9) \mathrm{Y}_n^m\left(\theta,\varphi\right)=\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{4\text{π}\left(n+m\right)!}}\mathrm{P}_n^m\left(\cos\theta\right)\rm{e}^{im\phi} (10) 式中： nmax表示最高阶数，一般取为15；n为阶数；m为次数；
{\boldsymbol{c}}^m_n 为系数矩阵；\mathrm{Y}_n^m 为球面谐波函数；\mathrm{P}_n^m 表示关联勒让德函数。Wei Deheng等人[24]提出了一种系数矩阵的重构算法。首先，定义球面描述符：
\begin{aligned}{d}_{n}=&{\|{\boldsymbol{C}}_{n}\|}_{{\rm{F}}}=\sqrt{\sum _{i\in \left(x,y,z\right)}\sum _{m=-n}^{n}{\|{c}_{i,n}^{m}\|}^{2}}=\\ &\quad\;\;\sqrt{\sum _{i\in \left(x,y,z\right)}\sum _{m=-n}^{n}{c}_{i,n}^{m}{c}_{i,n}^{{m}^{\mathrm{*}}}} \end{aligned} (11) 式中：Cn表示阶数为n的所有系数组成的矩阵；
{\|\cdot \|}_{{\rm{F}}} 表示Frobenius范数；*表示共轭转置。用幂函数拟合颗粒的球面描述符与阶数n的关系，得到参数α和β：
d_{\rm{\mathit{n}}}=\left\{\begin{array}{l}d_2\left(\dfrac{2}{n}\right)^{\alpha}\qquad n\in\left[\mathrm{2,8}\right] \\ d_9\left(\dfrac{9}{n}\right)^{\beta}\qquad n\in\left[\mathrm{9,15}\right]\end{array}\right. (12) 随后，分析球面描述符主成分，定义4个旋转不变球面谐波因子：
\left\{\begin{array}{l}I_\text{E}=\dfrac{p_2}{p_1} \\ I_\text{F}=\dfrac{p_3}{p_2}\end{array}\ \ \ \ \ \quad\left(p_1 > p_2 > p_3\right)\right. (13) \left\{\begin{array}{l}\dfrac{{d}_{2-8}}{{d}_{1}}=\displaystyle\sum _{n=2}^{8}\dfrac{{d}_{n}}{{d}_{1}}\\ \dfrac{{d}_{9-15}}{{d}_{1}}=\displaystyle\sum _{n=9}^{15}\dfrac{{d}_{n}}{{d}_{1}}\end{array}\right. (14) 式中：IE为延伸指数；IF为平面指数；d2-8/d1表征颗粒圆度；d9-15/d1表征颗粒的细观形状特征。
最后，根据旋转不变球面谐波因子分布规律重构系数矩阵，使用IE和IF计算C1：
{\mathit{C}}_{1}=-\sqrt{\frac{{\text{π}} }{6}}\left(\begin{array}{ccc}-1& 0& 1\\ I_\text{E} \text{i}& 0& I_\text{E} \text{i}\\ 0& \sqrt{2}I_\text{E} I_\text{F}& 0\end{array}\right) (15) 根据d2-8/d1和d9-15/d1的值结合式（12）计算d2～d15，并假设3个维度上的球面描述符相等，构造相应的系数矩阵C2～C15：
{\boldsymbol{C}}_{n}=\left(\begin{array}{ccccc}{c}_{x,n}^{-n}& \cdots & {c}_{x,n}^{0}& \cdots & {c}_{x,n}^{n}\\ {c}_{y,n}^{-n}& \cdots & {c}_{y,n}^{0}& \cdots & {c}_{y,n}^{n}\\ {c}_{z,n}^{-n}& \cdots & {c}_{z,n}^{0}& \cdots & {c}_{z,n}^{n}\end{array}\right)=\left(\begin{array}{ccccc}{k}_{x}\left({\alpha }_{2n}-{\alpha }_{2n+1}\text{i}\right)& \cdots & {k}_{x}{\alpha }_{1}& \cdots & {k}_{x}\left({\alpha }_{2n}+{\alpha }_{2n+1}\text{i}\right)\\ {k}_{y}\left({\beta }_{2n}-{\beta }_{2n+1}\text{i}\right)& \cdots & {k}_{y}{\beta }_{1}& \cdots & {k}_{y}\left({\beta }_{2n}+{\beta }_{2n+1}\text{i}\right)\\ {k}_{z}\left({\varepsilon }_{2n}-{\varepsilon }_{2n+1}\text{i}\right)& \cdots & {k}_{z}{\varepsilon }_{1}& \cdots & {k}_{z}\left({\varepsilon }_{2n}+{\varepsilon }_{2n+1}\text{i}\right)\end{array}\right) (16) \,其中\qquad \left\{\begin{array}{l}{k}_{x}={d}_{x,n}^{2}/\left({\alpha }_{1}^{2}+2{\displaystyle\sum }_{i=2}^{2n+1}{\alpha }_{i}^{2}\right)\\ {k}_{y}={d}_{y,n}^{2}/\left({\beta }_{1}^{2}+2{\displaystyle\sum }_{i=2}^{2n+1}{\beta }_{i}^{2}\right)\\ {k}_{z}={d}_{z,n}^{2}/\left({\varepsilon }_{1}^{2}+2{\displaystyle\sum }_{i=2}^{2n+1}{\varepsilon }_{i}^{2}\right)\end{array}\right. (17) 式中：αi, βi, εi为0～1之间符合均匀分布的随机数，kx, ky和kz为归一化因子。
1.3 离散元法数字岩心构建流程及评价方法
基于离散元法重构高应力影响的数字岩心构建方法为：1）建立常温常压数字岩心离散元模型，先使用1.2节方法构造符合实际形状特征的数据库，并转化为Clump模板库，再根据粒径分布从模板库中随机选取颗粒建模，并判断孔隙度是否达到实际岩心的孔隙度，不断调整每种Clump的比例，直至达到真实孔隙度（见图2），最后对比重构结果与实际岩心的两点相关函数曲线和线性路径相关函数曲线；2）在模型的边界施加应力，并在PFC3D中模拟伺服加载至设定应力，模拟深层超深层储层岩石不同方向受到的应力；3）将离散元模型体素化为数字岩心，提取孔隙网络模型分析孔隙几何拓扑结构，计算孔隙度和渗透率。
![]()
图 2 常温常压数字岩心离散元建模流程Figure 2. Flowchart of digital core modeling at normal temperature and pressure by DEM2. 常温常压数字岩心离散元建模
Bentheim砂岩具有完善的CT扫描图像和力学实验数据，可以用来验证考虑高应力数字岩心构建方法的可行性。通过颗粒分割与提取、球面谐波分析重构颗粒形状，将岩心图像转化为颗粒，用于离散元建模，得到常温常压下数字岩心的离散元模型。
2.1 颗粒分割与提取
选用的Bentheim砂岩图像大小为891×891×1400，分辨率为4.901 46 μm。为消除边界效应，选取300×300×600体素的区域作为分析对象，孔隙度为22.95%，使用分水岭算法处理图像，分割相互接触的颗粒并保存每个颗粒的所有体素，共提取500个完整颗粒。
2.2 球面谐波分析重构颗粒形状
提取图像中每个颗粒所有体素后，首先对每个颗粒进行球面谐波展开，阶数设置为15，计算每个颗粒的球面谐波函数系数矩阵，并保存相应的图形文件和系数矩阵；接着，分析了颗粒的等效粒径、纵横比（IE和IF的平均值）、球度、凸度和圆度的分布曲线（见图3），为Clump模板库及常温常压数字岩心离散元建模提供依据。
![]()
图 3 样品粒径及颗粒形状分析结果Figure 3. Analysis results of particle size and particle shapes of the sample接着分析球面描述符与阶数之间的关系，拟合得到式（12）中的参数α和β。在双对数坐标系中绘制了颗粒1～500的球面描述符与球面谐波函数阶数（SH阶数）之间的关系曲线（见图4），用幂函数进行拟合，拟合结果为：α = 0.743 7，β = 1.006。
![]()
图 4 球面描述符与SH阶数的关系曲线Figure 4. Relationship between the spherical descriptor and SH degree根据形状分析结果，设置不同比例的IE、IF、d2-8/d1和d9-15/d1，结合参数α、β的拟合结果重构球面谐波函数系数矩阵C，得到符合实际形状分布特征的颗粒轮廓数据库，数据库大小为50。使用ODEC方法转化为PFC3D中的Clump颗粒，形成模板库用于后续建模。
2.3 常温常压数字岩心的离散元建模
通过以上步骤，将常温常压下的数字岩心转化为离散元颗粒，随后根据离散元建模方法，实现常温常压数字岩心从图像到离散元模型的转化。常温常压数字岩心的离散元模型如图5（a）所示，模型尺寸为1.5 mm×1.5 mm×1.5 mm，使用边长为3 μm的体素对模型进行体素化，转化为数字岩心，体素尺寸为500×500×500（见图5（b））。从图5（b）可以看出，孔隙度受到边界效应的影响，边界处的孔隙度大，内部孔隙度小。
![]()
图 5 常温常压数字岩心离散元建模结果Figure 5. Results of digital core modeling at normal temperature and pressure by DEM不同表征单元体（representative elementary volume，REV）边长与孔隙度间的关系如图6所示，可以发现REV边长为200～400时，孔隙度较分布平缓。综合以上2个方面，选择REV为300，体素范围100～400，此时模型孔隙度为22.88%，与真实岩心的相对误差为0.30%。
重构数字岩心REV内的两点相关函数（S2）和线性路径相关函数（L）曲线与真实岩心的对比结果如图7所示。
![]()
图 7 重构岩心与真实岩心的两点相关和线性路径相关函数曲线对比结果Figure 7. Comparison of two-point correlation and linear path correlation function curves between reconstructed and real cores从图7可以看出，距离r为0时，S2和L的数值为模型的孔隙度（重构数字岩心和真实岩心分别为22.88%和22.95%）；r趋近于无穷大时，S2接近于孔隙度的平方，L趋近于0。此外，重构数字岩心和真实岩心的2个相关函数曲线基本重合，表明重构的常温常压数字岩心离散元模型能够反映真实岩心的孔隙结构。
根据图3（a）的粒径分布曲线，将粒径等比例放大15倍，建立与力学试验相同尺寸（直径2 cm，高度4 cm）的离散元模型[32]，如图8（a）所示；通过调整平行黏结接触模型中的参数，使模拟得到的应力应变曲线与实际应力应变曲线相吻合，结果如图8（b）所示，此时模型的微观力学参数：有效模量和黏结模量为1.30 GPa，法向-切向刚度比和黏结法向-切向刚度比均为1.50，抗拉强度0.70 MPa，内聚力0.50 MPa，摩擦角55.0°。
![]()
图 8 颗粒间微观力学参数标定模型及偏应力-应变曲线Figure 8. Micromechanical parameter calibration model among particles and the deviational stress-strain curve3. 考虑应力作用的数字岩心构建与分析
首先，分别设置3种围压（50/70, 70/90和90/110 MPa）和轴向应力（120, 150和180 MPa），在PFC3D中设置不同应力组合模拟应力加载，得到受应力压实后的离散元模型；然后，使用体素化算法将离散元模型转化为数字岩心，采用最大球算法提取不同应力作用下数字岩心的孔隙网络，分析不同应力下数字岩心的孔隙度、孔隙半径、喉道半径、喉道长度和比欧拉示性数的变化，并计算了渗透率。
3.1 数字岩心构建结果
不同应力组合条件下的数字岩心构建结果如表1所示。从表1可以看出，围压相同的情况下，轴向应力越大，颗粒的位移量越大；轴向应力相同的情况下，颗粒的位移量会随围压增大而减小。颗粒的位移导致岩心孔隙几何拓扑结构的改变，最终造成不同应力下数字岩心孔隙度和渗透率的差异。
表 1 不同应力组合下的数字岩心构建结果Table 1. Digital core construction results under different stress combinationss水平应力 σz = 120 MPa σz = 150 MPa σz = 180 MPa σx = 50 MPa
σy = 70 MPa
σx = 70 MPa
σy = 90 MPa
σx = 90 MPa
σy = 110 MPa
3.2 数字岩心孔隙几何拓扑结构分析
相同围压、不同轴向应力下的数字岩心孔隙几何拓扑结构如图9所示。图9（a）是围压σx = 70 MPa、σy = 90 MPa，轴向应力不同时的孔隙半径分布曲线，可以看出在应力作用下，孔隙半径的最大值由4.27 μm减小至3.70 μm；由于孔隙闭合，孔隙半径小于0.40 μm孔隙的比例略有减小；孔隙半径0.40～2.00 μm范围内分布情况剧烈波动。图9（b）为喉道半径分布曲线，随着轴向应力增大，喉道半径的最大值由3.41 μm减小至3.06 μm，小于0.50 μm和大于1.50 μm喉道的比例减小；0.50～1.50 μm范围内喉道的比例波动剧烈，但整体比例增大。图9（c）为喉道长度分布曲线，随着轴向应力增大，分布曲线整体右移，最大喉道长度由17.56 μm增大至21.88 μm。图9（d）为比欧拉示性数曲线，其与直线y = 0的交点值表征数字岩心连通性，其值越大，表示连通性越好；应力增大，曲线与直线y = 0交点值减小，表明岩心的连通性变差。
![]()
图 9 相同围压、不同轴压下的数字岩心孔隙几何拓扑结构Figure 9. Pore geometry topology of digital core under same confining pressure and different axial pressures相同轴向应力、不同围压下的数字岩心孔隙几何拓扑结构如图10所示。图10（a）和（b）分别为轴向应力为150 MPa时不同围压下数字岩心的孔隙半径分布曲线和喉道半径分布曲线，与常温常压数字岩心相比，应力的作用会造成孔隙半径和喉道半径的最大值减小，但是围压大的模型，减小幅度越小，能够保留更多的小半径孔隙和喉道。图10（c）为不同围压下数字岩心的喉道长度分布曲线，可以发现喉道长度分布曲线随围压增加向右移动，分布范围变大。图10（d）为不同围压下数字岩心比欧拉示性数曲线，其与直线y = 0交点值随围压增大而减小，表明数字岩心的连通性变差。
![]()
图 10 相同轴向应力、不同围压下的数字岩心孔隙几何拓扑结构Figure 10. Pore geometry topology of digital core under same axial stress and different confining pressures因此，改变围压和轴向应力，均会造成数字岩心孔隙半径喉道半径剧烈变化，其最大值减小，小于体素大小的孔隙和喉道闭合；喉道长度随着围压和轴向应力增大而增大，导致孔隙之间的传导距离变长，从而造成比欧拉示性数曲线与直线y = 0交点值变小，即数字岩心连通性变差。
3.3 不同应力数字岩心孔隙度和渗透率演化分析
不同应力组合作用下数字岩心的孔隙度和渗透率柱状图如图11所示。从图11（a）可以看出，围压和轴向应力增大，均会造成孔隙度降低，最大降低幅度为6.21%；结合3.2节的孔隙和喉道尺寸分析结果，孔隙度降低的主要原因是应力作用下小于体素尺寸孔隙和喉道的闭合，以及大于体素尺寸孔隙和喉道半径的减小，导致孔隙和喉道的数量减少和小半径孔隙喉道的比例增加，最终造成数字岩心孔隙度降低。从图11（b）可以看出，围压和轴向应力增大，同样会导致数字岩心渗透率降低，常温常压下数字岩心模型的渗透率为31.50 mD，90/110/180 MPa应力组合作用下数字岩心的渗透率为25.96 mD，降幅为17.61%。
![]()
图 11 不同应力组合作用下数字岩心孔隙度和渗透率的分布Figure 11. Porosity and permeability distribution of digital core under different stress combinationsYang Yongfei等人[33]实验研究 表明，围压作用下岩石孔隙度的降低幅度为12.21%，渗透率的降低幅度为31.0%，即渗透率的应力敏感性远超孔隙度。模拟结果表明，孔隙度最大降幅为6.21%，渗透率为17.61%，符合渗透率应力敏感性大于孔隙度应力这一实验规律。Lu Jun等人[34]研究表明，真三轴应力下深部储层岩石的渗透率随应力增大的最大降幅（偏应力为60 MPa时）为18.57%，模拟结果为17.61%，相对误差为4.85%。本文模拟结果和文献中实验结果的对比表明，模拟结果能够很好地反映深部岩石物性在应力作用下的变化规律。
4. 结论与建议
1）基于离散元法，提出了一种考虑深层超深层油气藏高应力影响的数字岩心重构方法。该方法可以构建出任意应力下的数字岩心。离散元建模时，采用球面谐波分析方法，并考虑了颗粒的实际形状，提升了模型的准确性，构建的数字岩心更能反映实际岩石的孔隙结构。
2）以Bentheim砂岩为例，构建了不同应力作用下的数字岩心，孔隙几何拓扑结构分析表明，随着围压和轴向应力增大，颗粒向模型内部位移量增大，孔隙和喉道被压缩，喉道长度拉长，连通性变差，最终导致数字岩心的孔隙度和渗透率降低；对比文献中的实验结果，模拟结果能精确反映应力对孔隙度和渗透率的影响。
3）应用该方法重构数字岩心时，必须同时获取深部储层岩石在常温常压下的CT扫描图像和应力应变曲线，用于离散元建模与应力加载模拟，最终生成深层超深层油气藏的数字岩心。在这个过程中，岩心图像向颗粒的转化是关键步骤，对图像的分辨率要求极高。因此，建议继续研究更高效和普适的图像颗粒化方法。
-
![]()
图 1 低渗透地层中流体流速与压力梯度的关系
Figure 1. Relationship between fluid flow velocity and pressure gradient in low permeability formation
![]()
图 4 压裂注采井组压力和饱和度的计算结果
Figure 4. Calculation results of the pressure and saturation of fractured injection-production well group
![]()
图 5 注水井压裂对注采井组产量的影响
Figure 5. Effect of fracturing in water injection well on the production of injection-production well group
-
[1] 宋付权,刘慈群. 含启动压力梯度油藏的两相渗流分析[J]. 石油大学学报(自然科学版), 1999, 23(3): 47–50. doi: 10.3863/j.issn.1674-5086.1999.03.014 SONG Fuquan, LIU Ciqun. Analysis of two-phase fluid flow in low permeability reservoir with the threshold pressure gradient[J]. Journal of the University of Petroleum, China(Edition of Natural Science), 1999, 23(3): 47–50. doi: 10.3863/j.issn.1674-5086.1999.03.014
[2] 程时清, 陈明卓. 油水两相低速非达西渗流数值模拟[J]. 石油勘探与开发, 1998, 25(1): 41–43. doi: 10.3321/j.issn:1000-0747.1998.01.012 CHENG Shiqing, CHEN Mingzhuo. Numerical simulation of oil-water low-velocity non-Darcy flow[J]. Petroleum Exploration and Development, 1998, 25(1): 41–43. doi: 10.3321/j.issn:1000-0747.1998.01.012
[3] 周涌沂,彭仕宓,李允,等. 低速非达西渗流的全隐式模拟模型[J]. 石油勘探与开发, 2002, 29(2): 90–93. doi: 10.3321/j.issn:1000-0747.2002.02.024 ZHOU Yongyi, PENG Shimi, LI Yun, et al. Fully implicit simulation model for low-velocity non-Darcy flow[J]. Petroleum Exploration and Development, 2002, 29(2): 90–93. doi: 10.3321/j.issn:1000-0747.2002.02.024
[4] 尹芝林,孙文静,姚军. 动态渗透率三维油水两相低渗透油藏数值模拟[J]. 石油学报, 2011, 32(1): 117–121. doi: 10.3969/j.issn.1001-8719.2011.01.020 YIN Zhilin, SUN Wenjing, YAO Jun. Numerical simulation of the 3D oil-water phase dynamic permeability for low-permeability reservoirs[J]. Acta Petrolei Sinica, 2011, 32(1): 117–121. doi: 10.3969/j.issn.1001-8719.2011.01.020
[5] 吕广忠,鞠斌山,栾志安. 油藏水力压裂区域分解模拟算法[J]. 石油大学学报(自然科学版), 1998, 22(5): 61–63. LYU Guangzhong, JU Binshan, LUAN Zhian. Domain decomposition simulation method for hydraulic fracturing area of reservoir[J]. Journal of the University of Petroleum, China (Edition of Natural Science), 1998, 22(5): 61–63.
[6] 苏玉亮,王霞,李涛,等. 人工裂缝对低渗透油田开发的影响研究[J]. 钻采工艺, 2006, 29(4): 33–34. doi: 10.3969/j.issn.1000-7393.2006.04.011 SU Yuliang, WANG Xia, LI Tao, et al. Influence of created fracture on low permeability reservoir development[J]. Drilling & Production Technology, 2006, 29(4): 33–34. doi: 10.3969/j.issn.1000-7393.2006.04.011
[7] 何勇明,孙尚如,徐荣伍, 等. 低渗透油藏污染井压裂增产率预测模型及敏感性分析[J]. 中国石油大学学报(自然科学版), 2010, 34(3): 76–79. doi: 10.3969/j.issn.1673-5005.2010.03.016 HE Yongming, SUN Shangru, XU Rongwu, et al. Prediction model for fracturing incremental recovery of damaged well in low-permeability reservoir and sensitivity analysis[J]. Journal of China University of Petroleum(Edition of Natural Science), 2010, 34(3): 76–79. doi: 10.3969/j.issn.1673-5005.2010.03.016
[8] BELHAJ H A, AGHA K R, NOURI A M, et al. Numerical modeling of Forchheimer’s equation to describe darcy and non-Darcy flow in porous media[R]. SPE 80440, 2003.
[9] SOLIMAN M Y. Numerical model estimates fracture production increase[J]. Oil and Gas, 1986, 84(41): 70–74.
[10] 温庆志,张士诚,王秀宇,等. 支撑裂缝长期导流能力数值计算[J]. 石油钻采工艺, 2005, 27(4): 68–70. doi: 10.3969/j.issn.1000-7393.2005.04.020 WEN Qingzhi, ZHANG Shicheng, WANG Xiuyu, et al. Numerical calculation of long - term conductivity of propping fractures[J]. Oil Drilling & Production Technology, 2005, 27(4): 68–70. doi: 10.3969/j.issn.1000-7393.2005.04.020
[11] 任勇,郭建春,赵金洲,等. 压裂井裂缝导流能力研究[J]. 河南石油, 2005, 19(1): 46–48. doi: 10.3969/j.issn.1673-8217.2005.01.017 REN Yong, GUO Jianchun, ZHAO Jinzhou, et al. A study on flow conductivity of fractures in a fractured well[J]. Henan Petroleum, 2005, 19(1): 46–48. doi: 10.3969/j.issn.1673-8217.2005.01.017
[12] 胥元刚,张琪. 变裂缝导流能力下水力压裂整体优化设计方法[J]. 大庆石油地质与开发, 2000, 19(2): 40–43. doi: 10.3969/j.issn.1000-3754.2000.02.014 XU Yuangang, ZHANG Qi. Overall optimizing designation method for hydraulic fracturing under variable fracture diverting capacity[J]. Petroleum Geology & Oilfield Development in Daqing, 2000, 19(2): 40–43. doi: 10.3969/j.issn.1000-3754.2000.02.014
[13] 孔祥言．高等渗流力学[M]．合肥: 中国科学技术大学出版社, 1999: 76–77. KONG Xiangyan. Advanced seepage mechanics[M]. Hefei: Press of University of Science and Technology of China, 1999: 76–77.
[14] 李淑霞, 谷建伟．油藏数值模拟基础[M]．东营: 中国石油大学出版社, 2009: 97–101． LI Shuxia, GU Jianwei. Fundamentals of numerical reservoir simulation[M]. Dongying: China University of Petroleum Press, 2008: 97–101.
[15] 戴嘉尊, 邱建贤．微分方程数值解法[M]．南京: 东南大学出版社, 2002． DAI Jiazun, QIU Jianxian. Numerical solutions for differential equations[M]. Nanjing: Southeast University Press, 2002.
[16] 张建国, 雷光伦．油气层渗流力学[M]．东营: 石油大学出版社, 1998: 46–47． ZHANG Jianguo, LEI Guanglun. Seepage mechanics of oil and gas reservoir[M]. Dongying: Petroleum University Press, 1998: 46–47.
-
期刊类型引用(9)
1. 伍藏原，张迅，李杨，姚杰，冯其红. 塔中深层稠油减氧空气驱氧化特征及机理研究. 西南石油大学学报(自然科学版). 2024(04): 131-137 . 
百度学术
2. 佘治成，陈利新，徐三峰，肖云，张键. 缝洞型碳酸盐岩油藏注氮气提高采收率技术. 西南石油大学学报(自然科学版). 2024(04): 138-148 . 百度学术
3. 王林，岳鹏，安高诧. 塔河油田油田水除氧剂的筛选与应用. 内蒙古石油化工. 2024(10): 36-39 . 百度学术
4. 牟汉生，陆文明，曹长霄，宋兆杰，石军太，张洪. 深水浊积岩油藏提高采收率方法研究——以安哥拉X油藏为例. 石油钻探技术. 2021(02): 79-89 . 本站查看
5. 刘中云，李兆敏，赵海洋. 缝洞型碳酸盐岩油藏注氮气致稠机理研究. 石油钻探技术. 2021(05): 75-80 . 本站查看
6. 李海波，李永会，姜瑜，陈凤，吕世瑶，展宏洋. 多层砂砾岩油藏注烟道气提高采收率技术. 断块油气田. 2020(01): 104-108 . 百度学术
7. 张伟，海刚，张莹. 塔河油田碳酸盐岩缝洞型油藏气水复合驱技术. 石油钻探技术. 2020(01): 61-65 . 本站查看
8. 贾国振. 略谈油田氮气开采技术的运用. 化工管理. 2020(07): 57-58 . 百度学术
9. 李二党，韩作为，高祥瑞，马明宇，邱钧超. 不同注气介质驱替致密油藏微观孔隙动用特征研究. 石油钻探技术. 2020(05): 85-91 . 本站查看
其他类型引用(1)
下载:
计量
- 文章访问数: 1214
- HTML全文浏览量: 557
- PDF下载量: 52
- 被引次数: 10
