王 敏，李銘華

(1.長江勘測規(guī)劃設(shè)計研究有限責(zé)任公司，湖北 武漢 430010；2.南京市長江河道管理處，江蘇 南京 210011)

在大型水電和地?zé)崮茉撮_發(fā)，石油、天然氣、地下水開采及核廢料深層埋置等工程問題中，均涉及裂隙巖體滲流及水力耦合機(jī)理的研究，巖石裂隙細(xì)觀流動特性是上述機(jī)理研究的關(guān)鍵問題之一[1-4]。

巖體裂隙中的細(xì)觀流動滿足不可壓縮流體動量守恒方程，即Navier-Stokes(N-S)方程：

(1)

(2)

式中，Q為流體流量，e為平板裂隙開度，J為水力梯度。天然裂隙表面凹凸不平，裂隙結(jié)構(gòu)之間的開度也并非定值，因此立方定律只適用于裂隙開度變化較緩或局部區(qū)域，公式(2)也被稱為局部立方定律。

式(2)中的裂隙開度e為物理寬度，也稱機(jī)械開度或力學(xué)開度，為考慮裂隙粗糙度的影響，常用等效水力開度eh代替裂隙機(jī)械開度，不同學(xué)者通過大量室內(nèi)裂隙滲流試驗(yàn)提出了各自的機(jī)械開度與等效水力開度之間的換算公式[5-7]。為有效地建立機(jī)械開度與水力開度的關(guān)系，通常需要引入描述裂隙粗糙度參數(shù)，如裂隙粗糙度系數(shù)Jrc，為考慮裂隙剪切錯動對裂隙等效水力開度的影響，Olsson等[8]通過實(shí)驗(yàn)研究提出了在剪切錯動條件下，機(jī)械開度與水力開度的換算關(guān)系：

(3)

式中：ds為裂隙剪切位移，dsp為峰值剪應(yīng)力時的剪切位移，Jrc 0為裂隙的初始粗糙度系數(shù)，Jrc mob為剪切錯動后的粗糙度系數(shù)。

從式(3)中可以看出，剪切位移對裂隙等效水力開度影響較大，但式(3)未考慮裂隙流態(tài)的影響，且公式是根據(jù)實(shí)驗(yàn)數(shù)據(jù)擬合得到，未揭示剪切位移對裂隙滲流影響的機(jī)理。本文采用格子玻爾茲曼模型(Lattice Boltzmann Method，簡稱LBM)[9-10]模擬單裂隙細(xì)觀滲流過程，研究在剪切條件下，考慮不同流速(不同雷諾數(shù))對細(xì)觀流動特性的影響，研究結(jié)果將為建立機(jī)械開度和等效水力開度的宏觀規(guī)律提供基礎(chǔ)。

1 數(shù)值模擬方法

1.1 MRT-LBM格子模型

在格子玻爾茲曼模型中，多松弛模型(Multiple-Relaxation-Time，簡稱MRT)由于采用多個因子控制LBM中的碰撞過程，能有效克服單松弛模型帶來的數(shù)值穩(wěn)定性差且存在滲透率與黏性相關(guān)等不足，并且能計算中等雷諾數(shù)條件下的流動過程[11]。在MRT-LBM格子模型中，速度分布函數(shù)f滿足如下演化方程[12]：

fi(x+eiΔt,t+Δt)-fi(x,t)=-M-1SM[f(x,t)-

feq(x,t)]

(4)

fi(x+eiΔt,t+Δt)-fi(x,t)=-M-1S[m-meq]

(5)

式中:S為非負(fù)松弛參數(shù)矩陣，ei為i方向的權(quán)函數(shù)，m=Mf，M為空間轉(zhuǎn)化矩陣，feq為平衡態(tài)速度分布函數(shù)。

1.2 方法驗(yàn)證

采用Fortran 90語言編寫LBM計算程序[13]，并采用圓柱繞流模型驗(yàn)證程序的有效性和正確性。圓柱繞流模型的進(jìn)口采用速度邊界，采用Zou等[14]格式實(shí)現(xiàn)，出口為充分發(fā)展邊界，上下邊壁及中間圓柱為速度無滑移邊界，采用非平衡外推格式[10-11]。取雷諾數(shù)為10、20和40三個計算工況。

雷諾數(shù)Re為40時的流速分布如圖1所示，由圖可見，流速在x方向和y方向均具有較好的對稱性，且在圓柱邊壁流速為零，上述結(jié)果表明計算方法在曲面邊界處理方面具有良好的計算穩(wěn)定性。

圓柱繞流模型在不同雷諾數(shù)條件下的流線如圖2所示。由圖2可見，圓柱后方均出現(xiàn)穩(wěn)定的渦流，且漩渦面積隨著雷諾數(shù)的增大而增大。為驗(yàn)證本文計算程序的正確性，采用渦流的相對長度2Leddy/D和分離角度θs兩個無量綱參數(shù)定量描述渦流區(qū)域的幾何特性，其中Leddy為渦流區(qū)域最遠(yuǎn)點(diǎn)到圓柱的距離，并與Ding等[15]的計算結(jié)果進(jìn)行對比，對比結(jié)果如表1所示，在不同雷諾數(shù)條件下，本文的圓柱繞流模擬結(jié)果與前人的模擬結(jié)果基本一致，表明本文模擬方法及計算程序的正確性成立。

圖2 不同雷諾數(shù)條件下的流線圖

2 剪切條件下單裂隙流動特性分析

2.1 物理模型

為獲取真實(shí)裂隙結(jié)構(gòu)面模型，首先采用巴西劈裂法[16]對巖體試樣進(jìn)行人工劈裂，得到粗糙裂隙面，再利用三維非接觸式輪廓儀掃描裂隙面，由采集的裂隙面數(shù)據(jù)生成的三維形貌如圖3所示，裂隙面數(shù)據(jù)密度為0.2 mm，Z方向的數(shù)據(jù)精度為10 μm。

如圖3所示，選取掃描裂隙面上的輪廓線A-A′為作為二維裂隙模型的下輪廓線，并將A-A′輪廓線豎直平移1 mm形成上輪廓線，得到二維裂隙模型，記為Model 0，如圖4所示。將上輪廓線向水平方向分別平移1.0 mm、2.0 mm、3.0 mm和4.0 mm，可得到剪切位移ds分別為1.0 mm、2.0 mm、3.0 mm和4.0 mm條件下的二維裂隙幾何模型，分別記為Model 1—Model 4。對裂隙模型Model 0—Model 4的力學(xué)開度進(jìn)行統(tǒng)計分析，結(jié)果如表2所示，對于完全吻合的裂隙模型Model 0，裂隙上下兩條輪廓線的相關(guān)系數(shù)為1.0，隨著剪切位移的增加，吻合程度降低，相關(guān)系數(shù)逐漸減小，裂隙開度的均方差也逐漸增大，說明裂隙開度的變異性隨剪切位移的增大而逐漸增大。

圖3 裂隙試樣的三維形貌

圖4 單裂隙幾何模型

2.2 模擬條件

表2 裂隙模型力學(xué)開度統(tǒng)計

2.3 流速分布

各裂隙模型B-B′截面(位置見圖4)處流速隨剪切位移的演化如圖5所示，圖5分別給出了雷諾數(shù)為1和100兩個計算工況下的結(jié)果。從圖5中可以看出，在不同雷諾數(shù)條件下，B-B′截面均存在一定大小的y方向流速uy，且速度大小受開度變化的影響，可能存在正負(fù)兩個方向的uy速度；當(dāng)Re=1時，即低雷諾數(shù)條件下，流速ux分布為典型的拋物線分布，且不同剪切位移ds條件下的ux速度分布形式基本一致，x方向的流速的最大值與截面處的裂隙開度成反比；當(dāng)Re=100時，即高雷諾數(shù)條件下，ux流速在靠近裂隙邊壁處的低流速區(qū)域有所增加，同時在剪切位移ds=4.0 mm的條件下，ux速度在一定寬度范圍內(nèi)出現(xiàn)負(fù)方向的流速，說明裂隙內(nèi)出現(xiàn)明顯的回流現(xiàn)象。

為揭示在高雷諾數(shù)條件下裂隙的細(xì)觀流動特征，圖6給出了不同剪切位移條件下的流線圖。從圖6中可見，當(dāng)ds=0 mm時，裂隙內(nèi)流線與流動主通道基本平行，隨著剪切位移值ds逐漸增加，在裂隙內(nèi)逐漸出現(xiàn)了大小不一的渦流，主要出現(xiàn)在開度變化劇烈的區(qū)域，且渦流數(shù)量隨著剪切位移的增加而增加，渦流大小也逐漸增大。當(dāng)ds=4.0 mm時，裂隙內(nèi)的渦流增加較為明顯，局部區(qū)域的渦流逐漸占據(jù)裂隙滲流的主通道，說明在高雷諾數(shù)條件下，剪切位移帶來的裂隙開度變化極大地影響裂隙通道內(nèi)的流動狀態(tài)。

2.4 渦流的形成和演化

圖7給出裂隙模型Model 4(ds=4.0 mm)在不同雷諾數(shù)條件下局部區(qū)域內(nèi)(x=50 mm～60 mm)的流線圖。由圖7中可見，當(dāng)雷諾數(shù)較小時(Re=1)，裂隙內(nèi)無渦流產(chǎn)生，隨著雷諾數(shù)增大，渦流才開始出現(xiàn)；當(dāng)Re=50～100，渦流的數(shù)量開始增多，渦流的面積也逐漸增大，導(dǎo)致裂隙內(nèi)的有效過流區(qū)域減小；當(dāng)Re=100～150，渦流面積在增大的同時，一些較大的渦流會在裂隙邊壁分裂出較小的次生渦流，原來一些較小的渦流會逐漸合并形成范圍更大的渦流，說明此時裂隙內(nèi)的形成了較強(qiáng)的紊流流態(tài)。

圖5 速度隨剪切位移的演化圖(B-B′：x=55 mm)

在高雷諾數(shù)條件下(Re=100)，各裂隙模型(ds=1.0 mm～4.0 mm)在局部區(qū)域內(nèi)(x=50 mm～60 mm)的流線如圖8所示。由圖8可見，裂隙開度分布隨著剪切位移的增加而變化，逐漸形成一個開度先縮小再突然擴(kuò)張的流動通道，在剪切位移逐漸增加的過程中，首先在裂隙開度擴(kuò)張區(qū)域的下壁面附近形成一個較小的渦流(ds=1 mm)，隨之渦流面積逐漸增大，同時在裂隙開度擴(kuò)張區(qū)域的上壁面也形成渦流(ds=2 mm)，當(dāng)剪切位移進(jìn)一步增加時(ds=3 mm～4 mm)，上下壁面的渦流區(qū)域沿順流方向逐漸增加擴(kuò)大，并交替占據(jù)裂隙內(nèi)主要流動通道。

圖6 不同剪切位移條件下Re=100.0時流線圖

圖7 渦流的形成和發(fā)展

2.5 水力開度

受渦流的的影響，裂隙內(nèi)的有效過流面積減小，從而減弱裂隙的過流能力，導(dǎo)致裂隙的等效水力開度降低。各裂隙模型的等效水力開度與平均機(jī)械開度比值eh/e隨雷諾數(shù)Re的變化如圖9所示，從圖中可以看出，在雷諾數(shù)較小時，即Re<1～10，裂隙的等效水力開度基本維持初始的等效水力開度不變，當(dāng)雷諾數(shù)增大時，即Re>10，裂隙的等效水力開度開始逐漸變小，且減小趨勢隨著雷諾數(shù)的增大而增大。同時，對各剪切位移的裂隙模型，裂隙剪切位移越大，裂隙的初始水力開度越小，裂隙的等效水力開度隨雷諾數(shù)的變化趨勢也越明顯。

圖8 渦流隨剪切位移的變化

圖9 水力開度與平均機(jī)械開度比值隨Re的變化

3 結(jié) 論

本文采用格子玻爾茲曼方法模擬在剪切錯動條件下巖石裂隙的滲流過程，考慮不同雷諾數(shù)對細(xì)觀流動特性的影響，主要結(jié)論如下：

(1) 發(fā)展了裂隙細(xì)觀流動特性數(shù)值模擬的格子Boltzmann方法，并采用圓柱繞流模型驗(yàn)證數(shù)值模擬方法的正確性和穩(wěn)定性。

(2) 數(shù)值模擬結(jié)果表明剪切錯動在高雷諾數(shù)條件下可導(dǎo)致裂隙產(chǎn)生明顯的非線性流動現(xiàn)象，主要表現(xiàn)為隨著雷諾數(shù)的增加，在靠近裂隙邊壁處低流速區(qū)域有所增加，同時出現(xiàn)負(fù)方向的流速，說明裂隙內(nèi)出現(xiàn)明顯的渦流區(qū)域，且非線性流動隨著剪切位移的增加增強(qiáng)。

(3) 數(shù)值模擬結(jié)果表明在剪切過程中，隨著雷諾數(shù)的增大，渦流的數(shù)量逐漸變多，渦流的范圍逐漸擴(kuò)大，從而導(dǎo)致有效過流通道的減小，裂隙的水力開度減小，且減小的趨勢隨著雷諾數(shù)的增大而增大。

(4) 計算結(jié)果揭示了在剪切錯動和流速增大條件下裂隙渦流的形成和演化過程，及渦流造成裂隙有效對流區(qū)減小、過流能力降低的細(xì)觀機(jī)制，對深化巖土介質(zhì)滲流和傳熱過程的認(rèn)識具有重要意義。