隨機裂隙巖體滲流-應力耦合分析的塊體單元法

殷德勝，汪衛明，陳勝宏

（武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072）

1 引 言

1959年法國Malpasset拱壩的潰壩事故，導致裂隙巖體滲流與應力的耦合問題成為巖石水力學的熱門課題。由于該問題的邊界條件和材料性質的復雜性，很難進行物理模型試驗，人們普遍采用數值解法。

有限單元法是連續介質力學的代表，它有著完備的理論體系，在巖體滲流-應力耦合研究中的應用最為廣泛。Noorishad等[1]利用土體固結分析理論，并運用有限元法中的節理單元，研究了非連續介質中的固液兩相介質的耦合間題；Oda[2]將裂隙巖體滲透張量與彈性張量統一用裂隙網絡幾何張量表示，建立了較為嚴密的裂隙巖體滲流特性與應力特性之間的耦合關系；在國內，陳平和張有天[3]考慮了裂隙變形的非線性特性，以裂隙平均隙寬為計算參數，對立方定理作了修正，提出了滲流-應力耦合計算方法；王媛等[4]則建立的同時以滲流水壓力和位移為未知值的4自由度全耦合法，這降低了耦合分析的難度。由于已有多種成熟的商業軟件如ABAQUS、ADINA等能夠模擬多場耦合問題，這將進一步加強有限元法在巖土工程滲流-應力耦合這一領域中的應用。

對于復雜的工程，已有的研究成果大多基于二維情況，由于三維裂隙巖體的幾何復雜性，目前尚無一個強大的前處理器能夠輕松地對其進行網格離散，這是限制有限元法在該領域進一步應用的主要原因。學者們只能退而求其次，大多采用等效連續介質模型進行耦合研究。然而，裂隙中的水主要通過裂隙等結構面流動，其滲流與傳統的多孔介質滲流存在著本質上的差別。對裂隙巖體進行滲流分析，宜采用裂隙網絡模型或雙重介質模型[5]。

雙重介質模型除裂隙網格外，將巖塊視為滲透系數較小的連續滲透介質，且能考慮巖塊孔隙與巖體裂隙之間水量的交換，是更接近于實際的模型。柴軍瑞等[6]將各種裂隙和孔隙按規模和滲透性分為4級，建立起巖體滲流場與應力場耦合分析的多重裂隙網絡模型。綜合來說，同連續介質方法相比，雙重介質模型優勢明顯，但由于繼承了網格剖分的困難，實現難度相對更大。綜合起來，不連續介質力學分析方法結合裂隙網絡滲流模型則是個相對合理可行的思路。

不連續變形分析方法（discontinuous deformation analysis method，DDA）屬于不連續介質力學范疇，十分適合處理裂隙巖體滲流、應力的計算問題。張國新等[7]應用DDA法，對日本某事故隧道進行了二維滲流-應力耦合分析，得到了邊坡失穩坍塌的機制。應用DDA法進行三維滲流-應力的耦合分析，在計算中尚存在接觸難以準確確定、罰函數如何施加等諸多困難，目前的研究成果很少。

塊體單元法也屬于不連續介質力學范疇，是一種能對節理裂隙巖體進行滲流、應力分析的數值方法[8]。在塊體單元法的基礎上，本文考慮了節理裂隙的剛度和導水系數隨開度的變化關系，并將常規的塊體單元法結構計算模塊和滲流分析模塊整合成一個系統，應用兩場交叉迭代技術，完成了隨機裂隙巖體的滲流-應力耦合研究。該方法前、后處理簡單，計算過程也不復雜，能有效地應用于三維隨機裂隙巖體的耦合分析中。

2 塊體單元法原理簡介

巖體的變形主要都集中于節理、裂隙等結構面上，相比而言，塊體的變形可以忽略不計。在此假設條件下，塊體單元法認為塊體為剛體，只考慮結構面的變形和強度特性，以塊體單元形心位移為基本未知量，通過建立塊體單元的平衡方程、結構面變形和塊體單元位移的幾何相容方程、結構面的本構方程，綜合推導出塊體系統的整體平衡方程。

同時，斷層、夾層、節理和裂隙等結構面不僅是巖體中的軟弱面，且由于其滲透系數通常遠大于完整巖石，故它們也是主要的滲流通道。假定巖塊不透水，即在巖塊中滲透系數k =0，滲流只發生在裂隙面上，根據裂隙網絡交叉處的流量守恒和分片變分原理，建立了與塊體單元法相配套的隨機裂隙巖體的滲流分析方法[9]。

有關塊體單元法的應力、滲流分析的詳細過程請參閱文獻[8]。

3 耦合機制

巖體的變形大部分由結構面產生，結構面隙寬的微小改變將引起滲流量的重大變化。同時，滲流場的變化會引起水荷載的變化，從而導致應力場的重分布。因而，裂隙巖體的結構場和滲流場之間存在著重要的耦合關系。

然而，在應力場和滲流場相互影響的定量描述上，目前尚無定論。由于裂隙面上的法向應力引起裂隙變形相對顯著，為抓住主要矛盾，本文僅研究法向應力對裂隙變形特性、滲透特性的影響。

3.1 結構面變形對裂隙剛度系數的影響

結構面的相對變形影響著剛度系數，從而間接地影響應力場的分布。陳勝宏等[10]曾提出有關節理面變形和應力關系的“充填模型”，將節理裂隙視為一種均勻介質，并據此推導了有關節理面法向、切向剛度系數與應力之間簡單而實用的關系式。

以單裂隙巖體說明裂隙滲流與應力的耦合關系。設圖1中的裂隙面受法向應力σz和切向應力τzx,τzy作用，由于裂隙面的開度較小，故應變滿足：

式中： εx、 εy分別為x、y向的正應變， γxy、 γyx為xy平面上的剪應變。

圖1 單裂隙巖體Fig.1 Rock mass containing one fracture

“充填模型”假設裂隙面上，“充填介質”是各向同性體。據此，可得到裂隙面局部坐標下的彈性本構關系為

式中：λ、G為裂隙內充填物的拉梅常數。記e0為裂隙面應力為0時的開度，ux、 uy、 uz分別為x、y、z三向的相對位移，則裂隙開度滿足：

于是式(2)可以變成：

式中： kn、 ks分別為裂隙面的法向、切向的剛度系數，可由充填物的彈性參數和開度確定：

由于開度e與法向變形uz相關，而uz又與應力有關聯，故 kn、 ks都是應力的函數。

從式(4)中取第1行，可以得到：

結合公式(3)，上式可寫成：

將式(7)積分后得到：

把式(9)代入式(5)，即可得到法向、切向剛度系數與法向應力的關系：

記 kn0= λ+ 2G /e0,ks0= G/e0，用它來表示初始剛度，即裂隙在無應力狀態時的剛度系數，則式(10)可寫成：

式(11)表明了裂隙剛度系數與法向應力之間存在著指數關系。

3.2 結構面變形對導水系數的影響

立方定律是研究裂隙巖體滲流的基礎。然而，天然節理并非完全光滑平直，也不是純粹的孔隙，因而立方定律的適用條件有限。近幾十年來，許多考慮了節理面的粗糙性和凸起物高度分布情況的修正立方定律[11－12]相繼被提出，其主要思想就是用等效水力隙寬代替幾何平均隙寬。然而，等效水力隙寬和幾何平均隙寬是完全不同的兩個概念[13]，若將其直接用于耦合研究，不但理論說不通，且流量與隙寬的3次方成正比，會導致耦合分析中計算過于敏感和收斂困難。故立方定律在耦合研究中直接應用較少。

“充填模型” 同樣適用于結構面的滲流特性的研究。該模型中假設裂隙中填充物的滲透系數不變，而導水系數隨著開度的變化而變化，并結合式(9)，得到：

式中：kd為導水系數；kd0為初始導水系數，即裂隙面無應力時的導水系數，可通過試驗獲得，其他同式(10)。

式(11)、(12)能較好地模擬裂隙剛度系數和導水系數與應力之間的關系。本文運用這套關系式，結合塊體單元法進行了裂隙巖體滲流-應力的耦合研究。

3.3 耦合步驟

根據結構面變形對應力場、滲流場的影響關系式，使用Fortran語言，編制了塊體元滲流-應力的耦合分析程序。該程序中，充分考慮塊體元應力、滲流分析的共同點，將2個分析模塊整合在一個程序中，顯著減少了公共計算部分的計算量，且相比傳統的兩場迭代，不需要2個程序相互調用時讀寫0時文件的時間開銷，且降低了文件傳遞導致的累積誤差，耦合計算的效率和精度均較高。

程序中，耦合分析的具體實施步驟如下：

(1)由初始剛度系數和初始導水系數，結合邊界條件，計算得到初始應力場 {σ}(0)和初始滲流場{φ}(0)；(2)結合式(11)、(12)，根據{σ}(0)計算裂隙面新的法向、切向剛度系數和導水系數，并重新生成結構計算的剛度矩陣和滲流計算的導水矩陣；(3)由新的導水矩陣和滲流邊界條件，重新計算滲流場{φ}(i)；(4)由新的滲流場 {φ}(i)計算出滲透荷載，并加上其它荷載，重新計算應力場 {σ}(i)；(5)計算應力、水頭的增量： {Δ σ} ={σ}(i)-{σ}(i-1)和{Δφ}={φ}(i)-{φ}(i-1)，并取{Δσ}和{Δφ}的2范數。若和小于給定的誤差值，則迭代結束，否則重新執行步驟(2)，直至滲流場和應力場都滿足收斂精度為止。

4 算例驗證

薛孌鸞[14]曾用基于雙重介質模型的復合單元法研究過裂隙巖體的滲流-應力耦合問題。本文采用她文中的算例，從規律上對本文方法進行驗證。如圖2所示的長方體巖塊試件，坐標系在左下角，試件的x方向長為1.0 m，y方向寬為0.1 m，z方向高為0.2 m。在z =0.1 m處有一條初始開度為1 mm的水平貫穿性裂隙。由于塊體單元法假設塊體為剛體，本文采用卓家壽等提出的界面元法[15－16]，用自定義的虛擬結構面將試件等分成若干塊體，將變形累積在結構面上，以此來近似模擬。

對于均質材料，根據文[16]的建議，虛擬結構面的法向剛度、切向剛度按式(13)選取：

式中： h1、 h2為塊體形心到結構面的距離；E為彈性模量；μ為泊松比。如圖2所示，總共引入9條虛擬結構面，將試件分成20個塊體。另外，虛擬結構面的導水系數的取值為實際裂隙導水系數的1%。所有結構面的參數見表1：

取巖塊試件的左側水位為 10 m，右側水位為5 m，其它外表面為不透水邊界。考慮水壓力并在上表面正中間的 23、53兩點處各施加集中力 F =0.1 MN，不考慮試件的自重。在除上表面外的其余5個側面均施加法向約束。

圖2 算例塊體系統Fig.2 Block element system of the example

表1 裂隙面的參數Table 1 Parameters of fractures

整理耦合和不耦合兩種工況下的成果，圖3為結構面上的水頭分布，圖4為不耦合和耦合兩種工況下y =0.05 m斷面上的流速矢量圖。

圖3 結構面上的水頭分布Fig.3 Hydraulic distribution on the structure

圖4 流速矢量圖Fig.4 Hydraulic velocity vector diagram

從圖3中可以看出，不考慮耦合作用時，結構面上的水頭呈線性分布；考慮耦合作用后，集中力作用點的左側部位，結構面的水頭趨于上游水位，右側部位的水頭則趨于下游水位。對比圖4可以看出：當不考慮耦合時，整個裂隙面上的流速是相同的；考慮耦合后，集中力作用點底部壓應力較大區域的導水系數減小，滲透坡降增加，流速增加，其他區域的流速減小。結果顯示了本文耦合機制和程序的合理性。

5 工程算例

圖5是貴州省烏江沙沱水電站非溢流壩段的典型剖面的示意圖。其巖基上含有F88斷層，j4、j5、j6三條夾層，并有兩組傾角約為45°的正交裂隙。沿壩軸線方向取寬度為1 m的壩段參與計算，對于斷層和夾層，采用確定性方式模擬；對于裂隙，采用蒙特卡羅法隨機生成，并用塊體識別方法[9]找出所有塊體。塊體系統如圖6所示。

圖5 重力壩剖面Fig.5 Profile diagram of gravity dam

圖6 塊體系統Fig.6 Block element system

對于地基部分，斷層、夾層、裂隙等結構面的參數可以通過實驗的方式獲得；對于壩體部分，仍按卓家壽的方法，按5 m的間距獲得虛擬結構面的參數。帷幕采用文獻[17]的方法，即將帷幕塊體周邊結構面視為相對不透水面，其導水系數取為較小值，本文取值為1×10-8m/s。

斷層、夾層的材料參數見表 2，裂隙參數見表3，壩體和基巖的參數見表4。

表2 結構面參數Table 2 Parameters of discontinuities

表3 裂隙參數Table 3 Parameters of fractures

表4 基巖和壩體混凝土參數Table 4 Parameters of rock foundation and dam concrete

同樣，計算也分為不耦合和耦合兩種工況。上下游水位見圖5，考慮的荷載有壩體和基巖的自重、靜水壓力和滲透荷載作用。

整理了兩種工況的計算成果對比。圖7為裂隙網絡的水頭等勢線，圖8為水力坡降矢量圖。

圖7 裂隙網絡等水頭線（單位：m）Fig.7 Hydraulic potential isolines of fracture network (unit: m)

圖8 水力坡降矢量圖Fig.8 Hydraulic gradient vector diagrms

從圖 7、8中可以看出：考慮耦合作用以后，壩前地基的水力坡降逐漸減小，水頭普遍增大；壩后基礎的變化規律與壩前類似；壩體下方地基的水平向主導的裂隙上的壓應力增加，導水系數減小，坡降、流速均有所增加，進而導致建基面的揚壓力增加，對壩體的穩定不利。

6 結 論

（1）塊體單元法本身即是為解決節理、裂隙巖體的應力、滲流問題而提出，將其應用于裂隙巖體的滲流-應力耦合中，有著先天的優勢。而且能考慮三維巖體中含有的復雜隨機裂隙網絡，這是其他模型所不能考慮或很難考慮的；

（2）法向應力對裂隙的滲透特性影響很大。在高壓應力的區域，導水系數減小，滲透坡降增加，流速增大；其他區域流速減小。從算例中可以看出，重力壩分析中若不考慮耦合作用，則可能導致揚壓力的計算值偏小。對實際工程會產生不利影響。

然而，就目前的計算機發展水平，采用塊體單元法進行滲流-應力耦合分析的計算量仍顯偏大，需要進一步優化數據結構、改進積分算法，甚至采用分布式并行計算，使該方法真正意義上用于工程的三維分析中。

[1]NOORISHAD J, AYATOLLAHI M S, WITHERSPOON P A. A finite-element method for coupled stress and fluid flow analysis in fractured rock masses[J]. International Journal of Rock Mechanics and Mining Sciences, 1982:185－193.

[2]ODA M. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses[J].Water Resources Research, 1986, 22(13): 1845－1856.

[3]陳平, 張有天. 裂隙巖體滲流與應力耦合分析[J]. 巖石力學與工程學報, 1994, 13(4): 299－308.CHEN Ping, ZHANG You-tian. Coupling analysis of seepage/stress for joint rock[J]. Chinese Journal of Rock Mechanics and Engineering, 1994, 13(4): 299－308.

[4]王媛, 徐志英, 速寶玉. 裂隙巖體滲流與應力耦合分析的四自由度全耦合法[J]. 水利學報, 1998, (7): 55－59.WANG Yuan, XU Zhi-ying, SU Bao-yu. Four-freedom complete method for the seepage-stress coupled analysis in fissured rock masses[J]. Journal of Hydraulic Engineering, 1998, (7): 55－59.

[5]張有天. 巖石水力學與工程[M]. 北京: 中國水利水電出版社, 2005.

[6]柴軍瑞, 仵彥卿. 巖體滲流場與應力場耦合分析的多重裂隙網絡模型[J]. 巖石力學與工程學報, 2000, 19(6):712－717.CHAI Jun-rui, WU Yan-qing. Research on multiple-level fracture network model for coupled seepage and stress fields in rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 2000, 19(6): 712－717.

[7]張國新, 武曉峰. 裂隙滲流對巖石邊坡穩定的影響——滲流、變形耦合作用的DDA法[J]. 巖石力學與工程學報, 2003, 22(8): 1269－1275.ZHANG Guo-xin, WU Xiao-feng. Influence of seepage on the stability of rock slope coupling of seepage and deformation by DDA method[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(8): 1269－1275.

[8]陳勝宏. 計算巖體力學與工程[M]. 中國水利水電出版社, 2006.

[9]殷德勝, 汪衛明, 陳勝宏. 三維隨機裂隙巖體滲流分析的塊體單元法[J]. 巖土力學, 2009, 30(8): 1269－1275.YIN De-sheng, WANG Wei-ming, CHEN Sheng-hong.Block element method for seepage analysis in three dimensional random fracture network[J]. Rock and Soil Mechanics, 2009, 30(8): 1269－1275.

[10]陳勝宏, 王鴻儒, 熊文林. 節理面滲流性質的探討[J].武漢水利電力學院學報, 1989, 22(1): 53－60.CHEN Sheng-hong, WANG Hong-ru, XIONG Wen-lin.Study on the seepage characteristics of joint surface[J].Journal of Wuhan University of Hydraulic and Electric Engineering, 1989, 22(1): 53－60.

[11]LOMIZE G M. Flow in fractured rocks[M]. Moscow:Gesenergoizdat, 1951.

[12]LOUIS C. Rock mechanics[M]. New York: Elsevier Science, 1974.

[13]王媛. 單裂隙面滲流與應力的耦合特性[J]. 巖石力學與工程學報, 2002, 21(1): 83－87.WANG Yuan. Coupling characteristic of stress and fluid flow within a single fracture[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(1): 83－87.

[14]薛孌鸞, 陳勝宏. 巖石裂隙滲流與法向應力耦合的復合單元模型[J]. 巖石力學與工程學報, 2007, 26(增刊):2613－2619.XUE Luan-luan, CHEN Sheng-hong. Composite element model of seepage-normal stress coupling for rock fractures[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(Supp.): 2613－2619.

[15]卓家壽, 趙寧. 不連續介質靜、動力分析的剛體-彈簧元法[J]. 河海大學學報, 1993, 21(5): 34－43.ZHUO Jia-shou, ZHAO Ning. Piecewise rigid body-interface spring method for problems of discontinuous medium[J]. Journal of Hohai University,1993, 21(5): 34－43.

[16]方義琳, 卓家壽, 章青. 具有任意形狀單元離散模型的界面元法[J]. 工程力學, 1998, 15(2): 27－37.FANG Yi-lin, ZHUO Jia-shou, ZHANG Qing. Interface stress element method for distinct model with elements of arbitrary shape[J]. Engineering Mechanics, 1998, 15(2):27－37.

[17]汪衛明, 徐明毅, 陳勝宏. 復雜邊界條件下的巖體網絡滲流分析[J]. 巖石力學與工程學報, 2001, 20(4): 473－476.WANG Wei-ming, XU Ming-yi, CHEN Sheng-hong.Analysis method of network seepage for blocky rock mass with complicated boundary condition[J]. Chinese Journal of Rock Mechanics and Engineering, 2001,20(4): 473－476.