武夢爽，楊中華，向 珂，2，白鳳朋，3，張 為

（1.武漢大學水資源與水電工程科學國家重點實驗室，湖北武漢 430072；2.清華大學水沙科學與水利水電工程國家重點實驗室，北京 100084；3.長江水資源保護科學研究所，湖北武漢 430051）

0 引言

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

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

1 模型建立及驗證

1.1 模型建立

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

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

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

1.2 邊界條件及求解方法

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

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

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

1.3 計算結果分析

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

圖3 數(shù)模與試驗的流場對比Fig.3 Velocity vector comparison between experiment and numerical model

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

2 數(shù)值模擬分析

2.1 工況設置

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

表1 工況設置Tab.1 Case setting

2.2 計算結果分析

2.2.1 環(huán)流結構

圖5展示了不同高程平面（z=0.1hg；z=0.5hg；z=0.9hg）上的二維時均流線，圖中，，m∕s。可以看出，淹沒條件下壩田水流的三維性很強，流線特性沿水深發(fā)生變化。

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

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

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

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

2.2.2 流速分布

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

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

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

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

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

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

2.2.3 混合層厚度

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

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

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

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

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

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

2.2.4 交換系數(shù)

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

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

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

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

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

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

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

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

3 結論

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

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

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

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