王桂林，王潤秋, 孫 帆

(1.重慶大學(xué) 土木工程學(xué)院，重慶 400045； 2.庫區(qū)環(huán)境地質(zhì)災(zāi)害防治國家地方聯(lián)合工程研究中心，重慶 400045)

1 研究背景

巖體破壞事實上是由巖體內(nèi)部細觀裂紋的閉合—產(chǎn)生—相互貫通連接而形成的宏觀表現(xiàn)，因此研究微裂紋的行為特點，對于揭示巖體宏觀行為特征是很有益處的[1-3]。大多數(shù)學(xué)者認為，巖體微裂紋的發(fā)展演化多由巖體礦物顆粒及其之間的膠結(jié)所控制[4-6]。自然界中的巖體顆粒多為幾何形狀復(fù)雜的礦物晶粒，并且晶粒與晶粒之間存在相應(yīng)的膠結(jié)物質(zhì)將其連接在一起，當巖體受到外荷載作用后，礦物晶粒便可能發(fā)生晶粒內(nèi)部破壞的穿晶裂紋以及沿晶粒間膠結(jié)破壞的晶間裂紋[7-9]，如圖1所示。

圖1 礦物的穿晶裂紋及晶間裂紋[9]Fig.1 Grain-boundary cracks and transgranularcracks[9]

研究巖體晶粒細觀裂紋的發(fā)展，除了試驗手段，人們越來越傾向于利用更加便利的數(shù)值模擬方法[10-11]。離散元法[12](Discrete Element Method, DEM)便是其中一種在研究巖體細觀行為中較為常用的方法。目前以離散元法為主要計算核心的程序包括PFC2d& PFC3d、UDEC & 3DEC以及開源離散元程序Yade[13]、ESyS-Particle[14]和LIGGGHTS[15]。然而早期離散元模型卻無法考慮礦物結(jié)構(gòu)、晶粒破壞以及巖體中礦物顆粒不同屬性，且模擬結(jié)果與實際情況也存在較明顯的偏差。有的研究表明利用以顆粒流為代表的PFC程序等所得到的單軸抗拉強度與抗壓強度之間的比值相對于實際巖體中的比值偏大[16-17]，Cho等[17]還指出PFC的單軸抗拉強度標定后模擬所得的三軸抗壓強度較低。除此之外，PFC中所標定的某些參數(shù)(摩擦系數(shù)、接觸模量、平行鍵模量等)實際的物理意義不明確。為解決這些問題，Potyondy[18]最早提出PFC-GBM(Grain Based Model)這一概念并用于模擬剛性圓盤顆粒的行為，該方法可以用來模擬不同礦物的行為，但仍然存在以上PFC的缺點，同時對于具有多邊形幾何特征的礦物而言，無法準確描述其幾何特征，因此建立以塊體離散元(UDEC/3DEC)為主要核心的GBM便顯示出其獨特的優(yōu)勢。Voronoi tessellation是一種利用離散點來形成不規(guī)則多邊形或多面體的方法，學(xué)者們常常將其與離散元相結(jié)合來模擬巖體的細觀特征[19-24]。對于UDEC而言，利用Voronoi所生成的模型(GBM)是內(nèi)置于程序中的，所生成的二維多邊形無法考慮到實際礦物生長過程中晶粒取向以及解理等，因此三維塊體離散元3DEC更適合用來利用GBM來解決早期離散元模型無法解決的問題。對此本文將利用Neper軟件生成Voronoi模型，然后導(dǎo)入到3DEC中。

目前對于塊體離散元 GBM的研究相對比較少，多停留在參數(shù)標定階段，Wang等[25-26]以花崗巖為例分析了3種不同顆粒屬性的參數(shù)標定過程，并利用塊體離散元 GBM考慮了晶粒內(nèi)和晶粒間的破壞模型。Ghazvinian等[27]研究了三維狀態(tài)下圓柱形花崗巖的損傷等，同時提出了自己的標定方法。綜合以前的研究，不難發(fā)現(xiàn)，大多文獻所給出的標定程序多圍繞彈性模量與泊松比，對于峰值強度的細觀參數(shù)標定程序不足，同時承載板與試件之間的接觸參數(shù)也少有考慮?；诖?，本文將在前人研究基礎(chǔ)之上，同時考慮承載板接觸參數(shù)以及峰值強度，細化研究對象并提出相應(yīng)的標定程序，之后利用改進的標定方法，標定不同礦物組分的花崗巖，進一步分析花崗巖中的穿晶裂紋以及晶界裂紋，并通過與室內(nèi)試驗結(jié)果對比，驗證改進標定方法后的塊體離散元 GBM方法在巖體細觀演化中的準確性。

2 GBM巖樣宏細觀參數(shù)分析

2.1 塊體離散元 GBM巖樣

模型采用板狀試件，寬高為50 mm×100 mm，厚度為1 mm，晶粒隨機生成，晶粒平均尺寸為3 mm，一共900個晶粒，如圖2(a)所示。模型上下兩側(cè)為加載鋼板，板厚0.05 mm，加載鋼板垂直方向兩側(cè)均施加垂直方向的速度約束，模型XZ平面均進行Y方向的位移約束，使之成為平面應(yīng)變問題。單軸壓縮試驗中控制加載速度為0.1 m/s，直拉試驗中控制加載的速度為0.01 m/s。值得注意的是無論是單軸壓縮試驗還是直拉試驗?zāi)Ｐ途捎盟欧刂?，同時在模型內(nèi)部設(shè)較大局部阻尼值(0.90)以保證在整個試驗的過程中巖體均為準靜態(tài)過程，從而實現(xiàn)準靜態(tài)加載。加載過程中監(jiān)測其橫向應(yīng)變、軸向應(yīng)變以及軸向應(yīng)力。已有文獻研究表明當應(yīng)力監(jiān)測區(qū)域增大時，監(jiān)測區(qū)域大小的影響將會消失[27]。因此，橫向應(yīng)變的監(jiān)測點位于試樣左右兩側(cè)距試件邊緣0.01 m處，軸向應(yīng)變的監(jiān)測點位于兩側(cè)壓板接觸面0.02 m，軸向應(yīng)力為橫向及軸向應(yīng)變監(jiān)測點所包圍區(qū)域的平均值σzz，如圖2(b)所示。

圖2 塊體離散元GBM和監(jiān)測點位置Fig.2 Block discrete element GBM and the position ofmonitoring points

模型的力學(xué)特征受到材料本身屬性的影響，為減少計算時間，本文將晶粒設(shè)為彈性，同時粒間采用庫倫滑移準則，晶粒與承載板之間采用接觸摩擦滑移，參考已有的研究經(jīng)驗[25-27]，表1、表2及表3分別給出了上述3種屬性參數(shù)值。值得注意的是接觸屬性里殘余抗拉強度以及殘余黏聚力設(shè)置為0[28-29]。

表1 承載板及晶粒參數(shù)Table 1 Bearing plate and grain parameters

表2 承載板與晶粒之間的接觸參數(shù)Table 2 Contact parameters between bearing plateand grains

表3 晶粒之間的初始接觸參數(shù)Table 3 Initial contact parameters between grains

2.2 單軸壓縮下粒間接觸屬性的影響

2.2.1 承載板接觸參數(shù)的影響

承載板的接觸參數(shù)主要包括承載板與晶粒之間的切向接觸剛度ks_bj、法向接觸剛度kn_bj以及摩擦角φbj。圖3為單因素下不同接觸剛度及摩擦角的數(shù)值模擬試驗結(jié)果，其中圖3中的(a)、(c)、(e)為承載板接觸參數(shù)對試樣的彈性模量E及泊松比ν的影響，圖3中的(b)、(d)、(f)為承載板接觸參數(shù)對試樣的峰值強度σf的影響。由圖3中的(a)、(c)、(e)可知：改變ks_bj時，試樣的彈性模量基本不發(fā)生變化，泊松比將會隨其增大而逐漸減小，但總體減小的量并不大；改變kn_bj對彈性模量及泊松比基本沒有影響；改變摩擦角，彈性模量基本不發(fā)生變化，但泊松比先減小后增大，并在60°左右達到最小值。由圖3中的(b)、(d)、(f)可知，總體來說試樣的峰值強度受承載板與晶粒間接觸屬性的影響相對較小，偶有起伏，但最大值與最小值之間的差距不大。綜合試驗結(jié)果可以發(fā)現(xiàn)在本試驗條件下，承載板接觸參數(shù)對于試驗結(jié)果的干擾相對較小，因此保證了數(shù)值試驗所采取的監(jiān)測手段以及所記錄的數(shù)據(jù)結(jié)果的可靠性。

圖3 承載板與晶粒之間接觸參數(shù)對試樣力學(xué)參數(shù)的影響Fig.3 Influence of contact parameters between bearingplate and grains on mechanical parameters of specimen

2.2.2 粒間接觸剛度的影響

研究表明[25-27]，粒間接觸剛度將會影響試樣的彈性模量E以及泊松比ν，但很少人給出粒間接觸剛度對峰值強度σf的影響，同時大部分的文獻僅僅給出粒間法向接觸剛度kn_jj不變的情況下粒間接觸剛度比值對彈性模量及泊松比的影響，而忽略了粒間剪切剛度ks_jj不變時粒間剛度比的影響。圖4(a)為保持粒間剪切剛度ks_jj為1 000 GPa/m，粒間剛度比對試樣的彈性模量及泊松比的影響，圖4(b)則為其對峰值強度的影響。增大粒間剛度比，試樣的泊松比、彈性模量以及峰值強度均隨之而增大。然而相對于圖5(a)而言，當保持粒間kn_jj為20 000 GPa/m時，彈性模量卻隨粒間剛度比的增大而減小，泊松比則與圖4(a)有相同的趨勢，同時對比圖4(b)與圖5(b)中的峰值強度，圖5(b)中的峰值呈現(xiàn)先增大后減小的趨勢，達到最大值時所對應(yīng)的粒間接觸剛度比值為5左右。保持粒間接觸剛度比值不變，增大粒間法向接觸剛度，彈性模量、泊松比以及峰值強度均會隨之單調(diào)增加，如圖6(a)、圖6(b)所示。

圖4 ks_jj=1 000 GPa/m時粒間剛度比對試樣力學(xué)參數(shù)的影響Fig.4 Influence of ratio of stiffness between grains on themechanical parameters of specimens(ks_jj=1 000 GPa/m)

圖5 kn_jj=20 000 GPa/m時粒間剛度比對試樣力學(xué)參數(shù)的影響Fig.5 Influence of ratio of stiffness between grains on themechanical parameters of specimens(kn_jj=20 000 GPa/m)

圖6 kn_jj/ks_jj=5時粒間剛度比對試樣力學(xué)參數(shù)的影響Fig.6 Influence of ratio of stiffness between grains onthe mechanical parameters of specimens (kn_jj/ks_jj=5)

綜合以上分析，試樣的彈性模量、峰值強度以及泊松比中，僅泊松比隨粒間接觸剛度的變化呈現(xiàn)單調(diào)性，因此在標定粒間接觸剛度時，首先應(yīng)考慮泊松比。通過計算圖4(a)、圖5(a)、圖6(a)中泊松比的平均差可以發(fā)現(xiàn)，當保持粒間接觸剛度比值不變時，其平均差(0.027)相對于保持粒間法向接觸剛度不變時的平均差(0.046)以及粒間剪切剛度不變時的平均差(0.079)最小，其離散性較小，能夠保證在標定過程中泊松比僅在小范圍內(nèi)變化。故在調(diào)整粒間接觸剛度時應(yīng)先觀察初始計算后的彈性模量及泊松比，若兩者均大于或小于試驗值，先保持粒間接觸剪切剛度不變，然后調(diào)節(jié)粒間接觸法向接觸剛度，如圖4(a)所示，反之保持粒間法向接觸剛度不變，調(diào)整粒間剪切剛度，如圖5(a)所示。由此標定泊松比的細觀參數(shù)同時使彈性模量趨近于實際值，之后再保持粒間剛度比，調(diào)整粒間法向接觸剛度，使得模量快速到達實際值。

2.2.3 粒間摩擦角、黏聚力、抗拉強度的影響

圖7、圖8、圖9分別為粒間摩擦角φjj、黏聚力cjj以及抗拉強度σt_jj對試樣的力學(xué)參數(shù)彈性模量E、泊松比ν以及峰值強度σf的影響。

圖7 粒間摩擦角對試樣力學(xué)參數(shù)的影響Fig.7 Influence of contact friction angle betweengrains on mechanical parameters of specimen

圖8 粒間黏聚力對試樣力學(xué)參數(shù)的影響Fig.8 Influence of contact cohesion between grains onmechanical parameters of specimen

圖9 粒間抗拉強度對試樣力學(xué)參數(shù)的影響Fig.9 Influence of contact tensile strength betweengrains on mechanical parameters of specimen

由圖7(a)、圖8(a)、圖9(a)可知，粒間摩擦角、黏聚力以及抗拉強度對試樣的彈性模量、泊松比影響相對較小，彈性模量基本維持在31.3 GPa左右，而泊松比則基本維持在0.25。由圖7(b)、圖8(b)、圖9(b)可知，粒間摩擦角、黏聚力以及抗拉強度對試樣的峰值強度的影響中，粒間摩擦角的影響最為顯著，如圖7(b)所示；粒間黏聚力對峰值強度的影響次之，趨勢與粒間摩擦角的影響一樣，如圖8(b)所示；粒間抗拉強度對峰值強度的影響最小，隨著粒間抗拉強度的增加，試樣的峰值強度增加緩慢，如圖9(b)所示。根據(jù)以上分析結(jié)果，適宜用粒間摩擦角、黏聚力及抗拉強度來標定峰值強度，其平均差分別為107.88、20.09、3.33 MPa，因此數(shù)值試驗所得的峰值強度過高或者過低時，可采用略微減小或增大粒間摩擦角的做法，因為這樣可以快速改變峰值強度，之后再改變粒間黏聚力以及粒間抗拉強度使之趨于實際值。值得注意的是在調(diào)整粒間摩擦角時應(yīng)選擇較小值，因為當粒間摩擦角的值增大到一定程度時(>60°)泊松比會產(chǎn)生突變，如圖7(a)所示。

2.3 抗拉強度細觀參數(shù)的標定

本文利用直拉試驗進行抗拉強度細觀參數(shù)的標定，所采用的模型和壓縮試驗保持一致，為保證拉伸時試樣與承載板之間不發(fā)生分離，初始時承載板和晶粒之間的接觸屬性改用摩爾庫倫滑移，并將黏聚力設(shè)置為較大值。試驗結(jié)果表明試樣的抗拉強度(σt)主要受晶粒間黏聚力Cjj及粒間抗拉強度σt_jj影響，圖10給出了這種影響對試樣整體抗拉強度的作用效果。由圖10(a)可以看出，晶粒間黏聚力對試樣抗拉強度的影響較弱，抗拉強度基本不隨著粒間黏聚力的大小發(fā)生太大的變化，極差僅為0.328 MPa。試樣的抗拉強度對于粒間抗拉強度的變化最為敏感，如圖10(b)所示，當粒間的抗拉強度增大時，試樣的抗拉強度也急劇增加，但始終小于粒間抗拉強度，其值大小約為粒間抗拉強度的一半。因此在校準抗拉強度時，改變粒間抗拉強度是最合理的。

圖10 粒間黏聚力和粒間抗拉強度對試樣抗拉強度的影響Fig.10 Influence of contact cohesion and contacttensile strength between grains on tensile strengthof specimen

3 單軸拉壓下塊體離散元GBM參數(shù)標定方法

綜合以上單軸壓縮試驗以及直拉試驗，以實驗室所得的巖體彈性模量、泊松比、峰值強度以及抗拉強度為參照對象，對塊體離散元 GBM數(shù)值模擬方法中的細觀參數(shù)進行標定，根據(jù)上述獲得的影響規(guī)律，其標定方法如圖11所示。

圖11 塊體離散元GBM標定方法Fig.11 Discrete element GBM calibration method

標定方法說明：模型標定時首先應(yīng)考慮初始計算結(jié)果中的泊松比和彈性模量與實際泊松比和彈性模量之間的大小關(guān)系，改變粒間接觸剛度以校準泊松比同時使得彈性模量快速趨于實際值，減少標定時間。然后，保持泊松比不變，調(diào)節(jié)法向接觸剛度以使彈性模量符合實際彈性模量，之后保持接觸剛度不變，按照此時的峰值強度大小調(diào)節(jié)粒間黏聚力以及粒間摩擦角以調(diào)整峰值強度達到實際值，最后通過直拉試驗調(diào)節(jié)粒間接觸抗拉強度以匹配總體的抗拉強度。

4 算例驗證

4.1 花崗巖塊體離散元GBM

花崗巖主要由長石、石英、黑云母等礦物組成。在建立數(shù)值模型時，目前已有的數(shù)值方法多將其考慮為同一種均質(zhì)的礦物，并將不同參數(shù)進行平均處理；部分學(xué)者利用其他GBM方法盡管考慮到了礦物的差異性，但忽略了實際礦物生長過程中解理的影響。因此為了使本文采用的塊體離散元GBM更加貼合實際，模型中的礦物作以下處理：

(1)長石具有2組相互垂直或接近垂直的相交解理，發(fā)生破壞時其裂紋一般沿著解理[30]，為簡化計算，數(shù)值采用平面應(yīng)變模型，因此僅考慮其中的一組解理。

(2)石英完全無解理，受力破壞時毫無規(guī)則，如果沒有預(yù)先存在裂紋則很難破壞[30]，故在計算之前先將其礦物晶粒預(yù)制為6個不均等的微晶粒。

(3)黑云母完全解理，但是當晶間或穿晶裂紋發(fā)展到黑云母時，其裂紋擴展大多沿云母片之間的夾層，形成穿晶裂紋的情況較為少見[31]，在平面應(yīng)變模型下考慮為片狀黑云母，不發(fā)生破壞。

基于此，為了實現(xiàn)晶粒本身的破壞，本文根據(jù)每一種礦物的解理性質(zhì)對每一個晶粒進行再次切割，預(yù)制成含解理的礦物，如圖12所示。本文所采用的花崗巖成分如表4所示[32]，由于綠泥石和方解石的含量較低，對整個模型的影響較小，因此為簡化模型未考慮這2種礦物，并將其占比分配給長石，數(shù)值模型采用調(diào)整后的花崗巖成分，如表5所示。

圖12 花崗巖數(shù)值模型Fig.12 Numerical model of granite

表4 花崗巖礦物組成比例[32]Table 4 Components and individual proportionof granite[32]

表5 模型中的花崗巖礦物組成Table 5 Modified components and individual proportionof granite

根據(jù)以上結(jié)果建立的塊體離散元GBM如圖12所示。利用上述標定方法對模型的參數(shù)進行標定，試驗彈性模量為31 GPa[32]，標定后的彈性模量為29.8 GPa，兩者差距不大，而最終獲得的接觸參數(shù)如表6所示。值得注意的是若要得到晶粒內(nèi)的接觸參數(shù)和晶界中的接觸參數(shù)是很復(fù)雜的，因此為減少計算量，數(shù)值模型晶粒間解理的參數(shù)和晶界間的參數(shù)采用同一組參數(shù)。

表6 晶粒及粒間接觸的細觀參數(shù)Table 6 Mesoscopic parameters of grain and contact

4.2 花崗巖晶界裂紋及穿晶裂紋分析驗證

利用自編的fish程序，統(tǒng)計單軸壓縮下花崗巖晶界裂紋(晶體與晶體之間的已破壞的接觸)以及晶內(nèi)裂紋(晶粒內(nèi)部破壞的接觸)累計統(tǒng)計曲線如圖13所示。

圖13 花崗巖穿晶裂紋及晶界裂紋統(tǒng)計Fig.13 Statistics of granite’s transgranular cracksand grain boundary cracks

由圖13可知，在早期的加載階段晶界裂紋與穿晶裂紋均未出現(xiàn)，隨著加載的進行晶界裂紋與穿晶裂紋的數(shù)量均開始增加，并且其增量逐漸增大，但是晶界裂紋的數(shù)量始終多于穿晶裂紋的數(shù)量，這與室內(nèi)試驗的結(jié)論相似[33]。

4.3 不同礦物之間晶界裂紋分析驗證

Ghasemi等[34]經(jīng)過大量的細觀試驗指出，長石的彈性模量小于石英的彈性模量，因此石英-黑云母晶界裂紋數(shù)目大于長石-黑云母晶界裂紋數(shù)目，如圖14所示。然而針對類似的這種現(xiàn)象，大多仍停留在細觀試驗上，試驗過程復(fù)雜且難以掌控。為解決這一難題，室內(nèi)數(shù)值模擬便成了一種良好的手段，但是目前已有的數(shù)值方法少有考慮晶界裂紋，更不用說不同種礦物之間的晶界裂紋。基于此，本文利用GBM模擬不同礦物晶界面細觀裂紋開展行為。由于試樣中礦物的含量不一樣，為保證試驗結(jié)果的準確性，其他條件保持不變，將長石和石英的含量調(diào)換，并在同樣的情景下加載，統(tǒng)計結(jié)果僅考慮其裂紋的數(shù)目，忽略峰值以及彈性模量的影響，得到如圖15所示的石英-黑云母與長石-黑云母累計晶界裂紋之差的統(tǒng)計曲線。

圖14 石英-黑云母晶界裂紋及長石-黑云母晶界裂紋[34]Fig.14 Quartz-biotite grain-boundary cracks andfeldspar-biotite grain-boundary cracks[34]

圖15 石英-黑云母與長石-黑云母之間的差值Fig.15 Difference between quartz-biotite andfeldspar-biotite

同樣的加載條件下，花崗巖單軸壓縮下礦物晶界的破碎和礦物的彈性模量息息相關(guān)，如圖15所示：彈性階段，兩者之間的差值逐漸增加，但增加的程度較小，此階段石英的彈性模量較大，抵抗變形的能力較大，黑云母的彈性模量小，抵抗變形的能力較小，兩者在晶界處容易產(chǎn)生晶界裂紋，而長石的彈性模量較小，可以與黑云母一起變形，故在晶界處產(chǎn)生的晶界裂紋較少，但總的來說兩者之間的差值不大；彈塑性階段，兩者之間的差值出現(xiàn)下降，長石-黑云母之間晶界裂紋的增長量大于石英-黑云母長石之間的晶界裂紋，其原因在于當長石的含量占花崗巖主要成分時，花崗巖的總體彈性模量小于石英占主要成分時，花崗巖容易進入塑性，裂紋增加明顯，因此出現(xiàn)下降段；峰后花崗巖破壞階段，石英占主要成分時，最先達到峰值強度，石英-黑云母由于能量的突然釋放晶界裂紋呈現(xiàn)爆發(fā)式的增長，其值遠大于長石-黑云母之間的晶界裂紋，此時長石占主要成分的花崗巖還未達到峰值，兩者之間的差值增量逐漸增大，直到后者達到峰值，差值開始穩(wěn)定增長，達到差值峰值之后兩者之間的裂紋差逐漸減小?？偟膩碚f兩者之間的差值在整個加載計算過程中均>0，因此可以證明花崗巖單軸壓縮下石英-黑云母之間的晶界裂紋大于長石-黑云母之間的晶界裂紋，這與實際的室內(nèi)細觀試驗結(jié)果一致[34]。

5 結(jié) 論

(1)巖體單軸壓縮試驗以及直拉試驗3DEC-GBM數(shù)值模擬結(jié)果表明，巖體彈性模量以及泊松比的細觀參數(shù)標定主要取決于晶粒間的接觸剛度，巖體峰值強度主要與粒間黏聚力、抗拉強度以及摩擦角的有關(guān)，而巖體的抗拉強度受晶粒間抗拉強度的影響最為顯著。

(2)花崗巖單軸壓縮過程中在礦物內(nèi)部以及晶界處均會產(chǎn)生細觀裂紋。經(jīng)與試驗結(jié)果對比，呈現(xiàn)的規(guī)律基本一致，花崗巖內(nèi)晶界裂紋的數(shù)量明顯大于穿晶裂紋的數(shù)量。同時其晶界裂紋受礦物的彈性模量影響較大，長石的彈性模量小于石英的彈性模量，因此在加載的過程，其與黑云母之間的晶界裂紋數(shù)小于石英與黑云母之間的晶界裂紋數(shù)，并且兩者之間的差值大小與加載過程中花崗巖內(nèi)部的應(yīng)力狀態(tài)有關(guān)。

(3)由于標定過程中未考慮不同礦物之間的接觸以及礦物內(nèi)部的接觸，也未考慮圓柱體等試件，因此，推廣本文方法時需要對這些因素進一步研究。