基于FLAC3D的貴州馬達嶺煤炭采動型滑坡機理研究

潘網(wǎng)生，傅良同，李 鑫，張曉周，韓國棟

(1.黔南民族師范學院 旅游與資源環(huán)境學院，貴州 都勻 558000；2.長安大學 環(huán)境科學與工程學院，陜西 西安 710054)

貴州省喀斯特出露面積約11×104km2，約占全省國土面積61.9%[1-3]，其喀斯特地貌形態(tài)以錐狀為主，類型多樣、形成發(fā)育過程相對復雜，其形態(tài)結(jié)構(gòu)獨特、錐形明顯，單體高度為60～250m，坡度為45°～ 47°[4-6]。近幾年，貴州喀斯特區(qū)域滑坡地質(zhì)災(zāi)害頻發(fā)，生命與財產(chǎn)損失巨大，如關(guān)嶺縣山體滑坡(2010年6月28日)，死亡99人；福泉市道坪鎮(zhèn)英坪村大型山體滑坡(2014年8月27日)，死亡23人；大方縣理化鄉(xiāng)偏坡村山體滑坡(2016年7月1日)，死亡21人；玉屏縣田坪鎮(zhèn)巖屋口村山體滑坡(2017年7月7日)，死亡2人；納雍縣張家灣鎮(zhèn)大型山體滑坡(2017年8月28日)，死亡15人；三荔高速荔波段滑坡(2017年9月23日)，死亡3人，六盤水六枝特區(qū)關(guān)寨鎮(zhèn)大型山體滑坡(2018年7月9日)，因預警及時以致無人員傷亡。上述滑坡地質(zhì)災(zāi)害已經(jīng)成為制約貴州社會經(jīng)濟可持續(xù)發(fā)展的重要因素之一。由于喀斯特滑坡災(zāi)害機理研究可以揭示滑坡災(zāi)害成因、發(fā)展及演變規(guī)律，也可以為早期預警和地質(zhì)災(zāi)害防治提供科學依據(jù)和理論支撐。因此，研究煤炭采動型滑坡機理具有重要的實踐意義。

國內(nèi)學術(shù)界關(guān)于貴州喀斯特滑坡機理研究，總體概括為巖性和巖土結(jié)構(gòu)等(內(nèi)因)，降雨和不科學的人類工程活動等(外因)。分以下兩種類型：其一是由自然地質(zhì)沉積所形成的軟弱結(jié)構(gòu)層在地表滲流和地下水溶蝕作用下，其粘聚力和內(nèi)摩擦角減小，成為天然滑動面，導致軟弱結(jié)構(gòu)層上方巖土體靜力平衡被打破，進而失穩(wěn)致滑。如貴州關(guān)嶺“6.28” 喀斯特山體滑坡地質(zhì)災(zāi)害等[7-13]；其二是人為斜坡開挖、礦產(chǎn)資源開采等人類工程活動導致邊坡應(yīng)力平衡被打破，在降雨或震動條件下誘發(fā)山體滑坡。如貴州印江巖口大型喀斯特山體滑坡地質(zhì)災(zāi)害等[14-18]。

貴州都勻馬達嶺滑坡屬于第二種類型—人為開挖煤層導致采空區(qū)上方巖體滑坡，也稱煤炭采動型滑坡[19]。趙建軍[20]基于物理模型研究其地質(zhì)力學模式；李禹霏[21]通過自動化監(jiān)測研究其穩(wěn)定性；亓星[22]研究其滑坡碎屑流形成泥石流的特征；王玉川[23]對其開展巖石力學試驗。本文在地質(zhì)野外調(diào)查基礎(chǔ)上，擬采用有限差分方法研究采動滑坡的變形過程及特征，以揭示其滑坡機理和驅(qū)動模式，相關(guān)研究成果可以應(yīng)用于類似喀斯特礦山開采的安全管理和地質(zhì)災(zāi)害防治。

1 數(shù)據(jù)來源及研究方法

研究區(qū)域地質(zhì)資料來源于前人研究成果及貴州省區(qū)域地質(zhì)志(1987版、2017版)[24，25]，地形地貌數(shù)據(jù)來源于國家地理空間數(shù)據(jù)云服務(wù)平臺及野外現(xiàn)場勘察、測繪?；贑AD-ANSYS構(gòu)建原始三維邊坡模型，導入FLAC3D2600軟件后模擬煤層開采對邊坡穩(wěn)定性的影響。

FLAC3D由美國ITASCA公司研究開發(fā)，是一款基于拉格朗日算法有限差分的顯式數(shù)值分析模擬軟件。該軟件包含11種彈塑性材料本構(gòu)模型，非常適合模擬巖土體的屈服、塑性流動、軟化及大變形問題，可以快速求解場的力學方程，準確分析巖土邊坡的穩(wěn)定性[26]。拉格朗日差分計算原理是基于顯式差分求解偏微分方程[27，28]，即通過運動方程計算節(jié)點運動速率，根據(jù)運動時間Δt，計算節(jié)點(ni)時步位移，由節(jié)點位移推算單元(ei)之間的相對位移(單元ei由邊長lj，lj+1，…lj+k構(gòu)成，邊長lj由節(jié)點nk，nk+1和構(gòu)成)，進而求解單元應(yīng)變增量(△εi，j)。由單元材料本構(gòu)方程獲取應(yīng)變條件下的應(yīng)力增量(△σi，j)，在此基礎(chǔ)上求解各節(jié)點的不平衡力(即合力)，根據(jù)中心差分原理[29]，計算節(jié)點的速度與加速度，進而再一次獲得各節(jié)點的時步位移。按照上述步驟往復迭代，擴展到整個區(qū)域，直至各節(jié)點的不平衡力足夠小或者各節(jié)點的位移趨于平衡[30]。

2 研究區(qū)域概況及地質(zhì)模型

馬達嶺滑坡位于貴州省都勻市江洲鎮(zhèn)富溪村干壩馬達嶺溝，溝口地理坐標為26°10′18.9″N，107°17′23.4″E[22](圖1)。馬達嶺溝流域面積0.86km2，流域全長1.62km，最大高程1570m(北部山脊)，最小高程1150m(南部溝口)。溝道走向大致由北向南，中游呈“V”型，上、下游呈“U”型，溝道平均縱向比降約28.2%，上陡下緩。災(zāi)害點地處亞熱帶季風濕潤氣候，雨熱同期，4—9月為雨季，年平均降雨量約1400mm，且主要集中在6—8月[22]。災(zāi)害點地層屬濱海沉積相，地質(zhì)構(gòu)造相對穩(wěn)定，溝口上方有一條東西走向的正斷層。地層地質(zhì)時代劃歸早石炭世大塘早期，為祥擺組(C1x)。地層結(jié)構(gòu)由灰黃色的石英砂巖、黃褐色砂質(zhì)頁巖、灰黑色炭質(zhì)頁巖及黑色煤層組成，其中夾少量泥灰?guī)r。地層主要特征為砂頁巖夾煤層，且被認為是貴州石炭系出露最廣、唯一含煤的巖組。由Ⅰ-Ⅰ′地質(zhì)剖面(圖2)可知，巖層產(chǎn)狀280°∠18°，緩傾坡內(nèi)。煤層主要有3層，自下而上分別為M1、M2和M3，其中M1和M2煤層厚度1.9～2.4m，M3煤層厚度0.2～0.9m不等(圖3)。因煤炭開采，坡頂產(chǎn)生縱橫交錯、寬度不等的裂隙，發(fā)育形成兩組優(yōu)勢結(jié)構(gòu)面[20]，大致產(chǎn)狀分別為34°∠87°、118°∠86°，受其影響滑坡體與山體拉裂，形成“圈椅狀”滑坡后壁。滑坡位于馬達嶺坡頂，滑坡體后緣肩部高程約1573.6m，滑坡前緣坡腳高程約1481.5m，相對高差約92.1m(圖3)?；禄瑒又鞣较蚣s204.6 °，滑坡體沿該主方向長度約88m，平均寬度102m，平均厚度92m，滑坡體體積達8.26×105m3。此外，滑坡物堆積于馬達嶺溝中部，堆積體平均厚度約6m[23]。

圖1 滑坡區(qū)域平面圖(滑坡前)

圖2 Ⅰ-Ⅰ′剖面地層結(jié)構(gòu)

圖3 滑坡全貌圖

由于FLAC3D構(gòu)建不規(guī)則三維模型具有一定的難度，因此本文以Ⅰ-Ⅰ′剖面為基礎(chǔ)，在ANSYS軟件平臺上恢復原始三維邊坡概化模型，之后再導入FLAC3D(圖4)。該邊坡模型x軸方向長0～426.49m，y軸方向?qū)?～100m，z軸方向高1400.00～1577.09m，模型共產(chǎn)生381786個節(jié)點和365950個單元，分19個地層，分別為砂巖、炭質(zhì)頁巖和煤巖互層。巖性參考文獻[23][31]，本文數(shù)值模擬計算時假設(shè)各層砂巖的巖性一致，各層炭質(zhì)頁巖的巖性一致，各層煤巖的巖性也一致(表1)。由于FLAC3D關(guān)于巖體變形模量僅接受體積模量K和剪切模量G，因此本文設(shè)定砂巖、炭質(zhì)頁巖及煤巖均彈塑性材料，服從莫爾-庫倫屈服準則，并依據(jù)文獻[32]確定體積模量K、剪切模量G與彈性模量E及泊松系數(shù)μ的關(guān)系。

圖4 初始概化三維邊坡模型

3 煤層開挖過程的邊坡變形分析

為深入分析煤層開挖對邊坡穩(wěn)定性的影響，本文設(shè)計5種工況，分別探討開挖進深、開挖作業(yè)面寬度及開挖煤層的先后次序等對馬達嶺邊坡變形的影響。

表1 巖石物理力學參數(shù)

3.1 工況設(shè)計

1)工況一：開挖M1、M2煤層至x=300m，作業(yè)面寬度8m。通過FLAC3D模擬M1、M2煤層開挖過程，并監(jiān)測邊坡應(yīng)力應(yīng)變狀態(tài)。M1、M2煤層開挖最終各形成5個工作面，每個工作面寬8m，工作面之間以煤柱作為支撐，其中，M1煤層開挖至x=300m，開挖進深為245m，M2煤層也開挖至x=300m，開挖進深為187m，工作面y軸方向的坐標分別為：4～12m；24～32m；44～52m；64～72m；84～92m，開挖位置及尺寸如圖5(a)所示。模擬過程中，在邊坡體中選取有代表性的點位A(256，0，95)、B(280，0，130)和C(360，0，140)(見圖2)作為監(jiān)測點，以記錄節(jié)點豎向位移、節(jié)點合速率等參數(shù)隨運行步的變化，其在格網(wǎng)中的節(jié)點編號分別為A(GP4104)、B(GP1815)、C(GP1414)。

2)工況二：開挖M1、M2煤層至x=400m，作業(yè)面寬度8m。M1煤層開挖至x=400m，開挖進深345m，M2煤層也開挖至x=400m，開挖進深287m，其他條件同開挖進深至x=300m一致，開挖位置及尺寸如圖5(a)所示。

3)工況三：開挖M1、M2煤層至x=400m，作業(yè)面寬度增加至12～16m。M1煤層開挖至x=400m，開挖進深為345m，M2煤層也開挖至x=400m，開挖進深為287m，工作面12～16m不等，y軸方向坐標分別為：4～20，24～40，44～60，64～80，84～92，開挖位置及尺寸見圖5(b)各工作面之間以煤柱作為支撐。其他條件同工況二。

4)工況四：在工況二基礎(chǔ)上開挖M3煤層，開挖至x=300m。在工況一基礎(chǔ)上開挖M3煤層，開挖至x=300m，M3煤層工作面在y軸方向分別為：4～12m；24～32m；44～52m；64～72m；84～92m，開挖位置及尺寸如圖5(c)所示。

5)工況五：同時開挖M1、M2、M3煤層，其他開挖條件同工況四，開挖位置及尺寸如圖5(c)所示。

圖5 各工況煤層開挖位置及工作面寬度(m)

3.2 邊坡變形結(jié)果分析

3.2.1 開挖進深對邊坡變形的影響

根據(jù)工況一和工況二的邊坡變形計算結(jié)果(圖6)，當M1、M2煤層開挖進深由300m增加至400m，邊坡整體變形增加了約2倍。開挖之初，坡腳前緣部分最先產(chǎn)生塑性變形，但由于坡腳前緣煤層上方的巖體厚度較薄，較小的勢能和載荷使得坡腳變形相對較小。隨著開挖進深增加，邊坡變形越來越大，且變形最大區(qū)域逐步由坡腳前緣向坡頂和邊坡內(nèi)部轉(zhuǎn)移。

圖6 不同開挖進深下的邊坡合位移

根據(jù)邊坡內(nèi)部A、B、C節(jié)點監(jiān)測數(shù)據(jù)，節(jié)點產(chǎn)生位移的先后次序分別為A、B、C，且由于節(jié)點在y軸方向的水平位移最小，在x軸方向的水平位移次之，而在z軸豎向位移最大，因此本文選取豎向位移作為節(jié)點位移監(jiān)測指標，另選取合速率作為節(jié)點速率監(jiān)測指標。由工況一的邊坡變形和節(jié)點豎向位移、合速率計算結(jié)果(圖6、圖7)可知，當M1、M2煤層開挖進深由300m增加至400m，節(jié)點合速率增加了約2倍。節(jié)點豎向位移與節(jié)點合速率與開挖進深均呈正比例關(guān)系(圖7)，即隨開挖進深的增加而增大。坡頂部位合位移將沿Z軸負向推進，其位移大小在任一Y剖面上沿X軸正向呈遞減分布，即形成不同的合位移區(qū)間。因此，邊坡巖體存在合位移差異，導致邊坡垂向節(jié)理和裂隙形成并發(fā)育，且隨著開挖進深的增加而增大，逐步向邊坡內(nèi)部推進，直至下切到M1、M2和M3煤層，形成潛在的優(yōu)勢滑動結(jié)構(gòu)面。

圖7 不同開挖進深下的節(jié)點豎向位移與合速率監(jiān)測

3.2.2 開挖作業(yè)面寬度對邊坡變形的影響

根據(jù)工況三的邊坡合位移計算可知，當M1、M2煤層的作業(yè)面開挖寬度由8m增加至12～16m時(2倍)，邊坡合位移增加了約10倍(圖8)。

圖8 增加作業(yè)面開挖寬度下的邊坡合位移

根據(jù)邊坡內(nèi)部A、B、C節(jié)點豎向位移監(jiān)測數(shù)據(jù)記錄，增加作業(yè)面寬度后，A、B節(jié)點的豎向位移增大了約11倍，C節(jié)點的豎向位移增加了約12.3倍。此外，增加作業(yè)面寬度后(2倍)，峰值合速率增大約4倍(圖9)。上述數(shù)據(jù)說明，邊坡內(nèi)部節(jié)點豎向位移與作業(yè)面寬度呈正比關(guān)系，即作業(yè)面寬度越大，節(jié)點豎向位移就越大，即作業(yè)面寬度對邊坡變形的影響要遠大于開挖進深。

圖9 增加開挖寬度下的節(jié)點豎向位移與合速率監(jiān)測

由此表明，隨著M1、M2煤層開挖作業(yè)面寬度的增加，邊坡合位移、節(jié)點豎向位移和節(jié)點合速率均呈不同程度增加，且節(jié)點豎向位移增幅遠大于節(jié)點合速率增幅。

3.2.3 開挖煤層先后次序?qū)吰伦冃蔚挠绊?/p>

由工況四的邊坡變形計算(圖10a，b)可知，邊坡最終合位移由兩次開挖疊加導致，且開挖M3煤層所產(chǎn)生的合位移主要集中于M3煤層上方巖體，合位移約為單獨開挖M1、M2煤層合位移的10 %～20%，說明深部煤層開挖對邊坡變形的影響要遠大于淺層煤層開采。此外，根據(jù)工況五的邊坡合位移計算結(jié)果(圖10c)，同時開挖煤層M1、M2和M3(工況五)產(chǎn)生的最大變形區(qū)域主要分布于邊坡后緣中下部，其整體變形較分步驟開挖煤層產(chǎn)生的變形(工況四)稍大，節(jié)點豎向位移和節(jié)點合速率記錄也驗證此現(xiàn)象(圖11)?？傮w而言，開挖煤層先后次序?qū)吰伦冃未嬖谝欢ǖ挠绊?，且同時開挖各煤層對邊坡變形的影響程度較分步驟開挖各煤層要大，但相對于開挖進深和開挖作業(yè)面寬度而言，其對邊坡變形的影響程度又相對較小。

圖10 不同開挖次序情況下的邊坡合位移

圖11 不同開挖次序下的節(jié)點豎向位移與合速率監(jiān)測

4 基于數(shù)值分析的滑坡機理推定

本文在現(xiàn)場調(diào)查和定期的邊坡變形監(jiān)測(GPS)基礎(chǔ)上，開展基于FLAC3D的數(shù)值模擬分析，得到與現(xiàn)場調(diào)查、監(jiān)測相一致的結(jié)論。說明本文選取邊坡合位移、節(jié)點豎向位移和節(jié)點合速率等指標監(jiān)測邊坡變形過程具有可操作性和實際意義，為揭示邊坡滑動機理和本質(zhì)特征提供了重要的觀測和量化工具。由此，借助該工具推定馬達嶺滑坡機理和過程如下：

單獨開挖M1、M2煤層時，隨煤層開挖進深的推進，邊坡內(nèi)部因變形位移差異產(chǎn)生沿x軸方向的若干條垂向優(yōu)勢結(jié)構(gòu)面(e、f、g、h)(圖12)。該類優(yōu)勢結(jié)構(gòu)面受煤層開挖作業(yè)面寬度控制，隨作業(yè)面寬度的增加而拉裂擴張。其向上延伸至坡頂，形成若干條地裂縫，成為坡面雨水入滲邊坡的優(yōu)勢通道，向下則切割邊坡，直至M1、M2煤層。

圖12 滑坡機理示意圖

在此基礎(chǔ)上若進一步開挖M3煤層，邊坡上方前緣再次遭受采掘擾動。此時位于M3開挖進深范圍內(nèi)、前期所形成的垂向優(yōu)勢結(jié)構(gòu)面，將隨M3煤層作業(yè)面寬度的增加而進一步拉裂擴張。受M3煤層采空影響，M3煤層上方巖體沿z軸負向產(chǎn)生位移，同時M3下方巖體因采空卸荷而反彈，產(chǎn)生若干卸荷拉張裂隙面，即橫向優(yōu)勢結(jié)構(gòu)面(d的下部)。在滿足一定的降雨條件時，M3煤層采空區(qū)匯集大量裂隙水，一部分裂隙水沿M3煤層開挖工作面入口溢出，另一部分則沿橫向優(yōu)勢結(jié)構(gòu)面繼續(xù)下滲。當M3煤層上方巖體發(fā)生失穩(wěn)垮塌，在重力作用下沖擊M3煤層下方巖體時，橫向優(yōu)勢結(jié)構(gòu)面鎖固段被打通，優(yōu)勢結(jié)構(gòu)面(d)在水、巖等綜合作用下實現(xiàn)連續(xù)和貫通。此時，上方巖體沿優(yōu)勢結(jié)構(gòu)面在剖面D處臨空剪出，剪出巖體受重力作用拋向馬達嶺溝谷，遭受撞擊后形成大小不等塊體散落于溝谷中段，滑坡地質(zhì)災(zāi)害由此發(fā)生。

馬達嶺滑坡地質(zhì)災(zāi)害發(fā)生后，該邊坡災(zāi)害隱患依然存在，甚至更為嚴峻。一方面，拋落在溝谷中的大小巖塊堆積物在極端天氣條件下可能會再次發(fā)生滑動，威脅溝谷下游安全。據(jù)2018年6月22日強降雨后的現(xiàn)場觀測記錄數(shù)據(jù)，該原滑坡堆積區(qū)(D下方)邊緣產(chǎn)生5～10cm寬度不等的蠕滑裂縫。另一方面，開挖M1、M2過程中所產(chǎn)生的一系列垂向優(yōu)勢結(jié)構(gòu)面e、f、g、h等，已經(jīng)將邊坡切割為若干部分，這些被切割巖體之所以能處于相對穩(wěn)定狀態(tài)，是因為邊坡巖層反傾向于坡內(nèi)，并缺少足夠水文地質(zhì)條件。然而，該邊坡及后緣山體巖層產(chǎn)狀為NW282°∠18°，在水、巖長期相互作用下，山體垂向優(yōu)勢結(jié)構(gòu)面和巖層水平節(jié)理、裂隙會持續(xù)擴張。滿足一定降雨條件后，地表水將沿后緣山體垂向優(yōu)勢結(jié)構(gòu)面滲入巖體層狀節(jié)理、裂縫，引起巖層順層錯動，直接擠壓和推動前述邊坡的被切割巖體沿優(yōu)勢結(jié)構(gòu)面(如：g-M1)剪出，到時可能會形成規(guī)模更大的滑坡地質(zhì)災(zāi)害。

5 結(jié)論與討論

5.1 結(jié)論

1)將基于快速拉格朗日差分計算的FLAC3D軟件應(yīng)用于貴州省都勻市馬達嶺滑坡機理分析，通過記錄邊坡變形最大不平衡力、單元和節(jié)點的豎向位移、合速率等參數(shù)變化來揭示邊坡滑動機理和過程，具有明顯的技術(shù)比較優(yōu)勢，可以在西南喀斯特崩塌和滑坡等地質(zhì)災(zāi)害研究中進一步推廣和應(yīng)用。

2)煤層開挖進深控制邊坡垂向優(yōu)勢結(jié)構(gòu)面的空間分布，煤層開挖作業(yè)面寬度控制垂向優(yōu)勢結(jié)構(gòu)面擴張程度，后者對邊坡穩(wěn)定性的影響遠大于前者。

3)深層煤層開挖對邊坡變形的影響遠大于淺層煤層開挖。此外，煤層開挖的先后次序?qū)吰伦冃未嬖谝欢ǖ挠绊?，但相對于開挖進深和開挖作業(yè)面寬度而言，其影響較小。

4)馬達嶺滑坡的形成與發(fā)育可以分為兩個階段，第一階段培育了垂向優(yōu)勢結(jié)構(gòu)面(開挖M1、M2煤層)和橫向優(yōu)勢結(jié)構(gòu)面(開挖M3煤層)；第二個階段(M3上層巖體重力失穩(wěn)垮塌，水、巖相互作用)實現(xiàn)橫向、垂向優(yōu)勢結(jié)構(gòu)面連續(xù)和貫通，在水、巖等綜合作用下，巖體沿優(yōu)勢結(jié)構(gòu)面臨空剪出。

5.2 討論

1)貴州省煤炭資源廣泛分布，有些礦區(qū)就在旅游景區(qū)周邊，甚至位于景區(qū)必經(jīng)之路旁，對旅游景區(qū)和游客安全構(gòu)成極大威脅。正如馬達嶺坡頂西側(cè)山體正在開發(fā)農(nóng)業(yè)觀光項目，建議相關(guān)部門應(yīng)該針對此項風險進行專門評估。

2)開挖M1、M2煤層所培育的大量垂向優(yōu)勢結(jié)構(gòu)面與M1、M2煤層采空區(qū)交會，實現(xiàn)了對后緣邊坡的完全切割。不利的巖層產(chǎn)狀(NW282°∠18°)及大量節(jié)理裂隙水的存在，可能正在孕育一個規(guī)模更大的滑坡地質(zhì)災(zāi)害。為此，將在后期通過持續(xù)的調(diào)查和監(jiān)測來修正和完善上述結(jié)論。

3)基于FLAC3D的邊坡變形過程分析是否適用于含溶洞的喀斯特邊坡穩(wěn)定性分析，以及是否適用于喀斯特溶洞結(jié)構(gòu)安全評價，尚無相關(guān)研究報道，值得進一步深入研究。此外，F(xiàn)LAC3D在這一領(lǐng)域的應(yīng)用，其可靠性依賴于現(xiàn)場試驗技術(shù)和手段的先進性，特別是巖體微細觀節(jié)理、裂隙等空間探測技術(shù)的重大突破。