淹沒(méi)壩田水流特性和物質(zhì)交換規(guī)律研究

武夢(mèng)爽，楊中華，向 珂，2，白鳳朋，3，張 為

（1.武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430072；2.清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084；3.長(zhǎng)江水資源保護(hù)科學(xué)研究所，湖北武漢 430051）

0 引言

丁壩廣泛應(yīng)用于航道整治和河道生態(tài)修復(fù)工程中，具有增加河床深度，改善航道條件和水生生境的作用。近年來(lái)，隨著生態(tài)文明建設(shè)的加強(qiáng)，丁壩的生態(tài)效應(yīng)受到廣泛的關(guān)注，在大型河流整治工程中，利用淹沒(méi)丁壩群來(lái)恢復(fù)河流退化的魚(yú)類(lèi)棲息地，增強(qiáng)河流生態(tài)系統(tǒng)的多樣性［1］。在丁壩群中，相鄰丁壩間形成壩田，壩田內(nèi)流速較低，水體在壩田滯留。在壩田內(nèi)，懸浮顆粒物的局部停留時(shí)間增加，加速浮游植物的生長(zhǎng)和繁殖，為植被的生長(zhǎng)提供條件。壩田內(nèi)的水流特性是影響棲息地地貌特征、生物種群分布規(guī)律的關(guān)鍵因素，對(duì)適宜生境的形成有著重要作用［2-4］。國(guó)內(nèi)外學(xué)者圍繞壩田的水流結(jié)構(gòu)和物質(zhì)交換規(guī)律開(kāi)展了大量的研究。許多學(xué)者通過(guò)試驗(yàn)和數(shù)模研究發(fā)現(xiàn)非淹沒(méi)情況下壩田內(nèi)的流動(dòng)可以視為一種水平二維運(yùn)動(dòng)，其環(huán)流模式與壩田的體型密切相關(guān)［5-7］；并且通過(guò)對(duì)壩田和主渠之間的混合層水動(dòng)力特性的研究發(fā)現(xiàn)，混合層內(nèi)的擬序結(jié)構(gòu)控制著壩田和主渠間的物質(zhì)交換［8-10］?；诖?，部分學(xué)者開(kāi)展了對(duì)淹沒(méi)壩田的研究：McCoy 等［6，11］利用大渦模擬比較了淹沒(méi)和非淹沒(méi)情況下單個(gè)壩田的水流特性、渦結(jié)構(gòu)及物質(zhì)交換，發(fā)現(xiàn)淹沒(méi)情況下壩田及其周?chē)鷧^(qū)域的水流運(yùn)動(dòng)和渦系更加復(fù)雜，橫向混合層和垂向混合層的發(fā)展及其驅(qū)動(dòng)壩田和主渠進(jìn)行物質(zhì)交換的作用機(jī)制是相似的，但淹沒(méi)狀態(tài)下的交換速率更快；Yossef等［12］通過(guò)水槽試驗(yàn)重點(diǎn)關(guān)注了連續(xù)淹沒(méi)壩田二維水流結(jié)構(gòu)，并對(duì)壩田環(huán)流模式、混合層水動(dòng)力特性和紊動(dòng)特性的參數(shù)化方法進(jìn)行了研究。劉易莊等［13］通過(guò)水槽實(shí)驗(yàn)對(duì)淹沒(méi)單一壩田不同間距的水流結(jié)構(gòu)進(jìn)行了研究，結(jié)果表明丁壩間距對(duì)壩田內(nèi)渦旋中心位置及渦量分布有著顯著影響。

隨著生態(tài)保護(hù)需求的提升，非淹沒(méi)高丁壩在工程中的應(yīng)用越來(lái)越少，工程中常采用高淹沒(méi)度的低矮的潛壩［14］，在起到河勢(shì)調(diào)整作用的同時(shí)避免對(duì)河流生境造成過(guò)大的擾動(dòng)。已有關(guān)于淹沒(méi)壩田的研究集中于低淹沒(méi)狀態(tài)，淹沒(méi)度（水深∕丁壩高度）均在1.5 以內(nèi)，并且淹沒(méi)度對(duì)壩田水流特性和物質(zhì)交換過(guò)程的影響機(jī)制并不明晰，因此本文通過(guò)數(shù)值模擬研究了單一壩田不同淹沒(méi)度下壩田水流特性和物質(zhì)交換規(guī)律，為河道生態(tài)修復(fù)工程的設(shè)計(jì)提供一定的理論參考。

1 模型建立及驗(yàn)證

1.1 模型建立

本文采用劉易莊等［13］的試驗(yàn)為算例對(duì)模型進(jìn)行驗(yàn)證，試驗(yàn)水槽長(zhǎng)16 m，寬0.4 m，高0.5 m。驗(yàn)證工況參數(shù)為丁壩長(zhǎng)W=0.05 m，厚B=0.01 m，高h(yuǎn)g=0.18 m，丁壩間距L=0.15 m，所形成的壩田寬長(zhǎng)比為λ=W∕L=0.33，是航道整治工程中常用的布置方式［12，15］，上游丁壩設(shè)在x=0 m位置處。流量Q=57.6 m3∕h，進(jìn)口流速U=0.20 m∕s，水深H=0.20 m。根據(jù)試驗(yàn)工況參數(shù)最終確定計(jì)算域?yàn)? m×0.4 m×0.3 m（如圖1 所示），其中水和空氣的組分比為2∶1。網(wǎng)格布置如圖2 所示，采用結(jié)構(gòu)網(wǎng)格對(duì)計(jì)算域進(jìn)行劃分，并對(duì)壩田附近網(wǎng)格進(jìn)行局部加密，所用總網(wǎng)格數(shù)為346 536個(gè)。

圖1 數(shù)值幾何模型（單位：m）Fig.1 Numerical geometric model

圖2 計(jì)算網(wǎng)格布置圖Fig.2 Computational mesh

1.2 邊界條件及求解方法

整個(gè)計(jì)算區(qū)域設(shè)置初始水位0.2 m。邊界條件為進(jìn)口設(shè)空氣和水流兩個(gè)進(jìn)口，均設(shè)置為速度進(jìn)口，流速為0.2 m∕s；出口為自由出流，計(jì)算域頂部按剛蓋假定的方法處理，其他設(shè)為固壁邊界。

基于商業(yè)求解器Fluent（版本14.5）對(duì)單個(gè)壩田水流特性和物質(zhì)交換規(guī)律進(jìn)行研究。Fluent使用有限體積法對(duì)連續(xù)性方程和雷諾時(shí)均N-S 方程進(jìn)行離散，并采用RNG k-ε 模型來(lái)封閉方程組。與標(biāo)準(zhǔn)k-ε 模型相比，RNG k-ε 模型能夠更好地預(yù)測(cè)分離流、二次流和旋流等復(fù)雜流動(dòng)，近年來(lái)，許多學(xué)者采用RNG k-ε模型研究了丁壩附近的流動(dòng)［16-19］。

自由面的追蹤采用隱式VOF 方案，近壁流動(dòng)采用標(biāo)準(zhǔn)的壁面函數(shù)法。采用PISO算法作為壓力-速度耦合算法求解離散的動(dòng)量方程，壓力差值的離散采用PRESTO 格式，體積分?jǐn)?shù)、湍動(dòng)能及耗散率的離散采用二階迎風(fēng)格式［20］，收斂標(biāo)準(zhǔn)為計(jì)算各變量殘差小于0.000 1。

1.3 計(jì)算結(jié)果分析

試驗(yàn)采用PIV 系統(tǒng)測(cè)量了3 個(gè)不同相對(duì)水深平面的流場(chǎng)，圖3 給出了數(shù)模和試驗(yàn)結(jié)果的流場(chǎng)圖對(duì)比，可以看出，在z=1∕6hg平面和z=1∕2hg平面上，存在一個(gè)順時(shí)針渦旋，隨著相對(duì)水深的增大，渦旋尺度減小，在z=5∕6hg平面上，由于受到壩頂?shù)挠绊?，流態(tài)混亂，不易形成渦旋。數(shù)模計(jì)算的流速場(chǎng)與試驗(yàn)結(jié)果吻合較好。圖4 展示了數(shù)模結(jié)果和試驗(yàn)數(shù)據(jù)的流速對(duì)比；在z=1∕6hg平面上沿流向在壩田內(nèi)的上游、中游以及下游位置選取了對(duì)比位置（即線s1-s1、s2-s2及s3-s3），并對(duì)試驗(yàn)數(shù)據(jù)附上了10%的誤差線?？梢钥闯鰯?shù)模結(jié)果與試驗(yàn)數(shù)據(jù)具有一致性，但由于數(shù)值模型未考慮液體的表面張力，并且試驗(yàn)可能存在一定的測(cè)量誤差，因此部分區(qū)域流速對(duì)比誤差大于10%。結(jié)合流場(chǎng)和流速對(duì)比來(lái)看，數(shù)模的整體流場(chǎng)和流速分布均接近試驗(yàn)結(jié)果，說(shuō)明本文所建立的數(shù)值模型可以較好地模擬淹沒(méi)雙丁壩附近的流場(chǎng)。基于此模型，進(jìn)行本文工況的模擬。

圖3 數(shù)模與試驗(yàn)的流場(chǎng)對(duì)比Fig.3 Velocity vector comparison between experiment and numerical model

圖4 數(shù)模與試驗(yàn)的時(shí)均流向速度對(duì)比Fig.4 Time-averaged flow velocity comparison between experiment and numerical model

2 數(shù)值模擬分析

2.1 工況設(shè)置

淹沒(méi)度不同，水流運(yùn)動(dòng)的復(fù)雜性不同。已有研究的淹沒(méi)度較小，但實(shí)際工程中的淹沒(méi)丁壩常處于高淹沒(méi)狀態(tài)，因此本文選取了從低淹沒(méi)狀態(tài)到高淹沒(méi)狀態(tài)的四個(gè)淹沒(méi)度（分別為1.11、1.6、2、2.5）進(jìn)行研究。丁壩體型、布置形式及水流條件與劉易莊等［13］的試驗(yàn)保持一致，具體工況設(shè)置見(jiàn)表1。

表1 工況設(shè)置Tab.1 Case setting

2.2 計(jì)算結(jié)果分析

2.2.1 環(huán)流結(jié)構(gòu)

圖5展示了不同高程平面（z=0.1hg；z=0.5hg；z=0.9hg）上的二維時(shí)均流線，圖中，，m∕s?？梢钥闯觯蜎](méi)條件下壩田水流的三維性很強(qiáng)，流線特性沿水深發(fā)生變化。

圖5 x-y平面時(shí)均流線圖Fig.5 Mean velocity streamlines in x-y planes

在近底面（z=0.1hg），淹沒(méi)度為1.11 時(shí)，壩田內(nèi)存在著明顯的橫向環(huán)流，旋心在壩頭連接線上，并且靠近下游丁壩；當(dāng)淹沒(méi)度增大時(shí)，環(huán)流結(jié)構(gòu)發(fā)生改變，下游丁壩壩頭附近存在一個(gè)小尺度橫向環(huán)流。在一半丁壩高度平面（z=0.5hg），4 種淹沒(méi)度下均存在單環(huán)流結(jié)構(gòu)，環(huán)流幾乎占據(jù)了整個(gè)壩田，并且環(huán)流中心隨著淹沒(méi)度的增大向主渠方向靠近。在壩頂附近平面（z=0.9hg），淹沒(méi)度為1.11時(shí)，下游丁壩和側(cè)壁的角落存在一個(gè)橫向環(huán)流；淹沒(méi)度為1.6和2時(shí)，壩田內(nèi)存在較為相似的環(huán)流結(jié)構(gòu)，環(huán)流中心均靠近上游丁壩；淹沒(méi)度為2.5 時(shí)，壩田處于高淹沒(méi)狀態(tài)，壩頂附近平面的流線相對(duì)平行，不存在橫向環(huán)流。

圖6 展示了不同縱剖面（y=0.2W；y=0.5W；y=0.8W）上的二維時(shí)均流線。對(duì)比不同工況下同一剖面上的流線可知，隨著淹沒(méi)度的增大，環(huán)流占據(jù)壩田的空間越大。淹沒(méi)度為1.11 時(shí)，近壁面（y=0.2W）和一半丁壩長(zhǎng)度平面（y=0.5W）上的環(huán)流結(jié)構(gòu)基本一致，上游丁壩壩頂附近存在垂向環(huán)流，并且由于自由面的影響，下游丁壩后方的回流越過(guò)下游丁壩壩頂進(jìn)入壩田；但在交界面附近（y=0.8W），壩田內(nèi)呈現(xiàn)雙環(huán)流結(jié)構(gòu)，壩頂附近和壩田中下部均存在垂向環(huán)流，兩個(gè)環(huán)流尺度相當(dāng)。淹沒(méi)度為1.6和2時(shí)，流線特性沿橫向變化不大，壩田內(nèi)存在單環(huán)流結(jié)構(gòu)，環(huán)流中心靠近壩頭連線和下游丁壩。淹沒(méi)度為2.5 時(shí)，不同縱剖面上壩田內(nèi)的二維流線幾乎是一樣的，垂向環(huán)流基本占據(jù)了整個(gè)壩田。

圖6 x-z平面時(shí)均流線圖Fig.6 Mean velocity streamlines in x-z planes

2.2.2 流速分布

為了綜合反映壩田與主渠之間流速的橫向和垂向分布，沿流向截取5 個(gè)截面，位置選取如圖7 所示，參考Sukhodolov［5］的研究，引入整體平均的思想，對(duì)時(shí)均流速在空間上進(jìn)行平均得到整體平均流速。在每個(gè)截面上對(duì)時(shí)均流向速度沿壩田高度（即丁壩高）進(jìn)行平均得到整體平均流速的橫向分布；沿壩田寬度（即丁壩長(zhǎng)）進(jìn)行平均得到整體平均流速的垂向分布。

圖7 截面位置示意圖Fig.7 Sketch of the selected positions

圖8 給出了4 種工況下不同截面上整體平均流向速度的橫向分布。對(duì)比不同截面的流速分布，可知所有工況下，壩田平均流速和主渠平均流速之間的差值沿流向先增后減，在x=0.7L截面（即d-d截面）達(dá)到最大。對(duì)比不同淹沒(méi)度下的流速橫向分布可知，隨著淹沒(méi)度的升高，丁壩對(duì)水流的阻隔作用減弱，主渠平均流速減小，因此壩田內(nèi)平均流速和主渠平均流速之間的差值相應(yīng)地減小?？偟膩?lái)看，淹沒(méi)度對(duì)整體平均流速的橫向分布的影響主要體現(xiàn)在流速大小而不影響其分布形狀。

圖8 整體平均流速橫向分布Fig.8 Ensemble-averaged streamwise velocity along the y-direction

圖9 給出了4 種工況下不同截面上整體平均流向速度的垂向分布。當(dāng)淹沒(méi)度h/hg=1.1 時(shí)，由于自由水面距壩頂?shù)木嚯x較小，受到自由水面的限制，流速的垂向分布與其他工況下的明顯不同，壩田和主渠的流速差沿流向不斷縮小，并且在靠近下游丁壩的截面上，壩頂以上區(qū)域時(shí)均流速為負(fù)值，這是因?yàn)橄掠味魏蠓降牟糠只亓魉髟竭^(guò)壩頂進(jìn)入壩田［如圖6（a）、（b）所示］，這種現(xiàn)象隨著淹沒(méi)度的升高消失。當(dāng)淹沒(méi)度為1.6、2 和2.5 時(shí)，壩田內(nèi)的回流速度沿流向先增后減，在x=0.7L截面（即d-d截面）達(dá)到最大。隨著淹沒(méi)度的升高，自由水面與壩頂?shù)木嚯x增加，靠近自由表面的流速沿流向逐漸趨于均勻。

圖9 整體平均流速垂向分布Fig.9 Ensemble-averaged streamwise velocity along the z-direction

2.2.3 混合層厚度

由于壩田與主渠之間存在明顯的速度差異，在交界面位置會(huì)形成混合層；一般采用混合層厚度δ來(lái)描述環(huán)流與主流之間的大流速梯度范圍及其對(duì)應(yīng)的擴(kuò)散角。橫向混合層厚度δl按照以下公式定義［8］：

式中：Um(x)和Ug(x)分別為主渠和壩田流向坐標(biāo)x處的整體平均流向速度；Um(x)和Ug(x)對(duì)應(yīng)的y軸坐標(biāo)是速度梯度絕對(duì)值小到可以忽略的位置，參照Mignot 等［8］的研究，本文選取0.5 s-1作為流速梯度閾值為交界面上各流向位置x所對(duì)應(yīng)的最大速度橫向梯度。

類(lèi)比橫向混合層厚度的定義式，將垂向混合層厚度δv按照以下公式定義：

圖10 展示了4 個(gè)工況下橫向混合層和垂向混合層厚度沿流向發(fā)展。對(duì)于橫向混合層，當(dāng)x<0.5L時(shí)，混合層厚度整體上沿流向呈現(xiàn)線性增加的規(guī)律；當(dāng)x>0.5L時(shí)，混合層厚度沿流向的增長(zhǎng)率放緩。隨著淹沒(méi)度增大，壩田和主渠之間的速度差不斷減小，因此混合層厚度及其增長(zhǎng)斜率也減小。對(duì)于垂向混合層，需要說(shuō)明的是，當(dāng)淹沒(méi)度為1.1 時(shí)，自由液面限制了垂向混合層的發(fā)展，混合層厚度遠(yuǎn)小于其他3種工況，并且受到下游丁壩后方回流的影響，在壩田下游已經(jīng)不存在垂向混合層，因此圖中只展示了其他3 種工況下的垂向混合層厚度沿流向的演變。3個(gè)工況下，混合層厚度沿流向不斷增加；并且隨著淹沒(méi)度的升高，混合層厚度相應(yīng)地增加。

圖10 不同工況下混合層厚度沿流向的演變Fig.10 Streamwise evolution of the mixing layer thicknesses for different cases

淹沒(méi)度升高，混合層的水動(dòng)力特性發(fā)生改變?；旌蠈雍穸仍酱?，混合層內(nèi)擬序結(jié)構(gòu)典型尺寸越大，通過(guò)混合層的質(zhì)量交換效率越高［21］。在本文模擬工況參數(shù)下，隨著淹沒(méi)度的升高，橫向混合層厚度減小，垂向混合層厚度增加。低淹沒(méi)狀態(tài)下（h∕hg=1.11），橫向混合層的厚度大于垂向混合層，因此橫向混合層在物質(zhì)交換過(guò)程中處于主導(dǎo)地位，橫向混合層中的擬序結(jié)構(gòu)控制主渠和壩田之間的質(zhì)量和動(dòng)量交換。隨著淹沒(méi)度的升高，垂向混合層發(fā)展，在高淹沒(méi)狀態(tài)下（h∕hg=2.5），垂向混合層厚度大于橫向混合層，此時(shí)垂向混合層在物質(zhì)交換過(guò)程中處于主導(dǎo)地位。

2.2.4 交換系數(shù)

經(jīng)典的壩田物質(zhì)交換理論認(rèn)為，壩田與主渠之間的質(zhì)量交換效率與兩區(qū)域之間的物質(zhì)濃度差異成正比［7］，即：

式中：Cg為壩田內(nèi)的濃度；t為時(shí)間；Cm為主渠的濃度；K為質(zhì)量交換系數(shù)，s-1，是定量描述壩田與主渠之間混合過(guò)程的常用指標(biāo)。

對(duì)于有初始濃度的壩田，主槽水流的卷吸作用會(huì)將物質(zhì)帶出壩田，壩田濃度下降梯度與剩余濃度成正比，比例系數(shù)即為交換系數(shù)。因此在初始時(shí)刻在壩田引入濃度為C0的物質(zhì)，并假設(shè)主渠濃度為0，對(duì)式（3）進(jìn)行求解，可以得到：

式（4）說(shuō)明了交換系數(shù)K代表著壩田和主渠之間的交換速率。

基于已建立的水動(dòng)力模型，在水流充分發(fā)展后，利用Flunet中的UDS 功能在壩田內(nèi)瞬時(shí)均勻地引入標(biāo)量物質(zhì)，此時(shí)設(shè)定為初始時(shí)刻，壩田內(nèi)無(wú)量綱標(biāo)量初始物質(zhì)濃度為1，主渠濃度設(shè)定為0，監(jiān)測(cè)壩田內(nèi)空間平均濃度隨時(shí)間變化過(guò)程直至99%的物質(zhì)離開(kāi)壩田，繪制壩田濃度（體積平均濃度）隨時(shí)間變化曲線（對(duì)數(shù)坐標(biāo)），曲線的斜率即為質(zhì)量交換系數(shù)K［6，9］。

圖11 給出了壩田平均濃度隨時(shí)間變化曲線。由圖11 可知，壩田內(nèi)的濃度以指數(shù)型下降，壩田和主渠之間的交換存在兩個(gè)階段。當(dāng)壩田內(nèi)的濃度降至50%～40%，交換系數(shù)發(fā)生改變。這一現(xiàn)象與McCoy 等［11］的研究具有一致性。第一個(gè)階段物質(zhì)濃度快速下降，交換系數(shù)大，這由通過(guò)混合層的直接交換控制的；第二個(gè)階段濃度下降速率放緩，交換系數(shù)減小，這是由于物質(zhì)位于環(huán)流中心，此時(shí)混合層內(nèi)的大尺度結(jié)構(gòu)無(wú)法滲透到環(huán)流中心。

圖11 壩田平均濃度隨時(shí)間變化過(guò)程Fig.11 Variation process of volume average concentration in groyne field with time

對(duì)比4 個(gè)工況下的交換系數(shù)，可以發(fā)現(xiàn)第一個(gè)階段的交換系數(shù)K1隨著淹沒(méi)度的升高而增大。當(dāng)淹沒(méi)度為2.5 時(shí)，混合層整體厚度達(dá)到最大，混合層內(nèi)的渦體侵占了壩田內(nèi)更多的水體，因此淹沒(méi)度為2.5 時(shí)交換效率最高。第二個(gè)階段的交換系數(shù)K2隨著淹沒(méi)度的升高先增后減，在淹沒(méi)度為2時(shí)達(dá)到最大。

3 結(jié)論

通過(guò)數(shù)值模擬對(duì)不同淹沒(méi)度下壩田的環(huán)流結(jié)構(gòu)、流速分布、混合層厚度和交換系數(shù)進(jìn)行了研究，得到了淹沒(méi)度對(duì)壩田水流特性和物質(zhì)交換的影響。

（1）淹沒(méi)條件下壩田內(nèi)存在橫向環(huán)流和垂向環(huán)流；淹沒(méi)度對(duì)流速的橫向分布的影響主要體現(xiàn)在流速數(shù)值大小而不影響其分布形狀。

（2）橫向混合層的厚度隨著淹沒(méi)度的升高而減小，而垂向混合層的厚度隨著淹沒(méi)度的升高而增大。低淹沒(méi)狀態(tài)下，壩田主要通過(guò)橫向混合層與主渠進(jìn)行交換，高淹沒(méi)狀態(tài)下，壩田主要通過(guò)垂向混合層與主渠進(jìn)行交換。

（3）不同淹沒(méi)度下，壩田和主渠之間的交換均存在兩個(gè)階段，第一個(gè)階段交換速率大，交換系數(shù)隨著淹沒(méi)度升高而降低，第二個(gè)階段交換速率放緩，交換系數(shù)隨著淹沒(méi)度升高先增后減。