逯學志,黃章峰,2,羅紀生
(1.天津大學力學系,天津300072;2.空氣動力學國家重點實驗室,四川 綿陽621000)
直接數值模擬(Direct Numerical Simulation,DNS)可提供流場的高分辨率的全部信息,是研究湍流機制的主要工具之一,尤其是對超聲速和高超聲速邊界層流動。然而采用空間模式的直接數值模擬(Spatial Direct Numerical Simulation,SDNS),其入流條件應能反映上游的特性,對于湍流計算而言入流條件如何給是一大難題。本文介紹了一種新的提法,即將由時間模式直接數值模擬(Temporal Direct Numerical Simulation,TDNS)得到的充分發展湍流流場(只需要一個時刻的流場)以適當方式轉換成SDNS的入口流場。計算結果證實該方法可行且有效。
對于高超聲速飛行器,下表面往往有兩到三個壓縮折角,飛行器下部的氣流經過這些壓縮折角后進入發動機。如果折角處流動狀態為層流狀態,往往會在折角處發生流動分離,從而導致氣流偏離下表面,只有較少部分的氣體流入了發動機,嚴重影響了發動機的性能。然而在邊界層的前部,流動一般是層流,而且自然狀態下層流段可以很長。這時為使轉捩提前,需要采用加人工粗糙度的方法,使之在第一個壓縮折角前就發生強迫轉捩。
美國主要采用地面試驗和飛行試驗的驗證方法來研究各種轉捩裝置的效果[1-3],并通過傳熱系數是否快速上升來判斷人工粗糙度促使轉捩是否成功。但實際上,實驗和數值模擬結果表明,由于壁面粗糙引起的流向渦也會導致St上升[4-5],如圖1所示,其趨勢和由于轉捩導致的很相似。其中St的定義為:因此傳熱系數上升既可能是由于轉捩導致的,也可能是由于較強的流向渦引起的,不是強迫轉捩后達到湍流狀態的可靠判據。有必要對添加湍流發生器后下游流場是否真的達到了湍流狀態進行研究。


圖1 由于壁面粗糙引起的流向渦導致傳熱系數St曲線上升Fig.1 Increasing of St led by stream-wise vortex caused by the wall roughness
本文利用充分發展湍流的相似性特性和可維持特性來研究強迫轉捩發生的可能性。采用空間模式直接數值模擬(SDNS)的方法,在入口處添加包含湍流信息的擾動,如果在下游流場能夠恢復湍流狀態,說明強迫轉捩有可能發生,此時入口擾動與湍流信息存在一定的關系;如果在下游流場趨近層流狀態,說明即使是在入口添加了包含湍流信息的擾動也無法達到湍流狀態,強迫轉捩很有可能不會發生,此時需要考慮更下游的位置。
DNS的控制方程為無量綱的守恒型的N-S方程,不做任何假設。在結構化的網格上采用有限差分法進行空間離散,并采用高精度的緊致差分格式計算;時間上采用滿足TVD(Total Variation Diminished)特性的高階Runge-Kutta方法離散;為了使用迎風差分格式,對非線性項還根據Jacobian矩陣特征值的正負進行通量分裂。除了給出合理的初始條件外,還需要給定邊界條件。對于邊界層流而言,在流向和展向均為周期性邊界條件,而在法向,一個是固壁一個是出流條件。
采用時間模式直接數值模擬(TDNS)的好處是計算量小,可以一直從層流計算到充分發展湍流;其缺點是流向存在局部平行流假設。研究結果表明,對于零壓力梯度的高速平板邊界層,其雷諾數很大,空間模式直接數值模擬(SDNS)和TDNS的計算結果相比,在定性和定量上的差別都很小[6]。
與TDNS不同,空間模式直接數值模擬(SDNS)在流向還需要入流條件,而入流條件應能反映上游的流動特性。對于湍流計算而言,入流條件如何給是一大難題。
采用TDNS在某一時刻的結果,這時有確定的邊界層厚度,有相應的Reynolds數。利用充分發展湍流的相似性特性,對TDNS在某一時刻的結果在法向進行相似變換,就可以得到指定Reynolds數下的流場,這一流場的Reynolds數也就是SDNS入口處的Reynolds數。
轉換的步驟如下[7-8]:(1)采用 TDNS計算得到零壓力梯度的平板邊界層的充分發展湍流流場,取某一時刻的結果XT(x,y,z),其湍流邊界層厚度為δT;(2)采用數值方法求出具體問題的層流解Xc(x,y,z),其計算域入口的層流邊界層厚度為δC;(3)假定流體在計算域入口處之前剛剛完成轉捩,則在計算域入口已是湍流,其邊界層厚度變為kδδC,其中kδ取2到3.5;(4)將TDNS的結果在法向做相似變換XT(x,z),其中=ykδδC/δT,并分解為平均量與脈動量之和;(5)在湍流邊界層內采用TDNS的平均流場,在湍流邊界層外采用SDNS的層流解 Xc(x0,y,z),中間采用光滑函數進行連接,從而得到SDNS入口處的平均流場;(6)以適當方式將TDNS的脈動量X′T空間分布轉換成SDNS所需的入口流場的時間序列,從而得到SDNS入口處的脈動流場。
轉換的示意圖如圖2所示。將轉換關系寫成x=ct,其中x的原點在TDNS中流場的出口處,方向朝向上游。這樣做的基本概念是流場中x處的擾動以速度c向下游傳播,在t時刻將到達出口處,而這正好可以作為SDNS的入口擾動。當x的值大于TDNS的流向計算域長度后,可利用流向周期性邊界條件的性質,減去流向計算域長度即可。這里取c=0.92[7]。
該方法已經被證實是解決空間模式直接數值模擬(SDNS)入口條件的一種有效方法,并且對不同Mach數、Reynolds數、壁面條件、典型邊界層類型均可行[7-8]。

圖2 由TDNS得到的某時刻充分發展湍流流場轉換成SDNS的入口流場的示意圖Fig.2 Schematic diagram of the transformation from TDNS to SDNS
本文選取某高超聲速飛行器的對稱面作為研究對象,物理模型的頭部為圓弧,下表面有三個坡道。計算域如圖3所示。以來流速度U0=1789m/s和頭部半徑r=2.5mm為基本參量進行無量綱化,無量綱參數Ma=6,單位Re=5.0×106,攻角為0°。采用SDNS的方法計算得到二維層流解,其網格參數為:流向2300個點;法向201個點,頭部激波內部的點數不少于170個,靠近壁面的1個邊界層名義厚度內的點數不少于75個,距離壁面最近的網格尺度不超過0.002個無量綱長度。
截取二維SDNS中下表面某坡道內的層流流場,并在展向做擴展,作為三維SDNS的層流解,其中x方向沿著壁面向右,y方向垂直壁面向下,z方向垂直xy平面向外。三維SDNS中的入口位置是指入口壁面處在二維SDNS計算域的x坐標值(圖3(b)所示),流向長度指x方向的長度,展向長度指z方向的長度,為21個無量綱長度,有128個網格點。本文選取文獻[7]中的平板湍流邊界層TDNS的結果,結合SDNS入口處層流邊界層厚度,按照1.3節給出三維SDNS的入流條件。

圖3 計算域示意圖Fig.3 Schematic diagram of computational domain
在下表面某位置添加湍流發生器后,此湍流發生器為下游流體提供了有限幅值的擾動,有可能導致強迫轉捩。如果在下游某個位置發生了轉捩,那么在轉捩后的湍流流場應該能夠維持湍流狀態。

圖4 第一個坡道的計算結果Fig.4 Results in the first ramp
入口處的不同位置代表湍流發生器的不同位置,入口處加擾動的強度代表發生器不同的形狀;通過考察湍流的可維持特性來研究湍流發生器對流場的影響,如果在入口引入湍流信息后,流場很快能調整到湍流狀態,說明在計算域內通過加湍流發生器能夠達到湍流狀態,否則在計算域范圍內難以達到湍流狀態。本文考察了入口處不同位置和不同湍流度的情況,共9種情況,詳情見表1,其中最大脈動速度均方根是指三個方向脈動速度的平方和的均方根沿法向剖面的最大值。

表1 計算工況Table 1 Computational cases
圖4(a)和圖4(b)給出了第一個坡道內3個算例的壁面摩擦系數沿流向的分布、流向平均速度剖面在靠近出口處沿法向的分布。可以看出,情況A和B的壁面摩擦系數逐漸減小,并趨近于層流狀態的值,而且入口位置越靠前,越趨近層流解。從圖4(b)可以看出,三種情況的平均速度剖面不同程度地偏離了層流解,并且入口位置越靠前、入口湍流脈動幅值越小,平均速度剖面越不飽滿。這表明在第一坡道內引入較小幅值的擾動,流場會趨近于層流狀態。
情況C中的壁面摩擦系數在x=0.28m之前是快速減小,而在之后減小的幅度較小,并且其幅值接近湍流狀態的幅值。此外在情況A、B和C三種的流向平均速度剖面來看,情況C的流向平均速度剖面最飽滿,有可能是湍流狀態。
圖4(c)給出了第一個坡道內情況C中不同位置的湍流Mach數沿法向的分布,可以看出在x=0.28m之后的湍流Mach數幾乎重合。根據湍流的相似特性,說明此時流場達到了湍流狀態。這表明在第一坡道內引入較大幅值的擾動,流場有可能趨近湍流狀態,強迫轉捩有可能發生。
圖5(a)和圖5(b)給出了第二個坡道內3個算例的壁面摩擦系數沿流向的分布、流向平均速度剖面在靠近出口處沿法向的分布。可以看出,添加擾動的入口位置越靠后,壁面摩擦系數偏離層流解越多,平均速度剖面越飽滿。當擾動的脈動速度均方根值最大值不小于0.15時,下游的流場均趨近湍流狀態。
圖5(c)給出了第二個坡道內情況C中不同位置的湍流Mach數沿法向的分布,可以看出在x=0.54m之后的湍流Mach數也同樣具有相似性。這一點與圖5中情況C中x>0.54m時壁面摩擦系數較平緩減小的結論是一致的,說明此時流場達到了湍流狀態。

圖5 第二個坡道的計算結果Fig.5 Results in the second ramp

圖6 第三個坡道的計算結果Fig.6 Results in the third ramp
圖6 (a)和圖6(b)給出了第三個坡道內3個算例的壁面摩擦系數沿流向的分布、流向平均速度剖在靠近出口處沿法向方向的分布。可以看出在入口處添加湍流較大幅值擾動后,經過一個快速的調整段,壁面摩擦系數在x=0.67m開始趨近湍流狀態的值,但如果擾動的幅值較小,壁面摩擦系數任然會趨近于層流狀態的值。從流向平均速度剖面同樣也可以看出,情況A和B的剖面很飽滿,而情況C的剖面在靠近壁面處趨近于層流解。
圖6(c)給出了第三個坡道內情況B中經過Van Direst變換的平均速度剖面,可以看出在x=0.68m到x=0.70m的平均速度剖面幾乎重合,符合平均速度的相似性特性。并且平均速度剖面與線性律和對數律都吻合的很好,說明情況B確實達到了湍流狀態。
采用空間模式直接數值模擬的方法,仿真了帶有三個壓縮折角的高超聲速飛行器下表面的湍流發生器,分析了其引起強迫轉捩的可能性,發現:
(1)在第一個坡道內不會發生自然轉捩。
(2)在第一個坡道上,擾動脈動速度均方根值小于0.15時,流場趨近于層流狀態,不會發生強迫轉捩;如果大于0.225時,流場有趨近湍流狀態的趨勢,有可能發生強迫轉捩。
(3)在第二個和第三個坡道上,擾動脈動速度均方根值為0.15時,流場趨近湍流狀態,有可能發生強迫轉捩,但如果小于0.03時,流場有趨近層流狀態的趨勢,不會發生強迫轉捩。
[1] TIRTEY S C,CHAZOT O.Characterization of hypersonic roughness induced transition for the EXPERT flight experiment[R].AIAA 2009-7215.
[2] BERRY S,DARYABEIGI K.Preliminary analysis of boundary layer transition on X-43Aflight 2[A].52nd JANNAF Propulsion Meeting[C],2004.
[3] BERRY S,DARYABEIGI K,WURSTER K.Boundary layer transition on X-43A[R].AIAA 2008-3736.
[4] LIU J T.Nonlinear instability of developing stream-wise vortices with applications to boundary layer heat transfer intensification through an extended Reynolds analogy[J].Phil.Trans.R.Soc.A.,2008,366:2699-2716.
[5] MOMAYEZ L,DUPONT P,PEERHOSSAINI H.Some unexpected effects of wavelength and perturbation strength on heat transfer enhancement by Go¨rtler instability[J].International Journal of Heat and Mass Transfer,2004,47:3783-3795.
[6] 黃章峰,曹偉,周恒.超聲速平板邊界層轉捩中層流突變為湍流的機理--時間模式[J].中國科學(G 輯),2005,35(5):537-547.
[7] 黃章峰,周恒.湍流邊界層的空間模式DNS的入口條件[J].中國科學(G輯),2008,38(3):310-318.
[8] 董明,周恒.超聲速鈍錐湍流邊界層DNS入口邊界條件的研究[J].應用數學和力學,2008,29(8):893-904.