90°等寬明渠交匯水流污染物混合規(guī)律數(shù)值模擬

毛頌平，朱 海，王玲玲，華祖林，陸圣杰，薛曉鵬

(1.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098； 2.河海大學(xué)環(huán)境學(xué)院，江蘇 南京 210098；3.寧波市水利水電規(guī)劃設(shè)計(jì)研究院有限公司，浙江 寧波 315192)

交匯水流是天然河流及人工河道中普遍存在的水力現(xiàn)象。干支流交匯處水流紊動(dòng)劇烈，對(duì)污染物的混合特性有著重要影響。隨著水環(huán)境問(wèn)題的日益突出，國(guó)內(nèi)外學(xué)者對(duì)明渠交匯水流污染物混合規(guī)律進(jìn)行了大量研究。在物理試驗(yàn)研究方面，陳凱霖等[1]基于非等寬交匯口模型研究了污染物質(zhì)量濃度場(chǎng)的分布規(guī)律，探究了交匯角和污染物質(zhì)量濃度對(duì)其的影響；Ribeiro等[2]采用野外實(shí)測(cè)和水槽模型試驗(yàn)相結(jié)合的方式，研究發(fā)現(xiàn)干支流流量比對(duì)交匯口物質(zhì)運(yùn)輸有著重要影響，在流量比不變的情況下，拓寬支流寬度有利于污染物混合；Hua等[3]通過(guò)建立辮狀多支流河流物理模型研究污染物的混合特性，研究表明污染物的釋放位置是污染物遷移過(guò)程的主要影響因素；唐洪武等[4]通過(guò)實(shí)測(cè)發(fā)現(xiàn)懸沙質(zhì)量濃度與水溫更能反映長(zhǎng)江與鄱陽(yáng)湖交匯前后的水質(zhì)摻混特性。在數(shù)值模擬研究方面，劉雪蘭[5]應(yīng)用SMA二維模型模擬了重慶兩江段流場(chǎng)和水質(zhì)，發(fā)現(xiàn)在匯流口一側(cè)的河岸會(huì)形成高質(zhì)量濃度的污染帶；李娟等[6]建立二維水動(dòng)力水質(zhì)模型模擬汾河入黃口污染帶分布特性，研究表明橫向上擴(kuò)散起主導(dǎo)作用，縱向上對(duì)流起主導(dǎo)作用；米潭等[7]建立平面二維自由表面流MIKE 21 FM模型，研究了Y形交匯區(qū)污染物擴(kuò)散規(guī)律，得出了匯流比與污染物橫向、縱向質(zhì)量濃度梯度的關(guān)系；Lee等[8]建立二位對(duì)流彌散數(shù)學(xué)模型，研究了韓國(guó)漢江某分汊河道下游的水面形態(tài)與污染物輸運(yùn)過(guò)程；茅澤育等[9]建立了明渠交匯口水流流動(dòng)及污染物輸移模型，較好地模擬了污染混合區(qū)的大小；Biron等[10]分別采用試驗(yàn)?zāi)Ｐ秃吞烊缓恿髂M了交匯口污染物的混合，發(fā)現(xiàn)干支渠底部存在高差會(huì)加速污染物的混合；魏娟等[11]建立水氣兩相流數(shù)學(xué)模型，研究了不同匯流比工況下交匯口污染物輸移擴(kuò)散特性，得出匯流比越大，污染帶越狹長(zhǎng)的結(jié)論；Ramon等[12]通過(guò)三維數(shù)學(xué)模型分析了埃布羅河和塞格雷河交匯口污染物的混合速率；Tang等[13]采用三維數(shù)值模型研究發(fā)現(xiàn)污染物的混合受水流形態(tài)的強(qiáng)烈影響，而匯流比和河床形態(tài)又深刻影響著水流形態(tài)；顧莉等[14]通過(guò)建立U形交匯河道水氣兩相流模型研究發(fā)現(xiàn)污染物橫向離散系數(shù)沿程呈單峰結(jié)構(gòu)分布，而縱向離散系數(shù)呈雙峰結(jié)構(gòu)分布。

綜上可知，匯流比是影響交匯水流水動(dòng)力特性的重要因素，對(duì)交匯河口污染物輸運(yùn)及混合特性具有重要影響，目前針對(duì)90°等寬明渠交匯水流污染帶三維空間分布特征的定量研究仍有不足，因此本文基于三維水動(dòng)力-污染物耦合輸運(yùn)數(shù)學(xué)模型，著重研究了不同匯流比情況下明渠交匯口下游污染物三維空間分布特征，分析了污染物質(zhì)量濃度的混合特性，獲得了污染帶最大寬度與匯流比的定量關(guān)系。

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

1.1 控制方程

本文采用基于雷諾時(shí)均方程(RANS)的三維水氣兩相流模型，模擬研究區(qū)域污染物質(zhì)量濃度的分布。計(jì)算區(qū)域內(nèi)水流運(yùn)動(dòng)滿足質(zhì)量守恒方程及動(dòng)量守恒方程，污染物質(zhì)量濃度滿足質(zhì)量濃度守恒方程。本文采用雷諾應(yīng)力模型封閉控制質(zhì)量守恒方程和動(dòng)量守恒方程，根據(jù)時(shí)均化法則，直接建立并?；字Z應(yīng)力的輸運(yùn)方程，沒有采用渦黏性假設(shè)，可以有效地計(jì)算各向異性的湍流流場(chǎng)，對(duì)復(fù)雜流動(dòng)的模擬有著更強(qiáng)的優(yōu)勢(shì)。

雷諾應(yīng)力輸運(yùn)方程如下：

DT,ij+DL,ij+Pij+Πij+εij

(1)

該方程需聯(lián)合耗散率的輸運(yùn)方程進(jìn)行求解：

(2)

本文采用質(zhì)量濃度通量的二階矩模式封閉控制質(zhì)量濃度守恒方程：

(3)

1.2 數(shù)值方法及邊界條件

控制方程的離散方法采用有限體積法，控制方程組擴(kuò)散項(xiàng)的離散均采用中心差分格式，水動(dòng)力模塊控制方程對(duì)流項(xiàng)的離散采用QUICK格式，污染物輸運(yùn)模塊控制方程對(duì)流項(xiàng)的離散采用二階迎風(fēng)格式。采用PISO算法對(duì)壓力和速度場(chǎng)進(jìn)行耦合計(jì)算。

交匯水流的相互作用導(dǎo)致交匯口水面劇烈變化，而自由表面的模擬對(duì)于模型的可靠性與準(zhǔn)確性十分重要，本文采用VOF法模擬自由表面。計(jì)算區(qū)域底部及側(cè)壁采用無(wú)滑移邊界條件。主、支渠入口邊界采用速度進(jìn)口，流向垂直于進(jìn)口斷面。出口邊界采用壓力出口，下游出口水位保持與試驗(yàn)值一致，出口壓力滿足靜水壓強(qiáng)分布，液面相對(duì)壓強(qiáng)為0。進(jìn)出口邊界處紊動(dòng)能和耗散率的值均根據(jù)經(jīng)驗(yàn)公式[11]確定：k=0.003 75u2，ε=k1.5/(0.4H)，其中u為斷面平均速度，H為水深。

2 模型驗(yàn)證

2.1 水動(dòng)力模型驗(yàn)證

本文將水動(dòng)力模型數(shù)值模擬結(jié)果與Weber等[17]的90°等寬明渠交匯流物理試驗(yàn)數(shù)據(jù)作對(duì)比，驗(yàn)證工況下游出口尾水高度H0= 0.306 m，弗勞德數(shù)Fr=0.37，雷諾數(shù)Re=186 000，匯流比q=Qu/Qd= 0.25(Qu為主渠進(jìn)口流量，Qd為下游出口流量，Qd= 0.17 m3/s)。數(shù)學(xué)模型計(jì)算區(qū)域?yàn)榫匦螖嗝婷髑粎R口，主、支渠交匯角為90°，沿程等寬，底坡水平，采用順接方式連接，主、支渠適當(dāng)加長(zhǎng)是為了消除邊界條件的影響，使水流紊動(dòng)充分發(fā)展，計(jì)算區(qū)域如圖1所示。采用結(jié)構(gòu)化六面體網(wǎng)格剖分計(jì)算區(qū)域，對(duì)交匯口、近壁區(qū)域、水氣交界面和渠底網(wǎng)格進(jìn)行加密。經(jīng)過(guò)網(wǎng)格無(wú)關(guān)性分析，驗(yàn)證工況網(wǎng)格數(shù)為24.6萬(wàn)。

圖1 計(jì)算區(qū)域示意圖

數(shù)學(xué)模型弗勞德數(shù)、雷諾數(shù)及流量比均與物理試驗(yàn)保持一致，坐標(biāo)系以渠寬W進(jìn)行無(wú)量綱化處理，速度以下游出口平均流速U0進(jìn)行無(wú)量綱化處理。水面線模擬值與試驗(yàn)值對(duì)比結(jié)果表明，不同縱向斷面水面線高度與物理試驗(yàn)實(shí)測(cè)值吻合較好，交匯口附近不同測(cè)點(diǎn)的速度數(shù)值模擬值與物理試驗(yàn)實(shí)測(cè)值吻合良好，能夠反映流速變化的整體趨勢(shì)。經(jīng)過(guò)驗(yàn)證，該模型能夠較好地反映交匯口水流的水動(dòng)力特性。典型縱斷面水面線對(duì)比如圖2所示，典型斷面流速對(duì)比如圖3所示。

圖2 y/W=0.5縱斷面水面線對(duì)比

圖3 x/W=-1橫斷面各測(cè)線上的流速對(duì)比

2.2 污染物輸運(yùn)模型驗(yàn)證

為了驗(yàn)證數(shù)值模型對(duì)于污染物輸運(yùn)的可靠性，本文選取陳凱霖等[1,18]的物理試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比。研究對(duì)象主、支渠交匯角為90°，主渠寬度為W，支渠寬度為1/3W，主、支渠上游長(zhǎng)10W，主渠下游長(zhǎng)10W。計(jì)算區(qū)域采用六面體網(wǎng)格進(jìn)行剖分，總網(wǎng)格數(shù)為19.2萬(wàn)。計(jì)算工況Fr=0.10，Re=6 357.79，q=0.81，主流入口污染物質(zhì)量濃度為0 μg/L，支流入口污染物質(zhì)量濃度為2 000 μg/L，均與物理試驗(yàn)保持一致。圖4為驗(yàn)證工況下x/W=0.8和x/W=1.33橫斷面上y/W=0.5測(cè)線上流速實(shí)測(cè)值與計(jì)算值的對(duì)比，可見流速模擬值與測(cè)量值分布趨勢(shì)基本一致，由于本文采用的是雷諾時(shí)均模型，模擬獲得的流速分布曲線相對(duì)試驗(yàn)值更為光滑平順，總體誤差在10%以內(nèi)。圖5為驗(yàn)證工況下y/W=0.167和y/W=0.5縱斷面上z/W=0.427深度處污染物質(zhì)量濃度縱向分布計(jì)算值與實(shí)測(cè)值的對(duì)比，兩者分布趨勢(shì)吻合基本良好。

圖4 不同橫斷面y/W=0.5處流速對(duì)比

圖5 不同縱斷面z/W=0.427處污染物質(zhì)量濃度對(duì)比

3 不同匯流比條件下明渠交匯水流污染物三維分布規(guī)律

3.1 工況設(shè)置

為研究明渠交匯水流污染物質(zhì)量濃度的三維分布規(guī)律及匯流比對(duì)其的影響，采用驗(yàn)證后的模型對(duì)q=0.083、0.25、0.417、0.583、0.750、0.917等6種工況進(jìn)行計(jì)算，所選取的工況與Weber等[17]的90°等寬明渠交匯流物理試驗(yàn)選取工況一致。所有工況主、支渠等寬，交匯角均為90°，下游出口斷面Fr=0.37，Re=186 000，Qd=0.17 m3/s。主流入口污染物質(zhì)量濃度C=0 μg/L，支流入口污染物質(zhì)量濃度C0=1 μg/L。

3.2 污染物質(zhì)量濃度三維空間分布特征

圖6和圖7分別展示了匯流比為0.25和0.75時(shí)支流污染物進(jìn)入主渠后下游各橫斷面的質(zhì)量濃度分布。支流高質(zhì)量濃度污染物進(jìn)入主渠后，在主渠靠近支渠一側(cè)形成污染物高質(zhì)量濃度區(qū)域；污染物隨著支流的初始動(dòng)能沖向主渠右岸，沿主渠橫向斷面形成明顯的質(zhì)量濃度梯度；交匯口下游x/W>-3范圍內(nèi)污染物等質(zhì)量濃度線總體保持垂直，在垂向上質(zhì)量濃度變化不明顯；隨著干支流的進(jìn)一步混合，遠(yuǎn)離交匯口的下游區(qū)域(x/W<-3)等質(zhì)量濃度線逐漸發(fā)生傾斜和扭曲，污染物質(zhì)量濃度呈現(xiàn)明顯的三維特征。

圖6 q=0.25時(shí)交匯口下游典型橫斷面流速矢量及污染物質(zhì)量濃度分布云圖

圖7 q=0.75時(shí)交匯口下游典型橫斷面流速矢量及污染物質(zhì)量濃度分布云圖

模擬結(jié)果表明，不同匯流比條件下交匯口下游污染物空間分布特征具有顯著差異。不同匯流比條件下污染物橫向輸運(yùn)范圍不同，匯流比較小時(shí)，污染物橫向輸運(yùn)的范圍更大，匯流比較大時(shí)，污染物橫向輸運(yùn)的范圍較小，這主要是不同匯流比工況下支流進(jìn)入主渠的動(dòng)能不同導(dǎo)致的；不同水深處污染物質(zhì)量濃度橫向混合范圍也和匯流比相關(guān)，匯流比較小時(shí)，表層污染物橫向混合范圍大于底層污染物混合范圍(圖6)，污染物質(zhì)量濃度等值線呈D形分布；匯流比較大時(shí)，表層污染物橫向混合范圍小于底層污染物混合范圍(圖7)，污染物質(zhì)量濃度等值線呈L形分布。污染物質(zhì)量濃度分布特征與斷面環(huán)流緊密相關(guān)，在小匯流比工況下，由于在渠道中部存在順時(shí)針的斷面環(huán)流(圖6(c)(d))，靠近渠道左岸的污染物向表層水體輸運(yùn)，靠近渠道右岸的污染物向底層水體輸運(yùn)，等質(zhì)量濃度線在水體表層和底層向左岸移動(dòng)，在水體中部向右岸移動(dòng)；大匯流比工況下，由于在靠近渠道左岸處存在逆時(shí)針的斷面環(huán)流(圖7(c)(d))，靠近渠道左岸的污染物向底層水體輸運(yùn)，等質(zhì)量濃度線逐漸發(fā)生傾斜，逐漸形成L形的高質(zhì)量濃度區(qū)域。

3.3 污染物質(zhì)量濃度橫向分布范圍

污染帶的最大寬度表征著污染物橫向分布的最大范圍，對(duì)于研究污染物對(duì)水體的影響范圍有著重要意義。在一般情況下，當(dāng)邊緣點(diǎn)的質(zhì)量濃度為同一斷面上最大質(zhì)量濃度的5%時(shí)，該點(diǎn)被認(rèn)為是污染帶的邊界點(diǎn)[19]，因此污染帶最大寬度b定義為污染帶邊緣質(zhì)量濃度是斷面質(zhì)量濃度5%時(shí)所占有的最大寬度(如圖8(d)所示)。圖8為不同匯流比工況表層水體(z/W=0.3)污染物質(zhì)量濃度分布云圖，圖中紅色虛線為污染帶邊界。當(dāng)支流流量占比較大時(shí)(對(duì)應(yīng)q≤0.417)，表層水體污染帶總能在交匯口下游x/W=-2范圍內(nèi)影響到右岸(圖8(a)(b)(c))；對(duì)于支流流量占比較小時(shí)(q>0.417)，污染帶最大寬度b隨著匯流比的增大呈逐漸減小趨勢(shì)(圖8(d)(e)(f))。

圖8 不同匯流比條件下交匯口表層污染物質(zhì)量濃度分布云圖

數(shù)值計(jì)算結(jié)果表明，不同水深處污染帶寬度與匯流比具有不同的相關(guān)關(guān)系。圖9為各工況下不同水深處污染帶相對(duì)最大寬度b/W與匯流比的關(guān)系曲線。由圖9(a)可知，交匯口表層水體(z/W=0.3)的污染物在小匯流比工況下(如q<0.417)其橫向污染范圍已達(dá)到右岸，污染帶最大寬度不隨匯流比的變化而變化；在大匯流比工況(如q>0.417)，交匯口表層水體污染帶最大寬度隨著匯流比的增大而減小，呈現(xiàn)良好的一次函數(shù)關(guān)系。中下層水體污染帶最大寬度分布與表層水體不同，統(tǒng)計(jì)不同水深處各工況污染帶最大寬度，并作無(wú)量綱化處理，擬合的曲線表明b/W與匯流比在中下層水體中呈現(xiàn)相似的規(guī)律且有著良好的二次曲線關(guān)系，見圖9(b)。

圖9 各工況下不同水深污染帶b/W與q的關(guān)系曲線

3.4 污染物質(zhì)量濃度混合特性

通常用斷面污染物質(zhì)量濃度的不均勻系數(shù)δ及標(biāo)準(zhǔn)差σ的沿程分布評(píng)估污染物的混合速率[12]。對(duì)于一個(gè)給定的斷面，δ表征與下游完全混合后污染物質(zhì)量濃度的偏差，計(jì)算公式如下：

δ=(Ci-Cp)/Cp

(4)

式中：Ci為某斷面污染物平均質(zhì)量濃度計(jì)算值；Cp為污染物平均預(yù)測(cè)質(zhì)量濃度。δ越接近于0，污染物混合越充分。污染物的平均預(yù)測(cè)質(zhì)量濃度Cp，即充分混合之后河道內(nèi)污染物的質(zhì)量濃度可按下式計(jì)算：

Cp=(CmQm+CtQt)/(Qm+Qt)

(5)

式中：Cm和Ct分別為主流和支流的污染物質(zhì)量濃度；Qm和Qt分別為主流和支流的流量。

圖10為不同匯流比工況下交匯口下游各橫斷面δ的沿程分布。由圖10可見，各工況污染物不均勻系數(shù)的沿程變化呈現(xiàn)相似的變化規(guī)律：先增大，在x/W=-2斷面達(dá)到最大值后逐漸減小，減小的速率先快后慢，最后趨于穩(wěn)定。匯流比越大，δ的峰值越大，即斷面平均質(zhì)量濃度與平均預(yù)測(cè)質(zhì)量濃度的差距越大；匯流比越大，δ的下降速度也越大，即污染物混合的速率越快。

圖10 交匯口下游各橫斷面δ變化曲線

δ反映的是給定斷面污染物質(zhì)量濃度與完全混合后水體中污染物質(zhì)量濃度的偏差，而給定斷面污染物質(zhì)量濃度分布均勻程度通常用標(biāo)準(zhǔn)差σ來(lái)表征[12]。圖11為不同工況下交匯口下游主渠各橫斷面標(biāo)準(zhǔn)差σ的變化曲線。由圖11可知，干支流流量越接近(匯流比越接近0.5)，標(biāo)準(zhǔn)差越大，污染物混合越不充分。在大匯流比情況下(如q=0.583、0.750、0.917)，標(biāo)準(zhǔn)差呈現(xiàn)先迅速增大后緩慢減小的規(guī)律；在小匯流比工況下(如q=0.083、0.250、0.417)，標(biāo)準(zhǔn)差呈現(xiàn)先迅速增大后減小，再小幅增大，最后緩慢減小的規(guī)律。這是由于小匯流比工況下，主渠近支渠一側(cè)形成了明顯的分離區(qū)[20]。分離區(qū)的存在使得過(guò)流面積收縮，主流在經(jīng)過(guò)分離區(qū)時(shí)流速先增大后減小。因此，在x/W=-2斷面下游由于流速緩慢減小，污染物向下游輸運(yùn)的質(zhì)量濃度通量有所降低，從而減弱了污染物的混合程度，故標(biāo)準(zhǔn)差出現(xiàn)小幅增大。

圖11 交匯口下游各橫斷面σ變化曲線

4 結(jié) 論

a.交匯水流污染物分布三維特征明顯，匯流比對(duì)交匯水流污染物的分布有著重要影響，主要體現(xiàn)在污染物橫向分布上。匯流比較小時(shí)，污染物橫向輸運(yùn)的范圍較大，匯流比較大時(shí)，污染物橫向輸運(yùn)的范圍較小；隨著干支流的混合，交匯口下游斷面等質(zhì)量濃度線逐漸發(fā)生傾斜與扭曲。

b.污染物的分布與二次流密切相關(guān)。小匯流比工況下，渠道中部存在順時(shí)針的斷面環(huán)流，交匯口下游各斷面污染物質(zhì)量濃度總體呈中部質(zhì)量濃度高、表層及底層質(zhì)量濃度低的D形分布；大匯流比工況下，渠道中部存在逆時(shí)針的斷面環(huán)流，交匯口下游各斷面污染物質(zhì)量濃度總體呈底層質(zhì)量濃度高、中上層質(zhì)量濃度低的L形分布。

c.小匯流比工況下，表層水體污染帶最大寬度為主渠寬度，不隨匯流比的變化而改變，大匯流比工況下，污染帶最大寬度與匯流比呈現(xiàn)良好的一次函數(shù)關(guān)系；中下層水體污染帶最大寬度隨著匯流比的增大而減小，與匯流比呈二次曲線關(guān)系。

d.污染物的混合速率與混合程度均與匯流比密切相關(guān)。匯流比越大，污染物混合的速率越快，匯流比越接近0.5，污染物越不易充分混合；小匯流比情景下交匯口下游出現(xiàn)的大尺度分離區(qū)會(huì)在局部區(qū)域內(nèi)增加污染物質(zhì)量濃度分布的不均勻性。