, , , , ,

(珠江水利委員會 a.珠江水利科學(xué)研究院； b.水利部珠江河口動力學(xué)及伴生過程調(diào)控重點實驗室； c.防汛抗旱辦公室， 廣州 510611)

一、二維耦合數(shù)學(xué)模型在感潮河網(wǎng)洪水風(fēng)險圖編制中的應(yīng)用

楊莉玲a, b,宋利祥a,鄧軍濤c,徐爽c,胡曉張a,孫倩雯a

(珠江水利委員會 a.珠江水利科學(xué)研究院； b.水利部珠江河口動力學(xué)及伴生過程調(diào)控重點實驗室； c.防汛抗旱辦公室， 廣州 510611)

為了解決感潮河網(wǎng)徑流、潮汐交匯，動力復(fù)雜，潰決洪水難以用經(jīng)驗公式準(zhǔn)確概化的問題，建立直接以潰口為耦合斷面的河網(wǎng)一維、保護(hù)區(qū)二維側(cè)向耦合模型，將感潮河網(wǎng)與保護(hù)區(qū)一體化，避免環(huán)境因素及經(jīng)驗參數(shù)的不確定性帶來的潰口流量估算誤差。典型算例和中順大圍潰決洪水情景模擬表明：洪水自潰口集中噴射出后分散流向圍內(nèi)，流態(tài)受圍內(nèi)下墊面影響顯著，潰口水位、流量隨外江潮位漲落而起伏變化，由潰口流量過程線計算所得潰口水量與根據(jù)淹沒區(qū)各單元的面積和水深計算的圍內(nèi)總水量一致。模擬成果直接反映下墊面、水頭差、潮位漲落，及潰口流態(tài)對潰決洪水的綜合影響，反映出本耦合模型計算潰口流量及對潰決洪水模擬的合理性、可信性，具有較好的應(yīng)用前景。

潰決洪水；一、二維耦合數(shù)學(xué)模型；感潮河網(wǎng)；洪水風(fēng)險圖；側(cè)向聯(lián)解

1 研究背景

為推進(jìn)我國洪水風(fēng)險圖的編制工作，水利部2008年啟動了“全國洪水風(fēng)險圖編制項目(一期)”工作，2012年啟動(二期)項目，全面開展洪水風(fēng)險圖的編制與推廣應(yīng)用工作。

中順大圍地處中山市西北部、順德區(qū)南部，東臨東海水道及小欖水道，西臨西江干流及磨刀門水道，圍內(nèi)河網(wǎng)交錯、水系復(fù)雜，人口密度大，土地肥沃，經(jīng)濟(jì)發(fā)達(dá)，外江洪潮頂托后潰決風(fēng)險較大，在“全國洪水風(fēng)險圖編制試點項目(二期)——風(fēng)暴潮風(fēng)險分析技術(shù)開發(fā)與試點項目”研究中被確定為試點區(qū)域[1]。

感潮河網(wǎng)區(qū)的潰決洪水同時受徑流、潮汐影響，十分復(fù)雜，潰口流量難以簡單用經(jīng)驗公式概化，在進(jìn)行洪水潰決情景預(yù)測模擬時，更是難以準(zhǔn)確確定潰口水頭，應(yīng)用有一定難度[2-3]。本文基于一、二維側(cè)向耦合模型，直接以潰口斷面作為耦合聯(lián)解斷面，實時反映外江潮漲、潮落及圍內(nèi)下墊面對潰口流量的影響，以及潰口分流對外江潮位的影響作用。外江一維模型用于模擬上游徑流及口外潮波的傳遞、遭遇、頂托作用，模擬范圍大、速度快；二維區(qū)域采用三角形單元剖分，可精細(xì)概化各潰口及圍內(nèi)結(jié)構(gòu)物、小河涌等。

本文以典型算例及中順大圍為例對建立的一、二維側(cè)向耦合潰決洪水?dāng)?shù)值模型進(jìn)行了驗證和應(yīng)用模擬,結(jié)果表明本耦合模型計算潰口流量及對潰決洪水的模擬合理、可信，可推廣應(yīng)用。

2 洪水演進(jìn)數(shù)學(xué)模型簡介

2.1 一維洪水演進(jìn)數(shù)學(xué)模型

一維非恒定流可采用圣維南方程組表示[4-5]：

(1)

(2)

式中：x為空間坐標(biāo)；t為時間；Z為水位；B為河寬；Q為流量；q為源、匯項；A為過流面積；g為重力加速度；u為流速；β為斷面動量修正系數(shù)；R為水力半徑；c=R1/6/n,為謝才系數(shù),n為曼寧綜合糙率；ul為旁側(cè)入流流速在主流方向的分量。

河網(wǎng)中各河道之間通過汊點連接實現(xiàn)耦合，即：

Zi,1=Zi,2=…=Zi,k；

(3)

(4)

式中:下標(biāo)i，k表示第i個汊點的第k條河道;kin，kout分別代表流入、流出汊點的河道數(shù)。

采用Preissmann四點隱式差分格式進(jìn)行方程離散，運(yùn)用四級聯(lián)解法進(jìn)行模型求解[6]。

2.2 二維洪水演進(jìn)數(shù)學(xué)模型

2.2.1 基本方程

二維淺水非恒定流連續(xù)方程和動量方程表示為[7-8]

(5)

式中：U為守恒變量；E，F(xiàn)分別為x,y方向通量；S為源項；S0x，S0y分別為x，y方向的底坡；Sfx，Sfy分別為x，y方向的摩阻坡降；h為水深；u，v分別為x，y方向的垂線平均流速。

2.2.2 數(shù)值求解

采用守恒性較好的有限體積法進(jìn)行控制方程離散，采用三角形網(wǎng)格進(jìn)行二維模型區(qū)域單元剖分，變量定義于單元中心。采用結(jié)構(gòu)簡單且穩(wěn)定性較好的HLLC近似Riemann算子計算界面通量，數(shù)值求解的詳細(xì)過程可參閱文獻(xiàn)[9]。

圖1 一、二維模型側(cè)向耦合示意圖Fig.1 Sketch of side coupling between 1-D and 2-D models

2.3 一、 二維模型側(cè)向耦合

一維河網(wǎng)模型與二維模型的側(cè)向耦合通過耦合斷面(潰口)的水力因子傳遞實現(xiàn)[7]，具體過程為：耦合斷面上、下游一維斷面計算水位按距離插值后提供給耦合斷面網(wǎng)格結(jié)點，采用寬頂堰流公式計算耦合斷面過流量Qk(潰口流量)，Qk既為二維模型邊界，也為耦合斷面下游側(cè)一維斷面的源、匯項(即式(2)中q)，二維模型計算后返回耦合斷面結(jié)點水位給一維模型進(jìn)行下一步計算。一、二維模型耦合方式示意圖如圖1。

2.4 一、二維耦合模型驗證

圖2 計算區(qū)域平面布置Fig.2 Layout of the simulation area

采用被廣泛用于一、二維耦合模型驗證的典型算例[7]來驗證本耦合模型，計算區(qū)域如圖2所示。河道寬4 m、底高程為0 m；蓄滯洪區(qū)地面高程為1.80 m，潰口寬8 m、堤頂高程為2.0 m。入口為流量邊界，取4.0 m3/s，出口為水位邊界，取2.2 m；一維河道單元空間步長取40 m，潰口局部加密，糙率取0.02，二維模型采用三角網(wǎng)格剖分，糙率取0.05，起始水深為0 m。

經(jīng)模擬所得t=0, 10, 20 min的河道水面線計算結(jié)果如圖3(a)所示，文獻(xiàn)[7]的模擬結(jié)果如圖3(b)所示，二者所反映水面線變化情況基本一致：堤防潰決后，河道水面迅速下降，20 min后水面調(diào)整至穩(wěn)定狀態(tài)。對比t=0及t=20 min的水面線可得：潰口上游河段，相同邊界入流下，兩水面線比降亦相同；潰口附近，潰口下游相鄰斷面的流速小于上游相鄰斷面，為了維持總能量平衡，部分動能轉(zhuǎn)化為勢能，水位沿水流運(yùn)動方向呈上漲趨勢；潰口下游河段，在潰口分洪作用下水面變平緩?？梢姡疚哪Ｐ统晒呛侠砜煽康?。

圖3 河道水面線計算結(jié)果對比Fig.3 Comparison of the computed water surface line

本模型直接基于圣維南控制方程進(jìn)行求解計算，不但精細(xì)刻畫了一、二維模型區(qū)域水動力變化，而且避免了經(jīng)驗法公式系數(shù)在不同河道的適應(yīng)性問題。

3 模型應(yīng)用

本文以“0814”黑格比風(fēng)暴潮遭遇上游5 a一遇洪水條件下的中順大圍潰決洪水為例進(jìn)行模型應(yīng)用分析。

3.1 模型研究范圍及網(wǎng)格剖分

一維模型范圍包含珠江三角洲及上游干流，上邊界分別取至馬口、三水、老鴉崗(流溪河)、麒麟咀(增江)、博羅(東江)、石咀(潭江)站；下邊界取至八大口門水文控制站的大虎、南沙、馮馬廟、橫門、燈籠山、黃金、西炮臺、官沖，共布置斷面3 621個(間距200～700 m)，模擬河道總長約1 710 km。

二維區(qū)域采用三角形網(wǎng)格單元，共剖分13萬余個網(wǎng)格，東、西邊界分別取至小欖水道和磨刀門水道，北邊界取至均安水道，南邊界取至五桂山山脈附近，總面積約800 km2，含15個鎮(zhèn)區(qū)，20余條圍內(nèi)河涌。

一、二維耦合模型研究范圍及耦合斷面位置如圖4所示。

圖4 一、二維耦合模型研究范圍及耦合斷面位置Fig.4 Study area of the coupled 1-D and 2-D modelsand locations of the coupling sections

圖5 潰口高程變化過程Fig.5 Development process of breach elevation

3.2 計算條件及邊界

3.2.1 初始條件

一維模型給定恒定流場(Q=Q0,Z=Z0)，二維模型給定靜流場，陸域初始水位、流速均為0；河涌內(nèi)初始水位取常水位。

3.2.2 邊界條件

一維模型上邊界采用流量控制，下邊界采用潮位過程控制，具體由一、二維風(fēng)暴潮流聯(lián)解模型中提取[1]；二維模型非耦合斷面均采用固壁邊界條件。

3.2.3 潰口條件

采用矩形斷面概化。潰口寬度參考韓國《洪水風(fēng)險繪制導(dǎo)則》推薦的無潰痕經(jīng)驗公式，即

Bb=1.9(lgB)4.8+20 。

(6)

式中:Bb為潰口寬(m);B為河寬(m)。 本潰口附近水道寬約為282 m，推算潰口寬度約為160 m。

潰口發(fā)展過程中高程變化如圖5所示，在(t1-t0)時間內(nèi)，潰口由h2潰至h1處，各參數(shù)可參考河堤地質(zhì)進(jìn)行設(shè)置，本研究中h1取圍內(nèi)近潰口處高程，h2取河堤頂高，t1與t0間隔5 s，瞬時潰決。

4 計算成果及分析

4.1 潰口流量

堤防潰口處的流量直接關(guān)系到潰堤洪水的演進(jìn)過程[10]，圖6為潰口處的水位過程線以及潰口處流量數(shù)學(xué)模型與經(jīng)驗公式(見參考文獻(xiàn)[11])計算結(jié)果對比。由圖6可知，兩者峰值相當(dāng)(誤差3.40%)、過程線相似。在外江潮位升高、潰口潰決后，潰口流量突增，然后隨著外江潮位的繼續(xù)抬高，潰口內(nèi)外潮差增至最大，潰口流量也隨之達(dá)峰值;其后，潰口潮位跟隨外江潮位一起漲落，潰口流量也隨之起伏變化，至潰口潮位降至潰口高程以下后，潰口處不再分洪進(jìn)流，潰口流量基本降為0。此外，在潰口低水位時，圍內(nèi)洪水呈向潰口區(qū)倒流，對潰口入流形成頂托之勢，從而模型計算潰口流量略小于經(jīng)驗公式計算值。

圖6 潰口處流量數(shù)學(xué)模型與經(jīng)驗公式計算結(jié)果對比Fig.6 Comparison of embankment-break dischargebetween numerical model and empirical formula

4.2 水 量

潰決洪水水量統(tǒng)計方式有2種：一是潰口泄入水量，由潰口流量過程線計算所得；二是圍內(nèi)總水量，根據(jù)淹沒區(qū)各單元的面積和水深計算、累加而得。本研究中2種方式統(tǒng)計逐時水量見圖7所示，可見，2種方法計算的水量逐時過程基本一致，變化趨勢完全相同，但由于潰口區(qū)域小，潰口流量不斷減小過程中的逐時統(tǒng)計削平了洪峰，造成一定水量損失，從而使?jié)⒖谒柯孕∮趪鷥?nèi)淹沒區(qū)統(tǒng)計總水量。

圖7 潰口總進(jìn)水量與淹沒區(qū)總水量對比Fig.7 Comparison of total water amount between breach discharge and flooded area

4.3 潰口沿程水面線、流態(tài)、流速

由于潰口的落差大，流速急，計算容易出現(xiàn)不穩(wěn)定現(xiàn)象，使計算結(jié)果失真；因此，分析潰口附近水面線、流態(tài)圖的合理性是把握穩(wěn)定性的重要環(huán)節(jié)。

為了分析潰口周邊區(qū)域水面線合理性，由潰口向西至海洲水道、向南至珠三角環(huán)線高速分別布設(shè)2條水位采樣線，繪制潰決后1，6，20，40，60 h的沿程水面線，如圖8(a)、圖8(b)所示。

圖8 潰口西向和南向沿程各點水面線Fig.8 Water surface lines along the west and south directions of breach

從圖8中可以看出，總體上，各時刻西向及南向水面線均表現(xiàn)為自潰口向外水面逐漸下降的變化趨勢，在距潰口4 000 m左右，受地勢大幅下降的影響，洪水傳播至此后水面也呈下降趨勢，由圖8(a)中6 h水面線可明顯看出。當(dāng)洪水進(jìn)一步潰入，地勢低洼區(qū)域水位不斷上升，水面坡降減緩，從細(xì)節(jié)處合理地反映了潰決洪水淹沒時的水面變化情況。

圖9 潰口潰決后其附近區(qū)域流態(tài)(5 h)及流速采樣點布置Fig.9 Flow pattern in the adjacent of breach five hours after embankment-break and location of sampling points

圖9給出了潰口四周區(qū)域5 h的流態(tài)圖。

從圖9中可看出：潰決洪水以潰口為中心向西、向南及西北方向泄流；潰口及地勢低洼區(qū)水流湍急，流勢強(qiáng)勁，而地勢較高區(qū)域，則水流分散，流勢較弱。此外，洪水自潰口泄入圍內(nèi)后，迅速流向西南流板涌內(nèi)，導(dǎo)致流板涌北部水面抬高，涌內(nèi)地勢較低，流速也較大于兩岸區(qū)域，部分洪水經(jīng)流板涌外溢至河涌南岸區(qū)域。流場圖較合理地再現(xiàn)了潰口潰決后的洪水傳播情況。

為了進(jìn)一步分析潰口附近區(qū)域流速計算成果的合理性，在潰口局部區(qū)域布置了6個流速采樣點(位置見圖9)，提取各采樣點流速、水位過程線如圖10所示。從圖10中可以看出：潰決后，由潰口向圍內(nèi)水位迅速抬升，流速急劇增大，以潰口處流速為最大；進(jìn)入圍內(nèi)后隨著泄流斷面迅速擴(kuò)寬，水流流速不斷減小。5#點位于潰口西北，地勢略高于四周(高程為2.71 m)，當(dāng)外江潮位下降至2.71 m以下后，該區(qū)域露灘，水位不再和其他測點一樣隨外江潮位的變化而變化，流速也隨之降為0。

圖10 各采樣點流速和水位的過程線Fig.10 Process lines of flow velocity and water level at sampling points

4.4 洪水淹沒情況

中順大圍內(nèi)潰決洪水最大淹沒水深見圖11所示。由圖11可見：風(fēng)暴潮洪水潰堤后，洪水由潰口進(jìn)入圍區(qū)，順小欖鎮(zhèn)較低地勢區(qū)域淹沒潰口周邊區(qū)域，部分洪水經(jīng)小欖鎮(zhèn)匯入流板涌，并迅速灌入中順大圍內(nèi)其他相連河涌水系，在低洼區(qū)域又溢出河涌，潰決洪水最終漫延至東升鎮(zhèn)之西北區(qū)域，以潰口附近的低洼區(qū)淹沒水深為最大，最大淹沒水深接近1.5 m。

圖11 潰決洪水最大淹沒水深分布Fig．11 Distributionofmaximuminundateddepthofembankment?breakflood

5 結(jié) 論

本文在潰口內(nèi)外分別采用一維模型和三角網(wǎng)格剖分的二維模擬技術(shù)，建立以潰口為耦合斷面的一、二維耦合洪水演進(jìn)模型，能實時反映徑流、潮汐綜合影響下的潰口水量傳遞和變化。通過典型算例及中順大圍潰決情景模擬計算表明，模型計算潰口流量、流速、流態(tài)均較合理，能較好地模擬出潰決洪水在圍內(nèi)不同下墊面特征區(qū)域的演進(jìn)情況。

[1] 珠江水利委員會珠江水利科學(xué)研究院. 風(fēng)暴潮風(fēng)險分析技術(shù)開發(fā)與試點編制報告[R]. 廣州：珠江水利委員會珠江水利科學(xué)研究院，2013.

[2] 夏軍強(qiáng)，王光謙，LIN Bin-liang，等. 復(fù)雜邊界及實際地形上潰壩洪水流動過程模擬[J].水科學(xué)進(jìn)展，2010 ，21(3)：289-298.

[3] 董耀華，黃煜齡. 天然河道長河段一維非恒定流數(shù)模研究[J]. 長江科學(xué)院院報,1994,11(2)：10-17.

[4] 李大鳴，胡永文，王 笑，等. 城市暴雨瀝澇數(shù)學(xué)模型的研究與應(yīng)用[J]. 長江科學(xué)院院報，2016，33(10)： 46-50.

[5] 韋直林，崔占峰. 一維河網(wǎng)非恒定流計算程序的初步研究[J]. 人民長江，2001, 32(12)：30-32.

[6] 吳壽紅. 河網(wǎng)非恒定流四級解法[J]. 水利學(xué)報，1985，(8)：42-50.

[7] 汪定揚(yáng)，明道華，周建華. 三峽工程多級船閘中間引渠水力學(xué)問題的數(shù)值計算研究[J]. 長江科學(xué)院院報， 1989，6(4)：6-14.[8] 張細(xì)兵，金 琨，林木松.潮流河段溫排水影響的平面二維數(shù)值模擬[J].長江科學(xué)院院報，2006，23(3)：13-16.[9] 王智勇，陳永燦，朱德軍，等. 一維-二維耦合的河湖系統(tǒng)整體水動力模型[J]. 水科學(xué)進(jìn)展，2011，22(4)：516-522.

[10]姚志堅，陳文龍. 潰壩洪水演進(jìn)數(shù)學(xué)模型研究及其在大鏡山水庫的應(yīng)用[J].人民珠江，2008 ，(3)：7-9.

[11]李 煒.水利計算手冊[K].2版.北京：中國水利水電出版社，2006.

(編輯：姜小蘭)

Application of a 1-D and 2-D Coupled Mathematical Modelin the Mapping of Tidal River Flood Risk

YANG Li-ling1,2, SONG Li-xiang1,2, DENG Jun-tao3, XU Shuang3, HU Xiao-zhang1,SUN Qian-wen1
(1.Pearl River Hydraulic Research Institute,Pearl River Water Resources Commission, Guangzhou 510611, China; 2.Key Laboratory of the Pearl River Estuarine Dynamics and Associated Process Regulation，Ministry of Water Resources，Guangzhou 510611，China; 3.The Pearl River Basin Flood Control and Drought Relief Office, Pearl River Water Resources Commission, Guangzhou 510611, China)

As the dynamics of tidal river network is complex with runoff and tide confluence, breach discharge in tidal river network is too complex to be estimated by empirical formula accurately. In view of this, a hydrodynamic model coupling 1-D model of river network and 2-D model of reserve area is developed with the breach as the lateral linkage. The model integrates the tidal river network and the reserve area, hence avoiding the deviation of breach discharge caused by the uncertainty of environmental factors and empirical parameters. Simulation results of typical example and embankment-break flood scene show that 1) the jet-flow of breach disperses into the reserve area, and the flow pattern is affected by the underlying surface significantly; 2) water level and discharge of the breach change with tide fluctuation, and the calculated total water amount according to breach discharge agrees with that from flooded areadepth. The simulation results reflect the combined influence of underlying surface, water head, tide level and flow pattern on breach flood. The proposed model is rational and reliable in calculating breach discharge and simulating embankment-break flood.

embankment-break flood; 1-D and 2-D coupling model; tidal river network; mapping of flood risk; lateral linkage

2016-06-12；

：2016-09-19

國家自然科學(xué)基金青年基金項目(51409286)

楊莉玲(1976-)，女，湖北宜昌人，高級工程師，博士，主要從事水動力、水環(huán)境數(shù)值模擬相關(guān)研究，(電話)15913122415(電子信箱)yanglilingok@163.com。

10.11988/ckyyb.20160590

2017,34(9):36-40

TV122.4

：A

：1001-5485(2017)09-0036-05