基于空間相關(guān)性不同組構(gòu)干熱巖裂隙發(fā)育

杜廣盛，陳世江*，馬萬(wàn)里

(1.內(nèi)蒙古科技大學(xué)礦業(yè)與煤炭學(xué)院，包頭 014010；2.內(nèi)蒙古自治區(qū)礦業(yè)工程重點(diǎn)實(shí)驗(yàn)室，包頭 014010)

隨著能源需求的清潔性和可再生性不斷攀升，地?zé)崮?、天然氣等能源被進(jìn)一步開發(fā)出來(lái)。早在1974年，美、英、法、德、日等國(guó)家就聯(lián)合建立了研究中心。干熱巖(hot dry rock)作為地?zé)崮荛_采最具價(jià)值的部分，它的開采方式是通過(guò)地下注入高壓水使其在溫度和地應(yīng)力作用下致裂其上部低應(yīng)力巖石，并通過(guò)該裂隙通道從而與周圍巖石交換能量最終回收[1-3]；以往的研究表明干熱巖普遍存在于地下較深部的以花崗巖為主的高溫巖石，因此對(duì)花崗巖的裂隙發(fā)育規(guī)律的研究有著十分重要的意義。

在裂隙研究的過(guò)程中，通常認(rèn)為巖石材料的破壞是由巖石內(nèi)部微裂紋逐漸展開貫通，最終連接形成一條或多條較粗的宏觀裂紋[4-5]。眾多研究表明巖石的細(xì)觀結(jié)構(gòu)特征決定其宏觀力學(xué)響應(yīng)同時(shí)也影響其裂隙發(fā)育過(guò)程，巖石不同的力學(xué)特性就是其組構(gòu)非均質(zhì)性的宏觀體現(xiàn)[6]。部分學(xué)者對(duì)于礦物組分、顆粒尺度分布特征與巖石力學(xué)性能的關(guān)系進(jìn)行了研究[7-9]，結(jié)果表明礦物組分、顆粒尺度分布對(duì)巖石強(qiáng)度有明顯影響。對(duì)于礦物質(zhì)組分含量及分布的研究，大量學(xué)者分別從均質(zhì)、隨機(jī)非均質(zhì)、原巖非均質(zhì)等不同角度進(jìn)行分析研究[10-14]，其中采用原巖非均質(zhì)性研究能夠更好地揭示巖石破壞機(jī)理和裂隙發(fā)展規(guī)律。但該方法采用的圖像處理過(guò)程中存在一定的誤差，同時(shí)對(duì)圖像中巖石組構(gòu)分布狀態(tài)無(wú)法進(jìn)行準(zhǔn)確的定量。唐欣薇等[15]、Tang等[16]、黃文敏等[17]引入了空間相關(guān)性函數(shù)用來(lái)對(duì)巖石組分分布和各向異性進(jìn)行描述，實(shí)現(xiàn)了通過(guò)改變參數(shù)表示不同礦物聚合分布的巖石結(jié)構(gòu)狀態(tài)。

隨著巖石力學(xué)研究的不斷深入，人們?cè)趲r石力學(xué)理論分析、試驗(yàn)、數(shù)值模擬等方面有了很大的進(jìn)展，針對(duì)巖石結(jié)構(gòu)的研究，其中基于離散元理論的數(shù)值模型能夠較好地表征巖石間結(jié)構(gòu)破壞時(shí)的力學(xué)關(guān)系[18-21]。現(xiàn)通過(guò)引入空間相關(guān)性函數(shù)，生成可控制空間相關(guān)長(zhǎng)度進(jìn)而控制組構(gòu)分布狀態(tài)的花崗巖數(shù)字圖像；然后應(yīng)用數(shù)字圖像技術(shù)對(duì)生成的圖像進(jìn)行區(qū)分，一方面消除巖石圖片采集的誤差，另一方面使分析過(guò)程中使用的圖像組分分布可以進(jìn)行有明確的參數(shù)標(biāo)定；最后采用細(xì)觀顆粒流數(shù)值模擬方法對(duì)花崗巖進(jìn)行抗拉強(qiáng)度試驗(yàn)?zāi)M，觀察其破壞過(guò)程中裂隙發(fā)展、結(jié)構(gòu)狀態(tài)對(duì)巖石強(qiáng)度的影響，進(jìn)而從破壞過(guò)程中的細(xì)觀角度分析不同組分狀態(tài)與巖石試件強(qiáng)度的關(guān)系。為研究不同組分花崗巖破裂行為提供一種新的分析思路。

1 含組構(gòu)信息花崗巖數(shù)值模型的構(gòu)建

基于巖石的細(xì)觀層次對(duì)巖石內(nèi)部結(jié)構(gòu)破壞規(guī)律與強(qiáng)度的關(guān)系進(jìn)行分析，需要獲得數(shù)值分析中最重要的巖石力學(xué)細(xì)觀參數(shù)，首先采用相應(yīng)室內(nèi)試驗(yàn)的方式獲得部分巖石物理參數(shù)，然后對(duì)原巖結(jié)構(gòu)獲取并進(jìn)行模型構(gòu)建，進(jìn)而進(jìn)行力學(xué)試驗(yàn)與數(shù)值試驗(yàn)匹配細(xì)觀參數(shù)。最終獲得該巖石較為可靠的細(xì)觀參數(shù)。

1.1 試件細(xì)觀物理參數(shù)獲取

試驗(yàn)采用蘇州紐邁電子科技有限公司生產(chǎn)的miniMR-60核磁共振分析測(cè)試儀(圖1)對(duì)試件進(jìn)行孔隙度檢測(cè)[22]。5個(gè)試樣測(cè)得孔隙度分別為1.01%、0.974%、0.921%、0.982%、0.881%，平均孔隙度為0.953 6%。

圖1 miniMR-60核磁共振分析測(cè)試儀Fig.1 Nuclear magnetic resonance analysis tester miniMR-60

為獲取花崗巖顆粒尺寸及分布情況，采用電子顯微鏡對(duì)花崗巖表面進(jìn)行觀測(cè)，在每個(gè)試件上選取八九個(gè)觀測(cè)點(diǎn)進(jìn)行圖像獲取，對(duì)其中較為明顯的顆粒結(jié)構(gòu)進(jìn)行測(cè)量(圖2)，觀測(cè)結(jié)果表明石英、長(zhǎng)石、云母顆粒直徑分別為0.14～0.25、0.20～0.32、0.07～0.13 mm。

圖2 巖石顆粒顯微成像Fig.2 Microscopic imaging of rock particles

1.2 不同組構(gòu)分布矢量化模型構(gòu)建

采用數(shù)碼相機(jī)設(shè)置相同參數(shù)獲取花崗巖表面圖像并進(jìn)行調(diào)整。

將拍攝的圖像在MATLAB軟件中進(jìn)行數(shù)值圖像處理，首先進(jìn)行灰度處理，使得圖像變成以0～255表示不同灰度的圖像；然后進(jìn)行以灰度值0～96為云母、96～195為石英、196～255為長(zhǎng)石不同閾值組分劃分；最后采用中值去噪獲得含三種組構(gòu)的圖像。

劃分后采用canny邊緣檢測(cè)算法對(duì)輪廓進(jìn)行確定，將獲取到的邊緣圖像導(dǎo)入CAD進(jìn)行矢量化形成可被識(shí)別的矢量邊界模型。

數(shù)值圖像處理過(guò)程及組分區(qū)分如圖3所示。

圖3 試件組分位置結(jié)構(gòu)獲取過(guò)程Fig.3 Acquisition process of component position structure of test piece

1.3 花崗巖抗拉強(qiáng)度的測(cè)定

在花崗巖薄板上鉆取Φ50 mm×25 mm的圓形試樣。選取有明顯組構(gòu)特征的試件，對(duì)其表面進(jìn)行打磨，使平整度控制在±0.05 mm范圍內(nèi)，試驗(yàn)共加工了5個(gè)試件。在鉆取過(guò)程中因鉆進(jìn)速度及薄板平整度等原因，使得加工試件尺寸略有不同。

巖石巴西劈裂試驗(yàn)采用Wance106微機(jī)控制電液伺服萬(wàn)能測(cè)試機(jī)如圖4所示，其最大垂直加載力為10 T，試驗(yàn)采用壓力控制方式，加載速度設(shè)置為200 N/s，在達(dá)到極限強(qiáng)度時(shí)自動(dòng)卸載5 s，自動(dòng)記錄應(yīng)力及加載時(shí)間。將待測(cè)件圓盤直接放置在巴西劈裂載具之間，通過(guò)兩端彈簧升起對(duì)試件抬升防止試件偏移，保證加載方向。

圖4 花崗巖巴西劈裂試驗(yàn)Fig.4 Brazilian splitting test of granite

按照上述試驗(yàn)方案，進(jìn)行了5個(gè)試件的巴西劈裂試驗(yàn)，試驗(yàn)結(jié)果見(jiàn)表1，花崗巖試件的抗拉強(qiáng)度采用國(guó)際巖石力學(xué)協(xié)會(huì)推薦計(jì)算公式[式(1)]進(jìn)行計(jì)算。

表1 花崗巖試件尺寸及實(shí)驗(yàn)結(jié)果Table 1 Size and experimental results of granite specimen

(1)

式(1)中：σt為巖石的抗拉強(qiáng)度；P為巖石的峰值載荷；D為試樣的直徑；t為試樣的厚度。

1.4 數(shù)值模型的建立及參數(shù)匹配

顆粒流分析程序(particle flow code，PFC)采用鍵粒子模型(bond particle model，BPM)平行黏結(jié)模型。該模型可以在顆粒與顆粒之間設(shè)置黏結(jié)組件，當(dāng)黏結(jié)組件生效時(shí)可以傳遞力的作用關(guān)系，當(dāng)?shù)竭_(dá)設(shè)置黏結(jié)組件的極限強(qiáng)度時(shí)該組件失效，模型轉(zhuǎn)變?yōu)閘inear Model線性接觸模型。將上述步驟處理后形成花崗巖試件圖像矢量邊界導(dǎo)入PFC模擬軟件中，生成含有巖石特征的空間相似模型。參照以往花崗巖離散元數(shù)值模擬的細(xì)觀參數(shù)標(biāo)定[23-26]，結(jié)合試驗(yàn)測(cè)得部分物理力學(xué)參數(shù)，對(duì)模型進(jìn)行參數(shù)匹配，經(jīng)過(guò)反復(fù)校準(zhǔn)，獲得與室內(nèi)試驗(yàn)應(yīng)力曲線相吻合的數(shù)值試驗(yàn)如圖5所示。細(xì)觀參數(shù)標(biāo)定結(jié)果如表2、表3所示。

圖5 數(shù)值模型與室內(nèi)試驗(yàn)參數(shù)匹配曲線Fig.5 Matching curve of numerical model and indoor test parameters

表2 顆粒單元的細(xì)觀參數(shù)Table 2 Mesoscopic parameters of the particle unit

表3 接觸黏結(jié)模型的細(xì)觀參數(shù)Table 3 Mesoscopic parameters of the contact bond model

2 不同組構(gòu)花崗巖劈裂數(shù)值試驗(yàn)

2.1 不同組構(gòu)花崗巖圖像生成過(guò)程

干熱巖開采過(guò)程中需要對(duì)地下深部的花崗巖進(jìn)行一定程度的致裂，而在成巖過(guò)程中，各個(gè)礦物組分不是隨機(jī)的具有一定的相關(guān)性特征，即一種礦物材料周圍一定區(qū)域內(nèi)的礦物分布狀態(tài)有一定的聚集性，而不同礦物則保持一定的區(qū)分度。由于不同組分的力學(xué)性質(zhì)有所區(qū)別，這對(duì)巖石破裂時(shí)的裂隙發(fā)育有著重要影響。

引入空間相關(guān)性函數(shù)，通過(guò)改變相關(guān)參數(shù)，生成不同組構(gòu)特征的巖石圖像，具體生成步驟如下。

步驟1在MATLAB軟件中建立一個(gè)二維隨機(jī)數(shù)組A，該數(shù)組中各元素ai,j的值均勻分布在(0,1)內(nèi)。若巖石僅由兩種含量不同的物質(zhì)組成，則可將數(shù)組A中元素進(jìn)行二值化；當(dāng)0

步驟2為了生成礦物組構(gòu)聚集性特征明顯的巖石圖像，須對(duì)圖6(a)巖石圖像進(jìn)行處理，處理方法就是引入空間相關(guān)函數(shù)對(duì)數(shù)組A進(jìn)行變換，此空間相關(guān)函數(shù)式見(jiàn)式(2)，疊加方法見(jiàn)式(3)。變換后得數(shù)組B，具體計(jì)算方法參見(jiàn)文獻(xiàn)[16]。由數(shù)組B可生成礦物組構(gòu)聚集特征明顯的巖石圖像；式(2)中參數(shù)L控制巖石礦物組構(gòu)聚集程度。圖6(b)是L=6 時(shí)的巖石圖像。

(2)

式(2)中:L為空間相關(guān)長(zhǎng)度；d為兩個(gè)元素點(diǎn)之間的有效距離。

(3)

式(3)中：bi,j為數(shù)組B在(i,j)坐標(biāo)下的元素值；ai+p,j+q為與元素bi,j點(diǎn)距離為(p,q)時(shí)數(shù)組A中的元素值；n為計(jì)算區(qū)域取值尺寸，選取為45。

步驟3將求解出的數(shù)組B乘對(duì)應(yīng)系數(shù)轉(zhuǎn)換為數(shù)值范圍在0～255的矩陣，此時(shí)該矩陣所含數(shù)值范圍與圖像所呈現(xiàn)的灰度值保持一致，借助灰度值組分劃分的方法采用相同的參數(shù)進(jìn)行劃分，既0～96為云母；96～195為石英；196～255為長(zhǎng)石。圖6(c)是經(jīng)灰度劃分后的3種不同組分的分布圖。

圖6 巖石組構(gòu)圖像數(shù)值實(shí)現(xiàn)過(guò)程Fig.6 Numerical realization of rock fabric image

2.2 不同組構(gòu)花崗巖試件模型表征

當(dāng)只考慮巖石結(jié)晶形成時(shí)細(xì)觀組分的空間相關(guān)性時(shí)，在計(jì)算過(guò)程中選取空間相關(guān)長(zhǎng)度分別為L(zhǎng)=2、4、6、8，介質(zhì)含量x=0.45，此時(shí)不同的空間相關(guān)長(zhǎng)度計(jì)算結(jié)果如圖7(a)所示。針對(duì)兩相介質(zhì)的礦物成分含量比例的情況，當(dāng)成分含量不同時(shí)，考慮其含空間相關(guān)性的細(xì)觀結(jié)構(gòu)表征。以空間相關(guān)長(zhǎng)度L=6為例，取其中某一介質(zhì)使其為低強(qiáng)度組分，分別賦予含量x=0.25、0.35、0.45，得到不同礦物成分含量計(jì)算結(jié)果如圖7(b)所示，很明顯在選取不同礦物成分含量時(shí)，計(jì)算結(jié)果不同介質(zhì)比例也有明顯的不同。

2.3 試驗(yàn)結(jié)果及分析

將表中細(xì)觀參數(shù)賦予模型試件，對(duì)生成的含不同空間參數(shù)的結(jié)構(gòu)進(jìn)行巴西劈裂試驗(yàn)，每種空間參數(shù)生成了3組隨機(jī)圖像，對(duì)同一圖像結(jié)構(gòu)在保證顆粒尺寸在一定范圍內(nèi)取值隨機(jī)變量，得到不同的5個(gè)試驗(yàn)結(jié)果，如圖8所示。

從圖8數(shù)據(jù)和圖7(a)中可以看出，空間相關(guān)性系數(shù)L能夠影響礦物組分聚合的排布，L值升高礦物組分更趨于大塊聚集。通過(guò)對(duì)不同L值的巴西劈裂模擬試驗(yàn)，L=4、6、8試件極限載荷平均值分別為13.577、11.482 77、11.200 55 kN?？估瓘?qiáng)度隨L的增大而減小。表4中數(shù)據(jù)表明當(dāng)L=6，x=0.25、0.35、0.45時(shí)礦物組分含量有明顯的變化同時(shí)其所對(duì)應(yīng)的試件極限載荷平均值為15.208、13.881、11.482 kN。當(dāng)x所表示的低強(qiáng)度組分含量增大，試件強(qiáng)度明顯降低。

圖7 不同組構(gòu)花崗巖試件圖像Fig.7 Image of different fabric granite specimens

3 不同組構(gòu)試件裂隙發(fā)育的分形分析

大量試驗(yàn)表明流體通過(guò)巖石間裂隙時(shí)，裂隙的復(fù)雜程度和裂隙總長(zhǎng)度分布情況，會(huì)影響流體通過(guò)裂隙的速度和熱交換的效率。因此對(duì)受載狀態(tài)下花崗巖裂隙的發(fā)育狀態(tài)研究有著很大意義。

事實(shí)上花崗巖內(nèi)部不同組構(gòu)的物理力學(xué)性質(zhì)有所差異，其微裂隙發(fā)育狀態(tài)也有著明顯的不同，因此微裂隙發(fā)育狀態(tài)的分析是連接花崗巖組構(gòu)與強(qiáng)度之間關(guān)系的橋梁。

3.1 分形維數(shù)計(jì)算過(guò)程及原理

分形維數(shù)是一種能夠?qū)ξ⒘严斗植紶顟B(tài)描述的手段[27-28]。選取其中組構(gòu)關(guān)系較為明顯的圖片，從加載過(guò)程中破裂點(diǎn)產(chǎn)生的附近，每間隔500步獲取試件數(shù)值模擬加載過(guò)程中內(nèi)部破裂圖像并記錄其相應(yīng)的破裂點(diǎn)數(shù)目，利用MATLAB軟件對(duì)產(chǎn)生的破裂形貌圖像去噪和灰度處理并對(duì)其分形維數(shù)計(jì)算，圖9為20 000步L=4的計(jì)算結(jié)果圖。在計(jì)算分形維數(shù)時(shí)，采用計(jì)盒維數(shù)算法，即對(duì)圖像進(jìn)行2n(n=1，2，…)不同尺度的劃分，并對(duì)其中所含裂隙的進(jìn)行計(jì)數(shù)，計(jì)算公式見(jiàn)式(4)。本次計(jì)算圖像大小設(shè)置為512×512像素，最大網(wǎng)格劃分為256×256。各階段破裂點(diǎn)數(shù)據(jù)及相應(yīng)分維數(shù)據(jù)發(fā)展趨勢(shì)如圖10、圖11所示。

圖9 20 000步L=4試件分形維數(shù)計(jì)算結(jié)果Fig.9 Calculation results of fractal dimension of 20 000 step L =4 specimen

(4)

式(4)中：D為計(jì)盒維數(shù)；N(C)為在網(wǎng)格尺度為C時(shí)所覆蓋分形體的格子數(shù)；C為網(wǎng)格尺度。

3.2 不同空間相關(guān)長(zhǎng)度加載階段破裂特征分形分析

由圖10可知，破裂點(diǎn)以及相應(yīng)分形維數(shù)隨加載步長(zhǎng)增加總體呈現(xiàn)增長(zhǎng)趨勢(shì)。加載破裂初始階段，當(dāng)L較小時(shí)，組分分布離散程度較高。在加載過(guò)程中空間相關(guān)長(zhǎng)度L高的結(jié)構(gòu)產(chǎn)生的破裂點(diǎn)較多，與其相對(duì)應(yīng)的分維值也相對(duì)較高；相同加載步長(zhǎng)下L=4的結(jié)構(gòu)比L=8的結(jié)構(gòu)強(qiáng)度高。

圖10 不同空間相關(guān)系數(shù)破裂點(diǎn)及相應(yīng)分形維數(shù)發(fā)展趨勢(shì)Fig.10 Fracture points with different spatial correlation coefficients and development trend of corresponding fractal dimensions

從裂隙的發(fā)育狀態(tài)可以看出，當(dāng)L=4時(shí)，花崗巖中不同的成分混合在一起，不同成分間產(chǎn)生的微小結(jié)構(gòu)面也較為均勻的散布在整個(gè)試件中，該狀態(tài)下加載時(shí)產(chǎn)生的裂隙數(shù)也相對(duì)較小，與其相應(yīng)的分形維數(shù)也較小，不利于應(yīng)力集中現(xiàn)象的產(chǎn)生；當(dāng)L=8時(shí)，花崗巖中不同成分區(qū)分較為明顯，不同組分之間有較大、較明顯的軟弱結(jié)構(gòu)，此時(shí)試件加載過(guò)程中產(chǎn)生的裂隙數(shù)相對(duì)較大，其分維數(shù)也較高，更加有利于應(yīng)力集中現(xiàn)象的出現(xiàn)。當(dāng)加載至17 500步時(shí)，L為4、8時(shí)，破裂點(diǎn)數(shù)為1 587、3 009，此時(shí)分維值為1.240 8、1.232 1，說(shuō)明L=8時(shí)破壞點(diǎn)雖多，但其規(guī)律性較強(qiáng)，都是在加載軸線附近的較大軟弱結(jié)構(gòu)處產(chǎn)生的破壞，應(yīng)力相對(duì)集中。而L=4破壞點(diǎn)雖少，但其分散性較高，在加載到試件即將破壞的后期，有部分裂隙趨向于分散發(fā)展，使其應(yīng)力分散，從而破壞強(qiáng)度上升。

3.3 不同組分含量加載階段破裂特征分形分析

如圖11所示，不同組分含量試件，在加載過(guò)程中的破壞點(diǎn)數(shù)與破壞圖像分維數(shù)隨加載步長(zhǎng)整體呈現(xiàn)增長(zhǎng)趨勢(shì)，當(dāng)x=0.25時(shí)，云母含量較少，試件中大量成分被設(shè)置為長(zhǎng)石和石英，因此形成的低強(qiáng)度結(jié)構(gòu)面較少，破壞過(guò)程中裂隙產(chǎn)生數(shù)量降低，初始裂隙分維值也因此較低，強(qiáng)度較高；當(dāng)x=0.45時(shí)，云母含量較多，試件中不同組分間的結(jié)構(gòu)面形成較多，破壞過(guò)程中裂隙產(chǎn)生較多，裂隙初始分維值較高，強(qiáng)度較低。當(dāng)加載至18 000步時(shí)，不同組分含量x為0.25、0.45時(shí)，破裂點(diǎn)數(shù)為1 193、2 222，此時(shí)分維值為1.224 9、1.200 5，這說(shuō)明x=0.25試件中的破裂點(diǎn)更加均勻地分布在整個(gè)空間內(nèi)，裂隙趨向于分散發(fā)育，而x=0.45的試件中破裂點(diǎn)由于集中在較弱組構(gòu)附近使得裂隙發(fā)育趨于集中，應(yīng)力集中現(xiàn)象明顯，因此其強(qiáng)度更低。從上述試驗(yàn)數(shù)據(jù)可以看出花崗巖抗拉強(qiáng)度受到破壞點(diǎn)數(shù)和破壞點(diǎn)分維數(shù)的影響。

圖11 不同組分含量破裂點(diǎn)及相應(yīng)分形維數(shù)發(fā)展趨勢(shì)Fig.11 Fracture point of different component content and development trend of corresponding fractal dimension

4 結(jié)論

通過(guò)對(duì)含空間相關(guān)性結(jié)構(gòu)干熱巖主要成分花崗巖的數(shù)值試驗(yàn)，可以得出以下結(jié)論。

(1)從巖石組構(gòu)的分布狀態(tài)出發(fā)，通過(guò)引入空間相關(guān)性函數(shù)構(gòu)建出不同組分分布的花崗巖表面礦物分布圖像，能夠很好地表征此類巖石的細(xì)觀結(jié)構(gòu)。

(2)數(shù)值試驗(yàn)表明，礦物的分布和含量情況不同與巖石整體強(qiáng)度有很大關(guān)系，不同組分結(jié)構(gòu)的干熱巖在開采時(shí)也需要選取不同的致裂工藝。

(3)破壞點(diǎn)分形值分析結(jié)果表明，不同的空間相關(guān)性長(zhǎng)度和組分含量導(dǎo)致巴西盤破裂點(diǎn)的擴(kuò)展方式不同，一種為相對(duì)集中擴(kuò)散，一種為相對(duì)分散擴(kuò)散，破裂呈現(xiàn)分散發(fā)展的試件其峰值載荷相對(duì)較大。較為分散的裂隙能夠有效地進(jìn)行熱量傳遞，但液體的流動(dòng)能力減弱。