基于MatDEM的尾礦壩壩體失穩(wěn)離散元數(shù)值模擬

周 濤 李啟航 黃宜超

(江西理工大學資源與環(huán)境工程學院，江西 贛州 314000)

尾礦庫是金屬或非金屬礦山尾礦、工業(yè)廢渣的儲存場所，通常通過筑壩或者圍地構成。礦山的安全生產(chǎn)及周邊居民的生命財產(chǎn)安全與尾礦庫壩體的穩(wěn)定性密切相關。尾礦壩的安全在世界各國的礦山建設過程中都受到密切關注，各國都十分重視尾礦壩的建設與管理。

長期以來專家學者利用室內相似材料模擬實驗和數(shù)值模擬技術對尾礦壩開展了廣泛而深入的研究。尾礦壩潰壩的失穩(wěn)模式可分為洪水漫壩失事和壩體結構失事，壩體結構失事模式又可分為滲透破壞失事模式、壩體失穩(wěn)模式和地震險情失穩(wěn)模式[1]。張興凱和陶東良都采用了相似物理實驗法分析了尾礦庫漫頂潰壩的全過程[2,3]。李旭[4]研究了降雨誘發(fā)的尾礦壩潰壩，利用Slide軟件對降雨下的邊坡進行定量分析。敬小非等[5]基于云南拉拉銅礦小打鵝尾礦庫工程設計資料,設計了相似材料模型實驗，模擬了洪水作用下尾礦壩的垮塌和潰決機制，結果表明：尾礦堆積壩的浸潤線變化存在滯后性；在水位上升過程中，壩坡中部出水平方向的總應力增量較垂直總應力增大；洪水導致尾礦壩所受的滲透力、孔隙水壓力、重力增大，削弱壩體材料的抗剪強度，加大自身荷載，從而發(fā)生結構失穩(wěn)。

尹光志等[6]指出尾礦壩組成材料有明顯的顆粒物質特征，在數(shù)值模擬中，相對于連續(xù)介質方法，采用離散元法探索尾礦顆粒組成的結構體的破壞特征以及規(guī)律更加合適。眾所周知在數(shù)值模擬中，正確的計算參數(shù)對最終結果的正確性起著決定性作用。在尾礦庫潰壩數(shù)值模擬中，張千貴等[7]基于離散元軟件PFC2D，采用了粘結彈性接觸模型，通過不斷調整細觀參數(shù)并進行雙軸數(shù)值模擬實驗來擬合實際室內實驗的曲線，從而確定最終的細觀參數(shù)值。劉春等[8]利用其開發(fā)的MatDEM，進行了茂縣滑坡模擬，在實測野外巖土體宏觀力學參數(shù)后，輸入軟件進行自動材料訓練，軟件自動給出最符合宏觀力學性質的細觀參數(shù)。利用MatDEM的自動訓練功能可以免去繁瑣的試算測試過程，可以節(jié)省大量時間精力。王學良等[9]在連續(xù)介質軟件Massfloss中根據(jù)前人對滑面摩擦參數(shù)確定了相應參數(shù)，同時在PFC3D考慮尾礦壩為飽和顆粒的極端情況，在軟件中將離散元顆粒視為彈性接觸模型，不考慮其粘聚力的影響。張曉宇等[10]基于離散元軟件MatDEM模擬了黏性土失水開裂過程，在根據(jù)參數(shù)完成建模后，利用前人總結的公式實現(xiàn)了對計算參數(shù)的動態(tài)計算，數(shù)值模擬效果十分理想。此方法實際是根據(jù)實驗公式?jīng)Q定材料的參數(shù)，依賴于前人對公式總結的正確性。目前，通過數(shù)值模擬方法物理模型的手段的主要研究對象是洪水漫頂導致的潰壩模擬，對其他潰壩模式的實驗以及數(shù)值模擬研究很少。利用離散元法進行細觀模擬的研究也很少。本文基于離散元軟件MatDEM利用預設參數(shù)建立了尾礦壩離散元模型。在模擬中，模擬了尾礦壩逆流漸進式滑坡的動態(tài)過程。通過對比發(fā)現(xiàn)壩體破壞特征以及顆粒流動特征與前人室內實驗結果較為一致。

1 離散元接觸模型

1.1 單元間的接觸模型

離散元法通過堆積和膠結具有特定力學性質的顆粒來構建巖土體結構。在此基礎上通過時間步迭代算法來進行數(shù)值模擬。本文模擬的尾礦壩潰壩實驗采用最基本的線彈性接觸模型，兩顆粒間靠法向與切向的彈簧來模擬相互作用力。通過以下公式來表征顆粒間的彈簧在法向與切向上的破壞準則。

Fn=KnXn

(1)

Fn=KsXs

(2)

Fsmax=Fs0-μpFn

(3)

其中，Kn為法向剛度；Ks為顆粒間的切向剛度；Fs0為顆粒間的初始抗剪強度；μp為單個顆粒的摩擦系數(shù)；Fn為垂直于切向彈簧上的法向壓力。當式(1)中的Xn超過斷裂位移Xb時兩個顆粒間的法向彈簧斷裂；當單元間的切向力超過式(3)中的Fsmax時兩顆粒間的切向彈簧斷裂，顆粒開始滑動。宏觀表現(xiàn)就是斷開的位置連成了一條弧線巖土體出現(xiàn)滑移面。因此式(1)以及式(3)中的參數(shù)的正確標定對模擬結果起決定性作用。

1.2 兩個不同單元的連接

每個單元都有自身的剛度kn，兩顆粒接觸時實際上是兩根彈簧的串聯(lián)，式(1)中的Kn是兩個單元間的連接剛度，對于法向剛度為kn1和kn2的兩個單元，其等效剛度Kn為：

(4)

對于切向剛度為ks1和ks2的兩個單元，等效切向剛度Ks為：

(5)

MatDEM在完成材料賦值操作后，不同材料的顆粒間的連接剛度會根據(jù)公式自動調整。同一種材料不同粒徑的顆粒也會賦予不同的顆粒剛度，因此顆粒間的連接剛度也不同。

2 尾礦壩模型建立

尾礦壩模型的建立步驟主要有堆積建模、重力沉積、切割地層分組與材料賦值等。邊坡主要包含初期壩、尾礦壩主體等，分別給這兩個賦予表1中材料砂土、水質以完成建模。

表1 建模中采用的兩種土體材料的力學參數(shù)

2.1 隨機堆積模型

堆積建模目的是模擬真實世界的重力沉積過程。如圖1a)所示,建立5.5 m×1.28 m的長方形模型箱，設置分散系數(shù)，在模型箱內生成0.14 cm～0.78 cm之間的隨機單元。堆積模型中生成的不同粒徑大小自身有著不同的切向剛度與法向剛度，不同粒徑顆粒間的連接的剛度也不同。重力沉積的過程中，規(guī)則的排列的單元被賦予了隨機的初速度，單元運動一段時間后到達隨機位置上；單元顆粒在重力的作用下向下沉積以模擬顆粒自然沉積的過程；利用上壓力板產(chǎn)生的荷載來壓實模型以模擬顆粒在重力沉積后經(jīng)歷的壓實作用。

2.2 數(shù)字高程建模與材料賦值

MatDEM可以利用導入的坐標點生成曲線或者曲面并以此為基礎劃分地層，對于任何復雜的幾何模型在確定了各層面的坐標后便可以很方便的構建。如圖1b)所示，曲線描述了不同地層界面的坐標位置。在堆積模型建成的基礎上，將這些曲線的坐標位置導入MatDEM后便可以將堆積模型中多余的隨機顆粒切割去除，效果如圖1c)所示。在完成了地層切割過后，對各地層賦予相應的材料參數(shù)，如圖1c)中，深灰色地層被賦予了水的材料性質以模擬庫中水，黑色地層被賦予了砂土材料。切割完地層并且賦予不同地層材料后，顆粒間的剛度發(fā)生了改變，地層的受力因此發(fā)生改變。再次對模型進行受力平衡，為了防止模型破壞顆粒飛出，需要加強顆粒間連接使其不可斷裂，在重力平衡后再將顆粒間的連接設置為可斷裂的狀態(tài)。至此邊坡模型的建立完成。MatDEM材料的細觀參數(shù)通過宏觀參數(shù)按照公式自動轉化而成，詳細公式見文獻[11]。

3 尾礦壩飽和—非飽和模擬

3.1 滲透力與重力

實際工程中很難對每個孔隙的滲流情況進行研究，也是沒有必要的，設定浸潤線如圖1d)所示,浸潤線以上為非飽和區(qū)域，浸潤線以下為飽和區(qū)域。研究壩體滲流，不考慮滲流在土壤孔隙中流動路徑的迂回曲折，只考慮滲流的主要流向，并認為全部的滲流空間均被水所充滿。在模擬飽和砂滲流主要作以下假設：1)壩體內部滲流為層流，認為壩內符合達西定律。2)壩體內部滲流為漸變流(杜平假定)，認為滲流中任意過水斷面的水平流速和比降都是相等的。

基于以上假設通過自編程序為每個浸潤線以下的顆粒賦予滲透力。根據(jù)達西定律，在水平與垂直方向賦予顆粒體力以合成與浸潤線方程切線方向一致、大小符合達西定律的滲透力。根據(jù)每個顆粒的坐標，實現(xiàn)給每個顆粒豎直方向賦予垂直向上的力以模擬浸潤線下的浮力。

3.2 飽和砂性能的模擬

浸潤線下經(jīng)過庫水的浸泡，浸潤線以下的壩體變?yōu)榱孙柡蜕啊Ｉ巴溜柡秃箢w粒間的連接破壞要比非飽和容易，這是因為顆粒間充滿了液體，顆粒間的距離增大，使得顆粒間的連接更加容易斷裂，砂土的抗張拉能力大幅度下降。王海東等[12]研究了不同飽和度下砂土的抗剪強度，發(fā)現(xiàn)內摩擦角隨含水率的增加呈線性下降；粘聚力隨著飽和度增加經(jīng)過上升達到峰值后呈線性減小。如式(3)，內摩擦系數(shù)與初始抗剪力是影響切向彈簧可承受最大力的重要指標。

4 數(shù)值模擬結果

通過建模完成了如圖1a)所示的二維尾礦壩模型。利用該軟件模擬了浸潤線抬升后并在保持穩(wěn)定的情況下，飽和砂顆粒間有效應力減小，抗剪強度減小，最后在重力作用下發(fā)生剪切破壞與滑坡。記錄了尾礦壩壩體失穩(wěn)的位移，總體沉降以及邊坡破壞特征等結果。在不考慮降雨的情況下，邊坡最終破壞形式呈現(xiàn)典型的逆坡牽引式滑坡(見圖2c))，破壞形式與前人的室內實驗結果基本一致[13]。

如圖2a),圖2b),圖2c)所示分別表示了數(shù)值模擬過程中第5時步、第15時步以及第25時步的位移場，每一個步時在軟件中迭代計算了23 250次，模擬現(xiàn)實時間0.186 s，由于水位已經(jīng)上升至破壞時的水位且沒有考慮砂土力學性質逐漸減弱的過程，即主要研究壩體破壞的過程，因此破壞發(fā)生較快。坡體中不同顏色代表了不同的位移值，具體位移值見各圖的顏色條。如圖2a)所示，隨著壩體開始崩潰，坡角處位移相比坡體后緣與中部明顯顏色更深，這意味著坡腳的位移有著明顯的滑動趨勢，最大水平位移達到1.6 cm左右，壩體的后緣以及中部幾乎沒有移動；而在豎直方向的位移上，壩體后緣部分相比壩體中部與坡腳有著明顯的沉降跡象，最大沉降達到了0.8 cm左右，這是因為后緣部分浸潤線上方非飽和土占比相對于壩體中部和坡腳大，浸潤線以下的砂體在滲流飽和后其有效應力大大降低，在上部非飽和砂的重力作用下有非飽和砂土更容易被向下壓實。

當潰壩模擬進行到第15步時，坡腳的水平移動距離已經(jīng)達到了4 cm以上，壩體坡腳至壩體中部部位都有較大的水平位移，從坡腳開始向中部逐漸減小，而中部以上至壩體后緣部分水平位移相比中部以下不是很大，但坡頂最大水平位移也達到了0.5 cm左右；在垂直位移方向，如圖2b)所示坡頂?shù)某两滴灰埔呀?jīng)達到了2.5 cm左右，沉降從坡頂向坡腳處開始逐漸減小，在壩體中部沉降為1 cm左右，而坡腳附近的沉降幾乎為零。此時壩體已經(jīng)出現(xiàn)了明顯的裂縫，這是因為水平方向互動位移不均導致顆粒間法向更容易斷開，再加上壩體沉降不均顆粒間的切向連接也遭到了破壞。

當潰壩模擬進行到了第25步時，如圖2c)坡腳水平位移超過5.5 cm，壩體中部的位移也增長到了3 cm左右，隨著壩體中部水平位移的增大，原有的兩條裂縫間逐漸出現(xiàn)了一條新的張拉裂縫，而坡頂?shù)乃轿灰七_到了1 cm左右；如圖2c)所示，坡頂處的沉降位移已經(jīng)達到了3.5 cm左右。原有裂縫的塊體沉降位移差別不大為2 cm左右，但該區(qū)域仍然出現(xiàn)了新的裂縫，該裂縫的產(chǎn)生主要是由于水平位移變化不均引起的張拉裂縫。

5 結語

本模擬根據(jù)前人總結的實驗公式為不同飽和程度的砂土體自動實現(xiàn)了參數(shù)的更新。利用離散元軟件Matdem建立了尾礦模型，模擬了尾礦庫沉降滑坡的細觀過程，分析了產(chǎn)生裂隙的成因，得出以下結論:尾礦壩牽引坡腳處水平位移大，逆坡向上水平位移呈遞減趨勢。尾礦壩沉降位移在邊坡頂部最明顯，順坡向下沉降位移逐漸減小。壩頂與壩體中部的沉降差大導致裂縫產(chǎn)生，在沉降差別不大的區(qū)域，裂縫產(chǎn)生的主要原因是水平位移不均所導致的。從整個潰壩的過程來看，壩體潰壩有著明顯的牽引式滑坡特征，壩腳以及壩體中部水平位移明顯大于坡頂位移，牽引了壩頂?shù)囊苿樱两挡痪铀倭肆芽p的壩體顆粒的斷裂，裂縫的出現(xiàn)且變寬主要是由于水平位移不均導致的，但沉降不均同樣能導致裂縫的產(chǎn)生。從合位移圖可知浸潤線上部整體有下滑趨勢，坡腳的位移變形最大如若增加初期黏土壩將會大大阻止壩體破壞的發(fā)生。MatDEM具有建模簡便，計算快速的特點。配合高性能計算機，基于其進一步利用野外邊坡的實際數(shù)字高程數(shù)據(jù)、巖土體強度測試數(shù)據(jù)建立實際工程的離散元模型以研究和預測尾礦壩潰壩災害的影響范圍。同時可以通過本文的方法研究降雨帶來的泥石流災害的危害范圍。