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

白文龍

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

0 引言

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

1 研究區(qū)概況

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

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

圖1 水庫概況圖

2 研究方法

水庫流量具有某些自然河流和湖泊的特征，并且還具有獨(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)變換來實(shí)現(xiàn)相對(duì)分層。σ坐標(biāo)轉(zhuǎn)換方法將自由水表面和不規(guī)則底 部轉(zhuǎn)化為σ坐標(biāo)系中的表面和底部坐標(biāo)平面?！八睢睘?。此方法不僅使整個(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)格過渡和鋸齒形邊界[2]。

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

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

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

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

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

（1）初始條件

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

（2）邊界條件

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

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

（3）參數(shù)值

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

3 研究結(jié)果

3.1 總磷變化趨勢

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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