劉 璐, 曹 偉

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

0 引 言

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

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

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

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

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

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

1 數(shù)值方法

1.1 基本流的計算

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

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

(a) 速度

(b) 溫度

(c) 中性曲線

(d) N值曲線

Fig.1Comparisonofthebaseflowcomputedbytwodifferentcalculationmethods

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

(1)

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

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

qw=hx(Tw-Tr)

(2)

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

(3)

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

(4)

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

(5)

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

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

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

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

1.2 轉捩預測方法

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

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

(6)

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

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

(7)

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

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

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

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

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

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

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

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

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

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

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

圖6分別給出了壁面溫度分布條件下、馬赫數(shù)為4.5、6.0和7.0時、4個時刻二維的第一和第二模態(tài)的中性曲線比較結果?？梢钥闯?，不同時刻壁面溫度分布條件的中性曲線相同流向位置對應的不穩(wěn)定區(qū)間范圍也不同，隨著時間增加，不穩(wěn)定區(qū)間向低頻方向移動，越來越靠近絕熱一側。馬赫數(shù)越高，隨著時間增加，第一模態(tài)越早出現(xiàn)。時間越靠后，第一模態(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ū)隨時間增加而增大，在絕熱壁面溫度條件下的增長區(qū)最大。在等溫壁面溫度條件下，已很難找到第一模態(tài)增長率隨頻率的變化曲線；第二模態(tài)的增長區(qū)隨時間的增加而減小。在等溫壁面溫度條件下的增長區(qū)最大。無論何種壁面條件，第二模態(tài)的最大增長率都大于第一模態(tài)的最大增長率，并且絕熱壁面溫度條件下第二模態(tài)的最大增長率與第一模態(tài)的最大增長率之比最小，對于等溫壁則相反?？梢缘贸鋈缦陆Y論：第一模態(tài)在絕熱邊界層中最不穩(wěn)定；第二模態(tài)在等溫邊界層中最不穩(wěn)定；在壁面溫度分布條件下，時間越長，第一模態(tài)越不穩(wěn)定，第二模態(tài)越穩(wěn)定。

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

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

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

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

圖9 壁面溫度分布條件下不同時刻的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)主導轉捩,未加說明的皆由第二模態(tài)主導轉捩;“-”表示在計算域內N值達不到4。

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

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

(a) t=120s

(b) t=300s

Fig.10ComparisonofenvelopecurvesofNvaluesfortwobaseflows

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

3 結 論

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

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

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