王 琳，李炎隆，武鈺淼

(西安理工大學 西北旱區(qū)生態(tài)水利國家重點實驗室， 陜西 西安 710048)

堰塞湖是一種典型的地質(zhì)災害，具有壩體結構松散、幾何形態(tài)各異、潰決風險高、潰決致災后果嚴重等特點。近年來，我國堰塞湖呈現(xiàn)大規(guī)模、高頻率、群發(fā)性、風險持續(xù)增加的趨勢，是目前堰塞湖潰決數(shù)量最多的國家[1]。堰塞湖一旦潰決，極易形成災害鏈，威脅下游人民群眾的生命財產(chǎn)安全[2]。

據(jù)統(tǒng)計[3]，9%的堰塞湖在1小時內(nèi)潰決，34%在1天內(nèi)潰決，67%在1個月內(nèi)潰決。2000年4月，易貢堰塞湖[4]形成，潰決時峰值流量達到94 810 m3/s。2008年汶川特大地震造成257處堰塞湖，其中最具代表性的是唐家山堰塞湖[5]和小崗劍堰塞湖[6]，潰決洪峰分別達到6 500 m3/s和3 950 m3/s，使得近百萬人的生命財產(chǎn)安全受到威脅。2018年白格堰塞湖更是兩次堵塞金沙江上游干流[7]，潰決洪峰流量分別達到10 000 m3/s和31 000 m3/s[8-9]，導致下游在建的蘇洼龍水電站圍堰拆除，葉巴灘、拉哇、巴塘等在建水電站工程延期。

沖蝕是堰塞湖潰決的主要因素之一，定量分析沖蝕過程是揭示堰塞湖壩體土料沖蝕特性的重要挑戰(zhàn)[2]。土料的沖蝕特性以沖蝕流速(沖蝕臨界剪應力)和沖蝕率表征。水流對土體表面顆粒施加剪應力，當該剪應力超過臨界值時，土料發(fā)生沖蝕，對應的水流平均速度為沖蝕流速。沖蝕流速決定了土體沖蝕開始和終止的時間，是影響沖蝕率的重要參數(shù)。由于沖蝕率與剪應力的關系均呈現(xiàn)不同形式，進而導致潰決洪水分析計算結果極為不穩(wěn)定，有時甚至不收斂。以唐家山堰塞湖為例，作者分別選用Einstein-Brown、Englund-Hensen、Du Boys和Meyer Peter & Muller四種不同的沖蝕模型開展?jié)Q洪水分析，通過與實測數(shù)據(jù)對比[10]，發(fā)現(xiàn)不同沖蝕模型計算出的潰決流量結果差距可達2～3倍。

因此，許多不同的沖蝕特性測量系統(tǒng)被開發(fā)用以研究沖蝕特性。廣泛應用的有原位流槽法[11]、室內(nèi)水槽法[12]、旋轉圓柱測試(Rotation Cylinder Test)[13]、沖蝕裝置(Erosion Function Apparatus，EFA)[14]、孔洞沖蝕測試(Hole Erosion Test)[15]、噴射沖蝕測試(Jet Erosion Test)[16]等。Wahl[17]在對HET和JET設備開展沖蝕試驗對比時發(fā)現(xiàn)，JET設備的沖蝕率比HET高一個數(shù)量級，臨界剪應力卻低兩個或多個數(shù)量級?？梢姴捎貌煌瑴y量系統(tǒng)可得到不同的沖蝕模型，由于采用的系統(tǒng)不同，會產(chǎn)生較大差異。這些沖蝕測量系統(tǒng)沒有標準，無法模擬實際潰決洪水流速，較難測量含有大體積巖石碎塊的堰塞湖土料沖蝕特性。

為模擬堰塞湖潰決過程，本文擬提出一種新的土料沖蝕特性測量方法，并通過自行研制的圓筒型沖蝕試驗設備(CETA)測量土料沖蝕特性。與上述測量系統(tǒng)相比，自行開發(fā)的CETA具有以下優(yōu)點。首先，能提供7 m/s的最高流速，可滿足潰決的快速變化過程。其次，可對粗料(砂、礫石)進行試驗。第三，可針對最大粒徑10 cm的土料開展沖蝕試驗。本文將通過介紹圓筒型沖蝕試驗設備，并選用易貢堰塞湖土料開展不同粒徑的啟動流速和沖蝕率的試驗研究，建立土體沖蝕流速與沖蝕率之間的關系式，檢測土料沖蝕特性測量方法的可靠性；并通過確定剪應力計算公式，建立堰塞湖土料沖蝕率與沖蝕剪應力間的沖蝕模型。

1 圓筒型沖蝕試驗設備(CETA)結構與原理

1.1 結 構

CETA由結構系統(tǒng)、傳動系統(tǒng)、沖蝕試驗系統(tǒng)、清淤循環(huán)系統(tǒng)和控制及數(shù)據(jù)采集系統(tǒng)組成，見圖1。傳動系統(tǒng)包括電機和葉輪(見圖2、3)。在實驗中，電機旋轉并驅動葉輪，葉輪又驅動圓筒中的水流來沖刷土料樣品。調(diào)速電機額定電壓為220 V，最高轉速為1 400 r/min，由沉砂池、水泵和散水器組成的清淤循環(huán)系統(tǒng)可實現(xiàn)水循環(huán)。土顆粒被沖走后，隨水流入沉砂池，沉砂池中的水通過水泵和散水器回到筒中，而土顆粒留在沉砂池中。在圓筒的四周專門設置了3個由有機玻璃制成的觀察窗，可用于觀察和拍攝土體沖蝕過程，每個觀察窗的尺寸約為40 cm×28 cm。另外，考慮到圓筒渦流可能會影響沖蝕過程，所以在筒底部加裝了一個圓柱體的燈環(huán)，使筒底接近一個環(huán)形水槽，降低渦流的影響，提高試驗的觀測和拍攝效果。試樣可采用原狀土或重塑土，也可以是細粒土或粗粒土。此圓筒的內(nèi)徑為104 cm，圍繞中軸的照明燈環(huán)直徑為35.4 cm，燈環(huán)的高度為30 cm。

圖1 圓筒型沖蝕試驗設備Fig.1 Cylindrical rotating erosion apparatus

圖2 電機Fig.2 Motor

1.2 原 理

設備通過電機帶動葉輪轉動，葉輪轉動促使水流發(fā)生流動，水流啟動后即沖刷已放置在圓筒底部的試樣，一旦試樣出現(xiàn)運動，此時的速度即為沖蝕啟動流速。采集系統(tǒng)可記錄流速和試樣沖刷過程。設備的葉輪最大設計轉速為11 m/s，試驗水深范圍為0.4～0.7 m。細粒料樣品，如粘土、粉土和細砂，可放入土壤樣品箱(見圖4(a))，類似于EFA。在試驗過程中，用提升系統(tǒng)將試樣推進5 mm，記錄在一定水流速度下5 mm試樣的沖蝕時間，可用于模擬真實堰塞湖壩體中的粗礫料土樣(見圖4(b))。確定流速下的沖蝕率可通過計算固定時間內(nèi)的土體沖蝕量計算。

圖4 土樣的兩種不同布置方式Fig.4 Two different layout patterns of the soil sample

具體試驗程序可簡化為如下。

1) 土樣制備。根據(jù)要求的密度和含水量制備土樣，將制備好的土樣放置在可升降的樣品箱或圓筒底部。

2) 組裝圓筒和攝像頭。隨后，將水注入圓筒，高度為55 cm。

3) 打開相機并調(diào)整其位置。

4) 打開主電源，啟動設備進行預熱。

5) 在沖蝕開始之前，旋轉速度一直增加，并記錄轉速。

6) 通過采集系統(tǒng)測量并記錄沖蝕流速和沖蝕速率。沖蝕率可由沖刷高度(5 mm)或測定時間間隔內(nèi)的殘余土料換算而成。

7) 根據(jù)試驗要求，提高轉速并重復步驟6。

8) 測試結束后，關閉電源并清潔儀器。

2 沖蝕剪應力計算方法

2.1 葉輪轉速與水流速度關系的標定

圓筒型沖蝕試驗設備工作時，是電機帶動葉輪轉動，再通過水的介質(zhì)作用，將速度傳遞到底部的試樣上。通過水的傳遞作用，葉輪中心位置處的流速值與底部試樣位置上的流速值存在差異，需要率定?？紤]到圓筒設備是新裝置以及它的特殊性，選用流速儀對試驗設備進行率定。采用FP111型直讀式流速儀(見圖5)，通過測量沖蝕設備過水斷面的不同位置在不同葉輪轉速下的流速來對設備內(nèi)部的流場進行率定，建立不同葉輪轉速下設備內(nèi)部流速場的分布特征，確定實測流速與葉輪轉速的相關關系，圖6為流速測定的測點位置示意圖。

圖5 FP111直讀式流速儀Fig.5 FP111 reading type current meter

圖6 測點分布圖 Fig.6 Measuring point

在不一樣的測點保持共同的轉速，在每個測點記錄的時間為40 s左右，并記錄其在此時間內(nèi)的平均流速和最大流速。結果見表1、表2和表3，其中水深均為55 cm。

表1 1、4、7、10四個點的率定結果Tab.1 Calibration results of four points 1, 4, 7, 10

表2 2、5、8、11四個點的率定結果Tab.2 Calibration results of four points 2, 5, 8, 11

表3 3、6、9、12四個點的率定結果Tab.3 Calibration results of four points 3, 6, 9, 12

在率定過程中，發(fā)現(xiàn)無論水量多少，測點2(見圖6)處的試樣均為首先啟動，驗證了圓筒內(nèi)部流速場的不均勻性。水流沿筒內(nèi)葉輪的方向轉動，水速分布主要由相對于葉輪的位置來決定，又因為2點的試樣最先啟動，故取各轉速對應的2，5，8，11四個測點(見圖6)最大流速值的平均值，得到轉速與水速之間的換算關系。其換算關系式為：

y=0.679x-0.0103

(1)

式中：y為底部中心流速值(m·s-1)；x為葉輪中心轉速值(m·s-1)。

2.2 沖蝕剪應力計算方法

文獻[18]認為在無法直接測量沖蝕率和剪應力的情況下，最有效的替代方法是使用Moody圖(見圖7)估計其剪應力。

圖7 Moody圖Fig.7 Moody chart

本研究對圓筒型沖蝕試驗設備采用Briaud[14]的剪應力計算公式：

(2)

式中：τ為剪應力(N·m-2)；ρw為水的密度(1 000 kg·m-3)；V為平均速度；f為摩擦系數(shù)。

Reynolds數(shù)Re計算公式為：

(3)

式中：D為管道的水力直徑；v為水在25℃時的運動粘度(10-6m2·s-1)。水力直徑D為水力半徑R的4倍，水力半徑由流動面積除以潤濕周長來定義，此公式表示為：

D=2ab/(a+b)

(4)

則水力直徑為：

(5)

當Re>100 000時，直接在Moody圖中尋找對應的摩擦系數(shù)f。

當Re<100 000時，用：

(6)

此方程為Blasius方程，式(6)可計算摩擦系數(shù)，不必使用Moody圖查詢。

3 試驗結果討論

見圖8，選用易貢的兩種典型土料開展沖蝕試驗，D50=8 mm，D50=10 mm。樣品的面積控制為1 075 cm2和1 085 cm2，結果見圖9。試驗中土體天然密度分別為1.845 g/cm3和1.585 g/cm3，水深為55 cm。沖蝕率定義為單位時間內(nèi)土料樣品被沖走的高度，土料質(zhì)量控制在4.00 kg和3.72 kg。

圖8 易貢土料典型級配曲線Fig.8 Gain size distributions of two Yigong landslide dam materials

圖9 典型粒徑級配沖蝕率Fig.9 Erosion rate of particle size gradation

3.1 沖蝕率

Briaud[14]運用EFA設備對干凈粗砂開展了沖蝕率研究。此土體的D50=3.375 mm，其沖蝕率最大為12 000 mm/h=0.33 cm/s，見圖10，其沖蝕流速與沖蝕率呈雙曲線函數(shù)關系。采用圓筒型沖蝕試驗設備針對易貢的典型土料(D50=8 mm，D50=10 mm)開展沖蝕試驗，試驗結果見表4和表5。分別對以上數(shù)據(jù)運用R語言進行線性回歸，建立沖蝕流速和沖蝕率之間的關系式，見式(7)和式(8)。

表4 粒徑級配D50=8 mm沖蝕率Tab.4 Erosion rate of particle size gradation D50=8 mm

表5 粒徑級配D50=10 mm沖蝕率Tab.5 Erosion rate of particle size gradation D50=10 mm

圖10 砂土的沖蝕曲線Fig.10 Erosion curve for coarse sand

沖蝕流速與沖蝕率的關系，經(jīng)回歸后當D50=8 mm時：

(7)

當D50=10 mm時：

(8)

式中：v為沖蝕流速(m/s),y為沖蝕率(m/s)。

由圖10可知，沖蝕流速與沖蝕率之間呈雙曲線函數(shù)關系，試驗中得到?jīng)_蝕率的最大值為0.17 cm/s，試驗結果中的量級和圖10中關于干凈粗砂的量級一致，數(shù)值略小于文獻中0.33 cm/s的數(shù)值。之所以出現(xiàn)此情況是因為EFA試驗選用的材料為干凈粗砂，而易貢主要堆積區(qū)材料除了塊石、砂之外還包含一些粉塵及碎屑，所以試驗結果稍小于Briaud[14]的數(shù)值。

Briaud[14]提到，EFA設備中的相同平均速度為1 m/s，砂的沖蝕率約是粘土的1 000倍，表明不同土體的沖蝕率可能存在較大差異。文獻中提到干凈的砂和礫石的沖蝕率數(shù)值量級為104mm/h，與本文試驗結果一致。

綜上所述，試驗結果顯示沖蝕流速與沖蝕率呈雙曲線函數(shù)關系，且有很大的相關性，通過結果分析認為圓筒型沖蝕試驗設備的可靠性值得信賴，也驗證了沖蝕率試驗結果的可靠性，可用來開展沖蝕剪應力研究。

3.2 沖蝕剪應力

參照EFA計算公式、圓柱的剖面為長方形，計算濕周，參考公式(2)～(6)計算剪應力，并采用回歸方程尋找沖蝕流速、沖蝕率和沖蝕剪應力τ之間的關系。其沖蝕流速、沖蝕率和沖蝕剪應力之間關系見表6和表7。

表6 D50=8 mm時沖蝕速度與沖蝕剪應力的關系Tab.6 Relationship between velocity and critical shear stress in D50=8 mm

表7 D50=10 mm時沖蝕速度與沖蝕剪應力的關系Tab.7 Relationship between velocity and critical shear stress in D50=10 mm

沖蝕流速與沖蝕剪應力的關系，經(jīng)回歸后當D50=8 mm時指數(shù)曲線為如下：

(9)

當D50=10 mm時指數(shù)曲線為如下：

(10)

式中：τ為沖蝕剪應力(N·m-2)；y為速度(m·s-1)。

3.3 沖蝕率剪應力模型

本課題組[19]基于唐家山實測的沖蝕率，認為土體材料抵抗沖蝕時，不應有無限“強度”，建議采用雙曲線模型，其形式如下：

(11)

式中：v為扣除臨界剪應力后的剪應力：

v=k(τ-τc)

(12)

雙曲線模型的沖蝕率是具有物理意義的，認為土體包含一定的強度，當沖蝕達到一定程度后，將不再增加(見圖11)。即雙曲線有一當v接近無限值時的漸進線，即dz/dt的極值1/b。其中單位變換因子k為100，其含義為v等于0時曲線的斜率。

圖11 堰塞湖土料沖蝕剪應力與沖蝕率的關系圖Fig.11 Relationship between erosion shear stress and erosion rate of dam-lake soil

在所有潰決洪水分析方法中，研究人員多選用繼承了泥沙-水力學領域中開發(fā)的土體沖蝕模型，如Engelund-Hansen、Meyer-Peter-Mueller及Einstein-Brown，分別應用于計算機程序MIKE11模型、BREACH模型和BEED模型。這些模型多用來評價低流速河床上的泥沙輸運模式，高速潰決洪水過程中的土體沖蝕行為很少研究。由于流速的定義關系，沖蝕剪應力與沖蝕率有很大的相關性，較快的流速產(chǎn)生較大的沖蝕剪應力。上述大多數(shù)用于測量粘性土體和巖石的沖蝕率分析表達式包含大量參數(shù)。Trammell[20]認為基于這些變量參數(shù)準確的預測其相關性是不切實際的，而且評估大量參數(shù)導致成本過高。

Einstein[21]，Partheniades[22]，van Prooijen 和 Winterwerp[23]都認為沖蝕率可由與沖蝕剪應力相關的函數(shù)關系式表示。Slagle[24]還認為可運用類似EFA-SRICOS或Miller-Sheppard模型[24-25]方法，使用純經(jīng)驗公式模擬沖蝕率剪應力關系。因此，純粹的經(jīng)驗方法可用來擬合沖蝕率-沖蝕剪應力關系。van Prooijen和Winterwerp認為在低于原定義的臨界剪應力時仍可能發(fā)生沖蝕；在大于臨界應力時，該線性近似仍是有效的。因此，通過首先計算沖蝕率，再繪制其與剪應力的函數(shù)關系圖，建立沖蝕率與沖蝕剪應力的沖蝕模型。

Briaud[14]提出，在干凈的砂和砂礫中，影響τc的主要土體參數(shù)是由D50引起的。本文改變D50的大小，研究沖蝕率和沖蝕剪應力的關系。沖蝕率與沖蝕剪應力數(shù)值見表3和表4。通過對沖蝕率和沖蝕剪應力的關系回歸后公式見下。

D50=8 mm：

(13)

D50=10 mm：

(14)

式中：τ為剪應力(N·m-2)；y為速度(m·s-1)。

本課題組[19]基于唐家山實測的沖蝕率，提出了沖蝕率-剪應力雙曲線模型，認為土體材料抵抗沖蝕時，不應有無限“強度”。上述研究證實沖蝕率與剪應力之間呈雙曲線函數(shù)關系。當剪應力較小時，沖蝕率與剪應力存在較強的雙曲線函數(shù)關系。

4 結 論

土料沖蝕特性可以用水流引起的沖蝕剪應力和土料沖蝕函數(shù)(即沖蝕速率與剪應力曲線)表征。本文介紹了一種新型土料沖蝕測量系統(tǒng)，該測量系統(tǒng)包括一個圓筒型沖蝕試驗裝置(CETA)和一種計算剪切應力的方法。所研制的CETA具有空間小、耗水量少、流速高等優(yōu)點，并能對D50小于30 mm的土體顆粒進行測量。運用此系統(tǒng)開展了易貢土料的啟動流速、沖蝕率、沖蝕剪應力研究。

試驗結果表明：該裝置的沖蝕率結果可靠，與原有文獻的沖蝕率數(shù)值處于同一量級，與Briaud的實驗結果一致，可用來進行堰塞湖土料沖蝕率-剪應力分析；提出了堰塞湖土料沖蝕剪應力的精確計算方法，發(fā)現(xiàn)沖蝕率與剪應力呈現(xiàn)雙曲線關系，建立了堰塞湖土料雙曲線沖蝕模型，可應用于堰塞湖潰決分析模擬。