3D Voronoi等效晶質模型在巖石破壞細觀研究中的應用*

韓 森， 張欽禮

(1. 中南大學 資源與安全工程學院，湖南 長沙 410083；2. 湖南有色金屬研究院，湖南 長沙 410100)

0 引言

采礦、隧道、水利水電和邊坡等工程與巖石力學特性密切相關，巖石是由不同性質礦物組成的集合體，受到巖石內部礦物晶粒形狀和組分、晶粒間的結構面及微裂隙等微觀結構特性影響[1]導致巖石宏觀力學性質千差萬別，即使同樣的巖石其力學特性也往往表現為各向異性[2]。因此，從巖石微觀結構入手，建立精細化模型，開展巖石破壞特征的細觀研究，對于巖體工程優化設計，具有重要的工程應用價值。

離散元法是一種非連續介質方法，不需要建立復雜的宏觀本構方程，且其將巖石視為由顆粒和不連續面組成的離散體，這與巖石的微觀結構較相似[3]。因此，離散元法在巖石細觀力學特性等方面的研究取得了較快發展，Itasca公司推出的顆粒流程序(Particle Flow Code,PFC)商業軟件，是由Cundall等[4]在離散元理論基礎上發展而來的，能夠實現具有不同力學參數多相介質的建模，并可以方便的展現細觀力學機理。PFC通過建立顆粒之間接觸的本構關系對巖石的力學特性進行數值模擬，首先Cundall和Potyondy提出了黏結顆粒接觸模型和平行黏結接觸模型,并對Lac du Bonnet花崗巖的破壞特征開展了細觀研究[5-7]，但應用于單軸試驗數值模擬時往往會高估巖石抗拉強度，偏離了巖石的實際情況[8]，在此基礎上，簇黏結模型[9]、等效晶質模型[10-11]、平直節理模型[12]以及蔣明鏡微觀膠結模型[13]等被相繼提出，有效解決了巖石單軸抗拉強度過高問題。本文重點建立反映巖石微觀結構的精細化模型，開展巖石破壞特征的細觀研究，根據鄭達等[14]采用電鏡對巖石斷裂后的斷口掃描研究發現微觀斷裂方式主要為穿晶斷裂和沿晶斷裂，等效晶質模型采用平行黏結和光滑節理接觸模型表征巖石的晶粒和晶粒間結構面，能夠展現穿(沿)晶斷裂行為。周喻等[15]構建了二維顆粒流的等效晶質模型體，通過計算與試驗結果對比驗證了等效晶質模型的適宜性。然而，自然界中的巖石為三維形態，巖石中顆粒的形態、組分及顆粒間的結構面等均與二維形態存在較大差異，本文基于3D Voronoi圖，采用PFC3D構建了能夠反映巖石內部晶粒結構的等效晶質模型，分析了晶粒排列均勻化程度及大小等微觀結構的變化對模型宏觀力學性質的影響程度，模擬了板巖單軸壓縮和巴西試驗，并從細觀角度揭示了巖石的破壞機理。

1 3D Voronoi晶質模型構建

1.1 生成3D Voronoi晶粒網

Voronoi晶胞的幾何形狀與多晶粒的微觀結構相似，因此在計算機仿真中常采用Voronoi圖構建多晶粒，本文采用Rhino6.0生成3D Voronoi晶粒網模擬巖石的多晶粒結構，具體步驟為：在PFC3D中隨機生成“粗”顆粒球體，其球體半徑與巖石晶粒尺寸相當，導出球體中心坐標點，形成三維空間點集合，在Rhino6.0中構建由點集產生的3D Voronoi晶粒網。

1.2 構建巖石晶粒結構

構建不同礦物組成成分的巖石晶粒結構，需要以下3個步驟：1)根據巖石的晶粒大小生成尺寸相當的“粗”顆粒體模型，依照巖石的礦物組分情況分組，并記錄分組信息；2)按照上文表述的方法構建晶粒網，并將“粗”顆粒的分組信息逐一對映至每個晶粒結構單元中；3)建立“細”顆粒體模型，將晶粒網覆蓋至該顆粒體模型中，利用向量數量積方法判斷每個“細”顆粒所在的晶粒結構單元，并將分組信息映射至對應的“細”顆粒體中，最后“細”顆粒模型被賦予平行接觸模型，晶粒網被賦予光滑節理接觸模型，平行黏結模型表征巖石晶粒體、光滑節理接觸模型表征巖石晶粒間的結構面，形成了與巖石內部結構相似的結構體。圖1顯示了含晶粒結構的圓柱體巖石模型的構建過程(不同的顏色表示不同的顆粒組分)。

圖1 巖石等效晶質模型構建過程Fig.1 Construction process of rock equivalent crystal model

2 3D Voronoi晶粒結構的優化研究

本節分析了3D Voronoi晶粒排列均勻化程度及大小等微觀結構的變化對模型宏觀力學性質的影響程度，并對模型中的晶粒結構進行優化。

選用粒徑比為1∶1，1∶1.5，1∶2，1∶2.5及1∶3(半徑均值4.5 mm)呈均勻分布“粗”顆粒集構建長×寬×高為150 mm×150 mm×150 mm的晶粒網，其晶粒個數分別為7 073，7 885，6 322，5 959和5 458個。

按照上述方法，沿晶粒網xy平面不同位置隨機構建3個圓柱體(φ50 mm×100 mm)模型，并對平行黏結接觸和光滑節理接觸賦予表1中的細觀參數。

表1 顆粒及接觸模型細觀參數Table 1 Microscopic parameters of particle and contact model

對圓柱體模型進行單軸壓縮(拉伸)試驗，得到的宏觀力學參數與顆粒粒徑比的關系如圖2所示。在單軸壓縮試驗中，模型的抗壓強度受粒徑比的影響較小，但隨著粒徑比的增大其離散性隨之增大；模型的彈性模量隨粒徑比增大逐步增長，但增長趨勢逐漸減??；泊松比波動性較大，但除了粒徑比1∶1時巖樣的泊松比顯著低于其他水平外，泊松比整體受粒徑比的影響并不明顯。在單軸拉伸試驗中，粒徑比1∶1時模型抗拉強度明顯低于其他水平，除此在外，模型抗拉強度受粒徑比的影響較小，但隨著粒徑比的增大其離散性隨之增大。當1∶2.5時彈性模量趨于穩定，且抗壓(拉)強度的離散程度在控制范圍之內，因此，本文選取“粗”顆粒粒徑比為1∶2.5。

為了研究晶粒體大小對模型宏觀力學性質的影響，在上述所構建的晶質等效模型體中(粒徑比1∶2.5)，沿模型體中心位置分別取φ30 mm×60 mm,φ40 mm×80 mm,φ50 mm×100 mm和φ60 mm×120 mm這4種尺寸的圓柱體模型進行單軸壓縮試驗。不同尺寸模型軸向應力-應變曲線如圖3所示。隨著模型尺寸的增大，應力-應變曲線趨于平穩，晶粒尺寸的影響逐漸減小，結合計算效率，本文選取φ50 mm×100 mm的巖樣模型。

圖2 宏觀力學參數與顆粒粒徑比的關系Fig.2 Relationship between macro-mechanical parameters and particle size ratio

圖3 不同尺寸模型軸向應力—應變曲線Fig.3 Axial stress-strain curves of models with different sizes

3 巖石破壞特征細觀研究

為了驗證3D Voronoi晶粒模型適用性和觀察其破壞時的細觀特征，本文選取某黑鎢礦圍巖——粉砂質板巖制成尺寸為φ50 mm×100 mm和φ50 mm×40 mm的巖樣，通過室內單軸壓縮和巴西試驗得出宏觀物理力學參數值。利用掃描電鏡(SEM)對巖樣受壓破壞后的斷口觀察，如圖4所示,巖石斷裂口呈現出穿(沿)晶斷裂的組合形態，圖4中圈定區域內以穿晶斷裂為主，斷裂面不平整，多為撕裂痕跡，其他區域以沿晶斷裂為主，斷裂面相對平整光滑，穿晶斷裂與沿晶斷裂相互穿插呈無規則分布。

圖4 巖樣斷裂斷口斷裂特征Fig.4 Characteristics of rock sample’s rupture surface

將巖樣切割打磨形成5組巖樣，利用x-射線儀對巖樣礦物組成成分進行分析。

該板巖主要由石英、黑云母、長石、綠泥巖及少量礦物成分組成，巖樣礦物組成成分半定量分析結果，如表2所示。

表2 巖樣主要礦物成分分析結果Table 2 Analysis results of main mineral components in rock samples %

根據表2巖樣礦物成分分析結果，選擇表征石英、黑云母、長石、綠泥石占比為28%，45%，22%和5%，粒徑為2.57～6.43 mm(粒徑比1∶2.5)均勻分布的“粗”顆粒集構建的晶粒網，并覆蓋至尺寸為φ50 mm×100 mm 和φ50 mm×40 mm的細顆粒(最小顆粒半徑1.11 mm,粒徑比1∶1.66)圓柱體中，組成晶質等效模型。

利用上述構建的粉質板巖晶質等效模型進行單軸壓縮和巴西試驗模擬。設置伺服加壓(拉)速率0.3，采用試錯方法對PFC3D細觀參數進行標定，標定后的礦物顆粒及平行黏結模型細觀參數見表3，光滑節理接觸模型細觀參數沿用表1中的參數。最終得到了與室內試驗結果相近的宏觀力學參數，如表4所示。圖5～6為應力-軸應變曲線圖。從圖中可以看出，數值模擬的曲線斜率和峰值強度均與室內試驗的曲線接近。峰后均表現出宏觀脆性破壞特征,如圖7所示。從而驗證了3D Voronoi等效晶質模型的適用性。

表3 礦物顆粒及平行黏結接觸模型細觀參數Table 3 Microscopic parameters of mineral grain and parallel bond contact model

在單軸壓縮試驗數值模擬過程中，裂紋沿巖樣晶粒結構面隨機產生，并緩慢增長，該階段應力-應變呈現線性關系；當應力接近峰值的50%時，沿晶斷裂出現快速生長和擴展，同時穿晶斷裂點在加載板接觸面附近出現和緩慢增長；當應力接近峰值的90%時，沿晶斷裂貫通形成主裂紋，穿晶斷裂生長速度加快，主要分布在周邊結構面發生斷裂和錯動的晶粒內部，特別是力學性質較弱的黑云母和綠泥石晶粒內部穿晶斷裂較明顯，這是由局部應力不平衡而引起的，穿晶斷裂擴展和貫通導致了局部失穩，該階段應力與應變的非線性關系明顯；峰值后，應力陡然下降，該階段沿晶和穿晶斷裂增長速率最快，局部失穩演化為宏觀斷裂，斷裂面呈現沿晶和穿晶斷裂的組合形態。巖樣在單軸壓縮破壞過程中以拉伸斷裂為主,在宏觀上表現為脆性斷裂，圖8為單軸壓縮試驗不同階段的斷裂數統計圖。

表4 粉砂質板巖宏觀參數試驗值和模擬值Table 4 Tested and simulated macroscopic parameters of silty slate

圖5 單軸試驗軸向應力(裂紋增長)—應變曲線Fig.5 Axial stress (crack growth)- strain curves of uniaxial tests

圖6 巴西試驗徑向應力(裂紋增長)—應變曲線Fig.6 Radial stress (crack growth)-strain curves of Brazilian tests

圖7 巖樣宏觀破壞后的對照Fig.7 Contrast diagram of of rock samples after macro-failure

圖8 單軸試驗不同階段斷裂數統計圖Fig.8 Statistical chart of fracture number at different stages of uniaxial tests

巴西試驗數值模擬過程中，其裂紋的生成同樣經歷了穩定發展、峰值前的快速增長及峰值后增長速率最快3個階段，但沿晶斷裂首先從加載板接觸面附近生成，然后沿加載板垂直方向擴展和貫通形成主裂紋，隨著晶粒周邊結構面發生斷裂和錯動，誘發晶粒內部的穿晶斷裂，穿晶斷裂擴展和貫通導致了局部失穩，應力達到峰值，最終形成與加載板大致呈垂直狀的宏觀斷裂，巖樣在巴西劈裂破壞過程中同樣以拉伸斷裂為主,在宏觀上表現為脆性斷裂，圖9為巴西試驗不同階段的斷裂數統計圖。

圖9 巴西試驗不同階段斷裂數統計圖Fig.9 Statistical chart of fracture number at different stages of Brazilian tests

4 結論

1)基于3D Voronoi圖，結合PFC3D中顆粒平行黏結及光滑節理模型，構建了一種能夠反映巖石礦物晶粒結構三維精細化建模方法，為巖石破壞特征的細觀研究提供了一種數值分析方法。

2)分析了晶粒排列均勻化程度及大小等微觀結構的變化對模型宏觀力學性質的影響程度，隨晶粒不規則程度的提高，彈性模量逐步增大，單軸抗壓(拉)強度離散型越來越明顯。晶粒大小不變的情況下，隨著模型尺寸的增大，晶粒結構對模型的宏觀力學性質影響逐漸減小。

3)通過裂紋生長經歷的穩定發展、快速增長及峰值后增長速率最快3個階段可以大致表征室內試驗彈性變形階段、非彈性變形階段和峰值后宏觀破壞階段；巖樣在單軸壓縮和巴西劈裂破壞過程中以拉伸斷裂為主,在宏觀上表現為脆性斷裂；巖樣的宏觀斷裂呈現沿晶和穿晶斷裂的組合特征，彈性變形階段裂紋的萌生和擴展以沿晶斷裂為主，而穿晶斷裂的擴展和貫通往往導致局部失穩，表現為非彈性變形階段。