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

張雋研,王學(xué)德

(南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094)

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

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

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

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

1 計(jì)算方法

1.1 DSMC方法介紹

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

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

1.2 算例驗(yàn)證

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

圖1 圓柱繞流示意圖

表1 計(jì)算參數(shù)

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

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

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

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

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

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

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

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

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

圖4 空腔幾何形狀示意圖

2.2 網(wǎng)格無關(guān)性驗(yàn)證

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

表4 網(wǎng)格驗(yàn)證計(jì)算參數(shù)

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

(1)

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

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

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

2.3 不同長(zhǎng)深比(L/D)對(duì)空腔流場(chǎng)的影響

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

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

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

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

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

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

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

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

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

圖9 空腔壁面無量綱長(zhǎng)度示意圖

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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