裂隙巖體滲流-應力耦合等效滲流阻模型

陶 煜，劉衛群

（1.中國礦業大學 力學與建筑工程建學院，江蘇 徐州 221116；2.中國礦業大學 深部巖土力學與地下工程國家重點實驗室，江蘇 徐州 221116）

1 引 言

在巖土工程實踐中，經常會遇到水文地質條件較為復雜的環境，受到多類水體運移及裂隙滲流作用的影響[1]。為了考察區域應力作用下裂隙巖體中復雜的滲流狀態，需要對其內部裂隙發育以及這些裂隙在應力作用下的滲流規律進行分析，一些國內學者已對此進行了相關的研究[2－4]。大量關于采動巖體滲流的試驗和實測結果表明，其滲透特性隨應力場改變而發生顯著變化。由于應力場變化導致區域裂隙系統的發育程度和組織形態的改變，因此，滲透率的改變與裂隙系統的發育程度和組織形態變化有著密切的聯系。

為此，區域裂隙系統的發育程度和組織形態成為本文研究工作的切入點。圖1是選自文獻[2]的關于不同采厚工作面推進200 m采動裂隙演化過程及形態特征圖。

通過觀察不難發現，裂隙巖體主干裂隙的演化有一個明顯的特征。就是在一段區域內，主干裂隙是沿著一定的方向發育的，并且存在分叉的現象。由于主干裂隙在區域裂隙系統中分布較為明顯，其貫通程度、裂隙張開度、跨越區域都超過非主干裂隙，所以在裂隙滲流中起著主導作用。因此，考慮區域主干裂隙應力影響下的滲流規律對分析裂隙巖體滲流特性具有重要意義。

為了考察貫通區域主干裂隙系統的滲透能力，本文運用平行板滲流立方定律，結合彈性力學方法和模擬電路知識，建立了一種基于等效滲流阻的滲流-應力耦合模型。通過該模型可以進一步地分析這種貫通主干裂隙系統受應力作用的滲流演化規律。

圖1 不同采厚工作面推進200 m采動裂隙演化過程及形態特征[2]Fig.1 Evolutions and configuration characteristics of MIF in caving 200 m with different mining thicknesses [2]

2 裂隙巖體平行板模型

對厚度為a、隙寬為b的平行板裂隙均勻流，運用平行板滲流立方定律有

式中：q為流量；P為水壓；μ為水的動力黏滯系數，其值與溫度有關。

假若裂隙方向與流速為v，截面積為S的水流方向成α角度，則

式（2）是單裂隙平行板模型隙寬與單寬流量之間的關系[5]。巖石中實際的壁面是粗糙的，若將由光滑平行板縫隙層流流態導出的立方定律用于實際裂隙，則需要對立方定理進行修正。實際裂隙面起伏不平，有許多凸體，致使水流在同一水力梯度下，同一隙寬的粗糙裂隙比光滑裂隙通過的流量要小。Ломизе[6]進行了粗糙裂隙水力試驗，發現裂隙面的突起度 Δ與隙寬的比值對裂隙過流能力有很大影響。他將立方定理加入了一個粗糙度修正系數C，

很多學者根據試驗或理論提出了各自的粗糙度修正公式，主要有 Ломизе（1951）、Louis（1967），Neuzil(1981)、Tsang（1981）、Barton（1985）、周創兵（1996）等。王媛和速寶玉[7]已就該方面的文獻進行了綜述。

取垂直于x軸的裂隙截面寬度為bv，則裂隙寬b=bvcosα ，則式（3）改寫為

3 滲流等效流阻模型的建立

天然裂隙巖體中的裂隙分布多種多樣，一般可分為較為離散細微的裂隙網絡和較大尺度的主干裂隙系統。前者參與滲流過程程度相對較弱，且運用連續介質理論方法處理較為成熟。而后者在貫通程度、張開度、跨越區域都優于前者，參與滲流作用的比重遠大于前者。作為滲流作用主體，主干裂隙系統成為本文滲流分析的主要對象。

圖2選自Nick Barton的巖體節理（Rock Joints）一書[8]，展現了不同巖體的裂隙分布?？梢暂^為明顯地觀察出主干裂隙系統與離散的裂隙網絡。在裂隙巖體中，主干裂隙系統是由有限個單裂隙組合而成。而每一個單裂隙之間僅有如圖3所示的2種基本組合關系。

本文對這兩種組合關系下的滲流特性分別進行討論。

3.1 串聯關系

如圖4所示，根據b1裂隙的流量和b2裂隙的流量相等的特點，有

式中：qf為裂隙流量；Δp為水的壓力差；α1、α2分別為上下部分裂隙的厚度；b1、b2分別為上下部分裂隙的寬度；Δx1、Δx2為上下部分裂隙的長度；α1、α2分別為上下裂隙分別與流動方向之間的夾角。在圖4中α1、α2都為0°。

而 Δp=Δp1+Δp2（如圖4所示），所以

圖2 幾種巖體裂隙圖[8]Fig.2 Fractures in several kinds of rock masses[8]

圖3 裂隙基本組合關系圖Fig.3 Relationship of fracture combinations

圖4 串聯情形平行板裂隙滲流系統Fig.4 Fracture seepage system in series

而在模擬電路中，串聯公式為 U=IR1+IR2，由式（6）可以看出在裂隙滲流中存在類似電路串聯關系。

3.2 并聯關系

兩裂隙并聯，如圖5所示，根據裂隙系統總流量等于b1裂隙流量和b2裂隙流量之和，可以得出

式中：α1、α2分別為左右兩個裂隙分別與流動方向之間的夾角。在圖5中α1、α2都為0°。

圖5 并聯情形平行板裂隙滲流系統Fig.5 Fracture seepage system in parallel

所以，根據式（7）我們可以得到在裂隙滲流中存在類似電路并聯的關系。

因此，可以把區域主干裂隙中的流體滲流等效為模擬電路中的串并聯電路群，兩者之間的相似量如表1中給出。

表1 模擬電路問題與滲流問題中的相似量表Table1 Similar capacities between electronics and seepage problems

假設一個區域主干裂隙系統有i個層次，即i層串聯等效滲流阻，每個層次有j個并聯等效滲流阻，那么此區域主干裂隙系統內的等效滲流阻Rf為

式中：bij、Cij、α分別為裂隙系統第i層串聯部分里的第j個并聯裂隙隙寬、裂隙粗糙度修正系數和裂隙與水流方向夾角。

裂隙系統流量為

由于裂隙系統的滲透率k與等效滲流阻之間的關系為

式中：A為滲流截面積。

故此裂隙系統的滲透率為

本文稱這種由等效滲流阻構成的串并聯滲流模型為“等效滲流阻模型”。

4 裂隙巖體應力與滲透率之間的關系

巖體中的應力變化導致裂隙寬度的變化，進而導致巖體滲透性的改變。為了探究裂隙巖體區域裂隙系統應力與滲透率之間的關系，假設在同一裂隙網絡中同一層面上分布著n條裂隙，如圖6所示。

圖6 三向應力條件下裂隙巖體應力-滲透性計算模型Fig.6 3D stress-permeability model for fractured rock mass

沿 x方向的巖體總位移Δutx等于裂隙的位移Δ ufx與巖塊位移 Δurx之和，所以

位移用應變表示有

式中：Sx、bx分別為巖體的寬度和多條裂隙沿x 方向的寬度和。

沿x方向巖體應變εtx又可表示為

式中：Emx為裂隙巖體在x方向上的彈性模量；γ為泊松比；σx、σy、σz分別為x、y、z方向上的應力。

沿x方向巖塊應變εrx可表示為

式中：Er為巖塊的彈性模量。

由此，可得到多條裂隙的隙寬變化總量

如果這j條裂隙位于第i層串聯中，那么式（19）可以改為

還可以進一步得到應力影響下該裂隙系統的滲透流量為

以及滲透率與應力之間的關系為

以上為本文建立的基于等效滲流阻的應力耦合滲流模型，可以用于分析應力作用下含區域主干裂隙巖體滲流演化規律。

5 算 例

結合上文建立的應力耦合等效滲流阻模型，選擇含典型的“人”字形組合裂隙系統的 50 mm×50 mm×100 mm立方體巖石試樣進行滲流分析，通過 Matlab計算程序考察該組合裂隙系統試樣的滲流特性，分析對象如圖7所示。

圖7 50 mm×50 mm×100 mm立方體裂隙巖樣圖Fig.7 50 mm×50 mm×100 mm cubic fractured rock sample

表2 全局參數設定表（以國際單位制計）Table2 Global parameters (using SI units)

表3 各條裂隙參數賦值表（以國際單位制計）Table3 Parameters of each fracture (using SI units)

在眾多粗糙度修正公式中，由于裂隙面粗糙度修正系數修正法相對于隙寬密度分布修正法等其他方法而言較為簡便，并且其使用限制為裂隙面必須高度分布比較均勻，而本算例符合其使用條件。所以在計算中運用了Louis根據試驗總結出的粗糙度修正公式

設巖石中水溫為15 ℃，則 μ=1.14×10-3N? s/m2。當Er=1 GPa，γ=0.2，μ=1.14×10-3N? s/m2（具體各類參數如表2、3所示），且巖體中僅存在裂隙法向應力σx時，通過本文建立的裂隙等效滲流阻模型，可得出應力與等效滲流阻的關系曲線，如圖 8所示。滲流阻與應力關系曲線表明，隨著裂隙面法向應力的增加，裂隙系統等效滲流阻不斷增加，兩者呈冪函數關系遞增。也可以得到滲透率與應力之間的關系曲線，如圖9所示。圖中滲透率與應力關系曲線顯示，隨著裂隙法向應力的增加，裂隙系統的滲透率不斷減少，呈冪函數關系遞減。

圖8 滲流阻與應力的關系Fig.8 Seepage resistance-stress relationship

圖9 滲透率與應力的關系Fig.9 Permeability-stress relationship

6 模型的擴展性探討

等效滲流阻模型具有一定的可擴展性，有很多可以繼續開展工作的地方，例如將等效滲流阻模型與有限元相結合，可以將裂隙系統劃分為若干子裂隙系統，如圖10所示。其中A～F為子裂隙系統，而 P1～ P9為聯系這些子系統的裂隙出入口水壓。

各子裂隙系統在應力場中所受的應力不同，而在子裂隙系統內可以將應力簡化為相同。故子裂隙系統的尺寸與應力場的應力梯度分布相關，應力變化緩慢的地方可以選擇相對較大的子裂隙系統尺寸，反之則應該選擇較小的尺寸。經過進一步細分后，可將等效滲流阻模型與有限元分析方法結合起來，對大范圍含裂隙巖體滲流區域應力場中主干裂隙系統的滲流演化規律進行分析。

圖10 子裂隙系統滲流圖Fig.10 Seepage in subsystem of fractures

7 結 論

等效滲流阻模型具有以下優點：

（1）由于采用現實裂隙構造，無需求解復雜的偏微分方程，較隨機裂隙網絡等其他模型更為實用。

（2）在一個區域裂隙系統內我們定義了“滲流阻”的概念，不去分個考慮每個裂隙的應力狀態下的滲流規律，化繁為簡，利用區域滲流阻整體考慮它在整體應力影響下的滲流規律，分析問題更為方便有效。

（3）模型基于平行板滲流理論、模擬電路、彈性力學知識，具有普適性和可擴展性強的優點，更客觀地描述應力變化下的區域主干裂隙系統的滲流規律。

由于并聯裂隙情形采用等間距平行板模型簡化，當分叉裂隙走向不同時存在一定誤差，但沿高應力方向裂隙取向相近實際區域主干裂網中非常普遍，基本能滿足模型要求。

文中模型經進一步細化后，可為巖體滲流分析問題提供理論支持。

[1]ZHANG J,SHEN B. Coal mining under aquifers in China:A case study[J]. International Journal of Rock Mechanics and Mining Sciences,2004,41(4): 629－639.

[2]楊科,謝廣祥. 采動裂隙分布及其演化特征的采厚效應[J]. 煤炭學報,2008,33(10): 1092－1096.YANG Ke,XIE Guang-xiang. Caving thickness effects on distribution and evolution characteristics of mining induced fracture[J]. Journal of China Coal Society,2008,33(10): 1092－1096.

[3]陳占清,繆協興,劉衛群. 采動圍巖中參變滲流系統的穩定性分析[J]. 中南大學學報(自然科學版),2004,35(1): 129－132.CHEN Zhan-qing,MIAO Xie-xing,LIU Wei-qun.Analysis of stability of parametric system of seepage flow in wall rock affected by mining[J]. Journal of Central South University (Natural Science),2004,35(1): 129－132.

[4]張永波,靳鐘銘,劉秀英. 采動巖體裂隙分形相關規律的實驗研究[J]. 巖石力學與工程學報,2004,23(20):3426－3429.ZHANG Yong-bo,JIN Zhong-ming,LIU Xiu-ying.Testing study of fractal correlation law of cracks in mined rock mass[J]. Chinese Journal of Rock Mechanics and Engineering,2004,23(20): 3426－3429.

[5]周志芳. 裂隙介質水動力學原理[M]. 北京: 高等教育出版社,2007.

[6]Ломизе М. Филътрация в трещиноватых породах[M].Москва: Госэнергоиздат,. 1951.

[7]王媛,速寶玉. 單裂隙面滲流特性及等效水力隙寬[J].水科學進展,2002,13(1): 61－67.WANG Yuan,SU Bao-yu. Research on the behavior of fluid flow in a single fracture and its equivalent hydraulic aperture[J]. Advances in Water Science,2002,13(1): 61－67.

[8]BARTON N. Rock Joints[M]. Rotterdam: A. A. Balkema Publishers,1990.