三維不可壓縮流的12速多松弛格子Boltzmann模型*

胡嘉懿 張文歡? 柴振華 施保昌 汪一航

1) (寧波大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,浙江 315211)

2) (華中科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,湖北 430074)

為提高多松弛(MRT)格子Boltzmann模型的計(jì)算效率,運(yùn)用反演法提出了一個(gè)求解三維不可壓縮流的12速M(fèi)RT格子Boltzmann模型(iD3Q12 MRT模型).這個(gè)模型比通常使用的D3Q13 MRT模型具有更高的計(jì)算效率.在數(shù)值模擬部分我們把iD3Q12 MRT模型與可壓縮性較小的一個(gè)13速多松弛模型(He-Luo D3Q13 MRT模型)在精確性和穩(wěn)定性方面作比較.通過模擬不同的流動(dòng),包括壓力驅(qū)動(dòng)的穩(wěn)態(tài)泊肅葉流、周期變化的壓力驅(qū)動(dòng)的非穩(wěn)態(tài)脈動(dòng)流、頂蓋驅(qū)動(dòng)的方腔流,可以發(fā)現(xiàn)iD3Q12 MRT模型模擬以上三種流動(dòng)時(shí)得到的數(shù)值解與解析解或與已有的結(jié)果符合很好,這說明我們提出的iD3Q12 MRT模型是準(zhǔn)確的.在模擬穩(wěn)態(tài)的泊肅葉流時(shí),兩個(gè)模型計(jì)算的速度場(chǎng)的全局相對(duì)誤差完全相同,且兩個(gè)模型都具有二階的空間精度.在模擬非穩(wěn)態(tài)脈動(dòng)流時(shí),大多情況下是12速模型的計(jì)算誤差更小,但在脈動(dòng)流的最大壓降增大時(shí),iD3Q12 MRT模型先發(fā)散,這說明He-Luo D3Q13 MRT模型具有更好的穩(wěn)定性.在模擬不同雷諾數(shù)下的頂蓋驅(qū)動(dòng)的方腔流時(shí),He-Luo D3Q13 MRT模型也比iD3Q12 MRT模型更穩(wěn)定.

1 引 言

格子Boltzmann方法(LBM)起源于20世紀(jì)70年代提出的格子氣自動(dòng)機(jī)方法(LGA),它克服了LGA方法的一些缺陷,例如消除了統(tǒng)計(jì)噪聲且其對(duì)應(yīng)的宏觀方程滿足伽利略不變性[1—3].并且LBM自身具有良好的計(jì)算局部性、程序的簡(jiǎn)潔性和拓展性等優(yōu)點(diǎn).另一方面從理論上,LBM方法可以從連續(xù)Boltzmann方程得到[4—6].因此,LBM作為計(jì)算流體動(dòng)力學(xué)(CFD)中一種有效的介觀數(shù)值模擬方法受到了廣泛的關(guān)注.它被應(yīng)用在模擬一些復(fù)雜的流體,例如：多相流[7—9]、懸浮液[10]、磁流體[11,12]、多孔介質(zhì)中的流體[13],還被用于求解一些偏微分方程,例如Burges方程[14]、對(duì)流擴(kuò)散方程[15,16]、泊松方程[17]、分?jǐn)?shù)階擴(kuò)散方程[18].

在研究人員的努力和發(fā)展下,LBM出現(xiàn)了lattice Bhatnagar-Gross-Krook (LBGK)模型或稱為單松弛(SRT)模型[19—23]、熵模型[24]、雙松弛(TRT)模型[25,26]、多松弛(MRT)模型[27—34]等.其中LBGK模型因其形式的簡(jiǎn)潔性在各種復(fù)雜的流體傳輸問題研究中受到了廣泛的應(yīng)用,其中最具有代表性的是Qian等[20]提出的DdQq模型.但由于LBGK模型以單松弛時(shí)間近似為基礎(chǔ),其在穩(wěn)定性和精確性方面存在不足.幾乎同時(shí),法國學(xué)者d'Humières[27]提出了廣義的格子Boltzmann模型即多松弛格子Boltzmann模型(MRT).MRT模型和LBGK模型的不同之處主要體現(xiàn)在它們的碰撞項(xiàng),MRT模型具有數(shù)量最多的自由度,即可調(diào)松弛因子的數(shù)量最多,而不同的松弛因子可以最大限度地優(yōu)化模型性質(zhì),如穩(wěn)定性和精確性[30,35],因此MRT模型受到越來越多的重視.常見的MRT模型有二維的D2Q9 MRT模型[27],三維的D3Q15 MRT模型[28]、D3Q19 MRT模型[28],至今為止三維空間中速度方向最少的MRT模型是由法國學(xué)者d'Humières等[29]提出的D3Q13 MRT模型.這些模型都可以通過Chapman-Enskog(C-E)展開恢復(fù)到Navier-Stokes方程 (N-S方程).速度集合多的模型具有更好的穩(wěn)定性和各向同性[36,37],但同時(shí)在計(jì)算時(shí)間和存儲(chǔ)空間的消耗上會(huì)增加.為此有學(xué)者提出了擁有更少速度方向的MRT模型,如D2Q8 MRT模型[32]、D3Q14 MRT模型和D3Q18 MRT模型[33].這幾種MRT模型的構(gòu)造原理是基于Guo等[22]提出的不可壓DdQq LBGK模型的宏觀量計(jì)算與分布函數(shù) f0無關(guān),且在DdQq MRT模型[28]恢復(fù)宏觀方程的過程中沒有用到能量的平方項(xiàng).因此將DdQq MRT模型的離散速度集合舍棄0速度方向、模型的矩中能量平方項(xiàng)丟掉,并且將它的變換矩陣去掉第一列和能量平方項(xiàng)所對(duì)應(yīng)的那一行再正交化得到新的DdQ(q-1) MRT模型的變換矩陣.最后將構(gòu)造的變換矩陣乘以原LBGK模型中去掉0方向的平衡態(tài)分布函數(shù)就得到DdQ (q-1) MRT模型的矩平衡態(tài).

本文在D3Q13 MRT模型[29]的基礎(chǔ)上提出了求解不可壓縮N-S方程的三維12速M(fèi)RT模型(iD3Q12 MRT),和已有的D2Q8,D3Q14,D3Q18等MRT模型的構(gòu)造方法相比,最大的難點(diǎn)在于D3Q13 MRT模型沒有對(duì)應(yīng)的LBGK模型,因此就不能通過變換矩陣乘以平衡態(tài)分布函數(shù)的方法得到矩平衡態(tài).但可仿照D2Q8,D3Q14,D3Q18 MRT模型的構(gòu)造方法構(gòu)造離散速度集合和變換矩陣,舍棄了原13速M(fèi)RT模型的矩中能量e 的那一項(xiàng),并且在13速M(fèi)RT模型沒有對(duì)應(yīng)的LBGK模型的情況下運(yùn)用反演法,即用13速M(fèi)RT模型變換矩陣的逆乘它的矩平衡態(tài)得到形式上的“平衡態(tài)分布函數(shù)”,再通過一系列構(gòu)造變成含有12個(gè)元素的“平衡態(tài)分布函數(shù)”,最后將構(gòu)造出的12速M(fèi)RT模型的變換矩陣乘以這個(gè)“平衡態(tài)分布函數(shù)”就得到了矩平衡態(tài),使得新的iD3Q12 MRT模型可以在低馬赫數(shù)(Ma)條件下通過C-E展開恢復(fù)到不可壓N-S方程.需要注意的是,本文在模擬部分與iD3Q12 MRT模型作對(duì)比的是D3Q13 MRT模型的不可壓版本,這里稱之為He-Luo D3Q13 MRT模型.這是因?yàn)樗木仄胶鈶B(tài)是將D3Q13 MRT模型的矩平衡態(tài)表達(dá)式中與速度相乘的密度ρ變成 ρ0,這與He-Luo LBGK模型[21]的平衡態(tài)分布函數(shù)的構(gòu)造方法相同.這種變化忽略了平衡態(tài)分布函數(shù)中Ma 三次方的同階或高階無窮小量,從而減少了LB模型的可壓縮效應(yīng).因此在模擬部分,選取可壓縮誤差更小的He-Luo D3Q13 MRT模型與iD3Q12 MRT模型作對(duì)比.

總之,基于D3Q13 MRT模型,運(yùn)用反演法提出了一個(gè)可以在低Ma 假設(shè)下恢復(fù)到不可壓N-S方程的iD3Q12 MRT模型,它可能是目前三維MRT模型中離散速度方向個(gè)數(shù)最少的一個(gè)模型,因此iD3Q12 MRT模型在計(jì)算量和存儲(chǔ)量的需求上更小.通過一系列數(shù)值模擬,我們將iD3Q12與He-Luo D3Q13 MRT模型作對(duì)比,驗(yàn)證了我們提出的iD3Q12 MRT模型的有效性,并考察了該模型在精確性和穩(wěn)定性方面與He-Luo D3Q13 MRT模型的差異.

2 三維13速M(fèi)RT模型

D3Q13 MRT模型是d'Humières等[29]提出的.在現(xiàn)有的MRT模型中,它滿足伽利略不變性和各向同性,并且是能夠通過C-E展開恢復(fù)到Navier-Stokes方程的具有最少的離散速度的一個(gè)MRT模型.在假設(shè)空間步長(zhǎng) δx =1 和時(shí)間步長(zhǎng)δt =1,即粒子速度 c =δx /δt =1 的情形下,D3Q13 MRT模型選取的離散速度如下：

模型的演化方程為

fi(x,t)是沿速度 ci移動(dòng)的粒子的分布函數(shù),是平衡態(tài)分布函數(shù),Λij是 13× 13 階碰撞矩陣Λ的元素.由速度方向 ci生成了兩兩正交的向量 ek,k∈ 0, 1,···, 12 ,再由這13個(gè)正交向量ek定義D3Q13 MRT模型的變換矩陣 T,

作用于分布函數(shù)f 得到13個(gè)矩mk=ek·f =

選取的矩平衡態(tài) m(eq)為

通過C-E展開,并在低Ma 假設(shè)下,該模型能夠恢復(fù)到可壓縮的Navier-Stokes方程：

其中

分別表示剪切黏度和體黏度,Λν和是與剪切黏度ν相關(guān)的碰撞因子.

He-Luo D3Q13 MRT模型的矩平衡態(tài)選取如下：

通過C-E展開可以將He-Luo D3Q13 MRT模型恢復(fù)到如下形式的Navier-Stokes方程：

其中

在低Ma 假設(shè)和 T?L/cs(T和L分別表示特征時(shí)間和特征長(zhǎng)度)的條件下,He-Luo D3Q13 MRT模型可以進(jìn)一步恢復(fù)到不可壓N-S方程：

在He-Luo D3Q13 MRT模型中,宏觀量的計(jì)算格式如下：

在后面的模擬中,假設(shè) ρ0=1.

3 反演法：iD3Q12 MRT模型矩平衡態(tài)的構(gòu)造

將D3Q13 MRT模型變換矩陣的逆乘它的矩平衡態(tài)得到形式上的“平衡態(tài)分布函數(shù)”,如下所示：

由于He-Luo LBGK模型[21]的平衡態(tài)分布函數(shù)的構(gòu)造是基于Qian等[20]的LBGK模型,并將其平衡態(tài)分布函數(shù)中速度u 前的密度ρ替換為 ρ0.因此,將(12)式中速度u 前的ρ替換成 ρ0,得到He-Luo LBGK模型式的“平衡態(tài)分布函數(shù)”,如下所示：

由于不可壓LBGK模型[22]的平衡態(tài)分布函數(shù)可以看作是在He-Luo LBGK模型的基礎(chǔ)上,非0方向的平衡態(tài)分布函數(shù)中的ρ換成 p /,與速度相乘的 ρ0令成1,0方向的“平衡態(tài)分布函數(shù)”定義為ρ0減去其他方向的平衡態(tài)分布函數(shù)之和.按照這種方法,變化(13)式,則得到不可壓LBGK模型式的“平衡態(tài)分布函數(shù)”,如下所示：

去掉(14)式中的 f0,得到下面1—12方向上的“平衡態(tài)分布函數(shù)”

4 三維12速M(fèi)RT模型

根據(jù)d'Humières等[29]提出的最原始的D3Q13 MRT模型,我們構(gòu)造出了新的三維12速多松弛格子Boltzmann模型(iD3Q12 MRT模型).首先將13速模型的0速度方向舍棄得到了12速模型的離散速度方向的集合：

不失一般性,這里假設(shè)粒子速度 c =1.iD3Q12 MRT模型的演化方程為

fi(x,t)是沿速度 ci移動(dòng)的粒子的分布函數(shù),是平衡態(tài)分布函數(shù),Λij是 12×12 階碰撞矩陣Λ的元素.同時(shí),演化方程也可以寫作向量形式：

其中|f (x,t)〉=(f1(x,t),f2(x,t),···,f12(x,t))′是列向量,符號(hào)′代表轉(zhuǎn)置算子.構(gòu)造一個(gè)MRT模型,變換矩陣和矩平衡態(tài)的構(gòu)造至關(guān)重要.按照構(gòu)造D2Q8[32],D3Q14和D3Q18 MRT模型[33]的變換矩陣的方法構(gòu)造了iD3Q12 MRT模型的變換矩陣,在構(gòu)造的過程中,舍棄了13速M(fèi)RT模型的矩中能量那一項(xiàng),得到的iD3Q12 MRT模型的變換矩陣如下：

通過變換矩陣T將分布函數(shù)f 變換為矩陣m

Λc是與守恒矩相關(guān)的松弛因子,由于守恒矩對(duì)應(yīng)的平衡態(tài)是它本身,因此這些松弛因子的取值不影響粒子的碰撞過程,可以將它們?nèi)≈禐?.Λν,和Λt是與非守恒矩相匹配的松弛因子,為了保持模型的穩(wěn)定性,取值一般在(0,2).我們選取的矩平衡態(tài)如下：

p 是壓力,|u |2是速度 u =(ux,uy,uz) 模的平方.演化方程在矩空間的形式如下：

其中

若令 τ=1/Λν,則

在低Ma 條件下,方程(23)可以寫成如下形式：

在iD3Q12 MRT模型中,宏觀量的計(jì)算如下：

5 數(shù)值模擬

為了驗(yàn)證提出的iD3Q12 MRT模型的有效性,我們模擬了三維泊肅葉流、脈動(dòng)流與頂蓋驅(qū)動(dòng)的方腔流.值得注意的是,在模型的建立和推導(dǎo)過程中,假設(shè) c =δx /δt =1.但在實(shí)際的計(jì)算中,c 有時(shí)不取1.這種情況下,我們將宏觀速度u 、壓力p 進(jìn)行單位化 u′=u /c ,p′=p /c2,再將單位化后的量代入我們提出的模型進(jìn)行計(jì)算,計(jì)算結(jié)束后再將所得到的 u′,p′分別乘以c ,c2,就得到我們所要求的宏觀速度和壓力.同時(shí),剪切黏性系數(shù)按照確定.此外,碰撞矩陣中的松弛因子取定如下：Λc=1.0 ,Λν=1/τ,Λt=1.8.我們也給出了He-Luo D3Q13 MRT模型計(jì)算得到的有關(guān)結(jié)果.對(duì)于邊界條件,使用的是非平衡態(tài)外推法[38].程序的計(jì)算流程如下：

1) 初始化 輸入計(jì)算參數(shù),宏觀速度u 和壓力p 的初始值都設(shè)置為0,并初始化分布函數(shù)f ；

2)矩空間碰撞

3) 速度空間遷移

4) 邊界處理采用非平衡態(tài)外推法.

5.1 三維穩(wěn)態(tài)的泊肅葉流

如圖1所示,三維泊肅葉流的物理空間限制在一個(gè)長(zhǎng)方體通道中,其長(zhǎng)、寬、高分別為 0≤x ≤l ,-a ≤y ≤a ,-b ≤z ≤b ,l =2 ,a =b =0.5 ,原點(diǎn)O代表流體入口平面的中心.

邊界條件設(shè)置為

圖1 三維泊肅葉流示意圖Fig.1.The schematic of three-dimensional Poiseuille flow.

pin和 pout分別表示進(jìn)出口壓力,三維泊肅葉流具有穩(wěn)態(tài)的解析解[39]是通道內(nèi)的壓力梯度,ν是流體的剪切黏度.在模擬中,初始狀態(tài)是：

dp /dx =(pout-pin)/l

并且設(shè) pin=1.1 和 pout=1.0 ,模擬達(dá)到穩(wěn)定狀態(tài)的判定準(zhǔn)則是

我們也計(jì)算了iD3Q12 MRT模型在模擬泊肅葉流時(shí)的空間精度的階.為此,計(jì)算了不同空間步長(zhǎng)下的泊肅葉流的速度場(chǎng)的全局相對(duì)誤差 GREu,GREu的表達(dá)式為

公式中的n 和a 分別表示數(shù)值解和解析解,ux,uy和 uz是流體速度u 的三個(gè)分量,同樣的關(guān)于i 的求和遍及所有網(wǎng)格點(diǎn).在不同松弛參數(shù)和不同空間步長(zhǎng)下由iD3Q12 MRT和D3Q13 MRT模型計(jì)算的泊肅葉流的速度場(chǎng)的全局相對(duì)誤差 GREu列在表1中.值得注意的是,在同樣的參數(shù)下,這兩個(gè)模型的全局相對(duì)誤差是完全一樣的.

圖2 泊肅葉流數(shù)值解與解析解的對(duì)比 (a) 泊肅葉流在 x =1 截面處z 取不同的值時(shí)水平速度 ux隨y 變化的函數(shù)圖像；(b) 在截面 z =0 處y 取不同的值時(shí)壓力p 隨x 變化的函數(shù)圖像；直線：解析解；符號(hào)：數(shù)值解；松弛因子Λν=1.3Fig.2.Comparison between numerical and analytical solutions of Poiseuille flow：(a) The variation of ux with y for different locations of z at section x =1 for Poiseuille flow；(b) the variation of pressure with x for different locations of y at section z =0 for Poiseuille flow.Lines,analytical solutions；symbols,numerical results；the relaxation parameter Λν=1.3.

表1 iD3Q12 MRT和D3Q13 MRT模型在不同松弛因子 Λν和不同空間步長(zhǎng)下計(jì)算得到的泊肅葉流的速度場(chǎng)的全局相對(duì)誤差GREuTable 1.The GREu of velocity field for Poiseuille flow computed by iD3Q12 MRT and D3Q13 MRT models under different relaxation parameters and different lattice spacings.

圖3 不同的 Λν下,模擬泊肅葉流得到的速度場(chǎng)的全局相對(duì)誤差 GREu隨空間步長(zhǎng) δx 的變化,符號(hào)代表數(shù)值解,連線表示擬合直線Fig.3.The variation of GREu of velocity field with the lattice spacing δx at different Λν for Poiseuille flow.Symbols represent numerical solutions,lines represent fitting line.

假設(shè) GREu=a (δx)b(a＞ 0) ,然后我們得到ln(GREu)=b ln(δx)+ln(a).如果 b =2 ,就稱數(shù)值解具有二階精度.用表1中的數(shù)據(jù)進(jìn)行最小二乘法的線性擬合,線性擬合的圖像在圖3中顯示.圖中三條直線分別對(duì)應(yīng)著iD3Q12 MRT模型在Λν=0.8, 1.0和1.3下的擬合直線,它們對(duì)應(yīng)的斜率分別為1.94,1.91和1.88,都接近2.這說明我們提出的iD3Q12 MRT模型在模擬穩(wěn)態(tài)的泊肅葉流時(shí)能夠達(dá)到二階的空間精度.

5.2 三維非穩(wěn)態(tài)的脈動(dòng)流

模擬三維脈動(dòng)流是為了驗(yàn)證所提出的iD3Q12 MRT模型在模擬非穩(wěn)態(tài)流時(shí)的精度.脈動(dòng)流的流域與泊肅葉流的流域相同,如圖1所示.但脈動(dòng)流在管道進(jìn)出口兩端的壓力梯度是周期性變化的,壓力梯度為

其中G是振幅,ω是頻率.脈動(dòng)流的解析解是[40]

其中

關(guān)于i 的求和遍及所有網(wǎng)格點(diǎn).圖4顯示了在四個(gè)不同時(shí)刻 t =T/4,T/2, 3T/4 和T下,在直線 x =1 ,z =0上水平速度 ux隨y 變化的函數(shù)圖像,此時(shí)η=2.8285,λν=1.522.水平速度 ux由 Umax無量綱化,Umax=1.876× 10-2是泊肅葉流在進(jìn)出口壓差為 pout-pin=-G*l 時(shí)的最大水平速度.從圖4可以得出iD3Q12 MRT模型在模擬非穩(wěn)態(tài)流時(shí)得到的數(shù)值解與解析解符合較好.

表2給出了在 η=2.8285 且 τ=1/Λν、ν一定的條件下,用iD3Q12 MRT和D3Q13 MRT模型模擬脈動(dòng)流時(shí)得到的速度場(chǎng)的全局相對(duì)誤差GREu.表3顯示了由表2中的數(shù)據(jù)計(jì)算得到的相鄰的兩個(gè)空間步長(zhǎng)下的iD3Q12 MRT和D3Q13 MRT模型的空間精度的階這里f 表示小的空間步長(zhǎng),c 表示大的空間步長(zhǎng).從表2中可以看出,在四個(gè)不同時(shí)刻下由iD3Q12 MRT模型計(jì)算的速度場(chǎng)的全局相對(duì)誤差都比D3Q13 MRT模型計(jì)算的速度場(chǎng)的全局相對(duì)誤差要略小一些.從表3中可以發(fā)現(xiàn)D3Q13 MRT模型比iD3Q12 MRT模型的階稍微更接近于2.根據(jù)表2中的數(shù)據(jù),我們繪制了用iD3Q12 MRT模型計(jì)算得到的速度場(chǎng)的全局相對(duì)誤差 GREu隨空間步長(zhǎng) δx 變化的函數(shù)圖像,如圖5所示.在四個(gè)不同時(shí)刻t = T/4,T/2,3T/4和T下擬合直線的斜率分別是1.98,1.88,1.97和1.87,這說明iD3Q12 MRT模型在模擬非穩(wěn)態(tài)脈動(dòng)流時(shí)具有二階的空間精度.

圖4 在 η=2.8285 時(shí)脈動(dòng)流在 x =1 ,z =0 處水平速度ux隨y 變化的函數(shù).直線：解析解；符號(hào)：數(shù)值解Fig.4.The variation of horizontal velocity ux with y for pulsatile flow at the location x =1 ,z =0 ,η=2.8285.Line,analytical solutions；symbols,numerical solutions.

表2 在 η=2.8285 時(shí),不同空間步長(zhǎng)下用iD3Q12 MRT模型和D3Q13 MRT模型模擬脈動(dòng)流所得的不同時(shí)刻下的速度場(chǎng)的全局相對(duì)誤差GREuTable 2.The global relative errors of the velocity field at different times for pulsatile flow simulated by iD3Q12 MRT and D3Q13 MRT models at different lattice spacings,η=2.8285.

圖5 同一周期四個(gè)不同時(shí)刻下變量 GREu隨空間步長(zhǎng)δx 的變化Fig.5.The variation of GREu with the lattice spacing at four different times in a period for pulsatile flow.

當(dāng)最大的壓差 Δp 增大時(shí),LB模型計(jì)算的速度場(chǎng)的可壓縮性誤差會(huì)增大,在這種情況下比較iD3Q12 MRT和D3Q13 MRT模型計(jì)算的速度場(chǎng)的全局相對(duì)誤差 GREu能凸顯兩個(gè)模型在精度上的差異.表4和表5分別顯示了最大的壓差 Δp 增大時(shí)iD3Q12 MRT和D3Q13 MRT模型在不同的空間步長(zhǎng)和松弛時(shí)間τ下所計(jì)算得到的T時(shí)刻的速度場(chǎng)的全局相對(duì)誤差.從表4可以看出,由iD3Q12 MRT模型計(jì)算的脈動(dòng)流在T時(shí)刻的全局相對(duì)誤差總是比D3Q13 MRT模型要略小一些.這說明iD3Q12 MRT模型的計(jì)算精度比D3Q13 MRT模型的計(jì)算精度要高.隨著最大壓差 Δp 的增加iD3Q12 MRT模型的計(jì)算開始發(fā)散但D3Q13 MRT模型卻沒有發(fā)散,如 δx =1/60 且 Δp =0.08 時(shí).當(dāng)空間步長(zhǎng)變?yōu)?δx =1/80 時(shí)再繼續(xù)增大最大壓差Δp 同樣會(huì)出現(xiàn)類似的情形,如 Δp =0.10 和Δp =0.12,這說明D3Q13 MRT模型在穩(wěn)定性方面較iD3Q12 MRT模型更好.從表5的數(shù)據(jù)可以看出不同的τ下,iD3Q12 MRT模型計(jì)算得到的速度場(chǎng)的全局相對(duì)誤差在絕大部分情況下都比D3Q13 MRT模型的要略小一些,但也有個(gè)別反常情形,如 τ=0.70 且 Δp =0.01 或 Δp =0.04 時(shí).在表5中還發(fā)現(xiàn)和表4相似的情形,即在 τ=0.55 時(shí)增大最大壓差 Δp 出現(xiàn)了iD3Q12 MRT模型先發(fā)散的情形,如 Δp =0.03 時(shí).從表5中還可以發(fā)現(xiàn),當(dāng)τ更小時(shí),增加 Δp ,iD3Q12 MRT和D3Q13 MRT模型容易發(fā)散.因此,在 Δp 較大時(shí),適當(dāng)?shù)卦龃螃涌梢蕴嵘齼蓚€(gè)模型的穩(wěn)定性.

表3 相鄰空間步長(zhǎng)下的iD3Q12 MRT和D3Q13 MRT模型的空間精度的階Table 3.The orders of the spatial accuracy of iD3Q12 MRT and D3Q13 MRT models under adjacent spacings.

表4 在 τ=0.5667 ,η=4.3416 ,最大壓差 Δp 增大時(shí)不同的空間步長(zhǎng)下由iD3Q12 MRT和D3Q13 MRT模型模擬的脈動(dòng)流在時(shí)刻T下的速度場(chǎng)所計(jì)算的全局相對(duì)誤差 GREu,空白處表示計(jì)算發(fā)散Table 4.The global relative error calculated by the velocity field at time T of pulsatile flow simulated by the iD3Q12 MRT and D3Q13 MRT models under different lattice spacings.The maximal pressure drop Δp of the channel increases,τ=0.5567,η=4.3416 are fixed.The blank indicates that the computation is divergent.

表5 在 δx =1/20 時(shí),最大壓差 Δp 增大時(shí)不同的松弛時(shí)間τ下由iD3Q12 MRT和D3Q13 MRT模型模擬的脈動(dòng)流由T時(shí)刻的速度場(chǎng)計(jì)算得出的全局相對(duì)誤差 GREu,空白處表示計(jì)算發(fā)散Table 5.The global relative error of the velocity field at time T of the pulsatile flow simulated by the iD3Q12 MRT and D3Q13 MRT models under different relaxation time τ.The maximal pressure drop of the channel is increased and δx =1/20is fixed.The blank indicates that the computation is divergent.

5.3 三維頂蓋驅(qū)動(dòng)的方腔流

三維方腔流包含旋渦運(yùn)動(dòng),Ku等[41]采用偽普方法計(jì)算的結(jié)果常被用作新數(shù)值方法模擬方腔流精度的評(píng)判標(biāo)準(zhǔn),因此我們也采用Ku的計(jì)算結(jié)果來檢驗(yàn)iD3Q12 MRT模型的精度.

流動(dòng)是在一個(gè)立方體盒子中進(jìn)行的,如圖6所示.流體是由最頂端的蓋子以常速度 U0=1.0 移動(dòng)而被驅(qū)動(dòng).流動(dòng)滿足三維不可壓N-S方程.雷諾數(shù)可由 Re =U0L/ν計(jì)算,這里 L=1.0 是立方體盒子的長(zhǎng)度,ν是剪切黏度.在模擬中,初始狀態(tài)的速度和壓力都設(shè)置為0.在 z =0.5 的截面上,對(duì)稱的邊界條件設(shè)置為

在 y =1 的截面上 x =0 和 x =1 處 u =0.在y =1且 x =0 和 x =1 旁邊的點(diǎn)處分別設(shè)置ux=0.3和 ux=1.0.初始條件和邊界條件的設(shè)置與Ku等[41]的相同.此外,為了滿足低Ma 假設(shè),固定 c =10 ,c／=1 時(shí)MRT模型的計(jì)算見參考文獻(xiàn)[31].速度域的收斂準(zhǔn)則是：

圖6 三維頂蓋驅(qū)動(dòng)的方腔流示意圖Fig.6.The schematic of three-dimensional lid-driven cavity flow.

圖7 不同的雷諾數(shù)下模擬方腔流,在截面 z =0.5 處豎直和水平中心線的速度分布 (a) Re = 100；(b) Re =400；(c)Re =1000Fig.7.The velocity distribution in the vertical and horizontal center lines at section z =0.5 for cavity flows at different Re ：(a) Re =100；(b) Re =400；(c) Re =1000.

表6 不斷增大雷諾數(shù)比較iD3Q12 MRT和He-Luo D3Q13 MRT模型在模擬方腔流時(shí)的穩(wěn)定性.?代表收斂,收斂準(zhǔn)則是(39)式Table 6.Comparing the stability of iD3Q12 MRT and He-Luo D3Q13 MRT models for three-dimensional cavity flows when the Reynolds number is continuously increased.The tick represents convergence,the convergence criterion is formula (39).

對(duì)i 和k 的求和遍及所有的網(wǎng)格點(diǎn)和所有的速度方向.首先采用iD3Q12 MRT模型模擬了Re =100,400和1000時(shí)的流動(dòng),并將模擬結(jié)果與已有的Ku等[41]的結(jié)果作對(duì)比,如圖7所示.可以看出由iD3Q12 MRT模型計(jì)算的數(shù)值結(jié)果與Ku等[41]的結(jié)果符合得很好,因此我們提出的iD3Q12 MRT模型在模擬三維頂蓋驅(qū)動(dòng)的方腔流時(shí)是準(zhǔn)確的.

不斷地增大雷諾數(shù)時(shí)我們觀察并記錄了iD3Q12 MRT和He-Luo D3Q13 MRT模型在模擬方腔流時(shí)的斂散性,結(jié)果顯示在表6中.從中可以看到在增大雷諾數(shù)時(shí)iD3Q12 MRT模型的穩(wěn)定性比He-Luo D3Q13 MRT模型的穩(wěn)定性稍弱一些.

6 結(jié) 論

提出了求解三維不可壓縮流的12速多松弛格子Boltzmann模型,即iD3Q12 MRT模型,它可能是現(xiàn)有的能推導(dǎo)出不可壓N-S方程的三維MRT模型中離散速度方向最少的一個(gè)模型,因此原則上具有更高的計(jì)算效率.構(gòu)建iD3Q12 MRT模型的基本思想是構(gòu)造一個(gè)12×12階的正交變換矩陣,并選取合適的矩平衡態(tài)使新模型能通過C-E展開恢復(fù)到不可壓的N-S方程.用iD3Q12 MRT模型分別模擬了穩(wěn)態(tài)的泊肅葉流、非穩(wěn)態(tài)的脈動(dòng)流、頂蓋驅(qū)動(dòng)的方腔流,驗(yàn)證了模型的準(zhǔn)確性和穩(wěn)定性,并將其準(zhǔn)確性和穩(wěn)定性與He-Luo D3Q13 MRT模型作了對(duì)比.

對(duì)于泊肅葉流和脈動(dòng)流,iD3Q12 MRT模型的數(shù)值解和解析解符合得很好,這說明iD3Q12 MRT模型在模擬穩(wěn)態(tài)流和非穩(wěn)態(tài)流時(shí)都是準(zhǔn)確的.在不同的空間步長(zhǎng) Δx 和松弛因子 Λν下用iD3Q12 MRT模型和He-Luo D3Q13 MRT模型模擬泊肅葉流,發(fā)現(xiàn)兩個(gè)模型的全局相對(duì)誤差 GREu完全相同.用兩個(gè)模型模擬了脈動(dòng)流并計(jì)算了不同空間步長(zhǎng)Δx 和不同時(shí)刻下的全局相對(duì)誤差,發(fā)現(xiàn)iD3Q12 MRT模型的全局相對(duì)誤差 GREu比He-Luo D3Q13 MRT模型稍小一些.全局相對(duì)誤差隨空間步長(zhǎng)的變化表明兩個(gè)模型在模擬穩(wěn)態(tài)的泊肅葉流和非穩(wěn)態(tài)的脈動(dòng)流時(shí)都具有二階精度.還通過增大管道最大壓降Δp 并調(diào)節(jié)空間步長(zhǎng) Δx 或松弛時(shí)間τ的方式模擬了脈動(dòng)流,發(fā)現(xiàn)在大部分的參數(shù)下iD3Q12 MRT模型計(jì)算的流場(chǎng)的全局相對(duì)誤差 GREu比He-Luo D3Q13 MRT模型計(jì)算的流場(chǎng)的全局相對(duì)誤差要小一些,但隨著最大壓降的增大iD3Q12 MRT模型比D3Q13 MRT模型先發(fā)散.這說明在模擬非穩(wěn)態(tài)流時(shí)iD3Q12 MRT模型比D3Q13 MRT模型的準(zhǔn)確性稍高但在穩(wěn)定性方面弱一些.此外,用iD3Q12 MRT模型模擬了包含旋渦運(yùn)動(dòng)的三維頂蓋驅(qū)動(dòng)方腔流,發(fā)現(xiàn)iD3Q12 MRT模型的模擬結(jié)果與已有的Ku等的結(jié)果符合得很好,并且iD3Q12 MRT模型能模擬的最大雷諾數(shù)比He-Luo D3Q13 MRT模型的要稍小.

附錄 iD3Q12 MRT模型的Chapman-Enskog展開

通過C-E展開將iD3Q12 MRT模型恢復(fù)到三維不可壓的N-S方程.對(duì)于不可壓縮流有

這里Ma 代表馬赫數(shù),δp 和 δρ分別代表壓力和密度的脈動(dòng)量.

首先引入展開

這里 ε=δt ,Di≡?t+ci·?=?t+ciα·?α=?t+cix?x+ciy?y+ciz?z,運(yùn)用上述(A2)式展開,并將演化方程

做Taylor展開,按照ε的不同階可得：

這里Di0≡?t0+ci·?=?t+ciα·?α=?t0+cix?x+ciy?y+ciz?z,并且將方程(A4b)代入到方程(A4c)中,方程(A4)可以轉(zhuǎn)化到矩空間中：

值得注意的是,P1和 ux,uy,uz在低Ma 下是守恒量,這樣在表達(dá)式(A6)中都為0.展開方程(A5b)有

(A7)式的前4個(gè)方程是

從(A1a)與(A1b) 可以得到 O(δp)=O(M2) 和O(u)=O(M),有 ?t0(3p /2+u2/2)=O(M2).忽略 O(M2) 項(xiàng),方程(A8)變成

它是不可壓N-S方程的連續(xù)方程.從方程(A5c)中可得

由(A9)式 + ε×(A11)式可得

根據(jù)方程(A7),算出

使用方程(A9),得到

將方程(A15)代入(A12)式,可得

其中,

由(A10)式,省略方程(A16)中 ξ(?(?·u)) 和u?·u 兩項(xiàng),得

這就是不可壓N-S方程動(dòng)量方程的分量形式.至此,在低馬赫數(shù)假設(shè)下,已經(jīng)通過C-E展開將iD3Q12 MRT模型恢復(fù)到不可壓的N-S方程,它可以寫作向量的形式,見方程(25).