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

陳 松

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

1 概 述

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

2 溢洪道過(guò)流特性

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

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

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

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

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

3.1 控制方程

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

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

(1)

運(yùn)動(dòng)方程可以寫(xiě)成：

(2)

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

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

(3)

3.2 有限體積法(FVM)求解

針對(duì)溢洪道繞流問(wèn)題，建立基于顯式有限體積法(FVM)的數(shù)學(xué)模型。為了應(yīng)用FVM格式，首先可以方便地重新編寫(xiě)控制方程式(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ù)測(cè)-校正方案求解。由于這是一種顯式方法，計(jì)算時(shí)間步長(zhǎng)必須基于數(shù)值穩(wěn)定性考慮。

(10)

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

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

4 模型應(yīng)用

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

4.1 溢洪道物理模型

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

4.2 網(wǎng)格劃分

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

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

4.3 邊界和初始條件

4.3.1 溢洪道水流邊界條件

邊界條件指定了計(jì)算流域邊界上的流動(dòng)變量或其梯度。上游邊界可以設(shè)置在水庫(kù)斷面上，在該斷面上可知水庫(kù)水位和來(lái)水流量。該段應(yīng)遠(yuǎn)離溢洪道，以避免反射效應(yīng)?；谒俣确植嫉倪吔鐥l件為：

u=uo(x0,y)

(11)

(12)

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

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

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

(13)

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

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

(14)

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

4.3.2 溢洪道水流初始條件

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

4.4 結(jié)果和討論

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

圖2 溢洪道泄流壓力分布

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

圖4 溢洪道流速分布

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

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

5 結(jié) 論

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

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

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

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

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

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