基于離散元的巖石單軸壓縮實(shí)驗(yàn)?zāi)M及斷裂特性分析

李 瑋，孟凡海，凌 鑫，李思琪，陳子賀

（東北石油大學(xué)石油工程學(xué)院，黑龍江大慶163318）

單軸抗壓強(qiáng)度作為巖石的抗壓強(qiáng)度的重要物理學(xué)指標(biāo)，是研究巖石力學(xué)性質(zhì)的基礎(chǔ)。在巖石力學(xué)及相關(guān)的學(xué)科中被廣泛應(yīng)用。單軸壓縮變形破壞是一個(gè)復(fù)雜的過程，得到諸多研究者的重視。巖石單軸壓縮的應(yīng)力應(yīng)變過程一般分為4個(gè)階段,即巖體內(nèi)部裂隙壓密、彈性形變、塑性形變以及破壞[1]。應(yīng)力應(yīng)變曲線與裂紋數(shù)應(yīng)變曲線表明,巖石在壓縮過程中應(yīng)力與裂紋增長(zhǎng)具有一定的對(duì)應(yīng)關(guān)系，即在應(yīng)力峰值的30%~80%時(shí)，裂紋穩(wěn)定增長(zhǎng)；當(dāng)應(yīng)力達(dá)到峰值的70%~80%以后，裂紋開始加速增長(zhǎng)[2]。在巖石的破壞形式方面，拉張破壞為主的觀點(diǎn)得到了廣泛的認(rèn)同[3?6]。離散元應(yīng)用于巖石力學(xué)可將巖石顆粒間的連接看作是靠力鏈連接而成，在受力的狀態(tài)下力鏈因斷裂而發(fā)生彎曲，方向由柱狀逐漸變成環(huán)狀[7]。胡光輝等[8]考慮巖石的時(shí)效變形損傷過程，引入巖石細(xì)觀單元時(shí)效損傷的應(yīng)力腐蝕模型，建立了基于離散元方法的巖石時(shí)效變形損傷破裂模型，并通過單軸壓縮及單軸蠕變的室內(nèi)實(shí)驗(yàn)和數(shù)值模擬對(duì)比驗(yàn)證了所建立的時(shí)效變形損傷破裂模型的合理性。

隨著計(jì)算機(jī)的快速發(fā)展，應(yīng)用數(shù)值模擬對(duì)巖石的力學(xué)實(shí)驗(yàn)進(jìn)行仿真模擬變得越來(lái)越普遍，離散元軟件的發(fā)展使通過離散元方法模擬巖石的力學(xué)實(shí)驗(yàn)變得越來(lái)越簡(jiǎn)便。離散單元法能夠從宏微觀層面上對(duì)數(shù)值模擬結(jié)果進(jìn)行分析，因此，該方法在巖石的數(shù)值模擬中也得到迅速的推廣和應(yīng)用[9]。本文借助離散元軟件，采用Hertz?Mindlin黏結(jié)接觸模型模擬大慶三肇地區(qū)深層巖石在受到單軸壓力的作用下的斷裂特性，分析其受力性質(zhì)，并結(jié)合離散元理論分析基于力鏈的顆粒體系接觸應(yīng)力。與單軸壓縮室內(nèi)真實(shí)實(shí)驗(yàn)對(duì)比，通過離散元數(shù)值模擬的優(yōu)勢(shì)在于可以在微觀的角度觀察到巖石的破壞特性。

1 室內(nèi)單軸實(shí)驗(yàn)

1.1 巖石試樣制備及實(shí)驗(yàn)設(shè)備

本文采用20 t級(jí)單軸實(shí)驗(yàn)機(jī)進(jìn)行單軸壓縮實(shí)驗(yàn)，如圖1所示。實(shí)驗(yàn)所用的巖樣均取自大慶三肇地區(qū)深層巖石，并制備成直徑20 mm、高度40 mm的圓柱形試件，加工好的巖心試件如圖2所示。

圖1 單軸壓縮實(shí)驗(yàn)機(jī)Fig.1 Uniaxial compression testing machine

1.2 實(shí)驗(yàn)結(jié)果分析

實(shí)驗(yàn)前將巖石兩端面涂抹適量的黃油。選取加載速度為0.2 mm/s。選取平均100 s的加載時(shí)間，進(jìn)行12次實(shí)驗(yàn)，為獲取單軸抗壓強(qiáng)度與巖石基本參數(shù)，選取實(shí)驗(yàn)結(jié)果較合理的3、8、10、12的4組試樣進(jìn)行研究，結(jié)果如圖3所示。

由圖3可以看出，在達(dá)到壓縮力峰值前表現(xiàn)為彈性壓縮,在峰值附近有小部分的回彈現(xiàn)象，樣品應(yīng)力在峰值后迅速下降，塑形變形不明顯，表明該泥質(zhì)砂巖脆性較強(qiáng)。當(dāng)試件抗壓強(qiáng)度達(dá)到極限以后，軸向力?軸向變形曲線有一小段的下降緊接著出現(xiàn)一個(gè)峰值。這種現(xiàn)象一般是因?yàn)榘l(fā)生弱剪切面破壞，在接近試件上端面位置發(fā)生輕微破碎，并且直接造成試樣整體破壞，所以緊接著上壓板接觸到殘余的試樣，從而再次出現(xiàn)一個(gè)峰值。

圖2 實(shí)驗(yàn)用巖心試樣Fig.2 Core sample for experiment

圖3 壓縮力-軸向變形曲線Fig.3 Compression force axial deformation curve

2 離散元的基本原理與力學(xué)模型建立

2.1 離散元法基本原理

離散元法的基本原理是：將研究對(duì)象劃分為相互獨(dú)立的一個(gè)個(gè)單元，根據(jù)單元之間的相互作用以及牛頓運(yùn)動(dòng)定律，采用靜態(tài)松弛法或動(dòng)態(tài)松弛法等迭代方法進(jìn)行循環(huán)迭代計(jì)算，確定在每一個(gè)時(shí)間步長(zhǎng)所有單元的受力及位移，并更新所有單元的位置[10]。通過對(duì)每個(gè)單元的微觀運(yùn)動(dòng)對(duì)離散單元進(jìn)行跟蹤計(jì)算，從而得到整個(gè)研究對(duì)象在宏觀方面的運(yùn)動(dòng)規(guī)律。根據(jù)離散元基本原理的定義，單元間相互作用的規(guī)律通常被當(dāng)作是瞬態(tài)平衡的問題，定義只要對(duì)象內(nèi)部在多種作用力的條件下達(dá)到平衡狀態(tài)，就認(rèn)為其處于平衡狀態(tài)。離散元法的離散單元從幾何形狀上可分為顆粒元和塊體元兩大類。塊體元中以四面體元、六面體元為主要的離散單元，對(duì)于二維的離散元單體問題可以是任意多邊形元，但應(yīng)用的范圍比較單一，而且每個(gè)離散單元的單體都具有一個(gè)基本的節(jié)點(diǎn)(通常取形心點(diǎn))。在模擬過程中，主要采用球體元作為離散單元，對(duì)于二維問題的離散元問題一般采用圓盤形單元作為離散單元。

巖石作為非均質(zhì)物質(zhì)，在受力過程中具有復(fù)雜的變化，有限元通過畫網(wǎng)格的方式模擬巖石，在模擬巖石的受力破壞過程，不能準(zhǔn)確反映巖石的性質(zhì)。通過離散元的方法，將巖石離散成一個(gè)個(gè)顆粒，利用有限尺度的黏合劑黏結(jié)起來(lái)，這種黏結(jié)與實(shí)際巖石結(jié)構(gòu)相似，通過黏結(jié)阻止顆粒的切向和法向相對(duì)運(yùn)動(dòng)，當(dāng)達(dá)到最大法向和切向應(yīng)力時(shí)，這種黏結(jié)就會(huì)發(fā)生破壞，在破壞后顆粒作為剛性球體彼此相互作用，這種黏結(jié)模型與混凝土及巖石的結(jié)構(gòu)相似，特別適合模擬巖石結(jié)構(gòu)。

2.2 接觸模型建立

為了準(zhǔn)確反映巖石的性質(zhì)，應(yīng)用Hertz?Mindlin黏結(jié)接觸模型對(duì)巖石試樣進(jìn)行單軸壓縮實(shí)驗(yàn)與模擬仿真，通過與真實(shí)實(shí)驗(yàn)對(duì)比，對(duì)參數(shù)進(jìn)行校核。接觸模型的黏結(jié)約束產(chǎn)生的合力在兩球形相交界面的中心（見圖4），在特定的時(shí)刻顆粒被黏結(jié)起來(lái)，在此之前顆粒通過默認(rèn)的Hertz?Mindlin接觸模型產(chǎn)生相互作用[11]。

圖4 巖石顆粒接觸模型Fig.4 Contact model of rock particles

進(jìn)行單軸壓縮室內(nèi)實(shí)驗(yàn)材料物性參數(shù)見表1。根據(jù)模擬要求對(duì)材料間的接觸參數(shù)進(jìn)行設(shè)置，材料之間的靜摩擦因數(shù)0.25，滾動(dòng)摩擦因數(shù)0.80，恢復(fù)系數(shù)0.05；材料與鋼板之間靜摩擦因數(shù)0.50，滾動(dòng)摩擦因數(shù)0.50，恢復(fù)系數(shù)0.02。

表1 材料物性參數(shù)Table 1 Material parameters

2.3 單軸壓縮實(shí)驗(yàn)?zāi)Ｐ徒?/h3>

模擬巖樣如圖5所示，其中上下平板為鋼制材質(zhì)，上側(cè)鋼板以特定的速度進(jìn)行勻速運(yùn)動(dòng)模擬加壓板對(duì)模擬巖石樣本進(jìn)行壓縮，用以模擬單軸壓縮實(shí)驗(yàn)的壓力，其中下側(cè)模擬鋼板為固定狀態(tài)，用來(lái)保證模擬巖石在受到模擬壓力的作用下，下端被固定且無(wú)相對(duì)運(yùn)動(dòng)。

圖5 模擬巖樣示意Fig.5 Schematic diagram of simulated rock sample

選取4組較理想的實(shí)驗(yàn)結(jié)果取平均值，進(jìn)行虛擬實(shí)驗(yàn)標(biāo)定。在EDEM軟件中,隨機(jī)巖石顆粒半徑由R0與倍率Cf相乘獲得，隨機(jī)巖石顆粒接觸半徑由R0'與倍率Cf相乘獲得。巖石顆粒黏結(jié)參數(shù)如表2所示。顆粒開始進(jìn)行黏結(jié)的時(shí)間為0.12 s。仿真時(shí)，固定時(shí)間步設(shè)置為1×10-6，初始仿真總時(shí)間為0.40 s，可根據(jù)具體時(shí)間實(shí)時(shí)進(jìn)行調(diào)整。EDEM的顆粒工廠可以根據(jù)用戶需要進(jìn)行定義，將編譯的dll文件導(dǎo)入顆粒工廠中。生成的巖石單軸壓縮仿真模型（黏結(jié)鍵分布）如圖6所示，生成的顆粒數(shù)為4 880，顆粒黏結(jié)鍵數(shù)為16 000，平均每個(gè)顆粒與周圍3.4個(gè)顆粒生成黏結(jié)鍵，該模型黏結(jié)充分。

2.4 模擬實(shí)驗(yàn)結(jié)果

模擬實(shí)驗(yàn)結(jié)果如圖7-10所示。

表2 巖石顆粒黏結(jié)參數(shù)Table 2 Cohesive parameters of rock particles

圖6 黏結(jié)鍵分布Fig.6 Distribution of bond bond

圖7 模型受力破壞示意Fig.7 Schematic diagram of stress failure of model

圖8 模型受力力鏈斷裂示意Fig.8 Schematic diagram of stress chain fracture of model

通過顆粒流中的球狀模型，借助離散元軟件可以很好地模擬出巖石在單軸壓力作用下斷裂的過程，通過離散元軟件的后處理功能給不同速度的離散元顆粒按照速度屬性的不同進(jìn)行分類著色（見圖7、8），模擬單軸壓縮實(shí)驗(yàn)的模擬總時(shí)間為0.40 s，根據(jù)模擬結(jié)果繪制應(yīng)力時(shí)間曲線（見圖9），可以看出不同的變形特性下顆粒的受力具體值及在受到單軸壓力的作用下力的變化趨勢(shì)。

圖9 壓縮力與時(shí)間關(guān)系Fig.9 Compression force time diagram

圖10 黏結(jié)力鏈數(shù)目與時(shí)間關(guān)系Fig.10 Time graph of the number of cohesive force chains

巖石在承受單軸壓力過程中的第一個(gè)階段為巖石的裂縫壓密階段，真實(shí)單軸壓縮實(shí)驗(yàn)時(shí)僅通過實(shí)驗(yàn)結(jié)果曲線不能揭示真實(shí)的巖石顆粒壓密的變化情況，通過離散元軟件可以看出在0.05 s有少數(shù)幾個(gè)顆粒具有較大速度，即巖石有被破壞的趨勢(shì)。通過圖8（a）可以看出，此時(shí)僅有極少數(shù)力鏈有斷裂的現(xiàn)象而且力鏈的方向與軸向載荷的方向相同。第二階段為力鏈斷裂鍵的增長(zhǎng)階段，從圖8(b)可以看出，此時(shí)較多顆粒有較大的速度，較多力鏈發(fā)生斷裂，此時(shí)模擬巖石試樣即將發(fā)生斷裂，但未達(dá)到巖石屈服極限，因此應(yīng)力繼續(xù)增大。此時(shí)的應(yīng)力達(dá)到峰值的81%左右，在此階段斷裂力鏈的數(shù)目增長(zhǎng)較快，但增長(zhǎng)速度比較平穩(wěn)，應(yīng)力時(shí)間曲線斜率近似為固定值，說(shuō)明模擬巖石試樣脆性較高。第三階段為巖石試樣完全開裂階段，從圖8（c）可以看出，此時(shí)巖石將要斷裂，在斷裂位置有較多的顆粒具有較大的速度，斷裂位置的力鏈多數(shù)發(fā)生斷裂，應(yīng)力達(dá)到極限，應(yīng)力時(shí)間曲線斜率近似為直線，巖石試樣完全被破壞。第四階段為到達(dá)峰值后應(yīng)力急速下降過程，從圖8（d）可以看出，此時(shí)模擬巖石試樣完全斷裂，兩部分?jǐn)嗔盐恢镁哂邢鄬?duì)速度，斷裂力鏈數(shù)目不再增加（見圖10），宏觀上能直接觀察到試樣的破壞狀態(tài)?？梢耘c單軸壓縮室內(nèi)實(shí)驗(yàn)對(duì)比，結(jié)果如圖11所示。

圖11 真實(shí)實(shí)驗(yàn)與模擬實(shí)驗(yàn)對(duì)比Fig.11 Comparison curves between real and simulated experiments

由圖11可以看出，首先在縮容段，室內(nèi)實(shí)驗(yàn)曲線與模擬曲線完全重合，此時(shí)主要是將巖石微小的裂縫被壓實(shí)的過程，因此，按照本文方法建立實(shí)驗(yàn)?zāi)Ｐ驮谳^小的顆粒狀態(tài)下可以較準(zhǔn)確模擬出巖石空隙情況，在彈性形變過程中，兩個(gè)曲線并沒有完全重合，但是兩種曲線近似平行，與巖石的力學(xué)性能符合，對(duì)比塑性變形過程可以看出兩條曲線具有微小的差異，但是差異在可接受的范圍，在最后的擴(kuò)容階段，室內(nèi)實(shí)驗(yàn)具有一個(gè)峰值，在模擬結(jié)果中未有此種現(xiàn)象發(fā)生，這是由于室內(nèi)實(shí)驗(yàn)最初在上端面發(fā)生破碎，當(dāng)最上面發(fā)生破碎時(shí)試樣壓縮力迅速下降，而在模擬過程中巖石在試樣中間發(fā)生斷裂，且壓縮力的變化過程較慢。

3 結(jié) 論

以顆粒離散元法為理論基礎(chǔ)，采用Hertz?Mindlin黏結(jié)接觸模型作為被壓縮材料，通過對(duì)仿真物理模型的影響因素的分析，并與室內(nèi)實(shí)驗(yàn)的結(jié)果進(jìn)行對(duì)比，證明了模擬結(jié)果的可靠性。

（1）通過單軸壓縮室內(nèi)實(shí)驗(yàn)得到巖石材料參數(shù)、接觸參數(shù)，通過虛擬參數(shù)標(biāo)定模擬出具有和大慶三肇地區(qū)深層巖石相似力學(xué)性質(zhì)的巖石。

（2）通過與單軸壓縮真實(shí)壓縮實(shí)驗(yàn)的對(duì)比，即與真實(shí)單軸壓縮實(shí)驗(yàn)應(yīng)力應(yīng)變的曲線進(jìn)行擬合，在完全破壞前也具有回彈現(xiàn)象發(fā)生，模擬實(shí)驗(yàn)與真實(shí)實(shí)驗(yàn)的擬合度較高，可以通過EDEM模擬單軸壓縮真實(shí)實(shí)驗(yàn)。

（3）通過模擬典型巖石脆性破壞過程，得力鏈在受力破裂之前沒有較大的變形過程，破壞形式為直接斷裂，也就是說(shuō)這些巖石在荷載作用下沒有顯著覺察地變形就突然破壞，產(chǎn)生這種破壞的原因可能是巖石中裂隙的發(fā)生和發(fā)展的結(jié)果。