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

(1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速空氣動(dòng)力研究所, 四川綿陽 621000; 2. 國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室, 北京 100191)

引 言

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

(1)

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

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

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

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

1 控制方程及計(jì)算方法

1.1 控制方程

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

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

(2)

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

物理空間宏觀流動(dòng)參數(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)]

對(duì)于氮?dú)釴2而言,T*=91.5 K,δ=1/1.55.

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

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

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

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

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

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

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

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

(3)

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

1.3 適于非定常流動(dòng)計(jì)算的差分方法

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

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

(4)

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

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

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

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

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

2 Karman渦街形成及運(yùn)動(dòng)的數(shù)值模擬

2.1 問題描述

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

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

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

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

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

(a) t=0.014 s

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

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

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

(a) Upper-half part of the cylinder surface

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

表2 圓柱繞流計(jì)算得到的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], 但不同計(jì)算條件下相差不多. 本文數(shù)值計(jì)算圓柱Karman渦街的尺寸比b/a=0.283, 與 von Karman理論推導(dǎo)的結(jié)果非常接近(見表3), 進(jìn)一步驗(yàn)證了氣體動(dòng)理論統(tǒng)一算法模擬Karman渦街等非定常流動(dòng)問題是適用的.

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

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

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

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

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

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

(a) t=0

3 結(jié)論

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

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

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

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

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

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