基于氣體動(dòng)理學(xué)格式的平板繞流數(shù)值模擬

謝 晴，謝軍龍，郭曉亮，莊曉東

華中科技大學(xué)能源與動(dòng)力工程學(xué)院，湖北 武漢 430074

描述流體運(yùn)動(dòng)有兩種方法，第一種是從宏觀角度描述流體，例如Euler方程，Navier-Stokes方程等；另一種是從微觀角度，它采用粒子速度分布函數(shù)描述高維相空間的粒子數(shù)密度，如Boltzmann方程?；诮橛^氣體動(dòng)理論的氣體動(dòng)理學(xué)格式（Gas Kinetic Scheme，GKS）從微觀角度為氣體流動(dòng)數(shù)值計(jì)算提供了很好的方法。Xu Kun等通過構(gòu)造不同的初始?xì)怏w分布函數(shù)以及平衡態(tài)分布函數(shù)，建立了等價(jià)于求解NS方程的GKS-NS格式以及理論上能夠求解全流域流動(dòng)的統(tǒng)一算法（UGKS格式）［1-3］。Jin等［4］在格式中引入速度空間自適應(yīng)技術(shù)，加快了計(jì)算速度，減少了多尺度流問題中的內(nèi)存需求。為高效模擬可壓縮流動(dòng)，學(xué)者研究了高精度氣體動(dòng)理學(xué)格式［5-7］。另外，在非結(jié)構(gòu)化網(wǎng)格［8］、動(dòng)網(wǎng)格［9］、浸入邊界［10］中應(yīng)用氣體動(dòng)理學(xué)格式進(jìn)行數(shù)值模擬，得到了精確的模擬結(jié)果。氣體動(dòng)理學(xué)格式在計(jì)算過程中表現(xiàn)出穩(wěn)定性，使得氣體動(dòng)理學(xué)格式得到越來越多的應(yīng)用。

繞流是流體力學(xué)研究的重要內(nèi)容之一［11-13］，在工業(yè)領(lǐng)域涉及面廣、影響面大，有重要的研究價(jià)值。目前國內(nèi)外應(yīng)用氣體動(dòng)理學(xué)格式在繞流方面的研究較少［14］。本文將介紹氣體動(dòng)理學(xué)格式在平板繞流中的應(yīng)用，選取合適的邊界條件構(gòu)造模擬二維平板繞流的物理模型。主要探究平板在不同間距條件下不同的繞流流場特性，揭示氣體動(dòng)理學(xué)格式在平板繞流數(shù)值模擬中良好的應(yīng)用前景。

1 氣體動(dòng)理學(xué)格式的基本原理

1.1 控制方程

在無外力場作用下，二維氣體動(dòng)理學(xué)BGK-Shakhov方程為：

式（1）中：f是氣體分布函數(shù)；u是粒子切向速度；v是粒子法向速度。

式（2）中：τ是碰撞時(shí)間；μ是動(dòng)力粘性系數(shù)；p是壓力。

修正后的分布函數(shù)為：

式（3）中：Pr是普朗特?cái)?shù)；c是隨機(jī)速度矢量；q是熱流矢量；R是通用氣體常數(shù)；T是溫度。

平衡態(tài)分布函數(shù)為：

式（4）中：ρ是密度；U，V是宏觀速度。

式（5）中：K是內(nèi)部自由度。

式（6）中：m是分子質(zhì)量；k是玻爾茲曼常數(shù)。

內(nèi)部自由度K和比熱比γ之間的關(guān)系如下所示：

宏觀守恒量與氣體分布函數(shù)的關(guān)系為：

由于氣體分子在碰撞過程中滿足質(zhì)量、動(dòng)量和能量守恒，因此有守恒性約束條件要求：

采用有限體積法對(duì)公式（1）進(jìn)行離散。物理空間，時(shí)間空間和粒子速度空間被劃分為許多小的 控 制 體 。 在 二 維 空 間 中 ，Ωi，j=ΔxΔy ，Δx=xi+1/2，j-xi-1/2，j，Δy=yi+1/2，j-yi-1/2，j。 粒 子 速度空間被劃分為矩形網(wǎng)格進(jìn)行離散，間距為Δu和Δv 。（k，l）速度間隔的中心位于（uk，vl）=(kΔu，lΔv)。則tn時(shí)刻的平均氣體分布函數(shù)為：

氣體動(dòng)理學(xué)的關(guān)鍵正是利用上述通解得到單元界面上隨時(shí)間演化的氣體分布函數(shù)，從而計(jì)算出數(shù)值通量，得到下一時(shí)刻的物理量［15］。

1.2 計(jì)算域和邊界處理

計(jì)算域?yàn)殚L方形，長L=400 mm，寬W=100 mm。平板布置方式為單平板（布置1）、前后平板（布置2）和上下平板（布置3）3種情況，具體布置方式及計(jì)算域如圖1所示。布置1、2中平板l=30 mm，布置3中平板長度l=20 mm。

圖1 平板布置形式及計(jì)算域示意圖：（a）布置1，（b）布置2，（c）布置3Fig.1 Schematic diagrams of flat-plate layout and calculation domain：（a）arrangement 1，（b）arrangement 2，（c）arrangement 3

求解流動(dòng)問題需要給出適當(dāng)?shù)倪吔鐥l件，邊界條件的處理方式對(duì)計(jì)算結(jié)果的精度和計(jì)算的穩(wěn)定性有很大的影響，各種邊界條件處理如下：

1）滑移壁面邊界：在邊界法線上根據(jù)鏡面反射原則取值，即把和邊界節(jié)點(diǎn)相鄰節(jié)點(diǎn)的各流動(dòng)量矢量值，取相反值賦到計(jì)算區(qū)域外的虛擬網(wǎng)格節(jié)點(diǎn)上。而在邊界切線方向上直接取為計(jì)算域內(nèi)與虛擬網(wǎng)格節(jié)點(diǎn)相鄰節(jié)點(diǎn)流動(dòng)量的值。

2）無反射邊界：根據(jù)零梯度條件，可通過外插的方法由流場的內(nèi)點(diǎn)得到宏觀邊界條件。零梯度條件為：

3）速度入口邊界：在射流入口處施加速度邊界 u=U0，v=0。

平板繞流中，入口邊界為速度入口，出口邊界為無反射邊界，上下邊界為滑移壁面邊界條件。網(wǎng)格數(shù)為401×101的結(jié)構(gòu)化網(wǎng)格。

1.3 算法準(zhǔn)確性對(duì)比驗(yàn)證

選用典型的Poiseuille流即流場無平板時(shí)驗(yàn)證程序的準(zhǔn)確性。Poiseuille流的解析解為：

式（16）中：μ為流體的動(dòng)力粘性系數(shù)；L為管道寬度；G=-dp/dx為管道內(nèi)的壓力梯度。

當(dāng)y=L/2時(shí)速度取得最大值：

以最大速度做無量綱處理，則：

計(jì)算域尺寸為 0≤x≤2，0≤y≤1。采用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格數(shù)為128×64。上下邊界采用滑移壁面邊界條件；左邊界采用壓力進(jìn)口邊界條件，入口壓力為1.01 MPa；右邊界采用壓力出口邊界條件，出口壓力為1.00 MPa。

在流域內(nèi)取3個(gè)截面，section 1，section2和section3分別表示與流向垂直的1/4位置截面、中間位置截面和3/4位置截面，截面示意圖如圖2所示。3個(gè)截面流向速度U與解析解的相對(duì)誤差δ如圖3所示，從圖3中可以看出：管道內(nèi)流體的速度與解析解最大誤差小于1.8%，且誤差較大的地方出現(xiàn)在上下邊界處，管道中心誤差很小，故氣體動(dòng)理學(xué)數(shù)值計(jì)算結(jié)果與解析解基本重合，初步驗(yàn)證了程序的正確性。

圖2 Poiseuille流截面示意圖Fig.2 Schematic diagram of Poiseuille flow section

圖3 Poiseuille流不同截面流向速度與解析解相對(duì)誤差圖Fig.3 Relative error diagram of flow velocity and analytical solution for different cross sections of Poiseuille flow

2 平板繞流數(shù)值模擬

利用氣體動(dòng)理學(xué)格式的基本理論，計(jì)算二維平板繞流的流體流動(dòng)情況。假設(shè)流體從左端均勻流入。

2.1 單個(gè)平板繞流

對(duì)不同雷諾數(shù)的單板繞流流場進(jìn)行模擬。圖4給出了布置1在不同雷諾數(shù)下的渦線圖。由圖4可以看出：Re≤100時(shí)，流場處于對(duì)稱尾流區(qū)，沒有渦的脫落，繞平板分離后的流體在平板后形成對(duì)稱的漩渦，并且只對(duì)平板附近流場有影響。但隨著雷諾數(shù)的增加，渦不斷被拉長，逐漸失去對(duì)稱性，渦開始在平板兩側(cè)周期性的脫落，形成一系列的渦即著名的卡門渦街。

圖4 布置1在不同雷諾數(shù)下的渦線圖：（a）Re=50，（b）Re=100，（c）Re=300Fig.4 Vortex diagrams of arrangement 1 at different Reynolds numbers：（a）Re=50，（b）Re=100，（c）Re=300

2.2 前后平板繞流

對(duì)雷諾數(shù)Re=300，兩板間距為10 mm，30 mm和50 mm三種情況下前后平板繞流流場進(jìn)行了模擬。圖5和圖6分別給出了布置2在不同板間距下的渦線圖和渦量圖。由圖5和圖6可知，對(duì)于低雷諾數(shù)下前后排列的平板繞流，隨著板間距的逐漸增大，其流場結(jié)構(gòu)呈現(xiàn)以下特征：1）當(dāng)平板間距較小時(shí)，兩個(gè)平板離的很近，此時(shí)流場形態(tài)類似于一個(gè)平板，流場處于對(duì)稱尾流區(qū)，沒有渦的脫落，如圖5（a）和圖6（a）所示；2）當(dāng)間距稍微增大時(shí)，上游平板脫落的剪切層附著于下游平板，在兩平板間出現(xiàn)渦量場，但仍沒有渦的脫落，如圖5（b）和圖6（b）所示；3）當(dāng)間距很大時(shí)，上游平板的渦脫落后擊打在下游平板上，并在下游出現(xiàn)卡門渦街，如圖5（c）和圖6（c）所示。當(dāng)平板距離比較小時(shí)沒有卡門渦街的形成，當(dāng)間距較大時(shí)則有卡門渦街的形成。

2.3 上下平板繞流

對(duì)雷諾數(shù)Re=300，兩板間距為10 mm，20 mm和40 mm三種情況，上下板繞流流場進(jìn)行了模擬。圖7和圖8分別給出了布置3在不同板間距下的渦線圖和渦量圖。

圖5 布置2在不同板間距時(shí)的渦線圖：（a）D=10 mm，（b）D=30 mm，（c）D=50 mmFig.5 Vortex diagrams of arrangement 2 at different plate spacings：（a）D=10 mm，（b）D=30 mm，（c）D=50 mm

圖6 布置2在不同板間距時(shí)的渦量圖：（a）D=10 mm，（b）D=30 mm，（c）D=50 mmFig.6 Vorticity diagrams of arrangement 2 at different plate spacings：（a）D=10 mm，（b）D=30 mm，（c）D=50 mm

圖7 布置3在不同板間距時(shí)的渦線圖：（a）D=10 mm，（b）D=20 mm，（c）D=40 mmFig.7 Vortex diagrams of arrangement 3 at different plate spacings：（a）D=10 mm，（b）D=20 mm，（c）D=40 mm

圖8 布置3在不同板間距時(shí)的渦量圖：（a）D=10 mm，（b）D=20 mm，（c）D=40 mmFig.8 Vorticity diagrams of arrangement 3 at different plate spacings：（a）D=10 mm，（b）D=20 mm，（c）D=40 mm

由圖7和圖8可知，對(duì)于低雷諾數(shù)下，上下排列的平板繞流，隨著板間距的逐漸增大，其流場結(jié)構(gòu)呈現(xiàn)以下特征：1）當(dāng)兩平板間距較小時(shí)，每個(gè)平板邊緣處有渦的脫落，并在距平板較近的地方渦合并形成共同的渦，出現(xiàn)卡門渦街，如圖7（a）和圖8（a）所示；2）隨著板間距的增大，在上下平板后面附近區(qū)域各自形成自己的渦，在離平板比較遠(yuǎn)的地方兩列渦開始逐漸合并形成一列渦，如圖7（b）和圖8（b）所示；3）當(dāng)板間距很大時(shí)，在兩平板下游各自出現(xiàn)卡門渦街，兩板之間的影響較小，如圖 7（c）和圖8（c）所示。

3 結(jié) 語

本文基于氣體動(dòng)理學(xué)格式對(duì)平板繞流進(jìn)行了數(shù)值模擬，主要結(jié)論如下：

1）對(duì)于單平板繞流，當(dāng)雷諾數(shù)較小時(shí)，流場處于對(duì)稱尾流區(qū)，沒有渦的脫落；當(dāng)雷諾數(shù)較大時(shí)，在平板下游會(huì)形成卡門渦街。

2）對(duì)于前后平板繞流，左右平板距離較小時(shí)沒有形成卡門渦街，但當(dāng)平板距離較大時(shí)就會(huì)形成卡門渦街。

3）對(duì)于上下平板繞流，當(dāng)平板間距比較小時(shí)，相當(dāng)于一個(gè)平板形成的繞流；而當(dāng)兩個(gè)平板距離比較大時(shí)，在平板下游較近距離處形成各自的渦旋，并且在離平板比較遠(yuǎn)的地方兩列渦開始逐漸合并形成一列渦。