應(yīng)變軟化模型中內(nèi)聚力對模擬結(jié)果影響的數(shù)值模擬研究

杜學(xué)領(lǐng)

(貴州理工學(xué)院 礦業(yè)工程學(xué)院，貴州 貴陽 550003)

0 引 言

應(yīng)變軟化模型是巖體研究中常用的本構(gòu)模型，盡管該模型來源于摩爾庫倫模型，但在峰后特性方面，應(yīng)變軟化模型被認為是與真實巖體特性更加接近的模型。因此，在巖石力學(xué)及巖土、礦業(yè)工程等領(lǐng)域，應(yīng)變軟化模型應(yīng)用廣泛。

在FLAC3D數(shù)值模擬中，Itasca公司將摩爾庫倫模型設(shè)定為理想彈塑性模型，達到峰值強度后，一般將會保持穩(wěn)定；應(yīng)變軟化模型的峰前特性與摩爾庫倫模型幾乎沒有差異，但在峰后則表現(xiàn)為脆性突然破壞或階段性軟化破壞[1]。物理實驗中，往往能夠觀察到煤巖體的峰后軟化特性，如砂巖在三軸圍壓下的測試表明，應(yīng)變?yōu)?.005～0.04時達到峰值強度，而后進入峰后軟化階段[2]；化學(xué)腐蝕作用下，黑云片巖峰值載荷所對應(yīng)的應(yīng)變量為0.18～0.28，化學(xué)腐蝕作用使達到峰值載荷所需的應(yīng)變量增加[3]；不同加載速度下，巖石的峰值載荷應(yīng)變?yōu)?.28%～0.43%，而煤的應(yīng)變?yōu)?.38%～0.75%，甚至可達到1%[4]。李鵬飛等[5]通過對花崗巖的三軸實驗分析認為,巖石內(nèi)摩擦角會隨著塑性參數(shù)而變化，但殘余內(nèi)摩擦角與初始內(nèi)摩擦角值相差不大；陸銀龍等[6]通過軟弱泥巖的常規(guī)三軸試驗則認為，廣義內(nèi)摩擦角會隨著圍壓增大而顯著降低，但在峰后軟化階段基本保持不變；經(jīng)緯等[7]分析認為，巖石峰后內(nèi)摩擦角可近似看作是不變的；亦有研究指出含水率、粗顆粒含量、內(nèi)摩擦角等對斷層泥的軟化特性產(chǎn)生主要影響[8]?？梢?，目前對于內(nèi)摩擦角已有一定研究，但對于內(nèi)聚力的評價則較少。

工程方面，基于能量原理的分析認為，煤體的應(yīng)變軟化特性是沖擊地壓的必要條件，該研究列舉了煤的峰值強度對應(yīng)應(yīng)變約為0.005[9]；李鐵[10]研究認為，軟煤的反復(fù)軟化-硬化，使軟煤在深部高應(yīng)力條件也有可能發(fā)生沖擊失穩(wěn)，在沖擊傾向性評價中，動態(tài)破壞時間、沖擊能量指數(shù)、剩余能量指數(shù)等也與試件的峰后特性密切相關(guān)[11]；周勇等[12]研究認為，忽略應(yīng)變軟化特性對于工程而言有可能存在風(fēng)險；王凱等[13]也認為，采用應(yīng)變軟化模型要比摩爾庫倫模型更適用于深部工程。

上述物理實驗及相關(guān)研究表明，煤巖體的應(yīng)變軟化特性是客觀存在的，甚至?xí)こ淘斐芍匾绊?。但同時也說明，不同巖體達到峰值載荷的應(yīng)變量不同，相同巖體在不同加載速度、不同試樣處理條件下，峰值載荷所對應(yīng)的應(yīng)變量也不同。

數(shù)值模擬中主要通過設(shè)置內(nèi)聚力、剪脹角、內(nèi)摩擦角、抗拉強度、累積剪切塑性應(yīng)變、累積拉伸塑性應(yīng)變等參數(shù)實現(xiàn)應(yīng)變較化[1-14]。臨界抗拉強度與抗拉強度、內(nèi)聚力、內(nèi)摩擦角等關(guān)系密切，且峰后脆性特征與抗拉強度的設(shè)置只可二選一，而實際應(yīng)用中往往主要控制內(nèi)聚力、內(nèi)摩擦角、抗拉強度。FLAC3D采用混合離散法顯式動態(tài)有限差分解算，計算結(jié)果不僅與參數(shù)設(shè)置有關(guān)，而且還與整個模型的網(wǎng)格單元數(shù)量、不同時步的單元狀態(tài)和不同單元間的交互作用等有關(guān)。對于巖體的非線性問題，模擬結(jié)果往往與理論模型有所差異。以往研究偏重于工程視角、參數(shù)量級、因素分析等方面，大都采用莫爾庫倫模型[15-18]，卻未對比本構(gòu)模型差異對模擬結(jié)果的影響，另對應(yīng)變軟化模型的軟化參數(shù)設(shè)置的評價也略顯不足。本文以FLAC3D中的應(yīng)變軟化模型為基礎(chǔ)，著重探討內(nèi)聚力軟化對模擬結(jié)果的影響。

1 研究模型建立

根據(jù)摩爾庫倫破壞準則，F(xiàn)LAC3D中巖體破壞主要由以下條件控制[1]：

(1)

σ3=σt，

(2)

(3)

巖體有3種破壞條件：滿足式(1)時，由σ1，σ3，c，φ共同決定破壞準則；滿足式(2)時，最大主應(yīng)力達到巖體抗拉強度時巖體發(fā)生破壞；巖體的臨界拉伸強度為巖體抗拉強度或式(3)中的較小值時，巖體也會發(fā)生破壞。

由于巖體受到的主應(yīng)力由加載條件和圍巖網(wǎng)格單元的相互作用共同決定，而c，φ，σt在數(shù)值模擬中由人為指定，因此這3個參數(shù)對巖體破壞均會產(chǎn)生影響，應(yīng)變軟化參數(shù)的影響也會體現(xiàn)在網(wǎng)格單元的屬性變化上。本文以內(nèi)聚力作為主要研究對象，進一步探討其他參數(shù)對模擬結(jié)果的影響。

建立直徑50 mm×高100 mm的標準圓柱體試樣模型，數(shù)值模型切面圖如圖1所示。參照CHEN J等[19]在FLAC2D中的研究方法，非加載端限制其軸向位移，加載端采用位移控制條件，邊界位移速度為10-7m/step，取A，B，C3點，并利用位移差與長度比作為軸向應(yīng)變，利用加載端平面內(nèi)的不平衡力計算軸向應(yīng)力。為方便對比模擬結(jié)果，模型基本參數(shù)與文獻[19]中材料參數(shù)相同，即密度2 100 kg/m3，體積模量1.1 GPa，剪切模量1.2 GPa，內(nèi)聚力3.42 MPa，抗拉強度1 MPa。文獻[19]中材料的內(nèi)摩擦角僅為13.8°，該數(shù)值顯然與巖石物理實驗存在一定差距，本文將初始內(nèi)摩擦角設(shè)置為35°。在此基礎(chǔ)上，通過改變應(yīng)變軟化參數(shù)，觀察加載過程中不同模型的響應(yīng)，并以此評價軟化參數(shù)對模擬結(jié)果的影響。

圖1 數(shù)值模型切面示意圖

2 摩爾庫倫模型中內(nèi)聚力的影響

保持其他參數(shù)不變，在摩爾庫倫模型中按照25%遞減改變內(nèi)聚力，觀察不同內(nèi)聚力時模型的響應(yīng)情況。圖2為采用摩爾庫倫模型時，不同內(nèi)聚力下應(yīng)力-時步演化結(jié)果。

由圖2可知，采用摩爾庫倫模型的應(yīng)力演化具有以下特點：

圖2 摩爾庫倫模型中不同內(nèi)聚力下應(yīng)力-時步演化結(jié)果

(1)內(nèi)聚力越大，峰值強度越高。

(2)內(nèi)聚力與峰值強度呈線性相關(guān)。當(dāng)內(nèi)聚力以25%幅值降低時，試樣單軸抗壓強度的峰值強度也表現(xiàn)出線性降低。內(nèi)聚力分別為3.42，

2.57，1.71，0.86 MPa時，試樣的峰值強度分別約為13.2，9.9，6.6，3.4 MPa，下降幅值分別為3.3，3.3，3.2 MPa，考慮到小數(shù)位數(shù)的保留和峰值強度波動情況，可以近似認為，試樣峰值強度隨著內(nèi)聚力同步等幅度降低，而下降平均幅值3.3 MPa與最大的峰值強度13.2 MPa的比值也為25%，進一步說明若在模擬中保持其他參數(shù)不變而只改變內(nèi)聚力，試樣單軸抗壓強度會隨著內(nèi)聚力改變進行等幅度改變，內(nèi)聚力與峰值強度呈線性相關(guān)。

(3)摩爾庫倫模型具有彈塑性。內(nèi)聚力非0時，應(yīng)力-時步均表現(xiàn)出彈塑性特點，但與理想彈塑性模型相比有差異：其一是達到塑性穩(wěn)定狀態(tài)前存在一個應(yīng)力調(diào)整階段，一般在此階段出現(xiàn)最大應(yīng)力，但最大應(yīng)力與塑性階段相對穩(wěn)定的應(yīng)力值在數(shù)值上差別不大，可以近似認為二者相等；其二是塑性階段應(yīng)力值為波動穩(wěn)定狀態(tài)而非絕對穩(wěn)定狀態(tài)，它圍繞某一水平波動變化，振幅較小。

(4)內(nèi)聚力為0時的特殊性。當(dāng)網(wǎng)格單元內(nèi)聚力為0時，加載初期會產(chǎn)生一個較小的峰值應(yīng)力，約1.2 MPa，為最大峰值強度13.2 MPa的9%，但該應(yīng)力出現(xiàn)的范圍非常小，峰值點出現(xiàn)在前100時步以內(nèi)，第500時步時應(yīng)力就已經(jīng)衰減到0.000 5 MPa，此后均保持在非常小的水平。對于FLAC3D的網(wǎng)格單元，即便內(nèi)聚力為0，依然會存在力的相互作用，這與現(xiàn)實不完全相符；當(dāng)模擬中的內(nèi)聚力參數(shù)過小時，若采用應(yīng)變軟化模型時將內(nèi)聚力軟化至0，內(nèi)聚力為0時的應(yīng)力波動有可能會影響到模擬的準確性。

3 應(yīng)變軟化模型中內(nèi)聚力對模擬結(jié) 果影響

3.1 軟化標志點對模擬結(jié)果影響

分別設(shè)置應(yīng)變軟化標志點(SSP)為0.05，0.10，0.20，觀察不同標志點對應(yīng)不同內(nèi)聚力時的應(yīng)變軟化效果，軟化終點均為(1,0)。圖3為不同應(yīng)變軟化參數(shù)的模擬結(jié)果。由圖3可知：(1)僅改變應(yīng)變軟化參數(shù)時，不影響試樣在峰前的表現(xiàn)，不同軟化參數(shù)的峰前特性與摩爾庫倫模型表現(xiàn)一致，這與FLAC3D手冊中的峰后特性由塑性應(yīng)變控制相一致，且在摩爾庫倫模型中的峰后強度要高于應(yīng)變軟化后的強度；(2)其他條件相同時，峰后塑性應(yīng)變的標志點越小，峰后的應(yīng)變軟化效果越明顯;(3)其他條件相同時，峰后應(yīng)變軟化參數(shù)與峰前參數(shù)差值越大，峰后的應(yīng)變軟化效果越明顯；(4)采用應(yīng)變軟化模型后，峰后的應(yīng)力非

線性演化特征明顯，相同標志點時，峰后應(yīng)力跌落與內(nèi)聚力的衰減并非簡單的一元線性關(guān)系。相同軟化程度時，不同的軟化標志點與內(nèi)聚力的衰減也不是簡單的線性關(guān)系。標志點選取和軟化程度設(shè)置使峰后軟化特性更加復(fù)雜，以致不同的參數(shù)設(shè)置可能獲得類似的初期峰后表現(xiàn)；(5)應(yīng)變軟化值與軟化前的內(nèi)聚力差值較小時，不同軟化標志點時軟化效果均并不明顯。

3.2 軟化程度與摩爾庫倫模型比較

以軟化標志點0.05為例，將該組模擬運算時步增加至6萬，圖4為摩爾庫倫模型與應(yīng)變軟化模型的結(jié)果對比。由圖4可知，對于摩爾庫倫模型(MC)，達到應(yīng)力峰值后接近理想塑性狀態(tài)，應(yīng)力基本保持不變。對于應(yīng)變軟化的模型(SS)，峰后應(yīng)力表現(xiàn)為不斷降低，軟化程度越大，峰后初期應(yīng)力降低較快，但經(jīng)歷快速降低后，應(yīng)力會逐漸過渡到降低速度下降階段，此時峰后應(yīng)力緩慢降低。由圖4可以看出，峰后軟化特性與單純摩爾庫倫模型中內(nèi)聚力參數(shù)之間的關(guān)聯(lián)性并不顯著，如圖4所示，將內(nèi)聚力軟化至1.71 MPa時，與摩爾庫倫模型中內(nèi)聚力1.71 MPa相比，應(yīng)變軟化模型中的峰前依然呈現(xiàn)出內(nèi)聚力為3.42 MPa的特征，而在峰后也不會穩(wěn)定在1.71 MPa。由圖4可以看出，在應(yīng)變軟化模型中，內(nèi)聚力軟化至1.71，0.86 MPa，峰后應(yīng)力快速跌落，曲線快速跌落的終點值與摩爾庫倫模型保持穩(wěn)定的應(yīng)力值接近，而當(dāng)內(nèi)聚力軟化至2.57 MPa時，這種現(xiàn)象并不明顯。

當(dāng)內(nèi)聚力軟化程度足夠大、軟化標志點適當(dāng)小時，應(yīng)力快速跌落終點，與采用相同內(nèi)聚力的摩爾庫倫模型相當(dāng)。由圖3可以看出，軟化標志點較大或內(nèi)聚力軟化程度較小時，這種現(xiàn)象并不明顯。

產(chǎn)生上述現(xiàn)象的原因可能為：加載過程中巖體已經(jīng)積累一定應(yīng)力，應(yīng)力增加至某一時刻達到相應(yīng)的破壞條件，網(wǎng)格單元破壞并將屬性參數(shù)重置為軟化參數(shù)，該參數(shù)與采用摩爾庫倫模型的相同，軟化導(dǎo)致的參數(shù)重置和未破壞單元的應(yīng)力重分布使應(yīng)力快速跌落，且跌落終點與采用摩爾庫倫模型相同參數(shù)時相近。隨著加載不斷進行，破壞單元逐漸增加，軟化后的單元數(shù)量已處于主導(dǎo)地位，但達到軟化標志終點(1，0)的單元仍占少數(shù)，因此應(yīng)力跌落由快速跌落轉(zhuǎn)入緩慢跌落。對于軟化程度較低的參數(shù)，由于軟化后的內(nèi)聚力與軟化前差別不明顯，使峰后應(yīng)力表現(xiàn)為緩慢下降。圖3中將軟化標志點設(shè)為0.2時的情況與此類似，由于峰后加載過程中較多的網(wǎng)格單元并沒有達到軟化條件，使峰后相當(dāng)長的時步內(nèi)軟化效果并不明顯。

3.3 不同監(jiān)測點的應(yīng)變影響

文獻[19]提出利用軸向兩點位移差與其距離的比值計算軸向應(yīng)變，但這種方法并未評估測點是否會帶來評價誤差，本文選擇如圖1所示的A，B，C3點，并按照上述方法計算軸向應(yīng)變，以軟化參數(shù)(0.05，0.86)為例予以說明，圖5為該示例的應(yīng)變演化過程。由圖5可知，達到峰值點前，不同測點獲得的應(yīng)變演化幾乎一致，該部分所需運算時步與峰值點的時步一致，均為5 300時步左右。達到峰值點后，應(yīng)變演化開始出現(xiàn)明顯差異，包含測點A在內(nèi)的計算方法在峰后應(yīng)變呈現(xiàn)不同程度的增加趨勢，應(yīng)變計算的另一點距離測點A越近，峰后應(yīng)變增加的越快，即AB的應(yīng)變增速快于AC獲得的應(yīng)變增速，BC應(yīng)變則呈現(xiàn)逐漸下降趨勢，這種趨勢與AB,AC趨勢相反。說明即使對于同一模型，選用同樣的應(yīng)變軟化參數(shù)，選擇不同測點評價應(yīng)變變化時，其結(jié)果可能不同。這與不同模型中的破壞形態(tài)有關(guān)，圖6為內(nèi)聚力0.86時不同模型的位移結(jié)果。

圖5 應(yīng)變-時步演化過程

圖6 內(nèi)聚力為0.86 MPa時的位移結(jié)果

由圖6(a)可知，對于摩爾庫倫模型，在單軸壓縮條件下，以圓柱體中部垂直平面為對稱平面，位移近似呈對稱分布。當(dāng)選擇不同應(yīng)變軟化標志點時，位移演化趨勢開始發(fā)生變化，標志點選擇越小，破壞越接近45°角剪切破壞，圖6(d)為較為明顯的剪切破壞，圖6(c)表現(xiàn)出傾斜剪切破壞，圖6(b)不明顯。從位移變化量看，加載端位移大于非加載端的位移，軟化越明顯的模型其非對稱性也越明顯。對于圖6(c)～(d)的剪切破壞，當(dāng)測點位于不同的軸向位置時，獲得的應(yīng)變演化趨勢也不同。文獻[19]中選擇的兩測點位于圓柱體的上下兩部分，與本文的AC選擇類似，由圖5可知，AC在峰后的應(yīng)變增加速度低于AB，甚至在后期二者的應(yīng)變差值會越來越大?？赏ㄟ^選取多個測點求取平均的方式緩解這一問題，但測點應(yīng)與模型的破壞相一致，以避免出現(xiàn)BC在峰后直接應(yīng)變降低的情況。

應(yīng)變標志點的選擇與模型整體監(jiān)測獲得的應(yīng)變之間關(guān)聯(lián)性并不明顯。如圖5中選擇0.05為應(yīng)變標志點，AB，AC，BC應(yīng)變達到0.05的時步并不一致，甚至BC應(yīng)變不可能達到0.05。應(yīng)變標志點的選擇是相對的，與模型的實際應(yīng)變量無法一一對應(yīng)。

3.4 抗拉強度與內(nèi)聚力的協(xié)同軟化效果

保持內(nèi)聚力分別沿軟化路徑1：(0，3.42)→(0.05，0.86)→(1，0)、軟化路徑2：(0，3.42)→(0.05，2.54)→(1，0)不變，設(shè)置不同的抗拉強度軟化路徑，觀察其協(xié)同軟化效果。抗拉強度的軟化路徑分別為：方案1，同步軟化路徑，(0，1)→(0.05，0.25)→(1，0)；方案2，快速軟化路徑，(0，1)→(0.01，0.25)→(1，0)；方案3，慢速軟化路徑，(0，1)→(0.25，0.25)→(1，0)。由圖7可知，同時對抗拉強度與內(nèi)聚力進行軟化，并未發(fā)揮出協(xié)同軟化的效果，增加抗拉強度的軟化參數(shù)后，應(yīng)力-時步演化變化非常微小，可忽略不計。

圖7 抗拉強度與內(nèi)聚力的協(xié)同軟化模擬結(jié)果

由式(1)～(3)可知，從網(wǎng)格單元的破壞條件看，抗拉強度的軟化效果與式(2)中的最大主應(yīng)力有關(guān)，而與內(nèi)聚力關(guān)系并不大。當(dāng)模型中的最大主應(yīng)力達不到抗拉強度時，設(shè)置抗拉強度的軟化參數(shù)對模擬并無影響。即使最大主應(yīng)力達到抗拉強度，但若最大主應(yīng)力集中在一定區(qū)域且整體范圍較小時，設(shè)置抗拉強度的軟化參數(shù)對模擬的影響也是較小的。由圖8可知，設(shè)軟化前巖體的抗拉強度為1 MPa，達到巖體的峰值強度后，模型中較大范圍的網(wǎng)格單元最大主應(yīng)力低于1 MPa，因此，破壞準則式(2)在本文研究條件下發(fā)揮作用的范圍有限。結(jié)合最大主應(yīng)力和塑性區(qū)分布可以看出，該模型以剪切破壞為主，剪切破壞帶的角度約為45 °，剪切破壞帶的分布范圍與最大主應(yīng)力超過1 MPa的網(wǎng)格單元集中范圍相一致，同時，在該剪切破壞帶周圍零星分布有拉伸破壞，與前文分析一致。

3.5 關(guān)于內(nèi)摩擦角的討論

圖9 式(3)的結(jié)果圖

角的影響更顯著，也會使式(3)結(jié)果偏大，有可能造成無法觸發(fā)式(3)發(fā)揮作用。因此，對于軟巖而言，在FLAC3D中將內(nèi)聚力軟化即可達到相應(yīng)的軟化效果，無需對內(nèi)摩擦角和抗拉強度進行軟化。

3.6 沖擊傾向性評價討論

在評價沖擊傾向性時，動態(tài)破壞時間、沖擊能量指數(shù)、剩余能量指數(shù)均要利用巖體加載的峰后部分。常規(guī)的數(shù)值模擬中，摩爾庫倫模型峰后相對穩(wěn)定，而應(yīng)變軟化模型則對軟化參數(shù)具有較大的依賴性。因此對于確定的巖體參數(shù)，由于既定輸入?yún)?shù)并不能反映天然巖體的離散特性，在數(shù)值模擬中通過峰后能量反映沖擊傾向性很可能不夠準確。按照預(yù)期選定軟化參數(shù)時，已經(jīng)排除了天然離散性或人為限定峰后表現(xiàn)，二者都無法與天然巖體的離散特性相匹配。因此，采用常規(guī)數(shù)值模擬的本構(gòu)模型并利用峰后曲線計算耗散變形能、峰后破壞能密度，若以此評價沖擊傾向性，則其研究結(jié)論的可靠性值得懷疑。

4 結(jié) 論

(1)采用摩爾庫倫模型時，內(nèi)聚力大小與峰值強度呈線性相關(guān)。內(nèi)聚力為0時，模型中依然存在短時的力作用。

(2)僅改變應(yīng)變軟化參數(shù)時，不影響試樣的峰前表現(xiàn)，不同軟化參數(shù)時的峰前特性與摩爾庫倫模型表現(xiàn)一致，其他條件相同時，峰后塑性應(yīng)變的標志點越小，峰后的應(yīng)變軟化效果越明顯，峰后應(yīng)變軟化參數(shù)與峰前參數(shù)差值越大，峰后的應(yīng)變軟化效果越明顯，設(shè)置的應(yīng)變軟化值與軟化前的內(nèi)聚力差值較小時，在不同軟化標志點時軟化效果均并不明顯。

(3)當(dāng)內(nèi)聚力軟化程度足夠大、軟化標志點適當(dāng)小時，其應(yīng)力快速跌落的終點與采用相同內(nèi)聚力的摩爾庫倫模型相當(dāng)。但軟化標志點較大或內(nèi)聚力軟化程度較小時，這種現(xiàn)象均不明顯。

(4)達到峰值點之前，不同測點獲得的應(yīng)變演化幾乎完全一致，該部分所需運算時步與峰值點的時步具有一致性。但達到峰值點后，應(yīng)變演化開始出現(xiàn)明顯差異，對于同一模型、選用同樣的應(yīng)變軟化參數(shù)，選擇不同測點評價應(yīng)變變化時，其結(jié)果也不同。

(5)一般在FLAC3D中將軟巖內(nèi)聚力軟化，即可達到相應(yīng)的軟化效果。在常規(guī)數(shù)值模擬中通過峰后能量反映沖擊傾向性很可能不夠準確。