二維后向臺階層流傳熱特性的數(shù)值模擬

仲敏波，謝緯安，袁銀男，，喜冠南，

(1.江蘇大學(xué)能源與動力工程學(xué)院，江蘇鎮(zhèn)江212013;2.南通大學(xué)機(jī)械工程學(xué)院，江蘇南通226019)

在能源、動力、機(jī)械過程的許多管內(nèi)換熱設(shè)備中，通過改變流道截面，使流體發(fā)生分離，能達(dá)到強(qiáng)化換熱的目的［1］.后向臺階雖然幾何結(jié)構(gòu)簡單，但因?yàn)榉蛛x點(diǎn)固定，而且包含著分離流動的一系列特征，常被作為研究分離流動的基準(zhǔn)模型之一.

迄今，針對后向臺階流動已經(jīng)有了大量的試驗(yàn)［2-5］和數(shù)值研究［6-10］.B.F.Armaly 等［2］首次對同一模型同時進(jìn)行了試驗(yàn)和數(shù)值研究，通過對比發(fā)現(xiàn)層流狀態(tài)下的二維數(shù)值模擬結(jié)果和試驗(yàn)結(jié)果吻合，而隨著雷諾數(shù)(Re)增大，由于臺階側(cè)壁的影響，數(shù)值模擬結(jié)果偏離試驗(yàn)結(jié)果.之后，眾多的科研工作者對這一復(fù)雜的流動結(jié)構(gòu)進(jìn)行了研究.然而，現(xiàn)有的關(guān)于后向臺階的研究大部分只針對流動特性和相關(guān)的流場，只有少部分是針對后向臺階的傳熱特性.H.Chaudhary等［11］通過二維數(shù)值模擬研究了層流(Re=1～200)時，Re和Pr的變化對流動和傳熱特性的影響;孫正中等［12］研究了湍流時臺階下壁面溫度變化對通道內(nèi)部溫度分布的影響.筆者通過數(shù)值模擬方式研究層流狀態(tài)下，雷諾數(shù)和臺階幾何參數(shù)(臺階高度、平板間高度)的變化對臺階底面?zhèn)鳠崽匦缘挠绊?，并結(jié)合流場中流動結(jié)構(gòu)的變化對傳熱特性的分布規(guī)律進(jìn)行分析.這對指導(dǎo)層流狀態(tài)下管內(nèi)換熱設(shè)備的結(jié)構(gòu)優(yōu)化，提高換熱設(shè)備的換熱效率具有重要意義.

1 模型的建立和驗(yàn)證

1.1 控制方程

數(shù)值模擬研究流體流過二維后向臺階后的流動傳熱情況.以連續(xù)、不可壓縮流體(空氣)為研究對象，所采用的控制方程包括質(zhì)量守恒、動量守恒和能量守恒方程:

質(zhì)量守恒方程為

式中:u，v分別為x，y方向速度分量;ρ為流體密度.

動量守恒方程為

式中:p為壓力;μ為黏性系數(shù).

能量守恒方程為

式中:cp為定壓比熱;h為比焓;λ為熱導(dǎo)率.

1.2 計算區(qū)域和邊界條件

二維后向臺階流動的計算區(qū)域如圖1所示.

圖1 后向臺階計算區(qū)域

采用直角坐標(biāo)系建立數(shù)值模型，坐標(biāo)原點(diǎn)O位于臺階垂直壁面和下游底面的交點(diǎn)處，x方向?yàn)檠刂黧w流動方向，y方向?yàn)榇怪庇诹黧w流動方向.S為臺階高度，H為平板間高度，L0為臺階上游進(jìn)口長度，L為臺階下游出口長度，ER=H/(H-S)為臺階擴(kuò)張比.數(shù)值計算設(shè)定的初始條件:進(jìn)口邊界處的溫度均勻分布θin=10℃，進(jìn)口速度uin充分發(fā)展，滿足拋物線分布;出口邊界處的速度場和溫度場符合邊界層近似理論.壁面邊界:臺階所有壁面速度無滑移，臺階下游底面被等溫加熱θw=40℃，其余壁面絕熱.

1.3 網(wǎng)格劃分和模型驗(yàn)證

采用非均勻網(wǎng)格劃分計算區(qū)域，由于流動參數(shù)在臺階分離點(diǎn)處和臺階上、下平板壁面處的變化梯度較大，為了保證計算結(jié)果的準(zhǔn)確性，需要對上述區(qū)域附近的網(wǎng)格加密.對計算區(qū)域劃分網(wǎng)格時，保證臺階分離點(diǎn)處的網(wǎng)格間距為最小，保證流向和法向的首個、末個網(wǎng)格間距分別為相鄰網(wǎng)格間距的0.5倍.

數(shù)值模擬采用課題組自行編制的FORTRAN計算程序.求解控制方程時，首先采用基于交錯網(wǎng)格的有限容積法離散控制方程，離散過程中動量方程的對流項和擴(kuò)散項分別采用三階迎風(fēng)差分(QUICK算法)和二階中心差分格式，時間項采用隱式求解法.對離散后的方程求解時，每一個時間步采用ADI算法，壓力場的修正采用SIMPLE法則.

為了驗(yàn)證數(shù)值計算方法的可靠性，對相同工況下計算得到的流體再附著點(diǎn)位置xr與前人相關(guān)的試驗(yàn)和數(shù)值模擬結(jié)果進(jìn)行了對比，結(jié)果如圖2所示，按以上數(shù)值模型計算的結(jié)果與文獻(xiàn)［6］中ER=2.00的數(shù)值模擬數(shù)據(jù)表現(xiàn)一致，與文獻(xiàn)［2］中ER=1.94的試驗(yàn)數(shù)據(jù)之間存在一個差異(隨著Re的增大，數(shù)值模擬結(jié)果低于試驗(yàn)結(jié)果).通過驗(yàn)證表明文中所用數(shù)值計算方法是正確可靠的.

圖2 再附著點(diǎn)位置比較

2 計算結(jié)果與討論

2.1 臺階底面努塞爾數(shù)和表面摩擦因數(shù)的分布

文獻(xiàn)［2］以ER=1.94的臺階作為研究對象，以臺階進(jìn)口高度(Hin=H-S)的2倍(2Hin)作為特征高度，得到Re2Hin＜1 200時，流場處于層流狀態(tài).將這一范圍換算到以臺階高度(S)作為特征高度，對應(yīng)的層流Re數(shù)范圍為ReS＜565.所以，文中的研究范圍定為Re=100～400.

臺階擴(kuò)張比ER=2.00，臺階上游進(jìn)口長度L0=S，下游出口長度L=100S.臺階高度S作為特征量，相應(yīng)的Re和Nu分別定義為

式中ˉuin為流體進(jìn)口平均速度.

式中qw為熱流密度.

流場穩(wěn)定后臺階下游底面努塞爾數(shù)Nu的分布情況如圖3所示.從圖中可以觀察到以下特征:①Re一定時，沿著流體流動方向，Nu先迅速增大，到達(dá)峰值后又迅速減小，沿程只有1個峰值，在主回流區(qū)，Nu隨著距臺階距離的增加而快速增加，到達(dá)再附著點(diǎn)時，Nu達(dá)到峰值，換熱效率達(dá)到最大，在再附著點(diǎn)下游，Nu隨著距臺階距離的增加先快速減小，然后慢慢減小;②隨著Re的增大，Nu的變化范圍也越大，峰值越大且峰值位置隨著Re的增加向下游移動，Re從100增加到400，峰值位置由x/S=5處變化到x/S=12處;③Re≥300時，Nu變化曲線在下降段有個明顯的曲折點(diǎn)，且曲折點(diǎn)的位置穩(wěn)定在x/S=22處，在曲折點(diǎn)下游Nu受Re的影響不大.

圖3 ER=2.00臺階底面努塞爾數(shù)的分布

對應(yīng)工況下，流場穩(wěn)定后臺階底面摩擦因數(shù)Cf的變化規(guī)律如圖4所示.

圖4 ER=2.00臺階底面表面摩擦因數(shù)的分布

從圖4可以觀察到以下特征:①Re分別為100，200時，Cf沿著流動方向，先由0迅速減小，到達(dá)最小值后又迅速增大到0，然后繼續(xù)向正值方向增大到最大值，之后Cf的變化比較平緩;②Re≥300時，Cf到達(dá)最大值后又出現(xiàn)下降.Re越大，Cf下降也越多;③ 隨著Re的增大，Cf的變化范圍減小，再附著點(diǎn)的位置(Cf=0)向下游移動.

2.2 流場與溫度場

Re變化時ER=2.00的后向臺階流動達(dá)到穩(wěn)定后內(nèi)部流場的變化情況圖5所示，Re為100和200時，流場只在靠近臺階處存在一個主回流區(qū)域.當(dāng)Re≥300時，位于主回流區(qū)域下游靠近臺階上平板壁面又出現(xiàn)了1個二次回流區(qū)域.二次回流區(qū)域的出現(xiàn)，擠壓了流體的中間流道，使流體的速度增加，而壓力減小，對應(yīng)位置的Cf出現(xiàn)下降，如圖4所示.流線和速度矢量圖如圖5所示，標(biāo)“▲”處表示再附著點(diǎn)位置，Re越大，主回流區(qū)長度越長.

隨著Re增大溫度分布云圖的變化情況如圖6所示，圖中的曲線表示等溫線.流體在臺階角點(diǎn)發(fā)生分離，經(jīng)過一段距離后，流體和壁面重新接觸，在再附著點(diǎn)處(圖中標(biāo)“▲”的位置)，由于流體沖擊壁面，壁面溫度邊界層破壞，溫度變化梯度最大，圖3中底面Nu在對應(yīng)位置達(dá)到最大.

圖5 流線和速度矢量圖 (ER=2.00)

圖6 溫度分布云圖(ER=2.00)

當(dāng)Re≥300，臺階上壁面出現(xiàn)二次回流區(qū)后，底面Nu變化曲線在下降段的曲折點(diǎn)固定.圖7，8比較了Re分別為300，350，400時后向臺階流道剖面上的量綱一速度、溫度沿法向的分布.選擇6個流道剖面，其中，剖面x/S=0，5，10，16位于曲折點(diǎn)上游，剖面x/S=22位于曲折點(diǎn)處，剖面x/S=28位于曲折點(diǎn)下游.從圖7，8可以看出:由于進(jìn)口速度不同，曲折點(diǎn)上游各剖面上的速度和溫度沿法向的變化并不一致;在剖面x/S=22處，速度和溫度在靠近底面處的變化梯度表現(xiàn)一致;在剖面x/S=28處，速度和溫度沿法向的變化大致相同.

圖7 法向速度分布 (ER=2.00)

圖8 法向溫度分布 (ER=2.00)

2.3 平板間高度變化對臺階底面Nu分布影響

保持臺階高度S不變，改變平板間高度H，圖9，10分別表示ER=1.50，3.00 時臺階底面Nu隨Re的變化.

圖9 ER=1.50臺階底面努塞爾數(shù)的分布

圖10 ER=3.00臺階底面努塞爾數(shù)的分布

從圖9，10可以看到，ER一定時，Nu的變化規(guī)律和圖3相似，但隨著ER的增大，Re的變化對Nu峰值的影響也增大.ER=1.50，2.00，3.00 時，Re從100 增大到400，Nu的峰值分別增加了約 1.3，2.0，2.6倍.圖9中ER=1.50時，由于Re≤400時臺階下游上壁面還未出現(xiàn)二次回流區(qū)，Nu變化曲線在下降段未出現(xiàn)固定的曲折點(diǎn);圖10中ER=3.00時，由于Re≥300時臺階下游上壁面出現(xiàn)了二次回流區(qū)，Nu變化曲線出現(xiàn)了1個固定的曲折點(diǎn).Re≥300后，ER由 2.00增大到 3.00，該曲折點(diǎn)的位置由x/S=22處向上游移動至x/S=15處，但曲折點(diǎn)處Nu僅由1.38增加到1.41.所以臺階上壁面出現(xiàn)二次回流區(qū)域后，增大ER能提高曲折點(diǎn)上游的Nu，但對曲折點(diǎn)下游的Nu影響較小.

臺階進(jìn)口速度保持0.190 8 m·s-1不變，Re=200，S=15.0 mm，H=22.5，25.0，30.0，45.0，60.0 mm時，臺階底面Nu的分布如圖11所示.

圖11 相同S，不同H下底面努塞爾數(shù)的分布

從圖11可以看出:H由60.0減小到22.5，ER由1.33增加到3.00，Nu的變化范圍隨H減小而增大，Nu峰值的位置隨H減小向上游移動.

2.4 臺階高度變化對臺階底面Nu分布影響

保持平板間高度H不變，改變臺階高度S.以H作為特征量，Re和Nu分別記為ReH和NuH，定義為

臺階進(jìn)口速度保持0.190 8 m·s-1不變，ReH=400，H=30.0 mm，S=7.5，10.0，15.0，18.0，20.0 mm時，臺階底面NuH的分布如圖12所示，S由7.5增加到20.0，ER由 1.33 增加到 3.00.NuH的變化情況分2部分:S由7.5增加到15.0，ER由1.33增加到2.00，NuH的峰值基本不變，但峰值位置隨著S增加向下游移動;S由15.0增加到20.0，ER由2.00增加到3.00，NuH的峰值快速增加，但峰值的位置隨著S增加變化不大，位于x/H=4附近.

圖12 相同H，不同S下底面努塞爾數(shù)的分布

3 結(jié)論

1)數(shù)值模型計算得到的流體再附著點(diǎn)位置與現(xiàn)有的試驗(yàn)和計算結(jié)果吻合，說明計算模型可靠，計算結(jié)果正確有效.

2)隨著Re的增大，臺階上平板壁面出現(xiàn)二次回流區(qū)后，Nu變化曲線在下降段有一明顯的曲折點(diǎn).Re繼續(xù)增大，曲折點(diǎn)的位置穩(wěn)定，在曲折點(diǎn)下游Nu受Re的影響不大.臺階上壁面出現(xiàn)二次回流區(qū)域后，增大ER能提高曲折點(diǎn)上游的Nu，但對曲折點(diǎn)下游的Nu影響較小.

3)保持S不變，減小H，ER由1.33增加到3.00，Nu的變化范圍增大，Nu峰值的位置向上游移動.

4)保持H不變，增大S.ER由1.33增加到2.00，NuH峰值基本不變，但峰值的位置隨著ER增加向下游移動;ER由2.00增加到3.00，NuH峰值快速增加，但峰值的位置隨著ER增加基本不變.

References)

［1］蘇 磊，張 紅，丁雷江，等.單回路脈動熱管傳熱性能的試驗(yàn)［J］.江蘇大學(xué)學(xué)報:自然科學(xué)版，2012，33(1):49-54.Su Lei，Zhang Hong，Ding Leijiang，et al.Experiment on heat transfer characteristics of single loop pulsating heat pipe［J］.Journal of Jiangsu University:Natural Science Edition，2012，33(1):49-54.(in Chinese)

［2］Armaly B F，Durst F，Pereira J C F，et al.Experimental and theoretical investigation of backward-facing step flow［J］.J Fluid Mech，1983，127:473-496.

［3］齊鄂榮，黃明海，李 煒，等.二維后向臺階流流動特性的實(shí)驗(yàn)研究［J］.實(shí)驗(yàn)力學(xué)，2006，21(2):225-232.Qi Erong，Huang Minghai，Li Wei，et al.A experimental study on the 2D time-average flow over a backward facing step via PIV ［J］.Journal of Experimental Mechanics，2006，21(2):225-232.(in Chinese)

［4］Tihon J，Pěnkavová V，?imík H M，et al.The transitional backward-facing step flow in a water channel with variable expansion geometry［J］.Experimental Thermal and Fluid Science，2012，40:112-125.

［5］Chen Z，Yi S H，Tian L F，et al.Flow visualization of supersonic laminar flow over a backward-facing step via NPLS［J］.Shock Waves，2013，23:299-306.

［6］Iwai H，Nakabe K，Suzuki K.Flow and heat transfer characteristics of backward-facing step laminar flow in a rectangular duct［J］.International Journal of Heat and Mass Transfer，2000，43(3):457-471.

［7］Barri M，El Khoury G K，Andersson H I，et al.DNS of backward-facing step flow with fully turbulent inflow［J］.International Journal for Numerical Methods in Fluids，2010，64(7):777-792.

［8］Malamataris N A，L?hner R.The computation of the eddy along the upper wall in the three-dimensional flow over a backward-facing step［J］.International Journal for Numerical Methods in Fluids，2012，68(9):1102-1125.

［9］Malamataris N A.A numerical investigation of wall effects in three-dimensional，laminar flow over a backward facing step with a constant aspect and expansion ratio［J］.International Journal for Numerical Methods in Fluids，2013，71(9):1073-1102.

［10］肖 瀟，吳時強(qiáng)，樊新建，等.后臺階流的水動力特征［J］.江蘇大學(xué)學(xué)報:自然科學(xué)版，2014，35(2):149-153.Xiao Xiao，Wu Shiqiang，F(xiàn)an Xinjian，et al.Hydraulic characteristics of backward-facing step flow［J］.Journal of Jiangsu University:Natural Science Edition，2014，35(2):149-153.(in Chinese).

［11］Chaudhary H，Dhiman A.Two-dimensional laminar fluid flow and heat transfer over a backward-facing step:effects of reynolds and prandtl numbers［J］.Heat Transfer Research，2011，42(4):379-402.

［12］孫正中，蘇莫明，周 錚.平面不可壓后臺階湍流流動的數(shù)值模擬:壁面溫度的影響［J］.航空計算技術(shù)，2009，39(4):71-74.Sun Zhengzhong，Su Moming，Zhou Zheng.Numerical simulation of 2-D incompressible turbulent flow over a backward-facing step:effects of wall temperature［J］.Aeronautical Computing Technique，2009，39(4):71-74.(in Chinese)