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

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

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

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

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

1 研究方法

1.1 數(shù)值模型

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

連續(xù)性方程：

(1)

動量方程：

(2)

紊動能k方程:

(3)

耗散率ε方程：

(4)

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

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

(5)

1.2 模型驗證

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

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

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

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

2.1 地形建模

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

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

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

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

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

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

(6)

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

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

3 結(jié)果分析

3.1 流速變化分析

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

3.2 傳播時間分析

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

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

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

3.4 影響范圍分析

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

4 結(jié)論

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