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

1.吉林大學建設(shè)工程學院, 長春130026；2.青海省環(huán)境地質(zhì)重點實驗室, 西寧810007；3.青海省環(huán)境地質(zhì)勘查局, 西寧810007；4.青海省地質(zhì)環(huán)境保護與災(zāi)害防治工程技術(shù)研究中心, 西寧810007

0 引言

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

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

1 地質(zhì)概況及模型建立

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

1.1 物理模型建立

在幾何相似的多個系統(tǒng)中，進行相同性質(zhì)的物理過程，當幾何上的對應(yīng)點在相應(yīng)時間點上的表現(xiàn)具有一定的一致性時，稱這樣的物理過程為相似現(xiàn)象。具有相似現(xiàn)象的物理系統(tǒng)在幾何形態(tài)上具有相似性，遵循著同樣的物理定律。

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

幾何相似系數(shù)：

(1)

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

應(yīng)力相似系數(shù)。

(2)

式中：σp為原型應(yīng)力；σm為模型應(yīng)力。

容重相似系數(shù)：

(3)

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

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

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

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

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

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

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

1.2 數(shù)值模型建立

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

表2 數(shù)值模型各層參數(shù)Table 2 Layer parameters of numerical model

圖2 模型建立及網(wǎng)格剖分Fig.2 Modeling and meshing

2 模擬結(jié)果及分析

2.1 物理模型模擬結(jié)果及分析

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

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

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

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

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

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

2.2 數(shù)值模型模擬結(jié)果及分析

2.2.1 摩爾庫倫本構(gòu)模型下數(shù)值模擬分析

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

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

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

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

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

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

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

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

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

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

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

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

圖8 第五年I-I剖面剪應(yīng)力以及垂向應(yīng)力分布圖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 理論公式驗證

為進一步驗證數(shù)值模擬的準確性，將模擬的圍巖壓力結(jié)果與理論公式中的普式壓力拱理論計算結(jié)果進行橫向?qū)Ρ?。普式理論認為，松散巖體中作用于硐頂?shù)膰鷰r壓力由自重應(yīng)力提供，工程中為了方便計算，往往將硐頂?shù)淖畲髧鷰r壓力簡化為均布荷載而不計硐軸線的變化引起的圍巖應(yīng)力變化。因此，硐頂?shù)淖畲髧鷰r壓力公式：

(4)

側(cè)向壓力公式：

(5)

式中：f為堅固性系數(shù)；b、b1為拱的矢高；a為側(cè)壁穩(wěn)定時平衡拱的跨度；a1為自然平衡拱的最大跨度。

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

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

4 結(jié)論

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

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

(3)煤礦開采導(dǎo)致地表的大范圍沉陷，嚴重影響研究區(qū)人類活動和建筑物安全。采空區(qū)上方，豎直方向上，地表中央的豎向位移大，兩側(cè)的豎向位移?。凰椒较蛏?，采空區(qū)的影響范圍略大于煤層開采的范圍。