基于Lattice Boltzmann方法的页岩张性裂缝渗流特征研究

曲冠政，周德胜，彭 娇，Randy Doyle Hazlett，Siwei Wu

(1.西安石油大学，陕西 西安 710065；2.The University of Tulsa McDougall School of Petroleum Engineering，Tulsa，74104，USA)

0 引 言

页岩储层的开发是油气开发领域的热点问题[1-5]，通常采用水平井分簇压裂技术使压裂裂缝沟通诱导裂缝和天然裂缝，实现流体从微纳米级基质到井筒的有效渗流。裂缝渗流特征及空间分布决定了储层的渗流特征，明确裂缝中流体渗流特征是研究页岩储层渗流的基础。平行板模型忽略粗糙性影响，首先被用于描述裂缝结构[6-9]，随后意识到粗糙性对渗流的重要影响，Lomize、李士斌等[10-15]从各个角度修正或描述裂缝的粗糙性，但获取的均是裂缝整体渗流特征，对裂缝内渗流物理场变化及特征参数分布等的认识尚不完善。目前关于裂缝内渗流特征的描述，可见报道的只有鞠杨等[10]采用人工手段计算渗流速度变化，但其裂缝物理模型是基于W-M函数生成的，本质上并非三维裂缝结构。根据岩石的破裂机理，页岩储层裂缝可归纳为张性裂缝和剪切滑移裂缝[16]。因此，采用物理实验和数值模拟技术相结合的方式，明确张性裂缝结构对流体渗流的影响，为研究单裂缝结构及裂缝网络结构渗流特征提供参考。

1 张性裂缝结构的获取

裂缝结构的获取和描述是研究流体在裂缝中渗流的关键，裂缝的获取方法主要有巴西劈裂实验法和函数生成法。函数生成法中，以W-M函数的应用最为广泛，其优势在于构造二维裂缝，但无法构造三维裂缝结构。为保证研究的可参考性，拟采用巴西劈裂实验获取裂缝结构。采用巴西劈裂实验将取自美国Barnett盆地的页岩岩心劈裂成人工裂缝，并用三维轮廓仪收集裂缝面高度分布数据；在所获裂缝面上随机提取表观尺寸为2.6 mm×1.5 mm的区域用于研究分析，并在后续研究中将其划分为离散渗流单元体。

2 流体在张性裂缝流动的模拟方法

2.1 Lattice Boltzmann模型

Lattice Boltzmann方法将流体宏观运动离散为微观尺度上粒子的迁移和碰撞，时间和空间的离散通过粒子离散速度联系；在保证微观粒子运动特征参数密度、动量和能量等的统计结果与宏观特征参数保持一致的前提下，通过微观粒子间的碰撞和迁移运动来推演宏观流体运动[17-21]。在单相流的Lattice Boltzmann模型中，碰撞项采用BKG模型，离散速度模型中Qian等[17]提出的DdQm模型(d为空间维数，m为离散速度数)被广泛采用。因此，流体在三维裂缝中的流动采用LBM-BKG方法D3Q19模型进行模拟。

2.2 准确性验证

采用Lattice Boltzmann方法专业软件PowerFLOW进行流体渗流的模拟计算，该软件是目前唯一一款基于Lattice Boltzmann方法开发的计算流体力学的商业软件。陈耀松、Wang、Du和曲冠政等的研究成果证明了该专业软件计算的准确性和可靠性[22-25]。

2.3 模拟参数

取水作为模拟流体，温度为293.15 K，运动黏度为1.007×10-6m2/s，流体密度为998.203 kg/m3；出口压力设为大气压；格子单位下，运动黏度通过调整松弛时间来确定，模拟中运动黏度取1/30，格子单位下流体密度取0.22，格子速度为1/3，由雷诺数计算公式，可得初始格子速度为1/30。流动模拟方面设置入口流速和出口压力边界条件。不考虑流体压缩性的影响。模拟裂缝体表观尺寸为2.6 mm×1.5 mm，模拟中固定裂缝开度为0.1 mm。采用笛卡尔正交网格，为确保模拟精度，裂缝开度方向上分辨率取10。将纵向和裂缝面上分别设为光滑边界条件和标准反弹边界条件。

3 张性裂缝流体渗流LBM模拟结果

3.1 渗流参数

在具体表征粗糙裂缝结构中流体渗流时的主要参数包括衡量惯性力与黏性力的雷诺数，反映压力和惯性力的欧拉数、泊肃叶数(雷诺数的函数)等[15]。以上参数均是基于流体特征参数(密度、黏度、裂缝开度)及渗流场参数(压力和速度)计算出来的流动分析参数。由于流体密度、黏度及裂缝开度等均为定值，参数的变化主要取决于压力和速度的分布，因此，研究集中考察压力和速度的变化。

3.2 可视化渗流场

3.2.1 整体渗流特征

以渗流速度为0.4 m/s为例，对流体在粗糙裂缝内的渗流进行数值模拟实验研究。由数值模拟实验可知：当渗流速度为0.4 m/s时，整体上裂缝中流体压力沿渗流方向呈线性下降趋势，粗糙性并不影响压力的线性趋势；同一渗流纵向截面上渗流速度呈非均匀分布，呈现两端小中间大的趋势；裂缝面处渗流速度基本为0；粗糙性造成不同截面渗流速度差别很大，以截面中心处最大速度的差距尤为明显；沿渗流流线上各处速度值并不相等，因流线的迂曲性而有差异。

3.2.2 三维方向上速度分布

为观察渗流速度分布，将速度场三维剖分。在垂直于x方向选取的10条速度切片，同一切片上速度呈非均匀分布，若沿y方向上裂缝面粗糙度较小则速度分布相对均匀；若局部粗糙度较大则明显形成渗流高速区和低速区，其原理与多孔介质渗流中流体沿大孔道突进的原理一致；在垂直于y方向均匀选取6条速度切片，速度分布形态与垂直于x方向选取的速度切片上的速度分布情况基本一致，粗糙性造成速度分布的非均匀性；在垂直于z方向的速度切片上也较明显展示了速度分布的非均质性。综合分析认为：裂缝中渗流稳定后，尽管渗流场宏观速度为平均速度，但由于粗糙性影响，速度分布依据空间位置的不同而相异；粗糙度大的区域，渗流速度较小；粗糙度较小的区域，渗流速度较大。研究局部粗糙度较大的裂缝面处的流体渗流发现：裂缝面的粗糙性造成了流体渗流路径的迂曲性；在裂缝面处，速度矢量受粗糙性影响明显，在部分裂缝面形态显著变化的区域表现出渗流速度的杂乱无序性，产生回流区；回流区的产生增加了流体渗流的阻力。

3.2.3 渗流速度等值面分布

考察渗流速度为0.4 m/s时裂缝中渗流速度等值面的分布：在平均渗流速度附近(0.2～0.5 m/s)时，该区间流线间是相互平行的，只有部分粗糙性变化程度较剧烈区域由于出现回流导致渗流速度的非一致性；当渗流速度偏离平均速度程度较大时，粗糙性对渗流速度等值面分布的影响较大， 渗流速度小于0.1 m/s的区域基本对应于粗糙性最大的部分，而渗流速度大于0.6 m/s的区域基本对应于粗糙性变化较为平缓的区域。

3.3 截线上渗流特征

在渗流场中随机选取渗流截线，分析截线上的压力和速度分布情况。随机在z方向上高度为1.43 mm处的渗流平面上沿渗流方向和垂直渗流方向上任意选取3条水平截线(沿渗流方向上3条随机渗流截线y分别为0.55、0.15、0.10 mm；垂直于渗流方向上的3条随机截线x分别为0.75、1.70、2.00 mm)，分析其压力和速度分布情况可知：压力沿渗流方向总体呈斜率固定的线性下降趋势且截线的斜率保持一致，但压力曲线稍有波动并非呈纯粹的线性下降趋势；垂直于渗流方向，压力分布随坐标变化呈无规则波动，该模拟中，沿垂向压力波动在1 Pa以内，约占渗流方向压力变化的5%；沿渗流方向和垂直于渗流方向上，同一水平截线上渗流速度差别很大；粗糙性导致裂缝面高度分布的非均匀性，流线依托于裂缝面形态，裂缝形态导致渗流流线的迂曲性。

3.4 张性裂缝离散单元渗流特征

为明确张性裂缝粗糙性结构中不同区域渗流特征参数分布，将裂缝面沿渗流方向均分为5个区域，中心位置分别为0.15、0.45、0.75、1.05、1.35 mm，并用垂直于xy平面尺寸为0.2 mm×0.3 mm的65个长方形结构覆盖渗流场，考察离散单元内的渗流特征。图1为离散单元内渗流速度与位置的关系。由图1可知：图1a、b中最大渗流速度分别为0.45、2.30 m/s，分别位于2.70、0.15 mm处，最小渗流速度分别为0.31、1.54 m/s，分别位于2.70、1.35 mm处，渗流速度最大值与最小值的差距分别为31%、33%；同一渗流区域中，沿渗流方向各点的速度相异；不同渗流区域，同一渗流截面上不同位置的速度亦不同。因此，由于粗糙性影响渗流速度因空间位置而异，在垂直于渗流剖面上同样有压差存在，前面关于截面压力分布的研究也证明了该观点；模拟渗流速度为0.4、2.0 m/s时，速度分布趋势基本一致。因此，在分析中可选取某一渗流速度下的渗流特征，从差异辨别度的角度考虑选取模拟渗流速度为2.0 m/s进行渗流特征分析。

图1 微观渗流速度分布

图2为渗流区域各渗流速度矢量值。由图2可知：x方向渗流速度矢量值与渗流总速度矢量值基本接近，而y方向渗流速度矢量值最小，基本为0；z方向渗流速度矢量值在0附近波动，波动幅度较大，如中心位置y=0.75 mm处渗流速度分布图。分析显示：整个渗流场沿y方向压力略有波动，但y方向上压差的驱动作用不足以抵消粗糙性导致的渗流场迂曲性的阻力作用，因此，y方向上渗流速度值较小；模拟渗流方向沿x方向，模拟速度为2.0 m/s，同时压力梯度沿x方向保持恒定值，流体压能及动能完全能克服粗糙性导致的迂曲性阻力作用及裂缝面的摩擦阻力作用等；同时，由于渗流路径沿z方向的迂曲性必然会造成z方向上具有一定的渗流速度分量，该分量值与驱动压差、渗流速度及裂缝面形态等相关。

图2 渗流速度矢量图

常用的表征裂缝面粗糙性参数包括粗糙度、迂曲度、分形维数[26-30]，粗糙度和分形维数是从裂缝面高度分布的复杂程度考虑，迂曲度是从渗流路径考虑。考虑到数据点采集量，采用粗糙度定性描述裂缝面高度分布[29]。

中心位置分别为0.15、0.45、0.75、1.05、1.35 mm，所对应的沿渗流方向各区域粗糙度分别为0.038、0.056、0.055、0.037、0.077 mm。5个横向渗流区域渗流特征见图3。为考察垂直于渗流方向上的渗流特征，纵向划分为13个相互平行的渗流区域，每个区域包含5个渗流单元，渗流特征见图3b、d、f。由图3a、b可知：沿渗流方向各区域渗流速度各异，分布为1.90～2.10 m/s，按照渗流速度由大至小对5个区域进行排序，所对应区域的中心位置依次为0.15、1.05、0.45、0.75、1.35 mm，5个区域粗糙度依次为0.038、0.037、0.056、0.055、0.077 mm。渗流速度基本与粗糙度呈反相关关系，个别区域稍有出入。分析表明：离散单元渗流速度具有非一致性，相邻区域间由于速度差异而存在剪切力，同时裂缝面附近的回流区域与粗糙度及渗流路径变化率有关，回流区对渗流速度有阻力作用；垂直于渗流方向各截面渗流速度也不一致，在1.96～2.02 m/s内波动，波动范围比沿渗流方向波动程度小(图3b)。

由图3c、d可知：横向截面压力总体呈线性均匀分布；相对于纵向压力分布，横向压力梯度很小；沿渗流方向压差总体呈线性下降趋势，但相邻考察点间压差呈波动趋势，波动趋势与局部粗糙性有关；总之，2个方向压力总体均呈线性分布，渗流阻力与渗流路径迂曲程度有关，因此，压力梯度随粗糙度的增加而增加。由图3e、f可知：x方向渗流速度比渗流总速度略小，横向截面上的压力梯度造成横向上具有较小波动；y、z方向渗流速度分量基本在0附近；图3f中，沿渗流方向上，x方向渗流速度依然占据渗流速度绝对比例，各剖面y方向渗流速度基本在0附近，由于渗流路径的迂曲性造成z方向速度分量值为-0.50～0.50 m/s。

图3 渗流场内部区域渗流特征

4 结 论

(1) 整体上压力沿渗流方向呈线性分布，但受粗糙性影响压力下降趋势具有小幅波动性；受裂缝面粗糙性影响，在垂直于渗流方向上也存在呈线性分布的压降，但相比于渗流方向的压降很小，可忽略。

(2) 沿渗流方向任意截线上，压力均呈波动状线性分布且线性分布规律保持一致；垂直于渗流方向截线上，压力分布受裂缝面粗糙性影响具有随机性。

(3) 沿渗流方向和垂直于渗流方向上切片各处的渗流速度不同；同一渗流流线上不同位置处渗流速度略有差别；受裂缝面粗糙性变化影响，部分区域会出现速度的杂乱无序性分布，产生回流区；沿渗流方向和垂直于渗流方向上，任选一条平行于xy平面的截线，其各处速度值具有随机性，该分布特点与粗糙性造成的流线分布特点有关。

(4) 当速度值在平均速度附近时，渗流速度等值面间基本上相互平行；当渗流速度偏离平均速度值较大时，速度等值面分布与裂缝面形态密切相关，渗流速度较低的等值面基本分布于裂缝面形态较粗糙位置，渗流速度较高的等值面基本分布于裂缝面形态分布平缓的区域。

(5) 渗流单元体中的渗流速度各不相同。由于粗糙性造成渗流路径迂曲性的差异，渗流单元体中的渗流总速度并不相等。渗流路径迂曲性造成三维方向上均有渗流速度，但沿渗流方向速度分量占总速度分量的绝对比重。

[1] 蒋裕强，董大忠，漆麟，等.页岩气储层的基本特征及其评价[J].天然气工业，2010，30(10)：7-12.

[2] 侯冰，陈勉，李志猛，等.页岩储集层水力裂缝网络扩展规模评价方法[J].石油勘探与开发，2014，41(6)：763-768.

[3] 张士诚，郭天魁，周彤，等.天然页岩压裂裂缝扩展机理试验[J].石油学报，2014，35(3)：496-503.

[4] 吴奇，胥云，张守良，等.非常规油气藏体积改造技术核心理论与优化设计关键[J].石油学报，2014，35(4)：706-714.

[5] 程远方，李友志，时贤，等.页岩气体积压裂缝网模型分析及应用[J].天然气工业，2013，33(9)：53-59.

[6] WITHERSPOON P A，WANG J S Y，IWAI K，et al.Validity of cubic law for fluid flow in a deformable rock fracture[J].Water Resources Research，1979：16(6)：1016-1024.

[7] 温庆志，刘欣佳，黄波，等.水力压裂可视裂缝模拟系统的研制与应用[J].特种油气藏，2016，23(2)：136-139，158.

[8] 温庆志，胡蓝霄，翟恒立，等.滑溜水压裂裂缝内砂堤形成规律[J].特种油气藏，2013，20(3)：137-139，158.

[9] 李小龙，肖雯，王凯，等.支撑剂在清洁压裂液中的沉降规律[J].大庆石油地质与开发，2015，34(2)：95-98.

[10] LOMIZE G M.Water flow through jointed rock [M].Moscow：Gesenergoizdat，1951：122.

[11] STEPHEN R Brown.Transport of fluid and electric current through a single fracture[J].Journal of Geophysical Research，1989，94(B7)：9429-9438.

[12] ERDINC Eker，SERHAT Akin.Lattice Boltzmann simulation of fluid flow in synthetic fractures[J].Transport in Porous Media，2006，65(3)：363-384.

[13] 张程宾，陈永平，施明恒，等.端面粗糙度的分形特征及其对微通道内层流流动的影响[J].物理学报，2009，58(10)：7050-7056.

[14] 李士斌，陈波涛，张海军，等.清水压裂自支撑裂缝面闭合残留宽度数值模拟[J].石油学报，2010，31(4)：680-683.

[15] 鞠杨，张钦刚，杨永明，等.岩体粗糙单裂隙流体渗流机制的实验研究[J].中国科学 E辑：技术科学，2013，43(10)：1144-1154.

[16] 钟建华，刘圣鑫，马寅生，等.页岩宏观破裂模式与微观破裂机理[J].石油勘探与开发，2015，42(2)：242-249.

[17] QIAN Y H.Lattice BGK models for Navier-Stokes Equation [J].Europhysics Letters，1992，17(6)：479-484.

[18] 何雅玲，王勇，李庆，等 .格子Boltzmann方法的理论与应用[M].北京：科学出版社，2009：98-100.

[19] 郭照立，郑楚光.格子Boltzmann方法的原理及应用[M].北京：科学出版社，2009：121-123.

[20] ZHAO Liguo，Zhao T S.Lattice boltzmann model for incompressible flows through porous media[J].Physical review E，2002，66(3)：162-173.

[21] YOSHINO M，INAMURO T.Lattice boltzmann simulation for flow and heat/mass transfer problems in a three-dimensional porous structure[J].International journal for numberical methods in fluids，2003，43(2)：183-198.

[22] 陈耀松，单肖文，陈沪东，等.计算流体力学的新方向及其在工业上的应用[J].中国科学 E辑：技术科学，2007，37(9)：1107-1116.

[23] WANG Yiwei，WANG Yang，AN Yiran，et al.LBM calculation of high-speed train aerodynamics[J].Science China Technological Sciences，2008，38(11)：1795-1804.

[24] DU Tezhuan，LI Xiangqun，ZHANG Xialiang，et al.Lattice Boltzmann Method used for the aircraft characteristics computation at high angle of attack[J].Science China Technological Sciences，2010，53(8)：2068-2073.

[25] 曲冠政，曲占庆，HAZLETT R D，等.页岩拉张型微裂缝几何特征描述及渗透率计算[J].石油勘探与开发，2016，43(1)：115-120.

[26] BELEM T，HOMAND-ETIENNE F.Quantitative parameters for rock joint surface roughness[J].Rock Mechanics and Rock Engineering，2000，33(4)：217-242.

[27] ZHANG Zhenyu，JAN Nemcik，QIAO Qiuqiu，et al.A model for water flow through rock fractures based on friction factor [J].Rock Mechanics and Rock Engineering，2015，48(2)：559-571.

[28] 闫寒，张文明，胡开明，等.随机粗糙微通道内流动特性研究[J].物理学报，2013，62(17)：174701-174706.

[29] LAMAS L N.An experimental and analytical study of the roughness of granite joints[J].Hydraulic Engineering，2013，33(4)：240-248.

[30] MAKURAT A，GUITERREZ M.Fracture flow and fracture cross flow experiments[C].SPE36732，1996：1-9.