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

(上海理工大學 能源與動力工程學院, 上海 200093)

0 引 言

當流體中同時存在溫度梯度和濃度梯度時，就會引發(fā)溫度擴散和濃度擴散，并且兩者相互作用，進而發(fā)生自然對流現(xiàn)象。這種由兩者綜合而引起的自然對流被稱為雙擴散對流[1]。雙擴散對流現(xiàn)象是一種常見的自然現(xiàn)象，在海洋學、大氣物理學、化學動力學等領域中也會涉及到這種流動，例如海洋環(huán)流、化學反應和大氣擴散等方面[2-5]。

其中室內火災的排煙是一個典型的雙擴散對流問題。當室內發(fā)生火災，室內溫度急劇上升，并且伴隨大量排煙，室內空氣在溫度差和濃度差的相互作用下運動。國內外學者對火災物理現(xiàn)象的研究大都采用如CFAST[6-7]和FDS[8-9]等火災模擬軟件，以及大渦模擬方法(LES)和雷諾時均法(RANS)[10-11]。

對雙擴散問題國內外眾多學者進行了大量的研究。Al-Amiri等[12]基于Galerkin加權余量法對頂蓋驅動的正方形腔體內的雙擴散混合對流過程進行了研究，得出了在理查德森數(shù)Ri較低時，頂蓋驅動作用會使方腔內傳質和傳熱速率增加的結論。在此基礎之上，Zhuo等[13]利用格子Boltzmann方法模擬分析了頂蓋驅動深腔內的流動分岔，研究了長寬比對第一分岔問題的影響，得出了與時間相關的渦流隨雷諾數(shù)Re呈現(xiàn)周期性和非周期性變化。雍玉梅等[14]利用格子Boltzmann方法，在前人基礎上模擬了溶解的氧氣分子在整個正方形方腔內的擴散過程，同時為控制氣體傳質過程的給熱條件提供依據(jù)。Xu等[15-17]則采用有限元方法研究了內置發(fā)熱圓開口方腔內，不同的特性參數(shù)和發(fā)熱圓位置對含有發(fā)熱圓開口方腔內雙擴散自然對流的影響。

上述研究大都基于傅里葉導熱定律和菲克擴散效應進行分析；然而當流體中溫度和濃度梯度較大時，Soret和Dufour效應則不能忽略。Soret效應是指溫度場的不均勻導致的傳質現(xiàn)象，Dufour效應是指濃度場的不均勻導致的傳熱現(xiàn)象。Soret效應為正值時，表示溫度梯度會驅使傳質從高溫到低溫運動，負值時表示溫度梯度會驅使傳質從低溫到高溫運動。Dufour效應與之類似，不再贅述。

國內外學者對Soret和Dufour效應也做了相應研究。Nithyadevi和Yang[18]在考慮Soret和Dufour效應的情況下，采用有限體積法分析了水在部分加熱方腔內的雙擴散自然對流，并研究了不同加熱位置、瑞利數(shù)、浮升力比等對流場和傳熱傳質的影響。而Cheng[19]采用邊界層近似解法，在考慮Soret和Dufour效應的條件下，研究了在充滿多孔介質的傾斜波浪狀表面的自由對流邊界層，發(fā)現(xiàn)了隨著Soret數(shù)的增大無量綱總傳熱率降低，而隨著Dufour數(shù)的增大無量綱總傳熱率增加的規(guī)律。Wang等[20]則采用了SIMPLE算法，通過構造一個包含Soret和Dufour效應的水平方腔，研究了雙擴散對流在非均勻網格中的振蕩特性，得出了隨著浮升力比和瑞利數(shù)的增加，雙擴散對流會經過穩(wěn)態(tài)對流，周期振蕩，擬周期性震蕩，混沌流動最后再恢復到周期振蕩的發(fā)展過程。近些年來，更多的科研人員開始在磁流體力學(Magnetohydrodynamics，MHD)對流的研究中考慮Soret和Dufour效應[21-23]。

本文主要從室內火災以及污染物擴散現(xiàn)象出發(fā)，抽象出內置高濃度發(fā)熱圓方腔內雙擴散自然對流現(xiàn)象問題，利用格子Boltzmann方法對其進行研究。格子Boltzmann作為一種數(shù)值模擬方法，具有天生的并行特性以及邊界條件處理簡單、程序易于實施、演化過程清晰等特點[24]。本文考慮Soret和Dufour效應對方腔內部傳熱傳質現(xiàn)象的影響，得到了不同參數(shù)條件下方腔內部的流線、等溫線和等濃度線分布特性，以及發(fā)熱圓表面的平均努賽爾數(shù)Nuav和平均舍伍德數(shù)Shav的變化規(guī)律。

1 問題描述和計算模型

1.1 物理問題

本文以室內火災物理現(xiàn)象抽象如下模型：將室內簡化為一個方腔，室內火災點簡化為方腔內發(fā)熱圓。其簡化物理模型如圖1所示。方腔模型四周壁面均為低溫低濃度，方腔內部發(fā)熱圓為高溫高濃度，模型方腔的長寬均為L，方腔中心有一直徑為d(d=0.4L)的發(fā)熱圓，發(fā)熱圓表面的溫度為Th，濃度為Ch，四周壁面溫度為Tc，濃度為Cc，重力加速度為g，其中Th>Tc，Ch>Cc。

圖1 物理模型Fig.1 Physical model

1.2 數(shù)學模型

假設流體為不可壓且不考慮黏性熱耗散，采用Boussinesq假設[25]，浮升力項中的密度由下式給出：

ρ=ρ0[1-βT(T-T0)-βC(C-C0)](1)

其中：ρ0、T0和C0分別為參考密度、溫度和濃度；βT為溫度引起的體積膨脹系數(shù)，K-1；βC為濃度引起的體積膨脹系數(shù)，m3·kg-1。

基于以上假設，描述二維考慮Soret和Dufour效應的方腔內雙擴散自然對流的宏觀控制方程為：

其中：u和p分別為二維速度矢量和壓力；T和C分別為溫度和濃度；ν、α和D分別是運動黏性系數(shù)、熱擴散系數(shù)和質擴散系數(shù)；κTC和κCT分別代表Soret系數(shù)和Dufour系數(shù)；F為浮力項，即:

F=-g[1-βT(T-T0)-βC(C-C0)](6)

定義參考速度uR=α/L，引入下列無量綱變量：

其中：t、p、T和C分別為時間、壓力、溫度和濃度，u、v分別為直角坐標系下水平方向和豎直方向的速度分量。

在上述條件下，圖1所示的方腔內雙擴散自然對流問題在直角坐標系下的無量綱宏觀數(shù)學描述方程如下：

Ra·Pr(θ-Br·S)(10)

描述用到的6個無量綱準則數(shù)是普朗特數(shù)Pr、瑞利數(shù)Ra、浮升力比Br、路易斯數(shù)Le、Soret數(shù)Sr和Dufour數(shù)Df，其定義分別如下：

其中：ν、α、D和g分別表示運動粘度、熱擴散系數(shù)、質擴散系數(shù)和重力加速度；在模擬分析中，瑞利數(shù)Ra設為105，普朗特數(shù)Pr設為0.71，路易斯數(shù)Le設為2.0。

邊界條件設定如下：

發(fā)熱圓表面的無量綱速度均為0，無量綱溫度和濃度均為1,即：U=V=0,θ=1.0,S=1.0。

四周壁面的無量綱速度、溫度和濃度均為0,即：U=V=0,θ=0,S=0。

2 Lattice Boltzmann模型

本文基于Guo等[26]提出的不可壓縮CLBGK模型進行分析，該模型的基本思想是用三個獨立的LBGK方程分別模擬速度場、溫度場和濃度場，然后通過Boussinesq假設近似將其耦合起來。

速度場模擬采用D2Q9模型，其演化方程為：

fi(r+eiΔt,t+Δt)-fi(r,t)=

其中ei為離散速度，采用Qian等[27]提出的離散速度模型，表示為：

其中，ωi代表權重系數(shù)，取值如下所示：

演化方程中外力項Fi表示為：

+βc(C-C0)](19)

相應的流體宏觀速度和壓力的計算式為：

溫度場和濃度場模擬采用D2Q5模型，其演化方程分別為：

Ti(r+eiΔt,t+Δt)-Ti(r,t)=

Ci(r+eiΔt,t+Δt)-Ci(r,t)=

相應的宏觀溫度和濃度的計算式表示為：

熱擴散系數(shù)α和質擴散系數(shù)D表示為：

發(fā)熱圓表面的平均努賽爾數(shù)Nuav和平均舍伍德數(shù)Shav由下式確定：

下式中n表示發(fā)熱圓表面法線方向上的單位矢量，φ表示圓角弧度值。

本模型發(fā)熱圓曲面邊界處理采用Bouzidi等[30]提出的空間插值的邊界處理方法，方腔壁面的邊界處理格式采用Guo等[31]提出的非平衡態(tài)外推格式；網格數(shù)采用200×200，速度、溫度及濃度最大絕對殘差均小于10-7。

3 程序驗證

為了對程序的準確性進行驗證，本文采用該方法對豎直腔內雙擴散自然對流進行了模擬，并將本文計算結果與文獻[32]中所得結果進行對比。表1給出了本文與對比文獻[32]中方腔熱壁面的Nuav和Shav的計算結果和誤差，兩者基本吻合。

表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 模擬結果與分析

圖2為Br=5.0時，方腔內流場、溫度場和濃度場在不同Sr數(shù)和Df數(shù)下的分布圖。從圖中可以看出，方腔內的流線、等溫線和等濃度線均關于y軸對稱，這是因為流體所處的邊界條件以及所受的外力均關于y軸對稱，且在瑞利數(shù)Ra=1×105的條件下，內部流場的流動沒有振蕩等非線性現(xiàn)象的產生。從流線圖的分布可知，在方腔頂部左右兩側各形成一個渦結構。當Sr數(shù)和Df數(shù)繼續(xù)同時減小時，左右兩個渦開始向方腔頂部擴散，并且在方腔底部左右兩側又各自產生了一個二次渦結構，渦的強度和范圍不太明顯。從等溫線可以看出，Sr=Df=-0.4時，發(fā)熱圓上方出現(xiàn)了熱羽流，并且隨著Sr數(shù)和Df數(shù)的增大，發(fā)熱圓上方的熱羽流變得越來越弱，方腔底部兩側的等溫線逐漸向發(fā)熱圓靠近；發(fā)熱圓上方附近的等溫線的溫度梯度較下方小。發(fā)熱圓下方的等溫線則變得稀疏，溫度梯度變小，整體傳熱能力下降。觀察等濃度線的變化可以發(fā)現(xiàn)，隨著Sr數(shù)和Df數(shù)逐漸增大，發(fā)熱圓上方的質羽流越來越弱，而且發(fā)熱圓下方的質邊界層越來越厚，傳質能力逐漸變弱，兩側的質羽流向下方靠攏。由方程(11)和(12)可知，當Sr數(shù)和Df數(shù)減小時，方程(11)和(12)右側數(shù)值的減小，方程左側的數(shù)值亦會減小，即無量綱溫度和濃度隨時間和在x，y方向的偏導數(shù)減小，所以其傳質能力和傳熱能力會均有所減弱。

圖3為Br=-5.0時，方腔內流場、溫度場和濃度場在不同Sr數(shù)和Df數(shù)下的分布圖。從圖3中可以觀察可知，當浮升力比Br=-5.0時，流線、等溫度線和等濃度線的分布圖與Br=5.0時的分布圖基本呈反對稱分布，因為Br并非單純的浮力，而是質浮升力和熱浮升力的比值，當Br=-5.0時，表明質浮升力和熱浮升力方向相反，且質浮升力占主導地位，故兩者基本呈反對稱分布。從流線分布圖可知，在流場中左右兩側對稱各出現(xiàn)了一個渦結構且旋渦中心靠近方腔的底部；隨著Sr數(shù)和Df數(shù)的同時減小，渦結構逐漸向發(fā)熱圓頂部擴散；當Sr數(shù)和Df數(shù)小于-0.2時，除形成對稱渦結構外，還會在發(fā)熱圓的上方出現(xiàn)二次渦，且運動方向與主旋渦的運動方向相反。觀察等溫線和等濃度線分布圖可以看出，隨著Sr數(shù)和Df數(shù)的同時增大，發(fā)熱圓下方的熱羽流形態(tài)逐漸變弱，兩側的等溫線慢慢向方腔底部偏移，傳熱強度變??；發(fā)熱圓頂部的質羽流變化微弱，但底部的質羽流變得愈不明顯，傳質能力變小，但相對于傳熱能力的減弱程度較小。

圖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)可以看出，當浮升力比Br不變時，發(fā)熱圓表面Nuav隨著Sr數(shù)和Df數(shù)的同時增加而幾乎呈線性減小，說明傳熱能力隨著Sr數(shù)和Df的同時增加而減弱。這是由于Sr數(shù)增強，溫度場的不均勻對傳質的影響增強，傳熱能力增強；但Df數(shù)增加，濃度場的不均勻對傳熱的影響也在增強，傳熱能力減弱，且Dufour效應的作用效果較大，即對傳熱能力的減弱能力更大，因此Nuav總體趨勢減小，傳熱能力減弱。但觀察可以發(fā)現(xiàn)，當Br=-1時,發(fā)熱圓表面的Nuav隨著Sr數(shù)和Df數(shù)的增加下降幅度很小，幾乎呈穩(wěn)定趨勢，這是因為當Br=-1時，質浮升力和熱浮升力的大小相同且方向相反，故而相互抵消，此時主要靠自然對流進行換熱，因此隨著Sr數(shù)和Df數(shù)的同時增大，其傳熱能力變化不大。同樣，從圖4(b)中可以觀察到此趨勢。進而觀察圖4(b)可知，當Sr數(shù)和Df數(shù)一定時，隨著浮升力比Br的增加， 發(fā)熱圓表面Nuav值先減小后增大，說明傳熱能力先減小后增大。這是因為當Br=-5時，質浮升力和熱浮升力方向相反，但質浮升力占主導，會抵消一部分熱浮升力作用，隨著Br的逐漸增大，質浮升力和熱浮升力的抵消效果逐漸增大，且在Br=-1處為低值點且Nuav值相差不大，這是質浮升力和熱浮升力完全相互抵消所致，這時，傳熱能力最弱；當Br逐漸增加為正值時，質浮升力和熱浮升力的方向相同，并且相互共同作用，Nuav值又開始逐漸增加，傳熱能力逐漸增強。

圖5是在不同浮升力比Br、Soret數(shù)Sr和Dufour數(shù)Df下，發(fā)熱圓表面Shav的曲線變化規(guī)律。觀察圖5(a)可知，當浮升力比Br不變時，Shav隨著Sr數(shù)和Df數(shù)的同時增加而緩慢減小，因為Sr數(shù)的增加，溫度場的不均勻對傳質的影響增強，傳質能力減弱；但Df數(shù)增加，濃度場的不均勻對傳熱的影響也在增強，傳質能力增強，且Soret效應的作用效果較大，又因為條件中Le>1，傳質強度本身較大，故其傳質能力減小程度較小。但當Br=-1時，Shav隨著Sr數(shù)和Df數(shù)的增加反而有所增加；這是因為當質浮升力和熱浮升力相互抵消后，Le>1表示其傳質強度較大，所以Shav會有所增加。從圖5(b)可以看出，當Sr數(shù)和Df數(shù)一定時，隨著浮升力比Br的增加，發(fā)熱圓表面Shav值先減小后增大，說明傳質能力先減弱后增強，其同樣是因為Br的不同導致質浮升力和熱浮升力方向不同；且同樣在Br=-1時，Shav值最小，傳質能力最弱。

(a) Nuav變化曲線1

(b) Nuav變化曲線2

(a) Shav變化曲線1

(b) Shav變化曲線2

5 結 論

本文采用格子Boltzmann方法，在給定條件Ra=1×105，Pr=0.71，Le=2.0，并且考慮Soret和Dufour效應的情況下，對方腔內雙擴散自然對流現(xiàn)象進行了研究；為了確保方法的準確性，與已有文獻進行了對比驗證。在不同浮升力比Br(-5≤Br≤5)、Soret數(shù)Sr(-0.4≤Sr≤0.4)和Dufour數(shù)Df(-0.4≤Df≤0.4)條件下，本文的數(shù)值模擬得出如下結論：

(1) 在相同的浮升力比Br下，隨著Soret數(shù)Sr和Dufour數(shù)Df數(shù)值同時增大時，方腔內傳熱能力逐漸減弱；在相同Soret數(shù)Sr和Dufour數(shù)Df下，隨著Br數(shù)的增大，傳熱能力先減小后增大。

(2) 在相同的浮升力比Br下，隨著Soret數(shù)Sr和Dufour數(shù)Df數(shù)值同時增大時，方腔內傳熱能力逐漸減弱，但當Br=-1時，傳熱能力有所增加；在相同Soret數(shù)Sr和Dufour數(shù)Df下，隨著Br數(shù)的增大，傳質能力先減小后增大。

本文采用格子Boltzmann方法對內置高濃度發(fā)熱圓的方腔內雙擴散自然對流現(xiàn)象的研究可為室內污染物擴散過程中的熱質耦合作用機理研究提供相應的理論分析支撐。