撫順西露天礦開挖作用下邊坡變形破壞滑動規律研究

郭子鈺 劉向峰 王來貴 王 勇2

（1.遼寧工程技術大學力學與工程學院，遼寧 阜新 123000；2.中煤科工集團沈陽設計研究院有限公司，遼寧 沈陽 110015）

近年來，露天采礦規模迅速擴張，開采技術越來越先進，露天采場越來越深，邊坡的高度和擾動面積逐漸增大，引發的滑坡事故增多[1-2]，許多露天煤礦邊坡[3-4]都經歷過變形，當變形量累積到一定程度時，萌生裂紋，裂紋之間相互擴展連接，再當坡體被裂紋貫穿后，就在重力或其他外界因素擾動下發生失穩滑動。由于滑坡成因和誘因十分復雜，邊坡的變形、破壞、滑動過程也具有階段性，處于不同變形階段的邊坡，距離邊坡整體破壞的時間間隔也不同，所以滑坡時間預報和空間預報至今還是一個世界性的難題[5]。

目前，曹蘭柱等[6]應用數值模擬的方法，模擬了寶日希勒露天煤礦的應力應變和位移變化規律，結果表明南幫邊坡失穩時，邊坡失穩以剪切作用為主，拉伸破壞相伴發生。王玉凱等[7]采用模型試驗對比天然和飽和工況下的排土場變形破壞過程，應用數字散斑和測點追蹤等軟件，監測位移分布特征及演化規律，結果表明排土場變形破壞形態呈現滑塌—牽引—推移態勢。張志飛等[8]采用顆粒流程序（PFC）對反傾層狀巖質邊坡變形破壞過程進行模擬，得到了邊坡巖層的主要變形破壞方式為彎曲變形、折斷破壞，變形首先在坡頂發生，破壞從坡腳開始，邊坡在變形破壞過程中表現為上部受拉下部受壓的特征。張碩等[9]對坡體裂縫發育與演化及變形破壞模式進行模型試驗，結果表明坡體裂縫體系演化機制為臨空面附近產生張拉裂縫且向后擴展—坡體側翼產生剪張裂縫—后緣產生貫通性的張拉裂縫。楊天鴻等[10]通過位移監測手段在弱層流變的基礎之上研究阜新海州露天礦邊坡變形破壞機理，劃分出了滑移變形、潰屈變形、隆起變形和剪切變形。

上述學者通過數值模擬、物理模型、位移監測等手段對邊坡變形、破壞、滑動進行研究并取得了豐碩的研究成果，但是在實際工程當中，邊坡變形、破壞、滑動過程是一個連續不間斷的一個災害演化，因此在前人研究的基礎上本項目應用連續—非連續方法（CDEM），以撫順西露天礦為研究對象，研究露天礦南幫邊坡變形、破壞、滑動全過程，從南幫邊坡的變形孕育到裂紋破壞，再到最后失穩滑動進行系統的分析和討論，從而對南幫邊坡以及相似工程背景的露天礦邊坡提供數據參考。

1 CDEM方法

中科院力學研究所提出了基于連續介質力學的離散元方法，簡稱CDEM[11-13]。CDEM的數值模型可以分成塊體部分、界面部分，如圖1所示[14]，真實界面是塊體真實存在的不連續面，如斷層、節理等，其物理力學參數與塊體的參數一致。虛擬界面在塊體中并不存在，主要有兩個作用：①把相鄰塊體連接，傳遞力學信息，②為裂紋的發育擴展提供潛在的備選通道。本研究應用虛擬界面和塊體單元來描述邊坡裂紋的擴展以及滑動，真實界面轉換為塊體單元進行計算，在坡面區域將單元細化，最大限度提供裂紋的擴展路徑。

CDEM方法的核心計算公式為

式中，u?(t)、u?(t)和u(t)分別是單元內所有節點的加速度列陣、速度列陣和位移列陣；M、C、K和F（t）分別為單元質量矩陣、阻尼矩陣、剛度矩陣和節點外部荷載列陣。

CDEM采用顯式歐拉前向差分法進行問題的求解，整個計算過程中通過不平衡率表征系統受力的平衡程度，計算流程如圖2所示。

2 變形破壞滑動機理

根據文獻[15]可知，邊坡失穩滑動的機理分為變形、破壞以及滑動，在變形階段滿足彈性準則，變形分為法向變形和切向變形；在破壞階段滿足摩爾庫倫準則，分為拉張破壞和剪切破壞；在滑動階段滿足動力學方程，分為穩定滑動和非穩定滑動，即負反饋和正反饋系統。

變形階段：滿足τ＜σ·tan?+c，其中τ為弱面剪應力，σ為弱面法向應力，c為弱面的黏聚力，?為內摩擦角。

破壞階段：滿足σ3＞σt，則塊體發生拉伸破壞，裂紋的方向與最小主應力σ3垂直，其中σ3為最小主應力，σt為塊體的抗拉強度；滿足σ1≥σ3tan2(π /4+?/2)+2ctan(π /4+?/2)，則考慮剪切破壞，裂紋方向與σ3夾角為π/4+?/2，其中σ1為最大主應力。

2.1 單元變形演化規律

在變形階段時，邊坡單元為線彈性模型，根據變形機理發生法向和切向變形，對應到數值計算中就是節點位置坐標發生更改，本節討論單元變形的演化過程。結合CDEM方法，將邊坡劃分單元網格進行迭代計算，單元發生法向和切向變形，如圖3所示，從圖中可以明顯地看出塊體單元圖3（a）到塊體單元圖3（b）再到塊體單元圖3（c）的變化。塊體單元圖3（a）中共有5個小三角形單元，節點數為6個，塊體單元圖3（c）三角形單元也有5個，節點數也為6個，但是單元圖3（a）的三角形單元開始階段都是正三角形單元，受到重力作用，導致各個節點受力不同，最終使得單元形狀發生改變，變得扁長，最后形成單元圖3（c）的形態，塊體單元圖3（a）的重心向右下方向移動一小段距離。在變形階段時，邊坡內各個單元的平均應力均滿足τ＜σ·tan?+c，經過法向和切向變形，逐漸累積變形能量。總之，在變形階段地質體受到重力作用影響，會發生切向和法向變形，邊坡地質體會產生微小變形，映射到數值計算中塊體單元在變形階段，單元間是連續的，節點與節點之間不會發生分離，在單元內部都是共享節點，單元的節點由于受到外力作用，導致位置坐標發生更改，單元的形狀就會發生微小的改變。

2.2 單元破壞演化規律

如圖4所示，在變形階段計算平衡后即進入到破壞階段，在單元界面設置接觸彈簧（接觸面），對應的接觸彈簧滿足拉破壞準則或剪破壞準則。

滿足破壞準則時界面發生斷裂，也就是當變形量逐漸增加，達到一定閾值時，單元間界面發生破壞，邊界面上的接觸彈簧斷裂，形成裂紋，由變形到破壞，即單元由連續狀態變為非連續狀態。由于單元之間所受到的應力不同，導致接觸面破壞分成拉張破壞和剪切破壞，對應的裂紋也分成拉張裂紋和剪切裂紋。即單元間發生破壞，裂紋就會沿著單元的邊界面進行擴展，即虛擬界面擴展。

在單元變形之后繼續發展，當滿足拉張破壞條件時，單元的節點會發生分離，單元與單元之間出現縫隙如圖4（b），即單元破壞。單元破壞是一個時間點，在這個特定時間點內單元破壞，即接觸面滿足拉張破壞準則，單元由連續狀態變成非連續狀態。隨著外力的擾動，連續區域內的單元滿足變形準則，繼續發生法向和切向變形，導致非連續區域單元間的縫隙越來越大，隨著單元應力變化，節點的應力也在改變，當相鄰單元間再次滿足破壞準則時，接觸界面又再次發生斷裂，就再次萌生裂紋，當裂紋逐漸連接成一條大裂紋時，就標著破壞區域完全貫穿破壞（圖4（c））。如果繼續發展就有可能發生滑動，就進入到滑動階段。

圖5為沿著單元界面的剪切破壞演化，和拉張破壞一樣，當滿足剪切破壞條件時，單元發生剪切破壞如圖5（b）所示。在圖5（b）的左側開始剪開一個小裂紋，在小裂紋的周圍單元是連續狀態，滿足變形準則，發生法向和切向變形，變形量繼續增加，當下一個界面滿足破壞條件時，則界面彈簧斷裂，形成裂紋，以此類推實現剪切裂紋的擴展，形成圖5（c）中的形態。

2.3 單元滑動演化規律

對于滑動階段則是在破壞的基礎上延續，邊坡破壞產生的大裂紋上下貫穿坡體，形成滑體，受到重力或其他擾動作用，滑體沿著弱層面向下發生滑動。滑動階段的單元演化主要是以滑體在破壞路徑上整體滑動（即裂紋擴展方向），局部單元產生微小變形，如果再次達到界面破壞準則，則繼續開裂形成新的裂紋，新的裂紋繼續擴展連接，形成新的離散區域，繼續滑動。由于技術原因，在單元中無法監測阻尼力變化，故監測速度變化以等效阻尼力變化。

南幫滑移方向為x正方向，以x正方向為例，式可以看出，下滑力P(t)基本不再改變，如果a增大，阻尼力減小，就可以等效阻尼力隨著速度增大而減小，系統為正反饋系統。

圖6（a）為局部剪切破壞最終狀態，圖6（b）和圖6（c）為滑動演化的單元形態。在剪切破壞結束時，裂紋完全貫穿破壞區域，由圖6（a）到圖6（b）在主裂紋區域相鄰的界面都已經完全分離，界面上的節點都為非共享節點，受到重力作用，再經過迭代計算，連續區域發生法向和切向變形，非連續區域每個單元的下滑力大于阻尼力，導致沿著裂紋面方向滑動，如圖6（c）所示，在滑移的過程中伴隨著變形和進一步的破壞。單元滑動演化過程又細分為穩定滑動和非穩定滑動，圖6只是滑動演化，沒有體現出穩定或非穩定滑動，在實際計算當中，設置監測點監測速度變化來描述穩定或非穩定滑動，判斷依據根據速度曲線斜率來判斷穩定或非穩定滑動。

3 開挖作用下南幫邊坡數值模擬

3.1 南幫滑坡模式分析

國內許多學者對撫順西露天礦南幫進行滑坡分析[16-18]，結合撫順西露天礦的工程地質因素等，本研究提出南幫的滑坡模式主要是滑移—剪斷模式，如圖7所示。南幫現在形成一個大滑體，其潛在滑動面埋深較大且其傾角（30°左右）大于坡角（19°～27°），構成“隱伏型”順傾坡[19]。

從圖7中可知，巖層自上而下為凝灰巖區域、弱層區域、玄武巖區域以及花崗片麻巖區域，南幫滑坡大概分成五部分，分別為拉張變形區域、沉陷變形區域、滑移變形區域、隆起變形區域以及剪切變形區域。邊坡頂部區域受到重力作用，使得頂部發生沉降變形差，形成地裂縫，在拉張變形區域，凝灰巖上部區域受到邊坡下部區域的拉力，產生拉張變形，使其發生沉陷變形，凝灰巖中部區域沿著弱層向下滑移，在坑底處，受到北幫的影響，滑動產生了抑制作用故產生隆起變形，最后在坑底區域剪出。

3.2 計算工況及模型參數

結合撫順西露天礦現狀，建立E1600南幫計算模型，以左下角為原點，底邊長為2 050 m，高為990 m（頂部標高為+126 m），單元為2 827個，節點為1 488個，如圖8所示。

南幫地層自上而下分布依次為凝灰巖（分3次開挖）、弱層、玄武巖、弱層和花崗片麻巖。

開挖作用下邊界條件：考慮重力場，左右兩側邊界x方向固定約束，底部邊界x、y方向固定約束，上部為自由邊界。南幫從頂部向下開挖分3次，第一次開挖130 m，第二次開挖160 m，第三次開挖140 m。

南幫監測點設置如圖8中所示，共計6個。

南幫巖土的物理力學參數根據礦區試驗資料確定算例的計算參數，如表1所示。

?

3.3 南幫變形、破壞、滑動規律演化

圖9為開挖作用下，隨著時步增加，南幫邊坡的變形破壞滑動過程，細分為變形壓密、變形穩定、破壞、第二次破壞滑移、第三次破壞滑移5個階段，當滿足收斂平衡條件后，就進入到下一個階段，后續以這5個階段來描述南幫邊坡的變形破壞滑動過程。

3.3.1 位移變化曲線

（1）坡頂區域。由圖10可知，坡頂區域x方向和y方向位移的總體變化趨勢，靠近坡面位置隨著開挖深度增加而逐級增大。在變形壓密階段，由于N1的片麻巖強度高于N2處的弱層，故N2的x方向變形量大，而N1處基本不變；y方向N2處的沉降量也高于N1處的。在第一次開挖變形穩定階段，x和y方向的數值基本不變，這時主要是積蓄能量，在破壞時釋放。進入到破壞階段，由于產生裂紋，釋放能量，導致巖體會產生回彈效應，使得位移曲線產生波動，x方向和y方向的變形量回彈，后來在重力的作用下，繼續壓實。在第二次開挖時，由于上部巖體被挖空，下部會釋放壓力，又使得巖體回彈，當第三次開挖時，又反復一次。經過3次的開挖，導致坡頂區域N2處產生x方向2 m，y方向-2.6 m的位移值；N1處x方向基本不改變，y方向-1.6 m。

（2）坡中區域。由圖11可知，坡中區域經過3次的開挖發生大規模滑動，坡中區域的N3和N4位移曲線差別很大，曲線間隔為滑動距離。在變形壓密階段，x方向基本不動，主要是y方向發生沉降變形，此時積蓄能量；進入第一次開挖穩定變形階段，x和y方向在原來的基礎上穩定不變，當到達破壞時釋放能量；在破壞階段，有回彈現象，由于在坡中位置，故第一次開挖，對N3和N4的影響不大，基本和坡頂處的位移一致；在第二次開挖時，坡中位置完全臨空，N4更靠近坡面位置，這時N4的x和y方向位移迅速增加；到第三次開挖時，滑動得就更加明顯。

（3）坡底區域。由圖12可知，在坡底區域經過3次開挖，坡底區域N5和N6的x方向位移曲線差別很大，曲線間隔為滑動距離。其中N6為坡底，在整個過程x方向基本不動，N5為坡趾處，到第三次開挖才對其產生影響，x方向位移迅速增長；y方向的變化情況：在變形壓密階段，發生沉降變形，經過開挖，底部產生回彈效應，由于上部和中部的滑體經過3次的開挖已經產生滑動，使坡底處有大量的剪切裂紋，坡底也隆起形成底鼓。

3.3.2 速度變化曲線

坡中區域位移量最大，故監測N4點的x方向和y方向速度變化，來描述整體邊坡在開挖情況下的滑動情況。本研究通過速度曲線斜率來描述邊坡的滑動特征。算例的一時步需要的時間為1×10-4s，共計5.6 s。

由圖13（a）可知N4點x方向的速度變化情況。在壓密階段，尋求計算平衡狀態，故速度波動較大。第一次開挖穩定階段，速度的變化基本不大，這時在積蓄變形能量，進入到破壞階段，由于裂紋萌生，釋放能量，速度增加，阻尼力也增加，增加到一定值時曲線斜率＜0，速度又逐漸降低，屬于負反饋系統。到第二次開挖時，巖體挖空，導致N4區域完全臨空，速度又迅速增加，由于阻尼力的影響，最后速度不變，這時x方向處于臨界狀態。在進入到第三次開挖階段，滑體開始并未整體滑移，N4的x方向速度逐漸降低，隨著時間推移，裂紋越來越多，上下部完全貫穿，滑面形成，曲線斜率＞0，最后速度迅速增大，演化成滑坡災害，這時滑坡屬于正反饋系統，邊坡失穩，最后x方向位移12.5 m。

由圖13（b）可知N4點y方向的速度變化情況。在變形壓密階段，主要是發生沉陷變形，y方向的速度變化較大；壓密階段結束后，速度基本不再變化，一直到第一次開挖穩定，隨后進入到破壞階段。由于裂紋的影響，裂紋萌生釋放變形能，巖體向上回彈，y方向的速度向上也迅速增加，而受到阻尼力的影響速度又降低，這時在y方向也是屬于負反饋系統。第二次開挖時，由于裂紋繼續延伸擴展，繼續釋放能量，巖體又一次回彈一部分，速度也跟著增加，隨后在重力和阻尼力的作用下，速度又向下增加，此時還是屬于負反饋系統，阻尼力一直阻礙速度變化。進入到第三次開挖時，由于開始產生巖體裂紋并完全貫穿，故速度逐漸減低，但是隨著裂紋增多，破裂區域逐漸連接，滑面形成，出現滑坡災害，速度就一直增大，邊坡失穩，最后y方向的位移為-11.5 m。

3.3.3 位移變化云圖

分析不同時段的云圖變化來描述災變過程，云圖選取12 000時步、16 000時步、26 000時步、36 000時步和56 000時步分別對應變形壓密階段最后、第一次開挖穩定階段最后、破壞階段最后、第二次開挖破壞滑移階段最后和第三次開挖破壞滑移階段最后，如圖14所示。

圖14（a）為變形壓密最后階段，整體來看，整個地層處于壓實狀態，由于地層的巖層不同，故在凝灰巖區域x方向變形較大，變形值為0.971 m。圖14（b）為第一次開挖穩定狀態最后，此時經過開挖作用，變形已經基本穩定，把地層的凝灰巖區域挖空，使得玄武巖臨空，故在坡面位置x方向變形較大，變形值為0.677 m。圖14（c）為破壞階段最后，可以明顯地看出，裂紋布滿坡面位置，尤其邊坡頂部處有大量的拉張裂紋，并向下延伸，但并沒貫穿邊坡，x方向位移逐漸向坡中區域靠攏，這是坡中由凝灰巖層和弱層構成，力學參數較弱，故變形值較大，變形值為2.147 m。圖14（d）為第二次開挖滑移結束階段，明顯地看出，裂紋在破壞階段的基礎上繼續擴展連接，形成新的大裂紋，裂紋也逐漸的向坡底處延伸，從圖中也可以看出坡底形成大量的剪切裂紋，整體的位移情況也在延續破壞階段的位移，并逐漸形成滑體輪廓，在此階段坡中的位移值仍是最大，位移值為5.257 m。圖14（e）為第三次開挖滑移最后，可以看出整個滑體沿著弱層面向下滑移，坡頂、坡中和坡底已經布滿裂紋，最后邊坡失穩滑動，最大位移值為16.04 m還是在坡中，這是由于在第二次開挖結束時裂紋已經基本貫穿邊坡，經第三次開挖，使得邊坡中下部臨空，沒有巖體阻礙滑動，形成滑體后，迅速向中下滑動，但是由于北幫和底部巖層的影響，使其坡中區域向下滑移在底部推擠，底部也形成底鼓。

3.3.4 塑性區變化云圖

分析不同時段塑性區云圖來描述災變過程。故云圖選取16 000時步、26 000時步、36 000時步和56 000時步分別對應第一次開挖穩定階段最后、破壞階段最后、第二次開挖破壞滑移階段最后和第三次開挖破壞滑移階段最后。如圖15邊坡整體塑性區云圖所示。圖中色帶的不同值分別表示邊坡的破壞狀態：0表示坡體未破壞；1表示當前拉壞；2表示當前剪壞；4表示過去拉壞；8表示過去剪壞。

圖15（a）為第一次開挖變形穩定最后，塑性區已經出現，并沿著弱層區域向下擴展，在坡頂區域出現拉張破壞塑性區，但無裂紋產生，當前剪切破壞塑性區沿著弱層向下。圖15（b）為破壞階段最后，可以看出裂紋已經在剪切和拉張塑性區處產生，裂紋也隨著塑性區的擴展而向下擴展。圖15（c）為第二次開挖破壞滑移最后，可以看出當前剪切破壞基本貫穿整個邊坡，塑性區面積也增加明顯，裂紋也布滿坡面位置，坡底有大量的剪切裂紋，坡頂有大量的拉張裂紋，在第二次開挖時邊坡就已經基本破壞，但是由于下部還有巖層的存在，故沒有發生大規模滑移，如果再次開挖，裂紋完全貫穿邊坡就會發生滑坡災害。圖15（d）為第三次開挖破壞滑移最后，看出邊坡沿著弱層面向下滑移，滑體中部的滑移量最大，整個滑體最后推擠在坡底處，而坡底由于受到剪切裂紋和上部巖體的擠壓也形成底鼓。為了更直觀地分析開挖過程邊坡地表的破裂程度，定義地表破裂度為

式中，α為地表破裂度，S1為地表的塑性區面積，S2為地表的整體面積。通過對塑性區云圖的處理，將破裂度畫成曲線表示，如圖16所示。

由圖16可知，邊坡坡面的破裂度隨著時步的增加而增加。在變形壓密階段，此時為彈性階段，故沒有塑性區產生。但是在第一次開挖變形穩定時，塑性區開始出現，這就表明巖體開始出現塑性區域，但無裂紋產生，破裂度達到0.36。進入到破壞階段，破裂度再次上升，達到0.5，這時邊坡已經破壞一半，但是裂紋并沒有貫穿邊坡。隨后進入到第二次開挖，破裂度再次攀升，達到0.8，這時坡面已經幾乎完全破壞，坡面也已經被裂紋基本貫穿，此時邊坡非常危險。再次開挖，當裂紋完全貫穿，就有很大概率發生失穩滑動，進入到第三次開挖后，破裂度到達0.9以上，邊坡發生失穩。

4 結論

（1）南幫邊坡在變形壓密階段到第三次開挖階段前，均屬于負反饋系統，屬于穩定滑動。第三次開挖階段之后，頂部裂紋繼續增多，但不發生大規模滑移，但中部弱層和底部區域發生失穩滑動，中部和底部屬于正反饋系統，滑體沿著弱層面向下滑移，最后推擠在坡底，坡底受到上部的壓力，再加上剪切裂紋的存在和北幫的阻礙，使得坡底形成底鼓，南幫的滑坡模式為滑移—剪切式。

（2）變形階段，邊坡整體向下壓密變形，越靠近坡面，x方向和y方向的變形位移越大；破壞階段，邊坡頂部受拉，萌生拉張裂紋，沿著弱層面擴展，坡底受剪，隨著開采深度增加，裂紋數目增多，裂紋與裂紋之間相互連接，坡頂和坡底裂紋逐漸連接，最終貫穿邊坡模型；滑動階段，滑體沿著裂紋面發生失穩滑動，最后在坡底剪出。

（3）隨著開挖深度增加，南幫弱層出露越多，邊坡的臨空區域增多，引發滑坡幾率越高，當開挖深度達到標高為-164 m后，南幫邊坡出現大規模滑坡。