非均勻矩形網格的局部網格細化LBM 算法研究1)

安 博 孟欣雨 楊雙駿 桑為民,2)

* (西北工業大學航空學院,西安 710072)

? (中國空氣動力研究與發展中心結冰與防除冰重點實驗室,四川綿陽 621000)

** (翼型、葉柵空氣動力學重點實驗室,西安 710072)

引言

格子玻爾茲曼方法[1-3](LBM)在經歷幾十年的快速發展后已逐漸成為流體力學研究中的重要分支和熱點話題.同時在計算流體力學(CFD)領域中得到廣泛應用,比如不可壓縮流動[4]、可壓縮流動[5]、多孔介質流[6]、氣動聲學[7]、多相流[8]、燃燒[9]、滲流[10]等.由于其經典的碰撞遷移理論完美地契合均勻正方形網格的網格結構,因此在LBM 相關數值模擬中應用非常廣泛.雖然基于均勻正方形網格的單松弛LBM (傳統LBM)在數值模擬方面表現不俗,但由于其魯棒性和數值穩定性較差,尤其當松弛因子 τ=0.5+3UL/(Re?x)接近0.5 時,數值穩定性逐漸喪失[1-3],因此其應用范圍基本僅限于低雷諾數Re的數值模擬.為了應對逐漸增大的雷諾數,就必須減小網格間距 ?x,同時還能滿足物面處以及流動變化劇烈區域流動細節的捕捉.而均勻正方形網格的細化必然導致網格總量的驟然增加,使得計算效率大幅下降,且對計算設備的要求也急劇提高.極大地限制傳統LBM 的應用與發展.

作者在之前的研究中提出了基于樹網格的LBM算法[11-12],可以在保證計算精度,滿足局部網格細化要求的前提下,大幅提高傳統LBM 的計算效率.此外,眾多學者也開展了傳統LBM 局部網格細化的算法優化研究[13-27].比如有限差分LBM[13-14]、插值補充LBM[15]、多重網格LBM[16-18]、多區塊LBM[19-21]、樹網格LBM[22-25]以及均勻矩形網格LBM[26-27]等.Mei 等[13]以及Guo 等[14]借鑒傳統的求解N-S 方程的手段,將有限差分方法與LBM結合.由于控制方程,即格子玻爾茲曼方程不再通過碰撞遷移演化求解,轉而依靠差分格式逼近.所以網格結構的使用也變得較為靈活,隨著貼體網格的引入,在網格整體數量不變的情況下解決了網格物面處的局部加密的問題.但是貼體網格的使用會帶來復雜的坐標變換以及在求解過程中需要構造高精度差分格式,總體來講未能有效提升計算效率.同樣地,基于貼體網格,He 等[15]使用的插值補充格式,依然保留了LBM的精髓,即分布函數的碰撞遷移過程,然后通過插值來近似求解碰撞遷移后虛擬點的分布函數,但計算過程比較復雜且插值格式精度較低.Marvriplis[16]、Patil 等[17]和王興勇等[18]將傳統CFD 中的多重網格技術應用于LBM 的數值模擬中.其核心思想是通過粗細網格的信息交換,使得低頻誤差在粗網格快速降低,高頻誤差在細網格得到光滑處理,從而提升計算效率.總體來講,多區塊LBM[19-21]和樹網格LBM[22-25]算法基本是殊途同歸,雖然每個學者的實現手段不盡相同,但是核心思想是一致的.即對流場內部的不同區域采取局部加密.這些區域往往是流動變化較為劇烈的區域或者是物面流動區.Zhou[26-27]提出了基于均勻矩形網格的LBM 算法,不同于本文非均勻矩形網格LBM 算法,Zhou 的計算模型沒有插值過程,用重新構建的矩形離散速度模型以及對應的平衡態分布函數取代傳統的正方形離散速度模型(見圖1),可以保證分布函數以矩形方式完成格點之間的碰撞遷移.雖然該方法可以在一定程度上減少總體網格數量,但局部網格細化效果不佳且不具備普遍的適用性.以上國內外的相關研究對本文工作予以很大的啟發,尤其是插值補充LBM 算法[15]中虛擬點分布函數的插值求解過程,對本文基于非均勻矩形網格25 點拉格朗日插值LBM 算法的構建提供了思路.

圖1 LBGK-D2Q9 離散速度模型Fig.1 Discrete velocity model of LBGK-D2Q9

1 25 點插值LBM 計算模型

本文研究使用最常用的LBGK-D2Q9 模型[28].其平衡態分布函數的構造如下式所示

其中,ωα是離散速度模型中不同離散方向的權系數;ρ為流體密度;u為流體速度矢量;eα是離散速度,其大小和方向構造如下所示

其中,c=?x/?t=1 為格子速度;cs為格子聲速;?x為網格步長;?t為時間步長.圖1 介紹了LBGK-D2Q9計算模型的離散速度模型.詳細說明了分布函數fα的碰撞遷移方式.

式(3)定義了LBGK 計算模型中碰撞遷移理論的控制方程,其中r為空間向量;ξ 為速度向量;τ 為松弛時間;分別為分布函數和平衡態的分布函數.

下式定義了流場宏觀物理量,速度u和密度 ρ 與微觀粒子分布函數fα之間的聯系.每一次演化過后,流場的速度和密度都會通過下式得到更新

傳統的LBM 算法基于均勻正方形網格,根據其碰撞遷移理論,分布函數的遷移過程必須保證在網格格點上.因為在傳統的LBM 算法中,格子速度 ψ.如圖2 所示,點A和O為網格格點.點O處的流體粒子,其2 方向的分布函數經過一個時間步長 ψ 的演化會遷移至A點,遷移距離為網格間距 ω,即.其他方向的分布函數也遵循同樣的規則遷移至O點周圍的各點.由此可見均勻正方形網格天然匹配LBM 的演化機理.通常情況下,為了得到較為理想的數值模擬結果,尤其是捕捉物面處以及流動劇烈變化區的流動細節,往往需要極為細膩的網格.這就使得均勻正方形網格的網格總量急劇增加,極大地降低了傳統LBM 的計算效率和魯棒性,限制了LBM 的應用.因此,本文采用了非均勻矩形網格用于流場局部加密既,保證了物面處和流動變化劇烈區的流場細節捕捉,同時相較于均勻正方形網格,在保證計算精度的前提下又能提升LBM 的計算效率.

圖2 均勻直角網格的遷移演示Fig.2 Streaming on uniform Cartesian grid

由于本文的計算網格是非均勻的矩形網格,每一層網格的網格間距都不相同,傳統LBM 中基于格點和格點之間遷移條件無法滿足.

如圖3 所示,點A和O均為網格格點.由于點A所在的網格層和點O所在的網格層之間的網格間距不同,?xOA

圖3 非均勻矩形網格的遷移演示Fig.3 Streaming on non-uniform rectangular grid

圖4,以A33點5 方向的分布函數f5為例,介紹了本文25 點拉格朗日插值格式的構造過程.由圖可見點Aij(i=1,2,···,5;j=1,2,···,5)均為計算域的實際網格點,而點aij(i=1,2,···,5;j=1,2,···,5)為沿5 方向遷移后的虛擬點,并且包括點A33以及其周圍的24 個點都參與了插值格式的構建.網格沿橫(ξ1,ξ2,ξ3和 ξ4)、縱(η1,η2,η3和 η4)向均有4 層網格且網格間距均不相同.?=?t·c為流體粒子單位時間步長 ?t內移動的距離(c=1.0),同時也是最細網格間距(?=?xfinest).本文沿橫、縱向靠近物面的第一層最細網格為 ?xfinest=1/1024,見圖5(b)中計算域四角黑色區域.以5 方向的分布函數為例,下式給出了本文25 點拉格朗日插值格式

圖4 非均勻矩形網格的遷移演示(α=5)Fig.4 Streaming on discrete direction (α=5)

圖5 網格對比Fig.5 Comparison between different mesh

2 數值模擬

2.1 計算域構建及邊界條件

為了驗證本文基于25 點拉格朗日插值的非均勻矩形網格LBM 算法,選取經典的頂蓋驅動方腔內流為基準算例(Ma=0.173 2 且Re<10 000).本文數值模擬最細網格步長(?xfinest=1/1024)的選取借鑒了作者之前研究工作[29-32]中針對腔體內流的網格獨立性驗證.基于此網格尺度,發現當Re<10 000 時能夠保證y+<0.1,數值模擬結果可信且能夠捕捉邊界層附近的流場信息.如圖5 所示,展示了頂蓋驅動方腔內流的網格結構,圖5(a)和圖5(b)分別對應均勻正方形網格和非均勻矩形網格.流動信息采集點位于腔體中心P(x=L/2,y=L/2)處.

為了便于研究,本文非均勻矩形網格的建立借鑒了等比數列的控制函數,將計算域沿橫向和縱向劃分為4 個對稱的子區間.每個子區間網格間距由下式控制

其中,N為網格點個數,σ 是等比控制函數的公比,在本文中當N=500 時,σ 為0.005 213 64.?xfinest=?yfinest定義了貼近物面處第一層網格的網格間距(最細網格).

2.2 邊界條件處理

本文涉及的邊界均為平直邊界(方腔4 個邊).其中頂蓋為驅動邊界,其余三邊均為物面邊界.由于非平衡態外推格式[33]易于代碼編寫且非常適用于平直邊界點的整體處理.本文延續之前研究中已成熟應用的非平衡態外推格式,將邊界點上的分布函數分為平衡態和非平衡態兩部分.其中平衡態部分由平衡態分布函數的定義近似獲得,而非平衡態部分則用非平衡態外推法求解.

如圖6 所示點A,B和C為物面邊界點,D點為流場點.根據LBM 的演化原理可知在每次演化之前需要知道每個點的分布函數,對于B點,其分布函數可分為平衡態和非平衡態兩部分

圖6 平直物面邊界Fig.6 Straight wall boundary

以物面邊界點B為例,該點的平衡態分布函數可用該點的宏觀物理量來構造,如果該點的宏觀物理量未知,則由D點的相應值代替.而非平衡態分布函數則由D點的非平衡態分布函數來近似代替

3 計算結果與分析討論

3.1 定常流動數值模擬

本文以頂蓋驅動方腔內流為例,針對本文RLBML25 算法在定常(Re=1000 和Re=5000)和非定常狀態(Re=8800 和Re=9500)下開展數值模擬分析驗證研究.

表1 介紹了當Re=1000 時本文的計算結果和其他文獻[34-35]的數據對比.第1 列P 代表流場中渦漩的位置,LB 表示左下側次級渦;RB 表示右下側次級渦;PC 代表中心主渦(見圖7).對應每個渦結構給出了4 個流動參數,分別為渦心橫坐標(X)、縱坐標(Y)、流函數(ψ)和渦量(ω).由表1 可知,本文提出的RLBM-L25 算法其結果與其他文獻的數據吻合較好,可以輸出可信的數值模擬結果.

表1 不同計算結果對比 (Re=1000)Table 1 Different results comparison (Re=1000)

類似表1,表2 列出了當Re=5000 時的數據對比.數據結構和內容與表1 一致.再次證明RLBML25 算法對于定常流動數值模擬表現良好.

表2 不同計算結果對比 (Re=5000)Table 2 Different results comparison (Re=5000)

如圖7(a)和圖7(b)所示,展示了定常狀態下的流場拓撲結構.可以明顯觀察到經典方腔驅動內流的典型流動特征,與其他文獻[34-39]的數據完全一致.中心主渦覆蓋大部分計算域從而主導整個流場,隨著雷諾數的增加,次級渦逐漸增大并伴隨邊角處出現細碎渦的現象.圖7(c)展示了Re=1000 時流場的局部放大圖,左側是基于均勻正方形網格的傳統LBM 算法(SLBM)的計算結果,右側是本文RLBML25 的計算結果.通過對比可以看出,本文基于非均勻矩形網格的網格加密方法可以保證物面處和流動突變區的流動細節捕捉.圖7(d)和圖7(e)分別展示了Re=1000 時基于均勻正方形網格的經典單松弛LBM(SLBM)且網格分辨率為1024 和基于非均勻矩形網格插值LBM(RLBM-500)且網格分辨率為500 時的渦量圖對比,如圖所示,渦量圖吻合得較好.

圖8 描述了沿頂蓋分布的剪切力fs,總體來講,本文RLBM-L25 算法網格分辨率分別為500 和600 且雷諾數等于1000 時的計算結果與傳統LBM算法網格分辨率為1024 時的結果較為吻合,但是在頂蓋左右兩側附近仍存在一定量級的誤差且不會因為插值格式的改變而消除,因此作者認為這種情況的出現首先跟頂蓋上網格分布數量有關,其次跟網格局部細化的控制方程也有關.

圖8 頂蓋剪切力分布Fig.8 Shear force fs distribution along the lid

圖9 速度均方根誤差曲線Fig.9 History of root mean square error eRMS of velocityux

圖10 擾動衰減率曲線Fig.10 History of perturbation decay rate ε of velocity ux

為了進一步討論分析本文算法的表現,還開展了基于不同網格分辨率的研究.以傳統算法SLBM且網格分辨率為1024 時的計算結果作為對比.

由表3 可見,第1 列給出了不同的網格分辨率(240～620).第2 列介紹了當Re=1000 時,信息采集點水平速度分量ux的收斂值.第3 列展示了每個網格分辨率相對于SLBM-1024 算法計算結果的相對誤差 ?ux(%).由表可見,隨著網格分辨率的增加,相對誤差逐漸減小,這與均勻正方形網格數值模擬所表現出來的特性一致.當網格分辨率從500 增加至620 時,相對誤差從1.7%減小至0.99%.但是由于網格數量的增加會降低計算效率,同時1.7%的相對誤差小于2.0%,是本文可接受的誤差范圍.綜合考慮,本文所有的數值模擬都采用500 的網格分辨率.

表3 不同網格數量計算結果對比Table 3 Comparison of different resolutions

此外,本文還開展了不同插值格式之間的對比研究.同樣地我們以傳統算法SLBM 且網格分辨率為1024 時的計算結果作為對比.結合表4,可以看到第1 列列舉了不同的插值格式,分別為拉格朗日9 點插值(RLBM-L9)、牛頓9 點插值(RLBM-N9)、埃爾米特9 點插值(RLBM-H9)以及拉格朗日25 點插值(RLBM-L25).第2 列介紹了當Re=1000 時,信息采集點水平速度分量ux的收斂值.第3 列展示了每個插值格式相對于SLBM-1024 計算結果的相對誤差 ?ux(%).由表可見,本文RLBM-L25 算法相對于其他插值算法表現優異.

4.2 非定常流動分數值模擬

為了進一步探討本文算法(RLBM-L25)對于非定常流動的表現,本文開展了非定常周期性流動的數值模擬.選取兩個非定常周期性雷諾數8800 和9500 作為數值模擬計算條件.

表5 展示了不同網格分辨率(RLBM-L25-500和RLBM-L25-620)雷諾數為9500 時計算結果,對SLBM-1024 的計算結果.第1 列數據(Data)分別列舉了周期性速度曲線的最大值(uxmax)、最小值(uxmin)、平均值(uxmean)以及平均值相對SLBM-1024 計算結果的相對誤差 ?uxmean.通過表5 我們發現,對于非定常周期性流動,隨著網格分辨率的增大,相對誤差逐漸減小.相對于網格分辨率620,500 的分辨率產生較大的相對誤差,但仍然在誤差允許范圍內(?uxmean<5%).作者認為,這個網格分辨率可以保證非定常周期性流動的數值模擬精度.

表5 不同網格分辨率計算結果對比Table 5 Comparison of different resolutions

如圖11 所示,介紹了流場在一個完整周期(T=1/f,f=0.432 6)內某一時刻的瞬時流線圖.其表現出的流場典型特征與作者之前針對方腔內流的研究[41-42]完全一致.

圖11 某一時刻周期性解的瞬時流線圖Fig.11 Snapshots of periodic solutions at a given moment

圖12 介紹了雷諾數為9500 時,信息采集點基于不同插值格式的計算結果對比.主畫面展示了速度頻譜圖對比,可見不同的插值格式對于周期性流動的振蕩頻率(居中內嵌圖)的預測基本一致,且本文RLBM-L25 算法表現優異.從右側內嵌圖可以看出,在預測周期性流動的速度振幅和平均值時,與SLBM-1024 的計算結果存在差異.相較于其他插值格式,RLBM-L25 算法的表現最好,即便如此,與SLBM-1024 的計算結果的相對誤差為3.5%左右.作者認為導致這一誤差的可能原因有兩個.其一,與本文網格細化策略有關,比如網格分辨率,網格間距控制函數等.此外插值格式也是一個需要考慮的因素,作者認為更精確的插值格式能很大程度上消除或降低這一誤差.

圖12 某一時刻周期性解的瞬時流線圖Fig.12 Snapshots of periodic solutions at a given moment

5 結論

本文提出了一種基于25 點拉格朗日插值的非均勻矩形網格LBM 算法.針對本文算法以方腔頂蓋驅動內流為例開展了定常和非定常流動的數值模擬驗證研究,得出主要結論如下.

(1)本文算法對于定常流動數值模擬表現優異,網格細化策略可以保證物面處和流動突變區的流動細節的捕捉.

(2)本文算法對于非定常周期性流動數值模擬表現良好,同樣可以保證物面處和流動突變區的流動細節的捕捉.相對基于均勻正方形網格(1024×1024)傳統LBM 的計算結果,雖然存在誤差,但總體控制在可接受的范圍內.

(3)相較于傳統的基于均勻正方形網格的LBM 算法,本文算法在保證計算精度的前提下可以提升計算效率2～5 倍.

(4)針對本文非均勻矩形網格,當網格分辨率逐漸增加時,基于理論解的相對誤差逐漸減小,與均勻正方形網格LBM 算法的特性一致.

(5)本文拉格朗日25 點插值格式,相較于格式拉格朗日9 點插值、牛頓9 點插值和埃爾米特9 點插值,表現更好.

(6)為了消除或降低誤差,可以改良網格細化策略,如提高網格分辨率和優化網格間距控制函數等.此外,采用更為精確的插值格式.