邢博陽,楊俊彥,余 躍,張 斌

(1.上海交通大學 航空航天學院·上?！?00240;2.上海航天控制技術研究所·上?！?01109)

0 引 言

臨近空間通常是指海拔從20km到100km不等的區(qū)域,近年來其已成為世界軍事強國的重要研究方向。臨近空間的空氣密度很小,并且隨高度呈指數(shù)下降趨勢[1]。例如,當海拔為40km時,空氣密度約為海平面標準密度的1/300;當海拔為60km時,空氣密度約為海平面標準密度的1/4000;當海拔為80km時,空氣密度接近真空。近幾十年來,許多學者通過實驗或數(shù)值模擬研究了后臺階的特征,并發(fā)現(xiàn)了流場中的一些典型結構,如膨脹波、剪切層、回流區(qū)、重新附著點,以及重新附著沖擊等[2-5]。朱楊柱[6]發(fā)現(xiàn),隨著馬赫數(shù)的增加,再附著點位置逐漸向后移動。陳植[7]通過實驗獲得了馬赫數(shù)為3的超音速層流和后臺階湍流結構,結果表明回流區(qū)的長度約為臺階高度的7.5倍。Eaton[8]回顧了亞音速湍流后臺階的再附著點,并指出回流區(qū)的長度在6~8倍的臺階高度之間。工程中出現(xiàn)的后臺階構型的流場通常伴隨著噴流制冷現(xiàn)象,目前也有大量針對噴流冷卻的后臺階流場研究。張鋒[9]采用基于納米粒子的平面激光散射技術對不同靜壓比下馬赫數(shù)為3的超聲速冷卻氣膜流場進行了實驗研究。結果表明,在波系結構、噴流厚度及湍流化程度等方面,靜壓比對超聲速冷卻氣膜產(chǎn)生了顯著的影響。朱志斌[10]以光學窗口外冷噴流為研究背景,采用大渦模擬方法對后臺階外形切向噴流混合流場進行了研究,獲得了流場結構特征、時空演化規(guī)律及流場密度脈動特性。鄧放[11]采用高精度格式求解了二維Navier-Stokes方程,研究了超聲速射流與同向超聲速后臺階流動相互作用的流場基本結構及規(guī)律,討論和分析了流場中復雜的波系結構及其相互干擾的流動現(xiàn)象。

如上所述,已有大量文獻對亞音速和超音速的后臺階流場進行了實驗和數(shù)值研究。但是,臨近空間的高超音速流動很少受到關注。在超音速和高超音速飛行器的光學頭罩中,通常會出現(xiàn)向后臺階結構,研究近空間高超聲速后臺階的流場對于工程應用而言具有重要意義。

具體而言,本文既關注不同噴流冷卻形式的冷卻效果,又關注相應條件下流場的氣動光學效應。文章使用由課題組開發(fā)并已通過驗證[12-14]的高精度直接模擬蒙特卡洛算法(Direct Simulation Monte Carlo,DSMC)[15]計算臨近空間高超聲速的后臺階流場。此外,提出了一種基于計算網(wǎng)格和直接仿真的光線追蹤技術,用于計算氣動光學效應。結果表明,臨近空間高超聲速后臺階流場在高度大于30km時可忽略不計。本文同時比較了四種噴流冷卻方式的制冷效果和氣動光學效應。

1 數(shù)值方法

1.1 DSMC方法

本文采用DSMC方法計算不同條件下的后臺階流場。DSMC是一種基于Boltzmann方程的數(shù)值模擬方法,采用大量分子模型模擬物理分子的真實運動,以計算宏觀物理量,如速度、溫度、密度、壓力、質(zhì)量等。DSMC方法能夠準確地模擬從自由分子流動區(qū)域到接近連續(xù)流動區(qū)域的大多數(shù)空氣動力學問題,是解決稀薄氣體區(qū)域和過渡區(qū)域中的空氣動力學問題的唯一成功的數(shù)值方法[16]。分子碰撞過程的模擬基于可變硬球模型(Variable Hard Sphere,VHS),化學反應采用量子-動理學(Quantum-Kinetic,Q-K)模型,采用自適應子網(wǎng)格技術滿足碰撞要求,采用變時間步長、metis負載平衡分區(qū)等加速技巧提高程序的計算效率。本程序能夠模擬超聲速、高超聲速的單一組分氣體和多組分混合氣體,并可考慮高溫空氣離解的可壓黏性流動,尤其在稀薄領域具有明顯優(yōu)勢。DSMC方法的工作過程可以概括為如下步驟:

(1)讀取網(wǎng)格數(shù)據(jù),記錄邊界條件信息;

(2)初始化流場,計算實際分子的進入數(shù)量;

(3)模擬分子運動,并與邊界相互作用;

(4)索引所有模擬分子;

(5)按概率選擇模擬分子,并發(fā)生碰撞;

(6)對單元和壁面信息進行采樣,并重復3~6次,使其穩(wěn)定;

(7)輸出流場和壁面信息。

1.2 光線追蹤方法

為研究光束穿過混合層流場的氣動光學效應,本文采用了基于流場網(wǎng)格和直接模擬的光線追蹤技術。這種方法基于光沿著直線傳播的假設,可以通過光線追蹤模擬光粒子在非均勻介質(zhì)中的運動軌跡,進而得到一系列的氣動光學評價參數(shù)。

粒子在穿越網(wǎng)格時,折射或反射現(xiàn)象均有可能出現(xiàn),如圖1所示。粒子的反射或折射現(xiàn)象與粒子是從光密介質(zhì)進入光疏介質(zhì)還是從光疏介質(zhì)進入光密介質(zhì)有關。具體而言,粒子的反射或折射取決于粒子歷經(jīng)的相鄰兩個網(wǎng)格的折射率大小。

圖1 二維光線追蹤示意圖Fig.1 The ray tracing method in 2D plane

圖2所示的光線追蹤流程可概述為以下幾個步驟:(1)初始化計算,在邊界上布置虛擬光粒子;(2)在給定時間步長內(nèi)根據(jù)牛頓定律遷移粒子;(3)確定粒子是否離開原始網(wǎng)格。如是,轉至(4),否則繼續(xù)遷移;(4)確定粒子是否穿越了計算域邊界,如是,停止,否則計算網(wǎng)格面編號及穿越的位置;(5)根據(jù)折射率更新粒子速度,并繼續(xù)遷移粒子。

圖2 光線追蹤流程圖Fig.2 The computational procedure of ray tracing method

粒子在穿越整個計算流場后,沿其經(jīng)歷的路徑可以得到該粒子的光程(Optical Path Length,OPL),計算公式如式(1)所示。其中,c為光線積分路徑,n是光線所通過每一個網(wǎng)格的折射率,由網(wǎng)格內(nèi)密度換算而來。

光程差(Optical Path Difference,OPD)同樣是光學評價中的重要參數(shù),由式(2)得到,其中尖括號表示空間平均。

2 算例設置和網(wǎng)格無關性驗證

2.1 網(wǎng)格模型

后臺階結構流場(Backward-Facing Step,BFS)具有展向中心對稱的特點,并且展向寬度與臺階高度足夠大,因此本文中的流場模擬可以采用二維模型進行計算,流場計算區(qū)域如圖3所示。后臺階上下游壁面互相平行,后臺階高度h為5mm,上游壁面和臺階下游壁面的長度分別為60mm和100mm。流場上游的流場寬度為40mm,下游流場寬度為45mm。

圖3 后臺階流場示意圖及CFD計算區(qū)域Fig.3 Schematic diagram of the BFS and the computational area of CFD

對于噴流介質(zhì)的選取,從容易獲得且價格低廉等角度考慮,選取氮氣作為噴流介質(zhì)。垂直噴流出口設置在距離后臺階上游9mm位置處,噴縫沿著展向設置,寬度為2mm。噴流速度馬赫數(shù)為1。另外,沿著后臺階豎直立面距離上壁面2mm處設置噴縫寬度為2mm的水平噴流出口,噴流速度馬赫數(shù)為3。所有噴口均為方形。

2.2 算例設置

本文的主要目的是研究臨近空間范圍內(nèi)的后臺階流動問題。臨近空間范圍為20km~80km,文章中主要討論的計算范圍為30km~60km,后臺階流動不同噴流控制的算例以30km流場作為基準。分別采用A、B、C、D四種不同的噴流方式,研究不同噴流方式對后臺階流動控制的影響。

對于后臺階流場,左側和右側邊界是自由流邊界條件,上下游下方邊界是等溫壁面條件,壁面條件為無滑移等溫壁面。給定壁面溫度為400K,噴流溫度為300K,不同位置的噴流可以起到冷卻作用。對于所有的算例,來流條件為馬赫數(shù)6。來流空氣組分為78%的氮氣(N2)和22%的氧氣(O2),空氣的摩爾分子質(zhì)量為28.96(kg/mol),比熱比為1.4。其他詳細參數(shù)(如靜壓、靜溫等參數(shù))如表1所示。其中,H表示飛行高度,P0表示來流壓強,T0表示來流氣體溫度,U0表示遠場來流速度,Ma表示馬赫數(shù)。

表1 不同高度后臺階流場計算條件Tab.1 Boundary conditions for the basic BFS flowfield

不同噴流控制流場的角度設置如圖4所示。噴流方式A、B的噴口在豎直立面,噴流方式A的噴流方向為沿水平方向,噴流方式B的方向為斜向上45°。噴流方式C、D的噴口位置在上游壁面上,噴流方式C的方向為豎直向上,噴流方式D的方向為斜向上45°。不同噴流參數(shù)(如靜壓、靜溫等參數(shù))的設置如表2所示。其中,Pj表示噴流介質(zhì)的壓強,Tj表示噴流氣體的溫度,Uj表示冷卻噴流的速度。

圖4 不同噴流方式示意圖Fig.4 Schematic diagram of different film cooling modes

表2 不同噴流方式的噴流計算條件Tab.2 Computational conditions for the BFS of different film cooling modes

2.3 網(wǎng)格無關性驗證

DSMC方法是直接求解Boltzman方程的方法,采用模擬流體粒子運動的方式對流場進行仿真。DSMC方法對網(wǎng)格尺寸和采樣粒子的數(shù)量均有較高的要求,因此要對后臺階流場問題進行網(wǎng)格無關性驗證。坐標系建立在后臺階拐角處,x、y兩個方向分別為流向和豎直方向。

本文采用81121網(wǎng)格、202521網(wǎng)格和401501網(wǎng)格三種不同加密程度的網(wǎng)格進行計算。網(wǎng)格如圖5所示,其中橫縱坐標分別表示為x/h、y/h,h為無量綱化的后臺階高度。

圖5 后臺階流場計算網(wǎng)格Fig.5 Computational grid for 2D BFS

網(wǎng)格密度沿著后臺階壁面加密,接近壁面的網(wǎng)格尺度為20μm。應用DSMC方法計算后臺階流場需使用三種不同尺寸的網(wǎng)格。后臺階流場的高度為30km,x/h=5處的縱向截面速度分布和溫度分布如圖6所示。由圖6可以看出,三套網(wǎng)格數(shù)值仿真結果中的速度和溫度沒有明顯的差異。放大來看,中等網(wǎng)格和加密網(wǎng)格的結果更加接近,因此可以采用中等網(wǎng)格或加密網(wǎng)格進行流場模擬計算。為了減少計算量,下文所有數(shù)值仿真均采用中等加密的網(wǎng)格進行計算。

3 結果與討論

3.1 不同飛行高度后臺階流動特征

首先研究無噴流后臺階流場在30km~60km的流動特征,以下計算結果均是DSMC算法的時均結果。圖7是不同高度下后臺階流場的密度云圖。從圖7可以明顯看出,隨著飛行高度增加,流場密度不斷降低。高度越高,空氣越稀薄,這是由大氣密度隨高度變化本身決定的?？梢钥吹?在高度為30km、40km時,流場存在明顯的前緣激波、后臺階回流區(qū)、在附著激波等流場結構。同時可以看到,隨著高度增加,流場壓縮性減小,前緣激波角度減小,并且后臺階回流區(qū)的大小也隨著飛行高度的增加而不斷減小。當高度為50km、60km時,流場更加稀薄,導致在附著激波也不再明顯。從顯示流場的流線可以看到,激波角度的增大使得流動速度方向向上抬升了微小的角度。

圖6 三種網(wǎng)格在X/h=5處的流向速度(左)和溫度(右)的網(wǎng)格無關性驗證Fig.6 Comparison of stream wise velocity(left)and temperature(right)along the vertical line ofX/h=5 for three different grid densities

圖7 不同飛行高度后臺階流場的密度云圖和流線Fig.7 The density contour with streamlines for the cases of basic BFS at different altitudes

通過DSMC仿真流場可以得到后臺階下游壁面的熱流系數(shù)CQ,如圖8所示。由圖8可以看到,在高度為30km、40km時,后臺階下游壁面的熱流系數(shù)在回流區(qū)右邊界處出現(xiàn)熱流系數(shù)峰值。在高度為50km、60km時,后臺階流場的來流密度出現(xiàn)了顯著降低,熱流系數(shù)在下游壁面持續(xù)增加?？傮w而言,隨著飛行高度的增加,后臺階流場下游壁面的氣動熱問題愈發(fā)嚴重。

3.2 不同噴流方式后臺階流動特征

采用DSMC方法計算不同噴流方式的流場(來流均為Ma=6),可以得到流場的密度云圖,如圖9所示。從圖9可以看到,對于A、B兩種平行于壁面的冷卻噴流方式,噴流方式A的回流區(qū)偏小,噴流出口的流線沿著噴口向兩邊略微膨脹。噴流經(jīng)過回流區(qū)沿著壁面向下游流動,對壁面起到了冷卻作用。

圖8 不同高度下后臺階流場的熱流系數(shù)Fig.8 CQ for BFS at different altitude

對于噴流方式B,由于噴流方向為斜向上45°,噴流下方的回流區(qū)相比噴流方式A略微增大。同樣,噴流在越過回流區(qū)后貼近下游壁面,可對流動起到冷卻作用,同時與上游自然來流形成剪切層。

對于噴流方式C,噴口設置在距離后臺階豎立面邊緣9mm的位置,垂直噴流與來流形成了典型的噴流結構?？梢钥吹?豎立面左側流場中有分離區(qū),噴口兩側有回流區(qū)、馬赫盤、膨脹激波等復雜的流場結構。自由來流在遇到噴流形成的膨脹波后,來流速度發(fā)生了改變,流線出現(xiàn)了抬升,后臺階區(qū)域出現(xiàn)了狹長的回流區(qū)。流場密度最大值約為噴流方式A、C的4.5倍。

對于噴流方式D,噴流角度為斜向上45°,流場相比噴流方式C,噴流兩側的回流區(qū)消失,并且噴流右側的馬赫盤也出現(xiàn)了消失。膨脹波角度相比噴流方式C更小。同樣地,后臺階區(qū)域也出現(xiàn)了狹長的回流區(qū)。流場密度最大值小于噴流方式C,但同樣大于噴流方式A和B。

圖9 不同噴流方式下后臺階流場密度云圖和流線Fig.9 The density contour with streamlines for the cases of BFS

同樣地,求得后臺階下游壁面熱流系數(shù),可以看到四種噴流方式下的壁面熱流總體小于沒有噴流冷卻的標準流場。如圖10所示,噴流方式A、B的熱流系數(shù)接近。在壁面前段,熱流系數(shù)為負,說明熱傳遞方向相反,降溫效果好。噴流方式C的冷卻效果強于噴流方式D,并且兩種噴流方式在上游壁面設置噴口的冷卻方式的冷卻效果均弱于在后臺階豎立面設置噴口的A、B噴流方式。

圖10 不同噴流方式下后臺階流場熱流系數(shù)Fig.10 CQ of BFS with different film cooling modes

3.3 不同高度后臺階氣動光學效應

使用632nm的光束垂直于流場上方,在X/h=0~10平面布撒光粒子模擬光線照射。光束穿過流場最終到達后臺階下游壁面。在整個過程中,光束穿過了前緣激波、在附著激波、回流區(qū)等復雜的流場結構,如圖11所示。

圖11 光束布置位置示意圖Fig.11 Schematic diagram of beam position

采用光線追蹤算法,模擬光粒子在流場中的傳輸,得到光波的波面畸變。計算得到OPD,如圖12所示。

圖12 不同飛行高度下后臺階流場OPD變化Fig.12 OPD of BFS at different altitude

可以看到,對于無噴流后臺階流場而言,飛行高度越高,波面畸變越小。這是由于隨著飛行高度增加,流場整體密度變小,并且不同流場結構本身的密度梯度的減小對光束傳輸?shù)挠绊憸p弱。在高度為40km時,OPD的幅值是高度為30km時OPD幅值的1/5,并且隨著飛行高度的增加,氣動光學OPD的幅值更小。因此,可以認為飛行高度在30km以上時,流場的氣動光學效應嚴重減弱,可以忽略不計。

使用光線追蹤方法可以通過計算得到每一個光粒子穿過非均勻流場的最終位置與無流動介質(zhì)的最終位置之間的距離,如圖13所示。對所有光子的距離偏移進行平均,可以得到在該條件下圖像的平均偏移,該平均偏移可以用來作為衡量氣動光學效應的一個參數(shù)。

圖13 光線通過非均勻介質(zhì)的圖像偏移Fig.13 Image shift of beam passing through inhomogeneous medium

設入射角θi為0.573°,通過計算得到不同飛行高度下每個光粒子圖像偏移量的平均值,如圖14所示。由圖14可以明顯看出,隨著飛行高度的增加,后臺階流場的圖像偏移急劇減小。當飛行高度增加到30km以上時,就可以忽略氣動光學效應。

圖14 不同飛行高度下后臺階流場的圖像偏移量Fig.14 Image shift of BFS at different altitude

3.4 不同噴流方式后臺階氣動光學效應

同樣,對四種不同噴流方式的流場采用上文的光學參數(shù)進行氣動光學計算,并用30km無噴流流場計算結果作為參照,可以得到波面畸變,如圖15所示。

圖15 不同噴流方式下后臺階流場OPD變化Fig.15 OPD of BFS with different film cooling modes

可以看出,對于四種噴流冷卻方式,噴流方式C(即上游壁面垂直噴流)的OPD幅值最大,相比無噴流后臺階流場的幅值大一個量級。同時,其他三種噴流方式的波面畸變明顯大于30km無噴流后臺階流場。

同樣采用光線追蹤計算光束平均偏移,可以得到不同噴流方式流場對圖像偏折的影響,如圖16所示。由圖16可以看到,兩種上游壁面開孔噴流方式C、D的圖像偏移要比A、B兩種噴流方式的圖像偏移更大,并且噴流方式C的圖像偏移遠大于其他三種噴流方式。

圖16 不同噴流方式下的圖像偏移量Fig.16 Image shift of BFS with different film cooling modes

綜上所述,在四種噴流方式控制下的后臺階流場中,垂直噴流方式的氣動光學效應最差,水平噴流和斜45°角水平噴流的氣動光學效應相對較好,水平噴流方式A略好于噴流方式B。

4 結 論

對于不同飛行高度的后臺階流場,飛行高度越高,壁面的氣動熱效應越嚴重,但同時氣動光學效應也隨之減弱。當高度在30km以上時,可以忽略氣動光學效應。

考慮到對后臺階壁面的噴流冷卻,采取了四種不同的噴流方式。其中四種噴流方式相對于無噴流后臺階流場均能夠改善壁面溫度,起到冷卻作用。其中,水平噴流的噴流方式冷卻效果最好。在四種噴流方式中,垂直噴流方式的氣動光學效應最嚴重,水平噴流方式的氣動光學效應相對較弱,但均嚴重于無噴流流場。

綜上所述,對于后臺階流場的冷卻,建議采用水平噴流的冷卻方式,這種方式在冷卻下游壁面的同時不會帶來更為嚴重的氣動光學效應。

參考文獻(References)

[1] HEDIN A E.Extension of the MSIS thermosphere model into the middle and lower atmosphere[J].Journal of Geophysical Research:Space Physics,1991,96(2):1159-1172.

[2]LE H,MOIN P,KIM J.Direct numerical simulation of turbulent flow over a backward-facing step[J].Journal of Fluid Mechanics,1997,330:349-374.

[3]ARMALY B F,DURST F,PEREIRA J C F,et al.Experimental and theoretical investigation of backward-facing step flow[J].Journal of Fluid Mechanics,1983,127:473-496.

[4]BISWAS G,BREUER M,DURST F.Backward-facing step flows for various expansion ratios at low and moderate Reynolds numbers[J].Journal of Fluids Engineering,2004,126(3):362-374.

[5]SHAHI M,KOK J B W,POZARLIK A.On characteristics of a non-reacting and a reacting turbulent flow over a backward facing step(BFS)[J].International Communications in Heat and Mass Transfer,2015,61:16-25.

[6]朱楊柱,易仕和,孔小平,等.基于NPLS的超聲速后臺階流場精細結構及其非定常特性[J].物理學報,2014,63(13):271-276.ZHU Y Z,YI S H,KONG X P,et al.Fine structures and the unsteadiness characteristics of supersonic flow over backward facing step via NPLS[J].Acta Physica Sinica,2014,63(13):271-276(in Chinese).

[7]陳植,易仕和,何霖,等.基于NPLS的超聲速層流/湍流后臺階流動精細結構研究[J].科學通報,2011,56(36):3057-3063.CHEN Z,YI S H,HE L,et al.An experimental study on fine structures of supersonic laminar/turbulent flow over a backward-facing step based on NPLS[J].Chinese Science Bulletin,2011,56(36):3057-3063(in Chinese).

[8]EATON J K,JOHNSTON J P.A review of research on subsonic turbulent flow reattachment[J].AIAA Journal,1981,19(9):1093-1100.

[9]張鋒,易仕和,吳宇陽,等.基于NPLS的靜壓比對高超聲速光學頭罩冷卻氣膜影響作用研究[J].空氣動力學學報,2019,37(5):762-769.ZHANG F,YI S H,WU Y Y,et al.Ratio of static pressure influence on cooling film of hypersonic optical dome based on NPLS[J].Acta Aerodynamica Sinica,2019,37(5):762-769(in Chinese).

[10]朱志斌,程曉麗,潘宏祿.超聲速噴流混合流場大渦模擬[J].航空動力學報,2019,34(1):210-216.ZHU Z B,CHENG X L,PAN H L.Large eddy simulation of supersonic jet mixingflow[J].Journal of Aerospace Power,2019,34(1):210-216(in Chinese).

[11]鄧放,韓桂來,孟寶清,等.超聲速后臺階流動/射流相互作用的數(shù)值模擬[J].氣體物理,2018,3(6):41-50.DENG F,HAN G L,MENG B Q,et al.Numerical simulation of supersonic flow/jet interaction at the double backward-facing step[J].Physics of Gases,2018,3(6):41-50(in Chinese).

[12]LI L Y,REN W,ZHANG B.Improved vertex weight for parallelization of DSMC in the near-continuum regimes[J].Journal of Aeronautics,Astronautics and Aviation,2014,46(2):88-95.

[13]REN W,LIU H,JIN S.An asymptotic-preserving Monte Carlo method for the Boltzmann equation[J].Journal of Computational Physics,2014,276:380-404.

[14]REN W,LIU H.Nonequilibrium hypersonic flows simulations with asymptotic-preserving Monte Carlo methods[C]//Proceedings of 29th International Symposium on Rarefied Gas Dynamics,Xi