徐曉惠，姚再興

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

相應(yīng)的臨界應(yīng)力與臨界應(yīng)變分別是

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

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

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

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

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

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

1 煤柱模型

1.1 煤柱模型簡(jiǎn)化

為了計(jì)算方便，對(duì)煤柱作如下簡(jiǎn)化:忽略自重及頂?shù)装鍖?duì)煤柱兩端的橫向限制、煤柱的應(yīng)力狀態(tài)及應(yīng)變狀態(tài)在整個(gè)煤柱內(nèi)是均勻的、對(duì)煤柱的應(yīng)力及應(yīng)變分析簡(jiǎn)化為對(duì)煤柱內(nèi)任一點(diǎn)的應(yīng)力及應(yīng)變分析。圖1是承受等效壓力p的煤柱簡(jiǎn)化模型。

圖1 簡(jiǎn)化的煤柱模型

1.2 應(yīng)力特點(diǎn)

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

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

壓應(yīng)力p以壓為正。主應(yīng)力向量增量記作:

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

1.3 應(yīng)變特點(diǎn)

與3個(gè)主應(yīng)力及其方向分別對(duì)應(yīng)的是煤柱的3個(gè)主應(yīng)變。第一主應(yīng)變?chǔ)?反映煤柱側(cè)向膨脹;第三主應(yīng)變?chǔ)?反映煤柱的壓縮;煤柱的第二主應(yīng)變?chǔ)?始終為零。因此，主應(yīng)變向量增量記作

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

2 承載能力分析

2.1 彈性階段解

2.1.1 彈性本構(gòu)關(guān)系

在足夠小的壓力下，煤柱可以認(rèn)為處于彈性受力狀態(tài)，此時(shí)的彈性本構(gòu)關(guān)系如下:

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

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

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

2.1.2 初始屈服解

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

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

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

相應(yīng)的臨界應(yīng)力與臨界應(yīng)變分別是

2.2 塑性軟化分析

2.2.1 塑性本構(gòu)關(guān)系

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

材料進(jìn)入屈服階段的本構(gòu)關(guān)系滿足

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

Dp為塑性矩陣，定義為

取塑性體應(yīng)變?chǔ)萷作為內(nèi)變量κ時(shí)，則

參數(shù)A定義為

塑性體應(yīng)變?chǔ)萷的微分表達(dá)式為

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

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

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

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

黏聚力c對(duì)塑性體應(yīng)變?chǔ)萷的依賴關(guān)系是

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

圖2 黏聚力隨塑性體應(yīng)變變化的高斯曲線

ξ與c0之一，及φ0，c0，cm，θp0由實(shí)驗(yàn)給出，ξ與c0的關(guān)系是:

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

2.2.3 軟化限制的分析

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

圖3 A隨塑性體應(yīng)變的變化曲線

由于A＞0，要求在A的極小值A(chǔ)min時(shí)也能成立，即

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

(1)Amin對(duì)ξ的導(dǎo)數(shù)

說(shuō)明Amin關(guān)于ξ是增函數(shù)。

(2)Amin對(duì)cm的導(dǎo)數(shù)

說(shuō)明Amin關(guān)于cm是減函數(shù)。

(3)Amin對(duì)φ的導(dǎo)數(shù)

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

(4)Amin對(duì)E的導(dǎo)數(shù)

說(shuō)明Amin是對(duì)E的增函數(shù)。

(5)Amin對(duì)ν的導(dǎo)數(shù)

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

當(dāng)煤柱受力進(jìn)入塑性狀態(tài)時(shí)，本構(gòu)關(guān)系滿足公式 (20)，其中的應(yīng)力向量增量dσ與應(yīng)變向量增量dε分別由公式(3)和公式(4)確定。由于煤柱的軟化，不能事先給出煤柱的最大承壓力，可以按控制應(yīng)變方法進(jìn)行求解。因此，對(duì)于給定的狀態(tài)，彈塑性矩陣認(rèn)為是已知的，第三主應(yīng)變?cè)隽恳彩且阎模瑒t可以求解dε1，dσ2和dp。

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

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

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

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

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

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

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

3 煤柱承載能力實(shí)例分析

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

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

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

黏聚力隨內(nèi)變量的變化 (Gauss曲線)見(jiàn)圖4。黏聚力先增加后減小，對(duì)應(yīng)材料的強(qiáng)化和軟化。

圖4 黏聚力隨塑性體應(yīng)變的變化曲線

強(qiáng)度參數(shù)k隨內(nèi)變量的變化見(jiàn)圖5。強(qiáng)度參數(shù)k先增加后減小，對(duì)應(yīng)材料的強(qiáng)化和軟化。

圖5 k隨塑性體應(yīng)變的變化曲線

3.2 應(yīng)變強(qiáng)化及軟化曲線

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

圖6 煤柱壓力隨壓應(yīng)變的變化曲線

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

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

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

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

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

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

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

圖8 軟化有限元模擬的煤柱壓力隨壓應(yīng)變的變化曲線

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

5 結(jié)論與展望

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

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

(3)對(duì)有軟化特點(diǎn)的平面應(yīng)變問(wèn)題，對(duì)邊受均布載荷，可以通過(guò)數(shù)值積分，得到足夠精確的解，可以用作軟化彈塑性分析計(jì)算機(jī)程序的測(cè)試算例。

［1］朱建明，彭新坡，姚仰平，等.SMP準(zhǔn)則在計(jì)算煤柱極限強(qiáng)度中的應(yīng)用［J］.巖土力學(xué)，2010，31(9):2987－2990.

［2］潘 岳，張 勇，吳敏應(yīng)，等.非對(duì)稱開(kāi)采礦柱失穩(wěn)的突變理論分析［J］.巖石力學(xué)與工程學(xué)報(bào)，2006，25(S2).

［3］李江騰，曹 平.非對(duì)稱開(kāi)采時(shí)礦柱失穩(wěn)的尖點(diǎn)突變模型［J］.應(yīng)用數(shù)學(xué)和力學(xué)，2005，26(8):1003－1008.

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

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

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

［7］張少杰，王金安，魏現(xiàn)昊.煤柱過(guò)渡區(qū)的應(yīng)力遷移與沖擊礦壓顯現(xiàn)特征［J］.中國(guó)礦業(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.