沈鴻杰，徐力群，劉子茜

(河海大學(xué)水利水電學(xué)院，江蘇 南京 210098)

尾礦庫作為一種蘊(yùn)含著很高勢能的人工泥石流危險(xiǎn)源，當(dāng)發(fā)生潰壩事故時，大量下泄尾砂在極短時間內(nèi)沖向下游，威脅下游居民和生態(tài)環(huán)境[1]。據(jù)統(tǒng)計(jì)， 2001—2015年我國共發(fā)生尾礦庫潰壩事故65起，其中洪水漫頂占事故總數(shù)的36.9%[2]。洪水漫頂是導(dǎo)致尾礦庫發(fā)生潰壩的主要原因之一，較為準(zhǔn)確地預(yù)測其潰壩過程及影響范圍，是制定應(yīng)急預(yù)案的基本條件[3]。

尾礦庫的潰壩過程類似于散粒體土石壩及堰塞壩的潰壩，尾礦庫潰壩計(jì)算理論的研究多參考水庫潰壩實(shí)例[4]。近年來，不少學(xué)者基于尾礦庫壩體材料的特殊性開展了相關(guān)的尾礦庫潰壩模型研究。袁兵等[5]根據(jù)多個大壩潰決實(shí)測資料，建立了尾礦壩潰壩數(shù)學(xué)模型，提出了尾礦庫潰壩沙流對下游影響的預(yù)測方法。劉磊等[6]通過建立尾礦壩潰決物理模型對潰壩過程進(jìn)行觀測，依據(jù)試驗(yàn)結(jié)果建立了尾礦庫漫頂潰壩洪水預(yù)測數(shù)學(xué)模型。武立功等[7]采用不同粒徑的尾砂作為筑壩材料，研究了尾砂粒徑對尾礦庫漫頂潰壩的影響。

基于唐家山堰塞湖實(shí)測潰決資料，Chen等[8]開發(fā)了DB-IWHR模型，提出了以雙曲線形式代替指數(shù)形式的沖刷模型，引入了簡化Bishop法模擬潰口擴(kuò)展過程，并采用速度增量形式直接求解控制方程，避免了數(shù)值迭代。張強(qiáng)[9]運(yùn)用DB-IWHR模型對“10·10”和“11·3”白格堰塞湖潰決過程進(jìn)行了反演分析，結(jié)果表明模型可較為準(zhǔn)確地計(jì)算洪水演進(jìn)過程。張貴金等[10]采用DB-IWHR模型計(jì)算堆石壩潰壩洪水過程，利用MIKE21模擬洪水演進(jìn)過程，評估了極端條件下堆石壩潰壩對下游的淹沒風(fēng)險(xiǎn)。

雖然尾礦庫潰壩過程與土石壩相似，但尾礦庫建造過程和條件與土石壩不同，目前的潰壩潰口數(shù)值分析模型中很少考慮到尾礦庫的自身特點(diǎn)，直接采用現(xiàn)有土石壩模型進(jìn)行潰壩洪水演進(jìn)計(jì)算，會產(chǎn)生較大誤差。此外，據(jù)尾礦庫潰壩調(diào)研實(shí)際情況，大量尾礦庫發(fā)生潰壩后，總泄沙量約為實(shí)際庫容的1/3，極少出現(xiàn)全部庫容尾砂下泄的情況。為此，本文根據(jù)尾礦庫結(jié)構(gòu)和泄沙特點(diǎn)，提出了一種改進(jìn)模型，以合理模擬尾礦庫潰壩過程，準(zhǔn)確預(yù)測尾礦庫發(fā)生潰壩后對下游的影響，為尾礦庫潰壩的防災(zāi)減災(zāi)提供理論依據(jù)。

1 改進(jìn)DB-IWHR模型

1.1 尾礦庫結(jié)構(gòu)特點(diǎn)

尾礦庫筑壩方法可分為上游式筑壩法、下游式筑壩法、中線筑壩法、高濃度堆積法和水庫式堆積法5種。我國目前修筑的尾礦庫多采用上游式筑壩法，約占金屬礦山尾礦庫的95%。上游式筑壩法修筑的尾礦庫一般由初期壩和多級子壩構(gòu)成，初期壩一般為透水堆石壩，而后期堆積壩主要利用庫內(nèi)尾礦料填筑而成，隨著生產(chǎn)的進(jìn)行，尾礦壩不斷加高形成多級子壩。在結(jié)構(gòu)特點(diǎn)及潰壩分析方面，尾礦壩與土石壩相比具有明顯差別，主要表現(xiàn)在以下幾個方面：

a.尾礦壩的安全性普遍低于土石壩?，F(xiàn)階段我國的尾礦壩多采用上游式筑壩法，其筑壩方法簡單、生產(chǎn)管理方便，但存在抗地震液化能力差、排滲系統(tǒng)易淤堵導(dǎo)致浸潤線抬高等問題，影響尾礦庫安全性。

b.尾礦庫內(nèi)沉積尾礦料物理力學(xué)性質(zhì)的時空變化規(guī)律更為復(fù)雜。土石壩壩體通常采用分區(qū)設(shè)計(jì)，同一分區(qū)內(nèi)材料的級配、壓實(shí)度和物理力學(xué)特性基本一致；而尾礦料在水力沉積作用下尾礦顆粒的空間差異性和沉積尾礦在自重荷載和上覆荷載作用下的固結(jié)變形和強(qiáng)度變化產(chǎn)生的時間變異性，使其物理力學(xué)特性的時空變化規(guī)律與土石壩相比更為復(fù)雜[11]。

c.尾礦庫潰壩下泄流體性質(zhì)更為復(fù)雜。土石壩潰決后下泄至河道的流體主要為庫水，尾礦庫潰決后主要下泄物為水砂混合物，其流變性質(zhì)與泥石流具有一定相似性，目前對尾砂漿體流變性的研究很少，且尾砂漿體的流體參數(shù)選取和數(shù)值模型沒有統(tǒng)一的標(biāo)準(zhǔn)[12]。

1.2 改進(jìn)DB-IWHR模型的潰口流量計(jì)算

1.2.1潰口流量

潰口流量按寬頂堰公式[8]進(jìn)行計(jì)算，即：

式中：Q為潰口流量，m3/s；C為綜合流量系數(shù)，m1/2/s；B為潰口寬度，m；z為潰口底高程，m；H為水庫水位，m；mb為寬頂堰流量系數(shù)；mq為側(cè)收縮系數(shù)；g為重力加速度，m/s2。

1.2.2潰口沖刷

潰口沖刷數(shù)值分析采用雙曲線模型[8]，公式如下：

(2)

其中τ′=k(τ-τc)

式中：r為侵蝕率，mm/s；τ為剪應(yīng)力，Pa；a、b為沖刷侵蝕參數(shù)；k為在剪應(yīng)力范圍內(nèi)允許r接近其極值的單位變換因子；τc為臨界剪應(yīng)力，Pa。

臨界剪應(yīng)力是判斷潰壩潰口開始與終止發(fā)展的重要參數(shù)，DB-IWHR模型中一個重要假定就是壩體組成物質(zhì)較均勻。對于上游式筑壩法修筑的尾礦庫，初期透水壩一般為碾壓堆石壩，在現(xiàn)代施工條件下，孔隙率較小，穩(wěn)定性較好，不易被水流沖刷破壞，而尾礦堆積壩主要利用庫內(nèi)尾礦料填筑而成，粒徑較小，臨界剪應(yīng)力也較小，遇洪水漫頂時易被沖刷發(fā)生潰壩。因此式(2)中τc的取值不符合尾礦壩實(shí)際情況。

為使模型更適用于尾礦庫漫頂潰壩的實(shí)際情況，對DB-IWHR模型進(jìn)行改進(jìn)，即在不完全潰壩模式下，根據(jù)尾礦壩剖面進(jìn)行材料分區(qū)，并采用修正Shields曲線[13]計(jì)算臨界剪應(yīng)力τc。當(dāng)τ>τc時，潰口開始沖蝕；當(dāng)τ<τc時，潰口完全發(fā)展，沖蝕結(jié)束。臨界剪應(yīng)力τc的計(jì)算公式為

(3)

式中：D*為無量綱泥沙尺寸參數(shù)；Re*為顆粒雷諾數(shù)；θc為希爾茲數(shù)；γs為尾礦砂容重，N/m3；γ為尾礦砂和水混合物容重，N/m3；D為尾礦顆粒粒徑，m；υ為水的運(yùn)動黏度，m2/s。

根據(jù)曼寧公式，可以得到臨界剪應(yīng)力與臨界流速的關(guān)系：

(4)

式中：R為水力半徑，m；J為水力坡度；h為水深，m；vc臨界流速, m/s；n為曼寧系數(shù)，s/m1/3。

修正的Shields曲線是由錢寧等[14]在原始Shields曲線基礎(chǔ)上，在Re*<2的層流區(qū)域及Re*很大的紊流區(qū)域內(nèi)對θc的取值進(jìn)行修正所得。修正的Shields曲線根據(jù)流體和顆粒的特性引入了無量綱泥沙尺寸參數(shù)D*，避免了數(shù)值迭代，簡化了求解過程。

表1 改進(jìn)DB-IWHR潰口模型參數(shù)

1.2.3潰口擴(kuò)展

潰口的邊坡穩(wěn)定性分析采用簡化的畢肖普圓弧滑動面分析方法[8]。在實(shí)際應(yīng)用中將計(jì)算所得圖形假設(shè)為一系列梯形，采用雙曲線模型逐級計(jì)算潰口寬度和傾角[4]，潰口寬度計(jì)算公式為

(5)

式中：B0為潰口初始寬度，m；Δz為潰口底高程增量，m；β為潰口邊坡傾角，(°)。

1.2.4初始條件

潰口流量可以通過單位時間內(nèi)水庫庫容的損失來確定，根據(jù)質(zhì)量守恒定律，可得水量平衡方程：

(6)

式中：W為庫容，m3；t為時間，s；q為入庫流量，m3/s。

在單潰情況下，假設(shè)潰壩起始條件為

(7)

可得初始潰口寬度及初始潰口底高程，計(jì)算公式為

(8)

(9)

式中：z0為初始潰口底高程，m；H0為初始庫水位，m；q0為天然入流量，m3/s；m為跌落系數(shù)。

1.3 尾礦庫潰壩下泄沙流演進(jìn)模擬方法

通過上述模型可以模擬尾礦庫潰壩潰口發(fā)展情況，預(yù)測下泄流量過程，但無法定量分析潰壩對庫區(qū)下游的影響。因此需要根據(jù)下游庫區(qū)地形，選取合適的計(jì)算模型，模擬尾礦庫潰壩下泄沙流的演進(jìn)過程，為尾礦庫風(fēng)險(xiǎn)預(yù)警提供數(shù)值參考依據(jù)。

本文基于分段模型法模擬尾礦庫漫頂潰壩下泄沙流演進(jìn)過程，具體步驟為： ①根據(jù)尾礦庫相關(guān)參數(shù)，利用改進(jìn)DB-IWHR模型計(jì)算尾礦庫漫頂潰壩潰口演化情況和下泄流量過程線；②基于尾礦庫下游地形資料，確定模型范圍，建立非結(jié)構(gòu)化網(wǎng)格，并生成尾礦庫下游影響區(qū)地形模型；③設(shè)置邊界位置、初始條件、糙率等相關(guān)參數(shù)，將計(jì)算得到的下泄流量過程線作為二維水動力模型的入口邊界條件；④通過二維水動力模型計(jì)算求解，模擬尾礦庫潰壩下泄沙流向下游的演進(jìn)情況，并對模擬結(jié)果進(jìn)行分析，得到尾礦庫漫頂潰壩對下游的影響，評估尾礦庫的潰壩風(fēng)險(xiǎn)。

1.4 改進(jìn)模型潰口演化模擬精度驗(yàn)證

以唐家山堰塞壩的泄流過程模擬為例驗(yàn)證改進(jìn)DB-IWHR潰口模型的合理性，模型的參數(shù)選擇參考相關(guān)資料[15]，具體參數(shù)如表1所示，原始模型和改進(jìn)模型模擬結(jié)果對比如圖1所示。

圖1 唐家山堰塞壩潰壩原始模型與改進(jìn)模型模擬結(jié)果對比

通過下泄流量過程(圖1(a))可以得出，改進(jìn)模型與原始模型對洪峰流量和洪峰流量到達(dá)時間的模擬結(jié)果差別不大，改進(jìn)模型模擬的洪峰流量為7 423.33 m3/s，洪峰到達(dá)時間為5.7 h。在到達(dá)洪峰流量之后，原始模型模擬的流量過程線隨著時間變化與實(shí)測值誤差逐漸變大，從潰壩開始至14 h時，實(shí)測總下泄量為12 434萬m3，原始模型模擬的總下泄量為20 565萬m3，比實(shí)測值大65%，而改進(jìn)模型模擬的總下泄量為13 891萬m3，僅比實(shí)測值大12%，因此，模擬的下泄流量過程更為合理。通過庫水位過程(圖1(a))可以得出，對于實(shí)測數(shù)據(jù)，庫水位從742.5 m下降至719.48 m，原始模型庫水位模擬值下降至698.76 m，遠(yuǎn)小于實(shí)測值，改進(jìn)模型庫水位模擬值下降至719.10 m，誤差為0.5%，模擬精度明顯提高。

通過潰口底高程和潰口寬度變化過程(圖1(b))可以得出，原始模型模擬的潰口底高程從740 m下降至700 m，屬于完全潰壩，與實(shí)際明顯不符。改進(jìn)模型通過控制最終潰口底高程，在底高程下降至719 m時停止縱向沖刷。由于最終的潰口寬度是根據(jù)邊坡穩(wěn)定計(jì)算得出，因此原始模型和改進(jìn)模型模擬結(jié)果較為一致，與實(shí)測數(shù)據(jù)基本符合。當(dāng)潰口發(fā)展到一定深度時，其縱向下切和橫向擴(kuò)展均已停止，潰口不再發(fā)展，因此導(dǎo)致原始模型模擬值與實(shí)測值誤差較大。而改進(jìn)模型通過修正臨界剪應(yīng)力控制最終潰口形態(tài)，模擬結(jié)果更符合實(shí)際情況，采用改進(jìn)模型進(jìn)行潰壩沙流演進(jìn)模擬能更準(zhǔn)確地預(yù)測尾礦庫漫頂潰壩對下游的影響。

2 工程應(yīng)用

2.1 工程概況

某擬建尾礦庫位于湖北省宜都市丘陵地帶，庫區(qū)東、西、南三面環(huán)山，在北部地形狹窄處修建尾礦庫。初期壩為碾壓堆石壩，壩高58 m，壩頂高程340 m；堆積壩采用上游式筑壩法，最終堆積標(biāo)高為360 m，總壩高為78 m，壩頂軸線長約422 m，設(shè)計(jì)最終總庫容為1 222.4萬m3，尾礦庫等級為三等。堆積壩分為尾粉砂和尾粉土兩個區(qū)，尾礦庫材料分區(qū)見圖2。

圖2 尾礦庫剖面及材料分區(qū)

2.2 模型建立

2.2.1二維水動力模型

根據(jù)尾礦庫下游周邊環(huán)境情況確定模型模擬的范圍，依據(jù)地勢下降最快方向，選取自初期壩壩址向下游延伸2 km作為模型范圍。在初期壩下游165 m、660 m、1 000 m、1 275 m、1 555 m處共有5處居民區(qū)，分別定義為居民區(qū)1、居民區(qū)2、居民區(qū)3、居民區(qū)4和居民區(qū)5，居民區(qū)建筑物最低高程分別為239.4 m、207.2 m、194.1 m、189.5 m、和190.8 m。用模擬范圍內(nèi)的等高線地形圖對模擬區(qū)域進(jìn)行網(wǎng)格劃分和地形插值，本次建模依據(jù)為該尾礦庫下游溝谷真實(shí)地形，按1∶1比例建立地形模型，如圖3所示。

圖3 尾礦庫下游影響區(qū)地形

2.2.2潰口下泄流量過程

采用湖北省某尾礦庫設(shè)計(jì)資料作為潰口流量過程模擬的輸入?yún)?shù)，應(yīng)用改進(jìn)DB-IWHR潰口模型對潰口演化過程進(jìn)行分析，將得到的潰口流量變化過程作為后續(xù)模擬的邊界條件，根據(jù)GB 50863―2013《尾礦設(shè)施設(shè)計(jì)規(guī)范》及設(shè)計(jì)報(bào)告確定沖刷侵蝕參數(shù)。由于構(gòu)筑尾礦壩的尾礦砂粒徑較小，計(jì)算得到的臨界剪應(yīng)力遠(yuǎn)小于唐家山堰塞湖，因此尾礦庫的臨界流速小于唐家山堰塞湖模擬中的取值。模型具體參數(shù)見表1。

計(jì)算得出的潰壩水位和下泄流量過程如圖4所示。模擬初期潰口侵蝕較慢，潰口流量增長較為緩慢，水位變化?。患s15 min時，潰口發(fā)生崩塌，潰口流量迅速增大，直至29.5 min時，潰口流量達(dá)到峰值1 420 m3/s；隨后下泄流量逐漸減小， 60 min時，水位穩(wěn)定在340 m，潰口發(fā)展達(dá)到穩(wěn)定狀態(tài)，潰壩過程結(jié)束。整個潰壩過程共60 min，峰值發(fā)生在29.5 min。潰壩過程主要分為3個階段，即初始緩慢泄流、流量快速增長至峰值、下泄沙流逐漸落峰。

圖4 潰壩下泄流量過程及庫水位過程模擬結(jié)果

2.2.3模型參數(shù)設(shè)定及邊界條件

根據(jù)潰壩下泄流量過程確定模型模擬時間為120 min，重點(diǎn)分析潰壩下泄沙流對下游的影響。尾礦庫下游主要為農(nóng)田及灌木，根據(jù)天然河道糙率表，取平均糙率0.04。由于下游無河道，因此將初始水深設(shè)置為0。

圖5 潰壩下泄沙流演進(jìn)過程

邊界條件包括入口邊界和出口邊界，選取初期壩軸線位置作為模型的入口邊界，采用模擬得到的下泄流量過程線作為邊界條件；選取尾礦壩下游2 km處作為模型出口邊界并設(shè)置為自由出流。

2.3 模擬結(jié)果與敏感性分析

2.3.1潰壩下泄沙流演進(jìn)過程

不同時刻下泄沙流的淹沒范圍和淹沒深度如圖5所示。潰壩開始初期潰口流量較小，淹沒深度及淹沒范圍較小，在29.5 min時潰口流量達(dá)到峰值，潰壩沙流開始快速向下游演進(jìn)并迅速堆積，在41.5 min時下泄尾砂演進(jìn)至下游邊界處，由于潰口流量的下降及下泄沙流從下游出口不斷流出，潰壩淹沒深度逐漸下降。

2.3.2潰壩下泄沙流對下游影響分析

根據(jù)潰壩下泄尾砂演進(jìn)路徑，繪制各敏感點(diǎn)淹沒深度隨潰壩演進(jìn)的時間變化曲線，如圖6所示。從下游各敏感點(diǎn)淹沒深度隨演進(jìn)時間的變化曲線來看，隨著演進(jìn)時間的增加，淹沒深度逐漸增大，其間某一時刻流量達(dá)到最大，淹沒深度達(dá)到峰值之后逐漸減小，最后趨于穩(wěn)定。

圖6 各敏感點(diǎn)淹沒深度及流速隨演進(jìn)時間變化曲線

從潰壩沙流演進(jìn)淹沒模擬結(jié)果來看，潰壩沙流在壩址處流量大，流速快，淹沒深度大，破壞性強(qiáng)；距離壩址越遠(yuǎn)，潰壩沙流到達(dá)時間越晚，最大淹沒深度越小，最大流速也越小。對于距初期壩165 m的居民區(qū)1，其淹沒深度最大達(dá)到9.21 m，而對于距初期壩較遠(yuǎn)的居民區(qū)5，其最大淹沒深度降低至2.77 m，未達(dá)到建筑物最低高程。下泄尾砂大多堆積在壩址附近，對1 km外的主要建筑物影響較小，符合沙流演進(jìn)的一般規(guī)律。

2.3.3改進(jìn)模型敏感性分析

考慮潰口初始寬度的影響，分別取B0=3.29 m、6.58 m、13.16 m；考慮沖刷侵蝕參數(shù)的影響，分別取a=1.0、b=0.000 1，a=1.1、b=0.000 3，a=1.1、b=0.000 7，根據(jù)下泄流量過程曲線判斷模型對其的敏感性。參數(shù)敏感性分析曲線如圖7所示。

圖7 參數(shù)敏感性分析曲線

由圖7可知，初始潰口寬度對流量過程影響有限，隨著初始潰口寬度的變大，潰壩初期流量較大，但峰值流量較小，且峰值到達(dá)時間較晚。沖刷侵蝕參數(shù)對流量過程影響很大，這是由于沖刷模型是該模型的重要組成部分，而沖刷侵蝕參數(shù)直接影響沖蝕速率，所以模型對于沖刷侵蝕參數(shù)的敏感性較強(qiáng)。

3 結(jié) 論

a.采用修正的Shields曲線改進(jìn)DB-IWHR模型中臨界剪應(yīng)力的計(jì)算，提出了模擬尾礦庫漫頂潰壩的改進(jìn)DB-IWHR模型。利用唐家山堰塞湖潰決洪水過程對改進(jìn)的模型進(jìn)行驗(yàn)證，結(jié)果表明改進(jìn)DB-IWHR模型的潰口演化及下泄過程模擬結(jié)果更符合實(shí)測數(shù)據(jù)。

b.采用改進(jìn)模型結(jié)合二維水動力模型對湖北省某尾礦庫進(jìn)行洪水漫頂潰壩影響數(shù)值模擬，結(jié)果符合尾礦砂流演進(jìn)的一般規(guī)律，得到的下游淹沒范圍、下游敏感點(diǎn)淹沒深度、淹沒出現(xiàn)時間等關(guān)鍵數(shù)據(jù)可為防災(zāi)減災(zāi)提供理論支撐。

c.對改進(jìn)模型的初始潰口寬度、沖刷侵蝕參數(shù)進(jìn)行敏感性分析，結(jié)果表明初始潰口寬度對流量過程影響不大，但模型對于沖刷侵蝕參數(shù)的敏感性較強(qiáng)。