城市地區(qū)水庫潰壩洪水演進模擬

，，

(華南理工大學 亞熱帶建筑科學國家重點實驗室，廣州 510641)

1 研究背景

城市地區(qū)是生命和財產高度集中的區(qū)域，區(qū)內水庫對城市的供水、防洪、興利、生態(tài)等發(fā)揮著重要的作用。隨著我國城市的快速擴張與發(fā)展，一些原在郊外的水庫已經進入了城區(qū)，這類城市水庫一旦潰決，大量水體瞬間下泄將給下游高度集中的生命和財產帶來無法估量的損失。因此，對城市地區(qū)水庫潰壩進行研究，探討潰壩洪水演進特性，對制定應急預案和防洪減災具有重要意義。

城市水庫與山區(qū)水庫在洪水演進研究方面的主要區(qū)別在于下游地區(qū)的地形、地表類型的不同。城市水庫下游分布著各種建筑物，這些建筑物會對演進的洪水產生阻滯、分流、延緩等作用。與此同時，諸如混凝土路面、瀝青路面、草地、叢林帶等組成的城市特有多元化的地表類型更給洪水演進帶來許多不確定性。因此，相對下游地形和地表較為簡單的山區(qū)水庫，城市水庫潰壩洪水的演進則顯得十分復雜。國內外不少學者對潰壩洪水進行過研究，數值模擬是這些研究中的主流研究方法[1-3]，而數值模擬的關鍵步驟是對地形進行處理。對于復雜的城市地形，一些研究對建筑物采取高度概化處理[4-5]，模擬結果雖在大尺度范圍內得到一定保證，但難以反映潰壩洪水在市區(qū)小范圍的發(fā)展特性，地形的失真使得難以獲取更精確的水流發(fā)展信息。在城市地區(qū)，不同類型的地表有著截然不同的糙率，是一個糙率“場”，不能簡單地以1，2個糙率值來代表。針對這種情況，本文在深圳龍口水庫的潰壩模擬研究課題中，利用GIS技術對下游復雜城市地形進行處理并建立高仿真的數字地形模型，對下游不同地表類型賦予不同糙率值以建立多元糙率場，以MIKE21水動力模塊為平臺構建了二維潰壩洪水演進模型，采取分段模擬方法[6]模擬下游研究區(qū)域洪水演進的過程。根據模擬計算結果，分析了下游的流態(tài)變化、淹沒面積、特征點水深變化等洪水信息，為制定潰壩事故早期預警以及人員疏散方案提供了科學的參考依據。

2 MIKE21水動力模型

2.1 模型簡介

MIKE21是丹麥水力研究所(DHI)開發(fā)的系列水動力學軟件之一，屬于平面二維自由表面流模型，其在平面二維自由表面數值模擬方面具有強大的功能。水動力模型是MIKE21最核心的基礎模塊，能模擬由于各種作用力的作用而產生的水位及水流變化。當用戶為模型提供了地形、底部糙率、風場和水動力學邊界條件等數據后，模型能計算出每個網格的水位和水流變化。

2.2 控制方程

MIKE21水動力模型采取有限體積法進行模擬，其主要控制方程為：

(1)

(2)

(3)

式中：t為時間；x,y為笛卡爾坐標系坐標；η為水位；d為靜止水深；h=η+d為總水深；u,v分別為x,y方向上的速度分量；Pa為大氣壓；f為哥氏力系數，f=2ωsinφ,ω為地球自轉角速度，φ為當地緯度；g為重力加速度；ρ為水的密度;ρ0為參考水密度；sxx,sxy，xyy分別為輻射應力分量；S為源項；us,vs分別為源項x,y方向上水流流速分量；τs,τb分別為風應力和水底摩擦力；Tij為水平黏滯應力項，包括黏性力、紊流應力和水平對流，這些量均可根據沿水深平均的速度梯度用渦流黏性方程得出，即

(4)

(5)

3 城市地區(qū)潰壩洪水演進模型

3.1 模型構建步驟

當城市水庫的土石壩由于各種因素被破壞造成潰決，大量水體將下泄到下游城市地區(qū)。由于受建筑物、道路、河道等眾多復雜要素的影響，水流在城市地區(qū)演進時會產生諸多不確定性，其演進特性往往比地形簡單的山區(qū)洪水的演進更為復雜。 為了解下泄洪水在城市地區(qū)的演進特性和發(fā)展態(tài)勢，計算洪水淹沒范圍和淹沒水深，本文以MIKE21的水動力模塊為平臺，構建城市地區(qū)水庫潰壩洪水演進模型，模型構建步驟如下：

(1) 利用GIS對水庫下游復雜地形進行處理，生成高仿真度地形模型。

(2) 對下游不同地表類型賦予不同糙率值，利用GIS生成糙率場。

(3) 根據水庫庫容、水深和潰口寬度，通過經驗公式推求潰口流量過程。

(4) 以水動力模型為平臺，確定研究范圍后設置初始、邊界條件，劃分網格后輸入地形、糙率數據并插值，輸入流量過程，設定其他參數后進行模擬。

(5) 對模擬結果進行水量平衡、流態(tài)、淹沒范圍、水深、流速等洪水特性進行分析。

由于潰壩洪水洪峰量較大、歷時較短，城市排洪管網及水流下滲的水量相對較小，因此在本次研究中暫不考慮城市排洪管網及下滲的影響。

3.2 地形數字高程模型

城市地區(qū)地形的顯著特征是聳立著形式各異、高矮不同、規(guī)則或不規(guī)則的各種建筑物，這些建筑物結合著公路、街道、山地、河道等元素使得城市地區(qū)地形變得極為復雜。構建合理的地形模型是洪水潰決模擬的最基礎也是最重要一步，模型合理與否會直接決定著下游洪水的發(fā)展態(tài)勢和演進結果。本研究借助GIS技術對下游城市地形進行處理并生成地形數字高程模型(DEM)，具體方法如下：

(1)將CAD地形圖中帶有坐標及高程信息的高程點導入GIS中，利用克里金插值法生成數字高程模型，并用柵格格式顯示出來。

(2)對CAD格式中閉合建筑物輪廓線(沒有閉合的需先進行閉合)導入GIS中，使用線矢量生成面矢量，對這些面矢量的建筑物賦予高度值，將面矢量轉化為柵格格式并疊加到上述的數字高程模型(DEM)即生成城市地區(qū)的數字高程模型。

(3)對河道的處理與建筑物處理方法相似，在上述城市地區(qū)數字高程模型減去河岸到河底的高度差即可。

3.3 糙率場的構建

糙率是洪水演進的一個重要參數，其設置合理與否直接影響模擬結果。當水庫潰決時，洪水演進所淹沒的城市地區(qū)則看作河床，對于擁有復雜地面類型的河床不能簡單設置一個糙率值，需要根據流經河床地面類型特征確定不同的糙率值。本研究參考相關試驗[7]和文獻[8]，結合實際情況對不同地表類型賦予的不同的糙率值，利用GIS技術生成城市地區(qū)的糙率場，具體方法如下：

(1) 在CAD地形圖上，根據具體建筑群、草地、叢林帶、公園、公路、街道、河流等特征對研究區(qū)域劃分出不同形狀的板塊，將這些板塊導入GIS中生成面矢量。

(2)對生成面矢量賦予相應糙率值，并將面矢量轉化為柵格格式即可構建區(qū)域的糙率場。

3.4 潰口流量計算

土石壩潰口處瞬時流量的計算有多種方法，不同方法計算結果也略有不同。本文潰決方式采取瞬間潰決，根據文獻[9]，采用如下系列公式計算瞬潰的最大潰壩流量：

(6)

式中：Qmax為大壩潰決最大流量(m3/s)；g為重力加速度(m/s2)；B為壩長(m)；bm為最終潰口寬度(m)；H0為潰決水深(m)。

洪水流量過程采用下式：

(7)

式中：W為水庫庫容(m3)；t為洪水排水歷時(s)。

4 實例分析

4.1 龍口水庫概況

龍口水庫位于深圳市東部龍崗區(qū)龍城街道，是東深供水水源工程的調節(jié)水庫，始建于1994年并于1995年8月竣工，主要解決龍崗區(qū)中心城及龍崗、龍城街道的供水，兼顧防洪任務。水庫共有5座副壩、1座主壩，其中主壩最大壩高28 m，壩頂長510 m，壩頂寬為7.8 m，上游坡比1∶3.25，下游坡比1∶2.75～1∶3，排水棱體高度8 m，筑壩材料為土石材料，土質為粉質黏土，采取逐層均勻碾壓修筑。本次研究模擬主壩潰決后水流從水庫下泄到下游城區(qū)演進過程。

4.2 模型輸入數據處理與生成

潰壩洪水演進模型最重要的數據為地形、糙率以及流量過程數據。根據龍口水庫下游城市地形的特點，利用上述地形處理方法生成龍口水庫數字高程模型(DEM)，如圖1(a)所示。插值生成的3D地形模型(3DEM)如圖1(b)所示，此時疊加了建筑物以及減去河道的數字高程模型可以達到很高的仿真度，確保為模型提供準確的地形數據。

(a) 數字高程模型

(b) 3D高程模型

糙率場的生成：建筑區(qū)的糙率與樓房的疏密、形狀、相對位置等眾多因素有關，本研究對建筑群糙率值取0.2；城市中街道不同于平整的混凝土面，具有明顯的阻礙物，例如街燈、綠化、路牌等因素影響，因此瀝青和混凝土路面均取0.025；龍崗河道經過整治后河床較為平整，河岸堤防已均用較為光滑的面磚修砌，河道糙率統(tǒng)一采取0.03；對于下游的草地取0.04，叢林帶取0.08，其它面積較少的地表類型給予相應的糙率。山區(qū)或郊區(qū)水庫下游糙率設置一般較為單一，通常取0.03左右[10]，而城市地區(qū)糙率值相對多元化，取值范圍更廣，一般在0.025～0.200范圍內。根據龍口水庫下游糙率信息，基于GIS制作龍口水庫下游河床糙率場，如圖2所示。

圖2 龍口水庫主壩下游地區(qū)糙率分布

圖3 潰口流量過程線

為了計算簡單方便，根據相關文獻[11]，本文只研究1種工況：假設龍口水庫在正常水位為72.0 m時因戰(zhàn)爭或地震造成主壩全部瞬間潰決。潰口寬度為大壩長度510 m，水庫庫容924萬m3，根據式(6)計算得潰口最大流量為33 174 m3/s，根據式(7)得到流量過程線如圖3所示。由圖3可見，潰壩流量過程在初始時刻最大，隨時間逐漸減少，且流量變化梯度較大，潰壩過程歷時較短，約23 min庫容基本被排空。

4.3 城市潰壩洪水演進模型構建

利用MIKE21的水動力模型模擬位于城區(qū)的龍口水庫的潰壩洪水演進過程，建立的二維模型(可參考圖1(a))設置2個上邊界和1個下邊界共3個開邊界，閉邊界條件一般均取為無滑移邊界條件，即法向流速為0。上邊界1為龍口水庫大壩潰口，潰口寬度假設為510 m(全潰)，按流量邊界控制；上邊界2設置在龍崗河上游(梧桐河與大康河交界處)，按流量邊界控制，并賦予5 m3/s的流量。下游開邊界位于龍崗河下游進入惠州邊境處，按水位邊界控制。閉邊界根據水庫下游可能淹沒范圍劃分，一般以高山、高地和高聳建筑物作為不過水邊界，閉邊界不與外界進行水量交換。開邊界與閉邊界共同組成一個閉合、狹長的研究區(qū)域，面積約44.63 km2，具體邊界設置可見圖1(a)。根據龍口水庫下游邊界情況劃分網格，本次模擬設定為三角形網格。網格劃分后，節(jié)點數175 048個，網格總數347 354個，最大網格200 m2，最小網格26.78 m2，網格平均大小為128.5 m2。

打開模型界面后，輸入處理好的地形、糙率、潰口流量數據，時間步長需根據計算網格的大小、水深、流速及實際計算的收斂性和穩(wěn)定性確定。由于研究范圍較大，生成網格較多，時間步長設置為30 s，模擬時間為12 h，模擬潰決開始時刻為0時。調整相關參數，在既定的邊界條件和初始條件下進行模型調試，確定合適的參數后進行模擬計算及結果輸出。

4.4 模擬結果及分析

輸出模擬結果后，對下泄水量、流態(tài)、淹沒范圍、淹沒水深、個別特征點流速和水深變化進行分析，以了解龍口水庫潰壩洪水演進的特性。

4.4.1 水量平衡分析

為確保模型計算的合理性，對水量平衡的分析是有必要的，驗證方法是由潰口流量過程線計算出潰壩后某段時間的排水量，再統(tǒng)計下游淹沒區(qū)的模擬淹沒水量，分析兩者的誤差。由于歷時2.0～2.5 h洪水陸續(xù)流進龍崗河，因此只統(tǒng)計到第2.5 h。由表1可知，最大的相對誤差僅3.34%，相對誤差基本上小于5%，可認為2種統(tǒng)計方式的水量基本相等，在水量平衡角度上說明模型是合理的。

表1 下泄水量與淹沒水量對比分析

圖4 各時刻流態(tài)圖

4.4.2 流態(tài)分析

根據輸出的流速信息結合Teapot平臺合成流態(tài)圖，圖4為瞬潰后第1，3，5，10 min的流態(tài)圖。圖中橫縱坐標為高斯平面直角坐標(km)。潰決瞬間潰口處流速最大，可達到15 m/s，水流高速下泄到下游地區(qū)。隨著水流的不斷前進以及時間的推移，流速將逐漸減?。? min后，主流流速大于7 m/s，3 min后主流流速為5～7 m/s，5 min后主流流速3～5 m/s，10 min后主流流速1.5～3.0 m/s，20 min后主流流速1.0～1.5 m/s，而過了30 min主流水流流速基本小于1.0 m/s。盡管水流在復雜的城市地形演進，但系列流態(tài)圖表明流態(tài)整體規(guī)律性強，流向合理不紊亂，流速值均勻無突變，無發(fā)散跡象，遵循水往低處流的客觀規(guī)律，說明模型是合理的。

4.4.3 特征點水深變化分析

為了解城市地區(qū)內部水位變化情況，可以選取典型特征點進行分析。本次按照離潰口由近到遠的原則依次選取6個典型特征點：壩腳點(特征點1)、龍崗區(qū)體育中心(特征點2)、體育中心會所(特征點3)、大運文化園(特征點5)、龍格路黃閣路交叉口(特征點13)和龍翔大道如意路交叉口(特征點15)，其具體位置見圖1(b)。由圖5可知，距離潰口越近，淹沒水深越高，特征點1，2，3，5因受到附近高大建筑物的雍水作用水深迅速增長，隨后迅速消落；特征點13，15附近沒有建筑物，水深上漲及消退相對平緩，由此可見建筑物對洪水雍高起到重要作用。

圖5 典型特征點水深變化

4.4.4 最大淹沒水深分析

最大淹沒水深是指潰壩洪水波在演進過程中所達到的最大淹沒水位，是一個衡量洪水淹沒嚴重程度的指標，最大淹沒水深越大，破壞越嚴重；反之，破壞相對越小。提取洪水最大淹沒水深可以幫助了解城區(qū)被潰壩洪水肆虐的重災區(qū)分布。利用水動力模塊統(tǒng)計出各個網格最大淹沒水深及淹沒范圍如圖6所示。

圖6 龍口水庫潰壩洪水最大淹沒水深及淹沒范圍

由圖6可知，受建筑物和公路街道的影響，城市地區(qū)淹沒產生明顯脈絡分布，局部建筑物間的最大淹沒水深較大。產生這些特點主要原因是由于洪水流經的短時期內，洪水遇到不透水建筑物需繞道而行；相對建筑物，道路與街道的地勢則更平坦，受到阻力較小，路面高程相對兩旁低，容易引導洪水前進；當來勢洶洶的洪水進入城區(qū)后，街道或道路兩側高聳的建筑物側向阻擋使得洪水涌高，若演進再遇到街道逐漸變窄或建筑物阻擋則雍水效應更加明顯，從而使得街道或道路兩側局部地區(qū)最大淹沒水深加大。

4.4.5 淹沒區(qū)統(tǒng)計分析

對淹沒區(qū)水深進行統(tǒng)計，下游淹沒面積共10.736 km2，最大淹沒水深超過3 m的占56.86%。龍城街道淹沒面積為9.915 km2，占總受淹面積92.35%，其中荷坳、愛聯2個居委會受淹面積分別為3.436，5.731 km2，占總受淹面積32.00%，53.38%。結合圖6分析可知，荷坳和愛聯居委會為受災重災區(qū)，這些區(qū)域盡管建筑較少，但離水庫潰口較近，且地勢低洼，受到洪水影響的風險最大。

5 結 語

以龍口水庫為例，構建了城市地區(qū)水庫潰壩洪水演進二維模型，通過對水量平衡、流態(tài)合理性分析證明模型合理可行。模擬結果表明：

(1) 利用GIS技術生成的高仿真度城市數字高程模型(DEM)，能較好模擬潰壩水流發(fā)展特性。

(2) 根據城市地區(qū)不同的地表類型構建多元的糙率場，能較好解決城市多元糙率的問題。

(3) 若龍口水庫主壩發(fā)生潰壩，荷坳、愛聯居委會將會是重災區(qū)，有關部門可根據模擬結果制定潰壩事故早期預警以及人員疏散方案。

本研究還存在一些不足，例如為計算簡便，沒有考慮潰口發(fā)展特性而采取了瞬潰假設，忽略了建筑物內部儲水作用等，因此在日后研究中尚需進一步完善。

參考文獻：

[1] 劉玉玲，王玲玲，周孝德，等.二維潰壩洪水波傳播的高精度數值模擬[J].自然災害學報，2010，19(5)：164-169.(LIU Yu-ling, WANG ling-ling, ZHOU Xiao-de,etal.High-resolution Numerical Simulation of 2D Dam-break Flood Waves[J].Journal of Natural Disasters, 2010, 19(5):164-169.(in Chinese))

[2] GALLEGOS H A，SANDERS B F，SCHUBERT J E.Two-dimensional, High-resolution Modeling of Urban Dam-break Flooding: A Case Study of Baldwin Hills, California[J].Advances in Water Resources，2009，32(8)：1323-1335.

[3] SHAKIBAEINIA A，JIN Y C.A Mesh-free Particle Model for Simulation of Mobile-bed Dam Break[J].Advances in Water Resources，2011，34(6)：794-807.

[4] PETACCIA G，SOARES-FRAZO S，SAVI F，etal.Simplified Versus Detailed Two-dimensional Approaches to Transient Flow Modeling in Urban Areas[J].Journal of Hydraulic Engineering，2010，136(4)：262-266.

[5] 王曉玲，張愛麗，陳華鴻，等.三維潰壩洪水在復雜淹沒區(qū)域演進的數值模擬[J].水利學報，2012，43(9)：1025- 1033.(WANG Xiao-ling, ZHANG Ai-li, CHEN Hua-hong,etal.Three-dimensional Numerical Simulation of Dam-break Flood Routing in the Complex Inundation Areas[J].Journal of Hydraulic Engineering, 2012, 43(9):1025-1033.(in Chinese))

[6] 謝任之.潰壩水力學[M].濟南：山東科學技術出版社，1993.(XIE Ren-zhi.Dam-break Hydraulics[M].Jinan: Shandong Science and Technology Press, 1993.(in Chinese))

[7] 張 瑋，鐘春欣，應翰海.草皮護坡水力糙率實驗研究[J].水科學進展，2007，18(4)：483-489.(ZHANG Wei, ZHONG Chun-xin, YING Han-hai.Experimental Study on Hydraulic Roughness of Revetment with Grass Cover[J].Advances in Water Science, 2007, 18(4):483-489.(in Chinese))

[8] 曾 祥，黃國兵，段文剛.混凝土渠道糙率調研綜述[J].長江科學院院報，1999，16(6)：1-4.(ZENG Xiang, HUANG Guo-bing, DUAN Wen-gang.Summarization of Investigation on Roughness of Concrete Channels[J].Journal of Yangtze River Scientific Research Institute, 1999, 16(6):1-4.(in Chinese))

[9] 落全富，安莉娜.青山水庫潰壩洪水模擬計算[J].浙江水利科技，2010，(2)：17-19.(LUO Quan-fu, AN Li-na.Simulated Calculation of Dam Break Flood for Qingshan Reservoir[J].Zhejing Hydrotechnics, 2010，(2):17-19.(in Chinese))

[10] 魏文禮，沈永明，孫廣才，等.二維潰壩洪水波演進的數值模擬[J].水利學報，2003，(9)：43-46.(WEI Wen-li, SHEN Yong-ming, SUN Guang-cai,etal.Numerical Simulation of 2D Dam-break Flood Wave[J].Journal of Hydraulic Engineering，2003，(9):43-46.(in Chinese))

[11] 郭洪巍，吳蔥蔥.逐漸潰壩和瞬時潰壩的模擬研究與比較[J].東北水利水電，2000，18(2)：1-3.(GUO Hong-wei, WU Cong-cong.Comparison and Research between Gradual and Instantaneous Dam Break Mathematical Model[J].Water Resource & Hydropower of Northeast, 2000, 18(2):1-3.(in Chinese))