張小雅,任春平,楊 帆

(太原理工大學(xué) 水利科學(xué)與工程學(xué)院,太原 030024)

0 引 言

在中小河流上修筑液壓壩,不僅可以有效改善河流生態(tài),城市河段的景觀,進(jìn)而實(shí)現(xiàn)水系重建、河流生態(tài)修復(fù),而且可以提高水資源利用效率[1-3]。低水頭液壓壩結(jié)構(gòu)簡(jiǎn)單,施工方便,可調(diào)控度較高,在中小河流上得到較為廣泛的應(yīng)用。但同時(shí),這些低水頭壩也會(huì)改變河流水沙過程,特別是在梯級(jí)壩群河段,會(huì)對(duì)河道沖淤特性產(chǎn)生重要影響[4]。

汾河中游生態(tài)修復(fù)核心段始于汾河二壩,終于文峪河入汾河口,全長(zhǎng)82.2 km,共計(jì)修筑15座液壓壩,相鄰壩間距約5 km,每座液壓壩由許多能自由升降的面板組成,可實(shí)現(xiàn)升壩攔水、塌壩行洪的目的,每個(gè)面板長(zhǎng)6 m,高3～5 m,一般可將每座液壓壩分成主槽壩和邊灘壩,遇到大洪水時(shí)可全部塌壩泄洪,小洪水時(shí)可以根據(jù)需要調(diào)整壩的開啟。圖1給出了汾河中游某一座液壓壩不同運(yùn)行方式過流場(chǎng)景(圖1(a)—圖1(c))及液壓壩下游主槽擺動(dòng)對(duì)堤防護(hù)坡的侵蝕(圖1(d))。由圖1可見,通過調(diào)控液壓壩(主槽段、邊灘段)不同運(yùn)行狀態(tài),可在一定程度上實(shí)現(xiàn)對(duì)水流過流狀態(tài)的調(diào)控,但不同運(yùn)行狀態(tài)對(duì)河道沖淤特性有怎樣的調(diào)控影響,還有待進(jìn)一步深入研究。

圖1 液壓壩不同運(yùn)行情況過流場(chǎng)景及液壓壩下游主槽擺動(dòng)對(duì)堤防護(hù)坡的侵蝕

目前,已有的研究大多關(guān)注高水頭壩對(duì)河道沖淤特性的影響,關(guān)于可調(diào)控低水頭梯級(jí)壩群的研究還鮮有報(bào)道。河道中修筑擋水建筑物后,會(huì)改變水動(dòng)力條件,進(jìn)而影響泥沙的輸移[5],對(duì)該方面的研究多采用二維數(shù)值模擬方法[6-9],數(shù)值模型涉及水動(dòng)力過程[10-13]和泥沙輸移過程[14-18]。已有學(xué)者就汾河中游梯級(jí)液壓壩群對(duì)洪水演進(jìn)及床面形態(tài)影響進(jìn)行了初步研究。任春平等[19]基于Delft3D數(shù)值研究了汾河中游梯級(jí)液壓壩群不同運(yùn)行方案對(duì)洪水演進(jìn)的影響。Ni等[20]基于二維水動(dòng)力-泥沙-床面形態(tài)耦合模型研究了不同運(yùn)行方案下的床面形態(tài)演變,但僅僅考慮了黏性沙。針對(duì)上述問題,本文以汾河中游二壩到義棠段為研究河段,該河段包含15座梯級(jí)液壓壩,基于Delft3D FM軟件構(gòu)建二維水沙數(shù)學(xué)模型,同時(shí)考慮黏性沙和非黏性沙,研究在不同洪水頻率下,液壓壩群不同運(yùn)行方案對(duì)河道沖淤變化的影響。本研究將為汾河中游液壓壩群的運(yùn)行提供科學(xué)依據(jù),進(jìn)而為該段河流生態(tài)修復(fù)的成功實(shí)施奠定基礎(chǔ)。

1 Delft3D FM模型構(gòu)建及驗(yàn)證

本文采用1996年、2016年、1995年和1998年洪水過程下二壩和義棠水文站觀測(cè)所得水沙數(shù)據(jù),基于Delft3D FM軟件構(gòu)建二維水動(dòng)力-泥沙耦合數(shù)學(xué)模型,對(duì)汾河二壩—義棠段的水沙輸移過程進(jìn)行數(shù)值模擬及驗(yàn)證。

1.1 控制方程

Delft3D FM軟件包括水動(dòng)力、波浪、泥沙和水質(zhì)等模塊,可應(yīng)用于海岸、河口和河流等區(qū)域的工程研究。本文通過將D-FLOW和D-MOR模塊耦合構(gòu)建二維水沙數(shù)學(xué)模型[21]。

模型采用σ坐標(biāo)系,二維水動(dòng)力模型控制方程如下。

(1)連續(xù)方程為

(1)

ξ方向和η方向的動(dòng)量方程分別為：

(2)

(3)

式中：w為垂向流速;σ為比例垂直坐標(biāo);ρ0為水的密度;Pξ和Pη分別為ξ和η方向上的靜水壓力梯度;Fξ和Fη分別為ξ和η方向上由水平雷諾應(yīng)力引起的不平衡;Mξ和Mη分別為ξ和η方向上由外部導(dǎo)致的動(dòng)量;νV為垂直渦黏系數(shù)。

(2)泥沙輸移方程。黏性沙輸沙率的計(jì)算采用Partheniades-Krone公式,即：

(4)

(5)

(6)

式中：E和D分別為侵蝕通量和沉積通量;M為侵蝕系數(shù)(kg/m2/s);ρs為泥沙密度(kg/m3);τ為床面剪切應(yīng)力;τcr,e和τcr,d分別為臨界侵蝕和臨界沉積床面剪切應(yīng)力(Pa);ωs為泥沙沉降速度(m/s);c為含沙量;cb為接近底部計(jì)算層的平均泥沙濃度(kg/m3);Δzb為近底層泥沙厚度(m)。

非黏性沙輸沙率的計(jì)算采用Engelund-Hansen (1967)全沙公式,即

(7)

式中：Sb和Ss,eq分別為推移質(zhì)和懸移質(zhì)輸沙率;α為校準(zhǔn)系數(shù);q為流速;Δ為相對(duì)密度,Δ=(ρs-ρ0)/ρ0,ρs為泥沙密度(kg/m3);C為謝才系數(shù);D50為泥沙中值粒徑(mm);g為重力加速度。

泥沙輸運(yùn)基本方程為

(8)

床面高程變化方程為

(9)

式中：c為含沙量;εs,x和εs,y分別為x、y方向的泥沙擴(kuò)散系數(shù);zb為床面高程;Sx、Sy分別為x、y方向的輸沙分量。

1.2 模型網(wǎng)格劃分

模型計(jì)算域?yàn)榉诤佣巍x棠段,全長(zhǎng)82.2 km,河寬100～400 m,平均縱向坡度為0.03%,本文考慮二壩—義棠段間匯入的昌源河、龍鳳河、磁窯河和文峪河4條支流,沿程共設(shè)15座液壓壩,見圖2。

圖2 研究區(qū)域及15座液壓壩平面位置

圖3 計(jì)算網(wǎng)格

對(duì)研究區(qū)域進(jìn)行非結(jié)構(gòu)化三角形網(wǎng)格劃分,劃分網(wǎng)格尺寸為50 m,對(duì)河寬較窄處進(jìn)行局部加密,網(wǎng)格數(shù)量為18 855個(gè),網(wǎng)格節(jié)點(diǎn)的夾角余弦值<0.02,滿足正交性要求,見圖3。

1.3 邊界條件及參數(shù)設(shè)置

本文將二壩作為上游邊界,義棠作為下游邊界。對(duì)于水動(dòng)力場(chǎng),上游設(shè)為流量-時(shí)間序列,下游設(shè)為水位-時(shí)間序列,支流邊界采用流量-時(shí)間序列。對(duì)于泥沙沖淤過程,上下游均設(shè)為含沙量-時(shí)間序列。其中1996年和2016年洪水重現(xiàn)期分別為10、2 a一遇。1995年和1998年洪水過程為常發(fā)洪水,1996年二壩和義棠2個(gè)斷面含沙量較大,1998年斷面含沙量介于1995年和1996年之間,2016年由于汾河二庫的運(yùn)行,斷面含沙量極小,見圖4。

圖4 二壩、義棠斷面含沙量

本文在數(shù)值模擬過程中涉及泥沙、河道糙率及模擬時(shí)間步長(zhǎng)相關(guān)參數(shù),下文給出具體參數(shù)設(shè)置。模擬河段主要由細(xì)沙和輕沙壤土組成,來洪量較大的情況下,泥沙以推移質(zhì)和懸移質(zhì)形式運(yùn)動(dòng),因此采用黏性沙和非黏性沙的泥沙組分[22-23],具體參數(shù)見表1。

表1 泥沙模型參數(shù)

以上下游水位和流量過程實(shí)測(cè)資料為基礎(chǔ),率定主槽糙率為0.007,邊灘糙率為0.03。D-FLOW模塊使用顯式平流計(jì)算原理,根據(jù)CFL(Courant-Friedrichs-Levy)計(jì)算準(zhǔn)則自動(dòng)設(shè)置動(dòng)態(tài)時(shí)間步長(zhǎng)限制,每個(gè)時(shí)間步長(zhǎng)的最大庫朗數(shù)(Max Courant nr)設(shè)為0.7[21]。

1.4 模型驗(yàn)證

本文采用1996年和2016年洪水過程進(jìn)行水動(dòng)力模型驗(yàn)證,驗(yàn)證結(jié)果見圖5。由圖5可知,二壩斷面水位和義棠斷面流量的峰值和位置都與實(shí)測(cè)數(shù)據(jù)較為吻合。由于未考慮植被對(duì)局部糙率的影響,部分時(shí)段水位、流量的模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)有差異。總體來看,整個(gè)洪水過程的水位和流量的計(jì)算結(jié)果與實(shí)測(cè)數(shù)據(jù)吻合較好,說明構(gòu)建的二維水動(dòng)力模型精度較高,可以用于后續(xù)水沙輸移模擬。

圖5 1996年和2016年洪水驗(yàn)證結(jié)果

圖6 驗(yàn)證斷面位置

本文采用斷面含沙量和斷面形態(tài)進(jìn)行泥沙模型驗(yàn)證,各斷面位置見圖6,二壩—義棠方向(右上—左下)共18個(gè)斷面,依次為CS-A、CS-01—CS-16、CS-B。將2016年洪水過程下二壩下游斷面CS-A和義棠上游斷面CS-B含沙量的計(jì)算結(jié)果與實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比,如圖7所示。將2014年洪水過程下河道16個(gè)斷面形態(tài)變化的計(jì)算結(jié)果與實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比,因篇幅有限僅展示部分?jǐn)嗝娴尿?yàn)證圖,見圖8。

圖7 斷面含沙量實(shí)測(cè)與模擬結(jié)果比較

圖8 斷面形態(tài)實(shí)測(cè)結(jié)果與模擬結(jié)果比較

根據(jù)統(tǒng)計(jì)學(xué)方法[24]對(duì)模擬結(jié)果進(jìn)行評(píng)價(jià),計(jì)算公式為

(10)

計(jì)算得出,斷面含沙量的模型效率系數(shù)分別為0.87和0.63,評(píng)價(jià)結(jié)果為極好和較好,斷面形態(tài)的模型效率系數(shù)為0.45～0.73,評(píng)價(jià)結(jié)果為好至較好。由于未考慮植被對(duì)局部糙率的影響,部分時(shí)間段的斷面含沙量模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)有差異。軟件對(duì)地形數(shù)據(jù)的差值處理導(dǎo)致斷面形態(tài)較為均勻,部分位置的模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)有差異?？傮w來看,本文所構(gòu)建的二維水沙數(shù)學(xué)模型模擬精度較高,能夠較好地模擬出河道床面形態(tài)變化過程,可以用于液壓壩群對(duì)河道沖淤變化影響的數(shù)值研究。

1.5 液壓壩工況設(shè)置

本文研究的河段在15個(gè)斷面處筑有液壓壩,每個(gè)斷面處液壓壩可分為主槽壩和邊灘壩。為研究液壓壩群不同運(yùn)行方案對(duì)河道沖淤變化的影響規(guī)律,本文針對(duì)4種運(yùn)行方案進(jìn)行了研究,見表2。

表2 液壓壩運(yùn)行方案

2 液壓壩群對(duì)河道沖淤變化的影響

本節(jié)運(yùn)用上述二維水沙數(shù)學(xué)模型,對(duì)4場(chǎng)洪水情況下液壓壩群的不同運(yùn)行方案進(jìn)行模擬,研究液壓壩群對(duì)河道沖淤變化的影響。

2.1 水動(dòng)力場(chǎng)及床面高程變化

本文對(duì)14#液壓壩進(jìn)行水動(dòng)力及床面高程變化分析。表3給出了壩前水深分布,在1996年、2016年、1995年和1998年4場(chǎng)洪水情況下,無液壓壩運(yùn)行時(shí)(方案1),壩前水深分別約為1.0～1.4、0.7～1.2、0.6～1.1、0.5～0.8 m,壩后局部流速最大約為2.5、2.0、2.0、1.5 m/s。液壓壩全部運(yùn)行時(shí)(方案2)對(duì)水流起到攔蓄作用,壩前水深分別增加了約0.8～1.4、0.5～1.1、0.9～1.0、0.8～1.0 m,形成過壩溢流,壩后局部流速最大約為6.0、3.0、4.0、2.0 m/s,為無壩運(yùn)行時(shí)的1.3～2.4倍。僅有主槽壩運(yùn)行時(shí)(方案3),壩前水深分別增加了約0.3～0.4、0.1～0.2、0.2～0.4、0.1～0.2 m,壩后局部流速最大約為3.0、2.5、2.5、2.0 m/s,為無壩運(yùn)行時(shí)的1.2～1.3倍。僅有邊灘壩運(yùn)行時(shí)(方案4),14#壩的壩前水深分別增加了約0.1～0.2、0～0.1、0～0.1、0～0.1 m,水流流向主槽,邊灘壩后流速明顯減小,約為0.1～0.2 m/s,主槽流速與無壩運(yùn)行時(shí)相近。

表3 14#壩的壩前水深、壩后流速分布

表4 15座液壓壩前床面高程變化ΔZ

在1996年、2016年、1995年和1998年4場(chǎng)洪水情況下,由于液壓壩攔水?dāng)r沙的作用,壩前產(chǎn)生泥沙淤積,壩后局部流速較大,壩后形成床面沖刷(壩上下游各100 m處)。無液壓壩運(yùn)行時(shí),泥沙沖淤范圍約為-0.5～0.5、-0.5～0.25、-0.4～0.2、-0.4～0.2 m。液壓壩全部運(yùn)行時(shí),由于液壓壩群聯(lián)合作用,處于下游的14#壩前泥沙淤積厚度較小,泥沙沖淤范圍約為-0.5～0.6、-0.5～0.05,-0.5～0.1、-0.5～0.2 m。僅有主槽壩運(yùn)行時(shí),沖淤范圍約為-0.5～0.5、-0.5～0.2、-0.5～0.1、-0.5～0.2 m。僅有邊灘壩運(yùn)行時(shí),壩后由于流速較小,邊灘處產(chǎn)生淤積最大值約為0.5、0.3、0.15、0.2 m,主槽沖淤范圍與無壩運(yùn)行時(shí)相同。

2.2 15座壩液壓壩前泥沙沖淤變化

本節(jié)選取每座液壓壩上游100 m處的斷面作為參考,分析不同液壓壩運(yùn)行方案下15座液壓壩前的床面泥沙沖淤高程變化范圍(見表4),研究液壓壩群對(duì)河道沖淤變化的影響規(guī)律。

由表4可知,在1996年、1995年、1998年3場(chǎng)洪水過程下,相對(duì)于方案1(無液壓壩運(yùn)行)泥沙沖淤變化情況,方案2(液壓壩全部運(yùn)行)時(shí)1#—9#壩前泥沙淤積厚度增加0.05～1.07 m,沖刷范圍減小0.04～0.40 m,而10#—15#液壓壩處在河道下游,流速和含沙量相對(duì)減小,泥沙沖淤范圍減小0.05～0.60 m;方案3(僅有主槽壩運(yùn)行)中壩前泥沙淤積無明顯變化規(guī)律,沖刷范圍相對(duì)增加0.05～0.30 m;方案4(僅有邊灘壩運(yùn)行)中泥沙沖淤變化范圍與方案1近似。2016年洪水過程由于含沙量較小,液壓壩全部運(yùn)行時(shí),1#—3#壩前淤積厚度增加0.05～0.10 m,僅有主槽壩運(yùn)行時(shí)床面整體沖刷范圍增加0.1～0.4 m,僅有邊灘壩運(yùn)行時(shí)整體泥沙沖淤范圍與方案1相近。

2.3 沖淤總量分析

下面分析4場(chǎng)洪水情況下整個(gè)模擬河段內(nèi)總的泥沙淤積量,公式為[22]

ΔV=V二壩-V義棠。

(11)

圖9 4場(chǎng)洪水過程在不同運(yùn)行方案下的泥沙淤積

式中：ΔV為河道淤積總量;V二壩為二壩斷面入口泥沙量(萬m3);V義棠為義棠斷面出口泥沙量(萬m3)。

圖9給出了4場(chǎng)洪水情況下的泥沙淤積總量,可知河道中泥沙總體呈現(xiàn)淤積狀態(tài)。1996年洪水過程下,方案1河道泥沙淤積量約為242.3萬m3,與無壩運(yùn)行時(shí)相比,方案2—方案4泥沙淤積量分別多了約81.2萬、10.2萬、3.1萬m3。2016年洪水過程下,由于汾河二庫的運(yùn)行,河道內(nèi)含沙量較小,方案1河道泥沙淤積量約為4.9萬m3,與無壩運(yùn)行時(shí)相比,方案2—方案4分別多了約0.4萬、0.1萬、0.05萬m3。1995年洪水過程下,方案1河道泥沙淤積量約為62.3萬m3,與無壩運(yùn)行時(shí)相比,方案2—方案4分別多了約5.3萬、3.1萬、1.3萬m3。1998年洪水過程下,方案1河道淤積量約為22.1萬m3,與無壩運(yùn)行時(shí)相比,方案2—方案4分別多了約8.1萬、3.3萬、1.5萬m3。可知在4場(chǎng)洪水過程下,方案2(液壓壩全部運(yùn)行)泥沙淤積量最大,約為方案1的1.08～1.36倍,方案1(無液壓壩運(yùn)行)淤積量最小,方案3和方案4(僅有主槽壩運(yùn)行和僅有邊灘壩運(yùn)行)的淤積量值介于兩者之間,分別約為方案1的1.04～1.15倍、1.02～1.07倍。

3 結(jié) 論

本文采用1996年、2016年、1995年和1998年4場(chǎng)洪水的汛期水沙數(shù)據(jù),設(shè)計(jì)了4種液壓壩運(yùn)行方案,基于Delft3D FM軟件構(gòu)建二維水動(dòng)力-泥沙耦合數(shù)學(xué)模型,研究了汾河中游生態(tài)修復(fù)核心區(qū)15座液壓壩作用下的河道泥沙沖淤變化過程。

(1)液壓壩的運(yùn)行導(dǎo)致壩前后水動(dòng)力場(chǎng)和床面高程變化范圍增大。在4種液壓壩運(yùn)行方案中,液壓壩全部運(yùn)行時(shí)對(duì)水動(dòng)力場(chǎng)和床面高程變化影響最大,僅有邊灘壩運(yùn)行時(shí)的影響最小,汛期在滿足防洪要求時(shí),可考慮不對(duì)邊灘壩進(jìn)行塌壩處理。

(2)隨著洪水過程的含沙量增大,液壓壩群對(duì)河道沖淤變化影響范圍隨之增加,僅有邊灘壩運(yùn)行時(shí)影響較小。4場(chǎng)洪水情況下方案1—方案4中15座液壓壩前泥沙沖淤范圍分別約為-0.5～1.4、-0.3～1.9、-0.5～1.6、-0.5～1.5 m。

(3)河道泥沙總體呈現(xiàn)淤積狀態(tài),液壓壩全部運(yùn)行時(shí)泥沙淤積量約為無壩運(yùn)行時(shí)的1.08～1.36倍,僅有主槽壩運(yùn)行的淤積量約為無壩運(yùn)行時(shí)的1.04～1.15倍,僅有邊灘壩運(yùn)行時(shí)淤積量為無壩運(yùn)行時(shí)的1.02～1.07倍。