崔弘毅編 譯
(國家電力監管委員會大壩安全監察中心,浙江杭州310014)
4家研究機構參與了IMPACT項目泥沙運動專題的研究,分別為:天主教魯汶大學(UCL,比利時)、特倫托大學(UdT,意大利)、里斯本高級技術學院(IST,葡萄牙)和法國農業暨環境工程研究所(Cemagref,法國里昂)。
本文主要闡述研究團隊在該項研究涉及到的不同區域中的主要發現:近壩址和遠壩址潰壩洪水、對潰壩洪水的試驗模擬和數字模擬以及通過與試驗室數據和現場數據的對比進行基準試驗來校驗模型。
2.1.1 試驗室數據
試驗室試驗可使研究人員重點關注某些具體的過程,可以對其進行更詳細的分析研究。通過充分的測量,可得到精確的數據,數據用途有二:一是可以幫助更好地理解正在發生的現象的物理原理,二是為數字模型提供校驗基準。在天主教魯汶大學(潰壩洪水)和特倫托大學(均一泥石流)的試驗室中進行了一系列試驗,考慮了兩類性質:
(1)在近壩場,伴隨著潰壩波的演進,發生了快速而強烈的沖蝕。洪流顯示了強烈的自由表面的特性:波破碎發生在中間位置(在壩址附近),接近垂直的水和殘渣墻在波浪前沿掀起了泥沙底床,造成瞬時強烈的泥石流(見圖1)。但是,在潰壩波前沿,泥石流卻意外地和均勻波相似。因此,該項研究的第一項研究內容是在更深入地調查研究潰壩流條件下的性態之前,研究均一條件下的泥石流的特性。
(2)在遠壩場,泥沙運動依然很劇烈,但泥沙的動態作用減弱。因泥沙分散、兩岸沖蝕和殘渣沉積,山谷地貌發生了戲劇性的變化,見圖2。第三項研究內容就是研究遠壩場的性態。

圖1 近壩場的泥石流Fig.1 Near-field geomorphic flow(UCL)

圖2 斷續的壩段潰決導致的兩岸沖蝕Fig.2 Bank erosion resulting from intermittent block failure
2.1.2 現場數據
泥沙運動主題研究所利用的現場案例是Lake Ha!Ha!潰決事件。Lake Ha!Ha!位于加拿大魁北克沙格奈河的一條支流上,于1996年潰決(Brooks&Lawrence,1999年)。在該潰決事件中,發生了大量的泥沙運動,山谷地形地貌發生了巨大的改變:天然河道變位偏移,兩岸沖蝕引起河道中等或大幅變寬,受不可沖蝕的基巖的影響,洪流進一步偏離發生變道(洪流被撕裂)。該項研究的目的是收集與高度瞬變流相符,包括強烈泥沙轉移的現場數據。對相關文獻進行核查,以便發現最合適的案例,并有足夠可用的數據。
2.2.1 試驗室數據
可用的數據系列有:
(1)均一泥石流,用以調查研究作用力和速度分布(UdT)。
①均一材料(PVC顆粒);
②級配材料(PVC顆粒和沙)。
(2)潰壩洪水,用以調查研究近壩場的影響:沖刷和泥石流流阻的形成(UCL)。
①大壩上下游河床高程一樣(Spinewine&Zech,2002b);
②河床高程有一個抬升,水庫內的高程更高(Spinewine&Zech,2003)。
(3)梯形河谷中的潰壩洪流,用以調查研究遠壩場的影響。兩岸沖蝕和河道拓寬(UCL)(Le Grelle等,2003,2004)。
①均一材料(沙);
②級配材料(沙和粗礫石)。
2.2.2 現場數據
關于Lake Ha!Ha!潰壩事件,收集到了洪水前后的大量數據。數據處理工作由天主教魯汶大學(UCL)、臺灣大學、魁北克大學和加拿大地質調查局聯合進行,其中,加拿大地質調查局是數據所有者。整套數據包含完整的洪水前后的數字地面模型(30 km范圍)以及重新制作的堤壩潰決的出流過程線。
2.3.1 試驗室數據
試驗取得了高質量的數據系列,用于校驗數字模型。
2.3.2 現場數據
現場數據可幫助識別和理解由潰壩洪水引起的地形地貌變化的關鍵特征,以便涵括未來數字模型的必要方面。另外,通過現有數字模型得出的結果和觀察結果之間的比較,可評價現有數字模型的性態。
研究的目的是對調查過程進行數學描述。在將現有的洪流描述延伸去解釋試驗認定的過程之前,需對文獻進行復核。
3.2.1 均一泥石流
由于快速顆粒流和氣體之間存在某些物理相似性,將動力理論用于顆粒材料需耗費大量工作。所有模型都假設微粒間的相互作用是由瞬時碰撞引起的,這就意味著只需考慮二元或雙微粒碰撞。
Jenkins和Hanes(1998年)曾將動力理論用于薄層水流,其顆粒由碰撞互相作用支持,而不是由紊流的速度波動支持。顆粒壓力的本構關系被視為與Chapman和Cowling于1970年提出的稠密分子氣體的準彈性近似,其描述了顆粒間碰撞率集中度的變化。
假設顆粒浮重整個被顆粒碰撞接觸所支撐,則從試驗中可能推導出顆粒壓力σs和剪切應力τs。該項目框架中關于明確本構關系的最主要進步為解釋了附加質量效應,即通過公式(1)替換泥沙密度ρs:

式中,Cs為顆粒集中度。
3.2.2 近壩場潰壩洪水
3.2.2.1 2D-V水平集模型
考慮到潰壩洪水第一階段速度的垂直分量不可忽略,最初是想開發一個2D-V模型,可以表示出大壩潰決后最初的瞬間垂直中正面發生的情況。最合適的模型貌似是水平集方法。該方法基于的假設為:洪流被細分為性態近似相同的層次,由尖銳的界面分開。多種介質(空氣、水、泥沙)中行進的界面對應更高維函數Φ的零點水平集,被定義為到界面的符號距離(Sethian,1999年)。就渦度和流函數Ψ而言,可建立Navier-Stokes方程,從方程可得到速度場。根據速度場,則可得水平集方程,得到有用的定性結果,見圖3。

圖3 水平集方法原理Fig.3 Level-set method principle
3.2.2.2 雙層淺水1D模型
最初的進展由Capart于2000年提出。洪流由三個層次代表:(1)上層水,由清水組成,深為hw;(2)運動的泥沙層,厚度為hs;(3)固定的基巖層,基巖高程zb,作為上限。在最初的模型中(Capart,2000年),假設泥沙的集中度不變(Cs=Cb),且水和泥沙(hs)混合體的上部都與清水層一樣,以同樣均一的速度運動(us=uw)。根據這些假設,剪應力應該在垂線上是連續的。模型得出了解析解(Fraccarollo和Capart,2002年),盡管很靈巧,但不能用于實際的幾何構型。
關于模型的最主要進展之一(Spinewine,2003;Spinewine和Zech,2002a)是賦予集中度和三個層次間速度以新的自由度(Cs≠ Cb,us≠ uw),見圖4。
該描述中得出的方程由二階Godunov有限容積法解答,其中流量用LHLL Riemann解來計算(Fraccarollo等,2003年)。

圖4 近壩場洪水流數學描述的假設Fig.4Assumptionformathematicaldescriptionofnear-fieldflow
3.2.3 遠壩場潰壩洪水
3.2.3.1 二維模型
首先,開發了一個2D可拓模型來表示近壩場,包括岸坡沖蝕機理。這里主要總結其方法,詳細信息參見Spinewine等(2002年)和Capart&Young(2002年)。關鍵點在于認定分離的水流和液狀泥漿層各自獨立流動,控制方程就完全能處理岸坡料坍塌陷入水流的滑塌事件。一旦潰決發生,則潰決后的洪水就像其它水和泥沙的運動模式。
因此,需要液化標準來明確什么時候及什么地點,岸坡從固態轉換成液態介質。對此,假定以下基本機理:當局部邊坡超過臨界角度αc時,則發生壩段潰決;延伸的潰決表面被確定為錐形,中心在潰決處,以剩余角αr<αc向外傾斜。最后,假設椎體以上的泥沙材料在潰決時立即液化。為解釋觀察到的淹沒區域和未淹沒區域的區別,如圖5所示,定義了4個不同的安息角:αc,subm和φr,subm用于描述淹沒區域,αc,em和αr,em用于描述未淹沒區域。

圖5 2D地壓潰決算子穩定示意圖Fig.5 Stability diagram for the 2D geostatic failure operator
3.2.3.2 全部岸坡潰決的一維模型
第二個選來耦合以上岸坡沖蝕機理的模型是一個一維方案,它包含了一個流體動力有限容積法和一個單獨的泥沙運動路線。有限容積法,其開發目的是為應對復雜地形(Soares-Fraz?o和Zech,2002年),它解出了流體動力淺水方程,因在一個計算時間步長上的縱向泥沙運動(推移質)導致斷面地形發生的部分改變可以從泥沙相的Exner連續性方程得出。除單元格上下游面的泥沙流量外,岸坡潰決導致的側向泥沙流入量設為體積Vs,其在計算時間步長末端的斷面上將重新分布。
水位升高Δh淹沒岸坡,引起如圖6所示的棱柱形部分材料失穩,最終引發潰決。試驗中,岸坡原始角α小于水面上穩定角αs,em,但大于水面下的穩定角αs,subm。因此,當水位上升時,岸坡立即變得不穩定,當達到潰決角度時,則發生潰決。對應于淹沒工況和出露工況,潰決角度分別為αf,subm和αf,em(在實際情況中,潰決角度αf比穩定角度αs小一點)。斷面上的沖蝕體積Vs也重新分布。

圖6 因淹沒岸坡引起的岸坡潰決Fig.6 Bank failure triggered by the submergence of the bank
如圖7所示,沖蝕材料沉積在河道中,淹沒部分,對應水下的安息角,泥沙堆積角度為αr,subm,出露部分則穩定在角度αr,em(堆積過程發生后的水上安息角)。所有這些安息角都由試驗中所用材料決定,測量都在靜態和動態試驗中完成。

圖7 岸坡沖蝕材料堆積示意圖Fig.7 Deposition of the material eroded from the banks
最終,數字一維模型關鍵在于以去耦合的方式解決過程中的三個不同關鍵步驟:(1)水流的流體動力路線;(2)縱向泥沙運動和由此導致的沖蝕和沉積;(3)岸坡潰決和由此導致的斷面形態的變化。
3.2.3.3 局部岸坡潰決的一維模型
以上方案在理想工況中適應性良好,其中,理想工況是指斷面明確,例如是長方形或梯形,在其棱角方面,也只有有限的幾個凸峰。但對天然河道來說,其斷面復雜得多,不能用這么簡單的方式來描述。在描述天然河道時,推薦使用Schmautz和Aufleger于2002年提出的方法(本文不介紹)。
斷面剖面被切分成小塊,從山谷邊開始,對每一小塊的穩定性進行復核(如圖8所示)。若岸坡斷面AB局部比臨界值(穩定角αs)更傾斜,則岸坡部分會發生轉動,直至達到A’B’位置,對應安息角αr(淹沒或出露)。結果,這一新位置會加劇周邊區域的穩定性(例如BC已移至B’C位置)。需多次觀察整個剖面,直到所有區域都穩定。

圖8 局部岸坡潰決模型原理Fig.8 Principle of the local bank-failure model
對于縱向的泥沙運動,無論對全部或是局部岸坡潰決模型,以下規則都適用:(1)若發生沖蝕,根據 Meyer-Peter.Müller公式,假設運動與(τb-τb,c)3/2的局部值成比例,其中τb和τb,c分別為河床的實際剪切應力和臨界剪切應力;(2)若發生沉積,則假設泥沙沿河床均勻沉積(不是水平地),這在后面被定義為斷面因素,其坡度小于水下的安息角。
3.3.1 2D-V模型
事實很快表明,在現實案例中使用這么復雜的模型幾乎是不可能的。提出整個現象沿特定垂直面的表達的約束條件是不現實的,因為對大多數真實河谷來說,其在大壩附近變得非常窄,而在橫河向,水深變化極大。因此,2D-V方法只能作為通往完全3D方法途中的一個有趣步驟,而完全3D方法是建模者的遠期目標。
3.3.2 1D和2D-H模型
描述潰壩洪水條件下泥沙運動的方程組還未完全建立起來。因此,現階段不可能開發出用于商業用途的工具包來解決隨潰壩事件而產生的與嚴重的泥沙轉移有關的各種問題?,F有模型,包括其延伸的各種可能(如之前提到過的2D-H模型),在實際應用中都太慢太缺乏效率。
3.4.1 2D-V模型
2D-V方法只能作為通往完全3D方法途中的一個有趣步驟,而完全3D方法才是建模者以后的重要攻克目標。
3.4.2 1D和2D-H模型
所有比較的模型中,至少在平坦河床基準試驗中,都存在一個普遍的缺點,那就是它們推進前沿太快,這是因為沒考慮垂直速度分量,因此造成了錯失泥沙流動的初期階段的事實。另一明確的結論是,要重現飽和殘積物前沿的沖蝕性態很困難,特別是之后發生部分再堆積的沖蝕。
但是,與數年前可獲得的結果比較,這樣的建模過程是令人驚嘆的。取得該進展的部分原因是基于新的數字成像測量技術,它使對液相和泥沙運動層的速度場進行實時觀測成為可能。
二層一維方法及其在2D-H模型中的延伸看起來都是有趣且有前途的方法。未來的努力方向是對岸坡潰決算子的描述,包括岸坡潰決的引發,加入對土力學的考慮就可以實現。從數字模擬的角度看,還應努力使計算在合理的時間內完成。
共進行了3個基準階段的測試,不同的建模者可以在IMPACT項目中進行盲模試驗測試其模型。通過比較盲模結果及分析不同數學描述后的模擬結果,可以知曉每種建模方法的優點和局限性。對每條基準,要求建模者提供其數值模型的描述,從而可以深入分析結果。
4.2.1 流過原本平坦河床的潰壩洪水
Spinewine和Zech(2002b)給出了基準描述,有4家研究機構得出了結果,即參與IMPACT項目泥沙運動專題的小組:法國農業與環境工程研究所(Cemagref,法國里昂)、特倫托大學(UdT,意大利)、里斯本高級技術學院(IST,葡萄牙)和天主教魯汶大學(UCL,比利時)。圖9顯示了典型結果,其中可見計算所得的河床高程、運動的泥沙層高程和洪水高程與試驗數據的對比。
4.2.2 流過原本階梯狀河床的潰壩洪水

圖9 潰壩波流過原始平坦易蝕性河床的基準試驗結果和數值結果Fig.9 Experimental and numerical results from the benchmark on dam-break wave over an initially flat erodible bed
Spinewine和Zech(2003)給出了基準描述,有4家研究機構得出了結果,即參與IMPACT項目泥沙運動專題的小組:法國農業與環境工程研究所(Cemagref,法國里昂)、特倫托大學(UdT,意大利)、里斯本高級技術學院(IST,葡萄牙)和天主教魯汶大學(UCL,比利時)。階梯狀河床基準下的各模型結果與試驗數據的比較見圖10,圖中顯示了在給定時間的不同高程。


圖10 潰壩波流過原始階梯狀河床的基準試驗結果和數值結果Fig.10 Experimental and numerical results from the benchmark on dam-break wave over an initially stepped bed
4.2.3 原本棱柱形河谷中由潰壩洪水引發的岸坡沖蝕
Grelle等(2003年)給出了基準描述圖11為其基準,4家研究機構得出了結果(圖12),即參與IMPACT項目泥沙運動專題的小組:法國農業與環境工程研究所(Cemagref,法國里昂)、特倫托大學(UdT,意大利)、里斯本高級技術學院(IST,葡萄牙)和天主教魯汶大學(UCL,比利時)。

圖11 試驗測量Fig.11 Experimental measurement
4.3.1 流過原本平坦河床和原本階梯狀河床的潰壩洪水
關于前沿波速,特倫托大學(UdT)的結果利用了校準過程,包括用這些波速數據作為校準參數。與之形成對比的是,運動泥沙層估值偏低,原因是假定運動泥沙層的集中度與河床材料一樣,但天主教魯汶大學(UCL)和里斯本高級技術學院(IST)的模型情況并非如此,事實是,運動泥沙層的集中度在降低,以使顆粒可以運動。由前沿移動造成的沖蝕僅在天主教魯汶大學(UCL)和法國農業與環境工程研究所(Cemagref)模型中出現,盡管法國農業與環境工程研究所的簡單模型不能提供任何關于運動泥沙層的結果,但仍得出了激波過后的關于水體表面的有用預估值。在這一點上,對沖蝕和沉積的不對稱處理可解釋為天主教魯汶大學模型的成功之處。
4.3.2 原本棱柱形河谷中由潰壩洪水引發的岸坡沖蝕
法國農業與環境工程研究所(Cemagref)模型中僅有河床運動,不考慮斷面的典型增寬,而其顯然是河岸穩定標準的需要。在特倫托大學(UdT)模型中,河岸沒有重塑形,而是所有河岸材料都被轉移到了底部,明顯導致河床高程的估值偏高。在描述岸坡失事時,里斯本高級技術學院(IST)的模型奇怪地失效了,而河床加深會引發該機理的發生。天主教魯汶大學(UCL)所使用的模型中,對因岸坡垮塌而導致的材料沉積,明顯利用了安息角的定義。
必須注意的是,試驗中,原本岸坡角度比臨界角大,強調了該現象。對于更平緩的岸坡,形態效應不那么重要。
大壩或堤壩失事產生的洪水可引起多種形式的、嚴重的泥沙運動。某些案例中,洪水帶走的材料體積可與潰壩釋放的洪水達到同一數量級(達幾百萬立方米),因此,泥沙運動的風險是巨大的。
現有模型還不能對大多數不確定性相關因子的識別做出詳細而敏感的分析。在模擬洪水事件中,重點在強調對泥沙運動核算的重要性上。
強烈的沖蝕和泥沙轉移使河谷地形發生了巨大而快速的變化。反過來,地形變化也對波浪性態有強烈的影響,因此,到達時間和最大水位高程也受到影響,而這是風險評估和組織預警的兩個關鍵參數。這就意味著影響泥沙運動預測的不確定性因素可能最終影響整個預測過程。為證實這些影響,對波浪在固定、光滑的河床和波浪在有移動的泥沙上的傳播性態進行比較(見圖13)。

圖12 10 s后的斷面圖Fig.12 Cross section after 10 s

圖13 固定河床(點線)和移動河床(實線)上的潰壩波在t=1.5 s時的比較Fig.13 Comparison between dam-break wave on fixed(dotted lines)and mobile bed(solid lines)at t=1.5 s
從圖13中可見,泥沙的移動轉移了部分有效勢能,因此波的前進速度大大放緩,這對預警和下游人口的應急規劃來說是一利好。但波前沿后面的水深大幅增加,至少在近壩場有大幅增加,增大了危險區域和潰壩受災區人民的風險。
此類現象中的大多數過程是不確定的。獲取建模所需的數據通常也非常困難。河床底部的組成材料不均一,厚度也不明確,材料會隨時間嚴重離析,特別是細粒材料。大壩下游河谷的河床材料也不均一,其中包含的土和巖石以不可預知的組合方式排列著。對河床進行測量非常困難、成本太高而且乏味。
這就意味著,現在還沒有標準的做法可推薦用于考慮泥沙對潰壩波的影響。同樣意味著終端用戶必須清楚知道其模型在這方面存在著局限性。
模型是在對潰壩的理想化狀態基礎上進行模擬的。問題被表示成為一個豎直平面,且大壩被假設為潰決后立即消失,沒考慮橫向效應。在近壩場泥沙運動中,僅考慮了河谷河床材料,忽略了潰壩本身釋放出來的材料。模型發展現階段,移動河床建模還沒與潰壩建模耦合起來。模型對理想狀況來說非常有應用前景,但要用于真實案例中的情況,還有很大差距。
對遠壩場,關鍵是要重現河谷的演變過程,即因上游固體材料輸移及壩堤潰決引起的持續沖蝕和泥沙沉積過程。對此,可以模擬部分地形變化,尤其是局部的地形變化,但僅限于幾公里范圍內。因存在許多隨機現象,因此潰決發生后的一系列事件變得不可預測,從而形成類似“不確定樹”現象,很難掌控。
經過數據整理,真實的案例研究包括:(1)數據解釋;(2)將合理的假設應用到數字模型中:確定網格,一維還是二維?邊界條件情況如何,有時候并不能獲取所有邊界條件的物理數據;(3)解決問題(運行計算);(4)對結果進行評判性分析,評估價值和結果的“真實性”。
1996年7月,Lake Ha!Ha!水庫副壩潰決,產生了強烈的潰壩波,導致下游30 km長河道的地形發生了重大改變,在Ha!Ha!灣匯入Saguenay河(見圖14)。

圖14 Lake Ha!Ha!示意圖(Brooks,2003年)Fig.14 Lake Ha!Ha!(Brooks,2003)
事實上,可以觀察到由潰壩災難引起的巨大地形變化的所有典型特征:大范圍的沉積,以致于河流河道發生變化(見圖15(a)),河道大幅加寬(見圖15(b)),有時因河床基巖底坎和河岸的出現而發生阻塞(見圖15(c)),改變了河床的地形地貌,改變了河流的路線等。


圖15 Lake Ha!Ha!潰壩洪水波洗劫后的典型地形變化(Brooks,2003年)Fig.15 Typical morphological evolutions after the Lake Ha!Ha!dam-break wave(Brooks 2003)
加拿大地質調查局、魁北克大學、臺灣大學和天主教魯汶大學(Capart等,2003年)做了巨大的努力來對獲得的數據進行判讀,以獲得有用的數據系列,這或許是現實案例中可用來校驗模型的最好的數據系列之一。
圖16為發生大規模地形變化的周邊區域的河床剖面變化圖。原河床剖面主要為不可沖蝕的巖石區域,對應圖16中綠色的斜線段。因在河岸線的壓低處發生了漫頂,河流河道發生改變。河床巖石區域被繞開,引起嚴重的溯源沖刷。
雖然法國農業與環境工程研究所(Cemagref)的一維模型沒使用關于泥沙運動的復雜的描述(普通泥沙轉移的Exner方程),但結果得到的河床剖面變化方向正確,盡管存在數值不穩定性。其中存在一些差距的原因與巖石位置有關,且與這一實際情況有關:計算在2天后就停止了,而此時沖蝕過程還沒有完成。
圖17為與圖16相對應的相同河段的水位剖面圖。從圖中可見河床移動和地形變化的作用。水位高程數據來自兩種不同的計算,即法國農業與環境工程研究所(Cemagref)的模型,計算工況為移動河床和與之形成對比的固定河床法(El Kadi和Paquier,2004年)。在有些地方,因河床的移動性,導致水位上升或下降達5 m(例如圖17中21~22 km處)。

圖16 Ha!Ha!河潰壩前后的河床剖面,對比數字模型Fig.16 Bed profile of the Ha!Ha!River before and after the dam break.Comparison with numerical models

圖17 Ha!Ha!河接近洪峰時(7月20日19∶30)的水位高程剖面圖(法國農業與環境工程研究所(Cemagref)的模型:固定河床和移動河床間的對比(2種方法))Fig.17 Water profile of the Ha!Ha!River near the flood peak(July 20,19:30).Cemagref numerical model:comparison between fixed and mobile bed(2 approaches)
臺灣大學使用的擴散移流模型是二維的、各向異性的,但其基本假設卻很簡單,該模型對Lake Ha!Ha!這一典型案例的應用來說非常有效。根據當地具體坡降,用沿坡的擴散來表示地形變化。不要求明確的流體力學計算:在凹陷處,地下水面是水平的,而其它地方的深度為0。

此方法得出的結果令人印象深刻,至少在有些河段給人留下了深刻印象,比如上游河段(見圖18)。二維模型貌似可以捕捉到剛好潰決的壩體下游處的沖蝕區域和主要沉積區域。一維模型也可以預測某些沉積特征,并且,在一定程度上,還能得到沖刷區域的時間精確的結果。對有些河段的加寬也有定性的表示。
在Lake Ha!Ha!這樣的復雜案例中,僅有最簡單的方法能成功得出某些結果,復雜的模型也無法應付海量的必需的數據,這是非常有意思的現象。
[1]Final Technical Report of Investigation of Extreme Flood Processes And Uncertainty[R].2005.
[2]H.CAPART,B.SPINEWINE,D.L.YOUNG and et al.The 1996 Lake Ha!Ha!breakout flood,Quebec:test data for geomorphic flood routing methods[R].Revision Submitted to the Journal of Hydraulic Research.