充填節(jié)理特性對(duì)巖石靜態(tài)壓縮力學(xué)行為影響的顆粒流模擬*

宋浪， 柴少波，， 宋博陽(yáng)， 柴連增， 杜宇翔

1.長(zhǎng)安大學(xué)建筑工程學(xué)院，陜西 西安 710064

2.江漢大學(xué)爆破工程湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430056

實(shí)際工程巖體中隨機(jī)分布著大量的節(jié)理，而節(jié)理往往含有一定的充填物。相比于非充填節(jié)理，充填節(jié)理由于裂隙貫通性好、充填材料強(qiáng)度低等特點(diǎn)，其力學(xué)強(qiáng)度和變形特性更差。在外荷載作用下充填節(jié)理層容易產(chǎn)生法向及切向變形，進(jìn)而易引起節(jié)理巖質(zhì)邊坡或地下洞室工程失穩(wěn)破壞。因此，研究荷載作用下充填節(jié)理巖石的力學(xué)及變形特性可為節(jié)理巖體工程穩(wěn)定性分析和安全性評(píng)價(jià)提供指導(dǎo)。

近年來(lái)對(duì)充填節(jié)理巖石力學(xué)及變形特性的相關(guān)研究持續(xù)升溫（Li et al.，2014；史玲等，2012；Liu et al.，2017）。Li et al.（2009）利用改進(jìn)的SHPB裝置研究發(fā)現(xiàn)充填節(jié)理的含水量和厚度對(duì)充填節(jié)理巖石動(dòng)態(tài)應(yīng)力-應(yīng)變曲線有明顯的影響。劉婷婷等（2017）采用RMT-150C試驗(yàn)機(jī)進(jìn)行充填節(jié)理的直接剪切試驗(yàn)，研究了法向壓力、剪切速率和充填厚度對(duì)其剪切力學(xué)特性的影響。許江等（2021）研究了循環(huán)荷載作用下充填厚度對(duì)結(jié)構(gòu)面剪切力學(xué)特性和破壞模式的影響。肖維民等（2019，2020）研究了節(jié)理粗糙度和充填度對(duì)薄層充填巖石節(jié)理抗剪強(qiáng)度、剪脹特性和破壞特征的影響。Chai et al.（2020）研究了不同充填厚度及不同充填材料對(duì)節(jié)理巖石試樣靜態(tài)和動(dòng)態(tài)力學(xué)特性的影響，研究發(fā)現(xiàn)充填節(jié)理的靜態(tài)和動(dòng)態(tài)抗壓強(qiáng)度，隨節(jié)理充填物強(qiáng)度的增大而增大，隨節(jié)理厚度的增大而減小；靜態(tài)和動(dòng)態(tài)壓縮下充填節(jié)理的變形特性與充填物材料性質(zhì)相關(guān)。此外柴少波等（2020，2022）還進(jìn)一步研究了累積沖擊和干濕循環(huán)作用對(duì)充填節(jié)理巖石力學(xué)特性的影響規(guī)律，一定程度上揭示了劣化作用下充填節(jié)理巖石力學(xué)特性的損傷規(guī)律。

以上研究均采用室內(nèi)物理試驗(yàn)的方法，但由于室內(nèi)試驗(yàn)中人工制備充填節(jié)理巖石試樣成本高、養(yǎng)護(hù)周期長(zhǎng)且試樣制作精度存在差異，故室內(nèi)物理試驗(yàn)存在一定的難度。此外，近年來(lái)隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，基于離散元法的數(shù)值計(jì)算軟件PFC 從散體介質(zhì)的微觀力學(xué)特征出發(fā)，構(gòu)建的數(shù)值模型可以真實(shí)逼近巖石材料的力學(xué)特性，反映節(jié)理巖體的非線性大變形特征，在模擬節(jié)理巖體此類非均質(zhì)、非連續(xù)材料方面具有很大的優(yōu)勢(shì)（王衛(wèi)華等，2005），受到了許多研究人員的青睞。例如，Zhang et al.（2012；2013）利用PFC2D程序研究了單軸壓縮過(guò)程中裂隙傾角對(duì)單裂隙和雙裂隙巖石裂紋擴(kuò)展情況的影響。金愛(ài)兵等（2012）進(jìn)行了巖石的雙軸試驗(yàn)?zāi)M，從巖橋長(zhǎng)度、節(jié)理長(zhǎng)度和節(jié)理傾角三方面對(duì)巖石的破裂模式和力學(xué)特性進(jìn)行了分析。黃彥華等（2014，2015）采用PFC2D研究了不同圍壓作用下含共面雙裂隙脆性砂巖的強(qiáng)度變化情況和裂紋演化規(guī)律，此外還研究了節(jié)理數(shù)量、節(jié)理傾角和圍壓對(duì)含兩組交叉節(jié)理巖體強(qiáng)度及破壞特征的影響。周喻等（2012）從宏觀和細(xì)觀角度深入探討了節(jié)理在直剪試驗(yàn)過(guò)程中的力學(xué)演化規(guī)律和破壞機(jī)制。張國(guó)凱等（2016）進(jìn)行了單軸壓縮下節(jié)理厚度和剛度比對(duì)試件宏觀力學(xué)參數(shù)、能量及波傳播規(guī)律影響的數(shù)值模擬，但缺乏與室內(nèi)物理試驗(yàn)結(jié)果的對(duì)比驗(yàn)證。許萬(wàn)忠等（2018）進(jìn)行了直剪試驗(yàn)數(shù)值模擬，研究JRC 系數(shù)、充填強(qiáng)度、接觸面黏結(jié)強(qiáng)度等參數(shù)對(duì)充填節(jié)理剪切強(qiáng)度的影響。連蓮等（2019）研究了不同充填厚度下充填節(jié)理的剪切破壞行為和強(qiáng)度特性。綜上所述，現(xiàn)有文獻(xiàn)關(guān)于節(jié)理巖石的顆粒流數(shù)值模擬中研究對(duì)象多為非充填節(jié)理，針對(duì)充填節(jié)理巖石的數(shù)值模擬研究較少，且其中剪切力學(xué)特性的研究居多，而壓縮力學(xué)特性的研究較為匱乏。

鑒于此，本文開(kāi)展充填節(jié)理特性對(duì)巖石靜態(tài)壓縮力學(xué)行為影響的顆粒流模擬。首先在PFC2D程序中構(gòu)建充填節(jié)理巖石數(shù)值模型，根據(jù)充填節(jié)理巖石的室內(nèi)物理靜態(tài)單軸壓縮試驗(yàn)結(jié)果進(jìn)行細(xì)觀參數(shù)標(biāo)定及合理性驗(yàn)證，其次進(jìn)一步研究不同節(jié)理傾角、多條節(jié)理數(shù)以及不同節(jié)理厚度情況下充填節(jié)理巖石的力學(xué)參數(shù)變化規(guī)律、微裂紋擴(kuò)展和破壞形態(tài)情況。研究成果可豐富充填節(jié)理力學(xué)特性研究，并為復(fù)雜分布充填節(jié)理巖體工程穩(wěn)定性分析提供參考。

1 數(shù)值模型建立及驗(yàn)證

本文前期研究（柴少波等，2022）的室內(nèi)靜態(tài)單軸壓縮試驗(yàn)，所用充填節(jié)理巖石試樣直徑為50 mm、高為35 mm（圖1）。為對(duì)比需要，本文在PFC中建立相同尺寸的數(shù)值模型（圖2）。

圖1 充填節(jié)理巖石實(shí)物圖Fig.1 Photos of filled jointed rock samples

圖2 充填節(jié)理巖石顆粒流數(shù)值模型Fig.2 Particle flow numerical model of filled jointed rock

PFC 程序中內(nèi)置接觸模型有多種：線性模型、接觸黏結(jié)模型、平行黏結(jié)模型、平直節(jié)理模型和光滑節(jié)理模型等，其中平行黏結(jié)模型可以很好地模擬巖石材料的力學(xué)行為（黃彥華等，2014；LEE et al.，2011），本數(shù)值模型采用平行黏結(jié)模型。平行黏結(jié)模型包括兩種接觸界面：其一是無(wú)限小的線彈性界面，該界面可以承受摩擦力，但無(wú)法承受張力和傳遞力矩，該界面模型等效于線性模型；其二是有具體尺寸的線彈性黏結(jié)界面，可以傳遞力和力矩，稱之為平行黏結(jié)。當(dāng)?shù)诙N界面黏結(jié)時(shí)，能抵抗扭矩且表現(xiàn)為線彈性，當(dāng)所受的作用力大于強(qiáng)度極限時(shí)，黏結(jié)模型失效，從而無(wú)法傳遞荷載，退化為線性模型。此外，由于數(shù)值模型中顆粒粒徑大小和數(shù)量會(huì)對(duì)宏觀力學(xué)參數(shù)產(chǎn)生一定影響（周喻等，2011），所以需要考慮模型最小尺度上顆粒數(shù)RES（resolution）對(duì)宏觀參數(shù)的影響，官方Itasca建議RES＞20，并定義

式中L為數(shù)值模型的最小尺度，Rmax為顆粒最大粒徑，Rmin為顆粒最小粒徑。

考慮計(jì)算效率和文獻(xiàn)中顆粒粒徑的選?。ò⒈葼柕牡龋?018；田文嶺等，2017；黃達(dá)等，2013），本文設(shè)置中Rmin=0.3 mm，Rmax=0.5 mm，根據(jù)式（1）計(jì)算得到本數(shù)值模型RES=43.9，滿足RES＞20的要求。數(shù)值模型共生成了顆粒3 370個(gè)，接觸7 675個(gè)，孔隙率為1.2%，顆粒密度為2 700 kg/m3。為了準(zhǔn)確模擬室內(nèi)靜態(tài)單軸試驗(yàn)的加載過(guò)程，上下墻的速度設(shè)定為0.05 m/s（田文嶺等，2017；楊圣奇等，2022），當(dāng)應(yīng)力跌至峰值應(yīng)力的70%時(shí)停止加載。

細(xì)觀參數(shù)標(biāo)定是建立準(zhǔn)確的顆粒流數(shù)值模型的關(guān)鍵，微觀參數(shù)標(biāo)定是一個(gè)非常繁雜的過(guò)程，PFC 軟件推薦使用“試錯(cuò)法”進(jìn)行細(xì)觀參數(shù)標(biāo)定，“試錯(cuò)法”本質(zhì)上就是通過(guò)不斷調(diào)試細(xì)觀參數(shù)，使得數(shù)值模擬結(jié)果和實(shí)際物理試驗(yàn)結(jié)果接近，從而達(dá)到逼近宏觀試驗(yàn)的目的。在細(xì)觀參數(shù)標(biāo)定過(guò)程中需綜合考慮充填節(jié)理巖石試樣以及各自分層的力學(xué)參數(shù)，巖石層和充填節(jié)理層細(xì)觀參數(shù)的比例關(guān)系參考實(shí)際物理實(shí)驗(yàn)兩者力學(xué)參數(shù)的比例關(guān)系。此外由于實(shí)際工程巖體中的節(jié)理充填層是由巖石碎屑分化形成的顆粒與具有一定黏性的礦物膠結(jié)形成的混合物，在沉積作用過(guò)程中節(jié)理充填層與巖石會(huì)存在一定的黏結(jié)力，且黏結(jié)力的大小與節(jié)理充填層內(nèi)部的黏結(jié)力基本相同。故本數(shù)值模型中巖石層和充填節(jié)理層之間的接觸也采用平行黏結(jié)模型，其細(xì)觀參數(shù)與充填節(jié)理層組內(nèi)顆粒間接觸的細(xì)觀參數(shù)保持一致。本研究中充填節(jié)理巖石數(shù)值模擬與室內(nèi)物理試驗(yàn)的應(yīng)力-應(yīng)變曲線對(duì)比結(jié)果如圖3（a）所示?？梢钥闯觯瑪?shù)值模擬與室內(nèi)試驗(yàn)的單軸壓縮應(yīng)力-應(yīng)變曲線除了前期不同，后期相似度很高，這是由于PFC2D程序中數(shù)值模型在初始狀態(tài)時(shí)通過(guò)伺服機(jī)制使得模型已經(jīng)達(dá)到了密實(shí)狀態(tài)，故無(wú)法模擬壓縮過(guò)程中巖樣的初始?jí)好茈A段（黃彥華等，2015；張寶玉等，2021；于利強(qiáng)等，2021）.所以在進(jìn)行數(shù)值模擬結(jié)果和室內(nèi)試驗(yàn)結(jié)果比較分析時(shí)，忽略室內(nèi)試驗(yàn)壓密階段充填節(jié)理層產(chǎn)生的壓密變形量，主要對(duì)比峰值應(yīng)力和彈性階段時(shí)彈性模量的差異以及最終破壞形態(tài)。數(shù)值模擬與室內(nèi)物理試驗(yàn)壓縮破壞形態(tài)結(jié)果對(duì)比如圖3（b）所示，可發(fā)現(xiàn)數(shù)值模型的最終破壞形態(tài)與實(shí)際物理試驗(yàn)充填節(jié)理巖石的破壞形態(tài)較為吻合。

圖3 數(shù)值模擬與室內(nèi)試驗(yàn)結(jié)果對(duì)比圖Fig.3 Comparison of numerical simulation and laboratory test results

通過(guò)“試錯(cuò)法”標(biāo)定后的數(shù)值模型細(xì)觀參數(shù)見(jiàn)表1。進(jìn)一步對(duì)比數(shù)值模擬與室內(nèi)物理試驗(yàn)結(jié)果的峰值應(yīng)力和彈性模量數(shù)值，物理試驗(yàn)結(jié)果中峰值應(yīng)力和彈性模量分別為47.35和0.845 GPa；數(shù)值模擬結(jié)果中峰值應(yīng)力和彈性模量分別為45.26 和0.851 GPa。可以發(fā)現(xiàn)峰值應(yīng)力的相對(duì)誤差＜ 4.5%，彈性模量的相對(duì)誤差＜ 1%，數(shù)值結(jié)果吻合度良好，表示該數(shù)值模型能夠較為真實(shí)模擬實(shí)際中的充填節(jié)理巖石，同時(shí)也驗(yàn)證了本文數(shù)值模型和細(xì)觀參數(shù)選取的正確性和可靠性。

表1 充填節(jié)理巖石數(shù)值模型的細(xì)觀參數(shù)表Table 1 Mesoscopic parameter table of numerical model of filled jointed rock

為更好展示單軸壓縮過(guò)程中微裂紋的發(fā)育情況，每間隔10 000 時(shí)步（step）記錄數(shù)值模型在單軸壓縮過(guò)程中的微裂紋擴(kuò)展情況（圖4）。由圖4 可以看出，在時(shí)步為20 000 時(shí)充填節(jié)理層出現(xiàn)首個(gè)微裂紋；之后微裂紋不斷產(chǎn)生并且多出現(xiàn)于充填節(jié)理層的外側(cè)區(qū)域；隨著時(shí)步的增加，充填節(jié)理層的微裂紋不斷向中部區(qū)域擴(kuò)展，同時(shí)還可以看出充填節(jié)理層外側(cè)顆粒被擠壓移位并發(fā)生“逃逸”現(xiàn)象。隨著壓縮過(guò)程的持續(xù)進(jìn)行，充填節(jié)理層遍布微裂紋且兩側(cè)巖石開(kāi)始出現(xiàn)微裂紋，兩側(cè)巖石的微裂紋不斷擴(kuò)展和匯聚形成肉眼可見(jiàn)的裂縫，繼而裂縫不斷沿徑向發(fā)展，最終形成貫通裂縫導(dǎo)致充填節(jié)理巖石發(fā)生破壞。這與室內(nèi)試驗(yàn)觀察結(jié)果一致（圖5）。由圖5 可以看出，首先是充填節(jié)理層被壓縮，充填節(jié)理層出現(xiàn)一些細(xì)小裂紋并且充填節(jié)理邊緣部位向外隆起并伴有碎屑脫落。隨著荷載的施加，充填節(jié)理層裂紋不斷發(fā)展并且充填節(jié)理層達(dá)到壓密狀態(tài)，兩側(cè)巖石開(kāi)始被壓縮并出現(xiàn)微小裂紋。持續(xù)加載后兩側(cè)巖石裂紋將繼續(xù)擴(kuò)展、貫通，最終導(dǎo)致充填節(jié)理巖石失穩(wěn)破壞。通過(guò)比對(duì)單軸壓縮過(guò)程中數(shù)值模型和室內(nèi)物理試樣的破壞全過(guò)程可以看出PFC2D程序可以準(zhǔn)確地模擬充填節(jié)理巖石的破壞過(guò)程，進(jìn)一步驗(yàn)證了顆粒流數(shù)值模擬的可行性。

圖4 單軸壓縮全過(guò)程數(shù)值模型的裂紋發(fā)育情況Fig.4 Crack development in numerical model of uniaxial compression process

圖5 單軸壓縮全過(guò)程充填節(jié)理巖石試樣破壞情況Fig.5 Failure of filled jointed rock in the whole process of uniaxial compression

2 結(jié)果分析

2.1 不同節(jié)理傾角

對(duì)不同充填節(jié)理傾角（0°、15°、30°、45°、60°、75°和90°）下數(shù)值模型進(jìn)行單軸壓縮數(shù)值模擬，如圖6 為單軸壓縮過(guò)程的應(yīng)力-應(yīng)變曲線。從圖6中可以看出，不同節(jié)理傾角的充填節(jié)理巖石對(duì)應(yīng)的應(yīng)力-應(yīng)變曲線差異很大，隨著節(jié)理傾角的增大，應(yīng)力-應(yīng)變曲線先出現(xiàn)下移的趨勢(shì)，其中節(jié)理傾角為30°時(shí)應(yīng)力-應(yīng)變曲線處于最下方，隨著節(jié)理傾角的持續(xù)增大，應(yīng)力-應(yīng)變曲線又逐漸向上方移動(dòng)，節(jié)理傾角為90°時(shí)應(yīng)力-應(yīng)變曲線處于最上方。進(jìn)一步分析峰值應(yīng)力、彈性模量與節(jié)理傾角的關(guān)系，如圖7為不同節(jié)理傾角下峰值應(yīng)力和彈性模量變化曲線。可以發(fā)現(xiàn)峰值應(yīng)力的變化規(guī)律呈現(xiàn)出明顯的U 型分布，彈性模量的變化規(guī)律呈現(xiàn)出近似的U型分布，以往文獻(xiàn)研究（陳新等，2011；鄧華鋒等，2019）中在不同節(jié)理傾角下非充填的單節(jié)理或多節(jié)理巖石峰值應(yīng)力和彈性模同樣也呈現(xiàn)出U 型分布特征，表明充填節(jié)理的研究結(jié)果與非充填節(jié)理巖石較為一致。

圖6 不同節(jié)理傾角下數(shù)值模型的應(yīng)力-應(yīng)變曲線Fig.6 Stress-strain curves of numerical models under different joint angles

圖7 不同節(jié)理傾角下峰值應(yīng)力和彈性模量變化曲線Fig.7 Curve of peak stress and elastic modulus under different joint angles

此外，從圖7 中還可以看出，節(jié)理傾角為90°時(shí)，峰值應(yīng)力和彈性模量最高；節(jié)理傾角為30°時(shí)，峰值應(yīng)力和彈性模量最低。節(jié)理傾角90°時(shí)峰值應(yīng)力和彈性模量分別為55.3 和2.348 GPa，較0°傾角時(shí)分別提高了22.18%和153.84%。節(jié)理傾角為30°時(shí)峰值應(yīng)力和彈性模量分別為9.59 和0.539 GPa，較0°傾角時(shí)分別下降了78.81%和41.73%。

為研究壓縮過(guò)程中不同節(jié)理傾角下微裂紋的發(fā)育情況，提取不同節(jié)理傾角下微裂紋數(shù)量隨應(yīng)變的變化曲線，如圖8所示。結(jié)果發(fā)現(xiàn)不同節(jié)理傾角下微裂紋的變化曲線差異較大。當(dāng)充填節(jié)理處于水平位置時(shí)數(shù)值模型的微裂紋發(fā)育主要呈現(xiàn)出前期增長(zhǎng)較快、中期增長(zhǎng)緩慢、后期迅速增長(zhǎng)的3階段特點(diǎn)。對(duì)照?qǐng)D4 對(duì)所呈現(xiàn)的微裂紋3 階段發(fā)育特點(diǎn)的分析如下：前期裂紋增長(zhǎng)較快的主要原因是充填節(jié)理層強(qiáng)度及模量較低，在壓力作用下更易產(chǎn)生裂紋且低強(qiáng)度材料中裂紋擴(kuò)展較快；中期裂紋增長(zhǎng)緩慢的主要原因是充填節(jié)理層裂紋不再增加，裂紋產(chǎn)生位置逐步轉(zhuǎn)移至兩側(cè)巖石層；后期迅速增長(zhǎng)的主要原因是兩側(cè)巖石層裂紋在壓力作用下不斷匯聚和擴(kuò)展，充填節(jié)理試樣趨于失穩(wěn)破壞。

從圖8中還可以發(fā)現(xiàn)當(dāng)節(jié)理傾角改變時(shí)微裂紋發(fā)育情況也發(fā)生改變，當(dāng)節(jié)理傾角為15°時(shí)數(shù)值模型的微裂紋發(fā)育還近似呈現(xiàn)為3階段特點(diǎn)，但可以看出微裂紋變化曲線出現(xiàn)明顯向左偏移現(xiàn)象，表明同一應(yīng)變下節(jié)理傾角為15°數(shù)值模型的微裂紋數(shù)量更多；當(dāng)節(jié)理傾角增加至30°時(shí)數(shù)值模型的微裂紋發(fā)育不再呈現(xiàn)典型的3階段特點(diǎn)，而是呈現(xiàn)出前期增長(zhǎng)緩慢、后期增長(zhǎng)迅速的2階段特點(diǎn)，隨著節(jié)理傾角的持續(xù)增加，節(jié)理傾角為45°、60°和75°時(shí)數(shù)值模型的微裂紋發(fā)育又近似呈現(xiàn)出3 個(gè)階段特點(diǎn)，但微裂紋發(fā)育曲線歷程明顯縮短，而當(dāng)節(jié)理傾角為90°時(shí)，微裂紋發(fā)育情況與節(jié)理傾角為30°時(shí)類似，上述現(xiàn)象說(shuō)明節(jié)理傾角改變會(huì)嚴(yán)重影響充填節(jié)理巖石裂紋的發(fā)育。進(jìn)一步提取最終破壞時(shí)不同節(jié)理傾角下微裂紋總數(shù)量，如圖9所示。從圖9中可知充填節(jié)理巖石在不同節(jié)理傾角下最終破壞時(shí)的總裂紋數(shù)分別為1 030、990、839、774、585、865和901，可知隨著節(jié)理傾角的變化，最終破壞時(shí)總裂紋數(shù)呈現(xiàn)出先減小后增大的變化趨勢(shì)。

圖9 不同節(jié)理傾角下最終破壞時(shí)微裂紋總數(shù)量Fig.9 Total number of microcracks at final failure under different joint angles

引入起裂應(yīng)力進(jìn)行深入分析，將單軸壓縮過(guò)程中微裂紋數(shù)為峰值應(yīng)力時(shí)微裂紋數(shù)的1%對(duì)應(yīng)的軸向應(yīng)力定義為起裂應(yīng)力（Potyondy et al.，2004），如圖10 所示。進(jìn)一步可以得到不同節(jié)理傾角下起裂應(yīng)力的變化情況，如圖11 所示。從圖11 中可以看出，當(dāng)節(jié)理傾角不同時(shí)，充填節(jié)理巖石數(shù)值模型的起裂應(yīng)力分別為8.33、7.52、4.06、5.46、10.65、17.24 和30.61 MPa。充填節(jié)理巖石數(shù)值模型的起裂應(yīng)力隨節(jié)理傾角的增大出現(xiàn)先減小后增大的趨勢(shì)，也呈現(xiàn)出近似U 型分布的特點(diǎn)，其中起裂應(yīng)力最小時(shí)對(duì)應(yīng)的節(jié)理傾角為30°，起裂應(yīng)力最大時(shí)對(duì)應(yīng)的節(jié)理傾角為90°，也就表明當(dāng)節(jié)理傾角為30°時(shí)充填節(jié)理巖石最容易率先產(chǎn)生裂紋，而節(jié)理傾角為90°時(shí)充填節(jié)理巖石最不易產(chǎn)生裂紋。

圖10 起裂應(yīng)力示意圖Fig.10 Diagram of initiation stress

圖11 不同節(jié)理傾角下起裂應(yīng)力變化曲線Fig.11 Fracture initiation stress curve under different joint angles

進(jìn)一步通過(guò)圖12 觀察不同節(jié)理傾角下充填節(jié)理巖石最終破壞形態(tài)可以看出，節(jié)理傾角會(huì)顯著影響充填節(jié)理巖石的破壞形態(tài)。當(dāng)充填節(jié)理處于水平位置時(shí)，數(shù)值模型最終破壞形態(tài)表現(xiàn)為充填節(jié)理層發(fā)生很大的豎向壓縮并且兩側(cè)巖石出現(xiàn)貫通裂縫；當(dāng)節(jié)理傾角為15°時(shí)，充填節(jié)理層產(chǎn)生了較大的豎向壓縮且數(shù)值模型的上側(cè)巖石出現(xiàn)滑移現(xiàn)象，巖石層微裂紋很少，微裂紋基本集中在充填節(jié)理層；隨著節(jié)理傾角的增大，充填節(jié)理層壓縮量減小，當(dāng)節(jié)理傾角為60°時(shí)，充填節(jié)理巖石總裂紋數(shù)明顯減小，充填節(jié)理層與巖石層出現(xiàn)一定程度的脫離，巖石層出現(xiàn)一條貫通大裂縫最終導(dǎo)致試樣破壞；隨著節(jié)理傾角的繼續(xù)增大，巖石層的微裂紋數(shù)量持續(xù)增多，當(dāng)節(jié)理傾角為90°時(shí)，可以看出巖石層的破壞最為嚴(yán)重，由上文所述可知節(jié)理傾角為90°時(shí)數(shù)值模型的峰值應(yīng)力和起裂應(yīng)力最大，由此可以推斷出當(dāng)節(jié)理傾角較小時(shí)，豎向壓力作用主要由充填節(jié)理層和兩側(cè)巖石層共同承擔(dān)，而當(dāng)節(jié)理傾角較大時(shí)，豎向壓力作用主要由兩側(cè)巖石層承擔(dān)。

圖12 不同節(jié)理傾角下充填節(jié)理巖石最終破壞形態(tài)Fig.12 Final failure pattern of filled jointed rock under different joint angles

2.2 不同充填節(jié)理厚度

圖13 為5 種不同充填厚度（2.0、3.5、5.0、6.5和8.0 mm）的充填節(jié)理巖石單軸壓縮的應(yīng)力-應(yīng)變曲線?？梢钥闯?，應(yīng)力-應(yīng)變曲線變化趨勢(shì)較有規(guī)律，隨著充填節(jié)理厚度的增加，充填節(jié)理巖石的應(yīng)力-應(yīng)變曲線逐漸向下方移動(dòng)，其變化趨勢(shì)與試驗(yàn)結(jié)果一致（Chai et al.，2020）。充填節(jié)理厚度分別為6.5 和8.0 mm 時(shí)應(yīng)力-應(yīng)變曲線峰前出現(xiàn)波動(dòng)現(xiàn)象，這是由于充填節(jié)理厚度較大時(shí)在壓縮過(guò)程中充填層顆粒會(huì)發(fā)生錯(cuò)動(dòng)和向外擠壓情況。充填厚度改變時(shí)充填節(jié)理巖石的峰值應(yīng)力和彈性模量的變化規(guī)律如圖14所示。

圖13 不同充填厚度下數(shù)值模型的應(yīng)力-應(yīng)變曲線Fig.13 Stress-strain curves of numerical models under different filling thicknesses

圖14 不同充填厚度下峰值應(yīng)力和彈性模量變化曲線Fig.14 Curve of peak stress and elastic modulus under different filling thickness

從圖14 可知當(dāng)充填節(jié)理厚度不斷增加時(shí)，數(shù)值模型的峰值應(yīng)力分別為57.82、 51.58、45.26、32.08 和23.83 MPa，彈性模量分別為1.31、1.103、0.851、0.809和0.754 GPa。充填厚度增加至3.5mm時(shí)峰值應(yīng)力和彈性模量分別降低了10.79% 和15.8%；充填厚度增加至5.0 mm時(shí)峰值應(yīng)力和彈性模量分別降低了21.72%和35.04%；充填厚度增加至6.5 mm 時(shí)峰值應(yīng)力和彈性模量分別降低了44.52%和38.24%；充填厚度增加至8.0 mm 時(shí)峰值應(yīng)力和彈性模量分別降低了58.79%和42.44%?？梢园l(fā)現(xiàn)峰值應(yīng)力和彈性模量隨充填厚度的增加而逐漸降低，并且彈性模量降低速度逐漸變緩，進(jìn)一步通過(guò)擬合發(fā)現(xiàn)峰值應(yīng)力的降低趨勢(shì)符合線性函數(shù)分布，而彈性模量的降低趨勢(shì)符合多項(xiàng)式函數(shù)分布，相關(guān)系數(shù)均大于0.97，說(shuō)明一致性良好。

深入分析不同充填厚度情況下微裂紋的發(fā)育情況，如圖15 所示。從圖15 中可以看出不同厚度充填節(jié)理巖石數(shù)值模型微裂紋的發(fā)育過(guò)程也呈現(xiàn)典型的3 階段特點(diǎn)，但隨著充填節(jié)理厚度的增加，微裂紋數(shù)隨應(yīng)變的變化曲線逐漸向上側(cè)移動(dòng)，從微裂紋發(fā)育的第2階段可以明顯看出微裂紋數(shù)量從高到低依次為8.0、6.5、5.0、3.5、2.0 mm，且充填節(jié)理厚度較大時(shí)第2 階段的微裂紋數(shù)量明顯增多。從圖16 可知當(dāng)充填節(jié)理厚度不斷增加時(shí)，最終破壞時(shí)數(shù)值模型的總裂紋數(shù)分別為744、903、1 030、1 251 和1 399，最終破壞時(shí)數(shù)值模型的總裂紋數(shù)隨著充填節(jié)理厚度的增加而不斷增加，通過(guò)擬合發(fā)現(xiàn)總裂紋數(shù)的增長(zhǎng)趨勢(shì)符合線性函數(shù)分布，且一致性良好。進(jìn)一步計(jì)算不同充填厚度下充填節(jié)理巖石的起裂應(yīng)力數(shù)值，結(jié)果發(fā)現(xiàn)隨著充填節(jié)理厚度的增加時(shí)起裂應(yīng)力分別為10.12、7.52、8.19、6.02 和6.32 MPa，可以看出不同厚度充填節(jié)理巖石的起裂應(yīng)力基本處于6～10 MPa 范圍之內(nèi)，說(shuō)明充填節(jié)理厚度改變對(duì)起裂應(yīng)力影響不大。

圖15 不同充填厚度下微裂紋數(shù)量變化曲線Fig.15 Variation curve of the number of microcracks under different filling thicknesses

圖16 不同充填厚度下最終破壞時(shí)微裂紋總數(shù)量Fig.16 Total number of microcracks at final failure under different filling thicknesses

如圖17 為不同厚度充填節(jié)理巖石數(shù)值模型的最終破壞形態(tài)。可以發(fā)現(xiàn)隨著充填節(jié)理厚度的增加，最終破壞時(shí)充填節(jié)理層顆粒受擠壓向兩側(cè)外溢、脫落情況嚴(yán)重，同時(shí)還發(fā)現(xiàn)兩側(cè)巖石層的微裂紋較少，微裂紋多集中于充填節(jié)理層，說(shuō)明壓縮過(guò)程中充填節(jié)理層受損程度更嚴(yán)重，此外還可以看出不同厚度充填節(jié)理巖石數(shù)值模型最終破壞模式為張拉破壞，可以表明改變充填節(jié)理厚度對(duì)充填節(jié)理巖石壓縮破壞模式并無(wú)影響。

圖17 不同充填厚度下充填節(jié)理巖石最終破壞形態(tài)Fig.17 Final failure pattern of filled jointed rock under different filling thickness

2.3 一組平行分布充填節(jié)理

實(shí)際自然界巖體內(nèi)部中節(jié)理是成組出現(xiàn)的，為進(jìn)一步對(duì)比1組節(jié)理與單個(gè)節(jié)理情況下力學(xué)及變形特性的差異性，在數(shù)值模型建立過(guò)程中設(shè)置1組平行分布的充填節(jié)理（3 條節(jié)理），并且控制充填節(jié)理的總厚度5 mm 不變，其中充填節(jié)理層和巖石層的細(xì)觀參數(shù)與上文一致。以此來(lái)研究含1組平行充填節(jié)理情況下節(jié)理巖石的力學(xué)及變形特性，如圖18為不同節(jié)理傾角下含1組平行充填節(jié)理數(shù)值模型的應(yīng)力-應(yīng)變曲線。進(jìn)一步可以得到含1 組平行充填節(jié)理的巖石試樣在不同節(jié)理傾角下的峰值應(yīng)力和彈性模量，如圖19所示。

圖18 不同節(jié)理傾角下含1組平行充填節(jié)理的應(yīng)力-應(yīng)變曲線Fig.18 Stress-strain curves of a group of parallel filled joint under different joint angles

圖19 不同節(jié)理傾角下含1組平行充填節(jié)理峰值應(yīng)力和彈性模量變化曲線Fig.19 Curves of peak stress and elastic modulus of a group of parallel filled joint under different joint angles

由圖18 中可以看出，含1 組平行充填節(jié)理的巖石試樣的應(yīng)力-應(yīng)變曲線隨著節(jié)理傾角的增大出現(xiàn)先向左下方移動(dòng)后向左上方移動(dòng)的現(xiàn)象，與含單個(gè)充填節(jié)理的巖石應(yīng)力-應(yīng)變曲線隨節(jié)理傾角的變化趨勢(shì)相同。不同的是，含1組平行充填節(jié)理的巖石試樣在節(jié)理傾角為45°時(shí)的應(yīng)力-應(yīng)變曲線處于最下方。由圖19 可以發(fā)現(xiàn)，隨著節(jié)理傾角的增加，峰值應(yīng)力和彈性模量均出現(xiàn)先減小后增大的趨勢(shì)。相比于含單個(gè)充填節(jié)理的巖石，含1組平行充填節(jié)理的巖石在節(jié)理傾角為0°時(shí)峰值應(yīng)力差異較小；節(jié)理傾角為15°和30°時(shí)峰值應(yīng)力分別提高了38.3%和104.9%；節(jié)理傾角為45°、60°、75°和90°時(shí)峰值應(yīng)力分別降低了38.4%、62.3%、46.5%和23.3%；可以發(fā)現(xiàn)當(dāng)節(jié)理傾角小于45°時(shí)含1 組平行充填節(jié)理的巖石峰值應(yīng)力較高，而節(jié)理傾角大于45°時(shí)含單個(gè)充填節(jié)理的巖石峰值應(yīng)力較高，其主要原因是當(dāng)節(jié)理傾角較小時(shí)，如圖20 所示可以看出含單個(gè)充填節(jié)理的巖石數(shù)值模型顆粒間接觸力分布不均勻，充填節(jié)理層兩側(cè)位置處的接觸力較大，而含1組平行充填節(jié)理的巖石數(shù)值模型顆粒間接觸力分布較為均勻。所以相比于充填節(jié)理集中布設(shè)在巖石中部，充填節(jié)理分散分布于巖石試樣內(nèi)部可以使得壓力作用下充填節(jié)理巖石整體受力更為均勻，可以減小顆粒間接觸斷裂的可能性。而當(dāng)節(jié)理傾角較大時(shí)，由于多條充填節(jié)理的存在，在壓縮過(guò)程中充填層與巖石層易發(fā)生分離脫落現(xiàn)象，從而降低了充填節(jié)理巖石整體的強(qiáng)度及穩(wěn)定性。

圖20 充填節(jié)理巖石數(shù)值模型顆粒間接觸力分布圖Fig.20 Distribution of contact forces between particles in numerical model of filled jointed rock

進(jìn)一步分析不同節(jié)理傾角下含1組平行充填節(jié)理的微裂紋變化曲線，如圖21 所示。可以看出節(jié)理傾角為0°時(shí)微裂紋變化曲線同樣呈現(xiàn)出典型的3階段特點(diǎn)。隨著節(jié)理傾角的增加，微裂紋發(fā)育的變化情況與含單個(gè)充填節(jié)理的巖石變化情況類似。圖22 為不同節(jié)理傾角下含1 組平行充填節(jié)理最終破壞時(shí)微裂紋總數(shù)量，可以看出隨著節(jié)理傾角的增加，最終破壞時(shí)總裂紋數(shù)也呈現(xiàn)出先減小后增大的趨勢(shì)，通過(guò)對(duì)比含單個(gè)節(jié)理的巖石數(shù)值模型的總裂紋數(shù)發(fā)現(xiàn)兩者相差較小，說(shuō)明充填節(jié)理數(shù)對(duì)單軸壓縮過(guò)程中微裂紋情況影響較小。

圖22 不同節(jié)理傾角下含1組平行充填節(jié)理最終破壞時(shí)微裂紋總數(shù)量Fig.22 The total number of microcracks in a group of parallel filling joints under different joint angles

圖23 為不同節(jié)理傾角下含1 組平行充填節(jié)理的巖石最終破壞形態(tài)。由圖23 可以看出，不同節(jié)理傾角下數(shù)值模型最終破壞形態(tài)差異較大，節(jié)理傾角為0°和15°時(shí)，充填節(jié)理層之間的巖石層均出現(xiàn)了貫通裂縫，節(jié)理傾角為30°時(shí)最終破壞形態(tài)完整性較好，當(dāng)節(jié)理傾角大于45°后充填節(jié)理層出現(xiàn)分離剝落現(xiàn)象，尤其當(dāng)節(jié)理傾角為75°時(shí)分離剝落現(xiàn)象最為嚴(yán)重，充填節(jié)理層均發(fā)生了明顯的滑移破壞，可知當(dāng)節(jié)理傾角較大且?guī)r石內(nèi)部含有多條充填節(jié)理時(shí)整體穩(wěn)定性較差、破壞程度嚴(yán)重。

圖23 不同節(jié)理傾角下含1組平行充填節(jié)理最終破壞形態(tài)Fig.23 Final failure form of a group of parallel filling joints under different joint angles

3 結(jié) 論

本文利用PFC2D程序構(gòu)建了充填節(jié)理巖石數(shù)值模型，結(jié)合室內(nèi)靜態(tài)單軸壓縮試驗(yàn)結(jié)果驗(yàn)證了細(xì)觀參數(shù)選取的可靠性以及顆粒流模擬的可行性，分析了不同節(jié)理傾角、多條節(jié)理數(shù)以及不同節(jié)理厚度情況下充填節(jié)理巖石的力學(xué)參數(shù)變化規(guī)律、微裂紋擴(kuò)展和破壞形態(tài)情況，主要結(jié)論如下：

1）隨著充填節(jié)理傾角的改變，峰值應(yīng)力和彈性模量呈現(xiàn)出U 型分布，最終破壞時(shí)的總裂紋數(shù)和起裂應(yīng)力呈現(xiàn)出先減小后增大的變化趨勢(shì)，此外改變充填節(jié)理傾角會(huì)嚴(yán)重影響充填節(jié)理巖石裂紋的發(fā)育狀況和最終破壞形態(tài)。

2）充填節(jié)理巖石的峰值應(yīng)力和彈性模量隨充填厚度的增加而逐漸降低，且彈性模量降低速度逐漸變緩，而最終破壞時(shí)的總裂紋數(shù)隨著充填厚度的增加而不斷增加，進(jìn)一步擬合發(fā)現(xiàn)峰值應(yīng)力和總裂紋數(shù)的變化規(guī)律符合線性函數(shù)分布，彈性模量的變化規(guī)律符合多項(xiàng)式函數(shù)分布。

3）不同充填厚度的節(jié)理巖石微裂紋的發(fā)育過(guò)程均呈現(xiàn)典型的3階段特點(diǎn)，即微裂紋前期增長(zhǎng)較快、中期增長(zhǎng)緩慢、后期迅速增長(zhǎng)；但充填節(jié)理厚度改變對(duì)起裂應(yīng)力和破壞模式影響不大。

4）巖石內(nèi)部存在1組平行充填節(jié)理時(shí)，隨著節(jié)理傾角的增加，峰值應(yīng)力和彈性模量均出現(xiàn)先減小后增大的趨勢(shì)，最終破壞時(shí)總裂紋數(shù)也呈現(xiàn)出先減小后增大的趨勢(shì)。但當(dāng)節(jié)理傾角較大會(huì)發(fā)生明顯的滑移破壞，充填節(jié)理巖石整體穩(wěn)定性較差、破壞程度嚴(yán)重。

本文進(jìn)行的數(shù)值模擬研究是基于二維顆粒流程序，二維顆粒流數(shù)值模擬未能反映節(jié)理面的空間形態(tài)特性和三維應(yīng)力狀態(tài)等真實(shí)情況，存在一定的不足和缺陷，但基于二維顆粒流程序的數(shù)值模擬一定程度上可以反映實(shí)際規(guī)律。此外，實(shí)際工程中充填節(jié)理分布復(fù)雜且充填節(jié)理材料差異明顯，本文研究中暫未考慮相關(guān)因素，僅是對(duì)節(jié)理傾角、充填厚度等因素展開(kāi)分析，研究結(jié)果在一定程度上可為后續(xù)分析復(fù)雜分布節(jié)理巖體穩(wěn)定性提供參考。