格子-Boltzmann方法模擬霜結(jié)晶生長

龔建英 孫金絹 李國君

(西安交通大學(xué)能源與動力工程學(xué)院熱流科學(xué)與工程教育部重點實驗室 西安 710049)

格子-Boltzmann方法模擬霜結(jié)晶生長

龔建英 孫金絹 李國君

(西安交通大學(xué)能源與動力工程學(xué)院熱流科學(xué)與工程教育部重點實驗室 西安 710049)

為了研究冷壁面上霜結(jié)晶的生長過程以及各參數(shù)的變化規(guī)律，利用格子-Boltzmann方法(LBM)方法，建立一個二維介觀模型，實驗驗證了模型的可靠性，分析了霜層溫度變化和密度增長規(guī)律，并直觀模擬出霜結(jié)晶凝聚變化過程。結(jié)果表明:霜層表面溫度在早期階段迅速增加，但增加率隨著結(jié)霜時間增加而減小;霜層內(nèi)部溫度隨霜層厚度的增加呈線性增長;霜層平均密度隨結(jié)霜時間增加呈現(xiàn)出先慢后快的規(guī)律;隨著結(jié)霜過程的進行，由于越往上，霜晶體積分?jǐn)?shù)越小，導(dǎo)致霜層內(nèi)部密度隨霜層厚度的上升而減小。

格子-Boltzmann方法 冷壁面 霜結(jié)晶

1 引言

結(jié)霜現(xiàn)象的研究是涉及制冷、空調(diào)、航空、航天，農(nóng)業(yè)等多領(lǐng)域中具有重要意義的一項基礎(chǔ)性研究。霜的存在不僅降低系統(tǒng)的傳熱性能，增加壓力損失，而且會導(dǎo)致除霜的能耗和運行成本的增加。因此，研究霜結(jié)晶生長的物理機制一直受到國內(nèi)外學(xué)者的普遍關(guān)注。

近幾十年中，結(jié)霜現(xiàn)象的研究方法主要分為數(shù)學(xué)模擬和實驗方法。Lee等［1］把霜層看作多孔介質(zhì)，提出了一個不同的理論模型。Sahin 發(fā)現(xiàn)霜的密度不僅僅是影響霜層有效導(dǎo)熱系數(shù)的唯一參數(shù)。許多學(xué)者還通過實驗方法研究了霜的形成過程，如，Tao等［3］用顯微鏡拍攝結(jié)霜過程，并把結(jié)霜過程分為3個階段:結(jié)晶生長期，霜層生長期和霜層充分生長期。Kandula［4］用實驗方法分析了冷壁面層流下霜的生長。雖然已經(jīng)對霜層生長作了許多研究，霜結(jié)晶凝聚的物理機制和影響過程尚未完全明確，有必要采用更加先進有效的數(shù)學(xué)模擬方法結(jié)合實驗，揭示霜結(jié)晶凝聚本質(zhì)規(guī)律。

隨著數(shù)值計算的快速發(fā)展，格子-Boltzmann方法(LBM)是近些年出現(xiàn)的一種新穎的數(shù)值計算方法，該方法基于分子動理論，具有清晰的物理背景，在宏觀上是離散方法，微觀上是連續(xù)方法，因而被稱為介觀模擬方法。該方法近些年來已成功地模擬流體流動和復(fù)雜結(jié)構(gòu)的一個新工具［5］。與傳統(tǒng)模擬方法相比，由于格子-Boltzmann方法對于復(fù)雜結(jié)構(gòu)的適應(yīng)性受到了廣泛的關(guān)注，該方法在許多傳統(tǒng)計算方法難以勝任的領(lǐng)域，如多孔介質(zhì)、多相流等中取得了眾多成功的應(yīng)用［6］。

本文采用LBM對冷壁面霜晶生長過程進行數(shù)值模擬，并對模型的可靠性進行驗證。從介觀層面揭示了霜層表面溫度和密度的變化規(guī)律，并闡明機理，直觀模擬出霜晶體生長的二維變化圖像，并對其進行深入分析。

2 格子-Boltzmann方法

格子-Boltzmann方法經(jīng)過多年發(fā)展，很多模型被提出用以滿足不同的應(yīng)用要求。其中最基本最常用的模型，大致分為3部分:(1)等溫模型;(2)非等溫模型;(3)其它非標(biāo)準(zhǔn)模型。其中，本文中流場和溫度場均采用了等溫模型Bhatnagar-Gross-Krook(BGK)［7］近似法中的D2Q9模型。流體流動的格子-Boltzmann方程表達如下:

式中:fi是粒子速度分布函數(shù)，ei是離散速度，δt是時間步長，τ是無量綱松弛時間，與運動粘度v=cs2(τ-0.5)δt有關(guān)。平衡態(tài)分布函數(shù)定義為:

式中:聲速cs=δx RT/δt，δx 是網(wǎng)格步長，RT=1/3。權(quán)系數(shù)分別為 ω0=4/9;當(dāng) i=1、2、3、4 時，ωi=1/9;當(dāng) i=5、6、7、8 時，ωi=1/36。離散速度 ei可以表達為:

因此，模型的宏觀密度和速度可以用以下方程計算:

通過Chapman-Enskog展開方法，格子-Boltzmann方程可以推導(dǎo)出Navier-Stokes下的連續(xù)方程和動量方程［8］。

溫度場中的演化方程如下所示:

式中:Ti是溫度分布態(tài)函數(shù)，τc是無量綱松弛時間。溫度的平衡態(tài)函數(shù)T(eq)i表達如下:

權(quán)系數(shù)ωi和等溫聲速cs與速度場中介紹的一樣。模型的溫度T定義如下:T=∑i

Ti(8)

通過Chapman-Enskog展開方法，格子-Boltzmann方程可以推導(dǎo)出本模型下對應(yīng)的Navier-Stokes溫度宏觀方程［9］。

格子-Boltzmann方程的演化過程比較清晰，主要包括遷移和碰撞兩個過程，遷移和碰撞是兩個相對獨立的粒子運動過程。

粒子遷移:

粒子碰撞:

3 邊界條件

LBM中，邊界條件通過分布函數(shù)獲得?？諝鈦砹魉俣缺３譃楹愣ǖ某跏妓俣?，空氣上方、出流的格點速度與冷壁面格點速度均為0，空氣上方、來流和出流的格點溫度保持恒定濕空氣溫度，而冷壁面格點溫度保持與冷板溫度一致。對于冷壁面格點和固體格點均采用無滑移邊界條件中的反彈格式［7］，即當(dāng)粒子到達格點后被原路返回到流體內(nèi)部，且運動方向和入射方向相反。符合準(zhǔn)靜態(tài)假設(shè)，雖然邊界只提供了一階精度，但是該方法需要較少的計算時間，能有效地處理復(fù)雜結(jié)構(gòu)的無滑移條件。如圖1a所示，邊界格點(i，1)的未知分布函數(shù)可由式(13)獲得，即:

除了以上格點采用反彈格式，其余格點均采用了非平衡態(tài)外推格式［7］，這種格式數(shù)值穩(wěn)定性好，計算量小并且適用范圍廣泛。邊界格點上的分布函數(shù)主要分為平衡態(tài)和非平衡態(tài)兩部分，其中平衡態(tài)部分是由一個新構(gòu)造的平衡態(tài)分布來近似，非平衡態(tài)部分是由一個一階精度的外推方法來獲得。如圖1b所示，格點EOA位于邊界上，格點BCD位于流場內(nèi)，格點FGH位于流場外。邊界格點O的未知分布函數(shù)f2，5，6(xO，t)，可由式(14)獲得，即:

圖1 邊界格式Fig.1 Boundary scheme

4 模擬結(jié)果和討論

運用LB模型模擬介觀尺度下冷壁面上霜結(jié)晶生長過程及其各參數(shù)變化規(guī)律。采用的幾何結(jié)構(gòu)被劃分為800×400的二維網(wǎng)格。首先，給定各個物理變量的初始值，如空氣的濕度、速度、溫度、密度，冷壁面溫度，松弛時間，晶格大小，以及邊界條件等。在網(wǎng)格底部，分布有一定數(shù)量的隨機晶核。晶核的數(shù)量由霜晶體初始體積分?jǐn)?shù)決定。根據(jù)方程(2)和方程(7)，求出流場和溫度場初始的平衡態(tài)分布函數(shù)。初始的粒子速度分布函數(shù)為fi(x，t0)=feqi(x，t0)，溫度分布函數(shù)為Ti(x，t0)=Teqi(x，t0)。然后，通過成核理論［10］求得霜晶體質(zhì)量增加，并得到霜晶體新一步的體積分?jǐn)?shù)。而格點隨著霜晶體體積分?jǐn)?shù)增加，由流體相變?yōu)楣腆w相。最后，由式(1)—式(8)結(jié)合邊界條件，求得流場和溫度場的分布函數(shù)和濕空氣在每一個格點上的密度、速度和溫度的分布。

圖2是霜層表面溫度隨霜結(jié)晶生長的變化曲線。圓圈代表Lee的實驗數(shù)據(jù)，線代表模型計算出的結(jié)果。從圖中可以看出，模擬的結(jié)果與文獻中的實驗數(shù)據(jù)的變化趨勢是比較一致的，證明了模型的可靠性。

圖2 霜表面溫度隨時間的變化Fig.2 Effect of frost time on frost surface temperature

霜層表面溫度是一個重要的參數(shù)，對通過水分子凝聚在冷壁面造成的霜晶體的密度變化及形態(tài)影響很大，由圖2可知，在霜晶體生長初期階段，霜表面溫度變化迅速，但是隨著時間的推移，溫度的變化速率變得緩慢，并不斷趨近于三相點溫度。

圖3是在矩形域內(nèi)霜結(jié)晶生長二維幾何結(jié)構(gòu)的溫度分布圖。圖中縱坐標(biāo)代表的霜晶體縱向生長的高度，橫坐標(biāo)表示的是霜晶體的橫向生長。從圖中可以直觀地觀察到晶體的凝聚與溫度隨結(jié)霜時間的變化過程，底部靠近冷壁面的位置，霜晶溫度較低，越往上溫度越高，而由于周圍環(huán)境中溫度較高的濕空氣向內(nèi)擴散，冰晶空隙處溫度較高。圖4顯示了不同時刻霜層內(nèi)溫度隨厚度的變化曲線，從圖中可以看出在任一結(jié)霜時刻，霜層內(nèi)部溫度隨霜層厚度增大幾乎都呈線性增長。在相同霜層厚度處，隨時間的推移溫度不斷降低，但變化幅度很小。

圖3 霜結(jié)晶凝聚直觀圖Fig.3 Simulated temperature field at different times

圖4 霜內(nèi)部溫度不同時刻隨霜層厚度的變化Fig.4 Effect of frost thickness on frost internal temperatures

圖5給出了霜的平均密度隨結(jié)霜時間的變化曲線。由圖可見:霜的密度隨時間在不斷增加。在結(jié)晶初期，由于霜晶呈現(xiàn)較強的枝狀生長，表現(xiàn)在霜層厚度的增加，而密度增長較慢，所以初期曲線呈凹形。隨著霜結(jié)晶的進行，濕空氣分子不斷擴散和凝華，使密度后期不斷增大，所以后期曲線斜率較大，即增長速率較大。除此之外，霜內(nèi)部密度還隨著霜層的厚度而變化，如圖6所示。從圖中可以看出，同一時刻，隨著霜層厚度的增加，霜層密度在降低，這是由于越往上，霜晶體積分?jǐn)?shù)越小，而密度越接近空氣密度;在相同霜層同一厚度處，隨著時間的增加，密度越來越大，原因在于水蒸氣分子不斷沉積凝聚到霜晶體中，導(dǎo)致密度不斷增大。

圖5 霜的平均密度隨時間的變化Fig.5 Effect of frost time on frost average density

5 結(jié)論

圖6 霜內(nèi)部密度不同時刻隨霜層厚度的變化Fig.6 Effect of frost internal density onfrost layer thickness

建立了與傳統(tǒng)方法不同的格子-Boltzmann模型成功模擬霜表面溫度和密度隨霜結(jié)晶生長的變化過程。模擬結(jié)果與文獻中的實驗數(shù)據(jù)較為吻合，得出以下主要結(jié)論:

(1)霜層表面溫度在早期階段迅速增加，但增加率減小。此外，內(nèi)部溫度隨霜層厚度的上升呈線性曲線。

(2)直觀模擬出晶體的凝聚與溫度隨結(jié)霜時間的變化過程，霜層靠近冷壁面的位置，霜晶溫度較低，沿冷壁面法線方向越往上霜層溫度越高。

(3)霜層平均密度隨結(jié)霜時間增加呈現(xiàn)出先慢后快的規(guī)律;霜層內(nèi)部密度隨霜層厚度的上升而減小，原因在于隨著結(jié)霜過程的進行，越往上，霜晶體積分?jǐn)?shù)越小，導(dǎo)致密度減小。

1 Lee K S，Kim W S，Lee T H.A one-dimensional model for frost formation on a cold flat surface［J］.International Journal of Heat and Mass Transfer，1997，40(18):4359-4365.

2 Sahin A Z.Effective thermal conductivity of frost during the crystal growth period［J］.International Journal of Heat and Mass Transfer，2000，43(4):539-553.

3 Tao Y X，Besant R W，Rezkallah K S.A mathematical model for predicting the densification and growth of frost on a flat plate［J］.International Journal of Heat and Mass Transfer，1993，36(2):356-363.

4 Kandula M.Frost growth and densification on a flat surface in laminar flow with variable humidity［J］.International Communications in Heat and Mass Transfer，2012，39(8):1030-1034.

5 Pan Chongxun，Luo Li-Shi，Miller C T.An evaluation of lattice Boltzmann schemes for porous medium flow simulation［J］.Computers ＆Fluids，2006，35(8-9):898-909.

6 Guo Zhaoli，Zhao T S.A lattice Boltzmann model for convection heat transfer in porous media［J］.Numerical Heat Transfer，Part B:Fundamentals:An International Journal of Computation and Methodology，2005，47(2):157-177.

7 郭照立，鄭楚光.格子Boltzmann方法的原理及應(yīng)用［M］.北京:科學(xué)出版社，2009.

8 Hou Shuling，Zou Qisu，Chen Shiyi，et al.Simulation of cavity flow by the lattice Boltzmann method［J］.Journal of Computational Physics，1995，118(2):329-347.

9 Ergun S.Fluid flow through packed columns［J］.Chemical Engineering Progress，1952，48(2):89-94.

10 Becker R，Doring.W The kinetic treatment of nuclear formation in supersaturated vapors［J］.Annals of Physics-Leipzig，1935，24:719.

Frost crystal growth simulation by lattice Boltzmann method

Gong Jianying Sun Jinjuan Li Guojun

(MOE Key Laboratory of Thermo-Fluid Science and Engineering，School of Energy and Power Engineering，Xi’an Jiaotong University，Xi’an 710049，China)

In order to calculate the frost crystal growth process on a cold flat surface and study the variation laws of parameters，a two-dimensional mesoscopic model was proposed based on LBM.The predicted results were found to be in good agreement with the data reported in the literatures.This simulation was concerned with temperature and density variation of the frost.The results show that the surface temperature increases rapidly in the early-stage，but the increase rate decreases with frosting.In addition，the interior temperatures rise in a linear curve with increase of the frost layer thickness.The average frost layer density increases with increasing frosting time.It was also found that the nearer the frost layer gets to the cold flat surfaces the higher the frost density is.

lattice Boltzmann method;frost crystal growth;cold flat surfaces

2013-03-07;

2013-06-14

國家自然科學(xué)基金項目(51106013、50806059)、中國博士后科學(xué)基金(2012M511998)。

龔建英，女，36歲，博士、講師。

TB611、TK121

A文章編號:1000-6516(2013)04-0010-04