剛度差異對巖石破裂影響的顆粒流模擬研究

陳軍濤，李 昊，褚衍玉，張 毅，唐道增，朱 君，李文昕，王開宇

(1. 山東科技大學(xué) 能源與礦業(yè)工程學(xué)院，山東 青島 266590 ； 2.煤炭資源高效開采與潔凈利用國家重點(diǎn)實(shí)驗(yàn)室，北京 100013； 3. 山東能源臨沂礦業(yè)集團(tuán)有限責(zé)任公司 邱集煤礦，山東 德州 251105； 4. 中國安全生產(chǎn)科學(xué)研究院，北京 100012； 5.山東科技大學(xué) 電氣與自動化工程學(xué)院，山東 青島 266590)

巖土工程的失穩(wěn)破壞往往與巖體內(nèi)部裂隙的張開、閉合、擴(kuò)展及貫通密切相關(guān)，研究巖石破裂過程及其規(guī)律對巖土工程具有重大意義[1-2]。實(shí)際工程中很難采用探測手段觀察到巖石的破裂過程[3]，顆粒流模擬則成為研究巖土類材料裂紋擴(kuò)展破裂的重要手段，而宏觀參數(shù)與細(xì)觀參數(shù)的轉(zhuǎn)換問題是模擬準(zhǔn)確性的關(guān)鍵[4]。國內(nèi)外學(xué)者對顆粒流模擬進(jìn)行了大量的研究，李守巨等[5]提出一種基于巴西劈裂試驗(yàn)的混凝土細(xì)觀本構(gòu)模型參數(shù)估計(jì)方法；尹小濤等[6]研究了黏結(jié)內(nèi)聚力和摩擦角對巖石各項(xiàng)宏觀力學(xué)特征的影響；袁康等[7]利用離散元軟件分析了加載過程中錨桿和結(jié)構(gòu)面之間的相互作用機(jī)理以及錨桿內(nèi)部軸向和剪切應(yīng)力情況；鄧樹新等[8]利用Plackett-Burman試驗(yàn)設(shè)計(jì)并分析了細(xì)觀參數(shù)對宏觀響應(yīng)的敏感性，建立了宏觀力學(xué)指標(biāo)和細(xì)觀參數(shù)之間的線性關(guān)系；Yoon[9]采用中心合成設(shè)計(jì)(central composite design, CCD)標(biāo)定了接觸模型的細(xì)觀參數(shù)，研究了單軸壓縮強(qiáng)度和細(xì)觀參數(shù)之間的關(guān)系；Coetzee等[10]通過組合的模擬剪切試驗(yàn)和壓縮試驗(yàn)研究了顆粒間摩擦系數(shù)和黏結(jié)剛度；Stavrou等[11]采用離散元方法進(jìn)行巖體的尺度分析，提出了塊石強(qiáng)度與其體積和條件的關(guān)系。以上研究大多揭示了平行黏結(jié)模型細(xì)觀參數(shù)(如黏結(jié)內(nèi)聚力、摩擦角、黏結(jié)強(qiáng)度)與宏觀參數(shù)之間的關(guān)系，而剛度參數(shù)與巖石的宏觀破裂密不可分[12]，目前鮮有這方面的研究。

圖1 線性平行黏結(jié)模型[9]Fig. 1 Linear parallel bond mode

1 試驗(yàn)方案

1.1 線性平行黏結(jié)模型

假設(shè)顆粒有效模量等于黏結(jié)有效模量，則法向剛度與切向剛度的比值等于線性平行黏結(jié)法向剛度與線性平行黏結(jié)切向剛度的比值。此時(shí)，模型彈性模量可表示為[16]：

E/Ec=a+bln(kn/ks)。

(1)

式中：E為彈性模量，GPa；Ec為有效模量，GPa；kn/ks為剛度比；a、b為常數(shù)。顆粒有效模量與法向剛度的關(guān)系如圖2所示，表達(dá)式為[17]：

圖2 線性平行黏結(jié)ball-ball接觸[9]Fig. 2 Linear parallel bond contact of ball-ball

emod=knL/A，

(2)

(3)

(4)

(5)

式中：R(1)為接觸兩端某一個(gè)顆粒的半徑；R(2)為接觸兩端另一個(gè)顆粒的半徑；L為接觸兩端顆粒的中心距；r為較小顆粒的半徑；A為接觸的截面積；t為接觸界面的厚度。

線性平行黏結(jié)有效模量與線性平行黏結(jié)法向剛度之間的關(guān)系[18]為：

(6)

綜上分析，當(dāng)模型剛度比確定時(shí)，顆粒、黏結(jié)有效模量可以控制顆粒和黏結(jié)的法向剛度和切向剛度，并且與宏觀彈性模量相關(guān)。眾所周知，巖石材料是典型的非均質(zhì)體，而視為均質(zhì)體的巖石試驗(yàn)?zāi)M獲得的結(jié)論與實(shí)際有一定差距。因此，本研究通過離散元模擬試驗(yàn)，研究細(xì)觀參數(shù)顆粒有效模量、黏結(jié)有效模量對巖石裂紋擴(kuò)展和破裂過程的影響。

圖3 模型示意圖Fig. 3 Model diagram

1.2 模擬方案設(shè)計(jì)

建立如圖3所示的離散元模型，尺寸為50 mm×100 mm，線性平行黏結(jié)模型。模型細(xì)觀參數(shù)顆粒有效模量、黏結(jié)有效模量為定值，孔隙率為0.1，顆粒半徑0.75×10-3m，密度2 500 kg/m3，顆粒局部阻尼系數(shù)0.7。通過墻體建立顆粒生成范圍，在范圍內(nèi)隨機(jī)生成顆粒。

現(xiàn)定義剛度差異系數(shù)Ka，反映模型內(nèi)部顆粒剛度和膠結(jié)體剛度的差異程度。

(7)

當(dāng)Ka>1時(shí)，顆粒剛度大于膠結(jié)物剛度，Ka越大，顆粒與膠結(jié)體的剛度差別越大；當(dāng)0

對模型賦予線性平行黏結(jié)后，查看接觸Contactpb_state數(shù)值，確認(rèn)接觸模型賦值成功。模擬試驗(yàn)采用位移方式加載，軸向加載速率為0.2 mm/s。通過fish函數(shù)記錄應(yīng)力-應(yīng)變變化過程、拉剪裂紋數(shù)目及裂紋擴(kuò)展過程，峰值應(yīng)力一半處視為彈性階段，求得試樣彈性模量，試驗(yàn)結(jié)束后保存模型并記錄破裂形態(tài)。

2 顆粒有效模量對巖石宏觀參數(shù)及破裂過程的影響特征

2.1 應(yīng)力-應(yīng)變曲線

圖4為黏結(jié)有效模量一定時(shí)，顆粒有效模量分別為1、3、5、7、9和11 GPa時(shí)的單軸壓縮全應(yīng)力-應(yīng)變曲線。可以看出，應(yīng)力-應(yīng)變曲線的前半部分基本為直線段，屬于彈性變形階段，峰前曲線明顯偏離直線，偏離直線的起點(diǎn)視為屈服點(diǎn)，屈服點(diǎn)到峰值之間為非穩(wěn)定破裂發(fā)展階段，峰后試樣表現(xiàn)出明顯的脆性。模擬曲線并未表征試驗(yàn)具有明顯的裂隙壓密階段，與室內(nèi)試驗(yàn)得到的裂隙壓密、彈性變形、微彈性裂隙穩(wěn)定發(fā)展、非穩(wěn)定破裂發(fā)展等階段有差異，這主要是因?yàn)檐浖茈y準(zhǔn)確模擬巖石復(fù)雜的原生孔隙和裂隙網(wǎng)絡(luò)，而默認(rèn)為致密模型，故在模擬結(jié)果中沒有明顯的裂隙孔隙壓密的現(xiàn)象。

圖4 不同顆粒有效模量的應(yīng)力-應(yīng)變曲線(單軸壓縮)Fig. 4 Stress-strain curves of effective modulus of different particles (uniaxial compression)

從圖4中還可以看出：巖石試樣破裂前，彈性模量和單軸抗壓強(qiáng)度隨顆粒有效模量的增大而增大；軸向應(yīng)變值隨顆粒有效模量的增大而減小，呈負(fù)相關(guān)關(guān)系。當(dāng)顆粒有效模量較大時(shí)，模型更容易出現(xiàn)脆性破壞；當(dāng)顆粒有效模量較小時(shí)，模型更容易出現(xiàn)塑性破壞。

PFC2D模擬試驗(yàn)中，顆粒剛度通過法向剛度和切向剛度控制，顆粒剛度越大，顆粒越難以發(fā)生彈性變形。由式(2)可知，顆粒有效模量與法向剛度正相關(guān)，在黏結(jié)有效模量和剛度比確定的情況下，顆粒有效模量越大，法向剛度和切向剛度越大，同樣應(yīng)力條件下模型變形量越小。

黏結(jié)有效模量一定時(shí)，顆粒有效模量為1、3、5、7、9和11 GPa時(shí)，試樣拉伸應(yīng)力峰值、彈性模量和應(yīng)變峰值幾乎沒有影響，僅顆粒有效模量為11 GPa的試樣在應(yīng)力峰值后發(fā)生了小幅波動。由式(1)可知，顆粒有效模量與顆粒剛度成正相關(guān)，拉伸應(yīng)力主要作用在顆粒與顆粒之間的接觸上。因此，在細(xì)觀參數(shù)黏結(jié)有效模量確定的情況下，拉伸條件下試樣各宏觀參數(shù)幾乎不受顆粒有效模量變化的影響。

2.2 顆粒有效模量對巖石宏觀參數(shù)的影響分析

根據(jù)2.1的研究結(jié)果，改變顆粒有效模量時(shí)，對巖石壓縮過程中彈性模量、峰值應(yīng)變、單軸抗壓強(qiáng)度等宏觀參數(shù)的變化進(jìn)行分析。單軸抗壓強(qiáng)度和應(yīng)變隨顆粒有效模量的增大呈非線性變化趨勢，當(dāng)顆粒有效模量不超過7 GPa時(shí)對模型的單軸抗壓強(qiáng)度影響較大，大于7 GPa時(shí)模型的單軸抗壓強(qiáng)度趨于穩(wěn)定。當(dāng)顆粒有效模量不超過5 GPa時(shí)對模型的應(yīng)變影響較大，大于5 GPa時(shí)模型的應(yīng)變趨于穩(wěn)定。反映顆粒有效模量增大到一定數(shù)值(本試驗(yàn)為7 GPa)后對模型單軸抗壓強(qiáng)度、應(yīng)變影響越來越小。

圖5～7分別為壓縮試驗(yàn)過程中，不同顆粒有效模量的巖石彈性模量、峰值應(yīng)變、單軸抗壓強(qiáng)度的變化曲線，可以看出隨著顆粒有效模量的增大，彈性模量和單軸抗壓強(qiáng)度均表現(xiàn)出上升趨勢，峰值應(yīng)變則表現(xiàn)出下降趨勢。對于擬合結(jié)果的R2，彈性模量(0.99)>峰值應(yīng)變(0.87)>單軸抗壓強(qiáng)度(0.84)，表明宏觀彈性模量與細(xì)觀顆粒有效模量呈較好的線性關(guān)系；對于擬合結(jié)果的斜率，單軸抗壓強(qiáng)度(5.20)>彈性模量(0.71)>峰值應(yīng)變(0.26)，說明顆粒有效模量發(fā)生改變時(shí)，巖石的單軸抗壓強(qiáng)度變化最大。

2.3 剛度差異系數(shù)對巖石破裂的影響分析

圖8為顆粒有效模量影響下不同剛度差異系數(shù)值時(shí)的巖石的壓縮破裂形態(tài)，圖9為顆粒有效模量影響下巖石剪切裂紋數(shù)量、拉伸裂紋數(shù)量、裂紋總數(shù)隨剛度差異系數(shù)的變化曲線，圖10為顆粒有效模量影響下拉伸、剪切裂紋占裂紋總數(shù)的比例隨剛度差異系數(shù)的變化曲線。

從圖8中可以看出，總體來說，剛度差異系數(shù)越大，巖石破裂產(chǎn)生的裂紋越集中；剛度差異系數(shù)越小，巖石破裂產(chǎn)生的裂紋越分散，說明剛度差異系數(shù)的大小影響著巖石破裂過程中的裂隙分布形態(tài)。從圖9中可以看出，當(dāng)剛度差異系數(shù)較小(≤7)時(shí)，拉伸裂紋數(shù)量和裂紋總數(shù)變化不穩(wěn)定，當(dāng)剛度差異系數(shù)較大(>7)時(shí)，拉伸裂紋數(shù)量和裂紋總數(shù)變化趨于穩(wěn)定，而剪切裂紋數(shù)量較少，在31～63范圍內(nèi)波動。圖10表明，隨剛度差異系數(shù)增大，拉伸裂紋和剪切裂紋占比僅在剛度差異系數(shù)為3時(shí)有較大波動，剛度差異系數(shù)為其他值時(shí)拉伸裂紋總數(shù)保持在80%左右，剪切裂紋總數(shù)保持在20%左右，說明剛度差異系數(shù)增大對巖石的剪切拉伸兩類裂紋占比幾乎沒有影響。

圖5 彈性模量Fig. 5 Elastic modulus

圖6 峰值應(yīng)變Fig. 6 Peak strain

圖7 單軸抗壓強(qiáng)度Fig. 7 Uniaxial compressive strength

圖8 壓縮破裂形態(tài)Fig. 8 Compression fracture morphology

圖9 裂紋數(shù)量Fig. 9 Number of cracks

圖10 拉剪裂紋占比Fig. 10 Proportion of tensile and shear cracks

3 黏結(jié)有效模量對巖石宏觀參數(shù)及其破裂過程的影響特征

3.1 應(yīng)力-應(yīng)變曲線

圖11為顆粒有效模量一定時(shí)，黏結(jié)有效模量分別為1、3、5、7、9和11 GPa時(shí)的單軸壓縮全應(yīng)力-應(yīng)變曲線。從圖11中可以看出：巖石破裂前，曲線斜率隨黏結(jié)有效模量的增大而增大，即彈性模量隨黏結(jié)有效模量的增大而增大；單軸抗壓強(qiáng)度和軸向應(yīng)變隨黏結(jié)有效模量的增大而減小。當(dāng)黏結(jié)有效模量較大時(shí)，模型更容易出現(xiàn)脆性破壞，反之更容易出現(xiàn)塑性破壞，但是影響程度較小。

由式(1)可知，黏結(jié)有效模量與模型的黏結(jié)法向剛度成正相關(guān)關(guān)系，在顆粒有效模量確定的情況下，增大黏結(jié)有效模量相當(dāng)于增大了接觸彈性體的剛度，同樣應(yīng)力條件下彈性體的剛度越小模型變形量越小。

圖12為黏結(jié)有效模量分別為1、3、5、7、9、11 GPa時(shí)的拉伸應(yīng)力-應(yīng)變曲線，可以看出隨著黏結(jié)有效模量增大，巖石的彈性模量增大，相同應(yīng)力水平下應(yīng)變量減小，峰值應(yīng)力幾乎保持不變。

圖11 不同黏結(jié)有效模量單軸壓縮應(yīng)力-應(yīng)變曲線Fig. 11 Stress-strain curves of uniaxial compression with different bonding effective modulus

圖12 不同黏結(jié)有效模量拉伸應(yīng)力-應(yīng)變曲線Fig. 12 Tensile stress-strain curves of different bonding effective modulus

3.2 黏結(jié)有效模量對巖石宏觀力學(xué)參數(shù)的影響分析

圖13～15為顆粒有效模量一定時(shí)，單軸壓縮試驗(yàn)過程中試樣壓彈性模量、壓縮峰值應(yīng)變和單軸抗壓強(qiáng)度隨黏結(jié)有效模量的變化趨勢。圖16～18為顆粒有效模量一定時(shí)，單軸拉伸試驗(yàn)過程中試樣拉彈性模量、拉伸峰值應(yīng)變、抗拉強(qiáng)度隨黏結(jié)有效模量的變化趨勢。

圖13 壓彈性模量Fig. 13 Elastic modulus under compression

圖14 壓縮峰值應(yīng)變Fig. 14 Compression peak strain

圖15 單軸抗壓強(qiáng)度Fig. 15 Uniaxial compressive strength

圖16 拉彈性模量Fig. 16 Tensile modulus of elasticity

圖17 拉伸峰值應(yīng)變Fig. 17 Tensile peak strain

圖18 抗拉強(qiáng)度Fig. 18 Tensile strength

單軸抗壓強(qiáng)度、拉應(yīng)變、壓應(yīng)變隨黏結(jié)有效模量的增大呈非線性變化趨勢，黏結(jié)有效模量變化對模型抗拉強(qiáng)度幾乎沒有影響。當(dāng)黏結(jié)有效模量不超過5 GPa時(shí)對模型的拉、壓應(yīng)變影響較大，當(dāng)黏結(jié)有效模量大于5 GPa時(shí)模型的拉、壓應(yīng)變趨于穩(wěn)定；當(dāng)黏結(jié)有效模量不超過3 GPa時(shí)對模型的單軸抗壓強(qiáng)度影響較大，當(dāng)黏結(jié)有效模量大于3 GPa時(shí)模型的單軸抗壓強(qiáng)度變趨于穩(wěn)定。反映黏結(jié)有效模量增大到一定數(shù)值(本試驗(yàn)為5 GPa)后對模型單軸抗壓強(qiáng)度、拉應(yīng)變、壓應(yīng)變影響越來越小。

通過觀察擬合結(jié)果的R2，拉彈性模量(1)>壓彈性模量(0.995)>單軸抗壓強(qiáng)度(0.75)>單軸抗拉強(qiáng)度(0.73)>壓縮峰值應(yīng)變(0.72)>拉伸峰值應(yīng)變(0.66)，發(fā)現(xiàn)宏觀參數(shù)拉、壓彈性模量與細(xì)觀參數(shù)黏結(jié)有效模量表現(xiàn)出良好的線性關(guān)系。從擬合結(jié)果斜率來看，拉彈性模量(1.61)>壓彈性模量(1.15)>單軸抗壓強(qiáng)度(0.97)>壓縮峰值應(yīng)變(0.33)>拉伸峰值應(yīng)變(0.22)>單軸抗拉強(qiáng)度(0.04)，說明黏結(jié)有效模量增大對拉、壓彈性模量的影響最大。

3.3 剛度差異系數(shù)對巖石破裂的影響分析

圖19為黏結(jié)有效模量變化影響下不同剛度差異系數(shù)值時(shí)試件的壓縮破裂形態(tài)，圖20為黏結(jié)有效模量變化影響下剪切裂紋數(shù)量、拉伸裂紋數(shù)量和裂紋總數(shù)隨剛度差異系數(shù)變化曲線，圖21為黏結(jié)有效模量變化影響下拉伸、剪切裂紋數(shù)量占比隨剛度差異系數(shù)變化曲線。結(jié)果表明，剛度差異系數(shù)較小時(shí)(≤1/5)，產(chǎn)生的裂紋數(shù)量較多，且裂紋數(shù)量變化較大；剛度差異系數(shù)較大時(shí)(>1/5)裂紋數(shù)目較穩(wěn)定。隨著剛度差異系數(shù)的減小，拉、剪裂紋的占比上下波動，剪切裂紋保持在16%～31%，拉伸裂紋占比保持在69%～84%，拉伸裂紋是剪切裂紋的2～4倍，說明剛度差異系數(shù)減小時(shí)對巖石的裂紋演化及破裂有一定影響，但總體影響不大。

圖19 壓縮破裂形態(tài)Fig. 19 Compression fracture morphology

圖20 裂紋數(shù)量變化Fig. 20 Change of crack number

圖21 拉剪裂紋占比Fig. 21 Proportion of tensile and shear cracks

4 實(shí)例驗(yàn)證

以文獻(xiàn)[19]模擬結(jié)果為例，參考上述本研究結(jié)果，對單軸壓縮試驗(yàn)進(jìn)行參數(shù)優(yōu)化研究。圖22(a)為文獻(xiàn)[19]的室內(nèi)試驗(yàn)結(jié)果，圖22(b)為文獻(xiàn)[19]的模擬試驗(yàn)結(jié)果(Ka=1)，可以發(fā)現(xiàn)試件破壞模式以剪切為主，形成了與試件軸向呈一定夾角的宏觀剪切裂紋，而模擬試驗(yàn)結(jié)果中裂隙分布散亂，巖石裂隙沒有相互貫通形成宏觀的剪切裂隙，與室內(nèi)試驗(yàn)獲得的結(jié)果相差較大。

圖22 巖石破裂特征Fig. 22 Rock fracture characteristics

由前面研究可知，當(dāng)顆粒有效模量或黏結(jié)有效模量增大時(shí)，試樣越容易出現(xiàn)剪切破壞，且顆粒有效模量影響更大；剛度差異系數(shù)值增大時(shí)，裂紋數(shù)目減小。因此，保持其他參數(shù)不變，調(diào)整顆粒有效模量，增大剛度差異系數(shù)值，使得巖石裂紋集中在實(shí)際破裂區(qū)域，最終出現(xiàn)明顯的宏觀剪切裂隙，破壞形態(tài)與室內(nèi)巖石破裂吻合性較好，調(diào)整前后的模擬試驗(yàn)結(jié)果見圖22，細(xì)觀力學(xué)參數(shù)見表1，驗(yàn)證了前面研究結(jié)論的準(zhǔn)確性。

表1 調(diào)整前后的細(xì)觀力學(xué)參數(shù)Tab. 1 Meso mechanical parameters before and after adjustment

5 結(jié)論

1) 模型黏結(jié)有效模量恒定時(shí)，顆粒有效模量越大，宏觀參數(shù)壓彈性模量和單軸抗壓強(qiáng)度越大，壓縮應(yīng)變反而減小；顆粒有效模量對拉伸條件下的宏觀參數(shù)幾乎沒有影響；當(dāng)顆粒有效模量較大時(shí)，模型易出現(xiàn)脆性破壞，反之易出現(xiàn)塑性破壞；顆粒有效模量與壓彈性模量具有較好的線性關(guān)系。

2) 模型顆粒有效模量恒定時(shí)，黏結(jié)有效模量越大，宏觀參數(shù)壓彈性模量和拉彈性模量越大，單軸抗壓抗拉強(qiáng)度、壓縮和拉伸應(yīng)變越?。划?dāng)黏結(jié)有效模量較大時(shí)，模型易出現(xiàn)脆性破壞，反之易出現(xiàn)塑性破壞；黏結(jié)有效模量與壓、拉彈性模量都具有良好的線性關(guān)系。模擬試驗(yàn)時(shí)，應(yīng)優(yōu)先標(biāo)定黏結(jié)有效模量，以提高模擬結(jié)果的準(zhǔn)確性。

3) 剛度差異系數(shù)反映了模型內(nèi)部顆粒和膠結(jié)體剛度的差異程度。剛度差異系數(shù)越大，模型越難產(chǎn)生裂隙，穩(wěn)定性好；剛度差異系數(shù)越小，模型越易產(chǎn)生裂隙，穩(wěn)定性差。

4) 以室內(nèi)巖石試驗(yàn)為例，對細(xì)觀參數(shù)進(jìn)行了優(yōu)化，參數(shù)優(yōu)化后的試驗(yàn)結(jié)果與室內(nèi)巖石破裂形態(tài)吻合較好，驗(yàn)證了本研究結(jié)論的準(zhǔn)確性。