基于摩爾庫倫準(zhǔn)則的數(shù)值流形法巖體破壞模擬分析

高文偉

（延安大學(xué)建筑工程學(xué)院，陜西 延安 716000）

隨著工程建設(shè)難度不斷增加，巖質(zhì)邊坡工程、地基工程、隧道工程中巖體破裂失穩(wěn)都與巖石裂紋擴(kuò)展密切相關(guān)［1-3］。巖石裂紋擴(kuò)展過程本身是從連續(xù)體到不連續(xù)體轉(zhuǎn)化的過程，巖石裂紋擴(kuò)展過程研究對巖體裂紋模型與本構(gòu)關(guān)系的建立、巖石體變形斷裂機(jī)制的揭示以及巖體工程穩(wěn)定性評價都具有重要意義［4］。

近幾年，隨著數(shù)值計算方法的不斷發(fā)展和計算機(jī)能力的不斷提高，數(shù)值分析方法在巖石裂紋擴(kuò)展研究中得到越來越多的應(yīng)用［5-6］。其中，數(shù)值流形法在處理巖石非連續(xù)性方面表現(xiàn)出了極大的優(yōu)勢，并逐漸得到了重視［7-8］。NING等［9］運(yùn)用數(shù)值流形法對巖體裂紋擴(kuò)展進(jìn)行了研究，該方法通過建立數(shù)學(xué)覆蓋和物理覆蓋系統(tǒng)，實(shí)現(xiàn)了對連續(xù)問題和非連續(xù)問題的統(tǒng)一求解，使得巖體裂紋擴(kuò)展不受網(wǎng)格劃分的影響，在巖石裂紋問題上表現(xiàn)出了很好的優(yōu)勢和前景。MA等［10］結(jié)合斷裂力學(xué)裂紋尖端的研究理論，根據(jù)裂紋尖端的應(yīng)力強(qiáng)度因子，以最大周應(yīng)力準(zhǔn)則為判據(jù)，判斷裂紋擴(kuò)展問題，在巖石爆破、水力壓裂等方面取得了一定的研究成果。張國新等［11-12］研究了二階流形法，借鑒有限單元法的二階位移函數(shù)的構(gòu)造方法，用6節(jié)點(diǎn)三角形構(gòu)成數(shù)學(xué)網(wǎng)格，用常數(shù)覆蓋函數(shù)和二階權(quán)函數(shù)構(gòu)成二階流形元法的位移函數(shù)，通過改進(jìn)的程序得出了二階流行法能夠較高精度的分析解決一般結(jié)構(gòu)的變形問題和接觸問題。徐棟棟等［13］通過對裂紋尖端附近的網(wǎng)格加密即通過在裂紋尖端影響區(qū)域內(nèi)的物理片上增加用于模擬應(yīng)力奇異性的增強(qiáng)位移函數(shù)的方法，保證了裂紋擴(kuò)展能夠在細(xì)分的網(wǎng)格中進(jìn)行，減小了裂紋擴(kuò)展對網(wǎng)格的依賴，擴(kuò)展了數(shù)值流形法。本文采用摩爾庫倫破壞準(zhǔn)則，對數(shù)值流形方法進(jìn)行分析，然后建立數(shù)值分析模型，對巖體在裂紋擴(kuò)展模式進(jìn)行分析研究，為巖體破裂研究提供新的思路。

1 數(shù)值流形方法

1.1 覆蓋位移原理及函數(shù)表達(dá)

數(shù)值流形方法是把有限覆蓋系統(tǒng)連接在一起，用以覆蓋全部材料體可用局部覆蓋所定義的函數(shù)來計算的方法。數(shù)學(xué)網(wǎng)格和物理網(wǎng)格相互重疊形成物理覆蓋。數(shù)學(xué)網(wǎng)格只定義數(shù)值解的精，它由用戶選擇；物理網(wǎng)格通常由材料的裂縫、邊界和分區(qū)界面等材料的物理及幾何性質(zhì)所決定。兩套覆蓋系統(tǒng)在處理非連續(xù)性問題時只對物理覆蓋進(jìn)行更新，數(shù)學(xué)網(wǎng)格并不更新。

采用有限單元覆蓋的數(shù)值流形方法，每個物理覆蓋都具有獨(dú)立的位移函數(shù)，因此在裂紋兩側(cè)由獨(dú)立的位移函數(shù)描述。圖1A表示含裂紋尖端的材料，流形單元被裂紋切割分為2種情況：完全切割（單元a），部分切割（單元b和單元c）。單元a是由3個數(shù)學(xué)覆蓋表示Mi（i=1，2，3），3個數(shù)學(xué)覆蓋形成6個獨(dú)立的物理覆蓋和（i=1，2，3），單元a被裂紋切割成2個獨(dú)立的流形單元a1（裂紋以上部分）和單元a（2裂紋以下部分），a1由物理覆蓋（i=1，2，3）描述，a2由物理覆蓋（i=1，2，3）描述，裂紋的相對位移為式（1）計算得到。

圖1 含裂紋的NMM

其中，?(i)為各物理覆蓋的勢函數(shù)，為不同的多項(xiàng)式，因此裂紋兩側(cè)就會發(fā)生相對位移。單元c是由3個數(shù)學(xué)覆蓋表示Mi（i=1，2，3），受到裂紋的影響，3個數(shù)學(xué)覆蓋形成4個獨(dú)立的物理覆蓋與，單元c被裂紋切割成2個流形單元，單元c1（裂紋以上部分）和單元c2（裂紋以下部分），c1由物理覆蓋（i=1，2，3）描述，c2由物理覆蓋、描述，裂紋的相對位移為

其中，?(1)為各物理覆蓋M1的勢函數(shù)，為不同的多項(xiàng)式，因此裂紋兩側(cè)有部分相對位移，這也就是裂紋尖端的奇異性。

1.2 塊體破裂算法

在數(shù)值流形開裂算法中，首先求解各流形單元的應(yīng)力狀態(tài)，由于流形單元面積太小，不能準(zhǔn)確表示開裂區(qū)域即裂紋尖端周圍的應(yīng)力狀態(tài)，因此以物理覆蓋為基本單位進(jìn)行開裂判斷，求出各個物理覆蓋應(yīng)力狀態(tài)，根據(jù)開裂準(zhǔn)則判斷是否發(fā)生開裂，若物理覆蓋應(yīng)力狀態(tài)達(dá)到開裂條件，則覆蓋中的單元發(fā)生開裂，根據(jù)覆蓋中最大主應(yīng)力和最大剪應(yīng)力以及塊體周邊現(xiàn)有的裂紋尖端情況，確定其開裂性質(zhì)以及開裂發(fā)生的方向和位置，規(guī)定開裂面在一個時步內(nèi)貫穿一個流形單元，并考慮單元破裂過程中的能量損耗。

以摩爾庫倫準(zhǔn)則為判據(jù)判斷裂紋擴(kuò)展，記壓應(yīng)力為正。在每個物理覆蓋上，最大主應(yīng)力、最小主應(yīng)力分別為σ1和σ3，根據(jù)摩爾應(yīng)力圓理論，在與最小主應(yīng)力σ3方向呈β角的面上的法相應(yīng)力和剪應(yīng)力分別為σ和τ，其中

巖石強(qiáng)度曲線服從庫倫準(zhǔn)則：

其中，c，φ分別為巖石粘聚力和內(nèi)摩擦角。

將式（3）、（4）代入式（5）可得

最小主應(yīng)力滿足-σ3=T0時，T0為巖石的抗拉強(qiáng)度，巖石發(fā)生剪切破壞。

2 巖體破裂模擬

將該方法應(yīng)用于巖體模型破裂模擬以驗(yàn)證其適用性，建立試樣如圖2A所示，試樣尺寸為5 cm×10 cm（寬×高），試樣中心處設(shè)置1條與水平方向夾角為45°、長13 mm、寬0.5 mm的預(yù)設(shè)裂紋，對應(yīng)建立的數(shù)值流形網(wǎng)格模型如圖2B所示。為了驗(yàn)證模型的正確性，巖體模型結(jié)果與ZHANG［14］的試驗(yàn)結(jié)果進(jìn)行對比，試樣的參數(shù)見表1。試驗(yàn)采用位移控制伺服加載函數(shù)，當(dāng)出現(xiàn)峰值應(yīng)力后，軸向應(yīng)變達(dá)到10%時，試驗(yàn)終止。

圖2 巖體單軸壓縮模型:（A）二維試樣中包含的單個預(yù)設(shè)裂縫（α=45°，l=13 mm）;（B）與物理試驗(yàn)相同高度和寬度的數(shù)值流形模型

表1 模擬試樣的力學(xué)參數(shù)[10]

2.1 壓縮條件下巖體破裂裂隙擴(kuò)展

圖3是含有預(yù)設(shè)裂紋的試樣模擬結(jié)果，模型黑色表示裂紋擴(kuò)展形態(tài)，圖3A-E中分別顯示的計算時步為10步、100步、200步、300步和400步，可以看出預(yù)設(shè)裂紋尖端首先產(chǎn)生微裂紋（圖3B），表明裂紋尖端有應(yīng)力集中，隨著試樣應(yīng)變的增大，裂紋逐步向試樣上下擴(kuò)展（圖3C-E），與ZHANG［14］的試驗(yàn)結(jié)果（圖3F）比較可以看出，本次模擬的裂紋擴(kuò)展模式和形態(tài)與巖體實(shí)際擴(kuò)展一致，次生裂紋均在裂紋尖端進(jìn)行啟裂，并以曲線路徑延伸擴(kuò)展至試樣頂部和底部，模擬結(jié)果說明本文所引入的摩爾庫倫破壞準(zhǔn)則對于模擬巖石破裂具有較好的效果。

圖3 巖石開裂過程模擬結(jié)果:(A)-(E)數(shù)值模擬結(jié)果;(F)實(shí)驗(yàn)結(jié)果[10]

2.2 應(yīng)力與裂紋尖端應(yīng)變特性分析

為了更好的對比預(yù)設(shè)裂紋對巖體強(qiáng)度的影響，在裂隙巖體模擬過程中，記錄加載板的應(yīng)力和應(yīng)變參數(shù)，并與完整巖體的應(yīng)力應(yīng)變結(jié)果進(jìn)行對比，繪制試樣的應(yīng)力-軸向應(yīng)變的關(guān)系曲線，如圖4所示。可以看出，對于完整的巖體試樣，存在明顯的彈性階段，該階段的曲線平直，屈服階段極短，破壞階段曲線驟降，表明完整試件的彈性及脆性較為明顯，峰值應(yīng)力為139.4 MPa。而預(yù)制裂隙對巖石試樣的強(qiáng)度與變形特性產(chǎn)生了明顯的影響。加載起始階段應(yīng)力與完整巖體接近，但在彈性階段，由于預(yù)制裂紋的存在，應(yīng)力顯著降低，且應(yīng)力曲線有著明顯的波動，峰值應(yīng)力60.6 MPa，較完整巖體下降了56.5%。

圖4 巖體應(yīng)力應(yīng)變曲線

在加載過程中記錄了所有試件裂紋尖端所在物理覆蓋的平均應(yīng)變，結(jié)果如圖5所示。根據(jù)裂紋可以看出，在計算小于195時步前（A區(qū)），裂紋尖端平均應(yīng)變較大，而從195步到260步之間（B區(qū)），裂紋尖端平均應(yīng)變迅速降低，而后又區(qū)域穩(wěn)定（C區(qū)）。根據(jù)對裂紋的類型進(jìn)行了劃分，具體裂紋類型可分為張拉裂紋和剪切裂紋。而單軸壓縮巖體在垂直于最大主應(yīng)力方向主要為拉伸狀態(tài)，這也表明該類型預(yù)制裂紋起始時剪切裂紋最后演變?yōu)閺埨土鸭y。

圖5 裂紋尖端平均應(yīng)變

3 結(jié)論

1）理論分析了數(shù)值流形法中覆蓋位移原理及函數(shù)表達(dá)，將塊體單元破壞過程中的數(shù)學(xué)網(wǎng)格、物理網(wǎng)格以及物理覆蓋的變化特征進(jìn)行了系統(tǒng)描述，為引入破壞本構(gòu)模型奠定了基礎(chǔ)。

2）將摩爾庫倫準(zhǔn)則引入數(shù)值流形法中，規(guī)定覆蓋中最大主應(yīng)力、最大剪應(yīng)力以及塊體周邊現(xiàn)有的裂紋尖端情況，確定其開裂性質(zhì)以及開裂發(fā)生的方向和位置，開裂面在一個時步內(nèi)貫穿一個流形單元，并考慮單元破裂過程中的能量損耗。

3）通過單軸壓縮試驗(yàn)驗(yàn)證了該破壞本構(gòu)的使用性，模擬結(jié)果表明該模型能夠較為準(zhǔn)確地反映巖體的破壞特征。預(yù)制裂隙對巖體試樣的強(qiáng)度與變形特性產(chǎn)生了明顯的影響。起始階段應(yīng)力與完整巖體接近，但在彈性階段，應(yīng)力顯著降低，且應(yīng)力曲線有著明顯的波動，峰值應(yīng)力較完整巖體下降顯著。