適用于可壓縮流動的 -Re t-fRe轉捩模型

劉清揚，雷娟棉，劉周，石磊，周偉江

1. 北京理工大學 宇航學院，北京 100081 2. 中國航天空氣動力技術研究院，北京 100074

高超聲速飛行器研制以及與之相關的高超聲速空氣動力學是當前研究的重點，其中高超聲速邊界層轉捩研究具有十分迫切的工程需求。邊界層轉捩對高超聲速飛行器的氣動特性影響顯著，而低冗余、一體化的高超聲速飛行器設計需要較為準確的預測轉捩位置和轉捩區(qū)的發(fā)展，這不僅要求高超聲速轉捩理論的進步和實驗能力的提高，同時也需要工程實用化的高超聲速轉捩數(shù)值預測手段。

最近的十幾年里，基于局部流場變量構造，通過求解偏微分輸運方程的轉捩模型得到廣泛的研究和應用。其中最具有代表性的是- 轉捩模型，模型構造了基于間歇因子和轉捩動量雷諾數(shù) 的輸運方程，并通過大量實驗數(shù)據(jù)分析建立的轉捩區(qū)長度與臨界轉捩雷諾數(shù)的經(jīng)驗關系式，使得該模型對轉捩預測有了一定的可靠性，因此被很多CFD代碼和軟件采用，成為當前工程研究最流行的轉捩模型。

標準- 轉捩模型基于低速流動發(fā)展得到，因此將- 轉捩模型應用于超聲速和高超聲速流動時需要對模型進行修正。張毅鋒等針對可壓縮性和湍流普朗特數(shù)進行了修正，并成功驗證了MF-1飛行試驗數(shù)據(jù)。Krause、張曉東、夏陳超、Frauholz等各自發(fā)展了針對高超聲速流動的關聯(lián)函數(shù)。Wang等利用相似解改進了- 轉捩模型在高超聲速下的轉捩觸發(fā)判據(jù)。這些針對性的改進使得- 轉捩模型的應用范圍得到進一步的擴展。

除了對- 轉捩模型進行適應高速流動的擴展外，針對高超聲速流動的特點也發(fā)展了特有的轉捩模型。王亮和符松基于非湍流脈動概念提出了適用于高超聲速流動的--轉捩模型。周玲等改進了--轉捩模型，使之擴展到高超聲速橫流轉捩模擬。徐家寬等構建了-轉捩模型，其中是層流脈動渦黏性系數(shù)，較好地預測了由第一和第二Mack不穩(wěn)定模態(tài)誘發(fā)的超聲速和高超聲速轉捩。Shi等基于高超聲速靜風洞數(shù)據(jù)發(fā)展了工程轉捩模型，較好地預測了HiFiRE-5外形的高超聲速轉捩。

這些面向超聲速和高超聲速流動所修正和發(fā)展的轉捩模型已經(jīng)取得一定的成功。與其他模型相比，- 轉捩模型可以通過與試驗的關聯(lián)包括多種轉捩效應，具有高度的靈活性和可擴展性，因而在高速轉捩模擬中獲得更多的應用。但高速可壓縮流動轉捩的實驗數(shù)據(jù)比低速更為稀缺，且影響高速流動轉捩的因素更多，難以從數(shù)據(jù)中獲得一個相對統(tǒng)一的轉捩準則和關聯(lián)函數(shù)，已有的修正和擴展往往針對特定條件，難以適應到更廣泛的情況，而且與原始低速模型不能很好的相容。

本文研究基于廣泛應用的- 轉捩模型，針對已有的轉捩準則引入可壓縮性修正，利用基于參考溫度法的雷諾數(shù)可壓縮性比擬關系修正現(xiàn)有基于不可壓縮流動的轉捩關聯(lián)函數(shù)，并通過額外構造的輸運方程避免使用全局參數(shù)，從而建立考慮流動可壓縮性的- -轉捩模型，以較小的計算代價實現(xiàn)- 轉捩模型向可壓縮高速流動的擴展，并具備從低速流動到高速流動的無縫統(tǒng)一模擬能力。

1 γ-Reθ t-fRe轉捩模型發(fā)展

1.1 轉捩準則的可壓縮性修正

在- 轉捩模型中，當流動的動量厚度雷諾數(shù)大于轉捩臨界雷諾數(shù) 后，轉捩被觸發(fā)。對不可壓絕熱平板邊界層，邊界層中渦量雷諾數(shù)和動量厚度雷諾數(shù)有如下關系：

max()=2193·

(1)

式中：和的定義為

=，=

(2)

其中：是到壁面的距離；為密度；為黏性系數(shù)；是應變率；是速度；是動量厚度；下標e代表邊界層的外緣。

因此- 轉捩模型的轉捩觸發(fā)控制函數(shù)為

=(2193 )

(3)

即當邊界層內(nèi)的超過2.193倍的轉捩臨界雷諾數(shù) 后，>1，此時轉捩被觸發(fā)。該條件基于不可壓絕熱平板邊界層得到，并不適應于可壓縮流動，需要修改。對可壓縮平板邊界層，和之間的關系隨邊界層外緣馬赫數(shù)、溫度以及壁面溫度而變化。圖1(a)給出了不同、相同和下可壓縮平板邊界層內(nèi)和的關系曲線；圖1(b)給出了馬赫數(shù)=6時不同壁溫下可壓縮平板邊界內(nèi)的和關系曲線。圖中為邊界層厚度。

從圖1中可見，在亞聲速時，當壁溫和邊界層邊緣溫度相等，與的比值基本和不可壓絕熱平板的值相等；在超聲速和高超聲速時，與的比值要大于不可壓平板的值2.193；在高馬赫數(shù)下，壁溫越高，2.193·和max()的偏差越大。在真實的高超聲速飛行中，壁溫顯著高于來流溫度，因此修正是必須的。對可壓縮流動，邊界層內(nèi)和關系可以修正為

max()=·2193·

(4)

圖1 可壓縮平板邊界層內(nèi)的Rev和Reθ關系曲線Fig.1 Relations between Rev and Reθ in compressible plate boundary layer

式中：是修正函數(shù)，本文采用的表達形式為

lg()=·(lg())+

·lg()+

(5)

其中：系數(shù)、和是的函數(shù)，

(6)

其中：=min(max(0.4,), 8)。

如圖1(a)所示，當邊界層邊緣馬赫數(shù)小于0.4 時，邊界層內(nèi)和關系基本不變，通過該限制確保擬合函數(shù)的取值合理。

式(5)中的是參考溫度，使用Eckert用于層流超聲速流動的參考溫度定義：

=050×+022×+028×

(7)

式中：是絕熱壁面溫度，定義為

(8)

其中：是恢復因子，空氣通常取0.85。對等溫壁面邊界，為給定值；對絕熱壁面邊界，=。

從修正函數(shù)的構造過程可見，該修正函數(shù)同時考慮了馬赫數(shù)和壁溫的影響。圖2給出了不同和不同/下修正函數(shù)的曲線，從圖中可見，隨著變小，在/逐漸趨近于1時，的值趨近于1，恢復為不可壓絕熱壁面平板邊界層流動的結果。而對高超聲速飛行器繞流，/通常在2～5之間，此時修正函數(shù)的值在1.5～2.5之間。

圖2 不同Mae下修正函數(shù)Cfc隨Tw/Te的變化Fig.2 Curves of Cfc vs Tw/Te with different Mae

在獲得修正函數(shù)后，需要將其應用至轉捩判斷中。一種直接的思路就是在式(3)的分母上乘以的值，然而的計算需要用到邊界層外緣的和，獲得這些值對于當前基于非結構網(wǎng)格的并行CFD代碼是十分困難的。

結合- 轉捩模型的構造思想，的值并不顯式的應用于式(3)的觸發(fā)判斷，而是通過修改轉捩雷諾數(shù) 來實現(xiàn)。- 轉捩模型原始采用基于實驗的轉捩準則和基于湍流度和壓力梯度參數(shù)的函數(shù)，可表示為

=()·()

(9)

針對可壓縮流動，應用修正函數(shù)，式(9)被修改為

=()·()·

(10)

通過對轉捩準則 加入可壓縮性修正，基于不可壓流動的轉捩準則被直接擴展到可壓縮流動，后續(xù)的算例計算表明這樣的擴展是合理的。

1.2 雷諾數(shù)可壓縮性比擬關系fRe輸運方程

在基本- 轉捩模型應用于高超聲速轉捩的模擬中，一個常見的問題是計算得到的轉捩區(qū)發(fā)展太慢，這表明間歇因子方程中的生成源項太小。一種解決方式是修改關聯(lián)函數(shù)以和實驗結果匹配，但這種方法不具備普適性，僅適用于要求的流動條件。

參考溫度法是平板可壓縮流動和不可壓縮流動之間的橋梁。通過參考溫度法，可壓縮層流平板摩阻和熱流可以和不可壓縮層流平板摩阻和熱流聯(lián)系在一起，兩者之間的動量厚度雷諾數(shù)也有相應的比擬關系：

(11)

(12)

(13)

圖3給出了=300 K時不同/下的曲線。從圖中可見，對實際高超聲速飛行器繞流，的值大約在4～16之間；而對不可壓流動，的值恢復為1，關聯(lián)函數(shù)恢復為原始形式。

圖3 不同Mae下雷諾數(shù)比擬關系fRe隨Tw/Te的變化Fig.3 Relations of Reynolds number relation fRe and Tw/Te with different Mae

圖4 關聯(lián)函數(shù)Flength隨變化的曲線Fig.4 Curve of Flength changing with

(14)

源項的表達式為

(15)

從式(15)可見，在邊界層內(nèi)部，源項為0，因此邊界層內(nèi)部的均由邊界層邊緣的值輸運進入，采用當?shù)亓鲃幼兞揩@得的不會對真實的造成影響。

基本- 轉捩模型的 定義為

(16)

式中：是到壁面的距離；是計算獲得的邊界層厚度；函數(shù)用于確保 函數(shù)不會在尾跡區(qū)被激活。相關函數(shù)定義和參數(shù)的值可見文獻[4]。式中()的計算是 成功的關鍵。基本- 轉捩模型采用如下方式計算：

(17)

式中：是渦量。從基本- 轉捩模型所用到的定義來看，邊界層外緣的厚度是通過不可壓平板邊界層的理論解近似獲得，對于可壓縮邊界層，這些函數(shù)關系需要修正。

雖然可以構造一個更復雜的函數(shù)來修正上述關系，但實踐發(fā)現(xiàn)，將原始定義中的=7.5簡單修正為

(18)

即可很好地改善可壓縮流動時邊界層內(nèi)外的判斷。

圖5給出了=8、=300 K、=1 500 K，采用不同邊界層內(nèi)外判斷定義時平板層流的()值的比較。從圖中可見，修正的定義更

圖5 不同δBL定義獲得的平板層流(y/δcomp)4Fig.5 (y/δcomp)4 of laminar flow on plate with different δBL

保守，與原始定義相比，修正的定義確保了在邊界層內(nèi)()=0，從而使得 =1。

至此，通過額外的輸運方程，將基本- 轉捩模型擴展為可壓縮- -轉捩模型，使之可用于超聲速和高超聲速轉捩的模擬；不可壓縮流動時比擬關系退化為1，可壓縮- -轉捩模型自然恢復為基本- 轉捩模型，從而實現(xiàn)了從不可壓低速流動至可壓縮高速流動的統(tǒng)一轉捩計算。

1.3 初始條件和壁面邊界條件

和基本- 轉捩模型一樣，- -轉捩模型中的3個輸運方程在壁面處的法向通量為0。

2 γ-Reθ t-fRe轉捩模型驗證

2.1 低速平板轉捩流動

- -轉捩模型在低速不可壓時會恢復為基本的- 轉捩模型。本節(jié)針對不同低速平板轉捩流動，驗證- -轉捩模型與基本- 轉捩模型的兼容性。

表1給出了所采用的零壓力梯度平板流動的計算條件，其中下標inlet表示距平板前緣0.04 m入口處的值。圖6給出了不同模型計算得到的摩阻系數(shù)和實驗的對比，其中湍流結果采用SST湍流模型計算得到。

表1 低速平板轉捩流動計算條件[4]Table 1 Transition flow conditions for low speed plate[4]

從圖6中可見，對于不可壓低速流動條件，- -模型的結果和基本- 模型的結果幾乎完全一致，而且和實驗吻合得相當好，這表明在低速流動條件下，- -模型自動恢復為基本- 模型，兩者在低速是完全兼容的。這同時也表明，在式(6)中采用0.4作為的下限進行截斷是恰當?shù)摹?/p>

圖6 低速平板轉捩流動的計算摩阻和實驗對比Fig.6 Friction for low speed plate transition from computation and experiment

轉捩模擬沒有很好地獲得T3A平板在轉捩結束時的過沖現(xiàn)象，這可能和當前采用的基本湍流模型有關。本文與- 轉捩模型進行耦合的是兩方程SST湍流模型，文獻[24]將- 轉捩模型與雷諾應力湍流模型進行耦合，較好地模擬了T3A平板在轉捩結束處的過沖現(xiàn)象。

2.2 高超聲速平板轉捩流動

高超聲速平板轉捩流動對應的實驗是Mee在激波風洞中進行的，本文選取的計算條件見表2，其中下標inlet表示距平板前緣0.04 m入口處的值。壁面邊界條件為等溫壁，壁溫=300 K。

表2 高超聲速平板流動條件Table 2 Flow conditions for hypersonic speed plate

計算網(wǎng)格約13萬網(wǎng)格單元，為確保收斂和正確地捕獲轉捩位置，壁面均小于1。圖7給出了不同條件下計算獲得的壁面Stanton數(shù)分布和實驗的比較。在該例中Stanton數(shù)定義為

(19)

從圖7中可見，4種來流條件下，- -模型獲得的熱流分布和實驗吻合得相當好，而基本- 模型獲得的轉捩區(qū)過長。特別對于高焓條件的Case 4，基本- 模型未能給出明顯的轉捩，而- -模型正確地捕獲了轉捩。總體而言，- -模型有效改善了高超聲速平板轉捩流動的模擬。

2.3 半錐角5°尖錐Ma=3.5超聲速轉捩流動

半錐角為5°的尖錐計算條件取自文獻[26-27]等在低擾動風洞的試驗，具體條件見表3。從表中可見，通過調(diào)整來流壓力，可以研究不同雷諾數(shù)下的尖錐轉捩特性。

圖7 高超聲速平板轉捩Stanton數(shù)Fig.7 Stanton number for hypersonic plate transition

表3 半錐角5°尖錐Ma=3.5來流條件

試驗數(shù)據(jù)是從低擾動風洞中獲得的，故計算時來流湍流度取為0.1%，來流/=10，壁面邊界條件為絕熱壁面邊界。

由于迎角為0°，計算域為1/4錐。計算網(wǎng)格是六面體網(wǎng)格，共約400萬網(wǎng)格單元，流向221個網(wǎng)格點，其中181個點在錐體表面，法向301個網(wǎng)格點，如圖8所示。

圖8 半錐角5°尖錐計算網(wǎng)格Fig.8 Computational grid for 5° half angle sharp cone

實驗采用的是絕熱壁面，轉捩的影響通過恢復因子分布體現(xiàn)。恢復因子定義為

(20)

式中:為參考溫度，是過激波之后的流場溫度;是過激波之后的總溫;是實驗或計算得到的壁面溫度。

圖9給出了計算得到的恢復因子和實驗的比較。從圖中可見，基本- 模型和- -模型都計算出了轉捩，基本- 模型得到的轉捩位置更靠前，而- -模型的轉捩位置和轉捩區(qū)發(fā)展與實驗吻合得更好，不同來流條件下的轉捩雷諾數(shù)都約為7.5×10。此外基本模型在轉捩前恢復因子有一個非物理的下降過程，文獻[9，28]中也報道過類似的情況；而- -模型的恢復因子在轉捩前一直保持層流結果，沒有非物理的下降。

圖9 恢復因子沿尖錐表面的分布Fig.9 Distribution of recovery factors along cone surface

實驗結果在湍流區(qū)的恢復因子有一個較尖銳的峰值，然后沿流向逐步緩慢下降，計算沒有觀察到這樣的現(xiàn)象，文獻[9]認為采用變湍流Prandtl數(shù)的湍流模型可能會改善計算結果。

半錐角5°尖錐=3.5的超聲速流動轉捩模擬表明，相比于基本- 模型，- -模型明顯改善了轉捩模擬結果。

2.4 半錐角5°尖錐Ma=6.0高超聲速轉捩流動

選取Horvath等在=6常規(guī)風洞中針對半錐角5°尖錐的高超聲速流動轉捩試驗，用于計算的來流條件見表4。

表4 半錐角5°尖錐Ma=6.0來流條件

計算網(wǎng)格參見圖8。計算采用等溫壁，所有計算的壁溫=0.59。在計算中入口湍流度=3.0%，入口處的(/)=16，下標inlet表示距尖錐頂點0.01 m入口處的值。

圖10給出了不同模型計算得到的無量綱傳熱系數(shù)和實驗的對比。這里=/(-)，為壁面熱流，為焓，為由Fay-Riddell駐點熱流公式給出的參考值。

從圖10中可見，本文所發(fā)展的- -模型較好地捕獲了轉捩，在不同雷諾數(shù)下轉捩起始位置和轉捩區(qū)的發(fā)展都和實驗吻合得相當好。對于Case 1和Case 2，試驗的無量綱傳熱系數(shù)存在明顯的峰值，而計算中沒有獲得類似的結果，這個現(xiàn)象和2.3節(jié)中=3.5的5°半錐角尖錐的結果類似，有可能和湍流模型對可壓縮流動的適應性相關。

圖10 傳熱系數(shù)沿尖錐表面的分布Fig.10 Distribution of heat transfer coefficients along cone surface

圖10中的Case 1和Case 4也給出了基本- 模型、層流和全湍流的結果。從圖中可見，基本- 模型預測的轉捩位置更靠前，但轉捩區(qū)的強度遠小于- -模型，轉捩區(qū)長度更長。這是因為基本- 模型的轉捩準則缺乏可壓縮性修正，使得轉捩位置提前；同時基本- 模型用于間歇因子γ方程的F過小，使得方程的生成源項太小，產(chǎn)生了不正確的轉捩發(fā)展。總體而言，相比于基本- 模型，- -模型有效改善了半錐角5°尖錐=6.0的高超聲速流動轉捩模擬精度。

3 結 論

通過一系列不同流動條件的轉捩算例對本文所發(fā)展的- -轉捩模型進行了考核檢驗，計算結果表明：

1) 在低速流動條件下，- -模型自動恢復為基本- 模型，兩者的結果幾乎完全一致。

2) 針對轉捩準則的可壓縮性修正改善了超聲速和高超聲速流動的轉捩起始位置預測。

3) 利用雷諾數(shù)比擬關系對關聯(lián)函數(shù)的修正改善了可壓縮流動轉捩區(qū)發(fā)展的預測。

4)- -模型將基本的- 模型擴展到可壓縮流動，實現(xiàn)了從低速至高速的無縫統(tǒng)一模擬能力。