高超聲速壁湍流入口條件生成方法的比較

禹 旻,袁湘江,朱志斌

(中國航天空氣動力技術(shù)研究院,北京 100074)

高超聲速壁湍流入口條件生成方法的比較

禹 旻*,袁湘江,朱志斌

(中國航天空氣動力技術(shù)研究院,北京 100074)

給定恰當(dāng)?shù)娜肟跅l件是開展壁湍流數(shù)值模擬的關(guān)鍵問題。選取基于有限差分的直接數(shù)值模擬方法結(jié)合高精度的計算格式，詳細(xì)討論了自然轉(zhuǎn)捩、波紋壁面結(jié)構(gòu)導(dǎo)致的“Bypass”轉(zhuǎn)捩和利用時間發(fā)展湍流場進(jìn)行參數(shù)回收促發(fā)轉(zhuǎn)捩這幾種湍流入口生成方法在高超聲速條件下的可行性，分析了這幾種方法各自存在的優(yōu)點和不足之處。計算結(jié)果表明，Bypass轉(zhuǎn)捩和利用時間發(fā)展流場進(jìn)行參數(shù)回收方法相比自然轉(zhuǎn)捩能快速促發(fā)轉(zhuǎn)捩，但自然轉(zhuǎn)捩得到的湍流場品質(zhì)更好。該項研究為高馬赫數(shù)壁湍流數(shù)值模擬入口條件的選取提供參考依據(jù)。

高超聲速壁湍流；入口條件；自然轉(zhuǎn)捩；波紋壁面；流場參數(shù)回收

0 引 言

湍流的多樣性使得描述一個湍流流動入口是非常困難的。對于湍流的數(shù)值模擬，通常需要在入口給定復(fù)雜的湍流脈動生成方法來模擬真實的湍流場。常用的湍流入口邊界生成方法有自然轉(zhuǎn)捩方法、合成湍流方法和流場參數(shù)回收方法等等[1-3]。

自然轉(zhuǎn)捩方法最為簡單[4]，但計算量大，計算域包含了層流區(qū)、轉(zhuǎn)捩區(qū)和湍流區(qū)。流場參數(shù)回收方法是采用計算本身來產(chǎn)生入口所需的湍流脈動信息，它其中又包括利用時間發(fā)展湍流場方法(流向周期性方法)[5-7]和回收/調(diào)節(jié)方法[8-9]，該方法適合的研究對象是不關(guān)心轉(zhuǎn)捩區(qū)的充分發(fā)展湍流。在工程上，經(jīng)常使用在壁面上引入粗糙度、凸起或壁面變形等使流動和幾何參數(shù)發(fā)生改變的方法，對轉(zhuǎn)捩產(chǎn)生影響[10-11]。

由于壓縮性的影響，高馬赫數(shù)下邊界層流動更為穩(wěn)定，因此在高超聲速下生成湍流場就更為困難。一些常用的生成湍流場的方法在高馬赫數(shù)條件下不能達(dá)到預(yù)期效果。比如回收/調(diào)節(jié)方法，該方法把下游某處速度場、熱力學(xué)參量和各種湍流統(tǒng)計量按照一定的比例縮放后引入入口，但這樣引入的前提是Morkovin假說成立[9]，其實該方法在M>5時就不是很有效了。在實際數(shù)值計算中也發(fā)現(xiàn)，在低馬赫數(shù)情況下(M=3)使用該方法能產(chǎn)生湍流場，而在高馬赫數(shù)情況下(M=6)就不能激發(fā)出湍流場。

本文主要針對高馬赫數(shù)流動，對自然轉(zhuǎn)捩、波紋壁面結(jié)構(gòu)和利用時間發(fā)展流場的參數(shù)回收這幾種方法進(jìn)行了數(shù)值計算，并對結(jié)果進(jìn)行了分析和評價。

1 控制方程和數(shù)值方法

考慮一般曲線坐標(biāo)系(ξ,η,ζ)，在此坐標(biāo)系下三維Navier-Stokes方程形式為：

采用完全氣體假設(shè)，方程的具體表達(dá)式參見文獻(xiàn)[12]。方程(1)的數(shù)值求解采用有限差分方法，通量采用Sterge-Warming矢通量分裂形式，通量離散波紋壁計算時使用五階WENO格式[13]，其他算例使用五階迎風(fēng)格式；黏性項采用六階中心差分格式；時間方向采用具有TVD性質(zhì)的三步三階Runge-Kutta方法[14]推進(jìn)求解。

邊界處理：壁面采用無滑移條件和等溫壁條件；出口采用外推方法；展向采用周期邊界條件。

2 幾種可壓縮湍流入口條件生成方法

氣體參數(shù)選取26 km高空處值，來流馬赫數(shù)為6，壁溫為600 K。計算時高空來流速度、密度和溫度值作為相應(yīng)無量綱參考量，長度用L=1 m作為無量綱參考量。

2.1 自然轉(zhuǎn)捩方法

在計算域入口存在邊界層問題時，將計算域入口取在上游層流段，根據(jù)Orr-Sommerfeld方程得到不穩(wěn)定特征模態(tài)，疊加至層流剖面上，激勵擾動經(jīng)過線性、非線性增長轉(zhuǎn)捩至成為完全發(fā)展的湍流。

在高馬赫數(shù)情況下，最不穩(wěn)定的是第二模態(tài)二維擾動[15]。為促使擾動快速進(jìn)入非線性增長階段，在計算域入口處加入一個第二模態(tài)二維擾動和一對第一模態(tài)三維擾動，擾動的形式為：

表1 入口擾動參數(shù)(無量綱量，自然轉(zhuǎn)捩方法)Table 1 Parameters of the disturbances at the inlet (dimensionless quantities,natural transition)

圖1中黑色實線給出了自然轉(zhuǎn)捩方法得到的物面摩擦系數(shù)分布，虛線是層流和湍流壁面摩擦系數(shù)的理論值[16]。在x=0.55 m左右曲線開始抬升，轉(zhuǎn)捩開始，經(jīng)過一段距離后達(dá)到湍流狀態(tài)。圖2給出了計算到統(tǒng)計定常后速度梯度張量二階不變量Q2=5000的三維等值面，可以看到，計算域入口附近流動是以二維結(jié)構(gòu)為主，隨著向下游的推進(jìn)，開始有三維結(jié)構(gòu)出現(xiàn)，Λ渦形成，在流向x>0.9 m到達(dá)完全湍流區(qū)。圖3給出了流向為x=0.9 m處的平均速度剖面沿壁面法向的分布，實線為平均速度，三角形為Van Driest變換后的平均速度，虛線分別為壁面律(u+=y+)和對數(shù)律(u+=2.5ln(y+)+5.8)，平均速度剖面與對數(shù)率和壁面率符合得很好。在計算域入口處添加上述擾動，數(shù)值模擬獲得了平板邊界層從層流失穩(wěn)轉(zhuǎn)捩到湍流的整個發(fā)展過程。

2.2 壁面上引入波紋結(jié)構(gòu)

在結(jié)構(gòu)表面上引入粗糙單元或改變壁面局部結(jié)構(gòu)，都可使局部流場狀態(tài)發(fā)生改變，促使流向渦的產(chǎn)生，使流動剪切作用增強，起到促進(jìn)轉(zhuǎn)捩的效果。

在流向位置為x1=0.45 m至x2=0.5 m之間布置波紋壁面，波紋壁的分布函數(shù)為：

其中：高度h=2.45 mm(約為0.7倍當(dāng)?shù)剡吔鐚雍穸?，kx=5×2π/(x2-x1)，kz=2×2π/lz，lz為展向長度。

圖4給出了Q2=5000的三維等值面，在下游由于波紋壁面的干擾作用，使得流動剪切作用增強，流場變得不穩(wěn)定，層流快速轉(zhuǎn)捩到湍流，從圖4中可以看出x≈0.6 m左右開始流動結(jié)構(gòu)就比較飽滿了。圖1中實心方形表示的是波紋壁方法得到的物面摩擦系數(shù)曲線，由于波紋壁的存在，能引發(fā)壁面摩擦系數(shù)曲線迅速抬升，在經(jīng)過波紋壁面之后，曲線會緩慢下降，經(jīng)過一段距離后達(dá)到完全湍流狀態(tài)。圖5給出了x=0.68 m處的平均速度剖面，可以看到平均速度剖面同對數(shù)率和壁面率符合得較好。

2.3 利用時間發(fā)展湍流場進(jìn)行參數(shù)回收方法

采用時間發(fā)展湍流場作為空間發(fā)展流場的入口條件，預(yù)先模擬的時間發(fā)展湍流場流向采用周期性邊界條件，在進(jìn)行回收時只需用預(yù)先模擬的湍流場的某一時刻的流動信息即可，預(yù)存的數(shù)據(jù)量小。

首先，計算一個隨時間發(fā)展的湍流場。與2.1節(jié)計算自然轉(zhuǎn)捩時引入擾動的方法類似，初始時刻在整個計算域中引入一些擾動，擾動形式為：

具體參數(shù)見表2，此處ωi代表擾動的增長率。

表2 擾動的具體參數(shù)(無量綱量,時間發(fā)展的DNS)Table 2 Parameters of the disturbances at the inlet (dimensionless quantities,temporal DNS)

圖6給出了壁面摩擦系數(shù)曲線隨時間的分布?？梢钥吹絫=0.2 s左右壁面摩擦系數(shù)曲線開始抬升，也就意味著轉(zhuǎn)捩的開始，到t=0.5 s左右達(dá)到峰值進(jìn)入湍流區(qū)，到達(dá)湍流區(qū)后系數(shù)曲線的變化趨于平緩且有所下降。圖7為0.8 s時速度梯度張量二階不變量Q2=5000的三維等值面，此時已是充分發(fā)展湍流。取0.8 s時刻的瞬時流場數(shù)據(jù)作為空間發(fā)展湍流流場的入口條件，假定此瞬時湍流流場數(shù)據(jù)為UT(x,y,z)，則空間模式DNS入口條件可表示為：

這里通過時空變換關(guān)系，x=ct，將某一時刻隨時間發(fā)展的湍流場整個計算域的流場信息轉(zhuǎn)換為空間發(fā)展的湍流流場入口處的時間序列。黃章峰等[5]認(rèn)為，壁湍流是具有大尺度相干結(jié)構(gòu)的，而相干結(jié)構(gòu)對湍流的性質(zhì)起著決定性的作用，而相干結(jié)構(gòu)具有某些波的特點，其平均傳播速度約為0.92，在實際計算中也發(fā)現(xiàn)設(shè)定為此對流速度能得到較好的效果，因此取傳播速度c為0.9。

圖1中的空心圓形為回收方法得到的壁面摩擦系數(shù)沿流向的分布，可看到由于入口條件的人為給定，流動需要經(jīng)過一小段調(diào)整，之后就維持在完全湍流狀態(tài)。圖8為Q2=5000的三維等值面，計算域入口開始就是湍流狀態(tài)，而且在向下游的推進(jìn)中能維持住湍流流動。圖9給出了流向為0.55 m處平均速度剖面，同對數(shù)率和壁面率符合得較好。

3 可壓縮湍流入口條件生成方法評價

對第2節(jié)中計算出的湍流場進(jìn)行歸納總結(jié)，從生成湍流場品質(zhì)和能否快速生成湍流場兩方面來對這幾類入口條件給定方法進(jìn)行評價。首先，圖10分別給出了y+=40時3種方法生成的湍流場速度分布云圖，自然轉(zhuǎn)捩方法生成的湍流場最為均勻飽滿；引入波紋壁面結(jié)構(gòu)促發(fā)轉(zhuǎn)捩生成的湍流場不如自然轉(zhuǎn)捩方法得到的湍流場均勻、飽滿；同樣，使用隨時間發(fā)展湍流場進(jìn)行參數(shù)回收方法時得到的湍流場也不如自然轉(zhuǎn)捩方法得到的湍流場飽滿。

另從圖1中的壁面摩擦系數(shù)曲線對比圖，可以直觀地看到引入波紋壁面結(jié)構(gòu)方法相比自然轉(zhuǎn)捩方法能更為快速的生成湍流場，不需要經(jīng)歷很長的流向距離(從流動入口位置0.4 m到完全湍流處0.7 m左右)。而采用時間發(fā)展流場進(jìn)行回收作為入口條件的方法，從計算域入口開始就是湍流場，只要經(jīng)過一小段的調(diào)整(入口0.4 m到0.45 m左右)就能獲得一個比較真實的充分發(fā)展湍流，此種方法也是一種快速生成湍流場的方法。

采用時間發(fā)展流場進(jìn)行回收方法可以推廣到其它外形和不同計算工況。董明等[6]證實此方法是對不同馬赫數(shù)、雷諾數(shù)、壁面條件和典型邊界層類型均可行的有效的空間模式直接數(shù)值模擬入口條件生成方法，在實際數(shù)值計算中也發(fā)現(xiàn)，使用文獻(xiàn)中對平均密度和溫度以及脈動量進(jìn)行相應(yīng)修正的方法，可導(dǎo)出不同計算條件下的入口條件。在這里用2.3節(jié)計算得到的隨時間發(fā)展的湍流場作為入口，可以計算出來流馬赫數(shù)7、高空27 km、壁面溫度為700 K時的湍流場。圖11為Q2=5000的三維等值面，可以看到同樣從計算域入口開始就是湍流狀態(tài)。圖12為壁面摩擦系數(shù)沿流向的分布，可看到由于入口條件的影響，流動需要經(jīng)過一段比之前相同計算工況條件下稍長一些的調(diào)整距離之后就維持在完全湍流狀態(tài)?？梢?，采用隨時間發(fā)展湍流場進(jìn)行參數(shù)回收方法，可以推廣到更高的馬赫數(shù)。

4 結(jié) 論

數(shù)值模擬研究了自然轉(zhuǎn)捩、波紋壁誘發(fā)轉(zhuǎn)捩、利用時間發(fā)展邊界層湍流這三類促發(fā)湍流場生成的入口條件給定方法。研究表明，在高馬赫流動下這三種方法都可以實現(xiàn)湍流場的生成。其中自然轉(zhuǎn)捩方法處理起來最為簡單，得到的湍流場也最為均勻飽滿，但是計算域長，網(wǎng)格量大，不能快速生成湍流場；波紋壁面的引入，能很快誘導(dǎo)湍流生成，是一種能快速生成湍流的方法，但生成的湍流場不如自然轉(zhuǎn)捩方法得到的湍流場均勻飽滿；利用時間發(fā)展流場進(jìn)行參數(shù)回收這種湍流生成方法，從計算域入口就是湍流場，也是一種能有效生成湍流場的方法，該方法還有一個好處就是可由同一預(yù)先模擬的時間發(fā)展湍流場導(dǎo)出不同馬赫數(shù)、雷諾數(shù)、壁面條件所需要的湍流入口條件，這對于要計算很多參數(shù)下的湍流情況是很有用的。

[1]Lund T S,Wu X H,Squires K D.Generation of turbulent in flow data for spatially-developing boundary layer simulations[J].Journal of Computational Physics,1998,140(2):233-258.

[2]Sivasubramanian J,Fasel H F.Direct numerical simulation of transition in a sharp cone boundary layer at Mach 6:Fundamental breakdown[J].Journal of Fluid Mechanics,2015,768:175-218.

[3]Xu S,Martin M P.Assessment of inflow boundary conditions for compressible turbulent boundary layers[J].Physics of Fluids,2004,16(7):2623-2639.

[4]Batten P,Goldberg U,Chakravarthy S.Interfacing statistical turbulence closures with large-eddy simulation[J].AIAA Journal,2004,42(3):485-492.

[5]黃章峰,周恒.湍流邊界層的空間模式DNS的入口條件[J].中國科學(xué)(G輯:物理學(xué) 力學(xué) 天文學(xué)),2008,38(3):310-318.

[6]董明,周恒.超聲速鈍錐湍流邊界層DNS入口邊界條件的研究[J].應(yīng)用數(shù)學(xué)和力學(xué),2008,29(8):893-904.

[7]逯學(xué)志,黃章峰,羅紀(jì)生.湍流入流條件在強迫轉(zhuǎn)捩研究中的應(yīng)用[J].空氣動力學(xué)學(xué)報,2012,30(6):749-753.

[8]陳逖,劉衛(wèi)東,范曉檣,等.“回收/調(diào)節(jié)”方法在混合LES/RANS模擬方法中的應(yīng)用[J].航空動力學(xué)報,2011,26(6):1215-1222.

[9]Xiao X D,Edwards J R,Hassan H A.inflow boundary conditions for hybrid large eddy/Reynolds averaged Navier-Stokes simulations[J].AIAA Journal,2003,41(8):1481-1489.

[10]Schneider S P.Effects of roughness on hypersonic Boundary-Layer transition[R].AIAA 2007-305,2007.

[11]Berry S A.Discrete roughness transition for hypersonic flight vehicles[R].AIAA 2007-307,2007.

[12]劉建新.小攻角鈍錐高超聲速邊界層的擾動演化[D].天津:天津大學(xué),2010.

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

[14]Gottlieb S,Shu C W.Total variation diminishing Runge-Kutta schemes[J].Mathematics of computation,1998,67(221):73-85.

[15]Yu M,Luo J S.Nonlinear interaction mechanisms of disturbances in supersonic flat-plate boundary layers[J].Science China:Physics Mechanics and Astronomy,2014,57(11):2141-2151.

[16] White F M.Viscous fluid flow[M].McGraw-Hill Education,2005.

Assessmentofinflowboundaryconditionsforhypersonicwallboundedturbulentflows

YU Min*,YUAN Xiangjiang,ZHU Zhibin

(ChinaAcademyofAerospaceAerodynamics,Beijing100074,China)

It is critically important to generate turbulent fluctuating inflow boundary conditions for numerical simulation in solving hypersonic turbulent boundary layers.A high-order finite-difference method was applied to investigate the inflow generation methods.Turbulent inflow boundary conditions are broadly classified into natural transition methods,recycling/rescaling-based methods,and synthetic turbulence generators.To investigate the feasibility of inflow generation methods at high Mach number,the simulation of a Mach 6 boundary layer over an isothermal flat plate was carried out.The validity and limitations of different inflow boundary conditions were also discussed.Several feasible methods were studied in great detail including natural transition from laminar to turbulent flows,bypass transition induced by wavy wall configuration,and turbulence generation via reintroducing the flow field obtained from temporal developing turbulent flow.The limitation of the recycling/rescaling method was also proposed in hypersonic flow condition.By analyzing the flow filed and statistical results,we summarized and discussed the advantages and drawbacks about these inflow generation approaches,and further assessed these methods.Compared with natural transition method,bypass transition method and reintroducing method can promote transition to turbulent flow field efficiently.But turbulent field has better quality generated by natural transition route.The paper can be a reference of imposing an inflow boundary condition for wall bounded turbulent flows at high Mach number.

hypersonic wall bounded turbulent flows; inflow boundary condition; natural transition; wavy wall; reintroducing methods

0258-1825(2017)06-0772-05

O357.5

A

10.7638/kqdlxxb-2015.0177

2015-09-10;

2015-12-15

國家自然科學(xué)基金(11772316)

禹旻*(1986-),女,天津人,高級工程師,研究方向:湍流直接數(shù)值模擬,穩(wěn)定性分析.E-mail:yumin8688@126.com

禹旻,袁湘江,朱志斌.高超聲速壁湍流入口條件生成方法的比較[J].空氣動力學(xué)學(xué)報,2017,35(6):772-776.

10.7638/kqdlxxb-2015.0177 YU M,YUAN X J,ZHU Z B.Assessment of inflow boundary conditions for hypersonic wall bounded turbulent flows[J].Acta Aerodynamica Sinica,2017,35(6):772-776.