煤層開采后覆巖穩定性物理模型及數值模擬

張謙，張延軍，于子望，張通，張雨，羅銀飛

1.吉林大學建設工程學院, 長春130026；2.青海省環境地質重點實驗室, 西寧810007；3.青海省環境地質勘查局, 西寧810007；4.青海省地質環境保護與災害防治工程技術研究中心, 西寧810007

0 引言

煤炭是中國重要能源之一，與國民經濟的發展有著密切聯系，隨著中國經濟的不斷發展，煤炭作為主要能源其消耗量也在與日俱增[1,2]。煤層開采過程中引發的上覆巖層移動及地表沉陷等問題嚴重威脅著人民的生活和生產。因此對煤層開采后上覆巖層變化規律的研究具有重大意義。現階段，對上覆巖層變化規律的研究方法主要為相似材料模擬和數值模擬。由于巖土介質本構的復雜性及結構與巖土介質之間的相互作用，相似材料模擬的模擬結果相對貼近實際，一直是研究覆巖穩定性的重要手段[3--6]，但相似材料模擬存在著實驗費用高，場地約束大等問題[7--9]。采用數值模擬具有可以靈活地改變模型的形狀及參數，受客觀條件影響小，節約經費且能處理復雜的材料本構關系等優點，技術上與相似材料模擬優勢互補，結果上相互印證。

隨著FLAC軟件的不斷成熟，近年來其已經成為眾多學者研究地表沉降規律的一種重要手段。李培現等[10]用FLAC3D對老采空區的地基穩定性進行了研究。郭俊廷等[11]用 FLAC3D對非充分采動下地表移動變形值進行分析，并獲取了地表的移動參數。李想等[12]用 FLAC3D軟件求解地表移動帶范圍，優化了圈定移動帶的方法。周振亮等[13]用FLAC求解出條帶開采對地表變形值的影響。何標慶[14]采用FLAC軟件對大型采空區群的穩定性進行模擬。一些學者[15，16]將FLAC模擬地表移動的結果與理論公式計算結果及INSAR監測數據進行對比，進一步證明了FLAC模擬的準確性。筆者通過建立物理模型和采用FLAC3D軟件進行數值模擬兩種方法，將鄂爾多斯某煤礦的物理模型開挖結果與數值模擬結果對比分析，得出該煤礦下部煤層開采對上覆巖層穩定性的影響范圍及規律。

1 地質概況及模型建立

研究區位于鄂爾多斯高原東北部，地貌總體上波狀起伏，河流不甚發育，屬于典型的構造剝蝕高原地貌。地勢開闊呈波狀起伏，枝狀溝谷發育其中；南部地勢平坦，多被風積砂覆蓋，并有洼地、湖盆散列其中。

1.1 物理模型建立

在幾何相似的多個系統中，進行相同性質的物理過程，當幾何上的對應點在相應時間點上的表現具有一定的一致性時，稱這樣的物理過程為相似現象。具有相似現象的物理系統在幾何形態上具有相似性，遵循著同樣的物理定律。

相似系數是原型與模型之間相同物理量之比，通常用C表示。在煤層相似材料模擬中常用的相似系數有：

幾何相似系數：

(1)

式中：δp為原型位移；δm為模型位移；為Lp原型長度；Lm為模型長度。

應力相似系數。

(2)

式中：σp為原型應力；σm為模型應力。

容重相似系數：

(3)

式中：λp為原型容重；λm為模型容重。

為準確建立研究區模型，收集了研究區內的工程地質、環境地質、水文地質資料，根據研究區內的鉆孔資料，以及參考前人的資料[17]，最終確定物理模型的相似常數為：CL=200、Cλ=1.5、Cσ=300。原地層平均容重為λ1=26.32 KN/m3，因此模型的相似配比材料容重設計為λ2=17 KN/m3。由于原型橫斷面高290 m，最終模型整體尺寸設計為長×高×寬=1.2 m×1.45 m×0.12 m。

為簡化模型鋪設，將部分巖層簡化為均布荷載鋪設于模型的上部，由于上覆巖層的巖性變化引起下部巖層荷載分布的動態與非均布的特點，模型的上邊界加載遵循的原則：巖性均一的覆巖可簡化為均布荷載加載于模型的上邊界，覆巖中有突變巖層時，只有突變巖層的上部巖層可以簡化為均布荷載。按照此原則最終將原型簡化為十層模型。模型由上至下分布如表1。

根據相似材料性質、模型特點、實際地質情況及前人對采空區模型實驗的經驗[18]，物理模型采用石英砂、細砂做骨料，石灰、石膏作為膠結材料，由下至上逐層建立模型，通過控制膠結材料的使用量控制巖層強度，以每次2～3 cm的厚度進行模型裝填，由于石膏具有遇水極易凝結硬化的特點，加入濃度為1.0%的硼砂溶液作為緩凝劑，并在每層之間加入云母粉模擬實際巖層之間的原生裂隙，模型制作完成后進行一周的晾置，得到最終模型。

表1 相似模型地層分布Table 1 Stratigraphic distribution of similar model

模型建立完成后，采用在頂層表面布置千分表的方式對實驗數據進行監測。在模型最上方非等距布設千分表，靠近模型中軸處，千分表的布設間距較小，隨著遠離中軸，千分表布設距離增大。開采方法采用水平方向壁式采煤法，開采實際寬度200 m(模擬開采寬度100 cm)的煤層。

圖1 物理模型圖Fig.1 Physical model diagram

1.2 數值模型建立

通過研究區地質數據統計資料，確定地層標高及厚度；參考鉆孔資料確定煤層埋深，最終建立地質概念模型。通過工作區內鉆孔物理、力學測量結果，結合實驗室力學實驗最終確定數值模型參數。為方便模型建立與計算，將物理性質相似的非關鍵層簡化為同一層，采用Midas描述并建立研究區數值模型(圖2)，參數設置如表2。

表2 數值模型各層參數Table 2 Layer parameters of numerical model

圖2 模型建立及網格剖分Fig.2 Modeling and meshing

2 模擬結果及分析

2.1 物理模型模擬結果及分析

實際生產中采用壁式開挖法對3--1煤層進行開采，開挖煤層實際厚度為4 m(模型厚度2 cm)，實際開采寬度200 m(模型寬度100 cm)，實際開采時間為5 a。實驗中，由于模型開挖層厚度較小(2 cm)，因此采用小刀進行開挖，將距模型右邊界20 cm處設為右端點，在右端點處以10 cm/h的速度由右向左對3--1煤層進行開挖，一次開挖時長為2 h。每次開挖完成后，每隔1 h測量并記錄冒落帶和裂隙帶高度，記錄千分表監測數據，驗證模型是否達到穩定，若模型達到穩定，則進行下一次開挖。

圖3 物理模型開挖效果圖Fig.3 Physical model excavation effect diagram

實驗過程中，第三次開挖，即開挖長度為60 cm時，可以看到模型出現裂隙，在第三次開挖過程中伴隨部分坍塌與脫落現象的出現(圖3)。從采動巖體移動和離層裂隙發育情況可以看出，隨著工作面的推進，巖體裂隙自上而下逐步發展。在采動作用下，距離開挖層較近的層面首先開裂，隨著開采范圍的增大，距離開挖層較遠的巖層也開始產生裂隙。隨著工作面的進一步推進，在前一工作面推進距離條件下的采動巖體裂隙網格上又疊加了后一工作面推進距離條件下的采動巖體裂隙網絡，從而使采動巖體裂隙分布更加復雜化，不僅有新的裂隙產生，而且原有裂隙也發生了一系列的擴展、閉合以及張開。

實驗發現，隨著開挖的進行，模型的裂隙發育率逐漸增大，開挖0.5 m時，模型的裂隙發育率為23.11%，當開挖至1.0 m時，模型的裂隙發育率為30.41%。經測量，裂隙帶高度為32.5 cm，約為采高的16.25倍。煤層采出后，冒落帶僅出現在開采區域，冒落塊體以外的巖體不受影響。當工作面開采至1.0 m時(實際200 m)，覆巖冒落高度為7.5 cm(實際15 m)，是采出煤層厚度的3.75倍。冒落帶的左側垮落角36°，右側垮落角57°。實驗開采完成時，裂隙并未延伸到含水砂層，但是考慮地層原有裂隙的存在，3--1煤層開挖還是會對上方的含水砂層產生較大的影響。

圖4 冒落帶、裂隙帶、底板位移圖Fig.4 Displacement diagram of caving zone, fissure belt and floor slab

根據地表千分表制作底板位移圖，采用網格測量法，制作冒落帶、裂隙帶圖。由圖4可以看出，底板的最大位移近似等于采高，最大位移為2 cm，裂隙帶位移遠小于底板位移，最大位移為0.9 cm，冒落帶最大位移為0.51 cm。以上3個層均具中間位移大，兩端位移小的特征，近似呈以采空區中心線為對稱軸的類似“V”字型。

2.2 數值模型模擬結果及分析

2.2.1 摩爾庫倫本構模型下數值模擬分析

數值模擬的主要目的是計算煤層開挖對上覆巖層造成的影響，首先應消除重力作用對上覆巖層的影響。因此，在開挖前，需要對初始應力進行求解，將初始應力解出的位移和應力值清零后再進行開挖，以此確定開挖造成的地表形變。

達到初始應力過程中的最大不平衡力的變化趨勢如圖5所示，從圖5中可以看出，最大不平衡力先增加后逐漸變小，最后趨于穩定，最終最大不平衡力數值<1.0×10-5，說明模型已經達到初始平衡。

圖5 最大不平衡力變化圖Fig.5 Variation diagram of maximum unbalanced force

根據研究區勘察報告中的采礦區開采情況，本次數值模擬采用走向長壁式開挖法對3--1煤層進行開挖，設計年采煤量288萬t，平均采高2.02 m，工作面長度250 m，年進尺400 m，開挖走向由南向北開挖，計算五年內采空區對上覆巖層穩定性的影響。

采用FLAC3D對采空區開挖模擬的豎向位移如圖6所示。由圖6可以看出，煤層開采后，煤層上方巖層逐漸懸空，地下形成采空區，此時頂板巖層可能會出現冒落、坍塌的現象，遠離采空區，這種現象也逐漸變弱。巖層的彎曲、破環現象一般沿著巖層的層理法線方向發生，伴隨著這些現象發生的還有部分離層現象與斷裂的產生。隨著開采的進行，地層的表面也開始出現這些現象，此時地表形成下沉盆地。由圖6可以看出，開挖對地表的影響范圍要大于采空區實際尺寸，地表下沉形狀近似橢圓狀。并且這種現象隨著時間的推移，變得愈發明顯，在第五年時，地表沉降及位移達到了最大值。地表最大沉降量也隨時間推移逐漸增大，第一年地表最大沉降量為0.093 m，第二年為0.301 m,第三年為0.500 m，第四年為0.676 m，第五年為0.833 m。

由于自重應力的存在，在采煤過程中，巖層的初始平衡狀態遭到破壞，應力重新調整、分布，最終達到新的平衡狀態。FLAC3D應力分布云圖通過顏色深淺顯示出工作區的各個區域位移值的大小關系，越趨近于紅色的區域，應力集中現象越明顯。由剪應力分布圖與垂向應力分布圖中可以看出，在開采的第一年，上覆巖層中出現應力集中現象，隨著煤層開采的推進，在采空區及覆巖所在的區域均有剪應力形成，且伴隨著開采面積的增大，剪應力的范圍與應力集中現象也逐漸地增大，在第五年，最大剪應力已經達到3.63×106Pa。隨著上覆巖層應力集中現象愈發明顯，上覆巖層逐漸破壞，當應力增加到一定程度，上覆巖層出現裂隙，導致巖層失去隔水性，最終發生破壞。

對比五年塑性區變化可以看出，當開采面積增大，應力隨之增大，塑性區的范圍也逐年增加。從垂向應力可以看出，由于開挖，地下應力重分布，采空區及上覆巖層周圍出現應力集中現象，在采空區兩側頂點處這種現象尤為明顯。由應力狀態可以推出，在第五年時，由于開挖作用，上覆巖層底板在重力及覆巖的共同作用下，在法線方向產生彎曲與移動，變形達到一定程度后，頂板內部應力劇增，應力增大到超過巖石的極限抗拉強度時，巖石發生破壞，塑性區出現，上覆巖層也隨之出現斷裂、破碎。

2.2.2 修正劍橋模型下3--1煤層開挖模擬

采用修正劍橋本構模型，對研究區開展數值模擬。建立模型的尺寸、地層分層及開挖方式均與摩爾庫倫本構模型相同，同樣選擇開挖3--1煤層，工作面長度為250 m，年進尺400 m，開挖走向為由南向北，其模擬結果如圖7、8所示。

從圖9可以看出，采用修正劍橋模型和庫倫模型模擬得到的位移分布形態基本相似，但沉降量模擬結果上，二者差異較大，修正劍橋模型模擬至第五年時，地表最大沉降量為0.377 m。出現這種結果的原因可能是修正劍橋模型主要適用于土層的模擬而摩爾庫倫模型主要適用于巖層的模擬。

圖6 3--1煤層開采位移分布云圖Fig.6 Displacement distribution cloud map of 3--1 coal seam

圖7 第五年I-I剖面塑性區分布剖面圖Fig.7 Distribution profile of plastic zone in I-I profile in the fifth year

圖8 第五年I-I剖面剪應力以及垂向應力分布圖Fig.8 Shear stress and vertical stress distribution map of I-I section in the fifth year

圖9 修正劍橋模型下開挖3--1煤層第五年位移分布云圖Fig.9 Correction of the fifth year displacement distribution of 3--1 coal seam under Cambridge model

3 理論公式驗證

為進一步驗證數值模擬的準確性，將模擬的圍巖壓力結果與理論公式中的普式壓力拱理論計算結果進行橫向對比。普式理論認為，松散巖體中作用于硐頂的圍巖壓力由自重應力提供，工程中為了方便計算，往往將硐頂的最大圍巖壓力簡化為均布荷載而不計硐軸線的變化引起的圍巖應力變化。因此，硐頂的最大圍巖壓力公式：

(4)

側向壓力公式：

(5)

式中：f為堅固性系數；b、b1為拱的矢高；a為側壁穩定時平衡拱的跨度；a1為自然平衡拱的最大跨度。

圖10 煤層開挖中軸線上覆巖層垂向最大應力對比圖Fig.10 Vertical maximum stress diagram of overburden strata on axis of coal seam excavation

對比普式壓力拱理論計算結果與摩爾庫倫模型數值模擬計算結果發現，兩者計算結果中最大豎向應力的變化數量級為105Pa。兩者的模擬曲線整體擬合效果較好，對于形態上的差距，這可能是由于普式壓力拱理論在計算等速開采時，計算結果只受開采寬度的影響，因此其計算的硐頂最大豎向壓力表現為線性增加。而數值模擬計算的結果由于考慮了材料的彈塑性變形以及圍巖摩擦力的影響，其計算的最大豎向應力曲線則表現出一定的緩變性。兩者計算結果的對比進一步證明了采用摩爾庫倫模型進行數值模擬的結果更加準確。

4 結論

(1)在鄂爾多斯某煤礦，隨著煤層開采面積的增大，上覆巖層會產生較大的剪應力，塑性區(達到力學強度)的范圍也會逐年增加。

(2)修正劍橋模型和庫倫摩爾模型第五年的地表沉降量模擬結果差距較大，對比理論計算解，可見庫倫摩爾模型在模擬煤層開挖對上覆巖層的影響方面具有優勢。

(3)煤礦開采導致地表的大范圍沉陷，嚴重影響研究區人類活動和建筑物安全。采空區上方，豎直方向上，地表中央的豎向位移大，兩側的豎向位移小；水平方向上，采空區的影響范圍略大于煤層開采的范圍。