B-C轉捩模型的標定與應用*

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

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

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

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

1 B-C轉捩模型簡介

非結構軟件平臺中的B-C轉捩模型與文獻[13]中的略有不同，修正了部分量綱不匹配的問題。間歇因子γ是轉捩過程中流動狀態的一個物理性質，整個轉捩過程γ從層流狀態時的0變為全湍流狀態時的1。B-C轉捩模型的思想就是建立一個間歇因子γ的表征函數用于實現轉捩的數值模擬。S-A湍流模型的方程為：

(1)

B-C模型在原湍流模型生成項上乘上γ函數，從而達到在湍流計算中引入轉捩的目的。

(2)

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

(3)

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

(4)

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

(5)

目前的B-C轉捩模型使用的轉捩動量厚度雷諾數表達式是一個與湍流強度相關的經驗關聯式[3]，即為：

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

(6)

2 B-C轉捩模型標定過程

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

關聯函數的標定一般選取低速平板實驗中的S&K、T3A、T3A-和T3B。其中低湍流度的自然轉捩有S&K和T3A-實驗，高湍流度的旁路轉捩有T3A和T3B實驗。

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

圖1 平板計算網格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

關聯函數Reθt的標定工作具體步驟如下所示：

1) 采用原有的Reθt經驗關聯式進行計算、并記錄。將計算結果與實驗值進行比較，若誤差很小，則不再進行計算。

2) 對于步驟1)中誤差較大的算例，需將Reθt定為常數代入計算。如果與實驗值誤差較大，則換另一常數代入計算。如此循環，直到計算值與實驗值之間的誤差很小為止。記錄下最后的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中的數據得到新的經驗關聯式，如式(7)所示：

Reθt=-4.351Tu3+80.49Tu2-

502.705Tu+1 162

Tu≥0.027;Reθt≥20

(7)

原經驗關聯式曲線以及新經驗關聯式曲線，如圖2所示。

圖2 經驗關聯式與計算值的對比Fig.2 The comparison of empirical correlation and calculated value

標定后的B-C轉捩模型的低速平板算例的計算結果如圖3-6所示。

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

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

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

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

3 標定后的B-C轉捩模型在二維低速問題上的應用

3.1 低雷諾數算例

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

圖7 E387翼型算例的網格圖Fig.7 Mesh of E387 airfoil

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

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

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

3.2 中高雷諾數算例

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

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

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

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

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

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

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

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

4 結 論

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

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

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