一种通过耦合储层流动确定水平井压裂改造体积的方法

1.本发明涉及油气藏开发技术领域，特别涉及一种通过耦合储层流动确定水平井压裂改造体积的方法。


背景技术：

2.目前，体积压裂已经成为了页岩气资源高效开发中不可或缺的技术手段。压裂改造体积的大小直接决定着压裂改造效果和单井可采储量(eur)的评估，因此确定水平井压裂改造体积对提高油气资源的开发效果具有重要意义。
3.在页岩气井的压裂施工过程中，压裂液会不断经过基质和裂缝滤失进入储层，同时地层压力也会随之向储层深部扩散，这会导致裂缝周围储层的孔隙结构发生改变，从而扩大压裂改造体积。目前，现有的水平井压裂改造体积确定方法有很多，但这些方法都忽略了上述作用，导致压裂液在储层中的波及范围无法计算，储层压力分布、含水饱和度分布无法获取等等，因此最终获取的压裂改造体积结果不够准确。


技术实现要素：

4.针对上述问题，本发明旨在提供一种通过耦合储层流动确定水平井压裂改造体积的方法，将压裂过程中储层流体的流动纳入考虑范畴，在获取储层含水饱和度场和压裂液滤失情况的基础之上，分析确定压裂改造体积，获取更准确的压裂改造体积结果。
5.本发明的技术方案如下：
6.一种通过耦合储层流动确定水平井压裂改造体积的方法，包括以下步骤：
7.s1：获取目标井的储层地质模型，根据所述储层地质模型建立目标井的储层网格，并在所述储层网格中添加各簇裂缝的初始裂缝单元，所述初始裂缝单元被所述储层网格划分为多个裂缝单元；
8.所述储层网格中水力裂缝的端点总数为u，裂缝端点编号为e，则e＝1,2,
…
,u-1,u；e为奇数时表示裂缝上端点，e为偶数时表示裂缝下端点；设裂缝单元总数为n
l
，裂缝单元编号为l，则l＝1,2,
…
,n
l-1,n
l
；设第l个裂缝单元长度为ξ
l
，第l个裂缝单元内的流体压力为p
f,l
，第l个裂缝单元的开启时间为t
f,l
，各簇裂缝的初始裂缝单元的开启时间均为t0；
9.s2：基于边界元位移不连续法，计算t0时间下第e个裂缝端点处的应力强度因子；
10.s3：根据t0时间下第e个裂缝端点处的应力强度因子，判断t0时间下第e个裂缝端点处是否发生裂缝扩展：
11.若第e个裂缝端点处，裂缝不发生扩展，则重复步骤s2-s3，计算t0时间下第e 1个裂缝端点处的应力强度因子，并判断t0时间下第e 1个裂缝端点处是否发生裂缝扩展；
12.若第e个裂缝端点处，裂缝发生扩展，则第e个裂缝端点所在的裂缝沿第e个裂缝端点方向增加增加一个裂缝单元，且裂缝单元总数n
l
＝n
l
1，新增裂缝单元的流体压力与其相邻裂缝单元内的流体压力相同，新增裂缝单元的开启时间重复步骤s2-s3，以新增裂缝单元后的数据为基础，计算t0时间下第e 1个裂缝端点处的应力强度因子，并判断t0时
间下第e 1个裂缝端点处是否发生裂缝扩展；
13.当所有裂缝端点处的裂缝扩展情况判断结束后，即可获得t0时间下裂缝扩展之后的裂缝单元总数裂缝单元编号t0时间下第l个裂缝单元的开启时间t0时间下第l个裂缝单元内的流体压力
14.s4：以步骤s3获得的数据为基础，结合嵌入式离散裂缝，运用气-水两相双重介质渗流模型进行数值计算，获得压裂施工过程中t1时间下的水力裂缝单元内的流体压力、基质系统的压力分布、基质系统的含水饱和度分布、微裂缝系统的压力分布、以及微裂缝系统的含水饱和度分布；
15.s5：以步骤s4获得的数据为基础，重复步骤s2-s4，计算t2时间下的水力裂缝单元内的流体压力、基质系统的压力分布、基质系统的含水饱和度分布、微裂缝系统的压力分布、以及微裂缝系统的含水饱和度分布；
16.s6：重复步骤s5，直至时间达到压裂施工时间t
end
，获得t
end
时间下的水力裂缝单元数、基质系统的含水饱和度分布、微裂缝系统的含水饱和度分布；
17.s7：根据步骤s6获得的t
end
时间下的数据，计算压裂改造体积。
18.作为优选，步骤s1中，所述储层网格为x-y坐标系下网格数为ni×
nj的结构化网格，所述结构化网格中的x
i,j
和y
i,j
分别代表每个网格的长度和宽度，下角标i和j代表每一个网格在储层中的位置；每个网格中存在基质系统和微裂缝系统，以此构成双重介质模型；所述基质系统和所述微裂缝系统的初始压力均为原始地层压力，所述基质系统和所述微裂缝系统的含水饱和度均为原始地层含水饱和度。
19.作为优选，步骤s1中，在所述储层网格中添加各簇裂缝的初始裂缝单元时：每簇裂缝的初始裂缝单元方向为y轴方向，每簇裂缝的初始裂缝单元的长度为n个网格的长度，所述n为大于等于3的整数。
20.作为优选，步骤s2具体包括以下子步骤：
21.s21：假设每个裂缝单元内流体压力相同，则在t0时间下裂缝端点处的流体压力与裂缝端点所处的裂缝单元内的流体压力相等；
22.s22：根据t0时间下的裂缝单元位置，确定x-y坐标系下每簇裂缝的上下端点位置坐标，分别表示为(x1,y1),
…
,(xe,ye)；
23.s23：以第1个裂缝单元的中心为原点，第1个裂缝单元的延伸方向为方向，第1个裂缝单元的延伸方向的垂直方向为方向，建立局部坐标系，并将步骤s22中所有的裂缝端点位置坐标转换为所述局部坐标系下的位置坐标，分别表示为
24.s24：计算t0时间下第1、第2，
……
，第n
l
个裂缝单元在第e个裂缝端点处产生的法向不连续位移；
25.s25：将步骤s24计算得到的第e个裂缝端点处所有的法向不连续位移进行叠加，以叠加后的法向不连续位移为基础，采用裂缝尖端应力强度因子计算模型计算t0时间下第e个裂缝端点处的应力强度因子。
26.作为优选，步骤s24中，t0时间下第1个裂缝单元在第e个裂缝端点处形成的所述法
向不连续位移通过下式进行计算：
[0027][0028][0029][0030]
式中：为t0时间下第1个裂缝单元在第e个裂缝端点处产生的法向不连续位移，mm；为t0时间下第e个裂缝端点处的流体压力，mpa；σh为储层的最小水平主应力，mpa；a
xx(1,e)
、a
xy(1,e)
、a
yy(1,e)
、f
a(1,e)
、f
b(1,e)
、f
c(1,e)
均为计算过程中的中间函数；α1为第1个裂缝单元局部坐标系的轴与x-y坐标系的x轴之间的夹角；e为储层岩石的杨氏模量，gpa；υ为储层岩石的泊松比；为第1个裂缝单元局部坐标系下的第e个裂缝端点坐标；ξ1为第1个裂缝单元的长度，m；
[0031]
步骤s25中，t0时间下第2，
……
，第n
l
个裂缝单元在第e个裂缝端点处形成的法向不连续位移的计算方法与t0时间下第1个裂缝单元在第e个裂缝端点处形成的法向不连续位移的计算方法相同。
[0032]
作为优选，步骤s25中，所述裂缝尖端应力强度因子计算模型如下所示：
[0033]
[0034][0035]
式中：为t0时间下第e个裂缝端点处的应力强度因子，re为第e个裂缝端点处裂缝单元的半长，mm；为t0时间下第e个裂缝端点处的总法向不连续位移，mm；为t0时间下第l个裂缝单元在第e个裂缝端点处产生的法向不连续位移，mm。
[0036]
作为优选，步骤s3中，判断t0时间下第e个裂缝端点处是否发生裂缝扩展的具体判断方法如下：
[0037]
将t0时间下的第e个裂缝端点处的应力强度因子与储层岩石的断裂韧性k
ic
进行对比：若则在第e个裂缝端点处，裂缝不发生扩展；若则裂缝会在第e个裂缝端点处发生扩展。
[0038]
作为优选，步骤s4中，所述气-水两相双重介质渗流模型包括：
[0039]
(1)水力裂缝单元向微裂缝滤失模型：
[0040][0041][0042]
式中：q
f-fw
为水力裂缝单元与微裂缝网格间的滤失量，m3；c为滤失系数；t为当前时刻，s；tf为水力裂缝单元开启时间，且kf为微裂缝网格的渗透率，d；lf为水力裂缝单元长度，m，lf＝ξ
l
且h为储层厚度，m；μw为压裂液粘度，mpa
·
s；为微裂缝网格内的点到水力裂缝单元的平均法向距离，m；pf为水力裂缝单元流体压力，mpa，且p
fw
为微裂缝网格的水相压力，mpa；af为微裂缝网格的面积，m2；l
f-f
为微裂缝网格的面积单元到裂缝k的距离，m；
[0043]
(2)水力压裂中的单相流动：
[0044][0045][0046]req
＝0.14[(lf)2 (hf)2]
1/2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0047]
式中：β为单位转换系数；kf为水力裂缝的渗透率，d；bw为压裂液的体积系数；δ
well
为判断水力裂缝单元与井筒相交关系的系数，如果裂缝单元与井筒相交，则δ
well
＝1；不相交，则δ
well
＝0；q
fw
为水力裂缝单元与井筒之间的流量交换，m3；vf为水力裂缝单元的体积，m3；φf为水力裂缝的孔隙度，％；wf为水力裂缝单元宽度，m；p
wf
为井底流压，mpa；r
eq
为有效半径，m；r
well
为井筒半径，m；s为表皮系数；hf为水力裂缝单元高度，m；
[0048]
(3)压裂过程中的基质系统和微裂缝系统的渗流模型：
[0049][0050][0051][0052][0053]
p
mc
＝p
mg-p
mw
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0054]
p
fc
＝p
fg-p
fw
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0055]
式中：为哈密顿算子；k
frw
、k
frg
为微裂缝网格中液相和气相的相对渗透率；δf为微裂缝网格是否含有水力裂缝的判断参数，当微裂缝网格有水力裂缝穿过时，δf＝1；当微裂缝网格无水力裂缝穿过时，δf＝0；vf为微裂缝网格的体积，m3；km为基质渗透率，md；k
mrw
、k
mrg
为微裂缝网格中液相和气相的相对渗透率；p
mw
、p
mg
为基质网格中液相和气相的压力，mpa；φf、φm为微裂缝和基质的孔隙度，％；s
fw
、s
fg
为微裂缝网格中的液相和气相饱和度；μg为气体的粘度，mpa
·
s；bg为气体的体积系数；p
fg
为微裂缝网格的气相压力，mpa；s
mw
、s
mg
为基质网格中的液相和气相饱和度；p
fc
、p
mc
为微裂缝和基质的毛管力，mpa；
[0056]
(4)初始条件：
[0057][0058][0059][0060]
式中：pf(l)为渗流模型计算中第l个裂缝单元的流体压力，mpa；p
fw(i,j)
、p
mw(i,j)
为渗流模型计算中(i,j)网格位置处的微裂缝和基质系统液相压力，mpa；p
fwi,j,t0
、p
mwi,j,t0
为t0时间下(i,j)网格位置处的微裂缝和基质系统液相压力，mpa；s
fw(i,j)
、s
mw(i,j)
为渗流模型计算中(i,j)网格位置处的微裂缝和基质系统含水饱和度；s
fwi,j,t0
、s
mwi,j,t0
为t0时间下(i,j)网格位置处的微裂缝和基质系统含水饱和度；t
f(l)
为渗流模型计算中第l个裂缝单元的开启时间，s；t为当前时间，s。
[0061]
作为优选，步骤s7中，所述压裂改造体积通过下式进行计算：
[0062]
[0063][0064]
式中：v0为压裂改造体积，m3；α
1,i,j
为基质系统中压裂改造范围的判断参数，当为基质系统中压裂改造范围的判断参数，当时，(i,j)网格位置处的基质系统受到改造，α
1,i,j
＝1；当时，α
1,i,j
＝0；h为储层厚度，m；φm为基质的孔隙度，％；α
2,i,j
为微裂缝系统中压裂改造范围的判断参数，当时，(i,j)网格位置处的微裂缝系统受到改造，α
2,i,j
＝1；若时，α
2,i,j
＝0；φf为微裂缝的孔隙度，％；为压裂结束后的水力裂缝长度，m；为t
end
时间下水力裂缝单元数。
[0065]
本发明的有益效果是：
[0066]
本发明结合嵌入式离散裂缝模型、边界元位移不连续法和储层双重介质模型对裂缝扩展、压裂液滤失和储层流体流动过程进行模拟，然后依据压裂波及范围来确定改造体积，本发明在裂缝扩展过程中较好地耦合了储层流动，同时具有较高的计算效率，能够及时高效指导现场施工。通过本发明不仅可以获取压后裂缝参数，还可以耦合储层流动获取压裂全过程中的含水饱和度分布和压力分布演化；进而确定页岩气水平井的压裂改造体积，为页岩气井压裂效果评价、单井eur评估和产能模拟提供依据，促进页岩气资源的高效开发。
附图说明
[0067]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0068]
图1为本发明一个具体实施例中基质的含水饱和度分布结果示意图；
[0069]
图2为本发明一个具体实施例中微裂缝的含水饱和度分布结果示意图。
具体实施方式
[0070]
下面结合附图和实施例对本发明进一步说明。需要说明的是，在不冲突的情况下，本技术中的实施例及实施例中的技术特征可以相互结合。需要指出的是，除非另有指明，本技术使用的所有技术和科学术语具有与本技术所属技术领域的普通技术人员通常理解的相同含义。本发明公开使用的“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。
[0071]
本发明提供一种通过耦合储层流动确定水平井压裂改造体积的方法，包括以下步骤：
[0072]
s1：获取目标井的储层地质模型，根据所述储层地质模型建立目标井的储层网格，并在所述储层网格中添加各簇裂缝的初始裂缝单元，所述初始裂缝单元被所述储层网格划分为多个裂缝单元；
[0073]
所述储层网格中水力裂缝的端点总数为u，裂缝端点编号为e，则e＝1,2,
…
,u-1,u；e为奇数时表示裂缝上端点，e为偶数时表示裂缝下端点；设裂缝单元总数为n
l
，裂缝单元编号为l，则l＝1,2,
…
,n
l-1,n
l
；设第l个裂缝单元长度为ξ
l
，第l个裂缝单元内的流体压力
为p
f,l
，第l个裂缝单元的开启时间为t
f,l
，各簇裂缝的初始裂缝单元的开启时间均为t0；
[0074]
在一个具体的实施例中，所述储层网格为x-y坐标系下网格数为ni×
nj的结构化网格，所述结构化网格中的x
i,j
和y
i,j
分别代表每个网格的长度和宽度，下角标i和j代表每一个网格在储层中的位置；每个网格中存在基质系统和微裂缝系统，以此构成双重介质模型；所述基质系统和所述微裂缝系统的初始压力均为原始地层压力，所述基质系统和所述微裂缝系统的含水饱和度均为原始地层含水饱和度，即：
[0075][0076][0077]
式中：为t0时间下(i,j)网格位置处的微裂缝和基质系统液相压力，mpa；p0为原始地层压力，mpa；s
fwi,j,t0
、s
mwi,j,t0
为t0时间下(i,j)网格位置处的微裂缝和基质系统含水饱和度；s
w,0
为原始地层含水饱和度。
[0078]
基于已建立的储层网格，在水平井压裂施工位置对应的网格处添加初始裂缝单元。在所述储层网格中添加各簇裂缝的初始裂缝单元时：每簇裂缝的初始裂缝单元方向为y轴方向，每簇裂缝的初始裂缝单元的长度为n个网格的长度，所述n为大于等于3的整数。
[0079]
需要说明的是，当n为1或2时，计算精度达不到要求，n越大计算精度越高，为了保证计算效率，在具体应用本发明时，n可直接设置为3。
[0080]
s2：基于边界元位移不连续法，计算t0时间下第e个裂缝端点处的应力强度因子；具体包括以下子步骤：
[0081]
s21：假设每个裂缝单元内流体压力相同，则在t0时间下裂缝端点处的流体压力与裂缝端点所处的裂缝单元内的流体压力相等；即：
[0082][0083]
式中：为t0时间下第e个裂缝端点处的流体压力，mpa；为t0时间下第e个裂缝端点所处的裂缝单元内的流体压力，mpa。
[0084]
s22：根据t0时间下的裂缝单元位置，确定x-y坐标系下每簇裂缝的上下端点位置坐标，分别表示为(x1,y1),
…
,(xe,ye)；
[0085]
s23：以第1个裂缝单元的中心为原点，第1个裂缝单元的延伸方向为方向，第1个裂缝单元的延伸方向的垂直方向为方向，建立局部坐标系，并将步骤s22中所有的裂缝端点位置坐标转换为所述局部坐标系下的位置坐标，分别表示为
[0086]
s24：计算t0时间下第1、第2，
……
，第n
l
个裂缝单元在第e个裂缝端点处产生的法向不连续位移；
[0087]
在一个具体的实施例中，t0时间下第1个裂缝单元在第e个裂缝端点处形成的所述法向不连续位移通过下式进行计算：
[0088][0089][0090][0091]
式中：为t0时间下第1个裂缝单元在第e个裂缝端点处产生的法向不连续位移，mm；为t0时间下第e个裂缝端点处的流体压力，mpa；σh为储层的最小水平主应力，mpa；a
xx(1,e)
、a
xy(1,e)
、a
yy(1,e)
、f
a(1,e)
、f
b(1,e)
、f
c(1,e)
均为计算过程中的中间函数；α1为第1个裂缝单元局部坐标系的轴与x-y坐标系的x轴之间的夹角；e为储层岩石的杨氏模量，gpa；υ为储层岩石的泊松比；为第1个裂缝单元局部坐标系下的第e个裂缝端点坐标；ξ1为第1个裂缝单元的长度，m；
[0092]
t0时间下第2，
……
，第n
l
个裂缝单元在第e个裂缝端点处形成的法向不连续位移的计算方法与t0时间下第1个裂缝单元在第e个裂缝端点处形成的法向不连续位移的计算方法相同，即同样采用公式(1)-(3)计算t0时间下第2，
……
，第n
l
个裂缝单元在第e个裂缝端点处形成的法向不连续位移。
[0093]
s25：将步骤s24计算得到的第e个裂缝端点处所有的法向不连续位移进行叠加，以叠加后的法向不连续位移为基础，采用裂缝尖端应力强度因子计算模型计算t0时间下第e个裂缝端点处的应力强度因子。
[0094]
在一个具体的实施例中，所述裂缝尖端应力强度因子计算模型如下所示：
[0095][0096][0097]
式中：为t0时间下第e个裂缝端点处的应力强度因子，re为第e个裂缝端点处裂缝单元的半长，mm；为t0时间下第e个裂缝端点处的总法向不连续位移，mm；为t0时间下第l个裂缝单元在第e个裂缝端点处产生的法向不连续位移，mm。
[0098]
需要说明的是，现有技术中，计算应力强度因子的模型有很多，但是本发明需要实现裂缝扩展过程和储层流体流动的耦合，因此采用上述边界元位移不连续法进行应力强度因子的计算。该方法与嵌入式离散裂缝模型有着类似的裂缝离散思路，都是通过将水力裂缝离散为一个个裂缝单元来实现求解，因此能够实现裂缝扩展过程和储层流体流动的顺利耦合。
[0099]
s3：根据t0时间下第e个裂缝端点处的应力强度因子，判断t0时间下第e个裂缝端点处是否发生裂缝扩展：
[0100]
若第e个裂缝端点处，裂缝不发生扩展，则重复步骤s2-s3，计算t0时间下第e 1个裂缝端点处的应力强度因子，并判断t0时间下第e 1个裂缝端点处是否发生裂缝扩展；
[0101]
若第e个裂缝端点处，裂缝发生扩展，则第e个裂缝端点所在的裂缝沿第e个裂缝端点方向增加增加一个裂缝单元，且裂缝单元总数n
l
＝n
l
1，新增裂缝单元的流体压力与其相邻裂缝单元内的流体压力相同，新增裂缝单元的开启时间重复步骤s2-s3，以新增裂缝单元后的数据为基础，计算t0时间下第e 1个裂缝端点处的应力强度因子，并判断t0时间下第e 1个裂缝端点处是否发生裂缝扩展；
[0102]
当所有裂缝端点处的裂缝扩展情况判断结束后，即可获得t0时间下裂缝扩展之后的裂缝单元总数裂缝单元编号t0时间下第l个裂缝单元的开启时间t0时间下第l个裂缝单元内的流体压力
[0103]
在一个具体的实施例中，判断t0时间下第e个裂缝端点处是否发生裂缝扩展的具体判断方法如下：将t0时间下的第e个裂缝端点处的应力强度因子与储层岩石的断裂韧性k
ic
进行对比：若则在第e个裂缝端点处，裂缝不发生扩展；若则裂缝会在第e个裂缝端点处发生扩展。
[0104]
s4：以步骤s3获得的数据为基础，结合嵌入式离散裂缝，运用气-水两相双重介质渗流模型进行数值计算，获得压裂施工过程中t1时间下的水力裂缝单元内的流体压力、基质系统的压力分布、基质系统的含水饱和度分布、微裂缝系统的压力分布、以及微裂缝系统的含水饱和度分布；所述气-水两相双重介质渗流模型包括：
[0105]
(1)水力裂缝单元向微裂缝滤失模型：
[0106][0107]
[0108]
式中：q
f-fw
为水力裂缝单元与微裂缝网格间的滤失量，m3；c为滤失系数；t为当前时刻，s；tf为水力裂缝单元开启时间，且kf为微裂缝网格的渗透率，d；lf为水力裂缝单元长度，m，lf＝ξ
l
且h为储层厚度，m；μw为压裂液粘度，mpa
·
s；为微裂缝网格内的点到水力裂缝单元的平均法向距离，m；pf为水力裂缝单元流体压力，mpa，且p
fw
为微裂缝网格的水相压力，mpa；af为微裂缝网格的面积，m2；l
f-f
为微裂缝网格的面积单元到裂缝k的距离，m；
[0109]
(2)水力压裂中的单相流动：
[0110][0111][0112]req
＝0.14[(lf)2 (hf)2]
1/2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0113]
式中：β为单位转换系数；kf为水力裂缝的渗透率，d；bw为压裂液的体积系数；δ
well
为判断水力裂缝单元与井筒相交关系的系数，如果裂缝单元与井筒相交，则δ
well
＝1；不相交，则δ
well
＝0；q
fw
为水力裂缝单元与井筒之间的流量交换，m3；vf为水力裂缝单元的体积，m3；φf为水力裂缝的孔隙度，％；wf为水力裂缝单元宽度，m；p
wf
为井底流压，mpa；r
eq
为有效半径，m；r
well
为井筒半径，m；s为表皮系数；hf为水力裂缝单元高度，m；
[0114]
(3)压裂过程中的基质系统和微裂缝系统的渗流模型：
[0115][0116][0117][0118][0119]
p
mc
＝p
mg-p
mw
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0120]
p
fc
＝p
fg-p
fw
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0121]
式中：为哈密顿算子；k
frw
、k
frg
为微裂缝网格中液相和气相的相对渗透率；δf为微裂缝网格是否含有水力裂缝的判断参数，当微裂缝网格有水力裂缝穿过时，δf＝1；当微裂缝网格无水力裂缝穿过时，δf＝0；vf为微裂缝网格的体积，m3；km为基质渗透率，md；k
mrw
、k
mrg
为微裂缝网格中液相和气相的相对渗透率；p
mw
、p
mg
为基质网格中液相和气相的压力，mpa；φf、φm为微裂缝和基质的孔隙度，％；s
fw
、s
fg
为微裂缝网格中的液相和气相饱和度；μg为气体的粘度，mpa
·
s；bg为气体的体积系数；p
fg
为微裂缝网格的气相压力，mpa；s
mw
、s
mg
为基质网格中的液相和气相饱和度；p
fc
、p
mc
为微裂缝和基质的毛管力，mpa；
[0122]
(4)初始条件：
[0123][0124][0125][0126]
式中：pf(l)为渗流模型计算中第l个裂缝单元的流体压力，mpa；p
fw(i,j)
、p
mw(i,j)
为渗流模型计算中(i,j)网格位置处的微裂缝和基质系统液相压力，mpa；p
fwi,j,t0
、p
mwi,j,t0
为t0时间下(i,j)网格位置处的微裂缝和基质系统液相压力，mpa；s
fw(i,j)
、s
mw(i,j)
为渗流模型计算中(i,j)网格位置处的微裂缝和基质系统含水饱和度；s
fwi,j,t0
、s
mwi,j,t0
为t0时间下(i,j)网格位置处的微裂缝和基质系统含水饱和度；t
f(l)
为渗流模型计算中第l个裂缝单元的开启时间，s；t为当前时间，s。
[0127]
s5：以步骤s4获得的数据为基础，重复步骤s2-s4，计算t2时间下的水力裂缝单元内的流体压力、基质系统的压力分布、基质系统的含水饱和度分布、微裂缝系统的压力分布、以及微裂缝系统的含水饱和度分布；
[0128]
s6：重复步骤s5，直至时间达到压裂施工时间t
end
，获得t
end
时间下的水力裂缝单元数、基质系统的含水饱和度分布、微裂缝系统的含水饱和度分布；
[0129]
s7：根据步骤s6获得的t
end
时间下的数据，计算压裂改造体积；所述压裂改造体积通过下式进行计算：
[0130][0131][0132]
式中：v0为压裂改造体积，m3；α
1,i,j
为基质系统中压裂改造范围的判断参数，当为基质系统中压裂改造范围的判断参数，当时，(i,j)网格位置处的基质系统受到改造，α
1,i,j
＝1；当时，α
1,i,j
＝0；h为储层厚度，m；φm为基质的孔隙度，％；α
2,i,j
为微裂缝系统中压裂改造范围的判断参数，当时，(i,j)网格位置处的微裂缝系统受到改造，α
2,i,j
＝1；若时，α
2,i,j
＝0；φf为微裂缝的孔隙度，％；为压裂结束后的水力裂缝长度，m；为t
end
时间下水力裂缝单元数。
[0133]
在一个具体的实施例中，以川南地区某区块的页岩气井为例，采用本发明所述通过耦合储层流动确定水平井压裂改造体积的方法对其进行模拟计算，获得水力压裂施工后的压裂改造体积。具体的计算步骤如下：
[0134]
(1)依据储层情况在直角坐标系下建立结构化网格，每个网格的长度和宽度均设置为1m，整个储层被划分为180
×
350的结构化网格，即ni＝180，nj＝350。通过公式(22)～(23)用原始地层压力p0＝40mpa和原始地层含水饱和度s
w,0
＝0.3给每一个储层网格处的微
裂缝系统和基质系统赋值，即s
fwi,j,t0
＝s
mwi,j,t0
＝0.3。
[0135]
(2)在已建立的网格的基础上，添加6条簇间距为20m的初始水力裂缝，每条初始裂缝的长度为3m，被储层网格划分成3个裂缝单元。此时，水力裂缝端点总数为u＝12，水力裂缝单元的总数n
l
＝18，每个水力裂缝单元长度ξ
l
＝1m，每个裂缝单元内的流体压力p
f,1
＝p
f,2
＝
…
＝p
f,18
＝40mpa，每个裂缝单元的开启时间t
f,1
＝t
f,2
＝
…
＝t
f,18
＝0s。
[0136]
(3)基于已建立的储层网格和裂缝单元，通过公式(24)获得12个裂缝端点处的流体压力p
1,0
、p
2,0
、
…
p
12,0
。
[0137]
(4)通过公式(1)～(5)计算t＝0s时第1个裂缝端点处的应力强度因子k
i,1,0
，将其与页岩储层的断裂韧性进行对比，若k
i,1,0
＞k
ic
，则裂缝会在第1个裂缝端点处发生扩展，水力裂缝总网格数n
l
＝n
l
1，新增裂缝单元的流体压力用与其相邻的裂缝单元内的流体压力进行赋值，新增裂缝单元的开启时间为0s；若k
i,1,0
≤k
ic
，则裂缝在第1个裂缝端点处不扩展，裂缝长度不变。
[0138]
(5)继续计算第2～12个裂缝端点处的应力强度因子k
i,2,0
、k
i,3,0
…ki,12,0
，并将计算结果与k
ic
进行比较，得出每一个裂缝端点处的裂缝扩展情况，最终获得，t＝0s时裂缝扩展延伸之后的裂缝单元总数n
l,0
；每个裂缝单元的开启时间t
f,1,0
、每个裂缝单元内的流体压力p
f,1,0
、
[0139]
(6)通过公式(17)～(19)将步骤(1)-(5)中计算得到的参数转换为数值模拟模型的初始条件。
[0140]
(7)通过对公式(6)～(16)计算获得t＝1s时每个水力裂缝单元内的流体压力p
f,1,1
、每个基质网格的压力分布p
mwi,j,1
；每个基质网格的含水饱和度分布s
mwi,j,1
、每个微裂缝网格的压力分布p
fwi,j,1
、每个微裂缝网格的含水饱和度分布s
fwi,j,1
。在本实施例中，滤失系数c＝0.5，水力裂缝单元宽度wf＝0.004m。
[0141]
(8)以步骤(7)获得的参数为基础，重复步骤(2)～(7)，计算t＝2s时间下的储层基质、微裂缝和水力裂缝参数；如此循环，计算直到压裂施工完成时t＝t
end
结束，计算结束后，获得水力裂缝单元总数为储层基质不同位置处的含水饱和度微裂缝系统不同位置处的含水饱和度基质和微裂缝的含水饱和度分布结果分别如图1-2所示。
[0142]
(9)根据步骤(8)获得的参数，结合公式(21)计算获得6条水力裂缝的总长度
[0143]
(10)根据步骤(8)获得的参数，通过与初始含水饱和度的对比获得基质系统和微裂缝系统中压裂改造范围的判断参数α
1,i,j
和α
2,i,j
。
[0144]
(11)根据步骤(8)-(10)的参数，结合公式(20)计算获得压裂改造体积v0＝18841.44m3。
[0145]
本发明考虑了压裂过程中储层流体的流动，在获取储层含水饱和度场和压裂液滤失情况的基础之上对压裂改造体积进行计算，如此获得的计算结果更加符合实际，与现有技术相比，具有显著的进步。
[0146]
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质
对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。