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

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

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

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

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

1 工程概況

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

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

2 潰壩數(shù)值模擬

2.1 模擬軟件評(píng)述

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

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

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

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

2.2 潰壩模擬區(qū)域

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

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

2.3 初始條件假設(shè)

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

2.4 模型建立及邊界條件

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

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

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

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

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

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

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

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

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

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

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

3.1 潰壩水砂淹沒(méi)范圍

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

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

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

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

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

圖5 水砂對(duì)地面的壓力分布及監(jiān)測(cè)點(diǎ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ù)壓強(qiáng)公式,折合尾砂最大淹沒(méi)深度為7.2 m,位于下游500 m 處的溝谷內(nèi),潰壩下游其他位置水砂的淹沒(méi)高度為2～3 m。

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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