基于CFD模型的洪水條件下溢洪道水流特性模擬

陳興龍

(甘肅省水利水電勘測設計研究院有限責任公司，蘭州 730099)

0 引 言

自20世紀50年代初以來，國內一直在使用物理模型研究溢洪道上的水流行為[1]。工程師可以使用一系列水力設計圖來設計任何給定設計洪水位的溢洪道剖面[2]。CFD技術在溢洪道水流分析中的應用相對較新[3]。隨著計算機技術的進步，更高效的CFD程序可以在三維中求解Navier-Stokes方程，并且也有許多湍流模型可供選擇。

在溢洪道設計中，剖面設計應確保在最大洪水下水流過溢洪道結構時，不會造成不利影響，如頂部和下游的氣穴。理想情況下，溢洪道表面應剛好承受設計水頭下的大氣壓。當水庫水位低于該洪水位時，溢洪道上方的壓力將高于大氣壓力。當水庫水位遠高于設計水頭時，沿溢洪道頂部將產生負壓，該壓力可能由于氣穴而損壞溢洪道的混凝土面，并對包括閘門結構在內的其他部件產生不利影響。

國內的大多數大壩及其相關水工建筑物是在20世紀50年代和60年代早期設計的，以應對當時的洪水。此后，收集和處理了更可靠的長期水文數據。在許多情況下，修訂后的最大洪水等級有所增加。為了得出最佳補救設計，許多大壩需要考慮最具成本效益的方法，以調查最大洪水增加情況下溢洪道流量的行為。在過去的研究中，使用物理比例模型是唯一的調查方法。隨著計算機技術的發展，數值方法得以使用，并且在降低成本和減少分析時間方面更具有優勢。

在本次調查中，將分析各種洪水條件下的標準WES溢洪道剖面，研究二維和三維模型。為了驗證模型結果，將計算結果與公布的數據進行比較，結果顯示一致性良好，為今后研究中應用CFD提供了技術保障。

1 數值計算方法

本研究使用的CFD程序為FLOW-3D，通過有限差分法求解Navier-Stokes方程。該算法是一種基于SOLA方法的擴展方法，該方法由Hirt等在洛斯阿拉莫斯國家實驗室開發。流體體積(VOF)方法用于計算自由表面運動。

所有控制微分方程，如連續性和動量方程，均采用面積(2D)和體積(3D)孔隙度函數表示。該公式(分數面積/體積障礙物表示法)用于模擬復雜的幾何區域。

任何復雜的障礙物幾何體都可以使用偏好技術來表示。每個單元(柵格)中障礙物占據的體積部分(或2D中的面積)在分析開始時定義。還計算了每個單元中的流體分數，流體分數的連續性方程、動量方程或輸運方程是使用偏好函數來表示的，使用有限差分近似對每個方程進行離散。與某些有限元/體積或邊界擬合網格方法不同，這種網格技術不需要重新網格，并且在瞬態分析期間不會導致任何網格變形，因此可以很容易地應用精確的求解算法。

以一個時間增量推進解決方案的基本算法包括以下3個步驟：

步驟1：根據動量(Navier-Stokes)方程的顯式近似，使用所有平流、壓力和其他加速度的初始條件或以前的時間步長值計算速度。

步驟2：調整壓力以滿足連續性方程。

步驟3：更新流體自由表面或界面，根據流體體積給出新的流體配置。

2 二維溢洪道模型

2.1 模型網格及邊界條件

考慮無墩的WES溢洪道模型。之所以選擇這種方法進行分析，是因為測量結果不受任何三維效應的影響。因此，該模型接近于真實的二維流動問題。

溢洪道剖面的幾何形狀符合水力設計。它有一個垂直的上游面和一條曲線，由壩頂中心線前方的3個半徑(R=0.04Hd、R=0.20Hd和R=0.50Hd；Hd為設計水頭)定義。波峰中心線下游的剖面由以下方程式定義：

(1)

二維溢洪道頂部視圖見圖1。網格由X(水平)方向的95個單元和Z(垂直)方向的98個單元組成。縱橫比盡可能保持統一，特別是在敏感的區域，以提高求解精度和計算速度。對于該網格，X和Z方向的最大縱橫比分別為2.3和2.5。請注意，這里使用Z方向代替方程式(1)中定義的Y方向。

圖1 二維溢洪道頂部的網格特寫

本次調查的設計水頭取10 m。左邊界(上游)距離壩頂25 m，右邊界(下游)距離壩頂22 m。底部邊界位于壩頂下方18 m，頂部邊界位于壩頂上方14 m。假設以下邊界條件：

上游邊界：流速為零的靜水壓力；流體高度等于H。

下游邊界：流出邊界。

上游底部：無流動—被下方障礙物堵塞。

下游底部：流出邊界。

頂部邊界：對稱性—在這種情況下，由于重力，無影響。

設置初始條件，使水頭為H的流體體積位于溢洪道頂部。當達到穩定狀態時，在15 s的總時間內進行瞬態流分析，這是通過檢查系統的流速和動能等結果確定的。使用1 000 kg/m3的恒定水密度，這里假設水是不可壓縮的。在負Z方向施加9.81 m/s2的重力加速度。檢查3種不同的水頭工況(H/Hd=1.33、1.00和0.50)。

2.2 壓力分布

穩態時沿波峰的表面壓力分布見圖2，這些壓力取自溢洪道附近的單元。圖2中還繪制了測量數據。可以觀察到，計算結果給出了略高的負壓，但總體趨勢和量級與測量數據非常一致；還可以看到一些壓力振蕩，它們可能歸因于局部網格效應。

對于設計水頭工況(H/Hd=1.00)，水流沿溢洪道產生的壓力接近于零，即使在數值模擬中未引入曝氣。當壓頭高于設計壓頭時，出現負壓；當水頭低于設計水頭時，產生正壓。圖3顯示了穩定狀態下溢洪道頂部H/Hd=1.33的壓力等值線，可以觀察到頂部上方的負壓區域。

圖2 不同水頭頂部壓力分布計算值與測量結果的比較

圖3 溢洪道頂部上方的負壓分布(單位：Pa)(H/Hd=1.33)

2.3 流速與流量分布

計算假設壁面完全光滑，且流動無黏性，湍流的影響將是未來研究的主題。根據試驗結果，尖頂堰/溢洪道上的流量可表示為：

Q=CLH1.5

(2)

式中：Q為流量，m3/s；C為流量系數；L為堰頂有效長度，m；H為壩頂上方的實測水頭，不包括流速水頭。

其中：流量系數近似為：

C=3.27+0.40(H/h)

(3)

式中：h為堰高。

穩態下波峰上方的速度矢量見圖4，每種情況下的計算流量和平均速度(波峰中心線X方向的速度)見表1。為了便于比較，表1中還顯示了式(2)中的流量和相應的平均流速。可以看出，計算值比經驗結果高估了約10%～20%，這可能與分析中使用的無黏流條件有關。

圖4 溢洪道壩頂上方的速度矢量(單位：m/s)

表1 流量和平均流速的比較

3 三維溢洪道模型

現在考慮帶橋墩的WES溢洪道模型。橋墩的存在將產生三維效應，此分析使用了相當粗糙的網格。除了模型中包含橋墩和橋臺外，幾何結構與二維情況類似。對稱條件應用于兩個z-x邊界平面。上游z-y邊界面上保持恒定的靜水壓頭，允許水通過下游的z-y和x-y邊界平面流出。穩態下的模型視圖見圖5。

對結果的初步檢驗表明，分析結果與公布的數據有較好的一致性，橋墩附近的負壓高于其余部分的負壓。

4 溢洪道結構的應用

鑒于計算結果與二維分析中公布的數據之間的良好一致性，使用相同的技術對混凝土重力壩中央溢洪道在洪水位增加的情況下進行分析。首先建立一個物理模型，對溢洪道在不同洪水位下的特性進行研究；繪制物理模型試驗的測量值，以進行比較，沿溢洪道中心線計算的壩頂壓力分布見圖6。

圖5 顯示穩定狀態下流過溢洪道的含水量和流速分布

圖6 不同洪水位的峰值壓力分布CFD分析和物理模型試驗結果的比較

由圖6可以觀察到，獲得了相對較好的一致性。

5 結 論

本文基于CFD模型，從二維和三維兩個方面對溢洪道洪水特性進行了模擬分析，并將計算結果與模型試驗結果進行了驗證，取得了良好的一致性。但計算結果高估了流速，低估了溢洪道沿線的壓力分布。