蔣力帥，武泉森，李小裕，丁 楠

(1.山東科技大學(xué) 礦山災(zāi)害預(yù)防與控制國家重點實驗室培育基地，山東 青島 266590; 2.山東科技大學(xué) 礦業(yè)與安全工程學(xué)院，山東 青島 266590)

采動應(yīng)力與采空區(qū)壓實承載耦合分析方法研究

蔣力帥1,2，武泉森2，李小裕2，丁 楠2

(1.山東科技大學(xué) 礦山災(zāi)害預(yù)防與控制國家重點實驗室培育基地，山東 青島 266590; 2.山東科技大學(xué) 礦業(yè)與安全工程學(xué)院，山東 青島 266590)

采空區(qū)垮落帶在覆巖沉降作用下的應(yīng)力恢復(fù)對采動應(yīng)力場有顯著影響。基于圍巖應(yīng)力與采空區(qū)壓實承載的耦合作用，研究了垮落帶巖體壓實承載的力學(xué)特性及參數(shù)反演方法，分析了裂隙帶巖體的拉伸劣化力學(xué)特性，提出了采動應(yīng)力場與采空區(qū)壓實承載耦合分析方法，通過實例仿真模擬分析進(jìn)行了驗證。研究表明：垮落帶漸進(jìn)壓實承載的力學(xué)特性是指數(shù)型表達(dá)式，壓縮剛度隨壓縮應(yīng)變呈非線性增長。由雙屈服本構(gòu)模型反演得到與理論解擬合度較好的垮落帶力學(xué)參數(shù)，實現(xiàn)了垮落帶壓實承載特性的仿真模擬。裂隙巖體劣化系數(shù)表達(dá)式反映了模量劣化特性，應(yīng)變軟化本構(gòu)模型增加模量劣化算法組成巖體劣化本構(gòu)模型，實現(xiàn)了應(yīng)變軟化和裂隙巖體模量劣化的動態(tài)分析，較真實嚴(yán)謹(jǐn)?shù)姆从沉严稁r體的力學(xué)特性?；诳迓鋷c裂隙帶特征的仿真模擬，實現(xiàn)了采動應(yīng)力場與采空區(qū)覆巖“兩帶”特征的耦合分析。通過工作面支承壓力場演化特征的實例仿真分析，驗證了采動應(yīng)力場與采空區(qū)壓實承載耦合分析方法的可行性和正確性，提高了采動應(yīng)力場和圍巖采動響應(yīng)研究的嚴(yán)謹(jǐn)性和可靠性。

采動應(yīng)力;垮落帶壓實;裂隙帶劣化;數(shù)值模擬;耦合分析

長壁工作面垮落法管理頂板是最為普遍的采煤方法。隨著工作面推進(jìn)，采場圍巖應(yīng)力分布不斷調(diào)整及演化，采空區(qū)直接頂垮落、基本頂及其上部巖層破斷運移、高位巖層緩慢沉降等破壞及移動現(xiàn)象隨之出現(xiàn)，垮落的松散巖體充填采空區(qū)形成垮落帶，上覆巖層破斷運移形成裂隙帶，高位巖層直至地表形成彎曲下沉帶。工作面采空區(qū)冒落巖體在上覆巖層沉降作用下，由松散狀態(tài)逐漸壓實，支撐能力不斷提高，最終形成支撐上覆巖層的承載體。這一動態(tài)力學(xué)過程實質(zhì)上是圍巖應(yīng)力與覆巖運移及垮落帶逐漸壓實的耦合作用過程，圍巖應(yīng)力集中引起覆巖運移，覆巖運移又導(dǎo)致圍巖應(yīng)力分布演化。其中，支承壓力分布是采動應(yīng)力場的核心，是沿空巷道布置、回采巷道圍巖控制、沖擊地壓防治、煤與瓦斯突出防治等實際問題的理論基礎(chǔ)和重要依據(jù)[1-3]。

國內(nèi)外學(xué)者對采動應(yīng)力場的分布特征都十分關(guān)注，基于連續(xù)介質(zhì)的數(shù)值模擬方法被廣泛應(yīng)用。在數(shù)值模擬計算中，通常以開挖方式處理工作面采空區(qū)，并未考慮垮落帶的漸進(jìn)壓實特性及裂隙帶特征的影響;或者將采出空間或垮落帶賦予軟弱巖性，并不符合垮落帶壓實過程的力學(xué)特性;或者選擇適合垮落帶特性的本構(gòu)模型，但力學(xué)參數(shù)選取缺少依據(jù)。PAPPAS和C.MARK[4]指出采空區(qū)破碎松散的冒落巖體壓實力學(xué)特性對探究長壁工作面開采引起的圍巖擾動具有重要意義，并通過試驗研究了破碎松散巖體的壓實力學(xué)特性;張振南、馬占國等[5-7]通過松散巖塊的壓實試驗，分析了變形模量與軸向應(yīng)力、軸向應(yīng)變的關(guān)系;YAVUZ[8]以英國某煤礦的開采沉陷為背景，研究了采空區(qū)內(nèi)巖體壓實與應(yīng)力恢復(fù)主要與采高、埋深等參數(shù)有關(guān)，并采用FLAC數(shù)值計算軟件反演了采空區(qū)壓實力學(xué)特性;白慶升、屠世浩等[9]基于采空區(qū)壓實理論進(jìn)行了FLAC3D反演分析。近年來，一些學(xué)者采用連續(xù)介質(zhì)數(shù)值模擬方法，對采空區(qū)進(jìn)行等效模擬或近似等效處理。目前，相關(guān)研究主要集中在采空區(qū)內(nèi)冒落巖體的力學(xué)特性方面，數(shù)值模擬中未能正確反映垮落帶巖體的壓實力學(xué)特性，而采動影響下裂隙帶巖體的力學(xué)特征及其對采動應(yīng)力場的影響卻鮮有探索。

垮落帶松散巖體的壓實與承載力學(xué)特性，是采空區(qū)應(yīng)力恢復(fù)與分布規(guī)律的基礎(chǔ)，明顯影響著采動應(yīng)力場特別是支承壓力分布。覆巖裂隙帶范圍較大、裂隙縱橫，裂隙帶巖層的采動劣化對采動應(yīng)力場的影響是不應(yīng)忽略的。因此，在前人有關(guān)研究成果的基礎(chǔ)上，基于采動應(yīng)力場與覆巖破壞特征的耦合分析，綜合研究采動影響下垮落帶冒落巖體的壓實承載特性和裂隙帶巖體的拉伸劣化特性，通過對FLAC3D中本構(gòu)模型的二次開發(fā)，研究其在連續(xù)介質(zhì)數(shù)值模擬中的耦合再現(xiàn)技術(shù)，進(jìn)而提出較真實嚴(yán)謹(jǐn)?shù)牟蓜討?yīng)力場仿真研究方法，并結(jié)合工程實例驗證該方法的可行性及效果。研究工作對支承壓力分布、區(qū)段煤柱穩(wěn)定性、沿空巷道布置、回采巷道采動影響與控制設(shè)計、沖擊地壓危險評估、煤與瓦斯突出防治、巖層與地表移動等方面，具有重要的科學(xué)意義及應(yīng)用價值。

1 垮落帶壓實承載特性及參數(shù)反演方法

1.1 采動覆巖破壞高度

長壁工作面開采過程中，上覆巖層將不斷發(fā)生垮落及破斷運移，在采空區(qū)形成覆巖“三帶”結(jié)構(gòu)，其中垮落帶和裂隙帶直接影響著采動應(yīng)力分布和覆巖運移，“兩帶”高度及其力學(xué)特性是采動應(yīng)力場仿真中模擬極其重要的組成部分，對水體下采煤及下解放層開采等具有重要意義[10]。

采用鉆孔雙端堵水器、鉆孔沖洗液漏失量等方法，通過現(xiàn)場實測可得到具體工程條件下的覆巖“兩帶”高度。在大量現(xiàn)場實測分析的基礎(chǔ)上，國內(nèi)外學(xué)者研究提出了工作面垮落帶和裂隙帶高度的計算式。Palchik[11]結(jié)合烏克蘭頓涅茨克煤田的實踐經(jīng)驗與前人理論分析成果，認(rèn)為垮落帶高度與采高及碎脹系數(shù)有關(guān)，即：

式中，Hc為垮落帶高度;h為采高;b為碎脹系數(shù)。

BAI,KENDORSKI等[12]在中、美兩國大量現(xiàn)場實測數(shù)據(jù)的基礎(chǔ)上，通過統(tǒng)計回歸分析得到垮落帶和裂隙帶高度的回歸計算式，見式(2)。該計算方法包含不同巖體強度對“兩帶”高度的影響，精確性較高，因而得到廣泛應(yīng)用。

式中，σc為單軸抗壓強度;Hf為裂隙帶高度;c1,c2,c3,c4為巖層強度系數(shù)，見表1。

表1垮落帶和裂隙帶巖層強度系數(shù)
Table1Coefficientsforheightofcavingandfracturedzone

巖石類別σc/MPac1c2c3c4堅硬>402116122中硬20～4047191636軟弱<206232315

我國“三下”采煤規(guī)程中提供了覆巖“兩帶”高度的計算方法。金太等[13]綜合分析了兗州礦區(qū)的現(xiàn)場實測資料，得到緩傾斜厚煤層分層開采和綜放開采“兩帶”高度的變化規(guī)律與計算式。根據(jù)現(xiàn)場實測結(jié)果，或針對覆巖條件采用上述有關(guān)計算式，可得到比較可靠的“兩帶”高度，為建立數(shù)值計算模型提供依據(jù)。

1.2 垮落帶壓實承載過程的力學(xué)特性

隨著工作面推進(jìn)，直接頂板發(fā)生冒落及碎脹，充填采空區(qū)的垮落帶巖體最初呈現(xiàn)松散狀態(tài)，在裂隙帶及彎曲下沉帶的沉降荷載作用下，逐漸壓實同時支撐上覆運移巖層?？迓鋷r體的漸進(jìn)壓實過程及其承載力學(xué)特性，對采空區(qū)與煤體應(yīng)力狀態(tài)有著極其重要的影響。對松散破碎巖體壓實演變過程中的巖體力學(xué)響應(yīng)特性，國內(nèi)外學(xué)者進(jìn)行了理論和試驗研究[5-9,14]。由SALAMON[14]提出的破碎巖體壓縮過程的應(yīng)力-應(yīng)變關(guān)系見式(3)，被國內(nèi)外研究人員廣泛認(rèn)可并應(yīng)用[2,8,15-17]。

式中，σ為垮落帶巖體所受垂直應(yīng)力;E0為巖體初始正切模量;ε為垂直應(yīng)力作用下垮落帶巖體的應(yīng)變量;εm為巖體碎脹后最大應(yīng)變。

其中參數(shù)E0主要受巖體碎脹系數(shù)b及σc的影響。YAVUZ[8]基于前人大量不同巖石單軸壓縮試驗數(shù)據(jù)，通過三維回歸分析得到E0的表達(dá)式：

巖體碎脹后達(dá)到的最大應(yīng)變εm可由碎脹系數(shù)b得到

聯(lián)立式(1),(2)可得碎脹系數(shù)與巖體強度參數(shù)及采高的關(guān)系：

由式(3),(4),(5)得到垮落帶巖體壓實特性的表達(dá)式：

表達(dá)式(7)反映了垮落帶巖體漸進(jìn)壓實-支撐過程的力學(xué)特性。由式(7)及式(3)～(5)可知，垮落帶巖體壓實過程中的應(yīng)力-應(yīng)變關(guān)系由碎脹系數(shù)b和單軸抗壓強度σc確定，其中碎脹系數(shù)對壓實力學(xué)特性的影響更為顯著。

巖體單軸抗壓強度為30 MPa時碎脹系數(shù)b對垮落帶巖體剛度的影響如圖1所示?？梢?，垮落帶巖體的壓縮剛度隨著壓縮應(yīng)變的增長而明顯提高，巖體碎脹系數(shù)對碎脹后的最大壓縮應(yīng)變及壓縮剛度影響很大。高碎脹系數(shù)巖體的初期壓縮剛度上升趨勢平緩，后期上升幅度很大，壓縮應(yīng)變較大;低碎脹系數(shù)巖體的壓縮剛度上升較快，壓縮應(yīng)變較小。

圖1 碎脹系數(shù)對垮落帶巖體壓縮剛度的影響Fig.1 Effect of bulking factor on the compression stiffness of caved rock

碎脹系數(shù)為1.3和1.5時單軸抗壓強度對垮落帶巖體剛度的影響如圖2所示?？梢姡谒槊浵禂?shù)相同的前提下，高強度垮落巖體的壓縮剛度明顯大于低強度巖體，相同應(yīng)力狀態(tài)下高強度垮落巖體的應(yīng)變總是小于低強度巖體。需要注意的是，垮落巖體的體積永遠(yuǎn)無法壓實至初始完整體積，仍保留殘余碎脹。

圖2 單軸抗壓強度對垮落帶巖體壓縮剛度的影響Fig.2 Effect of uniaxial compressive strength on the compression stiffness of caved rock

上述研究得到的垮落帶巖體應(yīng)力-應(yīng)變關(guān)系式，揭示了垮落帶巖體壓實承載的力學(xué)特性，能夠可靠地反映出從松散低剛度狀態(tài)到壓實高剛度狀態(tài)的力學(xué)特性演化，表征了垮落帶巖體在覆巖壓力作用下逐漸壓緊密實的工程實際過程，定量表達(dá)了整體剛度隨體積減小到最終壓實的指數(shù)型增長規(guī)律，可表達(dá)垮落帶巖體壓縮、承載及支撐過程的力學(xué)行為。

1.3 垮落帶模型力學(xué)特性參數(shù)的反演方法

在數(shù)值仿真模擬研究中，正確再現(xiàn)采空區(qū)垮落巖體壓實-支撐過程的力學(xué)響應(yīng)，是真實模擬工作面開采過程中支承壓力分布、煤柱應(yīng)力狀態(tài)、巷道采動影響及其圍巖穩(wěn)定性的重要基礎(chǔ)。針對采空區(qū)垮落帶的數(shù)值模擬，國內(nèi)外學(xué)者進(jìn)行了大量探索，主要方法是采用彈性模型或雙屈服本構(gòu)模型(double-yield model)[8-9]充填采空區(qū)和垮落帶，其中雙屈服模型準(zhǔn)確體現(xiàn)出垮落帶巖體壓實過程中的力學(xué)特性，近年來被科研人員所認(rèn)可和運用。

將巖石力學(xué)試驗所得到的參數(shù)代入式(3)～(7)，可得出該工作面頂板條件下采空區(qū)垮落帶巖體碎脹后的最大應(yīng)變εm、碎脹系數(shù)b和壓實過程中的應(yīng)力-應(yīng)變關(guān)系。建立FLAC3D垮落碎脹巖體精細(xì)數(shù)值計算模型，采用雙屈服本構(gòu)模型進(jìn)行壓縮模擬試驗，監(jiān)測模型的應(yīng)力-應(yīng)變關(guān)系，與垮落帶應(yīng)力-應(yīng)變關(guān)系的理論解對比，采用試錯反演法可以得到擬合度較好的采空區(qū)模型力學(xué)參數(shù)。

2 裂隙帶巖體的拉伸劣化力學(xué)特性

巖體抗拉強度很低，且遠(yuǎn)小于抗壓強度及抗剪強度，在拉應(yīng)力作用下極易使節(jié)理與裂隙擴展甚至貫通，工作面圍巖特別是頂?shù)装鍘r層中會產(chǎn)生拉應(yīng)力及拉伸破裂。隨著工作面推進(jìn)及垮落帶碎脹巖體充填采空區(qū)，上覆懸露巖層在彎曲拉應(yīng)力作用下發(fā)生撓曲下沉及離層，在拉應(yīng)力超過巖層抗拉強度的部位，拉伸裂隙擴展甚至產(chǎn)生拉破裂，破斷運移后得到垮落帶碎脹巖體的支撐，形成拉伸裂隙及離層裂隙發(fā)育的裂隙帶。

國內(nèi)外學(xué)者紛紛指出巖石完整性、裂隙產(chǎn)狀和發(fā)育程度直接影響到巖體的力學(xué)性能，并分別提出了巖體力學(xué)參數(shù)(尤其是楊氏模量)與裂隙發(fā)育程度的關(guān)系[18-19]。由于采動裂隙帶的復(fù)雜性和隱蔽性，現(xiàn)場實測其裂隙發(fā)育程度難度較大、成本較高，對裂隙帶內(nèi)裂隙發(fā)育程度的準(zhǔn)確、定量描述非常困難，筆者基于國際普遍認(rèn)可和廣泛應(yīng)用的地質(zhì)強度指標(biāo)GSI(Geological Strength Index)體系，研究得到巖體拉伸破壞后的殘余楊氏模量與拉應(yīng)力作用下的裂隙發(fā)育程度GSIt的關(guān)系[20]，即工程巖體拉伸劣化系數(shù)的表達(dá)式：

式中，Er為巖體發(fā)生拉伸破壞后的殘余楊氏模量;Em為巖體初始楊氏模量;GSIt為巖體在拉應(yīng)力作用下的裂隙發(fā)育程度。

由式(8)可見，隨著巖體裂隙發(fā)育程度的增加，楊氏模量呈現(xiàn)快速下降的變化趨勢，具有顯著的模量劣化特性。裂隙帶巖體的力學(xué)特性進(jìn)入了峰后力學(xué)性質(zhì)的劣化階段，特別是裂隙發(fā)育導(dǎo)致的楊氏模量劣化。

工作面圍巖破壞的主要形式是剪破裂和拉伸破裂，圍巖峰后特性軟化不僅是剪破裂狀態(tài)下的應(yīng)變軟化，還有拉破裂狀態(tài)下的楊氏模量劣化?；贛ohr-Coulomb強度準(zhǔn)則的應(yīng)變軟化本構(gòu)模型，通過定義力學(xué)參數(shù)與塑性變形的負(fù)相關(guān)關(guān)系，較好地反映出巖體峰后力學(xué)性能軟化、并呈現(xiàn)殘余強度的力學(xué)行為。而拉破壞分析僅判別出單元的破壞狀態(tài)與破壞類型，尚未實現(xiàn)楊氏模量的劣化。由于FLAC3D中內(nèi)置的Mohr-Coulomb模型和應(yīng)變軟化模型中楊氏模量均維持恒定，不能準(zhǔn)確嚴(yán)謹(jǐn)?shù)姆从沉严栋l(fā)育、縱橫的裂隙帶巖體力學(xué)特性。為此，基于連續(xù)介質(zhì)有限差分?jǐn)?shù)值模擬軟件FLAC3D，采用FISH語言對內(nèi)置的應(yīng)變軟化本構(gòu)模型進(jìn)行二次開發(fā)，增加拉伸裂隙模量劣化算法，形成工程巖體劣化模型[20]。

工程巖體劣化模型的算法核心由巖體破壞識別算法和破裂劣化算法組成，該模型在常規(guī)顯式有限差分計算中，通過每間隔一定運算時步后遍歷三維模型的全部單元體，逐步漸進(jìn)地動態(tài)識別每個剪切破裂和拉伸破壞的單元，并對其進(jìn)行相應(yīng)力學(xué)參數(shù)劣化，其中黏聚力、內(nèi)摩擦角弱化遵循應(yīng)變軟化模型，楊氏模量劣化遵循拉伸劣化模型，單元力學(xué)參數(shù)更新后繼續(xù)迭代運算，循環(huán)上述過程直至模型達(dá)到平衡狀態(tài)[20]。

工程巖體劣化模型在動態(tài)識別圍巖破壞性態(tài)的基礎(chǔ)上，根據(jù)破壞類型進(jìn)行應(yīng)變軟化和楊氏模量劣化，通過力學(xué)特性等效的方法實現(xiàn)了巖體拉伸破裂后根據(jù)裂隙發(fā)育程度對楊氏模量的劣化，從而較真實嚴(yán)謹(jǐn)?shù)姆从沉严稁r體的力學(xué)特性，對覆巖裂隙帶模擬與采動應(yīng)力場演化研究具有重要意義。

3 采動應(yīng)力場與覆巖兩帶耦合分析

3.1 采動應(yīng)力場與覆巖兩帶耦合方法

基于上述采動覆巖“兩帶”破壞特征及其力學(xué)機制，考慮采空區(qū)垮落帶壓實承載力學(xué)特性與裂隙帶巖體破斷劣化特性的動態(tài)演化力學(xué)過程，研究提出采動應(yīng)力場與采動覆巖的耦合分析方法，如圖3所示。

圖3 采動應(yīng)力場與覆巖兩帶耦合分析方法Fig.3 Diagram of numerical simulation method on mining- induced stress field

根據(jù)工作面地質(zhì)條件和巖石物理力學(xué)性質(zhì)，計算或?qū)崪y得到垮落帶范圍，依據(jù)垮落帶巖體壓實過程的力學(xué)特性，建立垮落巖體精細(xì)數(shù)值模型，采用試錯反演得到符合垮落帶壓實過程應(yīng)力-應(yīng)變關(guān)系的模型力學(xué)參數(shù)，通過以上步驟實現(xiàn)采空區(qū)垮落帶壓實過程的力學(xué)特性模擬。計算或?qū)崪y得到裂隙帶范圍，根據(jù)裂隙帶巖體的拉伸劣化特性，實現(xiàn)上覆裂隙帶中裂隙發(fā)育過程的模擬。在連續(xù)介質(zhì)的數(shù)值模擬中，實現(xiàn)垮落帶壓實和裂隙帶力學(xué)特性的仿真，使采動應(yīng)力場分析與采空區(qū)覆巖力學(xué)特征相耦合，從而可以提高采動影響研究過程的嚴(yán)謹(jǐn)性和結(jié)果的可靠性。

3.2 采空區(qū)覆巖力學(xué)特性參數(shù)分析

為進(jìn)行采動應(yīng)力場與采空區(qū)覆巖“兩帶”力學(xué)特性的耦合分析，驗證上述方法的可行性，以河南能化集團(tuán)趙固二礦工程地質(zhì)條件為例，運用上述方法分析確定采空區(qū)覆巖力學(xué)特性參數(shù)。

趙固二礦處于焦作煤田東部，11030工作面開采二1煤層，煤層傾角0°～11°，平均厚度6.16 m，煤層結(jié)構(gòu)簡單，層位穩(wěn)定，屬近水平穩(wěn)定厚煤層。工作面采用大采高綜合機械化開采工藝，后退式全部垮落法管理頂板，開采深度平均652 m。直接頂以泥巖、砂質(zhì)泥巖為主，基本頂以粉砂巖、細(xì)粒砂巖為主，底板以砂質(zhì)泥巖為主。二1煤層和頂板的巖石試樣采自煤層頂板20 m和底板15 m范圍，采用RMT-150型伺服試驗機進(jìn)行巖石力學(xué)試驗，試驗測試結(jié)果見表2。

表2煤巖物理力學(xué)參數(shù)
Table2Physicalandmechanicalpropertiesofrocksandcoal

巖石層位巖性視密度/(kg·m-3)抗拉強度/MPa抗壓強度/MPa楊氏模量/GPa泊松比砂巖25521091836316023頂板砂質(zhì)泥巖2591347562163025泥巖277723438295029二1煤煤層1435072204281030底板砂質(zhì)泥巖262928141899027石灰?guī)r27541352996826022

11030工作面采高6 m，根據(jù)巖石力學(xué)參數(shù)試驗結(jié)果表2，由式(2)和表1，計算得到工作面采空區(qū)的垮落帶高度為12.7 m，裂隙帶高度為45.5 m。將表2中有關(guān)巖石力學(xué)參數(shù)代入計算式(5),(6)，得到二1煤層工作面采空區(qū)垮落帶巖體碎脹后的最大應(yīng)變εm為0.32，碎脹系數(shù)b為1.47。將εm及b兩個參數(shù)代入計算式(3),(4)，得到垮落帶巖體壓實過程的應(yīng)力-應(yīng)變關(guān)系如圖4所示。據(jù)此，采用FLAC3D數(shù)值計算模型對垮落帶巖體進(jìn)行壓縮模擬試驗，采用試錯反演法得到與理論解擬合較好的模型應(yīng)力-應(yīng)變關(guān)系如圖4所示，數(shù)值模擬反演得到的垮落帶雙屈服本構(gòu)模型參數(shù)見表3。

圖4 雙屈服模型力學(xué)特性的數(shù)值模擬反演Fig.4 Back analysis of stress-strain characteristics of double-yield model

K/MPaG/MPaγ/(kg·m-3)c/MPaφ/(°)σt/MPa19900100017000001300

由巖石力學(xué)參數(shù)表2，采用Hoek-Brown強度準(zhǔn)則計算得到巖體力學(xué)參數(shù)見表4。其中K為體積模量，G為剪切模量，c為黏聚力，φ為內(nèi)摩擦角，cr為峰后殘余內(nèi)聚力，εp為巖體強度變?yōu)闅堄嘀禃r的塑性應(yīng)變。數(shù)值計算模型的煤層和巖層采用表4的力學(xué)參數(shù)。工作面覆巖屬于中硬覆巖條件，通過分析研究[21]，裂隙帶巖體GSIt取值50。

表4模型巖體力學(xué)參數(shù)
Table4Rockmassmechanicalproperties

巖石層位巖性K/GPaG/GPac/MPaφ/(°)cr/MPaεp/%砂巖91593945039001頂板砂質(zhì)泥巖52313240032001泥巖24112135021001二1煤煤13061431014001底板砂質(zhì)泥巖72403437034001石灰?guī)r96654247042001

4 應(yīng)用實例

為驗證上述方法的正確性和可靠性，以趙固二礦11030工作面工程為背景建立模型算例，應(yīng)用上述模擬方法研究工作面采動應(yīng)力場演化規(guī)律，并與理論分析結(jié)果進(jìn)行對比印證。

4.1 數(shù)值模型的建立

根據(jù)對稱性原則，以11030工作面中線為對稱軸，建立三維數(shù)值模型如圖5(a)所示。模型長度350 m，其中工作面推進(jìn)長度250 m，前后各留50 m邊界;模型寬度180 m，工作面二分之一的長度90 m，側(cè)向煤柱寬度90 m;模型高120 m。模型頂部施加15 MPa的垂直應(yīng)力，X,Y方向施加的水平應(yīng)力分別為垂直應(yīng)力的0.8和1.2倍[22]，模型四周和底部采用位移限定邊界。隨著工作面推進(jìn)，垮落帶冒落巖體由雙屈服本構(gòu)模型填充，力學(xué)參數(shù)見表3。各巖層采用工程巖體劣化模型作為力學(xué)判據(jù)。

圖5 采動應(yīng)力場演化FLAC3D模型Fig.5 FLAC3D for mining induced stress field simulation

在煤層直接頂布置如圖5(b)中虛線所示的4條應(yīng)力測線，其中測線1,2,3沿推進(jìn)方向布置，分別位于工作面對稱軸、距工作面邊界45 m處和工作面邊界;測線4垂直于推進(jìn)方向布置，從采空區(qū)邊界貫穿邊界煤柱，工作面推進(jìn)200 m后記錄工作面前后和煤柱的垂直應(yīng)力分布。

4.2 工作面前后支承壓力演化規(guī)律

一次采高為3 m和6 m的工作面前后垂直應(yīng)力分布如圖6所示?？梢?，工作面開采促使圍巖內(nèi)部應(yīng)力重新調(diào)整和分布，在工作面前方煤壁前方產(chǎn)生超前支承壓力區(qū)，采空區(qū)垮落帶巖體在上覆巖層作用下逐漸壓實，出現(xiàn)卸壓與應(yīng)力恢復(fù)現(xiàn)象，隨著遠(yuǎn)離工作面，采空區(qū)內(nèi)垂直應(yīng)力逐漸增大。

圖6 工作面垂直應(yīng)力分布平面圖Fig.6 Plan view of vertical stress distribution induced by mining

工作面中部、距邊界45 m及工作面邊界位置(即測線1,2,3)的垂直應(yīng)力分布如圖7所示?？梢?，在工作面不同部位，前方煤壁和后方采空區(qū)的垂直應(yīng)力分布存在著一定差異，沿著工作面長度范圍(x=0～90 m)并非均勻分布。在整個工作面范圍內(nèi)，煤壁和采空區(qū)的垂直應(yīng)力峰值點均位于工作面中部，由中部向邊界逐漸降低。

圖7 工作面前后支承壓力分布Fig.7 Distribution of mining induced abutment pressure

由圖6,7可見，采高增大對前支承壓力分布有顯著影響，主要體現(xiàn)為煤壁應(yīng)力降低區(qū)擴大、支承壓力峰值降低、峰值影響區(qū)域前移、影響范圍擴大。采高增大對工作面中部前方煤壁影響最顯著，采高由3 m增加到6 m時，前支承壓力峰值由32.3 MPa降低至30.9 MPa，峰值位置前移8.4 m;在距邊界45 m處，前支承壓力峰值由30.6 MPa降低至29.4 MPa，峰值位置前移4.8 m;在工作面邊界，前支承壓力峰值僅由27.9 MPa降低至27.7 MPa，比中部峰值分別下降了4.4 MPa和3.2 MPa，峰值位置前移0.6 m。

由圖6,7反映出工作面后方采空區(qū)冒落帶壓實造成的應(yīng)力恢復(fù)。采空區(qū)應(yīng)力恢復(fù)在工作面長度方向有著顯著差異，采空區(qū)中部的應(yīng)力恢復(fù)較高，距離采空區(qū)邊界越近應(yīng)力恢復(fù)越低。采高的增大也有明顯影響，在采空區(qū)中部工作面后方110 m處，采高3 m時后支承壓力恢復(fù)到原巖應(yīng)力的88.7%，而采高6 m時僅恢復(fù)到79.3%，應(yīng)力恢復(fù)速度較慢。這是由于垮落帶和裂隙帶的高度與采高成正比(見式(1),(2))，大采高開采擴大了覆巖破壞范圍，使采動覆巖需要更長的時間才能經(jīng)過應(yīng)力重新調(diào)整達(dá)到穩(wěn)定。

在距邊界45 m的采空區(qū)部位，采高3和6 m時的采空區(qū)應(yīng)力恢復(fù)差異不大，后支承壓力升高速度基本相近，在工作面后方110 m處后支承壓力分別恢復(fù)到原巖應(yīng)力的73.7%和67.4%。而在采空區(qū)邊界位置，應(yīng)力恢復(fù)特征與采空區(qū)中部相反，采高為3 m時，后支承壓力在工作面后方110 m處恢復(fù)到原巖應(yīng)力的20.1%，而采高為6 m時恢復(fù)到30.0%，這一現(xiàn)象與其他學(xué)者的研究結(jié)果一致[8,23]。這是由于采空區(qū)邊界位于側(cè)向基本頂結(jié)構(gòu)下方，圍巖應(yīng)力調(diào)整過程中應(yīng)力主要向采空區(qū)中部壓實區(qū)和煤柱轉(zhuǎn)移，因此采空區(qū)邊界應(yīng)力恢復(fù)程度顯著低于采空區(qū)中部。由于大采高工作面垮落帶及裂隙帶高度的增加，采空區(qū)中部大范圍的覆巖破壞減緩了應(yīng)力恢復(fù)速度，使圍巖應(yīng)力調(diào)整時部分應(yīng)力由采空區(qū)中間向兩側(cè)轉(zhuǎn)移;同時采高的增大使側(cè)向基本頂結(jié)構(gòu)高度增加，使采空區(qū)邊界位置的覆巖壓力隨之增大。

4.3 側(cè)向煤柱集中應(yīng)力與力學(xué)狀態(tài)演化規(guī)律

工作面開采除了在前方煤壁形成前支承壓力，同時也在采空區(qū)側(cè)向煤柱上形成垂直應(yīng)力集中，煤柱邊緣區(qū)域在集中應(yīng)力作用下發(fā)生塑性破壞而卸載，集中應(yīng)力逐漸向煤柱深部轉(zhuǎn)移，塑性破壞區(qū)或應(yīng)力極限平衡區(qū)也逐步向煤柱深部擴展直至彈性應(yīng)力區(qū)邊界[24-25]。

一次采高為3和6 m時側(cè)向煤柱垂直應(yīng)力分布如圖8(a)所示，走向中部煤柱垂直應(yīng)力分布如圖8(b)所示。可見，采高對側(cè)向煤柱的應(yīng)力分布影響非常顯著。采高為3 m時，側(cè)向煤柱應(yīng)力峰值位置進(jìn)入煤柱內(nèi)部12.9 m處，即極限平衡區(qū)寬度12.9 m，應(yīng)力峰值34.1 MPa;采高為6 m時，側(cè)向煤柱上應(yīng)力集中區(qū)整體向煤柱深部移動，應(yīng)力極限平衡區(qū)寬度24.2 m，且應(yīng)力集中范圍擴大，應(yīng)力峰值為34.7 MPa。一次采出厚度由3 m增大至6 m時，側(cè)向煤柱中應(yīng)力極限平衡區(qū)寬度由12.9 m擴大到24.2 m，擴大系數(shù)為1.88。因此，采高大小對側(cè)向煤柱應(yīng)力極限平衡區(qū)寬度、集中應(yīng)力程度有著顯著影響，是沿空巷道合理煤柱寬度、穩(wěn)定性分析與支護(hù)設(shè)計的重要因素。

圖8 側(cè)向煤柱垂直應(yīng)力分布Fig.8 Distribution of vertical stress in pillar

煤柱應(yīng)力極限平衡區(qū)寬度是設(shè)計煤柱尺寸、計算頂板側(cè)向斷裂位置和確定沿空掘巷合理位置等工作的重要指標(biāo)。應(yīng)力極限平衡區(qū)寬度可由式(9)得到[26]：

式中，x0為應(yīng)力極限平衡區(qū)寬度;λ為側(cè)壓系數(shù);φ0為煤層界面的內(nèi)摩擦角;K為最大垂直應(yīng)力集中系數(shù);γ為上覆巖層體積力;H為埋深;c0為煤層界面黏聚力;Px為支架對煤幫的支護(hù)強度。

計算式(9)根據(jù)表3得到的參數(shù)取值，其中煤層與頂?shù)装褰唤缑娴牧W(xué)參數(shù)c0和φ0低于煤體的c和φ值，c0和φ0分別取0.7 MPa和26°。由此計算得到采高3 m時應(yīng)力極限平衡區(qū)寬度x0為12.2 m，采高6 m時x0為24.4 m。上述數(shù)值模擬仿真結(jié)果與理論計算基本吻合。

上述工作面前后及側(cè)向支承壓力分布及演化規(guī)律，與采空區(qū)垮落帶漸進(jìn)壓實承載力學(xué)特性、裂隙帶覆巖運移特征、采空區(qū)及煤柱形態(tài)特征具有耦合關(guān)系。由此表明，采用研究得到的垮落帶力學(xué)特性及其雙屈服模型反演參數(shù)、裂隙帶巖體劣化模型，實現(xiàn)了采動應(yīng)力場與采空區(qū)覆巖特性的耦合，采動應(yīng)力演化的仿真模擬是可靠的，該方法對區(qū)段煤柱穩(wěn)定性、回采巷道采動影響與控制設(shè)計、巖層與地表移動、保護(hù)層開采等工程問題，具有重要的科學(xué)意義及應(yīng)用價值。

5 結(jié) 論

(1)在綜合分析覆巖“兩帶”高度計算及影響因素的基礎(chǔ)上，研究得到垮落帶漸進(jìn)壓實承載過程的力學(xué)特性表達(dá)式，垮落帶巖體壓縮剛度隨壓縮應(yīng)變呈現(xiàn)指數(shù)型增長，碎脹系數(shù)和單軸抗壓強度對壓縮剛度及應(yīng)變有顯著影響，闡述了垮落帶模型力學(xué)特性參數(shù)的反演方法，正確反映出采空區(qū)垮落帶壓實—支撐過程的力學(xué)響應(yīng)。

(2)基于采空區(qū)垮落帶壓實承載過程的力學(xué)特性與裂隙帶巖體的破斷劣化特性，采用反演力學(xué)參數(shù)后的雙屈服模型和二次開發(fā)得到的工程巖體劣化模型，等效描述垮落帶和裂隙帶的力學(xué)特性，提出了采動應(yīng)力場與覆巖兩帶耦合方法，實現(xiàn)了采動應(yīng)力場與采空區(qū)覆巖“兩帶”特征耦合分析。

(3)通過趙固二礦11030工作面應(yīng)用實例分析，采用采動應(yīng)力場與覆巖兩帶耦合方法，研究了工作面采動應(yīng)力場演化規(guī)律，結(jié)果表明采高對超前支承壓力分布、采空區(qū)應(yīng)力恢復(fù)和側(cè)向煤柱集中應(yīng)力分布均有顯著影響。將計算結(jié)果與理論分析進(jìn)行對比，驗證了采動應(yīng)力場與覆巖兩帶耦合方法的可行性和正確性，提高了采動應(yīng)力場和圍巖采動響應(yīng)研究的嚴(yán)謹(jǐn)性和可靠性，對區(qū)段煤柱穩(wěn)定性、回采巷道采動影響與控制設(shè)計、巖層與地表移動、保護(hù)層開采等工程問題，具有重要的科學(xué)意義及應(yīng)用價值。

[1] 謝和平,周宏偉,王金安,等.FLAC 在煤礦開采沉陷預(yù)測中的應(yīng)用及對比分析[J].巖石力學(xué)與工程學(xué)報,1999,18(4):397-401. XIE Heping,ZHOU Hongwei,WANG Jin’an,et al.Application of FLAC to predict ground surface displacements due to coal extraction and its comparative analysis[J].Chinese Journal of Rock Mechanics and Engineering,1999,18(4):397-401

[2] BADR S A.Numerical analysis of coal yield pillars at deep longwall mines[D].Colorado:Colorado School of Mines,2004.

[3] 浦海，繆協(xié)興.綜放采場覆巖冒落與圍巖支承壓力動態(tài)分布規(guī)律的數(shù)值模擬[J].巖石力學(xué)與工程學(xué)報,2004,23(7):1122-1126. PU Hai,MIAO Xiexing.Numerical simulation of dynamic falling of overlying rocks and abutment pressure in surrounding rocks for fully-mechanized top-coal caving stope[J].Chinese Journal of Rock Mechanics and Engineering,1999,23(7):1122-1126.

[4] PAPPAS D M,MARK C.Behavior of simulated long-wall gob material[R].Washington:United States Department of the Interior Bureau of Mines,1993.

[5] 張振南,繆協(xié)興,葛修潤.松散巖塊壓實破碎規(guī)律的試驗研究[J].巖石力學(xué)與工程學(xué)報,2005,24(3):451-455. ZHANG Zhennan,MIAO Xiexing,GE Xiurun.Testing study on compaction breakage of loose rock blocks[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(3):451-455.

[6] 張振南,茅獻(xiàn)彪,郭廣禮.松散巖塊壓實變形模量的試驗研究[J].巖石力學(xué)與工程學(xué)報,2003,22(4):578-581. ZHANG Zhennan,MAO Xianbiao,GUO Guangli.Experimental study on deformational modulus of friable rock during compactionp[J].Chinese Journal of Rock Mechanics and Engineering,2003,22(4):578-581.

[7] 馬占國.采空區(qū)破碎巖體壓實和滲流特性研究[M].徐州:中國礦業(yè)大學(xué)出版社,2009.

[8] YAVUZ H.An estimation method for cover pressure re-establishment distance and pressure distribution in the goaf of longwall coal mines[J].International Journal of Rock Mechanics and Mining Sciences,2004,41(2):193-205.

[9] 白慶升,屠世浩,袁永，等.基于采空區(qū)壓實理論的采動響應(yīng)反演[J].中國礦業(yè)大學(xué)學(xué)報,2013,42(3):355-361,369. BAI Qingsheng,TU Shihao,YUAN Yong,et al.Back analysis of mining induced responses on the basis of goaf compaction theory[J].Journal of China University of Mining & Technology,2013,42(3):355-361,369.

[10] 錢鳴高,石平五.礦山壓力與巖層控制[M].徐州:中國礦業(yè)大學(xué)出版社,2003.

[11] PALCHIK V.Influence of physical characteristics of weak rock mass on height of caved zone over abandoned subsurface coal mines[J].Environmental Geology,2002,42(1):92-101.

[12] BAI M,KENDORSKI F S,VAN ROOSENDAAL D J.Chinese and North American high-extraction underground coal mining strata behavior and water protection experience and guidelines[A].Proceedings of the 14th International Conference on Ground Control in Mining[C].West Virginia Univ.,Morgantown,WV (United States),1995.

[13] 金太,席京德,蔣金泉,等.緩傾斜厚煤層礦山壓力與巖層控制[M].徐州:中國礦業(yè)大學(xué)出版社,2000.

[14] SALAMON M.Mechanism of caving in longwall coal mining[A].Rock mechanics contribution and challenges[C].Proceedings of the 31st US symposium of rock mechanics,Golden,Colorado,1990:161-168.

[15] WHITTLES D N,LOWNDES I S,KINGMAN S W,et al.Influence of geotechnical factors on gas flow experienced in a UK longwall coal mine panel[J].International Journal of Rock Mechanics and Mining Sciences,2006,43(3):369-387.

[16] JIANG Lishuai,ZHANG Peipeng,CHEN lianjun,et al.Numerical approach for goaf-side entry layout and yield pillar design in fractured ground conditions[J].Rock Mechanics and Rock Engineering,DOI:10.1007/s00603-017-1277-0.

[17] XUE Junhua,WANG Hanpeng,ZHOU Wei,et al.Experimental research on overlying strata movement and fracture evolution in pillarless stress-relief mining[J].International Journal of Coal Science & Technology,2015,2(1):38-45.

[18] HOEK E,CARRANZA-TORRES C,CORKUM B.Hoek-Brown failure criterion-2002 edition[J].Proceedings of NARMS-Tac,2002,1:267-273.

[19] MITRI H S,EDRISSI R,HENNING J G.Finite-element modeling of cable-bolted stopes in hard-rock underground mines[J].Transactions-Society For Mining Metallurgy and Exploration Incorporated,1995,298:1897-1902.

[20] JIANG Lishuai,SAINOKI Atsushi,MITRI Hani S,et al.Influence of fracture-induced weakening on coal mine gateroad stability[J].International Journal of Rock Mechanics and Mining Sciences,2016,88:307-317.

[21] 蔣力帥.工程巖體劣化與大采高沿空巷道圍巖控制原理研究[D].北京:中國礦業(yè)大學(xué)(北京),2016.

[22] YAN S,BAI J,WANG X,et al.An innovative approach for gateroad layout in highly gassy longwall top coal caving[J].International Journal of Rock Mechanics and Mining Sciences,2013,59:33-41.

[23] LI W,BAI J,PENG S,et al.Numerical modeling for yield pillar design:a case study[J].Rock Mechanics and Rock Engineering,2015,48(1):305-318.

[24] HU Zhenqi,CHEN Chao,WU Xiao,et al.Surface movement and deformation characteristics due to high-intensive coal mining in the windy and sandy region[J].International Journal of Coal Science & Technology,2016,3(3):339-348.

[25] HU Haifeng,LIAN Xugang.Subsidence rules of underground coal mines for different soil layer thickness:Lu’an Coal Base as an example,China[J].International Journal of Coal Science & Technology,2015,2(3):178-185.

[26] 侯朝炯.巷道圍巖控制[M].徐州：中國礦業(yè)大學(xué)出版社，2013：186-194.

Numericalsimulationoncouplingmethodbetweenmining-inducedstressandgoafcompression

JIANG Lishuai1,2，WU Quansen2，LI Xiaoyu2，DING Nan2

(1.StateKeyLaboratoryofMiningDisasterPreventionandControlCo-foundedbyShandongProvinceandtheMinistryofScienceandTechnology,ShandongUniversityofScienceandTechnology,Qingdao266590,China; 2.CollegeofMiningandSafetyEngineering,ShandongUniversityofScienceandTechnology,Qingdao266590,China)

The goaf stress recovery under the effect of overlying strata subsidence has notable effect on mining-induced stress field.Based on the coupling effect between ground stress and goaf compression,the mechanical characteristics and properties back-analysis method of caved rock are studied,the tension-weakening characteristics of rock mass in fractured zone are analyzed,and hereby the coupling simulation method between mining-induced stress and goaf compression is proposed and verified using a simulation case study.The results show that the mechanical behavior of the compressed caved zone follows an exponential growth,the compression stiffness nonlinearly increases with the compression strain.The mechanical properties of caved zone obtained from back-analysis agree with the analytical results,which contribute to the simulation of the mechanical behavior of goaf.The tension-weakening model,which is developed from strain-softening model by adding the tension-weakening algorithm,can describe the weakening of Young’s modulus,hereby achieve to describe the mechanical behavior of fractured rock mass rigorously.Based on the mechanical behavior simulation of caved and fractured zone,the proposed method can conduct the coupling numerical simulation of mining-induced stress field and the overlying “two-zone”.

mining-induced stress;caved zone compression;fractured zone weakening;numerical simulation;coupling analysis

10.13225/j.cnki.jccs.2016.1717

TD323;TD823

:A

:0253-9993(2017)08-1951-09

國家自然科學(xué)基金資助項目(51704182);山東省自然科學(xué)基金資助項目(ZR2017BEE050);山東省高等學(xué)?？萍加媱澷Y助項目(J17KA212)

蔣力帥(1989—)，男，江蘇南通人，講師。E-mail:jlsh1989@126.com

蔣力帥，武泉森，李小裕，等.采動應(yīng)力與采空區(qū)壓實承載耦合分析方法研究[J].煤炭學(xué)報,2017,42(8):1951-1959.

JIANG Lishuai，WU Quansen，LI Xiaoyu,et al.Numerical simulation on coupling method between mining-induced stress and goaf compression[J].Journal of China Coal Society,2017,42(8):1951-1959.doi:10.13225/j.cnki.jccs.2016.1717