李致宇,王明玉,趙建輝,3,王慧芳

(1．中國科學(xué)院大學(xué)資源與環(huán)境學(xué)院,北京 100049;2．中國科學(xué)院大學(xué)水系統(tǒng)安全研究中心,北京 100049;3．河南大學(xué) 計算機與信息工程學(xué)院,河南開封 475001)

三維裂隙滲流模擬中裂隙對相交判斷過程的優(yōu)化

李致宇1,2,王明玉1,2,趙建輝1,2,3,王慧芳1,2

(1．中國科學(xué)院大學(xué)資源與環(huán)境學(xué)院,北京 100049;2．中國科學(xué)院大學(xué)水系統(tǒng)安全研究中心,北京 100049;3．河南大學(xué) 計算機與信息工程學(xué)院,河南開封 475001)

離散裂隙網(wǎng)絡(luò)模型是模擬預(yù)測礦井突水事故有力工具。精確構(gòu)建裂隙網(wǎng)絡(luò)的幾何結(jié)構(gòu)需要判斷大量裂隙之間的拓撲關(guān)系,而這一過程會顯著影響數(shù)值模擬效率。引入定向包圍盒技術(shù)對傳統(tǒng)相交分析“兩步法”進行算法優(yōu)化,提出運算成本估算公式和基于假想裂隙面的參數(shù)選取方法。經(jīng)多組相交實驗表明,在大規(guī)模裂隙模擬中采用定向包圍盒的整體運算時間與范圍盒和軸向包圍盒相比可節(jié)約40%和20%。

三維裂隙滲流;裂隙相交分析;離散裂隙網(wǎng)絡(luò);定向包圍盒

煤礦突水事故與導(dǎo)水裂隙緊密相關(guān)[1-2],大量彼此相交的裂隙交織成網(wǎng)絡(luò)為水源提供了通道,威脅著生產(chǎn)安全。離散裂隙網(wǎng)絡(luò)(DFN)模型是對天然裂隙網(wǎng)絡(luò)的概化與“數(shù)字重現(xiàn)”,是精確刻畫裂隙介質(zhì)滲流的重要基礎(chǔ)。Baecher提出將裂隙概化為三維圓盤[3],該方法幾何描述簡單易于數(shù)值實現(xiàn),此后眾多學(xué)者基于此對裂隙網(wǎng)絡(luò)的構(gòu)建方法與連通性進行了研究[4-7]。筆者后續(xù)的研究也是基于圓盤模型。

構(gòu)建裂隙網(wǎng)絡(luò)涉及一系列復(fù)雜的工作流程,其中對裂隙進行相交判斷是一項基礎(chǔ)性工作。通過相交判斷可以尋找出所有兩兩相交的裂隙組成相交裂隙對,將相交裂隙對彼此連接構(gòu)成裂隙網(wǎng)絡(luò),再剔除水力無效通道形成滲流路徑,后續(xù)數(shù)值模擬工作均在滲流路徑上開展[8]。

相交判斷是對裂隙之間空間關(guān)系的定量描述。有學(xué)者總結(jié)了三維圓盤裂隙精確相交判斷方法[9-13]。基本思路是通過一系列繁瑣的解析幾何運算直接計算兩個圓盤的空間位置關(guān)系,稱“直接法”。“直接法”在小規(guī)模數(shù)值模擬中應(yīng)用效果良好。如果共有N個裂隙,相交判斷次數(shù)為C2N,即0．5N(N-1)次。但隨著模擬規(guī)模不斷擴大,相交判斷計算量陡然增加,其消耗運算機時所占比重已不可忽視,甚至因所需運算資源超出現(xiàn)有計算機硬件水平而無法計算。

于是有學(xué)者提出使用初步去除和詳細去除的“兩步法”逐步篩選對滲流有貢獻的裂隙[14]。筆者將“兩步法”中的第1步稱為“粗相交”判斷,將第2步稱為“精相交”判斷。“粗相交”判斷的目標是最大程度地篩除掉那些絕對不可能相交的裂隙,保留少數(shù)可能相交的裂隙。則在第2步中僅對剩余裂隙進行“精相交”判斷。圖1為“兩步法”的流程示意。由于進行一次“粗相交”判斷的運算成本遠小于一次“精相交”判斷,所以先用低代價算法縮小目標集合,再使用高代價算法進行精確分析可顯著節(jié)約整體計算量。

圖1 裂隙相交判斷“兩步法”流程Fig．1 Two steps for fractures intersection analysis

“粗相交”判斷借鑒了計算機圖形學(xué)與計算幾何學(xué)中的包圍盒技術(shù)[15-16]。范圍盒(VB)是一種最簡單的包圍盒,因其計算方便首先用于相交裂隙對判斷及塊體識別中[14]。隨后有學(xué)者將軸向包圍盒(AABB)引入[17],在此基礎(chǔ)上又與R樹空間索引相結(jié)合,進一步提高了運算性能[18-19]。

在前人研究基礎(chǔ)上,筆者提出使用定向包圍盒(OBB)優(yōu)化“粗相交”判斷過程,進一步為“精相交”減負,以此提升“兩步法”總體運算效率。

1 包圍盒在“粗相交”判斷中的作用

1.1 基本原理

“粗相交”判斷的對象并不是裂隙本身,其思想是將每個裂隙放進各自獨立封閉的包圍盒內(nèi),裂隙相交判斷也退化為兩個包圍盒之間的相交判斷(圖2)。包圍盒不相交則其內(nèi)部的裂隙肯定不相交,但包圍盒相交只能表示其內(nèi)部的裂隙可能相交,故還需對這些“粗相交裂隙對”進行“精相交”判斷。

圖2 3種圓盤裂隙包圍盒的對比Fig．2 Comparison of the three types of bounding box

一個裂隙的包圍盒可由其體對角線兩端點構(gòu)成的點對(P,Q)定義。如圖2(a)所示,(P1i,Q1i)為裂隙O1包圍盒體對角線兩端點。相應(yīng)地也可定義裂隙O2包圍盒為(P2i,Q2i)。其中i取1,2,3,分別表示各點的3個坐標分量x,y,z。當滿足以下關(guān)系式時說明兩個包圍盒相交:

如圖2存在O1,O2,O3三個空間相離的裂隙圓盤。圖2(a)在裂隙中心以圓盤半徑r構(gòu)建邊長為2r的正方體;圖2(b)緊緊包裹住內(nèi)部裂隙且各邊與坐標軸平行;圖2(c)緊緊包裹住內(nèi)部裂隙但各邊與某個自定義坐標系平行。

利用上述3種包圍盒進行“兩步法”相交判斷,結(jié)果見表1。

3種包圍盒經(jīng)過“兩步法”均得出了相同結(jié)果,但范圍盒判斷次數(shù)最多,定向包圍盒次數(shù)最少。這是因為定向包圍盒可以更加緊密地包裹住裂隙,自身所占據(jù)空間越小,與其他包圍盒相交機會也越小,由此提高了對相離裂隙的篩除效果?？啥x指標e反映篩除效果,稱包圍盒的有效性。

表1 相交判斷結(jié)果Table 1 Intersection analysis results

指標e的取值區(qū)間為0～1,e值接近1說明“粗相交”得出的“粗相交裂隙對”數(shù)目與“精相交”判斷的最終結(jié)果相近,反映“粗相交”判斷越有效果;反之若e接近0,說明“粗相交”結(jié)果中包含大量實際上并不相交的裂隙對,仍需經(jīng)過“精相交”再次判斷以最終剔除,反映“粗相交”并未達到預(yù)期效果。利用式(2)可計算圖2中各個包圍盒的e值,范圍盒為0,定向包圍盒為0,軸向包圍盒為1。

上述案例中沒有任何裂隙相交,但采用范圍盒仍需進行3次“粗相交”判斷和3次“精相交”判斷,相比傳統(tǒng)“直接法”僅需3次“精相交”判斷反而增加了計算負擔。這一極端案例僅為凸顯不同包圍盒的構(gòu)造特點,在實際數(shù)值模擬中相交裂隙是普遍存在的,其數(shù)量是影響相交判斷過程的重要因素,將在下節(jié)運算成本分析中予以討論。

1.2 運算成本估算

對包圍盒自身的計算也要耗費一定機時,這在一定程度上會抵消運算效率的提升?？捎萌缦鹿椒謩e對“直接法”和“兩步法”運算成本進行估算:

式中,Cdir和C2step分別為“直接法”和“兩步法”的總運算成本;N為裂隙總數(shù);Cjing和Ccu分別為單次“精相交”和單次“粗相交”判斷的運算成本;Csur為構(gòu)建單個包圍盒的運算成本。

經(jīng)算法分析,若計Ccu的運算量(式(1))為單位運算量1,則Cjing≈13,Csur至少為1。此外,設(shè)Njing為“精相交”判斷結(jié)果,即實際相交的裂隙對數(shù)目;則dint=Njing/[0．5N(N-1)]表示實際相交裂隙對與可能相交裂隙對數(shù)目最大值之比,體現(xiàn)裂隙連通性的部分特征。為消除不同裂隙總數(shù)N的影響,取C2step和Cdir的比值M表示“兩步法”相對于“直接法”的改進程度。

如圖3所示,橫軸為假設(shè)已知幾種裂隙網(wǎng)絡(luò)的相交裂隙對比例(dint),分別進行“直接法”和基于不同e值包圍盒的“兩步法”判斷,計算其所需運算量比值M,并取lg M為縱軸。位于Y=0的橫線為“兩步法”與“直接法”性能優(yōu)勢分割線,該線下方的點表示采用“兩步法”可節(jié)省的計算量;上方的點表示適宜采用“直接法”。從各條曲線變化趨勢可知,隨著裂隙網(wǎng)絡(luò)中相交裂隙對比例(dint)的增加,“兩步法”的性能優(yōu)勢逐漸減小,但采用較高e值包圍盒能更好地延緩這種衰減。所以可對包圍盒進行適當改造以提高其e值,達到進一步優(yōu)化“兩步法”的目的。

圖3 相交判斷運算成本估算Fig．3 Evaluation on intersection analysis cost

2 裂隙圓盤定向包圍盒的建立方法

2.1 空間圓盤的正投影性質(zhì)

設(shè)全局坐標系為O-XYZ,X軸正向的方位角為dx。存在一個裂隙圓盤Of,傾向α,傾角β,半徑r,圓心三維坐標(x0,y0,z0)?？捎嬎愠隽严秷A盤法向量n的各分量[20],即

當圓盤與投影平面平行,投影圖形是半徑為r的正圓;與投影面垂直時是一條長度為2r的線段;與投影平面成任意角時為一橢圓(圖4(a))。投影橢圓的大小和空間朝向因圓盤產(chǎn)狀而變化。投影橢圓的長、短半軸通過式(6)計算,即

其中,a為投影橢圓長半軸;b為短半軸;nprj為法向量的投影向量。3個坐標面上的法向量投影nprj十分方便計算,如將法向量n=(l,m,n)向XOY面投影,只需將第3個分量n設(shè)為0,即nprj-XOY=(l,m,0)。投影橢圓的朝向可由短軸所在直線向量(nprj)確定。

2.2 定向包圍盒的計算方法

2.2.1 自定義坐標系

圖4 裂隙投影橢圓與局部坐標系Fig．4 Fracture projection and local coordinate

軸向包圍盒依照全局坐標建立,而定向包圍盒依照自定義坐標建立。地質(zhì)學(xué)中使用傾向傾角描述裂隙的產(chǎn)狀,產(chǎn)狀定義了空間中一組相互平行的平面。如已知裂隙上一點則可確定裂隙延展平面。為了符合地質(zhì)學(xué)的習(xí)慣,筆者基于產(chǎn)狀描述自定義坐標系的建立過程。自定義坐標系是根據(jù)一個假想裂隙面而建立的。如圖5(a)所示,假設(shè)在全局坐標系O-XYZ中有一個假想裂隙fgue,設(shè)其圓心位于坐標原點O,傾向α,傾角 β。在該裂隙面內(nèi)沿傾向方向設(shè)為 OXc軸,沿走向方向設(shè)為OYc軸,垂直于XcOYc面設(shè)為OZc軸,建立自定義三維直角坐標系O-XcYcZc。

自定義坐標系選定后,需將裂隙圓盤在全局坐標系下的參數(shù)向自定義坐標系轉(zhuǎn)換,包括裂隙中心坐標(x0,y0,z0),傾向α,傾角β和半徑r。其中傾向傾角可由圓盤法向量n=(l,m,n)統(tǒng)一表示。坐標系的變化不影響半徑,其值仍為r。轉(zhuǎn)換后的參數(shù)添加下角標c以示區(qū)別。以中心點坐標為例[11]:

(1)從OZ軸正向看并作為旋轉(zhuǎn)軸,將XOY面逆時針旋轉(zhuǎn)φ1角得到X1OY1,使得OX1軸與裂隙傾向方向一致。

(2)從OY1軸正方向看并作為旋轉(zhuǎn)軸,將X1OZ面逆時針旋轉(zhuǎn)φ2角得到XcOZc,使OXc軸與裂隙傾向線重合。

式中,φ1=dx-α,φ2=β。

同理可算出裂隙圓盤在自定義坐標系下的法向量nc=(lc,mc,nc)。得到自定義坐標下的新參數(shù)后,可按2．1節(jié)所述方法進一步計算新的投影橢圓參數(shù):長半軸ac,短半軸bc和法向量的投影向量nprj_c。

2.2.2 包圍盒的構(gòu)建

圖5 利用假想裂隙面建立自定義坐標系與定向包圍盒Fig．5 Custom coordinate system and OBBs

如圖5(b)所示,定向包圍盒的后續(xù)構(gòu)建步驟就是在自定義坐標系O-XcYcZc下構(gòu)建軸向包圍盒,后文計算均默認在自定義坐標系下進行,符號不再添加下腳標c。包圍盒可由各個坐標面上投影橢圓的軸向外包矩形拼接而成。以XOY面投影橢圓為例,在XOY面上建立如圖4(b)所示局部二維直角坐標系O-X′Y′,坐標原點位于投影橢圓中心,X′軸正向沿橢圓的短軸,方向與nprj-XOY一致,Y′軸正向沿橢圓長軸向上。同時將自定義坐標系的X軸和Y軸平移到投影橢圓中心O,平行于OX軸和OY軸作投影橢圓的軸向最小外包矩形。則投影橢圓在局部二維直角坐標系O-X′Y′上的標準方程為

OX′軸逆時針旋轉(zhuǎn)θ角與OX軸重合,則局部二維坐標系與自定義坐標系的轉(zhuǎn)換關(guān)系為

其中,θ角要使用nprj-XOY的分量求反正切函數(shù)。注意由于反正切函數(shù)的值域為[-0．5π,0．5π],需要根據(jù)終邊所在象限調(diào)整角度。對于XOY面有

假設(shè)一條與Y軸平行的直線X=ΔX與投影橢圓相切,令判別式等于0可解出

這就是軸向最小外包矩形在X軸上相對于投影橢圓中心的坐標增量。同理設(shè)與X軸平行的直線Y=ΔY與投影橢圓相切解出:

則XOY面上軸向最小外包矩形的兩邊長度分別為2ΔX和2ΔY,即得到包圍盒的長與寬。同理可利用YOZ面或XOZ面的投影橢圓計算出包圍盒的高2ΔZ。則包圍盒的對角線端點坐標為

式(15)計算出的(P,Q)是針對自定義坐標系的。如圖2(c)所示,在全局坐標系下定向包圍盒各邊不一定與坐標軸平行,所以無法用體對角線端點表示。在全局坐標下定義定向包圍盒需要知道其8個角點坐標,可參照與2．2．1節(jié)完全相反的步驟將參數(shù)從自定義坐標系轉(zhuǎn)回到全局坐標系。

3 定向包圍盒的實驗效果

假想裂隙面是任意的,所以自定義坐標系也是任意的。選擇不同的自定義坐標系會影響定向包圍盒的使用效果。本章通過數(shù)值實驗研究3個問題:①對比范圍盒、軸向包圍盒和不同自定義坐標系下定向包圍盒的有效性指標e值;②總結(jié)自定義坐標系的優(yōu)選方法;③逐步擴大裂隙模擬規(guī)模,考察不同方法的運算成本變化。

實驗采用Monte Carlo法對一含有3組裂隙的裂隙系統(tǒng)在邊長10～50 m的立方體巖塊范圍內(nèi)進行隨機生成,分別使用范圍盒、軸向包圍盒和不同定向包圍盒進行“兩步法”相交判斷。表2為裂隙模擬參數(shù),總裂隙密度為1．5個/m3。

表2 裂隙系統(tǒng)模擬參數(shù)Table 2 Parameters of the fracture system

3.1 包圍盒有效性指標的變化規(guī)律

自定義坐標系使用上一章節(jié)提出的假想裂隙面法定義,傾向從0°～360°以45°遞增,傾角從0°～90°以15°遞增,其二者組合可產(chǎn)生56個假想裂隙面,對應(yīng)56種定向包圍盒。再加上范圍盒和軸向包圍盒共計58個。為避免單次隨機模擬可能存在的特殊性,實驗中對10 m長巖塊進行了10次隨機實現(xiàn),相應(yīng)生成了 10套裂隙數(shù)據(jù)。相交裂隙對數(shù)目最小值為168 275對,最大值為204 006對,相差3萬多對。但橫向?qū)Ρ雀鞔坞S機實現(xiàn),相同定向包圍盒e值波動比較平穩(wěn),如圖6所示,e值最大變異系數(shù)為0．021,出現(xiàn)在采用傾向90°傾角30°的定向包圍盒上。實驗結(jié)果表明,e值在不同隨機實現(xiàn)中是穩(wěn)定的,故可任取一次隨機實現(xiàn)對e值變化規(guī)律進行研究。

表3為第1次隨機實現(xiàn)的計算結(jié)果。因數(shù)據(jù)較多,表中僅列出了首尾部分數(shù)據(jù)和中間的一些比較重要的數(shù)據(jù)點?？煽闯?0 m長巖塊中共有181 229對裂隙相交。不同包圍盒的“粗相交”結(jié)果差異較大,反映在e值上:使用范圍盒的e值最小為0．243。使用傾向90°傾角45°的定向包圍盒的 e值最大為0．472,是范圍盒的1．9倍,是軸向包圍盒的1．3倍。軸向包圍盒與使用傾向0°傾角0°的定向包圍盒e值相同,這是因為這個假想裂隙面所對應(yīng)的自定義坐標系就是全局坐標系,所建立的包圍盒也自然一樣。

圖6 不同定向包圍盒e值變異系數(shù)(10 m巖體10次實現(xiàn))Fig．6 CVs of e value of different OBBs(10 times realizations on 10 m sample)

為進一步尋找規(guī)律,將e值根據(jù)假想裂隙面的產(chǎn)狀繪成柱狀圖(圖7)。橫坐標為傾向(0°～315°),每個傾向值對應(yīng)7個相鄰的柱圖,表示傾角(0°～90°),縱坐標是e值。從圖7可以看出,e值大于0．4的點依次出現(xiàn)在由如下傾向、傾角組合中:(0°,45°), (90°,30°),(90°,45°),(90°,60°),(180°,45°), (270°,30°),(270°,45°),(270°,60°)。且e值在某些點重復(fù)對稱,如(0°,45°)和(180°,45°),(90°, 30°)和(270°,60°),(90°,45°)和(270°,45°),(90°, 60°)和(270°,30°)的e值分別相同。

表3 第1次隨機實現(xiàn)的實驗結(jié)果(10 m巖體樣本) Table 3 Results of the first realization(10 m sample)

通過分析可發(fā)現(xiàn),e值較大點所對應(yīng)的傾向傾角與各組裂隙的優(yōu)勢產(chǎn)狀相同或相近。如(90°,45°)是第1組裂隙的優(yōu)勢產(chǎn)狀,(90°,30°)和(90°,60°)是圍繞其小幅擺動;(180°,45°)是第2組裂隙的優(yōu)勢產(chǎn)狀;(270°,45°)是第3組裂隙的優(yōu)勢產(chǎn)狀。

至于在某些點上e值出現(xiàn)重復(fù)對稱性,是因為這些點的傾向之間相差180°,傾角之間互為余角,在這些自定義坐標系上構(gòu)建的定向包圍盒也相同。

3.2 自定義坐標系的選擇

根據(jù)上述實驗結(jié)果,采用某組裂隙的優(yōu)勢產(chǎn)狀作為假想裂隙面建立自定義坐標系的定向包圍盒e值較高。這是因為每個裂隙都是圍繞其優(yōu)勢產(chǎn)狀小幅擺動,沿著優(yōu)勢產(chǎn)狀建立的定向包圍盒對內(nèi)部裂隙包裹更緊密,在“粗相交”中的篩除效果好。

當巖體中僅含有1組裂隙時,可以選擇其優(yōu)勢產(chǎn)狀建立自定義坐標系。但當含有多組裂隙時,任意選取某組優(yōu)勢產(chǎn)狀可以得到較高e值,但如何選擇出最高的e值,仍是個值得深入研究的問題。

3.3 總體運算時間測試

運算成本包括空間成本和時間成本兩部分,隨著計算機硬件的飛速發(fā)展,人們更加關(guān)注時間成本。測試以總體運算時間衡量不同包圍盒在相交判斷中的C2step,模擬范圍為 10～50 m 長巖體(裂隙數(shù)目1 500～187 500)。其中定向包圍盒利用3．2節(jié)實驗中e值最高的傾向90°傾角45°建立,并對這3種包圍盒建立R樹索引[18]。進行3次重復(fù)實驗并對總體運算時間取平均值,所采用的計算機軟硬件配置為Intel Core i7 860,4GB內(nèi)存、Windows7 64位、C++和Python。測試結(jié)果如圖8所示。

圖8 不同數(shù)量裂隙相交判斷總體運算時間對比Fig．8 Comparison of total calculation times under different numbers of fractures

從圖8可以看出,隨著裂隙數(shù)量增加每種包圍盒的總體運算時間近似成線性增長,定向包圍盒的總體運算時間均小于軸向包圍盒與范圍盒。當裂隙規(guī)模較小時(10 m長巖體,1 500條裂隙),3種包圍盒的總體運算時間差異不大,實際時間分別為范圍盒26 s,軸向包圍盒20 s,定向包圍盒17 s。當裂隙數(shù)量接近20萬時(樣本長50 m,裂隙數(shù)187 500),范圍盒用時約為3．3 h,軸向包圍盒約為2．0 h,定向包圍盒約為1．6 h。綜合各種不同規(guī)模的實驗,采用定向包圍盒的總體運算時間相比范圍盒和軸向包圍盒分別可節(jié)約40%和20%。運算效率的提高幅度比e值提高幅度略小,這是因為式(4)中的Csur在增加。

4 結(jié)論與展望

(1)包圍盒有效性指標e是對相交分析效率高低的理論度量,定向包圍盒的e值可達范圍盒的1．9倍,軸向包圍盒的1．3倍。

(2)數(shù)值實驗中采用定向包圍盒的總體運算時間相比范圍盒和軸向包圍盒可節(jié)約40%和20%。

(3)構(gòu)建包圍盒所引入的額外運算量是造成理論分析與數(shù)值實驗中效率提升存在差異的主要原因。

仍有幾個方面值得深入研究:

(1)實驗表明選擇某組裂隙的優(yōu)勢產(chǎn)狀構(gòu)建定向包圍盒其e值較高,但當裂隙組數(shù)、密度、半徑等參數(shù)均發(fā)生變化時,如何選擇出e值較大或最大的定向包圍盒仍是個問題。

(2)若使用橢圓模型或多邊形模型概化裂隙,由于幾何屬性更加復(fù)雜,包圍盒構(gòu)建成本必然增加,仍可采用式(4)估算運算成本,但實際改善效果值得通過實驗進一步研究。

(3)包圍盒技術(shù)同空間索引和并行計算技術(shù)的結(jié)合是進一步提高裂隙對相交判斷運算效率的有效途徑。

[1] 國家安全生產(chǎn)監(jiān)督管理總局,國家煤礦安全監(jiān)察局．煤礦防治水規(guī)定[M]．北京:煤炭工業(yè)出版社,2009．

[2] 武 強,崔芳鵬,趙蘇啟,等．礦井水害類型劃分及主要特征分析[J]．煤炭學(xué)報,2013,38(4):561-565．

Wu Qiang,Cui Fangpeng,Zhao Suqi,et al．Type classification and main characteristics of mine water disasters[J]．Journal of China Coal Society,2013,38(4):561-565．

[3] Baecher G B,Lanney N A,Einstein H H．Statistical description of rock properties and sampling[A]．The 18th U．S．Symposium on Rock Mechanics[C]．U．S．A．,1977:1-8．

[4] Long J,Gilmour P,Witherspoon P A．A model for steady fluid flow in random three-dimensional networks of disc-shaped fractures[J]．Water Resources Research,1985,21(8):1105-1115．

[5] Andersson J,Dverstorp B．Conditional simulations of fluid flow in three-dimensional networks of discrete fractures[J]．Water Resources Research,1987,23(10):1876-1886．

[6] Cacas M C,Ledoux E,De Marsily G,et al．Modeling fracture flow with a stochastic discrete fracture network:Calibration and validation:1．the flow model[J]．Water Resources Research,1990,26 (3):479-489．

[7] Xu C,Dowd P A,Mardia K V,et al．A connectivity index for discrete fracture networks[J]．Mathematical Geology,2006,38(5):611-634．

[8] 陳劍平,肖樹芳,王 清．隨機不連續(xù)面三維網(wǎng)絡(luò)計算機模擬原理[M]．長春:東北師范大學(xué)出版社,1995．

[9] 宋曉晨,徐衛(wèi)亞．裂隙巖體滲流模擬的三維離散裂隙網(wǎng)絡(luò)數(shù)值模型(I):裂隙網(wǎng)絡(luò)的隨機生成[J]．巖石力學(xué)與工程學(xué)報, 2004,23(12):2015-2020．

Song Xiaochen,Xu Weiya．Numerical model of three-dimensional discrete fracture network for seepage in fractured rocks(I):Generation of fracture network[J]．Chinese Journal of Rock Mechanics and Engineering,2004,23(12):2015-2020．

[10] 張有天．巖石水力學(xué)與工程[M]．北京:中國水利水電出版社, 2005．

[11] 于青春,薛果夫,陳德基．裂隙巖體一般塊體理論[M]．北京:中國水利水電出版社,2007．

[12] 李新強,楊松青,汪小剛．巖體隨機結(jié)構(gòu)面三維網(wǎng)絡(luò)的生成和可視化技術(shù)[J]．巖石力學(xué)與工程學(xué)報,2007,26(12):2564-2569．

Li Xinqiang,Yang Songqing,Wang Xiaogang．Generation and visualization technologies of three-dimensional network of rockmass stochastic structural plane[J]．Chinese Journal of Rock Mechanics and Engineering,2007,26(12):2564-2569．

[13] Alghalandis Y F,Xu C,Dowd P A．A general framework for fracture intersection analysis:Algorithms and practical applications[A]．Australian Geothermal Energy Conference 2011[C]．Australia, 2011:15-20．

[14] 于青春,武 雄,大西有三．非連續(xù)裂隙網(wǎng)絡(luò)管狀滲流模型及其校正[J]．巖石力學(xué)與工程學(xué)報,2006,25(7):1469-1474．

Yu Qingchun,Wu Xiong,Ohnishi Yuzo．Channel model for fluid flow in discrete fracture network and its modification[J]．Chinese Journal of Rock Mechanics and Engineering,2006,25(7):1469-1474．

[15] 孫家廣,胡事民．計算機圖形學(xué)基礎(chǔ)教程[M]．北京:清華大學(xué)出版社,2005．

[16] Berg M D,Cheong O,Kreveld M V．Computational geometry:Algo

rithms and applications third edition[M]．Berlin:Springer,2008．[17] 劉華梅,王明玉．三維裂隙網(wǎng)絡(luò)滲流路徑識別算法及其優(yōu)化[J]．中國科學(xué)院研究生院學(xué)報,2010,27(4):463-470．

Liu Huamei,Wang Mingyu．Algorithms and their optimization for identification of the connected flow paths for 3D fracture networks in fractured rocks[J]．Journal of the Graduate School of the Chinese Academy of Sciences,2010,27(4):463-470．

[18] Liu Huamei,Wang Mingyu,Song Xianfeng．A new approach for effectively determining fracture network connections in fractured rocks using R tree indexing[J]．Journal of Coal Science and Engineering(China),2011,17(4):401-407．

[19] Liu Huamei,Wang Mingyu,Song Xianfeng．Stepwise approach for identifying pairwiseracture intersections in 3D fracture networks [J]．Journal of Graduate University of Chinese Academy of Sciences,2013,30(1):24-32．

[20] 賈洪彪,唐輝明,劉佑榮,等．巖石結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬理論與工程應(yīng)用[M]．北京:科學(xué)出版社,2008．

Optimizing fracture intersection analysis procedure in 3D fracture network seepage simulation

LI Zhi-yu1,2,WANG Ming-yu1,2,ZHAO Jian-hui1,2,3,WANG Hui-fang1,2
(1.College of Resources and Environment,University of Chinese Academy of Sciences,Beijing 100049,China;2.Center for Water System Security,University of Chinese Academy of Sciences,Beijing 100049,China;3.College of Computer and Information Engineering,Henan University,Kaifeng 475001,China)

The discrete fracture network model is a powerful tool for simulation and prediction on water inrush accidents of coal mine．A considerable amount of topological analysis is required to construct a fracture network in detail, which often significantly impacts on the simulation efficiency．Authors proposed to further optimize the existing“twostepwise approach”intersection analysis by means of introducingan oriented bounding box,presenting a theoretical estimation on the calculation amount and applying a parameter selection approach based on an imaginary fracture plane．Numerical experiments show that,compared with traditional volume box and axis aligned bounding boxapproaches,the total calculation time of the proposed procedure can be reduced by 40%and 20%under a large scale simulation on fractures intersection analysis．

3D fracture network seepage;intersection analysis;discrete fracture network;oriented bounding box

TD745．2

A

0253-9993(2014)11-2286-07

2013-09-23 責任編輯:韓晉平

國家重點基礎(chǔ)研究發(fā)展計劃(973)資助項目(2010CB428801);國家科技重大專項資助項目(2011ZX05060-005);國家高技術(shù)研究發(fā)展計劃(863)資助項目(2011AA050105)

李致宇(1985—),男,北京人,博士研究生。E-mail:greatlzy@sina．com。通訊作者:王明玉(1961—),男,山東樂陵人,教授,博士生導(dǎo)師。E-mail:mwang@ucas．a(chǎn)c．cn

李致宇,王明玉,趙建輝,等．三維裂隙滲流模擬中裂隙對相交判斷過程的優(yōu)化[J]．煤炭學(xué)報,2014,39(11):2286-2292．

10．13225/j．cnki．jccs．2013．1385

Li Zhiyu,Wang Mingyu,Zhao Jianhui,et al．Optimizing fracture intersection analysis procedure in 3D fracture network seepage simulation [J]．Journal of China Coal Society,2014,39(11):2286-2292．doi:10．13225/j．cnki．jccs．2013．1385