基于LB方法平板流立方定律的数值模拟研究

朱益腾

摘要：格子Boltzmann方法模拟岩体裂隙渗流场，较传统的数值模拟具有天生的并行特性、边界条件处理简洁、程序易于实施、能得到清晰的物理图像等独特优势。通过JRC数值生成法建立起在一定范围内的节理粗糙度系数JRC的二维岩体吻合裂隙结构。基于格子Boltzmann方法（LBM），设置出入口为非平衡态外推格式、上下岩体裂隙表面为标准反弹格式的边界条件，建立了模拟岩体裂隙渗流模型。通过C++编程实现格子Boltzmann方法模拟并验证了经典的泊肃叶流动。压力差为0.002和隙宽比为400：21的粗糙裂隙内，岩体裂隙渗流流态是线性次立方水流，并且随节理裂隙粗糙度系数JRC值的增加，次立方值偏离数值3越明显。在JRC=14.1的情况下，裂隙渗流在不同压力差驱动和不同隙宽、不同压力差下，可表现为次立方渗流、立方渗流和超立方渗流。

Abstract： Lattice Boltzmann Method （LBM） has been used to simulate the seepage field of the fractured rock masss， which is born with the advantages of parallel characteristics， simple boundary conditions， easy to implement the program， clear physical image compared with the traditional numerical simulation. JRC numerical generation method is used to establish the fracture structure of two-dimensional rock mass with JRC of the joint roughness coefficient in a certain range. Based on the Lattice Boltzmann Method （LBM）， the the boundary conditions are that the inlet and outlet are set to the non-equilibrium extrapolation scheme， and the crack surface of the upper and lower rock mass is set to the standard rebound format. C++ programming can realize the Lattice Boltzmann Method simulation and verify the classic Poiseuille flow. In the rough fracture with a pressure difference of 0.002 and the width ratio of 400：21， the flow pattern of the rock mass seepage is linear sub-cubic flow， and with the increase of the JRC value， the deviation of the cubic value from the value of 3 is more obvious. In the case of JRC=14.1， the the fracture seepage can be represented by the sub-cubic seepage， cubic seepage and super-cubic seepage under different pressure difference and different fracture width.

关键词：格子Boltzmann方法（LBM）；数值模拟；吻合岩体单裂隙渗流；立方定律

Key words： Lattice Boltzmann Method （LBM）；numerical simulation；water flow in the single fractured rock mass with a close correlation；cubic law

中圖分类号：TU45 文献标识码：A 文章编号：1006-4311（2017）10-0094-04

0 引言

岩体是指在一定范围内的自然地质体，它经历了漫长的自然历史过程，经受了各种地质作用，并在地应力作用下在其内部形成了各种地质构造行迹，如节理、褶皱、劈理、层理等，正是由于这些岩体裂隙的存在，为地下水的存储和流动提供了场地，由此导致的地质灾害及渗流问题在裂隙型油气田的开发、深基坑和隧道的开挖等方面得到重视并逐渐变成一个热点问题。因此，研究岩体渗流问题具有非常大的工程意义和工程应用价值。

从20世纪70年代开始W.Wittke用有限元法计算裂隙中的水流、Krizek等用有限差分法计算了各种裂隙网络系统内的流势形态。相应的，已有的计算流体力学方法也可以分为微观方法、介观方法和宏观方法三类[1]。介观方法中，流体被视为离散成一系列的流体粒子。最常见的介观模拟方法是格子Boltzmann方法。本文基于格子Boltzmann方法，在平直裂隙面下验证了该方法的可靠性，同时还分析研究了岩体节理裂隙粗糙度系数JRC、隙宽、压力差与渗流特性间的相互联系。

1 立方定律基本理论

由于理想的平板在自然界岩体裂隙中是不存在的，岩体裂隙往往是粗糙凹凸不平的，所以在光滑平板下推得的立方定理在粗糙裂隙渗流中不完全适用，所以要在原来的公式中做一些修正。依照上式我们已经知道隙宽的指数n=3时称之立方定理，为了让n=3命名为立方相应，本文将n<3命名为次立方，n>3命名为超立方。前苏联学者Lomize[11]、Amadei及速宝玉[12]等用试验研究，认为处于层流状态的天然裂隙渗流依然满足立方定理，仅需要对系数进行修正。目前，在实际工程应用和科学研究中，使用该定理的比较多。

2 格子Boltzmann模型

格子Boltzmann方法（Lattice Boltzmann method）诞生至今已有28年，格子Boltzmann方法以介观层次为视角，建立将流体离散成流体粒子、物理区域离散成一系列的格子、时间离散成一系列的时步的模型。一个完整的格子Boltzmann模型一般由三部分组成：格子，即离散速度模型；平衡态分布函数；分布函数的演化方程。

2.1 格子Boltzmann方法的理论

本文选用的D2Q9模型中的离散速度如图2。

3 讨论

3.1 二维平板模型中的泊肃叶流的模拟验证

为了验证该模型的可靠性，本文用D2Q9做了一个简单的平板模型：长为200格子单位，宽为50格子单位。进口压力P-in=1.001；出口压力P-out=0.999；雷诺数Re=20；tau=1.2；niu=0.125。

为了更为直观地看到平板流中流场的信息，图4给出了计算收敛的数值解成像后的流场信息图。

为了使取值更加精确，本文取平板中间部位（x=100）上的横截面水平速度值和Poiseuille理论值进行对比，如图5所示，可以很明显地看到两者吻合度非常高。

3.2 压力差、裂隙开度与n值的关系

从图6可以看出当压力差值ΔP=0.31、0.41时，刚开始当隙宽为11时n值均大于3为超立方流，随着隙宽的增加n值逐渐减小由立方流转变为次立方渗流；当压力差ΔP=0.21、0.11时，n值随隙宽先增大后逐渐减小，但n值始终未超过3，所以均为次立方渗流；当压力差ΔP=0.01、0.001时，n值随隙宽的增加逐渐增加，最后逐渐趋于平缓，并且n值始终未超过3，渗流表现为次立方渗流特性。

4 结论

本文使用C++编程技术，数值模拟生成了吻合的岩体单裂隙粗糙面。在LB方法的基础上，采用D2Q9模型，模拟了平板流并且研究了在压力作用驱动下的单裂隙的渗流特性，得到了以下结论：

①LB方法模拟泊肃叶流数值解和理论解误差非常小，可以很好地模拟泊肃叶流。

②利用LB方法来模拟平板流是可靠的，LB方法适用于本文所要用来模拟的岩体裂隙渗流模型。

③可以将裂隙开度和压力差值作为次立方渗流和超立方渗流的判别标准，即当压力差值ΔP≥0.31且隙宽小于等于31时裂隙渗流为线性超立方渗流，当压力差值ΔP≥0.31且隙宽大于31时裂隙渗流为线性立方渗流或次立方渗流；当压力差值ΔP<0.31时裂隙渗流为线性次立方水流。

参考文献：

[1]王志良，申林方，李邵军.岩体裂隙面粗糙度对其渗流特性的影响研究[J].岩土工程学报，2016：1-7.

[2]盛金昌，速宝玉.裂隙岩体渗流应力耦合研究综述[J].岩土力学，1998，19（2）：92-98.

[3]LIU J， ELSWORTH D， BRADYB H. Linking stress-dependent effective porosity and hydraulic conductivity fields to RMR[J]. International Journal of Rock Mechaincs and Mining Sciences，1999，36（5）：581-589.

[4]AMERICAN， BARTON N. Review of a new shear strength criterion for rock joints[J]. Engineering Geology， 1973， 7（4）：287-332.

[5]李化，黄润秋.岩石结构面粗糙度系数JRC定量确定方法研究[J].岩石力学与工程学报，2014，33（2）：3489-3497.

[6]Hyun-Sic Jang， Seong-Seung Kang， Bo-An Jang. Determination of Joint Roughness Coefficients Using Roughness Parameters[J].Rock Mech Rock Eng （2014） 47：2061-2073.

[7]Tse R，Cruden DM.Estimating joint roughnesscoefficients[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts . 1979， 16（2）：303-307.

[8]Maerz NH， Franklin JA， Bennett CP.Joint roughness measurement using shadow profilometry[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts . 1990，27（5）： 329-343.

[9]Clarke，Keith C.Computation of the FRACTAL DIMENSION OF TOPOGRAPHIC SURFACES USING THE TRIANGULAR PRISM SURFACE AREA METHOD [J]. Computers and Geosciences. 1986（12）：713-722.

[10]Turk N，Greig M.J，Dearman W.R，Amin F，F. Characterization of rock jointsurfaces by fractal dimension[C]//The28th U.S. Symposium on Rock Mechanics （USRMS）. Rotterdam：A.A. Balkema，1987.

[11]Lomize G M. Flow in Fractured Rock （in Russian ）[M].Gosemergaizdat， Moscaw， 1951： 127-129.

[12]速寶玉，詹关礼，赵坚.仿天然岩体裂隙渗流的实验研究[J].岩土工程学报，1995，17 （5）：19-24.

[13]Barton N， Bandis S， Bakhtar K. Strength Deformation and Conductivity Coupling of Rock Joints[J]. Int. J. RockMech. Min.Sci. &Geomech. Abstr， 1985， 22 （3）：121.