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

陳 松

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

1 概 述

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

2 溢洪道過流特性

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

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

3 溢洪道水流數(shù)值模擬

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

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

3.1 控制方程

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

對于包括瞬變流在內(nèi)的一般馬赫數(shù)流，其連續(xù)性方程為：

(1)

運(yùn)動方程可以寫成：

(2)

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

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

(3)

3.2 有限體積法(FVM)求解

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

(4)

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

F=iE1+jE2

(5)

G=[puv]T

(6)

(7)

(8)

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

(9)

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

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

(10)

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

在此基礎(chǔ)上，建立溢洪道水流動力學(xué)數(shù)值模型。在所建立的模型中，從給定的初始條件出發(fā)，求解下一步所有網(wǎng)格點(diǎn)的速度和壓力控制方程，計算新的自由曲面，生成新的網(wǎng)格系統(tǒng)，根據(jù)顯式時間推進(jìn)法計算每一時間步的流量，直至得到穩(wěn)態(tài)解。

4 模型應(yīng)用

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

4.1 溢洪道物理模型

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

4.2 網(wǎng)格劃分

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

圖1 溢洪道網(wǎng)格劃分

4.3 邊界和初始條件

4.3.1 溢洪道水流邊界條件

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

u=uo(x0,y)

(11)

(12)

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

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

最難模擬的邊界是時變自由面位置。影響自由表面的因素很多，如風(fēng)應(yīng)力、水與空氣的熱交換、表面張力應(yīng)力等，通常采用運(yùn)動學(xué)和動力學(xué)條件對自由表面進(jìn)行模擬。運(yùn)動條件是基于自由面是一個物質(zhì)表面的觀點(diǎn)，它的形式是：

(13)

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

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

(14)

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

4.3.2 溢洪道水流初始條件

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

4.4 結(jié)果和討論

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

圖2 溢洪道泄流壓力分布

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

圖4 溢洪道流速分布

圖5 溢洪道水流的流線形態(tài)

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

5 結(jié) 論

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

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

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

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

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

數(shù)值模擬結(jié)果與物理模型計算結(jié)果吻合較好，說明該模型在溢洪道實(shí)例水力模擬中的適用性，可將該模型應(yīng)用于溢洪道泄洪水流特性的模擬，從而進(jìn)一步加強(qiáng)大壩的安全運(yùn)行。