李相南,陳祖煜

(中國水利水電科學研究院巖土工程研究所,北京 100048)

兩種潰壩模型在唐家山堰塞湖潰決模擬中的對比

李相南,陳祖煜

(中國水利水電科學研究院巖土工程研究所,北京 100048)

基于唐家山堰塞壩潰決的實測數據,使用BREACH模型及中國水科院DB-IWHR模型,對唐家山堰塞壩潰決過程進行了反演分析,并對兩種模型進行參數敏感性分析。研究結果表明：使用針對唐家山的參數,兩種模型均可較好地反演唐家山潰決洪水的過程，BREACH模型的下游坡比和孔隙率是對結果影響比較敏感的輸入參數；DB-IWHR模型在沖刷參數和下游水深計算兩方面做了改善,較好地解決了這一問題；DB-IWHR提供了一個Excel格式的計算軟件,數值分析穩(wěn)定性好,使用簡潔客觀,可供類似堰塞湖應急處置時參考。

堰塞湖;BREACH模型;DB-IWHR模型;潰壩洪水分析;唐家山

近年來,由地震、泥石流等自然災害形成的天然堰塞壩在庫水位滲透作用或漫頂沖刷作用下,極易發(fā)生失控性潰壩,危害極大。較為準確地預測潰壩洪水流量,是做好應急預案工作的基本條件。但是潰壩洪水分析的理論和相應的程序還遠不能滿足安全防控的需要。在2008年汶川地震中,模擬預測唐家山堰塞湖潰決洪峰流量時,某一模型預測的最大峰值流量為46 000 m3/s,而實測峰值流量僅有6 500 m3/s[1]。

目前對其模型研究主要包括統(tǒng)計模型、參數化模型和基于物理機理的模型[2]。在實際工程應用中,通常采用基于物理機制的數學參數模型,通過考慮潰口的水流過程、潰口的沖刷侵蝕過程以及潰口的擴展過程,結合三者之間的相互影響建立了如BREACH[3]、BEED[4]、MIKE11[5]、HEC-RAS[6]等堰塞壩潰決計算程序。

陳祖煜等[7]提出了基于物理機制的數學參數模型,在Microsoft Excel中開發(fā)了DB-IWHR模型,提出潰口侵蝕率的雙曲線沖刷模型,引入了邊坡穩(wěn)定計算模擬潰口展寬。該模型采用速度增量形式直接求解控制方程,避免了數值迭代,其成果與唐家山的洪峰流量和洪峰時間吻合,且參數敏感性問題得到了較好的解決,相應的成果也在紅石巖堰塞湖的應急處置中進行了應用[8]。由Fread于1988年提出的BREACH模型是應用較為廣泛的潰壩洪水分析程序,采用Smart修正的 Meyer-Peter & Mueller公式計算沖刷過程,曾成功運用于提堂壩(Teton Dam)、草坪壩(Lawn Dam)、曼塔羅河堰塞壩(Mantaro Landslide Dam)、靈湖堰塞壩(Spirit Lake Blockage)的模擬和預測中[3]。

鑒于潰壩洪水分析包含較大的不確定性,在實際工作中,工程技術人員偏向于采用多個模型和相應的程序,對計算結果進行綜合分析,在此基礎上進行工程決策。本文針對2008年唐家山堰塞湖潰決過程,通過對以上2個模型和相應程序對比分析,考察不同潰壩模型對洪峰流量和洪峰時間的反演結果,并通過敏感性分析,對各個模型的穩(wěn)定性和適用性進行了研究,可供類似堰塞湖應急處置時參考。

1 潰壩模型對比

BREACH和DB-IWHR模型均采用了寬頂堰公式作為潰壩洪水計算的基礎方程,在不同的土石壩或堰塞壩的潰決洪水分析中,兩個模型主要有以下幾個共同點和差異點。

1.1 水力學計算

對于潰口流量計算,BREACH模型和DB-IWHR均采用寬頂堰(圖1)流量公式:

Q=CB(H-z)3/2

(1)

式中:Q為流入潰口通道的流量;B為引流槽進口斷面的寬度;H為庫水位;z為引流槽進口處河床高程;C為流量系數,陳祖煜等[7]對唐家山實測資料反演得到的C值為1.43 m0.5/s。

圖1 泄流槽入口處水力關系

1.2 沖刷計算

對潰口沖刷過程,BREACH模型采用Smart修正的Meyer-Peter & Mueller公式[9]:

(2)

式中:Qs為單寬沖刷率;D為水力半徑;S為下游坡度的倒數;P為潰口的濕周;n為曼寧系數;D30、D90分別為含量30%和90%的顆粒粒徑;Ω與臨界剪應力有關,相關公式Fread[3]在其文章中做了詳細論述,在此不做贅述。

陳祖煜等[7]在分析唐家山堰塞湖潰決過程中,發(fā)現(xiàn)指數形式的侵蝕率模型在流速較低的范圍內,適用性較好,但當潰壩水流流速達到5～10 m/s時,潰壩洪峰計算結果受參數影響極為敏感,經常導致計算結果不收斂,因此在DB-IWHR模型中提出了雙曲型侵蝕率模型,其表達式為

(3)

其中v=k(τ-τc)

圖2 基于實測資料的雙曲線和Meyer-Peter & Muller侵蝕率模型

1.3 潰口擴展計算

BREACH模型假定潰口最初為矩形,當潰口的兩邊崩潰,最初的長方形通道就變成一個梯形通道,邊坡與豎直方向夾角為α。隨著邊坡的變形破壞,采用梯形潰口進行邊坡計算。當發(fā)生崩潰時,深度消減到達臨界深度,大壩內摩擦角為φ,黏聚力為c,容重為γ。潰口兩側邊坡的穩(wěn)定性與潰口下切深度相關,潰口下切至臨界深度Hk時,潰口兩側邊坡發(fā)生失穩(wěn)坍塌,初始矩形潰口變?yōu)樘菪螡⒖?即當潰口深度H>Hk時,邊坡失穩(wěn)坍塌,k依次取1、2、3進行計算得到潰口深度(參見圖3(a))。BREACH模型為

B0=Bry

(4)

(5)

θk-1=90-α

(6)

式中:y為潰口通道的流動深度;Br為線性參數,在漫頂潰壩時,Br取2。

圖3 潰口擴展示意圖

BREACH模型另一個控制潰口寬度的機制是土體的邊坡穩(wěn)定性。在巖土工程中常用的是圓弧形式的邊坡穩(wěn)定分析方法, DB-IWHR模型的潰口擴展過程,即采用簡化Bishop法[11]自動搜索臨界圓弧滑裂面,進行潰口擴展數值分析。通過搜索可能的滑裂面,找到臨界滑裂面,從而得到下切深度,文獻[12]詳細介紹了相關原理和計算方法(圖3(b))。

1.4 堰后水深計算

水流以寬頂堰形式泄流時,在堰口會發(fā)生一個跌落現(xiàn)象(圖1)。確定堰后水深h是潰壩計算中需要獨立解決的問題。BREACH模型假定根據潰口流量,由曼寧公式可以得到潰口內水深：

(7)

式中:n為曼寧系數；J為下游坡度。經驗表明,曼寧公式中的J對計算結果影響極大[13]。陳祖煜等[7]發(fā)現(xiàn),如果將式(1)代入式(7),可得

(8)

式(8)中H-z項的指數為0.9,DB-IWHR模型將其近似取為1,定義m為水位跌落系數:

(9)

則有h=m(H-z)

(10)

式中m具有明確的物理意義和取值范圍(一般在0.5～0.8之間),后文將作詳細的敏感性分析,并與輸入坡降J來確定h的方案相對比,說明計算的穩(wěn)定性得到了明顯提高。

1.5 數值計算方法

BREACH在數值計算中,是通過給定的初始時間t0和時間步長Δt,計算相應的水位增量ΔH,沖刷深度Δz和流速變化量ΔV。

DB-IWHR模型以水流流速增量ΔV代替時間Δt為步長,給定初始流速V0和流速增量ΔV,避免迭代數值計算。該模型在計算中有較好的穩(wěn)定性,不會出現(xiàn)數值不收斂的情況。

2 唐家山潰壩過程模擬

2.1 唐家山潰決概況

2008年5月12日,中國四川省汶川發(fā)生8.0級大地震,受強烈地震的影響,位于四川通口河右岸、距北川縣城6 km處的唐家山發(fā)生了巨大的滑坡,滑坡體下切堵塞了通口河,形成了地震中最嚴重的堰塞湖。根據資料[14]唐家山堰塞壩的相關參數如下，壩頂高程741 m,壩底高程669 m,壩頂長度310 m,順河向長803 m,橫向最大寬度611 m,庫容2.3億m3,入庫流量80 m3/s,啟動流速2.6 m/s ，平均粒徑2～8 mm,初始潰口寬58 m,初始潰口深度0.18 m,潰口邊坡1.5∶1,最終潰口寬度131.2 m,最終潰口深度20.8 m,洪峰流量6 480 m3/s,潰口發(fā)展時間6～15 h,唐家山堰塞壩潰后剖面如圖4所示。

圖4 唐家山堰塞壩潰后剖面[14]

2.2 不同模型計算

不同模型計算參數如表1和表2所示。對唐家山采用此兩種模型計算獲得的流量、庫水位、潰口寬度、潰口底高程如圖5所示。

表1 BREACH模型計算參數

表2 DB-IWHR模型計算參數

注：①B0基于泄流槽的幾何參數以及3 m深的水流;②庫容水位關系由公式W=a1(H-Hr)2+b1(H-Hr)+c1擬合估算 ;③侵蝕參數為雙曲線模型中的參數。

圖5 唐家山堰塞湖兩種模型模擬結果對比

運用BREACH模型對唐家山堰塞湖潰決進行反演計算,可得到洪峰流量為6 478 m3/s,潰決時間為12:32左右。對于潰口的形態(tài),初始底部寬度設定為14 m,已經達到其泄流槽開挖寬度,之后潰口以14 m為底邊長從長方形潰口逐漸發(fā)展為梯形潰口,潰壩起始階段,潰口的下切速度較快,當入庫流量保持不變時,潰口頂部寬度也會一直發(fā)展,并不斷下切。

運用DB-IWHR模型對唐家山進行分析計算,可以得到唐家山潰決的洪峰流量為7 610 m3/s,潰決時間為13:20左右,流量過程和水位過程與實測結果差距不大,但洪峰流量到達時間相差較大。在達到洪峰流量之前,模擬結果與實測相近,洪峰流量之后則有了較明顯的偏差。

通過模型反演計算可知,使用針對唐家山堰塞湖的參數,BREACH模型和DB-IWHR模型均可較好地模擬唐家山的潰決過程。

3 潰壩模型敏感性分析

3.1 BREACH模型

BREACH模型計算潰壩時,主要參數包括壩體幾何特征和壩體材料參數,由于天然堰塞壩土體材料復雜,在應急處置時,很難及時測得所有有效參數,只能通過現(xiàn)場地質勘探進行估測,在此選擇壩體材料相關參數進行敏感性分析。

3.1.1 壩體模型參數

計算中選擇ZU=2.75,ZD=11.1,現(xiàn)分別選取不同的ZU、ZD對模型進行分析。

針對不同ZU、ZD的結果分別如圖6、圖7所示?？梢悦黠@發(fā)現(xiàn),BREACH模型中ZU的變化對結果影響不大,而ZD對結果影響較大,即改變堰塞壩不同的下游坡比,得到的結果差異性很大。分析中選取了ZD分別為6、9、10、11.1、12(即J為1/6、1/9、1/10、1/11.1、1/12),洪峰流量從32 370 m3/s變化到5 288.8 m3/s,可見,BREACH模型結果明顯受曼寧公式中J的輸入值的影響。

圖6 BREACH模型ZU敏感性分析

圖7 BREACH模型ZD敏感性分析

3.1.2 壩體材料孔隙率

圖8 BREACH模型孔隙率敏感性分析

計算中分別將孔隙率取為0.1、0.2、0.5進行分析,得到流量過程線對比見圖8。從圖8可以明顯發(fā)現(xiàn),壩體材料孔隙率對洪水過程影響較大?？紫堵嗜?.1、0.2、0.5,隨著孔隙率的增大,洪峰流量逐漸增大,而洪峰到達時間逐漸變小,洪峰流量分別是5 707 m3/s、6 478 m3/s、10 712 m3/s,變化最大達60%,選取孔隙率準確與否對模擬結果也有著重要影響。

3.1.3 壩體材料內摩擦角

計算中分別將內摩擦角φ取為35°、30°、20°,得到流量過程線對比見圖9。

圖9 BREACH模型內摩擦角敏感性分析

從圖9中可以發(fā)現(xiàn),內摩擦角對潰壩流量過程影響不大。在到達洪峰流量之前,內摩擦角越小,潰口出流量增加越快,洪峰流量到達的時間越快,且洪峰流量越大;而在到達洪峰流量之后,內摩擦角越小,潰口流量減小得也越快。

3.1.4 壩體材料黏聚力

計算中分別將黏聚力c取為40 kPa、30 kPa、20 kPa,得到流量過程線對比見圖10。

圖10 BREACH模型黏聚力敏感性分析

從圖10可以發(fā)現(xiàn),黏聚力與內摩擦角對潰壩過程的影響類似,其對流量過程的敏感性也較小。

綜上所述,BREACH模型中材料的孔隙率和下游坡比ZD兩個輸入參數是影響計算成果合理性的關鍵問題。

3.2DB-IWHR模型

3.2.1 跌落系數

如前所述,DB-IWHR模型的一個主要改進是使用跌落系數m來取代J這個輸入值,在典型工況計算時取0.8,分別將跌落系數m取為0.8、0.75、0.67,結果對比如圖11所示。

圖11 DB-IWHR模型淹跌落系數敏感性分析

從圖11可知,跌落系數對洪峰流量到達時間幾乎沒有影響；對于洪峰流量,跌落系數取值從0.8變?yōu)?.67,洪峰流量僅增加了5.8%,利用跌落系數m來取代J來計算堰后水深具有明顯的優(yōu)越性。

3.2.2 沖刷臨界啟動速度

計算中分別將沖刷臨界啟動速度Vc取為2.7 m/s、2.5 m/s、3 m/s,結果對比如圖12所示。

圖12 DB-IWHR模型沖刷臨界啟動速度敏感性分析

從圖12(a)中可以看出來,對于流量過程,調整啟動流速,整體未發(fā)生明顯變化;對于洪峰流量到達時間來看,啟動流速從2.5 m/s、2.7 m/s到3 m/s變化,時間變化均不足10%,影響很小。

陳祖煜等[7]對DB-IWHR的強度參數c、φ進行了敏感性分析考察,結果表明模型中側向擴展的邊坡穩(wěn)定分析方法對剪切強度參數并不敏感。

DB-IWHR潰壩模型的主要參數對結果的影響都不敏感,且模型參數量少、簡單易用,可以快速給出估算結果,在潰決洪水演進分析中,以及資料匱乏的初步分析階段具有重要的實用價值。

3.3 數值計算穩(wěn)定性

BREACH模型在計算時,若下游坡度小于4時,程序計算不收斂,無法得到結果。DB-IWHR模型采用直接求解格式,不包含數值迭代,計算迅速、簡明。

4 結 論

a. 本文比較了都使用寬頂堰模型基本形式的BREACH和DB-IWHR兩種潰壩洪水計算模型。DB-IWHR模型實際上是在BREACH模型理論框架下發(fā)展起來的,主要在以下幾個方面做了改進:①采用雙曲線沖刷模型代替了Smart修正的 Meyer-Peter & Muller公式;②使用跌落系數取代曼寧公式計算堰后水深以解決坡降J對計算結果高度敏感的問題;③應用圓弧滑動模式較合理地模擬壩體潰決橫向擴展;④應用速度積分代替時間積分,回避了非線性迭代計算,計算過程可以用簡單的Excel表格實現(xiàn)。

b. 對唐家山堰塞湖潰決反演的計算表明,使用針對唐家山堰塞湖的參數,兩種模型均可較好地模擬唐家山堰塞湖的潰決過程。

c. 對兩種模型輸入參數的敏感性分析表明:①在BREACH模型中,孔隙率的影響比較敏感,在實際使用范圍內取值,變動范圍最大達到60%;而在DB-IWHR的雙曲線模型中,參數1/a、1/b需要通過大量沖刷實驗積累經驗進行取值。②堰后水深的確定,BREACH模型對坡降J極其敏感,DB-IWHR模型使用跌落系數m后,對結果影響很小。③橫向擴展過程中,BREACH模型和DB-IWHR模型中黏聚力c、內摩擦角φ對結果都不敏感。④對于數值計算穩(wěn)定性,BREACH模型在合理選取J的基礎上計算穩(wěn)定,DB-IWHR模型則不會出現(xiàn)不穩(wěn)定的情況。

[1] 劉寧.唐家山堰塞湖應急處置與減災管理工程[J].中國工程科學,2008,10(12):67-72.(LIU Ning. The emergency handling of Tangjiashan Barrier Lake and disaster reduction management project [J].Engineering Sciences, 2008,10(12):67-72.(in Chinese))

[2] LIU F,FU X D,WANG G Q,et al.Physically based simulation of dam breach development for Tangjiashan Quake Dam, China[J].Environmental Earth Sciences,2012,65(4):1081-1094.

[3] FREAD D L. BREACH, an erosion model for earthen dam failures [D]. Silver Spring:Hydrologic Research Laboratory, National Weather Service, NOAA, 1988.

[4] SINGH V. Dam breach modeling technology [M]. Netherlands:Springer, 1996.

[5] 周興波,陳祖煜,陳淑婧,等.基于MIKE11的堰塞壩潰決過程數值模擬[J].安全與環(huán)境學報,2014,14(6):23-27.(ZHOU Xingbo,CHEN Zuyu,CHEN Shujing, et al. Simulation of the barrier dam breaching process based on MIKE 11[J]. Journal of Safety and Environment,2014,14(6):23-27. (in Chinese))

[6] Hydrologic Engineering Center of USACE. HEC-RAS: River analysis system hydraulic reference mannual version 4.1 [D].Davis C A: US Army Corps of Engineers,2010.

[7] CHEN Z,MA L,YU S,et al. Back analysis of the draining process of the Tangjiashan Barrier Lake [J].Journal of Hydraulic Engineering,2014,141(4):1-14.

[8] 王琳,李守義,于沭,等. 紅石巖堰塞湖應急處置的關鍵技術[J].中國水利水電科學研究院學報,2015,13(4):284-289.(WANG Lin,LI Shouyi,YU Shu,et al. Key techniques for the emergency disposal of Hongshiyan Landslide Dam[J]. Journal of China Institute of Water Resources and Hydropower Research,2015,13(4):284-289. (in Chinese))

[9] MEYER-PETER E, MULLER R. Formulas for bed-load transport [C]

Proceedings of the 2nd Meeting of the International Association for Hydraulic Structures Research. Stockholm:IAHR, 1948:39-64.

[10] 周興波,陳祖煜,李守義,等.不同推移質輸沙模型在潰壩洪水模擬中的對比分析[J].應用基礎與工程科學學報,2015,23(6):1097-1108.(ZHOU Xingbo,CHEN Zuyu,LI Shouyi,et al. Comparison of sediment transport model in dam break simulation[J]. Journal of Basic Science and Engineering,2015,23(6):1097-1108. (in Chinese))

[11] 陳祖煜.土質邊坡穩(wěn)定分析:原理·方法·程序[M].北京:中國水利水電出版社,2003.

[12] WANG Lin, CHEN Zuyu, WANG Naixin,et al. Modeling lateral enlargement in dam breaches using slope stability analysis based on circular slip mode [J]. Engneering Geology,2016,209:70-81.

[13] 陳淑婧.土石壩漫頂潰決機理模型及數值模擬方法研究[D].北京:北京工業(yè)大學,2015.

[14] LIU N,CHEN Z,ZHANG J,et al. Draining the Tangjiashan Barrier Lake [J].Journal of Hydraulic Engineering, 2010,136(11):914-923.

Comparison of two models for back analysis of dam breach of Tangjiashan Barrier Lake

LI Xiangnan, CHEN Zuyu

(DepartmentofGeotechnicalEngineering,ChinaInstituteofWaterResourcesandHydropowerResearch,Beijing100048,China)

Based on measured data of the dam breach of Tangjiashan Barrier Lake, the process of the dam breach was analyzed using the BREACH model and the DB-IWHR model developed by the China Institute of Water Resources and Hydropower Research. A sensitivity analysis of the two models to different parameters was also performed. Results show that, using the parameters of the dam, the two models can both reflect the flooding process induced by the dam breach; the calculated results of the peak flow are more sensitive to the input parameters of the downstream slope and porosity in the BREACH model; the DB-IWHR model improves the calculation of the erosion parameter and downstream water depth, resulting in calculated results less sensitive to the input parameters; and the DB-IWHR model provides software in the Excel format, with a high stability in numerical analysis, which is transparent, self-explanatory, and friendly to practitioners.

barrier lake; BREACH model; DB-IWHR model; dam breach-induced flood analysis; Tangjiashan

國家重點基礎研究發(fā)展計劃(973計劃)(2013CB036402);國家自然科學基金青年基金(51309260)

李相南(1992—),男,碩士研究生,主要從事潰壩模擬及洪水計算研究。E-mail:lixn0555@163.com

陳祖煜(1943—),男,中國科學院院士,博士,主要從事水利工程和巖土工程研究。E-mail:chenzuyu@iwhr.com

10.3880/j.issn.1006-7647.2017.02.004

TV122+.4

：A

：1006-7647(2017)02-0020-07

2016-08-22 編輯：鄭孝宇)