劉嘉夫，齊 昕

（重慶水利電力職業(yè)技術學院水利工程學院，重慶402160）

1 前 言

近年來，隨著水電工程和地下空間利用等大型巖體工程的興建，巖體滲流問題日益得到關注。在長期地質作用下，巖體中普遍存在各種不同尺寸的裂隙，由于裂隙是各種流體在巖體中流動的主要通道，對地下水流動起著主要的控制作用，因此關于裂隙中的流體滲流分析是近年來巖體滲流分析的熱點之一。在早期的研究中，一般忽略裂隙上、下表面的粗糙性，將粗糙裂隙視為理想光滑平行板模型［1-4］，裂隙的滲透性可以用立方定理（cubic law）來描述［5］，即通過裂隙的滲流體積流量與裂隙隙寬的三次方成正比，由于忽略了真實裂隙的粗糙性，因此立方定理只能定性地描述裂隙間的滲流情況。若考慮裂隙上、下表面的粗糙性以及相互接觸部分的影響，裂隙中的流體流動將變得迂回曲折，從而使得真實的滲流體積流量與立方定理計算值之間存在偏差。鑒于此，很多學者通過在立方定理中引入經驗修正系數(shù)f來描述粗糙裂隙與理想光滑平行板模型之間的偏差，由于這些研究只是考慮裂隙表面的粗糙性對裂隙滲流的影響，沒有考慮裂隙隙寬分布的不均勻性對滲流的顯著影響，因此引入修正系數(shù)f不能完全反映裂隙中的滲流情況。

S.R.Brown［6］采用分形函數(shù)構造粗糙裂隙的上、下表面，然后將兩表面疊加到一起構成裂隙空腔，采用有限差分方法求解裂隙面上的滲流體積流量分布情況，通過滲流體積流量分布情況可以看到裂隙中的滲流呈現(xiàn)出明顯的曲折效應。法國數(shù)學家B.B.Mandelbrot在1973年首次提出了分維的設想，創(chuàng)造了“分形（Fractal）”這個新術語，后來其又提出了分形幾何，用來描述自然界不規(guī)則以及雜亂無章的現(xiàn)象和行為［7］。自然界的大多數(shù)復雜幾何形狀都具有分形特性，巖石裂隙表面形狀也同樣具有分形特性，可以用分形維數(shù)來描述巖石裂隙表面的粗糙性。Y.H.Lee等［8］應用分形幾何的尺碼方法測量了裂隙剖面的分形維數(shù)；N.H.Maerz等［9-10］根據(jù)經驗建立了 JRC值與分形維數(shù)之間的關系式；S.Murata 和 T.Saito［11］研究了裂隙面分形參數(shù)對曲折效應的影響規(guī)律。由于分形維數(shù)僅是表征裂隙表面形貌特征的參數(shù)，而裂隙滲流的曲折性不僅與裂隙的表面形貌相關，也與裂隙的三維空腔組合形貌密切相關，因此結合裂隙三維空腔組合形貌對裂隙滲流的曲折性進行定量描述是研究裂隙滲流曲折性的新思路。

綜上所述，為了在粗糙裂隙滲流計算公式中考慮裂隙滲流曲折效應的影響，在對粗糙裂隙試件進行室內滲流試驗和有限元數(shù)值方法求解的基礎上，結合裂隙的三維空腔組合形貌，提出了裂隙滲流曲折因子的計算方法；然后借助分析巖石滲透系數(shù)的等效溝槽模型，在現(xiàn)有的立方定理修正計算公式中引入了曲折因子，得到了考慮裂隙曲折效應的滲流計算新公式，在利用數(shù)值計算結果確定新公式中的待定常數(shù)后，根據(jù)室內滲流試驗結果對新公式進行了驗證。

2 考慮曲折效應的裂隙滲流計算公式的推導

裂隙中的滲流流線都是曲線，這與光滑平行板模型中的直線型流線分布是完全不同的，表明裂隙的實際流徑要大于從入水口到出水口的直線距離，這主要是由于粗糙裂隙隙寬分布不均勻引起裂隙滲流流徑呈現(xiàn)出曲折特性。光滑裂隙的滲流可以用立方定理來描述，其表達式為

式中：Qc為裂隙的滲流體積流量；b為裂隙隙寬；?p為水頭壓力梯度；C為與裂隙幾何尺寸和流體物理特性有關的常數(shù)。

對于粗糙裂隙，考慮裂隙表面粗糙性影響的立方定理修正經驗公式可以寫為

式中：Q為考慮裂隙表面粗糙性的滲流體積流量；f為考慮裂隙表面粗糙性的修正系數(shù)；＜b＞為粗糙裂隙平均隙寬。

對式（2）進行適當?shù)淖儞Q，有

由式（3）可以看出，考慮裂隙表面粗糙性影響對立方定理進行修正，其本質是考慮裂隙面起伏不平對粗糙裂隙平均隙寬進行折減，從而將粗糙裂隙簡化為隙寬 ＜b＞／f1／3的光滑平行板模型，以滿足立方定理的計算條件。

在式（2）中，關于考慮裂隙表面粗糙性的修正系數(shù) f，G.M.Lomize 等［12-15］提出了不同的表達式，將這些表達式歸納起來可用統(tǒng)一的表達式表示：

式中：A和B均為待定常數(shù)，可以根據(jù)滲流試驗結果確定；Δ為裂隙表面絕對粗糙度。

在式（3）中引入修正系數(shù)F得到

對于平行板模型，有

式中：k為巖石的滲透系數(shù)；kc為按立方定理進行計算所得的滲透系數(shù)；τ為滲流平均曲折因子；Ts為裂隙面曲折系數(shù)。

式（6）表示粗糙裂隙滲流與修正立方定理式（2）之間的偏差，與式（5）中的系數(shù)1／F等效，即有

將式（7）代入式（5）中，有

將式（4）代入式（8）中進行整理可得

式（9）即為考慮裂隙面粗糙性和滲流曲折效應的粗糙裂隙滲流計算公式，其中曲折因子τ和裂隙面曲折系數(shù)Ts是分別根據(jù)裂隙三維空腔組合形貌和裂隙表面形貌特征計算得到的，只要確定了待定常數(shù)A和B，就可以將式（9）應用于粗糙裂隙的滲流計算。

3 裂隙面量測與分形計算

3.1 裂隙面制備

試驗中，裂隙面采用劈裂巖石作為原型，首先用不易變形且強度很高的樹脂材料拷貝裂隙表面的形貌，然后基于該樹脂試件表面形貌制作石膏材料試件，用來模擬天然的巖石試件。類巖石材料的成分為石膏、水和延緩劑。試件分為上、下兩部分，長200 mm、寬100 mm，上、下部分各高50 mm，總高為100 mm。由于上、下裂隙面由新鮮巖塊劈裂而來，試件吻合較好，因此初始接觸可近似視為1.0。

3.2 裂隙表面形貌量測

試驗選取的5組粗糙表面巖石裂隙試件分別標記為J1、J2、J3、J4和J5，其中裂隙試件J1和J2的數(shù)據(jù)來源于 B.Li等的研究成果［16］，裂隙試件 J3、J4 和 J5 的數(shù)據(jù)通過進行新的室內物理模型試驗獲取。試件表面形貌用KEYENCE公司生產的激光測試儀量測，該儀器精度可達到±20μm，分辨率可達到10μm，掃描儀可識別X、Y方向平面定位坐標系統(tǒng)，并且可根據(jù)掃描過程中預設路徑自行移動，掃描數(shù)據(jù)通過PC機可實現(xiàn)實時收集與處理。本研究中，各巖石裂隙表面在X、Y軸方向上測點間距均為0.01 mm，得到裂隙試件三維表面形貌如圖1所示。可以看出，試件J1表面平滑，整體不存在較大凸起；試件J2局部存在某些較大凸起，但整體呈平滑狀；試件J3表面較粗糙，但整體沒有較大凸起；試件J4表面較平滑，整體沒有較大凸起，但存在比較多的小凸起；試件J5表面較粗糙，有較大凸起，同時也存在較多的小凸起。幾組裂隙試件在試驗過程中均施加1 MPa恒定法向應力，剪切位移達18 mm，且每剪切1 mm便在裂隙內施加與剪切方向一致的大小為0.1 m的水頭進行透水試驗。隨著剪切位移的增加，上、下裂隙面相互錯動導致滲流路徑減小，在水頭保持不變的情況下，兩端的水力梯度會逐漸增大，為了消除該影響，在計算中需降低與剪切位移相對應的水力梯度。

圖1 裂隙面形貌

3.3 裂隙分形描述與計算

根據(jù)裂隙試件表面形貌測試數(shù)據(jù)，計算其分形維數(shù)。分形幾何中，最常用的分形維數(shù)計算方法是尺碼法和覆蓋法。其中尺碼法可以直接估算出復雜曲線的分形維數(shù)，而對于粗糙表面則需要采用間接覆蓋的方法，三角形棱柱表面積法［17］和投影覆蓋法［18］是具有代表性的兩種計算方法，但計算過程中對面積的近似計算使得這兩種方法的計算結果存在一定的偏差。針對該問題，周宏偉等［19］提出了表面分維的立方體覆蓋法，采用三維立方體網(wǎng)格直接覆蓋粗糙表面，具體計算過程如下。

（1）在XOY平面上存在長度為δ的正方形網(wǎng)格，該網(wǎng)格對應的四角點高度依次記為 h（i，j） 、 h（i，j+1） 、 h（i+1，j） 以及 h（i+1，j+1） ，其中 i、j取值為1～n-1（n為每條邊上的測點個數(shù)）。

（2）用長度為δ的立方體覆蓋粗糙表面，對應區(qū)域內覆蓋粗糙表面所需的立方體總數(shù)Ni，j為

（3）計算整個區(qū)域上的立方體總數(shù)

（4）改變測量尺度，再次計算立方體總數(shù)N，由分形理論可知，整個區(qū)域上的立方體總數(shù)N和對應的測量尺度δ之間有如下關系式

式中：k為比例系數(shù)；D為分形維數(shù)，反映了復雜形體占有空間的有效性，是復雜形體不規(guī)則性的量度。

由式（12）可知

為了簡捷表示，分形維數(shù)D的計算公式可以寫成

式中：δ0、N0為常數(shù)，其中

通過編寫Matlab程序，可以計算得出粗糙表面試件 J1、J2、J3、J4、J5的立方體覆蓋數(shù)目與觀測尺度之間的關系，如圖2所示，圖中直線斜率的絕對值即為粗糙表面分形維數(shù)。本文立方體測量尺度為0.2～100.0 mm，當測量尺度較大， lg（δ／δ0） ＞-1.69 時，粗糙表面分形維數(shù)精確到2.000，即粗糙表面不表現(xiàn)出分形特性，只有當測量尺度較小， lg（δ／δ0） ＜-1.69 時，粗糙表面才表現(xiàn)出分形特性。

圖2 立方體覆蓋數(shù)目與觀測尺度之間的關系

4 粗糙裂隙非線性流動特性

4.1 粗糙裂隙飽和滲流試驗

整個試驗裝置包括裂隙模型、上下游緩水箱、增壓水泵、水槽、連接水管、閥門、壓力表等。考慮到高水力梯度下水流流量大這一特點，采用了整體自循環(huán)水流系統(tǒng)，試驗裝置如圖3所示。裂隙模型由上、下兩塊3 cm厚有機玻璃板疊合而成，有機玻璃板具有足夠的剛度以抵抗高水壓引起的變形。J1～J5裂隙開度分別為1.58、2.20、1.50、2.45、2.10mm，裂隙開度采用千分表量測，共測8個點取其平均值。自循環(huán)系統(tǒng)中，設計了增壓泵以提供上游高水頭，增壓泵揚程30m，在H＝25m時排水流量達到250cm3／min，能滿足試驗中供壓穩(wěn)定的要求。水泵與上游緩水箱之間安裝控制閥門來控制流量，以達到調節(jié)水力梯度的目的。為了進行高水力梯度下的水流試驗，在裂隙模型進出口處分別連接封閉式緩水箱，在水箱頂蓋上安裝排氣閥，通過排氣口加壓，對模型和連接水箱排氣。平行板之間以及平行板與緩水箱之間通過γ夾固定連接。在試驗水槽中設置三角堰以觀測流量，三角堰前設一水位測量管，通過水位測針測讀堰上水位，精度達到0.1mm。若三角堰堰口超出水位小于1cm，則直接采用體積法測流量，測讀3次，取平均值。

4.2 非線性分形模型的驗證

4.2.1 系數(shù) A 和 B 的確定

圖3 試驗裝置

在求解待定常數(shù)A和B時，首先根據(jù)裂隙平均隙寬按立方定理計算滲流體積流量Qc，將其與數(shù)值計算結果進行比較，得到兩者的偏差。然后根據(jù)裂隙三維表面形貌和三維空腔組合形貌計算裂隙滲流曲折因子、裂隙面曲折系數(shù)和相對粗糙度，在此基礎上對待定常數(shù)A和B進行擬合分析。根據(jù)裂隙試件J1、J2及J3形貌數(shù)據(jù)計算得到的曲折效應參數(shù)見表1，擬合得到常數(shù) A、B 分別為 0.22、2.0。

表1 裂隙試件J1、J2及J3的曲折效應參數(shù)

把 A＝0.22、B＝2.0 代入式（9）得到考慮裂隙面粗糙程度和滲流曲折效應的裂隙滲流計算公式為

在式（15）中包含兩類參數(shù)，τ和 ＜b＞為根據(jù)裂隙三維空腔組合形貌計算得到的參數(shù)；Ts和Δ為根據(jù)裂隙三維表面形貌計算得到的參數(shù)。只要根據(jù)裂隙的三維空腔組合形貌和三維表面形貌數(shù)據(jù)確定了上述參數(shù)，就可以根據(jù)式（15）計算裂隙在設定滲流邊界條件下的滲流體積流量。

4.2.2 非線性分形模型與其他模型的對比

本文根據(jù)裂隙試件J1、J2及J3的形貌數(shù)據(jù)和數(shù)值計算結果求解了裂隙滲流計算公式中的待定常數(shù)A和B，得到了考慮裂隙面粗糙程度和滲流曲折效應的裂隙滲流計算公式。為了驗證該公式的正確性，下面將根據(jù)裂隙試件J4、J5的形貌數(shù)據(jù)，采用式（15）計算裂隙試件J4、J5的滲流體積流量，并將其與滲流數(shù)值計算結果進行比較。根據(jù)裂隙試件J4、J5的三維表面形貌和三維空腔組合形貌計算得到的裂隙滲流曲折因子、裂隙面曲折系數(shù)和相對粗糙度見表2。把表2中的數(shù)據(jù)代入式（15）中就可以計算得到試件J4、J5的滲流體積流量預測值，將其與Reynolds方程下的裂隙滲流數(shù)值計算結果進行比較，如圖4所示。

表2 裂隙試件J4、J5的曲折效應參數(shù)

圖4 試件J4、J5的裂隙滲流數(shù)值計算結果

由圖4可知，采用式（15）計算的裂隙滲流體積流量隨平均水力梯度的變化趨勢與數(shù)值計算結果基本一致，且兩者的計算結果吻合較好，可見采用裂隙試件J1、J2及J3數(shù)值計算結果擬合得到的待定常數(shù)A和B具有很好的可信度，由此得到的裂隙滲流計算公式可以用于計算裂隙的滲流體積流量，能夠較好地反映裂隙中的真實滲流情況。

5 結 論

（1）根據(jù)Reynolds控制方程下的裂隙滲流值數(shù)計算結果，提出了定量描述裂隙滲流曲折效應的平均曲折因子τ的定義，并根據(jù)裂隙的三維空腔組合形貌提出了平均曲折因子的計算方法。

（2）在總結現(xiàn)有考慮裂隙表面粗糙性的立方定理修正公式的基礎上，根據(jù)巖石滲透性等效溝槽模型，在滲流計算公式中引入裂隙滲流曲折因子τ表征裂隙滲流曲折效應，得到了綜合考慮裂隙表面粗糙性和裂隙滲流曲折效應的滲流計算新公式。

（3）將裂隙滲流計算新公式計算結果與Reynolds方程下的裂隙滲流數(shù)值計算結果進行對比，結果表明，考慮裂隙表面粗糙性和裂隙滲流曲折效應的滲流計算新公式計算結果與Reynolds方程下的裂隙滲流數(shù)值計算結果吻合較好，能真實地反映裂隙中的滲流情況。同時，滲流計算新公式具有計算簡便、參數(shù)物理意義明確、計算結果可靠的特點，在工程中有較好的應用前景。