劉 躍管小榮徐 誠

(南京理工大學(xué) 機械工程學(xué)院,南京 210094)

0 引 言

廣泛存在于航空航天、風(fēng)工程等空氣動力學(xué)領(lǐng)域的鈍體繞流問題通常為無序、非線性、多尺度的湍流分離運動,在使用數(shù)值模擬方法(CFD)對其進(jìn)行分析計算時,湍流的求解成為影響CFD計算精度的瓶頸問題[1]。在諸多的湍流模型中,兼顧計算精度與資源消耗的RANS/LES混合方法受到人們的青睞[2],混合模型在近壁區(qū)及遠(yuǎn)場使用RANS模型,而在分離區(qū)及尾跡區(qū)切換為更合理的類LES濾波模型,極大地提高了計算精度,且資源消耗遠(yuǎn)小于LES僅略大于RANS。近年來,不同學(xué)者發(fā)展了種類繁多的混合方法,其中包括DES(Detached Eddy Simulation)方法[3]、SAS(Scale Adaptive Simulation)方 法[4-7]、PANS(Partially Averaged Navier-Stokes)方法[8]、PITM(Partially Integrated Transport Model)方法、LNS(Limited Numerical Scales)方法等。其中DES方法是最早出現(xiàn)的混合方法,經(jīng)過一系列改進(jìn)(DDES、IDDES)已發(fā)展較為成熟[9-13]；SAS模型為能夠解析寬頻非定常湍流脈動而對計算網(wǎng)格依賴較小的新一代非定常湍流預(yù)測方法(URANS),通過加入?yún)⒄债?dāng)?shù)亓鲃拥鸟T卡門(von Karman)長度尺度L vk,合理地釋放更多流場信息[14-15]；PANS模型同樣是最近發(fā)展的一種RANS-DNS橋接模型,模型構(gòu)造相對簡單,通過加入?；膭幽鼙壤齾?shù)f k調(diào)整湍流的?；縼碚{(diào)整湍流求解尺度[16]；DES模型的提出者Spalart在2000年時曾預(yù)測,2045年左右計算機的計算能力才能滿足利用LES解決90%的湍流問題,因此,在現(xiàn)階段,驗證與發(fā)展RANS/LES方法仍是非常必要的[17-18]。

鈍體繞流中的圓柱繞流計算模型簡單,同時又包含如湍流壁面、剪切層失穩(wěn)、流動分離、周期性渦脫落等豐富的流動現(xiàn)象,是考察湍流模型計算性能的標(biāo)準(zhǔn)算例。Re＝3900的圓柱繞流是典型的亞臨界雷諾數(shù)流動,流動結(jié)構(gòu)復(fù)雜且耗費資源較少,人們對該雷諾數(shù)的圓柱繞流進(jìn)行了大量的實驗[19-22]及數(shù)值研究[23-29],有可靠豐富的數(shù)據(jù)可供參考。本文選用SAS模型及兩種變f k構(gòu)造函數(shù)的PANS模型計算了Re＝3900的圓柱繞流流場,并重點從流場特征捕捉、氣動力特性計算、渦黏性控制等方面對SAS模型及PANS模型進(jìn)行了比較,通過對比標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格后的計算結(jié)果,比較了兩類混合方法對網(wǎng)格的適應(yīng)性,并探討了兩類混合模型的尺度函數(shù)對流場的調(diào)控機制。

1 計算模型及數(shù)值方法

1.1 控制方程

控制方程為三維N-S方程,在一般計算坐標(biāo)系(τ,ξ,η,ζ)下的 守恒形 式為:

式中:為守恒變量為無粘通量為黏性通量,各物理量的具體表述參見文獻(xiàn)[1]。

1.2 湍流模型

1.2.1 SST-k-ω模型

SST模型由標(biāo)準(zhǔn)k-ω模型與k-ε模型混合修改而成,該模型兼顧了k-ω模型的近壁性能及k-ε模型的遠(yuǎn)場準(zhǔn)確性,SST模型中湍動能k和比耗散率ω的輸運方程可表示如下:

式中模型參數(shù)σk、α、β、σω、σω2均由對應(yīng)的k-ω模型參數(shù)φ1和k-ε模型參數(shù)φ2通過一定函數(shù)計算得到,其對應(yīng)的函數(shù)關(guān)系為:

SST模型中湍動黏度μt可表示為:

式中S為應(yīng)變率張量S ij的模,SST模型中F1、F2函數(shù)構(gòu)造方式及各模型參數(shù)值參考文獻(xiàn)[30]。

1.2.2 SST-SAS模型

尺度自適應(yīng)SAS(Scale-Adaptive Simulation)模型是一種較新的創(chuàng)新型湍流模型,它通過在湍流模型中加入刻畫當(dāng)?shù)亓鲃油負(fù)涞膙on Karman長度尺度L vk來作為湍流模型的第二長度尺度[6]。基于SST模型的SAS模型構(gòu)造方法為在式(3)SST模型的ω方程中加入發(fā)揮自適應(yīng)作用的源項QSAS,其定義式為:

由式(6)可以看出:在SAS模型中出現(xiàn)了由一階及二階速度導(dǎo)數(shù)決定的長度尺度項L vk,該長度尺度基于當(dāng)?shù)亓鲃佣c網(wǎng)格尺度無關(guān),在湍流邊界層內(nèi),L vk能?；瘧T性子區(qū)所有湍流脈動,同時,在非穩(wěn)態(tài)區(qū)域L vk能根據(jù)當(dāng)?shù)鼐W(wǎng)格分辨的湍流渦動態(tài)調(diào)整RANS長度尺度。在SAS模型中,湍動黏度μt仍按照式(5)計算。

1.2.3 SST-PANS模型

本文中所用PANS模型同樣是基于SST模型修改而成[8],其原理是在SST模型中引入兩個模型控制參數(shù)f k和fω:

式中:f k為?；膭幽芘c總的湍動能之比,fω是模化比耗散率與總的比耗散率的比值。如果假定f k和fω是不隨時間和空間變化的常數(shù),則k和ω的輸運方程(2)、(3)變?yōu)?

將式(8)、(9)與式(2)、(3)對比可以發(fā)現(xiàn)在PANS模型中主要是模型參數(shù)發(fā)生了變化,PANS模型中模型參數(shù)σku、σωu、σω2u、β表示為:

通過以上變化,將SST模型轉(zhuǎn)變?yōu)镾ST-PANS模型,由k、ε、ω的關(guān)系可得式(10)中fω與f k和fε的關(guān)系,在SST-PANS模型中,fε一般取為1,這樣PANS中的參數(shù)可以表示為f k的函數(shù)。通過調(diào)節(jié)f k的值(1-0),改變模化的湍動能與比耗散率的比值,從而實現(xiàn)RANS-DNS的轉(zhuǎn)化。這時湍動黏度μtu變?yōu)?

在早期PANS模型的應(yīng)用中,f k通常被設(shè)置為空間常數(shù),通過計算發(fā)現(xiàn)對于整個流場很難確定合適統(tǒng)一的f k值。對于復(fù)雜的分離流動問題,計算時分離區(qū)與尾跡區(qū)需要較小的f k值(f k＜1)釋放湍流信息,而近壁區(qū)與遠(yuǎn)場區(qū)采用RANS模型(f k＝1)更合理,因此近年來人們逐漸嘗試用適當(dāng)?shù)暮瘮?shù)來計算f k分布,本文選取其中兩種相對簡單、計算相對穩(wěn)定的變f k函數(shù)來進(jìn)行對比。

(a)類DES型f k函數(shù)[31-32]。該函數(shù)參考DES模型中模型尺度的構(gòu)造方法,將f k表示為網(wǎng)格尺度Δ與湍流尺度L u的函數(shù):

式中模型常數(shù)CPANS取0.3。用該函數(shù)進(jìn)行計算時,在分離區(qū)與尾跡區(qū),網(wǎng)格尺度要小于湍流尺度,因此f k＜1；在壁面及遠(yuǎn)場區(qū)域,網(wǎng)格尺度大于湍流尺度,因此f k＝1。文中將使用類DES型f k分布函數(shù)的PANS模型稱為PANS(des)模型。

(b)Tanh型f k函數(shù)[33]。在該函數(shù)中,f k可表示為:

對比式(12)與式(13)可知,Tanh型f k函數(shù)由湍流尺度、網(wǎng)格尺度及tanh函數(shù)修改而成,tanh函數(shù)的取值范圍為(-1,1),這樣根據(jù)流動與網(wǎng)格關(guān)系得到不同位置的f k取值范圍為(0,1)。文中將使用Tanh型f k函數(shù)的PANS模型稱為PANS(tanh)模型。

2 計算模型及數(shù)值方法

2.1 計算域、網(wǎng)格及邊界設(shè)置

文中圓柱繞流計算雷諾數(shù)為Re＝3900(基于圓柱直徑D與遠(yuǎn)場來流速度U∞),圖1為計算域及邊界條件示意圖。為減小上游來流對計算的影響并使尾跡區(qū)流動充分發(fā)展,上游來流區(qū)域(x方向)及橫向區(qū)域(y方向)延展至20D,下游計算區(qū)域設(shè)置為30D。根據(jù)近尾跡區(qū)流向結(jié)構(gòu)的展向波長設(shè)置計算域展向長度(z方向)為πD[2]。圓柱體設(shè)置為無滑移壁面,上游入口流動各物理量由無擾動自由來流條件給定,下游出口設(shè)置為無反射邊界,計算域的上下界面設(shè)置為自由滑移壁面,計算域的前后邊界設(shè)置為周期邊界。文中坐標(biāo)系設(shè)置如圖1右下角所示,軸向從上游到下游方向為x正方向,坐標(biāo)系符合右手螺旋法則,橫向為y方向,展向為z方向,方位角從上游圓柱前駐點為0°開始,順時針旋轉(zhuǎn)為壁面方位角θ的正方向,下游背風(fēng)點為180°。

圖1 計算域及邊界條件Fig.1 Computational domain and boundary conditions

圖2為圓柱繞流的局部網(wǎng)格示意,計算網(wǎng)格為分塊結(jié)構(gòu)化網(wǎng)格。為了保證計算精度,在壁面附近及尾跡分離區(qū)域?qū)W(wǎng)格進(jìn)行適當(dāng)加密,壁面第一層網(wǎng)格高度為0.002D,保證壁面y＋＜1。為了探討不同模型的尺度函數(shù)對網(wǎng)格的敏感性,文中設(shè)置了兩套計算網(wǎng)格,其中標(biāo)準(zhǔn)網(wǎng)格保證近壁增長率為1.06,分離區(qū)網(wǎng)格固定為0.1D；加密網(wǎng)格為保證近壁增長率為1.04,分離區(qū)網(wǎng)格固定為0.05D。參考文獻(xiàn)[24,25]中的網(wǎng)格設(shè)置見表1。

圖2 圓柱繞流局部網(wǎng)格Fig.2 Computational grid for the cylinder

表1 計算網(wǎng)格Table 1 Computational grids

2.2 方程求解設(shè)置

文中對所求方程組空間離散采用有限體積法(FVM),對流項中速度項的離散采用三階限制(TVD)中心差分格式,湍流項離散采用二階限制中心差分格式,梯度項及黏性項的離散采用四階中心差分格式,壓力與速度耦合選用PISO算法。為了保證計算的收斂性及穩(wěn)定性,時間項離散采用二階隱式后退歐拉(Backward)格式,并采用自適應(yīng)時間步長方法[24],時間步長設(shè)置為庫朗數(shù)(Courant number)小于0.5,并且保證在一個渦脫落周期內(nèi)包含近500個時間步。

3 計算結(jié)果分析

3.1 數(shù)據(jù)采集及流場宏觀結(jié)構(gòu)

文中分別計算了SST、SAS、PANS(des)、PANS(tanh)模型的圓柱繞流流場。周期性渦脫落是該計算雷諾數(shù)下的顯著特征。為了獲得合理的統(tǒng)計數(shù)據(jù),根據(jù)Franke等[34]的研究結(jié)果,圓柱繞流的統(tǒng)計周期應(yīng)大于40個脫落周期。因此參照文獻(xiàn)[34,25],本文中總的計算時間約為60個渦脫落周期,在對非定常結(jié)果的統(tǒng)計時去掉前10個脫落周期,將后50個周期內(nèi)的統(tǒng)計數(shù)據(jù)作為有效數(shù)據(jù)進(jìn)行模型的對比。計算結(jié)果進(jìn)行時間與空間平均。圖3中給出了以標(biāo)準(zhǔn)網(wǎng)格計算結(jié)果為例的圓柱繞流升力系數(shù)(C l)、阻力系數(shù)(C d)的時程(U∞/D無量綱化)變化曲線?？梢钥闯鯮e＝3900時的升力幅值已不再穩(wěn)定,但仍呈現(xiàn)出周期性變化趨勢,同時阻力系數(shù)表現(xiàn)出強烈的脈動現(xiàn)象,并在一平均值附近上下擺動。

圖3 升力、阻力系數(shù)隨時間變化及升力頻譜Fig.3 Time history of lift C l,drag C d coef ficients

圖4中給出了由升力系數(shù)時程曲線經(jīng)過傅里葉變換(FFT)獲得的渦脫落頻率功率譜,由此可以計算渦脫落頻率f及斯特勞哈數(shù)S t(S t＝f D/U∞)。各模型阻力系數(shù)(C d)及斯特勞哈數(shù)S t見表2。

圖4 升力頻譜Fig.4 Spectra of lift

圖5為四種模型基于標(biāo)準(zhǔn)網(wǎng)格計算得到的圓柱繞流流場二維瞬時(t＝300)渦量分布及用軸向速度染色的三維Q分布,對比可知:由于RANS模型湍流黏性較大,抹去了很多重要的流場細(xì)節(jié),計算結(jié)果只捕捉到大尺度的脫落渦結(jié)構(gòu),Q圖也只先顯示出近似二維的流動結(jié)構(gòu),這說明RANS模型基本上沒有求解小尺度流動的能力；SAS模型及PANS模型均捕捉到豐富的流動結(jié)構(gòu),二維渦量及三維Q值均顯示出復(fù)雜的周期性渦脫落現(xiàn)象。分析三者的Q圖可清晰地發(fā)現(xiàn)附著于壁面的二維層流剪切層及分離區(qū)流向與展向相嵌套的三維無序渦體結(jié)構(gòu)。SAS模型與PANS模型均允許形成湍流譜,表現(xiàn)出預(yù)期的類似LES模型的求解小尺度運動的能力。同時,仔細(xì)對比SAS模型與PANS模型的結(jié)果可以看出,PANS模型能捕捉到比SAS模型更豐富的湍流結(jié)構(gòu)。

圖5 不同模型計算得到的流場瞬時渦量與Q值圖(t＝300)Fig.5 Instantaneous vorticity and Q of different models(t＝300)

由實驗及以往計算結(jié)果[23]可知,圓柱繞流壁面剪切層的發(fā)展對外界影響因素非常敏感,微小外部干擾很容易造成邊界層的提前分離,使得剪切層變短,進(jìn)而影響回流區(qū)及整個流場的測量或計算結(jié)果。圖6為用渦量等值線表示的標(biāo)準(zhǔn)網(wǎng)格條件下四種模型的瞬時(t＝300)剪切層,通過對比可知剪切層長度關(guān)系為:SAS＞PANS(des)＞PANS(tanh)＞SST。

圖6 不同模型瞬時分離剪切層示意圖(t＝300)Fig.6 Instantaneous separating shear layers of different models(t＝300)

3.2 回流區(qū)速度分布

圖7為以SAS模型為例的y＝0平面軸向速度時均圖,可以清晰地看出在圓柱近尾跡區(qū)域形成一個流速與來流相反的區(qū)域,即圓柱繞流的回流區(qū)。圖8為圓柱下游尾跡區(qū)中心線上時均軸向速度分布。由圖可知,圓柱背風(fēng)點處速度為0,速度為逐漸增大的負(fù)值,到達(dá)回流區(qū)中心時速度負(fù)值最大。背風(fēng)點與空間流動的再附點即流向速度恢復(fù)到0的點之間的距離即為回流區(qū)長度(L r/D)。四種模型均模擬出了圓柱后的速度分區(qū),但彼此間差別較大,特別是對回流區(qū)的模擬。圖8(b)圖是對回流區(qū)中心線時均軸向速度的放大圖。首先需要說明的是用空心圓圈表示的實驗數(shù)據(jù)是Lourenco&Shih(1993)用PIV測量獲得[20],用實心圓表示的實驗數(shù)據(jù)是Parnaudeau等(2008)同樣用PIV測量獲得[21],可以看到兩者實驗測得的回流區(qū)長度差別較大。最近幾年的文獻(xiàn)也出現(xiàn)了與Parnaudeau實驗接近或與Lourenco實驗接近的兩種數(shù)值計算結(jié)果,有意思的是Kravchenko等(2000)在用B-Spline伽遼金方法對Re＝3900進(jìn)行的LES數(shù)值計算中[23],最先發(fā)現(xiàn)了數(shù)值結(jié)果與Lourenco實驗數(shù)據(jù)在回流區(qū)與近尾跡區(qū)速度型的差別。Kravchenko計算的回流區(qū)長度大于Lourenco實驗值,其當(dāng)時推測可能有外加擾動導(dǎo)致Lourenco實驗中圓柱壁面剪切層過早分離,造成剪切層變短,進(jìn)而影響回流區(qū)長度。Parnaudeau實驗數(shù)據(jù)在回流區(qū)及速度型方面均驗證了Kravchenko的數(shù)值計算結(jié)果。由于實驗時間較晚,實驗方法更先進(jìn)、嚴(yán)謹(jǐn),本文中我們認(rèn)為Parnaudeau的實驗值是更可信的。分析圖8(b)可知,SST模型計算的圓柱回流區(qū)長度誤差較大,SAS與PANS計算的回流區(qū)長度關(guān)系為SAS＞PANS(des)＞PANS(tanh)。三種模型計算的回流區(qū)區(qū)域位于兩個實驗值之間,與Parnaudeau實驗值仍存在一定誤差,其中SAS模型與Parnaudeau實驗最接近。以上結(jié)論與圖6中得出的剪切層長度結(jié)論也是一致的。

圖7 y＝0平面軸向速度時均分布Fig.7 Mean streamwise velocity along the y＝0 surface

圖8 圓柱尾跡中心線上時均軸向速度分布Fig.8 Mean streamwise velocity along the wake center line

圖9～圖11為圓柱繞流下游不同截面處的時均軸向速度分布(圖例見圖8)。可以看到軸向速度流型從近尾跡區(qū)到遠(yuǎn)尾跡區(qū)由U型向V型轉(zhuǎn)變,結(jié)合前部分析可知,剪切層過早分離同樣會導(dǎo)致速度型的過早轉(zhuǎn)變,圖中兩種實驗值的對比也說明了這一點。圖9近壁x/D＝0.58時,四種模型的速度型均為U型,彼此差別不大,而到回流區(qū)x/D＝1.06時,四種模型出現(xiàn)差別,SST模型為完全的V型,SAS模型與兩種PANS模型同樣沒有計算出更加正確的U型,而是介于U型與V型之間,其中SAS模型與Parnaudeau實驗得到的U型更接近,同時PANS(des)要稍好于PANS(tanh)；在近尾跡回流區(qū)后部x/D＝1.54、2.02位置處,除了SST沒有計算出正確的速度型,SAS模型與兩種PANS模型的計算結(jié)果只是速度峰值與實驗值存在一定誤差,整體均與實驗值吻合較好；在遠(yuǎn)尾跡區(qū)x/D＝4、7、10位置處,四種模型計算的速度型均較好,速度峰值與實驗存在較小誤差。以上分析說明近尾跡區(qū)域特別是回流區(qū)流動較為復(fù)雜、敏感,計算控制更困難,同時遠(yuǎn)尾跡區(qū)受近壁區(qū)影響較小,使得遠(yuǎn)尾跡區(qū)計算誤差小于近尾跡區(qū)。

圖9 不同截面的平均軸向速度分布(x/D＝0.58,1.06)Fig.9 Mean streamwise velocity at different surfaces(x/D＝0.58,1.06)

圖10 不同截面的平均軸向速度分布(x/D＝1.54,2.02)Fig.10 Mean streamwise velocity at different surfaces(x/D＝1.54,2.02)

圖11 不同截面的平均軸向速度分布(x/D＝4,7,10)Fig.11 Mean streamwise velocity at different surfaces(x/D＝4,7,10)

圖12為圓柱繞流下游近尾跡回流區(qū)附近不同截面處的時均橫向速度分布(圖例見圖8)。由x/D＝1.06處的時均橫向速度與實驗值可以看到,SAS與兩種PANS模型得到的速度型均與Parnaudeau實驗值較吻合,與SST模型比較顯示出較高的計算精度。x/D＝1.54、2.02處的對比結(jié)果可以看出四種模型在速度型與峰值上均存在一定誤差,與軸向速度計算結(jié)果類似。

圖12 不同截面的平均橫向速度分布(x/D＝1.06,1.54,2.02)Fig.12 Mean cross-flow velocity at different surface(x/D＝1.06,1.54,2.02)

3.3 加密網(wǎng)格計算結(jié)果

根據(jù)混合模型的構(gòu)造原理可知,其湍流尺度求解能力是基于對渦黏性的控制作用。圖13為三種混合模型計算的圓柱繞流時均湍動黏度μt/μ(湍動黏度與來流黏度之比)分布,上部為標(biāo)準(zhǔn)網(wǎng)格計算結(jié)果,下部為加密網(wǎng)格計算結(jié)果。僅從數(shù)值來看,SAS模型與兩種PANS模型的湍動黏度分布關(guān)系為SAS＞PANS(des)＞PANS(tanh)。PANS模型所求湍動黏度要小于SAS模型,這也解釋了圖4中PANS模型對流場的捕捉要好于SAS模型。但需要注意的是,較低的湍動黏度不一定獲得較準(zhǔn)確的數(shù)值計算結(jié)果。由文中氣動力的比較可知,SAS模型要好于PANS模型。通過比較網(wǎng)格加密后的湍流黏度可知,兩種PANS模型在近壁網(wǎng)格密集區(qū)域出現(xiàn)了湍流黏度較低的區(qū)域,這是由于PANS模型的f k構(gòu)造函數(shù)仍是由當(dāng)?shù)赝牧鞒叨群途W(wǎng)格尺度來控制,在近壁網(wǎng)格較密的條件下,PANS模型顯現(xiàn)出類似DES模型對網(wǎng)格敏感的缺點,這時近壁RANS區(qū)域被破壞,而網(wǎng)格尺度還沒有小到能夠進(jìn)行大渦模擬的程度。

圖13 時均湍流黏度比(μt/μ)分布(標(biāo)準(zhǔn)網(wǎng)格：上,加密網(wǎng)格：下)Fig.13 Time-averaged distributions ofμt/μ(baseline grid：top,fine grid：bottom)

圖14為三種混合模型在不同網(wǎng)格條件下計算的局部(x/D＝1.06)時均軸向雷諾應(yīng)力分布(圖例見圖15)?？梢钥吹絊AS模型及兩種PANS模型計算的雷諾應(yīng)力與實驗值均有較大差距,這可能是由于渦黏性模型的內(nèi)在缺陷造成的。我們在此關(guān)注不同的模型在不同網(wǎng)格條件下的計算差異?？梢钥吹絊AS模型在標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格條件下均高估了雷諾應(yīng)力值,但兩種網(wǎng)格下的計算結(jié)果差別不大；兩種PANS模型對不同的網(wǎng)格表現(xiàn)出相似的規(guī)律,在加密網(wǎng)格時雷諾應(yīng)力均出現(xiàn)了明顯的降低,結(jié)合圖13中的黏度分布可知,由于近壁區(qū)域網(wǎng)格的加密,湍流黏度降低,導(dǎo)致?；睦字Z應(yīng)力減少,同時有沒有求解足夠的雷諾應(yīng)力,導(dǎo)致了雷諾應(yīng)力不足現(xiàn)象的發(fā)生,這也驗證了PANS構(gòu)造函數(shù)對網(wǎng)格較敏感的缺陷。同時兩種PANS模型函數(shù)比較可以看出,PANS(des)模型計算的雷諾應(yīng)力及氣動力分布要稍優(yōu)于PANS(tanh)模型,這也說明了PANS(tanh)模型構(gòu)造函數(shù)對網(wǎng)格的敏感程度要大于PANS(des)模型的構(gòu)造函數(shù)。

圖14 局部時均軸向雷諾應(yīng)力分布Fig.14 Time-averaged distributions of streamwise Reynolds stress

圖15為三種模型在不同網(wǎng)格條件下的時均軸向速度分布對比。結(jié)合前文分析可知,SAS模型網(wǎng)格獨立性較好,標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格下的計算結(jié)果沒有體現(xiàn)出太大差別,且由于對剪切層的計算更準(zhǔn)確,計算得到的回流區(qū)大小更接近實驗值。兩種PANS模型由于對網(wǎng)格較敏感,在加密網(wǎng)格時延緩了剪切層的轉(zhuǎn)捩,導(dǎo)致回流區(qū)變長。這在PANS(tanh)模型的計算結(jié)果中表現(xiàn)更明顯。

圖15 時均軸向速度分布Fig.15 Time-averaged distributions of streamwise velocity

3.4 積分統(tǒng)計量比較

表2為四種模型計算的圓柱繞流平均積分統(tǒng)計量與相關(guān)實驗值及數(shù)值結(jié)果對比。可以看到在與剪切層運動密切相關(guān)的參數(shù)中,普通網(wǎng)格下SAS與兩種PANS模型均取得了與實驗值較吻合的回流區(qū)長度(L r/D)計算結(jié)果,與構(gòu)造相對完善的LES-THE、SA-IDDES、V2-f-DES模型的計算結(jié)果也是相近的,其中SAS模型的計算結(jié)果與實驗值最接近；加密網(wǎng)格時,兩種PANS模型均過大地預(yù)測了回流區(qū)長度。在阻力系數(shù)(C d)及分離角(θSep)的計算中,在標(biāo)準(zhǔn)網(wǎng)格下三種混合模型計算誤差均較小,C d滿足實驗值的0.98±0.05區(qū)間,在加密網(wǎng)格時,SAS模型的計算結(jié)果變化不大,兩種PANS模型獲得了較小的分離角及阻力系數(shù)。在其余的比較項中,可以看到四種模型對圓柱尾跡旋渦脫落頻率的計算及壁面壓力的計算均較準(zhǔn)確,所計算的斯特勞哈數(shù)(S t)與背風(fēng)點壓力系數(shù)(-C pb)差別不大。對于SST模型,由于其構(gòu)造?；怂械耐牧髁?導(dǎo)致重要的流場細(xì)節(jié)被抹平,其計算的阻力系數(shù)、背風(fēng)點壓力系數(shù)、回流區(qū)長度、分離角均存在較大誤差。同時可以看到兩種網(wǎng)格條件下PANS(des)計算結(jié)果稍好于PANS(tanh)。比較文獻(xiàn)[25]中兩種f k為常數(shù)的PANS模型計算結(jié)果可以看出,f k＝0.1時,流場f k值已處于較低水平,但計算的回流區(qū)長度要小于本文變f k計算結(jié)果,而f k＝0.3時的計算結(jié)果誤差已較大,這說明了變f k在流場計算中更合理,同時也說明了流場對f k值是比較敏感的。

表2 圓柱繞流平均積分量對比Table 2 Summary of the mean integrated quantities for the cylinder flow

3.5 混合模型調(diào)控機制比較

由前述分析可知,SAS及PANS等渦黏性混合模型均能通過降低流場湍動黏度來求解更多湍流尺度,但SAS及PANS對求解尺度的調(diào)控機制是完全不同的。由前文模型構(gòu)造我們可知,SAS模型中引入了von Karman長度尺度L vk,通過當(dāng)?shù)赝牧鞒叨萀與L vk的比例關(guān)系來控制ω輸運方程中源項QSAS的生成量,因此SAS模型在分離區(qū)是通過L/L vk的調(diào)控作用來增大ω的值,從而降低湍動黏度的值；PANS模型是通過加入一個f k參數(shù)來調(diào)控流場可解/不可解湍流量的比例,f k值越小,可解尺度就釋放的越多,文中PANS模型使用了兩種變f k函數(shù),其核心是當(dāng)?shù)赝牧鞒叨扰c網(wǎng)格尺度的比例關(guān)系,兩者都是通過減小ω方程中耗散項的量,從而使求解的ω值增大,進(jìn)而使湍動黏度減小。

圖16為SAS模型計算的瞬時及時均L/L vk分布?？梢钥吹皆谖槽E分離區(qū)L/L vk的值較大,而在遠(yuǎn)場區(qū)L/L vk基本為0。圖17和圖18是兩種變f k的PANS模型計算得到的瞬時及時均f k分布,可以看到在分離區(qū)f k的值處在較低水平,而在遠(yuǎn)場范圍f k＝1湍流信息均被?；?。通過比較兩類模型中對渦黏度起關(guān)鍵調(diào)控作用的L/L vk與f k分布可以看出,立足于當(dāng)?shù)亓鲃印W(wǎng)格無依賴的SAS模型局部分辨率更高,可以反映出流場結(jié)構(gòu)的無序與豐富狀態(tài)。PANS模型中網(wǎng)格尺度的關(guān)聯(lián)性使得其局部分辨率降低,且計算結(jié)果對網(wǎng)格及計算初值設(shè)置非常敏感。從以上來看,SAS模型對流場渦黏性的調(diào)控機制是有優(yōu)勢的。

圖16 SAS模型計算的瞬時及時均L/L vk分布Fig.16 Time-averaged and instantaneous distributions of L/L vk for SAS

圖17 PANS(des)模型計算得到的瞬時及時均f k分布Fig.17 Time-averaged and instantaneous distributions of f k for PANS(des)

圖18 PANS(tanh)模型計算得到的瞬時及時均f k分布Fig.18 Time-averaged and instantaneous distributions of f k for PANS(tanh)

對比圖17和圖18,顯示了PANS模型的兩種函數(shù)計算的f k分布差別不大,在近尾跡區(qū)f k均處在較低水平,這有效地調(diào)整了流場的求解尺度。同時可以看到PANS(des)模型相比PANS(tanh)模型在近尾跡區(qū)的分辨率更高,其在回流區(qū)后部及遠(yuǎn)尾跡區(qū)的f k值要稍小于后者。分析兩種模型的構(gòu)造可知,PANS(des)模型f k函數(shù)構(gòu)造相對簡單,但其對流場的調(diào)控作用更加直接和靈敏,而PANS(tanh)模型中f k函數(shù)中借助tanh函數(shù),tanh函數(shù)整體從－1～1變化,這使得函數(shù)的整體靈敏度下降,從而影響f k的分布。同時圖13中可以看到PANS(des)的湍動黏度稍高于PANS(tanh),其獲得的氣動力分布也稍好于后者。這也說明了兩種變f k函數(shù)比較,類DES的f k調(diào)控函數(shù)更好一些,計算的湍流耗散也更合理。

4 結(jié) 論

本文以Re＝3900的圓柱繞流作為算例比較了SAS、PANS(des)、PANS(tanh)湍流模型在大分離流動中的計算能力,主要從流場特征捕捉、氣動力特性計算、流場渦黏性控制等方面對所用湍流模型進(jìn)行了比較,并通過不同的網(wǎng)格計算條件對比了不同模型調(diào)控函數(shù)的網(wǎng)格敏感性。計算結(jié)果表明:三種混合模型均具有求解三維小尺度渦運動的能力,并能較完整地展現(xiàn)出亞臨界雷諾數(shù)下圓柱繞流層流剪切層形成、發(fā)展、失穩(wěn)、分離,在下游近尾跡區(qū)形成回流區(qū)直至遠(yuǎn)尾跡區(qū)形成剪切分離區(qū)的過程。同時,通過對比計算結(jié)果得到以下結(jié)論:

1)在對湍流結(jié)構(gòu)的顯示方面,PANS模型相比SAS模型能捕捉到更豐富的湍流結(jié)構(gòu),兩類模型對流場氣動力的計算均存在一定誤差,回流區(qū)計算誤差要大于遠(yuǎn)尾跡區(qū),三種模型氣動力計算與實驗吻合關(guān)系為SAS＞PANS(des)＞PANS(tanh)。

2)SAS模型L vk尺度立足當(dāng)?shù)亓鲃?局部流場分辨度更高,調(diào)控作用更強；PANS模型對網(wǎng)格設(shè)置較敏感。兩類模型相比,SAS模型的流場信息調(diào)控模式更優(yōu)秀,獲得的流動渦黏性更合理。

3)兩種PANS模型中,類DES可變f k分布函數(shù)構(gòu)造相對簡單,所得f k分布更準(zhǔn)確；借助tanh函數(shù)的f k分布函數(shù)調(diào)節(jié)靈敏度較低,計算得到的尾跡區(qū)f k值偏低,對流場調(diào)控作用稍差。