土石壩潰壩洪水分析：原理和計算程序

陳祖煜，陳淑婧，王 琳，張 強

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

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

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

1 基本原理

1.1 潰壩水力學計算基礎

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

(1)

(2)

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

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

(3)

圖1 寬頂堰泄流示意圖

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

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

(4)

潰口流量可以通過單位時間內水庫庫容的損失來確定，因此：

(5)

式中：q為入庫流量。

1.2 沖刷模型

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

(6)

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

v=k(τ-τc)

(7)

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

表1 雙曲線沖刷模型參數參考值

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

(8)

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

1.3 潰口側向擴展模型

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

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

圖2 潰口的側向崩塌過程等效簡化

(9)

(10)

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

圖3 經驗系數m1的回歸方程

圖4 經驗系數m2的回歸方程

2 數值計算方法

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

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

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

3 潰壩洪水分析程序

3.1 程序界面

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

DB-IWHR 2018需要輸入表2中有關地形地貌、水力學和巖土力學參數。

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

表 2 DB-IWHR需要輸入的參數列表

4.2 成果分析

有關唐家山堰塞湖潰決洪水計算成果分析，已在文獻[1]中詳細討論?？傮w的結論是，表2中侵蝕參數a,b是影響洪峰計算值最敏感的參數，需要在實踐中不斷摸索，深化認識。

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

(11)

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

為分析相關參數對計算結果的影響，DB-IWHR模型計算中對上述曼寧公式做了進一步的推導，提出使用一個跌落系數m簡化潰口水深計算，具體推到公式如下：

(12)

(13)

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

圖5 數值分析方法的流程表

圖6 DB-IWHR2018程序界面

表3 BREACH 模型模擬唐家山堰塞湖潰決流量的輸入參數

圖7 BREACH模型中下游坡比J的敏感性分析結果(輸入參數見表3)

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

線性模型具體形式如下：

(14)

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

圖8 線性模型參數敏感性分析結果

4 結 論

本文介紹中國水利水電科學研究院近期對土石壩潰決洪水分析方法的改進以及相應的軟件開發工作。對數值分析方法的改進主要為以下方面：

(1)提出了一個描述土體沖刷的雙曲線模型。此模型的兩個參數具有物理意義。參數1/a代表起動區沖刷速率依線性關系相對剪應力的比例系數，可在實驗室內測定。參數1/b代表在高流速區沖刷速率逼近“可能最大”的漸近值。表2提供了參考值。隨著經驗的積累，這些參考值可不斷更新。這一模型可以使計算的穩定性大大提高，不會出現因輸入參數微小變化導致洪峰值急劇變動的現象。

(2)對于潰口側向崩塌擴展也提出了一個經驗的雙曲線模型?？梢愿鶕馏w的容重、凝聚力和摩擦角確定相應的參數m1和m2。

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

在上述改進的基礎上開發了潰壩洪水分析程序DB-IWHR具有以下特點：

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

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

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