一种自支撑裂缝导流能力的数值计算方法与流程

本发明涉及石油领域，尤其是页岩清水压裂过程中一种自支撑裂缝导流能力的数值计算方法。

背景技术：

水力压裂技术是低渗透油气藏增产改造的重要措施。水力压裂是利用地面高压泵组，以超过地层吸收能力的排量将压裂液泵入地层来产生裂缝，然后继续注入带有支撑剂(砂粒)的压裂液，使裂缝继续延伸并在其中充填支撑剂，当压裂液返排后，在地层压力作用下，支撑剂在裂缝中起到支撑裂缝的作用，阻止裂缝闭合，从而在地层中形成具有一定长度和导流能力的填砂裂缝。

清水压裂是水力压裂的一种形式，被广泛应用于页岩油气藏的增产改造中。它具体是指在压裂过程中不加入支撑剂(砂粒)，仅通过将低粘度的压裂液泵入地层来产生裂缝。地下岩石性质差异较大，压裂形成的裂缝表面一般凹凸不平，同时还会在剪切作用下发生错位，因此即使不加入支撑剂，裂缝表面的凸点之间也可以相互支撑，形成自支撑裂缝，自支撑裂缝在地层压力作用下仍能保持一定开启程度，并具有一定的导流能力，从而达到改善油气流动条件和油气井增产的目的。因此，自支撑裂缝的导流能力是评价水力压裂(清水压裂)成功与否的重要指标之一。

技术实现要素：

本发明的目的在于提供一种自支撑裂缝导流能力的数值计算方法，该方法将目标储层页岩加工为长方体岩样，并劈裂分为两个具有粗糙表面的岩板，将两个岩板的粗糙表面错位放置形成错位裂缝，通过对错位裂缝进行变形分析和数值建模，最终计算出裂缝导流能力，该方法原理可靠，操作简便，能够简单、准确地计算自支撑裂缝的导流能力，对水力压裂具有重要的指导意义。

为达到以上技术目的，本发明采用以下技术方案：对自支撑错位裂缝的受力形态进行建模，绘制裂缝的“应力-位移”图版，基于图版计算裂缝受力变形后的裂缝宽度数据，再基于格子玻尔兹曼方程计算裂缝中的流体流速分布及密度分布数据，进而得到裂缝中的流体流量和压力数据，最终利用达西公式计算出裂缝导流能力。

利用三维激光扫描仪按专利cn201510319382.9所述的方法对岩板裂缝面进行扫描，获取裂缝表面的三维数据，并建立岩板的三维坐标系：取岩板长度方向为横轴方向，宽度方向为纵轴方向，高度方向为竖轴方向，且横轴和纵轴方向所在的面为水平面。

一种自支撑裂缝导流能力的数值计算方法，依次包括以下步骤：

(a)对目标储层裂缝面进行单元离散分析，建立裂缝受力变形模型，绘制裂缝应力-位移图版；

(b)利用步骤(a)绘制的应力-位移图版确定目标储层裂缝在已知闭合应力σf作用下的裂缝变形量zf；

(c)对目标储层岩板表面进行激光扫描，基于步骤(b)计算裂缝宽度矩阵wf；

(d)根据格子玻尔兹曼方法，利用下式计算裂缝内流体粒子微团的流速、密度及压力数据：

式中：α—流体粒子微团的规定运动方向；α＝0，1，2……18；

—流体粒子微团在空间和时间维度上的分布函数；

—流体粒子微团在α方向上的离散速度；

δt—离散时间步长；

τ0—松弛时间，一般计算为τ0＝3μ+0.5，μ为流体粘度；

—流体粒子微团在空间位置为时间为t+δt时刻的分布函数；

—流体粒子微团在α方向上，空间位置为时间为t的平衡分布函数；

ρ—时间为t时，流体粒子微团密度；

—时间为t时，流体粒子微团流动速度；

p—时间为t时，流体粒子微团的压力；

(e)基于步骤(d)，利用下式计算裂缝内流体流量q，裂缝入口压力p1和裂缝出口压力p2：

式中：wfin—裂缝入口处的裂缝宽度矩阵；

ufx—流体粒子微团在裂缝长度方向上的流速分布矩阵；

x—对目标储层裂缝面进行单元离散分析时，离散单元的边长；

wfinj—矩阵wfin中的元素；

ufxj—矩阵ufx中的元素；

pinj—裂缝入口处的压力分布矩阵pin中的元素；

poutj—裂缝出口处的压力分布矩阵pout中的元素；

k—对裂缝面进行扫描后，纵轴方向共有k列个扫描数点；

(f)使用以下达西公式计算裂缝导流能力：

式中：fcd—裂缝导流能力；

μ—流体粘度；

l—裂缝(岩板)长度；

b—裂缝(岩板)宽度。

本发明中，所述步骤(a)对目标储层裂缝面进行单元离散分析，建立裂缝受力变形模型，绘制裂缝应力-位移图版，过程如下：

对扫描裂缝面三维数据进行坐标变换，使两个裂缝表面刚好接触，对上裂缝面施加向下的初始位移量z0，在竖轴方向上选定一个离散单元，离散单元截面为正方形，边长为x，上裂缝面高度z1，竖轴方向上所受的压力为δfz，压缩量为δz1，在竖轴方向所受的应力值为：

式中：σm为裂缝岩体的抗压强度；

mc为岩体的应力突变系数(本发明中，推荐页岩应力突变系数mc取值0.5)；

υ为裂缝岩体的泊松比。

从而求得上裂缝面所受的应力值σ：

式中：m—对裂缝面进行立方单元离散后，横轴方向共有m排离散单元；

n—对裂缝面进行立方单元离散后，纵轴方向共有n列离散单元；

δσi,j—横轴方向第i排，纵轴方向第j列的离散单元所受的应力值。

以z表示上裂缝面的向下位移量，以σ表示上裂缝面所受的应力值，通过对上裂缝面施加不同的向下位移量，不断重复上述步骤，直到获得一定数量的(σ，z)数据点，绘制σ-z图版，得到z＝f(σ)曲线。

本发明中，所述步骤(b)利用步骤(a)绘制的应力-位移图版确定目标储层裂缝在已知闭合应力σf作用下的裂缝变形量zf，过程如下：在步骤(a)中，上裂缝面的向下位移量即为实际裂缝的变形量，所以，通过σ-z图版，给定任意应力值σ，利用z＝f(σ)曲线即可得到相应的裂缝变形量z，从而可得到目标储层裂缝在已知闭合应力σf作用下的裂缝变形量zf。

本发明中，步骤(c)对目标储层岩板表面进行激光扫描，基于步骤(b)计算裂缝宽度矩阵wf，过程如下：对扫描裂缝面三维数据进行坐标变换，使两个裂缝表面刚好接触，下裂缝面表面每个点的竖轴数据组成下裂缝面高度矩阵zd，上裂缝面表面每个点的竖轴数据组成上裂缝面高度矩阵zu。通过步骤(b)已确定了在闭合应力σf作用下，裂缝变形量为zf，即上裂缝面的向下位移量为zf。则在σf作用下，上裂缝面高度矩阵变为z′u：

z′u＝(zu′ij)＝(zuij-zf)，(i＝1,2,3,...,h；j＝1,2,3,...,k)(3)

式中：z′uij—z′u矩阵中，横轴方向第i排，纵轴方向第j列的点的高度值；

zuij—zu矩阵中，横轴方向第i排，纵轴方向第j列的点的高度值。

获得裂缝形态变形矩阵zc：

zc＝z′u-zd(4)

从而得到裂缝宽度矩阵wf为：

式中：wfij为裂缝中横轴方向第i排，纵轴方向第j列的点的裂缝宽度值；

zcij为裂缝形态变形矩阵zc矩阵中横轴方向第i排，纵轴方向第j列的点的高度值。

本发明中，步骤(d)根据格子玻尔兹曼方法，计算裂缝内流体流速、密度及压力数据，过程如下：

玻尔兹曼(boltzmann)方程(杨大勇.格子玻尔兹曼方法[m].电子工业出版社.2015)是描述粒子微团速度分布函数(简称“分布函数”)在时间和空间上变化的守恒方程，裂缝中的流体流动满足离散格子boltzmann方程耦合bgk方程(简称lbgk方程)：

本发明选用常见的三维格子玻尔兹曼速度模型d3q19单速模型进行计算，该模型具有离散速度的对称性，粒子演化机制简单，是目前广泛使用的三维格子玻尔兹曼模型。

在d3q19单速模型条件下，式(6)中：

α—流体粒子微团的规定运动方向；α＝0，1，2……18；

—流体粒子微团在空间和时间维度上的分布函数；

—流体粒子微团在α方向上的离散速度；

δt—离散时间步长；

τ0—松弛时间，一般计算为τ0＝3μ+0.5，μ为流体粘度；

—流体粒子微团在空间位置为时间为t+δt时刻的分布函数；

—流体粒子微团在α方向上，空间位置为时间为t的平衡分布函数。

式(6)描述了流体粒子微团的微观运动，通过对式(6)进行求解可获得流体粒子微团在某一位置上的分布函数求得分布函数之后，可以利用分布函数，根据以下三个公式求得流体粒子微团的密度、流速和压力：

式中：ρ—时间为t时，流体粒子微团密度；

—时间为t时，流体粒子微团流动速度，

—时间为t时，流体粒子微团在裂缝长度方向上的流速；

—时间为t时，流体粒子微团在裂缝宽度方向上的流速；

p—时间为t时，流体粒子微团的压力。

要对式(6)进行求解，首先需要求得式(6)中的平衡分布函数其表统一表达式为：

式中：ωα—权重系数；

cs—格子声速。

在d3q19单速模型格子模型中：

式(12)中，c为声速，即340m/s。

式(8)与式(11)中，取在横轴方向的分向量代入式(8)计算得到流体粒子微团在裂缝长度方向上的流速同理，取在纵轴方向的分向量代入式(8)计算得到流体粒子微团在裂缝宽度方向上的流速

对裂缝表面进行网格划分，在裂缝长度方向上划分h(i＝1,2,3,...,h)个节点，裂缝宽度方向上划分k个节点，保证网格划分数量与三维激光扫描仪扫描后得到的网格数量一致。在每个网格的节点处即是一个流体微团，利用计算机编程在离散时间步长为δt条件下迭代计算。根据格子玻尔兹曼发展经验，设定初始状态(t＝0时刻)下，裂缝内流体密度ρt＝0为1g/cm³，此时的分布函数计算式为：

根据式(8)计算得到每个节点处流体微团流动速度并根据式(10)计算t＝0时刻的平衡分布函数：

t＝0起，每个流体微团开始运动，采用定流量边界条件(lian-pingwang,andchengpeng.latticeboltzmannsimulationofparticle-ladenturbulentchannelflow[j].computersandfluids.2015:1-11.)计算t＝δt时刻的流体微团分布函数：

利用结合式(7)、式(8)、式(9)计算得到t＝δt时刻的密度分布数据速度分布数据以及压力分布数据根据式(10)计算t＝δt时刻的平衡分布函数，以此开始第二次迭代。以此反复迭代计算，直到t＝nδt时刻(n为自然数)达到式(17)的收敛条件：

式中：—时间为t＝(n-1)δt时，流体粒子微团在裂缝长度方向上的流速；

—时间为t＝(n-1)δt时，流体粒子微团在裂缝宽度方向上的流速。

本发明中，所述步骤(e)计算裂缝内流体流量q，裂缝入口压力p1和裂缝出口压力p2，过程如下：

获得裂缝宽度矩阵wf后，在矩阵wf中取i＝1条件下的所有元素组成裂缝入口处的裂缝宽度矩阵wfin。基于达到收敛条件的流体微团分布函数可计算出网格数i＝1条件下，所有节点处流体粒子微团在裂缝长度方向上的流速分布数据并组成矩阵ufx，矩阵ufx即为裂缝入口处流体粒子微团在裂缝长度方向上的流速分布矩阵。则裂缝内流量q为：

式中：wfinj—矩阵wfin中的元素；

ufxj—矩阵ufx中的元素；

x—离散单元截面边长。

基于达到收敛条件的流体微团分布函数可计算出网格数i＝1条件下，所有节点处流体粒子微团的压力数据并组成矩阵pin，矩阵pin即为裂缝入口处流体粒子微团的压力分布矩阵。同理，基于达到收敛条件的流体微团分布函数可计算出网格数i＝h条件下，所有节点处流体粒子微团的压力数据并组成矩阵pout，矩阵pout即为裂缝出口处流体粒子微团的压力分布矩阵。采用求平均值方法得到裂缝入口处力p1和出口处压力p2：

式中：pinj—矩阵pin中的元素；

poutj—矩阵pout中的元素。

附图说明

图1为裂缝空间直角坐标系的任意一个横-竖轴剖面示意图。

图2为选定离散单元受力分析示意图。

具体实施方式

下面根据附图进一步说明本发明。

本发明中，对错位裂缝受力变形进行分析并建立模型，具体过程为：

如图1所示，对扫描裂缝面三维数据进行坐标变换，将下裂缝面原始数中所有点的高度值减去下裂缝面的最低点高度值，即：

zd＝(zdij)＝[z0ij-min(z0ij)],(i＝1,2,3,...,h；j＝1,2,3,...,k)(20)

上式中：zd—下裂缝面高度矩阵；

zdij—下裂缝面高度矩阵中，横轴方向第i排，纵轴方向第j列的点的高度值；

z0ij—原始扫描数据度矩阵中，横轴方向第i排，纵轴方向第j列的点的高度值；

h—对裂缝面进行扫描后，横轴方向共有h列个扫描数点；

k—对裂缝面进行扫描后，纵轴方向共有k列个扫描数点；

min(z0ij)—原始扫描数据高度矩阵中最低点的高度值。

同样，将上裂缝面原始数中所有点的高度值减去上裂缝面的最低点高度值，再将处理后上裂缝面数据点经过坐标变换按图1所示放入下裂缝面的坐标系中，使两个裂缝表面刚好接触。在图1所示的坐标系中，下裂缝面表面每个点的竖轴数据即组成下裂缝面高度矩阵zd，上裂缝面表面每个点的竖轴数据即组成上裂缝面高度矩阵zu。

保持下裂缝面不变，对上裂缝面施加向下的初始位移量z0后，上裂缝面高度矩阵变为z′u：

z′u＝(z′uij)＝(zuij-z0)，(i＝1,2,3,...,h；j＝1,2,3,...,k)(21)

上式中：z′uij—z′u矩阵中，横轴方向第i排，纵轴方向第j列的点的高度值；

zuij—zu矩阵中，横轴方向第i排，纵轴方向第j列的点的高度值；

利用式(4)计算获得裂缝形态变形矩阵zc

在矩阵zc中，当元素zcij(zc矩阵中，横轴方向第i排，纵轴方向第j列的点的高度值)大于零时，表示该点在位移量为z0时，两个裂缝面没有发生接触；当zcij小于或等于零时，表示该点在位移量为z0时两个裂缝面发生了接触变形，且压缩量为|zcij|。

定义上裂缝面与下裂缝面的总压缩量矩阵为zt，矩阵中的元素为ztij：

在竖轴方向上对裂缝面进行立方单元离散，选定一个离散单元，如图2所示，已知：上裂缝面高度z1，下裂缝面高度z2，并通过矩阵zt可得到上裂缝面与下裂缝面的总压缩量δz，离散单元截面为正方形，边长为x，裂缝岩体杨氏模量为e，泊松比为υ，抗压强度为σm。

对上裂缝面离散单元进行受力分析，假设为纯弹性应变，其竖轴方向上所受的压力为δfz，应力为δσ，压缩量为δz1，应变为εz1，根据胡克定律有：

假设上裂缝面离散单元横向变形为δx1，则变形后的离散单元在竖轴方向受力面积为δa，则有：

将公式(24)代入公式(23)中，得到上裂缝面离散单元的受力变形方程：

同理，假设下裂缝面离散单元横向变形为δx2，可推出下裂缝面离散单元的受力变形方程为：

设离散单元上裂缝面在横轴方向的应变为εx1，根据泊松比定义得：

由于上、下粗糙面为同一性质材料，因此有：

离散单元上裂缝面与下裂缝面的总压缩量为δz：

δz1+δz2＝δz(29)

由公式(25)、(26)、(28)、(29)可建立如下计算方程组：

上述方程组含五个方程，可求解得到五个未知数(δz1、δz2、δx1、δx2、δfz)。

由公式(28)可得到δx1与δz1关系以及δx2与δz2关系

将公式(31)代入公式(25)中

上式可整理为：

同理可得：

由式(34)、式(35)可知，求解出δz1或δz2即可计算获得微元受力δfz值。

因此，联立公式(34)、(35)、(29)可得δz1计算表达式：

上式可整理出公式(37)如下：

上式可使用牛顿迭代法进行数值求解得到选定离散单元中上裂缝面的压缩量δz1。迭代过程中由物理对象客观条件即计算值为正且变形量不超过计算边界，设定迭代起始值及迭代精度。本发明推荐迭代起始值为1，迭代精度10^-8。

至此，可利用公式(34)计算出选定离散单元中上裂缝面所受的压力δfz。

公式(37)和公式(34)的推导过程假设了离散单位为纯弹性应变，实际上，随着岩体所受应力的增加，当岩体达到抗压强度后，岩样受到瞬间破坏，应力不再遵循纯弹性应变规律，故需要进一步判断选定离散单元所受的应力是否达到抗压强度：

通过公式(37)和公式(34)可计算获得δfz、δx1，则选定离散单元在竖轴方向上所受的应力值δσ为：

将公式(31)代入公式(38)中得：

裂缝岩体的抗压强度为当σm时，当δσ＜σm时，岩体发生线弹性变形，应力值为当σ0≥σm时，页岩发生应力损伤破坏，应力值为σmmc。此处，mc为岩体的应力突变系数，该系数定义为：岩体受力超过岩体抗压强度σm后，岩体发生瞬间破裂失效，并保持某一残余应力值，失效后的残余应力值与抗压强度的比值为应力突变系数mc。故选定离散单元在竖轴方向所受的应力值为：

本发明在实验室采用大量页岩样本进行应力—位移测试，得到的应力突变系数分布在0.4-0.6范围内。因此本发明推荐页岩应力突变系数mc取值0.5。

可利用求和公式计算上裂缝面所受的应力值：

上式中：σ—上裂缝面所受的应力值；

m—对裂缝面进行立方单元离散后，横轴方向共有m排离散单元；

n—对裂缝面进行立方单元离散后，纵轴方向共有n列离散单元；

δσi,j—横轴方向第i排，纵轴方向第j列的离散单元所受的应力值；

此时，已经计算出了在一个向下的初始位移量z0条件下，上裂缝面所受的应力值σ，现给定一个位移量步长t，将初始位移量增加t，继续计算上裂缝面所受的应力值。

以z表示上裂缝面的位移量，以σ表示上裂缝面所受的应力值，不断重复上述步骤，直到获得一定数量的(σ，z)数据点，利用(σ，z)数据点绘制σ-z图版，即可得到z＝f(σ)曲线，在σ-z图版中，给定任意应力值σ，利用z＝f(σ)曲线得到相应的裂缝变形量z。

自支撑裂缝的宽度一般是毫米级，为了尽可能获取更多的数据点，本发明建议初始位移量取值0.1mm，位移量步长t取值0.1mm。

一种自支撑裂缝导流能力的数值计算方法，依次包括以下步骤：

(a)对目标储层裂缝面进行单元离散分析，建立裂缝受力变形模型，绘制裂缝应力-位移图版；

(b)利用步骤(a)绘制的应力-位移图版确定目标储层裂缝在已知闭合应力σf作用下的裂缝变形量zf；

(c)对目标储层岩板表面进行激光扫描，基于步骤(b)计算裂缝宽度矩阵wf；

(d)根据格子玻尔兹曼方法，计算裂缝内流体流速、密度及压力数据；

(e)基于步骤(d)，计算裂缝内流体流量q，裂缝入口压力p1和裂缝出口压力p2；

(f)使用以下公式计算裂缝导流能力：

式中：fcd—裂缝导流能力；

μ—流体粘度；

l—裂缝(岩板)长度；

b—裂缝(岩板)宽度。