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

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

1 研究背景

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

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

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

2 潰壩洪水?dāng)?shù)值模擬

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

(1)

(2)

式中：xi為坐標(biāo)；ui為流體平均速度；t為時(shí)間；P為壓力；gi為重力加速度；v為水流運(yùn)動(dòng)黏度；vt為渦流運(yùn)動(dòng)黏度。

由于潰壩水流計(jì)算域的尺度較大，計(jì)算歷時(shí)較長(zhǎng)，兼顧計(jì)算精度與計(jì)算資源消耗，本研究采用RNGk-ε湍流模型。其湍動(dòng)能k和耗散率ε的輸運(yùn)方程可表述為

(3)

(4)

(5)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

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

上游水深H=0.4 m條件下不同水深比對(duì)各測(cè)點(diǎn)水位變化的影響如圖5所示，(a)、(b)、(c)分別為對(duì)應(yīng)下游水深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條件下四個(gè)測(cè)點(diǎn)水位隨時(shí)間演化規(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

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

3.2 沖擊動(dòng)力學(xué)分析本研究中設(shè)置18個(gè)測(cè)點(diǎn)用于捕捉下游壁面不同位置壓力演化過(guò)程，并對(duì)不同上游水深和水深比對(duì)壓力峰值的影響展開分析，測(cè)點(diǎn)位置布置如圖4(b)所示。

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

圖6 H=0.8，α=0.1、0.3、0.5條件下各壓強(qiáng)測(cè)量點(diǎn)相對(duì)壓強(qiáng)演化曲線Fig.6 Relative pressure evolution curves of each pressure measurement point for H=0.8，α=0.1，0.3，0.5

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

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

圖7 H=0.8 m時(shí)各點(diǎn)峰值壓力分布Fig.7 Peak pressure distribution at each point at H=0.8 m

當(dāng)y=37.5 mm或112.5 mm時(shí)，峰值變化過(guò)程基本重合，如圖7(a)所示。受邊壁的影響，兩側(cè)水流流速小于中間水流流速，導(dǎo)致y=75 mm時(shí)峰值壓力略大于y=37.5 mm(y=112.5 mm)的峰值壓力。但由于整個(gè)模擬中計(jì)算域?qū)挾葍H為15 cm，中間速度與兩側(cè)速度也相差較小，所以中間測(cè)點(diǎn)峰值壓力只是略高于兩側(cè)的峰值壓強(qiáng)。垂直方向的峰值壓強(qiáng)分布呈現(xiàn)先減小再變大再減小的趨勢(shì)，這是因?yàn)橄掠魏哟灿谐跏妓睿诔跏妓钪畠?nèi)，沿垂直方向，由于存在靜水壓強(qiáng)，峰值壓強(qiáng)分布是先減小。當(dāng)超過(guò)下游初始水深時(shí)，隨著z值的變大，峰值壓強(qiáng)分布先變大再變小，這是由于沖擊壓力以及垂直向上的射流引起的。圖7(b)展示了壓力峰值沿y軸方向呈對(duì)稱分布，z=0.3 m(P13、P14、P15)、z=0.6 m(P16、P17、P18)是在下游初始水深之上，z=0.3 m的壓力值主要是由沖擊壓強(qiáng)以及靜水壓強(qiáng)形成，而其余各點(diǎn)的壓強(qiáng)值主要是由靜水壓強(qiáng)值形成。圖7(c)則展示了α=0.3～0.5下的壓強(qiáng)峰值分布，隨著α變大，各點(diǎn)的壓強(qiáng)峰值在逐漸變小。峰值壓強(qiáng)在垂直方向和水平方向都是呈非線性分布。

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

4 結(jié)論

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

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

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

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