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

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

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

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

1 室內單軸實驗

1.1 巖石試樣制備及實驗設備

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

圖1 單軸壓縮實驗機Fig.1 Uniaxial compression testing machine

1.2 實驗結果分析

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

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

圖2 實驗用巖心試樣Fig.2 Core sample for experiment

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

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

2.1 離散元法基本原理

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

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

2.2 接觸模型建立

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

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

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

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

2.3 單軸壓縮實驗模型建立

模擬巖樣如圖5所示，其中上下平板為鋼制材質，上側鋼板以特定的速度進行勻速運動模擬加壓板對模擬巖石樣本進行壓縮，用以模擬單軸壓縮實驗的壓力，其中下側模擬鋼板為固定狀態(tài)，用來保證模擬巖石在受到模擬壓力的作用下，下端被固定且無相對運動。

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

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

2.4 模擬實驗結果

模擬實驗結果如圖7-10所示。

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

圖6 黏結鍵分布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

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

圖9 壓縮力與時間關系Fig.9 Compression force time diagram

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

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

圖11 真實實驗與模擬實驗對比Fig.11 Comparison curves between real and simulated experiments

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

3 結 論

以顆粒離散元法為理論基礎，采用Hertz?Mindlin黏結接觸模型作為被壓縮材料，通過對仿真物理模型的影響因素的分析，并與室內實驗的結果進行對比，證明了模擬結果的可靠性。

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

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

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