魯 俊,王玲玲,曹小紅

(1.張家港市水利局,江蘇張家港 215600;2.河海大學(xué)水利水電學(xué)院,江蘇南京 210098)

橫流中射流及其溫度標(biāo)量輸運(yùn)大渦模擬

魯 俊1,王玲玲2,曹小紅1

(1.張家港市水利局,江蘇張家港 215600;2.河海大學(xué)水利水電學(xué)院,江蘇南京 210098)

將動(dòng)力擬序渦黏性亞格子應(yīng)力模型拓展到溫度標(biāo)量亞格子模型中,數(shù)值模擬了橫流條件下有、無溫度標(biāo)量場(chǎng)作用的射流,得到的橫流條件下浮力射流的溫度和速度分布與Anwar的試驗(yàn)值吻合一致。在此基礎(chǔ)上,分析了有、無溫度標(biāo)量場(chǎng)作用下射流回流區(qū)域大小和射流軌跡線特性,對(duì)比分析了回流區(qū)域內(nèi)渦心和分離點(diǎn)處湍動(dòng)能和耗散率、擬渦能以及邊界層處擬序結(jié)構(gòu)等湍流特性。計(jì)算結(jié)果表明:溫度場(chǎng)的作用使射流的回流區(qū)域增大,射流速度軌跡線高度增加,回流區(qū)域內(nèi)湍流的湍動(dòng)能增加,邊界處擬序結(jié)構(gòu)的周期性不如無溫度場(chǎng)時(shí)明顯。

射流;橫流;大渦模擬;溫度場(chǎng);擬序結(jié)構(gòu)

在實(shí)際工程中常遇到廢溫水向河流中排放等問題。從水力學(xué)角度來看,這類問題可以概化為橫流中的射流及其標(biāo)量輸運(yùn)問題,因此,對(duì)此進(jìn)行系統(tǒng)研究具有重要的理論意義和實(shí)際應(yīng)用價(jià)值。

許多研究者對(duì)此類問題做了很多試驗(yàn)及理論和數(shù)值模擬研究,如文獻(xiàn)[1-4]研究了射流軌跡線等宏觀分布特性,但未能給出瞬時(shí)流場(chǎng)分布特性;韓會(huì)玲等[5]采用雷諾時(shí)均方程數(shù)值分析了橫流中射流,受雷諾時(shí)均方程的局限,僅能對(duì)射流的軌跡線等宏觀特性進(jìn)行分析,未能給出射流發(fā)展的整個(gè)瞬時(shí)過程。近幾年來,隨著計(jì)算機(jī)運(yùn)算能力的快速提升和大渦模擬理論的成熟,越來越多的研究者開始用大渦模擬技術(shù)來研究這一問題,如Majander等[6]利用經(jīng)典Smagorinsky亞格子模型模擬了橫流中圓孔射流及其標(biāo)量輸運(yùn)問題,發(fā)現(xiàn)采用經(jīng)典Smagorinsky亞格子模型在靠近壁面處誤差較大;Wegner等[7]用大渦模擬技術(shù)分析了橫流中射流摻混圖像特性,但沒有定量分析射流湍流特性。

由于經(jīng)典Smagorinsky亞格子模型模擬射流存在較大誤差以及在壁面處采用壁面函數(shù)或者采用渦黏性衰減函數(shù)來模化渦黏系數(shù)等不足,文獻(xiàn)[8]利用湍流的擬序結(jié)構(gòu),構(gòu)造了新的亞格子應(yīng)力模型,其模型系數(shù)采用動(dòng)力模式來求解,并成功地應(yīng)用于復(fù)雜環(huán)境中射流數(shù)值模擬[9-10]。本文利用此亞格子應(yīng)力模型數(shù)值模擬橫流條件下的射流,并定量分析有無其標(biāo)量輸運(yùn)對(duì)射流湍流特性的影響。

1 數(shù)學(xué)模型的建立

1.1 控制方程

對(duì)不可壓縮的流體流動(dòng)及其溫度輸運(yùn),可用連續(xù)方程、動(dòng)量方程、能量方程和狀態(tài)方程來控制。在大渦模擬理論中,用空間濾波把求解的變量分為能直接求解的大尺度量和需要模化的小尺度量。在Bousineq的假設(shè)下,其在一般坐標(biāo)下的空間濾波后的數(shù)學(xué)控制方程[11]為

式中:ξk、ξl為空間中任意坐標(biāo)(k=1,2,3;l=1,2, 3);J-1為Jacobian矩陣;ui為直角坐標(biāo)下速度分量; p為靜水壓強(qiáng);υ為流體運(yùn)動(dòng)黏性系數(shù);t為時(shí)間; τij為亞格子尺度應(yīng)力;T為溫度;k為分子擴(kuò)散系數(shù),,其中Pr為分子普朗特?cái)?shù);q為標(biāo)量亞格i子輸運(yùn)項(xiàng);ρa(bǔ)為環(huán)境流體密度;ρ為射流密度;gi為重力加速度;xi為直角坐標(biāo)。

亞格子尺度應(yīng)力τij采用由文獻(xiàn)[8]提出的考慮擬序結(jié)構(gòu)作用的動(dòng)力擬序渦黏性亞格子模型?；?類似可把上述模型拓展到標(biāo)量亞格子輸運(yùn)項(xiàng)qi模化,其結(jié)果如下:

式中:Cs為Smagorinsky常數(shù);Δ為濾波尺度,定義為Δ=1/2;νt為湍動(dòng)擴(kuò)散系數(shù);為反對(duì)稱變形張量;為對(duì)稱變形張量;α為加權(quán)權(quán)重,取α= 0.5;Prt為湍動(dòng)普朗特?cái)?shù),同Smagorinsky常數(shù)Cs一樣,其值一般不是一個(gè)常數(shù),是個(gè)動(dòng)態(tài)變化的值,本文利用動(dòng)力模式的思想來確定。

用式(6)確定系數(shù)Cs在空間分布變異較大,從而導(dǎo)致νt在空間變化劇烈,這種劇烈的變化往往會(huì)導(dǎo)致數(shù)值計(jì)算的不穩(wěn)定。為此采用文獻(xiàn)[8]提出的加權(quán)時(shí)間平均法:

同理,湍動(dòng)普朗特?cái)?shù)采用加權(quán)時(shí)間平均法如下:

式中:λ為時(shí)間加權(quán)值,數(shù)值計(jì)算表明λ取值為0.65最佳;下標(biāo)(n-1)表示上一時(shí)刻。

1.2 數(shù)值求解方法

為了計(jì)算不規(guī)則的區(qū)域和自由表面,采用垂直σ坐標(biāo),把一般坐標(biāo)下的方程采用以下方程轉(zhuǎn)化為σ坐標(biāo):

式中:η為自由表面波高;h為靜止水深。當(dāng)z在-h ～η之間變化時(shí),σ在0～1之間變化。把式(9)進(jìn)行微分求導(dǎo)代入式(1)、式(2)和式(3),得到σ坐標(biāo)的表達(dá)式,詳細(xì)推導(dǎo)過程見文獻(xiàn)[8,10],數(shù)值方法采用分步法,動(dòng)量方程分為3個(gè)計(jì)算過程,即對(duì)流計(jì)算過程、擴(kuò)散計(jì)算過程和壓力傳播過程;能量方程分為兩個(gè)計(jì)算過程,即對(duì)流計(jì)算過程和擴(kuò)散計(jì)算過程。壓力的變化通過求解Poisson壓力方程來滿足連續(xù)方程得到,采用計(jì)算速度快、魯棒性好的CGSTAB方法求解。程序采用Fortran語言編寫。

1.3 邊界條件

a.進(jìn)口邊界條件分為橫流進(jìn)口條件和射流進(jìn)口條件。為了快速產(chǎn)生湍流,橫流和射流進(jìn)口速度場(chǎng)采用給定的平均流速再加由方位角產(chǎn)生的隨機(jī)擾動(dòng)給定脈動(dòng)速度;射流進(jìn)口溫度場(chǎng)采用給定的平均溫度再加由方位角產(chǎn)生的隨機(jī)擾動(dòng)給定脈動(dòng)溫度。

c.壁面條件。底部固壁采用不可滑移邊界條件;側(cè)向固壁采用可滑移邊界條件。

d.自由表面。采用Lagrange-Euler法求解以下標(biāo)高函數(shù)捕捉自由表面形狀:

2 數(shù)值模擬結(jié)果及分析

2.1 數(shù)學(xué)模型驗(yàn)證

應(yīng)用上述模型和計(jì)算方法對(duì)Anwar[13]試驗(yàn)成果進(jìn)行對(duì)比驗(yàn)證。試驗(yàn)水槽長(zhǎng)10 m,寬0.6 m,高1.2m。76℃的熱水從射流窄縫處垂直射向溫度為12℃的環(huán)境流體。射流孔口寬為0.01m,長(zhǎng)為0.6 m。環(huán)境流體采用流量為0.23m3/s泵循環(huán)。選取其射流比R=W0/Ua=2的浮力射流工況(簡(jiǎn)稱CR1工況)進(jìn)行驗(yàn)證,其中W0為射流的垂向速度,Ua為橫流速度。

計(jì)算區(qū)域長(zhǎng)取為6 m,為了減少寬度方向網(wǎng)格數(shù)量和采用可滑移邊界條件,計(jì)算寬度僅取為0.1 m,高為1.0 m。經(jīng)網(wǎng)格無關(guān)性要求驗(yàn)證后,在x(槽長(zhǎng))方向、y(槽寬)方向和z(高度)方向(即σ方向)分別采用309、11、85網(wǎng)格劃分,其中,在x方向射流口寬度d范圍內(nèi)均勻布置11個(gè)網(wǎng)格節(jié)點(diǎn),非均勻網(wǎng)格中最小網(wǎng)格尺度為0.002,最大網(wǎng)格尺度為0.083;y方向采用均勻網(wǎng)格;在z(垂線)方向采用不均勻網(wǎng)格劃分,最小網(wǎng)格尺度為0.003,最大網(wǎng)格尺度為0.032。時(shí)間步長(zhǎng)Δt=2.5×10-4s,滿足計(jì)算過程對(duì)流穩(wěn)定性條件和擴(kuò)散穩(wěn)定性條件,計(jì)算時(shí)間為40 s。計(jì)算在CPU主頻為3.0 G、內(nèi)存為2 G的計(jì)算機(jī)上計(jì)算約耗時(shí)300 h。

圖1為時(shí)間平均下的溫度計(jì)算值和試驗(yàn)值的對(duì)比。為了和文獻(xiàn)中[13]中坐標(biāo)取值一致,圖中相對(duì)長(zhǎng)度δ按文獻(xiàn)[13]方式定義為垂直斷面溫度T-Ta

圖1 平均溫度大小分布

圖2 平均速度大小分布

2.2 計(jì)算結(jié)果分析

為了進(jìn)一步對(duì)比分析溫度場(chǎng)對(duì)射流特性的影響,計(jì)算同樣工況下無溫度場(chǎng)的湍動(dòng)射流(簡(jiǎn)稱CR2工況)。

2.2.1 速度場(chǎng)、雷諾應(yīng)力場(chǎng)和溫度場(chǎng)分析

圖3 CR1工況垂直平面流線

圖4 CR1工況垂直平面x方向速度等值線(單位:m/s)

圖3 為CR1工況時(shí)間平均下浮力射流流線。從圖3可知,環(huán)境流體受到窄縫射流的阻滯而形成繞流,因此在射流的上邊緣和下邊緣形成不對(duì)稱的壓強(qiáng)分布,此不平衡力使射流發(fā)生彎曲,彎曲程度受到來流速度和射流速度及其溫度差產(chǎn)生浮力的綜合影響。另外,圖3還顯示在射流前后下端區(qū)域呈現(xiàn)出三角形狀角渦,同時(shí)在射流背水面呈現(xiàn)出大范圍的回流區(qū),x方向流速為負(fù)值,如圖4所示。同時(shí)在橫流主流區(qū)和回流區(qū)之間形成一個(gè)強(qiáng)的切應(yīng)力區(qū),如圖5所示。圖6為時(shí)間平均下溫度等值線,可以看出,在射流出口附近,由于射流的初始動(dòng)量較大,最大溫度在垂向上逐漸抬升,直至某一距離后,由于射流的動(dòng)量作用減小,最大溫度值受到橫流和回流區(qū)的負(fù)速度影響,逐漸降落,直至貼壁,并入侵回流區(qū)。這一現(xiàn)象也為韓會(huì)玲等[5]用RANS雙方程數(shù)值模擬得到。結(jié)合Anwar[13]試驗(yàn)成果發(fā)現(xiàn),隨著射流比的減小,射流背水面的回流區(qū)域減小,流線的彎曲程度越大,流動(dòng)越貼近壁面;同時(shí),橫流的對(duì)流速度作用增大,在射流出口處射流溫度的衰減也越快,環(huán)境流體與射流流體的溫度摻混稀釋能力也越強(qiáng)。對(duì)比圖3和圖7,可以明顯地得出溫度場(chǎng)中正浮力作用使回流區(qū)域增大的結(jié)論。

圖5 CR1工況垂直平面雷諾應(yīng)力等值線(單位:Pa)

圖6 CR1工況垂直平面溫度等值線(單位:℃)

圖7 CR2工況垂直平面流線

2.2.2 射流軌跡線、湍動(dòng)能、耗散率和擬渦能分析

為了進(jìn)一步研究溫度標(biāo)量場(chǎng)對(duì)射流軌跡線的影響,圖8給出了射流邊緣速度內(nèi)外軌跡線。由圖8可知,約在z/d=0～7的范圍內(nèi),浮力射流和湍動(dòng)射流的速度內(nèi)外軌跡線重合,說明在射流出口附近動(dòng)量作用遠(yuǎn)遠(yuǎn)大于浮力作用;越遠(yuǎn)離射流,浮力作用越明顯,射流的速度內(nèi)外軌跡線都向上移動(dòng),同時(shí)速度外軌跡線附近溫差大而形成的浮力作用比速度內(nèi)軌跡線附近形成的浮力作用強(qiáng),致使在溫度場(chǎng)浮力作用下,CR1工況的速度內(nèi)外軌跡線之間的距離比CR2工況的速度內(nèi)外軌跡線之間的距離大。圖9為CR1工況最大速度軌跡線、最大溫度軌跡線和CR2工況最大速度軌跡線,可以看出在溫度場(chǎng)的正浮力作用下最大速度軌跡線高于無溫度場(chǎng)浮力作用下最大速度軌跡線,最大速度軌跡線與最大溫度軌跡線不重合,最大溫度軌跡線低于最大速度軌跡線。值得一提的是約在z/d=0～7的范圍內(nèi),浮力射流和湍動(dòng)射流的最大速度軌跡線重合,這從另外一個(gè)角度說明在射流出口附近動(dòng)量作用遠(yuǎn)大于浮力作用。

圖8 射流速度內(nèi)外軌跡線

圖9 射流最大速度和最大溫度軌跡線

圖10 渦心和分離點(diǎn)處垂直斷面湍動(dòng)能

圖10 和圖11分別為CR1工況和CR2工況渦心和分離點(diǎn)處垂直斷面的湍動(dòng)能和耗散率。從圖10可知,無論是否受溫度影響,回流區(qū)域以上的湍動(dòng)能k沿z軸正向逐漸減小,符合湍流湍動(dòng)能逐漸衰減的物理機(jī)制。但在渦心處不受溫度影響的湍動(dòng)能明顯大于受溫度影響的湍動(dòng)能,在分離點(diǎn)處兩者卻相反。而渦心和分離點(diǎn)處不受溫度影響的耗散率ε均大于受溫度影響的耗散率(圖11),這說明符合溫度標(biāo)量場(chǎng)在回流區(qū)域內(nèi)可以延緩?fù)牧髂芰克p這一物理機(jī)制,這是由于溫度場(chǎng)的脈動(dòng)可以增加其流體湍動(dòng)能,而這一點(diǎn)也為渦心和分離點(diǎn)處垂直斷面擬渦能0.5Ω所證實(shí),如圖12所示。

圖11 渦心和分離點(diǎn)處垂直斷面耗散率

圖12 渦心和分離點(diǎn)處垂直斷面擬渦能

2.2.3 擬序結(jié)構(gòu)及能譜分析

數(shù)值模擬的瞬時(shí)過程表明在回流區(qū)的分離點(diǎn)處,有渦旋從回流區(qū)脫離并向下游發(fā)展,如圖13所示。從圖13可知,近壁有大的渦團(tuán)存在并逐漸向下游發(fā)展,最后脫離整個(gè)計(jì)算區(qū)域。同時(shí)圖13還顯示瞬時(shí)溫度標(biāo)量渦團(tuán)和速度渦旋存在一定協(xié)同性。為了進(jìn)一步研究其變化特性,在分離點(diǎn)的后面建立一個(gè)典型觀察點(diǎn)A,并研究其速度分量時(shí)間序列。點(diǎn)A位于x=1.32 m、y=0.05 m和z=0.04 m處,即z+=zuτ/ν=85,其中為切應(yīng)力。

圖13 CR1工況瞬時(shí)速度場(chǎng)和溫度等值線(單位:℃)

圖14 和圖15分別為CR1工況和CR2工況下點(diǎn)A處的速度分量時(shí)間序列和相應(yīng)速度分量能譜圖。低頻處能譜可以理解為平均速度的能譜,而高頻處能譜可以理解為脈動(dòng)速度的能譜。對(duì)比圖15中兩種不同工況下速度分量能譜,可以發(fā)現(xiàn)約在0.2 Hz以下范圍內(nèi),CR2工況中的平均速度能譜值要低于CR1工況下的平均速度能譜值,這說明有溫度梯度作用下平均速度所做的功要比無溫度梯度作用下平均速度所做的功要大。而約在1 Hz以上范圍內(nèi),CR2工況下的脈動(dòng)速度能譜值要低于CR1工況下的脈動(dòng)速度能譜值,這說明有溫度梯度下脈動(dòng)速度的能量耗散要比無溫度梯度下脈動(dòng)速度的能量耗散要小,從而可說明溫度主動(dòng)標(biāo)量場(chǎng)的作用可以增加湍流的湍動(dòng)能并減小湍流的湍流能耗散。另外,從圖15還可發(fā)現(xiàn)約在0.2 Hz以上,CR2工況下擬序結(jié)構(gòu)的能譜值迅速增大,這說明擬序結(jié)構(gòu)的作用在無溫度梯度作用下增強(qiáng),計(jì)算結(jié)果顯示在其工況下擬序結(jié)構(gòu)頻率范圍約在0.2～1 Hz。

圖15 點(diǎn)A處速度分量能譜

3 結(jié) 語

本文將動(dòng)力擬序結(jié)構(gòu)渦黏性亞格子應(yīng)力模型拓展到溫度標(biāo)量亞格子模型中,數(shù)值模擬了橫流條件下有、無溫度標(biāo)量場(chǎng)作用的射流。數(shù)值模擬結(jié)果表明模型在計(jì)算橫流條件下浮力射流得出的溫度和速度分布值與試驗(yàn)結(jié)果一致。在主動(dòng)溫度標(biāo)量場(chǎng)作用下,射流的回流區(qū)域增大,射流速度軌跡線高度增加。溫度場(chǎng)作用使回流區(qū)域內(nèi)湍流的湍動(dòng)能增加,使湍流的湍動(dòng)能耗散減小,從而使湍流總能量衰減減慢。速度分量能譜圖還顯示了溫度場(chǎng)的正浮力作用使擬序結(jié)構(gòu)的周期性不如無溫度場(chǎng)時(shí)明顯,同時(shí)溫度梯度作用使湍流的總動(dòng)能量增加和能量耗散減小。

[1]ANDREOPOULOS J,PRATURI A,RODI W.Experiments on vertical plane buoyant jets in shallow water[J].J Fluid Mech,1986,168:305-336.

[2]PLESNIAK M W,CUSANO D M.Scalar mixing in a confined rectangular jet in cross-flow[J].J Fluid Mech, 2005,524:1-45.

[3]PETERSON S D,PLESNIAK M W.Evolution of jets emanating from short holes into cross-flow[J].J Fluid Mech,2004,503:57-91.

[4]CHEN K S,HWANG J Y.Experimental study on the mixing of one and dual-line heated jets with a cold crossflow in a confined channel[J].AIAA J,1991 29:353-360.

[5]韓會(huì)玲,李煒.橫流中正浮力射流近區(qū)特性預(yù)報(bào)[J].河北農(nóng)業(yè)大學(xué)學(xué)報(bào),1997:20(3):112-119.(HAN Huiling,LI Wei.Prediction of characteristics for near field of positive buoyant jets in cross-flow[J].Journal of Agricultural University of Hebei,1997,20(3):112-119. (in Chinese))

[6]MAJANDER P,SIIKONEN T.Large-eddy simulation of a round jet in a cross-flow[J].Journal of Heat and Fluid Flow,2006,27:402-415

[7]WEGNER B,HUAI Y A S.Comparative study of turbulent mixing in jet in cross-flow configurations using LES[J]. International Journal of Heat and Fluid Flow,2004,25: 767-775.

[8]LU Jun,TANG hongwu,WANG Lingling.A novel dynamic coherent eddy model and application to LES of turbulent jet with free surface[J].Science in China:Series G, 2010,53(9):1671-1680.

[9]LU Jun,WANG Lingling,TANG Hongwu,et al.Large eddy simulation of vertical turbulent jets under JONSWAP waves[J].Acta Mechanica Sinica,2011,27(2):189-199.

[10]LU Jun,WANG Lingling,TANG Hongwu,et al.Numerical investigating of vertical turbulent jets in different types of waves[J].China Ocean Engineering,2010,24(4):611-626.

[11]SAUGAUT P.Large eddy simulation for incompressible flows[M].3rd ed.Singpore:Springer press,2006.

[12]LILLY D K.A proposed modification of the Germano subgrid-scale closure method[J].Physics of Fluids,1992, 4(3):633-635.

[13]ANWAR H O.Two-dimensional buoyant jet in a current [J].Journal of Engineering Mathematics,1973,7(4): 297-311.

Large eddy simulation of jets with and without temperature scalar in a current

//LU Jun1,WANG Lingling2,CAO Xiaohong1(1.Water Conservancy Bureau of Zhangjiagang,Zhangjiagang 215600,China;2.College of Water Conservancy and Hydropower,Hohai University,Nanjing 210098,China)

The jet in cross flow was simulated using a dynamic coherence eddy subgrid stress model,which was expanded into a subgrid heat flux model,under the action of temperature field.The results show that the temperature and velocity distributions of a buoyant jet in cross flow obtained from the numerical simulation are in good agreement with Anwar’s experimental data.The characteristics of backflow zones and the trajectory lines of the jets were investigated with and without temperature field.The turbulence kinetic energy and the turbulence eddy dissipation in the vortex core were compared with those of the separation point of backflow zones.The calculated results show that the areas of the backflow zones,the jet trajectory line height,and the turbulence kinetic energy increase under the action of temperature field.The period of the coherent structures near the boundary layer becomes unobvious under the action of temperature field.

jets;cross flow;large eddy simulation;scalar turbulence;coherent structures

10.3880/j.issn.10067647.2013.02.006

O358

A

10067647(2013)02002606

2012-05-15 編輯:熊水斌)

國(guó)家自然科學(xué)基金(51279050)

魯俊(1980—),男,安徽馬鞍山人,工程師,博士,主要從事湍流與環(huán)境水力學(xué)數(shù)值模擬研究。E-mail:55939818@qq.com