梁志宏，鄭江麗*，陳秀洪，劉 夏

(1.水利部珠江河口治理與保護(hù)重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510611；2.珠江水利委員會(huì)珠江水利科學(xué)研究院，廣東 廣州 510611)

當(dāng)前，隨著經(jīng)濟(jì)和科技的飛速發(fā)展，中國工業(yè)蓬勃發(fā)展，工業(yè)工程的應(yīng)用范圍逐漸增大，在實(shí)踐生產(chǎn)過程中也漸漸凸顯其優(yōu)越性[1]。與此同時(shí)，建設(shè)項(xiàng)目施工期、運(yùn)行期的生產(chǎn)、生活等用水對(duì)取水河道流量提出了新的要求。項(xiàng)目的取水可靠性直接涉及到項(xiàng)目投入運(yùn)行后的穩(wěn)定產(chǎn)出效率，在提高項(xiàng)目經(jīng)濟(jì)效益方面發(fā)揮著至關(guān)重要的作用[2]。

自1992年Warren等[3]首次提出MIKE21建模系統(tǒng)起，MIKE21已歷經(jīng)近30年的發(fā)展與推廣，目前已廣泛應(yīng)用于多個(gè)領(lǐng)域。紀(jì)忠華等[4]利用MIKE進(jìn)行洪水演進(jìn)模擬，通過設(shè)置邊界條件、驗(yàn)證并率定模型參數(shù)等，推斷出了洪水演進(jìn)規(guī)律并制定了相應(yīng)的防洪方案和策略，為水利工程的防洪工作提供了重要參考。Li等[5]基于氣象水文數(shù)據(jù)和基本水質(zhì)數(shù)據(jù)，利用MIKE21建立了模擬河流湖泊的水質(zhì)擴(kuò)散模型，反映了排污口污水排放對(duì)所在江段的水質(zhì)影響情況，可作為排污口設(shè)置論證的參考依據(jù)。Soltani等[6]運(yùn)用MIKE21進(jìn)行水沙模擬研究，將原有方案與優(yōu)化方案進(jìn)行對(duì)比，提出了最優(yōu)的泥沙防治方案。此外，MIKE21還被用于內(nèi)澇[7-8]、通航[9]、風(fēng)浪[10]、突發(fā)水污染事件[11]和水系整治[12]等模擬研究。

從上述研究中可以看出，經(jīng)過近30年的發(fā)展和推廣，MIKE21已廣泛應(yīng)用于水利行業(yè)的多個(gè)領(lǐng)域[13-20]，但是在建設(shè)項(xiàng)目工程取水可靠性分析方面還尚有不足。取水可靠性分析一般采用水庫調(diào)節(jié)計(jì)算或水資源供需平衡計(jì)算方法，但對(duì)于河網(wǎng)區(qū)或來水受下游回水影響的區(qū)域，采用水動(dòng)力模型分析，如姬戰(zhàn)生等[21]根據(jù)河網(wǎng)調(diào)蓄原理構(gòu)建了以南沙平原骨干河道為基礎(chǔ)的概化河網(wǎng)，利用Mike11建立了水動(dòng)力學(xué)模型，驗(yàn)證了偏不利條件下取水口的取水可靠性。劉樹鋒等[22]針對(duì)咸潮河段展開研究，建立了潭江一維河網(wǎng)水質(zhì)數(shù)學(xué)模型，并與實(shí)測(cè)資料進(jìn)行對(duì)比分析，驗(yàn)證了取水口的取水可靠性，并提出了特枯年份的取水風(fēng)險(xiǎn)規(guī)避措施。邢凱等[23]針對(duì)以巢湖湖區(qū)為水源的建設(shè)項(xiàng)目，明確其取水可靠性分析的主要思路和方法。

因此，本文以上小河附近規(guī)劃建設(shè)項(xiàng)目——梧州市千萬噸級(jí)冶金新材料產(chǎn)業(yè)基地為例，基于MIKE21，建立并求解水動(dòng)力模型，旨在針對(duì)上小河上、下游河道寬度變幅較大的情況，驗(yàn)證已有取水方案的合理性和可靠性，為項(xiàng)目的建設(shè)和運(yùn)行提供科學(xué)支撐。

1 研究區(qū)域概況與資料收集

1.1 研究區(qū)域概況

梧州市千萬噸冶金新材料(鋼鐵基地)臨港經(jīng)濟(jì)區(qū)片區(qū)位于臨港經(jīng)濟(jì)區(qū)南片區(qū)，在梧州市總體規(guī)劃藍(lán)圖上，是梧州市城市工業(yè)發(fā)展的集聚區(qū)域之一，其取水口位于上小河，中途有秋米河匯入后，最終匯入潯江。上小河又名須羅江。發(fā)源于岑溪市白板村大竹頂，自南向北流，縱貫新地鎮(zhèn)、林水鎮(zhèn)，至林水鄉(xiāng)河口自然村匯入潯江，是潯江右岸的一級(jí)支流，河口距長(zhǎng)洲水利樞紐2.1 km處，長(zhǎng)洲水庫2007年開始蓄水后，上小河河口段被長(zhǎng)洲水庫淹沒，成為庫區(qū)的一部分。

本文研究區(qū)域主要包括：①秋米河區(qū)域，秋米河水筒口至秋米河與上小河交匯處，長(zhǎng)度約2.38 km；②上小河區(qū)域，上邊界位于擬建工業(yè)取水口上游約550 m，下邊界至上小河與潯江交匯口，長(zhǎng)度約5.65 km；③潯江區(qū)域，主要根據(jù)現(xiàn)狀梧州市飲用水源保護(hù)區(qū)和長(zhǎng)洲水利樞紐地理位置確定，上邊界距離上小河與潯江交匯口約4 km，下邊界為長(zhǎng)洲水利樞紐，距離上小河與潯江交匯口約2 km。研究區(qū)域見圖1。

圖1 研究區(qū)域示意

1.2 研究資料

1.2.1資料收集

由于上小河距長(zhǎng)洲水利樞紐僅2.1 km，且長(zhǎng)洲水利樞紐多年水文資料先后經(jīng)過了水電部中南院、昆明院、貴陽院、東北院、廣西電力院、珠委以及所屬各省(區(qū))的水利、水電設(shè)計(jì)院等單位的多次審查和校核。因此，本文采用長(zhǎng)洲水利樞紐1946—2017年共72年徑流資料進(jìn)行研究。

1.2.2資料整理

a)徑流計(jì)算。采用數(shù)學(xué)期望公式進(jìn)行經(jīng)驗(yàn)頻率計(jì)算，并用矩法公式估算統(tǒng)計(jì)參數(shù)(多年平均流量均值、Cv、Cs)初始值，再用P-Ⅲ型理論頻率曲線適線，確定長(zhǎng)洲水利樞紐壩址多年平均流量為6 124 m3/s，Cv=0.21，Cs/Cv=2.0。

b)典型年選取。根據(jù)長(zhǎng)洲水利樞紐壩址1946—2017年72年年徑流系列，按照豐、平、枯水年選擇年來水總量接近且年內(nèi)水量分配對(duì)項(xiàng)目取水較為不利(年內(nèi)枯水期來水量小)的年份作為典型年。根據(jù)典型年的年來水量及年內(nèi)月來水過程，采用同倍比方法縮放后得到長(zhǎng)洲水利樞紐壩址各保證率來水量年內(nèi)月分配過程。典型年選取結(jié)果及年內(nèi)月分配過程見表1。

表1 典型年月分配過程 單位：m3/s

2 模型建立及方案設(shè)置

新建取水口位于上小河，中途有秋米河匯入后，最終匯入潯江。由于上下游河寬差別很大，傳統(tǒng)的矩形網(wǎng)格應(yīng)用困難，因此，本次計(jì)算模型采用MIKE21模型。該模型網(wǎng)格節(jié)點(diǎn)布設(shè)靈活，便于局部加密，能精確地?cái)M合彎曲河道邊界，具有算法可靠、計(jì)算穩(wěn)定、界面友好、前后處理方便等突出優(yōu)點(diǎn)，已得到廣泛應(yīng)用。結(jié)合本項(xiàng)目河寬、水系連通等實(shí)際情況，采用零維和二維模型相結(jié)合的技術(shù)手段進(jìn)行研究。

2.1 控制方程

連續(xù)性方程表達(dá)式如下[20]：

(1)

動(dòng)量方程表達(dá)式如下：

(2)

(3)

式中h——靜止水深，即水面到某一基準(zhǔn)面的距離；η——水面相對(duì)于未擾動(dòng)基面的高度，即水位；t——時(shí)間；u、v——x、y方向上的流速分量；us、vs——源項(xiàng)水流流速；g——重力加速度；f——柯氏力參數(shù)；ρ——水密度；ρ0——水的參照密度；sxx、sxy、syy——波浪輻射應(yīng)力分量；Pa——大氣壓力；τsx、τsy——風(fēng)應(yīng)力分量；τbx、τby——底部摩擦應(yīng)力分量；Txx、Txy、Tyy——黏性項(xiàng)分量；S——源匯項(xiàng)。

2.2 邊界條件及初始條件

2.2.1邊界條件

在本研究采用的數(shù)值模式中，需給定2種邊界條件，即閉邊界條件和開邊界條件。所謂開邊界條件即水域邊界條件，可以給定水位、流量或調(diào)和常數(shù)。所謂閉邊界條件即水陸交界條件，計(jì)算水域與陸地交界的固邊界上Γ2有：

(4)

模型上邊界流量平水期采用潯江、上小河和秋米河多年平均流量，枯水期從最不利原則考慮，采用潯江、上小河和秋米河多年最小流量；模型下邊界水位分別采用長(zhǎng)洲水利樞紐在平水期和枯水期的潯江控制水位。根據(jù)上小河實(shí)際開發(fā)情況，上游建有2座水電站、5座攔河壩(用于蓄水發(fā)電或灌溉)，平水期下游來流量有可能小于多年平均流量，因此本次新增一個(gè)平水期計(jì)算方案，上小河來流量按生態(tài)流量考慮，秋米河和潯江按平水期時(shí)的邊界條件考慮。水動(dòng)力模型計(jì)算邊界條件見表2。

表2 水動(dòng)力邊界條件

2.2.2初始條件

(5)

式中ζ0(x,y,t0)、u0(x,y,t0)、v0(x,y,t0)——初始時(shí)刻t0的已知值。

2.2.3活動(dòng)邊界處理

本模型采用干濕點(diǎn)判斷法處理潮灘活動(dòng)邊界，如果模型中的區(qū)域是處在干濕邊交替區(qū)，為了避免模型計(jì)算出現(xiàn)不穩(wěn)定性，設(shè)定一個(gè)干水深(drying depth)、淹沒深度(flooding water depth)和濕水深(wetting depth)。當(dāng)某一單元的水深小于濕水深時(shí)，在此單元上的水流計(jì)算會(huì)被相應(yīng)調(diào)整，即不計(jì)算動(dòng)量方程，僅計(jì)算連續(xù)方程。而當(dāng)水深小于干深度的時(shí)候會(huì)被凍結(jié)而不參與計(jì)算。淹沒深度是用來檢測(cè)網(wǎng)格單元是不是已經(jīng)被淹沒。

2.3 網(wǎng)格布置

二維模型采用三角形網(wǎng)格剖分計(jì)算區(qū)域，三角形網(wǎng)格節(jié)點(diǎn)數(shù)為4 534個(gè)，三角形網(wǎng)格數(shù)為7 950個(gè)。上小河取水口附近的網(wǎng)格分辨率為2 m，可見計(jì)算網(wǎng)格是非常精細(xì)的，能高度擬合岸線，準(zhǔn)確地對(duì)河道進(jìn)行概化。模型網(wǎng)格布置見圖2。

圖2 二維水動(dòng)力模型計(jì)算區(qū)域網(wǎng)格布置

2.4 地形資料

計(jì)算范圍內(nèi)河道水下地形采用各區(qū)域最新實(shí)測(cè)地形，潯江整體采用2015年實(shí)測(cè)河道地形資料，上小河和秋米河采用2020年1月最新實(shí)測(cè)資料。長(zhǎng)洲水利樞紐建成后，桂平至長(zhǎng)洲段成為河道型水庫，當(dāng)水庫為枯水期最低消落水位時(shí)，可回水到達(dá)桂平壩下。受水庫回水的影響，河段水面較天然情況有所抬高，因此該河段會(huì)由于水位的壅高而導(dǎo)致流速有所下降，局部灘段存在泥沙淤積的可能性。但由于長(zhǎng)洲水利樞紐為低水頭徑流式電站，汛期壩前水位降低至18.76 m運(yùn)行，基本恢復(fù)到天然狀態(tài)，且由于長(zhǎng)洲樞紐上游干、支流上水庫的攔截，長(zhǎng)洲樞紐實(shí)際推移質(zhì)產(chǎn)沙面積不大，入庫推移質(zhì)泥沙數(shù)量極少，絕大部分來沙(占全年的96%)集中在汛期，進(jìn)入庫區(qū)的推移質(zhì)泥沙很快就可推移到壩前，進(jìn)而排到下游，因此長(zhǎng)洲水庫泥沙淤積主要是懸移質(zhì)泥沙淤積。相關(guān)研究結(jié)果表明，建庫后泥沙淤積是不連續(xù)的，幾乎全部集中在壩前段和常年回水區(qū)的幾個(gè)深潭內(nèi)，回水末端沒有淤積。非汛期來沙僅占全年的4%，水流夾帶的泥沙要少得多，河床也不會(huì)發(fā)生大的淤積。因此認(rèn)為2015—2020年潯江河段地形基本不變。水動(dòng)力模型計(jì)算區(qū)域河道水下地形見圖3。

圖3 二維水動(dòng)力模型計(jì)算區(qū)域河道水下地形

2.5 模型率定和驗(yàn)證

二維模型的率定和驗(yàn)證采用潯江大橋下游河道右岸蒼梧水廠2017年6—9月實(shí)測(cè)水位數(shù)據(jù)，當(dāng)模型底床糙率取值在0.021～0.029范圍時(shí)，模擬效果較好，模型實(shí)測(cè)和計(jì)算值較為吻合，模擬效果見圖4。對(duì)水位的模擬值及實(shí)測(cè)值進(jìn)行分析，發(fā)現(xiàn)水位模擬值和實(shí)測(cè)值的相對(duì)誤差在0.31%～7.13%，水位Nash-Stucliffe系數(shù)均在0.9以上，表明模型設(shè)置參數(shù)能用于本次研究計(jì)算。

圖4 蒼梧水廠2017年6—9月實(shí)測(cè)與模擬水位

3 結(jié)果分析及對(duì)策討論

3.1 結(jié)果分析

梧州市臨港經(jīng)濟(jì)片區(qū)新建取水口取水需求為15萬m3/d，約1.74 m3/s，取水口處河道底高程約17.8 m，不同方案取水口斷面河道水位情況見表3。

表3 預(yù)測(cè)方案及取水口斷面水位

由模型計(jì)算結(jié)果可知本項(xiàng)目取水后取水口所在斷面水位在方案2下最低，說明平水期，潯江運(yùn)行水位較低時(shí)，若上小河來流量不足，本項(xiàng)目取水風(fēng)險(xiǎn)最大。因此針對(duì)方案2下本項(xiàng)目能否取水進(jìn)行詳細(xì)分析。

方案2下，上小河來流量?jī)H0.57 m3/s，本項(xiàng)目取水必須靠潯江來水補(bǔ)充。根據(jù)MIKE模型計(jì)算的水深分布(圖5)，本項(xiàng)目取水口至上小河潯江河口，有多個(gè)河段水深小于0.5 m，需要重點(diǎn)關(guān)注。

圖5 最不利方案模型計(jì)算水深分布

選取潯江向上小河補(bǔ)水可能阻水點(diǎn)(+t1、+t2、+t3、+t4、+t5)的水深變化(圖6)可知，+t2、+t3、+t4斷面的水深為0 m，因此本項(xiàng)目在方案2下靠潯江來水補(bǔ)充存在取水風(fēng)險(xiǎn)。

圖6 最不利方案下可能阻水點(diǎn)水深變化

3.2 風(fēng)險(xiǎn)規(guī)避對(duì)策

鋼鐵企業(yè)供水保證率為97%，本次提出的3種方案，方案1和方案3的來水條件下，取水口水量豐富，取水能夠得到保障；但方案2即平水期若上小河上游來流量較小，本項(xiàng)目取水有風(fēng)險(xiǎn)。為此，提出以下3個(gè)方面的建議。

3.2.1疏浚河道

平水期長(zhǎng)洲水利樞紐運(yùn)行水位為18.758 m，為保證該水位情況下長(zhǎng)洲水利樞紐的回水達(dá)到取水口，需疏浚取水口至上小河潯江匯合口河段，特別是方案2情況下MIKE模型計(jì)算水深低于0.5 m甚至為0 m的河段。

3.2.2下移取水點(diǎn)

為了錯(cuò)開取水口至上小河潯江匯合口河段中河底高程較高的河段，建議將取水口下移至秋米河與上小河交匯點(diǎn)北面的河段。從方案2情況下MIKE模型計(jì)算水深(圖5)可知，秋米河與上小河交匯點(diǎn)北面至上小河潯江匯合口河段水深均大于1.5 m，可滿足本項(xiàng)目取水需求。

3.2.3平水期上游最小流量

通過MIKE模型試算，上小河上游來流量達(dá)到2.5 m3/s時(shí)，本取水口水位與方案1的水位一致，本取水口取水才有保障。

4 結(jié)論與展望

本研究建立了潯江-秋米河-上小河MIKE21模型，模擬了3種水文情況下取水河段的水位情況，根據(jù)模擬結(jié)果分析了長(zhǎng)洲水利樞紐庫區(qū)范圍內(nèi)潯江支流上小河取水可靠性，計(jì)算結(jié)果能夠?yàn)楣こ毯侠碓O(shè)計(jì)、水行政主管部門決策提供依據(jù)。得出的主要結(jié)論如下。

a)由于本取水口位于長(zhǎng)洲水利樞紐庫區(qū)范圍，枯水年長(zhǎng)洲水利樞紐蓄水水位較高，本項(xiàng)目取水有保障。

b)平水年長(zhǎng)洲水利樞紐需要預(yù)泄洪水，運(yùn)行水位較低，若上小河上游來流量有保證，本項(xiàng)目可取到水。由于上小河上游建有多座攔河壩，平水期如4月份左右取水口處來流量較小時(shí)，本項(xiàng)目取水存在較大風(fēng)險(xiǎn)。

c)為了保障本項(xiàng)目的取水需求，應(yīng)對(duì)平水期河道流量無法滿足項(xiàng)目取水需求的問題，可采用疏浚河道、下移取水口等方法，以提高工程取水保證率，從而規(guī)避可能出現(xiàn)的取水風(fēng)險(xiǎn)。