庫區(qū)幾何要素對潰壩洪水特性影響的模擬研究

黃子奇， 張建豐， 李 濤， 胡 暉， 高澤海， 杜佰林

(西安理工大學(xué) 西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710048)

水壩潰決引起的洪水通常比暴雨引起的水文洪水更具破壞性，會(huì)對人民生命財(cái)產(chǎn)和國民經(jīng)濟(jì)造成極其嚴(yán)重的危害[1]。對潰壩洪水進(jìn)行研究有助于人們計(jì)算分析潰壩洪水的危害，采取適當(dāng)?shù)念A(yù)防措施。在其他科學(xué)研究中，潰壩波被用于模擬海嘯波、風(fēng)暴潮、涌潮等災(zāi)害，基于潰壩水流的緩蓄快放技術(shù)也被應(yīng)用于景觀水體水質(zhì)的保持[2-3]。因此，對潰壩洪水的研究意義重大。

對潰壩洪水的理論研究可歸結(jié)為求解Saint-Venant方程組的Riemann問題，以精確描述潰壩水流運(yùn)動(dòng)規(guī)律。Ritter[4]提出了無限庫區(qū)長度情況下干底矩形棱柱體河槽中瞬潰的解析解。Dressier[5]提出了考慮流動(dòng)阻力的Ritter解。伍超等[6]用Riemann方程推導(dǎo)出平底無阻力無限長度U形棱柱體河道潰壩波的簡化解。王立輝等[7]利用特征函數(shù)和Riemann函數(shù)建立了斜底河道瞬潰的解析解。這些研究都是控制方程在簡化條件下的解析解，實(shí)際應(yīng)用受到了限制，潰壩洪水的理論研究仍然非常重要。

物理模型實(shí)驗(yàn)?zāi)軌蛐：藵尾ɡ碚摻?，研究洪水演進(jìn)特性、其影響范圍以及流態(tài)，是研究潰壩洪水的重要方法。Lauber等[8-9]通過矩形水槽實(shí)驗(yàn)研究了下游無水情況下潰壩波在平底河道與斜底河道的演進(jìn)規(guī)律，得出了下游最大水深、流速和水深分布的規(guī)律等水力特性。Pilotti等[10]研究了庫區(qū)平面形狀對潰壩洪水的影響，總結(jié)了庫區(qū)形狀對下游水深的影響以及經(jīng)驗(yàn)公式的適用性。Hooshyaripor等[11]論證了水庫橫截面邊坡坡度對潰壩洪水洪峰流量、峰現(xiàn)時(shí)間以及下游水深隨時(shí)間演變的重要作用。劉文軍等[12]對上下游初始水深比為0～0.9時(shí)的瞬時(shí)潰壩水流進(jìn)行了系統(tǒng)研究，分析了瞬潰洪水的運(yùn)動(dòng)特性并對濕底情況下潰壩水流的相似性進(jìn)行了分析。目前的研究對庫區(qū)寬度、長度等幾何因素對潰壩洪水的影響關(guān)注較少，且關(guān)注點(diǎn)多集中在各斷面水深過程的分析。

隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，通過數(shù)學(xué)模型模擬分析復(fù)雜流場成為研究熱點(diǎn)，與物理模型相比，具有造價(jià)低、容易調(diào)整實(shí)驗(yàn)條件、能避免外界因素干擾等優(yōu)點(diǎn)[13]。BRDAM、BREACH、DAMBREAK、BEED等模型已用于簡單潰壩問題的計(jì)算[14]。夏軍強(qiáng)等[15]建立了基于無結(jié)構(gòu)三角網(wǎng)格下采用有限體積法求解的二維水動(dòng)力學(xué)模型，用于模擬潰壩洪水在復(fù)雜邊界及實(shí)際地形上的流動(dòng)過程。胡暉等[16]利用MIKE 3 FM與Flow-3D對潰壩洪水進(jìn)行數(shù)值模擬，與實(shí)測實(shí)驗(yàn)數(shù)據(jù)和解析解比較，研究了兩種軟件的適用性，得出了MIKE 3 FM模型適用于大型區(qū)域中實(shí)際潰壩問題的計(jì)算，F(xiàn)low-3D模型適用于精細(xì)計(jì)算。壩址是上、下游的連接點(diǎn)，使用數(shù)值解法研究潰壩洪水時(shí)被作為下游的邊界條件使用，準(zhǔn)確推求壩址斷面流量過程線非常重要，然而在實(shí)際工程中，壩址斷面流量過程線多通過概化算法獲得，一般默認(rèn)潰壩開啟的時(shí)刻即為壩址峰現(xiàn)時(shí)間。使用Flow-3D進(jìn)行數(shù)值模擬獲得的壩址斷面流量過程線包含流量增大至洪峰流量的過程，與概化算法相比有一定的優(yōu)勢。

現(xiàn)有研究缺乏針對庫區(qū)寬度、長度等幾何因素對潰壩洪水特性影響的關(guān)注，且實(shí)驗(yàn)分析較少涉及流量和波速。本文應(yīng)用Flow-3D構(gòu)建水槽的數(shù)值模型，提出了獲取物理水槽實(shí)驗(yàn)流量數(shù)據(jù)的方法，通過物理模型實(shí)驗(yàn)流量數(shù)據(jù)驗(yàn)證了數(shù)值模型的合理性，通過數(shù)值模型研究庫區(qū)蓄水深度hu、庫區(qū)寬度B、庫區(qū)長度L對壩址斷面流量過程、潰壩洪水水面線、潰壩洪水洪峰流量、峰現(xiàn)時(shí)間、潰壩波波速等特性的影響，對潰壩洪水的模型實(shí)驗(yàn)作出一定的補(bǔ)充。

1 實(shí)驗(yàn)工具和方法

物理模型實(shí)驗(yàn)在西安理工大學(xué)露天試驗(yàn)場的矩形棱柱體水槽內(nèi)進(jìn)行，實(shí)驗(yàn)時(shí)間是2018年5月中旬。水槽形式見圖1，水槽以有機(jī)玻璃制成，底坡為0.002 5，糙率為0.01，寬0.3 m，高0.3 m，長10.8 m，距上游邊界2 m處設(shè)置一塊平板閘門將水槽分為上游蓄水庫區(qū)與下游河道，即上游蓄水庫區(qū)長度為2 m。壩址(閘門)上游設(shè)置10個(gè)觀測斷面，壩址下游設(shè)置6個(gè)觀測斷面，壩址處設(shè)置1個(gè)觀測斷面，以壩址為原點(diǎn)，潰壩洪水向壩址下游傳播為正方向建立數(shù)軸，各監(jiān)測斷面的坐標(biāo)見表1。

表1 監(jiān)測斷面坐標(biāo)

進(jìn)行實(shí)驗(yàn)時(shí)，其無量綱化閘門提升時(shí)間為(g/hu)1/2t<21/2，滿足瞬間潰壩條件(g為重力加速度，取9.8 m/s2；hu為庫區(qū)蓄水深度，m)[8]。

實(shí)驗(yàn)水深數(shù)據(jù)采用錄像法獲得。使用高速相機(jī)對準(zhǔn)觀測斷面，調(diào)節(jié)焦距使觀測斷面的水深刻度與秒表在拍攝范圍內(nèi)清晰可見，錄像速度為30 fps，開始錄像后，在快速打開閘門的同時(shí)按下秒表計(jì)時(shí)按鍵，待觀測斷面處水位不再變化時(shí)停止錄像。按幀查看視頻即可獲得觀測斷面的水深隨時(shí)間變化的過程，按相同時(shí)刻匯總17個(gè)觀測斷面的水深即可獲得不同時(shí)刻的水面線變化情況。實(shí)驗(yàn)中對每個(gè)工況都重復(fù)3次進(jìn)行，水深數(shù)據(jù)取3次實(shí)驗(yàn)的平均值得到。

由于流量數(shù)據(jù)無法直接測得，本文根據(jù)相鄰兩時(shí)刻水面線與邊界圍成的水體體積變化進(jìn)行推求，計(jì)算公式為：

(1)

式中：Δt為時(shí)間間隔；V(t)為t時(shí)刻觀測斷面與水槽邊界、觀測斷面上游的水面圍成的水體體積；Q(t+Δt/2)為t+Δt/2時(shí)刻觀測斷面的流量。

2 數(shù)學(xué)模型

2.1 公式和數(shù)學(xué)符號

在潰壩洪水?dāng)?shù)值模擬中，一般不考慮水流壓縮性， Flow-3D采用如下方程控制流體運(yùn)動(dòng)。

連續(xù)性方程：

(2)

動(dòng)量方程：

(3)

紊流模型采用采用RNGk-ε模型，其k、ε方程見下。

k方程：

(4)

紊動(dòng)耗散率ε方程：

(5)

(6)

(7)

(8)

(9)

(10)

式中：ui、uj為流速分量，i、j= 1，2，3；xi、xj為坐標(biāo)分量；ρ為體積分?jǐn)?shù)平均密度；μ為流體動(dòng)力黏滯系數(shù)，μt為流體湍動(dòng)黏度；P為作用在流體微元上的壓力；k為紊動(dòng)能；Gk為由于流速梯度引起的紊動(dòng)能產(chǎn)生項(xiàng)；ε為紊動(dòng)能耗散率；η為紊動(dòng)時(shí)間與應(yīng)變尺度的比值；Eij為時(shí)均應(yīng)變率；本研究取Cμ、Cε1、Cε2、σk、σε、η0、β等經(jīng)驗(yàn)常數(shù)為軟件默認(rèn)值。

Flow-3D采用流體體積法(VOF)追蹤模型沿程自由液面變化。VOF的基本原理是通過研究網(wǎng)格單元中流體和網(wǎng)格體積比函數(shù)F來構(gòu)造和追蹤自由液面[17]。若F=1，則該單元全部為液體所占據(jù)；若F=0，則該單元為無液體單元；當(dāng)0

(11)

式中參變量含義同上。

Flow-3D采用有限差分法通過將模型和計(jì)算域獨(dú)立進(jìn)行網(wǎng)格生成，將方程離散從而轉(zhuǎn)化為代數(shù)方程組并進(jìn)行數(shù)值求解。其獨(dú)特的Favor網(wǎng)格處理技術(shù)對三維實(shí)體做六面體劃分，能夠消除和順滑不平整區(qū)域，可解決網(wǎng)格模型失真問題。

2.2 模型建立

本研究利用Flow-3D按照物理模型的形式構(gòu)建數(shù)值模型，具體參數(shù)見表2。利用工況1、2、3可對比分析庫區(qū)蓄水深度hu對潰壩洪水的影響，利用工況3、4、5可對比分析庫區(qū)寬度B對潰壩洪水的影響，利用工況3、6、7可對比分析庫區(qū)長度L對潰壩洪水的影響。工況4與工況6的庫容相同，工況5與工況7的庫容相同，對比這些工況可知庫區(qū)長度與寬度對潰壩洪水水力要素作用的差異。綜合考慮求解精度與收斂速度，模型的網(wǎng)格尺寸設(shè)置為5 mm。

表2 數(shù)值模型參數(shù)Tab.2 Parameters of the numerical model

2.3 模型驗(yàn)證

為了驗(yàn)證數(shù)值模型的準(zhǔn)確性，本文將物理模型實(shí)驗(yàn)與數(shù)值模擬結(jié)果進(jìn)行對比。選取壩址上游距離壩址1 m處的斷面、壩址斷面、壩址下游距離壩址3 m 處的斷面對比數(shù)學(xué)模型與物理模型的水深、流量隨時(shí)間變化的過程，結(jié)果見圖2。觀察得出，在庫區(qū)蓄水深度為9 cm的工況下，將閘門快速提起，河道水體的突然外泄，閘門附近水位迅速下降，形成洪水波向下游演進(jìn)。當(dāng)洪水波傳播至某一斷面，該斷面水深快速升高，一段時(shí)間后水位逐漸下降，F(xiàn)low-3D模擬的結(jié)果與觀測結(jié)果吻合。其中，圖2(c)中下游斷面0 s至3 s時(shí)間內(nèi)水深不為0，是因?yàn)殚l門有少量漏水。上游斷面、壩址斷面、下游斷面水深數(shù)據(jù)的相對誤差分別為6.66%、9.11%、9.88%，流量數(shù)據(jù)的相對誤差分別為9.85%、10.63%、10.23%，說明模型模擬精度基本滿足要求，本研究采用的數(shù)值模擬方法合理，可用于潰壩洪水的模擬研究。

圖2 各斷面水深、流量變化Fig.2 Depth and discharge change for each section

3 模擬結(jié)果分析

3.1 壩址斷面流量過程

圖3為各種工況下的壩址流量過程線。閘門瞬間移除時(shí)，產(chǎn)生引起水位上升的正波向壩址下游傳播、引起水位下降的負(fù)波向壩址上游傳播。由于此時(shí)上下游水位差最大，壩址斷面瞬間到達(dá)最大流量。隨著洪水向壩址下游演進(jìn)，壩址上游水位降低，流量隨之減小。待上游水位基本穩(wěn)定后，上下游水位差將在一段時(shí)間內(nèi)保持穩(wěn)定，流量也保持穩(wěn)定。一段時(shí)間后，上下游水位差不斷減小，流量隨之減小，直到上下游水位相等時(shí)，流量減小為零。

圖3 壩址斷面流量過程線(x=0)Fig.3 Discharge process in damsite(x=0)

工況1、2、3的峰現(xiàn)時(shí)間分別是0.12 s、0.14 s、0.16 s，洪峰流量分別為0.046 m3/s，0.071 m3/s，0.099 m3/s。0～5 s之間有一個(gè)流量基本不變的穩(wěn)定階段，其流量穩(wěn)定時(shí)間約為2.69 s、1.97 s、1.60 s，穩(wěn)定流量約為0.039 m3/s、0.060 m3/s、0.083 m3/s。可見隨著庫區(qū)蓄水深度的增加，壩址處洪峰流量增大，峰現(xiàn)時(shí)間推遲，流量穩(wěn)定持續(xù)時(shí)間減少，穩(wěn)定流量增大。

工況3、4、5的洪峰流量分別為0.099 m3/s，0.132 m3/s，0.165 m3/s，峰現(xiàn)時(shí)間均是0.16 s。流量穩(wěn)定時(shí)間約為1.60 s、1.59 s、1.61 s，穩(wěn)定流量約為0.083 m3/s、0.111 m3/s、0.139 m3/s。可見，隨著庫區(qū)寬度的增加，壩址處洪峰流量增大，峰現(xiàn)時(shí)間基本不變，流量穩(wěn)定持續(xù)時(shí)間基本不變，穩(wěn)定流量增大。

工況3、6、7的洪峰流量均為0.099 m3/s，峰現(xiàn)時(shí)間均為0.16 s，流量穩(wěn)定時(shí)間約為1.60 s、2.43 s、3.36 s，穩(wěn)定流量約為0.083 m3/s?？梢?，隨著庫區(qū)長度的增加，洪峰流量、峰現(xiàn)時(shí)間均不變，流量穩(wěn)定持續(xù)時(shí)間增長，穩(wěn)定流量不變。

根據(jù)Ritter解的壩址洪峰流量公式，壩址洪峰流量與庫區(qū)寬度、庫區(qū)蓄水深度成正比，與庫區(qū)長度無關(guān)，與本研究相符。庫區(qū)長度的變化對壩址斷面的峰現(xiàn)時(shí)間沒有影響，與穩(wěn)定流量的持續(xù)時(shí)間成正比。各種工況下，穩(wěn)定流量均為洪峰流量的8/9，此時(shí)對應(yīng)的水深為庫區(qū)蓄水深度的4/9，根據(jù)Ritter解，壩址斷面將會(huì)長久維持在這一水深，研究結(jié)果與Ritter解不同的主要原因是Ritter解的基本假定是庫區(qū)擁有無限長度，而本研究是基于有限長度的水槽實(shí)驗(yàn)，因此在一段時(shí)間后，流量與水深均會(huì)下降，這也說明了本研究的合理性[4]。

3.2 潰壩洪水水面線

圖4為不同工況下壩址下游水面線。隨著庫區(qū)蓄水深度的增加，潰壩波附加比降增大，各斷面最大水深逐漸增大(圖4(a))。隨著庫區(qū)寬度的增大，水面線基本重合，潰壩波附加比降和各斷面最大水深基本不變(圖4(b))。隨著庫區(qū)長度的增大，各斷面最大水深增大，庫區(qū)長度與達(dá)到最大水深時(shí)間成正比。潰壩洪水向下游演進(jìn)的初始階段，較小的庫區(qū)長度達(dá)到最大水深之前水面線基本重合，潰壩波附加比降基本不變。之后水面線逐漸下降，與其他工況的水面線偏離，其他工況的水面線仍然重合。隨著潰壩洪水向下游演進(jìn)，潰壩波附加比降逐漸變緩，沿程各斷面最大水深逐漸降低，這是由于河道具有調(diào)蓄作用，潰壩波在傳播的過程中將逐漸坦化。

圖4 潰壩洪水水面線Fig.4 Dam-break flood water surface profile

3.3 洪峰流量

圖5為各斷面洪峰流量。隨著庫區(qū)蓄水深度、庫區(qū)長度的增加，沿程各斷面洪峰流量Qm均增大；隨著庫區(qū)寬度的增大，各斷面洪峰流量均增大。

圖5 各斷面洪峰流量Fig.5 Peak discharge in each section

為了說明庫區(qū)蓄水深度、庫區(qū)寬度、庫區(qū)長度對潰壩洪水洪峰減小的影響，用各斷面洪峰流量與壩址斷面洪峰流量的比值來表示，即為α，各工況下各斷面α值見表3。由工況1～5的α值可知，隨著庫區(qū)蓄水深度與庫區(qū)寬度的增大，沿程各斷面α基本不變，即庫區(qū)蓄水深度與庫區(qū)寬度的變化對下游的洪峰減小沒有影響。對比工況3與工況7，庫區(qū)長度由3 m增大到5 m，沿程各斷面α平均增加了8.2%，說明庫區(qū)長度的增大均對潰壩洪水的洪峰減小有抑制作用。

表3 不同工況下各斷面α值Tab.3 α value of each section under different conditions

3.4 各斷面波前到達(dá)時(shí)間

圖6為各斷面波前到達(dá)時(shí)間。隨著庫區(qū)蓄水深度的增加，各斷面的波前到達(dá)時(shí)間t1均明顯提前。隨著庫區(qū)寬度或長度的增加，各斷面的波前到達(dá)時(shí)間基本不變。說明潰壩波波速只與庫區(qū)蓄水深度有關(guān)，Ritter的研究結(jié)果與此相似[4]。主要原因是，水體流動(dòng)是因?yàn)閱挝恢亓克w的位能與壓能轉(zhuǎn)化為動(dòng)能，位能和壓能只與水深有關(guān)，與庫區(qū)長度和寬度無關(guān)。

圖6 各斷面波前到達(dá)時(shí)間Fig.6 Wave front arrival time in each section

見圖3，隨著庫區(qū)蓄水深度的增加，壩址斷面的流量過程線的穩(wěn)定時(shí)間縮短，庫區(qū)寬度的變化對穩(wěn)定時(shí)間基本沒有影響。隨著庫區(qū)長度的增加，穩(wěn)定時(shí)間增長，主要原因是負(fù)波到達(dá)水槽上邊界，被反射再傳播至壩址斷面的過程流量保持穩(wěn)定，之后再逐漸減小，庫區(qū)長度越大則負(fù)波的傳播距離越長，流量保持穩(wěn)定時(shí)間延長。庫區(qū)蓄水深度增大使波速增大，導(dǎo)致流量穩(wěn)定時(shí)間縮短，庫區(qū)寬度對波速與傳播距離均沒有影響，因此穩(wěn)定時(shí)間基本不變。

3.5 各斷面峰現(xiàn)時(shí)間

圖7為沿程各斷面的峰現(xiàn)時(shí)間。各種工況下，隨著潰壩波的傳播，沿程各斷面峰現(xiàn)時(shí)間t2逐漸推遲。隨著庫區(qū)蓄水深度的增大，各斷面峰現(xiàn)時(shí)間均提前，隨著庫區(qū)長度的增大，各斷面峰現(xiàn)時(shí)間均推遲，庫區(qū)寬度的變化對各斷面峰現(xiàn)時(shí)間基本沒有影響。工況4、5的峰現(xiàn)時(shí)間明顯遲于工況6、7，即庫容相同時(shí)，庫區(qū)長度較大水庫各斷面的峰現(xiàn)時(shí)間較遲。

圖7 各斷面峰現(xiàn)時(shí)間Fig.7 Time to peak in each section

不同庫區(qū)長度工況下，同一時(shí)刻波前所在位置相同，潰壩波的附加比降i相同(圖4(c))。隨著庫區(qū)長度的增加，各斷面最大水深增大，即波峰高度增大(圖4(c))。潰壩波波前到波峰的水平距離增加(l=h/i)。庫區(qū)長度不同的工況下，沿程各斷面波前到達(dá)時(shí)間相同，峰現(xiàn)時(shí)間不同，主要原因是在波速v相等的情況下(圖6(c))，波峰到達(dá)時(shí)間將推遲(t=l/v)。

靠近壩址的斷面，由于洪水水頭大，傳播速度快，屬于高風(fēng)險(xiǎn)區(qū)。根據(jù)對各斷面洪峰流量的研究結(jié)果，庫區(qū)長度較大時(shí)，靠近壩址斷面的洪峰流量僅以較小的幅度增大。對于這些斷面，潰壩洪水的危害基本沒有加劇或小幅度地加劇，但峰現(xiàn)時(shí)間的推遲，使?jié)魏樗l(fā)生時(shí)，人民財(cái)產(chǎn)有了更充裕的轉(zhuǎn)移時(shí)間，這降低了潰壩洪水的危害，有利于潰壩風(fēng)險(xiǎn)的管理。

4 結(jié) 論

1) Flow-3D對潰壩洪水的模擬結(jié)果與理論研究相符合，與物理實(shí)驗(yàn)實(shí)測數(shù)據(jù)誤差較小，并且比人工推求洪水過程更加精細(xì)，能夠廣泛使用。

2) 增加庫區(qū)蓄水深度hu，壩址及沿程各斷面洪峰流量增大，峰現(xiàn)時(shí)間提前，潰壩洪水對下游的危害加劇，潰壩波波速增大。

3) 增加庫區(qū)寬度B，壩址及沿程斷面洪峰流量增大，各斷面峰現(xiàn)時(shí)間均不變，潰壩洪水對下游的危害不變，潰壩波波速基本不變。

4) 增加庫區(qū)長度L，壩址斷面洪峰流量不變，沿程各斷面洪峰流量以較小幅度增大，壩址斷面峰現(xiàn)時(shí)間不變，沿程各斷面峰現(xiàn)時(shí)間推遲，潰壩波波速不受影響。

綜合來看，通過增加庫區(qū)蓄水深度hu來增大庫容，會(huì)加劇潰壩洪水對下游的危害。庫區(qū)長度較大的水庫，有利于潰壩發(fā)生時(shí)下游人民財(cái)產(chǎn)的轉(zhuǎn)移，降低了潰壩洪水的危害。而庫區(qū)寬度基本不影響潰壩洪水對下游的危害。本研究對建造水庫的工程設(shè)計(jì)以及潰壩風(fēng)險(xiǎn)管理具有指導(dǎo)意義。