我国油气资源总量相对不足、分布不均匀、成藏条件复杂。成熟油田二次开采和非常规储层勘探开发支撑着我国油气总量的持续增长, 提高致密储层(油水混合、低孔低渗)油气采收率是挖掘我国油气资源潜力需要解决的关键问题之一。

致密储层开采面临的最大问题是如何通过大型水力压裂技术实现经济产量, 导致该问题的直接原因是储层油气流动性差, 根本原因在于压力作用下孔隙/裂缝中流体流动机理尚不清楚。一般认为, 致密气藏基质渗透率小于0.1mD(1mD≈0.987×10^-3μm²), 基质孔隙度小于10%^[1-2], 致密气砂岩的孔喉大小为0.03~2.00mm^[3], 致密油藏的孔喉直径通常小于1μm^[2]。实验观测和利用数字岩心技术得到的研究成果表明, 致密储层岩石孔隙/裂缝达到微纳米尺度, 并形成复杂网络结构。当有机物分子大小与孔隙尺寸相当时, 流体流动性依赖于孔隙空间大小和形状。此外, 粘稠原油、油水混合物等复杂流体的渗透率还与流体流变学特征相关^[4]。传统孔隙介质模型基于球形孔或者单根毛细管假设, 无法描述致密储层孔隙-裂缝网络的拓扑连通结构, 也不能用于解释流体流动的非线性现象。

提高储层采收率技术的理论基础是孔隙介质力学, 其核心为孔隙/裂缝结构的渗透率模型。传统孔隙介质流体流动模型主要包括单毛细管模型、毛细管束模型和毛细管网络模型。单毛细管模型假设孔隙空间是一根微型圆管道, 流体在圆管道中发生泊肃叶流动。毛细管流体动力学是一个经典的研究领域。HAGENBACH^[5]给出了单根管道中泊肃叶流动的解析表达式(哈根-泊肃叶公式); WILBERFORCE^[6]研究了毛细管流体流动粘性问题; WASHBURN^[7]研究了静压力下毛细流动现象; BIOT^[8-11]基于单毛细管流动模型, 进行了重要的孔隙介质动力学理论研究, 建立了完全饱和孔隙介质中的波传播方程, 预测了两类纵波和一类横波。近50年来, Biot理论不断发展, 并且在各类流体饱和孔隙介质中得到广泛应用。在毛细管尺度远小于波长的前提下, 可以近似地认为流体在岩石内部“均匀”分布并且全局流动, 由地震波场导致的流体粘性损耗也“均匀”一致。随后的研究发现, 当探测频率超过100kHz时, 岩石中全局流动引起的衰减较明显^[12]。因单毛细管模型不考虑孔隙空间的复杂几何和拓扑结构, 故不适用于低孔低渗储层岩石。

毛细管束模型假设孔隙空间由一束均匀圆管组成^[13], BARTLEY等^[14]研究了均匀分布的毛细管束中不混溶流体流动的解析计算方法, 并利用该方法计算了相对渗透率。WANG等^[15]推导了相互作用毛细管束模型的相对渗透率解析表达式, 通过调整模型参数估算出与产油历史和压降相匹配的相对渗透率曲线。毛细管束模型虽然适用于满足达西定律的静态渗透率预测, 但是不适用于刻画孔隙-裂缝之间纵横连通的拓扑结构和天然岩石中具有复杂连通关系的孔隙-裂缝结构。

利用岩心切片电子显微镜观测和数字岩心技术得到的最新研究成果显示, 在岩石内部普遍存在裂缝-孔隙的连通网络。孤立孔隙和裂缝模型由于未考虑孤立孔隙/裂缝之间流体连通通道产生的影响, 故对地震波场速度存在一定的影响, 但是这些孔缝对储层流体渗透率的影响可以忽略不计。平行排列的裂缝系统能够较好地反映地层各向异性属性^[16-17], 但是很难描述跨平行裂缝层的流体流动现象。因此, 相对于研究局部孔隙-裂缝相互作用的模型, 毛细管网络模型, 尤其是三维裂缝-孔隙网络模型是一种更合理的描述孔隙连通性、流体流动性和骨架弹性模量的模型。

FATT^[18]较早提出了孔隙介质简化网络模型, 并利用规则网络上的管道表示孔隙空间。BRYANT等^[19]提出了一个多孔介质的无序网络模型, 认为真实多孔系统中孔隙尺寸随机分布的假设可能无效。SEEBURGER等^[20]提出将二维孔隙空间网络模型用于砂岩的静渗透率预测, 该模型应用于砂岩和花岗岩的计算结果表明, 网络模型可以根据围压函数预测渗透率大小。基于流体力学基本方程, BERNABÉ^[21]建立了牛顿流体在刚性和弹性管网中的水传导率和波频散方程, 随后, BERNABÉ^[22]通过管道和管网对牛顿流体的输运特性进行了数值模拟。JOHNSON等^[23]得到了不同管道长度和半径的网络模型中水力传导率与频率的关系。BLUNT等^[24-25]研究了由细管形成的网络中相对渗透率与饱和度的关系, 证明了网络模型能够解释微观流动如何影响厘米尺度上的网络平均特性。考虑到孔隙尺度的润湿性变化, 三维孔隙网络也可用于模拟砂岩内部排水和注水过程^[26-27]。大量研究表明, 虽然局部非均匀孔隙和裂缝参数对孔隙流体会产生一定影响, 但在孔隙网络节点相互连通的“均质化”效应下, 其对储层岩石整体渗透率的影响十分有限。因此, 相对于研究局部孔隙与裂缝相互作用的模型来说, 更需要将跨微观-宏观尺度的三维裂缝/孔隙网络系统作为研究对象, 从整体角度考察裂缝网络对岩石孔隙度、渗透率、地震波频散衰减以及储层采收率等诸多方面的影响。XIONG等^[4]分析了任意三维孔隙网络中流体流变学特征和孔隙连通关系对岩石渗透率的影响, 发现渗透率变化受脉冲压力波场频率、流体流变学属性和孔隙连通性的影响, 流体在低频时主要在稀疏大孔隙网络中流动, 而在高频时则以稠密小孔隙网络为主要流动通道^[4]。麦克斯韦流体在特定频率附近出现渗透率峰值, 峰值位置与网络结构平均配位数(mean coordination number, MCN)相关, 这表明在特定频率范围下, 存在最优的网络连接关系, 在这种连接关系下流体流动最容易。

非常规储层具有裂缝/孔隙结构复杂、孔隙流体部分饱和、非达西流动等特点, 当孔隙空间达到微纳米尺度时, 传统孔隙介质模型的许多基本假设不适用。渗透性通常被认为是岩石本身的一种性质, 不随压力和内部流动的气体类型而变化。然而, 这种假设在气体滑移现象显著的情况下不成立。首先, 在致密储层页岩气勘探开发中, 考虑到纳米尺度毛细管对气体平均自由程的影响, 需要引入Knudsen扩散效应模型。由于Knudsen扩散效应的影响, 导致气体渗透率增加, 经Klinkenberg修正和Knudsen修正后, 常规流动方程在一定条件下仍然可以使用。其次, 在注液压裂过程中, 孔隙流体可能出现较高的流速, 引起非达西流动, 此时需要引入Forchheimer模型修正流体流动计算结果。国内外学者对致密储层非达西流动理论模型进行了深入研究^[28-29], 由于这些复杂流体流动模型较为特殊, 其理论基础和适用范围的研究十分关键, 对于推动理论模型在实际工作中的应用至关重要, 故本文针对微纳米孔隙复杂流体流动模型及其适用性进行研究。本文详细分析了微纳米孔隙条件下复杂流体流动的Forchheimer流动模型、Klinkenberg修正模型、对Knudsen模型修正得到的B-K模型和高阶修正模型, 通过对理论模型进行数值计算和结果对比分析, 得到了不同孔隙尺度、流速和压力下孔隙流体满足各种模型的条件, 有助于判断复杂条件下不同流动模型的适用性。研究发现, Knudsen模型适用于研究微纳米孔隙中较低气压下气体流动, 对于孔隙特征尺寸很小的低孔渗储层, 流体流速通常很低, Forchheimer模型的修正效应并不明显。

1 方法原理 1.1 微纳米孔隙流动类型

当流体在微纳米流道中流动时, 如果流道的特征尺寸(如流道直径)接近甚至小于流体分子间的平均自由程, 那么相较于流体分子之间的碰撞, 分子与流道固壁的碰撞概率极大地增加, 这种流体扩散方式称为Knudsen扩散。

在Knudsen扩散中, 流体分子间碰撞极大地降低, 分子的运动轨迹为流道表面不同碰撞点之间的直线运动。Knudsen扩散的条件是流动空间特征尺寸小于分子平均自由程。分子平均自由程与温度和压力有关, 较低压力下气体更加稀薄, 分子平均自由程更大。对于速度满足麦克斯韦分布的气体来说, 平均自由程可以表示为^[30]:

$ l = \frac{{{k_{\rm{B}}}T}}{{\sqrt 2 {\rm{ \mathsf{ π} }}d_{\rm{g}}^2p}} $

(1)

式中:l是平均自由程; k_B是波尔兹曼常数; T是温度; p是压力; d_g是气体分子直径。

表 1为不同气体分子的碰撞直径^[31-33], 表 2为室温下空气分子在不同压力下的平均自由程^[34]。Knudsen数(K_n)用于判断是否发生Knudsen流动, K_n定义为分子平均自由程l与流动空间特征尺度Λ之比:

$ {K_{\rm{n}}} = \frac{l}{\varLambda } $

(2)

根据Knudsen数可选择使用连续介质流体力学模型或者统计力学模型。当Knudsen数接近或大于1时, 连续介质模型的适用性大大降低。

对液体和稠密气体而言, 平均自由程与分子的大小相当, 因此对于微纳米孔隙中的致密油而言, 只有当孔隙特征尺寸降到纳米甚至纳米以下时, 才会出现明显的Knudsen流动效应。根据Knudsen数, 可以将孔隙介质中气体流动大致分为4种情况^[35-36], 对于每种流动情况下的Knudsen数范围, 研究人员提出的划分范围略有差异。ROY等^[36]提出连续介质无滑移流动时, K_n < 0.001;滑移流动时, 0.001 < K_n < 0.100, 无滑移流动和滑移流动均满足纳维叶-斯托克斯方程。ZIARANI等^[35]提出连续介质滑移流动时, 0.01 < K_n < 0.10。本文采用ZIARANI等^[35]提出的划分标准, 4种不同的流动情况分别介绍如下。

1.1.1 连续介质无滑移流动(K_n < 0.01)

在这种情况下, 孔隙介质中流体流动速度较慢, 粘性项作用大于惯性项, 流体流动满足层流假设, 流量和压力梯度之间遵循达西定律:

$ {q_{{\rm{Darcy }}}} = \frac{{{\rm{ \mathsf{ π} }}{r^4}}}{{8\mu }}\frac{{\Delta P}}{L} $

(3)

式中:q_Darcy是流体流量; r是孔喉半径; μ是流体动力学粘性系数; L是流道长度; ΔP是流道两端压力差。

当流体流速较大时, 惯性项的作用逐渐增加并打破了层流假设, 此时需要对达西定律进行修正, Forchheimer模型是一种重要的修正模型。

1.1.2 滑移流动(0.01≤K_n≤0.10)

在这种情况下, 气体分子在固体表面发生滑移运动, 若适当修正模型, 常规流动方程仍然可以使用, 但是在气体-固体边界上该方程存在速度和温度的不连续间断。修正模型主要包括Klinkenberg修正模型和Knudsen修正模型。通常认为, Knudsen修正模型更准确, 它包括B-K二阶修正模型^[37]和高阶修正模型^[38]。Klinkenberg修正模型更易于使用, 因此被广泛应用。当K_n=0.10或K_n=0.01时, 即处于连续介质滑移流动区间的上、下边界, 亦可以利用Klinkenberg修正模型和Knudsen修正模型描述流体运动流动情况。

1.1.3 瞬变流动(0.1＜K_n≤10.0)

在这种情况下, 滑移(连续介质)流动和扩散流动同时并存, 但是传统流体动力学方程失效。尽管经过修正的达西定律在一定程度上仍然可用, 但是其结果存疑, 比较可靠的方法是使用Knudsen扩散方程或带滑移边界条件的伯内特方程^[39]来描述流体运动。

1.1.4 Knudsen(自由分子)流动(K_n>10)

这种情况下, 连续介质流体力学方程完全失效, 需要使用基于扩散模型的流动方程。只有在压力极低和孔喉十分致密的情况下才会发生Knudsen流动, 此时, 低密度流体分子之间几乎不发生碰撞。

1.2 Forchheimer模型

雷诺数是用来判断粘性流体流动状态的无量纲数, 其物理含义可以理解为惯性力和粘性力的对比关系。当孔隙介质中的流体因雷诺数低而产生粘性流动时, 可以用达西定律描述流量(流速)与压力梯度之间的线性关系。当流体的雷诺数较高时, 流体的惯性效应超过了粘性效应, 导致了湍流现象, 这种情况下达西定律失效。1901年, FORCHHEIMER^[40]提出了一个达西定律的修正公式, 并引入了速度的平方项作为达西定律的修正项, 该公式被称为Forchheimer公式, 具体如下:

$ - \nabla P = \frac{{\mu v}}{\kappa } + \beta \rho {v^2} $

(4)

式中:P是压力; v是流速; ρ是密度; κ是渗透率; β是孔隙介质的非达西系数。当地层中流体的流速很高时, (4)式等号右端第2项对速度的二次修正影响大, 此时Forchheimer模型对流体流动的修正效果明显。

该模型中的系数β又被称为惯性阻力系数, 在实际工作中通常根据Forchheimer图(孔隙介质的视渗透率倒数(1/κ_app)随伪雷诺数(ρv/μ)的变化曲线)通过斜率拟合得到。视渗透率κ_app定义为:

$ \frac{1}{{{\kappa _{{\rm{app}}}}}} = - \frac{{\nabla P}}{{\mu v}} = \frac{1}{\kappa } + \beta \frac{{\rho v}}{\mu } $

(5)

由(5)式可以看出, 1/κ_app跟ρv/μ之间满足线性关系, β是线性关系的斜率。通常可以依据雷诺数和Forchheimer数的大小, 判断连续介质流体流动是否产生非达西流动, 雷诺数R_e和Forchheimer数F₀分别表示如下:

$ {{R_{\rm{e}}} = \frac{{\rho dv}}{\mu }} $

(6)

$ {{F_0} = \frac{{\kappa \beta \rho v}}{\mu }} $

(7)

式中:ρ是流体密度; d是流道特征长度; v是流速; μ是动力学粘性系数。可以看出, 如果将κβ视作特征长度, 那么F₀即为另一种形式的雷诺数。

根据流体类型和孔隙介质的性质, 并基于实验数据^[41-42], 人们发现了不同的流动形态, 并依据雷诺数加以分类。

1) 达西流动状态:粘性力占主导地位, 雷诺数满足R_e＜1, 当R_e≈1时, 孔隙介质固体表面逐渐形成边界层。

2) 惯性流动状态:惯性流动状态从R_e=1~10开始, 持续到R_e=150, 孔隙介质固体表面边界层更加明显, 流体出现“惯性核”。

3) 不稳定层流状态:当R_e=150~300时, 流体开始出现弱震荡, 湍流开始发展。

4) 高度不稳定流动和混乱流动:当R_e＞300时, 出现类似管道流动中的湍流现象, 涡流逐渐占据主导地位。

(4) 式等号右端第2项对速度的二次修正使流体达西流动模型变为Forchheimer流动模型。真实流体从达西流动到Forchheimer流动的转变约始于雷诺数R_e=1, 对于是否存在一个临界雷诺数, 目前尚存在不同观点。一种观点认为, 当达西流动的线性部分引起的压力下降所占比例低于95%时, 可以认为是达西流动状态的结束阶段, 此时雷诺数R_e≈4.3^[43]。数值计算结果显示, 在雷诺数R_e=2~4时, 发生了从达西流动到Forchheimer流动的转变^[44]。

对于理想圆柱形直管道来说, 可以采用上述方法计算雷诺数。但是, 对于真实致密储层岩石来说, 孔隙结构不满足圆柱形管道的理想假设, 因此, 对于含有复杂孔隙的岩石, 采用上述雷诺数计算方法会带来难以估计的误差。对于弯曲孔隙的毛细作用模型, COMITI等^[45]给出了含有复杂孔隙结构的雷诺数定义, 并得到一个形似Forchheimer公式的方程:

$ R_{\mathrm{e}}=\frac{4 \rho v \tau}{(1-\varphi) \mu a_{\mathrm{vd}}} $

(8)

$ \frac{{\Delta P}}{H} = 2\gamma {\tau ^2}\mu a_{{\rm{vd}}}^2\frac{{{{(1 - \varphi )}^2}}}{{{\varphi ^3}}}v + 0.096{\kern 1pt} {\kern 1pt} {\kern 1pt} 8{\tau ^3}{a_{{\rm{vd}}}}\frac{{1 - \varphi }}{{{\varphi ^3}}}\rho {v^2} $

(9)

式中:ΔP为孔隙流体压力差; H为ΔP的等效距离; φ是孔隙度; τ是弯曲度; a_vd是动力学比表面积; μ是流体动力学粘性系数; v是流速; γ是形状系数。因此, 对于含有复杂孔隙连通结构的真实储层岩石, 还可以使用COMITI等^[45]提出的公式计算雷诺数。

1.3 Klinkenberg修正模型

受Knudsen扩散效应的影响, 气体分子在固体表面流动时发生滑移, 引起渗透率的增加。因此, 在这种情况下, 孔隙介质的气体渗透率往往大于流体渗透率。考虑到Knudsen效应的影响, KLINKENBERG^[46]提出了气体渗透率的修正公式:

$ {\kappa _a} = {\kappa _\infty }\left( {1 + \frac{{{b_\kappa }}}{P}} \right) $

(10)

式中:κ_a是发生滑移现象时的表观渗透率; κ_∞是等效流体渗透率, 表示在压力非常大的情况下可以忽略气体边界滑移效应的渗透率; P是孔隙流体压力; b_κ是滑移系数。b_κ与平均自由程λ、孔喉半径r和孔隙流体压力P之间的关系可表示如下:

$ \frac{{{b_\kappa }}}{P} = \frac{{4c\lambda }}{r} $

(11)

1.4 Knudsen模型修正的B-K模型

BESKOK等^[37]提出了微管道中流体流动的二阶修正模型方程为:

$ q = \left[ {1 + \alpha {{\left( {{K_{\rm{n}}}} \right)}^2}} \right]\left( {1 + \frac{{4{K_{\rm{n}}}}}{{1 - b{K_{\rm{n}}}}}} \right)\frac{{{\rm{ \mathsf{ π} }}{r^4}}}{{8\mu }}\frac{{\Delta P}}{L} $

(12)

式中:q是压力作用下通过微管道的流量; r是流道半径; L是流道长度; K_n是Knudsen数; μ是流体动力学粘性系数; α是无量纲稀薄系数。α是K_n和参数b的函数, 可表示如下:

$ \alpha = \frac{{2{\alpha _0}{{\tan }^{ - 1}}\left( {{\alpha _1}{K_{\rm{n}}}^{{\beta _1}}} \right)}}{{\rm{ \mathsf{ π} }}} $

(13)

式中:α₁=4;β₁=0.4;α₀=α_Kn→∞=64/[3π(1－4/b)], 滑移流动时, b=-1, α₀=64/(15π)。CIVAN^[47]提出采用无量纲稀薄系数的Knudsen修正模型, 研究表明利用上述两种修正模型得到的结果无明显差异^[35]。

微管道中流体流动方程可表示为:

$ {q = {f_{\rm{c}}}{q_{{\rm{Darcy }}}}} $

(14)

$ {{f_{\rm{c}}} = \left( {1 + \alpha {K_{\rm{n}}}} \right)\left( {1 + \frac{{4{K_{\rm{n}}}}}{{1 - b{K_{\rm{n}}}}}} \right)} $

(15)

式中:f_c为Knudsen修正系数。

同样地, KLINKENBERG^[46]对气体渗透率κ_a的修正公式可表示为:

$ {\kappa _{\rm{a}}} = f_{\rm{c}}^{\rm{K}}{\kappa _\infty } = \left( {1 + 4c{K_{\rm{n}}}} \right){\kappa _\infty } $

(16)

式中:f_c^K为Klinkenberg修正系数, f_c^K=1+4cK_n。

1.5 Knudsen高阶修正模型

ZHU等^[38]提出了具有滑移效应的渗透率非线性方程:

$ {\kappa _{\rm{a}}} = {\kappa _\infty }\left( {1 + A{{\rm{e}}^{\frac{\alpha }{P}}}} \right) $

(17)

式中:A是与κ_∞成反比的常数; α是与Knudsen数有关的系数。对方程(17)进行泰勒展开得到^[38]:

$ {\kappa _{\rm{a}}} = {\kappa _\infty }\left\{ {1 + A\left[ {1 + \frac{\alpha }{P} + {{\left( {\frac{\alpha }{P}} \right)}^2}\frac{1}{{2!}} + \cdots } \right]} \right\} $

(18)

从(18)式可以看出, 公式(16)是其一阶近似。研究显示, 当孔隙介质固有渗透率小于1mD时, 低渗透率岩心中气体遵循公式(16), 当固有渗透率大于1mD时, 渗透率与压力倒数之间存在非线性关系, 需要利用修正模型对其进行修正^[38]。

2 测试分析 2.1 致密储层中Forchheimer模型和Knudsen修正模型使用条件

实际应用中, 我们需要判断致密储层中岩石孔隙流体处于怎样的流动状态, 是否需要对渗透率进行相应的模型修正。Forchheimer模型和Knudsen修正模型是两种重要的修正模型, 孔隙流体雷诺数和Knudsen数大小是判断模型是否需要修正的条件。

2.1.1 判断Forchheimer模型修正的雷诺数

雷诺数表示为R_e=ρdv/μ, 对于甲烷气体来说, 压力和温度对气体密度有重要影响, 因此密度计算需要考虑测试环境条件。根据理想气体状态方程, 可以得到密度表达式为:

$ \rho = \frac{{MP}}{{RT}} $

(19)

式中:M是甲烷摩尔质量, M=16.0425×10^-3kg/mol; P, T分别为压力和温度; R为理想气体常数, R=8.31441J/(mol·K)。

温度和压力对甲烷气体的粘性系数有重要影响, HEIDARYAN等^[48]通过实验测量, 给出了粘性系数与温度和压力的拟合关系式:

$ \mu = \frac{{{A_1} + {A_2}{p_{\rm{r}}} + {A_3}p_{\rm{r}}^2 + {A_4}p_{\rm{r}}^3 + {A_5}{T_{\rm{r}}} + {A_6}T_{\rm{r}}^2}}{{1 + {A_7}{p_{\rm{r}}} + {A_8}{T_{\rm{r}}} + {A_9}T_{\rm{r}}^2 + {A_{10}}T_{\rm{r}}^3}} $

(20)

式中:p_r=P/p_c, p_c为甲烷临界值, p_c=4.595MPa; T_r=T/T_c, 其中T_c为临界温度, T_c=190.55K^[49-50]。

我们计算了不同温度、压力、孔隙特征尺寸和流速条件下流体的雷诺数, 图 1为不同流道特征尺寸和压力条件下(温度300K)甲烷气体的雷诺数分布情况, 该图可以用于判断是否需要对致密储层Forchheimer模型进行修正。当岩心孔隙结构分布和压力条件已知时, 可以迅速得到岩石孔隙流体对应的雷诺数, 并根据雷诺数大小来判断是否需要进行Forchheimer模型修正。雷诺数是流体流速大小的线性函数, 因此, 在已知流速的情况下, 可以计算得到雷诺数数值。其它温度条件下, 雷诺数随压力和孔隙尺寸变化的二维分布情况可采用上述方法求得, 此处不再赘述。

2.1.2 判断是否需要滑移流动模型修正的Knudsen数

利用(1)式和(2)式, 可以计算出不同温度、压力和孔隙特征尺寸条件下的Knudsen数。基于Knudsen数的大小, 根据流体发生滑移流动、瞬变流动和自由分子流动的条件, 可以划分不同流体的流动情况。对于甲烷气体来说, 当给定岩石孔隙结构分布和测试压力时, 可以根据不同流道特征尺寸和压力条件下甲烷气体的流动情况(图 2), 判断孔隙气体是否发生了滑移流动, 并进行模型修正。

2.2 Knudsen修正模型适用性计算分析

达西定律在常规气藏中起主导作用。在致密储层中, 当孔隙和孔喉尺寸仅有几纳米时, 达西定律只考虑了粘性流体流动, 忽略了气体扩散和滑移运动。相较于常规气藏的其它储层岩体, 天然气页岩的孔隙和喉道较小, 因此, 在天然气页岩中, 分子与孔壁的碰撞比分子本身的碰撞更为频繁。在这种情况下, 气体滑移过程变得非常重要, 基质渗透率取决于流经页岩的气体类型以及气体压力, 此时必须考虑Knudsen流动的影响。在等温条件下, 因为气体的平均自由程与压力的倒数成正比, 与分子半径平方的倒数成正比((1)式), 所以压力越低, 分子越小, 平均自由程越大, Knudsen数随之增大。综上可知, 在孔隙和喉道较小的岩石中气体滑移效应更明显, 这导致了孔隙流量的增加, 进而提高了基质渗透率。

从不同流道特征尺寸下甲烷气体的稀薄系数和修正系数分布情况(常温常压条件下)(图 3)可以看出, 当流道特征半径从100nm减小至10nm时, 无量纲稀薄系数α为1.1~1.3(图 3a), 修正系数f_c为7~60(图 3b), 较大的无量纲稀薄系数带来更明显的修正效果, 此时孔隙流体处于瞬变流动区间, Knudsen流动的影响不可忽略。

由于孔隙介质中的液体平均自由程很小, 基本跟分子自身尺寸在同一量级, 因此Knudson数很小。从不同流道特征尺寸下水的稀薄系数和修正系数分布(图 4)可以看出, 当流道特征半径从100nm减小至10nm时, 修正系数在1附近变化(最大值为1.1), 此时孔隙流体处于达西流动或滑移流动中, Knudsen流动的影响基本可以忽略。

甘油分子半径约为0.435nm, 摩尔质量为92.09g/mol, 由分子平均自由程公式$l = M/(\sqrt 2 \cdot {\rm{ \mathsf{ π} }}d_{\rm{m}}^2\rho {N_{\rm{A}}})$(M是分子质量, d_m为分子直径, ρ为甘油密度(20℃下约1.2613g/cm³)^[51], N_A为阿伏伽德罗常数)和Knudsen数, 可以计算出不同孔隙尺寸下Knudsen数的分布和修正系数。由Knudsen模型修正的B-K模型公式可知, 流体稀薄系数的值反映了Knudsen模型修正的B-K模型修正项的影响(方程(12)), 因此, 从稀薄系数变化曲线可以直观判断流体是否需要进行Knudsen修正。

从图 5可以看出, 当流道特征半径从100nm减小至10nm时, 甘油稀薄系数进一步降低, 只有水的稀薄系数的1/2, 修正系数在1附近变化(最大值为1.02), Knudsen数也进一步降低至0.0004~0.0036, 此时孔隙流体处于无滑移流动状态, Knudsen流动的影响可以完全忽略。

实验研究显示, 在压力逐渐增大过程中, 不同孔隙流体的平均自由程收敛于一个接近0的极限值^[52]。LI等^[53]研究了3类致密砂岩的孔隙-孔喉结构与孔隙度、渗透率之间的关系, 第1类的致密砂岩包含较多的纳米孔, 第2类的致密砂岩纳米孔与微米孔含量基本相当, 第3类的致密砂岩含有较多微米孔。研究发现, 致密砂岩渗透率与纳米孔和微米孔的孔隙度正相关, 基于这3类致密砂岩实验数据, 可以得到不同压力下的致密砂岩渗透率和平均自由程随压力倒数的变化关系(图 6), 其中归一化渗透率和平均自由程分别表示为κ/κ_min, l/l_min, 其中κ_min和l_min是渗透率和平均自由程的最小值。表 3、表 4和表 5分别为3类致密砂岩的实验参数(孔喉半径、压力)和渗透率^[53]。

计算结果显示, 致密砂岩渗透率与压力相关, 这种关系还与岩石中孔隙-孔喉结构类型分布有关。渗透率随压力增加(压力倒数降低)逐渐降低, 并收敛至接近于0。第1类和第2类致密砂岩的渗透率与压力倒数具有一定的线性关系(图 6a), 这与LETHAM^[52]的实验数据结果一致。第2类和第3类致密砂岩渗透率变化曲线较为平缓, 第1类致密砂岩的收敛过程更加陡峭, 这表明在包含较多纳米孔隙的致密砂岩中, 压力对渗透率的影响更为显著, 其原因在于纳米孔隙通道中流体滑移效应更明显, 压力对平均自由程的影响更为显著。平均自由程计算结果显示, 平均自由程随着压力增大(压力倒数降低)而降低, 二者呈线性关系(图 6b), 同时第1类致密砂岩中平均自由程随压力的变化更为剧烈, 这表明包含更多纳米孔隙的致密砂岩中流体的滑移效应更明显。

图 7、图 8和图 9为3类致密砂岩甲烷气体的稀薄系数和修正系数分布情况, 可以看出, 第1类致密砂岩包含较多纳米孔隙, 最小流道特征半径为32nm, 与常压下的甲烷气体修正系数f_c相比(图 3b), 较高压力下甲烷气体修正系数数值较小, 无量纲稀薄系数α为0.2~0.5, 滑移修正的影响较小, Knudsen数为0.001~0.010, 依据ZIARANI等^[35]的划分标准, 孔隙流体仍然处于无滑移流动区间(图 7), 而参考ROY等^[36]提出的连续介质无滑移流动区间标准(K_n < 0.001), 孔隙流体已经发生滑移流动。基于第2类和第3类致密砂岩实验测量数据, 计算出的无量纲稀薄系数α为0.2~0.5, 依据Knudsen数也可推算出流体处于无滑移流动区间(图 8和图 9)。

数值计算和实验数据对比分析表明, 在实验条件下, 孔隙压力往往达到几十个甚至几百个大气压, 在如此高的压力下, 孔隙气体的密集性极高, 平均自由程很小, 即使是在包含较多纳米孔隙的储层岩石中, Knudsen数通常也小于0.010。由此可知, Knudsen模型适用于研究致密储层中的非稠密气体(低压)流动, 不适用于研究高压气体和液体(油、水等)流动占主导地位的情况。

2.3 Forchheimer模型适用性计算分析

计算水在不同流道特征尺寸介质中流动时的雷诺数(扩散速度不太大), 结果如图 10所示, 水密度为1000kg/m³, 粘性为0.001Pa·s, 流动速度变化范围是0.1~1.0m/s, 孔隙特征尺寸变化范围是20~100nm。

计算微纳米孔隙介质中甲烷气体在不同流道特征尺寸介质中的雷诺数(扩散速度不太大), 结果如图 11所示, 在温度300K和压力20MPa的情况下, 甲烷气体密度约为128.7kg/m³, 粘性为3.3×10^－5Pa·s, 流动速度变化范围是0.1~1.0m/s, 孔隙特征尺寸变化范围是20~100nm。

计算车用机油SAE 10(20℃)在不同流道特征尺寸介质中的雷诺数的分布情况(扩散速度不太大), 结果如图 12所示, 油密度为775kg/m³, 粘性为0.065Pa·s, 流动速度变化范围是0.1~1.0m/s, 孔隙特征尺寸变化范围是20~100nm。

从图 10、图 11和图 12可以看出, 在微纳米孔隙介质中, 在扩散速度不太大的情况下, 流动的水和油的雷诺数低于0.1, 远未达到Forchheimer流动状态, 仍旧处于粘性流体的层流流动状态。甲烷气体的雷诺数接近1, 也仍未达到典型的Forchheimer流动状态。由此可知, 在流动速度低于1m/s时, Forchheimer修正模型不适用于研究致密储层中的流体流动。

3 结论

本文分析了致密储层中微纳米孔隙的复杂流体模型及其适用范围。对不同孔隙结构和流体参数的模型研究发现, Knudsen修正模型适用于研究微纳米孔隙中较低气压下气体的流动修正, 对于致密储层中油、水等液体或致密气体来说, Knudsen流动的影响可以忽略; Forchheimer模型适用于研究孔隙介质流速较高的情况, 当流动速度小于1m/s时, 雷诺数较小, 流动惯性项引起的修正可以忽略。模型对比分析表明, Knudsen修正模型是流体不满足连续介质流体力学特征时使用的修正模型, 通常用于孔隙特征尺寸很小的致密储层, Forchheimer模型是在流体湍流现象严重时使用的修正模型, 其流体流动仍然符合连续介质流体力学特征。对于孔隙特征尺寸很小的低孔渗储层, 孔隙流体的流速通常很小, 湍流现象不常见, 此类储层中Forchheimer模型的修正效应不明显。