張拴羊, 徐洪濤, 梁天生, 婁 欽, 楊 茉
(上海理工大學(xué) 能源與動(dòng)力工程學(xué)院, 上海 200093)
當(dāng)流體中同時(shí)存在溫度梯度和濃度梯度時(shí),就會(huì)引發(fā)溫度擴(kuò)散和濃度擴(kuò)散,并且兩者相互作用,進(jìn)而發(fā)生自然對(duì)流現(xiàn)象。這種由兩者綜合而引起的自然對(duì)流被稱為雙擴(kuò)散對(duì)流[1]。雙擴(kuò)散對(duì)流現(xiàn)象是一種常見的自然現(xiàn)象,在海洋學(xué)、大氣物理學(xué)、化學(xué)動(dòng)力學(xué)等領(lǐng)域中也會(huì)涉及到這種流動(dòng),例如海洋環(huán)流、化學(xué)反應(yīng)和大氣擴(kuò)散等方面[2-5]。
其中室內(nèi)火災(zāi)的排煙是一個(gè)典型的雙擴(kuò)散對(duì)流問題。當(dāng)室內(nèi)發(fā)生火災(zāi),室內(nèi)溫度急劇上升,并且伴隨大量排煙,室內(nèi)空氣在溫度差和濃度差的相互作用下運(yùn)動(dòng)。國(guó)內(nèi)外學(xué)者對(duì)火災(zāi)物理現(xiàn)象的研究大都采用如CFAST[6-7]和FDS[8-9]等火災(zāi)模擬軟件,以及大渦模擬方法(LES)和雷諾時(shí)均法(RANS)[10-11]。
對(duì)雙擴(kuò)散問題國(guó)內(nèi)外眾多學(xué)者進(jìn)行了大量的研究。Al-Amiri等[12]基于Galerkin加權(quán)余量法對(duì)頂蓋驅(qū)動(dòng)的正方形腔體內(nèi)的雙擴(kuò)散混合對(duì)流過程進(jìn)行了研究,得出了在理查德森數(shù)Ri較低時(shí),頂蓋驅(qū)動(dòng)作用會(huì)使方腔內(nèi)傳質(zhì)和傳熱速率增加的結(jié)論。在此基礎(chǔ)之上,Zhuo等[13]利用格子Boltzmann方法模擬分析了頂蓋驅(qū)動(dòng)深腔內(nèi)的流動(dòng)分岔,研究了長(zhǎng)寬比對(duì)第一分岔問題的影響,得出了與時(shí)間相關(guān)的渦流隨雷諾數(shù)Re呈現(xiàn)周期性和非周期性變化。雍玉梅等[14]利用格子Boltzmann方法,在前人基礎(chǔ)上模擬了溶解的氧氣分子在整個(gè)正方形方腔內(nèi)的擴(kuò)散過程,同時(shí)為控制氣體傳質(zhì)過程的給熱條件提供依據(jù)。Xu等[15-17]則采用有限元方法研究了內(nèi)置發(fā)熱圓開口方腔內(nèi),不同的特性參數(shù)和發(fā)熱圓位置對(duì)含有發(fā)熱圓開口方腔內(nèi)雙擴(kuò)散自然對(duì)流的影響。
上述研究大都基于傅里葉導(dǎo)熱定律和菲克擴(kuò)散效應(yīng)進(jìn)行分析;然而當(dāng)流體中溫度和濃度梯度較大時(shí),Soret和Dufour效應(yīng)則不能忽略。Soret效應(yīng)是指溫度場(chǎng)的不均勻?qū)е碌膫髻|(zhì)現(xiàn)象,Dufour效應(yīng)是指濃度場(chǎng)的不均勻?qū)е碌膫鳠岈F(xiàn)象。Soret效應(yīng)為正值時(shí),表示溫度梯度會(huì)驅(qū)使傳質(zhì)從高溫到低溫運(yùn)動(dòng),負(fù)值時(shí)表示溫度梯度會(huì)驅(qū)使傳質(zhì)從低溫到高溫運(yùn)動(dòng)。Dufour效應(yīng)與之類似,不再贅述。
國(guó)內(nèi)外學(xué)者對(duì)Soret和Dufour效應(yīng)也做了相應(yīng)研究。Nithyadevi和Yang[18]在考慮Soret和Dufour效應(yīng)的情況下,采用有限體積法分析了水在部分加熱方腔內(nèi)的雙擴(kuò)散自然對(duì)流,并研究了不同加熱位置、瑞利數(shù)、浮升力比等對(duì)流場(chǎng)和傳熱傳質(zhì)的影響。而Cheng[19]采用邊界層近似解法,在考慮Soret和Dufour效應(yīng)的條件下,研究了在充滿多孔介質(zhì)的傾斜波浪狀表面的自由對(duì)流邊界層,發(fā)現(xiàn)了隨著Soret數(shù)的增大無量綱總傳熱率降低,而隨著Dufour數(shù)的增大無量綱總傳熱率增加的規(guī)律。Wang等[20]則采用了SIMPLE算法,通過構(gòu)造一個(gè)包含Soret和Dufour效應(yīng)的水平方腔,研究了雙擴(kuò)散對(duì)流在非均勻網(wǎng)格中的振蕩特性,得出了隨著浮升力比和瑞利數(shù)的增加,雙擴(kuò)散對(duì)流會(huì)經(jīng)過穩(wěn)態(tài)對(duì)流,周期振蕩,擬周期性震蕩,混沌流動(dòng)最后再恢復(fù)到周期振蕩的發(fā)展過程。近些年來,更多的科研人員開始在磁流體力學(xué)(Magnetohydrodynamics,MHD)對(duì)流的研究中考慮Soret和Dufour效應(yīng)[21-23]。
本文主要從室內(nèi)火災(zāi)以及污染物擴(kuò)散現(xiàn)象出發(fā),抽象出內(nèi)置高濃度發(fā)熱圓方腔內(nèi)雙擴(kuò)散自然對(duì)流現(xiàn)象問題,利用格子Boltzmann方法對(duì)其進(jìn)行研究。格子Boltzmann作為一種數(shù)值模擬方法,具有天生的并行特性以及邊界條件處理簡(jiǎn)單、程序易于實(shí)施、演化過程清晰等特點(diǎn)[24]。本文考慮Soret和Dufour效應(yīng)對(duì)方腔內(nèi)部傳熱傳質(zhì)現(xiàn)象的影響,得到了不同參數(shù)條件下方腔內(nèi)部的流線、等溫線和等濃度線分布特性,以及發(fā)熱圓表面的平均努賽爾數(shù)Nuav和平均舍伍德數(shù)Shav的變化規(guī)律。
本文以室內(nèi)火災(zāi)物理現(xiàn)象抽象如下模型:將室內(nèi)簡(jiǎn)化為一個(gè)方腔,室內(nèi)火災(zāi)點(diǎn)簡(jiǎn)化為方腔內(nèi)發(fā)熱圓。其簡(jiǎn)化物理模型如圖1所示。方腔模型四周壁面均為低溫低濃度,方腔內(nèi)部發(fā)熱圓為高溫高濃度,模型方腔的長(zhǎng)寬均為L(zhǎng),方腔中心有一直徑為d(d=0.4L)的發(fā)熱圓,發(fā)熱圓表面的溫度為Th,濃度為Ch,四周壁面溫度為Tc,濃度為Cc,重力加速度為g,其中Th>Tc,Ch>Cc。
圖1 物理模型Fig.1 Physical model
假設(shè)流體為不可壓且不考慮黏性熱耗散,采用Boussinesq假設(shè)[25],浮升力項(xiàng)中的密度由下式給出:
ρ=ρ0[1-βT(T-T0)-βC(C-C0)](1)
其中:ρ0、T0和C0分別為參考密度、溫度和濃度;βT為溫度引起的體積膨脹系數(shù),K-1;βC為濃度引起的體積膨脹系數(shù),m3·kg-1。
基于以上假設(shè),描述二維考慮Soret和Dufour效應(yīng)的方腔內(nèi)雙擴(kuò)散自然對(duì)流的宏觀控制方程為:
其中:u和p分別為二維速度矢量和壓力;T和C分別為溫度和濃度;ν、α和D分別是運(yùn)動(dòng)黏性系數(shù)、熱擴(kuò)散系數(shù)和質(zhì)擴(kuò)散系數(shù);κTC和κCT分別代表Soret系數(shù)和Dufour系數(shù);F為浮力項(xiàng),即:
F=-g[1-βT(T-T0)-βC(C-C0)](6)
定義參考速度uR=α/L,引入下列無量綱變量:
其中:t、p、T和C分別為時(shí)間、壓力、溫度和濃度,u、v分別為直角坐標(biāo)系下水平方向和豎直方向的速度分量。
在上述條件下,圖1所示的方腔內(nèi)雙擴(kuò)散自然對(duì)流問題在直角坐標(biāo)系下的無量綱宏觀數(shù)學(xué)描述方程如下:
Ra·Pr(θ-Br·S)(10)
描述用到的6個(gè)無量綱準(zhǔn)則數(shù)是普朗特?cái)?shù)Pr、瑞利數(shù)Ra、浮升力比Br、路易斯數(shù)Le、Soret數(shù)Sr和Dufour數(shù)Df,其定義分別如下:
其中:ν、α、D和g分別表示運(yùn)動(dòng)粘度、熱擴(kuò)散系數(shù)、質(zhì)擴(kuò)散系數(shù)和重力加速度;在模擬分析中,瑞利數(shù)Ra設(shè)為105,普朗特?cái)?shù)Pr設(shè)為0.71,路易斯數(shù)Le設(shè)為2.0。
邊界條件設(shè)定如下:
發(fā)熱圓表面的無量綱速度均為0,無量綱溫度和濃度均為1,即:U=V=0,θ=1.0,S=1.0。
四周壁面的無量綱速度、溫度和濃度均為0,即:U=V=0,θ=0,S=0。
本文基于Guo等[26]提出的不可壓縮CLBGK模型進(jìn)行分析,該模型的基本思想是用三個(gè)獨(dú)立的LBGK方程分別模擬速度場(chǎng)、溫度場(chǎng)和濃度場(chǎng),然后通過Boussinesq假設(shè)近似將其耦合起來。
速度場(chǎng)模擬采用D2Q9模型,其演化方程為:
fi(r+eiΔt,t+Δt)-fi(r,t)=
其中ei為離散速度,采用Qian等[27]提出的離散速度模型,表示為:
其中,ωi代表權(quán)重系數(shù),取值如下所示:
演化方程中外力項(xiàng)Fi表示為:
+βc(C-C0)](19)
相應(yīng)的流體宏觀速度和壓力的計(jì)算式為:
溫度場(chǎng)和濃度場(chǎng)模擬采用D2Q5模型,其演化方程分別為:
Ti(r+eiΔt,t+Δt)-Ti(r,t)=
Ci(r+eiΔt,t+Δt)-Ci(r,t)=
相應(yīng)的宏觀溫度和濃度的計(jì)算式表示為:
熱擴(kuò)散系數(shù)α和質(zhì)擴(kuò)散系數(shù)D表示為:
發(fā)熱圓表面的平均努賽爾數(shù)Nuav和平均舍伍德數(shù)Shav由下式確定:
下式中n表示發(fā)熱圓表面法線方向上的單位矢量,φ表示圓角弧度值。
本模型發(fā)熱圓曲面邊界處理采用Bouzidi等[30]提出的空間插值的邊界處理方法,方腔壁面的邊界處理格式采用Guo等[31]提出的非平衡態(tài)外推格式;網(wǎng)格數(shù)采用200×200,速度、溫度及濃度最大絕對(duì)殘差均小于10-7。
為了對(duì)程序的準(zhǔn)確性進(jìn)行驗(yàn)證,本文采用該方法對(duì)豎直腔內(nèi)雙擴(kuò)散自然對(duì)流進(jìn)行了模擬,并將本文計(jì)算結(jié)果與文獻(xiàn)[32]中所得結(jié)果進(jìn)行對(duì)比。表1給出了本文與對(duì)比文獻(xiàn)[32]中方腔熱壁面的Nuav和Shav的計(jì)算結(jié)果和誤差,兩者基本吻合。
表1 不同Br下的Nuav和Shav的比較 (Ra=1×105,Pr=1.0,Le=2.0,Sr=0.1,Df=0.1)Table 1 Comparisons of Nuav and Shav at different Br (Ra=1×105,Pr=1.0,Le=2.0,Sr=0.1,Df=0.1)
圖2為Br=5.0時(shí),方腔內(nèi)流場(chǎng)、溫度場(chǎng)和濃度場(chǎng)在不同Sr數(shù)和Df數(shù)下的分布圖。從圖中可以看出,方腔內(nèi)的流線、等溫線和等濃度線均關(guān)于y軸對(duì)稱,這是因?yàn)榱黧w所處的邊界條件以及所受的外力均關(guān)于y軸對(duì)稱,且在瑞利數(shù)Ra=1×105的條件下,內(nèi)部流場(chǎng)的流動(dòng)沒有振蕩等非線性現(xiàn)象的產(chǎn)生。從流線圖的分布可知,在方腔頂部左右兩側(cè)各形成一個(gè)渦結(jié)構(gòu)。當(dāng)Sr數(shù)和Df數(shù)繼續(xù)同時(shí)減小時(shí),左右兩個(gè)渦開始向方腔頂部擴(kuò)散,并且在方腔底部左右兩側(cè)又各自產(chǎn)生了一個(gè)二次渦結(jié)構(gòu),渦的強(qiáng)度和范圍不太明顯。從等溫線可以看出,Sr=Df=-0.4時(shí),發(fā)熱圓上方出現(xiàn)了熱羽流,并且隨著Sr數(shù)和Df數(shù)的增大,發(fā)熱圓上方的熱羽流變得越來越弱,方腔底部?jī)蓚?cè)的等溫線逐漸向發(fā)熱圓靠近;發(fā)熱圓上方附近的等溫線的溫度梯度較下方小。發(fā)熱圓下方的等溫線則變得稀疏,溫度梯度變小,整體傳熱能力下降。觀察等濃度線的變化可以發(fā)現(xiàn),隨著Sr數(shù)和Df數(shù)逐漸增大,發(fā)熱圓上方的質(zhì)羽流越來越弱,而且發(fā)熱圓下方的質(zhì)邊界層越來越厚,傳質(zhì)能力逐漸變?nèi)酰瑑蓚?cè)的質(zhì)羽流向下方靠攏。由方程(11)和(12)可知,當(dāng)Sr數(shù)和Df數(shù)減小時(shí),方程(11)和(12)右側(cè)數(shù)值的減小,方程左側(cè)的數(shù)值亦會(huì)減小,即無量綱溫度和濃度隨時(shí)間和在x,y方向的偏導(dǎo)數(shù)減小,所以其傳質(zhì)能力和傳熱能力會(huì)均有所減弱。
圖3為Br=-5.0時(shí),方腔內(nèi)流場(chǎng)、溫度場(chǎng)和濃度場(chǎng)在不同Sr數(shù)和Df數(shù)下的分布圖。從圖3中可以觀察可知,當(dāng)浮升力比Br=-5.0時(shí),流線、等溫度線和等濃度線的分布圖與Br=5.0時(shí)的分布圖基本呈反對(duì)稱分布,因?yàn)锽r并非單純的浮力,而是質(zhì)浮升力和熱浮升力的比值,當(dāng)Br=-5.0時(shí),表明質(zhì)浮升力和熱浮升力方向相反,且質(zhì)浮升力占主導(dǎo)地位,故兩者基本呈反對(duì)稱分布。從流線分布圖可知,在流場(chǎng)中左右兩側(cè)對(duì)稱各出現(xiàn)了一個(gè)渦結(jié)構(gòu)且旋渦中心靠近方腔的底部;隨著Sr數(shù)和Df數(shù)的同時(shí)減小,渦結(jié)構(gòu)逐漸向發(fā)熱圓頂部擴(kuò)散;當(dāng)Sr數(shù)和Df數(shù)小于-0.2時(shí),除形成對(duì)稱渦結(jié)構(gòu)外,還會(huì)在發(fā)熱圓的上方出現(xiàn)二次渦,且運(yùn)動(dòng)方向與主旋渦的運(yùn)動(dòng)方向相反。觀察等溫線和等濃度線分布圖可以看出,隨著Sr數(shù)和Df數(shù)的同時(shí)增大,發(fā)熱圓下方的熱羽流形態(tài)逐漸變?nèi)?,兩?cè)的等溫線慢慢向方腔底部偏移,傳熱強(qiáng)度變小;發(fā)熱圓頂部的質(zhì)羽流變化微弱,但底部的質(zhì)羽流變得愈不明顯,傳質(zhì)能力變小,但相對(duì)于傳熱能力的減弱程度較小。
圖2 不同Sr和Df數(shù)下的流線、等溫線和等濃度線分布(Ra=1×105,Pr=0.71,Le=2.0,Br=5.0)Fig.2 Distributions of streamlines, isotherms and isoconcentrations at different Sr and Df(Ra=1×105,Pr=0.71,Le=2.0,Br=5.0)
圖3 不同Sr和Df數(shù)下的流線、等溫線和等濃度線分布(Ra=1×105,Pr=0.71,Le=2.0,Br=-5.0)Fig.3 Distributions of streamlines, isotherms and isoconcentrations at different Sr and Df(Ra=1×105,Pr=0.71,Le=2.0,Br=-5.0)
圖4是在不同浮升力比Br、Soret數(shù)Sr和Dufour數(shù)Df下,發(fā)熱圓表面Nuav的曲線變化規(guī)律。從圖4(a)可以看出,當(dāng)浮升力比Br不變時(shí),發(fā)熱圓表面Nuav隨著Sr數(shù)和Df數(shù)的同時(shí)增加而幾乎呈線性減小,說明傳熱能力隨著Sr數(shù)和Df的同時(shí)增加而減弱。這是由于Sr數(shù)增強(qiáng),溫度場(chǎng)的不均勻?qū)髻|(zhì)的影響增強(qiáng),傳熱能力增強(qiáng);但Df數(shù)增加,濃度場(chǎng)的不均勻?qū)鳠岬挠绊懸苍谠鰪?qiáng),傳熱能力減弱,且Dufour效應(yīng)的作用效果較大,即對(duì)傳熱能力的減弱能力更大,因此Nuav總體趨勢(shì)減小,傳熱能力減弱。但觀察可以發(fā)現(xiàn),當(dāng)Br=-1時(shí),發(fā)熱圓表面的Nuav隨著Sr數(shù)和Df數(shù)的增加下降幅度很小,幾乎呈穩(wěn)定趨勢(shì),這是因?yàn)楫?dāng)Br=-1時(shí),質(zhì)浮升力和熱浮升力的大小相同且方向相反,故而相互抵消,此時(shí)主要靠自然對(duì)流進(jìn)行換熱,因此隨著Sr數(shù)和Df數(shù)的同時(shí)增大,其傳熱能力變化不大。同樣,從圖4(b)中可以觀察到此趨勢(shì)。進(jìn)而觀察圖4(b)可知,當(dāng)Sr數(shù)和Df數(shù)一定時(shí),隨著浮升力比Br的增加, 發(fā)熱圓表面Nuav值先減小后增大,說明傳熱能力先減小后增大。這是因?yàn)楫?dāng)Br=-5時(shí),質(zhì)浮升力和熱浮升力方向相反,但質(zhì)浮升力占主導(dǎo),會(huì)抵消一部分熱浮升力作用,隨著Br的逐漸增大,質(zhì)浮升力和熱浮升力的抵消效果逐漸增大,且在Br=-1處為低值點(diǎn)且Nuav值相差不大,這是質(zhì)浮升力和熱浮升力完全相互抵消所致,這時(shí),傳熱能力最弱;當(dāng)Br逐漸增加為正值時(shí),質(zhì)浮升力和熱浮升力的方向相同,并且相互共同作用,Nuav值又開始逐漸增加,傳熱能力逐漸增強(qiáng)。
圖5是在不同浮升力比Br、Soret數(shù)Sr和Dufour數(shù)Df下,發(fā)熱圓表面Shav的曲線變化規(guī)律。觀察圖5(a)可知,當(dāng)浮升力比Br不變時(shí),Shav隨著Sr數(shù)和Df數(shù)的同時(shí)增加而緩慢減小,因?yàn)镾r數(shù)的增加,溫度場(chǎng)的不均勻?qū)髻|(zhì)的影響增強(qiáng),傳質(zhì)能力減弱;但Df數(shù)增加,濃度場(chǎng)的不均勻?qū)鳠岬挠绊懸苍谠鰪?qiáng),傳質(zhì)能力增強(qiáng),且Soret效應(yīng)的作用效果較大,又因?yàn)闂l件中Le>1,傳質(zhì)強(qiáng)度本身較大,故其傳質(zhì)能力減小程度較小。但當(dāng)Br=-1時(shí),Shav隨著Sr數(shù)和Df數(shù)的增加反而有所增加;這是因?yàn)楫?dāng)質(zhì)浮升力和熱浮升力相互抵消后,Le>1表示其傳質(zhì)強(qiáng)度較大,所以Shav會(huì)有所增加。從圖5(b)可以看出,當(dāng)Sr數(shù)和Df數(shù)一定時(shí),隨著浮升力比Br的增加,發(fā)熱圓表面Shav值先減小后增大,說明傳質(zhì)能力先減弱后增強(qiáng),其同樣是因?yàn)锽r的不同導(dǎo)致質(zhì)浮升力和熱浮升力方向不同;且同樣在Br=-1時(shí),Shav值最小,傳質(zhì)能力最弱。
(a) Nuav變化曲線1
(b) Nuav變化曲線2
(a) Shav變化曲線1
(b) Shav變化曲線2
本文采用格子Boltzmann方法,在給定條件Ra=1×105,Pr=0.71,Le=2.0,并且考慮Soret和Dufour效應(yīng)的情況下,對(duì)方腔內(nèi)雙擴(kuò)散自然對(duì)流現(xiàn)象進(jìn)行了研究;為了確保方法的準(zhǔn)確性,與已有文獻(xiàn)進(jìn)行了對(duì)比驗(yàn)證。在不同浮升力比Br(-5≤Br≤5)、Soret數(shù)Sr(-0.4≤Sr≤0.4)和Dufour數(shù)Df(-0.4≤Df≤0.4)條件下,本文的數(shù)值模擬得出如下結(jié)論:
(1) 在相同的浮升力比Br下,隨著Soret數(shù)Sr和Dufour數(shù)Df數(shù)值同時(shí)增大時(shí),方腔內(nèi)傳熱能力逐漸減弱;在相同Soret數(shù)Sr和Dufour數(shù)Df下,隨著Br數(shù)的增大,傳熱能力先減小后增大。
(2) 在相同的浮升力比Br下,隨著Soret數(shù)Sr和Dufour數(shù)Df數(shù)值同時(shí)增大時(shí),方腔內(nèi)傳熱能力逐漸減弱,但當(dāng)Br=-1時(shí),傳熱能力有所增加;在相同Soret數(shù)Sr和Dufour數(shù)Df下,隨著Br數(shù)的增大,傳質(zhì)能力先減小后增大。
本文采用格子Boltzmann方法對(duì)內(nèi)置高濃度發(fā)熱圓的方腔內(nèi)雙擴(kuò)散自然對(duì)流現(xiàn)象的研究可為室內(nèi)污染物擴(kuò)散過程中的熱質(zhì)耦合作用機(jī)理研究提供相應(yīng)的理論分析支撐。