周昔東，何小瀧，袁 浩，孫 倩

(1.重慶交通大學 河海學院，重慶 400074；2.重慶交通大學 西南水運工程科學研究所，重慶 400016；3.重慶交通大學 土木工程學院，重慶 400074)

1 研究背景

在氣候變化大環(huán)境下，全球氣溫逐漸升高[1]，使得我國極端降雨頻發(fā)，集中高頻率降雨極可能導致水庫潰決。潰壩導致水庫庫容在短時間內下泄，引發(fā)的洪水對下游建筑、河道河岸、下游大壩壩體帶來巨大沖擊，對于存在梯級水庫的河道，極可能引發(fā)連鎖反應，導致梯級潰壩。潰壩洪水不僅會嚴重破壞下游河道[2]以及結構物，而且會對下游的生態(tài)環(huán)境造成負面影響[3]。

潰壩實驗可以直觀展示出潰壩波的演進、壁面沖擊壓力的演化等諸多規(guī)律。研究者針對潰壩水流的演進與對建筑物的沖擊開展了系統(tǒng)實驗研究，研究內容主要集中在建筑物附近局部流場、水位的變化[4-5]，干濕河床對潰壩演進的影響[6-8]，下游結構物受到的沖擊及壓力演化過程[9-12]。受限于實驗尺度和觀測手段，實驗研究并不能完全反饋出真實潰壩過程中潰壩波演進過程中的坦化和流速變化及下游建筑物所受到的真實沖擊壓力。而數(shù)值模擬則彌補了上述缺陷，可以獲得潰壩過程中的自由液面的破碎和流場的演化過程等更為詳細的流場信息。研究表明由于潰壩水流運動速度快，且在結構附近的自由水面有明顯的變化，垂向壓力對下游結構的影響不僅有靜水壓力還有沖擊壓力[13]。并進一步獲得了不同上游庫區(qū)水位條件下潰壩洪水對下游結構物附近水流流態(tài)和所受沖擊壓力的影響[14]。

目前針對潰壩水流實驗和數(shù)值模擬研究主要集中在下游無障礙性潰壩水流演進和對下游有結構物阻擋的潰壩水流演進與沖擊這兩方面。由于潰壩洪水作用在結構上的沖擊脈動壓力的持續(xù)時間較短，這種沖擊壓力在許多關于潰壩水流研究中可能被低估或甚至被忽視[15]，此沖擊壓力雖持續(xù)時間短，但破壞性大。由于當前研究中極少考慮潰壩波與潰壩下游壩體相互作用的動態(tài)沖擊過程，本文基于RNGk-ε湍流模型對三維潰壩洪水的演進及潰壩波與下游壩體的相互作用開展系統(tǒng)數(shù)值研究，旨在揭示在上下游不同的初始水深比α(H/h，H為上游水深，h為下游水深)條件下，下游壩體壁面受到的沖擊壓強峰值的變化規(guī)律及不同下游水深對潰壩洪水演進影響。

2 潰壩洪水數(shù)值模擬

2.1 計算方法簡介本研究基于FLOW-3D商業(yè)軟件進行數(shù)值模擬，F(xiàn)LOW-3D采用有限差分法對控制方程進行離散，且此次計算采用隱式方法對壓力速度進行求解，其中采用的廣義最小殘差(GMRES)迭代法在求解大型代數(shù)方程時具有良好的收斂性、高效的計算效率和計算精度。潰壩水流演進過程可通過Navier-Stokes(N-S)方程開展模擬，其控制方程包括連續(xù)方程和動量方程：

(1)

(2)

式中：xi為坐標；ui為流體平均速度；t為時間；P為壓力；gi為重力加速度；v為水流運動黏度；vt為渦流運動黏度。

由于潰壩水流計算域的尺度較大，計算歷時較長，兼顧計算精度與計算資源消耗，本研究采用RNGk-ε湍流模型。其湍動能k和耗散率ε的輸運方程可表述為

(3)

(4)

(5)

式中：Rε為剪切湍流附加項；G為由于平均速度梯度而產(chǎn)生湍流動能的產(chǎn)生率。本研究中湍流模型中各個系數(shù)取值為：Cμ=0.085，C1ε=1.42，C2ε=1.68，σε=σk=0.7194，β=0.012，η0=4.38[14]。

同時本研究中使用流體體積法(VOF)捕獲自由液面，即在整個計算域內追蹤每一計算單元的體積分數(shù)對于每個網(wǎng)格單元，函數(shù)F(x1，x2，x3，t)定義為流體所占體積與整個單元的體積之比，當F值為1時表示單元格體積完全被水占據(jù)，F(xiàn)值為0時表示單元格體積完全被空氣占據(jù)。F值在0和1之間的單元表示對應于水氣交接面。

2.2 模型驗證及網(wǎng)格質量分析圖1中L為計算域長度，D為水箱高度，H為水庫初始水位高度，H1、H2、H3、H4為水位測量位置，圖中尺寸單位均為毫米。

圖1 數(shù)值模擬計算域示意圖Fig.1 Schematic diagram of the calculation domain of numerical simulation

模擬結果使用h/H、P/(ρgH)、t/(H/g)1/2對水位和壓強以及洪水演進時間進行無量綱化，以便后續(xù)系統(tǒng)地開展數(shù)值結果與實驗數(shù)據(jù)進行定量對比，并采用均方根誤差(RMSE)來定量評估水位數(shù)值模擬和實驗測量之間的差異。由于本文重點關注潰壩水流對下游垂直壩體沖擊的壓強峰值，所以對數(shù)值結果和實驗數(shù)據(jù)壓強之間的差異采用峰值誤差進行定量分析。

不同網(wǎng)格的數(shù)值結果與給定位置H1、H2、H3和H4實驗數(shù)據(jù)水位的對比如圖2所示，數(shù)值模型對波浪演化的特征提供了合理的預測，水位數(shù)值結果在二次波到來前各位置的實驗數(shù)據(jù)基本一致，較準確的預測了二次波的到達時間及其在H1，H2，H3和H4位置的水位高度。受限于網(wǎng)格尺度，本研究中無法準確獲得潰壩波沖擊引起的自由表面破碎所帶來的水位波動，導致數(shù)值模擬對二次波后的水位預測值與實驗值差異較大，預測精度降低。水位對應的實驗數(shù)據(jù)對比如表1所示，誤差最大值出現(xiàn)在受沖擊射流引起的強回流波影響最大的H4 位置，其最大誤差小于14%。H2和H4位置的相對誤差相對較大，造成H2位置水位誤差較大其原因可能為水庫水體突然下泄，下泄的水體不穩(wěn)定以及閘門提升過程對水體產(chǎn)生影響；而造成H4位置水位誤差較大其原因可能是測點靠近下游壩體墻壁，潰壩洪水到達壩體會形成壅水回流及表面破碎，引發(fā)水位振蕩。

圖2 數(shù)值模擬和實驗的水位比較Fig.2 Comparison between numerical and experiment water level

表1 不同網(wǎng)格尺寸下水位均方根誤差

圖3展示了各點P1、P2、P3、P4的數(shù)值和實驗壓強值的對比。其中圖3(a)展示了100次實驗測試的數(shù)值壓強預測和壓強時間歷程的置信區(qū)間，分別以97.5%和2.5%百分位水平為上限和下限[12]。數(shù)值模擬得到的壓強隨時間變化的過程與實驗中所獲一致，較好的反演了實驗獲得的P1、P2和P3處的峰值壓力。其中P1、P2點壓力隨時間變化曲線與實驗數(shù)據(jù)一致。受潰壩波和自由表面破碎的影響，對P3、P4點數(shù)值模擬獲得的最大壓力略小于實驗值。表1展示了不同網(wǎng)格尺度下各測點的最大壓力誤差，在5 mm網(wǎng)格尺寸下，最大壓強誤差在P3點，兩者誤差不超過18%；在8 mm網(wǎng)格尺寸下，最大壓力誤差在P1點，兩者誤差不超過15%；而在10 mm網(wǎng)格尺寸下，最大壓力誤差同樣在P1點，兩者誤差不超過25%。同時數(shù)值模擬的壓強峰值比實驗的壓力峰值略微滯后，其可能與數(shù)值模擬糙率與試驗模型有差異和實驗隨機誤差相關。通過網(wǎng)格密度分析表明，平均網(wǎng)格尺度為8 mm，可滿足計算精度要求，同時可減小計算資源耗費。

圖3 不同網(wǎng)格密度下壓力值與實驗值對比Fig.3 Comparison of pressure values with experimental values at different grid sizes

表2為水槽H1位置水位在不同時刻下本研究與文獻[16]中數(shù)值模擬結果與實驗結果對比獲得的水位相對誤差。結果表明本研究與文獻[16]相比其相對誤差略大，其原因可能是由于本研究采用8 mm的網(wǎng)格尺度略大于文獻[16]采用4 mm的網(wǎng)格尺度，同時本研究未涉及閘門的開啟過程對潰壩洪水演進的影響。因此在后續(xù)模擬研究中將進一步開展網(wǎng)格尺度的研究，并結合考慮閘門開啟過程對潰壩水動力過程模擬的影響。

表2 H1位置水位在不同時刻本文與文獻[16]的數(shù)值模擬結果相對誤差對比

2.3 計算模型與數(shù)值模擬工況本研究中計算區(qū)域及概化示意圖如圖4(a)所示。計算域為1.61 m×0.15 m×1.5 m立方體，該計算域初始分為兩部分，包括閘門上游0.6 m長的水庫和下游1.01 m長的濕河床。在下游河床的尾部采用固壁邊界用以模擬下游建筑物。計算域其余初始和邊界條件與第2節(jié)設置相同。水庫上游水深為H，下游濕河床水深為h。為了分析這些上下游水深對潰壩水流演進以及沖擊壓力的影響，在數(shù)值模擬中考慮了3個初始上游水深和5個水深比(定義為α=h/H，其中h初始下游水深)，各模擬工況如表3所示。

圖4 數(shù)值模擬的幾何模型以及下游墻壁測點布置圖Fig.4 Geometric model of numerical simulation and downstream wall measurement point arrangement

表3 數(shù)值模擬的初始水深條件

3 結果與分析

3.1 運動學分析本研究主要考慮上下游不同水深對潰壩波演進和下游壩面沖擊的影響。選取x1=855 mm、x2=1105 mm、x3=1355 mm、x4=1605 mm斷面研究潰壩波的波高演進過程及潰壩波演進的平均速度u。

上游水深H=0.4 m條件下不同水深比對各測點水位變化的影響如圖5所示，(a)、(b)、(c)分別為對應下游水深0.04 m、0.12 m、0.2 m，且α分別為0.1、0.3、0.5。

圖5 H=0.4 m，α=0.1，0.3，0.5條件下四個測點水位隨時間演化規(guī)律Fig.5 Evolution rule of water level with time at four measurement points under H=0.4 m，α=0.1，0.3，0.5

當下游水深為0.04 m，對應水深比α=0.1時，各測點水深隨時間變化如圖5(a)所示，x1位置最大波高為0.157 m，水位峰值到達時間為0.219 s；x2和x3位置處最大波高分別為0.188 m和0.155 m，最大水位峰值對應時間則分別為0.360 s和0.470 s；x4位置最大波高為0.789 m，最大峰值時間為0.839 s。潰壩波導致下游x1至x3區(qū)間內水面線峰值呈現(xiàn)一個先增加隨后坦化的趨勢。而受到下游壩體壁面的影響，x4位置潰壩波接觸壁面后發(fā)生壅水，波高迅速增加。而對比圖5不同水深條件下，α=0.3條件下，其相對峰值(相對峰值為該位置水位峰值減去初始水深)達到最大值。隨著水深比增大，其潰壩波傳播速度和各點水位峰值均減小。

3.2 沖擊動力學分析本研究中設置18個測點用于捕捉下游壁面不同位置壓力演化過程，并對不同上游水深和水深比對壓力峰值的影響展開分析，測點位置布置如圖4(b)所示。

3.2.1 壓強變化 圖6展示了上游水深H=0.8 m，不同上下游水深比α條件下各壓力點捕捉到的壓強變化過程。在相同α時較低和較高位置的壓強隨時間變化是明顯不同的。根據(jù)壓強的主要形成原因，較低位置測點的壓強可分為兩部分：①主要的潰壩洪水波沖擊下游垂直壩體墻壁時，由于動能損失而產(chǎn)生的瞬時動壓，其大小隨沖擊速度增加而增加；②洪水波沖擊后的靜水壓強，其大小隨垂直壩體墻壁前方水深的增加而增加。

圖6 H=0.8，α=0.1、0.3、0.5條件下各壓強測量點相對壓強演化曲線Fig.6 Relative pressure evolution curves of each pressure measurement point for H=0.8，α=0.1，0.3，0.5

下游水深對洪水波的演進是有阻礙作用且可以緩解洪水波與下游壩體的沖擊。當α=0.1時瞬時沖擊壓強遠大于其他任何時間階段的壓強，可達到第一次沖擊后最大壓強的7.5～22.1倍；α=0.3時瞬時沖擊壓強是第一次沖擊后最大壓力的1.1～1.3倍；α=0.5時瞬時沖擊壓強與第一次沖擊后最大壓強基本一致?？梢缘贸鱿掠螇Ρ谏喜糠譁y點的壓強，主要是靜水壓強，類似于下部分測點的第二種壓強。

3.2.2 峰值壓強 圖7(a)、圖7(b)是H=0.8 m、α=0.3(h=240 mm)時沿z軸和y軸的峰值壓力變化，圖7(c)是H=0.8 m、α=0.3～0.5、y=75 mm時沿z軸(垂直方向)的峰值壓強變化。

圖7 H=0.8 m時各點峰值壓力分布Fig.7 Peak pressure distribution at each point at H=0.8 m

當y=37.5 mm或112.5 mm時，峰值變化過程基本重合，如圖7(a)所示。受邊壁的影響，兩側水流流速小于中間水流流速，導致y=75 mm時峰值壓力略大于y=37.5 mm(y=112.5 mm)的峰值壓力。但由于整個模擬中計算域寬度僅為15 cm，中間速度與兩側速度也相差較小，所以中間測點峰值壓力只是略高于兩側的峰值壓強。垂直方向的峰值壓強分布呈現(xiàn)先減小再變大再減小的趨勢，這是因為下游河床有初始水深，在初始水深之內，沿垂直方向，由于存在靜水壓強，峰值壓強分布是先減小。當超過下游初始水深時，隨著z值的變大，峰值壓強分布先變大再變小，這是由于沖擊壓力以及垂直向上的射流引起的。圖7(b)展示了壓力峰值沿y軸方向呈對稱分布，z=0.3 m(P13、P14、P15)、z=0.6 m(P16、P17、P18)是在下游初始水深之上，z=0.3 m的壓力值主要是由沖擊壓強以及靜水壓強形成，而其余各點的壓強值主要是由靜水壓強值形成。圖7(c)則展示了α=0.3～0.5下的壓強峰值分布，隨著α變大，各點的壓強峰值在逐漸變小。峰值壓強在垂直方向和水平方向都是呈非線性分布。

圖8 α=0.3時不同上游水深在不同水平距離上的沿垂直方向峰值壓強分布Fig.8 Distribution of peak pressure along the vertical direction at different upstream water depths at different horizontal distances at α=0.3

圖8則展示了在α=0.3時不同上游水深下，峰值壓強在不同水平(y軸)距離上的沿垂直方向分布。對H=0.4 m，h=0.12 m分析，z軸方向0.12 m以下位置(線H=0.4 m，h=0.12 m中空心點)的壓強主要是由靜水壓強影響，在0.12 m以上的第一個位置(線H=0.4 m，h=0.12 m中實心點)主要受沖擊壓強以及靜水壓強影響，會出現(xiàn)瞬間壓強變大，最大峰值壓強為12.61 kPa，隨后水流受壁面阻擋，發(fā)生水位壅高，以上點的壓強主要是靜水壓強形成，是逐漸減??；工況H=0.6 m，h=0.18 m、H=0.8 m，h=0.24 m與工況H=0.4 m，h=0.12 m規(guī)律類似，最大沖擊壓強存在下游初始水深之上第一個點，最大峰值壓強分別為22.5 kPa、19.4 kPa。圖8(b)和圖8(c)與圖8(a)有類似的規(guī)律，在此就不做贅述。結果表明最大峰值壓強主要出現(xiàn)在下游水深之上，由潰壩波引起的沖擊壓強產(chǎn)生。

4 結論

本文采用RNGk-ε湍流模型對潰壩誘發(fā)洪水演進及其對下游垂直壁面沖擊的影響開展了研究。著眼于上游水深和上下游水深比對洪水演進的影響，得到如下結論：

(1)潰壩波導致下游x1至x3測點區(qū)間內水面線峰值呈現(xiàn)一個先增加隨后坦化的趨勢。而受到下游壁面的影響，x4位置潰壩波接觸壁面后發(fā)生壅水，波高迅速增加。α≥0.3條件下，其相對峰值達到最大值，隨著水深比增大，其潰壩波傳播速度和各點水位峰值均減小，表明下游水深對洪水波的演進有阻礙作用。

(2)α=0.1時瞬時沖擊壓強遠大于其他任何時間階段的壓強，可達到第一次沖擊后最大壓強的7.5～22.1倍。隨著α的增大，瞬時沖擊壓強值減小且與第一次沖擊后最大壓強值的倍數(shù)減小，可以得出下游墻壁上部分測點的壓強，主要是靜水壓強，類似于下部分測點的第二種壓強。

(3)下游垂直墻壁的峰值壓強在垂直和水平方向上的分布都是非線性的，垂直方向的峰值壓強分布呈現(xiàn)先減小再變大再減小的趨勢，而水平方向呈對稱分布，且最大沖擊壓強出現(xiàn)在下游水深之上的第一個點，由沖擊壓強產(chǎn)生。