基于非結(jié)構(gòu)網(wǎng)格的分蓄洪區(qū)水沙演進(jìn)數(shù)學(xué)模型研究

張細(xì)兵，歐治華，崔占峰，孫 先

基于非結(jié)構(gòu)網(wǎng)格的分蓄洪區(qū)水沙演進(jìn)數(shù)學(xué)模型研究

張細(xì)兵1，2，歐治華3，崔占峰1，孫 先1

（1．長(zhǎng)江科學(xué)院河流研究所，武漢 430010；2．武漢大學(xué)水利水電學(xué)院，武漢 430072；3．荊州市長(zhǎng)江河道管理局洪湖分局，湖北洪湖 433200）

基于非結(jié)構(gòu)三角形網(wǎng)格，采用有限體積算法，建立了分蓄洪區(qū)水沙演進(jìn)數(shù)學(xué)模型。模型網(wǎng)格能很好適應(yīng)分蓄洪區(qū)復(fù)雜的邊界條件，應(yīng)用了DEM數(shù)字模型插值技術(shù)，能快速獲取計(jì)算域的地形數(shù)據(jù)，實(shí)現(xiàn)了洪水演進(jìn)過(guò)程中泥沙輸移和河床沖淤變形的同步模擬。模型中還對(duì)網(wǎng)格劃分、數(shù)值離散、動(dòng)邊界及干河床處理等關(guān)鍵技術(shù)進(jìn)行了研究，并提供了相應(yīng)的處理方法。該模型在荊江分洪區(qū)得到成功運(yùn)用，較好地復(fù)演了1954年分洪過(guò)程。洪水計(jì)算成果可為分蓄洪區(qū)防洪決策與搶險(xiǎn)提供參考，分洪閘后最大沖刷深度計(jì)算成果可為工程設(shè)計(jì)提供參考。

非結(jié)構(gòu)網(wǎng)格；分蓄洪區(qū)；洪水演進(jìn)；泥沙沖淤；數(shù)值模擬

1 概 述

河流在自然或人工潰堤情況下，常常會(huì)有水沙的急變運(yùn)動(dòng)發(fā)生，洪水在洪泛區(qū)內(nèi)的演進(jìn)將給洪泛區(qū)內(nèi)人民生命財(cái)產(chǎn)安全構(gòu)成巨大威脅。因此，正確模擬潰壩洪水在洪泛區(qū)內(nèi)的演進(jìn)過(guò)程，將給防洪搶險(xiǎn)和防洪決策提供重要的依據(jù)。由于問(wèn)題的復(fù)雜性，國(guó)內(nèi)外許多學(xué)者開展了分蓄洪區(qū)洪水演進(jìn)的研究。如武漢大學(xué)的余明輝、張小峰［1］采用矩形網(wǎng)格和有限體積算法，建立了分蓄洪區(qū)洪水演進(jìn)數(shù)學(xué)模型；王嘉松［2］、王志剛［3］基于矩形網(wǎng)格，分別采用有限差分算法和有限體積算法，對(duì)局部潰壩水流演進(jìn)問(wèn)題進(jìn)行了研究；Macchione和Morelli［4］提出了用于求解潰壩水流的MacComack格式。國(guó)外許多商業(yè)化的軟件也具備了分蓄洪區(qū)洪水演進(jìn)功能，如丹麥DHI水環(huán)境研究所研制的MIKE-Flood模型能進(jìn)行分蓄洪區(qū)洪水演進(jìn)的一、二維嵌套模擬；美國(guó)密西西比大學(xué)研制的CCHE-Flood能進(jìn)行潰壩洪水演進(jìn)的模擬等。但以上模型大都采用結(jié)構(gòu)化網(wǎng)格，一般未考慮河床沖淤變形問(wèn)題。

為此，筆者在前人研究基礎(chǔ)上，建立了分蓄洪區(qū)的水沙演進(jìn)模型。其特點(diǎn)為，采用非結(jié)構(gòu)三角形網(wǎng)格，能很好擬合復(fù)雜的水域邊界，并能自由對(duì)重點(diǎn)區(qū)域進(jìn)行加密處理；在模擬洪水演進(jìn)的同時(shí)，還能模擬泥沙輸移和河床沖淤變形，從而為分蓄洪區(qū)的洪水調(diào)度及防洪搶險(xiǎn)提供更為豐富和準(zhǔn)確的數(shù)據(jù)。

2 非結(jié)構(gòu)網(wǎng)格生成技術(shù)

復(fù)雜流動(dòng)問(wèn)題數(shù)值計(jì)算中采用的網(wǎng)格，根據(jù)拓?fù)浣Y(jié)構(gòu)大致分為結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格2大類。非結(jié)構(gòu)化網(wǎng)格是一種無(wú)規(guī)則隨機(jī)的網(wǎng)格結(jié)構(gòu)，其優(yōu)點(diǎn)為：復(fù)雜區(qū)域的網(wǎng)格生成自動(dòng)化程度高、貼體性好；流動(dòng)方程直接在物理空間上的非結(jié)構(gòu)化控制體上離散和求解，不需要坐標(biāo)變換和方程變換，物理空間即為計(jì)算空間；能夠較方便地求解復(fù)雜邊界的流動(dòng)且更容易實(shí)現(xiàn)網(wǎng)格自適應(yīng)處理。

非結(jié)構(gòu)化網(wǎng)格生成方法主要有2種：Delaunay三角形化方法和陣面推進(jìn)法。這2種方法均能實(shí)現(xiàn)網(wǎng)格生成的自動(dòng)化；能夠?qū)崿F(xiàn)局部加密；能夠加以推廣生成三維非結(jié)構(gòu)化網(wǎng)格。筆者基于陣面推進(jìn)法，引入了背景網(wǎng)格（Background Mesh）概念，提出了河道有限元網(wǎng)格自動(dòng)剖分方法［5］。

3 非結(jié)構(gòu)網(wǎng)格上的數(shù)值求解方法

3．1 基本方程及其離散

非結(jié)構(gòu)化網(wǎng)格的方法不需要進(jìn)行坐標(biāo)轉(zhuǎn)換，物理空間即為計(jì)算空間，因此可采用單一物理坐標(biāo)系描述流動(dòng)方程和計(jì)算方法。

3．1．1 基本方程

建立在靜水壓強(qiáng)假定下的平面二維笛卡爾坐標(biāo)系中的水沙運(yùn)動(dòng)方程可寫為：

連續(xù)性方程

水流運(yùn)動(dòng)方程

懸沙運(yùn)動(dòng)方程

河床變形方程

式中：H為水深（m）；u和v分別為x和y方向的流速（m／s），M＝uh，N＝vh；Z為水位（m）；u0，v0，S0分別為源項(xiàng)流速分量和含沙量；n為Manning糙率系數(shù)；vt為紊動(dòng)粘性系數(shù)；q為單位面積上水流的源匯強(qiáng)度（m／s）；ρ為水體密度（kg／m3）；S為含沙量；S*為挾沙力；ω為泥沙顆粒沉速（m／s）；γ′為泥沙干密度（kg／m3）；α為恢復(fù)飽和系數(shù)；ε為泥沙擴(kuò)散系數(shù)。

方程（1）－（4）可寫成以下統(tǒng)一形式，

式中：φ＝｛H，M，N，HS｝；xi＝｛x，y｝；Γφ為廣義擴(kuò)散系數(shù)；Sφ為源項(xiàng)。

3．1．2 方程離散

為保證滿足動(dòng)量方程守恒，基矢量一般采用總體直角坐標(biāo)（即對(duì)整個(gè)計(jì)算區(qū)域設(shè)定直角坐標(biāo)系）中的速度分量。本文的節(jié)點(diǎn)布置方法采用單元中心法，網(wǎng)格為同位網(wǎng)格。對(duì)于任意形狀的控制體，方程（6）的控制體積法的積分形式如下，

式中：V為控制體的體積（單元面積）；A為控制體的表面積（單元邊長(zhǎng)度）；A為控制容積界面的面積矢量（邊矢量）；其正向與外法線單位矢量一致（圖1）；u為速度矢量。

圖1 單元矢量布置示意Fig．1 Sketch of unit vector

將式（7）應(yīng)用與圖2所示的控制體，可得

這里，j為控制體中界面（邊）的角標(biāo)。

圖2 界面上平均值的確定Fig．2 Determ ination of the average value on the interface

當(dāng)對(duì)流項(xiàng)界面插值采用混合方案，混合因子取［0，1］之間值，并采用延遲修正求解；用中心差分計(jì)算節(jié)點(diǎn)上的梯度；源項(xiàng)采用線性化處理方案時(shí)，可得到如下的離散結(jié)果

懸沙運(yùn)動(dòng)方程（4）的離散與此完全一致，這里就不單獨(dú)寫出其離散表達(dá)式。

3．2 離散方程的求解步驟

非結(jié)構(gòu)化同位網(wǎng)格SIMPLE算法的計(jì)算步驟如下：

（1）給定初始速度場(chǎng)和水位（壓力）場(chǎng)；

（2）由給定的速度場(chǎng)和水位場(chǎng)計(jì)算離散動(dòng)量方程組的系數(shù)和源項(xiàng)；

（3）求解動(dòng)量方程組，得出速度場(chǎng)；

（4）計(jì)算水位校正方程中的系數(shù)和源項(xiàng)；（5）求解水位校正方程，得到校止水位；（6）修正水位和速度，獲得本層次上滿足連續(xù)方程的速度和水位場(chǎng)；

（7）反復(fù)迭代，直至收斂。

由以上可見，在非結(jié)構(gòu)化網(wǎng)格上實(shí)施SIMPLE算法的步驟與結(jié)構(gòu)化網(wǎng)格上一致，只不過(guò)各個(gè)步驟的具體實(shí)施方式有較大差別。另外，所有公式及說(shuō)明都是對(duì)速度矢量而言的，實(shí)際計(jì)算時(shí)還需要分別對(duì)其直角坐標(biāo)的分量一一進(jìn)行。

4 模型在荊江分洪區(qū)的應(yīng)用

荊江分洪區(qū)位于荊江南岸，是荊江地區(qū)防洪系統(tǒng)的主要組成部分，分洪區(qū)面積921 km2，1995年全區(qū)共有人口48．25萬(wàn)，有效蓄洪容量54億m3。荊江分洪工程包括進(jìn)洪閘、節(jié)制閘和分洪區(qū)圍堤。進(jìn)洪閘位于太平口東岸，設(shè)計(jì)進(jìn)洪流量8 000 m3／s。

荊江河段的安全泄量約為68 000 m3／s（包括向洞庭湖分流），其中沙市河段只能通過(guò)50 000 m3／s。超過(guò)這個(gè)標(biāo)準(zhǔn)，則需要考慮運(yùn)用荊江分洪區(qū)。荊江分蓄洪區(qū)的主要功能在于分泄荊江河段的超承洪水流量，降低荊江的洪水位，使洪水能夠安全通過(guò)，而不致超過(guò)荊江大堤的防御能力，以確保荊江兩岸廣大農(nóng)田和人民生命財(cái)產(chǎn)安全。同時(shí)由于分洪后荊江水位降低，四口分泄入湖的流量減少，并且由于荊江下泄量減少，岳陽(yáng)出湖流量也可能相對(duì)增加，因而也就間接地減輕了洞庭湖的負(fù)擔(dān)，在1954年的洪水情況下，這項(xiàng)工程曾發(fā)揮了“江湖兩利”的重大作用。

采用以上建立的模型，對(duì)1954年荊江分洪區(qū)分洪過(guò)程進(jìn)行了復(fù)演計(jì)算。

4．1 模型的建立及相關(guān)問(wèn)題處理

4．1．1 基于DEM模型的地形自動(dòng)剖分技術(shù)

計(jì)算網(wǎng)格地形根據(jù)美國(guó)SRTM 90 m DEM數(shù)字模型插值得到。計(jì)算區(qū)域?yàn)檎麄€(gè)荊江分蓄洪區(qū)，以北閘為進(jìn)流邊界，南閘為出流邊界，分蓄洪區(qū)的圍堤為其不過(guò)水邊界。計(jì)算網(wǎng)格采用三角網(wǎng)格，節(jié)點(diǎn)總數(shù)2 964個(gè)，單元總數(shù)5 614個(gè)，最長(zhǎng)邊長(zhǎng)1 396 m，最短邊長(zhǎng)119 m，如圖3所示。

圖3 荊江分洪區(qū)網(wǎng)格劃分Fig．3 Grid partition of Jingjiang flood-division area

4．1．2 動(dòng)邊界模擬

潰壩水流計(jì)算關(guān)鍵在于如何進(jìn)行干河床與動(dòng)邊界的處理。筆者在前人研究基礎(chǔ)上，提出了能適應(yīng)潰壩水流動(dòng)邊界的處理方法［6］。其主要思路是：首先給定一很小水深hmin，根據(jù)計(jì)算得到的節(jié)點(diǎn)水深來(lái)判斷節(jié)點(diǎn)的“干濕”情況，當(dāng)h＞hmin時(shí)為濕節(jié)點(diǎn)；否則為干節(jié)點(diǎn)，并令干節(jié)點(diǎn)h＝hmin。對(duì)于干節(jié)點(diǎn)，還需根據(jù)相鄰節(jié)點(diǎn)水位來(lái)判斷是否進(jìn)行“凍結(jié)”處理，若其周圍均為干節(jié)點(diǎn)，或其周圍有濕節(jié)點(diǎn)但其水位低于該干節(jié)點(diǎn)河床高程時(shí)，則對(duì)該干節(jié)點(diǎn)進(jìn)行凍結(jié)處理，即糙率取為極大值（如1010）；否則“不凍結(jié)”，節(jié)點(diǎn)糙率為正常值。

4．1．3 參系數(shù)取值

由于荊江分蓄洪區(qū)內(nèi)缺乏實(shí)測(cè)的洪水資料，因此，糙率參考相關(guān)資料確定：樹林0．070，旱地0．065，水田0．050，水面0．025。如果某網(wǎng)格內(nèi)含有多種地形，則按照各種地形糙率的加權(quán)平均值確定該網(wǎng)格的糙率。此外，經(jīng)驗(yàn)表明，糙率隨水深增加而減小，并趨于穩(wěn)定，據(jù)此規(guī)律確定洪水演進(jìn)計(jì)算中網(wǎng)格節(jié)點(diǎn)的糙率。

紊動(dòng)粘性系數(shù)采用υt＝αu*h公式計(jì)算，α為常數(shù)，取為0．5，u*為摩阻流速。

4．2 計(jì)算條件

1954年洪水具有流量大、歷時(shí)長(zhǎng)的特點(diǎn)，荊江工程曾先后3次開閘運(yùn)用：第1次分洪從7月22日2時(shí)20分至27日13時(shí)10分，分洪總量達(dá)23．5億m3；第2次分洪從7月29日6時(shí)15分至8月1日15時(shí)55分，分洪量約17．2億m3；第3次分洪從8月1日21時(shí)40分至22日7時(shí)50分，分洪量約81．9億m3（包括臘林洲扒口分洪在內(nèi)）。以上3次分洪總量為122．6億m3（8月22日7時(shí)50分以后開閘及臘林洲扒口流量未計(jì)在內(nèi)）。

模擬計(jì)算時(shí)假定閘門瞬間開啟，忽略閘門的開啟過(guò)程。閘門上游水位保持不變直至蓄滿為止。

4．3 洪水演進(jìn)計(jì)算成果分析

計(jì)算結(jié)果表明，分洪閘下泄流量先緩慢增大，后快速減小，最大下泄流量可達(dá)8 000 m3／s，分洪區(qū)蓄滿所需時(shí)間為11 d左右，分洪量約為60億m3。

圖4給出了荊江分蓄洪區(qū)分洪后3 h和24 h時(shí)流速矢量和淹沒(méi)水深套繪圖，由圖中可見在縮窄部位流速增大，水流首先占據(jù)低洼部位，然后水位逐步抬升。

圖4 不同時(shí)刻流速與水深套繪圖Fig．4 Flow velocity and water depth at different times

圖5 為各采樣點(diǎn)流速過(guò)程線圖。由圖可見，在整個(gè)分洪過(guò)程中，潰口處的P1采樣點(diǎn)流速呈逐漸減少趨勢(shì)，其它采樣點(diǎn)流速在分洪初期增大，后隨著水位的抬高流速逐漸減少，最大流速為1．2 m／s。同一時(shí)刻，離潰口越近的采樣點(diǎn)流速一般越大。

4．4 泥沙沖淤計(jì)算成果分析

圖6為各采樣點(diǎn)沖淤幅度過(guò)程線（圖中負(fù)值表示沖刷）。由于P1采樣點(diǎn)位于分洪口門附近，故此處河床處于強(qiáng)烈沖刷狀態(tài)，沖刷初期含沙量最大，為2．0 kg／m3，沖刷深度在分洪后3h內(nèi)增加最快，最大沖深為1．8 m；P2采樣點(diǎn)位于分蓄洪區(qū)內(nèi)束窄段的進(jìn)口處，流速較大，除了初期略有淤積外，其它時(shí)間河床一直處于沖刷狀態(tài)，其最大沖刷深度為1．2 m；采樣點(diǎn)P3和P4由于距離分洪口門較遠(yuǎn)，大部分泥沙在其前面河床已經(jīng)落淤，因此這些位置河床淤積幅度很小。

圖5 各采樣點(diǎn)流速過(guò)程線Fig．5 Velocity process curves at each monitoring point

圖6 各采樣點(diǎn)沖淤幅度過(guò)程線Fig．6 Erosion and deposition height process curves at each monitoring point

5 結(jié) 論

（1）本文建立的分蓄洪區(qū)水沙演進(jìn)數(shù)學(xué)模型采用了無(wú)結(jié)構(gòu)網(wǎng)格，能很好適應(yīng)分蓄洪區(qū)復(fù)雜的邊界條件；應(yīng)用DEM數(shù)字模型插值技術(shù)，能快速獲取計(jì)算域的地形高程；模型可同時(shí)模擬干河床的洪水演進(jìn)和河床的沖淤變形。

（2）模型在荊江分洪區(qū)得到初步應(yīng)用，復(fù)演了1954年分洪過(guò)程，得到主要成果為：最大分洪流量約為8 000 m3／s，蓄滿時(shí)分洪量約為60億m3，所需時(shí)間為11 d左右；分洪閘后最大流速、含沙量和沖深分別為1．2 m／s、2．0 kg／m3和1．8 m。計(jì)算結(jié)果表明模型能適應(yīng)荊江分洪區(qū)復(fù)雜的地形邊界條件，模擬效果較好，成果可為分蓄洪區(qū)洪水調(diào)度及防洪搶險(xiǎn)提供參考。

［1］ 余明輝，張小峰．平面二維潰堤水流泥沙數(shù)值模擬［J］．水科學(xué)進(jìn)展，2001，12（3）：286－290．（YU Ming-hui，ZHANG Xiao-feng．Horizontal 2-D Uneven Sedi-mentMathematicalModel in Dike Burst［J］．Advances in Water Science，2001，12（3）：286－290．（in Chinese））

［2］ 王嘉松，何友生．淺水流動(dòng)與污染物擴(kuò)散的高分辨率計(jì)算模型［J］．應(yīng)用數(shù)學(xué)和力學(xué)，2002，23（7）：661－666．（WANG Jia-song，HE You-sheng．High-Resolution Numerical Model for Shallow Water Flows and Pollutant Diffusions［J］．Applied Mathematics and Mechanics，2002，23（7）：661－666．（in Chinese））

［3］ 王志剛，汪德爟，賴錫軍，等．下游為干河床的潰壩水流數(shù)值模擬［J］．水利水運(yùn)工程學(xué)報(bào)，2003，（2）：18－23．（WANG Zhi-gang，WANG De-guan，LAIXi-jun，et al．Numerical Simulation of Dam-Break Current Flowing in Dry Bed［J］．Hydro-science and Engineering，2003，（2）：18－23．（in Chinese））

［4］ Macchione F，Morelli M A．Practical Aspects in Compa- ring Shock-Capturing Schemes for Dam-break Problems［J］．J．Hydraul．Eng．，2002，129：187－195．

［5］ 張細(xì)兵．河道有限元網(wǎng)格自動(dòng)剖分方法研究［J］．長(zhǎng)江科學(xué)院院報(bào)，2002，（3）：19－21．（ZHANG Xi-bing．Au-to-Generating Method of Triangulation Grids for River Reach［J］．Journal of Yangtze River Scientific Research Institute，2002，（3）：19－21．（in Chinese））

［6］ 張細(xì)兵，范北林．潰壩洪水演進(jìn)平面二維數(shù)學(xué)模型初步研究［J］．長(zhǎng)江科學(xué)院院報(bào)，2006，（6）：19－21．（ZHANG Xi-bing，F(xiàn)AN Bei-lin．Numerical Simulation for Two-Dimensional Dam-Break Flood Flow［J］．Journal of Yangtze River Scientific Research Institute，2006，（6）：19－21．（in Chinese） ）

（編輯：周曉雁）

Unstructured Grid M odel of Flow and Sediment Evolution in Flood Diversion and Storage Area

ZHANG Xi-bing1，2，OU Zhi-hua3，CUIZhan-feng1，SUN Xian1
（1．River Research Institute，Yangtze River Scientific Research Institute，Wuhan 430010 China；2．College ofWater Resources and Hydroelectric Engineering，Wuhan University，Wuhan 430072 China；3．Honghu Sub-bureau，Jingzhou River Management Bureau of Yangtze River，Honghu 433200 China）

Based on unstructured triangular grid and Finite Volume Method，a flow and sedimentevolutionmodel in flood diversion and storage area is founded in this paper．One advantage is that themodel grids can adapt to compli-cated boundary conditions well，and the terrain data of the computational domain are obtained quickly by using DEM numericalmodel interpolation technology．Another advantage is that the sediment transport and riverbed de-formation can be simulated synchronously．Key techniques including grid partition，numerical separation，treat-ments ofmovable boundary and dry riverbed are studied and the corresponding processingmethods are introduced．Furthermore，themodel has been applied successfully to the simulation of flood diversion in 1954 in Jingjiang flood diversion and storage area．The calculated flood results can serve as reference for flood control policy-making and flood rescue，and themaximum erosion depth behind the flood diversion gate can provide reference for engineering design．

unstructured grid；flood diversion and storage area；flood evolution；sediment erosion and deposition；nu- merical simulation

TV143

A

1001－5485（2011）04－0075－05

2011-02-17

水利部公益性行業(yè)專項(xiàng)經(jīng)費(fèi)項(xiàng)目（200901003）；公益院所基金項(xiàng)目（CKSF2011001）

張細(xì)兵（1976-），湖北鄂州人，高級(jí)工程師，博士研究生，主要從事河流數(shù)值模擬研究，（電話）027-82829871（電子信箱）jss9871＠vip．163．com。