張幫穩(wěn),吳保生,章若茵

(清華大學(xué)水沙科學(xué)與水利水電工程國家重點實驗室,北京 100084)

三峽水庫泥沙淤積問題不僅關(guān)系到水庫的使用壽命和功能發(fā)揮,而且對水庫及下游的生態(tài)環(huán)境具有重要的影響[1]。近年來長江上游梯級水庫的建設(shè)導(dǎo)致三峽入庫徑流量和泥沙量發(fā)生變化,嚴重影響了三峽水庫綜合效益的發(fā)揮,為此,三峽水庫在“蓄清排渾”運用方式的基礎(chǔ)上,開展了中小洪水調(diào)度和汛末提前蓄水調(diào)度,目的是提高發(fā)電和航運效益,減輕壩下游的防洪壓力[2]。但汛期中小洪水調(diào)度和汛末提前蓄水會抬高庫區(qū)平均水位,降低庫區(qū)流速,加重庫區(qū)泥沙淤積。為緩解中小洪水調(diào)度及汛末提前蓄水調(diào)度帶來的泥沙淤積影響,三峽水庫于2012年和2013年汛期利用洪峰和沙峰的異步運動特性開展了沙峰調(diào)度試驗,取得了較好的排沙效果[3- 4]。

研究者基于實測資料對三峽庫區(qū)洪峰和沙峰的運動特性進行了分析,發(fā)現(xiàn)三峽水庫蓄水后,庫水位抬高,庫區(qū)水流流速減慢,沙峰輸移時間較蓄水前大幅增加,沙峰多滯后于洪峰,且沙峰滯后于洪峰的時間沿程逐漸增加,有利于進行水庫的沙峰調(diào)度[5- 7]。但目前關(guān)于洪峰和沙峰異步運動特性多以實測資料進行分析,基于試驗研究及數(shù)值模擬開展相對較少,對洪峰和沙峰異步運動特性缺乏較為深入的認識。利用數(shù)值模擬方法對不同情況下洪峰和沙峰運動特性進行系統(tǒng)分析,可為應(yīng)對未來愈加復(fù)雜的水沙變化和地形變化情勢下水庫的沙峰排沙調(diào)度提供參考。水沙數(shù)學(xué)模型在研究水沙輸移、泥沙淤積及河床演變規(guī)律中有著重要的作用?；跀嗝嫫骄虼瓜蚱骄囊痪S和二維數(shù)值模型[8- 9]能夠?qū)こ趟硢栴}進行一定程度的研究,但是泥沙淤積、沖刷和輸移是一個非常復(fù)雜的水動力現(xiàn)象,特別是在復(fù)雜的自然地形,例如河漫灘、彎曲河段和沙波地形等產(chǎn)生二次流和流體分離等三維水動力特性,嚴重影響了數(shù)值模擬結(jié)果的準確性。隨著計算機科學(xué)技術(shù)和高性能計算機的發(fā)展,三維模型[10- 11]逐漸應(yīng)用到工程實際當中,目前大多數(shù)采用三維數(shù)值模擬方法針對局部河段或庫區(qū)的泥沙輸運和淤積進行研究,但對洪水傳播過程中洪峰和沙峰異步運動特性研究較少。主要原因是洪峰和沙峰的運動特性與研究河流或水庫的尺度(時間尺度和空間尺度)有關(guān),較短的距離和較短的時間難以顯示出洪峰和沙峰運動規(guī)律。三峽庫區(qū)地形較為復(fù)雜,總長超660 km,對此進行長距離洪水演進的三維數(shù)值模擬具有較大的難度。

掌握洪峰與沙峰異步運動特性、開展入庫泥沙實時監(jiān)測與預(yù)報,是進行沙峰排沙調(diào)度的基礎(chǔ)。為進一步提高沙峰預(yù)報精度,掌握更精準的沙峰調(diào)度時機、制訂更科學(xué)合理的沙峰排沙調(diào)度方案,本文采用三維水沙數(shù)值模型SCHISM(Semi- implicit Cross- scale Hydroscience Integrated System Model)對三峽庫區(qū)汛期洪水傳播過程中洪峰和沙峰運動進行模擬,初步分析三峽水庫不同蓄水位下洪峰和沙峰在萬縣—三峽大壩庫區(qū)的異步運動特性,可為三峽水庫的沙峰調(diào)度提供參考依據(jù)。

1 數(shù)值模擬方法

1.1 水流的基本方程

SCHISM水動力學(xué)模型[12- 13]的控制方程采用基于Reynolds時均N- S方程,滿足靜壓假定和Boussinesq渦黏性假定。在笛卡爾坐標系下,對于不可壓縮流體N- S方程的連續(xù)性方程為

(1)

動量守恒方程:

(2)

(3)

式中:x和y分別表示笛卡爾水平坐標,z為垂向坐標,向上為正；u,v,w分別表示3個方向的流速；t為時間；f為柯氏力系數(shù)；η為自由水面；zb為河床底高程；ρ0和ρ分別表示參考密度和混合流體的密度；g為重力加速度；Kmh和Kmv分別為水平與垂直渦黏性系數(shù),其中垂向渦黏性系數(shù)根據(jù)紊流模型進行封閉,水平渦黏性系數(shù)采用常數(shù)化處理；pa為自由水面大氣壓強。

自由水面采用水位函數(shù)法處理,對連續(xù)方程(1)沿水深方向積分,可得自由水面方程

(4)

模型在河床底部的動力學(xué)邊界條件由底床摩擦剪應(yīng)力和底層水體的雷諾應(yīng)力平衡給出

(5)

式中:tbx和tby分別為床面的x、y方向摩擦剪應(yīng)力。

對于垂向紊動渦黏性系數(shù)Kmv,利用紊流閉合模型GLS(Generic Length Scale)[14]進行求解,包括紊動能k方程和通用紊動長度ψ方程。GLS紊流模型的控制方程為:

(6)

(7)

(8)

(9)

式中:k為紊動能；ψ為通用紊動長度參數(shù)；μ為鹽度、溫度等物質(zhì)的垂向擴散系數(shù)；υk和υy分別表示紊動能和通用紊動長度的垂向擴散系數(shù)；ε為紊動耗散項；Fw為壁函數(shù),在k- ε模式中為1；σk、σψ、cψ1、cψ2和cψ3為模型系數(shù)；M2和N2分別表示由于剪切變形和密度分層而引起的紊動能產(chǎn)生項,通用紊動長度ψ和紊動耗散項ε作為紊流模型中的關(guān)鍵參量,其表達式如下:

ψ=(cμ)mknlp,ε=(cμ)mknlp

(10)

式中:l為紊動摻混長度；cμ為常數(shù)0.3。在GLS模型k-ε模式雙方程紊流模式的參數(shù)取值見表1。

表1 GLS模型k- ε模式雙方程紊流參數(shù)值Table 1Parameter values of the turbulence in k- ε equation based on GLS model

1.2 懸移質(zhì)泥沙輸移基本方程

在懸移質(zhì)泥沙動力學(xué)中,對流擴散理論將泥沙視為單一的連續(xù)介質(zhì),并假設(shè)懸移質(zhì)泥沙運動與水流運動在垂直方向上存在速度差,且等于泥沙顆粒的沉降速度?；趯α? 擴散理論的三維懸移質(zhì)泥沙輸運方程表達式為

(11)

式中:Ci為i組含沙量；ws,i為i組泥沙顆粒的沉降速度；Ksv為泥沙垂向擴散系數(shù),通常假定與水流紊動黏性系數(shù)呈倍數(shù)關(guān)系,可通過紊流模型求解或采用經(jīng)驗關(guān)系估計，Ksv=Kmv/σs,σs為Schmidt數(shù),通常取值在0.6～1.2之間；Ksh為泥沙水平擴散系數(shù),考慮到水平擴散的量級遠小于垂向擴散,通常忽略不計；Dq為泥沙的沉積通量;Eq為泥沙的沖刷通量。三峽水庫中細顆粒泥沙存在絮凝現(xiàn)象,在泥沙運動力學(xué)中,采用絮凝因數(shù)F反映絮凝對泥沙顆粒的沉速的影響[15],其表達式為

F=ωs,i/ω0

(12)

式中:ω0為泥沙顆粒靜水中的沉速。廣泛采用的絮凝因數(shù)F=αdβ,其中,α取值為0.013[16]或0.005 5[17],β取值為-1.9[15]或-1.02[16]。張地繼等[16]采用α=0.005 5,β=-1.02研究了萬縣—廟河庫區(qū)沙峰的衰減規(guī)律,本研究中采用相同的絮凝因數(shù)系數(shù)。

1.3 河床沖刷演變基本方程

根據(jù)泥沙質(zhì)量守恒方程,懸移質(zhì)泥沙輸移過程中可將床面邊界條件視為床面附近泥沙通量的處理,包括床面泥沙的沉積通量Db和沖刷通量Eb,其表達式分別為:

Db=ωs,ic1,i

(13)

Eb=E0,i(1-qi)(τf/τcr,i-1)τf>τcr,i

(14)

式中:c1,i為數(shù)值模擬中最底部一層網(wǎng)格的i組泥沙濃度；E0,i為經(jīng)驗沖刷率系數(shù),取決于局部床面泥沙顆粒條件,取值大小范圍為10-4～10-2kg/(m2·s-1)；qi為床面層i組泥沙體積分數(shù)；τcr,i為i組泥沙顆粒臨界起動剪切應(yīng)力；τf為水流底部的剪切應(yīng)力。河床泥沙的沖淤和懸移質(zhì)泥沙的凈輸移導(dǎo)致河床表面的地形變化,其公式為

(15)

式中:Δh為床面高程變化量;ρs為泥沙顆粒的密度。

2 數(shù)值模型的建立及驗證

為了研究三峽庫區(qū)汛期洪水傳播過程中洪峰和沙峰異步運動特性,本文選擇萬縣站—三峽大壩的庫區(qū)段作為研究對象(圖1)。若把寸灘或清溪場作為入口邊界條件,則河段過長,計算耗時過大,計算效率極低。萬縣至大壩長280 km,不僅長度適中,并且也是目前少有的長距離、大范圍實際水庫的三維數(shù)值模擬,難度較高；另一方面原因,萬縣—三峽大壩庫區(qū)之間僅有靠近大壩廟河站可以測量流量和含沙量,長度不足以模擬洪峰和沙峰的傳播特性,因此只有萬縣站滿足本文的研究要求,可作為入口條件。模擬中采用的資料分別來源于萬縣站、奉節(jié)站、巫山站、廟河站和茅坪站,如圖1(a)所示,其中水文站有萬縣站和廟河站,具有日平均水位、流量和含沙量信息,水位站有奉節(jié)站、巫山站和茅坪站。

圖1 研究河段的地形高程、水文站和網(wǎng)格設(shè)置Fig.1 Topographic elevation,hydrological station and grid setting of the study reach

2.1 模型網(wǎng)格設(shè)置

萬縣—三峽大壩的地形較為復(fù)雜,寬闊和狹窄的河段交替變化,水平方向上進行網(wǎng)格設(shè)置時主要混合使用2種網(wǎng)格類型,在寬闊河段混合采用四邊形和三角形網(wǎng)格,主河道采用四邊形網(wǎng)格,岸灘采用三角形網(wǎng)格,在狹窄河段采用三角形網(wǎng)格,網(wǎng)格尺度約為20～23 m,共得到495 662個網(wǎng)格節(jié)點和856 056個網(wǎng)格單元。計算初始地形條件由2011年實測地形插值得到,插值后的局部地形如圖1(b)所示。此外,為了更好地模擬河道至壩前水深變化幅度大的特點并且提高計算效率,垂向上采用分層LSC2坐標[17],最大水深處可達36層,最淺處為16層,平均約為19層,每層厚度約5～6 m。圖1(c)為局部河段橫斷面的垂向網(wǎng)格示意圖。為了數(shù)值模擬的穩(wěn)定性,時間步長設(shè)置30 s。整個計算模型在清華大學(xué)高性能計算集群上采用280個核進行計算,模擬時間為8 d。

2.2 模型邊界條件

本文以三峽庫區(qū)2013年汛期7月1日—8月1日作為模擬時間段,萬縣站的流量和含沙量過程作為模型的入口邊界條件,茅坪站的水位過程作為出口邊界條件(圖2)。三峽庫區(qū)泥沙主要以懸移質(zhì)輸移為主,萬縣站7月實測的泥沙粒徑分別為0.002 mm、0.004 mm、0.008 mm、0.016 mm、0.032 mm和0.064 mm,對應(yīng)的泥沙級配分別為10.3%、12.3%、20.9%、26.0%、18.5%和12.0%,考慮到泥沙0.002 mm和0.004 mm 的粒徑很小,統(tǒng)一歸為0.003 mm的泥沙進行計算,床沙依據(jù)大斷面實測的泥沙級配進行插值設(shè)置,并且選擇與懸移質(zhì)同樣的5組粒徑泥沙。數(shù)值模擬泥沙輸移過程中暫時沒有考慮溫度和鹽度的影響。

圖2 數(shù)值模型邊界條件Fig.2 Boundary conditions in numerical model

2.3 數(shù)值模型驗證

為驗證數(shù)值模型模擬洪水傳播和泥沙輸移過程的正確性,分別對比分析了實測與模擬的奉節(jié)站和廟河站的水位,以及廟河站的流量、平均流速和含沙量。

圖3對比了奉節(jié)站和廟河站水位的模擬值和實測值,可以看出數(shù)值模擬的結(jié)果和實測值吻合較好,能夠較好地反映洪水傳播過程中水位的變化。圖4對比了廟河站斷面平均的流量、流速和含沙量的模擬值和實測值,可以看出洪水流量模擬值的變化和實測值的變化過程基本一致,局部的差別可能是由于萬縣—三峽大壩之間支流入?yún)R和復(fù)雜的局部地形糙率的影響,經(jīng)過分析2013年7月支流平均流量的總和占萬縣站流量的1%,本文中暫時沒有考慮支流徑流對主河道洪水傳播的影響。平均含沙量的模擬值與廟河站實測時刻的含沙量點較為接近,和2013年水文年鑒中日均含沙量有局部的差別。這是因為水文年鑒中日均含沙量是基于實測某一時刻的含沙量回歸插值得到的[18],本身存在一定的誤差,也可能是限于萬縣—三峽大壩之間實測資料的限制,模擬區(qū)域局部的沖淤不能很好地反映出來。

圖3 水位模擬值和實測值的對比結(jié)果Fig.3 Comparison of water level between simulated and measured results

圖4 廟河站平均流量、流速和平均含沙量的模擬值和實測值對比Fig.4 Comparison of average flow discharge,velocity and sediment concentration between simulated and measured results at Miaohe station

圖5給出了廟河站斷面流速和含沙量分布的瞬時模擬結(jié)果。從圖中可以看出流速和含沙量分布與斷面的幾何形態(tài)有關(guān),斷面含沙量的分布不均勻,呈現(xiàn)出分層的現(xiàn)象。圖6給出了廟河站垂向流速和含沙量的測量值和模擬值的對比結(jié)果,測點的位置見圖5,測點之間間隔120 m。從圖中可以看出垂向流速的數(shù)值模擬結(jié)果和測量值吻合較好,垂向含沙量的模擬結(jié)果和測量值存在一定的誤差,數(shù)值模擬斷面兩側(cè)的測點底部含沙量偏小,可能與數(shù)值模型沒有考慮推移質(zhì)的泥沙運動和斷面幾何形態(tài)有關(guān)。三維水沙數(shù)值模型能夠進行斷面垂向流速和含沙量的分析,相對于一維和二維數(shù)值模擬的精度更高。

圖5 廟河站斷面流速和含沙量模擬結(jié)果的瞬時分布Fig.5 Instantaneous distribution of simulated flow velocity and sediment concentration at Miaohe station

圖6 廟河站斷面垂向流速和含沙量模擬值和實測值的對比結(jié)果Fig.6 Comparison of vertical flow velocity and sediment concentration between simulated and measured results at Miaohe station

3 三峽庫區(qū)壩前不同蓄水位下洪峰和沙峰異步運動特性

三峽工程于2003年6月進入圍堰蓄水期,壩前水位汛期按135 m、枯季139 m運行；2006年汛后初期蓄水后,壩前水位按汛期144 m、枯季156 m運行；自2008年汛末三峽水庫進行175 m試驗性蓄水以來,工程進入175 m試驗性蓄水期。為了更好地分析三峽水庫蓄水以來汛期不同蓄水位下洪峰和沙峰異步運動特性的規(guī)律,本文設(shè)置3個研究方案。方案2為2013年7月實測洪水期間壩前水位；方案1和方案3分別為2013年7月實測洪水期間壩前水位減去和加上10 m。圖7給出了3種研究方案的壩前水位變化?；跀?shù)值模擬的結(jié)果分別提取萬縣—三峽大壩之間重要水文站的流量和含沙量隨時間的變化,用以分析洪峰和沙峰沿程的異步運動特性。

圖7 不同方案的壩前水位變化Fig.7 Temporal variation in water level in front of the dam in different schemes

圖8給出了整個研究區(qū)域三維流場瞬時的模擬結(jié)果,從圖中可以看出上游庫區(qū)的水流流速較大,壩前段水流流速較小,主要由于壩前段水深較大。

圖8 研究區(qū)域三維流場瞬時的模擬結(jié)果Fig.8 Instantaneous results of the three- dimensional simulated flow field in the study zone

圖9分別給出了不同方案下重要水文站的流量和含沙量隨時間變化的模擬結(jié)果對比。從沿程各水文站流量的模擬結(jié)果可以看出,蓄水位的變化對洪峰傳播時間的影響較?。粡暮沉磕M結(jié)果可以看出,蓄水位的變化對沙峰傳播時間和大小的影響較大,隨著蓄水位增加,沙峰滯后洪峰的時間逐漸增加,其實質(zhì)是洪水在河道型水庫向下游傳播過程中,隨著水深的增加,流速降低導(dǎo)致沙峰傳播越來越慢,沙峰滯后洪峰的時間也越來越大；同時由于水流流速減小,水流的挾沙能力降低,泥沙沿程不斷地落淤導(dǎo)致造成沙峰坦化,沙峰的峰值逐漸減小。從奉節(jié)站、巫山站和巴東站的含沙量隨時間變化的曲線形態(tài)可以看出在巫山站—巴東站庫區(qū)間存在局部的泥沙沖刷。

圖9 3種不同方案下沿程各水文站流量和含沙量的模擬結(jié)果對比Fig.9 Temporal variation of flow discharge and sediment concentration at each hydrographic station under different schemes

為了定量地反映壩前蓄水位的變化對三峽庫區(qū)洪峰和沙峰異步運動特性的影響,表2分別給出了不同方案下沿程各水文站的洪峰和沙峰到達時間及沙峰滯后洪峰的時間,表明壩前蓄水位的變化對洪峰傳播時間的影響較小,對沙峰傳播時間的影響較大；并且沿程各水文站沙峰滯后于洪峰的時間隨著壩前蓄水位增加越來越大,即沙峰滯后洪峰的時間隨著水流流速的降低越來越大。方案2和方案3相對于方案1壩前沙峰滯后于洪峰的時間分別增加了6.4%和16%。

表2 不同方案下沿程各水文站洪峰和沙峰到達時間及滯后時間Table 2Arrival time and lag time of flood peak and sediment peak reaching at each hydrographic station under the different schemes

4 結(jié) 論

本文采用三維水沙數(shù)值模型SCHISM對萬縣—三峽大壩280 km的庫區(qū)進行了洪水傳播和泥沙輸移的大尺度數(shù)值模擬,主要對比分析了不同壩前蓄水位下沿程各水文站洪峰和沙峰異步運動特性,結(jié)果表明:

(1) 三維數(shù)值模型SCHISM能夠較好地模擬三峽庫區(qū)長距離的洪水傳播和泥沙輸移過程,與實測值驗證結(jié)果較好。

(2) 三峽水庫壩前蓄水位的變化對洪峰傳播時間的影響不明顯,對沙峰傳播時間的影響較為顯著；壩前水位的增加導(dǎo)致水流流速減小,沙峰傳播減慢,引起庫區(qū)主要水文站沙峰滯后于洪峰的時間越來越大。

(3) 在萬縣—三峽大壩庫區(qū)洪水傳播過程中,隨著水深的增加,水流流速減小,水流挾沙能力降低,泥沙不斷的落淤導(dǎo)致沙峰峰值減小。

致謝:本研究得到了長江水利委員會水文局許全喜教高、武漢大學(xué)張為副教授、鄭珊副教授的幫助和支持,數(shù)值模型在清華大學(xué)探索100集群上計算,特此致謝！