周正印，楊 楠

（天津市市政工程設(shè)計(jì)研究院，天津市 300392）

1 研究背景

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

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

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

2 研究方法

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

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

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

對(duì)于第q相流體，體積分?jǐn)?shù)連續(xù)方程為：

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

k-ε紊流模型的連續(xù)性方程、動(dòng)量方程、紊動(dòng)能k方程和ε方程如下：

連續(xù)性方程：

式中：t為時(shí)間，s；ui為速度分量，m/s；xi為坐標(biāo)分量；ρ為密度，kg/m3；μ 為分子粘性系數(shù)，N·m/s；P為修正壓力，Pa；ui為紊流粘性系數(shù)；CD為源項(xiàng)；k為紊動(dòng)動(dòng)能，m2/s2；ε 為紊動(dòng)耗散率，m2/s2；?k、?ε 分別為 k、ε 的紊流普朗特?cái)?shù)，無因次；C1ε、C2ε為經(jīng)驗(yàn)常數(shù)，無因次[10]。

2.2 模型的邊界條件

（1）進(jìn)口邊界條件：如果進(jìn)口條件已知，給定速度如下式，湍流參數(shù)和根據(jù)經(jīng)驗(yàn)公式計(jì)算得到；如果不知道準(zhǔn)確的來流條件，只知道來流的流量條件（Q為進(jìn)口流量），A為進(jìn)口面積，可以速度方向給定在結(jié)構(gòu)化網(wǎng)格方向上（如河道），并設(shè)m→為網(wǎng)格的方向矢量，則進(jìn)口速度給定如下：

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

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

（4）壓力邊界條件：在計(jì)算域的邊界上，壓力應(yīng)滿足Neumann條件，即第二類邊界條件[11]。為保證計(jì)算的穩(wěn)定性，在規(guī)定的某一內(nèi)點(diǎn)上，壓力為給定值，且定義為參考?jí)毫ref。

2.3 模型求解

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

3 模型驗(yàn)證

針對(duì)三維模型的特點(diǎn)，采用S.Soare Frazao與Y.Zech在急轉(zhuǎn)彎河道模型中進(jìn)行的試驗(yàn)結(jié)果作為驗(yàn)證的依據(jù)[14]。

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

圖1 潰壩實(shí)驗(yàn)?zāi)Ｐ褪疽鈭D(單位：m)

圖2 網(wǎng)格模型

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

4 工程應(yīng)用

4.1 工程簡(jiǎn)介

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

圖3 沿程水深模擬值與實(shí)驗(yàn)值對(duì)比(單位：m)

4.2 計(jì)算模型及網(wǎng)格劃分

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

根據(jù)該河道設(shè)計(jì)要求，確定該河道設(shè)計(jì)流量為10.40 m3，河道一設(shè)計(jì)流量為11.20 m3，河道二設(shè)計(jì)流量為22.10 m3，河道三設(shè)計(jì)流量為30.10 m3，水動(dòng)力計(jì)算采用進(jìn)口流速邊界，出口采用流速邊界，并保證總進(jìn)口流量等于總出口流量，河道基本達(dá)到穩(wěn)定狀態(tài)。

圖4 干流河道網(wǎng)格模型

為了能夠更好的分析河道水力條件，在河道干流設(shè)置了8個(gè)典型斷面監(jiān)測(cè)水流變量，分別在河道起始點(diǎn)、河道交叉處、河道中間、河道轉(zhuǎn)彎、河道末端等處，斷面詳細(xì)布置情況見圖5。

圖5 河道干流8個(gè)典型斷面分布

4.3 VOF動(dòng)態(tài)及水流流速分析

圖6為干流河道已基本達(dá)到穩(wěn)定狀態(tài)時(shí)，6個(gè)典型斷面的VOF分布情況。根據(jù)圖3模擬結(jié)果計(jì)算水深可得，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可以看出，隨著斷面從上游向下游推進(jìn)，河道水深逐漸變深，并且斷面水平面的波動(dòng)比上游劇烈。斷面1由于處于河道起點(diǎn)，且為盲端，沒有其他補(bǔ)充水源，水平面較為平穩(wěn)，EI斷面最小水深值為2.790 m，最大水深值為2.791 m；E2斷面位于本研究河道與河道1交叉處，河道斷面具有一定的波動(dòng)性，水深比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交叉口處，且該處河道轉(zhuǎn)彎，水流變化最為劇烈，水深增加較大，最小水深值為3.295 m，最大水深值為3.304 m；E6斷面位于河道順直段，水流較為平緩，水位較E5斷面變化不大，最小水深值為3.359 m，最大水深值為3.303 m；E7斷面位于彎道河段，河道發(fā)生較大的彎曲，慣性離心力的存在而使得自由水面的平衡狀態(tài)遭到破壞，水流在彎道前的一段距離內(nèi)，凸岸處自由水面略高于凹岸[15]，E7斷面附近最小水深值為3.415 m，最大水深值為3.426 m；E8斷面位于河道最下游，由于河道變窄，流量最大，因此該斷面處的流速也最大，斷面水位具有一定的變化，最小水深值為3.428 m，最大水深值為3.434 m。

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

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

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

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

5 結(jié) 論

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