細微顆粒在方腔自然對流中運動特性的格子Boltzmann模擬

李勇

（華中科技大學煤燃燒國家重點實驗室 湖北武漢 430074）

細微顆粒在方腔自然對流中運動特性的格子Boltzmann模擬

李勇

（華中科技大學煤燃燒國家重點實驗室 湖北武漢 430074）

為了了解顆粒物在氣體中的流動規律，為大氣的顆粒物污染提供理論上的指導，本文將不可壓雙分布熱格子模型與基于點源顆粒的拉格朗日跟蹤法相結合，數值研究了封閉方腔自然對流中的細微顆粒物的運動特性。分析了瑞利數和粒徑大小對顆粒群流動特性的影響。結果表明，瑞利數較小時，顆粒在流場中會形成一個準平衡的循環區域，此時瑞利數和粒徑大小都會對顆粒的運動特性產生影響；而當瑞利數增加到一定臨界值時，不會形成類似的準平衡的循環區域，此時粒徑大小對顆粒的運動影響幾乎可以忽略。

格子Boltzmann方法；氣固兩相流；方腔自然對流

1 引言

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

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

2 數值方法

2.1 氣相流動的格子Boltzmann方法

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

其中fi，Ti分別為速度場和溫度場的離散分布函數，對應的平衡態分布函數定義為：

f(ieq)的參數設定為

模型中的作用力項定義如下：

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

各類參數的選取參照文獻[10]中的設置。

2.2 顆粒運動的拉格朗日跟蹤方法

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

其中，坎寧安校正系數C與顆粒的布朗力FB根據文獻[11]確定，熱泳力根據文獻[4]中的半經驗公式計算。

對方程(6)和(7)在一個時間步長內積分，可以得到顆粒速度和位移的計算格式：

3 模擬結果及討論

運用上述LBE-拉格朗日跟蹤法，我們對細微顆粒在二維封

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

首先，我們模擬了Ra=1000、顆粒直徑為1μm時的工況。圖1給出了不同時刻顆粒在方腔中的運動狀況。從圖中可以發現，一部分顆粒在熱泳力等的作用下很快與冷壁面接觸而被吸附。部分顆粒運動到一定的區域內由于力的作用達到平衡而處于類穩定狀態，這部分顆粒會一直在該區域內運動而不會被冷壁面吸附，形成一個準平衡的循環區域。另外，也有一小部分顆粒由于重力的作用會沉積在下平板上。

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

我們現在考慮不同瑞利數條件下顆粒群在封閉方腔內的流動狀況。顆粒直徑仍設定為1μm。圖2給出了瑞利數從100增加到8×105顆粒群在t=50s時的位置分布。可以發現，當Ra較小時，顆粒最終均會形成一個準平衡的區域并在該區域內運動。并且隨著Ra的增大，該準平衡區域越來越分散。當瑞利數為時，顆粒在方腔內隨機分布而不會形成類似的準平衡區域。

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

圖3 dp=1um時不同瑞利數的顆粒保留率

為了更好地理解顆粒群的運動特性，我們統計了方腔中未被壁面吸附的顆粒數占顆粒總數的比率，即顆粒保留率RF(Remaining Fraction)。不同瑞利數Ra下的顆粒保留率RF隨時間的變化曲線如圖3所示。初始階段所有顆粒隨機分布在方腔內，在熱泳力等的作用下，處于前述準平衡區域之外的外部顆粒迅速向冷壁面運動并被吸附，顆粒保留率RF隨之減小；而運動到準平衡區域之內的內部顆粒將保持在該準平衡區域內運動而不會被吸附，RF最終達到一個穩定值。并且Ra越小，RF減小得越快，最終達到的穩定值也越小。這是因為Ra越小，顆粒達到受力平衡的準平衡區域越接近方腔中心，顆粒更加集中，因此更多的顆粒處于準平衡區域外部，它們在熱泳力等的

作用下向冷壁面運動并被吸附。隨著Ra增大，該準平衡區域逐漸分散和模糊，最終消失。對于圖中Ra最大的情況，由于方腔內沒有明顯的準平衡區域存在，顆粒將逐漸運動到壁面上并被吸附，因而顆粒保留率呈單調下降趨勢。

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

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

4 結語

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

[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]劉若雷，楊瑞昌，由長福等.溫度梯度場內可吸入顆粒物運動特性及熱泳沉積[J].化工學報,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方法的原理及應用 [M].北京∶科學出版社,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—），男，漢族，湖北省武漢市人，華中科技大學煤燃燒國家重點實驗室，碩士研究生，主要研究領域為氣固兩相流動的計算數值模擬。