基于數(shù)值方法的大型水庫水溫結(jié)構(gòu)影響因素分析

鄒 琳，劉 偉，蔣定國

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

0 引 言

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

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

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

1 計算方法

1.1 Delft 3D-Flow模塊計算原理

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

1.2 計算網(wǎng)格

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

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

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

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

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

2 計算方案與邊界條件

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表5 β值與F值對比表Tab.5 Table of β and F values

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

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

圖6 不同入流水溫周期下庫尾水溫歷時對比圖Fig.6 Comparison diagram of reservoir tail water temperature duration under different periodic inflow water temperature

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

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

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

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

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

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

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

4 結(jié) 論

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