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

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

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

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

1 引 言

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

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

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

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

2 三維13速MRT模型

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

模型的演化方程為

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

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

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

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

其中

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

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

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

其中

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

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

在后面的模擬中,假設 ρ0=1.

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

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

由于He-Luo LBGK模型[21]的平衡態(tài)分布函數(shù)的構造是基于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模型的基礎上,非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速MRT模型

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

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

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

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

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

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

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

其中

若令 τ=1/Λν,則

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

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

5 數(shù)值模擬

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

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

2)矩空間碰撞

3) 速度空間遷移

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

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

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

邊界條件設置為

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

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

dp /dx =(pout-pin)/l

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

我們也計算了iD3Q12 MRT模型在模擬泊肅葉流時的空間精度的階.為此,計算了不同空間步長下的泊肅葉流的速度場的全局相對誤差 GREu,GREu的表達式為

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

圖2 泊肅葉流數(shù)值解與解析解的對比 (a) 泊肅葉流在 x =1 截面處z 取不同的值時水平速度 ux隨y 變化的函數(shù)圖像；(b) 在截面 z =0 處y 取不同的值時壓力p 隨x 變化的函數(shù)圖像；直線：解析解；符號：數(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模型在不同松弛因子 Λν和不同空間步長下計算得到的泊肅葉流的速度場的全局相對誤差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 不同的 Λν下,模擬泊肅葉流得到的速度場的全局相對誤差 GREu隨空間步長 δx 的變化,符號代表數(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.

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

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

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

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

其中

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

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

圖4 在 η=2.8285 時脈動流在 x =1 ,z =0 處水平速度ux隨y 變化的函數(shù).直線：解析解；符號：數(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 時,不同空間步長下用iD3Q12 MRT模型和D3Q13 MRT模型模擬脈動流所得的不同時刻下的速度場的全局相對誤差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 同一周期四個不同時刻下變量 GREu隨空間步長δx 的變化Fig.5.The variation of GREu with the lattice spacing at four different times in a period for pulsatile flow.

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

表3 相鄰空間步長下的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 增大時不同的空間步長下由iD3Q12 MRT和D3Q13 MRT模型模擬的脈動流在時刻T下的速度場所計算的全局相對誤差 GREu,空白處表示計算發(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 時,最大壓差 Δp 增大時不同的松弛時間τ下由iD3Q12 MRT和D3Q13 MRT模型模擬的脈動流由T時刻的速度場計算得出的全局相對誤差 GREu,空白處表示計算發(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 三維頂蓋驅動的方腔流

三維方腔流包含旋渦運動,Ku等[41]采用偽普方法計算的結果常被用作新數(shù)值方法模擬方腔流精度的評判標準,因此我們也采用Ku的計算結果來檢驗iD3Q12 MRT模型的精度.

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

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

圖6 三維頂蓋驅動的方腔流示意圖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模型在模擬方腔流時的穩(wě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).

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

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

6 結 論

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

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

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

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

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

首先引入展開

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

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

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

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

(A7)式的前4個方程是

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

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

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

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

使用方程(A9),得到

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

其中,

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

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