基于FLUENT三維數(shù)值模擬的潰壩洪水流動特性分析

王興華，付成華，穆宵梟

(西華大學能源與動力工程學院，成都 610039)

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

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

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

1 基本方程和計算方法

1.1 基本方程

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

連續(xù)性方程：

▽u=0

(1)

雷諾平均N-S方程：

(2)

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

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

k方程：

(3)

ε方程：

(4)

(5)

(6)

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

1.2 VOF法追蹤自由表面

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

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

(7)

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

(8)

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

(9)

這里的1和2分別表示液體和氣體。流體體積分數(shù)α滿足對流輸運方程的解：

(10)

1.3 數(shù)值求解法

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

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

(2)校正1。對假設(shè)的壓力場進行修正，求解壓力修正方程得到p′，計算壓力修正量、速度修正量為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***，若修正后的壓力場相對應(yīng)的速度場能夠滿足連續(xù)性方程，則p、u、v為正確的壓力速度分量，否則將令p*=p，u*=u，v*=v，重復步驟(1)直至滿足條件。

2 潰壩洪水演進分析

2.1 下游有水的全潰壩

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

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

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

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

圖2 不同時間水位高度變化模擬值與實驗值對比Fig.2 Comparison of simulated and experimental values of water level height change at different time

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

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

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

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

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

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

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

圖7 不同時間水位高度變化模擬值與實驗值對比Fig.7 Comparison of simulated and experimental values of water level height change at different time

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

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

2.3 下游無水的局部潰壩

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

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

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

表1 監(jiān)測點位置坐標Tab.1 Position coordinates of monitoring points

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

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

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

3 結(jié) 論

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

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

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

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

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

□