模型試驗及數值模擬下尾礦庫潰壩尾砂流演變預測

羅昌泰， 李棟偉， 余國平， 張潮潮， 徐 斌

(1.東華理工大學 土木與建筑工程學院，南昌 330013； 2.南昌工程學院 土木與建筑工程學院，南昌 330099；3.中國瑞林工程技術股份有限公司，南昌 330038)

尾礦庫是用以堆存礦山進行礦石選別后排出的尾礦或其他工業廢渣的場所，是一個具有高勢能的人造泥石流危險源[1]。Lemphers[2]對世界范圍內3 500座尾礦壩進行了統計分析，發現其潰壩失事概率是水庫大壩的10倍以上?；跉v史事件，Davies等[3]發現尾礦庫的潰壩概率在1∶700～1∶10 000。尾礦庫潰決所形成的尾砂流將造成下游居民的傷亡和設施的破壞，并造成環境災害。例如，2008年9月8日山西省襄汾縣新塔礦業發生的特大尾礦庫潰壩事故，造成281人死亡，34人受傷，并造成環境災難[4]；2010年9月21日，紫金礦業尾礦庫潰壩事件，造成22人死亡，房屋倒塌523戶783幢，損壞房屋6 370間[5]。數十年來，這種事故時有發生。隨著尾礦壩事故的頻發，引起了國內外的高度關注。經統計，地震液化、洪水漫頂和滲流是導致尾礦庫潰壩的主要原因，其直接影響是尾礦壩失穩。Rico等[6]分析了全世界147個尾礦庫事故，并提出了尾礦壩失穩的主要原因包括高陡壩(10～45 m)、異常降雨和地震液化，記錄了218次尾礦庫事故，其中極端降雨是主要原因。Yin等[7]通過對玉溪銅礦尾礦庫壩體結構的改變(包括降低壩體高度和增加壩體排水管等工藝)進行物理模型試驗，驗證之后并將該技術工藝應用于實際工程中以加固壩體。

為了預防和減少尾礦庫潰壩的事故，許多國家通過立法來規范尾礦庫的建造和管理，規定在尾礦庫的設計、建設、運行和監測過程中，必須進行嚴格的評價，防止造成尾礦庫的潰壩事故[8]。自20世紀60年代以來，國內外對尾礦庫潰壩的形成機理進行了深入、系統的研究[9]。目前對潰壩問題的研究主要有3種方法：理論分析、數值模擬和物理模型試驗[10]。早期的研究主要以理論分析為主，但理論研究通常局限于比較簡單的情形[11]。隨著科學技術的發展和計算手段的進步，數值模擬是當今主流的研究手段，諸多學者通過對尾礦壩的安全性和環境演變進行數值模擬，表明在尾礦庫潰壩預測中，沖沙距離、泥沙深度和沖擊壓力是關鍵因素[12-13]。物理模型試驗在尾礦庫潰壩研究中常作為預測這3個關鍵因素的經驗方法，與數值模擬方法進行相互的驗證及補充。然而，不同的數值模擬方法采用不同的數學模型，其基本理論、計算公式和方程也不同，比如利用不同的數值計算方法或方案(如：有限元法、有限體積法、概念有限差分方法、離數字高程模型(digital elevation model，DEM)法、光滑粒子流體動力學(smoothed particle hydrodynamics，SPH )方法、淺水波方程法)以及不同的本構模型(如：Bingham，HBP，Voellmy，Cross，and Viscous)來分析尾礦庫復雜的潰決特性，這也導致了潰壩結論存在較大的差異。

潰壩形成尾砂流的過程中流體的黏度是增大的，隨著漿體黏度的增大，它的沖擊力與沖擊力標準差呈非線性增大的[14]，在尾礦庫潰壩研究中，掌握其潰決后尾砂流的流速、沉積深度和移動距離往往是研究的關鍵，這對潰決尾砂流的監測、預警和預防是至關重要的。對諸多尾礦庫潰壩事件的統計分析發現在坡度相對平緩(300‰以上)的沖溝中，即使在暴雨條件下，沖溝也會保持穩定[15]。局部強降雨引起的尾礦庫潰壩是由于暴雨過程中地下水位突然升高，孔隙水壓力過大而引起的，導致溝道內巖土體失穩加劇，滑脫加速。但部分尾礦在暴雨過程中，地表附近并沒有從排水孔中滲出地下水，監測結果也表明地表附近的地下水位不受暴雨影響[16-17]。

為了闡明尾礦庫潰壩尾砂流的運動規律，預測尾砂流的運動距離和作用范圍，眾多學者采用了Ramms、CFD、DAN3D、SPH、Mass-flow等數值模型來計算其運動規律。其中RAMMS軟件使用了雙參數Voellmy摩擦模型來表達流動碎片之間的摩擦運動，研究表明，Voellmy摩擦模型能夠較準確地模擬尾砂流的運動軌跡[18]。RAMMS軟件在標定Voellmy模型過程中通常需要參考記錄良好的歷史事件，并確定可以在后續分析中使用最佳的擬合參數集。此外，RAMMS軟件還可以對DEM或DSM(digital surface model)的結果進行融合，可修改地形數據和增加附加參數的數量，這方便對尾砂流的預測[19]。然而，這些研究尚未得到實際尾礦庫潰壩的驗證，尾礦庫潰壩不同于邊坡泥石流，其流道較長、體積較大、流動深度較大，而且可能攜帶有害物質。因此，有必要開展RAMMS軟件在尾礦庫潰壩中的驗證應用。

本文旨在預測擬建尾礦庫潰壩后的尾砂流對下游的影響情況，給潰壩風險評估提供科學的依據。為了模擬潰壩對下游的影響，如影響范圍或沉積厚度，本研究采用了RAMMS中的Voellmy摩擦模型進行數值計算，并且按照最危險、影響程度最大的工況進行模擬計算，即模擬尾礦壩堆積到最終高程時發生潰壩，目的是充分評估尾礦壩失穩時的最大潛在危險區和受災范圍。同時，為了驗證數值模擬的準確性，利用DEM數據按1∶150的比例搭建了尾礦庫的三維山體物理模型進行潰壩試驗，驗證數值模擬的準確性以及校驗其計算參數。最后，提出尾礦庫潰壩對下游環境的影響程度以及做出準確的風險評估。

1 尾礦庫原型

擬建尾礦庫位于福建省。尾礦庫設計的最大庫容量約為9.3×107m3，滿堆后壩體位置在高程660～825 m。庫區靠近東海，每年6—9月間臺風和暴雨頻繁。因此，在尾礦庫建設前必須充分考慮極端降雨條件下，尾礦庫潰壩對下游所能造成的破壞及最大的影響范圍，對可能造成的破壞后果進行風險評估，并依此制定治理和預防事故的對策措施。

擬建尾礦庫庫區地形包括低丘和丘陵?？傮w地勢西高東低。庫區周圍山坡坡度為25°～45°，山谷和山脊的走向是從西北到東南，山溝海拔為555～855 m。庫區以砂質頁巖、板巖和花崗巖為特征。尾礦庫庫區位于所在溝谷的最上游，下游溝谷蜿蜒曲折，自壩腳至下游3 km(時州村)之間的地形為多處山體收口與開闊地帶交替的“葫蘆串狀”溝谷，溝底落差相對較大。下游8 km(北坑及新前排)開始，溝谷地形完全開闊，地勢趨于平緩，村落分布集中。下游村落主要分布情況如圖1所示，有：上地村(距尾礦壩0.7 km)，時州村(3 km)、老土樓(3.35 km)、門口庵(3.54 km)、時州水庫(5.8 km)、時州水電站(7.12 km)、北坑村(7.98 km)、土樓村(9.62 km)和嵩溪鎮(11 km)。

圖1 尾礦庫及下游村落地形(谷歌地圖)

2 數值模擬及模型試驗原理

2.1 數值模擬計算原理

RAMMS方法的核心是顆粒流平均深度運動方程的有效二階數值解，可在三維模型上計算任意位置、任意時刻的流量的沉積高度和運動速度[20]。

在較大尺度工程問題上，常常是將數值計算模型與地理信息體統(geo-information system,GIS)相結合起來應用[21]。RAMMS軟件首先是建立數字高程模型(DEM)，該文件是數值模擬預處理中最重要的輸入數據，能夠準確反映尾礦庫庫區的重要地形特征。DEM產生過程為：①從尾礦庫區調查地貌數據中提取并轉換DEM數據，得到等高線地形圖；②對潰壩研究區的地形文件、三角網文件和柵格文件進行依次生成；③生成尾礦庫潰壩研究區域的DEM文件。根據壩體高度、尾礦量和附近地形，結合DEM數據，建立尾礦庫三維模型，如圖2所示。

圖2 勘察地貌數據的DEM生產過程

使用GIS繪圖功能可指定單個或多個塊體釋放區域，或者通過輸入水力線來指定作為時間軸的流動，或使用輸入水位曲線來指定流量作為時間的函數，為用戶定義模擬的信息，如流體的釋放區域、流動模式和最終狀態?？梢栽诘匦文Ｐ蜕席B加地圖和遙感圖像來促進模型與歷史事件數據的校準(例如流速、沉積深度和流動距離)。RAMMS仿真計算中摩擦阻力的參數分為干摩擦(μ)和黏滯/湍流系數(ξ)[22]。摩擦阻力S(Pa)的計算公式為

(1)

式中：ρ為密度;g為重力加速度；尾砂運動表面的正應力N為N=ρhgcos(φ)，φ為傾斜角，h為流動高度，U=(ux,uy)T，由x方向和y方向的流速組成[23]。速度的大小關系式為

(2)

式(2)可以用質量平衡方程在Voellmy流變模型中表示為

?tH+?x(HUx)+?y(HUy)=Q(x,y,t)

(3)

式中:H為流動高度；Q(x,y,t)為批量生產源項。如果Q>0，則稱為夾帶速率；當Q=0時，沒有物質侵蝕/沉積；當Q<0時，表示沉積速率。流體的平均深度平衡方程如下

Sgx-Sfx

(4)

Sgy-Sfy

(5)

式中：Cx和Cy分別為剖面系數；gz為垂直重力加速度；Sgx=gxH,Sfx=gyH分別為x方向和y方向的驅動、重力加速度；Sfx，Sfy為驅動摩擦，重力加速度在x和y方向分別表示為

(6)

在Voellmy模型中，垂直方向上的接觸關系可以定義為多相莫爾-庫侖關系，其中nUx為x方向速度方向單位向量，nUy為y方向的速度方向單位矢量，ka/p為土壓力系數，表示為

(7)

式中:φ為內摩擦角；Ka/p為土壓力主動/被動。主動區是指膨脹流動區?·U≥0，被動區是指壓縮區域?·U<0。此外，除了Rankine理論外，Savage等[24]還采用了其他方法來估算土壓力系數。

利用上述方法可推出Voellmy流變公式為

(8)

式中:變量分別用長度L、速度(gl)/2和時間(L/g)/2測量，以獲得統一的流變值;k，μ′和s分別為橫向土壓比、有效動摩擦因數和運動方向的各向異性集合；重力矢量是z=(0,0,-1)，將擾動系數定義ζ=ξ/g為無量綱。

通過對式(8)的整理，可以得到重力巖土體黏性流動中流體內阻的計算公式為

(9)

式中:N0為流動尾礦的屈服應力；μ為一個“硬化”參數。與標準的Mohr-Coulomb類型關系不同，該公式確保當N→0時S→0，當U→0時S為0。它增加了剪切應力，從而導致尾砂流停止較早，這取決于N0的值。

2.2 物理模型試驗原理及模型搭建

2.2.1 試驗原理

物理模型試驗已廣泛應用于尾礦庫潰壩規律的研究，以模擬按特定比例微縮的地貌特征下尾礦庫潰壩的情況。本研究搭建了1∶150比例的試驗模型，并以物理性質與原尾砂相似的爐渣為模型試驗材料。模型中設計布置了梅花形滲透點，用于監測和調節尾礦庫浸潤線的高度，由于研究的目的是分析和預測尾礦庫潰壩破壞的后果。因此需保證尾礦庫的飽和度，即浸潤線高度，以達到尾礦庫潰壩流量。以往的研究表明，滲透潰壩試驗與漫頂潰壩試驗的特性相似[25]。該物理模型雖然簡化了降雨入滲對尾礦庫表面的侵蝕作用，以注水的方式達到漫頂的目的，對尾礦庫潰壩成泥砂流過程中的發生、發展和淹沒范圍進行了真實的論證。同時可驗證RAMMS數值模擬結果的正確性。

尾礦庫物理模型試驗遵循相似性準則，縮制成模型進行試驗研究。以模型重演與原型相似的自然形態，進行觀測，取得數據，然后按照一定的相似性準則引申于原型，以獲得原型的實際現象和性質。模型和原型的物理現象保持相似所必須遵守的規則稱為相似準則。模型應滿足幾何相似、水流運動相似和動力相似。對于研究潰壩尾砂對下游淹沒影響的模型，是以淤積為主的模型，尾砂的揚動現象可以不考慮，尾砂起動相似準則也可以放寬，重點在于滿足尾砂沉降相似、水流挾砂相似及沖淤變形相似3種條件[26]。

(1)尾砂沉降相似條件

根據紊動擴散理論的三元恒定均勻水流中不平衡輸沙狀態下的泥沙擴散方程正態模型中尾砂運動過程沿程沉降的相似條件

(10)

式中：下標ω為尾砂沉速；下標u為水流速度；下標L為長度比尺。

代入Stokes沉降速度比尺關系可得到懸移質沙粒徑比尺

(11)

式中：λd為粒徑比尺；γs為砂容重；γ為水容重。

(2)水流挾砂相似條件

用水流處于飽和狀態時的臨界含砂量來表示水流挾砂力S*。若實際含砂量S≥S*，將沿程發生淤積。水流挾砂力相似，則保證了模型與原型水流含砂量飽和與否相似。此時，水流含砂量比尺(λs)與水流挾砂能力比尺(λs*)相等，即

λs=λs*

(12)

(3)沖淤變形相似條件

由泥沙連續性方程，可得

(13)

根據相似性準則，選取合格的模型砂，原型砂與模型砂的物理性質指標，如表1所示，并計算尾礦庫潰壩模型在水流相似條件下的水力要素相關比尺，表2所示。

表1 原型砂與模型砂的物理性質指標

表2 物理模型試驗相似比尺

2.2.2 試驗模型搭建

物理模型是為了驗證數值模擬的正確性，進而校核數值計算的參數，因此模型搭建范圍可適當選擇關鍵部位即可，大范圍的評價及試驗還需利用數值模擬，可大量節省試驗時間和經費。試驗尾礦庫物理模型采用正態模型，建造比例為1∶150，根據野外調查和DEM數據進行搭建，模型分為六部分：山體地形模型、下游標志地物模型、水量控制系統、浸潤線測量系統、水砂回收系統及監控監測系統，如圖3所示。

圖3 物理模型設置和尺寸

模型尺寸(長×寬×高)為35.0 m×21.0 m×3.3 m，包含從尾礦庫到時州村的實際影響范圍。為了使得模型庫內尾砂飽和，同時形成穩定的滲流，在尾礦庫底部設置了水泵和壓力計等系統和浸潤線控制和測量系統，可對尾礦庫的入滲壓力和滲入量進行測量并能控制浸潤線的高度，并利用一系列的實時監測設備和分析軟件收集尾砂流的波形和速度等數據。

在尾礦庫物理模型中，通過埋入水位監測儀和流量監測點觀測浸潤線的變化規律(見圖3(b)和圖3(c))。采用RIEGL型三維激光掃描儀(3DLS)、無人機(UAV)、5KF20型高清攝像系統(HDCS)對尾礦庫潰壩位移變化進行監測；利用表面流場實時測量系統(VDMS-LAB-04C型)監測潰口流量變化和流量峰值過程。尾礦庫潰壩監測系統布置圖，如圖4所示。

圖4 流場監測系統的布置

3 計算及試驗結果

3.1 潰壩尾砂流演進結果

3.1.1 數值模擬結果

RAMMS模擬與物理模型均是使用笛卡爾坐標系作為數字高程模型快速捕獲地形的數據，如圖5(a)所示?；赩oellmy-Fluid方法，采用無通道尾砂流進尾砂流模擬，將尾礦庫深度、高程和平均坡角定義為RAMMS釋放區域的初始參考信息。由圖5(b)可知，淺灰色“L”形區域為包含關鍵地形信息的計算區域，深灰色區域為潰壩運動的釋放區域，其高度為54 m，體積為1.74×107m3。當各節點元素的動量均小于最大值的5%時，仿真停止。

(a) 尾礦庫DEM模型

RAMMS碎屑流動態分析軟件基于Voellmy流動摩擦理論，假定碎屑流是單向模型，將所有的尾砂作為一個整體，無法區分出固相和液相物質，摩擦阻力由干摩擦參數(μ)和黏滯/湍流型參數(ξ)控制[27]。因此模擬RAMMS模擬尾礦庫潰壩最大的難點在于漫頂潰壩過程中尾礦碎屑流的成分十分復雜，摩擦參數的選擇對結果影響明顯。模型試驗選擇與原型砂力學性質相似、造價低、容重適中的火電廠爐渣作為模型砂。通過現場測量，結合流槽試驗及潰壩模型試驗中尾砂流的流量流速變化，考慮尾砂流的流動路徑和分布，最終確定RAMMS數值模擬的摩擦參數為μ=0.07，ξ=1 700。根據Voellmy方法的內存大小和計算效率，最終生成1 334 514個計算機流量單元和1 338 190個流體節點。通過多次物理模型預備試驗，校驗后的RAMMS數值仿真模擬參數,如表3所示。

表3 模擬計算參數設置及主要計算結果

利用上述方法，對擬建尾礦庫研究區的潛在失效為出發點進行了建模和計算。RAMMS數值結果如圖6所示，RAMMS模擬仿真計算的潰壩進程描述如表4所示。

表4 RAMMS模擬仿真計算的潰壩進程

(a) t=30 s

3.1.2 物理模型試驗結果

試驗工況為洪水漫頂導致尾礦壩潰壩，該工況下尾礦庫庫區水位處于臨界狀態即最高洪水位，此時干灘長度0.67 m(相對應原型100 m)，在該狀態下遭遇洪水最終導致漫頂潰壩發生。本次潰壩模型具有落差顯著，勢能大的特點，潰壩后庫內存水急流直沖而下，可挾帶走大量的尾砂。潰壩尾砂流具體的演進情況如圖7所示，演進的進程描述如表4所示。

(a) t=45 s

模型試驗潰壩進程以及潰壩尾砂流對下游村莊影響的進程如表5所示，其中監測時間為從尾砂流到達壩腳時起算。相對應的原型時間依據時間比尺進行換算(見表1)。

表5 模型試驗潰壩尾砂流進程

物理模型試驗中尾沙總體的潰泄量是采用RIEGL型三維激光掃描儀(3DLS)進行掃描測量,如圖8所示。測得影響范圍為0.558 km2，總潰泄量為1.3×107m3，與RAMMS統計的傾瀉量1.33×107m3對比，二者誤差僅為2.3%。

(a) 潰壩模型中泥砂沉積厚度

在物理模擬試驗中，大壩潰壩和隨之而來的尾砂流發展速度很快。水的注入使得水位上升和尾礦庫坍塌導致礦漿快速移動和膨脹[28]，尾礦壩上的這種快速破壞和泥漿的發展，也體現在了物理模型試驗中決口和漿液的發育是由極端降雨和水分滲透引起的[29]。

從兩者的潰壩進程來看，總體的演進規律相仿，總體的潰泄影響范圍以及潰泄的尾砂量基本一致，存在的差異是在時間節點上。數值仿真模擬計算實現潰泄的方法是設置一個滑坡體下滑過程，而物理模型試驗潰壩的所形成的尾砂流是洪水漫頂沖刷的發展過程，故此存在差異，但不影響試驗結果的準確性，所以在獲取潰泄沖刷的時間時，應以物理模型試驗數據為準，更符合實際的情況。

3.2 淹沒深度與潰泄量對比

為測量潰壩尾砂流沿下游的淹沒范圍和淹沒深度，將研究區5個關鍵位置的淹沒深度與模型試驗結果進行了對比,如圖9所示。

(a) 模型試驗測量截面

在尾礦庫庫區下游關鍵位置設置監測斷面，采用模型試驗實際測量和數值模擬結果分析，對比各監測斷面的尾礦碎屑流的最終流動高度、流量隨時間的變化過程及最終截面處的流量大小。通過模型試驗與數值模擬結果對比分析可知，1#、4#、5#截面處的尾砂流動高度基本一致，上地村區域(2#、3#截面)的尾砂流動高度差別較大，分析是與尾砂的密度有關，模型試驗為了使壩體發生失穩破壞，短時間內持續給庫區加水，會導致尾礦庫在極短的時間內浸潤線高度升高，導致庫區大部分尾砂趨于飽和，但少部分尾砂依然處于非飽和狀態，壩體內部碎屑流的密度及飽和度分布與實際情況略有不同，但不影響試驗結果的準確性。為了下游社區和居民的安全，在風險評估中應以兩者之間的最大值為建議值，對擬建尾礦庫進行安全評價。

將RAMMS數值模擬結果與模型試驗進行對比發現，兩種分析過程相輔相成，相互驗證,如表6所示。兩種研究方案在尾礦庫潰壩動態分析過程中優勢互補，互相驗證。兩者試驗結果與單一方法相比更準確深入，能夠更全面的分析尾礦庫潰壩碎屑流的演進及致災影響問題。

表6 模型試驗與數值模擬結果對比表

4 結 論

本研究采用RAMMS數值模擬計算方法和物理模型試驗手段，仿真模擬了擬建尾礦庫在最終壩高時發生潰壩后尾砂流在下游地形上的運動規律。經過模型試驗與數值計算結果的相互對比與印證，獲得了尾礦庫潰壩破壞時尾砂流的演變進程和潛在的影響范圍等數據。本次尾礦庫潰壩的仿真模擬對工程實踐具有很強的指導意義，同時在試驗過程中也獲得了一些研究經驗。

(1) 研究得出擬建尾礦庫如在最終壩高時發生潰壩，9 min尾砂流開始沖擊距壩腳下游0.7 km處的村莊，18 min將對村莊造成全面淹沒，最終影響范圍約為0.558 km2，尾砂的潰泄總量可達1.3×107m3。試驗結果說明潰壩危害大且基本不具備預警時間，同時結合潛在的危險范圍，建議在尾礦庫的設計、管理和風險評估過程中，合理的增加防災減災的措施。

(2) 通過模型試驗和數值計算的對比，尾砂流潰泄總量二者誤差僅為2.3%，沿線淹沒深度基本一致，說明經模型試驗參數校驗后的RAMMS計算結果是可靠且有效的。潰壩模型試驗成本昂貴、耗時長、難度大且模擬的范圍有限，利用參數校驗后的RAMMS方法則能彌補這些缺陷，同時說明該研究方法能較好地解決尾礦庫潰壩尾砂流穿越真實地形的演進問題。

(3) 尾礦庫潰壩模型試驗復雜繁瑣且難于操控，本次試驗中發現主要影響因素有——地形模型的精確程度、模型砂的相似程度、模型堆砂與原型排砂貼近程度等。本研究中物理模型與數值建模采用相同DEM基礎數據，保證了地形邊界條件的一致性；結合潰壩泥砂流運動特性選擇沉降、水流挾砂和沖淤變形3個關鍵相似條件，確保選配的模型砂運動規律最相似于原型砂；實際庫內排砂是長期的沉積過程，通過調節模型堆砂壓實度滿足與原型相似的力學性質，最大程度貼近于原型排砂。