3D Voronoi等效晶質(zhì)模型在巖石破壞細(xì)觀研究中的應(yīng)用*

韓 森， 張欽禮

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

0 引言

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

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

1 3D Voronoi晶質(zhì)模型構(gòu)建

1.1 生成3D Voronoi晶粒網(wǎng)

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

1.2 構(gòu)建巖石晶粒結(jié)構(gòu)

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

圖1 巖石等效晶質(zhì)模型構(gòu)建過程Fig.1 Construction process of rock equivalent crystal model

2 3D Voronoi晶粒結(jié)構(gòu)的優(yōu)化研究

本節(jié)分析了3D Voronoi晶粒排列均勻化程度及大小等微觀結(jié)構(gòu)的變化對(duì)模型宏觀力學(xué)性質(zhì)的影響程度，并對(duì)模型中的晶粒結(jié)構(gòu)進(jìn)行優(yōu)化。

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

按照上述方法，沿晶粒網(wǎng)xy平面不同位置隨機(jī)構(gòu)建3個(gè)圓柱體(φ50 mm×100 mm)模型，并對(duì)平行黏結(jié)接觸和光滑節(jié)理接觸賦予表1中的細(xì)觀參數(shù)。

表1 顆粒及接觸模型細(xì)觀參數(shù)Table 1 Microscopic parameters of particle and contact model

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

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

圖2 宏觀力學(xué)參數(shù)與顆粒粒徑比的關(guān)系Fig.2 Relationship between macro-mechanical parameters and particle size ratio

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

3 巖石破壞特征細(xì)觀研究

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

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

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

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

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

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

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

表3 礦物顆粒及平行黏結(jié)接觸模型細(xì)觀參數(shù)Table 3 Microscopic parameters of mineral grain and parallel bond contact model

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

表4 粉砂質(zhì)板巖宏觀參數(shù)試驗(yàn)值和模擬值Table 4 Tested and simulated macroscopic parameters of silty slate

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

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

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

圖8 單軸試驗(yàn)不同階段斷裂數(shù)統(tǒng)計(jì)圖Fig.8 Statistical chart of fracture number at different stages of uniaxial tests

巴西試驗(yàn)數(shù)值模擬過程中，其裂紋的生成同樣經(jīng)歷了穩(wěn)定發(fā)展、峰值前的快速增長及峰值后增長速率最快3個(gè)階段，但沿晶斷裂首先從加載板接觸面附近生成，然后沿加載板垂直方向擴(kuò)展和貫通形成主裂紋，隨著晶粒周邊結(jié)構(gòu)面發(fā)生斷裂和錯(cuò)動(dòng)，誘發(fā)晶粒內(nèi)部的穿晶斷裂，穿晶斷裂擴(kuò)展和貫通導(dǎo)致了局部失穩(wěn)，應(yīng)力達(dá)到峰值，最終形成與加載板大致呈垂直狀的宏觀斷裂，巖樣在巴西劈裂破壞過程中同樣以拉伸斷裂為主,在宏觀上表現(xiàn)為脆性斷裂，圖9為巴西試驗(yàn)不同階段的斷裂數(shù)統(tǒng)計(jì)圖。

圖9 巴西試驗(yàn)不同階段斷裂數(shù)統(tǒng)計(jì)圖Fig.9 Statistical chart of fracture number at different stages of Brazilian tests

4 結(jié)論

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

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

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