陳祖煜，陳淑婧，王 琳，張 強(qiáng)

(1. 中國水利水電科學(xué)研究院, 北京 100048；2. 全國市長研修學(xué)院, 北京 100029；3. 西安理工大學(xué) 水利水電學(xué)院, 陜西 西安 710048)

土石壩潰決洪水分析是預(yù)警和應(yīng)急搶險(xiǎn)工作中非常重要的基礎(chǔ)工作，潰壩流量預(yù)測的準(zhǔn)確與否可直接影響應(yīng)急處置方案的制定與實(shí)施。這一分析工作涉及對非恒定高速急變流、泥沙動(dòng)力學(xué)和巖土力學(xué)多領(lǐng)域耦合的模擬過程。過去50 a來，許多學(xué)者研究開發(fā)基于物理機(jī)制的潰壩模型。但是由于其本身包涵的復(fù)雜問題以及缺乏較多工程案例驗(yàn)證的原因，尚有待繼續(xù)努力，以形成一門完善、成熟的學(xué)科分枝。

中國水利水電科學(xué)研究院在唐家山堰塞湖成功引流除險(xiǎn)后，根據(jù)實(shí)測資料開展了反演分析。在此基礎(chǔ)上對現(xiàn)有的潰壩洪水分析模型進(jìn)行了改進(jìn)。相應(yīng)的研究成果發(fā)表于美國ASCE J. Hydraulic Engineering[1]。 隨后，開發(fā)了潰壩洪水分析DB-IWHR。這一程序在Excel界面上開發(fā)，利用其強(qiáng)大的計(jì)算功能并結(jié)合VBA自主進(jìn)行開發(fā)，形成具有一定創(chuàng)新性的計(jì)算潰壩洪水過程線的電子計(jì)算表格，是一個(gè)面向一線工程技術(shù)人員的實(shí)用軟件。課題組又對易貢、小崗劍、板橋、白格等具有實(shí)測資料的堰塞湖潰決實(shí)例進(jìn)行了分析[2-5]，已在我國水利水電行業(yè)獲得了初步應(yīng)用。鑒于相關(guān)的研究論文多為英文，本文綜合前期的研究工作的要點(diǎn)，并對一些細(xì)節(jié)做進(jìn)一步的說明，以供國內(nèi)用戶參考。

1 基本原理

1.1 潰壩水力學(xué)計(jì)算基礎(chǔ)

潰口斷面的流量可用寬頂堰公式計(jì)算。

(1)

(2)

式中：mb,mq分別為寬頂堰的流量系數(shù)和側(cè)收縮系數(shù)。C為綜合流量系數(shù)，理論值為1.7 m1/2/s[6]，可在1.3到1.7范圍內(nèi)選取[7]。H和z分別為庫水位和渠底高程，m。見圖1。

水流經(jīng)過堰口跌落后水深為h，DB-IWHR采用簡化處理方案。建議在0.8(0.6之間假定一個(gè)跌落系數(shù)：

(3)

圖1 寬頂堰泄流示意圖

式中：h為潰口水深，m?？梢酝ㄟ^敏感性分析考察不同的m值對計(jì)算結(jié)果的影響。經(jīng)驗(yàn)表明，m的假定值對洪峰的計(jì)算值影響很小。

潰口斷面流速V可按式(4)確定：

(4)

潰口流量可以通過單位時(shí)間內(nèi)水庫庫容的損失來確定，因此：

(5)

式中：q為入庫流量。

1.2 沖刷模型

DB-IWHR建議了一個(gè)雙曲線模型，其形式如式(6)：

(6)

式中：v為扣除臨界剪應(yīng)力后的剪應(yīng)力。

v=k(τ-τc)

(7)

雙曲線有一當(dāng)v接近無限值時(shí)的漸近線，其值為1/b，代表“可能最大沖刷速率”。此處，k取100，1/a表示v等于0時(shí)曲線的斜率。該模型基于結(jié)構(gòu)材料的理解而建立的，即土體材料抵抗侵蝕時(shí)，不應(yīng)有無限“強(qiáng)度”。 表1提供相應(yīng)的參考值。

表1 雙曲線沖刷模型參數(shù)參考值

計(jì)算沖刷速率需要確定剪應(yīng)力，可按式(8)：

(8)

式中：γ為水的重度；R為水力半徑；S為能坡。

1.3 潰口側(cè)向擴(kuò)展模型

潰口底部不斷地被刷深的過程中，兩側(cè)邊坡發(fā)生崩塌失穩(wěn)，側(cè)面不斷地?cái)U(kuò)大。以前的潰壩分析模型采用楔形體法[8-10]，這一做法過于簡單，水科院團(tuán)隊(duì)采用了巖土工程中已經(jīng)被廣泛接受的圓弧滑動(dòng)面分析方法，并且比較合理地模擬了由于坡角下切導(dǎo)致的巖坡崩塌過程。相關(guān)的分析成果形成如圖2(a)所示階梯式潰決模式。但是在實(shí)際操作時(shí)，需輸入每一步的圓心坐標(biāo)、半徑和渠底坐標(biāo)，過于繁雜。在實(shí)際應(yīng)用時(shí)，仍將計(jì)算獲得的圖形簡化為一系列梯形，圖2(b)示。

近期的研究表明，在可接受的精度范圍內(nèi)，可用一個(gè)雙曲線模型來計(jì)算梯形側(cè)面傾角增量Δβ，使用以下公式即可確定。

圖2 潰口的側(cè)向崩塌過程等效簡化

(9)

(10)

式中：1/m1和1/m2分別表示雙曲線的初始切線和漸近線。文獻(xiàn)[11]提供了相應(yīng)的圖表以確定此兩個(gè)參數(shù)。在更新DB-IWHR程序時(shí)，我們還進(jìn)一步擬合了如圖3和圖4所示的線性關(guān)系。根據(jù)不同壩體土料的材料重度γ和強(qiáng)度參數(shù)內(nèi)聚力c和內(nèi)摩擦角φ，通過內(nèi)插得到雙曲線模型參數(shù)m1和m2。經(jīng)驗(yàn)表明，凝聚力c對m2的數(shù)值影響不大。故圖4僅依據(jù)摩擦系數(shù)tanφ來確定m2。

圖3 經(jīng)驗(yàn)系數(shù)m1的回歸方程

圖4 經(jīng)驗(yàn)系數(shù)m2的回歸方程

2 數(shù)值計(jì)算方法

常規(guī)的計(jì)算方法[8-9,12]都是通過給定的初始時(shí)間t0和時(shí)間步長Δt，計(jì)算相應(yīng)的水位增量ΔH，沖刷深度Δz和流速變化量ΔV。相應(yīng)的算法是高度非線性的。計(jì)算程序包含復(fù)雜的迭代求解過程，常會(huì)遇見收斂困難。通過觀察發(fā)現(xiàn)，一旦給定流速V，可以直接求出相應(yīng)的ΔH、Δz和Δt。因此新的算法采用給定初始流速V0和流速增量ΔV。

DB-IWHR2018是在Microsoft Excel 2007中用VBA語言編寫的程序，該程序能快速計(jì)算出大壩潰決的峰值流量。程序可在下面網(wǎng)站下載：

http://www.geoeng.iwhr.com/ytgcyjs/czy/zlxz/kbfx/A54160106index_1.htm.

3 潰壩洪水分析程序

3.1 程序界面

DB-IWHR包括主界面“Calculation”和庫容水位關(guān)系計(jì)算“W-H curve”。近期，作者增加了一個(gè)引導(dǎo)界面，用于初學(xué)者在一系列提示下輸入十余個(gè)參數(shù)(多數(shù)可取默認(rèn)值)，迅速進(jìn)入主界面計(jì)算。

DB-IWHR 2018需要輸入表2中有關(guān)地形地貌、水力學(xué)和巖土力學(xué)參數(shù)。

圖6所示唐家山堰塞湖潰決洪水計(jì)算界面。

表 2 DB-IWHR需要輸入的參數(shù)列表

4.2 成果分析

有關(guān)唐家山堰塞湖潰決洪水計(jì)算成果分析，已在文獻(xiàn)[1]中詳細(xì)討論?？傮w的結(jié)論是，表2中侵蝕參數(shù)a,b是影響洪峰計(jì)算值最敏感的參數(shù)，需要在實(shí)踐中不斷摸索，深化認(rèn)識(shí)。

(1)對淹沒系數(shù)計(jì)算方法的改進(jìn)。Fread[8]在BREACH模型中建議潰口水深計(jì)算可按明渠均勻流計(jì)算，即：

(11)

式中：n為曼寧系數(shù)；B為潰口寬度；J為下游坡坡比。在使用此式和相關(guān)的Breach程序時(shí)，需要輸入下游坡坡比J的數(shù)值。圖7是輸入不同的J值獲得的洪水過程，可見，J值對計(jì)算成果極為敏感。下游坡比J從0.05升高到0.24，潰口峰值流量增大24倍。

為分析相關(guān)參數(shù)對計(jì)算結(jié)果的影響，DB-IWHR模型計(jì)算中對上述曼寧公式做了進(jìn)一步的推導(dǎo)，提出使用一個(gè)跌落系數(shù)m簡化潰口水深計(jì)算，具體推到公式如下：

(12)

(13)

由上式(13)可知系數(shù)m具有一定的物理意義，即為水流從水庫至潰口處水頭跌落，故定義為跌落系數(shù)。文獻(xiàn)[1]的敏感性分析也考察不同的m值對計(jì)算結(jié)果的影響，結(jié)果表明，m的假定值對洪峰的計(jì)算值影響很小。因此，DB-IWHR不使用式(13),而是建議直接為m提供一個(gè)默認(rèn)值，如0.8。用戶可對其作敏感性分析，考察不同的m是否確實(shí)對計(jì)算成果影響不大。

圖5 數(shù)值分析方法的流程表

圖6 DB-IWHR2018程序界面

表3 BREACH 模型模擬唐家山堰塞湖潰決流量的輸入?yún)?shù)

圖7 BREACH模型中下游坡比J的敏感性分析結(jié)果(輸入?yún)?shù)見表3)

(2)雙曲線模型和線性模型的比較。文獻(xiàn)[1]根據(jù)實(shí)測沖刷率數(shù)據(jù)，采用雙曲線沖刷模型進(jìn)行擬合，以計(jì)算唐家山堰塞湖在潰決過程中的沖刷速率，潰決流量的計(jì)算結(jié)果顯示雙曲線模型參數(shù)大范圍變動(dòng)對流量計(jì)算結(jié)果不很敏感，并且，也發(fā)現(xiàn)雙曲線模型對計(jì)算成果十分敏感。本文進(jìn)一步論證線性沖刷模型對計(jì)算成果同樣十分敏感，不宜直接在潰壩洪水分析中應(yīng)用。

線性模型具體形式如下：

(14)

式中：a1是線性沖刷系數(shù)。唐家山堰塞湖實(shí)測沖刷率的線性回歸結(jié)果如圖8(a)所示，對應(yīng)圖中a1分別取0.1,0.2,0.3時(shí)，潰決流量曲線如圖8(b)所示，峰值流量依次為4219 m3/s，9024 m3/s，14 723 m3/s。計(jì)算結(jié)果說明線性模型系數(shù)小幅度變化時(shí)對計(jì)算峰值流量的影響卻很大，因此若采用線性模型增加了合理選取參數(shù)的難度，從而降低了潰壩流量模型計(jì)算結(jié)果的可靠性。

圖8 線性模型參數(shù)敏感性分析結(jié)果

4 結(jié) 論

本文介紹中國水利水電科學(xué)研究院近期對土石壩潰決洪水分析方法的改進(jìn)以及相應(yīng)的軟件開發(fā)工作。對數(shù)值分析方法的改進(jìn)主要為以下方面：

(1)提出了一個(gè)描述土體沖刷的雙曲線模型。此模型的兩個(gè)參數(shù)具有物理意義。參數(shù)1/a代表起動(dòng)區(qū)沖刷速率依線性關(guān)系相對剪應(yīng)力的比例系數(shù)，可在實(shí)驗(yàn)室內(nèi)測定。參數(shù)1/b代表在高流速區(qū)沖刷速率逼近“可能最大”的漸近值。表2提供了參考值。隨著經(jīng)驗(yàn)的積累，這些參考值可不斷更新。這一模型可以使計(jì)算的穩(wěn)定性大大提高，不會(huì)出現(xiàn)因輸入?yún)?shù)微小變化導(dǎo)致洪峰值急劇變動(dòng)的現(xiàn)象。

(2)對于潰口側(cè)向崩塌擴(kuò)展也提出了一個(gè)經(jīng)驗(yàn)的雙曲線模型?？梢愿鶕?jù)土體的容重、凝聚力和摩擦角確定相應(yīng)的參數(shù)m1和m2。

(3)提出了一個(gè)以速度增量為基礎(chǔ)的數(shù)值積分模式，回避了非線性迭代，使全部計(jì)算線性化。

在上述改進(jìn)的基礎(chǔ)上開發(fā)了潰壩洪水分析程序DB-IWHR具有以下特點(diǎn)：

(1)程序是在Excel界面上開發(fā)的計(jì)算表格。使用者只需在引導(dǎo)表格輸入15個(gè)參數(shù)。其友好的界面使絕大多數(shù)用戶不需要接受過于繁復(fù)的培訓(xùn)，即可開展計(jì)算工作。

(2)鑒于其線性化的處理，計(jì)算過程快速，魯棒性強(qiáng)。

(3)全部計(jì)算過程開源，透明。便于校核和二次開發(fā)。