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

周 濤 李啟航 黃宜超

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

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

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

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

1 離散元接觸模型

1.1 單元間的接觸模型

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

Fn=KnXn

(1)

Fn=KsXs

(2)

Fsmax=Fs0-μpFn

(3)

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

1.2 兩個(gè)不同單元的連接

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

(4)

對(duì)于切向剛度為ks1和ks2的兩個(gè)單元，等效切向剛度Ks為：

(5)

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

2 尾礦壩模型建立

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

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

2.1 隨機(jī)堆積模型

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

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

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

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

3.1 滲透力與重力

實(shí)際工程中很難對(duì)每個(gè)孔隙的滲流情況進(jìn)行研究，也是沒(méi)有必要的，設(shè)定浸潤(rùn)線如圖1d)所示,浸潤(rùn)線以上為非飽和區(qū)域，浸潤(rùn)線以下為飽和區(qū)域。研究壩體滲流，不考慮滲流在土壤孔隙中流動(dòng)路徑的迂回曲折，只考慮滲流的主要流向，并認(rèn)為全部的滲流空間均被水所充滿(mǎn)。在模擬飽和砂滲流主要作以下假設(shè)：1)壩體內(nèi)部滲流為層流，認(rèn)為壩內(nèi)符合達(dá)西定律。2)壩體內(nèi)部滲流為漸變流(杜平假定)，認(rèn)為滲流中任意過(guò)水?dāng)嗝娴乃搅魉俸捅冉刀际窍嗟鹊摹?/p>

基于以上假設(shè)通過(guò)自編程序?yàn)槊總€(gè)浸潤(rùn)線以下的顆粒賦予滲透力。根據(jù)達(dá)西定律，在水平與垂直方向賦予顆粒體力以合成與浸潤(rùn)線方程切線方向一致、大小符合達(dá)西定律的滲透力。根據(jù)每個(gè)顆粒的坐標(biāo)，實(shí)現(xiàn)給每個(gè)顆粒豎直方向賦予垂直向上的力以模擬浸潤(rùn)線下的浮力。

3.2 飽和砂性能的模擬

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

4 數(shù)值模擬結(jié)果

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

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

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

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

5 結(jié)語(yǔ)

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