王興華，付成華，穆宵梟

(西華大學(xué)能源與動(dòng)力工程學(xué)院，成都 610039)

2018年7月23日晚，老撾東南部桑片-桑南內(nèi)水電站副壩發(fā)生潰壩；2019年1月25日，巴西東南部米納斯吉拉斯州布魯馬迪紐市發(fā)生潰壩事故，均造成重大的生命財(cái)產(chǎn)損失及環(huán)境破壞。潰壩問(wèn)題作為一種低頻率、高危害的社會(huì)致災(zāi)因素[1]，越來(lái)越受到人們的關(guān)注，潰壩水流數(shù)值模型也被廣泛研究。為了模擬潰壩水流，開(kāi)發(fā)了大量的一維、二維模型[2]，這些模型適用于確定潰壩水流的淹沒(méi)深度和到達(dá)時(shí)間，但由于對(duì)控制方程做了假設(shè)，例如靜水壓力分布、不明顯的垂直加速度和自由表面曲率等[3]，不能準(zhǔn)確地表現(xiàn)水流在空間上的演進(jìn)變化規(guī)律，不適應(yīng)于計(jì)算潰壩初期、潰壩波前和潰壩結(jié)構(gòu)物附近水流，模擬這些復(fù)雜流動(dòng)應(yīng)基于使用雷諾平均N-S方程的三維模型。

基于雷諾平均N-S方程的三維潰壩模型，主要耦合紊流模型和自由界面捕獲法，利用有限體積法求解控制方程。紊流模型主要有一方程模型、兩方程模型(標(biāo)準(zhǔn)k-ε、RNGk-ε、Realizablek-ε、標(biāo)準(zhǔn)k-ω、SSTk-ω)和大渦模擬，其中標(biāo)準(zhǔn)k-ε模型在工程中應(yīng)用最多，計(jì)算量適中，有較多的數(shù)據(jù)積累且相對(duì)精確。自由界面的捕獲問(wèn)題一直是三維模型的一個(gè)難點(diǎn)，目前主要有無(wú)網(wǎng)格法、移動(dòng)網(wǎng)格法和固定網(wǎng)格法三類。如光滑粒子流體動(dòng)力學(xué)法(SPH)[4]在不使用計(jì)算網(wǎng)格的情況下，用一組流體粒子代替流體運(yùn)動(dòng)，模擬自由表面運(yùn)動(dòng)；移動(dòng)粒子半隱式法(MPS)[5, 6]，著眼于跟蹤粒子的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)性質(zhì)變化過(guò)程，各粒子之間通過(guò)核函數(shù)相互聯(lián)系；由Hirt和Nichols[7]提出的VOF法，通過(guò)跟蹤計(jì)算網(wǎng)格中每個(gè)計(jì)算單元相對(duì)于其中一個(gè)流體相的體積分?jǐn)?shù)來(lái)捕獲界面位置，是研究多相流最常見(jiàn)的方法之一，已被大多數(shù)潰壩水流模型所采用。Stansby P K等[8]通過(guò)試驗(yàn)研究了高度不同的兩立方水體的潰壩過(guò)程。Lauber和Hager[9]在平底矩形水槽中對(duì)潰壩洪水在下游無(wú)水情況下的演進(jìn)規(guī)律進(jìn)行研究，得出其水深和流速分布規(guī)律。王嘉松等[10]應(yīng)用組合型TVD格式, 模擬了瞬間全潰潰壩波的傳播、反射和繞射。肖瀟等[11]對(duì)比分析了潰壩水流的模擬方法，結(jié)果表明利用FLUENT得出的計(jì)算結(jié)果更為準(zhǔn)確，水位變化情況更加符合實(shí)際情況。王曉玲等[2]建立耦合VOF法與k-ε紊流模型的三維潰壩洪水演進(jìn)數(shù)學(xué)模型，模擬了三維潰壩洪水在復(fù)雜淹沒(méi)區(qū)域的演進(jìn)。袁岳等[1]結(jié)合有限元法和有限體積法處理二維淺水方程，模擬了潰壩水流流經(jīng)有障礙物的下游河道時(shí)的傳播過(guò)程，得到各時(shí)間段的水位變化和水流演進(jìn)情況?？偟膩?lái)說(shuō)，三維數(shù)值模擬已逐漸成為現(xiàn)階段研究潰壩水流的熱門方法，但結(jié)合紊流模型和自由面追蹤法系統(tǒng)地研究不同工況模型下水位變化規(guī)律和水流演進(jìn)的還較少，同時(shí)由于實(shí)際潰壩問(wèn)題具有很大的不確定性，利用數(shù)值模擬開(kāi)展?jié)嗡髁鲃?dòng)特性研究仍十分必要。

在前人研究的基礎(chǔ)上，本文利用FLUENT軟件建立潰壩三維數(shù)值計(jì)算模型，將下游有水的全潰壩、下游無(wú)水設(shè)障礙物的全潰壩、下游無(wú)水的局部潰壩3種工況整合到一起，采用VOF法對(duì)潰壩水流進(jìn)行自由表面追蹤，結(jié)合k-ε紊流模型模擬潰壩水流在演進(jìn)過(guò)程中的水位變化情況以及負(fù)波對(duì)水流傳播的影響。

1 基本方程和計(jì)算方法

1.1 基本方程

潰壩水流數(shù)值計(jì)算模型依據(jù)的基本方程有連續(xù)性方程和雷諾平均N-S方程。

連續(xù)性方程：

▽u=0

(1)

雷諾平均N-S方程：

(2)

式中：μeff=μ+μt；u為流速；t為時(shí)間；p為壓力；ρ為密度；f為外力；μ為動(dòng)力黏度；μt為湍流黏度，其中外力為重力，因此f=ρg，g為重力加速度。

上述方程采用保守形式，在含水區(qū)域內(nèi)求解控制方程時(shí)，由于密度是連續(xù)且恒定的，因此密度不包含在動(dòng)量方程的時(shí)間相關(guān)和對(duì)流項(xiàng)中。盡管與對(duì)流項(xiàng)相比，特別是在初始階段，潰壩流中的擴(kuò)散可以忽略不計(jì)，但為了保持雷諾平均N-S方程的一般形式仍然保留了擴(kuò)散相，如式(2)所示。

k方程：

(3)

ε方程：

(4)

(5)

(6)

式中：xi，xj為坐標(biāo)分量；Cu為經(jīng)驗(yàn)常數(shù)，建議取值為0.09，σk=1，C1ε=1.44，C2ε=1.92，σε=1.3。

1.2 VOF法追蹤自由表面

VOF法是根據(jù)每個(gè)時(shí)刻在網(wǎng)格單元中流體的體積變化來(lái)追蹤自由表面，其基本思想是，對(duì)于每一個(gè)計(jì)算單元，都有一個(gè)特定的標(biāo)量，表示給定單元的填充程度，例如水，如果在某個(gè)單元格中，該值為0，則為空；如果等于1，則充滿水；如果該值介于1和0之間，則可以說(shuō)該單元格在氣液交界面，也就是說(shuō)水的體積分?jǐn)?shù)α被定義為一個(gè)單元中的水體積與給定單元的總體積之比。相應(yīng)的，1-α表示給定單元中第二相的體積分?jǐn)?shù)。在VOF法中，氣液混合物的物理參數(shù)ρ、μ、μt可由體積分?jǐn)?shù)作平均。

ρ=αρ1+(1-α)ρ2

(7)

μ=αμ1+(1-α)μ2

(8)

μt=αμt1+(1-α)μt2

(9)

這里的1和2分別表示液體和氣體。流體體積分?jǐn)?shù)α滿足對(duì)流輸運(yùn)方程的解：

(10)

1.3 數(shù)值求解法

本模型數(shù)值求解采用PISO算法，該算法是SIMPLE的擴(kuò)展，在每個(gè)迭代步中增加了動(dòng)量修正和網(wǎng)格畸變修正過(guò)程，使得到的壓強(qiáng)場(chǎng)更準(zhǔn)確，計(jì)算收斂也越快。PISO算法包括一個(gè)預(yù)測(cè)和兩個(gè)校正，其目的是利用預(yù)測(cè)校正更好地滿足連續(xù)性方程和N-S方程。PISO算法步驟如下：

(1)預(yù)測(cè)：假設(shè)一個(gè)壓力場(chǎng)p*，利用p*求解動(dòng)量離散方程，得出相應(yīng)的速度場(chǎng)u*、v*。

(2)校正1。對(duì)假設(shè)的壓力場(chǎng)進(jìn)行修正，求解壓力修正方程得到p′，計(jì)算壓力修正量、速度修正量為u′、v′，得到壓力修正值和速度修正值p**、u**、v**。

p**=p*+p′

(11)

u**=u*+u′

(12)

v**=v*+v′

(13)

(3)校正2。

p***=p**+p″；p″=p*+p′

(14)

u***=u**+u″；u″=u*+u′

(15)

v***=v**+v″；v″=v*+v′

(16)

式中：p***、u***、v***分別為壓力修正值和速度修正值；p″、u″、v″分別為壓力修正量和速度修正量，設(shè)置p=p***，u=u***，v=v***，若修正后的壓力場(chǎng)相對(duì)應(yīng)的速度場(chǎng)能夠滿足連續(xù)性方程，則p、u、v為正確的壓力速度分量，否則將令p*=p，u*=u，v*=v，重復(fù)步驟(1)直至滿足條件。

2 潰壩洪水演進(jìn)分析

2.1 下游有水的全潰壩

(1)網(wǎng)格模型。在長(zhǎng)15.24 m、寬和高均為0.4 m的水庫(kù)中，閘門位于下游9.76 m處，水庫(kù)初始水深h0=0.1 m，下游水深為0.045 m，三維幾何模型如圖1所示，將空間尺度進(jìn)行歸一化處理，h/h0=0處為閘門位置，h為水位高度，重力加速度的方向與z軸相反。模型劃分為六面體結(jié)構(gòu)網(wǎng)格，數(shù)量約為224 700個(gè)。

圖1 三維模型圖(單位：m)Fig.1 Three-dimensional model

(2)邊界條件。由于模型整體偏細(xì)長(zhǎng)，單精度求解器可能不能足夠精確地表達(dá)各尺度方向的節(jié)點(diǎn)信息，因此本模型采用精確度更高的雙精度求解器計(jì)算，殘差值收斂標(biāo)準(zhǔn)為低于10-5。水庫(kù)頂部設(shè)置成壓力入口，出口邊界設(shè)置為壓力出口，出口位置壓力參考為一標(biāo)準(zhǔn)大氣壓，回流設(shè)為0。水庫(kù)底部和左右壁面設(shè)置為固壁邊界，默認(rèn)為是無(wú)滑移光滑壁面，采用壁面函數(shù)法。水庫(kù)上游處水體在重力作用下流向下游，設(shè)為內(nèi)部面邊界，流體體積分?jǐn)?shù)設(shè)為1。計(jì)入體積力，設(shè)置z方向重力加速度為-9.81 m/s2。下面兩種模型邊界條件設(shè)置相同。

(3)計(jì)算結(jié)果分析。對(duì)濕河床潰壩初期進(jìn)行數(shù)值模擬計(jì)算，不同時(shí)間(t=0.22 s、t=0.32 s、t=0.52 s、t=0.76 s)水位高度計(jì)算值與實(shí)驗(yàn)值[12]比較結(jié)果如圖2所示，圖3和圖4表示潰壩過(guò)程中不同時(shí)間(t=0.22 s、t=0.32 s、t=0.52 s、t=0.76 s)水面的變化。計(jì)算結(jié)果表明，閘門拉開(kāi)瞬間，上下游交界處水面有明顯的波動(dòng)，且以波浪形式向下游演進(jìn)，形成的第一個(gè)波浪波頂高度在演進(jìn)過(guò)程中變化不大[13](如圖3所示)，同時(shí)，t=32 s后由于負(fù)波的影響，上游水面出現(xiàn)波動(dòng)，波動(dòng)幅度隨水流演進(jìn)過(guò)程減小且并呈穩(wěn)定趨勢(shì)，閘門處水面也逐漸趨于平穩(wěn)。

圖2 不同時(shí)間水位高度變化模擬值與實(shí)驗(yàn)值對(duì)比Fig.2 Comparison of simulated and experimental values of water level height change at different time

圖3 不同時(shí)間水位變化圖Fig.3 Water level change chart at different time

圖4 不同時(shí)間水面的二維變化Fig.4 Two-dimensional variations of water surface at different times

2.2 下游設(shè)障礙物的全潰壩

(1)網(wǎng)格模型。取一長(zhǎng)9.03 m，寬0.3 m，高0.34 m的水庫(kù)，閘門位于水庫(kù)下游4.65 m處，水庫(kù)初始水深h0=0.25 m，下游為干河床，考慮與實(shí)際情況的吻合，障礙物設(shè)置成梯形，下游距閘門1.53 m 處，三維幾何模型如圖5所示，梯形障礙物截面尺寸如圖6所示。模型劃分為六面體結(jié)構(gòu)網(wǎng)格，數(shù)量約為20.71 萬(wàn)個(gè)。

(2)計(jì)算結(jié)果分析。不同時(shí)間下水位變化的數(shù)值模擬結(jié)果與實(shí)驗(yàn)值[14]的比較如圖7所示。從圖7可以看出，在t=1.9 s、t=2.8 s和t=6.68 s時(shí)，水面較為平滑，而在t=3.68 s和t=4.74 s時(shí)，水面出現(xiàn)小波浪，水面變得不穩(wěn)定。圖8和圖9表示潰壩過(guò)程中不同時(shí)間(t=1.9 s、t=2.8 s、t=3.3 s、t=3.68 s、t=4.74 s、t=6.68 s)水面的二維、三維變化。t=0 s時(shí)，水庫(kù)中的水體在重力作用下迅速流向下游，上游水位開(kāi)始下降，下游水位隨之上升，下游設(shè)置的梯形障礙物形成不規(guī)則地形，水流流經(jīng)下游障礙物時(shí)，水面發(fā)生劇烈變化，圖8可以看出，當(dāng)潰壩波到達(dá)梯形障礙物時(shí)，水流形成水躍；在t=1.9 s時(shí)，水面完全覆蓋了梯形障礙物，由于水流受障礙物的阻礙影響，梯形障礙物左側(cè)水位開(kāi)始上升；在t=2.8 s時(shí)，梯形障礙物左側(cè)的水位上升到一定高度，一些撞擊障礙物的水波向上游反向流動(dòng)；t=3.3 s時(shí)，障礙物左側(cè)上升的水位達(dá)到最大值，潰壩波破碎，產(chǎn)生的負(fù)波向上游傳播；t=3.68 s、t=4.74 s、t=6.68 s的水面線也顯示了負(fù)波向上游傳播的情況。綜上可知，在重力作用下，潰壩水流受下游障礙物的影響較大，產(chǎn)生的負(fù)波干擾了水流演進(jìn)。

圖5 三維模型圖(單位：m)Fig.5 Three-dimensional model

圖6 梯形障礙物截面尺寸(單位：m)Fig.6 Cross-sectional dimensions of trapezoidal obstacles

圖7 不同時(shí)間水位高度變化模擬值與實(shí)驗(yàn)值對(duì)比Fig.7 Comparison of simulated and experimental values of water level height change at different time

圖8 不同時(shí)間梯形障礙物水面的二維變化Fig.8 Two-dimensional variation of water surface of trapezoidal obstacle at different time

圖9 不同時(shí)間梯形障礙物水面的三維變化Fig.9 Three-dimensional variation of trapezoidal obstacle surface at different time

2.3 下游無(wú)水的局部潰壩

(1)網(wǎng)格模型。為分析局部潰壩水流沿程的運(yùn)動(dòng)，取長(zhǎng)為3 m,寬為2 m,高為0.7 m的水庫(kù)，兩個(gè)對(duì)稱的不透水墻組成局部有孔的大壩，大壩位于水庫(kù)下游1 m處，距左、右岸均為0.8 m,墻壁之間的距離為0.4 m,水庫(kù)的初始水位為0.6 m，下游河床干燥，三維幾何模型如圖10所示。模型劃分為六面體結(jié)構(gòu)網(wǎng)格，數(shù)量約為48.16萬(wàn)個(gè)。

圖10 模型三維圖(單位：m)Fig.10 Three-dimensional view of the model

(2)計(jì)算結(jié)果分析。在水庫(kù)上部、中部、潰口處與下游依次設(shè)置了G8A、GC、G5A、G0共4個(gè)監(jiān)測(cè)點(diǎn)，監(jiān)測(cè)點(diǎn)位置坐標(biāo)見(jiàn)表1。4個(gè)監(jiān)測(cè)點(diǎn)處水位隨時(shí)間變化的數(shù)值模擬計(jì)算結(jié)果與實(shí)驗(yàn)值[15]的比較如圖11所示。由計(jì)算結(jié)果可見(jiàn)，t=0 s時(shí)，潰壩開(kāi)始，水庫(kù)中的水體向下游運(yùn)動(dòng)，觀測(cè)點(diǎn)G5A與GC兩點(diǎn)處的水位隨時(shí)間遞減，隨后大約在t=0.28 s潰壩水流到達(dá)監(jiān)測(cè)點(diǎn)G8A，G8A處水位迅速升高后，隨著時(shí)間的推移逐漸趨于平緩。大壩中心處的G0點(diǎn)水位呈現(xiàn)出一定的波動(dòng)。由不同時(shí)間潰壩水流的二維和三維變化圖12和圖13可見(jiàn)，潰壩水流從潰口處流出后沿底部向兩側(cè)擴(kuò)散，并逐漸填滿下游區(qū)域。

表1 監(jiān)測(cè)點(diǎn)位置坐標(biāo)Tab.1 Position coordinates of monitoring points

圖11 各監(jiān)測(cè)點(diǎn)水位隨時(shí)間變化模擬值與實(shí)驗(yàn)值對(duì)比Fig.11 Comparison of simulated and experimental values of water level changes with time at each monitoring point

圖12 不同時(shí)間潰壩水流的二維變化Fig.12 Two-dimensional variation of dam-break flow at different times

圖13 不同時(shí)間潰壩水流的三維變化Fig.13 Three-dimensional variation of dam-break flow at different times

3 結(jié) 論

采用FLUENT軟件建立了三維潰壩水流模型，計(jì)算分析潰壩水流在下游的演進(jìn)過(guò)程?？紤]到實(shí)際潰壩問(wèn)題中存在許多的不確定性，例如地形起伏、河床侵蝕、水土流失、潰壩過(guò)程中的漸進(jìn)性等，分別考慮下游有水的全潰壩、下游無(wú)水設(shè)障礙物的全潰壩、下游無(wú)水的局部潰壩3種情況，用前人實(shí)驗(yàn)結(jié)果驗(yàn)證了數(shù)值模型的準(zhǔn)確性。通過(guò)對(duì)潰壩水流演進(jìn)過(guò)程和水位變化情況分析可得：

(1)下游有水時(shí)發(fā)生瞬間全潰，上下游交界處水面波動(dòng)明顯，而且出現(xiàn)的第一個(gè)波浪波頂高度最高，且在演進(jìn)過(guò)程中高度變化不大。雖然受負(fù)波影響上游水面出現(xiàn)波動(dòng)，但是隨著水流演進(jìn)波動(dòng)幅度減小，水面趨于穩(wěn)定。

(2)當(dāng)下游為有障礙物的不規(guī)則地形時(shí)，潰壩水流撞擊障礙物形成向上游的反射波，對(duì)水流演進(jìn)有一定影響，同時(shí)受障礙物阻擋作用影響，障礙物左側(cè)水位異常上升，而在實(shí)際情況中這種異常的水位上升可能會(huì)造成嚴(yán)重的環(huán)境破壞和經(jīng)濟(jì)財(cái)產(chǎn)損失。

(3)大壩發(fā)生局部破壞后，上游水位迅速下降，位于潰口處的水位逐漸降低且出現(xiàn)小幅波動(dòng)，潰壩水流整體從底部向兩側(cè)擴(kuò)散。

(4)潰壩水流受外界干擾條件的影響明顯，實(shí)際水流行進(jìn)過(guò)程中需盡量避免干擾，保證水流的平順。

□