王宇 秦凈凈 職保平

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

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

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

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

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

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

1 模擬計(jì)算模型

1.1 控制方程

數(shù)值模型計(jì)算控制方程[5]包括連續(xù)性方程、動(dòng)量方程、紊動(dòng)能k方程、紊動(dòng)能耗散率ε方程、挾沙力計(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ù)。

動(dòng)量方程為式中：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為作用在流體微元上的壓力。

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

紊動(dòng)能耗散率ε方程為式中：σ_ε為湍流耗散率對(duì)應(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ì)挾沙力。

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

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

1.2 模型建立及邊界條件

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

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

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

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

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

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

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

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

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

3.1 模擬計(jì)算工況

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

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

3.2 引水方案分析對(duì)比

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

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

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

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

不同引水方案下黑崗口水庫(kù)水位及泥沙分布情況見(jiàn)表3。

4 結(jié)論

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

參考文獻(xiàn)：

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

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

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

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

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

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

[7]關(guān)大瑋，程香菊.基于VOF模型的溢流壩三維水流數(shù)值模擬[J].中國(guó)水運(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)動(dòng)模擬研究[M].北京：海洋出版社，2008：50-56.

[10]湯立群，陳界仁，陳國(guó)祥，等.白沙水庫(kù)-維泥沙數(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.