劉 輝，左建宇，蘇麗娟，程 樺，朱曉峻，張鵬飛，王保國

(1.安徽大學 資源與環(huán)境工程學院，安徽 合肥 230601；2.安徽省礦山生態(tài)修復工程實驗室，安徽 合肥 230601；3.安徽大學 數學科學學院，安徽 合肥 230601；4.蚌埠市勘測設計研究院，安徽 蚌埠 233040)

0 引 言

儲量豐富的煤炭資源作為我國最重要的基礎能源，一直占據我國能源消費總量一半以上，擔負著為經濟發(fā)展提供主動力的責任[1-4]。我國煤炭資源賦存狀態(tài)具有一定的地域特色，東部礦區(qū)較西部礦區(qū)普遍具有煤層埋藏深、松散層厚且潛水位高等地質結構特點，其復雜的地質結構下開采地表移動變形較常規(guī)地質條件更為劇烈，主要表現為：下沉量增大、下沉系數接近甚至超過1.0、下沉盆地影響范圍擴大、沉陷后易形成積水區(qū)。因此對礦區(qū)地表生態(tài)環(huán)境和群眾居住生活產生了嚴重影響[5-8]。

針對巨厚松散層條件下開采的特殊地表移動變形規(guī)律，國內外學者從現場分析、理論分析和實驗模擬等多角度對其機理進行了大量研究。研究表明巨厚松散層下開采地表移動變形的特殊性不僅僅是由上覆巖層綜合巖性較軟所導致的，其不可忽略的因素是煤層開采引起上覆巖層內含水層疏水滲流固結所帶來的附加沉降[9-12]。對于上覆巖層內含水層疏水滲流固結引起的附加沉降，我國學者實驗模擬證明了含水層失水固結機理，量化分析了含水層固結沉降量，建立了松散層底部含水層失水固結模型，取得了一定的成果[13-18]。但是以往的研究手段主要為理論分析和相似材料模擬，難以全面考慮到應力場和滲流場的耦合作用對巨厚含水松散層下開采地表移動變形的影響，從而導致結果的偏差。

以菏澤礦區(qū)某礦典型的巨厚含水松散層地質采礦條件為原型，在分析工作面實測數據的基礎上，采用FLAC3D數值模擬軟件，充分考慮應力場和滲流場耦合作用的影響，對巨厚含水松散層下煤層開采地表移動變形進行了流固耦合模擬，總結了采動過程中滲流固結沉降動態(tài)規(guī)律，對比分析了有無含水層情況下煤層開采地表移動變形差異，研究了松散層內含水層位置對地表移動變形的影響，為更準確預測及控制巨厚含水松散層下煤層開采引起的地表移動變形，提供一定的理論依據和技術參考。

1 研究區(qū)概況

菏澤礦區(qū)某礦位于山東省菏澤地區(qū)鄆城縣城南約10 km處，處于巨野煤田的中北部，其范圍東起田橋斷層及田橋支斷層，西至煤系地層底界露頭，南起3925000緯線，北至25勘探線，南北長約14 km，東西寬約11 km，地表總面積達到了222 km2。礦區(qū)地面標高+41.60～+45.38 m，淺層地下水埋深2～7 m，上覆新近系和第四系松散層平均厚度達到了590 m，局部超過700 m，屬于典型的巨厚含水松散層地質條件，礦區(qū)綜合柱狀如圖1所示。

圖1 礦區(qū)綜合柱狀Fig.1 Stratigraphic column of mining area

通過對鉆孔柱狀圖以及水文地質報告分析可將上覆巨厚松散層概括劃分為5層含水、隔水層，見表1。巨厚松散層內第四系含水的砂層與隔水的黏土層相間沉積，屬于富水性中等的松散孔隙含水層，新近系含水的砂層單層呈犬牙交錯狀相連，屬于富水性強的松散孔隙承壓含水層。

表1 礦區(qū)含水、隔水層劃分Table 1 Division of aquifers and aquiclude in mining area

某礦1308工作面走向長度約630 m，傾向長度約230 m，主采煤層為山西組3煤層，煤層平均厚度為3 m，平均傾角為3°，為近水平煤層，煤層采用長壁垮落法開采，采深為765 m，工作面上覆松散層平均厚度為582 m，煤層頂板主要由粗砂巖、中細砂巖、砂質泥巖等組成，具有開采深度大、傾角小、松散層厚、松散層內含有多層含水層的特點。1308工作面移動盆地主斷面上方采用十字線布設了地表移動觀測線，走向觀測線全長約1 550 m，累計布設47個觀測點，傾向觀測線全長約1 800 m，累計布設53個觀測點，采動期間以平均35 d一次的頻率進行觀測，至停采后共觀測12次。根據地表移動觀測數據繪制了地表動態(tài)下沉曲線(圖2)，地表移動變形參數如下：

工作面1308下沉系數q1.092水平移動系數b0.25主要影響角正切tan β1.60邊界角δ056.90°

圖2 工作面走向線動態(tài)下沉曲線Fig.2 Dynamic sink curve of strike line of working face

可以看出，巨厚含水松散層下開采地表移動變形與常規(guī)地質條件下采動存在一定的特殊性。相比較之下，巨厚含水松散層下開采具有移動變形初始期時間短、起動距短、初始期和活躍期下沉速度快、下沉系數大于1、主要影響角正切和邊界角小、下沉盆地影響范圍大、下沉盆地邊界處水平位移大于下沉量且下沉盆地衰退期長、邊緣收斂緩慢的特點。

2 模型設計與方案

2.1 模型設計

為研究巨厚含水松散層下開采地表移動變形規(guī)律，以菏澤礦區(qū)某礦1308工作面地質采礦條件為原型，采用FLAC3D數值模擬軟件建立了模型尺寸為3 000 m×1 500 m×824 m(長×寬×高)的三維數值模擬模型。模型單元尺寸在水平方向為25 m，垂直方向按照巖層厚度不同進行合理設置。模型采高為3 m，煤層設計為水平煤層，設計的工作面長度為1 200 m，寬度為400 m，工作面沿走向方向分為12次開采，開挖步長均為100 m，開采深度為765 m，其上覆巖層內松散層厚度為582 m。為避免模型過小所導致的邊界應力效應，在工作面走向和傾向兩側各留900 m和550 m的邊界保護煤柱，三維數值模擬模型如圖3所示。

圖3 數值模擬模型Fig.3 Numerical simulation model

模型的位移邊界條件設置為：約束模型底面各個方向的位移，模型頂面為自由面，前后左右四面約束x、y水平方向的位移，但可以發(fā)生z方向的移動。滲流模型設置為各巖層均為各向同性且均勻等效的連續(xù)孔隙介質，力學模型采用Mohr-Coulomb屈服準則，數值模擬巖層參數根據地質資料進行概化合并見表2。

表2 數值模擬巖層參數Table 2 Numerical simulation of rock parameters

2.2 模擬方案

為研究巨厚含水松散層下開采地表移動變形規(guī)律，明確巨厚含水松散層在下沉盆地形成中的作用以及含水層位置對地表移動變形的影響，分別制定了動態(tài)固結沉降、無含水層、含水層位置3種數值模擬方案。

1)動態(tài)固結沉降方案。根據礦區(qū)含水、隔水層劃分表(表1)在松散層內設置含水、隔水層，進行流固耦合模擬，研究巨厚含水松散層下采動過程中滲流固結沉降動態(tài)規(guī)律。

2)無含水層方案。不考慮巨厚含水松散層在采動中產生的流固耦合影響，進行純力學開挖模擬，通過對比分析，研究有無含水層情況下煤層開采地表移動變形差異。

3)含水層位置方案。在松散層內設置厚度為40 m的單層含水層，通過調整含水層位置標高，分別模擬含水層底部與基巖的距離為0、100、200、300、400 m的5種情況，研究巨厚松散層內含水層位置標高對地表移動變形的影響。

3 巨厚含水松散層對下沉盆地影響

為明確巨厚含水松散層內的多層含水層在煤層開采下沉盆地形成過程中的作用及影響，對動態(tài)固結沉降方案模擬結果進行歸納總結分析，并與無含水層方案模擬結果進行對比分析，研究了巨厚含水松散層滲流固結沉降動態(tài)規(guī)律，并對比分析了有無含水層情況下煤層開采地表移動變形差異。

3.1 巨厚含水松散層滲流固結沉降動態(tài)規(guī)律

動態(tài)固結沉降方案通過循環(huán)開關流體滲流模塊從而實現流固耦合模擬：在煤層開挖后，先關閉FLAC3D中的流體滲流分析模塊，迭代計算至平衡，得到模型在單力學場中的土體不排水變形量(開挖變形沉降量)，然后開啟流體滲流模塊，流固耦合計算土體在該開挖階段內的排水變形量(滲流固結沉降量)，流固耦合計算完成后再進行下一步開挖，并如此往復循環(huán)計算，直至工作面開挖結束。

根據FLAC3D流固耦合模擬機制，對動態(tài)固結沉降方案模擬計算結果進行分析，繪制了隨工作面推進受開挖變形和滲流固結耦合作用影響的地表最大下沉發(fā)展曲線圖，如圖4所示。從圖4可看出，當工作面推進至1 200 m開挖結束時，地表最大下沉為1.988 m，其中開挖變形沉降在地表總沉降中占主導地位，由開挖變形引起的沉降量為1.820 m，而滲流固結引起的沉降量為0.168 m，結合地表總沉降量可得出巨厚含水松散層滲流固結引起的沉降量占最終地表總沉降的8.5%。

圖4 地表最大下沉發(fā)展曲線Fig.4 Maximum surface subsidence development curve

為研究巨厚含水松散層滲流固結沉降動態(tài)規(guī)律，通過計算得到了每一步開挖過程中開挖變形沉降增量和滲流固結沉降增量，并繪制了隨工作面推進的地表沉降增量變化曲線，如圖5所示。通過對地表沉降增量變化曲線進行分析可將開挖變形和滲流固結沉降增量的變化隨工作面推進劃分為同步增長期、動態(tài)變化期、同步減緩期3個階段。

圖5 地表沉降增量變化曲線Fig.5 Incremental change curve of surface subsidence

1)同步增長期。在工作面開采初期，工作面推進至200 m階段內，開挖變形和滲流固結沉降增量呈現出同步增長的趨勢，地表開始出現明顯下沉，起動距約為1/4煤層埋深。

2)動態(tài)變化期。工作面推進至200～600 m，開挖變形和滲流固結沉降增量在整體上共同出現明顯的增大和減小波動，其中滲流固結沉降增量隨工作面推進的變化波動更為劇烈，與此同時地表出現劇烈下沉變形。

3)同步減緩期。隨著工作面的繼續(xù)推進至終采線1 200 m處，逐步接近走向充分采動時，開挖變形和滲流固結沉降增量呈現出同步減緩的趨勢，但其整體上仍存在較小的變化波動。

3.2 地表移動變形差異對比分析

為更好地研究巨厚含水松散層的多層含水層對下沉盆地的影響，對比分析了動態(tài)固結沉降方案與無含水層方案的模擬結果，并繪制了工作面推進至1 200 m時動態(tài)固結沉降方案與無含水層方案的工作面走向下沉曲線，如圖6所示。下沉曲線整體關于采空區(qū)中心對稱且曲線保持連續(xù)漸變，由于工作面已推進至1 200 m，走向達到充分采動，下沉曲線底部呈現平底狀，兩者均符合厚松散層條件下采動引起的地表下沉曲線分布，但仍存在一定的差異性。

圖6 工作面走向下沉曲線Fig.6 Strike subsidence curve of working face

通過對模擬結果進一步對比分析，得到了地表最大變形，見表3。由表3可看出：考慮松散層內含水層影響進行流固耦合計算后，動態(tài)固結沉降方案相較于無含水層方案，最大下沉值增大了0.220 m、最大傾斜增大了0.178 mm/m、最大水平移動增大了0.063 m，而最大曲率和最大水平變形分別減小了0.003 mm/m2和0.587 mm/m。

表3 地表最大變形Table 3 Numerical simulation of rock mechanics parameters

表4給出了動態(tài)固結沉降方案與無含水層方案的地表移動變形參數，通過對比分析可得：

表4 地表移動變形參數對比Table 4 Comparison table of surface movement deformation parameters

1)考慮松散層內含水層影響進行流固耦合計算后，動態(tài)固結沉降方案的下沉系數達到了1.096，相較于無含水層方案的0.974，下沉系數增大了12.5%；相比較之下，水平移動系數變化極小，動態(tài)固結沉降方案的水平移動系數僅比無含水層方案減小了1.0%。

2)對2種情況的工作面主斷面走向線下沉數據分析計算可得出，動態(tài)固結沉降方案和無含水層方案的主要影響半徑分別為677、641 m，相比較之下考慮含水層影響進行流固耦合計算后，主要影響半徑增大了5.62%。因2種模擬計算方案的工作面平均采深不變，動態(tài)固結沉降方案的主要影響角正切比無含水層方案減小了5.28%。

3)動態(tài)固結沉降方案和無含水層方案的邊界角分別為48.86°和50.22°，相比較之下動態(tài)固結沉降方案考慮含水層影響進行流固耦合計算后，邊界角減小了2.71%。

綜上所述，巨厚松散層內含水層受采動影響導致松散層內孔隙水壓力消散，有效應力增大，產生滲流固結沉降，引發(fā)下沉盆地在開挖變形沉降的基礎上再平衡[19-20]。從而導致地表下沉和水平移動變形增大，主要影響半徑增大，主要影響角正切和邊界角減小，下沉盆地的地表移動影響范圍擴大，但最大曲率和最大水平變形的減小說明下沉盆地的整體形狀在巨厚松散層內含水層疏水滲流固結沉降的作用下更為平緩[21-22]。

4 地表變形參數與含水層位置的關系

上文中對巨厚含水松散層內的多層含水層在下沉盆地形成過程中所產生的影響進行了研究，得知巨厚松散層內含水層受采動影響產生的疏水滲流固結現象會導致地表移動變形增大、地表移動影響范圍擴大。為進一步探索巨厚松散層內不同位置標高的含水層對地表移動變形的影響，通過對含水層位置方案模擬數據結果進行分析，研究了地表移動變形參數與巨厚松散層內含水層位置之間的關系。

為方便表述含水層位置方案中所設置的含水層在松散層內所處的位置，定義含水層位置s為含水層底部距基巖的距離與松散層厚度的比值。

4.1 下沉系數變化規(guī)律

根據模擬結果，繪制了下沉系數q與含水層位置s之間的關系，如圖7所示。

圖7 下沉系數和水平移動系數與含水層位置關系Fig.7 Correlation curve of subsidence factor and displacement factor with aquifer location

從圖7中可得：含水層底部與基巖之間的距離為0時，下沉系數為1.098。隨著含水層底部與基巖之間距離的增加，下沉系數呈現出線性增長趨勢，當含水層底部與基巖之間的距離增加至400 m時，下沉系數達到最大值1.158。對含水層底部與基巖間距為0～400 m的5種模擬結果的下沉系數與含水層位置進行相關性分析，得出下沉系數q與含水層位置s之間的相關性公式為：

q=0.087 7s+1.101 2

(1)

綜上可得，下沉系數與松散層內含水層位置之間存在正線性關系。

4.2 水平移動系數變化規(guī)律

根據模擬結果，繪制了水平移動系數b與含水層位置s之間的關系圖，如圖7所示。

從圖7中可得：5種模擬結果的水平移動系數為0.287～0.300，數值波動較小。當含水層底部與基巖之間的距離為0時，水平移動系數為0.295，隨著含水層與基巖之間距離的增加，水平移動系數呈現出先增大后減小的變化規(guī)律，在含水層底部與基巖間距為200 m時達到0.300的最大值。對含水層底部與基巖間距為0～400 m的5種模擬結果的水平移動系數與含水層位置進行相關性分析，得出水平移動系數b與含水層位置s之間的相關性公式為：

b=-0.073 1s2+0.040 5s+0.293 8

(2)

綜上可得，水平移動系數與松散層內含水層位置之間存在先增大后減小的二次函數關系。

4.3 主要影響角正切變化規(guī)律

根據模擬結果，繪制了主要影響角正切tanβ與含水層位置s之間的關系，如圖8所示。

圖8 主要影響角正切和邊界角與含水層位置關系Fig.8 Correlation curve of tangent of major influence tan β and limit angle δ0 with aquifer location

由圖8可得：含水層底部與基巖之間的距離為0時，主要影響角正切為1.115。隨著含水層底部與基巖之間距離的增加，主要影響角正切呈現出先增大后減小的變化規(guī)律，在含水層底部與基巖間距為300 m時達到1.145的最大值。對含水層底部與基巖間距為0～400 m的5種模擬結果的主要影響角正切與含水層位置進行相關性分析，得出主要影響角正切tanβ與含水層位置s之間的相關性公式為：

tanβ=-0.238 3s2+0.165 4s+1.112 3

(3)

綜上可得，主要影響角正切與松散層內含水層位置之間存在先增大后減小的二次函數關系。

4.4 邊界角變化規(guī)律

根據模擬結果，繪制了邊界角δ0與含水層位置s之間的關系圖，如圖8所示。

從圖8可得：含水層底部與基巖之間的距離為0時，邊界角為63.77°。當含水層位置與基巖不再相鄰時，邊界角急劇減小至58.54°，并隨著含水層與基巖之間距離的增加，呈現出明顯的線性增長趨勢，增大至59.37°。對含水層底部與基巖間距為100～400 m的4種模擬結果的邊界角與含水層位置進行相關性分析，得出邊界角δ0與含水層位置s之間的相關性公式為：

δ0=1.600 2s+58.322

(4)

綜上可得，邊界角與松散層內含水層位置之間存在先減小后增大的關系。

5 結 論

1)某礦巨厚含水松散層采礦條件下受采動引起的滲流固結沉降占地表總下沉量的8.5%；隨工作面推進，開挖變形和滲流固結沉降增量的變化可劃分為同步增長期、動態(tài)變化期、同步減緩期3個階段。

2)巨厚含水松散層采礦條件下，松散層內含水層受采動影響產生的疏水滲流固結現象會導致地表移動變形增大、地表移動影響范圍擴大，但最大曲率和最大水平變形的減小說明下沉盆地的整體形狀更為平緩。

3)地表變形參數與巨厚松散層內含水層位置之間存在相關性：隨著松散層內含水層位置的上移，下沉系數線性增大，水平移動系數和主要影響角正切先增大后減小，而邊界角先減小后增大。