高超聲速三維邊界層橫流轉捩的數值研究

韓宇峰馬紹賢蘇彩虹

(天津大學機械學院 高速空氣動力學研究室,天津 300072)

0 引 言

高超聲速三維邊界層的轉捩預測是飛行器設計中必須考慮的重要問題,特別是對長時間在大氣中飛行的高速飛行器。由于湍流狀態下壁面摩阻和熱流比層流下大得多,只有準確地預測出轉捩位置,才能準確計算阻力和表面熱流,適當設計熱防護措施。實際中飛行器飛行時常常有后掠角,這時在邊界層內與外流流線垂直的方向上由于壓力梯度會產生二次流動,稱為橫流。橫流失穩是觸發三維邊界層轉捩至湍流的主要原因。然而,目前對橫流失穩觸發轉捩發生的機理,特別是對高超聲速三維邊界層[1-4],還并不清楚。

關于橫流不穩定性的研究,早期主要針對低速(不可壓或亞聲速、個別低超聲速)后掠翼和后掠平板的邊界層。國際上最具代表性的是德宇航Bippes及其合作者、以及美國Saric及其合作者的工作,比較全面的綜述可見文獻[3-4]。國內也開展了不少關于橫流感受性[5]、橫流穩定性[6-8]以及橫流轉捩預測模型的研究[9-10]。目前,已經有了比較一致的認識,即有橫流的邊界層中可以存在兩種由橫流不穩定性導致的流動——定常橫流渦和低頻的橫流行進波。定常渦對壁面粗糙度非常敏感,很容易被激發。激發出的定常渦經歷一定程度的線性增長后,在下游會由于非線性作用幅值出現飽和。飽和的定常渦導致平均流發生修正,使得在修正的平均流上高頻擾動再次失穩,即出現二次失穩現象,促使轉捩發生[11]。橫流轉捩中定常渦和行進波的非線性飽和過程很長,在此過程中,擾動幅值變化非常緩慢,轉捩發生的位置與開始發生非線性飽和的位置有相當長的距離。因此,若仍采用目前普遍用于Mack模態轉捩的預測方法,即以首次不穩定波(橫流定常渦或行進波)的幅值為基礎作為轉捩判據是不可行的[11-14]。若以此為判據,任何在輸入擾動和計算擾動演化過程中引入的很小的誤差,也會給轉捩位置的預測帶來很大的影響。因此,橫流失穩導致的轉捩預測相比Mack模態失穩而言,更依賴于對轉捩機理的認識和理解。

本文針對高超聲速后掠鈍板的三維邊界層,采用直接數值模擬方法(DNS)和全局不穩定性(Biglobal)的分析方法,研究了橫流定常渦模態的首次失穩和二次失穩；采用DNS計算了高超聲速三維邊界層橫流失穩轉捩至湍流的過程,分析了邊界層轉捩過程的機理。

1 基本流場

1.1 計算工況

本文選取半無限長后掠鈍板為研究對象,其前緣為一半徑為3.5 cm的圓柱形,如圖1所示。x、y、z分別為流向、法向、展向。飛行馬赫數為6,來流迎角為0°,后掠角為45°。飛行高度為30 km,來流溫度為226.5 K,以頭部半徑為特征長度的雷諾數為7.91×104。

圖1 后掠鈍板模型示意圖Fig.1 Schematic diagram of the swept blunt plate

1.2 計算方法和邊界條件

直角坐標系下三維N-S方程的守恒型控制方程為:

其中:U為守恒型變量分別表示密度、流向速度、法向速度、展向速度和溫度,c v為定容比熱；E、F、G分別為流向、法向和展向的無黏通量；E v、F v、G v為相應的黏性通量。N-S方程中具體系數矩陣見文獻[20]。

黏性系數采用Sutherland公式計算,具體為μ＝下標∞表示來流量。熱傳導系數κ由給出,Pr＝0.72,c p∞為來流定壓比熱。

“中國和印度的經驗表明，強大國家和強大社會同時出現，隨著時間的流逝，而相互平衡，相互抵消，這樣才會有較好形式的自由”。③馬克思·韋伯認為，中國在秦朝就具有現代國家治理的要素，但是沒有形成強大社會相平衡，法治要素缺乏，所以難以形成現代民主化的國家構建；相反，印度則因宗教、種姓制度形成了強大的社會要素基礎，但難以形成強大的統一國家，也影響了現代國家的構建。英國基于個人主義、權利制衡基礎上的政治治理結構，是特定條件下形成的特定結果。

由于模型是二維的,展向無限長,因此?/?z＝0。但由于流動存在后掠角,因此是三維的,即速度分量包含u、v和w。

采用基于有限差分的直接數值模擬程序計算基本流場。對流項離散采用五階WENO格式,在邊界處需進行降階處理[22]。對于邊界點,采用三點二階的偏心差分格式；對于次邊界點,采用四點三階的偏心差分格式。黏性項離散采用六階中心差分格式,時間推進采用四步四階Runge-Kutta法。壁面采用無滑移和絕熱邊界條件,遠場條件給定自由來流。頭部對稱面處采用對稱邊界條件,出口采用線性外推邊界條件。本文采用的基本流程序代碼已成功地用于后掠鈍板的基本流計算中,可見文獻[21]。

由于橫流渦擾動和二次失穩擾動是三維的,因此,對于定常橫流渦的演化,以及非定常二次失穩模態的演化,均采用三維N-S方程。

擾動演化的計算中,入口給定具體的擾動波形式,出口采用嵌邊區[23],壁面采用無滑移邊界條件,外邊界采用遠場邊界條件,展向采用周期邊界條件。數值離散格式采用與基本流計算中相同的格式。

1.3 基本流場

計算域如圖2所示。流向計算域為100個頭部半徑(3.5m),法向計算域包含激波。流向和法向網格數分別為1001和301,CFL數取為0.4。

圖2 基本流計算域示意圖Fig.2 Schematic diagram of the computing domain

求解二維的N-S方程直至所有網格點上的變量達到定常狀態。將速度在邊界層外緣的勢流方向和垂直于勢流的方向進行分解,得到勢流方向速度分量u s和橫流速度分量u c的分布,如圖3所示。可以看到,流動中橫流的量級為10-2。從頭部到下游,橫流強度逐漸減弱。

2 首次失穩——橫流定常渦的演化

首先對基本流進行線性穩定性(LST)分析。在平行流假設下,小擾動可以寫成如下形式:

其中(y)為擾動的形函數,擾動在流向和展向和時間方向具有波的形式。考慮空間模式,即擾動沿空間增長,那么頻率ω為實數,α和β為復數,其虛部的相反數分別表示擾動波在流向和展向的增長率。將上式帶入到線性化的N-S方程的擾動形式中,并忽略基本流量的流向導數及法向速度,得到準平行流假設下的線性穩定性方程:

在壁面和遠場擾動為零。這樣方程結合齊次邊界條件,構成特征值問題。通過求解可以得到邊界層擾動波的色散關系和擾動的特征函數。

圖3 勢流方向速度(a)和橫流方向速度(b)Fig.3 Potential flow(a)and crossflow velocities(b)

橫流有兩種失穩模態,一種是定常的橫流渦,另一種是低頻的橫流行進波。在高空低背景擾動情況下,定常橫流渦被認為是三維邊界層中占據主導作用的失穩模態[3]。因此,在首次失穩中,我們暫時只關注定常橫流渦。本文選取β0＝3這一定常渦擾動,來計算其在流場中的演化情況。

由于基本流中引入的是三維擾動,因此首先將基本流在展向延拓成一個波長,即L z＝2π/β0,然后在計算域的入口引入擾動來研究首次失穩擾動演化,擾動形式為:

其中,A0代表入口擾動幅值,取為10－3(y)代表線性穩定性分析得到的β0擾動的特征函數。

用DNS計算這一擾動在流場中的演化,直至流場定常。圖4給出了定常橫流渦的流向速度云圖。可以看出在上游,流向速度在展向分布比較均勻,隨著擾動向下游演化,壁面附近的流體逐漸被卷起,形成渦結構。

將DNS得到的擾動場沿展向做傅里葉變換,可以得到不同展向波數的各階譜的速度擾動幅值沿流向的變化,如圖5所示。圖中(0,1)表示加入的定常基本波(β0＝3),(0,2)和(0,3)分別表示2倍展向波數和3倍展向波數的波。可以看出,在上游基本波的幅值快速增長,在x＝20處由于幅值較大開始非線性作用產生基本流修正(0,0)和二次諧波(0,2)。當基本流修正增長到幅值約為0.1左右時,(0,1)波幅值達到飽和,隨后基本流修正也開始出現飽和。再向下游演化,擾動幅值不再增大,擾動維持在接近飽和的狀態。

圖4 定常橫流渦流向速度云圖Fig.4 Contours of the streamwise velocity at different locations

圖5 擾動各階譜的流向速度幅值Fig.5 The modal amplitudes of the streamwise fluctuation velocity

由于定常渦的非線性作用會產生基本流修正,修正后基本流的失穩特征將會出現變化。因此我們將圖5中得到的基本流修正項(0,0)加入到基本流中,再次分析穩定性特征。圖6給出了新的基本流在不同流向位置的失穩區內增長率等值線圖。同時,作為對比,也給出了原始基本流的穩定性分析結果。可以看出,在x＝20的地方,新的基本流增長率分布與原始基本流一致,這是由于基本流修正的幅值還比較小,對穩定性影響不大。在下游x＝30的位置,相比原始基本流,增長率以及失穩區域明顯減小；繼續向下游,在x＝40處新的基本流幾乎不再失穩。說明由于基本流修正的作用,流場不穩定區域減小,這意味著忽略擾動展向的變化,只考慮(0,0)基本流修正,只會得到更穩定的基本流。但是轉捩的發生依賴于流場在更大范圍內的失穩,說明具有展向尺度的二次失穩擾動可能是促進轉捩的關鍵因素。

圖6 線性穩定性分析得到的增長率等值線圖Fig.6 Growth rate contours obtained by linear stability analysis(a)for the mean flow profile after primary instability and(b)the base flow profile

3 二次失穩及轉捩至湍流的過程

首先對首次失穩的流場進行二維全局穩定性分析。全局穩定性分析方法可以研究基本流在兩個方向是快變而在另一個方向是慢變的三維流場的穩定性問題。擾動可以寫成如下形式:

其中(y,z)為擾動的形函數,擾動在流向和時間方向具有波的形式。與LST類似,經過推導可以得到關于(y,z)的二維特征值問題。具體公式見文獻[21]。

對首次失穩后流場的某一個流向站位的剖面作全局穩定性分析,可以找到不穩定的模態,稱為二次失穩模態。

在x＝36處,即首次失穩波接近飽和(幅值約為0.2)位置處進行全局不穩定分析,可以找到兩個二次模態,分別對應于z模態和y模態[13]。圖7(a)(b)分別給出了z模態和y模態的特征函數對于z模態,其最大值的位置在法向y≈1.5,即二次失穩擾動發生在飽和渦結構的側部,這一類型的模態是由流場中展向速度剪切導致的。而對于y模態,其最大值的位置在法向y≈1.6,對應飽和渦結構的頂部位置,這一類型的模態是由流場中法向速度剪切導致的。圖8給出了x＝36處,時間模式的全局穩定性分析無量綱頻率為2的y模態(b)的特征函數

圖7 x＝36,無量綱頻率為5的z模態(a)和

圖8 x＝36,二次失穩模態的增長率Fig.8 Growth rates of the secondary instability mode at x＝36

以x＝36為計算域入口,引入全局穩定性分析得到的兩個不同頻率(ω＝5和ω＝4)的二次失穩模態。為了使轉捩盡快發生,它們的初始幅值都選為A0＝0.01。使用DNS計算二次失穩擾動演化,其中流向計算域為x＝36～96,展向計算域仍為一個基本波的波長。流向、法向和展向計算域內的網格點數分別為2401、151和61。

圖9給出了計算得到的流場瞬時速度等值面。可以看到在上游靠近入口的位置,二次失穩擾動幅值很小,流場的渦結構主要呈現首次失穩后接近飽和的定常渦結構。隨著擾動向下游演化,高頻的二次失穩擾動增長起來,飽和渦結構開始扭曲,繼續向下游演化,規則的渦結構破碎,轉捩完成,最終形成湍流。

圖9 瞬時流向速度等值面(0.55)圖Fig.9 Contours of the streamwise velocity at the value of 0.55

圖10給出了不同流向位置的二次失穩擾動速度云圖以及首次失穩橫流渦速度等值線圖。從圖中可以看出在x＜51的范圍內,二次失穩擾動的形狀與首次失穩橫流渦的速度等值線形狀相近,說明首次失穩流場的剪切導致流動出現二次失穩。而在x＞51的下游,明顯觀察到二次失穩擾動形狀逐漸變得不規則,說明流場中不同頻率(展向波數)的高次諧波增長起來,流場開始出現轉捩。

為了進一步分析流場的轉捩過程,首先分析二次失穩擾動各階譜的演化特征。圖11給出了傅里葉分析得到的不同頻率的擾動幅值沿流向的變化,圖中ω＝4和ω＝5的波表示入口加入的兩個基本波,無量綱頻率ω＝0、1、2、3的波表示由基本波相互非線性作用激發的擾動,更高頻率的擾動在圖中用灰線表示。可以看到,在x＜44的范圍內,入口加入的基本波快速增長,同時由于它們的幅值較大,非線性作用得到的ω＝0和ω＝1的擾動被激發并且快速增長。繼續向下游演化,其它頻率的擾動也被激發,并且它們都會經歷一個快速增長階段,最終在x＝55附近各階譜的幅值趨于穩定,其中頻率為0的譜幅值最大,并且頻率越大,幅值越小。

圖10 不同流向位置的二次失穩擾動速度云圖以及首次失穩橫流渦速度等值線圖Fig.10 Contours of the disturbance velocities for the secondary instability mode and the contour lines of the streamwise velocities for the cross flow vortex

圖11 傅里葉變換得到的各階譜的幅值Fig.11 Fourier modal of amplitude evolution

圖12給出了壁面摩擦系數(C f)沿流向的變化,其中紫線是層流解的結果,紅線是采用SST模型得到的湍流剖面的結果,而黃線是首次失穩后橫流渦的壁面摩擦系數,黑色實線是二次失穩后的壁面摩擦系數。可以看到首次失穩后壁面摩擦系數有一定的抬升,但是抬升量不大,并沒有達到湍流的壁面摩擦系數值。說明首次失穩幅值增大到一定程度對平均流產生修正,導致壁面摩阻上升,但是由于首次失穩波已經飽和,其本身無法再使得壁面剪切增大,使得C f進一步抬升,因此轉捩并未發生。而加入二次失穩擾動之后,壁面摩擦系數在40＜x＜60的范圍內快速抬升,并達到湍流剖面的摩擦系數,說明在x＝40下游,流場開始轉捩,并逐步演化為湍流。這是由于引入高頻的二次失穩波后,在非線性作用下,各階擾動經歷快速的增長(見圖11),不斷從平均流中吸取能量,從而進一步修正平均流使得C f更加快速抬升。這樣的結果,導致壁面摩擦系數超過了相應的湍流的估計值,因而產生了“過沖”[24]的現象。當各階擾動的幅值達到一定程度,不再繼續增長,“過沖”現象逐漸減弱或消失,C f曲線會逐漸趨向于充分發展湍流的估計值。同時在C f曲線快速抬升的階段,低頻擾動快速增長并且幅值最大,說明低頻擾動的增長主導平均流修正及C f曲線的抬升。

圖12 壁面摩擦系數沿流向的變化Fig.12 Variation of the wall friction coefficient

將展向計算域拓展為兩個基本波波長(2L z),其他條件保持不變,計算了二次失穩擾動的演化,圖12中黑色虛線給出了壁面摩擦系數(C f)沿流向的變化,與展向計算域為一個基本波波長(黑色實線)結果相比,二者幾乎一致,特別是在C f曲線抬升階段,二者完全重合,說明展向計算域對橫流渦轉捩過程影響不大。

圖13 平均速度剖面的分布Fig.13 Distribution of the mean velocity

圖13給出了經過Van Driest變換后的平均速度剖面的分布,圖中虛線分別表示壁面律(u＋＝y＋)和對數律(u＋＝2.5lny＋＋4)。可以看出在61＜x＜91的范圍內,平均速度剖面與對數率和壁面率符合得很好,而且隨著向下游演化,剖面越來越接近于對數律,說明剖面逐漸向完全湍流發展。

4 結 論

本文以馬赫數為6的后掠鈍板邊界層為研究對象,采用直接數值模擬方法研究了橫流定常渦從首次失穩、二次失穩到轉捩發生的過程。得到以下結論:

1)橫流定常渦的非線性作用引起平均流修正,可使壁面摩擦系數曲線有一定程度的抬升,但是橫流渦飽和后壁面摩擦系數不再增大,轉捩不會發生。

2)在橫流渦非線性飽和(幅值約為0.2)的基礎上,發生二次失穩,產生高頻的不穩定波。二次失穩波由于非線性作用產生的低頻諧波及定常渦迅速增長,促使壁面摩擦系數急劇抬升,同時首次失穩的飽和橫流渦結構破碎,最終促使轉捩發生。

致謝：此項工作得到天津大學周恒院士的關心和指導,天津大學趙磊博士提供了DNS程序及全局穩定性分析程序,并在國家超算天津中心“天河一號”完成計算,特此致謝。