尺度自適應的離散統(tǒng)一氣體動理學格式

許 丁 劉 欣 孫 祥

(1. 西安交通大學航天航空學院 機械結構強度與振動國家重點實驗室，西安 710049;2. 陜西省先進飛行器服役環(huán)境與控制重點實驗室，西安 710049)

1 引言

隨著人們對多尺度流動問題的深入研究，傳統(tǒng)基于連續(xù)介質假設求解Navier-Stokes 方程的方法顯示出更多的局限性，而基于氣體動理學理論的方法所具有的優(yōu)勢則更多的被人們重視.氣體動理學從微觀角度出發(fā)，引入氣體速度分布函數的概念并建立其滿足的Boltzmann 方程.目前基于Boltzmann 方程的氣體動理學格式應運而生，如格子Boltzmann 方法(Lattice Boltzmann Method, LBM)[1-3]、氣體動理學格式(Gas Kinetic Scheme, GKS)[4,5]、統(tǒng)一氣體動理學格式(Unified Gas Kinetic Scheme, UGKS)[6,7]、離散統(tǒng)一氣體動理學格式(Discrete Unified Gas Kinetic Scheme,DUGKS)[8,9]等.這些方法與基于Navier-Stokes 方程的求解方法相比，可以獲得更多的流動信息，目前已取得眾多成功應用[10-16].但是，這些動理學方法又各有一些局限與不足.例如，LBM 作為一種有限差分方法對低速不可壓縮流動問題有著很好的適用性，但其對均勻網格的需求限制了其在高Re數下邊界層流動問題中的應用.GKS 作為一種有限體積方法，定位于求解傳統(tǒng)Navier-Stokes 方程，不能應用于連續(xù)介質假設失效的較大Kn 數流動.UGKS 方法是一種適用于多尺度流動問題的統(tǒng)一方法，在處理多尺度復雜流動中取得了較好的結果，但該方法要同時追蹤宏觀守恒量和分布函數，并要計算平衡態(tài)分布函數在界面處的空間導數和時間導數從而涉及大量矩陣運算，計算成本較高.DUGKS 方法沿特征線離散Boltzmann 方程，可認為是UGKS 方法的一種簡化形式，相比UGKS 來說DUGKS 計算中只需追蹤分布函數，且不再需要計算平衡態(tài)分布函數的時空導數，降低了計算開銷.但是DUGKS 在處理碰撞項時采用碰撞項的顯式處理和隱式處理的算術平均，這種處理方法主要基于數值的角度，未曾從物理上深入考慮兩個時間尺度即計算時間步長和碰撞松弛時間之間的相對大小對粒子碰撞過程宏觀表現(xiàn)的影響.

本文在原始DUGKS 基礎上，將碰撞項的處理改成顯式和隱式的加權平均，其中權系數依賴于當地碰撞松弛時間和計算時間步長之比，即當地Kn 數.通過數學和物理上的分析發(fā)現(xiàn)，該權系數可以根據當地Kn 數做出自適應調節(jié)，從而使所得方法具有尺度自適應特性，適用于所有Kn 數下的流動.

2 Boltzmann-BGK 模型方程沿特征線離散的一般形式

碰撞項采用BGK 模型的Boltzmann 方程如下

其中f=f(xi,t,uj)為t時刻、位置xi處、速度為uj的粒子速度分布函數，τ為碰撞松弛時間，g為平衡態(tài)分布函數

其中ρ為流體密度，R為氣體常數，T為溫度，D為空間維度，Ui為宏觀速度.宏觀量可以通過分布函數取矩得到，如對于二維等溫流動問題

其中dΞ = dudv, ψ= (1,u,v)T為碰撞不變量.根據粒子碰撞過程滿足質量與動量守恒定律，碰撞項需在任意時刻、任意空間位置滿足如下相容關系

對于Boltzmann-BGK 方程，若松弛時間τ為常數，則t+s時刻的分布函數可形式上寫為[4]

其中t為初始時刻，s為時間發(fā)展長度，x′i=xi-ui(s-t′)為粒子運動軌跡，fini(xi-uj) =f(xi,t,uj)代表初始時刻t的非平衡態(tài)分布函數.盡管Boltzmann-BGK 方程(1)基于微觀觀點建立，其所蘊含的流動尺度卻不限于微觀動理學尺度，這一點可以從式(5)清楚看出：式(5)右邊第一項積分項反映粒子之間相互大量碰撞作用的積累效果，該項已被證明可以看作為宏觀連續(xù)流動動力學模型[17]，對應為宏觀動力學尺度；式(5)右邊第二項反映粒子的遷移運動，對應微觀動理學尺度.式(5)將宏觀動力學尺度和微觀動理學尺度有機地結合到一起，可以很好反映多尺度流動特性.但是式(5)并不能直接使用，因為右邊第一項涉及積分的平衡態(tài)g是未知的，需要滿足相容關系式(4).因此從這個角度上來說式(5)只是式(1)形式上的解.

另一方面，認識到Boltzmann-BGK 方程(1)自身已將宏觀動力學尺度和微觀動理學尺度有機地結合到一起，因此可以直接從Boltzmann-BGK 方程(1)出發(fā)來構造多尺度計算格式.這里將(1)式沿特征線進行離散，可得到

其中gini(xi,uj) =g(xi,t,uj)和Ωini(xi,uj) = Ω(xi,t,uj)分別代表初始時刻t的平衡態(tài)分布函數和碰撞項.上式(6)右端的碰撞項采用了顯式和隱式加權平均處理.α為隱式部分對應的權系數，其為碰撞松弛時間τ和時間發(fā)展長度s的比值τ/s的函數，具體形式待定.上式(6)可整理為

式(7)為Boltzmann-BGK 方程沿特征線離散的一般形式，其右端前兩項和平衡態(tài)分布函數相關反映粒子之間的碰撞作用，和式(5)右端第一項積分項對應；式(7)右端第三項和初始非平衡態(tài)分布函數相關反映粒子的遷移效應并和式(5)右端第二項相對應.通過對比式(5)和式(7)中的非平衡態(tài)分布函數相關項，并令?3=e-s/τ，可確定權系數α為

由上式(9)可以清楚看出權系數α數學上直接依賴于時間比值τ/s，物理上的討論將在后面第4 節(jié)給出.確定了α后將其帶入到(8)式，可得系數?1, ?2和?3具體形式

為了進一步分析式(7)的數學物理特性，這里將從兩種極限流動的角度來討論，首先將其整理為以下形式

對于Kn→0 的宏觀連續(xù)流動問題，物理上理解為在一個時間步長內粒子經歷的碰撞次數足夠多，即τ ?s，此時數學上給出β1→0, β2→τ/s，從而(11)式可化為

其中Dt=?t+uj?xj.眾所周知，當對Boltzmann-BGK 方程(1)進行Chapman-Enskog展開至O(τ)階時也得到(13)式[4]，其對應宏觀方程為Navier-Stokes 方程.由此說明在宏觀連續(xù)流動區(qū)域式(7)對應為Navier-Stokes 方程，此時流動尺度為宏觀動力學長度和時間尺度.

另一方面，對于Kn?0 的稀薄流動或微尺度流動問題，有τ ?s，此時數學上給出β1→1, β2→1，從而(11)式化為

上式是自由分子流運動的解，此時流動對應微觀動理學尺度，長度尺度為分子平均自由程尺度，時間尺度為碰撞松弛時間尺度.

由上述分析可知，Boltzmann-BGK 方程沿特征線離散的一般形式(7)將宏觀動力學尺度和微觀動理學尺度有機地結合到一起，其中引入的系數?1, ?2和?3依賴于比值τ/s且可根據當地流動狀態(tài)的不同進行自適應調節(jié)，兩個時間尺度即碰撞松弛時間尺度和時間步長的比值τ/s在調節(jié)中起到了關鍵作用.

3 尺度自適應的離散統(tǒng)一氣體動理學格式

本文在文獻現(xiàn)有離散統(tǒng)一氣體動理學格式DUGKS[8]基礎上發(fā)展出具有尺度自適應特性的離散統(tǒng)一氣體動理學格式(Scale Adaptive Discrete Unified Gas Kinetic Scheme,SADUGKS).二者主要差別體現(xiàn)在碰撞項的處理不同，在DUGKS 中碰撞項采用算術平均，權系數恒為常數αDUGKS≡1/2；而在SADUGKS 中，碰撞項采用加權平均處理，權系數由式(9)給出，可根據當地流動狀態(tài)的不同進行自適應調節(jié).

將Boltzmann-BGK 方程(1)在網格單元Vj及時間步長(tn,tn+Δt)內積分，得

對式(15)中通量的時間積分項采用中點公式處理

對式(15)中碰撞項的時間積分采用加權梯形公式

其中α(Δt)為式(9)給出的權系數

將式(17)和式(18)代入到式(15)中，可得

從式(20)可以看出該式右端包含了碰撞項和通量的隱式部分，因此，在使用該式進行迭代遞推之前必須先解決碰撞項和通量的隱式離散問題，下面對其分別進行處理.

其中t=tn為初始時刻，時間發(fā)展長度為r= Δt/2, xi,b為第m個網格單元界面位置，α(r)由式(9)給出.引入輔助分布函數ˉf定義

上式也等價于

接下來，處理碰撞項中的隱式部分Ωn+1，引入新的輔助分布函數～f定義如下

由相容關系式(4)同時結合定義式(30)可知

將式(30)代入到式(20)，整理得

再次指出的是在原始DUGKS 格式中，從式(20)至(26)，式(30)和(32)中出現(xiàn)的所有權系數α恒取為αDUGKS≡1/2，而SADUGKS 中權系數形式由式(9)給出.

需要說明的是上述在推導SADUGKS 時速度空間是連續(xù)的，實際在使用SADUGKS 時還需將速度空間離散，其離散方法與DUGKS 格式[8]相同.

關于邊界條件的處理分兩類，針對宏觀連續(xù)流動問題，SADUGKS 格式在固體壁面處采用了無滑移邊界條件；而對大Kn 數流動問題，則采用了完全漫反射邊界條件，具體實現(xiàn)與DUGKS 格式[8]相同，亦不再贅述.

4 尺度自適應的離散統(tǒng)一氣體動理學格式的性質討論

首先，SADUGKS 格式基于DUGKS 格式發(fā)展出來，SADUGKS 格式同樣具有漸近保持性[8,18]，這一點可以清楚地從式(13)看出.

其次，SADUGKS 格式具有尺度自適應特性，能夠根據當地流動狀態(tài)不同進行自適應調節(jié).SADUGKS 格式中，時間步長選取按照CFL 條件確定

其中Δ 為空間計算網格尺度，CFL 為CFL 數；|u|max為粒子最大離散速度，通常與聲速c同量級.

同時，碰撞松弛時間τ可近似為

式中l(wèi)s為分子平均自由程，ˉu為粒子平均運動速率與聲速c同量級.因此可以得到

其中Kn 為當地努森數，反映當地的流動狀態(tài)和流動尺度.在SADUGKS 格式中，用來計算界面通量的式(21)和碰撞項加權離散的式(18)中皆出現(xiàn)權系數α，而權系數α依賴于比值τ/Δt，即依賴于當地Kn 數，因此，α取值的大小會隨著當地流動狀態(tài)和尺度的不同而改變.

眾所周知Boltzmann-BGK 模型將氣體從當前非平衡狀態(tài)向平衡態(tài)的過渡看作為一個馳豫過程.對于Kn→0 的宏觀連續(xù)流動問題，從物理上分析，此時碰撞弛豫時間遠小于時間步長，τ ?Δt，即在一個時間步長Δt內粒子經歷足夠多次的碰撞，因此弛豫過程受初始時刻碰撞項Ωini(xi-uiΔt,uj)的影響很小，即粒子的碰撞歷史效應可忽略，當前的分布函數f(xi,t+Δt,uj)將只決定于當前的碰撞Ω(xi,t+Δt,uj)；而另一方面，從式(6)的數學表達上分析，此時α →1，式(6)簡化為

可見數學結果與上述物理分析一致.

對于Kn?0 的稀薄流動或微尺度流動問題，從物理上分析，碰撞弛豫時間遠大于時間步長，τ ?Δt，即在一個時間步長Δt內粒子幾乎不發(fā)生碰撞，從而每次發(fā)生的碰撞都對當前的分布函數f(xi,t+ Δt,uj)有著至關重要的影響，因此初始碰撞作用Ωini(xi-uiΔt,uj)和當前碰撞作用Ω(xi,t+Δt,uj)都要保留.另一方面，從數學上來看，此時α →1/2，式(6)可簡化為

可見此時數學結果也與上述物理分析相一致.

根據上述分析可知，式(9)給出的權系數α緊密依賴于時間比值τ/Δt，并可根據當地流動狀態(tài)和流動尺度的不同進行自適應調節(jié)，因此SADUGKS 是一種尺度自適應的格式，適用于從連續(xù)流動到滑移流動、過渡流動及自由分子流所有流動問題，在第5 節(jié)中該格式的應用也很好地體現(xiàn)出了這一點.

5 在低速不可壓縮流動問題中的應用

5.1 二維非定常Taylor 渦問題

為了驗證SADUGKS 格式的空間精度，本文對二維非定常Taylor 渦問題[8]進行模擬.該問題的解析解如下所示其中，平衡態(tài)分布函數g可用式(39)由初始解析解得到；邊界條件采用周期性邊界條件.離散速度及權系數由3 點Gauss-Hermite 公式[8]決定.

為分析格式的空間精度，本文采用了一系列不同數量(N×N)的均勻網格；為了消除時間步長對誤差的影響，在計算中將時間步長設置為固定小量(Δt= 2τ).統(tǒng)計在tc=ln(2)/(8νπ2)時刻流場速度的誤差L2, L2定義為

其中U 和Ue分別指數值解速度和解析解速度.

采用64×64 均勻網格得到的不同時刻的速度剖面，如圖1 所示，可以看到SADUGKS格式得到的數值解與解析解能夠很好地吻合，證明了SADUGKS 格式對于Taylor 渦問題的適用性；表1 展示了在tc時刻網格量與誤差的關系，從中可以看出SADUGKS 格式具有空間2 階精度.

表1 SADUGKS 格式模擬Taylor 渦問題時的誤差與空間精度

圖1 SADUGKS 格式模擬Taylor 渦問題不同時刻的速度剖面(64×64 均勻網格)

5.2 宏觀連續(xù)流動下方腔驅動流問題

二維方腔驅動流問題為計算流體力學中驗證數值算法的經典算例，隨著雷諾數變化，流動結構及旋渦等也隨之改變，因此能夠很好地驗證算法的粘性分辨率及旋渦捕捉能力.為了驗證SADUGKS 格式，本節(jié)對宏觀連續(xù)流動下的二維方腔驅動流問題進行了數值模擬，并與文獻[19]的結果以及LBM 方法的計算結果進行了對比.

方腔邊長為L，上邊界以U的速度勻速運動，其余邊界靜止，流動雷諾數為Re=UL/ν，其中ν為運動學粘性系數.本節(jié)在計算中，令L= 1.0, U= 0.1, RT= 1/3，故流動馬赫數較小，流動可作為不可壓縮流動.CFL 設置為0.5，離散速度及權系數由3 點Gauss-Hermite 公式[8]決定，以保證與LBM 采用的D2Q9 格子模型一致.LBM 方法中所有速度邊界均采用非平衡態(tài)外推邊界條件[20]，SADUGKS 格式中所有邊界均采用無滑移邊界條件[8].

雷諾數Re= 1000 時，不同計算網格下沿方腔中線的速度剖面，如圖2 所示.從圖2 左圖可知，在采用相同的50× 50 均勻網格的條件下，SADUGKS 格式得到了比LBM 方法更好的結果，特別是下壁面附近優(yōu)勢更為明顯，表明SADUGKS 格式對壁面附近粘性邊界層的分辨較好.此外，作為一種有限體積格式，SADUGKS 格式可以使用非均勻網格，通過在壁面處進行加密來進一步減小計算誤差及所需網格量.圖2 右圖展示了SADUGKS 格式采用在壁面處進行局部加密的30×30 非均勻網格、LBM 方法采用100×100 均勻網格的結果對比，可以看出兩種方法得到的結果都與文獻[19]的結果吻合較好.

圖2 方腔驅動流Re=1000 時SADUGKS 與LBM 結果對比

此外，本文采用SADUGKS 格式計算了更高雷諾數的方腔驅動流，以檢驗SADUGKS格式對高雷諾數流動中旋渦的捕捉能力及格式的穩(wěn)定性.采用30×30 均勻網格、與圖2 左圖中相同的30×30 非均勻網格及加密后的50×50 非均勻網格等三種網格，雷諾數選擇為Re= 3200,5000,7500 和10000，計算結果如圖3 所示.可以看出SADUGKS 格式的結果總體上與文獻[19]的結果吻合較好，且隨著網格不斷加密，速度峰值的捕捉更加精確，50×50 非均勻網格足以很好地反映高雷諾數下的流動狀態(tài).在穩(wěn)定性方面，可以看到即使采用較為稀疏的30×30 的均勻網格，SADUGKS 格式也能夠對上述所有高Re數流動得到物理上合理的結果，說明該格式具有較好的穩(wěn)定性.

圖3 不同網格不同Re 數下，SADUGKS 計算方腔驅動流所得速度剖面

最后，分別使用DUGKS 和SADUGKS 格式計算了雷諾數為Re=3200,5000,7500和10000 的方腔驅動流，為了消除非均勻網格帶來的影響，這里采用50×50 的均勻網格，兩種格式計算時邊界條件設置和CFL 數取法相同，兩種格式的計算結果見圖4.可以看出從整體上來說SADUGKS 和DUGKS 所得結果均與文獻[19]的結果吻合較好，但是在局部點如速度極值點附近，SADUGKS 與文獻[19]的結果更為貼近，并且隨著Re的增加，DUGKS 在極值點附近和文獻[19]的結果偏差增大.

圖4 不同Re 數下DUGKS 和SADUGKS 計算方腔驅動流所得速度剖面

5.3 宏觀連續(xù)流動下后臺階流動問題

另一個經典宏觀連續(xù)流動問題是后臺階流動問題，流動過程中包含有流動分離與再附.本節(jié)對后臺階流動進行了數值模擬，并與文獻[21]中的實驗結果進行對比.流動示意圖見圖5 所示，坐標原點取在臺階處.

圖5 后臺階流動區(qū)域示意圖

來流入口速度為拋物線分布

其中Umax為最大速度.在后臺階流動中，雷諾數定義為Re= (H-h)Umax/ν.本文計算中取Umax= 0.1, l= 4.0, H= 3.0, h= 2.0, L足夠長(L= 35)以保證出口處流動充分發(fā)展，可采用非平衡態(tài)外推邊界條件[20]；臺階及上壁面處均采用無滑移邊界條件[8].本文采用均勻網格(250×100)，數值模擬了Re= 50 及Re= 150 的后臺階流動.圖6 和圖7 分別給出了Re=50 及Re=150 在四個不同位置處的水平速度剖面結果.

從圖6 中可以看出，Re= 50 的流動在x= 1.5 處出現(xiàn)了回流，說明流動經過后臺階發(fā)生流動分離，產生了旋渦；在臺階下游處，隨著x逐漸增大，速度逐漸發(fā)展為拋物線分布，表明臺階的影響逐漸減弱.圖7 顯示了相似的變化規(guī)律.此外，Re= 150 時，在x= 4.0 處仍然存在回流，表明隨著Re增加，粘性降低，旋渦的尺寸逐漸增大.在上述兩種不同Re下，SADUGKS 格式得到的數值結果均能與實驗結果較好地吻合，再次說明SADUGKS 格式對粘性流動分離、再附及旋渦具有較好的捕捉能力.

圖6 后臺階流動不同位置處的速度剖面圖，SADUGKS 所得結果與實驗結果的對比

圖7 后臺階流動不同位置處的速度剖面圖，SADUGKS 所得結果與實驗結果的對比

5.4 不同Kn 數下的等溫Couette 流動問題

本節(jié)中將SADUGKS 格式用于求解不同Kn 數下等溫Couette 流動，以驗證格式對大Kn 數問題的適用性.流動描述如下：兩無限大平板互相平行，距離H，之間充滿靜止流體，上下板分別以Uw和-Uw的速度沿水平方向反向運動，由于壁面處的粘性剪切作用，板間流體會運動并最終達到穩(wěn)定狀態(tài).這一問題已被諸多學者用包括DSMC[22]、UGKS[7]及DUGKS[8]在內的多種方法研究過，可以很好地檢驗格式求解大Kn 數問題的正確性及粘性剪切力的計算.

本文中，物理空間在豎直方向采用40 個均勻網格進行離散，速度空間采用28×28 半平面Gauss-Hermite[8,23]積分離散.入口與出口采用周期性邊界條件，上下壁面采用完全漫反射邊界條件[8].研究了多個不同Kn 數下的流動，并與文獻[22]結果及文獻[24]利用線性化Boltzmann 方程LBE 計算所得的結果進行對比.

三種方法所得不同Kn 數下的速度剖面，如圖8(a)所示，由于流動的對稱性，圖中只給出了流場的上半部分.從圖中可以看出，在所有Kn 數下，SADUGKS 格式的結果與其他方法均能很好地吻合，且隨著Kn 數不斷增大，壁面處的速度滑移效應也越來越顯著，三種方法均能很好地反映出大Kn 數流動下的速度滑移效應.圖8(b)所示為壁面處的切應力τw隨Kn 數的變化曲線，考慮到較大Kn 數下連續(xù)介質假設不成立及牛頓粘性切應力公式失效，圖中壁面切應力采用氣體動理學理論計算，即

圖8 SADUGKS 格式計算不同Kn 數(k =Kn/2)下的Couette 流動

6 結論

本文首先提出了Boltzmann 方程沿特征線離散的一般形式，在該式中碰撞項的離散采用了顯式和隱式加權平均的方法，其中權系數α依賴于當地碰撞松弛時間和計算時間步長的比值.同時基于該式，在文獻現(xiàn)有離散統(tǒng)一氣體動理學格式(DUGKS)基礎上提出了具有尺度自適應的離散統(tǒng)一氣體動理學格式(SADUGKS).接下來，本文對SADUGKS 格式的尺度自適應特性進行了討論，發(fā)現(xiàn)在SADUGKS 格式中，權系數α與反映當地流動尺度的當地努森數Kn 直接相關，當流動尺度發(fā)生變化時，α隨之變化，使得SADUGKS 格式能夠進行自適應調節(jié)，因此可以處理從宏觀連續(xù)流動到自由分子流的多尺度流動問題.最后，本文將SADUGKS 格式應用到若干低速不可壓縮流動問題中，通過這些經典算例對該格式進行了檢驗.計算結果表明：基于有限體積法的SADUGKS 在網格選用上更加靈活，格式穩(wěn)定性較好，對高Re數下的流動邊界層具有較高的分辨能力；對于流動分離、再附和旋渦也有較好的捕捉能力；此外，能夠很好地描述從小Kn 數的宏觀連續(xù)流動一直到大Kn 數的自由分子流動.

目前使用SADUGKS 僅在若干經典低速不可壓縮流動算例中進行了驗證，后續(xù)工作將在更加復雜的多尺度流動中進一步進行檢驗.另一方面，本文的研究基于Boltzmann-BGK 模型方程，結果僅適用于低速不可壓縮流動，后續(xù)將嘗試將SADUGKS 推廣到Boltzmann-Shakhov 模型方程中，并應用于非等溫、高速可壓縮流動及伴有間斷現(xiàn)象的復雜流動中去.