黃 燦 劉青泉 王曉亮
(北京理工大學宇航學院力學系,北京 100081)
我國水電資源豐富,當前已經建成各種類型的水庫接近10 萬座[1],西南地區的許多大型河流,如金沙江、大渡河等都已經形成了梯級水庫布局[2],一旦發生梯級潰決,將誘發大型潰壩洪澇災害,對下游沿岸的生產和生活帶來極大的影響.另一方面,地震災害活動經常誘發滑坡形成堰塞壩和堰塞湖,其穩定性和安全性存在極大的風險,也將帶來梯級潰決風險,如2018 年的白格堰塞體[3].因此研究梯級潰決洪水演進過程及其機理對維護梯級水庫安全運行具有重要意義.
近年來梯級潰決洪水過程受到了較大的關注,研究集中在室內實驗和數值模擬兩方面.Xue 等[4]開展了比較系統的室內梯級潰決室內實驗,主要探討了壩間距、庫容比對下游水庫水流洪峰的影響.Niu等[5]在斜槽中堆積了梯形大壩,研究了梯級大壩潰決形式和潰壩水流的傳播過程.Chen 等[6]通過實驗研究了潰壩水流對下游壩體的沖擊力作用,指出針對下游不同庫容存在3 種沖擊模式,并給出了作用力擬合公式.大量的工作集中在數值模擬方面,主要可以分為水文方法、淺水動力學模擬,以及三維Navier-Stokes 方程模擬.趙雪瑩等[7]使用經驗公式研究了貴陽柏山?花溪梯級水庫發生梯級潰決的洪水演進,分析了梯級水庫群不同水位組合總計4 種情況,誘發潰壩洪水對下游的影響,并指出降低庫水位將有效減輕潰壩洪水對下游影響.賀同坤等[8]使用Mike11 軟件模擬了瀾滄江中游五個水庫組成的梯級庫的特大洪水演進過程,指出當發生一萬年一遇洪水時,漫灣和景洪兩個水庫的最高水位接近壩頂.劉慶紅等[9]模擬了特定河流的梯級水庫群的連鎖響應以及向下游演進的過程.黃衛等[10]建立了梯級大壩潰決洪水過程的淺水動力學模型,研究了不同組合情況下梯級潰決洪水演進過程,指出梯級潰決洪峰存在增強機制.隨后Cao 等[11]進一步使用水沙耦合模型研究指出梯級潰決誘發洪水將更具威脅,當前的單壩潰決預案需要改進.Wang 等[12]比較了淺水動力學模型和Boussinesq 模型模擬梯級潰決洪水的傳播的爬升過程,并指出Boussinesq 在相對較大水深處能給出更加精確的解.Zhou 等[13]使用VOF 方法和湍流模型研究了梯級潰決洪水過程.Luo 等[14]用近年來快速發展的無網格粒子方法(SPH)模擬了梯級潰決過程,將潰壩洪水效應分為高速沖擊區和傳壓區兩個過程.李仟[15]研究了梯級土石壩潰決洪水的演進和調洪過程,對下游洪水過程進行了風險分析.袁岳[16]用RANS 模型和標準k–ε 模型模擬了不同坡角,壩間距及庫容比對潰決過程的影響,比較了單級潰決和梯級潰決不同潰決模式下水流流動特點.趙丹[17]研究上寺水庫大壩潰決的流量過程,并對當前使用的常見模型模擬精度進行了比較.
綜上所述,當前對梯級潰決開展了一些機理性實驗研究和各種不同類型的數值模擬,主要研究潰壩洪水過程,大部分都指出了梯級潰決相對單壩潰決存在增強現象,但是對梯級潰決洪峰增強機制缺乏基本的物理和力學解釋,特別是缺少質量和能量轉化如何增強梯級潰決洪水過程的認識.本文使用淺水動力學方法研究梯級潰決洪水過程,從梯級潰決上游來水對下游水庫的質量和能量補充及分布方面入手,對梯級潰決洪峰增強機制進行一些探索,并給出了一個射流水塔單壩潰決模式等效梯級潰決洪水過程.
土石壩的潰決過程十分復雜,涉及水流過程、泥沙沖刷以及壩體坍塌及其復雜的流固耦合問題[18-22],本文主要研究壩體潰決誘發的洪水過程,因此對壩體坍塌采取瞬間潰決的模式進行處理.雖然潰壩過程在局部水流結構復雜,但是潰壩洪水演化的過程基本發生在很長的河道,其水平尺度相比于水流的垂向尺度大得多,因此滿足淺水動力學的基本假設[23-24].本節主要給出潰壩洪水的基本控制方程和數值方法,并通過幾個經典案例驗證數值模型的有效性.
由于潰壩洪水過程的水平尺度遠大于垂向尺度,并且本文對潰壩過程采取瞬間潰決的方式處理,因此潰壩洪水演化過程可以用淺水動力學方程描述

其中,h和u分別表示水深和水平方向的垂向平均速度,g和b則分別表示重力加速度和地形高程.本文考慮大尺度洪水演進過程和機理研究,因此忽略摩擦作用.
淺水動力學模型(1)和(2)是一組非線性雙曲型偏微分方程,其典型的特征是會出現水躍(激波),同時由于源項的作用,存在一個平衡解,即所謂Well-Balance 條件,另一方面實際問題經常遇到一側水深為零的情況,即所謂的干濕條件,因此這套控制方程的求解算法需要具備至少3 個特征,捕捉激波、滿足Well-Balance 條件,以及捕捉干濕邊界.過去十年在非線性雙曲守恒律的高分辨率數值格式方面取得了重要的進展[25-27],本文綜合考慮了這3 方面的需求,結合近期淺水動力學計算方法的進展,融合了Kurganov的Riemann 解[28-30]和Audusse[31-32]的單元插值模式,構造了一套淺水動力學模型的離散格式.
首先將控制方程寫成守恒形式,如式(3)所示

其中,U=(h,hu)T,F=(hu,hu2+gh2/2)T,s(U)=[0,?gh(?b/?x)]T.
采用一套滿足 Well-Balance 行為的有限體積法[31]對守恒方程(3)進行離散,得到以下的離散格式,如式(4)所示

其中Uj=(hj,qj)T為第j個單元的單元平均值,數值通量Fj+1/2的構造采用MUSCL 方法[33],如式(5)所示

這里Uj+(1/2)?=(hj+(1/2)?,hj+(1/2)?uj,r)T,Uj+(1/2)+=(hj+(1/2)+,hj+(1/2)+uj+1,l)T.
采用滿足 Well-Balance 特性的插值方式構造hj+(1/2)?和hj+(1/2)+,其中hj+(1/2)?=max(0,hj,r+bj,r?bj+1/2),hj+(1/2)+=max(0,hj+1,l+bj+1,l?bj+1/2),這里bj+1/2=max(bj,r,bj+1,l).所有的單元左右重構值均采用二階MUSCL 方法構造,并采用Minmod 限制器抑制虛假振蕩[34].
這里的Riemann 近似求解器可以采用任意形式的Riemann 求解器,本文考慮計算穩定性,采用了近年來新發展的Kurganov 近似解[30].源項部分包括兩項,其表達式分別如式(6)和式(7)所示

時間離散采取二階Runge-Kutta 方法.
本節采取兩個典型的例子驗證1.2 中算法的有效性,兩個例子分別是潰壩問題和Well-Balance 問題,將給出兩個問題的數值解和解析解的比較.
1.3.1 潰壩問題簡單波解
經典潰壩問題以一側為水,另一側無水為初始情況,瞬間抽除擋板后的水流運動演化過程.該過程存在解析解,稱為簡單波解[35].設計初始水深h0=10 m,x=0 為初始擋板位置,本例的解析解為h=(h0/9)(2 ?x/ct)2,u=2(c+x/t)/3,其中波速,流量Q=hu.
圖1 給出了時刻0.072 1 s 的水深h和流量Q的計算結果,可以發現計算結果和解析解基本吻合.

圖1 潰壩問題在時刻0.0721 s 的水深和流量分布,數值解和理論解的比較,H 為水位,Q 為流量Fig.1 Height and discharge distribution of dam break flow at 0.0721 s,the red dashed line and blue symbols are analytical and numerical solutions respectively,where H is water depth,and Q is discharge
1.3.2 Well-Balance 行為
地形變化情況的淺水動力學方程存在一個穩態解結構,即h+b=const,q=0,即在流量為0 的時候,自由表面為常數.這個結構在數值模型中需要保持,稱為Well-Balance 問題.本例給出了底部為一個高度為5 m 的近似拋物線凸起,自由水面高程為10 m 例子的數值模擬.圖2 給出了本例在時刻10.093 2 s 的數值解,可以發現自由水面高程保持常數,同時流量在凸起位置處為很小的值,約10?5m2/s 量級,這在后續的工程計算中是可以接受的.

圖2 Well-Balance 問題的驗證,H+b 為水位,Q 為流量,綠色虛線為底床高程Fig.2 Verification for Well-Balance behavior over a bump denoted by red dashed line,where H+b is free surface elevation,and Q is discharge
梯級水庫是當前水資源開發利用中的常見模式,上游壩體的潰決往往又會導致下游壩體的潰決,其流動過程更為復雜.針對這種梯級水庫開發模式,本文模擬了典型的兩個梯級壩連續潰決的洪水演進過程,重點分析梯級潰決引發的洪峰增強現象.
如圖3 所示,在坡度為α 的河道上,上下游存在兩個壩體,其間距為L,兩個壩體圍成了兩個水庫,其庫容高度分別為h1和h2,梯級潰決過程是,首先上游壩體突然潰決,水流往下游傳播,下游壩前水位超過壩體高度后下游壩體瞬間潰決,形成潰壩洪水往下游傳播.主要關注潰壩洪水的演進過程及其對下游的影響,因此在下游處設置兩個測點P1和P2,測點P1距離下游水庫的距離為L1,測點P2距離測點P1的距離為L2.

圖3 梯級潰決問題描述以及測點Fig.3 Cascade dam-break model and detection points P1 and P2
這里模擬一個梯級潰決的洪水演進過程,其中h1=h2=100 m,L=1000 m,L1=L2=1000 m,α=12?,圖4 給出了4 個典型時刻的水位高程分布.0 s 時刻隨著上游壩體的潰決,上游水庫水體向下傾瀉而出,向下游水庫傳播,約10 s 左右上游潰壩洪水傳到下游水庫,如圖4(a)所示;形成一個穩定的向下游水庫傳播的水躍(激波),如圖4(b)所示;水躍傳播到下游壩體時候下游壩體潰決,從而形成洪水往下游傳播,如圖4(c)所示,部分上游水體和下游水庫合在一起宛如一個水塔為下游洪水提供充足的質量和動量來源;最終如圖4(d)所示,洪水繼續往下游傳播,直至消亡.

圖4 典型梯級潰決洪水演進過程水位分布Fig.4 Water elevation distribution at four typical times of cascade dam break flow
為研究梯級潰決下游河道的洪峰演化和增強現象,選取兩個參考測點P1和P2,其距離下游壩體分別為1000 m 和2000 m,考察其洪水演進時間歷程,結果如圖5 所示.為比較,同時給出了另外3 種潰決模式,分別是上游和下游壩體單獨潰決模式,以及將上游水體體積加載在下游水庫的單壩潰決模式(等體積潰決模式).圖5(a) 表明梯級潰決模式在P1測點的最高水位能達到23.56 m,下游單壩潰決模式在P1測點的最高洪峰為10.33 m,上游單壩潰決模式在P1測點的最高水位為7.50 m,而等體積潰決模式在P1測點的最高水位為17.19 m.綜合來看,梯級潰決模式在P1誘發的洪峰大大超過單板潰決洪峰,還超過上游下游單個潰決模式的最高洪峰之和,也超過等體積潰決模式的最高洪峰,在P2測點也有類似的規律.這論證了梯級潰決模式存在洪峰增強現象,上游水庫不僅對下游水庫進行了質量補充,還進行了大量的能量補充,但是這種質量和能量如何進行補充的,又是如何分布的,是梯級潰決洪水演進和洪峰增強的關鍵.
第2 節的計算分析表明梯級潰決誘發洪水洪峰超過單個壩體單獨潰決洪峰的總和,也超過將上游體積放到下游后等效的所謂等體積潰決模式的洪峰,梯級潰決模式中,下游水庫實際上既從上游水庫的潰壩洪水繼承了部分質量,而且繼承了部分能量,但是質量和能量是如何補充和分布的應當是梯級潰決洪峰增強的關鍵因素.本節重點分析梯級潰決的質量和能量轉化過程,并給出幾個關鍵機制的分析.

圖5 四種潰決模式下游測點P1(a)和P2(b)的洪水位時間演化,其中Cascade Mode 為梯級潰決模式,Downstream Mode 和Upstream Mode 分別為下游和上游單壩潰決模式,Single Equivalent Mode 表示等體積潰決模式Fig.5 Flood elevation evolution at detection points P1 (a)and P2 (b)for four dam break modes:cascade dam break mode,single downstream dam break mode,single upstream dam break mode,and single downstream dam break with full water mode
圖4 的梯級潰決洪水演化過程表明上游壩體潰決對下游水庫既增大了質量又增大了能量,圖6 給出了下游壩體潰決瞬間的水深分布和速度分布,下游水庫形成了一個基本穩定的梯形水塔,其自由面高程基本為常數(圖4(b)),由于底床為直線分布,因此水深為線性分布如圖6(a) 所示,其速度由后部至前部呈線性分布如圖6(b) 所示; 而后部實際拖了很長的一段水流,其速度由尾部向前部逐漸增大,前后兩部分通過一個運動水躍(激波)連接.因此初步的分析表明上游潰壩洪水對下游潰壩的效應體現在兩個方面,其一是增大了其質量,其二是增大了其動量,余下的水體質量和動量滯留在后部.以上機制表明可以考慮將梯級潰決誘發洪水過程簡化成一個帶動量的初始水體的單壩潰決情況處理.

圖6 梯級潰決下游壩體潰決瞬間的高程(a)和速度分布(b)Fig.6 Depth(a)and velocity(b)distribution of flood when downstream dam breaks
由于下游壩體臨近潰決瞬間的結構為前部集中分布的水體和后部拖了很長距離的一段水體,前部宛如一個水塔,將梯級潰決過程等效成帶動量水體的單壩潰決模式,稱之為“水塔”機制,因此建立一個所謂“水塔”模型,即帶動量水塔單壩潰決模型,如圖7 所示.水塔的質量和速度和原始模型協調給定,如圖8 所示,其中自由表面高程為常數,速度分布也考慮成常數.“水塔”模型的潰決洪水過程如圖9 所示,將初始時刻協調后,可以發現測點P1和測點P2的洪峰高度均小于原始梯級潰決過程的在兩個測點的洪峰高程.因此“水塔”機制僅僅反映了一部分梯級潰決洪峰增強機理.

圖7 梯級潰決洪水演進的水塔模型Fig.7 Water tower model of cascade dam break flow

圖8 梯級潰決洪水演進水塔模型的初始水深(a)和速度分布(b)Fig.8 Initial water depth(a)and velocity(b)of water tower model for cascade dam break flow

圖9 梯級潰決下游測點P1 (a)和P2 (b)洪水演進過程的水塔模型預測,虛線為原始模型計算結果,實線為水塔模型預測結果Fig.9 Prediction of flood level at detection positions P1 (a)and P2 (b)of cascade dam break flow by water tower model,dashed line is shallow water flow model prediction,and solid line represents water tower model prediction
以上水塔模型低估了梯級潰決洪峰的能量,進一步經過計算分析,對于本例,前部水塔的質量占總質量的61.2%,動量則占總動量的37.7%,后部水體的動量占據比例非常大,類似一股射流作用,對前部水塔有很強的質量、動量補充作用,以維持洪峰的高度,稱之為梯級潰決的“射流?水塔”機制.因此建立如圖10 所示的“射流?水塔”模型,射流為一股帶速度分布的等水深射流,而水塔則和水塔模型一致.射流的速度分布采取兩種模型,如圖11 所示,第一種為常數分布,第二種為線性分布,分別稱為射流水塔模型和分布式射流水塔模型.兩種模型模擬的洪水演進過程在測點P1和測點P2的洪峰演化如圖12 所示.結果表明射流水塔模型的預測結果明顯偏大,而分布式射流水塔模型的洪峰過程更加接近實際洪峰過程,但是依然存在一定的偏差.由此可知,分布式射流水塔模型相比于射流水塔模型、水塔模型以及單壩潰決洪峰疊加模式更能反映梯級潰決的洪水演進過程.

圖10 梯級潰決洪水演進過程的射流水塔模型Fig.10 Jet-water tower model for cascade dam break flow

圖11 梯級潰決射流水塔模型的水深(a)和速度分布(b),其中藍線虛線為射流水塔模型(JWT)的水深和速度分布,而紅線實線為分布式射流水塔模型(DJWT)的水深和速度分布,黑色點劃線是原始模型水深和速度分布Fig.11 Initial water depth(a)and velocity(b)of jet-water tower model for cascade dam break flow,blue dashed line,red solid line and black dash dot line represent water depth and velocity of jet-water tower model(JWT),distributed jet-water tower model(DJWT),and original shallow water model,respectively
梯級潰決洪水演進過程存在“射流?水塔”機制,一方面上游潰壩洪水增大下游庫容的質量和動量,另一方面殘余的水流構成一個高速的射流在下游潰壩洪水演進過程中不斷補充質量和動量,維持洪峰高度.第3 節總結了一個反映這種機制的分布式射流水塔模型,本節應用該模型預測兩組接近實際尺度的梯級壩潰決洪峰,考察不同的庫容比和壩間距的影響,總共設計7 組算例,其中坡降tan α=0.1%,計算區域為800 km,壩間距在100 km 量級,在距離下游壩處100 km 和200 km 處分別選取測點P1和測點P2,分別考察洪峰水位和流量的演進過程,同時選取原始淺水動力學模型的計算結果作為比較對象,相關結果整理后如表1 和表2 所示.

圖12 梯級潰決下游測點P1 (a)和P2 (b)洪水演進過程的兩類射流水塔模型預測,虛線為原始模型計算結果,藍線虛線為射流水塔模型(JWT)預測結果,紅線實線為分布式射流水塔模型(DJWT)預測結果,黑色點劃線為淺水動力學模擬結果Fig.12 Prediction of flood level at detection positions P1 (a)and P2 (b)of cascade dam break flow by jet-water tower model(JWT)in blue dashed line and distributed jet-water tower model(DJWT)in red solid line,and black dash dot line is the shallow water model solution

表1 實際尺度梯級壩潰決洪峰水位預測Table 1 Peak flood elevation prediction of real-scale cascade dam break flow

表2 實際尺度梯級壩潰決洪峰流量預測Table 2 Peak flood discharge prediction of real-scale cascade dam break flow
總體而言,分布式射流水塔模型反映了梯級潰決的洪水峰值,其預測的洪峰值和原始模型預測的洪峰值基本吻合.對于間距為100 km 的兩個壩體,庫容比在0.6 ~1.2 的范圍內,下游100 km 處洪水峰值的預測值誤差在6%以內,200 km 處洪峰的預測值誤差小于7%; 對于庫容比為1.0,間距在100 ~250 km范圍內的梯級潰決過程,下游200 km 以內,分布式射流水塔模型預測的洪峰誤差在15%以內.對于下游100 km 處測點,由表2 可知洪峰流量的預測值和淺水動力學模型計算結果吻合較好,下游測點200 km的P2測點,誤差隨著壩間距的增大而增大,這主要是由于潰決模式和洪水增強機理存在差異所致.
針對梯級潰壩洪水演進問題,本文建立了一維淺水動力學方程和數值模型,詳細研究了梯級潰的洪峰增強機制,并建立了一個梯級潰決洪水演進的單壩潰決等效模型.主要結論如下:
(1)梯級潰決的洪水過程存在增強機制,下游洪峰不僅超過單壩潰決洪峰,而且超過兩個單壩潰決的峰值之和,以及所謂等體積潰決模式的洪峰.
(2) 梯級潰決的洪峰增強由射流水塔機制控制,上游潰決誘發的洪水不僅增大下游水庫的質量和動量,形成“水塔”,而且在尾部殘留一個動量較大的射流,穩定和維持下游潰決后“水塔”的質量和動量.
(3)以梯級潰決洪水演進增強的射流?水塔機制為基礎,建立了一個等效的射流水塔單壩潰決等效模型,預測結果基本能夠反映梯級潰決下游洪峰演化過程,對實際百公里量級間距的梯級潰決洪峰過程預測良好.
梯級潰壩洪水過程十分復雜,潰決模式和增強機制方面依然還有待進一步研究,如在可沖刷河床情況下梯級潰決的洪峰增強機制,此外本文提出的射流水塔模型是否存在解析/半解析的參數化模型.本文所建立的梯級潰決洪水過程的單壩潰決等效模型為實際梯級潰決洪峰過程提供了理論上的可能,但是如何將本模型和實際河道如金沙江及其支流的梯級壩關鍵參數對應進行安全評估還有很多實際因素需要考慮,如地形、河道截面和能量耗散.