基于一、二維耦合水動力模型的贛西聯圩潰堤洪水風險分析

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

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

贛西聯圩防洪保護區位于贛江主支下游西岸，北臨鄱陽湖，保護區內水系繁多，同時水文條件錯綜復雜。加上贛西聯圩設計標準偏低，保護區內地形平坦，地勢低洼，汛期同時受贛江和錦江洪水的頂托影響，圩內洪澇災害嚴重。近60年贛西聯圩曾多次發生潰決，給圩區人民生產、生活造成嚴重的經濟損失[1]。因此進行贛西聯圩洪水風險分析，對于防汛指揮部門做出風險決策，保護當地居民生命財產安全，具有十分現實意義。

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

1 模型構建原理

1.1 河道一維模型

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

(2)

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

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

(4)

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

(5)

式中：θ為權重系數，可取0.5<θ<1.0。

1.2 洪泛區二維模型

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

(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為曼寧系數。

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

(16)

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

1.3 模型耦合

將一維河網模型與二維洪水模型進行耦合，構建最終的洪水分析模型。選定側堰流公式來實現潰口上下游水流信息的交互，具體為：

(17)

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

(18)

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

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

2 贛西聯圩洪水風險分析

2.1 區域概況

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

2.2 地形概化

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

圖1 贛西聯圩防洪保護區示意圖Fig.1 Schematic diagram of Ganxi dyke flood control protected area

對贛西聯圩防洪保護區計算區域采用非結構不規則網格進行網格劃分，使網格的大小隨地形地勢和阻水建筑物的分布靈活確定，而且盡可能地將影響水流的阻水建筑物作為網格邊界，充分反映計算域的特征。但是，必要的時候對保護區內的一些典型的線性阻水建筑物，如堤防、公路等，經合理概化，并對網格適當加密，在二維地形中充分反映其特征。研究區域部分斷面概化和網格概化如圖2所示。

圖2 斷面網格概化圖Fig.2 Generalized cross-section and topography

2.3 潰口設置與邊界條件

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

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

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

2.4 糙率確定

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

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

2.5 洪水演進結果分析

表2 網格糙率取值表Tab.2 Grid roughness values

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

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

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

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

表4 洪水風險數據統計表Tab.4 Statistics of flood risk data

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

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

圖5 洪水風險示意圖Fig.5 Flood risk map

3 結 語

表5 不同水深下的淹沒數據統計表Tab.5 Statistics of submerged data in different depths

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

□

[1] 熊一榮, 張俊暉. 贛西聯圩存在問題及發展建議[J]. 現代農業科技,2011,(10):249-251.

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

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

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

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

[6] 宋利祥. 潰壩洪水數學模型及水動力學特性研究[D]. 武漢:華中科技大學, 2012.

[7] 王 煒. 二維洪水數值模擬在太浦河洪水風險圖編制中的應用[D]. 上海: 同濟大學, 2008.

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

[9] 李大鳴, 林 毅, 周志華. 蓄滯洪區洪水演進一、二維數值仿真及其在洼淀聯合調度中的應用[J]. 中國工程科學, 2010,12(3):82-88.

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

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

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