壁面溫度分布的高超聲速平板邊界層穩(wěn)定性分析及轉(zhuǎn)捩預(yù)測

劉 璐, 曹 偉

(天津大學(xué) 力學(xué)系, 天津 300072)

0 引 言

高超聲速飛行器是航空航天技術(shù)領(lǐng)域研究和發(fā)展的重點(diǎn)，有許多關(guān)鍵問題亟待解決，其中之一就是如何準(zhǔn)確預(yù)測高超聲速邊界層的轉(zhuǎn)捩位置[1-2]，因其直接影響到氣動力、氣動熱計(jì)算的準(zhǔn)確性，這是高超聲速飛行器設(shè)計(jì)的基礎(chǔ)。

在航空工業(yè)領(lǐng)域，eN方法被認(rèn)為是預(yù)測轉(zhuǎn)捩位置的最有效方法，其中，N值往往結(jié)合實(shí)驗(yàn)和飛行數(shù)據(jù)給出。但高超聲速風(fēng)洞實(shí)驗(yàn)存在諸多難點(diǎn)，且早期eN方法忽略了多種因素影響，如今正在逐步完善[3-4]。Chen等[5]對Ma=3.5條件下尖錐和平板模型的轉(zhuǎn)捩現(xiàn)象進(jìn)行了研究，發(fā)現(xiàn)絕熱壁條件下測量得到的N=10的轉(zhuǎn)捩預(yù)測位置與線性理論結(jié)果吻合良好。Malik[6]也用eN方法對Ma=3.5條件下尖錐轉(zhuǎn)捩位置進(jìn)行了研究，并將N值調(diào)整為5，這與一般公認(rèn)的N=9～11相差甚遠(yuǎn)。蘇彩虹等[7]指出了傳統(tǒng)eN方法的不足并進(jìn)行了改進(jìn)，建議N值取為3.5。

在高超聲速飛行器邊界層穩(wěn)定性分析以及轉(zhuǎn)捩預(yù)測相關(guān)研究中，壁面一般采用絕熱或等溫條件。在不同壁面條件下，賈文利[8]研究了考慮真實(shí)氣體效應(yīng)對高超聲速鈍錐邊界層穩(wěn)定性及轉(zhuǎn)捩位置的影響。研究發(fā)現(xiàn)：對于絕熱壁條件，轉(zhuǎn)捩位置隨馬赫數(shù)增大而推后；對于等溫壁條件，在相同壁面溫度下，轉(zhuǎn)捩位置隨馬赫數(shù)增大而提前；在相同馬赫數(shù)下，隨著壁面溫度降低，第一模態(tài)更加穩(wěn)定，且壁面溫度條件對穩(wěn)定性產(chǎn)生影響，進(jìn)而影響轉(zhuǎn)捩位置；壁面溫度冷卻對第二模態(tài)擾動具有不穩(wěn)定效應(yīng)[9]。1984年，Mack[10]指出：與不可壓縮流動邊界層相比，超聲速邊界層存在多個(gè)不穩(wěn)定模態(tài)；當(dāng)馬赫數(shù)大于4時(shí)，第二模態(tài)最不穩(wěn)定。在來流馬赫數(shù)Ma=8條件下，Stetson等[11-12]以半錐角7°的尖錐為模型，研究了壁面溫度等對錐體穩(wěn)定性的影響，得到了與Mack一致的結(jié)論；將通過水冷系統(tǒng)冷卻的尖錐與非冷卻尖錐的數(shù)據(jù)進(jìn)行對比，發(fā)現(xiàn)冷壁使轉(zhuǎn)捩雷諾數(shù)減小、轉(zhuǎn)捩位置提前。蘇彩虹等[13]在零迎角小鈍度球頭錐邊界層轉(zhuǎn)捩問題的研究中指出：雖然第二模態(tài)增長率遠(yuǎn)大于第一模態(tài)，但其增長區(qū)范圍短，并不總是對轉(zhuǎn)捩起決定性作用；壁面溫度條件對轉(zhuǎn)捩位置影響較大，對于等溫壁情況，轉(zhuǎn)捩由二維的第二模態(tài)主導(dǎo)，而絕熱壁情況則由三維的第一模態(tài)主導(dǎo)。曹偉[14]通過對平板邊界層轉(zhuǎn)捩的研究發(fā)現(xiàn)，無論等溫壁或絕熱壁，在來流馬赫數(shù)為8、10、12時(shí)，轉(zhuǎn)捩均由第二模態(tài)主導(dǎo)。因此，轉(zhuǎn)捩究竟由何種模態(tài)主導(dǎo)，仍有待深入研究。

Kendall[15]曾對超聲速和高超聲速邊界層轉(zhuǎn)捩進(jìn)行了系統(tǒng)研究，發(fā)現(xiàn)距離尖錐尖端由近及遠(yuǎn)的4個(gè)熱電偶測得溫度各不相同，表明實(shí)際壁面溫度既非等溫亦非絕熱。壁面溫度條件直接影響基本流剖面、擾動演化特性、流動穩(wěn)定性和轉(zhuǎn)捩位置，只有得到符合實(shí)際流動的壁面溫度，邊界層穩(wěn)定性分析和轉(zhuǎn)捩預(yù)測才能得到符合實(shí)際的結(jié)果。因此，壁面采取何種條件也是值得研究的問題。

實(shí)際飛行中，等溫壁面條件只適用于飛行初始階段。隨后，由于邊界層內(nèi)氣動加熱，壁面溫度逐漸升高，沿流向各站位溫度不再相等。長時(shí)間飛行后，壁面溫度趨于定常；但由于壁面存在傳熱，絕熱條件并不嚴(yán)格滿足，壁面溫度沿流向仍有一個(gè)分布。如何得到不同時(shí)刻壁面溫度分布，并將其作為壁面溫度條件用以預(yù)測轉(zhuǎn)捩，是相關(guān)研究的關(guān)鍵所在。

針對薄壁超聲速流動(Ma=3～10)，文獻(xiàn)[16]提出了對流換熱系數(shù)計(jì)算公式，并根據(jù)實(shí)際飛行器材料熱物性，得到了任一時(shí)刻壁面溫度分布。本文利用此方法得出壁面溫度分布，通過直接數(shù)值模擬計(jì)算出基本流場，并分析其穩(wěn)定性；采用eN方法預(yù)測轉(zhuǎn)捩發(fā)生位置，與等溫和絕熱壁面溫度條件的穩(wěn)定性分析以及轉(zhuǎn)捩預(yù)測結(jié)果進(jìn)行對比；研究了求取某一時(shí)刻壁面等效溫度作為等溫壁面條件的壁溫、并以之代替壁面溫度分布條件進(jìn)行穩(wěn)定性分析和轉(zhuǎn)捩預(yù)測的可行性。

1 數(shù)值方法

1.1 基本流的計(jì)算

其中，x、y分別為流向及法向坐標(biāo)，U為流向速度，α、ω分別為流向波數(shù)和無量綱圓頻率(下文簡稱頻率)。

采用兩種方法計(jì)算了Ma=4.5條件下、t=6.63s時(shí)刻的基本流：第一種方法是壁面采用對流換熱公式，通過直接數(shù)值模擬從初始時(shí)刻計(jì)算到t=6.63s時(shí)刻的基本流；第二種方法是求得t=6.63s時(shí)刻壁面溫度分布后，以之作為壁面邊界條件，通過直接數(shù)值模擬得到定常的基本流。圖1為兩種方法得到的基本流對比(x= 0.6m處的法向速度、溫度剖面、中性曲線和N值曲線的對比)，可以看出兩者吻合很好。由于第二種方法更為省時(shí)，本文采用該方法得到基本流并進(jìn)行穩(wěn)定性分析。其詳細(xì)計(jì)算方法為：

(a) 速度

(b) 溫度

(c) 中性曲線

(d) N值曲線

Fig.1Comparisonofthebaseflowcomputedbytwodifferentcalculationmethods

用埃克特[17]參考焓方法求得平板上不同流向位置的對流換熱系數(shù)：

(1)

其中，hx為局部對流換熱系數(shù)，Pr=0.72，Tw、Hw為壁面溫度、壁面焓，Tr、Hr為恢復(fù)溫度、恢復(fù)焓。Tr、Hr可由經(jīng)驗(yàn)公式獲得。ρ*、μ*為參考焓H*下的密度和粘性系數(shù)。

對于超聲速流動，固壁與氣體間的局部對流換熱滿足傳熱公式：

qw=hx(Tw-Tr)

(2)

qw為熱流密度。將式(1)代入式(2)，可得熱流密度為：

(3)

平板壁面溫度增加量ΔT用集總參數(shù)法求得：

(4)

壁面溫度Tw可由來流溫度與壁面溫度增加量對時(shí)間的積分相加得出[16]：

(5)

c、ρ為金屬平板的比熱容、密度，δ為平板厚度的1/2，三者取值分別為520J/(kg·K)、4.7×103kg/m3和1mm。t表示在流場中經(jīng)過的時(shí)間，S、m分別表示所計(jì)算平板模型的面積和質(zhì)量。

圖2給出了馬赫數(shù)為4.5、6.0和7.0時(shí)、4個(gè)不同時(shí)刻(120、300、400和600s)的壁面溫度分布，以及等溫壁面溫度(取來流溫度)和絕熱壁面溫度沿流向的分布曲線?？梢钥闯觯旱葴睾徒^熱的壁面溫度沿流向不變；4個(gè)不同時(shí)刻的壁面溫度，越靠近前緣溫度越高。在t=120s時(shí)刻，Ma=4.5時(shí)，前緣溫度接近絕熱溫度，Ma=6.0和7.0時(shí)，前緣溫度達(dá)到絕熱溫度；對于其他3個(gè)時(shí)刻，前緣溫度都已達(dá)到絕熱溫度。在4個(gè)時(shí)刻，溫度都沿流向逐漸降低，且隨時(shí)間增加，逐步向絕熱溫度分布線一側(cè)靠近。

將式(5)得到的壁面溫度作為直接數(shù)值模擬的壁面溫度條件，分別計(jì)算馬赫數(shù)為4.5、6.0和7.0時(shí)、不同時(shí)刻(120、300、400和600s)壁面溫度分布的基本流場。等溫和絕熱的基本流場可直接由數(shù)值模擬(或相似性解)方法求得。

圖2 不同壁面溫度條件下的溫度分布曲線

1.2 轉(zhuǎn)捩預(yù)測方法

小擾動的演化可以用線性穩(wěn)定性理論來研究，而eN方法就是基于線性穩(wěn)定性理論的轉(zhuǎn)捩預(yù)測方法[18-19]。根據(jù)線性穩(wěn)定性理論，小擾動部分可以寫為行波形式：

φ′=φ(y)ei(αx+βz-ωt)

(6)

其中，φ=(U,V,W，T,p)，V、W分別為法向及展向速度，φ′為對應(yīng)的擾動量，z為展向坐標(biāo)，β為展向波數(shù)。對于時(shí)間模式，α和β為實(shí)數(shù)，ω為復(fù)數(shù)，ω的虛部對應(yīng)于擾動幅值沿時(shí)間的增長率；對于空間模式，β和ω為實(shí)數(shù)，α為復(fù)數(shù)，α虛部的相反數(shù)為擾動幅值沿流向的增長率。

eN方法中的N為幅值放大因子，定義為：

(7)

其中，A為小擾動幅值，A0為擾動初始幅值。對于高超聲速情況，目前還沒有系統(tǒng)的實(shí)驗(yàn)給出大小合適的N值作為轉(zhuǎn)捩判據(jù)。因此，根據(jù)文獻(xiàn)[7]和[20]，取擾動初始幅值為0.3‰，以幅值達(dá)到1.5%作為轉(zhuǎn)捩開始位置，代入式(7)求得N≈4。

對于第二模態(tài)，其展向波數(shù)β=0，下文給出的N值曲線為不同頻率擾動波增長曲線的包絡(luò)線；對于第一模態(tài)考慮的是三維波，由于展向波數(shù)β不為0，因此，對不同展向波數(shù)、不同頻率的擾動波增長率進(jìn)行積分，找出這些擾動波最早達(dá)到e4的N值曲線。

2 壁溫沿流向分布對基本流、流動穩(wěn)定性及轉(zhuǎn)捩位置預(yù)測的影響

2.1 不同壁面溫度條件對基本流剖面的影響

圖3給出了Ma=6.0、x=4m處、3種不同壁面溫度條件(等溫、絕熱、壁面溫度分布)下的法向速度剖面?？梢钥闯觯涸诓煌诿鏈囟葪l件下，邊界層中的速度不同，其中，等溫壁面溫度條件下，邊界層中的速度最大，速度邊界層最薄；絕熱壁面溫度條件下，邊界層中的速度最小，速度邊界層最厚。壁面溫度分布條件下，邊界層中的速度居中；在計(jì)算的4個(gè)時(shí)刻，邊界層內(nèi)同一法向位置的速度隨時(shí)間增加而減小，邊界層厚度增大。

圖3 不同壁面溫度條件下的速度剖面

圖4給出了Ma=6.0、x=4m處、3種不同壁面溫度條件下的法向溫度剖面。與法向速度剖面類比，可以得到相似結(jié)論：絕熱壁溫度最高，等溫壁溫度最低；溫度分布條件邊界層內(nèi)同一法向位置的壁面溫度隨時(shí)間增加而升高，并逐漸接近絕熱壁面條件的溫度剖面。

Ma=4.5和7.0時(shí)的影響規(guī)律與此相同，從略。

圖4 不同壁面溫度條件下的溫度剖面

2.2 不同壁面溫度條件流動穩(wěn)定性及轉(zhuǎn)捩預(yù)測的結(jié)果比較

圖5分別給出了馬赫數(shù)為4.5、6.0、7.0時(shí)、不同壁面溫度條件下(壁面溫度分布的壁面條件以120s時(shí)刻代表，下文同)二維的第一和第二模態(tài)的中性曲線對比。圖中，橫坐標(biāo)為中性曲線上的點(diǎn)到平板前緣的距離，即流向坐標(biāo)；縱坐標(biāo)為不穩(wěn)定波的頻率(圖6同)?？梢钥闯?，在相同流向位置，基本流的壁面溫度條件不同，中性曲線對應(yīng)的不穩(wěn)定區(qū)間頻率范圍不同。其中，絕熱壁面中性曲線對應(yīng)的不穩(wěn)定區(qū)間頻率較小，不穩(wěn)定區(qū)出現(xiàn)的臨界位置離平板前緣距離較遠(yuǎn)；等溫壁面中性曲線對應(yīng)的不穩(wěn)定區(qū)間頻率較大，范圍也比較大，不穩(wěn)定區(qū)出現(xiàn)的臨界位置離平板前緣較近；壁面溫度分布條件的中性曲線對應(yīng)的不穩(wěn)定區(qū)間頻率在等溫和絕熱壁面溫度條件之間，不穩(wěn)定區(qū)出現(xiàn)的臨界位置也在等溫和絕熱壁面溫度條件之間。

圖6分別給出了壁面溫度分布條件下、馬赫數(shù)為4.5、6.0和7.0時(shí)、4個(gè)時(shí)刻二維的第一和第二模態(tài)的中性曲線比較結(jié)果。可以看出，不同時(shí)刻壁面溫度分布條件的中性曲線相同流向位置對應(yīng)的不穩(wěn)定區(qū)間范圍也不同，隨著時(shí)間增加，不穩(wěn)定區(qū)間向低頻方向移動，越來越靠近絕熱一側(cè)。馬赫數(shù)越高，隨著時(shí)間增加，第一模態(tài)越早出現(xiàn)。時(shí)間越靠后，第一模態(tài)越靠近平板前緣，不穩(wěn)定區(qū)間也越大。

由圖5、6可以看出，在絕熱壁面條件下，中性曲線第一模態(tài)的不穩(wěn)定區(qū)間比其他條件下更大，即壁面加熱使第一模態(tài)的不穩(wěn)定區(qū)間變大；隨馬赫數(shù)增大，第二模態(tài)不穩(wěn)定區(qū)間逐漸與第一模態(tài)接近，直至重疊；在等溫壁面條件下，沒有找到第一模態(tài)的中性曲線。

以馬赫數(shù)4.5的情況為例，給出不同壁面溫度條件下、x=0.6m處、三維的第一模態(tài)和二維的第二模態(tài)增長率隨頻率的變化曲線，如圖7所示。可以看出，第一模態(tài)增長區(qū)隨時(shí)間增加而增大，在絕熱壁面溫度條件下的增長區(qū)最大。在等溫壁面溫度條件下，已很難找到第一模態(tài)增長率隨頻率的變化曲線；第二模態(tài)的增長區(qū)隨時(shí)間的增加而減小。在等溫壁面溫度條件下的增長區(qū)最大。無論何種壁面條件，第二模態(tài)的最大增長率都大于第一模態(tài)的最大增長率，并且絕熱壁面溫度條件下第二模態(tài)的最大增長率與第一模態(tài)的最大增長率之比最小，對于等溫壁則相反。可以得出如下結(jié)論：第一模態(tài)在絕熱邊界層中最不穩(wěn)定；第二模態(tài)在等溫邊界層中最不穩(wěn)定；在壁面溫度分布條件下，時(shí)間越長，第一模態(tài)越不穩(wěn)定，第二模態(tài)越穩(wěn)定。

圖5 不同壁面溫度條件下的中性曲線

圖6 壁面溫度分布條件下不同時(shí)刻的中性曲線

(b) 第二模態(tài)

Fig.7Variationofgrowthratewithfrequencyunderdifferentwalltemperatureconditions

即壁面冷卻使得第二模態(tài)更加不穩(wěn)定、第一模態(tài)更加穩(wěn)定。Ma=6.0和7.0的規(guī)律與此相同，從略。

圖8、9分別展示了Ma=4.5、6.0和7.0時(shí)、不同壁面溫度條件下N值曲線對比，以及壁面溫度分布條件下不同時(shí)刻的N值曲線對比。同時(shí)，圖中給出了三維的第一模態(tài)幅值放大倍數(shù)大于等于e4的情況。在Ma=4.5、6.0和7.0時(shí)，等溫壁面溫度條件下的轉(zhuǎn)捩位置比絕熱和壁面溫度分布條件下的更靠近平板前緣，由第二模態(tài)主導(dǎo)轉(zhuǎn)捩，且馬赫數(shù)越大，位置越靠前；而在絕熱壁面溫度條件下，馬赫數(shù)越大，轉(zhuǎn)捩位置越靠后。Ma=4.5和6.0時(shí)，轉(zhuǎn)捩由第一模態(tài)主導(dǎo)，Ma=7.0時(shí)，轉(zhuǎn)捩由第二模態(tài)主導(dǎo)；在壁面溫度分布條件下，轉(zhuǎn)捩位置都由第二模態(tài)主導(dǎo)，相同馬赫數(shù)下，隨著時(shí)間增加，轉(zhuǎn)捩位置越來越靠后；同一時(shí)刻，隨著馬赫數(shù)增大，轉(zhuǎn)捩位置提前。Ma=7.0、不同時(shí)刻的N值曲線對比則比較特殊，300s之后的N值包絡(luò)線都出現(xiàn)在300s的N值曲線左側(cè)，轉(zhuǎn)捩位置也都在300s之前，且隨著時(shí)間增加、溫度升高，轉(zhuǎn)捩位置越來越靠前。分析認(rèn)為：這種情況與400s以及大于400s時(shí)第一、二模態(tài)的中性曲線的交叉重疊有關(guān)。表1列出了不同壁面溫度條件的平板邊界層轉(zhuǎn)捩位置。

圖8 不同壁面溫度條件下的N值曲線

圖9 壁面溫度分布條件下不同時(shí)刻的N值曲線

Isothermal120s300s400s600sAdiabaticMa=4.51.42.9---2.5(1st)Ma=6.01.32.52.93.23.62.7(1st)Ma=7.01.22.43.13.02.92.9注:“1st”表示由第一模態(tài)主導(dǎo)轉(zhuǎn)捩,未加說明的皆由第二模態(tài)主導(dǎo)轉(zhuǎn)捩;“-”表示在計(jì)算域內(nèi)N值達(dá)不到4。

從表1可以看出，第二模態(tài)主導(dǎo)等溫和壁面溫度分布邊界條件的轉(zhuǎn)捩位置，即壁面溫度低的第二模態(tài)更不穩(wěn)定；第一模態(tài)在絕熱壁面條件、Ma=4.5和6.0的情況下主導(dǎo)轉(zhuǎn)捩，即壁面溫度高的第一模態(tài)更不穩(wěn)定，且轉(zhuǎn)捩位置較第二模態(tài)明顯提前。

對某一時(shí)刻，若能通過求取等效溫度作為等溫條件的壁溫，直接求得相似解以代替壁面溫度分布的基本流，而后進(jìn)行穩(wěn)定性分析及轉(zhuǎn)捩位置預(yù)測，這樣就可以簡化計(jì)算、節(jié)省計(jì)算時(shí)間。為此，以Ma=7.0、120和300s時(shí)刻的平均溫度(將壁面溫度積分，再除以流向長度，得到平均溫度)為例，分別以3.86和6.05倍的來流溫度作為等溫條件的壁面溫度來計(jì)算的基本流，與120和300s壁面溫度分布壁面條件的基本流得到的N值包絡(luò)進(jìn)行比較，結(jié)果發(fā)現(xiàn)兩者相差較大，如圖10所示。

(a) t=120s

(b) t=300s

Fig.10ComparisonofenvelopecurvesofNvaluesfortwobaseflows

選取若干不同溫度作為等溫壁面條件的壁溫，直接求得相似解，將其作為基本流場，與采用溫度分布通過直接數(shù)值模擬計(jì)算出的基本流場，分別進(jìn)行穩(wěn)定性分析及轉(zhuǎn)捩位置預(yù)測。結(jié)果發(fā)現(xiàn)，這兩種方法得到的N值曲線有較大差異，難以找到一個(gè)可以與壁面溫度分布等效的壁溫。

3 結(jié) 論

(1) 對于Ma=4.5、6.0和7.0的超聲速、高超聲速平板邊界層而言，壁面溫度越接近絕熱壁溫度，第一模態(tài)越不穩(wěn)定，第二模態(tài)越穩(wěn)定。在Ma=4.5和6.0的絕熱壁面溫度條件下，主導(dǎo)轉(zhuǎn)捩的是第一模態(tài)，其余情況轉(zhuǎn)捩均由第二模態(tài)決定。

(2) 壁面邊界條件對轉(zhuǎn)捩位置預(yù)測的影響很大。等溫壁面條件下的不穩(wěn)定頻率最高，轉(zhuǎn)捩位置較其他情況最靠近平板前緣，且馬赫數(shù)越大轉(zhuǎn)捩位置越靠前；絕熱壁面條件下的不穩(wěn)定頻率最低，馬赫數(shù)越大，轉(zhuǎn)捩位置越靠后；壁面溫度分布條件下，平板的壁面溫度越接近絕熱壁溫度，不穩(wěn)定區(qū)間越向低頻方向移動，且轉(zhuǎn)捩位置更靠后。與等溫壁面溫度條件類似的是，馬赫數(shù)越大，轉(zhuǎn)捩位置越靠前。需要指出的是，在Ma=7.0條件下，在400s以及大于400s時(shí)，第一和二模態(tài)發(fā)生重疊，故不完全滿足上述規(guī)律。

(3) 對于壁面溫度分布的邊界條件，很難找到某一個(gè)溫度作為等溫壁溫來代替。根據(jù)現(xiàn)有計(jì)算結(jié)果，仍需利用不同時(shí)刻的壁面溫度分布作為壁面邊界條件，通過直接數(shù)值模擬來獲得基本流。