下游障礙物位置對潰壩洪水流動特性的影響分析

林曉惠，付成華，王興華，白 帥，羅天璽

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

0 引 言

受水文、地質、施工等不確定因素以及暴雨、洪水、地震等自然災害的影響，潰壩事故時有發生，潰壩帶來的損失也是非常慘痛的。隨著計算機技術的發展，數值模擬越來越多地應用到潰壩問題的研究中。王曉玲等[1]針對潰壩洪水在復雜淹沒區域的演進，建立了耦合VOF法的三維k-ε湍流數學模型；Marrone等[2]采用delta-SPH模型分析潰壩水流沖擊不同形狀障礙物時的沖擊壓力，結果有較高的精度；Biscarini等[3]比較了淺水流動和三維模擬在不同潰壩條件下的計算結果。由國內外研究現狀可知，三維數學模型能更精細地模擬自由水面變化強烈的潰壩水流運動過程[4]，是現階段研究潰壩洪水泛濫區域、洪水深度、流速、洪水波傳播時間的熱門方法。

由于氣候變化，洪水泛濫區頻繁引起潛在傷亡和破壞，水流與障礙物間的相互作用導致災害等級增加，盡管存在這種相關性，但是文獻相對較少。本文利用FLUENT軟件建立潰壩洪水三維數值計算模型，耦合VOF法和RNGk-ε模型求解RANS方程，并采用PISO算法求解數值計算，模擬潰壩洪水在演進過程中的流動特性，進而研究障礙物位置變化對潰壩洪水流動特性的影響。

1 數學模型及計算方法

1.1 數學模型

1981年，Hirt和Nichols[5]提出了完整的VOF法理論體系和實現方法，廣泛應用于多種不混溶流體流動過程的截面捕獲。通過求解流體相體積分數方程，完成對各相之間界面的追蹤，從而實現對自由表面的追蹤[6]。方程如下

(1)

式中，相體積分數F=F(x，y，z，t)，定義為離散網格中各相流體的體積與網格體積的比值，且F∈[0，1]。

RNGk-ε湍流模型是由Yakhot等[7]使用重新歸一化組(RNG)方法開發的。

湍流動能k輸運方程

(2)

耗散率ε輸運方程

(3)

1.2 計算方法

基于有限體積法離散控制方程，包括連續方程、運動方程和能量方程[8]。對于黏性不可壓縮均質流體，ρ、μ、k均視為常數，于是控制方程組為

連續方程

水利現代化繪河清湖晏藍圖——訪江蘇省淮安市水利局局長、黨委書記黃克清……………………………… 韋鳳年，江 芳，郭 純等(16.59)

?u=0

(4)

N-S方程

(5)

能量方程

(6)

式中，ρ為密度；t為時間；u為流速；g為重力加速度；p為壓力；μ為流體黏度；cV為比定容熱容；T為熱量；Φ為耗散函數；k為熱傳導系數；q為輻射或其他原因單位時間傳給單位質量流體的熱量。

對于離散后的代數方程，選用FLUENT壓力速度耦合算法中收斂性較好的PISO算法進行求解。PISO算法適合瞬態不可壓縮流體，是典型的兩步校正算法[9]，主要實施步驟包括預估步、第一校正步、第二校正步。

(1)預估。假設壓力場p*，利用p*求解動量離散方程，得出速度場u*、v*。

(2)校正1。修正p*，求解壓力修正方程得到p′，計算壓力修正量、速度修正量為u′、v′，得到壓力修正值和速度修正值p**、u**、v**。

(3)校正2。對壓力修正方程進行修正得到p″，計算壓力修正量、速度修正量為u″、v″，得到壓力修正值和速度修正值p***、u***、v***。設p=p***，u=u***，v=v***，若修正后壓力場對應的速度場能滿足連續性方程，則p、u、v為正確的壓力速度分量，否則令p*=p，u*=u，v*=v，繼續迭代。

2 數值模擬模型及計算條件

2.1 數值模擬模型建立

基于文獻[10]的試驗模型，在水庫下游距閘門1.167 m處布設一45°弧形障礙物，距離左右岸均為0.295 m，初始水位H=0.55 m，簡化模型如圖1a所示。利用UG建立模型、ICEM劃分網格，導入FLUENT進行求解。采用結構化網格，最大網格0.005 m，最小網格0.001 m，網格總數為122 021個。

圖1 計算模型及測點布設(單位：m)

2.2 邊界條件

(1)進口邊界條件。重點關注水體自由跌落對障礙物的沖擊作用，故進口設置為墻(Wall)，以防止水流從進口處流出。

(2)出口邊界條件。在模型頂面和模型出口處設置為壓力出口邊界條件，參考壓力為一個標準大氣壓值，其中出口處壓力值與水深有對應關系，壓力值采用用戶自定義函數UDF給出。

(3)固壁邊界條件。障礙物及模型底部設置為固壁邊界(Wall)，采用無滑移邊界。糙率取值為0.012。

2.3 模型驗證

圖2為4個壓力測點處的壓力模擬值與試驗值的對比。從圖2可以看出，模擬的壓力值隨時間變化曲線與試驗值吻合程度較好，壓力分布能夠完全描述實驗并且充分模擬壓力峰值，可以使用該模型計算分析障礙物位置對潰壩洪水流動特性的影響。

圖2 各測點處的壓力模擬值與試驗值對比

2.4 計算工況

障礙物為45°弧形障礙物，初始水位為H=0.55 m，設置3種工況進行對比分析，即分別設置閘門與障礙物之間的距離L為0.583 5、1.167 m和1.750 5 m。

3 模擬計算結果及分析

3.1 壓力分析

對3種工況下的泄洪情況進行模擬，得出各監測點壓力隨時間的變化，如圖3所示，各監測點最大壓力對比見表1。由圖3、表1可以看出：3種工況下P1點最大壓力對比情況為工況3<工況2<工況1，到達時間為工況1<工況2<工況3；P3點壓力隨時間變化趨勢與P1點大致相同；P5點由于挑流作用出現負壓；P7點壓力值變化不明顯。此外，工況3下各監測點處的壓力比工況1減少14%～42%，比工況2減少4%～28%。計算各個工況下洪水流過障礙物前后的壓力對比情況為：工況1壓力減小74%，工況2壓力減小73%，工況3壓力減小84%。

表1 監測點壓力最值對比

圖3 不同監測點壓力隨時間的變化

下游障礙物位置距離潰壩閘門越遠，障礙物上各監測點的壓力峰值、均值越小，前后變化值相對越大。障礙物各監測點所受壓力負值的峰值與距離沒有明顯的對應關系；障礙物各監測點壓力最大值出現時間與距離成正比，即距離越遠，壓力峰值出現時間越晚。

3.2 水面線分析

3種工況下水面線隨時間的變化如圖4及表2所示。本次分析將潰壩洪水演進過程分為水流剛到達障礙物底端階段、水流到達障礙物頂面階段、水流完全越過障礙物階段、障礙物下游兩端包圍的無水區域被淹沒階段、水流到達出口階段5個階段。

表2 不同工況下水流到達各階段的時間對比

圖4 水面線隨時間的變化

從圖4及表2可以看出，工況3的挑流時長比工況1增加0.35 s，比工況2增加0.3 s。隨著間距的增加，水流到達障礙物各監測點的時間逐漸延后，當間距增大到一定程度時，第四階段與第五階段重合，并且間距越大，水流的挑流作用越明顯，通過障礙物時向前躍升高度越高。

3.3 流速分析

潰壩洪水向下游演進的過程中，遇到障礙物干擾，由于弧形障礙物對水流有躍升作用，導致水流越過障礙物后水流能量減少，速度明顯降低，之后再流向出口。計算3種工況下障礙物前、障礙物后以及出口平均流速變化的對比，見表3。

表3 平均流速變化對比

對比3種工況下障礙物前后平均速度的變化情況：工況1減小19.3%，工況2減小36%，工況3減小23.7%。考慮到L的不同，故對比3種工況下洪水出口處平均速度的關系，即工況2的平均速度比工況1增加12.2%，工況3的平均速度比工況1減小33.7%。因此，下游障礙物對水流流速有較大影響，且遠端障礙物更能有效地減緩水流流速，消能作用更強。

3.4 影響分析

下游障礙物對潰壩洪水的影響較大，潰壩洪水下泄越過障礙物時會發生挑流現象，且這個挑流的時間、距離隨著障礙物與閘門的間距發生變化。在洪水泛濫區，不同位置的障礙物，對潰壩洪水演進過程的影響不同，消能阻水的效果也不同。從模擬結果可見，離閘門遠端的障礙物對洪水沖擊力的減緩作用大于近端的障礙物。對于實際應用來說，下游障礙物削減的水能遠大于模型數據。因此，為達到良好的消峰消能目的，應當考慮下游障礙物的合適位置。

4 結 語

在已有試驗模型基礎上，本文利用潰壩洪水三維數學模型，模擬計算下游設置45°弧形障礙物的3種工況，得到了潰壩洪水在演進過程中的流動特性，以及障礙物位置變化對潰壩洪水流動特性的影響，對下游消能減災和防洪設計具有一定的參考意義。后續將結合潰壩方式、潰口變化、地形地貌影響、障礙物形狀大小變化等進一步研究潰壩洪水演進特性。