基于Voronoi 模型的脆性巖石微觀結構數值模擬1)

邱 鑫 林 緬, 鄭思平 陳天宇

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

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

??(東北大學深部金屬礦山安全開采教育部重點實驗室，沈陽110819)

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

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

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

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

1 UDEC–Voronoi 模型

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

1.1 運動方程

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

式中，xi為顆粒坐標位置，θ為顆粒繞質心轉動的角度。

在UDEC 的計算過程中，每個時間步長都會產生新的離散塊體位置，從而產生新的塊體之間的接觸力，利用作用在離散塊體上的合力和合力矩計算各塊體的直線加速度和角加速度，離散塊體的速度和位移則通過對時間增量的積分來確定。這個計算過程不斷重復，直至達到離散塊體系統的平衡狀態或出現一個持續的失效結果。

1.2 本構定律

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

2 致密巖石的微觀力學行為模擬

圖2 離散塊體接觸本構[21]

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

表1 不同脆性指數巖樣地化參數

圖 3 構建的 Voronoi 模型

式中，FQtz+Fsp+Cb為石英、長石和方解石的礦物含量總和，Ftotal為總體礦物含量，該值為100%。

通過單軸壓縮模擬實驗與巖樣的單軸壓縮實驗比對(圖4 中橫軸正應變為軸向應變，負應變為環向應變；下同)，校準模型的細觀力學參數，最終得到模型的細觀力學參數如表2 和表3 所示。

圖4 巖樣和模型的單軸壓縮應力?應變曲線

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

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

表2 礦物彈性參數[23-24]

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

表3 不同脆性指數模型礦物間接觸力學參數

2.2 不同脆性巖石的裂紋數量變化

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

壓縮過程中微裂紋的數量呈現出分階段的變化特點，根據微裂紋數量的變化特征，選取裂紋數量變化曲線上切線斜率水平達到2500 的特征點：裂紋起裂點的應力為σci，表征巖石內部開始出現微裂紋；選取切線斜率水平達到7500 的特征點：裂紋貫通點的應力為σcc，表征巖石內部微裂紋開始貫通，微裂紋急劇增長；σp為巖石的峰值強度。

圖5 不同脆性巖石的單軸壓縮應力?應變曲線

圖6 巖石的破壞形態(紅線為拉伸裂紋，藍線為剪切裂紋)

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

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

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

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

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

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

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

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

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

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

式中，C為拉裂紋密度，Li為統計區域內第i條裂紋的長度，S為統計區域的面積。

表4 增加的數值試件的地化參數

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

圖12 不同脆性巖石的起裂點應力和裂紋密度

圖13 展示了不同脆性巖石的峰值應力和裂紋密度，從圖中可以看出，致密巖石的峰值應力仍然符合隨著脆性指數的增加而增加的規律，A1 試樣的峰值強度約為A4 試樣的1.3 倍。但裂紋密度隨著脆性指數的增加呈現先增加后減小的趨勢。經過數據擬合，峰值點裂紋密度和脆性指數符合

式中，C為裂紋密度，B為脆性指數。

圖13 不同脆性巖石的峰值點應力和裂紋密度

在實際的巖樣測試中，試樣的脆性指數往往是容易測試的，然后根據經驗公式 (5) 可以對試樣在峰值應力處的裂紋損傷水平進行估計。另外，從圖13 中可以看出，當脆性指數約為0.6 時，峰值點裂紋密度和脆性指數的變化曲線出現拐點，可知在脆性指數小于 0.6 時巖石的性質出現了改變，該脆性指數可以作為區分頁巖脆性的經驗值。

3 結論和討論

本研究介紹了UDEC–Voronoi 模型的構建方法，建立了致密巖石的微結構模型，進行了不同脆性巖石的微觀力學行為模擬。結果顯示，通過微結構參數的校準，UDEC–Voronoi 模型能夠很好地模擬巖石力學特性曲線，表征巖石的變形和強度特性；在模擬致密巖石的微觀結構，監測不同脆性巖石微觀裂紋的生成與擴展過程，探究致密巖石的微觀破裂機制方面，該模型有良好的應用潛力。通過數值模擬獲得的研究結論主要有：

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

(2)隨著巖石脆性指數的增加，巖石的起裂應力呈現增加趨勢，而起裂點的裂紋密度呈現減少趨勢，高脆性巖石是在較高的應力狀態下和較低的裂紋損傷狀態下達到裂紋起裂點。

(3) 巖石的峰值應力隨著脆性指數的增加而增加，峰值點裂紋密度和脆性指數符合一擬合公式，根據該公式可以在已知巖石的脆性指數的情況下對峰值應力處巖石的裂紋損傷程度進行估計。可以把脆性指數0.6 作為區分頁巖脆性的經驗值。

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

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

(2)由于計算機能力的限制，構建的巖石模型的粒徑并不等于實際的巖石顆粒尺寸，而巖石顆粒尺寸對巖石的破裂機制具有一定的影響。隨著計算機能力的發展，模擬實際巖石的顆粒尺寸將有可能實現。

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