改進DB-IWHR模型及其在尾礦庫潰壩影響分析中的應用

沈鴻杰，徐力群，劉子茜

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

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

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

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

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

1 改進DB-IWHR模型

1.1 尾礦庫結構特點

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

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

b.尾礦庫內沉積尾礦料物理力學性質的時空變化規律更為復雜。土石壩壩體通常采用分區設計，同一分區內材料的級配、壓實度和物理力學特性基本一致；而尾礦料在水力沉積作用下尾礦顆粒的空間差異性和沉積尾礦在自重荷載和上覆荷載作用下的固結變形和強度變化產生的時間變異性，使其物理力學特性的時空變化規律與土石壩相比更為復雜[11]。

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

1.2 改進DB-IWHR模型的潰口流量計算

1.2.1潰口流量

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

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

1.2.2潰口沖刷

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

(2)

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

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

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

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

(3)

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

根據曼寧公式，可以得到臨界剪應力與臨界流速的關系：

(4)

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

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

表1 改進DB-IWHR潰口模型參數

1.2.3潰口擴展

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

(5)

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

1.2.4初始條件

潰口流量可以通過單位時間內水庫庫容的損失來確定，根據質量守恒定律，可得水量平衡方程：

(6)

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

在單潰情況下，假設潰壩起始條件為

(7)

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

(8)

(9)

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

1.3 尾礦庫潰壩下泄沙流演進模擬方法

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

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

1.4 改進模型潰口演化模擬精度驗證

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

圖1 唐家山堰塞壩潰壩原始模型與改進模型模擬結果對比

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

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

2 工程應用

2.1 工程概況

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

圖2 尾礦庫剖面及材料分區

2.2 模型建立

2.2.1二維水動力模型

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

圖3 尾礦庫下游影響區地形

2.2.2潰口下泄流量過程

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

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

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

2.2.3模型參數設定及邊界條件

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

圖5 潰壩下泄沙流演進過程

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

2.3 模擬結果與敏感性分析

2.3.1潰壩下泄沙流演進過程

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

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

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

圖6 各敏感點淹沒深度及流速隨演進時間變化曲線

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

2.3.3改進模型敏感性分析

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

圖7 參數敏感性分析曲線

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

3 結 論

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

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

c.對改進模型的初始潰口寬度、沖刷侵蝕參數進行敏感性分析，結果表明初始潰口寬度對流量過程影響不大，但模型對于沖刷侵蝕參數的敏感性較強。