鄒 琳，劉 偉，蔣定國(guó)

（1.三峽大學(xué)水利與環(huán)境學(xué)院，湖北宜昌 443002；2.中國(guó)長(zhǎng)江三峽集團(tuán)有限公司，北京 100038）

0 引 言

水電作為我國(guó)的重要能源資源，對(duì)國(guó)家經(jīng)濟(jì)建設(shè)和社會(huì)發(fā)展具有重大作用。然而水電在實(shí)現(xiàn)發(fā)電效益、促進(jìn)節(jié)能減排的同時(shí)，對(duì)河流生態(tài)環(huán)境也產(chǎn)生了環(huán)境效應(yīng)［1，2］，如夏季下泄水溫明顯低于下游河道水溫，從而對(duì)農(nóng)作物產(chǎn)量［3］、魚類的生長(zhǎng)繁殖產(chǎn)生不利影響［4-6］。隨著水庫(kù)的建成蓄水，庫(kù)區(qū)內(nèi)水體流速減緩，與未建水庫(kù)前的天然河道相比，水溫過(guò)程發(fā)生明顯改變，而作為水環(huán)境評(píng)價(jià)的重要因子之一，水溫問(wèn)題目前在水利工程生態(tài)效應(yīng)和水庫(kù)運(yùn)行管理中是一個(gè)重要課題，也是備受關(guān)注的熱點(diǎn)問(wèn)題［7-9］。

目前，國(guó)內(nèi)外針對(duì)水庫(kù)水溫研究方法大都采用現(xiàn)場(chǎng)監(jiān)測(cè)和數(shù)值模擬［10-13］。學(xué)者們基于國(guó)內(nèi)外多個(gè)水庫(kù)的實(shí)測(cè)資料，提出了許多經(jīng)驗(yàn)性水溫估算方法如東勘院張大發(fā)法［14］、朱伯芳法［15］以及垂向水溫分布公式［16］等。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，數(shù)值模擬成為研究水庫(kù)水溫結(jié)構(gòu)及演變規(guī)律的重要工具。數(shù)學(xué)模型的研究從一維［17］開(kāi)始，但垂向一維數(shù)學(xué)模型適用范圍有一定限制，一般用于回水長(zhǎng)度短且深的水庫(kù)，在入流影響較大的水庫(kù)，垂向一維數(shù)學(xué)模型并不適用。為提高水庫(kù)水溫模擬精度，立面二維數(shù)學(xué)模型隨之產(chǎn)生，在國(guó)內(nèi)外得到了廣泛應(yīng)用［18，19］。由于水流運(yùn)動(dòng)具有三維特性，所以水溫問(wèn)題實(shí)際上是一個(gè)三維問(wèn)題，因此，近年來(lái)國(guó)內(nèi)外研究學(xué)者將研究方向轉(zhuǎn)向三維數(shù)學(xué)模型［20-22］。任華堂［23］等建立三維水溫?cái)?shù)值模型對(duì)阿海水庫(kù)進(jìn)行模擬，研究發(fā)現(xiàn)取水口位置對(duì)下泄水溫、均溫層位置產(chǎn)生影響。張士杰［24］等利用MIKE3 建立二灘水庫(kù)的三維水溫模型，研究了上游來(lái)水水溫對(duì)水庫(kù)水溫結(jié)構(gòu)的影響，結(jié)果表明，隨著上游來(lái)水溫度升高，水庫(kù)水溫相應(yīng)升高。黃廷林［25］等建立黑河金盆水庫(kù)的三維水溫模型，研究發(fā)現(xiàn)氣溫與金盆水庫(kù)水體熱分層的穩(wěn)定性有顯著的正相相關(guān)關(guān)系。

河道型水庫(kù)的水動(dòng)力條件、熱力學(xué)性質(zhì)復(fù)雜，自然條件下，水溫結(jié)構(gòu)的影響因素眾多、共同作用、相互耦合，僅根據(jù)所有影響因素共同作用后的溫度值進(jìn)行分析，難以探明各個(gè)因素的作用程度和范圍。故以西南地區(qū)某水庫(kù)實(shí)測(cè)資料概化成的模型水庫(kù)作為研究對(duì)象，采用數(shù)值模擬方法對(duì)重要影響因素的單獨(dú)作用情況進(jìn)行模擬，分析其對(duì)水庫(kù)水溫結(jié)構(gòu)的影響程度和范圍，為相關(guān)研究和工程設(shè)計(jì)提供資料參考。

1 計(jì)算方法

1.1 Delft 3D-Flow模塊計(jì)算原理

水庫(kù)主要在垂向上出現(xiàn)水溫分層，同時(shí)受發(fā)電、灌溉以及泄洪的影響，水庫(kù)水流具有明顯的三維特征，因此本文采用Delft 3D-Flow 模塊進(jìn)行三維模擬。Delft 3D 是由荷蘭代爾夫特研究所開(kāi)發(fā)，能模擬二維和三維的水體流動(dòng)、波浪運(yùn)動(dòng)、水質(zhì)演變等，其中水動(dòng)力（Flow）模塊是模型的核心。為了減少由于邊界不規(guī)則引入的離散化誤差，水平方向采用正交曲線坐標(biāo)系，垂向上采用σ坐標(biāo)系。模型在淺水和Boussinesq 假定下求解不可壓縮流體的納維-斯托克斯方程，紊流模型使用k-ε模型，熱輸運(yùn)采用三維的對(duì)流-擴(kuò)散方程，方程中利用源項(xiàng)和匯項(xiàng)來(lái)模擬熱量的進(jìn)出，在水氣界面考慮了風(fēng)、輻射、熱對(duì)流等因素作用，使得模型可以模擬復(fù)雜的熱輸運(yùn)過(guò)程［26，27］。模型得到了董菁水庫(kù)、香溪河流域、大寧河流域等實(shí)測(cè)資料驗(yàn)證［28，29］。

1.2 計(jì)算網(wǎng)格

根據(jù)西南地區(qū)某水庫(kù)實(shí)測(cè)資料概化成模型水庫(kù)，水庫(kù)長(zhǎng)40 km，庫(kù)容2.000 8億m3，模型從上游入口至800 m處為平底，入口河底高程98 m，從800 m 處開(kāi)始變坡，底部高程逐漸降低，至壩前處為斜坡段，壩前河底高程0 m，河道寬100 m。水庫(kù)劃分網(wǎng)格時(shí)采用正交曲線網(wǎng)格，網(wǎng)格的余弦值小于0.02，滿足計(jì)算要求，生成的網(wǎng)格數(shù)為75 000 個(gè)，沿河道方向網(wǎng)格數(shù)為1 000 個(gè)，橫向網(wǎng)格數(shù)為5個(gè)，垂向分15層。

圖1 模型水庫(kù)網(wǎng)格圖Fig.1 Model reservoir grid diagram

1.3 模型參數(shù)、率定及驗(yàn)證

采用西南某水庫(kù)2012年實(shí)測(cè)水溫?cái)?shù)據(jù)進(jìn)行對(duì)比驗(yàn)證，并對(duì)模型參數(shù)進(jìn)行率定，率定后的模擬結(jié)果見(jiàn)圖2，可以看出，計(jì)算得到的壩前水庫(kù)水溫值與實(shí)測(cè)值變化趨勢(shì)基本相同，絕對(duì)誤差最大值0.32 ℃，模擬結(jié)果較好，模擬精度可滿足要求。此時(shí)曼寧系數(shù)為0.037，水平渦流黏滯系數(shù)為1.0 m2/s，斯坦頓數(shù)為0.001 3，重力加速度為9.81 m/s2，水的密度為1 000 kg/m3。

圖2 壩前水溫驗(yàn)證結(jié)果Fig.2 Water temperature verification results in front of dam

2 計(jì)算方案與邊界條件

水庫(kù)水溫的影響因素有水文水力因素、氣象因素以及水庫(kù)幾何特征等，現(xiàn)有研究表明入流條件、氣象條件以及取水口位置對(duì)水溫分布影響顯著［24］，因此，重點(diǎn)分析入流邊界和水氣界面交換對(duì)水庫(kù)水溫結(jié)構(gòu)的影響，同時(shí)考慮取水口位置的影響。模型上游入流邊界采用流量邊界，下游設(shè)置取水口，為簡(jiǎn)化計(jì)算，入流邊界主要考慮入流流量和入流水溫，水氣界面考慮氣溫的周期性變化，其余因素如云量、相對(duì)濕度、輻射量、風(fēng)速等參考西南地區(qū)實(shí)際情況進(jìn)行輸入。

入流流量是水庫(kù)水溫結(jié)構(gòu)的重要影響因素，為分析不同入流流量對(duì)水庫(kù)水溫結(jié)構(gòu)的影響程度和范圍，考慮取水口位置的不同，設(shè)置不同的入流流量進(jìn)行模擬分析，具體方案見(jiàn)表1。水庫(kù)初始水溫設(shè)為15 ℃，初始水位設(shè)為100 m，采用流量入流邊界；出庫(kù)流量設(shè)定值與入庫(kù)流量相同，入流水溫設(shè)為20 ℃，取水口高7 m，寬20 m。

表1 入流流量方案Tab.1 Incoming flow scheme

入流水溫的周期性變化也會(huì)對(duì)水庫(kù)水溫結(jié)構(gòu)產(chǎn)生影響，為分析不同周期的入流水溫對(duì)水庫(kù)水溫結(jié)構(gòu)的影響程度和范圍，設(shè)置不同周期的入流水溫進(jìn)行模擬分析，具體方案見(jiàn)表2。水庫(kù)初始水溫設(shè)為15 ℃，初始水位設(shè)為100 m，采用流量入流邊界，出庫(kù)流量設(shè)定值與入庫(kù)流量相同，為100 m3/s，來(lái)流水溫簡(jiǎn)化為正弦變化，均值15 ℃，振幅5 ℃，即入流水溫在10～20 ℃之間進(jìn)行周期性變化，周期分別為12、24、168 h，取水口高7 m，寬20 m，取水口底板高程75 m。

表2 周期性入流水溫方案Tab.2 Periodic inflow water temperature scheme

不同時(shí)間尺度的水氣界面邊界條件對(duì)水庫(kù)水溫結(jié)構(gòu)的影響效果不同，主要考慮氣溫，其他因素像云量、相對(duì)濕度、輻射量、風(fēng)速等參考西南地區(qū)實(shí)際情況進(jìn)行輸入，為分析不同周期的氣溫對(duì)水庫(kù)水溫結(jié)構(gòu)的影響程度和范圍，設(shè)置不同周期的氣溫進(jìn)行模擬分析，具體方案見(jiàn)表3。水庫(kù)初始水溫設(shè)為15 ℃，初始水位設(shè)為100 m，采用流量入流邊界，出庫(kù)流量設(shè)定值與入庫(kù)流量相同，為100 m3/s，入流水溫為15 ℃，氣溫設(shè)為正弦變化，均值15 ℃，振幅5 ℃，即氣溫在10 ℃至20 ℃之間進(jìn)行周期性變化，周期分別為7、30、90 d，取水口高7 m，寬20 m，取水口底板高程75 m。

表3 周期性氣溫方案Tab.3 Periodic temperature scheme

3 模擬結(jié)果與分析

3.1 入流流量對(duì)水庫(kù)水溫結(jié)構(gòu)影響分析

水庫(kù)初始水溫設(shè)為15 ℃，來(lái)流水溫設(shè)為20 ℃，可以通過(guò)水庫(kù)水溫的變化過(guò)程來(lái)分析入流流量對(duì)水庫(kù)水溫結(jié)構(gòu)的影響。圖3 表示在入流流量為100 m3/s 的情況下，取水口底板高程為85 m、75 m 時(shí)的水庫(kù)水溫歷時(shí)圖。由圖可知，在同一來(lái)流流量下，隨著水深的增加，水庫(kù)水溫達(dá)到與入流水溫同溫所需要的時(shí)間增長(zhǎng)，同時(shí)對(duì)比圖3（a）和（b）還可以看出，以取水口底板高程為分界線，在分界線以上，水庫(kù)水溫升溫速度較快，而在分界線以下，水體的升溫速度明顯減緩。這是受取水口位置影響，取水口底板高程以上水體，流速大，摻混強(qiáng)，因而可在一定時(shí)間內(nèi)使這部分水體水溫趨于均勻，而取水口底板高程以下水體，流速小，摻混弱，升溫極其緩慢。

圖3 入流流量為100 m3/s時(shí)水庫(kù)水溫歷時(shí)圖Fig.3 Diagram of water temperature duration when inflow flow is 100 m3/s

圖4 比較了當(dāng)入流流量分別為50、100、150、200 和300 m3/s時(shí)，在水庫(kù)壩前斷面處取水口底板高程位置，水庫(kù)水體的升溫過(guò)程。由圖4 可知，在水庫(kù)同一位置處，隨著入流流量的增大，水溫升溫所需要的時(shí)間減少。這是因?yàn)殡S著入口斷面處入流流量的增大，水體流速增大，置換作用增強(qiáng)，因而可以加速水體混合時(shí)間。

圖4 不同入流流量下取水口底板高程處水溫升溫過(guò)程對(duì)比圖Fig.4 Comparison diagram of water temperature warming process at intake floor elevation under different inflow flows

為分析入流流量、水庫(kù)水溫變化時(shí)間和水庫(kù)庫(kù)容間的關(guān)系，擬采用公式（1）對(duì)其進(jìn)行分析。

式中：q為入流流量，m3/s；t為入流流量作用時(shí)間，s；V為入流流量影響的庫(kù)容，m3。

由于入流流量對(duì)取水口底板高程以下水體影響不大，僅考慮取水口底板高程以上水體，V=V0-V1，其中V0為水庫(kù)全部庫(kù)容，V1為取水口底板高程以下庫(kù)容。

表4 為根據(jù)公式（1）計(jì)算所得的值，其中t0為水庫(kù)受入流影響水溫開(kāi)始發(fā)生變化的時(shí)間，以水溫變化率大于10-4（℃/s）的點(diǎn)為準(zhǔn)，此時(shí)計(jì)算所得的為F0，t1為水庫(kù)受入流影響水溫變化至穩(wěn)定狀態(tài)所需時(shí)間，以水溫變化率小于10-4（℃/s）的點(diǎn)為準(zhǔn)，此時(shí)計(jì)算所得的為F1，可見(jiàn)入流流量、水庫(kù)水溫變化時(shí)間和水庫(kù)庫(kù)容存在對(duì)應(yīng)關(guān)系。在遇到F小于F0的洪量，對(duì)水庫(kù)水溫結(jié)構(gòu)一般沒(méi)有影響，在遇到F大于F0且小于F1的洪量，水庫(kù)水溫結(jié)構(gòu)發(fā)生變化，在遇到F等于F1的洪量，水庫(kù)水溫結(jié)構(gòu)變化至另一穩(wěn)定狀態(tài)，此后的來(lái)水不會(huì)再對(duì)水溫結(jié)構(gòu)產(chǎn)生影響。根據(jù)F值，在已知水庫(kù)取水口位置和入流流量的情況下，可以大致預(yù)估水庫(kù)水溫結(jié)構(gòu)。

表4 入流流量、水庫(kù)水溫變化時(shí)間和庫(kù)容關(guān)系表Tab.4 Table of relation between inflow discharge，water temperature change time and storage capacity

傳統(tǒng)的判別水庫(kù)水溫結(jié)構(gòu)的β法定義為一次洪水總量比總庫(kù)容，分層型水庫(kù)在遇到β大于1 的洪量，會(huì)變成臨時(shí)混合型，在遇到β小于0.5 洪量，對(duì)水庫(kù)水溫結(jié)構(gòu)一般沒(méi)有影響，在遇到β在0.5 至1 之間的洪量，對(duì)水庫(kù)水溫結(jié)構(gòu)的影響一般介于兩者之間。表5 以入流流量100 m3/s 為例，在水庫(kù)取水口底板高程為85、75 m 的情況下，β值和F值的對(duì)比表。由表5 可知，在取水口底板高程為75 m 時(shí)，當(dāng)t=300 h，洪量為1.08 億m3時(shí)，β值為0.54，F(xiàn)值為1.50，β法和F值法均表明此次洪量已經(jīng)對(duì)壩前水庫(kù)水溫結(jié)構(gòu)產(chǎn)生影響；當(dāng)t=496 h，洪量達(dá)到1.785 6 億m3時(shí)，β值為0.89，表明壩前水溫結(jié)構(gòu)已經(jīng)接近混合，而F值為2.48，表明來(lái)水對(duì)水庫(kù)水溫結(jié)構(gòu)影響已經(jīng)結(jié)束，此后的來(lái)水不會(huì)再對(duì)水溫結(jié)構(gòu)產(chǎn)生影響；當(dāng)t=700 h，洪量達(dá)到2.52 億m3時(shí)，β值為1.26，表明壩前水溫結(jié)構(gòu)已經(jīng)混合。由圖5 可知，F(xiàn)值判斷更為準(zhǔn)確。

圖5 不同洪量下水庫(kù)水溫分布圖Fig.5 Distribution diagram of reservoir water temperature at different flood volumes

表5 β值與F值對(duì)比表Tab.5 Table of β and F values

3.2 周期性入流水溫對(duì)水庫(kù)水溫結(jié)構(gòu)影響分析

圖6 表示入流水溫周期分為12、24、168 h（7 d）時(shí)，庫(kù)尾水溫歷時(shí)對(duì)比圖。取水口底板高程75 m，由圖6可知，當(dāng)入流水溫周期為12 h，選取的庫(kù)尾斷面處總水深約為25 m，表中層受來(lái)流水溫周期變化基本一致，底層振幅略有衰減且升溫過(guò)程發(fā)生比較明顯變形；當(dāng)入流水溫周期為24 h，選取的庫(kù)尾斷面處總水深約為26 m，表中層受來(lái)流水溫周期變化基本一致，底層升溫過(guò)程發(fā)生比較明顯變形；當(dāng)入流水溫周期為7 d，選取的庫(kù)尾斷面處總水深約為30 m，表中底層受入流水溫周期變化基本一致。隨著入流水溫周期的增大，影響的庫(kù)尾水深范圍變大，使得小于等于此水深的水庫(kù)庫(kù)尾范圍內(nèi)整體變化趨勢(shì)與入流水溫變化趨勢(shì)基本一致，但仍可以看出入流水溫對(duì)庫(kù)尾的最大影響水深仍在取水口底板高程分界線附近。

圖6 不同入流水溫周期下庫(kù)尾水溫歷時(shí)對(duì)比圖Fig.6 Comparison diagram of reservoir tail water temperature duration under different periodic inflow water temperature

表6列出了不同入流水溫周期在壩前的水溫波動(dòng)范圍。由表可知，不論入流周期如何變化，隨著水深的增加，水溫逐漸降低，這是因?yàn)楫?dāng)入流水溫以周期性進(jìn)行變化時(shí)，高溫水從表層流入水庫(kù)，低溫水逐漸下沉。水溫波動(dòng)范圍的差值先增大后減小，大約在取水口底板高程分界線處差值達(dá)到最大。當(dāng)處于同一水深，隨著入流周期的增大，水溫波動(dòng)范圍也隨之增大。

表6 不同入流水溫周期下水溫波動(dòng)范圍表Tab.6 Water temperature fluctuation range under different periodic inflow water temperature

3.3 周期性氣溫邊界條件對(duì)水庫(kù)水溫結(jié)構(gòu)影響分析

由于水面是和空氣相接觸的，水氣界面需要進(jìn)行熱交換，氣溫在一定程度上會(huì)對(duì)水庫(kù)水溫產(chǎn)生影響。周期性氣溫邊界條件對(duì)水庫(kù)水溫結(jié)構(gòu)影響分析考慮了氣溫的周期性變化，其余因素如云量、相對(duì)濕度、輻射量、風(fēng)速等參考西南地區(qū)實(shí)際情況進(jìn)行輸入。圖7～9表示氣溫周期為7、30、90 d時(shí)，選取壩前斷面的水庫(kù)水溫隨時(shí)間變化圖。由圖可知，以水溫變化大于0.2 ℃的區(qū)域設(shè)為氣溫影響區(qū)域，當(dāng)氣溫變幅一致時(shí)，隨著氣溫周期的增加，對(duì)水庫(kù)水溫的影響范圍逐漸增大，當(dāng)氣溫周期為7 d時(shí)，經(jīng)過(guò)一個(gè)周期后，影響水深約為20 m，表現(xiàn)為影響區(qū)域以內(nèi)水溫升高，當(dāng)氣溫周期為30 d 時(shí)，經(jīng)過(guò)一個(gè)周期后，影響水深約為25 m，表現(xiàn)為影響區(qū)域以內(nèi)水溫升高，當(dāng)氣溫周期為90 d時(shí)，經(jīng)過(guò)一個(gè)周期后，影響整個(gè)影響水庫(kù)水深，表現(xiàn)為隨著水深的增加，水溫逐漸降低。同時(shí)，隨著氣溫周期的增加，表層水溫的波動(dòng)范圍增加，當(dāng)氣溫周期為7 d 時(shí)，表層水溫波動(dòng)范圍為15.0～15.4 ℃，當(dāng)氣溫周期為30 d 時(shí)，表層水溫波動(dòng)范圍為15.0～16.0 ℃，當(dāng)氣溫周期為90 d 時(shí)，表層水溫波動(dòng)范圍為15.0～17.1 ℃。

圖7 氣溫周期為7 d時(shí)壩前水庫(kù)水溫隨時(shí)間變化圖Fig.7 Diagram of water temperature variation with time in front of dam when temperature periodicity is 7 d

圖8 氣溫周期為30 d時(shí)壩前水庫(kù)水溫隨時(shí)間變化圖Fig.8 Diagram of water temperature variation with time in front of dam when temperature periodicity is 30 d

圖9 氣溫周期為90 d時(shí)壩前水庫(kù)水溫隨時(shí)間變化圖Fig.9 Diagram of water temperature variation with time in front of dam when temperature periodicity is 90 d

4 結(jié) 論

針對(duì)概化的模型水庫(kù)，采用數(shù)值方法分析入流流量、入流水溫和氣溫對(duì)大型水庫(kù)水溫結(jié)構(gòu)的影響程度和范圍。研究表明，受取水口位置影響，取水口底板高程以上水體，流速大，摻混強(qiáng)，而取水口底板高程以下水體，流速小，摻混弱，因此對(duì)于大型深水庫(kù)，入流流量、入流水溫對(duì)水庫(kù)水溫結(jié)構(gòu)的影響區(qū)域與取水口高程有一定相關(guān)性；隨著入流流量的增大，水庫(kù)水溫結(jié)構(gòu)變化至另一穩(wěn)定狀態(tài)所需要的時(shí)間減少；隨著入流水溫周期的增大，影響的水深范圍增大，但最大影響水深仍在以取水口底板高程為分界線附近；隨著氣溫周期的增大，對(duì)水庫(kù)水溫的影響范圍逐漸增大，經(jīng)過(guò)一個(gè)周期后，當(dāng)氣溫周期為7 d 時(shí)，影響水深約為20 m，當(dāng)氣溫周期為30 d時(shí)，影響水深約為25 m，當(dāng)氣溫周期為90 d 時(shí)，可影響整個(gè)影響水庫(kù)水深；同時(shí)改進(jìn)后的F值法，僅考慮取水口底板高程以上水體，相比較于β法，可更為準(zhǔn)確地判斷壩前水溫結(jié)構(gòu)的變化。