基于空間相關性不同組構干熱巖裂隙發育

杜廣盛，陳世江*，馬萬里

(1.內蒙古科技大學礦業與煤炭學院，包頭 014010；2.內蒙古自治區礦業工程重點實驗室，包頭 014010)

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

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

隨著巖石力學研究的不斷深入，人們在巖石力學理論分析、試驗、數值模擬等方面有了很大的進展，針對巖石結構的研究，其中基于離散元理論的數值模型能夠較好地表征巖石間結構破壞時的力學關系[18-21]。現通過引入空間相關性函數，生成可控制空間相關長度進而控制組構分布狀態的花崗巖數字圖像；然后應用數字圖像技術對生成的圖像進行區分，一方面消除巖石圖片采集的誤差，另一方面使分析過程中使用的圖像組分分布可以進行有明確的參數標定；最后采用細觀顆粒流數值模擬方法對花崗巖進行抗拉強度試驗模擬，觀察其破壞過程中裂隙發展、結構狀態對巖石強度的影響，進而從破壞過程中的細觀角度分析不同組分狀態與巖石試件強度的關系。為研究不同組分花崗巖破裂行為提供一種新的分析思路。

1 含組構信息花崗巖數值模型的構建

基于巖石的細觀層次對巖石內部結構破壞規律與強度的關系進行分析，需要獲得數值分析中最重要的巖石力學細觀參數，首先采用相應室內試驗的方式獲得部分巖石物理參數，然后對原巖結構獲取并進行模型構建，進而進行力學試驗與數值試驗匹配細觀參數。最終獲得該巖石較為可靠的細觀參數。

1.1 試件細觀物理參數獲取

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

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

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

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

1.2 不同組構分布矢量化模型構建

采用數碼相機設置相同參數獲取花崗巖表面圖像并進行調整。

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

劃分后采用canny邊緣檢測算法對輪廓進行確定，將獲取到的邊緣圖像導入CAD進行矢量化形成可被識別的矢量邊界模型。

數值圖像處理過程及組分區分如圖3所示。

圖3 試件組分位置結構獲取過程Fig.3 Acquisition process of component position structure of test piece

1.3 花崗巖抗拉強度的測定

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

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

圖4 花崗巖巴西劈裂試驗Fig.4 Brazilian splitting test of granite

按照上述試驗方案，進行了5個試件的巴西劈裂試驗，試驗結果見表1，花崗巖試件的抗拉強度采用國際巖石力學協會推薦計算公式[式(1)]進行計算。

表1 花崗巖試件尺寸及實驗結果Table 1 Size and experimental results of granite specimen

(1)

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

1.4 數值模型的建立及參數匹配

顆粒流分析程序(particle flow code，PFC)采用鍵粒子模型(bond particle model，BPM)平行黏結模型。該模型可以在顆粒與顆粒之間設置黏結組件，當黏結組件生效時可以傳遞力的作用關系，當到達設置黏結組件的極限強度時該組件失效，模型轉變為linear Model線性接觸模型。將上述步驟處理后形成花崗巖試件圖像矢量邊界導入PFC模擬軟件中，生成含有巖石特征的空間相似模型。參照以往花崗巖離散元數值模擬的細觀參數標定[23-26]，結合試驗測得部分物理力學參數，對模型進行參數匹配，經過反復校準，獲得與室內試驗應力曲線相吻合的數值試驗如圖5所示。細觀參數標定結果如表2、表3所示。

圖5 數值模型與室內試驗參數匹配曲線Fig.5 Matching curve of numerical model and indoor test parameters

表2 顆粒單元的細觀參數Table 2 Mesoscopic parameters of the particle unit

表3 接觸黏結模型的細觀參數Table 3 Mesoscopic parameters of the contact bond model

2 不同組構花崗巖劈裂數值試驗

2.1 不同組構花崗巖圖像生成過程

干熱巖開采過程中需要對地下深部的花崗巖進行一定程度的致裂，而在成巖過程中，各個礦物組分不是隨機的具有一定的相關性特征，即一種礦物材料周圍一定區域內的礦物分布狀態有一定的聚集性，而不同礦物則保持一定的區分度。由于不同組分的力學性質有所區別，這對巖石破裂時的裂隙發育有著重要影響。

引入空間相關性函數，通過改變相關參數，生成不同組構特征的巖石圖像，具體生成步驟如下。

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

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

(2)

式(2)中:L為空間相關長度；d為兩個元素點之間的有效距離。

(3)

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

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

圖6 巖石組構圖像數值實現過程Fig.6 Numerical realization of rock fabric image

2.2 不同組構花崗巖試件模型表征

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

2.3 試驗結果及分析

將表中細觀參數賦予模型試件，對生成的含不同空間參數的結構進行巴西劈裂試驗，每種空間參數生成了3組隨機圖像，對同一圖像結構在保證顆粒尺寸在一定范圍內取值隨機變量，得到不同的5個試驗結果，如圖8所示。

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

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

3 不同組構試件裂隙發育的分形分析

大量試驗表明流體通過巖石間裂隙時，裂隙的復雜程度和裂隙總長度分布情況，會影響流體通過裂隙的速度和熱交換的效率。因此對受載狀態下花崗巖裂隙的發育狀態研究有著很大意義。

事實上花崗巖內部不同組構的物理力學性質有所差異，其微裂隙發育狀態也有著明顯的不同，因此微裂隙發育狀態的分析是連接花崗巖組構與強度之間關系的橋梁。

3.1 分形維數計算過程及原理

分形維數是一種能夠對微裂隙分布狀態描述的手段[27-28]。選取其中組構關系較為明顯的圖片，從加載過程中破裂點產生的附近，每間隔500步獲取試件數值模擬加載過程中內部破裂圖像并記錄其相應的破裂點數目，利用MATLAB軟件對產生的破裂形貌圖像去噪和灰度處理并對其分形維數計算，圖9為20 000步L=4的計算結果圖。在計算分形維數時，采用計盒維數算法，即對圖像進行2n(n=1，2，…)不同尺度的劃分，并對其中所含裂隙的進行計數，計算公式見式(4)。本次計算圖像大小設置為512×512像素，最大網格劃分為256×256。各階段破裂點數據及相應分維數據發展趨勢如圖10、圖11所示。

圖9 20 000步L=4試件分形維數計算結果Fig.9 Calculation results of fractal dimension of 20 000 step L =4 specimen

(4)

式(4)中：D為計盒維數；N(C)為在網格尺度為C時所覆蓋分形體的格子數；C為網格尺度。

3.2 不同空間相關長度加載階段破裂特征分形分析

由圖10可知，破裂點以及相應分形維數隨加載步長增加總體呈現增長趨勢。加載破裂初始階段，當L較小時，組分分布離散程度較高。在加載過程中空間相關長度L高的結構產生的破裂點較多，與其相對應的分維值也相對較高；相同加載步長下L=4的結構比L=8的結構強度高。

圖10 不同空間相關系數破裂點及相應分形維數發展趨勢Fig.10 Fracture points with different spatial correlation coefficients and development trend of corresponding fractal dimensions

從裂隙的發育狀態可以看出，當L=4時，花崗巖中不同的成分混合在一起，不同成分間產生的微小結構面也較為均勻的散布在整個試件中，該狀態下加載時產生的裂隙數也相對較小，與其相應的分形維數也較小，不利于應力集中現象的產生；當L=8時，花崗巖中不同成分區分較為明顯，不同組分之間有較大、較明顯的軟弱結構，此時試件加載過程中產生的裂隙數相對較大，其分維數也較高，更加有利于應力集中現象的出現。當加載至17 500步時，L為4、8時，破裂點數為1 587、3 009，此時分維值為1.240 8、1.232 1，說明L=8時破壞點雖多，但其規律性較強，都是在加載軸線附近的較大軟弱結構處產生的破壞，應力相對集中。而L=4破壞點雖少，但其分散性較高，在加載到試件即將破壞的后期，有部分裂隙趨向于分散發展，使其應力分散，從而破壞強度上升。

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

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

圖11 不同組分含量破裂點及相應分形維數發展趨勢Fig.11 Fracture point of different component content and development trend of corresponding fractal dimension

4 結論

通過對含空間相關性結構干熱巖主要成分花崗巖的數值試驗，可以得出以下結論。

(1)從巖石組構的分布狀態出發，通過引入空間相關性函數構建出不同組分分布的花崗巖表面礦物分布圖像，能夠很好地表征此類巖石的細觀結構。

(2)數值試驗表明，礦物的分布和含量情況不同與巖石整體強度有很大關系，不同組分結構的干熱巖在開采時也需要選取不同的致裂工藝。

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