一种考虑多重因素影响的页岩气表观渗透率计算方法与流程

本发明涉及一种考虑多重因素影响的页岩气表观渗透率计算方法，属于非常规油气开发
技术领域：
。
背景技术：
：页岩微纳米孔道中气体流动能力是目前非常规油气领域研究热点。吸附气和自由气共存于页岩孔隙中。气体在页岩中的质量传输有以下几种机制：粘性流、knudsen扩散、表面扩散。同种气体分子与分子之间碰撞产生粘性流，分子与壁面碰撞产生knudsen扩散，吸附在孔隙壁面的气体分子沿孔隙表面蠕动产生表面扩散。由于连续性假设在微纳米孔道中不适用，一般采用气体分子运动自由程与孔隙半径的比值knudsen数判定页岩中的气体运移模式，不同knudsen数对应气体运移机制不同。目前，普遍认为无机质孔隙如黏土孔隙具有大量吸附气是建立在固－气界面吸附理论基础上，其对页岩气藏不适应。由于黏土矿物具有极强的亲水性，实际储层黏土孔隙表面吸附水膜，并且孔隙中存在一定毛细水及可动水。在黏土孔隙表面存在水膜情况下，大多数甲烷分子表现为聚集在孔隙中心，或者吸附在水膜上，而并非吸附在黏土表面，这就在很大程度上降低了黏土对甲烷的吸附能力。相比于有机质对甲烷强吸附能力，无机质孔隙吸附能力可以忽略。beskok和civan等人建立了基于knudsen数的气体表观渗透率模型，javadpour和wu等人建立了基于多重气体运移机制叠加的气体表观渗透率模型。但目前基于单管分析的页岩气表观渗透率模型忽略了页岩地层的非均质性，并且没有考虑有机质孔隙与无机质孔隙在气体传输、空间分布、孔隙半径分布上的差异。此外，目前研究页岩气体流动并没有考虑气体在微纳米孔道中的相态变化。技术实现要素：针对现有技术的不足，本发明提供了一种考虑多重因素影响的页岩气表观渗透率计算方法；本发明主要针对目前页岩地层气体渗透率预测问题，充分考虑页岩地层特征、气体性质、地层应力作用，建立了一种考虑多重因素影响的页岩气表观渗透率计算方法。术语解释：1、蒙特卡洛方法：蒙特卡洛方法是二十世纪四十年代中期提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数(通常为伪随机数)来解决很多计算问题的方法。蒙特·卡洛方法在计算物理、金融、经济学等领域应用广泛。2、孔隙相：三维数据体中所有数字0占据的空间位置。3、页岩三维数字岩心：基于页岩二维剖光扫描电镜图像采用多点地质统计学方法三维重构得到的三维数据体。本发明的技术方案为：一种考虑多重因素影响的页岩气表观渗透率计算方法，包括步骤如下：(1)采用多点地质统计学方法，基于页岩二维剖光扫描电镜图像，重建页岩三维数字岩心；(2)计算页岩三维数字岩心的孔隙度、迂曲度、多孔介质特征长度；(3)计算页岩有机质和无机质中的气体渗透率。根据本发明优选的，步骤(1)，包括步骤如下：a、将所选用的页岩二维剖光扫描电镜图像四个顶点中左上角顶点作为原点，规定x方向为页岩二维剖光扫描电镜图像短边方向，y方向为页岩二维剖光扫描电镜图像长边方向，从原点开始截取x方向的n个像素，y方向n个像素的正方形图像，作为要研究的训练图像，构建查找树，n的大小不超过页岩二维剖光扫描电镜图像短边方向的像素个数；由于在二维图像中对孔隙相每个位置进行连续的测量计算量较大，因此，选取页岩二维剖光扫描电镜图像上若干孔隙相分布密集区域的像素点作为采样点，把采样点作为初始条件数据，随机分配到临近的网格节点上，网格节点是指页岩二维剖光扫描电镜图像上每个像素点的位置；b、对于任一待模拟网格节点u，确定在给定n个条件数据值s(uα)的情况下，属性s(u)取k个状态值中任一个状态值的条件概率如式(ⅰ)所示：式(ⅰ)中，sk指属性s(u)中的任一个状态值、skα指α个条件数据值下属性s(u)取k个状态值中任一个状态值条件数据事件；dn为由n个条件数据值s(uα)联合构成的数据事件，式(ⅰ)中的分母为某个模式出现的概率,α＝1,…,n；分子为某个模式出现和待模拟网格节点u取sk的情况同时出现的概率；对于所有的未模拟的网格节点，基于条件概率分布函数，使用蒙特卡洛方法提取模拟值，统计得到的模拟值个数，若获得的模拟值个数为n×n个，则将获得的n×n个模拟值作为新生成的二维图像，若获得的模拟值个数小于n×n个，则将得到的模拟值作为下一步的条件数据，继续模拟，直到获得的模拟值个数为n×n个；c、将步骤b新生成的二维图像放置于页岩二维剖光扫描电镜图像上方,得到垂直方向上高度为2层二维图像的三维图像,将步骤b新生成的二维图像作为下一步的训练图像，重复步骤a-b，直至获得基于训练图像得到的新生成的二维图像，继续放置于已构建的三维图像的上方，构建出垂直方向上高度为3层二维图像的三维图像；d、不断重复步骤a-c，直到获得n-1层生成的二维图像，二维图像中的每一个像素都对应着三维空间的相应的体素，加上步骤a选用的的页岩二维剖光扫描电镜图像，共有n层二维图像，构成了x、y、z三个方向大小为n×n×n的多孔介质三维图像，该多孔介质三维图像成为页岩三维数字岩心。根据本发明优选的，步骤(2)，计算页岩三维数字岩心的孔隙度、迂曲度、多孔介质特征长度，是指：气体在三维多孔介质中运移时努森数kn表示为：λ为气体平均分子自由程；lc为多孔介质特征长度，基于弦长计算多孔介质特征长度，弦长的几何含义是指：用一条直线切割多孔介质，分布在孔隙相中不连续的直线；多孔介质特征长度lc通过弦长计算，如式(ⅰ)所示：式(ⅰ)中，z表示孔隙相中弦长长度，p(z)表示三维多孔介质中切割出的弦长中弦长长度为z的概率；根据弦长的物理意义，表征在岩心数据体上即为连续的0数据。通过基于页岩三维数字岩心模拟，迂曲度为：沿程路径与绝对路径之比；沿程路径为三维数字岩心中沿着孔隙相(数字0)与骨架相(数字1)的边界气体流过的长度；绝对路径指三维数字岩心边长；孔隙度为：孔隙相与页岩三维数字岩心之比。根据本发明优选的，步骤(3)，计算页岩有机质和无机质中的气体渗透率，包括：计算地层有效应力对多孔介质特征长度的影响，如式(ⅱ)所示：式(ⅱ)中，lc_stress是指考虑有效应力后的多孔介质特征长度；lc0为原始地层有效应力σeff0下的多孔介质特征长度，σeff为变化后的地层有效应力，p1为孔隙闭合压力，m为地层胶结指数；页岩中吸附气量采用langmuir等温吸附表示，计算吸附气在孔隙壁面的覆盖度θ，如式(ⅲ)所示：式(ⅲ)中，pl是指langmuir压力，单位为mpa；z是指气体压缩因子，无因次；θ是指真实气体孔隙壁面覆盖度，无因次；po是指地层压力，单位为mpa；计算考虑吸附气表面覆盖校正后的多孔介质特征长度lc_or，如式(ⅳ)所示：lc_or＝lc_stress-dmθ(ⅳ)式(ⅳ)中，dm是指气体分子直径，单位为m；计算微纳米孔道气体临界压力临界温度变化，通过式(ⅴ)、(vi)进行计算：式(ⅴ)、(vi)中，tc是指气体临界温度，单位k；pc是指气体临界压力，单位mpa；r是指理想气体常数，单位j/(mol·k)；σ是指lennard-jones作用尺寸，单位m；a是指甲烷范德华方程常数，单位mpa·kmol-2；b是指甲烷范德华方程常数，单位m6·mpa·kmol-2；ε是指lennard-jones势能参数，单位j-1；有效多孔介质特征长度lceff在无机质介质中为考虑有效应力后的多孔介质特征长度lc_stress，在有机质介质中为考虑吸附气表面覆盖校正后的多孔介质特征长度lc_or；计算气体压缩因子z，如式(ⅶ)、式(ⅷ)、式(ⅸ)所示：ppr是指拟压力，无因次；tpr是指拟温度，无因次；t是指地层温度，单位为k；计算粘度，如式(ⅹ)、式(ⅺ)、式(ⅻ)、式(xiii)、式(xiv)所示：μ＝(1×10-4)kexp(xρy)(ⅹ)y＝2.447-0.2224x(xiv)式(ⅹ)、式(ⅺ)、式(ⅻ)、式(xiii)、式(xiv)中，μ是指气体粘度，单位为pa·s,x、y、k为计算气体粘度公式中的相关参数，受地层温度和气体摩尔质量影响。m是指气体摩尔质量，单位为g/mol。计算真实气体平均分子自由程λ，如式(xv)所示：计算在有机质介质中的多孔介质平均knudsen数kn，如式(xvi)所示：计算在无机质介质中的多孔介质平均knudsen数kn，如式(xvii)所示：计算自由气单个孔隙气体流量qfree，如式(xviii)所示：式(xviii)中，无因次气体稀薄系数根据达西公式，自由气单个孔隙渗透率吸附气通过表面扩散方式进行流动，因此，计算真实气体表面扩散流量，如式(xix)所示：对于有机质孔隙，不同页岩岩石干酪根类型不同，气体吸附能力不同。孔隙材质对气体流动的影响可通过在模型中通过最大吸附浓度作为模型参数进行输入来进行考虑。计算最大吸附浓度cmax，通过实验室等温吸附曲线获得的langmuir体积vl进行计算，如式(xx)所示：式(xx)中，εks是指有机质占总岩石骨架体积比例；camax是指实验室测量得到的最大吸附气密度m3/mol,考虑气体覆盖度的影响，计算表面扩散系数ds，如式(xxi)所示：式(xxi)中，△h是指覆盖度“0”时的等量吸附热/(j·mol-1)；κm是指表面气体分子前进速度常数/m/s；κb是指表面气体分子阻塞速度常数，表面气体分子前进过程中由于其他气体分子的存在而退回原地的速度常数/m/s；κ是指表面气体分子阻塞系数,无因次；当κm＞κb，不论运动前方的位置是否已被分子占据，都将发生表面扩散；当κb＞κm，气体分子被堵塞，表面扩散停止。根据达西公式，计算吸附气单管渗透率ksurface，如式(xxii)所示：有机质孔隙中气体以吸附气和自由气两种形式流动，考虑孔隙度和迂曲度影响，计算页岩有机质介质渗透率kor，如式(xxiii)所示：无机质孔隙中气体主要以自由气流动，考虑孔隙度和迂曲度影响，计算页岩无机质介质渗透率kin，如式(xxiv)所示：本发明的有益效果为：1、本发明节省计算时间的同时保证了计算精度。传统高精度表观渗透率计算需要基于三维页岩数字岩心数据体直接模拟气体在孔隙相中的流动来计算不同温度压力下的表观渗透率，计算时间长，一般需要数周的时间，并且需要巨大的高性能计算机集群作为计算资源支撑。本发明通过引入多孔介质特征长度来准确反映三维页岩数字岩心数据体孔隙相结构对气体流动的影响，仅需要几分钟时间进行计算，大大简化了气体渗透计算过程的同时保证了计算的准确性。2、目前气体表观渗透率模型无法计算同时考虑有机质含量，有机质孔隙与无机质孔隙共存条件下的气体表观渗透率，并且忽略了有机质孔隙与无机质孔隙在气体传输、空间分布、孔隙半径分布上的差异。本发明综合考虑了有机质孔隙结构、无机质孔隙结构、有机质含量对有机质孔隙中自由气和吸附气流动的影响、多重气体运移机制、微纳米孔道中气体性质变化，能够准确预测页岩气实际生产条件气藏压力下降过程中的表观渗透率变化，为高效开发提供指导。3、本发明原始输入参数的页岩扫描电镜图像相比如目前气体表观渗透率模型需要的实验室测量获得的孔隙尺寸分布曲线以及三维纳米ct扫描页岩获得三维页岩数字岩心数据相比易获得、时间短、费用低、适合大规模推广应用。附图说明图1为本发明考虑多重因素影响的页岩气表观渗透率计算方法流程示意图；图2为实施例2中页岩无机质扫描电镜图像示意图；图3为实施例2中页岩无机质三维数字岩心示意图；图4为实施例2中页岩有机质扫描电镜图像示意图；图5为实施例2中页岩有机质三维数字岩心示意图；图6为弦长示意图；图7为三维页岩数字岩心切片弦长计算方法示意图；图8为页岩无机质孔隙表观渗透率在地层温度400k下随压力变化的示意图；图9为页岩有机质孔隙表观渗透率在地层温度400k下随压力变化的示意图。具体实施方式下面结合说明书附图和实施例对本发明作进一步限定，但不限于此。实施例1一种考虑多重因素影响的页岩气表观渗透率计算方法，包括步骤如下：(1)采用多点地质统计学方法，基于页岩二维剖光扫描电镜图像，重建页岩三维数字岩心；(2)计算页岩三维数字岩心的孔隙度、迂曲度、多孔介质特征长度；(3)计算页岩有机质和无机质中的气体渗透率。实施例2根据实施例1所述的一种考虑多重因素影响的页岩气表观渗透率计算方法，如图1所示。根据图2中页岩无机质扫描电镜图像(像素点个数：400×400,像素点分辨率：50nm)数值重构无机质三维数字岩心，具体步骤如下：a、将所选用的页岩无机质扫描电镜图像四个顶点中左上角顶点作为原点，从原点开始截取x方向的400个像素，y方向400个像素的正方形图像，作为要研究的训练图像，构建查找树；选取页岩无机质扫描电镜图像上若干孔隙相分布密集区域的像素点作为采样点，把采样点作为初始条件数据，随机分配到临近的网格节点上，网格节点是指页岩二维剖光扫描电镜图像上每个像素点的位置；b、对于任一待模拟网格节点u，确定在给定n个条件数据值s(uα)的情况下，属性s(u)取k个状态值中任一个状态值的条件概率如式(ⅰ)所示：式(ⅰ)中，sk指属性s(u)中的任一个状态值、skα指α个条件数据值下属性s(u)取k个状态值中任一个状态值条件数据事件；dn为由n个条件数据值s(uα)联合构成的数据事件，式(ⅰ)中的分母为某个模式出现的概率,α＝1,…,n；分子为某个模式出现和待模拟网格节点u取sk的情况同时出现的概率；对于所有的未模拟的网格节点，基于条件概率分布函数，使用蒙特卡洛方法提取模拟值，统计得到的模拟值个数，若获得的模拟值个数为400×400个，则将获得的400×400个模拟值作为新生成的二维图像，若获得的模拟值个数小于400×400个，则将得到的模拟值作为下一步的条件数据，继续模拟，直到获得的模拟值个数为400×400个；c、将步骤b新生成的二维图像放置于页岩无机质扫描电镜图像上方,得到垂直方向上高度为2层二维图像的三维图像,将步骤b新生成的二维图像作为下一步的训练图像，重复步骤a-b，直至获得基于训练图像得到的新生成的二维图像，继续放置于已构建的三维图像的上方，构建出垂直方向上高度为3层二维图像的三维图像；d、不断重复步骤a-c，直到获得399层生成的二维图像，二维图像中的每一个像素都对应着三维空间的相应的体素，加上步骤a选用的的页岩无机质扫描电镜图像，共有400层二维图像，构成了x、y、z三个方向大小为400×400×400的多孔介质三维图像，该多孔介质三维图像成为页岩无机质三维数字岩心，如图3所示，物理尺寸为20μm×20μm×20μm。同样的方法基于图4中页岩有机质扫描电镜图像(像素点个数：300×300,像素点分辨率：3.36nm)数值重构有机质三维数字岩心，如图5所示，物理尺寸为1.008μm×1.008μm×1.008μm。气体在三维多孔介质中运移时努森数kn表示为：λ为气体平均分子自由程；lc为多孔介质特征长度，基于弦长计算多孔介质特征长度，弦长的几何含义是指：用一条直线切割多孔介质，分布在孔隙相中不连续的直线。如图6所示闭合曲面表示固体骨架，分散在固体骨架间的直线算成为弦长。多孔介质特征长度lc通过弦长计算，如图7所示，如式(ⅰ)所示：式(ⅰ)中，z表示孔隙相中弦长长度，p(z)表示三维多孔介质中切割出的弦长中弦长长度为z的概率；根据弦长的物理意义，表征在页岩数字岩心数据体上即为连续的0数据。通过基于页岩无机质数字岩心和有机质数字岩心分别计算，得到页岩无机质多孔介质特征长度为131.4nm,页岩有机质多孔介质特征长度为9.82nm。通过基于页岩三维数字岩心模拟，迂曲度为：沿程路径与绝对路径之比；沿程路径为三维数字岩心中沿着孔隙相(数字0)与骨架相(数字1)的边界气体流过的长度；绝对路径指三维数字岩心边长；计算得到页岩无机质三维数字岩心迂曲度为2.3，有机质三维数字岩心迂曲度为1.8。孔隙度为：孔隙相与页岩三维数字岩心之比。计算得到页岩无机质三维数字岩心孔隙度为0.12，有机质三维数字岩心孔隙度为0.07。计算地层有效应力对多孔介质特征长度的影响，如式(ⅱ)所示：式(ⅱ)中，lc_stress是指考虑有效应力后的多孔介质特征长度；lc0为原始地层有效应力σeff0下的多孔介质特征长度，σeff为变化后的地层有效应力，p1为孔隙闭合压力，m为地层胶结指数；页岩中吸附气量采用langmuir等温吸附表示，计算吸附气在孔隙壁面的覆盖度θ，如式(ⅲ)所示：式(ⅲ)中，pl是指langmuir压力，单位为mpa；z是指气体压缩因子，无因次；θ是指真实气体孔隙壁面覆盖度，无因次；po是指地层压力，单位为mpa；计算考虑吸附气表面覆盖校正后的多孔介质特征长度lc_or，如式(ⅳ)所示：lc_or＝lc_stress-dmθ(ⅳ)式(ⅳ)中，dm是指气体分子直径，单位为m；计算微纳米孔道气体临界压力临界温度变化，通过式(ⅴ)、(vi)进行计算：式(ⅴ)、(vi)中，tc是指气体临界温度，单位k；pc是指气体临界压力，单位mpa；r是指理想气体常数，单位j/(mol·k)；σ是指lennard-jones作用尺寸，单位m；a是指甲烷范德华方程常数，单位mpa·kmol-2；b是指甲烷范德华方程常数，单位m6·mpa·kmol-2；ε是指lennard-jones势能参数，单位j-1；有效多孔介质特征长度lceff在无机质介质中为考虑有效应力后的多孔介质特征长度lc_stress，在有机质介质中为考虑吸附气表面覆盖校正后的多孔介质特征长度lc_or；计算气体压缩因子z，如式(ⅶ)、式(ⅷ)、式(ⅸ)所示：ppr是指拟压力，无因次；tpr是指拟温度，无因次；t是指地层温度，单位为k；计算粘度，如式(ⅹ)、式(ⅺ)、式(ⅻ)、式(xiii)、式(xiv)所示：μ＝(1×10-4)kexp(xρy)(ⅹ)y＝2.447-0.2224x(xiv)式(ⅹ)、式(ⅺ)、式(ⅻ)、式(xiii)、式(xiv)中，μ是指气体粘度，单位为pa·s,x、y、k为计算气体粘度公式中的相关参数，受地层温度和气体摩尔质量影响。m是指气体摩尔质量，单位为g/mol。计算真实气体平均分子自由程λ，如式(xv)所示：计算在有机质介质中的多孔介质平均knudsen数kn，如式(xvi)所示：计算在无机质介质中的多孔介质平均knudsen数kn，如式(xvii)所示：计算自由气单个孔隙气体流量qfree，如式(xviii)所示：式(xviii)中，无因次气体稀薄系数根据达西公式，自由气单个孔隙渗透率吸附气通过表面扩散方式进行流动，因此，计算真实气体表面扩散流量，如式(xix)所示：对于有机质孔隙，不同页岩岩石干酪根类型不同，气体吸附能力不同。孔隙材质对气体流动的影响可通过在模型中通过最大吸附浓度作为模型参数进行输入来进行考虑。计算最大吸附浓度cmax，通过实验室等温吸附曲线获得的langmuir体积vl进行计算，如式(xx)所示：式(xx)中，εks是指有机质占总岩石骨架体积比例；camax是指实验室测量得到的最大吸附气密度m3/mol,考虑气体覆盖度的影响，计算表面扩散系数ds，如式(xxi)所示：式(xxi)中，△h是指覆盖度“0”时的等量吸附热/(j·mol-1)；κm是指表面气体分子前进速度常数/m/s；κb是指表面气体分子阻塞速度常数，表面气体分子前进过程中由于其他气体分子的存在而退回原地的速度常数/m/s；κ是指表面气体分子阻塞系数,无因次；当κm＞κb，不论运动前方的位置是否已被分子占据，都将发生表面扩散；当κb＞κm，气体分子被堵塞，表面扩散停止。根据达西公式，计算吸附气单管渗透率ksurface，如式(xxii)所示：有机质孔隙中气体以吸附气和自由气两种形式流动，考虑孔隙度和迂曲度影响，页岩有机质介质气体表观渗透率kor，如式(xxiii)所示：无机质孔隙中气体主要以自由气流动，考虑孔隙度和迂曲度影响，页岩无机质介质气体表观渗透率kin，如式(xxiv)所示：根据式(xxiii)和式(xxiv)分别计算页岩无机质数字岩心、有机质数字岩心表观渗透率，模型发明输入参数见表1。表1计算得到的不同地层压力下页岩无机质数字岩心表观渗透率见图8，页岩有机质数字岩心表观渗透率见图9。从图8和图9可以看出，随着地层压力升高，页岩无机质介质气体表观渗透率和有机质介质气体表观渗透率均有所降低。页岩无机质介质气体表观渗透率变化范围在2.3×10-5μm2到2.5×10-5μm2之间，页岩有机质介质气体表观渗透率变化范围在7×10-7μm2到1.25×10-6μm2之间。为验证本发明的效果及相比于目前方法的优势，将典型地层压力40mpa下本发明计算得到的无机质介质气体表观渗透率计算结果与直接基于图3页岩无机质数字岩心气体流动模拟(以格子波尔兹曼方法为例)得到的无机质介质气体表观渗透率计算结果进行对比。格子波尔兹曼方法计算气体表观渗透率技术参考文献为[zhaoj,yaoj,zhangm,zhangl,yangy,sunh,ans,lia.studyofgasflowcharacteristicsintightporousmediawithamicroscalelatticeboltzmannmodel.scientificreports.2016sep2；6:32393.]。计算结果见表2，表2为典型地层压力40mpa下本发明气体表观渗透率计算结果与直接基于三维页岩数字岩心数据体流动模拟结果(以格子波尔兹曼方法为例)对比表。表2页岩有机质数字岩心页岩无机质数字岩心基于三维页岩数字岩心数据体流动模拟计算时间t/s8.1×1051.92×106直接流体模拟气体表观渗透率计算结果k/μm27.33×10-72.42×10-5本发明计算时间t/s405960本发明气体表观渗透率计算结果k/μm27.41×10-72.35×10-5从表2可以看出本发明大大简化了气体渗透计算时间的同时保证了计算的准确性。当前第1页12