大采深綜放開采覆巖移動規律離散元數值模擬研究

余學義，王昭舜，楊 云，毛旭魏

( 1.西安科技大學 能源學院，陜西 西安 710054；2.西安科技大學 教育部西部礦井開采與災害防治重點實驗室，陜西 西安 710054 )

實現科學采礦已然成為新時期煤礦開采的必然要求[1-2]。近年來，隨著煤礦機械化程度的不斷提高，不少煤礦已采用大采深綜放開采的生產方式，高強度的資源開發必然帶來開采沉陷等破壞。故而大采深綜放開采地表移動變形規律的研究就顯得尤為重要，而此項研究的首要條件就是了解上覆巖層的移動變形規律。

早在20世紀初，便有專家通過簡單的力學原理解釋了一些礦山壓力現象，并初步形成了礦山壓力假設。隨著壓力拱假說，懸臂假說和鉸接巖石假說的發展，礦山壓力理論已經形成了一個系統。在19世紀80年代，錢鳴高[3-4]等提出控制全部或部分上覆地層的地層稱為關鍵地層；茅獻彪[5]等研究了采場覆蓋層中關鍵層的破壞規律和關鍵層負荷的變化規律；許家林[6-12]等提出了基于關鍵層控制的局部充填開采技術以及影響深部開采關鍵層對地表沉陷影響的一些因素；黃慶享[13-15]等研究了厚沙層地質條件下關鍵層上部的載荷分布，傳遞規律以及破斷規律；余學義[16-18]等對大采高覆巖災變規律進行了研究，分析了影響下沉的采礦因素，研究得出大采深綜放開采地表移動變形規律；李化敏[19]等研究了特厚煤層綜放開采中的放煤理論。

眾多學者經過長期研究，還針對近距離厚煤層組工作面覆巖破壞規律、特厚煤層不同開采階段“弱-弱”結構覆巖破壞高度、厚煤層綜放開采支架圍巖關系、頂板運動、窄煤柱圍巖控制及含夾矸特厚煤層綜放工作面頂煤破碎機理等開展了大量研究[20-21]。

到目前為止，眾多學者已經進行了大量研究工作，但是大采深綜放開采條件下，隨著工作面開采數量的增加，地表從極不充分-非充分-充分采動過程中覆巖的演化規律對地表移動變形的影響鮮有深入研究。本文基于陳家溝煤礦特殊的地質條件，研究厚埋深大采高綜放多工作面開采過程中的覆巖演化規律以及對地表移動變形的影響。

1 試驗礦井地質條件

通過多個鉆孔探測，采區地層自上而下有表土層的第四系全新統，總厚度為268.75 m的上第三系甘肅群，涵蓋直羅組與延安組的中侏羅統，下侏羅統富縣組與上三疊統延安群。詳情如圖1，表1所示。

陳家溝煤礦主采5煤層，工作面平均寬120 m，煤層傾角6°，分層開采高度12 m。平均開采深度為506 m，最大開采深度達到580 m，即大采深( 煤層埋深500～1 000 m )煤礦[22]。

圖1 巖層柱狀圖Fig.1 Histogram of rock formation

表1 陳家溝煤礦井田地層Table 1 Field layer profiles of Chenjiagou Mine

根據鉆孔取芯分析，按巖性分類，上覆巖為中硬巖性，且泥巖厚度占較大比例。由于泥巖本身具有很好的防水性，其已成為覆巖中良好的隔水層。

2 地表移動觀測結果分析

工作面測點布置如圖2所示。8512工作面累計回采長度2 235 m，走向觀測線共計35個有效工作測點，點號分別為Z1～Z35，傾向觀測線共計35個工作測點，點號分別為Q1～Q35；8513工作面累計推進1 710 m，走向觀測線共計118個工作測點，點號分別為zx1～zx118，傾向觀測線共計40個工作測點，點號分別為qx1～qx40。

8512工作面自2013-11—2015-11觀測結束共計觀測15次。8513工作面從2016-05-16開始，截止到2019-01-13，累計觀測23次。8512工作面走向最大下沉達到1.124 m，8513工作面走向最大下沉達到2.226 m，由于2個工作面走向均屬超充分采動，故將傾向地表下沉作為研究重點，現對2個工作面傾向的觀測數據進行分析。

圖2 工作面測點布置Fig.2 Layout of measuring points in working face

8512工作面傾向開采寬度120 m，開采寬深比小于1/4，屬于極不充分采動，根據多次觀測數據，繪制傾向觀測點下沉曲線如圖3所示，未出現下沉盆地。8513工作面根據觀測數據，繪制觀測點下沉曲線如圖4所示。傾向累計工作面寬度為246 m，屬于非充分采動。

圖3 8512工作面傾向下沉曲線Fig.3 Sinking curves of 8512 working face

圖4 8513工作面傾向下沉曲線Fig.4 Sinking curves of 8513 working face

通過分析觀測數據可知，在8512首采工作面開采后，下沉變化過程呈規律性增大且下沉量均不大。根據8513工作面觀測結果發現，在qx24點、qx27點處下沉發生突變，而其正處于8512工作面的采空區上覆巖層，說明8513工作面的開采對8512工作面采空區產生擾動，導致8512工作面上覆巖層發生活化，造成下沉量增加。在地表發育過程中由于采深較大，當回采結束后，采空區內圍巖應力處于新的平衡，處于穩定階段；當再次開采，擾動范圍內的工作面會打破這種應力平衡，已達到平衡的采空區會被活化，從而導致二次沉降，這與觀測數據吻合，故而在大采深開采過程中，地表達到穩態的周期會較長。

由此得出：8512工作面開采后由于是單一工作面開采，此時地表處于極不充分采動狀態，受關鍵層支承上部壓力的影響，傳遞到地表表現出沉降量較小，下沉不明顯，沉降影響范圍也較小。在8513工作面開采后，地表下沉逐漸增大，且部分區域產生擾動，此時地表處于非充分采動狀態。故地表由極不充分采動向非充分采動發育時，下沉會有突變，且下沉速度突變，即是關鍵層斷裂失穩所致。

此外，利用地表移動盆地的主斷面上實測下沉曲線求得地表移動角參數。具體結果詳見表2。

表2 實測角參數Table 2 Measured angular parameters

3 覆巖結構控制條件分析

3.1 關鍵層分析

該礦區煤層覆蓋層的綜合巖性為中硬型。較硬的地層主要分布在侏羅紀地層。上第三系甘肅群地層一般較軟，結構疏松，抗壓強度低，遇水易開裂和風化；第四系地層薄而軟；特厚煤層中煤5層的頂部為砂巖，屬于較厚的硬巖地層；在中侏羅統第二和第三層中，有幾層較厚且堅硬的巖性，這些地層在覆蓋層的運動和變形中起著控制作用。使用固定支撐梁的力學模型來分析關鍵層的應力狀態。

設其撓曲方程為

則滿足邊界條件：

利用伽遼金法通過解算，可得固支梁的最大撓度為

式中，I為慣性矩；E為巖體的彈性模量，GPa；q為巖體上的載荷，kN/m2。

3.2 關鍵層的確定

上覆巖層的荷載q系由巖梁自重及其上覆巖層共同作用，因此計算各巖層承載的荷載可按式( 4 )計算，詳細參數見表3。

表3 鉆孔部分巖層力學參數Table 3 Mechanical parameters of some rock strata in borehole

式中，E1，E2，…，En為各層巖層的彈性模量，MPa；h1，h2，…，hn為各層巖層的厚度，m；γ1，γ2，…，γn為各層巖層的容重，kN/m2；( qn)1為第1層巖層所受上方第2層至第n層的載荷，kPa。

計算結果見表4。

表4 上覆巖層承載載荷Table 4 Bearing load of overlying strata

計算得出堅硬巖層為第3層礫巖、第6層砂質泥巖、第11層細粒砂巖和第19層粉砂巖，并計算破斷距綜合得出第19層為關鍵層1，第11層為關鍵層2，關鍵層3為第3層礫巖。

計算結果表明，8512工作面開采后關鍵層2以下的巖層均斷裂破壞，關鍵層2起著控制覆巖沉陷的作用，但是關鍵層2以上巖層在移動過程中還會出現離層結構，到達地表的沉陷量小于關鍵層2的撓曲量。

4 覆巖移動演化規律模擬

4.1 模擬試驗數據分析

數值模擬采用3DEC塊體離散元數值模擬軟件，模擬基于陳家溝煤礦八采區地質采礦條件的多工作面開采覆巖移動變形規律。

結合地質資料和巖石力學試驗結果，得出各巖層物理力學參數并采用莫爾-庫侖塑性模型來進行計算,得到

式中，K為體積模量；E為楊氏模量；G為剪切模量；v為泊松比。

本模型巖層節理面材料模型采用庫侖滑動模型( 完全彈塑性 )，各節理面物理力學參數見表5。

表5 各節理面力學參數Table 5 Mechanical parameters of each joint surface

續 表

4.2 試驗模型建立

為了研究陳家溝煤礦大采深綜放開采覆巖移動變形規律，采用3DEC三維計算模型，模型具有6個自由面，水平方向和底部邊界采取薄層的方法進行速度為零的固定。模型尺寸為1 000 m×30 m×570 m，本文模擬單一煤層開采，煤層埋深525 m，煤層為水平煤層，煤層厚度26 m，設計4個開采工作面，分別為8512，8513，8511，8514工作面，工作面寬度均為120 m，留設區段煤柱6 m，各單元在運動時一般不發生彈跳，這是由于運動時動能轉化成熱能耗散掉的緣故，因此煤巖塊的運動是不可逆的過程，為避免巖塊在平衡位置振動，應加入阻尼來吸收系統的動能，本文采用了自適應阻尼。三維數值計算模型如圖5所示。

圖5 3DEC三維計算模型Fig.5 3D computing model of 3DEC

4.3 模擬結果分析

4.3.1 851 2工作面模擬結果分析

在模擬過程中，在4個工作面的覆巖和地表運動穩定后，可以獲得覆巖移動變形圖、地表沉降云圖和測量線數據。通過覆巖移動變形圖和云圖，可以獲得覆巖和地表沉降的分布，地表的最大沉降點，地表沉降的范圍和裂隙破裂帶的范圍；地表沉降曲線可以通過對地表測量線數據的分析得出。

關鍵層1初次破斷以及周期破斷時塊體的位置變化，如圖6所示，隨著8512工作面的推進，垮落高度呈上升趨勢，開采至10 m時，垮落高度6 m；開采至20，30 m 時，垮落高度分別增加到14.50 m 和23.40 m。工作面持續推進，直接頂繼續垮落，由于巖層間的互相擠壓，當開采至70 m時，關鍵層1初次破斷，此時監測到地表下沉0.54 m；繼續開采至27.10 m，失穩現象嚴重，第1次周期來壓，由圖6可知，下部巖層產生了約5 m的超前斷裂，此時監測到地表下沉0.77 m。

圖6 關鍵層1首次破斷Fig.6 First break of the key stratum 1

模擬8512工作面開采結束，如圖7所示，開采尺寸為120 m，平均采高9 m，工作面寬深比為0.24，此時地表處于極不充分采動狀態。覆巖中關鍵層1發生破斷，破斷距分別為70.0 m和27.1 m。從煤層頂板上方向上33.4 m為覆巖垮落帶，再向上76.9 m的范圍為裂縫帶，即裂縫帶發育高度為110.3 m，裂采比為12.2，裂縫帶向上直達地表范圍為覆巖的彎曲下沉帶。由于開采范圍小，開采對覆巖擾動破壞影響到關鍵層2下部，關鍵層2以下的巖層均斷裂破壞，關鍵層2表現為微量的撓曲變形，但是關鍵層2以上巖層在移動過程中還會出現離層結構，其上覆巖土層表現為彎曲下沉。根據模擬觀測結果表明在開采一個工作面后，覆巖移動變形規律與關鍵層2撓曲形態吻合，再次證明，地表處于極不充分采動狀態是由于關鍵層起控制作用。

圖7 8512工作面開采結束Fig.7 End of 8512 working face mining

4.3.2 851 3工作面模擬結果分析

8513工作面在開采過程中，關鍵層2初次破斷以及周期破斷時塊體的位置變化如圖8所示，工作面在持續推進過程中，上覆巖層向上發育迅速。當開采至146.6 m時，關鍵層2底部初次破斷，整體表現為局部破斷的撓曲變形，此時地表下沉1.18 m；繼續推進70.2 m，累計推進213.8 m時，關鍵層2周期破斷，此時監測地表下沉1.22 m。工作面初次來壓是70.0 m，周期來壓是21.7 m，這與現場實測及理論計算結果均基本吻合。

圖8 關鍵層2首次破斷Fig.8 First break of the key stratum 2

模擬8513工作面開采結束，如圖9所示，開采尺寸為240 m，平均采高12 m，工作面寬深比為0.48，地表處于非充分采動狀態。覆巖移動破壞向上傳遞致使關鍵層2中下部發生破斷，關鍵層2初次破斷距與周期破斷距分別為146.6 m和70.2 m。覆巖垮落帶高度增至54.8 m，再向上100 m的范圍為裂縫帶，即裂縫帶發育高度為154.8 m，其裂采比為12.9，關鍵層整體表現為局部破斷的撓曲變形。裂縫帶向上379.1 m直達地表為彎曲下沉帶范圍。由于8512，8513工作面之間區段煤柱的支承作用，在覆巖一定高度范圍內，破斷的巖層呈現“波浪形”分布，由于關鍵層2在自重和上覆巖層載荷作用下產生變形，與下部巖層運動相反，因此在地表處于非充分采動狀態下，關鍵層2與其下部巖層之間產生了離層，關鍵層2略有失穩。此處亦證實了在地表移動觀察過程中下沉值發生突變的原因。

圖9 851 3工作面開采結束Fig.9 End of 8513 working face mining

4.3.3 模擬與實測數據分析

( 1 ) 模擬與實測數據對比分析

8512 工作面開采后，實測地表最大下沉為1.124 m，數值模擬結果顯示地表最大下沉為1.156 m，下沉系數為0.128，與實測實時下沉率0.124基本吻合。8513工作面開采結束后，由于關鍵層2的控制作用，下沉發育至地表后下沉量2.2 m，如圖10所示，與現場實測的數據2.226 m基本吻合。模擬下沉量誤差70 mm左右。

圖10 模擬開采與實測數據對比Fig.10 Comparison of simulated mining and measured data

( 2 ) 模擬與實測數據回歸分析

統計學中的假設檢驗和回歸分析法可以研究2個變量之間是否存在相關關系及數量規律，并使用數學公式的方法來表達變量之間的關系[23]。本文應用該方法分析，模擬開采2個工作面后得到的下沉數據與該工作面在現場實測的數據，驗證其是否存在相關性，并檢驗模擬開采數據對實測數據有無顯著性影響。

對觀測值與模擬值進行回歸分析，回歸模型見表6，并由方差分析法分析2組數據，判斷系數R2值均接近最大值1，表明2變量間的線性關系密切，2組變量即為正線性相關，如圖11所示；F值均小于F臨界值，即在顯著水平0.05下，數值模擬數據對實測下沉數據無顯著性影響，說明模擬開采下沉數據與地表觀測下沉數據基本吻合。也表明在數值模擬中參數的選取是合理的。

表6 模擬結果與實測結果回歸模型匯總Table 6 Summary of simulation and measured regression models

圖11 模擬數據與實測數據回歸曲線Fig.11 Regression curves of simulated data and measured data

由前2個工作面模擬開采數據與實測值的對比分析可知，數值模擬在參數的選取以及模擬地質情況上是合理的、可取的，故而進行第3個工作面及第4個工作面開采，從而獲取覆巖及地表移動變形的相應參數是必要的、科學的。

4.3.4 多工作面模擬結果分析

模擬8511工作面開采完成，如圖12所示，開采規模達到360 m，開采高度為12 m，工作面寬深比為0.72，地表仍處于非充分采動狀態。覆巖中關鍵層2完全斷裂，斷裂距離為70.2 m，關鍵層2上方有2層超厚的泥巖，會抑制導水裂縫帶發育，破裂不會遍及整個巖體，而是裂縫僅發展到一定高度。覆巖垮落帶的高度為54.8 m，由于厚砂質泥巖的抑制作用，且表現為同步彎曲，故砂質泥巖以上為彎曲下沉帶，導水裂縫帶的發育高度為196 m( 開采高度的16.3倍 )，彎曲下沉帶為向上336.7 m。關鍵層2和上覆的泥巖相隔較大。關鍵層2破裂后，由于自重的原因逐漸壓實，該斷面煤柱的支護影響范圍大于8513工作面開采結束時的支護影響范圍，最大地面沉降為3.8 m。

圖12 8511工作面開采結束Fig.12 End of 8511 working face mining

圖13 8514工作面開采結束Fig.13 End of 8514 working face mining

模擬8514工作面開采完畢，如圖13所示，開采規模達到480 m，開采高度為12 m，工作面寬深比為0.96，此時，采區已被充分開采，基本達到充分采動狀態。“三帶”發育高度與8511工作面開采完成后基本一致，只是下沉盆地范圍增大。覆巖中關鍵層2與其下伏巖層之間的離層閉合，這是由于：一方面，開采尺寸的繼續增大，關鍵層的鉸接結構范圍在橫向擴展的同時也伴隨著整體撓曲下沉；另一方面，關鍵層下伏巖層破斷后持續壓實壓密到一定程度后，不再發生繼續下沉，隨著工作面的推進，覆巖發育轉向橫向發育，破壞范圍擴大，下沉盆地也隨著擴大。當關鍵層的鉸接結構撓曲下沉量達到下伏巖層的最大允許下沉值時，離層不再存在。此時，在充分采動狀態下，關鍵層2及其覆巖的壓力通過巖層傳遞至采空區，所有這些都由該段煤柱承擔，煤柱應力達到強度極限從而產生塑性變形，失去了煤柱的支承作用，完全被煤屑充填密實后，關鍵層3對地表起絕對控制作用，上覆巖層的運動基本達到了穩定狀態，模擬結果表明最大地面沉降為5.7 m。

如圖14可知，隨著開采工作面的增加，除下沉值下沉曲線逐步顯現下沉盆地，邊界角、移動角也隨著開采工作面的增加而減小，下沉盆地在擴大。下沉曲線由不對稱分布逐漸變為對稱分布，最大下沉點在采空區中央。

圖14 模擬開采后各工作面下沉曲線Fig.14 Simulated subsidence curves of each working face after mining

綜上所述：分析4個工作面開采模擬結果可知，隨著開采工作面數量的增加，開采寬度增加，寬深比增大，上覆巖層也在逐步發生變化，8512工作面開采后，寬深比為0.24，屬于極不充分采動狀態，由于關鍵層的控制作用，地表沉陷小；8513工作面開采后，寬深比為0.48，地表從極不充分采動變為非充分采動，由于關鍵層失穩，地表急劇下沉，覆巖內部結構失穩；開采8511工作面后，寬深比為0.72，地表仍處于非充分采動狀態，沉降值仍在增加，關鍵層完全破裂，由于厚泥巖的抑制作用，裂縫只發展到一定深度；8514工作面開采后，寬深比為0.96，關鍵層2的鉸接結構范圍已經完全密實不再沉降，覆巖的發育轉向橫向發展，破壞范圍擴大，從而沉降盆地擴大，此時，地表基本達到充分采動狀態，關鍵層3未產生破斷，對地表起絕對控制作用，最大下沉值達到5.7 m。由8512與8513工作面對比來看，不同采高情況下，下沉量的變化較大，覆巖破壞程度較大，垮落帶與導水裂縫帶在采高12 m的條件下發育較大，比采高9 m條件下均有所增加。采高越大，覆巖垮落帶高度增加，導水裂縫帶向上發育尺寸基本不變。

4.4 對比分析

從開采至今已完成2個工作面的開采任務，地表狀態也由最初的極不充分采動發育至如今的非充分采動，未來地表將達到充分采動狀態。故而現將地表發育過程中模擬值與實測值進行對比說明。

由表2，7可知，地表在達到穩態后，地表下沉值達到最大，極不充分狀態下由于覆巖關鍵層控制，地表下沉量較小；非充分采動狀態下覆巖關鍵層破斷失穩，地表下沉速度加速，地表產生較寬的臺階裂縫；充分采動狀態下，覆巖由縱向發育變為橫向發育，地表表現為開采沉陷區擴大，導水裂縫帶范圍增大。模擬得出達到充分采動時，下沉系數為0.475，最大下沉值為5.7 m；由地表移動角量參數可知，地表在發育過程中，隨著下沉盆地的擴大，移動角及邊界角均有約3°的減小，主要影響范圍角也相應有較大幅度的減小，這表明下沉盆地范圍擴大較明顯，影響范圍擴大；地表由極不充分采動向非充分采動發育過程中，角度變動較大，由非充分采動向充分采動發育過程中，角度變化較小，最大下沉角逐步接近90°。

表7 數值模擬角參數Table 7 Numerical simulation of angular parameters

5 結 論

( 1 ) 由地表移動觀測數據分析，8512工作面開采后地表處于極不充分采動，8513工作面開采后地表處于非充分采動，地表由極不充分采動向非充分采動發育時，8512工作面采空區活化，產生二次沉降，造成下沉量突變，且下沉速度突變，即由關鍵層斷裂失穩所致。

( 2 ) 應用關鍵層理論結合3DEC數值模擬方法，模擬多工作面開采關鍵層斷裂失穩與地表移動的關系。結果表明：工作面寬深比為0.24時，地表極不充分采動狀態下由于覆巖中關鍵層控制垮落，地表表現為沉降較小，裂縫寬度較窄，地表破壞范圍較小；工作面寬深比為0.48時，地表處于非充分采動狀態下覆巖關鍵層由離層—失穩—破斷，地表也相應的表現為沉降突變、裂縫寬度增加并產生臺階裂縫且破壞范圍增加；工作面寬深比為0.72時，裂采比為16.3倍采高，地表仍處于非充分采動狀態，工作面寬深比為0.96時，基本達到充分采動狀態下，覆巖由縱向發育變為橫向發育，地表表現為開采沉陷區擴大，主要影響范圍角減小，影響范圍增大。

( 3 ) 模擬結果得出地表達到充分采動時，下沉系數為0.475，最大下沉值為5.7 m，大采深綜放開采采高越大，垮落帶高度越大，導水裂縫帶向上發育尺寸基本不變，同時應用統計學回歸分析與方差分析方法將地表移動觀測數據與數值模擬數據做假設檢驗，確定了模擬參數的可靠性。為該礦實現科學開采提供了理論依據。