, ,

(長江科學院 水力學研究所， 武漢 430010)

1 研究背景

研究土石壩的潰決過程及其壩下游洪水演進，對于我國處置潰壩突發(fā)性洪水災害，提升我國水利安全，具有十分重大的意義[1]。

土石壩的潰決涉及到水力學、土力學以及泥沙動力學等多學科，受壩體材料、壩體型式與尺度、水庫庫容及上下游水位等多種因素的影響。從20世紀60年代開始，世界各國相繼開始了對潰壩問題的研究，逐步形成了如“跌坎沖刷(headcut erosion)”、“切溝侵蝕(gully erosion)”等能較為合理地揭示潰壩過程的理論，為潰壩洪水預測奠定了良好理論基礎并提供了豐富的驗證資料[2-5]。

由于潰壩洪水的巨大危害性，世界各國均開發(fā)了計算模型用于預測與預報，如基于潰口線性發(fā)展假設的美國國家氣象局DAMBRK模型與SMPDBK模型、基于泥沙輸移理論并考慮潰口邊坡坍塌的Beed與Breach等模型[4]。近年來潰壩洪水平面二維數(shù)值模型的研究取得了長足的進展。Fracacrollo[6]建立了概化的室內(nèi)潰壩物理模型,開展了試驗研究并進行了數(shù)值模擬，潰口采用瞬潰假定；崔丹等[7]基于有限體積法與具有總變差減小特性MacCormack格式建立潰壩洪水演進數(shù)值模型，模型中潰口的發(fā)展采用了線性假定；夏軍強等[8]建立了基于無結構三角網(wǎng)格下采用有限體積法求解的二維水動力學模型，潰口采用瞬潰假定；楊志等[9]建立了黑河金盆水庫大壩潰口近區(qū)二維數(shù)值模型和下游地區(qū)潰壩洪水演進耦合數(shù)學模型，模型中潰口的發(fā)展及潰口流量采用公式計算；岳志遠等[10]建立二維水沙耦合數(shù)學模型，采用修正的上揚通量公式來描述泥沙運動和滑坡堰塞體潰決過程，將潰壩過程與洪水演進耦合模擬，相比前述模型有了較大的進展，但未考慮潰口發(fā)展過程中潰口邊坡的坍塌。

潰口流量過程是壩下洪水預測的核心參數(shù)，其精度不夠可能導致預警過度或不足。潰口發(fā)展過程是決定潰口流量過程的關鍵因素，它與潰口水力要素相互耦合，具有很強的時變特性與空間分布的不均性。合理地描述壩體潰決過程中潰口“剪剝”式(或“跌坎”式)發(fā)展模式、壩面雙螺旋流淘刷引起潰口橫向展寬以及邊坡失穩(wěn)坍塌的現(xiàn)象，需要詳細的潰口區(qū)流場隨時間變化的信息，如水深分布、流速分布等。上述計算模型由于其一維性很難體現(xiàn)潰口發(fā)展的空間分布不均勻性，而基于水深積分的平面二維數(shù)值模型雖不能反映潰口區(qū)水流的三維特性，但能提供的水位(水深)場、平面流場及對不規(guī)則潰口形態(tài)的描述能力等，能較好地滿足“跌坎”移動速度、壩體下切速度以及邊坡失穩(wěn)坍塌計算所需的水力參數(shù)，從而實現(xiàn)對潰口復雜沖蝕過程的模擬。相比一維數(shù)值模型僅采用潰口斷面平均流速來計算潰口發(fā)展過程有著本質(zhì)的區(qū)別，因此構建壩體潰決過程與潰壩洪水演進耦合的數(shù)值模型是潰壩洪水預測的發(fā)展趨勢。此外，潰口區(qū)間斷波的捕捉、洪水演進過程中引起的干濕邊界的處理、模型的守恒性等亦是潰壩數(shù)學模型中需要解決的重點問題。

本文采用雙曲線型沖蝕速率表達式描述壩體沖蝕、簡化Bishop法搜索臨界滑裂面描述潰口邊坡坍塌、具有總變差不增特性的MacCormack有限體積法離散控制方程，建立了壩體潰決過程與潰壩水流完全耦合維數(shù)值模型，相比采用經(jīng)驗公式計算潰口發(fā)展[9]、瞬潰假定[6、8]、線性潰決假定[7]和不考慮潰口邊坡坍塌[10]的模型，潰口流量過程的合理性及壩下洪水預測精度均可得到較大的提高。

2 數(shù)學模型

2.1 基本方程

(1)

(2)

將控制方程統(tǒng)一寫成守恒型式，即

(3)

式中：U為守恒性變量，U=(h,hu,hv)；F和T分別代表對流項與擴散項；S代表源項。

2.2 潰口沖蝕

對于黏性和非黏性土體的沖蝕率與切應力之間關系，國內(nèi)外學者已進行了大量的研究， Gauche基于水槽試驗在2010年提出了針對非黏性土的指數(shù)形式的沖蝕速率表達式,即

(4)

(5)

(6)

式中：γ為水體重度(9.8 kN/m3)；J為水力坡度；R′為水力半徑，當渠道寬深比足夠大時，可近似取水深h;V為水流流速。

2.3 潰口橫向展寬

潰決過程中，當側邊坡的重力以及由于滲透造成的滲透力大于土體的摩擦力和黏聚力時，潰口側邊坡將變得不穩(wěn)定而坍塌?，F(xiàn)在大多數(shù)的模型均采用直線型的滑楔體模型，如Breach模型、Beed模型、Osman模型均屬于此類。本文潰口展寬采用式(7)所示的簡化Bishop法搜索臨界滑裂面模型。

(7)

式中：ci為第i條土條塊的黏聚力(kN/m2)；bi為第i條土條塊的寬度(m)；Wi為第i條土條塊的重力(kN)；φi為第i條土條塊的內(nèi)摩擦角(°)；θi為第i條土條塊與水平面的夾角(°)；rw為水體重度(kN/m3)；hti為浸潤線以下土條高度(m)；FB為土坡安全系數(shù)。

3 求解方法

3.1 離散格式

由于潰口區(qū)水流非恒定性強，水位梯度大，需要引入較大的數(shù)值黏性才能保持格式的數(shù)值穩(wěn)定，而數(shù)值黏性的增大會導致潰壩波峰值被抹平，降低模擬精度。為了提高計算的穩(wěn)定性模擬精度，基于MacCormack有限體積法，引入間斷波捕捉格式，在其校正步實施通量修正。目前TVD,ENO,AUSM等格式得到了廣泛的應用，引入由Yee修正后的Harten TVD格式[12-13]，即：

ΔtS(t+Δt) ；

(8)

(9)

式中：目標p,c分別代表預測步和校正步;F代表界面法向通量；A代表界面法向面積；Δt為時間步長；ΔV為控制體體積;i為網(wǎng)格邊數(shù)。

式(9)中

(11)

ψ(ak+1/2+γk+1/2)αk+1/2。

(12)

式中:ak+1/2為右特征值；αk+1/2為特征列向量元素，利用函數(shù)ψ對αk+1/2進行熵修正,ψ的表達式為

式中δ為一極小正數(shù)。式(12)中γk+1/2的表達式為

(14)

式(12)和式(14)中限制函數(shù)gk采用minmod限制器，即

(15)

3.2 網(wǎng)格劃分與干濕邊界

采用非規(guī)則局部不連續(xù)四邊形結構網(wǎng)格劃分計算區(qū)域[7]，使模型具備對復雜地形的適應能力，同時可減少無效計算網(wǎng)格，提高模型的計算效率。針對潰壩洪水波在傳播過程中，引起計算網(wǎng)格劇烈干濕變化，采用了與數(shù)值解法相匹配的基于單元屬性與基于界面屬性的動邊界處理技術[7]。

后來，宴姝去英國名校利茲大學念書。英國的博物館資源豐富，上學的3年間，宴姝最大的印象就是跟著老師跑了各種博物館、畫廊。

4 工程算例

4.1 工程概況

洈水水庫位于湖北省西南部松滋洈水鎮(zhèn)，與湖南省澧縣接壤，壩址處于洈水流域中游洈水鎮(zhèn)，距松滋城區(qū)約40 km。洈水水庫主壩為土壩，最大壩頂長1 640.00 m，最大壩高42.95 m，擋水壓力大，一旦發(fā)生潰決，潰壩洪水將對下游形成嚴重威脅。

4.2 成果分析

4.2.1 潰口沖蝕與發(fā)展過程

圖1為沿潰口中心線沖蝕縱剖面過程圖，可以看出，潰口形成1～2 min，水流主要沖刷壩軸線下游約55 m、高程71 m以上的壩面； 3 min以后，壩軸線下游約35 m處附近的壩面沖刷加劇，約5 min沖蝕至65 m高程，隨后快速向上、下游發(fā)展。

圖1 潰口中心線沖蝕過程縱剖面Fig.1 Longitudinal process lines of erosion along breach center line

圖2為潰口沖蝕與發(fā)展平面過程。由沖蝕發(fā)展過程可知，0～4 min壩面潰口逐步展寬，同一時刻從壩頂至壩腳潰口寬度沿程減??；6 min左右，壩頂至1/2壩高范圍的潰口發(fā)展較快，其寬度相對其上、下游較大；8 min左右，上、下游壩面的潰口寬度迅速增大，至10 min左右形成方向相對的“喇叭”型潰口，且沿潰口中心線，壩體均已潰至65 m的壩腳高程；15～35 min潰口向左右兩岸展寬，至39 min潰口基本穩(wěn)定。

圖2 潰口沖蝕與發(fā)展過程平面圖Fig.2 Plane diagrams of breach erosion and developing process

4.2.2 水面過程與穩(wěn)定性

潰口區(qū)0～20 min瞬時水面見圖3(圖中灰色為水面)，可以看出，隨著潰口的逐漸展寬，潰口下泄流量迅速增大，4～8 min，潰壩洪水快速淹沒壩下低洼河槽，隨即淹沒河漫灘；10 min左右，洪水上溯至兩溢洪道出口處。

圖3 潰口區(qū)瞬時水面圖Fig.3 Instantaneous water surface images of breach zone

觀察潰壩洪水的演進可以看出，模型較好地模擬了蓄洪區(qū)在洪水到達時陸域迅速被淹沒的過程，以及洪峰過后地面逐漸出露的過程，表明模型的動邊界技術是合理可靠的。

4.2.3 斷面流量過程

95.77 m庫水位條件下主壩潰決后潰口及壩下游各斷面流量過程線見圖4。從圖4可以看出，主壩潰決后潰口處流量迅速上升并達到39 074 m3/s的峰值流量，在峰值附近短時振蕩，隨后逐漸減小。壩下不同里程5個斷面監(jiān)測到的洪峰流量沿程依次減小，至壩下22 km處的金雞山斷面，洪峰流量減小至6 166 m3/s。各斷面流量過程線較為光滑，除峰值附近短時低頻振蕩外，未出現(xiàn)高頻振蕩，表明數(shù)學模型的穩(wěn)定性較強。

圖4 壩下游各斷面流量過程Fig.4 Process lines of discharge in different sections in the downstream of dam

4.2.4 守恒性分析

95.77 m庫水位條件下，庫區(qū)的初始水量8.347 7×107m3，模擬計算3 h后統(tǒng)計計算區(qū)域內(nèi)的水量為8.339 8×107m3，水量誤差約-1‰。

5 結 論

本文建立的耦合數(shù)值模型合理地模擬了壩體的潰決過程與壩下洪水演進過程，潰口的沖蝕與發(fā)展符合土石壩體潰決的一般規(guī)律；模型在潰口附近的急緩流過渡區(qū)表現(xiàn)出對間斷波良好的捕捉能力，在較大水位梯度區(qū)表現(xiàn)了良好的穩(wěn)定性；采用的動邊界處理技術對干濕劇烈變化區(qū)域具有較強的模擬能力；模型的水量誤差約1‰，守恒性良好。

參考文獻：

[1] 劉 寧.21世紀中國水壩安全管理、退役與建設的若干問題[J].中國水利,2004，(23):27-30.

[2] HANSON G J, ROBINSON K M, COOK K R. Prediction of Headcut Migration Using a Deterministic Approach[J].Transactions of ASAE, 2001, 44(3):525-531.

[3] TEMPLE D M, MOORE J S. Headcut Advance Prediction for Earth Spillways[J]. Transactions of the ASAE, 1997, 40(3):557-562.

[4] 朱勇輝,廖鴻志,吳中如. 國外土壩潰壩模擬綜述[J].長江科學院院報,2003，20(2):26-29.

[5] 陳生水, 方緒順, 鐘啟明, 等. 土石壩漫頂潰壩過程離心模型試驗與數(shù)值模擬[J]. 巖土工程學報, 2014,(5):922-932

[6] FRACCAROLLO L, TORO E F. Experimental and Numerical Assessment of the Shallow Water Model for Two-dimensional Dam-break Type Problems[J]. Journal of Hydraulic Research, 1995, 33(6):843-864.

[7] 崔 丹,槐文信,姜治兵. 潰壩洪水演進的數(shù)值模擬[J]. 華中科技大學學報(自然科學版), 2012，(7):34-37.

[8] 夏軍強,王光謙,LIN Bin-liang，等.復雜邊界及實際地形上潰壩洪水流動過程模擬[J].水科學進展,2010,21(3):289-298.

[9] 楊 志，馮民權. 潰口近區(qū)二維數(shù)值模擬與潰壩洪水演進耦合[J]. 水利水運工程學報, 2015,(1):8-19.

[10] 岳志遠,曹志先,閆 軍.滑坡體潰決洪水數(shù)學模型研究[J].水動力學研究與進展,2008,23(5):492-500.

[11] 李相南,陳祖煜. 兩種潰壩模型在唐家山堰塞湖潰決模擬中的對比[J].水利水電科技進展, 2017，(2):20-26.

[12] YEE H C. Explicit and Implicit Multidimensional Compact High-resolution Shock-capturing Methods: Formulation[J]. Journal of Computational Physics, 1997,131(1):216-232.

[13] YEE H C, SJ?GREE B. Adaptive Filtering and Limiting in Compact High Order Methods for Multiscale Gas Dynamics and MHD Systems[J]. Computers & Fluids, 2008, 37(5):593-619.