溢洪道水流的二維有限體積法數值模擬

陳 松

(水利部新疆維吾爾自治區水利水電勘測設計研究院， 烏魯木齊 830000)

1 概 述

溢洪道和其他泄洪建筑物設計的主要目的是安全地將洪水從大壩輸送到下游河道，并防止大壩漫頂，特定類型溢洪道的選擇和設計基于項目的特定目的、水文、泄洪要求、地形、地質、大壩安全和項目經濟。由于提供高效輸送水力和結構堅固的溢洪道對大壩的安全以及下游河流的生命和財產非常重要，因此水力模型被廣泛用于解釋水力現象的特性[1]。在綜合考慮泄流能力、壓力、水面線、消能布置等水力設計因素的情況下，對溢洪道進行水力高效設計，而且溢洪道寬度對大壩的經濟分析和安全運行起著重要作用。目前，世界范圍內大壩安全的主要問題是提供足夠的溢洪道泄洪容量，但許多水壩設計和建造于20世紀60和70年代，水文信息有限。當時設計的許多水壩的溢洪道容量不足，可能會遇到潛在的問題，如在洪水條件增加的情況下，溢洪道頂部產生過大的負壓，從而威脅大壩安全，因此溢洪道水動力學知識對溢洪道的安全設計具有重要意義[2]。本文針對溢洪道水流的流體力學問題，借助于數值模型，對溢洪道的流速分布、壓力分布和流量特性進行模擬研究。

2 溢洪道過流特性

傳統的溢洪道又稱為反弧型溢洪道，主要由4部分組成：進水段、控制段、泄槽和消能防沖段。在壩頂上游，水流處于亞臨界流狀態(或逐漸變化的)；在壩頂之后，由于溢洪道為陡峭斜坡面，使得水流由亞臨界狀態變為超臨界狀態。溢洪道水流本質上是具有明顯曲率且快速變化的水流，在頂部的水流中同時發生兩個現象：沿剖面形成湍流邊界層并且逐漸增厚，以及主流速度的逐漸增加和深度的逐漸減小。

溢洪道泄洪時由于斜坡面陡峭，水流流速過高，從而產生水流摻氣。由于夾帶空氣而導致的流動膨脹增加了自由表面，從而間接影響設置側壁高度和選擇徑向門的耳軸高度。對于溢洪道水流，垂直加速度非常重要，流速和壓力會沿水流方向以及垂直方向變化。溢洪道水流基本上是在頂部附近快速變化的水流，在垂直方向上流線明顯彎曲。由于溢洪道實體表面的曲線性質以及壩頂下游的陡坡，與實體邊界的抗剪強度相比，垂直加速度在流動中起主導作用。

3 溢洪道水流數值模擬

利用高性能計算機和更高效的CFD(computational fluid dynamics)程序，可以在合理的時間和成本內對水工結構的性能進行數值研究。以往研究人員大多采用有限差分法模擬溢洪道水流。由于有限差分法對曲線體幾何的適用性較弱，它在任意曲線體邊界的重力驅動自由表面流動中的應用受到很大限制。而有限體積法、有限元法、邊界元法等對曲面實體邊界具有良好的適應性，得到了廣泛的應用。在高速可壓縮流動模擬中，守恒對于準確捕捉沖擊波和其他流動不連續性是非常重要的，也是有限體積法的優點[3-5]。

本文基于空間弱可壓縮流動方程，建立溢洪道繞流的二維有限體積法數值模擬模型。對于在垂直方向上變化迅速且在橫向(如溢洪道上的水流)上顯示出可忽略的水流變化進行建模，可以使用二維垂直模型。該模型的簡化形式可作為一個寬度平均模型。

3.1 控制方程

自然界中典型的自由表面流具有較大的雷諾茲數，其特征是大尺度湍流混合。這些流動的馬赫數很小，通常被視為不可壓縮流動。在快速變化的條件下，一階壓縮性是不可忽略的。實際上，流體是可壓縮的，壓力總是與密度有關。因此，即使對于所謂的不可壓縮流型，也可以使用可壓縮形式的方程。為了使問題更易于數值求解，引入可壓縮性。Song和Yuan發展了適用于小馬赫數和大雷諾數實際流動的弱可壓縮流體力學方程[6]。

對于包括瞬變流在內的一般馬赫數流，其連續性方程為：

(1)

運動方程可以寫成：

(2)

式中：p為壓力，kPa；ρo為流體密度，kg/m3，ao為聲速，m/s；V為速度矢量，m/s。

在上述方程中，τij是剪應力張量，假設重力作用于Y方向。滿足式(1)和式(2)的流動稱為弱可壓縮流動。忽略運動方程式(2)中的黏性項，得到的運動方程稱為歐拉方程：

(3)

3.2 有限體積法(FVM)求解

針對溢洪道繞流問題，建立基于顯式有限體積法(FVM)的數學模型。為了應用FVM格式，首先可以方便地重新編寫控制方程式(1)和式(3)，其保守形式如下：

(4)

其中：G為流量變量(p，u，v)；F為流量矢量。

F=iE1+jE2

(5)

G=[puv]T

(6)

(7)

(8)

方程式(4)可以在任意有限體積上積分，應用散度定理，將體積積分轉化為曲面積分。把它平均到現有的體積：

(9)

其中：G為G的體積平均值，m3；V為體積，m3；n為單位法向量；s為控制體積的表面積，m2。

方程式(9)由兩步顯式預測-校正方案求解。由于這是一種顯式方法，計算時間步長必須基于數值穩定性考慮。

(10)

其中：s為控制體積的表面積，m2。在滿足方程式(10)條件的整個區域上，計算的Δt的最小值作為計算時間步長。

在此基礎上，建立溢洪道水流動力學數值模型。在所建立的模型中，從給定的初始條件出發，求解下一步所有網格點的速度和壓力控制方程，計算新的自由曲面，生成新的網格系統，根據顯式時間推進法計算每一時間步的流量，直至得到穩態解。

4 模型應用

本文以某攔洪壩溢洪道為例，利用所建立的數值模型對其進行水動力模擬。溢洪道設計洪水流量為55 m3/s。數值模型中模擬的部分原型包括：部分水庫(長58.61 m，深39.144 m)；上游溢洪道壩頂剖面長8.61 m，下游反弧壩頂剖面長27.862 m。溢洪道部分的反弧壩頂剖面考慮到下游切點。

4.1 溢洪道物理模型

大壩溢洪道的物理模型研究是以1∶50的幾何相似比例建造的二維截面模型上進行的。為了觀察溢洪道上游和下游的水流狀況，包括消能工的性能，在玻璃水槽中復制了模型。一個完整的跨度和兩邊的兩個半跨度被復制。沿跨距中心在溢洪道表面安裝直徑為5 mm的壓力計，用于測量壓力；流量測量采用1.5 m寬尖頂堰；水位測量采用最小計數為0.1 mm的指針式儀表，從而對溢洪道模型的泄流能力、溢洪道表面的壓力和水面線進行水力模型研究，最后將數值模擬結果與大壩溢洪道的物理模型研究結果進行比較。通過物理模型研究，主要從使用功能方面對溢洪道結構進行優化設計，以保證工程的技術的可行性。

4.2 網格劃分

建立貼體網格系統進行有限體積法建模。圖1為網格系統，將流場離散為18 400(230×80=18 400)個四邊形單元網格。

圖1 溢洪道網格劃分

4.3 邊界和初始條件

4.3.1 溢洪道水流邊界條件

邊界條件指定了計算流域邊界上的流動變量或其梯度。上游邊界可以設置在水庫斷面上，在該斷面上可知水庫水位和來水流量。該段應遠離溢洪道，以避免反射效應。基于速度分布的邊界條件為：

u=uo(x0,y)

(11)

(12)

其中：x為流動方向；y為垂直方向；u、v為X和Y方向的速度分量；uo為上游邊界X方向的速度分量，uo值由給定流量和水深確定。對于壩頂形狀效應的研究，由于溢洪道下游斜坡的水流是超臨界的，下游條件對上游水流沒有影響。下游段可選擇在流量充分發展的斜坡段上，以便假定速度和壓力梯度為零。

固體邊界條件有3種：完全滑移邊界條件、部分滑移邊界條件和無滑移邊界條件。全滑移是指內網格處的切向速度等于固體表面上的切向速度；無滑移是指固體表面上的切向速度為零；對于部分滑移情況，應使用壁面函數。3種邊界條件的選擇取決于網格大小和邊界層厚度的相對大小。在溢洪道和水庫的固體邊界處，采用全滑移條件。固體邊界沒有流動。

最難模擬的邊界是時變自由面位置。影響自由表面的因素很多，如風應力、水與空氣的熱交換、表面張力應力等，通常采用運動學和動力學條件對自由表面進行模擬。運動條件是基于自由面是一個物質表面的觀點，它的形式是：

(13)

式中：Zf為沿法向的自由表面位移；u為沿流動方向的速度矢量；vf為垂直方向的自由表面速度，m/s。

該方程在迭代過程中，采用Mac-Cormack格式來修正新的自由面位置。忽略表面張力效應的動態邊界條件是零應力條件。對于該模型，自由表面零應力的動態邊界條件簡化如下：

(14)

其中：n為垂直于自由曲面的單位向量。

4.3.2 溢洪道水流初始條件

溢洪道的初始流量為66 236 m3/s，其設計最大洪水流量為88 315 m3/s，相應的排放強度為143.99 m3/(s·m)。

4.4 結果和討論

經過60 000個時間步，期間模擬結果不時地被檢驗，直到得到一個穩定的狀態解。整個計算過程在1.0 GHz的處理器上進行了約12 h。并將物理模型上觀測到的測壓壓力、水面線與模擬壓力、水面線進行了比較，計算值與實驗值吻合較好。圖2為收斂后的壓力分布(以等高線的形式)，同時也給出不同水位對應的壓力水頭。圖3為解決方案收斂后固體表面上的壓力分布，并與物理模型結果進行比較，結果顯示模型觀測壓力與模擬計算壓力結果相吻合。從圖2和圖3可以很明顯地看出，在流動緩慢且穩定的水庫部分，靜水壓力占主導地位。從圖4和圖5可以看出，在數值模型中上游壩頂剖面上形成一個溫和的分離帶，由于重力作用水流勢能不斷轉化為動能，速度整體變大。由速度矢量圖可以看出，溢洪道表面流速前半段增速較快，后半段增速較慢；而溢洪道底部水流前半段增速較慢，后半段增速較快，在溢洪道的后半段同一切面的流速相同。

圖2 溢洪道泄流壓力分布

圖3溢洪道面壓力分布與物理模型計算結果的比較

圖4 溢洪道流速分布

圖5 溢洪道水流的流線形態

流量系數的計算值和實驗值分別為0.72和0.69，流量系數計算值比實測值高4%。這可能是由于在寬均二維數值模型中，忽略了黏性項和紊流，假定了良好的直線流動，沒有由于橋墩和橋臺的端部收縮而造成的損失。數值模型顯示了溢洪道壩頂附近的臨界流段位置，其中無量綱水力參數弗勞德數是統一的，表明流態由次臨界變為超臨界流。

5 結 論

本文利用弱可壓縮流動方程，建立了基于有限體積法的溢洪道壩頂斷面水力特性數值模型。對選定溢洪道的流速分布、壓力分布和泄流特性進行了估算，并與已有的物理模型資料進行了比較。主要結論如下：

1) 流量系數的計算值和實驗值分別為0.72和0.69，比物理模型的實驗值高4%。計算值較高的原因可能是忽略了黏性項和湍流，這通常會在水流中產生阻尼效應，并且在二維數值模型中沒有由于橋墩和橋臺引起的端部收縮造成的損失。

2) 從描繪壓力等值線和流線的圖中可以很明顯看出，水庫部分存在靜水壓力，溢洪道部分由于快速變化的加速流量而受到非靜水壓力的影響。

3) 數值模型顯示了溢洪道頂部附近的臨界流量段(弗勞德數=1)的位置，這與模型研究的既定事實非常吻合。

4) 上游壩頂剖面未能正確引導水流通過壩頂，因此在數值模型中，上游壩頂剖面上形成一個溫和的分離區。

數值模擬結果與物理模型計算結果吻合較好，說明該模型在溢洪道實例水力模擬中的適用性，可將該模型應用于溢洪道泄洪水流特性的模擬，從而進一步加強大壩的安全運行。