萬兵

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

0 引 言

對于超聲速飛行器的設計，層流-湍流轉捩是必須要考慮的問題，也是流體力學研究的熱點。在超聲速飛行條件下，飛行器鈍頭體附近產生弓形激波。當氣流穿過前緣附近激波角變化較大的弓形激波時，激波后的熵值變化較大，形成一層橫向熵梯度較大的區(qū)域，稱之為熵層。氣體向下游流動，將包圍著飛行器表面，從而增加表面氣動熱流[1]。因此研究熵層的流動穩(wěn)定性，對于飛行器氣動熱防護的設計具有實際意義。

早期國外對鈍體模型的實驗研究發(fā)現(xiàn)，熵層對邊界層穩(wěn)定性的影響具有重要作用。Stetson[2]通過實驗研究發(fā)現(xiàn)，在邊界層外存在增長著的擾動，這表明熵層存在不穩(wěn)定性。Lysenko[3]通過鈍板實驗也測量出熵層不穩(wěn)定模態(tài)，同時指出這個不穩(wěn)定模態(tài)影響著邊界層內的第一模態(tài)。為了證實Stetson實驗等關于熵層不穩(wěn)定性的存在，Dietz & Hein[4]利用數(shù)值方法計算鈍板模型，指出熵層有較小的增長率，其擾動演化與文中[4]Aerodynamisches Institute的實驗結果吻合。Fedorov[5]根據(jù)Yakura[6]漸近法分析得到了熵層基本流方程表達式，在此基礎上進行線性穩(wěn)定性分析，得出同樣的結果。

早在實驗[3, 7-8]研究鈍度對邊界層轉捩的影響時發(fā)現(xiàn)，在小鈍度范圍內增加鈍度，轉捩位置往下游移動；當鈍度增加到一定大小時，轉捩位置反而提前，這種現(xiàn)象稱之為“transition reversal”。然而，Rosenboom & Hein[9]和Lei & Zhong[10]在計算N值時并沒有發(fā)現(xiàn)這種現(xiàn)象，且鈍度不斷增加后，N值越小。Balacumar[11]和Kara et al[12]研究不同鈍度的鈍板、鈍楔和鈍錐的慢聲波感受性時，得到的感受性系數(shù)亦是隨著鈍度的增加而始終減小。Fedorov[13]認為這種“transition reversal”的現(xiàn)象可能與熵層不穩(wěn)定性有關。究竟是否存在“transition reversal”的現(xiàn)象是值得研究的。

本文將利用直接數(shù)值模擬計算馬赫數(shù)為4.5、前緣為圓頭的0°迎角鈍板，并在此基礎上進行線性穩(wěn)定性分析，研究超聲速繞鈍板的熵層不穩(wěn)定性。

1 控制方程和計算方法

1.1 控制方程

0°迎角超聲速繞鈍板流動是二維流動問題，直角坐標系(x,y)下二維可壓縮流動N-S的無量綱守恒型方程為：

(1)

其中，U=[ρ,ρu,ρv,ρes]T,

(2)

其中，

圖1 坐標系和計算模型示意圖Fig.1 Schematic diagram of coordinatesystem and computation model

1.2 計算參數(shù)和邊界條件

本文選擇馬赫數(shù)為4.5、前緣為圓頭的0°迎角鈍板，可參考圖1的示意圖。氣體為理想氣體，黏度和熱傳導系數(shù)均采用Sutherland公式，流動條件見表1。

表1 超聲速繞鈍板流動條件Table 1 Conditions of supersonic flow over a blunt plate

對于直接數(shù)值模擬，計算模型可分為基本流計算和擾動演化。不同的計算模型邊界條件也將不同。在對稱面處，兩個模型均采用對稱邊界條件；在出口處，各變量由上游三個點進行二次外推；在壁面處，基本流計算時采用無滑移絕熱壁，擾動演化時速度脈動為0，溫度脈動的法向導數(shù)為0；在遠場，基本流計算時上邊界等于來流值，擾動演化時上邊界各變量的脈動值均為0。

1.3 數(shù)值方法

對于直接數(shù)值模擬，基本流計算模型切換到擾動演化計算模型時，為了避免切換過程中引入數(shù)值擾動，需要保證通量分裂方法和時間離散格式相同，以及邊界層和熵層區(qū)域的空間離散格式相同。在這兩個計算模型下，通量分裂方法均選擇AUSM plus分裂方法，時間離散格式均選擇三步Runge-Kutta格式。但對于空間離散，控制方程的對流項所選擇的計算格式將不同。計算基本流時，對流項在激波附近采用五階WENO格式(Jiang & Shu[14])，其他區(qū)域包括邊界層和熵層區(qū)域采用五階迎風格式；而進行擾動演化時，對流項統(tǒng)一采用五階迎風格式。

本文的計算無量綱長度為1500，即7.5 m。在進行基本流計算時，計算域的流向網格和法向網格均為不均勻網格。對于法向網格分布，壁面附近網格較密，沿法向向上逐漸變稀，并且保證邊界層區(qū)域內約包含101個法向點，熵層中廣義拐點(廣義拐點下文將介紹)附近的區(qū)域內約包含31個法向點。對于流向網格分布，頭部區(qū)域較密，沿下游逐漸變稀。這樣計算域的網格點數(shù)為1001×301。在進行擾動演化時，法向網格分布與上述一致，但流向網格為均勻分布，并保證一個擾動波波長大小的區(qū)域內至少包含21個流向點。這樣計算域的網格點數(shù)為3001×301。數(shù)值計算程序由本課題組提供，其正確性可參見文獻[15]。

2 計算結果

2.1 熵層和邊界層

圖2給出直接數(shù)值模擬得到的溫度云圖，同時給出激波、熵層外緣和邊界層外緣的位置曲線，圖2(a)給出的是流向區(qū)域為x=0～100的云圖，圖2(b)給出的是流向區(qū)域為x=100～1000的云圖，其中區(qū)域為x=400～430的云圖是紅色矩形的放大圖，粉色曲線表示激波位置，黃色曲線表示熵層外緣位置，黑色實線表示邊界層外緣位置。邊界層外緣與壁面之間為邊界層區(qū)域，邊界層外緣與熵層外緣之間為熵層區(qū)域。

圖2 溫度云圖Fig.2 Contour of temperature

邊界層外緣位置可通過總焓h0參數(shù)確定。根據(jù)Kimmel & Klein[16]的經驗，總焓由壁面開始過峰值后等于來流總焓的1.005倍，則為該流向位置的邊界層外緣位置。例如，在x=200處，如圖3(a)所示，給出總焓h0沿垂直壁面法向y的變化，黑點即為邊界層外緣位置。對于熵層外緣的位置，可根據(jù)熵變的大小來確定，Laible[17]選擇的經驗方法是基于當?shù)丶げê箪刈儲的大小來判斷，這樣在確定每個位置的熵層外緣時需計算該位置激波后的熵變，增加了工作量。本文采用新的簡單的經驗方法，是基于某一固定值來確定。即，熵層外緣滿足沿法向的熵變Δs等于駐點處熵變的0.01倍。某一位置的無量綱熵變Δs為：

圖3(b)中給出熵變Δs沿垂直壁面法向y的變化，并根據(jù)上述方法確定出熵層外緣的位置。如此鈍板流場中的熵層區(qū)域和邊界層區(qū)域也就可以確定了，如圖2所示。需要說明的是，在圖2的x=0～1000流向范圍內，沒有明顯看出熵層外緣與邊界層外緣相交的趨勢，換句話說，熵層外緣與邊界層外緣夾角很小，熵層區(qū)域的特征仍然顯著，熵層還沒被“熵吞”。在這段區(qū)域內，熵層外緣與邊界層外緣近乎平行，熵層流場在下游滿足平行流假設。

圖3 總焓和熵變以及溫度和速度沿垂直壁面法向的變化Fig.3 Total enthalpy, entropy, temperature and stream-wisevelocity along the wall-normal direction

2.2 熵層不穩(wěn)定性

圖4給出x=200的流向速度和ρ?u/?y曲線，可以看出ρ?u/?y曲線出現(xiàn)兩個極大值，故存在兩個廣義拐點，分別位于邊界層和熵層，其中熵層區(qū)域的廣義拐點記為GIP-Entropy。根據(jù)無粘不穩(wěn)定理論[18]，除了邊界層中的不穩(wěn)定模態(tài)，在基本流場中的熵層區(qū)域也有無粘不穩(wěn)定模態(tài)。

圖4 流向速度和ρ?u/?y沿法向y的變化Fig.4 Stream-wise velocity and ρ?u/?yalong the wall-normal direction

為驗證理論，針對熵層的廣義拐點，進行線性穩(wěn)定性理論(LST)分析。線性穩(wěn)定性理論能夠描述擾動的參數(shù)特征。它基于小擾動、平行流假設，忽略擾動方程的非線性項和基本量沿流向的變化，小擾動q′形式如下：

(3)

圖5給出二維擾動的空間增長率-αi及相速度隨圓頻率的變化，可以看出中性解的上支界點的無量綱圓頻率約為0.13，對應無量綱相速度約為0.92，與圖4中廣義拐點GIP-Entropy對應位置的流向速度相等，這表明該模態(tài)為熵層引起的不穩(wěn)定模態(tài)。圖6給出該不穩(wěn)定模態(tài)的特征函數(shù)幅值分布，可以發(fā)現(xiàn)特征函數(shù)幅值的峰值對應的無量綱法向位置約為5.16，與圖4中熵層中廣義拐點的位置相等，進一步證明了熵層中存在不穩(wěn)定模態(tài)。

圖5 空間增長率及相速度隨圓頻率的變化Fig.5 Spatial rate and phase velocity with circular frequency

圖6 流向速度特征函數(shù)幅值分布Fig.6 Amplitude distribution of stream-wisevelocity eigenfunction

圖7給出x=200處不同圓頻率下無量綱流向速度和溫度的特征函數(shù)幅值分布，七個圓頻率表示在圖5中。圖7 (a) 可以看出，較小圓頻率時流向速度特征函數(shù)峰值在熵層廣義拐點上，隨著圓頻率的增加，該峰值幅值逐漸減小，同時更靠近壁面的位置的幅值開始增加，并超過原峰值。可以發(fā)現(xiàn)，此時的峰值位置是圖4中Boundary-layer edge和GIP-Entropy之間的極小值對應的位置。而圖7(b)可以看出，不同頻率下溫度特征函數(shù)幅值的峰值位置均在熵層廣義拐點GIP-Entropy上。

(a) 流向速度

(b) 溫度

圖8和圖9給出x=200的中性曲線及βr=0、 0.04、 0.08、 0.12時空間增長率隨圓頻率的變化，可以看出熵層存在著三維不穩(wěn)定波，但二維不穩(wěn)定波的空間增長率最大，因此在研究熵層擾動對邊界層的影響時，選擇熵層二維不穩(wěn)定波。

圖8 中性曲線Fig.8 Neutral curve

圖9 不同展向波數(shù)下空間增長率隨圓頻率的變化Fig.9 Spatial rate along circular frequency withdifferent span-wise wave numbers

2.3 熵層擾動的演化

圖10給出熵層二維中性曲線沿流向的變化，可以看出，中性曲線是半開放半封閉的。相比于邊界層內中性曲線的左封閉右開放的形狀，熵層的形狀則相反，即左開放右封閉。換句話說，熵層擾動在上游是不穩(wěn)定的，而到下游則是穩(wěn)定的。圖11給出無量綱圓頻率為0.07的空間增長率沿流向的變化，可以看出空間增長率沿流向單調遞減，且約在x=1000之后空間增長率為負，擾動開始衰減。另外注意到增長率為萬分之一的數(shù)量級，擾動幅值變化很小。

圖10 熵層二維中性曲線沿流向的變化Fig.10 Two-dimensional neutral curve of entropy-layermode along stream-wise direction

圖11 空間增長率沿流向的變化Fig.11 Spatial growth rate along stream-wise direction

為研究熵層擾動沿流向的演化過程，在x=400處加入通過穩(wěn)定性分析得到的熵層擾動，進行擾動的數(shù)值模擬。圖12給出密度脈動的云圖，可以看出擾動主要分布于熵層區(qū)域，邊界層內擾動未見到明顯變化。

圖13給出數(shù)值模擬的速度脈動的幅值曲線沿流向的變化，同時還給出線性PSE和LST的結果?？梢园l(fā)現(xiàn)約在x=1000之前，擾動幅值緩慢上升，在臨界位置之后緩慢下降。而且線性PSE(線性PSE的程序驗證參見文獻[19])和LST的結果也與數(shù)值結果吻合得很好，由于LST對方程組作了平行流假設，說明了熵層內流動在下游具有很好的平行性。

圖12 密度脈動的云圖Fig.12 Contour of density fluctuation

圖13 速度脈動的幅值曲線沿流向的變化Fig.13 Amplitude of velocity fluctuation alongstream-wise direction

3 結 論

本文以前緣為圓頭的0°迎角鈍板為模型，通過直接數(shù)值模擬(DNS)計算基本流場，在此基礎上利用線性穩(wěn)定性(LST)分析研究了熵層中的不穩(wěn)定性，同時根據(jù)經驗方法確定了熵層區(qū)域和邊界層區(qū)域，得到結論如下：

1) 熵層中存在廣義拐點，通過穩(wěn)定性分析表明熵層中存在不穩(wěn)定擾動。另外，熵層存在三維不穩(wěn)定波，但二維不穩(wěn)定波增長率最大。

2) 擾動圓頻率較小時流向速度特征函數(shù)峰值在熵層廣義拐點上，隨著圓頻率的增加，該位置的峰值大小逐漸減小，同時邊界層外緣與熵層廣義拐點之間的某一位置的幅值開始增加，并超過原峰值大小。

3) 熵層不穩(wěn)定模態(tài)集中在前緣附近，且空間增長率很小。利用直接數(shù)值方法模擬熵層擾動的演化過程，擾動在前緣附近幅值緩慢上升，之后緩慢下降。在距離前緣較遠的范圍內DNS結果與LST結果吻合得很好，表明熵層流場在下游具有很好的平行性。

需要指出的是，究竟是否存在“transition reversal”的現(xiàn)象還需要進一步研究。

致謝:本文所使用的DNS程序由趙磊博士提供，PSE程序由張存波博士提供，在這里特別感謝。

參考文獻：

[1]Zhang Hanxin.The problem of entropy layer in the hypersonic flow[J].Acta Aeronautica et Astronautica Sinica, 1965, (2): 20-43.(in Chinese)張涵信.高超聲速運動中的熵層問題[J].航空學報, 1965, (2): 20-43.

[2]Stetson K F, Thompson E R, Donaldson J C, et al.Laminar boundary layer stability experiments on a cone at Mach 8, Part 2: blunt cone[C]//22nd Aerospace Sciences Meeting, Reno, Nevada, USA: AIAA, 1984, 84-0006.

[3]Lysenko V I.Influence of the entropy layer on the stability of a supersonic shock layer and transition of the laminar boundary layer to turbulence[J].Journal of Applied Mechanics and Technical Physics, 1990, 31(6): 868-873.

[4]Dietz G, Hein S.Entropy-layer instabilities over a blunted flat plate in supersonic flow[J].Physics of Fluids, 1999, 11(1999): 7-9.

[5]Fedorov A, Tumin A.Evolution of disturbances in entropy layer on blunted plate in supersonic flow[J].AIAA Journal, 2004, 42(1): 89-94.

[6]Yakura J K.Theory of entropy layers and nose bluntness in hypersonic flow[J].Hypersonic Flow Research, 1962, 7(1): 421-470.

[7]Stetson K F, Rushton G H.Shock tunnel investigation of boundary-layer transition at M Equals 5.5[J].AIAA Journal, 1967, 5(5): 899-906.

[8]Softley E J, Graber B C, Zempel R E.Experimental observation of transition of the hypersonic boundary layer[J].AIAA Journal, 1969, 7(2): 257-263.

[9]Rosenboom I, Hein S, Dallmann U.Influence of nose bluntness on boundary-layer instabilities in hypersonic cone flows[C]//Proc.30th Fluid Dynamics Conference and Exhibit.Norfolk, VA, USA: AIAA, 1999, 99-3591.

[10]Lei J, Zhong X.Linear stability analysis of nose bluntness effects on hypersonic boundary layer transition[J].Journal of Spacecraft and Rockets, 2012, 49(1): 24-37.

[11]Balakumar P.Stability of supersonic boundary layers over blunt wedges[C]//Proc.36th Fluid Dynamics Conference and Exhibit.San Francisco, California, USA: AIAA, 2006, 06-3053.

[12]Kara K, Balakumar P, Kandil O.Effects of nose bluntness on stability of hypersonic boundary layers over a blunt cone[C]//Proc.37th Fluid Dynamics Conference and Exhibit.Miami, FL, USA: AIAA, 2007, 07-4492.

[13]Fedorov A V.Instability of the entropy layer on a blunt plate in supersonic gas flow[J].Journal of Applied Mechanics and Technical Physics, 1990, 31(5): 722-728.

[14]Jiang G S,Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics, 1996, 126(1): 202-228.

[15]Zhao Lei.Study on instability of stationarycrossflow vortices in hypersonic swept blunt plate boundary layers[D].Tianjin: Tianjin University, 2016.(in Chinese)趙磊.高超聲速后掠鈍板邊界層橫流定常渦失穩(wěn)的研究[D].天津: 天津大學, 2016.

[16]Kimmel R L, Klein M A, Schwoerke S N.Three-dimensional hypersonic laminar boundary-layer computations for transition

experiment design[J].Journal of Spacecraft & Rockets, 1993, 34(4): 409-415.

[17]Laible A C.Numerical investigation of boundary-layer transition for cones at Mach 3.5 and 6.0[D].Tucson, Arizona, USA: University of Arizona, 2011.

[18]Lees L, Lin C C.Investigation of the stability of the laminar boundary layer in a compressible fluid[M].Kitty Hawk, North Carolina, USA: NACA, 1946.

[19]Zhang Cunbo.Research on nonlinear mode interactions relating to supersonic boundary layer transition[D].Tianjin: Tianjin University, 2017.(in Chinese)張存波.與超聲速邊界層轉捩有關的擾動非線性作用研究[D].天津: 天津大學, 2017.