適用于可壓縮流動的 -Re t-fRe轉(zhuǎn)捩模型

劉清揚(yáng)，雷娟棉，劉周，石磊，周偉江

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

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

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

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

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

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

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

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

1.1 轉(zhuǎn)捩準(zhǔn)則的可壓縮性修正

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

max()=2193·

(1)

式中：和的定義為

=，=

(2)

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

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

=(2193 )

(3)

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

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

max()=·2193·

(4)

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

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

lg()=·(lg())+

·lg()+

(5)

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

(6)

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

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

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

=050×+022×+028×

(7)

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

(8)

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

從修正函數(shù)的構(gòu)造過程可見，該修正函數(shù)同時考慮了馬赫數(shù)和壁溫的影響。圖2給出了不同和不同/下修正函數(shù)的曲線，從圖中可見，隨著變小，在/逐漸趨近于1時，的值趨近于1，恢復(fù)為不可壓絕熱壁面平板邊界層流動的結(jié)果。而對高超聲速飛行器繞流，/通常在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ù)后，需要將其應(yīng)用至轉(zhuǎn)捩判斷中。一種直接的思路就是在式(3)的分母上乘以的值，然而的計算需要用到邊界層外緣的和，獲得這些值對于當(dāng)前基于非結(jié)構(gòu)網(wǎng)格的并行CFD代碼是十分困難的。

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

=()·()

(9)

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

=()·()·

(10)

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

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

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

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

(11)

(12)

(13)

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

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

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

(14)

源項的表達(dá)式為

(15)

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

基本- 轉(zhuǎn)捩模型的 定義為

(16)

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

(17)

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

雖然可以構(gòu)造一個更復(fù)雜的函數(shù)來修正上述關(guān)系，但實踐發(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。

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

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

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

2 γ-Reθ t-fRe轉(zhuǎn)捩模型驗證

2.1 低速平板轉(zhuǎn)捩流動

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

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

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

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

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

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

2.2 高超聲速平板轉(zhuǎn)捩流動

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

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

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

(19)

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

2.3 半錐角5°尖錐Ma=3.5超聲速轉(zhuǎn)捩流動

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

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

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

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

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

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

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

(20)

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

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

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

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

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

2.4 半錐角5°尖錐Ma=6.0高超聲速轉(zhuǎn)捩流動

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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