彈塑性軟化煤柱的承載能力及其有限元模擬

徐曉惠，姚再興

(1.遼寧工程技術大學礦業(yè)學院，遼寧阜新123000;2.煤炭科學研究總院礦山安全技術研究分院，北京100013; 3.北京大學工學院，北京100080)

相應的臨界應力與臨界應變分別是

在煤炭資源的地下開采中，為承受上覆巖層的壓力，通常會形成長條形煤柱 (也叫礦柱)。煤柱的寬度不僅影響煤炭資源的采出率，也直接關系到巷道或采空區(qū)的安全。

傳統(tǒng)的煤柱設計是以強度理論為依據(jù)的，只要認為煤柱某點進入屈服狀態(tài)，便認為煤柱被破壞。顯然，這種設計偏于保守。朱建明等人［1］推導了在SMP(Spatially mobilized plane)準則條件下煤柱的極限強度計算公式，認為煤柱極限強度計算應該考慮中主應力的影響，否則偏低。

突變理論的引入，使煤柱的設計納入頂?shù)装褰M成的系統(tǒng)。潘岳等［2］引入突變理論分析非對稱開采礦柱的失穩(wěn)。李江騰等［3］提出了非對稱開采時礦柱失穩(wěn)尖點突變模型。王連國等［4］基于尖點突變模型研究了礦柱失穩(wěn)機理，提出了礦柱失穩(wěn)判據(jù)。

王學濱［5］對礦柱的受力與穩(wěn)定性進行了模擬，并考慮了不均勻性及礦柱端部的限制。張曉君［6］的模擬再現(xiàn)了采空區(qū)從變形到破壞的全過程。張少杰等［7］采用數(shù)值模擬、現(xiàn)場微震和應力監(jiān)測，研究了工作面回采時，煤體內(nèi)的應力分布及遷移規(guī)律，揭示了工作面的沖擊礦壓顯現(xiàn)特征。

對于長條形的煤柱，在煤柱走向方向上各物理力學參數(shù)及受力和變形都沒有變化，符合平面應變的條件。因此可以抽象為平面應變問題。假設煤柱材料的破壞滿足Drucker－Prager破壞準則［8］。

本文把應力應變的非線性看作由材料的內(nèi)變量決定，以經(jīng)典彈塑性理論為基礎推導了它們的關系。在積分表達時，采用了數(shù)值積分，并用具有軟化特點的彈塑性有限元進行了模擬。

1 煤柱模型

1.1 煤柱模型簡化

為了計算方便，對煤柱作如下簡化:忽略自重及頂?shù)装鍖γ褐鶅啥说臋M向限制、煤柱的應力狀態(tài)及應變狀態(tài)在整個煤柱內(nèi)是均勻的、對煤柱的應力及應變分析簡化為對煤柱內(nèi)任一點的應力及應變分析。圖1是承受等效壓力p的煤柱簡化模型。

圖1 簡化的煤柱模型

1.2 應力特點

按照拉應力為正，且3個主應力由大到小的規(guī)定，即σ1≥σ2≥σ3，第一主應力σ1為水平向指向采空區(qū)，滿足下式:

第二主應力σ2是煤柱的走向方向，為壓應力;第三主應力σ3為煤柱橫截面的壓應力p，即

壓應力p以壓為正。主應力向量增量記作:

其中，t表示轉(zhuǎn)置。

1.3 應變特點

與3個主應力及其方向分別對應的是煤柱的3個主應變。第一主應變ε1反映煤柱側(cè)向膨脹;第三主應變ε3反映煤柱的壓縮;煤柱的第二主應變ε2始終為零。因此，主應變向量增量記作

其中，t表示轉(zhuǎn)置。

2 承載能力分析

2.1 彈性階段解

2.1.1 彈性本構關系

在足夠小的壓力下，煤柱可以認為處于彈性受力狀態(tài)，此時的彈性本構關系如下:

式中，D為彈性矩陣;σ，ε分別是應力向量與應變向量。

考慮到煤柱的受力與變形特點，可以解得彈性狀態(tài)下的應力與應變關系如下:

式中，E是彈性模量;ν是泊松比。

2.1.2 初始屈服解

當煤柱從彈性狀態(tài)進入臨界塑性狀態(tài)時，應力滿足初始屈服面方程。把式(1)，(2)，(8)代入如下的Drucker－Prager屈服面方程

可求得煤柱進入屈服狀態(tài)時臨界壓力ps為

式中，J2是應力偏量的第二不變量;I1是應力的第一不變量，式中的強度參數(shù)α，k在初始內(nèi)變量處取得。I1，J2在主應力空間的表達式如下:

相應的臨界應力與臨界應變分別是

2.2 塑性軟化分析

2.2.1 塑性本構關系

與Drucker－Prager屈服面方程相應的屈服面函數(shù)記作

材料進入屈服階段的本構關系滿足

式中Dep表示彈塑性矩陣，定義為

Dp為塑性矩陣，定義為

取塑性體應變θp作為內(nèi)變量κ時，則

參數(shù)A定義為

塑性體應變θp的微分表達式為

式中，k'和α'分別是k和α對內(nèi)變量κ的導數(shù);C是彈性矩陣D的逆矩陣;G是剪切模量 (剛性模量)，G=E/(2(1+ν));K是體積模量，K= E/(3(1－2ν))。

2.2.2 強度參數(shù)的變化模型

通常實驗室給出的是Mohr－Coulomb形式的強度參數(shù)，即黏聚力c和內(nèi)摩擦角φ。通過這兩個參數(shù)再等效給出Drucker－Prager的強度參數(shù)α，k。這里采用Drucker－Prager圓與Coulomb六邊形外頂點重合的等效方法，即:

內(nèi)摩擦角φ(弧度，對應的角度表示為φ0)保持常數(shù)，滿足下列關系式:

黏聚力c對塑性體應變θp的依賴關系是

式中，cr表示殘余黏聚力;cm為凸出黏聚力，它是最大黏聚力cM與殘余黏聚力cr的差;為最黏塑性體應變，塑性體應變達到該值，黏聚力達到最大值cM;ξ稱為胖度，胖度越大高斯函數(shù)曲線越緩; c0表示初始黏聚力，其變化曲線如圖2所示。

圖2 黏聚力隨塑性體應變變化的高斯曲線

ξ與c0之一，及φ0，c0，cm，θp0由實驗給出，ξ與c0的關系是:

式 (24)中的α'為零，k'的表達式如下:

2.2.3 軟化限制的分析

在公式 (22)塑性矩陣Dp的定義中，參數(shù)A包含的材料彈塑性的重要性質(zhì)。強度參數(shù)按式(27)和式 (28)給出時，k'＞0反映的是材料的強化;k'=0反映的是材料的理想塑性;k'＜0反映的是材料的軟化。在材料的彈塑性模擬，特別是軟化模擬時，把A控制在合理的范圍內(nèi)非常重要。A隨內(nèi)變量的變化而變化，在時，A取得極大值;在時，A取得極小值。圖3是按實例分析中給定參數(shù)后A的變化情況。

圖3 A隨塑性體應變的變化曲線

由于A＞0，要求在A的極小值Amin時也能成立，即

需要注意，式中α是φ的函數(shù)。在數(shù)值計算中，Amin的合理性關系到計算能否進行下去或是否收斂的重要性。以下分析影響Amin取值的各個因素。

(1)Amin對ξ的導數(shù)

說明Amin關于ξ是增函數(shù)。

(2)Amin對cm的導數(shù)

說明Amin關于cm是減函數(shù)。

(3)Amin對φ的導數(shù)

式中，cm的數(shù)量級是106，K的數(shù)量級是109，ξ的數(shù)量級是 10－3。若成立，則關于 φ是增函數(shù);若 cm成立，則關于φ是減函數(shù)。

(4)Amin對E的導數(shù)

說明Amin是對E的增函數(shù)。

(5)Amin對ν的導數(shù)

2.3 塑性軟化的數(shù)值解

當煤柱受力進入塑性狀態(tài)時，本構關系滿足公式 (20)，其中的應力向量增量dσ與應變向量增量dε分別由公式(3)和公式(4)確定。由于煤柱的軟化，不能事先給出煤柱的最大承壓力，可以按控制應變方法進行求解。因此，對于給定的狀態(tài)，彈塑性矩陣認為是已知的，第三主應變增量也是已知的，則可以求解dε1，dσ2和dp。

由于彈塑性矩陣是內(nèi)變量的函數(shù)，與公式(20)對應的積分表達不易給出。數(shù)值積分的計算過程如下:

(1)首先依據(jù)公式(10)到公式(18)給出初始屈服狀態(tài)的解，此時塑性體應變是零。

(2)依據(jù)公式 (21)計算彈塑性矩陣。

(3)任意給出第三主應變增量Δε3＜0(其絕對值越小，精度越高);依據(jù)與公式(20)相應的增量式Δσ=DepΔε，求出相應的增量Δε1，Δσ2和Δp。

(4)依據(jù)與公式 (25)對應的增量式Δθp=［1，1，1］(Δε－CΔσ)確定塑性體應變的增量，進而求出塑性體應變。

(5)如果還需要繼續(xù)，轉(zhuǎn)到第 (2)步。

采用這種方法，事實上給出一個彈塑性軟化分析的精確解。這不僅給出煤柱承載能力，而且為軟化彈塑性分析計算機程序給出一個測試算例。

3 煤柱承載能力實例分析

3.1 輸入?yún)?shù)

殘余黏聚力與凸出黏聚力的比cr/cm=0.35;峰值黏聚力內(nèi)變量θp0=0.01;內(nèi)摩擦角φ0=15°;胖度ξ=0.01;凸出黏聚力cm=1MPa。

對于如上給出的參數(shù)，A隨內(nèi)變量的變化見圖3。可見，在這一算例中，A的最小值也是大于0的。

黏聚力隨內(nèi)變量的變化 (Gauss曲線)見圖4。黏聚力先增加后減小，對應材料的強化和軟化。

圖4 黏聚力隨塑性體應變的變化曲線

強度參數(shù)k隨內(nèi)變量的變化見圖5。強度參數(shù)k先增加后減小，對應材料的強化和軟化。

圖5 k隨塑性體應變的變化曲線

3.2 應變強化及軟化曲線

按數(shù)值積分給出煤柱壓力隨壓應變變化的曲線，見圖6。

圖6 煤柱壓力隨壓應變的變化曲線

可見，煤柱在2.2MPa的壓力下開始進入屈服狀態(tài);在約5.0MPa的壓力下達到最大承載能力;煤柱的殘余承載能力約1.4MPa。

4 彈塑性軟化煤柱壓縮有限元模擬

4.1 軟化材料彈塑性有限元原理

與傳統(tǒng)有限元程序相比，這一程序引入了弧長法，從而在延拓算法中不需要預先指定壓力的增加或是減小，這就有可能使煤柱壓力達到最大值時，煤柱的壓縮繼續(xù)增加，而煤柱的壓力減小，煤柱進入不穩(wěn)定平衡狀態(tài)。

4.2 計算模型及控制參數(shù)

根據(jù)對稱性，取煤柱的1/4建立有限元幾何模型，如圖7。左側(cè)約束水平向位移;下側(cè)約束豎向位移。上側(cè)加均布壓力。劃分三角形網(wǎng)格，網(wǎng)絡共計2690個，節(jié)點共計1416個。模擬得出的煤柱壓力與煤柱壓應變關系見圖8。

圖7 煤柱模擬的有限元模型

圖8 軟化有限元模擬的煤柱壓力隨壓應變的變化曲線

可見，煤柱在2.5MPa的壓力下開始進入屈服狀態(tài);在約4.8MPa的壓力下達到最大承載能力;煤柱的殘余承載能力約2.2MPa。

5 結論與展望

(1)在壓力作用下，具有軟化特性的煤柱開始表現(xiàn)出彈性性質(zhì)，隨后進入塑性狀態(tài)。在強化階段，承載的壓力仍能增加;達到峰值承載能力后，如果要使煤柱保持平衡，必須減小壓力，煤柱進入軟化階段，煤柱處于不穩(wěn)定的平衡狀態(tài)。

(2)煤柱按彈性極限或殘余承載力設計是非常保守的，在經(jīng)濟上不合理，在能取得較準確的物理力學參數(shù)的基礎上按峰值承載力設計可充分發(fā)揮煤柱的承載能力，但要取足夠大的安全系數(shù)。

(3)對有軟化特點的平面應變問題，對邊受均布載荷，可以通過數(shù)值積分，得到足夠精確的解，可以用作軟化彈塑性分析計算機程序的測試算例。

［1］朱建明，彭新坡，姚仰平，等.SMP準則在計算煤柱極限強度中的應用［J］.巖土力學，2010，31(9):2987－2990.

［2］潘 岳，張 勇，吳敏應，等.非對稱開采礦柱失穩(wěn)的突變理論分析［J］.巖石力學與工程學報，2006，25(S2).

［3］李江騰，曹 平.非對稱開采時礦柱失穩(wěn)的尖點突變模型［J］.應用數(shù)學和力學，2005，26(8):1003－1008.

［4］王連國，繆協(xié)興.基于尖點突變模型的礦柱失穩(wěn)機理研究［J］.采礦與安全工程學報，2006，23(2):137－140.

［5］王學濱.具有初始隨機材料缺陷的礦柱漸進破壞模擬［J］.中國礦業(yè)大學學報，2008，37(2):196－200.

［6］張曉君.礦柱及圍巖對采空區(qū)破壞影響的數(shù)值模擬研究［J］.采礦與安全工程學報，2006，23(1):123－126.

［7］張少杰，王金安，魏現(xiàn)昊.煤柱過渡區(qū)的應力遷移與沖擊礦壓顯現(xiàn)特征［J］.中國礦業(yè)，2013，22(5):73－78.

［8］Drucker，D.C.，Prager，W..Soil mechanics and plastic analysis for limit design［J］.Quarterly of Applied Mathematics，1952，10 (2):157－165.