胡良才，郭大平，李哲輝，李玉雷

(中核第四研究設(shè)計(jì)工程有限公司，河北 石家莊 050021)

尾礦庫(kù)是礦山的三大控制性工程之一，用于貯存水冶廠選冶礦石后剩余的大量尾礦[1]，受筑壩材料特性及缺乏管理等因素影響，尾礦庫(kù)發(fā)生潰壩事故的概率遠(yuǎn)遠(yuǎn)高于傳統(tǒng)水庫(kù)[2]。尾礦庫(kù)一旦潰壩，產(chǎn)生的泥石流將快速?zèng)_泄向尾礦庫(kù)下游，破壞農(nóng)田以及交通設(shè)施等，給下游的居民造成嚴(yán)重的生命及財(cái)產(chǎn)損失[3]1。

國(guó)內(nèi)外諸多學(xué)者采用經(jīng)驗(yàn)公式法、模型試驗(yàn)、數(shù)值模擬等手段對(duì)尾礦庫(kù)潰壩問(wèn)題進(jìn)行了研究。陳殿強(qiáng)等[4]針對(duì)鳳城市某尾礦庫(kù)，采用經(jīng)驗(yàn)公式計(jì)算得到尾礦庫(kù)潰口寬度、潰壩壩址最大流量、潰壩最大流量沿程演進(jìn)、砂流傳播時(shí)間等，并根據(jù)計(jì)算結(jié)果提出了尾礦庫(kù)預(yù)防及應(yīng)急措施。連國(guó)璽等[5]采用經(jīng)驗(yàn)公式計(jì)算得到某鈾尾礦庫(kù)潰壩影響范圍。尹光志等[6]采用模型試驗(yàn)方法，研究了不同高度下尾礦壩瞬時(shí)全部潰壩情況下泥漿的流態(tài)、演進(jìn)規(guī)律以及動(dòng)力特性，研究了泥漿在尾礦壩下游流動(dòng)過(guò)程中的沖擊力強(qiáng)度、淹沒(méi)范圍及演進(jìn)規(guī)律。張力霆等[7]利用自主研發(fā)的尾礦庫(kù)模型試驗(yàn)平臺(tái)，研究了尾礦庫(kù)潰壩破壞形式。劉磊等[8]采用物理模型試驗(yàn)方法，對(duì)尾礦壩漫頂潰壩洪水及潰口變化過(guò)程進(jìn)行了研究。林江等[9]在潰壩洪水?dāng)?shù)學(xué)模型的基礎(chǔ)上，加入非恒定流輸沙方程，建立能同時(shí)模擬水流和泥沙運(yùn)動(dòng)的二維潰壩水流泥沙數(shù)學(xué)模型，并對(duì)云南某尾礦庫(kù)潰壩后的洪水演進(jìn)及礦砂在下游的淤積過(guò)程進(jìn)行了二維數(shù)值模擬。金佳旭等[10]采用非牛頓流體模型中的Bingham模型，針對(duì)遼寧某尾礦庫(kù)潰壩后的砂流演進(jìn)過(guò)程進(jìn)行了數(shù)值模擬，重點(diǎn)分析了潰壩后泥石流流態(tài)變化、速度矢量變化及最終堆積形態(tài)?；糇犀|[11]采用強(qiáng)度折減法分析某尾礦庫(kù)壩體穩(wěn)定性，計(jì)算得出壩體潰壩范圍，利用計(jì)算流體力學(xué)軟件對(duì)尾礦庫(kù)潰壩進(jìn)行了數(shù)值模擬，得出了潰壩泥石流淹沒(méi)范圍、流動(dòng)速度和淹沒(méi)深度等重要指標(biāo)。Marina Pirulli等[12]采用數(shù)值模擬方法，研究了意大利Stava Valley尾礦庫(kù)潰壩泥石流流動(dòng)特性及影響范圍。陳俊等[13-14]采用數(shù)值模擬方法研究了某銅礦擬建尾礦庫(kù)潰壩情況，通過(guò)設(shè)置監(jiān)測(cè)點(diǎn)，從潰壩泥石流流場(chǎng)、水深和尾砂淤積方面分析潰壩泥石流對(duì)下游村莊的影響。

經(jīng)驗(yàn)公式法不能考慮地形影響，計(jì)算結(jié)果往往與實(shí)際情況存在較大誤差；模型試驗(yàn)法費(fèi)用高、試驗(yàn)周期長(zhǎng)，同時(shí)存在有縮尺效應(yīng)，無(wú)法準(zhǔn)確研究諸多工程中關(guān)心的實(shí)際問(wèn)題。數(shù)值模擬法費(fèi)用低、計(jì)算周期較短，且可以不需考慮縮尺效應(yīng)的影響，直接對(duì)工程原型進(jìn)行模擬研究，可以獲得大量詳實(shí)的資料，還可以進(jìn)行方案比較選擇最優(yōu)設(shè)計(jì)或試驗(yàn)方案，指導(dǎo)模型試驗(yàn)的進(jìn)行，節(jié)省大量試驗(yàn)時(shí)間[3]5。隨著計(jì)算機(jī)硬件及數(shù)值模擬理論的發(fā)展，數(shù)值模擬技術(shù)越來(lái)越多的應(yīng)用于尾礦庫(kù)的潰壩研究。筆者針對(duì)某鈾尾礦庫(kù)，建立包括尾礦庫(kù)壩體及下游三維精細(xì)地形的數(shù)值模型，進(jìn)行尾礦庫(kù)潰壩三維數(shù)值模擬研究。

1 研究方法

1.1 數(shù)值模型

RNGk-ε紊流數(shù)值模型可更好地模擬流體流線彎曲程度較大的流動(dòng)[15]，因此采用RNGk-ε紊流數(shù)值模型來(lái)封閉基本方程組，對(duì)鈾尾礦庫(kù)潰壩進(jìn)行三維數(shù)值模擬研究，控制方程[16]如下。

連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

紊動(dòng)能k方程:

(3)

耗散率ε方程：

(4)

采用有限體積法來(lái)離散控制方程，VOF(Volume of Fluid)方法[18]追蹤水流的自由表面。VOF方法定義一個(gè)體積函數(shù)F，若單元體內(nèi)充滿流體，F(xiàn)=1；若單元體內(nèi)無(wú)流體，F(xiàn)=0；單元體內(nèi)同時(shí)含有流體以及空氣，0

相間界面的追蹤方程可表示為

(5)

1.2 模型驗(yàn)證

采用Martin等[19]的水槽試驗(yàn)驗(yàn)證潰壩數(shù)值模型。根據(jù)數(shù)值模擬軟件特點(diǎn)，建立與試驗(yàn)?zāi)Ｐ统叽缦嗤娜S數(shù)值模型，水槽試驗(yàn)及數(shù)值模型示意圖如圖1所示。

模型計(jì)算區(qū)域長(zhǎng)度l=12.7 cm，寬度b=5.72 cm，高度h=6.72 cm。采用結(jié)構(gòu)化網(wǎng)格，順?biāo)鞣较驗(yàn)閤方向，重力方向?yàn)閦方向，網(wǎng)格的單元尺寸Δx=0.025 cm，Δy=Δz=0.1 cm，設(shè)置計(jì)算軟件自動(dòng)調(diào)整時(shí)間步長(zhǎng)大小，確定初始的時(shí)間步長(zhǎng)Δt=0.001 s。采用VOF方法追蹤流體自由表面，RNGk-ε紊流數(shù)值模型封閉控制方程。數(shù)值模型上游、兩側(cè)以及底部的邊界條件設(shè)定為法向速度為零、無(wú)滑移的固體邊壁；模型出口為自由出流，設(shè)置為壓力出口邊界條件；模型頂部給定為壓力入口邊界條件。確定初始水體長(zhǎng)l0=2.86 cm，寬b0=5.72 cm，高h(yuǎn)0=5.72 cm。水槽試驗(yàn)與數(shù)值模擬水流波前同時(shí)間關(guān)系曲線如圖2所示。

由圖2可見(jiàn)，對(duì)于水流與時(shí)間關(guān)系的模擬，數(shù)值模擬與模型試驗(yàn)結(jié)果規(guī)律一致，誤差較小，說(shuō)明采取的數(shù)值模型真實(shí)可靠。

2 尾礦庫(kù)潰壩數(shù)值模擬

2.1 地形建模

采用三維建模軟件CATIA建立尾礦庫(kù)庫(kù)區(qū)及下游河道三維地形模型：1)使用Autodesk Civil 3D軟件，根據(jù)實(shí)測(cè)地形圖建立地形曲面，并提取庫(kù)區(qū)及下游河道地形坐標(biāo)及高程，導(dǎo)出點(diǎn)云數(shù)據(jù)；2)進(jìn)入CATIA的Digitized Shape Editor模塊，導(dǎo)入點(diǎn)云數(shù)據(jù)，過(guò)濾、修補(bǔ)點(diǎn)云數(shù)據(jù)(圖3a)；3)根據(jù)導(dǎo)入的點(diǎn)云數(shù)據(jù)生成mesh網(wǎng)格面，分析、檢查、修補(bǔ)網(wǎng)格面(圖3b)；4)進(jìn)入Quick Surface Reconstruction模塊，生成實(shí)體地表面(圖3c)；5)進(jìn)入Part模塊，將地形邊界投影到基準(zhǔn)面上作為草圖邊界，拉伸草圖至曲面，形成三維地形實(shí)體(圖3d)。

2.2 尾礦庫(kù)三維數(shù)值模型

某鈾尾礦庫(kù)主要由初期壩、尾礦堆積壩、副壩、排洪設(shè)施等組成。初期壩高25 m，壩頂標(biāo)高125 m，尾礦堆積壩高27 m，總壩高52 m。初期壩上游壩坡坡度為1∶2.5，下游壩坡坡度為1∶2～1∶2.5；初期壩內(nèi)庫(kù)容堆滿后，后期采用尾礦砂堆積筑壩，尾礦堆積壩坡度為1∶3.5[20]?；趯?shí)測(cè)的1∶1 000比例地形圖，建立鈾尾礦庫(kù)局部潰壩三維模型，模型長(zhǎng)4 000 m，寬3 000 m。模型計(jì)算區(qū)域包括尾礦庫(kù)庫(kù)區(qū)、壩體及尾礦壩下游5 km范圍內(nèi)河道。為了避免因計(jì)算過(guò)程中發(fā)生封頂現(xiàn)象而導(dǎo)致計(jì)算失敗，設(shè)定庫(kù)區(qū)最大模擬高程為160 m，高于尾礦庫(kù)現(xiàn)狀最大壩高[17]140。

假定尾礦庫(kù)潰壩產(chǎn)生的泥石流為介于流體和散粒體之間的一種特殊的運(yùn)動(dòng)形式，可采用流體流動(dòng)的能量方程和連續(xù)方程近似描述，數(shù)值模型采用如下基本假定：1)尾礦庫(kù)庫(kù)區(qū)基巖及周?chē)襟w為不透水邊界；2)潰壩泥石流演進(jìn)過(guò)程中，不考慮單個(gè)顆粒的體積變形；3)潰壩泥石流為各向同性連續(xù)介質(zhì)流體，其流動(dòng)符合賓漢流動(dòng)模型的規(guī)律；4)不考慮潰壩泥石流對(duì)下游河道的沖刷。

尾礦庫(kù)潰壩從形成決口到基本形成穩(wěn)定的潰決斷面，時(shí)間過(guò)程非常短暫，為安全考慮，尾礦庫(kù)潰壩按瞬時(shí)潰壩處理。根據(jù)某鈾尾礦庫(kù)壩體穩(wěn)定性分析計(jì)算得到的壩體滑動(dòng)面位置[21]，確定尾礦壩滑移破壞區(qū)，即為局部潰壩范圍。對(duì)校核工況下，尾礦庫(kù)瞬時(shí)局部潰壩情況進(jìn)行三維數(shù)值模擬。根據(jù)黃河水利委員會(huì)水利科學(xué)研究院實(shí)際資料分析求得潰壩決口平均寬度b，計(jì)算公式為

b=k(W1/2B1/2H0)1/2，

(6)

式中：b—潰壩決口平均寬度，m；W—潰壩時(shí)庫(kù)容，萬(wàn)m3；B—壩頂長(zhǎng)度，m；H0—壩前水深(水頭)，m；k—與壩體土質(zhì)有關(guān)的系數(shù)，對(duì)黏土k值約為0.65，壤土k值約為1.3。計(jì)算得到尾礦庫(kù)壩體潰口寬度為102.9 m。

采用結(jié)構(gòu)化網(wǎng)格對(duì)數(shù)值模型進(jìn)行網(wǎng)格劃分，自尾礦庫(kù)壩體向下游方向?yàn)閤方向，重力方向?yàn)閦方向。經(jīng)過(guò)多次試算后確定計(jì)算網(wǎng)格尺寸為Δx=Δy=6.0 m，Δz=2.0 m，網(wǎng)格總數(shù)約1 800萬(wàn)。設(shè)定初始時(shí)間步長(zhǎng)Δt=0.001 s，設(shè)置計(jì)算軟件可以自動(dòng)調(diào)整時(shí)間步長(zhǎng)大小。采用VOF方法追蹤流體自由表面，RNGk-ε紊流數(shù)值模型來(lái)封閉紊流基本方程組。尾礦庫(kù)庫(kù)區(qū)模型頂面設(shè)置為壓力入口邊界，壓力為101.325 kPa；庫(kù)底及下游河道設(shè)定為法向速度零、無(wú)滑移的固體邊壁，采用標(biāo)準(zhǔn)壁面函數(shù)法處理近壁的黏性底層；下游河道出口方向?yàn)樽杂沙隽?，設(shè)置為壓力出口邊界，壓力為101.325 kPa。尾礦庫(kù)下游為水田、灌木叢等，下墊面情況較為復(fù)雜，根據(jù)國(guó)內(nèi)外相關(guān)研究成果，下游河道曼寧系數(shù)n取0.03。計(jì)算初始條件為庫(kù)區(qū)內(nèi)壩頂標(biāo)高以下水的體積分?jǐn)?shù)為1，即庫(kù)區(qū)被水體充滿。在尾礦庫(kù)壩址處至壩下游5 km處，每1 km設(shè)置1個(gè)監(jiān)測(cè)點(diǎn)，共計(jì)6個(gè)監(jiān)測(cè)點(diǎn)，以監(jiān)測(cè)潰壩泥石流流速、傳播時(shí)間、泥石流深等水力要素變化情況。

3 結(jié)果分析

3.1 流速變化分析

圖4為尾礦庫(kù)壩址處流速變化過(guò)程線。由圖可見(jiàn)，尾礦庫(kù)局部潰壩情況下，壩址處潰壩泥石流流速隨時(shí)間增加迅速增大至峰值，之后逐漸下降，最終趨于零。潰壩5 s后壩址處流速達(dá)到峰值，最大流速約為14.12 m/s。

3.2 傳播時(shí)間分析

尾礦庫(kù)潰壩后，潰壩泥石流快速?zèng)_泄向下游，從洪水到達(dá)至最大洪峰到達(dá)時(shí)，兩者時(shí)間間隔僅為數(shù)秒至數(shù)十秒，且越靠近尾礦庫(kù)，間隔時(shí)間越短，從而留給下游民眾安全撤離的時(shí)間越短。因此，鈾尾礦庫(kù)運(yùn)行管理過(guò)程中，應(yīng)加強(qiáng)監(jiān)測(cè)，提前預(yù)警，為居民安全撤離爭(zhēng)取寶貴的時(shí)間。局部潰壩情況下，尾礦庫(kù)壩址及下游各特征斷面洪水起漲時(shí)間、最大洪峰傳播時(shí)間、洪水消落時(shí)間見(jiàn)圖5。

3.3 泥石流深結(jié)果分析

尾礦庫(kù)壩址及下游各特征斷面監(jiān)測(cè)點(diǎn)泥石流深變化見(jiàn)圖6。潰壩事故發(fā)生后，壩址處及下游各特征斷面潰壩泥石流深隨時(shí)間增加迅速增大，之后逐漸下降。在潰壩40 s后，壩址處泥石流深達(dá)到最大值，約為12.31 m。

3.4 影響范圍分析

根據(jù)數(shù)值模擬計(jì)算得到的各特征斷面最大泥石流深，結(jié)合尾礦庫(kù)實(shí)測(cè)地形圖，繪制局部潰壩泥石流影響范圍見(jiàn)圖7。在尾礦庫(kù)下游的潭元等幾個(gè)村莊設(shè)置有監(jiān)測(cè)點(diǎn)，各監(jiān)測(cè)點(diǎn)均未監(jiān)測(cè)到潰壩泥石流，說(shuō)明局部潰壩情況下，潰壩泥石流未直接沖擊影響到這幾個(gè)村落。

4 結(jié)論

基于計(jì)算流體力學(xué)軟件Flow-3D建立潰壩三維數(shù)值模型，采用水槽試驗(yàn)驗(yàn)證數(shù)值模型，結(jié)果表明建立的數(shù)值模型能很好的模擬水流運(yùn)動(dòng)情況。采用三維設(shè)計(jì)軟件CATIA建立某鈾尾礦庫(kù)壩體三維模型及下游河道三維精細(xì)地形，基于VOF方法、RNGk-ε紊流模型，建立某鈾尾礦庫(kù)潰壩三維數(shù)值模型，對(duì)尾礦庫(kù)局部潰壩進(jìn)行的三維數(shù)值模擬表明：尾礦庫(kù)潰壩后，潰壩泥石流流速、深度等迅速增大到峰值，之后逐漸平穩(wěn)下降；局部潰壩情況下，尾礦庫(kù)潰壩泥石流不會(huì)直接沖擊影響下游各村莊；RNGk-ε紊流模型可用于對(duì)尾礦庫(kù)潰壩泥石流水力特性進(jìn)行三維數(shù)值模擬，研究成果可為工程設(shè)計(jì)、運(yùn)行管理提供很好的借鑒。