基于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.