某水利工程入庫(kù)口三維水動(dòng)力特征模擬分析

白文龍

（鞍山龍逸水利水電建筑工程有限公司，遼寧 鞍山 114100）

0 引言

大伙房水庫(kù)位于撫順市東部，是遼沈地區(qū)重要的生活用水水源地，水質(zhì)質(zhì)量對(duì)下游人民群眾飲水安全十分重要。水庫(kù)上游入庫(kù)河流主要途徑山區(qū)農(nóng)田，防止農(nóng)業(yè)面源污染是防止大伙房水庫(kù)富營(yíng)養(yǎng)化，保證水質(zhì)的重要途徑，主要包括防止農(nóng)藥污染、化肥污染和畜禽養(yǎng)殖業(yè)污染。對(duì)大伙房水庫(kù)入庫(kù)河流之一的社河建立進(jìn)水流量影響的三維水動(dòng)力數(shù)學(xué)模型，對(duì)區(qū)域水環(huán)境具有一定的意義。

1 研究區(qū)概況

大伙房水庫(kù)位于遼寧省東北部，撫順市東部。該壩址位于遼寧省撫順市東州，距撫順市18 km，距沈陽(yáng)市68 km。大伙房水庫(kù)是中華人民共和國(guó)成立后由水利部組織的多功能水利工程。蓄水灌溉，下游工業(yè)用水和生活用水，養(yǎng)殖用地和景觀用地。這也是我國(guó)第一個(gè)自行設(shè)計(jì)建設(shè)的大型水利工程。它的樞紐項(xiàng)目包括溢洪道，水壩，輸水等，并且在水庫(kù)區(qū)域還建有取水建筑物。

農(nóng)業(yè)面源污染是大伙房水庫(kù)富營(yíng)養(yǎng)化的重要因素，主要包括農(nóng)藥污染，化肥污染和畜禽養(yǎng)殖業(yè)污染。社河流經(jīng)撫順縣部分鄉(xiāng)鎮(zhèn)，流入大伙房水庫(kù)，造成水質(zhì)污染。

圖1 水庫(kù)概況圖

2 研究方法

水庫(kù)流量具有某些自然河流和湖泊的特征，并且還具有獨(dú)特的流量特征。為了充分理解污染物在水流中的運(yùn)動(dòng)規(guī)律，需要對(duì)水流的三維水動(dòng)力特性和污染物的運(yùn)動(dòng)進(jìn)行仿真分析。提出了三維水動(dòng)力和水質(zhì)控制方程。垂直使用s坐標(biāo)變換，并使用有限體積法離散和求解控制方程。

2.1 三維流體動(dòng)力學(xué)數(shù)學(xué)模型

在垂直方向上使用s坐標(biāo)變換來(lái)實(shí)現(xiàn)相對(duì)分層。σ坐標(biāo)轉(zhuǎn)換方法將自由水表面和不規(guī)則底 部轉(zhuǎn)化為σ坐標(biāo)系中的表面和底部坐標(biāo)平面。“水深”為1。此方法不僅使整個(gè)計(jì)算的水域在垂直方向上具有相同的網(wǎng)格數(shù)，而且可以任意分層以確保在淺水中具有較高的垂直分辨率。另一方面，就數(shù)值方法而言，σ坐標(biāo)系中方程的離散解要容易得多。因此，垂直s坐標(biāo)變換被廣泛應(yīng)用于水流和輸運(yùn)的三維數(shù)值模擬中。圖2是σ坐標(biāo)和z坐標(biāo)之間的關(guān)系的示意圖[1]。

圖2 σ坐標(biāo)變換示意圖

2.2 網(wǎng)格剖分

有限體積法中使用的非結(jié)構(gòu)化網(wǎng)格通常由三角形或四邊形網(wǎng)格組成。目前，優(yōu)選使用四邊形網(wǎng)格，因?yàn)楫?dāng)網(wǎng)格節(jié)點(diǎn)的數(shù)量相同時(shí)，三角形網(wǎng)格單元的數(shù)量和網(wǎng)格接口的數(shù)量通常是四邊形網(wǎng)格的兩倍或更多，因此計(jì)算量較大，三角網(wǎng)格計(jì)算精度和穩(wěn)定性差；相反，四邊形網(wǎng)格具有良好的計(jì)算穩(wěn)定性，計(jì)算結(jié)果的準(zhǔn)確性，并且二階粘度更易于處理。三角形網(wǎng)格和四邊形網(wǎng)格可以結(jié)合使用，以四邊形為主，三角形為補(bǔ)充。三角網(wǎng)格主要用于局部地形變化，粗細(xì)網(wǎng)格過(guò)渡和鋸齒形邊界[2]。

計(jì)算區(qū)域在平面上劃分為三角形元素，而柱狀元素用于垂直劃分（如圖3所示），計(jì)算單位是空間中的三角棱鏡。垂直分為N層，每層的高度為DSl，l=1，N，所有變量都定義在三棱柱的中心。

圖3 網(wǎng)格劃分示意圖

庫(kù)區(qū)長(zhǎng)35 km，最寬處4 km，最窄處0.3 km，最大水深34 m，平均水深12 m。渾河、蘇子河和社河是三大水庫(kù)河。庫(kù)區(qū)底部的地形如圖4所示。計(jì)算區(qū)被一個(gè)非結(jié)構(gòu)化的三角網(wǎng)格劃分，并在社河庫(kù)區(qū)進(jìn)行了適當(dāng)?shù)募用堋Ｆ矫嫔峡偣膊贾昧? 619個(gè)節(jié)點(diǎn)和4 711個(gè)網(wǎng)格元素。網(wǎng)格側(cè)約250 m長(zhǎng)，加密部分約100 m長(zhǎng)。垂直方向分為等距離的各層，并分為10層。網(wǎng)格分布如圖4所示。

圖4 庫(kù)區(qū)計(jì)算區(qū)域網(wǎng)格布置

2.3 定解條件及參數(shù)取值

（1）初始條件

將初始流速設(shè)為0，初始水位為131 m。

（2）邊界條件

岸邊界采用不可訪問(wèn)的邊界條件，水流的法向速度為0，V×n =0，其中V為速度矢量，n為岸邊界的法向矢量。

社河流量變化很大，因此社河進(jìn)入水庫(kù)的邊界是流量的時(shí)間序列。出水邊界采用水位邊界條件，正常高水位129.85 m作為出水邊界條件。

（3）參數(shù)值

用Smagorinsky模型計(jì)算水平渦流粘度系數(shù)，cs值為0.28；用Kolmogorov-Prandl求解垂直渦粘度系數(shù)，其中k和e通過(guò)求解ke方程獲得。在沒(méi)有相關(guān)粗糙度數(shù)據(jù)的情況下，粗糙度通常為0.002～0.01m，根據(jù)相關(guān)研究，其值為0.01 m。

3 研究結(jié)果

3.1 總磷變化趨勢(shì)

8月4日至8月8日入流和總磷濃度隨時(shí)間的變化曲線見(jiàn)圖5。自8月4日降雨以來(lái)，入流量隨時(shí)間先增加，達(dá)到479 m3/s后，流速下降。儲(chǔ)存中總磷濃度的總體趨勢(shì)也先上升然后下降，但是在此過(guò)程中會(huì)反復(fù)出現(xiàn)曲折。

圖5 社河入庫(kù)流量、總磷隨時(shí)間變化曲線

3.2 庫(kù)區(qū)三維水動(dòng)力特征時(shí)間演變分析

本節(jié)主要將入水量作為水動(dòng)力模擬的上限條件，將相應(yīng)的壩址出口水位作為模擬的下限條件，并以平均風(fēng)速2.3 m/s和主導(dǎo)風(fēng)向（東北風(fēng)）為基礎(chǔ)。濕季為自由水面條件，建立一個(gè)考慮到來(lái)水流影響的水庫(kù)區(qū)域三維水動(dòng)力數(shù)學(xué)模型，分析庫(kù)區(qū)三維水動(dòng)力特征[3]。

當(dāng)模型計(jì)算達(dá)到穩(wěn)定狀態(tài)時(shí)，將研究區(qū)域中自由水表層（表層）的流動(dòng)模式用于分析。模擬了12 h，20 h，36 h，50 h，65 h，80 h 時(shí)在儲(chǔ)層區(qū)域中的地表流場(chǎng)分布。

根據(jù)不同時(shí)期庫(kù)區(qū)地表流場(chǎng)分布：當(dāng)模擬時(shí)間為12 h、20 h和36 h時(shí)，社河水庫(kù)入庫(kù)流量分別為478 m3/s，248 m3/s，153 m3/s。由于流入儲(chǔ)層的射流很大，因此湍流對(duì)儲(chǔ)層流場(chǎng)的影響要大于風(fēng)力流對(duì)儲(chǔ)層流場(chǎng)的影響。社河進(jìn)水口的水流結(jié)構(gòu)呈現(xiàn)出由霍夫和蓬勃流控制的水動(dòng)力特性以及顯著的下游流量特性。當(dāng)模擬時(shí)間為50 h、65 h、80 h時(shí)，社河的入水量分別為86 m3/s，48 m3/s和31 m3/s。流入水庫(kù)的射流逐漸減小，而風(fēng)速保持不變。社河水庫(kù)入口附近的水流速度降低，水庫(kù)中的水流速度增加。大伙房水庫(kù)房進(jìn)水口的水流結(jié)構(gòu)顯示出以風(fēng)為主導(dǎo)的水動(dòng)力特性。受盛行的風(fēng)向東北風(fēng)的影響，水庫(kù)中的一部分水將流向社河進(jìn)口，形成順時(shí)針的大循環(huán)，水從社河流入水庫(kù)。

當(dāng)模擬時(shí)間為 12 h、20 h、36 h 時(shí)，吞吐流對(duì)水庫(kù)流場(chǎng)的影響大于風(fēng)生流對(duì)水庫(kù)流場(chǎng)的影響，大伙房水庫(kù)社河入庫(kù)口水域水流結(jié)構(gòu)呈現(xiàn)受吞吐流支配的水動(dòng)力特征和顯著的順岸流特征。當(dāng)模擬時(shí)間為 50 h、65 h、80 h 時(shí)，風(fēng)生流對(duì)水庫(kù)流場(chǎng)的影響大于吞吐流對(duì)水庫(kù)流場(chǎng)的影響，大伙房水庫(kù)社河入庫(kù)口水域水流結(jié)構(gòu)呈現(xiàn)以風(fēng)生流為支配的水動(dòng)力特征。

3.3 庫(kù)區(qū)水動(dòng)力三維空間分布特征

當(dāng)模型計(jì)算達(dá)到穩(wěn)定狀態(tài)時(shí)，分析研究區(qū)域內(nèi)自由水表面層（表層），層4（中間層）和床層表層（底部層）的流動(dòng)模式。

比較不同時(shí)間的表層，中間層和底層的流程圖：每一時(shí)刻底層流場(chǎng)的速度大于中間層流場(chǎng)的速度[4，5]。這是因?yàn)轱L(fēng)應(yīng)力方向和水流入水庫(kù)壩址的流向是不同的，并且地表流場(chǎng)受風(fēng)向的影響很大。

在模擬時(shí)間12 h和20 h，各層的水流結(jié)構(gòu)基本相同，且方向均沿社河入庫(kù)水庫(kù)至壩址方向，水流運(yùn)動(dòng)表現(xiàn)出良好的下游特性流。但是，在仿真時(shí)間36 h、50 h、65 h和80 h中，不同層的流動(dòng)結(jié)構(gòu)是不同的。此時(shí)，地表流場(chǎng)顯示了西北海岸的陸上流動(dòng)特征和東北海岸的近海流動(dòng)特征。中流場(chǎng)顯示了西北岸和東北岸的沿海水流特征。西北海岸的底部流場(chǎng)顯示了近海流的特征，并顯示了東北海岸的陸上流的特征。此外，在模擬時(shí)間為50 h、65 h和80 h時(shí)，底部和中間流場(chǎng)的順時(shí)針循環(huán)比表層的順時(shí)針循環(huán)更明顯。

每一時(shí)刻底層流場(chǎng)的速度都大于中間層流場(chǎng)的速度和表面層流場(chǎng)的速度。在12 h和20 h的模擬時(shí)間內(nèi)，不同層的流動(dòng)結(jié)構(gòu)基本相同，其方向?yàn)檠厣绾铀畮?kù)方向的壩址方向，模擬時(shí)間為36 h、50 h、65 h和80 h具有不同的水流結(jié)構(gòu)層。此外，在模擬時(shí)間為50 h、65 h和80 h時(shí)，底部和中間流場(chǎng)的順時(shí)針循環(huán)比表層的順時(shí)針循環(huán)更明顯。

4 結(jié)論

以大伙房水庫(kù)社河入庫(kù)口附近水域?yàn)檠芯繀^(qū)域，建立了以σ為坐標(biāo)系的三維水動(dòng)力數(shù)學(xué)模型。采用非結(jié)構(gòu)三角網(wǎng)格劃分計(jì)算區(qū)域，采用有限體積法離散求解控制方程，建立了考慮來(lái)水流量影響的庫(kù)區(qū)三維水動(dòng)力數(shù)學(xué)模型。得出的結(jié)論是，大伙房水庫(kù)社河入水口在12 h、20 h和36 h呈現(xiàn)出受吞流和顯著下游流量特性支配的水動(dòng)力特性。在50 h、65 h和80 h時(shí)，大伙房水庫(kù)大棚進(jìn)水口的水流結(jié)構(gòu)呈現(xiàn)出受風(fēng)流支配的水動(dòng)力特性。36 h、50 h、65 h和80 h不同層的流動(dòng)結(jié)構(gòu)不同。