PANS方法在雙圓柱繞流數(shù)值模擬中的性能分析

杜若凡，閻超，羅大海

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191)

湍流問題由于其自身的復(fù)雜特性，已經(jīng)成為制約計(jì)算流體動力學(xué)(Computational Fluid Dynamics，CFD)精確計(jì)算的瓶頸問題，極大地影響著CFD 預(yù)測的精度［1］.在當(dāng)前條件下，RANS/LES 混合方法由于其經(jīng)濟(jì)性與有效性，相較于其他方法有較為明顯的優(yōu)勢.近年來，研究者們發(fā)展了種類繁多的混合方法，包括DES(Detached-Eddy Simulation)方法［2］、SAS(Scale Adaptive Simulation)方法［3］、PANS(Partially Averaged Navier-Stokes)方法［4］以及 LNS(Limited Numerical Scales)方法［5］等.其中，PANS方法基于RANS方程進(jìn)行修改，相比RANS方法對分離流動的模擬更為有效，但求解的尺度又遠(yuǎn)少于 LES求解的尺度［6］，PANS方法在單圓柱繞流［7］和跨聲速空腔流動［8］等數(shù)值模擬中均得到了成功應(yīng)用.

相比于上面提到的單圓柱繞流和跨聲速空腔流動，雙圓柱繞流問題更為復(fù)雜，包含豐富的流動現(xiàn)象，如湍流邊界層、流動分離、自由剪切層失穩(wěn)、渦脫落以及渦與圓柱壁面相互作用等，是考察湍流模擬方法的優(yōu)秀算例.在 AIAA BANC［9］(Benchmark problems for Airframe Noise Computations)會議及歐盟 ATAAC［10］(Advanced Turbulence simulation for Aerodynamic Application Challenges)項(xiàng)目中，雙圓柱繞流均是標(biāo)準(zhǔn)算例之一.該問題可視作對飛機(jī)起落架外形的簡化，而起落架是飛機(jī)起降過程中噪聲的主要來源，故該流動問題的研究對認(rèn)識起落架繞流流場和分析噪聲產(chǎn)生機(jī)理具有重要意義.

美國Langley研究中心對雙圓柱繞流進(jìn)行了多次實(shí)驗(yàn)［11-13］，獲取了豐富可靠的實(shí)驗(yàn)數(shù)據(jù).文獻(xiàn)［9，14］對參與第一屆BANC研討會的眾多組織機(jī)構(gòu)的結(jié)果進(jìn)行了匯總分析，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，總結(jié)了 DDES(Delayed Detached Eddy Simulation)、FSM(Flow Simulation Methodology)、LES(Large Eddy Simulation)等不同湍流模擬方法對雙圓柱繞流的模擬能力，發(fā)現(xiàn)不同方法得到的結(jié)果差異明顯.文獻(xiàn)［15］采用多種RANS/LES混合方法對雙圓柱繞流進(jìn)行模擬，對不同方法進(jìn)行對比分析.文獻(xiàn)［16］采用 DES(Detached-Eddy Simulation)和IDDES(Improved Delayed Detached-Eddy Simulation)等方法借助雙圓柱繞流探討數(shù)值耗散對計(jì)算結(jié)果的影響，取得了一些有意義的成果.在此選取該流動問題作為研究對象，對本文的PANS方法進(jìn)行考察與評估.

1 計(jì)算方法

PANS方法以RANS方程為模板，本文采用的PANS方法基于兩方程SST(Shear Stress Transport)模式，原始 SST 模式方程［17］為

式中:fk為?；膭幽躪u與總的湍動能k的比值;fω為模化比耗散率ωu與總的比耗散率ω的比值.假定fk和fω是不隨時(shí)間和空間變化的常數(shù)，從而SST模式方程轉(zhuǎn)化為

式中:修改后的模型常數(shù)分別為

?；臏u黏性系數(shù)μu和湍動能生成項(xiàng)Pku分別為

式中:ρ為密度;F2為混合函數(shù);S為應(yīng)變率;a1為常數(shù);各參數(shù)詳見文獻(xiàn)［18］.

由此，通過將SST模式中的常數(shù)替換為與fk和fω有關(guān)的函數(shù)，較為便捷地得到了PANS方法模型方程.PANS方法中，通過調(diào)節(jié)fk和fω的值，決定?；耐膭幽?耗散率的比例，進(jìn)而確定求解的湍流尺度.對于高雷諾數(shù)流動，一般取 fω=1/fk，fk可選取不同的值［7］.然而對于復(fù)雜流動問題模擬，尤其是針對高雷諾數(shù)流動的近壁區(qū)，在計(jì)算網(wǎng)格密度達(dá)不到大渦模擬的要求時(shí)，采用RANS處理是更為合理的，此時(shí)fk在近壁區(qū)取值為1，PANS方法退化為RANS，而在其他流動區(qū)域(分離區(qū)或尾跡區(qū))fk取值小于1，以求解部分的湍流尺度.因此，PANS方法中，fk可變的做法是更為合適的.但在方程(4)和方程(5)的求導(dǎo)過程中，可變fk的引入會增加額外的fk對于時(shí)間和空間的導(dǎo)數(shù)項(xiàng)(或稱之為交換誤差)，目前的PANS方法忽略了交換誤差對計(jì)算結(jié)果的影響.

總之，對fk的取值方式可以分為兩種，一種是讓fk在整個(gè)計(jì)算域內(nèi)取一個(gè)定值，本文中稱為統(tǒng)一fk方法;另一種是讓fk在不同區(qū)域取不同的值，本文稱為可變fk方法，本文選取的具體fk分布函數(shù)［18］為

式中:lu為當(dāng)?shù)氐耐牧髁鲃映叨?Δ為3個(gè)方向網(wǎng)格尺度Δx、Δy和Δz的最大值;CPANS為常數(shù)，此處取0.3.可以看出，可變fk方法中fk的取值取決于當(dāng)?shù)赝牧鞒叨群途W(wǎng)格尺度的比值，這一做法與SST DES模型中FDES函數(shù)的構(gòu)造非常相似.近壁區(qū)和遠(yuǎn)場的湍流特征尺度很小，與網(wǎng)格尺度相比為小值，此時(shí)fk分布函數(shù)保證近壁區(qū)和遠(yuǎn)場處fk=1，PANS模型退化為RANS方程;在尾跡區(qū)，湍流特征尺度大于當(dāng)?shù)鼐W(wǎng)格尺度，此時(shí)fk＜1，較小的fk使得PANS方法可以直接求解更多的湍流尺度，模型表現(xiàn)出更低的渦黏性，可以更為準(zhǔn)確地捕捉尾跡區(qū)的流動結(jié)構(gòu).

2 計(jì)算結(jié)果與分析

雙圓柱模型由兩個(gè)直徑相同的圓柱體沿來流方向排列而成，如圖1所示，兩圓柱體中心線之間的距離L=3.7D，D為圓柱直徑，θ為圓柱方位角.坐標(biāo)原點(diǎn)位于前圓柱圓心.自由來流馬赫數(shù)Ma=0.1274，基于圓柱直徑的雷諾數(shù)ReD=1.66×105.

圖1 雙圓柱構(gòu)型的xOy截面示意圖Fig.1 Schematic diagram of xOy plane for tandem cylinders configuration

結(jié)合課題組以往計(jì)算經(jīng)驗(yàn)［19-20］，進(jìn)行網(wǎng)格生成與計(jì)算方法選擇.計(jì)算網(wǎng)格如圖2所示，在文獻(xiàn)［10］提供的網(wǎng)格基礎(chǔ)上進(jìn)行修改得到，計(jì)算域流向×法向×展向?yàn)?0D×20D×3D，網(wǎng)格總量約為625萬，壁面第一層法向網(wǎng)格高度保證y+＜1，且在兩圓柱中間區(qū)域以及后圓柱尾跡區(qū)域網(wǎng)格基本保持各向同性(即3個(gè)方向網(wǎng)格間距相同).對流項(xiàng)離散采用5階WENO的Roe格式，黏性項(xiàng)采用4階中心差分.壁面設(shè)置為無滑移絕熱壁，展向邊界設(shè)為周期性邊界條件.時(shí)間步長取 Δt=5.9×10-5s，所得結(jié)果進(jìn)行時(shí)間與展向平均.

圖2 雙圓柱構(gòu)型計(jì)算網(wǎng)格示意圖Fig.2 Schematic diagram of computational grid for tandem cylinders configuration

首先，結(jié)合瞬時(shí)流場，觀察雙圓柱繞流流動結(jié)構(gòu)，分析各方法所得流場存在的差異，定性比較各方法的不同.

圖3 不同方法瞬時(shí)展向渦量云圖Fig.3 Contours of instantaneous spanwise vortices from different methods

圖3為瞬時(shí)展向渦量云圖，以圖3(a)為例觀察流動結(jié)構(gòu).來流經(jīng)過前圓柱后拖出剪切層，剪切層逐漸發(fā)展，之后迅速失穩(wěn)，二維渦結(jié)構(gòu)逐漸發(fā)展為流向渦與展向渦相嵌套的三維渦結(jié)構(gòu)，并出現(xiàn)復(fù)雜的周期性渦脫落.脫落的渦在兩圓柱間相互作用，拉伸、變形、破碎，并進(jìn)一步撞擊到后圓柱壁面.渦撞擊到后圓柱壁面后流動更加紊亂，與后圓柱剪切層相互作用，在尾跡區(qū)出現(xiàn)復(fù)雜的分離流動.對比圖3中各結(jié)果，對于統(tǒng)一fk方法，隨著fk取值增大，瞬時(shí)展向渦量云圖中小尺度渦結(jié)構(gòu)逐漸減少，且前圓柱剪切層長度存在差異:fk=0.1時(shí)，兩圓柱間流動結(jié)構(gòu)豐富、精細(xì);fk=0.3時(shí)，剪切層長度明顯減小，觀察到的渦結(jié)構(gòu)明顯減少;fk=0.5時(shí)，所得結(jié)果與URANS結(jié)果類似，幾乎看不到任何小尺度流動結(jié)構(gòu)，完全不同于實(shí)驗(yàn)觀測的結(jié)果.這是因?yàn)椋谳^大fk值下，過多的湍流脈動被模化，流場細(xì)節(jié)被嚴(yán)重抹平，與真實(shí)流動情況不符.對于可變fk方法，其求解湍流尺度的豐富程度介于fk=0.1與fk=0.3的統(tǒng)一方法之間.將本文PANS方法結(jié)果與SST DES 方 法 及 圖 4［13］BART PIV(Particle Image Velocimetny)實(shí)驗(yàn)結(jié)果進(jìn)行對比，除fk=0.5的統(tǒng)一方法外，各方法都可以表現(xiàn)出求解小尺度流動結(jié)構(gòu)的能力.

圖5為前、后圓柱壁面平均壓力系數(shù)Cp分布，θ以上游駐點(diǎn)為0°，順時(shí)針方向?yàn)檎鐖D1所示.前、后圓柱壓力系數(shù)分布特點(diǎn)相似，均在90°位置出現(xiàn)壓力峰值，圓柱背風(fēng)面分離產(chǎn)生壓力平臺區(qū)，壓力系數(shù)保持在-0.5左右.圖5～圖7中“EXP.1”來自文獻(xiàn)［12］實(shí)驗(yàn)結(jié)果，“EXP.2”來自文獻(xiàn)［13］實(shí)驗(yàn)結(jié)果.對前圓柱，fk=0.1的統(tǒng)一方法所模擬的壓力峰值(絕對值)最小，分離提前;fk=0.3及fk=0.5時(shí)所得壓力峰值高于實(shí)驗(yàn)值，分離推遲，分析認(rèn)為這是由于其模擬的渦黏性偏高造成的.可變fk方法與實(shí)驗(yàn)結(jié)果及SST DES計(jì)算結(jié)果均吻合較好，差異不大.對后圓柱，不同方法所得結(jié)果的規(guī)律性與前圓柱一致，仍為可變fk方法和SST DES方法與實(shí)驗(yàn)結(jié)果最為接近.其中，可變 fk方法與“EXP.1”結(jié)果更為接近，SST DES方法更接近“EXP.2”的結(jié)果.

圖4 BART PIV實(shí)驗(yàn)結(jié)果［13］Fig.4 Experimental results from BART PIV［13］

圖5 雙圓柱壁面平均壓力系數(shù)分布Fig.5 Mean pressure coefficient distribution of surfaces of both cylinders

圖6給出了不同方法流場中心線上平均流向速度分布的對比，橫坐標(biāo)為x坐標(biāo)與D的比值，縱坐標(biāo)為流向速度u與來流速度u∞的比值.流動經(jīng)過前、后圓柱均出現(xiàn)分離區(qū)，對應(yīng)平均流向速度小于0的區(qū)域.對比各方法，計(jì)算得到的分離區(qū)大小差異明顯.對于統(tǒng)一fk方法，隨著fk取值增大，所得分離區(qū)逐漸變小.可變fk方法與SST DES和實(shí)驗(yàn)結(jié)果均吻合較好，僅存在細(xì)微差異:在圓柱間區(qū)域，可變fk方法所得分離區(qū)略大于實(shí)驗(yàn)結(jié)果和SST DES結(jié)果;在后圓柱尾跡區(qū)，可變fk方法所得分離區(qū)稍小.需要注意的是，圖6(b)中，“EXP.1”為后圓柱有轉(zhuǎn)捩帶裝置的實(shí)驗(yàn)結(jié)果，可以看出轉(zhuǎn)捩帶對流向速度的實(shí)驗(yàn)結(jié)果有明顯影響，可變fk方法和SST DES的計(jì)算結(jié)果與帶有轉(zhuǎn)捩帶的實(shí)驗(yàn)結(jié)果更為吻合.

圖6 中心線流向平均速度分布Fig.6 Mean streamwise velocity distribution at centerline

圖7 中心線二維湍動能分布Fig.7 2D turbulence kinetic energy distribution at centerline

圖8為可變方法計(jì)算的fk分布.在圓柱壁面附近的邊界層區(qū)域及遠(yuǎn)場，fk被限制為1，保證該區(qū)域由RANS主導(dǎo);在圓柱間及尾跡區(qū)等存在豐富渦結(jié)構(gòu)的區(qū)域，fk取值大致在0.15～0.25.需要說明的是，fk在上游圓柱剪切層外由小值快速變化到1，該區(qū)域不由模型黏性主導(dǎo)，此處fk的快速變化對PANS模型黏性影響很小，從PANS模型計(jì)算結(jié)果與實(shí)驗(yàn)及DES結(jié)果的對比看來，fk快速變化導(dǎo)致的交換誤差對PANS模型的正確求解影響不大.與統(tǒng)一fk方法相比，可變fk方法體現(xiàn)出較為明顯的優(yōu)勢.分析認(rèn)為，統(tǒng)一方法中fk難以取到最優(yōu)值，尤其對于復(fù)雜流場，統(tǒng)一的fk取值難以滿足不同區(qū)域的計(jì)算需要;可變fk方法考慮當(dāng)?shù)鼐W(wǎng)格尺度與流動信息，兼顧不同區(qū)域的計(jì)算要求，對近壁區(qū)的流動合理地以RANS處理，同時(shí)通過調(diào)節(jié)fk的分布保證分離的剪切層和圓柱尾跡區(qū)更多的湍流尺度得以求解，適用于復(fù)雜流動模擬.

圖8 可變fkPANS方法計(jì)算中fk的分布Fig.8 Distribution of fkfor variable fkPANS method

3 結(jié)論

1)對于統(tǒng)一fk方法，模型參數(shù)fk對計(jì)算結(jié)果有顯著影響.隨著fk取值增大，PANS方法?；嗟耐牧鞒叨?，展向流動受到抑制，小尺度渦減少，流動結(jié)構(gòu)趨于粗糙，流動更多呈現(xiàn)出二維特征.與實(shí)驗(yàn)結(jié)果相比，fk=0.3、0.5所得結(jié)果流動分離推遲，分離區(qū)偏小，湍流脈動偏弱.而fk=0.1所得結(jié)果分離提前，分離區(qū)偏大，脈動過強(qiáng)，與實(shí)驗(yàn)同樣存在明顯差異.故對于雙圓柱繞流，統(tǒng)一fk方法適應(yīng)性不足，難以處理該復(fù)雜流動問題，無法給出令人滿意的結(jié)果.

2)可變fk方法將模型參數(shù)fk與當(dāng)?shù)亓鲃有畔⒑途W(wǎng)格尺度信息關(guān)聯(lián)起來，構(gòu)造fk在計(jì)算域內(nèi)的分布函數(shù)，自動調(diào)節(jié)fk在不同流動區(qū)域的取值.計(jì)算結(jié)果表明，可變fk方法具有求解小尺度湍流結(jié)構(gòu)的能力，所得渦結(jié)構(gòu)豐富，流場正確合理.通過流場平均量及脈動量的對比分析，可變fk方法與SST DES計(jì)算結(jié)果及實(shí)驗(yàn)結(jié)果均吻合較好，相較于統(tǒng)一fk方法體現(xiàn)出了較為明顯的優(yōu)勢，是一種具有工程應(yīng)用前景的新型RANS/LES混合方法.

3)下一步研究工作中，考慮改進(jìn)及構(gòu)造新型的fk分布函數(shù)，以更為準(zhǔn)確地模擬復(fù)雜湍流問題，并通過對方法更加深入地考察評估，將其進(jìn)一步推向工程應(yīng)用.

References)

［1］閻超，于劍，徐晶磊，等.CFD模擬方法的發(fā)展成就與展望［J］.力學(xué)進(jìn)展，2011，41(5):562-589.Yan C，Yu J，Xu J L，et al.On the achievements and prospects for the methods of computational fluid dynamics［J］.Advances in Mechanics，2011，41(5):562-589(in Chinese).

［2］ Spalart P R，Jou W H，Strelets M，et al.Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach［C］∥Proceedings of 1st AFOSR International Conference on DNS/LES，Advances in DNS/LES.Columbus:Greyden Press，1997:137-147.

［3］ Menter F R，Kuntz M，Bender R.A scale-adaptive simulation model forturbulent flow predictions，AIAA-2003-0767［R］.Reston:AIAA，2003.

［4］ Girimaji S S，Srinivasan R，Jeong E.PANS turbulence model for seamless transition between RANS and LES:Fixed point analysis and preliminary results［C］∥Proceedings of ASME/JSME 2003 4th Joint Fluids Summer Engineering Conference.New York:ASME，2003，2:1901-1909.

［5］ Batten P，Goldberg U，Chakravarthy S.Reconstructed sub-grid methods for acoustics predictions at all Reynolds numbers，AIAA-2002-2511［R］.Reston:AIAA，2002.

［6］ Frendi A，Tosh A，Girimaji S S.Flow past a backward-facing step:Comparison of PANS，DES and URANS results with experiments［J］.International Journal for Computational Methods in Engineeing Science and Mechanics，2007，8(1):23-38.

［7］ Lakshmipathy S，Reyes D A，Girimaji S S.Partially averaged Navier-Stokes method:Modeling and simulation of low Reynolds number effects in flow past a circular cylinder，AIAA-2011-3107［R］.Reston:AIAA，2011.

［8］ Basu D，Hamed A，Das K.Assessment of partially averaged Navier-Stokes(PANS)multiscale model in transonic turbulent separated flows［C］∥Proceedings of FEDSM，5th Joint ASME/JSME Fluids Engineering Conference.New York:ASME，2007，2:1451-1459.

［9］ Lockard D P.Summary of the tandem cylinder solutions from the benchmark problems for airframe noise computations-I workshop，AIAA-2011-0353［R］.Reston:AIAA，2011.

［10］ Schwamborn D.Results and lessons learned from the EU-project ATAAC［C］∥ERCOFTAC Symposium Unsteady Separation in Fluid-Structure Interaction.Mykonos，Greece:［s.n.］2013，6:17-21.

［11］ Jenkins L N，Khorrami M R，Choudhari M M，et al.Characterization of unsteady flow structures around tandem cylinders for component interaction studies in airframe noise，AIAA-2005-2812［R］.Reston:AIAA，2005.

［12］ Hutcheson F V，Brooks T F.Noise radiation from single and multiple rod configurations，AIAA-2006-2629［R］.Reston:AIAA，2006.

［13］ Jenkins L N，Neuhart D H，Mcginley C B，et al.Measurements of unsteady wake interference between tandem cylinders，AIAA-2006-3202［R］.Reston:AIAA，2006.

［14］ Lockard D P，Khorrami M R，Choudhari M M，et al.Tandem cylinder noise predictions，AIAA-2007-3450［R］.Reston:AIAA，2007.

［15］ Weinmann M，Sandberg R D，Doolan C.Flow and noise prediction for a tandem cylinder configuration using novel hybrid RANS/LES approaches，AIAA-2010-3787［R］.Reston:AIAA，2010.

［16］ Xiao Z X，Liu J，Huang J B，et al.Numerical dissipation effects on massive separation around tandem cylinders［J］.AIAA Journal，2012，50(5):1119-1136.

［17］ Menter F R.Two equation eddy-viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32(8):1598-1605.

［18］ Luo D H，Yan C，Liu H K，et al.Comparative assessment of PANS and DES for simulation of flow past a circular cylinder［J］.Journal of Wind Engineering ＆ Industrial Aerodynamics，2014，134:65-77.

［19］劉佳，閻超，趙瑞，等.DES方法對返回艙氣動熱環(huán)境數(shù)值模擬［J］.北京航空航天大學(xué)學(xué)報(bào)，2013，39(5):590-594.Liu J，Yan C，Zhao R，et al.Simulation of re-entry capsule thermodynamics environment by DES method［J］.Journal of Beijing University of Aeronautics and Astronautics，2013，39(5):590-594(in Chinese).

［20］李君哲，閻超，柯倫，等.氣動熱CFD計(jì)算的格式效應(yīng)研究［J］.北京航空航天大學(xué)學(xué)報(bào)，2003，29(11):1022-1025.Li J Z，Yan C，Ke L，et al.Research on scheme effect of computational fluid dynamics in aerothermal［J］.Journal of Beijing University of Aeronautics and Astronautics，2003，29(11):1022-1025(in Chinese).