壩下深部煤層開采壩體移動變形規律的數值模擬

張 峰，秦洪巖，題正義

（1.山西工程技術學院 采礦工程系，山西 陽泉 045000；2.遼寧工程技術大學 礦業學院，遼寧 阜新 123000；3.華北科技學院，北京 101601）

庫壩下深部煤層開采，既要保證煤層開采的安全性，也要保證壩體的安全運營，掌握煤層開采過程中壩體的移動變形機理，才能最大限度地對庫壩下煤炭資源進行合理開發利用[1]。但煤層埋藏深度大、地表庫壩土體結構小，煤層埋藏深度與壩體高度比值較大，常規數值模擬、相似材料模擬和理論計算等方法難以推算壩體的移動變形規律。采用數值模擬方法，按照邊長5 m 至少要劃分300 多萬個單元網格，根據求解時間與網格數大致呈N4/3的正比關系[2-3]，完成一次模型的開挖計算大約需要6.8 d，完成整個模擬需要的時間之長是非常驚人的。利用相似材料按照1∶200 的相似比進行模擬計算，建立至壩體的全尺寸模型過大，現有裝置滿足不了[4-6]。依據煤層上覆巖層的礦壓顯現規律推斷壩體的移動變形，由于煤層埋藏深度過大，頂板巖層結構復雜，推算結果也難以保證其準確度[7-8]。為此，提出將地表和壩體結構的移動變形分開研究，定量地分析地表移動變形下壩體結構的移動變形，獲取壩體結構隨工作面不同推進長度下的移動變形規律，建立壩下深部煤層開采壩體結構局部放大的時空序列反演機制，能夠為庫壩下煤層開采壩體移動變形機理的分析和壩體的安全性評價提供理論指導依據。

1 分析方法的構建

1.1 方法的提出

壩下煤層開采后，采空區覆巖的移動變形自下而上，首先波及到壩體周圍的地表，然后帶動壩體一起移動變形。將此傳遞過程分解為2 個部分[9-12]，一是工作面開采引起的覆巖及地表的移動變形，得到隨工作面推進形成不同尺寸的下沉盆地；二是逐漸擴大的下沉盆地帶動壩體的移動變形，得到隨下沉盆地不斷擴大后壩體不同部位的移動變形。如果能夠舍去工作面開采和下沉盆地形成這一過程，按照下沉盆地形成的過程直接開挖盆地，從而觀測分析壩體的移動變形和破壞，將會大大縮小模型尺寸、減少運算時間和提高小目標的分析精度，壩下開采覆巖及地表移動變形過程分解如圖1。

1.2 模型構建及實施步驟

由圖1 可以看出，如果直接在地表進行開挖，壩體的移動變形會立刻發生相應的移動變形，甚至破壞。這與深部工作面開采引起的緩慢、連續變形不相符的，為了使模擬能夠再現煤層開采引起的壩體變形破壞，在地表以下保留一定厚度的緩沖巖土層，根據采空區覆巖變性破壞的“三帶”結構，留設的巖土層厚度至少要大于垮落帶的高度，才能保證地表不切冒。

圖1 壩下開采覆巖及地表移動變形過程分解示意圖Fig.1 Decomposition diagrams of mining overburden and surface movement and deformation process under the dam

用于挖掘的盆地選擇由數值模擬或地表移動變形預計得出的下沉盆地。開挖步距可分成2 部分：一是在壩體下沉之前設置較大的步距；二是在壩體下沉時設置較小步距。簡化模型的開挖順序如圖2。

圖2 簡化模型的開挖順序Fig.2 Excavation sequence of simplified model

簡化模擬分析的實施步驟如下：

1）依據概率積分法和數值模擬計算出煤層開采后地表的最終下沉盆地，選擇最大盆地作為開挖盆地的大小。模型的長寬高依據下沉盆地的影響范圍進行確定。

2）按照計算結果中不同推進長度時地表下沉曲線影響范圍確定開挖范圍，對模型中的下部巖土進行開挖，隨時觀察地表的移動變形。

3）針對模擬結果中不同推進距離下的地表下沉曲線，逐次進行開挖、觀測、記錄壩體移動變形規律，直至壩體移動變形穩定。

4）根據各階段壩體的移動變形狀態，總結分析隨著工作面推進、地表下沉變化時壩體的移動變形規律。

2 數值模擬

根據局部放大分析方法的內涵，結合地表移動變形預計結果和數值模擬得到的工作面推進過程中地表的下沉盆地，通過淺部開挖下沉盆地，模擬壩體受采動影響下的變形破壞規律，檢驗分析方法的可行性。

2.1 模擬方案

依據地表移動變形預計影響壩體走向的最大范圍為570 m，確定模型的尺寸為600 m×2 300 m×87 m。壩體下覆留設30 m 厚的巖土層，中部留設開挖體，開挖體下覆留設40 m 左右巖層用于消除底板對開挖體的影響。

模型開挖尺寸依據工作面與壩體中心位置之間的距離劃分2 種：一種是當間距大于200 m 時，此時開采還沒有對壩體產生影響，為節約運算時間，選擇工作面每推進100 m 時地表移動變形范圍作為開挖步長；另一種是間距在200 m 以內時，為深入觀測壩體移動變形，選擇工作面每推進5 m 時地表移動變形范圍作為開挖步長。

利用概率積分法預計壩下工作面開采后地表任意點的下沉值，將任意點的下沉值對應的三維空間坐標導入到Midas 建模軟件中[13-14]，繪制出開采壩下工作面地表移動變形的下沉盆地，按照確定尺寸將下沉盆地從模型中分割出來，在開挖體中分割出工作面方向和推進方向均為曲面的開挖步距，然后按照工作面的推進方向對開挖體依次進行分割，劃分網格。巖石力學基礎參數見表1。

表1 巖石力學基礎參數Table 1 Basic parameters of rock mechanics

在大壩走向方向的壩體中心線上布置114 個觀測點，測點間距為5 m，測線長570 m，觀測壩體的移動變形情況。

2.2 壩體變形破壞規律及特征

觀測不同開挖步距下壩體的位移變化、應力變化和塑性區變化[15-16]，分析壩體的移動變形規律，檢驗分析方法的可行性。

2.2.1 位移分析

根據工作面不同推進長度下的計算結果，提取壩體測線中點隨著工作面開采的監測數據，繪制的壩體隨開采變化的下沉曲線如圖3。

圖3 壩體中心位置的下沉曲線Fig.3 Settlement curves at the center of dam body

由圖3 可以看出，工作面推進470 m 時，壩體開始受采動影響，此時壩體中心位置的最大下沉值為0.092 m；推進至1 800 m 時，壩體測線中心位置達到該條件下的最大下沉值7.17 m。

2.2.2 應力分析

由于土體發生變形破壞后其塑性破壞區域無法恢復至原始狀態，為了更好地說明壩體的變形破壞過程，選擇壩體第1 次達到全部塑性破壞時作為重點分析對象，詳細說明壩體的移動變形過程，根據位移變化分析結果，選擇工作面推進至500、600、700、800、850、900、910、920、930、2 000 m 的應力變化云圖進行分析，壩體中心位置的應力變化云圖如圖4（圖中左側為迎水坡，右側為背水坡）。

由圖4 可以看出：

圖4 壩體中心位置的應力變化云圖Fig.4 Diagrams of stress change at the center of dam body

1）在工作面推進500 m 時，壩體迎水坡側干砌石上端出現約1 m 長的壓應力集中區（負為壓應力，正為拉應力），最大值為-0.38 MPa，其余區域的應力值基本保持不變。

2）當工作面推進至600 m 時，壩體背水坡側壩腳處出現應力集中，最大應力值為0.054 MPa。迎水坡側壩肩位置（干砌石段）應力集中區域沿斜坡向下轉移，最大應力值為-0.32 MPa，其它區域的應力值均呈現增大的變化，且為壓應力。

3）當工作面推進至700 m 時，壩體背水坡側壩腳處的應力集中區域依然存在，且應力值由0.054 MPa 增大至0.099 MPa。迎水坡側（干砌石段）應力集中區域基本保持不變，其上覆土體的應力值為正值。其它區域的應力值繼續增大。

4）當工作面推進至800 m 時，背水坡側壩腳處應力集中區域增大，應力值繼續增大（0.378 MPa）。迎水坡（干砌石段）應力集中區域的應力值由-0.316 MPa 增大至-0.693 MPa，應力集中區域也呈現增大的變化趨勢，但迎水坡側壩腳位置（壩基土體）出現了小范圍的應力集中現象，最大應力值為0.01 MPa，其它區域的應力值繼續增大。

5）當工作面推進至850 m 時，背水坡側壩腳處應力集中區域減小，其上覆的應力集中區域減小，總體呈現拉應力區域增大，壓應力區域減小。迎水坡側（干砌石段）應力集中區域增大，應力值由-0.693 MPa 增大至-1.13 MPa，壩腳位置（壩基土體）應力集中區域減小，應力值由0.01 MPa 增大至0.80 MPa，壩頂位置出現小范圍的應力集中區域，最大應力值為-0.87 MPa。其它區域的應力值逐漸減小，均降低至-0.2 MPa 以下。

6）當工作面推進至900 m 時，背水坡側壩腳處應力集中區域增大，其上覆的應力集中區域減小，應力值由拉應力轉化為壓應力。迎水坡側（干砌石段）應力集中區域繼續增大，最大應力值增大至1.7 MPa，壩腳位置的應力集中區域也逐漸增大。其它區域的應力值繼續減小，且為壓應力。

7）當工作面推進至910 m 時，背水坡側壩腳處應力集中區域基本保持不變。迎水坡側（干砌石段）應力集中區域基本不變，處于拉應力狀態，壩腳位置的應力集中區域和應力值均在減小，總體呈現為拉應力狀態。壩頂位置的應力集中區域增大，應力值由-0.87 MPa 增大至1.1 MPa。壩體中心區域出現了壓應力轉化為拉應力變化，最大應力值由-0.95 MPa轉為0.059 MPa。其它區域的應力值繼續減小，且為壓應力。

8）當工作面推進至920 m 時，背水坡側壩腳的應力集中區域消失，應力值降低為-0.021 MPa，上覆的應力集中區域消失，保持在壓應力狀態。迎水坡側（干砌石段）應力集中區域減小，應力值增大至1.9 MPa，壩腳位置的應力集中區域減小，應力值也在減小，轉變為壓應力狀態。壩頂位置的應力集中區域和應力值基本保持不變。壩體中部位置拉應力區域擴展至壩體頂部，且隨著高度的增大，應力值逐漸減小，幾乎為0。其它區域的應力值也逐漸減小，且為壓應力。

9）當工作面從920 m 推進至2 000 m 的過程中，背水坡側、迎水坡側和壩體中部位置土體的應力值呈現1 個緩慢增大的過程，應力值均為負值（壓應力），最大應力值為-0.25 MPa，說明土體在經歷壓縮變形過程。

由上述壩體各部位的應力變化可以得出：壩體經歷了從拉伸變形破壞到壓縮-還原的交替變化過程[17-18]，只是各部位出現拉伸-壓縮轉變過程的時間節點不同，受采動影響開始階段背水坡側壩腳先經歷拉伸變形破壞，其次為迎水坡側壩體，最后為壩體中部位置。經歷由拉伸轉變為壓縮的過程也是由背水坡側壩體開始，其次為迎水坡側壩體，最后為壩體中部位置。

2.2.3 塑性破壞分析

根據壩體應力變化云圖對應的塑性破壞云圖如圖5，由此分析壩體的變形破壞規律及特征。

由圖5 可以看出：

圖5 壩體中心位置的塑性破壞云圖Fig.5 Diagrams of plastic failure at the center of the dam

2）當工作面推進600 m 時，背水坡側壩腳出現小范圍內的拉伸破壞區域（長1.7 m，高1 m），說明此區域內土體受下部采空區側拉伸變形較大（應力值由壓應力轉換為拉應力）。壩體迎水坡側和壩體中部位置均未出現變形破壞。

3）當工作面推進至700 m 時，壩體背水坡側變形破壞區域增大，說明此處壩體受到下部表土層的下沉，對其產生的拉應力增大所致，拉伸破壞區域長2.5 m，高1.3 m。壩體迎水坡和壩體中部位置均未出現變形破壞。

4）當工作面推進至800 m 時，壩體背水坡側產生拉伸破壞的區域繼續增大（長4.8 m，高2.9 m）。迎水坡側壩腳出現小范圍的變形破壞區域（長1.6 m，高1 m），說明此處壩體隨著背水坡側壩腳位置的拉伸破壞，壩體中部位置和迎水坡側壩腳向背水坡側傾斜變形，迎水坡側壩腳上覆沒有限制其移動變形的實體結構，其受到的拉應力大于土體的自身強度，產生拉伸破壞。壩體中部位置未出現變形破壞區域。

5）當工作面推進至850 m 時，壩體背水坡側的拉伸破壞區域繼續增大（長17.8 m，高4.5 m），此處壩體受到下部表土層的拉伸作用繼續增大所致。迎水坡側壩體變形破壞區域出現急劇增大現象，說明此處壩體受到背水坡側下沉變形，拉伸幅度加大，迎水坡側土體受到的拉應力也在增大，拉伸破壞區域增大。

6）當工作面推進至900 m 時，壩體背水坡側的拉伸破壞區域增大速度減緩，是由于背水坡側壩體受到下部表土層拉應力減弱所致，背水坡側壩體由拉應力轉化為壓應力。迎水坡側壩腳受到背水坡側的拉應力減弱，變形破壞區域增大緩慢。而壩體中部位置開始受到下部表土層的拉應力，出現拉伸破壞，此時壩體中部位置（底部）僅有1 m 左右的土體沒有被拉伸破壞區域影響。

7）當工作面推進至910 m 時，壩體背水坡側壩體繼續保持壓應力狀態，壩體迎水坡側受到的下部表土層拉應力減弱，拉伸破壞區域增長緩慢。壩體中部位置拉伸破壞區域增大，且逐漸向上發展，此時壩體的塑性破壞區域已貫穿壩基，變形破壞區域很有可能成為滲流的主要通道。

8）當工作面推進至920 m 時，壩體迎水坡側壩腳位置開始由拉伸變形破壞轉變為壓縮變形破壞。壩體中部位置的拉伸破壞區域繼續增大，直至工作面推進至930 m 時，整個壩體處于拉伸破壞區域。壩體中部位置也開始由拉伸變形轉向壓縮變形狀態。

隨著工作面開挖結束，壩體一直保持著拉伸破壞狀態。雖然土體發生拉伸破壞后，會出現壓縮、還原的過程，但軟件只能記錄出現過變形破壞，塑性破壞區域無法顯示還原過程，因此，在工作面繼續向前推進的過程中，壩體整體結構均出現了拉伸破壞，僅僅依靠壩體結構的自重使壩體恢復到原始結構狀態，需要經歷很漫長的歷史過程。

青辰想說些什么安慰，卻又不知道說什么才好。他能夠想象到當時的戰況有多慘烈，那日在天葬場，他已經見識了唐飛霄的手段。

3 壩體移動變形觀測

利用topconGPS 水準儀和全站儀對工作面開采過程中壩體的移動變形進行觀測[19-20]，從壩體中心位置向庫區方向在壩頂布設1 條觀測線，設置20 個測點，測點間距為30 m，測線全長600 m，壩下開采壩體移動變形測點布置如圖6。設計工作面推進至600 m 之前，工作面每推進50 m 觀測1 次壩體的移動變形值；工作面推進600 m 至1 800 m 過程中，每推進10 m 觀測1 次；工作面推進1 800 ～2 000 m過程中，每推進50 m 觀測1 次。

圖6 壩下開采壩體移動變形測點布置圖Fig.6 Layout of measuring points of movement and deformation of mining dam body under the dam

選取壩體中部1#、2#和3#測點的觀測結果，分析壩體受采動影響的移動變形規律，工作面推進過程中3 個測點的下沉值見表2。

表2 工作面推進過程中3 個測點的下沉值Table 2 Subsidence values of 3 measuring points during the advancement of working face

依據表2 中工作面不同推進長度下壩體1#、2#和3#測點的下沉值，繪制的壩體上3 個測點受采動影響下的下沉變化曲線如圖7。

由圖7 中測點下沉變化曲線可知：1#、2#和3#測點受采動影響的下沉變化趨勢基本相同，受采動的影響起始點在600 m 左右，在工作面推進至1 500 m 以后各測點處壩體的下沉基本平穩；工作面推至1 800 m 左右時，其中1#測點和2#測點處壩體達到最大下沉，最大下沉值為7.077 m，3#測點處壩體的最大下沉值為7.032 m。

圖7 受采動影響下壩體上3 個測點的下沉變化曲線Fig.7 Subsidence variation curves of three measuring points on the dam body under the influence of mining

實測結果與數值模擬結果對比如圖8。

圖8 實測結果與數值模擬結果對比圖Fig.8 Comparison between measured results and numerical simulation results

從圖8 壩體中部位置的下沉情況對比可知：模擬得到壩體中心位置的下沉過程變化曲線與實測結果中1#測點變化曲線基本吻合，模擬得到壩體最終的下沉值為7.17 m，與壩體真實的最大下沉值7.077 m 相比，絕對誤差為0.093 m，誤差較小。同時，上述壩體的移動變形過程經歷了實測工作面采動后地表下沉曲線的各個位置，符合地表移動變形特征，也證明了局部方法分析方法具有一定科學性和合理性。

4 結 語

1）提出了用淺部開挖替代煤層開采的局部放大分析方法。

2）壩體經歷了從拉伸變形破壞到壓縮-還原的交替變化過程，只是各部位出現拉伸-壓縮轉變過程的時間節點不同，受采動影響開始階段背水坡側壩腳先經歷拉伸變形破壞，其次為迎水坡側壩體，最后為壩體中部位置。經歷由拉伸轉變為壓縮的過程也是由背水坡側壩體開始，其次為迎水坡側壩體，最后為壩體中部位置。

3）數值模擬得到壩體的最大下沉值為7.17 m，與實測結果7.077 m 相比，絕對誤差為0.093 m，誤差較小，驗證了構建的壩下開采壩體變形破壞特征簡化分析方法的科學性和合理性。