丁志雄，李 娜，許小華，李德龍

（1.中國(guó)水利水電科學(xué)研究院，北京 100038； 2.水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038；3.江西省水利科學(xué)研究院，江西 南昌 330029）

1 研究背景

河流堤防是抵御洪水侵害的一項(xiàng)重要水利工程設(shè)施，堤防一旦發(fā)生因洪水導(dǎo)致的漫頂或潰決，其造成的影響或損害可能極其嚴(yán)重，而尤以堤防潰決造成的災(zāi)害更具有突發(fā)性和不可預(yù)測(cè)性。為了減輕潰堤洪水的災(zāi)害影響，除了轉(zhuǎn)移受災(zāi)群眾和重要財(cái)物外，往往會(huì)采取及時(shí)、快速封堵修復(fù)或排水的防洪搶險(xiǎn)措施。近年來(lái)，國(guó)內(nèi)外許多學(xué)者研究了潰堤洪水的數(shù)值模擬，其主要是利用二維水動(dòng)力模型模擬潰堤洪水在二維平面上的演進(jìn)情況［1-6］。然而，對(duì)于潰堤發(fā)生后可能的采取封堵修復(fù)或排水的搶險(xiǎn)措施考慮的情況少。本文利用適用性較寬泛的二維非恒定流完全圣維南方程組建立潰堤洪水模型，并利用基于共軛單元離散的有限體積法進(jìn)行模型求解。以江西撫河2010年唱?jiǎng)P堤潰堤洪水為例，考慮當(dāng)時(shí)實(shí)際的潰堤封堵修復(fù)過(guò)程和潰堤分流的實(shí)際情況，進(jìn)行潰堤洪水的模擬反演，并與現(xiàn)場(chǎng)采集的實(shí)際洪痕進(jìn)行對(duì)比驗(yàn)證，以期為潰堤洪水的防洪搶險(xiǎn)提供參考和借鑒。

2 區(qū)域概況

唱?jiǎng)P堤位于江西省第二大河撫河的中下游右岸，撫州市臨川區(qū)東北部，涉及湖南、羅湖、唱?jiǎng)P、羅針、云山5個(gè)鄉(xiāng)（鎮(zhèn)），為一條封閉圩堤，圩堤全長(zhǎng)81.8 km，保護(hù)耕地0.93萬(wàn)hm2，保護(hù)人口19.9萬(wàn)人。

2010年6月?lián)岷影l(fā)生大洪水，全線超警戒水位。6月21日18時(shí)30分，唱?jiǎng)P堤在靈山何家段（樁號(hào)32+923～33+270）發(fā)生決口（以下簡(jiǎn)稱“何家段潰口”），起初潰口寬度5 m，后迅速發(fā)展到60 m，到6月22日7時(shí)30分潰口寬度擴(kuò)至347 m，2010年6月22日上午拍攝的潰口處照片如圖1所示［7］，圖中標(biāo)識(shí)了潰口寬和潰口入流方向。堤防潰決造成受災(zāi)鄉(xiāng)鎮(zhèn)4個(gè)、受災(zāi)村41個(gè)，被淹區(qū)平均水深2.5 m至4 m，其中羅針鎮(zhèn)、唱?jiǎng)P鎮(zhèn)受災(zāi)最嚴(yán)重，整個(gè)受淹區(qū)域人口約10萬(wàn)人。

2010年6月23日6時(shí)30分左右，江西撫州唱?jiǎng)P堤內(nèi)的洪水在羅針鎮(zhèn)長(zhǎng)湖村附近再次沖開一個(gè)新缺口（以下簡(jiǎn)稱“長(zhǎng)湖村潰口”），缺口寬度約150 m。洪水一部分泄入撫河，一部分倒灌東鄉(xiāng)河。同時(shí)6月23日武警水電部隊(duì)開始進(jìn)行靈山何家段潰口的封堵，至6月27日18時(shí)靈山何家段潰口勝利合龍［8］。因此6月27日18時(shí)后，唱?jiǎng)P圩堤內(nèi)何家段潰口不再有水量流入。

唱?jiǎng)P堤潰堤洪水具有一定的代表性，首先堤防潰決具有較長(zhǎng)的發(fā)展過(guò)程；其次潰決后有人工堵口修復(fù)工作的介入；還有下游自然沖開堤防進(jìn)行洪水排出等現(xiàn)象的發(fā)生。因此研究清楚潰堤洪水的發(fā)展演進(jìn)過(guò)程，人工堵口修復(fù)對(duì)洪水的影響及其保護(hù)區(qū)進(jìn)水后的排水輸干等問(wèn)題，對(duì)于堤防潰決后的防洪、搶險(xiǎn)、救災(zāi)以及災(zāi)后恢復(fù)重建等都具有非常重要的參考價(jià)值。

圖1 唱?jiǎng)P堤靈山何家決口處實(shí)拍照片（2010年6月22日上午拍攝）

3 潰堤洪水模擬模型及其求解算法

3.1 潰堤洪水模擬模型中國(guó)水利水電科學(xué)研究院于1990年代開發(fā)了一套基于水動(dòng)力學(xué)的二維洪水模擬仿真模型（IWHR-FRAS）［9］，該模型主要采用擴(kuò)散波方程對(duì)洪水進(jìn)行模擬演算，模型較好地適用于城市的暴雨內(nèi)澇及較小區(qū)域的潰堤洪水模擬，近年來(lái)結(jié)合全國(guó)重點(diǎn)地區(qū)洪水風(fēng)險(xiǎn)圖的編制，對(duì)該模型作了進(jìn)一步的改進(jìn)，針對(duì)大范圍的潰堤洪水模擬，模型采用適用性較寬泛的二維非恒定流的完全圣維南方程組，考慮了水流的慣性項(xiàng)和摩阻項(xiàng)，具體如下：

質(zhì)量守恒方程（連續(xù)方程）：

動(dòng)量守恒方程（運(yùn)動(dòng)方程）：

式中：H為水位；h為水深；u、v分別為流速在x、y方向的分量；t為時(shí)刻；q為源匯項(xiàng)；g為重力加速度；ce為渦黏系數(shù)；cf為摩阻系數(shù)，由曼寧糙率確定；fk為科氏力系數(shù)。

3.2 模型求解算法二維非恒定流圣維南方程的數(shù)值求解方法主要有有限差分法（Finite Difference Method，F(xiàn)DM）、有限單元法（Finite Element Method，F(xiàn)EM）以及有限體積法（Finite Volume Method，F(xiàn)VM）等。有限體積法是近年來(lái)興起的一種用于不可壓縮流體計(jì)算的數(shù)值求解方法，由于該方法從物理規(guī)律出發(fā)，每一離散方程都是有限大小體積上的某物理量的守恒表達(dá)式，在推導(dǎo)過(guò)程中物理概念清晰，并可以保證離散方程的守恒性，同時(shí)該方法能像有限單元法一樣適用于不規(guī)則網(wǎng)格和復(fù)雜邊界情況，但處理效率要遠(yuǎn)高于有限單元法，所以在當(dāng)今流體力學(xué)數(shù)值模擬計(jì)算中有著越來(lái)越廣泛的應(yīng)用。這里主要利用基于共軛單元離散的有限體積法進(jìn)行潰堤洪水的完全圣維南方程組的求解，共軛離散單元如圖2所示，共軛單元及節(jié)點(diǎn)不是真實(shí)存在的，由離散網(wǎng)格單元根據(jù)單元的共軛關(guān)系自動(dòng)計(jì)算得到，用于輔助計(jì)算離散網(wǎng)格單元的水位和流量。

圖2 共軛離散單元示例(實(shí)線和圓點(diǎn)為離散單元和節(jié)點(diǎn)，虛線和方點(diǎn)為共軛單元及節(jié)點(diǎn))

將上述完全圣維南方程組的連續(xù)方程寫為矢量形式如下［4］：

式中：V=（u，v）為流速矢量；為偏微分算子。

運(yùn)動(dòng)方程的矢量形式如下：

根據(jù)高斯散度定理，單元的通量與單元水體體積的變化保持恒定［5］。對(duì)于連續(xù)方程（式（4））在空間共軛離散以及時(shí)間離散后每個(gè)單元的水量近似滿足如下關(guān)系式：

式中：Ω（H）為單元不同水位下的水體體積，由單元形態(tài)特征計(jì)算成水位-水體體積關(guān)系查算表；im為離散單元的邊數(shù)；ni第i邊元（單元的一條邊）的法向量；Ai（Ht）為第i邊元的過(guò)水面積，由邊元形態(tài)特征計(jì)算成水位-過(guò)水面積關(guān)系查算表；Q為區(qū)域降雨、蒸發(fā)或入滲量。

對(duì)于運(yùn)動(dòng)方程（式（5）），在空間共軛離散以及時(shí)間離散后每個(gè)節(jié)點(diǎn)流速近似滿足如下關(guān)系式［6］：式中：由廣義重心坐標(biāo)插值得到，包括流速大小和方向；jm為節(jié)點(diǎn)相關(guān)單元個(gè)數(shù)；cj為節(jié)點(diǎn)第j個(gè)相關(guān)單元與所在共軛網(wǎng)格單元形狀關(guān)系系數(shù)；ci為節(jié)點(diǎn)第j個(gè)相關(guān)單元所有節(jié)點(diǎn)的相關(guān)單元與所在共軛網(wǎng)格單元形狀關(guān)系系數(shù)。

在設(shè)定邊界及初始條件下，得到每個(gè)網(wǎng)格的水位及每個(gè)節(jié)點(diǎn)（每個(gè)節(jié)點(diǎn)對(duì)應(yīng)一個(gè)邊元）的流速vt及Ht，從而由式（6）、式（7）可以聯(lián)立求解得到vt+1及Ht+1，并能夠保證在每個(gè)離散網(wǎng)格和時(shí)間步長(zhǎng)上水流的近似質(zhì)量守恒和近似動(dòng)量守恒。

4 唱?jiǎng)P堤潰堤洪水模擬模型建立

唱?jiǎng)P堤洪水分析模型的模擬計(jì)算面積為86.51 km2，概化分為14 130個(gè)不規(guī)則網(wǎng)格，如圖3所示。其中，包括邊元數(shù)28 240條、節(jié)點(diǎn)數(shù)14 758個(gè)。在網(wǎng)格剖分過(guò)程中，以唱?jiǎng)P堤圩區(qū)的計(jì)算范圍作為邊界約束條件，將京福高速公路、316國(guó)道等阻水建筑物作為內(nèi)部約束條件，并以每個(gè)網(wǎng)格面積不超過(guò)0.05 km2為標(biāo)準(zhǔn)進(jìn)行網(wǎng)格剖分，對(duì)于特殊地形地物，如高速公路等進(jìn)行了局部加密處理。網(wǎng)格的糙率按土地利用類型進(jìn)行糙率取值，每個(gè)網(wǎng)格可能包含多種土地利用類型，則用面積加權(quán)平均進(jìn)行求取網(wǎng)格的糙率。

圖3 唱?jiǎng)P堤潰決模擬分析二維模型

圖4 唱?jiǎng)P堤潰口及撫河上下游水文站位置

5 唱?jiǎng)P堤潰堤洪水的模擬反演及結(jié)果分析

5.1 潰堤洪水模擬反演唱?jiǎng)P堤潰堤洪水計(jì)算主要考慮堤防潰決的發(fā)展過(guò)程以及何家段潰口的修復(fù)過(guò)程，潰口處采用水位邊界，模擬時(shí)間步長(zhǎng)為10 s。

（1）潰口發(fā)展過(guò)程。根據(jù)唱?jiǎng)P堤靈山何家段2010年洪水實(shí)際潰決情況，2010年6月21日18時(shí)30分發(fā)生潰決，初始潰口寬度為5 m，最終潰口寬度為347 m，潰口發(fā)展歷時(shí)約13 h。建立潰口的線性發(fā)展方程如下：

式中：Y為潰口寬度；a為潰口初始寬度（5 m）；X為潰口歷時(shí)，min。

隨著唱?jiǎng)P堤地內(nèi)蓄水量的增加6月23日6時(shí)30分左右唱?jiǎng)P堤羅針鎮(zhèn)長(zhǎng)湖村發(fā)生潰決，潰口寬度約150 m，潰口歷時(shí)約5 h，因此建立潰口的線性發(fā)展方程如下：

式中：Y為潰口寬度；X為潰口歷時(shí)，min。

（2）潰口處水位。唱?jiǎng)P圩堤潰口位置位于撫河的廖家灣水文站和李家渡水文站之間（如圖4）。2010年6月唱?jiǎng)P堤發(fā)生潰決洪水時(shí)，這兩個(gè)水文站均有較詳細(xì)的水位記錄，該水位資料全面反映了唱?jiǎng)P堤潰口上下游水位的變化過(guò)程，因此可以用這兩個(gè)站點(diǎn)的實(shí)測(cè)水位資料，按距離的線性關(guān)系插值計(jì)算出唱?jiǎng)P堤潰口處的水位，作為模型的水位邊界輸入。

（3）何家段潰口的修復(fù)過(guò)程。2010年6月23日武警水電部隊(duì)奉命封堵靈山何家段決口，至6月27日18時(shí)靈山何家段決口勝利合龍，因此何家段潰口的修復(fù)從潰口完成后約24 h開始，潰口修復(fù)經(jīng)歷的時(shí)間為107 h，建立何家段潰口開始修復(fù)后的寬度變化方程如下：

式中：Y為開始修復(fù)后的潰口寬度；a為潰口寬度（347 m）；X為潰口修復(fù)歷時(shí)，min。

5.2 反演結(jié)果分析

（1）潰口處入流過(guò)程。模型計(jì)算的唱?jiǎng)P堤何家段潰決入流流量變化過(guò)程如圖5示，由圖可見潰口處入流流量在潰口發(fā)生后約10 h內(nèi)迅速增加，最大達(dá)到1544.475m3/s。最初潰口流量迅速增加的過(guò)程與實(shí)際潰口的發(fā)展過(guò)程基本一致。從何家段潰口處堤內(nèi)外水位變化過(guò)程看（圖6），潰口發(fā)生后，堤外水位一直均高于堤內(nèi)水位，直至6月27日以后堤外水位才降低到比堤內(nèi)水位低，與該潰口處的完成封堵時(shí)間基本一致。

圖5 何家段潰口處入流流量過(guò)程

圖6 何家段潰口處堤內(nèi)外水位變化過(guò)程

（2）出口流量過(guò)程。模型計(jì)算的唱?jiǎng)P堤長(zhǎng)湖村潰決出流流量變化過(guò)程如圖7所示，由圖可見長(zhǎng)湖村潰口處出流流量在6月23日6時(shí)發(fā)生潰決后，出流流量逐漸增加，變化趨勢(shì)與長(zhǎng)湖村潰口處堤內(nèi)外水位的變化過(guò)程基本一致（圖8）。

圖7 長(zhǎng)湖村潰口處出流流量變化過(guò)程

圖8 長(zhǎng)湖村潰口處堤內(nèi)外水位變化過(guò)程

同時(shí)由圖8可以看到，長(zhǎng)湖村潰口處堤內(nèi)（模型網(wǎng)格編號(hào)565）在何家段潰堤后最初8小時(shí)還沒有水，8 h后何家段潰堤洪水到達(dá)，水位逐漸增加，但仍低于長(zhǎng)湖村潰口處堤外水位（東鄉(xiāng)河水位）；到6月22日15時(shí)，長(zhǎng)湖村潰口處堤內(nèi)外水位持平，達(dá)31.31 m；6月23日6時(shí)30分長(zhǎng)湖村潰口，此時(shí)長(zhǎng)湖村堤內(nèi)外水位差0.96 m，唱?jiǎng)P堤圩區(qū)洪水通過(guò)該潰口直接泄入東鄉(xiāng)河直至排入撫河。

隨著長(zhǎng)湖村潰口的發(fā)生，圩區(qū)內(nèi)水位逐漸降低，到6月24日23時(shí)長(zhǎng)湖村潰口處水位基本要接近堤外（東鄉(xiāng)河）水位，但由于后續(xù)洪峰的到來(lái)，何家段潰口入流量的增加，長(zhǎng)湖村潰口處堤內(nèi)外水位差又有所增加，可見模型基本反映了潰口的入流、出流以及圩區(qū)內(nèi)的水位變化情況。

圖9 圩區(qū)內(nèi)進(jìn)出水量及區(qū)內(nèi)剩余水量變化過(guò)程曲線

圖10 洪痕觀測(cè)點(diǎn)模型計(jì)算水深變化過(guò)程

圖11 現(xiàn)場(chǎng)考察GPS記錄路線圖及考察點(diǎn)位

圖12 洪痕水深現(xiàn)場(chǎng)照片

（3）唱?jiǎng)P堤圩區(qū)內(nèi)進(jìn)出水量變化過(guò)程分析。唱?jiǎng)P堤圩區(qū)自何家段潰堤開始至6月27日18時(shí)完成何家段潰堤的封堵，由模型計(jì)算得圩區(qū)內(nèi)的進(jìn)出累積水量及區(qū)內(nèi)蓄水量的變化過(guò)程如圖9所示。由圖9可見唱?jiǎng)P堤圩區(qū)內(nèi)自6月21日16時(shí)何家段潰決開始，總進(jìn)水量持續(xù)增加，到6月27日達(dá)到2.34億m3；圩區(qū)內(nèi)自6月23日6時(shí)長(zhǎng)湖村潰決開始由圩區(qū)排出的水量也持續(xù)增加，到6月27日18時(shí)外排水量達(dá)到2.02億m3；唱?jiǎng)P堤圩區(qū)內(nèi)蓄水量自6月21日16時(shí)何家段潰決開始水量持續(xù)增加，在6月23日6時(shí)以前蓄水量增加較快，6月23日6時(shí)以后受長(zhǎng)湖村潰決自然排水影響，唱?jiǎng)P堤圩區(qū)內(nèi)蓄水量增加逐漸放緩，到6月23日7時(shí)，唱?jiǎng)P堤圩區(qū)內(nèi)蓄水量達(dá)到最大，為1.14億m3，之后圩區(qū)內(nèi)蓄水量開始減小，到6月27日18時(shí)圩區(qū)內(nèi)蓄水量減小到0.32億m3。由于長(zhǎng)湖村潰口處地面高程約為28.35 m，圩區(qū)內(nèi)大部分區(qū)域均高于該高程，因此圩區(qū)內(nèi)的水均可通過(guò)該潰口自然排干。

圖13 潰堤后洪水模擬演進(jìn)不同時(shí)刻淹沒范圍變化

（4）洪痕點(diǎn)水深變化過(guò)程。2015年1月19日到唱?jiǎng)P堤現(xiàn)場(chǎng)進(jìn)行了考察，察看了2010年洪水的何家段潰口位置、長(zhǎng)湖村潰口位置以及仍保留有2010年堤防潰決洪水時(shí)的洪痕位置（對(duì)應(yīng)模型網(wǎng)格編號(hào)6644）。模型計(jì)算洪痕處的水深變化過(guò)程如圖10所示。從圖可見，何家段潰口發(fā)生后約1個(gè)小時(shí)，洪水到達(dá)洪痕觀測(cè)點(diǎn)，并且水深迅速增加到最大水深1.85 m左右。

現(xiàn)場(chǎng)考察行進(jìn)GPS記錄的路線及考察點(diǎn)位如圖11所示。洪痕位置的水深淹沒狀況如現(xiàn)場(chǎng)照片圖12所示，從現(xiàn)場(chǎng)實(shí)測(cè)對(duì)比可見，最高洪痕處距地面約1.95 m?？梢娔Ｐ陀?jì)算結(jié)果與現(xiàn)場(chǎng)考察觀測(cè)到的洪痕水深非常接近，誤差約為0.1 m。

（5）淹沒范圍變化過(guò)程。根據(jù)唱?jiǎng)P堤潰決洪水模擬演進(jìn)計(jì)算結(jié)果，繪制不同時(shí)刻的洪水淹沒范圍及水深分布如圖13。由圖可見潰決發(fā)生后12 h（2010年6月22日4時(shí)），唱?jiǎng)P堤圩區(qū)內(nèi)大部分區(qū)域已被淹，演進(jìn)48 h（2010年6月22日16時(shí)）后唱?jiǎng)P堤圩區(qū)內(nèi)淹沒范圍進(jìn)一步擴(kuò)大，洪水演進(jìn)到第72 h（2010年6月23日16時(shí)）受長(zhǎng)湖村所在圩堤段潰決（2010年6月23日6時(shí)）影響，唱?jiǎng)P堤圩區(qū)內(nèi)淹沒范圍已經(jīng)有所減少，洪水演進(jìn)到第144 h（2010年6月27日18時(shí)），受何家段潰決口門的封堵成功以及長(zhǎng)湖村所在圩堤潰決的自然排水，圩區(qū)內(nèi)淹沒范圍已大大減少。

6 結(jié)論

2010年6月江西撫河唱?jiǎng)P堤潰決洪水是近年來(lái)發(fā)生的一次較大的潰堤洪水災(zāi)害，洪水影響人口十多萬(wàn)人，潰堤有較長(zhǎng)的發(fā)展過(guò)程，潰堤發(fā)生后既有人工修復(fù)的介入，又有堤防自然潰決排水過(guò)程的發(fā)生，因此該堤防潰決具有較突出的典型性。

利用改進(jìn)的中國(guó)水利水電科學(xué)研究院洪澇仿真模型（IWHR-FRAS），針對(duì)潰堤洪水建立了模擬唱?jiǎng)P堤潰堤洪水的數(shù)學(xué)模型，進(jìn)行了潰堤及洪水發(fā)展過(guò)程的模擬反演。結(jié)合潰堤后現(xiàn)場(chǎng)考察的實(shí)測(cè)資料，對(duì)模型反演結(jié)果進(jìn)行了分析，結(jié)果表明：所建立的唱?jiǎng)P堤潰堤洪水模擬模型可以較好地模擬反演潰堤、修復(fù)及排水等因素影響下洪水的發(fā)展變化過(guò)程；計(jì)算的潰口流量過(guò)程與堤內(nèi)外的水位變化過(guò)程相一致；潰堤后唱?jiǎng)P堤保護(hù)區(qū)內(nèi)的水量變化過(guò)程及淹沒范圍變化比較合理，與實(shí)際情況也相吻合；計(jì)算的洪痕處的水深與現(xiàn)場(chǎng)觀測(cè)的洪痕水深比較接近，誤差只有10 cm，符合洪水模擬的一般技術(shù)規(guī)范要求。因此利用該模型，建立潰堤洪水的實(shí)時(shí)模擬分析系統(tǒng)，結(jié)合人工修復(fù)干預(yù)下的實(shí)時(shí)潰堤洪水進(jìn)行模擬演算，則可以為潰堤發(fā)生后的防洪、搶險(xiǎn)、救災(zāi)等工作的開展提供重要的技術(shù)參考。