充填節理特性對巖石靜態壓縮力學行為影響的顆粒流模擬*

宋浪， 柴少波，， 宋博陽， 柴連增， 杜宇翔

1.長安大學建筑工程學院，陜西 西安 710064

2.江漢大學爆破工程湖北省重點實驗室，湖北 武漢 430056

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

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

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

鑒于此，本文開展充填節理特性對巖石靜態壓縮力學行為影響的顆粒流模擬。首先在PFC2D程序中構建充填節理巖石數值模型，根據充填節理巖石的室內物理靜態單軸壓縮試驗結果進行細觀參數標定及合理性驗證，其次進一步研究不同節理傾角、多條節理數以及不同節理厚度情況下充填節理巖石的力學參數變化規律、微裂紋擴展和破壞形態情況。研究成果可豐富充填節理力學特性研究，并為復雜分布充填節理巖體工程穩定性分析提供參考。

1 數值模型建立及驗證

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

圖1 充填節理巖石實物圖Fig.1 Photos of filled jointed rock samples

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

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

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

考慮計算效率和文獻中顆粒粒徑的選取（阿比爾的等，2018；田文嶺等，2017；黃達等，2013），本文設置中Rmin=0.3 mm，Rmax=0.5 mm，根據式（1）計算得到本數值模型RES=43.9，滿足RES＞20的要求。數值模型共生成了顆粒3 370個，接觸7 675個，孔隙率為1.2%，顆粒密度為2 700 kg/m3。為了準確模擬室內靜態單軸試驗的加載過程，上下墻的速度設定為0.05 m/s（田文嶺等，2017；楊圣奇等，2022），當應力跌至峰值應力的70%時停止加載。

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

圖3 數值模擬與室內試驗結果對比圖Fig.3 Comparison of numerical simulation and laboratory test results

通過“試錯法”標定后的數值模型細觀參數見表1。進一步對比數值模擬與室內物理試驗結果的峰值應力和彈性模量數值，物理試驗結果中峰值應力和彈性模量分別為47.35和0.845 GPa；數值模擬結果中峰值應力和彈性模量分別為45.26 和0.851 GPa。可以發現峰值應力的相對誤差＜ 4.5%，彈性模量的相對誤差＜ 1%，數值結果吻合度良好，表示該數值模型能夠較為真實模擬實際中的充填節理巖石，同時也驗證了本文數值模型和細觀參數選取的正確性和可靠性。

表1 充填節理巖石數值模型的細觀參數表Table 1 Mesoscopic parameter table of numerical model of filled jointed rock

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

圖4 單軸壓縮全過程數值模型的裂紋發育情況Fig.4 Crack development in numerical model of uniaxial compression process

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

2 結果分析

2.1 不同節理傾角

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

圖6 不同節理傾角下數值模型的應力-應變曲線Fig.6 Stress-strain curves of numerical models under different joint angles

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

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

為研究壓縮過程中不同節理傾角下微裂紋的發育情況，提取不同節理傾角下微裂紋數量隨應變的變化曲線，如圖8所示。結果發現不同節理傾角下微裂紋的變化曲線差異較大。當充填節理處于水平位置時數值模型的微裂紋發育主要呈現出前期增長較快、中期增長緩慢、后期迅速增長的3階段特點。對照圖4 對所呈現的微裂紋3 階段發育特點的分析如下：前期裂紋增長較快的主要原因是充填節理層強度及模量較低，在壓力作用下更易產生裂紋且低強度材料中裂紋擴展較快；中期裂紋增長緩慢的主要原因是充填節理層裂紋不再增加，裂紋產生位置逐步轉移至兩側巖石層；后期迅速增長的主要原因是兩側巖石層裂紋在壓力作用下不斷匯聚和擴展，充填節理試樣趨于失穩破壞。

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

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

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

圖10 起裂應力示意圖Fig.10 Diagram of initiation stress

圖11 不同節理傾角下起裂應力變化曲線Fig.11 Fracture initiation stress curve under different joint angles

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

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

2.2 不同充填節理厚度

圖13 為5 種不同充填厚度（2.0、3.5、5.0、6.5和8.0 mm）的充填節理巖石單軸壓縮的應力-應變曲線。可以看出，應力-應變曲線變化趨勢較有規律，隨著充填節理厚度的增加，充填節理巖石的應力-應變曲線逐漸向下方移動，其變化趨勢與試驗結果一致（Chai et al.，2020）。充填節理厚度分別為6.5 和8.0 mm 時應力-應變曲線峰前出現波動現象，這是由于充填節理厚度較大時在壓縮過程中充填層顆粒會發生錯動和向外擠壓情況。充填厚度改變時充填節理巖石的峰值應力和彈性模量的變化規律如圖14所示。

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

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

從圖14 可知當充填節理厚度不斷增加時，數值模型的峰值應力分別為57.82、 51.58、45.26、32.08 和23.83 MPa，彈性模量分別為1.31、1.103、0.851、0.809和0.754 GPa。充填厚度增加至3.5mm時峰值應力和彈性模量分別降低了10.79% 和15.8%；充填厚度增加至5.0 mm時峰值應力和彈性模量分別降低了21.72%和35.04%；充填厚度增加至6.5 mm 時峰值應力和彈性模量分別降低了44.52%和38.24%；充填厚度增加至8.0 mm 時峰值應力和彈性模量分別降低了58.79%和42.44%。可以發現峰值應力和彈性模量隨充填厚度的增加而逐漸降低，并且彈性模量降低速度逐漸變緩，進一步通過擬合發現峰值應力的降低趨勢符合線性函數分布，而彈性模量的降低趨勢符合多項式函數分布，相關系數均大于0.97，說明一致性良好。

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

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

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

如圖17 為不同厚度充填節理巖石數值模型的最終破壞形態。可以發現隨著充填節理厚度的增加，最終破壞時充填節理層顆粒受擠壓向兩側外溢、脫落情況嚴重，同時還發現兩側巖石層的微裂紋較少，微裂紋多集中于充填節理層，說明壓縮過程中充填節理層受損程度更嚴重，此外還可以看出不同厚度充填節理巖石數值模型最終破壞模式為張拉破壞，可以表明改變充填節理厚度對充填節理巖石壓縮破壞模式并無影響。

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

2.3 一組平行分布充填節理

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

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

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

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

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

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

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

圖23 為不同節理傾角下含1 組平行充填節理的巖石最終破壞形態。由圖23 可以看出，不同節理傾角下數值模型最終破壞形態差異較大，節理傾角為0°和15°時，充填節理層之間的巖石層均出現了貫通裂縫，節理傾角為30°時最終破壞形態完整性較好，當節理傾角大于45°后充填節理層出現分離剝落現象，尤其當節理傾角為75°時分離剝落現象最為嚴重，充填節理層均發生了明顯的滑移破壞，可知當節理傾角較大且巖石內部含有多條充填節理時整體穩定性較差、破壞程度嚴重。

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

3 結 論

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

1）隨著充填節理傾角的改變，峰值應力和彈性模量呈現出U 型分布，最終破壞時的總裂紋數和起裂應力呈現出先減小后增大的變化趨勢，此外改變充填節理傾角會嚴重影響充填節理巖石裂紋的發育狀況和最終破壞形態。

2）充填節理巖石的峰值應力和彈性模量隨充填厚度的增加而逐漸降低，且彈性模量降低速度逐漸變緩，而最終破壞時的總裂紋數隨著充填厚度的增加而不斷增加，進一步擬合發現峰值應力和總裂紋數的變化規律符合線性函數分布，彈性模量的變化規律符合多項式函數分布。

3）不同充填厚度的節理巖石微裂紋的發育過程均呈現典型的3階段特點，即微裂紋前期增長較快、中期增長緩慢、后期迅速增長；但充填節理厚度改變對起裂應力和破壞模式影響不大。

4）巖石內部存在1組平行充填節理時，隨著節理傾角的增加，峰值應力和彈性模量均出現先減小后增大的趨勢，最終破壞時總裂紋數也呈現出先減小后增大的趨勢。但當節理傾角較大會發生明顯的滑移破壞，充填節理巖石整體穩定性較差、破壞程度嚴重。

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