熱格子Boltzmann方法分析及應用

陳 杰， 錢躍竑

(上海大學上海市應用數學和力學研究所，上海200072)

格子Boltzmann方法(lattice Boltzmann method，LBM)是近20年發展成熟起來的一種數值計算方法.LBM基于氣體動理論，通過分布函數的演化獲得宏觀信息.作為一種簡單且能處理復雜流動問題的有效數值方法［1-2］，LBM具有良好的數值穩定性、天然的并行性、簡單的邊界處理等優點，自出現之日起就被廣泛用于多孔介質流［3］、多相流［4］、反應擴散系統［5］等諸多領域.早期的LBM只應用于等溫流動(或無熱流動)的模擬，但是基于這種方法具備處理復雜問題的能力以及解決傳熱問題的需要，研究者一直在不斷地探索研究熱格子Boltzmann模型，已形成了一些經過數值驗證具有模擬熱流動能力的熱LBM［6-10］，并應用于多孔介質流動與傳熱、燃燒及化學反應流、湍流等問題.本研究簡述了不同熱格子Boltzmann模型的基本理論，并通過數值分析對比了不同熱格子Boltzmann模型的計算結果及數值特性，進而用于多孔介質流動傳熱問題中.

1 等溫LBM基本原理

LBM中除時間、空間被離散之外，無限維的粒子速度空間也都被離散成有限的速度序列.在標準LBM模型中，物理空間被離散成正方形(體)格子，流體粒子在格點x上碰撞并按離散速度E=［e0，e1，…，eq－1］遷移到x+eiδt格點.fi(x，t)定義為t時刻在格點x上速度為ei的粒子密度，滿足如下的格子Boltzmann方程：

式中，cs為格子聲速，Wi為不同速度粒子的權重.本研究在數值模擬中均采用D2Q9模型.

2 熱格子Boltzmann模型

現有的熱格子Boltzmann模型通常可以分為兩大類：第一類是流場溫度場耦合統一求解的模型，如多速格子Boltzmann模型(multi-speed LBM，MSLBM)、熵格子Boltzmann方法(entropic LBM，ELBM);另一類則是對流場與溫度場分別求解，如被動標量格子Boltzmann模型(passive scalar LBM，PSLBM)、雙分布函數(double-distribution-function，DDF)模型，以及其他與傳統計算流體動力學(computational fluid dynamics，CFD)結合的混合方法，如混合熱格子Boltzmann方法(hybrid-thermal LBM，HTLBM).

2.1 多速格子Boltzmann模型(MSLBM)

多速格子Boltzmann模型是等溫LBM模型的直接推廣，其密度、速度、內能等均由速度分布函數的各階速度矩得到.Qian［6］基于等溫LBGK模型，提出了D1Q5，D2Q13，D3Q21，D3Q25熱力學LBGK模型.在這些模型中，除了要滿足等溫模型的守恒條件外，還應滿足能量守恒和平衡態熱通量為0的條件：

平衡態分布函數是Maxwell分布的截斷形式：

式中，Ap，Bp，Dp為待定參數，由滿足的守恒條件確定.平衡態包含了速度的三階項，離散速度也在D2Q9的基礎上在主坐標軸上增加了4個速度.Qian［6］采用此模型對一維激波管、二維 Rayleigh-Benard對流進行了模擬，證明了該模型的有效性.

MSLBM具有良好的物理基礎，宏觀方程絕對耦合，已成功模擬了一些傳熱現象，但只能模擬狹窄的溫度范圍和較小的Ma數，存在穩定性問題，限制了該模型的廣泛應用.

2.2 熵格子Boltzmann方法(ELBM)

熵格子Boltzmann方法考慮了H定理，通過在守恒約束下最小化波爾茲曼H函數求解平衡態分布函數，由此得出的正定的分布函數保證了模型的穩定性和準確性［11］.Prasianakis等［10］將ELBM拓展到熱流動問題的求解中，證實了該方法的有效性，本研究參照此方法.

H函數定義為

Prasianakis等［12］采用在ELBM中加入高階量的補償算法，較大地提高了基于D2Q9標準格子的ELBM可模擬的溫差和Ma數，但是模型實施較為復雜.

2.3 雙分布函數模型

雙分布函數模型，即存在兩個分布函數：密度分布函數和內能(溫度或總能)分布函數，其中密度分布函數用于模擬速度場，而內能(溫度或總能)分布函數則用來模擬溫度場.溫度、內能或總能分布函數均通過不同的方式構造，但其演化都獨立于密度分布函數.

2.3.1 被動標量格子Boltzmann模型(PSLBM)

被動標量格子Boltzmann模型基于如下原理：在忽略壓力做的功和粘性熱耗散的情況下，溫度可以看作是隨流體運動的一個標量，遵循對流擴散方程.由于此方程與組分濃度場的控制方程一樣，于是Shan［7］提出使用兩組分模型模擬單組分熱流動問題：組分1模擬流體的運動;組分2模擬被動的溫度場.平衡態密度函數為

2.3.2 內能雙分布函數模型

內能雙分布函數模型最早由He等［8］提出，其速度場仍用密度分布函數演化模擬，溫度場則由內能分布函數模擬.該模型的基本思想是通過對連續Boltzmann方程進行特殊的離散得到等溫LBM，如果進行同樣的操作，則熱LBM可以由離散內能的演化方程得到.

根據內能的定義ρε=∫(ξ－u)2/2f dξ，引入內能分布函數g(r，ξ，t)=(ξ－u)2f/2，并引入新的碰撞模型，得到內能分布函數滿足的演化方程：

式中，q=(ξ－u)·［?tu+(ξ·)u］.然后對演化方程離散，得到可用于數值計算的離散的分布演化方程，具體的離散過程詳見文獻［8］.

相比于PSLBM，內能DDF的構造更具有物理基礎，并包含了粘性熱耗散和可壓縮功.相比于MSLBM，DDF模型具有更好的數值穩定性，Pr數不受限制，因此被廣泛用于各種近似不可壓流體流動與傳熱問題.

2.4 混合熱格子Boltzmann模型(HTLBM)

HTLBM是指使用 LBM解速度場，使用傳統CFD解溫度場，并通過一定的方式相互影響.這種方法利用了LBM能簡單處理復雜流動問題的優勢以及傳統CFD在傳熱問題上的成熟技術，可以處理一些僅僅使用傳統CFD較難解決的復雜流動傳熱問題.最初，Lallemand等［13］將多速多松弛模型和有限差分法(finite difference method，FDM)相結合，提出了混合模型，速度場用多松弛LBM求解，溫度場采用FDM求解.

本研究采用有限容積法(finite volume method，FVM)與LBM相結合的混合方法，即采用如下的FVM求解能量守恒方程：

式中，S為廣義源項，包括壓力做的功和粘性熱耗散.

速度場與溫度場的耦合通過在LBM中添加溫度相關的外力項以及在FVM中添加廣義源項S來實現.此外，普朗特數、比熱容等熱物性以及隨溫度變化的輸運系數可以實現相應的調節.本研究中FVM與LBM采用同一套網格系統，FVM采用絕對穩定且具有與LBM相同精度的二階迎風格式(second-order upwind scheme，SUS).

PSLBM，DDF以及HTLBM這類模型的一個關鍵之處在于流場與溫度場之間的耦合，其模型往往不滿足氣體完全狀態方程，溫度場對速度場的影響只是通過施加一個外力來實現.如Guo等［9］針對Boussinesq方程組，通過在密度分布函數演化方程中增加一個外力項以實現溫度對流場的影響.Filippova等［14］基于HTLBM研究了小Ma數下高溫燃燒，用溫度場修正密度場以滿足狀態方程.

3 計算結果及分析

為了進一步對比各類模型，本研究采用ELBM，PSLBM，內能DDF模型以及HTLBM，對熱Couette流、封閉方腔自然對流和多孔介質內非等溫流動等問題進行了模擬對比.

3.1 熱Couette流模擬

考慮兩平板間熱Couette流，上平板以速度U向右運動，下板靜止，且上下平板分別保持恒溫Th，Tc，且Th＞Tc.橫截面溫度廓線的解析形式為

式中，H為平板間距離，Pr=ν/χ為普朗特數，χ為熱擴散系數，Ec=U2/［Cp(Th－Tc)］為埃克特數.

熱Couette流中不考慮流體可壓縮性的影響，而粘性耗散效應明顯，因而分別運用ELBM，內能DDF模型和HTLBM對該問題進行了模擬，網格數均為64×64.模擬中Re=UH/ν=20，計算結果如圖1所示.固定Pr=4，Ec分別為1，10和20的無量綱溫度廓線，散點為不同方法的計算值，曲線為解析解公式(10).由圖可見，三種模型都成功模擬了粘性耗散效應，且與解析解吻合得很好.

本工作進一步研究了三種模型的計算效率問題.圖2給出了溫度殘差隨CPU時間的變化曲線，可見ELBM和HTLBM明顯優于內能DDF模型.

3.2 封閉方腔自然對流模擬

封閉方腔尺寸為H(正方形邊長)，左右壁面分別保持恒溫Th，Tc，且Th＞Tc，上下壁面絕熱，四壁面速度均為無滑移邊界.方腔內充滿均質空氣，考慮向下的重力.描述自然對流的無量綱參數Ra數定義為

圖1 熱Couette流溫度廓線Fig.1 Temperature variation of the thermal Couette flow

圖2 熱Couette流溫度殘差變化曲線Fig.2 Temperature residuals variation of the thermal Couette flow

式中，β為熱膨脹系數.物性滿足Boussinesq假設，這里通過施加外力G=－β(T－T0)g實現溫度場對速度場的影響.

在方腔自然對流中，可壓縮效應以及粘性耗散效應可忽略不計.從模型分析可以看出，PSLBM在這種情況下與DDF模型類似，而ELBM邊界實施較為復雜.因此，本研究分別采用不包含粘性耗散效應的PSLBM和HTLBM對該問題進行了模擬，模擬中Pr=0.71，Ra數分別為104，105和106.圖3和圖4分別為HTLBM在不同Ra數下流動穩定后得到的流線、等溫線，與以往的數值及實驗結果一致.由圖3可見，隨著Ra數的增大，方腔中心的近似圓形的渦逐漸變成橢圓形，進而分裂成兩個渦.當Ra= 106時，兩個渦分別向左右壁面移動，在中心出現了第三個渦.由圖4可見，隨著Ra數的增大，豎直的等溫線逐漸變得水平，主導的傳熱機理由導熱變為對流.

為了進一步定量考核，本研究計算了努塞爾數Nu和平均努塞爾數 Numean.表1給出了熱壁面的Numean、最大Nu數Numax及相應位置的yNumax、水平中心線上最大速度vmax及相應的位置x、垂直中心線上最大速度umax以及相應的位置y.HTLBM和PSLBM求解的結果與Barakos等［15］的基準解一致.

同樣，本研究對HTLBM和PSLBM的計算效率進行了對比，圖5所示為兩種方法模擬自然方腔對流Ra=105時，速度殘差隨CPU時間的變化曲線.可以明顯看出，兩種方法中殘差均呈現震蕩下降趨勢，且HTLBM收斂快于PSLBM，HTLBM殘差收斂到10－7以下時的耗時為PSLBM的57%.

圖3 方腔自然對流不同Ra數的流線Fig.3 Predicted streamlines of natural convection

圖4 方腔自然對流不同Ra數的等溫線Fig.4 Predicted temperature profiles of natural convection

表1 數值解與基準解對比Table 1 Comparison of numerical results between thermal models and benchmarks

圖5 方腔自然對流速度殘差變化曲線Fig.5 Velocity residuals variation of the natural convection

3.3 多孔介質非等溫流動模擬

多孔介質內部結構十分復雜，其流動傳熱現象也相當復雜.格子Boltzmann方法在模擬孔隙內的流體運動時可以方便地使用反彈格式處理復雜流場，因此，該方法在孔隙尺度模擬多孔介質內部復雜流動上有明顯的優勢及較高的計算率.對于多孔介質內流動與傳熱的問題，以往使用比較廣泛的是PSLBM和內能DDF模型.本研究將HTLBM用于多孔介質流動與傳熱分析中，并與PSLBM進行了對比.

本研究分析了分形多孔介質中的自然對流，分形結構采用Sierpinski地毯，依次對分形等級N=2和3的Sierpinski情況進行了模擬.無量綱控制參數Pr=0.71，Ra數分別為104，105和106，固體區域溫度保持線性溫度分布.圖6為采用HTLBM計算N= 2分形結構內自然對流得到的流線圖，圖7為相應的等溫線.由圖可見，模擬結果與PSLBM一致，隨Ra數的逐步增大，傳熱機理由導熱主導變化為對流主導.圖8為N=3，Ra=106時的流線圖及等溫線.由圖可見，固體的增多明顯地抑制了對流作用.

同樣對HTLBM在計算效率的問題上和PSLBM進行了對比.圖9為Ra=106時兩種方法模擬N=2分形結構時的速度殘差曲線，此時HTLBM耗時為PSLBM的76%，仍具有優勢.

圖6 多孔介質方腔自然對流流線(N=2)Fig.6 Predicted streamlines of porous cavity(N=2)

圖7 多孔介質方腔自然對流等溫線(N=2)Fig.7 Predicted temperature profiles of porous cavity(N=2)

圖8 多孔介質方腔自然對流流線及等溫線(N=3)Fig.8 Predicted streamlines and temperature profiles of porous cavity(N=3)

4 結論

本研究簡要介紹了幾種熱格子Boltzmann模型(MSLBM，ELBM，PSLBM，內能DDF模型及HTLBM)，并運用不同熱格子模型求解了兩個典型算例以及多孔介質流動傳熱問題，得到如下結論.

圖9 多孔方腔自然速度殘差變化曲線Fig.9 Velocity residuals variation of porous cavity

(1)速度場溫度場耦合求解的模型還需要進一步發展才能被廣泛應用.

(2)相比于PSLBM和DDF模型，HTLBM在保證計算精度的前提下，具有較高的計算效率.

(3)數值模擬驗證了HTLBM在處理多孔介質復雜結構時可行、有效，且比PSLBM的效率高.

［1］ QIANY H，D’HUMIERESD，LALLEMANDP.Lattice BGK models for Navier-Stokes equation ［J］.Europhysics Letters，1992，17(6)：479-484.

［2］ QIANY H，SUCCIS，ORSZAGS A.Recent advances in lattice Boltzmann computing［M］∥ DIETRICH S.Annual reviews of computational physicsⅢ.New Jersey：World Scientific Publishing Company，1995：195-224.

［3］ ZHAOC Y，DAIL N，TANGG H，et al.Numerical study of natural convection in porous media(metals) using lattice Boltzmann method (LBM) ［J］.International Journal of Heat and Fluid Flow，2010，31 (5)：925-934.

［4］ 嚴永華，石自媛，楊帆.液滴撞擊液膜噴濺過程的LBM模擬［J］.上海大學學報：自然科學版，2008，14(4)：399-404.

［5］ 李青，徐旭峰，周美蓮.三維斑圖形成的格子Boltzmann方法模擬［J］.上海大學學報：自然科學版，2007，13(5)：516-518.

［6］ QIANY H.Simulating thermohydrodynamics with lattice BGK models［J］.Journal of Scientific Computing，1993，8(3)：231-242.

［7］ SHANX.Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method［J］.Physical Review E，1997，55(3)：2780-2788.

［8］ HEX，CHENS，DOOLENG D.A novel thermal model for the latticeBoltzmann method in incompressible limit［J］.Journal of Computational Physics，1998，146 (1)：282-300.

［9］ GUOZ，ZHENGC，SHIB，et al.Thermal lattice Boltzmann equation for low Mach number flows：Decoupling model［J］.Physical Review E，2007，75 (3)：036704.

［10］ PRASIANAKISN I，CHIKATAMALAS S，KARLINI V，et al.Entropic lattice Boltzmann method for simulation of thermal flows［J］.Mathematics and Computers in Simulation，2006，72(2)：179-183.

［11］ ANSUMALIS，KARLINI V，OTTINGERH C.Minimal entropic kinetic models for hydrodynamics ［J］.Europhysics Letters，2003，63(6)：798-804.

［12］ PRASIANAKISN I，KARLINI V.Lattice Boltzmann method for simulation of compressible flows on standard lattices［J］.Physical Review E，2008，78(1)：016704.

［13］ LALLEMANDP，LUO L S.Theoryofthelattice Boltzmann method：Acoustic and thermal properties in two and three dimensions［J］.Physical Review E，2003，68(3)：036706.

［14］ FILLIPPOVAO，HANELlD.A novellatticeBGK approach for low Mach number combustion［J］.Journal of Computational Physics，2000，158(2)：139-160.

［15］ BARAKOSG，MITSOULISE，ASSIMACOPOULOSD.Natural convection flow in a square cavity revisited：Laminar and turbulent models with wall functions［J］.International Journal for Numerical Methods in Fluids，1994，18(7)：695-719.