引言

孔隙网络模型是预测多孔介质流动特性的有效手段，通常采用球体表征多孔介质的孔隙，等截面柱体表示喉道，并将孔隙空间拓扑结构抽象为孔隙中轴线，直观体现孔隙空间和喉道的连通情况和分布情况^[1-3]。

孔隙网络模型由Fatt^[4-6]于1956年提出，Fatt通过随机设置网格连线的半径建立了简单的二维随机网络模型，并预测了毛管压力和相对渗透率。1977年，Chatzis和Dullien^[7]在Fatt模型的基础上建立了三维随机网络模型，其研究结果表明三维模型的空间表征能力更强。1990年，Jerauld和Salter^[8]通过建立规则孔隙网络模型，研究了孔隙空间拓扑结构以及孔隙、喉道的空间关联性。2005和2006年，王克文等^[9-10]分别用Weibull方法建立了随机孔隙网络模型，并分析了孔喉结构参数及润湿性对电阻率的影响。2006年，陈民锋等^[11]结合逾渗理论，采用Weibull截断分布作为孔喉特征分布函数，建立了三维孔隙网络模型并分析了微观孔隙结构参数（如孔喉比、配位数、形状因子）对水驱油特征的影响。2007年，姚军等^[12]通过控制孔喉尺寸分布和配位数，建立了三维随机网络模型，并提高了模型的连通性。侯健等^[13]于2011年在随机孔隙网络模型的基础上分析了微观参数对宏观参数的影响。1992年，Bryant等^[14-16]1992年，首次进行了真实孔隙网络模型的研究。随后，Oren和Bakke^[17-19]通过模拟岩石的沉积、压实和成岩的过程建立了基于过程法的孔隙网络模型。而随着X射线CT扫描技术的应用，得以通过算法提取出真实岩芯的孔隙网络模型^[20-22]。随机孔隙网络模型相对真实孔隙网络模型对岩芯围观孔喉参数反映差，但基于真实岩芯的微观孔喉参数建立的随机孔隙网络模型也可反映出岩芯的孔喉特征。

本文利用基于真实页岩气藏岩芯的微观孔喉参数建立了随机孔隙网络模型，并模拟了大于临界流速的不同流动速率下渗流路径。

1 流动控制方程

流体在毛管束中流动流量与流动压降成正比

$ {q_{ij}} = {g_{ij}}\Delta {p_{ij}} $

(1)

式中：${q_{ij}}$—节点$i$和$j$之间的流量，cm$^3$/s；

${g_{ij}}$—节点$i$和$j$之间的导流系数，cm$^3$/(s$\cdot$MPa)；

$\Delta {p_{ij}}$—节点$i$和$j$之间的压差，MPa。

页岩气储层的孔隙尺度为微纳米级，需考虑“滑脱效应”，Javadpour等人^[23]给出了一个理论的滑脱因子$F$来修正管流中的滑脱速度。

$ F= 1{\rm{ + }}{\left( {{{10}^3} \times \dfrac{{8{\rm{R}}T}}{{\mathsf{π} M}}} \right)^{0.5}}\dfrac{\mu }{{ p {r_{ij}}}}\left( {\dfrac{2}{\alpha } - 1} \right) $

(2)

式中：$F$—无因次系数；

R—气体常数，R=8.314 J/(mol$\cdot$K)；

$T$—温度，K；

$M$—摩尔质量，kg/kmol；

$\mu$—气体黏度，Pa$\cdot$s；

$p$—平均压力，MPa；

${{r_{ij}}}$—节点$i$和$j$之间的喉道平均半径，cm；

$\alpha$—切向动量调节系数，0.8。

气体通过纳米孔隙的总质量流量是Knudsen扩散和压力差作用的结果，Knudsen扩散系数${D_{\rm{K}}}$的定义如下^[24-25]

$ {D_{\rm{K}}} = \dfrac{{2{r_{ij}}}}{3} {\dfrac{{8{\rm{R}}T}}{{\mathsf{π} M}}} $

(3)

因此，导流系数满足方程^[23-24]

$ {g_{ij}} = \left[ {\dfrac{{2{r_{ij}}M}}{{3000{\rm{R}}T{\rho _{{\rm{avg}}}}}}{{\left( { \dfrac{{8000{\rm{R}}T}}{{\mathsf{π} M}}} \right)}^{0.5}} + F\dfrac{{r_{ij}^2}}{{8\mu }}} \right]\dfrac{{\mathsf{π} r_{ij}^2}}{{{l_{ij}}}} + \\\hspace{5em}{D_{\rm{K}}}{c_{\rm{g}}}{\rho _{{\rm{avg}}}} $

(4)

式中：$\rho _{\rm{avg}}$—气体平均密度，kg/m$^3$；

$l_{ij}$—节点$i$和$j$之间的喉道长度，cm；

$c _{\rm{g}}$—气体压缩系数，Pa$^{-1}$。

由质量守恒方程可知流经任一节点的流出流量等于流入流量，即总流量为0

$ \sum\limits_{i, j} {{q_{ij}}} = 0 $

(5)

而流量与流动压降成正比，则流动条件为

$ {q_{ij}} = \left\{ \begin{array}{l} {g_{ij}}({p_i} - {p_j} + {p_{\rm{c}}}), {p_i} - {p_j} + {p_{\rm{c}}} > 0\\ 0, {p_i} - {p_j} + {p_{\rm{c}}} < 0 \end{array} \right. $

(6)

式中：

$p_i$—节点$i$处的压力，MPa；

$p_j$—节点$j$处的压力，MPa；

$p_{\rm{c}}$—节点$i$和$j$之间的毛管力，MPa。

毛管力$p_{\rm{c}}$可根据Young-Laplace方程得到

$ {p_{\rm{c}}} = \dfrac{{2\sigma \cos\theta }}{r} $

(7)

式中：

$\sigma$—界面张力，N/m；

$\theta $—流体接触角，(°)(强水湿系统的接触角为180°)；

$r$—喉道半径，cm。

两相界面处，如图 1中节点5的流动方程

$ \left( {{g_{35}} + {g_{45}} + {g_{56}} + {g_{57}}} \right){p_5} - {g_{35}}{p_3} - {g_{45}}{p_4} -\\\hspace{3em} {g_{56}}{p_6} - {g_{57}}{p_7} = {C_{\rm{s}}}{g_{56}}{p_{{\rm{c}}56}} + {C_{\rm{s}}}{g_{57}}{p_{{\rm{c}}57}} $

(8)

式中：

$g_{mn}$—节点$m$和$n$($m$和$n$为节点代号，如式(8)的3、4、5、6、7等)之间的导流系数，cm$^3$/(s$\cdot$MPa)；

$p_m$—节点$m$处的压力，MPa；

$p_{{\rm{c}}mn}$—节点$m$和$n$之间的毛管压力，MPa；

$C_{\rm{s}}$—与入口处压力$p_{\rm{inlet}}$和喉道长度$l$相关的函数，表达式为

$ {C_{\rm{s}}} = \dfrac{{2\sigma }}{{{p_{{\rm{inlet}}}}l}} $

(9)

式中：

$p_{{\rm{inlet}}}$—入口处的压力，MPa；

$l$—喉道长度，cm。

然后结合Young-Laplace方程及压降$\Delta p$即可算出各节点的流量，根据达西公式则可算出岩芯的渗透率$K$

$ K = \dfrac{{q\mu l}}{{A\Delta p}} $

(10)

式中：$K$—渗透率，mD；

$q$—流量，cm$^3$/s；

$A$—岩芯横截面积，cm$^2$；

$\Delta p$—两端的的压差，MPa。

2 孔隙网络建模

本文使用的孔隙网络模型为基于逾渗理论并根据真实岩芯数据建立的随机网络模型，其非均质性主要通过连通概率$P$对节点配位数$z$的影响，以及孔喉半径的分布差异实现，即根据给定的孔喉半径尺寸和变异系数对每个结点的赋值，建模流程如图 2。

2.1 确定节点间距和模型大小

节点间距的赋值即通过计算真实岩芯的喉道长度获得。建立孔隙网络模型大小的原则为越大越好，但孔隙网络规模越大需要的计算机计算能力越强。与数字岩芯中的代表体积元(REV)原理类似，孔隙网络模型的规模达到一定尺度后其微观结构统计参数即保持稳定，因此，只需超过此阈值，建立的随机网络即可具有与真实网络相同的孔喉特征。而Jerauld等人的研究表明，当规则网络的节点数达到一定时(逾渗网络节点数需大于20$\times$20$\times$20)，随机网络模型即能反映真实岩芯的孔隙拓扑结构^{[8, 28-31]}。

2.2 确定晶格类型和配位数赋值

本文使用配位数$z_{\max}$=8的体心立方(BCC)模型。通过设定连通概率$P$对喉道半径随机赋值，即对$P$的键(喉道)赋予半径分布范围内的随机半径值，其余键的半径值则为0。其中连通概率$P$等于实际连通配位数$z$/$z_{\max}$，表达式如下

$ z = {z_{\max }}P $

(11)

式中：

$z_{\max }$—模型最大配位数，无因次；

$z$—模型连通配位数，无因次；

$P$—连通概率，%。

2.3 喉道半径赋值

考虑储层具有强非均值性，用储层岩石孔隙半径变异系数来描述其非均质性，本文网络模型采用随机函数来配置喉道半径，来实现非均质性孔隙网络模拟。喉道半径生成公式为

$ r\!=\!{{\rm{e}}^{{\rm{rand}}()\% \dfrac{{{\mathop{\rm int}} (50\lg {r_{\max }})\!-\!{\mathop{\rm int}} (50\lg {r_{\min }})}}{{50}}}} +\\\hspace{3em} \lg {r_{\min }} $

(12)

式中：

e—欧拉数，常数；

rand()%—随机函数；

$r_{\max }$—最大喉道半径，µm；

$r_{\min }$—最小喉道半径，µm。

孔隙网络模型中任一喉道半径都在式(12)中的$r_{\max }$到$r_{\min }$之间。

3 流动模拟结果及分析

本文模拟的页岩气藏储层岩石的原始数据为：孔喉半径平均值约为0.35 µm，喉道长度均值为4.63 µm，配位数均值为4.36，孔隙度为3.11%。

根据上一节提到的最低节点数要求，本文采用的模型节点数规模为70$\times$30。连通概率$P$和实际配位数$z$的对应取值如表 1所示，因此根据气藏的原始配位数4.36，本文使用的概率$P$值为54.50。

常用的喉道半径分布函数形式通常为两种：对数均匀分布和均匀分布。对数均匀分布使用的变异系数一般为0.05、0.30、0.55、0.80和1.05，均匀分布使用的变异系数一般为0.05、0.30和0.55。本文所采用的半径值为0.35 µm，半径$r$与变异系数$C_{\rm{V}}$的对应关系满足表 2；对数均匀分布和均匀分布的变异系数、$r_{\max }$和$r_{\min }$取值见表 3。

本文完成了随机孔隙网络建模和水驱气流动模拟，建立的孔隙网络模型的孔隙度为2.96%，实测孔隙度为3.11%；误差为4.82%，证明建立的孔隙网络模型是可靠的。

根据临界渗流理论，若岩石孔隙喉道空间的非均质性较强，流体倾向于在阻力较小的大孔喉通道(临界路径)中流动。水驱气过程中，较高流动速度(大于临界流速$v_{\rm{c}}$时)会破坏气/水界面稳定性，造成水侵现象，同时渗流通道水侵后可能形成的绕流、卡断和水锁现象，水相滞留则会造成气相运移通道的堵塞，造成节点渗流通道的改变，进而在网络模型中对渗流路径造成改变。当渗流速度小于临界流速时，水气两相存在明显相界面，相界面随时间推移向前推移，形成一种拟稳定的相界面的现象。

临界流速$v_{\rm{c}}$的影响因素较多，包含孔喉半径分布、喉道长度、毛管力梯度和渗透率梯度等，通常由拟合流速与渗透率曲线获得，在本文生成的孔隙网络模型中水驱气的临界流速约为0.74$\times$10$^{-6}$ mL/s。

3.1 模拟结果

在70$\times$30的随机BCC网络模型中模拟渗流速度大于临界渗流速度$v_{\rm{c}}$时水相驱替气相过程。渗流速度分别为1.1$\times$10$^{-4}$，1.1$\times$10$^{-5}$，5.5$\times$10$^{-6}$及1.1$\times$10$^{-6}$ mL/s；孔喉半径变异系数值为0.8，最大孔喉半径0.665 µm，最小孔喉半径0.035 µm，孔喉半径按照对数均布分布；连通概率为54.5%，对应配位数$z$为4.36。图 3分别对应以上不同流速的渗流情况。

在驱替速度大于临界流速的情况下，随着渗流速度的增大，水的波及系数总体呈降低的趋势，但流动路径的变化并不统一：(1) 图 3d与图 3c相比，随着流速增加，位于图中上部的优势流通通道依然存在，而部分中下部的初始未流通的孔隙在速度增加的情况下被打开。(2)在图 3c与图 3b的对比中可发现：随着流速的继续增加，上部低流速未流通的孔隙贯通成为了优势流通通道，下部通道的优势则在弱化。(3)而通过图 3b与图 3a的对比可发现，路径随压力的变化趋势发生了变化，压力的进一步增大让中部的孔喉得以重新打开并形成优势通道，上部和下部的通道则均在收缩。

对比图 3中4个流动结果图可发现，流动路径随渗流速度的变化并非一成不变：在渗流速度大于临界流速时，流速的增加总体上会强化优势流通通道，使得渗流路径变少，但同时可能会打开低渗流速度下未流通的路径，因此优势流动通道的发展趋势会发生变化，在不同的流速增加状态下会形成不同的优势流动通道。

渗流速度小于临界流速$v_{\rm{c}}$的水驱气过程中，主要影响因素为毛管力(压降梯度的影响较小)，气水两相流体的界面会趋于稳定，呈活塞式推进，波及效率比大于临界流速时高很多。

3.2 渗流路径的分形特征

水相驱替气相过程中，流速较高时，气水两相流体的界面会失稳，有利于驱替水相穿透原饱和流体，建立一组高速的渗流通道，非湿相流体抑制湿相流体的产出，渗流主要发生在两条优势渗流路径，路径与分形曲线相似(以$v$=1.1$\times$10$^{-4}$ mL/s为例，如图 4)，本文通过建立分形图形进行了对比验证，如下。

(1) 生成分形图形

分形图形常用迭代函数系统(Iterated Function System，IFS)生成，核心为利用整体和局部之间的自相似结构，通过函数迭代变换和压缩产生分形图形。其迭代过程如下：

将$\omega_1$，$\omega_2$，…，$\omega_N$个压缩变换同时作用于原始图形$E_0$上，产生$\omega_1$($E_0$)，$\omega_2$($E_0$)，…，$\omega_N$($E_0$)，生成输出图形$E_1$，表示为

$ {E_1} = \boldsymbol{W}\left( {{E_0}} \right) = {\omega _1}\left( {{E_0}} \right)\&{\omega _2}\left( {{E_0}} \right)\& \cdots {\omega _N}\left( {{E_0}} \right) $

(13)

式中：

$E_0$—原始图形，无因次；

$E_1$—第一次变换后生成的图形，无因次；

$\omega _1$—第一个压缩变换，无因次；

$\omega _2$—第二个压缩变换，无因次；

$\omega _N$—第$N$个压缩变换，无因次；

$\boldsymbol{W}$—对原始图形分别进行$\omega _1$，$\omega _2$，…，$\omega _N$个压缩变换后形成的集合，无因次；

&—逻辑“与”。

然后进行第二次迭代，将$E_1$作为输入图形，并再次将以上变换同时作用于$E_1$之上并得到$E_2$

$ {E_2}\!=\!\boldsymbol{W}\left( {{E_1}} \right)\!=\!{\omega _1}\left( {{E_1}} \right)\&{\omega _2}\left( {{E_1}} \right)\& \cdots {\omega _N}\left( {{E_1}} \right) $

(14)

式中：

$E_2$—第二次变换后生成的图形，无因次。

重复以上步骤直到收敛或生成目标图形

$ {E_{k + 1}} = \boldsymbol{W}\left( {{E_k}} \right) ~~~~~k = 0, 1, 2, \cdots $

(15)

式中：

$E_k$—第$k$次变换后生成的图形，无因次；

$E_{k + 1}$—第$k+1$次变换后生成的图形，无因次。

本文将树形渗流路径上的大枝作为原始图形$E_0$，利用上述迭代过程生成了树形分形曲线(图 4，图 5)。

(2) 自相似性质验证

分形结构具有任意小的比例细节，传统几何学难以对整体和局部的关系进行表征，同时局部和整体存在自相似。可通过迭代的方式生成分形曲线对比验证图形是否具有分形特征。本文利用Matlab程序生成分形曲线，分析渗流路径和树形分形曲线的关系，树形分形曲线的迭代过程如下图。

树形分形曲线与网络模型渗流路径比较(以$v$=1.1$\times$10$^{-4}$ mL/s和$v$=5.5$\times$10$^{-6}$ mL/s为例)如图 6，图 7。分形相似包括大尺度上的相似与细节上的相似，如图 5，图 6和图 7所示，非稳定界面的渗流路径与树形分形曲线总体吻合度较好，符合大尺度上的相似性；渗流速度正相关于压降梯度，而随着压降梯度的增大，树形路线曲线小枝始终位于压降梯度方向上的分形曲线双节点上，说明渗流路径和分形曲线具有细节上的相似性；因此，渗流路径符合分形的特征，具有自相似性，即渗流路径图为分形图形。而流动路径分形曲线的小枝总位于压降梯度方向上的双节点上，也说明了增大流速是造成渗流界面失稳的主要因素。

4 结论

（1）在毛管束流动中，渗流速度与驱替压差呈正相关，高流速的驱替过程中，压降梯度对渗流路径起主要影响作用，高渗流速度会破坏气/水界面稳定性，更利于水侵而降低产气量；而低渗流速度下的气/水界面趋于稳定，对气藏的动用程度更高。

（2）在渗流速度高于临界流速的情况下，流动路径随渗流速度的变化总体呈现出渗流速度越高流动路径越少的趋势，优势流动通道会随着流速的增加得以加强，但并非一成不变，流速的增加可能会将低速下未流通的孔隙打开，进而形成优势流动通道，因此导致优势流动通道的位置发生变化。渗流速度小于临界流速时，气水两相流体的界面会趋于稳定，呈活塞式推进，波及效率比大于临界流速时高。

（3）流体渗流速度高时，两相渗流界面失稳的渗流路径呈树形多支状，通过分形曲线自相似性验证和迭代过程模拟发现失稳渗流路径具分形特征，流动路径分形曲线为压降梯度方向上的双节点，也说明了增大流速是造成渗流界面失稳的主要因素。