平衡分布正值性對格子Boltzmann方法數值表現影響分析*

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

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

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

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

1 引 言

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

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

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

2 pGHQ速度離散方案

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

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

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

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

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

而對應的權重系數計算公式則為

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

其中 bi是網格常數c的多項式,引入次方項 ξi與Hermite多項式的關系式,

則(5)式可以表達為

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

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

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

對應的格子速度為

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

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

表2 格子Boltzmann離散模型.為方便展示,表格羅列的是用于對應構造高階模型的一維模型參數.實際計算中需根據(9)式和(10)式對表格中的參數進行張量構造.另外需注意的是表格中羅列的是網格常數c,其與格子聲速 cs 存在如下換算關系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離散模型.為方便展示,表格羅列的是用于對應構造高階模型的一維模型參數.實際計算中需根據(9)式和(10)式對表格中的參數進行張量構造.另外需注意的是表格中羅列的是網格常數c,其與格子聲速 cs 存在如下換算關系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

其中格子速度,權重系數以及格子聲速 cs的取值與D2H2模型相同.

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

3 數值分析

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

其解為

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

圖1 格子Boltzmann離散模型正值性分析 虛線框表示模型正值區域,其涵蓋區域為 [-Upos,Upos]×[-Upos,Upos],模型 Upos 標注在圖中模型名稱之后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.

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

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

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

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

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

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

表3 算例設計Table 3.Case design for lattice Boltzmann method.

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

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

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

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

式中 ua和 va為理論結果.

圖2 各格子Boltzmann模型在Taylor-Green渦設計計算工況下初始時刻計算流域 [0,2π]×[0,2π]內 κ 值色階圖 圖中列對應格子Boltzmann模型,行則對應表3中計算工況; 每個色階圖上方的數值 Δu 值,即正值區域值 upos 設計計算工況 uLB,0 差值;為方便觀察正值性隨計算工況遭到破壞的變化過程,圖中對 κ 值做了截斷處理,僅區別 [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]π流域內,不同計算工況下理論解(Theoretical)與各格子Boltzmann模型的速度平方和 U2 等高線圖; 圖中列對應格子Boltzmann模型,行則對應表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渦設計計算工況下,各格子Boltzmann模型計算誤差隨時間演變,圖中標號a,b,c,d分別對應表3各設計計算工況; 為方便表示,圖中橫坐標單位為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.

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

圖5 Taylor-Green渦設計計算工況下,Δu 值為正值的各格子Boltzmann模型數值模擬誤差對比.圖中實線,虛線和虛點線分別為工況a,b,c (見表3)結果,離散模型則用顏色和符號標注.為方便表示,圖中橫坐標單位為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.

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

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

圖6 Taylor-Green渦數值模擬誤差與 Δu 值相關性分析 (a)對比了 Δu 接近的D2H4模型a,b,c工況與D2H5模型b,c,d工況計算誤差; (b)對比了 Δu 在 -0.30 附近所有數值模擬計算誤差.所有圖中均保留了D2H3-2模型在a,b,c,d工況下的計算誤差作為參照.曲線上標注值為該數值模擬的 Δu 值.圖中a,b,c,d工況分別用實線,虛線,虛點線及點線標注,而離散模型則用顏色和符號標注.為方便表示,圖中橫坐標單位為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渦數值模擬誤差與初始 κ 值相關性分析.圖中對比了在 Δu 值為負情況下,初始 κ 值偏離 1 幅值可忽略的算例計算誤差,包括模型D2H3-1工況c,模型D2H3-2工況d,模型D2H4工況b,模型D2H5工況c,d和模型D2H6工況b,c.圖中均保留了D2H3-2模型在a,b,c,d工況下的計算誤差作為參照.曲線上標注值為該數值模擬的 Δu 值.圖中a,b,c,d工況分別用實線,虛線,虛點線及點線標注,而離散模型則用顏色和符號標注.為方便表示,圖中橫坐標單位為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.

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

4 結 論

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

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

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

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

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

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

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