張 亞,彭樂文
(河海大學(xué)地球科學(xué)與工程學(xué)院,江蘇 南京 211100)
Monte-Carlo隨機(jī)模擬技術(shù)是通過一定的隨機(jī)數(shù)生成方法,生成服從某一隨機(jī)變量的概率分布形式的隨機(jī)數(shù)序列。由于裂隙的分布具有很大的隨機(jī)性、不確定性,可以把裂隙幾何參數(shù)當(dāng)作一種隨機(jī)變量,并且大量研究表明裂隙幾何參數(shù)的分布服從一定概率分布形式,因此可以將Monte-Carlo隨機(jī)模擬技術(shù)運用到裂隙巖體中來研究裂隙分布特征以及裂隙網(wǎng)絡(luò)模型。離散裂隙網(wǎng)絡(luò)(DFN,discrete fracture network)將巖體切割成很多由相鄰結(jié)構(gòu)面接觸而組合的塊體,塊體的變形及運動均由結(jié)構(gòu)面控制,并且裂隙網(wǎng)絡(luò)為地下水流唯一的滲流通道。對裂隙巖體來說,離散裂隙網(wǎng)絡(luò)模型比等效連續(xù)介質(zhì)模型更能刻畫地下水實際的滲流情況[1]。
隨著計算機(jī)技術(shù)的發(fā)展,隨機(jī)裂隙網(wǎng)絡(luò)的模擬也成為當(dāng)今研究的熱點之一。王晉麗等[2]利用Monte-Carlo隨機(jī)模擬技術(shù)生成二維裂隙網(wǎng)絡(luò)滲流模型,研究了裂隙網(wǎng)絡(luò)中水流在不同邊界條件下的流動特征;李新強(qiáng)等[3]根據(jù)結(jié)構(gòu)面的分布特征利用計算機(jī)模擬技術(shù)建立了裂隙網(wǎng)絡(luò)模型,并對該模型進(jìn)行了等效巖體滲透張量的求解;張彥洪等[4]通過人工模擬方法建立隨機(jī)裂隙網(wǎng)絡(luò)模型,研究了裂隙開度和滲透壓力對裂隙巖體滲流應(yīng)力耦合特性的影響。在有限的地質(zhì)測量基礎(chǔ)上,馮學(xué)敏等[5]運用Monte-Carlo方法建立了與實際巖體裂隙同分布的隨機(jī)裂隙網(wǎng)絡(luò);敖雪菲等[6]在研究壩基裂隙巖體灌漿時,基于Monte-Carlo方法建立了裂隙巖體三維網(wǎng)絡(luò)模型,結(jié)果表明該法能有效反應(yīng)裂隙真實特征并獲得貼近實際的灌漿結(jié)果;劉日成等[7]通過隨機(jī)模擬技術(shù)建立了2種DFN模型,考慮2種邊界條件研究了DFN的非線性滲流特征;劉泉聲等[8]基于裂隙網(wǎng)絡(luò)模擬技術(shù)Monte-Carlo法,通過Fish語言編制DFN模型生成程序DFN-GEN,研究了應(yīng)力對等效滲透系數(shù)的影響。DFN建模技術(shù)應(yīng)用領(lǐng)域比較廣泛,不僅在裂隙巖體滲流方面得到應(yīng)用,而且郎曉玲等[9]、劉廣峰等[10]也將DFN模型建模技術(shù)運用到了油田開發(fā)等方面。
根據(jù)大量的實踐經(jīng)驗,可以發(fā)現(xiàn)裂隙分布非常復(fù)雜,并且裂隙幾何參數(shù)難以確定,但許多研究表明在裂隙網(wǎng)絡(luò)模擬中假設(shè)結(jié)構(gòu)面為薄圓盤形狀是比較符合實際情況的[11]。因此,研究利用Monte-Carlo隨機(jī)模擬技術(shù)建立裂隙為薄圓盤形狀的三維裂隙網(wǎng)絡(luò)模型,并對模型進(jìn)行驗證、統(tǒng)計窗裂隙交線分析以及裂隙網(wǎng)絡(luò)連通率計算。
Monte-Carlo隨機(jī)模擬方法的基本原理是利用[0,1]區(qū)間標(biāo)準(zhǔn)均分隨機(jī)數(shù),根據(jù)由已知裂隙幾何參數(shù)的密度分布函數(shù)求得的抽樣公式來獲得服從給定裂隙幾何參數(shù)分布形式的隨機(jī)變量。
目前,在工程應(yīng)用中多運用數(shù)學(xué)方法通過計算機(jī)編程來產(chǎn)生所要的隨機(jī)數(shù)。同余法是產(chǎn)生標(biāo)準(zhǔn)均勻分布隨機(jī)數(shù)最為常用的數(shù)學(xué)方法之一,它包括乘法線性同余法、混合線性同余法、二次同余法和加法同余法(見表1),其中研究使用的是混合線性同余法。
表1 常用的同余法
隨機(jī)數(shù)檢驗是對利用隨機(jī)數(shù)生成器產(chǎn)生的偽隨機(jī)數(shù)序列與真正的[0,1]均勻分布隨機(jī)數(shù)序列相似程度的檢驗。常用的檢驗方法包括:參數(shù)檢驗、均勻性檢驗、獨立性檢驗等,其中均勻性檢驗又包括χ2檢驗和K-S檢驗。統(tǒng)計檢驗不能完全肯定或否定某一假設(shè),最好能夠多做幾種統(tǒng)計檢驗,以便有較大的把握保證使用的隨機(jī)數(shù)列有較好的統(tǒng)計性質(zhì)。研究所用的隨機(jī)數(shù)檢驗方法如下:
(1)
該統(tǒng)計量公式漸近地服從χ2(m-1) ,通過查χ2分布表,即可對頻率差異性做顯著性檢驗。
(2) 獨立性檢驗 獨立性檢驗是檢驗一組隨機(jī)數(shù)相互之間是否存在相關(guān)性。設(shè)u1,u2,…,ui為一組按先后順序排列的隨機(jī)數(shù)序列,前后距離為k的隨機(jī)數(shù)相關(guān)系數(shù)為
(2)
因此,統(tǒng)計量為
(3)
當(dāng)在n充分大時,Vk近似服從N(0,1)分布。
在Monte-Carlo隨機(jī)模擬中,隨機(jī)變量是通過在標(biāo)準(zhǔn)均勻分布隨機(jī)數(shù)的基礎(chǔ)上經(jīng)過一定的概率分布函數(shù)公式變換處理得到的,并且服從不同分布形式,這一過程稱為隨機(jī)變量抽樣。抽樣的方法主要包括:直接抽樣法(又稱反函數(shù)法)、復(fù)合抽樣法、舍選抽樣法和變換抽樣法等,其中直接抽樣法最為常用和有效。
直接抽樣法的基本思路為:設(shè)隨機(jī)變量t服從分布p(t),累積分布函數(shù)P(t),則P(t)的值域為[0,1],所以可以均勻分布隨機(jī)數(shù)ui作為P(ti)的函數(shù)值。根據(jù)累積分布函數(shù)與分布函數(shù)的關(guān)系,在P(ti)和ui之間建立t與u對應(yīng)關(guān)系:
(4)
求式(4)的反函數(shù)得
t=p-1(u),
(5)
將p(t)代入式(4),再將式(4)代入求反函數(shù)之后的式(5),即可得到抽樣表達(dá)式。幾種常見分布形式的抽樣表達(dá)式如表2所列。
表2 常見分布的抽樣表達(dá)式
雖然結(jié)構(gòu)面的發(fā)育具有隨機(jī)性,但是同組結(jié)構(gòu)面在巖體內(nèi)分布具有一定的規(guī)律可循,即結(jié)構(gòu)面具有等距性和韻律性。結(jié)構(gòu)面的分布規(guī)律不僅為研究巖體結(jié)構(gòu)力學(xué)和巖體水力學(xué)模型提供了依據(jù),還有助于研究結(jié)構(gòu)面的力學(xué)特性。
結(jié)構(gòu)面的等距性是指同組結(jié)構(gòu)面在一定尺度上、一定范圍內(nèi)具有相同的間距。無論是大尺度還是小尺度,這種等距現(xiàn)象都存在。結(jié)構(gòu)面的等距性主要受以下4種因素的影響:巖石力學(xué)性質(zhì)、應(yīng)力強(qiáng)度、構(gòu)造作用和巖層厚度。除了等距性,結(jié)構(gòu)面還表現(xiàn)出韻律性,即裂隙發(fā)育密集帶和稀疏帶交替相見出現(xiàn),且?guī)c帶之間具有一定間距。結(jié)構(gòu)面韻律性分布的形成機(jī)制一直是學(xué)者研究的熱點,目前公認(rèn)的2種理論模式為:飽和模式與傳播模式。
裂隙巖體中力學(xué)特征及水力特征受裂隙控制,而裂隙的空間分布又可以通過幾何特征來表征。結(jié)構(gòu)面的幾何特征主要是指結(jié)構(gòu)面的形態(tài)以及結(jié)構(gòu)面幾何參數(shù),其中結(jié)構(gòu)面幾何參數(shù)又包括方位、規(guī)模、張開度和粗糙度等。
(1) 結(jié)構(gòu)面形狀 結(jié)構(gòu)面形狀是指在忽略厚度的情況下,在延伸平面上的幾個形狀。由于地質(zhì)條件的復(fù)雜性,在不同形成條件下、結(jié)構(gòu)面的形狀完全不同,很難用一個統(tǒng)一的形狀來概括,因此許多學(xué)者提出了不同的概念模型,主要包括:正交模型(見圖1)、圓盤模型(見圖2)和多邊形模型等。
圖1 結(jié)構(gòu)面正交模型Fig.1 Structural surface orthogonal model
圖2 結(jié)構(gòu)面圓盤模型Fig.2 Structural surface disc model
(2) 結(jié)構(gòu)面參數(shù) 結(jié)構(gòu)面的幾何參數(shù)一般包括:產(chǎn)狀、間距、跡長、粗糙度、張開度、充填物和連通性等,其中間距是指相鄰結(jié)構(gòu)面的垂直距離,用來反映結(jié)構(gòu)面發(fā)育的密集程度,另外根據(jù)所測的平均間距可將巖體進(jìn)行描述(見表3);跡長是指結(jié)構(gòu)面與露頭面的交線,用來描述和評價結(jié)構(gòu)面的連續(xù)性,國際巖石力學(xué)學(xué)會提出的分級標(biāo)準(zhǔn)見表4;張開度是指裂隙兩壁間的垂直距離,是表征滲透性重要的幾何參數(shù)之一。
表3 結(jié)構(gòu)面間距分級表
表4 結(jié)構(gòu)面連續(xù)性分級表
裂隙發(fā)育具有很大的隨機(jī)性,即裂隙幾何參數(shù)具有隨機(jī)性,屬于隨機(jī)變量,因此可以用相應(yīng)的概率分布來描述。目前,常見的裂隙幾何參數(shù)的概率分布函數(shù)包括:均勻分布、負(fù)指數(shù)分布、正態(tài)分布和對數(shù)正態(tài)分布等[12]。
根據(jù)大量的野外資料表明,結(jié)構(gòu)面的產(chǎn)狀、跡長、間距及張開度等幾何要素服從于某種隨機(jī)分布,其中均勻分布一般適用于產(chǎn)狀的分布形式,負(fù)指數(shù)分布一般適用于跡長、間距和張開度分布形式,如圖3所示;正態(tài)分布一般適用于產(chǎn)狀和跡長分布形式,如圖4所示;對數(shù)正態(tài)分布一般適用于跡長、間距和張開度分布形式。常見的幾何要素概率分布規(guī)律見表5[13]。
圖3 間距的負(fù)指數(shù)分布形式Fig.3 Negative exponential distribution of spacing
野外巖體裂隙分布是復(fù)雜的,大量的工程實踐證明,裂隙的分布特征服從一定的規(guī)律。因此,可以通過野外露頭、平洞等裂隙的統(tǒng)計分析,建立裂隙分布的概率模型,并通過Monte-Carlo隨機(jī)模擬方法,在計算機(jī)上建立反映服從該概率分布規(guī)律的三維裂隙網(wǎng)絡(luò)模型。裂隙網(wǎng)絡(luò)模擬過程是根據(jù)現(xiàn)場裂隙統(tǒng)計資料獲得裂隙幾何參數(shù)的概率分布規(guī)律,應(yīng)用Monte-Carlo隨機(jī)模擬方法,對已知的密度函數(shù)進(jìn)行隨機(jī)抽樣,進(jìn)而得到與實際分布函數(shù)相應(yīng)的隨機(jī)變量,進(jìn)而推算出每條裂隙的空間位置坐標(biāo),并結(jié)合裂隙幾何參數(shù)在計算機(jī)上生成形象直觀的三維裂隙網(wǎng)絡(luò)。
圖4 傾向的正態(tài)分布形式Fig.4 Normal distribution of tendency
實際上,隨機(jī)裂隙網(wǎng)絡(luò)建模的過程就是根據(jù)露頭得到有限的裂隙信息推求服從這些概率分布特征的巖體內(nèi)部以及更大區(qū)域的裂隙分布,建立能夠反映現(xiàn)場實際情況的裂隙網(wǎng)絡(luò)模型[12]?;镜娜S裂隙網(wǎng)絡(luò)模擬流程如圖5所示。
根據(jù)文獻(xiàn)[1]中某水電站壩址區(qū)裂隙的幾何參數(shù)統(tǒng)計資料,并結(jié)合Monte-Carlo模擬技術(shù)建立三維裂隙網(wǎng)絡(luò)模型。裂隙的幾何參數(shù)如表6所列,其優(yōu)勢方向為北東向和北西向。對裂隙的幾何分布進(jìn)行擬合,部分?jǐn)M合結(jié)果如圖6、圖7所示。
表5 結(jié)構(gòu)面幾何要素經(jīng)驗概率分布
圖5 三維裂隙網(wǎng)絡(luò)模擬流程Fig.5 3D crack network simulation flowchart
圖6 間距(NE)頻數(shù)直方圖(負(fù)指數(shù)分布)Fig.6 Frequency histogram (Negative exponential distribution) of crack spacing (NE)
圖7 傾向(NE)頻數(shù)直方圖(正態(tài)分布)Fig.7 Frequency histogram (Normal distribution) of tendency (NE)
設(shè)置裂隙生成區(qū)域為200×200×200 m3,正北向為X軸,正東向為Y軸,垂直向為Z軸。由于該區(qū)域發(fā)育2組優(yōu)勢裂隙,故將NE向、NW向分布命名為DFN01和DFN02。2組裂隙的數(shù)量大小由結(jié)構(gòu)面體密度控制,而求得所需要生成的2組裂隙數(shù)量的方法可參考文獻(xiàn)[14]中所述方法,主要思路如下:
由線密度λd和結(jié)構(gòu)面間距d的定義可知兩者互為倒數(shù),則有
(6)
(7)
表6 某水電站壩址區(qū)裂隙幾何參數(shù)統(tǒng)計
(8)
因此,模擬區(qū)域V內(nèi)的結(jié)構(gòu)面數(shù)量n為
(9)
經(jīng)計算可知,NE向裂隙組DFN01所需生成裂隙數(shù)量為307條,NW向DFN02所需61條,由此生成的三維裂隙網(wǎng)絡(luò)模型如圖8所示,其中DFN01標(biāo)記為藍(lán)色,DFN02標(biāo)記為綠色。
圖8 三維裂隙網(wǎng)絡(luò)模型Fig.8 3D crack network model
(1) 模型的驗證 根據(jù)裂隙網(wǎng)絡(luò)模型求得所生成DFN的幾何參數(shù)信息,結(jié)合現(xiàn)場裂隙統(tǒng)計得到的裂隙幾何特征對通過Monte-Carlo隨機(jī)模擬技術(shù)建立的三維裂隙網(wǎng)絡(luò)進(jìn)行驗證。模型驗證一般可以考慮對比模型數(shù)據(jù)與實測數(shù)據(jù)之間結(jié)構(gòu)面平均傾角傾向、平均跡長、線密度等參數(shù)。在模型中設(shè)置4條測線來計算所建DFN的線密度,如圖9所示。對比模型數(shù)據(jù)與實際數(shù)據(jù)之間線密度的差異性,并計算相對誤差百分比,結(jié)果見表7。
圖9 測線及統(tǒng)計窗布置示意圖Fig.9 Measuring line and statistical window layout map
由表7裂隙網(wǎng)絡(luò)線密度模擬值與實測值對比來看,相對誤差百分比均≤10%,認(rèn)為該模型模擬結(jié)果比較貼近實際情況。
(2) 模型分析 研究對模型的分析主要包括統(tǒng)計窗裂隙交線分析和裂隙連通率計算分析。統(tǒng)計窗設(shè)置如圖9中紅色平面所示,統(tǒng)計窗坐標(biāo)分別為(-200,-200,0)(-200,200,0)(200,200,0)(200,-200,0)。統(tǒng)計窗裂隙交際線如圖10所示。統(tǒng)計窗上藍(lán)色裂隙交線為DFN01與統(tǒng)計窗的交線,綠色為DFN02與統(tǒng)計窗的交線。
由圖10可知,DFN01和統(tǒng)計窗相交而成的交線與DFN02和統(tǒng)計窗交線的夾角約為70°,這與表6中所示的2組裂隙走向的差值大致一致。DFN裂隙連通率主要是通過裂隙間的交線的長度來反映的,即交線長度越長,裂隙連通率越大,裂隙連通率計算如圖11所示。由圖11可知,大部分裂隙交線的走向為NE向,這是因為DFN01數(shù)量較多,約為DFN02數(shù)量的5倍;裂隙連通率較大的交線是2組優(yōu)勢裂隙中半徑較大的裂隙相交而形成。
表7 裂隙網(wǎng)絡(luò)線密度模擬值與實測值對比
圖10 統(tǒng)計窗上裂隙交線分布Fig.10 Distribution of the crack intersection on the statistical window
圖11 裂隙連通率計算Fig.11 Crack connectivity calculation
研究采用Monte-Carlo隨機(jī)模擬方法建立了三維裂隙網(wǎng)絡(luò)模型,主要通過混合線性同余法生成隨機(jī)數(shù),并且對隨機(jī)數(shù)進(jìn)行χ2檢驗和獨立性檢驗,結(jié)合實測裂隙統(tǒng)計與分析得到的裂隙幾何參數(shù)概率分布的抽樣公式對隨機(jī)變量進(jìn)行抽樣,得到與實際分布函數(shù)相應(yīng)的隨機(jī)變量,并且假設(shè)結(jié)構(gòu)面為薄圓盤形狀,進(jìn)而建立DFN模型。
在裂隙網(wǎng)絡(luò)模擬之后,需要進(jìn)行模型的驗證,通過對比實測數(shù)據(jù)與模擬數(shù)據(jù)之間的線密度參數(shù)的差異性來檢驗?zāi)P偷目煽啃?驗證結(jié)果表明模型擬合度較高。建立模型之后,除了對模型統(tǒng)計窗裂隙交線分析和裂隙連通率計算分析,還可以近似地求得RQD、滲透系數(shù)張量等參數(shù)。利用Monte-Carlo隨機(jī)模擬技術(shù)對三維裂隙網(wǎng)絡(luò)模型的研究,不僅為裂隙巖體的建模提供了一種很好的思路,并且為后續(xù)對裂隙巖體滲流的研究工作提供了技術(shù)支持。