邱 鑫 林 緬, 鄭思平 陳天宇

?(中國科學(xué)院力學(xué)研究所，北京100190)

?(中國科學(xué)院大學(xué)，北京100049)

??(東北大學(xué)深部金屬礦山安全開采教育部重點(diǎn)實(shí)驗(yàn)室，沈陽110819)

巖石顆粒尺度的微觀結(jié)構(gòu)控制著巖石的微觀力學(xué)行為，從而控制著巖石的宏觀力學(xué)反應(yīng)。在顆粒尺度上，巖石的晶粒形狀、晶粒礦物剛度、礦物接觸分布和其他缺陷等導(dǎo)致了巖石微觀結(jié)構(gòu)的非均質(zhì)性[1-3]，這將引起內(nèi)部應(yīng)力分布的不均一性，影響斷裂萌生和斷裂擴(kuò)展的力學(xué)行為[4-5]，從而控制巖石在載荷作用下的破壞、變形能力和裂紋模式。然而，對(duì)于非均質(zhì)性如何影響脆性巖石破壞機(jī)制的理解仍不完整[6]。

隨著計(jì)算機(jī)計(jì)算能力的不斷提高，數(shù)值模擬技術(shù)不斷發(fā)展并提供了模擬巖石微觀結(jié)構(gòu)和研究其巖石力學(xué)特性的可能性。目前在模擬巖石的顆粒結(jié)構(gòu)方面，主要有以下四種方法：(1) 顆粒流 PFC 實(shí)現(xiàn)的圓盤狀顆粒[7-8]；(2) RFPA2D 代碼實(shí)現(xiàn)的正方形顆粒[9]；(3) 離散元 UDEC[10]以及混合有限離散代碼 Y-Geo[11]和 ELFEN[12]實(shí)現(xiàn)的三角形顆粒；(4) UDEC 實(shí)現(xiàn)的多邊形顆粒[13-14]。從巖石實(shí)際的微觀結(jié)構(gòu)觀察來看，多邊形顆粒結(jié)構(gòu)更加接近于真實(shí)的情況。在UDEC 中，泰森多邊形劃分方法被應(yīng)用于生成巖石的多邊形顆粒結(jié)構(gòu)，即UDEC–Voronoi模型。在 UDEC–Voronoi 模型中，多邊形隨機(jī)形狀的巖石顆粒與顆粒間接觸嵌合自鎖，接觸面的力學(xué)特性由法向剛度、切向剛度、內(nèi)聚力、抗拉強(qiáng)度和摩擦角等參數(shù)決定。

UDEC–Voronoi 模型被開發(fā)后，許多學(xué)者對(duì)其進(jìn)行了研究和應(yīng)用。Fabjan 等[15]對(duì)該模型進(jìn)行了廣泛的微觀參數(shù)敏感性分析，提供了有效巖石模型校準(zhǔn)的指南。同時(shí)，通過與連續(xù)介質(zhì)模型的比較，提出UDEC–Voronoi 模型可以更好地表征巖石的膨脹、裂縫形態(tài)和峰后行為。而孫博等[16]則對(duì)比了Trigon模型和Voronoi 模型在微觀參數(shù)敏感性方面的差異，展現(xiàn)了多邊形顆粒結(jié)構(gòu)模型和三角形顆粒結(jié)構(gòu)模型在表征巖石力學(xué)特征時(shí)存在的差異。目前，UDEC–Voronoi 模型已成功模擬了砂巖、火成巖、油頁巖等巖石的微觀結(jié)構(gòu)，并在巖石的裂紋起裂與擴(kuò)展規(guī)律、壓縮強(qiáng)度和變形破壞特征等方面的研究中得到了應(yīng)用[17-20]。

本研究介紹了UDEC–Voronoi 模型的構(gòu)建方法，構(gòu)建了致密巖石的微結(jié)構(gòu)模型。該模型可以模擬致密巖石的單軸壓縮實(shí)驗(yàn)，還可以監(jiān)測(cè)不同脆性致密巖石在壓縮過程中的裂紋數(shù)量、裂紋分布以及裂紋密度在壓縮過程中的變化特征。這些特征有助于我們更好地理解致密巖石的微觀結(jié)構(gòu)破壞機(jī)制，從而更好地理解致密巖石的宏觀破壞機(jī)制。

1 UDEC–Voronoi 模型

本研究采用的是 UDEC–Voronoi 建模方法，該方法利用 UDEC 內(nèi)建的 Voronoi 鑲嵌模式的自動(dòng)生成器，將一個(gè)特定區(qū)域的模型細(xì)分為隨機(jī)大小的多邊形。通過這種方式，多邊形可以代表顆粒、顆粒組合或完整巖石中隨機(jī)擾動(dòng)樣本的其他缺陷[3,10](圖 1)。

1.1 運(yùn)動(dòng)方程

根據(jù)式(1) 中的新速度確定新的顆粒位置為

式中，xi為顆粒坐標(biāo)位置，θ為顆粒繞質(zhì)心轉(zhuǎn)動(dòng)的角度。

在UDEC 的計(jì)算過程中，每個(gè)時(shí)間步長(zhǎng)都會(huì)產(chǎn)生新的離散塊體位置，從而產(chǎn)生新的塊體之間的接觸力，利用作用在離散塊體上的合力和合力矩計(jì)算各塊體的直線加速度和角加速度，離散塊體的速度和位移則通過對(duì)時(shí)間增量的積分來確定。這個(gè)計(jì)算過程不斷重復(fù)，直至達(dá)到離散塊體系統(tǒng)的平衡狀態(tài)或出現(xiàn)一個(gè)持續(xù)的失效結(jié)果。

1.2 本構(gòu)定律

在低應(yīng)力條件下巖石內(nèi)部的破裂主要為沿晶破裂，且Voronoi 模型中顆粒不會(huì)發(fā)生破壞，破壞只沿著顆粒邊界發(fā)生，因此在本研究中顆粒本構(gòu)模型選用各向同性彈性本構(gòu)模型。各向同性彈性本構(gòu)模型具有卸載時(shí)可逆變形的特征，應(yīng)力應(yīng)變定律是線性的，與路徑無關(guān)。顆粒模型的變形特征由體積模量和剪切模量表征。接觸的本構(gòu)模型選取不含殘余強(qiáng)度的庫侖滑移?面接觸模型，如圖2 所示。接觸模型的變形特征由其法向剛度kn和切向剛度ks表征，強(qiáng)度特征由其摩擦角φ、黏聚力c和抗拉強(qiáng)度T來表征。圖 2 中 ?τs為接觸面切向應(yīng)力改變量，為切向位移增量的彈性部分，?us為總切向位移增量。

2 致密巖石的微觀力學(xué)行為模擬

圖2 離散塊體接觸本構(gòu)[21]

本研究選取了 4 類典型的頁巖巖樣進(jìn)行能量色散 X 射線光譜學(xué)分析，獲取巖樣的礦物組分，并采用 Jin 等[22]提出的基于礦物成分的脆性指數(shù)(式(3)) 計(jì)算了試樣的脆性指數(shù)，整理成表1。接著按礦物組成隨機(jī)分布4 組礦物顆粒，分別代表致密巖石的4 種礦物組分并賦予不同的彈性參數(shù)，不同礦物之間的接觸被分別賦予不同的微觀力學(xué)參數(shù)，受計(jì)算機(jī)能力的限制，構(gòu)建的 Voronoi 模型的平均粒徑為 1.5 mm (圖 3)。

表1 不同脆性指數(shù)巖樣地化參數(shù)

圖 3 構(gòu)建的 Voronoi 模型

式中，F(xiàn)Qtz+Fsp+Cb為石英、長(zhǎng)石和方解石的礦物含量總和，F(xiàn)total為總體礦物含量，該值為100%。

通過單軸壓縮模擬實(shí)驗(yàn)與巖樣的單軸壓縮實(shí)驗(yàn)比對(duì)(圖4 中橫軸正應(yīng)變?yōu)檩S向應(yīng)變，負(fù)應(yīng)變?yōu)榄h(huán)向應(yīng)變；下同)，校準(zhǔn)模型的細(xì)觀力學(xué)參數(shù)，最終得到模型的細(xì)觀力學(xué)參數(shù)如表2 和表3 所示。

圖4 巖樣和模型的單軸壓縮應(yīng)力?應(yīng)變曲線

2.1 不同脆性巖石的壓縮破壞特征

本研究接著構(gòu)建了不同脆性巖石的單軸壓縮模擬實(shí)驗(yàn)，4 個(gè)試樣分別采取相同的加載條件，獲取的應(yīng)力?應(yīng)變曲線如圖5 所示。隨著脆性的降低，巖石的單軸抗壓強(qiáng)度從200 MPa 下降至150 MPa，彈性模量從34.5 GPa 減少至8.87 GPa，同時(shí)顯現(xiàn)出更強(qiáng)的延性特征，峰前應(yīng)力?應(yīng)變曲線也從彈性類型轉(zhuǎn)變?yōu)閺?塑性類型，應(yīng)力?應(yīng)變曲線的向下彎曲度不斷提升，巖石的應(yīng)變值也顯著增加了。

表2 礦物彈性參數(shù)[23-24]

4 種模型的破壞形態(tài)分別如圖6(a)～圖6(d) 所示。巖石出現(xiàn)的破壞主要為拉伸破壞，局部也可見一些剪切破壞。相較于A1 和 A2 高脆性巖石，A3 和A4 低脆性巖石破壞時(shí)部分區(qū)域出現(xiàn)滑移趨勢(shì)，且試樣的環(huán)向變形較大，側(cè)幫和底部出現(xiàn)了較大程度的滑移和脫落。

表3 不同脆性指數(shù)模型礦物間接觸力學(xué)參數(shù)

2.2 不同脆性巖石的裂紋數(shù)量變化

為了模擬不同脆性巖石在壓縮過程中的裂紋擴(kuò)展情況，本研究利用 FISH 編寫了監(jiān)測(cè) Voronoi 模型內(nèi)部裂紋變化信息的功能，并選取不同的參量進(jìn)行分析。由于模型在壓縮過程中主要發(fā)生拉伸破壞，故主要監(jiān)測(cè)拉伸裂紋的數(shù)量變化。在軸向力不斷施加的過程中，巖石的微裂紋數(shù)量變化情況如圖 7～圖 10 所示。

壓縮過程中微裂紋的數(shù)量呈現(xiàn)出分階段的變化特點(diǎn)，根據(jù)微裂紋數(shù)量的變化特征，選取裂紋數(shù)量變化曲線上切線斜率水平達(dá)到2500 的特征點(diǎn)：裂紋起裂點(diǎn)的應(yīng)力為σci，表征巖石內(nèi)部開始出現(xiàn)微裂紋；選取切線斜率水平達(dá)到7500 的特征點(diǎn)：裂紋貫通點(diǎn)的應(yīng)力為σcc，表征巖石內(nèi)部微裂紋開始貫通，微裂紋急劇增長(zhǎng)；σp為巖石的峰值強(qiáng)度。

圖5 不同脆性巖石的單軸壓縮應(yīng)力?應(yīng)變曲線

圖6 巖石的破壞形態(tài)(紅線為拉伸裂紋，藍(lán)線為剪切裂紋)

圖7 A1 試樣的微拉伸裂紋變化

圖8 A2 試樣的微拉伸裂紋變化

圖9 A3 試樣的微拉伸裂紋變化

圖10 A4 試樣的微拉伸裂紋變化

不同脆性巖石的裂紋起裂點(diǎn)和裂紋貫通點(diǎn)顯現(xiàn)出不同的特征。首先分析裂紋起裂點(diǎn)。以各試樣的裂紋起裂點(diǎn)應(yīng)力與峰值強(qiáng)度的百分比值作為衡量標(biāo)準(zhǔn)，脆性巖石的裂紋起裂點(diǎn)應(yīng)力水平較低脆性巖石高約10%σp，裂紋數(shù)量則較低脆性巖石更少。其次分析裂紋貫通點(diǎn)。脆性巖石在峰前階段存在明顯的裂紋貫通點(diǎn)，而低脆性巖石微裂紋數(shù)量保持較為緩慢的增長(zhǎng)趨勢(shì)直到試樣破壞，其裂紋貫通點(diǎn)與峰值應(yīng)力點(diǎn)幾乎重合，不存在明顯的裂紋貫通點(diǎn)。此外，在峰前階段，脆性巖石的微裂紋數(shù)量增加速度要明顯低于低脆性巖石；到達(dá)應(yīng)力峰值點(diǎn)時(shí)，低脆性巖石的微裂紋數(shù)量大約是高脆性巖石的1.5 倍，表明低脆性巖石需要大量的裂紋累積才能形成宏觀貫通裂紋破壞。

2.3 不同脆性巖石的裂紋分布變化

為了直觀地觀察壓縮破壞過程中微裂紋的擴(kuò)展過程，繪制了巖石模型在不同加載階段的微裂紋分布圖(圖11)，需要說明的是，圖中的裂紋是該階段的某一時(shí)刻的裂紋，所以存在加載初期的剪切裂紋隨著加載的進(jìn)行向拉伸裂紋轉(zhuǎn)化的情況。可以發(fā)現(xiàn)，微拉伸裂紋傾向于與試樣加載方向平行，而微剪切裂紋則與試樣加載方向成不同角度。在裂紋起裂階段，微拉伸裂紋和剪切裂紋的數(shù)量均較少，且裂紋的長(zhǎng)度較短；到達(dá)裂紋貫通階段后，裂紋的長(zhǎng)度得到了顯著的增長(zhǎng)，裂紋開始相互貫通，短裂紋的數(shù)量明顯減少，此時(shí)剪切裂紋為主導(dǎo)型裂紋；到達(dá)峰值應(yīng)力階段后，微剪切裂紋的比例開始降低，巖石開始出現(xiàn)破壞，此時(shí)拉伸裂紋成為主導(dǎo)型裂紋；到達(dá)峰后階段，拉伸裂紋進(jìn)一步增加，在拉伸、剪切裂紋共同的作用下，試樣形成了更大尺度的破壞。

圖11 A2 試樣各階段裂紋分布(紅線為拉伸裂紋；藍(lán)線為剪切裂紋)

2.4 不同脆性巖石的裂紋密度變化

為了更好地理解巖石破裂與裂紋擴(kuò)展之間的聯(lián)系，根據(jù) Bristow 提出的有關(guān)裂紋密度的定義[25](式4)，將裂紋密度作為衡量巖石內(nèi)部裂紋擴(kuò)展程度的指標(biāo)。為了更好地探究不同脆性巖石的裂紋密度變化與起裂應(yīng)力、峰值應(yīng)力之間的規(guī)律，增加了脆性指數(shù)分別為 0.86，0.74 和 0.5 的數(shù)值試件，其地化參數(shù)如表4 所示，增加的數(shù)值試件的構(gòu)建過程與A1～A4 試樣完全一致。

式中，C為拉裂紋密度，Li為統(tǒng)計(jì)區(qū)域內(nèi)第i條裂紋的長(zhǎng)度，S為統(tǒng)計(jì)區(qū)域的面積。

表4 增加的數(shù)值試件的地化參數(shù)

圖12 展示了不同脆性巖石的起裂應(yīng)力和起裂點(diǎn)裂紋密度，從圖中可以看出，隨著巖石脆性指數(shù)的增加，巖石的起裂應(yīng)力呈現(xiàn)增加趨勢(shì)，起裂點(diǎn)的裂紋密度呈現(xiàn)減少趨勢(shì)，高脆性的A4 試樣的起裂點(diǎn)裂紋密度約為低脆性的A1 試樣的2 倍，而A1 試樣的起裂應(yīng)力約為A4 試樣的1.7 倍，可見高脆性巖石是在較高的應(yīng)力狀態(tài)下和較低的裂紋損傷狀態(tài)下達(dá)到裂紋起裂點(diǎn)。換句話說，相較于低脆性巖石，高脆性巖石會(huì)產(chǎn)生更多高能量的裂縫[26]。

圖12 不同脆性巖石的起裂點(diǎn)應(yīng)力和裂紋密度

圖13 展示了不同脆性巖石的峰值應(yīng)力和裂紋密度，從圖中可以看出，致密巖石的峰值應(yīng)力仍然符合隨著脆性指數(shù)的增加而增加的規(guī)律，A1 試樣的峰值強(qiáng)度約為A4 試樣的1.3 倍。但裂紋密度隨著脆性指數(shù)的增加呈現(xiàn)先增加后減小的趨勢(shì)。經(jīng)過數(shù)據(jù)擬合，峰值點(diǎn)裂紋密度和脆性指數(shù)符合

式中，C為裂紋密度，B為脆性指數(shù)。

圖13 不同脆性巖石的峰值點(diǎn)應(yīng)力和裂紋密度

在實(shí)際的巖樣測(cè)試中，試樣的脆性指數(shù)往往是容易測(cè)試的，然后根據(jù)經(jīng)驗(yàn)公式 (5) 可以對(duì)試樣在峰值應(yīng)力處的裂紋損傷水平進(jìn)行估計(jì)。另外，從圖13 中可以看出，當(dāng)脆性指數(shù)約為0.6 時(shí)，峰值點(diǎn)裂紋密度和脆性指數(shù)的變化曲線出現(xiàn)拐點(diǎn)，可知在脆性指數(shù)小于 0.6 時(shí)巖石的性質(zhì)出現(xiàn)了改變，該脆性指數(shù)可以作為區(qū)分頁巖脆性的經(jīng)驗(yàn)值。

3 結(jié)論和討論

本研究介紹了UDEC–Voronoi 模型的構(gòu)建方法，建立了致密巖石的微結(jié)構(gòu)模型，進(jìn)行了不同脆性巖石的微觀力學(xué)行為模擬。結(jié)果顯示，通過微結(jié)構(gòu)參數(shù)的校準(zhǔn)，UDEC–Voronoi 模型能夠很好地模擬巖石力學(xué)特性曲線，表征巖石的變形和強(qiáng)度特性；在模擬致密巖石的微觀結(jié)構(gòu)，監(jiān)測(cè)不同脆性巖石微觀裂紋的生成與擴(kuò)展過程，探究致密巖石的微觀破裂機(jī)制方面，該模型有良好的應(yīng)用潛力。通過數(shù)值模擬獲得的研究結(jié)論主要有：

(1)在壓縮過程中，致密巖石內(nèi)部的拉伸裂紋逐漸增加，最終成為主導(dǎo)性裂紋，導(dǎo)致巖石的破壞形態(tài)多為拉伸破壞。高脆性巖石的單軸抗壓強(qiáng)度要更高，具有更高的彈性模量；而低脆性巖石的塑性特征更加明顯，峰前應(yīng)力?應(yīng)變曲線從彈性類型轉(zhuǎn)換為彈?塑性類型。

(2)隨著巖石脆性指數(shù)的增加，巖石的起裂應(yīng)力呈現(xiàn)增加趨勢(shì)，而起裂點(diǎn)的裂紋密度呈現(xiàn)減少趨勢(shì)，高脆性巖石是在較高的應(yīng)力狀態(tài)下和較低的裂紋損傷狀態(tài)下達(dá)到裂紋起裂點(diǎn)。

(3) 巖石的峰值應(yīng)力隨著脆性指數(shù)的增加而增加，峰值點(diǎn)裂紋密度和脆性指數(shù)符合一擬合公式，根據(jù)該公式可以在已知巖石的脆性指數(shù)的情況下對(duì)峰值應(yīng)力處巖石的裂紋損傷程度進(jìn)行估計(jì)?？梢园汛嘈灾笖?shù)0.6 作為區(qū)分頁巖脆性的經(jīng)驗(yàn)值。

必須指出，本研究介紹的 UDEC–Voronoi 模型還存在一些限制，主要有以下幾點(diǎn)：

(1) 構(gòu)建的巖石微結(jié)構(gòu)模型中含有石英、長(zhǎng)石、伊利石和方解石四種礦物顆粒，在計(jì)算模擬中只考慮了致密巖石的一般性特征，并未考慮巖石學(xué)、礦物學(xué)和微觀結(jié)構(gòu)等細(xì)微特征，且未考慮四種礦物顆粒大小的不同。

(2)由于計(jì)算機(jī)能力的限制，構(gòu)建的巖石模型的粒徑并不等于實(shí)際的巖石顆粒尺寸，而巖石顆粒尺寸對(duì)巖石的破裂機(jī)制具有一定的影響。隨著計(jì)算機(jī)能力的發(fā)展，模擬實(shí)際巖石的顆粒尺寸將有可能實(shí)現(xiàn)。

(3) UDEC–Voronoi 模型為二維模型，未來的研究將致力于把Voronoi 模型擴(kuò)展到三維的巖石模型中，這將使得構(gòu)建的微結(jié)構(gòu)模型更加接近于巖石材料的實(shí)際狀況。