Soret和Dufour效應(yīng)對(duì)方腔內(nèi)雙擴(kuò)散自然對(duì)流影響的格子Boltzmann模擬

張拴羊, 徐洪濤, 梁天生, 婁 欽, 楊 茉

(上海理工大學(xué) 能源與動(dòng)力工程學(xué)院, 上海 200093)

0 引 言

當(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ī)律。

1 問題描述和計(jì)算模型

1.1 物理問題

本文以室內(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

1.2 數(shù)學(xué)模型

假設(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。

2 Lattice Boltzmann模型

本文基于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。

3 程序驗(yàn)證

為了對(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)

4 模擬結(jié)果與分析

圖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

5 結(jié) 論

本文采用格子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)的理論分析支撐。