李滔,李骞,胡勇,彭先,冯曦,朱占美,赵梓寒
(1.中国石油西南油气田公司勘探开发研究院,成都 610041;2.清华大学工程力学系,北京 100084;3.中国石油西南油气田公司,成都 610051)
岩石中存在性质不同的微裂缝,且相互连通形成复杂的微裂缝网络[1-6]。微裂缝网络不同程度地提高了岩石渗流能力,使得相同孔隙度下的岩石渗透性差异大,储集层渗透性与孔隙度相关性较差,加大了储集层储量可动性和渗流能力研究的难度[7-8],因此需要进一步提出评价微裂缝网络整体性质的新方法。
岩石微裂缝主要通过岩石薄片、铸体薄片和扫描电镜等技术识别[2]。前人用裂缝长度[9]、等效裂缝开度[10]、裂缝表面粗糙度[11-12]、裂缝曲折率[12-13]和倾角[9]等参数表征单一微裂缝的基本特征,但镜下观测结果难以形成对微裂缝网络的整体认识。随着微CT技术和数字岩心技术的发展,三维数字岩心可以直观展示岩石内部复杂的微裂缝网络形态,定性认识微裂缝网络性质。然而定量评价微裂缝网络性质难度大。由于不同微裂缝的倾角往往不同,难以使用统一的倾角来描述微裂缝网络整体走向[4-6]。裂缝网络连通性常被定义为每条裂缝上的平均交点数[9,14],或者为连通部分裂缝网络的密度占总体裂缝网络密度的比例[15],但一个局部连通的密集裂缝网络,其整体渗透率不一定大[16]。所以目前缺乏定量表征微裂缝网络性质的有效手段[2,9,16-17]。
基于有限差分的数值模拟方法[18]适用于宏观尺度渗流模拟,但难以满足微裂缝网络渗流模拟需求。微观数值模拟方法主要包括孔隙网络模拟法[19]、计算流体动力学方法(CFD)[20]和格子-玻尔兹曼方法(LBM)[17]。其中,基于微观动力学的LBM具有清晰的微观粒子背景,便于处理不规则几何边界,适用于研究不规则微裂缝与孔隙的耦合渗流[17,21-23]。Gao等[24]和李滔等[25]采用LBM开展多孔介质渗流模拟,证明了其可靠性。多松弛格子-玻尔兹曼方法(MRT-LBM)克服了单松弛格子-玻尔兹曼方法(SRT-LBM)模拟精度受弛豫时间影响的不足,有效提高了模拟结果的准确性[26]。
本文结合对岩石微裂缝和孔隙结构的认识,提出定量评价微裂缝网络连通性和走向的方法,采用数值算法构建一系列裂缝性多孔介质和裂缝-孔隙性多孔介质模型,运用 MRT-LBM模拟模型中流体渗流,以揭示微裂缝网络对岩石渗透率的影响机理。
四川盆地致密砂岩和碳酸盐岩的岩心薄片、扫描电镜和微CT扫描图像如图1所示。其中,图1a—图1b为致密砂岩薄片鉴定图像,图 1c—图 1d为碳酸盐岩薄片鉴定图像,图1e—图1h为致密砂岩扫描电镜图像,图1i—图1j为分辨率为0.62 μm的致密砂岩微CT扫描图像,图1k—图1l为分辨率为14.10 μm的碳酸盐岩微CT扫描图像。致密砂岩岩样孔隙度4.1%~7.3%,渗透率 0.05×10−3~0.14×10−3μm2;碳酸盐岩岩样孔隙度3.2%~7.5%,渗透率 0.03×10−3~1.50×10−3μm2。如图所示,致密砂岩和碳酸盐岩中均发育形态复杂的微裂缝,包括构造缝、粒间缝、层间缝、溶蚀缝等;部分微裂缝相交,形成微裂缝网络。致密砂岩微裂缝开度为0.002~0.050 mm;碳酸盐岩微裂缝开度为 0.002~0.100 mm,与王璐等[27]的认识一致。
图1 致密砂岩和碳酸盐岩中的微裂缝
微观图像和三维数字岩心直观展示了岩石中的微裂缝网络形态,可获取微裂缝开度和数量等信息。但受限于观测尺度,部分微裂缝延伸至薄片范围外,难以准确获取长度[2]。而裂缝长度影响微裂缝网络的连通性,对岩石渗流起着重要作用;同时受分辨率限制,三维数字岩心未能捕获岩石中大部分基质孔隙。若直接基于数字岩心开展模拟,会高估岩石微裂缝网络连通性、低估基质孔隙渗流能力。
高压压汞实验可准确测量岩石微观孔隙性质,以四川盆地某气藏储集层岩石为例,测得其基质孔隙半径主要为0.1~0.7 μm,近似满足正态分布(见图2)。本文结合该气藏储集层岩石数字岩心和高压压汞实验结果,采用数值算法构建裂缝性多孔介质和裂缝-孔隙性多孔介质模型。其中,裂缝性多孔介质模型(简称裂缝性模型)含有贯穿性微裂缝网络(忽略基质孔隙影响),裂缝-孔隙性多孔介质模型(简称裂缝-孔隙性模型)含有孔隙和未贯穿性微裂缝网络。裂缝-孔隙性模型构建过程如下:①设定多孔介质区域,将整个区域均设为固相;②结合高压压汞实验结果,设定多孔介质基质孔隙孔径满足图 2所示正态分布,确定不同尺寸孔隙所占节点数量,在多孔介质内随机指定孔隙的半径;当生成孔隙的节点数量达到预期,停止孔隙生长,确保孔隙连通(见图3a);③基于微观图像和数字岩心分析结果,设置微裂缝开度和长度,并将其随机分布于多孔介质;④设置不同裂缝数量,重复步骤①—③,即可生成含不同微裂缝网络的裂缝-孔隙性模型(见图3b)。只需执行步骤③和④即可生成不同裂缝性模型。
图2 四川盆地某气藏储集层岩石基质孔隙半径分布
图3 多孔介质模型
岩石渗流能力只与相互连通的孔隙和微裂缝有关,采用四连通区域标记算法[23]进一步去除多孔介质中的孤立孔隙和裂缝。基于形态学原理[28],按一定线条间隔沿多孔介质x、y方向分别插入多条直线,单位长度直线所遇到平均裂缝数的比值定义为裂缝网络走向因子(A):
图4为6个由10条长400格子,开度20格子的微裂缝组成的微裂缝网络示意图。将线条间隔设置为5格子,沿多孔介质x和y方向插入直线得到nx和ny。运用(1)式计算得到6个微裂缝网络的A值为0.05~17.61,不同走向的微裂缝网络A取值不同。当倾向于y方向的微裂缝数量多于倾向于x方向时,A值小于1;当倾向于x方向的微裂缝数量多于倾向于y方向时,A值大于1;当倾向于x方向的微裂缝数量等于y方向时,A值等于 1。当微裂缝走向均为y方向时,A值最小(A=0.05);当微裂缝走向均为x方向时,A值最大(A=17.61)。随着倾向于x方向微裂缝数量的增加,A值增大,因此A值可以定量表征微裂缝网络整体走向。
图4 不同走向的微裂缝网络示意图
微裂缝网络在多孔介质中形成新的渗流通道,微裂缝间可能直接连通,也可能通过基质孔隙间接连通。以往裂缝网络连通性的评价方法未能兼顾裂缝网络局部与多孔介质整体的关系,适用性较差[15]。本文将某一流动路径上与基质孔隙相连通的微裂缝水平长度之和的最大值(Lf,max)与多孔介质长度(L)的比值定义为微裂缝网络连通系数(f),表征流体渗流路径上微裂缝的最大占比,f取值0~1,f值越大,多孔介质中微裂缝网络整体连通性越好。当f<1时,微裂缝网络处于多孔介质内部,为裂缝-孔隙性模型;当f=1时,微裂缝网络连通多孔介质出入口端,为裂缝性模型。运用数值算法筛选多孔介质中与孔隙连通的裂缝,并判断裂缝间连通情况,便可求得Lf,max(见图 5,红色圆圈表示选择的起始点和终止点)。通过(2)式计算得到图5a和图5b中f分别为0.58和0.70。
图5 微裂缝网络连通系数计算示意图
MRT-LBM由Boltzmann方程直接推导,可表示为[23]:
M为分布函数fi的变换矩阵[29],(3)式可写为:
二维9个离散速度方向(D2Q9)模型中平衡态分布函数为:
离散速度ci定义为:c0=(0,0),c1=-c3=(c,0),c2=-c4=(0,c),c5=-c7=(c,c),c6=-c8=(-c,c)。对角矩阵S定义为:
本文选取sρ=sj=1.0,se=1.19,sε=1.4,sq=8(2-sv)/(8-sv)[29],sv与运动黏度有关。宏观流体的密度、速度和运动黏度如(7)式所示。
多孔介质孔隙表面不规则,可直接应用经典的LBM边界条件[28]。本文中多孔介质出入口端均采用压力边界[25,30],上下边界采用标准反弹边界。
采用Fortran语言编写MRT-LBM模拟程序,并验证其模拟多孔介质流体渗流的准确性。
2.2.1 裂缝性模型流体渗流模拟
结合储集层岩石微裂缝的开度、数量、长度等数据,采用数值算法构建 28个裂缝性模型,通过MRT-LBM模拟裂缝开度分别为30,60,120 μm时多孔介质中的流体(水)连续流动。裂缝性模型构建参数及模拟条件如表 1所示,模拟中流体流动均为层流(Re值较小)。模拟偏差用(8)式计算,模拟偏差小于1×10−6时,模拟结果不再变化,结束迭代过程[30-31]。
表1 裂缝性模型构建参数及MRT-LBM模拟参数
图6所示为微裂缝开度为60 μm的裂缝性模型的模拟结果,可见裂缝网络显著影响多孔介质中流体渗流路径以及速度场分布。基于模拟结果,进一步删除模型中未参与流动的微裂缝,得到裂缝性模型的连通孔隙度(即多孔介质中相互连通孔隙的体积与多孔介质总体积之比)。将线条间隔设置为5个格子,采用(1)式计算裂缝性模型的裂缝网络走向因子。裂缝性模型的裂缝网络走向因子为 0.99~2.70,连通孔隙度为2.3%~5.0%,二者无明显相关性(见图7)。
图6 开度60 μm的裂缝性模型模拟结果
图7 裂缝性模型连通孔隙度与裂缝网络走向因子的关系
2.2.2 裂缝-孔隙性模型流体渗流模拟
结合对储集层岩石性质的认识,设定裂缝-孔隙性模型的基质孔隙度为 5%,且孔隙半径满足正态分布(见图2),微裂缝孔隙度为5%。考虑到岩石中微裂缝尺度普遍大于基质孔隙,将微裂缝开度分别设为2 μm和4 μm。微裂缝分布于裂缝-孔隙性模型内部,随机设定其走向和长度,微裂缝起止点均与基质孔隙连通。模型详细构建参数见表2。为考虑随机性影响,不同开度下均生成20个裂缝-孔隙性模型,共计40个。
表2 裂缝-孔隙性模型构建参数及MRT-LBM模拟参数
结合(2)式计算裂缝-孔隙性模型的裂缝网络连通系数,采用D2Q9 MRT-LBM模拟其中流体(水)连续流动,模拟条件如表2所示。图8对比了基质模型和裂缝网络连通系数分别为 0.36,0.70,0.80时裂缝-孔隙性模型(裂缝开度2 μm)的模拟结果。可以看出,裂缝网络显著影响裂缝-孔隙性模型中流体渗流路径以及速度场分布,部分渗流通道中的流体渗流速度显著增加(见图 8c、图 8d),优势通道效应显著。随后,删除裂缝-孔隙性模型中未参与流动的微裂缝,将线条间隔设置为5个格子,采用(1)式计算裂缝-孔隙性模型的裂缝网络走向因子,得到模型的连通孔隙度与裂缝网络走向因子和裂缝网络连通系数的关系(见图9、图10)。当裂缝开度为2 μm时,裂缝-孔隙性模型的连通孔隙度为8.4%~9.5%,裂缝网络走向因子为0.99~1.58,裂缝网络连通系数为0.11~0.80;裂缝开度为4 μm时,裂缝-孔隙性模型的连通孔隙度为8.6%~9.2%,裂缝网络走向因子为0.89~3.63,裂缝网络连通系数为0.34~0.70。裂缝-孔隙性模型的连通孔隙度与裂缝网络走向因子和裂缝网络连通系数均无明显相关性。
图8 基质模型和裂缝-孔隙性模型模拟结果
图9 裂缝-孔隙性模型连通孔隙度与裂缝网络走向因子的关系
图10 裂缝-孔隙性模型连通孔隙度与裂缝网络连通系数的关系
通常用渗透率表征多孔介质的渗流能力,而迂曲度反映了多孔介质中流体渗流路径的曲折程度,是影响多孔介质渗透性的重要参数[25]。基于 MRT-LBM 模拟结果,采用(9)式和(10)式分别计算裂缝性模型和裂缝-孔隙性模型的渗透率和迂曲度[28],发现两个模型的渗透率与连通孔隙度均无明显相关性(见图 11、图 12)。裂缝-孔隙性模型中的基质具有相同的孔隙度和孔径分布,说明其渗透率的变化主要受微裂缝网络的影响。因此需要进一步分析多孔介质渗透率与裂缝网络性质的关系,探究微裂缝网络影响多孔介质渗透性的机理。
图11 裂缝性模型渗透率与连通孔隙度的关系
图12 裂缝-孔隙性模型渗透率与连通孔隙度的关系
裂缝性模型渗透率与裂缝网络走向因子呈较好的正相关性(见图13),当裂缝网络走向因子由1.0增加至2.7时,不同裂缝开度下的裂缝性模型渗透率增加幅度均超过160%;而其迂曲度与裂缝网络走向因子呈负相关性(见图14),二者的多项式拟合效果最好,拟合相关系数平方达 0.953(见表 3)。裂缝-孔隙性模型的渗透率和迂曲度与裂缝网络走向因子均无明显相关性(见图 15、图 16)。进一步分析裂缝性模型出口端速度剖面,发现微裂缝网络走向因子越大,多孔介质出口端流体渗流速度越大(见图17)。从微观角度分析,裂缝网络走向因子越大,裂缝性模型的迂曲度越小,渗流路径越平直,渗透率就越大。但对于裂缝-孔隙性模型(f<1),基质孔隙与微裂缝尺寸普遍存在数量级差异,基质孔隙是制约裂缝-孔隙性模型渗流能力的关键因素[32]。因此,微裂缝网络走向对裂缝性模型渗流能力的影响显著大于裂缝-孔隙性模型。
图13 裂缝性模型渗透率与裂缝网络走向因子的关系
图14 裂缝性模型迂曲度与裂缝网络走向因子的关系
表3 裂缝性模型迂曲度-裂缝网络走向因子拟合函数及裂缝-孔隙性模型渗透率-裂缝网络连通系数拟合函数
图15 裂缝-孔隙性模型渗透率与裂缝网络走向因子的关系
图16 裂缝-孔隙性模型迂曲度与裂缝网络走向因子的关系
图17 裂缝性模型速度剖面与裂缝网络走向因子的关系
裂缝-孔隙性模型的渗透率与裂缝网络连通系数呈较好的正相关关系(见图18)。当裂缝开度为2 μm且裂缝网络连通系数由0.11增至0.80时,裂缝-孔隙性模型的渗透率增加幅度达61.83%;当裂缝开度为4 μm且裂缝网络连通系数由 0.34增至 0.70时,裂缝-孔隙性模型的渗透率增加幅度达48.43%。裂缝-孔隙性模型渗透率与裂缝网络连通系数的指数式拟合效果最好(以开度2 μm为例),拟合相关系数平方为0.832(见表 3)。
图18 裂缝-孔隙性模型渗透率与裂缝网络连通系数的关系
随后分析了模型(f=0为基质模型,0 此外,部分裂缝-孔隙性模型中优势通道效应显著,通道二在裂缝网络连通系数为0.70和0.80时,以及通道四在裂缝网络连通系数为0.36时,其流量占比均超过0.5(见图19)。总体来看裂缝网络连通系数越大,多孔介质中越易出现优势通道效应(f≥0.7时尤为显著)。优势通道效应存在时,多孔介质中其余通道的流量占比减小,而通过这些通道的流体流量大多增加,少数情况下也存在个别通道传输流体流量减少的现象。 图19 模型出口端各通道流量占比与裂缝网络连通系数的关系 裂缝开度为60 μm和120 μm时的裂缝性模型渗透率分别约为裂缝开度为30 μm时裂缝性模型的4倍和16倍,其连通孔隙度、渗透率、迂曲度与裂缝网络走向因子的相关性均与裂缝开度为30 μm的裂缝性模型一致(见图 13、图 14)。即对于裂缝性模型,其渗流能力主要取决于裂缝网络走向和裂缝开度,而裂缝开度的影响更显著。裂缝-孔隙性模型的基质孔隙特征长度为0.344 μm[31],当裂缝开度与基质特征长度的比值由5.814增加至11.628时(微裂缝开度由2 μm增加至4 μm),裂缝-孔隙性模型的渗透率变化幅度小于 10%(见图 18)。所以,裂缝开度对裂缝-孔隙性模型渗流能力的影响远小于裂缝网络连通性。 目前气藏储量可动性评价方法主要包括实验法和数值模拟法,综合考虑经济因素、实验结果和现场测试数据求取储集层可动储量孔隙度下限的方法取得了较好效果。但由前文可知,含微裂缝网络多孔介质的渗流能力与孔隙度无明显相关性,因此采用孔隙度划分储集层储量可动下限的方法未考虑微裂缝网络对岩石渗流能力的影响。以本文得到的多孔介质渗透率与微裂缝网络连通性指数关系式为例,基质孔隙度和裂缝开度相同时,当裂缝网络连通系数由0增至1,多孔介质渗透性增加120.98%,储量可动性差异显著。因此针对不同气藏,需要结合微观实验分析岩石微裂缝性质(裂缝开度、数量和连通性等),评价微裂缝网络对岩石渗透性的影响程度。对于微裂缝网络影响较弱的气藏,直接使用实验法等得到的孔隙度下限评价储量动用情况;而对于微裂缝网络影响显著且水体能量较弱的气藏,根据影响程度不同幅度地下调实验法得到的储量动用孔隙度下限,提升储集层评价的准确性;对于微裂缝网络影响显著且水体能量较强的气藏,岩石中存在水驱时,优势通道效应的存在会导致波及效率下降[33],加剧水封气现象,降低储集层储量可动性。 准确获取储集层渗透性是气井产能评价的关键。试井分析得到的近井区和远井区储集层渗透性通常差异较大,非均质性显著。由本研究可知,近井区和远井区储集层渗流能力的差异可能是由不同的裂缝网络走向导致的。裂缝开度和数量相同时,当裂缝网络走向因子由1.00增加至2.64,裂缝性多孔介质的迂曲度减小幅度超过25%,渗透率增加幅度超过160%。考虑裂缝网络走向的影响可进一步提升对气井产能的认识程度。 裂缝网络走向因子和连通系数可较好地定量表征多孔介质中微裂缝网络的整体走向和连通性。随着裂缝网络走向因子的增加,裂缝性多孔介质的迂曲度显著减小,渗透率逐渐增加,且其迂曲度与裂缝网络走向因子满足多项式关系式;随着裂缝网络连通系数的增加,裂缝-孔隙性多孔介质容易出现优势通道效应,渗透率增大,且其渗透率与裂缝网络连通系数满足指数关系式。本研究结果可用于计算不同裂缝发育情况下岩石的渗透率张量,有助于提升评价储集层储量可动性和气井产能的准确性。 符号注释: A——裂缝网络走向因子,无因次;c——格子速度,m/s;ci——离散速度,m/s;cs——声速,m/s;ei——格子速度,m/s;f——裂缝网络连通系数,无因次;fi——离散分布函数,kg/m3;feq,i——平衡态分布函数,kg/m3;i——离散方向;K0——模型渗透率,10−3μm2;Ke——裂缝-孔隙性模型基质渗透率,10−3μm2;L——多孔介质长度,m;Lfn——流动路径上单条微裂缝水平方向长度,m;Lf,max——流动路径上微裂缝网络水平方向长度之和的最大值,m;mi——矩空间的函数;meq,i——矩空间的平衡态函数;M——转换矩阵;n——渗流通道上的裂缝总数;nx,ny——单位长度直线沿多孔介质x和y方向所遇到的平均裂缝数;R——相关系数,无因次;Re——雷诺数,无因次;se,sj,sq,sv,sε,sρ——能量矩、动量矩、能量通量矩、应力张量矩、能量的平方矩、密度矩,无因次;S——对角矩阵;t——时间,s;u——流体速度,m/s;u(x,y)——多孔介质(x,y)处流体速度,m/s;u(x,y,t)——t时刻多孔介质(x,y)处流体速度,m/s;ux(x,y)——多孔介质(x,y)处x方向的流体速度,m/s;x——x方向平均流体速度,m/s;wi——权系数,无因次;x——多孔介质的位置向量;x,y——直角坐标系的2个方向;δt——时间步长,s;δE——模拟结果的相对偏差,无因次;▽p——压力梯度,Pa/m;μ——动力黏度,Pa·s;ν——运动黏度,m2/s;ρ——流体密度,kg/m3;τ——迂曲度,无因次;Ω——碰撞矩阵;φ——多孔介质总孔隙度,%;φeff——多孔介质连通孔隙度,%。3.3 微裂缝开度对多孔介质渗流能力的影响
4 讨论
4.1 储量可动性和储集层评价
4.2 气井产能评价
5 结论