趙寧寧,李二強,馮吉利
(1.中國礦業(yè)大學(xué)(北京) 深部巖土力學(xué)與地下工程國家重點實驗室,北京,100083;2.中國礦業(yè)大學(xué)(北京) 力學(xué)與建筑工程學(xué)院,北京,100083;3.洛陽理工學(xué)院 土木工程學(xué)院,河南 洛陽,471000)
頁巖、片麻巖等沉積巖及板巖、煤巖等變質(zhì)巖中通常存在著發(fā)育良好的層理面,這些規(guī)模不等的層理將完整巖體切割成橫觀各向同性巖體,由此引起的巖石力學(xué)特性差異對工程巖體的安全與穩(wěn)定性有重要影響[1-4]。由于層狀巖體具有明顯的橫觀各向同性,巖體強度不僅僅與自身強度有關(guān),結(jié)構(gòu)面的分布及性質(zhì)也會對巖體造成一定的影響,且破壞劈裂形式往往與結(jié)構(gòu)面的產(chǎn)狀相關(guān)。為了研究層狀巖體特性及破裂形態(tài),人們采用離散單元法(discrete element model,簡稱DEM)將剛性圓盤和球體作為離散顆粒,通過相互作用來模擬組件的力學(xué)行為。隨著離散元的發(fā)展及計算機性能的逐漸提升,計算效率顯得尤為重要。已有的離散元軟件如顆粒流程序在二維計算上效率較高,對于三維大規(guī)模模型,計算效率相對較低?;陔x散元框架MUSEN,采用GPU加速,可以實現(xiàn)大規(guī)模的計算。在MUSEN中,通過固體鍵將顆粒連接起來以模擬層狀巖體在準(zhǔn)靜態(tài)壓縮過程中的破壞行為,固體鍵充當(dāng)顆粒間黏聚力的作用。因此,開展層狀巖石抗壓力學(xué)特性試驗及基于離散元框架MUSEN提高數(shù)值模擬效率對指導(dǎo)大規(guī)模工程實踐具有重要的意義。
學(xué)者們針對BPM理論及模擬進(jìn)行了大量研究?;趶椥约僭O(shè),在均質(zhì)球形顆粒中,許多模型可以描述其應(yīng)力場。HERTZ[5]計算了接觸區(qū)域的壓力分布、接觸力以及接觸面的半徑。HUBER[6]提出的本構(gòu)模型不僅可以描述半空間顆粒接觸區(qū)域的應(yīng)力分布,而且可以應(yīng)用于計算球體的接觸體積[7]。HAHN[8]計算了整個球體的空間應(yīng)力分布。ANTONYUK 等[9]概述了求解顆粒間接觸力-位移的本構(gòu)方程。顆粒的力學(xué)破壞行為由微觀黏接機制決定,已有的研究證實了顆粒鍵的微觀及宏觀相關(guān)性[10-12]。RUMPF 模型[13]描述了顆粒的抗拉強度,該模型考慮了隨機模型中單個球體顆粒在接觸點處的黏性力。BIKA 等[14]推廣了RUMPF 模型[13]和KENDALL[12]模型,并將其推廣到多孔顆粒結(jié)構(gòu)中,發(fā)現(xiàn)具有高黏性液體黏合劑和高飽和度的顆粒,其抗拉強度依賴于毛細(xì)管壓力[11]。DELENNE 等[15]采用莫爾-庫侖破壞準(zhǔn)則對耦合圓柱桿在不同的應(yīng)力或加載模式(拉伸、壓縮和剪切)下的變形和破壞行為進(jìn)行了建模和模擬。在數(shù)值模擬方面,主要采用DEM 對顆粒破裂進(jìn)行研究。基于離散元模擬計算單個顆粒在離散時間步中的運動。最初,離散單元法是為模擬理想的球形顆粒群而開發(fā)的,然而,該方法已被推廣應(yīng)用于模擬非球形物體或原生顆粒組成的團聚體。為了重現(xiàn)團聚體的形狀和內(nèi)部結(jié)構(gòu),可以使用鍵合顆粒模型(BPM)方法[16-17]。BPM 將顆粒表示為一組由理想彈性或黏彈性固體鍵連接的初始顆粒[16,18-23],可以在顆粒級尺度上獲得作用在顆粒上的外力和顆粒內(nèi)部應(yīng)力的詳細(xì)數(shù)據(jù),而這些數(shù)據(jù)不能從實驗測量中獲得。與研究彈性顆粒復(fù)合材料力學(xué)性能的BPM 方法類似,另一種方法則將材料視為一組彈性彈簧。然而,與BPM 方法相反,該方法不能有效地描述材料的斷裂。BPM 方法還可以用于分析在實驗中很難分離的結(jié)構(gòu)參數(shù)的影響。例如,在BPM 中可以很容易地隔離實心橋強度的影響,而在實驗中使用不同的黏結(jié)劑可能改變顆粒的內(nèi)部結(jié)構(gòu),致使實心橋強度對團聚體強度的影響變得復(fù)雜。在試驗方面,眾多學(xué)者對層狀巖體的抗拉、抗壓力學(xué)特性開展了相關(guān)研究。MOUSAVI NEZHAD 等[24-26]通過對層狀頁巖和片巖等開展間接拉伸試驗、有限元及離散元數(shù)值分析,總結(jié)了層理及基質(zhì)對抗拉及抗壓強度的影響,發(fā)現(xiàn)層理對抗拉強度和劈裂破壞具有顯著影響。針對層狀板巖,劉運思等[27-30]通過單軸及巴西劈裂試驗揭示了橫觀各向同性強度特征和破壞機制;馬天壽等[31]通過巴西劈裂以及聲發(fā)射檢測試驗,明確頁巖橫向與縱向彈性模量、抗拉強度及破壞模式各向異性特征。
BPM是DEM的擴展,已得到廣泛應(yīng)用。在離散元法中,每一對顆??梢酝ㄟ^固體鍵連接。POTYONDY 等[21,23]通過定義化學(xué)鍵的本構(gòu)關(guān)系,可以捕獲不同的宏觀效應(yīng),如材料硬化或軟化。近年來,BPM 已被越來越多用于巖石材料的二維和三維模擬數(shù)值分析,但較少應(yīng)用于層狀炭質(zhì)板巖中。當(dāng)前在巖石層理的數(shù)值模擬研究中存在以下幾個問題:1) 模擬層理巖石破壞時多為二維模擬,無法展現(xiàn)其內(nèi)部結(jié)構(gòu)破裂特性,而三維離散元模擬軟件如離散元軟件PFC 對于大規(guī)模模型的計算效率較低。2) 構(gòu)建層理巖體模型時,其層理結(jié)構(gòu)大多數(shù)給定很薄的厚度,為一條直線(二維)或平滑平面(三維),而真實層理結(jié)構(gòu)面應(yīng)為鋸齒狀結(jié)構(gòu)。因此,需要采用合適的層理模型及數(shù)值方法進(jìn)行模擬。本文作者首先構(gòu)建不同層理傾角的層理模型;其次,研究微觀與宏觀彈性模量的關(guān)系并探究法向抗拉與切向抗剪強度對單軸抗壓強度的影響,分析不同模型對層狀板巖力學(xué)響應(yīng)的影響。最后,分別對傾角為0°,30°,45°,60°和90°層理進(jìn)行數(shù)值模擬并與試驗結(jié)果進(jìn)行對比,并對層狀板巖靜力破壞的主要裂紋萌生以及力學(xué)響應(yīng)進(jìn)行分析。
試驗樣品取自陜西木寨嶺隧道施工隧道掌子面,取樣點巖體呈薄層狀構(gòu)造,所取巖塊賦存層理顯著,其層理存在規(guī)模不等的差異,經(jīng)粗略測量層間距為4.5~22.5 mm。該試樣密度為2.688 g/cm3,經(jīng)X 射線衍射測定試樣主要由石英與黏土礦物組成。層狀炭質(zhì)板巖基本力學(xué)參數(shù)測定如下:水平層理炭質(zhì)板巖單軸抗壓強度約為49.2 MPa,彈性模量E約為7.0 GPa,泊松比μ為0.20;而垂直層理炭質(zhì)板巖單軸抗壓強度約為52.5 MPa,彈性模量E′約為7.8 GPa,泊松比μ′為0.23。
圖1所示為層狀炭質(zhì)板巖單軸試樣示意圖及試樣實物圖,底面直徑為50 mm、高為100 mm。試樣加工工藝主要包括鉆取巖芯、切割柱型圓盤及平整試樣。為減小水分對巖體特性的影響,將制備的試樣在通風(fēng)條件下靜置60 d。通過烘干法計算可得:試樣制備完成時,其含水率為1.63%~1.95%,通風(fēng)靜置后,含水率為0.21%~0.24%。假設(shè)層理面與水平面夾角為層理傾角,本試驗中,取0°,30°,45°,60°及90°這5種傾角的層狀炭質(zhì)板巖試樣。
圖1 不同層理傾角炭質(zhì)板巖單軸壓縮試樣Fig.1 Uniaxial compression samples of slate under different bedding angles
此次試驗采用200 t 大噸位巖石加載系統(tǒng),如圖2所示,該系統(tǒng)由主機、加載系統(tǒng)、測量控制系統(tǒng)、計算機控制和數(shù)據(jù)處理系統(tǒng)組成,試驗加載時采用位移控制,加載速率為0.1 mm/min。
圖2 炭質(zhì)板巖單軸試驗系統(tǒng)Fig.2 Uniaxial compressive setup for carbonaceous slate
在BPM中,每1對初始顆粒可以通過1個或多個固體鍵連接。每個顆粒和每個固體鍵都可以設(shè)定獨特的幾何和材料參數(shù)以構(gòu)建具有復(fù)雜結(jié)構(gòu)的非均質(zhì)材料。本文采用的是軟球離散單元,顆粒間或顆粒與邊界重疊可以看作是局部材料變形。本文研究3種不同類型的交互作用:1) 初始顆粒之間相互作用即固體鍵連接。2) 顆粒與邊界幾何之間的相互作用。3) 固體鍵用來描述相鄰顆粒間相互作用的力、力矩和扭矩。由于炭質(zhì)板巖具有近乎理想的彈性特性,顆粒-顆粒接觸采用線彈性接觸模型,顆粒-邊界面接觸采用Hertz-Mindlin 模型[22]。對于顆粒間的相互作用,分別在法向和切向?qū)αM(jìn)行分解。
式中:Fn,t為法線方向上的合力;Fn,d為法線方向上的阻尼力;kn為顆粒法向抗拉強度;e為2 個顆粒間恢復(fù)系數(shù);rn為單位向量,且方向與法向力Fn方向一致;Vn為兩顆粒重疊體積;υrel,n為2個顆粒間相對速度;M*為質(zhì)量分別為m1和m2的2個顆粒的等效質(zhì)量。
式中:Ft,t和Ft,d分別為切線方向上的合力和阻尼力;kt為顆粒切向抗剪強度;Ft,p為前一次迭代的切向力矢量;ut為在當(dāng)前時間步長切向重疊的增量;υrel,t為相互作用顆粒的相對速度。此外,F(xiàn)t受到與滑動摩擦有關(guān)的顆粒間摩擦因數(shù)μsl的限制:
若不滿足式(8),則根據(jù)式(9)修正式(6)中Ft:
顆粒與顆粒間合力為法向力矢量與切向力矢量之和:
初始顆粒之間的連接鍵為具有初始長度和半徑的理想圓柱形實體固體鍵,半徑不能超過接觸對最小半徑,且固體鍵在法向和切向附加了阻尼力[16,23]。使用固體鍵模型的問題之一是在模擬過程中幾乎無阻尼振蕩。為了減少這種影響,增加人工阻尼器。一般來說,作用在固體鍵上的應(yīng)力是相應(yīng)顆粒之間相互作用的結(jié)果。為了模擬材料的斷裂過程,分析各黏結(jié)劑中的應(yīng)力,并與材料的抗拉/抗壓強度σmax和切向強度τmax進(jìn)行比較。若滿足式(11)或式(12)中的1 個條件,則固體鍵斷開并從模型中移除。
式中:Ft,b和Fn,b分別為固體鍵兩端的軸力和剪力;固體鍵矢量合力可表示為F2=Ft,b+Fn,b(其中Ft,b和Fn,b分別為固體鍵法向力與切向力矢量),顆粒-顆粒間矢量合力與固體鍵矢量合力之和為顆粒的總力,即F=F1+F2,且根據(jù)牛頓第二運動定理即可求出顆粒的加速度;Mt,b和Mn,b分別為作用在鍵上的扭轉(zhuǎn)力矩和彎矩;Ab和Rb分別為橫切面面積和半徑;Jb和IT,b分別為扭轉(zhuǎn)慣量和轉(zhuǎn)動慣量。
圖3 所示分別為0°,30°,45°,60°和90°的層理模型。圖4 所示為30°層理模型的邊界條件。其中,圖4(a)所示為顆粒與固體鍵模型,圖4(b)所示為固體鍵(綠色代表層理固體鍵,紅色代表基質(zhì)固體鍵),圖4(c)所示為顆粒模型。模型中顆粒直徑為1 mm,孔隙率為0.38,顆粒數(shù)為300 975個,固體鍵單元數(shù)為3 141 792 個。構(gòu)建模型方法如下:1) 首先在MUSEN中形成底面直徑為50 mm、高為100 mm的幾何形狀。2) 根據(jù)直徑及孔隙率在幾何內(nèi)部形成相應(yīng)的顆粒及固體鍵,在形成顆粒過程中可以使用GPU 加速。3) 將形成的顆粒及固體鍵以文本格式提取出來,通過程序編制,將等間距的層理固體鍵單元挑選出來,最終形成如圖3所示的層理。該層理模型構(gòu)建方式有以下2個優(yōu)點:首先,使用同一種模型模擬5種層理,其內(nèi)部構(gòu)造相同。其次,從圖4(b)可以看出,層理固體鍵在模型中呈鋸齒狀,與巖體材料構(gòu)造較接近。圖4所示為層理30°時的邊界條件,上下為加載板,為了與試驗條件保持一致,將下板固定,上板以0.1 mm/min的速度向下移動。
圖3 不同傾角層狀板巖物理模型Fig.3 Physical model of layered slate with different inclination angle
圖4 層理傾角為30°時的邊界條件Fig.4 Boundary conditions when the bedding angle is 30°
3.1.1 微觀與宏觀彈性模量校正
以均質(zhì)板巖顆粒模型為標(biāo)本,微觀彈性模量與宏觀彈性模量關(guān)系曲線如圖5所示。其中,微觀彈性模量代表數(shù)值模擬時顆粒的彈性模量,宏觀彈性模量代表根據(jù)荷載-位移曲線求出斜率并進(jìn)行換算得出的試驗彈性模量。數(shù)值模擬中,微觀彈性模量的選取至關(guān)重要。為了使試件宏觀彈性模量與微觀彈性模量對應(yīng),開展模擬試驗對其進(jìn)行校正。具體方法為:在幾何尺寸相同的條件下,隨機構(gòu)建模型1~6共6種顆粒模型,保證每組模型顆粒直徑及孔隙率一致,并對每組模型取不同的微觀彈性模量,分別為2,8,12和16 GPa,之后,可以根據(jù)每組顆粒模型(對應(yīng)4個彈性模量)的應(yīng)力-應(yīng)變曲線求出宏觀彈性模量,其中模擬邊界條件與圖4 中的一致,下部鋼板保持不變,頂板以0.1 mm/min的速度向下移動。
通過圖5可以看出,微觀模量與宏觀彈性模量大致呈線性關(guān)系,因此,在測得試驗彈性模量的情況下,可根據(jù)插值計算求取微觀彈性模量。對于此次單軸壓縮試驗彈性模量約為7.0 GPa,通過參數(shù)校正可得微觀彈性模量為6.35 GPa。
圖5 微觀模量與宏觀模量關(guān)系曲線Fig.5 Relationship between the micro-modulus and the macro-modulus
3.1.2 固體鍵拉剪強度對巖體力學(xué)特性的影響及加載速度的選取
下面分析固體鍵切向抗剪強度與法向抗拉強度對巖體力學(xué)特性的影響。在相同顆粒模型條件下(如表1 所示),剪切強度約比抗拉強度大2 倍。當(dāng)?shù)缺仍黾訌姸葧r(其中模擬邊界條件與圖4 中的一致),下部鋼板保持不變,頂板以0.1 mm/min 的速度向下移動,最終的荷載-位移曲線如圖6所示,其中τ表示切向剪切強度。從圖6 可以看出,隨著固體鍵拉、切強度逐漸增加,對應(yīng)的峰值荷載及位移均增加,說明隨著拉、切強度逐漸增加,顆粒間的黏聚力增強,致使峰值荷載增加;另一方面,需選取較小的加載使模型處于準(zhǔn)靜態(tài)加載,以消除動能對結(jié)果的影響。當(dāng)加載速度過大時,荷載-位移曲線會出現(xiàn)間斷性的波浪式上升,而隨著速度逐漸減小,逐漸呈現(xiàn)1條光滑的曲線。經(jīng)過測試,當(dāng)速度取0.1 mm/min 時,模型大致處于準(zhǔn)靜態(tài)加載狀態(tài)。固體鍵的拉、切強度對荷載峰值影響較大,因此,基于上述原則且通過試錯法選取合適的法向抗拉和切向抗剪強度,以增強數(shù)值模擬結(jié)果與試驗結(jié)果的吻合度。
表1 基質(zhì)與層理固體鍵單元參數(shù)Table 1 Matrix and layered solid bond element parameters
圖6 同一種模型不同切向強度(法向強度)下荷載-位移曲線Fig.6 Load-displacement curves of different tangential stiffness(normal stiffness) under the same model
在數(shù)值模擬中,邊界條件與圖4中的一致,下部鋼板保持不變,頂板以0.1 mm/min 的速度向下移動,根據(jù)隨機性原理,形成5種不同顆粒分布模型(分別編號為模型Ⅰ~Ⅴ),分析不同的顆粒模型對巖體力學(xué)響應(yīng)的影響。圖7所示為單軸壓縮下的均質(zhì)模型破壞過程示意圖。5種不同顆粒模型條件下對應(yīng)的荷載-位移曲線如圖8 所示。從圖8 可以看出:盡管顆粒模型不同,但是最終力學(xué)響應(yīng)曲線大致相近。因此,為了減小模型的影響,使用同一種顆粒模型進(jìn)行數(shù)值計算分析。
圖7 單軸壓縮下模型劈裂破壞過程示意圖Fig.7 Schematic diagrams of model splitting failure process under uniaxial compression
圖8 相同參數(shù)及邊界條件下不同顆粒模型荷載-位移曲線Fig.8 Load-displacement curves of different physical models under the same parameters and boundary conditions
對于此次單軸壓縮試驗,樣品取自陜西木寨嶺隧道掌子面,在層理模型切向強度與法向強度的選取上,對于巖石材料,其抗剪強度大于抗拉強度,因此,切向強度取值大于法向強度,且基質(zhì)的黏結(jié)性比層理黏結(jié)性高,即層理固體鍵強度大于基質(zhì)固體鍵強度?;谝陨显瓌t且通過試錯法選取法向和切向強度,以增強數(shù)值模擬結(jié)果與試驗結(jié)果的吻合度(見表1)。
圖9 所示為0°,30°,45°,60°及90°層理模型下試驗與模擬荷載-位移曲線對比結(jié)果。從圖9 可以看出:1) 0°,30°,45°,60°及90°層理模型的峰值荷載依次為98.374,81.457,78.673,80.751 和102.595 kN,隨著傾角逐漸增加,峰值荷載大致呈現(xiàn)“U”型分布。2) 模擬曲線與試驗曲線均呈現(xiàn)出線彈性、瞬時劈裂跌落的發(fā)展態(tài)勢,試樣呈現(xiàn)為顯著脆性劈裂破壞特征。3) 模擬與試驗荷載-位移曲線大致吻合。
圖9 不同層理下試驗與模擬荷載-位移曲線Fig.9 Load-displacement curves of test and simulation under different beddings
圖10 所示為0°,30°,45°,60°及90°層理條件下,試樣與模型破壞過程示意圖。由圖10可見:
1) 在0°層理破壞過程的起始階段,層理固體鍵斷裂,且垂直于層理的基質(zhì)固體鍵也逐漸斷裂;隨著頂板壓力的積累,基質(zhì)固體鍵斷裂,并伴隨少量的層理固體鍵斷裂,最終發(fā)生劈裂拉伸破壞,見圖10(a)。
2) 在30°層理破壞過程的起始階段,兩端層理結(jié)構(gòu)面發(fā)生剪切滑移,隨著上頂板壓力逐漸增加,中部結(jié)構(gòu)面逐漸被破壞,且沿層理結(jié)構(gòu)面法向發(fā)生劈裂破壞,最終呈現(xiàn)基質(zhì)結(jié)構(gòu)面與層理結(jié)構(gòu)面交錯劈裂破壞,見圖10(b)。
3) 在45°與60°層理破壞過程的起始階段,裂紋隨著上端或下端發(fā)生層理結(jié)構(gòu)面滑移剪切破壞,隨著荷載逐漸增加,剪力逐漸增加,中間層理結(jié)構(gòu)面逐漸被破壞,最終沿著層理結(jié)構(gòu)面形成貫通的裂紋,破壞模式為剪切滑移破壞,見圖10(c)和(d)。
4) 在90°層理破壞過程的起始階段,巖體兩端基質(zhì)結(jié)構(gòu)面與層理結(jié)構(gòu)面幾乎同時發(fā)生斷裂,隨著時間延長,裂紋沿著層理結(jié)構(gòu)發(fā)生劈裂拉伸破壞,最終形成1條或多條貫通裂紋,見圖10(e)。
圖10 不同層理條件下試樣與模型破壞過程對比Fig.10 Comparison of model and samples failure process under different bedding conditions
5) 模型最終破壞形態(tài)與試驗結(jié)果較接近,且數(shù)值模型能夠更好地再現(xiàn)劈裂破壞過程。
圖11 所示為0°層理模型內(nèi)部固體鍵單元破壞過程。從圖11 可以看出:當(dāng)固體鍵顏色變?yōu)榧t色時,根據(jù)式(11)和(12),此位置的鍵將會從模型中移除。從圖11(a)~(c)可以發(fā)現(xiàn):層理固體鍵先斷裂,隨著時間的延長,沿著層理法向方向的固體鍵力逐漸增加,最終形成1條貫穿裂紋,模型發(fā)生劈裂破壞,如圖11(d)~(f)所示。
圖11 0°層理模型內(nèi)部固體鍵單元破壞過程Fig.11 Failure processes of the internal solid bond in the 0° layering model
圖12所示為30°層理模型內(nèi)部固體鍵單元破壞過程。從圖12 可以看出:在初始階段,模型沿著棱角層理破裂;隨著時間延長,向內(nèi)部層理發(fā)生滑移斷裂,且沿著層理法向方向,固體鍵發(fā)生部分?jǐn)嗔眩罱K模型以滑移劈裂模態(tài)發(fā)生失穩(wěn)破壞。
圖12 30°層理模型內(nèi)部固體鍵單元破壞過程Fig.12 Failure processes of the internal solid bond in the 30° layering model
1) 通過對不同宏觀彈性模量及微觀彈性模量進(jìn)行校正,得出其呈線性關(guān)系,通過插值法可以將試驗彈性模量轉(zhuǎn)換為微觀彈性模量進(jìn)行數(shù)值計算。
2) 各組層狀板巖荷載-位移曲線均具有線彈性-瞬時劈裂跌落的發(fā)展趨勢,試樣均呈現(xiàn)為顯著脆性劈裂破壞特征。
3) 當(dāng)層理傾角分別為0°,30°,45°,60°和90°時,層狀板巖峰值荷載依次為98.374,81.457,78.673,80.751 和102.595 kN,隨著傾角逐漸增加,層狀板巖峰值荷載大致呈現(xiàn)“U”型分布。
4) 隨層理傾角變化,層狀板巖劈裂破壞形態(tài)復(fù)雜多變,但主要呈現(xiàn)為沿層理面的滑剪破壞以及受基質(zhì)和層理綜合影響的拉剪混合破壞模式。
5) 基于固體鍵合顆粒模型的層狀板巖力學(xué)響應(yīng)及劈裂破壞形態(tài)與試驗結(jié)果具有較好的一致性,同時該方法還可再現(xiàn)不同層理傾角下層狀板巖的劈裂演化過程。