平衡分布正值性對格子Boltzmann方法數(shù)值表現(xiàn)影響分析*

葉歡鋒 金頔 匡波 楊燕華2)

1) (上海交通大學(xué)核科學(xué)與工程學(xué)院,上海 200240)

2) (國家電投中央研究院核電軟件技術(shù)中心,北京 100029)

在傳統(tǒng)中,格子Boltzmann離散模型的數(shù)值表現(xiàn)影響因素分析主要關(guān)注于平衡分布的矩精度,而模型的正值性(作為粒子分布描述,分布函數(shù)需恒為正)則僅作為模型的附屬屬性,用于計算時工況約束.隨著部分高斯-厄米特求積公式離散方案的提出,模型正值性被發(fā)現(xiàn)是獨立于矩精度的模型屬性,可以通過格子速度調(diào)整.研究人員推測平衡分布正值性對格子Boltzmann方法的數(shù)值表現(xiàn)存在顯著影響,可以通過改善平衡分布正值性改善模型數(shù)值表現(xiàn).相比提升模型矩精度方案,正值性改善方案具有計算量優(yōu)勢.然而鑒于高階模型邊界處理的缺失,相關(guān)推測并未得到具體數(shù)值計算證實.本文采用周期邊界的Taylor-Green渦算例,回避了邊界處理問題,詳細(xì)分析了正值性對數(shù)值表現(xiàn)的影響,包括平衡分布正值區(qū)域內(nèi)計算精度穩(wěn)定性以及模型平衡分布正值區(qū)域大小對計算影響,并與矩精度影響進(jìn)行對比.計算結(jié)果顯示,模型正值區(qū)域內(nèi)計算精度并不恒定,隨著工況靠近正值區(qū)域上界,計算精度下降,但總體上均具備較好的精算精度.模型數(shù)值表現(xiàn)同時受到矩精度與模型正值性影響,矩精度影響主要體現(xiàn)在模型是否滿足Galilean不變性上,對于滿足Galilean不變性模型,其數(shù)值表現(xiàn)則取決于模型正值性.基于此,本文認(rèn)為通過改善模型正值性提升格子Boltzmann方法數(shù)值表現(xiàn)是切實可行的方案,并推薦基于滿足Galilean不變性條件下選擇具有最寬正值區(qū)域的模型,而不必執(zhí)著于模型矩精度.另外從本文的數(shù)值結(jié)果來看,高階模型模型的數(shù)值表現(xiàn)均好于經(jīng)典D2Q9模型,特別是D2H3-2模型,是文中涉及模型的最優(yōu)者,值得進(jìn)一步深入研究.總體而言,通過數(shù)值分析首次系統(tǒng)性梳理了模型正值性對計算影響,并與矩精度進(jìn)行對比分析.本文證實了正值性對計算的影響,為離散模型選擇和改進(jìn)提供了新的方向.

1 引 言

最近幾十年,格子Boltzmann方法在計算流體領(lǐng)域受到了廣泛的關(guān)注.作為一種介觀尺度計算方案,格子Boltzmann方法既能提供相較于宏觀計算方法更為深入的視角,又能避免純粒子尺度計算方法所需要的大量計算資源,因此在諸多傳統(tǒng)方法難以適用的復(fù)雜流體領(lǐng)域得到了廣泛應(yīng)用,如多組分多相流[1-4],微納米尺度流[5-7]、多孔介質(zhì)流[6,8]、粒子懸浮流[9,10]、血液流[11,12]、湍流[13-15]等.格子Boltzmann方法描述的是一群虛構(gòu)粒子在流場內(nèi)遷移和碰撞的過程.流場內(nèi)粒子信息通過粒子分布函數(shù)來描述,而流場的宏觀信息通過對粒子分布函數(shù)做相應(yīng)的速度矩積分得到.格子Boltzmann方法的控制方程可以通過對BGK近似下的Boltzmann方程沿特征線積分得到[16-18],而具體的算法實現(xiàn)則還需對得到的特征線積分在速度空間上進(jìn)行離散[17,18],如常見的速度離散模型D2Q9等.然而低階矩精度速度離散模型在還原宏觀控制方程式的過程中會引入誤差,如D2Q9模型還原動量守恒方程時會引入額外黏性項,導(dǎo)致方程不滿足Galilean不變性[19-21].為解決該問題,研究者們提出了高階速度離散模型,如Shan[22,23]和Shim[24]的工作.然而由于方案實現(xiàn)的復(fù)雜性,高階速度模型數(shù)量非常有限,導(dǎo)致對于高階模型數(shù)值表現(xiàn)的理解集中于模型的矩精度.

最近Ye等[25]提出了部分高斯-厄米特求積公式(partial Gauss-Hermite quadrature,pGHQ),顯著簡化了離散模型的構(gòu)建.Ye等[25]用該方案系統(tǒng)性地搜索了矩精度介于{3-7}離散速度在[-10,10]范圍內(nèi)的可行離散模型,結(jié)果顯示可行模型的數(shù)量非常豐富,遠(yuǎn)超預(yù)期,尤其是存在大量精度相同但是離散速度不同的模型.這些模型從宏觀角度分析是等價的,即用Chapman-Enskog展開可以還原出相同的宏觀方程.但是在微觀上,模型有著明顯的區(qū)別,特別是平衡分布的正值區(qū)域(所有分布為正值的宏觀速度區(qū)域).為方便討論,本文用正值性指代分布正值相關(guān)性質(zhì)包括正值區(qū)域等,且在后續(xù)討論中不區(qū)分兩者使用.相對于傳統(tǒng)理解中作為模型用于計算時工況的約束條件,在pGHQ方案下,平衡分布正值性可以看作是獨立于矩精度的模型屬性,可以通過離散速度調(diào)整.鑒于負(fù)值平衡分布違背格子Boltzmann方法的物理原則,即粒子分布不存在負(fù)值,Ye等[25]認(rèn)為離散模型的正值性將對模型的數(shù)值表現(xiàn)有直接影響,可以通過選擇模型離散速度來改善計算.對比通過提升模型矩精度改善格子Boltzmann方法數(shù)值表現(xiàn)的傳統(tǒng)構(gòu)想,調(diào)整模型離散速度設(shè)想在計算量上存在優(yōu)勢.然而鑒于高階模型的具體實現(xiàn)需要設(shè)計對應(yīng)的邊界處理,Ye等[25]文章中并沒給出具體數(shù)值驗證.本文選擇可以采用周期性邊界處理的Taylor-Green渦算例,回避了高階離散模型邊界處理匱乏的問題.通過詳細(xì)的數(shù)值實驗,文章分析了平衡分布正值性對離散模型數(shù)值穩(wěn)定性和精度的影響,并與模型矩精度影響進(jìn)行對比,驗證通過提升模型正值性來改善格子Boltzmann方法數(shù)值表現(xiàn)的可行性.

本文首先回顧Ye等[25]提出的pGHQ離散模型構(gòu)造方案,并選擇了部分離散模型用作具體計算.然后本文用選擇的模型計算Taylor-Green渦,分析平衡分布正值性對格子Boltzmann方法計算穩(wěn)定性和精度的影響,并與模型矩精度影響進(jìn)行對比,驗證通過提升模型正值性來改善格子Boltzmann方法數(shù)值表現(xiàn)的可行性.最后總結(jié)本文工作.

2 pGHQ速度離散方案

pGHQ速度離散方法基于Hermite展開理論[26],即平衡分布可以通過對Maxwell分布以Hermite多項式形式展開構(gòu)造.以一維為例,對于m階精度的平衡分布,其分布函數(shù)可以通過如下方式構(gòu)造：

式中 ρ 和u分別為密度和宏觀流速,vα和 wα為離散格子速度和對應(yīng)的權(quán)重系數(shù),在Hermite展開理論中被稱為無量綱網(wǎng)格常數(shù)c,而 Hi(x) 則為i階物理Hermite多項式,表1羅列了前5階物理Hermite多項式的具體表達(dá).無量綱網(wǎng)格常數(shù)c,離散格子速度 vα和對應(yīng)的權(quán)重系數(shù) wα是通過求解對應(yīng)的 0-2m 階數(shù)值積分方程得到,

式中 ξα=vαc.方程(2)是高階非線性方程組,為減少其求解難度,Shan[22,23]和Shim[24]預(yù)先定義了對稱格子速度集 {0,±v1,···} ,使得方程組轉(zhuǎn)變?yōu)閮H包含無量綱網(wǎng)格常數(shù)c和權(quán)重系數(shù) {wα} 未知量的方程,并假定m階矩精度平衡分布構(gòu)造需要2m-1個格子速度.至此離散速度模型構(gòu)造轉(zhuǎn)變?yōu)闃?gòu)造 {0,±v1,···,±vm-1} 并求解其對應(yīng)方程組(2).盡管Shan和Shim等通過預(yù)定義格子速度集{0,±v1,···,±vm-1}顯著簡化了離散速度模型構(gòu)造,但是求解代入格子速度集后的方程組(2)依舊繁瑣而復(fù)雜,無法大規(guī)模驗證可行的格子速度集.因此現(xiàn)有的高階可行離散模型數(shù)量十分有限,這阻礙了對于格子Boltzmann速度離散模型的進(jìn)一步研究.另外對于模型矩精度與格子速度數(shù)量之間的關(guān)系也沒經(jīng)過系統(tǒng)而深入地分析.因此在Shan和Shim的研究中,其重點更偏重于突破高階離散模型的缺失,并未形成高階模型系統(tǒng)性認(rèn)識.這使得對于高階模型的認(rèn)知依舊局限于傳統(tǒng)“模型階數(shù)決定模型數(shù)值表現(xiàn)”的猜測,即模型階數(shù)越高數(shù)值計算表現(xiàn)越好.

表1 前五階物理Hermite多項式Table 1.The first five physicists′ Hermite polynomials.

Ye等[25]通過改造Gauss-Hermite求積公式,給出了方程(2)的等價求積公式,即pGHQ求積公式.根據(jù)pGHQ求積公式,互不相同的q節(jié)點ξα和其權(quán)重系數(shù) wα若要滿足方程組(2),只需保證其節(jié)點多項式 Wq(ξ) 基于Hermite多項式表述時,僅涉及 2m-q+1 到q階之間的Hermite多項式,

而對應(yīng)的權(quán)重系數(shù)計算公式則為

完整的推導(dǎo)可參考文獻(xiàn)[25].不同于(2)式同時包含節(jié)點 ξα和權(quán)重系數(shù) wα,(3)式僅為節(jié)點 ξα的函數(shù).因此當(dāng)給定格子速度集 {vα} 時,可以根據(jù)(3)式直接得到網(wǎng)格常數(shù)c的約束方程,其具體流程是將格子速度集代入其節(jié)點多項式并展開多項式,

其中 bi是網(wǎng)格常數(shù)c的多項式,引入次方項 ξi與Hermite多項式的關(guān)系式,

則(5)式可以表達(dá)為

而根據(jù)pGHQ (3)式,要使得格子速度滿足m階矩精度離散模型構(gòu)造,其節(jié)點多項式基于Hermite多項式表述時,應(yīng)當(dāng)僅涉及 2m-q+1 到q階之間的Hermite多項式,即 0 到 2m-q 階Hermite多項式系數(shù)為 0 ,

而系數(shù) ai和 bi一樣,是僅包含網(wǎng)格常數(shù)c的多項式.因此(8)式事實上是給定格子速度集 {vα} 的網(wǎng)格常數(shù)c約束方程組.求解(8)式,若方程組有解,則將求解的網(wǎng)格常數(shù)以及格子速度代入(4)式,即可得到對應(yīng)的權(quán)重系數(shù).將網(wǎng)格常數(shù),格子速度以及權(quán)重系數(shù)代入(1)式,即可得到m階矩精度離散模型.

在Hermite展開理論中,高維模型可以通過對一維模型做張量乘積得到.以二維為例,其平衡函數(shù)構(gòu)造即為

對應(yīng)的格子速度為

Shan[22,23]和Shim[24]的離散方案,pGHQ速度離散方法顯著簡化了離散模型構(gòu)建過程.其算法的簡潔和通用性讓系統(tǒng)性研究格子速度集成為可能.Ye等[25]用該方案系統(tǒng)性地搜索了矩精度分別為{3-7}，離散速度在 [-10,10]范圍內(nèi)的可行離散模型,得到了大量等矩精度的離散模型.這些模型矩精度相同,但是正值性存在顯著差別.鑒于負(fù)值平衡分布違背格子Boltzmann方法的物理原則,即粒子分布不存在負(fù)值,Ye等[25]認(rèn)為模型正值性會影響模型的數(shù)值表現(xiàn),并推測可以通過選擇格子速度集改善模型正值性的方式改善格子Boltzmann方法的計算穩(wěn)定性和精度.

為驗證該推論的可行性,本文選擇了表2所列的離散模型,模型以維度n和速度矩精度m來表示,即DnHm.文章主要從如下幾個方面驗證正值性影響：1)平衡分布正值性對于數(shù)值表現(xiàn)影響,是否在正值區(qū)域內(nèi),模型都能保持計算的穩(wěn)定性和相應(yīng)的精度; 2) 平衡分布正值性與模型矩精度對數(shù)值表現(xiàn)影響對比; 3)通過調(diào)整平衡分布正值性改善數(shù)值表現(xiàn)可行性.鑒于非對稱網(wǎng)格對于平衡分布的正值性影響主要在于平移正值區(qū)域,其在理論上更適合具有主流速度的計算問題.而本文采用的Taylor-Green渦流場是對稱性,其主流速度為0,因此并不適合針對于非對稱網(wǎng)格的分析.而從本文研究目標(biāo)來看,單純采用對稱離散模型足夠確認(rèn)相關(guān)結(jié)果.因此本文在模型選擇上并不涉及非對稱離散模型.同時為驗證Hermite展開理論的實用性,本文還額外采用了經(jīng)典D2Q9模型作為比較,

表2 格子Boltzmann離散模型.為方便展示,表格羅列的是用于對應(yīng)構(gòu)造高階模型的一維模型參數(shù).實際計算中需根據(jù)(9)式和(10)式對表格中的參數(shù)進(jìn)行張量構(gòu)造.另外需注意的是表格中羅列的是網(wǎng)格常數(shù)c,其與格子聲速 cs 存在如下?lián)Q算關(guān)系c=1/csTable 2.Discrete model description of lattice Boltzmann method.For the seeking of space saving,the parameters illustrated in this table are of the corresponding unidimensional models.In numerical implementations,they require tensor product illuminated in Eq.(9) and Eq.(10).It should be noted that the table lists the lattice constant c instead of the lattice sonic speed cs ,which can be expressed as c=1/cs.

表2 格子Boltzmann離散模型.為方便展示,表格羅列的是用于對應(yīng)構(gòu)造高階模型的一維模型參數(shù).實際計算中需根據(jù)(9)式和(10)式對表格中的參數(shù)進(jìn)行張量構(gòu)造.另外需注意的是表格中羅列的是網(wǎng)格常數(shù)c,其與格子聲速 cs 存在如下?lián)Q算關(guān)系c=1/csTable 2.Discrete model description of lattice Boltzmann method.For the seeking of space saving,the parameters illustrated in this table are of the corresponding unidimensional models.In numerical implementations,they require tensor product illuminated in Eq.(9) and Eq.(10).It should be noted that the table lists the lattice constant c instead of the lattice sonic speed cs ,which can be expressed as c=1/cs.

Model nameimages/BZ_182_878_603_916_628.pngDiscrete velocityset { } Lattice constant cimages/BZ_182_1811_603_1857_628.png{images/BZ_182_1546_658_1768_692.png ,1.6667images/BZ_182_1886_658_1987_692.png }D2H3-1 Weights { }D2H2{ }images/BZ_182_711_662_790_696.pngimages/BZ_182_1070_658_1262_691.png{images/BZ_182_1408_720_1630_754.png ,images/BZ_182_1644_720_1877_758.png ,images/BZ_182_1892_720_2125_758.png }D2H3-2{ }images/BZ_182_679_724_821_758.pngimages/BZ_182_1058_720_1274_753.png{ }images/BZ_182_681_787_819_820.pngimages/BZ_182_1058_782_1274_816.png{images/BZ_182_1408_783_1630_817.png ,images/BZ_182_1644_783_1877_820.png ,images/BZ_182_1892_783_2125_820.png }D2H4{images/BZ_182_1412_837_1633_871.png ,images/BZ_182_1648_837_1881_874.png ,,}D2H5{ }images/BZ_182_650_862_850_895.pngimages/BZ_182_1058_857_1274_891.pngimages/BZ_182_1912_837_2128_870.pngimages/BZ_182_1650_879_1866_912.png{images/BZ_182_1412_925_1633_959.png ,images/BZ_182_1648_925_1881_962.png ,,,images/BZ_182_1757_967_1990_1004.png }D2H6{ }images/BZ_182_619_950_881_983.pngimages/BZ_182_1058_945_1274_979.pngimages/BZ_182_1912_925_2128_958.pngimages/BZ_182_1526_967_1742_1000.png{ }images/BZ_182_590_1042_911_1075.pngimages/BZ_182_1058_1037_1274_1071.png{images/BZ_182_1412_1013_1633_1047.png ,images/BZ_182_1648_1013_1881_1050.png ,,,,images/BZ_182_1881_1055_2114_1092.png }images/BZ_182_1912_1013_2128_1046.pngimages/BZ_182_1402_1055_1619_1088.pngimages/BZ_182_1650_1055_1866_1088.png

其中格子速度,權(quán)重系數(shù)以及格子聲速 cs的取值與D2H2模型相同.

對于DnHm類型的模型,其正值區(qū)域是比較容易分析的,即只需分析對應(yīng)一維模型的所有分布保持正值性的流體速度區(qū)域,高維模型的正值區(qū)域即為該區(qū)域的張量乘積結(jié)果.對于本文涉及的{0,±v1,···,±vm-1}類型一維模型,其正值區(qū)域具有對稱性.經(jīng)過張量乘積后,該對稱性在高維下依舊得到保留,因此本文以一維正值區(qū)域上限來指代模型正值區(qū)域,符號記為 Upos.以本文D2H2模型對應(yīng)的一維模型為例,其平衡分布在 [-0.82,0.82]之間均為正值,超過該范圍后則分布出現(xiàn)負(fù)值.而D2H2模型是一維模型的張量乘積結(jié)果,因此流體速度在[-0.82,0.82]×[-0.82,0.82]區(qū)間內(nèi),D2H2模型均為正值,即D2H2模型正值區(qū)域為[ -0.82,0.82]×[-0.82,0.82],其 Upos=0.82.通過同樣分析可以得到本文其他DnHm型的正值區(qū)域.比較特殊的是D2Q9模型,不同于DnHm類型模型的正值區(qū)域為一維張量乘積結(jié)果,D2Q9模型下不同方向的流體速度是相互關(guān)聯(lián)的.為方便分析,取D2Q9模型正值區(qū)域中最大上滿足D2Q9正值性的最大 Upos值.本文涉及模型的正值區(qū)域如圖1所示.

3 數(shù)值分析

兩維Taylor-Green渦作為等溫不可壓Navier-Stokes方程組的有限幾組精確解,是流體計算中常用的經(jīng)典算例,其控制方程為：

其解為

式中 ρ0為流體常物性密度; p0為參考壓力; u0為Taylor-Green渦特征速度; kx,ky為周期因子;Td為衰減特征時間，其計算式為1/Td=為方便討論,本文將周期因子恒定選擇為 kx=ky=1.根據(jù)兩維Taylor-Green渦的解,當(dāng)流場初始速度和壓力分布滿足：

圖1 格子Boltzmann離散模型正值性分析 虛線框表示模型正值區(qū)域,其涵蓋區(qū)域為 [-Upos,Upos]×[-Upos,Upos],模型 Upos 標(biāo)注在圖中模型名稱之后Fig.1.Positivity analysis of lattice Boltzmann discrete models.The dash squares denote the positive areas of each lattice Boltzmann discrete model,which are constructed as [-Upos,Upos]×[-Upos,Upos].The Upos values of each model are annotated after the model names.

其速度場和壓力場分別以指數(shù) exp(-t/Td) 和exp(-2t/Td)衰減,且速度場和壓力場在流場內(nèi)分別以 2π 和 π 為周期成周期分布.Taylor-Green渦解十分有標(biāo)識性,統(tǒng)一的指數(shù)衰減形式非常適合用于測試算法的穩(wěn)定性.而周期分布避免了邊界處理,適用于驗證目前邊界處理缺乏的高階離散模型.

本文用格子Boltzmann計算了流域為[0,2π]×[0,2π]的Taylor-Green渦演變.流體的運動黏性ν為 0.1m2/s ,特征速度 u0為 1m/s.網(wǎng)格密度為200×200.模擬采用單松弛因子格子Boltzmann方程,

其中 Δt 和 ΔL 分別為計算時間步長和網(wǎng)格長度,ν為流體運動黏性,cs為格子聲速.流場的初始分布為：

其中 uLB,0為 u0在格子Boltzmann體系下的值,uLB,0=u0Δt/ΔL.不同于初始流場引入,初始壓力場的引入需經(jīng)過轉(zhuǎn)換處理.本文根據(jù)壓力和密度關(guān)系式通過設(shè)置初始密度場來引入壓力初始條件,

為避免初始密度出現(xiàn)負(fù)值,本文統(tǒng)一取 ρp0=10 ,ρ0=1.從Taylor-Green渦的設(shè)置可以看出,要驗證平衡分布正值性影響,可以通過調(diào)整 uLB,0來實現(xiàn),即通過調(diào)整計算步長來改變 u0在格子Boltzmann體系下的值,人為控制流場初始時刻平衡分布的正值性.根據(jù)模型正值區(qū)域分析結(jié)果(圖1),本文設(shè)計了如表3中所列的計算工況.

根據(jù)Taylor-Green渦的初始條件設(shè)置,隨著uLB,0增大,格子Boltzmann的初始分布將開始出現(xiàn)負(fù)值,而具體的負(fù)值情況則取決于平衡分布模型.為方便分析和展現(xiàn)具體分布正值情況,本文設(shè)計了如下正值指標(biāo)：

表3 算例設(shè)計Table 3.Case design for lattice Boltzmann method.

當(dāng)分布中存在負(fù)值時,κ 值將大于1,而不存在負(fù)值時則 κ 值恒為1.κ 偏離1的幅值可以看作分布正值性破壞的嚴(yán)重程度.該指標(biāo)通過對 ρ 歸一,消除了Taylor-Green渦流場內(nèi)不同密度對正值性的影響.圖2給出了設(shè)計計算工況下,各模型初始分布的 κ 值情況.圖中額外標(biāo)注了對應(yīng)工況下模型正值區(qū)域值 upos與設(shè)計 uLB,0的差值,記作 Δu ,

用以輔助理解設(shè)計計算工況對模型正值區(qū)域的偏離程度.Δu 值為正,表示計算工況完全位于模型正值區(qū)域內(nèi),Δu 值為負(fù)則表示計算工況超過正值區(qū)域,部分網(wǎng)格出現(xiàn)負(fù)值分布.另外需注意的是,圖2對 κ 值的范圍做了截斷處理,即僅區(qū)別[1.0,2.0]范圍內(nèi)變化,以方便觀察分布正值性隨設(shè)計工況而遭到破壞的變化過程.κ 值圖顯示,當(dāng) Δu 值為負(fù)時,分布的正值性遭到破壞,κ 值大于1,如D2H2模型在b,c,d工況情況下所示.但是對于給定負(fù)值 Δu ,平衡分布正值性遭到的破壞程度取決于具體模型.以D2H2和D2H6模型為例,D2H2在 Δu 值為 -0.18 時,初始 κ 分布出現(xiàn)了明顯偏離1值; 而D2H6模型在 Δu 值為 -1.02 時,初始 κ 分布偏離 1 值的幅度依舊小于前述情況.

對于模型數(shù)值表現(xiàn)分析,本文主要考量計算的流場分布以及誤差演化.流場分布分析通過比較Taylor-Green渦初始流場強(qiáng)度衰減一半時(t=3.4657359s)模型計算結(jié)果與理論解的速度平方和U2等高線圖實現(xiàn),

鑒于速度平方和 U2在流域內(nèi)成周期性分布,本文僅截取周期區(qū)域[0.3775,0.6225]π×[0.3775,0.6225]π內(nèi)圖像(見圖3).誤差演化分析則關(guān)注模型計算結(jié)果與理論解誤差在 t=3.4657359s 時間內(nèi)隨時間的演變(圖4),其計算公式為

式中 ua和 va為理論結(jié)果.

圖2 各格子Boltzmann模型在Taylor-Green渦設(shè)計計算工況下初始時刻計算流域 [0,2π]×[0,2π]內(nèi) κ 值色階圖 圖中列對應(yīng)格子Boltzmann模型,行則對應(yīng)表3中計算工況; 每個色階圖上方的數(shù)值 Δu 值,即正值區(qū)域值 upos 設(shè)計計算工況 uLB,0 差值;為方便觀察正值性隨計算工況遭到破壞的變化過程,圖中對 κ 值做了截斷處理,僅區(qū)別 [1.0,2.0]范圍變化Fig.2.The initial κ filled contours of designed Taylor-Green vortex simulations for each lattice Boltzmann model in the fluid domain [0,2π]×[0,2π].Each column is a simulation set of a lattice Boltzmann model with different designed cases.Each row is a set of a simulation case with different lattice Boltzmann models.The values in the contour panels are the values of Δu ,i.e.the difference between the upos and the designed uLB,0.For the sake of identifying the positivity violation process of a lattice Boltzmann model,the filled contours truncate the range of κ ,limiting into [1.0,2.0].

圖3 Taylor-Green渦半衰時刻 t=3.4657359s ,[0.3775,0.6225]π×[0.3775,0.6225]π流域內(nèi),不同計算工況下理論解(Theoretical)與各格子Boltzmann模型的速度平方和 U2 等高線圖; 圖中列對應(yīng)格子Boltzmann模型,行則對應(yīng)表3中計算工況Fig.3.The U2 contour lines of Taylor-Green vortex in fluid domain [0.3775,0.6225]π×[0.3775,0.6225]πat half-value decay time t=3.4657359s,including the theoretical solution (denoted as “Theoretical”) and the numerical results of lattice Boltzmann models.Each column is a set of simulations under the designed cases.Each row is a set of simulations under selected lattice Boltzmann models including the theoretical solution.

圖4 各Taylor-Green渦設(shè)計計算工況下,各格子Boltzmann模型計算誤差隨時間演變,圖中標(biāo)號a,b,c,d分別對應(yīng)表3各設(shè)計計算工況; 為方便表示,圖中橫坐標(biāo)單位為0.173286795sFig.4.The error evolution of each simulation.The panel a,b,c,d renders each designed case in Table 3 respectively.For the sake of rendering,the unit of the time axis is scaled as 0.173286795 s.

從計算結(jié)果來看,對于每個DnHm模型,其數(shù)值精度隨著 uLB,0增大而變壞,無論是計算結(jié)果等高線與理論解的相似程度(圖3),還是誤差演化(圖4)都體現(xiàn)了該趨勢.但與模型正值性一受到破壞計算即發(fā)散的D2H2模型不同的是,高階DnHm模型對于模型正值性破壞并未體現(xiàn)出間斷式的數(shù)值表現(xiàn)變化.對于單個計算工況,DnHm系列模型除不滿足Galilean不變性的D2H2模型外,計算精度隨 upos減小而變壞.而D2H2模型精度則差于所有高階模型,包括正值區(qū)域更小的D2H4模型.對比D2H2模型與其他高階模型的數(shù)值表現(xiàn),模型的Galilean不變性對于模型的計算穩(wěn)定性具有顯著影響.從這一角度而言,矩精度對模型數(shù)值表現(xiàn)影響主要體現(xiàn)在模型是否滿足Galilean不變性上.對于滿足Galilean不變性的模型,矩精度影響不再顯著,模型數(shù)值表現(xiàn)主要受正值性影響.

圖5 Taylor-Green渦設(shè)計計算工況下,Δu 值為正值的各格子Boltzmann模型數(shù)值模擬誤差對比.圖中實線,虛線和虛點線分別為工況a,b,c (見表3)結(jié)果,離散模型則用顏色和符號標(biāo)注.為方便表示,圖中橫坐標(biāo)單位為0.173286795 sFig.5.The error comparison of Taylor-Green vortex simulations with postive value of Δu.The designed configurations of Taylor-Green vortex are denoted with line styles,in which the solid,dashed and dash-dotted line indicates case a,b and c in Table 3 respectively,meanwhile the lattice Boltzmann models are labeled with line colors and markers.For the sake of rendering,the unit of the time axis is scaled as 0.173286795 s.

分析模型正值區(qū)域內(nèi)計算表現(xiàn),即 Δu 大于0的數(shù)值模擬,所有計算結(jié)果等高線與理論解均保持較好的符合,而對比誤差(圖5),計算結(jié)果未明顯表現(xiàn)出模型正值區(qū)域與精度之間跨工況相關(guān)性.具體而言,在模型正值區(qū)域內(nèi),具有較寬正值區(qū)域的模型并不對計算恒定保持優(yōu)勢精度,以D2H3-2模型工況b為例,其計算差于所有其他3階及3階以上DnHm模型工況a計算結(jié)果.進(jìn)一步而言,在模型正值區(qū)域內(nèi),其計算精度并不是恒定的,計算工況越靠近正值區(qū)域上限,模型精度越低.因此模型正值區(qū)域優(yōu)勢主要體現(xiàn)在同一計算工況對比.

在數(shù)值計算中,獨立的精度指標(biāo)是非常實用的指示參量,即對于任意計算工況,只需保證該精度指標(biāo)一致,則計算精度一致.鑒于正值區(qū)域的精度優(yōu)勢主要體現(xiàn)在同工況下模型之間的對比,無法表征數(shù)值模擬的確切精度,本文嘗試探索其他參量作為精度指標(biāo)的可能性.在Taylor-Green渦算例中,文章設(shè)計了 Δu 值和 κ 值來表征模型在對應(yīng)計算工況下的正值性情況.本文探索這兩指標(biāo)作為衡量計算精度的可能性.文章首先檢驗 Δu 值,以D2H4模型工況a,b,c與D2H5模型工況b,c,d為例,兩兩分別具有近似的 Δu 值.若格子Boltzmann方法的數(shù)值表現(xiàn)與 Δu 值強(qiáng)相關(guān),則D2H4模型與D2H5模型對應(yīng)算例精度將呈現(xiàn)明顯的相關(guān)性.從圖6 (a)來看,對比D2H3-2模型隨工況變化而導(dǎo)致精度變化程度,對應(yīng)D2H4模型與D2H5模型結(jié)果之間的精度差異幅度依舊明顯,并不足以體現(xiàn)格子Boltzmann方法計算精度與 Δu 值相關(guān)性.另外取 Δu 值在 (-0.30±0.6) 附近計算結(jié)果 (圖6 (b)),同樣對比D2H3-2模型隨工況變化而導(dǎo)致精度變化程度,選取的計算精度并未收斂于足夠識別Δu值影響的值域附近.從上述分析可以看出,Δu 值并不具備衡量格子Boltzmann方法數(shù)值精度的能力.對于 κ 值,本文選取了模型正值性已遭受破壞但初始 κ 值尚未體現(xiàn)顯著偏離的計算結(jié)果 (圖7).類似地,對比D2H3-2模型隨工況變化而導(dǎo)致精度變化程度,誤差對比顯示,選取算例并未顯示顯著的收斂性,因此 κ 值也不具備衡量格子Boltzmann方法數(shù)值表現(xiàn)的能力.

圖6 Taylor-Green渦數(shù)值模擬誤差與 Δu 值相關(guān)性分析 (a)對比了 Δu 接近的D2H4模型a,b,c工況與D2H5模型b,c,d工況計算誤差; (b)對比了 Δu 在 -0.30 附近所有數(shù)值模擬計算誤差.所有圖中均保留了D2H3-2模型在a,b,c,d工況下的計算誤差作為參照.曲線上標(biāo)注值為該數(shù)值模擬的 Δu 值.圖中a,b,c,d工況分別用實線,虛線,虛點線及點線標(biāo)注,而離散模型則用顏色和符號標(biāo)注.為方便表示,圖中橫坐標(biāo)單位為0.173286795 sFig.6.The numerical performance of Taylor-Green vortex simulations vs the value of Δu.Panel a plots the numerical errors of model D2H4 under case a,b,c and model D2H5 under case b,c,d,which possess close value of Δu respectively.Panel b plots the numerical errors of simulations with a value of Δu around -0.30.The numbers labeled on the curves are their values of Δu.The simulation configurations are denoted with line style,in which solid,dashed,dash-dotted and dotted line indicates case a,b,c and d respectively,meanwhile the lattice Boltzmann models are labeled with line colors and markers.In all panels,the results of model D2H3-2 with four designed configurations are plotted for a reference.The designed configurations of Taylor-Green vortex are denoted with line styles,in which the solid,dashed,dash-dotted and dotted line indicates case a,b,c and d in Table 3 respectively,meanwhile the lattice Boltzmann models are labeled with line colors and markers.For the sake of rendering,the unit of the time axis is scaled as 0.173286795 s.

圖7 Taylor-Green渦數(shù)值模擬誤差與初始 κ 值相關(guān)性分析.圖中對比了在 Δu 值為負(fù)情況下,初始 κ 值偏離 1 幅值可忽略的算例計算誤差,包括模型D2H3-1工況c,模型D2H3-2工況d,模型D2H4工況b,模型D2H5工況c,d和模型D2H6工況b,c.圖中均保留了D2H3-2模型在a,b,c,d工況下的計算誤差作為參照.曲線上標(biāo)注值為該數(shù)值模擬的 Δu 值.圖中a,b,c,d工況分別用實線,虛線,虛點線及點線標(biāo)注,而離散模型則用顏色和符號標(biāo)注.為方便表示,圖中橫坐標(biāo)單位為0.173286795 sFig.7.The numerical performance of Taylor-Green vortex simulations vs initial κ.The error evolutions of numerical simulations with negative values of Δu but negligible departure of initial κ from 1 are plotted,including model D2H3-1 under case c,model D2H3-2 under case d,model D2H4 under case b,model D2H5 under case c,d and model D2H6 under case b,c.The results of model D2H3-2 with four designed configurations are also plotted for a reference.The designed configurations of Taylor-Green vortex are denoted with line styles,in which the solid,dashed and dashdotted line indicates case a,b and c in Table 3 respectively,meanwhile the lattice Boltzmann models are labeled with line colors and markers.For the sake of rendering,the unit of the time axis is scaled as 0.173286795 s.

對比經(jīng)典的D2Q9模型,在正值區(qū)域內(nèi),D2H2模型結(jié)果與之驚人的一致,如工況a下兩者的等高線圖和誤差演變曲線所示.但在穩(wěn)定性上,D2H2模型則略差于D2Q9模型.而本文涉及的所有高階模型(矩精度大于等于3)計算穩(wěn)定性和精度都優(yōu)于D2Q9模型.從本文數(shù)值計算結(jié)果來看,D2H3-2模型具有最佳的計算穩(wěn)定和精度,同時對比其他高階DnHm系列模型,其計算量也具有明顯優(yōu)勢.因此從數(shù)值計算來說,D2H3-2模型值得進(jìn)一步研究.

4 結(jié) 論

本文用Taylor-Green渦回避了高階DnHm系列離散模型邊界處理缺乏的現(xiàn)狀,數(shù)值驗證了Ye等[25]提出的離散模型的精度受到模型的正值性影響的推測.文章選擇了矩精度為{2-6}的最緊湊對稱離散模型,以及具有較寬正值區(qū)間的D2H3-2模型作為分析對象,并與經(jīng)典D2Q9模型進(jìn)行比較.根據(jù)涉及離散模型的正值區(qū)域,文章通過調(diào)整時間步長設(shè)計了橫跨各模型正值區(qū)域與非正值區(qū)域的四個計算工況.分析各模型對四個計算工況數(shù)值模擬的結(jié)果,文章得出了如下結(jié)論.

1) DnHm系列在正值區(qū)域內(nèi),均能維持不錯的計算結(jié)果,但是并不能保證統(tǒng)一精度,即隨著計算工況靠近正值區(qū)域上限計算精度會隨之下降.在非正值區(qū)域計算時,模型的表現(xiàn)受模型是否滿足Galilean不變性影響.具備Galilean不變性的模型,其在非正值區(qū)域上也具有一定的計算穩(wěn)定性,并未出現(xiàn)間斷式數(shù)值表現(xiàn)變化.而不具備Galilean不變性的模型,其計算精度在非正值區(qū)域存在間斷式變化,即非正值區(qū)域計算直接發(fā)散.

2)對于DnHm系列模型,矩精度的影響主要體現(xiàn)在是否滿足Galilean不變性.對于不滿足Galilean不變性模型,除了其在非正值區(qū)域數(shù)值精度存在間斷式變化外,其正值區(qū)域內(nèi)精度也差于正值區(qū)域更小的高階模型.而對于滿足Galilean不變性的模型,矩精度影響不再顯著,其數(shù)值表現(xiàn)主要取決于正值區(qū)域.

3)對于滿足Galilean不變性的模型,在相同計算工況下,模型精度與模型正值區(qū)域正相關(guān),即模型正值區(qū)域越寬計算結(jié)果精度越好.鑒于模型正值性僅適合分析相同計算工況模型的相對數(shù)值精度表現(xiàn),本文探索了 Δu 值與 κ 值用于衡量格子Boltzmann方法數(shù)值表現(xiàn)的可能性,從分析結(jié)果來看,兩者并不具備表征能力.

4) D2H2模型與經(jīng)典的D2Q9模型有著非常一致的精度表現(xiàn),盡管穩(wěn)定性略有不如.而高階DnHm模型在穩(wěn)定性和計算精度上都優(yōu)于D2Q9模型.特別是D2H3-2模型,在所有本文涉及的DnHm模型中,具有最佳的穩(wěn)定性和精度,是值得深入研究的模型.

平衡分布正值性對計算的影響可以從粒子分布需保持正值得到解釋：離散模型對于Maxwell分布全局正值性的破壞導(dǎo)致模型僅在有限范圍有效,該范圍影響著模型的計算表現(xiàn).需要指出的是,該有限范圍并不等同于模型正值區(qū)域,其原因是離散模型在正值區(qū)域內(nèi)也僅為Maxwell分布的近似.因此計算中,模型正值區(qū)域內(nèi)的精度并不一致.而正值性與矩精度影響對比分析顯示,三階Hermite多項式展開近似即能很好近似Maxwell分布函數(shù)形狀.另外,本文對于離散模型的分析并不考慮模型間Ma數(shù)一致性問題,即計算工況在不同模型下Ma數(shù)并不相同.對于需要嚴(yán)格保持Ma數(shù)一致性的計算問題,相關(guān)結(jié)論還需進(jìn)一步驗證分析,本文在此不過多展開.

總結(jié)而言,本文通過數(shù)值分析首次系統(tǒng)性梳理了離散模型正值性對計算影響,并與矩精度進(jìn)行對比分析.文章證實了正值性對計算的影響,確認(rèn)通過提升模型正值性改善計算的可行性,為pGHQ離散方案中可行方案篩選與未來模型進(jìn)一步改進(jìn)提供了方向.