基于氣體動理論的二維Karman渦街?jǐn)?shù)值模擬

吳俊林, 李志輝, 蔣新宇, 彭傲平

(1. 中國空氣動力研究與發(fā)展中心超高速空氣動力研究所, 四川綿陽 621000; 2. 國家計算流體力學(xué)實驗室, 北京 100191)

引 言

Boltzmann方程[1]是描述氣體分子運動狀態(tài)和運動規(guī)律的確定論基本方程, 主要描述單原子氣體或簡單氣體的分子輸運特性. 不考慮外力場影響的Boltzmann方程如公式(1)所示

(1)

其中,f為氣體分子的速度分布函數(shù)；t，r，ξ分別為時間、 位置空間坐標(biāo)、 分子速度；ξr為碰撞分子的相對運動速度；σ為微分碰撞截面, 與分子之間的作用力模型相關(guān)；Ω為兩個分子碰撞后相對速度的方向；R3表示ξ1的三維速度空間積分.

Boltzmann方程通過對氣體分子速度分布函數(shù)矩積分得到流動的宏觀參數(shù), 可描述從稀薄流到連續(xù)流跨流域氣體流動由非平衡態(tài)向平衡態(tài)演化的過程[2]. 然而, 直接采用理論分析或數(shù)值方法研究Boltzmann方程很難實現(xiàn), 困難在于[3-4]: (1) 幾率密度分布函數(shù)f具有至少7個自變量, 維度太高; (2) 碰撞項高度非線性, 且具有高維積分特點, 并與分子碰撞模型相關(guān), 極其復(fù)雜. 雖然Aqarwal等[5]、 Kolobov等[6]在求解Boltzmann方程或廣義Boltzmann方程(generalized Boltzmann equation， GBE)方面取得了一些進(jìn)展, 但也并未實現(xiàn)對碰撞項的直接求解, 并且由于求解過程復(fù)雜, 計算代價很大, 導(dǎo)致很難應(yīng)用到工程中.

正是由于Boltzmann方程具有高維度、復(fù)雜碰撞積分項的特點, 對于跨流域非定常問題, 特別是稀薄過渡區(qū)的低速非定常流動, 其數(shù)值求解應(yīng)該考慮基于Boltzmann的運動模型方程[7-12]. 近年, Polikarpov等[8]通過求解非定常Shakhov模型方程得到了氣流穿過方形切口的瞬態(tài)流動, 并說明定常流場建立的時間與背壓壓比和氣體稀薄度相關(guān). Lihna-ropoulos等[9]通過數(shù)值求解時間依賴的BGK運動方程分析了圓柱管道中的氣體啟動過程. Chigullapalli等[7]求解Boltzmann-ESBGK 模型方程發(fā)展了一種三維非定常稀薄流動求解器. 實踐表明, 基于Boltzmann模型方程的數(shù)值求解方法是模擬稀薄流到連續(xù)流跨流域非定常流動的有效手段.

為了求解Boltzmann模型方程以得到跨流域數(shù)值模擬方法, 李志輝等在速度空間采用離散速度坐標(biāo)法, 位置空間應(yīng)用計算流體力學(xué)有限差分方法, 建立起從稀薄流到連續(xù)流的跨流域氣體動理論統(tǒng)一算法(gas-kinetic unified algarithm, GKUA)[13-14], 并在定常流動研究中得到了較好的應(yīng)用[15-16]. 實際上, 統(tǒng)一算法的求解過程是通過長時間非定常輸運模擬得到最終的定常流動狀態(tài), 因此可通過技術(shù)改進(jìn)得到整個非定常氣體流動演化過程. 前期已分析了一維、二維跨流域非定常流動特征與稀薄效應(yīng)對非定常流場影響機(jī)制[17], 并用于求解分析稀薄平面噴流擴(kuò)散進(jìn)入真空環(huán)境的非定常過程[18]. 結(jié)果表明, 基于GKUA的非定常流動求解器對于高真空自由分子流、稀薄流、過渡流等流區(qū)的超聲速非定常流動具有較好適用性. 因此, 本文對連續(xù)流區(qū)的經(jīng)典Karman渦街問題進(jìn)行模擬, 以確認(rèn)該跨流域非定常流動模擬算法對于連續(xù)流區(qū)低速流動的適應(yīng)性.

1 控制方程及計算方法

1.1 控制方程

GKUA能夠求解各種模型方程[19], 包括BGK[20], ES-BGK[21]， Shakhov[22]或者其他模型. 對于氮氣或者空氣, 室溫狀態(tài)下分子的轉(zhuǎn)動能就已完全激發(fā)[23], 因此采用考慮轉(zhuǎn)動自由度影響的雙原子氣體動理論模型方程能夠更好地描述氮氣流動. 此外, 眾多學(xué)者已經(jīng)對高溫條件下多原子氣體分子的振動能開展了有意義的分子動力學(xué)建模分析和數(shù)值計算工作[24-27], 后續(xù)如果研究高速、高溫氣體動力學(xué)問題時應(yīng)該加以考慮.

基于對轉(zhuǎn)動自由度松弛變化特性的Rykov模型研究[28-30], 采用轉(zhuǎn)動慣量描述氣體分子的自旋運動, 將分子總角動量守恒作為一個碰撞不變量, 在求解Boltzmann模型方程的統(tǒng)一算法理論框架[13-16]下, 構(gòu)建考慮轉(zhuǎn)動非平衡效應(yīng)的Boltzmann模型方程

(2)

式中, 下標(biāo)i表示空間維度方向; 角標(biāo)t， r分別表示彈性碰撞和非彈性碰撞, 如νt和νr是彈性碰撞頻率和非彈性碰撞頻率;μ為黏性系數(shù);Z為非彈性碰撞松弛因子；xi為位置空間坐標(biāo)；ξi為氣體分子速度分量；c為氣體分子熱運動速度；n，T，P分別為氣體分子數(shù)密度、氣體流動溫度和壓力；k為Boltzmann常數(shù)；m為分子質(zhì)量；δ對于分子間相互作用規(guī)律來說是一個常數(shù)δ=μt/(mnD),D是自擴(kuò)散系數(shù).

物理空間宏觀流動參數(shù)可以由f0和f1在速度空間上積分[13,28]求得

kTr=ε,Pt=nkTt,P=nkT

黏性系數(shù)μ和非彈性碰撞松弛因子Z可定義為如下的形式[29-30]

μt=μ(Tt)=μ(T*)t2/3/φ(t),t=Tt/T*

φ(t)=0.767+0.233t-1/6exp[-1.17(t-1)]

對于氮氣N2而言,T*=91.5 K,δ=1/1.55.

1.2 無量綱化和速度空間的離散

采用氣體動理論數(shù)值計算方法可直接捕捉速度分布函數(shù)隨時間的演化. 空間坐標(biāo)、時間、數(shù)密度、流動速度、溫度、黏性的無量綱參考量為[13-16]

而能量(包括平動能和轉(zhuǎn)動能)的無量綱參考量以及分布函數(shù)f0和f1的無量綱參考量分別為

mn∞(2RT∞)3/2,n∞(2RT∞)-3/2

mn∞RT∞(2RT∞)-3/2

對于二維流動, 為了對速度空間進(jìn)行降維處理, 在對速度分布函數(shù)方程(2)應(yīng)用離散速度坐標(biāo)法之前, 需要引入3個約化速度分布函數(shù), 對模型方程和基于速度空間的宏觀流動矩積分進(jìn)行約化處理, 這樣可以減少對計算機(jī)內(nèi)存的需求, 提高計算效率.

得到關(guān)于二維約化速度分布函數(shù)gi(t,x,y,ξx,ξy),i=1,2,3的方程. 在此基礎(chǔ)上, 使用離散速度坐標(biāo)法對速度分量ξx，ξy進(jìn)行數(shù)值離散.

應(yīng)用離散速度坐標(biāo)方法將物理空間的坐標(biāo)(x,y)轉(zhuǎn)換到計算平面(ζ,η), 模型方程轉(zhuǎn)化成為每個離散速度坐標(biāo)點上包含非線性源項的雙曲型守恒方程

(3)

在離散速度坐標(biāo)法中, 采用合適的積分規(guī)則對分布函數(shù)在速度空間進(jìn)行矩積分, 其中積分所依賴的分布函數(shù)僅在某些離散速度坐標(biāo)點上有確定值. 對分布函數(shù)在全速度域上的積分可基于這些離散速度點上的分布函數(shù)值得到. 求解方式是采用Gauss-Hermite數(shù)值積分法, 將全域積分轉(zhuǎn)變?yōu)榛跈?quán)重因子的求和. 基于此, 構(gòu)造某些自適應(yīng)選取的離散速度坐標(biāo)點(ξxσ,ξyε), 以及對速度空間這些點上的分布函數(shù)值進(jìn)行合理迭代求解至關(guān)重要, 這將直接影響數(shù)值模擬算法的收斂性和適用性. 值得注意的是, 這些離散速度坐標(biāo)點確定的速度分布函數(shù)值必須保持正定性.

1.3 適于非定常流動計算的差分方法

對于非定常流場的數(shù)值計算, 除了差分格式須滿足耗散控制、色散控制和激波控制條件以保證計算過程中的穩(wěn)定性、不產(chǎn)生虛假波動以及良好地捕捉激波外, 還應(yīng)特別滿足保頻譜原則, 并且要求計算網(wǎng)格和邊界計算方法要與內(nèi)點計算匹配協(xié)調(diào)[31-32]. 一般采用固定時間步長的方法進(jìn)行顯式的時間推進(jìn), 但時間步長必須滿足流場計算的格式穩(wěn)定性條件.

基于非定常時間分裂方法[4,13-16,19]將控制方程(3)分裂成為包含非線性源項的碰撞松弛方程和計算平面(ζ,η)上的對流運動方程. 采用3階WENO格式[19,33]離散速度分布函數(shù)方程的對流項. 源項的碰撞松弛積分采用3階Runge-Kutta 方法求解. 考慮到在真實氣體流動中對流運動和碰撞松弛具有相似的過程, 在數(shù)值計算中需要前向時間步長和后向時間步長的耦合迭代. 基于此, 在時間推進(jìn)上采用基于非定常時間分裂法的3階顯式Runge-Kutta[34]時間推進(jìn), 得到在各個離散速度坐標(biāo)點處數(shù)值求解約化速度分布函數(shù)方程的氣體運動論耦合迭代數(shù)值格式

(4)

其中, 下標(biāo)s表示源項,ζ,η為計算平面坐標(biāo). (4)式表示在Δt時間步長之后的分布函數(shù)值通過半時間步長(Δt/2)的兩次松弛演化得到, 這樣可以提高該差分離散格式的時間精度.

方程(4)的每一項對應(yīng)于非線性方程(3)的數(shù)值計算, 其中Ls,Lη和Lζ分別表示碰撞松弛源項以及兩個計算平面坐標(biāo)方向的對流運動項的差分方程.

引入3階WENO格式[35]對計算平面進(jìn)行掃描. 對于時間項, 可以使用Runge-Kutta 法(RK-3)來提高精度.

數(shù)值計算過程中的時間步長Δt由格式穩(wěn)定條件約束

其中, CFL為時間步長調(diào)節(jié)系數(shù), 一般CFL<1.

2 Karman渦街形成及運動的數(shù)值模擬

2.1 問題描述

一定條件下的定常來流繞過鈍體時, 物體兩側(cè)可能會周期性地脫離出旋轉(zhuǎn)方向相反、排列規(guī)則的雙列線渦, 經(jīng)過非線性作用后會形成Karman渦街. 以圓柱繞流為例, 低Reynolds數(shù)時圓柱后方有一對尾渦; 隨著Re的增加, 尾跡出現(xiàn)震蕩; 繼續(xù)增加Re, 近尾跡中渦交替脫落, 進(jìn)入下游, 形成持續(xù)較長距離的穩(wěn)定渦列, 即Karman渦街[36]. 本文選取Re=1 200 時作為典型計算條件, 對應(yīng)的計算狀態(tài)如表1所示, 氣體為氮氣.

表1 低速Karman渦街算例的典型計算狀態(tài)Table 1 Classical computing condition for low-speed Karman vortex street case

2.2 圓柱尾跡區(qū)的Karman渦街

圓柱直徑d=1 m, 數(shù)值計算過程中流場網(wǎng)格關(guān)于x軸對稱(如圖1所示), 物理空間的網(wǎng)格量為(301×201), 分子速度空間的離散速度點為(16×16), 關(guān)于速度空間的矩積分采用Gauss-Hermite積分法. 圓柱尾跡區(qū)Karman渦街的形成、發(fā)展、運動及耗散過程如圖2所示, 圖中給出的是不同時刻流場壓力等值線云圖. 初始時刻流場完全處于來流狀態(tài), 壁面擾動從t=0時刻開始傳播. 對于低速流動, 擾動向全流場傳播. 在圖2(a)所示的t=0.014 s 時刻, 物面邊界擾動已經(jīng)傳播到離圓柱較遠(yuǎn)的距離, 圓柱周圍流場基本成形. 此時能夠明顯看到圓柱頭部區(qū)域的高壓、 圓柱上下兩側(cè)的低壓以及圓柱后方的擾動區(qū)都是上下對稱的. 到圖2(b)所示的t=0.024 s時刻, 圓柱后端尾跡區(qū)出現(xiàn)了上下對稱的兩個低壓旋渦結(jié)構(gòu). 可以看到, 低壓旋渦結(jié)構(gòu)的產(chǎn)生是由逆壓梯度導(dǎo)致的流動分離. 分離出的旋渦持續(xù)向下游流動, 并逐漸發(fā)展、擴(kuò)大. 到圖2(c)所示的t=0.057 s時刻, 尾跡區(qū)兩個渦在向下游流動的同時相互接觸、互相影響, 但此時兩個旋渦依然是對稱結(jié)構(gòu). 而在圖2(d)的t=0.569 s時刻, 尾跡區(qū)兩個旋渦結(jié)構(gòu)“對稱破缺”, 上下兩個渦交替“擴(kuò)大—收縮”, 并繼續(xù)向下游流動. 這種對稱性被打破的現(xiàn)象來源于流場的非穩(wěn)態(tài)演化, 與外界干擾、網(wǎng)格結(jié)構(gòu)等并無關(guān)系. 隨著旋渦結(jié)構(gòu)的持續(xù)運動、發(fā)展、相互作用, 到圖2(e)的t=0.687 s時刻形成了“擬序結(jié)構(gòu)”的尾跡區(qū)旋渦序列, 并最終形成了圖2(f), (g)所示“Karman渦街”非定常流動狀態(tài). 其中, 圖2(f)和(g)分別是半個周期的 Karman渦街運動耗散過程, 旋渦結(jié)構(gòu)呈現(xiàn)上下交替演化現(xiàn)象, 從圖中可以看到Karman渦街脫落、發(fā)展、運動、耗散的整個過程.

圖1 低速圓柱繞流的物理空間網(wǎng)格Fig. 1 Grid system in physical space for low-speed flow past a cylinder

(a) t=0.014 s

圖3給出了圓柱繞流尾跡區(qū)在半個周期內(nèi)的Karman渦街交替脫落現(xiàn)象. 從流場的溫度等值線云圖中能夠更明顯地看到圓柱上下頂點稍靠后的位置交替“甩出”渦旋結(jié)構(gòu), 及其渦脫離的整個過程. 有意思的是, 甩出的渦旋結(jié)構(gòu)為低壓、低溫, 而另一側(cè)對應(yīng)被壓縮的區(qū)域則是高壓、高溫. 產(chǎn)生的Karman渦街在離圓柱尾部4d距離以后逐漸被耗散掉. 圖4 詳細(xì)繪出了半個周期內(nèi)的流線圖, 展示了圓柱尾跡區(qū)下半部甩出旋渦到上半部甩出旋渦的交替演化過程.

圖3 圓柱繞流尾跡區(qū)半個周期內(nèi)的Karman渦街交替脫落現(xiàn)象(溫度)Fig. 3 Alternately dropping phenomenon of the Karman vortex street in half period at the wake-region of flow past a cylinder(temperature contours)

圖4 圓柱繞流尾跡區(qū)半個周期內(nèi)的流線圖Fig. 4 Streamlines in half period at the wake-region of flow past a cylinder

圖5給出了本文GKUA得到圓柱上下表面的壓力分布與連續(xù)流區(qū)N-S方程數(shù)值計算結(jié)果的比較, 其中N-S方程計算中分別采用了層流模型和兩方程k-ε湍流模型(RANS). 可以看到, GKUA得到的表面壓力分布規(guī)律與連續(xù)流區(qū)經(jīng)典N-S方程基本一致. 對于該低速流動, 迎風(fēng)面壓力值略高于背風(fēng)區(qū), 且GKUA得到的迎風(fēng)面壓力分布與層流模型和湍流模型都吻合較好. 背風(fēng)面壓力值的差別相對來說要大一些, 但也在可控范圍內(nèi)(最大偏差低于2.3%). 同時可以看出, 上下表面壓力分布的不對稱性正體現(xiàn)出了流動的非定常特性, 這跟尾跡區(qū)渦的交替脫落緊密相關(guān).

(a) Upper-half part of the cylinder surface

該狀態(tài)下圓柱繞流數(shù)值計算得到的Strouhal數(shù)與文獻(xiàn)[37]中的實驗值比較如表2所示. 可以看到, 本文數(shù)值計算得到的表征Karman渦街非定常流動特征的Strouhal數(shù)與實驗值和公認(rèn)理論值0.21都是非常接近的. 考慮到數(shù)據(jù)量的問題, 在數(shù)值計算過程中并沒有保存任意時刻的流場結(jié)果, 而是隔一定時間保存一次, 導(dǎo)致Karman渦街的周期估算值存在一定偏差. 鑒于此, 本次圓柱繞流尾跡區(qū)Karman渦街的數(shù)值模擬結(jié)果是可信的.

表2 圓柱繞流計算得到的Strouhal數(shù)與文獻(xiàn)比較Table 2 Comparison of the Strouhal numbers from simulation and experiment

Karman渦街的尺寸比定義為b/a, 該參數(shù)可以用來描述渦列, 其中a為同渦列中相鄰旋渦之間的距離,b為兩列旋渦的間隔(見圖6). von Karman[38]給出渦列穩(wěn)定的必要條件為尺寸比b/a≈0.281. 雖然后人指出Karman渦街穩(wěn)定的尺寸比并不嚴(yán)格等于 0.281[39], 但不同計算條件下相差不多. 本文數(shù)值計算圓柱Karman渦街的尺寸比b/a=0.283, 與 von Karman理論推導(dǎo)的結(jié)果非常接近(見表3), 進(jìn)一步驗證了氣體動理論統(tǒng)一算法模擬Karman渦街等非定常流動問題是適用的.

圖6 Karman渦街尺寸比示意圖Fig. 6 Sketch of the length ratios for Karman vortex street

表3 圓柱Karman渦街尺寸比的計算值與理論值比較Table 3 Comparison of the length ratios from simulation and theory analysis

2.3 非對稱尖劈產(chǎn)生的渦街

圓柱繞流的Karman渦街是由流動的對稱性破缺引起的, 而非對稱尖劈繞流則由于流場本身的非對稱性, 誘導(dǎo)產(chǎn)生的Karman渦街在空間上是非對稱的, 但在時間上呈現(xiàn)出周期變化的規(guī)律. 非對稱尖劈的物理空間流場網(wǎng)格如圖7所示, 規(guī)模為58 403個網(wǎng)格點, 分子速度空間采用(32×32)的Gauss-Hermite積分法. 圖8給出了一個周期內(nèi)非對稱尖劈誘導(dǎo)渦街的流場壓力等值線云圖, 其中t=0表示該周期的起始時刻, 并不是初始時刻. 可以看到, 上下尖劈誘導(dǎo)產(chǎn)生的渦系結(jié)構(gòu)在向下游發(fā)展時是非常規(guī)整的, 形成了相互干擾又彼此獨立的擬序結(jié)構(gòu). 在圖8(a)該周期開始時刻, 上部尖劈出現(xiàn)誘導(dǎo)渦分離現(xiàn)象; 到圖8(b)Tp/4時刻上部誘導(dǎo)渦發(fā)育完全, 脫離尖劈向下游運動; 到圖8(c)Tp/2時刻下部尖劈開始產(chǎn)生誘導(dǎo)渦; 到圖8(d)3Tp/4時刻下部誘導(dǎo)渦發(fā)育完全, 也向下游運動; 到圖8(e)t=Tp時刻上部尖劈又產(chǎn)生新的誘導(dǎo)渦. 如此循環(huán)往復(fù), 構(gòu)成了非對稱尖劈誘導(dǎo)Karman渦街的非定常流動過程.

圖9給出了一個周期內(nèi)t=Tp/2,t=Tp兩個時刻的非對稱尖劈繞流下游的渦系結(jié)構(gòu)發(fā)展, 其中展示的是流線圖. 可以看到, 非對稱尖劈下游渦系結(jié)構(gòu)比圓柱繞流下游的擬序結(jié)構(gòu)要復(fù)雜得多, 前者上下尖角誘導(dǎo)產(chǎn)生的渦核大小略有差別, 且相隔半個周期的渦系結(jié)構(gòu)并不關(guān)于y軸對稱. 這些差異是由尖劈的非對稱性引起的.

圖7 非對稱尖劈繞流的物理空間網(wǎng)格Fig. 7 Grid system in physical space for flow past an asymmetric wedge

(a) t=0

3 結(jié)論

本文在氣體動理論統(tǒng)一算法框架下, 應(yīng)用離散速度坐標(biāo)法和3階WENO格式, 以及3階Runge-Katta時間推進(jìn)法, 數(shù)值求解考慮轉(zhuǎn)動自由度影響的Rykov模型方程, 得到了能模擬雙原子氣體跨流域非定常流動問題的一種確定論數(shù)值計算方法. 對經(jīng)典的二維Karman渦街非定常流動現(xiàn)象的模擬, 說明了該跨流域非定常模擬方法能夠適用于連續(xù)流區(qū)的低速流動, 并復(fù)現(xiàn)了關(guān)于Karman渦街流動現(xiàn)象的一些規(guī)律和理論:

(1)低速圓柱繞流的Karman渦街是由于流動演化的對稱性破缺而產(chǎn)生周期性變化規(guī)律, 與圓柱外形、流場網(wǎng)格等無關(guān);

(2)圓柱繞流的Karman渦街不僅在時間上有周期性變化規(guī)律, 而且圓柱上下誘導(dǎo)產(chǎn)生的擬序結(jié)構(gòu)在空間上也具有對稱性, 這種空間對稱性并不體現(xiàn)在同一時刻, 而是體現(xiàn)在半個周期上的空間對稱;

(3)圓柱繞流的誘導(dǎo)渦產(chǎn)生必然是發(fā)生在圓柱上下頂點靠后的逆壓梯度區(qū);

(4)非對稱尖劈誘導(dǎo)產(chǎn)生的Karman渦街并不具有空間對稱性, 但依然表現(xiàn)出時間上的周期性變化, 且其渦系結(jié)構(gòu)更加復(fù)雜.

致謝本工作得到973計劃(2014CB744100), 國家自然科學(xué)基金(11902339, 11325212)資助.