基于非結構網格的分蓄洪區水沙演進數學模型研究

張細兵，歐治華，崔占峰，孫 先

基于非結構網格的分蓄洪區水沙演進數學模型研究

張細兵1，2，歐治華3，崔占峰1，孫 先1

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

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

非結構網格；分蓄洪區；洪水演進；泥沙沖淤；數值模擬

1 概 述

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

為此，筆者在前人研究基礎上，建立了分蓄洪區的水沙演進模型。其特點為，采用非結構三角形網格，能很好擬合復雜的水域邊界，并能自由對重點區域進行加密處理；在模擬洪水演進的同時，還能模擬泥沙輸移和河床沖淤變形，從而為分蓄洪區的洪水調度及防洪搶險提供更為豐富和準確的數據。

2 非結構網格生成技術

復雜流動問題數值計算中采用的網格，根據拓撲結構大致分為結構化網格和非結構化網格2大類。非結構化網格是一種無規則隨機的網格結構，其優點為：復雜區域的網格生成自動化程度高、貼體性好；流動方程直接在物理空間上的非結構化控制體上離散和求解，不需要坐標變換和方程變換，物理空間即為計算空間；能夠較方便地求解復雜邊界的流動且更容易實現網格自適應處理。

非結構化網格生成方法主要有2種：Delaunay三角形化方法和陣面推進法。這2種方法均能實現網格生成的自動化；能夠實現局部加密；能夠加以推廣生成三維非結構化網格。筆者基于陣面推進法，引入了背景網格（Background Mesh）概念，提出了河道有限元網格自動剖分方法［5］。

3 非結構網格上的數值求解方法

3．1 基本方程及其離散

非結構化網格的方法不需要進行坐標轉換，物理空間即為計算空間，因此可采用單一物理坐標系描述流動方程和計算方法。

3．1．1 基本方程

建立在靜水壓強假定下的平面二維笛卡爾坐標系中的水沙運動方程可寫為：

連續性方程

水流運動方程

懸沙運動方程

河床變形方程

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

方程（1）－（4）可寫成以下統一形式，

式中：φ＝｛H，M，N，HS｝；xi＝｛x，y｝；Γφ為廣義擴散系數；Sφ為源項。

3．1．2 方程離散

為保證滿足動量方程守恒，基矢量一般采用總體直角坐標（即對整個計算區域設定直角坐標系）中的速度分量。本文的節點布置方法采用單元中心法，網格為同位網格。對于任意形狀的控制體，方程（6）的控制體積法的積分形式如下，

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

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

將式（7）應用與圖2所示的控制體，可得

這里，j為控制體中界面（邊）的角標。

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

當對流項界面插值采用混合方案，混合因子取［0，1］之間值，并采用延遲修正求解；用中心差分計算節點上的梯度；源項采用線性化處理方案時，可得到如下的離散結果

懸沙運動方程（4）的離散與此完全一致，這里就不單獨寫出其離散表達式。

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

非結構化同位網格SIMPLE算法的計算步驟如下：

（1）給定初始速度場和水位（壓力）場；

（2）由給定的速度場和水位場計算離散動量方程組的系數和源項；

（3）求解動量方程組，得出速度場；

（4）計算水位校正方程中的系數和源項；（5）求解水位校正方程，得到校止水位；（6）修正水位和速度，獲得本層次上滿足連續方程的速度和水位場；

（7）反復迭代，直至收斂。

由以上可見，在非結構化網格上實施SIMPLE算法的步驟與結構化網格上一致，只不過各個步驟的具體實施方式有較大差別。另外，所有公式及說明都是對速度矢量而言的，實際計算時還需要分別對其直角坐標的分量一一進行。

4 模型在荊江分洪區的應用

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

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

采用以上建立的模型，對1954年荊江分洪區分洪過程進行了復演計算。

4．1 模型的建立及相關問題處理

4．1．1 基于DEM模型的地形自動剖分技術

計算網格地形根據美國SRTM 90 m DEM數字模型插值得到。計算區域為整個荊江分蓄洪區，以北閘為進流邊界，南閘為出流邊界，分蓄洪區的圍堤為其不過水邊界。計算網格采用三角網格，節點總數2 964個，單元總數5 614個，最長邊長1 396 m，最短邊長119 m，如圖3所示。

圖3 荊江分洪區網格劃分Fig．3 Grid partition of Jingjiang flood-division area

4．1．2 動邊界模擬

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

4．1．3 參系數取值

由于荊江分蓄洪區內缺乏實測的洪水資料，因此，糙率參考相關資料確定：樹林0．070，旱地0．065，水田0．050，水面0．025。如果某網格內含有多種地形，則按照各種地形糙率的加權平均值確定該網格的糙率。此外，經驗表明，糙率隨水深增加而減小，并趨于穩定，據此規律確定洪水演進計算中網格節點的糙率。

紊動粘性系數采用υt＝αu*h公式計算，α為常數，取為0．5，u*為摩阻流速。

4．2 計算條件

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

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

4．3 洪水演進計算成果分析

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

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

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

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

4．4 泥沙沖淤計算成果分析

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

圖5 各采樣點流速過程線Fig．5 Velocity process curves at each monitoring point

圖6 各采樣點沖淤幅度過程線Fig．6 Erosion and deposition height process curves at each monitoring point

5 結 論

（1）本文建立的分蓄洪區水沙演進數學模型采用了無結構網格，能很好適應分蓄洪區復雜的邊界條件；應用DEM數字模型插值技術，能快速獲取計算域的地形高程；模型可同時模擬干河床的洪水演進和河床的沖淤變形。

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

［1］ 余明輝，張小峰．平面二維潰堤水流泥沙數值模擬［J］．水科學進展，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］ 王嘉松，何友生．淺水流動與污染物擴散的高分辨率計算模型［J］．應用數學和力學，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］ 王志剛，汪德爟，賴錫軍，等．下游為干河床的潰壩水流數值模擬［J］．水利水運工程學報，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］ 張細兵．河道有限元網格自動剖分方法研究［J］．長江科學院院報，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］ 張細兵，范北林．潰壩洪水演進平面二維數學模型初步研究［J］．長江科學院院報，2006，（6）：19－21．（ZHANG Xi-bing，FAN 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

水利部公益性行業專項經費項目（200901003）；公益院所基金項目（CKSF2011001）

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