臨近空間高超聲速空腔幾何特性的DSMC研究

張雋研,王學德

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

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

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

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

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

1 計算方法

1.1 DSMC方法介紹

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

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

1.2 算例驗證

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

圖1 圓柱繞流示意圖

表1 計算參數

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

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

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

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

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

2.1 模擬參數的設置

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

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

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

圖4 空腔幾何形狀示意圖

2.2 網格無關性驗證

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

表4 網格驗證計算參數

定義一個無量綱參數Ω來研究網格尺度對計算結果的影響,其表達式為

(1)

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

圖5 3種不同網格尺度的熱流密度對比

圖6 3種不同尺度網格的局部放大示意圖

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

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

圖7 60 km處不同長深比的馬赫數等值線及再循環區域示意圖

當長深比小于4時,空腔內的再循環區域僅有一個,隨著邊界層的擠壓效應,空腔內的再循環區域逐漸地一分為二。空腔內當L/D≥6時,剪切層逐漸觸及到空腔底板,使再循環區域分為2個獨立的循環區域,此時稱空腔由開式變為閉式。

文獻[14]中提出,長深比通常用于對空腔進行分類,以長深比L/D=10為空腔變形極限值,當L/D>10時空腔由開式變為閉式。根據本文的計算結果,在臨近空間高度為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時,空腔變為閉式。因此,一個重要的發現是,用于區分空腔類型的標準(即L/D)不是恒定的,而是隨自由流條件而變化。將本文的計算結果與文獻[15]的結果比較后可以得出,對于具有更高海拔和更大速度的自由流動,L/D的標準會變得更低。本節旨在找出不同高度及不同馬赫數時使空腔由開式變為閉式的具體長深比(L/D)的值。

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

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

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

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

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

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

圖10 空腔底部面壁S2的表面傳熱系數

2.4 不同物面傾角對空腔內流場的影響

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

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

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

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

空腔前壁面不同傾角的表面傳熱系數如圖14所示。

圖14 不同壁面傾角時壁面S1的表面傳熱系數

當θ≥60°時,該壁面處表面傳熱系數的大小由空腔左側的再循環區域決定,因此θ=75°算例的表面傳熱系數大于θ=60°算例;當θ≤45°時,該壁面處表面傳熱系數的大小由自由來流決定。與θ=45°的算例相比,θ=30°算例的空腔前壁面受到自由來流的影響更大,因此該算例下空腔前壁面的表面傳熱系數更大。

3 結論

本文研究了不同高度以及不同馬赫數下不同的長深比L/D對空腔內在循壞區域的影響,以及在固定的長深比下,不同的物面傾角θ對空腔內在循環區域的影響。主要得出以下結論:

①找出了不同高度和馬赫數下空腔內流動結構發生形變的L/D的臨界值,發現高度和速度的增加可以使該臨界值變小,且L/D=4為該臨界值的極限值。

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

③空腔底部壁面與后壁面的表面傳熱系數隨長深比L/D的增大而增大。當壁面傾角θ≥60°時,空腔前壁面的表面傳熱系數隨θ的增大而增大;當壁面傾角θ≤45°時,空腔前壁面的表面傳熱系數隨θ的增大而減小。