橫流中射流及其溫度標量輸運大渦模擬

魯 俊,王玲玲,曹小紅

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

橫流中射流及其溫度標量輸運大渦模擬

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

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

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

射流;橫流;大渦模擬;溫度場;擬序結構

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

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

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

1 數學模型的建立

1.1 控制方程

對不可壓縮的流體流動及其溫度輸運,可用連續(xù)方程、動量方程、能量方程和狀態(tài)方程來控制。在大渦模擬理論中,用空間濾波把求解的變量分為能直接求解的大尺度量和需要?；男〕叨攘?。在Bousineq的假設下,其在一般坐標下的空間濾波后的數學控制方程[11]為

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

亞格子尺度應力τij采用由文獻[8]提出的考慮擬序結構作用的動力擬序渦黏性亞格子模型模化,類似可把上述模型拓展到標量亞格子輸運項qi模化,其結果如下:

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

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

同理,湍動普朗特數采用加權時間平均法如下:

式中:λ為時間加權值,數值計算表明λ取值為0.65最佳;下標(n-1)表示上一時刻。

1.2 數值求解方法

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

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

1.3 邊界條件

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

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

d.自由表面。采用Lagrange-Euler法求解以下標高函數捕捉自由表面形狀:

2 數值模擬結果及分析

2.1 數學模型驗證

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

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

圖1為時間平均下的溫度計算值和試驗值的對比。為了和文獻中[13]中坐標取值一致,圖中相對長度δ按文獻[13]方式定義為垂直斷面溫度T-Ta

圖1 平均溫度大小分布

圖2 平均速度大小分布

2.2 計算結果分析

為了進一步對比分析溫度場對射流特性的影響,計算同樣工況下無溫度場的湍動射流(簡稱CR2工況)。

2.2.1 速度場、雷諾應力場和溫度場分析

圖3 CR1工況垂直平面流線

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

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

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

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

圖7 CR2工況垂直平面流線

2.2.2 射流軌跡線、湍動能、耗散率和擬渦能分析

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

圖8 射流速度內外軌跡線

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

圖10 渦心和分離點處垂直斷面湍動能

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

圖11 渦心和分離點處垂直斷面耗散率

圖12 渦心和分離點處垂直斷面擬渦能

2.2.3 擬序結構及能譜分析

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

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

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

圖15 點A處速度分量能譜

3 結 語

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

[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]韓會玲,李煒.橫流中正浮力射流近區(qū)特性預報[J].河北農業(yè)大學學報,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 編輯:熊水斌)

國家自然科學基金(51279050)

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