馬萬里，陳世江

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

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

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

以上研究豐富和發(fā)展了巖石非均性表征方法和巖石破裂過程中裂隙演化、能量釋放特征，巖石的非均質(zhì)性同樣對水力壓裂、爆破等有著重大影響，建立了基于MC準(zhǔn)則的應(yīng)變軟化本構(gòu)模型，通過FLAC內(nèi)嵌的FISH語言及Matlab編程，實(shí)現(xiàn)了巖石的非均質(zhì)賦參及彈性模量、內(nèi)聚力、內(nèi)摩擦角隨塑性應(yīng)變?nèi)趸接懥瞬煌|(zhì)度巖石破裂過程及不同裂隙處理方式的優(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分布，且無相關(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)變軟化過程。其內(nèi)置的應(yīng)變軟化模型是基于單元的塑性剪切應(yīng)變對單元的內(nèi)聚力、內(nèi)摩擦角、膨脹角的參數(shù)進(jìn)行弱化，對拉破壞單元的適用性較弱，因此，本文在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ù)值之和一半對應(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時表示巖石材料處于彈性階段，εp大于0時表示巖石進(jìn)入塑性階段，此時其內(nèi)聚力、內(nèi)摩擦角、彈性模量等參數(shù)將發(fā)生弱化。

若以內(nèi)聚力的弱化程度對損傷變量進(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)行對比。由于未進(jìn)行實(shí)驗(yàn)室試驗(yàn)，為了驗(yàn)證模型的可行性，本文以單參量軸向應(yīng)變(F)作為微元強(qiáng)度值衡量巖石損傷程度，并和楊圣奇等[21]試驗(yàn)數(shù)據(jù)進(jìn)行對比，對比結(jié)果及參數(shù)取值見圖1。由圖1可知，理論模型無法體現(xiàn)巖石壓縮過程中的初始壓密階段，彈性階段匹配性較高，同時屈服破壞階段應(yīng)力跌落貼合度較高。巖石單軸壓縮表現(xiàn)出較強(qiáng)的彈脆性，均值度m取值較大，為了理論曲線和實(shí)驗(yàn)室曲線較為接近，同樣可能導(dǎo)致模型參數(shù)的失真，但總體上驗(yàn)證了模型的可行性和正確性，表明本文建立模型具有較強(qiáng)的適用性，可以反映巖石全應(yīng)力-應(yīng)變整個過程。

圖1 實(shí)測曲線-理論曲線對比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ì)性對巖石破裂特征影響及破裂特征描述，本小節(jié)建立數(shù)值計(jì)算模型尺寸為100 mm×50 mm×50 mm立方體模型(圖2)，模型一共25 000個網(wǎng)格，上下表明施加位移速度邊界5e-7m/step。初始參數(shù)見表1，同時計(jì)算過程中對彈性單元數(shù)量、能量進(jìn)行監(jiān)測。

圖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)變軟化求解過程

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

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

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

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

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

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

圖5 巖石破裂過程破裂點(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ì)材料無法體現(xiàn)的。同時隨著均質(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對應(yīng)的抗壓強(qiáng)度分別為9.43 MPa、11.27 MPa、12.21 MPa，和均質(zhì)度呈對數(shù)函數(shù)關(guān)系。

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

隨著均質(zhì)度增大，聲發(fā)射事件數(shù)分布形態(tài)向倒V型轉(zhuǎn)變(圖9)，即峰前聲發(fā)射事件數(shù)逐漸增大，到達(dá)峰值時突然大幅度提高；隨著均質(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)，由之帶來的問題則是無法顯示巖石的破裂程度，如果僅從塑性區(qū)對某關(guān)鍵層位巖層穩(wěn)定性判定將會帶來誤差，并且相對保守。

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

(6)

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

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

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

4 結(jié) 論

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

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

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