Soret和Dufour效應對方腔內雙擴散自然對流影響的格子Boltzmann模擬

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

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

0 引 言

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

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

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

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

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

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

1 問題描述和計算模型

1.1 物理問題

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

圖1 物理模型Fig.1 Physical model

1.2 數學模型

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

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

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

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

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

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

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

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

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

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

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

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

邊界條件設定如下：

發熱圓表面的無量綱速度均為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代表權重系數，取值如下所示：

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

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

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

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

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

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

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

熱擴散系數α和質擴散系數D表示為：

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

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

本模型發熱圓曲面邊界處理采用Bouzidi等[30]提出的空間插值的邊界處理方法，方腔壁面的邊界處理格式采用Guo等[31]提出的非平衡態外推格式；網格數采用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數和Df數下的分布圖。從圖中可以看出，方腔內的流線、等溫線和等濃度線均關于y軸對稱，這是因為流體所處的邊界條件以及所受的外力均關于y軸對稱，且在瑞利數Ra=1×105的條件下，內部流場的流動沒有振蕩等非線性現象的產生。從流線圖的分布可知，在方腔頂部左右兩側各形成一個渦結構。當Sr數和Df數繼續同時減小時，左右兩個渦開始向方腔頂部擴散，并且在方腔底部左右兩側又各自產生了一個二次渦結構，渦的強度和范圍不太明顯。從等溫線可以看出，Sr=Df=-0.4時，發熱圓上方出現了熱羽流，并且隨著Sr數和Df數的增大，發熱圓上方的熱羽流變得越來越弱，方腔底部兩側的等溫線逐漸向發熱圓靠近；發熱圓上方附近的等溫線的溫度梯度較下方小。發熱圓下方的等溫線則變得稀疏，溫度梯度變小，整體傳熱能力下降。觀察等濃度線的變化可以發現，隨著Sr數和Df數逐漸增大，發熱圓上方的質羽流越來越弱，而且發熱圓下方的質邊界層越來越厚，傳質能力逐漸變弱，兩側的質羽流向下方靠攏。由方程(11)和(12)可知，當Sr數和Df數減小時，方程(11)和(12)右側數值的減小，方程左側的數值亦會減小，即無量綱溫度和濃度隨時間和在x，y方向的偏導數減小，所以其傳質能力和傳熱能力會均有所減弱。

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

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

圖5是在不同浮升力比Br、Soret數Sr和Dufour數Df下，發熱圓表面Shav的曲線變化規律。觀察圖5(a)可知，當浮升力比Br不變時，Shav隨著Sr數和Df數的同時增加而緩慢減小，因為Sr數的增加，溫度場的不均勻對傳質的影響增強，傳質能力減弱；但Df數增加，濃度場的不均勻對傳熱的影響也在增強，傳質能力增強，且Soret效應的作用效果較大，又因為條件中Le>1，傳質強度本身較大，故其傳質能力減小程度較小。但當Br=-1時，Shav隨著Sr數和Df數的增加反而有所增加；這是因為當質浮升力和熱浮升力相互抵消后，Le>1表示其傳質強度較大，所以Shav會有所增加。從圖5(b)可以看出，當Sr數和Df數一定時，隨著浮升力比Br的增加，發熱圓表面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效應的情況下，對方腔內雙擴散自然對流現象進行了研究；為了確保方法的準確性，與已有文獻進行了對比驗證。在不同浮升力比Br(-5≤Br≤5)、Soret數Sr(-0.4≤Sr≤0.4)和Dufour數Df(-0.4≤Df≤0.4)條件下，本文的數值模擬得出如下結論：

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

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

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