細(xì)微顆粒在方腔自然對(duì)流中運(yùn)動(dòng)特性的格子Boltzmann模擬

李勇

（華中科技大學(xué)煤燃燒國家重點(diǎn)實(shí)驗(yàn)室 湖北武漢 430074）

細(xì)微顆粒在方腔自然對(duì)流中運(yùn)動(dòng)特性的格子Boltzmann模擬

李勇

（華中科技大學(xué)煤燃燒國家重點(diǎn)實(shí)驗(yàn)室 湖北武漢 430074）

為了了解顆粒物在氣體中的流動(dòng)規(guī)律，為大氣的顆粒物污染提供理論上的指導(dǎo)，本文將不可壓雙分布熱格子模型與基于點(diǎn)源顆粒的拉格朗日跟蹤法相結(jié)合，數(shù)值研究了封閉方腔自然對(duì)流中的細(xì)微顆粒物的運(yùn)動(dòng)特性。分析了瑞利數(shù)和粒徑大小對(duì)顆粒群流動(dòng)特性的影響。結(jié)果表明，瑞利數(shù)較小時(shí)，顆粒在流場(chǎng)中會(huì)形成一個(gè)準(zhǔn)平衡的循環(huán)區(qū)域，此時(shí)瑞利數(shù)和粒徑大小都會(huì)對(duì)顆粒的運(yùn)動(dòng)特性產(chǎn)生影響；而當(dāng)瑞利數(shù)增加到一定臨界值時(shí)，不會(huì)形成類似的準(zhǔn)平衡的循環(huán)區(qū)域，此時(shí)粒徑大小對(duì)顆粒的運(yùn)動(dòng)影響幾乎可以忽略。

格子Boltzmann方法；氣固兩相流；方腔自然對(duì)流

1 引言

鑒于近年來頻繁出現(xiàn)的霧霾天氣，對(duì)于細(xì)微顆粒物對(duì)大氣的污染治理及控制越來越緊迫，因而研究顆粒物在流體中的運(yùn)動(dòng)規(guī)律具有重大意義。顆粒的熱泳沉積被視作一種有效的除塵手段，一些學(xué)者對(duì)這類問題進(jìn)行過研究。1998年，Thakurta等人應(yīng)用DNS對(duì)細(xì)微顆粒在管道流中的熱泳沉積進(jìn)行了數(shù)值研究[1]。Tsai等人2001年發(fā)展了一個(gè)簡(jiǎn)單的關(guān)系式來評(píng)估熱泳作用對(duì)氣溶膠顆粒沉積作用的影響[2]。2009年劉若雷等人研究了可吸入顆粒物在溫度梯度場(chǎng)內(nèi)的運(yùn)動(dòng)特性，得到了PM2.5熱泳沉積的半經(jīng)驗(yàn)公式[3]。同年，Akbar等人對(duì)細(xì)微顆粒在方腔自然對(duì)流中的運(yùn)動(dòng)進(jìn)行了模擬[4]，結(jié)果表明熱泳力和布朗擴(kuò)散效應(yīng)對(duì)顆粒的運(yùn)動(dòng)起主要作用。然而，他們的結(jié)果并沒有指出當(dāng)瑞利數(shù)和粒徑大小變化時(shí)這兩種作用力對(duì)顆粒運(yùn)動(dòng)特性的影響。因此，本文將對(duì)顆粒群在方腔自然對(duì)流中的運(yùn)動(dòng)進(jìn)行數(shù)值模擬，并研究瑞利數(shù)和粒徑大小對(duì)顆粒運(yùn)動(dòng)特性的影響。

隨著大型計(jì)算機(jī)性能的不斷提高，數(shù)值模擬方法被迅速應(yīng)用到氣固兩相流領(lǐng)域的研究中[5-8]。相對(duì)于傳統(tǒng)的數(shù)值計(jì)算方法而言，格子Boltzmann方法（Lattice Boltzmann Method，LBM）具有演化過程清晰、編程容易和適合并行計(jì)算等優(yōu)勢(shì)。因此LBM被廣泛應(yīng)用氣固兩相流的數(shù)值模擬[9]。郭照立等人2002年提出的不可壓雙分布熱格子模型（DDF-LBE）具有計(jì)算簡(jiǎn)單、數(shù)值穩(wěn)定性良好等優(yōu)點(diǎn)[10]。因此本文使用上述熱格子LBM模型，并與拉格朗日跟蹤法結(jié)合，模擬顆粒群在方腔自然對(duì)流中的流動(dòng)，研究瑞利數(shù)與粒徑大小對(duì)顆粒群運(yùn)動(dòng)特性的影響，希望能為細(xì)微顆粒的污染提供一些理論上的指導(dǎo)。

2 數(shù)值方法

2.1 氣相流動(dòng)的格子Boltzmann方法

本文考慮二維問題，采用郭照立等人2002年提出的不可壓雙分布（DDF-LBE）熱格子模型模擬流場(chǎng)。該模型用兩組演化方程分別模擬速度場(chǎng)和溫度場(chǎng)，其演化方程如下：

其中fi，Ti分別為速度場(chǎng)和溫度場(chǎng)的離散分布函數(shù)，對(duì)應(yīng)的平衡態(tài)分布函數(shù)定義為：

f(ieq)的參數(shù)設(shè)定為

模型中的作用力項(xiàng)定義如下：

流體的宏觀速度、壓力、溫度分別計(jì)算如下：

各類參數(shù)的選取參照文獻(xiàn)[10]中的設(shè)置。

2.2 顆粒運(yùn)動(dòng)的拉格朗日跟蹤方法

對(duì)于固體顆粒的運(yùn)動(dòng)模擬，采用基于點(diǎn)源顆粒的拉格朗日跟蹤方法，主要考慮重力，曳力，布朗力和熱泳力對(duì)顆粒的作用，并且忽略顆粒對(duì)氣相的作用以及顆粒間的相互作用。在上述假設(shè)下，顆粒的運(yùn)動(dòng)方程可以表示為：

其中，坎寧安校正系數(shù)C與顆粒的布朗力FB根據(jù)文獻(xiàn)[11]確定，熱泳力根據(jù)文獻(xiàn)[4]中的半經(jīng)驗(yàn)公式計(jì)算。

對(duì)方程(6)和(7)在一個(gè)時(shí)間步長(zhǎng)內(nèi)積分，可以得到顆粒速度和位移的計(jì)算格式：

3 模擬結(jié)果及討論

運(yùn)用上述LBE-拉格朗日跟蹤法，我們對(duì)細(xì)微顆粒在二維封

閉方腔自然對(duì)流中的流動(dòng)進(jìn)行了模擬。模擬時(shí)所采用網(wǎng)格大小是256×256，封閉方腔的上下平板均是絕熱壁面，左平板為高溫壁面，右平板為低溫壁面。方腔自然對(duì)流的瑞利數(shù)其中L是特征長(zhǎng)度，v和D分別為流體的運(yùn)動(dòng)粘度和熱擴(kuò)散率，β為熱體積膨脹系數(shù)。起始狀態(tài)顆粒以初速度為零隨機(jī)分散在方腔內(nèi)，顆粒數(shù)目為2000。計(jì)算過程中，假設(shè)顆粒和壁面發(fā)生接觸即被壁面吸附。

首先，我們模擬了Ra=1000、顆粒直徑為1μm時(shí)的工況。圖1給出了不同時(shí)刻顆粒在方腔中的運(yùn)動(dòng)狀況。從圖中可以發(fā)現(xiàn)，一部分顆粒在熱泳力等的作用下很快與冷壁面接觸而被吸附。部分顆粒運(yùn)動(dòng)到一定的區(qū)域內(nèi)由于力的作用達(dá)到平衡而處于類穩(wěn)定狀態(tài)，這部分顆粒會(huì)一直在該區(qū)域內(nèi)運(yùn)動(dòng)而不會(huì)被冷壁面吸附，形成一個(gè)準(zhǔn)平衡的循環(huán)區(qū)域。另外，也有一小部分顆粒由于重力的作用會(huì)沉積在下平板上。

圖1 Ra=1000時(shí)1um顆粒在不同時(shí)刻的位置分布: (a)0s,(b)5s,(c)25s,(d)50s

我們現(xiàn)在考慮不同瑞利數(shù)條件下顆粒群在封閉方腔內(nèi)的流動(dòng)狀況。顆粒直徑仍設(shè)定為1μm。圖2給出了瑞利數(shù)從100增加到8×105顆粒群在t=50s時(shí)的位置分布?？梢园l(fā)現(xiàn)，當(dāng)Ra較小時(shí)，顆粒最終均會(huì)形成一個(gè)準(zhǔn)平衡的區(qū)域并在該區(qū)域內(nèi)運(yùn)動(dòng)。并且隨著Ra的增大，該準(zhǔn)平衡區(qū)域越來越分散。當(dāng)瑞利數(shù)為時(shí)，顆粒在方腔內(nèi)隨機(jī)分布而不會(huì)形成類似的準(zhǔn)平衡區(qū)域。

圖2 不同瑞利數(shù)Ra下顆粒群的位置分布: (a)100,(b)103,(c)104,(d)8×105

圖3 dp=1um時(shí)不同瑞利數(shù)的顆粒保留率

為了更好地理解顆粒群的運(yùn)動(dòng)特性，我們統(tǒng)計(jì)了方腔中未被壁面吸附的顆粒數(shù)占顆?？倲?shù)的比率，即顆粒保留率RF(Remaining Fraction)。不同瑞利數(shù)Ra下的顆粒保留率RF隨時(shí)間的變化曲線如圖3所示。初始階段所有顆粒隨機(jī)分布在方腔內(nèi)，在熱泳力等的作用下，處于前述準(zhǔn)平衡區(qū)域之外的外部顆粒迅速向冷壁面運(yùn)動(dòng)并被吸附，顆粒保留率RF隨之減??；而運(yùn)動(dòng)到準(zhǔn)平衡區(qū)域之內(nèi)的內(nèi)部顆粒將保持在該準(zhǔn)平衡區(qū)域內(nèi)運(yùn)動(dòng)而不會(huì)被吸附，RF最終達(dá)到一個(gè)穩(wěn)定值。并且Ra越小，RF減小得越快，最終達(dá)到的穩(wěn)定值也越小。這是因?yàn)镽a越小，顆粒達(dá)到受力平衡的準(zhǔn)平衡區(qū)域越接近方腔中心，顆粒更加集中，因此更多的顆粒處于準(zhǔn)平衡區(qū)域外部，它們?cè)跓嵊玖Φ鹊?/p>

作用下向冷壁面運(yùn)動(dòng)并被吸附。隨著Ra增大，該準(zhǔn)平衡區(qū)域逐漸分散和模糊，最終消失。對(duì)于圖中Ra最大的情況，由于方腔內(nèi)沒有明顯的準(zhǔn)平衡區(qū)域存在，顆粒將逐漸運(yùn)動(dòng)到壁面上并被吸附，因而顆粒保留率呈單調(diào)下降趨勢(shì)。

同時(shí)，考慮粒徑對(duì)流動(dòng)過程中顆粒保留率RF變化曲線的影響。我們統(tǒng)計(jì)了粒徑分別為1μm,200nm,50nm時(shí)不同Ra下的顆粒保留率RF隨時(shí)間變化，如圖4所示。由圖4中可知，Ra較小時(shí)，粒徑大小對(duì)顆粒保留率的影響比較明顯。粒徑越小，流動(dòng)過程中顆粒保留率RF減小得越快，并且最終達(dá)到的穩(wěn)定值也越小。這是因?yàn)閷?duì)于細(xì)微顆粒，布朗運(yùn)動(dòng)的劇烈程度較大，布朗擴(kuò)散的加劇使顆粒的沉積速度增加。而隨著瑞利數(shù)增加，熱泳力的作用越來越強(qiáng)烈，粒徑大小的影響則逐漸減弱；當(dāng)瑞利數(shù)增加到足夠大時(shí)，熱泳力對(duì)顆粒的影響起主導(dǎo)作用，粒徑大小的影響幾乎可以忽略，如圖4(d)。

圖4 顆粒直徑對(duì)顆粒保留率的影響: Ra=(a)100,(b)103,(c)104,(d)8×105

4 結(jié)語

為了研究顆粒物在氣體中的運(yùn)動(dòng)規(guī)律，為顆粒污染物的治理提供一些理論指導(dǎo)。本文將不可壓雙分布（DDF-LBE）熱格子模型與拉格朗日跟蹤法相結(jié)合，對(duì)顆粒群在方腔自然對(duì)流里的流動(dòng)進(jìn)行了數(shù)值研究。結(jié)果表明，當(dāng)瑞利數(shù)較小時(shí)，顆粒在流場(chǎng)中會(huì)形成一個(gè)準(zhǔn)平衡的循環(huán)區(qū)域，在這個(gè)區(qū)域內(nèi)的顆粒保持準(zhǔn)平衡狀態(tài)，不會(huì)被壁面吸附。隨著瑞利數(shù)的增大，該準(zhǔn)平衡區(qū)域逐漸分散、模糊，最終消失，此時(shí)顆粒逐漸向壁面運(yùn)動(dòng)并被吸附。當(dāng)瑞利數(shù)較小時(shí)，粒徑大小對(duì)顆粒群在流體中的流動(dòng)特性影響較為明顯。粒徑越小，在整個(gè)流動(dòng)時(shí)間段里的顆粒保留率減小得越快，最后達(dá)到的穩(wěn)定值也越小；隨著瑞利數(shù)增大，粒徑大小的影響逐漸減弱；當(dāng)瑞利數(shù)逐漸增大到某一個(gè)臨界值時(shí)，粒徑大小對(duì)顆粒運(yùn)動(dòng)特性的影響基本可以忽略。

[1] Thakurta D.G., Chen M., Mclaughlin J.B. etc. Thermophoretic depositionof small particles in a direct numerical simulation of turbulentchannel flow [J]. International Journal of Heat and Mass Transfer,1998,41,4167- 4182.

[2] Tsai R., Liang L.J. Correlation for thermophoretic deposition ofaerosol particles onto cold plates [J].Journal of Aerosol Science,2001,32,473- 487.

[3]劉若雷，楊瑞昌，由長(zhǎng)福等.溫度梯度場(chǎng)內(nèi)可吸入顆粒物運(yùn)動(dòng)特性及熱泳沉積[J].化工學(xué)報(bào),2009,60,1623-1628.

[4] M.K.Akbar, M.Rahman,S.M. Ghiaasiaan. Particle transport in asmall square enclosure in laminar natural convection [J]. Aerosol Science,2009,40,747- 761.

[5] U. Lantermann, D. Haanel. Particle Monte Carlo and lattice- Boltzmannmethods for simulations of gas- particle flows [J]. Comput. Fluids,2007, 36∶ 407- 422.

[6] I. Alexander. Effects of the variation of mass loading and particledensity in gas- solid particle flow in pipes [J].Powder Technol,2009, 193∶176-181.

[7] Balachandar S, Eaton JK. Turbulent dispersed multiphase flow[J].Annual Reviewof Fluid Mechanics,2010,42∶ 111- 133.

[8] Durst F, Milojevic D, Schonung B. Eulerian and Lagrangian predictionsof particulate two- phase flows∶ A numerical study [J].AppliedMathematical Modelling, 1984,8(2)∶101- 115.

[9]郭照立,鄭楚光.格子Boltzmann方法的原理及應(yīng)用 [M].北京∶科學(xué)出版社,2009.

[10] Z.L.Guo, B.C.Shi, and C.G.Zheng. A coupled lattice BGK modelfor the Boussinesq equations[J]. Int.T.Numer.Fluids,2002,39, 325- 342.

[11] Friedlander, S. K. Smoke dust and haze. 2nd ed [M]. Oxford UniversityPress,2000.

[12] M.Honrrmann and M.Peric. Finite volume multigrid prediction oflaminar natural convection∶Benchmark solutions [J]. Int. J. Numer.Method.Fluids,1990,11,189- 207.

李勇（1989—），男，漢族，湖北省武漢市人，華中科技大學(xué)煤燃燒國家重點(diǎn)實(shí)驗(yàn)室，碩士研究生，主要研究領(lǐng)域?yàn)闅夤虄上嗔鲃?dòng)的計(jì)算數(shù)值模擬。