剪切條件下單裂隙細觀流動特性研究

王 敏，李銘華

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

在大型水電和地熱能源開發，石油、天然氣、地下水開采及核廢料深層埋置等工程問題中，均涉及裂隙巖體滲流及水力耦合機理的研究，巖石裂隙細觀流動特性是上述機理研究的關鍵問題之一[1-4]。

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

(1)

(2)

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

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

(3)

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

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

1 數值模擬方法

1.1 MRT-LBM格子模型

在格子玻爾茲曼模型中，多松弛模型(Multiple-Relaxation-Time，簡稱MRT)由于采用多個因子控制LBM中的碰撞過程，能有效克服單松弛模型帶來的數值穩定性差且存在滲透率與黏性相關等不足，并且能計算中等雷諾數條件下的流動過程[11]。在MRT-LBM格子模型中，速度分布函數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為非負松弛參數矩陣，ei為i方向的權函數，m=Mf，M為空間轉化矩陣，feq為平衡態速度分布函數。

1.2 方法驗證

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

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

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

圖2 不同雷諾數條件下的流線圖

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

2.1 物理模型

為獲取真實裂隙結構面模型，首先采用巴西劈裂法[16]對巖體試樣進行人工劈裂，得到粗糙裂隙面，再利用三維非接觸式輪廓儀掃描裂隙面，由采集的裂隙面數據生成的三維形貌如圖3所示，裂隙面數據密度為0.2 mm，Z方向的數據精度為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的力學開度進行統計分析，結果如表2所示，對于完全吻合的裂隙模型Model 0，裂隙上下兩條輪廓線的相關系數為1.0，隨著剪切位移的增加，吻合程度降低，相關系數逐漸減小，裂隙開度的均方差也逐漸增大，說明裂隙開度的變異性隨剪切位移的增大而逐漸增大。

圖3 裂隙試樣的三維形貌

圖4 單裂隙幾何模型

2.2 模擬條件

表2 裂隙模型力學開度統計

2.3 流速分布

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

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

2.4 渦流的形成和演化

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

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

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

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

圖7 渦流的形成和發展

2.5 水力開度

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

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

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

3 結 論

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

(1) 發展了裂隙細觀流動特性數值模擬的格子Boltzmann方法，并采用圓柱繞流模型驗證數值模擬方法的正確性和穩定性。

(2) 數值模擬結果表明剪切錯動在高雷諾數條件下可導致裂隙產生明顯的非線性流動現象，主要表現為隨著雷諾數的增加，在靠近裂隙邊壁處低流速區域有所增加，同時出現負方向的流速，說明裂隙內出現明顯的渦流區域，且非線性流動隨著剪切位移的增加增強。

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

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