劉嘉欣，閻志坤，鐘啟明,3，陳亮，單熠博

(1. 河海大學(xué)土木與交通學(xué)院，江蘇 南京，210098；2. 南京水利科學(xué)研究院巖土工程研究所，江蘇 南京，210024；3. 水利部土石壩破壞機理與防控技術(shù)重點實驗室，江蘇 南京，210029)

我國采礦業(yè)在帶來巨大經(jīng)濟效益的同時，也制造了大量廢料，其中最主要的是尾礦。目前對尾礦的主流處置方法是筑壩存儲，存儲形成的結(jié)構(gòu)物統(tǒng)稱為尾礦庫或尾礦壩。據(jù)統(tǒng)計，全國共有14 217座尾礦庫，其中鐵礦的尾礦庫最多[1]。尾礦庫運行狀況與周圍環(huán)境及居民生活關(guān)系密切，隨著社會經(jīng)濟的發(fā)展，世界各國也越發(fā)重視尾礦庫的建設(shè)和安全管理。據(jù)統(tǒng)計[2]，從2001 年至2015年，我國共發(fā)生尾礦庫事故99 起，其中由洪水漫頂引發(fā)的潰壩事故占潰壩總數(shù)的35%。2008年9月8日，山西省襄汾縣發(fā)生特別重大潰壩事故，造成277 人死亡、4 人失蹤、33 人受傷，直接經(jīng)濟損失達9 619.2 萬元[3]；2010 年9 月21 日，廣 東 省 茂名市信宜紫金礦業(yè)有限公司發(fā)生尾礦庫潰壩事故，造成22 人死亡，523 戶房屋倒塌，815 戶房屋受損[4]；2019年1月25日，巴西米納斯吉拉斯州布魯馬迪紐市發(fā)生鐵礦石尾礦庫潰壩事故，造成165人死亡、160 人失蹤[5]。目前，對尾礦庫潰壩的研究大多借鑒土石壩潰壩的研究成果[6]，但尾礦庫的壩料特性和構(gòu)筑方式均與土石壩的不同，尾礦庫潰壩時的潰口下泄物固相含量較高，且其失事概率也比水庫大壩的失事概率高10 倍以上[7]。現(xiàn)有的針對尾礦庫的模型試驗研究主要集中在尾礦庫失穩(wěn)和潰壩后泥石流演進等方面[8-12]，有關(guān)漫頂潰壩洪水流量變化過程和潰口發(fā)展過程的研究較少。數(shù)值模擬也主要集中在對壩體穩(wěn)定性的研究[13]和潰壩事故后流體運動規(guī)律的研究[14-16]，對潰壩過程的精細(xì)化模擬方法的研究較少。

綜上所述，我國尾礦庫存在數(shù)量多、事故頻發(fā)的現(xiàn)狀，且人們對尾礦庫的潰壩機理的認(rèn)識不夠清晰，考慮流固耦合的尾礦庫潰壩過程數(shù)值模擬方法較少。為保護環(huán)境、保障人民生命財產(chǎn)安全，有必要深入研究尾礦庫的漫頂潰壩機理，建立能夠有效模擬潰壩過程的精細(xì)化數(shù)值模擬方法，提高潰口流量過程和潰口發(fā)展過程的預(yù)測精度，以期為科學(xué)制定應(yīng)急預(yù)案提供理論與技術(shù)支撐。

1 尾礦庫漫頂潰壩機理

物理模型試驗是揭示尾礦庫潰壩機理的有效方式，但尾礦庫潰壩過程的模擬涉及材料力學(xué)、巖土力學(xué)、流體力學(xué)等多個學(xué)科，需綜合考慮多個物理量的相似關(guān)系，在技術(shù)上存在一定難度，故有關(guān)尾礦庫潰壩機理與潰壩過程的模型試驗研究仍處于起步階段。張紅武等[17]認(rèn)為對于尾礦庫潰壩模型，設(shè)計時應(yīng)考慮水流重力相似、水流阻力相似、水流挾砂相似、尾礦料懸移相似、河床變形相似及尾礦料起動相似，并提出了模型砂的選擇及模型試驗方法。武立功等[18]按1∶100的縮尺比探究不同粒徑尾砂筑壩的尾礦庫發(fā)生漫頂潰決時的潰決機理。趙懷剛等[19]將某尾礦庫及下游河道整體按1∶400 比例進行縮尺，研究山地復(fù)雜地形條件下尾礦庫潰決過程及潰后泥石流的演進規(guī)律。但現(xiàn)有的模型試驗縮尺比例較小，存在重力不相似等問題[20]，使得試驗結(jié)果與真實潰決機理存在偏差。

離心模型試驗通過高速旋轉(zhuǎn)的離心機所產(chǎn)生的離心力大幅提高模型的應(yīng)力水平，使模型試驗結(jié)果與原型結(jié)果更為接近，因此，離心模型試驗在巖土工程中應(yīng)用較多。陳生水等[21-22]基于南京水利科學(xué)研究院研制的潰壩離心模型試驗系統(tǒng)開展?jié)坞x心模型試驗，揭示了尾礦庫的漫頂潰壩機理和潰口流量過程線，并得到了潰壩過程中的潰口演化規(guī)律。

1.1 模型試驗設(shè)計

1.1.1 潰壩離心模型試驗系統(tǒng)

離心機高速旋轉(zhuǎn)所產(chǎn)生的N倍重力加速度場，使得土工構(gòu)筑物幾何尺寸縮小N倍后，模型和原型的應(yīng)力水平相同。模型與原型之間通過相似準(zhǔn)則進行連接，相似準(zhǔn)則是用模型試驗結(jié)果反推原型力學(xué)行為的理論依據(jù)。潰壩是水土耦合作用過程，需要依據(jù)能量守恒定律建立壩料和潰壩水流的離心試驗相似準(zhǔn)則，本試驗中模型參數(shù)與原型參數(shù)相似比見表1。

表1 潰壩離心模型試驗中模型參數(shù)與原型參數(shù)相似比Table 1 Similarity criterion between model and prototype parameters in centrifugal model test of dam breach

以南京水利科學(xué)研究院大型土工離心機(NHRI-400 GT)為基礎(chǔ)構(gòu)建試驗系統(tǒng)，該系統(tǒng)主要由離心機大流量水流控制系統(tǒng)、潰壩試驗專用模型箱以及數(shù)據(jù)采集與視頻監(jiān)控系統(tǒng)組成(見圖1[23])。

圖1 NHRI-400 GT潰壩離心模型試驗系統(tǒng)示意圖[23]Fig.1 Schematic diagram of NHRI-400 GT dam breach centrifugal model test system[23]

離心機大流量水流控制系統(tǒng)由貯水箱、接水環(huán)、伺服水閥、流量計和管路等組成(見圖2)。該系統(tǒng)具有如下特點：1)在離心機大臂上設(shè)置接水環(huán)，解決了固定件與旋轉(zhuǎn)件之間的磨損和滲漏問題，保證了高重力場下接水環(huán)內(nèi)部水流的受力均勻，實現(xiàn)了水流從地面環(huán)境到高重力場之間的平穩(wěn)過渡；2)研發(fā)伺服水閥流量控制裝置可實現(xiàn)潰壩過程中入庫流量的準(zhǔn)確控制，實測流量與目標(biāo)流量最大相對誤差控制在1%以內(nèi)。

圖2 離心機大流量水流控制系統(tǒng)Fig.2 Large flow control system for centrifuge

潰壩試驗專用模型箱的特點在于將薄壁量水堰內(nèi)嵌于模型箱端板，實現(xiàn)了潰口流量的精確測量。模型箱內(nèi)部有效長×寬×高為1 200 mm×400 mm×800 mm，一側(cè)設(shè)置有機玻璃，薄壁堰口下游面與上游直立面成30°夾角，堰口寬度為280 mm，堰頂至模型箱底的距離為200 mm(見圖3)。為了在超重力環(huán)境中測量堰前水頭，在模型箱中放置水壓力傳感器(見圖4)。采用動態(tài)數(shù)據(jù)采集與傳輸技術(shù)以及視頻實時監(jiān)測系統(tǒng)和基于粒子圖像測量技術(shù)(PIV)的圖像測量系統(tǒng)，構(gòu)建數(shù)據(jù)采集傳輸與圖像監(jiān)控系統(tǒng)，以實現(xiàn)潰壩過程的多角度監(jiān)測和攝錄。

圖3 薄壁量水堰尺寸Fig.3 Size of thin-wall water measuring weir

圖4 水壓力傳感器Fig.4 Water pressure sensors

1.1.2 模型試驗參數(shù)

本次離心模型試驗所用尾礦料取自馬鋼集團姑山礦業(yè)公司青山林場尾礦庫，通過篩分試驗測得模型尾礦料中各粒組質(zhì)量分?jǐn)?shù)，見表2(其中d為粒徑)。模型試驗壩料級配曲線如圖5 所示，模型試驗尾礦料的基本物理力學(xué)特性參數(shù)見表3。

表3 離心模型試驗尾礦料基本特性參數(shù)Table 3 Basic characteristic parameters of tailings for centrifugal model test

圖5 模型試驗壩料級配曲線Fig.5 Gradation curve of model test dam material

表2 模型試驗壩料各粒組質(zhì)量分?jǐn)?shù)Table 2 Each granule group of dam material in model tests %

尾礦庫不存在顯著的上游壩坡，且?guī)靺^(qū)內(nèi)尾礦料自然沉積時間較短，固結(jié)較弱，以水砂混合為主，庫容量小?？紤]到尾礦庫的材料、結(jié)構(gòu)特點以及試驗研究內(nèi)容，模型主要考慮干灘面、壩頂和下游壩坡部分。

離心模型試驗主要研究漫頂水流對尾礦庫壩頂和下游坡潰口的沖刷影響，測量潰口發(fā)展過程和潰口流量變化過程相關(guān)參數(shù)等。離心模型試驗設(shè)計A，B 和C 這3 組模型，尾礦庫試驗壩參數(shù)見表4。

表4 尾礦庫潰壩離心模型試驗壩參數(shù)Table 4 Parameter values of tailings dams breaching for centrifugal model tests

模型A 不設(shè)初始潰口，用來觀察潰口的形成和發(fā)展過程；模型B在壩頂靠近模型箱有機玻璃一側(cè)開挖一個頂寬為2 cm、底寬為1 cm、深度為1 cm 的倒梯形凹槽作為初始潰口，用來觀測潰口的垂向發(fā)展過程；模型C在壩頂中間位置設(shè)置初始倒梯形潰口，初始潰口頂寬為4 cm、底寬為2 cm、深度為1 cm，觀測潰口橫向擴展過程。在模型上、下游設(shè)置圓開孔不銹鋼板透水支撐邊界。模型試驗壩干密度為1.80 g/cm3，分層填筑，填筑完成后在上游庫區(qū)內(nèi)加水，保持水位在干灘面高度并靜置24 h形成初始滲流場。下面選取模型B進行試驗結(jié)果分析。

1.2 潰口演化規(guī)律

模型B成模后引導(dǎo)漫頂水流從右側(cè)初始潰口處的壩頂開始沖刷，試驗在50g(g為重力加速度)離心加速度條件下開展，通過對試驗視頻錄像資料進行分幀處理，得到尾礦庫潰口發(fā)展過程，如圖6所示，尾礦庫剖面形態(tài)演化過程如圖7所示。由圖7可以看出：潰口的最終深度為23.81 cm。

圖7 尾礦庫潰壩過程剖面形態(tài)演化示意圖Fig.7 Profile morphologic evolution diagram of tailings dam breach process

根據(jù)尾礦庫潰壩模型試驗觀測資料，可將潰壩過程分為如下3個階段。

1)下游坡沖蝕(見圖6(a)和6(b))。水流漫過壩頂后，沿最易沖蝕的路徑運動，沖蝕下游坡面，形成許多細(xì)小的初始沖槽，沖槽內(nèi)水流較小，潰口發(fā)展也相對較慢。2)溯源沖蝕(見圖6(c)和6(d))。初始階段的沖槽發(fā)展成較大的沖槽，并向上游發(fā)展至壩頂，沖槽內(nèi)同時存在垂直下切和橫向擴展現(xiàn)象，但下切速度要大于橫向擴展的速度，以下切發(fā)育為主。3)沿程沖蝕(見圖6(e)和6(f))。潰口水流主要作用于沖槽兩側(cè)邊壁，側(cè)壁連續(xù)沖刷和邊坡坍塌交替發(fā)生，坍塌體逐漸被水流帶走，當(dāng)水流無法繼續(xù)沖蝕壩料時潰壩過程停止，3組模型結(jié)果均顯示潰壩結(jié)束時無殘留壩高。綜合試驗?zāi)Ｐ虯和C的結(jié)果可以看出：尾礦庫潰口發(fā)展過程主要包括由水流沖刷引起的潰口快速下切和橫向沖刷以及潰口邊坡失穩(wěn)坍塌引起的間歇性橫向擴展。

圖6 尾礦庫潰壩離心模型試驗過程Fig.6 Centrifugal model test processes of tailings dam breach

1.3 潰口流量變化過程

在潰壩過程的初始階段，壩前水位上升，約10 s 后達到最高水位44.56 cm。由于潰口迅速發(fā)展，潰壩水流流量持續(xù)增加，庫水位迅速下降。潰壩開始30 s左右，水位下降速度減緩，這主要是因為庫水位已退至干灘面尾部，從模型C的試驗結(jié)果可以看出，此階段潰口呈“漏斗狀”。40 s 后潰口發(fā)展到干灘面尾部，潰口流量迅速增加，庫水位又迅速下降。此后，庫水位下降速度隨潰口發(fā)展逐步變緩，過程曲線呈“鋸齒狀”，這主要是由于潰口邊坡大塊滑坡導(dǎo)致潰壩水流的流量降低，外部入流作用導(dǎo)致水位短時間內(nèi)抬升。離心模型試驗潰口流量變化如圖8所示。

圖8 尾礦庫漫頂潰壩潰口流量變化Fig.8 Breach flow changes of tailings dam due to overtopping failure

由圖8可以看出：潰壩初始階段潰口流量快速增大，在潰壩32 s 時達到峰值流量8.75 L/s，此階段潰口快速下切，這主要是因為匯水面積不斷增大。在達到峰值流量之后的30 s內(nèi)，潰口流量迅速下降，上一階段較大的潰口流量使水位迅速下降至干灘面尾部底床附近，庫中沒有足夠的水流迅速補充到潰口處，因此潰口流量迅速下降。

通過尾礦庫潰壩離心模型試驗可以看出，在水流作用下，尾礦庫與均質(zhì)土壩漫頂潰壩過程存在一定的相似之處，但由于尾礦庫結(jié)構(gòu)、材料的獨特性，在潰壩水流作用下形成高濃度水流，潰口下切速度非?？欤茐哪芰Ω鼜?，潰口發(fā)展規(guī)律和流量過程更復(fù)雜。因此，針對尾礦庫漫頂潰壩的數(shù)值模擬要重點考慮水流運動與尾礦料輸移的相互作用。

2 數(shù)學(xué)模型

2.1 數(shù)學(xué)模型控制方程

基于靜水壓力假設(shè)的二維淺水方程可作為描述潰口流量的控制方程[27-29]。但由于尾礦庫潰口下泄物中的固相含量較高，需要考慮挾砂水流的沖蝕對潰口形態(tài)演化的影響。因此，在淺水方程中引入挾砂水流的濃度，進而可考慮水流、壩料和潰口底床變化的耦合作用，控制方程如下：

式中：Φ為守恒變量，Φ=[h,uh,vh,hCt]T；F(Φ)和G(Φ) 為 單 元 界 面 的 法 向 通 量，F(xiàn)(Φ)=[uh,u2h,uvh,huCt]T，G(Φ)=[vh,uvh,v2h,hvCt]T；S(Φ)為源項，S(Φ)=[s1,s2,s3,s4]T；s1，s2，s3和s4可分別表示為

其中，h為潰口水深；u和v分別為x和y方向的平均流速；Ct為挾砂水流質(zhì)量濃度；pm為潰口處壩體孔隙率；U為深度平均流速；mb為潰口底床的坡度系數(shù)，mb=[1+(?zb/?x)2+(?zb/?y)2]1/2；qt*為壩料沖蝕率；L為壩料沖蝕從非平衡到平衡狀態(tài)的適應(yīng)長度；ρw為水的密度；ρs為壩料的密度；g為重力加速度。

2.2 數(shù)學(xué)模型求解

2.2.1 控制方程離散

模型的控制方程采用矩形網(wǎng)格上的Godunov型格式有限體積法計算。式(1)的離散化方程為

由于控制方程中含有考慮流固耦合和河床變化影響的源項，為方便起見，將控制方程分為齊次偏微分方程和常微分方程。

2.2.2 通量計算

本研究通過前瞻性隊列研究，在 435例乙型肝炎病毒感染患者中進一步探討和驗證了 2型糖尿病與乙型肝炎病毒相關(guān)肝癌發(fā)病的關(guān)系，并提供了 2型糖尿病與乙型肝炎病毒相關(guān)肝癌發(fā)病風(fēng)險的前瞻性研究證據(jù)。但本研究也存在一定的局限性，如沒有考慮糖尿病病程、糖尿病用藥、抗病毒治療等因素對研究結(jié)果的影響，因此尚需納入相關(guān)因素進一步開展更大樣本量的研究。總之，在本隊列人群中，2型糖尿病和乙型肝炎病毒相關(guān)肝癌發(fā)病有關(guān)，2型糖尿病增加了乙型肝炎病毒感染患者的肝癌發(fā)病風(fēng)險。

有限體積算法的關(guān)鍵是如何計算通過界面的流體的流量通量和動量通量以及尾礦料流量通量。有限體積法不存在守恒誤差，其主要誤差來自網(wǎng)格單元界面上通量計算誤差。對于齊次偏微分方程，采用HLLC黎曼求解器計算流體質(zhì)量和動量的網(wǎng)格界面通量。

圖9 二維有限體積網(wǎng)格Fig.9 Two dimensional finite volume mesh

尾礦料流量通量不采用黎曼求解器，界面上尾礦料通量為水流通量與尾礦料質(zhì)量濃度的乘積。相鄰2個計算單元界面只有一個水流通量，而界面左右兩側(cè)存在2種尾礦料質(zhì)量濃度，尾礦料的界面通量可以通過中波速SM與0的關(guān)系來確定[31]：

2.2.3 數(shù)值重構(gòu)

為了在空間上達到二階精度，選擇一種分段線性數(shù)值重構(gòu)的方法(MUSCL 格式)提高解的空間精度[32]。為了避免出現(xiàn)大梯度下的假振蕩現(xiàn)象，采用全變分衰減(TVD)格式來描述振蕩[33]，利用坡度限制器Δ?i來控制變量或通量的梯度。

2.2.4 源項離散

由于模型控制方程分解為齊次偏微分方程和常微分方程，所以源項采用直接離散處理。源項離散方式對數(shù)值解的精度有很大的影響，模型控制方程中的源項主要包括水面梯度項、潰口底床摩擦阻力項、濃度梯度項等。

潰口底床坡度在靜水條件下會產(chǎn)生虛假流動現(xiàn)象，故采用改進的表面梯度法(SGM)計算水面梯度項[35]。

為了得到穩(wěn)定的計算結(jié)果，采用半隱式格式計算潰口底床摩擦阻力項[27,36]。

此外，采用中心差分格式處理濃度梯度項，以獲得二階精度[27,37]。

2.2.5 尾礦料沖蝕率

由式(1)可知：影響沖蝕過程的關(guān)鍵輸入?yún)?shù)為L和qt*。L可通過下式[38]確定：

式中：Lb為推移質(zhì)自適應(yīng)長度；α為懸移質(zhì)自適應(yīng)系數(shù)。

尾礦料沖蝕率qt*包括懸移質(zhì)沖蝕率qs*和推移質(zhì)沖蝕率qb*這2 個部分，懸移質(zhì)沖蝕率可通過下式[39]計算：

式中：Φs為量綱一懸移質(zhì)沖蝕率；γs和γ分別為尾礦料和混合物的容重；d50為尾礦料平均粒徑；τ為流體剪應(yīng)力，τ=γRJ；R為水力半徑；J為能坡；τc為尾礦料的臨界起動剪應(yīng)力。

同理，推移質(zhì)沖蝕率可通過下式[39]計算：

式中：Φb為量綱一推移質(zhì)沖蝕率；τb為潰口床面的流 體 剪 應(yīng) 力，τb=ρgN2mbU2/h1/3；N為 曼 寧 糙 率系數(shù)。

當(dāng)壩體坡度較陡時，需要考慮重力分量對尾礦料起動條件和沖蝕率的影響。采用考慮重力影響的方法對τb和τc進行修正[40]。

式中：τbm和τcm分別為潰口床面的流體剪應(yīng)力修正值和尾砂料臨界剪應(yīng)力修正值；φ為尾礦料的休止角；θ為壩體邊坡坡角；cosωz為平面流動方向上河床高程的梯度。

2.2.6 潰口邊坡失穩(wěn)

隨著尾礦庫潰壩的持續(xù)，潰口邊坡可能發(fā)生失穩(wěn)，模型采用矩形網(wǎng)格上二維邊坡滑動的數(shù)值計算方法[31,37]，即當(dāng)相鄰網(wǎng)格中心點連線的坡角大于尾礦料休止角時，邊坡將發(fā)生滑坡，形成潰口底床，其方程表達式如下：

式中：zbi和zbi+1分別為計算點i和i+1的潰口底床高程；Δzbi和Δzbi+1分別為滑坡造成的計算點i和i+1的潰口底床高程變化量；Δli為計算點i和i+1 之間的距離；當(dāng)計算點i至i+1 為上坡方向時，tanφ取正值，當(dāng)計算點i至i+1 為下坡方向時，tanφ取負(fù)值。

若8 個相鄰網(wǎng)格中有1 個以上滿足滑坡條件，則將這些網(wǎng)格分為上坡網(wǎng)格和下坡網(wǎng)格。相鄰網(wǎng)格的床面高程變化可用下式表示：

式中：Nup和Ndown分別為上坡和下坡網(wǎng)格數(shù)。

同時，中心網(wǎng)格與8個相鄰單元上滑坡尾礦料總質(zhì)量應(yīng)滿足質(zhì)量守恒定律：

由滑坡條件和質(zhì)量守恒定律可以得到中心網(wǎng)格底高程變化量的表達式：

2.2.7 模型穩(wěn)定性判定

在計算迭代步驟中，若潰口底床高程變化量Δzb過大，則無法保證尾礦料沖蝕率計算結(jié)果的準(zhǔn)確性[41]；經(jīng)調(diào)試發(fā)現(xiàn)Δzb應(yīng)小于水深的10%。此外，為了保證計算效率，采用動態(tài)時間步長方法，即在每次迭代前，模型通過式(17)為迭代計算選擇最優(yōu)時間步長。

2.2.8 數(shù)值模擬流程

尾礦庫漫頂潰決過程數(shù)值模擬計算流程圖如圖10所示。

圖10 尾礦庫漫頂潰決過程數(shù)值模擬流程圖Fig.10 Numerical simulation flow chart of tailings dam overtopping breach process

3 漫頂潰決過程數(shù)學(xué)模型驗證

為了驗證基于尾礦庫潰壩機理的數(shù)學(xué)模型和數(shù)值計算方法的合理性，首先選擇雙側(cè)擴寬水槽中動床模型[42]進行試驗，該模型試驗常被學(xué)者們作為基準(zhǔn)試驗來驗證流固耦合模型的準(zhǔn)確性[27,43]。隨后，選擇尾礦庫漫頂潰壩離心模型B進行試驗，測試流固耦合模型的合理性。

3.1 雙側(cè)擴寬水槽中動床模型試驗

本文選用文獻[42]中的算例分析動床條件下潰壩水流引起的底床變化。水槽試驗參數(shù)見表5，試驗平面布置示意圖如圖11 所示。試驗通過快速抽起1 m寬的閘門模擬潰壩過程，閘門兩側(cè)為不透水塊，上、下游寬度均為1.0 m。水槽底床預(yù)先鋪設(shè)飽和砂，鋪設(shè)范圍從潰口上游1 m至潰口下游9 m。以水槽定床為參考面，蓄水池初始水位為0.47 m。在潰口下游處預(yù)先布置8 個超聲探頭測量水位變化，測量點U1～U8具體坐標(biāo)值見表6。

表5 水槽試驗參數(shù)Table 5 Parameters for flume experiment

圖11 水槽試驗平面布置示意圖[42]Fig.11 Schematic layout diagram of flume experiment[42]

表6 水位測量點坐標(biāo)Table 6 Coordinates of water level measuring points

各監(jiān)測點水位計算結(jié)果與實測結(jié)果對比如圖12所示，最終水位實測值與計算值對比見表7。由圖12 和表7 可以看出：數(shù)值計算模型較好地預(yù)測了各監(jiān)測點水位變化。監(jiān)測點U2和U3距離潰口較近，水流出現(xiàn)振蕩趨勢，結(jié)合底床變化來看，可能是由底床高程的劇烈變化所致。監(jiān)測點U1和U4在水槽擴寬的位置，在初始階段，水位數(shù)值計算結(jié)果明顯低于實測值，監(jiān)測點U5和U8的水位變化較為平穩(wěn)。

表7 監(jiān)測點最終水位測量值與計算值對比Table 7 Comparison of measured and calculated values of final water level at each monitoring point

圖12 各監(jiān)測點水位實測值與計算值變化過程對比Fig.12 Comparison of measured and calculated water level variations at each monitoring point

數(shù)值計算與實測的最終底床厚度對比如圖13所示。由圖13 可以看出：接近潰口的區(qū)域表現(xiàn)出較為明顯的沖蝕特征，數(shù)值計算結(jié)果顯示底床泥砂被挾砂水流全部帶走，底床泥沙厚度接近于0 m，實測結(jié)果也在0.015 m 以下。下游水槽壁都表現(xiàn)出明顯的沉積特征，當(dāng)潰壩波傳到水槽兩側(cè)，沖擊側(cè)壁形成水躍，水流強度降低，在下游和水側(cè)兩側(cè)形成沉積特征，底床泥沙厚度為0.125～0.145 m，高于原始高度0.085 m，與實測結(jié)果相似?？紤]到潰壩水流輸砂過程的復(fù)雜性，可認(rèn)為本文的模型較好地反映了挾砂水流作用下的沖蝕與沉積特征。

圖13 實測與計算所得最終底床形狀對比Fig.13 Comparison of measured and calculated morphologies of final bottom bed

3.2 尾礦庫漫頂潰壩離心模型試驗

為進一步驗證模型計算的準(zhǔn)確性，選取尾礦庫漫頂潰壩離心模型B進行試驗驗證，模型B尾礦料特性和試驗壩參數(shù)分別見表3 和表4，按照表1中的相似準(zhǔn)則，將模型試驗實測值轉(zhuǎn)化為原型實測值。表8所示為不同參數(shù)實測值與計算值的對比結(jié)果。由表8可以看出不同參數(shù)相對誤差均在±5%以內(nèi)。

表8 尾礦庫漫頂潰壩離心模型B不同參數(shù)實測值與計算值對比Table 8 Comparison of measured and calculated values of different parameters in centrifugal model B for overtopping breach of tailings dam

圖14 所示為尾礦庫潰壩潰口流量實測值與計算值對比，圖15 所示為壩體縱斷面潰口發(fā)展過程計算結(jié)果。由圖14和圖15可以看出：在潰壩開始初期，潰口發(fā)展緩慢，此階段主要為下游坡沖蝕；隨著潰壩的發(fā)展，溯源沖刷發(fā)展至壩頂潰口的位置，潰口在挾砂水流和邊坡滑坡的交替作用下快速發(fā)展，下切速度快于橫向發(fā)展的速度，此時段水位迅速下降；后續(xù)潰口流量逐漸減小，潰口發(fā)展速率減緩。由圖14和圖15可以看出：數(shù)值模擬得到的潰壩過程與離心模型試驗結(jié)果基本吻合，驗證了模型的合理性。

圖14 潰口流量實測值與計算值對比Fig.14 Comparison of measured and calculated breach flow hydrographs

圖15 潰口縱剖面形態(tài)變化計算結(jié)果Fig.15 Calculated results of morphology evolution in the longitudinal section

4 結(jié)論

1)按照潰口發(fā)展規(guī)律和潰壩水動力條件，可將尾礦庫漫頂潰壩過程劃分為下游坡沖蝕、溯源沖蝕、沿程沖蝕3個階段。

2)基于潰壩水流的挾砂能力，提出一個可模擬尾礦庫漫頂潰壩過程中流固耦合現(xiàn)象的數(shù)學(xué)模型，該模型綜合考慮了潰壩過程中的水砂混合物水動力條件、壩料在陡坡上的非平衡沖蝕和潰口邊坡失穩(wěn)坍塌等尾礦庫潰壩特點，并采用基于變時間步長的有限體積法對控制方程進行求解。

3)流固耦合模型具有良好的適應(yīng)性，可準(zhǔn)確捕捉潰口的演化規(guī)律和潰口流量的突變特征，且計算結(jié)果與實測結(jié)果基本吻合，驗證了數(shù)學(xué)模型和數(shù)值計算方法的合理性。