獻(xiàn)縣泛區(qū)洪水演進(jìn)數(shù)學(xué)模型研究

李大鳴，孟 政，李彥卿，陳 碩，卜世龍，王騰飛

(1.天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072； 2.生態(tài)環(huán)境部海河流域北海海域生態(tài)環(huán)境監(jiān)督管理局 生態(tài)環(huán)境監(jiān)測與科學(xué)研究中心，天津 300170)

1 研究背景

近年來，極端天氣頻繁出現(xiàn)，尤其以暴雨最為顯著。而為防止暴雨洪水對下游一些重點(diǎn)城市及設(shè)施產(chǎn)生危害，往往會在河道兩岸低洼地帶設(shè)置蓄滯洪區(qū)，讓多余的洪水暫存在這里。隨著近些年來蓄滯洪區(qū)內(nèi)居民區(qū)及農(nóng)田的發(fā)展，其下墊面已經(jīng)發(fā)生了很大的變化。利用當(dāng)前資料模擬蓄滯洪區(qū)的洪水演進(jìn)，對蓄滯洪區(qū)風(fēng)險(xiǎn)圖的制作和蓄滯洪區(qū)的建設(shè)很重要。國外對洪水演進(jìn)大都采用馬斯京根法(Muskingum)：Farahani等[1]利用四參數(shù)非線性Muskingum和Shark算法，提出了一種新的洪水演進(jìn)方法；Saeed等[2]利用三參數(shù)的Muskingum模型和改進(jìn)的Bat算法進(jìn)行河段洪水演進(jìn)；Retsinis等[3]采用水文及水力學(xué)方法在棱形渠道中進(jìn)行了動態(tài)洪水演進(jìn)模擬。國內(nèi)較多地采用水文水力學(xué)耦合方法對洪水演進(jìn)進(jìn)行研究：蘆云峰[4]利用水文水動力耦合的洪水演進(jìn)模型探討了水庫調(diào)度對洪水水流結(jié)構(gòu)和斷面過流能力的影響；于富強(qiáng)等[5]耦合水文模型和二維水動力模型對泉州梅溪流域進(jìn)行了洪水演進(jìn)模擬；于琛等[6]基于有限元方法(DEM)分析的洪水淹沒算法，模擬了資料不完備狀態(tài)下的洪水；董柏良等[7]利用物理試驗(yàn)?zāi)Ｐ停治隽私ㄖ锩芏群途G化帶布置對洪水演進(jìn)的影響；程坤等[8]用MIKE21軟件模擬了西藏藏木水電站潰壩洪水演進(jìn)；魏博文等[9]對Mike11一維河道模型和Mike21二維洪泛區(qū)模型進(jìn)行了動態(tài)耦合，模擬了丘陵地區(qū)叉網(wǎng)式河流段的洪水；吳儀等[10]基于HEC-RAC模型的二維模塊高效地模擬了胖頭泡蓄滯洪區(qū)的洪水；王棟等[11]基于EFDC水動力模塊建立了清河下游洪水演進(jìn)水動力模型；李大鳴等[12]利用一二維耦合的洪水演進(jìn)數(shù)學(xué)模型模擬了西三洼洪水演進(jìn)，并制作了洪水風(fēng)險(xiǎn)圖。

前人已經(jīng)利用水文、水力學(xué)或水文水力學(xué)耦合方法對洪水進(jìn)行了模擬，取得了不錯(cuò)的模擬結(jié)果。本文在模型建立過程中充分考慮鐵路、堤防、公路等的影響，將網(wǎng)格通道的類型劃分得較為細(xì)致，使其結(jié)果更加貼合實(shí)際情況，這為防洪規(guī)劃與決策提供了更加可信的依據(jù)。

2 數(shù)學(xué)模型理論基礎(chǔ)

2.1 基本控制方程

連續(xù)方程：

(1)

動量方程：

(2)

(3)

式中：H為水深；Z為水位，Z=Z0+H，Z0為底高程；q為源匯項(xiàng)；u、v分別為x、y方向的平均流速；M、N分別為x、y方向上的單寬流量，且M=Hu，N=Hv；n為糙率；g為重力加速度。

2.2 有限體積離散

本文采用有限體積法對控制方程進(jìn)行離散，中點(diǎn)的通量可用中心格式(如取相鄰兩格子形心處通量的平均)或逆風(fēng)格式確定。

2.2.1 連續(xù)方程的離散

將方程(1)改寫成矢量形式，將其在控制體內(nèi)進(jìn)行積分，則方程(1)可離散為

(4)

2.2.2 動量方程的離散

考慮到蓄滯洪區(qū)不同的地形情況及植被、公路、鐵路等建筑物的影響，本模型將控制體的通道按不同情況進(jìn)行分類，概化為地面型通道、河道型通道和其他通道進(jìn)行處理，并給所有通道附加特征信息，以相應(yīng)的水力學(xué)公式進(jìn)行通量計(jì)算。

(1)地面型通道，即通道兩側(cè)單元為陸地地面，且通道上沒有堤防等阻水建筑物?？紤]到滯洪區(qū)內(nèi)的地形起伏不大，地面洪水演進(jìn)主要受到重力和阻力的作用，忽略加速度項(xiàng)，即只保留方程(2)和方程(3)中的相應(yīng)項(xiàng)，利用差分方法離散得到地面型通道的動量離散方程，即

(5)

(2)河道型通道，即通道兩側(cè)網(wǎng)格均為河道型網(wǎng)格，動量方程中保留局地加速度項(xiàng)、重力項(xiàng)和阻力項(xiàng)，即保留方程(2)和方程(3)中的相應(yīng)項(xiàng)，利用差分方法離散得到河道型通道的動量離散方程，即

(6)

(3)其他類型通道，即對于滯洪區(qū)內(nèi)寬度較小的河流，不便于將其劃分成獨(dú)立的單元網(wǎng)格，也不能忽略不計(jì)，具體的計(jì)算方法如下：

①特殊通道與兩側(cè)網(wǎng)格之間的流量計(jì)算，采用寬頂堰流公式，即

(7)

式中：σs為淹沒系數(shù)；m為流量系數(shù)；b為通道寬度。

②動量方程離散同河道型通道離散方法，沿河道的單寬流量Qs采用河道型通道式(6)計(jì)算。

③采用連續(xù)方程來計(jì)算特殊單元的水位，即

(8)

(4)連續(xù)堤或缺口堤：當(dāng)計(jì)算區(qū)域有許多高于地面的阻水建筑物，如堤防、鐵路、公路等時(shí)，其流量采用寬頂堰溢流公式來計(jì)算，離散后得到

(9)

2.3 一維河道模型理論

明渠非恒定流方程采用水位-流量關(guān)系表示，即

(10)

式中：Z為水位；t為時(shí)間坐標(biāo)；B為?；瘜挾?；Q為流量；x為沿程斷面坐標(biāo)；q為旁側(cè)出流單寬流量；A為過流面積；C為謝才系數(shù)；R為水力半徑。

一維模型中斷面所在位置視為單元中心并計(jì)算守恒變量，斷面上下游Δx范圍內(nèi)的守恒量全部由斷面的守恒量代替，如圖1所示。斷面間距的中點(diǎn)作為相鄰單元的交界面，并在此處計(jì)算數(shù)值通量。

圖1 一維模型離散格式Fig.1 Discretization scheme of one-dimensional model

根據(jù)圖1所示，一維模型控制方程的離散形式可以表示為

(11)

式中：U1為一維模型的守恒向量；F1為一維模型的通量向量；S1為一維模型的源項(xiàng)向量。其中：

(12)

式中：ql表示源匯項(xiàng)，通常為旁側(cè)入流；So為底坡坡降；Sf為摩阻坡降；I1為靜水壓力項(xiàng)；I2為沿縱向上斷面寬度變化引起的壓力。

其中：

(13)

式中：b(x，δ)表示深度δ位置的過水?dāng)嗝鎸挾龋籬為斷面水深；zb為斷面底高程。

3 獻(xiàn)縣泛區(qū)洪水演進(jìn)模型的建立

3.1 研究區(qū)域概況

獻(xiàn)縣泛區(qū)位于河北省衡水市和滄州市交界處，在衡水饒陽縣與滄州獻(xiàn)縣之間。獻(xiàn)縣泛區(qū)內(nèi)主要有滹沱河下游河道穿過，泛區(qū)西起饒陽縣大齊村西肅臨公路，東至獻(xiàn)縣樞紐，南北以滹沱河南大堤、北大堤為界。獻(xiàn)縣泛區(qū)主要承納并宣泄滹沱河及滹滏區(qū)間的洪瀝水，獻(xiàn)縣泛區(qū)內(nèi)地勢低洼，從西南向東北傾斜，坡降在1/4 000～1/6 000之間，最低點(diǎn)張村鄉(xiāng)12.5 m(1985國家高程基準(zhǔn)，下同)，最高點(diǎn)臨河鄉(xiāng)19.1 m。獻(xiàn)縣泛區(qū)東西長為22.5 km，南北寬12.5 km，泛區(qū)面積約為331 km2。泛區(qū)由滹沱河為深槽的行洪道分割，分為南泛區(qū)和北泛區(qū)。獻(xiàn)縣泛區(qū)平面示意圖見圖2。

圖2 獻(xiàn)縣泛區(qū)平面示意圖Fig.2 Map of flood plain in Xianxian County

為研究超標(biāo)準(zhǔn)洪水演進(jìn)過程，選擇圖3為超標(biāo)準(zhǔn)洪水的研究范圍。研究范圍東起獻(xiàn)縣樞紐，西到西李莊；南迄石黃高速公路，北至滹沱河北大堤。東西長度約為67 km，南北寬度約為35 km。研究范圍包含滹沱河西李莊至獻(xiàn)縣河段及兩岸行洪河床、獻(xiàn)縣北泛區(qū)、獻(xiàn)縣南泛區(qū)、滹沱河南大堤以南與石黃高速公路以北的低洼地區(qū)，面積約為1 775 km2。在研究范圍內(nèi)超標(biāo)洪水分洪口門有7個(gè)，分別是大齊、東草蘆、張弛、耿各莊、故城、羅屯和東呈干分洪口門。各分洪口門的主要分洪控制條件見表1。

在表1中，1至4號扒口按順序啟用，控制水位為扒口附近的滹沱河河道中水位；5至7號扒口啟用控制水位為姚莊處滹沱河河道斷面水位?？刂屏髁恐饕渣S壁莊水庫下泄流量為標(biāo)準(zhǔn)。

表1 分洪控制條件Table 1 Control conditions for flood diversion

3.2 模型計(jì)算網(wǎng)格的剖分

計(jì)算區(qū)域網(wǎng)格剖分：模型范圍內(nèi)滹沱河從西李莊至獻(xiàn)縣樞紐大約75 km，劃分為75個(gè)斷面，平均每千米一個(gè)斷面，河道平均糙率為0.024。其他區(qū)域網(wǎng)格剖分情況如表2所示，計(jì)算模型網(wǎng)格分布如圖3所示。

圖3 計(jì)算模型網(wǎng)格Fig.3 Computational grids

表2 區(qū)域網(wǎng)格剖分情況Table 2 Regional meshing subdivision

3.3 模型的邊界條件

模型的邊界條件：滹沱河入流邊界主要是黃壁莊水庫的下泄水量。10 a一遇下泄水量為18.58億m3，最大流量為800 m3/s；20 a一遇下泄水量為25.66億m3，最大流量為2 990 m3/s；50 a一遇下泄水量為34.46億m3，最大流量為3 000 m3/s。模型邊界滹沱河洪水入流過程見圖4。

圖4 滹沱河各頻率洪水入流過程Fig.4 Process of flood inflow of varied frequency of Hutuo River

泛區(qū)退水口門控制工程是獻(xiàn)縣樞紐，建于1967年，位于獻(xiàn)縣城西北5 km，是滹沱河、滏陽河和滏陽新河三河匯流處，該樞紐由子牙河節(jié)制閘、子牙新河主槽進(jìn)洪閘、灘地進(jìn)洪堰3座建筑物組成?？刂屏饔蛎娣e為45 500 km2。子牙河節(jié)制閘共3孔，設(shè)計(jì)流量為600 m3/s，子牙新河進(jìn)洪閘共6孔，設(shè)計(jì)流量為943 m3/s，灘地溢洪堰設(shè)計(jì)流量為5 050 m3/s。獻(xiàn)縣樞紐水位-流量關(guān)系見圖5。

圖5 獻(xiàn)縣樞紐水位-流量關(guān)系Fig.5 Stage-discharge relation at Xianxian County

3.4 模型的驗(yàn)證

采用滹沱河現(xiàn)狀行洪能力的成果表對模型計(jì)算的水位進(jìn)行驗(yàn)證，如表3所示。對比發(fā)現(xiàn)在相同流量情況下，同位置處水位的絕對誤差都在0.1 m以下，說明模型較為準(zhǔn)確，可以用來進(jìn)行獻(xiàn)縣蓄滯洪區(qū)的洪水演進(jìn)。

表3 同流量水位對比Table 3 Comparison of water level with equal flow rate

4 洪水演進(jìn)模擬結(jié)果

4.1 10 a一遇洪水演進(jìn)模擬

在黃壁莊水庫以10 a一遇洪水下泄時(shí)，洪峰最大流量為800 m3/s。根據(jù)衡水市分洪口門運(yùn)用方案，不啟用獻(xiàn)縣泛區(qū)，通過滹沱河河道行洪。在河道無水條件下，10 a一遇洪水沿河道演進(jìn)，最大洪峰流量為800 m3/s，洪峰流量持續(xù)時(shí)間約為569 h，河道下游邊界最大泄洪量也為800 m3/s。在進(jìn)洪量與泄洪量持平的情況下，河道水面線保持不變。圖6為河道主要斷面的水位歷時(shí)變化，在進(jìn)洪流量與泄洪流量持平的行洪時(shí)間段水位不變。洪水前鋒大約在18 h后到達(dá)姚莊，大約35 h后到達(dá)獻(xiàn)縣樞紐。圖7為沿河道水面線與堤高比較，說明堤高滿足10 a一遇洪水行洪條件。

圖6 滹沱河10 a一遇洪水主要斷面的水位歷時(shí)變化Fig.6 Time-histories of water level of main sections of Hutuo River in the presence of 10-year event flood

圖7 滹沱河10 a一遇洪水沿程高水位與堤高比較Fig.7 Comparison between high water level and dike height of Hutuo River in the presence of 10-year event flood

4.2 重現(xiàn)期20、50 a一遇洪水演進(jìn)模擬

在黃壁莊水庫以20、50 a一遇洪水下泄時(shí)，根據(jù)衡水市分洪口門運(yùn)用方案，依次運(yùn)用大齊、東草蘆、張弛和耿各莊口門分洪，啟用獻(xiàn)縣泛區(qū)。除分洪口門在達(dá)到運(yùn)用條件時(shí)啟用外，當(dāng)水位高于滹沱河堤防時(shí)，會形成過流堤防。20、50 a一遇洪水的洪峰最大流量、3 000 m3/s以上流量持續(xù)時(shí)間、分洪口啟用時(shí)間、洪峰前導(dǎo)波出現(xiàn)時(shí)間以及南北泛區(qū)水位達(dá)到最大值時(shí)間等具體情況如表4所示。在河道無水條件下，20 a一遇洪水在85 h前沿河道演進(jìn)，最大洪峰流量為3 299 m3/s，洪峰流量持續(xù)時(shí)間約為7 h，50 a 一遇洪水在77 h前沿河道演進(jìn)，河道下游邊界最大泄洪量為800 m3/s。

表4 20、50 a一遇洪水情況Table 4 Situations of floods with 20-year and 50-year return periods

圖8為20、50 a一遇時(shí)河道主要斷面的水位歷時(shí)變化，在進(jìn)洪量較大時(shí)，河道上水位增高；在大齊、東草蘆、張馳和耿各莊分洪后水位明顯回落；在進(jìn)洪量減小后，大齊、東草蘆、張馳水位再次明顯回落，耿各莊水位回落不明顯。

圖8 河道主要斷面的水位歷時(shí)變化Fig.8 Time-histories of water level in major sections of river channel

表5為滹沱河主要位置的堤高與10、20、50 a一遇洪水高水位比較，10 a一遇時(shí)，左右堤超高均為2 m以上；20、50 a一遇時(shí)，姚莊下游洪水高水位超出堤高，且在50 a一遇情況下，姚莊下游水位超出堤高較大。圖9為沿河道水面線與堤高比較，說明姚莊下游堤高不滿足20、50 a一遇洪水行洪條件，右堤首先會出現(xiàn)漫堤，左堤在靠近獻(xiàn)縣樞紐附近會出現(xiàn)漫堤。

表5 滹沱河主要位置的堤高與10、20、50 a一遇洪水高水位比較Table 5 Comparison between the dike height at major positions of Hutuo River and the high water level in floods with 10-year，20-year，and 50-year return periods

圖9 沿河道水面線與堤高比較Fig.9 Comparison between river surface line and dike height

5 結(jié) 論

本文以一二維非恒定流方程作為基本方程，利用有限體積法對方程離散求解，在充分考慮了地形的狀況下，更為細(xì)致地概化了控制體間的通道，建立洪水演進(jìn)模型。最終研究成果為：在10 a一遇洪水情況下，不啟用泛區(qū)蓄洪，最大進(jìn)泄洪量都為800 m3/s 時(shí)，河道水面線保持不變；主要位置堤高均比洪水高水位高2 m以上，說明河堤滿足行洪要求。重現(xiàn)期為20 a一遇洪水時(shí)，在86 h后依次運(yùn)用大齊、東草蘆、張弛和耿各莊口門分洪，啟用獻(xiàn)縣泛區(qū)；大約在300 h獻(xiàn)縣泛區(qū)水位達(dá)到最大值；在20 a一遇高水位時(shí)，姚莊下游河道水位超過堤高，右堤首先會出現(xiàn)漫堤，左堤在靠近獻(xiàn)縣樞紐附近出現(xiàn)漫堤。重現(xiàn)期為50 a一遇洪水時(shí)，大約在280 h獻(xiàn)縣泛區(qū)水位就達(dá)到最大值；姚莊下游水位超出堤高較大，說明下游堤高不滿足50 a一遇洪水行洪條件?？梢愿鶕?jù)本文研究成果，繪制出不同重現(xiàn)期的風(fēng)險(xiǎn)圖；同時(shí)結(jié)合社會經(jīng)濟(jì)資料，進(jìn)一步評估洪水災(zāi)害造成的經(jīng)濟(jì)損失。