耦合VOF法的河道三維水流特性數值模擬研究

周正印，楊 楠

（天津市市政工程設計研究院，天津市 300392）

1 研究背景

隨著我國經濟的深入發展，國家對生態環境的保護越來越重視，污染防治是2018年我國提出的三大攻堅戰之一，被放在了前所未有的重要位置上。但是長期以來，我國污染治理進程相對滯后，河湖空間管理和控制能力相對薄弱，水流不暢、水岸雜亂、水景不佳、水體黑臭等水環境問題日趨惡化，內澇等水安全問題日益凸顯，成為經濟社會和生態環境協調發展的制約因素。城市中小河道水體黑臭是點源污染、面源污染、內源污染等綜合因素導致的結果[1]，其中喪失生態功能的水體，往往流動性降低或完全消失，導致水體復氧能力衰退，形成缺氧、厭氧的水環境，引發水質惡化，加劇水體黑臭。目前針對黑臭水體的治理方法主要是包括截污、清淤、生態、曝氣、活水等一系列方法在內的綜合處理[2]，而其中河道水動力條件是實施上述治理措施的依據，因此對河道水流進行數值模擬，改善河道水動力條件，可為黑臭水體治理提供前期支持。

目前，對河道水流特性進行數值模擬的研究起步較早，研究較多。在國外，Renato[3]采用二維高精度淺水方程對意大利2014年Secchia河的洪水事件進行了數值模擬，為達到較快的計算速度，采用了并行計算方法；Siavash[4]以伊朗Karun河的實測數據為基礎，采用二維數值模型，將固定的曼寧糙率改為隨河道條件相應動態變化的糙率，提高了河道水流數值模擬的精度；Ao[5]考慮水質，采用三維數值模型對河道水流進行了數值模擬，為河道水質管理提供了技術支持。在國內，張瑋[6]針對目前大多河流存在著平面順直、斷面單一的問題，以奧地利某順直河道段低水生態修復工程為例，利用二維水流數學模型，研究修復工程的平面布局；彭畢帥[7]采用二維淺水方程對長江下游黑沙洲水道在不同流量下灘脊線加高、岸線崩退、左槽封堵3中工況下的水流變化進行了數值模擬；王甲榮[8]建立基于有限體積法的平面二維水動力模型對海河流域魯北平原典型河流彎道進行數值模擬，分析不同工況下河床、河道邊坡的沖淤變化，選出適宜河流彎道生態治理的河床形態。

綜上所述，目前國內外對河流水流特性的數值模型多集中于一、二維，忽略了對河流自由水面問題的研究，并對水流豎直方向上的變化進行了平均化處理，缺乏對河流自由水面波動的模擬，簡化了實際河流的復雜彎曲流態與水流細節，降低了模擬的精度。針對上述問題，本文建立了三維高精度河道地形模型，考慮河道水流的三維非穩態運動，建立耦合流體體積函數（volume of fluid，VOF）法與紊流模型的三維數值模型，研究了河道干流段水流三維運動規律，得到各斷面水流的水深、流速等特征值，為進一步改善河流水動力條件及水生態措施的實施提供理論依據和決策支持。

2 研究方法

2.1 耦合VOF法的三維k-ε模型

流體力學中淺水流動模型（即圣維南方程組）是對實際流動的一種簡化和概化的數學模型，無法精確、穩定地模擬河道水流問題[9]。因此，本文提出一種耦合VOF法和k-ε雙方程紊流模型的數值模擬方法，使之能夠適用于求解更廣泛的具有復雜邊界形態和流動特征的自由表面流動。

VOF法允許較陡的自由表面和非單一表面，對河流水質遷移轉化自由水面數值模擬采用VOF法可精確計算各計算斷面水位、自由水面形狀、污染物分布及流速等變量，并以可視化形式表達。VOF法通過定義一個流體體積函數F，F是空間和時間的函數，即 F（x，y，z，t）。在離散網格內，F 取值是網格內各相流體的體積與能夠被流體通過的空間體積的比值。在每一個單元內，變量或代表一相，或代表多相，用aq表示單元內第q相流體的體積分數。在每一個控制體中，所有的相體積分數之和為1，即：

對于第q相流體，體積分數連續方程為：

式中：u、v、w 分別為主場沿 x、y、z方向的分速度；Saq為該相的質量源，無源情況下為零。

k-ε紊流模型的連續性方程、動量方程、紊動能k方程和ε方程如下：

連續性方程：

式中：t為時間，s；ui為速度分量，m/s；xi為坐標分量；ρ為密度，kg/m3；μ 為分子粘性系數，N·m/s；P為修正壓力，Pa；ui為紊流粘性系數；CD為源項；k為紊動動能，m2/s2；ε 為紊動耗散率，m2/s2；?k、?ε 分別為 k、ε 的紊流普朗特數，無因次；C1ε、C2ε為經驗常數，無因次[10]。

2.2 模型的邊界條件

（1）進口邊界條件：如果進口條件已知，給定速度如下式，湍流參數和根據經驗公式計算得到；如果不知道準確的來流條件，只知道來流的流量條件（Q為進口流量），A為進口面積，可以速度方向給定在結構化網格方向上（如河道），并設m→為網格的方向矢量，則進口速度給定如下：

（2）出口邊界條件：通常在計算域的出口，各速度分量(u，v，w)以及k和ε均取為第二類齊次邊界條件，即（n→代表出口斷面的法向）：

（3）固壁邊界條件：在壁面上采用無滑移條件，即：u=v=w=0。為計算近壁區的紊流，采用壁面函數法。

（4）壓力邊界條件：在計算域的邊界上，壓力應滿足Neumann條件，即第二類邊界條件[11]。為保證計算的穩定性，在規定的某一內點上，壓力為給定值，且定義為參考壓力Pref。

2.3 模型求解

采用有限體積法離散控制方程，差分格式采用QUICK差分格式[12]，并采用基于SIMPLE算法改進的壓力的隱式算子分割算法（Pressure Implicit Split Operator，PISO 算法）對方程進行求解[13]。

3 模型驗證

針對三維模型的特點，采用S.Soare Frazao與Y.Zech在急轉彎河道模型中進行的試驗結果作為驗證的依據[14]。

試驗在一個有90°急彎的模型河道中進行。上游設置一個水箱用于儲水，平面大小為2.44 m×2.44 m×2.39 m，河道為矩形斷面，寬0.495 m，玻璃壁面；上游河道長3.92 m，急彎下游河道長2.92 m，河道床面高出水庫底面0.33 m。末端自由出流，水箱中的初始水位高出河道床面0.25 m，實驗布置見圖1。模型采用貼體網格技術對模型進行劃分，網格單元的大小為x×y×z=0.01 m×0.01 m×0.01 m，共劃分出88 130個網格，計算網格見圖2。

圖1 潰壩實驗模型示意圖(單位：m)

圖2 網格模型

實驗前河道為干河床，試驗時將水箱閘門快速提起，模擬河水在河道中的流動，通過玻璃邊壁處的攝影捕捉某些斷面的水位變化過程，圖3為沿程水深模擬值與實驗值對比。將模擬得到的t=3 s與t=7 s時下有河道中水深與實驗值對比，模擬結果與實驗結果吻合，最大誤差為9.7%，最小誤差為3.1%，平均誤差為6.5%。

4 工程應用

4.1 工程簡介

阜陽市位于安徽省西北部，淮河以北，是淮海經濟區重要組成部分。近年來，阜陽市的經濟建設取得了長足的發展，但是污染治理進程相對滯后，由于合流制管道排污、初期雨水污染、水體流動性差、污染物淤積等因素，導致部分河道水體黑臭。為改善阜陽市城區人居環境，提升城市品位，阜陽市決定啟動阜陽市城區水系綜合整治（含黑臭水體治理）項目。本文以其中某典型黑臭水體河道為研究對象，應用建立的耦合VOF法和雙方程的紊流數值模型對河道水流進行數值模擬，為進一步的黑臭水體綜合治理措施提供科學依據和決策支持。

圖3 沿程水深模擬值與實驗值對比(單位：m)

4.2 計算模型及網格劃分

該河道長5.2 km，河水自西向東流動，上游為盲端，下游接入另一條河道，中間與3條城區河道交叉，加上沿河雨水及合流制排口，為該河道的主要來水。依照河道平面CAD布置圖及斷面圖，將包含河底形狀及高程等實際地形數據的格式文件導入到CFD軟件中，實現了真實地形資料在CFD軟件計算網格模型中的精確表達，彌補了以往CFD建模中通過坐標繪制網格模型而使網格精確性不足的局限。采用三角形網格、無結構貼體網格劃分技術和局部加密網格生成技術建立干流河道的計算網格模型，河道網格單元總數為712718個。河道流域及河道干流的三維網格模型的全局圖見圖4。由于河道長寬高比例大，因此在全局顯示的情況下將使得河道在Y、Z方向網格幾乎無法清楚顯示，河道幾乎成為線狀，但河道各彎道及轉彎仍可見。

根據該河道設計要求，確定該河道設計流量為10.40 m3，河道一設計流量為11.20 m3，河道二設計流量為22.10 m3，河道三設計流量為30.10 m3，水動力計算采用進口流速邊界，出口采用流速邊界，并保證總進口流量等于總出口流量，河道基本達到穩定狀態。

圖4 干流河道網格模型

為了能夠更好的分析河道水力條件，在河道干流設置了8個典型斷面監測水流變量，分別在河道起始點、河道交叉處、河道中間、河道轉彎、河道末端等處，斷面詳細布置情況見圖5。

圖5 河道干流8個典型斷面分布

4.3 VOF動態及水流流速分析

圖6為干流河道已基本達到穩定狀態時，6個典型斷面的VOF分布情況。根據圖3模擬結果計算水深可得，E1、E2、E3、E4、E5、E6、E7 和 E8 斷面的平均水深值分別為2.79 m、2.98 m、3.00 m、3.02 m、3.30 m、3.36 m、3.42 m和3.43 m，斷面處的水深分布見圖7。

由圖7可以看出，隨著斷面從上游向下游推進，河道水深逐漸變深，并且斷面水平面的波動比上游劇烈。斷面1由于處于河道起點，且為盲端，沒有其他補充水源，水平面較為平穩，EI斷面最小水深值為2.790 m，最大水深值為2.791 m；E2斷面位于本研究河道與河道1交叉處，河道斷面具有一定的波動性，水深比E1斷面大，最小水深值為2.978 m，最大水深值為2.983 m；E3斷面位于本研究河道與河道2交叉處，水深比E2斷面大，最小水深值為2.997 m，最大水深值為3.003 m；E4斷面位于河道順直段，水流較為平緩，水位變化不大，最小水深值為3.019 m，最大水深值為3.020 m；E5斷面位于本研究河道與和河道3交叉口處，且該處河道轉彎，水流變化最為劇烈，水深增加較大，最小水深值為3.295 m，最大水深值為3.304 m；E6斷面位于河道順直段，水流較為平緩，水位較E5斷面變化不大，最小水深值為3.359 m，最大水深值為3.303 m；E7斷面位于彎道河段，河道發生較大的彎曲，慣性離心力的存在而使得自由水面的平衡狀態遭到破壞，水流在彎道前的一段距離內，凸岸處自由水面略高于凹岸[15]，E7斷面附近最小水深值為3.415 m，最大水深值為3.426 m；E8斷面位于河道最下游，由于河道變窄，流量最大，因此該斷面處的流速也最大，斷面水位具有一定的變化，最小水深值為3.428 m，最大水深值為3.434 m。

圖6 干流河道典型斷面VOF分布

本研究根據河道起點、河道交叉點、河道順直、河道急彎四個特點選取6個典型斷面分析河道的水流速度分布，為河道的水生態分析提供指導，見圖8。F1斷面位于河道起點，河道斷面較寬，寬度約29 m，由于該斷面處于河道起點，無外源水補充，可以發現河道水流流速比較慢，且起點為盲端，不受外界水流擾動，水流流態非常穩定，表現為低流速下的旋流，水體平均流速為0.24 m/s；F2斷面位于本研究河道與河道1的交叉口處，受河道補充水源的影響，該處水流表現為旋流、繞流、渦流、反射等紊流現象，表現為復雜的三維流動，水體平均流速為0.71 m/s；F3斷面位于本研究河道與河道2的交叉口處，水流流態與F2斷面類似，只是由于受河道2補水的影響，河道流量加大，水體流速變大，該處水體平均流速為0.83 m/s；F4斷面位于本研究河道與河道3的交叉口處，該處河道邊界條件非常復雜，一方面該處為兩條河道的交匯處，另一方面該處為急彎斷面，本研究河道的水流流向為自南向北，而河道3的水流流向為自北向南，所以兩條河道的流向相反，由圖8（F4）處的水流分析可以發現，河道水流整體還是按照本研究河道的流向流動，即自南向北，本研究河道有一部分水流流向河道3的下游，成為河道3的補充水源，而河道3的上游水流則主要匯入本研究河道，該處平均水流速度為0.84 m/s；F5斷面位于本研究河道順直段，河道中央水流較少受到外部邊界條件影響，水流流態較為穩定，但是由于河道寬度不規則變化，在河岸處河水表現為較明顯的紊流現象，并且受河岸阻力影響，水流流速較慢，該處平均水流速度為0.69 m/s；F6斷面位于本研究河道急彎處，水流受河道轉向的影響，在河岸處水體表現為明顯的紊流現象，而中央水流則隨河道轉彎，該處平均水流速度為0.76 m/s。

圖7 干流河道典型斷面水深分布（單位：m）

圖8 干流河道典型斷面的流速分布（單位：m/s）

5 結 論

流水不腐，是減緩甚至基本消除河流黑臭的關鍵因素，因此對河道水動力特性進行數值模擬研究對消除河流黑臭具有重要意義。針對目前河流水動力數值模擬的研究多集中于一、二維，忽略了河道水流在豎直方向上的運動，無法精確模擬河流自由水面形狀的復雜性和整體位置的不斷變化的問題，本文采用CFD技術，考慮水流自由水面以及豎直方向的運動，建立了基于VOF法的河流水力三維紊流數學模型，并采用S.Soare Frazao與Y.Zech在急轉彎河道模型中進行的試驗結果對本模型進行了驗證。以阜陽市某典型黑臭水體河道為研究對象，運用無結構貼體網格劃分技術和局部加密網格生成技術建立了河道的計算網格模型，對河段水動力特性進行了三維數值模擬研究。結果表明：在河道受其他河流補水后，水深逐漸增大，河道最小水深為2.79 m，最大水深3.43 m；在河道起點盲段，水流為緩慢的旋流；在河道交叉以及急彎處，水流表現為旋流、繞流、渦流、反射等復雜的三維紊流現象，并且水流速度較順直河段大，河道最小流速為0.24 m/s，最大流速為0.84 m/s。因此通過本文建立的三維紊流數學模型，實現了對自由水面、水深、流速的精確模擬，為河流水質的模擬提供基礎，為河道水生態治理與恢復提供科學的決策支持。