U形彎道中湍流對流換熱的數值模擬

徐文潔,　賴煥新

(華東理工大學承壓系統與安全教育部重點實驗室,上海 200237)

U形彎道中湍流對流換熱的數值模擬

徐文潔,賴煥新

(華東理工大學承壓系統與安全教育部重點實驗室,上海 200237)

摘要：網格獨立解是湍流對流換熱數值研究的前提,而邊界條件的恰當處理則是模擬結果準確性的關鍵。結合方形截面U形彎道中的湍流對流換熱,采用SST k-ω模型對彎道內流場和溫度場進行計算,重點分析了速度場和溫度場的網格獨立解。同時結合第二類邊界條件(熱流密度條件),對比分析其兩種處理方式對數值解準確性的影響。結果表明,由于不可壓縮流動中速度場和溫度場之間較弱的耦合關系,獲得網格獨立解所需的網格數目存在差別,溫度場對網格密度要求更高。熱流密度邊界條件的兩種處理方式的比較表明,Kader(1981)的代數方法在網格粗糙的情況下仍然能夠較好地吻合實驗值,具有更好的通用性和網格適應性。

關鍵詞：網格獨立解; 湍流; 對流換熱; 邊界條件

U形彎道中的對流傳熱現象廣泛存在于熱能動力工程設備中,如換熱器盤管內的流動,其中流動與換熱效率對熱能動力工程設備的效率、阻力性能等有很大的影響。為此,前人針對U形彎道內的對流換熱做了很多工作。Chang等[1]采用激光多普勒儀測量了U形彎道中發展流動的平均速度和雷諾數。Johnson和Launder[2]在Chang等的基礎上,對U形彎道內的傳熱系數和溫度場分布進行了詳細的研究。Chang[3]用實驗研究了U形彎道中有旋流的傳熱特征。數值模擬方面,Iacovides等[4]研究發現二階矩差分模型(DSM)的預測結果比常用的二階矩代數模型(ASM)更接近實驗值,并提出了一種引入無壁面反射項的新DSM模型,其能顯著改善計算結果。Suga[5-6]改進了廣義梯度擴散假設模型(GGDH),使其對傳熱的計算結果更為可靠。Etemad等[7-9]對直管和U形彎道進行了數值模擬,并選用了不同的模型進行比較分析,發現U形彎道入口段長度和上游邊界層厚度對流場和溫度場的分布極為重要。Shin等[10]提出了一種橢圓概念湍流熱流密度模型,并對U形彎道內的湍流流動和傳熱進行了數值模擬。Nobari等[11]研究了不同旋轉角度對U形彎道內傳熱效率的影響。Moon等[12]以加強對流傳熱和減少摩擦損失為目標,研究了不同結構的旋轉U形彎道內的流動和傳熱。Salameh等[13]選用了不同的湍流模型并結合不同的壁面邊界條件數值模擬了U形彎道內的湍流流動及傳熱,結果發現結合了標準壁面函數法的重整化的k-ε模型和可實現的k-ε模型所預測得到的結果與實驗值最為吻合。Shen等[14]研究了旋轉對帶肋條U形彎道內對流傳熱的影響。由于U形彎道內流動和傳熱現象的復雜性,數值模擬中湍流模型的通用性、算法和網格的精度、邊界條件的恰當處理等,都將對其流場與溫度場的準確預測造成重要的影響。

網格獨立解是湍流對流換熱數值模擬的前提和基礎。對于湍流強制對流換熱,當流體的密度和黏度不受溫度影響時,速度場對溫度分布起決定性作用,但溫度場對速度場分布的反作用則較小甚至可以忽略,兩者為弱耦合關系。此時,速度場和溫度場不一定能同時達到網格獨立解。前人的上述關于U形彎道湍流對流換熱的數值研究中,并未見到同時對速度場和溫度場進行網格無關性的檢查。另一方面,邊界條件的處理也是計算結果準確性的關鍵。當能量方程為第二或第三類邊界條件時,需要對未知的邊界節點進行離散,其中需要采用一定的插值法對邊界點溫度外推計算,而此時這類溫度邊界外推的方法對溫度全場結果的準確性具有十分重要的影響。

基于此,本文結合具有詳細測量數據的方形截面180°U形彎道[9]中的湍流對流換熱,采用SSTk-ω模型對彎道內流場和溫度場進行計算,重點分析速度場和溫度場的網格獨立解;同時,結合第二類溫度邊界條件(熱流密度條件),對比分析兩種壁面溫度處理法對傳熱結果準確性的影響。

彎管幾何結構如圖1(a)所示，具體尺寸如圖1(b)所示。橫截面邊長D=88.9 mm,半徑比Rc/D=3.357;入口的平均速度Ub=8.484 2 m/s;2z/D=0表示彎道對稱中心面;2z/D=1表示彎道側壁面。其中z方向為左壁指向右壁,z=0表示左壁,z=D表示右壁。定義量綱為一徑向坐標R*=(R-Ri)/(Ro-Ri),其中R*=0表示彎道內壁,R*=1表示彎道外壁,Ro表示外徑，Ri表示內徑。

圖1　正方形截面180°U形管結構及坐標系示意圖

1流動模型與計算方法

1.1湍流計算模型

本文計算采用Menter在2003年修正版本的SSTk-ω模型[15],其湍流動能k和湍動頻率ω方程分別為

(1)

(2)

式中:ui和xi分別代表直角坐標系中各方向的速度和坐標;ν和νt分別表示流體黏度和湍流黏度;Pk表示湍流動能生成率。閾值函數

(3)

(4)

閾值函數F2定義為

(5)

對湍流動能生成率進行修正以防止在流動停滯區造成湍流堆積:

(6)

k和ω方程中所涉及到的常數如α、β、σk和σω,均按照φ=φ1F1+φ2(1-F1) 進行計算。其中,φ1和φ2分別表示k-ω模型和標準k-ε模型中相應的常數。β*=0.09,α1=5/9,α2=0.44,β1=3/40,β2=0.082 8,σk1=0.85,σk2=1.0,σω1=0.5,σω2=0.856,σt=1.0。

計算時采用課題組自有的高性能代碼,利用有限體積離散格式和非交錯網格SIMPLE算法求解三維不可壓雷諾時均Navier-Stokes方程。計算采用的離散格式為:對流項采用三階SMART格式,壓力和擴散項采用二階中心差分格式。利用經典的Rhie-Chow動量插值方法[16]來解決速度-壓力之間的耦合問題。代數方程的求解采用了交替方向迭代法(Alternation Direction Iteration,ADI)[17]。本文使用的數值方法和代碼,在前期發表的論文中已有詳細的介紹,請參考文獻[18-20]。

1.2邊界條件的處理

根據文獻[2]和文獻[9]中的流體條件,流動介質為空氣,忽略其浮升力。進口流體按照Re=56 800和模型尺寸計算得到,平均速度Ub=8.484 2 m/s。出口給定充分發展的邊界條件,即假設流動參數沿流向的一階偏導數為零。

計算流體的初始溫度為273 K,4個固體壁面均設為恒定熱流密度邊界條件,熱流密度qw=250 W/m2。入口面的湍動頻率ωin=10Ub/D。壁面采用無滑移條件,速度設為0。壁面處k的邊界條件取為k=0;壁面處的湍動速率ωw=10×6ν/β1y2。

當壁面熱流密度恒定時,需要對第1層內節點進行離散,再選用一定的插值法計算得到壁面溫度。本文采用了兩種方法,第1種方法為近壁層流假設方法(Laminar Hypothesis,LH)。實際上在非常靠近壁面的地區,湍流脈動動能強烈地被衰減,而耗散率則達到其最大值,因此分子黏性的作用會變得顯著起來[17]。此時,能量方程在邊界上的導熱系數就是流體的分子導熱系數λ,其壁面溫度為

(7)

其中:TP表示第1層內節點的流體溫度;yP表示第1層內節點到壁面的距離。第2種為Kader[21]提出的代數方法,其中壁面量綱為一的溫度表達式為

(8)

其中y+表示第1層內節點的量綱為一的距離。

(9)

(10)

將量綱為一溫度定義為

(11)

其中特征溫度Tτ由熱流密度qw、密度ρ、比熱容cp和剪切速率uτ計算得到;Tw和TP分別表示壁面溫度和第1層內節點溫度,剪切速率為

(12)

因此壁面溫度:

(13)

2計算結果與分析

2.1網格無關性檢查

本文首先對速度場和溫度場進行了網格無關性檢查,網格方案如表1所示。不同網格方案下的主流速度分布及溫度分布見圖2與圖3。

表1　網格節點布置方案

從圖2和圖3中可以看出,方案1與方案2的計算結果相差較大。進一步加密網格后,對于速度場,方案3與方案2的計算結果已趨于一致,只在部分區域有微小差異。因此,可以認為方案2的網格密度可為SSTk-ω模型求解速度場提供網格獨立解;然而從圖3中可以觀察到,方案2與方案3的溫度計算值差異較大,即方案2的網格密度并不能為求解溫度場提供網格獨立解。直到網格進一步加密,方案3與方案4的計算結果趨于一致,并與實驗值基本吻合。因此認為方案3的網格密度才能為求解溫度場提供網絡獨立解。對比速度場和溫度場的網格無關性檢查可以發現,兩者得到網格獨立解所需的網格數目不同,溫度場的計算對網格密度要求更高。因此,在研究對象包含對流換熱問題時,單一對速度場進行網格無關性檢查是不充分的。

圖2　不同網格方案下的主流速度分布圖

圖3　不同網格方案下的溫度分布圖

2.2主流速度

圖4給出了沿流向發展各截面的主流速度計算值與實驗值的分布圖。從圖中可以看到,隨著流體進入彎道,二次流導致流體速度峰值逐漸向外壁偏移。在彎曲段中心處(θ=90°),速度曲線中明顯新出現了一個“凹陷”現象,這個低速區域的形成,是對稱面附近的二次流向上、下壁面劇烈運動導致的。從圖4可以看到,計算結果與實驗結果吻合良好,隨著流體向彎道后半段流動,這個“凹陷”現象逐漸減弱,在彎曲段θ=130°,“凹陷”退化成一個小波紋。總體來看,盡管U形管內流體流動情況復雜,但計算結果始終與實驗結果良好地吻合。

2.3二次流速度

圖5示出了各截面內的二次流速度計算值與實驗值的比較。在θ=45°處,除了在靠近側壁附近區域(2z/D=0.75),出現了反向二次流,其他位置均為正向二次流。說明在θ=45°處,二次流沿著對稱中心面由內壁流向外壁,而沿著側壁由外壁向內壁返回。在θ=90°截面處,靠近側壁附近區域(2z/D=0.75),實驗值顯示,所有位置處的二次流均由內壁流向外壁,二次返回流僅存在于小于D/8寬度的窄小范圍內,而計算結果在D/8位置處捕捉到二次返回流。由于在該位置處計算結果與實驗結果有較大偏差,為了驗證計算的可靠性,本文用商業軟件Fluent對流場進行模擬,采用同樣的湍流模型和網格數,結果如圖5所示。可以看出,Fluent的結果與本文計算趨勢一致,同樣在D/8位置處出現了二次返回流。與實驗結果的比較表明,本文計算結果整體上優于Fluent的模擬結果。

圖4　各截面的主流速度分布圖

2.4溫度場

本文采用了兩種方法對壁面溫度進行處理,其計算結果分別如圖6～圖9所示。為了與網格分布方案3的計算結果形成對比,并比較兩種壁面溫度處理法在較為粗糙網格條件下的差別,首先在網格分布方案2(60×60×170)的條件下,對比兩種壁面溫度處理法對溫度場計算結果準確性的影響;并且在網格較為粗糙時,改變近壁網格密度(即改變y+),其結果變化也將更為明顯,有利于更好地分析討論。從圖6～圖9中可以看到,在網格分布方案2下,近壁層流假設法和Kader代數方法計算得到的溫度分布曲線完全重合,即當第1層內節點位于近壁黏性底層時,兩種處理法沒有差異,并且在近壁區域,計算結果與實驗結果基本吻合。但在遠離壁面區域,計算結果與實驗結果有較大的偏差。為了進一步研究兩種處理法對溫度場計算結果的預測能力,在保證網格總數不變的前提下,改變網格分布,使y+≈12(原y+≈0.7)。此時兩種處理法有了明顯不同的表現:近壁層流假設法計算得到的溫度曲線與實驗結果偏離更大,而Kader代數法更貼近了實驗結果,在對稱面附近區域這種差異尤為明顯。

圖5　各截面的二次流速度分布圖

再選用網格分布方案3(80×80×170,y+≈0.7)對溫度場進行數值計算。同樣地,兩種處理法的計算曲線完全重合,并且與網格分布方案2對比可以發現,無論在近壁區域還是對稱面區域,其計算結果均與實驗值吻合良好,再次說明了溫度場的計算對網格密度要求較高,速度場達到網格獨立解并不是溫度場計算的充分條件。

2.5Nusselt數

Nusselt數是流體力學中的一個量綱為一數,表示對流換熱的強烈程度,是傳熱研究領域的重要衡量指標。其定義為

(14)

其中:qw代表熱流密度;DH表示彎道橫截面水力直徑;TW表示壁面溫度;TB代表按體積計算的流體平均溫度;λ代表計算流體的導熱系數。

圖10給出了兩種壁面溫度處理法在不同網格分布方案的情況下得到的Nu數與實驗值的比較。從圖中可以看到,與溫度場相同,在y+≈0.7時,兩種壁面溫度處理法的計算結果完全重合;但與溫度場結果分析不同的是,當y+≈12時,兩種壁面溫度處理法的計算結果均與實驗值有較大差異。說明在研究Nu數時,要求壁面網格需要達到一定的密度。另一方面,在同樣y+≈0.7的情況下,比較兩種網格分布方案(80×80和60×60)可以發現,幾乎在任意截面、任意位置處,80×80網格分布方案下其計算結果都比60×60網格分布方案下的計算結果更貼近實驗值,再次說明了傳熱問題的研究對網格密度的要求更高。

圖6　各截面的溫度分布圖(θ=45°)

圖7　各截面的溫度分布圖(θ=90°)

圖8　各截面的溫度分布圖(θ=135°)

圖9　各截面的溫度分布圖(θ=180°)

圖10　彎曲段各截面Nu數分布圖

3結論

本文針對熱能動力工程中廣泛出現的U型彎道對流換熱問題,結合具有經典實驗數據的實例,研究了速度場和溫度場網格無關解所需網格數不同的問題,同時比較了熱流邊界條件下壁面溫度外推的兩種方式,主要結論如下:

(1)湍流對流換熱計算中,溫度場的計算對網格密度要求更高,單一地只檢查速度場的網格無關性并不充分。

(2)彎道內部溫度場分布的計算的準確性更多地取決于整體網格密度,而不僅僅是壁面處網格密度;當計算資源受到限制,網格較為粗糙時,與近壁層流假設方法相比,Kader代數方法與實驗值吻合明顯更好,具有更好的通用性和網格適應性。

(3)雖然Kader代數方法在網格較粗糙時可獲得整體上與實驗吻合的溫度場分布,但由于壁面處溫度梯度較大,此時Nu數的計算仍應盡可能采用y+值小的網格。

參考文獻：

[1]CHANG S M,HUMPHREY J A C,MODAVI A.Turbulent-flow in a strongly curved u-bend and downstream tangent of square cross-sections[J].PhysicoChemical Hydrodynamics,1983,4(3):243-269.

[2]JOHNSON R W,LAUNDER B E.Local Nusselt number and temperature field in turbulent flow through a heated square-sectioned U-bend[J].International Journal of Heat and Fluid Flow,1985,6(3):171-180.

[3]CHANG T.An investigation of heat transfer characteristics of swirling flow in a 180° circular section bend with uniform heat flux[J].KSME International Journal,2003,17(10):1520-1532.

[4]IACOVIDES H,LAUNDER B E,LI H.Application of a reflection-free DSM to turbulent flow and heat transfer in a square-sectioned U-bend[J].Experimental Thermal & Fluid Science,1996,13(96):419-429.

[5]SUGA K.Predicting turbulence and heat transfer in 3-D curved ducts by near-wall second moment closures[J].International Journal of Heat & Mass Transfer,2003,46(1):161-173.

[6]SUGA K,HORINOUCHI N,NAGAOKA M.Application of a higher order GGDH heat flux model to three-dimensional turbulent U-bend duct heat transfer[J].Transactions of the Asme Serie C Journal of Heat Transfer,2003,125(1):200-203.

[7]ETEMAD S,SUNDéN B.Numerical analysis of turbulent convective heat transfer processes in a square-sectioned U-bend duct[C]//15th Australasian Fluid Mechanics Conference,Sydney,Australia:The University of Sydney,2004:13-17.

[8]ETEMAD S,SUNDéN B.Analysis of developing and fully developed turbulent flow and heat transfer in a square-sectioned U-bend duct[C]//Proceedings of the ASME Summer Heat Transfer Conference.San Francisco,CA:[s.n.],2005:15-22.

[9]ETEMAD S,SUNDéN B,DAUNIUS O.Turbulent flow and heat transfer in a square-sectioned U-bend[J].Progress in Computational Fluid Dynamics,an International Journal,2006,6(1):89-100.

[10]SHIN J,CHOI Y,AN J.Numerical analysis of turbulent flow and heat transfer in a square sectioned U-bend duct by elliptic-blending second moment closure[J].Journal of Mechanical Science & Technology,2007,21(2):360-371.

[11]NOBARI M R H,NOUSHA A,DAMANGIR E.A numerical investigation of flow and heat transfer in rotating U-shaped square ducts[J].International Journal of Thermal Sciences,2009,48(3):590-601.

[12]MOON M A,KIM K Y.Multi-objective optimization of a guide vane in the turning region of a rotating U-duct to enhance heat transfer performance[J].Heat & Mass Transfer,2012,48(11):1941-1954.

[13]SALAMEH T,SUNDEN B.A numerical investigation of heat transfer in a smooth bend part of a U-duct[J].International Journal of Numerical Methods for Heat & Fluid Flow,2013,24(1):137-147.

[14]SHEN Z,ZHANG D,XIE Y.Heat transfer enhancement in gas turbine U-shaped channel with rib tabulators under rotational effect[J].Energy Procedia,2014,61:2501-2504.

[15]MENTER F R,KUNTZ M,LANGTRY R.Ten years of industrial experience with the SST turbulence model[C]//Proceedings of the Fourth International Symposium on Turbulence,Heat and Mass Transfer.Antalya,Turkey:Begell House,2003:625-632.

[16]RHIE C M,CHOW W L.Numerical study of the turbulent flow past an isolate airfoil with trailing edge separation[J].AIAA Journal,1982,21:1525-1532.

[17]陶文銓.數值傳熱學的近代進展[M].北京：科學出版社,2000.

[18]LAI H X,YAN Y Y,GENTLE C R.Calculation procedure for conjugate viscous flows about and inside single bubbles[J].Numerical Heat Transfer,Part B,Fundamentals:An International Journal of Computation and Methodology,2003,43(3):241-265.

[19]LAI H X,XING G L,TU S T,etal.A pressure-correction procedure with high-order schemes[J].International Journal of Numerical Methods for Heat and Fluid flow,2011,21(3):331-352.

[20]馬坤,楊金龍,賴煥新.大曲率彎管中湍流的計算與模型考證[J].華東理工大學學報(自然科學版),2013,39(5):617-624.

[21]KADER B A.Temperature and concentration profiles in fully turbulent boundary layers[J].International Journal of Heat & Mass Transfer,1981,24(9):1541-1544.

Numerical Simulation of Turbulent Convective Heat Transfer in a U-Bend

XU Wen-jie,LAI Huan-xin

(Key Laboratory of Pressurized Systems and Safety,Ministry of Education,East China University of Science and Technology,Shanghai 200237,China)

Abstract:The grid-independent solution is a premise of numerically studying turbulent convective heat transfer,while the appropriate treatment of boundary conditions is critical to the accuracy of numerical results.In this paper,turbulent heat transfer in a squared U-bend is studied with the SST k-ω model.The flow and temperature fields are calculated,and the grid-independent solutions of velocity and temperature fields are analyzed emphatically.Under the second type of boundary condition (heat flux condition),the influence of two treatments on the accuracy of numerical solutions is contrastively analyzed.Because the coupling of velocity and temperature fields is weak,the mesh size required for grid-independent solutions is different for the two fields.A finer mesh is needed for the grid-independent solution of the temperature field.The comparative analysis of the two treatments of the heat flux condition demonstrates that the algebraic formulation of Kader (1981) is in good agreement with experimental data even though the mesh is relatively coarse,so it has better versatility and grid adaptability.

Key words:grid-independent solution; turbulent; convective heat transfer; boundary condition

收稿日期：2015-06-10

基金項目：國家自然科學基金(51176048)

作者簡介：徐文潔(1991-),女,碩士生,主要從事彎道對流換熱的數值模擬研究。E-mail:cassie_xuxu@163.com 通信聯系人：賴煥新,E-mail:hlai@ecust.edu.cn

文章編號：1006-3080(2016)02-0281-10

DOI:10.14135/j.cnki.1006-3080.2016.02.021

中圖分類號：O357.5

文獻標志碼：A