基于數(shù)值模擬的南湖流場(chǎng)優(yōu)化與水質(zhì)改善研究

蔡 一，尚福強(qiáng)，唐莉華，孔 彬，王 碩，王可平，馬進(jìn)龍

(1.清華大學(xué)水利水電工程系，北京 100084； 2.中電建路橋集團(tuán)有限公司，北京 100048)

很多城市景觀湖泊由于結(jié)構(gòu)相對(duì)封閉、缺乏有效水量交換，導(dǎo)致水體流速較慢，自凈化能力弱；同時(shí)，為滿足城市景觀需求，湖岸多設(shè)計(jì)成曲折多彎型，容易形成流動(dòng)性更差的“死水”區(qū)域。由于湖體缺乏有效的水量置換手段，導(dǎo)致各類污染物、營養(yǎng)物質(zhì)等易于沉積，進(jìn)而導(dǎo)致水環(huán)境惡化[1-2]。

針對(duì)城市景觀湖泊水體流動(dòng)性差的特點(diǎn)，人工水循環(huán)是一種行之有效的湖泊水環(huán)境治理手段[3]。有關(guān)研究表明，湖泊型景觀水體主要由季節(jié)性降雨與出入流控制湖區(qū)流量大小，且流量決定著其他水力因子的分布與變化特征[4]。通過人為增設(shè)取水口延長引水時(shí)間，加大引水流量[4-5]，可以有效地提高湖泊水體交換能力，促進(jìn)水體循環(huán)流動(dòng)，改善城市景觀湖泊水動(dòng)力條件，實(shí)現(xiàn)流速的提高與流場(chǎng)分布的優(yōu)化，加速污染物的遷移與轉(zhuǎn)化[6]。因此，可以通過人工方式對(duì)流場(chǎng)進(jìn)行改變，從而達(dá)到對(duì)水體環(huán)境質(zhì)量的改善效果[7-9]。

對(duì)湖泊流場(chǎng)分布的研究方法主要包括野外觀測(cè)法、物理試驗(yàn)法和數(shù)值模擬法等3類[7]，其中，數(shù)值模擬法通過數(shù)值計(jì)算方法求解數(shù)學(xué)模型，能較為精確地再現(xiàn)湖泊水體流動(dòng)，是應(yīng)用最為廣泛的研究方法，如EFDC(environmental fluid dynamics code)、WASP(water quality analysis simulation program)和MIKE系列軟件等都是較為成熟的數(shù)值模擬模型，在國內(nèi)應(yīng)用較多[10-11]。由于一般湖泊寬淺的形態(tài)特點(diǎn)，目前主要采用平面二維模型進(jìn)行湖泊流場(chǎng)模擬，如美國陸軍工程兵團(tuán)和揚(yáng)·伯明翰大學(xué)研發(fā)的SMS-RMA2模型、荷蘭代爾夫特理工大學(xué)研發(fā)的Delft2D、丹麥水利研究所研發(fā)的MIKE21等軟件。此類機(jī)理性模型的理論基礎(chǔ)堅(jiān)實(shí)、平臺(tái)工具成熟，能夠?qū)Ω鞣N工況下水域的水動(dòng)力-水質(zhì)情況進(jìn)行全方位的模擬，在實(shí)際工作中更具適用性。其中，MIKE21軟件經(jīng)過多年來的完善發(fā)展，已經(jīng)被廣泛應(yīng)用于包括湖泊、水庫在內(nèi)的各類水體的模擬研究中[9,11-13]，為指導(dǎo)相關(guān)湖泊的水環(huán)境治理提供了科學(xué)工具。

本文以典型城市景觀湖泊山東鄆城南湖為例，基于MIKE21軟件進(jìn)行湖區(qū)流場(chǎng)和濃度場(chǎng)的模擬，提出多個(gè)水源熱泵回水口設(shè)置方案，以提高湖區(qū)水體的流動(dòng)性，改善湖區(qū)水質(zhì)，并通過分析不同方案對(duì)流場(chǎng)的改善效果，尋求最佳設(shè)置方案。

1 研究區(qū)概況

鄆城縣南湖公園位于山東省菏澤市西南部平原區(qū)，東臨德商路、西臨勝利街、南臨東溪路、北臨裕民路，水面面積約27.2 hm2，平均水深約為3.6 m。南湖承接了周邊約2.2 km2范圍的雨水徑流，通過圖1所示的5個(gè)雨水管進(jìn)入湖中，是湖區(qū)主要的污染來源。同時(shí)，南湖北岸的南湖賓館空調(diào)采用水源熱泵系統(tǒng)，通過湖區(qū)北岸的取水口將湖水引入熱泵系統(tǒng)用于空調(diào)熱量交換，然后通過回水管退水到湖中，設(shè)計(jì)引水量800 m3/h，退水升溫不超過1℃。由于湖中水量遠(yuǎn)大于引退水量，本文暫不考慮水溫變化對(duì)流場(chǎng)和濃度場(chǎng)的影響。取水口由并列的2個(gè)1.5 m×2 m矩形口組成，回水口由6個(gè)直徑為0.478 m的圓管束組成。取水口及回水口位置如圖1所示。南湖為人工城市景觀湖泊，因此湖盆平坦且湖岸坡度較大。所在的鄆城縣地處華北平原，屬于暖溫帶半濕潤東亞季風(fēng)大陸性氣候，四季分明，多年年均降水量為695 mm，年均蒸發(fā)量為1 860 mm。

南湖是典型的城市封閉景觀水體，水體流動(dòng)性差，自凈化能力弱，易受污染且不容易自我修復(fù)。根據(jù)2018年6月進(jìn)行的水質(zhì)采樣和檢驗(yàn)數(shù)據(jù)，按照GB 3838—2002《國家地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》，湖區(qū)水質(zhì)主要指標(biāo)如溶解氧、氨氮、化學(xué)需氧量的質(zhì)量濃度都在Ⅳ～Ⅴ類水范圍內(nèi)，葉綠素a含量處于富營養(yǎng)化水平，未達(dá)到一般景觀水域要求，受到當(dāng)?shù)毓芾聿块T的關(guān)注。從空間分布來看，湖區(qū)中央及湖區(qū)南部的水樣水質(zhì)較差。

圖1 雨水管及水源熱泵取水和回水口位置Fig.1 Locations of rainwater pipes and inlet and outlet of water-source heat pump

2 模型構(gòu)建和參數(shù)確定

2.1 模型構(gòu)建

基于MIKE21軟件建立南湖湖區(qū)二維平面水動(dòng)力-水質(zhì)耦合模型，對(duì)湖區(qū)流場(chǎng)及濃度場(chǎng)分布進(jìn)行模擬計(jì)算。建模主要采用以下模塊：

a.水動(dòng)力模塊。水動(dòng)力模塊是MIKE21軟件的基礎(chǔ)，能夠模擬研究區(qū)域的水位和流速變化，主要應(yīng)用于平面二維自由表面流，模塊能夠靈活地識(shí)別干濕單元。水動(dòng)力模塊主要包含了連續(xù)性方程和動(dòng)量方程[14-16]。

b.水質(zhì)模塊。水質(zhì)模塊包括對(duì)流擴(kuò)散模塊和水質(zhì)反應(yīng)模塊Ecolab。對(duì)流擴(kuò)散模塊可以模擬物質(zhì)在水體中的對(duì)流和擴(kuò)散過程；Ecolab模塊通過設(shè)置變化量的反應(yīng)方程實(shí)現(xiàn)對(duì)于復(fù)雜水質(zhì)情況的模擬[17-18]，同時(shí)也可以考慮各物質(zhì)之間的相互轉(zhuǎn)換關(guān)系，主要用于全氮素、葉綠素、大腸桿菌濃度等復(fù)雜反應(yīng)體系的模擬。

利用網(wǎng)格編輯器對(duì)南湖進(jìn)行三角形網(wǎng)格劃分，基本單元尺寸為30 m，并對(duì)彎曲岸邊局部進(jìn)行加密處理，如圖2所示。計(jì)算區(qū)域網(wǎng)格總數(shù)為8 496，最小單元面積為8.0 m2。

利用MIKE21軟件中的水動(dòng)力模塊與水質(zhì)反應(yīng)模塊Ecolab，依次添加降水、蒸發(fā)、風(fēng)場(chǎng)等氣象數(shù)據(jù)并設(shè)置相關(guān)參數(shù)。氣象數(shù)據(jù)來源于國家氣象信息中心《中國地面氣候資料日值數(shù)據(jù)集(V3.0)》，本文采用2018年逐日實(shí)測(cè)氣象數(shù)據(jù)序列，計(jì)算分析湖區(qū)的流場(chǎng)分布。模型的邊界輸入條件為5個(gè)雨水管的水量、水質(zhì)和水源熱泵的抽退水量，初始條件為2018年1月的水位和水質(zhì)數(shù)據(jù)。

圖2 南湖水下地形和計(jì)算網(wǎng)格Fig.2 Underwater topography and model mesh of Nanhu Lake

2.2 模型參數(shù)確定及驗(yàn)證

水動(dòng)力模塊需要考慮的基本參數(shù)包括起始時(shí)間、時(shí)間步數(shù)、主時(shí)間步長、CFL(Courant-Friedrich Levy)收斂條件數(shù)、底床摩擦力、科氏力、風(fēng)場(chǎng)、冰蓋、引潮勢(shì)、降水量、蒸發(fā)量、波浪輻射應(yīng)力等?？紤]到南湖面積不大、水深較小，水體流動(dòng)較為緩慢以及南湖周邊單位用水情況等特點(diǎn)，模型設(shè)置時(shí)忽略科氏力、冰蓋、引潮力、波浪、熱效應(yīng)等的影響。為保證模型運(yùn)行穩(wěn)定，采用CFL收斂條件數(shù)作為判斷計(jì)算的收斂條件，參考文獻(xiàn)[15]，設(shè)置為0.8，代入式(1)，可計(jì)算得到相應(yīng)的主時(shí)間步長，即最大允許時(shí)間步長為12 s。本文采用10 s作為計(jì)算步長。

(1)

式中：C為CFL收斂條件數(shù)；g為重力加速度，m/s2；h為水深，m；u、v分別為x、y方向的流速，m/s；Δt為時(shí)間步長，s；Δx、Δy分別為x、y方向的網(wǎng)格長度，m。

參考文獻(xiàn)[19-21]選取Smagorinsky系數(shù)表征渦流粘度，其值為0.28，選取曼寧系數(shù)表示底床摩擦，其值為32 m1/3/s。風(fēng)場(chǎng)對(duì)上層水流影響較大[22]，風(fēng)摩擦力設(shè)定為風(fēng)速的函數(shù)，水面風(fēng)應(yīng)力的大小取決于水面風(fēng)速的強(qiáng)弱，可由經(jīng)驗(yàn)公式得出：

(2)

水質(zhì)計(jì)算模塊參數(shù)眾多，特別是Ecolab模塊中的很多參數(shù)無法通過試驗(yàn)測(cè)得。因此，綜合參考相關(guān)文獻(xiàn)設(shè)置初始參數(shù)，并通過控制點(diǎn)模擬，與實(shí)測(cè)結(jié)果對(duì)比，對(duì)參數(shù)進(jìn)行調(diào)整，最終確定合適的參數(shù)。

南湖為新建人工湖泊，缺乏足夠的實(shí)測(cè)數(shù)據(jù)，因此，僅通過水量平衡方法對(duì)流場(chǎng)模擬結(jié)果進(jìn)行合理性分析?；?018年模擬結(jié)果，湖區(qū)直接降水量為17.9萬m3，管道入湖的雨水徑流量為70.3萬m3，湖面蒸發(fā)量為41.1萬m3，湖底滲漏量為50.1萬m3。根據(jù)水量平衡，2018年南湖水量減少了3.1萬m3。結(jié)合南湖地形高程情況，湖水位平均下降約0.12 m。根據(jù)現(xiàn)場(chǎng)有關(guān)資料，南湖2018年的水位實(shí)際下降約0.15 m，與模型所得結(jié)果比較接近。采樣點(diǎn)水質(zhì)指標(biāo)模擬值與實(shí)測(cè)值對(duì)比見表1，相對(duì)誤差均在15%以內(nèi)，因此，可以采用該模型對(duì)不同方案的流場(chǎng)和濃度場(chǎng)進(jìn)行模擬與分析。

表1 南湖水質(zhì)指標(biāo)的實(shí)測(cè)值與模擬值對(duì)比Table 1 Comparison of observed and simulated values of water quality indicators in Nanhu Lake

3 計(jì)算方案和結(jié)果分析

3.1 計(jì)算方案

南湖岸邊的南湖賓館空調(diào)系統(tǒng)采用水源熱泵，每天抽取湖水約1.9萬m3，并通過回水口退回湖中，形成湖水的循環(huán)利用。取水口及原設(shè)計(jì)的回水口位置見圖1。為了提高湖水的流動(dòng)性，擬在不同位置增設(shè)回水口，利用MIKE21軟件進(jìn)行不同方案的流場(chǎng)模擬和分析。

根據(jù)對(duì)現(xiàn)狀流場(chǎng)的模擬結(jié)果，分別增設(shè)2個(gè)、3個(gè)和4個(gè)回水口，具體方案和相應(yīng)回水口的位置及各方案回水口流量見表2和圖3，各方案的退水總量不變。

表2 不同方案及設(shè)計(jì)流量Table 2 Design schemes of outlets and discharges

3.2 結(jié)果分析

根據(jù)模擬結(jié)果，采用整個(gè)流場(chǎng)分析與代表點(diǎn)分析相結(jié)合的方法，對(duì)不同方案下的流場(chǎng)進(jìn)行比較。由于降水會(huì)擾動(dòng)湖面，導(dǎo)致湖區(qū)流場(chǎng)變化，為避免降水對(duì)流場(chǎng)的擾動(dòng)，選擇前1周內(nèi)未發(fā)生明顯降水的2018年3月16日作為典型日來考察湖區(qū)的流場(chǎng)分布。原方案下的流場(chǎng)分布情況如圖4所示，不同方案改善后的流場(chǎng)分布情況如圖5所示。

(a) 方案1

圖4 原方案下2018年3月16日0:00的流場(chǎng)分布Fig.4 Fluid distribution at 0:00 of March 16, 2018 in original scheme

(a) 方案1

由圖4可以看出，原方案的湖中心存在大片低速區(qū)，且流速長期維持在0.006 m/s以下，流場(chǎng)狀況不佳。由圖5可以看出，增設(shè)回水口之后，其附近一定時(shí)間內(nèi)均發(fā)生了流場(chǎng)的變化，特別是原來接近“靜水”的湖心區(qū)域流速明顯增加，增幅超過80%；在各增加排水口方案下，低流速區(qū)域的面積均有一定程度的減小，流場(chǎng)狀況變好；對(duì)比方案1和方案4、方案2與方案3都可以看出，隨著回水口的數(shù)量增加，湖泊中的低流速區(qū)域明顯變小。參考已有研究結(jié)果[23]認(rèn)為，流速增加有利于水體的混合，并能提高水面的復(fù)氧能力，從而增強(qiáng)湖泊的自凈能力，即流速越大越有利于湖泊水環(huán)境的改善。本研究中由于回水口的出水打破了原有的“死水”環(huán)境，使得湖水能夠?qū)崿F(xiàn)更快、更大范圍內(nèi)的流動(dòng)，有利于湖水的充分混合。

流場(chǎng)的改變?cè)诒苊馕廴疚锞植烤奂耐瑫r(shí)更有利于提高湖泊的自凈能力。為了驗(yàn)證流場(chǎng)改變對(duì)水環(huán)境的改善效果，選取流速變化較明顯的典型位置(圖4中的1號(hào)和2號(hào)點(diǎn))，分析不同方案下溶解氧和氨氮月均質(zhì)量濃度變化的模擬結(jié)果(圖6和圖7)。

(a) 溶解氧

由模擬結(jié)果可以看出，兩個(gè)典型位置的水質(zhì)變化規(guī)律相似，增設(shè)回水口數(shù)量能夠在一定程度上改善南湖的局部水質(zhì)。增加回水口方案下溶解氧質(zhì)量濃度有所升高，變化幅度在5%以內(nèi)；與溶解氧相比，增加回水口方案對(duì)氨氮的削減效果更為顯著。

(a) 溶解氧

其中，冬春季由于入湖污染負(fù)荷少，水質(zhì)改善效果尤為明顯，在各方案下，1號(hào)和2號(hào)點(diǎn)氨氮質(zhì)量濃度春季降幅普遍超過30%，顯著改善了附近水域的水質(zhì)。由此可見，不同回水口方案對(duì)流場(chǎng)的改變，不僅表現(xiàn)在流速變化上，還表現(xiàn)在對(duì)湖水的混合作用。流速增大可以提高污染物的降解速率，同時(shí)由混合作用帶來的稀釋效果，也會(huì)提升水體自凈能力。因此氨氮質(zhì)量濃度的變化幅度比流速和溶解氧質(zhì)量濃度的增幅更為顯著。

對(duì)比不同方案的結(jié)果表明，各方案之間的差異明顯，特別是在春季的作用效果差異較大。其中，方案3對(duì)氨氮的削減效果最好，1號(hào)點(diǎn)在方案3下的年均氨氮質(zhì)量濃度較原方案降低了61.2%，方案1、方案2和方案4下的年均氨氮質(zhì)量濃度較原方案分別降低了19.0%、34.6%和10.1%；2號(hào)點(diǎn)在方案3下的年均氨氮質(zhì)量濃度較原方案降低了58.2%，而方案1、方案2和方案4下的年均氨氮質(zhì)量濃度較原方案分別降低了9.0%、38.8%和11.5%?？梢钥闯觯桨?對(duì)南湖水質(zhì)改善最為有利。

4 結(jié) 論

a.基于MIKE21軟件的二維平面水動(dòng)力-水質(zhì)耦合模型，模擬得到不同方案下的流場(chǎng)和濃度場(chǎng)，結(jié)果表明，回水口對(duì)流場(chǎng)的影響主要受其布設(shè)位置與數(shù)量的影響，增設(shè)回水口改善了湖區(qū)的流場(chǎng)分布，改變了封閉湖體存在大面積低流速區(qū)域的現(xiàn)狀，對(duì)改善湖區(qū)水質(zhì)有較好的效果。

b.根據(jù)模型模擬結(jié)果，各方案對(duì)湖區(qū)水質(zhì)改善的效果不同，以氨氮為分析指標(biāo)，方案3(即設(shè)置4個(gè)回水口方案)的水質(zhì)改善效果最好，與原方案相比，兩個(gè)觀測(cè)點(diǎn)氨氮質(zhì)量濃度的平均降幅為 59.7%，有利于南湖的水質(zhì)改善。

c.模型參數(shù)參考了相關(guān)研究成果進(jìn)行設(shè)定，并根據(jù)水量平衡分析對(duì)其合理性進(jìn)行了評(píng)估，今后可結(jié)合現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)進(jìn)一步對(duì)模型進(jìn)行完善，以提高其模擬精度。