潰堤山洪淹沒風險評估水動力耦合模型及應用

田福昌， 張興源， 苑希民

(天津大學 水利工程仿真與安全國家重點實驗室， 天津 300072)

1 研究背景

多年來，我國一直高度重視山洪災害防范與預警工作，山區潰堤洪水風險分析評估是其中一項重要研究內容。利用水動力數值模型模擬計算潰堤山洪演進過程，可為山洪預警和風險調度提供較為準確的風險信息。山區溝道潰堤洪水淹沒過程主要包括溝道與保護區洪水演進以及兩者之間洪水的動態交互，通過建立溝道與保護區洪水聯合計算的一、二維水動力耦合模型，采用側向耦合銜接方式實現溝道與保護區之間水量與動量的實時動態傳遞，能夠較為準確地分析評估潰堤山洪淹沒風險。

近年來，眾多學者針對一、二維水動力耦合模型開展了大量研究工作。陳文龍等[1]建立了基于側向聯解的一維-二維耦合水動力模型，并對郁江中下游河段防洪保護區潰堤及漫堤洪水演進進行了模擬；田福昌[2]建立了考慮漫、潰堤位置處流態變化的河道-泛區二維動態耦合水動力數值模型，并對黃河寧夏段洪水漫灘運動以及青銅峽防洪保護區洪水演進過程進行了模擬；付成威等[3]利用建立的一、二維水動力耦合模型模擬了谷堆圩保護區潰堤洪水的演進過程；苑希民等[4]借鑒全二維氣相色譜理論建立了漫潰堤聯算的全二維水動力模型，并成功用于黃河寧蒙段河道及兩岸灌區的漫潰堤洪水模擬計算；陳俊鴻等[5]建立了基于多種優化方法的一、二維耦合水動力模型，模擬了贛西聯圩防洪保護區潰堤洪水演進過程；Dutta等[6]建立了有限差分法的一、二維水動力耦合模型，并對湄公河漫溢洪水運動過程進行了模擬；徐祖信等[7]建立了基于有限元法的一、二維水動力耦合模型，并應用于黃浦江平原感潮河網的流場模擬；Cai Xin等[8]與Xie Zuotao等[9]基于元胞自動機理論建立了洪水演進計算模型并模擬了荊江洞庭湖洪水流動過程；陳智洋等[10]建立了橫陽支江一、二維水動力耦合數值模型，計算了橫陽支江超標準洪水漫堤、潰堤等工況下的洪水演進情況；李長躍[11]、苑希民等[12]建立了考慮潰堤分流與暴雨等多源洪水條件下的一、二維水動力耦合模型，并將其應用于茨南淝左片防洪保護區多源洪水運動的耦合模擬。

本文建立模擬山洪溝潰堤洪水演進過程的一、二維水動力耦合模型，將計算區內道路、灌渠堤防等概化為寬頂堰并作線性處理，以模擬線狀地物附近水流流態及其對洪水演進的影響。采用非結構化網格剖分計算區域，沿線狀地物縮小網格尺度并加密剖分計算網格，實現線狀地物等特殊邊界附近網格地形及水流運動的動態耦聯，在此基礎上利用干濕邊界理論對模型進行優化處理，建立具有真實地形的山區溝道潰堤洪水淹沒計算模型，并將其應用于寧夏大武口溝潰堤洪水淹沒風險的分析評估。

2 一、二維水動力耦合模型基本原理

2.1 一維水動力模型

一維水動力模型采用的基本方程如下：

(1)

(2)

式中：B為過流橫斷面寬度，m；qs為溝道旁側流量，m3/s；QS為溝道總流量，m3/s；A為溝道過水斷面面積，m2；g為重力加速度，m/s2；t為時間，s；Z為溝道水位，m；C為謝才系數；s為距離坐標；R為水力半徑；i為溝道底坡降。

利用Abbott六點隱式格式[13]離散上述方程組，該離散格式在每一個網格節點按順序交替計算水位和流量，因此能夠在相當大的Courant數下保持計算穩定，可采用較大計算步長以節省計算時間。

2.2 二維水動力模型

二維水動力模型采用的基本方程如下[14]：

連續方程：

(3)

動量方程：

(4)

(5)

式中：H為水深，m；Z為水位，m；qc為連續方程中的源匯項；M與N分別為x和y方向的垂向平均單寬流量，m2/s；u和v分別為垂向平均流速在x與y方向的分量，m/s；g為重力加速度，m/s2；n為曼寧糙率系數。

2.3 一、二維水動力模型耦聯方式

采用側向連接方式實現溝道一維模型與保護區二維模型的耦合銜接。潰堤洪水流態與寬頂堰流較為接近，故潰口流量Qb采用寬頂堰公式[3,6]計算:

(6)

h1=max(Z1,Z2)-Zb，h2=max(min(Z1,Z2)-Zb,0)

式中:Z1、Z2分別為一維、二維模型在耦合界面處的水位，m；Zb為耦合界面的底高程，m；lb為潰口寬度，m。

3 研究實例

3.1 研究區域概況

大武口溝位于石嘴山市北側，為賀蘭山東麓最大的山洪溝。大武口溝流域面積576 km2，溝長50 km，平均比降11.50‰，較大支流有北叉溝、塔塔溝、八號泉溝、大燈溝等。溝道范圍內地勢北高南低、西高東低，海拔1 350～1 400 m。溝道右岸為山地陡岸，左岸地勢較平緩，相對高差較小，兩岸植被稀疏，巖石裸露，溝道砂礫等堆積嚴重。

大武口溝為季節性河流，洪水主要由暴雨產生，洪水特性主要表現為年際變化大、季節性特征明顯(多出現在7-8月份)、來勢兇猛、暴發頻繁、洪峰流量隨匯水面積的增大而緩慢增加、洪峰陡漲陡落、峰高量小、峰型尖瘦，洪水歷時一般為6～12 h。本文計算區域為大武口溝城區段(自上游石大公路橋至下游平汝鐵路橋)及其兩岸保護區，總面積為395.19 km2。

圖1 大武口溝模型計算區域示意圖

3.2 潰堤山洪淹沒水動力耦合模型建立

3.2.1 地形概化

(1)溝道斷面設置。溝道斷面是一維水動力模型的重要基礎數據，大武口溝城區段計算范圍為石大公路橋至平汝鐵路橋(大武口攔洪庫入庫斷面)，溝道長度8.02 km，根據溝道實際寬度及其蜿蜒曲折情況，對溝道斷面進行內插加密處理，共設置160個計算斷面，斷面間距均值為500 m。

(2)網格剖分。非結構化網格具有很強的邊界適應能力，能夠對任意形狀和聯通區域進行網格剖分，便于控制網格密度，易于修改和調整，更容易獲得高質量網格地形。因此本文采用非結構化三角形網格單元剖分研究區域，并在興民村揚水渠、平汝鐵路、S301、第二農場渠沿線以及區域邊界進行網格加密處理，計算區域共劃分網格31.20×104個，最小網格面積50 m2，最大網格面積3 500 m2，計算總面積395.19 km2。

3.2.2 參數設定

(1)干濕邊界。模型計算過程中，設定干水深Hdry和濕水深Hwet的分界可以消除不穩定性并提高計算效率，利用干濕水深與網格淹沒水深(H)進行比較，通過網格淹沒水深與干濕水深的關系來判別網格的通量(動量通量和質量通量)計算。若網格H>Hwet，該網格為濕網格，質量通量和動量通量同時參與計算；若網格H

(2)糙率和計算步長。根據《寧夏石嘴山市大武口區大武口溝河流治理工程初步設計報告》、《洪水風險圖編制技術細則(試行)》及《水力計算手冊(第二版)》，確定大武口溝城區段溝道綜合糙率為0.034，堤外平面計算區綜合糙率值為0.06。綜合考慮模型穩定及運算效率等多種因素[15]，設定大武口溝城區段溝道一維水動力模型計算迭代步長為1 s，計算區二維水動力模型最大計算迭代步長為10 s，最小計算迭代步長為0.01 s。為實現一維模型和二維模型固定時間步長內的動態耦合，耦合模型計算時間步長為1 s。

(3)潰口參數。大武口溝城區段堤防設計防洪標準為50年一遇，考慮計算溝道洪水潰堤分流淹沒風險的最大化，按超標準洪水量級選取，故本文確定大武口溝城區段洪水分析標準為100年一遇。考慮大武口溝城區段保護區可能遭遇最大風險，綜合河勢地形、地質狀況、工程狀況和歷史出險情況，確定潰口位置位于平汝鐵路上游1 320 m處，并結合現場調查和防汛專家咨詢意見等確定潰口寬度為60 m，堤防瞬間潰決，潰堤時水位達到設計洪水位。

3.2.3 邊界條件 大武口溝城區段一維水動力模型上邊界條件為石大公路橋斷面遭遇100年一遇洪水，如圖2所示；下邊界為大武口攔洪庫入庫斷面，控制條件為該斷面水位-流量關系，如圖3所示；采用側向銜接概化方式計算潰口分流流量過程，以實現固定時間步長內溝道與保護區洪水的實時動態交互及耦合計算。

考慮到計算區域內道路、灌渠渠堤等線狀地物對洪水演進過程干擾性較強，影響區域流場及淹沒風險分布情況。本文考慮將研究區域內興民村揚水渠、第二農場渠、平汝鐵路和S301等線狀地物作為模型內邊界并概化為寬頂堰，即當洪水水位未達到線狀地物頂部高程時，線狀地物起阻水作用，區域不過水；當洪水水位超過線狀地物頂部高程時，水流漫溢通過。同時考慮線狀地物沿程缺口及橋梁等主要過流形式，將線狀地物在橋梁或缺口處斷開，斷開寬度即為過水區域斷面寬度。

3.3 潰堤洪水淹沒計算結果與風險分析

3.3.1 潰口分流洪水過程分析 根據所建大武口溝城區段潰堤山洪一、二維水動力耦合模型，模擬計算大武口溝100年一遇洪水演進過程以及潰口分流洪水在保護區內的淹沒運動過程，根據洪水計算結果提取平汝鐵路橋上游處潰口(上游臨近斷面)溝道流量過程(如圖4所示)和水位過程(如圖5所示)，潰口側向分流流量過程如圖6所示。

平汝鐵路橋上游堤防潰決時潰口處斷面洪水水位等于設計水位，同時考慮大武口溝堤防現狀防洪標準較高及抗洪搶險能力等多種因素，設定潰堤分流截止時潰口處溝道水位等于地面高程(潰口分流積水后地面高程)。由圖4和圖5分析可知，潰口處溝道流量為845 m3/s時堤防開始潰決，此時溝道洪水水位為1 108.03 m，潰堤分流截止時溝道洪水水位為1 107.49 m，潰口分流歷時為7 h。由圖6分析可知，起潰時刻潰口分流流量為35.74 m3/s，分流洪峰流量為197.10 m3/s，經統計分流總量為179.98×104m3。

3.3.2 潰堤洪水演進及風險分析 洪水由平汝鐵路橋上游潰口處側向分流進入二維平面計算區域，由所建大武口溝城區段潰堤山洪一、二維水動力耦合模型計算得到不同時刻洪水淹沒水深分布情況，如圖7所示。

圖2 石大公路橋斷面100年一遇設計洪水過程 圖3 大武口攔洪庫入庫斷面水位-流量關系圖

圖4 潰口上游臨近斷面溝道流量過程 圖5 洪水平汝鐵路橋上游處溝道水位變化過程 圖6 平汝鐵路橋上游潰口側向分流流量過程

圖7 洪水演進淹沒水深分布圖

由圖7分析可知：受區域地形分布影響，潰堤洪水進入計算區后向東南方向演進，受到平汝鐵路和興民村揚水渠阻擋；洪水演進2 h，興民村揚水渠與平汝鐵路之間區域逐漸被淹沒，洪水向第二農場渠方向演進，計算區內最大水深3.02 m，淹沒面積1.12 km2；洪水演進4 h，洪鋒北側到達S301與平汝鐵路之間中心區域，計算區內最大水深3.21 m，淹沒面積2.49 km2；洪水演進8 h，受S301阻擋，中心區域大面積積水，計算區內最大水深2.71 m，淹沒面積6.03 km2；洪水演進16 h時已基本趨于穩定狀態，受興民村揚水渠、平汝鐵路和S301阻水作用，淹沒區最大積水面積7.18 km2，最大積水深度為3.32 m(位于平汝鐵路路基附近)。考慮洪水沖刷作用影響，大武口溝沿程橋梁、護坡護岸工程、水文測站、計算區內渠堤及道路路基等將遭受嚴重損害，應提前做好汛前檢查與應急搶險防御工作。

4 結 論

(1) 建立一、二維水動力耦合模型模擬山區溝道潰堤洪水演進過程。為了準確反映計算區線狀地物對水流運動的影響，采用非結構化網格剖分與局部網格加密處理技術對計算區域進行了地形概化，利用干濕邊界理論優化模型，并將淹沒區具有擋水作用的線狀構筑物概化為寬頂堰，實現了淹沒區特殊邊界與非結構化網格的無縫耦聯，以更好地擬合復雜地形，建立具有真實地形的山區溝道潰口分流洪水淹沒計算模型。

(2) 將模型應用于寧夏大武口溝保護區潰堤山洪風險分析評估，準確高效地模擬了潰堤洪水演進過程和風險分布特征，明確了主要淹沒區域與防洪重點關注區，計算結果可為大武口溝防汛指揮與應急搶險提供支持與參考。

(3) 所建大武口溝潰堤山洪淹沒風險評估水動力耦合模型，可推廣應用于溝道工程規劃設計與防洪評價、洪災評估與風險區劃等領域。由于未有歷史大場次山洪潰堤事件的詳細記載資料，本文尚無法實現模型更為精準的率定與驗證，后續可根據歷史及現場資料搜集情況對其進行補充完善，以最大限度提高模型計算精準度，為防汛指揮與應急決策等提供更為準確的風險信息。