田福昌， 張興源， 苑希民

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

1 研究背景

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

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

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

2 一、二維水動(dòng)力耦合模型基本原理

2.1 一維水動(dòng)力模型

一維水動(dòng)力模型采用的基本方程如下：

(1)

(2)

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

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

2.2 二維水動(dòng)力模型

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

連續(xù)方程：

(3)

動(dòng)量方程：

(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 一、二維水動(dòng)力模型耦聯(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，溝長(zhǎng)50 km，平均比降11.50‰，較大支流有北叉溝、塔塔溝、八號(hào)泉溝、大燈溝等。溝道范圍內(nèi)地勢(shì)北高南低、西高東低，海拔1 350～1 400 m。溝道右岸為山地陡岸，左岸地勢(shì)較平緩，相對(duì)高差較小，兩岸植被稀疏，巖石裸露，溝道砂礫等堆積嚴(yán)重。

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

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

3.2 潰堤山洪淹沒(méi)水動(dòng)力耦合模型建立

3.2.1 地形概化

(1)溝道斷面設(shè)置。溝道斷面是一維水動(dòng)力模型的重要基礎(chǔ)數(shù)據(jù)，大武口溝城區(qū)段計(jì)算范圍為石大公路橋至平汝鐵路橋(大武口攔洪庫(kù)入庫(kù)斷面)，溝道長(zhǎng)度8.02 km，根據(jù)溝道實(shí)際寬度及其蜿蜒曲折情況，對(duì)溝道斷面進(jìn)行內(nèi)插加密處理，共設(shè)置160個(gè)計(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)場(chǎng)渠沿線以及區(qū)域邊界進(jìn)行網(wǎng)格加密處理，計(jì)算區(qū)域共劃分網(wǎng)格31.20×104個(gè)，最小網(wǎng)格面積50 m2，最大網(wǎng)格面積3 500 m2，計(jì)算總面積395.19 km2。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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