非均勻矩形網(wǎng)格的局部網(wǎng)格細(xì)化LBM 算法研究1)

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

* (西北工業(yè)大學(xué)航空學(xué)院,西安 710072)

? (中國空氣動(dòng)力研究與發(fā)展中心結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室,四川綿陽 621000)

** (翼型、葉柵空氣動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室,西安 710072)

引言

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

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

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

1 25 點(diǎn)插值LBM 計(jì)算模型

本文研究使用最常用的LBGK-D2Q9 模型[28].其平衡態(tài)分布函數(shù)的構(gòu)造如下式所示

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

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

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

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

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

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

由于本文的計(jì)算網(wǎng)格是非均勻的矩形網(wǎng)格,每一層網(wǎng)格的網(wǎng)格間距都不相同,傳統(tǒng)LBM 中基于格點(diǎn)和格點(diǎn)之間遷移條件無法滿足.

如圖3 所示,點(diǎn)A和O均為網(wǎng)格格點(diǎn).由于點(diǎn)A所在的網(wǎng)格層和點(diǎn)O所在的網(wǎng)格層之間的網(wǎng)格間距不同,?xOA

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

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

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

圖5 網(wǎng)格對(duì)比Fig.5 Comparison between different mesh

2 數(shù)值模擬

2.1 計(jì)算域構(gòu)建及邊界條件

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

為了便于研究,本文非均勻矩形網(wǎng)格的建立借鑒了等比數(shù)列的控制函數(shù),將計(jì)算域沿橫向和縱向劃分為4 個(gè)對(duì)稱的子區(qū)間.每個(gè)子區(qū)間網(wǎng)格間距由下式控制

其中,N為網(wǎng)格點(diǎn)個(gè)數(shù),σ 是等比控制函數(shù)的公比,在本文中當(dāng)N=500 時(shí),σ 為0.005 213 64.?xfinest=?yfinest定義了貼近物面處第一層網(wǎng)格的網(wǎng)格間距(最細(xì)網(wǎng)格).

2.2 邊界條件處理

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

如圖6 所示點(diǎn)A,B和C為物面邊界點(diǎn),D點(diǎn)為流場(chǎng)點(diǎn).根據(jù)LBM 的演化原理可知在每次演化之前需要知道每個(gè)點(diǎn)的分布函數(shù),對(duì)于B點(diǎn),其分布函數(shù)可分為平衡態(tài)和非平衡態(tài)兩部分

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

以物面邊界點(diǎn)B為例,該點(diǎn)的平衡態(tài)分布函數(shù)可用該點(diǎn)的宏觀物理量來構(gòu)造,如果該點(diǎn)的宏觀物理量未知,則由D點(diǎn)的相應(yīng)值代替.而非平衡態(tài)分布函數(shù)則由D點(diǎn)的非平衡態(tài)分布函數(shù)來近似代替

3 計(jì)算結(jié)果與分析討論

3.1 定常流動(dòng)數(shù)值模擬

本文以頂蓋驅(qū)動(dòng)方腔內(nèi)流為例,針對(duì)本文RLBML25 算法在定常(Re=1000 和Re=5000)和非定常狀態(tài)(Re=8800 和Re=9500)下開展數(shù)值模擬分析驗(yàn)證研究.

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

表1 不同計(jì)算結(jié)果對(duì)比 (Re=1000)Table 1 Different results comparison (Re=1000)

類似表1,表2 列出了當(dāng)Re=5000 時(shí)的數(shù)據(jù)對(duì)比.數(shù)據(jù)結(jié)構(gòu)和內(nèi)容與表1 一致.再次證明RLBML25 算法對(duì)于定常流動(dòng)數(shù)值模擬表現(xiàn)良好.

表2 不同計(jì)算結(jié)果對(duì)比 (Re=5000)Table 2 Different results comparison (Re=5000)

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

圖8 描述了沿頂蓋分布的剪切力fs,總體來講,本文RLBM-L25 算法網(wǎng)格分辨率分別為500 和600 且雷諾數(shù)等于1000 時(shí)的計(jì)算結(jié)果與傳統(tǒng)LBM算法網(wǎng)格分辨率為1024 時(shí)的結(jié)果較為吻合,但是在頂蓋左右兩側(cè)附近仍存在一定量級(jí)的誤差且不會(huì)因?yàn)椴逯蹈袷降母淖兌?因此作者認(rèn)為這種情況的出現(xiàn)首先跟頂蓋上網(wǎng)格分布數(shù)量有關(guān),其次跟網(wǎng)格局部細(xì)化的控制方程也有關(guān).

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

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

圖10 擾動(dòng)衰減率曲線Fig.10 History of perturbation decay rate ε of velocity ux

為了進(jìn)一步討論分析本文算法的表現(xiàn),還開展了基于不同網(wǎng)格分辨率的研究.以傳統(tǒng)算法SLBM且網(wǎng)格分辨率為1024 時(shí)的計(jì)算結(jié)果作為對(duì)比.

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

表3 不同網(wǎng)格數(shù)量計(jì)算結(jié)果對(duì)比Table 3 Comparison of different resolutions

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

4.2 非定常流動(dòng)分?jǐn)?shù)值模擬

為了進(jìn)一步探討本文算法(RLBM-L25)對(duì)于非定常流動(dòng)的表現(xiàn),本文開展了非定常周期性流動(dòng)的數(shù)值模擬.選取兩個(gè)非定常周期性雷諾數(shù)8800 和9500 作為數(shù)值模擬計(jì)算條件.

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

表5 不同網(wǎng)格分辨率計(jì)算結(jié)果對(duì)比Table 5 Comparison of different resolutions

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

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

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

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

5 結(jié)論

本文提出了一種基于25 點(diǎn)拉格朗日插值的非均勻矩形網(wǎng)格LBM 算法.針對(duì)本文算法以方腔頂蓋驅(qū)動(dòng)內(nèi)流為例開展了定常和非定常流動(dòng)的數(shù)值模擬驗(yàn)證研究,得出主要結(jié)論如下.

(1)本文算法對(duì)于定常流動(dòng)數(shù)值模擬表現(xiàn)優(yōu)異,網(wǎng)格細(xì)化策略可以保證物面處和流動(dòng)突變區(qū)的流動(dòng)細(xì)節(jié)的捕捉.

(2)本文算法對(duì)于非定常周期性流動(dòng)數(shù)值模擬表現(xiàn)良好,同樣可以保證物面處和流動(dòng)突變區(qū)的流動(dòng)細(xì)節(jié)的捕捉.相對(duì)基于均勻正方形網(wǎng)格(1024×1024)傳統(tǒng)LBM 的計(jì)算結(jié)果,雖然存在誤差,但總體控制在可接受的范圍內(nèi).

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

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

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

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