本发明涉及一种复杂裂缝网络渗透率计算方法，属于非常规油气开发
技术领域：
。
背景技术：
：21世纪经济快速发展，能源需求不断攀升，煤、石油、天然气等传统能源的消耗也逐步增加，能源压力日益增大。在能源需求的新形势下，作为常规油气能源的有效补充，页岩气以非常规天然气的形式进入人们的视野。近年来，页岩气凭借其巨大的资源潜力与经济价值而备受关注，并且随着水平井钻井技术和体积压裂技术的发展，日益成为油气资源勘探开发的重要领域和目标。页岩气储集层在成岩及烃类演化过程中，容易在有机质边缘或薄弱面处形成规则、不规则天然裂缝，大量发育的天然裂缝对页岩气质量传输影响显著。而在页岩气开采过程中，都会用到压裂施工，由于页岩储层本身具有较为发育的天然裂缝，并且页岩由石英等脆性矿物组成，在压裂后改造区域会形成人工裂缝与天然裂缝相互交叉的复杂裂缝网络。为了能够准确模拟出页岩气压裂施工后的产能，关键在于建立符合现场实际的页岩气复杂裂缝网络，并且准确计算出其渗透率。而如何准确获取复杂裂缝网络渗透率是一个亟待解决的难题。目前国内外有大量学者开展了关于裂缝网络建模的相关研究，对于流体在裂缝中的流动主要有连续体模型、离散裂缝网络模型以及两者的混合模型。针对页岩压裂后形成复杂裂缝问题的特殊性，通常选择离散裂缝网络模型进行研究。分形维数在分析岩石裂缝中应用较为广泛，通常情况下分形维数量化了曲线或者曲面在一定范围内填充空间的程度，可以更接近实际地来描述裂缝的分布及其形态特征。目前，在离散裂缝网络建模过程中对页岩压裂后裂缝形态的表征有所欠缺，大多数的裂缝形态主要是根据经验参数取得分形维数的大小。同样地，在复杂裂缝网络渗透率计算时，目前存在不能反映真实储层情况，不能同时考虑张、剪切裂缝以及围压等因素的缺点。本发明模型通过分形维数生成的裂缝网络更加接近实际地层条件下裂缝的情况，且考虑裂缝密度、裂缝开度、围压，以及张开裂缝和剪切裂缝同时存在的多种地层条件下的因素，通过在在裂缝复杂网络生成的基础上，利用改进的立方定律计算出渗透率也因此更加接近实际地层渗透率。技术实现要素：本发明提供了一种复杂裂缝网络渗透率计算方法，其目的在于，解决现有技术中存在的上述问题。本发明的技术方案如下：一种复杂裂缝网络渗透率计算方法，其特征在于，包括以下步骤：s1、从页岩储层段钻取岩心，并将钻取的岩心制成直径为2.5cm、长度为5.0cm的标准岩心；s2、对岩心端面进行ct扫描，扫描后图像进行灰度处理，识别出裂缝条数、位置及形态特征；s3、基于s2中的图像，运用盒子法来获取表征岩石复杂裂缝网络的分形维数；s4、基于s3中的分形维数，利用蒙特卡洛随机建模理论建立二维离散裂缝网络模型；s5、基于二维离散裂缝网络模型，通过改进的立方定律求解复杂裂缝网络渗透率。进一步地，所述步骤s2中对岩心进行ct扫描即利用形貌扫描仪拍下岩心端面的裂缝细节；利用编制的程序，将岩石裂缝面间距图进行灰度处理，并将灰度图处理成二值图，完成数字化处理。进一步地，所述步骤s3的盒子法是用正方形格子去覆盖裂缝间距图；给定盒子尺码δ，可以计算出覆盖住裂缝间距所需的总盒子数目n；运用盒子法求得的分形维数表达式为：式中：d为分形维数，无量纲；n为盒子数目，无量纲；δ为盒子尺码，无量纲；i为计数符号。进一步的，所述步骤s4的二维离散裂缝网络模型的建模方法为蒙特卡洛随机建模方法，所述二维离散裂缝网络模型的生成由以下步骤组成：a、模拟区域生成；b、裂缝数目确定；c、裂缝位置表征；d、裂缝几何参数；e、裂缝网络模型生成。进一步的，所述步骤s5中复杂裂缝网络渗透率包括正常张开型裂缝渗透率、剪切型裂缝渗透率，计算渗透率过程中，由张力形成的张开型和由剪应力形成的剪切型裂缝渗透率分别进行计算，最后将正应力和剪应力诱导产生裂缝所具有渗透率进行叠加，所述正应力和剪应力诱导产生裂缝所具有渗透率的计算公式如下：式中：kx为复杂缝网在x方向上的渗透率，10-3μm2；ky为复杂缝网在y方向上的渗透率，10-3μm2；fx为在x方向上正应力产生的裂缝等效频率，1/m；fy为在y方向上正应力产生的裂缝的等效频率，1/m；br为裂缝的初始平均开度，μm；bmax为受压下裂缝最大开度，μm；αx为x方向上的压力系数，无量纲；αy为y方向上的压力系数，无量纲；σx为x方向上正应力大小，mpa；σy为y方向上正应力大小，mpa；βx为在xy下的函数斜率，无量纲；βy为在yx下的函数斜率，无量纲；fdx为在x方向上剪应力产生的裂缝等效频率，1/m；fdy为在y方向上剪应力产生的裂缝等效频率，1/m；dmax为裂缝发生剪切破坏后的最大开度，μm；k—x与y方向应力比，无量纲；γx——x方向剪切压力系数，无量纲；γy——y方向剪切压力系数，无量纲；kc——裂缝变形的临界压力比，无量纲；kc为裂缝变形的临界压力比，无量纲。本发明的有益效果为：本发明通过在钻取地层实际岩心的基础上，通过扫描获取分形维数，并根据分形维数利用蒙特卡洛随机建模理论生成裂缝网络模拟地层复杂缝网情况；在此基础上，通过改进的立方定律定量计算复杂缝网渗透率。该方法充分考虑了裂缝密度、裂缝开度、围压，以及张开裂缝和剪切裂缝同时存在等因素，通过该方法计算出的复杂缝网渗透率也更加接近储层真实情况。附图说明为了更清楚地说明本发明实施方式的技术方案，下面将对实施方式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它相关的附图。图1本发明提供的流程图；图2本发明提供的1号岩心扫描图；图3本发明提供的2号岩心扫描图；图4本发明提供的3号岩心扫描图；图5本发明提供的模型在x方向上渗透率与岩心实验测试渗透率的对比图；图6本发明提供的模型在y方向上渗透率与岩心实验测试渗透率的对比图；图7本发明提供的模型计算得到的渗透率与岩心实验测试渗透率的对比图；图8本发明提供的模型计算得到的不同水平主应力下，x方向上渗透率的变化图；图9本发明提供的模型计算得到的不同水平主应力下，y方向上渗透率的变化图；图10本发明提供的模型计算得到的不同水平主应力下，复杂裂缝网络渗透率的变化图；图11本发明提供的模型计算得到的不同分形维数下，复杂裂缝网络渗透率的变化图；图12本发明提供的模型计算得到的不同面密度，复杂裂缝网络渗透率的变化图；图13本发明提供的模型计算得到的不同裂缝开度，复杂裂缝网络渗透率的变化图。具体实施方式下面结合附图和实施例对本发明进一步说明。为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施方式，都属于本发明保护的范围。如图1所示，本发明提供的一种复杂裂缝网络渗透率计算方法，包括以下步骤：s1、从页岩储层段钻取岩心，并将钻取的岩心制成直径为2.5cm、长度为5.0cm的标准岩心；s2、对岩心端面进行ct扫描，扫描后图像进行灰度处理，识别出裂缝条数、位置及形态特征；s3、基于s2中的图像，运用盒子法来获取表征岩石复杂裂缝网络的分形维数；s4、基于s3中的分形维数，利用蒙特卡洛随机建模理论建立二维离散裂缝网络模型；s5、基于二维离散裂缝网络模型，通过改进的立方定律求解复杂裂缝网络渗透率。在本发明提供的实施例中，所述步骤s2中对岩心进行ct扫描即利用形貌扫描仪拍下岩心端面的裂缝细节；利用编制的程序，将岩石裂缝面间距图进行灰度处理，并将灰度图处理成二值图，完成数字化处理。在本发明提供的实施例中，所述步骤s3的分形维数在分析岩石裂缝中应用较为广泛，可以更接近实际来描述裂缝的分布及其形态特征。盒子法是用正方形的格子去覆盖岩石裂缝间距图，盒子大小是变化的；给定盒子尺码δ，可以数出覆盖住岩石裂缝间距所需的总盒子数目n；假设第i步覆盖使用δi×δi的格子，所需盒子数目为ni(δi)，在第i+1步需要的格子(δi+1)×(δi+1)，则需要盒子数目ni+1(δi+1)。得到任意两个尺度下所需盒子数之比与尺码之比的关系如下：推广到：n＝aδ-d(2)得到分形维数表达式为：式中：d为分形维数，无量纲；n为盒子数目，无量纲；δ为盒子尺码，无量纲；i为计数符号。在本发明提供的实施例中，所述步骤s4的二维离散裂缝网络模型的建模方法为随机建模方法，所述二维离散裂缝网络模型的生成由以下步骤组成：a、模拟区域生成：模拟区域可以根据地层中裂缝系统分布范围或者根据研究需要来选择。裂缝网络模拟时，一般先要生成裂中心点。考虑到有裂缝只有部分在模拟区域内，中心点位置可能在模拟区域外。因此，为减少边界效应的影响，应扩大模拟区域。为了便于区别，扩大后的模拟区域称为生成域，而原来模拟区域称为分析区。为保证模拟结果可靠，生成域应该足够大，原则上每一个侧面距分析域的距离应大于裂缝可能出现的最大半径。b、裂缝数目确定：模拟区域的大小确定后，就可以根据裂缝的密度来确定该区域内每组裂缝的数目。若迹长服从指数分布，裂缝的面密度可由下式确定。c、裂缝空间位置表征：对于二维离散裂缝网络模型，裂缝的位置可以由中心点及迹长，或者裂缝的两个端点完全确定，即式中：θ为裂缝与x轴的夹角，°；l——裂缝长度，m；x，y——生成裂缝坐标随机变量，m；x0，y0——生成裂缝中心点坐标随机变量，m。d、裂缝几何参数：a、按均匀分布生成裂缝中心点坐标随机变量x0，y0，做二维模拟时只需生成其中的两个变量；b、生成裂缝的倾向及倾角，生成二维离散裂缝模型时，根据研究目的来选择其中一个变量；c、生成裂缝的迹长；d、生成裂缝宽度，通常把每组裂缝的宽度设为一常量。e、裂缝网络模型生成及简化。当各组裂缝的数目、中心点坐标、裂缝产状、裂缝半径及缝宽确定后，生成域内裂缝也就完全确定，生成域裂缝生成后，根据分析几何边界来生成裂缝网络模型。在本发明提供的实施例中，所述步骤s5中复杂裂缝网络渗透率包括正常张开型裂缝渗透率、剪切型裂缝渗透率，计算渗透率过程中，正应力和剪应力分别进行计算，最后将正应力和剪应力诱导产生裂缝所具有渗透率进行叠加，表达式如下：kx＝knx+kdx(5)ky＝kny+kdy(6)式中：kx——在x方向上复杂缝网的渗透率，10-3μm2；ky——在y方向上复杂缝网的渗透率，10-3μm2；knx——在x方向上由于裂缝正常闭合后的渗透率，10-3μm2；kny——在y方向上由于裂缝正常闭合后的渗透率，10-3μm2；kdx——在x方向上由于裂缝发生剪切破坏产生的渗透率，10-3μm2；kdy——在y方向上由于裂缝发生剪切破坏产生的渗透率，10-3μm2。在离散裂缝网络中，由于单条裂缝并不能反映出不同长度不同开度的裂缝，因此，在发明中使用等效裂缝频率(fx，fy)和等效开度(bx,by)的概念，使用这种形式来相对接近的表述离散裂缝网络中的复杂裂缝变形。通过等效频率和等效开度，正应力产生的张开或者闭合裂缝在x和y方向上(knx,kny)的渗透率可以通过下式计算所得：等效开度是指用一个值来近似代表整个离散裂缝网络所有裂缝的宽度，然后根据式(7、8)的立方定律来表达离散裂缝网络的水力特性。其中裂缝的等效频率可以通过数值实验和初始开度反推得到：式中：bi为裂缝的初始平均开度，μm；裂缝的开度与地应力之间的关系，可以由如下公式表示：b＝br+bm＝br+bmaxexp(-α'σ)(11)式中：br为裂缝的初始平均开度，μm；bmax为受压下裂缝最大开度，μm；α'为正常闭合下正常剪切下的函数斜率，无量纲；σ为正应力大小，mpa。根据式(11)可以从一维的公式推导到二维条件，在x方向上和在y方向上裂缝开度的大小：bx＝br+bmax＝br+bmaxexp{-(αxσx+αyσy)}(12)式中：αx为x方向上应力系数，无量纲；αy为y方向上应力系数，无量纲；σx为x方向上正应力大小，mpa；σy为y方向上正应力大小，mpa。同理：by＝br+bmax＝br+bmaxexp{-(βxσx+βyσy)}(13)式中：βx为正常闭合下正常剪切下的函数斜率，无量纲；βy为正常闭合下正常剪切下的函数斜率，无量纲；σx为x方向上正应力大小，mpa；σy为y方向上正应力大小，mpa。接下来计算剪切应力条件下裂缝渗透率，计算公式如下：式中：fdx为在x方向上剪切裂缝的等效频率，1/m；fdy为在y方向上剪切裂缝的等效频率，1/m；dx为在x方向上剪切裂缝的等效开度，μm；dy为在y方向上剪切裂缝的等效开度，μm。由于只有在临界压力附近的裂缝才会剪切破坏，因此只有部分裂缝会出现剪切变形，其它大部分裂缝都不会出现剪切变形。所以发生剪切破坏的裂缝产生的渗透率会小于之前计算所得的正常闭合裂缝产生的渗透率。由裂缝剪切变形产生的开度大小可以通过下式计算：对于k<kc：dx＝0(16)dy＝0(17)对于k≥kc：dx＝dmax[1-exp{-γx(k-kc)}](18)dy＝dmax[1-exp{-γy(k-kc)}](19)式中：kc为裂缝变形的临界压力比，无量纲；dmax为裂缝发生剪切破坏后的最大开度，μm。导致裂缝破坏的临界压力比和临界方向值可以通过库伦破坏准则来计算：式中：φ为岩石内摩擦角，°；φf为岩石临界破坏角，°。最后，x和y方向上的渗透率(kx,ky)大小即可通过叠加正常闭合和剪切变形的裂缝渗透率得到，如下：式中：kx为裂缝在x方向上的渗透率，10-3μm2；ky为裂缝在y方向上的渗透率，10-3μm2；fx为在x方向上正应力产生的裂缝的等效频率，1/m；fy为在y方向上正应力产生的裂缝的等效频率，1/m；br为裂缝的初始平均开度，μm；bmax为受压下裂缝最大开度，μm；αx为x方向上的压力系数，无量纲；αy为y方向上的压力系数，无量纲；σx为x方向上正应力大小，mpa；σy为y方向上正应力大小，mpa；βx为在xy下的函数斜率，无量纲；βy为在yx下的函数斜率，无量纲；fdx为在x方向上剪应力产生的裂缝的等效频率，1/m；fdy为在y方向上剪应力产生的裂缝等效频率，1/m；dmax为裂缝发生剪切破坏后的最大开度，μm；kc为裂缝变形的临界压力比，无量纲。实施例：如图2～4所示，选取三块岩心，首先对岩心进行ct扫描，对扫描图像进行灰度处理，识别出裂缝在岩石表面的分布位置及形态特征，并根据本发明模型获得的裂缝长度、裂缝密度和分形维数等参数建立二维离散裂缝网络。根据扫描结果获取的分形维数等参数见表1。表1计算基础参数ⅰ为了进一步说明本发明方法的适用性，进一步使用表2的数据进行了敏感性分析。表2计算基础参数ⅱ参数单位数值分形维数-1.41最大裂缝长度rmaxm150最小裂缝长度rminm1裂缝区域宽度m100裂缝区域长度m100裂缝面密度条/m25bx上的压力系数(αx，αy)-0.035，0.132by上的压力系数(βx，βy)-0.113，0.069dx,dy的剪切应力系数γx，γy-1.88，1.95剩余开度(br)*μm50最大机械变形开度(bmax)*μm25最大剪切开度(dmax)*μm50裂缝初始开度(bi)*μm20摩擦角°24.9临界压力比-2.66水平主应力σxmpa1～40水平主应力σympa1～40注：*实验测试数据如图5所示，根据式(22)求得x方向上渗透率大小，得到模型在x方向上渗透率与岩心实验测试渗透率的对比。如图6所示，根据式(23)求得y方向上渗透率大小，得到模型在y方向上渗透率与岩心实验测试渗透率的对比。如图7所示，根据式(24)计算得到的渗透率与岩心实验测试渗透率的对比。如图8所示，根据式(22)求得x方向上渗透率大小，得到模型在x方向上渗透率随x与y方向应力比值增加的变化情况。如图9所示，根据式(23)求得y方向上渗透率大小，得到模型在y方向上渗透率随x与y方向应力比值增加的变化情况。如图10所示，根据式(24)求得复杂裂缝网络渗透率大小，得到该区域的平均渗透率随压力的变化情况。如图11所示，复杂裂缝网络渗透率计算模型，得到不同分形维数下，复杂裂缝网络渗透率的变化情况；如图12所示，复杂裂缝网络渗透率计算模型，得到不同面密度下，复杂裂缝网络渗透率的变化情况；如图13所示，复杂裂缝网络渗透率计算模型，得到不同裂缝开度下，复杂裂缝网络渗透率的变化情况；首先通过对比x与y方向渗透率与实验测试所得岩心渗透率的值进行对比可以看出，岩心1计算得到的渗透率为30.5×10-3μm2，实验测得的渗透率为31.45×10-3μm2。对于岩心2和岩心3，计算得出的渗透率大小分别为55.7×10-3μm2和11.2×0-3μm2，实验测得的渗透率分别为50.06×10-3μm2和13.13×10-3μm2。通过对比可以发现，本发明方法计算所得渗透率与实验室测试得到的渗透率误差较小。从图10可以看出，在水平主应力不等比例增加的条件下，保持x方向上的水平主应力不变，然后增大y方向上水平主应力，应力的比值k从0.5增加到5。岩石的渗透率在k值增大的情况下，出现了先减小后增大的趋势，裂缝渗透率先从2.3×10-3μm2减小到了0.7×10-3μm2。当k值大于2.5以后，渗透率又增加到了2.1×10-3μm2；从图9可以看出，在y轴上渗透率的变化没有x轴方向明显，原因在于x轴上的应力在不断增大，y方向上的渗透率增加受限。通过对比不同分形维数可以发现，离散裂缝网络的分形维数从1.15增加到1.45的过程中，离散裂缝网络的渗透率在不断的增加(图11)；通过对比面密度可以发现随着区域内裂缝密度的增大，渗透率的增加非常明显，说明在区域内裂缝密度的增加对渗透率的贡献非常大，能够显著的改变区域内渗透率的大小。当分形维数从1.15增加到1.45时，在裂缝面密度等于5条/m2时，缝网渗透率从0.15×10-3μm2增加到0.92×10-3μm2；在裂缝面密度等于10条/m2时，缝网渗透率从0.97×10-3μm2增加到1.73×10-3μm2；在裂缝面密度等于15条/m2时，缝网渗透率从1.75×10-3μm2增加到2.71×10-3μm2(图12)；通过对比裂缝开度可以发现，当裂缝开度为30μm时,随着裂缝密度从2条/m2增加到16条/m2，复杂缝网渗透率从0.1×10-3μm2增加到了0.81×10-3μm2；当裂缝开度为120μm时,随着裂缝密度从2条/m2增加到16条/m2，渗透率随着面密度从2×10-3μm2增加到了16×10-3μm2。本发明提供的一种复杂裂缝网络渗透率测试方法的有益效果为：本发明模型利用对实际岩心进行扫描，通过分形维数生成的裂缝网络，并且考虑了裂缝密度、裂缝开度、围压多种地层实际因素，使得本发明模型更加接近实际地层条件下裂缝的形态分布，通过该模型计算出渗透率也因此更加接近实际地层渗透率。以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。当前第1页12