黃 帆， 姚 池， 周創(chuàng)兵， 楊建華

(南昌大學(xué) 建筑工程學(xué)院， 江西 南昌 330031)

在水利水電工程、礦山工程和核廢料填埋工程中，面臨著大量的裂隙巖體滲流問題。對裂隙巖體滲流問題認(rèn)識不清造成過許多重大工程事故，如1959年法國Malpasset拱壩由于壩肩揚(yáng)壓力異常發(fā)生潰壩，造成500余人的死亡和300億法郎的經(jīng)濟(jì)損失。裂隙巖體在漫長的地質(zhì)演化過程中，發(fā)展成為多裂隙切割而成的“不連續(xù)體”，滲流主要發(fā)生在相互連通的裂隙網(wǎng)絡(luò)中。因此，研究裂隙網(wǎng)絡(luò)滲流、開發(fā)更符合實(shí)際情況的裂隙網(wǎng)絡(luò)滲流模型，具有重要的理論和現(xiàn)實(shí)意義。

巖體裂隙系統(tǒng)的幾何參數(shù)通常通過對野外巖石露頭進(jìn)行統(tǒng)計(jì)調(diào)查獲得。由于采樣范圍有限，因而獲得的裂隙幾何參數(shù)的分布模型包含了很大的不確定性，為了減少這種不確定性，使之更好地反映真實(shí)情況，通常做法是采用離散裂隙網(wǎng)絡(luò)模型(DFN)[1]來擬合裂隙幾何參數(shù)的統(tǒng)計(jì)分布。

目前運(yùn)用蒙特卡羅法模擬裂隙網(wǎng)絡(luò)方面，Dowd等[2]介紹了一種基于高斯結(jié)構(gòu)的新方法，模擬過程中考慮了不同裂隙組之間的空間相關(guān)性：Xu等[3]給出了三維裂隙網(wǎng)絡(luò)模擬技術(shù)；李愛華等[4]通過Fortran程序編制了隨機(jī)裂隙網(wǎng)絡(luò)模擬程序。目前，離散裂隙網(wǎng)絡(luò)模型滲流分析方面的研究有待進(jìn)一步深入。黃勇等[5]應(yīng)用蒙特卡羅法隨機(jī)生成裂隙網(wǎng)絡(luò)系統(tǒng)，基于交叉裂隙網(wǎng)絡(luò)的流量守恒定理，推導(dǎo)了隨機(jī)裂隙網(wǎng)絡(luò)的滲流模型；王晉麗等[6]根據(jù)裂隙幾何參數(shù)分布的統(tǒng)計(jì)特性，基于蒙特卡羅模擬技術(shù)建立了二維裂隙網(wǎng)絡(luò)滲流模型，在不同邊界條件下分析了二維裂隙網(wǎng)絡(luò)中的流體流動(dòng)。Liu等[7]提出了一種代表巖體裂隙網(wǎng)絡(luò)幾何特征的分形模型，將分形特征與裂隙網(wǎng)絡(luò)的等效滲透率聯(lián)系起來，使用蒙特卡羅法生成包含分形維數(shù)的裂隙網(wǎng)絡(luò)模型，其尺寸遵循一定冪律分布。

然而，在以往的DFN建模工作中裂隙各幾何參數(shù)，如裂隙位置、跡長、方向和水力開度等這些參數(shù)通常被認(rèn)為是不相關(guān)的。上述研究在模擬隨機(jī)裂隙網(wǎng)絡(luò)時(shí)也都沒有考慮裂隙跡長和開度的相關(guān)性，但一些研究人員通過對現(xiàn)場測繪結(jié)果的觀察推斷裂隙的長度和孔徑之間存在一定相關(guān)關(guān)系[8-11]，這種相關(guān)性對裂隙網(wǎng)絡(luò)的生成和滲流分析具有重大影響，甚至導(dǎo)致分析結(jié)果出現(xiàn)較大誤差。鑒于此，本文在傳統(tǒng)離散裂隙網(wǎng)絡(luò)模型基礎(chǔ)上，引入裂隙跡線長度和開度的相關(guān)性函數(shù)，建立一種具有參數(shù)相關(guān)性的裂隙網(wǎng)絡(luò)模型，模擬出更符合實(shí)際的裂隙場，對隨機(jī)裂隙網(wǎng)絡(luò)模型的發(fā)展能夠有所補(bǔ)充，并基于此分析了兩種開度條件下裂隙網(wǎng)絡(luò)的滲流，對裂隙網(wǎng)絡(luò)滲流分析具有一定指導(dǎo)意義。

1 裂隙的幾何參數(shù)和假設(shè)

1.1 裂隙位置

在DFN數(shù)值模擬中裂隙的位置通常被假定為服從泊松分布[12]，為了反映裂隙的位置，采用生成服從泊松分布的裂隙中心點(diǎn)來表達(dá)，其泊松過程通過基于遞歸算法生成的隨機(jī)數(shù)來實(shí)現(xiàn)，遞推公式如下：

Ri+1=27Ri-INT(27Ri)

(1)

式中：Ri為[0，1]之間隨機(jī)分布的均勻數(shù)，INT(x)為取整函數(shù)，即取內(nèi)部函數(shù)的x整數(shù)部分，初始值R0由線性同余法生成。如果裂隙中心點(diǎn)生成空間定義二維直角坐標(biāo)系，且范圍在[xa1，xa2]和[ya1，ya2]之內(nèi)，根據(jù)直角坐標(biāo)系的定義，可以生成裂隙中心點(diǎn)坐標(biāo)(xi，yi)，表達(dá)式為：

xi=xa1+Ri+1(xa1-xa2)

(2)

圖1 泊松過程確定的裂隙中心位置(案例尺寸30 m×30 m)Fig.1 Locations of fracture center determined by a Poisson process (example at 30 m×30 m scale)

yi=yb1+Ri(yb1-yb2)

(3)

圖2 不同分形維數(shù)D對裂隙跡線長度的影響Fig.2 Effects of different fractal dimension D on fracture trace length

該算法產(chǎn)生裂隙的中心坐標(biāo)應(yīng)用于所有正方形DFN模型。圖1為尺寸30 m×30 m的DFN數(shù)值模型中裂隙中心點(diǎn)的位置坐標(biāo)。圖中裂隙密度為0.45/m2。

1.2 裂隙跡長

在野外露頭調(diào)查過程中，通過分析一些較大規(guī)模裂隙的幾何特征，文獻(xiàn)[13]給出了裂隙跡長按冪律分布的累積概率密度函數(shù)：

(4)

式中：L為裂隙的跡線長度，Lmax和Lmin為最大和最小裂隙跡長，取值范圍為0.5～250 m[13]，F(xiàn)是區(qū)間為[0，1]的隨機(jī)均勻分布；D是擬合的分形維數(shù)[14]，一般取2.2±0.2，圖2顯示了隨著分形維數(shù)D取不同的值，裂隙跡長累積概率曲線的變化，D取值為2，取值1作為比較。

1.3 裂隙走向

根據(jù)文獻(xiàn)[12]，裂隙走向遵循Fisher分布，這在離散裂隙網(wǎng)絡(luò)數(shù)值建模中被普遍采用，其平均方位角的偏差角累積概率密度函數(shù)由下式給出：

(5)

圖3 不同F(xiàn)isher常數(shù)K對裂隙方位角的偏差角θ的影響Fig.3 Influences of different Fisher constant K on deviation angle θ of fracture

式中：K為Fisher常數(shù)；F是區(qū)間為[0，1]的隨機(jī)均勻分布。由圖3可見，偏差角θ隨K變大而減小，當(dāng)K=10時(shí)接近75 %的裂隙走向偏差角θ<0.5 rad，說明區(qū)域內(nèi)某一組裂隙整體走向變異性不大，而當(dāng)K=1和K=5時(shí)，裂隙走向偏差角θ普遍較大，有的甚至達(dá)到3 rad，與原走向相反。

1.4 裂隙開度

假定通過單個(gè)裂隙的流體運(yùn)動(dòng)遵循立方定律，可以從實(shí)驗(yàn)室剪切流試驗(yàn)間接計(jì)算出裂隙開度[15]：

(6)

式中:h為裂隙開度;Q為流量(m3/s);w為流道寬度(m);υ為流體運(yùn)動(dòng)黏滯系數(shù)(m2/s);i為水力梯度;g為重力加速度(m/s2)。為了簡化計(jì)算，將裂隙開度當(dāng)作恒定值，根據(jù)標(biāo)準(zhǔn)液壓試驗(yàn)結(jié)果[16-17]，裂隙開度服從對數(shù)正態(tài)分布，其密度函數(shù)為:

(7)

式中:a和b為對數(shù)正態(tài)分布的一階矩和二階矩，在使用野外露頭調(diào)查中的裂隙幾何特征數(shù)據(jù)進(jìn)行DNF數(shù)值模擬時(shí)，截?cái)嗔严钝E長和開度的分布函數(shù)是不能避免的。本文研究中，假設(shè)裂隙開度截?cái)喾植己瘮?shù)為:

(8)

式中:h為裂隙開度，ha和hb分別為裂隙開度的上下限，分別取為1和200 μm。當(dāng)開度值小于ha和大于hb時(shí)，截?cái)喾植己瘮?shù)的值為0。對數(shù)正態(tài)分布的截?cái)嗬鄯e分布函數(shù)為:

(9)

經(jīng)過代數(shù)運(yùn)算，上式可寫為:

FLNTR(h)=(g(h)-g(ha))/(g(hb)-g(ha))

(10)

(11)

式中:erfinv(x)是逆誤差函數(shù);F是[0，1]之間的隨機(jī)數(shù)。

1.5 裂隙開度與跡長的關(guān)系

出于簡單考慮，裂隙開度與跡線長度通常被假定為不相關(guān)。但國外學(xué)者在對黏土拉伸變形過程中的斷裂模式進(jìn)行研究時(shí)，發(fā)現(xiàn)了裂隙的長度和開度的相關(guān)性特征[8-11]，并提出了裂隙開度與裂隙跡線長度的冪律相關(guān)函數(shù)：

h=alβ

(12)

式中：β為特征指數(shù)，其取值為0.5～2.0。通過研究裂隙開度與長度的相關(guān)性[11]，指出隨著裂隙長度和位移的增加,裂隙的中值寬度趨于增加，裂隙跡長的截?cái)嗬鄯e概率密度函數(shù)可定義為:

(13)

式中：Lmax和Lmin為最大和最小裂隙跡長，根據(jù)文獻(xiàn)[18]，將式(13)代入式(11)中，得到了一個(gè)基于裂隙開度對數(shù)正態(tài)分布函數(shù)和跡長冪律分布函數(shù)的相關(guān)函數(shù)：

(14)

圖4 隨機(jī)裂隙網(wǎng)絡(luò)Fig.4 Random fracture network

圖5 裂隙跡長和走向的玫瑰圖Fig.5 Fracture trace length and rose diagrams of direction

式中：D為裂隙的分形維數(shù)，ha和hb分別為裂隙開度大小的上下限(同式(8)中取值一致)，初始值由式(6)給出。因此，a和b可以直接作為輸入?yún)?shù)，然后通過式(14)計(jì)算出裂隙開度。

2 隨機(jī)裂隙網(wǎng)絡(luò)數(shù)值模擬

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

本文研究的裂隙數(shù)據(jù)源(表1)基于湖北省黃石市金山店鐵礦項(xiàng)目[19]。該礦采用分段崩落法生產(chǎn)，長期開采后在地表形成多個(gè)大面積的塌陷區(qū)，由于塌陷區(qū)相互交叉覆蓋，造成大量裂隙貫穿地表，裂隙分布很具有代表性。為了更直觀地觀察兩種裂隙生成方法的不同，基于表1中裂隙的初始幾何參數(shù)生成隨機(jī)裂隙網(wǎng)絡(luò)，為了消除尺寸效應(yīng)的影響，采用了統(tǒng)一30 m×30 m的邊界尺寸。圖4為文獻(xiàn)[19]中傳統(tǒng)裂隙生成方法和本文所提出的裂隙生成方法生成的裂隙，可見，傳統(tǒng)的隨機(jī)裂隙網(wǎng)絡(luò)分布更均勻，這是因?yàn)楹笳咄ㄟ^基于遞歸算法生成的隨機(jī)數(shù)來實(shí)現(xiàn)裂隙的位置，即前一個(gè)裂隙位置對后一個(gè)裂隙位置的生成有影響，而不像傳統(tǒng)裂隙生成方法直接采用隨機(jī)均勻分布來實(shí)現(xiàn)，相應(yīng)的跡長和走向的玫瑰圖如圖5所示，圖中顯示本文方法生成的裂隙網(wǎng)絡(luò)中同一裂隙組內(nèi)的各裂隙走向變異性更大，變化區(qū)間尺度大概為30°～35°，而傳統(tǒng)方法生成的裂隙網(wǎng)絡(luò)中同一裂隙組內(nèi)的各裂隙走向變異性更小，變化區(qū)間尺度為25°～30°，因此，裂隙生成方法較傳統(tǒng)方法生成了走向變異性更大的裂隙網(wǎng)絡(luò)，這是因?yàn)橛?jì)算裂隙走向時(shí)對結(jié)果進(jìn)行偏差角估計(jì)而造成的，從而能夠更好地模擬實(shí)際地質(zhì)裂隙隨機(jī)網(wǎng)絡(luò)[12]。

3 與真實(shí)裂隙場結(jié)果對比

圖6 二維真實(shí)裂隙場Fig.6 Two-dimensional real fracture

選取的真實(shí)裂隙場是通過對裂隙進(jìn)行現(xiàn)場露頭調(diào)查，統(tǒng)計(jì)分析裂隙幾何參數(shù)的特性，最后用布爾模擬方法隨機(jī)生成的一個(gè)二維裂隙場[20](圖6)。為了方便比較，真實(shí)場的寬度、長度均取為30 m，總裂隙率約為15%。與真實(shí)場進(jìn)行比較時(shí)假設(shè)所有裂隙的張開度均相同。根據(jù)表1所示統(tǒng)計(jì)數(shù)據(jù)，該裂隙真實(shí)場的裂隙走向大致可以分為3組，分別為: 40° ～70°，190°～240°以及 290°～330°，而裂隙跡長介于0～12 m。將本文所用方法的模擬結(jié)果與真實(shí)場統(tǒng)計(jì)特征進(jìn)行比較，對比如圖7，可以看出模擬結(jié)果與真實(shí)場標(biāo)準(zhǔn)曲線之間擬合效果很好，誤差在5% 以內(nèi)。因此，這種隨機(jī)方法能很好地模擬真實(shí)裂隙場，并可應(yīng)用于后續(xù)的裂隙網(wǎng)絡(luò)滲流分析。

圖7 裂隙統(tǒng)計(jì)特征結(jié)果對比Fig.7 Comparison between statistical characteristics of fracture

圖8 二維真實(shí)裂隙場Fig.8 Two-dimensional real fracture

4 裂隙網(wǎng)絡(luò)穩(wěn)定滲流分析

在離散裂隙網(wǎng)絡(luò)中，巖體中的滲流分布取決于隨機(jī)網(wǎng)絡(luò)中裂隙的幾何參數(shù)和連通率，根據(jù)裂隙滲流規(guī)律，可知無交叉或者只和1條裂隙交叉的裂隙對滲流沒有影響，因此不能作為巖體中的滲流通道，所以在進(jìn)行裂隙網(wǎng)絡(luò)滲流分析前要?jiǎng)h減和截?cái)嗄切]有作用的裂隙，提取出隨機(jī)網(wǎng)絡(luò)中的連通裂隙，圖8顯示了經(jīng)過凈化處理后形成的連通裂隙網(wǎng)絡(luò)。

假定巖塊不透水，流體只在連通裂隙中流動(dòng)。在二維直角坐標(biāo)系中，對裂隙網(wǎng)絡(luò)中任一點(diǎn)的總水頭，可以定義[21]為：

φ=z+p/rw

(15)

式中：z為垂直坐標(biāo)；p為裂隙中水的壓力；rw為水的重度，工程中一般取值為1×104N/m3。本次研究不考慮裂隙交叉點(diǎn)及裂隙表面粗糙度的影響，假定巖體中單個(gè)裂隙的滲流是一維運(yùn)動(dòng)。根據(jù)Darcy定律，裂隙ij內(nèi)的流速可表示為

uij=-kij?φ/?L

(16)

式中：kij為裂隙ij的滲透系數(shù)；L為裂隙跡長。則利用局部立方定律，得到kij-hij函數(shù)關(guān)系式

kij=hij2g/(12υ)

(17)

式中：hij為裂隙隙寬；g為重力加速度，一般取值為9.8 m/s2；υ為水的運(yùn)動(dòng)黏滯系數(shù)，一般取值范圍為0.72×10-6～1.516×10-6m2/s，υ取25°下的近似值0.897×10-6m2/s。單裂隙內(nèi)水的流動(dòng)應(yīng)滿足連續(xù)性方程：

?uij/?L=0

(18)

根據(jù)質(zhì)量守恒原理，裂隙交叉點(diǎn)總流量的代數(shù)和為0，即

∑hijuij=0

(19)

根據(jù)表1生成DNF模型，采用作者之前開發(fā)的裂隙網(wǎng)絡(luò)滲流計(jì)算程序[21]，在裂隙網(wǎng)絡(luò)模型基礎(chǔ)上對巖體滲流進(jìn)行研究。由于在作者之前的研究中對該計(jì)算程序已有驗(yàn)證，因此直接采用該方法，而不對其計(jì)算合理性進(jìn)行重復(fù)驗(yàn)證。假定模型右側(cè)邊界水頭為40 m，左側(cè)邊界水頭為10 m，頂部邊界和底部邊界均不透水。為了研究裂隙開度與跡長的相關(guān)性對巖體裂隙滲透特性的影響，比較兩種情況：(1)裂隙的隙寬根據(jù)相關(guān)性計(jì)算式(14)得出；(2)所有裂隙的隙寬為常數(shù) 65 μm。

圖9顯示了兩種不同裂隙開度條件下裂隙網(wǎng)絡(luò)內(nèi)部流量分布，圖10給出了裂隙網(wǎng)絡(luò)模型下游邊界裂隙出口的流量分布，可以看出隙寬均勻分布時(shí)裂隙網(wǎng)絡(luò)模型內(nèi)部及下游出口的流量分布比較均勻，且滲透率較大，而對于裂隙隙寬根據(jù)相關(guān)性計(jì)算式(14)得出的情況，裂隙網(wǎng)絡(luò)模型內(nèi)部及下游出口的流量分布差異性較大，優(yōu)勢流路徑較為突顯，滲流主要由開度較大和跡長較長的裂隙所控制，細(xì)小裂隙的流量較小，導(dǎo)水效果很差，圖11給出了裂隙網(wǎng)絡(luò)內(nèi)部水壓力的分布，裂隙網(wǎng)絡(luò)左下角的水壓力最大，向網(wǎng)絡(luò)模型右上方向逐步遞減，符合滲流水壓力規(guī)律，證明了程序的有效性。

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

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

圖11 裂隙網(wǎng)絡(luò)內(nèi)部水壓力分布Fig.11 Water pressure distribution in DFN model

5 結(jié) 語

基于蒙特卡羅法，引入裂隙跡線長度和開度的相關(guān)性函數(shù)，建立了一種新的二維裂隙網(wǎng)絡(luò)生成算法，并對生成的裂隙網(wǎng)絡(luò)進(jìn)行了滲流分析。通過與不考慮相關(guān)性的傳統(tǒng)隨機(jī)裂隙網(wǎng)絡(luò)方法比較，得出以下結(jié)論：

(1)隨機(jī)裂隙網(wǎng)絡(luò)生成通過基于遞歸算法生成的隨機(jī)數(shù)來實(shí)現(xiàn)裂隙的位置，取代傳統(tǒng)裂隙生成方法中的隨機(jī)均勻分布，同時(shí)對裂隙走向進(jìn)行偏差角估算，使得裂隙網(wǎng)絡(luò)分布更符合實(shí)際。

(2)基于黏土拉伸模式中裂隙開度和跡長的冪律相關(guān)函數(shù)，同時(shí)結(jié)合裂隙跡長的累積概率密度函數(shù)，提出了裂隙開度和跡長的相關(guān)函數(shù)，使得隨機(jī)裂隙網(wǎng)絡(luò)生成更符合物理規(guī)律，對隨機(jī)裂隙網(wǎng)絡(luò)模型的發(fā)展能夠有所補(bǔ)充。

(3)通過將模擬結(jié)果與真實(shí)場統(tǒng)計(jì)特征進(jìn)行比較，發(fā)現(xiàn)模擬結(jié)果與真實(shí)場標(biāo)準(zhǔn)曲線之間擬合效果很好，誤差在5% 以內(nèi)?？梢姳痉椒芎芎玫貙?shí)現(xiàn)對真實(shí)裂隙場的模擬，驗(yàn)證了改進(jìn)算法的有效性。

(4)對復(fù)雜裂隙網(wǎng)絡(luò)滲流過程的分析表明，考慮開度和跡長的相關(guān)性能反映裂隙網(wǎng)絡(luò)滲流的結(jié)構(gòu)特征，從而證明了所提方法的合理性。

參 考 文 獻(xiàn)：

[1] CUNDALL P A. A computer model for simulating progressive, large scale movement in blocky rock systems[C]∥Symp. ISRM, Nancy, France, 1971: 129- 136.

[2] DOWD P A, XU C, MARDIA K V, et al. A comparison of methods for the stochastic simulation of rock fractures[J]. Mathematical Geology, 2007, 39(7): 697- 714.

[3] XU C, DOWD P. A new computer code for discrete fracture network modelling[J]. Computers & Geosciences, 2010, 36(3): 292- 301.

[4] 李愛華, 朱江. 基于二維裂隙網(wǎng)絡(luò)模擬的巖塊搜索與巖層追蹤方法[J]. 水利水運(yùn)工程學(xué)報(bào), 2014(6): 65- 70. (LI Aihua, ZHU Jiang. A rock block auto-searching and rock layer auto-tracing method based on 2D fracture network simulation[J]. Hydro-Science and Engineering, 2014(6): 65- 70. (in Chinese))

[5] 黃勇, 周志芳. 巖體滲流模擬的二維隨機(jī)裂隙網(wǎng)絡(luò)模型[J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版), 2004, 32(1): 91- 94. (HUANG Yong, ZHOU Zhifang. 2- D stochastic fracture network model for simulation of seepage through fissured rocks[J]. Journal of Hohai University(Natural Science), 2004, 32(1): 91- 94. (in Chinese))

[6] 王晉麗, 陳喜, 黃遠(yuǎn)洋, 等. 巖體裂隙網(wǎng)絡(luò)隨機(jī)生成及連通性研究[J]. 水文地質(zhì)工程地質(zhì), 2013, 40(2): 30- 35. (WANG Jinli, CHEN Xi, HUANG Yuanyang, et al. A study of stochastic generation and connectivity of fracture network in rock mass[J]. Hydrogeology and Engineering Geology, 2013, 40(2): 30- 35. (in Chinese))

[7] LIU R, JIANG Y, LI B, et al. A fractal model for characterizing fluid flow in fractured rock masses based on randomly distributed rock fracture networks[J]. Computers and Geotechnics, 2015, 65: 45- 55.

[8] HATTON C G, MAIN I G, MEREDITH P G. Non-universal scaling of fracture length and opening displacement[J]. Nature, 1994, 367: 160- 162.

[9] RENSHAW C E, PARK J C. Effect of mechanical interactions on the scaling of fracture length and aperture[J]. Nature, 1997, 386(6624): 482.

[10] VERMILYE J M, SCHOLZ C H. Relation between vein length and aperture[J]. Journal of Structural Geology, 1995, 17(3): 423- 434.

[11] STONE D. Sub-surface fracture maps predicted from borehole data: An example from the Eye-Dashwa pluton, Atikokan, Canada[C]∥International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts. Pergamon, 1984, 21(4): 183- 194.

[12] PRIEST S D. Discontinuity analysis for rock engineering[M]. Springer Scienceand Business Media, 2012.

[13] ADLER P M, THOVERT J F. Fractures and fracture networks[M]. Springer Scienceand Business Media, 1999.

[14] 周創(chuàng)兵, 熊文林. 不連續(xù)面的分形維數(shù)及其在滲流分析中的應(yīng)用[J]. 水文地質(zhì)工程地質(zhì), 1996(6): 1- 6. (ZHOU Chuangbing, XIONG Wenlin. Fractal dimension of discontinuous surface and its application in seepage analysis[J]. Progress in Hydrogeology, 1996(6): 1- 6. (in Chinese))

[15] MIN K B, JING L, STEPHANSSON O. Determining the equivalent permeability tensor for fractured rock masses using a stochastic REV approach: method and application to the field data from Sellafield, UK[J]. Hydrogeology Journal, 2004, 12(5): 497- 510.

[16] DVERSTORP B, ADNERSSON J. Application of the discrete fracture network concept with field data: Possibilities of model calibration and validation[J]. Water Resources Research, 1989, 25(3): 540- 550.

[17] CACAS M C, LEDOUX E, MARSILY G, et al. Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model[J]. Water Resources Research, 1990, 26(3): 491- 500.

[18] DE DREUZY J R, DAVY P, BOUR O. Hydraulic properties of two-dimensional random fracture networks following a power law length distribution 2. Permeability of networks based on lognormal distribution of apertures[J]. Water Resources Research, 2001, 37(8): 2079- 2095.

[19] 劉海軍. 基于蒙特卡羅法的巖體裂隙網(wǎng)絡(luò)模型及滲透張量的研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 2011. (LIU Haijun. Studies on discrete fracture network modeling based on Monte-carlo and permeability tensor of rock masses[D]. Harbin: Harbin Institute of Technology, 2011. (in Chinese))

[20] 張弛, 吳劍鋒, 陳干, 等. 裂隙網(wǎng)絡(luò)生成的隨機(jī)模擬研究[J]. 水文地質(zhì)工程地質(zhì), 2015, 42(4): 12- 17. (ZHANG Chi, WU Jianfeng, CHEN Gan, et al. A comparative study of stochastic modeling for generating fracture networks[J]. Hydrogeology and Engineering Geology, 2015, 42(4): 12- 17. (in Chinese))

[21] 姚池, 姜清輝, 葉祖洋, 等. 裂隙網(wǎng)絡(luò)無壓滲流分析的初流量法[J]. 巖土力學(xué), 2012, 33(6): 1896- 1903. (YAO Chi, JIANG Qinghui, YE Zuyang, et al. Initial flow method for unconfined seepage problems of fracture networks[J]. Rock and Soil Mechanics, 2012, 33(6): 1896- 1903. (in Chinese))