孫鴻昌，郝 喆，楊青潮

（1.遼寧大學環(huán)境學院,遼寧 沈陽 110036；2.長安大學公路學院,陜西 西安 710064）

我國很多尾礦庫已進入中老年期，接近或達到設計庫容，尾礦壩也接近設計高度[1]；且多為上游式尾礦庫，規(guī)模小、數(shù)量多，下游多分布人類生產(chǎn)生活區(qū)而形成“頭頂庫”[2]。受庫水位影響，上游式尾礦庫的浸潤線不斷變化，導致壩體飽和區(qū)范圍改變，使得壩體強度受到嚴重影響，一旦發(fā)生潰壩，下游人民生命財產(chǎn)安全將受到嚴重損害[3-4]。因此，開展不同尾礦庫水位下的尾礦壩滲流破壞及潰壩淹沒范圍的研究具有重要現(xiàn)實意義。

國內(nèi)外學者在邊坡穩(wěn)定性及潰壩分析方面取得了一定的研究進展。劉克輝等[5]利用物理模型試驗裝置開展溝槽在不同彎曲度下的潰壩模擬研究，得出溝槽彎曲度與潰壩泥漿淹沒深度呈非線性增加；郝喆等[6]結(jié)合Ansys 和FLAC3D軟件，開展尾礦庫加高增容過程中壩體滲流的穩(wěn)定性分析及研究；陳勝亮等[7]利用多場耦合分析法，探討了邊坡非飽和滲流及其形變量等變化規(guī)律，得出邊坡應變軟化被非飽和滲流的時空分布所控制，易使飽和-非飽和分界帶形成應力集中；皇甫凱龍等[8]應用二維潰壩數(shù)值模型分析計算潰壩水砂的流動狀態(tài)，得出尾礦壩潰決水砂與水壩具有顯著不同的流動狀態(tài)；王興華等[9]通過Fluent 數(shù)值模擬并考慮有水全潰壩、無水全潰壩、無水局部潰壩等3 種工況，研究分析尾礦庫潰壩洪水的流動特性；Zandarín 等[10]采用流固耦合有限元方法模擬了潰壩誘因，指出毛細水是影響滲透性較差壩體穩(wěn)定的關鍵因素；Yuan 等[11]利用Fluent 軟件從水砂流速、沖擊力和淹沒深度等方面分析了尾砂流的演進特征，得出了下游各位置的水砂分布；孔祥云等[12]通過自制潰壩試驗裝置研究分析尾礦庫水位變動工況下的壩體位移形變情況；姚志堅等[13]利用“等效糙率”方法代替下游構(gòu)筑物，并將該方法應用于二維潰壩模擬中；李強等[14]針對實際尾礦庫，將流固耦合與強度折減相結(jié)合，開展尾礦壩滑帶形成過程及滑帶空間分布研究；侯圣山等[15]采用FLO-2D 模型計算了降雨條件下泥石流的流動狀態(tài)及運動特征，并劃分了泥石流災害危險性分區(qū)圖。

以上針對尾礦壩穩(wěn)定性與潰壩的模擬分析多是獨立進行的，多采用實際地形概化后的二維模型，缺乏三維尾礦壩穩(wěn)定性與精細地形下潰壩過程的綜合分析，也缺乏下游構(gòu)筑物對水流滯水作用的可靠處理方法。本文以遼寧省某上游式尾礦庫為例，開展正常水位、洪水位、漫頂水位等3 種工況下尾礦壩穩(wěn)定性研究，采用等效糙率方法將下游構(gòu)筑物與地形對潰壩泥砂的阻水作用進行替代模擬，同時實現(xiàn)尾礦壩穩(wěn)定性分析與高精度地形下尾礦庫水砂兩相流三維潰壩模擬相結(jié)合，開展尾礦壩穩(wěn)定性-潰壩的綜合剖析及水砂在下游地形流動特點分析，為尾礦庫安全運行、災害分析和應對策略提供科學依據(jù)。

1 計算原理

1.1 非飽和流固耦合分析

尾礦壩內(nèi)由于浸潤線的存在，使得浸潤線以下形成飽和區(qū)，以上形成非飽和區(qū)，而庫內(nèi)水位高度變化引起的浸潤線變化也可認為是飽和-非飽和滲流過程[16]。FLAC3D中非飽和流固耦合微分方程由平衡方程、連續(xù)方程、基質(zhì)吸力方程組成[17]，平衡方程為：

式中：D——剛度張量；

ε(r)——應變張量；

r——位移向量；

χ——非飽和土參數(shù)；

pw——孔隙水壓力/Pa；

pnw——孔隙氣壓力/Pa；

γ——土體重度/（N·m-3）。

將達西定律與質(zhì)量守恒方程聯(lián)立得連續(xù)性方程：

式中：ρw——水密度/（kg·m-3）；

kw、knw——水、氣滲透系數(shù)/（m·s-1）；

γw、γnw——水、氣重度/（N·m-3）；

Z——位置水頭/m；

Sr——飽和度/%；

n——孔隙率；

H——Henry 系數(shù)。

非飽和尾礦土的基質(zhì)吸力與土體含水率有關，土-水特征曲線方程為：

式中：Sr0——初始飽和度/%；

n0——初始孔隙率。

通過尾礦庫流固耦合求解，計算出其應力場及孔壓場的分布，加入摩爾-庫倫模型，結(jié)合強度折減法，實現(xiàn)滲流與強度折減同步進行，從而計算尾礦壩穩(wěn)定性。

1.2 潰壩泥砂分析

1.2.1 VOF 模型

尾礦庫潰壩可看作潰決水砂與大氣多相流動，初始階段及潰壩階段水砂均與大氣存在接觸，而VOF 模型在多相流體模擬分析中效果最佳[18]。VOF法在計算域內(nèi)對模型中的每一相引入一個變量F，其為在網(wǎng)格內(nèi)各相流體的體積與流體空間總體積的比值，水氣界面的追蹤可通過求解F的輸運方程[19]：

式中：F=F(x,t)——空間和時間函數(shù)；

x——坐標向量，其分量為xi；

ui——速度分量。

1.2.2k-ε湍流模型

尾礦庫潰壩水砂往往以湍流的形式向下游快速無規(guī)則流動，而Fluent 提供了多種湍流模型 。本次計算參考劉學炎[20]的潰壩模擬研究成果，選取應用廣泛且適用于多相流模擬的標準k-ε湍流模型。

k-ε模型是基于能量方程和擴散速率方程得出的一種經(jīng)驗模型，其湍流動能k和耗散率ε控制方程[21]分別為：

k方程：

ε方程：

式中：G——紊動能生成率/（m·s-1）；

ρ——流體密度/（kg·m-3）；

μ——分子動力黏性系數(shù)；

μt——紊動黏性系數(shù)；

k——紊動動能/（m2·s-2）；

ε——紊動耗散率/（m2·s-2）；

σk——k的紊流普朗特數(shù)；

C1ε、C2ε——經(jīng)驗常數(shù)。

1.2.3 黏滯系數(shù)

潰壩水砂黏滯系數(shù)的影響因素眾多，本文采用水砂濃度得出的黏滯系數(shù)計算方程[22]：

式中：η——水砂動力黏度/（Pa·s）；

μ0——水的黏滯系數(shù)，常溫下 μ0為0.99；

Cv——水砂濃度/%。

經(jīng)勘察得水砂濃度為67%，則η為15.84 Pa·s。

1.2.4 最大泄砂流量

洪水漫頂造成尾礦庫從開始潰壩到穩(wěn)定的時間短暫，可按瞬潰處理[23]。瞬時垂向局部潰壩時，存在高度h1的殘壩，壩址最大泄砂流量[24]為：

式中：QM——壩址泄砂流量/（m3·s-1）；

B——潰壩前水面寬度/m；

H0——潰壩前水深/m；

g——重力加速度/（m·s-2）；

h1——殘壩高/m。

1.2.5 下游潰壩水砂流量

潰壩水砂在不同位置處的流量受距離和地形影響顯著[25]，而本次潰壩模擬既有山谷狹窄地形，又有農(nóng)田平緩地形，因此對下游重要構(gòu)筑物位置的流量分析尤為重要。下游水砂流量計算方程[26]：

式中：QLM——下游L處最大流量/（m3·s-1）；

W——庫容/m3；

L——距壩址距離/m；

υk——經(jīng)驗系數(shù)，山區(qū)河道取7.15，平原河道取3.13。

2 工程概況及相關參數(shù)

2.1 模擬區(qū)域概況

某上游式尾礦庫位于遼寧省，尾礦庫及下游地形見圖1。壩體由初期壩及后期堆積壩組成，加高增容后總壩高約102 m，壩頂長203 m，總庫容約16.95×106m3，壩體平均外坡比為1∶4。模擬區(qū)域東西方向長1.8 km，南北方向平均寬650 m。該尾礦庫下游為辦公區(qū)、選礦場、礦料堆積場、農(nóng)田、溫室以及磚廠。根據(jù)《尾礦堆積壩巖土工程技術規(guī)范》（GB 50547—2010）[27]，此尾礦庫等級標準為二等庫，尾礦壩級別為2 級。

圖1 尾礦庫及下游地形Fig.1 Topography of the tailings pond and the area downstream

2.2 尾礦壩物理力學參數(shù)

將現(xiàn)場采的原狀尾礦土試樣進行直剪速率為0.02 mm/min 的慢剪試驗，得到黏聚力及內(nèi)摩擦角。經(jīng)篩分的樣品烘干后，以含水率為5%的梯度制備尾細砂、尾粉砂、尾粉質(zhì)黏土樣品至飽和，其飽和含水率分別為35%、33%、30%。通過接觸濾紙法對樣品進行基質(zhì)吸力測定，經(jīng)VG 模型[28]擬合得到土-水特征曲線（圖2）。結(jié)合尾礦庫現(xiàn)場勘察結(jié)果，得到尾礦壩基本物理參數(shù)（表1）。

圖2 尾礦壩各土層土壤土-水特征曲線Fig.2 Soil water characteristic curve of each layer of the tailings dam

表1 尾礦壩物理力學參數(shù)Table 1 Physical and mechanical parameters of the tailings dam

3 庫水位變化對穩(wěn)定性的影響

根據(jù)前期礦山地質(zhì)鉆探資料，結(jié)合本次現(xiàn)場調(diào)查結(jié)果，建立尾礦庫穩(wěn)定性計算模型（圖3）。尾礦庫穩(wěn)定性模擬計算通過FLAC3D軟件實現(xiàn)，通過賦值孔隙水壓力將坡表及初期壩的邊界條件設為透水邊界，在進行滲流計算時通過摩爾-庫倫本構(gòu)模型對賦予物理力學參數(shù)的土體進行強度折減。為了更好地對數(shù)值模擬結(jié)果進行定量分析，對模型設置A、B、C、D 等4 個關鍵點（圖3）進行位移形變監(jiān)測。

圖3 尾礦庫模型及監(jiān)測點位置Fig.3 Tailings pond model and location of detection points

3.1 庫水位變化對滲流的影響

根據(jù)《尾礦設施設計規(guī)范》（GB 50863—2013）[29]中的上游式尾礦庫最小安全超高（洪水位與壩頂之間高差）和干灘長度（壩頂與庫內(nèi)水邊線的水平距離）限值（表2），可知研究區(qū)尾礦庫洪水位的安全超高和干灘長度分別為1.0 m 和100 m。結(jié)合實際尾礦庫的現(xiàn)場勘察結(jié)果，得到尾礦庫正常水位運行情況下的安全超高為1.3 m、干灘長度122 m。基于此，將庫內(nèi)水位高度劃分為漫頂水位（安全超高0 m）、洪水位（安全超高1 m）及正常水位（安全超高1.3 m）。通過計算得出各水位下的浸潤線形態(tài)（圖4）。

表2 上游式尾礦壩的最小安全超高和最小干灘長度Table 2 Minimum safe superelevation and minimum dry beach length of the upstream tailings dam

圖4 不同水位浸潤線位置Fig.4 Location of phreatic line at different water levels

由圖4 可知，浸潤線的埋深隨尾礦庫水位的升高而變小，由正常水位提升至洪水位，浸潤線埋深下降5～8 m，干灘長度縮短24 m；由洪水位升高至壩頂，浸潤線埋深下降5～6 m，干灘長度縮短為0 m。而漫頂水位工況下，浸潤線埋深較小，壩頂處穩(wěn)態(tài)飽和區(qū)沿坡面向下運移8 m 左右，逐步到達1 級臺階處，壩頂部分達到飽和狀態(tài)，土體的黏聚力迅速下降，極易造成壩體局部失穩(wěn)從而發(fā)生潰壩事故。

3.2 庫水位變化對剪切帶及穩(wěn)定性系數(shù)的影響

通過滲流場與應力場耦合，利用強度折減法來分析現(xiàn)場正常水位、洪水位、漫頂水位等3 種工況下壩體的穩(wěn)定性。為便于結(jié)果分析，截取尾礦壩中線位置處的穩(wěn)定性計算結(jié)果剖面，圖5 為3 種水位工況的剪切帶。

圖5 不同水位剪切帶Fig.5 Shear zones with different water levels

由圖5 可知，庫水位對剪切帶及尾礦壩穩(wěn)定性有顯著影響。剪切應變率是剪切帶隨壩體變形破壞的速率，隨水位上漲，壩體剪切應變率增大1 個數(shù)量級，表明水位變化的影響使剪切應變加快。正常水位時，剪切應變率最大值為2.95×10-6，出現(xiàn)在3 級臺階處。從正常水位升高到洪水位時，剪切應變率增大為5.78×10-5，尾礦壩穩(wěn)定系數(shù)由1.80 下降至1.32。剪切帶在縱向上不斷向壩體內(nèi)部延伸，剪切帶從最初10～12 m 深發(fā)展到15～18 m 深；在橫向上從壩頂貫穿到3 級臺階發(fā)展到從壩頂貫穿到4 級臺階，水平上增加了20 m。而水位達到壩頂時，剪切應變率進一步增大為3.32×10-4，出現(xiàn)在1 級臺階處，尾礦壩穩(wěn)定系數(shù)由1.32 下降到1.18。剪切帶急劇縮短，只在1 級臺階處貫穿。究其原因，洪水漫頂時，浸潤線沿坡面向下運移使得壩頂處尾礦土含水率遠大于其他位置土壤，達到局部飽和。

3.3 庫水位變化對壩體變形的影響

由于水位變化使得壩體內(nèi)部應力發(fā)生改變，導致壩體受應力變化發(fā)生形變，形變進一步影響尾礦壩土體結(jié)構(gòu)，使其穩(wěn)定性明顯下降。因此在計算完成后提取監(jiān)測點A、B、C、D 的位移變化曲線圖（圖6），對3 種工況下壩體變形作出定量化分析。

由圖6 可知，各水位工況下尾礦壩形變量差異顯著，總體呈先急增后平穩(wěn)趨勢。計算初期，各監(jiān)測點位移數(shù)值迅速增加，在設置大變形情況下壩體水平變化明顯，在2 000～4 000 步之間，位移曲線產(chǎn)生拐點；在計算中期時，壩體形變量在0.02～0.05 m 之間上下小幅度波動，此階段壩體未產(chǎn)生較大形變；進入計算后期，計算不斷收斂，位移曲線不再發(fā)生變化。

圖6 各監(jiān)測點位移變化曲線圖Fig.6 Displacement change curve of each monitoring point

正常水位時，監(jiān)測點A、B、C 監(jiān)測到壩體形變量，變化量值為B＞A＞C，分別為0.19，0.17，0.09 m。洪水位時，監(jiān)測點A、B、C、D 監(jiān)測到壩體形變量，變化量值為C＞B＞A＞D，分別為0.26，0.24，0.22，0.19 m。漫頂水位時，監(jiān)測點A、B、C 監(jiān)測到壩體形變量，其變化量值為A＞B＞C，分別為0.54，0.23，0.15 m。綜上，尾礦壩位移形變最大值位置應處于潛在滑面中線與坡面交點偏下處，與剪切應變率最大位置略有差異。究其原因，剪切應變率最大位置處表示壩體失穩(wěn)時優(yōu)先破壞的位置，而臨空面處于破壞處之上，在向下滑動過程中位移變化較大。

4 潰壩模擬

上述穩(wěn)定性分析結(jié)果表明，漫頂水位時壩頂產(chǎn)生局部飽和，潰壩風險最大，以此工況作為尾礦庫潰壩模擬的基礎。漫頂水位時尾礦壩潛在滑帶由壩頂貫穿到1 級臺階，此時穩(wěn)定系數(shù)最低，僅為1.18。規(guī)范未提供漫頂穩(wěn)定系數(shù)，洪水漫頂出現(xiàn)的概率完全高于地震荷載特殊工況。由于穩(wěn)定系數(shù)1.18 略高于特殊工況穩(wěn)定系數(shù)1.15，因此可以認為尾礦壩處于欠穩(wěn)定狀態(tài)。根據(jù)圖5，將壩頂?shù)? 級臺階的高度作為尾礦庫潰壩模擬潰口的高度，其值為13.2 m。將FLAC3D的計算結(jié)果調(diào)入到Ansys-Fluent 中，進行尾礦庫潰壩水砂在下游流動狀態(tài)的模擬計算。

4.1 模型建立及邊界條件

采用Rhino 和Fluent 軟件進行建模及分析計算。在Rhino 中通過尾礦庫的等高線建立復雜潰壩模擬地形圖。將建立好的模型導入Fluent 軟件中進行計算，使?jié)螖?shù)值模擬結(jié)果與實際情況更相符。在Fluent 中設置模型邊界條件，將模型上部空氣界面設置為壓力入口，四周界面為壓力出口，山體地形和壩體為墻體[30]。水砂受地形摩擦力影響，而不同地形對水砂的摩擦影響不同，本文通過等效糙率法來控制地形對水砂的摩擦影響[31]，糙率值的取值見表3。

表3 地形對水砂摩擦影響的糙率值Table 3 Influence of topography on the water sand friction roughness value

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

4.2.1 潰壩水砂淹沒范圍

根據(jù)Fluent 計算結(jié)果，得出尾礦庫發(fā)生潰壩后潰決水砂在不同時刻流動演進情況的體積分數(shù)（Volume fraction）云圖（圖7）。由圖7 可知，t為0 s 時，尾礦庫內(nèi)水砂和堆積壩處于相對靜止狀態(tài)，隨著時間推移，堆積壩壩頂出現(xiàn)潰口瞬間潰決。0～100 s 時，潰口附近尾砂由于水流拖曳力作用沖刷形成初期水砂混合流，后方尾砂隨之開始滑動；100 ～300 s 時，尾礦庫內(nèi)水砂在重力勢能的作用下快速流出；300～700 s 為水砂流在山谷內(nèi)的匯流演進階段，300 ～500 s 時潰壩水砂淹沒下游辦公區(qū)，潰壩水砂的前端來到礦料堆積區(qū)，水砂速度在高差明顯的山谷內(nèi)將會進一步變大；700 s 以后，潰壩水砂在農(nóng)田等地勢平坦處向四周發(fā)散流動，700～900 s 時，前端水砂流到低勢能區(qū)，水砂淹沒范圍更加廣泛，但水砂伴隨尾砂料之間的碰撞摩擦其速度將會降低；1 300 s 之后，水砂趨于停止。

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

4.2.2 潰壩水砂速度演進

潰壩水砂在山谷及農(nóng)田兩處的流動狀態(tài)分別為匯聚流動和發(fā)散流動，山谷及農(nóng)田的水砂速度云圖見圖8。由圖8 可知，500 s 時，山谷內(nèi)地勢陡峭，云圖顏色由堆積壩上的藍色轉(zhuǎn)為山谷內(nèi)的黃綠色，水砂流速進一步增大。900 s 時，農(nóng)田處地勢平緩，云圖顏色由山谷內(nèi)的深紅色轉(zhuǎn)為黃綠色，水砂發(fā)散流動，流速降低。潰壩水砂流速整體狀態(tài)呈“增大—增大—減小”趨勢，主要表現(xiàn)為在堆積壩受重力作用從0 逐漸增大，流經(jīng)山谷受地形地勢影響流速進一步增大，到達農(nóng)田時水砂發(fā)散流動，流速開始降低，山谷內(nèi)辦公區(qū)附近流速最大。

圖8 潰壩水砂速率云圖Fig.8 Cloud chart of water and sand velocity of dam break

4.2.3 下游水砂流量

此上游式尾礦庫潰壩水砂最終淹沒范圍東西方向長1 362 m，南北方向平均寬440 m，淹沒面積約為0.6 km2。初期壩東方向的辦公區(qū)、礦料堆積場及下游農(nóng)田處于主要淹沒區(qū)內(nèi)，而選礦廠房和料場地勢較高，基本不受潰壩影響。由式（9）（10）計算可知，此尾礦庫最大泄砂流量為2 905.43 m3/s，下游470 m 處辦公區(qū)出現(xiàn)最大水砂流量為2 872.64 m3/s，下游1 020 m 處農(nóng)田出現(xiàn)最大水砂流量為2 651.62 m3/s。辦公區(qū)與農(nóng)田相距550 m，最大水砂量相差221.02 m3/s，由式（10）可知，隨著L的增大，QLM逐漸減小。因此，在尾礦庫選址時，距下游重要建筑物距離應作為優(yōu)先考慮因素。

5 結(jié)論

（1）尾礦壩浸潤線埋深隨庫水位升高而不斷變小，漫頂水位時，壩頂浸潤線沿坡面向下運移，使得壩頂部分達到局部飽和，土體黏聚力急劇下滑，造成優(yōu)先滑動破壞現(xiàn)象。

（2）從正常水位到洪水位時，剪切帶在縱向上不斷向壩體內(nèi)部延伸；在橫向上從壩頂貫穿到3 級臺階發(fā)展到從壩頂貫穿到4 級臺階，洪水漫頂時，浸潤線沿坡面向下運移使得壩頂局部達到飽和、優(yōu)先發(fā)生破壞，剪切帶急劇縮短。

（3）由FLAC3D穩(wěn)定性分析軟件計算出失穩(wěn)時的潛在滑面，以此作為潰口，利用Rhino 建立尾礦庫下游精細地形，將穩(wěn)定性分析結(jié)果作為Fluent 中潰壩模擬初始條件并進行計算，得出水砂的流動狀態(tài)：在山谷狹窄地形下，水砂具有流速快，水砂流量大，淹沒范圍小的匯集流動；在農(nóng)田等平坦地區(qū)，水砂具有流速慢，水砂流量小，淹沒范圍廣的發(fā)散流動。

（4）通過設置地表糙率，能夠有效降低模擬結(jié)果與實際情況的偏差；發(fā)生潰壩后，潰口水砂在重力勢能和堆積壩體的摩擦共同作用下迅速向下游流動，由于山谷內(nèi)高差影響，此時水砂流速最大；隨后，潰口的水位不斷降低，受下游地勢平緩和地表面的粗糙程度影響，水砂不會向下游產(chǎn)生較大區(qū)域的流動。