熱化學非平衡高超聲速平板邊界層線性穩定性分析

陳賢亮, 符 松

(清華大學 航天航空學院, 北京 100084)

0 引 言

層流向湍流的轉捩會使壁面熱流和摩阻提高最多一個量級，因此對轉捩位置的準確預測對于高速飛行器設計具有重要意義。但是轉捩現象是極其復雜的，一個重要表征是轉捩過程強烈地依賴于邊界條件和外部環境。Morkovin等[1]根據來流環境擾動量級的不同，給出了如圖1所示的五條不同轉捩路徑，其中環境擾動最小的路徑被稱為自然轉捩，適用于背景湍流度較低的高空大氣環境。對于自然轉捩，外部擾動經過感受性機制在邊界層內形成擾動波，而后擾動首先經歷線性模態增長，再由于非線性引發模態間相互作用，最終經過更復雜的非線性過程發展為湍流。由于高超聲速邊界層流動轉捩問題的重要性和復雜性，相關研究是近期空氣動力學的熱點之一，文獻[2]和文獻[3]對此進行了綜述。

圖1 不同來流擾動幅值下轉捩路徑示意圖[1]Fig.1 Paths to boundary layer transition with different freestream disturbance amplitudes [1]

高超聲速流動的另一個重要特點是隨著馬赫數升高而快速上升的流場溫度。量熱完全氣體下高馬赫時正激波后溫度近似與馬赫數平方成正比，極高的溫度會激發空氣分子的振動能和引起空氣的化學分解甚至電離，導致量熱完全氣體假設的失效，這就是所謂的高溫(真實)氣體效應[4]。圖2給出了標準大氣壓下空氣發生不同熱化學過程的溫度范圍。基于統計熱力學，研究者們建立了對熱化學平衡態的描述方法，并通過系列擬合式實現了對熱化學平衡流場的求解[5-6]。

然而研究者們很快發現熱化學過程的時間尺度在高超聲速高焓流動下常與流動時間尺度相當，例如圖3給出了不同溫度和密度下氧氣熱化學過程的特征時間，而高超聲速流動的速度恰在km/s量級， 因此流場內的熱化學過程并不能處處達到平衡態，需要額外考慮振動能、電子能與組分質量的守恒方程，來建立對熱化學非平衡流場的模擬[7]。熱化學過程對層流場的主要影響可定性描述為：通過能量松弛和化學反應降低流場溫度，使得流場密度升高、激波層變薄；在相同馬赫數下，來流溫度和密度越高，則熱化學過程的時間尺度越小，流場越趨向于平衡態。但是熱化學過程對穩定性和轉捩的影響難以系統回答，這一方面是由于熱化學過程所涉及的尺度和參數眾多，與擾動場有復雜的相互作用，另一方面是由于較惡劣的氣動環境導致相關實驗較為缺乏，因而高溫流場模型的可靠性還有待進一步驗證。

1991年，Malik等[8]率先分析了“真實氣體效應”對平板邊界層穩定性的影響，他們比較了熱化學平衡氣體與量熱完全氣體的第二模態增長率，發現真實氣體效應會使第二模態更不穩定，這與他們預期的“真實氣體效應對不穩定模態的影響類似冷壁效果”一致，但是最不穩定波的頻率卻與預期相反地下降了。Stuckert等[9]則發現平衡流假設下“熱”壁面也會出現量熱完全氣體冷壁時才有的“超聲速模態”。天津大學的樊宇[10]、萬兵兵[11]等也對平衡流假設下馬赫數、高度、輸運模型等對穩定性和轉捩的影響進行了研究。除此之外，研究者們進一步分析了熱化學非平衡流動下邊界層的穩定性特征，發現與平衡流時有顯著不同[9,12]，從而更加關注時間尺度在流動發展和擾動傳播中的作用。同時，研究者們發現熱化學非平衡效應對邊界層穩定性的影響并不單一，而是與氣體模型、黏性模型、邊界條件等多種因素密切相關[13-15]，例如熱化學過程引起的流場溫度降低和邊界層變薄有使得第二模態增長率增加的趨勢，但能量松弛引起的脹壓黏性又會降低擾動增長率(尤其是對CO2等極性氣體)。另一方面，從穩定性方程數學求解的角度，Federov等[16]提出的高速邊界層穩定性分析新框架為研究模態演化與相互作用提供了新思路。在此基礎上，Bitter等[17]詳細分析了振動非平衡效應對馬赫數5左右極冷壁平板邊界層的影響，發現振動非平衡效應有助于在極冷壁條件下形成超聲速不穩定模態。除了線性穩定性分析之外，相關直接數值模擬[18-20]、拋物化擾動方程[21-22]等研究也正在開展。同時，T5[23]、LENS[24]、HIEST[25]、 VKI-Longshot[26]等高焓風洞也進行了針對考慮高溫氣體效應的流動穩定性和轉捩的驗證與分析工作。

圖2 標準大氣壓下空氣發生振動能激發和 化學反應的溫度范圍[4]Fig.2 Temperature ranges of air for vibrational and chemical processes under normal atmosphere [4]

(a) 振動松弛特征時間

(b) 化學分解特征時間圖3 不同溫度和密度下氧氣振動松弛和 化學分解特征時間Fig.3 Characteristic time of oxygen versus temperature and density for thermal and chemical processes

綜合來看，熱化學非平衡效應對高超聲速邊界層的穩定性和轉捩過程有重要影響，但相關研究依然是有限和不充分的。本文旨在利用線性穩定性框架，分析熱化學非平衡平板邊界層離散譜模態的演化特性，并嘗試探究不同熱化學過程如何對不穩定模態的增長率和頻率產生影響。

1 物理模型和數值方法

1.1 熱化學非平衡N-S方程

本文的計算工質為空氣，不考慮電離取經典5組分模型(N2、O2、NO、N和O)。增加振動能守恒和組分質量守恒方程，得到描述熱化學非平衡流動的N-S方程如式(1)所示，其中引入了Park[27]的雙溫度模型，即使用熱力學平動/轉動溫度T和振動溫度Tv這兩個溫度來描述流場。

(1)

式(1)中U表示守恒變量，F和Fv分別表示無黏和黏性通量，S表示源項。各項具體寫為：

hs=cptr,sT+ev,s+Δh°f,s

(2)

黏性應力采用牛頓流體本構和Stokes假設，能量擴散采用Fourier導熱定律，質量擴散采用Fick定律。混合物的黏性相關系數使用Wilke規則[28]對各組分進行加權，其中組分黏性系數由Blottner[29]給出，組分熱擴散系數采用Eucken關系[30]，組分質量擴散系數假設恒定的Schmidt數Sc=0.5。

振動能的非平衡松弛過程引入Landau-Teller方程來描述：

其中松弛特征時間τs基于Millikan和White[31]的半經驗擬合。

化學反應過程采用有限速率化學反應模型。所用5組分5反應模型寫為：

R1: N2+M? 2N+M

R2: O2+M? 2O+M

R3: NO+M? N+O+M

R4: N2+O ? NO+N

R5: NO+O ? O2+N

其中R1～R3是分解反應，R4和R5是交換反應，M表示催化劑，可取5組分中任意一個。化學反應平衡常數由Gibbs自由能得到，其中自由能擬合式由文獻[32]給出，反應速率常數使用Arrhenius公式，其中擬合系數采用Park模型[33]。

振動溫度的邊界條件與溫度邊界相一致，取為絕熱或等溫。組分邊界取決于壁面的催化條件，本文使用無催化壁。

將N-S方程(1)寫為算子形式以便后續分析：

q=[ρ,u,v,w,T,Ys,Tv],s∈[2,5]

(3)

1.2 邊界層基本流計算

穩定性分析需要高精度邊界層基本流，這可以通過求解上節N-S方程得到。但N-S方程的求解時間一般較長，為此嘗試求解邊界層方程以提高求解效率。對來流0°迎角的平板幾何，在沒有流動分離的情況下，借鑒量熱完全氣體的Lees-Dorodnitsyn相似性變換來求解邊界層方程：

利用連續性方程消去法向速度后，得到變換后的φ=[U,T,Ys,Tv]滿足下列形式的方程：

上式中非平衡源項S的存在使得即使對平板也不存在相似性解。針對于此，本文采用流向推進解[34]，在η和ξ方向分別用Chebyshev配置點法和三階差分進行離散，然后利用牛頓法進行隱式迭代。此方法求解效率很高，在每個流向位置所需迭代一般不超過6次。

1.3 模態穩定性計算

代入式(3)后得到擾動方程：

(4)

上式可整理為矩陣形式：

其中F、A、B、C、D和H為均為10×10的矩陣，且只與基本流有關；N為非線性項，在線性分析時被忽略。引入當地平行流假設與行波解假設：

模態擾動求解是一個特征值問題，矩陣規模為(19Ny)2，其中Ny是法向網格點數。本文的法向導數離散采用Chebyshev配置點法。

2 結果與分析

2.1 計算結果驗證

首先驗證基本流。選取馬赫10熱化學非平衡絕熱平板算例，來流參數與Malik[8]相同：靜溫350 K單位雷諾數6.6×106/m，N2占比78%。將邊界層方程解(BLS)與N-S方程解進行比較，其中N-S解所用的組分數目、非平衡模型和計算參數均與BLS相同。比較兩者在x=0.6 m處的剖面如圖4所示。可以看出，兩者的速度和振動溫度剖面重合得很好，而溫度和組分剖面有2%左右的偏差，這主要是由于邊界層方程解忽略了流向壓力梯度。總體來看，使用上述流向推進法求解邊界層方程有很高的計算效率和精度，適用于穩定性分析計算。

再驗證穩定性計算結果。首先選取Bitter[17]的熱非平衡、化學凍結的極冷壁平板算例，來流條件為馬赫數4.5、靜溫1500 K，靜壓10 kPa，壁溫為300 K。比較流向位置Reδ=2000處的模態演化過程如圖5所示。由于冷壁效果，此時F模態成為了不穩定模態，且在ω>0.48后成為超聲速模態,引起了增長率曲線的拐折。

(a) 溫度/速度剖面

(b) 組分剖面圖4 馬赫數10熱化學非平衡絕熱平板 邊界層解與N-S方程解比較Fig.4 Comparison between BLS and N-S solvers for Mach 10 non-equilibrium adiabatic flat-plate

圖5 馬赫數4.5熱非平衡極冷壁平板的增長率和相速度驗證

Fig.5 Verification of growth rate and phase velocity for Mach 4.5 thermal non-equilibrium flat-plate with highly cooled wall

最后驗證高馬赫數下引起較充分化學反應的算例。選取馬赫數10絕熱平板，單位雷諾數6.6×106/m、流向位置0.6 m。分別與Miró[35](來流靜溫600 K)的化學非平衡結果，和Malik[8](來流靜溫350 K)的熱化學平衡結果進行比較。與Miró比較增長率曲線如圖6所示，與Malik比較增長率曲線和第二模態擾動形函數如圖7所示。可以看出，本文穩定性計算結果同樣能與文獻較好地符合。平衡氣體算例的增長率在高頻時與文獻稍有偏離，這可能是由于黏性模型的不一致引起。

圖6 馬赫數10化學非平衡絕熱平板增長率驗證Fig.6 Verification of growth rate for Mach 10 chemical non-equilibrium adiabatic flat-plate

(a) 增長率曲線

(b) 第二模態擾動形函數圖7 馬赫數10熱化學平衡絕熱平板增長率和形函數驗證Fig.7 Verification of growth rate and shape function for Mach 10 equilibrium adiabatic flat-plate

2.2 邊界層離散譜的演化

關于離散譜的演化過程，Federov等[16]已針對量熱完全氣體進行了詳細敘述，以下檢驗熱化學非平衡邊界層是否存在相似演化過程。圖8(a)給出了馬赫數10非平衡絕熱平板邊界層(計算參數與“基本流驗證”處相同)的二維擾動全局譜。該全局譜是由離散模態和由離散點近似的連續譜共同組成。以相速度cr為分類依據，圖中連續譜出現了快聲波、慢聲波和熵/渦波等四個分支，這與非平衡N-S方程無黏通量線化矩陣的特征值分布相一致。除了連續譜，圖中還存在從連續譜中脫落的單個離散譜模態，典型之一是此時位于不穩定區的傳統意義上的第二模態，以紅色箭頭在圖中標出。

(a) ω=0.07的二維擾動全局譜

圖8(b)給出了不同離散模態的增長率和相速度隨頻率ω的變化曲線。隨著頻率從0開始逐漸升高，連續譜中會有一系列模態發生“cut-in”，即從連續譜中脫落成為離散譜。從快聲波分支中脫落的離散模態按照脫落順序依次稱為F1、F2、F3、…模態等。慢聲波中也會脫落一個離散模態，稱為S模態。以下重點關注同步過程，“同步”指的是兩個模態的相速度相互接近時發生的過程，這會導致模態增長率的間斷或突增/降。具體來說，本算例中，F和S模態首先分別與快慢聲波發生同步，即相速度分別達到1±1/Ma。隨著F1模態相速度逐漸降低，其在ω=0.06附近與熵/渦波同步、導致F1模態增長率間斷(如圖中①處虛線圈部分所示)，隨后在ω=0.07附近與S模態同步，導致S模態相速度曲線的扭曲和增長率的突增(如圖中②處虛線圈部分所示)，這樣產生的不穩定S模態即是所謂第二模態。F2及更高模態也經歷了類似過程，但同步強度已小于F1。在冷壁條件下，F與S模態的同步過程會出現分叉，即F模態也可能發展為不穩定模態，正如圖5所示。由上述分析，熱化學非平衡平板邊界層存在與量熱完全氣體類似的模態演化和同步過程。

2.3 熱化學過程對不穩定模態的影響

以下比較熱化學凍結氣體(相當于量熱完全氣體)、非平衡氣體和平衡氣體的穩定性特征，仍選用馬赫數10絕熱平板算例，比較溫度剖面和增長率曲線如圖9所示(其中T和C分別表示熱和化學過程，F、N和E分別表示凍結、非平衡和平衡態)。

(a) 溫度邊界層剖面

(b) 增長率曲線圖9 馬赫數10絕熱平板在熱化學凍結、非平衡和平衡 假設下邊界層剖面和增長率曲線比較

Fig.9 Comparison of profiles and growth rates for Mach 10 adiabatic flat-plate with frozen, non-equilibrium and equilibrium flow

隨著熱化學過程由凍結到非平衡再到平衡，邊界層溫度和厚度逐漸降低，使得非平衡和平衡態氣體的第二模態增長率和對應頻率相比量熱完全氣體都有增加，這與量熱完全氣體降低壁溫時的效果類似。但從非平衡態向平衡態的過渡卻出現相反趨勢：平衡態的第二模態增長率峰值沒有繼續增加、與非平衡態相當，而且對應頻率反而小于非平衡態。這說明熱化學過程除了降低流場溫度外，還有其他影響穩定性的途徑。

2.3.1 非平衡源項對基本流的影響

給出熱化學凍結、非平衡和平衡假設下基本流的壁面處溫度T、振動溫度Tv和組分質量分數Ys隨流向的變化如圖10所示。作為非平衡流動的參照，圖中以黑線表示的量熱完全氣體的振動能未被激發、化學反應也未發生，而以紅線表示的平衡流的壁溫和壁面組分沿流向也不變化。以藍線表示的非平衡流在x=0處為凍結狀態，而后向下游逐漸趨向平衡態。

(a) 壁面溫度和振動溫度

(b) 壁面 O2和O組分圖10 馬赫數10絕熱平板在熱化學凍結、非平衡和 平衡假設下壁面物理量隨流向變化曲線

Fig.10 Streamwise distribution of wall quantities for Mach 10 non-equilibrium adiabatic flat-plate with different disturbance assumptions

引入無量綱數Damk?hler數來描述不同過程時間尺度的相對關系：

(6)

2.3.2 非平衡源項的擾動的影響

Mack[36]詳細分析了量熱完全氣體的模態穩定性，指出第二及更高模態主要源于無黏框架。由于高馬赫數時第二模態成為主導模態，因此對熱化學非平衡的無黏擾動方程進行分析，給出方程如式(7)所示。由于Mack模態在二維時增長率最大，因此展向波數β取為0。

(7)

由于氧氣相比氮氣在更低的溫度時激發振動能和發生化學分解，因此對空氣考慮相同溫度下熱化學過程最顯著的情況YO2,∞=1，由式(7)解出振動能和組分質量的擾動形函數：

(8)

(9)

作為驗證，在完整黏性框架下分別計算包括和忽略非平衡源項擾動的穩定性方程，即分別求解：

(10)

文獻[37]將式(8)中第二模態不穩定波的頻率與邊界層厚度進行了關聯，得到以下結果：

(11)

對比式(6)，可見Dadist相比Dabase沿流向衰減更慢，所以擾動引起的非平衡效應的影響域相比基本流更靠下游，或者說擾動相比基本流更趨向于熱化學凍結態，因此有些論文中假設“基本流為非平衡態、擾動為平衡態”是不符合物理過程的。

綜上所述，在大氣來流條件和空氣5組分模型下，對由第二模態主導的高超聲速二維邊界層，擾動方程中新增的非平衡源項的擾動對穩定性的影響很小，非平衡效應主要是通過改變邊界層剖面來間接影響穩定性的。

(a) 振動溫度擾動幅值

(b) 增長率曲線

圖11 包括和忽略非平衡源項的擾動對馬赫數10熱化學非平衡絕熱平板的振動溫度擾動和增長率的影響

Fig.11 Comparison of shape function and growth rate between those with and without the disturbance of source term for Mach 10 non-equilibrium adiabatic flat-plate

2.3.3 聲速的影響

由上述分析，聲速是影響Mack模態的重要參數，但不同氣體模型假設的聲速計算式是不同的。由于非平衡源項的擾動對穩定性影響較小，因此只需比較熱化學凍結態、熱平衡化學凍結態和熱化學平衡態這三種模型。聲速表示壓力小擾動在等熵條件下的傳播速度，定義式為：

(12)

對熱化學凍結態(即量熱完全氣體)，聲速只與溫度有關：

(13)

對熱平衡和化學平衡態，等熵過程同步引起了振動溫度和組分的變化，推導出熱平衡化學凍結和熱化學平衡態這兩種假設下的聲速表達式如下所示：

以下進行數值結果分析。為控制變量，在同樣的基本流(只能是熱化學平衡態基本流)上傳播上述三種不同假設的小擾動，比較三種聲速和對應的模態增長率曲線如圖12所示。如綠色箭頭所示，本算例下熱化學凍結擾動的聲速最大，熱化學平衡擾動的聲速最小。與此相對應，聲速越小則第二和第三模態的最大增長率越大，最不穩定波對應的頻率越小。

(a) 聲速比較

(b) 增長率曲線圖12 馬赫數10熱化學平衡絕熱平板在 不同擾動假設下的聲速和增長率比較

Fig.12 Comparison of the speeds of sound and growth rates for Mach 10 equilibrium adiabatic flat-plate with different disturbance assumptions

由此，圖9中平衡氣體第二模態的增長率和頻率的反轉趨勢可以被解釋：相比于非平衡態，雖然平衡態氣體的流場溫度和邊界層厚度進一步降低，引起第二模態最大增長率和對應頻率增加，但同時平衡態假設下等熵傳播的擾動還引起了振動能和組分的變化，使得聲速降低，從而引起第二模態最大增長率的升高和對應頻率降低。由于這兩個作用的增長率峰值對應的頻率不同，因此相互疊加后造成平衡氣體第二模態的最大增長率與非平衡氣體持平、對應的頻率相比非平衡氣體出現降低。

3 結 論

本文利用空氣5組分(N2、O2、NO、N和O)模型，開展了考慮熱化學非平衡效應的高超聲速平板邊界層的線性穩定性研究。針對非平衡邊界層基本流不存在相似性解的問題，本文實現和驗證了求解邊界層非相似性解的流向推進方法，相比直接求解N-S方程有很高效率。同時，針對熱化學凍結、非平衡和平衡等不同氣體模型，本文驗證了所發展的線性穩定性計算程序的可靠性。

以此為基礎，本文重點研究了熱化學過程對模態穩定性的影響，并探究了離散譜模態的同步和演化過程。計算結果表明，熱化學非平衡邊界層中存在與量熱完全氣體類似的離散譜演化和模態同步過程。對于由第二模態主導的高超聲速二維邊界層，本文的計算和分析表明：第一，擾動所引起的非平衡效應的影響域相比基本流更靠近下游，即擾動相比基本流更趨向于熱化學凍結態；第二，擾動方程中新出現的非平衡源項的擾動項對穩定性影響很小，非平衡過程主要是通過改變基本流剖面來間接影響穩定性；第三，對考慮熱化學過程的邊界層，聲速同樣是影響第二及更高模態(Mack模態)的重要參數，熱化學平衡態假設引起的聲速計算式的變化能夠解釋邊界層溫度和厚度降低時第二模態頻率反而降低的非常規趨勢。