李金羽，王光進(jìn)，崔 博，趙 冰，艾嘯韜

（1.昆明理工大學(xué) 國土資源工程學(xué)院，云南 昆明 650093；2.云南省中-德藍(lán)色礦山與特殊地下空間開發(fā)利用重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650093）

0 引言

隨著使用年限的增加，尾礦庫最初的設(shè)計庫量達(dá)到了飽和，因此大部分尾礦庫都進(jìn)行了擴(kuò)容。擴(kuò)容后的尾礦庫高度更高，所帶來的危險就更大。眾所周知，尾礦庫潰壩事故由于原因多變、機(jī)理復(fù)雜、突發(fā)性強(qiáng)、破壞力大，往往造成慘重人員傷亡、巨額財產(chǎn)損失以及難以修復(fù)的環(huán)境污染[2]。一旦潰壩，將會帶來災(zāi)難性的后果。據(jù)不完全統(tǒng)計，截止到2017年，僅云南一省便有尾礦庫637座，按等別可分為：二等庫10座、三等庫51座、四等庫203座、五等庫373座；按安全度可分為：正常庫608座、病庫29座。其中有“頭頂庫”72座（下游1 km距離內(nèi)有居民或重要設(shè)施的尾礦庫）。從數(shù)據(jù)上看，我國現(xiàn)存尾礦庫的數(shù)量不可小覷。因此，對尾礦庫進(jìn)行潰壩模擬對可能發(fā)生的危險提前做好相應(yīng)的措施，具有十分重要的現(xiàn)實(shí)意義。

目前對于FLO-2D軟件的應(yīng)用大部分局限于泥石流的研究，少數(shù)應(yīng)用于尾礦庫潰壩數(shù)值模擬中。李書娜[3]引入FLO-2D模擬軟件，對區(qū)域尾礦壩潰壩事故進(jìn)行了24 h災(zāi)害發(fā)展情況模擬。模擬過程的實(shí)時動態(tài)展現(xiàn)了研究區(qū)域各處的流動速度及流體深度隨時間變化情況，為尾礦庫潰壩安全評價提供了一種定量的方法。阮德修[4]為了研究潰決型泥石流危險性，探討不同降雨頻率下暴發(fā)泥石流的結(jié)果，利用FLO-2D軟件對華溪溝在不同降雨頻率下潰決與非潰決條件進(jìn)行了數(shù)值模擬。許志發(fā)[5]等研究3種下游河道坡降情況下尾砂流的演進(jìn)和沉積規(guī)律，并引入FLO-2D與室內(nèi)潰壩模型試驗(yàn)進(jìn)行參照。本文引入工程實(shí)例為對象，通過數(shù)值模擬，研究不同地形和粗糙度下尾礦庫潰后尾砂流演進(jìn)規(guī)律，最終堆深和影響范圍。

1 潰壩數(shù)值模擬

1.1 FLO-2D軟件簡介

FLO-2D采用一維變量模式與二維漫地流模式仿真一維渠道及二維漫地流之流況，并計算發(fā)生溢堤時，堤外渠道與堤內(nèi)地表間水流之互動機(jī)制，能有效模擬本計劃區(qū)之淹水情形。其渠道演算方式為利用連續(xù)方程式和動力波模式以有限差分方法進(jìn)行模擬，能模擬流體流動的速度及深度，并推估合理的淹沒面積，作為災(zāi)害影響范圍預(yù)測之工具。其數(shù)值模型所采用的控制方程如式（1）～式（3）[6-7]，流變方程如式（4）～式（6）[8]：

式中：h為流深，m；I為水力坡降，%；μ、ν分別為水平和垂直切向速度，m/s。Sfx及Sfy分別為x和y切向摩擦坡降，%；Sox及Soy分別為x、y切向河道坡降，%，g 為重力加速度，m/s2；Sf、Sy、Sν及Std為摩擦坡降、屈服坡降、黏性坡降和紊流-分散坡降；τy為屈服應(yīng)力，Pa；γm為比重 Nm/m3；K 為層流阻力系數(shù)；η 為黏滯系數(shù)；n 為曼寧系數(shù)；Cν為體積濃度 mol/L；α1、β1為黏滯系數(shù)的參數(shù)；α2、β2為屈服應(yīng)力的參數(shù)。

1.2 尾礦庫模型的建立

（1）尾礦庫概況

工況1尾礦庫計總壩高74 m，堆積壩坡的外坡總坡比為1∶3.6，壩頂平臺標(biāo)高1790 m。尾礦庫匯水面積0.568 km2，庫區(qū)占地約0.3641 km2，全庫容820.16萬m3，有效庫容739.94萬m3，屬三等庫。

工況2尾礦庫設(shè)計總高55 m，壩頂平臺標(biāo)高為1730 m，庫內(nèi)流域內(nèi)總匯水面積約1.690 km2，全庫容為983萬m3，有效庫容836萬m3，屬于三等庫，平均外坡比為 1∶3.6。

利用FLO-2D軟件完成了兩座尾礦庫，在設(shè)計壩頂平臺標(biāo)高時（工況1標(biāo)高為1790 m，工況2標(biāo)高為1730 m）全庫容和三分之一庫容潰砂量的潰壩數(shù)值模擬分析，確定了尾礦庫潰壩對下游村莊的影響范圍和程度，獲取不同情況下潰壩后尾礦壩下游溝谷的淹沒范圍、沖擊范圍。同時還通過修改曼寧系數(shù)改變尾礦庫下游粗糙度進(jìn)行對比比較。

由于套筒內(nèi)灌漿料的本構(gòu)關(guān)系尚無成熟理論模型，該項(xiàng)研究灌漿料采用彈塑性模型，泊松比為0.2。通過定義材料的極限彈性應(yīng)力和楊氏模量來實(shí)現(xiàn)彈性階段的應(yīng)力-應(yīng)變關(guān)系，而非彈性階段的應(yīng)力-應(yīng)變關(guān)系則采用規(guī)范[9]提供的混凝土應(yīng)力-應(yīng)變關(guān)系來確定。為了驗(yàn)證模型參數(shù)的準(zhǔn)確與否，對C80混凝土的模型參數(shù)進(jìn)行了ABAQUS驗(yàn)證。C80混凝土，泊松比為0.2，彈性模量為33 GPa，如表2所示。

（2）建立計算模型

數(shù)字地形模型是FLO-2D模擬的基礎(chǔ)資料。依據(jù)研究區(qū)1∶10000等高線的地形數(shù)據(jù)，在BIGMAP中將所需要的計算區(qū)域下載，再由FLO-2D軟件進(jìn)行識別。模擬前，需要將計算區(qū)域劃分為網(wǎng)格。網(wǎng)格單元的尺寸越小，模型的分辨率會越高，同時模擬運(yùn)算的時間也會隨之增加。因此我們需要在網(wǎng)格單元大小（網(wǎng)格單元數(shù)）、所需的模型分辨率和模擬運(yùn)行時間尋找一個平衡點(diǎn)。在嘗試了25 m×25 m、50 m×50 m和100 m×100 m三種不同的計算網(wǎng)格后，本文最終采用了50 m×50 m的計算網(wǎng)格。按照這種計算，網(wǎng)格工況1共劃分2077個計算網(wǎng)格單元，工況2共劃分1103個計算網(wǎng)格單元。并對網(wǎng)格進(jìn)行高程點(diǎn)的插值計算。模擬地形如圖1～圖2所示。

圖1 工況1建模Fig.1 Modeling of working condition 1

圖2 工況2建模Fig.2 Modeling of working condition 2

（3）計算參數(shù)的選取

模擬所需要的參數(shù)的取值均是由現(xiàn)場測量出的基本物理力學(xué)數(shù)據(jù)，經(jīng)1.1節(jié)所列出的方程計算后得出。

通過查找當(dāng)?shù)厮氖謨?，在模擬時將當(dāng)?shù)?4小時暴雨均值數(shù)據(jù)導(dǎo)入流量過程線中。兩種工況下庫容量設(shè)置相同，以工況1容量為準(zhǔn)，即820.16萬m3。針對工況1在標(biāo)高1790 m，工況2在標(biāo)高1730 m時，全庫容以及三分之一庫容的潰砂量及影響范圍進(jìn)行模擬。相關(guān)物理力學(xué)計算參數(shù)的取值見表1。

表1 物理力學(xué)計算參數(shù)的取值表Tab.1 Valuetableofphysicalandmechanicalcalculationparameters

1.3 潰壩模擬分析

1.3.1 同一尾礦庫不同庫容條件下潰壩模擬

使用Mapper軟件將FLO-2D模擬潰壩后的潰壩尾砂流沖擊高度和沉積深度圖與原始地形圖相結(jié)合，得到曼寧系數(shù)0.04下兩種工況在不同庫容下的尾砂演進(jìn)過程和沖擊高度圖（圖3）。

由圖 3（a）、圖 3（b）對比可得，潰壩尾砂流的影響范圍和淹沒深度，在同一尾礦庫不同庫容條件下存在一定的差異，即隨著庫容的增大，其潰壩后的破壞力就越強(qiáng)，下泄尾砂流的流量和流速也隨著增大，所淹沒的范圍跟沉積深度也會隨著庫容的增加而變大。同樣，對同樣庫容下不同尾礦庫的對比可知，其規(guī)律同樣適用，如圖 3（c）、圖 3（d）所示。

由圖4、圖5可知，工況1潰壩后尾砂主要沉積在上游300～580 m和下游1000～1100 m地勢低洼點(diǎn)，工況2主要沉積在上游380～500 m和下游1300～1400 m地勢平緩區(qū)。工況1沖擊高度與沉積深度全過程差值較大；工況2前中期差值較大，后期較為平緩；尾砂流最大沖擊高度隨著下游河道坡降的抬升而逐漸降低，最深淤積點(diǎn)也隨著往下游發(fā)展。

結(jié)合圖 4 和圖 5，由圖 3（a）、圖 3（c）對比可得，不同地形條件下即使庫容相同，它們的潰后尾砂流沉積位置也大不相同。由于兩種工況下的匯水面積不同，本文使用FLO-2D軟件人為修改了工況2的匯水面積，使得兩種工況下的匯水面積保持一致。從圖中可以看出工況1與工況2的地形走勢均為狹長形，其中工況1的坡降為2.08%，工況2的坡降為4.7%。由于兩種工況下的坡降不同導(dǎo)致工況1的沖擊高度和沉積深度比工況2的小，工況1的淹沒范圍大于工況2。此外，雖然尾砂在地形的拐彎處會出現(xiàn)沉積現(xiàn)象，但是其總體趨勢仍然是沿著溝壑流出。對圖3（b）、圖3（d）類比得出的規(guī)律也驗(yàn)證了這一結(jié)果。

圖3 曼寧系數(shù)0.04下兩種工況在不同庫容下的尾砂沖擊高度圖Fig.3 Impact height diagram of tailings under two working conditions and different storage capacity(Manning coefficient is 0.04)

圖4 工況1全庫容潰砂量時縱剖面的沖擊高度與沉積深度分析結(jié)果（1790 m高程）Fig.4 Working condition 1:Analysis results of impact height and deposition depth from the perspective of longitudinal section in case the total storage capacity is crushed(1,790 m elevation)

圖5 工況2全庫容潰砂量時縱剖面的沖擊高度與沉積深度分析結(jié)果（1730 m高程）Fig.5 Working condition 2:Analysis results of impact height and deposition depth from the perspective of longitudinal section in case the total storage capacity is crushed(1,730 m elevation)

1.3.2 同一尾礦庫不同曼寧系數(shù)條件下潰壩

通過修改曼寧系數(shù)改變尾礦庫的下游粗糙程度，得到了曼寧系數(shù)0.04和0.50下兩種工況在不同庫容下的尾砂演進(jìn)過程和沖擊高度圖，見圖6～圖7。

由圖 6（a）、圖 6（b）對比可以得到，當(dāng)曼寧系數(shù)增大時，尾礦庫潰壩影響范圍基本沒有變化，但是潰后下游的尾砂流的沖擊高度和尾沉積深度顯著提高，加大了破壞力，對下游生命財產(chǎn)的威脅更大。在工況2條件下的圖7（a）、圖7（b）對比結(jié)果也證實(shí)了這一規(guī)律。

圖7 工況2在三分之一庫容下不同曼寧系數(shù)下沖擊高度對比圖Fig.7 Working condition 2:Comparison of impact heights under different manning coefficients and at one-third storage capacity

2 結(jié)論

通過以上模擬與計算可以得到以下結(jié)論：

（1）利用FLO-2D軟件對尾礦庫潰壩進(jìn)行建模、計算，能夠直觀的得到尾礦庫潰后淹沒范圍，從而對下游村莊、工廠等建筑物的保護(hù)方法提供了科學(xué)依據(jù)。

（2）庫容對潰壩泥石流演進(jìn)過程起著至關(guān)重要的作用，在同一尾礦庫條件下，庫容越大，對下游的沖擊和影響范圍就越大，所造成的損失也就越大。

（3）不同地形的尾礦庫即使庫容相同，其潰后的尾砂沉積深度和淹沒高度也是不一樣的；相似地形下，庫容一樣的尾礦庫，坡降越大，潰后的尾砂流沉積深度就越深，沖擊高度越高，對下游的威脅就越大。

（4）粗糙度對模擬結(jié)果有一定的影響，曼寧系數(shù)越大，下游粗糙度就越大，對潰后的泥石流阻礙性就越大，造成的尾砂流沉積深度就越大。