B-C轉(zhuǎn)捩模型的標(biāo)定與應(yīng)用*

余秋陽，鄧小剛，包蕓，王光學(xué)，王靖宇，張懷寶

(1.中山大學(xué)工學(xué)院，廣東 廣州 510006；2.國防科技大學(xué)空天科學(xué)學(xué)院，湖南 長沙 410073；3.中山大學(xué)物理學(xué)院，廣東 廣州 510006)

轉(zhuǎn)捩現(xiàn)象的研究對飛行器設(shè)計非常重要。層流邊界層和湍流邊界層結(jié)構(gòu)的不同會導(dǎo)致很多氣動力參數(shù)在轉(zhuǎn)捩后發(fā)生很大的變化。20世紀(jì)末，湍流模式理論顯示出它在工程應(yīng)用上的優(yōu)勢，很多模型如S-A模型[1]、k-ω模型[2]等涌現(xiàn)出來，全湍流的計算也得到了很多令人振奮的結(jié)果。但與此同時，轉(zhuǎn)捩現(xiàn)象的數(shù)值模擬卻沒有得到很好的發(fā)展，原因可以歸結(jié)為轉(zhuǎn)捩現(xiàn)象本身的復(fù)雜性[3]：① 在于物理過程的復(fù)雜性。轉(zhuǎn)捩的類型包括自然轉(zhuǎn)捩、旁路轉(zhuǎn)捩、分離誘導(dǎo)轉(zhuǎn)捩等，同時湍流邊界層在強順壓梯度下還可能再層流化，很難從一兩個物理量的變化得到全部的物理信息。② 在于數(shù)學(xué)表達(dá)的復(fù)雜性。轉(zhuǎn)捩流動中線性和非線性效應(yīng)的相互影響導(dǎo)致幾乎不可能用單一數(shù)學(xué)模型來刻畫這一物理現(xiàn)象。

C轉(zhuǎn)捩模型由Onur Bas和Samet等在2013年提出[13]。它與γ-Reθ模型一樣基于當(dāng)?shù)仃P(guān)聯(lián)。B-C模型使用函數(shù)公式而不是輸運方程去激發(fā)轉(zhuǎn)捩過程。它耦合于S-A湍流模型上使用，因其同樣包含經(jīng)驗關(guān)聯(lián)式，也可以根據(jù)不同的轉(zhuǎn)捩關(guān)聯(lián)進行標(biāo)定。本文詳細(xì)介紹了B-C轉(zhuǎn)捩模型，利用零壓力梯度平板轉(zhuǎn)捩實驗數(shù)據(jù)得到一組關(guān)聯(lián)函數(shù)；并針對E387、NLF(1)-416翼型做了驗證計算，將得到的結(jié)果與風(fēng)洞實驗結(jié)果、全湍流計算結(jié)果進行了比較，獲得了有益的結(jié)論。

1 B-C轉(zhuǎn)捩模型簡介

非結(jié)構(gòu)軟件平臺中的B-C轉(zhuǎn)捩模型與文獻[13]中的略有不同，修正了部分量綱不匹配的問題。間歇因子γ是轉(zhuǎn)捩過程中流動狀態(tài)的一個物理性質(zhì)，整個轉(zhuǎn)捩過程γ從層流狀態(tài)時的0變?yōu)槿牧鳡顟B(tài)時的1。B-C轉(zhuǎn)捩模型的思想就是建立一個間歇因子γ的表征函數(shù)用于實現(xiàn)轉(zhuǎn)捩的數(shù)值模擬。S-A湍流模型的方程為：

(1)

B-C模型在原湍流模型生成項上乘上γ函數(shù)，從而達(dá)到在湍流計算中引入轉(zhuǎn)捩的目的。

(2)

這樣的轉(zhuǎn)捩模型不會改變計算過程中的方程個數(shù)，形式上也較簡單，這是它優(yōu)于輸運方程轉(zhuǎn)捩模型之處。其中，γ函數(shù)的具體形式為：

(3)

(3)式由兩部分組成：第一部分將渦雷諾數(shù)與轉(zhuǎn)捩起始動量厚度雷諾數(shù)關(guān)聯(lián)起來，起到判斷轉(zhuǎn)捩位置的作用；第二部分在近壁面區(qū)起作用。近壁面區(qū)渦雷諾數(shù)等于0，以至于第二部分也會等于0。第一部分使用的關(guān)聯(lián)方法來源于文獻[4]，其具體形式為：

(4)

(3)式中參數(shù)δ1=0.002，δ2=5.0，且

(5)

目前的B-C轉(zhuǎn)捩模型使用的轉(zhuǎn)捩動量厚度雷諾數(shù)表達(dá)式是一個與湍流強度相關(guān)的經(jīng)驗關(guān)聯(lián)式[3]，即為：

Reθt=803.73(Tu+0.606 7)-1.027

(6)

2 B-C轉(zhuǎn)捩模型標(biāo)定過程

B-C轉(zhuǎn)捩模型的標(biāo)定是指利用T3系列低速平板實驗的實驗數(shù)據(jù)對模型中的經(jīng)驗關(guān)聯(lián)式進行標(biāo)定。B-C模型中的經(jīng)驗關(guān)聯(lián)式只有Reθt的經(jīng)驗關(guān)聯(lián)，原Reθt經(jīng)驗關(guān)聯(lián)在B-C模型中不能很好地模擬轉(zhuǎn)捩位置，與實驗結(jié)果有所偏差，而，工程上更關(guān)心的是轉(zhuǎn)捩位置，因此本文的標(biāo)定工作主要依據(jù)轉(zhuǎn)捩位置。本文采用與原有Reθt類似的樣條擬合曲線，描述對Reθt的標(biāo)定。

關(guān)聯(lián)函數(shù)的標(biāo)定一般選取低速平板實驗中的S&K、T3A、T3A-和T3B。其中低湍流度的自然轉(zhuǎn)捩有S&K和T3A-實驗，高湍流度的旁路轉(zhuǎn)捩有T3A和T3B實驗。

標(biāo)定過程中使用的網(wǎng)格如圖1所示。網(wǎng)格數(shù)量為348×119，共41 412個矩形單元。網(wǎng)格在壁面處加密以保證邊界層內(nèi)的正確模擬。另外，在x=0處也加密了網(wǎng)格。標(biāo)定計算在非結(jié)構(gòu)軟件平臺上進行，對流項離散采用二階精度的Roe格式，并使用預(yù)處理技術(shù)以適應(yīng)低速流動。低速平板實驗來流條件見表1, 表中FSTI為來流湍流強度。

圖1 平板計算網(wǎng)格Fig.1 Computational mesh for the flat plate

算例U/(m·s-1)μ/(Pa·s)ρ/(kg·m-3)FSTI/%S&K50.11.8×10-51.20.18T3A-19.81.8×10-51.20.874T3A5.41.8×10-51.23.5T3B9.41.8×10-51.26.5

關(guān)聯(lián)函數(shù)Reθt的標(biāo)定工作具體步驟如下所示：

1) 采用原有的Reθt經(jīng)驗關(guān)聯(lián)式進行計算、并記錄。將計算結(jié)果與實驗值進行比較，若誤差很小，則不再進行計算。

2) 對于步驟1)中誤差較大的算例，需將Reθt定為常數(shù)代入計算。如果與實驗值誤差較大，則換另一常數(shù)代入計算。如此循環(huán)，直到計算值與實驗值之間的誤差很小為止。記錄下最后的Reθt值。

3)記錄步驟1)和步驟2)的所有Reθt值。新Reθt值與原Reθt值比較如表2所示。

表2 平板算例的Tu和ReθtTable 2 Tu and Reθt of the flat plate

擬合表2中的數(shù)據(jù)得到新的經(jīng)驗關(guān)聯(lián)式，如式(7)所示：

Reθt=-4.351Tu3+80.49Tu2-

502.705Tu+1 162

Tu≥0.027;Reθt≥20

(7)

原經(jīng)驗關(guān)聯(lián)式曲線以及新經(jīng)驗關(guān)聯(lián)式曲線，如圖2所示。

圖2 經(jīng)驗關(guān)聯(lián)式與計算值的對比Fig.2 The comparison of empirical correlation and calculated value

標(biāo)定后的B-C轉(zhuǎn)捩模型的低速平板算例的計算結(jié)果如圖3-6所示。

圖3 S&K平板算例的計算結(jié)果Fig.3 Calculation result of S&K flat plate

圖4 T3A平板算例的計算結(jié)果Fig.4 Calculation result of T3A flat plate

圖5 T3A-平板算例的計算結(jié)果Fig.5 Calculation result of T3A- flat plate

圖6 T3B平板算例的計算結(jié)果Fig.6 Calculation result of T3B flat plate

3 標(biāo)定后的B-C轉(zhuǎn)捩模型在二維低速問題上的應(yīng)用

3.1 低雷諾數(shù)算例

Eppler E387翼型的轉(zhuǎn)捩仿真是典型的低雷諾數(shù)翼型轉(zhuǎn)捩算例。低雷諾數(shù)下的翼型繞流問題十分復(fù)雜，通常會有分離泡[14]，影響轉(zhuǎn)捩的因素也較多。本文選取E387翼型在雷諾數(shù)為2×105和3×105時的流動，對本文的標(biāo)定模型進行驗證計算。算例采用O網(wǎng)格構(gòu)型，共129 870個單元，如圖7所示。網(wǎng)格在翼面附近加密，第一層網(wǎng)格控制在10-6量級以保證y+<1。

圖7 E387翼型算例的網(wǎng)格圖Fig.7 Mesh of E387 airfoil

本算例來流湍流強度FSTI=0.1%，對流項離散采用二階精度的Roe格式，時間離散采用歐拉隱式格式，離散方程采用LU-SGS方法求解。為了更好的模擬低速情況，加快算法的收斂，本文采用了預(yù)處理技術(shù)。本文根據(jù)氣動力參數(shù)的計算結(jié)果，繪制了升力系數(shù)曲線和極曲線，分析了升阻力系數(shù)的變化趨勢。最后，將轉(zhuǎn)捩模型的計算結(jié)果與全湍流模型計算結(jié)果、實驗結(jié)果進行了對比分析，實驗結(jié)果取自文獻[13]。當(dāng)計算攻角α從-3°增大到12°時，Re分別為2×105和3×105的升力系數(shù)、極曲線如圖8-9所示。

從圖中可以看出，在不同雷諾數(shù)下，轉(zhuǎn)捩模型的計算結(jié)果都明顯優(yōu)于全湍流模型。全湍流模型可以很好地模擬升力系數(shù)，但阻力系數(shù)卻與實驗結(jié)果相差較大，可見轉(zhuǎn)捩現(xiàn)象的影響在氣動力參數(shù)計算上的重要性。另外，還可以發(fā)現(xiàn)不管是升力變化還是極曲線的計算結(jié)果與實驗結(jié)果都很吻合。隨著攻角的不斷增大，超過失速角后，可能是由于流動分離過大，轉(zhuǎn)捩模型的計算結(jié)果會出現(xiàn)偏差。

圖9 Re=3×105時的升力系數(shù)和極曲線Fig.9 Lift coefficient and polar curve at Re=3×105

3.2 中高雷諾數(shù)算例

NASA/LANGLEY NLF(1)-416翼型被廣泛應(yīng)用于轉(zhuǎn)捩模型的測試計算[13、15]。NLF(1)-416算例的實驗結(jié)果取自文獻[13]。本文計算了Re=1×106和2×106兩種工況，來流湍流強度為0.03%。本算例的網(wǎng)格采用C型結(jié)構(gòu)網(wǎng)格，共189 945個單元，如圖10所示。翼面上分布了1 048個網(wǎng)格點，網(wǎng)格在翼面附近和尾跡區(qū)加密，翼面上第一層網(wǎng)格高度設(shè)為1×10-5以保證y+<1。

圖10 NLF(1)-416翼型算例的網(wǎng)格圖Fig.10 Mesh of NLF(1)-416 airfoil

3.2.1 氣動力參數(shù) 本算例采用與E387翼型一樣的計算方法與離散格式。當(dāng)計算攻角α從-6°增大到12°時，Re分別為1×106和2×106的升力系數(shù)、極曲線如圖11-12所示。

從圖11-12可以看出， 對于在不同雷諾數(shù)工況下的NLF(1)-416翼型來說，轉(zhuǎn)捩模型的計算結(jié)果同樣都明顯優(yōu)于全湍流模型，更加符合實際流動。從升力系數(shù)和極曲線的變化趨勢來看，轉(zhuǎn)捩模型的計算結(jié)果與實驗結(jié)果大致相同。

3.2.2 摩阻系數(shù)分布 6°攻角下摩阻系數(shù)計算結(jié)果如圖13-14所示。

圖13-14將轉(zhuǎn)捩模型的計算結(jié)果與全湍流模型進行了對比。根據(jù)摩阻系數(shù)變化趨勢可以判斷流動是否發(fā)生了層流分離。當(dāng)摩阻系數(shù)沿弦向先降為負(fù)值時，說明出現(xiàn)了回流即發(fā)生層流分離；然后，摩阻系數(shù)迅速增長為正值，說明分離導(dǎo)致轉(zhuǎn)捩并且湍流又附著到物面上。從圖13中可以看出：在Re=1×106的工況下，上翼面在翼型弦長32.7%處發(fā)生了分離，在40%處再附著；下翼面在在翼型弦長61%處發(fā)生分離，在68.5%處再附著。從圖14處可以看出：在Re=2×106的工況下，上翼面沒有發(fā)生流動分離現(xiàn)象，下翼面分離點位于翼型弦長61%處，再附著點位于65.7%處。由摩阻系數(shù)計算結(jié)果可以看出轉(zhuǎn)捩模型的計算結(jié)果可以體現(xiàn)層流分離泡在中高雷諾數(shù)翼型轉(zhuǎn)捩過程中的影響。

圖11 Re=1×106時的升力系數(shù)和極曲線Fig.11 Lift coefficient and polar curve at Re=1×106

圖12 Re=2×106時的升力系數(shù)和極曲線Fig.12 Lift coefficient and polar curve at Re=2×106

4 結(jié) 論

本文確定了一組經(jīng)驗關(guān)聯(lián)函數(shù)，并且在低速平板算例中的摩阻系數(shù)分布與實驗結(jié)果吻合得很好，說明了本文的標(biāo)定工作是成功的；在二維低速問題中對比了標(biāo)定后的B-C轉(zhuǎn)捩模型和全湍流模型，選擇了中低雷諾數(shù)兩個算例。計算結(jié)果表明標(biāo)定后的B-C轉(zhuǎn)捩模型可以很好地捕捉轉(zhuǎn)捩過程，中等攻角下計算出的升力阻力系數(shù)與實驗結(jié)果吻合得較好。標(biāo)定后的B-C轉(zhuǎn)捩模型還可以反映層流分離泡在中高雷諾數(shù)翼型轉(zhuǎn)捩過程中的影響。這些都說明了B-C轉(zhuǎn)捩模型的合理性。

圖13 Re=1×106時的摩阻系數(shù)計算結(jié)果Fig.13 Frictions coefficient around the airfoil at Re=1×106

圖14 Re=2×106時摩阻系數(shù)計算結(jié)果Fig.14 Frictions coefficient around the airfoil at Re=1×106