汴西湖引水對黑崗口水庫的影響研究

王宇 秦凈凈 職保平

摘要：為了研究汴西湖建成后黑崗口水庫泥沙分布變化的相關(guān)問題，通過對汴西湖從黑崗口水庫引水?dāng)?shù)據(jù)的分析，結(jié)合黑崗口-汴西湖地形數(shù)據(jù)，依據(jù)插值和擬合數(shù)學(xué)理論，利用Flow-3D軟件進(jìn)行建模計(jì)算得到黑崗口水庫的水沙分布，計(jì)算結(jié)果表明修建汴西湖后引水閘位置處引水?dāng)嗝媸湛s，造成一定量的淤積。若不采取相應(yīng)的措施，將會造成引水效率下降，從而使黑崗口水庫和汴西湖的連接干渠倒虹吸喪失正常輸水功能。同時(shí)，綜合對比不同引水方案，得出了根據(jù)實(shí)際降水量的時(shí)間分布來調(diào)節(jié)輸水量的優(yōu)化方案。

關(guān)鍵詞：數(shù)值模擬;水沙分布;引水;汴西湖;黑崗口水庫

中圖分類號：TV145 文獻(xiàn)標(biāo)志碼：A

水庫中的泥沙沉積問題是泥沙研究中的一項(xiàng)重要內(nèi)容。很多學(xué)者利用物理模型對該問題進(jìn)行了研究，給出了諸多水庫中水沙的分布變化規(guī)律。隨著泥沙研究的發(fā)展及數(shù)值模擬方法的不斷完善，利用計(jì)算機(jī)數(shù)值模擬軟件可以更為快捷方便地對該類常見問題進(jìn)行研究。Flow-3D軟件自推出后，其基于CFD的特有VOF（Volume of Fluid）技術(shù)在數(shù)值模擬計(jì)算的準(zhǔn)確度上可取得極好的效果。張紅藝等[1-3]通過對高含沙水庫和河道的數(shù)值模擬得出了泥沙分布相關(guān)規(guī)律;假冬冬等[4]通過三維數(shù)學(xué)模型系統(tǒng)分析了水沙條件變化對河型河勢的影響及水庫泥沙沖淤分布情況。

黑崗口水庫位于河南省開封市西北部，自黃河引流蓄水。黑崗口引黃灌區(qū)調(diào)蓄水庫（又名汴西湖）工程位于開封西部新區(qū)，是銜接開封新老城區(qū)的湖泊-森林生態(tài)景觀帶。該工程從黑崗口水庫引流黃河水，經(jīng)西干渠倒虹吸地下輸水通道直達(dá)汴西湖。汴西湖自2010年開始修建，至2015年一期工程完工正式蓄水。

為了分析汴西湖部分修建完工后對引水源頭——黑崗口水庫水沙分布的影響，本文在前人研究的基礎(chǔ)之上，以黑崗口水庫和汴西湖兩者地形數(shù)據(jù)為基礎(chǔ)，建立了黑崗口水庫一汴西湖模型，同時(shí)以實(shí)測水沙資料作為計(jì)算邊界條件進(jìn)行了三維數(shù)值模擬，根據(jù)模擬結(jié)果對修建汴西湖后黑崗口水庫的水沙沉積和分布等問題進(jìn)行了分析研究，以期為水庫的運(yùn)行和管理工作提供依據(jù)。

1 模擬計(jì)算模型

1.1 控制方程

數(shù)值模型計(jì)算控制方程[5]包括連續(xù)性方程、動量方程、紊動能k方程、紊動能耗散率ε方程、挾沙力計(jì)算方程。根據(jù)計(jì)算水力學(xué)相關(guān)理論[6-8]可知連續(xù)性方程為式中：u、v、w分別為x、y、z三個(gè)方向上的流速分量;A_x、A_y、A_z分別為x、y、z三個(gè)方向的面積分?jǐn)?shù)。

動量方程為式中：t為時(shí)間;G_x、G_y、G_z分別為x、y、z三個(gè)方向的重力加速度：f_x、f_y、f_z為x、y、z三個(gè)方向的黏滯力;V_F為體積分?jǐn)?shù);ρ為流體密度;p為作用在流體微元上的壓力。

紊動能k方程為式中：k為紊動能;ε為紊動能耗散率;μ為紊動能黏性系數(shù)聲;μ_t=ρC_uk²/ε，C_y=0.0845;G_k為紊動能k的產(chǎn)生項(xiàng)，σ_k為湍動能對應(yīng)的Prandtl數(shù)，為1.39;x_i、x_j為方向變量。

紊動能耗散率ε方程為式中：σ_ε為湍流耗散率對應(yīng)的Prandtl數(shù)，為1.39;C^*_ε1、C_ε2為經(jīng)驗(yàn)常數(shù)。

挾沙力計(jì)算方程[9]為式中：S^*為沖沙質(zhì)挾沙力;S_W為沖瀉質(zhì)含沙量;S^*_w為沖瀉質(zhì)挾沙能力;s_i^*為分組挾沙力;S^*_b為床沙質(zhì)挾沙力。

在沖淤過程中，床沙級配在不斷調(diào)整，本模型床沙級配調(diào)整方程采用吳偉民、李義天模式：

P_bu=[ΔZ_i+（E_m-ΔZ）P_obi]/E_m（7）式中：P_obi、P_bi分別為時(shí)段初、時(shí)段末的床沙級配;E_m為床沙可動層厚度，其大小與河床沖淤狀態(tài)、沖淤強(qiáng)度及沖淤歷時(shí)有關(guān)，單向淤積時(shí)E_m=ΔZ（Z為河床泥沙厚度），單向沖刷時(shí)E_m的限制條件是保證床面有足夠的泥沙補(bǔ)償。

1.2 模型建立及邊界條件

在利用數(shù)值模擬方法求解控制方程時(shí)，基本思路是對控制方程進(jìn)行離散，通過求解得到的離散方程組來達(dá)到最終需要的結(jié)果[10-12]。而在空間域上離散控制方程則需要使用網(wǎng)格。目前的網(wǎng)格主要分為結(jié)構(gòu)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格兩大類，在三維建模中，常用的網(wǎng)格單元有四面體、六面體、三棱體等不同形式，在整個(gè)計(jì)算域上，網(wǎng)格通過節(jié)點(diǎn)聯(lián)系在一起。在本文計(jì)算中，由于計(jì)算模型為三維模型（見圖1），因此網(wǎng)格劃分采用直角坐標(biāo)系六面體結(jié)構(gòu)網(wǎng)格，其中計(jì)算區(qū)域在x-y面上布置，z為豎直方向。

利用黑崗口水庫和汴西湖分別建立模型，整體邊界為計(jì)算范圍。模型采用自適應(yīng)邊界較好的正交網(wǎng)格進(jìn)行模擬。網(wǎng)格均勻分布，模型計(jì)算時(shí)上游開邊界由實(shí)測引水流量控制，含沙量采用典型代表值。網(wǎng)格計(jì)算采用自適應(yīng)邊界較好的三維正交網(wǎng)格進(jìn)行模擬。

進(jìn)口邊界條件按照實(shí)際工況進(jìn)行確定。上邊界為水流與空氣交界面，采用VOF模型。床面及壁面則采用標(biāo)準(zhǔn)壁面函數(shù)法計(jì)算。下游出口邊界條件由黑崗口水庫實(shí)測水位控制。

模擬初始條件按照水庫正常蓄水位進(jìn)行確定，為了更好地推算未來黑崗口水庫的水沙變化趨勢，計(jì)算時(shí)長設(shè)定為10，計(jì)算步長設(shè)置為24h。

2 模型驗(yàn)證計(jì)算

根據(jù)現(xiàn)狀黑崗口水庫建立三維數(shù)值模型，并根據(jù)實(shí)際邊界條件和初始條件來驗(yàn)證該模型的可靠性。初始狀態(tài)按照黑崗口水庫設(shè)計(jì)說明書進(jìn)行設(shè)置，進(jìn)出口條件按照開封黃河河務(wù)局提供的涵閘水文資料中黑崗口閘多年（1999-2007年）平均引水流量（Q=2.76m³/s）進(jìn)行設(shè)定，實(shí)際情況黑崗口水庫流速分布見圖2，2016年黑崗口水庫航拍圖見圖3。

將按照實(shí)際情況模擬出的黑崗口水庫流速分布圖與航拍圖進(jìn)行初步對比，可以看出引水口至湖心島之間的水流流速較大且沖刷嚴(yán)重。由于航拍圖中無法看出實(shí)測數(shù)據(jù)，因此選取代表點(diǎn)（A、B、C）處的水深和泥沙淤積深度模擬結(jié)果與實(shí)測數(shù)據(jù)進(jìn)行對比，見表1。

根據(jù)對比可以確定該模型較為精確且計(jì)算結(jié)果處于可接受誤差范圍之內(nèi)，因此可以利用該模型進(jìn)行下一步的研究分析。

3 模擬計(jì)算及分析

3.1 模擬計(jì)算工況

汴西湖建成后每年經(jīng)西干渠輸水5387萬m³，將對黑崗口水庫的水沙條件造成一定影響。根據(jù)黑崗口水庫引水閘引水流量及汴西湖引水流量設(shè)定不同的工況進(jìn)行仿真計(jì)算，進(jìn)口含沙量采用多年平均值。水沙計(jì)算工況見表2。

水沙仿真計(jì)算結(jié)果見圖4。由圖4可知：在修建汴西湖后，由于開啟了分水口，因此水流向西南方向偏移，整體含沙量分布也向汴西湖引水方向轉(zhuǎn)移;原有分布中，湖中島東側(cè)、沉沙廊道北側(cè)沉沙聚集區(qū)域沉沙量減小，泥沙遷移至西干渠進(jìn)水口與沉沙廊道之間，以及湖中島與對應(yīng)的南岸之間;黑崗口水庫南岸沿線形成淤積;西干渠輸水渠道采用倒虹吸形式進(jìn)行輸水，輸水口底板高程較高，汴西湖輸水后，進(jìn)口端具有一定的沉沙量，端口處及兩岸（寬2～5m）產(chǎn)生一定的淤積（淤積深度0.5～1.0m）。

3.2 引水方案分析對比

汴西湖修建后，由于大部分水量需要由黑崗口水庫供給，同時(shí)汴西湖規(guī)劃中有分水口引水，因此需要對汴西湖從黑崗口水庫引水的方案進(jìn)行對比。目前可供選取的方案主要有以下幾個(gè)。

方案一：每月平均從黑崗口水庫引水，引水流量保持在2m³/s左右（工況3、工況4）。該方案下不考慮黑崗口水庫從黃河的引水量及開封市的降水量，模擬發(fā)現(xiàn)雖然引水量保持不變，但沒有考慮汴西湖的蒸發(fā)量及未來分水口的引水量，因此整個(gè)汴西湖及黑崗口水庫的整體水量在全年中并不穩(wěn)定。豐水期由于降水量較大，因此黑崗口水庫水位上升，而引水流量相對較小，會造成黑崗口水庫部分區(qū)域泥沙淤積嚴(yán)重，可以從圖4（c）、圖4（d）中看出，黑崗口引水閘后至湖心島之間的泥沙淤積情況較為嚴(yán)重，若干年后在淤積作用下，湖心島范圍會逐漸擴(kuò)大，從而影響黑崗口水庫中的水流，也會影響東側(cè)供水渠的引水。

方案二：不考慮時(shí)間等要素，只根據(jù)汴西湖水量變化引水，引水流量在0～5m³/s之間調(diào)節(jié)。該方案沒有考慮豐水及枯水期影響，僅僅根據(jù)汴西湖水量調(diào)節(jié)，可能會造成部分極端工況，如在豐水期由于汴西湖水量可能較大，因此引水量減小甚至不引水（工況2），此時(shí)黑崗口水庫水量會有較大變化，水位明顯上升，同時(shí)水庫中水流流速較小，泥沙淤積會出現(xiàn)極端情況。根據(jù)模擬推算可以得出，該種極端情況下黑崗口水庫整體泥沙淤積會明顯增大。另外，若枯水期汴西湖水量較小，僅僅考慮補(bǔ)水而將連接汴西湖和黑崗口水庫的西干渠閘門全開大流量（5m³/s）引水（工況5），則會造成大量泥沙向水庫西南處引水口推移，泥沙淤積主要集中在引水渠口周圍，從而使西干渠喪失輸水功能。

方案三：考慮豐水期和枯水期因素調(diào)節(jié)引水流量。開封市年內(nèi)降水集中在5-8月，因此黑崗口水庫在該時(shí)間段內(nèi)水位升高，此時(shí)應(yīng)體現(xiàn)出汴西湖的分流作用，將連接黑崗口水庫與汴西湖的西干渠閘門全開，保證汴西湖最大引水流量（工況6）。該種工況下會對黑崗口水庫的泥沙進(jìn)行重新調(diào)節(jié)，局部淤積可能有所減弱。開封市枯水期為11月—翌年1月，為保證黑崗口水庫的水量，此時(shí)西干渠引水量應(yīng)根據(jù)情況減?。ㄒ髁?～2m³/s）（工況1、工況3），維持汴西湖水位在最低水位之上。在該工況下，由于水庫水量及引水流速減小，因此泥沙淤積分布較為廣泛，應(yīng)做好黑崗口水庫的挖沙清淤工作。

不同引水方案下黑崗口水庫水位及泥沙分布情況見表3。

4 結(jié)論

汴西湖作為黑崗口水庫的調(diào)蓄水工程，經(jīng)過建立數(shù)值模型的分析計(jì)算與推演，同時(shí)綜合不同引水方案的對比，得出優(yōu)化方案即根據(jù)實(shí)際降水量的時(shí)間分布來調(diào)節(jié)輸水量，豐水期引水、枯水期調(diào)蓄將對黑崗口水庫的泥沙淤積及分布起到優(yōu)化作用，在滿足景觀需要的同時(shí)也可以起到防洪等作用。在黑崗口水庫整體的泥沙分布方面，未修建汴西湖前黑崗口水庫引水過程中由于黃河水中含沙量較高，因此在黑崗口引水渠道出口至水庫中心島之間產(chǎn)生大量泥沙淤積，通過模擬計(jì)算及與實(shí)測數(shù)據(jù)進(jìn)行對比分析，得到淤積厚度達(dá)到3m。修建汴西湖后，由于西南處開設(shè)一流量較大的引水口，在往汴西湖輸水的過程中，泥沙分布整體向西南方向偏移500～700m，分水口附近泥沙厚度達(dá)到1.5～2.0m，因此應(yīng)按照實(shí)際情況及模擬結(jié)果對清淤工作進(jìn)行一定的調(diào)整。同時(shí)，在西干渠引水閘位置處由于引水?dāng)嗝媸湛s，因此造成泥沙淤積厚度達(dá)0.5m，距離水面僅有0.2～0.6m，目前實(shí)際工程中設(shè)置有排沙網(wǎng)，可以起到一定的阻礙泥沙的作用，但通過模擬推算后發(fā)現(xiàn)，該位置若不采取相應(yīng)的清淤排沙措施，則會造成引水效率下降，若干年后連接黑崗口水庫與汴西湖西干渠的倒虹吸將會因泥沙淤積問題而喪失正常輸水功能。

參考文獻(xiàn)：

[1]張紅藝，楊明，張俊華，等.高含沙水庫泥沙運(yùn)動數(shù)學(xué)模型的研究及應(yīng)用[J].水利學(xué)報(bào)，2001，32（11）：20-25.

[2]楊程，李春光，景何仿，等.寧夏沙坡頭河段二維水沙運(yùn)動數(shù)值模擬研究[J].水力發(fā)電，2013，39（5）：8-11，23.

[3]馮淑萍，楊程.沙坡頭庫區(qū)三維水沙運(yùn)移數(shù)值模擬研究[J].水力發(fā)電，2014，40（7）：22-25，56.

[4]假冬冬，邵學(xué)軍，周建銀，等.水沙條件變化對河型河勢影響的三維數(shù)值模擬研究[J].水力發(fā)電學(xué)報(bào)，2014，33（5）：108-113.

[5]王福軍.計(jì)算流體動力學(xué)分析[M].北京：清華大學(xué)出版社，2004：58-76.

[6]耿敬，馬世領(lǐng)，李剛，等.龍頭橋水庫溢洪道流場數(shù)值模擬及優(yōu)化研究[J].水力發(fā)電學(xué)報(bào)，2017，36（10）：74-83.

[7]關(guān)大瑋，程香菊.基于VOF模型的溢流壩三維水流數(shù)值模擬[J].中國水運(yùn)（下半月），2015，15（6）：51-55，59.

[8]CHUNG R S，KIM J，GE W，et al.Reducing Erosion ofEarthen Levees Using Engineered Flood Wall Surface[J].Journal of Geotechnical and Geoenvironmental Engineering，2011，137（10）：874-881.

[9]楊明遠(yuǎn)，嚴(yán)以新，孔俊.珠江口水流泥沙運(yùn)動模擬研究[M].北京：海洋出版社，2008：50-56.

[10]湯立群，陳界仁，陳國祥，等.白沙水庫-維泥沙數(shù)學(xué)模型及電廠取水方案計(jì)算[J].河海大學(xué)學(xué)報(bào)，1998（5）：37-41.

[11]MAJID A Rhinder，AURELIEN Babarit，LIONEL Gentaz，et al.Potential Time Domain Model with Viscous Correctionand CFD Analysis of a Generic Surging Floating Wave En-ergy Converter[J].International Journal of Marine Energy，2015，10：70-96.

[12]JI Un，JANG Fun kyung，KIM Gwonhan.Numerical Mod-eling of Sedimentation Control Scenarios in the ApprocachHannel of the Nakdong River Estuary Barrage，South Korea[J].International Journal of Sediment Research，2016，31（3）：257-263.