弱混合條件下支流入?yún)R口水力特性數(shù)值模擬

陳雪妮，楊中華，尹 煒，吳 菲，白鳳朋

（1. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072； 2. 長江水資源保護科學(xué)研究所，湖北 武漢 430051； 3. 長江水利委員會湖庫水源地面源污染生態(tài)調(diào)控重點實驗室，湖北 武漢 430051； 4. 湖北工業(yè)大學(xué)土木建筑與環(huán)境學(xué)院，湖北 武漢 430068）

0 引 言

明渠交匯現(xiàn)象普遍存在于天然河流系統(tǒng)中［1］，交匯區(qū)域水動力特性復(fù)雜多變，干流受到支流頂托作用在交匯口上游發(fā)生壅水現(xiàn)象、流速降低，下游靠近支流側(cè)形成分離區(qū)，對污染物有滯留作用［2，3］；當支流來流量較小時，支流入?yún)R口形成大面積回流區(qū)域，造成水體滯留、泥沙沉降和污染物聚集［4］。因此明渠交匯口是影響水環(huán)境保護和航運整治的關(guān)鍵部位，交匯區(qū)水力特性的研究一直是水利工程設(shè)計和運行管理中的重點和難點。

國內(nèi)外不少學(xué)者針對干支流交匯區(qū)域水流特性及其變化規(guī)律展開了研究。Taylor［5］最早于1944 年對明渠交匯問題開展了試驗研究，明確提出了交匯口研究的必要性及重要性，為后續(xù)研究建立了基礎(chǔ)。此后，Best［6］通過水槽模型試驗提出將交匯河段劃分為滯留區(qū)、偏離區(qū)、分離區(qū)、加速區(qū)、剪切層和流動復(fù)原區(qū)六個區(qū)域。Yang［7］分別采用標準k-ε、RNGk-ε和RSM模型分析了不同湍流模型對分離區(qū)尺寸和二次流的影響，發(fā)現(xiàn)RNGk-ε模型模擬效果最好。Webber［8］等基于ADV 對90°交匯口三維水流結(jié)構(gòu)進行了詳細測量，為交匯口的數(shù)值模擬研究提供了大量實驗數(shù)據(jù)作為驗證支撐。茅澤育［9］通過模型試驗和理論分析發(fā)現(xiàn)交匯區(qū)域斷面存在二次環(huán)流，且分離區(qū)尺寸與匯流比有關(guān)。目前交匯區(qū)的研究成果多側(cè)重于交匯角和匯流比對匯流口及干流河段的影響作用［10-12］，近年來部分學(xué)者開始將干支流交匯段的研究重心向支流河道轉(zhuǎn)移，張琦［13］通過數(shù)值模擬對明渠交匯段的小型支流河道展開了研究，重點分析了不同彎曲形態(tài)和不同匯流比影響下的支流河口水面形態(tài)和泄流能力。

當支流來流量極小時，水體和污染物無法下泄進入干流河道，在適宜的光照和溫度條件下易造成庫區(qū)水體富營養(yǎng)化。近年來隨著對庫區(qū)支流富營養(yǎng)化問題研究的深入，水體流速、表層環(huán)流和紊動混合等水動力特征越來越多的成為研究富營養(yǎng)化問題的重要因子［14，15］。本文將這種支流入?yún)R不影響干流流場的交匯條件定義為弱混合條件，弱混合條件下交匯河段的交匯角和匯流比對支流入?yún)R口水力特性的影響機制尚不明晰。陸建宇［16］對河道型水庫庫區(qū)干支流交匯段進行概化，通過建立三維數(shù)學(xué)模型研究了天然河道尺度下不同交匯角對支流入?yún)R口紊動能分布和床面切應(yīng)力的作用機制，發(fā)現(xiàn)弱混合條件下紊動能峰值區(qū)域隨交匯角增加逐漸向支流右岸上層移動，交匯角＜90°時紊動能變化較小，交匯角＞90°時紊動能增長幅度較大。陸對弱混合條件下交匯角對支流紊動特性的影響較為細致，而匯流比對紊動能分布的影響以及交匯角和匯流比兩個因素影響的對比有待進一步研究。

基于此，本文以弱混合條件下交匯河段支流入?yún)R口水力特性為研究對象，采用室內(nèi)水槽試驗數(shù)據(jù)驗證后的RNGk-ε模型，對不同交匯角和不同匯流比的支流入?yún)R口水動力特性和紊動特性進行分析和研究，以期闡明交匯角和匯流比對支流入?yún)R口水力特性的影響。

1 數(shù)學(xué)模型

1.1 控制方程

連續(xù)性方程：

式中：ui為xi方向的時均速度；i=1，2，3。

動量方程：

式中：ρ為水的密度；p為時均壓強；t為時間；γi為水的動力黏性系數(shù)；fi為外部體積力。

1.2 RNG k-ε紊流模型

RNGk-ε模型對紊動黏度進行修正，從而引入平均流動中的漩流運動狀況，可以很好地反映應(yīng)變率高及流線彎旋狀態(tài)明顯的流體形態(tài)［17］。RNGk-ε模型的k及ε方程如下：

k方程：

ε方程：

式中：αk、αε為湍流普朗特數(shù)，均取1.39；μeff=μ+μt，，Cμ= 0.084 5；為紊動能平均速度梯度產(chǎn)生項；Gb為浮力產(chǎn)生項，不可壓縮流體Gb= 0；?；?shù)C1ε= 1.42，C2ε= 1.92，C3ε= 1.0；，η0=4.377，β= 0.012，；Sk、Sε為源項。

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

2.1 試驗?zāi)Ｐ图傲繙y系統(tǒng)

試驗在武漢大學(xué)水力學(xué)實驗室的矩形玻璃水槽中進行。干流水槽長8 m，寬0.3 m，高0.5 m，底坡為1/1 000，進口采用水泵供水并通過多道消浪柵平穩(wěn)水流，在水槽上游采用三角形薄壁量水堰測量流量。支流水槽長1.8 m，寬0.1 m，高0.45 m，支流采用管道供水，進口處鋪設(shè)鵝卵石平穩(wěn)水流，流量控制系統(tǒng)采用精度為0.5 級的電磁流量計。水槽出口處設(shè)有尾門以控制水位。

流速采用粒子圖像測速儀（Particle Image Velocimetry，簡稱PIV）進行測量，PIV 拍攝范圍為18 cm×14 cm，水槽可視化測量區(qū)域共80 cm，A～D 為相機拍攝位置，x-y平面測量Y1～Y4四個斷面，x-z平面測量S1～S3三個斷面，相機位置和斷面位置如圖1所示。

圖1 水槽試驗 PIV拍攝斷面布置示意圖（單位：cm）Fig.1 PIV shooting section layout diagram of experiment

2.2 計算區(qū)域及網(wǎng)格劃分

圖2 所示為交匯角α=90°工況的計算區(qū)域，原點位于交匯點上游，x、y、z軸正方向如圖所示。根據(jù)試驗水槽進行建模，干流寬0.3 m，水深h=0.15 m，支流長1.8 m，寬0.1 m，交匯口距干流入口3.2 m。選用六面體結(jié)構(gòu)化網(wǎng)格對計算區(qū)域進行剖分（圖3），網(wǎng)格數(shù)為345 600（α=90°）～398 340（α=150°），橫向最大網(wǎng)格尺寸為0.01 m，縱向為0.02 m，垂向為0.01 m。對交匯口附近網(wǎng)格進行加密。

圖2 計算區(qū)域示意圖（單位：cm）Fig.2 Computational area

圖3 計算網(wǎng)格示意圖Fig.3 Computational mesh

2.3 數(shù)值方法、邊界條件及計算工況

控制方程的離散方法采用有限體積法，速度和壓力解耦采用PISO 算法［18］，其他項離散均采用QUICK 格式，以RNGk-ε模型封閉雷諾時均方程中的紊動黏性系數(shù)。干、支流進口邊界采用速度進口，流向垂直于進口斷面，干流入口流速Um=0.22 m/s，支流入口流速va隨工況不同有所變化；出口邊界條件采用壓力出口，滿足靜壓假定；進出口的紊動能和耗散率邊界值均按照經(jīng)驗公式的計算結(jié)果初始化；計算區(qū)域底部及側(cè)壁采用無滑移邊界條件；根據(jù)水槽模型試驗觀察，低匯流比條件下支流水面比降較小，故模型頂部采用剛蓋假定。

為了研究交匯角和匯流比對支流流場的影響，設(shè)計6 種計算工況，具體工況設(shè)置如表1所示。

表1 數(shù)值模擬計算工況Tab.1 Case setting of numerical simulation

2.4 模型驗證

選取水槽模型試驗工況中交匯角α=90°、匯流比RQ=0.25（干支流流量分別為Q干流=10 L/s、Q支流=2.5 L/s）的測量結(jié)果對模型結(jié)果進行驗證，其中模型邊界條件與相應(yīng)試驗邊界參數(shù)一致。以支流寬度W為特征長度對橫縱坐標進行無量綱化得到x*、y*，無量綱化水深z*=z/h，將時均流速U除以干流入口流速Um得到無量綱化U*，下文中出現(xiàn)的物理量若無特別說明均經(jīng)過上述量綱歸一化處理。圖4 所示為交匯口附近區(qū)域近水面處（z*=0.93）的x-y平面流場矢量圖與實測結(jié)果的對比結(jié)果，可以看出支流入?yún)R后干流交匯區(qū)產(chǎn)生水流偏轉(zhuǎn)現(xiàn)象，交匯區(qū)下游形成分離區(qū)，數(shù)值模擬可以較好反映出交匯口附近流速偏轉(zhuǎn)情況、分離區(qū)范圍及二次流現(xiàn)象。圖5為不同位置處的無量綱化流速分量Ux、Uz沿垂向分布情況與水槽試驗實測值Vx、Vz的對比結(jié)果，可以看出模擬結(jié)果與實測資料吻合較好。

圖4 近水面處流場矢量圖與實測結(jié)果的對比結(jié)果Fig.4 Comparison of velocity vector near the water surface with the measured results

圖5 不同位置時均流速垂向分布與實測數(shù)據(jù)的對比結(jié)果Fig.5 Comparison of the vertical distribution of time-averaged flow velocity at different positions with the measured results

3 支流入?yún)R口水力特性

3.1 支流入?yún)R口水流結(jié)構(gòu)分區(qū)

當天然河道中支流來流較小時，支流水體由于干流的阻礙作用在入?yún)R口右岸產(chǎn)生大范圍回流，將產(chǎn)生回流的區(qū)域定義為回水區(qū)，回水區(qū)的存在造成水面收縮從而導(dǎo)致左岸水流加速形成加速區(qū)［圖6（a）］。對于回水區(qū)尺寸的界定，如圖6（b）所示，參考分離區(qū)尺寸界定的等值線法［8，19］，即順支流向流速為零的等值線為回水區(qū)邊界，零等值線區(qū)域的最大長度和最大寬度定義為回水區(qū)的長度L和寬度B。

圖6 支流入?yún)R口水流結(jié)構(gòu)分區(qū)Fig.6 Tributary confluence flow structure zones

3.1.1 不同水深斷面支流回水區(qū)形狀

交匯角α=90°、匯流比RQ=0.005 條件下回水區(qū)長度L、寬度B在不同水深斷面的分布規(guī)律如圖7 所示，L和B沿水深方向整體呈現(xiàn)出先減小后增加的變化規(guī)律。近底面附近（z*＜0.2）回水區(qū)形狀為“半葫蘆”型（圖8），L和B為全水深斷面最大值，說明近底面處入?yún)R口回流范圍較大，原因可能為此時支流流速較?。ń咏?），受到干流阻礙作用后更易發(fā)生流向偏轉(zhuǎn)。隨著水深斷面上移，L和B都逐漸減小；回水區(qū)形狀和寬度在z*=0.4～0.85 范圍內(nèi)基本不變，表明中上層水體回水區(qū)和加速區(qū)維持動態(tài)平衡；最后在近水面處（z*＞0.85）L和B都出現(xiàn)較大幅度增長，回水區(qū)形狀轉(zhuǎn)變?yōu)椤般^”型，這是由于回水區(qū)在表層位置流速較大的緣故。

圖7 支流回水區(qū)形狀參數(shù)沿垂線分布Fig.7 Vertical distribution of tributary backwater zone shape parameters

圖8 支流回水區(qū)形狀示意圖Fig.8 Schematic diagram of the tributary backwater zone

3.1.2 交匯角和匯流比影響下入?yún)R口分區(qū)特性

圖9 為不同交匯條件下支流入?yún)R口近水面處（z*=0.93）x-y平面流速云圖。α=30°時，回水區(qū)范圍從左岸向支流上游延伸至右岸，回流流速峰值出現(xiàn)在交匯口右岸處；α=90°時，右岸回水區(qū)和左岸加速區(qū)分區(qū)明顯且形狀規(guī)則，流速峰值較30°時有所減??；當交匯角增加至150°時，回流范圍基本覆蓋了整個支流入?yún)R口，回水區(qū)長度L長達8.66W，加速區(qū)緊貼左岸呈細長三角形形狀。隨著支流來流量增加［圖9（b）、（d）、（e）、（f）］，匯流比增大，回水區(qū)和加速區(qū)范圍逐漸縮小，流速峰值也隨之減小，當RQ增加至0.01 時，回水區(qū)和加速區(qū)消失，支流入?yún)R基本不再受到干流的影響，表層水體恢復(fù)至正常流動。

圖9 不同工況下入?yún)R口x-y平面流速云圖Fig.9 Flow velocity clouds under different cases in x-y plane

3.2 匯流比影響下入?yún)R口縱向時均流速分布

圖10 為α=90°時不同匯流比下不同位置處縱向時均流速沿垂線分布圖（其中v/va是以支流進口流速va為基準的無量綱化縱向流速）?？梢钥闯?，4種匯流比下不同橫截面處左岸流速分布基本遵循對數(shù)分布規(guī)律，且匯流比越小加速區(qū)越明顯；隨著x增加，支流右岸流速逐漸減小并產(chǎn)生回流現(xiàn)象。近交匯口（y*=-0.2）截面上的回流現(xiàn)象僅出現(xiàn)在近右岸壁面處。沿支流向上游至y*=-1.0 截面處，加速區(qū)內(nèi)流速增加，RQ≤0.008 工況下回水區(qū)范圍沿橫向擴大至x*=0.8處，RQ=0.01時右岸流速分布恢復(fù)至對數(shù)分布，表明y*=-1.0 截面上水體基本不再受到干流影響。至y*=-2橫截面，除RQ=0.005 時仍有小范圍回流，其余工況下支流均恢復(fù)正常流動。

圖10 不同匯流比下不同位置縱向流速沿垂線分布Fig.10 Vertical distribution of longitudinal velocity with different convergence ratios

3.3 交匯角和匯流比影響下入?yún)R口紊動特性

水流的紊動特性是影響藻類生長的重要水動力因子，而紊動能是衡量交匯區(qū)紊動強度的重要參數(shù)，本節(jié)對弱混合條件下支流入?yún)R口紊動能分布進行了分析研究。綜合3.1 節(jié)和3.2 節(jié)可知，l*＜0.2 時入?yún)R口分區(qū)不明顯（其中l(wèi)*=l/W，為支流斷面至交匯口的無量綱化距離），l*＞0.2 時A5 工況下入?yún)R口流速恢復(fù)正常，故選取斷面M（l*=0.2）分析交匯口水流的紊動能分布（圖11）。α=30°時，斷面M右岸水面處有極小的紊動能，水流摻混均發(fā)生在l*＜0.2 近交匯口范圍內(nèi)。隨著交匯角增加，紊動能峰值區(qū)域逐漸向支流上游移動，且峰值區(qū)域面積和峰值大小顯著增大。在不同匯流比條件下，RQ=0.005 時紊動能峰值出現(xiàn)在斷面M右側(cè)中部，峰值大小約為3.29×10-4m2/s2，表明干支流摻混主要發(fā)生在支流右岸。隨著匯流比增加，紊動能峰值區(qū)域逐漸減小且向水面移動。RQ=0.01 時，整個斷面紊動能都較小，說明此時干支流摻混區(qū)域已逐漸脫離支流。

圖11 不同工況下斷面M紊動能分布Fig.11 Turbulent kinetic energy distribution of cross-section M under different cases

圖12 為不同交匯條件下斷面M縱向流速分布，正值區(qū)域（紅色）表示支流水體入?yún)R干流，負值區(qū)域（藍色）表示水體倒灌支流。交匯角α=30°時，回流主要發(fā)生在水面處；隨著交匯角增加，回流區(qū)逐漸轉(zhuǎn)移至右岸全水深斷面；α=150°時，回流區(qū)占據(jù)右側(cè)大部分區(qū)域，干支流摻混最為劇烈，產(chǎn)生紊動能峰值區(qū)域［圖11（c）］。如圖所示，4種匯流比工況下均表現(xiàn)為左側(cè)入?yún)R右側(cè)倒灌，回流流速峰值均出現(xiàn)在右岸壁面，表明倒灌水體從右岸進入支流后與來流水體形成摻混，導(dǎo)致附近流速梯度增加形成紊動能峰值區(qū)；隨著匯流比增加，負值區(qū)域范圍減小且逐漸上移，表明水體倒灌也隨之向右岸上層發(fā)展。

圖12 不同工況下斷面M的縱向流速分布Fig.12 Longitudinal velocity distribution of cross-section M under different cases

4 結(jié) 論

通過數(shù)值模擬對不同交匯角和不同匯流比下支流入?yún)R口的回水區(qū)范圍、紊動能和流速分布進行了研究，得到了低匯流比下干流對支流入?yún)R口水流特性和紊動特性的影響隨交匯角和匯流比的變化情況。

（1） 弱混合條件下支流入?yún)R口右岸產(chǎn)生回流形成回水區(qū)，左側(cè)水體擠壓形成加速區(qū)。

（2） 支流回水區(qū)長度和形狀變化主要受到交匯角和匯流比兩種因素影響，不同交匯角條件下回水區(qū)長度大小為：150°＞30°＞90°；交匯角α=90°時回水區(qū)長度隨匯流比RQ的增加而減小，RQ=0.01時回水區(qū)基本消失。

（3） 回流流速峰值出現(xiàn)在近右岸壁面處，匯流比RQ越小，回水區(qū)最大寬度的橫截面位置距交匯口越遠；交匯角α=90°時支流入?yún)R口受影響區(qū)域在y*＞-2范圍內(nèi)。

（4） 在不同交匯角和匯流比條件下均存在紊動能峰值區(qū)域，區(qū)域面積及峰值大小與匯流比RQ成反比，與交匯角α成正比；交匯角對紊動能的影響較匯流比更大。