SST-DES在小分離流動(dòng)數(shù)值模擬中的改進(jìn)

王翔宇 李 棟

(西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力學(xué)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室，西安710072)

1 問題提出

分離渦模擬［1］(DES，Detached Eddy Simulation)是當(dāng)前模擬復(fù)雜航空航天工程問題中分離流動(dòng)的一種主要方法.作為一種RANS/LES的混合方法，DES將湍流模型中的長(zhǎng)度尺度表征為一種壁面距離函數(shù)控制器，邊界層內(nèi)流動(dòng)仍受雷諾平均方法 (RANS，ReynoldsAveragedNavier Stokes)湍流模型主導(dǎo)，并未顯著提高計(jì)算量，而在大尺度分離區(qū)或遠(yuǎn)場(chǎng)簡(jiǎn)化為亞格子模型大大降低了湍流黏性，起到類似大渦模擬方法(LES，Large Eddy Simulation)隱式濾波的效果.早期的DES的一個(gè)主要問題是模型應(yīng)力損耗(MSD，Modeled Reynolds Stress Depletion)導(dǎo)致的網(wǎng)格誘導(dǎo)分離現(xiàn)象(GIS，Grid Induced Separation)，邊界層內(nèi)受到LES方法干擾湍流黏性降低，在小分離流動(dòng)中出現(xiàn)明顯的非物理分離解［2］.為此，文獻(xiàn)［3］提出了延遲分離渦模擬(DDES，Delayed Detached Eddy Simulation)，通過引入控制函數(shù)保證邊界層內(nèi)的 RANS求解，延遲了 RANS的作用范圍.

然而很多時(shí)候除了壁面上的分離點(diǎn)，對(duì)邊界層下游的自由剪切層(如翼型的尾跡區(qū))流動(dòng)信息同樣關(guān)注.從這一點(diǎn)來看，DDES對(duì)邊界層內(nèi)RANS的延遲是以犧牲自由剪切層的LES可解性為代價(jià)的——延遲得到的較大湍流黏性可能從RANS/LES交界面開始向LES區(qū)作用，干擾了LES的分辨能力，出現(xiàn)了所謂的灰區(qū)［4-5］(grey area)，如圖1所示.

圖1 理想的DES模擬效果和實(shí)際的DES模擬效果［6］Fig.1 Comparison between the ideal situation and actual situation in DES［6］

從流動(dòng)特性考慮，對(duì)大分離流動(dòng)，漩渦的劇烈生成發(fā)展衰減耗散的過程同樣是湍動(dòng)能不斷發(fā)展的過程，這個(gè)貫穿整個(gè)流場(chǎng)的過程將流動(dòng)發(fā)展初期的可解信息不斷傳播循環(huán)到整個(gè)LES區(qū)，形成一個(gè)自反饋的機(jī)制使得灰區(qū)變得模糊不清甚至忽略不計(jì)［4］.但對(duì)小分離流動(dòng)而言，較弱的漩渦信息傳遞使得這種反饋機(jī)制幾乎不復(fù)存在，分離越弱，RANS/LES交界區(qū)域受上游邊界層RANS影響越大，灰區(qū)效應(yīng)越明顯.基于以上，本文旨在發(fā)展一種能夠更為準(zhǔn)確表征LES區(qū)域的改進(jìn)SSTDES方法，降低小分離流動(dòng)中灰區(qū)對(duì)自由剪切層的不利影響.

2 SST-DES在小分離流動(dòng)中的改進(jìn)

其中，μt和μl分別為湍流黏性系數(shù)和層流黏性系數(shù);Ω為當(dāng)?shù)亓鲃?dòng)的應(yīng)變率張量;表征RANS的?；牧鏖L(zhǎng)度尺度;Δ為某個(gè)體現(xiàn)LES特性的選定的網(wǎng)格長(zhǎng)度濾波尺度;fDDES為DDES中采用的延遲函數(shù)——本文采用已經(jīng)得到廣泛應(yīng)用的、SST模型中自帶的 f2函數(shù)［8-9］.式(1)～式(3)中其他參數(shù)均為標(biāo)準(zhǔn)k-ω SST湍流模型給定，可參考文獻(xiàn)［7］.

而Yoshizawa提出的 k方程亞格子模式［10］(k-epuation sub-grid model)的表達(dá)式為

其中Cs=0.05，Cd=1.0均為該亞格子模式中給定的常數(shù).觀察SST-DES的表達(dá)式，在遠(yuǎn)離壁面時(shí)僅變更式(2)的耗散項(xiàng)，轉(zhuǎn)換為類似式(7)k方程亞格子模式的方程構(gòu)型，起到類似大渦模擬的作用——但此時(shí)湍流黏性系數(shù)仍由SST中k-ε湍流模型表達(dá)式控制，即式(1)為

最終的湍流黏性仍受較大的RANS長(zhǎng)度尺度控制，混淆了此時(shí)式(5)對(duì)湍動(dòng)能方程的影響.從另一個(gè)角度看，當(dāng)從RANS區(qū)域轉(zhuǎn)變?yōu)長(zhǎng)ES區(qū)域時(shí)，即使湍動(dòng)能k能夠迅速減小，由于ω方程并未變化，上游壁面附近得到的較低湍動(dòng)能耗散比率ω會(huì)在交界處顯著抑制LES區(qū)域的解析能力，延緩?fù)牧鞯陌l(fā)展.為此有必要將式(8)改進(jìn)為標(biāo)準(zhǔn)的k方程亞格子模式中的式(6)構(gòu)型:

為防止出現(xiàn)類似DES中直接引入亞格子模式導(dǎo)致的網(wǎng)格誘導(dǎo)分離現(xiàn)象，借鑒DDES的思路，通過引用混合函數(shù)的方法延遲RANS作用區(qū)域.最終式(1)可寫為

這樣，在遠(yuǎn)場(chǎng)(f2=0)將原 SST-DES中的“類”k方程亞格子模式轉(zhuǎn)化為其標(biāo)準(zhǔn)形式，阻斷了之前變量ω在RANS和LES之間的信息傳遞.進(jìn)一步分析，當(dāng)SST-DES起到LES作用時(shí)，將式(9)代入式(11)，近似認(rèn)為f2=0，化簡(jiǎn)后得到湍流黏性系數(shù)的對(duì)流-擴(kuò)散方程:

限于篇幅，上式中復(fù)雜的擴(kuò)散項(xiàng)表達(dá)式在此略去.觀察式(12)，一方面由于LES濾波尺度Δ出現(xiàn)在了生成項(xiàng)，而生成項(xiàng)又在源項(xiàng)中處于主導(dǎo)地位，因此Δ的表達(dá)形式比其在一般的DES中更為重要.本文選擇k方程亞格子模式中默認(rèn)的Δ:

代替原 SST-DES 中的 Δ =max(Δx，Δy，Δz)［11-12］.

另一方面，當(dāng)式(10)處于生成耗散平衡時(shí):

化簡(jiǎn)得到

而SA-DES中采用的Smagorinksy亞格子模式為

其中C=0.182.顯然最終湍流黏性系數(shù)的表達(dá)式(15)與Smagorinksy亞格子模式方程(16)構(gòu)型完全相同.同時(shí)，將k方程亞格子模型中的Cs和Cd代入式(15)后發(fā)現(xiàn)與式(16)中的常數(shù)并不相等，需要對(duì)二者取值重新標(biāo)定.通過對(duì)各項(xiàng)同性湍流模擬結(jié)果進(jìn)行比較，本文算例中選取Cs=0.1，Cd=1.0，保證了式(15)右端的常系數(shù)約為Smagorinksy亞格子模式中給定的0.182.

3 計(jì)算模型和結(jié)果分析

本文仿照Wei和Ning［13］在提出改進(jìn)SA-DES時(shí)采用的AS239翼型最大升力點(diǎn)繞流進(jìn)行算例驗(yàn)證.實(shí)驗(yàn)顯示在來流條件為Ma=0.15，Re=2×106時(shí)，當(dāng)迎角為13.3°，AS239翼型達(dá)到最大升力點(diǎn)，分離發(fā)生在翼型尾部且分離渦很?。?4-15］.計(jì)算網(wǎng)格在二維流向切面為323×121，壁面法向第1層網(wǎng)格距離約為1×10-5，y+≈1.翼型展向在0.4c長(zhǎng)度內(nèi)拉伸30層.選用基于LU-SGS的雙時(shí)間步長(zhǎng)推進(jìn)方法，物理時(shí)間步為0.01c/U∞.分別采用SST-DES，SST-DDES和上文提出的改進(jìn)方法(簡(jiǎn)記為DES，DDES和I-DES)進(jìn)行數(shù)值模擬.

觀察圖2，相比 DES，DDES通過延遲 RANS提高了壁面附近湍流黏性，但同時(shí)也使得尾跡區(qū)湍流黏性大幅提高，幾乎達(dá)到了和壁面附近相同的量級(jí)，可能進(jìn)一步加劇灰區(qū)的影響.I-DES的結(jié)果在邊界層與DDES類似，而在尾跡區(qū)則幾乎“截?cái)唷绷薉DES過大的湍流黏性，促進(jìn)了湍流的快速發(fā)展.進(jìn)一步從渦量云圖可以看出，由于網(wǎng)格誘導(dǎo)分離DES在翼型中部即出現(xiàn)了明顯的大分離漩渦并不斷向下游發(fā)展.DDES雖然避免了在翼型表面的提前分離但延遲RANS使得尾跡區(qū)漩渦發(fā)展受到抑制變得模糊不清.相比之下，I-DES在翼型表面渦量分布與DDES較為類似而在后緣的自由剪切層的漩渦運(yùn)動(dòng)更加清晰.

圖3和圖4分別給出了摩擦阻力系數(shù)分布和翼型表面壓力系數(shù).相比于實(shí)驗(yàn)得到分離點(diǎn)約0.82c，DES僅為0.7c左右，翼型尾部邊界層內(nèi)模化應(yīng)力嚴(yán)重不足，而DDES和I-DES分離點(diǎn)均在0.8c～0.85c之間，也就是說I-DES這種混合函數(shù)構(gòu)造的形式起到了類似DDES改善DES中網(wǎng)格誘導(dǎo)分離的效果.對(duì)表面壓力分布而言，在尾部分離區(qū)I-DES與實(shí)驗(yàn)數(shù)據(jù)幾乎完全吻合.

圖5給出了尾跡區(qū)前部3個(gè)站位處的流向速度分布曲線.

圖2 DES，DDES和I-DES的湍流黏性系數(shù)和渦量對(duì)比Fig.2 Comparison between DES，DDES and I-DES in turbulent viscosity and vorticity magnitude

圖3 摩擦阻力系數(shù)曲線Fig.3 Skin friction coefficient along the suction side

圖4 翼型表面壓力系數(shù)曲線Fig.4 Pressure coefficient along the aerofoil

圖5 翼型尾跡區(qū)流向速度分布比較Fig.5 Comparison of streamwise velocity profiles in the wake region

由圖5不難看出，受上游RANS區(qū)影響，DES和DDES在尾跡區(qū)特別是靠近在翼型的位置湍流黏性偏大，模化剪切應(yīng)力偏大，直觀來看流向“黏滯”強(qiáng)烈，速度型梯度比I-DES大很多——由于延遲RANS，灰區(qū)影響更強(qiáng)，DDES在尾跡區(qū)流向比DES更偏離實(shí)驗(yàn)結(jié)果，而這種偏離隨著站位的靠后、湍流的發(fā)展而不斷減小.顯然，在小分離流尾跡區(qū)I-DES取得了最靠近實(shí)驗(yàn)值的速度分布結(jié)果.

4 結(jié)論

1)對(duì)小分離流動(dòng)而言，漩渦傳遞湍動(dòng)信息能力較弱，SST-DES中下游LES區(qū)域缺乏上游附體流區(qū)域解析湍流信息，受RANS區(qū)域較大湍流黏性干擾嚴(yán)重，在RANS/LES交界附近湍流發(fā)展受到抑制，即出現(xiàn)了所謂灰區(qū)效應(yīng);

2)通過對(duì)原始SST-DES方程的分析，降低了其中類亞格子模式對(duì)湍流黏性的多參數(shù)影響，促進(jìn)了自由剪切層湍流的發(fā)展，用標(biāo)準(zhǔn)k方程亞格子模式將SST-DES重新構(gòu)造，得到了一種改進(jìn)的SST-DES方法;

3)通過對(duì)AS239翼型最大升力點(diǎn)小分離流動(dòng)的數(shù)值模擬，改進(jìn)后的SST-DES方法在翼型表面具有類似DDES延遲RANS的特性，避免出現(xiàn)DES中網(wǎng)格誘導(dǎo)分離現(xiàn)象，同時(shí)在尾跡區(qū)大大降低了湍流黏性，明顯減輕了灰區(qū)的影響，驗(yàn)證了該改進(jìn)方法的正確性.

References)

［1］ Strelets M.Detached eddy simulations of massively separated flows［R］.AIAA-2001-0879，2001

［2］ Spalart P R.Young-person's guide to detached-eddy simulation grids［R］.NASA-TM-2001-211032，2001

［3］ Spalart P R，Deck S，Shur M L，et al.A new version of detachededdy simulation，resistant to ambiguous grid densities［J］.Theoretical and Computational Fluid Dynamics，2006，20(3):181-195

［4］ Spalart P R.Detached-eddy simulation［J］.Annual Review of Fluid Mechanics，2009，41:181-202

［5］ Menter F，Kuntz M.Adaption of eddy-viscosity turbulence models to unsteady separated flow behind vehicles［J］.The Aerodynamics of Heavy Vehicles:Trucks，Buses，and Trains，2004，19:339-352

［6］ Michel U，Eschricht D，Greschner B，et al.Advanced DES methods and their application to aeroacoustics［C］//Progress in HybridRANS-LES Modelling.Berlin:Springer-Verlag，2010:59-76

［7］ Mockett C.A comprehensive study of detached-eddy simulation［D］.Germany:Berlin Institute of Technology，2009

［8］ Rochm M，Krta M，Meson S.Zonal hybrid RANS-LES method for static and oscillating airfoils and wings［R］.AIAA-2006-1256，2006

［9］ Lynch C E，Smithy M J.Hybrid RANS-LES turbulence models on unstructured grids［R］.AIAA-2008-3854，2008

［10］ Yoshizawa A，Horiuti K.A statistically-derived subgrid scale kinetic energy model for the large-eddy simulation of turbulent flows［J］.Journal of the Physical Society of Japan，1985，54(8):2834-2839

［11］ Deck S.Zonal-detached eddy simulation of the flow around a high-lift configuration［J］.AIAA Journal，2005，43(11):2372-2384

［12］ Deck S.Recent improvements of the zonal detached eddy simulation(ZDES)formulation［J］.Theoretical and Computational Fluid Dynamics，2012，26(6):523-550

［13］ Wei W，Ning Q.Balancing destruction and production in S-A model-based hybrid RANS-LES for flow around an aerofoil with mild separation［C］//Progress in Hybrid RANS-LES Modelling.Berlin:Springer-Verlag，2010:379-388

［14］ Schmidtf S，Thiele R.Detached eddy simulation of flow around A-airfoil［J］.Flow，Turbulenceand Combustion，2003，71(1-4):261-278

［15］ Gleyzes C，Capbern P.Experimental study of two AIRBUS/ONERA airfoils in near stall conditions.Part I:boundary layers［J］.Aerospace Science and Technology，2003，7(6):439-449