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

華中科技大學能源與動力工程學院，湖北 武漢 430074

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

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

1 氣體動理學格式的基本原理

1.1 控制方程

在無外力場作用下，二維氣體動理學BGK-Shakhov方程為：

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

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

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

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

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

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

式（5）中：K是內部自由度。

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

內部自由度K和比熱比γ之間的關系如下所示：

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

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

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

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

1.2 計算域和邊界處理

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

圖1 平板布置形式及計算域示意圖：（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

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

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

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

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

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

1.3 算法準確性對比驗證

選用典型的Poiseuille流即流場無平板時驗證程序的準確性。Poiseuille流的解析解為：

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

當y=L/2時速度取得最大值：

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

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

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

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

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

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

利用氣體動理學格式的基本理論，計算二維平板繞流的流體流動情況。假設流體從左端均勻流入。

2.1 單個平板繞流

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

圖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 前后平板繞流

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

2.3 上下平板繞流

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

圖5 布置2在不同板間距時的渦線圖：（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在不同板間距時的渦量圖：（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在不同板間距時的渦線圖：（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在不同板間距時的渦量圖：（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可知，對于低雷諾數(shù)下，上下排列的平板繞流，隨著板間距的逐漸增大，其流場結構呈現(xiàn)以下特征：1）當兩平板間距較小時，每個平板邊緣處有渦的脫落，并在距平板較近的地方渦合并形成共同的渦，出現(xiàn)卡門渦街，如圖7（a）和圖8（a）所示；2）隨著板間距的增大，在上下平板后面附近區(qū)域各自形成自己的渦，在離平板比較遠的地方兩列渦開始逐漸合并形成一列渦，如圖7（b）和圖8（b）所示；3）當板間距很大時，在兩平板下游各自出現(xiàn)卡門渦街，兩板之間的影響較小，如圖 7（c）和圖8（c）所示。

3 結 語

本文基于氣體動理學格式對平板繞流進行了數(shù)值模擬，主要結論如下：

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

2）對于前后平板繞流，左右平板距離較小時沒有形成卡門渦街，但當平板距離較大時就會形成卡門渦街。

3）對于上下平板繞流，當平板間距比較小時，相當于一個平板形成的繞流；而當兩個平板距離比較大時，在平板下游較近距離處形成各自的渦旋，并且在離平板比較遠的地方兩列渦開始逐漸合并形成一列渦。