尹則高,王 盛,孫 燦

(1.中國海洋大學(xué)工程學(xué)院,山東 青島 266100；2.浙江禹貢信息科技有限公司,浙江 杭州 311265)

射流是高濃度鹽水和廢水排放等工業(yè)領(lǐng)域常見的物理現(xiàn)象。根據(jù)射流流體與環(huán)境流體物理特性(密度、溫度等)的差異,一般可將射流分為純射流和浮力射流[1]。純射流流體與環(huán)境流體物理特性一致,動力特性相對簡單,對其研究較早且較完善。前人通過理論、實驗和數(shù)值模擬的方法,得到了靜水條件下純射流的時均流速沿程分布、射流核長度、射流擴(kuò)展率等相關(guān)經(jīng)驗公式,并得出了各斷面的流速分布自相似且符合高斯分布的結(jié)論；與靜水條件相比,波浪作用下純射流與環(huán)境水體間相互作用明顯增強(qiáng),射流擴(kuò)展寬度和紊動強(qiáng)度大幅增加,射流稀釋能力顯著提高[2-7]。然而,自然條件下射流流體與周圍環(huán)境水體常存在密度差,形成復(fù)雜的浮力射流。靜水條件下,諸多學(xué)者對浮力射流時均、紊動特性進(jìn)行了實驗研究和理論分析,得到了射流最大上升高度、穩(wěn)態(tài)高度、撞擊位置和各特征位置稀釋度等經(jīng)驗公式[8-13]。波浪條件下,其時間尺度與射流時間尺度相近,波浪引起水質(zhì)點往復(fù)運動會對浮力射流特性產(chǎn)生顯著影響[14-15]。

本文建立規(guī)則波作用下垂直圓管負(fù)浮力射流的三維數(shù)學(xué)模型,著重分析波參數(shù)對射流中軸速度衰減和射流速度半寬的影響規(guī)律,探討不同波-射流動量比下負(fù)浮力射流的振蕩擴(kuò)散機(jī)制。建立量綱一平均速度衰減系數(shù)、平均特征半寬經(jīng)驗公式,以期為工程實際提供參考。

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

1.1 控制方程

三維波浪運動的控制方程為RANS形式。RNGK-ε模型考慮了紊流漩渦和浮力效應(yīng)的影響,對負(fù)浮力射流的紊流特征具有較高的模擬精度。本研究中,射流流體與環(huán)境水體存在一定密度差,對常密度運動流體的連續(xù)性方程進(jìn)行改寫,得到考慮變密度影響的連續(xù)性方程:

(1)

動量方程:

(2)

紊動動能(K)方程:

(3)

紊動耗散率(ε)方程:

(4)

式中:αK和αε為紊流普朗特數(shù)；ueff為有效動力黏度系數(shù)；GK為由層流平均速度梯度而產(chǎn)生的紊流動能；Gb為浮力產(chǎn)生的紊流動能；YM為可壓縮修正系項；C1ε、C2ε和C3ε為經(jīng)驗常數(shù),分別取1.44、1.92和0.20。

利用VOF方法追蹤水氣界面的實時變化。

(5)

式中:F=F(x1,x2,x3,t)為水體體積分?jǐn)?shù)。當(dāng)F=1時,網(wǎng)格充滿液體；當(dāng)F=0時,網(wǎng)格充滿空氣；當(dāng)0

1.2 邊界條件

在造波邊界定義自由水面高度和水質(zhì)點的速度來生成規(guī)則波,其波面高度變化(η)可以用下式表示:

(6)

水質(zhì)點水平速度:

(7)

水質(zhì)點垂直速度:

(8)

式中:H為波高,m；T為波周期,s；k為波數(shù),m-1；d為靜水深,m；ω為波圓頻率,s-1。

數(shù)值水槽末端采用孔隙率為0.7的多孔介質(zhì)材料,減少波浪反射。利用質(zhì)量動量法模擬負(fù)浮力射流源項；通過定義射流源的位置、形狀和方向,求解器在源上自動生成這些粒子。使用隱式粒子-流體耦合法計算從時間步長N到N+1的粒子速度,近似的動量方程為:

(9)

(10)

求解方程(9)和(10),得

(11)

2 物理模型試驗及數(shù)學(xué)模型驗證

2.1 物理模型試驗

為驗證所建立數(shù)學(xué)模型的可靠性,在中國海洋大學(xué)水動力學(xué)實驗室中的波浪水槽進(jìn)行了試驗,水槽長25 m,寬1.0 m,高1.2 m,上述長寬高分別對應(yīng)x,y,z方向；水槽前端為推板造波機(jī),末端為多孔介質(zhì)消波材料。射流裝置布置在水槽中部,垂直噴口內(nèi)徑D=0.015 m,距離水槽底部h=0.10 m,射流出口中心坐標(biāo)為(0,0,0)。水箱內(nèi)鹽水密度ρj=1 050 kg/m3(Dm=1.23×10-9m2/s),經(jīng)進(jìn)水管垂直射入密度ρa(bǔ)=1 000 kg/m3的清水中。利用進(jìn)水管上安裝的CN-8300型水泵提供并調(diào)節(jié)射流流量,利用進(jìn)水管上安裝的HR-LWGB-10型渦輪流量計監(jiān)測流量。采用小威龍聲學(xué)多普勒流速儀(ADV)測量波浪作用下射流縱剖面上的三維速度,采用波高儀觀測波面過程。本文試驗和數(shù)學(xué)模型靜水深d=0.6 m。物理模型試驗及數(shù)學(xué)模型示意如圖1所示。

圖1 物理模型試驗及數(shù)學(xué)模型示意Fig.1 Sketch of experimental and mathematical model

2.2 網(wǎng)格獨立性驗證

本文利用FLOW-3D軟件建立并求解上述數(shù)學(xué)模型。網(wǎng)格劃分中,選用理查森外推法(Richardson extrapolation)計算離散誤差,并用網(wǎng)格收斂指數(shù)(IGC)方法測試網(wǎng)格獨立性。這里選取了3種不同尺寸的遠(yuǎn)場網(wǎng)格,S1=0.01 m,S2=0.02 m,S3=0.04 m；為了精確捕捉射流及波面附近區(qū)域的水動力特征,分別對上述3種遠(yuǎn)場網(wǎng)格劃分下的波面和射流附近局部區(qū)域進(jìn)行加密,加密網(wǎng)格尺寸分別為0.002 5 m,0.005 m和0.01 m。表1給出了波高H=0.10 m、波周期T=1.5 s、射流口直徑D=0.015 m、射流口速度wj=1.57 m/s、t=17.8 s時距離造波端5 m處波面水質(zhì)點速度、波面高度及計算誤差?？梢钥闯?網(wǎng)格1和網(wǎng)格2的離散誤差較小,分別為0.7%和0.2%。圖2(a)為距離造波端5 m 處的波面高度時程曲線?？梢钥闯鼍W(wǎng)格3與其他2種網(wǎng)格劃分的計算結(jié)果出現(xiàn)了相對較大的誤差,但是網(wǎng)格1和網(wǎng)格2之間的結(jié)果誤差較小。圖2(b)給出了t=20 s時沿水槽方向波面高度的空間分布,分別計算網(wǎng)格1、網(wǎng)格2和網(wǎng)格3劃分條件下x=8 m和x=10.99 m處的波高,得到單位距離內(nèi)相對波高衰減值分別為0.2%、0.31%和0.5%,相對較小,說明3種尺寸的網(wǎng)格劃分不會對計算波面高度產(chǎn)生大的影響。綜合考慮計算精度和時間的平衡,選擇網(wǎng)格2參與后續(xù)數(shù)值計算。

表1 3種網(wǎng)格下的誤差大小Table 1 Discretization error using three grids

圖2 3種尺寸網(wǎng)格的計算波面高度曲線Fig.2 Computational wave elevation using three grids

2.3 三維波浪水槽驗證

為驗證建立的數(shù)值波浪水槽準(zhǔn)確性,首先對無射流條件下波浪水槽進(jìn)行數(shù)值模擬計算。圖3給出了波高H=0.10 m、波周期T=1.5 s條件下,距離波浪水槽造波端9 m處穩(wěn)定后的規(guī)則波波面時程曲線和理論值的對比?？梢钥闯?數(shù)值計算波面與理論波面變化趨勢比較吻合,相對誤差均在5%以內(nèi)。數(shù)值計算中波高出現(xiàn)輕微衰減的現(xiàn)象,這是由于數(shù)學(xué)模型水體存在一定的黏性損耗,本文建立的三維數(shù)值波浪水槽是可靠的。

圖3 波面時程曲線驗證Fig.3 Wave surface elevation history validation

2.4 負(fù)浮力射流數(shù)值模擬的驗證

圖4給出了靜水和規(guī)則波(H=0.10 m,T=1.5 s,D=0.015 m,wj=1.57 m/s)條件下射流中軸(經(jīng)過射流口中心點的垂線)上的量綱一平均垂向速度(w0/wj)計算值與試驗值,其中利用相位平均法得到射流中軸上某一點的平均垂向速度w0。發(fā)現(xiàn)隨著射流距離的增加,平均垂向速度相應(yīng)降低,計算值與試驗數(shù)據(jù)吻合較好,說明本文建立的規(guī)則波作用下負(fù)浮力射流模型合理可靠,可以進(jìn)行后續(xù)的計算工作。

圖4 射流中軸上量綱一平均垂向速度驗證Fig.4 Verification of the dimensionless average vertical velocity at the jet axis

3 計算結(jié)果與分析

基于驗證后的數(shù)學(xué)模型,對不同規(guī)則波參數(shù)下負(fù)浮力射流的水動力行為進(jìn)行計算,其中射流口直徑D=0.015 m,射流口速度wj=1.57 m/s,ρj=1 050 kg/m3,ρa(bǔ)=1 000 kg/m3。工況1為靜水無波浪條件,其他工況波高和波周期組合分別為:工況2(0.06 m,1.5 s)、工況3(0.08 m,1.5 s)、工況4(0.10 m,1.5 s)、工況5(0.12 m,1.5 s)、工況6(0.10 m，1.2 s)、工況7(0.10 m,1.8 s)、工況8(0.10 m，2.1 s)。

3.1 射流濃度場和速度場分布

圖5給出了T=1.5 s、H=0.10 m工況4射流口中心所在縱剖面(y=0)上4個典型時刻的濃度場和速度場分布?？梢钥闯?上述4個時刻,遠(yuǎn)離波面的射流口附近流向基本呈現(xiàn)向上狀態(tài),且流速較大。規(guī)則波上跨零點時刻(t=T/4),受波面上升水質(zhì)點影響,距離射流口較遠(yuǎn)處的流向基本呈現(xiàn)向上狀態(tài),入射波帶動射流流體向上運動。波峰時刻(t=T/2),波面水質(zhì)點運動方向與波浪傳播方向一致,距離射流口較遠(yuǎn)處的流速基本偏向正x方向；從入射波上跨零點到波峰期間,波面附近水質(zhì)點運動帶動射流軸向正x方向擺動。下跨零點時刻(t=3T/4),受波面下降水質(zhì)點影響,距離射流口較遠(yuǎn)處的流向基本呈現(xiàn)垂直向下狀態(tài),帶動射流頂端云團(tuán)向下運動。波谷時刻(t=T),波面水質(zhì)點運動方向與波浪傳播方向相反,距離射流口較遠(yuǎn)的流向基本偏向負(fù)x方向；從下跨零點到波谷期間,波面附近水質(zhì)點流向的改變帶動射流軸向負(fù)x方向擺動?？傮w來看,波浪運動帶動了射流軸的擺動,加快了射流流體與環(huán)境水體的摻混,有助于射流流體稀釋和擴(kuò)散。

圖5 y=0縱剖面上典型時刻水體速度場和濃度場分布Fig.5 Velocity field and concentration field distribution at typical times under y=0 longitudinal section

3.2 波高對射流特性的影響

圖6(a)給出了T=1.5 s時,不同波高條件下射流中軸上量綱一平均垂向速度(w0/wj)與量綱一射流距離(z/D)的關(guān)系曲線?？梢钥闯?隨著波高增加,w0/wj相應(yīng)降低,射流中軸上平均垂向速度w0衰減加快。從圖中還可以看出,z/D<10時,w0/wj衰減較快；z/D>10時,射流初始動量的影響相對較小,負(fù)浮力和波浪的影響逐漸增加,其衰減速度趨緩,可以用經(jīng)驗公式(12)[15]來描述。

(12)

式中:B和n是衰減系數(shù)和指數(shù),分別表示射流中軸上量綱一平均垂向速度衰減的早晚和衰減率。隨著B值增加,射流中軸上任一點平均垂向速度衰減程度降低,即衰減越晚；隨著n值增加,射流中軸上任兩點的平均垂向速度差值減小,即衰減率減小。

圖6(b)給出了圖6(a)工況下z/D>10范圍內(nèi)的w0/wj數(shù)據(jù)；為了與純水射流試驗進(jìn)行對比,這里引用了文獻(xiàn)[16]在規(guī)則波作用下觀測的試驗數(shù)據(jù)。基于公式(12),對文獻(xiàn)[16]試驗數(shù)據(jù)和本文計算值擬合得到的衰減指數(shù)n分別為-1.07和-1.65～-1.75,本文的n值變化不大,說明波高對平均垂向速度沿射流中軸的衰減率影響不大。由于負(fù)浮力加快了射流軸上垂向速度衰減,本文w0/wj小于文獻(xiàn)[16]的試驗數(shù)據(jù)。另外,相應(yīng)擬合線相關(guān)系數(shù)較高,說明w0/wj與z/D的關(guān)聯(lián)度較好。為分析方便,對本文計算值進(jìn)行擬合時n取-1.65和-1.75的均值-1.70,得到0.06 m、0.08 m、0.10 m和0.12 m波高條件下對應(yīng)的B值分別為19.41、16.42、13.39和10.41,上述B值迅速減小,差值也較大,說明波高對平均垂向速度沿射流中軸衰減的早晚有較大影響；波高越大,w0/wj衰減越早。

圖6 不同波高條件下射流中軸上量綱一平均垂向速度分布Fig.6 Dimensionless average vertical velocity distribution along the jet axis under various wave heights

圖7給出了T=1.5 s時,距離射流口不同距離水平剖面的縱軸線上量綱一平均垂向速度(w/w0)分布,其中w為某一點的平均垂向速度?？梢钥闯?距離射流口較近時,其分布基本符合高斯分布,具有一定的自相似性。隨著遠(yuǎn)離射流口,w/w0分布范圍逐漸變寬。波高從H=0.08 m增加到0.10 m和0.12 m時,z值較大剖面上的w/w0開始呈現(xiàn)雙峰分布,這是由于較強(qiáng)的波浪振蕩引起射流中軸較大幅度的擺動引起的。隨著波高增加,射流擴(kuò)展率相應(yīng)增加,開始出現(xiàn)雙峰的水平剖面也更加靠近射流口。z值較大的水平剖面上,遠(yuǎn)離射流中軸的w/w0出現(xiàn)了負(fù)值,這可能是此處距離射流口較遠(yuǎn),射流初始動量對其影響相對較小,而射流鹽水密度大于清水造成的負(fù)浮力影響較大造成的。

圖7 距離射流口不同距離的水平剖面縱軸線上量綱一平均垂向速度分布Fig.7 Dimensionless average vertical velocities along the longitudinal line at horizontal sections with different distances from the jet exit

圖8為T=1.5 s、不同波高條件下量綱一射流特征速度半寬(bw/D)與量綱一垂向距離(z/D)的關(guān)系,這里bw定義為水平剖面縱軸線上w0/wj=0.5所對應(yīng)的2個x值之差的一半?？梢钥闯?隨著z/D的增加,靜水和波浪條件下的bw/D均相應(yīng)增加；靜水條件下的bw/D與z/D基本呈現(xiàn)線性關(guān)系。z/D<5時,靜水和波浪條件下的bw基本一致,原因是在這段區(qū)域內(nèi)射流主要受射流初始動量影響,波浪和負(fù)浮力對其影響較小；z/D>5時,隨著波高增加,bw和射流擴(kuò)展率相應(yīng)增加,這是由于射流初始動量影響沿程衰減,波浪和負(fù)浮力對射流影響逐漸增加造成的。

圖8 不同波高條件下射流特征速度半寬與量綱一垂向距離關(guān)系Fig.8 Half width of jet velocity relationships with dimensionless vertical distance under various wave heights

3.3 波周期對射流特性的影響

圖9(a)為H=0.10 m、不同波周期條件下射流中軸上量綱一平均垂向速度(w0/wj)的變化?？梢钥闯?隨著波周期增加,w0/wj相應(yīng)降低,射流中軸上垂向速度衰減越明顯；z/D<10時,w0/wj衰減較快；z/D>10時,其衰減速度趨緩。圖9(b)給出了圖9(a)工況下z/D>10范圍內(nèi)的w0/wj數(shù)據(jù),為了與純水射流試驗進(jìn)行對比,同樣引用了文獻(xiàn)[16]在規(guī)則波作用下觀測的試驗數(shù)據(jù)?；诠?12),分別對文獻(xiàn)[16]試驗數(shù)據(jù)和本文計算值進(jìn)行擬合,得到n分別為-1.07和-1.61～-1.80,本文的n值變化不大,說明波周期對平均垂向速度沿射流中軸的衰減率影響不大。另外,相應(yīng)擬合線相關(guān)系數(shù)較高,說明w0/wj與z/D的關(guān)聯(lián)度較好,這也與圖6(b)給出的趨勢相吻合。為分析方便,對本文計算值進(jìn)行擬合時n取-1.61和-1.80的近似均值-1.70,得到波周期由1.2 s增加到1.5 s、1.8 s和2.1 s時,所對應(yīng)的衰減系數(shù)B值分別從17.34降低到13.39、10.13和8.34,差別也較大?？梢钥闯?波周期對平均垂向速度沿射流中軸衰減的早晚有較大影響，波周期越大,w0/wj衰減越早。

圖9 不同波周期條件下射流中軸上量綱一平均垂向速度分布Fig.9 Dimensionless average vertical velocity distribution along the jet axis under various wave periods

圖10為H=0.10 m、不同波周期條件下射流特征速度半寬沿射流中軸相對距離(z/D)的變化?？梢钥闯?隨著z/D的增加,靜水和波浪條件下的bw均相應(yīng)增加；靜水條件下的bw/D與z/D基本呈現(xiàn)線性關(guān)系。z/D<5時,靜水和波浪條件下的bw基本一致,原因是在這段區(qū)域內(nèi)射流主要受射流初始動量影響,波浪和負(fù)浮力對其影響較小。z/D>5時,隨著波周期增加,bw相應(yīng)增加；波周期T=2.1 s時,z/D=7時bw急劇增加,在z/D=13.33處其增幅開始趨于平緩,這是由于波周期較大的工況下射流軸發(fā)生了明顯偏折。

圖10 不同波周期條件下射流特征速度半寬與量綱一垂向距離關(guān)系Fig.10 Characteristic velocity half-width relationship with the dimensionless vertical distance under various wave periods

3.4 波-射流動量比與射流特性的關(guān)系

Mori和Chang[8]引入波-射流動量比(RM)參數(shù)研究波浪環(huán)境下非浮力射流的振蕩特性:

(13)

式中:A為波幅,m。

本研究中射流流體和環(huán)境水體之間存在著一定密度差,這里引入環(huán)境流體與射流流體的密度比ρa(bǔ)/ρj,對公式(13)進(jìn)行改寫,得到考慮浮力效應(yīng)的波-射流動量比(RMB):

(14)

圖11為各計算工況下衰減系數(shù)(B)與RMB的關(guān)系和擬合直線(其擬合公式形式為B=C1RMB+C2,C1和C2為系數(shù))。從圖11(a)可以看出,波周期一定情況下,隨著波高或波幅增加,RMB增加,B值相應(yīng)降低,兩者呈負(fù)相關(guān)關(guān)系。從圖11(b)可以看出,波高一定情況下,隨著波周期增加,RMB和B值均相應(yīng)減小,兩者呈正相關(guān)關(guān)系?？梢?隨著波高和波周期增加,平均垂向速度沿射流中軸衰減的越早,稀釋和擴(kuò)散效果也越好。

圖11 衰減系數(shù)與波-射流動量比的關(guān)系Fig.11 Relationship between decay coefficient and wave-to-jet momentum ratio

圖12為各組工況下射流平均特征速度半寬(b/D)與RMB的關(guān)系曲線和擬合直線(其擬合公式形式為b/D=D1RMB+D2,D1和D2為系數(shù)),b由公式(15)給出,其中Zmax和Zmin分別為Z方向計算域的最大和最小值。從圖12(a)可以看出,波周期一定情況下,隨著波高或波幅增加,RMB和b/D值均相應(yīng)增加,兩者呈正相關(guān)關(guān)系。從圖12(b)可以看出,波高一定情況下,隨著波周期增加,RMB減小,b/D值相應(yīng)增加,兩者呈負(fù)相關(guān)關(guān)系?？梢?波浪對射流的擴(kuò)展率有較明顯的影響；隨著波高和波周期增加,射流平均特征速度半寬增加,有利于射流鹽水的稀釋和擴(kuò)散。

圖12 平均特征速度半寬與波-射流動量比的關(guān)系Fig.12 Relationship between average characteristic velocity half-width and wave-to-jet momentum ratio

(15)

4 結(jié) 論

本文通過物理模型試驗和數(shù)值模擬相結(jié)合的方法,開展了規(guī)則波作用下垂直圓管負(fù)浮力射流特性研究,探討了波參數(shù)、波射流動量比與射流特性之間的關(guān)系。得到以下結(jié)論:

(1)規(guī)則波作用下垂直圓管負(fù)浮力射流的數(shù)值模擬計算結(jié)果與試驗結(jié)果較吻合,說明建立的數(shù)學(xué)模型是合理可靠的。

(2)波高和周期對平均垂向速度沿射流中軸的衰減指數(shù)影響不大,而對衰減系數(shù)影響較大。隨著波高或波周期的增加,射流特征速度半寬增加,射流中軸上垂向速度的衰減系數(shù)減小。波浪促進(jìn)了射流的振蕩,有助于射流流體與環(huán)境水體的混合,提高了射流流體的擴(kuò)散稀釋效率。

(3)波周期一定、波高變化條件下,射流中軸上量綱一平均垂向速度衰減系數(shù)與考慮浮力效應(yīng)的波-射流動量比呈負(fù)相關(guān)關(guān)系,平均特征速度半寬與考慮浮力效應(yīng)的波-射流動量比呈正相關(guān)關(guān)系。波高一定、波周期變化條件下,射流中軸上量綱一平均垂向速度衰減系數(shù)與考慮浮力效應(yīng)的波-射流動量比呈正相關(guān)關(guān)系,平均特征速度半寬與考慮浮力效應(yīng)的波-射流動量比呈負(fù)相關(guān)關(guān)系。另外,建立了考慮浮力效應(yīng)的波-射流動量比表征的衰減系數(shù)和平均特征速度半寬公式,可對高濃度鹽水排放高度及擴(kuò)展范圍進(jìn)行相關(guān)預(yù)測分析。

需要說明的是,本研究選取的數(shù)值計算尺度都是以室內(nèi)試驗條件為基礎(chǔ)的,入射波也選取相對簡單的規(guī)則波,與實際海況存在較大區(qū)別。下一步需要對實際復(fù)雜海況條件下的負(fù)浮力射流特性進(jìn)行深入研究。