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

殷德勝，汪衛(wèi)明，陳勝宏

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

1 引 言

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

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

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

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

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

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

2 塊體單元法原理簡介

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

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

有關(guān)塊體單元法的應(yīng)力、滲流分析的詳細過程請參閱文獻[8]。

3 耦合機制

巖體的變形大部分由結(jié)構(gòu)面產(chǎn)生，結(jié)構(gòu)面隙寬的微小改變將引起滲流量的重大變化。同時，滲流場的變化會引起水荷載的變化，從而導(dǎo)致應(yīng)力場的重分布。因而，裂隙巖體的結(jié)構(gòu)場和滲流場之間存在著重要的耦合關(guān)系。

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

3.1 結(jié)構(gòu)面變形對裂隙剛度系數(shù)的影響

結(jié)構(gòu)面的相對變形影響著剛度系數(shù)，從而間接地影響應(yīng)力場的分布。陳勝宏等[10]曾提出有關(guān)節(jié)理面變形和應(yīng)力關(guān)系的“充填模型”，將節(jié)理裂隙視為一種均勻介質(zhì)，并據(jù)此推導(dǎo)了有關(guān)節(jié)理面法向、切向剛度系數(shù)與應(yīng)力之間簡單而實用的關(guān)系式。

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

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

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

“充填模型”假設(shè)裂隙面上，“充填介質(zhì)”是各向同性體。據(jù)此，可得到裂隙面局部坐標(biāo)下的彈性本構(gòu)關(guān)系為

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

于是式(2)可以變成：

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

由于開度e與法向變形uz相關(guān)，而uz又與應(yīng)力有關(guān)聯(lián)，故 kn、 ks都是應(yīng)力的函數(shù)。

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

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

將式(7)積分后得到：

把式(9)代入式(5)，即可得到法向、切向剛度系數(shù)與法向應(yīng)力的關(guān)系：

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

式(11)表明了裂隙剛度系數(shù)與法向應(yīng)力之間存在著指數(shù)關(guān)系。

3.2 結(jié)構(gòu)面變形對導(dǎo)水系數(shù)的影響

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

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

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

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

3.3 耦合步驟

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

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

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

4 算例驗證

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

對于均質(zhì)材料，根據(jù)文[16]的建議，虛擬結(jié)構(gòu)面的法向剛度、切向剛度按式(13)選?。?/p>

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

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

圖2 算例塊體系統(tǒng)Fig.2 Block element system of the example

表1 裂隙面的參數(shù)Table 1 Parameters of fractures

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

圖3 結(jié)構(gòu)面上的水頭分布Fig.3 Hydraulic distribution on the structure

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

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

5 工程算例

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

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

圖6 塊體系統(tǒng)Fig.6 Block element system

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

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

表2 結(jié)構(gòu)面參數(shù)Table 2 Parameters of discontinuities

表3 裂隙參數(shù)Table 3 Parameters of fractures

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

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

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

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

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

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

6 結(jié) 論

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

（2）法向應(yīng)力對裂隙的滲透特性影響很大。在高壓應(yīng)力的區(qū)域，導(dǎo)水系數(shù)減小，滲透坡降增加，流速增大；其他區(qū)域流速減小。從算例中可以看出，重力壩分析中若不考慮耦合作用，則可能導(dǎo)致?lián)P壓力的計算值偏小。對實際工程會產(chǎn)生不利影響。

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

[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]陳平, 張有天. 裂隙巖體滲流與應(yīng)力耦合分析[J]. 巖石力學(xué)與工程學(xué)報, 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]王媛, 徐志英, 速寶玉. 裂隙巖體滲流與應(yīng)力耦合分析的四自由度全耦合法[J]. 水利學(xué)報, 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]張有天. 巖石水力學(xué)與工程[M]. 北京: 中國水利水電出版社, 2005.

[6]柴軍瑞, 仵彥卿. 巖體滲流場與應(yīng)力場耦合分析的多重裂隙網(wǎng)絡(luò)模型[J]. 巖石力學(xué)與工程學(xué)報, 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]張國新, 武曉峰. 裂隙滲流對巖石邊坡穩(wěn)定的影響——滲流、變形耦合作用的DDA法[J]. 巖石力學(xué)與工程學(xué)報, 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]陳勝宏. 計算巖體力學(xué)與工程[M]. 中國水利水電出版社, 2006.

[9]殷德勝, 汪衛(wèi)明, 陳勝宏. 三維隨機裂隙巖體滲流分析的塊體單元法[J]. 巖土力學(xué), 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]陳勝宏, 王鴻儒, 熊文林. 節(jié)理面滲流性質(zhì)的探討[J].武漢水利電力學(xué)院學(xué)報, 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]王媛. 單裂隙面滲流與應(yīng)力的耦合特性[J]. 巖石力學(xué)與工程學(xué)報, 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]薛孌鸞, 陳勝宏. 巖石裂隙滲流與法向應(yīng)力耦合的復(fù)合單元模型[J]. 巖石力學(xué)與工程學(xué)報, 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]卓家壽, 趙寧. 不連續(xù)介質(zhì)靜、動力分析的剛體-彈簧元法[J]. 河海大學(xué)學(xué)報, 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]. 工程力學(xué), 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]汪衛(wèi)明, 徐明毅, 陳勝宏. 復(fù)雜邊界條件下的巖體網(wǎng)絡(luò)滲流分析[J]. 巖石力學(xué)與工程學(xué)報, 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.