田福昌， 張興源， 苑希民

(天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室， 天津 300072)

1 研究背景

多年來，我國一直高度重視山洪災(zāi)害防范與預(yù)警工作，山區(qū)潰堤洪水風(fēng)險分析評估是其中一項(xiàng)重要研究內(nèi)容。利用水動力數(shù)值模型模擬計(jì)算潰堤山洪演進(jìn)過程，可為山洪預(yù)警和風(fēng)險調(diào)度提供較為準(zhǔn)確的風(fēng)險信息。山區(qū)溝道潰堤洪水淹沒過程主要包括溝道與保護(hù)區(qū)洪水演進(jìn)以及兩者之間洪水的動態(tài)交互，通過建立溝道與保護(hù)區(qū)洪水聯(lián)合計(jì)算的一、二維水動力耦合模型，采用側(cè)向耦合銜接方式實(shí)現(xiàn)溝道與保護(hù)區(qū)之間水量與動量的實(shí)時動態(tài)傳遞，能夠較為準(zhǔn)確地分析評估潰堤山洪淹沒風(fēng)險。

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

本文建立模擬山洪溝潰堤洪水演進(jìn)過程的一、二維水動力耦合模型，將計(jì)算區(qū)內(nèi)道路、灌渠堤防等概化為寬頂堰并作線性處理，以模擬線狀地物附近水流流態(tài)及其對洪水演進(jìn)的影響。采用非結(jié)構(gòu)化網(wǎng)格剖分計(jì)算區(qū)域，沿線狀地物縮小網(wǎng)格尺度并加密剖分計(jì)算網(wǎng)格，實(shí)現(xiàn)線狀地物等特殊邊界附近網(wǎng)格地形及水流運(yùn)動的動態(tài)耦聯(lián)，在此基礎(chǔ)上利用干濕邊界理論對模型進(jìn)行優(yōu)化處理，建立具有真實(shí)地形的山區(qū)溝道潰堤洪水淹沒計(jì)算模型，并將其應(yīng)用于寧夏大武口溝潰堤洪水淹沒風(fēng)險的分析評估。

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

2.1 一維水動力模型

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

(1)

(2)

式中：B為過流橫斷面寬度，m；qs為溝道旁側(cè)流量，m3/s；QS為溝道總流量，m3/s；A為溝道過水?dāng)嗝婷娣e，m2；g為重力加速度，m/s2；t為時間，s；Z為溝道水位，m；C為謝才系數(shù)；s為距離坐標(biāo)；R為水力半徑；i為溝道底坡降。

利用Abbott六點(diǎn)隱式格式[13]離散上述方程組，該離散格式在每一個網(wǎng)格節(jié)點(diǎn)按順序交替計(jì)算水位和流量，因此能夠在相當(dāng)大的Courant數(shù)下保持計(jì)算穩(wěn)定，可采用較大計(jì)算步長以節(jié)省計(jì)算時間。

2.2 二維水動力模型

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

連續(xù)方程：

(3)

動量方程：

(4)

(5)

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

2.3 一、二維水動力模型耦聯(lián)方式

采用側(cè)向連接方式實(shí)現(xiàn)溝道一維模型與保護(hù)區(qū)二維模型的耦合銜接。潰堤洪水流態(tài)與寬頂堰流較為接近，故潰口流量Qb采用寬頂堰公式[3,6]計(jì)算:

(6)

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

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

3 研究實(shí)例

3.1 研究區(qū)域概況

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

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

圖1 大武口溝模型計(jì)算區(qū)域示意圖

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

3.2.1 地形概化

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

(2)網(wǎng)格剖分。非結(jié)構(gòu)化網(wǎng)格具有很強(qiáng)的邊界適應(yīng)能力，能夠?qū)θ我庑螤詈吐?lián)通區(qū)域進(jìn)行網(wǎng)格剖分，便于控制網(wǎng)格密度，易于修改和調(diào)整，更容易獲得高質(zhì)量網(wǎng)格地形。因此本文采用非結(jié)構(gòu)化三角形網(wǎng)格單元剖分研究區(qū)域，并在興民村揚(yáng)水渠、平汝鐵路、S301、第二農(nóng)場渠沿線以及區(qū)域邊界進(jìn)行網(wǎng)格加密處理，計(jì)算區(qū)域共劃分網(wǎng)格31.20×104個，最小網(wǎng)格面積50 m2，最大網(wǎng)格面積3 500 m2，計(jì)算總面積395.19 km2。

3.2.2 參數(shù)設(shè)定

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

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

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

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

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

3.3 潰堤洪水淹沒計(jì)算結(jié)果與風(fēng)險分析

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

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

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

圖2 石大公路橋斷面100年一遇設(shè)計(jì)洪水過程 圖3 大武口攔洪庫入庫斷面水位-流量關(guān)系圖

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

圖7 洪水演進(jìn)淹沒水深分布圖

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

4 結(jié) 論

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

(2) 將模型應(yīng)用于寧夏大武口溝保護(hù)區(qū)潰堤山洪風(fēng)險分析評估，準(zhǔn)確高效地模擬了潰堤洪水演進(jìn)過程和風(fēng)險分布特征，明確了主要淹沒區(qū)域與防洪重點(diǎn)關(guān)注區(qū)，計(jì)算結(jié)果可為大武口溝防汛指揮與應(yīng)急搶險提供支持與參考。

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