超聲速繞鈍板熵層不穩(wěn)定性的研究

萬(wàn)兵

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

0 引 言

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

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

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

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

1 控制方程和計(jì)算方法

1.1 控制方程

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

(1)

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

(2)

其中，

圖1 坐標(biāo)系和計(jì)算模型示意圖Fig.1 Schematic diagram of coordinatesystem and computation model

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

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

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

對(duì)于直接數(shù)值模擬，計(jì)算模型可分為基本流計(jì)算和擾動(dòng)演化。不同的計(jì)算模型邊界條件也將不同。在對(duì)稱面處，兩個(gè)模型均采用對(duì)稱邊界條件；在出口處，各變量由上游三個(gè)點(diǎn)進(jìn)行二次外推；在壁面處，基本流計(jì)算時(shí)采用無(wú)滑移絕熱壁，擾動(dòng)演化時(shí)速度脈動(dòng)為0，溫度脈動(dòng)的法向?qū)?shù)為0；在遠(yuǎn)場(chǎng)，基本流計(jì)算時(shí)上邊界等于來(lái)流值，擾動(dòng)演化時(shí)上邊界各變量的脈動(dòng)值均為0。

1.3 數(shù)值方法

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

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

2 計(jì)算結(jié)果

2.1 熵層和邊界層

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

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

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

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

圖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)兩個(gè)極大值，故存在兩個(gè)廣義拐點(diǎn)，分別位于邊界層和熵層，其中熵層區(qū)域的廣義拐點(diǎn)記為GIP-Entropy。根據(jù)無(wú)粘不穩(wěn)定理論[18]，除了邊界層中的不穩(wěn)定模態(tài)，在基本流場(chǎng)中的熵層區(qū)域也有無(wú)粘不穩(wěn)定模態(tài)。

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

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

(3)

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

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

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

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

(a) 流向速度

(b) 溫度

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

圖8 中性曲線Fig.8 Neutral curve

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

2.3 熵層擾動(dòng)的演化

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

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

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

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

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

圖12 密度脈動(dòng)的云圖Fig.12 Contour of density fluctuation

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

3 結(jié) 論

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

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

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

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

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

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

參考文獻(xiàn)：

[1]Zhang Hanxin.The problem of entropy layer in the hypersonic flow[J].Acta Aeronautica et Astronautica Sinica, 1965, (2): 20-43.(in Chinese)張涵信.高超聲速運(yùn)動(dòng)中的熵層問(wèn)題[J].航空學(xué)報(bào), 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].天津: 天津大學(xué), 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)張存波.與超聲速邊界層轉(zhuǎn)捩有關(guān)的擾動(dòng)非線性作用研究[D].天津: 天津大學(xué), 2017.