基于Rhino 和Fluent 耦合的尾礦庫潰壩數值模擬分析

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

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

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

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

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

1 工程概況

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

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

2 潰壩數值模擬

2.1 模擬軟件評述

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

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

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

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

2.2 潰壩模擬區域

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

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

2.3 初始條件假設

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

2.4 模型建立及邊界條件

Fluent 軟件的建模能力有限,而Rhino 軟件中可以通過等高線建立復雜潰壩模擬的地形圖。為此,本研究利用Rhino 建模和Fluent 計算進行耦合,使潰壩數值模擬結果更符合實際情況。將構建的模型導入Fluent 軟件中進行計算,在Fluent 中設置計算模型和邊界條件,具體如下:

(1)讀入網格模型,檢查網格模型質量,若有體積量為負的情況,需要重新定義網格。該模型的網格數量為350 萬。地表形態及主要網格見圖3。

圖3 地表形態及主要網格Fig.3 Surface morphology and main grid

(2)模型設置。選取Transient 非穩定瞬態求解,計算模型選取VOF 多相流模型,相數為2,k-ε湍流模型。

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

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

(5)流場初始化。選取壓力入口作為流場初始化條件,初始流場基本處于靜止狀態,將尾砂部分設置為第二相即尾砂材料,使其體積分數為1。計算在自重作用下的潰壩流動情況。

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

3 模擬結果與分析

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

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

3.1 潰壩水砂淹沒范圍

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

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

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

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

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

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

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

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

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

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

3.3 水砂速度演進情況

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

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

潰壩水砂端頭速度變化曲線如圖8所示。分析圖8 可知:0～300 s 時,速度迅速增加,300 s 時潰壩水砂端頭速度達到2.9 m/s;300～600 s 時,速度持續增加,在600 s 時端頭速度達到最大值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 結 論

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

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

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

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

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

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