劉春友 李作旭 王連平 ,?, ,2)

* (南方科技大學(xué)力學(xué)與航空航天工程系,廣東深圳 518055)

? (南方科技大學(xué)復(fù)雜流動及軟物質(zhì)研究中心,廣東深圳 518055)

** (廣東省湍流基礎(chǔ)研究與應(yīng)用重點(diǎn)實(shí)驗(yàn)室,廣東深圳 518055)

引言

在實(shí)際的工程應(yīng)用中,各種各樣的流體力學(xué)問題隨處可見,而如何準(zhǔn)確且高效的描述其對應(yīng)的流動特征一直是工業(yè)界和科學(xué)界重點(diǎn)關(guān)注的問題.近半個(gè)世紀(jì)以來,計(jì)算機(jī)技術(shù)的迅猛發(fā)展使得大規(guī)模數(shù)值求解復(fù)雜問題成為可能,而數(shù)值模擬作為一種有效的研究手段也得到了更加廣泛的關(guān)注和應(yīng)用.不同于直接數(shù)值求解納維-斯托克斯(Navier-Stokes,N-S)方程的傳統(tǒng)計(jì)算流體力學(xué)方法,介觀計(jì)算流體動力學(xué)方法是直接數(shù)值求解基于分子動力學(xué)理論的Boltzmann 方程,而格子玻爾茲曼方法(lattice Boltzmann method,LBM)則是近30 年來出現(xiàn)的一種較為高效的介觀計(jì)算流體動力學(xué)方法.它較高的計(jì)算效率和靈活性使其可以適用于各種高度復(fù)雜的物理現(xiàn)象,如湍流[1]、燃燒[2]、多相流[3]、磁流體[4]和粒子懸浮[5-6]等問題.由于其數(shù)值耗散較低,故LBM也常被應(yīng)用于氣動聲學(xué)模擬中[7-8].同時(shí),LBM 也憑借演化過程中“碰撞-遷移”的局部性,使得其非常適合于大規(guī)模的并行計(jì)算[9-10].除此之外,許多學(xué)者為了充分利用LBM 的優(yōu)點(diǎn),還通過反設(shè)計(jì)的思想在原始的Boltzmann 方程中添加一些特定的源項(xiàng),以使得修改后的LBM 可以滿足特定需求[11-12],而這也使得LBM 方法的適用范圍得到了極大的擴(kuò)展.

流體的流動特點(diǎn)一般是隨著Reynolds 數(shù)的增加,會產(chǎn)生越來越精細(xì)的動力學(xué)結(jié)構(gòu).而工業(yè)領(lǐng)域里涉及到的流動問題大多是高Reynolds 數(shù)的湍流,通過直接數(shù)值計(jì)算來解析所有的湍流結(jié)構(gòu)往往需要較高的計(jì)算成本.因此,為了在有限的計(jì)算資源下獲得較為準(zhǔn)確的結(jié)果,往往需要提高局部分辨率來解析這些局部流域內(nèi)包含的精細(xì)流動結(jié)構(gòu).但由于LBM算法在時(shí)空離散時(shí)是沿著特征線進(jìn)行的離散,故在給定離散速度集后LBM 中的時(shí)間步長和空間步長便通過離散速度相互耦合,這一特點(diǎn)導(dǎo)致標(biāo)準(zhǔn)的LBM只能使用均勻的直角網(wǎng)格.

為了克服標(biāo)準(zhǔn)LBM 只能使用均勻網(wǎng)格這一缺陷,Filippova 等[13]提出了基于LBM 的局部網(wǎng)格加密(FH)方法.該方法首先需要使用均勻的直角粗網(wǎng)格(背景網(wǎng)格)覆蓋整個(gè)流場,然后在需要進(jìn)行加密的區(qū)域上對網(wǎng)格進(jìn)行加密,不同網(wǎng)格分辨率區(qū)域下仍使用標(biāo)準(zhǔn)的LBM 演化關(guān)系且使用同一離散速度集.而粗細(xì)網(wǎng)格之間的耦合則通過流場物性參數(shù)以及分布函數(shù)間的轉(zhuǎn)換實(shí)現(xiàn),如可以通過保證粗細(xì)網(wǎng)格間運(yùn)動黏度系數(shù)一致來得到粗細(xì)網(wǎng)格間格子松弛系數(shù)的關(guān)系,然后通過保證粗細(xì)網(wǎng)格間宏觀量一致得到分布函數(shù)的轉(zhuǎn)換關(guān)系.文中以二維圓柱擾流為測試算例,在使用局部網(wǎng)格加密技術(shù)后得到了較好的數(shù)值結(jié)果,這也驗(yàn)證了網(wǎng)格局部加密的有效性.局部網(wǎng)格加密算法的提出對于標(biāo)準(zhǔn)LBM 的意義是重大的,這意味著其可以通過局部加密來節(jié)省大量的計(jì)算成本,也使得LBM 被大規(guī)模應(yīng)用于工業(yè)界成為了可能.

與FH 方法不同,Lin 等[14]建議不對粗細(xì)網(wǎng)格間分布函數(shù)的非平衡部分做任何縮放,而是直接交換粗細(xì)網(wǎng)格的分布函數(shù).Dupuis 等[15]比較了FH 方法和Lin 的方法后發(fā)現(xiàn),當(dāng)考慮非平衡部分的轉(zhuǎn)換時(shí),可以得到更精確的數(shù)值結(jié)果.而Touil 等[16]更是直接指出,Lin 的方法實(shí)際上是混淆了連續(xù)分布函數(shù)和離散分布函數(shù)之間的區(qū)別,其直接在粗細(xì)網(wǎng)格間交換分布函數(shù)會導(dǎo)致該方法的精度下降.此外,Rohde等[17]提出了一種基于有限體積思想且可以適用于任何碰撞模型的局部網(wǎng)格加密算法,但其也犯了與Lin 相同的錯(cuò)誤,即不對分布函數(shù)做任何轉(zhuǎn)換,這種處理方式在模擬高Reynolds 數(shù)下的流動時(shí),很容易在網(wǎng)格加密界面處出現(xiàn)間斷和非物理震蕩.而許多學(xué)者[18-21]的相關(guān)研究也表明,在局部網(wǎng)格加密算法中,粗細(xì)網(wǎng)格間非平衡態(tài)分布函數(shù)的轉(zhuǎn)換是必要且不可或缺的.

眾所周知,LBM 的靈活性很大程度上來自于可以通過Chapman-Enskog (CE)展開[22]靈活設(shè)計(jì)源項(xiàng).因此,在考慮源項(xiàng)的前提下,建立粗細(xì)網(wǎng)格間分布函數(shù)的轉(zhuǎn)換關(guān)系,對于擴(kuò)展局部網(wǎng)格加密算法在LBM中的應(yīng)用范圍至關(guān)重要.Huang 等[23]基于CE 展開,提出了一種在多松弛碰撞模型(multiple relaxation time,MRT)[24]下利用多層網(wǎng)格技術(shù)計(jì)算被動標(biāo)量的模型.該方法在考慮外力項(xiàng)存在的情況下,在矩空間下執(zhí)行粗細(xì)網(wǎng)格間的數(shù)據(jù)轉(zhuǎn)換.Liou 等[25]同樣也基于CE 展開,提出了一種在BGK (Bhatnagar-Gross-Krook)模型[26]下,同時(shí)考慮外力和標(biāo)量源項(xiàng)的加密方法.值得注意的是,這些方法都是在 CE 展開的基礎(chǔ)上導(dǎo)出的,其推導(dǎo)過程較為復(fù)雜,且在推導(dǎo)中對分布函數(shù)的非平衡態(tài)使用了一階CE 近似,而這有可能限制局部網(wǎng)格加密算法在高階LBM[27]中的應(yīng)用.本文的主要目標(biāo)是,在考慮任意源項(xiàng)存在的前提下,提出一套規(guī)范且簡潔的粗細(xì)網(wǎng)格在重合點(diǎn)處的數(shù)據(jù)轉(zhuǎn)換關(guān)系推導(dǎo)方法.

本文的其余部分按如下順序展開.在第 1 節(jié),梳理BGK 模型和MRT 模型下LBM 的時(shí)空離散過程,尤其是指出在時(shí)空離散過程中離散分布函數(shù)變量的引入,從而為后文局部網(wǎng)格加密算法中的數(shù)據(jù)轉(zhuǎn)換進(jìn)行必要的鋪墊;同時(shí),梳理BGK 模型和MRT 模型下非平衡態(tài)分布函數(shù)的初始化方法,為后文數(shù)值驗(yàn)證中分布函數(shù)的初始化提供依據(jù).第 2 節(jié)則是憑借第 1 節(jié)中的理論基礎(chǔ),分別詳細(xì)給出BGK 模型和MRT 模型下粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系的推導(dǎo)過程.在第 3 節(jié)中,通過對強(qiáng)迫泰勒-格林渦流動和平板泊肅葉流中對流-擴(kuò)散問題進(jìn)行數(shù)值模擬,以驗(yàn)證第 2 節(jié)中轉(zhuǎn)換關(guān)系對復(fù)雜源項(xiàng)的適應(yīng)性;通過對頂蓋驅(qū)動方腔流動進(jìn)行數(shù)值模擬來展示局部網(wǎng)格加密技術(shù)在處理復(fù)雜流動問題方面的優(yōu)勢;通過對一維剪切波問題的模擬,來觀察局部網(wǎng)格加密技術(shù)對LBM 在數(shù)值耗散方面的影響.最后,在第 4 節(jié)給出本文結(jié)論.

1 理論基礎(chǔ)

速度離散后的Boltzmann 方程(discrete velocity Boltzmann equation,DVBE)可以表示如下[28]

這里,ξα表示第 α 個(gè)離散速度點(diǎn)對應(yīng)的離散速度;fα(x,t) 代表t時(shí)刻在位置x處以速度 ξα移動的粒子對應(yīng)的分布函數(shù);?α代表碰撞項(xiàng);Sα為可以代表包括體力項(xiàng)在內(nèi)的任意源項(xiàng).以二維問題中被廣泛使用的D2Q9 速度集為例,如圖1 所示,其離散速度點(diǎn)可以表示如下[28]

圖1 D2Q9 離散速度模型示意圖Fig.1 Schematic diagram of D2Q9 discrete velocity model

不同離散速度點(diǎn)對應(yīng)的權(quán)重為

需要注意的是,本文使用的D2Q9 離散速度集長度未經(jīng)縮放,即格子聲速:Cs=1.

在標(biāo)準(zhǔn)的LBM 中,時(shí)空離散是沿著特征線進(jìn)行的離散,這種離散方法導(dǎo)致了時(shí)間步長和空間步長通過離散的速度集進(jìn)行了耦合.這就使得整個(gè)計(jì)算域都必須使用均勻的直角網(wǎng)格,以保證流體粒子在一個(gè)時(shí)間步長 ?t后可以精確地移動到鄰近的網(wǎng)格節(jié)點(diǎn)上,這也是LBM 需要使用 on-lattice 速度集的原因.而對于其他的介觀方法,如:基于有限體積法對DVBE 經(jīng)行離散的DUGKS[29]等方法,這些方法中時(shí)間步長和空間步長是解耦的,故其可以使用非均勻網(wǎng)格,且既可以使用on-lattice 的速度集也可以使用off-lattice 的速度集.接下來,我們將給出標(biāo)準(zhǔn)LBM 中對DVBE 的時(shí)空離散過程直接對DVBE 左右兩側(cè)沿特征線積分,可得如下結(jié)果

從上式不難看出,分布函數(shù)fα從t時(shí)刻到t+?t時(shí)刻,時(shí)間上推進(jìn)了 ?t而空間上推進(jìn)的長度為?x=ξα?t,這也正是LBM 中時(shí)空耦合的原因.同時(shí)在推導(dǎo)上式的過程中,我們對DVBE 等號右邊項(xiàng)使用了梯形積分公式來近似,其中分別代表?α(x+ξα?t,t+?t),Sα(x+ξα?t,t+?t).不難發(fā)現(xiàn)式 (4)右邊項(xiàng)中包含了t+?t時(shí)刻的未知量,而引入新的變量gα(x,t) 可以消除這種隱式[30],該離散分布函數(shù)變量的表達(dá)式如下

現(xiàn)在可以使用離散分布函數(shù)變量gα對略去高階小量的式 (4) 重新表示

在給定碰撞項(xiàng) ?α以及源項(xiàng)Sα的具體形式后,結(jié)合式 (5) 對式 (6) 簡單整理后便可以得到對應(yīng)碰撞模型下的LBM 演化式.下面以BGK 碰撞模型和MRT 碰撞模型為例,給出兩種模型下LBM 演化式的詳細(xì)推導(dǎo)過程.

1.1 BGK 碰撞模型

速度離散后的BGK 碰撞模型可以直接表示為

式中,q為離散速度集中離散速度點(diǎn)的個(gè)數(shù).而H(n)(ξα) 代表以 ξα為自變量的第n階Hermite 多項(xiàng)式[31],前兩階的Hermite 多項(xiàng)式可以表示如下

其中,式 (8) 中的密度 ρ 和速度u可以由fα以及離散速度 ξα表示如下

而根據(jù)式 (11)和式(12) 以及Hermite 多項(xiàng)式的具體表達(dá)式可知,密度 ρ、速度u以及黏性應(yīng)力張量σ正好分別對應(yīng)連續(xù)分布函數(shù)變量fα的第 0 階、第1階Hermite 矩以及非平衡態(tài)分布函數(shù)的 第2階Hermite 矩

將式 (7)代入到式(5)后,可以反解出

式中,τν代表格子松弛系數(shù),其與網(wǎng)格分辨率有關(guān)

由于時(shí)空離散前后平衡態(tài)部分并沒有發(fā)生任何改變,故為了統(tǒng)一符號,取將式 (14) 和 式(15) 代入式 (7),可以得到BGK 碰撞模型下完全由新變量gα表示的碰撞項(xiàng)?α

將式 (16) 代入式 (6),便可以得到使用BGK 碰撞模型的標(biāo)準(zhǔn)LBM 演化方程

而通過Hermite 展開得到的力項(xiàng)表達(dá)式 (18),也與Guo 等[32]通過反設(shè)計(jì)得到的結(jié)果完全一致.

除此之外,由于式 (17) 的解gα≠fα,故由離散分布函數(shù)變量gα計(jì)算宏觀量的表達(dá)式需要重新給出.通過式 (14),可以得到連續(xù)分布函數(shù)變量fα對應(yīng)的各階Hermite 矩a(n)與離散分布函數(shù)變量gα的各階Hermite 矩間的關(guān)系

進(jìn)而可以通過式 (13)、式(20) 以及式 (24) 得到當(dāng)源項(xiàng)為體力項(xiàng)時(shí),由離散分布函數(shù)變量gα表示的密度 ρ、速度u以及由離散非平衡態(tài)分布函數(shù)變量表示的黏性應(yīng)力張量σ

最后,在經(jīng)過合適的初始化后,式 (17) 便可以顯式推進(jìn).而當(dāng)初始流場中包含非平衡態(tài)時(shí),正確的初始化分布函數(shù)對于一些非定常問題是至關(guān)重要的.其中平衡態(tài)分布函數(shù)可以直接根據(jù)式 (8) 通過初始場的密度以及速度給出,而非平衡態(tài)分布函數(shù)的初始化可以借由Hermite 展開給出[33].對于等溫的流動問題,分布函數(shù)的非平衡態(tài)部分可以直接使用其前兩階Hermite 展開來近似,即

將式 (24) 代入到上式中,便可以得到初始流場在BGK 模型下且考慮源項(xiàng)為體力項(xiàng)時(shí),非平衡態(tài)分布函數(shù)的具體表達(dá)形式

由于我們的目的是保證初始時(shí)刻分布函數(shù)所對應(yīng)的密度、速度以及黏性應(yīng)力張量與初始條件一致,故這里只是展到二階.

以上便是從DVBE 到LBM 的詳細(xì)推導(dǎo)過程.其中注意的是,在離散過程中引入了離散分布函數(shù)變量gα,而且核心演化關(guān)系式 (17) 實(shí)際求解的也是離散分布函數(shù)變量gα,而不是時(shí)空離散前的連續(xù)分布函數(shù)變量fα.此外,從式 (5) 不難看出,離散分布函數(shù)變量gα與時(shí)間步長 ?t相關(guān),進(jìn)而也與網(wǎng)格分辨率有關(guān).故當(dāng)使用式 (17) 作為計(jì)算流域上的演化公式時(shí),即使忽略時(shí)空離散誤差,同一時(shí)空點(diǎn)下不同網(wǎng)格分辨率對應(yīng)的gα仍然是不一致的,所以粗細(xì)網(wǎng)格重合點(diǎn)上的gα并不能直接交換,而這也是局部網(wǎng)格加密算法中不同網(wǎng)格分辨率在重合點(diǎn)處離散分布函數(shù)gα需要轉(zhuǎn)換的原因[16].

1.2 MRT 碰撞模型

MRT 模型是d'Humieres[24]在1992 年提出的一種碰撞模型,該模型可以被視為廣義的BGK 模型.相比于只有一個(gè)松弛參數(shù)的原始BGK 模型,MRT模型中則有著多個(gè)松弛參數(shù)用以描述不同物理量趨向平衡態(tài)的過程.除了少數(shù)幾個(gè)有著明確物理含義的松弛參數(shù)外,MRT 模型還有若干個(gè)可調(diào)參數(shù),相關(guān)文獻(xiàn)[34-35]在對MRT 模型進(jìn)行系統(tǒng)研究后指出,對這些參數(shù)進(jìn)行合理的調(diào)節(jié)可以有效改善數(shù)值穩(wěn)定性以及減小邊界條件的誤差.

MRT 模型對應(yīng)的速度離散后的碰撞項(xiàng)可以表示如下[36]

其中 Λ 代表物理松弛系數(shù)矩陣,其與網(wǎng)格分辨率無關(guān);而的表達(dá)式仍然使用式 (8).上式中使用了與D'Humieres 相同的記號,以“〈·| ”,“|·〉 ”來分別表示行向量以及列向量,即

這里上標(biāo)“T ”代表轉(zhuǎn)置操作符.將式 (28) 代入到式(5) 的矩陣形式中后,可以反解出由離散分布函數(shù)變量 |g〉 所表示的連續(xù)分布函數(shù)變量 |f〉 如下

其中I為單位矩陣;上標(biāo)“-1 ”代表矩陣求逆操作符.進(jìn)而可以得到MRT 模型下由離散分布函數(shù)變量|g〉表示的碰撞項(xiàng)

將式 (30) 代入到式 (6) 的矩陣形式中,可以得到

若記

則將上式代入到式 (31) 中后,便可以得到MRT 模型下關(guān)于離散分布函數(shù)變量 |g〉 的演化方程

其中,M是一個(gè)q×q的正交變換矩陣;M-1為M的逆矩陣,且滿足以下關(guān)系

式中,|m〉 代表轉(zhuǎn)換矩陣M下分布函數(shù) |g〉 對應(yīng)的一組矩.對式 (33) 兩邊同時(shí)左乘轉(zhuǎn)換矩陣M后便可得到MRT 模型下關(guān)于 |m〉 的演化方程[23]

這里 |ms(x,t)〉=M|S(x,t)〉.而轉(zhuǎn)換矩陣M與速度集的選取有關(guān),以D2Q9 速度集為例,其對應(yīng)的轉(zhuǎn)換矩陣M為[34]

其中,sν與運(yùn)動黏度系數(shù) ν,以及se與體積黏度系數(shù)ζ的關(guān)系可以由CE 分析給出[37]

與BGK 模型類似,MRT 模型下演化關(guān)系式(33) 的解 |g〉≠|(zhì)f〉,故也需要給出MRT 模型下,由離散分布函數(shù)變量 |g〉 計(jì)算宏觀量的表達(dá)式.通過引入Hermite 多項(xiàng)式矩陣H和Hermite 轉(zhuǎn)換矩陣MH,可以根據(jù)式 (29) 得到MRT 模型下,連續(xù)分布函數(shù)變量 |f〉 對應(yīng)的各階Hermite 矩 |a〉 與離散分布函數(shù)變量 |g〉 對應(yīng)的各階Hermite 矩 |ag〉 間的關(guān)系;以及連續(xù)非平衡態(tài)分布函數(shù)變量 |fneq〉 對應(yīng)的各階Hermite矩 |aneq〉 與離散非平衡態(tài)分布函數(shù)變量 |gneq〉 對應(yīng)的各階Hermite 矩間的關(guān)系如下(各變量的具體表達(dá)式以及詳細(xì)推導(dǎo)過程見附錄A)

進(jìn)而我們可以通過式 (13)、式(19)、式(39 a)以及附錄A 中MH得到MRT 模型下當(dāng)源項(xiàng)為體力項(xiàng)時(shí),由離散分布函數(shù)變量gα表示的密度 ρ、速度u

上述結(jié)果也與Chen 等[37]通過CE 分析得到結(jié)果保持一致.

當(dāng)初始流場中包含非平衡態(tài)時(shí),MRT 模型下的分布函數(shù)也需要正確的初始化.同樣可以使用Hermite 展開來初始化MRT 模型下的分布函數(shù),其中平衡態(tài)分布函數(shù)的初始化與BGK 模型一致,而離散非平衡態(tài)分布函數(shù)的初始化,則需要將附錄A 中 |gneq〉 的前6 個(gè)矩代入到式 (26),進(jìn)而可以得到MRT 模型下當(dāng)源項(xiàng)為體力項(xiàng)時(shí),初始流場中的具體表達(dá)式

上式中各參數(shù)的具體形式可以參見附錄A.

此外,由于矩陣M和矩陣初始化后就已確定,且其在式 (33) 演化過程中不再改變,故可以在初始化時(shí)將計(jì)算后單獨(dú)存儲起來以避免在執(zhí)行碰撞步驟時(shí)重復(fù)計(jì)算,從而節(jié)省計(jì)算量.且直接使用式 (33) 去執(zhí)行分布函數(shù)的碰撞步驟相比于使用式(34) 和式 (35) 更加節(jié)省時(shí)間,經(jīng)我們程序測試后發(fā)現(xiàn)CPU 消耗時(shí)間減少 25% 左右.

2 局部網(wǎng)格加密中的數(shù)據(jù)轉(zhuǎn)換

如前所述,同一點(diǎn)處不同網(wǎng)格分辨率下的離散分布函數(shù)變量gα之間的轉(zhuǎn)換是必要且不能忽略的.使用局部網(wǎng)格加密算法需要保證,當(dāng)粗細(xì)網(wǎng)格均使用同一離散速度集且在忽略時(shí)空離散誤差后,粗細(xì)網(wǎng)格區(qū)域在同一時(shí)空點(diǎn)下數(shù)值求解后得到的連續(xù)分布函數(shù)變量fα必須是一致的,而且不同分辨率下同一點(diǎn)處的物性參數(shù)如:物理松弛系數(shù) τ 以及運(yùn)動黏度系數(shù) ν 等也必須保持一致.

LBM 在對式 (1) 做時(shí)空離散時(shí)為了消除隱式,引入了與分辨率有關(guān)的離散分布函數(shù)變量gα以及相應(yīng)的格子松弛系數(shù) τν(BGK 模型)和(MRT 模型)

根據(jù)上式可知,不同網(wǎng)格分辨率下的同一時(shí)空點(diǎn)上,格子松弛系數(shù) τν和以及實(shí)際求解得到的離散分布函數(shù)變量gα和 |g〉 是不一致的,故粗細(xì)網(wǎng)格在交換這些量之前需要經(jīng)過合適的轉(zhuǎn)換.與之相反的是,物理松弛系數(shù) τ 以及 Λ-1是不隨網(wǎng)格分辨率變化的;而當(dāng)考慮粗細(xì)網(wǎng)格區(qū)域均使用同一離散速度集時(shí),則數(shù)值求解得到的連續(xù)分布函數(shù)變量fα和 |f〉,在不考慮時(shí)空離散誤差的前提下,也是不隨網(wǎng)格分辨率變化的.故接下來將以此為橋梁,建立起B(yǎng)GK模型和MRT 模型下,不同網(wǎng)格分辨率間的數(shù)據(jù)轉(zhuǎn)換關(guān)系.

下文中,與粗網(wǎng)格和細(xì)網(wǎng)格區(qū)域相關(guān)的任何量分別用上標(biāo)“c”和“f”表示.設(shè)局部加密時(shí)粗細(xì)網(wǎng)格的空間步長比例為

局部網(wǎng)格加密算法中不同網(wǎng)格分辨率下的區(qū)域仍然使用的是同一離散速度集,故在考慮對流尺度下,粗細(xì)網(wǎng)格的時(shí)間步長也有同樣的比例關(guān)系

故為了保證粗細(xì)網(wǎng)格上物理時(shí)間的統(tǒng)一推進(jìn),粗網(wǎng)格上運(yùn)行一步則細(xì)網(wǎng)格上需要運(yùn)行n步.故n常設(shè)為大于 1 的整數(shù).

2.1 BGK 模型下的局部網(wǎng)格加密

為了將我們的方法與現(xiàn)有推導(dǎo)過程更好地對比,在給出新的方法前,先梳理一下Liou 等[25]在考慮源項(xiàng)時(shí),分布函數(shù)轉(zhuǎn)換關(guān)系的推導(dǎo)過程.

首先對式 (17) 使用泰勒展開

考慮CE 展開如下

其中 ? 為與克努森數(shù)Kn同階的小量.將式 (47) 代入到式 (46) 中,便可以得到與 ?0,?1,?2同階的方程

而式 (48b) 等號左邊項(xiàng)均為宏觀量與離散速度ξα的函數(shù),故當(dāng)粗細(xì)網(wǎng)格區(qū)域使用同一離散速度集時(shí),其與網(wǎng)格分辨率無關(guān),故如下關(guān)系式成立

而格子松弛系數(shù) τν的轉(zhuǎn)換則是通過保證動力黏度系數(shù) μ=p(τν-0.5)?t一致的基礎(chǔ)上推導(dǎo)出來的,即

以上便是Liou 等在考慮源項(xiàng)存在時(shí),給出的BGK 模型下粗細(xì)網(wǎng)格間離散非平衡態(tài)分布函數(shù)變量以及格子松弛系數(shù) τν轉(zhuǎn)換關(guān)系的推導(dǎo)過程,而轉(zhuǎn)換關(guān)系的推導(dǎo)過程本質(zhì)上是利用CE 展開尋找與網(wǎng)格分辨率無關(guān)的量.不難發(fā)現(xiàn)整個(gè)推導(dǎo)過程是較為復(fù)雜的,而且使用了來近似對于恢復(fù)N-S 方程而言,這樣的近似是足夠的,但對于需要捕捉稀薄效應(yīng)的高階LBM 方法[38]而言,一階近似顯然是不夠的.下面將給出一種更加簡潔且直觀的推導(dǎo)方法.

如前所述,BGK 模型下的物理松弛系數(shù) τ 和 連續(xù)分布函數(shù)變量fα與網(wǎng)格分辨率無關(guān),故由式 (15)和式 (14)可知如下關(guān)系是成立的

對上式簡單整理后,便可以得到不同網(wǎng)格分辨率在重合點(diǎn)處的數(shù)據(jù)轉(zhuǎn)換關(guān)系如下

將式 (52b) 二邊減去平衡態(tài)并利用式 (52a),便可以得到與式 (50) 一致的轉(zhuǎn)換關(guān)系式.

相比之下,我們的推導(dǎo)過程中除忽略時(shí)空離散誤差外,并沒有做任何額外的假設(shè)(沒有對非平衡態(tài)分布函數(shù)做近似假設(shè)),故上述轉(zhuǎn)換關(guān)系同樣保證了粗細(xì)網(wǎng)格重合點(diǎn)處高階非平衡態(tài)分布函數(shù)的一致.然而,Liou 等在使用的假設(shè)后,也得到了與我們相同的結(jié)果.而這主要是由于各階非平衡態(tài)分布函數(shù)間的遞推關(guān)系導(dǎo)致的(詳細(xì)推導(dǎo)過程見附錄 B)

上述推導(dǎo)過程不僅解釋了我們的推導(dǎo)結(jié)果與前人結(jié)果一致的原因,同時(shí)也為之前學(xué)者對非平衡態(tài)分布函數(shù)做一階近似提供了理論依據(jù),而且上述推導(dǎo)過程也表明了局部網(wǎng)格加密算法同樣可以被擴(kuò)展到高階LBM 中.

2.2 MRT 模型下的局部網(wǎng)格加密

下面將繼續(xù)使用新方法給出MRT 模型下,不同網(wǎng)格分辨率在重合點(diǎn)處離散分布函數(shù)變量 |gc〉 和|gf〉 以及格子松弛系數(shù)矩陣間的轉(zhuǎn)換關(guān)系.與上一節(jié)BGK 模型下的推導(dǎo)過程類似,不同分辨率下的物理松弛系數(shù)矩陣 Λ-1應(yīng)保持一致,故由附錄A 可得

而由CE 分析可知,sq和s?并不會直接出現(xiàn)在N-S 方程中,故可以被視為自由參數(shù),但其對數(shù)值穩(wěn)定性及邊界條件誤差有一定影響[34].Chen 等[37]認(rèn)為這兩個(gè)參數(shù)在不同網(wǎng)格分辨率下可以不需要轉(zhuǎn)換,即但本文中為了保持轉(zhuǎn)換關(guān)系的一致性,仍然采用式 (58) 作為粗細(xì)網(wǎng)格間s?和sq的轉(zhuǎn)換關(guān)系式.

同樣,直接從式 (29) 出發(fā),基于連續(xù)分布函數(shù)變量 |f〉 在不同網(wǎng)格分辨率下的重合點(diǎn)處保持一致,便可以得到MRT 模型中粗細(xì)網(wǎng)格在時(shí)空重合點(diǎn)處的離散分布函數(shù)變量 |g〉 的轉(zhuǎn)換關(guān)系

且該結(jié)果與Huang 等[23]的結(jié)果也完全一致.

綜上,我們以保證不同網(wǎng)格分辨間連續(xù)分布函數(shù)變量fα,|f〉 以及物理松弛系數(shù) τ,Λ-1一致為基礎(chǔ),以更加直觀且簡潔的方法分別推導(dǎo)了BGK 模型和MRT 模型下,不同網(wǎng)格分辨率在時(shí)空重合點(diǎn)處離散分布函數(shù)變量gα,|g〉 和格子松弛系數(shù) τν,間的轉(zhuǎn)換關(guān)系,其相比于之前學(xué)者[23,25]的推導(dǎo)過程更加簡潔且直觀.

3 局部網(wǎng)格加密算法的數(shù)值驗(yàn)證

本節(jié)數(shù)值模擬中將采用與Yu 等[18]相同的網(wǎng)格排布方式,其示意圖如圖2 所示.同時(shí),為了避免局部網(wǎng)格加密算法中插值誤差對數(shù)值模擬的影響,我們使用與Gendre 等[39]相同的高階插值公式.而對于加密界面處缺失的離散分布函數(shù),采用Chen 等[40]提出的只使用空間插值的方法.

圖2 粗細(xì)網(wǎng)格排布示意圖.灰色區(qū)域?yàn)榇旨?xì)網(wǎng)格的重疊區(qū)域Fig.2 Schematic diagram of the coarse and fine grid arrangement.The grey area is the overlap between the coarse and fine grids

下面將通過使用局部網(wǎng)格加密技術(shù)的LBM來對強(qiáng)迫泰勒-格林渦流動、平板泊肅葉流中對流-擴(kuò)散問題、頂蓋驅(qū)動方腔流動和一維剪切波問題進(jìn)行模擬,在前兩個(gè)算例中均包含非均勻且非定常的源項(xiàng),故以此來驗(yàn)證粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系對復(fù)雜源項(xiàng)的適應(yīng)性;而有著較為復(fù)雜流動現(xiàn)象的頂蓋驅(qū)動方腔流動則可以較好地驗(yàn)證局部網(wǎng)格加密技術(shù)的有效性;一維剪切波問題的模擬則是用來檢驗(yàn)局部網(wǎng)格加密技術(shù)對原本的LBM 算法在耗散方面的影響.

3.1 強(qiáng)迫泰勒-格林渦流動

在標(biāo)準(zhǔn)的 Taylor-Green vortex 流動中沒有外力的參與,因此流動會隨時(shí)間單調(diào)衰減.為了驗(yàn)證源項(xiàng)存在時(shí),粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系的正確性,我們以Min 等[41]使用的二維強(qiáng)迫泰勒-格林渦流動作為驗(yàn)證算例.該流動中不僅有體力的參與,且該體力非均勻且非定常,故可以較為充分地驗(yàn)證BGK 模型以及MRT 模型下轉(zhuǎn)換關(guān)系的正確性.

該流動問題對應(yīng)的控制方程為

流動的物理區(qū)域被設(shè)為:(x,y)∈(0:Lx,0:Ly),4 周為周期邊界條件,考慮單位質(zhì)量受到的體力為b=k2ν(1-Q)u時(shí),則該流動對應(yīng)的理論解如下

式中,U0代表初始速度的幅值;kx=2π/Lx和ky=2π/Ly分別代表x和y方向上的波數(shù)P0代表參考壓力(可以為任意值).為了簡單,我們采用kx=ky,Lx=Ly=L=1 以及P0=1.此外,Q代表衰減因子,可以通過調(diào)節(jié)參數(shù)Q來控制流動的衰減速度.

(1)Q=1 時(shí),代表沒有外力,各項(xiàng)宏觀量按標(biāo)準(zhǔn)的泰勒-格林渦流動衰減;

(2)Q=0.5 時(shí),代表各項(xiàng)宏觀量的指數(shù)衰減率僅為體力不存在時(shí)的一半;

(3)Q=0 時(shí),體力輸入的能量與黏性耗散平衡,流場中各項(xiàng)宏觀量與初始流場保持一致且不隨時(shí)間變化;

(4)Q=-0.5 時(shí),體力輸入的能量大于黏性耗散,流場中的壓力和速度等宏觀量的絕對值將隨著時(shí)間的增加而增加.

使用t=0 時(shí)的解作為流場各宏觀量的初始值,并根據(jù)式 (8) 對平衡態(tài)分布函數(shù)進(jìn)行初始化;且由于該流動問題對應(yīng)的初始流場中包含非平衡態(tài)部分,故需要根據(jù)式 (27) 和式 (42) 分別對BGK 模型和MRT 模型下的離散非平衡態(tài)分布函數(shù)初始化;而兩者之和便可以作為初始流場對應(yīng)的離散分布函數(shù)變量gα.

粗細(xì)網(wǎng)格上均采用D2Q9 離散速度集.初始速度幅值U0對應(yīng)的Mach 數(shù)為 0.01;Reynolds 數(shù)Re=U0L/ν=100;網(wǎng)格數(shù):Nx=Ny=100.加密區(qū)域范圍為 (x/L,y/L)∈(0.4:0.6,0.1:0.4),加密比例為2.MRT 模型下各松弛系數(shù)如下

在細(xì)網(wǎng)格中,MRT 模型下的各松弛系數(shù)直接使用式 (57) 及 式 (58) 轉(zhuǎn)換而來.此外還需注意的是原來在計(jì)算速度u時(shí),有如下公式(考慮源項(xiàng)為體力項(xiàng))

但本算例中b=k2ν(1-Q)u,其大小與u直接相關(guān),故需要將b代入到式 (63) 中來顯式得到u的計(jì)算公式

為了充分驗(yàn)證粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系的正確性,首先分別繪制了BGK 模型和MRT 模型在Q=0.5 以及tU0/L=1 時(shí),相對壓力P-P0、水平方向速度ux以及黏性應(yīng)力 σxx的分布云圖.數(shù)值計(jì)算結(jié)果如圖3 所示,圖中的黑色方框?yàn)榫W(wǎng)格加密的界面,其所包圍的區(qū)域?yàn)榧用軈^(qū)域.從圖中可以看出,各宏觀量在加密界面兩側(cè)均光滑無間斷.

圖3 強(qiáng)迫泰勒-格林渦各宏觀量的分布云圖(Q=0.5, tU0/L=1,黑色矩形框內(nèi)的區(qū)域?yàn)榧用軈^(qū)域)Fig.3 The contours of each macroscopic quantity of the forced Taylor-Green vortex (Q=0.5, tU0/L=1,the area inside the black rectangular box is the refinement region)

此外,還分別計(jì)算了當(dāng)Q為 0.5,0.0 和 -0.5 時(shí),在 (x/L,y/L)=(0.4,0.25) 點(diǎn)(粗細(xì)網(wǎng)格重合點(diǎn))處,相對壓力P-P0、水平方向速度ux以及黏性應(yīng)力σxx隨時(shí)間的變化,其數(shù)值結(jié)果如圖4 所示.在每張子圖中,繪制了在Q取不同的值時(shí),BGK 模型以及MRT 模型下數(shù)值結(jié)果與理論解果的對比.其中實(shí)線代表由式 (61) 給出的理論解,圓圈和倒三角分別代表BGK 模型和MRT 模型下的數(shù)值結(jié)果,而不同的顏色代表Q取不同的值.觀察圖4 后發(fā)現(xiàn),兩種碰撞模型的數(shù)值結(jié)果都與理論解吻合得非常好.

圖4 強(qiáng)迫泰勒-格林渦中在 (x/L,y/L)=(0.4,0.25) 點(diǎn)處各宏觀量隨時(shí)間的變化.(a)～(c)分別為BGK 模型和MRT 模型下的相對壓力P-P0、水平方向速度 ux 和黏性應(yīng)力 σxx 與理論結(jié)果的對比Fig.4 The variation of each macroscopic quantity with time at the point (x/L,y/L)=(0.4,0.25) in the forced Taylor-Green vortex.(a)～(c) are comparisons between the relative pressure P-P0,velocity in the x direction ux and viscous stress σxx under BGK model and MRT model,respectively,and the theoretical results

最后,為了進(jìn)一步觀察粗細(xì)網(wǎng)格界面兩側(cè)宏觀量的連續(xù)性,我們還繪制了當(dāng)Q為 0.5,0.0 和-0.5以及tU0/L=1 時(shí),在y/L=0.25 處BGK 模型和MRT模型下各宏觀量隨x/L變化的剖線圖,同時(shí)也與理論解進(jìn)行了對比.其數(shù)值結(jié)果如圖5 所示,圖中兩條黑色虛直線間的區(qū)域?yàn)榧用軈^(qū)域.從圖中可以看出,兩種碰撞模型都取得了非常好的數(shù)值結(jié)果,且各宏觀量在加密界面兩側(cè)光滑連續(xù).

圖5 強(qiáng)迫泰勒-格林渦中各宏觀量的剖線圖(tU0/L=1, y/L=0.25,兩條黑色虛線間的區(qū)域?yàn)榧用軈^(qū)域).(a)～(c)分別為BGK 模型和MRT 模型下相對壓力 P-P0、水平方向速度 ux 和黏性應(yīng)力 σxx 與理論結(jié)果的對比Fig.5 Profile diagram of each macroscopic quantity of forced Taylor-Green vortex (tU0/L=1, y/L=0.25,the area between the two black dashed lines is refined).(a)～(c) are comparisons between the relative pressure P-P0,velocity in the x direction ux and viscous stress σxx under BGK model and MRT model,respectively,and the theoretical results

BGK 模型和MRT 模型下良好的數(shù)值結(jié)果也印證了第 2 節(jié)中粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系的正確性.接下來將以平板泊肅葉流中對流-擴(kuò)散問題為驗(yàn)證算例,進(jìn)一步測試BGK 模型下粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系對復(fù)雜源項(xiàng)的適應(yīng)性.

3.2 平板泊肅葉流中對流-擴(kuò)散問題

如前所述,合理設(shè)計(jì)LBM 中的源項(xiàng)S可以極大地?cái)U(kuò)展該方法的適用范圍,許多學(xué)者通過CE 分析反設(shè)計(jì)源項(xiàng)S的形式,以使得修改后的LBM 可以滿足特定要求[11].甚至在很多情況下,LBM 被當(dāng)作一款求解一些特殊偏微分方程的高效通用求解器[42].其中,對流-擴(kuò)散方程作為一種可以描述由于對流和擴(kuò)散過程而引起的傳遞熱量、質(zhì)量或其他物理量輸運(yùn)現(xiàn)象[43]的特殊偏微分方程,被廣泛應(yīng)用于求解各種對流-擴(kuò)散問題.該方程可以被表示如下

式中,?(x,t) 代表一個(gè)標(biāo)量,u(x,t) 是由N-S 方程控制的速度,D為擴(kuò)散系數(shù).

對于如何使用LBM 求解對流-擴(kuò)散方程,之前的學(xué)者們利用CE 分析給出了許多不同的方案[44-48].這里使用Chai 等[48]提出的通過CE 分析反設(shè)計(jì)源項(xiàng)進(jìn)而恢復(fù)對流-擴(kuò)散方程的方案,該方法的標(biāo)量場中分布函數(shù)的演化方程如下

這里 ?α對應(yīng)標(biāo)量場中的分布函數(shù),而Rα則是為了在恢復(fù)式 (65) 時(shí)減小額外計(jì)算誤差所添加的源項(xiàng),且和Rα分別表示如下

這里 ?? 可以直接使用如下公式計(jì)算

接下來將用上述方案去求解平板泊肅葉流中對流-擴(kuò)散問題.在計(jì)算這個(gè)問題時(shí),需要兩套分布函數(shù):gα用于計(jì)算速度場中的密度 ρ 以及速度u;?α用于計(jì)算標(biāo)量場中的 ?.

該問題對應(yīng)的計(jì)算簡圖如圖6 所示.其中上下板處的灰色區(qū)域?yàn)榧用軈^(qū)域,加密范圍為:y/H∈{(0.0:0.1)∪(0.9:1.0)},加密比例為 2.其他參數(shù)為:定常時(shí)中心速度U0對應(yīng)的Mach 數(shù)為 0.01;下板標(biāo)量:?b=0.0,上板標(biāo)量:?t=0.1;Reynolds 數(shù)Re=HU0/ν=10,H為槽道寬度;Peclet 數(shù)Pe=HU0/D=10;網(wǎng)格數(shù)Nx=6,Ny=20,粗細(xì)網(wǎng)格均采用D2Q9 速度集以及BGK 碰撞模型.該流動采用水平方向的體力驅(qū)動,單位質(zhì)量所受的體力為:bx=8νU0/H2.左右為周期邊界條件,上下邊界均使用非平衡態(tài)外推邊界條件[49],其中邊界點(diǎn)處密度的計(jì)算使用非平衡態(tài)反彈邊界條件[50]中密度的計(jì)算方法.

圖6 平板泊肅葉流中對流-擴(kuò)散問題示意圖Fig.6 Schematic diagram of diffusion problem in a planar Poiseuille flow

該流動問題中速度ux(y,t) 所滿足的控制方程為

對應(yīng)的解析解可以通過分離變量法求得,其級數(shù)形式可以表示為

式中k=2n+1.此外,標(biāo)量 ? (y,t) 所滿足的控制方程可以表示如下

其對應(yīng)的解析解也同樣可以通過分離變量法給出,且級數(shù)形式表示如下

其中J=-D??+?u代表擴(kuò)散通量與對流通量之和,而Jx,Jy分別為x,y方向上的分量.

我們使用局部網(wǎng)格加密技術(shù),分別計(jì)算了在t?=tν/H2=0.04,0.1,1.0 時(shí),水平方向的速度ux、標(biāo)量 ?、通量Jx和Jy并與通過式 (71)和式(73) 得到的理論結(jié)果進(jìn)行了對比.數(shù)值結(jié)果如圖7 所示,圖中的黑色虛線代表粗細(xì)網(wǎng)格的界面,實(shí)線代表理論解,圓圈代表數(shù)值解,不同的顏色代表不同的時(shí)刻.從圖中可以看出,不同時(shí)刻下的數(shù)值結(jié)果都與理論解吻合得非常好,而且在加密界面的兩側(cè)各宏觀量均光滑且連續(xù).良好的數(shù)值結(jié)果也進(jìn)一步驗(yàn)證了第 2節(jié)中粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系的正確性,同時(shí)也驗(yàn)證了該轉(zhuǎn)換關(guān)系對復(fù)雜源項(xiàng)的良好適應(yīng)性.

圖7 平板泊肅葉流中對流-擴(kuò)散問題(黑色虛線為加密界面).(a)～(d)分別為BGK 模型下水平方向速度 ux、標(biāo)量 ?、通量 Jx 和 Jy 與理論結(jié)果的對比ux,scalar ?,flux Jy and Jx compared with theoretical results under the BGK modelFig.7 Diffusion problem in a planar Poiseuille flow (the black dashed lines are the refinement interface).(a)～(d) represent the velocity in x direction

3.3 頂蓋驅(qū)動方腔流動

頂蓋驅(qū)動方腔流作為經(jīng)典的基準(zhǔn)算例已被廣泛用作數(shù)值方法的測試案例.為了對數(shù)值結(jié)果進(jìn)行合理的評價(jià),使用Tamer 等[51]的數(shù)值結(jié)果作為參考數(shù)據(jù)并與之進(jìn)行對比.使用LBM 算法去計(jì)算頂蓋驅(qū)動方腔流時(shí),如果使用的邊界條件較為粗糙且網(wǎng)格分辨率較低時(shí),頂蓋左右兩個(gè)角點(diǎn)處往往會產(chǎn)生明顯的壓力“噪聲”.

下面使用局部加密算法對該流動進(jìn)行數(shù)值模擬,用以進(jìn)一步展示局部加密算法的優(yōu)勢.本次模擬使用BGK 模型以及D2Q9 速度集,雷諾數(shù)Re=UwL/ν=1000,其中L代表頂蓋長度,Uw為頂蓋速度其對應(yīng)的馬赫數(shù)為 0.1,ν 為運(yùn)動黏度系數(shù).邊界條件:左右壁面以及下壁面均使用修正反彈格式,頂蓋處使用運(yùn)動邊界對應(yīng)的反彈格式[52],即

圖8 頂蓋驅(qū)動方腔流4 個(gè)角點(diǎn)示意圖,其中黑色粗實(shí)線代表壁面Fig.8 Diagram of four corner points of the square cavity flow driven by the top lid,where the thick black line represents the wall surface

其中 (xb,yb) 代表右上角點(diǎn)對應(yīng)的坐標(biāo),而 2,3,6 方向均使用靜止邊界的反彈格式.左上角點(diǎn)也是同樣處理.

對整個(gè)頂蓋附近實(shí)施加密,加密范圍為:y/L∈(0.7:1.0),加密比例為 2.不同分辨率下壓力分布云圖如圖9 所示.其中,不僅對比了使用局部加密(UCG-L)和使用均勻粗網(wǎng)格(分辨率與局部加密中的粗網(wǎng)格分辨率一致,UCG)以及均勻細(xì)網(wǎng)格(分辨率與局部加密中的細(xì)網(wǎng)格分辨率一致,UFG)下的結(jié)果,同時(shí)還添加了與UCG-L 達(dá)到穩(wěn)態(tài)時(shí)理論計(jì)算耗時(shí)相同的均勻網(wǎng)格(EQ)下的結(jié)果,用來對比使用局部網(wǎng)格加密算法的非均勻網(wǎng)格與相同計(jì)算量下的均勻網(wǎng)格兩者的數(shù)值結(jié)果.從圖9 中不難看出,采用局部加密后原來的“噪聲”得到了極大的消除,而且其對“噪聲”的抑制效果比相同計(jì)算耗時(shí)下的均勻網(wǎng)格更好.這也側(cè)面驗(yàn)證了局部網(wǎng)格加密方法在消除局部“噪聲”的有效性.值得一提的是,Dong 等[35]通過詳細(xì)的分析指出,這里的壓力噪聲主要是邊界條件隱藏誤差導(dǎo)致的,且可以通過利用雙松弛時(shí)間碰撞模型優(yōu)化加以有效消除,其文章中的數(shù)值結(jié)果也證明了這一點(diǎn).

圖9 Re=1000 時(shí)頂蓋驅(qū)動方腔流的壓力分布云圖對比,(a)～(d)分別為UCG,UFG,EQ 和UCG-L 對應(yīng)的壓力分布云圖 ((d)中的黑色直線為加密界面,界面以上為加密區(qū)域)Fig.9 When Re=1000,the contours of pressure of the lid-driven cavity flow are compared.(a)～(d) are the contours of pressure corresponding to UCG,UFG,EQ and UCG-L,respectively (the black line in (d) is the grid refinement interface,and the area above the black line is the refinement area)

圖10 中繪制了使用不同分辨率時(shí)沿垂直方向穿過幾何中心的水平速度ux隨y/L的變化關(guān)系以及沿水平方向穿過幾何中心的垂直速度uy隨x/L的變化關(guān)系.而且從圖10(b) 以及圖10(d) 中可以很清楚地看到即使是遠(yuǎn)離加密區(qū)域UCG-L 的數(shù)值結(jié)果仍要好于與其所需計(jì)算量相同的均勻網(wǎng)格下的結(jié)果,而且在使用更少計(jì)算資源的情況下還取得了與UFG 非常相近的數(shù)值結(jié)果,這就表明在局部流域執(zhí)行網(wǎng)格加密可以很好地減小局部數(shù)值誤差對全流域的影響.

圖10 Re=1000 時(shí)頂蓋驅(qū)動方腔流中無量綱速度的剖線圖,UCG,UFG,EQ 和UCG-L 與參考數(shù)據(jù)的對比.(a) 沿垂直直線穿過幾何中心的ux/Uw,(b)為其局部放大,(c) 沿水平直線穿過幾何中心的 uy/Uw,(d)為其局部放大Fig.10 Profile diagram of dimensionless velocity in a lid-driven cavity flow when Re=1000,comparison of UCG,UFG,EQ and UCG-L with reference data.(a) ux/Uw in a vertical straight line through the center of the geometry,(b) local amplification of (a),(c) uy/Uw in a horizontal straight line through the center of the geometry,(d) local amplification of (c)

在表1 中,我們比較了3 種情況下每 1000 個(gè)迭代步的 CPU 實(shí)際消耗時(shí)間.如果假設(shè)在 UCG 上運(yùn)行到穩(wěn)態(tài)所需的 CPU 時(shí)間是t,若考慮對流時(shí)間尺度以及二維情況下則UFG 下運(yùn)行到相同狀態(tài)下,理論上所需的 CPU 時(shí)間是 8t,而上述加密的UCG-L只需要 3.4t,表1 中的測量結(jié)果也基本支持上述論斷.值得注意的是,在本算例中UCG-L 所需的計(jì)算時(shí)間還不及 UFG 的一半,但其數(shù)值結(jié)果卻與 UFG非常相近.有理由相信在某些特殊算例中通過合理布置加密區(qū)域,局部網(wǎng)格加密可以實(shí)現(xiàn)更高的加速比.

表1 不同分辨率下程序演化1000 步時(shí)CPU 消耗時(shí)間對比Table 1 Comparison of CPU consumption time when the program evolves 1000 steps at different resolutions

3.4 一維剪切波問題

為了觀察局部網(wǎng)格加密技術(shù)對LBM 算法在數(shù)值耗散方面的影響,我們選擇了一維剪切波問題[53]作為研究算例,其示意圖如圖11 所示.該問題對應(yīng)的控制方程及初始化流場如下

圖11 一維剪切波問題示意圖Fig.11 Schematic diagram of one-dimensional shear wave problem

式中,u0代表初速流場中速度的幅值.由于物理黏性的存在,速度ux會隨著時(shí)間逐漸衰減,其對應(yīng)的理論解為

針對該問題數(shù)值模擬的各參數(shù)和設(shè)置如下:上下使用周期邊界條件;粗細(xì)網(wǎng)格上均采用BGK 碰撞模型和D2Q9 速度集;非平衡態(tài)分布函數(shù)的初始化仍然使用式 (27).初始速度幅值u0對應(yīng)的Mach數(shù)為 0.1;Reynolds 數(shù)Re=Hu0/ν=10,H為槽道寬度;網(wǎng)格數(shù)Nx=4,Ny=20.分別計(jì)算了加密區(qū)域范圍為:(x/L,y/L)∈(0.0:1.0,0.1:0.25),(x/L,y/L)∈(0.0:1.0,0.1:0.4),(x/L,y/L)∈(0.0:1.0,0.1:0.6) 和(x/L,y/L)∈(0.0:1.0,0.4:0.6),且加密比例均為2 的情況下,UCG-L 和UCG 下的數(shù)值耗散 νN與物理黏性ν之比,其中數(shù)值黏性可以通過反解式 (78) 求得

分別繪制出UCG-L 和UCG 下的 νN/ν 隨著y/H變化的剖線圖,其數(shù)值結(jié)果如圖12 所示.圖12 展示了在不同區(qū)域使用局部網(wǎng)格加密后,在給定時(shí)刻(t?=tν/H2=4π2/ln2) 下的數(shù)值黏性與物理黏性之比(紅色線)以及均勻粗網(wǎng)格下的數(shù)值黏性與物理黏性之比(藍(lán)色線),圖中的黑色虛線代表粗細(xì)網(wǎng)格的分界線.從圖中可以發(fā)現(xiàn),不同加密區(qū)域下,加密界面兩側(cè)的數(shù)值黏性相比于均勻網(wǎng)格下,表現(xiàn)出了不一樣的結(jié)果,既有增大的情況也有降低的情況.而其中數(shù)值黏性小于物理黏性現(xiàn)象的出現(xiàn)可能與加密界面附近的非物理震蕩[19,54-55]有著密切的關(guān)系,故還需要一系列更加系統(tǒng)的理論研究將局部網(wǎng)格加密技術(shù)引起的數(shù)值耗散和數(shù)值色散與加密區(qū)域附近的非物理震蕩聯(lián)系起來,而這也將作為我們后續(xù)工作的重點(diǎn).

圖12 一維剪切波問題在不同加密區(qū)域下均勻網(wǎng)格和非均勻網(wǎng)格對應(yīng)的數(shù)值黏性和物理黏性之比的分布情況(黑色虛線為加密界面).(a)～(d)分別為不同加密區(qū)域下UCG-L 和UCG 對應(yīng)的數(shù)值黏性和物理黏性之比隨 y/H 變化的剖線圖Fig.12 For one-dimensional shear wave problem,the distribution of the ratio of numerical viscosity and physical viscosity corresponding to uniform and non-uniform grids in different refinement regions (the black dashed line is the refinement interface).(a)～(d) are the cross-section plots of the ratio of numerical viscosity and physical viscosity corresponding to UCG-L and UCG with y/H in different refinement regions,respectively

圖12 一維剪切波問題在不同加密區(qū)域下均勻網(wǎng)格和非均勻網(wǎng)格對應(yīng)的數(shù)值黏性和物理黏性之比的分布情況(黑色虛線為加密界面).(a)～(d)分別為不同加密區(qū)域下UCG-L 和UCG 對應(yīng)的數(shù)值黏性和物理黏性之比隨 y/H 變化的剖線圖 (續(xù))Fig.12 For one-dimensional shear wave problem,the distribution of the ratio of numerical viscosity and physical viscosity corresponding to uniform and non-uniform grids in different refinement regions (the black dashed line is the refinement interface).(a)～(d) are the cross-section plots of the ratio of numerical viscosity and physical viscosity corresponding to UCG-L and UCG with y/H in different refinement regions,respectively (continued)

4 總結(jié)

本文在考慮任意源項(xiàng)存在的前提下,提出了一套規(guī)范且簡潔的粗細(xì)網(wǎng)格在重合點(diǎn)處的數(shù)據(jù)轉(zhuǎn)換關(guān)系推導(dǎo)方法,該方法既可以用于BGK 模型也可以適用于MRT 模型,同時(shí)也容易被推廣到其他碰撞模型下的數(shù)據(jù)轉(zhuǎn)換.該推導(dǎo)方法從保證粗細(xì)網(wǎng)格重合點(diǎn)處對應(yīng)的連續(xù)分布函數(shù)一致出發(fā),除忽略離散誤差外,并不依賴于CE 展開以及非平衡態(tài)近似.這也說明,局部網(wǎng)格加密算法不僅適用于連續(xù)流,而且也適用于高階非平衡態(tài)流動問題.此外,我們還從理論上證明了保證粗細(xì)網(wǎng)格間非平衡態(tài)部分的一階CE 近似一致,便可以保證整個(gè)非平衡態(tài)部分的一致,這也給之前學(xué)者對非平衡態(tài)部分做一階CE 近似提供了理論依據(jù).同時(shí),我們還通過引入新的Hermite 轉(zhuǎn)換矩陣MH直接構(gòu)建了MRT 模型下,時(shí)空離散前后分布函數(shù)對應(yīng)的各階Hermite 矩間的關(guān)系,從而將BGK 模型下基于Hermite 展開的初始化方法推廣到了MRT 模型下.

為了驗(yàn)證包含復(fù)雜源項(xiàng)的粗細(xì)網(wǎng)格間數(shù)據(jù)轉(zhuǎn)換關(guān)系,考慮了兩個(gè)具體算例.第1 個(gè)算例是外力作用下的不可壓強(qiáng)迫泰勒-格林渦流動,這個(gè)問題涉及非均勻非定常流場,并且由源項(xiàng)表示的外力也為非均勻場.初始流場對應(yīng)的分布函數(shù)采用Hermite 展開來準(zhǔn)確構(gòu)建.流場內(nèi)部的部分區(qū)域采用細(xì)網(wǎng)格.計(jì)算結(jié)果顯示粗細(xì)網(wǎng)格邊界處壓力、速度、應(yīng)力場均光滑銜接,這些宏觀量隨時(shí)間和空間的演化剖面都與解析解完全吻合,并且BGK 碰撞模型和MRT 碰撞模型下的計(jì)算結(jié)果也基本完全一致,證明了我們推導(dǎo)的基于兩種碰撞模型轉(zhuǎn)換關(guān)系的正確性.第2 個(gè)算例是平板泊肅葉流驅(qū)動的非定常對流-擴(kuò)散問題,這里標(biāo)量場對應(yīng)的Boltzmann 方程需要引進(jìn)源項(xiàng),才能準(zhǔn)確恢復(fù)宏觀非定常對流-擴(kuò)散方程.在上下兩個(gè)壁面附近,我們采用細(xì)網(wǎng)格局部加密.本問題的標(biāo)量場在流向有對流通量,在垂直壁面方向有擴(kuò)散通量,并且它們分布都不均勻,且隨時(shí)間演化.通過分離變量法,給出了標(biāo)量場的非定常解析解.計(jì)算結(jié)果顯示速度場、標(biāo)量場和標(biāo)量通量的兩個(gè)分量在粗細(xì)網(wǎng)格界面都光滑過渡,并與解析解完全吻合,進(jìn)一步驗(yàn)證了一般源項(xiàng)下轉(zhuǎn)換關(guān)系的正確性.此外,我們還使用局部網(wǎng)格加密技術(shù)對經(jīng)典的頂蓋驅(qū)動方腔流動問題進(jìn)行了數(shù)值模擬,數(shù)值結(jié)果表明局部網(wǎng)格加密技術(shù)在處理局部壓力噪聲方面有著明顯的優(yōu)勢,而且通過對比使用局部加密的非均勻網(wǎng)格和不使用局部加密的均勻網(wǎng)格的計(jì)算耗時(shí)以及宏觀速度的剖線圖,可以較為清楚地看到局部網(wǎng)格加密技術(shù)在處理復(fù)雜流動問題時(shí)的優(yōu)勢.最后,為了觀察局部網(wǎng)格加密技術(shù)對LBM 算法在數(shù)值耗散方面的影響,我們通過在一維剪切波中的不同區(qū)域使用局部網(wǎng)格加密,并將其數(shù)值黏性與均勻粗網(wǎng)格下的數(shù)值黏性進(jìn)行了對比,發(fā)現(xiàn)由局部加密引起的數(shù)值黏性與加密區(qū)域的選取有很大的關(guān)系,且在部分加密區(qū)域下還出現(xiàn)了數(shù)值黏性小于物理黏性的情況,一系列更加精細(xì)的算例和理論分析需要被考慮來闡述該現(xiàn)象的發(fā)生.

在上述的數(shù)值驗(yàn)證中,我們的算例均為二維算例,并沒有涉及更加復(fù)雜的三維流動問題,這主要是由于LBM 中的局部網(wǎng)格加密技術(shù)經(jīng)過20 多年的發(fā)展已經(jīng)趨于成熟,相關(guān)的文獻(xiàn)中[25,37,56]均已報(bào)道了大量使用局部網(wǎng)格加密技術(shù)計(jì)算復(fù)雜三維流動問題的測試結(jié)果,相關(guān)的數(shù)值結(jié)果均展示了局部網(wǎng)格加密技術(shù)在三維復(fù)雜流動問題中的優(yōu)勢及適應(yīng)性.本文的主要創(chuàng)新點(diǎn)在于從粗細(xì)網(wǎng)格間分布函數(shù)需要轉(zhuǎn)換的本質(zhì)出發(fā),構(gòu)建了一套直觀且簡潔的包含任意源項(xiàng)下分布函數(shù)轉(zhuǎn)換關(guān)系的推導(dǎo)過程,且我們的推導(dǎo)結(jié)果與前人的結(jié)果完全一致,這一點(diǎn)在 2.1 中也給出了解釋.而如上所述,LBM 中的局部網(wǎng)格加密技術(shù)已經(jīng)被廣泛應(yīng)用到各種復(fù)雜流動問題的計(jì)算中,且其使用的分布函數(shù)轉(zhuǎn)換關(guān)系與我們的也完全一致,故本文中并沒有對更加復(fù)雜的三維流動進(jìn)行模擬,而是借用二維情況下的強(qiáng)迫泰勒-格林渦流動和平板泊肅葉流中的對流-擴(kuò)散問題,著重考察了分布函數(shù)轉(zhuǎn)換關(guān)系對復(fù)雜源項(xiàng)的適應(yīng)性,以及通過對頂蓋驅(qū)動方腔流動的模擬初步展示局部網(wǎng)格加密技術(shù)在處理復(fù)雜流動問題方面的優(yōu)勢.

同時(shí)需要注意的是,粗細(xì)網(wǎng)格重合點(diǎn)處的分布函數(shù)轉(zhuǎn)換關(guān)系式 (48)和式(49) 中的粗細(xì)網(wǎng)格分辨率之比n理論上可以采用任意正整數(shù),但本文粗細(xì)網(wǎng)格間分布函數(shù)的轉(zhuǎn)換關(guān)系是在忽略時(shí)空離散誤差的前提下得出的,而實(shí)際情況下由于離散誤差的存在,不同網(wǎng)格分辨率下得到的數(shù)值結(jié)果中所包含的離散誤差也是不一樣的,故粗細(xì)網(wǎng)格界面處的數(shù)值結(jié)果間往往會產(chǎn)生一定的“數(shù)值間斷”,且該“數(shù)值間斷”隨著加密比例的增加而增加,對于一些較為復(fù)雜的流動,這往往會在粗細(xì)網(wǎng)格界面產(chǎn)生非物理的震蕩,故上面的算例中局部網(wǎng)格的加密倍數(shù)均為 2,而且絕大多數(shù)文獻(xiàn)的網(wǎng)格加密比例在計(jì)算中往往也被限定為 2.

必須指出的是,本文的分析并沒有考慮粗細(xì)網(wǎng)格間不同數(shù)值離散誤差有可能在網(wǎng)格加密界面處產(chǎn)生的數(shù)值震蕩.這個(gè)問題在已有文獻(xiàn)中有一些討論[19,57],但需要進(jìn)一步做理論分析并找到消除數(shù)值震蕩的有效方法.其中,不同碰撞模型對粗細(xì)網(wǎng)格界面的數(shù)值穩(wěn)定性有著不同的影響[19],但也需要從理論上進(jìn)一步做分析.如果網(wǎng)格加密邊界與物理壁面相交,那么物理壁面的反彈格式局部誤差[35]與粗細(xì)網(wǎng)格界面誤差交匯,此時(shí)往往會在壁面附近產(chǎn)生一定的數(shù)值震蕩[54],該問題的理論分析會更具挑戰(zhàn)性,但對實(shí)際應(yīng)用問題卻非常重要.這些都是需要進(jìn)一步研究的問題.

此外,測試和分析局部網(wǎng)格加密算法在數(shù)值色散、數(shù)值耗散方面的表現(xiàn)對于該技術(shù)的實(shí)際應(yīng)用有著非常重要的意義,尤其是有助于加密界面處非物理震蕩問題[19,54-55]的解決.在LBM 的局部網(wǎng)格加密算法中,粗細(xì)網(wǎng)格上仍然使用標(biāo)準(zhǔn)的LBM 算法,而不同的是需要在粗細(xì)網(wǎng)格交界面處進(jìn)行合理的數(shù)據(jù)轉(zhuǎn)換.其中,數(shù)據(jù)轉(zhuǎn)換時(shí)往往涉及未知分布函數(shù)的構(gòu)建,目前構(gòu)建方式可以分為有限差分方法和有限體積方法兩類,而不同的構(gòu)建方式對局部網(wǎng)格加密算法的數(shù)值耗散和數(shù)值色散都有不同的影響.故研究局部網(wǎng)格加密算法對LBM 在數(shù)值色散和數(shù)值耗散方面的影響,不僅需要詳細(xì)對比不同的構(gòu)建格式在數(shù)值耗散和數(shù)值色散方面的表現(xiàn),還需要從理論上分析不同構(gòu)建格式間結(jié)果差異的原因,以及數(shù)值耗散和數(shù)值色散對加密界面處非物理震蕩的影響,而這也將是我們下一步的工作重點(diǎn).考慮到該問題的復(fù)雜性,故我們將在后續(xù)文章中嘗試對該問題做較為完整和詳細(xì)的介紹.

附錄A MRT 模型下時(shí)空離散前后的分布函數(shù)對應(yīng)的各階Hermite 間的關(guān)系

首先,需要對式(29)兩邊同時(shí)左乘Hermite多項(xiàng)式矩陣H

為了讓上式可以直接反映連續(xù)分布函數(shù)變量 |f〉 的各階Hermite 矩與離散分布函數(shù)變量 |g〉 各階Hermite 矩間的關(guān)系,上式中的Hermite 多項(xiàng)式矩陣H可以直接由Hermite 多項(xiàng)式組成,以D2Q9 速度集為例

上式中各階Hermite 多項(xiàng)式的具體形式可以表示如下[58]

此時(shí),式 (A1) 中的H|f〉 可以表示如下

進(jìn)而式 (A1)可以表示為

上式直接給出了連續(xù)分布函數(shù)變量 |f〉 和離散分布函數(shù)變量 |g〉 對應(yīng)的Hermite 矩間的關(guān)系.在對式 (32) 兩邊取逆后,可以得到

若記Hermite 轉(zhuǎn)換矩陣MH為如下形式

結(jié)合式 (A7) 與式 (32) 后,式 (A6) 可以被進(jìn)一步簡化如下

D2Q9 速度集下的Hermite 轉(zhuǎn)換矩陣MH可以顯式表示為

上述Hermite 轉(zhuǎn)換矩陣MH中各參數(shù)表示的具體表達(dá)式如下

當(dāng)考慮源項(xiàng)為體力項(xiàng)時(shí),|gneq〉 的前6 個(gè)矩的具體形式如下

此外,對比式 (24) 和式 (A16) 后可以發(fā)現(xiàn):當(dāng)MRT 模型中與體積黏度系數(shù)有關(guān)的格子松弛參數(shù)se被設(shè)為se=sν時(shí),MRT 模型中的密度 ρ、速度u、黏性應(yīng)力 σ 以及初始化離散分布函數(shù)gα(Hermite 展開到兩階) 的計(jì)算表達(dá)式與BGK 模型下完全一致.

附錄B 各階連續(xù)非平衡態(tài)分布函數(shù)的遞推關(guān)系

直接從BGK 模型下的DVBE 出發(fā)

進(jìn)行簡單移項(xiàng)后

對式 (B2)使用若干次麥克斯韋迭代

上式也可以寫作如下形式

歸納可得如下關(guān)系