內(nèi)置發(fā)熱圓的方腔內(nèi)雙擴散混合對流數(shù)值模擬研究

肖瑞雪， 李貝貝， 徐洪濤， 楊 茉

（上海理工大學能源與動力工程學院，上海 200093）

在自然界和工程應(yīng)用中，許多過程同時受到熱擴散和質(zhì)擴散共同作用，如噴霧和閃蒸過程、霧化液體燃料的燃燒、旋風蒸發(fā)和干燥過程、化工和食品加工中的脫水及晶體生長等［1－6］，這些問題可歸結(jié)為雙擴散混合對流問題.國內(nèi)外許多學者對雙擴散混合對流的問題進行了模擬研究.陳文等［7］研究了室內(nèi)雙擴散混合對流輸運過程，分析了送風強度、熱源強度及污染源強度對室內(nèi)空氣流動、熱量與污染物輸運過程的影響，并給出合理的控制室內(nèi)溫度和污染源傳播的途徑.Teamah等［8］對頂蓋驅(qū)動的矩形腔內(nèi)雙擴散混合對流進行了研究，分析了頂蓋向左右兩個方向移動的情況，得出了兩種移動情況下隨著理查德數(shù)Ri的減小，傳熱和傳質(zhì)都減弱等結(jié)論.Al－Amiri等［9］對頂蓋驅(qū)動腔體內(nèi)雙擴散過程進行了研究，指出低Ri時，頂蓋驅(qū)動作用使得腔體內(nèi)的傳熱和傳質(zhì)效應(yīng)增強.Basak等［10］研究了底部均勻和不均勻加熱兩種情況下方腔內(nèi)的混合對流情況，分析了不同雷諾數(shù)Re、普朗特數(shù)Pr以及格拉曉夫數(shù)Gr的影響，并對努賽爾數(shù)Nu的變化情況進行了分析.Nayak等［11］對頂蓋驅(qū)動的立方腔內(nèi)雙擴散混合對流進行了研究，頂蓋以恒定速度移動，并且溫度和濃度梯度在腔體外部，分析了二維和三維情況下熱壁面平均努賽爾數(shù)Nuav與平均舍伍德數(shù)Shav隨著Re的變化情況.Kandaswamy等［12］分析了密度最大值對變壁溫多孔介質(zhì)腔內(nèi)雙擴散自然對流的影響.Sivasankaran等［13］研究了多孔腔內(nèi)不定常密度流體的雙擴散自然對流，給出了速度矢量圖、等溫線圖、流線圖以及等濃度線分布圖，并分析了Nuav及Shav的變化情況.Nishimura等［14］研究了壁面具有水平溫度和濃度梯度的矩形腔內(nèi)振蕩雙擴散自然對流，討論了浮升力比對雙擴散對流的影響以及模型的震蕩性問題.Chamkha等［15］對壁面具有相反的溫度和濃度梯度的矩形腔體內(nèi)磁流體雙擴散自然對流進行了研究.Ching等［16］采用有限元法模擬分析了在Le一定時，Br，Ri以及移動壁的方向?qū)φ乔粌?nèi)傳熱和傳質(zhì)的影響.Chamkha等［17］研究了充滿均勻多孔介質(zhì)的傾斜長方形腔體內(nèi)雙擴散不穩(wěn)定層流，對腔體內(nèi)振蕩流動進行了預測并且觀察到了衰減振蕩.Sun等［18］建立了一個微可壓縮流模型研究封閉腔體內(nèi)二元混合理想氣體的雙擴散對流，分析了變密度對充滿二元混合理想氣體的垂直腔內(nèi)自然對流的影響.Hasanuzzama等［19］研究了Le對三角形腔體內(nèi)傳熱和傳質(zhì)的影響.

從上述文獻綜述可知，目前大部分研究集中于頂蓋驅(qū)動或側(cè)壁移動情況下的雙擴散混合對流，較少研究內(nèi)部放置加熱圓的方腔內(nèi)雙擴散混合對流問題，然而這是一個典型的污染物處理、熱舒適性研究以及煙氣控制等的簡化模型.因此，本文對內(nèi)部放置加熱圓的方腔內(nèi)二元雙擴散混合對流現(xiàn)象進行了數(shù)值模擬研究，分析了Ri，Le，Br對方腔內(nèi)雙擴散混合對流的影響，給出了不同情況下的等溫度線、等濃度線以及流線圖，并分析了加熱圓表面的Nuav和Shav的變化規(guī)律.

1 計算模型和方法

1.1 物理模型

如圖1所示，方腔內(nèi)置一加熱圓，加熱圓位于方腔中心.方腔的長寬均為L，圓的直徑為d，其中d＝0.4L.方腔設(shè)有開口，開口的大小為0.1L.流體從左側(cè)底部流進，右側(cè)頂部流出，進口速度為ui，溫度為Ti，濃度為ci，圓表面的溫度為Th，濃度為ch，方腔其余壁面絕熱且濃度梯度為零，重力加速度為g，其中Th＞Ti，ch＞ci.方腔內(nèi)流體Pr＝0.7.

圖1 物理模型Fig.1 Schematic diagram of thephysical model

1.2 數(shù)學描述

在進行方腔內(nèi)對流換熱和傳質(zhì)的數(shù)值計算時，為便于處理由于溫差和濃度差而引起浮升力項，常采用Boussinesq假設(shè)［20］：

其中，ρ為密度；ρ0為與Τ0相對應(yīng)的流體密度，kg／m3；βt為溫度引起的體積膨脹系數(shù)為濃度引起的體積膨脹系數(shù)

定義下列無量綱參數(shù)

式中，u，v分別為直角坐標系下水平方向和豎直方向的速度分量；T，p，c分別為溫度、壓力和濃度；α，D，υ，g，為熱擴散系數(shù)、質(zhì)擴散系數(shù)、運動黏度、重力加速度.則圖1所示的方腔內(nèi)雙擴散混合對流問題在直角坐標系下的無量綱數(shù)學描述方程為

無量綱邊界設(shè)定為

式中，n為表面法向.圓表面的局部Nu和局部Sh，可表示為

式中，h，hm為圓表面的局部換熱系數(shù)和傳質(zhì)系數(shù).圓表面的Nuav和Shav可表示為

1.3 計算方法

控制方程采用有限體積法［21］求解，方程采用二階迎風格式離散.數(shù)值計算采用SIMPLE（semiimplicit method for pressure linked equations）算法進行壓力速度耦合計算.為確保計算的準確性，進行了網(wǎng)格獨立性驗證，網(wǎng)格數(shù)分別取50×50，60×60，70×70，80×80，90×90和100×100，如圖2所示.在網(wǎng)格數(shù)為100×100的情況下，數(shù)值解基本穩(wěn)定.因此，最終選定網(wǎng)格數(shù)為100×100.數(shù)值計算時，連續(xù)性方程、動量方程和傳質(zhì)方程的計算殘差取10－3，能量方程的計算殘差取10－6.

圖2 不同網(wǎng)格下的NuavFig.2 Grid independency check for the average Nusselt number around the heated cylinderat Ri＝0.01，Br＝10and Le＝1

為了確保計算方法的準確性，本文對頂蓋驅(qū)動的矩形腔內(nèi)雙擴散混合對流進行了研究，并將計算所得結(jié)果與文獻［8］中的結(jié)果作了比較，如圖3和圖4所示.從圖中可知目前的數(shù)值方法得出的結(jié)果與文獻中的基本吻合，計算方法可行.

圖3 等濃度線、流線和等溫度線分布圖（Pr＝0.7，Ri＝1，Br＝1，Le＝10）Fig.3 Isoconcentrations，streamlines and isotherms at Pr＝0.7，Ri＝1，Br＝1and Le＝10

圖4 本文與文獻［8］中Shav的比較Fig.4 Comparison of average Sherwood numbers

2 模擬結(jié)果分析

2.1 路易斯數(shù)Le的影響

Le是確定質(zhì)量擴散和熱擴散相對關(guān)系的重要參數(shù)之一，本文研究了Le在0.1與10范圍內(nèi)，方腔內(nèi)的雙擴散混合對流情況，圖5給出了Ri＝1及 Br＝1時不同Le下等濃度線、流線和等溫線分布情況.

圖5 不同Le下等濃度線、流線以及等溫線分布情況（Pr＝0.7，Ri＝1，Br＝1）Fig.5 Effects of Lewis number on isoconcentrations，streamlines and isotherm at Pr＝0.7，Ri＝1 and Br＝1

從圖5可以看出流線和等溫線的分布基本不隨Le的變化而變化，而在發(fā)熱圓表面的質(zhì)邊界層隨著Le數(shù)的增加逐漸變薄，使得傳質(zhì)速率增加.Le＝0.1時，等濃度線均勻分布在方腔左下角進口附近，方腔右上角濃度梯度幾乎為零.隨著Le的增大，等濃度線由進口逐漸向發(fā)熱圓集中.當Le＝10時，等濃度線大多集中在發(fā)熱圓周圍及出口附近，進口的濃度梯度幾乎為零，加熱圓上部的質(zhì)羽流形狀明顯.從圖5流線分布圖可以看出，在進口附近有一個漩渦，流線部分從加熱圓左側(cè)繞過，部分從加熱圓右側(cè)流過，這是因為Ri＝1時自然對流與強制對流作用相當.從圖5等溫線分布圖可以看出等溫線在發(fā)熱圓和兩側(cè)壁之間分層分布，在發(fā)熱圓左下角等溫線分布較密集.

圖6 不同Le下Nuav和Shav隨Ri的變化情況Fig.6 Average Nusselt and Sherwood numbers

圖6給出了Br＝1時，不同Le下平均努賽爾數(shù)和平均舍伍德數(shù)隨Ri的變化情況.從圖6（a）可以看出同一Ri下，Nuav幾乎不隨Le的變化而變化.Nuav隨著Ri的增加而逐漸減小，這是因為隨著Ri的增加，強制對流作用減弱，從而減弱了整體換熱效果.從圖6（b）可以看出當Le≤1時，Shav隨著Ri的增加而逐漸減??；當Le≥5，Shav隨著Ri的增加先增加而后減小，因為低Le下Ri的影響占主導地位.從圖6（b）還可以看出Shav隨著Le的增加而增加，這是因為Le的增加使得傳質(zhì)效果增強.

2.2 浮升力比Br的影響

浮升力比Br也是影響方腔內(nèi)雙擴散混合對流的重要因素之一，不同的Br使得熱邊界層增厚或減薄.圖7給出了Ri＝1和Le＝0.1時不同Br下等濃度線、流線和等溫線分布圖.Le＝0.1表明質(zhì)擴散比熱擴散強.從圖中可以看出，在Le＝0.1時，等濃度線幾乎不隨Br的變化而變化，并且在進口附近均勻傾斜分層分布，圖中可以觀察到，數(shù)值為0.95的等濃度線隨著Br的增大，逆時針偏轉(zhuǎn).Br＝－10時，質(zhì)浮升力占主導地位，且與熱浮升力的方向相反.從圖中可以看出，流線全部從發(fā)熱圓左側(cè)流過.隨著Br的增加，大部分流線從發(fā)熱圓左上角向右下角轉(zhuǎn)移.當Br≥0時，流線在進口附近形成一個漩渦，并隨著Br的增加逐漸增大，而等溫線幾乎不發(fā)生變化.

圖7 不同Br下等濃度線、流線和等溫線分布圖（Pr＝0.7，Ri＝1，Le＝0.1）Fig.7 Effects of buoyancy ratio on isoconcentrations，streamlines and isotherms at Pr＝0.7，Ri＝1 and Le＝0.1

圖8給出了強制對流Ri＝0.01時，不同Le下平均努賽爾數(shù)Nuav和平均舍伍德數(shù)Shav隨Br的變化情況.浮升力比的增大將增強自然對流作用，減弱強制對流作用.從圖8（a）可以看出在Le＝0.1時，Nuav不隨Br的增大而變化，即浮升力比對低路易斯數(shù)的影響較小.其它情況下，Nuav隨著Br的增大而逐漸減小.當Br從－10變化到10時，Nuav減小的速率隨著Le的增大而增大.當Br＝0時，從方程（4）可以看出，傳質(zhì)對流動沒有影響，因而不同Le下Nuav相同.當Br＞0時，Y方向的動量方程中源項增大，Le的增大增強了傳質(zhì)速率，平均努賽爾數(shù)Nuav隨著Le的增大而減小.相反，當Br＜0時，平均努賽爾數(shù)Nuav隨著Le的增大而增大.同理，浮升力比的增大增強了自然對流作用，減弱了強制對流作用，從圖8（b）可以很明顯地看出，在0.1≤Le≤1時，Shav幾乎不隨Br的變化而變化，只是因為此范圍內(nèi)的質(zhì)傳遞速率較大.當Le＝5，10和15時，Shav受Br的影響較明顯且隨Br的增大而減小.同一Br下，Shav隨著Le的增大而增大.

圖8 不同Le下Nuav和Shav隨Br的變化情況Fig.8 Average Nusselt and Sherwood number

為了更準確地了解加熱圓不同位置處的傳熱傳質(zhì)情況，更好地反映Br對雙擴散混合對流的影響.圖9給出了在Le＝1和Ri＝1時，不同Br下發(fā)熱圓表面局部Sh的變化情況.從圖中可以看出，當Br≥0時，局部Sh的最大值在φ＝240°處.當Br＝－10時，Sh的最大值在90°＜φ＜120°之間.

圖9 發(fā)熱圓表面局部Sh（Le＝1，Ri＝1）Fig.9 Local Sherwood number around the heated cylinder at Le＝1and Ri＝1

3 結(jié) 論

對內(nèi)置發(fā)熱圓的方腔內(nèi)雙擴散混合對流問題進行了數(shù)值模擬分析，得出以下結(jié)論：

a.當Ri＝1和Br＝1時，流線和等溫線的分布基本不隨Le的變化而變化，等濃度線隨著Le的增大由進口逐漸向發(fā)熱圓和出口集中.當Br＝1時，Nuav幾乎不隨Le的變化而變化，隨著Ri的增加而逐漸減小，Shav隨著Le的增加而增加.當Br＝1和Le≤1時，Shav隨著Ri的增加而逐漸減?。划擝r＝1和Le≥5，Shav隨著Ri的增加先增加而后減小.

b.當Ri＝1和Le＝0.1時，等濃度線幾乎不隨Br的變化而變化，并且在進口附近呈傾斜分層分布.隨著Br的增加大部分流線從發(fā)熱圓左上角向右下角轉(zhuǎn)移.當Br≥0時，流線在進口附近形成一個漩渦，并隨著Br的增加逐漸增大，而等溫線幾乎不發(fā)生變化.

c.當Ri＝0.01時，Le＝0.1時的Nuav幾乎不隨著Br的變化而變化，其它情況下Nuav隨著Br的增大而逐漸減小，當Br從－10變化到10時，Nuav減小的速率隨著Le的增大而增大.當Br＞0時，Nuav隨著Le的增大而減小，當Br＜0時，Nuav隨著Le的增大而增大.在0.1≤Le≤1時，Shav幾乎不隨Br的變化而變化，只是因為此范圍內(nèi)的質(zhì)傳遞速率較大.當Le＝5，10和15時，Shav受Br的影響較明顯且隨Br的增大而減小.同一Br下，Shav隨著Le的增大而增大.

［1］ Marcos J D，Izquierdo M，Parra D.Solar space heating and cooling for Spanish housing：potential energy savings and emissions reduction［J］.Solar Energy，2011，85（11）：2622－2641.

［2］ Serrano－Arellano J，Xamán J，álvarez G.Optimum ventilation based on the ventilation effectiveness for temperature and CO2distribution in ventilated cavities［J］.International Journal of Heat and Mass Transfer，2013，62：9－21.

［3］ Ghaddar N，Ghali K，Chehaitly S.Assessing thermal comfort of active people in transitional spaces in presence of air movement［J］.Energy and Buildings，2011，43（10）：2832－2842.

［4］ Vijaya Venkata Raman S，Iniyan S，Goic R.A review of solar drying technologies［J］.Renewable and Sustainable Energy Reviews，2012，16（5）：2652－2670.

［5］ Yang M，Ma N.Free convection in a liquidencapsulated molten semiconductor in a vertical magnetic field［J］.International Journal of Heat and Mass Transfer，2005，48（19／20）：4010－4018.

［6］ Hasan M，Mujumdar A S.Simultaneous mass and heat transfer undermixed convection along a vertical cone［J］.International Journal of Energy Research，1985，9（2）：129－140.

［7］ 陳文，趙福云，湯廣發(fā)，等.耦合墻體擴散的室內(nèi)雙擴散混合對流輸運過程［J］.暖通空調(diào)，2006，36（8）：12－18.

［8］ Teamah M A，El－Maghlany W M.Numerical simulation of double－diffusive mixed convective flow in rectangular enclosure with insulated moving lid［J］.International Journal of Thermal Sciences，2010，49：1625－1638.

［9］ Al－Amiri A M，Khanafer K M，Pop I.Numerical simulation of combined thermal and mass transport in a square lid－driven cavity［J］.International Journal of Thermal Sciences，2007，46（7）：662－671.

［10］ Basak T，Roy S，Sharma P K，et al.Analysis of mixedconvection flows within a square cavitywithuniform and non－uniform heating of bottom wall［J］.International Journal of Thermal Sciences，2009，48（5）：891－912.

［11］ Nayak A K，Bhattacharyya S.Double－diffusive convection in a cubical lid－driven cavitywith opposing temperature and concentration gradients［J］.Theoretical and Computational Fluid Dynamics，2012，26（6）：565－581.

［12］ Kandaswamy P，Eswaramurthi M，Lee J.Density maximum effect on double－diffusive natural convection in a porous cavity with variable wall temperature［J］.Transport in Porous Media，2008，73（2）：195－210.

［13］ Sivasankaran S，Kandaswamy P，Ng CO.Double diffusive convection of anomalous density fluidsin a porous cavity［J］.Transport in Porous Media，2008，71（2）：133－145.

［14］ Nishimura T，Wakamatsu M，Morega A M.Oscillatory double－diffusive convection in arectangular enclosure with combined horizontaltemperature and concentration gradients［J］.International Journal of Heat and Mass Transfer，1998，41（11）：1601－1611.

［15］ Chamkha A J，Al－Naser H.Hydromagnetic doublediffusive convection in arectangular enclosure with opposing temperature and concentration gradients［J］.International Journal of Heat and Mass Transfer，2002，45（12）：2465－2483.

［16］ Ching Y C，?ztop H F，Rahman M M，et al.Finite element simulation of mixed convection heat and mass transfer in a right triangular enclosure［J］.International Communications in Heat and MassTransfer，2012，39（5）：689－696.

［17］ Chamkha A J，Al－Naser H.Double－diffusive convection in an inclined porous enclosure with opposing temperature and concentration gradients［J］.International Journal of Thermal Sciences，2001，40（3）：227－244.

［18］ Sun H，Lauriat G，Sun D L，et al.Transient doublediffusive convection in an enclosure with large density variations［J］.International Journal of Heat and Mass Transfer，2010，53（4）：615－625.

［19］ Hasanuzzaman M，Rahman M M，?ztop H F，et al.Effects of Lewis number on heat and mass transfer in a triangular cavity［J］.International Communications in Heat and Mass Transfer，2012，39（8）：1213－1219.

［20］ Gray D D，Giorgin A.The validity of the boussinesq approximation for liquids and gases［J］.International Journal of Heat and Mass Transfer，1976，19（5）：545－551.

［21］ Patankar S V.Numerical heat transfer and fluid flow［M］.New York：McGraw－Hill，1980.