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

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

0 引 言

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

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

1 數(shù)學模型及計算方法

1.1 數(shù)學模型

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

(1)

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

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

湍流動能k輸運方程

(2)

耗散率ε輸運方程

(3)

1.2 計算方法

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

連續(xù)方程

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

?u=0

(4)

N-S方程

(5)

能量方程

(6)

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

對于離散后的代數(shù)方程，選用FLUENT壓力速度耦合算法中收斂性較好的PISO算法進行求解。PISO算法適合瞬態(tài)不可壓縮流體，是典型的兩步校正算法[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***，若修正后壓力場對應的速度場能滿足連續(xù)性方程，則p、u、v為正確的壓力速度分量，否則令p*=p，u*=u，v*=v，繼續(xù)迭代。

2 數(shù)值模擬模型及計算條件

2.1 數(shù)值模擬模型建立

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

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

2.2 邊界條件

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

(2)出口邊界條件。在模型頂面和模型出口處設置為壓力出口邊界條件，參考壓力為一個標準大氣壓值，其中出口處壓力值與水深有對應關系，壓力值采用用戶自定義函數(shù)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種工況下的泄洪情況進行模擬，得出各監(jiān)測點壓力隨時間的變化，如圖3所示，各監(jiān)測點最大壓力對比見表1。由圖3、表1可以看出：3種工況下P1點最大壓力對比情況為工況3<工況2<工況1，到達時間為工況1<工況2<工況3；P3點壓力隨時間變化趨勢與P1點大致相同；P5點由于挑流作用出現(xiàn)負壓；P7點壓力值變化不明顯。此外，工況3下各監(jiān)測點處的壓力比工況1減少14%～42%，比工況2減少4%～28%。計算各個工況下洪水流過障礙物前后的壓力對比情況為：工況1壓力減小74%，工況2壓力減小73%，工況3壓力減小84%。

表1 監(jiān)測點壓力最值對比

圖3 不同監(jiān)測點壓力隨時間的變化

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

3.2 水面線分析

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

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

圖4 水面線隨時間的變化

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

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 影響分析

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

4 結 語

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