時間步長對格子玻爾茲曼法模擬室內(nèi)氣流精度的影響

韓夢濤

（華中科技大學建筑與城市規(guī)劃學院，湖北 武漢 430074）

近年，基于格子玻爾茲曼法的大渦模擬（lattice Boltzmann method-based large-eddy simulation，LBM-LES）開始應用于建筑［1-3］和城市風環(huán)境［4］模擬。與當前風環(huán)境主流的有限體積法（finite volume method，F(xiàn)VM）在宏觀尺度上求解物理量不同，LBM用虛擬的、包含有限種速度模式的微觀分布函數(shù)表示流體粒子的集合，并通過分布函數(shù)的碰撞和遷移來模擬流體運動［5］。與基于FVM 的大渦模擬（FVM-LES，即風環(huán)境模擬的主流LES 方法）［6-7］相比，LBM-LES 算法簡單，邊界條件易于實現(xiàn)，且在LES 計算中無需求解壓力泊松方程［8］，計算速度更快，在復雜湍流風環(huán)境模擬中具有較大潛力［9-10］。

LBM的控制方程是格子玻爾茲曼方程，其中的關鍵項是碰撞算子。碰撞算子的形式?jīng)Q定了待求解流體的性質(zhì)。BGK（Bhatnagar-Gross-Krook）近似模型［11］是最常用的碰撞算子，能以簡單的形式獲得足夠精度的模擬結果。既往研究表明，BGK算子在模擬環(huán)境風問題，尤其是室內(nèi)氣流中取得了較好的成果。Elhadidi等［1］利用含BGK算子的LBM模擬了較粗網(wǎng)格條件下室內(nèi)空間的氣流分布并與FVM 進行了比較。Han 等［2］則系統(tǒng)總結了使用LBM-LES 模擬室內(nèi)空氣流動時不同計算條件對模擬精度的影響，并與FVM-LES 進行了詳細比較。上述研究表明含BGK算子的LBM-LES可有效模擬室內(nèi)空氣流動，并可取得與FVM-LES類似的模擬結果。

同時理論分析表明，利用含BGK 算子的LBMLES 可推導出低馬赫數(shù)（Mach number，M）流體的連續(xù)性方程和納維-斯托克斯（Navier-Stokes，N-S）方程，但推導出的N-S 方程與風環(huán)境模擬中常用的方程形式有所差異［12-14］，從而產(chǎn)生模擬誤差。該誤差與模擬時間步長δt的設定緊密相關，不恰當?shù)摩膖設置可能導致計算結果出現(xiàn)壓縮性誤差或數(shù)值振蕩。目前，在應用LBM-LES模擬風環(huán)境問題時，δt的設置引起的誤差問題尚未引起足夠重視和充分討論。為了在應用LBM-LES 模擬風環(huán)境時對如何設置δt提供依據(jù)，本文梳理了既往研究，從理論方面系統(tǒng)總結了δt的設置可能引起的壓縮性誤差或數(shù)值振蕩，并以室內(nèi)空氣流動案例為例，定量討論了δt造成壓縮性誤差或數(shù)值振蕩的程度。

1 含BGK的LBM中與時間步長相關的誤差理論分析與總結

1.1 含BGK的格子玻爾茲曼方程及LBM-LES 方法的理論回顧

本節(jié)首先對含BGK 碰撞算子的格子玻爾茲曼方程及LBM-LES 方法的理論進行簡要回顧，以便后文對誤差進行理論分析。該方程如式（1）所示。

式中：fa為a方向分布函數(shù)；ea為a方向上fa的離散速度；r和t分別為fa所在位置向量和時間為fa的平衡函數(shù)；δt為離散時間步長；τ為fa的松弛時間。

格子玻爾茲曼法方程在微觀層面描述了流體粒子的分布函數(shù)隨時間發(fā)展的演化過程。當分布函數(shù)確定后，流體速度u、密度ρ及壓力p等宏觀物理量可通過式（2）求得。其中，es為格子聲速，在三維問題中的值為，其他參數(shù)含義同前。

在面對高雷諾數(shù)Re 的湍流問題時，可基于LBM 開展LES 計算（LBM-LES）。根據(jù)LES 理論，流體的總粘性νtot由分子粘性ν及亞格子粘性νsgs共同構成（即νtot=ν+νsgs）［15］。同時，基于LBM 理論，流體的總粘性νtot與總松弛時間τtot存在如式（3）所示的關系：

基于式（3）可用總松弛時間τtot替換式（1）中的松弛時間τ以開展LBM-LES 計算。與傳統(tǒng)FVMLES 相同，只需采取合適的LES 亞網(wǎng)格模型計算亞網(wǎng)格粘性νsgs，即可開展LBM-LES計算。

1.2 實際物理量到格子玻爾茲曼單位的轉換

在式（1）中，以真實物理量（通常包含量綱）所度量的流體原型問題首先被映射到格子玻爾茲曼單位的物理量（通常是量綱一的）進行模擬，模擬完成后再將其映射回真實物理量以輸出結果。故模擬第一步是確定適當?shù)霓D換參數(shù)。在無外力等溫流體問題中，LBM主要關注流體粘性ν、速度u、壓力p及位置r參數(shù)，在進行上述物理量的轉換時只需網(wǎng)格分辨率δx和時間步長δt兩個轉換參數(shù)。各物理量轉換關系如式（4）所示，其中上標ph 和lb 分別表示真實物理量及其對應的格子波爾茲曼單位物理量。

其中，網(wǎng)格分辨率δx通常根據(jù)湍流復雜度、模擬所需精度及計算量共同確定，與傳統(tǒng)FVM-LES 基本相同。需要注意的是，LBM中采用的是均勻正立方體網(wǎng)格，無法如FVM-LES 一樣在局部復雜湍流區(qū)域加密網(wǎng)格，故應注意使用測試網(wǎng)格獨立性等方式來確定網(wǎng)格尺寸。這在既往研究中已得到多次驗證［2，4］。而時間步長δt的確定則與FVM-LES有較大差異。過大或過小的δt都可能導致明顯的精度誤差，這將是本文接下來的討論重點。

1.3 壓縮性誤差

利用BGK算子，可從格子玻爾茲曼方程中推導出形如式（5）的N-S方程［14-15］：

式（5）中：O(?2)和O(M3)分別為與?2和M3相關的高階省略項；?為努森數(shù)（Knudsen number），與馬赫數(shù)M成正相關，即O(?2)～O(M2)［12］。式（5）相比標準的N-S 方程多了與M相關的附加項。其中M是以格子波爾茲曼單位定義的，如式（6）所示：

式中：|uph|為局部物理速度uph的大小。

式（5）同時表明，推導出的N-S方程以可壓縮形式存在（即無法消去各項中的密度ρ），意味著LBMLES在處理不可壓縮流體問題時本質(zhì)上是一種偽可壓縮方法，并會產(chǎn)生所謂的“壓縮性誤差”［12］。雖然嚴格來說不可壓的流體并不存在，但在處理諸如建筑與城市風環(huán)境等低流速問題時，F(xiàn)VM-LES 通常采用不可壓縮的N-S 方程。由此可見LBM-LES 與FVM-LES 在計算時針對是否可壓縮的處理方法上有較大區(qū)別。既往研究［12，15］表明，LBM 的壓縮性誤差包括與密度梯度相關的誤差，以及上述與M相關附加項導致的誤差。

式（6）表明，即使在局部風速與網(wǎng)格分辨率δx一定時，較大的時間步長δt會增大M，從而增大使式（5）中與M及?相關項的值（即壓縮性誤差）。故在模擬不可壓縮湍流問題時，與傳統(tǒng)FVM-LES相比，不恰當?shù)摩膖可能導致M增大，從而使得模擬結果產(chǎn)生壓縮性誤差。Skordos［16］曾嘗試用LBM模擬層流狀態(tài)下的二維泰勒渦旋流和剪切流，發(fā)現(xiàn)渦旋流的模擬值和解析值之間的誤差隨著M的減小而減小，并最終實現(xiàn)穩(wěn)定收斂。而在剪切流中，隨著M的減小，模擬誤差呈先減小而后增大的趨勢。Reider 等人［12］從理論上推導了壓縮性誤差，并證實在Re=100且周期為2π的衰減泰勒渦流模擬結果準確性隨著M的降低而提高。

應當注意，BGK 算子中的壓縮性誤差與FVMLES中庫朗數(shù)C的不正確設置引起的誤差并非同一概念，盡管C也是由δx和δt間的取值關系造成。在FVM-LES中，在處理低M不可壓縮流體時，通常建議選擇適當?shù)臅r間步長δt以將C控制在小于1（即C= ||uph＜1），否則將造成結果誤差甚至模擬發(fā)散。而在LBM-LES 中，保證模擬穩(wěn)定性的一個必要條件是M＜1，即|ulb|＜1 3 ≈0.577，否則模擬將直接發(fā)散。故若要保證LBM-LES 模擬正常穩(wěn)定進行，則必有C= ||ulbδlbt δlbx＜1 始終成立（因為LBM規(guī)定了δlbt=1和δlbx=1）。

1.4 超松弛與數(shù)值振蕩

從式（1）中可明顯看出BGK算子的本質(zhì)是表現(xiàn)分布函數(shù)fa(r，t)以一定的速率向平衡分布狀態(tài)feqa(r，t)的演化過程，即松弛過程。該式可改寫為如式（7）的形式：

應當注意，τ不可小于因為根據(jù)式（3），流體的粘性不可為負。Krüger 等人［17］研究了初始條件為f0feq0=1.1、且feq0為恒定值條件下的BGK 算子，得到了如圖1 所示的f0與feq0的關系，對應了上述的三種形態(tài)模式。

圖1 BGK 算子中的亞松弛、全松弛及超松弛算例（重繪自Krüger等［17］）Fig.1 Example of under, full, and over relaxation in BGK collision(reproduced from Krüger et al[17])

圖1 表明理想的碰撞過程是亞松弛或全松弛，即fa平滑地或直接向feqa演化。在實際模擬中，全松弛難以達到，因為不可能經(jīng)過一步就完成模擬。而在超松弛中fa將圍繞feqa振動并以指數(shù)幅度衰減，最后達到feqa。但τ過小可能導致振動過于劇烈，fa無法達到feqa，最終得到錯誤結果。

綜合式（3）、式（4）可得到如式（8）的關系：

式（8）表明，由于流體粘性νphtot通常為定值，故當確定網(wǎng)格分辨率δx后，τtot的大小與δt正相關。δt的設定影響τtot的大小，從而決定了碰撞過程的松弛模式。理想狀態(tài)下τ應不小于1，據(jù)式（8）可知需要或然而這在風環(huán)境模擬中較難滿足。由于空氣動粘性系數(shù)極?。〝?shù)量級為10-5m2·s-1），即使在LES 計算中加上亞網(wǎng)格粘性νsgs也很難大于例如，在風環(huán)境模擬中，由于庫朗數(shù)和計算量的雙重限制，很難在設置δx=1 m 的同時使得δt=10-5s；也無法使用δx=10-3m 來匹配δt=10-1s。盡管在風速較小的局部區(qū)域可能滿足但建筑或城市尺度的風環(huán)境中，大部分區(qū)域（特別是所關注的區(qū)域）主要以高雷諾數(shù)湍流為主，故風環(huán)境模擬中BGK 算子主要表現(xiàn)為超松弛模式。通常這種超松弛模式下的振動類似湍流脈動，然而不恰當?shù)木W(wǎng)格分辨率δx和時間步長δt之間的關系會導致類似圖1中τ=0.51 所示的劇烈振動，最終導致模擬結果產(chǎn)生數(shù)值振蕩。尤其是在相同的網(wǎng)格分辨率δx下，過小的δt會導致νlb過小，從而導致嚴重的數(shù)值振蕩。這與傳統(tǒng)的FVM-LES具有極大的不同。

綜合上述分析可知，在LBM-LES 進行風環(huán)境模擬時，當δx確定后，過大的δt可能導致與M相關的項產(chǎn)生較大的壓縮性誤差，而δt過小則使松弛時間τ減小，可能造成松弛碰撞算子產(chǎn)生數(shù)值振蕩。這是LBM-LES 相比傳統(tǒng)FVM-LES 在模擬設置上的一個重要區(qū)別。FVM-LES 中，只要離散時間步長滿足C＜1 則不會對模擬結果產(chǎn)生顯著的影響。既往研究已經(jīng)討論了層流狀態(tài)下二維流動中LBMLES 的壓縮性誤差［12］，而以風環(huán)境模擬為代表的湍流狀態(tài)下不同時間步長對模擬結果的影響尚未充分討論。下節(jié)將以室內(nèi)氣流為例定量討論這一問題。該案例邊界條件相對簡單、純粹、易于控制，便于進行定量研究。

2 等溫室內(nèi)氣流模擬案例

本文采用國際能源機構推薦的標準等溫室內(nèi)氣流案例（IEA-Annex 20［18］），對含BGK 算子的LBMLES進行壓縮性誤差研究。該案例形狀及取樣位置如圖2 所示，房間特征參數(shù)為L H=3，h H=0.056，r H=0.16。其中L，W和H分別代表房間長度、進深和高度，且H=3.0 m。h和r分別代表氣流入口及出口高度。由房間高度及氣流入口速度定義的Re≈89 000。模擬參數(shù)詳表1。本文中，含標準Smagorinsky 亞 網(wǎng) 格 模 型［19］及BGK 碰 撞 算 子 的LBM-LES 應用于本模擬。Han 等［2］的既往研究已經(jīng)表明該亞網(wǎng)格模型和BGK 算子能較好地適應室內(nèi)氣流模擬，取得滿意的模擬精度。

表1 模擬參數(shù)及相關邊界條件Tab.1 Simulation parameters and boundary conditions

圖2 等溫室內(nèi)氣流案例的幾何形狀及取樣點分布（修改自IEA-Annex 20［18］）Fig.2 Geometry and sampling points distribution of isothermal indoor airflow case (modified from IEA-Annex 20[18])

之前Han 等［2］已經(jīng)討論了LBM-LES 的網(wǎng)格分辨率對模擬精度的影響，本文基于其研究結果僅選取δx=1 75H及1 150H兩種網(wǎng)格，討論不同時間步長δt對模擬精度的影響。具體工況設置如表2所示。工況名稱規(guī)定為“XaTb”，表示網(wǎng)格分辨率及時間步長分別為δx=H a及δt=1bs。

表2 工況設置Tab.2 Case settings

3 模擬結果與討論

3.1 平均與湍流脈動風速結果及精度分析

圖3 顯示了uˉ（時間平均風速的x方向分量）和基于時間平均的脈動風速標準差的x方向分量）的模擬結果。所有結果均用入口風速Uin進行量綱歸一化。圖中添加了Nielsen 等［21］的實驗數(shù)據(jù)用于驗證模擬的準確性。

圖3 部分區(qū)域與的量綱一化模擬結果與實驗結果對比Fig.3 Comparison of simulation and experiment results of normalized andin some regions

除精度最低的X75T50，幾乎所有工況模擬結果均能再現(xiàn)uˉ和的分布趨勢，且模擬精度有隨著δt的減小而提高的趨勢。X75T800和X75T1600中，和均出現(xiàn)了輕微的空間數(shù)值振蕩。而在X150工況組中，X150T200 的和都達到了最佳精度。隨著δt的減小，精度并未提高，反而出現(xiàn)了明顯的數(shù)值振蕩，從而降低模擬精度。

本文采用式（9）所定義的L2誤差范數(shù)［17］εq定量評估模擬精度。式中qEXP(r)和qLBM(r)分別代表實驗和模擬中位置r處的物理量q。L2 誤差范數(shù)考慮了Nielson 實驗數(shù)據(jù)的所有點。εq越小代表模擬與實驗之間的誤差越小，從而模擬精度越高。

圖4顯示了不同δt時所有工況的L2誤差范數(shù)變化曲線。隨著δt從1/50 s 減小到1/200 s，X75 工況組中和的誤差均減小，隨后則顯著增大。可以推測，δt從1/50 s 減小到1/200 s 時的精度提升可能是由于壓縮性誤差的減小；而δt＜1/200 s時精度降低應該是由于超松弛導致，因為在這些工況中觀察到了和的數(shù)值振蕩。對于X150 工況組，δt=1/200 s 時，模擬精度最高，而后隨著δt的減小的模擬精度迅速衰減，這也應歸因于超松弛引起的振蕩。3.2和3.3節(jié)將深入討論這些推測。

圖4 不同δt的模擬精度L2誤差范數(shù)變化曲線Fig.4 L2 error norms of different δt values

3.2 壓縮性誤差的討論

X75T50、X75T100 和X75T200 這三個工況的模擬精度有較大差異，但其風速模擬結果并未顯示出明顯的數(shù)值振蕩，且三者的網(wǎng)格設置完全一致僅δt不同。這表明它們的精度差異極有可能是由于不同δt導致的壓縮性誤差所引起，于是本節(jié)選擇這三個工況分析壓縮性誤差。

圖5顯示了三個工況中各區(qū)域的時間平均密度ρˉ與初始值ρ0相比的相對偏差()-ρ0ρ0。如1.3節(jié)所述，LBM-LES 是一種偽可壓縮模擬方法，即使在模擬同一個不可壓縮問題時，不同的δt設置亦會導致密度變化顯著。X75T50的δt較大，導致ρˉ明顯偏離了初始值，尤以入口附近區(qū)域更為明顯。

圖6顯示了所有工況中整個模擬域的空間平均密度相對于初始值的偏差。該偏差反映了LBMLES 計算中密度的壓縮程度。隨著δt減小一半，空間平均密度的差異幾乎呈指數(shù)衰減。當密度偏差小于0.5 %后逐漸達到穩(wěn)定。此時的密度接近初始值，可忽略壓縮性的影響。

圖6 所有工況中全空間平均密度與初始值的相對偏差變化Fig.6 Deviation between spatial-averaged density and initial value of all cases

垂直截面上量綱一化uˉ及M的計算結果如圖7所示。在X75T50 中，來自入口的氣流有遠離天花板向下的趨勢，導致該區(qū)域的模擬結果精度較差。據(jù)圖7b所示，該區(qū)域M大于該工況其他區(qū)域及其他工況的M，并超過了0.3。這與Krüger 等［17］的研究一致。該研究建議LBM-LES模擬場中的由ulb定義的M不應大于0.3，否則將產(chǎn)生較顯著的壓縮性誤差。同時，圖5 顯示該區(qū)域密度變化劇烈。這再次表明較大的密度梯度將會導致顯著的壓縮性誤差。隨著M的降低，入口處的氣流趨于水平，表明壓縮性誤差可通過降低M得到一定的補償。X75T50的其他區(qū)域或其他工況中M均小于0.3，故可忽略uˉ的壓縮性誤差。同時，X150工況組中M＜0.3始終成立，表明X150 工況組的壓縮性誤差均不明顯，因而X150 工況組中模擬精度并未隨著δt的降低而提高。

圖5 部分工況的中央垂直平面量綱一化時間平均密度偏差( -ρ0 )ρ0分布Fig.5 Deviation of time-averaged density(-ρ0 )ρ0 at central vertical section in some cases

圖7 部分工況中垂直截面上量綱一化uˉ及M分布Fig.7 Vertical distribution of normalized uˉand M of some cases

以上討論證實，在使用LBM-LES 求解室內(nèi)湍流時，過大的M會導致速度場產(chǎn)生明顯的壓縮性誤差。通過調(diào)整δt可減小M，以補償誤差。為了減少壓縮性誤差造成的影響，應盡量保證流場中最大風速區(qū)域的M＜0.3，即值得注意的是，M是由ulb定義的格子玻爾茲曼單位的參數(shù)。即使模擬問題原型相同，也可以通過使用不同的δt改變M，這與物理場中由真實速度定義的M不同。

3.3 超松弛導致的數(shù)值振蕩討論

3.1 節(jié) 圖3 顯 示，在X75T400、X75T800、X75T1600工況和X150工況組中的大多數(shù)工況中均觀察到明顯數(shù)值振蕩，這可能由于超松弛碰撞模式造成，本節(jié)將對此進行分析。圖8 顯示了所有工況的f0在點(x，y，z)=(H，H/2，H/2)處當流場達到充分發(fā)展后某個1 s 周期（第1 080～1 081 s）內(nèi)的松弛狀況。選擇點(x，y，z)=(H，H/2，H/2)是因為在該點處X75T50、X75T100 和X75T200 中沒有明顯的振蕩，而在其他工況中發(fā)生振蕩，即存在一個振蕩發(fā)生的臨界狀況。圖中f0(t)/f0eq(t)表明某一t時刻f0與f0eq間的大小關系

圖8 所有工況在(x，y，z)=(H，-H/2，H/2)位置處1 s內(nèi)f0的超松弛現(xiàn)象Fig.8 Over relaxation phenomenon in one second at(x,y,z)=(H,-H/2,H/2)of all cases

在所有工況中，f0圍繞f0eq來回擺動，表明所有工況都是超松弛碰撞模式，而不是理想的亞松弛模式，此時f0并不向f0eq呈指數(shù)衰減。這一結果證實了在湍流中，超松弛碰撞模式比亞松弛更常見。在X75 工況組中，f0擺動的“頻率”隨δt的減小而增大。同時擺動“幅度”隨δt減小而減小。從X75T50 到X75T200，擺動過程似乎沒有形成數(shù)值振蕩，而應該是由湍流脈動導致。然而在X75T400、X75T800 和X75T1600 中，擺動過程演化為可見的高頻振動，與速度場的模擬結果發(fā)生數(shù)值振蕩的狀況一致。類似地，在X150組中，當δt降低到一定程度時，分布函數(shù)形成了高頻振動，最終形成了速度場發(fā)生數(shù)值振蕩，這在X150T1600至X150T3200之間尤為明顯。

由此可見，當δt減小到一定程度時，超松弛碰撞模式所對應的湍流脈動會最終演化成分布函數(shù)的高頻振動，并最終導致宏觀速度場的數(shù)值振蕩。然而，超松弛碰撞模式何時演化為高頻振動則較為復雜，其與局部湍流的流動模式、網(wǎng)格尺寸、流體性質(zhì)及碰撞算子等都有很大關系，并非只與δt線性相關，故較難判斷數(shù)值振蕩的臨界δt。一般建議在消除壓縮性誤差的前提下盡量增大δt以避免數(shù)值振蕩。

4 結論

本文分析并總結了采用含BGK 碰撞算子的LBM-LES模擬風環(huán)境問題時，時間步長δt對模擬精度的影響，并以室內(nèi)氣流模擬案例對其進行了定量討論。主要結論如下：

（1）LBM-LES 是一種偽可壓縮方法，在處理不可壓縮問題時模擬域中的密度在模擬過程中會發(fā)生變化。過大的δt會使得密度變化較大，導致速度場產(chǎn)生壓縮性誤差。較小的δt可以減小壓縮性誤差。

（2）在模擬湍流時，BGK 碰撞算子通常表現(xiàn)為超松弛碰撞模式，即分布函數(shù)表現(xiàn)為一定程度的擺動。過小的δt會導致擺動會演化成高頻振動，最終使得速度場發(fā)生數(shù)值振蕩。該現(xiàn)象在網(wǎng)格分辨率相對較高時更容易產(chǎn)生。

（3）在確定網(wǎng)格分辨率（如網(wǎng)格獨立性測試）后，δt應足夠小以滿足最大風速區(qū)域M＜0.3（即δt≤以減小壓縮性誤差。在此基礎上盡量采用較大的δt以防止數(shù)值振蕩的發(fā)生。

本文僅粗略確定了δt的取值上限，今后的工作將著重于如何確定δt的合理范圍，并建立δt與其他物理量之間的定量關系。此外，對應于風環(huán)境中高Re 問題的模擬，通常采用比BGK 更為復雜、魯棒性更高的碰撞算子（如MRT、cumulant LBM 等），δt的變化對這些碰撞算子的影響也應予以進一步研究。