石 莎，范子武，張 銘，費香波，烏景秀

(南京水利科學研究院，江蘇 南京 210029)

浯溪口水利樞紐位于景德鎮(zhèn)市蛟潭鎮(zhèn)境內(nèi)，距景德鎮(zhèn)市區(qū)40 km，是昌江干流中游一座以防洪為主，兼顧供水、發(fā)電等的綜合利用工程．水庫總庫容為4．27億m3，大壩壩長538．4 m，最大壩高45．6 m，正常蓄水位56 m，設(shè)計洪水位62．3 m，校核洪水位64．3 m．工程等別為二等，非溢流壩、溢流壩、河床式廠房(擋水部分)等主要建筑物的級別為2級，次要建筑物為3級，臨時建筑物為4級．混凝土壩防洪標準采用100年一遇洪水設(shè)計，1 000年一遇洪水校核;土壩防洪標準采用100年一遇洪水設(shè)計，2 000年一遇洪水校核．

浯溪口水利樞紐工程壩址下游的主要防洪保護對象是景德鎮(zhèn)市．景德鎮(zhèn)市2010年城市建設(shè)用地面積72．84 km2，總?cè)丝?6．21萬，工業(yè)總產(chǎn)值369．61億元．浯溪口水利樞紐建成后，通過浯溪口水庫與景德鎮(zhèn)市城防工程聯(lián)合調(diào)度，可以使景德鎮(zhèn)市區(qū)的防洪標準提高到50年一遇設(shè)防標準．但是，水庫大壩自身安全性所導(dǎo)致的潰壩洪水將給下游地區(qū)帶來潛在風險［1-2］．因此對浯溪口水利樞紐潰壩洪水演進的模擬計算對于其下游防洪，特別是景德鎮(zhèn)市的防洪管理決策具有十分重要的意義．

基于潰壩水流的破壞性和潰壩后果的嚴重性，世界各國高度重視大壩的安全問題．自19世紀以來，各國學者分別從理論分析、物理模型試驗、數(shù)值模擬、潰壩洪水演進模型研究、歷史資料統(tǒng)計分析等方面對潰壩問題進行了系統(tǒng)研究［3］．其中，物理模型試驗和數(shù)值模擬是近幾十年來應(yīng)用較多的研究方式．然而，潰壩洪水的物理模型試驗耗資高，且操作上具有一定的技術(shù)難度，因而目前數(shù)值模擬在潰壩洪水演進方面是較常用的研究手段．潰壩的速度和程度決定了潰壩洪水的大小和傳播速度，從而影響了潰壩洪水的影響范圍和程度［4］．潰壩水流的流態(tài)十分復(fù)雜，在洪水演進過程中，通常伴隨著急流、緩流和臨界流的情況［5］．因而構(gòu)建合理的潰壩洪水演進數(shù)學模型，計算潰壩洪水的淹沒范圍、水深和到達時間等要素，能為建立潰壩風險圖和評估潰壩損失提供科學依據(jù)，從而對水庫和堤防失事影響做出定量估算，并合理確定水庫或堤壩防洪設(shè)計標準及避險措施［6］．

TELEMAC-2D是一款基于有限單元法求解二維淺水方程的開源軟件，可以用于研究風暴潮、河流海岸動力學、波的傳播、泥沙和污染物輸移等，其網(wǎng)格劃分為非結(jié)構(gòu)化三角形網(wǎng)格，對于重要區(qū)域可進行局部加密．TELEMAC-2D在西歐、加拿大和美國等地區(qū)與國家河流模擬中應(yīng)用較為廣泛［7-12］．此外，加拿大學者D．Lam等在所研究的多模型集成決策支持系統(tǒng)中應(yīng)用TELEMAC-2D模擬湖流部分［13］．德國學者將TELEMAC-2D應(yīng)用于洪水污染研究、海嘯對德國沿岸城市影響研究以及河口支流水流流態(tài)研究等［14-16］．西班牙學者應(yīng)用TELEMAC-2D研究了沉水植物對水流阻力及渦黏性的影響［17］．荷蘭學者應(yīng)用TELEMAC-2D研究了Belgium的Ijzer河口的水動力學及泥沙輸移［18］．法國學者應(yīng)用TELEMAC-2D對魚道水力學特性進行了研究［19］．英國學者將TELEMAC-2D作為有限單元法代表性軟件與有限差分法進行了對比性研究［20］．TELEMAC-2D在中國的應(yīng)用則比較少，可查文獻有中國學者Zhang Cheng應(yīng)用TELEMAC-2D對云南洱海的二維湖流進行研究［21］，童朝鋒等應(yīng)用TELEMAC-2D系統(tǒng)建立長江口鹽水數(shù)值模型計算了不同徑流、潮汐條件下，崇啟大橋橋墩區(qū)域的潮位和鹽度變化過程［22］．

本文應(yīng)用TELEMAC-2D建立浯溪口下游洪水演進二維數(shù)學模型，對浯溪口壩址下游河道洪水演進進行模擬，采用河道水面線和河道斷面水位流量過程線驗證模型合理性，并且將幾組不同重現(xiàn)期的潰壩流量作為上游邊界條件，對浯溪口水庫下游二維潰壩洪水演進進行了數(shù)值模擬，計算了昌江河道不同斷面的流量水位變化過程、洪水是否漫堤和洪水到達的時間等要素，為洪水風險評估提供了科學依據(jù)．

1 數(shù)學模型

TELEMAC系統(tǒng)基于有限元數(shù)值求解程序，包含了明渠流的數(shù)字繪圖等普通功能．研究范圍覆蓋了水動力學、波、輸移質(zhì)、地下水和水質(zhì)等．TELEMAC-2D作為TELEMAC系統(tǒng)中的一個計算模塊，其計算基于Saint-Venant方程組，由Navier-Stokers方程組垂向積分再取平均所得．方程采用有限單元法進行求解，其中對于非線性項采用了相應(yīng)的假定和近似法［23］．應(yīng)用TELEMAC-2D所建立的浯溪口水利樞紐下游二維潰壩洪水演進數(shù)值模型如下．

1．1 基本方程

基本方程為垂向平均的Navier-Stokers方程組．

式中:h為水深;u，v分別為x，y方向的速度;g為重力加速度;νt為動力擴散系數(shù);z為自由表面高程;t為時間;Sh為水流源或匯;Sx，Sy分別為動力方程的源或匯;div(*→)為散度;▽→(*)為梯度．

1．2 邊界條件確定及參數(shù)選取

TELEMAC-2D中采用關(guān)鍵字定義計算參數(shù)、邊界條件和初始條件等．對于浯溪口潰壩洪水二維演進數(shù)學模型，邊界條件類型數(shù)據(jù)存儲于文件名為*．lci的邊界條件文件中，4個整型:LIHBOR，LIUBOR，LIVBOR和LITBOR，它們分別有0～6共7個取值，每個取值代表不同的邊界類型．此文件由JENAT軟件生成．TELEMAC-2D的求解是基于特征曲線法，應(yīng)用中要求在流體邊界的每個點都給出變量h，u，v的具體情況，TELEMAC-2D中的原則是:(1)急流入口，給定流速和水深;(2)緩流入口，給定流速和自由水深;(3)急流出口，自由流速和自由水深;(4)緩流出口，自由流速和給定水深．在浯溪口二維潰壩洪水數(shù)學模型中，上邊界為流量過程邊界，LIHBOR=4(自由水深的開邊界)，LIUBOR/LIVBOR=5(給定流量過程的開邊界)，此處不考慮污染物情況，因此LITBOR不進行賦值．下邊界為水位邊界，LIHBOR=5(給定水深的開邊界)，LIUBOR/LIVBOR=4(自由流速的開邊界)，同樣不考慮污染物情況，因此LITBOR不進行賦值．

而邊界條件中的具體數(shù)據(jù)，分別存儲于文件名為*．lbf的文件(上邊界為流量過程邊界:不同重現(xiàn)期下浯溪口潰壩洪水流量過程)中和文件名為*．qz的文件(下邊界為水位過程:昌江下游斷面的水位-流量關(guān)系特性曲線)中．計算時分別在配置文件(*．str)中取關(guān)鍵字LIQUID BOUNDARIES FILE=*．lbf，STAGEDISCHARGE CURVESFILE=*．qz和STAGE-DISCHARGE CURVES=0(或1)，即可調(diào)用文件中準備好的流體邊界數(shù)據(jù)．關(guān)于固壁邊界關(guān)鍵字取值為LAW OF BONTTOM FRICTION=2，表示采用Chezy公式，F(xiàn)RICTION COEFFICIENT=50．初始條件對應(yīng)關(guān)鍵字取為INITIAL DEPTH=0．01計算時步采用變步長，初始時步給定為TIME STEP=10，采用變時步關(guān)鍵字為VARIABLE TIME-STEP=YES，總計算時長根據(jù)潰壩流量過程定，關(guān)鍵字取值為DURATION=500 000(計算時步所有取值單位均為s)．

1．3 計算區(qū)域確定及網(wǎng)格劃分

如圖1所示，計算區(qū)域以60 m等高線為界，模型上游入口為浯溪口大壩壩址處，下游出口距上游入口62 km．渡峰坑位于景德鎮(zhèn)市西郊墾殖場莊下村，景德鎮(zhèn)城區(qū)在渡峰坑上游約3．7 km處．

本文采用JANET軟件進行網(wǎng)格劃分和網(wǎng)格質(zhì)量檢驗．全區(qū)域為三角形單元，河道局部加密，滿足河道橫斷面至少5個計算節(jié)點的要求．整個區(qū)域總計28 542個節(jié)點，56 428個單元．其中，最大長度673．6 m，最小長度21．5 m;最大面積158 848 m2，最小面積252 m2;最大角度123°，最小角度19°(見圖1)．考慮到需要對各節(jié)點的高程數(shù)據(jù)進行插值，根據(jù)已知的等高線密度確定網(wǎng)格密度;JANET對網(wǎng)格質(zhì)量檢驗合格，一個工況的計算耗時約2．5 h．

圖1 模型區(qū)域及模型網(wǎng)格Fig．1 Model area and model mesh

2 模型驗證

模型驗證分為2個部分:一是以測量斷面同期的實測水位為邊界條件，對昌江縱斷面水面線與實測水面線進行對比驗證;二是對渡峰坑斷面的水位流量關(guān)系曲線進行驗證．

2．1 水面線驗證

2009年3月至5月，景德鎮(zhèn)市水文局在浯溪口壩址至渡峰坑的昌江段進行了斷面測量，同時獲得了對應(yīng)時期的斷面水位資料．本文以斷面測量同期的實測水位為邊界控制，計算得到昌江縱斷面水面線，與實測水面線進行對比(見圖2)，浯溪口壩址位于106 km處，下游延續(xù)到渡峰坑約62 km處．從圖中可見，計算水面線和實測水面線之間有良好的擬合關(guān)系．上游偏差較大，是由于上游約100 km處有一攔水壩，對過水有短期的影響，造成實測水面線有比較明顯的水頭差．下游67 km左右有個水位陡降的過程，這是由于此處(渡峰坑上游不遠處)有一支流匯入，對上游來水有一定頂托作用．但由于本次驗證采用水位邊界，流量較小，當流量較大時，昌江整體水位抬高，攔水壩的阻水作用將大幅減?。畯尿炞C結(jié)果看，浯溪口二維數(shù)值模擬所采用的基礎(chǔ)數(shù)據(jù)，包括地形、斷面均合理可靠．

2．2 渡峰坑斷面水位流量關(guān)系曲線驗證

圖2 昌江水面線對比Fig．2 Water surface profiles of the Chang River

渡峰坑水文站設(shè)立于1941年，流量資料最為齊全．因此模型采用渡峰坑斷面水位-流量關(guān)系進行驗證．上游邊界為流量邊界，下游取出口處水位-流量關(guān)系為下邊界，進行模型計算．從計算結(jié)果中提取渡峰坑斷面各個時刻的水位和流量數(shù)據(jù)，以流量為橫軸、水位為縱軸繪制水位-流量關(guān)系，與渡峰坑斷面的實測水位流量過程進行對比(見圖3)．由圖3可見，計算所得的水位-流量關(guān)系離散點均勻分布于實測水位流量關(guān)系點附近，并且?guī)缀鯇崪y水位-流量關(guān)系線包括其中．證明下游邊界所取的水位流量關(guān)系是合理可靠的．由以上驗證結(jié)果可知，TELEMAC-2D系統(tǒng)應(yīng)用于浯溪口潰壩洪水演進二維數(shù)值模擬計算是合理有效的．

圖3 渡峰坑斷面水位-流量關(guān)系散點分布Fig．3 Water level-discharge curves

3 模型應(yīng)用

3．1 設(shè)計洪水組合

由設(shè)計暴雨推求設(shè)計洪水的方法得到浯溪口水庫的入庫設(shè)計洪水，再與水庫不同特征水位(設(shè)計洪水位62．3 m、校核洪水位64．3 m)進行組合，見表1．當水位超過設(shè)計洪水位，水庫按泄洪能力泄流．采用改進的BREACH模型，考慮不同入庫洪水重現(xiàn)期和特征水位，模擬計算出各種組合下浯溪口土壩壩段漫頂潰決潰口處流量過程，見圖4，以此作為下游潰壩洪水演進數(shù)值模擬的上游邊界條件．

表1 模型計算組合Tab．1 Calculation combinations of the model

圖4 不同工況下潰壩洪水過程Fig．4 Dam-break flood routing under different conditions

3．2 各種組合的計算結(jié)果分析

3．2．1 淹沒特征值 表2為12種組合情況下大壩下游洪水演進特征值．由表可見，設(shè)計水位條件下的最大流速和最大水深均較校核洪水位時大，這是由于校核洪水位的潰壩發(fā)生時間較早，入庫洪水還未達到洪峰流量時，潰壩洪水已經(jīng)發(fā)展至洪峰，而入庫洪水洪峰傳遞至壩址位置時，潰壩下泄流量已經(jīng)錯過洪峰值，此時的入庫洪水流量與潰壩下泄流量進行疊加所產(chǎn)生的最大洪峰流量反而小于二者在設(shè)計洪水位時所對應(yīng)的疊加洪峰流量．由圖4中不同工況下潰壩洪水過程也可看出，潰壩洪水洪峰值會隨著洪水重現(xiàn)期的增大而增大并且在發(fā)生時間上有所提前．此外，5 000年一遇校核洪水位和設(shè)計洪水位對應(yīng)的最大淹沒水深、最大流速和最大流量等特征數(shù)據(jù)均小于2 000年一遇．原因是當遭遇5 000年一遇特大洪水時，入庫洪水與潰壩洪水的錯峰更加明顯，使得二者的洪量疊加所產(chǎn)生的洪峰小于遭遇2 000年一遇洪水時二者的疊加值．

表2 洪水演進特征參數(shù)Tab．2 Flood routing characteristic parameters

3．2．2 淹沒范圍 圖5為組合3、9情況下浯溪口土壩壩段潰壩洪水的演進過程，黑色部分表示未被洪水淹沒的陸地，白色部分為水淹部分，白色與黑色的遞進表示水深逐漸減小．由圖可見，下游淹沒范圍較大，景德鎮(zhèn)城區(qū)也在淹沒范圍內(nèi)，且積水較多．此外，由圖5還能看出，隨著時間的推移，洪峰流量由上游向下游傳播，潰壩1 h至20 h小時過程中，潰壩水流流量越來越大．

圖5 兩種組合下淹沒范圍對比Fig．5 Comparison of the flooded downstream areas under two conditions

4 結(jié)語

本文應(yīng)用TELEMAC-2D對浯溪口下游二維潰壩洪水演進進行了數(shù)值模擬．通過對昌江縱斷面水面線和渡峰坑斷面的水位流量關(guān)系曲線驗證計算，表明模型所選取的參數(shù)是合理的，TELEMAC-2D系統(tǒng)能較好地模擬二維潰壩洪水演進過程．利用該模型模擬計算了浯溪口水利樞紐建成后，遭遇100年一遇、200年一遇、500年一遇、1 000年一遇、2 000年一遇、5 000年一遇的入庫洪水流量情況下，漫頂潰壩水流洪水演進的基本參數(shù)，結(jié)果顯示:

(1)由于入庫洪水與潰壩洪水疊加作用，時間上產(chǎn)生的錯峰，使得校核洪水情況下特征參數(shù)小于設(shè)計洪水下的特征參數(shù);洪水重現(xiàn)期高達5 000年一遇時，由于潰壩水流與入庫水流的疊加作用，錯峰明顯，以至其相關(guān)特征數(shù)據(jù)小于2 000年一遇．

(2)對于同一設(shè)計水位情況下，洪水重現(xiàn)期越長，潰壩洪水在下游演進過程中的最大流速、最大流量、最大淹沒水深和最大淹沒面積都越來越大．

(3)500年一遇來流情況產(chǎn)生潰壩的下游洪水最大水深為30．87 m，此時昌江局部地區(qū)出現(xiàn)了洪水漫堤現(xiàn)象，因此沿昌江河堤的防洪能力，特別是景德鎮(zhèn)城區(qū)的堤防防洪能力需要適當提高．

本文現(xiàn)階段主要研究了漫頂潰壩洪水演進，為浯溪口下游防洪風險管理提供了一定的可靠的基礎(chǔ)數(shù)據(jù)支持．若進一步研究滲透潰壩洪水演進，并將其結(jié)果與漫頂潰壩洪水進行對比，對防洪預(yù)報預(yù)警與洪水風險管理的意義更為重大．

［1］龍曉飛，高龍華．茜坑水庫潰壩洪水數(shù)值模擬研究［J］．人民珠江，2011(2):42-44．(LONG Xiao-fei，GAO Long-hua．Dam-break flood routing simulation of Qiankeng reservoir［J］．Pearl River，2011(2):42-44．(in Chinese))

［2］張利民，徐耀，賈金生．國外潰壩數(shù)據(jù)庫［J］．中國防汛抗旱，2007(增刊):2-7．(ZHANG Li-min，XU Yao，JIA Jinsheng．Foreign dam failure database［J］．China Flood and Drought Management，2007(Suppl):2-7．(in Chinese))

［3］李云，李君．潰壩洪水模型實驗研究綜述［J］．水科學進展，2009，20(2):304-310．(LI Yun，LI Jun．Review of experimental study on dam-break［J］．Advances in Water Science，2009，20(2):304-310．(in Chinese))

［4］朱勇輝，廖鴻志，吳中如．國外土壩潰壩模擬綜述［J］．長江科學院院報，2003，20(2):26-29．(ZHU Yong-hui，LIAO Hong-zhi，WU Zhong-ru．Review on oversea earth-dam-break modeling［J］．Journal of Yangtze River Scientific Research Institute，2003，20(2):26-29．(in Chinese))

［5］張大偉，李丹勛，王興奎．基于非結(jié)構(gòu)網(wǎng)格的潰壩水流干濕變化過程數(shù)值模擬［J］．水力發(fā)電學報，2008，27(5):98-102．(ZHANG Da-wei，LI Dan-xun，WANG Xing-kui．Numerical modeling of dam-break water flow with wetting and drying change based on unstructured grids［J］．Journal of Hydroelectric Engineering，2008，27(5):98-102．(in Chinese))

［6］史宏達，劉臻．潰壩水流數(shù)值模擬研究進展［J］．水科學進展，2006，17(1):129-135．(SHI Hong-da，LIU Zhen．Review and progress of research in numerical simulation of dam-break water flow［J］．Advances in Water Science，2006，17(1):129-135．(in Chinese))

［7］BRIèRE C，ABADIE S，BRETEL P，et al．Assessment of TELEMAC system performances，a hydrodynamic case study of Anglet，F(xiàn)rance［J］．Coastal Engineering，2007，54(4):345-356．

［8］ASARO G，PARIS E．The effects induced by a new embankment at the confluence between two rivers:TELEMAC results compared with a physical model［J］．Hydrological Processes，2000，14(13):2345-2353．

［9］WOLF J．Coastal flooding:impacts of coupled wave-surge-tide models［J］．Natural Hazards，2009，49(2):241-260．

［10］SYME W J，PINNELL M G，WICKSJM．Modelling flood inundation of urban areas in the UK using 2D/1D hydraulic models［C］∥8th National Conference on Hydraulics in Water Engineering，Australia．2004．

［11］JONESJE，DAVIESA M．Storm surge computations for the west coast of Britain using a finite element model(TELEMAC)［J］．Ocean Dynamics，2008，58(5-6):337-363．

［12］DI BALDASSARRE G，CASTELLARIN A，MONTANARI A，et al．Probability-weighted hazard maps for comparing different flood risk management strategies:a case study［J］．Natural Hazards，2009，50(3):479-496．

［13］LAM D，LEONL，HAMILTON S，et al．Multi-model integration in decision support system:a technical user interface approach for watershed and lake management scenarios［J］．Environmental Modelling＆Software，2004，19(3):317-324．

［14］BüTTNER O，SCHULZ M，MATTHIES M，et al．Flood and pollutant dispersal simulation in urban areas［C］∥Conference Proceedings of the International Conference“Science and Information Technologies for Sustainable Management of Aquatic Ecosystems”，2009:12-16．

［15］LEHFELDT R，MILBRADT P，PLüSSA，et al．Propagation of a tsunami-wave in the North Sea［J］．Die Küste，2007，72:105-123．

［16］MALCHEREK A．Application of TELEMAC-2D in a narrow estuarine tributary［J］．Hydrological Processes，2000，14(13):2293-2300．

［17］VIONNET C A，TASSI P A，MARTIN VIDE JP．Estimates of flow resistance and eddy viscosity coefficients for 2 D modelling on vegetated floodplains［J］．Hydrological Processes，2004，18(15):2907-2926．

［18］GIARDINO A，LBRAHIM E，ADAM S，et al．Hydrodynamics and cohesive sediment transport in the Ijzer Estuary Blgium:Case Study［J］．Journal of Waterway，Port，Coastal and Ocean Engineering，2009，135(4):176-184．

［19］CHORDA J，MAUBOURGUET M M，ROUX H，et al．Two-dimensional free surface flow numerical model for vertical slot fishways［J］．Journal of Hydraulic Research，2010，48(2):141-151．

［20］JONES J E，DAVIES A M．An intercomparison between finite difference and finite element(TELEMAC)approaches to modelling west coast of Britain tides［J］．Ocean Dynamics，2005，55(3-4):178-198．

［21］ZHANG Cheng，GUILBAUD C．2D hydrodynamic modelling for Erhai Lake applying TELEMAC modelling system［EB/OL］．http:∥www．paper．edu．cn．2008．

［22］童朝鋒，劉豐陽，邵宇陽，等．長江口北支崇啟大橋處潮位和鹽度過程研究［J］．水道港口，2012，33(4):291-298．(TONG Chao-feng，LIU Feng-yang，SHAO Yu-yang，et al．A modeling study on tidal and salinity process at Chongqi bridge cross section located in north branch of Yangtze Estuary［J］．Journal of Waterway and Harbor，2012，33(4):291-298．(in Chinese))

［23］HERVOUET JM．Hydrodynamics of free surface flows-modelling with the finite element method［M］．West Sussex:John Wiley＆Sons Ltd，2007．