新疆射月溝水庫潰壩過程數值模擬

馬 黎，鐘啟明，楊 蒙，吳 迪，李大成，梅勝堯

(1.中國電建集團貴陽勘測設計研究院有限公司，貴州 貴陽 550081； 2.南京水利科學研究院，江蘇 南京210029； 3.水利部水庫大壩安全重點實驗室，江蘇 南京 210029)

近年來，受全球氣候變化影響，我國極端暴雨事件呈增加趨勢[1]。2000 年1 月至2021 年12 月，我國共發生99 起潰壩事故，多起潰壩事故發生在干旱和半干旱地區，其中新疆和內蒙古自治區分別發生9 起和6 起由極端暴雨導致的潰壩事故[2-3]。如2018 年7 月31 日，位于新疆自治區哈密市的射月溝水庫大壩因壩頂溢流發生潰決[4]。合理預測潰壩過程對潰壩洪水致災后果的評價具有重要意義，采用數值分析方法模擬潰口流量過程和潰口尺寸演化過程是預測潰壩過程的有效手段。目前常用的潰壩過程數學模型可分為經驗模型和基于潰決機理的模型[5-7]。經驗模型一般通過考慮壩體形態、水庫特征、壩料特性等信息，基于回歸分析或機器學習方法建立經驗計算式，獲取潰壩時的潰口峰值流量、最終潰口尺寸和潰壩歷時等潰壩參數[5]。但經驗模型無法獲取潰壩時的潰口流量過程和潰口發展過程?；谕潦瘔螡Q機理的潰壩過程數學模型[8-14]，一般采用寬頂堰流或孔流公式模擬漫頂或滲透破壞潰壩時的潰口流量，并利用各種沖蝕公式來描述壩料的輸移，一些模型還可以考慮潰口邊坡的失穩[13-14]，通過按時間步長迭代的數值計算方法獲取各時間步長的潰口流量和潰口尺寸。

本文通過對射月溝水庫潰壩案例的調研，總結潰壩水流的沖蝕特性和潰口發展過程。在此基礎上建立基于壩體潰決機理的簡化漫頂潰壩過程數學模型，反演射月溝水庫的潰壩過程，并與國外典型的NWS BREACH 模型的計算結果進行對比分析。

1 射月溝水庫潰壩事件調研

射月溝流域地處哈爾里克山南坡，由6 條支流匯流形成射月溝水庫（圖1）[15]。射月溝水庫集水面積約408 km2，徑流補給主要由季節性冰雪融水、夏季降雨和泉水組成，年平均降雨強度約為150 mm，年平均徑流量為1.58×107m3。

圖1 射月溝流域水文概況Fig.1 Hydrological survey of Sheyuegou Basin

射月溝水庫大壩距離二宮村頭宮隊2 km，距離哈密市115 km。大壩為瀝青混凝土心墻砂礫石壩，3 層瀝青混凝土心墻厚度均為0.4 m，壩頂高程為1 496.65 m，防浪墻頂高程為1 497.65 m；大壩最大壩高41.15 m，壩頂長403 m、寬6 m，上、下游坡比（垂直/水平）分別為1∶2.25 和1∶2.00[16]，射月溝水庫大壩典型斷面見圖2。

圖2 射月溝水庫大壩典型斷面示意Fig.2 Schematic diagram of the typical section of Sheyuegou reservoir dam

1.1 2018 年“7·31”洪水產匯流過程

2018 年7 月31 日01:00 至13:00，射月溝流域發生特大暴雨，小堡、射月溝、秦城3 個氣象站12 h的降水強度分別為116、105 和79 mm（3 個氣象站位置見圖1），均打破了當地的極值記錄[17]。2018 年7 月31 日22:00，新疆自治區水文局基于坡面面積法，通過現場洪水痕跡調查，得出射月溝水庫入庫洪峰流量達到1 848 m3/s，遠超設計標準，導致水庫發生漫頂潰壩。趙懸濤等[15]通過模擬獲取了射月溝水庫2018 年7 月31 日的降雨-入庫流量過程（圖3），入庫峰值流量相對誤差僅為3.6%，在本文的潰壩過程模擬中采用此入庫洪水流量過程線。

圖3 射月溝水庫2018 年7 月31 日降雨-入庫流量過程線Fig.3 The rainfall-inflow hydrograph of Sheyuegou reservoir on July 31, 2018

1.2 射月溝水庫漫頂潰壩過程

2018 年7 月31 日8:00 至11:40，一名目擊者用手機錄制了射月溝水庫潰壩過程的多段視頻。雖然視頻是臨時拍攝的，但獲取的信息對研究人員直觀了解壩體沖蝕特性和潰口形態演變過程具有重要意義。

根據現場視頻，可將射月溝水庫潰壩過程分為4 個階段：（1）初始沖坑出現。水流越過壩頂后表現為全斷面溢流，由于壩頂處水流流速較小，頂部未出現明顯沖蝕；受重力作用，水流在下游坡行進時流速逐漸增大，導致下游坡某一位置的土體發生沖蝕，形成初期沖坑（圖4(a)）。（2）溯源沖蝕。在漫頂水流的作用下，沖坑以上下游坡逐漸變陡，直至接近垂直；隨后發生溯源沖蝕，在此過程中下游坡面幾乎保持垂直，沖蝕不斷向上游發展（圖4(b)～(c)）。（3）初始潰口出現。當溯源沖蝕發展至上游水庫后，由于壩頂高程突然下降，初始潰口出現，漫頂水流匯聚到潰口處，其他部位的漫頂沖蝕結束（圖4(d)）。（4）潰口擴展。隨著漫溢水頭的逐漸增大，潰口在橫斷面和縱斷面的沖蝕逐漸加快；由于潰口形態的快速變化，潰口邊坡發生失穩；當潰口處水流的剪應力小于壩料的臨界啟動剪應力時，沖蝕過程停止（圖4(e)～(f)）。

圖4 射月溝水庫潰壩過程照片Fig.4 Photos of the Sheyuegou dam breach process

根據潰壩視頻和現場調查，大壩內部的瀝青混凝土心墻無法抑制潰口的發展，其原因主要有兩個方面：一是瀝青混凝土心墻均位于壩體上游部位（圖2），主要起防滲作用，對初期潰口的發展影響微弱；二是溯源沖蝕發展至瀝青混凝土心墻位置時，潰口處的漫溢水頭已達到較大值，僅有40 cm 厚度的瀝青混凝土心墻無法阻止潰口的發展。因此在數值模擬中不考慮心墻。

2 基于潰決機理的簡化漫頂潰壩過程數學模型

基于射月溝水庫漫頂潰壩過程的視頻和現場調查資料，建立基于物理機制的簡化潰壩數學模型。由于壩體潰決時間相對較短，滲流的影響可以忽略不計。模型主要包括4 個部分：潰口流量計算、初始沖坑位置確定、縱斷面潰口發展和橫斷面潰口發展。

2.1 潰口流量計算

水量平衡是潰壩模擬中必須遵循的基本原理，利用水量平衡關系描述水庫水位變化關系：

式中：As為與庫水位相對應的庫面面積；zs為庫水位；t為時間；Qin為入庫流量；Qb為潰口流量；Qs為溢洪道泄流量；QT為閘門泄流量。

潰口流量可采用寬頂堰流量方程[18]表示：

式中：c1和c2分別為流量系數，其中c1= 1.7 m0.5/s，c2= 1.3 m0.5/s[18]；b為潰口底寬；H=zs－zb，其中zb為潰口底高程；m為潰口邊坡坡比（垂直/水平）的倒數；ksm為尾水修正系數[18-19]。

2.2 初始沖坑位置確定

根據相關文獻[20-21]，土石壩潰壩時初始沖坑的位置計算如下：

式中：ln為從下游坡頂順坡至初始沖坑的距離；β0為大壩下游坡角；dn為下游坡法向水深；Frn為初始沖坑處的弗勞德數；Un為水流在初始沖坑處及越過沖坑后的法向平均流速，Un=C(Rnsinβ0)1/2=C(Hsinβ0)1/2，其中C為謝才系數，Rn為水力半徑；g為重力加速度；Bw為壩頂處漫頂水流的頂部寬度；Bn為下游壩坡處漫頂水流的頂部寬度；當漫頂水流越過初始沖坑后，dn=Qb/(UnBn)。由于潰壩第一階段為全斷面溢流，因此假設Bw=Bn。

2.3 縱斷面潰口發展

根據射月溝潰壩案例的視頻記錄，在數學模型中可將溯源沖蝕過程分為兩個階段：（1）初始沖坑出現后，水流繼續掏刷沖坑，下游壩坡坡角逐漸變大至基本垂直，形成陡坎；（2）陡坎保持直立狀向上游發生運移，直至進入水庫。

潰口的底床沖蝕率采用基于剪應力的沖蝕公式[22]表示：

式中：E為底床沖蝕率；kd為壩料沖蝕系數，可通過試驗測量[23]或經驗公式計算[24]，即kd= 10(ρw/ρd)exp，其中ρw為水的密度取1 000 kg/m3，ρd為壩料的干密度，cr為壩料的黏粒含量；τb為潰壩水流的剪應力，可通過曼寧公式計算[10]；τc為壩料顆粒的臨界起動剪應力，可由希爾茲曲線[25]獲取。

對于下游壩坡變陡階段，坡角增量計算如下：

式中：dβ/dt為下游坡角變化率；E1為沖坑處的底床沖蝕率；E0為壩頂處的底床沖蝕率。

對于陡坎運移階段，利用時間平均的沖蝕率公式來計算溯源沖蝕率[26]：

式中：dx/dt為陡坎運移速率；CT為陡坎運移系數；q為潰口單寬流量；He為陡坎高度。

隨著陡坎的運移，陡坎上游側的壩體在縱斷面方向逐漸變薄，在水動力作用下可能沿水平面破壞，可通過比較壩體在水平破壞面上的驅動力和抵抗力來評價陡坎上游側壩體的穩定性。驅動力由上游水壓力和陡坎底部的水流剪應力組成，而抵抗力則包括水平破壞面和陡坎兩側的黏聚力和摩擦力及尾水引起的水壓力。當驅動力大于抵抗力時，陡坎上游側壩體會發生失穩破壞，并假定失穩的壩體物質在瞬間被漫頂水流沖走，至此溯源沖蝕發展至水庫。

2.4 橫斷面潰口發展

當溯源沖蝕發展至水庫后，漫頂水頭突然增大，壩體橫斷面出現初始潰口，根據射月溝潰壩視頻假設初始潰口的形狀為矩形。隨著水流的不斷侵蝕，初始潰口迅速加深加寬，當潰口深度達到臨界高度時，潰口邊坡沿平面發生失穩破壞，初始潰口由矩形變為倒梯形，此后漫頂水流繼續沖蝕潰口直至潰壩結束。

假設漫頂水流對潰口底部和側壁具有同等的沖蝕能力，則潰口頂寬增量為：

式中：ΔBt為潰口頂寬增量；nloc為潰口位置參數（nloc= 1 代表單側沖蝕，nloc= 2 代表雙側沖蝕）；α為潰口邊坡坡角。

潰口底寬增量表示為：

式中：ΔBb為潰口底寬增量。

利用極限平衡法評價潰口邊坡的穩定性，當失穩塊體自重引起的驅動力大于失穩塊體自重和土體黏聚力引起的抵抗力時，潰口邊坡發生失穩破壞。視頻記錄和其他模型試驗表明，邊坡失穩時的破壞面近乎平面，并假定失穩土體被潰壩水流瞬間沖走。潰口邊坡失穩的驅動力Fd和抵抗力Fr分別為：

式中：W為失穩塊體重力；θ為失穩后的潰口邊坡坡角；γs為邊坡土體重度；Hs為潰口邊坡的臨界高度；φ為土體內摩擦角；C為土體黏聚力。

采用按時間步長迭代的數值計算方法模擬潰壩過程中的水土耦合效應，每個時間步長可輸出潰口流量、潰口尺寸、庫水位等潰壩參數。

3 射月溝水庫潰壩案例數值模擬

采用上述建立的數學模型對射月溝水庫潰壩過程進行反演分析，通過視頻記錄和現場實測數據，驗證模型的合理性，并對關鍵輸入參數進行敏感性分析。

3.1 輸入參數

從拍攝的視頻中可以看到，射月溝水庫大壩于7 月31 日9:10 開始漫頂溢流，并在當日11:00 左右潰壩結束。模擬持續時間為3 h，時間步長為1 s。根據視頻記錄，將潰壩計算時的初始水位設置為1 497.75 m，超過防浪墻頂高程0.1 m，入庫流量過程如圖3 所示。壩料沖蝕系數（kd）和陡坎運移系數（CT）是潰壩模擬的重要指標。根據設計報告，ρd=2 240 kg/m3，cr= 0.04，由前文經驗公式可得，kd=3.0 cm3/N/s。另外，基于潰壩過程視頻確定的陡坎移動時間，推導出CT= 0.001 /m1/6/s2/3。表1 給出了射月溝水庫和壩體的物理力學指標[4,16]。

表1 模型輸入參數Tab.1 Input parameters of the numerical model

3.2 計算結果分析

由于初始漫頂溢流的水位為0.1 m，由式（4）可以得出，初始沖坑位于下游頂部順坡向下8.93 m 處。計算結果顯示，陡坎在潰壩后0.90 h（對應時間為10:04）進入水庫，現場觀測時間為10:00，計算結果與實測數據吻合較好。另外，根據計算結果，在潰壩后1.38 h（對應時間為10:33）出現潰口峰值流量，這也與實測數據10:25 出現峰值流量相吻合。

表2 給出了射月溝水庫潰壩案例的計算和實測潰壩參數，包括潰口峰值流量（Qp）、潰口最終頂寬（Btf）、潰口最終底寬（Bbf）、潰口最終深度（Df）和潰口峰值流量出現時刻（Tp）。從表2 可以看出，潰口峰值流量、潰口最終頂寬和底寬、潰口峰值流量出現時刻的計算誤差均較小。由于沒有潰口流量過程和潰口形態發展過程的實測數據，僅有推求的峰值流量和實測最終潰口尺寸，圖5 和圖6 主要展示了潰口流量過程和潰口形態發展過程的計算值。另外，從圖5 和圖6 的變化趨勢看出，壩體橫斷面和縱斷面潰口形態的演變趨勢與現場視頻記錄相吻合，這表明本文所提出的數值模擬方法能夠較為準確地描述射月溝水庫的潰壩過程。

表2 射月溝水庫計算與實測潰壩參數值對比Tab.2 Comparison of numerical and measured breaching parameters of Sheyuegou reservior

圖5 潰口流量過程計算值與潰口峰值流量實測值Fig.5 Numerical and measured results of breach flow discharge

圖6 潰口形態發展過程計算值與潰口最終尺寸實測值Fig.6 Numerical and calculated results of breach development

3.3 參數敏感性分析

在射月溝水庫潰壩過程的數值模擬中，壩料沖蝕系數（kd）和陡坎運移系數（CT）是具有高度不確定性的關鍵輸入參數。kd對潰壩過程中的沖蝕率有較大影響，是決定潰壩過程和潰口尺寸的關鍵參數；CT直接決定陡坎的運移速率，也對潰壩過程具有重要影響。通過將數學模型中的相關參數（kd和CT）分別乘以0.5 和2.0，重新計算射月溝水庫的潰壩過程，分析模型分別對這2 個參數的敏感性。表3 和表4 分別顯示了在不同kd和CT條件下，模型計算出的潰口峰值流量、潰口最終尺寸（長度、寬度和深度）和潰口峰值流量出現時刻的變化情況。圖7 和圖8 還分別比較了不同kd和CT條件下的潰口流量過程。

表3 不同壩料沖蝕系數（kd）對潰壩過程的影響Tab.3 Influence of different dam material erodibility coefficients (kd) on dam breach process

表4 不同陡坎運移系數（CT）對潰壩過程的影響Tab.4 Influence of different headcut migration coefficients (CT) on dam breach process

圖7 不同壩料沖蝕系數kd 條件下的潰口流量過程Fig.7 Breach hydrographs under different kd

圖8 不同陡坎運移系數CT 條件下的潰口流量過程Fig.8 Breach hydrographs under different CT

對于沖蝕系數（kd），各潰壩參數的總體敏感性由強到弱依次為：潰口峰值流量，潰口峰值流量出現時刻，潰口最終尺寸。當沖蝕系數增大1 倍或縮小至50%時，潰口峰值流量將分別產生57.5%和?41.7%的變化，而潰口峰值流量對潰壩致災后果評價具有重要意義，因此在選取參數時應注重其合理性。對于陡坎運移系數（CT），由于原始模型中CT的取值較小，當CT縮小至50%時，對計算結果的影響較??；當CT增大1 倍時，各潰壩參數的敏感性有所增強，尤其是潰口峰值流量出現時間，較原始模型提前了23 min。各潰壩參數的總體敏感性由強到弱依次為：潰口峰值流量出現時刻、潰口最終尺寸、潰口峰值流量?？傮w而言，kd主要影響潰口峰值流量和潰口峰值流量出現時刻，而CT主要影響潰口峰值流量出現時刻。

4 與典型潰壩模型比較

選擇美國NWS BREACH 模型[8]模擬射月溝水庫的潰壩過程，然后通過計算結果比較其與本文模型的優缺點。選擇NWS BREACH 模型的原因在于，該模型也是基于潰決機理的簡化漫頂潰壩過程數學模型，且目前在國際上應用較廣泛。表5 列舉了NWS BREACH 模型和本文模型的特點。

表5 NWS BREACH 模型和本文模型的特點總結Tab.5 Summary of characteristics of NWS BREACH model and the present model

NWS BREACH 模型采用修正Meyer-Peter-Mueller 公式[27]模擬壩料的沖蝕過程，其中起控制作用的參數是d90/d30（d90和d30分別為小于該粒徑的顆粒質量占總質量30%和90%對應的顆粒粒徑），反映了壩料顆粒的分布特征。根據射月溝壩料的顆分曲線[4]，d90/d30約為45。模型其他輸入參數可參考表1。

表6 對比了NWS BREACH 模型和本文模型預測的潰壩參數結果及相對誤差。圖9 為2 個模型預測的射月溝水庫潰壩過程中的潰口流量過程線。

表6 不同基于潰決機理的簡化潰壩模型計算結果對比Tab.6 Comparison of numerical results of different simplified physically-based dam breach models

圖9 不同基于潰決機理的簡化潰壩模型預測的潰口流量過程Fig.9 Predicted breach hydrographs of different simplified physically-based dam breach models

由表6 和圖9 可以看出，NWS BREACH 模型預測的潰壩參數和潰口流量過程與本文模型相比表現出更大的相對誤差。由于壩料的沖蝕特性對潰壩過程至關重要，而NWS BREACH 模型中決定沖蝕系數的關鍵參數是d90/d30，射月溝水庫的堆石材料表現出寬級配的特征，且d90/d30= 45，導致壩料的沖蝕率較大，因而預測得到的潰口峰值流量偏大，峰值流量出現時刻提前，潰壩歷時縮短。雖然NWS BREACH模型預測的潰口流量過程線呈現出較大的誤差，但由于潰壩歷時較短，潰口形態后期無明顯擴展，因此潰口最終形態誤差較小。通過射月溝水庫潰壩實例的比較發現，本文模型的預測精度總體高于NWS BREACH 模型。

5 結 語

根據射月溝水庫潰壩過程的視頻記錄和現場調查，本文建立了基于潰決機理的描述水庫潰壩過程的簡化數學模型。利用射月溝水庫潰壩案例的實測資料對模型的準確性進行驗證，并通過與其他典型潰壩模型預測結果的對比，證明了本文模型預測結果的合理性。主要結論如下：

（1）射月溝水庫的潰壩機制主要表現為縱斷面方向初始沖坑的形成與溯源沖蝕，以及橫斷面方向潰口的連續下切、側向沖蝕和潰口邊坡的間歇性失穩坍塌。此外，隨著沖蝕的不斷加劇，縱斷面方向大壩下游坡角逐漸減小，而橫斷面潰口的邊坡角基本保持不變，直至發生潰口邊坡失穩。

（2）參數敏感性分析表明，壩料沖蝕系數對潰壩過程有顯著影響，潰口峰值流量和潰口峰值流量出現時刻對壩料沖蝕系數的敏感性大于潰口最終尺寸，且潰口峰值流量對壩料沖蝕系數的敏感性最高。另外，潰口峰值流量出現時刻比其他潰壩參數對陡坎移動系數更為敏感。因此，在進行潰壩模擬時，壩料沖蝕系數和陡坎移動系數的選取至關重要。

（3）因壩料沖蝕率公式的選擇和對縱、橫斷面潰口發展過程計算方法的過度簡化，目前廣泛應用的NWS BREACH 未能合理描述射月溝潰壩過程的特點。但由于對土石壩潰決過程的認識不足，且測量數據存在較大的不確定性，因此僅通過單個潰壩案例難以全面展示各模型的功能。本文對模型預測能力的對比僅作為參考，未來仍需對各模型的精度開展進一步驗證。