喻海軍,吳俊峰,吳濱濱,孫宏圖,馬建明,高 佳
(1.中國(guó)水利水電科學(xué)研究院,北京 100038;2.水利部防洪抗旱減災(zāi)工程技術(shù)研究中心,北京 100038;3.河北省水利水電第二勘測(cè)設(shè)計(jì)研究院,河北 石家莊 050021;4.河北省防汛抗旱指揮部辦公室,河北 石家莊 050011)
蓄滯洪區(qū)是我國(guó)江河防洪體系中的重要組成部分,是保障重點(diǎn)城市和地區(qū)防洪安全,減輕洪水災(zāi)害的有效措施[1]。蓄滯洪區(qū)洪水?dāng)?shù)值模擬對(duì)蓄滯洪區(qū)洪水管理、洪水損失評(píng)估以及流域洪水調(diào)度決策等都具有重要的意義[2-4]。近年來(lái),國(guó)內(nèi)外學(xué)者在蓄滯洪區(qū)洪水?dāng)?shù)值模擬方面做了大量研究,取得了豐碩成果。王秀杰等[5]建立了河道與防洪保護(hù)區(qū)的全二維動(dòng)態(tài)耦合模型,實(shí)現(xiàn)了寬淺河道與防洪區(qū)的無(wú)縫銜接,較好地模擬了水流流態(tài)的變化。李大鳴等[6]采用一二維耦合水動(dòng)力學(xué)模型實(shí)現(xiàn)了大清河五洼滯洪區(qū)的實(shí)時(shí)洪水調(diào)度分析。Vorogushyn等[7]采用一維河道、潰堤和二維模型實(shí)現(xiàn)了潰堤洪水在蓄滯洪區(qū)內(nèi)的演進(jìn)模擬,對(duì)德國(guó)Elbe河蓄滯洪區(qū)內(nèi)的洪災(zāi)損失進(jìn)行了評(píng)估。一般來(lái)說(shuō),大型蓄滯洪區(qū)的洪水模擬一般需要構(gòu)建一二維耦合模型[6-7],但往往由于面積大、河流眾多等原因?qū)е履P蜁r(shí)效性比較差,給實(shí)時(shí)預(yù)報(bào)分析帶來(lái)一定的困難,因而十分有必要研究相對(duì)精準(zhǔn)且高效的數(shù)值模擬方法。本文采用基于二維網(wǎng)格邊元上設(shè)置河道的方法構(gòu)建一、二維耦合模型,并選擇海河流域最大的蓄滯洪區(qū)大陸澤寧晉泊蓄滯洪區(qū)為典型應(yīng)用區(qū)域,對(duì)本文建立的一、二維耦合模型進(jìn)行了驗(yàn)證與應(yīng)用,并分析了下滲因素對(duì)蓄滯洪區(qū)行洪的影響。
2.1 一維模型河道洪水采用一維水動(dòng)力學(xué)模型模擬,控制方程為一維圣維南方程組:
式中:q為旁側(cè)入流;Q為流量;B為水面寬度;Z為水位;A為過(guò)水面積;α為動(dòng)量校正系數(shù);Sf為摩阻比降。
采用基于有限體積法的Godunov格式對(duì)上述一維圣維南方程組進(jìn)行離散[8]。
2.2 二維模型地表洪水演進(jìn)采用二維水動(dòng)力學(xué)模型,控制方程采用守恒型的二維淺水方程:
式中:h、u、v、b分別為水深、x和y方向流速、底高程;Sb是底坡項(xiàng);Sox、Soy分別是x和y方向的底坡比降;Sf是摩阻項(xiàng);Sfx、Sfy分別是x和y方面的摩阻比降。
采用基于有限體積法的Godunov格式對(duì)上述二維淺水方程組進(jìn)行離散求解,其中Riemann問(wèn)題采用Roe格式的近似解進(jìn)行計(jì)算,底坡源項(xiàng)采用特征分級(jí)離散,保證模型的守恒性,阻力源項(xiàng)采用隱式離散提高模型的穩(wěn)定性,采用MUSCL空間重構(gòu)和兩步Runge-Kutta法使得模型具有時(shí)間和空間二階精度,所有變量都定義在單元中心[9-10]。
2.3 模型耦合一維河道和二維地表水動(dòng)力模型的水流交互采用邊元耦合的方式(如圖1所示),即在二維網(wǎng)格邊元上設(shè)置和建立一維河網(wǎng)模型,網(wǎng)格劃分時(shí)不考慮河道寬度,因而河道寬度不影響網(wǎng)格單元?jiǎng)澐?,可以極大的提高網(wǎng)格質(zhì)量,從而提高模型演算效率。相對(duì)于傳統(tǒng)分別構(gòu)建一二維模型,河道范圍內(nèi)不劃分網(wǎng)格或者全部采用二維模型計(jì)算,河道處進(jìn)行網(wǎng)格加密等方法,采用邊元設(shè)置河道來(lái)耦合的方式無(wú)論從網(wǎng)格劃分便捷性還是計(jì)算效率來(lái)看都得到了極大的提升,特別適用于在計(jì)算范圍內(nèi)存在較多對(duì)洪水演進(jìn)有影響的中小河道的情形。
圖1 二維網(wǎng)格邊元設(shè)置河道示意
圖2 河道與地面水流交互示意
如圖2所示,假設(shè)某時(shí)刻河道與地表網(wǎng)格水量通過(guò)側(cè)向交換的方式交換的流量為Ql,采用堰流公式近似計(jì)算交換流量的方法如下[11]:
其中:hmax和hmin分別采用下式計(jì)算:
式中:Zr、Zc分別為堰上、下游水位,分別取河道和二維網(wǎng)格單元的水位值;Ze為堰的高程,一般取河堤岸的高程;be為堰的寬度,一般取單元格與河道相連邊的邊長(zhǎng)。
3.1 區(qū)域概況本文選擇典型應(yīng)用區(qū)域大陸澤及寧晉泊蓄滯洪區(qū)位于河北省南部(如圖3所示),海河流域子牙河水系支流滏陽(yáng)河系中游,河北省邢臺(tái)市東北部,是全國(guó)第三大蓄滯洪區(qū),也是海河流域第一大蓄滯洪區(qū)和關(guān)鍵防洪工程,總面積約2041 km2,總滯洪量35.52億m3,涉及邢臺(tái)市9個(gè)縣區(qū),135.27萬(wàn)人。上游河流多處于太行山迎風(fēng)區(qū),源短、坡陡、流急,洪水突發(fā)性、致災(zāi)性強(qiáng)。蓄滯洪區(qū)內(nèi)分布有眾多中小河流及堤壩,且堰閘、分區(qū)滯洪等調(diào)度規(guī)則復(fù)雜,導(dǎo)致一維河道與二維滯洪區(qū)洪水演進(jìn)過(guò)程復(fù)雜。大陸澤、寧晉泊蓄滯洪區(qū)按50年一遇洪水設(shè)計(jì),從建國(guó)后至今60多年的蓄滯洪區(qū)運(yùn)用情況來(lái)看,大陸澤實(shí)際運(yùn)用4次,最高蓄洪水位33.41 m,最大蓄洪水量10.40億m3,蓄洪歷時(shí)57 d;寧晉泊實(shí)際運(yùn)用3次,最高蓄洪水位30.70 m,最大蓄洪水量29.27億m3,蓄洪歷時(shí)60 d。
3.2 模型構(gòu)建
(1)網(wǎng)格劃分。采用非結(jié)構(gòu)四邊形網(wǎng)格對(duì)整個(gè)研究區(qū)域進(jìn)行網(wǎng)格剖分,將區(qū)域內(nèi)高速公路、鐵路和河流作為內(nèi)部約束邊界,最終共劃分成約兩萬(wàn)個(gè)網(wǎng)格。在網(wǎng)格插值時(shí),將堤防、高速公路、鐵路設(shè)置成線形阻水建筑物在網(wǎng)格邊元上予以考慮,將測(cè)量的實(shí)際高程作為阻水建筑物的高度,采用堰流公式計(jì)算建筑物兩側(cè)的過(guò)流水量。河道糙率取0.025,其它區(qū)域糙率取值范圍在0.040~0.050之間?;诰W(wǎng)格邊元設(shè)置的河道分布如圖4(a)所示,網(wǎng)格局部放大如圖4(b)所示。
圖3 大陸澤寧晉泊蓄滯洪區(qū)位置示意
圖4 研究區(qū)域網(wǎng)格示意
(2)邊界條件。研究區(qū)域內(nèi)一共包含滏陽(yáng)河、留壘河、洺河、南澧河、順?biāo)樱ㄆ呃锖樱?、牛尾河、馬河(白馬河、小馬河和李陽(yáng)河)、泜河、午河(泲河)、北沙河(槐河)、小漳河和洨河等12條中小河道。在實(shí)際模擬計(jì)算時(shí),河道上游入流采用流量邊界條件(典型洪水過(guò)程如圖5所示),下游艾辛莊樞紐的出流采用水位流量關(guān)系邊界(如圖6所示)。
(3)分區(qū)滯洪。大陸澤及寧晉泊蓄滯洪區(qū)被滏陽(yáng)河右堤和小漳河右堤劃成3個(gè)分區(qū)(如圖7所示),在實(shí)際調(diào)度運(yùn)行中,當(dāng)分區(qū)1的水位達(dá)到29.4 m時(shí),由小南堤分洪口門(mén)向分區(qū)二分洪,當(dāng)分區(qū)二的水位超過(guò)29.4 m時(shí)由小漳河右堤分洪口門(mén)向分區(qū)三分洪。模型構(gòu)建時(shí)考慮了蓄滯洪區(qū)分區(qū)滯洪的調(diào)度規(guī)則,兩個(gè)分洪口門(mén)則采用潰口的方式進(jìn)行概化考慮,達(dá)到設(shè)定的條件(h>29.4 m),自動(dòng)進(jìn)行分洪。
圖5 典型洪水過(guò)程(洺河20年和50年一遇設(shè)計(jì)洪水過(guò)程)
圖6 艾辛莊樞紐水位流量關(guān)系
圖7 分區(qū)滯洪示意
(4)模型驗(yàn)證。由于缺乏蓄滯洪區(qū)運(yùn)用的實(shí)測(cè)數(shù)據(jù),無(wú)法對(duì)構(gòu)建的模型進(jìn)行較為準(zhǔn)確的率定驗(yàn)證,因而將模型計(jì)算結(jié)果與《子牙河系防洪規(guī)劃》(2008年)以及《大陸澤、寧晉泊蓄滯洪區(qū)建設(shè)與管理可行性研究報(bào)告》(2014年)中的模擬結(jié)果進(jìn)行比對(duì)分析,以對(duì)模型結(jié)果合理性進(jìn)行檢驗(yàn)。現(xiàn)狀條件下,計(jì)算50年一遇洪水,采用位于蓄滯洪區(qū)上游的環(huán)水村(位于邢臺(tái)市任縣)以及下游艾辛莊兩處(位置示意見(jiàn)圖7)的最高洪水位進(jìn)行對(duì)比,結(jié)果如表1所示。總體上,本次構(gòu)建的模型與其它研究結(jié)果有一定差別,但水位差總體控制在0.15 m以?xún)?nèi),表明模型結(jié)果是可信的,差別主要來(lái)源可能是模型算法和概化方法不完全一樣。
3.3 模型應(yīng)用大陸澤及寧晉泊蓄滯洪區(qū)位于海河流域,屬于干旱半干旱地區(qū),近幾十年由于地下水長(zhǎng)期處于超采狀態(tài),地下水位持續(xù)下降,導(dǎo)致河道和地表行洪時(shí)下滲比較強(qiáng)烈,對(duì)洪水演進(jìn)(如洪峰、洪水淹沒(méi)時(shí)間等)有著重要的影響。由于以往針對(duì)大陸澤及寧晉泊蓄滯洪區(qū)的研究成果多沒(méi)有考慮下滲的影響[3],從洪水調(diào)度管理的角度來(lái)看,不考慮下滲屬于偏安全的做法,存在一定的合理性,但由于蓄滯洪區(qū)一旦啟用,蓄水時(shí)間長(zhǎng),下滲量就相對(duì)比較大,其影響不應(yīng)被忽略。本文著重對(duì)比分析了考慮和不考慮下滲情況下,蓄滯洪區(qū)洪水演進(jìn)情況的區(qū)別,下滲率取值參考程亮等[12]在海河流域的研究成果,穩(wěn)定下滲率取值12.3 mm/h。
圖8給出了考慮和不考慮下滲兩種情況下遭遇50年一遇洪水淹沒(méi)范圍。從圖8中可以看出,在t=50 h時(shí)上下游已經(jīng)有部分洪水從河道中溢出,但此時(shí)淹沒(méi)深度均較小,隨著時(shí)間增長(zhǎng),淹沒(méi)范圍逐漸增大??傮w上看,不考慮下滲時(shí),50年一遇洪水需要啟用全部的三個(gè)分洪區(qū),而考慮下滲時(shí)則只需啟用第一個(gè)分洪區(qū),差別比較明顯。分析同一時(shí)刻考慮和不考慮下滲兩種情況,可以明顯看出同一時(shí)刻不考慮下滲淹沒(méi)范圍和淹沒(méi)深度均大于考慮下滲時(shí),而考慮長(zhǎng)歷時(shí)洪水下滲對(duì)洪水淹沒(méi)范圍和淹沒(méi)深度都有較大影響。
表1 模擬結(jié)果對(duì)比
圖8 50年一遇不同時(shí)刻淹沒(méi)圖
圖9給出了20年一遇和50年一遇環(huán)水村附近洪水位變化過(guò)程。從圖9可以看出,考慮下滲時(shí)的洪峰值明顯小于不考慮下滲時(shí),數(shù)值上分別減少了0.42 m和0.41 m,另外,考慮下滲時(shí)的退水時(shí)間與不考慮下滲時(shí)相比分別減少了195 h和296 h。從整個(gè)蓄滯洪區(qū)來(lái)看,以遭遇50年一遇洪水為例,不考慮下滲時(shí)蓄滯洪區(qū)需要約60 d的時(shí)間將洪水基本排盡,考慮下滲時(shí),則只需要約36 d時(shí)間,減少了近24 d時(shí)間。
圖9 環(huán)水村洪水變化過(guò)程
(1)基于在二維邊元設(shè)置河道方法,建立了大陸澤及寧晉泊蓄滯洪區(qū)一二維耦合的地表和河道洪水?dāng)?shù)值模型,在與其它模型計(jì)算結(jié)果進(jìn)行對(duì)比驗(yàn)證時(shí),水位差總體控制在0.15 m以?xún)?nèi),表明該模型結(jié)果是可信的,在此基礎(chǔ)上對(duì)20年一遇和50年一遇洪水進(jìn)行了模擬與分析。應(yīng)用結(jié)果表明基于二維邊元設(shè)置河道的一二維耦合模型能夠很好的應(yīng)用于中小河流較多的蓄滯洪區(qū),具有廣闊的應(yīng)用前景。
(2)對(duì)比分析地表下滲對(duì)蓄滯洪區(qū)洪水的影響可見(jiàn):考慮下滲作用時(shí),蓄滯洪區(qū)50年一遇洪水只需要啟用1個(gè)分區(qū),而不考慮時(shí)則需要啟用3個(gè);考慮下滲時(shí)退水時(shí)間與不考慮下滲時(shí)相比也相應(yīng)減少了近24 d。表明干旱半干旱地區(qū)下滲對(duì)蓄滯洪區(qū)行洪具有較大的影響,在蓄滯洪區(qū)實(shí)際調(diào)度應(yīng)用時(shí),應(yīng)該考慮下滲因素對(duì)洪量減少的影響,以較好發(fā)揮蓄滯洪區(qū)的效用。