杜一鳴 高正紅 舒博文 邱福生,1) 宋辰星

* (沈陽航空航天大學(xué)航空宇航學(xué)院,沈陽 110136)

? (西北工業(yè)大學(xué)航空學(xué)院,西安 710072)

** (中國航發(fā)沈陽發(fā)動機(jī)研究所,沈陽 110015)

引言

基于雷諾平均N-S (Reynolds-averaged Navier-Stokes,RANS) 框架的渦黏性湍流模式(eddyviscosity model,EVM)經(jīng)過40 多年的發(fā)展對于邊界層附著流動、中等壓力梯度和小分離流動以及自由剪切流動都能得到較理想的模擬結(jié)果[1-2],在廣泛的工程實踐中得到了比較充分的驗證.雖然計算機(jī)技術(shù)推動了直接數(shù)值模擬(DNS) 和大渦模擬(LES)等方法的發(fā)展,但RANS 方法被認(rèn)為仍是未來幾十年工程應(yīng)用中的主要模擬手段[3-4].然而由于渦黏性模式本身存在的基礎(chǔ)理論、模式構(gòu)造和計算方法等方面的問題,導(dǎo)致其對于一些工程常見的分離和剪切流動難以給出滿足氣動設(shè)計需求的模擬結(jié)果[3-4].Kato 和Launder[5]就曾指出,基于傳統(tǒng)渦黏性假設(shè)的湍流模式實際在除了簡單剪切流動中的表現(xiàn)都不好.

雖然輸運(yùn)方程在一定程度上體現(xiàn)了湍流的輸運(yùn)特性,但EVM 中渦黏性系數(shù)僅采用湍流量定義的方式(即 νT=k/ω) 實際只有在湍流平衡流動中才成立[1].平衡流(equilibrium flow)[6]指的是湍動能生成與耗散達(dá)到局部平衡,且相對于對流和輸運(yùn)足夠大的流動狀態(tài).在這種流動區(qū)域中,流動狀態(tài)(例如雷諾應(yīng)力等)完全由當(dāng)?shù)貤l件決定,而基本不存在上游歷史效應(yīng).邊界層內(nèi)層就基本具有這個性質(zhì).實際上,對于滿足局部平衡的湍流,包括Prantdl 混合長度理論[7]在內(nèi)的大多數(shù)湍流模式都能得到較好的模擬結(jié)果.但整個邊界層內(nèi)并不完全滿足局部平衡性,即湍流非平衡流動(non-equilibrium flow),尤其是當(dāng)存在逆壓梯度時,此時采用傳統(tǒng)EVM 的渦黏性系數(shù)定義方式就會得到比實際更大的雷諾應(yīng)力預(yù)測結(jié)果,進(jìn)而導(dǎo)致分離較小,最典型的流動就是激波?邊界層分離(shock wave-boundary layer separation).

Coakley[8]于1983 年通過修正渦黏性系數(shù)的方式首先將“雷諾應(yīng)力輸運(yùn)”(Reynolds-stress transport)的思想引入兩方程湍流模式.1985 年出現(xiàn)的Johnson-King (J-K)零方程模式[9]通過限制速度虧損律層內(nèi)的主雷諾應(yīng)力符合Bradshaw 假設(shè)[10-11],降低了渦黏性系數(shù)從而提高了模式逆壓梯度流動的預(yù)測能力.這種方法向方程中引入了“遲滯效應(yīng)”[12],相當(dāng)于考慮了雷諾應(yīng)力的輸運(yùn)特性.后來,Menter[13]在1992 年提出k-ω切應(yīng)力輸運(yùn)(shear-stresstransport,SST)兩方程模式時參考J-K 模式的處理方式引入了Bradshaw 假設(shè),基本解決了諸如RAE2822超臨界翼型等二維激波?邊界層分離的預(yù)測問題.然而對于較強(qiáng)的三維激波?邊界層分離流動,目前的EVM 均無法進(jìn)行正確模擬,其激波和分離預(yù)測位置向上游移動,嚴(yán)重偏離實驗結(jié)果.以較為經(jīng)典的ONERA M6 機(jī)翼[14]為例,跨聲速α=6.06°狀態(tài)下激波相比于小攻角狀態(tài)更強(qiáng),同時展向更大范圍內(nèi)都存在邊界層分離,包括k-ωSST 和Spalart-Allmaras(S-A)[15]模式在內(nèi)的EVM 都無法得到合理的計算結(jié)果[4].近年來的計算分析[4,16-19]表明,采用雷諾應(yīng)力模式(Reynolds-stress model,RSM)才能夠?qū)ι鲜龉r進(jìn)行較準(zhǔn)確的預(yù)測.RSM 直接求解各雷諾應(yīng)力分量的輸運(yùn)方程,不存在雷諾應(yīng)力各向同性和非平衡流動模擬問題[19].但此類二階封閉模式建模難度較大,且目前數(shù)值穩(wěn)定性和收斂性仍未達(dá)到工程應(yīng)用的水平[3-4].

由于湍流非平衡流動的復(fù)雜性,目前針對渦黏性模式這一問題的研究較少,尚未形成廣泛有效的修正策略.Li 等[20-21]的研究具有代表性,他們在開發(fā)渦黏性湍流模式結(jié)冰翼分離模擬能力時也考慮了非平衡性的影響,并以湍動能生成和耗散之比作為開關(guān)參數(shù)針對ω方程的耗散項進(jìn)行了修正,取得了較好的改進(jìn)效果,但并未進(jìn)行激波分離流動方面的驗證討論.Rumsey 和Vatsa[22]針對ONERA M6 機(jī)翼α=5.06°和α=6.06°狀態(tài)對包括J-K 零方程模式和S-A 一方程模式在內(nèi)的幾種湍流模式進(jìn)行了網(wǎng)格和格式耗散的收斂性研究,發(fā)現(xiàn)采用耗散較大的數(shù)值方法(例如矢通量分裂格式FVS)以及稀網(wǎng)格能夠得到收斂且與實驗吻合較好的計算結(jié)果.一旦降低數(shù)值格式耗散水平或加密網(wǎng)格,計算就會得到不穩(wěn)定的大分離,激波位置大幅向前緣移動.Frink[23]采用S-A 模式在稀網(wǎng)格上也得到了ONERA M6 機(jī)翼α=6.06°狀態(tài)y/b=90%截面較好的壓力分布模擬結(jié)果,但他認(rèn)為這只是“某種巧合”,因為由稀網(wǎng)格引入的截斷誤差還很大,同時“翼尖附近復(fù)雜的流動結(jié)構(gòu)可能已經(jīng)超過了S-A 模式的模擬能力”.然而上述研究并未對EVM 三維激波?邊界層分離模擬失效的具體原因作出分析.

從機(jī)理角度出發(fā),實際上,對于由逆壓梯度導(dǎo)致的邊界層分離流動,邊界層內(nèi)的雷諾應(yīng)力水平基本決定了其抵抗逆壓梯度的能力以及分離流動的特性,因此湍流模式對雷諾應(yīng)力的求解直接影響其激波-邊界層分離流動的模擬精度.在渦黏性假設(shè)的框架下,雷諾應(yīng)力本構(gòu)關(guān)系和渦黏性系數(shù)定義方式就是兩個突破點.

以目前較常用的k-ωSST 兩方程湍流模式為修正對象,首先通過典型二維和三維算例對模式的激波?邊界層分離模擬能力進(jìn)行評估,分析三維激波?邊界層分離的物理特性和模式失效機(jī)理,并討論現(xiàn)有非線性雷諾應(yīng)力本構(gòu)關(guān)系對于此類問題的改進(jìn)效果.在此基礎(chǔ)上,提出兩種針對渦黏性系數(shù)的修正方法,并通過與原始模式和RSM 的雷諾應(yīng)力特性對比、網(wǎng)格收斂性分析等確認(rèn)修正方法進(jìn)行三維激波?邊界層分離模擬的合理性和有效性.

1 模擬能力評估

1.1 k-ω SST 模式與Bradshaw 假設(shè)

相比于k-ε模式,k-ω模式更適合邊界層流動模擬,對逆壓梯度主導(dǎo)的附著流動和小分離流動有較高的模擬精度,但對湍流變量的遠(yuǎn)場邊界條件比較敏感[1].Menter[13]基于k-ω模式和k-ε模式各自的特點,巧妙地將兩種模式相融合,并在邊界層內(nèi)外通過開關(guān)函數(shù)(blending/switch function)進(jìn)行切換.直接將兩種模式形式和參數(shù)進(jìn)行混合的稱為基準(zhǔn)模式(baseline model,BSL).Menter 認(rèn)為EVM 與RSM 的主要區(qū)別在于EVM 僅考慮了湍動能和耗散率的輸運(yùn)特性,而未直接考慮雷諾應(yīng)力本身的輸運(yùn)效應(yīng).而對于非平衡流動,渦黏性系數(shù)的定義缺陷將導(dǎo)致雷諾應(yīng)力不再符合其輸運(yùn)特性,進(jìn)而得到更大的雷諾應(yīng)力預(yù)測結(jié)果,使分離預(yù)測偏小.針對這一問題,k-ωSST 模式在BSL 模式的基礎(chǔ)上引入了Bradshaw 假設(shè)(不可壓縮形式)[10-11]

即對于邊界層、尾跡和混合層流動,主雷諾應(yīng)力分量與湍動能k的比值為常數(shù).實驗測量表明,常數(shù)a1基本保持在0.3 左右,稱為Bradshaw 常數(shù)[10]或Townsend 常數(shù)[11].Bradshaw 假設(shè)被公認(rèn)能夠較準(zhǔn)確地反映主雷諾應(yīng)力分量在對數(shù)律層和速度虧損律層(尾跡區(qū))中的特征.根據(jù)湍流邊界層分層結(jié)構(gòu),對數(shù)律層和速度虧損律層是受勢流區(qū)壓力梯度影響最大的區(qū)域,通常湍流模式在對數(shù)律層的表現(xiàn)對于準(zhǔn)確預(yù)測中等壓力梯度流動非常重要,而模式對于邊界層尾跡(速度虧損律層)的預(yù)測精度決定了模式對強(qiáng)逆壓梯度流動的預(yù)測能力[1,24].

當(dāng)湍動能生成項大于耗散項時有Sˉ >a1ω,上式就會保證雷諾應(yīng)力不超過 ρa(bǔ)1k.其中,開關(guān)函數(shù)F2用于限制Bradshaw 假設(shè)僅在邊界層和尾跡內(nèi)起作用.相比之下,k-ωBSL 模式僅采用 νT=k/ω 的原始定義方式.加入如式(3)的渦黏性系數(shù)限制后,相當(dāng)于強(qiáng)制雷諾應(yīng)力在逆壓梯度下滿足Bradshaw假設(shè),此時雷諾應(yīng)力不會被過度預(yù)測,能夠顯著提升模式的二維激波?邊界層分離模擬能力.圖1 給出了采用不同湍流模式及圖2 網(wǎng)格(y+<2)[25]得到的RAE2822超臨界翼型[26]Ma=0.75,Re=6.2 × 106,α=2.72°(AGARD Case 10)狀態(tài)[27-28]激波?邊界層分離的計算壓力系數(shù)和摩擦系數(shù)(下表面Cf取相反數(shù))對比.可以發(fā)現(xiàn),SST 模式與RSM 計算結(jié)果基本一致,激波位置和強(qiáng)度與實驗數(shù)據(jù)吻合較好,BSL 模式模擬的激波位置靠后,且分離泡較小.

圖1 計算壓力系數(shù)和摩擦系數(shù)分布與實驗的對比Fig.1 Comparison of computed pressure and skin friction coefficients distribution with experiment

圖2 RAE2822 翼型計算網(wǎng)格Fig.2 Computational grid of RAE2822 airfoil

本文RANS 計算均采用NASA CFL3 D 結(jié)構(gòu)化網(wǎng)格求解器[27]完成.該求解器在有限體積框架下求解三維可壓縮非定常RANS 方程,有多種湍流模式可供選擇,計算效率和計算精度滿足工程應(yīng)用需求.所采用的k-ωSST 模式為2003 年Menter 等[29]提出的改進(jìn)版本.此外,RSM 計算結(jié)果將用作對比分析,采用的是Togiti 和Eisfeld[30]于2015 年提出的SSG/LRR-g模式,其中g(shù)≡ (1/ω)1/2為長度尺度相關(guān)量.

1.2 三維激波-邊界層分離算例驗證

ONERA M6 機(jī)翼[14](圖3)包含了復(fù)雜的激波?邊界層分離流動特征,同時外形簡單且實驗數(shù)據(jù)比較完善,是CFD 研究的經(jīng)典算例.其中Ma=0.84,α=3.06°狀態(tài)對應(yīng)附著流動,機(jī)翼上表面的λ 形激波結(jié)構(gòu)是其主要特征.隨著攻角的增加,上表面分離區(qū)逐漸增大,當(dāng)α=6.06°時,外翼段分離已經(jīng)非常明顯.本文選取的計算狀態(tài)在表1 中給出.計算網(wǎng)格規(guī)模為510 萬(圖4,y +≈ 1).

圖4 ONERA M6 計算網(wǎng)格Fig.4 Computational grid of ONERA M6 wing

表1 ONERA M6 機(jī)翼計算狀態(tài)Table 1 Computation conditions of ONERA M6 wing

圖3 ONERA M6 機(jī)翼形狀與模式參數(shù)Fig.3 Shape and parameter of ONERA M6 wing

圖5 給出了ONERA M6 機(jī)翼α=3.06°和α=4.08°狀態(tài)(分別對應(yīng)原始實驗Case 2308 和Case 2563)典型截面計算壓力分布與實驗的對比.可以看出,對于分離不強(qiáng)的小攻角流動狀態(tài),k-ωSST 模式具有比較理想的預(yù)測精度,計算結(jié)果與實驗吻合很好,且與BSL 模式和RSM 的計算結(jié)果差別較小.

圖5 小攻角狀態(tài)計算壓力分布與實驗的對比Fig.5 Comparison of computed pressure distribution with experiment at small angles of attack

然而,對于更大攻角下的強(qiáng)激波?大分離流動,目前大多數(shù)渦黏性模式都無法勝任.圖6 和圖7 給出了α=6.06°狀態(tài)不同湍流模式計算的上表面流動形態(tài)及典型截面壓力分布.可以發(fā)現(xiàn),k-ωSST 模式得到的激波位置從y/b=65%開始大幅向前緣移動,與RSM 的計算結(jié)果和實驗結(jié)果的偏離程度不斷增大[17,23].在α=5.06°狀態(tài)下,這種情況發(fā)生在y/b=80%之后[23].

圖6 表面壓力分布與極限流線對比Fig.6 Comparison of surface pressure distribution and limited streamline pattern

圖7 α=6.06°狀態(tài)計算壓力分布與實驗的對比(非線性模式與RSM、SST 和BSL 模式)Fig.7 Comparison of computed pressure distribution with experiment at α=6.06° (nonlinear models with RSM,SST and BSL models)

2 失效機(jī)理分析

為進(jìn)一步確認(rèn)造成渦黏性湍流模式三維激波?邊界層分離模擬失效的原因,對RSM 和k-ωSST 兩種模式進(jìn)行網(wǎng)格收斂性分析,自研網(wǎng)格序列包括微量網(wǎng)格(T,網(wǎng)格量0.74 百萬)、稀網(wǎng)格(C,網(wǎng)格量1.35 百萬)、中等網(wǎng)格(M,網(wǎng)格量3.3 百萬)、密網(wǎng)格(F,網(wǎng)格量5.1 百萬,圖4)和特密網(wǎng)格(EF,網(wǎng)格量11.0 百萬);物面法向首層網(wǎng)格高度均為y +≈ 1,這參考了Frink[23]的研究發(fā)現(xiàn),即法向網(wǎng)格密度相比于表面網(wǎng)格密度對結(jié)果的影響很小.從圖8 中可以發(fā)現(xiàn),RSM 在所有截面都表現(xiàn)出了較好的網(wǎng)格收斂趨勢,激波模擬結(jié)果隨網(wǎng)格加密逐漸陡峭且更加接近實驗結(jié)果(但在EF 網(wǎng)格上收斂性較差,未完成計算).對于k-ωSST 模式,內(nèi)翼段(y/b<65%)也表現(xiàn)出了正常的網(wǎng)格收斂趨勢,但在外翼段卻表現(xiàn)出隨著網(wǎng)格加密激波位置前移、逐漸偏離實驗結(jié)果的趨勢,這與Rumsey 和Vatsa[22]的結(jié)論一致.稀網(wǎng)格的計算結(jié)果也表明,某些看似正確的數(shù)值模擬結(jié)果其實是具有誤導(dǎo)性的,因為“實際上整個計算的截斷誤差還未通過加密網(wǎng)格或降低數(shù)值格式耗散得到充分降低”[22].考慮到外翼段激波位置即是分離點,造成上述預(yù)測失準(zhǔn)的一種合理解釋是:外翼段激波?邊界層分離流動本身是由分離主導(dǎo)的,對于傳統(tǒng)EVM,雷諾應(yīng)力預(yù)測偏小導(dǎo)致分離靠前;分離后物面壓力恢復(fù)嚴(yán)重受阻進(jìn)而導(dǎo)致激波的產(chǎn)生,也即分離點直接決定了激波位置.因此,正確預(yù)測三維激波?邊界層分離的關(guān)鍵首先在于正確的雷諾應(yīng)力求解,而非網(wǎng)格或數(shù)值格式.而雷諾應(yīng)力求解又依賴模式的建模方法,具體來說就是雷諾應(yīng)力和渦黏性系數(shù)的定義.

圖8 α=6.06°狀態(tài)網(wǎng)格收斂性對比Fig.8 Comparison of grid convergence at α=6.06°

從雷諾應(yīng)力定義的角度考慮,基于渦黏性假設(shè)的線性雷諾應(yīng)力本構(gòu)關(guān)系是雷諾應(yīng)力求解的基礎(chǔ).完整的不可壓縮流動雷諾應(yīng)力本構(gòu)關(guān)系(constitutive relation)具有如下形式

計算穩(wěn)定性和精度較好.為了發(fā)展更簡潔魯棒的方法,Spalart[35]于2000 年針對SA 模式提出了“二次本構(gòu)關(guān)系”(quadratic constitutive relation,QCR)方法,在線性雷諾應(yīng)力本構(gòu)關(guān)系的基礎(chǔ)上增加了非線性項,即

對各向異性流動模擬有明顯改善,原理上適用于所有渦黏性模式.上述變量定義可參見文獻(xiàn)[34-35].由于S-A 模式本身不求解湍動能,在使用渦黏性假設(shè)時無法計入其中的湍動能項(即式(4)第三項).于是,2013 年Mani 等[36]針對S-A-QCR 方法引入了湍動能項的近似計算方法.2020 年Rumsey 等[37]通過近似項系數(shù)調(diào)整進(jìn)一步改進(jìn)了原始S-A-QCR修正方法.而基于湍動能方程的兩方程模式能夠直接求解湍動能,因此在構(gòu)造k-ωSST-QCR 方法時,后來加入的湍動能近似項可以不必考慮,直接采用原始形式的雷諾應(yīng)力本構(gòu)關(guān)系式(4)即可.

分別采用Rumsey-Gatskik-ωEASM 和k-ωSST-QCR 方法對ONERA M6 機(jī)翼α=6.06°狀態(tài)進(jìn)行模擬,計算結(jié)果見圖7.可以發(fā)現(xiàn),k-ωEASM 雖然起到了一定延緩分離的作用,但這種作用靠近外翼段逐漸消失,翼尖附近的計算結(jié)果與線性模式基本相同.而k-ωQCR 方法在該算例中基本未表現(xiàn)出明顯效果.從式(5)和式(6)的構(gòu)造形式可以看出,非線性雷諾應(yīng)力本構(gòu)關(guān)系考慮的是流動中旋轉(zhuǎn)、剪切等引起的非線性效應(yīng),因此適用于模擬雷諾應(yīng)力各向異性主導(dǎo)的流動(特別是二次流動),也能夠提升分離流動中雷諾應(yīng)力分量的計算精度[38],但對于改善湍流非平衡流動模擬能力并沒有實質(zhì)作用.

3 修正方法構(gòu)造

根據(jù)上述分析和討論,本節(jié)從渦黏性系數(shù)的定義角度出發(fā),分別基于Bradshaw 假設(shè)和湍流長度尺度提出兩種針對k-ωSST 模式的三維激波?邊界層分離模擬修正方法.

3.1 基于Bradshaw 假設(shè)的修正方法

從圖6 和圖7 的對比結(jié)果看,未考慮Bradshaw假設(shè)的BSL 模式反而得到了與RSM 接近的表面壓力分布和分離狀態(tài),在分離明顯的外翼段截面壓力分布也與實驗接近.由于Bradshaw 假設(shè)是利用Klebanoff 零壓力梯度平板邊界層試驗數(shù)據(jù)[10]提出的,因此只針對附著流動適用;隨后的研究表明,其在弱逆壓梯度和小分離流動中也基本可用[39],但可能不適用于其他流動狀態(tài)[24].本文的一種合理的推斷是,在三維強(qiáng)逆壓梯度作用下,各雷諾應(yīng)力分量的量級接近,因此Bradshaw 假設(shè)成立的基礎(chǔ)——“主雷諾應(yīng)力分量”的概念不復(fù)存在(即標(biāo)量近似形式不再適用),此時不能再假設(shè)僅有某個分量是主要的而忽略其他分量的作用.多個量級相當(dāng)?shù)睦字Z應(yīng)力分量會使Bradshaw 假設(shè)對雷諾應(yīng)力的限制(式(1))不再準(zhǔn)確(過度限制了雷諾應(yīng)力的生成),導(dǎo)致SST 模式在三維較強(qiáng)的激波分離流動中預(yù)測雷諾應(yīng)力嚴(yán)重不足,造成分離和激波位置靠前.但從圖7 還可以看出,BSL 模式在外翼段的模擬精度并不高,同時對于分離再附后的壓力分布計算精度較差,因此單純通過屏蔽Bradshaw 假設(shè)并不能解決問題,但嘗試對Bradshaw 假設(shè)進(jìn)行相應(yīng)修正不乏是一種可行措施.

針對高超聲速中的激波?邊界層干擾(shock wave-boundary layer interaction,SWBLI)模擬,近些年已有一些針對Bradshaw 假設(shè)的修正嘗試.Georgiadis 等[39-40]對Bradshaw 常數(shù)a1的實驗和理論基礎(chǔ)進(jìn)行了研究,并嘗試提高常數(shù)a1至0.34,0.355 甚至0.37,對多個超聲速SWBLI 問題進(jìn)行了模擬.結(jié)果表明,這種方式對提高SWBLI 問題模擬精度很有幫助.文獻(xiàn)[41]在針對SWBLI 問題的湍流模式不確定性分析中,也將常數(shù)a1的取樣范圍定為了0.31～ 0.40.根據(jù)Bradshaw 假設(shè)的作用機(jī)理,提高a1的目的在于放寬對雷諾應(yīng)力生成的限制,進(jìn)而提高雷諾應(yīng)力水平并推遲分離.

通過測試對比,將Bradshaw 常數(shù)a1調(diào)整為0.355 并對ONERA M6 機(jī)翼α=6.06°狀態(tài)進(jìn)行模擬,計算結(jié)果如圖9 和圖10 所示.可以發(fā)現(xiàn),采用Bradshaw假設(shè)修正方法(記為“SST-a1=0.355”)得到的截面壓力分布與RSM 接近,同時與實驗數(shù)據(jù)相比,修正方法在外翼段(即y/b=90%,95%及99%截面)得到了精度優(yōu)于RSM 和BSL 模式的激波模擬結(jié)果,此外激波?分離再附后壓力分布的模擬準(zhǔn)確性也與RSM 相當(dāng),且相比于BSL 模式,與實驗數(shù)據(jù)更加吻合.需要說明的是,由于Bradshaw 假設(shè)是一種基于實驗觀測的半經(jīng)驗關(guān)系式,直接對其進(jìn)行經(jīng)驗化修正是否符合物理需要進(jìn)一步考慮.

圖9 α=6.06°狀態(tài)計算壓力分布與實驗的對比(修正方法與RSM、SST 和BSL 模式對比)Fig.9 Comparison of computed pressure distribution with experiment at α=6.06° (correction methods with RSM,SST and BSL models)

圖10 表面壓力分布與極限流線對比(修正方法)Fig.10 Comparison of surface pressure distribution and limited streamline pattern (correction methods)

3.2 基于湍流長度尺度的修正方法

無論是RANS 方法還是LES、DNS 方法,尺度與能量的關(guān)系都是最重要的理論基礎(chǔ).從機(jī)理角度看,湍流的(特征)長度尺度對應(yīng)湍流邊界層的大尺度結(jié)構(gòu),支配湍流邊界層的特性;同時尺度也是湍流能量的另一種體現(xiàn),根據(jù)湍流渦拉伸和能量級聯(lián)理論,尺度越大則能量也越大.RANS 湍流建模實際都可歸結(jié)到對湍流長度尺度的建模問題上.但RANS方法并不存在解析尺度,因此在湍流模式中,長度尺度實際是一種“平均尺度”的概念.對于k-ω模式,假設(shè)ω和渦黏性系數(shù) νT都是湍動能k的函數(shù),根據(jù)量綱分析有

式中C為常數(shù).可以發(fā)現(xiàn),長度尺度l是通過湍流單位耗散率ω定義的(Wilcox 稱ω方程為“尺度決定方程”(scale determining equation)[8],即單位耗散率ω為湍流長度尺度相關(guān)量),ω方程決定了建模尺度也就決定了模式所能預(yù)測的湍動能(雷諾應(yīng)力)大小.換句話說,誰(湍流模式)能夠更精準(zhǔn)地對湍流特征長度尺度進(jìn)行建模,誰就能更準(zhǔn)確地把握湍流的主要特征,得到更理想的模擬結(jié)果.針對激波?邊界層分離模擬問題,可以考慮對ω方程進(jìn)行修正,進(jìn)而提高渦黏性系數(shù)和雷諾應(yīng)力的建模精度.很多湍流模式的問題如旋轉(zhuǎn)/曲率問題等[38]也都是以此為出發(fā)點進(jìn)行定位、分析和修正的.

針對長度尺度最著名的研究當(dāng)屬混合長度理論[7].該理論雖因未考慮流動歷史效應(yīng)已經(jīng)很少使用,但其對湍流大尺度的研究和目前的RANS 湍流/轉(zhuǎn)捩建模方法具有指導(dǎo)意義.Prandtl 類比分子動理論,認(rèn)為湍流脈動速度(類比于分子熱運(yùn)動的平均速度)應(yīng)大致取決于當(dāng)?shù)氐钠骄俣忍荻扰c流體微元跳動距離(類比于分子自由程長度)的乘積,例如對于二維平行剪切流,沿y向跳動的流體微團(tuán)具有如下形式的脈動速度

其中l(wèi)m稱為混合長度,它相當(dāng)于流體微元的拉格朗日自由程,這個距離小于平均流動變化的尺度(即在這個距離內(nèi)跳動基本不喪失其原有速度),從概念上講可以將它視為對湍流大尺度的一種度量.隨后,馮·卡門和Van Driest 相繼提出并完善了混合長度的計算方式[33],其中馮·卡門提出湍流邊界層中當(dāng)?shù)孛}動長度尺度基本與壁面距離成正比(由于壁面對于湍流脈動起限制作用),即

其中κ稱為卡門常數(shù),約等于0.41.這說明湍流大尺度主要集中在遠(yuǎn)離壁面的邊界層對數(shù)律層(對數(shù)率層中黏性效應(yīng)也基本可以忽略)和速度虧損率層,這也正是上文提到的邊界層內(nèi)受壓力梯度影響最大的區(qū)域,應(yīng)該作為修正對象.Kolmogorov[42]稱l為“external (length) scale of turbulence”,其中的“external” 指的就是邊界層外層.

根據(jù)上文的分析,原始湍流模式在三維激波-邊界層分離流動中雷諾應(yīng)力預(yù)測不足,即渦黏性系數(shù)或湍動能預(yù)測偏小,應(yīng)提高建模尺度.由于單位耗散率ω與長度尺度l成反比(式(7)),為了提高建模長度尺度應(yīng)該考慮降低ω;而降低ω最直接的方式就是增大ω方程的耗散項,這也是旋轉(zhuǎn)/曲率修正等采用的思路.根據(jù)上文的討論,對于強(qiáng)逆壓梯度下的非平衡流動,湍動能的生成與耗散不再相當(dāng),因此兩者之比是比較合適的指標(biāo)函數(shù),即

然而并不能在湍動能生成項一旦大于耗散項時就啟動修正,畢竟兩者在附著流動和弱逆壓梯度下的比值也并不是一直相等的,因此需要對上述比值設(shè)置一個開啟閾值,相當(dāng)于在逆壓梯度大于一定程度時才進(jìn)行修正.根據(jù)Menter[13]的研究,湍動能生成率與耗散率的最大比值即便在復(fù)雜剪切層中也僅在2 左右.通過對不同攻角下ONERA M6 機(jī)翼流動的測試并與RSM 對比,發(fā)現(xiàn)該閾值取為1.5 能夠基本滿足預(yù)測需求.故構(gòu)造如下形式的ω方程耗散項修正函數(shù)

縱觀三季度的冰箱線上市場，三門、多門、對開門冰箱呈現(xiàn)出不同幅度的增長，以多門勢頭最盛，漲幅高達(dá)61.9%；而線下除多門冰箱有7.5%的增長外，其他品類均出現(xiàn)不同程度的下跌?？傮w來看，多門冰箱表現(xiàn)搶眼，已成為消費者的主流選擇，且上升趨勢十分明顯；對開門冰箱上升勢頭有所放緩，但依然是不少消費者的最終選擇；三門冰箱以高性價比也獲得了越來越多的線上客戶青睞，但在對消費體驗要求更高的線下市場，已經(jīng)少有昔日風(fēng)光；兩門冰箱則逐漸式微，尤其在線下市場，消費者的目光更容易被中高端產(chǎn)品吸引。

其中α的作用類似一個放大系數(shù),即在Pk/Dk大于1.5 后控制修正量的大小.修正系數(shù)FPD與原始ω方程的耗散項相乘就得到新的耗散項

采用上述修正方法對ONERA M6 機(jī)翼α=6.06°狀態(tài)進(jìn)行模擬可以得到如圖11 所示的結(jié)果.以外翼段RSM 計算結(jié)果為參照,放大系數(shù)α經(jīng)過測試取為4 較為合適.圖12 給出了采用上述修正方法得到的y/b=90%截面湍動能生成與耗散之比云圖,可以發(fā)現(xiàn)在翼型前緣上游的大部分區(qū)域兩者之比均大于2.壓力分布及局部Pk/Dk云圖顯示(圖13),前緣上游沿流向存在較大范圍逆壓梯度,Pk/Dk很快增大至2 以上,這也是跨聲速翼型/機(jī)翼流動較為明顯的特征,同時也是非平衡性修正發(fā)揮作用的主要區(qū)域[21].由圖13 前緣附近云圖可以看出,Pk/Dk在壁面附近很快下降至1 左右;而在分離之前,從邊界層中間位置(對數(shù)律層)開始,Pk/Dk又開始增長并在分離之前達(dá)到2 以上,此時修正方法也會發(fā)揮作用.可以發(fā)現(xiàn),若將上述開啟閾值(目前為1.5)減小,則修正將在前緣及分離前更上游的位置開啟,使得ω降低、建模尺度提高的區(qū)域增大,進(jìn)而造成雷諾應(yīng)力增大,分離和激波位置將后移.與原始k-ωSST 模式模擬結(jié)果相比(圖11),加入FPD修正項后推遲了外翼段的分離,計算精度基本與RSM 相當(dāng).但相比于Bradshaw 假設(shè)修正方法(即“SST-a1=0.355”)的計算結(jié)果,FPD使得湍動能增長過度,激波和分離位置被大大推后,同時分離再附后壓力分布模擬結(jié)果并不理想.分析造成上述結(jié)果的原因:雖然采用湍動能生成項與耗散項之比能夠識別前緣附近及邊界層內(nèi)非平衡性的增長,但這仍然沒有考慮到長度尺度的實際建模不準(zhǔn)問題.也就是說,“湍動能生成項與耗散項之比大于1.5”只是依據(jù)計算結(jié)果的判斷,此時計算結(jié)果是否準(zhǔn)確、建模長度尺度是否足夠還應(yīng)該進(jìn)行判斷,以決定是否開啟FPD.這就需要對長度尺度進(jìn)行量化.

圖11 采用FPD 修正的壓力分布與其他結(jié)果的對比Fig.11 Comparison of FPD-fixed pressure distribution with other results

圖12 采用FPD 修正的截面湍動能生成/耗散之比Fig.12 Sectional contour of Pk/Dk using FPD correction

圖13 y/b=90%截面前緣及分離前湍動能生成/耗散之比Fig.13 Pk/Dk contour at leading edge and before separation of y/b=90% section

在S-A 模式[15]的建模過程中,利用了一個長度尺度的比例定義,即

其中,y為最近壁面的距離,為S-A 模式中定義的一個量,它等于渦量加上一個修正項.數(shù)值計算表明,在對數(shù)律層r≈ 1,向更外層發(fā)展這個數(shù)值開始減少,在更內(nèi)層由于 νT很小,這個值逐漸接近于0.可以從雷諾應(yīng)力的角度理解r的變化.由于壁面的阻尼作用,雷諾應(yīng)力在壁面附近(黏性底層)幾乎為0,而在黏性底層之上雷諾應(yīng)力和渦黏性系數(shù) νT將快速增大,并在對數(shù)律層某處達(dá)到最大,此時分子黏性的作用基本可以忽略.因此可以認(rèn)為對數(shù)律層所包含的長度尺度和能量是最大的,而 νT可以作為衡量尺度變化的變量.根據(jù)量綱分析,可以采用應(yīng)變不變量或渦量進(jìn)行量綱調(diào)整,若采用S-A模式定義的,就會得到一個具有長度量綱的量,即.同時上文已經(jīng)指出,混合長度lm=κy也可以作為湍流邊界層大尺度的度量,因此在對數(shù)律層附近兩者應(yīng)近似相等,即

但由于是采用計算值定義的,而混合長度lm=κy在計算網(wǎng)格確定時就已經(jīng)確定了,因此式(14) 可以理解為“計算”長度尺度(分子) 與“應(yīng)有”長度尺度(分母 κy)的比.如果建模及計算的渦黏性系數(shù)不足,則上式中的計算“長度尺度”就會小于“應(yīng)有長度尺度”,比例系數(shù)r不再等于1.因此可以利用這個關(guān)系,當(dāng)r小于某個值后才開啟式(12)中的FPD項.采用QCR 方法[35]中用于歸一化的參考量(下式根號項)代替,最終得到的比例系數(shù)為

借鑒粗糙壁面修正中Hellsten-Laine 開關(guān)函數(shù)[43]的形式構(gòu)造“長度尺度開關(guān)函數(shù)”,采用雙曲正切函數(shù)對是否開啟FPD進(jìn)行控制

其中系數(shù)β和n是類比Hellsten-Laine 開關(guān)函數(shù)[43]引入的:β的作用可以理解為當(dāng)長度尺度比例系數(shù)r小于1/β時開啟FPD(此時tanh 項接近于0,Fscale≈1),而冪指數(shù)n的作用在于使Fscale在r=1/β左右進(jìn)行0-1 快速切換.于是新的修正函數(shù)FPD為

與1 相比取最大值的操作意味著不存在耗散項的減小.以O(shè)NERA M6 機(jī)翼α=6.06°狀態(tài)外翼段RSM 計算結(jié)果為參照,通過測試將系數(shù)β和n分別取為9 和3.圖14 給出了不同修正方式下的y/b=90%截面湍流量云圖,可以發(fā)現(xiàn),經(jīng)過Fscale限制后的FPD(式(17))較未限制時(式(11))對單位耗散率的抑制作用減弱,分離位置更靠前.

圖14 對修正函數(shù)FPD 進(jìn)行限制前后的湍流量計算結(jié)果對比Fig.14 Comparison of the computed turbulence variables before and after limiting the modified function FPD

圖9 中也顯示了采用長度尺度修正方法(記為“SST-length scale”)得到的截面壓力分布,可以發(fā)現(xiàn),length scale 修正方法與Bradshaw 假設(shè)修正方法得到的結(jié)果比較接近,初步證明了所提出的兩種修正方法的有效性.

圖10 進(jìn)一步給出了兩種修正方法計算的表面流動形態(tài),其中壓力分布接近.不同的是,兩種修正方法得到的外翼段分離區(qū)較RSM 的預(yù)測結(jié)果更大(RSM 外翼段有明顯的分離再附),且length scale 修正方法在內(nèi)翼段第二道激波處也預(yù)測到了一定的分離再附流動,這些現(xiàn)象是否存在還需要實驗等方法的進(jìn)一步驗證.

4 驗證與討論

4.1 網(wǎng)格收斂性分析

上文的網(wǎng)格收斂性分析已表明,ONERA M6 機(jī)翼大攻角流動是由分離主導(dǎo)的,無法準(zhǔn)確求解是由于k-ωSST 模式未能準(zhǔn)確描述流動的主要特性,這一點無法通過網(wǎng)格加密或提高格式精度彌補(bǔ).采用所提出的兩種修正方法后(圖15),外翼段表現(xiàn)出了合理的網(wǎng)格收斂特性(即隨著網(wǎng)格加密,激波分辨率提高),表明兩種修正方法均能正確反映該問題的流動特征,截斷誤差能夠隨網(wǎng)格加密不斷減小,這進(jìn)一步驗證了兩種修正方法的合理性.從圖15 中還可以看出,length scale 修正方法對表面網(wǎng)格密度的敏感性較低,稀網(wǎng)格規(guī)模就能基本滿足模擬精度要求.

圖15 α=6.06°狀態(tài)兩種修正方法的網(wǎng)格收斂性對比Fig.15 Grid convergence comparison of two corrections at α=6.06°

4.2 其他攻角狀態(tài)驗證

為了進(jìn)一步確認(rèn)兩種修正模擬方法的正確性,對ONERA M6 機(jī)翼另外兩個攻角狀態(tài)進(jìn)行數(shù)值模擬,典型截面的計算結(jié)果如圖16 所示.可以看出,α=4.08°狀態(tài)各個模式和修正方法的模擬結(jié)果接近;而對于α=5.06°狀態(tài),兩種修正模擬方法的計算結(jié)果與RSM 基本相同,相比于原始模式改進(jìn)明顯.

圖16 兩種修正方法應(yīng)用于其他攻角狀態(tài)的計算壓力分布對比Fig.16 Computed pressure distribution comparison of two corrections at other angles of attack

圖16 兩種修正方法應(yīng)用于其他攻角狀態(tài)的計算壓力分布對比(續(xù))Fig.16 Computed pressure distribution comparison of two corrections at other angles of attack (continued)

4.3 雷諾應(yīng)力特性分析

下面分析三維激波?邊界層分離流動的雷諾應(yīng)力特性.圖17(a)為采用DNS 方法計算的三維零壓力梯度湍流平板邊界層雷諾應(yīng)力分布,其中u,v,w分別代表流向、法向和展向速度.可以發(fā)現(xiàn),近壁區(qū)雷諾應(yīng)力分量之間存在比較明顯的差異,其中流向雷諾正應(yīng)力和流向?法向雷諾切應(yīng)力明顯大于其余應(yīng)力分量,這時主雷諾應(yīng)力這一概念是成立的.圖17(b)為ONERA M6 機(jī)翼α=6.06°狀態(tài)y/b=90%站位前緣x/c=8%處邊界層內(nèi)的雷諾應(yīng)力分布.由于原始SST 模式預(yù)測的分離非常靠前,已經(jīng)覆蓋了該位置,因此雷諾應(yīng)力分量都很大(點劃線).以Bradshaw 假設(shè)修正方法(虛線)為例,與RSM (實線)計算的雷諾應(yīng)力分布對比可以發(fā)現(xiàn),法向和展向雷諾正應(yīng)力的相對大小較平板邊界層明顯增大,特別是展向雷諾正應(yīng)力增長明顯;同時,除流向?法向切應(yīng)力外,其余兩個切應(yīng)力分量也增長明顯,三者基本具有相同量級,均不能被忽略.這也證明了在三維激波?邊界層分離流動中,主雷諾應(yīng)力分量并不唯一,對各向雷諾正應(yīng)力和切應(yīng)力都需要進(jìn)行比較準(zhǔn)確(量級至少要正確)的評估才能最終得到正確的分離和壓力分布計算結(jié)果.

從圖17 中還可以發(fā)現(xiàn),Bradshaw 假設(shè)修正方法得到的雷諾應(yīng)力分量之間的相對大小與RSM 的計算結(jié)果是一致的,但絕對大小還是有區(qū)別的.整體來看,與RSM 相比,渦黏性模式修正方法計算的雷諾正應(yīng)力較大,切應(yīng)力較小.

圖17 不同流動中雷諾應(yīng)力分量的對比Fig.17 Comparison of Reynolds-stress component for different flows

4.4 湍流平板邊界層壁面律驗證

最后,對湍流平板邊界層進(jìn)行確認(rèn)計算,驗證上述修正方法是否會對附著流動模擬造成影響.參考Wieghardt 等[44]湍流平板實驗的參數(shù)和數(shù)據(jù),計算網(wǎng)格及邊界條件如圖18 所示.圖19 給出了不同方法得到的壁面律曲線與實驗結(jié)果及光滑平板湍流邊界層壁面律理論曲線的對比.可以發(fā)現(xiàn),原始k-ωSST 模式及兩種修正方法在線性底層和對數(shù)律層均得到了一致且與實驗和理論曲線吻合的計算結(jié)果,說明修正方法并未對模式原有的附著流動模擬能力造成不利影響,具有一定的通用性.

圖18 湍流平板邊界層計算網(wǎng)格Fig.18 Computational grid of turbulent boundary layer of flat plate

圖19 不同方法的壁面律曲線與實驗數(shù)據(jù)和理論值的對比Fig.19 Computational wall-function law by different methods comparing with experimental and theoretical data

5 結(jié)論

針對常用渦黏性湍流模式在三維激波分離流動中的模擬能力不足,基于k-ωSST 模式提出了兩種修正方法,均取得了滿意的改進(jìn)效果,為渦黏性湍流模式應(yīng)用于工程中的此類流動提供了思路.本文主要結(jié)論如下.

(1)渦黏性湍流模式中的渦黏性系數(shù)定義方式不適用于非平衡流動,k-ωSST 模式為此引入的Bradshaw 假設(shè)考慮了雷諾應(yīng)力的輸運(yùn)特性,使模式具備了較好的二維激波分離流動模擬能力.

(2)三維激波分離流動不存在主雷諾應(yīng)力分量(存在多個量級相當(dāng)?shù)睦字Z應(yīng)力分量),使得Bradshaw假設(shè)限制了雷諾應(yīng)力的生成從而失效.同時,目前已有的非線性雷諾應(yīng)力本構(gòu)關(guān)系并不能從本質(zhì)上提高模擬能力,原因是非線性修正主要針對旋轉(zhuǎn)、剪切引起的二次效應(yīng).

(3)基于提高Bradshaw 常數(shù)和基于長度尺度的修正方法均使k-ωSST 模式表現(xiàn)出了合理的網(wǎng)格收斂特性(其中后者的網(wǎng)格敏感性較低),對于ONERA M6 機(jī)翼大攻角狀態(tài)和其他攻角狀態(tài)也能夠進(jìn)行準(zhǔn)確模擬,驗證了本文修正方法的合理性和有效性.

由于本文針對三維激波分離流動提出的修正方法與k-ωSST 模式引入Bradshaw 假設(shè)的初衷相反,因此兩種修正方法應(yīng)用于二維激波分離流動會得到稍向后移的激波分離位置.這涉及到方法的二維/三維適用性問題,需要對修正函數(shù)等進(jìn)行更細(xì)節(jié)的考慮.此外,由于長度尺度修正方法的模式化、經(jīng)驗化意味更強(qiáng),本文所給定的修正參數(shù)的通用性還有待進(jìn)一步驗證和改進(jìn).最后,若處理更高馬赫數(shù)的問題,壓縮性修正也是需要考慮的問題.