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

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

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

** (廣東省湍流基礎研究與應用重點實驗室,廣東深圳 518055)

引言

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

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

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

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

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

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

1 理論基礎

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

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

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

不同離散速度點對應的權重為

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

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

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

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

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

1.1 BGK 碰撞模型

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

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

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

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

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

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

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

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

而通過Hermite 展開得到的力項表達式 (18),也與Guo 等[32]通過反設計得到的結果完全一致.

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

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

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

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

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

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

1.2 MRT 碰撞模型

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

MRT 模型對應的速度離散后的碰撞項可以表示如下[36]

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

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

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

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

若記

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

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

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

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

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

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

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

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

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

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

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

2 局部網格加密中的數(shù)據(jù)轉換

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

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

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

下文中,與粗網格和細網格區(qū)域相關的任何量分別用上標“c”和“f”表示.設局部加密時粗細網格的空間步長比例為

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

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

2.1 BGK 模型下的局部網格加密

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

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

考慮CE 展開如下

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

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

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

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

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

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

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

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

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

2.2 MRT 模型下的局部網格加密

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

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

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

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

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

3 局部網格加密算法的數(shù)值驗證

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

圖2 粗細網格排布示意圖.灰色區(qū)域為粗細網格的重疊區(qū)域Fig.2 Schematic diagram of the coarse and fine grid arrangement.The grey area is the overlap between the coarse and fine grids

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

3.1 強迫泰勒-格林渦流動

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

該流動問題對應的控制方程為

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

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

(1)Q=1 時,代表沒有外力,各項宏觀量按標準的泰勒-格林渦流動衰減;

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

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

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

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

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

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

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

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

圖3 強迫泰勒-格林渦各宏觀量的分布云圖(Q=0.5, tU0/L=1,黑色矩形框內的區(qū)域為加密區(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)

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

圖4 強迫泰勒-格林渦中在 (x/L,y/L)=(0.4,0.25) 點處各宏觀量隨時間的變化.(a)～(c)分別為BGK 模型和MRT 模型下的相對壓力P-P0、水平方向速度 ux 和黏性應力 σxx 與理論結果的對比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

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

圖5 強迫泰勒-格林渦中各宏觀量的剖線圖(tU0/L=1, y/L=0.25,兩條黑色虛線間的區(qū)域為加密區(qū)域).(a)～(c)分別為BGK 模型和MRT 模型下相對壓力 P-P0、水平方向速度 ux 和黏性應力 σxx 與理論結果的對比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ù)值結果也印證了第 2 節(jié)中粗細網格間數(shù)據(jù)轉換關系的正確性.接下來將以平板泊肅葉流中對流-擴散問題為驗證算例,進一步測試BGK 模型下粗細網格間數(shù)據(jù)轉換關系對復雜源項的適應性.

3.2 平板泊肅葉流中對流-擴散問題

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

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

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

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

這里 ?? 可以直接使用如下公式計算

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

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

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

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

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

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

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

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

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

圖7 平板泊肅葉流中對流-擴散問題(黑色虛線為加密界面).(a)～(d)分別為BGK 模型下水平方向速度 ux、標量 ?、通量 Jx 和 Jy 與理論結果的對比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 頂蓋驅動方腔流動

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

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

圖8 頂蓋驅動方腔流4 個角點示意圖,其中黑色粗實線代表壁面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) 代表右上角點對應的坐標,而 2,3,6 方向均使用靜止邊界的反彈格式.左上角點也是同樣處理.

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

圖9 Re=1000 時頂蓋驅動方腔流的壓力分布云圖對比,(a)～(d)分別為UCG,UFG,EQ 和UCG-L 對應的壓力分布云圖 ((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 中繪制了使用不同分辨率時沿垂直方向穿過幾何中心的水平速度ux隨y/L的變化關系以及沿水平方向穿過幾何中心的垂直速度uy隨x/L的變化關系.而且從圖10(b) 以及圖10(d) 中可以很清楚地看到即使是遠離加密區(qū)域UCG-L 的數(shù)值結果仍要好于與其所需計算量相同的均勻網格下的結果,而且在使用更少計算資源的情況下還取得了與UFG 非常相近的數(shù)值結果,這就表明在局部流域執(zhí)行網格加密可以很好地減小局部數(shù)值誤差對全流域的影響.

圖10 Re=1000 時頂蓋驅動方腔流中無量綱速度的剖線圖,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 個迭代步的 CPU 實際消耗時間.如果假設在 UCG 上運行到穩(wěn)態(tài)所需的 CPU 時間是t,若考慮對流時間尺度以及二維情況下則UFG 下運行到相同狀態(tài)下,理論上所需的 CPU 時間是 8t,而上述加密的UCG-L只需要 3.4t,表1 中的測量結果也基本支持上述論斷.值得注意的是,在本算例中UCG-L 所需的計算時間還不及 UFG 的一半,但其數(shù)值結果卻與 UFG非常相近.有理由相信在某些特殊算例中通過合理布置加密區(qū)域,局部網格加密可以實現(xiàn)更高的加速比.

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

3.4 一維剪切波問題

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

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

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

針對該問題數(shù)值模擬的各參數(shù)和設置如下:上下使用周期邊界條件;粗細網格上均采用BGK 碰撞模型和D2Q9 速度集;非平衡態(tài)分布函數(shù)的初始化仍然使用式 (27).初始速度幅值u0對應的Mach數(shù)為 0.1;Reynolds 數(shù)Re=Hu0/ν=10,H為槽道寬度;網格數(shù)Nx=4,Ny=20.分別計算了加密區(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ù)值結果如圖12 所示.圖12 展示了在不同區(qū)域使用局部網格加密后,在給定時刻(t?=tν/H2=4π2/ln2) 下的數(shù)值黏性與物理黏性之比(紅色線)以及均勻粗網格下的數(shù)值黏性與物理黏性之比(藍色線),圖中的黑色虛線代表粗細網格的分界線.從圖中可以發(fā)現(xiàn),不同加密區(qū)域下,加密界面兩側的數(shù)值黏性相比于均勻網格下,表現(xiàn)出了不一樣的結果,既有增大的情況也有降低的情況.而其中數(shù)值黏性小于物理黏性現(xiàn)象的出現(xiàn)可能與加密界面附近的非物理震蕩[19,54-55]有著密切的關系,故還需要一系列更加系統(tǒng)的理論研究將局部網格加密技術引起的數(shù)值耗散和數(shù)值色散與加密區(qū)域附近的非物理震蕩聯(lián)系起來,而這也將作為我們后續(xù)工作的重點.

圖12 一維剪切波問題在不同加密區(qū)域下均勻網格和非均勻網格對應的數(shù)值黏性和物理黏性之比的分布情況(黑色虛線為加密界面).(a)～(d)分別為不同加密區(qū)域下UCG-L 和UCG 對應的數(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ū)域下均勻網格和非均勻網格對應的數(shù)值黏性和物理黏性之比的分布情況(黑色虛線為加密界面).(a)～(d)分別為不同加密區(qū)域下UCG-L 和UCG 對應的數(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 總結

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

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

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

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

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

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

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

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

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

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

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

進而式 (A1)可以表示為

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

若記Hermite 轉換矩陣MH為如下形式

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

D2Q9 速度集下的Hermite 轉換矩陣MH可以顯式表示為

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

當考慮源項為體力項時,|gneq〉 的前6 個矩的具體形式如下

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

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

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

進行簡單移項后

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

上式也可以寫作如下形式

歸納可得如下關系