淤地壩淤積對漫頂潰壩洪水影響數值模擬研究

張兆安， 侯精明， 劉占衍， 馬利平， 李占斌， 李小綱

(1.西安理工大學 省部共建西北旱區生態水利國家重點實驗室， 陜西 西安 710048；2.河北省邯鄲水文水資源勘測局， 河北 邯鄲 056006)

1 研究背景

為了治理黃土高原嚴重的水土流失問題，淤地壩作為最有效的防治水土流失的溝道工程措施[1]，在我國黃河流域被大量建造并使用[2]，它能夠防洪攔沙、淤地增產、充分利用水沙資源和改善生態環境[3]，通過抬高溝底侵蝕基準面，提高了流域斜坡穩定性，減少了溝坡重力侵蝕發生的可能性[4]。但在實際運行中，隨著淤地壩淤積高度變高，漫頂潰壩問題變得越來越頻發。由于氣候變化，我國極端降雨導致的致災暴雨事件增加[5]，僅在1994年，陜北就有7347座淤地壩在暴雨中被摧毀[6]。因此，研究淤地壩潰壩洪水的影響具有十分重大的意義。

近年來我國西北淤地壩事故頻發，國內關于淤地壩的研究日漸增多，王保勤[7]分析了滲水在淤地壩災害中的危害，在淤地壩設計施工等方面提出建議；魏霞等[6]認為從長時間尺度、大范圍來看水毀不是件壞事，淤地壩災害防治研究對淤地壩的建設有著重要意義，肯定了對淤地壩水毀潰壩的研究；許五弟等[8]按照不同的暴雨洪水和庫容狀態對淤地壩進行了模擬分析，驗證了淤地壩潰壩預報預警地理信息模型應用的可行性，但沒有直觀體現淤地壩的潰壩規律；王乃欣等[9]利用虛擬現實軟件，在地形不變的情況下改變初始水位，對淤地壩潰決進行仿真研究，使淤地壩潰壩結果更加直觀，但并未研究淤積高度對于淤地壩潰壩的影響；高海東等[10]用SINMAP模型定量評估了淤地壩淤積高度對斜坡穩定性的影響，研究了淤地壩淤積至不同高度時的斜坡穩定性，為淤地壩建設提供參考。但以上對于淤地壩災害的研究都缺乏對淤地壩淤積高度與其潰壩洪水關系的研究。

為研究不同運行周期淤積高度改變下潰壩洪水特征的變化規律，本文采用耦合了潰口演變過程的二維水動力模型，模擬了淤積高度分別為20%、40%、60%、80%壩高情形下的淤地壩漫頂潰決及下游洪水演進過程，分析了淤積高度對潰壩洪水峰值的量化影響。研究成果可用于指導淤地壩運行維護及防洪減災工作。

2 研究方法

本文采用考慮潰口演變過程的DB-IWHR模型計算潰口演變下的潰口處流量過程，并使用潰壩洪水演進模型(GAST模型)計算下游洪水演進過程，潰壩過程模擬流程見圖1。

圖1 潰壩過程模擬流程圖

GAST模型即顯卡加速的地表水及其伴隨輸移過程模型(GPU Accelerated Surface Water Flow and Transport Model (GAST))，該模型基于Godunov的有限體積法求解二維淺水方程，該方法可以穩定地解決不連續問題，嚴格保持物質守恒，應用二階算法提高了模擬的計算效率和精度。它還耦合并模擬模型中淺水流動的一些其他物理化學反應，如沉積物和污染物運輸。該模型的另一個特征就是使用GPU計算來提高計算性能。

2.1 潰壩洪水演進模型控制方程

平面二維淺水方程(SWEs)又稱圣維南方程，被廣泛用于模擬自然河流、湖泊和沿海地區的淺水流動，同時它還可以描述潰壩問題。二維淺水方程可以寫成：

(1)

其中變量矢量、不同方向上的通量矢量以及源項矢量可以表示如下：

(2)

式中：q為變量矢量，包括水深h， m；qx和qy為x、y方向的單寬流量，m3/s；F和G分別為x、y方向的通量矢量;g為重力加速度，m3/s；u和v分別為x、y方向的流速，m3/s；S為源項矢量；zb為河床底面高程，m；Cf=gn2/h1/3為謝才系數，m1/2/s。

2.2 潰壩洪水演進模型數值方法

GAST模型采用Godunov格式的有限體積法對平面二維淺水方程進行求解。選用具有二階時空精度的MUSCL型格式以初始水深來進行空間二階精度的變量重組，以重組后的數值作為計算變量，使數值變量守恒。在控制單元內，界面上的物質與動量通量通過HLLC近似黎曼求解器計算[11]；通過靜水重構來修正干濕邊界處負水深；底坡源項用底坡通量法進行處理。摩阻源項(摩阻力)使用較為穩定的半隱式法計算，以此來提高穩定性[12]，并采用二階顯式Runge Kutta[13]方法進行時間步長的推進。

本模型模擬潰壩過程時，采用GPU并行計算技術實現高速運算[14]，以此來提高模型計算效率，但模型未考慮潰壩時的土動力過程。

2.3 潰口演變模型

本文使用DB-IWHR模型(DB模型)來模擬考慮潰口演變的情況下的潰口流量過程。DB模型由中國水利水電科學研究院開發，考慮了潰口斷面的流量會受到水流在垂直方向沖刷影響，也會受到側向擴展影響的問題[15]。由于下游河床很低，潰口出口水流按自由出流考慮，這時潰口斷面流量可用寬頂堰公式計算[16]：

Q=CB(H-z)3/2

(3)

式中：Q為流量，m3/s；B為斷面寬度，m；H為庫水位，m；z為進口處河床高程，m；C為綜合流量系數，經驗取值為1.4 m1/2/s[16]。

關于沖刷計算，DB-IWHR模型提出一種雙曲型侵蝕率模型，認為水流沖刷土石料時，土石料抵抗沖刷侵蝕的能力并非是無限制的[15]。關于潰口擴展計算，模型采用簡化Bishop[17]法自動搜索臨界圓弧滑裂面，對潰口擴展進行分析，相關原理和計算方法在文獻[18]中有詳細說明，本文不再贅述。

該耦合模型在文獻[19]中模擬計算實際潰壩洪水已經得到驗證。

3 淤積高度對潰壩洪水影響模擬

為研究不同運行周期淤積高度改變下潰壩洪水特征的變化規律，本文采用耦合了潰口演變過程的二維水動力模型，分別對矩形溝道、王茂溝某淤地壩模擬其淤積高度為20%、40%、60%、80%壩高情形下淤地壩漫頂潰決及下游洪水演進過程，其中矩形溝道及王茂溝某淤地壩壩高皆為10 m，并結合實際工程淤地壩壩高，將模擬壩高提升至20、30 m進行模擬，對所得結論進行一般性地拓展。

3.1 矩形溝道淤地壩潰決洪水過程模擬

為了研究淤積高度對潰壩洪水的影響，本文設計一個1 000 m×200 m的矩形溝道，以左下角作為零點，在x正方向795～800 m處為淤地壩(圖2)。

圖2 矩形溝道示意圖(單位：m)

3.1.1 矩形溝道潰壩方案 本文設定淤地壩壩高10 m，壩寬和壩長為5和200 m。壩體下游是一個坡面，坡度取1∶0.6。在漫頂的條件下，對不同淤積高度的淤地壩進行潰壩模擬，并在下游截取代表斷面(圖2中以S1、S2表示)進行研究。

3.1.2 DB-IWHR模型參數設置 用DB-IWHR模型建立潰口關系并計算潰口流量，再通過GAST模型計算下游河道洪水演進過程，設定曼寧系數為0.025[20]，下游為自由出流開邊界，其余為閉邊界，參考DB模型使用說明書及相關文獻[21]，設定參數如表1所示。

表1 DB-IWHR模型參數設定

設定入庫流量Qin=0，保持初始水位H0=10 m始終不變。庫容—水位關系滿足以下方程：

dw=(Qn-Qin)dt

(4)

dH=dw/(200×200)

(5)

式中： dw為水量，m3； dH為水位，m; dt為時間步長，s;Qn為流量，m3/s。

3.1.3 矩形溝道潰壩模擬結果 圖3和4分別為矩形溝道在庫區淤積高度不同的情況下，下游代表斷面的潰壩洪水流量過程線，以及兩斷面處洪峰流量與淤積高度擬合曲線(此處取均方根誤差最小的二次多項式擬合)。由圖3、4可以看出，洪峰流量隨淤積高度的增加而下降，且下降趨勢逐漸變緩。其中，S1處二次多項式均方根誤差為1.92 m3/s，相關系數為0.999；S2處二次多項式均方根誤差為2.87 m3/s，相關系數為0.999。

表2為不同淤積高度下代表斷面的洪峰流量及洪峰流量削減率表，淤積高度至20%時，兩斷面的洪峰流量削減率增加至44.4%左右；淤積高度至40%時，兩斷面洪峰流量削減率增加至72.0%左右，故洪峰流量隨淤積高度的增加而減小，洪峰流量削減率隨淤積高度增大而增大。

3.2 王茂溝某淤地壩潰決過程模擬

3.2.1 研究區域概況 王茂溝位于陜西省榆林市綏德縣，是韭園溝流域的一個分支，流域面積5.97 km2，其中溝間地占58.4%；溝谷地占41.6%，主溝長3.75 m[22]。自1953年開始綜合治理后，建造淤地壩40多座，已攔泥超150×104m3，在溝內已經形成較完整的壩系[23]。

根據激光雷達掃描技術獲取生成的1 m精度數字地形高程文件，選擇一處典型小支溝，潰口完好，垂直于下游水流流向截取兩個代表斷面D1、D2，研究區支溝地形及代表斷面位置如圖5所示。

表2 代表斷面洪峰流量及洪峰流量削減率

3.2.2 王茂溝潰壩洪水過程模擬結果 對王茂溝淤地壩地形及相關參數做類似矩形河道的處理，下游為自由出流開邊界，其余為閉邊界，在壩后的庫區內，設定初始水位始終為10 m水位，不隨淤積高度發生改變。模擬使用了高精度實際地形，經過模擬可以得到三維展示圖(圖6)及計算結果。

圖7、8分別為王茂溝淤地壩在庫區淤積高度不同的情況下，下游代表斷面的潰壩洪水流量過程線，以及代表斷面處洪峰流量與淤積高度擬合曲線，得到洪峰流量與淤積高度可用二次多項式進行擬合。由圖8可看出，洪峰流量隨淤積高度的增加而下降，且下降趨勢逐漸變緩。其中，D1處二次多項式均方根誤差7.359 m3/s，相關系數為0.995；D2處二次多項式均方根誤差6.686 m3/s，相關系數為0.996。表3為代表斷面不同淤積高度下的洪峰流量及洪峰流量削減率，淤積高度至20%時，兩斷面的洪峰流量削減率增加至48.2%左右；淤積高度至40%時，兩斷面洪峰流量削減率增加至74.0%左右，故可得出實際淤地壩地形與矩形溝道相似的結論。

表3 代表斷面洪峰流量及洪峰流量削減率表

圖3 代表斷面處潰壩洪水流量過程線

圖4 代表斷面處洪峰流量與淤積高度擬合曲線

圖5 研究區域地形及代表斷面位置圖

圖6 潰壩洪水演進三維展示圖

3.3 不同壩高下的淤積程度對潰壩洪水過程影響模擬

由以上分析得出，當淤地壩壩高為10m時，在各個算例中均得到統一的潰壩洪水演變規律，故本文再次計算理想溝道地形條件下壩高為20和30 m的情形，淤積高度每次增加量為淤地壩壩高的20%，結果如下：

圖9、10分別為不同壩高的情況下，下游代表斷面D1的潰壩洪水流量過程線，以及代表斷面D1處洪峰流量與淤積高度擬合曲線。可得出與10 m壩高時類似的結論，壩高20 m時二次多項式均方根誤差為4.641 m3/s，相關系數為0.999；壩高30 m時二次多項式均方根誤差為16.435 m3/s，相關系數為0.998。

表4為不同壩高的情況下，代表斷面D1不同淤積高度下的洪峰流量及洪峰流量削減率，淤積高度增至壩高的20%時，兩種壩高的洪峰流量削減率平均增加至34.7%左右；淤積高度至40%時，兩種壩高洪峰流量削減率平均增加至64.4%左右。得到結論與壩高為10 m時類似，即洪峰流量隨淤積高度的增加而減小，洪峰流量削減率隨淤積高度增大而增大。

圖7 代表斷面處潰壩洪水流量過程線

圖8 代表斷面處洪峰流量與淤積高度擬合曲線

圖9 不同壩高代表斷面D1處潰壩流量過程線

圖10 不同壩高代表斷面D1處洪峰流量與淤積高度擬合曲線

表4 代表斷面D1洪峰流量及洪峰流量削減率表

4 結 論

為研究不同運行周期淤地壩淤積高度改變下漫頂潰壩洪水特征的變化規律，本文對不同淤積高度下的淤地壩潰壩洪水進行數值模擬研究，結論如下：

(1)淤地壩潰壩洪峰流量隨淤積高度增加而減小，并可用二次多項式擬合，相關系數均在0.99左右，擬合精度較高，研究成果對淤地壩潰壩災害預報有一定意義。

(2)在淤地壩淤積初期和中期，淤積高度對于潰壩洪水影響較大，淤積高度至壩高20%和40%時，平均洪峰流量削減率分別達40.50%和68.71%，淤地壩淤積對于降低潰壩洪水災害影響效果顯著。

(3)本文采用考慮了潰口演變的水動力模型，分別模擬計算了在不同地形、不同壩高情形下不同淤積高度淤地壩漫頂潰決及下游洪水演進過程，分析了淤地壩淤積對潰壩洪水的影響，各算例得到的規律一致，建議在淤地壩建壩完成初期(淤地厚度在壩高20%左右)加強對于潰壩災害的預防措施。

本文研究了不同運行周期淤積高度改變下潰壩洪水特征的變化規律，為淤地壩建設、維護及防洪減災提供了理論依據。本次研究未考慮泥沙輸移，水庫來流等因素對潰壩洪水的影響以及管涌潰壩的情形，將在未來研究中加入未考慮因素，使結果更加合理。