0 引言

近年来，随着常规油气资源储量的减少以及油气勘探开发技术的提高，我国正持续加大以页岩气、致密气为代表的非常规油气资源的勘探开发力度^[1-3]。不同于常规油气藏，页岩气藏富含微纳米级孔隙，孔渗条件差，基质孔隙度通常 < 12%，渗透率通常 < 1 mD^[4]。目前，页岩气藏的商业开发主要依赖于通过大型水力压裂来尽可能沟通储层中发育的天然微裂缝。所以，开展气体在页岩储层中的渗流机理研究，对于页岩气藏的渗流理论研究及工程实践都具有重要意义。

由于页岩气藏独特的储层物性，页岩气的运移过程呈现出从分子尺度到宏观尺度的多尺度的特性。页岩气的流动机理包括了自由气在人工裂缝中的达西流动和在基质微孔隙中的连续流、滑脱流以及过渡流，吸附气在微孔隙表面吸附/解吸以及扩散^[5]。通常对非常规裂缝性气藏渗流规律的研究手段主要有室内实验以及数值模拟方法^[6-7]。由于非常规气藏超低孔、低渗的储层特性，对实验设备要求较高，流体的渗流过程也难以进行直观地观测，所以，通过数值方法模拟流体在多孔介质中的渗流方式得到广泛的应用。目前已有大量的用于表征气体在纳米级孔隙中运移的模型。此类模型主要分为2类：一类是基于Hagen-Poiseuille模型对其边界条件进行修正，从而表征页岩气的多种运移机制^[8-10]；另一类是按一定的权重系数对多种运移机制进行耦合并考虑气体分子间的相互作用建立的真实气体传输模型^[11]。此类模型仅能用于描述页岩气在单一纳米孔隙中的运移规律，尚不能直接用于研究页岩气在天然微裂缝发育的复杂多孔介质中的运移特征。

格子Boltzmann方法(LBM)以其易于处理复杂几何边界问题、编程较简单等优势在近年来被大量用于非常规气藏微观渗流模拟研究^[12-14]。LBM模拟分为孔隙尺度以及表征单元(REV)尺度，其中孔隙尺度主要用于模拟流体在微纳米尺度孔隙中的流动，而不能直接用于模拟流体在天然裂缝发育的多孔介质中的流动过程^[15]。Guo等^[16]通过在广义格子Boltzmann方程中添加修正项来实现孔隙尺度向REV尺度的升级。此后，该方法被大量运用于非常规油气藏渗流模拟。Zhao等^[17]运用REV尺度的格子Boltzmann方法，模拟了气相在重构多孔介质中的渗流过程，并对多孔介质渗透率进行了计算。Zhao等^[15]应用格子Boltzmann方法对页岩气在基质中的渗流过程进行了REV尺度的模拟，并对页岩气渗流规律的影响因素进行了分析。张烈辉等^[18]基于页岩扫描电镜图像运用REV尺度的格子Boltzmann方法模拟了页岩气在单一天然微裂缝中的渗流。Zhang等^[19]基于页岩扫描电镜建立的二值图像模拟了水相在富含有机质页岩多孔介质中的渗流过程，并通过在二值图像上规则排布天然裂缝的方式分析了天然裂缝的存在对页岩气藏渗流规律的影响。

本文首先重构天然微裂缝发育的页岩储层多孔介质模型，然后对Zhao等^[15]以及Liu等^[20]提出的格子Boltzmann模型围绕页岩气藏特性对演化方程以及边界条件进行针对性修正，并将其升级到REV尺度，使其能够表征页岩气在裂缝性多孔介质中滑脱、吸附/解吸以及气体扩散等运移机制；应用修正模型模拟气体在天然微裂缝发育的多孔介质中的流动过程；最后基于模拟结果对气体渗流特性以及渗透率影响因素进行分析，以期为页岩气藏的渗流规律研究提供新思路。

1 页岩气藏表征单元体尺度微观渗流模型 1.1 物理模型

本文建立的页岩气藏渗流模拟流动域如图 1所示。整个渗流过程考虑二维流动，流场上、下均为不渗透边界；整个渗流过程为从入口处通过压力驱动从左至右的流动；气体渗流过程的速度分布用格子Boltzmann方法中的二维九速模型(D2 Q9)进行表征。由于模型的求解过程复杂，为了提高求解效率，对模型作出了如下假设：①流场中仅考虑甲烷气体的单相流动；②渗流过程为等温渗流过程；③对于天然微裂缝的物性参数仅考虑裂缝长度、裂缝宽度以及裂缝密度。

1.2 修正REV尺度的格子Boltzmann模型 1.2.1 外力项修正

本文使用LBM来表征气体在上述物理模型中的渗流过程。该方法是将Boltzmann方程进行离散，把流体分子运动分为迁移和碰撞等2个过程^[21]。广义格子Boltzmann演化方程为

$ {f_i}(x + {e_i}\Delta t,t + \Delta t) - {f_i}(x,t) = - \frac{1}{\tau }[{f_i}(x,t) - f_i^{eq}(x,t)] + \Delta t{F_i}(x,t),\quad i = 0,1, \cdots ,N $

(1)

式中：x为格子空间格点位置；t为格子空间的某一时刻；f_i(x，t)为t时刻格子空间中x格点沿i方向的气体密度分布函数；Δt为格子空间的时间步长；τ为驰豫时间；F_i(x，t)为格子空间x格点沿i方向的外力项函数；e_i为格子空间的离散速度。

e_i在D2 Q9模型中气体分子的运动方向分为以下9个方向。

$ {e_i} = \left\{ {\begin{array}{*{20}{l}} {(0,0),}&{i = 1}\\ {\left[ {{\rm{cos}}\frac{{(i - 1)\pi }}{2},{\rm{sin}}\frac{{(i - 1)\pi }}{2}} \right],}&{i = 2 \backsim 5}\\ {\left[ {{\rm{cos}}\frac{{(2i - 9)\pi }}{4},{\rm{sin}}\frac{{(2i - 9)\pi }}{4}} \right],}&{i = 6 \backsim 9} \end{array}} \right. $

(2)

式中：i为方向编号。

在式(1)中，f_i^eq为平衡态分布函数，本文模型考虑多孔介质的非均质性，其表达式如下^[22]：

$ f_i^{eq}(x,t) = {\omega _i}\rho \left[ {1 + \frac{{u{e_i}}}{{c_{\rm{s}}^2}} + \frac{{{u^2}\left( {e_i^2 - c_{\rm{s}}^2I} \right)}}{{2\varphi c_{\rm{s}}^4}}} \right] $

(3)

式中：ω_i为加权系数，取值为ω₀ = 4/9，ω_{1 - 4} = 1/9，ω_{5 - 8} = 1/36；c_s为格子空间声速，${c_{\rm{s}}} = \frac{{\sqrt 3 c}}{3}$，c一般取值1，故${c_{\rm{s}}} = \frac{{\sqrt 3 }}{3}$；I为单位矩阵；ρ为格子空间流体密度；u为格子空间流体速度；φ为多孔介质孔隙度，%。

驰豫时间对于LBM模型至关重要，由于在页岩气藏中，气体的流动特性由努森数(Kn)控制，所以对于驰豫时间必须进行相应的修正，其表达式^[13]为

$ {\tau _{\rm{e}}} = \frac{1}{2} + \sqrt {\frac{6}{\pi }} \frac{L}{{\Delta {\kern 1pt} x}}\frac{{Kn{\kern 1pt} \varphi (Kn){{(1 + 0.5b{\kern 1pt} \rho {\kern 1pt} \chi )}^2}}}{\chi } $

(4)

式中：L为流场特征长度，m；Δx格子尺寸，m；φ(Kn)为考虑努森层的修正函数，φ(Kn) = $\frac{1}{{1 + 2Kn}}$；χ为修正系数。

χ的计算式为

$ \begin{array}{l} \chi = 1 + \frac{5}{8}b\rho + 0.286{\kern 1pt} {\kern 1pt} {\kern 1pt} 9{(b\rho )^2} + 0.110{\kern 1pt} {\kern 1pt} {\kern 1pt} 3{(b\rho )^3} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.038{\kern 1pt} {\kern 1pt} {\kern 1pt} 6{(b{\kern 1pt} \rho )^4} \end{array} $

(5)

式中：b = 2 π d³/(3 m)。其中，d为甲烷分子直径，m；m为甲烷分子质量，kg。

气体的动力黏度也是影响渗流过程的重要参数，其在格子空间的计算式为

$ \nu = c_{\rm{s}}^2\left( {{\tau _{\rm{e}}} - \frac{1}{2}} \right)\Delta {\kern 1pt} t $

(6)

式(1)中的格子空间外力项F_i考虑孔隙表面对于气体分子的作用，其离散形式为

$ \begin{array}{*{20}{l}} {{F_i} = {\omega _i}{\kern 1pt} \rho \left( {1 - \frac{1}{{2{\tau _{\rm{e}}}}}} \right)\left[ {\frac{{{e_i}{\kern 1pt} F}}{{c_{\rm{s}}^2}} + \frac{{u{\kern 1pt} F(e_i^2 - c_{\rm{s}}^2)}}{{\varphi {\kern 1pt} c_{\rm{s}}^4}}} \right] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} G{\kern 1pt} \varPhi (x,t)\sum\limits_i {{\omega _i}} (x + {e_i}\Delta t,t){e_i}} \end{array} $

(7)

式中：Φ(x, t)为势函数；F为格子空间气体分子所受总外力项。

F的计算式为

$ F = - \frac{{\varphi {\kern 1pt} \nu }}{K}u - \frac{{\varphi {\kern 1pt} {F_\varepsilon }}}{{\sqrt K }}|u|u + \varphi {\kern 1pt} G $

(8)

式中：G为格子空间气体分子所受外部力，主要为基质及裂缝表面对气体分子的吸附力；F_ε为几何形状因子；K为固有渗透率，m²。

F_ε与K的计算式^[23]为

$ {F_\varepsilon } = \frac{{1.75}}{{\sqrt {150{\kern 1pt} {\varphi ^3}} }},K = \frac{{{\varphi ^3}d_{\rm{p}}^2}}{{150{{(1 - \varphi )}^2}}} $

(9)

式中：d_p为平均孔隙直径，m。

气体的宏观密度以及速度可以通过下式进行计算

$ {\rho = \sum\limits_{i = 0}^8 {{f_i}} (x,t)} $

(10)

$ {\rho {\kern 1pt} u = \sum\limits_{i = 0}^8 {{e_i}} {f_i}(x,t) + \frac{{\Delta {\kern 1pt} t}}{2}\rho {\kern 1pt} F} $

(11)

由于式(8)中也包含流体流速u，所以要对式(8)进行求解，必须联立式(8)及式(11)

$ u = \frac{V}{{{d_0} + \sqrt {d_0^2 + {d_1}|V|} }} $

(12)

V为临时变量，其计算式为

$ V = \frac{{\sum\limits_{i = 0}^8 {{e_i}} {f_i}(x,t)}}{\rho } + \frac{{\Delta {\kern 1pt} t}}{2}\varphi {\kern 1pt} G $

(13)

联立式(12)及式(13)则可以得到气体宏观速度u的计算式为

$ u = \frac{{\frac{{\sum\limits_{i = 0}^8 {{e_i}} {f_i}(x,t)}}{\rho } + \frac{{\Delta {\kern 1pt} t}}{2}\varphi G}}{{{d_0} + \sqrt {d_0^2 + {d_1}\left| {\frac{{\sum\limits_{i = 0}^8 {{e_i}} {f_i}(x,t)}}{\rho } + \frac{{\Delta {\kern 1pt} t}}{2}\varphi {\kern 1pt} G} \right|} }} $

(14)

d₀及d₁的完整表达式为

$ {d_0} = \frac{1}{2}\left( {1 + \frac{{\Delta {\kern 1pt} t}}{2}\frac{{\varphi {\kern 1pt} v}}{K}} \right),{d_1} = \frac{{\Delta {\kern 1pt} t}}{2}\frac{{\varphi {\kern 1pt} {F_\varepsilon }}}{{\sqrt K }} $

(15)

1.2.2 边界条件修正

本文模型考虑了页岩气藏中气体吸附/解吸及扩散效应，所以必须对传统的格子Boltzmann方法的边界条件进行修正。岩石表面的气体滑移流动速度的表达式^[24]为

$ {u_{\rm{s}}} = (1 - \alpha ){u_{\rm{g}}} + \alpha {\kern 1pt} {u_{\rm{w}}} $

(16)

式中：u_g为离边界层最近的气流速度，m/s；u_w为假定的固体表面速度，m/s；α为比例系数。

表征的气体在孔隙及裂缝表面的吸附情况，可由下式进行计算

$ \alpha = \frac{{{C_{\rm{ \mathsf{ μ} }}}}}{{{C_{{\rm{ \mathsf{ μ} }}{\kern 1pt} {\rm{max}}}}}} = \frac{{\frac{P}{{{P_{\rm{L}}}}}}}{{1 + \frac{P}{{{P_{\rm{L}}}}}}} $

(17)

式中：C_μ及C_μmax分别为吸附气浓度及最大吸附浓度，g/m³；P_L为兰格缪尔压力，MPa。

式(16)中的u_w为假定的边界滑移速度，在页岩气藏须要考虑气体的吸附以及表面扩散作用，其表达式为

$ {u_{\rm{w}}} = \frac{{{D_{\rm{S}}}}}{{{D_{\rm{K}}}{\rho _{\rm{a}}}}}\left[ {{\rho _{\rm{g}}}{u_{\rm{g}}}\frac{{k{\rho _{{\rm{ \mathsf{ μ} }}{\kern 1pt} {\rm{max}}}}}}{{{{(1 + k\rho )}^2}}}} \right] $

(18)

式中：D_K为气体分子扩散系数，m²/s；D_S为气体表面扩散系数，m²/s。

在格子空间s可由下式进行计算^[25]：

$ {D_{\rm{S}}} = C_{\rm{S}}^2({\tau _{\rm{e}}} - 0.5)\Delta {\kern 1pt} t $

(19)

2 结果与讨论 2.1 模型验证

Guo等^[23]提出了对REV尺度的格子Boltzmann模型的验证方法，后来该方法被大量应用于REV尺度格子Boltzmann模型准确性的验证^{[15, 19-20]}。本文也选取该方法对所建模型的准确性进行验证。

首先设置模型上、下不渗透平板间的多孔介质孔隙度为20%，流体在模型入口处受到恒定的压力驱动，其流动过程满足Brinkman达西方程

$ \frac{{{\nu _{\rm{e}}}}}{\varphi }\frac{{{\partial ^2}u}}{{\partial {\kern 1pt} {y^2}}} - \frac{\nu }{K}u + G = 0 $

(20)

其解析解为

$ u = \frac{{G{\kern 1pt} K}}{\nu }\left\{ {1 - \frac{{{\rm{cosh}}\left[ {\frac{{a(y - H)}}{2}} \right]}}{{{\rm{cosh}}\left( {\frac{{a{\kern 1pt} H}}{2}} \right)}}} \right\} $

(21)

式中：$a = \sqrt {\nu \varphi /K{\nu _{\rm{e}}}} $；ν，ν_e分别是流体的黏度和有效黏度，本文取ν/ν_e = 1，m²/s；L，H分别为流场的长度和宽度，m。

格子空间中的一些必要参数取值为：网格数L_x × L_y = 100 × 100，流体密度ρ = 1，黏度ν = 1/6，驱动压差ΔP = 1。图 2为不同达西数下的解析解与本文所建立模型的数值解的对比。从图 2可看出，本文模型的计算结果与经典解析解吻合，证明了本文所建立的REV尺度格子Boltzmann模型(LB模型)具有较高的准确性。

2.2 算例及参数敏感性分析 2.2.1 基本算例

本文模型模拟所用的基本参数及其格子空间尺度向真实物理空间尺度的转换关系如表 1所列。在LBM中主要通过保持格子空间和物理空间流动特征数不变的准则来实现尺度的转换。将表 1中的参数带入本文所建立的LB模型中，然后通过编程求解可以得到气体在重构页岩基质多孔介质中的渗流动态。图 3展示了在不同时间步长下的气体流速分布情况。从图 3可以看出，当天然裂缝较发育时，大量气体沿天然裂缝窜进，天然裂缝沟通程度越显著时气流速度越大。以上结论证实了天然裂缝发育及储层中天然裂缝的沟通作用对于页岩气藏开发十分重要。

2.2.2 天然微裂缝物性参数分析

为了进一步明确天然微裂缝对于页岩气藏气体渗流特性的影响，图 4展示了不同裂缝密度下(ρ_f= 0.02，0.04，0.08)的气体流速分布情况。在3种模拟条件下，仅改变裂缝的线密度、裂缝宽度W_f =(2~ 4)×10^-6 m，裂缝长度l_f =(1~4)×10^-4 m、裂缝倾角θ = 0°~60°，其余流体参数保持表 1中的数值不变。从图 4可以看出，当裂缝密度增大时，气体流动通道增加，渗流能力显著增强。

图 5为不同天然裂缝密度条件下，所建立裂缝性多孔介质的表观渗透率与基质固有渗透率比值的变化情况。从图 5可以看出，随着天然裂缝密度的增加，表观渗透率将显著增大。

同样仅改变裂缝开度，保持天然裂缝ρ_f = 0.04及l_f=(1~4)×10^-4 m不变，将裂缝参数以及流体参数带入LB模型，计算所得表观渗透率的增大情况如图 6所示。从图 6可以看出，随着裂缝开度的增加，表观渗透率增量会随之增大。这是由于天然裂缝开度增加后，气体在天然裂缝的流动通道增大，从而引起气体表观渗透率增大。

仅改变裂缝长度，保持其他裂缝参数为ρ_f= 0.04及W_f =(2~4)×10^-6 m不变，同样将上述裂缝参数及流体物性参数带入本文LB模型进行计算，得到表观渗透率增大倍数情况如图 7所示。从图 7可以看出，随着裂缝长度的增加，表观渗透率的增大倍数也随之增加。对比裂缝宽度对于表观渗透率增大量的贡献，裂缝长度对于表观渗透率的增大作用更为显著。这是由于随着裂缝长度的增加，基质向裂缝的窜流会加剧，从而引起表观渗透率的大幅上升。

2.2.3 气体渗流特征参数分析

在页岩气藏中，气体运移规律不同于常规储层。为了进一步明确页岩储层中气体多重运移的规律，有必要对气体渗流特征参数进行敏感性分析。

当储层致密，气体渗流通道狭窄时，气体分子自由程将接近于其流动通道的特征长度，此时Kn对于气体流动的影响将加剧。气体在页岩储层中的流动可由Kn的大小将其划分为4种流态^[26]。其中Kn ≤ 0.001为连续介质流；0.001 < Kn ≤ 0.100为滑脱流；0.1 < Kn ≤ 10.0为过渡流；Kn > 10为自由分子流。为了能够定量分析Kn对于气体流动特征的影响，笔者模拟了Kn为0.007，0.069，0.687条件下的气体渗流过程。

图 8表明，随着Kn的增大，气体流动过程中克氏效应越显著，气体流速也随之增大。为了能够直观观测Kn对于气体渗流速度的影响，图 9展示了不同Kn下的气体流速剖面。从图 9可看出，当Kn持续增大后，裂缝壁面的滑脱效应更加显著，使得气体沿天然微裂缝窜进的速度大幅提升；特别是当Kn > 0.1时，流态进入过渡流后，气体流速增幅显著变大。图 10为不同Kn下的裂缝性多孔介质表观渗透率随出口压力的变化情况。从图 10可以看出，随着Kn的增大，表观渗透率也随之增大；特别是当Kn > 0.1时，表观渗透率增幅明显增大。

页岩气藏中气体的吸附/解吸对于气体渗流有巨大的影响。为了能够直观表征气体吸附/解吸作用对于页岩气藏表观渗透率的影响，本文计算了不同吸附气量下多孔介质的表观渗透率。图 11为不同吸附气浓度下表观渗透率的变化情况。从图 11可以看出，随着吸附气量的增加，多孔介质的表观渗透率随之下降。出现该结果的主要原因是由于在页岩有机质中，吸附气吸附在孔隙表面会占据部分游离气的流动通道，使基质的有效孔隙减小，从而影响气体的渗流^[27]。

目前，大量研究表明^[28-30]，页岩气藏有机质中吸附气的表面扩散作用是影响页岩气渗流过程的关键因素。在本文模型中，气体表面扩散作用的强弱是通过表面扩散系数D_s来进行表征。图 12为不同表面扩散系数比值下(D_s/D_k = 0，0.01，0.1，1)，页岩气藏多孔介质表观渗透率随压力的变化情况。从图 12可以看出，在压力较低的情况下，考虑气体表面扩散效应的表观渗透率数值较不考虑表面扩散效应的情况能够提升2~5倍。这样的模拟结果表明气体表面扩散效应对于表观渗透率的贡献作用是不可忽略的。

在页岩有机质中，气体表面扩散效应并不能无限变化，其强弱受制于吸附气量的多少。所以，表面扩散系数以及吸附气量是影响页岩气表面扩散效应的2个关键参数。图 13为考虑不同表面扩散效应下，多孔介质表观渗透率随吸附气量的变化情况。从图 13可以看出，表面扩散效应对表观渗透率的贡献会随着吸附气量的增大而增加，吸附气量越大时，不同表面扩散系数对表观渗透率的影响也愈发显著。

3 结论

(1) 重构了天然微裂缝发育的页岩储层多孔介质，考虑页岩气渗流过程中的滑脱、吸附/解吸、表面扩散效应等运移机制对广义格子Boltzmann模型进行了修正，并基于经典REV尺度格子Boltzmann模型验证方法对重构模型的正确性进行了验证。

(2) 应用重构模型对页岩气渗流过程进行了模拟。结果表明，当储层天然微裂缝较发育时，绝大部分气体沿天然微裂缝窜进，天然微裂缝成为页岩气渗流的主要通道。通过对天然微裂缝物性参数进行敏感性分析发现，裂缝密度对于页岩气渗流的促进作用最为显著；随着裂缝密度增大3~4倍，储层表观渗透率可增大10倍以上；裂缝长度及裂缝开度对表观渗透率的贡献则较小。

(3) 参数敏感性分析表明，努森数(Kn)是影响页岩气渗流特征的主要参数，随着Kn的增大，克氏效应显著，气体流速随之增大，特别当Kn > 0.1时，表观渗透率增幅明显增大；表观渗透率会随着吸附气量的增大而减小，特别当储层压力较低时，该现象更为显著；表面扩散效应也是影响页岩气渗流的主要因素之一，同等条件下考虑气体表面扩散效应的表观渗透率较忽略该因素可提升2~5倍，但提升作用受制于吸附气量的多少。