孫鴻昌 郝 喆 侯永莉 尹亮亮 張 穎

(1.中煤科工集團(tuán)沈陽研究院有限公司,遼寧 撫順 113122;2.煤礦安全技術(shù)國家重點實驗室,遼寧 撫順 113122;3.遼寧大學(xué)環(huán)境學(xué)院,遼寧 沈陽 110036;4.遼寧有色勘察研究院有限責(zé)任公司,遼寧 沈陽 110033)

尾礦庫是高勢能的泥石流源頭[1],我國尾礦庫具有體積小、數(shù)量眾多的特點,加之我國人口多、分布廣[2],尾礦庫的建造場地往往很難避開人類居住區(qū),導(dǎo)致多數(shù)尾礦庫下游存在人類活動[3],形成了“頭頂庫”。根據(jù)相關(guān)統(tǒng)計,截至2018年底,我國仍然存在“頭頂庫”1 425 座,一旦發(fā)生潰壩,將嚴(yán)重威脅下游居民和重要工農(nóng)業(yè)、交通設(shè)施及環(huán)境安全[4]。近年來“頭頂庫”潰壩事故屢見不鮮,例如2008年9月,山西省新塔礦業(yè)有限公司發(fā)生尾礦庫潰壩事故,波及下游500 m 左右的礦區(qū)辦公樓、集貿(mào)市場和部分民宅,房屋摧毀若干,造成276 人死亡;2015年11月,湖南郴州云錫礦業(yè)尾礦庫因排泄設(shè)施損壞,導(dǎo)致4 人失聯(lián),造成7.9 億的經(jīng)濟損失以及嚴(yán)重的環(huán)境污染[5]。為此,一些學(xué)者開展了尾礦庫潰壩模擬研究,李全明等[6]以潰壩泥漿的動力學(xué)過程作為研究重點,得到了尾礦庫潰壩水砂流動的計算方法;郝喆等[7]通過ANSYS 和FLCA3D相結(jié)合,對尾礦庫加高增容過程中的壩體滲流進(jìn)行了穩(wěn)定性分析;陳殿強等[8]采用大壩潰決的經(jīng)驗公式,分析了某鐵礦尾礦庫潰壩時的尾砂淹沒范圍及深度等;王興華等[9]通過Fluent 數(shù)值模擬研究了尾礦庫潰壩洪水不同時刻的流動特性;HANSON 等[10]通過多次大規(guī)模逐漸潰壩模擬試驗,得出漫頂引起的堤壩潰決的形成過程和時間可以極大地影響水庫放水的速率;WANG 等[11]對潰壩后下泄洪水的演變進(jìn)行了數(shù)值模擬。

根據(jù)現(xiàn)有成果分析可知,當(dāng)前的尾礦庫潰壩模擬大多采用基于數(shù)學(xué)模型的理論分析方式,或者進(jìn)行實際地形簡化后的二維數(shù)值模擬思路。潰壩的三維數(shù)值模擬在水庫大壩中應(yīng)用較多,在尾礦庫中應(yīng)用較少,而尾礦庫潰壩后形成水砂兩相流,其流動狀態(tài)與水庫的洪水單相流有很大區(qū)別。前人多采用單一數(shù)值模擬軟件進(jìn)行穩(wěn)定性計算或潰壩流動狀態(tài)分析,在精細(xì)地形描述和高精度建模方面存在不足。潰壩模擬中多將流動狀態(tài)等效為帶側(cè)板的明渠流,缺乏三維真實地形下的流動過程模擬。因此,采用專業(yè)的3D造型軟件和具有強大計算能力的數(shù)值模擬軟件進(jìn)行耦合分析,實現(xiàn)高精度復(fù)雜地形下的尾礦庫水砂兩相流三維潰壩模擬,具有重要意義。

作為各級應(yīng)急管理部門正在推行的非煤礦山安全管理新方向,尾礦庫潰壩數(shù)值模擬技術(shù)標(biāo)準(zhǔn)尚未建立,分析方法、建模精度、計算內(nèi)容和成果要求等也未明確。本研究以遼寧省某“頭頂庫”為例,采用VOF多相流模型與k-ε湍流模型理論,利用Rhino 和Fluent 軟件相耦合,考慮下游存在溝谷、農(nóng)田、村莊、河流和公路等狀態(tài)的復(fù)雜地形條件,開展尾礦庫加高后潰壩的三維數(shù)值模擬,研究漫頂潰壩后的尾砂流動規(guī)律,獲得尾砂在不同時刻對下游的淹沒范圍、堆積情況和流動狀態(tài),對尾礦庫防災(zāi)減災(zāi)和環(huán)境保護(hù)具有一定的參考意義。所提出的耦合建模分析思路,可為類似“頭頂庫”潰壩分析和標(biāo)準(zhǔn)建立提供一定的借鑒。

1 工程概況

某尾礦庫位于遼寧省鳳城市,尾礦庫平面如圖1所示。該尾礦庫初期壩設(shè)計為堆石壩,壩底標(biāo)高222.70 m,壩頂標(biāo)高232.00 m,相對高度9.30 m。主壩與副壩采用上游法筑壩,內(nèi)坡比為1 ∶2,外坡比為1 ∶1.5。該尾礦庫為山谷型,設(shè)計加高增容后的堆積壩頂標(biāo)高為285.00 m,堆積后總壩高約48.00 m。西側(cè)堆積壩體為主壩平均外坡比為1 ∶4.5,南側(cè)副壩外坡比為1 ∶2.4,東側(cè)壩體為擋水堆石壩體,外坡比約1 ∶2.0,庫容約1 300 萬m3。該尾礦庫下游為農(nóng)田、公路以及農(nóng)村,副壩南側(cè)為廠區(qū),綜上,根據(jù)《尾礦堆積壩巖土工程技術(shù)規(guī)范》(GB 50547—2010)[12]的規(guī)定,按現(xiàn)狀壩高、庫容等綜合考慮,該尾礦庫等級為Ⅲ等庫,尾礦壩級別為3 級。

圖1 遼寧某尾礦庫平面(單位:m)Fig.1 Plan of a tailings pond in Liaoning Province

2 潰壩數(shù)值模擬

2.1 模擬軟件評述

本研究數(shù)值模擬軟件采用Rhino 建模工具[13]和Fluent 流體分析軟件。

(1)Rhino 軟件具備較傳統(tǒng)建模思路更為出色的NURBS 建模方式[14],可以更好地控制模型表面的曲線度,創(chuàng)建出更逼真、生動的復(fù)雜模型和效果圖,實現(xiàn)高精度的復(fù)雜地形真實模擬。

(2)Fluent 軟件具有較高的計算速度、穩(wěn)定性和精度,強大的后處理功能,可以模擬各種類型流體的真實流動狀態(tài),但軟件自帶的建模功能有很大的局限性[15]。

(3)Rhino 軟件和Fluent 軟件相耦合,可有效解決目前數(shù)值模擬建模能力差、地形精度低、模型太過理性化、缺乏真實性的問題,能夠有效體現(xiàn)出下游復(fù)雜地形對潰壩水砂流動的影響。

2.2 潰壩模擬區(qū)域

在Rhino 軟件中建立尾礦庫加高增容擴建工程終期條件下的潰壩三維計算模型,模型內(nèi)包括初期壩、堆積壩及所在區(qū)域的實際地形。研究區(qū)域長度為3.1 km,最大寬度1.8 km,研究范圍約4.65 km2,研究區(qū)域內(nèi)主要有青城子村、新式街、松樹溝、雙下等4個村莊以及潰壩下游農(nóng)田。尾礦庫建模范圍如圖2所示。

圖2 尾礦庫建模范圍Fig.2 Modeling scope of tailings pond

2.3 初始條件假設(shè)

考慮到尾礦庫生產(chǎn)時間不斷加長,排洪設(shè)施有可能發(fā)生尾砂下陷、尾砂沉積造成堵塞,或者排水斜槽發(fā)生坍塌,預(yù)制管接口滲漏等情況[16]。在這種極端情況下排洪能力下降或消失,庫內(nèi)水位升高,浸潤線埋深抬升過高,從而導(dǎo)致尾礦庫發(fā)生潰壩[17]。為簡化計算并進(jìn)一步提高模擬結(jié)果與實際情況的契合度,本研究做出兩點假設(shè):① 假設(shè)500 a 一遇洪水過程超過10 h,庫內(nèi)水位持續(xù)上漲達(dá)到壩頂,洪水位上漲的時間過程在數(shù)值模擬中忽略;② 假設(shè)庫內(nèi)排洪系統(tǒng)失效,計算過程中的庫內(nèi)水砂主要從潰口下泄,且潰壩模擬求解中無其他外來水量疊加計算。

2.4 模型建立及邊界條件

Fluent 軟件的建模能力有限,而Rhino 軟件中可以通過等高線建立復(fù)雜潰壩模擬的地形圖。為此,本研究利用Rhino 建模和Fluent 計算進(jìn)行耦合,使?jié)螖?shù)值模擬結(jié)果更符合實際情況。將構(gòu)建的模型導(dǎo)入Fluent 軟件中進(jìn)行計算,在Fluent 中設(shè)置計算模型和邊界條件,具體如下:

(1)讀入網(wǎng)格模型,檢查網(wǎng)格模型質(zhì)量,若有體積量為負(fù)的情況,需要重新定義網(wǎng)格。該模型的網(wǎng)格數(shù)量為350 萬。地表形態(tài)及主要網(wǎng)格見圖3。

圖3 地表形態(tài)及主要網(wǎng)格Fig.3 Surface morphology and main grid

(2)模型設(shè)置。選取Transient 非穩(wěn)定瞬態(tài)求解,計算模型選取VOF 多相流模型,相數(shù)為2,k-ε湍流模型。

(3)設(shè)置材料參數(shù)。材料為潰壩尾砂和空氣兩種,由于潰壩常伴隨著暴雨等條件的發(fā)生,潰壩時的尾砂常處于飽和狀態(tài),所以根據(jù)實際調(diào)查和試驗結(jié)果,選取尾砂的平均密度為1 750 kg/m3,黏度系數(shù)為15 Pa·s,空氣選取標(biāo)準(zhǔn)大氣壓下的空氣參數(shù)。將空氣設(shè)置為主相,潰壩尾砂設(shè)置為次相。

(4)邊界條件。將模型上部空氣界面設(shè)置為壓力入口,四周界面為壓力出口,山體地形和壩體為墻體。

(5)流場初始化。選取壓力入口作為流場初始化條件,初始流場基本處于靜止?fàn)顟B(tài),將尾砂部分設(shè)置為第二相即尾砂材料,使其體積分?jǐn)?shù)為1。計算在自重作用下的潰壩流動情況。

(6)迭代求解。設(shè)置時間步為0.01 s,每次迭代循環(huán)步數(shù)20 s,迭代總步數(shù)18 萬步。

3 模擬結(jié)果與分析

本研究采用VOF 計算模型與k-ε湍流模型理論相結(jié)合的計算方法[18]。VOF 法在計算域內(nèi),引入了相體積分?jǐn)?shù)F,即在網(wǎng)格內(nèi)各相流體的體積與流體空間總體積的比值[19];選取適用于剪切流動強、邊界層分離等復(fù)雜流動情況[20]的k-ε湍流模型。

尾礦庫發(fā)生潰壩后,庫區(qū)內(nèi)水砂向下游傾瀉而下,水砂的淹沒范圍能夠直觀地反映出不同時刻水砂對下游的淹沒情況;水砂對地面的壓力能夠計算出水砂的淹沒深度;下泄水砂端頭速度大小能夠反映出水砂對下游構(gòu)筑物沖擊力的大小。因此根據(jù)模擬結(jié)果,對潰壩后水砂的淹沒范圍、壓力、速度進(jìn)行進(jìn)一步分析。

3.1 潰壩水砂淹沒范圍

水砂流的演進(jìn)情況可表示尾礦庫潰壩的全部過程。根據(jù)分析結(jié)果,提取了尾礦庫發(fā)生潰壩過程以及潰決水砂流中不同時刻流動演進(jìn)情況的體積分?jǐn)?shù)(Volume fraction)云圖,如圖4所示。

圖4 不同時刻潰壩水砂流演進(jìn)云圖Fig.4 Nephogram of water and sand flow evolution of dam break at different times

由圖4 可知:當(dāng)尾礦庫排洪系統(tǒng)失效同時遭遇500 a 一遇洪水發(fā)生潰壩,壩體破壞狀態(tài)表現(xiàn)為庫內(nèi)水砂漫過主壩壩頂,從壩頂傾瀉而下,隨后主壩頂破碎,受洪水拖曳力的作用,壩坡尾砂受沖刷下泄,庫內(nèi)洪水誘發(fā)尾砂流態(tài)下泄,演化運動過程位移變化明顯。t=0 時,尾礦堆積壩體處于相對靜止?fàn)顟B(tài),隨著時間不斷增長,堆積壩壩頂潰口持續(xù)發(fā)展,在0～100 s 時,潰口附近尾砂由于水流拖曳力作用沖刷形成初期水砂混合流,后方尾砂隨之開始滑動,臨界流動面向后方擴展;在100～300 s 時,潰口后方尾礦堆積體逐漸滑出,在后方滑出水砂的流動推動下,壩坡尾砂出現(xiàn)沖刷溝壑;400 s 左右前端砂流淹沒初期壩,此后階段為水砂流在初期壩下游的演進(jìn)階段;在600～900 s 時,水砂流在初期壩外坡腳至所在溝谷口范圍流動運移;在900～1 500 s 時,潰壩的勢能轉(zhuǎn)化為動能過程中伴隨尾砂料之間的碰撞摩擦導(dǎo)致能量消耗,前端水砂流到達(dá)低勢能區(qū),在下游村莊和公路所在平地減速緩慢演進(jìn)淤積;t=1 800 s 之后尾礦壩體滑動逐漸達(dá)到新的平衡狀態(tài),壩體垮塌運動趨于停止,此時水砂到達(dá)下游的1 035 m 位置。此次潰壩水砂量約297萬m3,水砂覆蓋下游面積約0.98 km2。

3.2 潰壩水砂壓力(淹沒高度)情況

尾礦庫下游地形條件的復(fù)雜情況對潰壩的影響程度較大[21]。若下游高差大,相應(yīng)的潰出距離及影響范圍更大;若尾礦庫下游區(qū)域為平原開闊區(qū)域,在縱坡度較緩的情況下,相對其他下游地形條件,影響距離較小。因此,根據(jù)分析結(jié)果,提取了1 800 s 尾礦庫潰決后水砂對地面壓力云圖并選取距初期壩外壩腳250、500、750、1 000 m 4 處作為水砂壓力監(jiān)測點,如圖5所示。

圖5 水砂對地面的壓力分布及監(jiān)測點位置Fig.5 Pressure distribution of water and sand on the ground and the position of monitoring points

由圖5 可知:大部分尾砂在距離初期壩外坡腳500～600 m 區(qū)域內(nèi)淤積。根據(jù)壓強公式,折合尾砂最大淹沒深度為7.2 m,位于下游500 m 處的溝谷內(nèi),潰壩下游其他位置水砂的淹沒高度為2～3 m。

本次潰壩模擬泥石流漸進(jìn)式泄向下游地區(qū)。根據(jù)監(jiān)測點處的壓力變化規(guī)律,繪制了水砂壓力隨時間的變曲線,如圖6所示。

圖6 與壩腳不同距離處的壓力—時間曲線Fig.6 Pressure-time curves at different distances from the foot of the dam

由圖6 可知:各位置處水砂壓力在不同時刻出現(xiàn)峰值,反映了潰壩水砂在下游逐步淹沒的進(jìn)程。水砂流到監(jiān)測點位置后,壓力以約1.1 kPa/s(合0.064 m/s)的高梯度增加到最大值,然后迅速減低,隨著時間增加,監(jiān)測點壓力衰減速度越來越緩慢,最后趨于穩(wěn)定不再變化。下游750 m 處的農(nóng)田地勢較低,壓力略高于250 m 和1 000 m 處,下游500 m 處地勢低洼且有大量尾砂在此堆積,所以水砂壓力最大。隨后由于流距增長以及地面阻力的共同作用,水砂壓力開始減弱,最后趨于穩(wěn)定。

3.3 水砂速度演進(jìn)情況

潰壩后不同時刻的水砂速度云圖如圖7所示。由圖7 可知:t=300 s 時,潰壩水砂到達(dá)初期壩底部,由于潰口與壩底存在48 m 的高差,潰口水砂在重力勢能和堆積壩體摩擦的共同作用下從0 迅速增大;隨后由于潰口增大,水砂流量加大;t=600 s 時,水砂流經(jīng)下游500 m 處,水砂流速進(jìn)一步增加,快速向下游流動;隨著時間不斷增長,堆積壩壩頂潰口持續(xù)發(fā)展,受下游地勢平緩和地表面粗糙程度的影響,水砂的整體流速降低,尾砂在500 m 處堆積;此后由于尾礦庫水砂重力勢能降低;t=1 800 s 時,潰口與庫區(qū)內(nèi)水位逐步持平,水砂的流速非常小,不會向下游產(chǎn)生較大區(qū)域的流動。

圖7 不同時刻潰壩水砂速度云圖Fig.7 Water and sand velocity nephogram of dam break at different times

潰壩水砂端頭速度變化曲線如圖8所示。分析圖8 可知:0～300 s 時,速度迅速增加,300 s 時潰壩水砂端頭速度達(dá)到2.9 m/s;300～600 s 時,速度持續(xù)增加,在600 s 時端頭速度達(dá)到最大值3.7 m/s,此時水砂的沖擊力最大;之后水砂流速逐漸降低,到1 800 s 時流速趨于0 停止不動。

圖8 潰壩水砂端頭速度—時間曲線Fig.8 Velocity-time cures of water and sand at the end of dam break

4 結(jié) 論

針對精細(xì)真實三維地形下尾礦庫潰壩模擬分析存在的不足,通過Rhino 建立高精度地形圖,與Flunet 模擬軟件相耦合,分析了潰壩水砂受真實地形影響的流動狀態(tài),得出以下結(jié)論:

(1)采用VOF 多相流計算模型與k-ε湍流模型理論,利用Rhino 建模和Fluent 計算相耦合,可以有效實現(xiàn)高精度地形下的水砂兩相流潰壩三維數(shù)值模擬,獲得尾砂在不同時刻對下游的淹沒范圍、堆積特征和流動狀態(tài)。

(2)“頭頂庫”發(fā)生潰壩后,在下游地勢較平緩和谷口較寬闊的情況下,發(fā)生潰壩后水砂的淹沒范圍廣泛,其下游水砂淹沒形狀近似為三角形:距尾礦壩較近位置淹沒范圍寬,水砂前端受溝谷地形影響淹沒范圍窄。尾礦庫潰壩后的潰決持續(xù)時間長,水砂的淹沒范圍達(dá)到下游的1 km 左右。

(3)潰壩水砂漸進(jìn)式泄向下游地區(qū),潰壩水砂在下游各處壓力都具有先增大后減小,最后趨于穩(wěn)定的特點。大部分尾砂會淤積在潰壩下游溝谷內(nèi),導(dǎo)致壓強高于其他地區(qū)。

(4)發(fā)生潰壩后,潰口水砂在重力勢能和堆積壩體摩擦的共同作用下從0 迅速增大,在某一時刻,流速達(dá)到峰值,此時水砂沖擊力最大,對下游破壞力最強;隨后,潰口水位不斷降低,受下游地勢平緩和地表面粗糙程度的影響,水砂流速緩慢降低不會向下游產(chǎn)生較大區(qū)域的流動。

(5)在汛期來臨前,建議拓寬加深下游河道,增強溝谷的輸水性和堆積能力,降低水砂流通性,減緩潰壩水砂流速和淹沒范圍,防止對下游產(chǎn)生較大的沖擊力和危害;在下游建立攔擋導(dǎo)流設(shè)施,降低水砂的淹沒深度,延長水砂到達(dá)下游村莊的時間,保障居民生命財產(chǎn)安全。根據(jù)本研究尾礦庫潰壩水砂淹沒范圍的模擬結(jié)果,建議在下游淹沒范圍的1.5 倍區(qū)域內(nèi)禁止建設(shè)廠礦和村落。