戎貴文，袁 岳，肖柏青，向英奇

(1．安徽理工大學(xué)地球與環(huán)境學(xué)院，安徽淮南 232001;2．水文水資源與水利工程科學(xué)國家重點實驗室，江蘇南京 210098)

病險水壩會對人民生命和財產(chǎn)安全構(gòu)成潛在的威脅，潰壩洪水多年來一直是水利工程師和學(xué)者們重點關(guān)注的問題。潰壩水流具有突發(fā)性和災(zāi)害性，當(dāng)潰壩發(fā)生時，準(zhǔn)確預(yù)測潰壩洪水的演進(jìn)過程對保護(hù)下游人民生命和財產(chǎn)安全具有重要意義。潰壩水流向下游演進(jìn)過程中如果遇到變化的斷面(如局部斷面收縮)，水流流態(tài)會發(fā)生較大的變化。河道中的堤壩、樹木、碎石以及其他一些建筑物都可以造成河道的局部斷面收縮。在這種情況下，大壩下游各處隨時間變化的水位、負(fù)波的形成與傳播、正波和負(fù)波的速度等都是影響準(zhǔn)確預(yù)測潰壩洪水演進(jìn)過程的重要參數(shù)。

對于潰壩研究來說，由于難以獲得隨時間變化的速度、水位、流量以及波前、波峰的移動速度等水動力要素的現(xiàn)場測量數(shù)據(jù)，因此，人們通常在實驗室中進(jìn)行模型試驗研究。Dressler［1］在水平矩形河道中對潰壩水流進(jìn)行了試驗研究。Bell等［2］研究了直線型河道和彎曲河道內(nèi)的二維潰壩水流。Bellos等［3］在干-濕河床上試驗研究了分汊河道內(nèi)的潰壩水流。Lauber等［4］在矩形光滑的干水槽內(nèi)進(jìn)行了潰壩水流的試驗研究，重點研究了潰壩水流正波和負(fù)波的形成和傳播。雖然這些試驗研究得到了一些關(guān)于潰壩水流的數(shù)據(jù)，但是由于試驗環(huán)境的假定、試驗?zāi)Ｐ偷拇笮?、測量裝置的影響等使得試驗結(jié)果的應(yīng)用有一定的局限性。隨著測量技術(shù)的提高以及計算機的高度發(fā)展，潰壩洪水研究結(jié)果變得更加準(zhǔn)確和直觀。Aureli等［5］采用試驗手段和MacCormack計算格式數(shù)值模擬了局部收縮河道內(nèi)的潰壩水流。戎貴文等［6］采用MacCormack預(yù)測-校正技術(shù)的隱式數(shù)值格式模擬了潰壩水流。符傳君等［7］建立了潰壩水流運動的三維數(shù)學(xué)模型，使用k-ε模型對水流方程進(jìn)行封閉，并對急彎河道進(jìn)行了數(shù)值模擬計算。Liang［8］提出了一種求解不規(guī)則地形下的二維淺水水流數(shù)值模擬技術(shù)，模擬了局部收縮河道內(nèi)的潰壩水流。王嘉松等［9］用TVD格式對瞬間全潰潰壩波的傳播、反射和繞射進(jìn)行了數(shù)值模擬。Chang等［10］采用光滑粒子流體動力學(xué)(SPH)方法求解淺水方程，研究了潰壩水流經(jīng)過有三角形障礙物河道時的流動情況。Kocaman等［11］利用新的攝像技術(shù)試驗研究了局部收縮對潰壩水流的影響。Lindsey等［12］利用試驗和大渦模型(LES)數(shù)值模擬的方法研究了二維潰壩水流，得到了潰壩水流的一些水位變化剖面圖。Jaan等［13］為研究潰壩水流遇到障礙物的運動過程，采用表面梯度逆風(fēng)方法(SGUM)求解淺水方程(SWEs)，并與采用不可壓縮光滑粒子流體動力學(xué)方法(ISPH)的模擬結(jié)果進(jìn)行了比較。王曉玲等［14］建立耦合流體體積函數(shù)(VOF)法的三維標(biāo)準(zhǔn)k-ε紊流模型，模擬了三維潰壩洪水在復(fù)雜淹沒區(qū)域的演進(jìn)過程，并繪制了潰壩洪水風(fēng)險圖。雖然國內(nèi)外對潰壩水流的研究已有很多，但是對三維局部斷面收縮地形內(nèi)潰壩水流流動的研究卻很少，而且對該種地形下負(fù)波的形成與傳播規(guī)律的研究更為罕見。

本文采用能夠更好地處理高應(yīng)變率及流線彎曲程度較大流動的RNG k-ε紊流模型［15-16］和具有傾斜網(wǎng)格校正能力、允許使用大時間步長、適合過渡流動計算的壓力隱式算子分割(PISO)算法［17］，研究了三維局部斷面收縮情況下潰壩水流的演進(jìn)過程，分析了負(fù)波的形成和傳播過程，得到了負(fù)波形成和傳播的規(guī)律以及潰壩水流水位變化特征。

1 數(shù)學(xué)模型及數(shù)值求解方法

1．1 控制方程

式中:ρ為流體的密度，kg/m3;t為時間，s;xi、xj為坐標(biāo)分量(i，j=1，2，3)，m;ui、uj為速度分量(i，j=1，2，3)，m/s;p為修正壓力，Pa;μ為分子動力黏性系數(shù)，N·m/s;k為紊動能，m2/s2;ε為紊動耗散率，m2/s2;σk、σε分別為紊動能k和紊動耗散率ε對應(yīng)的普朗特數(shù);μt為紊動黏度;G 為紊動能生成項;cμ、C1ε、C2ε為經(jīng)驗常數(shù)，C1ε的表達(dá)式如下:

式中:Sij為平均應(yīng)力張量。根據(jù)Yakhot等［15］的研究結(jié)果，σk、σε、cμ、C2ε的取值分別為 0. 7179、0. 7179、0.0845和1.68。

1．2 數(shù)值求解方法

采用有限體積法離散控制方程，對控制方程中的時間項使用隱式歐拉格式離散，對Navier-Stokes方程中的對流項采用TVD格式離散［18］，與此同時使用一階迎風(fēng)格式離散紊流模型控制方程中的對流項，其余各項采用二階中心差分進(jìn)行離散。物理變量采用交錯網(wǎng)格布置，在每個控制單元中對微分方程進(jìn)行積分，再將積分方程線性化，把控制方程離散為可以數(shù)值求解的代數(shù)方程，采用PISO算法對方程組進(jìn)行求解［17］。該算法由 Issa［19］于 1986 年提出，適用于非穩(wěn)態(tài)問題。PISO算法包括1個預(yù)測步驟和2個校正步驟，可以看作是在SIMPLE算法的基礎(chǔ)上增加了1個校正步驟，通過一系列的預(yù)測-校正來求解每一時刻的流動物理量［20］。

1．3 自由表面的處理

潰壩水流自由表面采用VOF方法處理［21］。與其他處理自由表面的方法，如剛蓋假定法、MAC法和LINC法相比，VOF方法用體積分?jǐn)?shù)代替了標(biāo)記點，只用一個函數(shù)就可以描述自由表面的各種復(fù)雜變化，大幅度節(jié)約了計算機內(nèi)存和計算時間，而且計算精度也較高［16］。VOF方法的基本原理是通過研究網(wǎng)格單元中第q相流體體積和網(wǎng)格體積比函數(shù)αq來確定自由面，追蹤流體的變化，而非追蹤自由液面上質(zhì)點的運動。若αq=1，說明該單元全部為第q相流體所占據(jù);若αq=0，則該單元中沒有第q相流體;若0＜αq＜1，表明該單元包含了第q相流體和一相或者其他多相流體的界面。ρ和μ由體積分?jǐn)?shù)加權(quán)平均得出，計算公式如下:

式中:下標(biāo)w和a分別代表水和氣體;αw為水的體積分?jǐn)?shù)。

2 數(shù)值模擬

2．1 物理模型

本文所采用物理模型是Kocaman等［11］的水槽模型(圖1)，假定模型表面光滑。該模型整體上是一個六面體，長8.9 m，寬0.30 m，高0.34 m。用一塊擋板將模型分為上下游兩個部分，上游為有水區(qū)域，代表水庫庫區(qū)，下游為無水區(qū)域。初始時刻水槽上游水位為0.25 m，下游無水。在下游1.52 m處河道兩側(cè)有一個對稱的三棱柱體，形成了一個局部斷面收縮的區(qū)域，最窄處只有0.1 m。

圖1 水槽模型示意圖 (單位:cm)

2．2 網(wǎng)格劃分及邊界條件

2．2．1 網(wǎng)格劃分

在笛卡兒坐標(biāo)系下對上述物理模型進(jìn)行網(wǎng)格剖分，x、y、z分別代表物理模型的長、高、寬方向，由于坐標(biāo)原點的選擇對文中的研究內(nèi)容沒有影響，理論上可以選擇在任意位置，因此將坐標(biāo)原點(0，0，0)選在擋板底部中心處。由于該模型是對稱的，對其進(jìn)行數(shù)值模擬計算時只需取計算模型的一半，所以在笛卡兒坐標(biāo)系下只對模型的一半進(jìn)行網(wǎng)格劃分，劃分的網(wǎng)格為均勻結(jié)構(gòu)性網(wǎng)格，每個網(wǎng)格的大小為Δx=Δy=Δz=8 mm。

2．2．2 邊界條件

由于模型是對稱的，將模型對稱面的邊界條件設(shè)為鏡面對稱;模型的下游出口邊界條件設(shè)為自由出流;由于上游沒有來水條件，將上游邊界設(shè)為壁面邊界;該試驗是在常規(guī)大氣壓下進(jìn)行的，故將模型的頂部邊界條件設(shè)為鏡面對稱;其他壁面的邊界都設(shè)為壁面邊界。

2．3 模型驗證

在初始時刻，上游長4.65m、寬0.15m、高0.25m的區(qū)域內(nèi)充滿水，其他區(qū)域無水。在模型軸線處設(shè)置了6個測點，測點位置見圖2。試驗?zāi)M了這6個測點處水位隨時間的變化規(guī)律，通過與Kocaman等［11］的試驗數(shù)據(jù)的比較，驗證文中所采用模型的可行性。

圖2 測量水位的位置點 (單位:cm)

圖3顯示了6個測點處水位隨時間變化的模擬值與測量值，圖3表明該模型可以準(zhǔn)確地模擬局部斷面收縮地形內(nèi)的潰壩水流，驗證了該模型的可行性。

2．4 模擬結(jié)果

2．4．1 負(fù)波的形成與傳播

將擋板突然移去后，模擬的潰壩水流迅速流向下游。當(dāng)潰壩水流到達(dá)局部收縮的斷面時，一部分水流經(jīng)過該收縮斷面后繼續(xù)向下游流動;另一部分水流擁堵在局部收縮段，抬高了這部分的水位，形成了一個波峰。當(dāng)波峰達(dá)到一定高度時就會向上游移動形成負(fù)波，從圖4可以清楚地看出形成負(fù)波的過程。隨著擋板的突然移去，水流在t=0.9 s時到達(dá)收縮段，部分水流通過收縮段繼續(xù)流向下游。在t=0～1.2s時，水流漸漸擁堵在收縮段中，水位不斷升高。從t=1.5s開始收縮段上游形成了強烈的水躍，從收縮段到突擴段水位迅速降低。在t=2.1 s時，水流到達(dá)收縮段的最窄斷面，此時水位迅速上升，部分水流反向上游流動。在t=3.0 s之前，收縮段的水位一直在升高，自由液面的形狀為凸形。與此同時，流向突擴段的水流自由液面呈現(xiàn)出凹形。

圖3 各測點處水位隨時間變化的測量值與模擬值

圖4 負(fù)波的形成與傳播過程

潰壩水流由于通過局部收縮的斷面，形成了負(fù)波。按時間順序?qū)⒇?fù)波的變化過程分為兩個階段，即負(fù)波的形成階段和負(fù)波的傳播階段。從圖4可以看出，t=0.9～3 s為負(fù)波的形成階段，t=3～5 s為負(fù)波的傳播階段。筆者從能量的轉(zhuǎn)化角度來解釋負(fù)波的形成和傳播過程。由于收縮斷面的存在，當(dāng)水流的動能不足以使水流通過斷面時，水流就會擁堵在收縮口處將動能轉(zhuǎn)化為勢能。這個能量轉(zhuǎn)化的過程導(dǎo)致了收縮口上游的水位逐漸升高，水流的流態(tài)也相應(yīng)地發(fā)生變化，從上游的急流變?yōu)榫徚鳌８鶕?jù)水力學(xué)的知識，在自然條件下當(dāng)水流流態(tài)由急流轉(zhuǎn)變?yōu)榫徚鲿r必然發(fā)生水躍現(xiàn)象，這就很好地解釋了從t=1.5 s起開始形成水躍的原因。

圖5 測點處水流速度隨時間的變化過程

圖5為測點處水流速度、負(fù)波速度隨時間的變化過程，圖5(a)反映了潰壩水流向下游傳播時，上述6個測點處流速的變化情況，圖5(b)反映了局部收縮斷面處負(fù)波速度隨時間的變化情況。由于取潰壩水流正常流動方向為正方向，所以向上游傳播的負(fù)波速度為負(fù)值，圖5(b)中以－ν波表示負(fù)波速度。從圖5(a)可以看出，在t=1～3 s時，收縮斷面處的測點P4、P5、P6處向下游流動的速度大幅度減小，這主要是由于上游大量水流到達(dá)最窄斷面后水流不能及時流出，擁堵在收縮斷面附近。圖5(b)中，在t=2～3 s時，隨著負(fù)波的移動，P4、P5、P6 處的速度增大，反映了負(fù)波的波前速度隨著負(fù)波的移動逐漸增加。在t=3～5 s時，向下游流動的潰壩水流速度與向上游移動的負(fù)波速度逐漸保持不變，此時負(fù)波正向上游傳播，并與向下游傳播的水流處于能量平衡狀態(tài)，這與圖4中3～5 s時負(fù)波的自由液面漸漸呈現(xiàn)水平狀態(tài)保持一致。

2．4．2 隨時間變化的潰壩水流水位

圖6為各測點處潰壩水流的水位變化過程。從圖6(a)可以看出，隨著擋板的突然移去，P1處的水位迅速降低，而隨著水流的運動，P2處水位隨后迅速升高，直到t=1.6s左右。在t=1.6～5.0s之間，P2處的水位基本保持恒定，水位約為0.116 m。Ritter［22］研究得出此處的水位 h=4/9h0(h0為來流水位)，P2處模擬得到的水位與Ritter的研究結(jié)果吻合。當(dāng)t=5.0 s時，收縮段產(chǎn)生的負(fù)波影響了P2處的水位，P2處的水位隨之迅速升高。同樣的情形也可以從P1、P3處的水位變化看出。P1、P2、P3處的水位隨時間的變化，間接地反映了負(fù)波的傳播過程。從圖6(b)中可以看出，t=0.9 s左右時水流到達(dá)收縮段的最窄部分，P4、P5、P6處的水位隨后迅速上升。在t=3 s時P5處的水位達(dá)到最大值，約為0.20m，隨后逐漸降低。圖6(b)恰好說明了t=0～3 s為負(fù)波形成的時間，t=3 s后負(fù)波開始傳播，P4、P5、P6處的水位保持在一定高度，直至水流越過收縮段最窄處，隨后這3個測點處的水位逐漸降低。

圖6 各測點處潰壩水流的水位變化過程

3 結(jié)論

本文利用三維RNG k-ε紊流模型對局部斷面收縮情況下的三維潰壩水流進(jìn)行了數(shù)值模擬，給出了水槽內(nèi)負(fù)波的形成與傳播過程、潰壩水流在測點處的水位和速度變化，分析了局部收縮對潰壩水流水力特性的影響。通過與已有試驗數(shù)據(jù)的比較，驗證了該模型的可行性，表明該模型可以準(zhǔn)確地模擬局部斷面收縮地形內(nèi)的潰壩水流。采用該紊流模型得到了三維局部斷面收縮水槽內(nèi)潰壩水流的水位變化過程以及負(fù)波的形成和傳播規(guī)律，可為研究實際情況下的潰壩水流提供參考。

但本文在研究潰壩水流時，假設(shè)大壩的潰壩形式為瞬間全潰，下游河道壁面光滑，同時潰壩水流中只含有氣液兩相。與全潰壩相比，部分潰壩液面變化復(fù)雜，并且上游形成的波形也不同;在相同的時間內(nèi)潰口的大小對潰壩波推進(jìn)的距離有影響，特別是潰壩初期潰口的大小對推進(jìn)距離影響較大。實際的潰壩過程中，大壩并非瞬間全潰，而是部分逐漸潰決，而且由于水流的沖刷作用，水體中含有泥沙以及其他雜物，所以該模型還存在一定的局限性。為了更準(zhǔn)確地模擬真實的潰壩水流，可以在模擬時加入泥沙，假設(shè)壁面有一定的粗糙度等，也可以將計算區(qū)域內(nèi)的網(wǎng)格剖分得更精細(xì)，模擬的準(zhǔn)確度將會更高，但計算時間和計算成本也會相應(yīng)增加。

［1］DRESSLER R F．Comparison of theories and experiments for the hydraulic dam-break wave［J］．Int Ass Sci Hydrol Publi，1954，38(3):319-328．

［2］BELL S W，ELLIOT R C，CHAUDHRY M H．Experimental results of two dimensional dam-break flows［J］．Journal of Hydraulic Research，1992，30(2):225-252．

［3］BELLOS C V，SOULIS J V，SAKKAS J G．Experimental investigation of two dimensional dam-break induced flows［J］．Journal of Hydraulic Research，1992，30(1):47-63．

［4］LAUBER G，HAGER W H． Experiments to dam-break wave: horizontal channel ［J］． Journal of Hydraulic Research， 1998， 36 ( 3) : 291-308．

［5］AURELI F，MIGNOSA P，TOMIROTTI M．Numerical simulation and experimental verification of dam-break flows with shock［J］．Journal of Hydraulic Research，2000，38(3):197-206．

［6］戎貴文，魏文禮，嚴(yán)建軍．潰壩洪水演進(jìn)數(shù)值模擬研究［J］．黑龍江水專學(xué)報，2008，35(3):12-14．(RONG Guiwen，WEI Wenli，YAN Jianjun．Numerical simulation of 2D dam-break flood waves［J］．Journal of Heilongjiang HydraulicEngineering，2008，35(3):12-14．(in Chinese))

［7］符傳君，練繼建． 潰壩水流在復(fù)雜河道中傳播的三維數(shù) 值模擬［J］． 水利學(xué)報，2007，38 ( 10) : 1151-1157． ( FU Chuanjun，LIAN Jijian． Three-dimensional numerical simulation of dam-break flow in complicated river sections［J］． Journal of Hydraulic Engineering，2007，38 ( 10 ) :1151-1157． ( in Chinese) )

［8］LIANG Qiuhua．Simulation of shallow flows in non-uniform open channels［J］．Journal of Fluid Mechanics，ASME，2008，130(1):51-59．

［9］王嘉松，倪漢根，金生．瞬間全潰潰壩波的傳播、反射和繞射的數(shù)值模擬［J］．水動力學(xué)研究與進(jìn)展:A輯，2000，15(1):1-7．(WANG Jiasong，NI Hangen，JIN Sheng．Simulation of propagation，reflection and diffraction of bores caused by sudden and full destruction of a dam［J］．Journal of Hydrodynamics:Ser A，2000，15(1):1-7．(in Chinese))

［10］CHANG Tsangjung，KAO Hongming，CHANG Kaohua，et al．Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydrodynamics［J］．Journal of Hydrology，2011，408(1):78-90．

［11］KOCAMAN S，OZMEN C H．The effect of lateral channel contraction on dam break flows:laboratory experiment［J］．Journal of Hydrology，2012(432/433):145-153．

［12］LINDSEY A L，JASIM I M，HANIF C．Experimental and numerical investigations of two-dimensional dam-break flows［J］．Journal of Hydraulic Engineering，2012，139(6):569-579．

［13］JAAN H P，SHAO Songdong，HUANG Yuefei，et al． Evaluations of SWEs and SPH numerical modelling techniques for dam break flows ［J ］． Engineering Applications of Computational Fluid Mechanics，2013，7( 4) : 544-563．

［14］王曉玲，張愛麗，陳華鴻，等．三維潰壩洪水在復(fù)雜淹沒區(qū)域演進(jìn)的數(shù)值模擬［J］．水利學(xué)報，2012，43(9):1025-1033．(WANG Xiaolin，ZHANG Aili，CHEN Huahong，et al．Three-dimensional numerical simulation of dam-break flood routing in the complex inundation areas［J］．Journal of Hydraulic Engineering，2012，43(9):1025-1033．(in Chinese))．

［15］YAKHOT V，ORZAG S A．Renormalization group analysis of turbulence:basic theory［J］．J Scient Comput，1986，1(1):3-51．

［16］戎貴文，魏文禮，劉玉玲，等．涌潮作用下丁壩附近水流運動特性的數(shù)值模擬研究［J］．水利學(xué)報，2012，43(3):296-301．(RONG Guiwen，WEI Wenli，LIU Yulin，et al．Study on flow characteristics near spur dikes under tidal bore［J］．Journal of Hydraulic Engineering，2012，43(3):296-301．(in Chinese))

［17］魏文禮，王德意．計算水力學(xué)理論及應(yīng)用［M］．西安:陜西科學(xué)技術(shù)出版社，2001．

［18］劉玉玲，王玲玲，周孝德，等．二維潰壩洪水波傳播的高精度數(shù)值模擬［J］．自然災(zāi)害學(xué)報，2010，19(5):164-169．(LIU Yuling，WANG Lingling，ZHOU Xiaode，et al．High-resolution numerical simulation of 2D dam-break flood waves［J］．Journal of Natural Disasters，2010，19(5):164-169．(in Chinese))

［19］ISSA R I．Solution of the implicitly discretized fluid flow equations by operator-splitting［J］． Journal of Computational Physics，1986，62:40-65．

［20］李人憲．有限體積法基礎(chǔ)［M］．北京:國防工業(yè)出版社，2008．

［21］HIRT C W，NICHOLS B D．Volume of fluid(VOF)method for the dynamics of free boundary［J］．Journal of Computational Physics，1981，39:201-225．

［22］RITTER A．The propagation of water waves［J］．Journal of German Engineers Association，1892，36(33):947-954．