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

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

三維裂隙滲流模擬中裂隙對(duì)相交判斷過程的優(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é) 計(jì)算機(jī)與信息工程學(xué)院,河南開封 475001)

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

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

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

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

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

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

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

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

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

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

1.1 基本原理

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

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

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

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

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

3種包圍盒經(jīng)過“兩步法”均得出了相同結(jié)果,但范圍盒判斷次數(shù)最多,定向包圍盒次數(shù)最少。這是因?yàn)槎ㄏ虬鼑锌梢愿泳o密地包裹住裂隙,自身所占據(jù)空間越小,與其他包圍盒相交機(jī)會(huì)也越小,由此提高了對(duì)相離裂隙的篩除效果?？啥x指標(biāo)e反映篩除效果,稱包圍盒的有效性。

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

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

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

1.2 運(yùn)算成本估算

對(duì)包圍盒自身的計(jì)算也要耗費(fèi)一定機(jī)時(shí),這在一定程度上會(huì)抵消運(yùn)算效率的提升。可用如下公式分別對(duì)“直接法”和“兩步法”運(yùn)算成本進(jìn)行估算:

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

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

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

圖3 相交判斷運(yùn)算成本估算Fig．3 Evaluation on intersection analysis cost

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

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

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

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

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

2.2 定向包圍盒的計(jì)算方法

2.2.1 自定義坐標(biāo)系

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 定向包圍盒的實(shí)驗(yàn)效果

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

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

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

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

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

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

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

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

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

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

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

3.2 自定義坐標(biāo)系的選擇

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

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

3.3 總體運(yùn)算時(shí)間測試

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

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

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

4 結(jié)論與展望

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

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

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

仍有幾個(gè)方面值得深入研究:

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

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

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

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

[2] 武 強(qiáng),崔芳鵬,趙蘇啟,等．礦井水害類型劃分及主要特征分析[J]．煤炭學(xué)報(bào),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] 陳劍平,肖樹芳,王 清．隨機(jī)不連續(xù)面三維網(wǎng)絡(luò)計(jì)算機(jī)模擬原理[M]．長春:東北師范大學(xué)出版社,1995．

[9] 宋曉晨,徐衛(wèi)亞．裂隙巖體滲流模擬的三維離散裂隙網(wǎng)絡(luò)數(shù)值模型(I):裂隙網(wǎng)絡(luò)的隨機(jī)生成[J]．巖石力學(xué)與工程學(xué)報(bào), 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] 李新強(qiáng),楊松青,汪小剛．巖體隨機(jī)結(jié)構(gòu)面三維網(wǎng)絡(luò)的生成和可視化技術(shù)[J]．巖石力學(xué)與工程學(xué)報(bào),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é)報(bào),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] 孫家廣,胡事民．計(jì)算機(jī)圖形學(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ò)滲流路徑識(shí)別算法及其優(yōu)化[J]．中國科學(xué)院研究生院學(xué)報(bào),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 責(zé)任編輯:韓晉平

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

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

李致宇,王明玉,趙建輝,等．三維裂隙滲流模擬中裂隙對(duì)相交判斷過程的優(yōu)化[J]．煤炭學(xué)報(bào),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