﻿ 孔隙网络模拟渗流速度对页岩气开采的影响
 西南石油大学学报(自然科学版)  2019, Vol. 41 Issue (6): 19-27

Influence of Seepage Velocity on Shale Gas Exploration by Pore Network Simulation
GUO Xiao , YANG Kairui, YANG Yurui, JIA Haowei
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China
Abstract: Pore network model is an effective method to predict the flow characteristics of porous media. Based on characteristic values of shale gas reservoir, a pore network was established. And gas-water phase flow was simulated in that pore network. Seepage image was drawn by solving the pressure matrix equations of the model. By simulating flow under different velocities, the water penetration and seepage path were analyzed. As the results shows, the seepage velocity was positively correlated with the displacement pressure difference. When seepage velocity is higher than critical velocity vc, pressure drop gradient has a main effect on the seepage path, which conforms to characteristics of fractal curve, and seepage velocity has an effect on stability of oil-gas interface and destroys equilibrium of gas/water interface thus facilitates water penetration and reduce gas production when it's high. There will be a stable gas/water interface and a large amount of gas production when seepage velocity is lower than vc.
Keywords: pore network model     shale gas     flow simulation     seepage velocity     seepage path     fractal theory

1 流动控制方程

 ${q_{ij}} = {g_{ij}}\Delta {p_{ij}}$ (1)

$\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\cdotK)； T—温度，K； M—摩尔质量，kg/kmol； \mu—气体黏度，Pa\cdots； p—平均压力，MPa； {{r_{ij}}}—节点i$$j$之间的喉道平均半径，cm；

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

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

 ${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)

$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}}} = \dfrac{{2\sigma \cos\theta }}{r}$ (7)

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

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

$r$—喉道半径，cm。

 图1 流动模拟示意图 Fig. 1 Schematic diagram of flow simulation
 $\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 孔隙网络建模流程图 Fig. 2 Flow chart of pore network modeling 2.1 确定节点间距和模型大小 节点间距的赋值即通过计算真实岩芯的喉道长度获得。建立孔隙网络模型大小的原则为越大越好，但孔隙网络规模越大需要的计算机计算能力越强。与数字岩芯中的代表体积元(REV)原理类似，孔隙网络模型的规模达到一定尺度后其微观结构统计参数即保持稳定，因此，只需超过此阈值，建立的随机网络即可具有与真实网络相同的孔喉特征。而Jerauld等人的研究表明，当规则网络的节点数达到一定时(逾渗网络节点数需大于20\times20\times20)，随机网络模型即能反映真实岩芯的孔隙拓扑结构[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 流动模拟结果及分析

3.1 模拟结果

 图3 不同渗流速度下的水相流动路径 Fig. 3 Flow path of water phase at different seepage velocity

3.2 渗流路径的分形特征

 图4 失稳界面的渗流路径 Fig. 4 Flow path of water with unstable phase-interface

(1) 生成分形图形

$\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_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$次变换后生成的图形，无因次。

 图5 树形分形曲线的迭代过程 Fig. 5 Iteration process of tree curve generation

(2) 自相似性质验证

 图6 树形分形曲线与$v$=1.1$\times$10$^{-4}$ mL/s时失稳两相渗流路径比较 Fig. 6 Comparison between fractal tree curve and two-phase unstable flow path when $v$=1.1$\times$10$^{-4}$ mL/s
 图7 树形分形曲线与$v$=5.5$\times$10$^{-6}$ mL/s时失稳两相渗流路径比较 Fig. 7 Comparison between fractal tree curve and two-phase unstable flow path when $v$=5.5$\times$10$^{-6}$ mL/s
4 结论

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

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

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

 [1] SHIN H, LINDQUIST W B, SAHAGIAN D L, et al. Analysis of the vesicular structure of basalts[J]. Computers & Geosciences, 2005, 31(4): 473-487. doi: 10.1016/j.cageo.-2004.10.013 [2] JIANG Z, WU K, COUPLES G, et al. Efficient extraction of networks from three-dimensional porous media[J]. Water Resources Research, 2007, 43(12): 2578-2584. doi: 10.1029/2006WR005780 [3] CHRIS P. Distance-ordered homotopic thinning:A skeletonization algorithm for 3D digital images[J]. Computer Vision and Image Understanding, 1998, 72(3): 404-413. doi: 10.1006/cviu.1998.0680 [4] FATT I. The network model of porous media Ⅰ. Capillary characteristics[J]. Petroleum Transactions, AIME, 1956, 207: 144-159. [5] FATT I. The network model of porous media Ⅱ. Dynamic properties of a single size tube network[J]. Petroleum Transactions, AIME, 1956b, 207: 160-163. [6] FATT I. The network model of porous media Ⅲ. Dynamic properties of networks with tube radius distribution[J]. Petroleum Transactions, AIME, 1956b, 207: 164-181. [7] CHARTZIS I, DULLIEN F A L. Modeling pore structure by 2-D and 3-D networks with application to sandstones[J]. Journal of Canadian Petroleum Technology, 1997, 16: 97-108. doi: 10.2118/77-01-09 [8] JERAULD G R, SALTER S J. The effects of porestructure on hysteresis in relative permeability and capillary pressure:Pore-level modeling[J]. Transport in Porous Media, 1990, 5(2): 103-151. doi: 10.1007/BF00144600 [9] 王克文, 孙建孟, 关继腾. 两相流下岩石电性的三维网络模型模拟[J]. 测井技术, 2005, 29(4): 289-292. WANG Kewen, SUN Jianmeng, GUAN Jiteng. Three dimensional network modeling of electrical properties of rock with two phase fluids[J]. Well Logging Technology, 2005, 29(4): 289-292. doi: 10.3969/j.issn.1004-1338.-2005.04.004 [10] 王克文, 孙建孟, 关继腾, 等. 聚合物驱后微观剩余油分布的网络模型模拟[J]. 中国石油大学学报(自然科学版), 2006, 30(1): 72-76. WANG Kewen, SUN Jianmeng, GUAN Jiteng, et al. Network model modeling of microcosmic remaining oil distribution after polymer flooding[J]. Journal of China University of Petroleum, 2006, 30(1): 72-76. doi: 10.3321/j.issn:1000-5870.2006.01.015 [11] 陈民锋, 姜汉桥. 基于孔隙网络模型的微观水驱油驱替特征变化规律研究[J]. 石油天然气学报, 2006, 28(5): 91-95. CHEN Minfeng, JIANG Hanqiao. Law of characteristic variations of microscopic water displacement based on pore net work model[J]. Journal of Oil and Gas Technology, 2006, 28(5): 91-95. doi: 10.3969/j.issn.1000-9752.-2006.05.027 [12] 姚军, 陶军, 李爱芬. 利用三维随机网络模型研究油水两相流动[J]. 石油学报, 2007, 28(2): 94-97, 101. YAO Jun, TAO Jun, LI Aifen. Research on oil water two phase flow using 3D random network model[J]. Acta Petrolei Sinica, 2007, 28(2): 94-97, 101. doi: 10.3321/j.issn:0253-2697.2007.02.018 [13] 侯健, 高达, 李振泉, 等. 储层微观参数对宏观参数影响的网络模拟研究[J]. 计算力学学报, 2011, 28(1): 78-83. HOU Jian, GAO Da, LI Zhenquan, et al. Network modeling of the influence of reservoir microscopic parameters on macroscopic parameters[J]. Chinese Journal of Computational Mechanics, 2011, 28(1): 78-83. [14] BRYANT S L, KING P R, MELLOR D W. Network model evaluation of permeability and spatial correlation in areal random sphere packing[J]. Transport in Porous Media, 1993, 11(1): 53-70. doi: 10.1007/BF00614635 [15] BRYANT S L, MELLOR D W, CADE C A. Physically representative network models of transport in porous media[J]. AIChE Journal, 1993b, 39(3): 387-396. doi: 10.1002/aic.690390303 [16] BRYANT S L, BLUNT M. Prediction of relative permeability in simple porous media[J]. Physical Review A, 1992, 46(4): 2004-2012. doi: 10.1103/physreva.46.2004 [17] OREN P E, BAKKE S, ARNTZEN O J. Extending predictive capabilities to network models[C]. SPE 52052-PA, 1998. doi:10.2118/52052-PA [18] OREN P E. Process based reconstruction of sandstones and prediction of transport properties[J]. Transport in Porous Media, 2002, 46(2-3): 311-343. doi: 10.1023/A:-1015031122338 [19] ØREN P E, BAKKE S. Reconstruction of Berea sandstone and pore-scale modeling of wettability effects[J]. Journal of Petroleum Science and Engineering, 2003, 39(3-4): 177-199. doi: 10.1016/S0920-4105(03)00062-7 [20] ZHAO H Q, MACDONALD I F, KWIECIEN M J. Multiorientation scanning:A necessity in the identification of pore necks in porous media by 3-D computer reconstruction from serial section data[J]. Journal of Colloid and Interface Science, 1994, 162(2): 390-401. doi: 10.1006/jcis.-1994.1053 [21] BALDWIN C A, SEDERMAN A J, MANTLE M D, et al. Determination and characterization of the structure of a pore space from 3D volume images[J]. Journal of Colloid and Interface Science, 1996, 181(1): 79-92. doi: 10.1006/jcis.1996.0358 [22] SHEPPARD A P, SOK R M, AVERDUNK H. Improved pore network extraction methods[C]//Proceedings of International Symposium of the Society of Core Analysts, Toronto, Canada, 2005. [23] JAVADPOUR F, FISHER D, UNSWORTH M. Nanoscale gas flow in shale gas sediments[J]. Petroleum Society of Canada, 2007, 46(10): 55. doi: 10.2118/07-10-06 [24] SOEDER D J. Porosity and permeability of eastern Devonian gas shale[J]. SPE Formation Evaluation, 1988, 3: 116-124. doi: 10.2118/15213-PA [25] SAKHAEE P A, STEVEN L B. Gas permeability of Shale[C]. SPE 146944, 2011. doi:10.2118/146944-MS [26] MEHMANI A, PRODANOVIĆ M. The application of sorption hysteresis in nano-petrophysics using multiscale multiphysics network models[J]. International Journal of Coal Geology, 2014, 128-129: 96-108. doi: 10.1016/j.coal.2014.03.008 [27] MEHMANI A, PRODANOVIĆ M, JAVADPOUR F. Multiscale, multiphysics network modeling of shale matrix gas flows[J]. Transport in Porous Media, 2013, 99(2): 377-390. doi: 10.1007/s11242-013-0191-5 [28] DONG Hu, BLUNT M J. Pore-network extraction from micro-computerized-tomography images[J]. Physical Review E, 2009, 80(80): 036307. doi: 10.1103/PhysRevE.-80.036307 [29] LIU Xuefeng, SUN Jianmeng, WANG Haitao. Reconstruction of 3D digital cores using a hybrid method[J]. Applied Geophysics, 2009, 6(2): 105-112. doi: 10.1007/s11770-009-0017-y [30] JAVADPOUR F. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone)[J]. Journal of Canadian Petroleum Technology, 2009, 48(8): 16-21. doi: 10.2118/09-08-16-DA [31] YALE D P, NUR A. Network modeling of flow, storage, and deformation in porous rocks[J]. Seg Technical Program Expanded Abstracts, 1949, 4(1): 437-437. doi: 10.1190/1.1892843