磨盤山水庫(kù)溫度場(chǎng)的數(shù)值模擬

伍悅濱，徐 瑩，田 禹，張海龍

(1.哈爾濱工業(yè)大學(xué)城市水資源與水環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室，哈爾濱150090，ybwu@hit.edu.cn; 2.哈爾濱工業(yè)大學(xué)土木工程學(xué)院，哈爾濱150090;3.哈爾濱供水工程有限責(zé)任公司，哈爾濱150076)

水庫(kù)蓄水帶來了防洪、發(fā)電、灌溉，供水等綜合效益，同時(shí)也引起了河流水文、泥沙、生態(tài)等各方面的環(huán)境影響，隨著社會(huì)經(jīng)濟(jì)的發(fā)展，大量工業(yè)廢水和城市污水未經(jīng)處理排入河流、水庫(kù)，超過了水體的自凈能力，河流，水庫(kù)污染日益嚴(yán)重，影響著人類的生存和發(fā)展.水溫是水質(zhì)的一個(gè)重要影響因素，因此進(jìn)行水庫(kù)水溫分析和預(yù)測(cè)，是水庫(kù)環(huán)境評(píng)價(jià)、規(guī)劃、治理和管理的基礎(chǔ).

磨盤山水庫(kù)位于拉林河干流上游五常市沙河子鄉(xiāng)沈家營(yíng)村的上游1.8 km附近，拉林河是松花江右岸的一級(jí)支流，流域總面積19200 km2.除供水外，磨盤山水庫(kù)兼有防洪、灌溉及其他綜合利用效益.水庫(kù)的建成將改變庫(kù)區(qū)及大壩下游河段的水文和水環(huán)境狀況，進(jìn)而對(duì)水庫(kù)下游河段的灌溉和生活用水產(chǎn)生重大影響.

1 數(shù)學(xué)模型

EFDC(Environmental Fluid Dynamics Code)是由美國(guó)國(guó)家環(huán)保署資助開發(fā)，用于模擬河流、湖泊、水庫(kù)、海灣、濕地和河口等水環(huán)境系統(tǒng)中流場(chǎng)、溫度場(chǎng)、濃度場(chǎng)及生態(tài)過程的三維數(shù)值計(jì)算模型［1－2］.對(duì)水庫(kù)水溫進(jìn)行模擬的數(shù)學(xué)模型主要采用Mellor等［3－5］提出的控制方程.

(1)連續(xù)方程

(2)曲線正交坐標(biāo)中x向和y向動(dòng)量方程

(3)溫度

式中:u和v分別是沿著正交曲線水平方向x向和y向的流速分量;mx和my是平面坐標(biāo)轉(zhuǎn)換因子，m=mxmy.

其中，z為無量綱垂向坐標(biāo)，w*為物理垂向速度;水深H=h+ζ，h代表實(shí)際垂向的底高程;ζ代表水位;動(dòng)量方程中f表示科里奧利(Coriolis)參數(shù); Av為垂向紊動(dòng)粘性系數(shù);Qu和Qv表示動(dòng)量方程中的源匯項(xiàng);θ為水溫;在溫度方程中的Qθ項(xiàng)是外源輸入輸出項(xiàng).連續(xù)方程中的QSS和QSW表示河床與水體之間的泥沙及水的交換通量，其正值表示通量由河床入水體中，反之為負(fù)值.三角函數(shù)δ(0)表示這類通量進(jìn)入水體的最底層.

2 數(shù)值計(jì)算及分析

磨盤山水庫(kù)為典型的河道型水庫(kù)，干流回水長(zhǎng)約13 km.地形總趨勢(shì)為東南高、西北低，呈不對(duì)稱的“U”字型，地面坡降約為3.8%.磨盤山水庫(kù)壩高45 m，正常蓄水位318 m，相應(yīng)的水庫(kù)總庫(kù)容355.92×106m3，調(diào)節(jié)庫(kù)容為292.95×106m3.汛期限制水位為 317.00 m，重復(fù)利用庫(kù)容為27.27×106m3.設(shè)計(jì)洪水位為318.77 m，相應(yīng)的設(shè)計(jì)總庫(kù)容為375.98×106m3.校核洪水位為320.89 m，相應(yīng)的校核總庫(kù)容為443.02×106m3.

2.1 網(wǎng)格劃分

庫(kù)區(qū)地形資料是進(jìn)行水庫(kù)模擬計(jì)算的基礎(chǔ)，首先根據(jù)水庫(kù)地形高程圖，按照模型要求劃分計(jì)算網(wǎng)格.計(jì)算中時(shí)間步長(zhǎng)為10 s，運(yùn)行周期為1 d，即86 400 s，模型運(yùn)行時(shí)間設(shè)定為30個(gè)周期.控制方程的數(shù)值解法為有限體積法和有限差分法的結(jié)合，對(duì)于模擬對(duì)象，沿流動(dòng)方向采用正交曲線網(wǎng)格離散，而在垂直方向上則采用分層方法進(jìn)行離散，進(jìn)而引入了層厚的概念.根據(jù)磨盤山水庫(kù)的地理特點(diǎn)，將其水平方向上劃分1 052個(gè)矩形網(wǎng)格，每個(gè)單元格大小為85 m×65 m.考慮溫度分層效應(yīng)，垂直方向上分為8層，占水深比例分別為12.5%，12.5%，20%，20%，20%，5%，5%，5%.水庫(kù)概化及計(jì)算區(qū)域單元格見圖1.

圖1 水庫(kù)概化及計(jì)算區(qū)域單元格

2.2 初始條件

拉林河流域內(nèi)現(xiàn)有沈家營(yíng)、五常等13處水文站，根據(jù)磨盤山水庫(kù)工程地點(diǎn)和資料情況，選用的主要分析測(cè)站是拉林河的沈家營(yíng)測(cè)站.測(cè)站觀測(cè)項(xiàng)目有降水、蒸發(fā)、氣溫、氣壓、日照、溫度和風(fēng)速等以及磨盤山庫(kù)區(qū)的日氣溫變化.初始水位水深根據(jù)地形圖得到，初始水溫根據(jù)水庫(kù)管理區(qū)測(cè)量資料以及《水利水電工程水文計(jì)算規(guī)范》(SDJ214－83)(試行)條例中給出.地面坡降為3.8%，初始計(jì)算流速設(shè)為零.

2.3 邊界條件

根據(jù)庫(kù)區(qū)管理站測(cè)量的2007年入庫(kù)流量和出庫(kù)流量作為上、下游邊界條件，入庫(kù)流量過程線見圖2，出庫(kù)流量過程線見圖3.

圖2 入流流量過程線

圖3 出流流量過程線

2.4 模型驗(yàn)證

通過對(duì)磨盤山水庫(kù)枯水期、平水期和豐水期的模擬，得到了水庫(kù)出流溫度變化曲線，并與實(shí)際出流溫度變化進(jìn)行了比較，結(jié)果如圖4～6所示.模擬結(jié)果與實(shí)測(cè)值基本吻合，說明計(jì)算模型及方法的可信性.

圖4 2007年4月壩前水溫實(shí)測(cè)與模擬值對(duì)比

圖5 2007年6月壩前水溫實(shí)測(cè)與模擬值對(duì)比

圖6 2007年9月壩前水溫實(shí)測(cè)與模擬值對(duì)比

2.5 模擬結(jié)果及分析

圖7 壩前斷面的溫度分層

本文分別模擬了2007年枯水期(四月)、平水期(六月)、豐水期(九月)磨盤山水庫(kù)的水溫變化規(guī)律.壩前水溫在枯水期、平水期、豐水期的溫度場(chǎng)計(jì)算結(jié)果分別見圖7.水溫在4月份為升溫期，入流水溫和表面熱通量逐漸增加，上層水體水溫增長(zhǎng)迅速;在4月末，表層水溫已達(dá)到10℃，但是靠近底部的水溫幾乎沒有變化，仍保持低于6℃的低溫，在水面附近出現(xiàn)了較大的溫度梯度，形成了較明顯的水溫分層.6月份為高溫期，入流水溫和氣溫都達(dá)到最高.在表面很薄的水體中，由于水氣熱交換形成了表層溫躍層［6－7］，表層水溫保持在15℃作用，底部的水溫緩慢上升了2～3℃，但仍然保持低于8℃的溫度，在整個(gè)垂向斷面上的溫差達(dá)到了接近10℃的水平，在這一時(shí)期，水體上部分保持穩(wěn)定的溫度分層現(xiàn)象.9月份為降溫期，氣溫和入流水溫逐漸下降，水體向大氣散失熱量，冷水下沉使表層溫躍層消失，水庫(kù)在整體垂直方向形成水溫分層，并開始出現(xiàn)層與層之間的垂向摻混，即秋季翻滾現(xiàn)象［8－10］.

水庫(kù)的這種分層結(jié)構(gòu)對(duì)水質(zhì)具有重要的影響.對(duì)于磨盤山水庫(kù)而言，由于其溫度分層作用，不同深度的水質(zhì)狀況將有所不同.在枯水期，水庫(kù)水溫混合基本均勻，導(dǎo)致水中的營(yíng)養(yǎng)物質(zhì)沿垂向均勻混合.又加之冬季水中生化反應(yīng)緩慢，消耗較少的營(yíng)養(yǎng)物質(zhì)，至春季時(shí)沉積了大量營(yíng)養(yǎng)物質(zhì)，為藻類快速生長(zhǎng)提供了物質(zhì)基礎(chǔ).平水期處于夏季，水庫(kù)上部分水體水溫分層較為明顯.此時(shí)溫躍層在垂直方向上具有“屏障”作用［11－13］，水中的營(yíng)養(yǎng)物質(zhì)在垂直方向上的交換量被大大的削減.在豐水期，汛期來臨，洪水挾帶著沖起的河底淤沙，涌進(jìn)水庫(kù)，使大量的亞磷酸鹽進(jìn)入水流.水庫(kù)水溫的穩(wěn)定分層被破壞后，“屏障”消失，水庫(kù)底層的沉淀物和有害物質(zhì)會(huì)翻滾到表層，對(duì)水體造成污染.

圖8是枯水期、平水期和豐水期在水庫(kù)輸水洞斷面的模擬水溫分布剖面圖.磨盤山水庫(kù)的輸水洞位于右岸，進(jìn)口位于右岸壩軸線上游約250 m處.水體自輸水洞進(jìn)入長(zhǎng)直管線，向凈水廠輸送原水.由模擬結(jié)果可知，水庫(kù)在枯水期和平水期沿深度方向上水溫分布穩(wěn)定，水溫較低，可以直接選取水面以下10 m的水庫(kù)水作為飲用水水源;在豐水期，水溫結(jié)構(gòu)發(fā)生變化.隨著氣溫降低，表層水溫也下降，相應(yīng)的水密度增大，導(dǎo)致表層水開始下沉，底層水上升，在整個(gè)深度范圍內(nèi)表層和底層發(fā)生了混合，這種秋季翻滾現(xiàn)象可能會(huì)引發(fā)水質(zhì)問題［14－15］.

圖8 輸水洞斷面水溫分布

3 結(jié)論

1)通過實(shí)測(cè)值與模擬值的對(duì)比驗(yàn)證，應(yīng)用EFDC建立的水庫(kù)水動(dòng)力學(xué)模型，能夠很好的模擬水庫(kù)中水溫的時(shí)空分布.因此基于該模型進(jìn)行水體中流場(chǎng)、溫度場(chǎng)、濃度場(chǎng)及生態(tài)過程的數(shù)值模擬是進(jìn)行水庫(kù)環(huán)境保護(hù)、管理及運(yùn)行的有力工具.

2)磨盤山水庫(kù)屬于分層型水庫(kù)，水溫出現(xiàn)季節(jié)性分層現(xiàn)象，庫(kù)表分層現(xiàn)象最為明顯，至庫(kù)底分層減弱，庫(kù)底水溫變化平緩.枯水期垂向分層現(xiàn)象不明顯，平水期形成表層溫躍層，對(duì)于垂向傳質(zhì)具有一定程度的屏障作用，豐水期出現(xiàn)秋季翻滾，對(duì)水質(zhì)影響較大.

3)水庫(kù)取水口深度應(yīng)該隨季節(jié)調(diào)整.在枯水期和平水期沿深度方向上水溫分布穩(wěn)定，大部分水體水溫較低，水深10 m以下均可布置取水口;而在豐水期出現(xiàn)的秋季翻滾現(xiàn)象，形成水體垂向的自循環(huán)流動(dòng)，可能引發(fā)水質(zhì)問題.因此，豐水期至冬季時(shí)段，取水口位置應(yīng)有所調(diào)整，并建議凈水廠強(qiáng)化水處理工藝.

［1］HAMRICKJ M.A Three-dimensional environmental fluid dynamics computer code:Theoretical and computational aspects［M］.Gloucester，Massachusetts:Virgina Institute of Marine Science，the College of William and Mary，1992:1－63.

［2］李云生，劉偉江，吳悅穎，等.美國(guó)水質(zhì)模型研究進(jìn)展綜述［J］.水利水電技術(shù)，2006(2):68－73.

［3］MELLOR G L，BLUMBERG A F.Modeling vertical and horizontal diffusivities with the sigma coordinate system［J］.Mon Wea Rev，1985，20(5):851－875.

［4］MELLOR G L，YAMADA T.Development of turbulence closure model for geophysical fluid problems［J］.Rev Geophysical Space Phys，1982，20(6):851－875.

［5］ARAKAWA A，LAMBV R.Computational design of the basic dynamical processes of the UCLA general circulation model［M］//CHANG J.Methods in Computational Physics.New York:Academic Press，1977:173－265.

［6］梁瑞峰，李嘉，李克峰，等.立面二維水庫(kù)水溫模型的并行研究［J］.四川大學(xué)學(xué)報(bào)，2009，41(1):30－33.

［7］張士杰，彭文啟.二維水庫(kù)水溫結(jié)構(gòu)及影響因素研究［J］.水利學(xué)報(bào)，2009，40(10):1254－1259.

［8］OUGOLNITSKY G A，USOV A B.Water quality control in watercourses［J］.Water Resources，2003，30(2): 226－232.

［9］申滿斌，陳永燦，劉昭偉，等.岸邊排放污染物濃度場(chǎng)三維渾水水質(zhì)模型研究［J］.水力發(fā)電學(xué)報(bào)，2005，4(3):93－98.

［10］郭磊，高學(xué)平，張晨，等.北大港水庫(kù)水質(zhì)模擬及分析［J］.長(zhǎng)江流域資源與環(huán)境，2007，16(1):11－16.

［11］WOOL T A，DAVIE S R，HUGO N R.Development of three-dimensional hydrodynamic and water quality models to support TMDL decision process for the Neuse River Estuary，North Carolina［J］.Journal of Water Resources Planning and Management，2003，129(4):295－306.

［12］JIN K R，HAMRICK J H，TISDALE T.Application of a three-dimensional hydrodynamic model for lake okeechobee［J］. Journal of Hydranlic Engineering ASCE，2001，126:13－14.

［13］陳景秋，趙萬星，季振剛.重慶兩江匯流水動(dòng)力模型［J］.水動(dòng)力學(xué)研究與進(jìn)展，2005，12:829－835.

［14］JOHNSON B H，KIM K W，HEATH R E，et al.Validation of three-dimensional hydrodynamic model of chesapeake bay［J］.J Hyd Engrg，2003，119:2－20.

［15］李錦繡，杜斌，孫以三.水動(dòng)力條件對(duì)富營(yíng)養(yǎng)化影響規(guī)律探討［J］.水利水電技術(shù)，2005(5):23－28