丁壩挑角及長度對回流區(qū)流動特性影響的數(shù)值模擬

吳伊平，張 巖，崔鵬義，黃遠東，張紅武

（1.上海理工大學(xué) 環(huán)境與建筑學(xué)院，上海 200093；2.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084）

丁壩作為一種航道整治、防洪護坡的建筑物，廣泛應(yīng)用于我國江河和海岸線[1-2]。丁壩在岸邊迎托水流，把水流挑離河岸保護堤岸；約束水流，控制河流寬度，穩(wěn)固制導(dǎo)線；挑離水流進入主流加快流速，加速河流對底部河床的沖刷，增加行槽深度，提升河流輸沙能力[3]。丁壩的存在使得河流在丁壩附近的水流形態(tài)變得復(fù)雜。水流流經(jīng)丁壩時，在壩頭產(chǎn)生邊界層分離，形成復(fù)雜的三維旋渦結(jié)構(gòu)；在丁壩附近產(chǎn)生局部沖刷和淤積[4]。丁壩上游產(chǎn)生的壅水現(xiàn)象使得水面向丁壩對岸傾斜，而在丁壩斷面后水面普遍下凹。水流在丁壩壩頭處產(chǎn)生流動分離，大部分水流進入主流，少部分回流到丁壩后端，形成回流區(qū)[5]。

對回流區(qū)流動特性的探究是丁壩繞流問題的研究熱點，回流區(qū)域的大小對橫斷面上流速重新分布關(guān)系密切，正確預(yù)測回流區(qū)大小，了解丁壩的掩護范圍在實際工程中具有重要意義[6]。眾多學(xué)者通過水槽試驗、理論分析和數(shù)值模擬等對回流區(qū)域流動特性進行探究。李國斌等[7]從二維水深平均水流運動方程出發(fā)，推導(dǎo)出回流區(qū)長度和最大回流寬度計算公式，并與實測值進行了比較，但由于天然河道的復(fù)雜性，其公式的實際應(yīng)用仍需進一步驗證。陳稚聰?shù)萚8]利用三維激光顆粒動態(tài)分析儀（3D-PDA）對丁壩回流區(qū)的流速進行詳細測量，并根據(jù)流速分布和流量連續(xù)性原理對丁壩回流區(qū)進行詳細劃分，進一步研究縱向回流尺度和流量的變化規(guī)律。隨著計算機硬件及數(shù)值計算方法的快速發(fā)展，相比較水槽試驗較高的成本，數(shù)值模擬方法以其較低的成本、詳盡的預(yù)測數(shù)據(jù)，以及不斷提高的計算精度成為當前科學(xué)研究的重要工具。曹曉萌等[9]利用數(shù)值模擬和水槽試驗相結(jié)合探究丁壩群對河流系統(tǒng)3 個作用尺度的劃分準則。鄭艷等[3]利用RNGk-ε湍流模型，對非淹沒丁壩在寬淺河道中形成的流場進行了二維數(shù)值模擬，結(jié)果表明丁壩的相對長度和河寬束窄系數(shù)呈線性相關(guān)。

然而丁壩挑角及長度對回流區(qū)流動特性影響的研究還相對較少，實際流動中水流流向和丁壩不可能是理想的垂直狀態(tài)（即丁壩挑角為90 °），所以本文從丁壩挑角出發(fā)，結(jié)合丁壩長度變化，利用數(shù)值模擬的方法探究丁壩挑角和長度對回流區(qū)大小和水流紊動特性的影響。

1 數(shù)值模型及驗證

1.1 數(shù)值模型和方法

基于曹曉萌等[9]的水槽試驗?zāi)Ｐ停⒈狙芯康臄?shù)值模型，如圖1 所示?？紤]寬淺河道，X方向為水流方向，設(shè)計長度為20 m；Y方向為河寬方向，設(shè)計長度為1.2 m。壩體為非淹沒丁壩，其固定在模型下側(cè)，距離進口斷面6 m 處。丁壩投影長度b1（丁壩在垂直于水流方向的平面上的投影）分別為 0.12，0.2，0.3 m，在相同投影長度下分別設(shè)置挑角為60 °，90 °和120 °。設(shè)計進口斷面體積流量為0.061 8 m3/s，河道入口水深0.3 m。

圖1 數(shù)值模型二維平面圖Fig.1 Numerical model two-dimensional plan

根據(jù)圖1 給定的數(shù)值模型，設(shè)計計算區(qū)域大小為長20 m、寬1.2 m 和高0.3 m。利用Gambie2.4.6進行網(wǎng)格劃分，采用混合網(wǎng)格和局部加密技術(shù)對計算區(qū)域進行網(wǎng)格劃分。圖2（a）為其中一種工況下的平面二維網(wǎng)格示意圖（采用混合網(wǎng)格）。為更好地探究丁壩附近流場變化，對丁壩附近網(wǎng)格進行局部加密，如圖2（b）所示。最小網(wǎng)格尺寸在壩頭處為0.001 25，相鄰網(wǎng)格按照1∶1.04～ 1∶1.08 增長率遞增，各工況網(wǎng)格總數(shù)均在50～60 萬。

圖2 網(wǎng)格劃分示意圖Fig.2 Diagram of mesh generation

進口為速度入口，速度大小為0.172 m/s（由體積流量和設(shè)計水深確定）；出口為自由出流；模型的底面、側(cè)面和壩體表面都設(shè)置為速度無滑移的壁面邊界條件；上表面采用剛蓋假定，對應(yīng)頂部邊界設(shè)置為對稱邊界。采用Fluent 作為計算軟件，選用壓力基求解器，壓強和速度耦合采用Simple算法，對流項離散選用二階迎風格式。

1.2 控制方程

采用三維穩(wěn)態(tài)Navier-Stokes 方程來描述非淹沒丁壩繞流，其控制方程如下。

連續(xù)方程

動量方程

式中：xi，xj分別為i，j方向 上 的坐標；ui，uj分別為流體時均速度在i，j方向上的分量；ρ為水的密度；p為流體壓強；μ為分子粘滯系數(shù)，μt為湍流黏性系數(shù)；gi為重力加速度在i方向上的分量。

采用標準k-ε湍流模型[10]對方程組進行封閉，湍動能k和湍動能耗散率ε的輸運方程分別為

湍流黏性系數(shù)

湍動能產(chǎn)生項

式中：Cμ，σk，σε，C1ε，C2ε均為經(jīng)驗常數(shù)，其取值分別為0.09，1.0，1.3，1.44，1.92。

1.3 模型驗證

采用文獻[9]的試驗數(shù)據(jù)對本研究數(shù)值模型進行驗證。試驗中正挑丁壩長度為0.3 m，采用超聲波多普勒測速儀（ADV）對壩前0.2 m 和壩后0.4 m 處X方向上的流速Ux進行測量。驗證結(jié)果如圖3 所示，可見壩前壩后速度的數(shù)值模擬和試驗測量結(jié)果在數(shù)值和分布趨勢上都十分吻合，表明數(shù)值模型采用的相關(guān)參數(shù)可以合理地模擬非淹沒丁壩的繞流運動，因此采用上述模型設(shè)置對本文其他工況進行研究。

圖3 模型驗證Fig.3 Model validation

2 結(jié)果分析

2.1 流場和回流尺度

9 種工況下計算所得水面處速度云圖如圖4所示。從圖中可以看出，由于丁壩的挑流作用，當水流流經(jīng)丁壩時，主流收縮，流速V明顯增大，在丁壩下游形成明顯的順時針旋渦。從研究的不同工況可知：a.同一丁壩挑角下，隨著丁壩長度的增加，主流收縮程度增大，流經(jīng)丁壩的水流速度增大，并且達到最大收縮斷面的距離增加，沿程擴散的速度也增大，主流沿程擴散的距離增長；而在流線圖中水流流經(jīng)丁壩，在丁壩后形成不同尺度的順時針旋渦。b.同一投影長度下，隨著丁壩挑角的增加，壩后回流區(qū)產(chǎn)生的旋渦渦心向下游移動，但與壩頭之間的水平距離基本保持不變。c.同一投影長度下，不同挑角的丁壩，流速分布云圖和流線圖分布規(guī)律基本一致。

圖4 水面速度云圖及流線圖（丁壩長度-角度）Fig.4 Water surface velocity contours and streamlines

為研究丁壩挑角及投影長度對壩后回流區(qū)尺寸的影響，表1 列出了所研究工況壩后回流區(qū)的長度和寬度?？傮w來看，所有工況最大回流區(qū)長度均在1～ 3.5 m 之間，相同角度下，隨著投影長度的增加，回流區(qū)長度也增加，這和其他學(xué)者試驗測量結(jié)果一致。表1 給出了最大回流區(qū)長度L與丁壩投影長度b1 的比值L/b1，可以得出：丁壩壩后的回流區(qū)長度與丁壩的投影長度有關(guān)，且回流區(qū)長度一般是丁壩投影長度的11 倍；另外，隨著丁壩投影長度的增加，L/b1 基本不變。研究表明，丁壩相對投影長度b1/B與最大相對回流長度L/B成正相關(guān)，其線性關(guān)系為：L/B=1.910b1/B+1.226（R2=0.86，R為相關(guān)系數(shù)）。

另外，壩后最大回流寬度同樣隨著丁壩投影長度的增加而增大，且與丁壩的投影長度相關(guān)。一般壩后最大回流寬度是丁壩投影長度的1.7～2 倍左右，且隨著長度的增加，倍數(shù)略微減小。分析表明，丁壩的相對投影長度b1/B與最大相對回流寬度W/B成正相關(guān)，其線性關(guān)系為：W/B=0.274b1/B+0.206（R2=0.91）。

表1 所研究工況壩后回流區(qū)尺寸Tab.1 Dimension of recirculation zone behind spur dike in this study

2.2 湍動能分布

為研究丁壩設(shè)置對所產(chǎn)生的湍流混合作用強弱的影響，選取幾個不同丁壩、不同水深處的代表工況，對其湍動能分布云圖進行分析，如圖5所示。從流線圖可以看出，湍動能主要分布在回流渦區(qū)的下游，湍動能最大值出現(xiàn)在回流區(qū)末端與主流交匯處。水流流經(jīng)丁壩，在壩頭位置開始產(chǎn)生流動分離，湍動能產(chǎn)生并沿著流動分離區(qū)域向兩側(cè)及下游逐漸衰減。圖5（a）～（c）為3 種挑角下投影長度為0.3 m的丁壩水面的湍動能分布，可以看出湍動能分布區(qū)域在不同丁壩挑角下基本一致，回流區(qū)湍動能最大值與回流區(qū)渦心的相對位置基本保持不變；壩頭處隨著角度的增加湍動能最大 值從60 °的0.01 m2/s2減小為120 °的0.005 m2/s2，這是由于丁壩與來流流線方向的角度不同導(dǎo)致丁壩對主流挑流作用的強弱不同引起的。當丁壩角度小于等于90 °時，丁壩對主流流動起到抑制作用，此時流動分離更為劇烈，因此湍動能變大；當丁壩挑角大于90 °時，丁壩承接主流，隨著角度的增大，流動分離現(xiàn)象變?nèi)?，因此湍動能減小。圖5（b）、（d）、（e）為90 °丁壩、不同投影長度下、水面處的湍動能分布，可以明顯看出投影長度直接影響湍動能分布區(qū)域，投影長度越長，其湍動能影響范圍越大；且隨著投影長度的增加，湍動能最大值也增加，這是由于投影長度越長，壩頭處水流流速越大，流動分離越為劇烈，因此產(chǎn)生的湍動能越大。圖5（c）、（f）、（g）為120 °丁壩投影長度0.3 m 下、不同水深處的湍動能分布，可以看出，隨著深度的增加，湍動能分布規(guī)律相似，但湍動能的平均值減小，代表相同湍動能等值線范圍也減小。

圖5 湍動能分布圖（投影長度-角度-水深）Fig.5 Contours of turbulent kinetic energy (projection length-angle-water depth)

2.3 湍流黏度分析

湍流黏度μt是表征能量耗散度的量[11]，其表達式如式（5）所示，它是空間坐標的函數(shù)，表征流動的狀態(tài)，而非物性參數(shù)。提取9 種工況下不同深度的湍流黏度分布并作圖分析，選取具有代表性的工況如圖6 進行說明。從圖6 可知，湍流黏度分布區(qū)域與湍動能類似，不過湍流黏度最大值的位置靠近河岸。

圖6（a）～（c）為丁壩投影長度0.3 m、不同挑角情況下水面處的湍流黏度分布。相同投影長度下，不同角度回流區(qū)的湍流黏度影響區(qū)域基本一致，但分布規(guī)律存在差別。從云圖可以看出，隨著角度的增大，相同的等值線的范圍擴大，代表相同湍流黏度影響范圍增大，并且湍流黏度的最大值向回流區(qū)河岸靠近。另外，研究發(fā)現(xiàn)，不同投影長度下，其回流區(qū)湍流分布情況基本一致，研究選取90 °丁壩進行分析，如圖6（b）、（d）、（e）所示，從云圖可以明顯看出湍流黏度主要分布在回流區(qū)及其下游；并且當投影長度增加時，湍流黏度影響范圍增大，相應(yīng)湍流黏度最大值也隨之增大。圖6（c）、（f）、（g）為投影長度0.3 m、挑角120 °下不同水深處的湍流黏度分布，可以看出，不同水深處湍流黏度分布規(guī)律基本一致，在水面處（即水深為0）平均值最大，隨著水深的增加，相同大小湍流黏度的范圍逐漸縮小，但其最大值位置與渦心的相對位置基本保持不變。另外，研究發(fā)現(xiàn)，不同挑角、不同長度的丁壩沿水深方向的湍流黏度分布規(guī)律基本一致。

圖6 湍流黏度分布圖（投影長度-角度-水深）Fig.6 Contours of turbulent viscosity (projection length-angle-water depth)

3 結(jié)論

采用試驗驗驗證的數(shù)值模型研究了丁壩挑角及長度對回流區(qū)流動特性的影響，可以得到以下結(jié)論：

a.丁壩壩后的最大相對回流區(qū)長度L/b1 及相對寬度W/B與丁壩相對投影長度b1/B成正相關(guān)，線性關(guān)系為L/B=1.910b1/B+1.226（R2=0.86）及W/B=0.274b1/B+0.206（R2=0.91）。

b.湍動能最大值出現(xiàn)在壩后旋渦右上方的流動分離處，其距離壩頭的相對位置基本保持不變。具體表現(xiàn)為當投影長度為0.12 m 時，其距離壩頭約為0.9 m；投影長度為0.2 m 時，其距離壩頭約為1.6 m；投影長度為0.3 m 時，其距離壩頭約為2.4 m。湍動能產(chǎn)生后沿著流動分離區(qū)域向兩側(cè)及下游逐漸衰減。壩后湍動能的分布與丁壩的投影長度有關(guān)。

c.湍動能和湍流黏度分布主要集中在回流區(qū)末端和河流下游。丁壩的投影長度直接影響湍動能和湍流黏度的大小和范圍；丁壩挑角對壩頭湍動能影響較大，挑角越小壩頭湍動能越大；同一丁壩情況下湍動能和湍流黏度隨著水深的增加其分布規(guī)律不變，其值隨著水深的增加而小幅減小。