：
1.本发明涉及石油天然气开发领域，具体涉及一种页岩油藏启动压力梯度的计算方法。
背景技术：2.近年来，世界对油气资源的需求增长迅速，油气资源的开发开采逐渐从常规油气储层转向为非常规油气储层。页岩油藏分布广泛，但储层基质普遍存在启动压力梯度，因此亟需对页岩油藏的启动压力梯度展开理论研究。
3.页岩油储层岩石矿物以石英、长石、粘土矿物、碳酸盐矿物为主，脆性矿物含量高，有机质发育，属于低孔低渗储层，孔隙结构复杂多样，储集空间为无机孔、有机孔和裂缝，其中无机孔包括粘土孔和脆性矿物孔，页岩油在各类孔隙结构中不会立即流动，必须当驱动压力梯度达到阈值时流动才会发生，此时的阈值就被称为启动压力梯度。目前页岩油藏基质启动压力梯度研究主要以实验为主，通过对不同岩心进行实验测定，拟合回归得到含水饱和度与启动压力梯度的经验公式，单一将启动压力梯度的成因归结于储层中的束缚水，无法揭示启动压力梯度的产生机理和影响因素；虽然部分理论模型从流体自身受力平衡出发，建立了启动压力梯度的机理模型，但没有综合考虑基质不同孔隙类型(有机孔、无机孔)中流体的流动特点。
技术实现要素：4.本发明的目的是为了克服现有方法中的不足，提供一种页岩油藏启动压力梯度的计算方法，克服了现有技术中没有考虑基质不同孔隙类型(有机孔、无机孔)中流体流动的缺陷。
5.本发明的目的是通过以下技术方案来实现的。
6.本技术提供了一种页岩油藏启动压力梯度的计算方法，包括步骤：
7.(1)以页岩油藏基质单毛管为研究对象，引入有效滑移长度表征孔隙类型，考虑真实粘度的影响，建立单毛管的页岩油流动受力平衡方程；
8.(2)基于步骤(1)中建立的受力平衡方程，获取单毛管的流速方程以及流量方程，通过分形理论进行尺度升级，获取得到岩心尺度的流量方程；
9.(3)根据岩心尺度的水相流量方程，令流量为0，计算得到对应的页岩油藏启动压力梯度。
10.进一步地，有效滑移长度表达式为
[0011][0012]
式中，l
se
为有效滑移长度；l
st
为真实滑移长度；l
sa
为表观滑移长度；μ
b
为体相流体粘度；μ
d
为有效粘度；λ为毛管直径；
[0013][0014]
式中；c为液体常数；θ为润湿角。
[0015]
进一步地，有效滑移长度为正值表示孔隙类型为无机孔，有效滑移长度为负值表示孔隙类型为有机孔。
[0016]
进一步地，有效粘度获取方法为：
[0017]
考虑体相流体粘度和核心层
‑
边界层界面面积的影响，对体相流体粘度和边界层流体粘度进行加权获取有效粘度：
[0018]
计算式为
[0019][0020]
其中：
[0021]
a
id
＝π[(λ/2)2‑
(λ/2
‑
δ)2]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0022]
a
td
＝π(λ/2)2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0023]
式中，μ
d
为有效粘度；μ
nw
为边界层流体粘度；a
id
为边界层流体区域所占面积；a
td
为毛细管总横截面积；μ
b
为体相流体粘度；δ为边界层厚度。
[0024]
进一步地，边界层流体粘度可表示为：
[0025][0026]
进一步地，体相流体粘度可表示为：
[0027][0028]
式中：t为储层温度，k。
[0029]
进一步地，受力平衡方程为：
[0030][0031]
式中，
△
p为驱动压差；r为毛细管半径；l
se
为有效滑移长度；η0为极限剪应力；μ
d
为流体有效粘度；v为流体流动速度；l
t
为毛细管实际长度。
[0032]
进一步地，所述单毛管的流速方程表达式为
[0033][0034]
式中，dt为迂曲度分形维数；l0为毛细管直线距离；
[0035]
[0036]
式中，r为滑移边界距离。
[0037]
进一步地，单毛管的流量方程表达式为
[0038][0039]
进一步地，岩心尺度的流量方程为
[0040][0041]
式中，n为分段积分的分段数；下标max为最大值；下标min为最小值；i为第i段；d
f
为孔隙分形维数。
[0042]
进一步地，页岩油藏启动压力梯度表达式为：
[0043][0044]
其中
[0045][0046][0047]
式中：g为页岩油藏启动压力梯度。
[0048]
与现有技术相比，本发明带来的有益技术效果包括：
[0049]
(1)本发明解决了页岩油藏储层孔类型多样以及启动压力梯度计算困难等特点。
[0050]
(2)本发明根据页岩储层基质同时存在有机孔、无机孔以及各自的润湿性不同，引入有效滑移长度来描述不同润湿性对页岩油流动的影响，充分考虑页岩油自剪切流动、有效滑移长度以及真实粘度的影响，计算出页岩油藏基质启动压力梯度。
[0051]
(3)本发明对单毛管采用分段积分的方式获取流量，能够解决粘度随单毛管长度方向动态变化的问题。
附图说明
[0052]
图1为无机孔启动压力梯度与润湿角的变化关系曲线图。
[0053]
图2为有机孔启动压力梯度与含水饱和度的变化关系曲线图。
[0054]
图3为单毛管启动压力梯度与孔隙度的变化关系曲线图。
[0055]
图4为单毛管启动压力梯度与孔隙直径的变化关系曲线图。
[0056]
图5为单毛管启动压力梯度与滑移长度的变化关系曲线图。
[0057]
图6为基质尺度启动压力梯度与孔隙度的变化关系曲线图。
[0058]
图7为基质尺度启动压力梯度与最大孔径的变化关系曲线图。
[0059]
图8为基质尺度启动压力梯度与滑移长度的变化关系曲线图。
具体实施方式
[0060]
以下将结合实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成相应技术效果的实现过程能充分理解并据以实施。本技术实施例以及实施例中的各个步骤，在不相冲突前提下可以相互结合，所形成的技术方案均在本发明的保护范围之内。
[0061]
另外，本技术中给出了一些实施例以及相应步骤，但是在某些情况下，可以以不同于本技术中采用的理论模型来建立最终的启动压力计算模型，例如可以采用现有技术中的其他单毛管流量计算方法以及粘度描述方程等。
[0062]
计算分析页岩油藏基质孔隙有效滑移长度和真实粘度，基于综合滑移边界条件，采用真实滑移方程可以有效表征流体分子在毛细管壁面上的滑移，分析页岩毛细孔的有效滑移长度以及真实粘度的推导过程；
[0063]
(1.1)有效滑移长度
[0064]
为表征毛管中流体流动的有效滑移，引入真实滑移长度和有效滑移长度。
[0065]
真实滑移长度可表示为：
[0066][0067]
式中，l
st
为流体真实滑移长度，nm；c为液体常数，取0.41；θ为润湿角，
°
。；
[0068]
有效滑移长度可表示为：
[0069][0070]
式中，l
se
为有效滑移长度，nm；l
st
为真实滑移长度，nm；l
sa
为表观滑移长度，nm；μ
b
为体相流体粘度，mpa
·
s；λ为毛细管直径，nm；μ
d
为有效粘度，mpa
·
s。
[0071]
(1.2)有效粘度
[0072]
有效粘度上取决于体相流体粘度和核心层
‑
边界层界面面积，对体相流体粘度和边界层流体粘度进行加权平均可得：
[0073][0074]
其中：
[0075]
a
id
＝π[(λ2)2‑
(λ/2
‑
δ)2]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0076]
a
td
＝π(λ/2)2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0077]
式中，μ
d
为有效粘度，mpa
·
s；μ
nw
为边界层流体粘度，mpa
·
s；a
id
为边界层流体区域
所占面积，nm2；δ为边界层厚度，nm，取0.7nm；μ
b
为体相流体粘度，mpa
·
s；a
td
为毛细管总横截面积，nm2。
[0078]
边界层流体粘度可表示为：
[0079][0080]
体相流体粘度可表示为：
[0081][0082]
式中：t为储层温度，k。
[0083]
(2)基于极限剪应力理论，可以建立页岩油藏基质毛细管流动方程和启动压力梯度模型。
[0084]
根据牛顿内摩擦定律，同时引入有效滑移长度并建立页岩油流动的受力平衡方程：
[0085][0086]
其中可用迂曲度分形来描述真实毛细管长度，表达式如下：
[0087][0088][0089]
特别地，对于毛细管流动过程而言，有效滑移长度相当于增加或减少了毛细管有效流动半径，这与毛细管的润湿性有关，根据不同孔隙类型润湿性的特点，有效滑移长度为正值表示此时的孔隙为无机孔，相反则为有机孔。
[0090]
式中，
△
p为驱动压差，mpa；r为毛细管半径，nm；l
se
为有效滑移长度，nm；η0为极限剪应力，mpa；μ
d
为流体有效粘度，mpa
·
s；v为流体流动速度，nm/s；l
t
为毛细管实际长度，nm；l
t
为毛细管实际距离，nm；l0为毛细管直线距离，nm；dt为迂曲度分形维数，1<dt<2；τ
av
为毛细管平均迂曲度，无量纲；λav为毛细管平均直径，nm。
[0091]
滑移边界条件如下：
[0092][0093]
将式(9)、(10)、(11)代入式(8)积分可得毛细管流速分布为：
[0094][0095]
其中：
[0096][0097]
对单毛管整个流动区域进行积分可得体积流量为：
[0098][0099]
当流体刚发生流动时，即此时的驱动压力梯度达到阈值，此时的压力梯度为启动压力梯度，令q(r)＝0，可得毛细管启动压力梯度为：
[0100][0101]
其中：
[0102][0103]
a3＝d
t
(d
t
+1)(d
t
+2)(d
t
+3)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0104]
式中：g为单毛细管启动压力梯度，mpa/m。
[0105]
计算页岩油藏基质启动压力梯度，通过分形理论尺度升级，同时通过分段积分抵消真实粘度的动态变化，具体为：
[0106]
(2.1)分形理论
[0107]
根据分形理论，储层孔隙直径≥λ的单元体横截面毛管总数n为：
[0108][0109]
对式(17)中λ求导可得：
[0110][0111]
设λmin/λmax＝β，则df可以表示为：
[0112][0113]
式中：d
f
为孔隙分形维数，无因次；d为欧几里得维数，d＝2；为孔隙度，％。
[0114]
单元体横截面可视为不同直径毛细管组成，则可求出单元体横截面积如下：
[0115][0116]
(2.2)分段积分
[0117]
考虑有效粘度随管径的变化，为了方便求解，将页岩油储层所有毛细管直径离散
化为j个微小单元，第i段(r
min,i
≤r
i
≤r
max,i
)的单毛细管流量可以写成下式：
[0118][0119]
对式(21)在(r
min,i
≤r
i
≤r
max,i
)上积分得第i段的总流量为：
[0120][0121]
其中
[0122][0123]
a
id,i
＝π[r
av,i2
‑
(r
av,i
‑
δ
i
)2]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(24)
[0124]
a
td,i
＝πr
av,i2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0125][0126]
m
i
＝r
min,i
/r
max,i
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0127]
式中：q
i
为第i段毛细管总流量，nm3/s；μ
d,i
(r
av,i
)为第i段毛细管平均粘度，mpa
·
s；a
id,i
为第i段毛细管边界层流体区域所占面积，nm2；a
td,i
为第i段毛细管总横截面积，nm2；r
av,i
为第i段毛细管平均半径，nm；m
i
为第i段毛细管最小最大孔隙半径比，无因次；δ
i
为第i段毛细管边界层厚度，nm。
[0128]
当然在计算过程中，本领域技术人员也可忽略粘度随毛管长度方向的变化情况，例如在原油粘度变化不大的情况下，本领域技术人员可以设置粘度为某一定值，从而忽略上述过程中的粘度变化描述过程。
[0129]
(3)根据岩心尺度的水相流量方程，令流量为0，计算得到对应的页岩油藏启动压力梯度。
[0130]
通过代数叠加每个微小单元流量，得到总流量如下：
[0131][0132]
式中：n为分段积分的分段数，无因次。
[0133]
将式(22)代入到式(27)可得分段积分后的总流量方程为：
[0134][0135]
在式(28)中令q＝0可得页岩油藏基质启动压力梯度如下：
[0136][0137]
其中
[0138][0139][0140]
式中：g为基质启动压力梯度，mpa/m。
[0141]
各项基本参数如表1所示：
[0142]
表1各项基本参数
[0143][0144][0145]
(2)计算结果
[0146]
由图1可知：随着润湿角的增加，无机孔单毛管启动压力梯度逐渐增加，孔隙直径为500nm时，当润湿角从3
°
增加到90
°
，无机孔单毛管启动压力梯度从0.08mpa/m增加到65.04mpa/m；同时在同一润湿角下，孔隙直径越大，无机孔单毛管启动压力梯度越小，不同
孔隙直径下无机孔单毛管启动压力梯度之间的差异在润湿角较大的时候比较大。
[0147]
由图2可知：有机孔单毛管启动压力梯度随着润湿角的增加而大幅增加，这是因为润湿角越大，储层流体与毛细管壁之间的相互作用越强，有效滑移效应越弱，故启动压力梯度越高。
[0148]
由图3可知：孔隙度对单毛管启动压力梯度的影响分为两个阶段，在孔隙度达到3％以前，单毛管启动压力梯度随孔隙度的增加而快速下降，孔隙度超过3％以后，随着孔隙度的增加，启动压力梯度下降速度逐渐减缓；同时极限剪切应力对启动压力梯度也有影响，极限剪切应力为0时，单毛管启动压力梯度为0，在相同条件下，极限剪切应力越大，单毛管启动压力梯度越大。
[0149]
由图4可知：随着孔隙直径的增加，单毛管启动压力梯度先快速下降再缓慢降低，极限剪切应力为0.2pa时，当孔隙直径从20nm增加到200nm，单毛管启动压力梯度下降了95％，当孔隙直径从200nm增加到600nm，单毛管启动压力梯度下降了68％。
[0150]
由图5可知：随着有效滑移长度的增加，单毛管启动压力梯度逐渐下降，最大孔隙直径为500nm时，当有效滑移长度从0.1nm增加到3nm时，单毛管启动压力梯度从8.78mpa/m下降到5.98mpa/m，下降幅度为32％。
[0151]
由图6可知：孔隙度对基质启动压力梯度的影响同样分为两个阶段，在孔隙度达到3％以前，基质启动压力梯度随孔隙度的增加而迅速下降，孔隙度超过3％以后，随着孔隙度的增加，启动压力梯度下降速度逐渐减缓；总体而言，基质启动压力梯度下降速度比单毛管快，同时，相同条件下基质启动压力梯度比单毛管启动压力梯度大。
[0152]
由图7可知：随着最大孔隙直径的增加，基质启动压力梯度先快速下降再缓慢降低，极限剪切应力为0.2pa时，当最大孔隙直径从20nm增加到200nm，基质启动压力梯度下降了94％，当最大孔隙直径从200nm增加到600nm，基质启动压力梯度下降了77％。
[0153]
由图8可知：随着有效滑移长度的增加，基质启动压力梯度逐渐下降，最大孔隙直径为500nm时，当有效滑移长度从0.1nm增加到3nm时，基质启动压力梯度从20.16mpa/m下降到8.11mpa/m，下降幅度达60％。
[0154]
以上通过实施例对本发明进行具体描述，有必要在此指出的是，本实施例仅是本发明的优选实施例，并非对本发明作任何限制，也并非局限于本文所披露的形式，不应看作是对其他实施例的排除。而本领域人员所进行的改动和简单变化不脱离本发明技术思想和范围，则均属于本发明技术方案的保护范围内。