非均質(zhì)巖石單軸壓縮細(xì)觀模擬及破裂特征描述

馬萬(wàn)里，陳世江

(內(nèi)蒙古科技大學(xué)礦業(yè)研究院，內(nèi)蒙古 包頭 014010)

巖石的非均質(zhì)性向來(lái)都是國(guó)內(nèi)外學(xué)者研究的重點(diǎn)，不同學(xué)者采用不同數(shù)值計(jì)算軟件對(duì)非均質(zhì)巖石的破裂特征進(jìn)行研究[1-2]。巖石的破裂特征對(duì)沖擊地壓、頂板冒落、巷道支護(hù)具有重要意義。

為了描述巖石的非均質(zhì)性，不同學(xué)者采用不同的方式進(jìn)行刻畫(huà)，唐春安等[3]通過(guò)假設(shè)巖石參數(shù)服從Weibul方式對(duì)巖石彈性模量、抗壓抗拉強(qiáng)度進(jìn)行賦值；田志等[4]通過(guò)數(shù)字巖石物理的形式對(duì)地層沉積過(guò)程中的非均質(zhì)性進(jìn)行表征；羅榮等[5]通過(guò)蒙特卡洛法對(duì)巖石非均質(zhì)性進(jìn)行描述；陳沙等[6]以圖像處理為媒介，利用Matlab灰度識(shí)別獲取巖石的截面特征，實(shí)現(xiàn)了巖石的非均質(zhì)表征。隨著計(jì)算機(jī)分析手段的進(jìn)化，CT掃描三維重構(gòu)技術(shù)得到了極大的發(fā)展，不同研究方向?qū)W者以CT掃描技術(shù)為基礎(chǔ)，研究了巖石的細(xì)觀破壞、滲流、結(jié)構(gòu)變化等，獲得了較多有益成果，極大發(fā)展了巖石力學(xué)[7-9]。巖石破裂方面：謝和平等[10]以能量耗散角度對(duì)巖石破裂過(guò)程不同階段進(jìn)行了精細(xì)刻畫(huà)；程昊等[11]采用C++對(duì)FLAC進(jìn)行二次開(kāi)發(fā)，實(shí)現(xiàn)了巖石破裂過(guò)程中聲發(fā)射參量導(dǎo)出；朱澤奇等[12]分析了花崗巖破裂及聲發(fā)射特性，明確了非均質(zhì)對(duì)巖石局部應(yīng)變產(chǎn)生重要影響；霍丙杰等[13]以FLAC數(shù)值計(jì)算軟件分析了房式采空區(qū)煤柱塑性區(qū)分布規(guī)律，明確了煤柱拉剪復(fù)合破壞機(jī)制。為了更好地描述巖石每個(gè)單元的破碎程度，基于單元損傷、單元安全度、破壞接近度等描述巖石應(yīng)力狀態(tài)及三維重構(gòu)裂隙越來(lái)越來(lái)受到人們的重視[14-15]，朱萬(wàn)成等[16]以彈性損傷建立了非均質(zhì)巖石應(yīng)力-損傷本構(gòu)模型；王峰[17]以MC準(zhǔn)則為基礎(chǔ)，定義了以單元應(yīng)力狀態(tài)點(diǎn)距包絡(luò)線距離為基礎(chǔ)的破壞接近度概念。

以上研究豐富和發(fā)展了巖石非均性表征方法和巖石破裂過(guò)程中裂隙演化、能量釋放特征，巖石的非均質(zhì)性同樣對(duì)水力壓裂、爆破等有著重大影響，建立了基于MC準(zhǔn)則的應(yīng)變軟化本構(gòu)模型，通過(guò)FLAC內(nèi)嵌的FISH語(yǔ)言及Matlab編程，實(shí)現(xiàn)了巖石的非均質(zhì)賦參及彈性模量、內(nèi)聚力、內(nèi)摩擦角隨塑性應(yīng)變?nèi)趸接懥瞬煌|(zhì)度巖石破裂過(guò)程及不同裂隙處理方式的優(yōu)缺點(diǎn)，為今后非均質(zhì)爆破試驗(yàn)及數(shù)值計(jì)算提供基礎(chǔ)性工作。

1 非均質(zhì)巖石應(yīng)變軟化本構(gòu)模型建立及驗(yàn)證

1.1 巖石非均質(zhì)性的描述

煤巖體物理力學(xué)性質(zhì)和顆粒分布、缺陷、節(jié)理等分布表現(xiàn)出較強(qiáng)的差異性，局部應(yīng)力、應(yīng)變導(dǎo)致不同位置應(yīng)力分布不同，假設(shè)其彈性模量、內(nèi)聚力、內(nèi)摩擦角、抗拉強(qiáng)度等參數(shù)服從weibul分布，且無(wú)相關(guān)性，可表示為式(1)[18]。

(1)

式中：ε*為微元強(qiáng)度；m為均質(zhì)度系數(shù)，表征巖石的非均質(zhì)程度，其值越大，表明巖石均質(zhì)性較高；M為某一力學(xué)參數(shù)的均值，表征了巖石某一力學(xué)性能強(qiáng)度。

1.2 基于MC準(zhǔn)則的應(yīng)變軟化本構(gòu)模型

FLAC3D內(nèi)置MC本構(gòu)模型為理想彈塑性模型，即進(jìn)入塑性階段后應(yīng)力隨著應(yīng)變的增加保持不變，很難描述巖石的應(yīng)變軟化過(guò)程。其內(nèi)置的應(yīng)變軟化模型是基于單元的塑性剪切應(yīng)變對(duì)單元的內(nèi)聚力、內(nèi)摩擦角、膨脹角的參數(shù)進(jìn)行弱化，對(duì)拉破壞單元的適用性較弱，因此，本文在FLAC內(nèi)置應(yīng)變軟化模型基礎(chǔ)上進(jìn)行修正，建立以塑性應(yīng)變?yōu)閰?shù)弱化依據(jù)參量，則其內(nèi)聚力、內(nèi)摩擦角、彈性模量與塑性應(yīng)變關(guān)系可表示為式(2)。

(2)

式中：c0、φ0、E0和cr、φr、Er分別為初始和殘余內(nèi)聚力、內(nèi)摩擦角、彈性模量；εp*為初始和殘余彈性模量數(shù)值之和一半對(duì)應(yīng)的塑性應(yīng)變。

因此，修正后的MC準(zhǔn)則為式(3)[19]。

(3)

式中：Nφ=(1+sinφ)/(1-sinφ)；φ為摩擦角；c為內(nèi)聚力；σt為抗拉強(qiáng)度；c(εp)為隨塑性應(yīng)變變化的內(nèi)聚力，εp為0時(shí)表示巖石材料處于彈性階段，εp大于0時(shí)表示巖石進(jìn)入塑性階段，此時(shí)其內(nèi)聚力、內(nèi)摩擦角、彈性模量等參數(shù)將發(fā)生弱化。

若以內(nèi)聚力的弱化程度對(duì)損傷變量進(jìn)行描述，損傷變量定義為式(4)。

(4)

根據(jù)連續(xù)介質(zhì)損傷力學(xué)理論，其本構(gòu)關(guān)系可以表示為式(5)[20]。

σ=Eε(1-D)

(5)

結(jié)合式(1)～式(5)得到基于MC準(zhǔn)則的非均質(zhì)巖石應(yīng)變軟化(應(yīng)力-損傷)本構(gòu)模型。

1.3 模型驗(yàn)證

本文建立的應(yīng)變軟化(應(yīng)力-損傷)模型為多參數(shù)弱化本構(gòu)模型，和FLAC內(nèi)置應(yīng)變軟化模型具有一定的相似性，為了驗(yàn)證模型的合理性及正確性，通常需要實(shí)驗(yàn)室試驗(yàn)反演獲取模型各參數(shù)數(shù)值，然后將實(shí)驗(yàn)室和FLAC數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比。由于未進(jìn)行實(shí)驗(yàn)室試驗(yàn)，為了驗(yàn)證模型的可行性，本文以單參量軸向應(yīng)變(F)作為微元強(qiáng)度值衡量巖石損傷程度，并和楊圣奇等[21]試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，對(duì)比結(jié)果及參數(shù)取值見(jiàn)圖1。由圖1可知，理論模型無(wú)法體現(xiàn)巖石壓縮過(guò)程中的初始?jí)好茈A段，彈性階段匹配性較高，同時(shí)屈服破壞階段應(yīng)力跌落貼合度較高。巖石單軸壓縮表現(xiàn)出較強(qiáng)的彈脆性，均值度m取值較大，為了理論曲線和實(shí)驗(yàn)室曲線較為接近，同樣可能導(dǎo)致模型參數(shù)的失真，但總體上驗(yàn)證了模型的可行性和正確性，表明本文建立模型具有較強(qiáng)的適用性，可以反映巖石全應(yīng)力-應(yīng)變整個(gè)過(guò)程。

圖1 實(shí)測(cè)曲線-理論曲線對(duì)比Fig.1 Comparison of measured curve and theoretical curve

2 非均質(zhì)巖石破裂特征數(shù)值計(jì)算

2.1 數(shù)值計(jì)算模型建立及邊界條件

為了驗(yàn)證將上文建立的非均質(zhì)巖石應(yīng)變軟化本構(gòu)模型引入FLAC，分析巖石非均質(zhì)性對(duì)巖石破裂特征影響及破裂特征描述，本小節(jié)建立數(shù)值計(jì)算模型尺寸為100 mm×50 mm×50 mm立方體模型(圖2)，模型一共25 000個(gè)網(wǎng)格，上下表明施加位移速度邊界5e-7m/step。初始參數(shù)見(jiàn)表1，同時(shí)計(jì)算過(guò)程中對(duì)彈性單元數(shù)量、能量進(jìn)行監(jiān)測(cè)。

圖2 數(shù)值計(jì)算模型Fig.2 Numerical calculation model

表1 單軸抗壓有限差分法模擬數(shù)據(jù)Table 1 Simulation data of uniaxial compressionfinite difference method

2.2 非均質(zhì)性及應(yīng)變軟化求解過(guò)程

為了表征巖石的非均質(zhì)性，本文通過(guò)Matlab編程生成weibul函數(shù)txt文件，數(shù)據(jù)點(diǎn)和單元網(wǎng)格數(shù)量保持一致，通過(guò)FLAC讀取文件的形式實(shí)現(xiàn)對(duì)不同參數(shù)的對(duì)應(yīng)賦參，同時(shí)由于本文不考慮數(shù)據(jù)點(diǎn)之間的關(guān)聯(lián)性，因此采用多個(gè)weibul數(shù)據(jù)文檔進(jìn)行賦值。具體FISH編程思路為：將模型每個(gè)單元力學(xué)參數(shù)值和weibul函數(shù)txt文件數(shù)值進(jìn)行一一對(duì)應(yīng)，同時(shí)利用parse功能實(shí)現(xiàn)數(shù)據(jù)的分隔，并采用額外變量顯示，均質(zhì)度為7的內(nèi)聚力和內(nèi)摩擦角分布特征見(jiàn)圖3。

圖3 巖石材料非均質(zhì)描述Fig.3 Description of heterogeneity of rock material

建立FISH函數(shù)對(duì)單軸壓縮過(guò)程中單元內(nèi)聚力、內(nèi)摩擦角、彈性模量等參數(shù)與塑性應(yīng)變關(guān)系進(jìn)行描述，計(jì)算一步求解每個(gè)單元的塑性應(yīng)變，并根據(jù)式(2)對(duì)模型各參數(shù)進(jìn)行修改，并監(jiān)測(cè)單元的彈塑性狀態(tài)、數(shù)量及彈性能，再計(jì)算一步并不斷循環(huán)，最終實(shí)現(xiàn)非均質(zhì)巖石單軸壓縮破裂全過(guò)程數(shù)值計(jì)算，具體求解過(guò)程見(jiàn)圖4。

圖4 數(shù)值計(jì)算求解流程Fig.4 Numerical calculation flow

2.3 單軸壓縮全應(yīng)力-應(yīng)變各階段

以均質(zhì)度m=7為例分析單軸壓縮過(guò)程中能量及裂隙演化全過(guò)程(圖5)，單軸壓縮初期，巖石內(nèi)部缺陷較大的點(diǎn)即力學(xué)參數(shù)較弱的點(diǎn)發(fā)生破裂，但整體數(shù)量較小，呈隨機(jī)分布，隨著加載步數(shù)增大，破裂點(diǎn)數(shù)量不斷增大，巖石仍處于彈性階段，此時(shí)，破裂點(diǎn)開(kāi)始出現(xiàn)匯集、融合等現(xiàn)象，這些為宏觀裂隙的出現(xiàn)提供了必要條件。巖石進(jìn)入屈服階段后，開(kāi)始出現(xiàn)大量拉伸破裂點(diǎn)，一旦進(jìn)入破壞階段，大量單元狀態(tài)由正在剪切破壞狀態(tài)向過(guò)去剪切狀態(tài)轉(zhuǎn)變，破裂點(diǎn)從隨機(jī)分布的無(wú)序狀態(tài)向局部剪切帶形成的有序狀態(tài)演變，這和混沌理論較為一致，最終形成“X”型剪切帶。巖石彈性能演化和整個(gè)應(yīng)力應(yīng)變過(guò)程貼合度較高(圖6)，且在應(yīng)力峰值時(shí)到達(dá)最大值，能量演化曲線呈現(xiàn)凹型演化和實(shí)驗(yàn)室應(yīng)力應(yīng)變曲線較為類似。若以巖石破裂點(diǎn)作為聲發(fā)射參量，可以看出隨著加載步數(shù)的增加巖石的總破裂點(diǎn)基本恒定，且聲發(fā)射參量在巖石峰值時(shí)稍靠右側(cè)達(dá)到最大值，之后由于巖石大面積卸載，導(dǎo)致聲發(fā)射參量不斷減小，這個(gè)過(guò)程為巖石內(nèi)部顆粒重新裝配過(guò)程，本文稱之為平靜期(圖7)。

圖5 巖石破裂過(guò)程破裂點(diǎn)分布Fig.5 Distribution of fracture points in the process of rock fracture

圖6 應(yīng)力能量曲線Fig.6 Stress energy curve

圖7 聲發(fā)射曲線Fig.7 Acoustic emission curve

3 討 論

3.1 不同均質(zhì)度能量演化及聲發(fā)射參量

不同均質(zhì)度條件下巖石應(yīng)力應(yīng)變曲線如圖8所示。巖石的彈性模量、峰值強(qiáng)度和均質(zhì)度呈正相關(guān)，即巖石的材料參數(shù)分布越統(tǒng)一，高強(qiáng)度材料參數(shù)單元數(shù)量越多，巖石的宏觀強(qiáng)度越高，這是均質(zhì)材料無(wú)法體現(xiàn)的。同時(shí)隨著均質(zhì)度的增加巖石的脆性程度逐漸增加，表現(xiàn)為峰后階段應(yīng)力跌落速度越快，這和內(nèi)部顆粒的膠結(jié)程度、顆粒強(qiáng)度有較大關(guān)系。均質(zhì)

圖8 不同均質(zhì)度應(yīng)力應(yīng)變及彈模分布Fig.8 Stress strain and elastic modulus distribution of different homogeneity

度為3、5、7對(duì)應(yīng)的抗壓強(qiáng)度分別為9.43 MPa、11.27 MPa、12.21 MPa，和均質(zhì)度呈對(duì)數(shù)函數(shù)關(guān)系。

從不同均質(zhì)度單元破壞數(shù)量及聲發(fā)射參量可以看出，巖石試樣壓縮初始階段，均質(zhì)度較小的試件已有部分破裂單元，隨著均質(zhì)度的增大，破裂點(diǎn)出現(xiàn)位置對(duì)應(yīng)應(yīng)變逐漸增加，表明了高均質(zhì)材料具有較強(qiáng)的穩(wěn)定性。但從單元破壞總數(shù)量上來(lái)看，單元總破壞數(shù)量和均質(zhì)度并無(wú)直接關(guān)聯(lián)。不同均質(zhì)度對(duì)應(yīng)的聲發(fā)射特征差異性較大，均質(zhì)度較小時(shí)，聲發(fā)射事件數(shù)基本呈n型分布，僅在屈服階段出現(xiàn)小幅度提升；

隨著均質(zhì)度增大，聲發(fā)射事件數(shù)分布形態(tài)向倒V型轉(zhuǎn)變(圖9)，即峰前聲發(fā)射事件數(shù)逐漸增大，到達(dá)峰值時(shí)突然大幅度提高；隨著均質(zhì)度增大聲發(fā)射頻率及強(qiáng)度表現(xiàn)出不同的形式，這和前人總結(jié)的聲發(fā)射特征具有高度的一致性，即群震型、前震-主震-余震型和主震型3種典型模式[22]。

圖9 不同均質(zhì)度單元破壞及聲發(fā)射參量Fig.9 Failure and AE parameters of different homogeneity units

3.2 巖石破裂形態(tài)描述

為了更好地表征巖石破裂形態(tài)，不同學(xué)者以不同的思路進(jìn)行處理，通常我們采用FLAC自帶的塑性區(qū)顯示方式，認(rèn)為巖石進(jìn)入塑性區(qū)即進(jìn)入破壞狀態(tài)，由之帶來(lái)的問(wèn)題則是無(wú)法顯示巖石的破裂程度，如果僅從塑性區(qū)對(duì)某關(guān)鍵層位巖層穩(wěn)定性判定將會(huì)帶來(lái)誤差，并且相對(duì)保守。

本文則通過(guò)對(duì)比以塑性區(qū)、破壞接近度、損傷值(本文定義)等方式對(duì)破裂點(diǎn)的描述，分析各種形式的優(yōu)缺點(diǎn)，以期獲得更合理的表征形式。其中，關(guān)于塑性區(qū)的定義本文不再贅述，破壞接近度是空間應(yīng)力的一點(diǎn)沿最不利應(yīng)力路徑到屈服面的距離與相應(yīng)的最穩(wěn)定參考點(diǎn)在相同羅德角方向上沿最不利應(yīng)力路徑到屈服面的距離之比，見(jiàn)式(6)。

(6)

式中：I1、J2為應(yīng)力張量的第一不變量、第二不變量；θσ為羅德角。

本文同樣定義了損傷值表達(dá)式，已在前文進(jìn)行說(shuō)明，不同形式描述巖石破裂特征見(jiàn)圖10。由圖10可知，三者都可以反映出巖石的破裂特征，其中，塑性區(qū)表征程度較為一般，它僅反映了巖石破裂點(diǎn)分布，而損傷值則將剪切帶刻畫(huà)的更為細(xì)致，損傷值為1時(shí)表示巖石已完全破壞，不同損傷值表明了巖石所處的應(yīng)力狀態(tài)；損傷值為0時(shí)表示巖石處于彈性狀態(tài)。破壞接近度則更為細(xì)致，它將單元的應(yīng)力狀態(tài)統(tǒng)一起來(lái)進(jìn)行體現(xiàn)，從彈性、塑性、破壞狀態(tài)對(duì)每個(gè)單元的應(yīng)力狀態(tài)進(jìn)行刻畫(huà)，破壞接近度小于1時(shí)表明巖石處于彈性狀態(tài)，但其彌補(bǔ)了損傷值的不足，即使單元處于彈性狀態(tài)，它仍可對(duì)較危險(xiǎn)單元進(jìn)行劃分；破壞接近度大于2時(shí)表示單元處于破裂狀態(tài)，可以根據(jù)其值判定巖石較為危險(xiǎn)位置。綜上所述，基于破壞接近度的單元接近度評(píng)判法則有著較強(qiáng)的優(yōu)越性。

圖10 不同破裂表征形式Fig.10 Different fracture characterization forms

4 結(jié) 論

1) 本文建立了非均質(zhì)巖石應(yīng)力-損傷本構(gòu)模型，并采用前人數(shù)據(jù)驗(yàn)證了模型的正確性。通過(guò)Matlab軟件和FISH語(yǔ)言實(shí)現(xiàn)了FLAC非均質(zhì)巖石單軸壓縮破裂數(shù)值計(jì)算。

2) 不同均質(zhì)度單軸壓縮數(shù)值計(jì)算表明：?jiǎn)屋S抗壓強(qiáng)度和彈性模量和均質(zhì)度呈對(duì)數(shù)函數(shù)關(guān)系，隨著均質(zhì)度系數(shù)的增大，巖石由彈塑性向彈脆性演變。

3) 巖石的破裂形態(tài)由雜亂無(wú)章向局部有序演變，對(duì)不同描述破裂特征方式進(jìn)行了討論，確定了單元破壞接近度的優(yōu)越性。