基于Voronoi塊體模型的砂巖力學特性研究

孫 浩, 陳帥軍, 金愛兵, 李 海

(1. 北京科技大學 金屬礦山高效開采與安全教育部重點實驗室, 北京 100083; 2. 北京科技大學 土木與資源工程學院, 北京 100083)

隧道、邊坡、橋梁和硐室等各種工程施工中遇到的巖石具有不同的力學行為,而針對巖石力學行為的研究是解決各種工程問題的基礎.不同種類巖石力學行為的差異不僅體現在其礦物成分和膠結方式,巖石中各類微孔隙和礦物結晶亦嚴重影響巖石的力學特性[1].

物理試驗和數值模擬一直是巖石力學特性研究的常用手段.物理試驗中,真實巖樣[2]和相似材料試樣[3]被廣泛應用于巖石力學性質及結構面性質的研究中;而最新興起的3D打印試樣技術[4]為巖石力學行為研究注入了新的活力.與數值模擬相比,物理試驗具有能反映研究對象物理力學特性的優點,但也存在控制變量困難、可重復性差以及細觀機理研究困難等弊端.數值模擬研究中常用的方法主要包括:有限單元法(finite element method,FEM)[5]、有限差分法(finite difference method,FDM)[6]、邊界元法(boundary element method,BEM)[7]和離散單元法(discrete element method,DEM)[1]等.其中,DEM因其能夠直觀反映試樣內裂紋擴展及貫通過程而被廣泛應用.目前,應用于試樣尺度巖石力學細觀研究的離散元軟件有YADE[8],ESyS-Particle[9]和PFC(particle flow code)[10]等,均采用球形或球形簇作為基質材料;UDEC(universal distinct element code)[11]軟件雖然可生成不規則剛性塊體,但更適用于邊坡和地下工程開挖等大變形問題研究.PFC是由Cundall等[12]提出的用于模擬無黏結顆粒材料的細觀離散元數值軟件,而在Potyondy等[13]提出黏結模型后,PFC開始逐漸應用于各類巖石力學研究.

利用PFC軟件進行試樣尺度巖石力學數值試驗研究一般分為以下三個步驟:① 建模,建立基于顆粒體系的巖石試樣模型;② 接觸模型選擇及細觀參數匹配,選擇接觸模型,并對墻體、顆粒自身參數以及顆粒-顆粒、顆粒-墻體間接觸模型參數的匹配研究;③ 數值試驗,開展單軸壓縮、三軸壓縮和直接剪切等各類巖石力學數值試驗研究.上述三步驟模擬的準確性均會影響所研究結果的可靠性.目前,針對步驟②的接觸模型選擇及細觀參數確定問題,眾多學者已分別采用平行黏結模型(parallel bond model,PBM)[10]、平節理模型(flat-joint model,FJM)[2]和伯格斯模型(Burger’s model)[14]等接觸模型開展了各類巖石的物理力學特性研究;針對步驟③,相關學者也已開展了PFC單軸壓縮試驗[15]、三軸壓縮與卸荷試驗[2]、巴西劈裂試驗[16]以及巖石動荷載試驗[17]等各類力學試驗的數值模擬研究.本文以山西呂梁地區典型砂巖試樣為依據,從數值模型建立(步驟①)和接觸模型選擇及細觀參數匹配(步驟②)兩個方面對現行數值模擬方法進行優化.

首先,在數值模型構建方面(步驟①),國內外絕大多數學者均采用球形顆粒或球形簇進行巖石力學問題研究.而相關室內試驗研究已表明:礦物組成、粒徑尺寸、顆粒形狀等因素均對巖石強度具有顯著影響[18].含球形和球形簇顆粒模型無法精確模擬巖石結構中的礦物晶粒形狀,致使目前基于PFC的巖石力學問題研究中所構建的數值模型準確性較差.Potyondy等[19]提出了基于PFC2D的細觀晶粒模型(grain-based model,GBM),研究巖石晶粒的非均勻性對其強度和裂紋擴展模式的影響,但其本質還是采用球形顆粒團簇的方式實現巖石晶粒的模擬,且晶粒破裂后仍是一個個獨立的球形顆粒,球形顆粒間較小的接觸面積與真實巖石晶粒間接觸差異較大,同時,球形顆粒團簇又大幅增加了模型內顆粒數目,降低了計算效率.目前,數值模擬中建模階段(步驟①)所得數值模型與真實巖石試樣誤差較大,且接觸模型選擇及細觀參數匹配(步驟②)的準確性對試樣力學性質具有重要影響.因此,有必要對巖石力學PFC模擬的建模階段開展優化研究.

在接觸模型選擇及細觀參數匹配(步驟②)方面,與連續介質方法直接輸入室內試驗所得宏觀物理力學參數不同,PFC中宏-細觀參數之間無定量關系.因此,宏-細觀參數的匹配研究是PFC模擬的前提.研究表明: PBM可以模擬巖石的彈性模量、泊松比和無側限抗壓強度[20].但球形顆粒所提供的較小配位數導致球形顆粒間自鎖效應偏低,且基于球形顆粒的PBM不能提供合適的轉動阻力[21].因此,模擬所得無側限抗壓強度與抗拉強度的比值(即壓拉比)僅能在3～7范圍內波動[13],小于真實巖石的壓拉比(10～20)[22].Scholtès等[23]認為較低壓拉比的數值模型在模擬同時存在壓縮和拉伸應力路徑的問題時誤差較大.

綜上所述,在采用球形顆粒和PBM進行巖石力學相關研究時存在兩個明顯缺陷:① 球形顆粒與巖石材料的晶粒形狀相差較大;② 無側限抗壓強度與抗拉強度之比(壓拉比)較低.基于此,本文以山西呂梁地區典型砂巖為例,通過室內單軸壓縮試驗和掃描電鏡試驗獲得砂巖的應力-應變關系及微觀晶粒圖像,基于Voronoi剖分技術構建剛性塊體數值試樣,并利用軟化黏結模型(soft-bond model,SBM)開展無側限壓縮和直接拉伸數值模擬研究,并與傳統基于球形顆粒和PBM的數值模擬方法對比,分析剛性塊體數值試樣和SBM在模擬典型砂巖晶粒形狀不規則性和大壓拉比等方面的優勢,研究方法和結果能夠提高細觀離散元軟件PFC在砂巖力學行為數值模擬研究中的適用性與可靠性,為其他晶粒型結構巖石在多應力路徑下的模擬研究提供借鑒.

1 PFC中的Voronoi剖分

1.1 3D Voronoi剖分原理

3D Voronoi剖分是將三維空間內的某個區域劃分為一系列的單元體[24].給定歐幾里德空間內有限個點P={p1,p2,…,pN},并且P∈R3,空間中的每一個點均隸屬于一個Voronoi單元體Vi,Vi定義如下:

Vi={xi∈‖x-xi‖≤‖x-xj‖,?i≠j} .

(1)

其中,‖x-xi‖為歐幾里德距離[25].空間中各點隨機分布,每一個Voronoi單元體均可視作空間中隨機點的影響范圍,而每一個Voronoi單元體均由隨機點P生成.上述算法計算所得每個Voronoi單元體的分界面為兩個空間相鄰點連線的中垂面,因此并未考慮Voronoi單元體的尺寸變化.通過引入Voronoi多邊形的概念,對Voronoi單元體進行正則化,即可生成不同粒徑分布的Voronoi單元體[26].

1.2 PFC中剛性塊體模型構建

在PFC中利用Voronoi剖分原理生成標準試樣的步驟如圖1所示.首先生成圓柱形墻體;其次在圓柱形墻體內生成一定孔隙率和粒徑分布的球形顆粒,并進行靜力平衡計算;最后使用Voronoi剖分命令,以球心坐標為隨機點,以球半徑為歐幾里德距離進行Voronoi剖分和正則化,得到數值試樣.采用TESCAN VEGA3掃描電子顯微鏡(圖2)獲得的砂巖試樣SEM圖像如圖3a所示,球形顆粒如圖3b所示,基于Voronoi剖分生成試樣的局部放大圖如圖3b所示.從圖3中可以看出,通過Voronoi剖分所建立的塊體模型顆粒大小不一、形狀各異,相較于球形顆粒,更符合砂巖中礦物晶粒分布.

圖1 基于Voronoi剖分的PFC剛性塊體建模過程

圖2 TESCAN VEGA3掃描電子顯微鏡

圖3 砂巖晶粒結構和兩種數值試樣顆粒局部圖

研究表明[21],球形顆粒間的配位數較小,導致球形顆粒間的自鎖效應偏低,而基于Voronoi剖分原理建立的剛性塊體則有效解決了這一難題.以典型四面體剛性塊體顆粒為例,如圖4所示,在理想情況下,四面體剛性塊體中4個面、6條邊和4個頂點均能生成接觸,且每條邊和每個頂點又能同時生成多個接觸,而球形顆粒僅能生成4個接觸(根據Voronoi剖分原理,生成Voronoi多面體的面數等于原始球形顆粒的接觸個數),剛性塊體間的配位數高于球形顆粒,導致剛性塊體間的自鎖效應明顯提高.基于Voronoi塊體的建模方法是對PFC傳統建模方法(Ball,GBM以及Clump)的有力補充,塊體之間的運動和相互作用力均遵循牛頓第二定律.

圖4 四面體剛性塊體

2 基于SBM的剛性塊體試樣數值模擬

2.1 數值模型構建

本次研究基于SBM建立剛性塊體試樣模型.在以往研究中,球形和球形簇一直作為模擬巖石的基質材料,通過對墻體施加一定速度的方式對整個試樣施加壓力,模擬室內試驗中無側限壓縮過程.由于剛性塊體之間緊密接觸,孔隙率較小,軟件的算法誤差會造成顆粒的局部凸起(圖5),使試件的上、下表面不平整,而凸起的塊體在用墻體或其他平整面作為加載端時會形成局部的應力集中,造成試樣上下端部的局部破壞,使模擬結果產生較大誤差.因此,本次數值試驗通過建立高度為110 mm的數值試樣,設置上下5 mm部分剛性塊體分別為上下加載端[27].由于加載端和試樣一次建模成型,試樣和加載端結合較好,通過上下加載端勻速運動對試樣施加壓力的方式模擬無側限壓縮過程.該方法充分考慮了試樣和加載端的接觸條件,消除了由于上下端面的不平整帶來的誤差.剛性塊體數值模型試樣如圖6所示.

圖5 數值試樣中剛性塊體的局部凸起

圖6 剛性塊體數值模型試樣

2.2 SBM原理簡介

SBM是PFC中開發的接觸黏結模型,可以應用于黏結和未黏結的顆粒[27],本文主要介紹和使用黏結顆粒.軟化黏結模型黏結部分的接觸行為如圖7所示,其中ζ為軟化系數,γ為軟化抗拉強度系數,σc為抗拉強度,kn為法向剛度,kt為切向剛度,c,φ為內聚力和內摩擦角.

圖7 軟化黏結模型顆粒間黏結狀態和流變組件[27]

默認情況下(即軟化系數ζ=0), SBM類似于PBM,而當軟化抗拉強度系數γ=1.0時,黏結鍵的軟化行為也將受到抑制.在ζ≠0和γ≠1.0時黏結鍵會發生拉伸軟化.SBM的軟化行為只對拉伸響應起作用,而對剪切響應無效,因此,在抗剪強度超過極限之后黏結鍵失效破裂.一旦黏結鍵的拉力超過了黏結鍵的抗拉強度,黏結鍵間的拉力不會像PBM中那樣自動設置為0,而是不斷減小,當拉力小于γ和峰值抗拉強度的乘積時,黏結鍵破裂失效.軟化速率由軟化系數ζ決定,如果ζ=1,拉伸軟化以與拉力達到峰值抗拉強度相同的絕對速率進行(速率符號相反)[27].

2.3 基于剛性塊體的SBM細觀參數確定

以室內試驗力學參數為對照數據,通過宏-細觀參數匹配研究,獲得可靠的軟化黏結接觸模型PFC細觀參數組合.室內試驗砂巖巖樣取自山西呂梁地區,加工規格為φ50 mm×100 mm,采用YAW-600微機控制電液伺服剛性壓力試驗機(圖8)進行試驗獲得砂巖試樣應力-應變關系.通過構建基于球形顆粒和剛性塊體模型的數值試樣,對比探究剛性塊體模型模擬砂巖力學特性的可行性.

圖8 YAW-600單軸壓力試驗機

利用PFC3D軟件建立高為110 mm(上下各5 mm作為加載端),直徑φ=50 mm的三維數值模型,Voronoi剖分前球體顆粒直徑為2.0～3.32 mm,生成剛性塊體總數目為16 142,剛性塊體間采用SBM.與此同時,建立高度為110 mm(上下各5 mm作為加載端),直徑φ=50 mm的標準圓柱形試樣,試樣內球形顆粒尺寸為2.0～3.32 mm,生成球形顆粒總數目為14 959,顆粒間采用傳統的PBM.通過不斷調整顆粒間黏結的細觀參數,匹配室內試驗中砂巖試樣的彈性模量、峰值強度、峰值應變和泊松比等參數.數值模擬軸向加載速率為0.05 m/s,終止加載條件為峰值強度的80%.確定的SBM和PBM細觀參數如表1所示.

表1 SBM和PBM細觀參數匯總表

真實巖石試樣中存在眾多微孔隙,而數值模擬試樣中缺少可壓縮的微孔隙,導致離散元數值模擬方法仍無法模擬真實巖石中的微孔隙壓密階段[28].因此,數值模擬試驗中試樣的彈性模量略小于室內試驗中的彈性模量.球形顆粒和PBM已經在巖石力學領域廣泛應用,并已被證實可以模擬砂巖試樣的應力-應變曲線，兩種數值模型的應力-應變曲線如圖9所示.從圖9中可以看出,本文采用的SBM和剛性塊體在模擬砂巖應力-應變曲線方面與PBM和球形顆粒并無差異,即證明了本文中采用SBM和剛性塊體數值模擬方法的可靠性.數值模擬與室內試驗所得宏觀力學參數如表2所示.由表2可知，數值模擬所得試樣的峰值強度、峰值應變、彈性模量和泊松比等宏觀力學參數均與室內試驗所得結果較為接近.因此,認為采用基于剛性塊體試樣的SBM可以用于砂巖力學行為的數值模擬試驗中.

表2 室內試驗和數值模擬所得試樣力學參數

圖9 室內試驗和數值模擬所得試樣單軸壓縮應力-應變關系

2.4 兩種數值模型試樣無側限壓縮破壞模式對比分析

本文從數值模型構建和接觸模型選擇兩個方面對數值模擬過程進行優化,采用的剛性塊體模型與真實砂巖顆粒形狀較為接近,剛性塊體間的細觀接觸行為與真實顆粒間的接觸行為近似,且顆粒間配位數較大,顆粒自鎖效應較大.因此,宏觀數值試樣的破壞模式與真實砂巖試樣的破壞模式較為接近.圖10為真實砂巖試樣和兩種數值試樣的破壞模式:剛性塊體試樣破壞模式為經典巖石力學中常見的剪切破壞,更符合真實砂巖的破裂模式;而球形顆粒模型易產生局部破壞,從而為巖石破裂模式的匹配帶來一定困難,這也是學者在利用PFC3D進行巖石力學研究中對完整巖石破裂模式匹配較少的重要原因.因此,從試樣破裂模式可以明顯發現剛性塊體模型優于傳統球形顆粒模型.

圖10 無側限壓縮試樣破壞圖

2.5 兩種數值模型試樣拉伸試驗模擬結果對比分析

本文所采用的軟化黏結接觸模型通過引入軟化系數ζ來改善黏結鍵的力學行為,軟化系數不僅影響無側限壓縮數值模擬結果,亦對試樣抗拉強度影響較大.圖11為采用剛性塊體與SBM數值試樣和球形顆粒與PBM數值試樣的直接拉伸的應力-應變曲線,采用球形顆粒與PBM的數值試樣抗拉強度為10.53 MPa,采用剛性塊體與SBM的數值試樣抗拉強度為4.46 MPa,從表2可知,球形顆粒與PBM數值試樣和剛性塊體與SBM數值試樣的無側限抗壓強度分別為53.23 MPa 和52.78 MPa.無側限抗壓強度相等情況下,基于SBM的剛性塊體試樣抗拉強度較低.

圖11 兩種數值試樣抗拉應力-應變曲線

通過上述對比分析,從力學參數和破壞模式等角度探究了基于剛性塊體試樣的軟化黏結接觸模型的優越性,為后續巖石力學數值模擬研究提供一種新型建模方法.

3 基于SBM的剛性塊體數值模型研究

3.1 基于SBM模型的巖石壓拉比分析

通常情況下,巖石的壓拉比在10～20之間[22],表3為常見的三種典型巖石的無側限抗壓強度和抗拉強度.為了獲得兩種模型試樣的抗拉強度,分別針對三種典型巖石進行兩種模型數值模擬拉伸試驗,每種巖石試樣進行2組拉伸數值模擬,共計6組數值拉伸試驗.

表3 三種巖石單軸抗壓、抗拉強度[29-31]

傳統的PBM模型所能達到的壓拉比僅僅為3～7[13].SBM模型通過在PBM模型基礎上引入軟化系數,可以明顯提高數值試樣的壓拉比.表4為三種巖石PBM和SBM模型的無側限抗壓強度和抗拉強度數值模擬數據.從表4中可以看出,SBM相較于PBM來說,壓拉比明顯增大,更趨近于典型巖石的壓拉比,從而提高模擬結果的準確性.

表4 三種巖石的壓拉比數值模擬結果

3.2 SBM和PBM壓拉比影響機制分析

巖石類材料的抗壓強度遠大于其抗拉強度,通常情況下,典型巖石的壓拉比在10～20之間[22].對于PBM而言,其黏結鍵的抗拉和抗剪強度在黏結鍵被移除之前始終保持不變,其黏結鍵的抗拉力學行為如圖12a所示(抗剪強度類似).隨著顆粒間拉力達到黏結鍵的抗拉強度,黏結鍵破裂并被移除,此時顆粒間表現為無黏結行為.對于采用PBM模型模擬典型巖石的無側限壓縮和抗拉數值試驗而言,顆粒間的黏結強度參數(主要是指黏結抗拉強度和黏結抗剪強度)決定數值模擬試樣的無側限壓縮峰值強度和抗拉強度,這就導致在使用PBM來模擬典型巖石時,壓拉比僅能隨著黏結強度參數的變化在3～7范圍內波動[13],不能準確模擬真實典型巖石的壓拉比特性.

對于SBM來說,通過引入軟化系數,使黏結鍵在作用力超過其抗拉強度時并未直接移除,而是進入一個軟化狀態(圖12b).這個軟化狀態使粒間拉力超過黏結抗拉強度時黏結鍵仍有部分承載力,其可顯著增強試樣的抗拉強度和無側限抗壓強度,使較小的黏結強度參數即可模擬獲得較大的無側限抗壓強度和抗拉強度.

圖12 兩種黏結鍵的力學行為[32]

無側限壓縮中接觸黏結抗拉強度和接觸黏結抗剪強度相結合共同支撐試樣,如圖13所示,而抗拉試驗中僅有接觸黏結抗拉強度決定試樣強度.當PBM試樣接觸黏結抗剪強度過高時,在無側限壓縮過程中,試樣的破壞將以拉伸破壞為主,并不能無限提高試樣的抗壓強度,因此,PBM試樣的壓拉比始終維持在一個范圍內.SBM模型則引入軟化系數來增強試樣抗拉和抗壓強度,但抗壓強度增加幅度較抗拉強度大,進而使試樣的壓拉比增大.

圖13 PBM和SBM試樣壓拉比的影響示意圖

由于目前尚未見有關軟化黏結接觸模型SBM的研究,且本文并未對細觀參數間的復雜作用機理進行深入探究,因此,壓拉比的范圍較小.本文后續將結合典型巖石的無側限壓縮和抗拉試驗,深入研究SBM細觀參數對試樣力學性質的影響機制,進一步優化試樣壓拉比.

3.3 討 論

對于PBM模型從細觀角度而言,其黏結鍵間作用力在超過黏結鍵的抗拉強度之后會立即跌落至0,導致黏結失效.對于PBM模型在反映試樣宏觀力學特性中,則是達到峰值強度之后應力的迅速跌落,難以模擬部分軟巖應力-應變關系的峰后延性行為.而對于SBM模型而言,在未激活軟化之前,其表現為PBM的性質,因此可以模擬典型砂巖應力-應變行為;激活軟化后,在黏結鍵間作用力超過其抗拉強度之后,黏結鍵間作用力并未直接跌落至0,而是以一個恒定的速率跌落至某一數值,軟化期間顆粒間應力持續存在,反映到宏觀試樣性質上,在試樣應力達到峰值強度之后不會造成應力的突然跌落,因此,可模擬部分軟巖應力-應變關系的峰后延性行為.

針對李桂臣等[32]研究中所采用的兩個泥巖試樣,采用本文中的基于SBM的剛性塊體模型和基于PBM的球形顆粒模型進行數值模擬,所得應力應變關系如圖14所示.從圖14中可以看出,軟化系數的引入可以有效減少試樣的峰后脆性破壞模式,從而使試樣從峰后脆性向峰后延性轉變.

圖14 泥巖物理試驗和數值模擬應力應變關系

脆性、延性的材料特性最初屬于金屬材料力學范疇,隨后擴展到巖石材料力學中.脆性、延性性質的差異更多地體現在細觀結構構造的差異中,同時脆性、延性的定義和界定也不僅是峰后,發生破壞或者達到峰值強度前的變形量相差顯著[33].

4 結 論

1) 相較于球形顆粒試樣,基于Voronoi剖分生成的剛性塊體試樣顆粒形狀各異,更符合砂巖等含結晶結構巖石中礦物的不規則晶粒形狀.

2) 相較于平行黏結接觸模型僅能模擬出3～7的壓拉比,軟化黏結接觸模型能夠模擬出砂巖更大的壓拉比,更能真實表征砂巖在多應力路徑下的力學特性.

3) 相較于傳統的球形顆粒數值模擬方法,剛性塊體間配位數的增加提高了顆粒間的自鎖效應,使顆粒間的接觸行為更符合真實砂巖晶粒間的力學行為,進而能模擬出砂巖典型的剪切破壞模式.