俞增鑫，漆文邦，2，張 強(qiáng)，吳倩倩

(1.四川大學(xué)水利水電學(xué)院，四川 成都 610065；2.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065；3.南昌工程學(xué)院瑤湖學(xué)院，江西 南昌 330099)

我國是世界上建成水庫大壩數(shù)量最多的國家，據(jù)中國大壩協(xié)會統(tǒng)計，目前我國已建成水庫大壩98795座，總庫容9035億m3，其中中、小型水庫占94%以上[1]。這些水庫大壩的安全不僅關(guān)系興利更是關(guān)系到下游區(qū)域的財產(chǎn)和生命安全，因此有必要對一些大壩開展?jié)窝芯糠治觥：芏鄬W(xué)者對潰壩水流數(shù)值模擬進(jìn)行深入研究并取得了大量的研究成果。在一維模擬研究方面，Harten[2]提出了一種高分辨率的TVD格式并將其與歐拉方程結(jié)合起來，可以在復(fù)雜的條件下捕捉到激波從而得出精度較高的解析解；張永祥[3]對一維潰壩洪水波進(jìn)行數(shù)值模擬，發(fā)現(xiàn)運(yùn)用時空守恒法可以很好地對潰壩洪水波進(jìn)行模擬，還有對潰壩洪水波有較強(qiáng)地捕捉性能、具有較高精度以及可以顯示潰壩洪水波的運(yùn)動特性等特點(diǎn)；在二維模擬研究方面，朱琳[4]將Level-Set的方法引入到二維潰壩問題的研究上來，將高階迎風(fēng)格式和二階投影的方法結(jié)合起來對二維N-S方程進(jìn)行求解，采用CIP方法來解Level-Set方程組，計算求得潰壩運(yùn)動界面不同時刻的變化情況及其相應(yīng)的速度分布；劉玉玲[5]建立了高分辨率的二維淺水方程的數(shù)學(xué)模型，運(yùn)用二階TVD格式確定有限體積交界面處的物理通量，模擬矩形斷面明渠大壩瞬潰情況下的洪水波的運(yùn)動特性和演進(jìn)過程。

但很多研究都是從單一維度的理想狀態(tài)為基礎(chǔ)進(jìn)行研究，因此可能存在一定的局限，如：一維模型不能模擬全部流域的情況，不能直觀地反映實(shí)際洪水演進(jìn)情況；二維模型穩(wěn)定性不好，無法顯示垂直方向的流態(tài)參數(shù)以及對時間步長的選擇要求較高等。MIKE Flood模塊可以將一維、二維耦合模擬計算，可減少單一維度的局限性，提高計算精度，可以很好地模擬潰壩洪水的演進(jìn)情況，但探討入庫洪水、潰決歷時和終潰高度的組合對大型水庫工程潰壩洪水演進(jìn)的影響較少?；诖耍疚囊阅乘畮旃こ虨槔?，綜合考慮洪水大小、潰壩歷時和終潰高度對洪水演進(jìn)的影響，模擬分析下游河道和城鎮(zhèn)區(qū)域的水力特征情況，考慮工程最不利工況的淹沒情況。

1 模型介紹

1.1 一維模型控制方程

一維非恒定水動力學(xué)模型的基本方程可以用在服從不可壓縮、均質(zhì)流體、河床坡度小、流態(tài)為一維的水流等假定的圣維南方程組來描述，其控制方程為：

式中：x為計算點(diǎn)的流程，m；Q為過流流量，m3/s；Z為水位，m；g為重力加速度，m/s2；A為過水?dāng)嗝婷娣e，m2；q為側(cè)向單位長度的入流流量，m2/s；t為計算點(diǎn)時間坐標(biāo)，s；R為水力半徑；為動量修正系數(shù)；n為粗糙系數(shù)。

方程組的求解采用Abbott-Ioncscu六點(diǎn)隱式有限差分格式來解[6]。這種格式的優(yōu)點(diǎn)是能讓水位和流量交替計算，這樣可以使得計算的穩(wěn)定性更高，提高計算速度和精度等。

1.2 二維模型

平面二維淺水水流控制方程如下：

式中：h是水深，m；x、y、z是空間和時間坐標(biāo)；u和v為和方向的流速分量，m/s；M、N分別為在x和y方向單寬流量的分量，m2/s；M=hu，N=hv；g為重力加速度，m/s2；zb為河床高程，m；n為糙率系數(shù)；γt為紊動粘滯系數(shù)。

此次研究在對模型控制方程的計算過程中方程離散時空間采用有限體積法，運(yùn)用守恒格式對水流連續(xù)方程進(jìn)行離散，保證計算域內(nèi)水量守恒，時間采用蛙跳法。模型采用非結(jié)構(gòu)化三角形網(wǎng)格，三角形網(wǎng)格在地形條件和邊界條件更復(fù)雜的情況下可以更好地模擬邊界情況，在地形變化的反映程度上更好[7]，這樣可以便于計算機(jī)對研究區(qū)域的劃分，所需數(shù)據(jù)簡單，計算速度更快，效果較矩形網(wǎng)格好。

1.3 耦合模型

MIKE Flood是MIKE系列軟件中的一個模塊，可以將MIKE11的一維模型和MIKE21的二維模型耦合連接起來進(jìn)行洪水演進(jìn)模擬。一維河道模型和二維地形模型的耦合連接方式有標(biāo)準(zhǔn)連接和側(cè)向連接，本文采用側(cè)向連接[8]，側(cè)向連接的連接方式見圖1。

圖1 側(cè)向連接示意圖

2 模型構(gòu)建

2.1 工程概況

某水庫位是一個具有多年調(diào)節(jié)性能的以發(fā)電為主，兼顧防洪和其他綜合效益的大型水利水電工程。該水電站壩址以上流域面積3330 km2，干流流域占比93%，平均坡降0.253%，河流長度187 km，河面寬度為100 m～300 m。流域較狹窄，長105 km，平均寬度約為31.7 km。大壩為鋼筋混凝土面板堆石壩，壩頂高程為171 m，最大壩高為162 m，壩頂長為507 m。水庫校核洪水位169.15 m，總庫容41.9億m3；正常蓄水位160.0 m，相應(yīng)庫容35.2億m3；防洪高水位161.5 m；汛限水位156.5 m；死水位120.0 m，調(diào)節(jié)庫容21.26億m3。壩址及斷面區(qū)域在概化流域內(nèi)的相對位置見圖2。

圖2 壩址及斷面相對位置示意圖

2.2 建模參數(shù)

2.2.1 糙率

根據(jù)工程相關(guān)資料，流域內(nèi)的糙率變化范圍為0.04～0.065之間，參考該工程施工圍堰潰堰洪水專題研究報告中的取值，取平均糙率為0.04。

2.2.2 潰口形狀

將潰口形狀概化為等腰梯形，潰口的尺寸大小由最終潰口底寬、潰口形狀參數(shù)和最終潰口底部高程決定。參考相關(guān)文獻(xiàn)資料[9]，潰口模型中的平均潰口寬度取值范圍的經(jīng)驗(yàn)值為，壩前水面寬，近似取壩寬。最終潰口底寬b與平均潰口寬度關(guān)系[10]為：

式中：hd是潰口以上水深，m；形狀參數(shù)取值范圍一般為0≤m≤2，與壩體材料有關(guān)，本文取定m=2；最終潰底高程hm一般為水庫的底部高程或者泄水建筑物的底部高程，考慮前期施工強(qiáng)度較大，設(shè)定一期澆筑高程98 m和泄水口高程49 m為終潰高程。假定潰決終止時最終潰口底寬不變，根據(jù)上式，終潰高程為98 m時，取b=260 m；潰決到泄水孔高程時，取b=160 m。

2.2.3 潰口流量

模型采用逐漸潰的形式對潰壩水流進(jìn)行模擬，對于潰口流量的計算，假定潰口從小孔開始線性發(fā)展破壞，洪水流量參考文獻(xiàn)[11]按下式進(jìn)行計算：

Qmax=A[2g(H-Hp)/(1+fL/D)]0.5

式中：Qmax為潰決最大流量，m3/s；H為庫區(qū)水位，m；A為斷面面積，m2；Hp為潰口高度，m；f為摩擦系數(shù)，與壩體材料有關(guān)；L為潰口沿水流方向長度，m；D為寬度，m。

潰口流量過程按以下公式進(jìn)行計算：

Q=Qmax(Qmax/5Vt-1)

2.2.4 設(shè)計洪水頻率及潰壩歷時

考慮在潰壩模擬過程中不同頻率的來流量和不同的潰壩歷時對結(jié)果的影響，結(jié)合相關(guān)數(shù)據(jù)資料，設(shè)定來流洪水頻率為200 a一遇、500 a一遇，潰壩歷時假設(shè)為1 h和3 h。

2.2.5 起潰水位

大壩壩頂高程為171 m，設(shè)定大壩壩前初始水位為汛限水位156.5 m。對于起潰水位計算問題，外國學(xué)者Singh[12]研究分析多個潰壩事故，總結(jié)得出漫頂潰決的土石壩的洪水漫頂高程在0.15 m～0.61 m之間，考慮到電站為混凝土面板堆石壩且工程規(guī)模較大，防洪能力較一般土石壩強(qiáng)，參考相關(guān)文獻(xiàn)[13]，取潰壩水位為171.61 m。

2.3 計算方案

考慮不同因素對潰壩水流模擬結(jié)果的影響，設(shè)定計算方案見表1。

表1 計算方案參數(shù)表

計算方案中的所有工況都采用模型耦合計算，對于潰壩洪水的演進(jìn)情況，河道采用一維模型模擬，下游淹沒區(qū)域采用二維模型模擬，河道和開闊地區(qū)的流量和水位變化情況均可通過耦合模型的計算得出。

3 計算及結(jié)果分析

3.1 潰口流量分析

根據(jù)計算方案設(shè)定模型參數(shù)進(jìn)行模擬計算，計算得出四個工況下潰口的流量過程見圖3。從圖3中可以看出，潰決歷時和終潰高度對潰口流量影響較大：工況1和工況2，潰壩歷時分別為1 h和3 h，潰口流量對應(yīng)為260215.88 m3/s和225372.22 m3/s；工況3和工況4，終潰高度分別為98 m和49 m，潰口流量對應(yīng)為268521.06 m3/s和411805.19 m3/s。在該模型中，洪水大小對潰口流量的影響較小，考慮是此工程規(guī)模較大，這兩個工況下庫區(qū)下泄庫容遠(yuǎn)大于上游來流洪水量，因此削弱了上游洪水對潰口流量的影響。

圖3 各工況下潰口流量過程

從圖中可以看出，設(shè)計洪水大小對大壩潰決時間影響較大：500 a一遇洪水工況下大壩的潰決時間會比200 a一遇工況提前約1 h。且在工況1和工況2中，因?yàn)闈螝v時不同，歷時長工況下潰口流量的峰現(xiàn)時間推遲約2 h。

3.2 沿程水位及流量分析

選取壩址及河道沿程幾個重要城鎮(zhèn)所在斷面作為特征斷面，模擬計算得出的各特征斷面的最高水位和最大流量情況見圖4。

從圖4可知，潰壩洪水對下游的流量和水位的影響是沿程減小的；在該模型中，終潰高度對下游區(qū)域的流量和水位影響最大，潰決歷時次之：工況2中，青田縣和梅岱的洪峰流量為178003.4 m3/s和128959.4 m3/s，最高水位為38.55 m和12.86 m；工況4中兩地流量分別為283822.5 m3/s和189252.6 m3/s，最高水位為46.99 m和17.89 m。

(a)洪峰流量沿程變化

(b)最高水位沿程變化

3.3 淹沒范圍及水深分析

為更好研究下游區(qū)域淹沒情況，選取支流匯水口上游船寮、下游油竹和上戊三地研究潰壩洪水對城鎮(zhèn)的淹沒影響。三地在各工況下入口的流量過程見圖5(時間均以大壩潰決時間為基準(zhǔn))。

(a)船寮地區(qū)

(b)油竹地區(qū)

(c)上戊地區(qū)

從圖5中可以看出：船寮地區(qū)比油竹和上戊地區(qū)高程高，但其入洪時間早，洪峰流量大，考慮其地理位置離大壩近，潰壩后從支流雍到船寮入口的洪水量較大；在工況1和工況3中，三地入口處流量過程差異不大，各地的入洪時間和峰現(xiàn)時間相差不大，且洪峰流量較接近；而在工況4下各地的入洪時間以及洪峰流量相較于工況2均有較大差異：工況4中，船寮地區(qū)入洪時間為2.2 h，洪峰流量為21020.9 m3/s；油竹地區(qū)入洪時間為2.5 h，洪峰流量為12623.5 m3/s；上戊地區(qū)入洪時間為3.3 h，洪峰流量為9703.22 m3/s。在工況2中，相較于工況4，船寮地區(qū)入洪時間為3.6 h，推遲約1.4 h，洪峰流量為6813.01 m3/s，相差14207.89 m3/s；油竹地區(qū)入洪時間為4.2 h，推遲約1.7 h，洪峰流量為4009.07 m3/s，相差8614.43 m3/s；上戊地區(qū)入洪時間為5.4 h，推遲約1.8 h，洪峰流量為3320.2 m3/s，相差6382.9 m3/s。

根據(jù)前面分析可知，終潰高度越低，潰壩歷時越短，洪水流量越大，對下游區(qū)域的影響程度越大。根據(jù)工程最不利原則，本文只分析上述研究區(qū)域在工況4情況下的區(qū)域淹沒情況。各地在工況4下的最大淹沒情況見圖6。結(jié)合水深數(shù)據(jù)，船寮地區(qū)入洪約2.2 h后達(dá)到最高水位，最大水深約出現(xiàn)在區(qū)域中部，為24.3 m，平均水深約為15 m；油竹地區(qū)入洪約2.7 h后達(dá)到最高水位，最大水深出現(xiàn)在入口附近，為22.2 m，平均水深約為10.5 m；上戊地區(qū)入洪約2.9 h后達(dá)到最高水位，最大水深出現(xiàn)在入口附近，為9.47 m，平均水深約為5 m。

(a)船寮地區(qū)

(b)油竹地區(qū)

(c)上戊地區(qū)

4 結(jié)語

本文用MIKE Flood耦合模型對潰壩洪水的演進(jìn)情況進(jìn)行模擬計算，分析了潰壩歷時、終潰高度等因素對下游水位和流量的影響，研究下游重要城鎮(zhèn)的流量情況和淹沒情況。

(1)終潰高度變化造成的潰口流量以及下游區(qū)域流量和水位情況較其他因素明顯增加；大壩潰決歷時長，潰口洪峰流量峰現(xiàn)時間較歷時短推遲，下游區(qū)域流量和水位增量均減小。

(2)大壩潰決后，船寮地區(qū)離潰口最近，其入洪時間較另外兩地分別提前約30 min和1 h，受淹情況最嚴(yán)重，最大平均淹沒水深約15 m，應(yīng)該加強(qiáng)該地的避險轉(zhuǎn)移工作，其他兩地也應(yīng)做好避險和洪水排導(dǎo)的工作。