鐘啟明,陳生水,趙聯(lián)楨,任 強,曹 偉

(1.南京水利科學(xué)研究院巖土工程研究所,江蘇南京 210029;

2.南京水利科學(xué)研究院水文水資源與水利工程國家重點實驗室,江蘇南京 210029;

3.南京林業(yè)大學(xué)土木工程學(xué)院,江蘇南京 210037)

地震、降雨等自然災(zāi)害往往引發(fā)大規(guī)?；潞蜕襟w崩塌,使受災(zāi)區(qū)域河流阻塞,形成天然堰塞壩。堰塞壩的工作條件、壩體幾何特征以及壩體物質(zhì)組成和內(nèi)部結(jié)構(gòu)都與人工土石壩存在明顯差別。首先在工作條件上,因沒有溢洪道和其他泄洪設(shè)施,上游持續(xù)來水使壩前水位快速上漲,堰塞壩很容易發(fā)生漫頂潰決;其次從壩體形態(tài)來看,堰塞壩壩頂一般凹凸不平,存在天然凹槽,且體型多呈不規(guī)則形狀,沿河流方向堆積范圍明顯超過狹窄河谷的寬度;再者從壩體結(jié)構(gòu)來講,因壩體材料沒有經(jīng)過人工篩選和施工碾壓,大部分堰塞壩結(jié)構(gòu)疏松,不均勻性強,壩體材料多為礫石、粗顆粒,粗細混雜,粒徑級配范圍變化大[1-2]。調(diào)查表明[3],90%的堰塞壩都在其形成1a內(nèi)潰決。一旦發(fā)生潰決,所產(chǎn)生的洪水或泥石流將對下游產(chǎn)生嚴重危害。因此,有必要研究堰塞壩漫頂潰決機理及潰決過程,為正確評估堰塞壩潰決致災(zāi)后果、科學(xué)制定堰塞壩潰決防洪應(yīng)急預(yù)案提供技術(shù)支撐。

1 堰塞壩潰決機理和潰決過程分析

目前關(guān)于堰塞壩潰決機理和潰決過程的研究工作主要參考人工土石壩的成果,對于如何合理模擬寬級配堰塞壩材料在高強度水流侵蝕條件下的非平衡沖刷特性的研究工作很少涉及。王光謙等[4]提出在堰塞壩潰壩模型中考慮潰口水力沖刷和重力坍塌擴展,但是該模型使用的沖蝕公式是在恒定均勻流條件下建立的一般平衡沖蝕公式;黃金池[5]利用水槽試驗結(jié)果,將高強度水流非平衡沖刷公式引入到潰口擴展模型,并將潰口垂直下切、橫向擴展、壩坡溯源沖刷3種過程結(jié)合起來建立堰塞壩潰壩數(shù)學(xué)模型,但是該模型的沖蝕公式主要基于室內(nèi)平面水槽試驗,泥沙粒徑范圍為0.18～13.00mm,無法合理反映堰塞壩邊坡上礫石、塊石等大粒徑滑坡堆積物的起動問題,且沖蝕公式中的參數(shù)不易獲取。

要正確模擬堰塞壩漫頂潰決過程,需要從堰塞壩材料粒徑范圍寬的特點入手,同時考慮壩體密實度、壩體材料強度、下游坡度、漫頂水流流速的影響。筆者基于上述考慮,根據(jù)堰塞壩壩體材料的寬級配特性,建立描述堰塞壩漫頂潰決潰口發(fā)展規(guī)律與流量過程的數(shù)值模型,并利用該模型對唐家山堰塞壩的潰決過程進行數(shù)值模擬,驗證了該模型和計算方法的合理性。

2 堰塞壩漫頂潰決過程數(shù)值模型

2.1 寬級配壩體顆粒臨界起動流速

堰塞壩漫頂潰決過程中水流的沖蝕過程以及非均勻巖土顆粒的起動具有非恒定的特性。假設(shè)水流作用下忽略顆粒間的黏性特征,首先根據(jù)堰塞壩壩體土石料顆粒的受力情況求得顆粒起動的臨界流速(圖1)。

圖1 土體顆粒在壩坡上的受力示意圖Fig.1 Forces acting on a soil particle along dam slope

對土體代表顆粒1而言,考慮水流作用下巖土體顆粒所受的力一般有重力W、水流拖曳力Fd[6]、上舉力Fl[6],其表達式分別為

式中：ρs——土顆粒的密度;ρw——水的密度;d50——土顆粒的平均粒徑;g——重力加速度;Cd——拖曳力系數(shù),一般取0.4[7];v——水流流速;Cl——上舉力系數(shù),一般取0.1[7]。

需要指出的是,堰塞壩體土石料級配寬泛,最大粒徑和最小粒徑相差很大。為此,引入與水流方向垂直的附加作用力R(圖1)來考慮粗顆粒對細顆粒的阻攔、遮蔽作用以及細顆粒對粗顆粒的包圍、填實作用。由于目前對R的研究很少,尚無精確理論公式,可近似假定R與顆粒的平均剪力成比例[8-11],即

其中

式中：R——附加作用力;φ——比例系數(shù)與顆粒面積系數(shù)的乘積;τs——不均勻顆粒的平均剪力;Km——無因次系數(shù);M——緊密系數(shù),代表顆粒組成的密實程度,反映顆粒組成的不均勻程度和細顆粒填充孔隙的程度[9];Cu——不均勻系數(shù)。

由式(4)可得

其中

式中K為系數(shù),可以根據(jù)不均勻顆粒起動流速試驗資料及天然河流實測資料確定,在0.785～1.727范圍內(nèi)變化,此處取1.3[9]。

如圖1所示,土顆粒起動時受到的摩擦力可表示為

式中：Ff——土顆粒受到的摩擦力;φ——土顆粒間的內(nèi)摩擦角;θ——壩坡坡角;c——土體的黏聚力。通過受力分析可知,土顆粒1起動的臨界條件為

將式(1)(2)(3)(5)(6)代入式(7)可得土顆粒在壩坡上的臨界起動流速：

2.2 臨界起動流速公式驗證

唐家山堰塞壩潰決過程中,現(xiàn)場測量得到壩體顆粒的起動流速為2.4m/s[12],因此,可利用該起動流速來驗證本文提出的臨界起動流速公式的合理性。根據(jù)相關(guān)資料[13],計算參數(shù)取值為：d50=0.03m,c=60kPa,M=0.75,θ=23°,φ=45°,g=9.8m/s2,ρs=2.6×103kg/m3。將上述參數(shù)代入式(8),可計算得出臨界起動流速vc=2.12m/s(如不考慮附加力作用,臨界起動流速為1.35m/s),與實測結(jié)果基本一致,驗證了寬級配顆粒起動流速的合理性。

2.3 壩頂溢流沖刷速率

當土體顆粒在水流作用下起動后,在潰壩水流的作用下,壩頂潰口和下游壩坡發(fā)生沖蝕,由于壩體材料粒徑級配范圍變化大,因此在分析不同土體陡水槽沖蝕試驗結(jié)果[14-16]的基礎(chǔ)上,選擇d90與d30為代表粒徑,筆者建議了一個計算堰塞壩土石料沖蝕率的經(jīng)驗公式：

其中

式中：Qs——水流沖蝕率;d90——小于某粒徑土石料質(zhì)量分數(shù)為90%所對應(yīng)的顆粒粒徑;d30——小于某粒徑土石料質(zhì)量分數(shù)為30%所對應(yīng)的顆粒粒徑;B——潰口寬度;v*——摩阻流速;vb——潰口底流速;ˉv——水流平均流速;H——堰塞湖水位高程;Hc——潰口底部高程;J——水力梯度;N——潰口處糙率;Qb——潰口流量。

當漫頂潰壩發(fā)生后,水流沿著初始潰口沖蝕下游壩坡,潰口流量可采用以下寬頂堰公式計算：

式中m為流量系數(shù),此處可取0.5[17]。

2.4 沖蝕公式驗證

建議的壩體材料沖蝕公式可反映筑壩材料顆粒級配、密實度、強度、滲透通道傾角和壩坡、摩阻流速及水流速度對沖蝕量的影響,公式中各物理量意義明確且易于測定。為了驗證所建議公式的合理性,運用該公式對4種不同土體的單寬沖蝕量進行計算分析,計算結(jié)果與陡水槽沖蝕試驗結(jié)果[14-16]的對比見圖2。從圖2可以看出,計算結(jié)果位于試驗結(jié)果的變化區(qū)域內(nèi),這表明從工程應(yīng)用角度來看,所建議的壩體材料沖蝕公式是可以接受的。

圖2 水流流速與單寬沖蝕率的關(guān)系Fig.2 Relationship between flow velocity and erosion rate

2.5 潰口發(fā)展

土石壩漫頂破壞潰口發(fā)展狀況見圖3。壩體初始潰口形狀并非一直保持不變,而是隨漫頂水流流速的增大和壩體沖蝕的加劇不斷加深加寬,時間段增量Δti內(nèi)水流下切深度增量為

時間段Δt內(nèi)水流下切深度增量為

式中j為Δt時段內(nèi)的時間步長數(shù)。

如果不考慮潰口邊坡失穩(wěn)和坍塌引起的潰口橫向擴展,潰口底部的沖蝕速率應(yīng)和潰口邊坡的沖蝕速率大體相等,因此假定潰口的深度和寬度以同樣的速率發(fā)展[18],則水流對壩體潰口兩側(cè)的直接沖刷形成的潰口寬度增量ΔB可表達為

時間段Δt內(nèi)水庫水位變化量為

式中：Qin——入庫流量;Sa——庫水位為H時的水庫面積。

潰口受到水流的連續(xù)沖蝕發(fā)生垂向下切和橫向擴展,邊坡也隨水流沖蝕變得越來越陡,當垂向下切深度達到臨界深度時,潰口邊坡發(fā)生間歇性失穩(wěn)坍塌(圖3)。臨界深度可采用極限平衡方法導(dǎo)出[18]：

式中βk為潰口邊坡臨界坡角。

3 數(shù)值計算方法

可通過迭代計算來對上述潰壩數(shù)學(xué)模型進行求解,以模擬堰塞壩漫頂破壞潰口發(fā)展過程,得出潰口流量過程線。為證明所建議模型的有效性和合理性,針對堰塞壩漫頂引起的潰壩問題,進行潰壩洪水過程計算,計算實例為唐家山堰塞壩,具體求解過程見圖4,圖中Hbm為潰口最終底高程,Bu為潰口最終頂寬,Bm為潰口最終底寬。

圖3 土石壩漫頂破壞潰口發(fā)展示意圖Fig.3 Breach development of earth-rocKdam due to overtopping failure

4 模型驗證

“5.12”汶川大地震堵江形成的堰塞壩中以唐家山堰塞壩規(guī)模最大,對下游威脅也最大。壩體位于北川縣城上游4.6km涪江支流通口河峽谷中,其堰塞壩體積約為2037萬m3,按我國堰塞湖壩體等級劃分,唐家山堰塞壩屬于極高危險、潰決損失極嚴重的Ⅰ級風(fēng)險堰塞壩[19]。

堰塞體的基本形狀為寬頂堰,上游坡度18°～22°,下游坡度35°～40°。順河向堰頂寬度為150m左右,堰塞體除表部少量為碎石土外,主體為基巖滑坡形成的巨石、孤塊石,其中下部巖性為弱風(fēng)化～微風(fēng)化,仍保持原地層層序關(guān)系,破裂巖體結(jié)構(gòu)密實,強度高。壩體最大壩高124.4m,順河向底長803.4m,頂長300.0m左右,寬611.8m,體積約為2.037×107m3,壩底高程為669.5 m,上游坡度約為20°,下游平均坡比為1∶1.28[19]。經(jīng)過應(yīng)急搶險,截至2008年6月1日,堰塞壩頂開挖出1條梯形泄流槽,斷面形狀為梯形,兩側(cè)邊坡為1∶1.5。泄流渠堰頂高程為740.4m,底寬8 m,深13m,總長695m,上游平緩段縱坡為0.6%,下游陡坡段縱坡分別為24%和16%[19]。

由于泄流槽經(jīng)2次爆破,6月10日才開始泄流,所以選取6月10日0：00為起始計算時間。根據(jù)上述基本條件,計算時間步長取Δt=0.02h,計算時間設(shè)定為20 h。其他計算參數(shù)如下：堰塞體高度為124.4m,堰塞體橫河向頂寬為611.8m,堰塞體順河向長度為300.0m,下游平均坡比為1∶1.28,初始潰口寬度為8.0m,初始潰口深度為13m,上游入庫流量為0m3/s,堰塞體平均粒徑為0.03m,粒徑比d90/d30為30,堰塞體孔隙率為0.4,堰塞體材料平均密度為2.6×103kg/m3,內(nèi)摩擦角為45°,堰塞體下游平均坡角為23°,堰塞體黏聚力為60kPa,緊密系數(shù)為0.75,流量系數(shù)為0.5。

計算得出唐家山堰塞壩潰壩洪峰流量出現(xiàn)在壩體潰決后11.8h,洪峰流量為6137.4m3/s(實測最大值為6500.0m3/s[19]),見圖5。圖6為潰口發(fā)展過程,潰口頂寬202.0m,底寬120.7 m,深度22.0m(實測潰口呈上寬下窄的“倒梯形”,其開口寬145.0～225.0m,底寬100.0～145.0m,坡高10.0～60.0m[19])。由圖5～6可知,計算結(jié)果與實測資料接近,表明筆者所建議的堰塞壩潰壩數(shù)值模型和數(shù)值計算方法是合理的。

圖4 堰塞壩漫頂潰決數(shù)值模擬流程Fig.4 Flow chart of numerical simulation of barrier dam breach due to overtopping failure

圖5 潰口流量過程對比Fig.5 Comparison of calculated and observed discharge hydrographs at breach

圖6 潰口發(fā)展過程Fig.6 Breach development

5 結(jié) 語

根據(jù)堰塞壩壩體土石料的寬級配特性,建立描述堰塞壩漫頂潰決過程中潰口發(fā)展規(guī)律與流量過程的數(shù)值模型。該模型引入與水流方向垂直的附加作用力來考慮粗顆粒對細顆粒的阻攔、遮蔽作用以及細顆粒對粗顆粒的包圍、填實作用。基于不同土體陡水槽沖蝕試驗結(jié)果,針對寬級配堰塞壩材料在高強度水流侵蝕條件下的沖蝕特性,提出了一個能反映壩體材料顆粒級配、密實度、強度、壩體坡度、水流速度及摩阻流速對沖蝕量影響的土體沖蝕公式,公式中各物理量意義明確且易于測定。利用該模型對唐家山堰塞壩的潰壩過程進行模擬計算,得出的潰口發(fā)展規(guī)律與潰壩洪水流量過程與實測結(jié)果接近,驗證了該模型和計算方法的合理性,從而為正確評估堰塞壩潰決致災(zāi)后果、科學(xué)制定堰塞壩潰決防洪應(yīng)急預(yù)案、減輕因堰塞壩潰決造成的損失提供了一種新的有效手段。

[1]任強,陳生水,鐘啟明,等.堰塞壩的形成機理與潰決風(fēng)險[J].水利水電科技進展,2011,31(5)：30-34.(REN Qiang,CHEN Shengshui,ZHONG Qiming,et al.Formation mechanism and breaching failure risKof barrier dams[J].Advances in Science and Technology of Water Resources,2011,31(5)：30-34.(in Chinese))

[2]鄧明楓,陳寧生,胡桂勝,等.松散及弱固結(jié)堰塞體潰壩形式與流量過程[J].水利水電科技進展,2011,31(1)：11-14.(DENGMingfeng,CHEN Ningsheng,HUGuisheng,et al.Outburst patterns and discharge process for weakly consolidated and loose barrier dams[J].Advances in Science and Technology of Water Resources,2011,31(1)：11-14.(in Chinese))

[3]COATA JE,SCHUSTER RL.The formation and failure of natural dams[J].Geological Society of America Bulletin,1988,100(1)：1054-1068.

[4]王光謙,鐘德鈺,張紅武,等.汶川地震唐家山堰塞湖泄流過程的數(shù)值模擬[J].科學(xué)通報,2008,53(24)：3127-3133.(WANG Guangqian,ZHONG Deyu,ZHANGHongwu,et al.Numerical simulation of dischargeprocess for Tangjiashanbarrier lake due to Wenchuan earthquake[J].Chinese Science Bulletin,2008,53(24)：3127-3133.(in Chinese))

[5]黃金池.堰塞壩漫頂潰口流量變化過程的數(shù)值模擬[J].水利學(xué)報,2008,39(10)：1235-1240.(HUANG Jinchi.Numerical modeling of flow through breach of landslide dams[J].Journal of Hydraulic Engineering,2008,39(10)：1235-1240.(in Chinese))

[6]錢寧,萬兆惠.泥沙運動力學(xué)[M].北京：科學(xué)出版社,2003.

[7]韓其為,何明民.泥沙起動規(guī)律及起動流速[M].北京：科學(xué)出版社,1999.

[8]李榮,王迎春.非均勻沙起動規(guī)律研究[J].泥沙研究,1999(2)：27-32.(LI Rong,WANG Yingchun.Study on laws of threshold motion for non-uniform sediment[J].Journal of Sediment Research,1999(2)：27-32.(in Chinese))

[9]秦榮昱,王崇浩.河流推移質(zhì)運動理論及應(yīng)用[M].北京：中國鐵道出版社,1996.

[10]王協(xié)康,敖汝莊.泥沙起動條件及機理的非線性研究[J].長江科學(xué)院院報,1999,16(4)：39-41.(WANG Xiekang,AO Ruzhuang.Nonlinear study on sediment incipient motion condition and its mechanism[J].Journal of Yangtze River Scientific Research Institute,1999,16(4)：39-41.(in Chinese))

[11]楊具瑞,方鐸,何文社,等.非均勻沙起動規(guī)律研究[J].水利學(xué)報,2002,33(10)：82-86.(YANG Jurui,FANG Duo,HE Wenshe,et al.Study on laws of incipientmotionfor non-uniform sediment[J].Journal of Hydraulic Engineering,2002,33(10)：82-86.(in Chinese))

[12]LIU Ning,CHEN Zuyu,ZHANG Jianxin,et al.Draining the Tangjiashan barrier lake[J].Journal of Hydraulic Engineering,ASCE,2010,136(11)：914-923.

[13]胡卸文,黃潤秋,施裕兵,等.唐家山滑坡堵江機制及堰塞壩潰壩模式分析[J].巖石力學(xué)與工程學(xué)報,2009,28(1)：181-189.(HUXiewen,HUANG Runqiu,SHI Yubing,et al.Analysis of blocking river mechanism of Tangjiashan landslide an dambreaking mode of its barrier dam[J].Chinese Journal of RocKMechanics and Engineering,2009,28(1)：181-189.(in Chinese))

[14]MEYER-PETER E,MULLER R.Formulas for bed load transport[C]//Proceedings of Second Meeting International Association for Hydraulic Research.Stockholm：IAHR,1948：39-64.

[15]SMARTG M.Sediment transport formula for steep channels[J].Journal of Hydraulic Engineering,ASCE,1984,110(3)：267-276.

[16]范家驊,陳裕泰,金德春,等.懸移質(zhì)挾沙能力水槽試驗研究[J].水利水運工程學(xué)報,2011(1)：1-16.(FAN Jiahua,CHEN Yutai,JIN Dechun,et al.Experimental studies on carrying capacity of suspended load[J].Hydro-Science and Engineering,2011(1)：1-16.(in Chinese))

[17]李家星,趙振興.水力學(xué)[M].南京：河海大學(xué)出版社,2001.

[18]SINGH V P.Dam breach modeling technology[M].Dordrecht,Netherland：Kluwer Academic Publisher,1996.

[19]劉寧,張建新,林偉,等.汶川地震唐家山堰塞引流除險工程及潰壩洪水演進過程[J].中國科學(xué)E輯：技術(shù)科學(xué),2009,39(8)：1359-1366.(LIU Ning,ZHANG Jianxin,LIN Wei,et al.Drainage reinforcement project and dam breaKflood process of Tangjiashanbarrier due to Wenchuan earthquake[J].Science in China E：Technological Sciences,2009,39(8)：1359-1366.(in Chinese))