姚 池，姜清輝， ，葉祖洋，周創(chuàng)兵，

（1. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2. 武漢大學(xué) 土木建筑工程學(xué)院，武漢 430072）

1 引 言

當(dāng)前，描述裂隙巖體滲流的數(shù)學(xué)模型主要有3種方法，即等效連續(xù)介質(zhì)模型、離散介質(zhì)模型和雙重介質(zhì)模型。等效連續(xù)介質(zhì)模型將裂隙巖體等效為連續(xù)體[1－3]，假定地下水在巖體中的流動服從Darcy 定律，然后運用經(jīng)典的連續(xù)滲流理論進(jìn)行分析，目前應(yīng)用較為廣泛，但該模型不能很好地反映空間結(jié)構(gòu)面組合形成的裂隙網(wǎng)絡(luò)的導(dǎo)水作用及優(yōu)勢流效應(yīng)，在某些情況下數(shù)值模擬結(jié)果與實際出入很大。離散介質(zhì)模型假定巖塊不透水[4－6]，滲流只發(fā)生在巖體裂隙內(nèi)，地下水通過相互連通的裂隙網(wǎng)絡(luò)進(jìn)行傳導(dǎo)，巖體滲透特性受裂隙的跡長、走向、傾角和隙寬等參數(shù)控制，這種方法抓住了巖體中滲流主要受裂隙等結(jié)構(gòu)面控制這一本質(zhì)特征，具有廣泛的應(yīng)用前景，其主要缺點在于實際工程中裂隙網(wǎng)絡(luò)資料不易獲取，因此，常通過Monto-Carlo 法隨機(jī)生成具有相同統(tǒng)計規(guī)律的裂隙網(wǎng)絡(luò)來代替真實的裂隙網(wǎng)絡(luò)。雙重介質(zhì)模型同時考慮了裂隙的強導(dǎo)水性、巖塊的儲水作用及巖塊和裂隙的水量交換[7－8]，但該模型對裂隙系統(tǒng)的幾何形態(tài)和空間分布存在較嚴(yán)格的假設(shè)，從而限制了其應(yīng)用。

由于巖體本身賦存環(huán)境的復(fù)雜性，裂隙的發(fā)育極不規(guī)則，裂隙網(wǎng)絡(luò)的幾何構(gòu)成往往非常復(fù)雜，導(dǎo)致滲流在巖體中的運動顯示出明顯的非均勻性和各向異性，很多情況下并不存在所謂的 REV（representative elementary volume）。此時，如用等效連續(xù)介質(zhì)方法進(jìn)行滲流分析實際上是不合適的，應(yīng)該采用離散裂隙網(wǎng)絡(luò)的模擬方法?？傮w上，相對基于多孔連續(xù)介質(zhì)滲流分析方法，目前基于離散裂隙網(wǎng)絡(luò)模型滲流問題研究相對較少，尤其對有自由面的裂隙網(wǎng)絡(luò)滲流問題，國內(nèi)外只有少數(shù)的幾篇論文專門對其進(jìn)行了探討。Jing 等[9]利用DDA（discontinuous deformation analysis）形成的裂隙網(wǎng)絡(luò)，采用剩余流量法求解裂隙網(wǎng)絡(luò)自由面。王恩 志[10]借鑒初流量法的思想，提出了一種計算裂隙網(wǎng)絡(luò)滲透自由面的迭代算法。柴軍瑞等[11]提出含自由面巖體裂隙網(wǎng)絡(luò)非線性滲流分析的傳導(dǎo)矩陣調(diào)整法。由于上述方法主要借鑒求解連續(xù)介質(zhì)滲流自由面的一些常用的直覺化求解方法，這些方法不具有嚴(yán)格的理論的基礎(chǔ)，且因為沒有解決溢出點的奇異性，計算過程中往往會出現(xiàn)數(shù)值的不穩(wěn)定，對于含有大規(guī)模的復(fù)雜裂隙網(wǎng)絡(luò)則難以收斂。

變分不等式方法是近年來求解連續(xù)介質(zhì)無壓滲流分析理論上的一個重要進(jìn)展[12－13]，本文把單條裂隙概化為平行板模型，通過將Darcy 定律延拓到整個計算區(qū)域，定義了在整個區(qū)域上的非線性邊值問題，并將潛在溢出面邊界條件歸納為Signorini型邊界條件，建立與PDE 提法等價的變分不等式(variational inequality， VI)提法，從而在理論上解決了出滲點的定位問題。通過采用有限元方法的線單元模擬裂隙，在有限元計算格式中，引入連續(xù)的罰Heaviside 函數(shù)，顯著改善了離散裂隙網(wǎng)絡(luò)滲流自由面迭代過程中的數(shù)值不穩(wěn)定性。

2 有自由面裂隙網(wǎng)絡(luò)滲流問題的描述

對于一個任意的二維裂隙網(wǎng)絡(luò)，裂隙網(wǎng)絡(luò)的連通性及裂隙的跡長、水力開度、位置等參數(shù)決定了地下水滲流在巖體中的分布，由于孤立的裂隙、只與一條裂隙相交的枝狀裂隙及裂隙端部的獨頭不能作為滲流通道，對流場沒有貢獻(xiàn)，因此，在進(jìn)行滲流分析前需對隨機(jī)生成的裂隙進(jìn)行必要的刪減和截斷，經(jīng)過凈化處理形成連通的裂隙網(wǎng)絡(luò)，凈化處理前后裂隙網(wǎng)絡(luò)的對比可參見圖1。

不失一般性，以圖1 所示裂隙巖體邊坡滲流問題為例，對于有自由面的無壓滲流問題，區(qū)域Ω（ Ω=Ωw∪ Ωd）上的滲流實際上僅在自由面Γf以下的濕區(qū)Ωw中運動。在多孔連續(xù)介質(zhì)滲流理論中，自由面是一個由計算域上壓力水頭為0 且法向流速為0 的點組成的連續(xù)曲面。而對于裂隙巖體而言，當(dāng)忽略巖塊滲流時，裂隙巖體的滲流是沿裂隙網(wǎng)絡(luò)定向進(jìn)行的，自由面在裂隙處表現(xiàn)為潛水面，依次連接裂隙處的潛水面位置可形成裂隙巖體滲流的自由面邊界（見圖1）。

圖1 裂隙巖體邊坡滲流示意圖 Fig.1 Schematic diagram of seepage flow through a fractured rock slope

假定邊坡上游水位面位于A 點，下游水位面位于D 點，BC 為不透水邊界，E 為溢出點，AE 為由裂隙處的潛水位構(gòu)成的自由面。假定巖塊不透水，地下水僅在自由面以下的裂隙網(wǎng)絡(luò)濕區(qū) wΩ 內(nèi)產(chǎn)生流動。

裂隙網(wǎng)絡(luò)濕區(qū)內(nèi)任一點的總水頭φ 可以定義為

式中：z 為垂直坐標(biāo)分量；p 為裂隙水壓力；wγ 為水的重度。

圖2 裂隙ij 的局部坐標(biāo)系 Fig.2 Local coordinate system of the fracture ij

把巖體中的裂隙概化為光滑平行板模型，在二維裂隙網(wǎng)絡(luò)中，單個裂隙上的滲流是一維運動。如圖1 裂隙ij 的局部坐標(biāo)所示，對任意裂隙ij 建立局部坐標(biāo)系，根據(jù)Darcy 定律，在裂隙ij 內(nèi)的流速可表示為

式中：ijk 為裂隙ij 的滲透系數(shù)。根據(jù)立方定理可得，

式中：g 為重力加速度；ijb 為裂隙ij 的隙寬；υ 為水的運動黏滯系數(shù)。

水在裂隙ij 內(nèi)的流動應(yīng)滿足連續(xù)性方程：

可知，ijv 在裂隙ij 上為常量。設(shè)有m 條裂隙匯聚于節(jié)點i，如圖3 所示。由于節(jié)點i 本身不具備儲存流體的能力，依據(jù)質(zhì)量守恒原理，節(jié)點總流量的代數(shù)和為0，即

圖3 通過節(jié)點i 的滲流量 Fig.3 Schematic diagram for flow distribution at fracture intersection.

水在裂隙網(wǎng)絡(luò)濕區(qū)內(nèi)的運動應(yīng)滿足連續(xù)性方程（4）和（5）以及下列邊界條件： （1）水頭邊界條件

（2）流量邊界條件

式中：Γq為隔水邊界（規(guī)定流出為負(fù)）。

（3）自由面邊界條件

式中：p 為裂隙ij 與自由面的交點。

（4）出滲面邊界條件

至此，給出了帶自由面裂隙網(wǎng)絡(luò)滲流問題的基本方程和邊界條件，自由面fΓ 的位置事先未知，因此，確定其位置成了帶自由面滲流分析的主要內(nèi)容。

3 裂隙網(wǎng)絡(luò)滲流全域上邊值問題的PDE(partial differential equation)提法

為了建立全區(qū)域Ω 上的邊值問題，借鑒初流量法的思想，將Darcy 定律延拓到整個區(qū)域

式中：0v 為初流速，初流速0v 的引入是為了抵消干區(qū) dΩ 實際不存在的流速，其具體形式為

其中

顯然，(10)式中的 vij仍滿足連續(xù)性方程。為了回避Ωw的不確定性，將構(gòu)造一個定義在全域上的邊值問題，該邊值問題的基本方程由連續(xù)性方程、質(zhì)量守恒方程和廣義Darcy 定律構(gòu)成，水頭邊界條件和流量邊界條件保持不變。剩下的邊界AGFD 為Signorini 型互補邊界條件，

綜上所述，定義在整個裂隙網(wǎng)絡(luò)上邊值問題的提法是：求函數(shù)φ，使其滿足控制方程（4）和（5）、外部邊界條件（6）、（7）和式(13)以及內(nèi)部自由面邊界條件式（8）。

4 裂隙網(wǎng)絡(luò)滲流全域上邊值問題的VI 提法

上述PDE 提法中含有未知的內(nèi)部自由邊界fΓ和潛在的溢出面Signorini型互補條件，數(shù)值求解中選取其試探函數(shù)非常困難。通過建立與PDE 提法等價的變分不等式(VI)提法，使得自由邊界fΓ 以及潛在溢出邊界上的互補條件轉(zhuǎn)化為自然邊界條件，從而大大減小數(shù)值求解的難度，并改善數(shù)值計算的收斂性和穩(wěn)定性。

定義試探函數(shù)集VIΦ 為

則在數(shù)學(xué)上與上述PDE 提法完全等價的變分不等式(VI)提法可表述為：在ΦVI中求一函數(shù)φ，使得對?ψ ∈ΦVI，都有如下不等式成立

其中， ξ (φ ，ψ - φ)和 ζφ(ψ - φ)的形式為

為了完成兩種提法的等價性證明，首先利用分部積分展開式 ξ (φ ，ψ - φ ) - ζφ(ψ - φ)：

式中：n 為整個區(qū)域所有裂隙節(jié)點的總數(shù)；im 為匯聚于i 節(jié)點的裂隙條數(shù)。

4.1 PDE 提法?VI 提法

如果φ 是PDE 提法的解，將連續(xù)性方程式（4）和（5），外部邊界條件（6）、（7）和（13）代入式（18），并注意到條件式（14），得

因此，φ 也是VI 提法的解。

4.2 VI 提法?PDE 提法

假定φ 是變分(VI)提法的解， ?ψ ∈ΦVI都滿足式（18）。分別取ψ = φ + θ1和ψ = φ - θ1，其中θ1為在所有裂隙節(jié)點上等于0 的任意函數(shù)，可以得到連續(xù)性方程為

因此，式（18）可以簡化為

在式（21）中分別取ψ = φ + θ2，ψ = φ - θ2，其中θ2為在 Γφ、 Γq和 Γs邊界上裂隙節(jié)點等于0 的任意函數(shù)，可以得到內(nèi)部任意節(jié)點上的質(zhì)量守恒方程為

因而式（21）可以簡化為

對式（23）分別取ψ = φ + θ3，ψ = φ - θ3，其中 θ3為在 Γφ和 Γs邊界上裂隙節(jié)點等于0 的任意函數(shù)，可以得到 Γq上裂隙節(jié)點的流量邊界條件為

因此，式（23）可以進(jìn)一步簡化為

考慮到上式等號右邊的第二項為一常數(shù)，要使不等式對于 ?ψ ∈ΦIV都成立，則等號右邊的第一項必須非負(fù)，即

考慮到式（14），很明顯可以得到

在式（25）中當(dāng)然也可以取 zψ= ，如此就有

由式（27）、（28）以及在sΓ 上iizφ ≤ ，可得sΓ 上的互補條件為

為導(dǎo)出內(nèi)部自由面fΓ 邊界條件，將式（18）中的Ω 表示成 wΩ 和 dΩ 上的積分之和為

利用上面已經(jīng)導(dǎo)出的條件，可得對于?ψ ∈ΦIV，

分別取ψ = φ + θ4和ψ = φ - θ4，θ4為在 Γφ和 Γs邊界上裂隙節(jié)點等于0 的任意函數(shù)，如此可得：

至此，已經(jīng)導(dǎo)出了PDE 提法中的全部方程及邊界條件。

5 有限元數(shù)值求解格式

5.1 自適應(yīng)罰Heaviside 函數(shù)

盡管可以對式（15）進(jìn)行離散來求其數(shù)值解，但由于其中的Heaviside 函數(shù)的不連續(xù)性，將導(dǎo)致數(shù)值不穩(wěn)定性，主要原因是Heaviside 函數(shù)為非連續(xù)的階梯式函數(shù)，在對單元積分時會造成數(shù)值的跳躍。為避免這種數(shù)值跳躍，采用自適應(yīng)罰Heaviside 函數(shù)來代替Heaviside 函數(shù)，其具體表達(dá)式為

式中：2λ 為自由面附近由干區(qū)到濕區(qū)的過渡層厚度，在計算中可以取為穿過自由面裂隙平均長度的1/10；? 為自適應(yīng)放大系數(shù)，為獲得較好的收斂性，可適當(dāng)調(diào)整其大小。

5.2 有限元求解格式

將裂隙網(wǎng)絡(luò)中的單條裂隙看成是有限元中的線單元，單元節(jié)點即為裂隙的兩個端點，對裂隙ij，其上任意一點的水頭可以用節(jié)點i、j 的水頭表示為

式中： Ni、 Nj為形函數(shù)， Ni= 1 -l /lij， Nj= l /lij。

經(jīng)過有限元近似后得到的變分不等式提法為：尋找一水頭向量 φr+1∈I，對于 ?ψ ∈，都有以下不等式成立，

其中

式中：n 表示整個區(qū)域所有裂隙節(jié)點的總數(shù)；r 和r+ 1表示迭代次數(shù)。求解離散形式的變分不等式的算法有很多，本文選用Zheng 所建議的算法[14]。

5.3 迭代過程

Signorini 型變分不等式方法求解有自由面的滲流問題時，溢出邊界的判斷和自由面的迭代是同時進(jìn)行的。為表述方便，將節(jié)點總集N 劃分為4 個子集：

式中：inN 為整個裂隙網(wǎng)絡(luò)區(qū)域內(nèi)部節(jié)點集合。

基于上述離散型VI 提法的無壓滲流分析方法的有限元迭代算法如下：

步驟1：構(gòu)造總體滲透傳導(dǎo)矩陣K ，設(shè)對邊界Γs上的 ?i，都有 i ∈。求解方程 Kφ1=0，得到φ1。

步驟2：令 1r=

（3）如果 |φr+1- φr||＜ ε|| φr||，則轉(zhuǎn)入步驟3，否則令 r = r+ 1，轉(zhuǎn)入（1）。

步驟3：令φ =φr，計算結(jié)束。

6 算 例

基于本文所建立的裂隙網(wǎng)絡(luò)無壓滲流分析的變分不等式提法，采用Visual C++語言和面向?qū)ο蟮某绦蛟O(shè)計方法，研制開發(fā)了相應(yīng)的計算程序FracSeepage，并用多個算例對程序進(jìn)行了考察，驗證了程序的有效性和魯棒性。

6.1 均質(zhì)矩形壩的滲流分析

如圖4 所示，有一均質(zhì)矩形壩，壩高為12 m，寬為10 m，底部邊界隔水，上游水位為10 m，下游水位為2 m。采用間距為0.2 m，有著常水力開度的兩組正交裂隙形成的裂隙網(wǎng)絡(luò)來模擬均質(zhì)矩形壩的滲流，則裂隙的等效隙寬可以表示為

式中：k 為壩體材料的滲透系數(shù)，B 為裂隙間距。計算中選取均質(zhì)壩體的滲透系數(shù) k = 4.13× 10-3mm/s，則裂隙的等效隙寬 b= 0.1 mm。

對于該算例，圖4 給出了自由面的經(jīng)驗解及裂隙網(wǎng)絡(luò)滲流數(shù)值計算結(jié)果，其中經(jīng)驗解為

對于滲流溢出點，經(jīng)驗解為z=4.47 m，數(shù)值解為z=4.27 m。表明數(shù)值解與經(jīng)驗解較為吻合，但數(shù)值解的自由面降落略大于經(jīng)驗解。

圖4 均質(zhì)矩形壩滲流自由面位置的比較 Fig.4 Comparison of locations of the free surface from different approaches

根據(jù)Dupuit 假定，均質(zhì)壩體滲流的總流量可以表示為

式中：1h 和2h 分別表示上下游水位，L 表示壩體寬度。

采用 Dupuit 公式得到壩體滲流量為1.98 × 10-5m2/s，基于裂隙網(wǎng)絡(luò)滲流分析的數(shù)值解為 1.95 × 10-5m2/s，兩者誤差僅為1.79%，表明數(shù)值解的計算結(jié)果是合理可靠的。

6.2 復(fù)雜裂隙網(wǎng)絡(luò)滲流分析

為了研究裂隙巖體的滲透特性和水力耦合特性， 國 際 合 作 項 目 Decovalex[15]第 5 期（Decovalex-2011）的Task 提供了一個含有7 797裂隙、29 814 個節(jié)點的復(fù)雜裂隙網(wǎng)絡(luò)模型即DFN（discrete fracture network），如圖5 所示。產(chǎn)生該DFN 模型的裂隙參數(shù)來自英國Sellafield 地區(qū)的現(xiàn)場地質(zhì)測量[16]，發(fā)育有4 組裂隙，裂隙產(chǎn)狀服從Fisher 分布，裂隙跡長服從負(fù)指數(shù)分布，表1 給出了裂隙的幾何參數(shù)及統(tǒng)計分布規(guī)律。

圖5 裂隙網(wǎng)絡(luò)模型 Fig.5 Geometry of fracture system in the DFN model

表1 裂隙幾何參數(shù) Table 1 Fracture parameters used for DFN generation

本文采用該DFN 模型進(jìn)行裂隙網(wǎng)絡(luò)無壓滲流分析，假定模型左側(cè)水頭邊界為10 m，右側(cè)水頭邊界為20 m，底部為不透水邊界。為了分析裂隙的水力開度對裂隙網(wǎng)絡(luò)滲流的影響，計算考慮了兩種情況：（1）所有裂隙的隙寬為常數(shù)65 um；（2）裂隙的隙寬服從對數(shù)正態(tài)分布且裂隙跡長與隙寬相關(guān)，兩者的關(guān)系函數(shù)為[17]

式中：minl 和maxl 為最小和最大的裂隙跡長；D 為分形維數(shù)；g 為截斷誤差函數(shù)；imb 和inb 為裂隙隙寬的上、下限。

圖6給出了DFN模型在兩種不同水力開度分布條件下裂隙網(wǎng)絡(luò)內(nèi)的流量分布及自由面位置，圖7給出了模型左邊界裂隙出口的流量分布。從圖6、7可以看出，對于所有裂隙隙寬為常水力開度的情況，裂隙網(wǎng)絡(luò)內(nèi)部以及裂隙在下游邊界出口的流量分布相對比較均勻，而且自由面形狀變化平緩；對于裂隙隙寬服從對數(shù)正態(tài)分布且與跡長相關(guān)的情況，裂隙網(wǎng)絡(luò)內(nèi)部及出口的流量分布主要由隙寬較大、跡長較長的裂隙所控制，優(yōu)勢流的效應(yīng)比較明顯，充分反映了長大裂隙的導(dǎo)水作用。同時，受裂隙隙寬隨機(jī)分布影響，自由面形狀在局部起伏較大。圖8還給出了裂隙網(wǎng)絡(luò)內(nèi)部水頭等勢線的分布，可以看出，水頭自裂隙網(wǎng)絡(luò)右側(cè)向左側(cè)隨整體水流方向逐步遞減，與連續(xù)介質(zhì)滲流規(guī)律基本一致。

圖6 裂隙網(wǎng)絡(luò)內(nèi)部流量分布 Fig.6 Flow rates distribution in the DFN model

圖7 下游邊界出口裂隙流量分布 Fig.7 Flow rates distribution in each fracture intersecting the left vertical boundary of the model

圖8 裂隙網(wǎng)絡(luò)水頭等勢線分布 Fig.8 The water head contours in unit of m

7 結(jié) 論

（1）引入初流速來抵消在裂隙網(wǎng)絡(luò)干區(qū)實際不存在的流速，將Darcy 定律擴(kuò)展到整個區(qū)域，定義了在整個區(qū)域上的非線性邊值問題，并將潛在溢出面邊界條件歸納為Signorini 型邊界條件，建立等價的變分不等式提法，從而在理論上解決了出滲點的定位問題。

（2）通過使用連續(xù)型的自適應(yīng)Heaviside 函數(shù) 代替階梯型的Heaviside 函數(shù)，避免了在自由面附近由干區(qū)到濕區(qū)的過渡區(qū)域進(jìn)行積分時形成的數(shù)值跳躍，從而顯著提高了在精確定位自由面和溢出點的數(shù)值求解過程中的穩(wěn)定性和程序的魯棒性。

（3）來自均質(zhì)壩體和Decovalex Task 的復(fù)雜裂隙網(wǎng)絡(luò)算例驗證了本文算法的有效性，同時表明，該算法對復(fù)雜裂隙網(wǎng)絡(luò)具有良好的適應(yīng)性，并能反映裂隙網(wǎng)絡(luò)滲流的結(jié)構(gòu)特征。

[1] 張有天， 陳平， 王鐳. 有自由面滲流分析的初流量法[J]. 水利學(xué)報， 1988，8(1)： 18－26. ZHANG You-tian， CHEN Ping， WANG Lei. Initial flow method for seepage analysis with free surface[J]. Journal of Hydraulic Engineering， 1988， 8(1)： 18－26.

[2] 毛昶熙. 滲流計算分析與控制[M]. 北京： 中國水利水電出版社， 2003.

[3] LONG J C S， REMER J S， WILSON C R， et al. Porous media equivalents for networks of discontinuous fractures[J]. Water Resources Research， 1982， 18(3)： 645－658.

[4] 毛昶熙， 陳平， 李祖貽， 等. 裂隙巖體滲流計算方法研究[J]. 巖土工程學(xué)報， 1991， 13(6)： 1－10. MAO Chang-xi， CHEN Ping， LI Zu-yi， et al. Research on fractured rock mass seepage calculation method[J]. Chinese Journal of Geotechnical Engineering， 1991， 13(6)： 1－10.

[5] LONG J C S， GILMOUR P， WITHERSPOON P A. A method for steady fluid flow in random three-dimensional networks of disc-shaped fractures[J]. Water Resources Research， 1985， 21(8)： 1105－1115.

[6] 張奇華， 鄔愛清. 三維任意裂隙網(wǎng)絡(luò)滲流及其解法[J]. 巖石力學(xué)與工程學(xué)報， 2010， 29(4)： 720－730. ZHANG Qi-hua， WU Ai-qing. Three-dimensional arbitrary fracture network seepage model and solution[J]. Chinese Journal of Rock Mechanics and Engineering， 2010， 29(4)： 720－730.

[7] 王嬡， 速寶玉. 裂隙巖體滲流模型綜述[J]. 水科學(xué)進(jìn)展， 1996， 7(3)： 276－282. WANG Yuan， SU Bao-yu. Comment on the models of seepage flow in fractured rock masses[J]. Advances in Water Science， 1996， 7(3)： 276－282.

[8] BARENBLATT G I. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks[J]. Journal of Applied Mathematical Mechanics， 1960， 24 (5)： 1286－1303.

[9] JING L， MA Y， FANG Z. Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis (DDA) method[J]. International Journal of Rock Mechanics and Mining Sciences， 2001， 38(3)： 343－355.

[10] 王恩志. 剖面二維裂隙網(wǎng)絡(luò)滲流計算方法[J]. 水文地質(zhì)工程地質(zhì)， 1993， 20(4)： 27－29. WANG En-zhi. Seepage calculation method in fissure networks on vertical section[J]. Hydrogeology & Engineering Geology， 1993， 20(4)： 27－29.

[11] 柴軍瑞， 仵彥卿. 巖體裂隙網(wǎng)絡(luò)滲流自由面的確定方法[J]. 工程勘察， 2000， 1： 23－24. CHAI Jun-rui， WU Yan-qing. The method for determination of the position of free surface in fractured rock masses[J]. Geotechnical Investigation & Surveying， 2000， 1： 23－24.

[12] ZHENG H. A new formulation of Signorini’s type for seepage problems with free surfaces[J]. International Journal for Numerical Methods in Engineering， 2005， 64(1)： 1－16.

[13] CHEN Y， ZHOU C， ZHENG H. A numerical solution to seepage problems with complex drainage systems[J]. Computers and Geotechnics， 2008， 35(3)： 383－393.

[14] ZHENG T S， LI L， XU Q Y. An iterative method for the discrete problems of a class of elliptical variational inequalities[J]. Applied Mathematics and Mechanics， 1995， 16(4)： 351－358.

[15] NIRE X. Evaluation of heterogeneity and scaling of fractures in the Borrowdale Volcanic Group in the Sellafield area[R]. UK： [s. n.]， 1997.

[16] JING L， TSANG C F， STEPHANSSON O. DECOVALEX－An international cooperative research project on mathematical models of coupled THM processes for safety analysis of radioactive waste repositories[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts， 1995， 32(5)： 389－398.

[17] ALIREZA B， LANRU J. Hydraulic properties of fractured rock masses with correlated fracture length and aperture[J]. International Journal of Rock Mechanics and Mining Sciences， 2007， 44(5)： 704－719.