魏晶晶，李連俠，廖華勝，劉德偉，李秋林

(1.水力學與山區(qū)河流開發(fā)保護國家重點實驗室，成都 610065;2.密歇根州立大學，美國 密歇根州 48824)

0 前 言

我國是世界上的筑壩大國，目前已建水庫8.7萬余座，其中土石壩占我國大壩總數(shù)的93%[1]。據(jù)權威數(shù)據(jù)統(tǒng)計，由于滲流破壞而造成土石壩失事的比例占到35%左右[2]，因此滲流是影響大壩安全的重要因素，研究滲流對于工程實踐具有重要的現(xiàn)實價值和指導作用。

傳統(tǒng)的滲流計算方法主要是實驗法，但實驗法花費大，時間長，實驗結果容易受實驗條件及周圍環(huán)境的影響，因此它的應用受到了一定的限制。隨著計算機科學的發(fā)展，數(shù)值模擬方法由于可以有效地降低成本并且較為方便的模擬復雜的流動現(xiàn)象等優(yōu)勢受到了國內(nèi)外學者的廣泛關注，并取得了豐碩的成果。鮑呈蒼[3]利用SEEP/W軟件(基于FEM)通過計算大壩浸潤線、等勢線及滲透坡降等，對滲流對大壩邊坡穩(wěn)定的影響進行了研究。張嘎等[4]采用固定網(wǎng)格有限元分析法對溪洛渡水電站上游圍堰五種防滲結構進行滲流分析，驗證了該方法可用于形狀復雜、各分區(qū)滲透性差異較大的土石壩。陳玉茹[5]采用基于二維拉格朗日有限差分法的二維有限元分析法，分析了滲流量、滲流通道對滲流路徑的影響。歐陽君等[6]將土石壩滲流的有限元計算與壩坡穩(wěn)定分析的強度折減有限元法相結合對土石壩邊坡穩(wěn)定進行分析，驗證了該方法的可行性。楊林德[7]等將飽和土體視為均質(zhì)、連續(xù)的各向異性彈塑性多孔介質(zhì)，根據(jù)虛位移原理推導出飽和土體內(nèi)各向異性滲流直接耦合的有限元法計算公式。付成華等[8]利用有限元法對土石壩進行滲流分析，得出了庫水位下降過程中按穩(wěn)定滲流和非穩(wěn)定滲流所得壩體浸潤線之間的平均相對誤差，提出了由組合參數(shù)判定非穩(wěn)定滲流壩體浸潤線簡化計算條件的3組判據(jù)。倪小東、王媛[9]將顆粒流程序運用于管涌發(fā)展階段的研究，從散粒體組成的多孔介質(zhì)的細觀力學特征出發(fā)，運用離散單元法研究了土體管涌機理，并得出管涌破壞的影響因素。Xu Y.Q.及Lam. L等[10-12]用飽和-非飽和滲流模型分析了一座土石壩斷面的滲透穩(wěn)定性。高朝軍[13]等采用邊界元通過移動網(wǎng)格法求解自由面的位置，在迭代過程中采用二次插值的方法來確定溢出點。孫偉建等[14]在虛單元法的基礎上，提出了一種滲流自由面的曲線擬合法，選取適應能力強的三角形單元進行有限元分析，無須調(diào)整自由面單元。華靜等[15]分析探討了有限元法應用于土石壩滲流計算存在的滲流自由面迭代、單元格剖分不當?shù)葐栴}并提出了相應的解決方法。

目前這些傳統(tǒng)的有限元、有限差分、邊界元等數(shù)值模擬方法發(fā)展已經(jīng)較為成熟，但由于土體復雜的孔隙結構以及流體和介質(zhì)之間復雜的相互作用使得這些方法都存在一定的問題，宏觀層面的研究只能知道滲流的統(tǒng)計平均特性，不能了解多孔介質(zhì)內(nèi)部微觀效應，這些方法也不適于處理復雜邊界條件，程序編制復雜，計算量大時會影響數(shù)值計算的穩(wěn)定性，因此，探索高效準確的數(shù)值方法對于滲流研究的進一步發(fā)展十分重要。

格子玻爾茲曼方法(Lattice Boltzmann Methods, LBM)是一種介于宏觀連續(xù)介質(zhì)模型和微觀分子動力學模型之間的介觀方法,既克服了微觀分子動力學模型所存在的時間尺度和空間尺度的局限性，同時又克服了宏觀分子動力學模型引入連續(xù)性假設所帶來的計算精度問題[16]。相比于傳統(tǒng)的CFD方法，LBM能夠模擬微觀角度多孔介質(zhì)的非連續(xù)流動，具有物理背景清晰，編程簡單，能夠處理復雜邊界條件，算法簡單且數(shù)值精度高等優(yōu)勢，并成功在流體力學、量子力學、傳熱學等諸多領域得到應用。Hazlett[17]利用LBM方法模擬研究了孔隙介質(zhì)多相流；姚軍等[18]人采用數(shù)字巖心和LBM方法分析了巖石滲透率變化；Degruyter等人[19]利用同步加速CT成像與LBM方法研究了火山浮石的孔隙結構特征并分析了孔隙率與滲透率之間的關系。格子Boltzmann方法在滲流問題的研究中也得到了成功的應用，樊火等[20]從微觀角度對裂隙滲流的機理進行研究，用MRT-LBM模擬滲流運動，提出計算滲透率的截面滲透率加權法，進而對兩種天然裂隙進行滲流模擬，并預估其滲透率。田智威等[21]采用LBM對含裂隙煤體瓦斯?jié)B流規(guī)律進行數(shù)值模擬，定量研究了滲流速度和流量隨著裂隙寬度增加的變化規(guī)律。

目前格子Boltzmann方法應用在滲流方面的研究主要集中在裂隙滲流方面，很少直接應用于壩體浸潤線的模擬上，而自由面的確定是滲流研究的重要環(huán)節(jié)，只有確定了自由面，才能確定滲流區(qū)的范圍，從而才能得到準確的滲流場。本文嘗試將格子Boltzmann方法用于模擬壩體滲流浸潤線，并與傳統(tǒng)方法模擬結果進行對比，以論證其可行性和可靠性，以期為常規(guī)滲流模擬提供一種新的簡單高效的數(shù)值計算方法。

1 滲流和LBM基本理論

1.1 滲流基本方程

二維滲流的連續(xù)性方程可由下式來描述，它是研究地下水運動的基本方程：

(1)

式中：H為水頭；Kx和Ky分別為x和y方向的滲透系數(shù)(本次研究中滲透系數(shù)取各項同性，即Kx=Ky=K)；α為含水層壓縮系數(shù)；β為水的體積壓縮系數(shù)；n為孔隙度；t為時間；μs為給水度；ρ為水的密度；g為重力加速度。

1.2 二維格子玻爾茲曼模型

系統(tǒng)不受外力時，玻爾茲曼運輸(LBM)方程可以表示為：

▽f=Ω

(2)

式中：f=f(r,c,t)表示t時刻速度為c～c+dc且處于r～r+dr范圍內(nèi)的分子數(shù)，Ω=Ω(f)表示分布函數(shù)由初始狀態(tài)到最終狀態(tài)的變化率，稱為碰撞算子，常用式(3)所示的BGKW近似來模擬。

(3)

式中：τ稱為無量綱弛豫時間；feq為離散速度空間的局部平衡態(tài)分布函數(shù)。

格子玻爾茲曼模型有很多種，采用圖1所示的二維九速(D2Q9)模型，則LBM方程(2)可離散為如下方程：

fi(r+ciΔt,t+Δt)-fi(r,t)=

(4)

式中：r為位置矢量；Δt表示時間步長；ci為D2Q9模型[22]從中心節(jié)點出發(fā)的9個速度向量，分別為c1=(0,0)、c2(1,0)、c3(0,1)、c4(-1,0)、c5(0,-1)、c6(1,1)、c7(-1,1)、c8(-1,-1)、c9(1,1)。

圖1 D2Q9模型Fig.1 D2Q9 model

對于D2Q9模型，feqi(r,t)平衡分布函數(shù)可用下式表示:

feqi=ρwi

(5)

權函數(shù)為：

(6)

為了將格子玻爾茲曼方程恢復到宏觀表達式，可得宏觀滲流方程(1)中的水頭與分布函數(shù)的對應關系：

(7)

相應的滲流速度表達式為：

(8)

宏觀的壓力可以通過能量守恒方程得到，即：

(9)

通過Chapman-Enskog展開，可以建立介觀尺度和宏觀尺度的聯(lián)系，即松弛時間τ與滲透系數(shù)K的關系表示為：

(10)

式中：Δx為格子步長。

1.3 邊界條件

在格子玻爾茲曼方法中，邊界條件起著至關重要的作用，它會對數(shù)值計算的精度、效率和穩(wěn)定性產(chǎn)生極大的影響。本文所采用的邊界條件為動力學格式和反彈格式[23]。

第一類邊界條件采用通量守恒原理。例如在左邊界，對于給定的φ值，存在：feq5-f5+feq7-f7=w5φwall+w7φwall，f7可以由遷移過程獲得，那么f5=w5φwall+w7φwall-f7，同理，也可以計算出其他的未知遷移函數(shù)。

(11)

1.4 模型驗證

根據(jù)以上理論和模型，采用Fortran語言編寫了基于D2Q9的二維LBM計算程序。為了驗證該程序的準確性，采用一維熱傳導模型進行驗證。模型長L=4.0。左邊界x=0處保持恒溫T1=6，右邊界x=L處保持恒溫T2=0。

這一問題的解析解為如下形式：

(12)

根據(jù)上述邊界條件編寫代碼，擴散系數(shù)α=0.25，分別取t=1 min和t=10 h，進行試算。圖2為計算得到的不同時刻的總水頭等勢線，與解析解結果吻合良好，說明采用D2Q9模型編制的程序準確，模型精度較高。

圖2 不同時刻水頭分布沿x方向?qū)Ρ菷ig.2 Water head distribution in different time along the x direction

2 LBM在土石壩滲流計算中的應用

2.1 帶排水井矩形壩體滲流計算

2.1.1 物理模型

選取的矩形壩模型如圖3所示，寬40.0 m，上游水深為60.0 m,下游水深為0 m。為簡化邊界條件處理，認為壩體高也為60.0 m。D點處設置一井，井的流量為0.01 m3/s，格子大小 ，因為格子玻爾茲曼方法不需要考慮穩(wěn)定性問題，所以為簡單起見，計算時間步長設為1 s。

圖3 矩形壩體模型Fig.3 Rectangular dam model

進行邊界處理時，上邊界AB按照反射格式處理。左邊界AC給定了邊界函數(shù)值，即H=60 m。下邊界CD和右邊界BD按照不透水情況處理。D點處因為有井，所以需要進行特殊處理，本次計算中按源項處理[24]。

2.1.2 計算結果

將編制的基于D2Q9模型的LBM計算程序用于上述矩形壩體滲流問題中，可得到不同時刻的滲流場，圖4為時間t=1 000 min時的總水頭和壓力水頭等勢線計算結果，與有限元軟件(SEEP/W)計算得到的結果對比情況見圖5。圖5中虛線為LBM計算結果，實線為SEEP/W的計算結果。

根據(jù)圖4和圖5，可以發(fā)現(xiàn)壩體內(nèi)水頭和壓力等值線沿程呈下降趨勢，浸潤線(壓力為0)沿程降落迅速，未在下游壩坡出逸，由于下游排水井的影響，流場方向在排水井附近呈輻射狀，與實際物理現(xiàn)象符合。LBM和SEEP/W計算出的總水頭和壓力水頭等勢線吻合良好，誤差極小，證明了LBM在矩形土石壩滲流模擬中具有較高的可靠性，而LBM方法為顯式方法，無須有限元方法的多次迭代求解，計算效率更高。

圖4 t=1 000 min時LBM模擬總水頭和壓力水頭等勢線結果Fig.4 t=1 000 min LBM simulate total head and pressure head results

圖5 LBM(虛線)和SEEP/W(實線)模擬總水頭和壓力水頭結果比較Fig.5 LBM(dotted line) and SEEP/W(solid line) simulate total head and pressure head results

2.2 梯形壩體滲流計算

2.2.1 物理模型

為了驗證LBM在非規(guī)則計算區(qū)域中的適應性，選取了如圖6所示梯形壩模型，高12 m，長52 m，上下游壩坡均為1∶2，壩基鋪設長為12 m的水平排水褥墊。為簡化邊界條件，取上游水位與壩高相等，即12 m，下游水位為0 m，初始時內(nèi)壩體內(nèi)及下游無水。

由于LBM采用的為正方向網(wǎng)格，在進行模擬時，需要對計算區(qū)域邊界進行識別，區(qū)域Ⅰ采取第一類邊界條件處理方式，區(qū)域Ⅱ采取第二類邊界條件處理方式。

圖6 梯形壩體模型Fig.6 Trapezoidal dam model

2.2.2 計算結果

不同時刻得到的壩體內(nèi)水頭和滲透壓力等值線見圖7～圖10?？梢钥闯?，模擬初期，壩體內(nèi)水頭等值線非常集中，降落迅速，壓力等值線和水頭等值線在水頭陡降區(qū)接近平行，隨著模擬時間的增加，水頭等值線和壓力等值線也趨于正交，浸潤線逐漸抬升，到5 h時左右，壩體內(nèi)滲流場趨于穩(wěn)定。

圖7 t=1 min時LBM模擬總水頭和壓力水頭等勢線結果Fig.7 t=1min LBM simulate total head and pressure head results

圖8 t=20 min時LBM模擬總水頭和壓力水頭等勢線結果Fig.8 t=20 min LBM simulate total head and pressure head results

圖9 t=5 h時LBM模擬總水頭和壓力水頭等勢線結果Fig.9 t=5 hours LBM simulate total head and pressure head results

圖10 t=1 min, t=20 min, t=5 h 的浸潤線Fig.10 t=1 min, t=20 min, t=5 h saturation line

圖11為SEEP/W計算的穩(wěn)定滲流場計算結果，圖12顯示了LBM和SEEP/W計算得到水頭等值線結果對比情況，虛線為LBM計算結果，實線為SEEP/W計算結果。LBM和SEEP/W計算得出的滲流場吻合良，說明了本文采用的LBM模型在計算不規(guī)則邊界時采用的處理方法可行，計算準確性有保證。

圖11 SEEP/W模擬總水頭和壓力水頭等勢線結果Fig.11 SEEP/W simulate total head(a)and pressure head(b)results

圖12 LBM(虛線)和SEEP/W(實線)模擬的壓力等值線結果Fig.12 LBM(dotted line) and SEEP/W(solid line) simulate pressure contour results

3 結 論

格子玻爾茲曼方法已經(jīng)成為流體力學計算中的一個備受矚目的新方法，該方法可以適用于復雜的邊界條件和描述復雜系統(tǒng)的演化過程，并且容易編程，計算效率高。本文運用LBM分別模擬了具有規(guī)則邊界的矩形壩體和不規(guī)則邊界的梯形壩體的滲流計算，準確得到了壩體內(nèi)的滲流場及浸潤線等，與傳統(tǒng)有限元方法得到計算結果相一致，表明了LBM在壩體滲流方面具有較高的可行性和準確性，與傳統(tǒng)方法相比，LBM邊界處理簡單方便，計算效率較高，程序編制容易且易于推廣到三維問題中，是一種值得在滲流計算中推廣的方法。

□

[1] 孫繼昌.中國的水庫大壩安全管理[J].中國水利,2008,(20):10-14.

[2] 張 斌. 土石壩滲流研究的現(xiàn)狀與發(fā)展[J]. 河南水利與南水北調(diào),2016,(5):111-112.

[3] 鮑呈蒼,基于SEEP/W分析土石壩滲流穩(wěn)定[J],中國水運,2011,11(6):147-148.

[4] 張 嘎,張建民,江春波. 溪洛渡水電站上游圍堰滲流分析及防滲型式比較[J].水力發(fā)電學報,2002,(S1):82-87.

[5] 陳玉茹. 基于土石壩滲流和壩坡穩(wěn)定的二維有限元分析[J]. 浙江水利水電學院學報,2014,(4):13-15.

[6] 歐陽君,徐千軍,侍克斌,等. 土石壩邊坡穩(wěn)定性分析的溫控參數(shù)折減有限元法[J]. 巖土力學,2011,(8):2 549-2 554.

[7] 楊林德,楊志錫.各向異性飽和土體的滲流耦合分析和數(shù)值模擬[J].巖石力學與工程學報,2002,21(10):1 447-1 451.

[8] 付成華,陳敏林,周洪波. 基于有限元法的土石壩非穩(wěn)定滲流壩體浸潤線簡化計算條件分析[J]. 長江科學院院報,2004,(4):60-62.

[9] 倪小東,王 媛. 運用離散單元法研究土體管涌機理初探[C]∥ 水工滲流研究與應用進展,2006:37-42.

[10] XuY Q, UnamiK, KawachiT. Optimal hydraulic design of earth dam cross section using saturated-unsaturated seepage flow model [J]. Advances in Water Resources, 2003,26.

[11] Xu YQ,Unami K, Kawachi T. Optimal hydraulic design of earth dam cross section using saturated-unsaturated seepage flow model [J].Advances in Water resources,2003,26(1).

[12] Lam L. Fredlund D G. Saturated-Unsaturated Transient Finite Element Seepage Model for Geotechnical [J]. Advances in Water Resources, 1984,7(3).

[13] 高朝軍,阿布都熱西提·阿布都外力. 基于邊界元法的滲流自由面計算[J]. 水資源與水工程學報,2008,(6):32-35.

[14] 孫偉建,侯興民,李遠東,等. 基于虛單元法求解滲流自由面的曲線擬合法[J]. 水力發(fā)電,2016,(11):50-67.

[15] 華 靜,楊華舒. 土石壩滲流計算中的有限元應用研究[J]. 中國水能及電氣化, 2012,(7):15-18.

[16] 張 霞,李鳳濱,盛金昌,等. 格子Boltzmann方法模擬土體滲流場研究[J]. 水電能源科學,2012,(10):61-64.

[17] Hazlett R. Simulation of capillary-dominated displacements in microtomo grapphic images of reservoir rocks [J]. Trans Porous Media, 1995,20:21-35.

[18] 姚 軍,趙秀才,衣艷靜,等.儲層巖石微觀結構性質(zhì)的分析方法[J].中國石油大學學報(自然科學版)，2007,31:80-86.

[19] Degruyter W, Burgisser A, Bachmann O, et al. Synchrotron X-ray microtomography and lattice Boltzmann simulations of gas flow through volcanic pumices [J].Geosphere,2010,(6):470-481.

[20] 樊 火,鄭 宏,基于MRT-LBM 的分形裂隙網(wǎng)絡滲流數(shù)值模擬[J].中國科學：技術科學，2013,43(12):1 338-1 345.

[21] 田智威,譚云亮,劉兆霞.含裂隙煤體瓦斯?jié)B流規(guī)律的LBM數(shù)值模擬[J].煤炭學報,2013,(8):1 376-1 380.

[22] 尚 智,陳 碩,趙 鈞. Lattice Boltzmann方法不同邊界條件的混用[J]. 上海交通大學學報,2006,(10):1 800-1 805.

[23] 聶德明,林建忠. 格子Boltzmann方法中的邊界條件[J]. 計算物理,2004,(1):21-26.

[24] 黃俊濤,張 力,雍穩(wěn)安,等. 格子Boltzmann方法解擴散方程的復雜邊界條件研究[J]. 應用數(shù)學和力學,2014,(3):305-312.