, ,

(長(zhǎng)江科學(xué)院 水力學(xué)研究所， 武漢 430010)

1 研究背景

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

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

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

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

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

2 數(shù)學(xué)模型

2.1 基本方程

(1)

(2)

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

(3)

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

2.2 潰口沖蝕

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

(4)

(5)

(6)

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

2.3 潰口橫向展寬

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

(7)

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

3 求解方法

3.1 離散格式

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

ΔtS(t+Δt) ；

(8)

(9)

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

式(9)中

(11)

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

(12)

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

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

(14)

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

(15)

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

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

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

4 工程算例

4.1 工程概況

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

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壩面潰口逐步展寬，同一時(shí)刻從壩頂至壩腳潰口寬度沿程減小；6 min左右，壩頂至1/2壩高范圍的潰口發(fā)展較快，其寬度相對(duì)其上、下游較大；8 min左右，上、下游壩面的潰口寬度迅速增大，至10 min左右形成方向相對(duì)的“喇叭”型潰口，且沿潰口中心線，壩體均已潰至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瞬時(shí)水面見圖3(圖中灰色為水面)，可以看出，隨著潰口的逐漸展寬，潰口下泄流量迅速增大，4～8 min，潰壩洪水快速淹沒壩下低洼河槽，隨即淹沒河漫灘；10 min左右，洪水上溯至兩溢洪道出口處。

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

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

4.2.3 斷面流量過程

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

圖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，模擬計(jì)算3 h后統(tǒng)計(jì)計(jì)算區(qū)域內(nèi)的水量為8.339 8×107m3，水量誤差約-1‰。

5 結(jié) 論

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

參考文獻(xiàn)：

[1] 劉 寧.21世紀(jì)中國(guó)水壩安全管理、退役與建設(shè)的若干問題[J].中國(guó)水利,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] 朱勇輝,廖鴻志,吳中如. 國(guó)外土壩潰壩模擬綜述[J].長(zhǎng)江科學(xué)院院報(bào),2003，20(2):26-29.

[5] 陳生水, 方緒順, 鐘啟明, 等. 土石壩漫頂潰壩過程離心模型試驗(yàn)與數(shù)值模擬[J]. 巖土工程學(xué)報(bào), 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] 崔 丹,槐文信,姜治兵. 潰壩洪水演進(jìn)的數(shù)值模擬[J]. 華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012，(7):34-37.

[8] 夏軍強(qiáng),王光謙,LIN Bin-liang，等.復(fù)雜邊界及實(shí)際地形上潰壩洪水流動(dòng)過程模擬[J].水科學(xué)進(jìn)展,2010,21(3):289-298.

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

[10] 岳志遠(yuǎn),曹志先,閆 軍.滑坡體潰決洪水?dāng)?shù)學(xué)模型研究[J].水動(dòng)力學(xué)研究與進(jìn)展,2008,23(5):492-500.

[11] 李相南,陳祖煜. 兩種潰壩模型在唐家山堰塞湖潰決模擬中的對(duì)比[J].水利水電科技進(jìn)展, 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.