基于一、二維耦合水動力模型的贛西聯(lián)圩潰堤洪水風(fēng)險分析

陳俊鴻，劉小龍，王 崗，彭思韋,，張慶梓,，陳煉鋼

(1. 武漢大學(xué)，武漢 430072；2. 東南大學(xué)，南京 210096；3. 河海大學(xué)，南京 210098；4. 南京水利科學(xué)研究院，南京 210029)

贛西聯(lián)圩防洪保護區(qū)位于贛江主支下游西岸，北臨鄱陽湖，保護區(qū)內(nèi)水系繁多，同時水文條件錯綜復(fù)雜。加上贛西聯(lián)圩設(shè)計標(biāo)準(zhǔn)偏低，保護區(qū)內(nèi)地形平坦，地勢低洼，汛期同時受贛江和錦江洪水的頂托影響，圩內(nèi)洪澇災(zāi)害嚴(yán)重。近60年贛西聯(lián)圩曾多次發(fā)生潰決，給圩區(qū)人民生產(chǎn)、生活造成嚴(yán)重的經(jīng)濟損失[1]。因此進行贛西聯(lián)圩洪水風(fēng)險分析，對于防汛指揮部門做出風(fēng)險決策，保護當(dāng)?shù)鼐用裆敭a(chǎn)安全，具有十分現(xiàn)實意義。

近年來，隨著計算機性能的不斷提高和數(shù)值計算技術(shù)的高速發(fā)展，數(shù)值模擬逐漸成為研究洪水演進的主要手段，主要可分為水文學(xué)方法和水動力學(xué)方法。水文學(xué)方法是以水量平衡方程和槽蓄關(guān)系為基礎(chǔ)，反映水流在河道中的推移和坦化[2,3]。該類方法求解簡單方便，但參數(shù)的率定需大量的實測資料，模型適用性差。而水動力學(xué)方法物理意義明確，適用性強，按參數(shù)的空間分布特性，其數(shù)學(xué)模型可分為一維、二維和三維模型。對于大型流域而言，由于其主要的空間尺度是縱向，一般從整體上采用一維非恒定流數(shù)學(xué)模型模擬水流的演進和運動規(guī)律[4]。而下游洪泛區(qū)多處于平原地區(qū)，洪水傳播具有較強的平面二維特征，因此采用平面二維非恒定流數(shù)學(xué)模型模擬洪水在下游洪泛區(qū)的演進過程[5-7]。在此基礎(chǔ)上，將一維模型和二維模型耦合后，模擬一維河道洪水潰堤后，再模擬洪水在洪泛區(qū)詳細的二維洪水演進過程，具有更高的精度和明顯的優(yōu)越性，已被廣泛應(yīng)用于潰堤洪水模擬中[8-11]。本文將建立一、二維耦合洪水演進模型，一維河道和二維洪泛區(qū)水流運動分別采用四點隱式Preissmann格式以及中心格式的有限體積法進行離散求解，模擬極端洪水過程獲取相關(guān)洪水風(fēng)險要素并進行洪水風(fēng)險分析與評估。

1 模型構(gòu)建原理

1.1 河道一維模型

一維河道(河網(wǎng))的洪水運動用Saint-Venant方程組描述：

(2)

式中：x、t分別為河長、時間；A為過水面積；Q為斷面流量；z為斷面水位；α為動量修正系數(shù)；k為流量模數(shù)；q為旁側(cè)入流；vx為入流沿水流方向的速度分量。

上、下游邊界的控制條件一般采用水位過程控制、流量過程控制、流量與水位關(guān)系控制等形式。對Saint-Venant方程式(1)和式(2)采用四點隱式的Preissmann格式進行離散求解：

(4)

f=(fni+1+fni+fni+1+fni)/4

(5)

式中：θ為權(quán)重系數(shù)，可取0.5<θ<1.0。

1.2 洪泛區(qū)二維模型

對于寬淺河道和洪泛區(qū)內(nèi)大范圍的潰堤洪水運動，可以采用淺水方程來描述，控制方程的守恒形式如下:

(6)

(7)

F=Ei+Gj

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

式中：h為水深；t為時間變量；u、v分別為x、y方向的流速；Sox、Soy分別為x、y方向的坡底源項；Sfx、Sfy分別為x、y方向的摩阻項；g為重力加速度；z為底高程；n為曼寧系數(shù)。

本文采用了中心格式的有限體積法對控制方程式(6)進行離散，控制體取為任意的三角形單元，應(yīng)用格林公式可得如下形式：

(16)

式中：Ai為i單元的計算面積；Γi為i控制體的邊界條件；n為邊界的外法線方向單位向量；F為界面通量。

1.3 模型耦合

將一維河網(wǎng)模型與二維洪水模型進行耦合，構(gòu)建最終的洪水分析模型。選定側(cè)堰流公式來實現(xiàn)潰口上下游水流信息的交互，具體為：

(17)

式中：Q為側(cè)堰流量；H為側(cè)堰首端河渠斷面水深；b為側(cè)堰堰寬；Cd為流量系數(shù)。當(dāng)建立潰堤水流模型時，可根據(jù)湍流方程推求流量系數(shù)的近似方程：

(18)

式中：m為模塊限制系數(shù)；d為平均水深；DX為水流源匯距離；n為曼寧系數(shù)。

將一維模型和二維模型耦合后，模型系統(tǒng)可以自動分析分洪對河道水位的影響，以及河道水位下降反過來對洪泛區(qū)分洪量的影響，提高了模擬精度和可靠性。

2 贛西聯(lián)圩洪水風(fēng)險分析

2.1 區(qū)域概況

贛西聯(lián)圩防洪保護區(qū)(見圖1)位于贛江主支下游西岸，北臨鄱陽湖，東隔主支與廿四聯(lián)圩相望，西南為濱湖崗地，圩堤南起新建縣樵舍，北過高棠分支，西過青山，至大唐西崗地，圩區(qū)內(nèi)設(shè)有方洲斜塘蓄滯洪區(qū)。保護區(qū)內(nèi)主要河流為鐵河，貫穿整個贛西聯(lián)圩防洪保護區(qū)，西部有螞蟻河支流。保護區(qū)地勢西高東低，西邊為低山丘陵，東邊為贛江大堤。外洪堤線總長48.002 km，圩堤內(nèi)有樵舍、象山、鐵河、大塘、金橋5個鄉(xiāng)鎮(zhèn)，圩內(nèi)總集水面積241.08 km2，保護面積114.58 km2，保護耕地約9 900 hm2，保護人口12.06萬人。贛西聯(lián)圩所處區(qū)域?qū)賮啛釒Ъ撅L(fēng)濕潤氣候區(qū)，氣候溫和，雨量充沛，降水量年際變化較大、年內(nèi)分配不均，多年平均降雨量為1 530.3 mm，4-6月降水量占全年48.3%。多年平均蒸發(fā)量為954.9 mm，實測多年平均氣溫17.2 ℃，多年平均無霜期約250 d，多年平均風(fēng)速2.95 m/s，多年平均濕度80%。

2.2 地形概化

保護區(qū)內(nèi)河道斷面是模型計算的最基本的單元，斷面數(shù)據(jù)的準(zhǔn)確性直接影響到模型計算結(jié)果的精確程度。對有實測大斷面資料的河道，直接將實測的斷面數(shù)據(jù)，包括起點距、相應(yīng)點高程、糙率以及斷面間距等信息導(dǎo)入模型，自動生成斷面及其之間的河道連接。對無實測斷面的河道，采用地形文件來生成河道大斷面。根據(jù)地形數(shù)據(jù)的具體情況，以0.5 km左右為斷面間距，由上游向下游逐一在地面模型上切割出大斷面。斷面的切取遵循3個原則：一是要覆蓋河道兩岸的堤防，如果沒有堤防，則將斷面延伸到兩岸高地；二是斷面主河槽的方向要垂直于河道，而對于漫灘地區(qū)，則需要垂直于具體的水流流向；三避開高程線相對稀疏的地區(qū)。

圖1 贛西聯(lián)圩防洪保護區(qū)示意圖Fig.1 Schematic diagram of Ganxi dyke flood control protected area

對贛西聯(lián)圩防洪保護區(qū)計算區(qū)域采用非結(jié)構(gòu)不規(guī)則網(wǎng)格進行網(wǎng)格劃分，使網(wǎng)格的大小隨地形地勢和阻水建筑物的分布靈活確定，而且盡可能地將影響水流的阻水建筑物作為網(wǎng)格邊界，充分反映計算域的特征。但是，必要的時候?qū)ΡＷo區(qū)內(nèi)的一些典型的線性阻水建筑物，如堤防、公路等，經(jīng)合理概化，并對網(wǎng)格適當(dāng)加密，在二維地形中充分反映其特征。研究區(qū)域部分?jǐn)嗝娓呕途W(wǎng)格概化如圖2所示。

圖2 斷面網(wǎng)格概化圖Fig.2 Generalized cross-section and topography

2.3 潰口設(shè)置與邊界條件

根據(jù)歷史潰口及險工險段分析，結(jié)合目前河道堤防存在的險工隱患段與保護對象重要性，選擇易發(fā)生潰口的位置為：永建洲、李家垾(見圖1)。所選潰口位置分別位于贛西聯(lián)圩防洪保護區(qū)上、下位置，能夠反應(yīng)不同潰口位置對洪水淹沒時空分布的影響。贛江上游邊界為外洲站設(shè)計或?qū)嶋H流量過程，下游邊界為鄱陽湖水位過程。利用外洲1949-2013年65 a實測洪峰流量資料，再加上1901、1924年2 a的調(diào)查特大洪水資料構(gòu)成頻率分析所需的洪峰流量系列，采用年最大值法選擇實測洪峰流量系列，以此為基礎(chǔ)進行洪峰流量頻率分析，選擇2010年洪水為典型洪水，按照洪峰控制的同倍比放大法得到外洲站20年一遇設(shè)計洪水，洪水過程線如圖3所示。本文根據(jù)贛西聯(lián)圩防洪保護區(qū)堤防防洪標(biāo)準(zhǔn)及歷史上曾發(fā)生的較大洪水，選擇計算方案為永建洲、李家垾贛江20年一遇設(shè)計洪水湖口2010年實測水位過程，具體方案名稱見表1。

圖3 外洲站設(shè)計洪水過程線Fig.3 Designed flood hydrograph in Waizhou Station

方案編號潰口名稱計算方案1永建洲贛江20年一遇設(shè)計洪水湖口2010年實測水位過程2李家垾贛江20年一遇設(shè)計洪水湖口2010年實測水位過程

2.4 糙率確定

一維水動力學(xué)模型參數(shù)主要為河道糙率。糙率是表征河道底部、岸坡影響水流阻力的綜合系數(shù)，是水力計算的重要靈敏參數(shù)，也是水動力數(shù)學(xué)模型中最重要的參數(shù)[12]。本文采用2010年6月17日-2010年8月31日昌邑、樓前、滁槎等3個站計算水位與實測水位的對比結(jié)果，對河道糙率進行率定。經(jīng)過調(diào)算，確定贛西聯(lián)圩防洪保護區(qū)贛江一維河道模型糙率取值0.018～0.05之間，平均糙率為0.032。

二維水動力學(xué)模型參數(shù)主要為二維洪泛區(qū)的糙率。由于缺乏實測資料，本文根據(jù)計算區(qū)域地形、地貌的實際情況，結(jié)合以往規(guī)劃設(shè)計資料和經(jīng)驗值分析確定保護區(qū)糙率。計算區(qū)域地表采用不同下墊面不同的糙率分區(qū)，盡可能反映洪水演進的真實情況，如表2所示。

2.5 洪水演進結(jié)果分析

表2 網(wǎng)格糙率取值表Tab.2 Grid roughness values

設(shè)定險工險段永建洲、李家垾發(fā)生瞬時潰堤，潰口寬度分別為322和247 m。對各計算方案的潰口水量平衡進行了分析，具體驗證見表3?？傮w而言，模型對于水量平衡計算的相對誤差均控制在1%以內(nèi)，說明該模型模擬潰堤洪水的精度較高。以方案1為例，從開始模擬時刻2010年6月27日0∶00到6月28日12∶00之間沒有發(fā)生潰決，贛江水位正常升高，保護區(qū)水位等于地面高程，潰口流量為零；在28日12∶00，李家垾潰口發(fā)生潰決，洪水從贛江涌入保護區(qū)，潰口流量由零瞬時達到最大，贛江水位逐漸降低，保護區(qū)內(nèi)水位增高；潰決開始到30日18∶00這一時段內(nèi)，贛江水位高于保護區(qū)內(nèi)水位，洪水從贛江流入保護區(qū)，但水位差逐漸減小，潰口流量逐漸減?。辉?0日18∶00，保護區(qū)內(nèi)水位與贛江水位齊平，此時，潰口流量降為零。潰口流量過程線見圖4。

表3 潰口水量平衡驗證表Tab.3 Water balance verification of breaches

圖4 潰口流量過程線Fig.4 Discharge change process of dyke breaches

根據(jù)潰口洪水分析計算結(jié)果及區(qū)域社會經(jīng)濟情況，并結(jié)合保護區(qū)地物分布，統(tǒng)計各方案的洪水風(fēng)險數(shù)據(jù)，具體結(jié)果見表4?？梢园l(fā)現(xiàn)：在贛江20年一遇鄱陽湖2010年實測水位過程中，上、下游水位較低，潰口流量小，潰口最大流量小于900 m3/s，淹沒面積小，洪水前鋒到達時間小于3 h，淹沒歷時不足24 h，其中永建洲20年一遇方案淹沒面積僅8.46 km2。防洪保護區(qū)內(nèi)的淹沒分布來看，基本遵從地形高低的原則，地勢低洼區(qū)以及洪水主流區(qū)淹沒水深均大于3 m。

表4 洪水風(fēng)險數(shù)據(jù)統(tǒng)計表Tab.4 Statistics of flood risk data

綜合分析贛西聯(lián)圩防洪保護區(qū)洪水風(fēng)險狀況，可以得到以下結(jié)論：①贛西聯(lián)圩防洪保護區(qū)錦江洪水損失大小與潰口位置有關(guān)，從高到低依次為李家垾潰口洪水風(fēng)險、永建洲潰口洪水風(fēng)險，原因在于李家垾潰口位于贛江上游，上游水位較下游更高，由堰流公式可知李家垾潰口流量要大于永建洲潰口流量，同時保護區(qū)水量變化量也更大。②根據(jù)計算成果，突發(fā)潰堤災(zāi)難時，洪水到達時間短，淹沒水深大，洪水主流區(qū)避洪轉(zhuǎn)移的機會較小，應(yīng)分別向兩側(cè)高地、高速公路方向盡快轉(zhuǎn)移，盡可能減少生命財產(chǎn)損失。

最終將計算成果以洪水風(fēng)險圖的形式表達，并附上不同淹沒水深等級時的淹沒數(shù)據(jù)統(tǒng)計，其中方案1部分成果如圖5和表5所示。研究成果對防汛部門制定救災(zāi)決策方案具有重大參考價值。

圖5 洪水風(fēng)險示意圖Fig.5 Flood risk map

3 結(jié) 語

表5 不同水深下的淹沒數(shù)據(jù)統(tǒng)計表Tab.5 Statistics of submerged data in different depths

本文對贛西聯(lián)圩防洪保護區(qū)建立了潰堤洪水一、二維耦合水動力模型，利用糙率分區(qū)、網(wǎng)格加密等技術(shù)優(yōu)化了模型，提高了模型精度。以永建洲、李家垾贛江20年一遇設(shè)計洪水湖口2010年實測水位過程為基礎(chǔ)，對計算區(qū)域進行了模擬，結(jié)果表明，此模型能滿足水量平衡，具有較高的計算精度。獲取洪水風(fēng)險數(shù)據(jù)，結(jié)合區(qū)域社會經(jīng)濟情況做了洪水風(fēng)險分析。與傳統(tǒng)模型相比，該模型具有很好的多功能性和適用性，可應(yīng)用于潰堤及其任意組合方案的洪水演進模擬并具有很好的穩(wěn)定性。目前的潰壩計算主要以瞬間全潰為基礎(chǔ)，其計算簡單方便，所需資料少，成本低，但精度不高，如何考慮潰口發(fā)展過程仍然是目前模型計算中的重要技術(shù)問題。

□

[1] 熊一榮, 張俊暉. 贛西聯(lián)圩存在問題及發(fā)展建議[J]. 現(xiàn)代農(nóng)業(yè)科技,2011,(10):249-251.

[2] 李 麗. 分布式水力模型的匯流演算研究[D]. 南京: 南京大學(xué),2007.

[3] 芮孝芳. Muskingum法及其分段連續(xù)演算的若干理論探討[J]. 水科學(xué)進展,2002,13(6):682-688.

[4] 白玉川, 黃本勝. 河網(wǎng)非恒定流數(shù)值模擬的研究進展[J]. 水利學(xué)報,2000,(12):43-47.

[5] 苑希民, 田福昌, 王麗娜. 漫潰堤洪水聯(lián)算全二維水動力模型及應(yīng)用[J]. 水科學(xué)進展,2015,26(1):83-90.

[6] 宋利祥. 潰壩洪水?dāng)?shù)學(xué)模型及水動力學(xué)特性研究[D]. 武漢:華中科技大學(xué), 2012.

[7] 王 煒. 二維洪水?dāng)?shù)值模擬在太浦河洪水風(fēng)險圖編制中的應(yīng)用[D]. 上海: 同濟大學(xué), 2008.

[8] 張大偉, 李丹勛, 陳稚聰, 等. 潰堤洪水的一維、二維耦合水動力模型及應(yīng)用[J]. 水力發(fā)電學(xué)報, 2010,29(2):149-154.

[9] 李大鳴, 林 毅, 周志華. 蓄滯洪區(qū)洪水演進一、二維數(shù)值仿真及其在洼淀聯(lián)合調(diào)度中的應(yīng)用[J]. 中國工程科學(xué), 2010,12(3):82-88.

[10] 姜曉明, 李丹勛, 王興奎. 基于黎曼近似解的潰堤洪水一維-二維耦合數(shù)學(xué)模型[J]. 水科學(xué)進展, 2012,23(2):214-221.

[11] 付成威, 苑希民, 楊 敏. 實時動態(tài)耦合模型及其在洪水風(fēng)險圖中的應(yīng)用[J].水利水運工程學(xué)報, 2013,(5):32-38.

[12] 張小琴, 包為民, 梁文清, 等. 河道糙率問題研究進展[J]. 水力發(fā)電, 2008,34(6):98-100.