張雋研,王學德

(南京理工大學 能源與動力工程學院,江蘇 南京 210094)

臨近空間高超聲速飛行器的氣動布局與外形設計一直是近年來研究的熱點問題。由于臨近空間大氣的特殊性,在該空域飛行的飛行器,其結(jié)構強度及熱力學防護主要受到氣動力和氣動熱的影響,因此,其外形布局的設計尤為重要。在實際應用中,在飛行器的表面往往需要增加一定的空腔結(jié)構來滿足和強化具體的應用需求,如在高超聲速飛行器鈍頭體的前緣點引入空腔結(jié)構,可以使飛行器頭部空腔的外緣形成“冷卻環(huán)”,有效降低了局部熱流[1];高超聲速火箭撬試驗中,為了方便試驗體的安裝與固定,火箭撬下方的滑軌會與試驗體之間形成空腔結(jié)構[2];用于增強高超聲速飛行器機動性的空氣舵,其底部也設計為便于安裝的類似空腔結(jié)構的縫隙[3];制造公差、非相似材料的不均勻碰撞或者傳感器的安裝,導致飛行器的表面不可避免地出現(xiàn)縫隙或缺陷這種類空腔結(jié)構,使飛行器在飛行過程中出現(xiàn)安全隱患[4]。以上應用都基于空腔這一基本的幾何外形。為了更好地解決帶空腔高超聲速飛行器空腔幾何結(jié)構的設計,需要對臨近空間空腔幾何特性及腔內(nèi)流動問題進行研究。

空腔流動問題一直是空氣動力學研究的經(jīng)典問題之一,高速自由流流過空腔時,會形成2種基本的流動結(jié)構,即空腔內(nèi)的剪切層及低速的再循環(huán)區(qū)域[5-6]。根據(jù)剪切層對空腔內(nèi)部再循環(huán)區(qū)域的影響,空腔分為閉式和開式兩大類[7]。開式空腔和閉式空腔內(nèi)流體的流動特性并不相同。對于開式空腔,剪切層和邊界層并未觸及腔體底部,高速自由流對剪切層的沖擊使其產(chǎn)生自持振蕩,這使得空腔周圍的氣動噪聲增大,給腔內(nèi)的部署設施以及空腔自身結(jié)構的安全帶來隱患;對于閉式空腔,邊界層和剪切層逐漸觸及腔體底部,空腔內(nèi)再循環(huán)區(qū)域一分為二,使得空腔底部的壓強增大,對空腔內(nèi)部組件的安全分離產(chǎn)生不利的影響[8]。

國內(nèi)外許多學者對連續(xù)流區(qū)域的空腔進行了很多數(shù)值模擬和實驗研究工作,但是對于稀薄流區(qū)域的空腔研究較少。靳旭紅等[9]研究了70 km,75 km,80 km,90 km 4個不同高度下稀薄流域高超聲速縫隙流動問題,并分析了稀薄氣體效應對縫隙內(nèi)部流動結(jié)構和壁面熱流的影響。Palharini等[4,10]首先對比了高超聲速稀薄流動條件下,二維空腔與三維空腔內(nèi)部的傳熱系數(shù)、壓力系數(shù)與摩阻系數(shù),結(jié)果表明腔體后緣y方向上二維算例的傳熱系數(shù)比三維算例的值大幾倍,而腔體底部x方向的值吻合良好。隨后其又研究了70 km高度處、馬赫數(shù)為25時不同長深比對空腔壁面空氣動力學特性的影響,發(fā)現(xiàn)長深比為4時空腔內(nèi)部流動結(jié)構發(fā)生改變。Guo[11]采用結(jié)構網(wǎng)格,研究了60 km高度處、馬赫數(shù)為8時不同長深比以及不同高度的前、后壁面對空腔內(nèi)再循環(huán)區(qū)域以及剪切層的影響。

至今,學者們對稀薄流區(qū)域的不同高度以及不同馬赫數(shù)下空腔由開式變?yōu)殚]式的具體長深比數(shù)值還沒有給出定量的研究計算,對稀薄流區(qū)域的空腔的研究,給出的高度跨度以及馬赫數(shù)跨度還不夠廣。因此,本文從稀薄流臨近空間出發(fā),以空腔作為研究對象,研究高度為50～90 km,馬赫數(shù)為5～20。由于直角坐標網(wǎng)格在處理傾斜壁面時需要對單元格進行切割處理,對物面的處理相對復雜,基于此,本文利用非結(jié)構網(wǎng)格良好的貼體性以及生成的便利性,研究了不同傾斜角下的空腔前壁面對空腔內(nèi)再循環(huán)區(qū)域的影響。

1 計算方法

1.1 DSMC方法介紹

DSMC 方法是直接研究分子的微觀運動的統(tǒng)計學方法,克服了解玻爾茲曼方程的復雜性。該方法通過大量的模擬分子來仿真真實分子的微觀運動情況。其本質(zhì)是在時間步長Δt內(nèi),通過運動子程序和碰撞子程序?qū)⒎肿拥倪\動和分子的碰撞這2個耦合的過程進行解耦,通過采樣子程序采集到每個流場單元的宏觀流動數(shù)據(jù),最后按單元編號將數(shù)據(jù)輸出。DSMC方法中,網(wǎng)格起到以下2個作用[12]:一是空間離散化宏觀流場;二是有利于選取分子碰撞對。

本文計算網(wǎng)格采用非結(jié)構貼體網(wǎng)格,分子模型采用變徑硬球(VHS),能量交換模型采用Larsen-Borgnakke模型,物面采用恒溫壁面邊界條件,反射模型采用Cercignani-Lampis-Lord(CLL)完全漫反射模型。

1.2 算例驗證

本文采用圓柱繞流模擬算例來驗證程序的有效性,圓柱繞流的示意圖如圖1所示,其中Kn為克努森數(shù),用于表征空氣的稀薄程度。計算條件取自文獻[13],模擬氣體為氬氣,模擬分子個數(shù)為151萬,計算網(wǎng)格數(shù)為61 336,程序循環(huán)步數(shù)為20 000。來流速度v、來流溫度T、來流密度ρ、物面溫度Tsurf、馬赫數(shù)以及時間步長Δt的設置見表1。

圖1 圓柱繞流示意圖

表1 計算參數(shù)

v/(m·s-1)T/KMaρ/(kg·m-3)Tsurf/KΔt/μs2 634.1200102.818×10-55000.01

圖2給出了按表1參數(shù)條件運行的密度、溫度和壓力的等值線圖。圖3給出了本文圓柱繞流的壓力系數(shù)Cp、摩阻系數(shù)Cf和表面?zhèn)鳠嵯禂?shù)Ch計算結(jié)果與文獻[13]的比較,圖中的橫坐標Φ為圓柱表面的切線與x軸正方向所成的角度。由圖可以看出,本文的計算結(jié)果與文獻值吻合較好,驗證了本文程序采用非結(jié)構網(wǎng)格計算的正確性。

圖2 圓柱繞流的密度、溫度、壓力等值線圖

圖3 圓柱繞流計算結(jié)果與文獻[13]結(jié)果的比較

2 空腔的幾何特性研究與分析

2.1 模擬參數(shù)的設置

定義空腔底部長度L與深度D的比值為長深比。本文取L/D=1～10,長深比間隔為1,以及前壁面傾斜角θ=30°,45°,60°,75°,研究空腔內(nèi)再循環(huán)區(qū)域的流動特征,模擬高度H=50 km,60 km,70 km,80 km,90 km,來流速度Ma=5,10,15,20。不同高度下的來流條件見表2,表中,n為分子的數(shù)密度。由于高度的不同,空氣的稀薄程度也不盡相同,網(wǎng)格尺度約為分子自由程的1/3,時間步長與來流速度的乘積要小于網(wǎng)格尺度,每個網(wǎng)格中的模擬分子數(shù)約為20～30。表3給出了L/D=6時不同高度算例中網(wǎng)格數(shù)K、分子數(shù)N與時間步長Δt的具體設置。本文研究計算的空腔模型如圖4所示。

表2 不同高度高超聲速來流條件

表3 不同高度,L/D=6算例的計算參數(shù)

圖4 空腔幾何形狀示意圖

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

DSMC方法的計算精度主要取決于網(wǎng)格密度、采樣次數(shù)、每個網(wǎng)格單元中的粒子和時間步長。根據(jù)經(jīng)驗,采樣時間、每個網(wǎng)格單元中的粒子數(shù)量和時間步長都設置為適當?shù)闹?僅研究網(wǎng)格尺度的影響。時間步長Δt、網(wǎng)格數(shù)K以及分子數(shù)N的設置見表4。

表4 網(wǎng)格驗證計算參數(shù)

定義一個無量綱參數(shù)Ω來研究網(wǎng)格尺度對計算結(jié)果的影響,其表達式為

(1)

式中:a為修正系數(shù),這里取值為1.009;s為網(wǎng)格單元的平均線性尺寸,即相鄰節(jié)點之間的平均距離;λ為來流分子的平均自由程;Kn為來流的克努森數(shù);l為流場的特征長度。以50 km處L/D=6,Ma=10,取空腔前沿X=-0.3 m至X=0 m處物面的體積熱流密度Q的算例為例,計算3套Ω值分別為1,0.67和0.33的網(wǎng)格,來驗證3種尺度網(wǎng)格的獨立性,計算結(jié)果如圖5所示。3種不同尺度網(wǎng)格的局部放大圖如圖6所示?？梢钥闯?每個網(wǎng)格尺度的熱流密度很好地相互吻合,Ω=0.33和Ω=0.67尺度的網(wǎng)格計算結(jié)果幾乎重合。為了減小計算損耗,選擇中等網(wǎng)格作為計算網(wǎng)格處理下文中的各種工況。

圖5 3種不同網(wǎng)格尺度的熱流密度對比

圖6 3種不同尺度網(wǎng)格的局部放大示意圖

2.3 不同長深比(L/D)對空腔流場的影響

首先以60 km的高度為例,空腔的長深比為1～10,間隔為1。圖7給出了高度為60 km時不同長深比下馬赫數(shù)的等值線圖,圖中以流線的形式標識了空腔中的再循環(huán)區(qū)域以及剪切層中的流動。可以看到,當長深比從1逐漸增加到10時,邊界層也經(jīng)歷著顯著的變化,剪切層區(qū)域為圖中的淺灰色區(qū)域,該區(qū)域通過流線標識,如圖7的圖例所示。

圖7 60 km處不同長深比的馬赫數(shù)等值線及再循環(huán)區(qū)域示意圖

當長深比小于4時,空腔內(nèi)的再循環(huán)區(qū)域僅有一個,隨著邊界層的擠壓效應,空腔內(nèi)的再循環(huán)區(qū)域逐漸地一分為二。空腔內(nèi)當L/D≥6時,剪切層逐漸觸及到空腔底板,使再循環(huán)區(qū)域分為2個獨立的循環(huán)區(qū)域,此時稱空腔由開式變?yōu)殚]式。

文獻[14]中提出,長深比通常用于對空腔進行分類,以長深比L/D=10為空腔變形極限值,當L/D>10時空腔由開式變?yōu)殚]式。根據(jù)本文的計算結(jié)果,在臨近空間高度為60 km處,當L/D>6時,空腔呈閉式。也就是說,文獻[14]中提到的用于分類空腔類型的標準不再適用于稀薄高超聲速流動的情況。文獻[15]中,空腔深度保持在3 mm,長深比分別為1,2,3,4,5,模擬所用的自由流動條件為:來流速度7 600 m/s,高度80 km,在這種情況下,即使L/D=3,空腔依然會變形,而L/D=4,5時,空腔變?yōu)殚]式。因此,一個重要的發(fā)現(xiàn)是,用于區(qū)分空腔類型的標準(即L/D)不是恒定的,而是隨自由流條件而變化。將本文的計算結(jié)果與文獻[15]的結(jié)果比較后可以得出,對于具有更高海拔和更大速度的自由流動,L/D的標準會變得更低。本節(jié)旨在找出不同高度及不同馬赫數(shù)時使空腔由開式變?yōu)殚]式的具體長深比(L/D)的值。

本文還計算了其他高度及馬赫數(shù)時空腔由開式變?yōu)殚]式時的臨界值,具體結(jié)果見表5。由于50 km和60 km處Ma=5的工況下空腔發(fā)生類型轉(zhuǎn)變的長深比值均大于10,超出了本文的計算范圍,因此表5中并未給出具體數(shù)值。60～90 km處Ma=10,15,20的閉式空腔流動示意圖如圖8所示。

表5 空腔由開式變閉式的L/D臨界值

從表5中可以看出,隨著高度的增加空腔由開式變?yōu)殚]式的長深比的值在逐漸減小,并且隨著馬赫數(shù)的增加,空腔發(fā)生變形的長深比值也呈減小趨勢,并且高度對空腔發(fā)生形變的影響大于馬赫數(shù)的影響。當Ma=5時,不同高度下空腔發(fā)生變形的長深比臨界值的間隔很大,該間隔隨高度的增加而減小,并且,該間隔隨馬赫數(shù)的增加也呈減小趨勢。H=50 km時,馬赫數(shù)對空腔發(fā)生變形的長深比的數(shù)值的間隔影響也較大,長深比L/D的數(shù)值間隔隨著馬赫數(shù)的增加而減小,并且,隨著高度的增加,長深比的間隔也呈減小趨勢。H=90 km時長深比不再隨馬赫數(shù)的變化而變化,當Ma=20,H>70 km時,L/D的值也不再因高度的變化而發(fā)生改變,因此,鄰近空間區(qū)域內(nèi)空腔由開式變閉式的長深比臨界值的極限為4。

圖8 不同工況下空腔由開式變閉式的流線圖

將空腔內(nèi)部壁面的長度進行歸一化處理,定義前壁面為S1,其無量綱量長度為L1;底部壁面為S2,其無量綱長度為L2;后壁面為S3,其無量綱長度為L3,如圖9所示。60 km處Ma=10的空腔后壁面處表面?zhèn)鳠嵯禂?shù)隨長深比的變化趨勢分別如圖10和圖11所示,從圖中可以看出,隨著長深比的增大,空腔壁面的表面?zhèn)鳠嵯禂?shù)也呈增大趨勢,并且該趨勢與文獻[4]中的結(jié)果吻合。隨著無量綱長度的增加,長深比大于7時表面?zhèn)鳠嵯禂?shù)呈先增大后減小趨勢,原因在于長深比過大時會在空腔后方產(chǎn)生再附激波,對流場產(chǎn)生一定影響,并且分子撞擊后壁面會反射使得近壁面處速度降低,也會使熱流密度下降。

圖9 空腔壁面無量綱長度示意圖

圖10 空腔底部面壁S2的表面?zhèn)鳠嵯禂?shù)

2.4 不同物面傾角對空腔內(nèi)流場的影響

以60 km處,Ma=10,L/D=7的工況為例,通過觀察空腔內(nèi)的再循環(huán)區(qū)域,研究物面傾角對其流動特性產(chǎn)生的影響。

不同傾斜角下馬赫數(shù)的等值線圖如圖12所示,物面傾角分別設置為θ=30°,45°,60°,75°。從圖中可以看出,剪切層對左側(cè)空腔底部的擠壓效應隨左側(cè)物面傾斜角的增大而減小。傾角θ=30°時,邊界層幾乎完全觸及空腔左側(cè)底部,使左側(cè)的再循環(huán)區(qū)域完全消失,而右側(cè)再循環(huán)區(qū)域幾乎不受影響。當θ>60°時,左側(cè)出現(xiàn)較小的再循環(huán)區(qū)域。隨著傾角的增大,左側(cè)再循環(huán)區(qū)域也逐漸擴大,因此在保證空腔的空氣動力學性能的前提下,減小壁面傾角可以有效地減少再循環(huán)區(qū)域?qū)涨坏挠绊?。以本文研究?nèi)容為例,空腔壁面傾角為30°為宜。將前、后壁面的傾角同時設置為30°,馬赫數(shù)等值線圖如圖13所示,由圖可見,空腔內(nèi)流動穩(wěn)定,沒有出現(xiàn)再循環(huán)區(qū)域。

圖12 不同壁面傾角下的馬赫數(shù)等值線圖

圖13 前、后壁面傾角均為30°時的馬赫數(shù)等值線圖

空腔前壁面不同傾角的表面?zhèn)鳠嵯禂?shù)如圖14所示。

圖14 不同壁面傾角時壁面S1的表面?zhèn)鳠嵯禂?shù)

當θ≥60°時,該壁面處表面?zhèn)鳠嵯禂?shù)的大小由空腔左側(cè)的再循環(huán)區(qū)域決定,因此θ=75°算例的表面?zhèn)鳠嵯禂?shù)大于θ=60°算例;當θ≤45°時,該壁面處表面?zhèn)鳠嵯禂?shù)的大小由自由來流決定。與θ=45°的算例相比,θ=30°算例的空腔前壁面受到自由來流的影響更大,因此該算例下空腔前壁面的表面?zhèn)鳠嵯禂?shù)更大。

3 結(jié)論

本文研究了不同高度以及不同馬赫數(shù)下不同的長深比L/D對空腔內(nèi)在循壞區(qū)域的影響,以及在固定的長深比下,不同的物面傾角θ對空腔內(nèi)在循環(huán)區(qū)域的影響。主要得出以下結(jié)論:

①找出了不同高度和馬赫數(shù)下空腔內(nèi)流動結(jié)構發(fā)生形變的L/D的臨界值,發(fā)現(xiàn)高度和速度的增加可以使該臨界值變小,且L/D=4為該臨界值的極限值。

②減小物面的傾角可以使剪切層更易觸及空腔底部,有效地消除空腔內(nèi)的再循環(huán)區(qū)域,30°以下的角度為消除再循環(huán)區(qū)域且防止再次生成再循環(huán)區(qū)域的適合角度。

③空腔底部壁面與后壁面的表面?zhèn)鳠嵯禂?shù)隨長深比L/D的增大而增大。當壁面傾角θ≥60°時,空腔前壁面的表面?zhèn)鳠嵯禂?shù)隨θ的增大而增大;當壁面傾角θ≤45°時,空腔前壁面的表面?zhèn)鳠嵯禂?shù)隨θ的增大而減小。