寧 偉 李 瑜 吳蒙蒙 楊 超

（1.洛陽(yáng)欒川鉬業(yè)集團(tuán)股份有限公司；2.欒川縣環(huán)境保護(hù)局；3.河海大學(xué)水利水電學(xué)院）

隨著我國(guó)基礎(chǔ)設(shè)施建設(shè)行業(yè)的大力發(fā)展，對(duì)各種金屬和非金屬礦產(chǎn)資源的需求也迅速增加。在礦產(chǎn)資源的開發(fā)和利用中，不可避免地產(chǎn)生大量以漿狀形式排出的廢棄礦渣，一般的礦山企業(yè)都會(huì)選擇建設(shè)尾礦庫(kù)堆放，來(lái)減少尾礦渣給環(huán)境帶來(lái)的污染。尾礦庫(kù)是一個(gè)重大的危險(xiǎn)源，一旦發(fā)生潰壩，產(chǎn)生的后果相當(dāng)嚴(yán)重［1］。因此，進(jìn)行尾礦庫(kù)潰壩數(shù)值模擬研究，分析尾砂演進(jìn)與泥石流到達(dá)時(shí)間特性，可以為評(píng)價(jià)尾礦庫(kù)建設(shè)方案的安全性和合理性提供依據(jù)。

尾礦庫(kù)潰壩是一個(gè)復(fù)雜的水沙耦合過(guò)程，其研究成果豐碩。尹光志等［2］運(yùn)用數(shù)值模擬和室內(nèi)模擬的方法，對(duì)云南玉溪某尾礦庫(kù)進(jìn)行潰壩砂流演進(jìn)規(guī)律模擬，探討泥沙在下游的淹沒范圍。袁兵等［3］通過(guò)研究幾個(gè)尾礦庫(kù)的潰決案例，考慮尾礦渣的物理性質(zhì)和流動(dòng)過(guò)程中的變形來(lái)建立數(shù)學(xué)模型，對(duì)尾礦庫(kù)潰壩災(zāi)害進(jìn)行研究。林江等［4］等利用非恒定流輸砂方程模擬潰壩沙流的演進(jìn)過(guò)程，建立了一種新的模型來(lái)同時(shí)模擬下泄水流和尾礦砂的運(yùn)動(dòng)，并將其運(yùn)用在云南某工程中。DL Breach模型在國(guó)內(nèi)土石壩潰壩模擬分析中的應(yīng)用愈來(lái)愈廣泛，但鮮有人將其運(yùn)用到尾礦庫(kù)潰壩上。周清勇等［5］將Breach模型成功應(yīng)用于七一水庫(kù)風(fēng)險(xiǎn)分析；楊德瑋等［6］將Breach模型與MIKE模型進(jìn)行耦合分析，成功應(yīng)用于鄭州常莊水庫(kù)；李相南等［7］同樣對(duì)比分析了Breach模型與DBIWHR模型在唐家山堰塞湖潰決中的反演分析，認(rèn)為2種模型均可較好地反演唐家山潰決洪水的過(guò)程。

本研究基于DL Breach模型，利用數(shù)值模擬的方法，探索沙流運(yùn)動(dòng)的規(guī)律，計(jì)算尾礦庫(kù)下游的淹沒范圍、淹沒深度，預(yù)測(cè)尾礦壩潰決可能引起的災(zāi)害，對(duì)礦區(qū)居民的安全保障和礦區(qū)環(huán)境的保護(hù)有十分重要的指導(dǎo)意義。

1 潰壩數(shù)學(xué)模型

1.1 模型分類與選擇

參考土石壩的潰壩模型，尾礦庫(kù)潰壩的數(shù)學(xué)模型也可以分為3類：第一類是基于參數(shù)的模型，該模型通過(guò)利用已有的尾礦庫(kù)工程的潰壩數(shù)據(jù)，根據(jù)經(jīng)驗(yàn)確定尾礦庫(kù)潰壩的重要參數(shù)；第二類基于尾礦庫(kù)潰壩機(jī)理的簡(jiǎn)化數(shù)學(xué)模型，該模型利用水動(dòng)力學(xué)相關(guān)的推導(dǎo)，建立了一個(gè)簡(jiǎn)化的數(shù)學(xué)模型并用來(lái)模擬潰壩沙流的下泄流量過(guò)程線和潰壩沙流演進(jìn)過(guò)程；第三類是基于尾礦庫(kù)潰壩機(jī)理的精細(xì)數(shù)學(xué)模型，該模型利用水動(dòng)力壩料沖蝕方程來(lái)建立潰壩的三維數(shù)學(xué)模型，可以更加精細(xì)地模擬尾礦庫(kù)的潰壩過(guò)程［8］。DL Breach模型屬于第二類潰壩模型，該模型基于數(shù)學(xué)機(jī)理，假定潰口形狀為梯形，預(yù)測(cè)潰口發(fā)展趨勢(shì)和潰壩導(dǎo)致的泄流水位過(guò)程線，通過(guò)輸入大壩幾何尺寸、大壩寬高值和材料內(nèi)外屬性值等重要參數(shù)，來(lái)求得含有壩體潰決的各類相關(guān)參數(shù)以及大壩潰口流量過(guò)程曲線。本研究以澀草湖尾礦庫(kù)為例，基于DL Breach模型獲取尾礦庫(kù)潰壩的下泄沙流流量過(guò)程線，模擬下泄沙流的演進(jìn)過(guò)程。

1.2 DL Breach模型

DL Breach模型是一個(gè)基于數(shù)學(xué)機(jī)理的預(yù)測(cè)潰口特征（尺寸、形成時(shí)間）和潰決土壩引起的泄流水位過(guò)程線的數(shù)學(xué)模型，該模型原理基于水力學(xué)、泥沙運(yùn)動(dòng)、土力學(xué)、大壩的幾何屬性和材料屬性及水庫(kù)屬性［9］。DL Breach模型將漫頂潰決分為2個(gè)階段：第一階段為強(qiáng)烈決口或侵蝕階段，決口水流為超臨界，受上游控制。第二階段為一般決口或進(jìn)口演變階段，其水流為亞臨界流，受上游或上下游共同控制。在第一階段，使用堰流方程計(jì)算決口流量；在第二階段，使用Keulegan公式計(jì)算決口流量。該模型將漫頂潰壩模式下的潰口橫截面近似為梯形，潰口縱截面近似為平頂。該模型已經(jīng)通過(guò)使用50組大壩潰決數(shù)據(jù)和4組實(shí)驗(yàn)案例數(shù)據(jù)來(lái)對(duì)比，計(jì)算的峰值潰壩流量、潰口寬度等參數(shù)與實(shí)測(cè)數(shù)據(jù)基本一致。因此將DL Breach模型應(yīng)用于尾礦庫(kù)潰壩下泄沙流的計(jì)算中，同樣具有很重要的參考意義。

1.2.1 潰口流量計(jì)算

潰口的流量計(jì)算可以分為2個(gè)階段，第一個(gè)階段可以采用寬頂堰公式來(lái)計(jì)算流量：

式中，H=z s-z b，m，z s是上游水頭，m，z b是潰口的底部高程，m；b為潰口的底部寬度，m；m為下游潰口邊坡的坡度；系數(shù)c1=1.7，c2=1.3；k sm為堰流尾水淹沒修正系數(shù)；z t是尾水位的高程，m。

潰決的第二階段采用Keulegan公式計(jì)算，該公式是簡(jiǎn)化的恒定非均勻流能量方程，可以考慮河道的收縮或擴(kuò)張所引起的水頭損失：

式中，L是下泄沙流方向的距離，m；Q是潰口處的流量，m3/s；A是尾礦庫(kù)潰口的面積，m2；R是水力半徑，m；B是潰口處的水面寬度，m；n是曼寧系數(shù)；g是重力加速度，m/s2；λe n和λex分別是下泄沙流入口和出口處的水頭損失，m；ρ是水的密度，kg/m3；ρa(bǔ)是空氣密度，kg/m3；C d是風(fēng)阻系數(shù)；U w in是風(fēng)速，m/s；θwin是潰口中心線與風(fēng)向的夾角，（°）。

1.2.2 輸砂模型

尾礦庫(kù)潰壩的下泄沙流中存在著懸移質(zhì)與河床質(zhì)泥沙，總的輸砂率qb可以表示為

式中，R為河槽水力半徑，m；ωs為泥沙流速，m/s；U為水流速，m/s；q*為懸移質(zhì)輸砂率；q b*為河床質(zhì)泥沙單位時(shí)間單寬輸砂率，m2/s；d為泥沙粒徑，mm；τc=為泥沙的飽和容重，N/m3，γ為泥沙的天然容重n為河床的曼寧糙率系數(shù)，τb為河床剪應(yīng)力，Pa。

1.2.3 潰口擴(kuò)展過(guò)程

尾礦庫(kù)潰口邊坡的穩(wěn)定程度與尾礦砂的性質(zhì)、堆積高度等因素密切相關(guān)，可以采用滑楔法來(lái)模擬潰口邊坡的發(fā)展（圖1）。

2 潰壩沙流運(yùn)動(dòng)計(jì)算原理

2.1 尾砂流變模型

尾礦庫(kù)是一種特殊的水工建筑物，當(dāng)庫(kù)容較大、堆積壩較高時(shí)，庫(kù)內(nèi)尾砂具有較高的勢(shì)能，一旦潰決水流挾帶尾礦砂從潰口處沖出。尾礦庫(kù)下泄的沙流是一種散粒體與水復(fù)合材料，組分復(fù)雜多樣。從尾礦庫(kù)排放的尾砂造成的下泄沙流性質(zhì)類似泥石流，潰壩初期由于水量較多，其性質(zhì)更接近稀性泥石流。由泥石流引起的土壤流動(dòng)可以被認(rèn)為是流體和散粒體之間特殊的運(yùn)動(dòng)形式。它可以用與流體流動(dòng)類似的動(dòng)態(tài)方程和連續(xù)方程來(lái)描述。當(dāng)沙流流動(dòng)時(shí)具有初始屈服應(yīng)力，表現(xiàn)出非牛頓流體的特性，根據(jù)工程經(jīng)驗(yàn)，尾砂下泄沙流符合Bingham流動(dòng)模式，Bingham理想塑性體的切應(yīng)力與變形速率成線性函數(shù)關(guān)系，但并未經(jīng)過(guò)原點(diǎn)，存在屈服應(yīng)力。不同流體的剪切應(yīng)力與變形速率關(guān)系如圖2所示。Bingham流體在運(yùn)動(dòng)的過(guò)程中，密度、黏度以及屈服應(yīng)力等均不隨時(shí)間和空間而改變，因此在計(jì)算過(guò)程中，把下泄沙流假設(shè)為Bingham流體，其運(yùn)動(dòng)控制方程僅需滿足連續(xù)性定理、動(dòng)量守恒和能量守恒，可以簡(jiǎn)化運(yùn)動(dòng)控制方程和數(shù)值模擬的過(guò)程。同時(shí)假設(shè)尾礦漿在流動(dòng)過(guò)程中，尾礦砂與水充分混合，兩者之間沒有相對(duì)運(yùn)動(dòng)速度，這也符合下泄沙流的一般特征。

Bingham方程是根據(jù)Bingham塑性流體的流變曲線寫出的，其形式為

式中，τ為剪切應(yīng)力，Pa；τ0為屈服應(yīng)力，Pa；ηp為塑性黏度為速度梯度，s-1。

符合Bingham方程的流體稱為Bingham流體，塑性流體也稱Bingham流體。應(yīng)力τ0和塑性黏度ηp是公式中2個(gè)主要的流變參數(shù)。屈服應(yīng)力是指Bingham流體為實(shí)現(xiàn)層流流動(dòng)所需要克服的動(dòng)態(tài)結(jié)構(gòu)力，而塑性黏度則反映了層流狀態(tài)下的Bingham流體在處于動(dòng)平衡狀態(tài)時(shí)內(nèi)摩擦力的大小。

Bingham流體的表觀黏度可表示為

2.2 控制方程

流動(dòng)控制方程建立在基于Boussinesq和流體靜壓假定的不可壓縮雷諾平均N-S方程的解決方案的基礎(chǔ)上?？紤]到上述針對(duì)尾礦壩潰壩沙流的基本假定，可近似成下述平面深度平均方程組來(lái)描述。

連續(xù)方程：

動(dòng)量方程：

式中，t為時(shí)間，s；x,y為笛卡爾坐標(biāo)系坐標(biāo)，m；η為水面高度，m；d為靜水深，m；h為總水深，m；f為科里奧利參數(shù)；g為重力加速度，m/s2；ρ為密度，kg/m3；s xx，s xy，s yx和s yy為輻射應(yīng)力張量的分量，Pa；ρ0為水的密度，kg/m3；S為點(diǎn)源的流量大小，m/s；uˉ，vˉ為x，y方向上的速度分量，m/s。

3 工程實(shí)例

3.1 概況

澀草湖尾礦庫(kù)位于三聯(lián)村澀草湖下游的鉆天道與大干溝內(nèi)，初期壩為碾壓堆石壩，設(shè)計(jì)總庫(kù)容5 742.5萬(wàn)m3，尾礦庫(kù)總壩高為199.0 m，平均堆積邊坡1∶5.0。鉆天道在初期壩下游約550 m處，原居住4處村宅，現(xiàn)房屋已拆除，人員已搬遷。初期壩下游1 680，1 960，2 560 m處有3座橋，分別定義為1#橋、2#橋、3#橋，這3座橋附近分布有17戶村宅。

3.2 潰壩計(jì)算參數(shù)確定

澀草湖尾礦庫(kù)初期壩上、下游邊坡為1∶2.0，壩高79 m，壩頂寬4 m，壩軸線長(zhǎng)約165 m，后期堆積壩平均堆積邊坡1∶5.0，尾礦堆高120 m，總壩高199 m，設(shè)計(jì)總庫(kù)容5 742.5萬(wàn)m3，有效庫(kù)容4 708.8萬(wàn)m3。

根據(jù)提供的工程地質(zhì)勘察資料和設(shè)計(jì)資料，結(jié)合現(xiàn)狀尾礦庫(kù)滲流場(chǎng)反演分析成果，在尾礦庫(kù)剖面上，可以將其劃分為5個(gè)區(qū)域，具體分區(qū)和參數(shù)如圖3和表1所示。

?

3.3 邊界條件及初始條件

邊界條件包括上游邊界和下游邊界。上游的邊界條件為下泄沙流的流量過(guò)程線，由DL Breach模型算出；下游邊界條件為開邊界，位置在距離初期壩壩頂3 400 m處，其余模型邊界作為陸地邊界，取法線方向速度分量為零。

尾礦庫(kù)下游河道中水流的流量很小，水深很淺，因此忽略其影響，認(rèn)為初始水深為0 m，尾礦庫(kù)下游的初始條件取天然的地形高程。

3.4 計(jì)算工況與流量過(guò)程

根據(jù)相關(guān)國(guó)家規(guī)范，并參考類似工程經(jīng)驗(yàn)，澀草湖尾礦庫(kù)可能發(fā)生的潰壩情況為洪水漫頂潰壩工況。漫頂潰壩是大壩潰壩的最主要形式，潰口流量過(guò)程線由DL Breach模型計(jì)算得出，由圖4可知，DL Breach模型得出的潰口最大流量為246.33 m3/s。

3.5 潰壩影響預(yù)測(cè)

3.5.1 淹沒范圍預(yù)測(cè)

本模型范圍內(nèi)沿程溝谷狹長(zhǎng)且較平緩，斷面基本呈V字形，溝谷地表多樹木，兩側(cè)山體高。由DL Breach模型計(jì)算出的潰口流量過(guò)程線可知，整個(gè)潰壩流程為3 600 s，選取幾個(gè)典型時(shí)刻：T=600，1 200，1 800，2 800 s，下泄沙流不同時(shí)刻的淹沒范圍見圖5。

3.5.2 下游淹沒深度

選取下游區(qū)域有代表性的幾個(gè)敏感點(diǎn)，計(jì)算下泄沙流在這些地方的淹沒深度，并繪制了如圖6的潰壩沙流淹沒深度隨時(shí)間變化的曲線。

從計(jì)算結(jié)果來(lái)看，尾砂會(huì)積聚在下游鉆天道內(nèi)，下游空地、溪溝、灌木等會(huì)被淹沒；鉆天道附近居民已經(jīng)搬遷，一直到初期壩下游1.5 km以內(nèi)都沒有重要建筑物，該范圍內(nèi)沙流淤積面積大，擴(kuò)散范圍廣。截滲壩和攔沙壩上游（鉆天道溝谷內(nèi)）會(huì)有大量的泥石流淤積，而其下游淹沒深度小，樹木、植被也對(duì)泥石流的演進(jìn)具有明顯的阻礙作用，危害較小。

在整個(gè)沙流運(yùn)動(dòng)的過(guò)程中，尾砂在下游的堆積體的深度也隨著距離的增大而減小，距離初期壩越遠(yuǎn)，淹沒深度越小。潰壩后水深總體上先增大后減小，大致可分為3個(gè)階段：增長(zhǎng)階段、穩(wěn)定階段、衰減階段，這也反映了沙流運(yùn)動(dòng)的過(guò)程。尾砂行進(jìn)過(guò)程中，堆積深度由0 m開始逐漸增加，增大至最大深度后，繼續(xù)向下游流動(dòng)，此后堆積深度逐漸減小，并趨于平緩，這符合沙流演進(jìn)的一般規(guī)律。鉆天道原居民點(diǎn)處的淹沒深度最大可達(dá)到2.49 m，根據(jù)計(jì)算結(jié)果可以發(fā)現(xiàn)，距離初期壩越近，尾礦砂的堆積深度越深，而不會(huì)隨著下泄沙流的運(yùn)動(dòng)往下游推進(jìn)從而堆積得更高，這也與尾砂這一特殊材料的性質(zhì)有關(guān)。對(duì)于漫頂潰壩，洪水漫過(guò)壩頂沖刷壩體表面，這種情況下尾礦壩潰口緩慢發(fā)展，按照DL Breach模型計(jì)算出的流量過(guò)程曲線，下泄沙流的流量先增大再減小，數(shù)值模擬計(jì)算結(jié)果表明，該模型不僅適用于普通的土石壩，也適用于尾礦壩這一特殊的水工建筑物。利用DL Breach模型來(lái)預(yù)測(cè)尾礦庫(kù)潰壩后下泄沙流的演進(jìn)過(guò)程是合理可靠的。

4 結(jié)語(yǔ)

（1）結(jié)合具體的工程案例，利用數(shù)值模擬的方法，計(jì)算某尾礦庫(kù)潰壩后下泄沙流的演進(jìn)過(guò)程，得到了下游淹沒范圍、各敏感點(diǎn)淹沒深度以及流量、潰壩沙流抵達(dá)時(shí)間等關(guān)鍵數(shù)據(jù)，符合尾礦沙流運(yùn)動(dòng)的一般規(guī)律。

（2）實(shí)例計(jì)算結(jié)果：潰壩下泄沙流的堆積深度總體上先增大后減小，增大到最大深度后，沙流繼續(xù)向下游流動(dòng)，此后堆積深度慢慢減??；下泄沙流在鉆天道溝谷內(nèi)淹沒深度較大，淤積較多，最大堆積深度達(dá)2.48 m，但此處沒有民宅，因此不涉及下游群眾的安問(wèn)題；隨著與初期壩的距離逐漸變大，下游的3#橋處的最大淹沒深度不超過(guò)1 m，沙流對(duì)下游主要建筑物影響較小，整體安全風(fēng)險(xiǎn)可控，不會(huì)威脅到村民的生命安全。

（3）根據(jù)利用DL Breach模型進(jìn)行計(jì)算的尾礦庫(kù)潰壩數(shù)值模擬，可以預(yù)測(cè)尾礦壩潰決可能引起的災(zāi)害，對(duì)尾礦庫(kù)的安全運(yùn)行具有指導(dǎo)意義。