淤地壩淤積對(duì)漫頂潰壩洪水影響數(shù)值模擬研究

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

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710048；2.河北省邯鄲水文水資源勘測(cè)局， 河北 邯鄲 056006)

1 研究背景

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

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

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

2 研究方法

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

圖1 潰壩過(guò)程模擬流程圖

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

2.1 潰壩洪水演進(jìn)模型控制方程

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

(1)

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

(2)

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

2.2 潰壩洪水演進(jìn)模型數(shù)值方法

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

本模型模擬潰壩過(guò)程時(shí)，采用GPU并行計(jì)算技術(shù)實(shí)現(xiàn)高速運(yùn)算[14]，以此來(lái)提高模型計(jì)算效率，但模型未考慮潰壩時(shí)的土動(dòng)力過(guò)程。

2.3 潰口演變模型

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

Q=CB(H-z)3/2

(3)

式中：Q為流量，m3/s；B為斷面寬度，m；H為庫(kù)水位，m；z為進(jìn)口處河床高程，m；C為綜合流量系數(shù)，經(jīng)驗(yàn)取值為1.4 m1/2/s[16]。

關(guān)于沖刷計(jì)算，DB-IWHR模型提出一種雙曲型侵蝕率模型，認(rèn)為水流沖刷土石料時(shí)，土石料抵抗沖刷侵蝕的能力并非是無(wú)限制的[15]。關(guān)于潰口擴(kuò)展計(jì)算，模型采用簡(jiǎn)化Bishop[17]法自動(dòng)搜索臨界圓弧滑裂面，對(duì)潰口擴(kuò)展進(jìn)行分析，相關(guān)原理和計(jì)算方法在文獻(xiàn)[18]中有詳細(xì)說(shuō)明，本文不再贅述。

該耦合模型在文獻(xiàn)[19]中模擬計(jì)算實(shí)際潰壩洪水已經(jīng)得到驗(yàn)證。

3 淤積高度對(duì)潰壩洪水影響模擬

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

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

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

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

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

3.1.2 DB-IWHR模型參數(shù)設(shè)置 用DB-IWHR模型建立潰口關(guān)系并計(jì)算潰口流量，再通過(guò)GAST模型計(jì)算下游河道洪水演進(jìn)過(guò)程，設(shè)定曼寧系數(shù)為0.025[20]，下游為自由出流開(kāi)邊界，其余為閉邊界，參考DB模型使用說(shuō)明書(shū)及相關(guān)文獻(xiàn)[21]，設(shè)定參數(shù)如表1所示。

表1 DB-IWHR模型參數(shù)設(shè)定

設(shè)定入庫(kù)流量Qin=0，保持初始水位H0=10 m始終不變。庫(kù)容—水位關(guān)系滿足以下方程：

dw=(Qn-Qin)dt

(4)

dH=dw/(200×200)

(5)

式中： dw為水量，m3； dH為水位，m; dt為時(shí)間步長(zhǎng)，s;Qn為流量，m3/s。

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

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

3.2 王茂溝某淤地壩潰決過(guò)程模擬

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

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

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

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

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

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

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

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

圖5 研究區(qū)域地形及代表斷面位置圖

圖6 潰壩洪水演進(jìn)三維展示圖

3.3 不同壩高下的淤積程度對(duì)潰壩洪水過(guò)程影響模擬

由以上分析得出，當(dāng)淤地壩壩高為10m時(shí)，在各個(gè)算例中均得到統(tǒng)一的潰壩洪水演變規(guī)律，故本文再次計(jì)算理想溝道地形條件下壩高為20和30 m的情形，淤積高度每次增加量為淤地壩壩高的20%，結(jié)果如下：

圖9、10分別為不同壩高的情況下，下游代表斷面D1的潰壩洪水流量過(guò)程線，以及代表斷面D1處洪峰流量與淤積高度擬合曲線?？傻贸雠c10 m壩高時(shí)類似的結(jié)論，壩高20 m時(shí)二次多項(xiàng)式均方根誤差為4.641 m3/s，相關(guān)系數(shù)為0.999；壩高30 m時(shí)二次多項(xiàng)式均方根誤差為16.435 m3/s，相關(guān)系數(shù)為0.998。

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

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

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

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

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

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

4 結(jié) 論

為研究不同運(yùn)行周期淤地壩淤積高度改變下漫頂潰壩洪水特征的變化規(guī)律，本文對(duì)不同淤積高度下的淤地壩潰壩洪水進(jìn)行數(shù)值模擬研究，結(jié)論如下：

(1)淤地壩潰壩洪峰流量隨淤積高度增加而減小，并可用二次多項(xiàng)式擬合，相關(guān)系數(shù)均在0.99左右，擬合精度較高，研究成果對(duì)淤地壩潰壩災(zāi)害預(yù)報(bào)有一定意義。

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

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

本文研究了不同運(yùn)行周期淤積高度改變下潰壩洪水特征的變化規(guī)律，為淤地壩建設(shè)、維護(hù)及防洪減災(zāi)提供了理論依據(jù)。本次研究未考慮泥沙輸移，水庫(kù)來(lái)流等因素對(duì)潰壩洪水的影響以及管涌潰壩的情形，將在未來(lái)研究中加入未考慮因素，使結(jié)果更加合理。