精細(xì)地形下鈾尾礦庫潰壩三維數(shù)值模擬研究

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

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

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

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

經(jīng)驗(yàn)公式法不能考慮地形影響，計(jì)算結(jié)果往往與實(shí)際情況存在較大誤差；模型試驗(yàn)法費(fèi)用高、試驗(yàn)周期長，同時(shí)存在有縮尺效應(yīng)，無法準(zhǔn)確研究諸多工程中關(guān)心的實(shí)際問題。數(shù)值模擬法費(fèi)用低、計(jì)算周期較短，且可以不需考慮縮尺效應(yīng)的影響，直接對工程原型進(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ù)越來越多的應(yīng)用于尾礦庫的潰壩研究。筆者針對某鈾尾礦庫，建立包括尾礦庫壩體及下游三維精細(xì)地形的數(shù)值模型，進(jìn)行尾礦庫潰壩三維數(shù)值模擬研究。

1 研究方法

1.1 數(shù)值模型

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

連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

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

(3)

耗散率ε方程：

(4)

采用有限體積法來離散控制方程，VOF(Volume of Fluid)方法[18]追蹤水流的自由表面。VOF方法定義一個(gè)體積函數(shù)F，若單元體內(nèi)充滿流體，F(xiàn)=1；若單元體內(nèi)無流體，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ū)域長度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í)間步長大小，確定初始的時(shí)間步長Δt=0.001 s。采用VOF方法追蹤流體自由表面，RNGk-ε紊流數(shù)值模型封閉控制方程。數(shù)值模型上游、兩側(cè)以及底部的邊界條件設(shè)定為法向速度為零、無滑移的固體邊壁；模型出口為自由出流，設(shè)置為壓力出口邊界條件；模型頂部給定為壓力入口邊界條件。確定初始水體長l0=2.86 cm，寬b0=5.72 cm，高h(yuǎn)0=5.72 cm。水槽試驗(yàn)與數(shù)值模擬水流波前同時(shí)間關(guān)系曲線如圖2所示。

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

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

2.1 地形建模

采用三維建模軟件CATIA建立尾礦庫庫區(qū)及下游河道三維地形模型：1)使用Autodesk Civil 3D軟件，根據(jù)實(shí)測地形圖建立地形曲面，并提取庫區(qū)及下游河道地形坐標(biāo)及高程，導(dǎo)出點(diǎn)云數(shù)據(jù)；2)進(jìn)入CATIA的Digitized Shape Editor模塊，導(dǎo)入點(diǎn)云數(shù)據(jù)，過濾、修補(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 尾礦庫三維數(shù)值模型

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

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

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

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

(6)

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

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

3 結(jié)果分析

3.1 流速變化分析

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

3.2 傳播時(shí)間分析

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

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

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

3.4 影響范圍分析

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

4 結(jié)論

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