彭浩，朱光源，張邵祥，蔣水華，蔡木良，黃勁松

(1.南昌大學工程建設學院，江西 南昌 330031;2.國網(wǎng)江西省電力有限公司電力科學研究院，江西 南昌 310058)

尾礦庫是礦山企業(yè)安全生產(chǎn)過程中的重要基礎設施，但同時也是礦山企業(yè)最大的危險源。一旦尾礦庫發(fā)生潰壩事故，高勢能下泄尾砂流將對下游人民生命財產(chǎn)安全和生態(tài)環(huán)境造成巨大影響[1]。洪水漫頂潰壩是尾礦庫較為常見的潰壩模式。據(jù)統(tǒng)計，2001—2013年我國發(fā)生尾礦庫潰壩事故共計65起，其中由洪水漫頂引起的潰壩事故占潰壩事故總數(shù)的36.9%[2]。因此，研究尾礦庫漫頂潰壩的發(fā)展過程，分析下泄尾砂流的演進規(guī)律，科學預測尾礦庫失事后果和制定科學合理的避險轉移方案，具有非常重要的工程應用價值。

隨著計算機技術的日益成熟，數(shù)值模擬方法因成本低、耗時少等優(yōu)勢越來越多地被應用于解決尾礦庫潰壩過程仿真難題。鄭欣[3]利用Fluent軟件模擬了金山尾礦庫潰壩尾砂演進過程，得到了下泄砂流的影響范圍。阮德修等[4]結合Flo2D軟件和3DMine數(shù)字礦山軟件實現(xiàn)了湖南某尾礦庫潰壩過程數(shù)值模擬。Mahdi等[5]發(fā)展了考慮下泄尾砂流非牛頓性質和黏塑性關系的尾礦庫潰壩數(shù)值模擬方法，實現(xiàn)了對加拿大Alberta尾礦庫潰壩尾砂演進過程的合理預測。

盡管這些研究工作為解決尾礦庫潰壩模擬問題提供了有效途徑，但是目前的方法通常將潰壩尾砂流視作物理性質恒定不變的非牛頓流體，并沒有考慮尾礦庫在發(fā)生漫頂潰壩過程中水對尾砂的沖刷作用以及下泄尾砂流含沙量的時變特征。為此，本文依托江西省永平銅礦燕倉尾礦庫5號副壩，利用Civil-3D軟件建立尾礦庫及下游區(qū)域三維幾何模型，采用RNGk-ε湍流模型和泥沙沖刷模型構建尾礦庫漫頂潰壩水動力學數(shù)值模型，對尾礦庫漫頂潰壩過程進行精細化模擬，預測潰壩淹沒范圍及下游村莊的淹沒水深、泥沙淤積厚度、流速等關鍵信息。在此基礎上，制定避險轉移方案，為尾礦庫潰壩風險管理與決策提供理論依據(jù)。

1 尾礦庫漫頂潰壩數(shù)值模擬

本文在Flow-3D軟件平臺上，采用RNGk-ε湍流模型和泥沙沖刷模型構建尾礦庫漫頂潰壩水動力學數(shù)值模型[6-8]。

1.1 水流控制方程

在笛卡爾坐標系中，潰壩水流可視做不可壓縮流體，采用連續(xù)方程和Navier-Stokes動量方程描述流體運動[9]，其中連續(xù)方程為

(1)

Navier-Stokes動量方程[9]為

(2)

式中：u，v，w分別為流體在x，y，z3個方向的流速分量；Ax，Ay，Az分別為在計算單元中各截面流體所占面積與截面總面積的比值；VF為計算單元中的流體體積分數(shù)；p為壓強；ρ為流體密度；gz，gy，gz分別為3個方向的重力加速度；fx，fy，fz分別為3個方向的黏滯力加速度。

1.2 RNG k-ε湍流模型

RNGk-ε湍流模型由Yakhot等[10]于1986年提出。該模型基于嚴格的統(tǒng)計技術，在傳統(tǒng)k-ε湍流模型的基礎上改進而來，能夠有效地解決低雷諾數(shù)和具有強剪切區(qū)的湍流流動問題，計算精度高，適用范圍廣。其控制方程包括紊動能k方程和耗散率ε方程如下：

(3)

(4)

式中：k為紊動能；ε為紊動能耗散率；μ為流體動力黏滯系數(shù)；μt為流體湍動黏度；Gk為紊動能k的產(chǎn)生項；αk，αε，C1ε，C2ε為無量綱系數(shù)，根據(jù)文獻[6]，這些參數(shù)取值為αk=αε=1.39，C1ε=1.42，C2ε=1.68。

1.3 泥沙沖刷模型

泥沙沖刷模型[9,11]通過以下4個方面的計算來實現(xiàn)對泥沙沖刷、淤積和運移的精細化模擬：(1)泥沙起動計算；(2)泥沙沉降計算；(3)懸浮泥沙隨流體運移計算；(4)推移質泥沙運移計算。

泥沙起動指的是水流的拖拽使得沉積在河床上的泥沙被水流帶走，成為懸浮泥沙；沉降指的是水流中的懸浮泥沙在重力作用下沉積。泥沙的挾帶和沉降需依據(jù)水流速度而定，當水流速度大于泥沙起動速度時，泥沙將懸浮在水流中運動；當水流速度小于泥沙沉降速度時，泥沙將逐漸沉積或沿河床推移運動。

1.4 邊界條件

數(shù)值模型所涉及的邊界條件主要包括入流邊界、壁面邊界、出流邊界和對稱邊界四大類。

(1) 入流邊界。主要包括壓力邊界、流速邊界、流量邊界。由于水深和壓力具有相關性，壓力邊界主要約束了邊界處的壓力和水深；流速邊界和流量邊界主要通過約束邊界處的流速和流量來控制入流條件。

(2) 壁面邊界。Flow-3D中常用固壁邊界來模擬壁面約束，在邊界處應用無滑移條件以及垂直于邊界的零速度條件來模擬在流體運動問題中的固體邊壁約束。

(3) 出流邊界。出流邊界包括連續(xù)邊界和自由出流邊界。連續(xù)邊界流速梯度和壓力梯度為0，即流體流經(jīng)邊界時各個物理參數(shù)恒定不變，該邊界對計算結果干擾較大。自由出流邊界假設邊界符合Sommerfeld Radiation條件，流體在流經(jīng)該邊界時的擾動和波動是平順的，較為符合實際情況。

(4) 對稱邊界。對稱邊界條件在邊界處應用零梯度條件以及垂直于邊界的零速度條件。對于具有對稱性的計算模型，可運用對稱邊界減小計算量。另外，若在建模時劃分了多個網(wǎng)格區(qū)域，則網(wǎng)格區(qū)域的邊界重合處需要設置為對稱邊界條件以進行網(wǎng)格區(qū)域間的數(shù)據(jù)交換，實現(xiàn)網(wǎng)格區(qū)域的連接。

2 模型驗證

為了驗證所建立的尾礦庫漫頂潰壩水動力學數(shù)值模型的有效性，首先選取文獻[12]中漫頂潰壩水槽試驗為典型算例進行模型驗證。

2.1 計算模型

該試驗在一長17 m，寬1.2 m的水槽中進行，如圖1所示[12]。槽頂距地面2.4 m，槽深1.5 m，坡度為0。水槽上端為一個5.0 m×5.0 m蓄水池，池深2.4 m。試驗中壩高設置為1.2 m，頂寬0.2 m，上游壩坡為1:1.25，下游壩坡為1:1.2。為保證下游坡面沖刷均勻，試驗開始前在壩頂布置一塊高0.1 m的擋板，之后供水水箱持續(xù)向蓄水池內供水。當壩前水位與擋板上緣平齊時，迅速提起擋板，潰壩開始。潰壩試驗過程中供水水箱的供水流量恒定為20 L·s-1。

圖1 水槽模型試驗平面布置圖[12]

根據(jù)漫頂潰壩水槽試驗建立三維數(shù)值模型并對其進行網(wǎng)格劃分，數(shù)值模型尺寸與室內試驗裝置尺寸相同，如圖2所示。整個壩體由無黏性砂組成，中值粒徑為3.77 mm，密度為2 650 kg·m-3，臨界希爾茲系數(shù)采用Soulsby-Whitehouse公式計算[11]。為了模擬室內潰壩試驗初始條件，上游水深設為1.3 m(高出壩頂0.1 m)。

圖2 水槽模型網(wǎng)格剖分

采用正六面體結構化網(wǎng)格對模型進行網(wǎng)格剖分。根據(jù)該模型幾何特征，設置了兩個網(wǎng)格區(qū)域，這兩個網(wǎng)格區(qū)域的網(wǎng)格尺寸均為0.1 m×0.1 m×0.1m。網(wǎng)格區(qū)域2的上游側設置為流量邊界(Q)，用于模擬試驗過程中的上游來水；網(wǎng)格區(qū)域1和區(qū)域2重合邊界處定義為對稱邊界(S)，實現(xiàn)網(wǎng)格區(qū)域連接；模型上方設置為壓力邊界(P)，模擬大氣壓條件，壓力設置為101.325 kPa；網(wǎng)格區(qū)域1的下游側設置為自由出流邊界(O)；其余部分定義為無滑移固壁邊界(W)，如表1所示。為與文獻[12]中物理試驗時間保持一致，數(shù)值模擬時間也設為250 s。

表1 水槽模型邊界條件

2.2 潰口流量分析

圖3給出了水槽模型試驗和數(shù)值模擬獲得的潰口流量qV隨時間t的變化關系曲線。可知，潰壩發(fā)生以后，潰口流量短時間內迅速增大，到達峰值后開始下降，并逐漸趨近于上游來流流量，這顯然符合漫頂潰壩潰口流量發(fā)展的一般規(guī)律[13]。本文通過數(shù)值模擬得到的各個時刻潰口流量變化規(guī)律與實測值基本吻合，表明本文所建立的尾礦庫漫頂潰壩水動力學數(shù)值模型有效、可行。

t/s

3 工程應用

3.1 工程概況

本文以江西省永平銅礦燕倉尾礦庫5號副壩[7]為例，進行漫頂潰壩尾礦演進過程模擬。該副壩初期壩為均質黏土壩，軸線底部高程為129.5 m，初期壩頂高程為135.0 m，設計外坡比為1:2，設計內坡比為1:2.5，壩頂寬為3 m；堆積壩采用上游法進行分散管放礦堆筑，最終堆積壩頂高程為150 m。正常蓄水位為146.5 m，設計外坡比為1:5，設計干灘坡度為1:167。尾礦壩下游約2 km范圍內有兩個村莊以及多條公路和大面積農(nóng)田，如圖4所示。

圖4 尾礦庫下游衛(wèi)星圖

3.2 潰壩模型

根據(jù)尾礦庫實際尺寸以及下游地形圖，基于地形等高線和高程點信息，采用Civil-3D軟件建立尾礦庫及下游2 km范圍內地形和房屋的實體幾何模型。為提高計算效率，本數(shù)值模擬洪水重現(xiàn)期取規(guī)范上限1 000年，根據(jù)《江西省暴雨洪水查算手冊》推薦方法，計算出千年一遇洪水總量為6.15×105m3，與庫區(qū)內能容納水量(6.19×105m3)相差不大。故可在模擬開始時預先在庫區(qū)內充滿水體，保持庫水位與壩頂高程相等。尾礦庫遭遇洪水漫頂時，壩頂最薄弱部位一般先形成主潰口，但是數(shù)值模擬無法模擬壩頂最薄弱部位。故本文按照文獻[8]的做法，在尾礦壩中部位置預先設置一個寬30 m和深3 m的矩形斷面初始潰口，將漫頂過程逐漸導向最危險工況發(fā)展。

為了分析潰壩尾砂流演進過程及其對下游村莊的影響，在尾礦庫下游設置了若干監(jiān)測點以獲取下游關鍵區(qū)域的淹沒信息。如圖5所示，在新巖前村內設置監(jiān)測點1，其高程為107.12 m；因塘棣源村面積較大，房屋分布較廣，設置了2個監(jiān)測點，分別是位于村口的監(jiān)測點2和位于村內的監(jiān)測點3，其高程分別為90，92.69 m。

圖5 尾礦庫及下游地形三維模型

為了保證數(shù)值模擬精度和效率，將研究區(qū)域劃分為兩個網(wǎng)格區(qū)域，如圖6所示。其計算網(wǎng)格尺寸均為3 m×3 m×3 m，總有效網(wǎng)格數(shù)約為450萬。模型上方設置為壓力邊界(P)，模擬大氣壓條件，壓力值為101.325 kPa；網(wǎng)格區(qū)域1和區(qū)域2重合邊界處定義為對稱邊界(S)，實現(xiàn)網(wǎng)格區(qū)域連接；模型下游出口處設置為自由出流邊界(O)；其余部分有山體阻擋或者為地面，定義為無滑移固壁邊界(W)。具體邊界條件設置見表2。此外，泥沙模型參數(shù)設置如下：平均粒徑為0.112 mm，密度為2 940 kg·m-3，臨界希爾茲系數(shù)為0.05；挾帶系數(shù)為0.01，推移質系數(shù)為8，休止角為34.7°。經(jīng)多次試算確定漫頂潰壩數(shù)值模擬時間為3 200 s。

表2 尾礦庫三維模型邊界條件設定

圖6 三維數(shù)值模型網(wǎng)格剖分

3.3 潰壩尾砂流演進過程及流速

為分析潰壩尾砂流演進過程及流速變化規(guī)律，圖7給出了不同時刻的潰壩尾砂流流速變化云圖。

漫頂發(fā)生后100 s，見圖7(a)，潰口水流通過初始潰口向下游演進并逐漸沖刷壩體，形成的挾砂水流演進至尾礦庫下游約200 m處，潰口流速約為7 m·s-1。潰壩發(fā)生后500 s，見圖7(b)，潰口繼續(xù)擴寬并向上游發(fā)展，潰壩尾砂流已經(jīng)演進至尾礦庫下游約1 000 m處，其流速約為5 m·s-1，即將對塘棣源村造成影響。相比之下，新巖前村地勢較高，只有少量潰壩尾砂流越過山地，進入新巖前村。潰壩發(fā)生后1 000 s，見圖7(c)，在水流的沖刷作用下潰口持續(xù)向上游發(fā)展，此時潰壩尾砂流一部分已經(jīng)進入新巖前村內，村內尾砂流流速在0～3 m·s-1范圍內變化，但是由于新巖前村內建筑物地勢較高，潰壩尾砂流并沒有對其產(chǎn)生影響。另一方面，潰壩尾砂流已經(jīng)演進至塘棣源村，靠近房屋處尾砂流流速約為5 m·s-1。潰壩發(fā)生后1 500 s，見圖7(d)，由于庫水較少，庫內尾砂的沖刷和侵蝕作用較小，因而下泄尾砂流也較少，但是溝谷中和新巖前村內的尾砂流仍然持續(xù)向下游演進，導致下游地勢較低的塘棣源村的尾砂流流速并未明顯減弱，維持在4 m·s-1左右。潰壩發(fā)生后2 000 s，見圖7(e)，庫水已經(jīng)下泄完畢，下游區(qū)域僅剩塘棣源村前還有潰壩尾砂流以較低的流速演進。潰壩發(fā)生后3 200 s，見圖7(f)，尾礦庫下游區(qū)域已無尾砂流演進，大量泥沙已經(jīng)沿程淤積。

(a) 100 s

3.4 淹沒深度

圖8給出了3個監(jiān)測點的淹沒水深h隨時間t的變化關系曲線?？芍瑵挝采傲髟诼敽蠹s750 s到達新巖前村口監(jiān)測點1處，約1 190 s達到峰值，峰值水深為3.42 m；潰壩尾砂流在漫頂后約570 s到達位于塘棣源村口的監(jiān)測點2處，約1 330 s達到峰值，峰值水深為4.56 m；潰壩尾砂流在漫頂后約860 s到達位于塘棣源村內的監(jiān)測點3處，約970 s達到峰值，峰值水深為1.77 m。由于潰壩尾砂流在沿下游溝谷演進過程中遭遇了山體地形阻攔被分為了兩股(見圖7)，其中一股在分岔口跨過地形直接演進至塘棣源村，導致布置于該村的監(jiān)測點2和監(jiān)測點3處的水位上升，隨著該股尾砂流持續(xù)向下游演進，監(jiān)測點2和監(jiān)測點3的水位逐漸達到局部峰值后開始下降；另一股尾砂流越過地形演進至新巖前村，隨后繼續(xù)演進至塘棣源村，到達塘棣源村的時間較晚，這便造成了塘棣源村監(jiān)測點2的水位第2次上升至局部峰值而后再次下降。此外，監(jiān)測點3由于地勢較高，當?shù)?股尾砂流消散后，第2股尾砂流還未到達，故該監(jiān)測點水深在潰壩后1 000 s驟減至0。

t/s

結合布置于村莊內3個監(jiān)測點的高程和峰值水深信息可知，新巖前村內底部高程低于110.5 m的房屋建筑物將會受到潰壩尾砂流影響；而塘棣源村內底部高程低于94.5 m的房屋建筑物將會受到潰壩尾砂流影響。

3.5 最大淹沒范圍及尾砂淤積情況

由上述潰壩尾砂流演進分析可知，新巖前村和塘棣源村出現(xiàn)最大淹沒的時刻分別是在漫頂發(fā)生后1 160，1 280 s，見圖9。潰壩尾砂流進入新巖前村后主要聚集于村內低洼處，但村內房屋總體地勢較高，房屋最低高程為112.12 m，高于110.5 m，故該村房屋均不受潰壩尾砂流的影響。塘棣源村內部分房屋低于94.5 m，存在被尾砂流沖擊和淹沒的可能性，故該村受到潰壩的影響較大。

(a) 新巖前村

下泄尾砂總量約為4.82×105m3，大量淤積于尾礦庫下游溝谷和兩個村莊中，掩埋或損壞村內的房屋、道路等設施。圖10為潰壩結束后下游尾砂淤積云圖。可以清楚看出，潰壩后的泥沙淤積分布情況。尾砂在新巖前村淤積較輕，并且由于該村房屋建筑物地勢較高，淤積尾砂并未對其造成影響。大量尾砂淤積于塘棣源村以及塘棣源村前農(nóng)田處，地勢低洼處尾砂淤積厚度達到10 m左右，部分地勢較低的房屋建筑物將受到較為嚴重的淤積尾砂的影響，甚至存在被尾砂完全掩埋的可能性。

圖10 尾砂淤積云圖

4 避險轉移方案

由上述尾砂流演進分析結果可知，永平銅礦燕倉尾礦庫發(fā)生漫頂潰壩后，新巖前村房屋均不受潰壩尾砂流的影響，而塘棣源村內高程低于94.5 m的房屋建筑物將受到尾砂流沖擊或掩埋，故需要提前制定受災人員的疏散、撤離和避險轉移方案。轉移方案的制定步驟主要如下：(1)根據(jù)區(qū)域地形、交通以及事故影響范圍等信息，選取合適的應急避險場所中心點；(2)根據(jù)需轉移人數(shù)確定應急避險場所，可以按步驟(1)確定的中心點為圓心規(guī)劃一個平面圓形區(qū)域作為應急避險場所[14]；(3)計算受災人員到達應急避險場所的時間。

應急避險場所中心點的選取一方面應遵循安全性的原則，即此地不會被潰壩尾砂流沖擊或掩埋。在能夠保證本身安全的前提下，盡可能考慮能夠容納較多的人口，并遠離災害源，選擇地勢相對較高的地方；另一方面應考慮可達性，即受災人員能否快速到達事先安置好的應急避險場所。在遵循上述兩個原則的前提下，還應確保所選地點盡可能靠近道路，便于人員迅速撤離至應急避險場所。

經(jīng)統(tǒng)計，塘棣源村共計15戶61人將受到尾砂流威脅，根據(jù)《避險轉移圖編制技術要求(試行)》[15]規(guī)定，應急避險場所人均面積不小于3 m2，由此可確定應急避險場所面積應不小于183 m2。故應急避險圓形場所半徑擬定為8 m，總面積約為200 m2，可以滿足受災人員安置要求。

塘棣源村內道路縱橫，撤離最短路線總路程與受災房屋到應急避險場所中心點的直線距離相差不大，故選取受災房屋到應急避險場所中心點的直線作為轉移路線，以此來估算人員撤離時間。人員撤離速度取4 m·s-1[14]。經(jīng)計算，各受災居民的轉移距離在181.26～289.32 m范圍內，轉移時間在45.32～72.33 s內。由3.4節(jié)可知，潰壩尾砂流在漫頂發(fā)生后570 s到達塘棣源村。這表明，按照上述方法制定的避險轉移方案(圖11)，可在確保預警工作有效的前提下，塘棣源村可能受災的居民有充足的時間完成疏散、撤離和避險轉移。

圖11 塘棣源村居民避險轉移方案示意圖

5 結論

(1) 從尾砂流流速、流深、淹沒范圍、泥沙淤積厚度等方面調查了尾砂流對下游居民、村莊公路和農(nóng)田等造成的影響。燕倉尾礦庫漫頂潰壩下泄尾砂共計4.82×105m3，大部分尾砂沿程淤積，其中一部分會尾砂淹沒村莊公路和農(nóng)田。其中，新巖前村淹沒高度為110.5 m，但因新巖前村房屋建筑物高程均高于110.5 m，故該村基本沒有受到潰壩尾砂流影響；相比之下，塘棣源村淹沒高度為94.5 m，因該村部分房屋高程低于94.5 m，故受到了潰壩尾砂流沖擊或淤埋等嚴重影響。

(2)在潰壩尾砂流演進模擬的基礎上制定了塘棣源村居民的避險轉移方案。潰壩尾砂流到達塘棣源村時間為570 s，該村居民轉移時間為45.32～72.33 s。在保證預警工作有效的前提下，該村居民有充足的時間按照制定的方案完成避險轉移。