劉 躍管小榮徐 誠
(南京理工大學(xué) 機(jī)械工程學(xué)院,南京 210094)
廣泛存在于航空航天、風(fēng)工程等空氣動(dòng)力學(xué)領(lǐng)域的鈍體繞流問題通常為無序、非線性、多尺度的湍流分離運(yùn)動(dòng),在使用數(shù)值模擬方法(CFD)對(duì)其進(jìn)行分析計(jì)算時(shí),湍流的求解成為影響CFD計(jì)算精度的瓶頸問題[1]。在諸多的湍流模型中,兼顧計(jì)算精度與資源消耗的RANS/LES混合方法受到人們的青睞[2],混合模型在近壁區(qū)及遠(yuǎn)場(chǎng)使用RANS模型,而在分離區(qū)及尾跡區(qū)切換為更合理的類LES濾波模型,極大地提高了計(jì)算精度,且資源消耗遠(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模型為能夠解析寬頻非定常湍流脈動(dòng)而對(duì)計(jì)算網(wǎng)格依賴較小的新一代非定常湍流預(yù)測(cè)方法(URANS),通過加入?yún)⒄债?dāng)?shù)亓鲃?dòng)的馮卡門(von Karman)長(zhǎng)度尺度L vk,合理地釋放更多流場(chǎng)信息[14-15];PANS模型同樣是最近發(fā)展的一種RANS-DNS橋接模型,模型構(gòu)造相對(duì)簡(jiǎn)單,通過加入?;膭?dòng)能比例參數(shù)f k調(diào)整湍流的?;縼碚{(diào)整湍流求解尺度[16];DES模型的提出者Spalart在2000年時(shí)曾預(yù)測(cè),2045年左右計(jì)算機(jī)的計(jì)算能力才能滿足利用LES解決90%的湍流問題,因此,在現(xiàn)階段,驗(yàn)證與發(fā)展RANS/LES方法仍是非常必要的[17-18]。
鈍體繞流中的圓柱繞流計(jì)算模型簡(jiǎn)單,同時(shí)又包含如湍流壁面、剪切層失穩(wěn)、流動(dòng)分離、周期性渦脫落等豐富的流動(dòng)現(xiàn)象,是考察湍流模型計(jì)算性能的標(biāo)準(zhǔn)算例。Re=3900的圓柱繞流是典型的亞臨界雷諾數(shù)流動(dòng),流動(dòng)結(jié)構(gòu)復(fù)雜且耗費(fèi)資源較少,人們對(duì)該雷諾數(shù)的圓柱繞流進(jìn)行了大量的實(shí)驗(yàn)[19-22]及數(shù)值研究[23-29],有可靠豐富的數(shù)據(jù)可供參考。本文選用SAS模型及兩種變f k構(gòu)造函數(shù)的PANS模型計(jì)算了Re=3900的圓柱繞流流場(chǎng),并重點(diǎn)從流場(chǎng)特征捕捉、氣動(dòng)力特性計(jì)算、渦黏性控制等方面對(duì)SAS模型及PANS模型進(jìn)行了比較,通過對(duì)比標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格后的計(jì)算結(jié)果,比較了兩類混合方法對(duì)網(wǎng)格的適應(yīng)性,并探討了兩類混合模型的尺度函數(shù)對(duì)流場(chǎng)的調(diào)控機(jī)制。
控制方程為三維N-S方程,在一般計(jì)算坐標(biāo)系(τ,ξ,η,ζ)下的 守恒形 式為:
式中:為守恒變量為無粘通量為黏性通量,各物理量的具體表述參見文獻(xiàn)[1]。
1.2.1 SST-k-ω模型
SST模型由標(biāo)準(zhǔn)k-ω模型與k-ε模型混合修改而成,該模型兼顧了k-ω模型的近壁性能及k-ε模型的遠(yuǎn)場(chǎng)準(zhǔn)確性,SST模型中湍動(dòng)能k和比耗散率ω的輸運(yùn)方程可表示如下:
式中模型參數(shù)σk、α、β、σω、σω2均由對(duì)應(yīng)的k-ω模型參數(shù)φ1和k-ε模型參數(shù)φ2通過一定函數(shù)計(jì)算得到,其對(duì)應(yīng)的函數(shù)關(guān)系為:
SST模型中湍動(dòng)黏度μ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ù)亓鲃?dòng)拓?fù)涞膙on Karman長(zhǎng)度尺度L vk來作為湍流模型的第二長(zhǎng)度尺度[6]。基于SST模型的SAS模型構(gòu)造方法為在式(3)SST模型的ω方程中加入發(fā)揮自適應(yīng)作用的源項(xiàng)QSAS,其定義式為:
由式(6)可以看出:在SAS模型中出現(xiàn)了由一階及二階速度導(dǎo)數(shù)決定的長(zhǎng)度尺度項(xiàng)L vk,該長(zhǎng)度尺度基于當(dāng)?shù)亓鲃?dòng)而與網(wǎng)格尺度無關(guān),在湍流邊界層內(nèi),L vk能?;瘧T性子區(qū)所有湍流脈動(dòng),同時(shí),在非穩(wěn)態(tài)區(qū)域L vk能根據(jù)當(dāng)?shù)鼐W(wǎng)格分辨的湍流渦動(dòng)態(tài)調(diào)整RANS長(zhǎng)度尺度。在SAS模型中,湍動(dòng)黏度μt仍按照式(5)計(jì)算。
1.2.3 SST-PANS模型
本文中所用PANS模型同樣是基于SST模型修改而成[8],其原理是在SST模型中引入兩個(gè)模型控制參數(shù)f k和fω:
式中:f k為模化湍動(dòng)能與總的湍動(dòng)能之比,fω是?;群纳⒙逝c總的比耗散率的比值。如果假定f k和fω是不隨時(shí)間和空間變化的常數(shù),則k和ω的輸運(yùn)方程(2)、(3)變?yōu)?
將式(8)、(9)與式(2)、(3)對(duì)比可以發(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),改變模化的湍動(dòng)能與比耗散率的比值,從而實(shí)現(xiàn)RANS-DNS的轉(zhuǎn)化。這時(shí)湍動(dòng)黏度μtu變?yōu)?
在早期PANS模型的應(yīng)用中,f k通常被設(shè)置為空間常數(shù),通過計(jì)算發(fā)現(xiàn)對(duì)于整個(gè)流場(chǎng)很難確定合適統(tǒng)一的f k值。對(duì)于復(fù)雜的分離流動(dòng)問題,計(jì)算時(shí)分離區(qū)與尾跡區(qū)需要較小的f k值(f k<1)釋放湍流信息,而近壁區(qū)與遠(yuǎn)場(chǎng)區(qū)采用RANS模型(f k=1)更合理,因此近年來人們逐漸嘗試用適當(dāng)?shù)暮瘮?shù)來計(jì)算f k分布,本文選取其中兩種相對(duì)簡(jiǎn)單、計(jì)算相對(duì)穩(wěn)定的變f k函數(shù)來進(jìn)行對(duì)比。
(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)行計(jì)算時(shí),在分離區(qū)與尾跡區(qū),網(wǎng)格尺度要小于湍流尺度,因此f k<1;在壁面及遠(yuǎn)場(chǎng)區(qū)域,網(wǎng)格尺度大于湍流尺度,因此f k=1。文中將使用類DES型f k分布函數(shù)的PANS模型稱為PANS(des)模型。
(b)Tanh型f k函數(shù)[33]。在該函數(shù)中,f k可表示為:
對(duì)比式(12)與式(13)可知,Tanh型f k函數(shù)由湍流尺度、網(wǎng)格尺度及tanh函數(shù)修改而成,tanh函數(shù)的取值范圍為(-1,1),這樣根據(jù)流動(dòng)與網(wǎng)格關(guān)系得到不同位置的f k取值范圍為(0,1)。文中將使用Tanh型f k函數(shù)的PANS模型稱為PANS(tanh)模型。
文中圓柱繞流計(jì)算雷諾數(shù)為Re=3900(基于圓柱直徑D與遠(yuǎn)場(chǎng)來流速度U∞),圖1為計(jì)算域及邊界條件示意圖。為減小上游來流對(duì)計(jì)算的影響并使尾跡區(qū)流動(dòng)充分發(fā)展,上游來流區(qū)域(x方向)及橫向區(qū)域(y方向)延展至20D,下游計(jì)算區(qū)域設(shè)置為30D。根據(jù)近尾跡區(qū)流向結(jié)構(gòu)的展向波長(zhǎng)設(shè)置計(jì)算域展向長(zhǎng)度(z方向)為πD[2]。圓柱體設(shè)置為無滑移壁面,上游入口流動(dòng)各物理量由無擾動(dòng)自由來流條件給定,下游出口設(shè)置為無反射邊界,計(jì)算域的上下界面設(shè)置為自由滑移壁面,計(jì)算域的前后邊界設(shè)置為周期邊界。文中坐標(biāo)系設(shè)置如圖1右下角所示,軸向從上游到下游方向?yàn)閤正方向,坐標(biāo)系符合右手螺旋法則,橫向?yàn)閥方向,展向?yàn)閦方向,方位角從上游圓柱前駐點(diǎn)為0°開始,順時(shí)針旋轉(zhuǎn)為壁面方位角θ的正方向,下游背風(fēng)點(diǎn)為180°。
圖1 計(jì)算域及邊界條件Fig.1 Computational domain and boundary conditions
圖2為圓柱繞流的局部網(wǎng)格示意,計(jì)算網(wǎng)格為分塊結(jié)構(gòu)化網(wǎng)格。為了保證計(jì)算精度,在壁面附近及尾跡分離區(qū)域?qū)W(wǎng)格進(jìn)行適當(dāng)加密,壁面第一層網(wǎng)格高度為0.002D,保證壁面y+<1。為了探討不同模型的尺度函數(shù)對(duì)網(wǎng)格的敏感性,文中設(shè)置了兩套計(jì)算網(wǎng)格,其中標(biāo)準(zhǔn)網(wǎng)格保證近壁增長(zhǎng)率為1.06,分離區(qū)網(wǎng)格固定為0.1D;加密網(wǎng)格為保證近壁增長(zhǎ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 計(jì)算網(wǎng)格Table 1 Computational grids
文中對(duì)所求方程組空間離散采用有限體積法(FVM),對(duì)流項(xiàng)中速度項(xiàng)的離散采用三階限制(TVD)中心差分格式,湍流項(xiàng)離散采用二階限制中心差分格式,梯度項(xiàng)及黏性項(xiàng)的離散采用四階中心差分格式,壓力與速度耦合選用PISO算法。為了保證計(jì)算的收斂性及穩(wěn)定性,時(shí)間項(xiàng)離散采用二階隱式后退歐拉(Backward)格式,并采用自適應(yīng)時(shí)間步長(zhǎng)方法[24],時(shí)間步長(zhǎng)設(shè)置為庫朗數(shù)(Courant number)小于0.5,并且保證在一個(gè)渦脫落周期內(nèi)包含近500個(gè)時(shí)間步。
文中分別計(jì)算了SST、SAS、PANS(des)、PANS(tanh)模型的圓柱繞流流場(chǎng)。周期性渦脫落是該計(jì)算雷諾數(shù)下的顯著特征。為了獲得合理的統(tǒng)計(jì)數(shù)據(jù),根據(jù)Franke等[34]的研究結(jié)果,圓柱繞流的統(tǒng)計(jì)周期應(yīng)大于40個(gè)脫落周期。因此參照文獻(xiàn)[34,25],本文中總的計(jì)算時(shí)間約為60個(gè)渦脫落周期,在對(duì)非定常結(jié)果的統(tǒng)計(jì)時(shí)去掉前10個(gè)脫落周期,將后50個(gè)周期內(nèi)的統(tǒng)計(jì)數(shù)據(jù)作為有效數(shù)據(jù)進(jìn)行模型的對(duì)比。計(jì)算結(jié)果進(jìn)行時(shí)間與空間平均。圖3中給出了以標(biāo)準(zhǔn)網(wǎng)格計(jì)算結(jié)果為例的圓柱繞流升力系數(shù)(C l)、阻力系數(shù)(C d)的時(shí)程(U∞/D無量綱化)變化曲線??梢钥闯鯮e=3900時(shí)的升力幅值已不再穩(wěn)定,但仍呈現(xiàn)出周期性變化趨勢(shì),同時(shí)阻力系數(shù)表現(xiàn)出強(qiáng)烈的脈動(dòng)現(xiàn)象,并在一平均值附近上下擺動(dòng)。
圖3 升力、阻力系數(shù)隨時(shí)間變化及升力頻譜Fig.3 Time history of lift C l,drag C d coef ficients
圖4中給出了由升力系數(shù)時(shí)程曲線經(jīng)過傅里葉變換(FFT)獲得的渦脫落頻率功率譜,由此可以計(jì)算渦脫落頻率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)格計(jì)算得到的圓柱繞流流場(chǎng)二維瞬時(shí)(t=300)渦量分布及用軸向速度染色的三維Q分布,對(duì)比可知:由于RANS模型湍流黏性較大,抹去了很多重要的流場(chǎng)細(xì)節(jié),計(jì)算結(jié)果只捕捉到大尺度的脫落渦結(jié)構(gòu),Q圖也只先顯示出近似二維的流動(dòng)結(jié)構(gòu),這說明RANS模型基本上沒有求解小尺度流動(dòng)的能力;SAS模型及PANS模型均捕捉到豐富的流動(dòng)結(jié)構(gòu),二維渦量及三維Q值均顯示出復(fù)雜的周期性渦脫落現(xiàn)象。分析三者的Q圖可清晰地發(fā)現(xiàn)附著于壁面的二維層流剪切層及分離區(qū)流向與展向相嵌套的三維無序渦體結(jié)構(gòu)。SAS模型與PANS模型均允許形成湍流譜,表現(xiàn)出預(yù)期的類似LES模型的求解小尺度運(yùn)動(dòng)的能力。同時(shí),仔細(xì)對(duì)比SAS模型與PANS模型的結(jié)果可以看出,PANS模型能捕捉到比SAS模型更豐富的湍流結(jié)構(gòu)。
圖5 不同模型計(jì)算得到的流場(chǎng)瞬時(shí)渦量與Q值圖(t=300)Fig.5 Instantaneous vorticity and Q of different models(t=300)
由實(shí)驗(yàn)及以往計(jì)算結(jié)果[23]可知,圓柱繞流壁面剪切層的發(fā)展對(duì)外界影響因素非常敏感,微小外部干擾很容易造成邊界層的提前分離,使得剪切層變短,進(jìn)而影響回流區(qū)及整個(gè)流場(chǎng)的測(cè)量或計(jì)算結(jié)果。圖6為用渦量等值線表示的標(biāo)準(zhǔn)網(wǎng)格條件下四種模型的瞬時(shí)(t=300)剪切層,通過對(duì)比可知剪切層長(zhǎng)度關(guān)系為:SAS>PANS(des)>PANS(tanh)>SST。
圖6 不同模型瞬時(shí)分離剪切層示意圖(t=300)Fig.6 Instantaneous separating shear layers of different models(t=300)
圖7為以SAS模型為例的y=0平面軸向速度時(shí)均圖,可以清晰地看出在圓柱近尾跡區(qū)域形成一個(gè)流速與來流相反的區(qū)域,即圓柱繞流的回流區(qū)。圖8為圓柱下游尾跡區(qū)中心線上時(shí)均軸向速度分布。由圖可知,圓柱背風(fēng)點(diǎn)處速度為0,速度為逐漸增大的負(fù)值,到達(dá)回流區(qū)中心時(shí)速度負(fù)值最大。背風(fēng)點(diǎn)與空間流動(dòng)的再附點(diǎn)即流向速度恢復(fù)到0的點(diǎn)之間的距離即為回流區(qū)長(zhǎng)度(L r/D)。四種模型均模擬出了圓柱后的速度分區(qū),但彼此間差別較大,特別是對(duì)回流區(qū)的模擬。圖8(b)圖是對(duì)回流區(qū)中心線時(shí)均軸向速度的放大圖。首先需要說明的是用空心圓圈表示的實(shí)驗(yàn)數(shù)據(jù)是Lourenco&Shih(1993)用PIV測(cè)量獲得[20],用實(shí)心圓表示的實(shí)驗(yàn)數(shù)據(jù)是Parnaudeau等(2008)同樣用PIV測(cè)量獲得[21],可以看到兩者實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度差別較大。最近幾年的文獻(xiàn)也出現(xiàn)了與Parnaudeau實(shí)驗(yàn)接近或與Lourenco實(shí)驗(yàn)接近的兩種數(shù)值計(jì)算結(jié)果,有意思的是Kravchenko等(2000)在用B-Spline伽遼金方法對(duì)Re=3900進(jìn)行的LES數(shù)值計(jì)算中[23],最先發(fā)現(xiàn)了數(shù)值結(jié)果與Lourenco實(shí)驗(yàn)數(shù)據(jù)在回流區(qū)與近尾跡區(qū)速度型的差別。Kravchenko計(jì)算的回流區(qū)長(zhǎng)度大于Lourenco實(shí)驗(yàn)值,其當(dāng)時(shí)推測(cè)可能有外加擾動(dòng)導(dǎo)致Lourenco實(shí)驗(yàn)中圓柱壁面剪切層過早分離,造成剪切層變短,進(jìn)而影響回流區(qū)長(zhǎng)度。Parnaudeau實(shí)驗(yàn)數(shù)據(jù)在回流區(qū)及速度型方面均驗(yàn)證了Kravchenko的數(shù)值計(jì)算結(jié)果。由于實(shí)驗(yàn)時(shí)間較晚,實(shí)驗(yàn)方法更先進(jìn)、嚴(yán)謹(jǐn),本文中我們認(rèn)為Parnaudeau的實(shí)驗(yàn)值是更可信的。分析圖8(b)可知,SST模型計(jì)算的圓柱回流區(qū)長(zhǎng)度誤差較大,SAS與PANS計(jì)算的回流區(qū)長(zhǎng)度關(guān)系為SAS>PANS(des)>PANS(tanh)。三種模型計(jì)算的回流區(qū)區(qū)域位于兩個(gè)實(shí)驗(yàn)值之間,與Parnaudeau實(shí)驗(yàn)值仍存在一定誤差,其中SAS模型與Parnaudeau實(shí)驗(yàn)最接近。以上結(jié)論與圖6中得出的剪切層長(zhǎng)度結(jié)論也是一致的。
圖7 y=0平面軸向速度時(shí)均分布Fig.7 Mean streamwise velocity along the y=0 surface
圖8 圓柱尾跡中心線上時(shí)均軸向速度分布Fig.8 Mean streamwise velocity along the wake center line
圖9~圖11為圓柱繞流下游不同截面處的時(shí)均軸向速度分布(圖例見圖8)??梢钥吹捷S向速度流型從近尾跡區(qū)到遠(yuǎn)尾跡區(qū)由U型向V型轉(zhuǎn)變,結(jié)合前部分析可知,剪切層過早分離同樣會(huì)導(dǎo)致速度型的過早轉(zhuǎn)變,圖中兩種實(shí)驗(yàn)值的對(duì)比也說明了這一點(diǎn)。圖9近壁x/D=0.58時(shí),四種模型的速度型均為U型,彼此差別不大,而到回流區(qū)x/D=1.06時(shí),四種模型出現(xiàn)差別,SST模型為完全的V型,SAS模型與兩種PANS模型同樣沒有計(jì)算出更加正確的U型,而是介于U型與V型之間,其中SAS模型與Parnaudeau實(shí)驗(yàn)得到的U型更接近,同時(shí)PANS(des)要稍好于PANS(tanh);在近尾跡回流區(qū)后部x/D=1.54、2.02位置處,除了SST沒有計(jì)算出正確的速度型,SAS模型與兩種PANS模型的計(jì)算結(jié)果只是速度峰值與實(shí)驗(yàn)值存在一定誤差,整體均與實(shí)驗(yàn)值吻合較好;在遠(yuǎn)尾跡區(qū)x/D=4、7、10位置處,四種模型計(jì)算的速度型均較好,速度峰值與實(shí)驗(yàn)存在較小誤差。以上分析說明近尾跡區(qū)域特別是回流區(qū)流動(dòng)較為復(fù)雜、敏感,計(jì)算控制更困難,同時(shí)遠(yuǎn)尾跡區(qū)受近壁區(qū)影響較小,使得遠(yuǎn)尾跡區(qū)計(jì)算誤差小于近尾跡區(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ū)附近不同截面處的時(shí)均橫向速度分布(圖例見圖8)。由x/D=1.06處的時(shí)均橫向速度與實(shí)驗(yàn)值可以看到,SAS與兩種PANS模型得到的速度型均與Parnaudeau實(shí)驗(yàn)值較吻合,與SST模型比較顯示出較高的計(jì)算精度。x/D=1.54、2.02處的對(duì)比結(jié)果可以看出四種模型在速度型與峰值上均存在一定誤差,與軸向速度計(jì)算結(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)
根據(jù)混合模型的構(gòu)造原理可知,其湍流尺度求解能力是基于對(duì)渦黏性的控制作用。圖13為三種混合模型計(jì)算的圓柱繞流時(shí)均湍動(dòng)黏度μt/μ(湍動(dòng)黏度與來流黏度之比)分布,上部為標(biāo)準(zhǔn)網(wǎng)格計(jì)算結(jié)果,下部為加密網(wǎng)格計(jì)算結(jié)果。僅從數(shù)值來看,SAS模型與兩種PANS模型的湍動(dòng)黏度分布關(guān)系為SAS>PANS(des)>PANS(tanh)。PANS模型所求湍動(dòng)黏度要小于SAS模型,這也解釋了圖4中PANS模型對(duì)流場(chǎng)的捕捉要好于SAS模型。但需要注意的是,較低的湍動(dòng)黏度不一定獲得較準(zhǔn)確的數(shù)值計(jì)算結(jié)果。由文中氣動(dòng)力的比較可知,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模型對(duì)網(wǎng)格敏感的缺點(diǎn),這時(shí)近壁RANS區(qū)域被破壞,而網(wǎng)格尺度還沒有小到能夠進(jìn)行大渦模擬的程度。
圖13 時(shí)均湍流黏度比(μ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)格條件下計(jì)算的局部(x/D=1.06)時(shí)均軸向雷諾應(yīng)力分布(圖例見圖15)??梢钥吹絊AS模型及兩種PANS模型計(jì)算的雷諾應(yīng)力與實(shí)驗(yàn)值均有較大差距,這可能是由于渦黏性模型的內(nèi)在缺陷造成的。我們?cè)诖岁P(guān)注不同的模型在不同網(wǎng)格條件下的計(jì)算差異??梢钥吹絊AS模型在標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格條件下均高估了雷諾應(yīng)力值,但兩種網(wǎng)格下的計(jì)算結(jié)果差別不大;兩種PANS模型對(duì)不同的網(wǎng)格表現(xiàn)出相似的規(guī)律,在加密網(wǎng)格時(shí)雷諾應(yīng)力均出現(xiàn)了明顯的降低,結(jié)合圖13中的黏度分布可知,由于近壁區(qū)域網(wǎng)格的加密,湍流黏度降低,導(dǎo)致?;睦字Z應(yīng)力減少,同時(shí)有沒有求解足夠的雷諾應(yīng)力,導(dǎo)致了雷諾應(yīng)力不足現(xiàn)象的發(fā)生,這也驗(yàn)證了PANS構(gòu)造函數(shù)對(duì)網(wǎng)格較敏感的缺陷。同時(shí)兩種PANS模型函數(shù)比較可以看出,PANS(des)模型計(jì)算的雷諾應(yīng)力及氣動(dòng)力分布要稍優(yōu)于PANS(tanh)模型,這也說明了PANS(tanh)模型構(gòu)造函數(shù)對(duì)網(wǎng)格的敏感程度要大于PANS(des)模型的構(gòu)造函數(shù)。
圖14 局部時(shí)均軸向雷諾應(yīng)力分布Fig.14 Time-averaged distributions of streamwise Reynolds stress
圖15為三種模型在不同網(wǎng)格條件下的時(shí)均軸向速度分布對(duì)比。結(jié)合前文分析可知,SAS模型網(wǎng)格獨(dú)立性較好,標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格下的計(jì)算結(jié)果沒有體現(xiàn)出太大差別,且由于對(duì)剪切層的計(jì)算更準(zhǔn)確,計(jì)算得到的回流區(qū)大小更接近實(shí)驗(yàn)值。兩種PANS模型由于對(duì)網(wǎng)格較敏感,在加密網(wǎng)格時(shí)延緩了剪切層的轉(zhuǎn)捩,導(dǎo)致回流區(qū)變長(zhǎng)。這在PANS(tanh)模型的計(jì)算結(jié)果中表現(xiàn)更明顯。
圖15 時(shí)均軸向速度分布Fig.15 Time-averaged distributions of streamwise velocity
表2為四種模型計(jì)算的圓柱繞流平均積分統(tǒng)計(jì)量與相關(guān)實(shí)驗(yàn)值及數(shù)值結(jié)果對(duì)比??梢钥吹皆谂c剪切層運(yùn)動(dòng)密切相關(guān)的參數(shù)中,普通網(wǎng)格下SAS與兩種PANS模型均取得了與實(shí)驗(yàn)值較吻合的回流區(qū)長(zhǎng)度(L r/D)計(jì)算結(jié)果,與構(gòu)造相對(duì)完善的LES-THE、SA-IDDES、V2-f-DES模型的計(jì)算結(jié)果也是相近的,其中SAS模型的計(jì)算結(jié)果與實(shí)驗(yàn)值最接近;加密網(wǎng)格時(shí),兩種PANS模型均過大地預(yù)測(cè)了回流區(qū)長(zhǎng)度。在阻力系數(shù)(C d)及分離角(θSep)的計(jì)算中,在標(biāo)準(zhǔn)網(wǎng)格下三種混合模型計(jì)算誤差均較小,C d滿足實(shí)驗(yàn)值的0.98±0.05區(qū)間,在加密網(wǎng)格時(shí),SAS模型的計(jì)算結(jié)果變化不大,兩種PANS模型獲得了較小的分離角及阻力系數(shù)。在其余的比較項(xiàng)中,可以看到四種模型對(duì)圓柱尾跡旋渦脫落頻率的計(jì)算及壁面壓力的計(jì)算均較準(zhǔn)確,所計(jì)算的斯特勞哈數(shù)(S t)與背風(fēng)點(diǎn)壓力系數(shù)(-C pb)差別不大。對(duì)于SST模型,由于其構(gòu)造?;怂械耐牧髁?導(dǎo)致重要的流場(chǎng)細(xì)節(jié)被抹平,其計(jì)算的阻力系數(shù)、背風(fēng)點(diǎn)壓力系數(shù)、回流區(qū)長(zhǎng)度、分離角均存在較大誤差。同時(shí)可以看到兩種網(wǎng)格條件下PANS(des)計(jì)算結(jié)果稍好于PANS(tanh)。比較文獻(xiàn)[25]中兩種f k為常數(shù)的PANS模型計(jì)算結(jié)果可以看出,f k=0.1時(shí),流場(chǎng)f k值已處于較低水平,但計(jì)算的回流區(qū)長(zhǎng)度要小于本文變f k計(jì)算結(jié)果,而f k=0.3時(shí)的計(jì)算結(jié)果誤差已較大,這說明了變f k在流場(chǎng)計(jì)算中更合理,同時(shí)也說明了流場(chǎng)對(duì)f k值是比較敏感的。
表2 圓柱繞流平均積分量對(duì)比Table 2 Summary of the mean integrated quantities for the cylinder flow
由前述分析可知,SAS及PANS等渦黏性混合模型均能通過降低流場(chǎng)湍動(dòng)黏度來求解更多湍流尺度,但SAS及PANS對(duì)求解尺度的調(diào)控機(jī)制是完全不同的。由前文模型構(gòu)造我們可知,SAS模型中引入了von Karman長(zhǎng)度尺度L vk,通過當(dāng)?shù)赝牧鞒叨萀與L vk的比例關(guān)系來控制ω輸運(yùn)方程中源項(xiàng)QSAS的生成量,因此SAS模型在分離區(qū)是通過L/L vk的調(diào)控作用來增大ω的值,從而降低湍動(dòng)黏度的值;PANS模型是通過加入一個(gè)f k參數(shù)來調(diào)控流場(chǎng)可解/不可解湍流量的比例,f k值越小,可解尺度就釋放的越多,文中PANS模型使用了兩種變f k函數(shù),其核心是當(dāng)?shù)赝牧鞒叨扰c網(wǎng)格尺度的比例關(guān)系,兩者都是通過減小ω方程中耗散項(xiàng)的量,從而使求解的ω值增大,進(jìn)而使湍動(dòng)黏度減小。
圖16為SAS模型計(jì)算的瞬時(shí)及時(shí)均L/L vk分布??梢钥吹皆谖槽E分離區(qū)L/L vk的值較大,而在遠(yuǎn)場(chǎng)區(qū)L/L vk基本為0。圖17和圖18是兩種變f k的PANS模型計(jì)算得到的瞬時(shí)及時(shí)均f k分布,可以看到在分離區(qū)f k的值處在較低水平,而在遠(yuǎn)場(chǎng)范圍f k=1湍流信息均被模化。通過比較兩類模型中對(duì)渦黏度起關(guān)鍵調(diào)控作用的L/L vk與f k分布可以看出,立足于當(dāng)?shù)亓鲃?dòng)、對(duì)網(wǎng)格無依賴的SAS模型局部分辨率更高,可以反映出流場(chǎng)結(jié)構(gòu)的無序與豐富狀態(tài)。PANS模型中網(wǎng)格尺度的關(guān)聯(lián)性使得其局部分辨率降低,且計(jì)算結(jié)果對(duì)網(wǎng)格及計(jì)算初值設(shè)置非常敏感。從以上來看,SAS模型對(duì)流場(chǎng)渦黏性的調(diào)控機(jī)制是有優(yōu)勢(shì)的。
圖16 SAS模型計(jì)算的瞬時(shí)及時(shí)均L/L vk分布Fig.16 Time-averaged and instantaneous distributions of L/L vk for SAS
圖17 PANS(des)模型計(jì)算得到的瞬時(shí)及時(shí)均f k分布Fig.17 Time-averaged and instantaneous distributions of f k for PANS(des)
圖18 PANS(tanh)模型計(jì)算得到的瞬時(shí)及時(shí)均f k分布Fig.18 Time-averaged and instantaneous distributions of f k for PANS(tanh)
對(duì)比圖17和圖18,顯示了PANS模型的兩種函數(shù)計(jì)算的f k分布差別不大,在近尾跡區(qū)f k均處在較低水平,這有效地調(diào)整了流場(chǎng)的求解尺度。同時(shí)可以看到PANS(des)模型相比PANS(tanh)模型在近尾跡區(qū)的分辨率更高,其在回流區(qū)后部及遠(yuǎn)尾跡區(qū)的f k值要稍小于后者。分析兩種模型的構(gòu)造可知,PANS(des)模型f k函數(shù)構(gòu)造相對(duì)簡(jiǎn)單,但其對(duì)流場(chǎng)的調(diào)控作用更加直接和靈敏,而PANS(tanh)模型中f k函數(shù)中借助tanh函數(shù),tanh函數(shù)整體從-1~1變化,這使得函數(shù)的整體靈敏度下降,從而影響f k的分布。同時(shí)圖13中可以看到PANS(des)的湍動(dòng)黏度稍高于PANS(tanh),其獲得的氣動(dòng)力分布也稍好于后者。這也說明了兩種變f k函數(shù)比較,類DES的f k調(diào)控函數(shù)更好一些,計(jì)算的湍流耗散也更合理。
本文以Re=3900的圓柱繞流作為算例比較了SAS、PANS(des)、PANS(tanh)湍流模型在大分離流動(dòng)中的計(jì)算能力,主要從流場(chǎng)特征捕捉、氣動(dòng)力特性計(jì)算、流場(chǎng)渦黏性控制等方面對(duì)所用湍流模型進(jìn)行了比較,并通過不同的網(wǎng)格計(jì)算條件對(duì)比了不同模型調(diào)控函數(shù)的網(wǎng)格敏感性。計(jì)算結(jié)果表明:三種混合模型均具有求解三維小尺度渦運(yùn)動(dòng)的能力,并能較完整地展現(xiàn)出亞臨界雷諾數(shù)下圓柱繞流層流剪切層形成、發(fā)展、失穩(wěn)、分離,在下游近尾跡區(qū)形成回流區(qū)直至遠(yuǎn)尾跡區(qū)形成剪切分離區(qū)的過程。同時(shí),通過對(duì)比計(jì)算結(jié)果得到以下結(jié)論:
1)在對(duì)湍流結(jié)構(gòu)的顯示方面,PANS模型相比SAS模型能捕捉到更豐富的湍流結(jié)構(gòu),兩類模型對(duì)流場(chǎng)氣動(dòng)力的計(jì)算均存在一定誤差,回流區(qū)計(jì)算誤差要大于遠(yuǎn)尾跡區(qū),三種模型氣動(dòng)力計(jì)算與實(shí)驗(yàn)吻合關(guān)系為SAS>PANS(des)>PANS(tanh)。
2)SAS模型L vk尺度立足當(dāng)?shù)亓鲃?dòng),局部流場(chǎng)分辨度更高,調(diào)控作用更強(qiáng);PANS模型對(duì)網(wǎng)格設(shè)置較敏感。兩類模型相比,SAS模型的流場(chǎng)信息調(diào)控模式更優(yōu)秀,獲得的流動(dòng)渦黏性更合理。
3)兩種PANS模型中,類DES可變f k分布函數(shù)構(gòu)造相對(duì)簡(jiǎn)單,所得f k分布更準(zhǔn)確;借助tanh函數(shù)的f k分布函數(shù)調(diào)節(jié)靈敏度較低,計(jì)算得到的尾跡區(qū)f k值偏低,對(duì)流場(chǎng)調(diào)控作用稍差。